基于多層K-means在森林點(diǎn)云中的單木識(shí)別算法

谷志新,裴方睿

(東北林業(yè)大學(xué) 信息與計(jì)算機(jī)工程學(xué)院,哈爾濱150040)

森林是世界上最重要的生態(tài)環(huán)境之一,對(duì)改善全球氣候及生態(tài)環(huán)境起著重要的作用[1-2],而森林結(jié)構(gòu)參數(shù)的準(zhǔn)確測(cè)算對(duì)于林業(yè)資源的保護(hù)、精細(xì)管理與研究有著重要的意義[3-6]。森林結(jié)構(gòu)的測(cè)算可以依靠彩色圖像數(shù)據(jù)[7]或激光雷達(dá)[8]采集的點(diǎn)云數(shù)據(jù),相對(duì)于彩色圖像,點(diǎn)云數(shù)據(jù)可以更好地描述物體的幾何信息。

激光雷達(dá)以激光束照射目標(biāo),物體表面將收到的激光束返射回傳感器,通過(guò)回波脈沖和發(fā)射激光束之間的延遲得到采樣點(diǎn)與傳感器之間的距離,繼而可以實(shí)現(xiàn)對(duì)森林結(jié)構(gòu)參數(shù)的估測(cè)[9-13]。單木識(shí)別是估測(cè)單木結(jié)構(gòu)參數(shù)的前提,現(xiàn)在已有多種算法能夠成功提取森林點(diǎn)云中的單木。根據(jù)點(diǎn)云數(shù)據(jù)處理方式的不同,可將已有的單木識(shí)別算法分為二維圖像搜索處理和三維空間位置處理兩類。其中,基于冠層高度模型(Canopy Height Model,CHM)識(shí)別單木是二維圖像搜索處理數(shù)據(jù)中使用比較廣泛的方法,森林冠層高度模型是一個(gè)表達(dá)植被距離地面高度的模型[14]。劉魯霞等[15]對(duì)點(diǎn)云數(shù)據(jù)進(jìn)行分層?xùn)鸥窕苫叶葓D,通過(guò)Hough變換算法識(shí)別柵格圖像圖中的近似圓來(lái)確定單木位置；杜意鴻等[16]通過(guò)將分水嶺算法分割后的油松冠狀多邊形與森林冠層高度模型重疊,搜索局部最大值,濾除樹冠半徑小于指定閾值的偽樹冠頂點(diǎn),得到最終的樹冠頂點(diǎn)。但是,將點(diǎn)云數(shù)據(jù)轉(zhuǎn)換為二維后,無(wú)法識(shí)別樹干點(diǎn)與非樹干點(diǎn),造成Hough變換會(huì)在非樹干點(diǎn)處識(shí)別出圓,因此,將點(diǎn)云轉(zhuǎn)為二維數(shù)據(jù)的算法受到噪點(diǎn)的影響更大。三維空間位置處理相比較二維圖像搜索處理更好地保存了點(diǎn)云的幾何信息,對(duì)點(diǎn)云信息提取更精確[17]。Li等[18]通過(guò)利用點(diǎn)之間的相對(duì)間距以及通過(guò)定義適當(dāng)?shù)拈g隔閾值,從點(diǎn)云中依次識(shí)別樹干;孫拱等[19]先擬合圓柱作為聚類中心,再計(jì)算點(diǎn)與聚類中心的歐式距離決定是否為同一樹木的點(diǎn)。與二維圖像搜索相比,基于三維空間位置處理數(shù)據(jù)的方式具有更好的抗噪性,更適合處理點(diǎn)云,但是在點(diǎn)云聚類的過(guò)程中,需要對(duì)點(diǎn)云數(shù)據(jù)逐點(diǎn)分析,聚類的閾值難以選取,導(dǎo)致樹干點(diǎn)的識(shí)別精度較低,并且算法效率較低。

本研究以孟家崗林場(chǎng)落葉松人工林作為研究對(duì)象,在激光雷達(dá)掃描點(diǎn)云數(shù)據(jù)的基礎(chǔ)上提出了一種多層K-means聚類的方式,提高原有算法的單木識(shí)別精度以及提高算法的效率。

1 研究區(qū)域概況與數(shù)據(jù)預(yù)處理

1.1 研究區(qū)域概況

研究區(qū)位于樺南縣北部的孟家崗林場(chǎng),地理坐標(biāo)為北緯46°20′16″～46°30′50″,東經(jīng)130°32′42″～130°52′36″,隸屬佳木斯市林業(yè)局,行政區(qū)域?qū)贅迥峡h境內(nèi)。林場(chǎng)地處完達(dá)山西麓余脈,以低山丘陵為主,坡度較為平緩,大部分坡度在10～20°之間。森林總面積約16 000km2,人工林始建于1956年,以落葉松(Larixgmelinii)、樟子松(Pinussylvestris)等針葉人工林為主,面積約占整個(gè)林分面積的2/3,其它為天然次生林,約占1/3,森林覆蓋率達(dá)81.7%。孟家崗林場(chǎng)衛(wèi)星地圖如圖1所示。

圖1 孟家崗林場(chǎng)衛(wèi)星地圖

1.2 數(shù)據(jù)采集與數(shù)據(jù)概要

1.2.1背包激光雷達(dá)采集數(shù)據(jù)

本研究于2021年5月,使用背包激光雷達(dá)對(duì)實(shí)驗(yàn)樣地進(jìn)行掃描。當(dāng)日天氣晴,適合對(duì)數(shù)據(jù)進(jìn)行采集。采集過(guò)程中以樣地邊緣一處作為起點(diǎn),采集人員持設(shè)備以“S”型對(duì)樣地進(jìn)行掃描,最后從邊緣處回到起點(diǎn)使掃描路線閉合。掃描所使用的設(shè)備參數(shù)如表1所示。

表1 背包激光雷達(dá)參數(shù)

1.2.2數(shù)據(jù)概要

實(shí)驗(yàn)樣地是以落葉松構(gòu)成的單層林,占地面積約760m2,樣地內(nèi)包含649株落葉松,郁閉度為0.6。經(jīng)過(guò)背包雷達(dá)掃描后原始數(shù)據(jù)包含14 007 171個(gè)點(diǎn),通過(guò)計(jì)算各點(diǎn)的平均最小距離值來(lái)計(jì)算點(diǎn)云的密度,計(jì)算公式如下。

(1)

1.3 數(shù)據(jù)預(yù)處理

在森林點(diǎn)云數(shù)據(jù)中包含大量與樹干無(wú)關(guān)的噪點(diǎn),過(guò)多的噪點(diǎn)會(huì)明顯影響算法的效率與正確性,因此,需要對(duì)數(shù)據(jù)進(jìn)行去噪操作。

1.3.1地面噪點(diǎn)的去除

地面噪點(diǎn)是森林點(diǎn)云噪點(diǎn)中占比最多的噪點(diǎn)類型,本文采用RANSAC[20-21]算法進(jìn)行地面點(diǎn)云數(shù)據(jù)的去除。首先,隨機(jī)選取確定模型參數(shù)所需的最小點(diǎn)數(shù)求解模型參數(shù)；然后,去除大于預(yù)定閾值的點(diǎn)。如果小于閾值的點(diǎn)占集合中總點(diǎn)數(shù)的比例接近預(yù)定的閾值,則更新模型,重復(fù)N次保留最優(yōu)結(jié)果。

1.3.2離群點(diǎn)的去除

激光雷達(dá)采集過(guò)程中會(huì)由于特殊原因采集到非樹干且非地面的噪點(diǎn),如飛行物等,這類噪點(diǎn)的密度遠(yuǎn)小于樹干點(diǎn),本文采用半徑離群值去除該類噪點(diǎn),刪除在給定球形范圍中密度小于閾值的點(diǎn)。

2 單木識(shí)別方法

2.1 點(diǎn)云水平分層

點(diǎn)云水平分層是為了解決森林點(diǎn)云數(shù)據(jù)中存在兩個(gè)問(wèn)題:

第一,點(diǎn)云數(shù)據(jù)的數(shù)據(jù)量過(guò)于龐大,而聚類中每次迭代都需要重新計(jì)算聚類中心,并且將點(diǎn)云數(shù)據(jù)從2類聚合為1類的時(shí)間也遠(yuǎn)大于從4類聚合為2類,以單木聚類為例,每次迭代的時(shí)間如圖2所示。

圖2 算法收斂時(shí)間與迭代次數(shù)的關(guān)系

其中,測(cè)試所用的點(diǎn)云數(shù)據(jù)包含9 418個(gè)點(diǎn),由圖2可以看出,隨著迭代次數(shù)的增加,每次算法的收斂時(shí)間也會(huì)增加。

第二,聚類中心的計(jì)算依賴于局部最大值,只對(duì)整體提取局部最大值會(huì)導(dǎo)致單木點(diǎn)云的過(guò)聚類與漏聚類。對(duì)點(diǎn)云進(jìn)行分層計(jì)算,可以分別計(jì)算出各切片層的局部最大值,提高聚類的正確率,各個(gè)切片層之間的聚類結(jié)果也不會(huì)互相影響。當(dāng)各層聚類完成后,以上下兩層點(diǎn)云簇中心位置最近為合并條件,第一層與第二層合并,重新計(jì)算點(diǎn)云簇中心,再與第三層完成合并,重復(fù)此操作直至全部切片合并完成。

切片的層數(shù)可根據(jù)物體的形態(tài),點(diǎn)云間的密度決定[16]。

2.2 多層K-means聚類識(shí)別單木

K-means是一種基于距離評(píng)價(jià)兩個(gè)對(duì)象的相異度算法,傳統(tǒng)的K-means算法僅依靠?jī)蓚€(gè)對(duì)象的距離來(lái)判斷相異度,這種方法在森林點(diǎn)云數(shù)據(jù)的處理中存在問(wèn)題:森林中樹木之間存在大量的交叉,無(wú)法準(zhǔn)確識(shí)別樹木之間的交叉部分,為了解決這個(gè)問(wèn)題,本研究在傳統(tǒng)的K-means算法對(duì)相異度的評(píng)價(jià)上進(jìn)行了改進(jìn)。

2.2.1相異度M的計(jì)算

相異度是判斷點(diǎn)云簇之間是否可以合并的標(biāo)準(zhǔn),當(dāng)兩點(diǎn)云簇相異度低時(shí),被視為同一物體的點(diǎn)云,在大多數(shù)使用K-means進(jìn)行單木識(shí)別的研究中,點(diǎn)云簇之間的相異度計(jì)算只依靠點(diǎn)云簇之間的歐式距離[17-19],這樣造成了屬于多株單木的點(diǎn)存在交叉時(shí)錯(cuò)誤聚類,本研究在計(jì)算相異度中引入了方向項(xiàng)作為另一個(gè)影響因子[22],這在一定程度上降低了異類點(diǎn)云交叉造成的錯(cuò)誤聚類。

1)距離因子的計(jì)算。計(jì)算兩個(gè)集合間的距離作為兩集合之間的距離因子S,具體公式為:

S=d(pi,pj)

(2)

式中:pi,pj為兩集合間歐式距離(Euclidean Distance)最小的點(diǎn),d(pi,pj)為兩點(diǎn)之間歐式距離,在三維空間中,歐式距離公式如下:

d(pi,pj)=sqrt((xi-xj)2+(yi-yj)2+(zi-zj)2)

(3)

式中:xi,xj為pi,pj的橫軸坐標(biāo);yi,yj為pi,pj的縱軸坐標(biāo);zi,zj為pi,pj豎軸坐標(biāo)。

2)方向因子的計(jì)算。以法向量為基礎(chǔ)計(jì)算方向因子T,本研究中,法向量的計(jì)算基于K近鄰估計(jì),由某一點(diǎn)及其歐式距離最短的k個(gè)點(diǎn)估算平面,過(guò)該點(diǎn)作平面的垂線即為該點(diǎn)的法向量,k的數(shù)值由點(diǎn)云的密度決定。方向因子的計(jì)算如公式4所示:

(4)

式中:vi,vj分別為點(diǎn)pi,pj的法向量。

pi,pj的選擇與計(jì)算距離因子中的點(diǎn)相同,vi,vj夾角計(jì)算公式如公式5所示:

(5)

式中:xi,xj為vi,vj的橫軸坐標(biāo)；yi,yj為vi,vj的縱軸坐標(biāo)；zi,zj為vi,vj豎軸坐標(biāo)。

最終,兩點(diǎn)云集合相異度M計(jì)算公式如公式6所示:

M=S+ε×T

(6)

式中:ε為待定系數(shù),用于方向因子與距離因子的平衡,當(dāng)相異度M小于預(yù)定閾值時(shí),點(diǎn)云集合視為相似點(diǎn)云簇,可以合并。

2.2.2點(diǎn)云局部最大值的計(jì)算

本研究區(qū)域中,以落葉松為主,在形態(tài)上接近圓錐體,對(duì)原始點(diǎn)云數(shù)據(jù)去噪后樹干位置大致位于點(diǎn)云垂直投影的中心位置,在點(diǎn)云分層后,通過(guò)imregionalmax函數(shù)可以判斷各層次在中心位置上的局部最大值點(diǎn),其原理是找出當(dāng)前區(qū)域內(nèi)連通性最強(qiáng)的點(diǎn),即找出在中心位置上的點(diǎn),并且該點(diǎn)在近鄰的點(diǎn)中連通性最強(qiáng)。

2.2.3多層K-means聚類流程

多層K-means聚類步驟如下:1)點(diǎn)云橫向切片；2)分層輸入點(diǎn)云集合；3)將點(diǎn)云集合中的每一個(gè)點(diǎn)初始化為一個(gè)集合,即P={p1,p2,p3,…,pn}初始化為P={{p1},{p2},{p3},…,{pn}}；4)初始化空集合C,視為已合并區(qū),P為待合并區(qū),計(jì)算P中相鄰集合的相異度M,將P中相異度小于合并閾值的集合取出合并成新的集合放入集合C中,直至P為空時(shí)該步驟結(jié)束；5)將C中的元素重新放入P中重新計(jì)算相異度,合并后放入C中；6)重復(fù)“4)”-“5)”步驟直至函數(shù)收斂,其中函數(shù)收斂的條件為最后一次聚類輸出的C與上一次聚類的輸出C完全相等；7)不同層次點(diǎn)云合并。流程如圖3所示。

圖3 多層K-means流程

3 實(shí)驗(yàn)與分析

本文使用Python與開源庫(kù)Open3D進(jìn)行實(shí)驗(yàn),點(diǎn)的法向量以15個(gè)近鄰點(diǎn)計(jì)算,相異度的計(jì)算中,距離因子與方向因子的平衡系數(shù)ε的值為0.3。

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

本文使用RANSAC算法和半徑離群值去除地面噪點(diǎn)以及其他非地面且非樹干的噪點(diǎn),根據(jù)點(diǎn)云形態(tài),半徑離群值去噪的方式為在半徑為2單位距離的球體中,點(diǎn)云密度小于20的點(diǎn)為離群點(diǎn)。原始點(diǎn)云、去除地面噪點(diǎn)、去除離群噪點(diǎn)分別為圖4—圖6所示。

在圖5中,黑色部分為地面點(diǎn)云,紅色為非地面點(diǎn)云,在圖6中,紅色為未聚類的單木點(diǎn)云,藍(lán)色部分為去除的離群點(diǎn)。

圖4 原始點(diǎn)云

圖5 去除地面噪點(diǎn)

圖6 去除離群點(diǎn)

3.2 多層K-means聚類結(jié)果

改變聚類的層數(shù)會(huì)對(duì)識(shí)別單木的精度產(chǎn)生影響,如果分層較多會(huì)使單層點(diǎn)云密度小,聚類后的單木形態(tài)不準(zhǔn)確,如果層數(shù)少則會(huì)影響局部最大值的位置,當(dāng)一個(gè)點(diǎn)云簇距離該層的局部最大值較遠(yuǎn)時(shí)，會(huì)造成相異度較大,導(dǎo)致本屬于該樹木的點(diǎn)云無(wú)法聚類。 以單株樹木聚類為例,使用不同層數(shù)進(jìn)行對(duì)比實(shí)驗(yàn),其中測(cè)試使用的單木點(diǎn)云包含9 418個(gè)點(diǎn)。測(cè)試結(jié)果如表2所示。

表2 不同層數(shù)對(duì)識(shí)別單木精度的影響

由表2可知,當(dāng)層數(shù)為5時(shí),單木的精度最高。

同時(shí),不同的層數(shù)對(duì)算法的識(shí)別率也有影響,不同層數(shù)在原始數(shù)據(jù)去噪后的單木識(shí)別率如表3所示。

由表3可知,當(dāng)層數(shù)為7時(shí),算法的識(shí)別率最高,但是在考慮單木聚類精度的情況下,本文選擇將點(diǎn)云分為5層進(jìn)行實(shí)驗(yàn)。

在計(jì)算結(jié)束后,將所有聚類結(jié)果展示在三維空間中,聚類的效果如圖7所示。

表3 不同層數(shù)對(duì)樹木識(shí)別率的影響

圖7 點(diǎn)云聚類結(jié)果

實(shí)驗(yàn)樣地共有649株樹木,經(jīng)過(guò)與原始點(diǎn)云的比對(duì),其中,未能被識(shí)別的單木有32株,將不同單木識(shí)別為同一單木的情況有4類,共有59株單木識(shí)別失敗。未能被識(shí)別的單木與不同單木識(shí)別為同一種單木的案例如圖8—圖9所示。最終將被識(shí)別樹木與實(shí)驗(yàn)樣地調(diào)查的數(shù)據(jù)比對(duì),未發(fā)現(xiàn)將其他物體識(shí)別為樹木的情況。

圖8 未能被識(shí)別的單木

圖9 將多株單木識(shí)別為同一單木

結(jié)果的統(tǒng)計(jì)如表4所示。

表4 實(shí)驗(yàn)結(jié)果的統(tǒng)計(jì)

經(jīng)過(guò)觀察原始點(diǎn)云,未能識(shí)別的單木主要原因在于被其他樹木遮擋嚴(yán)重而導(dǎo)致點(diǎn)云稀疏,經(jīng)過(guò)半徑離群值去噪后點(diǎn)云數(shù)量過(guò)少,點(diǎn)與點(diǎn)間距離過(guò)大,導(dǎo)致無(wú)法被聚類。而被聚類成同一單木的多株單木，是由于樹木之間交叉過(guò)于密集造成的。

3.3 多層K-means與單層K-means對(duì)比

在保持點(diǎn)云數(shù)據(jù)點(diǎn)的數(shù)量一致、算法參數(shù)一致的情況下,分別使用多層K-means與單層K-means進(jìn)行聚類,算法的收斂時(shí)間如圖10所示。

圖10 多層K-means與單層K-means收斂時(shí)間對(duì)比

在同一數(shù)據(jù)集下,多層K-means與單層K-means的單木識(shí)別識(shí)別率對(duì)比如表5所示。

表5 同一數(shù)據(jù)集中多層K-means與單層K-means的效果對(duì)比

由圖10和表5可知,多層K-means從收斂時(shí)間與聚類效果都優(yōu)于傳統(tǒng)的單層K-means。

4 結(jié)束語(yǔ)

以點(diǎn)云數(shù)據(jù)為基礎(chǔ)獲取森林信息,是對(duì)森林結(jié)構(gòu)準(zhǔn)確測(cè)算的重要手段之一[23],準(zhǔn)確的單木識(shí)別是這些研究的基礎(chǔ),本算法在傳統(tǒng)K-means聚類的基礎(chǔ)上進(jìn)行了分層處理,添加了K-means聚類條件,降低算法的收斂時(shí)間,提高單木點(diǎn)云聚類的準(zhǔn)確率,得出結(jié)論如下:

1)對(duì)點(diǎn)云進(jìn)行切片操作可以提升聚類的準(zhǔn)確性,降低算法的收斂時(shí)間,切片操作大幅度減少了單次聚類所需要計(jì)算的點(diǎn),計(jì)算后只需對(duì)不同層的同一類點(diǎn)云拼接就可以得到單木的點(diǎn)云。由于聚類算法的特性,聚類中的一次錯(cuò)誤聚類就會(huì)導(dǎo)致聚類中心的計(jì)算錯(cuò)誤,從而導(dǎo)致了后續(xù)繼續(xù)依照錯(cuò)誤的中心進(jìn)行聚類,導(dǎo)致聚類錯(cuò)誤的點(diǎn)數(shù)增加,而對(duì)點(diǎn)云進(jìn)行分層,可以分別計(jì)算不同層次的點(diǎn)云中心,避免了不同層次間聚類中心互相影響。

2)在相異度的計(jì)算中,使用以法向量為基礎(chǔ)計(jì)算方向因子可以更充分地考慮點(diǎn)云中不同點(diǎn)云簇間的關(guān)系,在一定程度上解決了不同單木之間交叉的問(wèn)題。

3)本聚類算法不進(jìn)行圓柱擬合等操作,因?yàn)閷?duì)于較為稀疏的點(diǎn)云數(shù)據(jù)來(lái)說(shuō),點(diǎn)云的缺失會(huì)導(dǎo)致樹干部分無(wú)法擬合圓柱而導(dǎo)致漏檢測(cè)的情況,本算法逐點(diǎn)進(jìn)行計(jì)算,確保所有點(diǎn)均被考慮,不被樹干形態(tài)所影響。

4)森林結(jié)構(gòu)是森林經(jīng)營(yíng)管理的重要依據(jù),精準(zhǔn)的森林結(jié)構(gòu)參數(shù)可以使相關(guān)從業(yè)人員判斷該森林中單木的情況,從而提高林分質(zhì)量,充分發(fā)揮林地生產(chǎn)力。