曾偉生,孫鄉(xiāng)楠,王六如,王 威,蒲 瑩

(國家林業(yè)和草原局調(diào)查規(guī)劃設(shè)計院,北京 100714)

激光雷達(dá)是一項新興的主動遙感技術(shù),它通過發(fā)射激光來測定傳感器與目標(biāo)物之間的距離[1-3]。激光雷達(dá)分為對空和對地兩大類,而對地激光雷達(dá)又可分為星載、機載、車載、定位4類,其中最為常用的是機載激光雷達(dá)系統(tǒng)[2,4]。機載激光雷達(dá)通過獲取較大范圍內(nèi)森林的三維掃描數(shù)據(jù),可用來估計蓄積量等主要森林參數(shù)[1,2,4]。由于激光雷達(dá)在估計林木高度和林分結(jié)構(gòu)方面具有獨特優(yōu)勢,自20世紀(jì)80年代中期[5]首次成功用于估計森林蓄積量以后,基于激光雷達(dá)開展森林蓄積量估計方面的研究,已經(jīng)取得了很多成果[6-14],國外已經(jīng)開始將其應(yīng)用于森林資源清查[15-18]。我國開展激光雷達(dá)技術(shù)應(yīng)用于林業(yè)方面的研究相對較晚[2-4],目前,利用機載激光雷達(dá)技術(shù)估計森林蓄積量仍基本處于研究階段[10,11,14,19],尚未見到其應(yīng)用于森林資源調(diào)查實踐的報道。

森林蓄積量是森林資源調(diào)查的核心指標(biāo),已作為兩項森林資源約束性指標(biāo)之一納入生態(tài)文明建設(shè)評價考核指標(biāo)體系。激光雷達(dá)技術(shù)的出現(xiàn),為森林資源調(diào)查領(lǐng)域的森林蓄積量估計開辟了一條新的途徑。在國家機構(gòu)改革的新形勢下,如何充分發(fā)揮激光雷達(dá)技術(shù)的優(yōu)勢,為新時期森林資源調(diào)查監(jiān)測體系的發(fā)展提供技術(shù)支撐,是當(dāng)前面臨的一項重要課題。本研究將利用東北林區(qū)10種森林類型的1 947個樣地的地面實測蓄積量數(shù)據(jù)和機載激光雷達(dá)數(shù)據(jù),通過利用多元線性回歸、非線性回歸及啞變量建模方法[20-22],建立基于相同激光雷達(dá)變量和不同森林類型參數(shù)的蓄積量模型,既為森林資源調(diào)查提供計量依據(jù),也為規(guī)范森林蓄積量建模與評價方法提供科學(xué)參考。

1 數(shù)據(jù)與方法

1.1 數(shù)據(jù)資料

本研究所用數(shù)據(jù)來自中國國土勘測規(guī)劃院的招投標(biāo)項目“主要樹種航空林分材積表編制”(編號:GXTC-A-19070081),包括地面樣地調(diào)查數(shù)據(jù)和機載激光雷達(dá)數(shù)據(jù)兩部分。

1)地面樣地調(diào)查數(shù)據(jù)。該數(shù)據(jù)涉及東北林區(qū)的云冷杉林(Piceaspp.&Abiesspp.)、落葉松林(Larixspp.)、樟子松林(Pinussylvestrisvar.mongolica)、紅松林(Pinuskoraiensis)、櫟樹林(Quercusspp.)、樺樹林(Betulaspp.)、楊樹林(Populusspp.)、榆樹林(Ulmusspp.)、椴樹林(Tiliaspp.)和水胡黃林(Fraxinusmandshurica,Juglansmandshurica&Phellodendronamurense)等10種主要森林類型,共計2 020個樣地,其中每個類型不少于200個。調(diào)查范圍覆蓋內(nèi)蒙古大興安嶺、黑龍江大興安嶺、伊春、松花江、牡丹江、長白山等6個林區(qū)的12個片區(qū),調(diào)查時間為2019年9月至11月。樣地為600m2圓形樣地,中心點用RTK技術(shù)定位,定位精度在1m以內(nèi)。除每株樣木測量胸徑以外,還測量了15株不同徑階的樣木樹高,以此為基礎(chǔ)建立樹高-胸徑回歸模型,推算每株樣木的樹高,并依據(jù)部頒二元立木材積表計算樣木材積,從而得到樣地的每公頃蓄積量,作為建模的目標(biāo)變量(Y,m3/hm2)。

2)機載激光雷達(dá)數(shù)據(jù)。該數(shù)據(jù)范圍與地面樣地調(diào)查相同,獲取時間為2019年9月至10月,數(shù)據(jù)點云密度大于4個點/m2。在對原始數(shù)據(jù)進(jìn)行預(yù)處理的基礎(chǔ)上,利用LiDAR360軟件進(jìn)行激光雷達(dá)點云數(shù)據(jù)的分類、平差,以及數(shù)字高程模型(DEM)、數(shù)字表面模型(DNS)和冠層高度模型(CHM)的處理,最后提取反映林分高度、密度和結(jié)構(gòu)信息的98個變量作為建模的解釋變量,其中高度變量56個(X01—X56),強度變量42個(X57—X98)。

在建模之前,需要對2 020個樣地中的異常數(shù)據(jù)作剔除處理。由于涉及的解釋變量個數(shù)達(dá)到98個之多,很難逐一繪制目標(biāo)變量與解釋變量之間的散點圖來判定和剔除異常數(shù)據(jù),本研究首先采用全部數(shù)據(jù)建模,再根據(jù)主要解釋變量的殘差圖來判定和剔除異常數(shù)據(jù),但剔除異常數(shù)據(jù)的比例原則上要控制在5%以內(nèi),最后參與建模的樣地數(shù)為1 947個。表1為10種森林類型參與建模樣地的每公頃蓄積量的特征值。

表1 建模樣地每公頃蓄積量的特征值

1.2 建模方法

基于前述10種森林類型1 947個樣地的機載激光雷達(dá)數(shù)據(jù)和地面實測數(shù)據(jù),首先建立多元線性和非線性回歸模型,并通過對比分析,研究確定合適的基礎(chǔ)回歸模型；然后利用啞變量建模方法,建立基于相同激光雷達(dá)變量和不同森林類型參數(shù)的蓄積量標(biāo)準(zhǔn)模型,并編制林分蓄積量表。

1.2.1多元線性回歸

采用逐步回歸方法確定有顯著相關(guān)的激光雷達(dá)信息變量,擬合如下多元線性回歸模型:

Y=a0+a1X1+a2X2+…+akXk+ε

(1)

式中:Y為目標(biāo)變量(每公頃蓄積量,m3/hm2)；X1,X2,…,Xk為解釋變量；a0,a1,…,ak為模型參數(shù),參數(shù)的t值原則上應(yīng)該大于2,否則視為無統(tǒng)計學(xué)意義(即與0無顯著差異),從模型中剔除；ε為誤差項,假定其服從均值為0的正態(tài)分布。

1.2.2非線性回歸

在線性模型(1)基礎(chǔ)上,再建立非線性回歸模型:

Y=b0X1b1X2b2…Xkbk+ε

(2)

式中:b0,b1,…,bk為非線性模型參數(shù),其他符號的意義同模型(1)；模型(2)的參數(shù)求解采用非線性回歸估計方法。

通過對多元線性模型與非線性模型的相關(guān)評價指標(biāo)進(jìn)行對比分析,確定與森林蓄積量相關(guān)最緊密的激光雷達(dá)信息變量,并提出統(tǒng)一形式的模型作為基礎(chǔ)回歸模型。

1.2.3啞變量模型

在確定蓄積量基礎(chǔ)模型后,采用啞變量建模方法[20-22],通過引入啞變量代表不同的森林類型,建立基于相同激光雷達(dá)變量和不同森林類型參數(shù)的蓄積量標(biāo)準(zhǔn)模型。以二元非線性模型為例,其表達(dá)式如下:

Y=(∑b0i×Si)×X1(∑b1i×Si)×X2(∑b2i×Si)+ε

(3)

式中:Si為反映不同森林類型的啞變量(i=1,2,…,10)；b0i,b1i,b2i為不同森林類型的參數(shù)；其他符號同前。模型(3)的參數(shù)求解方法同模型(2)。

1.2.4模型評價

用于模型評價的指標(biāo)包括以下6項:確定系數(shù)R2、估計值的標(biāo)準(zhǔn)差SEE、總體相對誤差TRE、平均系統(tǒng)誤差A(yù)SE、平均預(yù)估誤差MPE和平均百分標(biāo)準(zhǔn)誤差MPSE[22-23]。其中MPE和MPSE的計算公式如下:

(4)

(5)

對建立的蓄積量回歸模型,計算以上6項指標(biāo)值,根據(jù)指標(biāo)大小進(jìn)行模型評價。從實用性角度考慮,一般要求模型的TRE和ASE均在±5%以內(nèi),MPE小于5%,MPSE小于20%。另外,殘差圖也是評價模型的重要參考依據(jù)。一個好的模型,殘差應(yīng)當(dāng)呈隨機分布。

2 結(jié)果與分析

利用10種森林類型1 947個樣地的蓄積量數(shù)據(jù)和激光雷達(dá)數(shù)據(jù),分別擬合多元線性模型(1)和非線性模型(2)?；谄紤],此處僅列出模型的評價指標(biāo)如表2所示,每個模型的解釋變量及相應(yīng)的參數(shù)估計值從略。

表2 多元線性和非線性蓄積量模型的評價指標(biāo)

從表2中的解釋變量個數(shù)k可知,非線性模型相比線性模型要簡單一些,多數(shù)森林類型僅用2個變量就能得到滿意的結(jié)果。從表3的評價指標(biāo)對比來看,總體相對誤差TRE幾乎無差異；平均預(yù)估誤差MPE差異也不大,且都達(dá)到了預(yù)期的5%以內(nèi)；確定系數(shù)R2和估計值的標(biāo)準(zhǔn)差SEE線性模型好于非線性模型；但平均系統(tǒng)誤差A(yù)SE和平均百分標(biāo)準(zhǔn)誤差MPSE則大部分非線性模型好于線性模型,尤其是云冷杉林、紅松林、水胡黃林3種類型的ASE和MPSE,其差異更為突出。究其原因,主要是線性模型的估計值出現(xiàn)極小值甚至負(fù)值導(dǎo)致的。綜合來看,非線性模型(2)的擬合效果要好于線性模型(1),特別是可以避免出現(xiàn)負(fù)值。因此,應(yīng)該優(yōu)先考慮采用非線性模型。

關(guān)于非線性模型(2),10種森林類型中有6種僅包含2個解釋變量,其中作用最大的變量是反映激光雷達(dá)點云高度的變量(X24,X26,X35,X36或X44),第二個變量是反映點云強度的變量(X81,X82,X86,X90或X92)。如果將云冷杉林、落葉松林、紅松林和水胡黃林的非線性模型從三元或四元模型簡化為二元模型,其確定系數(shù)R2分別從0.771,0.731,0.846和0.466減少至0.761,0.672,0.821和0.403,平均預(yù)估誤差MPE分別從3.65%,3.71%,2.66%和2.82%增加至3.72%,4.09%,2.85%和2.98%。從模型的實用性考慮,這樣的變化幅度是可以接受的。因此,從模型的簡約性和實用性出發(fā),最終選定二元非線性模型作為基礎(chǔ)回歸模型。

再進(jìn)一步分析二元模型中包含的點云高度變量和強度變量,一個為顯著正相關(guān)變量,另一個為顯著負(fù)相關(guān)變量。將全部10種森林類型的2 020個樣地數(shù)據(jù)合在一起進(jìn)行相關(guān)分析后發(fā)現(xiàn),正相關(guān)最大的變量集中分布在3個區(qū)段:X35—X38(點云高度平均值、二次冪平均值、三次冪平均值、中位數(shù)),X23—X26(點云高度50%,60%,70%,75%分位數(shù)),X04—X07(點云累積高度20%,25%,30%,40%分位數(shù))；負(fù)相關(guān)最大的變量也分布在3個區(qū)段:X90—X91(點云強度平均值、中位數(shù)),X79—X83(點云強度50%,60%,70%,75%,80%分位數(shù)),X46(點云高度偏斜度)(圖1)。這些變量應(yīng)該是需要重點關(guān)注的激光雷達(dá)變量,尤其是排在正負(fù)相關(guān)兩端的變量更是核心變量。

注圖中29個實心柱形對應(yīng)的變量為10種森林類型的線性和非線性模型所涉及的29個解釋變量

經(jīng)建模比較,用正負(fù)相關(guān)最大的2個變量X35和X82(相關(guān)系數(shù)分別為0.709和-0.420)建立的二元模型,其確定系數(shù)反而低于用X35和X90(相關(guān)系數(shù)-0.407)組合的二元模型。由于X35和X90是最符合預(yù)期的激光雷達(dá)點云信息變量,本文把基于這2個變量的二元蓄積量模型定義為標(biāo)準(zhǔn)模型。剔除部分異常樣地后,利用1 947個樣地數(shù)據(jù)擬合基于X35和X90的二元模型,其確定系數(shù)達(dá)到0.675,僅低于表2中云冷杉林、落葉松林和紅松林非線性蓄積量模型的確定系數(shù),高于其他7種森林類型非線性模型的確定系數(shù)。與標(biāo)準(zhǔn)模型相對應(yīng)的啞變量模型為:

Y=(∑b0i×Si)×X35(∑b1i×Si)×X90(∑b2i×Si)+ε

(6)

式中:X35為點云平均高度,X90為點云平均強度,其他符號同模型(3)。啞變量模型(6)的確定系數(shù)達(dá)到0.792,10種森林類型的參數(shù)估計值及模型評價指標(biāo)如表3所示。

表3 啞變量模型(6)的參數(shù)估計值和模型評價指標(biāo)

從表3的各項評價指標(biāo)可知,采用標(biāo)準(zhǔn)模型估計10種森林類型的蓄積量,其擬合效果與表2中的非線性模型(2)差異不大。但是,表2中不同森林類型的非線性模型(2)不僅變量個數(shù)有差異,參數(shù)值也相差懸殊；而表3中的標(biāo)準(zhǔn)模型(6),不同森林類型都是基于相同的2個變量X35和X90(僅水胡黃林例外),每個變量對應(yīng)的參數(shù)值大小也比較一致。所有模型的平均預(yù)估誤差MPE均在5%以內(nèi),總體相對誤差TRE和平均系統(tǒng)誤差A(yù)SE都在±3%以內(nèi)且接近于0,平均百分標(biāo)準(zhǔn)誤差MPSE為20%左右,但模型的確定系數(shù)R2還是存在較大差異,僅紅松林達(dá)到0.8以上,云冷杉林、櫟樹林、楊樹林在0.7～0.8之間,落葉松林、樺樹林、榆樹林在0.6～0.7之間,樟子松林、椴樹林和水胡黃林在0.6以下。尤其是水胡黃林,其標(biāo)準(zhǔn)模型的確定系數(shù)僅為0.111,改用另外2個變量X13和X62后,模型的確定系數(shù)也僅為0.440。究其原因,除了水胡黃林這個森林類型樹種組成較為復(fù)雜、混交程度也大于其他類型外,還有一個重要原因就是樣本的結(jié)構(gòu)不夠理想,這主要也是取樣范圍受到了限制,難以找到合適的樣地。表1中的數(shù)據(jù)也表明水胡黃林的蓄積量變動系數(shù)最小,最小值和最大值之間的分布區(qū)間最窄。

另外,表3中10個森林類型蓄積量模型的殘差隨解釋變量的分布基本是隨機的,圖2和圖3是云冷杉林和落葉松林蓄積量模型的相對殘差隨解釋變量X35的分布情況,其他8個森林類型的殘差圖也類似,此處略。

圖2 云冷杉林蓄積量的相對殘差分布 圖3 落葉松林蓄積量的相對殘差分布

根據(jù)表3中10種森林類型標(biāo)準(zhǔn)蓄積量模型的參數(shù)估計值,就可以編制林分材積表,表4為基于激光雷達(dá)數(shù)據(jù)的云冷杉林二元林分材積表。

表4 云冷杉林二元林分材積表

3 討論與結(jié)論

機載激光雷達(dá)作為一種現(xiàn)代先進(jìn)技術(shù)手段,為準(zhǔn)確估計森林蓄積量,提高森林資源調(diào)查工作效率提供了可能[1,15-17]。本文結(jié)合我國現(xiàn)階段對森林資源調(diào)查監(jiān)測的新要求,從森林蓄積量這一核心因子入手,基于1 947個樣地的地面實測數(shù)據(jù)和機載激光雷達(dá)數(shù)據(jù),利用啞變量建模方法,建立了東北林區(qū)10種主要森林類型的蓄積量估計模型,模型的確定系數(shù)R2在0.44～0.81之間,與Hollaus等[9]、González-Ferreiro等[12]和Bottalicoa等[13]建立的蓄積量模型的確定系數(shù)基本一致,略低于Holmgren等[8]研建的2個蓄積量模型的確定系數(shù)(分別為0.82和0.90),顯著高于范鳳云[10]利用0.39個點/m2的激光雷達(dá)數(shù)據(jù)所建刺槐林和油松林蓄積量模型,其相關(guān)系數(shù)R分別為0.696和0.453。森林蓄積量模型R2的高低,主要取決于模型中采用的激光雷達(dá)數(shù)據(jù)變量。激光雷達(dá)變量包括不同分位數(shù)高度、密度和最大回波高度等,但不同的研究所選擇的變量不盡相同[4]。本研究提取了反映激光點云高度和強度的98個變量來估計蓄積量,最終選擇的模型解釋變量為點云高度平均值X35和點云強度平均值X90,這既是正負(fù)相關(guān)最大的2個激光雷達(dá)變量,也是最符合預(yù)期的2個變量,所以我們將基于這2個變量的蓄積量模型定義為標(biāo)準(zhǔn)模型。Bottalicoa等[13]提取了49個激光雷達(dá)變量作為森林蓄積量的解釋變量,結(jié)果發(fā)現(xiàn)點云高度平均值和三次冪平均值與蓄積量的相關(guān)最高,與本研究的結(jié)果是基本一致的。

森林蓄積量模型的實用性,主要取決于2項誤差指標(biāo)MPE和MPSE的大小,前者反映對總體蓄積量的平均估計誤差,后者反映對林分水平蓄積量的平均估計誤差。本研究所建10個森林類型的蓄積量模型,MPE在2.88%～4.42%之間,均未超過5%；MPSE在16.76%～25.52%之間,其中紅松林、榆樹林、椴樹林、水胡黃林小于20%,另外6個森林類型除樟子松林略超25%以外,其他均小于25%?！渡仲Y源規(guī)劃設(shè)計調(diào)查技術(shù)規(guī)程》[24]對小班調(diào)查蓄積量的精度等級分A,B,C三級,要求相對誤差分別不超過15%,20%,25%。因此,10個森林類型的蓄積量模型有4個可以滿足B級的精度要求,其他6個可以滿足或基本滿足C級的精度要求。采用傳統(tǒng)調(diào)查方法是非常費工費力的,而基于激光雷達(dá)數(shù)據(jù)構(gòu)建蓄積量模型,只需通過信息提取就能得到符合精度要求的蓄積量估計值。因此,利用機載激光雷達(dá)開展森林資源規(guī)劃設(shè)計調(diào)查,在技術(shù)上是完全可行的。

根據(jù)本研究的相關(guān)結(jié)果,可以得出以下結(jié)論:1)利用機載激光雷達(dá)技術(shù)估計森林蓄積量,非線性模型總體上要優(yōu)于線性模型，從森林資源調(diào)查的實用性考慮,二元與多元非線性模型之間無實質(zhì)性差異；2)與森林蓄積量正負(fù)相關(guān)最大的2個激光雷達(dá)變量是點云高度平均值和點云強度平均值,以這2個變量為基礎(chǔ)建立的標(biāo)準(zhǔn)蓄積量模型,可以作為森林蓄積量估計的通用模型；3)基于啞變量建模方法,建立具有不同森林類型參數(shù)和相同激光雷達(dá)變量的蓄積量模型,是林業(yè)建模實踐中值得推廣的一種可行做法；4)本文所建東北地區(qū)10種森林類型的蓄積量回歸模型,其預(yù)估精度達(dá)到《森林資源規(guī)劃設(shè)計調(diào)查技術(shù)規(guī)程》的要求,可以在實踐中推廣應(yīng)用。