李 哲,宋青峰,朱新廣,胡 勇*,鞏彩蘭,卜弘毅
(1中國科學院紅外探測與成像技術重點實驗室,中國科學院上海技術物理研究所,上海 200083;2中國科學院大學,北京 100049;3中國科學院分子植物科學卓越創新中心,上海 200032)
玉米是全世界的主要糧食作物之一[1],提高玉米產量能夠有效解決糧食安全問題。在當今大數據時代,植物表型組學在育種中的作用開始得到重視[2],該學科通過研究植物個體間表型參數差異,能夠有針對性地對農作物進行改良育種,提高作物產量[3]。植物表型組學的發展瓶頸在于如何獲取大量植株表型參數。
傳統的表型參數獲取主要依靠人工使用直尺和量角器等工具測量,存在效率低、主觀誤差大、有破壞性等缺點。為解決這些缺點,三維表型參數提取方法[4-6]成為研究熱點。目前三維信息提取類研究主要通過深度相機掃描獲取三維點云,如Kinect法[7]。此方法能夠直接獲得高精度的植株點云,操作簡單,但存在所需設備昂貴、使用受天氣影響、數據量龐大等問題。隨著相機硬件設備和圖像算法的發展,基于機器視覺重建三維結構方法[8]因設備廉價、數據提取速度快、魯棒性好等優點,得到廣泛應用。但玉米葉片具有薄、存在翻折、葉片邊緣為不規則波浪形、互相遮擋等特點,目前基于機器視覺提取玉米三維表型特征的方法存在準確率低、耗時長、需要一定手工操作、不能完全自動化等問題。
針對上述問題,本研究擬提出一種高精度、高通量和自動化的玉米生長動態定量化測量方法,以期實現玉米植株三維表型參數以及生長動態自動化測量。
田間試驗在中國科學院上海生命科學研究院植物生理生態研究所實驗基地(30°94'N,121°13'E)進行。玉米選用自交系A619和W64A,這兩個玉米自交系株型結構差別主要在于W64A的葉片數比A619多、葉片長度比A619短,兩者的莖稈高度和葉片寬度無明顯差異。每個玉米自交系種植15盆,10盆用于測量,5盆用于備份。盆栽種植采用白色20 L塑料桶,每桶播種2粒已催芽的種子,間距7 cm,播種深度3 cm,長至三葉期后每桶僅留下1株長勢均勻的玉米。植株放置于戶外生長。播種后30 d開始測量,每隔7 d測一次,測量持續到玉米抽穗,共測量4次。
測量時將植株置于室內,采用瞬時成像設備[9]拍照。瞬時成像設備由64臺18—55 mm鏡頭的佳能EOS 1300D型相機、拍攝支架和分布式數據采集控制器組成(圖1),每個支架8臺相機,共8個支架圍繞植株拍攝。瞬時拍攝系統能在亞秒級時間輸出64個不同角度的拍攝圖像,避免拍攝同一植株過程中因光線變化或植株抖動等因素導致的三維重建誤差。相鄰兩個角度的拍攝圖像重疊率達75%以上,保證三維重建獲取點云的準確率。

圖1 瞬時成像設備結構圖Fig.1 Instantaneous imaging equipment structure
采用張正友[10]的方法標定相機參數,內參標定使用棋盤格,外參標定使用標定架。標定架如圖2所示,標定架一共包含160個靶標,其中149個黑色靶標用于圖像拼合,11個白色靶標用于確定坐標與尺寸。標定每臺相機參數后,移動世界坐標系,使得XY軸方向如圖2c所示,Z軸沿標定架方向,原點位于地面上方25 cm處,即玉米植株基部。以靶標之間距離為標準調整世界坐標系的比例,使得世界坐標系下的長度與真實世界長度相同。

圖2 標定架結構圖Fig.2 Calibration frame structure
測量當天,先拍攝一組背景圖像,再依次把每棵植株放到標定的世界坐標系原點位置拍攝,每張拍攝圖像去除對應背景圖像可以得到無背景的植株圖像,如圖3所示。
拍攝得到64個角度的無背景植株圖像后,基于機器視覺獲取玉米植株三維點云。首先利用尺度不變特征變換從每幅圖像中提取特征點;然后使用快速近似最近鄰搜索庫實現兩兩圖像間描述子向量匹配,并選擇匹配點最多的兩幅圖計算初始點云;之后每次選取其余圖像的特征點中與目前已知點云匹配數量最高的圖像加入,利用隨機抽樣一致算法根據匹配點求新增圖像的投影矩陣,并將新計算的三維點加入已知點云;不斷重復上一步過程,直至所有圖像都被添加,再利用光束平差法進行優化得到稀疏點云;最后通過多視覺聚簇法將稀疏點聚類到不同的影像集,每個影像集分別基于面片的三維多視角立體視覺算法通過法線和位置約束對稀疏點云形成的種子面片進行擴散,利用灰度一致性和幾何一致性約束對面片過濾,最終得到稠密三維點云。
去除背景可以減少背景噪聲點云和錯誤匹配點云,提高點云精度,但點云重建結果仍存在少部分背景噪聲,可以利用世界坐標系位置特殊性再次去除背景噪聲。玉米植株根部位于世界坐標系原點而且世界坐標系尺度與真實世界一致,設置植株大致的長寬高范圍,即植株點云的坐標值范圍,可以獲取無背景噪聲植株點云(圖3 d)。

圖3 去除背景以及三維重建結果Fig.3 Background removal and 3D reconstruction results
基于拉普拉斯的骨架提取算法[11-12]得到玉米點云骨架的主要步驟為:(1)構造單環鄰域;(2)幾何收縮,獲取零體積的網格或點云;(3)拓撲細化,得到一維曲線;(4)中心化處理,得到骨架。該算法魯棒性強,但提取的玉米骨架誤差大,不符合玉米植株特性。本研究通過葉尖和葉基部重定位骨架提高參數精度。
葉尖提取優化:葉尖點缺失導致計算葉長值小于實測值,通過葉尖提取優化提高葉長精度。優化步驟如下:1)根據骨架點之間鄰接關系劃分葉片骨架和莖稈骨架;2)依據式(1)得到第i個葉片點云LPi;3)在LPi中找與距離最大的100個點作為葉尖備選點Ti。是LSi最后一個點,即第i個葉片基部點。是LSi第一個點,即第i個葉片葉尖點。(4)對于Ti中第j個點(j=1,2,3,...,100),它與30個鄰近點nk(k=1,2,3,...,30)的向量和為,在滿足式(2)的所有點中選取與距離最遠的點,作為葉尖點,添加到第i個葉片骨架中。

式中LSi是第i個葉片骨架點集合,P是濾波點云集合。

葉尖提取優化即在原始骨架上新增一個骨架點作為葉尖點,圖4a中黑色的點為新增葉尖點,紅色點為原始骨架。新增葉尖點后,骨架長度隨之增加,而且從圖4b中可以看出新增葉尖點全部精準定位在植株真實葉尖位置,優化后葉片骨架長度彌補了原始骨架長度缺失,更加接近葉長的真實值。

圖4 葉尖提取優化結果示意圖Fig.4 Optimization results of leaf tip extraction
莖稈骨架優化:經觀察發現,玉米莖稈骨架在葉片基部附近向葉片方向傾斜,針對這個問題,利用玉米莖稈基本呈直線的先驗知識,優化莖稈骨架。優化步驟如下:(1)計算莖稈骨架點集SS與均值點,ss∈SS的誤差和矩陣;(2)使用奇異值分解[13]計算e的特征向量,e最大特征值對應的特征向量就是莖稈的方向向量;(3)莖稈擬合直線為,將每個莖稈骨架點z值帶入l重新計算位置。
葉基部重定位:葉片基部與莖稈相連,但莖稈骨架優化后,葉基部可能不在莖稈上,需要重新定位。優化過程偽碼如圖5所示。

圖5 葉基部重定位偽代碼Fig.5 Pseudocode for leaf base relocation
莖稈骨架優化前后,莖稈骨架點位置如圖6a和6b所示,優化結果符合玉米植株實際情況。葉基部重定位前后,葉基部點位置如圖6c和6d所示,優化后的葉基部點全部落在莖稈擬合直線上,還原了葉基部點的真實位置,葉基部點精確定位是求取高精度葉夾角和葉長等參數的基礎。

圖6 葉基部高度重定位Fig.6 Relocation of leaf base height
株高、葉基部高度和葉長提取:將植株基部到植株最上部展開葉的最高點的高度差作為株高,由于植株基部在Z=0平面上,所以株高就是所有點云坐標z值的最大值;葉基部高度是植株基部到葉片著生點位置的高度差,葉基部重定位后已知每個葉片的著生點坐標,又由于植株基部在Z=0平面上,葉基部高度即為骨架著生點坐標z值;葉長為葉片的相鄰骨架點之間連接線段長度的總和。
葉最大寬度提取:做一個過第i個葉片骨架中點,垂直于中點切線方向的垂面。由于點云P稀疏,垂面Si與點云的交點較少甚至可能沒有,所以應該將垂面Si左右移動0.5 cm,對夾在中間的點集Wi做直線擬合,如圖7a所示;將Wi所有點中x值最小和最大的兩點和值帶入擬合直線得到和,第i個葉片的葉最大寬度如式(3)所示。

葉夾角提取:規定莖稈與葉片1∕4長度位置的中心點為葉夾角,如圖7b所示。其中點o是葉片骨架著生點,點m是莖稈擬合直線上一點,點n是葉片1∕4長度位置截線的中點,截線獲取方式與葉最大寬度獲取方式相同。葉夾角按公式(4)計算。

最小包圍盒體積提取:對植株點云P進行主成分分析,并將植株點云轉到主成分方向上構成新的點云集合P'。最小包圍盒是求取能夠完整包裹物體且棱長平行于坐標軸的最小長方體,如圖7c所示,將點云在XYZ軸的范圍相乘即可得到最小包圍盒體積,如式(5)所示。


圖7 植株表型參數提取示意圖Fig.7 Plant phenotypic parameters extraction
如表1所示,兩個玉米自交系在骨架優化后葉長和葉基部高度精度都有明顯提升。這是因為葉尖提取優化和葉基部重定位分別精確定位了每個葉片的葉尖和葉著生點,所以自交系A619和W64A的葉長均方根誤差下降了50%左右。葉基部高度精度提升是由于葉基部重定位,得到了葉基部點精確位置。同時可以看到,骨架優化對葉最大寬度參數準確率也有少量提升,這是因為人工測量的葉最大寬度是葉片一半長度處的展開寬度,與葉長相關,葉最大寬度精度隨著葉長精度提高也稍有提升。

表1 骨架提取算法改進前后表型參數誤差對比Table 1 Comparison of phenotypic parameter errors before and after improvement of skeleton extraction algorithm
參數提取結果中葉基部高度略大于真實測量值,而葉長小于真實測量值,這是由于人工測量葉基部高度和葉片長度是從葉片背面測量,而參數提取算法是從莖稈和葉片分離較明顯的位置測量得到的,即葉基部包裹莖稈那一部分是誤差的主要來源。
A619和W64A兩個玉米自交系葉長、葉最大寬度和葉基部高度的測量值和計算值的比較如圖8所示。葉長、葉最大寬度和葉基部高度的R2均在0.9以上,表明即使是不同的玉米自交系,此方法也能夠較準確地反映葉長、葉最大寬度寬和基部高度屬性,展現了較好的魯棒性。

圖8 A619和W64A的葉長、葉最大寬度和葉基部高度的測量值和計算值Fig.8 Measured and calculated values for leaf length,maximum leaf width and leaf base height of A619 and W64A
如圖9所示,出苗后21—42 d,自交系A619株高從42 cm生長到138 cm,自交系W64A株高從40 cm生長到134 cm。A619和W64A在出苗后28—35 d快速增高,35 d后株高基本不變。

圖9 A619和W64A株高動態變化Fig.9 Dynamic changes of plant height of A619 and W64A
圖10展示了自交系A619和W64A葉長、葉最大寬度和葉基部高度的變化曲線。由于玉米底部葉片會在生長過程中逐漸衰老脫落,本研究選取第4—7片葉子(L4—L7)記錄動態變化曲線,葉片按照葉基部高度從下到上排序,序號越小葉齡越大。從誤差棒可以看出,W64A的動態范圍較大,每棵植株間差異大。從變化曲線可以看出,A619在出苗后28 d,第4、5、6片葉子完全展開,之后隨著植株生長,葉片發黃萎縮,葉片長度和葉最大寬度略微減小;第7片葉子在整個測量時間一直在生長狀態。W64A生長狀態與A619基本一致,不同在于第5和第6片葉子直到出苗后第35天才完全展開。

圖10 A619和W64A的葉長、葉最大寬度和葉基部高度變化曲線Fig.10 Variation curves of leaf length,maximum leaf width and leaf base height of A619 and W64A
植株點云骨架在器官分割以及參數求解中有較好的效果,近兩年越來越多的研究人員將植株骨架用于三維表型參數提取。2018年,溫維亮等[14]利用接觸式三維數字化儀獲取玉米各器官的植株骨架結構,并制定了參數獲取的規范流程。該方法需要人工使用探測筆沿著玉米每個器官劃過,通過定位跟蹤得到骨架,參數獲取方法耗時耗力,而且葉片易受力位移造成誤差。本研究只需使用瞬時成像設備拍攝圖像,經由算法在10 min內獲取玉米骨架并輸出參數,避免人為誤差,實現快速和全自動化測量。2019年,蘇偉等[15]提出基于地面激光掃描獲取大田玉米表型參數的方法,該方法提取表型參數方法與本研究相同,都是通過點云得到三維骨架,然后利用骨架信息獲取三維表型參數。但該方法激光掃描獲取點云掃描時間長、點云數據量龐大,而且需要手動完成點云配準和裁剪,不能實現完全自動化;而本研究采用的基于機器視覺的三維重建只需要多角度圖像就可以實現自動、快速的點云獲取。2020年,Wu等[16]提出了一種低成本的表型分析平臺,該方法通過計算機視覺得到植株三維模型,通過骨架提取算法得到骨架,并利用骨架提取表型參數,該方法株高、葉寬和葉面積的均方根誤差分別為2.96、1.27和75.03。本研究同樣采用機器視覺和骨架提取算法,但是通過改進植株骨架,提升了獲取表型參數精度。本研究中玉米葉長、葉最大寬度、葉基部高度的均方根誤差分別為4.64、0.47和3.99,對于相同的參數葉寬,均方根誤差相對于該方法少50%以上。
本研究采用瞬時成像設備,實現64個角度圍繞植株瞬時成像,且相鄰相機圖像重疊度達75%以上,保證了重建精度;利用基于拉普拉斯的骨架提取算法獲取植株點云骨架,并對葉尖和葉基部骨架進行改進,準確分割葉片,提高了參數提取精度;從圖像直接得到三維表型參數,實現了完全自動化。
本研究提出了一種高精度、高通量、自動化的玉米表型參數提取算法,該算法主要改進了植株骨架葉尖和葉基部的識別過程。利用改進的算法獲得的玉米葉長、葉最大寬度和葉基部高度的平均絕對百分比誤差分別下降5.89%、0.04%和0.86%。基于本研究方法提取的表型參數決定系數均大于0.9,證明算法能夠準確地反映葉長、寬和基部高度屬性,并用于全自動化參數提取,為高通量玉米表型分析提供了有效的方法。