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

1.中國(guó)林業(yè)科學(xué)研究院 資源信息研究所,北京 100091;2.國(guó)家林業(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)測(cè)中的應(yīng)用日益廣泛（Hyypp? 等，2008）。機(jī)載激光雷達(dá)ALS（Airborne Laser Scanning）可以精確刻畫(huà)冠層垂直結(jié)構(gòu)，提取林分和單木尺度的森林參數(shù)（Sa?kov 等，2017；李增元 等，2016）。ALS點(diǎn)云數(shù)據(jù)的單木分割對(duì)林業(yè)科學(xué)研究和生產(chǎn)實(shí)踐具有重要意義，基于準(zhǔn)確的分割結(jié)果，可以從中獲取描述單木空間結(jié)構(gòu)特征和生物化學(xué)組分屬性的單木因子，為評(píng)價(jià)森林生長(zhǎng)和生態(tài)功能、評(píng)估森林破壞程度和監(jiān)測(cè)森林再生情況提供數(shù)據(jù)支持（Chen等，2007）。

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

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

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

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

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

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

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

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

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

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

本研究采用中國(guó)林業(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ī)，平均絕對(duì)飛行高度為1000 m，相對(duì)飛行高度約650 m。ALS 點(diǎn)云采集于LiCHy 系統(tǒng)，點(diǎn)云數(shù)據(jù)密度約為20 pts/m2。對(duì)ALS點(diǎn)云的處理包括基于體素的孤立點(diǎn)去除、基于不規(guī)則加密三角網(wǎng)TIN （Triangulated Irregular Network）的地面點(diǎn)分類(lèi)、高度歸一化處理和地面點(diǎn)去除，此點(diǎn)云數(shù)據(jù)作為原始點(diǎn)云數(shù)據(jù)。

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

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

3 原理與方法

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

圖2 基于Nystr?m-譜聚類(lèi)算法的單木分割流程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 是一種快速有效的聚類(lèi)算法，本文使用mean shift 方法來(lái)完成體素化過(guò)程，并將mean shift的帶寬參數(shù)設(shè)置得較小，以實(shí)現(xiàn)過(guò)分割來(lái)達(dá)到體素化目的（帶寬參數(shù)的詳細(xì)設(shè)置見(jiàn)3.5 節(jié)中關(guān)于參數(shù)配置部分）。與使用給定分辨率的立方體來(lái)構(gòu)建規(guī)則體素點(diǎn)云的方法不同，在使用mean shift 方法構(gòu)造的體素空間中，每個(gè)體素點(diǎn)由mean shift 的聚類(lèi)結(jié)果表示，體素位置由類(lèi)內(nèi)中心點(diǎn)坐標(biāo)確定，體素權(quán)重等于類(lèi)內(nèi)的點(diǎn)數(shù)。后續(xù)的分割過(guò)程將基于體素空間完成，可以通過(guò)mean shift 聚類(lèi)結(jié)果的類(lèi)別標(biāo)簽將原始點(diǎn)云中的每個(gè)點(diǎn)對(duì)應(yīng)到體素點(diǎn)云中，從而將體素空間中得到的分類(lèi)結(jié)果映射回原始點(diǎn)云，最終得到單木的分割結(jié)果。同時(shí)，可以通過(guò)體素權(quán)重將原始點(diǎn)云的空間分布信息表征在體素空間中，即，體素權(quán)重越大表明該局部區(qū)域的原始點(diǎn)密度越高。在實(shí)驗(yàn)過(guò)程中，使用mean shift聚類(lèi)算法得到的體素點(diǎn)云數(shù)據(jù)量較原始點(diǎn)云約減少了10 倍，大大減輕了算法在后續(xù)過(guò)程中的計(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)沒(méi)有規(guī)則的空間位置索引，故無(wú)法像傳統(tǒng)立方體體素化方法那樣快速訪(fǎng)問(wèn)k-近鄰，本文采用KD-Tree數(shù)據(jù)結(jié)構(gòu)來(lái)彌補(bǔ)這一不足。

3.3 Nystr?m方法

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

3.3.1 Nystr?m理論

對(duì)一有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)間的相似度矩陣，且有對(duì)角化形式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，且其對(duì)角化為S=USΛS，其中，US∈Rn×n且列元素為為S的特征向量，ΛS是對(duì)角元素為S對(duì)應(yīng)特征值的對(duì)角矩陣。那么=VΛSVT即為W的近似特征分解，且有

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

式中，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ù)該過(guò)程直到采樣點(diǎn)集的大小滿(mǎn)足要求。

圖3 MSSS采樣方法過(guò)程Fig.3 Flowchart of MSSS sampling method

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

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

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

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

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

3.5 參數(shù)配置

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

對(duì)于mean shift 體素化中的關(guān)鍵參數(shù)核函數(shù)帶寬，由于體素化處理是要?jiǎng)?chuàng)建局部小單元（即體素），因此帶寬不宜過(guò)大，以生成大量細(xì)碎的聚類(lèi)結(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)體素化過(guò)程。這里，除帶寬參數(shù)外，MeanShift 中的其他參數(shù)均采用默認(rèn)值。

對(duì)于式（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 精度評(píng)價(jià)

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

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

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

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

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

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

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

4 結(jié)果與分析

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

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

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

4.2 分割精度評(píng)價(jià)

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

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

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

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

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

4.3 不同分割算法的對(duì)比

4.3.1 分割精度對(duì)比

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

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

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

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

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

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

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

4.3.2 計(jì)算效率對(duì)比

除了分割精度對(duì)比，計(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ù)增加的變化趨勢(shì)，對(duì)樣地按照點(diǎn)數(shù)由小到大排序，樣地號(hào)及點(diǎn)數(shù)如x軸所示。

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

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

5 結(jié) 論

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

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

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