結(jié)合激光雷達和三維性狀分析的田間小麥產(chǎn)量分級研究

施凌天,丁國輝,夏云鵬,葛玉峰,周濟,3*

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

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

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

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

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

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

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

1 材料與方法

1.1 田間試驗設(shè)計

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

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

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

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

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

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

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

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

1.5 點云的分割和校準算法設(shè)計

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

圖2 基于背包式激光雷達生成的大田尺度冠層高度模型(CHM)對小區(qū)作物點云分割及校準的算法流程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 大田點云分割為減少重復(fù)獲取小區(qū)尺度的地理坐標,本算法結(jié)合三維點云與CHM完成了對大田尺度的三維點云的小區(qū)分割:首先,使用Gaussian函數(shù)對代表田塊的CHM進行高斯濾波,在模糊小區(qū)細節(jié)的基礎(chǔ)上,減少細節(jié)對小區(qū)輪廓檢測的干擾;其次,使用threshold_sauvola函數(shù)初步將地面和非地面區(qū)域分開(圖2-b上部),以此繪制田塊中各個小區(qū)的邊緣輪廓;最后,使用hough_line函數(shù)［21］(設(shè)置檢測角度間隔90°,圖2-b),通過直線檢測小區(qū)邊緣。

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

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

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

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

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

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

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

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

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

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

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

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

(1)

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

1.8 人工檢測和算法評估

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

2 結(jié)果與分析

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

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

圖5 各生長相關(guān)性狀變化曲線及其關(guān)鍵采樣時期與小區(qū)產(chǎn)量及人工測量數(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)鍵生育時期,以完成產(chǎn)量分級,將人工測量的小區(qū)產(chǎn)量與各生育時期的冠層性狀進行線性回歸分析,用R2來評估關(guān)鍵生育時期與小區(qū)產(chǎn)量的相關(guān)性(圖5-b)。拔節(jié)期(150 DAS)小麥生長迅速,不同品種間冠層結(jié)構(gòu)差異較大,造成大部分小區(qū)三維冠層表面積較小,部分品種的三維冠層表面積較大。所以,該階段的三維冠層表面積數(shù)據(jù)分布較為分散(圖5-b;第3行,150 DAS)。此后的生育階段,小麥冠層逐漸緊密且結(jié)構(gòu)差異降低,三維冠層表面積數(shù)據(jù)分布逐漸緊密(圖5-b;第3行,170～197 DAS)。人工測量的單位面積產(chǎn)量與4個冠層性狀的相關(guān)性均在185 DAS(灌漿期,其中N0處理的小麥處于灌漿中后期,N180、N270處理的小麥處于灌漿中期;圖5-b中紅色虛線框)達到最大值,包括冠層高度(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)量等級的關(guān)鍵生育時期。

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

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

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

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

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

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

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

3 結(jié)論

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

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

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

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