孫拱 辛穎 閆燦

摘?要：為快速提取森林單木樹高，提高單木樹高的提取精度，本研究提出一種改進的 K-means 聚類算法提取單木樹高。即利用內(nèi)蒙古依根地區(qū)的野外實測數(shù)據(jù)和機載 LiDAR 點云數(shù)據(jù)，以各單木實測坐標為該算法的聚類中心，對機載 LiDAR 點云數(shù)據(jù)進行圓柱擬合聚類，將圓柱體內(nèi)垂直方向最高的點視為該單木樹冠最高點，利用提取的樹高與野外實測樹高進行相關(guān)性分析。分析發(fā)現(xiàn)，基于改進的 K-means 算法提取的樹高與實測樹高相關(guān)系數(shù) R2為 0.852。結(jié)果說明改進的 K-means 算法能夠很好的用于機載 LiDAR 點云數(shù)據(jù)單木樹高的提取，降低了傳統(tǒng)測量方法的局限性，挖掘了機載LiDAR點云數(shù)據(jù)的應(yīng)用潛力。

關(guān)鍵詞：機載 LiDAR;改進 K-means 算法;樹高

中圖分類號：S758文獻標識碼：A???文章編號：1006-8023（2019）06-0031-05

Estimating Single Tree Height Based on Improved K-means Algorithm

SUN Gong， XIN Ying*， YAN Can

（College of Engineering and Technology， Northeast Forestry University， Harbin 150040）

Abstract：In order to quickly extract the height of forest single tree and improve the extraction precision of single tree height， this study proposed an improved K-means clustering algorithm to extract single tree height. Using the field observation data and airborne LiDAR point cloud data in the Yigen area of Inner Mongolia， the measured coordinates of each single tree were used as the clustering center of the algorithm， and the cylindrical lie point cloud data was clustered and clustered in the vertical direction. The highest vertical point in the cylinder was regarded as the highest point of the single tree canopy， and the correlation between the extracted tree height and the field measured tree height was analyzed. The analysis found that the correlation coefficient R2 between the tree height and the measured tree based on the improved K-means algorithm was 0.852. The results show that the improved K-means algorithm can be used to extract the single-tree height of airborne LiDAR point cloud data， which reduces the limitations of traditional measurement methods and explores the application potential of airborne LiDAR point cloud data.

Keywords：Airborne LiDAR; improved K-means algorithm; tree height

0?引言

激光雷達（Light Detection and Ranging， LiDAR）是一種主動遙感技術(shù)，目前在森林資源調(diào)查、數(shù)字城市和基礎(chǔ)測繪等方面有廣泛應(yīng)用[1]。通過該技術(shù)能夠直接提取樹高，作為估測其他參數(shù)的基礎(chǔ)數(shù)據(jù)其在測量中存在難度大、耗時長和精度低等問題，故精確獲取單株木的樹高對反演其他森林結(jié)構(gòu)參數(shù)及減少人工工作量具有現(xiàn)實意義[2-6]，而 LiDAR 作為目前實現(xiàn)對地觀測最新技術(shù)之一，對高度的探測獨具優(yōu)勢[7]。

從機載 LiDAR 數(shù)據(jù)中提取森林結(jié)構(gòu)參數(shù)技術(shù)在林業(yè)中應(yīng)用廣泛，而如何從大量點云數(shù)據(jù)中準確提取單木參數(shù)信息更是研究的重點[8]。李增元等[9]提出目前單木分割方法主要分為 3 類：點云歸一化單木法、體元單木法和冠層高度模型（Canopy Height Model， CHM）單木法; Wang 等[10]介紹了一種基于 LiDAR 點云的全自動三維單樹建模的全過程;Brandtberg 等[11]利用機載 LiDAR 點云數(shù)據(jù)提取出 CHM，基于 CHM 柵格圖像采用局部最大值方法確定樹冠頂點;Reitberger 等[12]提出的 Normalized Cut-基于圖論的分割思想;這些研究表明基于機載 LiDAR 數(shù)據(jù)可以實現(xiàn)單木樹高的提取，進而有利于反演其他森林結(jié)構(gòu)參數(shù)并提高森林經(jīng)營的效益。Morsdorf 等[13]以數(shù)字表面模型的局部最大值為聚類中心，通過 K-means 算法聚類實現(xiàn)單木分割，其樹高估測精度 RMS=0.600;趙晨陽等[14]采用 K-means 算法以單木樹冠頂點作為初始聚類中心，經(jīng)過 4 次迭代提取單木樹高，其估測精度為 R2=0.892;基于 K-means 算法的單木樹高估測研究表明，在聚類迭代的過程中機載點云數(shù)據(jù)量大的特點會導致聚類過程不穩(wěn)定且速度較慢，同時由于 K-means 算法較適合對密度分布均勻且大致呈球形的數(shù)據(jù)進行聚類分析，這與機載 LiDAR 數(shù)據(jù)并不相符，故即使按合適比例將點云相對高程值縮放使其分布大體呈球形，仍會影響聚類的精度。而本研究提出改進的 K-means 算法，即以野外實測坐標（x， y）為聚類中心擬合單木圓柱后進行聚類，將聚類結(jié)果限制在圓柱范圍內(nèi)，不需要多次迭代，使聚類速度提高，且由于該聚類過程是投影在二維平面進行，避免了原始 K-means 算法適合于球形空間聚類的缺陷，降低聚類過程中產(chǎn)生的誤差。

本研究以內(nèi)蒙古依根地區(qū)為研究區(qū)，對機載激光雷達數(shù)據(jù)進行去噪、高程歸一化后，利用改進的K-means算法對機載LiDAR點云進行聚類并提取樹高，最后進行相關(guān)性分析。以實現(xiàn)避免樹高測量過程中的資源浪費和提高機載激光雷達獲取單木樹高精度的目的。

1?研究區(qū)概況

本研究區(qū)是位于內(nèi)蒙古境內(nèi)北部、內(nèi)蒙古大草原的上庫力農(nóng)場，經(jīng)度為 120°36′50.48″～120°52′56.53″E，緯度為 50°21′11.08″～50°24′32.00″N，地處內(nèi)蒙古自治區(qū)呼倫貝爾市額爾古納市的上庫力鄉(xiāng)，如圖 1 所示。該農(nóng)場為寒溫帶大陸性季風氣候，海拔為 600～700 m，具有較復(fù)雜的山岳地形地貌特征。山脈丘陵陰坡廣泛分布著以白樺（Betula platyphylla Suk）為主的天然次生林，混生樹種主要包括樟子松（Pinus sylvestris var.mongolica）和落葉松（Larix gmelinii）等，林下灌木層主要由筐柳（Salix linearisipularis）和石棒繡線菊（Spiraea media）等組成。

2?材料與方法

2.1?樣地布設(shè)及數(shù)據(jù)獲取

2.1.1?激光雷達數(shù)據(jù)

本研究所用數(shù)據(jù)是于 2012 年 8 月 26 日進行飛行試驗所獲取的 22 個東西方向條帶數(shù)據(jù)，所用觀測儀器為 Leica RCD105 相機，其鏡頭為 65 mm，飛行搭載所用雷達系統(tǒng)為機載激光掃描系統(tǒng) Leica ALS60，飛行的相對高度為 1 300 m，其波長為 1 064 nm，掃描方式是線性掃描，光斑大小為 0.22 m，點云密度為每平方米 2～4 個脈沖，圖像分辨率為 20 cm，圖像重疊率達 50%以上。研究區(qū)數(shù)據(jù)采集過程如圖 2 所示。

2.1.2?地面數(shù)據(jù)

為驗證飛行數(shù)據(jù)，在 2012 年 8 月，在研究區(qū)中均勻布設(shè)了 33 個半徑為 10 m 的圓形樣方，利用 Trimble GeoXT6000GPS 定位儀對各樣方中心點，樣方各角點及樣方內(nèi)所有單木進行差分定位，并結(jié)合基站數(shù)據(jù)對定位坐標進行差分運算，使差分定位精度在厘米級，滿足激光雷達提取單木位置精度要求，不會出現(xiàn)實測樹木位置偏離的現(xiàn)象。采用 Vertex IV 超聲波樹木測高測距儀測量單株樹高及枝下高，在 30 cm 范圍內(nèi)其高度分辨率達 0.1 m。此次野外調(diào)查共測得單木 1 500 株，其中白樺 1 238 株，落葉松 244 株，樟子松 18 株。表1為地面實測樹高的基本統(tǒng)計量。

2.2?數(shù)據(jù)預(yù)處理

2.2.1?機載 LiDAR 點云去噪

機載激光雷達中將數(shù)據(jù)采集過程中由于飛行物等返射回的極高點信號和低于地面水平的物體（如井壁等）反射回的極低點信號統(tǒng)稱為噪聲點。本研究采用 Isolated points 算法判斷孤立點，即在假定以某點為中心點在其 25 cm 半徑內(nèi)無其他點的前提下可判斷該點為噪聲點，最后對未能剔除的噪聲點采取手工剔除的方法[15]。

2.2.2?機載 LiDAR 點云分類

對去噪后的數(shù)據(jù)進行點云濾波處理，即分離地面點和非地面點。本研究采用不規(guī)則三角網(wǎng)格法（Triangulated Irregular Network，TIN），TIN 是建立數(shù)字高程模型（DEM）和數(shù)字表面模型（DSM）的重要方法之一， TIN 濾波原理是根據(jù)點云間高程差值是否滿足臨界條件和對滿足條件的鄰近點進行分類，即首先選擇一部分最低點作為種子點建立最初的不規(guī)則三角網(wǎng)模型，然后連續(xù)加入滿足臨界條件的點作為種子點，建立新的不規(guī)則三角網(wǎng)模型，進行迭代，直到?jīng)]有新的種子點加入時終止，最終完全分類出地面點和非地面點[16]。

2.2.3?高程歸一化

將分類后 LiDAR 地面點云數(shù)據(jù)按適當比例的單元劃分成柵格圖像，獲取柵格內(nèi)點云數(shù)據(jù)的高程均值，并生成數(shù)字高程模型（digital elevation model， DEM），然后非地面點云數(shù)據(jù)的高程減去 DEM 的高程， 得到植被點云數(shù)據(jù)的高程值，此過程即高程歸一化[17]。

2.3?改進的 K-means 算法

原始 K-means 算法是常用的聚類算法之一，是一種基于距離的迭代式算法，但它會出現(xiàn)聚類結(jié)果陷入局部最優(yōu)解及因采用誤差平方和準則函數(shù)出現(xiàn)的將大聚類簇分割開的情況，且在聚類中心一定的情況下，將該聚類算法用于提取單木樹高，并不能體現(xiàn)出多次迭代的優(yōu)勢。

為了改善上述問題并減少迭代次數(shù)，提高聚類速度同時降低聚類誤差。本研究提出改進的 K-means 算法，利用依根試驗區(qū)野外實測數(shù)據(jù)和激光數(shù)據(jù)進行實驗，算法步驟如下：

（1）確定初始聚類中心。即從 n 個樣本數(shù)據(jù)中選取k個特征數(shù)據(jù)作為中心點，本研究以野外調(diào)查中定位儀獲取的各單木坐標（x，y）為初始聚類中心。

（2）擬合單木圓柱。將初始聚類中心作為圓柱底面圓的中心點，根據(jù)實際樹高情況確定圓柱高為 25 m（過低會導致樹高被低估，過高則不能避免異常點），經(jīng)實驗比較發(fā)現(xiàn)底面圓半徑為 r=1.5 m 時樹高估測的效果最好，在此基礎(chǔ)上擬合單木圓柱。

（3）計算歐式距離并進行聚類。擬合單木圓柱后將所有激光雷達數(shù)據(jù)投影到底面圓所在的二維平面并計算各個點到底面圓中心點的歐氏距離，距離不大于設(shè)定半徑的所有點歸為一類，按此方法依次聚類出所有單木，并獲取聚類后的單木點云數(shù)據(jù)。其中，歐式距離的計算公式：

D=SQRT （（x1-x2）2 +（y1-y2）2） ?。 （1）

式中：D表示歐式距離;x1，x2表示兩點的橫坐標值;y1，y2是兩點的縱坐標值。

（4）提取估測樹高。將聚類后各圓柱體內(nèi) Z 值最高的點視為單木樹冠頂點，其相對高程值視為樹高。

2.4?相關(guān)性分析

本研究基于改進 K-means 算法得到估測樹高后與實測樹高進行相關(guān)性分析。其中估測樹高與實測樹高的相關(guān)程度用 R2 評價，其值越大，說明相關(guān)程度越高[18]。

3?結(jié)果與分析

3.1?點云數(shù)據(jù)預(yù)處理結(jié)果

對激光點云進行預(yù)處理，首先對點云進行去噪處理，本研究采用 Terrosolid 軟件對點云進行去噪處理，采用 Isolated points 算法判斷孤立點，假設(shè)當前點為中心的 25 m 半徑的球形空間內(nèi)沒有其它點的情況下，判斷當前點為噪聲點，最后在軟件中手動去噪。去噪后進行分類處理，本研究采用不規(guī)則三角網(wǎng)格法，分類出地面點和非地面點。接著對點云進行歸一化處理，得到預(yù)處理后的點云數(shù)據(jù)，結(jié)果如圖 3 所示。

3.2?改進 K-means 算法聚類結(jié)果

將預(yù)處理后的 las 點云數(shù)據(jù)按樣方導入到 matlab 中運行得到樣本數(shù)據(jù)，將野外實測數(shù)據(jù)按各樣方整理得到聚類中心數(shù)據(jù)，將樣本數(shù)據(jù)和聚類中心數(shù)據(jù)分別導入，按照改進的 K-means 算法運行依次得到各單木估測樹高和聚類后點云數(shù)據(jù)，運行過程中生成擬合單木圓柱體，以第 8 個樣方為例，結(jié)果如圖 4 所示。

3.3?相關(guān)性分析結(jié)果

經(jīng)改進的 K-means 算法聚類后得到單木估測樹高，將實測數(shù)據(jù)與激光雷達的估測數(shù)據(jù)進行相關(guān)程度分析，將結(jié)果記錄見表 2。

4?結(jié)論與討論

目前基于機載激光雷達技術(shù)獲取林分平均樹高技術(shù)已成熟，但從識別單木到獲取森林單木信息仍有很大的進步空間，故本研究基于改進 K-means 算法在已知單木坐標（x，y）前提下進行提取單木樹高實驗，結(jié)合實驗過程可得到以下結(jié)論。

（1）本研究使用野外數(shù)據(jù)時要注意將經(jīng)緯度轉(zhuǎn)化為大地坐標，與點云數(shù)據(jù)值對應(yīng)，且野外實測數(shù)據(jù)中有少數(shù)的枯死木和叢生木、樹冠枯死木，以及低矮灌木在整理數(shù)據(jù)的過程中應(yīng)將其作為噪聲點剔除，避免產(chǎn)生較大誤差，聚類過程中需擬合單木圓柱，圓柱底面圓半徑會影響到樹高估測的精度，半徑設(shè)定較小會使估測結(jié)果偏低，半徑較大會錯誤聚類，這些因素都影響到本研究的精度。

（2）本次研究表明估測樹高與實測樹高相關(guān)系數(shù) R2 為 0.852，Morsdorf 等[13]基于 K-means 算法實現(xiàn)單木分割，其樹高估測精度 RMSE=0.600;趙晨陽等[13]采用 K-means 算法估測樹高精度為 R2=0.892;說明基于改進的 K-means 算法和原始 K-means 算法都適用于單木樹高提取，但本研究沒有取得較高的精度主要有兩個原因，首先半徑的設(shè)定不靈活，因為天然林分中樹木生長不均一，有高大樹木的也有低矮的樹木，其對應(yīng)的半徑設(shè)定也應(yīng)該有區(qū)別，如設(shè)定動態(tài)半徑（高于 10 m 的樹木r=1.5，低于 10 m 的 r=1.2）樹高估測的精度應(yīng)有所提高;其次本研究不需要已知的z坐標，只利用 x、y 坐標就可以實現(xiàn)樹高的估測，與原始 K-means 算法相比減少了高程參數(shù)，提高了聚類速度，精度沒有得到較大提高。

（3）目前雖然機載激光雷達技術(shù)已經(jīng)廣泛應(yīng)用于林業(yè)的研究，但實驗的高成本及人力的大量輸出仍是面臨的主要問題，本研究可以實現(xiàn)在已知 x、y坐標的情況下較準確的獲取樹高而不需人力測量，可以避免測量過程中的大量人力輸出，同時提高機載激光雷達獲取單木樹高的精度，且基于單木水平的樹高估測相對于整體林分平均樹高的估測更符合實際生產(chǎn)意義。

（4）隨激光雷達技術(shù)的不斷發(fā)展，大范圍獲取林分樹高的技術(shù)已逐漸成熟，但單木分割的相關(guān)研究表明，在提取樹高的過程中目前的方法都有各自的適用范圍，不能夠推廣，且單木水平的森林結(jié)構(gòu)參數(shù)相對于林分更有實用價值，故為了能獲取更高的精度，降低成本，提出新的普遍適用的聚類算法應(yīng)是未來發(fā)展的重心。

【參?考?文?獻】

[1]吳嬌嬌，張亞紅，楊凱博，等.機載激光雷達在林業(yè)中的應(yīng)用[J].安徽農(nóng)業(yè)科學，2016，44（35）：209-212.

WU J J， ZHANG Y H， YANG K B， et al. Application of airborn LiDAR in forestry[J]. Journal of Anhui Agricultural Sciences， 2016， 44（35）： 209-212.

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

XING Y Q， YAO S T， LI M Y， et al. Estimation algorithm of forest aboveground biomass based on airborne full-waveform LiDAR data[J]. Forest Engineering， 2017， 33（4）：21-26.

[3]王曉林，郭斌.柞樹樹高與胸徑相關(guān)關(guān)系的研究[J].森林工程，2012，28（6）：18-21.

WANG X L， GUO B. Study on the correlation between height and diameter at breast height for Quercus mongolica[J]. Forest Engineering， 2012， 28（6）：18-21.

[4]李善堯.樟子松天然林樹高-胸徑模型及胸徑分布規(guī)律分析[J].林業(yè)科技，2018，43（2）：10-13.

LI S R. Analysis of height-diameter relationship and diameter distribution of pinus sylvestris var. mongolica nature forest[J].Forestry Science & Technology，2018，43（2）：10-13.

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

YOU H T， XING Y Q， WANG M， et al. The progress on estimating forest biomass using small footprint LiDAR data[J]. Forest Engineering， 2014， 30（3）：39-42.

[6]邢艷秋，尤號田，霍達，等.小光斑激光雷達數(shù)據(jù)估測森林樹高研究進展[J].世界林業(yè)研究，2014，27（2）：29-34.

XING Y Q， YOU H T， HUO D， et al. Research progress in estimating forest tree height using small footprint lidar data[J]. World Forestry Research， 2014，27 （2）： 29-34.

[7]周淑芳.基于機載LiDAR與航空像片的單木樹高提取研究[D].哈爾濱：東北林業(yè)大學，2007.

ZHOU S F. Research on single tree height extraction based on airborne LiDAR and aerial photograph[D]. Harbin： Northeast Forestry University， 2007.

[8]劉群.基于小光斑機載LiDAR數(shù)據(jù)的單木三維分割[D].北京：北京林業(yè)大學，2016.

LIU Q. Single tree 3D segmentation based on small spot airborne LiDAR data[D]. Beijing： Beijing Forestry University， 2016.

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

LI Z Y， LIU Q W， PANG Y. Review on forest parameters inversion using LiDAR[J]. Journal of Remote Sensing， 2016， 20（5）： 1138-1150.

[10]WANG Y S， WEINACKER H， KOCH B， et al. Lidar point cloud based fully automatic 3D single tree modelling in forest and evaluations of the procedure[J]. The International Archives of the Photogrammetry， Remote Sensing and Spatial Information Sciences， 2008：45-52.

[11]BRANDTBERG T， WARNER T A， LANDENBERGER R E， et al. Detection and analysis of individual leaf-off tree crowns in small footprint， high sampling density lidar data from the eastern deciduous forest in North America[J]. Remote Sensing of Environment， 2003， 85（3）： 290-303.

[12]REITBERGER J， SCHNORR 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.

[13]MORSDORF F， MEIER E， ALLGOWER B， et al. Clustering in airborne laser scanning raw data for segmentation of single trees[J]. International Archives of the Photogrammetry Remote Sensing & Spatial Information Sciences， 2003： 27-33.

[14]趙晨陽，邢艷秋，霍達，等.基于機載LiDAR落葉松樹冠幾何形狀三維重建關(guān)鍵技術(shù)研究[J].西北林學院學報，2015，30（2）：186-190.

ZHAO C Y， XING Y Q， HUO D， et al. Key technologies of three-dimensional geometric reconstruction of larch crown shape based on LiDAR data[J]. Journal of Northwest Forestry University， 2015， 30（2）： 186-190.

[15]WANG H， ZHANG Y， LI P， et al. A method of deriving DEM from airborne LiDAR points cloud data[C]. Proceedings of the Urban Remote Sensing Event， 2013：13-16.

[16]毛建華，何挺，曾齊紅，等.基于TIN的LIDAR點云過濾算法[J].激光雜志，2007，28（6）：36-38.

MAO J H， HE T， ZENG Q H， et al. A filtering algorithm for LIDAR point cloud with TIN[J]. Laser Journal， 2007， 28（6）： 36-38.

[17]焦義濤，邢艷秋，霍達，等.基于機載LiDAR點云估測林分的平均樹高[J].西北林學院學報，2015，30（3）：170-174.

JIAO Y T， XING Y Q， HUO D， et al. Study on mean canopy height estimation from airborne LiDAR point cloud data[J]. Journal of Northwest Forestry University， 2015， 30（3）： 170-174.

[18]CHATTERJEE S， PRICE B. Regression analysis by example[M]. New York： J. Wiley， 2003.