基于機載LiDAR的偽大光斑波形數(shù)據(jù)反演林分平均高

閆 燦，邢艷秋，高金萍，辛 穎，田 昕

（1.東北林業(yè)大學 森林作業(yè)與環(huán)境研究中心，黑龍江 哈爾濱 150040；2.國家林業(yè)與草原局調(diào)查規(guī)劃設計院， 北京 100714；3.黑龍江省測繪科學研究所，黑龍江 哈爾濱 150040）

準確及時地估算森林結構參數(shù)對于森林資源管理以及碳循環(huán)和生物多樣性研究具有重要意義[1-2]。其中，林分平均高指樣方內(nèi)所有林木的加權平均樹高，作為估計其他參數(shù)的基礎，對森林生物量/儲蓄量估算及動態(tài)變化的研究至關重要[3]。利用多種遙感技術手段獲取林分平均高可以有效進行森林資源估測[4-7]，而激光雷達技術（Light detection and ranging,LiDAR）在林業(yè)研究中得到了廣泛應用。

激光雷達系統(tǒng)可以獲取相互獨立的點云數(shù)據(jù)和波形數(shù)據(jù)，波形數(shù)據(jù)是未經(jīng)處理的原始回波序列，點云數(shù)據(jù)由硬件系統(tǒng)直接產(chǎn)生[8]。當前，已有很多研究基于小光斑激光雷達（直徑0.2～0.3 m）實現(xiàn)了林分平均高的提取，Qin 等[9]利用小光斑激光雷達獲取的多次回波點云數(shù)據(jù)估算冠層高度；穆喜云等[10]基于冠層高度模型（Canopy height model,CHM）使用四分位數(shù)法估計林分平均高；焦義濤[11]通過設定植被高度平均值閾值的方式估測林分平均高。分析發(fā)現(xiàn)，基于小光斑離散點云數(shù)據(jù)可以有效估測森林林分平均高，但點云數(shù)據(jù)在估測樹高方面存在著一些問題，如：點云分類算法、點云密度等，都會限制估測精度的提高。而波形數(shù)據(jù)相較于點云可以提高垂直結構信息含量，但小光斑由于光斑直徑小的原因，造成了樹頂無法被準確識別的問題，進而會影響樹高的估計[12]。

針對上述問題，本研究提出以小光斑機載LiDAR 波形數(shù)據(jù)為基礎，在波形數(shù)據(jù)高斯分解預處理的基礎上，將高斯分量在時間軸上按脈沖能量比例以加權平均的方法進行疊加，最終形成一條類似大光斑LiDAR 波形的數(shù)據(jù)，即偽大光斑波形數(shù)據(jù)，旨在降低小光斑LiDAR 因光斑直徑過小帶來的影響，增加冠頂有效信息的獲取，更加客觀地反映樣地內(nèi)整體情況。在此基礎上，分別基于小光斑LiDAR 波形數(shù)據(jù)和偽大光斑LiDAR 波形數(shù)據(jù)提取特征參數(shù)，并結合野外樣方實測平均樹高建立LiDAR 林分平均樹高反演模型，進而進行比較分析。

1 研究區(qū)概況與數(shù)據(jù)來源

1.1 研究區(qū)概況

研究區(qū)為位于內(nèi)蒙古自治區(qū)呼倫貝爾市額爾古納上庫力農(nóng)場的依根農(nóng)林交錯區(qū)，地理坐標 位 置 為120°36′50.48″～120°52′56.53″ E，50°21′11.08″～50°24′32.00″ N（圖1）。研究區(qū)地貌主要以山丘地形為主，氣候為寒溫帶大陸性季風氣候，境內(nèi)主要河流為根河，主要生長著天然次生白樺林。

圖1 研究區(qū)樣地布設方案Fig.1 Sample layout scheme of study area

1.2 激光雷達數(shù)據(jù)

所用研究區(qū)數(shù)據(jù)采集于2012年8月26日，由運5 飛機搭載Leica ALS60 激光雷達掃描系統(tǒng)，坐標系采用UTM 投影坐標系和WGS84 大地坐標系，飛機飛行高度約為1 300 m。此次飛行獲取的東西方向23 個條帶數(shù)據(jù)覆蓋整個研究區(qū)，激光雷達數(shù)據(jù)采集過程中飛行軌跡如圖2所示。

圖2 研究區(qū)飛行軌跡示意Fig.2 Diagram of flight trajectories

1.3 野外數(shù)據(jù)采集

為了使地面實測數(shù)據(jù)與激光雷達數(shù)據(jù)保持一致性，對研究區(qū)林分于2012年8―9月進行了野外調(diào)查。此次實地調(diào)查共設置了10 m×10 m 的方形樣地80 個，借助Trimble GeoXT6000 GPS 對樣地的中心點進行定位，定位精度達50 cm，結合Vertex IV 超聲波測高測距儀、胸徑尺和皮尺分別測量每木樹高、枝下高、胸徑及冠幅。野外數(shù)據(jù)測量結果如表1所示。

表1 野外數(shù)據(jù)基本統(tǒng)計量Table1 Basic statistics of measured data

基于Lorey's 平均樹高方法計算林分平均樹高，公式如式（1）所示。

式（1）中，hi為第i棵樹的樹高（m），hL為林分的平均樹高（m），n為樣地樹木的株數(shù)，gi為第i棵樹木的胸高斷面積（m2）。

2 研究方法

2.1 激光雷達數(shù)據(jù)預處理

原始波形數(shù)據(jù)讀取后采用高斯分解方法進行預處理，主要包括原始波形平滑去噪、高斯分量拐點計算及左右拐點確定、初始特征參量估計和高斯分量的波形擬合過程。本研究通過四分位數(shù)方法進行波形去噪，即設定原始波形振幅的某一四分位數(shù)為閾值以剔除噪聲。拐點是依據(jù)高斯函數(shù)二階導數(shù)的情況來判斷。在確定左右拐點后，根據(jù)公式計算高斯分量的初始特征參數(shù)，包括位置、半波寬和振幅。波形擬合則采用非線性最小二乘法。

2.2 偽大光斑波形擬合

2.2.1 權數(shù)獲得

利用DIGSILENT依次搭建了代表①、②、③、④四類負荷分布狀況的典型單電源輻射狀配網(wǎng)拓撲模型如圖1至圖4所示。

預處理得到每條原始波形數(shù)據(jù)對應的高斯分量波形與高斯函數(shù)十分接近，其中，每一高斯分量波形與對應采樣時間的積分即為脈沖強度值（圖3），圖中I2為返回信號的積分，物理意義為脈沖能量。

圖3 全波形中的三維點與屬性Fig.3 3D points and attributes derived from a waveform

以高斯函數(shù)為模型進行模擬，高斯函數(shù)的積分值對應曲線橫軸上一定區(qū)間的面積，即脈沖能量值，其計算公式如式（2）所示。

式（2）中，A、t、σ分別為高斯分量中對應特征參數(shù)振幅、位置和半波寬，由波形預處理得到。

由高斯函數(shù)積分公式可知，高斯分量的面積與振幅和半波寬之積成正比，由此計算樣地內(nèi)各高斯分量脈沖能量占總脈沖能量的比例，并將其視為各高斯分量特征參數(shù)對應權數(shù)，計算公式如式（3）所示。

式（3）中，θi為樣地內(nèi)第i條高斯分量波形脈沖能量在總脈沖能量中的權數(shù)，Ai為樣地內(nèi)第i條高斯分量的振幅值，σi為樣地內(nèi)第i條高斯分量的半波寬，n為樣地內(nèi)回波波形分量的個數(shù)。

2.2.2 偽大光斑高斯函數(shù)模型的建立

在得到樣方內(nèi)各高斯分量波形脈沖能量占總脈沖能量權數(shù)的基礎上，分別計算疊合后高斯函數(shù)的位置、振幅及半波寬，形成偽大光斑波形對應高斯函數(shù)模型。計算公式如式（4）所示。

式（4）中，n為樣地內(nèi)回波波形分量的個數(shù)，f(x)表示偽大光斑波形對應高斯函數(shù)模型。

2.3 參數(shù)提取

表2 回波波形參數(shù)說明Table2 Description of echo waveform parameters

圖4 大光斑波形參數(shù)示意Fig.4 Large footprint waveform parameter diagram

2.4 LiDAR 林分平均樹高反演模型建立

森林垂直結構參數(shù)與波形數(shù)據(jù)關系密切，本研究提取波形長度（Extent）、波形后緣長度（TraEdg）和波形前緣長度（LeaEdg）等與平均樹高相關的回波波形參數(shù)作為自變量，以野外實測平均樹高H作為因變量擬建立平均樹高反演模型。綜合考慮研究區(qū)地形與森林類型特點，從野外采集的80 個樣方中隨機選取60 個地形平坦的提取提取波形參數(shù)，并將剩余20 個樣方用于模型精度評價。繪制波形參數(shù)與實測平均樹高的散點圖，觀察其相關性，進而確定平均樹高回歸模型。

觀察圖5各參數(shù)與實測平均樹高散點圖可以發(fā)現(xiàn)波形后緣長度（TraEdg）和波形前緣長度（LeaEdg）與實測平均樹高相關性很弱，波形長度（Extent）與實測平均樹高線性關系較為明顯。大光斑波形數(shù)據(jù)進行平均樹高反演時，常用波形長度進行估算，本研究結合波形長度（Extent）建立林分平均樹高反演模型如式（5）。

式（5）中：a為參數(shù)對應系數(shù)，b為常量，H為野外實測Lorey's 平均樹高，Extent 為波形長度。

2.5 模型評價指標

林分平均樹高估測模型擬合相關性采用決定系數(shù)R2來評價，其值越大，說明建模精度越好。模型預測精度用P評價，值越大模型預測精度越好，計算方式如式（6）。

式（6）中，P為預測精度，H為野外實測Lorey's平均樹高，為林分平均樹高反演模型估算值。

圖5 實測平均樹高與波形參數(shù)散點Fig.5 Measured average tree height and waveform parameters scatter plots

3 結果與分析

3.1 波形數(shù)據(jù)高斯分解

原始波形數(shù)據(jù)經(jīng)高斯分解預處理后，結果如圖6所示，其中圖6(a)為原始波形數(shù)據(jù)經(jīng)去噪后檢測出的拐點位置及左右拐點的判斷，圖6(b)為高斯分量的最小二乘波形擬合結果。

圖6所示預處理結果表明，研究區(qū)機載LiDAR 全波形數(shù)據(jù)經(jīng)高斯分解處理流程得到很好的預處理效果。

圖6 機載LiDAR 全波形數(shù)據(jù)預處理Fig.6 Preprocessing of airborne LiDAR full-waveform data

3.2 偽大光斑波形擬合

基于高斯分解方法對樣地內(nèi)所有波形進行分解處理，并將分解后所有高斯分量繪制出來，如圖7(a)所示為一個10 m×10 m 樣地內(nèi)的全部波形高斯分量。然后將高斯分量在時間軸上按脈沖能量比例以加權平均的方法進行疊加，擬合形成一條偽大光斑波形數(shù)據(jù)，即實現(xiàn)了樣地內(nèi)所有波形高斯分量的疊加，進而獲得波形參數(shù)，結果如圖7(b)所示。

圖7 偽大光斑波形擬合Fig.7 Pseudo large footprint waveform fitting

激光雷達系統(tǒng)記錄波形是按時間先后順序依次對光斑內(nèi)各次反射信號數(shù)字化采樣而成，激光雷達穿透首次回波對應的冠層到達末次回波對應的地面，包含森林垂直結構信息。由此結合圖7(b)可以獲得偽大光斑波形長度（Extent）。

3.3 模型構建及精度評價

基于偽大光斑波形長度（Extent），建立平均樹高回歸模型（表3），并結合式（6）對模型進行精度評價，結果如圖8所示。

表3 模型及精度評價Table3 Model and accuracy evaluation

圖8 驗證樣方平均樹高預測值與實測值Fig.8 Verify sample average tree height predictions and measured values

從模型擬合結果來看，基于偽大光斑所建模型反演精度優(yōu)于小光斑波形反演精度，本研究的研究區(qū)地形相對平坦，研究對象為天然次生白樺林，樹干修直，枝葉堆疊較少，通透性好，使得末次回波脈沖能夠順利到達地面而較少地被阻擋，故通過波形長度能夠很好地反演林分平均高,且具有較高的可靠性。

分析發(fā)現(xiàn)，實測平均樹高值與預測值很接近，但也偶有偏差較大的情況，猜測這種現(xiàn)象是因為本研究是對小光斑激光雷達波形數(shù)據(jù)先波形分解后再形成的偽大光斑波形數(shù)據(jù)，此過程會有信息的缺失現(xiàn)象，主要是由原始小光斑波形數(shù)據(jù)預處理（如去噪、平滑）造成且無法避免。具體表現(xiàn)為去噪時噪聲閾值過大則會導致有效信息缺失，閾值較小則會使有效波形數(shù)據(jù)包含過多噪聲。平滑時，高斯濾波器選用高斯函數(shù)較脈沖寬度小，使估測平均樹高略低于實測平均樹高。但總體而言，基于偽大光斑波形能夠準確反演林分平均高，且精度較波形反演模型有所提高，本研究提出的將樣地內(nèi)所有高斯分量特征參數(shù)按脈沖能量比加權平均得到偽大光斑波形數(shù)據(jù)的疊加方法切實可行。

4 結論與討論

為了探索小光斑LiDAR 波形數(shù)據(jù)在估測森林冠層結構參數(shù)中的應用潛能，本研究提出形成偽大光斑波形數(shù)據(jù)的方法，基于小光斑波形數(shù)據(jù)和偽大光斑波形數(shù)據(jù)分別建立林分平均高反演模型，并比較分析。結果表明，偽大光斑模型反演精度高于波形反演模型（R2=0.61，P=90.65%）。

研究林分平均高提取方式發(fā)現(xiàn)，已有很多運用成熟的遙感技術手段。如：Kodani 等[13]使用LiDAR 在常綠闊葉林林分估計了平均株高，決定系數(shù)為0.79；Lin 等[14]建立了平均冠層高度的隨機森林回歸估計模型，總體平均精度93.17%；胡凱龍等[15]基于GLAS 大光斑波形數(shù)據(jù)實現(xiàn)林分平均高的提取，相關系數(shù)0.66；Fang 等[16]基于GLAS 大光斑波形數(shù)據(jù)估測山地森林冠層高度，決定系數(shù)為82.8%。分析可知，不同遙感技術手段估測林分平均高度的精度較本研究有高有低，研究方法也各有優(yōu)劣，在應用時需根據(jù)實際數(shù)據(jù)情況選擇合適的方法。本研究實現(xiàn)了小光斑機載LiDAR 波形數(shù)據(jù)估測林分平均樹高，將小光斑波形數(shù)據(jù)擬合成偽大光斑，既彌補了小光斑因光斑過小而易錯失樹頂信息缺陷，又發(fā)揮了小光斑波形數(shù)據(jù)精確刻畫林木特征的優(yōu)勢，在保證精度的同時為定量估測森林冠層結構參數(shù)提供了新的方法和思路。

本研究研究對象為天然次生白樺林，樣地樹種單一，研究區(qū)域地形結構相對平坦，對于其他地形區(qū)域樹種的適用性有待驗證。同時基于小光斑機載LiDAR 全波形數(shù)據(jù)估測林分平均樹高時，對波形數(shù)據(jù)處理效果會在一定程度上影響波形參數(shù)的準確提取，而建立林分平均樹高估測模型所選取的模型參數(shù)波峰長度（PeakLeg）、波形前緣長度（LeaEdg）和波形后緣長度（TraEdg）可直接影響估算結果的精度。脈沖回波信號由于受到大氣、航高航速及系統(tǒng)等因素的影響，不可避免地產(chǎn)生大量的噪聲，且對于回波信號較弱的波形，有效信號受到噪聲的影響往往十分嚴重，不能有效分離噪聲將會影響波形分解的精度甚至產(chǎn)生錯誤的分解結果。本研究基于高斯分解的方法完成了全波形數(shù)據(jù)的處理流程，取得了很好的效果，但在采用四分位數(shù)方法去噪的過程中，認為對噪聲的估計偏大，導致有效波形信息缺失，使模型參數(shù)波形前緣長度（LeaEdg）和波形后緣長度（TraEdg）偏小，影響了平均樹高的估測精度。在波形平滑處理中，選擇了比發(fā)射脈沖寬度稍小一點的高斯函數(shù)進行濾波，使模型參數(shù)波形前緣長度（LeaEdg）偏小，低估了林分平均高度。因此，未來應著眼于波形數(shù)據(jù)處理方法的改進研究，并嘗試擴大研究范圍，將該方法應用于復雜地形人工林進行檢驗，以期在挖掘小光斑激光雷達應用于森林結構參數(shù)估測時能夠進一步提高估算精度。