王偉偉，龐勇，杜黎明，張鐘軍，梁曉軍

1.中國林業(yè)科學(xué)研究院 資源信息研究所,北京 100091;2.國家林業(yè)和草原局 林業(yè)遙感與信息技術(shù)重點(diǎn)實(shí)驗(yàn)室,北京 100091;3.北京師范大學(xué) 人工智能學(xué)院,北京 100875

1 引 言

激光雷達(dá)LiDAR（Light Detection And Ranging）可以獲得高精度的森林空間結(jié)構(gòu)信息，在森林資源調(diào)查監(jiān)測中的應(yīng)用日益廣泛（Hyypp? 等，2008）。機(jī)載激光雷達(dá)ALS（Airborne Laser Scanning）可以精確刻畫冠層垂直結(jié)構(gòu)，提取林分和單木尺度的森林參數(shù)（Sa?kov 等，2017；李增元 等，2016）。ALS點(diǎn)云數(shù)據(jù)的單木分割對林業(yè)科學(xué)研究和生產(chǎn)實(shí)踐具有重要意義，基于準(zhǔn)確的分割結(jié)果，可以從中獲取描述單木空間結(jié)構(gòu)特征和生物化學(xué)組分屬性的單木因子，為評價(jià)森林生長和生態(tài)功能、評估森林破壞程度和監(jiān)測森林再生情況提供數(shù)據(jù)支持（Chen等，2007）。

對ALS 點(diǎn)云數(shù)據(jù)的單木分割算法可以大致分為二維方法和三維方法（Lindberg 和Holmgren，2017）。二維方法通過表示冠層上部輪廓的CHM（Canopy Height Model） 或DSM （Digital Surface Model）柵格數(shù)據(jù)，將其中的局部最大值識(shí)別為樹頂來分割單木。搜索局部最大值時(shí)，可以選擇固定或可變的窗口尺寸（Popescu 等，2002；Zimble等，2003）?；诠趯由媳砻娴木植孔畲笾?，常用的單木分割方法包括分水嶺方法（Chen 等，2006；Zhao 等，2014；李旺等，2015），區(qū)域增長法（Solberg等，2006；Pang等，2008；Zhen等，2015），以及其他基于二維圖像分割的算法（Brandtberg，2007；李平昊等，2018）。盡管這些方法對于樹頂明確且樹冠輪廓清晰的高層單木具有很好的分割效果，但使用的二維柵格數(shù)據(jù)損失了幾乎全部的林下信息，對于許多低矮單木缺乏識(shí)別能力。

與二維方法相比，更好地利用完整LiDAR 點(diǎn)云數(shù)據(jù)的三維分割方法在單木分割問題中受到越來越普遍的關(guān)注。聚類方法是一種常用的三維方法，其主要包括K-means 聚類、層次聚類和局部最大值聚類等（Morsdorf 等，2004；Lee 等，2010；Li 等，2012）。Gupta 等（2010）比較了幾種改進(jìn)的K-means聚類和層次聚類方法，發(fā)現(xiàn)在K-means算法中使用局部最大值作為初始種子點(diǎn)并縮小點(diǎn)云高度取得了更好的結(jié)果，同時(shí)層次聚類方法沒有得到令人滿意的結(jié)果。Ayrey 等（2017）提出了一種對點(diǎn)云數(shù)據(jù)分層使用K-means 的層疊方法，對樹冠覆蓋實(shí)現(xiàn)了較高的估計(jì)精度。這類算法對于大樹和單一簡單林分可以產(chǎn)生合理的分割效果，但是算法對于多層復(fù)雜林分的應(yīng)用存在局限性（Williams等，2020）。

基于譜圖理論的譜聚類方法在分割問題中表現(xiàn)優(yōu)異，其對數(shù)據(jù)分布沒有限制，且無需假定初始種子點(diǎn)。譜聚類方法通過構(gòu)造數(shù)據(jù)點(diǎn)間的相似度矩陣并進(jìn)行特征分解，在特征空間中實(shí)現(xiàn)數(shù)據(jù)的聚類分割，然后將結(jié)果映射回原始的數(shù)據(jù)空間。譜聚類方法在LiDAR 點(diǎn)云數(shù)據(jù)處理中的挑戰(zhàn)在于高計(jì)算復(fù)雜度，尤其是相似度矩陣的特征分解步驟。Nystr?m 近似是對譜聚類方法的一種有效簡化（Fowlkes 等，2004），其使用采樣技術(shù)和數(shù)值逼近理論，通過插值權(quán)重將采樣點(diǎn)的特征向量擴(kuò)展到非采樣點(diǎn)，從而快速地得到原始相似度矩陣的近似特征值和特征向量，減少算法的計(jì)算復(fù)雜度。由于不同的采樣結(jié)果會(huì)對原始相似度矩陣產(chǎn)生不同的近似，因此采樣方法對Nystr?m 方法的精確性至關(guān)重要。已有很多研究提出了不同的采樣方法（Zhang 和You，2011；Zeng 等，2014；Cohen 等，2015）。Bouneffouf 和Birol（2015）指出，可以通過最小化新的和已有的采樣點(diǎn)之間相似度的平方和來最大化采樣點(diǎn)間相似度矩陣行列式，從而最小化近似誤差?；谶@一概念，Bouneffouf 和Birol（2015）提出了考慮采樣點(diǎn)方差和相似度的逐步迭代MSSS（Minimum Sum of Squared Similarities）方法，該方法與其他采樣方法在不同數(shù)據(jù)集的對比實(shí)驗(yàn)中展示出很好的精度結(jié)果和優(yōu)異的采樣性能，同時(shí)，理論分析中也進(jìn)一步證實(shí)了該采樣方法的優(yōu)越性（Bouneffouf和Birol，2016）。

此外，考慮到三維點(diǎn)云數(shù)據(jù)量產(chǎn)生的算法效率限制和存儲(chǔ)空間困境，體素化方法成為許多三維方法對原始點(diǎn)云數(shù)據(jù)處理的有效方式。體素化方法將原始點(diǎn)云投影到體素空間中，進(jìn)而在體素空間中進(jìn)行單木分割。Wang 等（2016）根據(jù)點(diǎn)密度使用給定分辨率的立方體來構(gòu)建體素點(diǎn)云，然后根據(jù)體素的空間鄰域關(guān)系從冠層中檢測樹頂。Reitberger 等（2009）對體素化點(diǎn)云使用歸一化切割NCut（Normalized Cut）方法（Shi和Malik，2000；王濮等，2019）進(jìn)行迭代二分實(shí)現(xiàn)了對單木的三維分割，但是其高的計(jì)算復(fù)雜度限制了對大規(guī)模數(shù)據(jù)的普及和推廣。

與傳統(tǒng)規(guī)則立方體的體素化方法相比，超體素方法（李平昊等，2018）提供了一種更為靈活有效的思路來構(gòu)建更為緊致的體素空間?？梢赃x擇合適的分割方法對原始點(diǎn)云進(jìn)行過分割，將獲取的過分割結(jié)果視為體素點(diǎn)，從而構(gòu)造整個(gè)體素空間。本文算法參考超體素的思路，使用均值漂移（mean shift）算法實(shí)現(xiàn)對LiDAR 點(diǎn)云的體素化，有效地減小后續(xù)計(jì)算過程中的數(shù)據(jù)量。作為一種快速有效的聚類算法，mean shift 算法不需要對數(shù)據(jù)分布或聚類數(shù)目進(jìn)行假設(shè)，通過迭代地將每個(gè)搜索點(diǎn)移向偏移均值點(diǎn)處來對點(diǎn)進(jìn)行分組。在此過程中，搜索點(diǎn)總是沿著局部密度遞增的方向移動(dòng)到密度極大值點(diǎn)處，最終構(gòu)造出合理緊致的體素空間。對構(gòu)造的體素空間，本文算法使用Nystr?m 方法實(shí)現(xiàn)譜聚類算法中特征分解步驟的快速近似，得到最終的單木分割結(jié)果。其中，關(guān)鍵的采樣過程通過MSSS 方法完成，以快速獲取最佳采樣結(jié)果。本文算法旨在通過一系列優(yōu)化手段提升譜聚類算法的計(jì)算效率，使其優(yōu)越的分割性能在機(jī)載LiDAR 點(diǎn)云數(shù)據(jù)的單木分割問題中得以發(fā)揮，為大規(guī)模數(shù)據(jù)應(yīng)用提供可能。

2 實(shí)驗(yàn)區(qū)與數(shù)據(jù)

2.1 實(shí)驗(yàn)區(qū)基本情況

研究區(qū)位于中國黑龍江省佳木斯市的孟家崗林場（46°20′N—46°30′N，130°32′E—130°52′E），其地理位置如圖1所示。孟家崗林場中的人工林面積達(dá)1.5 萬ha，占總林地面積的76.7%。人工林的90%由落葉松（Larix olgensisHenry），紅松（Pinus koraiensisSieb.et Zucc.），樟子松（Pinus sylvestrisL.var.mongolica Litv.）和云杉（Picea asperataMast.）組成。

圖1 孟家崗林場位置及不同株數(shù)密度的樣地分布Fig.1 Study area and the distribution of field plots with different stem densities

2.2 LiDAR和地面實(shí)測數(shù)據(jù)

本研究采用中國林業(yè)科學(xué)研究院機(jī)載遙感系統(tǒng)CAF-LiCHy 于2017年5月31日—6月15日 采 集的LiDAR 數(shù)據(jù)。其中LiDAR 傳感器為Riegl 的LMS-Q680i（Pang 等，2016），航攝平臺(tái)為運(yùn)-5 多用途飛機(jī)，平均絕對飛行高度為1000 m，相對飛行高度約650 m。ALS 點(diǎn)云采集于LiCHy 系統(tǒng)，點(diǎn)云數(shù)據(jù)密度約為20 pts/m2。對ALS點(diǎn)云的處理包括基于體素的孤立點(diǎn)去除、基于不規(guī)則加密三角網(wǎng)TIN （Triangulated Irregular Network）的地面點(diǎn)分類、高度歸一化處理和地面點(diǎn)去除，此點(diǎn)云數(shù)據(jù)作為原始點(diǎn)云數(shù)據(jù)。

本研究選擇了11 個(gè)落葉松樣地?cái)?shù)據(jù)進(jìn)行單木分割實(shí)驗(yàn)以評估算法性能，其位置分布如圖1 所示。地面調(diào)查數(shù)據(jù)于2017年6月—7月獲得，表1對樣地的關(guān)鍵參數(shù)進(jìn)行了總結(jié)。胸徑DBH（Diameter at Breast Height）使用圍尺測量得到，且只記錄DBH≥5 cm 的單木參數(shù)。樹高使用超聲波樹木測高儀測得，11 個(gè)樣地的平均樹高從13—25 m不等，且同一個(gè)樣地內(nèi)樹高差異較小。冠幅使用鋼卷尺測量東西及南北兩個(gè)方向，并計(jì)算其平均值作為平均冠幅。樣地中心地理位置由全球衛(wèi)星導(dǎo)航系統(tǒng)GNSS（Global Navigation Satellite System）手持接收機(jī)測量獲取，同時(shí)該接收機(jī)使用參考接收機(jī)獲取的實(shí)時(shí)差分信號(hào)進(jìn)行校正，單木位置通過GNSS 手持接收機(jī)、全站儀、地基激光雷達(dá)和皮尺等測量數(shù)據(jù)的空間校正得到，定位精度在1 m 之內(nèi)（梁曉軍 等，2020）。此外，落葉松人工林樣地單木總體分布相對密集，但樣地之間存在一定株數(shù)密度差異，據(jù)此將11 個(gè)樣地按照株數(shù)密度分為高（H）、中（M）和低（L）等3 類，用以進(jìn)一步檢驗(yàn)算法對不同株數(shù)密度樣地的分割性能。最終，去除樣地內(nèi)的枯死木記錄，得到的單木株數(shù)為991株。

表1 樣地?cái)?shù)據(jù)參數(shù)Table 1 Characteristics of field plots

3 原理與方法

本文基于譜聚類方法的思路框架，并通過引入基于均值漂移（mean shift）的體素化和Nystr?m方法解決其計(jì)算效率低的問題，圖2顯示了本文算法的一般工作流程。傳統(tǒng)譜聚類算法直接對原始數(shù)據(jù)構(gòu)造相似圖，對相似度矩陣進(jìn)行特征分解，在特征空間使用K-means 方法進(jìn)行聚類，最終將結(jié)果映射回原始數(shù)據(jù)空間得到分割結(jié)果。與之相比，本文算法有以下兩個(gè)關(guān)鍵的改進(jìn)之處：

圖2 基于Nystr?m-譜聚類算法的單木分割流程Fig.2 Flowchart of individual tree segmentation using Nystr?m-based spectral clustering algorithm

（1）使用mean shift 方法將原始的LiDAR 點(diǎn)云數(shù)據(jù)轉(zhuǎn)換到體素空間以合理壓縮數(shù)據(jù)量；

（2）使用Nystr?m 方法快速計(jì)算相似度矩陣的近似特征向量和特征值。其中的關(guān)鍵采樣步驟由考慮采樣點(diǎn)方差和相似度的逐步迭代MSSS 采樣方法實(shí)現(xiàn)。

3.1 基于mean shift的體素化

Mean shift 是一種快速有效的聚類算法，本文使用mean shift 方法來完成體素化過程，并將mean shift的帶寬參數(shù)設(shè)置得較小，以實(shí)現(xiàn)過分割來達(dá)到體素化目的（帶寬參數(shù)的詳細(xì)設(shè)置見3.5 節(jié)中關(guān)于參數(shù)配置部分）。與使用給定分辨率的立方體來構(gòu)建規(guī)則體素點(diǎn)云的方法不同，在使用mean shift 方法構(gòu)造的體素空間中，每個(gè)體素點(diǎn)由mean shift 的聚類結(jié)果表示，體素位置由類內(nèi)中心點(diǎn)坐標(biāo)確定，體素權(quán)重等于類內(nèi)的點(diǎn)數(shù)。后續(xù)的分割過程將基于體素空間完成，可以通過mean shift 聚類結(jié)果的類別標(biāo)簽將原始點(diǎn)云中的每個(gè)點(diǎn)對應(yīng)到體素點(diǎn)云中，從而將體素空間中得到的分類結(jié)果映射回原始點(diǎn)云，最終得到單木的分割結(jié)果。同時(shí)，可以通過體素權(quán)重將原始點(diǎn)云的空間分布信息表征在體素空間中，即，體素權(quán)重越大表明該局部區(qū)域的原始點(diǎn)密度越高。在實(shí)驗(yàn)過程中，使用mean shift聚類算法得到的體素點(diǎn)云數(shù)據(jù)量較原始點(diǎn)云約減少了10 倍，大大減輕了算法在后續(xù)過程中的計(jì)算負(fù)擔(dān)。

3.2 相似度函數(shù)與相似圖

參考Von Luxburg（2007）的工作，本文使用常用的高斯相似度函數(shù)構(gòu)造k-近鄰圖，其很容易構(gòu)造稀疏鄰接矩陣，且不易受到不良參數(shù)選擇的影響?？紤]單木的形態(tài)和體素的屬性，本文采用的高斯相似度函數(shù)構(gòu)造如下：

式中，s(xi，xj)表示兩個(gè)體素xi和xj之間的相似度；KNN 表示k-近鄰；和別是兩個(gè)體素xi和xj之間的水平和垂直歐氏距離，其被分配不同的比例因子σxy和σz以適用于單木的橢球形狀；兩個(gè)加權(quán)因子ni和nj是體素xi和xj的權(quán)重，以維持體素空間與原始點(diǎn)云的一致性。此外，考慮到基于mean shift的體素化產(chǎn)生的體素點(diǎn)沒有規(guī)則的空間位置索引，故無法像傳統(tǒng)立方體體素化方法那樣快速訪問k-近鄰，本文采用KD-Tree數(shù)據(jù)結(jié)構(gòu)來彌補(bǔ)這一不足。

3.3 Nystr?m方法

在體素空間中構(gòu)建完成k-近鄰圖即可得到體素點(diǎn)間的相似度矩陣，本文使用Nystr?m 方法快速計(jì)算相似度矩陣的近似特征向量和特征值。

3.3.1 Nystr?m理論

對一有N個(gè)體素點(diǎn)的數(shù)據(jù)集，構(gòu)造其k-近鄰圖得到相似度矩陣W∈RN×N，Nystr?m 方法將所有的點(diǎn)分為n個(gè)采樣點(diǎn)和m個(gè)剩余點(diǎn)(m=N-n)，進(jìn)而將相似度矩陣W分塊為

式中，A∈Rn×n表示采樣點(diǎn)間的相似度矩陣，且有對角化形式A=UΛUT；B∈Rn×m表示采樣點(diǎn)和剩余點(diǎn)之間的相似度矩陣；C∈Rm×m表示剩余點(diǎn)間的相似度矩陣。

Fowlkes 等（2004）給出了W的近似正交特征向量的構(gòu)造方法。若矩陣A正定，定義S=A+A-12BBTA-12，且其對角化為S=USΛS，其中，US∈Rn×n且列元素為為S的特征向量，ΛS是對角元素為S對應(yīng)特征值的對角矩陣。那么=VΛSVT即為W的近似特征分解，且有

通過這種方式，特征分解的計(jì)算復(fù)雜度由對W的O(N3)減小為對S的O(n3)。此外，對于歸一化譜聚類算法，還需要實(shí)現(xiàn)相似度矩陣的歸一化，即其中，為相似度矩陣的度矩陣，是對角元素的對角矩陣。由于有

式中，1表示單位列向量，ar，br∈Rn分別表示矩陣A和B的行和向量，bc∈Rm是矩陣B的列和向量。那么就可以將矩陣A和B歸一化為

3.3.2 MSSS采樣

Nystr?m 方法的關(guān)鍵在于采樣方法，本文使用考慮采樣點(diǎn)方差和相似度的逐步迭代MSSS 采樣方法（Bouneffouf 和Birol，2015）得到最佳采樣點(diǎn)。MSSS方法綜合了增量采樣IS（Incremental Sampling）和最小相似度采樣MSS（Minimum Similarity Sampling）。如圖3所示，該方法首先從完整數(shù)據(jù)集中隨機(jī)選擇兩個(gè)點(diǎn)到采樣點(diǎn)集中；然后從剩余點(diǎn)集中隨機(jī)取出一定比例的點(diǎn)，計(jì)算這些點(diǎn)與采樣點(diǎn)集中采樣點(diǎn)之間的相似度平方和，并將最小值點(diǎn)作為新的采樣點(diǎn)放入采樣點(diǎn)集中；重復(fù)該過程直到采樣點(diǎn)集的大小滿足要求。

圖3 MSSS采樣方法過程Fig.3 Flowchart of MSSS sampling method

根據(jù)Bouneffouf 和Birol（2015）的評估建議和模擬實(shí)驗(yàn)，將每次從剩余點(diǎn)集中隨機(jī)取出的子集T占剩余點(diǎn)集Y的比例設(shè)置為10%，期望的采樣點(diǎn)集合X的大小也確定為原始體素點(diǎn)集合的10%。

3.4 單木分割及參數(shù)提取

在計(jì)算相似度矩陣的特征值和特征向量之后，最終的聚類數(shù)目k可以通過特征值間隔啟發(fā)式（eigengap heuristic）來進(jìn)行選擇（Von Luxburg，2007）。令λ1，…，λn表示的特征值，即式（5）中ΛS的對角元素，目標(biāo)是選擇k值使得特征值λ1，…，λk的值都非常小，而λk+1值相對較大，即在第k個(gè)和第k+1 個(gè)特征值之間存在間隔，|λk+1-λk|比較大，該間隔表明數(shù)據(jù)集包含k個(gè)類。

接著，在特征空間中使用K-means 方法得到分割為k類的聚類結(jié)果。在特征空間中，K-means算法的直接分割對象是前k個(gè)特征向量構(gòu)成的矩陣的歸一化行元素，其得到的labels 值索引與體素點(diǎn)一一對應(yīng)，通過labels 值的索引就可以將結(jié)果映射回體素空間。同時(shí)，體素點(diǎn)與原始點(diǎn)云通過mean shift的聚類標(biāo)簽一一對應(yīng)，可以通過標(biāo)簽索引將分割結(jié)果從體素空間映射回原始點(diǎn)云空間，最終得到單木的聚類點(diǎn)云。

對于得到的單木點(diǎn)云聚類，單木參數(shù)直接從三維點(diǎn)云信息中獲取。在本研究中，單木位置和樹高由最高點(diǎn)的空間坐標(biāo)來確定，即(xtree，ytree，htree)=(xHighest，yHighest，zHighest)。

3.5 參數(shù)配置

在基于Nystr?m 的譜聚類算法中，有3 個(gè)參數(shù)對算法的單木分割效果有關(guān)鍵影響，即mean shift體素化方法中的核函數(shù)帶寬、高斯相似度函數(shù)的寬度參數(shù)σxy和σz。

對于mean shift 體素化中的關(guān)鍵參數(shù)核函數(shù)帶寬，由于體素化處理是要?jiǎng)?chuàng)建局部小單元（即體素），因此帶寬不宜過大，以生成大量細(xì)碎的聚類結(jié)果構(gòu)成體素空間。在本文算法中，核函數(shù)帶寬值由Python sklearn.cluster 模塊中的estimate_bandwidth 函數(shù)計(jì)算得到。其中的關(guān)鍵參數(shù)quantile的取值表示進(jìn)行近鄰搜索時(shí)的近鄰占樣本的比例，本文將點(diǎn)云密度視為平均水平下的局部近鄰值的近似表示，以根據(jù)點(diǎn)云數(shù)據(jù)特點(diǎn)控制體素化尺度，在保持體素空間合理性的同時(shí)壓縮后續(xù)計(jì)算數(shù)據(jù)量。由此，使用點(diǎn)云密度和點(diǎn)云數(shù)據(jù)點(diǎn)數(shù)的比值作為quantile 的參數(shù)值，將計(jì)算得到的帶寬值帶入sklearn.cluster 模塊中的MeanShift 函數(shù)以實(shí)現(xiàn)體素化過程。這里，除帶寬參數(shù)外，MeanShift 中的其他參數(shù)均采用默認(rèn)值。

對于式（1）中高斯相似度函數(shù)的寬度參數(shù)σxy和σz，即水平和垂直歐氏距離的比例因子，考慮到單木的形狀近似橢球，高度約為冠幅直徑的6倍左右，將σz設(shè)為σxy的6 倍，以使水平和垂直距離達(dá)到同一個(gè)尺度。根據(jù)實(shí)驗(yàn)經(jīng)驗(yàn)，取σxy= 3.16 m，由此得到σz取值為18.96 m。

3.6 精度評價(jià)

單木分割結(jié)果的準(zhǔn)確性從檢測率和單木參數(shù)精度兩個(gè)方面進(jìn)行評估。

在算法的單木檢測率方面，參考Eysn等（2015）和Wang 等（2016）的思路，對每株外業(yè)實(shí)測參考單木在一定空間范圍內(nèi)尋找與之匹配的分割單木。整個(gè)匹配搜索過程按照實(shí)測參考單木樹高遞減順序進(jìn)行，對每株實(shí)測參考單木，搜索其冠幅范圍內(nèi)滿足表2 中樹高差的分割單木作為最佳匹配的候選，這里的樹高差準(zhǔn)則是參考Eysn 等（2015）的設(shè)計(jì)。

當(dāng)多株分割單木滿足要求成為候選時(shí)，匹配過程從近到遠(yuǎn)進(jìn)行。如果更遠(yuǎn)的候選表現(xiàn)出更小的樹高差，且其距當(dāng)前實(shí)測參考的距離比較近候選的距離增量在2.5 m 之內(nèi)，那么更遠(yuǎn)的候選成為更佳匹配。重復(fù)此過程，直到檢查完所有實(shí)測參考單木。

圖4以一對匹配單木為例，展示了上述匹配方法的過程，其中的數(shù)值表示樹高。根據(jù)實(shí)測參考單木樹高為20 m 及表2 的準(zhǔn)則可得，樹高差準(zhǔn)則為ΔH＜4 m。最終，高21 m的分割單木成為高20 m的實(shí)測參考單木的最佳匹配，二者樹高差ΔH=1 m，距離為2.2 m。

表2 匹配方法的樹高準(zhǔn)則Table 2 Height criterion for individual tree matching

根據(jù)上述匹配方法，分割結(jié)果可被分為TP（True Positive）、FN（False Negative）和FP（False Positive）3 種，分別代表匹配、漏檢和過檢，且滿足TP + FN = Nreference和TP + FP = Ndetection，這里，Nreference和Ndetection分別表示實(shí)測參考單木和算法分割單木的株數(shù)。提取率（Rextr）、匹配率（Rmatch）、漏檢率（ROm）和過檢率（RCom）可以通過式（7）—（10）計(jì)算得出。

式中，Rextr為分割株數(shù)與實(shí)測株數(shù)的比值，故其取值可能大于100%，即為存在相對較多的過檢。此外，用RMSextr，RMSmatch，RMSOm和RMSCom分別表示多個(gè)樣地的提取率、匹配率、漏檢率和過檢率的均方根，以查看算法的整體效果。對于單木參數(shù)，計(jì)算所有匹配結(jié)果的樹高R2和RMSE 值，以評估參數(shù)的準(zhǔn)確性。

4 結(jié)果與分析

4.1 基于Nystr?m譜聚類算法的單木分割結(jié)果

圖5 展示了本文提出的基于Nystr?m 的譜聚類算法的分割結(jié)果示例，分別從頂視圖、斜視圖、側(cè)視圖展示了分割后的單木三維點(diǎn)云，相鄰單木的點(diǎn)云用不同顏色渲染顯示?？梢钥闯觯惴▽LS點(diǎn)云分割為合理的符合單木形態(tài)的結(jié)果。

圖5 基于Nystr?m譜聚類算法的分割示例Fig.5 Examples of segmentation results using Nystr?m-based spectral clustering algorithm form

4.2 分割精度評價(jià)

圖6 中展示了本文提出的基于Nystr?m 的譜聚類算法對研究區(qū)11 個(gè)樣地的分割結(jié)果的提取率和匹配率，并將匹配率按不同株數(shù)密度分不同灰度顯示?？梢钥闯?，本文算法對所有樣地的匹配率大致隨株數(shù)密度的減小而提高。從整體水平看（表3），算法提取了104%的參考單木，其中65%是正確匹配，這對于漏檢和過檢是一個(gè)較好的平衡。對于單木參數(shù)的提取精度，匹配結(jié)果被驗(yàn)證為具有較高的樹高精度，R2值和RMSE 值分別為0.88和1.57 m（圖7）。

圖6 基于Nystr?m-譜聚類算法的樣地單木檢測率Fig.6 Detection rates for individual tree segmentation using Nystr?m-based spectral clustering algorithm in field plots

圖7 匹配結(jié)果的樹高精度Fig.7 Height accuracy of the matched trees

為進(jìn)一步檢驗(yàn)算法的分割性能，分別考慮對不同株數(shù)密度和樹高層單木的識(shí)別情況。如表3所示，在株數(shù)密度方面，隨著株數(shù)密度的增加，RMSmatch從72%降低至61%。但是，在高株數(shù)密度樣地中發(fā)現(xiàn)了最佳提取率（RMSextr.=90%），比整體提取率低13%。低株數(shù)密度樣地的過檢率和漏檢率都是最少的，RMSCom和RMSOm分別為30%和28%。對于中等株數(shù)密度樣地而言，匹配率處于居中的水平，但是提取率和過檢率均較高，RMSextr.和RMSCom分別為106%和40%。在對不同樹高層的分割結(jié)果中（表5），本文的算法對20—25 m 和大于25 m 樹高層的單木匹配率分別為77%和78%，對于15—20 m 和小于15 m 單木的匹配率分別為45%和51%，該結(jié)果在表5 的對比中屬于較優(yōu)表現(xiàn)。

表3 基于Nystr?m-譜聚類算法的分割精度Table 3 Segmentation accuracy of Nystr?m-based spectral clustering algorithm

4.3 不同分割算法的對比

4.3.1 分割精度對比

為進(jìn)一步驗(yàn)證本文算法的單木分割性能，將本文提出的基于Nystr?m的譜聚類算法NSC（Nystr?mbased Spectral Clustering）與其他算法進(jìn)行對比，分別考慮K-means 算法和譜聚類算法SC（Spectral clustering）。K-means 算法作為目前ALS 點(diǎn)云數(shù)據(jù)單木分割問題中常用的聚類方法，與其進(jìn)行對比可以驗(yàn)證本文算法在基于聚類方法的單木分割算法中的表現(xiàn)；譜聚類算法作為本文算法的理論基礎(chǔ)，與其進(jìn)行對比可以查看本文算法對譜聚類進(jìn)行一系列改進(jìn)的有效性。對于算法參數(shù)設(shè)置，參考Gupta 等（2010）的工作，在K-means 算法中縮小點(diǎn)云高度以取得更好的結(jié)果，并通過多次實(shí)驗(yàn)選取最佳分割結(jié)果；對于SC 算法，由于其與本文算法的理論相似，故參考NSC 的最優(yōu)參數(shù)選擇方法對其進(jìn)行參數(shù)設(shè)置。

圖8 展示了不同算法在研究區(qū)的11 個(gè)樣地中的單木分割精度。這里，為進(jìn)行比較，對K-means算法和譜聚類算法均使用本文算法選擇的分割單木株數(shù)，故在同一樣地內(nèi)3 種算法的RMSextr.值一致?？梢钥闯?，譜聚類算法在所有樣地中均取得了最高的匹配率，本文算法的表現(xiàn)緊隨其后。同時(shí)，3 種算法對于9 號(hào)、10 號(hào)和11 號(hào)3 個(gè)低株數(shù)密度樣地的單木分割結(jié)果的匹配率較為接近，但K-means算法對中高株數(shù)密度樣地的單木識(shí)別精度較低。

圖8 不同分割算法的單木檢測率Fig.8 Detection rates for different segmentation algorithms

表4中對不同分割算法在不同株數(shù)密度樣地中的分割結(jié)果進(jìn)行了定量對比，其中，RMSmatch表示整體匹配率，而RMSmatch_H、RMSmatch_M、RMSmatch_L分別表示對高、中、低株數(shù)密度樣地匹配率的均方根結(jié)果。譜聚類算法在所有匹配結(jié)果中均取得了最優(yōu)的精度，整體匹配率為70%，且對不同株數(shù)密度的分割精度差距較小，對高、中、低株數(shù)密度樣地的匹配率分別達(dá)到68%、69%和75%，表明其相對穩(wěn)定的出色性能。K-means算法在整體匹配率上比譜聚類算法低7%，且分割精度隨株數(shù)密度的增加與譜聚類算法的差距逐漸增大，對低、中、高株數(shù)密度樣地的匹配率分別較譜聚類算法下降了5%、8%和14%。本文的算法基于譜聚類理論，通過體素化過程和Nystr?m 方法優(yōu)化計(jì)算性能（關(guān)于計(jì)算性能的討論可見4.3.2節(jié)），對譜聚類算法的分割精度產(chǎn)生了少許的損失，整體匹配率為65%，較譜聚類算法下降了5%，對高、中、低株數(shù)密度樣地的匹配率分別為61%、63%和72%，較譜聚類算法分別下降了7%、6%和3%，但仍優(yōu)于K-means 算法。對于單木參數(shù)，3 種方法得到的匹配樹高精度大致相同，R2值均在0.88 左右，RMSE值均為1.56 m左右。

表4 不同分割算法對不同株密度樣地的分割結(jié)果對比Table 4 Comparison of segmentation results by stem densities for different segmentation algorithms

表5比較了不同算法對不同樹高層的單木分割結(jié)果。同樣地，譜聚類算法在所有樹高層中均取得了最高的匹配率，對最高兩層大于25 m 和20—25 m 的匹配率分別達(dá)到83%和81%，對15—20 m和小于15 m 樹高層的匹配率分別為58%和51%。K-means算法對不同樹高層的匹配率均低于譜聚類算法，在最高兩層的匹配率均下降5%，而最低兩層則下降了10%左右。本文算法較譜聚類算法的匹配精度有少許下降，對最高兩層的匹配率下降4%—5%，最低兩層下降6%—7%，但仍優(yōu)于K-means算法。

表5 不同分割算法對不同樹高層的分割結(jié)果對比Table 5 Comparison of segmentation results by height layers for different segmentation algorithms

4.3.2 計(jì)算效率對比

除了分割精度對比，計(jì)算效率是算法性能的另一個(gè)重要衡量。圖9 比較了3 種算法在研究區(qū)的11 個(gè)樣地中的計(jì)算時(shí)間，算法均使用Python 3.6 編程實(shí)現(xiàn)，于64 位8 核3.6 GHz 主 頻的Windows 操作系統(tǒng)下運(yùn)行。為查看計(jì)算時(shí)間隨點(diǎn)云點(diǎn)數(shù)增加的變化趨勢，對樣地按照點(diǎn)數(shù)由小到大排序，樣地號(hào)及點(diǎn)數(shù)如x軸所示。

譜聚類算法雖然在分割精度上取得了最佳的結(jié)果（4.4.1 節(jié)），但是計(jì)算效率十分低下，由圖9可見，當(dāng)點(diǎn)數(shù)在30000 左右時(shí)，計(jì)算時(shí)間約為2000 s。本文的算法對11 個(gè)樣地的平均計(jì)算時(shí)間為17 s，是3 種算法中計(jì)算效率最高的，其計(jì)算效率最快約提高到譜聚類算法的96 倍（11 號(hào)樣地）和K-means算法的3倍（2號(hào)樣地）。其中，本文算法在5 號(hào)樣地的計(jì)算時(shí)間略高于K-means 算法，但整體水平仍優(yōu)于K-means 算法（平均計(jì)算時(shí)間30 s）。同時(shí)，由圖9 可以看出，隨著點(diǎn)數(shù)的增加，譜聚類算法和K-means 算法的計(jì)算時(shí)間都呈較為明顯的上升趨勢，而本文算法的計(jì)算時(shí)間對點(diǎn)數(shù)增加的反應(yīng)則較為平緩?？梢姡疚乃惴ㄔ谟?jì)算效率上比傳統(tǒng)譜聚類和K-means算法具有優(yōu)勢。

圖9 不同單木分割算法的計(jì)算時(shí)間Fig.9 Computing time for different segmentation algorithms

5 結(jié) 論

針對現(xiàn)有機(jī)載LiDAR 點(diǎn)云單木分割方法的不足，本文提出了基于Nystr?m 的譜聚類算法完成單木分割，并從獲取的單木聚類點(diǎn)云中提取關(guān)鍵單木因子。該算法基于譜聚類方法，同時(shí)引入了mean shift體素化和Nystr?m 方法，在保持譜聚類算法優(yōu)越表現(xiàn)的同時(shí)，減小了譜聚類算法的空間和時(shí)間復(fù)雜度。對落葉松人工林?jǐn)?shù)據(jù)單木分割的整體匹配率為65%，對低、中、高株數(shù)密度樣地的匹配率分別達(dá)到72%、63%、61%，對20—25 m和大于25 m 樹高層的單木匹配率為77%和78%。與現(xiàn)有方法相比，本文的算法在分割實(shí)驗(yàn)中整體表現(xiàn)良好，單木識(shí)別率與譜聚類算法差別較小，且優(yōu)于K-means 算法。在算法效率方面，本文的算法對11 個(gè)樣地的平均計(jì)算效率最高，最快約提高到譜聚類算法的96 倍和K-means 算法的3 倍。提取的單木樹高具有較高的精度，為機(jī)載LiDAR 點(diǎn)云數(shù)據(jù)的單木分割提供了一個(gè)可行的方案。

本文算法的改進(jìn)大大提高了譜聚類方法的計(jì)算效率，改進(jìn)后的算法更適合于大規(guī)模機(jī)載LiDAR 點(diǎn)云數(shù)據(jù)的單木分割實(shí)驗(yàn)。由于譜聚類算法本身對數(shù)據(jù)分布的適應(yīng)性很強(qiáng)，本文方法是對譜聚類算法的近似優(yōu)化，理論上可適用于多種林分情況（單一或混合樹種、不同株密度水平等）?？紤]到本文所選研究數(shù)據(jù)為落葉松人工林?jǐn)?shù)據(jù)，算法對其他林分情景的適用性還有待進(jìn)一步探索。

志 謝感謝中國林業(yè)科學(xué)研究院林業(yè)研究所、東北林業(yè)大學(xué)林學(xué)院、黑龍江省孟家崗林場在外業(yè)調(diào)查方面提供的幫助。