瞿 帥，張曉麗，朱程浩，霍朗寧，劉會玲

(北京林業(yè)大學 精準林業(yè)北京市重點實驗室，北京 100083)

森林是陸地生態(tài)系統(tǒng)的主體，森林及其變化對陸地生物圈及其他地表過程有著重要影響[1]。森林資源的可持續(xù)開發(fā)利用依賴于有效的森林資源調查和經營，而傳統(tǒng)的森林資源調查方式需要投入大量的人力和物力，調查周期長，數據更新速度遠遠不能滿足當今信息化時代對森林資源調查數據的實時更新和獲取需求[2-4]。

隨著遙感技術的發(fā)展，森林資源調查擁有了強大的技術支撐[5]。近年來，光學遙感技術已廣泛應用于森林參數(如森林蓄積量、葉面積指數、生物量等)的提取及反演[6-11]，但是光學遙感受多種因素的影響，如云、混合像元等[12],而且由于光學遙感不具備三維空間探測的能力，從而限制了光學遙感在森林結構參數提取及反演方面的應用。

激光雷達是一種主動方式的遙感探測技術，具有較強的穿透能力，能夠用來探測森林空間結構及林下地面信息，具有光學遙感無可比擬的優(yōu)勢[13],在森林資源調查及經營管理方面具有重要的應用價值[14]。

現(xiàn)階段利用機載激光雷達進行森林資源調查的工作已經展開，但對于獲取的大數據量的激光雷達點云數據面向林業(yè)需求的處理工作往往比較滯后，如何有效利用機載激光雷達數據及時獲取成果以保障森林資源調查工作的如期完成是一個亟待解決的問題。因此，本研究將結合實際生產需要，利用機載激光雷達數據結合地理信息系統(tǒng)技術進行森林資源調查系統(tǒng)的設計與試驗，以期實現(xiàn)自動化、高精度的森林資源調查，滿足森林資源數據庫及時更新的需求。

1 系統(tǒng)分析與設計

1.1 功能需求分析

林業(yè)資源調查獲取的機載激光雷達點云數據一般數據量較大，系統(tǒng)需要在足夠大的內存下，支持千萬級別數量的點云數據載入，并進行流暢的數據渲染顯示；能夠加載影像；支持不同格式點云數據的轉換；可以實現(xiàn)點云去噪并進行高精度的點云濾波處理，生成數字高程模型；從單木尺度實現(xiàn)單木定位、單木樹高、冠幅的提取，從林分尺度實現(xiàn)株樹密度、林分平均樹高的提取。

1.2 系統(tǒng)結構設計

本系統(tǒng)以QT為圖形界面框架，結合VTK三維圖形引擎庫與GDAL影像處理庫，在Visual Studio 2013集成開發(fā)環(huán)境下采用C++語言開發(fā)。系統(tǒng)結構設計使用3層結構，分別是應用層、邏輯層、數據層。

應用層也可以稱為表現(xiàn)層，主要表現(xiàn)為人機交互，是對用戶公開的部分；邏輯層主要實現(xiàn)數據組織管理、各種數據處理算法；數據層主要實現(xiàn)不同數據格式的數據與本系統(tǒng)進行數據交換的接口。本系統(tǒng)所使用的數據以文件形式提供，系統(tǒng)為不同類型及格式的數據設計相應的數據解析接口，并封裝成數據操作類庫，方便本系統(tǒng)對外部數據進行操作管理及分析，并具備較高的可擴展性。3層結構間的關系具體表現(xiàn)為：應用層調用邏輯層實現(xiàn)的相關數據處理算法完成應用層的功能實現(xiàn)，邏輯層通過數據操作類庫完成對數據層相關數據的組織及管理。系統(tǒng)結構示意圖如圖1所示。

圖1 系統(tǒng)結構

機載激光雷達森林資源調查系統(tǒng)主要分為數據輸入、地形信息提取、林木參數提取、成果輸出4大功能模塊。系統(tǒng)主界面如圖2所示。

1.3 核心模塊實現(xiàn)

1.3.1 地形信息提取模塊 機載LiDAR能夠穿透植被遮擋，直接獲取地面點信息，因而成為數字高程模型(DEM)獲取的首選方式，而原始點云數據中還包括建筑、植被等非地面信息，要獲取地面點，就必須采取一定的算法將非地面點從原始點云中剔除，這個過程稱為點云濾波[15]。常規(guī)點云濾波算法主要包括數學形態(tài)學濾波[16-17]、逐步加密濾波[18-20]等算法。針對林區(qū)地形特征較為復雜的特點，本系統(tǒng)中實現(xiàn)了一種漸近不規(guī)則三角網加密濾波算法，算法流程如圖3所示。

圖2 系統(tǒng)主界面

圖3 算法流程

算法基本步驟：

1)對原始點云數據進行中值濾波去噪，去除由于系統(tǒng)誤差或者飛鳥激光反射點等噪聲點云，避免對后續(xù)步驟算法的影響。

2)點云格網化，設定一格網尺寸，根據點云中空間點的X、Y坐標將去噪后點云劃分入規(guī)則格網內，并將每一規(guī)則格網內高程最低點作為初始地面種子點。

3)將初始地面種子點作為構建初始不規(guī)則三角網的種子點，在構建初始不規(guī)則三角網過程中，引入坡度閾值判斷，即在構建三角網時先進行相鄰構網節(jié)點的坡度計算，如果小于坡度閾值(此處坡度閾值應設定較大值)則該點加入構網，否則標記為非地面點，以后步驟中不再參與運算，最終獲得初始不規(guī)則三角網。

4)在初始不規(guī)則三角網基礎上，對未劃入地面點云的點進行坡度及高差閾值的判斷，如果滿足條件則劃入地面點，對地面不規(guī)則三角網進行迭代加密。

5)重復步驟4，當不再有新的地面點加入三角網時，停止迭代，得到滿足條件的三角網點云，即為最終地面點云。

在軟件實現(xiàn)方面，本系統(tǒng)中將此算法封裝成一個算法類，格網尺寸、坡度閾值、高差閾值等參數可在軟件濾波算法界面上進行設定，作為算法傳入參數，執(zhí)行濾波最終得到地面點云。為了根據獲得的地面離散點云得到數字高程模型，系統(tǒng)采用GDAL影像處理庫中離散三維點反距離加權插值算法將地面點云插值形成DEM影像，以標準GeoTIFF格式輸出。插值函數為GDALGridCreate，插值算法設置為GGA_InverseDistanceToAPower。

1.3.2 林木參數提取模塊 此模塊主要包含單木參數提取及林分參數提取功能，單木參數主要包含單木樹高、冠幅以及單木位置，林分參數主要包括株樹密度及林分平均高。

1.3.2.1 單木參數提取 機載LiDAR點云數據已成功用于森林垂直結構參數及水平結構參數的估測[21]，可以直接提取單木樹高及冠幅。在點云濾波處理的基礎上可以得到林區(qū)點云數據包含的地面點云和非地面點云，林區(qū)非地面點云基本相當于植被點云，對原始植被點云進行高程歸一化處理，則輸出點云中每個空間點的高程即為該點相對于地面的絕對高度，然后利用歸一化點云進行單木點云分割(本系統(tǒng)中實現(xiàn)了一種基于歸一化植被點云高度分層的K-Means聚類點云分割算法)，最終依據單木點云分割結果進行樹高和冠幅的提取。算法流程如圖4所示。

圖4 算法流程

算法基本步驟：

1)點云歸一化：將點云濾波得到的植被點云數據以水平X、Y坐標為依據，根據生成的DEM影像分辨率尺寸劃分到不同格網，每個落在格網內的植被點高程減去對應DEM格網的值，即得到歸一化后的植被點云，此時植被點云中點的高度就相當于植被高度。

2)單木點云分割：按照歸一化后植被點云高度分布情況，將植被點云分成若干層，在每一層設定一鄰域檢測窗口尺寸，進行高度局部最大值檢測，對每一層的植被點云以局部最大值點為初始聚類中心進行三維空間點K-Means聚類，在每層點云聚類完成之后，從最上層開始，判斷相鄰上下2層各聚類點云中各聚類中心間水平距離是否小于設定的聚類中心距離閾值，如果小于閾值則合并上下2層中對應的聚類點云，直到所有層比較合并完畢，最終得到單木分割點云。

3)參數提取：對得到的單木分割點云，統(tǒng)計每個單木點云的最高點高度，將其作為提取的樹高，并將該點的水平坐標作為單木定位坐標，計算每個單木分割點云投影到水平面上的凸包面積，以該面積為樹冠近圓形投影面積，計算圓直徑作為單木平均冠幅。

在本系統(tǒng)中將點云歸一化、單木點云分割以及樹高、冠幅提取分步進行，針對每個子功能封裝為一個算法類，并且設置了不同的軟件界面窗口，用于設定各子功能參數，最終實現(xiàn)了單木樹高及冠幅參數的提取，提取結果以表格形式輸出，記錄了每棵樹的位置信息、樹高、平均冠幅。

1.3.2.2 林分參數提取

1)株樹密度：在實際林業(yè)調查中常使用圓形樣地法[22-23]來確定局部區(qū)域的林木株樹密度，該方法可以提高調查效率。本研究中以方形樣地為調查區(qū)域，根據樣地單木分割結果，可以統(tǒng)計出樣地的林木株樹，通過式(1)計算出樣地的株樹密度值。

(1)

式中,N為株樹密度(株/hm2)；n是樣地內樹木株數；S是樣地面積(hm2)。

2)林分平均高：本系統(tǒng)中實現(xiàn)的林分平均高提取方法，以龐勇[24]等的研究為依據，樣地歸一化植被點云數據上四分位數處的高度與實測樹高相關性高。

先選定一定數量的樣地調查數據作為訓練樣本，計算樣地歸一化植被點云上四分位高度，具體計算為：對樣地范圍內的歸一化植被點云根據高度進行排序，然后計算總高度的上四分位處的高度[20]，即為樣地歸一化植被點云上四分位高度，再將其與實測林分平均高建立線性回歸方程，然后可以用預留檢驗樣本進行精度驗證。最終根據樣地點云及回歸方程實現(xiàn)待測樣地的林分平均高提取。

其中，實測林分平均樹高計算采用斷面積加權計算方法，計算公式如下：

(2)

式中，H為林分平均高，hi為第i棵樹的樹高，gi為第i棵樹的胸高斷面積，k為林分株數。

2 試驗區(qū)案例研究

2.1 試驗區(qū)概況

試驗區(qū)位于甘肅省張掖市大野口森林區(qū)域(97°20′-103°50′E，36°40′-39°50′N)，位于祁連山自然保護區(qū)，海拔2 700～3 200 m。樣地內樹種為天然青海云杉純林，分布茂密，林齡組成結構主要為成熟林。林下植被種類單一，地表覆蓋物主要為苔蘚。

2.2 數據采集及預處理

2.2.1 LiDAR數據獲取 本次試驗區(qū)的機載激光雷達數據由RIEGL公司研制的LiteMapper 5600系統(tǒng)進行測量，配置的激光掃描儀為Riegl LMS-Q560，工作波長1 550 nm，光斑平均直徑約為0.38 m。飛機的相對航高約700 m，整個測區(qū)點云的平均密度為3.43個/m2，地表點的定位精度為水平方向0.1 m、垂直方向0.03 m。LiDAR點云數據采用WGS84坐標系和UTM投影(北半球6度第48帶)。

2.2.2 樣地調查數據 設立一邊長100 m正方形范圍的大尺寸樣地，為方便測樹調查，將其劃分成16個子樣地，在試驗區(qū)已知的大地控制點上建立GPS基準站，使用靜態(tài)差分GPS對大樣地四角點及各子樣地的中心坐標和四角點位置坐標進行精確定位，以子樣地為單位進行每木測樹調查，數據主要包括2個部分：一是林木實測數據，包含實測樣地內單木胸徑、樹高以及單木冠幅(南北、東西)等。二是地面定位點數據，采用全站儀對樣地內所有掛牌林木的地面位置進行精確測定。利用大樣地及子樣地角點定位坐標，使用系統(tǒng)點云裁剪功能獲取大樣地及每個子樣地的原始點云數據。大樣地點云數據如圖5所示，以高程渲染顯示。

圖5 點云數據

2.3 地形信息精度驗證

本次研究對系統(tǒng)的數字高程模型提取值進行了精度驗證，以16塊子樣地中心位置實測高程為真值，將該系統(tǒng)生成的數字高程模型對應位置的高程值與實測數據作了對比(表1)。對比數據分析可得，系統(tǒng)高程提取值的最大偏差為1.102 m，平均偏差為0.737 m。

2.4 林木參數精度驗證

為驗證軟件系統(tǒng)提取的單木樹高及冠幅，將該試驗區(qū)實測數據與軟件提取結果進行了對比。根據統(tǒng)計，大樣地共調查1 465株樹，去除枯立木、死木和缺少單木位置的單木，同時對于融合的樹進行拆分，共獲得1 445株樹，大樣地單木分割結果共提取1 223株樹(圖6)。根據統(tǒng)計輸出單木提取位置信息依據最近鄰匹配原則，在提取單木中選取距離實測單木位置最近單木作為其匹配單木，將匹配單木實測樹高與提取值進行對比，可得單木樹高平均相對誤差為13.66%；將實測單木的東西、南北向冠幅平均值與提取單木樹冠進行對比，可得單木冠幅的平均相對誤差為14.18%(表2)。

表1 軟件提取高程與實地測量值比較

圖6 單木分割結果

對于各子樣地中株樹密度的精度驗證，根據大樣地單木分割結果以及16個子樣地中各匹配單木，將其與16個子樣地實地調查結果進行對比，可得軟件提取的株樹密度與實地調查株樹密度值平均相對誤差為5.83%；而在軟件提取林分平均高參數時，選取1～10號共10塊子樣地作為訓練樣本，構建樣地歸一化植被點云上四分位高度與實測林分平均高的線性回歸方程為Y=0.727 7X+5.659 9，R2為0.752 9，其中Y為實測林分平均高，X為樣地歸一化植被點云上四分位高度，利用剩余6個子樣地點云結合回歸方程計算林分平均高并與實測數據進行對比，可得軟件提取的林分平均高與實測林分平均高的相對平均誤差為9.86%(表3)。

表2 單木參數提取與實地測量值比較

表3 林分參數提取與實地測量值比較

3 結論與討論

基于機載激光雷達的森林資源調查系統(tǒng)充分利用激光雷達可以探測森林垂直空間結構的能力，體現(xiàn)出了其在林業(yè)調查領域的優(yōu)勢。與傳統(tǒng)調查方式相比，該系統(tǒng)可以快速、高效的獲取調查區(qū)域的地形信息及相關林木參數信息。以甘肅省張掖市大野口關灘森林站超級樣地試驗區(qū)為例，進行機載激光雷達森林資源調查系統(tǒng)的驗證研究，該系統(tǒng)可以有效獲取試驗區(qū)域地形信息，提取單木樹高、冠幅的平均相對誤差分別為13.66%、14.18%，提取株樹密度、林分平均高的平均相對誤差分別為5.83%、9.86%。

針對本系統(tǒng)進行林木參數提取方面，由于樣區(qū)內點云密度不高，去除樣區(qū)地面點之后，實際林木點云更加稀疏，而且激光點不能夠保證掃描到樹頂，因而會造成利用單木分割算法提取計算的單木樹高普遍偏低，同樣也會造成單木樹冠提取結果偏低，而當有高大樹木遮擋鄰近較矮樹木時，單木分割時，可能會將其劃分為一棵樹木，有可能會造成個別提取的單木冠幅大于實測值。對于林分參數，株樹密度的精度與單木分割的結果有直接的關系，而以樣地歸一化植被點云上四分位高度結合回歸方程計算林分平均高的精度對比單木樹高提取精度有所提高，與激光點云不一定能夠掃描到樹頂從而影響單木樹高提取精度的分析相印證。

現(xiàn)階段，本系統(tǒng)尚有諸多不足之處，僅僅針對林分結構單一的針葉純林進行了林木參數提取及驗證，對于不同植被類型以及林分結構較為復雜的情況，系統(tǒng)中設計實現(xiàn)的單木分割算法具有一定的局限性，林木參數提取的精度會有較大程度的降低。因此后期考慮針對不同植被類型，將機載激光雷達系統(tǒng)搭載的光學相機獲取的高分辨率光學影像與激光雷達數據有機結合，特別是針對闊葉林，單純利用激光點云數據難以進行較高精度的林木參數提取，考慮以光學影像進行單木樹冠分割，然后以單木樹冠分割結果限定單木激光點云的劃分，再進一步利用限定后的點云進行林木參數的提取，從而提高林木參數的提取精度，這將是下一步研究的重點方向。