基于三維點云的大豆植株器官分割及表型分析

肖奕同， 劉帥， 侯晨連， 劉琦， 李富忠， 張吳平

（山西農(nóng)業(yè)大學軟件學院，山西 太谷 030801）

植物表型是指在基因和環(huán)境相互作用下產(chǎn)生的能夠直觀反映植物生長狀態(tài)的特征及性狀[1]。植物表型參數(shù)的分析與育種息息相關[2]，然而表型參數(shù)的測量結果受人員、環(huán)境和設備等因素影響較大，極大限制了大規(guī)模遺傳育種篩選的效率[3]。因此，開發(fā)高效、準確的表型參數(shù)測量方法對提高植物表型測量效率、促進植物育種至關重要?；诙S圖像的植物表型測量方法因設備要求低、易于使用等優(yōu)勢逐漸應用于農(nóng)業(yè)工程領域[4-6]。但是，二維圖像在表示成簇植株的冠層時有較大缺陷[7]，冠層結構中被遮擋的部位在圖像中完全缺失。三維點云不僅具有色彩信息，同時彌補了二維圖像深度信息不足的缺陷，在很大程度上避免了因葉片粘連、遮擋等問題而難以準確測量表型的問題。史蒲娟等[8]基于從不同角度拍攝的圖像序列，對油菜植株進行三維重建，提取了株高、葉柄長及葉長等表型參數(shù)。李楊先等[9]采集不同視角的二維圖像并結合運動恢復結構算法實現(xiàn)簸箕柳植株稀疏點云的重建，提取了株高、基徑、葉片面積、分枝數(shù)和分枝角等表型參數(shù)。近年來，隨著對植物表型分析需求的日益增多及三維點云表型研究的深入，對表型參數(shù)測量的精準度也有了更高的要求[10]，精準的表型分析需要依靠植株器官水平或具體部位的信息[11-13]。已有研究均是基于完整的植株點云數(shù)據(jù)進行表型分析，沒有實現(xiàn)植株器官尺度的分割。朱超等[14]基于點云骨架和最優(yōu)傳輸距離實現(xiàn)了玉米點云莖葉分割，分割算法平均總體準確率達到0.92。陽旭等[15]基于棉花植株的多時序點云數(shù)據(jù)，采用基于隨機采樣一致的區(qū)域生長算法完成冠層提取及分割。然而，相較于玉米及棉花植株，大豆植株葉片成簇且葉片間相互粘連現(xiàn)象嚴重，通過骨架提取[14]及柱體擬合[15]的方法均難以快速實現(xiàn)大豆植株莖葉分割。對植株點云進行高效、精確的器官分割既是表型分析的前提也是植物表型研究的重大難題[16]，為了解決多分枝作物葉片相互粘連帶來的器官分割及高通量表型測量困難問題，本研究基于大豆植株的整體稠密三維點云，利用點云法向量、曲率差異，結合基于歐式聚類的區(qū)域生長算法實現(xiàn)大豆植株器官分割，并提取葉面積、葉寬、葉長、葉傾角、莖粗等表型參數(shù)，為多分枝作物莖葉分割、單葉分割以及表型參數(shù)的測量提供技術支持。

1 材料與方法

1.1 數(shù)據(jù)獲取

本研究選取10 株分枝期盆栽大豆為研究對象。將植株置于光照均勻且無風的場景，利用RGB 相機（GoPro，HERO8 Black）以半球形狀環(huán)繞大豆植株進行多視角拍攝（圖1）。根據(jù)株型差異，采集30~50 張多視角圖像組成圖像序列。圖像序列中相鄰兩幅圖像間相互重疊區(qū)域達到70%以上。最后，分別測量各植株花盆口外徑。

圖1 相機位置參照Fig. 1 Camera position reference

1.2 點云重建與濾波

基于SFM（structure from motion）+MVS（multi view stereo）方法，利用Agisoft Photoscan 軟件重構大豆植株稀疏點云（圖2A）及稠密點云（圖2B）。基于人工測量所得的真實花盆口外徑，利用軟件中的標尺工具實現(xiàn)植株點云尺度變換。

圖2 重建的大豆植株三維點云Fig. 2 Reconstructed 3D point clouds of soybean plant

在獲取大豆植株點云數(shù)據(jù)時，由于設備精度、環(huán)境因素等影響，點云數(shù)據(jù)中不可避免地出現(xiàn)一些離植株點云主體較遠的離群點和噪聲點[17]。除了這些隨機誤差產(chǎn)生的噪聲點之外，受作物本身生長限制，如葉片交叉、葉片共面、葉片重疊等因素的影響，點云數(shù)據(jù)中往往存在一些貼近植株邊緣的冗余點。為了更好地進行器官分割、表型分析等后續(xù)處理，基于PCL（point cloud library）[18]對植株點云進行濾波處理，包括以下幾個步驟：①直通濾波過濾掉地面及花盆；②顏色濾波剔除與植株本身顏色差異較大的點，例如土壤及葉片邊緣的雜色冗余點；③統(tǒng)計濾波剔除離群點；④降采樣降低點云密度，減少計算時間。經(jīng)過濾波處理后的大豆植株點云如圖3所示。

圖3 濾波處理后的大豆植株點云Fig. 3 Point clouds of soybean plant after filtering

1.3 器官分割

1.3.1莖葉分割 大豆植株葉片相對扁平，同一點在合理尺度范圍的法線方向變化不大，但其莖稈近似于圓柱形，不同尺度法向量會產(chǎn)生較大差異，若法線差大于設定閾值，則認為該點屬于莖稈，反之該點屬于葉片部分。因此，可以利用不同尺度下法向量特征的差異性提取出大豆植株莖稈部分。利用法線微分差異（difference of normals，DoN）[19]實現(xiàn)莖葉分割，具體步驟為：①計算輸入點云數(shù)據(jù)的平均距離（Rm），作為較大支持半徑（Rl）以及較小支持半徑（Rs）的取值依據(jù)；②對于輸入點云數(shù)據(jù)中的每個點（P），分別利用較大的支持半徑（Rl）與較小支持半徑（Rs）計算其法向量，根據(jù)式（1）進行差分運算得到該點差分法向量；③過濾向量域，分割出大豆植株莖稈點云；④基于歐式距離聚類，合并閾值內(nèi)的冠層點云。部分葉片上的點受成像精度的影響以及閾值的選取可能被錯誤地歸為莖稈（圖4C）而被移除。同時，一些孤立的莖部，特別是靠近葉片的部分，可能被錯誤地歸類在葉片點云（圖4B）中。為了提高后期表型參數(shù)測量的精確度，利用歐式聚類移除多余的莖部點云（圖4D）并且回填漏掉的葉片點云（圖4E），得到較為完整的冠層結構（圖4F）。

圖4 大豆植株莖葉分割Fig. 4 Stem and leaf segmentation of soybean plant

式中，P表示點云數(shù)據(jù)表示大半徑及大半徑下的法向量表示小半徑及小半徑下的法向量，?n?表示差分法向量。

1.3.2單葉分割 由于葉片點云密度不均勻且存在遮擋、粘連等情況，去除莖稈的大豆植株冠層葉片點云仍需進一步分割才能實現(xiàn)高效準確的表型參數(shù)分析。單葉分割的難點在于如何識別并斷開葉片粘連部位。對于每個點，根據(jù)式（2）計算其協(xié)方差矩陣(C)的特征向量和特征值，根據(jù)式（3）進一步運算估計得到曲率（δ）。

式中，Pi表示點云數(shù)據(jù)，k表示Pi鄰近點的數(shù)目表示鄰近點的三維質心，T表示點云矩陣的轉置，λj是協(xié)方差矩陣的第j個特征值是第j個特征向量；λ0、λ1、λ2為每個點云數(shù)據(jù)協(xié)方差矩陣的3個特征值。

葉片粘連部位的點分布相對混亂且具有較大的曲率值（圖5A），剔除曲率值大于某個閾值的點，粘連葉片能順利地相互分離。去除葉片粘連點可能會導致葉片邊緣結構缺損，從而影響分割之后得到的單片葉片（圖5B）。因此，本研究利用歐式聚類回填部分粘連點，綜合考慮曲率、法向量夾角以及點間距離，采用改進的區(qū)域生長算法實現(xiàn)冠層葉片分割，得到較完整的單片葉片（圖5C）。

圖5 大豆植株冠層葉片分割Fig. 5 Leaves segmentation of soybean canopy

1.4 表型參數(shù)提取

1.4.1葉面積 基于拉普拉斯平滑（laplace smoothing）的德洛內(nèi)三角剖分（delaunay triangulation）法（簡稱LDT 法）計算大豆葉片葉面積。拉普拉斯平滑是一種網(wǎng)格平滑算法，其原理是在一定的松弛度（r）下將每個點用其鄰域點的中心來代替（式4），通過不斷迭代，得到較為光滑的網(wǎng)格。隨著平滑過程中迭代次數(shù)及松弛系數(shù)增大，平滑程度越大，網(wǎng)格會不斷向中心收縮。為了避免過度收縮造成的誤差，本研究在r=0.2、iteration=20的情況下對葉片進行平滑處理。平滑后利用三角剖分將散亂的點云數(shù)據(jù)剖分為一系列三角網(wǎng)格（圖6E），通過計算所有三角網(wǎng)格面積之和（式5）來代替葉面積。

圖6 大豆表型參數(shù)測量Fig. 6 Measurement of soybean phenotypic parameters

式中，x′表示點的新坐標位置，x表示點的當前坐標位置，j表示當前點的鄰域點個數(shù)，xi表示第i個鄰域點的坐標位置，S表示葉片點云面積，n表示三角網(wǎng)格個數(shù)，Si表示單個三角網(wǎng)格面積。

1.4.2葉寬 分割后的大豆葉片仍處于植株坐標系下（圖6A），此時通過每個坐標軸最大點和最小點建立起來的包圍盒稱為軸對齊包圍盒（axisaligned bounding box，AABB），其空間冗余大，對葉片的表型特征描述誤差就大。為了更加貼近葉片真實的形狀，基于主成分分析構建葉片有向包圍盒（oriented bounding box，OBB）（圖6B），提取包圍盒的寬代表葉片的寬度（圖6C）。具體步驟為：①利用主成分分析獲得葉片點云的質心及3 個主方向，將點云轉換至葉片坐標系下；②基于獲得的主方向和質心，將輸入的葉片點云轉換至原點，建立變換到原點的葉片點云包圍盒；③計算y方向上最大點和最小點間的距離。

1.4.3葉長 本研究采用提取點集以擬合葉片中脈的方法，能夠提取葉片點云上近似于中脈的最短曲線（圖6F）。具體步驟為：①在葉片坐標系下，找到葉片點云空間中距離最遠的2 個點并分別定義為起點和終點；②將起點作為基點，基于最鄰近（K-NearestNeighbor，KNN）算法搜索其鄰近點，找出鄰近點中距基點與終點距離之和最小的點作為下個基點；③重復搜索，直到基點與終點重合，此時找出的所有基點構成的曲線長度就是葉片長度。

1.4.4葉傾角 葉傾角是葉片腹面的法線與天頂軸的夾角，實際上它也是葉面與XOY平面的夾角。在得到葉片有向包圍盒后，進一步將其逆變換至原始植株坐標系下，此時包圍盒Z軸與天頂軸之間的夾角即為葉傾角（圖6C）。

1.4.5莖粗 本研究統(tǒng)一將花盆口上方5 cm 處的平面上莖稈點云歐式距離最遠的2 個點定義為最大點與最小點，將2 點水平方向上的差值作為大豆植株的莖粗。利用直通濾波截取目標莖稈點云，莖粗測量如圖6D所示。

1.4.6表型參數(shù)真實值測量 針對葉片表型參數(shù)，隨機選取10株大豆植株中的30片葉片用來測試所提出的表型參數(shù)測量方法。葉面積、葉寬、葉長和葉傾角均為30個樣本，莖粗為10個樣本。對涉及到的表型參數(shù)真實值作以下約定。

葉面積：將利用標定法[20]對獲取的葉片二維圖像進行像素轉換得到的值作為葉面積真實值。

葉寬：將用卷尺測量平鋪葉片的最大寬度距離值作為葉寬真實值。

葉長：將用卷尺測量平鋪葉片葉枕到葉尖的距離值作為葉長真實值。

葉傾角：將葉片腹面與水平線間夾角值作為葉傾角真實值。

莖粗：將用游標卡尺測量花盆口平面上方5 cm處寬度距離值作為莖粗真實值。

1.5 評價指標

莖葉分割的主要目的在于分離大豆植株莖稈與冠層葉片，以分割后冠層葉片點云數(shù)量評價其分割效果。冠層葉片點云數(shù)量用Nc表示，分割前點云總數(shù)用N表示，基于以上統(tǒng)計根據(jù)式（6）計算分割后冠層葉片點云分割率（Rc）。

單葉分割的主要目的在于分離大豆植株冠層中粘連的葉片，以分割后葉片聚類點云總數(shù)評價其分割效果。葉片聚類點云總數(shù)用Nl表示，分割前點云總數(shù)即冠層葉片點云數(shù)Nc，基于以上統(tǒng)計根據(jù)式（7）計算分割后單葉點云分割率（Rl）。

采用線性回歸分析方法評價表型參數(shù)（葉面積、葉寬、葉長、葉傾角和莖粗）測量結果，以本研究算法測量值與人工實測值之間的平均相對誤差（mean relative error，MRE）、均方根誤差（root mean square error，RMSE）和決定系數(shù)R2為指標，量化各參數(shù)的精度與誤差。

2 結果與分析

2.1 器官分割結果

2.1.1莖葉分割結果 研究表明，當Rs>Rm且Rl≤2Rs的情況下才能體現(xiàn)不同尺度下法向量差異。由圖7 可知，當尺度比例較小時（Rl=2Rs），冠層葉片與莖稈差異不明顯，無法分離出冠層葉片；尺度比例較大時（Rl=5Rs），大量葉片點云錯誤歸類為莖稈點云；當Rl=3Rs時，葉片部分與莖稈的差異最為明顯。在該條件下，當閾值較小時（threshold=0.03），冠層葉片點云過分割；當閾值較大時（threshold=0.10），冠層葉片點云欠分割。為了能夠最大程度地去除莖稈并且保留冠層葉片部分，試驗均選取閾值為0.08的分割結果進行后續(xù)處理。

圖7 不同參數(shù)下莖葉分割結果Fig. 7 Results of stem and leaf segmentation under different parameters

為了驗證莖葉分割方法的準確性和有效性，將本研究提出的方法（DoN）與隨機采樣一致性法（RANSAC）[15]進行比較。表1 為不同方法下植株莖葉分割結果統(tǒng)計，可知本研究方法冠層葉片點云分割率整體高于隨機采樣一致性法，冠層葉片點云平均分割率分別為84.24%和78.34%。

表1 植株莖葉分割結果統(tǒng)計Table 1 Statistics of stem and leaf segmentation results of plants

2.1.2單葉分割結果 由圖8 可知，當閾值設置過大（number of neighbours = 60，smoothness threshold=15/180π，curvature threshold=1）時，葉片不能完全分割并且有大量點歸類于未分類的紅色點集中（圖8A）；隨著閾值減?。╪umber of neighbours=30， smoothness threshold=12/180π，curvature threshold=0.05），大量未分類的點能夠重新歸類，葉片分割效果最好（圖8B），在該條件下基于歐式聚類回填部分點云進而得到完整的葉片點云（圖8C）。表2為植株葉片分割結果統(tǒng)計，利用本研究算法實現(xiàn)葉片分割后，葉片單葉點云分割率均達到95%以上，單葉點云平均分割率為96.39%。

表2 植株葉片分割結果統(tǒng)計Table 2 Statistics of leaf segmentation results of plants

圖8 不同參數(shù)下葉片分割結果Fig. 8 Results of leaf segmentation under different parameters

2.2 表型參數(shù)測量結果評價分析

為了驗證葉面積計算方法的準確性和有效性，將本研究提出的方法（LDT）與投影法（projection）[21]和貪婪投影三角化法（GPT）[22]進行比較。將大豆葉片葉面積測量值與真實值進行擬合，圖9 顯示了3 種葉面積測量方法的回歸比較，其中本研究方法的回歸線最接近標準回歸線y=x，回歸方程為y=0.971 1x+0.271 7，決定系數(shù)R2=0.987 9, 算法測量值與真實值接近，表明本研究方法可以應用于大豆葉片葉面積的測量。通過表3 可知，在測量葉片面積時，投影法平均相對誤差最大，決定系數(shù)R2和均方根誤差分別為0.978 7和1.251 1 cm2；貪婪投影三角化法平均相對誤差低于投影法但高于本研究方法，決定系數(shù)R2和均方根誤差分別為0.979 1 和0.763 3 cm2；本研究提出的方法平均相對誤差最小，決定系數(shù)R2和均方根誤差分別為0.987 9 和0.541 7 cm2，表明本研究方法整體優(yōu)于其他2種方法。

表3 不同方法測量葉面積結果評價Table 3 Result evaluation of leaf area by different methods

圖9 大豆葉面積不同方法測量值與真實值比較Fig. 9 Comparison of soybean leaf area by different method

通過表4 可知，葉寬、葉長及莖粗算法測量值與真實值的平均相對誤差分別為2.68%、2.87%和3.99%，決定系數(shù)R2分別為0.961 3、0.962 6 和0.963 4，均方根誤差分別為0.142 2、0.175 5 和0.047 5 cm。將葉寬、葉長及莖粗算法測量值與真實值進行擬合，擬合結果如圖10A、B和C所示，回歸方程分別為y=0.984 8x+0.086 3、y=0.897 0x+0.424 1 和y=0.830 2x+0.215 5。算法測量值與真實值接近，表明本研究方法可以應用于大豆葉片葉寬、葉長及莖粗的測量。

表4 葉寬、葉長及莖粗測量結果評價Table 4 Result evaluation of leaf width， leaf length and stem diameter

圖10 大豆表型參數(shù)測量值與真實值比較Fig. 10 Comparison of soybean phenotypic parameters extracted with actual values

葉傾角算法測量值與真實值的平均相對誤差為7.22%，決定系數(shù)R2和均方根誤差分別為0.931 1 和3.279 6°。將葉傾角算法測量值與真實值進行擬合，擬合結果如圖10D 所示，回歸方程為y=0.986 0x+0.442 2。葉傾角算法測量值與真實值接近，但其平均相對誤差相較于其余表型參數(shù)而言偏高，表明葉傾角測量算法在應用于大豆植株葉傾角測量方面仍有較大改進空間。

2.3 時效統(tǒng)計分析

本研究采用Intel Core i7-5930K 處理器以及NVIDIA Quadro P4000 GPU 進行試驗，不同處理階段時效統(tǒng)計結果見表5。每株大豆植株點云重建平均用時為24.17 min，濾波處理、莖葉分割、單葉分割以及表型參數(shù)提取平均用時分別為7.35、14.63、8.42 和46.72 s。點云重建受場景以及特征點數(shù)量的影響耗時較長，在匹配特征點時通過縮小目標區(qū)域能夠進一步提高重建效率。

表5 時效分析結果Table 5 Results of processing time analysis

3 討論

基于三維點云的植株表型參數(shù)測量為表型數(shù)據(jù)獲取提供了有效的解決方案。相對于傳統(tǒng)的常規(guī)測量方法，基于三維點云的測量技術具有非接觸、速度快等諸多優(yōu)點。相較于二維圖像，三維點云數(shù)據(jù)包含的深度信息能夠使植株與無關場景天然解耦，很大程度上避免了植株本身形變造成的誤差。已有研究探討了基于植株三維點云測量表型參數(shù)的實踐[8-9,13-15]，并取得了良好的效果，但仍需要大量人為干預。本研究提出的器官分割及表型分析方法與已有研究[14-15]相比，在一定程度上解決了葉片相互粘連造成的莖葉分割及高通量表型測量困難等問題，提高了表型分析的效率和準確度。通過法向量差異實現(xiàn)莖葉分割，并結合點云曲率特征對冠層葉片實施了基于歐式聚類的區(qū)域生長分割，以上器官尺度的分割使植株表型參數(shù)（葉面積、葉寬、葉長、葉傾角及莖粗）測量更加自動、快速、可靠。結果顯示，器官分割后冠層葉片點云平均分割率為84.24%，單葉點云分割率均高于95%，葉面積、葉寬、葉長、葉傾角和莖粗參數(shù)測量值與人工實測值的決定系數(shù)分別為0.987 9、0.961 3、0.962 6、0.931 1 和0.963 4，均方根誤差分別為0.541 7 cm2、0.141 2 cm、 0.175 5 cm、3.279 6°和0.047 5 cm。每株大豆植株點云重建和濾波處理、莖葉分割、單葉分割以及表型參數(shù)提取平均用時分別為24.17 min 和7.35、14.63、8.42 和46.72 s。該結果表明，本研究方法能夠快速地實現(xiàn)大豆植株器官分割，為多分枝作物高通量自動化表型參數(shù)測量提供高效的技術支持。

本研究方法可以類推應用于其他冠層密度大，葉片較為平坦，葉片呈近圓形、橢圓形的多分枝植株。在后續(xù)的研究中，一方面應考慮冠層密度小但葉片卷曲，葉片呈披針形的植株，提高器官分割及表型測量算法應用于不同形態(tài)植株上的精度。另一方面，針對大豆開花結莢期分枝和葉片較多、葉片遮擋嚴重的情況，未來應根據(jù)植株自相似性考慮從部分到整體的三維重建與表型分析。