曾偉生 孫鄉(xiāng)楠 王六如 王 威 蒲 瑩
(國家林業(yè)和草原局調(diào)查規(guī)劃設(shè)計(jì)院 北京 100714)
激光雷達(dá)(light detection and ranging,LiDAR)即光探測與測量,是一種通過發(fā)射激光來測定傳感器與目標(biāo)物之間距離的主動遙感技術(shù)(趙峰等,2008; 李增元等,2016; 龐勇等,2019)。按探測目標(biāo)不同,激光雷達(dá)可分為對空和對地兩大類,而對地激光雷達(dá)按傳感器搭載平臺不同,又分為星載、機(jī)載、車載和定位4類,其中機(jī)載激光雷達(dá)最為常用(劉魯霞,2014; 李增元等,2016)。機(jī)載激光雷達(dá)能夠獲取較大范圍森林的三維掃描數(shù)據(jù),可用于蓄積量等主要森林參數(shù)的估計(jì)(趙峰等,2008; 劉魯霞,2014; 李增元等,2016)。
在國外,激光雷達(dá)用于估計(jì)森林蓄積量始于20世紀(jì)80年代中期(MacLeanetal.,1986),經(jīng)過30多年的研究和實(shí)踐,激光雷達(dá)技術(shù)取得了長足發(fā)展,尤其是基于機(jī)載激光雷達(dá)數(shù)據(jù)估計(jì)森林蓄積量的研究已經(jīng)獲得了很多成果(N?sset,1997; Holmgrenetal.,2003; Holmgren, 2004; Hollausetal.,2009; Gonzlez-Ferreiroetal.,2012; Bouvieretal.,2015; Bottalicoaetal.,2017),為遙感技術(shù)在森林資源調(diào)查領(lǐng)域的應(yīng)用開辟了新的途徑(Hyypp?etal.,2012; Penneretal.,2013; Whiteetal.,2016; Z?rneretal.,2018)。而在我國,激光雷達(dá)在林業(yè)方面的應(yīng)用相對較晚(劉魯霞,2014; 李增元等,2016; 龐勇等,2019),目前基于機(jī)載激光雷達(dá)估計(jì)森林蓄積量基本還處于研究階段,尚未見到具有良好實(shí)用性的研究成果(付甜,2010; 范鳳云,2010; 楊明星等,2019)。
森林蓄積量是納入《綠色發(fā)展指標(biāo)體系》和《生態(tài)文明建設(shè)考核目標(biāo)體系》的森林資源約束性指標(biāo)之一,也是《自然資源調(diào)查監(jiān)測體系構(gòu)建總體方案》中特別強(qiáng)調(diào)的重要森林資源指標(biāo)。如何充分發(fā)揮激光雷達(dá)技術(shù)優(yōu)勢,服務(wù)于新時期森林資源調(diào)查監(jiān)測工作,是當(dāng)前面臨的一項(xiàng)重要課題。具體到基于機(jī)載激光雷達(dá)數(shù)據(jù)估計(jì)森林蓄積量,盡管國內(nèi)外已經(jīng)取得了一些成果,但研建的蓄積量模型僅適用于特定區(qū)域,不具有普適性。從森林資源調(diào)查的實(shí)用性角度出發(fā),研究建立變量相同、結(jié)構(gòu)穩(wěn)定、在較大地域范圍內(nèi)具有普適性或通用性的模型,是當(dāng)前需要重點(diǎn)解決的問題(Bouvieretal.,2015)。鑒于此,本研究利用東北林區(qū)落葉松(Larixspp.)林、紅松(Pinuskoraiensis)林、樺樹(Betulaspp.)林、楊樹(Populusspp.)林 4種森林類型790個樣地的激光雷達(dá)數(shù)據(jù)和地面實(shí)測蓄積量數(shù)據(jù),首先采用多元線性回歸和非線性回歸方法,分別建立基于機(jī)載激光雷達(dá)數(shù)據(jù)的森林蓄積量回歸模型,通過對比分析,確定具有相同變量和統(tǒng)一結(jié)構(gòu)形式的普適性模型; 然后采用啞變量建模方法(Zengetal.,2011; 鄭冬梅等,2013; Zeng,2015),建立基于相同激光雷達(dá)變量的不同森林類型蓄積量模型,以期為規(guī)范森林蓄積量建模與評價提供科學(xué)參考,為森林資源調(diào)查提供計(jì)量依據(jù)。
東北林區(qū)地處黑龍江、吉林和內(nèi)蒙古 3 省(區(qū)),包括大興安嶺、小興安嶺、完達(dá)山、張廣才嶺、長白山等山系,地理位置119°36′—134°05′E、41°37′—53°33′N,年降水量400~1 000 mm。東北林區(qū)是我國森林資源主要集中分布區(qū)之一,針葉林主要有落葉松林、樟子松(Pinussylvestrisvar.mongolica)林、紅松林和云(Piceaasperata)冷(Abiesfabri)杉林,闊葉林主要有樺樹林、櫟樹(Quercusspp.)林、楊樹林、榆樹(Ulmusspp.)林、椴樹(Tiliaspp.)林和水胡黃林[包括水曲柳(Fraxinusmandshurica)、胡桃楸(Juglansmandshurica)、黃菠蘿(Phellodendronamurense)],另外還有一些混交林類型。本研究從針葉林和闊葉林中各選取2種主要森林類型作為研究對象,其中針葉林選取落葉松林和紅松林,闊葉林選取樺樹林和楊樹林,未考慮結(jié)構(gòu)較復(fù)雜的混交林。
1.2.1 地面樣地調(diào)查數(shù)據(jù) 地面樣地調(diào)查數(shù)據(jù)來自821塊樣地,其中落葉松林202塊、紅松林200塊、樺樹林205塊、楊樹林214塊。調(diào)查范圍覆蓋大興安嶺、小興安嶺、完達(dá)山、張廣才嶺、長白山等林區(qū),并選擇具有典型代表性的12個片區(qū),調(diào)查時間為2019年9—11月。樣地為600 m2圓形樣地,中心點(diǎn)用RTK技術(shù)定位,定位精度在1 m以內(nèi)。除每株樣木測量胸徑外,還測量15株不同徑階的樣木樹高,以此為基礎(chǔ)建立樹高-胸徑回歸模型,推算每株樣木樹高,并依據(jù)部頒二元立木材積表計(jì)算樣木材積,從而得到樣地的蓄積量,作為建模的目標(biāo)變量(Y,m3·hm-2)。
1.2.2 機(jī)載激光雷達(dá)數(shù)據(jù) 機(jī)載激光雷達(dá)數(shù)據(jù)的獲取范圍與地面樣地調(diào)查相同,獲取時間為2019年9—10月。采用RIEGL-VUX-1UAV型激光掃描系統(tǒng),其基本參數(shù)為:精度10 mm,最大測距范圍920 m,最大發(fā)射頻率550 kHz,最大有效測量速率每秒500 000 次; 數(shù)據(jù)點(diǎn)云密度每平方米大于4個點(diǎn),現(xiàn)地定位精度0.11 m。在對原始數(shù)據(jù)進(jìn)行預(yù)處理的基礎(chǔ)上,利用商業(yè)化專用軟件LiDAR360進(jìn)行激光雷達(dá)點(diǎn)云數(shù)據(jù)的分類、平差以及數(shù)字高程模型(DEM)、數(shù)字表面模型(DSM)和冠層高度模型(CHM)的處理,最后提取反映林分高度、密度和結(jié)構(gòu)信息的98個變量作為建模的解釋變量,其中高度變量56個(X01~X56)、強(qiáng)度變量42個(X57~X98)(為省篇幅,具體含義不詳列,后面會對用到的變量給出解釋)。
1.2.3 異常數(shù)據(jù)處理 在建模前,對821塊樣地中的異常數(shù)據(jù)進(jìn)行剔除處理。由于涉及的解釋變量個數(shù)達(dá)98個之多,很難逐一繪制目標(biāo)變量與解釋變量之間的散點(diǎn)圖來判定和剔除異常數(shù)據(jù),因此本研究首先采用全部數(shù)據(jù)建模,然后根據(jù)主要解釋變量的殘差圖來判定和剔除異常數(shù)據(jù),并確保剔除異常數(shù)據(jù)比例不超過5%。最后參與建模的樣地數(shù)為790塊,其中落葉松林197塊、紅松林197塊、樺樹林202塊、楊樹林194塊。表1為4種森林類型參與建模樣地的林分年齡、平均樹高和蓄積量的變動范圍。
表1 建模樣地主要林分因子的變動范圍Tab.1 The ranges of main stand variables for modeling plots
利用4種森林類型790塊樣地的機(jī)載激光雷達(dá)數(shù)據(jù)和地面實(shí)測蓄積量數(shù)據(jù),基于從特殊到普遍、再從普遍到特殊的建模思路,首先采用多元線性回歸和非線性回歸方法,分別建立線性和非線性回歸模型,通過對比分析,確定具有相同變量和統(tǒng)一結(jié)構(gòu)形式的普適性模型; 然后采用啞變量建模方法,建立基于相同激光雷達(dá)變量的不同森林類型蓄積量模型。
1.3.1 線性回歸模型 采用逐步回歸方法確定有顯著相關(guān)的解釋變量(付甜,2010; Bottalicoaetal.,2017),擬合如下多元線性回歸模型:
Y=a0+a1X1+a2X2+ … +akXk+ε。
(1)
式中:Y為蓄積量(m3·hm-2);X1、X2、…、Xk為從98個機(jī)載激光雷達(dá)信息變量中挑選出來的解釋變量;a0、a1、…、ak為模型參數(shù),參數(shù)估計(jì)值的t值原則上應(yīng)大于2,否則視為無統(tǒng)計(jì)學(xué)意義,應(yīng)從模型中剔除;ε為誤差項(xiàng),假定其服從均值為0的正態(tài)分布。
1.3.2 非線性回歸模型 在線性模型(1)基礎(chǔ)上,建立如下非線性回歸模型(Bouvieretal.,2015):
(2)
式中:b0、b1、…、bk為非線性模型參數(shù)。
模型(2)參數(shù)求解采用非線性加權(quán)回歸估計(jì)方法,以消除異方差的影響。
通過對非線性回歸模型和線性回歸模型的擬合結(jié)果及其評價指標(biāo)進(jìn)行對比分析,確定影響森林蓄積量估計(jì)的主要激光雷達(dá)信息變量,并提出反映普遍規(guī)律的模型作為基礎(chǔ)模型。
1.3.3 啞變量模型 確定普適性基礎(chǔ)模型后,由于模型具有相同變量和統(tǒng)一結(jié)構(gòu)形式,不同森林類型的蓄積量模型只是參數(shù)估計(jì)值存在差異,因此采用啞變量建模方法(Zengetal.,2011; 鄭冬梅等,2013; Zeng,2015),同時建立多種森林類型的蓄積量模型。以二元非線性模型為例,其表達(dá)式如下:
Y=(∑b0i·Si)·X1(∑b1i·Si)·X2(∑b2i·Si)+ε。
(3)
式中:Si為反映不同森林類型的啞變量(i=1、2、3、4);b0i、b1i、b2i為不同森林類型參數(shù)。
啞變量的賦值方法為:對于落葉松林樣地,S1=1,S2=S3=S4=0; 對于紅松林樣地,S2=1,S1=S3=S4=0; 對于樺樹林樣地,S3=1,S1=S2=S4=0; 對于楊樹林樣地,S4=1,S1=S2=S3=0。模型(3)參數(shù)求解方法同模型(2)。
1.3.4 模型評價 采用確定系數(shù)(R2)、估計(jì)值的標(biāo)準(zhǔn)誤差(standard error of the estimate,SEE)、總體相對誤差(total relative error,TRE)、平均系統(tǒng)誤差(mean system error,MSE)、平均預(yù)估誤差(mean predictive error,MPE)和平均百分標(biāo)準(zhǔn)誤差(mean percentage standard error,MPSE)對模型進(jìn)行評價(曾偉生等,2011; Zeng,2015),其中MPE和MPSE的計(jì)算公式如下:
(4)
(5)
對于建立的蓄積量回歸模型,計(jì)算以上6項(xiàng)指標(biāo),根據(jù)指標(biāo)大小進(jìn)行模型評價。另外,殘差圖也是評價模型的重要參考依據(jù),一般要求殘差呈隨機(jī)分布。
多元線性模型(1)和非線性模型(2)的解釋變量個數(shù)(k)及6項(xiàng)評價指標(biāo)見表2。
表2 線性和非線性蓄積量模型的評價指標(biāo)Tab.2 Evaluation indices of linear and nonlinear stand volume models
由表2可知,2種針葉林類型的非線性模型相比線性模型要簡單一些,非線性模型的解釋變量只有3或4個,線性模型的解釋變量達(dá)4或6個; 2種闊葉林類型的線性和非線性模型,都只包含2個解釋變量。從評價指標(biāo)對比來看,2種闊葉林類型均為非線性模型好于線性模型,2種針葉林類型由于線性模型的解釋變量個數(shù)較多,其R2、SEE和MPE略好于非線性模型,TRE幾乎無差異,但MSE和MPSE則為非線性模型明顯好于線性模型,且這2項(xiàng)指標(biāo)幾乎所有線性模型均出現(xiàn)異常情況,究其原因是估計(jì)值出現(xiàn)極小值甚至負(fù)值導(dǎo)致的。綜合來看,非線性模型(2)的擬合效果要好于線性模型(1),因此優(yōu)先考慮采用非線性模型,其具體擬合結(jié)果見表3。
表3 非線性蓄積量模型(2)的擬合結(jié)果①Tab.3 Fitting results of nonlinear stand volume model(2)
關(guān)于非線性模型(2),4種森林類型中有2種僅包含2個解釋變量,一是點(diǎn)云高度變量(X35或X37),二是點(diǎn)云強(qiáng)度變量(X80或X82)。如果將落葉松林和紅松林的非線性模型從三元和四元模型簡化為二元模型,其R2分別從0.727、0.818減至0.676、0.805,MPE分別從3.95%、2.89%增至4.29%、2.97%。從模型簡約性和實(shí)用性考慮,這樣的變化幅度是可以接受的,因此選定二元非線性模型作為統(tǒng)一形式的蓄積量模型。
進(jìn)一步分析二元模型中包含的點(diǎn)云高度變量和強(qiáng)度變量,一個為顯著正相關(guān)變量(X35、X37或X38),一個為顯著負(fù)相關(guān)變量(X80、X81或X82)。將4種森林類型790塊樣地數(shù)據(jù)合在一起進(jìn)行相關(guān)分析發(fā)現(xiàn),正相關(guān)最大的變量為X35(點(diǎn)云平均高度),負(fù)相關(guān)最大的變量為X90(點(diǎn)云平均強(qiáng)度),表4為正負(fù)相關(guān)最大的各10個變量。
表4 與蓄積量正負(fù)相關(guān)最大的20個激光雷達(dá)變量Tab.4 20 LiDAR variables having the greatest positive and negative relation with stand volume
從表4可以看出,正相關(guān)最大的10個變量集中分布在3個區(qū)段:X35~X37(點(diǎn)云高度平均值、二次冪平均值、三次冪平均值)、X23~X26(點(diǎn)云高度50%、60%、70%、75%分位數(shù))、X06~X08(點(diǎn)云累積高度30%、40%、50%分位數(shù)); 負(fù)相關(guān)最大的10個變量也集中分布在3個區(qū)段:X90~X91(點(diǎn)云強(qiáng)度平均值、中位數(shù))、X77~X83(點(diǎn)云強(qiáng)度30%、40%、50%、60%、70%、75%、80%分位數(shù))、X46(點(diǎn)云高度偏斜度)。尤其需要重點(diǎn)關(guān)注的激光雷達(dá)變量,是位于正負(fù)相關(guān)2個頂端的變量X35和X90。利用790塊樣地數(shù)據(jù)擬合基于X35和X90的二元模型,其R2達(dá)0.790,高于落葉松林、樺樹林和楊樹林3類蓄積量模型的R2,僅低于紅松林模型的R2。由于X35和X90是最符合預(yù)期的激光雷達(dá)點(diǎn)云信息變量,故本研究將基于這2個變量的二元蓄積量模型定義為標(biāo)準(zhǔn)模型,其相應(yīng)的啞變量模型為:
Y=(∑b0i·Si)·X35(∑b1i·Si)·X90(∑b2i·Si)+ε。
(6)
取1/X35為權(quán)函數(shù),采用加權(quán)回歸方法擬合啞變量模型(6),其R2達(dá)0.866,4種森林類型的參數(shù)估計(jì)值和模型評價指標(biāo)見表5。
表5 啞變量模型(6)的參數(shù)估計(jì)值和模型評價指標(biāo)Tab.5 Parameter estimates and evaluation indices of dummy model(6)
由表5可知,采用變量相同、結(jié)構(gòu)統(tǒng)一的標(biāo)準(zhǔn)二元模型估計(jì)4種森林類型蓄積量,其效果與表2中模型(2)差異不大,模型參數(shù)值大小也較一致,說明該模型穩(wěn)定可靠,具有普適性; 而表3中不同森林類型的非線性模型(2)不僅變量個數(shù)有差異,而且參數(shù)值也相差懸殊。另外,表5中4種森林類型蓄積量模型殘差隨解釋變量分布基本是隨機(jī)的,圖1、2所示為落葉松林和紅松林蓄積量模型相對殘差隨解釋變量X35的分布情況,圖3、4所示為落葉松林和紅松林加權(quán)回歸模型殘差隨解釋變量X35的分布情況,反映了消除異方差后的效果(樺樹林和楊樹林模型殘差也相似,此處略)。
圖1 落葉松林蓄積量模型的相對殘差分布Fig. 1 Relative residual distribution of stand volume model for Larix spp.
圖2 紅松林蓄積量模型的相對殘差分布Fig. 2 Relative residual distribution of stand volume model for Pinus koraiensis
圖3 落葉松林加權(quán)回歸模型的殘差分布Fig. 3 Residual distribution of weighted regression model for Larix spp.
圖4 紅松林加權(quán)回歸模型的殘差分布Fig. 4 Residual distribution of weighted regression model for Pinus koraiensis
激光雷達(dá)作為一種主動遙感技術(shù),為高效估計(jì)森林蓄積量、提高森林資源調(diào)查監(jiān)測工作效率提供了可能(趙峰等,2008; Hyypp?etal.,2012; Penneretal.,2013; Whiteetal.,2016)。利用激光雷達(dá)技術(shù)估計(jì)森林蓄積量的研究最早可追溯到20世紀(jì)80年代(MacLeanetal.,1986),我國將激光雷達(dá)技術(shù)應(yīng)用于森林蓄積量、生物量方面的估計(jì)相對較晚,盡管近10年來已取得了一些研究成果(付甜,2010; 范鳳云,2010; 龐勇等,2012; 劉清旺等,2016; 楊明星等,2019),但具有良好實(shí)用性的成果不多,也尚未形成基于激光雷達(dá)技術(shù)估計(jì)森林蓄積量的應(yīng)用技術(shù)規(guī)范。本研究結(jié)合我國現(xiàn)階段的信息需求,從森林資源調(diào)查監(jiān)測最關(guān)注的蓄積量因子入手,提出了變量相同、結(jié)構(gòu)形式統(tǒng)一的普適性模型,并基于啞變量建模方法,建立了東北林區(qū)4種主要森林類型的蓄積量估計(jì)模型,模型的確定系數(shù)(R2)在0.68~0.81之間,略低于Holmgren(2004)研建的2個蓄積量模型(R2分別為0.82和0.90),與Hollaus等(2009)建立的2個蓄積量模型(R2分別為0.79和0.81)基本一致,與Bottalicoa等(2017)建立的針葉林和櫟樹林蓄積量模型(R2在0.69~0.83之間)也較接近。當(dāng)然,模型的擬合精度與森林結(jié)構(gòu)、研究區(qū)大小也是密切相關(guān)的。
森林蓄積量模型R2高低,在很大程度上還取決于從激光雷達(dá)數(shù)據(jù)中提取的信息變量。激光雷達(dá)信息變量包括不同分位數(shù)高度、不同分位數(shù)密度和最大回波高度等,但不同研究選擇的變量不盡相同(劉魯霞等,2014)。本研究提取反映點(diǎn)云高度和強(qiáng)度的98個變量建立蓄積量回歸模型,進(jìn)入多元線性模型的變量2~6個,進(jìn)入多元非線性模型的變量2~4個。從分析可知,反映林分高度的變量是最重要的,而在眾多高度變量中,體現(xiàn)點(diǎn)云平均高度的變量與蓄積量相關(guān)性最高; 點(diǎn)云強(qiáng)度變量在一定程度上反映了林分的疏密狀況和水平結(jié)構(gòu)特征,與蓄積量大小呈負(fù)相關(guān),即林分高度相同時,株數(shù)越多林分越密,蓄積量就越大,點(diǎn)云強(qiáng)度就越??; 株數(shù)越少林分越疏,點(diǎn)云強(qiáng)度就越大。從實(shí)用性角度出發(fā),非線性二元模型與多元模型沒有實(shí)質(zhì)性差異,而二元模型采用的解釋變量為X35(點(diǎn)云高度平均值)和X90(點(diǎn)云強(qiáng)度平均值),這既是全部98個變量中正負(fù)相關(guān)最大的2個變量,也是最符合預(yù)期的2個變量,所以本研究將基于這2個變量的二元蓄積量模型定義為標(biāo)準(zhǔn)模型。范鳳云(2010)基于激光點(diǎn)云數(shù)據(jù)提取了林分高度平均值、最大值、最小值、標(biāo)準(zhǔn)差4個變量,建立油松(Pinustabulaeformis)和刺槐(Robiniapseudoacacia)的林分蓄積量回歸模型,其唯一的解釋變量就是林分高度平均值; Bouvier等(2015)建立了基于4個變量的通用性模型,其中影響最大的變量是平均冠層高度,與點(diǎn)云平均高度類似; Bottalicoa等(2017)提取了49個變量作為森林蓄積量等主要參數(shù)的解釋變量,所建森林主要參數(shù)回歸模型的變量在1~4個之間,其中與蓄積量相關(guān)性最大的變量是點(diǎn)云高度平均值和三次冪平均值,與本研究結(jié)果基本一致。
森林蓄積量模型的實(shí)用性主要取決于平均預(yù)估誤差(MPE)和平均百分標(biāo)準(zhǔn)誤差(MPSE)的大小,前者衡量對總體蓄積量的估計(jì)誤差,后者衡量對樣地或林分水平蓄積量的估計(jì)誤差。本研究建立的4種森林類型蓄積量模型,MPE在2.89%~4.25%之間,均未超過5%; MPSE在18.13%~24.43%之間,其中紅松林小于20%,其他3種森林類型小于25%。《森林資源規(guī)劃設(shè)計(jì)調(diào)查技術(shù)規(guī)程》(國家質(zhì)量監(jiān)督檢驗(yàn)檢疫總局等,2011)對小班調(diào)查公頃蓄積量的精度分A、B、C三級,要求誤差分別不超過15%、20%和25%??梢姡t松林蓄積量模型可以滿足B級精度要求,其他3種森林類型蓄積量模型可以滿足C級精度要求。采用傳統(tǒng)調(diào)查方法要達(dá)到上述技術(shù)規(guī)程的精度要求是非常費(fèi)工費(fèi)力的,而利用激光雷達(dá)技術(shù)構(gòu)建的蓄積量回歸模型,只需提取相關(guān)信息變量就能得到有一定精度保證的蓄積量估計(jì)值,因此利用激光雷達(dá)技術(shù)開展森林資源規(guī)劃設(shè)計(jì)調(diào)查在技術(shù)上是可行的。
1) 基于機(jī)載激光雷達(dá)數(shù)據(jù)估計(jì)森林蓄積量,非線性模型優(yōu)于線性模型; 從森林資源調(diào)查監(jiān)測的實(shí)用性考慮,非線性多元模型與二元模型無實(shí)質(zhì)性差異。
2) 激光點(diǎn)云高度平均值和強(qiáng)度平均值是與蓄積量正負(fù)相關(guān)最大的2個變量,以此為基礎(chǔ)建立的二元蓄積量模型具有普適性,可作為森林蓄積量估計(jì)的標(biāo)準(zhǔn)模型。
3) 通過引入啞變量代表不同森林類型,可同時建立基于相同激光雷達(dá)變量和統(tǒng)一結(jié)構(gòu)形式的多種森林類型蓄積量模型,是實(shí)踐中值得推廣的一種可行方法。
4) 本研究建立的東北林區(qū)落葉松林、紅松林、樺樹林和楊樹林蓄積量模型,其預(yù)估精度均達(dá)到森林資源調(diào)查相關(guān)技術(shù)規(guī)定要求,可在實(shí)踐中推廣應(yīng)用。