施凌天,丁國(guó)輝,夏云鵬,葛玉峰,周濟(jì),3*

(1.南京農(nóng)業(yè)大學(xué)前沿交叉研究院/人工智能學(xué)院/作物表型組學(xué)交叉研究中心,江蘇 南京 210095;2.美國(guó)內(nèi)布拉斯加大學(xué)林肯分校工程學(xué)院/農(nóng)業(yè)與自然資源研究所,林肯 203 LW,美國(guó);3.英國(guó)國(guó)立農(nóng)業(yè)植物研究所/劍橋作物研究中心,劍橋 CB3 0LE,英國(guó))

小麥?zhǔn)鞘澜缟现饕募Z食作物之一,全球約有35%～40%的人口以小麥為主糧,維持其產(chǎn)量對(duì)世界糧食安全和經(jīng)濟(jì)發(fā)展均十分重要［1］。小麥產(chǎn)量是多種產(chǎn)量要素共同作用的結(jié)果,可通過關(guān)鍵性狀(如穗數(shù)、冠層覆蓋度和三維葉面積等)表征,且與冠層相關(guān)［2］。小麥產(chǎn)量中大部分貢獻(xiàn)來自頂部冠層的光合作用［3］,冠層表型會(huì)對(duì)最終產(chǎn)量產(chǎn)生較大影響,因此分析產(chǎn)量與關(guān)鍵產(chǎn)量性狀之間的關(guān)聯(lián)對(duì)小麥高產(chǎn)育種、栽培和農(nóng)業(yè)生產(chǎn)都具有重要意義。研究與產(chǎn)量相關(guān)的性狀還可為尋找不同環(huán)境條件下導(dǎo)致產(chǎn)量差異的原因提供依據(jù),為選育不同氣候下農(nóng)藝性狀更優(yōu)良的小麥品種,進(jìn)而為推廣高產(chǎn)栽培技術(shù)提供精準(zhǔn)表型分型技術(shù)和大數(shù)據(jù)支撐［4］。

近年來,隨著田間作物表型組研究的不斷發(fā)展［5］,高精度、非破壞性的表型采集技術(shù)和圖像處理、人工智能分析方法得到廣泛應(yīng)用［6］,這些技術(shù)方法在揭示植物生長(zhǎng)和產(chǎn)量等性狀方面發(fā)揮了關(guān)鍵作用。例如Bauer等［7］使用載人輕型飛機(jī)獲取超大規(guī)模生菜產(chǎn)區(qū)的航空影像,再結(jié)合計(jì)算機(jī)視覺和深度學(xué)習(xí)分析種植區(qū)內(nèi)生菜冠層的歸一化植被指數(shù)反射率,以此完成生菜計(jì)數(shù)和大小分級(jí),進(jìn)而完成對(duì)市場(chǎng)銷量(marketable yield)的預(yù)測(cè);Fernandez-Gallego等［8］使用手持熱像儀獲取植株冠層熱成像,再采用對(duì)比度增強(qiáng)和濾波技術(shù)計(jì)算麥穗數(shù)量,以此對(duì)小麥產(chǎn)量進(jìn)行預(yù)測(cè);Khan等［9］結(jié)合地面移動(dòng)平臺(tái)(mobile ground platform,MGP)和無(wú)人機(jī)進(jìn)行RGB(red-green-blue,即紅綠藍(lán)可見光)圖像采集,分析作物冠層特征,并使用冠層性狀解析春小麥在不同施肥處理下的表型差異。以上這些方法基于植物冠層對(duì)產(chǎn)量相關(guān)性狀進(jìn)行了采集和分析,針對(duì)冠層性狀和產(chǎn)量表型等進(jìn)行了有益的嘗試。然而,這些方法一般基于二維RGB、多光譜或熱輻射圖像,無(wú)法獲取完整且高精度的作物冠層結(jié)構(gòu)信息,因此很難對(duì)產(chǎn)量形成(yield formation)時(shí)的冠層形態(tài)變化進(jìn)行更進(jìn)一步的三維結(jié)構(gòu)重現(xiàn)和動(dòng)態(tài)性狀分析。

與其他采集方法不同,光探測(cè)和測(cè)距(light detection and ranging,LiDAR,即激光雷達(dá))具有穿透能力強(qiáng)、掃描范圍廣、靈活性高及對(duì)光線變化敏感性弱等優(yōu)點(diǎn),可快速獲取對(duì)象物體的三維結(jié)構(gòu)和形態(tài)特征［10］。根據(jù)LiDAR的應(yīng)用場(chǎng)景,其采集方式可大致分為3種類型:1)能生成高分辨率模型的定點(diǎn)地面激光雷達(dá)(terrestrial LiDAR),如Greaves等［11］利用Riegl VZ-400激光雷達(dá)采集灌木生物量和葉面積;2)可進(jìn)行高通量檢測(cè)的機(jī)載激光雷達(dá)(airborne LiDAR),如ten Harkel等［12］使用無(wú)人機(jī)搭載VUX-SYS激光雷達(dá)采集馬鈴薯、甜菜和冬小麥的冠層數(shù)據(jù),進(jìn)而估算生物量和作物高度;3)可靈活獲取作物三維形態(tài)的移動(dòng)激光雷達(dá)(mobile LiDAR),如Saeys等［13］將LMS 200和LMS 400安裝在聯(lián)合收割機(jī)上采集小麥密度和單位體積;Zhu等［14］結(jié)合背包式LiDAR和CropQuant-3D軟件對(duì)不同施氮水平下的小麥三維形態(tài)建立了基于氮素利用率的品種分類和鑒選流程。

以上LiDAR設(shè)備的使用為田間作物表型采集提供了可靠方法,但在針對(duì)冠層三維表型的采集中仍存在各類問題［15］。例如:1)地面激光雷達(dá)的數(shù)據(jù)采集效率較差,在測(cè)量占地面積大且地形復(fù)雜的田間試驗(yàn)時(shí)難以獲取所有作物的形態(tài);2)機(jī)載激光雷達(dá)的精度和兼容性較低,且飛行時(shí)產(chǎn)生的氣流還會(huì)影響作物冠層的形態(tài),造成采集誤差;3)移動(dòng)激光雷達(dá)中,手持式LiDAR雖易于攜帶使用,但構(gòu)建的三維模型質(zhì)量較差,很難完成大規(guī)模、高質(zhì)量的田間作物測(cè)繪;4)車載LiDAR操作要求和維護(hù)成本較高,且較難進(jìn)入某些地形進(jìn)行數(shù)據(jù)采集;5)背包式激光雷達(dá)雖然可在一定程度上克服上述限制,能快速靈活地在不同地點(diǎn)獲取相對(duì)高精度的作物三維形態(tài)數(shù)據(jù),但現(xiàn)階段其使用仍停留在單個(gè)生育時(shí)期的尺度［14］。

對(duì)作物三維點(diǎn)云的表型分析算法也是限制田間產(chǎn)量研究的瓶頸之一。例如:Madec等［16］結(jié)合激光雷達(dá)和無(wú)人機(jī)對(duì)灌溉和水分脅迫條件下的小麥表型進(jìn)行測(cè)量,基于株高對(duì)地上部生物量和產(chǎn)量進(jìn)行預(yù)測(cè),但該方法在對(duì)密集小區(qū)分析時(shí)準(zhǔn)確性低且單一性狀很難全面描述復(fù)雜田間環(huán)境下的產(chǎn)量變化;Jimenez-Berni等［17］結(jié)合SICK LMS 400激光雷達(dá)和Phenomobile Lite移動(dòng)平臺(tái)采集田間小麥表型數(shù)據(jù),并使用三維體素指數(shù)(3D voxel index,3DVI)和三維剖面指數(shù)(3D profile index,3DPI)估算地上生物量,然而該方法無(wú)法解釋灌漿期生物量變化,對(duì)冠層結(jié)構(gòu)的分析精度較低,并明顯低估了地上部作物生物量。值得指出的是,目前作物三維表型分析大都依靠單個(gè)性狀進(jìn)行評(píng)估,造成了算法對(duì)不同生育時(shí)期的分析結(jié)果偏差大、試驗(yàn)結(jié)果重復(fù)性低等問題［18］。

針對(duì)上述國(guó)內(nèi)外研究中的問題,本研究借鑒了Zhu等［14］通過背包式激光雷達(dá)構(gòu)建田間小麥三維點(diǎn)云和氮素響應(yīng)性狀分析的算法流程,優(yōu)化后提出了一套通過配準(zhǔn)修復(fù)錯(cuò)誤三維點(diǎn)云分類,再結(jié)合半自動(dòng)性狀處理進(jìn)行多生育時(shí)期、冠層性狀分析和田間產(chǎn)量分級(jí)的算法流程。為驗(yàn)證分析結(jié)果,本研究通過人工統(tǒng)計(jì)的產(chǎn)量數(shù)據(jù)和相關(guān)性分析對(duì)4個(gè)冠層形態(tài)性狀進(jìn)行驗(yàn)證,并綜合分析了影響產(chǎn)量分級(jí)的表型特征,為自動(dòng)化三維性狀提取和大規(guī)模田間產(chǎn)量分級(jí)提供了新思路和解析算法。

1 材料與方法

1.1 田間試驗(yàn)設(shè)計(jì)

根據(jù)人工產(chǎn)量測(cè)量結(jié)果,以54個(gè)中國(guó)長(zhǎng)江中下游的冬小麥品種(如鎮(zhèn)麥系、寧麥系和揚(yáng)麥系等)為試驗(yàn)材料,開展小麥氮素響應(yīng)和產(chǎn)量分級(jí)評(píng)估田間試驗(yàn),試驗(yàn)在鎮(zhèn)江市農(nóng)業(yè)科技創(chuàng)新中心進(jìn)行。采用裂區(qū)試驗(yàn)設(shè)計(jì),以3種氮素處理(0、180和270 kg·hm-2)為主區(qū)(9個(gè)田塊),54個(gè)小麥品種為副區(qū)(即田塊中的小區(qū)),設(shè)置3次重復(fù),共486個(gè)小區(qū)(圖1-a)。氮肥以基肥和拔節(jié)肥按照3∶2(質(zhì)量比)的比例施用,磷肥(P2O5)和鉀肥(K2O)均按照120 kg·hm-2作為底肥在播種前一次性施入。小麥播種密度為每公頃240萬(wàn)株基本苗,試驗(yàn)小區(qū)長(zhǎng)3 m,寬2 m,行間距25 cm。

圖1 使用背包式激光雷達(dá)對(duì)3種施氮水平下的54個(gè)小麥品種(3個(gè)重復(fù),共約2 916 m2)進(jìn)行三維點(diǎn)云采集和預(yù)處理Fig.1 The 3D crop phenotyping of wheat field experiments using backpack LiDAR together with data pre-processing, which has collected 3D point clouds from 54 wheat varieties(three biological replicates,approximately 2 916 m2) under three levels of nitrogen fertilization treatmentsN0、N180和N270:施氮量為0、180和270 kg·hm-2 Nitrogen application rates of 0,180 and 270 kg·hm-2. 下同 The same as follows.

1.2 表型數(shù)據(jù)采集

本研究所使用的田間表型采集設(shè)備是背包式LiDAR ROBIN(GeoSLAM公司,英國(guó);圖1-b)。該設(shè)備集成了RIEGK VUX-1 UAV激光雷達(dá)和緊湊型慣性導(dǎo)航系統(tǒng)(inertial navigation system),可測(cè)量精度為5 mm,田間測(cè)量準(zhǔn)確度約為10 mm,測(cè)量范圍為3～1 000 m。田間數(shù)據(jù)采集的基本流程包括:設(shè)備組裝和校準(zhǔn),按照預(yù)定采集路徑以1.5 m·s-1步行進(jìn)行表型采集(圖1-a);按行進(jìn)計(jì)劃完成3組路線的測(cè)繪(即紅、黃、藍(lán)色標(biāo)注路徑,圖1-a),再完成路徑校準(zhǔn);在田間完成作物測(cè)繪數(shù)據(jù)保存。為防止附近基站干擾LiDAR造成地理坐標(biāo)偏移,在離試驗(yàn)田較遠(yuǎn)的固定地點(diǎn)架設(shè)RTK(Real-time kinematic)基站。為保證數(shù)據(jù)質(zhì)量,采集時(shí)分3列行進(jìn),對(duì)單列田塊(即豎向的3個(gè)田塊,圖1-a)以S型路線繞行,平均單列采集時(shí)間約10～15 min,對(duì)整個(gè)田間試驗(yàn)采集耗時(shí)為1.5 h左右。在2021—2022年小麥生長(zhǎng)季中,對(duì)關(guān)鍵生育時(shí)期(從拔節(jié)到成熟期)的田間小麥分7次采集,分別為4月4日、4月10日、4月18日、4月24日、4月29日、5月9日和5月21日。期間不同小麥品種冠層特征變化和差異均較為明顯,符合本研究對(duì)冠層性狀表型差異的試驗(yàn)要求。

1.3 算法的軟件實(shí)現(xiàn)和測(cè)試

本研究針對(duì)三維冠層性狀研發(fā)了一系列分析算法流程,為提高下文中算法描述和軟件實(shí)現(xiàn)的可讀性和簡(jiǎn)潔性,本節(jié)集中介紹了實(shí)現(xiàn)算法使用的開源科學(xué)軟件開發(fā)庫(kù)。軟件實(shí)現(xiàn)通過編程語(yǔ)言Python(3.7)在Windows 10計(jì)算機(jī)上完成。其中,在處理三維點(diǎn)云數(shù)據(jù)時(shí)使用了Laspy(1.7.0)、Whitebox(1.3.0)、GDAL(3.1.4)、Rasterio(1.1.8)和Open3D(0.10.0.1)庫(kù);在處理二維圖像時(shí)使用了scikit-image(0.18.3)和OpenCV-Python(3.4.16.59)庫(kù);在對(duì)不同數(shù)據(jù)結(jié)構(gòu)(如數(shù)組和柵格)處理時(shí)使用了Numpy(1.19.4)、Pandas(1.3.4)和Scipy(1.7.1)庫(kù)［6］。通過PyCharm集成開發(fā)環(huán)境(integrated development environment,IDE)完成了源代碼編寫和結(jié)果測(cè)試。下文提及的各類函數(shù)和數(shù)據(jù)處理均通過上述科學(xué)軟件開發(fā)庫(kù)實(shí)現(xiàn)。

算法的軟件運(yùn)行與測(cè)試由2部分組成:1)田塊尺度,對(duì)單個(gè)田塊的點(diǎn)云文件［421～688 MB,共(1.5～2.6)×107個(gè)三維空間點(diǎn)］進(jìn)行小區(qū)分割,平均用時(shí)18～30 min,占用內(nèi)存15～16 G;2)小區(qū)尺度,對(duì)單個(gè)田塊內(nèi)54個(gè)小區(qū)的點(diǎn)云文件［351～586 MB,單個(gè)小區(qū)(0.01～0.12)×107個(gè)三維點(diǎn)］批量進(jìn)行數(shù)據(jù)優(yōu)化和性狀分析,平均用時(shí)為6～10 min,占用內(nèi)存150～640 MB。除數(shù)據(jù)預(yù)處理和小區(qū)點(diǎn)云分割部分需人工選擇相應(yīng)的解決方案,其余流程均可自動(dòng)化完成。

1.4 數(shù)據(jù)預(yù)處理流程

表型數(shù)據(jù)采集完成后,通過激光雷達(dá)ROBIN的配套處理軟件(即MMPrcoess、AEROoffice和grafNav)進(jìn)行數(shù)據(jù)預(yù)處理(圖1-b)。表型數(shù)據(jù)包括優(yōu)化采集時(shí)的路徑數(shù)據(jù),去除信號(hào)斷點(diǎn)和干擾數(shù)據(jù),配合RTK獲取的精確校準(zhǔn)地理坐標(biāo)信息的全球?qū)Ш叫l(wèi)星系統(tǒng)(global navigation satellite system,GNSS)數(shù)據(jù),將生成的初始三維點(diǎn)云數(shù)據(jù)保存為.las文件類型(圖1-c上部)。然后,使用TerraSolid系列軟件中的TerraMatch和TerraScan等模塊對(duì)載入的初始點(diǎn)云進(jìn)行區(qū)域劃分(即分割成9個(gè)田塊)、游離噪點(diǎn)去除和點(diǎn)云層次初步分類等點(diǎn)云數(shù)據(jù)預(yù)處理步驟(圖1-c下部)。

針對(duì)預(yù)處理獲得的三維點(diǎn)云(圖1-d上部),使用spatial.cKDTree函數(shù)對(duì)單個(gè)田塊點(diǎn)云建立索引,將漂浮于主要點(diǎn)云周邊且大于最大噪點(diǎn)判定距離［19］的點(diǎn),使用公式Avd.+k*Std. (式中,某一空間點(diǎn)與其相鄰點(diǎn)之間的平均距離為Avd.,其標(biāo)準(zhǔn)差為Std.,掃描距離倍數(shù)為k)判定為離群點(diǎn)。通過判定離群點(diǎn)定義三維噪點(diǎn),去除后完成降噪處理［14］。然后,使用height_above_ground函數(shù)對(duì)去噪點(diǎn)云進(jìn)行濾波,將不平整的地面點(diǎn)云歸一化,得到接近真實(shí)田間小麥生長(zhǎng)的三維點(diǎn)云文件。為方便對(duì)小區(qū)點(diǎn)云的批量分割,本研究特別將點(diǎn)云文件從三維轉(zhuǎn)換成二維柵格(raster)數(shù)據(jù),以此生成冠層高度模型(canopy height model,CHM)。在設(shè)置好柵格模型分辨率等參數(shù)后,通過lidar_tin_gridding函數(shù)生成每個(gè)田塊的CHM,并完成數(shù)據(jù)預(yù)處理流程(圖1-d下部)。

1.5 點(diǎn)云的分割和校準(zhǔn)算法設(shè)計(jì)

基于預(yù)處理獲取的點(diǎn)云數(shù)據(jù),設(shè)計(jì)了一種半自動(dòng)化的點(diǎn)云分割、校準(zhǔn)和性狀提取流程(圖2-a),包括:小區(qū)點(diǎn)云分割及重新定位、基于小區(qū)的點(diǎn)云分類、方向校準(zhǔn)和冠層性狀提取等步驟。首先,根據(jù)CHM轉(zhuǎn)換為地理坐標(biāo),完成小區(qū)范圍內(nèi)的點(diǎn)云分割;其次,使用Harris等［20］的角點(diǎn)檢測(cè),自動(dòng)化獲取小區(qū)點(diǎn)云的4個(gè)頂點(diǎn);最后,對(duì)CHM的方向校準(zhǔn),以便后續(xù)的自動(dòng)化性狀分析。在處理過程中,因?yàn)辄c(diǎn)云初步分類與田間實(shí)際情況存在偏差,算法還增加了點(diǎn)云分類修正環(huán)節(jié)。

圖2 基于背包式激光雷達(dá)生成的大田尺度冠層高度模型(CHM)對(duì)小區(qū)作物點(diǎn)云分割及校準(zhǔn)的算法流程Fig.2 The algorithmic workflow of plot-based segmentation and 3D point clouds calibration based on field-level canopy height model(CHM)acquired by the backpack LiDAR

1.5.1 大田點(diǎn)云分割為減少重復(fù)獲取小區(qū)尺度的地理坐標(biāo),本算法結(jié)合三維點(diǎn)云與CHM完成了對(duì)大田尺度的三維點(diǎn)云的小區(qū)分割:首先,使用Gaussian函數(shù)對(duì)代表田塊的CHM進(jìn)行高斯濾波,在模糊小區(qū)細(xì)節(jié)的基礎(chǔ)上,減少細(xì)節(jié)對(duì)小區(qū)輪廓檢測(cè)的干擾;其次,使用threshold_sauvola函數(shù)初步將地面和非地面區(qū)域分開(圖2-b上部),以此繪制田塊中各個(gè)小區(qū)的邊緣輪廓;最后,使用hough_line函數(shù)［21］(設(shè)置檢測(cè)角度間隔90°,圖2-b),通過直線檢測(cè)小區(qū)邊緣。

1.5.2 修正錯(cuò)誤邊界小區(qū)邊緣檢測(cè)中會(huì)出現(xiàn)一些錯(cuò)誤,如檢測(cè)直線過多、直線斜率錯(cuò)誤和部分邊界缺失等(圖2-c)。需通過算法改進(jìn)檢測(cè)結(jié)果:1)針對(duì)檢測(cè)直線過多問題,以間距閾值(即相鄰2行或2列間隔的0.75倍)為篩選條件,刪除小于此閾值的小區(qū)行列線;2)針對(duì)直線斜率錯(cuò)誤問題,計(jì)算所有檢測(cè)直線斜率的中位數(shù),以此刪除錯(cuò)誤檢測(cè)線;3)針對(duì)部分直線無(wú)法檢測(cè)問題(如圖2-b紅色圓圈位置),取最接近的直線坐標(biāo),以間距閾值分析空缺直線的代表坐標(biāo)點(diǎn),再結(jié)合行(列)斜率補(bǔ)充缺失直線和最邊緣直線。

1.5.3 分割小區(qū)點(diǎn)云完整獲取小區(qū)邊界后(圖2-d上部),利用小區(qū)行(列)線提取各小區(qū)的4個(gè)頂點(diǎn)坐標(biāo)(圖2-d下部)。算法使用transform函數(shù)獲取仿射變換(affine transformation)矩陣,將這些頂點(diǎn)的像素坐標(biāo)(pixel-based coordinates)轉(zhuǎn)換為對(duì)應(yīng)地理坐標(biāo)(geographic coordinates)。算法批量輸出與各小區(qū)頂點(diǎn)的地理坐標(biāo)對(duì)應(yīng)的Shapefile文件,即ESRI Shapefile,最后使用clip_lidar_to_polygon函數(shù)完成對(duì)田塊中每個(gè)小區(qū)的三維點(diǎn)云的分割和保存。

1.5.4 小區(qū)點(diǎn)云的方向校準(zhǔn)分割各小區(qū)點(diǎn)云后,自動(dòng)校準(zhǔn)小區(qū)CHM的方向,以便后續(xù)自動(dòng)化性狀處理。首先使用cornerHarris函數(shù)提取出小區(qū)圖像的拐點(diǎn),結(jié)合圖像邊緣值篩選出4個(gè)頂點(diǎn)(圖2-e);然后使用getPerspectiveTransform函數(shù)對(duì)圖像進(jìn)行二維透視變換,通過迭代完成對(duì)象田塊內(nèi)所有54個(gè)小麥小區(qū)的方向校準(zhǔn)。

1.6 小區(qū)點(diǎn)云的重分類

基于小區(qū)的三維點(diǎn)云會(huì)存在分類異常情況,例如由于植物遮掩或采集數(shù)據(jù)時(shí)的氣候情況,某些小區(qū)中間部分會(huì)呈現(xiàn)出點(diǎn)云凹陷或密度不均等問題。基于這些異常小區(qū)完成的冠層性狀分析(如高度或冠層結(jié)構(gòu))會(huì)出現(xiàn)與實(shí)際數(shù)據(jù)偏差較大的情況。因此,針對(duì)小區(qū)點(diǎn)云分類偏差構(gòu)建了分類修正算法(圖3-a),包括分類三維點(diǎn)云、矯正點(diǎn)云分類、修正點(diǎn)云和批量處理等步驟。

圖3 通過小區(qū)作物空間特征分析完成小麥作物的三維點(diǎn)云分類重校準(zhǔn)Fig.3 The re-classification of 3D point clouds based on wheat plants’ spatial features at the plot level

1.6.1 點(diǎn)云分類錯(cuò)誤在預(yù)處理過程中,TerraSolid會(huì)對(duì)某些三維點(diǎn)云錯(cuò)誤分類,其原因包括點(diǎn)云稀疏或小區(qū)中心區(qū)域三維空間點(diǎn)較少,造成地面點(diǎn)云部分過少。因此,濾波處理會(huì)將冠層部分點(diǎn)云沉降并誤分類為地面點(diǎn)。此外,由于地面起伏,一些接近地面部分的小麥點(diǎn)云也會(huì)被誤分類為地面點(diǎn)(ground level points),一些地面點(diǎn)云也會(huì)被分類為低植被點(diǎn)(low vegetation points,圖3-b)。

1.6.2 修正三維點(diǎn)云為保證性狀分析的準(zhǔn)確性,構(gòu)建了自動(dòng)化流程,完成對(duì)去噪后小區(qū)三維點(diǎn)云重新校準(zhǔn)分類,包括:1)按照TerraSolid預(yù)處理的分類,即未分類點(diǎn)(unclassified)、地面點(diǎn)(ground)、低植被點(diǎn)(low vegetation)、中植被點(diǎn)(medium vegetation)和高植被點(diǎn)(high vegetation),將點(diǎn)云重新進(jìn)行歸類;2)根據(jù)點(diǎn)云的空間特征信息,查找錯(cuò)誤分類(即小麥最高高度約為1.5 m,因此可去除高植被點(diǎn)分類;中植被點(diǎn)一般正確;低植被點(diǎn)、地面點(diǎn)和未分類點(diǎn)需重新歸類,圖3-c);3)根據(jù)小麥各階段高度設(shè)置閾值,對(duì)低植被點(diǎn)和地面點(diǎn)重新劃定分類,即基于大田的平均地面高度統(tǒng)一劃定該田塊內(nèi)的地面高度閾值,在比較不同時(shí)期地面高度后,使用大于99%地面三維點(diǎn)的高度為閾值(H99th),將地面點(diǎn)和植被點(diǎn)重新劃分(圖3-d);4)完成地面點(diǎn)分類后,對(duì)未分類點(diǎn)、地面點(diǎn)和低植被點(diǎn)重新進(jìn)行歸類,即低于閾值的三維點(diǎn)為地面點(diǎn),高于閾值的為低植被點(diǎn);5)以此流程,基于小區(qū)對(duì)所有小麥點(diǎn)云進(jìn)行批量修正,完成濾波處理和CHM提取(圖3-e)。利用俯視圖對(duì)小區(qū)點(diǎn)云切割后,小麥冠層會(huì)遮擋住冠層或高植被點(diǎn)下方的三維點(diǎn)云,因此對(duì)小區(qū)切割時(shí)都會(huì)保存小麥下部點(diǎn)云(即低植被點(diǎn))或地面點(diǎn),以保證修正三維點(diǎn)云。此外,大田中部小麥冠層下方點(diǎn)云稀少并不會(huì)影響后續(xù)基于冠層的性狀分析結(jié)果。

1.7 基于冠層的形態(tài)性狀分析

對(duì)校正后的點(diǎn)云和對(duì)應(yīng)的CHM可提取冠層性狀,如作物高度(crop height)、冠層覆蓋度(canopy coverage)、三維冠層表面積(3D canopy surface area)和三維冠層指數(shù)(3D canopy index,3DCI)［17］。對(duì)于作物冠層高度,本研究選取了大于小區(qū)90%高度的三維點(diǎn)［22-23］(H90th),并將這些三維點(diǎn)的高度平均值作為小區(qū)冠層平均高度,以減少小區(qū)植株高度差異造成的誤差(圖4-a、b)。對(duì)于冠層性狀,以最高高度的60%為閾值高度提取冠層部分點(diǎn)云(圖4-c);然后基于冠層點(diǎn)云俯視圖(圖4-d),使用threshold_local函數(shù)得到對(duì)應(yīng)小區(qū)冠層區(qū)域的二值圖(圖4-g),以此計(jì)算小區(qū)尺度的覆蓋度(0～1,其值越接近1,表明小區(qū)冠層覆蓋度越高)。對(duì)3DCI性狀的計(jì)算,通過二維傅里葉變換［24］(discrete Fourier transform,DFT)將獲取的小區(qū)圖像轉(zhuǎn)換成頻域(頻譜圖),并將零頻率分量移至頻譜中心(圖4-e)。然后,使用高斯擬合對(duì)頻譜圖中高度值的頻率和振幅進(jìn)行擬合并繪制曲線(圖4-h),使用trapz函數(shù)獲取曲線下面積,歸一化處理后獲得3DCI。如圖4-h所示,‘寧麥13’在3種氮素處理下的3DCI曲線差異明顯,表明3DCI值越小,冠層越密集植株高度差異越小,反之則植株高度及冠層結(jié)構(gòu)差異大。

圖4 基于小區(qū)的小麥關(guān)鍵生長(zhǎng)性狀分析算法流程Fig.4 The algorithmic workflow of key growth-related traits of wheat measured at the plot level

對(duì)于三維冠層表面積性狀,使用voxel_down_sample函數(shù)對(duì)冠層點(diǎn)云進(jìn)行體素化處理(圖4-f),降低點(diǎn)云密度后完成關(guān)鍵三維空間點(diǎn)的標(biāo)定;再使用creat_from_point_cloud_alpha_shape函數(shù)重建冠層三角網(wǎng)格模型(圖4-i);最后使用get_surface_area函數(shù)得到小區(qū)內(nèi)小麥冠層的三維表面積。由于不同小區(qū)點(diǎn)云數(shù)量的不同,點(diǎn)云稀少的小區(qū)冠層模型會(huì)出現(xiàn)不完整情況,所以通過式(1)計(jì)算了三維冠層表面積。

(1)

式中:SCS表示三維冠層表面積性狀;SCM表示體素化后冠層模型的表面積;SP表示小區(qū)面積;k為不完整的冠層模型表面積系數(shù),其取值范圍為0.6～0.7(即接近1個(gè)標(biāo)準(zhǔn)差的范圍)。

1.8 人工檢測(cè)和算法評(píng)估

人工數(shù)據(jù)集(ground truth data)由2部分組成:1)通過人工統(tǒng)計(jì)獲取田間不同品種的株高、單位面積穗數(shù)(SN m-2)和單位面積產(chǎn)量。對(duì)田間獲取的人工數(shù)據(jù)采取標(biāo)準(zhǔn)化處理,包括:交叉計(jì)算、歷史數(shù)據(jù)驗(yàn)證和截尾均值去除離群值;2)使用ImageJ軟件測(cè)量小區(qū)的冠層高度和冠層覆蓋度等。將算法獲取的作物高度、冠層覆蓋度、三維冠層表面積和3DCI等性狀與人工統(tǒng)計(jì)的產(chǎn)量和冠層性狀數(shù)據(jù)進(jìn)行相關(guān)性分析,使用決定系數(shù)R2、均方根誤差(root mean square error,RMSE)和P值等綜合評(píng)價(jià)算法獲取性狀的可靠性。

2 結(jié)果與分析

2.1 冠層性狀的變化及與人工數(shù)據(jù)的量化分析

2.1.1 冠層性狀的動(dòng)態(tài)變化及產(chǎn)量分級(jí)的關(guān)鍵時(shí)期通過7個(gè)時(shí)間點(diǎn)的4個(gè)冠層性狀(即作物高度、冠層覆蓋度、三維冠層表面積和3DCI),分析了從拔節(jié)到灌漿后期(播種后150～197 d;day after sowing,DAS)54個(gè)小麥品種的冠層性狀變化過程(圖5-a)。在3種氮素處理下,冠層性狀變化趨勢(shì)相似:作物高度和冠層覆蓋度呈上升趨勢(shì),在185 DAS到達(dá)最高值;三維冠層表面積和3DCI持續(xù)增加,在185 DAS后上升變緩。4個(gè)目標(biāo)性狀在高氮(270 kg·hm-2,N270)和中氮(180 kg·hm-2,N180)處理的變化曲線接近,與低氮處理(0 kg·hm-2,N0)的表型變化差異較大。

圖5 各生長(zhǎng)相關(guān)性狀變化曲線及其關(guān)鍵采樣時(shí)期與小區(qū)產(chǎn)量及人工測(cè)量數(shù)據(jù)的相關(guān)性分析Fig.5 The curves of key growth-related traits and their key sampling date together with the correlation analysis between these traits and grain production as well as manual measurements

為了確定與最終產(chǎn)量相關(guān)性較高的關(guān)鍵生育時(shí)期,以完成產(chǎn)量分級(jí),將人工測(cè)量的小區(qū)產(chǎn)量與各生育時(shí)期的冠層性狀進(jìn)行線性回歸分析,用R2來評(píng)估關(guān)鍵生育時(shí)期與小區(qū)產(chǎn)量的相關(guān)性(圖5-b)。拔節(jié)期(150 DAS)小麥生長(zhǎng)迅速,不同品種間冠層結(jié)構(gòu)差異較大,造成大部分小區(qū)三維冠層表面積較小,部分品種的三維冠層表面積較大。所以,該階段的三維冠層表面積數(shù)據(jù)分布較為分散(圖5-b;第3行,150 DAS)。此后的生育階段,小麥冠層逐漸緊密且結(jié)構(gòu)差異降低,三維冠層表面積數(shù)據(jù)分布逐漸緊密(圖5-b;第3行,170～197 DAS)。人工測(cè)量的單位面積產(chǎn)量與4個(gè)冠層性狀的相關(guān)性均在185 DAS(灌漿期,其中N0處理的小麥處于灌漿中后期,N180、N270處理的小麥處于灌漿中期;圖5-b中紅色虛線框)達(dá)到最大值,包括冠層高度(R2=0.836 2,RMSE=0.097 6)、冠層覆蓋度(R2=0.783 9,RMSE=0.221 4)、三維冠層表面積(R2=0.771 8,RMSE=0.151 7)和3DCI(R2=0.738 3,RMSE=0.451 5),且所有P值小于0.001。所以,灌漿期被用于在田間判定產(chǎn)量等級(jí)的關(guān)鍵生育時(shí)期。

此外,為了解不同氮素處理的產(chǎn)量差異,以此量化分析高產(chǎn)小麥的冠層性狀特征。本研究還以不同顏色的散點(diǎn)對(duì)不同氮素處理的小麥單位面積產(chǎn)量進(jìn)行劃分:低氮處理(藍(lán)色散點(diǎn))的小麥單位面積產(chǎn)量較低,且與中、高氮處理基本沒有交集;中氮(綠色散點(diǎn))和高氮處理(紅色散點(diǎn))的單位面積產(chǎn)量分布類似,因此選擇將中、高氮處理的小麥品種合并,以此進(jìn)行后續(xù)的高產(chǎn)性狀分類研究。

2.1.2 關(guān)鍵生育時(shí)期生長(zhǎng)性狀與人工數(shù)據(jù)對(duì)比在進(jìn)行產(chǎn)量分級(jí)前,為了確保算法計(jì)算小麥冠層性狀的可靠性,在185 DAS選取算法獲取的4個(gè)冠層性狀與人工測(cè)量數(shù)據(jù)進(jìn)行線性回歸分析(圖5-c):算法得出的作物高度和人工測(cè)量數(shù)據(jù)的決定系數(shù)R2為0.866(RMSE=5.66 cm,P<0.001),算法得出的冠層覆蓋度和人工測(cè)量數(shù)據(jù)的R2為0.899 3(RMSE=0.057 4,P<0.001),算法得出的三維冠層表面積和人工測(cè)量數(shù)據(jù)的R2為0.836 4(RMSE=0.170 3,P<0.001),算法得出的3DCI和人工測(cè)量穗數(shù)的R2為0.769 5(RMSE=0.265 5,P<0.001)。相關(guān)性分析表明,本研究提出的算法分析獲得的冠層性狀與人工測(cè)量的性狀和產(chǎn)量數(shù)據(jù)均具有顯著相關(guān)性,因此能用于量化描述分蘗期后不同時(shí)期小麥冠層性狀和不同品種的產(chǎn)量差異。

2.2 不同小麥品種的冠層性狀差異

2.2.1 基于產(chǎn)量分級(jí)的小麥性狀差異根據(jù)小麥冠層性狀和單位面積產(chǎn)量在N180和N270處理下相似度高的特點(diǎn)(圖5-a、b),將N180和N270處理合并后對(duì)54個(gè)小麥品種完成了產(chǎn)量分級(jí)(圖6-a)。首先,根據(jù)曹寒冰等［25］通過人工統(tǒng)計(jì)的單位面積產(chǎn)量,將小麥品種劃分為3個(gè)產(chǎn)量等級(jí),基于箱線圖也對(duì)產(chǎn)量進(jìn)行了3個(gè)等級(jí)的劃分:Ⅰ級(jí)(紅色,>75%,共14個(gè)品種)、Ⅱ級(jí)(綠色,25%～75%,27個(gè)品種)、Ⅲ級(jí)(藍(lán)色,<25%,12個(gè)品種,1個(gè)離散點(diǎn))。其次,基于3個(gè)產(chǎn)量等級(jí),通過雷達(dá)圖對(duì)比了其所屬品種在185 DAS時(shí)的冠層性狀差異(圖6-b):作物高度分布范圍為75～95 cm,冠層覆蓋度分布范圍為0.55～0.65,三維冠層表面積分布范圍為9～11 m2,3DCI分布范圍為0.65～0.75。以上結(jié)合不同產(chǎn)量等級(jí)對(duì)性狀差異進(jìn)行分析,可更好了解不同品種產(chǎn)量分布情況,為了解影響產(chǎn)量的關(guān)鍵因素提供支撐。其中,產(chǎn)量Ⅰ級(jí)(代表高產(chǎn)小麥品種)的高度、覆蓋度和冠層表面積數(shù)值較大,3DCI較小(即小麥冠層緊密且三維結(jié)構(gòu)差異小);產(chǎn)量Ⅱ級(jí)的品種4個(gè)冠層性狀數(shù)值對(duì)比其他2個(gè)等級(jí)適中;產(chǎn)量Ⅲ級(jí)(代表相對(duì)低產(chǎn)品種)的高度、覆蓋度和冠層表面積數(shù)值較小,3DCI較大(即小麥冠層松散且高度差異較大)。

圖6 根據(jù)產(chǎn)量分級(jí)和性狀聚類分析對(duì)N180/N270氮素處理下53個(gè)品種的4個(gè)冠層性狀的綜合分析Fig.6 A comprehensive comparison of canopy-related traits under N180/N270 N treatments according to clusters of 53 wheat varieties’ yield-production levels and four canopy-level traits

2.2.2 基于冠層性狀聚類的小麥性狀差異除了基于產(chǎn)量分級(jí)解析性狀差異,本研究還利用各品種小麥在185 DAS時(shí)的4個(gè)冠層性狀,使用K-means完成了聚類分析。根據(jù)性狀差異將53個(gè)小麥品種(同2.2.1節(jié)類似,去除了1個(gè)離散品種)分為3類:第1類群體有23個(gè)品種;第2類群體有19個(gè)品種;第3類群體有11個(gè)品種(圖6-c)。聚類分析中PC1(principal components)的貢獻(xiàn)率為74.35%,主要包括作物高度(特征值為0.52)、冠層覆蓋度(特征值為0.50)和三維冠層表面積(特征值為0.54);PC2的貢獻(xiàn)率為14.22%。根據(jù)聚類獲得的3個(gè)性狀分類(紅、綠、藍(lán)色,圖6-d),其高度、冠層覆蓋度、三維冠層表面積和3DCI的分布特征和范圍與基于產(chǎn)量分級(jí)獲得的雷達(dá)圖高度相似。

2.2.3 產(chǎn)量分級(jí)品種和性狀聚類品種的對(duì)比由于通過產(chǎn)量分級(jí)和性狀聚類獲取的冠層性狀特征分布具有很高的相似性,進(jìn)一步對(duì)比了產(chǎn)量Ⅰ級(jí)(紅色;高產(chǎn)小麥)對(duì)應(yīng)品種與性狀聚類第2類(綠色)對(duì)應(yīng)小麥品種。結(jié)果(表1)顯示,2種分類對(duì)應(yīng)的品種重合度為85.7%(12個(gè)品種),基本包含了在中、高氮處理下產(chǎn)量表現(xiàn)較優(yōu)異的已發(fā)表的小麥品種［26-27］。其余2類對(duì)應(yīng)的品種重合度也達(dá)到了66.7%(圖6-e)。綜上所述,在關(guān)鍵生育時(shí)期通過激光雷達(dá)獲取的三維點(diǎn)云分析冠層性狀,再使用具有高產(chǎn)特征的表型組(圖6-b、d)能夠完成對(duì)不同小麥品種的田間產(chǎn)量分級(jí),該方法具有很強(qiáng)的可行性和可操作性。

表1 在185 DAS時(shí)產(chǎn)量分級(jí)和性狀聚類對(duì)應(yīng)的相同品種Table 1 Same wheat varieties identified by yield-classification and trait clustering at 185 DAS

3 結(jié)論

1)利用背包式激光雷達(dá)在田間試驗(yàn)中采集了不同生育時(shí)期的54個(gè)小麥品種(3個(gè)重復(fù))的三維點(diǎn)云,基于前人研究研發(fā)了大田點(diǎn)云矯正、精準(zhǔn)小區(qū)點(diǎn)云分割、修正和提取冠層區(qū)域點(diǎn)云等三維點(diǎn)云分析流程,在小區(qū)尺度對(duì)不同小麥品種的多個(gè)冠層性狀進(jìn)行了基于時(shí)序的分析。該表型分析流程通過精準(zhǔn)分割小區(qū)點(diǎn)云和點(diǎn)云重分類,減少了田間地形和點(diǎn)云稀疏對(duì)自動(dòng)化提取冠層性狀的影響。在批量使用前,本研究提出的方法經(jīng)過對(duì)比篩選(如基于聚類或密度的去噪算法、基于坡度或插值的濾波算法等),能夠更簡(jiǎn)單、高效地完成三維點(diǎn)云處理與冠層性狀分析,與其他算法獲得的結(jié)果對(duì)比后獲得更接近田間真實(shí)情況的分析結(jié)果。但是由于數(shù)據(jù)預(yù)處理和小區(qū)邊界檢測(cè)仍需人工介入,后續(xù)研究將進(jìn)一步對(duì)這兩方面作出改進(jìn),以期構(gòu)建針對(duì)關(guān)鍵小麥冠層性狀的大規(guī)模、自動(dòng)化提取平臺(tái)。

2)對(duì)所試小麥品種完成了從拔節(jié)到灌漿后期的冠層性狀分析,即在3種氮素處理下作物高度、冠層覆蓋度、三維冠層表面積和3DCI的變化趨勢(shì)。結(jié)果表明,灌漿期(185 DAS)的單位面積產(chǎn)量與4個(gè)冠層性狀的相關(guān)性達(dá)到最大值,因此該時(shí)期可作為在田間對(duì)不同小麥品種進(jìn)行產(chǎn)量分級(jí)的關(guān)鍵生育時(shí)期。本研究中使用的小麥品種較多,品種間生育時(shí)期存在差異,而且灌漿期各階段對(duì)預(yù)測(cè)小麥產(chǎn)量的影響不同(灌漿初期判斷小麥的產(chǎn)量潛力,灌漿中期估計(jì)小麥的產(chǎn)量情況,灌漿后期觀察小麥的產(chǎn)量表現(xiàn))［28］。后續(xù)研究將綜合使用灌漿期不同階段的觀察結(jié)果,更全面評(píng)估小麥產(chǎn)量。通過對(duì)冠層性狀與人工測(cè)量數(shù)據(jù)的線性回歸分析,算法獲得的冠層性狀與人工測(cè)量數(shù)據(jù)具有顯著相關(guān)性,驗(yàn)證了本研究構(gòu)建的性狀分析算法的可靠性。然而,本方法的特征提取精度和應(yīng)用場(chǎng)景仍有提升和擴(kuò)展的空間,特別是針對(duì)其他谷類作物的測(cè)試仍處在嘗試階段,后續(xù)研究將在本研究的分析流程上進(jìn)行改進(jìn)。

3)通過產(chǎn)量分級(jí)和性狀聚類獲取的冠層性狀特征有很高的相似性。例如,產(chǎn)量Ⅰ級(jí)(即高產(chǎn)品種)與第2類(性狀聚類)對(duì)應(yīng)的品種重合度為85.7%;2類對(duì)應(yīng)品種的高度、覆蓋度和冠層表面積較大,3DCI較小(即小麥冠層緊密且三維結(jié)構(gòu)差異小),具有明顯的高產(chǎn)性狀特征。結(jié)果證明,基于背包式激光雷達(dá)采集并分析三維冠層性狀,在關(guān)鍵生育時(shí)期對(duì)不同小麥品種進(jìn)行產(chǎn)量分級(jí)具有一定的可行性和可操作性。

4)在關(guān)鍵生育時(shí)期對(duì)小麥產(chǎn)量分級(jí)能為遺傳育種中大規(guī)模鑒選優(yōu)質(zhì)高產(chǎn)小麥基因型,進(jìn)而解析作物在不同環(huán)境下產(chǎn)量的形成規(guī)律、降低高產(chǎn)性狀分析難度以及為品種選育提供量化依據(jù)。在實(shí)際栽培和農(nóng)業(yè)生產(chǎn)中,精確預(yù)測(cè)產(chǎn)量也可為田間動(dòng)態(tài)改進(jìn)栽培措施、精準(zhǔn)準(zhǔn)備收獲期的各項(xiàng)工作(如農(nóng)機(jī)、貨運(yùn)和倉(cāng)儲(chǔ)安排等)提供依據(jù)。特別是本研究提出的基于三維冠層性狀的分級(jí)方法能在作物生殖生長(zhǎng)(reproductive phase)早期構(gòu)建不同品種的產(chǎn)量等級(jí)劃分,在小麥遺傳育種、精準(zhǔn)栽培和農(nóng)業(yè)生產(chǎn)中均有較大的應(yīng)用價(jià)值。