梁秀英 周風燃 陳 歡 梁 博 許錫晨 楊萬能
(1.華中農業大學工學院, 武漢 430070; 2.華中農業大學作物遺傳改良國家重點實驗室, 武漢 430070;3.華中農業大學植物科學技術學院, 武漢 430070)
玉米是世界三大糧食作物之一[1-4]。在玉米生長過程中,株高、莖粗、葉面積等形態學參數決定植株對養分的競爭和獲取,進而影響玉米的產量[5-7]。葉投影面積、葉夾角等參數直接影響玉米植株對光能的吸收[8],葉片數直接反映了整株玉米的地上部生物量[9],單株最小包圍盒體積反映了玉米植株的空間緊湊度,單片葉最小包圍盒體積、葉周長、葉投影寬度、葉投影長度等參數直接影響葉片的空間分布[10],而玉米葉片的空間分布狀態變化能夠較大程度地影響玉米植株的光能利用率,進而影響玉米的生物量和產量[11]。因此,及時準確獲取玉米植株的株高、莖粗、葉面積、葉片數、葉夾角等形態學參數是研究玉米植株的必要條件。傳統的玉米植株形態學參數測量主要依靠人工或二維圖像測量。人工測量速度慢、成本高、主觀誤差較大[12-13],并且對植株有損傷,所以不適合玉米植株的連續測量[14];基于二維圖像的測量往往由于植株間的遮擋導致精度不高,甚至由于維度的限制不適合用于測量葉面積等形態學參數[15]。
近幾年,為獲得更高的測量精度,基于三維重建的方法被用于目標性狀提取[16-20]。運動恢復結構(SfM)是一種基于多視圖幾何基本原理的三維重建技術,具備自校正、受環境約束較少、室內室外均能使用等優勢,因此,在三維重建中使用非常普遍[21-23]。馮林等[24]基于無人機采集的傾斜影像與部署的地面控制點,采用SfM算法構建了高精度侵蝕溝表面模型,對其建模精度與數字高程模型、正射影像等進行分析,并與傳統正射航圖建模結果進行比較,證明該方法在侵蝕溝的高精度建模與檢測方面具有顯著優勢。陳鈺等[25]用簡單便攜的非測量智能手機獲取礫巖露頭的圖像序列,基于SfM算法進行三維重建、生成點云數據,結果表明,該方法可以快速準確地對地質體進行重建,為野外地質調查及施工等工作提供一種新的思路。在植物表型測量方面,基于SfM算法的植株三維重建具有重建精度較高,能夠動態、無損、大批量重建研究對象等優點,為植株性狀提取和表型參數測量提供了一種新的方法[26-28]。張慧春等[29]構建了一套機器視覺系統,用于采集擬南芥植株的二維圖像,利用SfM算法生成三維點云數據,測量擬南芥的葉片寬度、長度、主莖長度、葉片面積和葉片夾角等表型參數,與傳統人工接觸式測量值相比,該系統測量的擬南芥葉片寬度、長度、主莖長度、葉片面積、葉片間夾角的平均相對誤差分別為9.830%、10.100%、1.070%、4.090%和4.370%,效率高、速度快,可滿足擬南芥形態表型分析的需要。HUI等[30]基于RGB相機獲取茄子、青椒和黃瓜3種植物不同生長時期的圖像序列,利用SfM算法重建三維點云模型,基于得到的點云數據,提取3種植物的葉片數、葉長、葉寬和葉面積等表型參數,并與利用激光掃描儀得到的三維點云重建結果進行對比,結果表明,基于圖像序列的三維模型與對應的三維激光掃描模型差別較小,其三維重建精度較高。HE等[31]搭建了一套多視角的立體成像系統,用于采集草莓果實360°圖像數據,并基于SfM算法重建三維模型,利用自定義的軟件進行三維點云分割,提取草莓果實的高度、長度、寬度、體積、花萼大小、顏色和瘦果數7個性狀參數,與人工實測數據對比表明,提取的7個性狀參數相關性良好、準確性較高。
本文基于SfM算法重建玉米植株三維點云模型,提出一種自動分割、提取玉米單株、莖稈、葉片的方法。使用前期研制的履帶式拍照小車在戶外自動獲取多視角玉米植株二維圖像,采用SfM算法進行三維重建,生成點云數據,提出玉米三維點云自動分割和性狀提取方法,以實現戶外玉米株高、單株最小包圍盒體積、莖粗、葉片數、葉周長、葉面積、單片葉最小包圍盒體積、葉投影面積、葉投影寬度、葉投影長度、葉夾角等參數的動態、準確、無損測量。
本實驗共選擇3個玉米品種(晉單73、大豐30、鄭單958),每個品種有30株(分3個重復種植,每個重復單獨種植一行,一行每個品種各10株),共90株玉米植株樣本,玉米行間距為90 cm,每行內兩株玉米的初始距離為30 cm。為了數據后期的處理方便,實驗樣本采用盆栽種植,避免植株之間的遮擋。用于種植玉米的土壤先在太陽下曬干,然后將曬干的土壤初步碾碎,過6 mm網篩去除土壤中的石頭和雜草,保證土壤的均一性,最后將尿素、磷酸二氫鉀、氯化鉀肥料按照4∶3∶1的配比攪拌均勻。在玉米生長過程中提供正常的水肥管理,待玉米播種60 d后(玉米抽穗期)開始進行玉米株高、莖粗、葉面積等參數的動態無損測量及人工對比驗證。
由于天氣等外在原因,最終用于實驗的玉米植株樣本共80株。人工測量玉米植株株高、莖粗、葉面積和地上部生物量的具體步驟為:① 3名工作人員分別用塑料長尺以盆沿為零刻度線基準測量同一株玉米植株高度,取3名工作人員讀數的平均值作為該玉米植株的最終人工實測株高。② 3名工作人員用游標卡尺(0~200 mm開式,精度0.01 mm)測量同一株玉米植株的相同部位(距根部10 cm處)并讀取數值,取3名人員讀數的平均值作為該玉米植株的人工實測莖粗。③用葉面積儀(YMJ-C/YMJC型)測量玉米葉片面積。④用電子天平(LT2002T型)稱量玉米植株的地上部生物量。
為全自動地在戶外獲取玉米植株二維圖像,本實驗使用前期研制的能按規定路線自動行走的履帶式拍照小車采集不同視角下的玉米植株圖像。小車主要由計算機系統、圖像采集系統、可編程邏輯控制器(Programmable logic controller,PLC)系統、磁導航系統、供電系統、履帶式底盤結構等組成,如圖1所示。為達到比較好的重建效果,要求相鄰兩幅圖像的重疊率在75%以上[13],這要求相機的移動速度(即履帶式拍照小車的速度)不能太快,通過測試選取小車速度為0.1 m/s。相鄰兩幅圖像的拍攝間隔為5~6 cm。圖像采集系統由兩個可見光相機組成,分別用來采集側視圖像和俯視圖像,俯視相機鏡頭與水平方向夾角為90°,與地面的垂直距離為150 cm,側視相機鏡頭與水平方向夾角為30°,與地面的垂直距離為110 cm,相機型號均為Canon EOS 77D,鏡頭為EF-S18-200mm f/3.5-5.6 IS,分辨率為6 000像素×4 000像素,圖像格式為jpg;計算機內存為3.5 GB,主頻3.5 GHz。相機獲取圖像經過USB線傳送到計算機,最終獲取的俯視圖和側視圖各700幅。

圖1 履帶式小車戶外工作現場圖Fig.1 Work scene of autonomous crawler in outdoor1.磁條 2.背景布 3.俯視相機 4.太陽能板 5.側視相機 6.小車框架 7.控制箱 8.磁導航傳感器
實驗所用軟件為Microsoft Visual Studio 軟件、Robot Operating System下的跨平臺開源C++點云庫(Point cloud library,PCL)和三維重建開源軟件Visual SFM,性狀提取系統流程圖如圖2所示。履帶式小車在戶外以0.1 m/s的速度連續工作4 h獲取玉米植株二維圖像;將獲取的圖像導入到Visual SFM軟件進行三維重建生成點云數據;將點云數據進行下采樣、去噪、坐標校正等預處理;根據預處理后的點云數據,使用直通濾波、條件歐氏聚類、直徑擬合等算法實現自動去除背景點云、分割單株、提取莖稈和葉片等功能;基于分割出來的單株、莖稈和葉片,利用距離最值遍歷、三角面片化、條件歐氏聚類和隨機采樣一致性等算法自動測量株高、單株最小包圍盒體積、莖粗、單片葉最小包圍盒體積、葉投影寬度、葉投影長度、葉周長、葉面積、葉投影面積、葉片數和葉夾角等參數;與人工測量的株高、莖粗、葉面積作對比,來評估SfM重建玉米植株的精度和魯棒性。

圖2 性狀提取系統流程圖Fig.2 Flowchart of traits extraction system
玉米植株三維模型重建是在Visual SFM軟件中完成,重建的基本原理是利用大量重疊度較高的圖像通過SfM算法重建物體三維模型[32]。SfM算法的核心主要包括特征點提取、立體匹配、姿態估計和參數優化等[33-35]。具體重建過程如下:
(1)特征點提取與匹配:通過尺度不變局部特征描述子(Scale-invariant feature transform,SIFT)提取圖像特征,用K維空間二叉樹(kd-tree)模型計算兩幅圖像特征點之間的歐氏距離來進行特征點的立體匹配。
(2)稀疏點云重建:用集束調整算法(Bundle adjustment,BA)減少因圖像增多而累計的誤差,并估計每幅圖像的相機參數,生成稀疏點云。
(3)稠密點云重建:用多視角立體集群算法(Cluster multi-view stereo,CMVS)對圖像集進行聚簇,以減少重建過程的數據量,提高運算速度和重建精度,然后用多視角拼接算法(Patch-based multi-view stereo,PMVS)通過匹配、膨脹、過濾,在局部光度一致性和全局可見性約束下完成稠密點云的重建[36]。
將履帶式拍照小車獲取的700幅圖像導入到Visual SFM軟件進行三維重建并生成三維點云數據。圖3為玉米植株三維重建過程。

圖3 玉米植株三維重建過程Fig.3 Three-dimensional reconstruction process of maize plants
由于玉米三維重建后點云數據很大,夾雜著許多噪聲背景點云,且玉米三維點云模型與實際植株在標準三維空間方向、尺度上不一致,因此需要進行下采樣、降噪、坐標校正等預處理。
1.3.1點云數據下采樣
由于用Visual SFM軟件重建出來的原始玉米三維點云數據很稠密,為了節約算法運算時間,采用下采樣方法來精簡點云數據。下采樣是在不改變點云外輪廓的前提下大幅度減少點云的數量,從而提高算法運算效率。下采樣的方法采用三維體素柵格法[37-38],通過將輸入的點云數據創建三維體素柵格,然后在每個體素內,用體素中所有點的重心近似代替該體素內的所有其他點,這樣所有體素內的點最終都由該體素的重心所代替,從而達到精簡點云數據、減少運算時間、提高程序運行速度的目的。如圖4b所示,點的數量降至點云原圖(圖4a)的30%,但單株的外輪廓幾乎沒有變化,不影響后面性狀參數的提取。

圖4 玉米點云下采樣和去噪效果Fig.4 Effects of downsampling and denoising
1.3.2點云數據降噪
由于受RGB相機精度、實驗操作者經驗、外界環境等因素的影響,點云數據中不可避免地會出現一些噪聲點和局外點。這些噪聲點和局外點對性狀提取、特征匹配和曲面重建都有不良影響[39]。本研究采用統計方法去除噪聲點和局外點[40-41]。對于下采樣后點云數據中的每個點,計算其與k近鄰的平均距離(k近鄰的值設置為50),通過假設得到的結果服從高斯分布,具有一個平均值(μ)和標準差(σ),如果所有這些點的平均距離大于閾值Tthreshold,將被視為異常值,從數據集中刪除,計算公式如下
Tthreshold=σC
(1)
式中C——用戶所使用的常量,本研究取1
本研究也根據三維點云的顏色信息來去除噪聲點,計算所有點的超綠分量值(ExG)[42],將ExG小于0.36的點視為噪聲點,將其從三維點云數據中去除,計算公式如下
(2)
式中R、G、B——紅、綠、藍顏色分量
ExG——超綠分量值
玉米三維點云去噪效果如圖4c所示,其中點的數量降至下采樣后單株點云數量的95%,與圖4b相比,減少的點不多但都為噪聲點(圖4b中的紅色橢圓框)。
1.3.3點云數據坐標校正
為了精確地提取玉米植株性狀,需進行玉米三維點云的坐標校正。為了獲取玉米三維點云模型和真實玉米植株的比例,需進行坐標比例校正,以花盆為基準計算出坐標比例,計算公式如下
(3)
式中hreal——花盆真實高度
hreconstructed——花盆重建模型的高度
r——坐標比例
由于SfM重建的三維玉米點云模型的主方向與PCL點云庫中可視化下的三維坐標軸方向存在一定的偏差,所以需要使用旋轉與平移矩陣進行坐標轉換校正,將點云模型的質心點轉換到坐標系原點,轉換的主要步驟如下:①計算預處理后的三維點云模型質心點的三維坐標。②基于主成分分析法(Principal component analysis,PCA)[43]將點云模型進行質心化,然后計算點云模型的協方差,求解協方差矩陣的特征值和特征向量,特征向量即為玉米三維點云的主方向。③根據三維點云的質心坐標和主方向創建一個相對于空間坐標系原點的平移和旋轉矩陣。④通過創建的旋轉和平移矩陣,將原始點云轉換為以原點為中心的新點云,空間坐標軸為新點云主方向的位置。
點云坐標轉換公式如下
TA=MTTO
(4)
式中MT——平移旋轉矩陣
TO——原始點云數據
TA——轉換之后的點云數據
在可視化PCL點云數據庫下,新位置的點云模型主方向與玉米實際株高的生長方向不一致,為了便于后續點云數據的處理,需要再次使用平移旋轉矩陣使玉米植株生長方向與坐標軸z的正方向保持一致。點云坐標轉換校正流程如圖5所示。圖中,紅軸和藍軸分別表示x軸和z軸,黃色十字交點為坐標原點。

圖5 點云坐標轉換校正流程圖Fig.5 Flowchart of point cloud coordinate conversion correction
玉米植株三維點云分割主要分為分割玉米植株與花盆及以下點云、單株分割、莖稈提取以及葉片提取4部分。處理過程如圖6所示,具體點云處理步驟如下:①預處理后的原始玉米點云如圖6a所示。通過直通濾波[44],以花盆真實高度為基準將目標點云分割為玉米植株和花盆及以下點云兩部分,分割效果如圖6b所示。圖6中不同顏色代表不同類,綠色代表玉米植株點云,藍色則代表花盆及以下點云。②將分割出來的玉米植株采用條件歐氏聚類算法[45-46]分割出單株玉米植株,分割效果如圖6c所示。圖6c中不同顏色代表不同的玉米植株點云。③提取出單株玉米點云,如圖6d所示,基于隨機采樣一致性算法(Random sample consensus,RANSAC)[47]進行圓柱擬合分割出玉米莖稈,分割效果如圖6e所示。同時提取出非莖稈點云,如圖6f所示,將提取的非莖稈點云再次采用統計法去除噪聲點和局外點,去除后的效果如圖6g所示,被去除的點如圖6f中的紅色橢圓框所示,把去除噪聲點和局外點的非莖稈點云再次采用條件歐氏聚類算法分割出單片葉片,分割效果如圖6h所示,圖6h中不同的顏色代表不同的類。

圖6 玉米三維點云分割流程Fig.6 Flowchart of three-dimensional point cloud segmentation

圖7 玉米性狀參數提取流程Fig.7 Traits parameters extraction process of maizes
基于分割得到的玉米點云,計算玉米株高、單株最小包圍盒體積、莖粗、葉片數、葉周長、葉面積、單片葉最小包圍盒體積、葉投影面積、葉投影寬度、葉投影長度、葉夾角等參數,性狀參數提取流程如圖7所示,具體參數計算方法如下:
(1)玉米株高
本研究中玉米單株的最高點與盆沿平面在z軸方向的差值默認為玉米的株高。提取出單株玉米植株點云(圖7a),遍歷所有點,由于校正后單株玉米的生長方向與z軸的正方向一致,且玉米的底部平面(盆沿平面)為oxy平面,所以要找到單株玉米植株z坐標的最大值和盆沿平面的z坐標值(當前為0),其差值的絕對值即為單株玉米植株的高度(圖7b)。
(2)單株玉米植株最小包圍盒體積
提取出單株玉米植株,校正前的單株玉米植株如圖7c所示,黃色線條所組成的長方體為包圍盒,校正后的單株玉米植株已轉換至其主方向上(圖7d所示,黃色線條所組成的長方體即為最小包圍盒)。要找到校正后單株玉米點云最大的x、y、z坐標值和最小的x、y、z坐標值,得到8個頂點,這8個頂點連接起來所形成的長方體體積即為單株最小包圍盒體積。
(3)玉米莖粗
本研究中盆沿平面往上10 cm的平面內莖稈點云坐標最大值與最小值的差值默認為整株玉米的莖粗。該莖粗為本研究的定義,人工測量莖粗時也是按照這個標準進行測量。提取出莖稈點云(圖7m),首先用算法自動截取z=10 cm處平面內的點云數據(圖7n),遍歷該平面內的所有點,找到該平面內坐標的最大值和最小值,它們的差值絕對值即為莖粗(圖7o)。
(4)葉片數
提取去除噪聲點和局外點的非莖稈點云(圖7e),進行條件歐氏聚類,分割出單片葉片(圖7f,不同的顏色代表不同的類),其中聚類出的不同類的數量即為葉片數。
(5)葉周長
分割出三維玉米葉片(圖7g),采用貪婪投影三角算法[48]進行三角面片化,面片化后的三維葉片模型由若干個空間小三角面片組成,其中每個小三角面片的3條邊至少被使用一次,而最外圍的邊即葉片的邊界僅被某個小三角面片使用一次,不會被其它三角面片所共用,所以找出所有僅被使用一次的邊,計算該邊的長度并求和即為葉周長(圖7h,紅色邊界線的長度為葉周長)。
(6)玉米葉面積
分割出玉米葉片,首先采用貪婪投影三角算法進行三角面片化,面片化后的葉片模型由若干個空間三角面片組成,然后通過海倫公式計算單個空間三角面片的面積,最后通過面積求和公式計算單個葉片的面積,計算公式如下
(5)
(6)
式中pi——面片化三角形周長的一半
ai、bi、ci——面片化三角形各邊邊長
n——總面片數i——面片索引序號
Si——單個空間三角形面片的面積
S3D——單片葉片總面積
(7)玉米單片葉最小包圍盒體積
分割出玉米單片葉片后,將葉片點云轉換至主方向位置,轉換前的玉米單片葉片如圖7i所示,黃色線條所組成的長方體為包圍盒,轉換后的玉米單片葉片如圖7j所示,黃色線條所組成的長方體即為最小包圍盒。找到轉換后單片葉點云最大的x、y、z坐標值和最小的x、y、z坐標值,得到8個頂點,這8個頂點連接起來所形成的長方體體積即為單片葉最小包圍盒體積。
(8)玉米葉投影面積
將分割出的玉米葉片向oxy平面投影(圖7p),生成對應的投影葉片點云(圖7q),首先采用貪婪投影三角算法對投影葉片點云進行三角面片化,面片化后的投影葉片模型由若干個平面三角面片組成,然后通過海倫公式計算單個平面三角面片的面積,最后通過面積求和公式計算單片葉片的葉投影面積,計算公式如下
(7)
(8)
式中pj——面片化三角形周長的一半
aj、bj、cj——面片化三角形各邊長
m——總面片數j——面片索引序號
Sj——單個平面三角形面片投影面積
S2D——單片葉片總投影面積
(9)葉投影寬度
將分割出的玉米葉片向oxy平面投影,生成對應的投影葉片點云,把生成的投影葉片點云轉換至主方向位置(圖7r),計算寬度方向坐標最大值和最小值,其差值的絕對值默認為本研究中的葉投影寬度(圖7r)。
(10)葉投影長度
將分割出的玉米葉片向oxy平面投影,生成對應的投影葉片點云,把生成的投影葉片點云轉換至主方向位置,計算長度方向坐標最大值和最小值,其差值的絕對值默認為本研究中的葉投影長度(圖7r)。
(11)葉夾角
分割出單片玉米葉片后,由于整片玉米葉片彎曲幅度比較大且每片玉米葉片彎曲的幅度不一致,將其用隨機采樣一致性算法進行平面擬合,求取所得到的平面與莖稈的一個方向向量的夾角,該方法所求取的葉夾角可信度比較低,因此本研究統一截取葉片點云的一部分進行葉夾角求取。截取的規則為:以靠近莖稈的一端為起始端,在向葉片延伸的方向上截取某個閾值內的葉片點云(圖7k所示,綠色為完整的單片葉片點云,紅色為截取的點云,藍色為莖稈所在位置),將截取的葉片點云進行平面擬合,擬合出的平面與莖稈的一個方向向量夾角默認為葉夾角(如圖7l所示,紅色為截取的點云,黃色平面為截取的點云所擬合出來的平面,藍色為莖稈所在位置)。
將自動測量的玉米植株株高、莖粗、葉面積與人工測量值進行對比來評價本文方法的精度。測量結果精度用平均絕對百分比誤差(Mean absolute percentage error,MAPE)、均方根誤差(Root mean square error,RMSE)和決定系數R2(Determination coefficient)來評估。
本實驗對抽穗期的80株玉米樣本用SfM進行三維重建,自動分割玉米的單株、莖稈和葉片點云,并進行測量分析。將自動測量的株高、莖粗、葉面積和葉投影面積與人工測量值進行對比,結果如圖8所示。

圖8 玉米植株性狀的人工測量值與算法測量值的對比圖Fig.8 Comparison plots of algorithm and manual measurements of maize plant traits
由圖8a可看出,R2=0.970,MAPE為3.163%,RMSE為3.557 cm,算法測量株高的準確率為96.837%。由圖8b可知,R2=0.842,MAPE為4.760%,RMSE為1.540 mm,算法測量的莖粗準確性為95.240%。由圖8c可看出,R2=0.901, MAPE為19.102%, RMSE為48.163 cm2,算法測量的葉面積準確率為80.898%。由圖8d可看出,R2=0.728,MAPE 為69.321%, RMSE 為95.493 cm2。相比圖8c,投影葉面積測量值與葉面積測量值相差較大,說明葉面積測量方法比單個方向(投影至oxy平面)投影葉面積測量方法更接近實際玉米葉面積。圖8結果表明,本文方法具有較高的準確性,算法測量值與人工測量值具有較好的一致性。
單株最小包圍盒體積、葉片數、單片葉最小包圍盒體積、葉周長、葉投影寬度、葉投影長度和葉夾角7個性狀參數的測量結果如圖9所示。在單株性狀參數上,晉單73品種玉米的單株最小包圍盒體積和葉片數總體上高于大豐30和鄭單958品種,而在單片葉性狀上,3個品種的5個性狀間(單片葉最小包圍盒體積、葉周長、葉投影寬度、葉投影長度和葉夾角)沒有非常顯著的差異。
將80株玉米植株分為低地上部生物量和高地上部生物量兩類,通過SPSS軟件將單株玉米植株的株高、最小包圍盒體積、莖粗和葉片數4個性狀做單因素方差分析,顯著性差異P值分別為0.000 3、0.000 4、0.317 0和0.241 5。當P≤0.01時,說明因素對試驗指標的影響非常顯著;當0.01
0.05時,說明因素對試驗指標的影響不顯著。單株玉米植株的株高、最小包圍盒體積、莖粗和葉片數區分低地上部生物量玉米植株和高地上部生物量玉米植株的顯著性統計結果如圖10所示。由此可知,單株玉米植株的株高、最小包圍盒體積可以顯著區分低地上部生物量玉米植株和高地上部生物量玉米植株,而莖粗和葉片數不能區分低地上部生物量玉米植株和高地上部生物量玉米植株。
針對傳統的玉米植株性狀測量方法存在主觀性強、勞動強度大、有損傷等問題,以80株抽穗期玉米植株為研究對象,研究了基于運動恢復結構的玉米植株三維重建方法和性狀提取方法。結果表明,與人工測量值相比,株高、莖粗、葉面積的平均絕對百分比誤差(MAPE)分別為3.163%、4.760%、19.102%,均方根誤差(RMSE)分別為3.557 cm、1.540 mm、48.163 cm2,決定系數(R2)分別為0.970、0.842、0.901。同時,實驗驗證了株高和單株最小包圍盒體積可以顯著區分低地上部生物量玉米植株和高地上部生物量玉米植株。

圖9 玉米葉片性狀參數測量結果Fig.9 Measurement results of trait parameters in maize leaves

圖10 單株性狀參數顯著性差異分析Fig.10 Analysis of significant differences in trait parameters of single plant