基于機載全波形LiDAR數據的森林地上生物量估測算法研究

邢艷秋，姚松濤，李夢穎，謝杰，閆燦

(東北林業大學 森林作業與環境研究中心，哈爾濱 150040)

基于機載全波形LiDAR數據的森林地上生物量估測算法研究

邢艷秋，姚松濤，李夢穎，謝杰，閆燦

(東北林業大學 森林作業與環境研究中心，哈爾濱 150040)

為提高森林地上生物量(AGB)的估測精度，本研究以白樺林為研究對象，以機載全波形激光雷達(LiDAR)數據為研究數據，首先提出了機載全波形LiDAR數據讀取與全波形特征信息提取的相關算法，然后結合具體算法的實現提取出每條全波形數據對應的各波形分量的能量信息，進而依據波形能量信息計算出每個樣地的全波形激光穿透指數(FiLPI)，之后通過全波形激光穿透指數FiLPI建立其與對應樣地實測森林AGB的統計回歸模型，同時將森林AGB的估測值與森林AGB的實測值進行對比。結果表明全波形激光穿透指數FiLPI與森林AGB具有很好的相關性(R2為0.885，RMSE為0.095)，并且森林AGB的估測值與實測值之間誤差的波動較小，提高了森林AGB的估測精度，彌補和提供了機載全波形LiDAR數據估測森林AGB的方法和思路。

機載激光雷達；全波形數據；波形數據讀取；波形特征信息提取；白樺林；地上生物量

0 引言

在林業調查中，森林地上生物量(Aboveground Biomass，AGB)通常指的是單位面積內植被的絕干物質量[1-2]。傳統森林AGB的測定方法耗時耗力耗材且更新較慢。激光雷達(Light Detection And Ranging，LiDAR)作為一種主動式遙感技術，廣泛應用于森林經營管理與生態系統研究等領域[3]。按回波模式，LiDAR系統通常可分為離散回波和全波形兩種，其中，離散回波LiDAR通過記錄冠層的單個或多個離散回波信號以獲得不同密度的點云數據；全波形LiDAR則通過記錄返回脈沖的全部能量來得到亞米級森林垂直剖面。LiDAR遙感技術估測森林AGB是非破壞性的，有利于生態環境的保護。

大多研究學者基于離散回波LiDAR點云數據對森林AGB進行了大量的估測研究，Bortolot等[4]提出了一種適用于不同林型和不同區域的森林AGB估測的可視化算法，該算法以LiDAR點云數據提取的林木參數和外業實測數據為訓練樣本同AGB建立回歸關系，進而估測其它區域的森林AGB，結果表明森林AGB的外業實測值同估測值的相關性在0.59～0.82之間，能夠成功估測不同區域的森林AGB。劉清旺等[5]研究了離散回波LiDAR點云數據估測單株樹AGB，由點云數據生成的冠層高度模型結合優化的單株樹冠特征識別算法估測相關的結構參數，結果表明，由LIDAR點云數據得到的單株樹估測參數與實測參數顯著相關(R2=0.729)，可以估測單株樹AGB。王濤等[6]研究表明，森林AGB同葉面積指數等參數相關，然而有關機載全波形LiDAR數據估測森林AGB的研究相對較少。

為研究機載全波形LiDAR數據在森林AGB估測研究中的應用潛能，本研究提出了全波形數據讀取和波形特征信息提取算法，進而利用波形能量信息參數進行森林AGB估測研究，彌補機載全波形LiDAR數據估測森林AGB的方法和思路。

1 研究區與數據獲取

1.1 研究區概況

本研究選擇內蒙古自治區額爾古納市東南部的依根農林交錯區作為研究區域，地理位置坐標為120°36′50.48″E～120°52′56.53″E，50°21′11.08″N～ 50°24′32.00″N。植物資源以天然次生白樺林(Betulaplatyphylla)為主，混生樹種包括落葉松(Larixgmelinii)、樟子松(Pinussylvestrisvar.mongolica)等，林下灌木層主要由石棒繡線菊(Spiraeamedia)、筐柳(Salixlinearistipularis)等組成。

1.2 外業樣地數據

本研究于2012年8月至9月對研究區天然次生白樺林進行了相應的地面調查。綜合考慮林地的森林類型及生長狀況，一共布設了80個邊長為10 m的方形樣地，測量了樹高、胸徑、冠幅等單木參數。采用Vertex IV超聲波測高測距儀測量了每木樹高，30 m范圍內其高度分辨率為0.1 m；采用胸徑尺測量每木胸徑(1.30 m樹高處周長)，并利用皮尺測量東西和南北方向的單木冠幅。根據研究需求，借助Trimble GeoXT6000手持GPS定位儀對樣地中心點進行定位，定位精度達50 cm。外業樣地數據統計結果見表1。

表1 外業樣地數據基本統計量

1.3 機載全波形LiDAR數據

本研究機載全波形LiDAR數據采集時間為2012年8月26日，天氣狀況良好，采用的飛機型號為運-5，相對飛行高度1300 m，搭載的LiDAR掃描系統為Leica ALS60，采用的大地坐標系統為WGS84，投影坐標系統為UTM。Leica ALS60系統最高掃描頻率為200 kHz，最大視場角為75°。每個激光脈沖在地面上所形成的光斑直徑大約為0.2～0.3 m，全波形數據是返回脈沖被接收望遠鏡接收后由采樣器對其進行數字化采樣而成，記錄了激光光斑內目標物對應的完整波形數據[7-8]。

2 數據處理與研究方法

2.1 樣地森林地上生物量計算

王正文等[9]對位于內蒙古自治區呼倫貝爾市額爾古納右旗境內白樺林生物量進行了評估，并確定了白樺林AGB相對生長模型如公式(1)所示，模型相關系數達到0.985。由于本研究區與前者地理位置相鄰且樹種相同，因此采用公式(1)計算樣地內單木AGB。

(1)

則樣地AGB平均值為公式(2)所示：

(2)

式中：Wi表示第i株單木生物量，kg；Di表示第i株單木胸徑，cm；Hi表示第i株單木樹高，m；Wm表示樣地平均地上生物量，kg；n表示樣地內單木株數。

2.2 全波形數據讀取

全波形文件的讀取是分析處理波形數據的首要前提。本研究采集的全波形數據經過坐標轉換等處理后，以LAS1.3格式存儲在文件中，LAS1.3增加了對全波形數據的支持，以二進制形式存儲，記錄了接收到的完整波形數據，包括采樣間隔和實時記錄的方位等。通過訪問一個激光點數據，能夠以映射字段訪問到與該點關聯的返回脈沖波形，讀取到對應的波形信息。計算機系統訪問內存速度比訪問硬盤速度快得多，在讀取數據時如果少訪問硬盤而更多地訪問內存，那么讀取速度將會快速提升。基于這種思路，本研究使用IDL編程語言借助內存映射方法[10]實現了研究區LeicaALS60LAS1.3全波形數據的快速讀取。

2.3 分析處理與波形特征信息提取

本研究將原始全波形數據近似地看作不同目標物后向散射高斯波形的疊加，是噪聲與一系列高斯分量之和在時間軸上的疊加，即原始全波形數據是等時間間隔采樣的序列{tm，f(tm)}，其中{tm，m=1，…，n}表示采樣時間，{f(tm)，m=1，…，n}表示采樣時間tm對應的振幅值，符合混合高斯分布。使用公式(3)來表示由一系列簡單高斯分量組成的原始波形數據。

(3)

2.3.1 噪聲去除

每條完整全波形數據對應著19.2 m的沿脈沖距離，在實際工作環境中若不考慮異常噪聲點，則可穿透目標物的分布范圍通常要小于這個長度[11]。據此可選某個四分位數振幅值作為去噪閾值，從而剔除整個波形數據中的噪聲振幅值。四分位數方法適合于具有脈沖特征信號數據的平滑去噪，四分位數位置的計算如公式(4)所示。

(4)

2.3.2 波形分量個數確定

全波形數據高斯分解的準確性和效率性與初始參數的選取有很大關系。本研究依據獲得的數據以離散點數據的形式給出，提出了離散型函數求解拐點。根據拐點特性，拐點前后曲線的上升(或下降)趨勢不一樣，若點P為該曲線的拐點，則其肯定滿足公式(5)。

(kAB-kBP)(kPC-kPD)<0。

(5)

式中：A、B表示點P前任意兩點，C、D表示點P后任意兩點，kAB，kBP，kPC，kCD分別表示直線AB、直線BP、直線PC、直線CD的斜率。該算法尤其適用于離散點數據橫坐標為遞增等間距的情況，采集的LeicaALS60全波形數據本身是遞增等間隔采樣的數據。所以在離散點數據遞增等間隔采集時，若點(tm，f(tm))為波形的拐點，則滿足的判斷條件由公式(5)變換為了簡化公式(6)。

G(m)=(f(tm)+f(tm-2)-f(tm-1))(f(tm+2)+f(tm)-2f(tm+1))<0。

(6)

式中：m=3，4，5，…，N-3，N表示全波形數據采樣點個數，通過全波形數據的讀取結果可知N為128。全波形數據可能存在多個拐點，當算法開始時，需要將所有的序列點{tm，f(tm)}依次代入公式(6)，計算出G(m)，當G(m)小于0時算法即可終止并存儲記錄m時的點(tm，f(tm))。由于高斯分布關于直線t=μ對稱，基于此將高斯分布的兩個拐點分別稱為左拐點和右拐點，自然地，一個高斯分布由其一個左拐點(tleft=μ-σ，f(tleft))及其一個右拐點(tright=μ-σ，f(tright))唯一確定，進而確定該波形數據中的高斯分量個數[12-13]。

2.3.3 高斯分量初始特征參數估計

為解決實際波形數據中的左右拐點劃分，本研究在去除背景噪聲后將檢測出的拐點分別與其右邊第一個點連成直線，如果直線斜率大于零，則該拐點為左拐點；如果直線斜率小于零，則該拐點為右拐點。從而可將判斷出的所有拐點分為左、右拐點。

小學時期正是學生成長的重要時期。在這個時期，學生的身心都在發展過程中，是對于學生進行思想培養、促進思維成長的重要階段。在過去的教學當中，課堂上的主體是教師，學生知識被動的進行學習，學生更多是死記硬背，使得學生的思想固化，阻礙了學生多種思維的發展，以至于影響了學生未來的成就。因此在現今的小學教學當中，應該提出學生的主體性地位，培養學生多方面能力，促進學生的發展。

(7)

在確定了一個高斯分量的左拐點與其右拐點后，可以反推出其位置tm和半波寬σm的初始估計值，同時，將其左拐點和右拐點之間點的最大振幅值作為其峰值Am的初始估計值，它們的初始估計值如公式(7)所示。據此可以得到每個高斯分量的初始參數(Am、tm、σm)，從而完成該波形數據高斯分量初始特征參數估計步驟。

2.3.4 波形能量信息計算

高斯分量的能量信息是其波形擬合曲線與采樣時間軸所形成的面積，因而其計算方法如公式(8)所示。

(8)

式中：Em表示第m個高斯分量的能量值；Am表示通過波形擬合計算出的第m個高斯分量的波峰振幅值；tm表示通過波形擬合計算出的第m個高斯分量的位置；σm表示通過波形擬合計算出的第m個高斯分量的標準差(半寬)。

2.4 全波形數據激光穿透指數

本研究采用統計模型借助全波形激光穿透指數LPIfi對森林地上生物量進行分析和驗證。全波形數據激光穿透指數反映了激光脈沖的穿透能力，由以上內容得知，本研究第i塊樣地內的全波形激光穿透指數FiLPI恰好如公式(9)所示。

(9)

式中：i=1，2…80，N表示樣地內個數；Eci和Egi分別表示第i塊樣地內的全波形激光脈沖中的冠層返回脈沖能量和地面返回脈沖能量，需要通過波形能量信息計算公式(8)計算得到。

從外業采集的80個樣地數據中隨機選取60個作為建模樣本，余下20個作為檢驗樣本分別進行森林地上生物量回歸模型的建模和驗證。建模精度用決定系數R2評價，其值越大，說明模型構建結果越好；模型預測精度用RMSE(Root Mean Square Error，RMSE)進行評價，其值越小，說明模型預測效果越好。

3 結果與分析

3.1 樣地生物量計算

根據外業樣地實測數據按照2.1中的方法得到樣地AGB統計結果見表2。

表2 樣地生物量統計結果

3.2 全波形數據讀取結果與顯示

本研究按照上述方法借助IDL8.3平臺實現了研究區機載全波形LiDAR數據的讀取，并將其以波形的方式顯示出來，結果如圖1所示(某一條全波形數據)。此外，不僅順利讀取了全波形數據，還同時獲取了與全波形數據相關的參數，如采樣的起始時間、結束時間、采樣間隔、采樣個數和坐標位置參數等，這些數據為后續的波形分析處理與波形特征信息提取提供基礎參數信息。

圖1 編譯程序讀取并顯示的某條全波形數據Fig.1 A waveform record of this compiled program

由結果可知波形采樣數為128，相較于傳統離散回波LiDAR，全波形LiDAR不僅能記錄目標物的三維坐標信息[14-15]，還能存儲整個返回波形。由于讀取全波形數據文件需要專業軟件，本研究基于采集的LAS1.3機載全波形LiDAR數據，實現了全波形數據的快速讀取及可視化，為后續的波形分析處理與波形特征信息提取提供基礎。

3.3 全波形數據分析處理與波形特征信息提取結果

本研究根據后續工作的需要，通過上述闡述的方法對研究區樣地內的Leica ALS60 LAS1.3格式的全波形數據進行了處理。共有80個邊長為10 m的方形樣地，記錄了約 8.2 萬條的全波形數據，其中大約99.93%的都能得到很好處理，其余的根據局部最大值法進行了處理。某條原始全波形數據的各步驟處理結果如圖2～4所示。

圖2 某條原始全波形數據的四分位數法噪聲估計Fig.2 Noise estimation of the quartiles of the raw full-waveform data

圖3 拐點計算結果與原始全波形數據的比較Fig.3 Comparison between result of inflection point calculation and the raw full-waveform data

圖4 高斯分量波形擬合結果與原始全波形數據的比較Fig.4 Comparison between component waveform fitting and the raw full-waveform data

在全波形掃描器的情況下，波形分解還可提供脈沖寬度和脈沖強度作為屬性，其中脈沖強度是返回信號的積分，因此該強度值物理上是脈沖能量。相比之下，傳統離散回波LiDAR系統僅提供返回信號的三維坐標，并且在許多情況下也提供強度，然而一般情況下對該強度值是如何在這些系統中得到的知之甚少，通常它是返回脈沖的最大振幅值，其被錯誤地稱為強度[14]。

通過波形數據的高斯分解，可以提取每個波形分量的振幅、寬度、位置等重要的波形特征信息，進而計算每個波形分量的能量信息，即圖5中所示的起波點與止波點之間圍城的面積，且其計算方法如公式(8)所示。

圖5 Leica ALS60全波形數據及波形特征信息Fig.5 Leica ALS60 full-waveform data and waveform parameters

3.4 模型構建結果與精度分析

經統計分析得知，通過機載全波形LiDAR數據求出的全波形激光穿透指數FiLPI與實測AGB有較好的相關性(R2=0.885，RMSE=0.095)，模型構建結果如圖6所示，建模結果比較理想且實測AGB與全波形激光穿透指數FiLPI具有明顯的負相關關系：AGB實測值隨著FiLPI的增大而減小，即當森林地上生物量越大或森林越茂密時激光脈沖的穿透能力反而越小，這與森林樣地實際情況基本吻合。

圖6 FiLPI與AGB實測值的統計回歸關系Fig.6 Statistical Regression Relation between FiLPI and Measured AGB

通過驗證樣本進行預測，將計算出的AGB的估測值與AGB的實測值進行對比，可以直觀地得到估測值與實測值之間的誤差，如圖7所示，可知AGB實測值與AGB估測值之間存在一定誤差，但出現較大誤差的情況極少。

通過計算得到驗證樣地的RMSE為21.12 kg，且誤差波動較小，這說明該方法能夠有效地進行森林AGB的估測。因此，基于本研究的機載全波形LiDAR數據的相應算法與研究方法可以準確地估測森林AGB，彌補了全波形數據估測森林AGB的方法和思路。

圖7 實測AGB與預測AGB的預測誤差Fig.7 The prediction error of the field-measured AGB against the predicted AGB

4 結論與討論

本研究基于機載LiDAR全波形數據，通過全波形數據相關處理算法的實現，能夠準確地提取每個波形分量以及相應的波形特征信息，在實際處理中表現出較好的適應性。此外在每條全波形數據處理過程中，存儲記錄了其中每個波形分量的波峰振幅值、位置和波寬等波形特征信息(參數)，尤其是首末次波形分量的，并計算存儲對應的能量信息。進而提取出全波形激光穿透指數FiLPI估測參數，并通過統計模型估測了森林樣地AGB，表明本文基于機載全波形LiDAR數據的相應算法研究可以準確地估測森林AGB，彌補了全波形數據估測森林AGB的方法和思路。

由于本研究的林型為闊葉林，樹種為白樺，在林型和樹種的多樣性方面比較單一。因此，在下一步工作中，希望將本研究中的算法應用到不同類型的林區，以使機載全波形LiDAR數據森林AGB估測算法研究能夠得到推廣應用驗證，以提高相關算法研究及統計模型的理論性和可推廣性。

[1]尤號田，邢艷秋，王萌，等.小光斑激光雷達數據估測森林生物量研究進展[J].森林工程，2014，30(3)：39-42.

[2]李妍，徐新良，張超.中國喬木林碳儲量變化研究[J].森林工程，2015，31(4)：50-55.

[3]李增元，劉清旺，龐勇.激光雷達森林參數反演研究進展[J].遙感學報，2016，20(5)：1138-1150.

[4]Bortolot Z J，Wynne R H.Estimating Forest Biomass Using Small Footprint LiDAR Data：An Individual Tree-Based Approach that Incorporates Training Data[J].Isprs Journal of Photogrammetry & Remote Sensing，2005，59(6)：342-360.

[5]劉清旺，李增元，陳爾學，等.機載LIDAR點云數據估測單株木生物量[J].高技術通訊，2010，20(7)：765-770.

[6]王濤，龔建華，張利輝，等.基于機載激光雷達點云數據提取林木參數方法研究[J].測繪科學，2010，35(6)：47-49.

[7]王素元，馬洪超，王杰棟，等.基于分組LM算法的全波形LiDAR高斯分解[J].測繪與空間地理信息，2016(7)：144-147.

[8]邱賽，邢艷秋，田靜，等.星載LiDAR與HJ-1A/HSI高光譜數據聯合估測區域森林冠層高度[J].林業科學，2016，52(5)：142-149.[9]王正文，王德利，臧傳來，等.大興安嶺次生林白樺對林下日陰菅及其它主要草本植物的影響[J].生態學報，2001，21(8)：1301-1307.

[10]楊群.在VB利用內存映射文件讀取大數據量原始遙感影像圖[J].計算機工程與應用，2003，39(6)：135-137.

[11]盧昊，龐勇，徐光彩，等.機載激光雷達全波形數據與系統點云差異的定量分析[J].武漢大學學報信息科學版，2015，40(5)：588-593.

[12]段乙好，張愛武，劉詔，等.一種用于機載LiDAR波形數據高斯分解的高斯拐點匹配法[C]//全國地理信息科學博士生學術論壇，2014-11-21，2014：189-197.

[13]Mallet C，Bretar F.Full-Waveform Topographic Lidar：State-of-the-Art[J].Isprs Journal of Photogrammetry & Remote Sensing，2009，64(1)：1-16.

[14]Reitberger J，Schn?rr C，Krzystek P，et al.3D Segmentation of Single Trees Exploiting Full Waveform LIDAR data[J].Isprs Journal of Photogrammetry & Remote Sensing，2009，64(6)：561-574.

[15]Wagner W，Ullrich A，Ducic V，et al.Gaussian Decomposition and Calibration of a Novel Small-Footprint Full-Waveform Digitising Airborne Laser Scanner[J].Isprs Journal of Photogrammetry & Remote Sensing，2006，60(2)：100-112.

Estimation Algorithm of Forest Aboveground Biomass Based on Airborne Full-waveform LiDAR Data

Xing Yanqiu，Yao Songtao，Li Mengying，Xie Jie，Yan Can

(Center for Forest Operational Environment，Northeast Forestry University，Harbin 150040)

In order to improve the estimation accuracy of forest aboveground biomass(AGB)，this paper takes birch forest as the research object and airborne full-waveform light detection and ranging(LiDAR)data as the research data.First，the algorithms of LiDAR data reading and full-waveform information extraction were proposed.Then the waveform energy corresponding to each sample of full-waveform data was extracted and the full-waveform laser penetration indexFiLPIwas calculated.After that，the statistical regression model of AGB was established by using the indexFiLPI.Finally，the estimated value of forest AGB was compared with the measured value.Results showed that the indexFiLPIhas a good correlation with forest AGB(R2is 0.885，RMSEis 0.095)，and the error between the estimated value and the measured value of forest AGB was small，which improved the estimation accuracy of forest AGB and provided the method and idea of forest AGB estimated by airborne full-waveform LiDAR data.

Airborne LiDAR；full-waveform data；waveform data reading；waveform feature information extraction；birch forest；aboveground biomass

2016-04-18

國家林業局林業公益性行業科技專項基金(201504319)；國家重點基礎研究發展計劃項目(2013CB733404)

邢艷秋，博士，教授。研究方向：森工管理與林業信息工程。E-mail：yanqiuxing@nefu.edu.cn

邢艷秋，姚松濤，李夢穎，等.基于機載全波形LiDAR數據的森林地上生物量估測算法研究[J].森林工程，2017，33(4)：21-26.

S 758.5

A

1001-005X(2017)04-0021-06