基于機載激光雷達的國外松片林蓄積量反演研究

馬 驤，姜文龍，閆保銀，金 悅，崔 立

(1.常州市金壇區林業技術指導站，江蘇 金壇 213200； 2.南京國圖信息產業有限公司，江蘇 南京 210037)

森林是地球生態系統的重要組成部分，在給人類帶來林副產品產生經濟效益的同時，還具備涵養水源、調節區域氣候條件，實現碳平衡、碳循環等功能，產生生態效益。在IPCC發布的《全球溫控1.5 ℃特別報告》中明確指出，“林業是實現1.5 ℃溫控目標的重要路徑之一，通過林業固碳除碳等技術，可實現溫室氣體的大幅減少，在促進森林生態系統恢復，提高基于生態系統和社區的適應性以及濕地管理方面潛力巨大。”在此背景下，如何快速準確地獲取森林資源信息，實現動態監測，為森林保護、森林經營等提供依據，成為當前森林資源管理的重要任務之一。

激光雷達(LiDAR)的出現為森林資源調查方式的轉變提供契機。作為一種主動遙感技術，具有精度高、成本低等優勢，早在20世紀80年代已用于林業研究[1]。1984年，Nelson等[2]通過研究激光雷達與樹木冠層郁閉度關系，指出兩者間相關性顯著；MacLean等[3]通過研究樹冠垂直剖面與森林蓄積量相關性，指出2者的對數呈線性相關。進入21世紀后，越來越多學者也相繼證明激光雷達對于估測樹高、胸高斷面積、生物量等森林參數準確性較高[4-6]。Morsdorf等[7]采用林木冠幅分割算法對瑞士國家森林公園針葉林進行分割，再生成種子點進行聚類分析，獲得樹冠直徑后按圓面積公式反演了林木胸徑；劉清旺等[8-9]基于機載激光雷達技術，使用多元線性回歸、隨機森林模型等方法，對森林平均高、生物量、郁閉度等參數進行提取分析，其反演結果較光學遙感而言精度更高；Zhao等[10]采用2種尺度不變模型估測生物量，結果表明該2種方法R2范圍在0.80—0.95之間，可靠性較高。本文以江蘇省常州市金壇區碳匯計量山區監測點中的國外松片林作為研究對象，進行機載激光雷達蓄積量反演試驗，評價這一技術在江蘇省蘇南地區國外松林的適用性，以及在蓄積量反演方面的可能性，為今后森林資源調查新技術應用提供參考。

1 研究區概況與數據獲取

1.1 研究區概況

研究區位于江蘇省常州市金壇區4 km×4 km碳匯計量山區監測點范圍內，國家2000大地坐標系(China Geodetic Coordinate System 2000, CGCS2000)范圍坐標：西北至(20 718 118.46，3 517 999.08)，西南至(20 718 118.46，3 513 999.07)，東南至(20 722 118.48，3 513 999.07)，東北至(2 0722 118.48，3 517 999.08)。區內地貌以丘陵為主，為北亞熱帶季風區，年均氣溫15.3 ℃四季分明，雨量充沛，年降水量1 063.5 mm。日照充足，日照率46%，無霜期228 d，年平均濕度為78%。

研究區面積為16 km2，根據2018年碳匯計量監測成果，區內森林總蓄積量為48 486.04 m3，分布樹種主要有濕地松(PinuselliottiiEngelmann)、杉木(Cunninghamialanceolata)、香樟(Cinnamomumcamphora)、女貞(Ligustrumlucidum)、柳杉(Cryptomeriafortunei)等。本次研究對象位于4 km×4 km碳匯計量山區監測點(見圖1)，77，250號國外松小班范圍內。其中77號小班用于激光雷達森林參數反演，250號小班用于反演數據驗證。

圖1 試驗區范圍

1.2 LiDAR數據獲取

分別于2021年4月16日、5月12日，利用激光雷達在4 km×4 km碳匯計量山區監測點，77，250號小班范圍內，進行機載激光雷達進行試驗數據采集。使用大疆經緯M300 RTK作為飛行平臺，飛行航速為6 m/s，總覆蓋面積約為10 hm2，飛行時天氣晴朗，少量高云分布，對于LiDAR飛行影響較少。

在本次飛行中，激光雷達系統采用大疆禪思L1(DJIL1)，集成Livox激光雷達模塊、高精度慣導、測繪相機、三軸云臺等模塊。飛行時設置航高80 m，旁向重疊度50%。鏡頭角度-90°，激光掃描方式為三回波重復掃描，類似于傳統線掃激光雷達，能獲取更均勻、精度更高的掃描結果，也更穩定。

1.3 地面樣地數據獲取

2021年5月對77，250號小班進行外業調查。

本次地面實測調查分為大樣地與小樣地。大樣地以77號小班范圍為界，在其中按照20 m×20 m劃分小樣地，共計34塊。在樣地范圍內記錄樹種、胸徑、樹高等參數，所得數據將以小樣地為單位進行林分參數反演。

在250號小班內，按照《森林資源規劃設計調查技術規程》[11]中小班測樹因子調查的相關要求，深入小班內部，選擇具有代表性的調查點，以調查點為中心布設樣地并每木檢尺，記錄樹種、胸徑、樹高、郁閉度等測樹因子，按照《江蘇省主要樹種一元材積表》計算公頃蓄積，所得數據用于林分參數反演的數據可靠性驗證。

圖2 樣地設置

2 研究方法

2.1 樣地數據處理

根據樣地內每木檢尺獲得的胸徑，對照《江蘇省主要樹種一元材積表》，查詢樹種為“馬1”一列，以及胸徑對應的徑階，填寫單株立木的材積，最后再統一求和，獲得樣地內總蓄積量。計算公式為

M樹種=V樹種/ ΣAi；

M全小班=∑(M樹種×A)。

式中，A：小班面積；Ai:第i塊小樣地面積；V樹種：某樹種各樣方蓄積之和；M樹種：某樹種小班每公頃蓄積；M全小班：全小班蓄積。

對于77號小班，其范圍內的34塊小樣地需逐塊計算平均胸徑，用于后期反演模型的構建，250號小班需求得小班平均胸徑、每公頃蓄積量，用于反演結果驗證。

2.2 LiDAR數據處理

在林分反演模型構建之前，首先需將點云去噪，其次對點云進行地面點分類，以此分離地面點與地物點，然后分別生成數字地面模型(DTM)、數字表面模型(DSM)和冠層高度模型(CHM)，并基于CHM生成種子點，在點云歸一化處理(NPC)，消除地形影響后，最后進行單木分割。本次對于LiDAR數據的處理均在LiDAR 360軟件中完成。

地面點分類是點云數據處理的基礎操作，本文采用改進的漸進加密三角網濾波算法(IPTD)[12]分類地面點。其原理是：對于不含建筑物的點云，以默認值作為格網大小，取格網內最低點作為起始種子點，構建初始三角網。遍歷所有待分類點后，查詢各點水平投影所落入的三角形，并計算點到三角形的距離d及點到三角形3個頂點與三角形所在平面所成角度的最大值，分別與迭代距離與迭代角度進行比較，如果小于對應閾值，則將此點判定為地面點，并加入三角網中。重復此過程，直至所有地面點分類完畢。

地面點分類后，需經插值處理生成DSM，DTM。LiDAR 360軟件中提供了反距離權重插值(IDW)、克里金插值(Kriging)和不規則三角網插值(TIN)方法。本文采用TIN插值法，分別生成DSM和DTM，兩者作差后生成CHM。利用生成的DTM高程對點云作歸一化處理，消除地形影響。設置地面高度、最小樹高再進行單木分割。

2.3 基于單木分割的反演模型構建

胸徑是計算林木蓄積量的重要參數之一，若能準確獲取胸徑值，則可根據地方相應的一元材積表求算蓄積量。

樹高、冠幅等參數可通過機載激光雷達直接測得，而胸徑則需通過建立回歸方程的方式間接獲取[13]。本文擬參考劉清旺[14]、何祺勝等[15]的研究方法，以34塊小樣地為單位，統計各樣地內LiDAR估測的單木平均樹高、平均冠幅等參數，作為自變量，實測平均胸徑作為因變量，建立實測胸徑與LiDAR估測參數的一元回歸方程，間接估測胸徑。共設置3組：

(1)實測胸徑與LiDAR樹高

(2)實測胸徑與LiDAR冠幅

(3)實測胸徑的自然對數與LiDAR樹高的自然對數

一元回歸方程建立后，根據相關系數，確定最優估測方程。

2.4 反演模型可靠性評價

反演模型采用擬合優度(R2)、均方根誤差(RMSE)和相對均方根誤差(rRMSE)進行可靠性評價。

R2表示模型的擬合優度，一般R2值越高，表明該模型的擬合度就越高。計算公式如下：

RMSE表示實際值與預測值之間的均方根,誤差越小，則方程擬合效果越好，計算公式為

rRMSE為相對均方根誤差，其值應當越小越好，計算公式為

3 結果與分析

3.1 樣地實測數據分析

77號小班的LiDAR數據經數據處理后共獲得單木309株，與實測株數381株相比單木識別率達81%。在小班內按20 m×20 m大小劃分小樣地34塊，以小樣地為單位，分別統計LiDAR估測數據(含樹高、冠幅、冠幅面積和冠幅體積)和實測胸徑的平均數、中位數、最大值、最小值以及標準差，結果如表1所示。

表1 34塊小樣地林分因子統計

在LiDAR估測數據中，平均樹高為18.7 m，其最大值為21.1 m，最小值為14.2 m，標準差達到1.5；平均冠幅為7.1 m，其最大值為12.8 m，最小值為3.9 m，標準差2.0；平均冠幅面積與體積標準差分別為24.6和128.3，數據離散程度較高。

在實測胸徑方面，平均胸徑為29.4 cm，最大值為33.8 cm，最小值為20 cm，標準差為2.6。

3.2 胸徑反演模型構建

建模前需對34塊小樣地數據進行樣本異常值診斷和處理。經分析后發現，第16號小樣地的LiDAR估測平均樹高較低，其原因是由于該樣地范圍內林下分布植被，導致單木分割時錯將林下植被作為國外松進行分割，降低了樣地平均樹高值。因此需重新設置最小樹高，將林下植被進行過濾，修正平均樹高后進行胸徑反演。

基于LiDAR估測與樣地實測數據結果，按照基于單木分割的反演模型構建中的方法，設置3組反演模型：

(1)實測胸徑與LiDAR樹高

(2)實測胸徑與LiDAR冠幅

(3)實測胸徑的自然對數與LiDAR樹高的自然對數

反演結果如圖3、表2所示。

表2 實測胸徑與估測胸徑回歸方程

在第(1)組模型中，以LiDAR估測樹高作為自變量，實測胸徑作為因變量建立反演模型。結果顯示：獲得方程y=26.7 lnx-48.7，擬合優度R2=0.741 1，均方根誤差(RMSE)=0.004。

在第(2)組模型中，以LiDAR估測冠幅作為自變量，實測胸徑作為因變量建立反演模型。結果顯示：獲得方程y=- 0.233 9x2+ 3.142 3x+ 19.797，擬合優度R2=0.603 1，均方根誤差(RMSE)=0.008。

在第(3)組模型中，首先將LiDAR估測樹高與實測胸徑分別進行自然對數轉換，再以LiDAR估測樹高的自然對數作為自變量，實測胸徑的自然對數為因變量，建立反演模型。結果顯示：獲得方程y=- 2.184 9x2+13.569x-17.605，擬合優度R2=0.806 3，均方根誤差(RMSE)=0.001。

圖3 3組胸徑反演模型

綜上所述，在3組反演模型中，分別將LiDAR估測樹高和實測胸徑取自然對數后，再建立反演模型的效果最佳，此時擬合優度(R2)達到最高，均方根誤差(RMSE)最小。因此，結合本文樣地實測數據分析，當LiDAR估測樹高范圍在14—22 m時，推薦使用第3組方程y=- 2.184 9x2+13.569x-17.605作為胸徑的最優估測方程。

3.3 LiDAR估測蓄積量驗證

250號小班面積為3.269 6 hm2，在其內部選擇2個具有代表性的點位布設樣地，并每木檢尺，測得該小班平均胸徑為25.8 cm，平均樹高12.2 m，公頃蓄積量為183.6 m3/hm2。

250號小班經LiDAR飛行后，獲得估測平均樹高為12.7 m；將該小班單木估測樹高代入胸徑反演模型y=- 2.184 9x2+13.569x-17.605進行胸徑反算，得LiDAR估測平均胸徑為26.6 cm。查《江蘇省主要樹種一元材積表》中“馬1”列，根據LiDAR估測的單木胸徑統計對應徑階的材積，得估測公頃蓄積量為157.2 m3/hm2。

根據《森林資源規劃設計調查技術規程(GB/T 26424-2010)》“13 調查精度”中允許誤差的要求，主要小班調查因子允許誤差分經營單位性質、小班森林類別等分為A，B，C 3個等級。250號小班位于金壇區薛埠鎮東進村金沙灣高爾夫球場內，屬金壇區重點林區，因此小班調查因子允許誤差適用等級“A”。

經LiDAR估測數據與實測數據對比(如表3)，平均胸徑誤差為3.1%，平均樹高誤差為4.1%，公頃蓄積量誤差為14.4%，均在《森林資源規劃設計調查技術規程(GB/T 26424-2010)》中允許誤差范圍內(平均胸徑、平均樹高、每公頃蓄積量誤差需控制在5%，5%，15%范圍以內)，符合精度要求。

表3 LiDAR估測與實測數據對比

4 結論與討論

4.1 LiDAR估測受多因素制約

LiDAR技術較傳統的人工調查而言可節約大量人力物力，提升工作效率[16]。然而，LiDAR技術在實際應用時易受多因素制約，影響精度。制約因素主要包括：LiDAR采樣密度、航帶匹配誤差和數據處理方法[17]。

在本研究中，選擇的250，77號小班均為國外松片林，周圍無建筑物遮擋，從單木結構角度來說確保了樹冠上部的回波信號，減少了信號的滅失，同時點密度為100個/m2，較好地反映了試驗區現狀；在航帶匹配方面，對LiDAR飛行小班均進行了航帶平差處理，保證了定位精度；在數據處理方面，由于試驗區地貌均為丘陵，地形趨于平緩，因此在地面點分離方面采用改進的漸進加密三角網濾波算法(IPTD)，將點云分為地面點與植被點，然后分別生成DTM，DSM，CHM和NPC，對點云進行單木分割。總體而言本次LiDAR估測盡量減少了試驗誤差，估測精度較高。

4.2 LiDAR估測蓄積量可能性探討

一般而言，林分蓄積量與生物量隨平均胸徑的增加而增加[18]。研究表明，LiDAR估測樹高與實測樹高關系密切，而樹高也與胸徑、蓄積量有較高的相關性[19-20]。

本研究通過建立3組胸徑反演模型，結合擬合優度(R2)挑選出優選方程，反算了250號小班的林木胸徑。最后查《江蘇省主要樹種一元材積表》，得公頃蓄積量為157.2 m3/hm2。雖然其結果符合精度要求，但與實測公頃蓄積量183.6 m3/hm2仍有較大差距。分析原因，筆者認為部分區域樹冠較為茂密，同時林下植被較為復雜，單木分割時對于樹冠邊緣難以區分，最終導致獲得的單木數量較實際偏少，估測公頃蓄積量較實測數據偏小。

4.3 LiDAR技術應用綜合展望

本次采用LiDAR技術對金壇區國外松片林小班進行了估測，并獲得了良好的效果，對于森林資源調查領域新技術應用提供了技術參考。

目前對于機載激光雷達用于森林資源調查更偏向于理論研究，實際應用較少。其主要原因在于LiDAR飛行時客觀因素制約較多，影響數據的誤差，而且飛行成本偏高，若進行區縣范圍內的全域飛行，將投入大量物力。此外，飛行效果視林分結構而定，若林分結構、樹種結構、林下植被較為簡單，則LiDAR估測精度較高，若林分結構復雜區域，其精度則會降低[16]。

綜上所述，LiDAR技術是未來發展的大趨勢，前景應用廣闊。在精度符合要求的情況下，可在樣地水平成果的基礎上大范圍應用，節約人力物力，提高工作效率。