基于動態K閾值的蘋果葉片點云聚類與生長參數提取

劉 剛 張偉潔 郭彩玲

(1.中國農業大學現代精細農業系統集成研究教育部重點實驗室， 北京 100083；2.中國農業大學農業農村部農業信息獲取技術重點實驗室， 北京 100083； 3.唐山學院機電工程系， 唐山 063000)

0 引言

葉冠是果樹光合作用的主要場所，其空間形態結構直接影響冠內的光照分布、光合生產力及枝葉器官的生理狀態，進而影響果實品質和產量[1-2]。葉片作為葉冠的主要組成部分，獲取其空間分布規律及生長參數對獲取精確的葉冠結構具有重要的研究價值，也是進一步研究冠層光照分布計算的基礎[3]，可為果樹的自動化修形剪枝提供技術支持[4-6]。

葉片的識別方法多局限于二維空間中，雖然存在運算速度快等優點[7-8]，但難以進一步提取生長參數。地面三維激光掃描儀在農業領域的應用[9-12]使得快速獲取果樹冠層三維空間信息成為可能，進而可實現冠層葉片的自動分割及生長參數的提取，能在較大程度上減少人力，提高作業效率[13-15]。

國內外研究人員對點云分割方法進行了大量研究，主要分為6類[16]：基于區域邊界的分割方法適用于邊界明顯、分布規則的數據[17-18]；基于區域增長的分割方法適用于有明顯分布差異的點云[19]；基于圖的分割方法和基于模型匹配的分割方法對規則點云有很好的聚類效果[20]；基于機器學習的分割方法對訓練樣本的要求較高[21]；基于聚類的分割方法采用特征聚類，根據不同的點云分布特征，選擇適當的特征檢測方法[22-23]。果樹葉片空間分布相對密集，枝葉連續生長難以區分，因此利用特征聚類實現葉片分割較為合理。

葉片生長參數提取研究中，吳升等[24]通過計算葉片點云法向量與XOY坐標平面的夾角，實現點云在XOY坐標平面的正投影，進而將葉片幾何參數求解降維到二維空間上；鄭一力等[25]采用主成分分析與線性評判分析相結合的方法，實現了葉片的多特征降維計算。降維法計算葉片生長參數簡單快速，但忽略了葉片的三維空間特性。針對以上問題，本文提出一種葉片點云聚類分割算法，融合局部凹凸性算法(Locally convex connected patches，LCCP)并改進K-means算法，動態獲取K值，進而研究葉片的生長參數提取方法。

1 材料與方法

1.1 點云數據獲取方法

隨機選取中國農業大學蘋果樹采摘機器人試驗基地(北京市昌平區南口鎮辛力莊村)的1棵富士蘋果樹作為研究對象，樹高3.1 m，樹齡4年，采用美國Trimble 公司TX8型地面三維激光掃描儀(具體參數如表1所示)進行點云掃描。由于試驗果樹周圍還生長有其他果樹，遮擋較多，為了獲得高質量的點云，數據獲取試驗采用5站掃描，如圖1所示。其中的s1～s5為5個掃描點位，掃描儀的工作位置與被測果樹樹干的距離為3.5～5 m，c1～c6為用于配準定位的6個標靶球(標準直徑100 mm)，分布于試驗果樹周圍，保證在每站掃描能識別3個以上標靶球，以便多站點云的配準，掃描實地場景如圖2a所示。本試驗不考慮環境(如風速、溫濕度、光照度等)的影響。

表1 Trimble TX8型掃描儀產品參數Tab.1 Product parameters of Trimble TX8

圖1 掃描儀布置示意圖Fig.1 Scanner layout diagram

圖2 點云數據采集Fig.2 Point cloud acquisition

利用基于標靶球提取的算法[26]配準多站點云，點云數據獲取過程中由于外界環境干擾和設備自身的誤差，產生大量的噪點，采用前期研究成果[27]——基于距離均值計算的方法，剔除噪點。點云密度大且分布不均勻，采用均勻網格法完成抽稀處理，最終所得的蘋果樹冠層點云數據如圖2b所示，圖2a中白色小球為標靶球。

1.2 點云處理方法

本研究算法由C++語言編寫實現，點云處理庫為點云庫(Point cloud library，PCL)。獲取的點云數據以三維坐標(x，y，z)的形式存儲，數據量龐大，約有115萬點。為了簡化果樹葉片點云的搜索和存儲，創建新的葉片點云PCD存儲格式，每條記錄存儲一個葉片節點，包括的屬性有葉片中心點Li(xi,yi,zi)、葉片法向量Ni、葉長Lleni和葉寬Lwidi。

本文研究方法流程主要分為4部分(圖3)：數據獲取、點云預處理、葉片點云聚類分割和葉片生長參數提取。數據獲取試驗需在晴朗無風條件下進行；點云預處理主要包括點云配準、去噪和抽稀；采用基于動態K閾值的算法，融合LCCP算法和改進的K-means算法，隨機選取冠層的一個枝條，作為研究對象，進行葉片點云聚類分割；葉片生長參數提取方法采用PCA方法和邊界提取方法，分別實現葉片生長角度和葉長葉寬參數值的獲取。

圖3 葉片點云分割及生長參數提取流程圖Fig.3 Flow chart of apple tree point cloud segmentation and growth parameters extraction

1.2.1葉片點云的聚類方法

葉片點云聚類分割，采用基于動態K閾值的算法，融合LCCP算法，并采用改進的K-means算法，實現K閾值的動態獲取。該算法主要包含以下過程：超體聚類、LCCP聚類[28]、K閾值獲取和K-means分割。

超體聚類是一種過分割方法，過分割涉及的參數有體素距離Rv、晶核距離Rs。八叉樹初始化點云，過分割生成獨立的超體素塊pi=(xi,ni,Ei)，其中pi為超體素索引、xi為體素中心、ni為體素法向量、Ei為鄰接體素邊集合，計算鄰接體素相似度完成超體聚類

(1)

式中D——體素合并的概率

wc、ws、wn——顏色、空間、法向量的信息權重

Dc、Ds、Dn——顏色、空間、法向量容差

以超體聚類所得體素集合為研究對象，建立體素索引向量，采用CC(Extended convexity criterion)和SC(Sanity criterion)方法判定鄰接體素凹凸性關系完成LCCP聚類。

(2)

(3)

β=∠(n1,n2)

(4)

式中 CCb(pi,pj)——鄰接體素塊凹凸性

n1、n2——鄰接體素塊法向量

βT——法向量角度偏差閾值

i、j——體素索引號

當鄰接體素面不連接時，CC的判定準則不適用，此時需要使用SC準則對鄰接體素關系進行判定，SC判定依據為

(5)

其中θ(pi,pj)=min(∠(d,s),∠(d,-s))=

min(∠(d,s),180°-∠(d,s))

(6)

(7)

s=n1×n2

(8)

式中 SCb(pi,pj)——鄰接體素凹凸性

conv(pi,pj)=CCb(pi,pj)∧SCb(pi,pj)

(9)

式中 conv(pi,pj)——pi、pj的凹凸性

動態K閾值的獲取主要分為兩個步驟：一方面計算K=4時，K-means聚類后各點集三維質心間歐氏距離；另一方面計算葉片中心點到擬合直線的距離均值，作為距離閾值。動態獲取K閾值的具體方法如下：

(1)根據LCCP所得類簇，設定K=4，進行K-means聚類分割，并計算4個點群的中心點Ckm(xm,ym,zm)。

(2)計算4個中心點之間的6組歐氏距離(其中m、n為類簇中心點編號)

(10)

(3)手動截取枝條葉片點云，計算所有葉片點云中心點Cs(xs,ys,zs)，根據所得點云中心點擬合直線L，計算葉片中心點到直線L的距離均值，作為距離閾值Dt，此閾值為統計結果，不需要重復計算。

(4)將步驟(2)所得的中心點歐氏距離，與距離閾值Dt進行比較，大于、小于距離閾值的個數分別為a、b，根據a與b的大小關系，確定此點群需要分割的點群個數K，其中K值的選定準則為

K=4(a=6)a+12(b

(11)

式中 ?·」——向下取整

計算每個類簇的K值，再次進行K-means聚類，得到單葉片和枝干點云。

1.2.2葉片的生長參數提取方法

將葉片點云三維空間信息作為輸入參數，利用PCA法構建協方差矩陣，計算葉片點云的法向量，即計算由葉片點云三維坐標組成的協方差矩陣的特征值和特征向量。構建的協方差矩陣為

(12)

Avj=λjvj(j=0,1,2)

(13)

式中n0——葉片點云總數

Pi——葉片點坐標

λj——第j個特征值

vj——λj對應的特征向量

計算所得的最大特征值max_eigenvalue對應的特征向量max_eigenvector()作為葉片法平面方向向量，可根據計算所得結果，進一步計算葉方位角、葉傾角等參數。

圖4為葉片生長參數求解示意圖，根據前期研究成果，提取葉片邊界點云，計算邊界點中點與葉片中心點的最大距離，作為葉片的葉長。如圖4b所示，N為計算所得葉片生長方向向量，O為三維質心，l1和l2為葉片長軸端點，與l1、O、l2所在平面近似垂直的最大歐氏距離為葉寬，w1和w2為葉片寬軸端點。

圖4 葉片點云邊界提取及參數求解示意圖Fig.4 Sketch of leaf points boundary extraction and parameter solution

圖5 傳統方法聚類效果圖Fig.5 Sketches of segmentation by traditional methods

2 結果與分析

2.1 葉片點云聚類效果與精度

葉片的聚類分割方法主要考慮了葉片點云的分布特征。葉片與枝干連接且分布密集，導致葉片與葉片相互區分困難、葉片與枝干相互分割困難，若采用傳統的聚類算法則會有較多不足，圖5采用3種傳統的點云聚類方法，對同一單枝點云數據進行處理，圖中可見，采用DBSCAN算法無法將連續的枝、葉點云分割開；另外在點云獲取過程中由于遮擋等原因，會導致局部點云缺失，若僅利用K-means算法或Region-growing算法，會出現大量的單葉片被分割成兩個到多個聚類的情況，導致所得單葉片點云不完整。

如上分析，需要找到一種既能成功分割單葉片點云，也能保證其完整性的方法。如圖6所示，采用超體聚類和LCCP算法可以顧及點云局部凹凸性，將整體枝條點云分割為多簇，即將整體長枝條分割為較小的具有多個完整葉片的點云類簇，如圖6b中的黃色立方體框為點云類簇的包圍盒，此步驟能較好地保證葉片的完整性，避免一個葉片被分割到不同的類簇中。

圖6 超體聚類與LCCP聚類效果圖Fig.6 Sketches of super-voxel and LCCP clustering

圖7 動態K閾值分割結果Fig.7 Sketches of dynamic-K-threshold segmentation

超體聚類中涉及的參數設置主要有Rv=10.0,Rs=0.3，由于點云無顏色信息，故ws=1.0,wn=10.0,Ds=1.0,Dn=1.0,wc=0,Dc=0。由圖6可看出，葉片分割最大的難點是將葉片點云與連接的枝干點云進行分割。LCCP聚類結果能夠保證葉片點云的完整性，但不能直接將葉片、枝干分割開來，而K-means算法將距離作為相似性度量的標準，葉片、枝干點云間具有明顯的區域分布性，根據不同的LCCP聚類結果，設定不同的K值完成K-means聚類，即將葉片點云進行有效分割。因此，設定K-means算法中的K=4，計算所得4組聚類中兩兩的中心點距離，與距離閾值進行比較，根據動態K閾值判定準則可以動態獲取K值，再次利用K-means方法聚類，得到枝葉分離的點云聚類結果。主要包含葉-葉(圖7b)、枝-葉(圖7c)、葉-枝-葉(圖7d)以及葉-枝-葉-葉(圖7e)4種情況。圖8a為枝條葉片最終分割結果，圖8b為最終蘋果樹冠層葉片分割結果。

圖7為不同K閾值所得分割結果，圖7a為動態K=1的分割結果，只包含1個葉片的點云數據；圖7b、7c為K=2的分割結果，圖7d為K=3的分割結果，圖7e為K=4的分割結果，圖7所包含的結果只代表不同的動態K值的情況，不包含所有的葉片分割結果。識別葉片和枝干點云，需要借助點云可視化的方法。為了驗證本文提出的葉片分割方法的有效性，將本文所得冠層葉片分割結果與多種方法所得結果進行比較，如表2所示。

分析表2數據可得，DBSCAN方法所得聚類葉片數量較少，是由于單個聚類點集中往往會包含1個以上葉片。相反地，K-means方法和Region-growing方法所得葉片相對較多，是由于點云中若存在少量空洞或局部缺失則會導致數據不連續，算法執行時會將此區域聚類為獨立的點集。本文提出的葉片聚類方法，葉片識別的正確率達到90%以上，另外，動態K閾值的獲取方法也減少了人工干預，較大提高了運算效率。

2.2 生長參數提取效果與精度

葉片點云作為輸入參數，利用PCA法構建協方差矩陣并計算葉片點云法向量，作為葉片生長方向向量；利用邊界提取方法獲取葉片邊界點云，計算邊界點與葉片中心點的歐氏距離，從而獲取葉長、葉寬的參數值。為進一步驗證本文提出的生長參數提取方法的可行性，隨機選取10組葉片，每組中包含20個葉片，定義平均葉長、葉寬作為評價指標，分別記錄利用實測方法、投影法、本文提出的方法所得結果，如表3所示。可以看出，利用本文提出的生長參數提取方法，葉片生長角度誤差小于5°的葉片數量占葉片總數量的90%以上，葉長、葉寬參數值的平均準確率均達到91%以上，相較于降維法所得結果，本文提出的方法能夠進一步提高提取精度。

圖8 葉片點云分割效果圖Fig.8 Sketches of leaf points segmentation

聚類方法葉片分割數量/個正確率/%本文方法110290.77K-means139485.18DBSCAN100782.95Region-growing143881.55實測值1214

3 結論

(1)根據蘋果樹葉片點云分布特征，提出了基于動態K閾值的葉片聚類方法，實現了冠層枝葉點云的分離，進而完成了葉片生長參數的提取。基于動態K閾值的算法，克服了傳統算法中人工干預過多的缺點，可以有效提高生長參數的提取精度。

表3 葉片生長參數提取的精度對比Tab.3 Comparison of accuracy of leaf growth parameters extraction

(2)將本文方法與傳統聚類方法的試驗效果進行對比，結果顯示，本文方法能夠正確識別出的葉片數占冠層總葉片數的90%以上；本文提出的葉片生長參數提取方法所得結果，生長角度誤差小于5°的在90%以上，葉長、葉寬平均準確率均在91%以上，與投影法提取準確率比較表明，生長參數的提取精度得到提高。

(3)本文提取的葉片生長參數種類較少，未實現枝葉點云的自動化分類識別，在后續的研究中，需要借助可視化的方法對此進行完善和實現。