聯合LiDAR和多光譜數據森林地上生物量反演研究

巨一琳，姬永杰，黃繼茂，張王菲

(1.西南林業大學林學院，云南 昆明 650224;2.西南林業大學地理與生態旅游學院，云南 昆明 650224)

森林是陸地生態系統的主體，在生物圈的生物地球化學過程中發揮重要的“緩沖器”和“調節器”的功能，對改善生態環境、維護生態平衡有著重要的作用[1]。森林地上生物量(above ground biomass, AGB)作為森林生產力的重要評價指標，是評估森林碳收支的重要參數，是系統發揮其他生態功能的物質基礎，也是陸地生態系統進行碳循環的重要內容[2-3]。森林地上生物量的傳統調查方法是通過現地測量森林的胸徑及樹高等因子，使用異速生長方程或蓄積量-生物量模型計算森林地上生物量?；谫Y源清查的人工調查方法可準確估算森林地上生物量，但其需要花費大量時間、人力和物力，還會對森林造成一定的干擾與損傷，且只適于較小區域的森林AGB估算。不同于傳統資源調查方法，使用遙感技術立足不同傳感器獲取影像，可快速準確、連續、動態、長期地估算森林AGB，同時可在大區域范圍內實現森林AGB的估算?；谶b感技術的森林AGB反演有光學遙感、激光雷達、微波遙感等方式，在不造成森林生態損傷的前提下可極大提高對森林AGB反演的效率和精度[4]。

多光譜遙感具有豐富的光譜特征及植被指數，但其穿透性差，難以獲取森林垂直結構特征，并且估算森林生物量時容易出現易飽和、敏感性差的現象[5-6]。激光雷達(light delection and ranging, LiDAR)具有穿透植被葉冠，并能直接獲取目標的三維坐標，在地形、林木高度和植被空間結構的探測方面具有極大優越性[7]。將多光譜遙感與LiDAR數據相結合，進行森林AGB反演是當前研究的熱點之一。Popescu等[8]融合機載LiDAR數據和多光譜數據，結合地面調查數據建立森林生物量反演模型，結果顯示兩種數據的融合比單一數據反演AGB精度更高。徐婷等[9]借助LiDAR和Landsat8 OLI多光譜數據采用線性多元逐步回歸法分別建立森林生物量反演模型，結果表明兩種數據建立的綜合模型反演精度更高。卜帆[10]結合機載LiDAR和高光譜數據進行生物量反演，結果表明兩種數據結合時采用支持向量機建立森林生物量反演模型精度較高，決定系數(R2)達到0.7。

利用遙感數據進行森林AGB反演通常采用參數模型和非參數模型兩類方法，多元線性逐步回歸為參數模型方法中反演森林AGB既常用又有效的算法。非參數模型中的KNN-FIFS(k-nearest neighbor with fast iterative features selection, KNN-FIFS)算法是韓宗濤[11]2017年基于馬氏距離K-NN算法，提出以迭代方式優選遙感特征，并將特征因子優化組合的反演模型。韓宗濤等[12]以Landsat8 OLI和SAR數據為數據源，采用KNN-FIFS算法對大興安嶺根河研究區森林AGB進行估測，結果表明KNN-FIFS算法精度顯著優于多元線性逐步回歸法(SMLR)及K-NN算法，R2達到0.77。張少偉等[13]采用KNN-FIFS算法結合2009年和2014年兩期主、被動遙感數據，進行特征組合及快速優化構建估測模型，結果表明KNN-FIFS算法可大幅度提升高維度、多模式遙感特征優選效率。

本研究以LiDAR數據和Landsat8 OLI數據為數據源，結合地面實測數據，進行數據源的不同組合，并篩選優化特征，再以多元線性回歸及KNN-FIFS兩種算法建模，對研究區典型的寒溫帶落葉松林進行森林AGB反演，最后，對各自模型的適宜性及精度及準確性進行對比分析，從而為推動大興安嶺地區典型森林AGB動態監測提供參考。

1 材料與方法

1.1 研究區概況

研究區位于內蒙古根河市大興安嶺森林生態系統國家野外科學觀測研究站(121°30′～121°31′E，50°49′～50°51′N)，面積為102 km2，高程范圍為810～1 116 m，為我國目前緯度最高的森林生態系統野外科學觀測站。研究區屬寒溫帶大陸季風氣候，年均氣溫-5.3 ℃，氣溫年較差47.4 ℃，為高緯度多年凍土、寒溫帶森林生態系統的典型區域。該生態站森林植被覆蓋度達75%以上，優勢樹種有興安落葉松(Larixgmelinii)、樟子松(Pinussylvestrisvar.mongolica)、闊葉樹白樺(Betulaplatyphylla)和山楊(Populusdavidiana)等。研究區地理位置見圖1。

1.2 研究區數據獲取

1.2.1 機載LiDAR數據獲取

本研究使用的機載LiDAR數據以“運-5”為機載平臺，載有Leica機載雷達系統，獲取時間為2012年8—9月，平均飛行高度2 700 m，共32個航帶，激光脈沖發射頻率為100～200 kHz，掃射角度為±35°。傳感器記錄了首次和末次激光返回信息，記錄點云數據的三維坐標信息、高程值、強度信息以及返回類型等信息[14]。

1.2.2 Landsat8 OLI數據獲取

本研究使用的多光譜數據為Landsat8 OLI數據中的B1—B7波段，從地理空間數據云官方平臺獲取?？紤]到2013年8—9月研究區光學遙感影像云層覆蓋比較嚴重，另選取成像時間為2013年10月19日、影像無云層覆蓋、航帶號為123/24的一景數據，具體參數見表1。

表1 Landsat8 OLI陸地成像儀參數

1.2.3 樣地調查數據獲取

研究區樣地調查數據分別獲取于2012年37塊30 m×30 m固定樣地和2013年18塊45 m×45 m臨時樣地，共55塊。樣地位置采用差分GPS定位，精度控制在1 m以內。樣地進行每木調查，因子為胸徑、樹高、枝下高、冠幅、相對坐標，同時記錄樣地林分特征，灌木、草本的種類以及高度等參數。樣地優勢樹種包括興安落葉松和白樺兩種，樣地內兩類樹單木生物量的計算采用式(1)和式(2)得到，樣地生物量采用陳傳國等[15]的冪指數異數生長方程式(3)得到。

(1)

(2)

W樣地=∑W1+W2+……+Wn。

(3)

式中：W為森林生物量，DBH為實測胸徑。

1.3 遙感數據處理

1.3.1 LiDAR數據預處理

LiDAR數據處理前期部分由數據供應商完成轉換、拼接及質量檢測；后期部分主要借助于數字綠土LiDAR360軟件。主要包括以下步驟：點云數據定義坐標系、點云去噪、提取地面點、提取植被點、提取DEM以及點云歸一化。歸一化后得到的冠層高度即去除地形高度影響得到的真實樹高，最后將點云特征提取到樣地尺度上。

1.3.2 Landsat8 OLI數據預處理

本研究獲取的Landsat8 OLI數據已完成地形的幾何校正，因此在預處理影像時不用做幾何校正。使用ENVI5.3軟件進行輻射定標和大氣校正，通過輻射定標可以實現DN值與輻射亮度、反射率值等物理量之間的轉化；大氣校正消除大氣和光照等因素對地物反射的影響從而獲取其真實的地表反射率。

1.4 遙感影像特征提取

機載LiDAR第1次回波信息魯棒性較強，將第1次脈沖回波形成的參數應用在森林AGB反演結果中效果較好[4,16]。

Landsat8 OLI影像提取4類變量因子，分別為單波段因子、植被指數、主成分分析因子、紋理因子。單波段因子包括B1—B7，植被指數包括歸一化植被指數(INDVI)、差值植被指數(IDVI)、比值植被指數(IRVI)、土壤調節植被指數(ISAVI)、增強植被指數(IEVI)、有效葉面積指數(ISLAVI)[4]、地表反照率(RAlbedo)、B4/RAlbedo、B547、B65、B74等變量。主成分分析因子提取了前3個波段(PCa1—PCa3)；紋理因子提取了B1—B7波段的8種紋理信息，包括均值(Me)、方差(variance，Var)、對比度(contrast，Con)、信息熵(entropy，En)、協同性(homogeneity、Homo)、相異性(dissimilarity,Dis)、二階矩(second order moment，Sec)、相關性(correlation，Cor)[17]。LiDAR和Landsat8 OLI變量含義及公式見表2。

表2 LiDAR和Landsat8 OLI特征變量匯總

1.5 變量篩選和建模分析

1.5.1 偏最小二乘法及變量投影重要性分析

偏最小二乘回歸是對多元線性回歸模型的一種擴展，是在普通多元回歸的基礎上增加主成分分析、典型相關分析的思想[17], 很好地解決了只用一個線性模型來描述獨立變量Y與預測變量組X之間的關系。

變量投影重要性(VIP)是基于偏最小二乘法基礎上進行變量篩選的過程，用來判斷單個自變量在解釋因變量的重要性，如果自變量對因變量的重要性越大，則VIP值越大；如果各自變量對y的解釋作用一樣，所有的VIP值均為1；有學者建議用0.8作為臨界值來區分重要變量與不重要變量[18]，其計算公式為：

(4)

式中：k為自變量個數，ch為相關自變量提取的主成分，r(y,ch)為因變量和主成分的相關系數，表示主成分對y的解釋能力，whj為自變量在主成分上的權重。

1.5.2 線性多元逐步回歸

線性多元逐步回歸是一種線性多元回歸模型進行變量篩選的算法[19]，其回歸模型為：

Y=β0+β1x1+β2x2+……+βmxm+ε。

(5)

式中：Y是生物量的值；x1、x2……xm是預測變量，β0為常數項，β1、β2……βm為回歸方程的系數；ε是隨機誤差值。

1.5.3 KNN-FIFS方法

本研究采用韓宗濤提出的KNN-FIFS方法，是在K-NN方法的基礎上進行特征因子優化組合，提高森林參數反演效率[11]。KNN-FIFS組合優化基本原理如下：

1)設樣地數為n，特征數為m，由樣地數據和遙感特征提取訓練數據，即樣地對應特征所在像元的值；

2)初始化最優特征子集為空集，因此最優模型均方根誤差理論上為最大值；基于K-NN法, 依次利用特征建立森林AGB反演模型，得到最優特征子集個數的K-NN反演模型及每個模型對應的均方根誤差(RMSE)；

3)得到的最優RMSE, 即RMSE最小值, 設該值為σ(RMSEb), 將研究區RMSE最大值設為RMSE0[σ(RMSE)]。若σ(RMSEb)<σ(RMSE0)則將σ(RMSEb)賦給σ(RMSE0), 并將RMSEb對應的特征子集賦給特征子集，反之迭代結束。

KNN-FIFS算法在迭代運行中，當距離度量標準確定時,K值將影響反演結果，而K值則受樣地信息以及參數等因素影響[12]。其中K值默認為1～11，光譜特征信息提取窗口為1～11。

1.5.4 精度驗證評價方法

反演精度評價采用決定系數(R2)和均方根誤差[RMSE，式中記為σ(RMSE)]，以及相對均方根誤差[RMSEr，式中記為σ(RMSEr)]來檢驗，R2越接近1，代表著模型精度越高；精度驗證采用留一交叉驗證方法。RMSE值越低，說明回歸模型更準確[20]。計算公式如下：

(6)

(7)

(8)

2 結果與分析

本研究先采用偏最小二乘法對不同的遙感數據源組合方式提取的參數進行信息重要性排序，再依據其重要性排序確定最優參數組合，最后使用線性多元逐步回歸和KNN-FIFS算法構建森林AGB反演模型。

2.1 遙感影像參數優選

采用偏最小二乘法回歸算法進行參數優選，依據變量投影的重要性(VIP)做變量篩選，選擇0.8以上的變量后續分析建模。將LiDAR和Landsat8 OLI所提取的參數作為自變量輸入SIMCA 14.1中進行偏最小二乘法回歸算法的建模和VIP值排序，篩選VIP>0.8的影像參數進行后續森林AGB反演模型構建。單一LiDAR遙感數據時篩選得出18個參數；單一Landsat8 OLI遙感數據時篩選得出38個參數；LiDAR、Landsat8 OLI兩種遙感數據組合時篩選得出33個參數(圖2)。

2.2 森林AGB反演模型構建

將在單一LiDAR、單一Landsat8 OLI及LiDAR & Landsat8 OLI聯合遙感數據3種組合形式下的優選參數及地面實測樣地數據，分別作為自變量和因變量輸入線性多元逐步回歸和KNN-FIFS模型，使用R2和RMSE分析各模型反演的結果，并分析各模型精度和適宜性。

2.2.1 線性多元逐步回歸模型建立

在多源遙感數據分別是單一LiDAR、單一Landat8 OLI、LiDAR & Landat8 OLI 3種組合時，將使用偏最小二乘法篩選而得的18、38、33個參數輸入線性多元逐步回歸模型，最后進一步篩選得出以下特征和形成最優模型。模型結果及最優參數見表3和圖3。

表3 多元逐步回歸參數方程

2.2.2 KNN-FIFS模型建立

在兩種遙感數據源及3種組合方式下，將以偏最小二乘算法篩選出的參數輸入KNN-FIFS模型。KNN-FIFS算法的最佳模型參數見表4，進一步特征篩選結果見圖4。在單一Landat8 OLI、LiDAR & Landat8 OLI兩種數據組合方式下，均篩選得出B65參數,說明該參數對于森林AGB相關性較大且較為敏感。這是因為植被在紅光波段有著很強的吸收特性,同時在近紅外波段有著很強的反射特性，故選擇近紅外和紅光波段進行植被遙感監測效果最佳。另在兩種數據組合最后篩選參數中分別包含了冠層密度和高度分位參數，說明這兩種體現森林水平和垂直結構因子的參數，在反演森林AGB中發揮著關鍵性的作用。

表4 KNN-FIFS特征組合

2.3 森林地上生物量反演結果對比

由模型評價結果(表5)可知模型反演精度如下：①單一LiDAR數據的線性多元逐步回歸模型和KNN-FIFS模型反演精度相差不大，R2、RMSE以及RMSEr僅相差0.02、2.05 t/hm2和1.57%，總體上線性多元回歸模型的精度更好，對于數據的適宜性更好。②單一Landsat8 OLI數據的線性多元逐步回歸模型和KNN-FIFS模型反演精度相差較大，R2、RMSE以及RMSEr已相差0.36、9.64 t/hm2以及15.42%，KNN-FIFS模型明顯具有較高的精度及對Landsat8 OLI數據的適宜性。③LiDAR & Landsat8 OLI數據組合的兩種模型反演精度均較高，分別為線性多元逐步回歸模型和KNN-FIFS模型在各種數據組合方式的最高。兩種模型的反演精度也相差較小，R2、RMSE和RMSEr相差只有0.04、2.99 t/hm2和2.36%。

表5 模型評價結果

通過以上反演結果分析表明：①采用LiDAR和Landsat8 OLI兩種數據相結合，構建線性多元逐步回歸和KNN-FIFS模型均可達到最高的反演精度，機載激光雷達和多光譜遙感數據的多源數據結合可發揮各自數據優勢，能表征森林水平和垂直結構，最大程度反演森林的三維結構，進而提高森林AGB反演精度和準確度。②采用參數優選技術可有效提高森林AGB反演精度。首先采用偏最小二乘算法對兩種數據源3種組合方式時總體參數篩選，后又在兩種模型構建的同時進行了特征優選，兩次的參數優選過程去除了和森林AGB相關性低、敏感性差的參數。③在單一LiDAR和LiDAR&Landsat8 OLI兩種組合下，針對該研究區典型寒溫帶落葉松林森林AGB，使用兩種模型的反演精度幾乎相差不大，說明模型在該數據組合及森林類型模型適宜性相似。而在單一Landsat8 OLI時，KNN-FIFS反演精度明顯高于線性多元回歸模型，KNN-FIFS模型的適宜性更強。

由散點圖(圖5)分析可知，兩種模型在兩種數據源的3種組合方式下的估測值均出現不同程度的高估和低估現象，其中單一LiDAR數據源時，兩種模型反演的結果均出現明顯的低估現象，而LiDAR & Landsat8 OLI組合下高估與低估現象卻不明顯。關于估測值與實測值的線性擬合性上，森林AGB值約80 t/hm2以下時，兩種模型在兩種數據3種組合方式下總體上均表現出較好的線性擬合性，實測值與估測值之間相關性較高；但隨著森林AGB數值的增加，擬合性表現得較為離散。這是因為機載LiDAR對森林的穿透能力，能夠提供準確的樹高和直觀的冠層垂直結構信息，結合反應森林水平結構的光譜信息，表現得在單一機載LiDAR數據和兩種數據結合反演森林AGB中精度為較高，線性擬合性較好。而在單一Landsat8 OLI數據時，整體擬合性均低于LiDAR擬合結果，尤其在線性多元逐步回歸模型反演時表現最差，低于KNN-FIFS反演效果。這是因為在KNN-FIFS擬合方法中，估測參數只與相鄰的k個樣本有關，因此樣地點的加權求值可以減少樣本不平衡問題，能更好地描述各參數與影像之間的非線性關系[12]。對于整體擬合結果較差，這與多光譜數據具有豐富的光譜特征及植被指數，但穿透性差，難以獲取對森林垂直結構特征參數，對森林AGB反演具有較大的結構性局限性有關。

2.4 森林地上生物量空間分布

基于3種數據源采用線性多元逐步回歸算法和KNN-FIFS算法進行單一反演以及聯合反演得到航飛區森林地上生物量制圖(圖6)。

從圖 6a—6c可以看出，基于多元逐步回歸算法的LiDAR數據的制圖結果中，生物量分布<40、≥40～70以及≥160～200 t/hm2的森林區域各占30%，其余占10%；采用Landsat8 OLI數據的制圖結果中生物量≥40～70 t/hm2的區域約占研究區的40%，≥100～130 t/hm2的區域約占研究區的40%，<40 t/hm2的區域約占研究區的20%；采用融合LiDAR & Landsat8 OLI數據的制圖結果中生物量≥40～70 t/hm2的區域約占研究區的50%，≥130～160 t/hm2的區域約占研究區的30%，<40 t/hm2的區域約占研究區的20%。另外結合實際土地利用現狀(圖6b)可發現Landsat8 OLI數據的森林地上生物量結果中的高值區出現在了林區主干道附近，有明顯的錯估現象。

從圖6d—6f中可以看出，基于KNN-FIFS算法的LiDAR數據的制圖結果中，生物量分布<40 t/hm2、≥40～70 t/hm2的區域各占研究區的40%，其余占20%；采用Landsat8 OLI數據的制圖結果中生物量≥40～70 t/hm2的區域占50%，<40 t/hm2的區域約占研究區的20%，≥130～160 t/hm2約占30%；采用融合數據的制圖結果中生物量≥40～70 t/hm2約占研究區的50%，<40 t/hm2約占研究區的30%，其余占20%。結合地面調查信息來看，基于KNN-FIFS算法反演出的生物量制圖與調查信息基本相符，<40 t/hm2的區域和≥40～70 t/hm2的區域占比較大。另外結合根河生態站航飛區森林地上生物量高值區主要分布在航飛區北部以及山脊線附近，低值區主要分布在海拔較低的林區主干道。航飛區域森林地上生物量的空間分布趨勢與實際的地貌特征相符。航飛區北部處于山坡的陰坡處，可以較好地存留水分，因此森林地上生物量較高；在航飛區域的林區主干道附近，由于修建主干道這種工程可以改變土層結構以及條件，因此類似的人為干擾破壞植被活動會導致森林地上生物量偏低。

3 討 論

1)森林地上生物量的監測對區域生態植被保護及生態系統平衡都有著重要的意義，相比采用傳統人工地面調查森林AGB方法，使用遙感技術可快速獲取大區域尺度、實時及經濟的森林AGB分布情況。本研究以LiDAR和Landsat8 OLI影像為數據源結合地面實測數據，對不同組合形式下的遙感參數進行篩選，并采用線性多元逐步回歸和KNN-FIFS算法，對航飛區森林AGB進行反演及分析研究。

研究結果中采用單一LiDAR數據的線性多元回歸和KNN-FIFS模型反演結果的R2為0.76和0.74。龐勇等[21]將機載LiDAR數據劃分為針葉林、闊葉林和針闊混交林3個類型進行森林AGB反演，決定系數可達0.82以上，優于筆者的建模結果。本研究采用的機載LiDAR數據雖能準確獲取森林垂直結構信息，但仍受真實地形起伏的影響，使歸一化的點云高度產生誤差，又未考慮森林類型的影響，在一定程度上影響了本研究在單一LiDAR數據情況下的精度。

Landsat8 OLI數據的兩種模型的反演精度R2為0.24和0.60，反演精度結果相比LiDAR數據結果較低。Foody等[22]研究表明體現光學特征的光譜信息在不同的地物之間存在較大差異，地形存在起伏較大區域或郁閉度較低區域，光譜信息受到地形等影響較大導致不能準確表述植被特征，使得森林AGB數值產生10～15 t/hm2的反演誤差。羅洪斌等[4]同樣結合Landsat8 OLI和LiDAR兩種數據采用偏最小二乘法對橡膠林地上生物量進行遙感反演中，采用單一Landsat8 OLI數據建模精度為0.56，同樣比單一LiDAR數據時建模精度低。上述學者在采用單一Landsat8 OLI數據進行森林AGB反演時精度與本研究結果大致相當，多光譜數據反演精度低的原因多與光譜信息不能反映森林垂直結構、受地形影響較大及采用數據分辨率低等密切相關。

本研究融合LiDAR & Landsat8 OLI兩種數據源采用兩種模型反演森林AGB，反演精度R2分別0.84和0.80。胡凱龍等[23]利用光學數據的紋理特征和機載激光雷達的點云特征進行森林地上生物量反演分析，在不區分森林類型的情況下，建模精度低于本研究的建模結果，分別為0.73和0.79。徐婷等[9]基于機載LiDAR和光學數據進行建模，分析在不同森林類型的情況下估算綜合森林地上生物量，得到的R2可以達到0.88，證實了將典型森林類型劃分為不同森林類型會提高建模精度。本研究在未考慮森林類型的情況下，聯合兩種數據源建?？傮w上提高了反演的精度，這說明LiDAR數據和多光譜數據聯合、參數組及特征優選對我國寒溫帶森林AGB反演精度提高是行之有效的，另采用兩種模型也是適宜的。

2)本研究通過對LiDAR數據和Landsat8 OLI數據提取出的特征變量，采用偏最小二乘法進行特征優選，再通過線性多元逐步回歸算法和KNN-FIFS算法反演森林地上生物量形成結論如下：LiDAR和Landsat8 OLI兩種數據源的融合相比單一數據源對森林AGB的反演精度更高；在兩種數據3種組合方式下進行特征的優選，采用偏最小二乘法算法進行篩選對反演精度的提高是有效的；線性多元回歸及KNN-FIFS兩個模型在不同組合下，對本研究區的森林類型反演是適宜的。

在未來的研究中可將光譜分辨率高的高光譜數據納入，加大不同數據的優勢互補，深刻反饋森林三維結構的細化特征；還可嘗試參數及非參數其他模型的算法研究，尋找最優模型。另在加大特征提高的同時，應深入其他篩選算法，去除冗余特征降維提取相關性最大特征；對于熱帶森林等不同森林類型也應加大樹種、模型適宜性研究及精度和飽和點提高的研究。

根據本研究結果看，采用偏最小二乘法進行特征變量的篩選，可以為后續的森林地上生物量建模提高精度。從特征變量VIP值排序上看，LiDAR百分位高度變量特征在排序中較靠前，反映了其在森林地上生物量有一定優勢。此外，在紋理因子的提取中，可以考慮設置窗口的大小以便更好地進行特征選擇。