本發(fā)明屬于機載激光雷達森林遙感應用領域，具體涉及一種基于機載激光點云聚集關系的單木分割方法。

背景技術：

機載激光雷達(airbornelightdetectionandranging，a-lidar)是通過搭載在航空飛機上的傳感器所發(fā)出的激光來測定傳感器與目標物之間距離的一種有效且主動的遙感技術，不僅可以精確、快速的獲取地球表面幾何信息(如三維地理坐標)、回波地物的輻射信息，而且可以應用于大尺度上，再加上其具有很強的穿透森林植被的能力,為森林結構參數(shù)的提取以及森林生態(tài)系統(tǒng)功能參數(shù)的反演提供了全新的技術手段。

在林業(yè)應用方面，利用機載激光雷達數(shù)據(jù)進行研究大致包含森林參數(shù)反演(樹高、冠幅、葉面積指數(shù)等)、單木分割、林木判別以及森林3d模型重建。林木判別是從原始激光點云中識別出植被回波點的方法；而單木分割算法主要用來解決的問題是在林木判別后得到的植被回波點中分割出屬于特定單棵樹木的點云。在實際野外測量中，比如lai值、樹高、冠幅等實測參數(shù)往往都是針對單棵樹木進行人工測量的，而在利用a-lidar或者其他遙感數(shù)據(jù)反演森林參數(shù)過程中，基本上都是以這些單木實測數(shù)據(jù)源作為驗證數(shù)據(jù)，而且單木分割的位置精度、點云數(shù)據(jù)的密度等也影響著森林參數(shù)反演的精度,因而單木分割這個步驟往往是不能忽視的。

機載激光雷達數(shù)據(jù)不僅提供了高精度的地理位置數(shù)據(jù)，并且采用“自上而下”的方式掃描植被冠層，得到的回波數(shù)據(jù)中包含描述森林冠層的垂直分布特性的信息，使得在森林樣方下的單木分割成為可能。

目前國際上存在的單木分割方法中，最常用的方法是根據(jù)植被冠層的冠層高度模型(canopyheightmodel，chm)進行單木分割，大致過程是它首先通過建立的數(shù)字表面模型(digitalsurfacemodel，dsm)與數(shù)字高程模型(digitalelevationmodel，dem)相減得到chm，然后從chm中定位出局部最大值，并自定義一個范圍來擴展至被識別單棵樹的冠幅，以此冠幅為基礎提取出該棵樹。使用chm方法進行單木分割時，由于建立的chm常常會使用插值點代替真實的三維空間點，會造成不可逆的信息丟失，再加上自定義的窗口大小對chm進行處理，窗口大小與樹高和冠幅大小相聯(lián)系，因而在利用chm進行分割時往往造成錯誤分割以及不適當?shù)男螤罘指睢?/p>

而在三維空間上直接對原始點云進行分割發(fā)展的越來越好，比如h.lee等人(2010)，根據(jù)局部最大值的點與點之間的距離提出可適應性聚類方法直接對植被點云數(shù)據(jù)進行聚類，使其單木分割的總體精度達到了95.1％。lu等人(2014)基于強度和激光點云數(shù)據(jù)的在空間分布的3d結構，提出了先識別樹干的自低向上的方式分割闊葉林，結果表明該方法能從原始點云中檢測出84％的樹木來，而在檢測出的樹木中進行分割的正確率達到了97％。li等人(2012)使用小光斑離散回波激光雷達點云數(shù)據(jù)，通過考慮兩棵樹之間的相對間隔，以及形狀指數(shù)來判斷激光點，結果表明該方法能從原始點云中檢測出86％的樹木來，正確分割率達到了94％。

技術實現(xiàn)要素：

針對上述存在問題或不足，本發(fā)明提供了一種基于機載激光點云聚集關系的單木分割方法，基于原始點云數(shù)據(jù)，根據(jù)單棵樹內(nèi)的機載激光點云在三維空間上的聚集關系以及鄰近的樹與樹之間的點云在三維空間上的聚集關系進行單木分割，提供了單木分割的新路徑。

本發(fā)明具體技術方案如下：

步驟1、林木判別：

由于單木分割算法只是針對植被點云數(shù)據(jù)的，把回波點分類為植被回波點、地面回波點和其余回波點。

步驟2、局部最大值法初始分割：

在未編碼(默認編碼為0)的原始點云中尋找回波點類別屬于植被的且坐標z值最大的點，并記錄該點的坐標序號，同時把該點編碼設為k(k>＝1)，用半徑為r(實測目標樹木數(shù)據(jù)優(yōu)勢植被的最小冠幅)的圓標識該區(qū)域，并且把該區(qū)域內(nèi)的所有植被回波點都設為編碼k；然后在原始點云中去除已經(jīng)編碼了的點，樹編碼號變?yōu)閗+1，以此循環(huán)下去，直到原始點云中所有植被點的編碼都不為0。

步驟3、刪除錯誤編碼樹：

當同一編碼樹的點云數(shù)量小于激光點云密度或兩顆樹之間的間隔小于3r，修改該樹的編碼為0。

步驟4、聚類：

a、使用步驟3得到編碼數(shù)據(jù)，依據(jù)樹編碼號選取任意兩顆樹并確定兩樹之間的間隔i。

b、if判斷，如果i小于maxr，maxr為實測目標樹木數(shù)據(jù)優(yōu)勢植被的最大冠幅半徑,則依次執(zhí)行①、②、③、④。

①、由于在樹冠邊緣的不規(guī)則，高樹相對矮樹樹高的邊緣植被點與其頂點形成的垂直夾角θ是不同的；首先比較兩樹的高矮，以低樹的中心點為原點，不考慮高程坐標z，建立笛卡爾直角坐標系，在高樹頂點與低樹頂點連線上的點云作為高樹的邊緣植被點(如圖4所示)。②、求最小高樹邊緣植被激光點與高樹頂點形成的閾值正切值mintanθ(原理如圖3),即以mintanθ為閾值tanθ。

③、記錄與mintanθ比較的激光點：

為了減少聚類計算時間以及矮樹邊緣異值點對聚類結果的影響：記錄兩樹頂點連線為長邊長度及方向，矮樹樹冠幅為寬邊的平行四邊形所屬激光點point2peek，且該點的編碼不為高樹編碼的點(如圖4所示)。

④、根據(jù)多木判別原理進行聚類：首先判斷正切值，如果步驟③中記錄的激光點對應的正切值大于mintanθ，判斷該點編碼對應的矮樹為高樹的一部分，并把該矮樹對應的點編碼改為高樹的編碼，將平行四邊形內(nèi)所屬激光點遍歷(原理如圖3)。

⑤、根據(jù)單木歸屬原理進行聚類：如果point2peek對應的正切值小于mintanθ，那么判斷該點為高樹的冠幅點，修改該點的編碼為高樹的編碼，遍歷直到修改完該棵樹的所有回波點(原理如圖2)。

d、重復執(zhí)行a、b步驟，直到遍歷完所有組樹。

步驟5：最近鄰點法處理，把未被編碼的激光點分到距離最近的已編碼聚類的樹中。

步驟④和步驟⑤的聚類原理

單木歸屬原理：在同一棵樹中，假設在高度為h時,邊緣植被點與頂點形成的垂直夾角為θ，邊緣植被點以內(nèi)的點與頂點形成的夾角為θ`；在高度低于h時，冠幅內(nèi)的點與頂點形成的夾角為a、b(如圖2所示)。在同一高度h，θ`大于θ角，在高度小于h時，a、b小于θ角。若已知高度大于h的點屬于該棵樹的冠層點，那么可以求出在高度為h時邊緣點與頂點形成的最小角度θ；根據(jù)高度小于h的其余點和頂點形成的夾角a、b與θ的關系，可以判斷其余點與該樹的所屬關系。

多木判別原理：在一組樹中，與矮樹頂點同一高度上，高樹的邊緣植被點與其頂點形成的垂直夾角為θ；高樹頂點與矮樹頂點形成的夾角為a；矮樹其余點與高樹頂點形成的夾角為b(如圖3所示)。如果在矮樹中的角a小于θ，且存在部分角b小于θ，那么矮樹不劃歸為高樹冠層部分；否則把矮樹劃歸為高樹。

本發(fā)明首先使用terrascan軟件進行林木判別，把原始機載激光雷達數(shù)據(jù)分類為地面回波點、植被回波點和其余回波點；然后只針對植被回波點利用局部最值法進行初始的單木判別，同時在此過程中刪除錯誤判別的樹；再利用單棵樹內(nèi)的點云以及鄰近的樹與樹之間點云在三維空間上的聚集關系進行聚類處理，從而得到較為精確的分割結果；最后跟據(jù)最近鄰點法把未分類的回波點劃分到離其最近的單木中，從而提取出整片區(qū)域的樹木來，其流程如圖1所示。

本發(fā)明是通過利用機載激光雷達獲取森林植被冠層三維激光點云數(shù)據(jù)，并根據(jù)這些三維點云數(shù)據(jù)在單棵樹內(nèi)的聚集關系以及鄰近的樹與樹之間的聚集關系，經(jīng)歷林木判別、局部最大值法初始分割、刪除錯誤的分割單樹、循環(huán)聚類和最近鄰點法處理，建立基于原始三維激光點云數(shù)據(jù)的從森林中分割出單棵樹的單木分割方法。

綜上所述，本發(fā)明相比chm單木分割方法，能避免信息丟失、減少錯誤分割和不適當?shù)男螤罘指睿指钚Ч麅?yōu)異，并且計算時間也得到了大幅度的縮減，能快速、準確的對森林激光點云數(shù)據(jù)分割出屬于特定單棵樹的回波激光點來，有助于單樹的參數(shù)反演(如葉面積指數(shù)、樹高、生物量等)、單樹的三維結構特征重建以及對森林植被參數(shù)反演的驗證。

附圖說明

圖1為本發(fā)明的流程示意圖；

圖2為單木歸屬原理圖；

圖3為多木判別原理圖；

圖4為聚類過程中的確定兩棵樹的相對方向以及取值范圍；

圖5為森林樣方林木判別結果圖；

圖6為本發(fā)明的操作界面圖；

圖7為初始識別和刪除錯誤編碼樹后的結果圖；

圖8為單木分割最終結果圖；

圖9為使用實測樣方內(nèi)的林木棵數(shù)來驗證分割精度的結果圖；

圖10為使用實測樹高來驗證分割精度的結果圖。

具體實施方式

以下通過一個樣方的實施并結合附圖對本發(fā)明作進一步詳細說明：

開發(fā)環(huán)境microsoftvisualstudio2010、開發(fā)平臺框架.netframework3.5、三維可視化引擎directx11、編程語言c#。

步驟1，以內(nèi)蒙古呼倫貝爾市西北部，額爾古納市東南部的上庫力農(nóng)場為研究區(qū)域，研究區(qū)域地理位置坐標為東經(jīng)120°36′50.48″～120°52′56.53″，北緯50°21′11.08″～50°24′32.00″，并選取實地測量的30個樣方作為單木分割驗證區(qū)，其優(yōu)勢植被為白樺樹。利用leicalas60掃描系統(tǒng)獲取該區(qū)域的機載激光雷達數(shù)據(jù)，航飛時間在2012年8月16日至2012年9月25日，系統(tǒng)具體參數(shù)參見表1。并基于terrascan軟件把原始激光回波點分類為植被、地面、其余回波點，分類結果如圖5。以下步驟為對樣方編號為yf_1，中心坐標為(120.770998°，50.357450°)的單木分割詳細步驟。

表1機載激光雷達系統(tǒng)參數(shù)

步驟2，在原始激光點云中尋找局部最大高程值點，對該點進行編碼并以優(yōu)勢植被最小冠幅的近似圓標識該區(qū)域(樣方實測最小冠幅為1m)，以此循環(huán)下去，直到所有的原始點云都被編碼，最后得到初始識別的樹木棵樹為197(圖6所示)。采用上述局部最值法進行單木分割，由于只是使用優(yōu)勢植被的最小冠幅r來標識各個區(qū)域，而r并不能代表所有樹的冠幅，那么得到的分割結果將會存在多分現(xiàn)象、錯誤識別現(xiàn)象以及形狀不規(guī)范性現(xiàn)象等。為了提高分割精度，需要其他的方法對初始分割結果進行后處理。

步驟3，在初始識別得到的結果中，如果所識別的樹包含的激光回波點數(shù)小于5個點(這個點數(shù)閾值根據(jù)點云密度和最小冠幅面積而設定)，就把屬于該棵樹的所有編碼設為未編碼狀態(tài)，編碼值設為0。當兩棵樹之間的距離小于1.5米時，取消矮樹的所有編碼，最終得到修改后的樹木棵樹為45棵，如圖6、圖7所示。

步驟4，在實測中該樣方單棵樹的最大冠幅為4.2m，這里把最大半徑閾值設為2.1m，根據(jù)激光點云在單樹或者相鄰樹的聚集情況進行聚類處理，最后得到的分割棵樹為38棵，而實測棵樹為33棵。聚類結果如圖7所示。

步驟5，根據(jù)距離把未被編碼的回波點分到最近的已編碼已聚類的樹中，如圖7所示。

步驟6，精度檢驗。分別從單木分割后得到的樣方內(nèi)樹木棵樹和分割后得到的點云數(shù)據(jù)進行單木樹高的反演兩個方面對分割結果的精度進行了檢驗。如圖8、圖9所示。在圖8中，對實地測量的30個樣方進行分割，分割后得到的樹木棵樹的精度分別是，絕對系數(shù)r²為0.7398，均方根誤差rmse為6.845棵。在圖9中，對分割后得到數(shù)據(jù)進行了單木提取以及單木樹高反演，樹高反演精度分別是，絕對系數(shù)0.7912，均方根誤差1.894米。

該實例根據(jù)本發(fā)明提出的方法，對白樺樹樣方的機載激光雷達點云數(shù)據(jù)進行分析，依照技術方案所述，實現(xiàn)了對30個樣方的單木分割，并對分割后的結果進行了兩種方法的精度檢驗，精度檢驗結果分別為(r²＝0.7398，rmse＝6.845棵)、(r²＝0.7912，rmse＝1.894米)，都取得了較好的精度，充分說明該方法的可適用性。由檢驗結果可知，在一般情況下，分割后的單木棵樹比實測的稍多，是因為在分割時，如果單棵樹存在比較大的分支，那么單棵樹內(nèi)的回波點云聚集情況將會改變?yōu)橄噜彉涞狞c云聚集情況，造成多分，尤其是多枝丫的樹出現(xiàn)多分的頻率會更高。樹高反演的均方根誤差稍微顯大，是因為樹高反演的數(shù)據(jù)來源是經(jīng)過單木分割后得到的點云數(shù)據(jù)，單木分割中本身就存在一定的誤差(如錯誤分割、少分、多分等)，這樣就導致單木分割中的誤差傳遞到單木樹高反演中，進一步帶來偏差。該分割算法更適用于不存在大的分支的樹種(針葉林)或者樹木冠層形狀比較規(guī)范的闊葉林等。