基于機載激光雷達(dá)的亞熱帶森林生物量估算模型研究

方志良，陳 奇，任 引，王艷軍

(1. 湖南科技大學(xué)地理空間信息技術(shù)國家地方聯(lián)合工程實驗室，湖南 湘潭 411100；2. 夏威夷大學(xué)馬諾阿分校地理與環(huán)境系，美國 檀香山 96822；3. 中國科學(xué)院城市環(huán)境研究所城市環(huán)境與健康重點實驗室，福建 廈門 361021)

1 研究背景

森林生物量(above-ground biomass, AGB)是森林生態(tài)系統(tǒng)發(fā)揮其生態(tài)功能的物質(zhì)基礎(chǔ)，是森林固碳能力的重要標(biāo)志，在全球碳循環(huán)中扮演著重要角色[1]。對森林生物量和碳貯量現(xiàn)狀及其變化速率的準(zhǔn)確估計，有助于了解森林在區(qū)域和全球碳循環(huán)中的作用，能夠為應(yīng)對氣候變化提供關(guān)鍵信息。特別是我國亞熱帶地區(qū)，水熱條件優(yōu)越，生物多樣性豐富，在中國乃至全球碳循環(huán)中均起到關(guān)鍵作用。近年來，隨著遙感技術(shù)的飛速發(fā)展，利用光學(xué)、雷達(dá)和激光雷達(dá)(light detection and ranging, LiDAR)數(shù)據(jù)進行了大量的森林生物量研究[2-5]。特別是機載LiDAR能夠獲得被動光學(xué)遙感所不能反映的垂直結(jié)構(gòu)信息，且這些結(jié)構(gòu)參數(shù)與森林生物量有著很好的相關(guān)關(guān)系，進而提高森林生物量遙感的估算精度[6]。

利用機載激光雷達(dá)技術(shù)對森林生物量的準(zhǔn)確估算是了解全球碳循環(huán)和緩解氣候變化的重要手段，而異速生長方程及森林生物量估算模型是其中的關(guān)鍵。機載LiDAR數(shù)據(jù)并不能直接計算生物量，首先需要實地野外樣地調(diào)查數(shù)據(jù)結(jié)合異速生長方程對單木生物量進行計算，其次將單木生物量累加得到樣地生物量，然后與LiDAR特征變量構(gòu)建生物量估算模型，最后利用該模型對整個LiDAR覆蓋區(qū)域進行森林生物量的估算。大量研究表明，不同的異速生長方程會較大程度地影響最終生物量的估算精度[7]。如Chen等[8]研究了美國森林清查中的3種不同尺度的異速生長方程在美國太平洋西北部3個研究點之間的差異和關(guān)系，發(fā)現(xiàn)物種的組成和樹木大小均會是引起3種異速生長方程在生物量估算上的差異性。Zhao等[9]利用區(qū)域尺度和國家尺度的兩種異速生長方程比較了8個模型的回歸性能，提出異速生長方程的可用性和不確定性對森林生物量估算的影響。目前，森林生物量估算研究大多是根據(jù)研究區(qū)域選擇已有的異速生長方程或者根據(jù)解析木數(shù)據(jù)擬合的異速生長方程，較少研究不同異速生長方程的選擇對森林生物量估算模型的影響。

近年來，周國逸等[10]對文獻中900多套異速生長方程進行綜合分析，建立了較統(tǒng)一且實用的分省(自治區(qū)、直轄市)分樹種生物量方程(簡稱I)、全國優(yōu)勢種(組)生物量方程(簡稱II)、分省混合種(組)生物量方程(簡稱III)、中國森林主要群系生物量方程(簡稱IV)共4套異速生長方程。但是，4種異速生長方程的差異性及其在亞熱帶區(qū)域適用性研究尚缺乏。從LiDAR數(shù)據(jù)中提取生物量的方法已很成熟。在傳統(tǒng)的統(tǒng)計方法方面，研究人員采用了一系列參數(shù)或非參數(shù)統(tǒng)計模型來進行生物量的量化，包括簡單線性回歸、逐步線性回歸、線性混合模型方法[11-12]。機器學(xué)習(xí)方法較適合處理高維問題，已廣泛應(yīng)用于森林生物量估算，如人工神經(jīng)網(wǎng)絡(luò)、支持向量回歸(support vector regression, SVR)以及隨機森林(random forests, RF)等[13-14]。但需要考慮不同異速生長方程的約束條件，對比分析基于機載LiDAR點云的森林生物量估算模型。

為了更好地理解基于機載激光雷達(dá)數(shù)據(jù)的亞熱帶森林生物量估算模型及其受異速生長方程的影響，本研究旨在比較樣地尺度上的不同異速生長方程的生物量估算差異性，研究基于機載LiDAR數(shù)據(jù)的亞熱帶森林生物量估算模型及其誤差驗證分析。

2 研究區(qū)域與數(shù)據(jù)

2.1 研究區(qū)域概況

研究區(qū)域為福建省龍巖市北部地區(qū)(圖1)，屬亞熱帶海洋性季風(fēng)氣候，以山地為主，地形比較復(fù)雜，境內(nèi)山嶺沿東北—西南走向，大體呈平行分布，自然資源豐富，植被群落類型復(fù)雜多樣，主要喬木樹種有閩楠(Phoebebournei)、絲栗栲(Castanopsisfargesii)、栲樹(CastanopsisfargesiiFranch)、木荷(SchimasuperbaGardn)、楠木(PhoebezhennanS. Lee)、杉木(Cunninghamialanceolata)、馬尾松(PinusmassonianaLamb)等。

圖1 研究區(qū)域范圍

2.2 地面調(diào)查數(shù)據(jù)

本研究區(qū)域包括96個邊長為25.82 m的方形樣地，應(yīng)用中海達(dá)定位儀準(zhǔn)確記錄樣地空間位置。通過2013年現(xiàn)場實地調(diào)查，在每塊樣地中測量胸徑(diameter breast high,DBH)大于5 cm的各類樹木信息，主要包括樹種、胸徑、平均樹高、樹的狀態(tài)(存活、砍伐、枯死等)等，同時獲取了同區(qū)域的239棵解析木數(shù)據(jù)。

2.3 機載LiDAR數(shù)據(jù)

機載LiDAR數(shù)據(jù)于2013年獲取，采用的激光雷達(dá)系統(tǒng)是ALS70-HP，數(shù)碼相機系統(tǒng)是LEICA RCD30，相對航高約200～4 000 m，相對地面飛行速度約250 km/h。該機載LiDAR系統(tǒng)中，激光器工作波長為1 064 nm，脈沖重復(fù)頻率為50 kHz，激光束發(fā)散角為0.15 mrad。獲取的機載LiDAR點云數(shù)據(jù)平均密度約6點/m2，小范圍精細(xì)飛行區(qū)達(dá)9點/m2，地表定位平面精度為0.2 m，高程精度為0.3 m，同步獲取的CCD航空影像數(shù)據(jù)空間分辨率為0.2 m。

3 研究方法

3.1 基于異速生長方程的樹木單體生物量估算

使用4種異速生長方程來計算單木生物量，其來源于周國逸等[12]編譯的生物量方程。根據(jù)該文獻，樹木單體生物量(AGB)計算為4個生物量成分的總和，公式為：

AGB=Btrunk+Bbranch+Bleaf+Broot

(1)

式中：Btrunk表示樹干生物量；Bbranch表示樹枝生物量；Bleaf表示樹葉生物量；Broot表示樹根生物量。

由于外業(yè)采集數(shù)據(jù)只包含樹的胸徑，而不包括單棵樹高信息。因此，只選取該方法中基于胸徑計算生物量的異速生長方程，公式為：

B=a×DBHb

(2)

式中：B為樹木器官生物量(kg)；DBH為胸徑(cm)；a，b為系數(shù)。

選取4種異速生長方程，其中前兩種按樹種類型劃分，后兩種按森林類型劃分。 I：分省(自治區(qū)、直轄市)分樹種，樹種包括杉木、馬尾松、木荷、其他硬闊類4類。II：全國優(yōu)勢種(組)，樹種包括杉木、馬尾松、木荷、其他硬闊類4類。III：分省混合種(組)，森林類型包括針葉林、闊葉林、針闊混交林3類。IV：中國森林主要群系，森林類型包括針葉林、闊葉林、針闊混交林，并在每種森林類型的基礎(chǔ)上按胸徑大小劃分為3種，總共9類。

3.2 基于機載LiDAR點云的森林生物量估算模型

機載LiDAR點云數(shù)據(jù)的地面濾波和處理采用的是TIFFS軟件工具[15]和Terrasolid軟件輔助處理，主要包括：(1)應(yīng)用數(shù)學(xué)形態(tài)學(xué)方法對所有激光雷達(dá)數(shù)據(jù)進行處理，提取窗口數(shù)據(jù)，搜尋窗口內(nèi)的最低點，接著進行開運算和閉運算的處理，將激光點云數(shù)據(jù)分類為地面點和非地面點，進一步運用高程閾值法，從非地面點中分離出建筑、電力線等非植被點，剩余的即為植被點；(2)將地面點進行插值生成數(shù)字高程模型(DEM)，然后根據(jù)植被點的高度坐標(biāo)與相應(yīng)DEM進行歸一化處理，計算每個點的樹冠高度，即得到冠層高度模型；(3)基于冠層高度模型及其三維結(jié)構(gòu)，計算高度均值、方差、峰度、分位數(shù)、偏度、各垂直高度百分比數(shù)等15個特征變量，即生成樣地尺度的生物量預(yù)測變量特征數(shù)據(jù)集[16]。

根據(jù)機載LiDAR點云計算得到的特征數(shù)據(jù)集結(jié)合樣地估算生物量，設(shè)計和構(gòu)建了多種生物量估算模型，包括逐步回歸、支持向量回歸(SVR)、隨機森林(RF)和BP神經(jīng)網(wǎng)絡(luò)。

逐步回歸模型是從一個自變量開始，根據(jù)自變量對因變量作用的顯著程度，從大到小地依次逐個引入回歸方程。當(dāng)引入的自變量由于后面變量的引入而變得不顯著時，要將其剔除掉。引入一個自變量或從回歸方程中剔除一個自變量，為逐步回歸的一步。對每一步均要進行F值檢驗，以確保每次引入新的有統(tǒng)計意義，即方差貢獻顯著的變量前回歸方程中只包含對因變量作用顯著的變量。這個過程反復(fù)進行，直至既無不顯著的變量從回歸方程中剔除，又無顯著變量可引入回歸方程時為止。公式為：

W=α0+α1x1+α2x2+…+αnxn+ε

(3)

式中：W為根據(jù)地面實測數(shù)據(jù)計算的樣地生物量；α0，α1，α2，…，αn為模型系數(shù)；x1，x2，…，xn為基于機載激光雷達(dá)數(shù)據(jù)提取的特征變量；n為參與逐步回歸的特征變量的總數(shù)；ε為正態(tài)分布誤差項[ε～N(0,σ2)]。

機載LiDAR特征變量與生物量之間的關(guān)系通常是非線性的。支持向量回歸(SVR)通過使用核函數(shù)將原始輸入特征空間映射到一個新的高維特征空間[17]，將問題線性化。采用臺灣大學(xué)林智仁博士等開發(fā)的開源程序包LibSVM3.1工具，首先對樣本數(shù)據(jù)進行歸一化處理，將數(shù)據(jù)映射到[0-1]之間，再利用常見的核函數(shù)，通過格網(wǎng)搜索法確定模型參數(shù)。隨機森林RF的回歸過程通常涉及到兩個關(guān)鍵參數(shù)，一個為決策樹數(shù)量，即重抽樣的次數(shù)，另一個為隨機特征數(shù)量，即用來分割節(jié)點的最大輸入變量個數(shù)。首先根據(jù)均方根誤差(RMSE)最小的原則確定決策樹數(shù)量，然后利用袋外數(shù)據(jù)對所有樣本的特征變量的重要性進行計算，最后按特征變量的重要性依次引入模型，選擇預(yù)測能力強、相關(guān)性好的特征變量組合。BP神經(jīng)網(wǎng)絡(luò)將特征變量作為神經(jīng)網(wǎng)絡(luò)輸入變量，樣地生物量作為神經(jīng)網(wǎng)絡(luò)輸出變量，在仿真訓(xùn)練過程中不斷調(diào)試參數(shù)，最終構(gòu)建合適的神經(jīng)網(wǎng)絡(luò)模型進行生物量估算。

為了比較回歸模型的性能，采用五折交叉驗證方法進行模型擬合精度評估，主要指標(biāo)包括確定系數(shù)(R2)和均方根誤差(RMSE)。

4 結(jié)果分析

4.1 基于異速生長方程的生物量估算對比分析

4.1.1單木尺度的生物量估算對比分析

異速生長方程誤差的計算要求對單木進行破壞性測量。而在本研究中，利用解析木數(shù)據(jù)(包含樹干、樹枝、樹葉、樹根及整棵樹的生物量)對異速生長方程加以驗證，表1總結(jié)了4種異速生長方程驗證結(jié)果。

表1 4種異速生長方程估算生物量與實測生物量平均值比較

平均而言，分省分樹種方法的生物量最高，中國森林主要群系方法最低，但4種異速生長方程得到的估算生物量均低于實測生物量。分省分樹種和全國優(yōu)勢種方法在生物量估算上更為準(zhǔn)確，而另外兩種方法，特別是中國森林主要群系方法在生物量估算上嚴(yán)重偏低。這可能與分省分樹種和全國優(yōu)勢種方法在樹種上做了區(qū)分有關(guān)，對各個樹種生物量估算更精確。圖2顯示了解析木數(shù)據(jù)與4種異速生長方程估算生物量之間的關(guān)系。

圖2 4種異速生長方程單木生物量估算與實測值對比

對低生物量的樹木，4種方法估算值都較為接近實測值，但在高生物量的樹木中均出現(xiàn)估算偏低的情況。這可能在建立異速生長方程時，與解析木的選取有關(guān)。在實際野外調(diào)查中，大徑級解析木會比較少，這時由小徑級的解析木數(shù)據(jù)擬合的生物量方程估算大徑級的生物量，可能會導(dǎo)致估算結(jié)果的偏高或者偏低，而在本研究中出現(xiàn)了偏低的情況。

4.1.2樣地尺度的生物量估算對比分析

表2統(tǒng)計了4種異速生長方程在樣地尺度上的生物量密度，樣地數(shù)選擇96個。

表2 4種異速生長方程在樣地尺度上生物量密度

在樣地平均生物量中，分省分樹種方法最高(80.8 t/hm2)，中國森林主要群系方法最低(76.8 t/hm2)，這些平均生物量差別較小，最大差值僅5%。但在樣地最大生物量密度統(tǒng)計上，4種方法存在較大差異，最大差值為23%。對4種異速生長方程估算的樣地生物量密度進行兩兩比較(圖3)。對于低生物量的樣地，任意兩個異速生長方程生物量估算比較接近，而對于大生物量樣地，不同異速生長方程生物量估算有明顯差異。這與Zhao等[9]得出的結(jié)論相似。

圖3 4種不同異速生長方程在樣地尺度上生物量比較

4.2 基于機載LiDAR的森林生物量估算模型

基于5折交叉驗證實驗的所有模型R2和RMSE如表3所示。

表3 基于不同異速生長方程和不同統(tǒng)計方法的激光雷達(dá)生物量模型的交叉驗證統(tǒng)計

總體來說，基于中國森林主要群系的生物量估算模型擬合度最好，后面依次是分省分樹種、分省混合種、全國優(yōu)勢種。把所有生物量估算模型的結(jié)果取平均值，這4種方法對應(yīng)的R2依次為0.614、0.596、0.612和0.678。相對應(yīng)的RMSE分別為37.8 t/hm2、39.9 t/hm2、39.1 t/hm2和32.1 t/hm2。在誤差統(tǒng)計中，RMSE僅用于描述激光雷達(dá)生物量估算模型誤差，其假設(shè)異速生長方程無誤差。然而在4.1中的分析表明，不同異速生長方程的誤差是不同的，特別是中國森林主要群系異速生長方程存在較大誤差，即使其生物量估算模型誤差較小，但其在生物量估算上的誤差也可能很大。

4.3 基于機載LiDAR的森林生物量估算模型差異性分析

在測試的4種不同生物量估算模型統(tǒng)計方法(1種回歸模型，3種機器學(xué)習(xí)方法)中。RF模型擬合結(jié)果最好，BP神經(jīng)網(wǎng)絡(luò)擬合結(jié)果最差。為了比較生物量估算模型預(yù)測生物量與異速生長方程估算生物量之間的關(guān)系，以分省分樹種方法為例進行異速生長方程計算的生物量模型擬合和驗證(圖4)。

圖4 散點圖顯示預(yù)測生物量與參考生物量對比

通過比較散點圖中的模型線發(fā)現(xiàn)，4種生物量估算模型的共同點：大的生物量值往往被低估，而小的生物量值通常被高估。這與Chen[8]發(fā)現(xiàn)有共同之處，RF方法用于生物量回歸時，由于該模型預(yù)測的是訓(xùn)練數(shù)據(jù)值的平均值，會導(dǎo)致低估高生物量區(qū)域和高估低生物量區(qū)域。但不局限于RF方法，其他的線性回歸以及機器學(xué)習(xí)方法均可得出相同的結(jié)論。通過對分省分樹種方法中生物量大于100 t/hm2的32個樣地的計算，估算生物量的最大值和平均值分別為261 t/hm2、145 t/hm2，而SR、SVR、BP神經(jīng)網(wǎng)絡(luò)以及RF預(yù)測生物量的最大值依次為214 t/hm2、166 t/hm2、238 t/hm2、183 t/hm2，相應(yīng)的平均值依次為125 t/hm2、116 t/hm2、129 t/hm2、131 t/hm2，在最大值和平均值水平上，預(yù)測生物量均低于估算生物量。高生物量樣地生物量被低估是一個普遍現(xiàn)象，Gao等[18]也在利用生物量估算模型對高生物量樣地進行生物量估算時發(fā)現(xiàn)，多種模型出現(xiàn)生物量低估的情況，特別是RF和k-近鄰法模型對生物量嚴(yán)重低估，導(dǎo)致無法正確預(yù)測生物量。

研究結(jié)果表明，在亞熱帶研究區(qū)域，BP神經(jīng)網(wǎng)絡(luò)模型的回歸性能結(jié)果較差(相關(guān)系數(shù)比SR模型低0.175)，隨即森林RF的結(jié)果較好，和Oliveira等[19]研究結(jié)果一致。但與Chen[8]和Li等[20]研究發(fā)現(xiàn)的RF效果較差形成對比。由此得出，不同森林類型、樹種組成以及使用的異速生長方程等不同，會導(dǎo)致同樣的生物量模型的結(jié)果差異。另外，研究結(jié)果揭示了不同移速生長方程及其對應(yīng)的機載LiDAR點云生物量估算模型，二者在生物量估算中的不一致性。進而，在森林生物量估算誤差的定量評估中，應(yīng)使用誤差傳播理論[21]全面考慮整個過程的誤差，包括野外測量誤差、異速生長方程誤差、遙感數(shù)據(jù)誤差、生物量估算模型誤差等[22-24]。

5 結(jié)論

選取福建省龍巖市的亞熱帶森林為研究對象，構(gòu)建多種異速生長方程和生物量估算模型，利用機載LiDAR數(shù)據(jù)和樣地調(diào)查數(shù)據(jù)對樣地上的生物量進行監(jiān)測，結(jié)果表明：在異速生長方程的選擇上，分省分樹種的方法在樣地生物量估算上更為準(zhǔn)確，中國森林主要群系的方法效果最差；在基于機載激光雷達(dá)的森林生物量估算模型中，隨機森林的方法精度最高，BP神經(jīng)網(wǎng)絡(luò)的方法精度最低。

對生物量估算模型的直接比較(如R2)而不考慮異速生長方程誤差的差異是不合適的，未來需要更多地研究森林生物量估算過程中的誤差，了解相關(guān)誤差的組成及其傳播方式。