基于動(dòng)態(tài)K閾值的蘋(píng)果葉片點(diǎn)云聚類與生長(zhǎng)參數(shù)提取

劉 剛 張偉潔 郭彩玲

(1.中國(guó)農(nóng)業(yè)大學(xué)現(xiàn)代精細(xì)農(nóng)業(yè)系統(tǒng)集成研究教育部重點(diǎn)實(shí)驗(yàn)室， 北京 100083；2.中國(guó)農(nóng)業(yè)大學(xué)農(nóng)業(yè)農(nóng)村部農(nóng)業(yè)信息獲取技術(shù)重點(diǎn)實(shí)驗(yàn)室， 北京 100083； 3.唐山學(xué)院機(jī)電工程系， 唐山 063000)

0 引言

葉冠是果樹(shù)光合作用的主要場(chǎng)所，其空間形態(tài)結(jié)構(gòu)直接影響冠內(nèi)的光照分布、光合生產(chǎn)力及枝葉器官的生理狀態(tài)，進(jìn)而影響果實(shí)品質(zhì)和產(chǎn)量[1-2]。葉片作為葉冠的主要組成部分，獲取其空間分布規(guī)律及生長(zhǎng)參數(shù)對(duì)獲取精確的葉冠結(jié)構(gòu)具有重要的研究?jī)r(jià)值，也是進(jìn)一步研究冠層光照分布計(jì)算的基礎(chǔ)[3]，可為果樹(shù)的自動(dòng)化修形剪枝提供技術(shù)支持[4-6]。

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

國(guó)內(nèi)外研究人員對(duì)點(diǎn)云分割方法進(jìn)行了大量研究，主要分為6類[16]：基于區(qū)域邊界的分割方法適用于邊界明顯、分布規(guī)則的數(shù)據(jù)[17-18]；基于區(qū)域增長(zhǎng)的分割方法適用于有明顯分布差異的點(diǎn)云[19]；基于圖的分割方法和基于模型匹配的分割方法對(duì)規(guī)則點(diǎn)云有很好的聚類效果[20]；基于機(jī)器學(xué)習(xí)的分割方法對(duì)訓(xùn)練樣本的要求較高[21]；基于聚類的分割方法采用特征聚類，根據(jù)不同的點(diǎn)云分布特征，選擇適當(dāng)?shù)奶卣鳈z測(cè)方法[22-23]。果樹(shù)葉片空間分布相對(duì)密集，枝葉連續(xù)生長(zhǎng)難以區(qū)分，因此利用特征聚類實(shí)現(xiàn)葉片分割較為合理。

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

1 材料與方法

1.1 點(diǎn)云數(shù)據(jù)獲取方法

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

表1 Trimble TX8型掃描儀產(chǎn)品參數(shù)Tab.1 Product parameters of Trimble TX8

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

圖2 點(diǎn)云數(shù)據(jù)采集Fig.2 Point cloud acquisition

利用基于標(biāo)靶球提取的算法[26]配準(zhǔn)多站點(diǎn)云，點(diǎn)云數(shù)據(jù)獲取過(guò)程中由于外界環(huán)境干擾和設(shè)備自身的誤差，產(chǎn)生大量的噪點(diǎn)，采用前期研究成果[27]——基于距離均值計(jì)算的方法，剔除噪點(diǎn)。點(diǎn)云密度大且分布不均勻，采用均勻網(wǎng)格法完成抽稀處理，最終所得的蘋(píng)果樹(shù)冠層點(diǎn)云數(shù)據(jù)如圖2b所示，圖2a中白色小球?yàn)闃?biāo)靶球。

1.2 點(diǎn)云處理方法

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

本文研究方法流程主要分為4部分(圖3)：數(shù)據(jù)獲取、點(diǎn)云預(yù)處理、葉片點(diǎn)云聚類分割和葉片生長(zhǎng)參數(shù)提取。數(shù)據(jù)獲取試驗(yàn)需在晴朗無(wú)風(fēng)條件下進(jìn)行；點(diǎn)云預(yù)處理主要包括點(diǎn)云配準(zhǔn)、去噪和抽??；采用基于動(dòng)態(tài)K閾值的算法，融合LCCP算法和改進(jìn)的K-means算法，隨機(jī)選取冠層的一個(gè)枝條，作為研究對(duì)象，進(jìn)行葉片點(diǎn)云聚類分割；葉片生長(zhǎng)參數(shù)提取方法采用PCA方法和邊界提取方法，分別實(shí)現(xiàn)葉片生長(zhǎng)角度和葉長(zhǎng)葉寬參數(shù)值的獲取。

圖3 葉片點(diǎn)云分割及生長(zhǎng)參數(shù)提取流程圖Fig.3 Flow chart of apple tree point cloud segmentation and growth parameters extraction

1.2.1葉片點(diǎn)云的聚類方法

葉片點(diǎn)云聚類分割，采用基于動(dòng)態(tài)K閾值的算法，融合LCCP算法，并采用改進(jìn)的K-means算法，實(shí)現(xiàn)K閾值的動(dòng)態(tài)獲取。該算法主要包含以下過(guò)程：超體聚類、LCCP聚類[28]、K閾值獲取和K-means分割。

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

(1)

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

wc、ws、wn——顏色、空間、法向量的信息權(quán)重

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

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

(2)

(3)

β=∠(n1,n2)

(4)

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

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

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

i、j——體素索引號(hào)

當(dāng)鄰接體素面不連接時(shí)，CC的判定準(zhǔn)則不適用，此時(shí)需要使用SC準(zhǔn)則對(duì)鄰接體素關(guān)系進(jìn)行判定，SC判定依據(jù)為

(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的凹凸性

動(dòng)態(tài)K閾值的獲取主要分為兩個(gè)步驟：一方面計(jì)算K=4時(shí)，K-means聚類后各點(diǎn)集三維質(zhì)心間歐氏距離；另一方面計(jì)算葉片中心點(diǎn)到擬合直線的距離均值，作為距離閾值。動(dòng)態(tài)獲取K閾值的具體方法如下：

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

(2)計(jì)算4個(gè)中心點(diǎn)之間的6組歐氏距離(其中m、n為類簇中心點(diǎn)編號(hào))

(10)

(3)手動(dòng)截取枝條葉片點(diǎn)云，計(jì)算所有葉片點(diǎn)云中心點(diǎn)Cs(xs,ys,zs)，根據(jù)所得點(diǎn)云中心點(diǎn)擬合直線L，計(jì)算葉片中心點(diǎn)到直線L的距離均值，作為距離閾值Dt，此閾值為統(tǒng)計(jì)結(jié)果，不需要重復(fù)計(jì)算。

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

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

(11)

式中 ?·」——向下取整

計(jì)算每個(gè)類簇的K值，再次進(jìn)行K-means聚類，得到單葉片和枝干點(diǎn)云。

1.2.2葉片的生長(zhǎng)參數(shù)提取方法

將葉片點(diǎn)云三維空間信息作為輸入?yún)?shù)，利用PCA法構(gòu)建協(xié)方差矩陣，計(jì)算葉片點(diǎn)云的法向量，即計(jì)算由葉片點(diǎn)云三維坐標(biāo)組成的協(xié)方差矩陣的特征值和特征向量。構(gòu)建的協(xié)方差矩陣為

(12)

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

(13)

式中n0——葉片點(diǎn)云總數(shù)

Pi——葉片點(diǎn)坐標(biāo)

λj——第j個(gè)特征值

vj——λj對(duì)應(yīng)的特征向量

計(jì)算所得的最大特征值max_eigenvalue對(duì)應(yīng)的特征向量max_eigenvector()作為葉片法平面方向向量，可根據(jù)計(jì)算所得結(jié)果，進(jìn)一步計(jì)算葉方位角、葉傾角等參數(shù)。

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

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

圖5 傳統(tǒng)方法聚類效果圖Fig.5 Sketches of segmentation by traditional methods

2 結(jié)果與分析

2.1 葉片點(diǎn)云聚類效果與精度

葉片的聚類分割方法主要考慮了葉片點(diǎn)云的分布特征。葉片與枝干連接且分布密集，導(dǎo)致葉片與葉片相互區(qū)分困難、葉片與枝干相互分割困難，若采用傳統(tǒng)的聚類算法則會(huì)有較多不足，圖5采用3種傳統(tǒng)的點(diǎn)云聚類方法，對(duì)同一單枝點(diǎn)云數(shù)據(jù)進(jìn)行處理，圖中可見(jiàn)，采用DBSCAN算法無(wú)法將連續(xù)的枝、葉點(diǎn)云分割開(kāi)；另外在點(diǎn)云獲取過(guò)程中由于遮擋等原因，會(huì)導(dǎo)致局部點(diǎn)云缺失，若僅利用K-means算法或Region-growing算法，會(huì)出現(xiàn)大量的單葉片被分割成兩個(gè)到多個(gè)聚類的情況，導(dǎo)致所得單葉片點(diǎn)云不完整。

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

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

圖7 動(dòng)態(tài)K閾值分割結(jié)果Fig.7 Sketches of dynamic-K-threshold segmentation

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

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

分析表2數(shù)據(jù)可得，DBSCAN方法所得聚類葉片數(shù)量較少，是由于單個(gè)聚類點(diǎn)集中往往會(huì)包含1個(gè)以上葉片。相反地，K-means方法和Region-growing方法所得葉片相對(duì)較多，是由于點(diǎn)云中若存在少量空洞或局部缺失則會(huì)導(dǎo)致數(shù)據(jù)不連續(xù)，算法執(zhí)行時(shí)會(huì)將此區(qū)域聚類為獨(dú)立的點(diǎn)集。本文提出的葉片聚類方法，葉片識(shí)別的正確率達(dá)到90%以上，另外，動(dòng)態(tài)K閾值的獲取方法也減少了人工干預(yù)，較大提高了運(yùn)算效率。

2.2 生長(zhǎng)參數(shù)提取效果與精度

葉片點(diǎn)云作為輸入?yún)?shù)，利用PCA法構(gòu)建協(xié)方差矩陣并計(jì)算葉片點(diǎn)云法向量，作為葉片生長(zhǎng)方向向量；利用邊界提取方法獲取葉片邊界點(diǎn)云，計(jì)算邊界點(diǎn)與葉片中心點(diǎn)的歐氏距離，從而獲取葉長(zhǎng)、葉寬的參數(shù)值。為進(jìn)一步驗(yàn)證本文提出的生長(zhǎng)參數(shù)提取方法的可行性，隨機(jī)選取10組葉片，每組中包含20個(gè)葉片，定義平均葉長(zhǎng)、葉寬作為評(píng)價(jià)指標(biāo)，分別記錄利用實(shí)測(cè)方法、投影法、本文提出的方法所得結(jié)果，如表3所示。可以看出，利用本文提出的生長(zhǎng)參數(shù)提取方法，葉片生長(zhǎng)角度誤差小于5°的葉片數(shù)量占葉片總數(shù)量的90%以上，葉長(zhǎng)、葉寬參數(shù)值的平均準(zhǔn)確率均達(dá)到91%以上，相較于降維法所得結(jié)果，本文提出的方法能夠進(jìn)一步提高提取精度。

圖8 葉片點(diǎn)云分割效果圖Fig.8 Sketches of leaf points segmentation

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

3 結(jié)論

(1)根據(jù)蘋(píng)果樹(shù)葉片點(diǎn)云分布特征，提出了基于動(dòng)態(tài)K閾值的葉片聚類方法，實(shí)現(xiàn)了冠層枝葉點(diǎn)云的分離，進(jìn)而完成了葉片生長(zhǎng)參數(shù)的提取?；趧?dòng)態(tài)K閾值的算法，克服了傳統(tǒng)算法中人工干預(yù)過(guò)多的缺點(diǎn)，可以有效提高生長(zhǎng)參數(shù)的提取精度。

表3 葉片生長(zhǎng)參數(shù)提取的精度對(duì)比Tab.3 Comparison of accuracy of leaf growth parameters extraction

(2)將本文方法與傳統(tǒng)聚類方法的試驗(yàn)效果進(jìn)行對(duì)比，結(jié)果顯示，本文方法能夠正確識(shí)別出的葉片數(shù)占冠層總?cè)~片數(shù)的90%以上；本文提出的葉片生長(zhǎng)參數(shù)提取方法所得結(jié)果，生長(zhǎng)角度誤差小于5°的在90%以上，葉長(zhǎng)、葉寬平均準(zhǔn)確率均在91%以上，與投影法提取準(zhǔn)確率比較表明，生長(zhǎng)參數(shù)的提取精度得到提高。

(3)本文提取的葉片生長(zhǎng)參數(shù)種類較少，未實(shí)現(xiàn)枝葉點(diǎn)云的自動(dòng)化分類識(shí)別，在后續(xù)的研究中，需要借助可視化的方法對(duì)此進(jìn)行完善和實(shí)現(xiàn)。