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

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

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

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

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

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

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

1.1 功能需求分析

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

1.2 系統(tǒng)結(jié)構(gòu)設(shè)計

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

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

圖1 系統(tǒng)結(jié)構(gòu)

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

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

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

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

圖3 算法流程

算法基本步驟：

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

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

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

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

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

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

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

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

圖4 算法流程

算法基本步驟：

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

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

3)參數(shù)提?。簩Φ玫降膯文痉指铧c云，統(tǒng)計每個單木點云的最高點高度，將其作為提取的樹高，并將該點的水平坐標(biāo)作為單木定位坐標(biāo)，計算每個單木分割點云投影到水平面上的凸包面積，以該面積為樹冠近圓形投影面積，計算圓直徑作為單木平均冠幅。

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

1.3.2.2 林分參數(shù)提取

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

(1)

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

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

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

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

(2)

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

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

2.1 試驗區(qū)概況

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

2.2 數(shù)據(jù)采集及預(yù)處理

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

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

圖5 點云數(shù)據(jù)

2.3 地形信息精度驗證

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

2.4 林木參數(shù)精度驗證

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

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

圖6 單木分割結(jié)果

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

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

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

3 結(jié)論與討論

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

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

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