朱冰琳 李 敏 劉扶桑 賈奧博 毛 秀 郭 焱
(中國農(nóng)業(yè)大學(xué)土地科學(xué)與技術(shù)學(xué)院, 北京 100193)
高效選育高產(chǎn)的作物新品種是解決世界糧食安全問題的重要途徑[1]。影響現(xiàn)代育種的關(guān)鍵因素之一是能否快速、精確地對大量育種品系進行全面的表型評價,從中篩選出具有期望表型的品系[2-3]。目前,基因標記輔助選擇、基因組編輯[1]、高通量基因測序等技術(shù)已得到廣泛應(yīng)用,利用這些技術(shù)可以高效地挖掘調(diào)控作物重要性狀的基因信息[4]。然而,作物表型的研究技術(shù)和方法卻無法與高通量基因測序技術(shù)的發(fā)展相匹配[5-6]。
近年來興起的高通量獲取作物表型的方法或“表型組學(xué)”,被認為可以顯著提高植物表型信息的獲取效率[7]。目前開發(fā)的表型平臺大多應(yīng)用于溫室等特定環(huán)境中[8-10]。田間表型平臺也逐漸開發(fā)出來,但往往體積龐大、價格昂貴[11],且存在較大的應(yīng)用局限。如星載表型平臺較難應(yīng)用于面積較小的育種實驗小區(qū),不適合于高頻率的動態(tài)監(jiān)測[7,12];機載平臺的花費較大,需要較多的人員維護,且在低空獲取高分辨率圖像時對冠層的擾動較大[13];地面表型平臺能夠同時獲取多光譜、高光譜及熱紅外等圖像信息,但通常為特定株高范圍的作物設(shè)計,其行走受限于田間地面狀態(tài),靈活性和便攜性較差[7]。
目前,基于RGB相機獲取植株多視角圖像序列后,利用運動恢復(fù)結(jié)構(gòu) (Structure from motion, SFM)算法可重建植株三維結(jié)構(gòu)。該方法已在溫室中廣泛應(yīng)用,如監(jiān)測溫室內(nèi)小麥、茄子、黃瓜、辣椒等生長過程[14-15],對大田玉米也有初步的嘗試[16],但重建范圍較小。無人機(Unmanned arial vehicle, UAV)平臺搭載RGB相機,以數(shù)十米或更高的飛行高度獲取冠層圖像序列,可以實現(xiàn)大田尺度的冠層點云重建[17],并可進一步計算株高[18]、植被指數(shù)(Normalized difference vegetation index, NDVI)[19-20]、葉面積指數(shù)(Leaf area index, LAI)[21]、植被覆蓋度[17]等群體冠層參數(shù)。然而,實施精準農(nóng)業(yè)需要有作物冠層結(jié)構(gòu)精確信息的支持[22],基于理想株型選育作物新品種[23]也需要精確的冠層結(jié)構(gòu)信息[24]。隨著無人機的微型化發(fā)展,其對冠層的擾動變小,可以通過超低空飛行獲取更清晰的冠層圖像,使快速獲取高精度的冠層結(jié)構(gòu)信息成為可能。
通過點云數(shù)據(jù)網(wǎng)格化可建立三維網(wǎng)格的冠層結(jié)構(gòu)模型。文獻[25]利用Geomagic Studio軟件實現(xiàn)了蘋果樹葉片的網(wǎng)格重建;文獻[26]對黃瓜、西瓜、甜瓜等植株的葉片利用Delaunay三角剖分的方法生成初始網(wǎng)格曲面,再對網(wǎng)格曲面進行優(yōu)化處理,實現(xiàn)了冠層結(jié)構(gòu)建模;文獻[27]采用改進的Crust方法實現(xiàn)了樹干的網(wǎng)格重建。這些研究大多在溫室內(nèi)進行,且葉片形態(tài)較為簡單,受遮擋較少;在室外選擇的研究對象(如樹干),其形態(tài)也不易受到風(fēng)的影響,重建效果較好。目前,基于點云數(shù)據(jù)對大田作物群體冠層結(jié)構(gòu)建模的研究較少。
本文利用超微小型無人機超低空航拍獲取大田不同生育期植株個體及群體的多視角圖像序列,基于偽極點-Crust算法構(gòu)建其冠層結(jié)構(gòu)模型,并基于大田實測數(shù)據(jù)對重建冠層結(jié)構(gòu)模型的精度進行評估,以期為高效、高精度構(gòu)建大田作物冠層結(jié)構(gòu)模型、獲取特定的作物表型參數(shù)提供新途徑。
玉米大田實驗于2019年在中國農(nóng)業(yè)大學(xué)吉林梨樹實驗站(43°16′N,124°26′E)進行。供試品種為先玉335,種植株距為20 cm,行距為50 cm。共3個小區(qū),小區(qū)長、寬為24、9 m(圖1a)。種前施加基肥,用量為純N 80 kg/hm2、P2O5120 kg/hm2、K2O 100 kg/hm2,采用鏵式犁和旋耕機將肥料旋進土壤耕層中。在玉米第6葉及第11葉展開時追施40 kg/hm2純N。大田出苗日為5月19日。
出苗后43 d(苗期),選取一個先玉335小區(qū)的玉米群體為目標,隨機選擇相鄰的4個條帶,并從中選取18棵植株(圖1b),分別對每株上、中、下部位各葉片用標簽紙標記。為確定幾何換算系數(shù),在目標群體中心的地面上放置邊長為40 cm、啞光材質(zhì)的正方形校正板。在多云且風(fēng)速小的天氣條件下,采用搭載哈蘇相機的Mavic 2 Pro型(大疆創(chuàng)新科技有限公司,中國)超微小型無人機對目標群體進行航拍。航拍時采用圓周型航線飛行2圈,飛行高度及半徑均設(shè)置為5 m。相機鏡頭與水平面的夾角為45°,飛行速度為0.8 km/h,相鄰圖像的拍攝間隔為2 s,圖像重疊度為80%。每次航拍獲取160~180幅圖像(表1)。在出苗后112 d(成熟期)進行了同樣方式的航拍,由于此時植株中下部相互遮蔽嚴重,因而選取一株、一行內(nèi)相鄰2、3、4株玉米為航拍目標(圖1b),航拍前去除目標植株周邊的植株,并標記葉片、放置校正板。每次航拍后,立即測量目標植株標記葉片的葉長、葉最大寬及株高(植株基部與其最高點間的垂直距離)。

表1 無人機平臺、飛行參數(shù)及圖像信息Tab.1 UAV platform, flight parameters and image information
使用Metashape Professional (Agisoft, RU)軟件對航拍獲取的多角度圖像序列(圖2a)進行冠層點云重建。獲得點云后(圖2b、2c),手動將選定的目標植株分割出來,利用支持向量機(Support vector machine, SVM)分類器根據(jù)顏色信息進行點云分類(圖2d),刪除后得到僅包含植被顏色信息的點云(圖2e)。在成熟期時,由于地面有部分與植被顏色相近的落葉,采用隨機抽樣最大似然估計算法(Maximum likelihood estimation sample consensus, MLESAC)[28]對地面落葉點云進行平面擬合(圖2f),去除后得到植株點云(圖2g)。
由于點云數(shù)據(jù)來源于大田航拍圖像,因而不可避免地帶有噪聲點。為消除噪聲點對冠層結(jié)構(gòu)構(gòu)建的干擾,本文采用基于空間單元格搜索k-近鄰點的方法剔除噪聲點。設(shè)P={p1,p2, …,pn}是重建冠層上所有點的集合,P中與pi歐氏距離最小的k個點稱為pi的k-近鄰,記為N(pi)。該算法的核心思想是對點云數(shù)據(jù)的最大包圍盒進行一定步長(本文設(shè)定為0.1 m)的空間劃分(圖2h),通過計算一個點p所在的子包圍盒及其在空間上相鄰的26個子包圍盒(圖2i)中所有點的歐氏距離,由小到大進行排序,選取前k個距離,即為該點與k個近鄰點間歐氏距離的集合。通過k-近鄰的歐氏距離平均值D(pi)來衡量此點的局部鄰域關(guān)系,根據(jù)點集P的整體鄰域關(guān)系確定距離閾值。如果D(pi)在設(shè)定閾值的范圍外,則認為是噪聲點[25](圖2j),刪除后得到最終的植株點云(圖2k)。
Crust算法為基于點云重建三維網(wǎng)格結(jié)構(gòu)的經(jīng)典算法之一。采用該算法時,先對點云中所有點進行一次三維Delaunay剖分,得到樣本點集的Voronoi頂點集。在此基礎(chǔ)上計算各點的正負極點,將得到的極點集與已有點云點集的并集再次進行三維Delaunay剖分[29], 采用傳統(tǒng)Crust算法計算量很大,效率較低。為了提高剖分的效率,本文在構(gòu)建冠層結(jié)構(gòu)模型時,在去噪后的點云中添加少量的偽極點,利用偽極點代替正負極點構(gòu)建原始點集,從而較大程度地減少剖分的計算量,將其稱之為偽極點-Crust算法[27]。
執(zhí)行偽極點-Crust算法時,首先確定樣本點云的最大包圍盒,并以包圍盒X、Y、Z方向上距離最大值的M倍(本文將M設(shè)置為0.2)作為步長向外拓展得到更大的四面體包絡(luò)面[27]。分別在包絡(luò)面上對各邊N等分處(本文N取4)均勻插值形成6N2個(即96個)偽極點集合(圖3a,以1行4株玉米點云群體偽極點的建立為例)。然后對偽極點與已有點云點集的并集進行三維Delaunay剖分,得到初始四面體集合To。該集合中存在較多非冠層結(jié)構(gòu)的三角面片,需通過以下步驟進行篩選。
(1)冠層結(jié)構(gòu)四面體的篩選遍歷集合To,判斷各四面體的頂點中是否含有偽極點。
如果包含偽極點,則將其標記為“非冠層結(jié)構(gòu)”。然后,獲取To中所有的三角面片并剔除重復(fù)的三角面片得到集合T′o。遍歷T′o,判斷每個三角面片與To中各初始四面體的映射關(guān)系,標記不在集合邊界的所有三角面片。為評估非邊界三角面片所在的2個相鄰四面體間的相似度,采用了由這2個相鄰四面體各自外接球半徑ri1、ri2之間夾角θ的余弦值表征的交集因子F(Intersection factor) (圖3b)
(1)
其中
F(i)∈[-1,1] (i=1, 2, …,j)
式中j——所有非邊界三角面片所在的四面體的個數(shù)
ci1、ci2——2個相鄰四面體外接球的球心坐標
設(shè)置F的閾值Ft為0.95,F(xiàn)大于Ft表明θ小,即2個外接球重合的部分很大(圖3b);F小于-Ft表明θ大,即2個外接球重合的部分很小(圖3b)。遍歷具有相同非邊界三角面片的相鄰四面體,通過標記為“非冠層結(jié)構(gòu)”屬性的四面體及F與Ft的大小關(guān)系判斷另一個四面體的屬性。如果遍歷一次后未能標記所有四面體的屬性,則將閾值逐步降低,直至每個四面體的屬性標記完畢。將屬性為“非冠層結(jié)構(gòu)”的四面體剔除,得到屬于冠層結(jié)構(gòu)的四面體集合Tp。
(2)冠層三角面片的提取
從得到的冠層四面體集合Tp中取出所有三角面片,計算各三角面片3個頂點Z坐標之和,以其和最大的三角面片為起始面片,確定其法向量的方向,然后計算其鄰接三角面片的法向量,將法向量變化最小的三角面片作為其相鄰的三角面片,直至遍歷所有三角面片。
由上述篩選過程得到的三角面片集構(gòu)建初始冠層結(jié)構(gòu)模型(圖3c)。由于在大田環(huán)境下獲取的點云密度不均一,利用該方法構(gòu)建的初始冠層結(jié)構(gòu)模型中存在大量的邊長畸長的異常面片(圖3c)。為此,計算了所構(gòu)建模型中所有三角面片最大邊的長度,統(tǒng)計確定邊長的閾值后,剔除最大邊長度大于閾值的三角面片(圖3d)。為消除所得到的器官曲面不光滑問題,將面片的頂點沿其法向量方向根據(jù)平均曲率進行移動平滑[30],得到最終的冠層結(jié)構(gòu)模型(圖3e)。
根據(jù)校正板點云的計算邊長及實際邊長確定長度換算系數(shù)后,基于所構(gòu)建模型計算目標植株的株高;利用Geomagic Studio (Geomagic,美國)手動提取標記葉片,計算其葉長及最大葉寬;對組成葉片的每個小三角面片面積進行累加,獲得標記葉片的葉面積。基于大田標記葉的葉長和葉最大寬的測量值計算參比葉面積為
Am=αLmWm
(2)
式中Lm——標記葉的葉長
Wm——最大葉寬
α——經(jīng)驗系數(shù),取0.75
通過比較實際測量的株高、葉長、葉最大寬、參比葉面積與由冠層結(jié)構(gòu)模型得到的計算值,以決定系數(shù)(R2)、均方根誤差(Root mean square error,RMSE)、相對均方根誤差(Root mean square relative error,rRMSE)、平均誤差(Mean error,ME)來評價基于UAV航拍圖像重建的冠層結(jié)構(gòu)模型的精度。
在大田完成一個目標冠層的航拍需要約5 min。利用Metashape Professional軟件對各圖像序列處理后,得到目標群體的冠層點云(圖4a、4d,成熟期以1行4株群體為例,下同)。經(jīng)SVM分類器對點云分類并基于MLESAC算法對地面落葉點云進行平面擬合后,去除非植被點云,得到目標群體的植被點云(圖4b、4e)。對所得到的植被點云采用空間單元格搜索k-近鄰點的方法確定噪聲點。去噪后,不同生育期的冠層點云均能很好地保留冠層信息(圖4c、4f)。
基于偽極點-Crust法進行了冠層結(jié)構(gòu)模型的初步構(gòu)建,所得到的冠層模型植株間及葉片間存在大量異常面片(圖5a、 5d,黑色矩形框內(nèi)為指定區(qū)域放大后的葉片形態(tài))。通過對每個三角面片最大邊長進行篩選,剔除異常面片,得到相對精確的冠層結(jié)構(gòu)模型(圖5b、5e)。最后對冠層表面進行平滑處理,獲得最終的冠層結(jié)構(gòu)模型(圖5c、 5f)。冠層面片顏色由構(gòu)成每個三角面片頂點的顏色插值得到。
為評估基于無人機航拍圖像重建的冠層結(jié)構(gòu)模型的精度,比較了株高、葉長、葉最大寬、葉面積的模型計算值與大田實測值(圖6,圖中n為樣本數(shù))。結(jié)果表明,苗期冠層結(jié)構(gòu)模型的計算株高、葉長、葉最大寬、葉面積與大田實測值的R2均在0.90以上;RMSE、rRMSE和ME的值均很小(表2);成熟期重建的冠層結(jié)構(gòu)模型的株高、葉長、葉最大寬計算值與大田實測值間的R2均在0.90以上;RMSE、rRMSE和ME的值均很小(表2);計算葉面積與參比葉面積的一致性稍差,相應(yīng)評估參數(shù)R2、RMSE、rRMSE、ME分別為0.76、87.7 cm2、0.187、-27.1 cm2(表2)。

表2 冠層結(jié)構(gòu)模型精度評價參數(shù)Tab.2 Parameters for evaluating accuracy of canopy structure models
農(nóng)田冠層輻射利用效率是影響作物產(chǎn)量的重要因素。將冠層結(jié)構(gòu)模型與冠層光分布模型、光合模型相耦合,可實現(xiàn)作物新株型的精確評估及大田栽培模式的優(yōu)化,從而有利于作物高產(chǎn)。進行此研究的前提是獲取精確的大田冠層結(jié)構(gòu)信息。采用立體視角或多視角圖像方法[15,31]、三維數(shù)字化方法[32]等雖然可以較精確地獲取大田作物冠層結(jié)構(gòu)信息,但多費時、費力,通常只能獲取較小群體的結(jié)構(gòu)信息,使得所構(gòu)建的冠層結(jié)構(gòu)模型不具有很好的代表性。故迫切需要開發(fā)高通量、高精度獲取大田作物冠層結(jié)構(gòu)的方法。
采用無人機平臺獲取大田作物信息,具有高通量的優(yōu)勢。由于以往所采用的無人機在超低空飛行時,會明顯擾動冠層,從而難以獲得高精度的航拍圖像序列。本研究采用超微小型無人機,其攜帶的RGB相機分辨率高,且在距冠層數(shù)米高度處獲取冠層圖像時對冠層的擾動較小。獲取半徑為5 m的圓形區(qū)域內(nèi)數(shù)百株植株的冠層圖像僅需數(shù)分鐘,且操作簡單,采用該方法可無損、快速地獲取苗期玉米群體圖像。在生長中后期,植株間的遮蔽逐漸加大,使得無法獲得植株中下部的結(jié)構(gòu)信息。為此,在航拍前去除了目標植株周圍的植株。結(jié)果表明,以單株或一行多株為目標進行航拍,均可獲取相對完整的冠層結(jié)構(gòu),器官尺度參數(shù)的提取精度較高。本研究采用的偽極點-Crust算法,通過添加少量的偽極點代替?zhèn)鹘y(tǒng)Crust算法中原始點云對應(yīng)的正負極點,可大幅度減少極點的數(shù)量。以1行4株玉米為例,原有的正負極點數(shù)為1 392 490個,采用偽極點-Crust算法后降為96個,計算效率較高。實現(xiàn)5 m×5 m范圍內(nèi)的苗期植株群體結(jié)構(gòu)模型的構(gòu)建僅需10 min左右,對于成熟期4株玉米小群體的模型構(gòu)建僅需1 min左右。未來可以在生育中后期探索更大目標群體的航拍方法,以提高獲取作物生長中后期冠層結(jié)構(gòu)信息的效率。該方法還可應(yīng)用于植株葉片相互遮蔽不是很嚴重的其他作物,如高粱、生長中前期的大豆等,以及部分果樹、蔬菜等,對這些作物的應(yīng)用還需要系統(tǒng)地評估。
在玉米生長的中后期,葉片相互遮蔽嚴重,使得葉片間的噪聲點較多,采用搜索k-近鄰方法去除噪聲點的效果不夠理想。需要進一步優(yōu)化去噪算法,以減少噪聲點對構(gòu)建精確冠層結(jié)構(gòu)模型的干擾。目前所構(gòu)建的冠層模型,還存在因點云缺失導(dǎo)致所構(gòu)建的葉片形態(tài)不完整的問題,需要采用算法予以改進。基于構(gòu)建的精確冠層結(jié)構(gòu)模型,可進一步計算葉片角度[33]并進行大田冠層空間光分布的模擬分析,為評估不同作物株型、篩選高產(chǎn)作物品種提供支撐。
采用超微小型無人機平臺獲取了大田玉米冠層前期和后期的航拍圖像,實現(xiàn)了冠層點云的快速重建。基于SVM分類器提取了目標冠層的三維點云,并實現(xiàn)了冠層噪聲點的去除。采用偽極點-Crust算法初步構(gòu)建了冠層結(jié)構(gòu)模型,通過剔除異常的三角面片并進行移動平滑,得到了精度較高的玉米冠層結(jié)構(gòu)模型。基于所構(gòu)建模型提取的株高、葉長、最大葉寬及葉面積與大田測量值有較好的一致性。本研究所采用的圖像獲取方法、點云預(yù)處理流程及冠層結(jié)構(gòu)建模的算法效率較高,并可提取較精確的表型參數(shù),為高效、精確獲取大田作物冠層結(jié)構(gòu)信息提供了新的途徑。