基于骨架和最優(yōu)傳輸距離的玉米點(diǎn)云莖葉分割和表型提取

朱 超，苗 騰，許童羽，李 娜，鄧寒冰，周云成

（1. 沈陽(yáng)農(nóng)業(yè)大學(xué)信息與電氣工程學(xué)院，沈陽(yáng) 110866；2. 遼寧省農(nóng)業(yè)信息化工程技術(shù)研究中心，沈陽(yáng) 110866）

0 引 言

玉米（Zea mays L.）是世界上重要的糧食作物，其產(chǎn)量對(duì)保障全球糧食供應(yīng)至關(guān)重要。表型分析可輔助玉米育種，提高玉米產(chǎn)量[1]。傳統(tǒng)表型參數(shù)獲取方法主要是人工測(cè)量，效率低、主觀性強(qiáng)、破壞性強(qiáng)[2]。近年來(lái)，計(jì)算機(jī)視覺技術(shù)的發(fā)展為作物表型參數(shù)的自動(dòng)測(cè)量提供了條件[3-4]。

基于二維數(shù)字圖像的視覺技術(shù)率先被應(yīng)用于植物表型參數(shù)的提取[5-7]。Granier等[8]設(shè)計(jì)了Phonopsis成像系統(tǒng)，通過(guò)圖像提取葉片和根的生長(zhǎng)信息來(lái)監(jiān)測(cè)植物對(duì)土壤水分的響應(yīng)。Pan等[9]使用高分辨率的數(shù)碼相機(jī)圖像來(lái)估計(jì)作物覆蓋率。Tanabata等[10]利用圖像分析方法識(shí)別水稻種子輪廓并計(jì)算種子長(zhǎng)、寬、周長(zhǎng)等性狀參數(shù)。

隨著三維激光掃描儀[11]、飛行時(shí)間相機(jī)[12]以及激光雷達(dá)[13]等三維傳感器的普及，基于三維點(diǎn)云的視覺技術(shù)被用于作物表型研究，如作物識(shí)別[14]、長(zhǎng)勢(shì)監(jiān)測(cè)[15]、冠層重建[16-17]以及各類表型參數(shù)的測(cè)量[18-21]等，這些研究證明了從三維空間提取植物表型具有更高的穩(wěn)健性和準(zhǔn)確性[22-23]。為了提取器官尺度的表型參數(shù)，需要在植株點(diǎn)云上進(jìn)行器官分割。Wu等[24]提取玉米植株的三維骨架，并在骨架上實(shí)現(xiàn)器官分割，進(jìn)而提取葉、莖尺度的表型參數(shù)。Elnashef等[25]提出了一種基于張量的三維點(diǎn)云分類與分割算法，利用一階和二階張量分析實(shí)現(xiàn)莖、葉點(diǎn)云分類，再采用基于密度的空間聚類算法[26]對(duì)葉片進(jìn)行分割，但這種聚類方法不能處理頂部新葉點(diǎn)云緊密挨在一起的情況。Jin等[27]開發(fā)了基于中值歸一化向量的區(qū)域生長(zhǎng)莖葉分割算法，該算法首先手動(dòng)標(biāo)記出莖底部的點(diǎn)云，然后通過(guò)中值歸一化向量計(jì)算區(qū)域生長(zhǎng)方向，從葉基到葉尖依次對(duì)每個(gè)葉片進(jìn)行分割。該方法簡(jiǎn)單有效，對(duì)中下部葉片的分割效果很好，但對(duì)于新葉，容易出現(xiàn)誤分割現(xiàn)象。

目前，玉米植株點(diǎn)云的莖葉分割方法對(duì)于完全展開的葉片是有效的，但是對(duì)于間距小或相互包裹在一起的新生葉，分割效果仍然較差。針對(duì)該問(wèn)題，本研究提出一種基于骨架和最優(yōu)傳輸距離的莖葉點(diǎn)云分割方法，嘗試自動(dòng)對(duì)植株頂部的新生葉進(jìn)行準(zhǔn)確分割，為玉米高通量表型檢測(cè)、三維幾何重建提供技術(shù)手段。

1 材料與方法

1.1 點(diǎn)云數(shù)據(jù)獲取與預(yù)處理

田間試驗(yàn)于2019年5月至7月在沈陽(yáng)農(nóng)業(yè)大學(xué)玉米試驗(yàn)田進(jìn)行。選取種植行距60 cm、株距25 cm的苗期、拔節(jié)期玉米為研究對(duì)象，使用FreeScan X3手持激光掃描儀進(jìn)行三維點(diǎn)云獲取，該激光掃描儀的主要參數(shù)如表1所示。

數(shù)據(jù)采集時(shí)，為防止玉米葉片萎蔫，將植株移植到盆中后立即澆水，隨即將玉米轉(zhuǎn)移到室內(nèi)進(jìn)行點(diǎn)云數(shù)據(jù)采集。由于FreeScan X3掃描儀需要借助激光感光片進(jìn)行三維定位，因此數(shù)據(jù)采集時(shí)，將粘貼有感光片的多個(gè)支架圍在植株四周，保證激光線能同時(shí)掃描到植株和感光片。為提高數(shù)據(jù)獲取質(zhì)量，采集時(shí)將掃描儀與玉米植株之間的距離控制在300 mm之內(nèi)。整個(gè)數(shù)據(jù)獲取過(guò)程如圖 1a所示。掃描得到的原始點(diǎn)云數(shù)據(jù)中除了植株外，還包含了部分支架和花盆的點(diǎn)云（圖1b）。使用三維點(diǎn)云處理軟件CloudCompare將支架點(diǎn)和花盆點(diǎn)云手動(dòng)去除，并進(jìn)行下采樣處理，得到最終的玉米植株點(diǎn)云（圖1c）。下采樣后，每個(gè)植株點(diǎn)云包含5 000～8 000個(gè)點(diǎn)，該數(shù)量的點(diǎn)云即可保留玉米植株的形態(tài)特征，也能對(duì)骨架進(jìn)行高效提取。

表1 FreeScan X3激光掃描儀主要參數(shù)Table 1 Main parameters of FreeScan X3 laser scanner

1.2 莖葉分割方法

本研究提出的莖葉分割方法以植株散亂點(diǎn)云為輸入，通過(guò)骨架提取、骨架分解、點(diǎn)云坐標(biāo)系變換、粗分割、精細(xì)分割和表型參數(shù)提取6個(gè)主要步驟，實(shí)現(xiàn)對(duì)莖葉器官的準(zhǔn)確分割和表型參數(shù)的自動(dòng)測(cè)量。

1.2.1 骨架提取

采用拉普拉斯骨架提取算法[28]生成玉米植株骨架。給定一個(gè)玉米植株點(diǎn)云集合P，首先找到P中每個(gè)點(diǎn)的K近鄰點(diǎn)計(jì)算拉普拉斯坐標(biāo)，并利用拉普拉斯收縮操作將P收縮成零體積點(diǎn)云集合C，K為用戶設(shè)定的最近鄰參數(shù)；其次，使用最遠(yuǎn)距離采樣策略從C中采樣得到骨架頂點(diǎn)集合U；然后將U中所有相鄰的骨架頂點(diǎn)相連，任意2個(gè)相鄰頂點(diǎn)均會(huì)生成1條邊，所有邊形成的邊集合用符號(hào)E′表示；最后，使用邊塌縮操作從邊集合E′中去除冗余邊，形成最終的邊集合E。由骨架頂點(diǎn)集合U和邊集合E表示的無(wú)向圖G=(U,E)，即為提取出的植株骨架。該骨架提取算法還為U中的每個(gè)骨架頂點(diǎn)保存了一個(gè)對(duì)應(yīng)的點(diǎn)云集合。使用ui表示U中第i個(gè)（i＞0）骨架頂點(diǎn)，其對(duì)應(yīng)的點(diǎn)云集合用Ci表示，ui由Ci中的點(diǎn)云通過(guò)拉普拉斯收縮和最遠(yuǎn)距離采樣操作獲得。每個(gè)對(duì)應(yīng)點(diǎn)云集合Ci都是植株點(diǎn)云集合P的子集，且C1∪C2∪…CN=P，N為骨架頂點(diǎn)的個(gè)數(shù)。拉普拉斯骨架提取算法生成的玉米植株骨架如圖2所示，整個(gè)骨架由邊和骨架頂點(diǎn)構(gòu)成。稱1條邊上的2個(gè)頂點(diǎn)互為鄰居，根據(jù)頂點(diǎn)鄰居數(shù)量的不同，可將骨架頂點(diǎn)分為3類：

1）根頂點(diǎn)：僅有1個(gè)鄰居的骨架頂點(diǎn)。

2）連接頂點(diǎn)：有2個(gè)以上鄰居的骨架頂點(diǎn)。

3）分支頂點(diǎn)：只有2個(gè)鄰居的骨架頂點(diǎn)。

拉普拉斯骨架提取算法中的最近鄰參數(shù)K會(huì)顯著影響植株骨架的提取結(jié)果。K值越大，骨架越光滑，抗噪聲能力越強(qiáng)，但對(duì)頂部新生葉局部骨架的提取能力減弱。但過(guò)小的K值會(huì)增加冗余邊，需要更多次的邊塌縮操作來(lái)進(jìn)行刪除，從而減慢骨架提取算法的處理速度。通過(guò)試驗(yàn)，本研究明確了對(duì)參數(shù)K的2種設(shè)定方案。第一種方案針對(duì)頂部葉片間距較大的植株，K參數(shù)設(shè)為16；第二種針對(duì)頂部葉片間距較小的植株，K參數(shù)設(shè)為6。2種K參數(shù)設(shè)定方案對(duì)玉米植株骨架提取效果的影響如圖3所示。植株頂部2個(gè)葉片間距較小時(shí)，如果采用第一種K參數(shù)設(shè)定方案，則其中1個(gè)葉片的骨架丟失（圖3a），而采用第二種方案，新葉上半部的局部骨架被正確提?。▓D3b）。植株頂部葉片間距較大時(shí)，無(wú)論采用哪種方案，頂部葉片骨架均能正確提取（圖3c和圖3d）。

1.2.2 骨架分解

得到植株骨架之后，需將其分解為莖子骨架和葉片子骨架。從圖2中的植株骨架拓?fù)浣Y(jié)構(gòu)可以看出，在根頂點(diǎn)中，存在唯一根頂點(diǎn)由莖點(diǎn)云生成（稱為莖根頂點(diǎn)），其余根頂點(diǎn)均由葉尖點(diǎn)云生成（稱為葉根頂點(diǎn)）。只需找到莖根頂點(diǎn)，即可較容易地識(shí)別出莖子骨架，進(jìn)而提取葉片子骨架，因此莖根頂點(diǎn)是后續(xù)骨架分解的基礎(chǔ)。

很多研究利用深度相機(jī)獲取植株點(diǎn)云數(shù)據(jù)。在數(shù)據(jù)獲取時(shí)可通過(guò)調(diào)整相機(jī)位置使植株點(diǎn)云所在坐標(biāo)系的某個(gè)坐標(biāo)軸與植株生長(zhǎng)方向重合，因此通過(guò)坐標(biāo)值即可較容易地確定莖點(diǎn)云的位置[29-30]。但是手持激光三維掃描儀獲取的點(diǎn)云數(shù)據(jù)屬于散亂點(diǎn)云，點(diǎn)云與所在空間坐標(biāo)系間的關(guān)系是隨機(jī)的，因此植株生長(zhǎng)方向并未與某個(gè)坐標(biāo)軸重合，這增大了莖根頂點(diǎn)的識(shí)別難度。本研究利用葉片形態(tài)特征對(duì)莖根頂點(diǎn)進(jìn)行識(shí)別，具體步驟如下：

1）提取根連子骨架

采用迪杰斯特拉算法[31]提取每個(gè)根頂點(diǎn)到其拓?fù)渚嚯x最近的連接頂點(diǎn)間的最短路徑，將最短路徑形成的局部骨架稱為根連子骨架。隨后對(duì)每個(gè)根連子骨架中的頂點(diǎn)進(jìn)行重排序。排序時(shí)，以根頂點(diǎn)作為子骨架的第一個(gè)頂點(diǎn)，并按邊的連接順序依次排列剩余頂點(diǎn)。

2）計(jì)算根連子骨架的單位方向向量

對(duì)于排序后的任一根連子骨架，選擇中間位置的頂點(diǎn)和末端頂點(diǎn)，兩者的空間坐標(biāo)分別用符號(hào)ν1和ν2表示。利用上述2個(gè)頂點(diǎn)計(jì)算根連子骨架的單位方向向量空間坐標(biāo)V，計(jì)算如式（1）所示

3）利用葉片形態(tài)特征識(shí)別莖根頂點(diǎn)

觀察發(fā)現(xiàn)，正常玉米葉片中間到葉基的區(qū)域均朝向天空方向，因此葉片根連子骨架的方向向量也會(huì)朝向天空，相反莖根連子骨架的方向向量會(huì)朝向地面，2類子骨架方向向量的示意圖如圖4所示。其中，2個(gè)葉片根連子骨架的方向向量間的夾角會(huì)小于π/2；而莖根連子骨架的方向向量與葉片根連子骨架的方向向量的夾角會(huì)大于π/2。因此通過(guò)計(jì)算不同根連子骨架方向向量之間的夾角即可識(shí)別莖根頂點(diǎn)，具體計(jì)算方法如下。

用符號(hào)Vτ表示第τ個(gè)根連子骨架的單位方向向量，用符號(hào)Vφ表示第φ個(gè)根連子骨架的單位方向向量，則2個(gè)子骨架之間的夾角（θτφ，rad）可通過(guò)式（2）計(jì)算。

利用式（2）計(jì)算每個(gè)根連子骨架與其他所有根連子骨架的夾角，并求其平均值，平均值大于π/2的子骨架即為莖根連子骨架，該子骨架中的根頂點(diǎn)即為莖根頂點(diǎn)，其余根連子骨架的根頂點(diǎn)均為葉根頂點(diǎn)。

識(shí)別出根頂點(diǎn)的類型后，將植物骨架分解成多個(gè)葉片子骨架和1個(gè)莖子骨架。在植物骨架無(wú)向圖G中，葉片根連子骨架直接作為葉片子骨架；采用迪杰斯特拉算法[31]提取莖根頂點(diǎn)到其最遠(yuǎn)連接頂點(diǎn)的最短路徑作為莖子骨架，并對(duì)莖子骨架的頂點(diǎn)進(jìn)行排序，以莖根頂點(diǎn)作為子骨架的第1個(gè)頂點(diǎn)，并按邊的連接順序依次排列剩余頂點(diǎn)。

1.2.3 點(diǎn)云坐標(biāo)系變換

三維掃描儀獲取的散亂點(diǎn)云數(shù)據(jù)與所在空間坐標(biāo)系間的關(guān)系是隨機(jī)的，不利于通過(guò)空間坐標(biāo)值來(lái)判斷點(diǎn)在植株上的相對(duì)位置，因此本研究構(gòu)建了一個(gè)更符合玉米形態(tài)特征的植株坐標(biāo)系，并將點(diǎn)云坐標(biāo)從原始坐標(biāo)系轉(zhuǎn)換到植株坐標(biāo)系中。利用骨架分解得到的莖子骨架構(gòu)建植株坐標(biāo)系，具體包含以下3個(gè)步驟：

1）采用隨機(jī)抽樣一致算法[32]對(duì)莖子骨架中的骨架頂點(diǎn)進(jìn)行單位方向向量擬合，該向量用符號(hào)L1表示，L1的負(fù)方向指向莖根頂點(diǎn)。

2）將植株點(diǎn)云投影到與L1垂直的二維平面獲得二維投影點(diǎn)，采用主成分分析算法[33]計(jì)算投影點(diǎn)的第一主成分向量和第二主成分向量，分別用符號(hào)L2和L3表示。

3）構(gòu)建植株坐標(biāo)系：莖根頂點(diǎn)作為植株坐標(biāo)系的原點(diǎn)O；方向向量L2、L3和L1分別作為植株坐標(biāo)系的X軸、Y軸和Z軸。

得到植株坐標(biāo)系后，將植株點(diǎn)云和骨架頂點(diǎn)的坐標(biāo)值從原始坐標(biāo)系轉(zhuǎn)換到植株坐標(biāo)系。對(duì)于任意三維點(diǎn)p，設(shè)其在原始坐標(biāo)系下的三維坐標(biāo)為（p(x)’，p(y)’，p(z)’），在植株坐標(biāo)系下的坐標(biāo)為（p(x)，p(y)，p(z)），2組坐標(biāo)值的轉(zhuǎn)換可通過(guò)式（3）計(jì)算得到。

式中L1(x)、L1(y)、L1(z)分別表示原始坐標(biāo)系下，方向向量L1的X、Y、Z軸坐標(biāo)值；同理，L2(x)、L2(y)、L2(z)表示L2在原始坐標(biāo)系下的X、Y、Z軸坐標(biāo)值；L3(x)、L3(y)、L3(z)表示L3在原始坐標(biāo)系下的X、Y、Z軸坐標(biāo)值；O(x)、O(y)、O(z)表示植株坐標(biāo)系原點(diǎn)O在原始坐標(biāo)系下的X、Y、Z軸坐標(biāo)值。

在植株坐標(biāo)系下，可通過(guò)Z軸坐標(biāo)來(lái)明確一個(gè)點(diǎn)在植株上的高度，Z坐標(biāo)越大，該點(diǎn)在植株上所處的位置越高。本研究后續(xù)的分割算法均在植株坐標(biāo)系中進(jìn)行。

1.2.4 粗分割

為了便于描述后續(xù)點(diǎn)云分割過(guò)程，定義?u表示未分割的點(diǎn)云集合；定義A0表示已分割的莖點(diǎn)云集合；定義Ak(1≤k≤m）表示已分割的葉片點(diǎn)云集合，其中m表示葉片子骨架的數(shù)量。將莖點(diǎn)云集合和葉片點(diǎn)云集合統(tǒng)稱為器官點(diǎn)云集合。

利用分解出的莖、葉子骨架進(jìn)行粗分割。在骨架提取時(shí)，已為每個(gè)骨架頂點(diǎn)ui均保存了一個(gè)對(duì)應(yīng)點(diǎn)云集合Ci，借助這些對(duì)應(yīng)點(diǎn)云集合可以實(shí)現(xiàn)莖葉器官的粗分割。粗分割的具體步驟如下：

1）葉片點(diǎn)云粗分割

對(duì)于每個(gè)葉片子骨架，提取該子骨架中的所有骨架頂點(diǎn)，并將這些骨架頂點(diǎn)的對(duì)應(yīng)點(diǎn)云集合合并，形成一個(gè)葉片點(diǎn)云集合。每個(gè)葉片子骨架均會(huì)分割出一個(gè)葉片點(diǎn)云集合，因此當(dāng)葉片子骨架個(gè)數(shù)為m時(shí)，最終分割得到的葉片點(diǎn)云集合數(shù)量也為m。

2）莖點(diǎn)云粗分割

將莖子骨架中骨架頂點(diǎn)的對(duì)應(yīng)點(diǎn)云集合合并，形成唯一的莖點(diǎn)云集合A0。

莖葉粗分割結(jié)果如圖5a所示，其中，大量葉基部的點(diǎn)云被誤分割到莖中，因此需要進(jìn)行精細(xì)分割將這些誤分割點(diǎn)從莖中移除，并分割到正確的葉片中。

1.2.5 精細(xì)分割

首先通過(guò)3個(gè)約束操作從莖點(diǎn)云集合A0中提取非莖點(diǎn)云，具體步驟如下：

1）高度約束

高度約束要求莖點(diǎn)云高度不能過(guò)高，超過(guò)一定高度的點(diǎn)大概率屬于葉片。在莖子骨架中，按照骨架頂點(diǎn)排列順序，倒序遍歷頂點(diǎn)，找到第二個(gè)出現(xiàn)的連接頂點(diǎn)，將該連接頂點(diǎn)的Z坐標(biāo)值作為閾值，與A0中點(diǎn)云的Z坐標(biāo)進(jìn)行比較，高于閾值的點(diǎn)從A0中移除，并加入未分割點(diǎn)云集合?u中。

2）半徑約束

半徑約束要求莖點(diǎn)云與莖中軸的距離不能過(guò)大，超過(guò)一定距離的點(diǎn)大概率屬于葉片。首先計(jì)算A0集合中所有點(diǎn)Z坐標(biāo)的中值z(mì)1；再?gòu)腁0中選出Z坐標(biāo)小于z1/2的點(diǎn)擬合莖中軸，并計(jì)算這些點(diǎn)到莖中軸距離的平均值rc。如果A0中某個(gè)點(diǎn)到莖中軸的距離大于αrc（α為用戶設(shè)定的半徑參數(shù)，α＞0），則將其從A0中移除，并加入集合?u中。

3）距離約束

距離約束要求莖點(diǎn)云彼此之間的距離不能過(guò)大，超過(guò)一定范圍的點(diǎn)大概率屬于葉片。首先采用K-dTree搜索方法找到A0中每個(gè)點(diǎn)的32個(gè)空間鄰近點(diǎn)，并計(jì)算該點(diǎn)到其所有鄰近點(diǎn)的平均距離；隨后估算所有平均距離的平均值d和方差s。如果某個(gè)點(diǎn)到其32個(gè)空間鄰近點(diǎn)的平均距離大于d+s，則將該點(diǎn)從A0中移除，并加入集合?u中。

通過(guò)上述3個(gè)約束操作，將A0中可能屬于葉片的點(diǎn)轉(zhuǎn)移到?u中，轉(zhuǎn)移之后的結(jié)果如圖5b所示。

接下來(lái)，采用基于最優(yōu)傳輸距離的莖葉分割算法，將?u中的未分割點(diǎn)依次分割到對(duì)應(yīng)的器官點(diǎn)云集合中，實(shí)現(xiàn)精細(xì)分割。莖葉分割算法的具體步驟如下：

1）計(jì)算植株點(diǎn)云的最優(yōu)傳輸距離矩陣

本研究采用最優(yōu)傳輸距離作為兩點(diǎn)之間的距離度量。假設(shè)植株點(diǎn)云P共有n個(gè)點(diǎn)，ph(1≤h≤n)表示P中第h個(gè)點(diǎn)的三維坐標(biāo)，pl(1≤l≤n)表示P中第l個(gè)點(diǎn)的三維坐標(biāo)。ph與pl之間的最優(yōu)傳輸距離dhl可通過(guò)基于Sinkhorn散度的最優(yōu)傳輸算法[34]計(jì)算得到。該算法構(gòu)建的約束方程如式（4）所示。

式中ε為用戶設(shè)定的光滑參數(shù)。

通過(guò)Sinkhorn固定迭代法[35]對(duì)式（4）進(jìn)行求解，得到一個(gè)最優(yōu)傳輸距離矩陣D，該矩陣為n×n方陣，矩陣中的元素為2個(gè)點(diǎn)之間最優(yōu)傳輸距離的倒數(shù)，例如，D中第h行、第l列的數(shù)值即為點(diǎn)ph和點(diǎn)pl之間的最優(yōu)傳輸距離的倒數(shù)1/dhl。因此，得到矩陣D后，即可獲得任意2個(gè)點(diǎn)之間的最優(yōu)傳輸距離。光滑參數(shù)ε的設(shè)定會(huì)影響矩陣D的結(jié)果，ε越大，D矩陣對(duì)角線上的元素值越接近1/n，其余位置元素值越接近0；ε越小，D矩陣對(duì)角線上的元素值越小于1/n，其余位置元素值越大于0。

2）未分割點(diǎn)云排序

為提高新生葉的分割精度，本研究采用從植株上部到下部的順序?qū)?u中的未分割點(diǎn)依次進(jìn)行分割。將?u中的點(diǎn)根據(jù)Z坐標(biāo)值按從大到小的順序進(jìn)行排序，之后依次取出，執(zhí)行步驟3。

3）點(diǎn)分割

從?u中取出1三維點(diǎn)p，計(jì)算其到每個(gè)器官點(diǎn)云集合的距離。三維點(diǎn)p到一個(gè)器官點(diǎn)云集合的距離被定義為p與集合中所有點(diǎn)的最優(yōu)傳輸距離的最小值。比較點(diǎn)p與所有器官點(diǎn)云集合之間的距離，將點(diǎn)p加入到距離最小的器官集合中。循環(huán)執(zhí)行步驟3，直到?u集合為空。精細(xì)分割的結(jié)果如圖5c所示。

本研究采用了最優(yōu)傳輸距離作為兩點(diǎn)之間的距離度量。相比于歐氏距離，在本研究的分割流程中，使用最優(yōu)傳輸距離能夠更好地處理頂部大葉包裹小葉的情況。采用2種距離的頂部葉片分割效果對(duì)比如圖6所示，其中，在含有大葉包裹小葉情況的點(diǎn)云區(qū)域，采用最優(yōu)傳輸距離能夠?qū)⒆钔獠康拇笕~點(diǎn)云正確分割出來(lái)，而歐式距離會(huì)將大葉點(diǎn)云誤分割為小葉點(diǎn)云。

1.2.6 表型參數(shù)提取

利用分割出的葉、莖器官點(diǎn)云提取表型參數(shù)，表型參數(shù)提取方法如圖7所示。

株高參數(shù)（cm）提取值為植株點(diǎn)云Z坐標(biāo)最大值和最小值的差值；冠幅參數(shù)（cm）提取值為植株點(diǎn)云X坐標(biāo)最大值和最小值的差值。

利用莖點(diǎn)云提取莖高（cm）和莖粗（cm）參數(shù)。首先使用隨機(jī)抽樣一致算法擬合莖點(diǎn)云的中軸線段，該線段的長(zhǎng)度作為莖高參數(shù)；之后計(jì)算每個(gè)莖點(diǎn)云到中軸的距離，所有距離的中值的2倍作為莖粗參數(shù)。

利用葉片點(diǎn)云提取葉長(zhǎng)（cm）和葉寬（cm）參數(shù)。首先采用主成分分析方法計(jì)算葉片點(diǎn)云的3個(gè)主成分軸。其中第1主成分軸的2個(gè)端點(diǎn)之間的最短路徑稱為葉長(zhǎng)路徑，將該路徑的長(zhǎng)度作為提取的葉長(zhǎng)參數(shù)。為了提取葉寬參數(shù)，首先沿葉長(zhǎng)路徑將葉點(diǎn)云切分成11個(gè)區(qū)間；然后在中間區(qū)間的點(diǎn)云中找到沿第2主成分軸的2個(gè)端點(diǎn)，并提取兩者之間的最短路徑T1，同時(shí)在這些點(diǎn)云中找到第3主分量軸的2個(gè)端點(diǎn)，并提取兩者之間的最短路徑T2，將T1和T2中的最長(zhǎng)路徑作為葉寬路徑，其長(zhǎng)度即為葉寬參數(shù)。

1.3 評(píng)價(jià)指標(biāo)

本研究選取30株不同形態(tài)的玉米植株點(diǎn)云對(duì)莖葉分割和表型提取方法進(jìn)行定量評(píng)估。

使用精確度（Precision）、召回率（Recall）、微F1分?jǐn)?shù)（Micro F1-score）和總體準(zhǔn)確率（Overall accuracy）對(duì)單個(gè)植株點(diǎn)云的莖葉分割精度進(jìn)行評(píng)估[27]。利用CloudCompare軟件手動(dòng)分割植株點(diǎn)云獲得分割真實(shí)值。假設(shè)1個(gè)玉米植株點(diǎn)云共有n1個(gè)點(diǎn)、M1個(gè)器官，對(duì)于其中第λ個(gè)器官點(diǎn)云集合，用TPλ表示屬于該器官的點(diǎn)云被正確分割為該器官的點(diǎn)數(shù)量，用FPλ表示不屬于該器官的點(diǎn)云被誤分割成該器官的點(diǎn)數(shù)量，用FNλ表示屬于該器官的點(diǎn)云被誤分割為其他器官的點(diǎn)數(shù)量，則該器官點(diǎn)云的器官精確度Pλ和器官召回率Rλ可通過(guò)式（5）和式（6）計(jì)算得出。

玉米植株莖葉分割的精確度、召回率、微F1分?jǐn)?shù)和可通過(guò)式（7）～式（10）進(jìn)行計(jì)算。

采用線性回歸分析評(píng)估表型參數(shù)的人工實(shí)測(cè)值與本研究方法提取值之間的關(guān)系，并使用決定系數(shù)R2（Coefficient of determination）和均方根誤差（Root Mean Square Error，RMSE）進(jìn)行定量化評(píng)估，計(jì)算如式（11）和式（12）所示

式中M2為表型數(shù)據(jù)樣本個(gè)數(shù)；η為表型數(shù)據(jù)樣本標(biāo)號(hào)；qη為表型參數(shù)提取值；ˉqη為表型參數(shù)提取值的平均值；Qη為參數(shù)人工實(shí)測(cè)值；ˉQη為參數(shù)實(shí)測(cè)值的平均值。

2 結(jié)果與討論

2.1 莖葉點(diǎn)云分割

算法在配置為2.2 GHZ CPU、DDR 32 G內(nèi)存的筆記本工作站上進(jìn)行了測(cè)試。算法處理8 000個(gè)點(diǎn)左右的玉米植株，平均處理時(shí)間在16 s，其中70%的時(shí)間用于骨架提取。

2.1.1 算法參數(shù)

影響本研究分割效果的參數(shù)主要為半徑參數(shù)α和光滑參數(shù)ε。

半徑參數(shù)α影響莖點(diǎn)云分割結(jié)果，該值越大，分割得到的莖越粗越高；α值越小，莖越細(xì)越短（圖8）。通過(guò)多次試驗(yàn)發(fā)現(xiàn)，針對(duì)苗期植株點(diǎn)云，α=1.0能夠取得較好的效果；針對(duì)拔節(jié)期植株，α=0.8效果較好。

光滑參數(shù)ε會(huì)影響最優(yōu)傳輸距離矩陣的計(jì)算，進(jìn)而影響精細(xì)分割的結(jié)果。在本研究的精細(xì)分割流程中，相比于歐式距離，采用最優(yōu)傳輸距離可以更好地處理大葉包裹小葉的情況。但隨著光滑參數(shù)ε的增大，采用最優(yōu)傳輸距離的分割結(jié)果會(huì)逐漸趨向于采用歐式距離的結(jié)果，當(dāng)ε＞50后，兩者效果基本相同。通過(guò)試驗(yàn)發(fā)現(xiàn)，在本研究中ε=5時(shí)分割結(jié)果最好。不同光滑參數(shù)ε值下的分割效果如圖9所示，當(dāng)ε=5時(shí)，算法對(duì)頂部新葉的分割結(jié)果要遠(yuǎn)遠(yuǎn)好于ε=50和ε=100的情況；ε=50和ε=100時(shí)的分割結(jié)果基本相同，這是因?yàn)楫?dāng)ε較大時(shí)，采用最優(yōu)傳輸距離的效果已經(jīng)等同于歐式距離。

2.1.2 莖葉分割精度

將本研究莖葉分割方法與現(xiàn)有常用的點(diǎn)云分割方法進(jìn)行對(duì)比，包括區(qū)域增長(zhǎng)法[36]、法線差分法[37]和基于近似凹凸性分析的方法[38]。使用其他方法對(duì)玉米植株點(diǎn)云進(jìn)行分割時(shí)，針對(duì)每個(gè)植株數(shù)據(jù)，均采用窮舉策略來(lái)調(diào)整這些方法中的參數(shù)值，選取總體精度值最大的結(jié)果作為分割結(jié)果參與評(píng)估。除了對(duì)比上述3個(gè)方法外，還在本研究的莖葉分割流程中，將最優(yōu)傳輸距離替換為歐式距離，形成基于骨架和歐式距離的方法參與結(jié)果評(píng)估。5種方法的分割精度對(duì)比結(jié)果如表2所示，本研究提出的基于骨架和最優(yōu)傳輸距離的方法結(jié)果最好，30個(gè)樣本的平均精確率、平均召回率、平均微F1分?jǐn)?shù)和平均總體精度分別為0.967、0.961、0.964和0.967，略高于本研究提出的基于骨架和歐式距離的方法，遠(yuǎn)遠(yuǎn)高于另外3種方法的各類指標(biāo)。結(jié)果說(shuō)明，本研究構(gòu)建的莖葉分割方法流程更加適合進(jìn)行玉米植株點(diǎn)云分割，而最優(yōu)傳輸距離作為點(diǎn)的距離度量要好于歐式距離。

表2 5種方法分割精度的對(duì)比Table 2 Comparison of segmentation accuracy of five methods

本研究分割結(jié)果的可視化示例圖如圖10所示，分別給出了總體精度最高、最低，微F1分?jǐn)?shù)最高、最低，葉片數(shù)量最多、最少，植株高度最高、最低的可視化效果，從可視化結(jié)果也可以看出，本研究分割的結(jié)果與真實(shí)值非常相近。

一些點(diǎn)云莖葉分割算法在去除莖點(diǎn)云之后，采用基于歐式距離的聚類方法[39]或基于密度的空間聚類方法[25]進(jìn)行葉片分割。測(cè)試了本研究提出的基于骨架和最優(yōu)傳輸距離的莖葉分割方法和上述2種聚類方法對(duì)新葉的分割效果。測(cè)試時(shí)，將植株的完全展開葉點(diǎn)云刪除，只保留最上部的4片新葉點(diǎn)云和莖點(diǎn)云，以便更加清晰地對(duì)新葉分割結(jié)果進(jìn)行評(píng)估。由于2種聚類方法不能分割莖，因此在測(cè)試時(shí)刪除莖點(diǎn)云，只對(duì)葉片進(jìn)行分割；但運(yùn)行基于骨架和最優(yōu)傳輸距離的莖葉分割方法時(shí)，則保留莖點(diǎn)云。3種方法對(duì)新葉分割精度的對(duì)比如表3所示，基于骨架和最優(yōu)傳輸距離的莖葉分割方法對(duì)新葉點(diǎn)云的分割能力要優(yōu)于其他2種聚類方法，其平均微F1分?jǐn)?shù)、平均精確率、平均召回率和平均總體準(zhǔn)確率達(dá)到0.956、0.957、0.957和0.964。由3種方法對(duì)新葉分割結(jié)果的可視化圖可知，與手動(dòng)分割的真實(shí)值相比（圖11a），本研究提出的基于骨架和最優(yōu)傳輸距離的莖葉分割方法較好地分割出了4個(gè)新葉（圖11b），而基于歐式距離的聚類方法（圖 11c）和基于密度的空間聚類方法（圖11d）均將最外側(cè)的2個(gè)葉片誤分割成同一個(gè)葉片，這是由于2個(gè)葉片的葉基部點(diǎn)云緊密相連導(dǎo)致的。

2.2 表型參數(shù)提取

各表型參數(shù)提取值與人工實(shí)測(cè)值的對(duì)比如圖12所示，對(duì)比結(jié)果表明，株高提取值和實(shí)測(cè)值之間具有較強(qiáng)的相關(guān)性，決定系數(shù)R2和均方根誤差RMSE分別為0.99和1.71 cm；冠幅的R2和RMSE分別為0.99和3.44 cm；莖高的R2和RMSE分別為0.96和7.07 cm；莖粗的R2和RMSE分別為0.97和0.41 cm；葉寬的R2和RMSE分別為0.93和0.85 cm；葉長(zhǎng)的R2和RMSE分別為0.96和5.28 cm。

株高和冠幅2個(gè)表型參數(shù)的提取結(jié)果與人工實(shí)測(cè)值非常接近，反映出本研究構(gòu)建的植株坐標(biāo)系適合用來(lái)表示植物的形態(tài)結(jié)構(gòu)。

表3 3種方法新葉分割精度的對(duì)比Table 3 Comparison of new leaf segmentation accuracy of three methods

葉寬參數(shù)的誤差主要由下述2個(gè)因素造成：1）如果葉寬路徑的2個(gè)端點(diǎn)未沿葉脈對(duì)稱，會(huì)導(dǎo)致葉寬路徑與葉脈不垂直，造成提取的葉寬參數(shù)偏大；2）本研究使用葉片中部的點(diǎn)云提取葉寬，但某些葉片最寬的區(qū)域并不在葉片中部，導(dǎo)致提取的葉寬參數(shù)偏小。

葉長(zhǎng)的誤差是由于提取的葉長(zhǎng)路徑與葉脈不重合導(dǎo)致的。下述3個(gè)因素會(huì)導(dǎo)致該問(wèn)題：1）如果葉片的葉尖區(qū)域下垂嚴(yán)重，其第一主成分軸的端點(diǎn)就不會(huì)落在葉尖區(qū)域，導(dǎo)致葉長(zhǎng)路徑中缺失葉尖部分的葉脈點(diǎn)云；2）本研究方法對(duì)部分葉片的葉領(lǐng)區(qū)域分割精度低，導(dǎo)致葉長(zhǎng)路徑中缺失葉領(lǐng)點(diǎn)；3）如果莖葉分割算法將莖上的點(diǎn)誤分割為葉片點(diǎn)云，則葉長(zhǎng)路徑中會(huì)包含莖點(diǎn)云。前2個(gè)因素會(huì)導(dǎo)致提取的葉長(zhǎng)參數(shù)偏低，而第3個(gè)因素會(huì)使提取的葉長(zhǎng)參數(shù)偏大。

莖高的誤差來(lái)自于新葉葉基與莖連接區(qū)域點(diǎn)云的誤分割。莖粗的誤差主要來(lái)自以下2個(gè)因素：1）本研究以莖點(diǎn)云到莖中軸距離的中值作為莖半徑，這種計(jì)算方法將莖的橫截面看作是圓，但一些玉米植株的莖橫截面更接近橢圓，導(dǎo)致莖粗的提取值偏低；2）如果莖葉分割算法將完全展開葉的葉領(lǐng)點(diǎn)云誤分割為莖點(diǎn)云，會(huì)導(dǎo)致莖粗提取值偏高。

總體上看，雖然提取的表型參數(shù)存在高估或低估，但與手工測(cè)量值具有非常高的相關(guān)性，因此一定程度上也驗(yàn)證了本研究表型參數(shù)提取方法的穩(wěn)定性和實(shí)用性。

3 結(jié) 論

本研究提出了基于點(diǎn)云骨架和最優(yōu)傳輸距離的玉米莖葉自動(dòng)分割方法。該方法先通過(guò)植株點(diǎn)云骨架的提取和分解實(shí)現(xiàn)對(duì)器官位置、類別的識(shí)別，再采用基于最優(yōu)傳輸距離的精細(xì)分割算法實(shí)現(xiàn)對(duì)莖葉器官的精確分割，最后利用分割出的莖葉點(diǎn)云自動(dòng)提取株高、冠幅、莖高、莖粗、葉長(zhǎng)和葉寬6種表型參數(shù)。通過(guò)分析試驗(yàn)結(jié)果，得出以下結(jié)論：

1）本研究可以對(duì)苗期、拔節(jié)期玉米植株進(jìn)行精準(zhǔn)的莖葉分割，平均精確率、召回率、微F1分?jǐn)?shù)和總體精準(zhǔn)度分別為0.967、0.961、0.964和0.967，高于區(qū)域增長(zhǎng)法、法線差分法和基于近似凹凸性分析的分割方法。

2）本研究對(duì)新生葉的分割效果較好，平均微F1分?jǐn)?shù)、平均精確率、平均召回率和平均總體精度為0.956、0.957、0.957和0.964，高于基于歐式距離的聚類方法和基于密度的空間聚類方法。

3）本研究可自動(dòng)提取6個(gè)表型參數(shù)，株高參數(shù)的決定系數(shù)R2（Coefficient of determination）和均方根誤差（Root Mean Square Error，RMSE）分別為0.99和1.71 cm；冠幅的R2和RMSE分別為0.99和3.44 cm；莖高的R2和RMSE分別為0.96和7.07 cm；莖粗的R2和RMSE分別為0.97和0.41 cm；葉寬的R2和RMSE分別為0.93和0.85 cm；葉長(zhǎng)的R2和RMSE分別為0.96和5.28 cm。

本研究為玉米表型高通量測(cè)量、玉米植株三維重建提供了一種有效技術(shù)手段。下一步將研究對(duì)長(zhǎng)有雄穗、果穗的玉米植株點(diǎn)云進(jìn)行分割，以實(shí)現(xiàn)對(duì)玉米全生育期表型信息的自動(dòng)提取。