基于作物生長模型與機(jī)器學(xué)習(xí)算法的區(qū)域冬小麥估產(chǎn)

馬戰(zhàn)林 文 楓 周穎杰 魯春陽 薛華柱 李長春

(1.河南城建學(xué)院測繪與城市空間信息學(xué)院, 平頂山 467001; 2.河南城建學(xué)院管理學(xué)院, 平頂山 467001;3.河南理工大學(xué)測繪與國土信息工程學(xué)院, 焦作 454003)

0 引言

河南省為中國小麥主產(chǎn)地之一,冬小麥產(chǎn)量約占全國小麥產(chǎn)量的25%。精準(zhǔn)、高效、實時的區(qū)域冬小麥產(chǎn)量預(yù)測,不僅能夠為河南省各級政府、相關(guān)部門制定農(nóng)業(yè)農(nóng)村政策提供技術(shù)支持,同時對保障中國糧食安全、農(nóng)業(yè)可持續(xù)發(fā)展、種植結(jié)構(gòu)優(yōu)化和數(shù)字農(nóng)業(yè)建設(shè)也有重要意義[1]。

當(dāng)前,作物估產(chǎn)常用的方法有人工區(qū)域調(diào)查法、作物生長模型及機(jī)器學(xué)習(xí)算法3類。人工區(qū)域調(diào)查方法工作量大且成本較高,難以實現(xiàn)大范圍、高頻的區(qū)域尺度作物估產(chǎn),但其優(yōu)勢是獲取的點源信息準(zhǔn)確度高。作物生長模型雖融入光、溫、水、土壤等條件為環(huán)境驅(qū)動變量,對作物生長的全過程、產(chǎn)量形成機(jī)理具有較好的科學(xué)解釋性,但僅在單點尺度上建立和實現(xiàn)作物生長發(fā)育動態(tài)模擬。遙感對地觀測的優(yōu)勢在于提供大空間、面元尺度的作物物候、生長冠層信息。實現(xiàn)遙感信息與作物生長模型的耦合,則能夠利用空間上連續(xù)、時間上動態(tài)變化的衛(wèi)星觀測數(shù)據(jù)得到模型中較難獲取的區(qū)域級尺度參數(shù),以校正、優(yōu)化作物生長模型與產(chǎn)量形成過程的相關(guān)參數(shù),使模擬的產(chǎn)量結(jié)果更為準(zhǔn)確[2]。結(jié)合遙感數(shù)據(jù)與作物生長模型實現(xiàn)區(qū)域估產(chǎn)的方法主要包含驅(qū)動法和同化法。相比于驅(qū)動法,同化法已成為近年來的研究熱點。

同化法以改進(jìn)或調(diào)整作物生長模型中與作物生長總生物量、糧食產(chǎn)量形成密切相關(guān)的輸入?yún)?shù)或初始條件,從而縮小作物生長模型模擬的和遙感數(shù)據(jù)反演的狀態(tài)變量間的差距,達(dá)到準(zhǔn)確估測作物生長模型參數(shù)和提高作物估產(chǎn)精度的目的。針對作物模型敏感性參數(shù)分析,常用的敏感度分析方法主要有Morris、FAST、Sobol、EFAST等方法[3]。對于同化方法,當(dāng)前同化法主要有擴(kuò)展卡爾曼濾波(Extended Kalman filter,EKF)、粒子濾波(Particle filter,PF)、集合卡爾曼濾波(Ensemble Kalman filter,EnKF)等順序同化算法和基于最優(yōu)控制的全局?jǐn)M合法,如四維變分(4DVAR)[4]、三維變分(3DVAR)及VW-4DenSRF[5](Variable time window and four-dimensional extension)等方法。其中,EnKF、4DVAR同化法最為常用。相關(guān)研究表明,EnKF在相對較小的研究區(qū)域上能夠取得較高的產(chǎn)量模擬精度[6]。如謝毅等[7]對比4DVAR和EnKF不同同化算法模擬關(guān)中平原冬小麥產(chǎn)量精度,結(jié)果表明EnKF算法在同化LAI時的產(chǎn)量估算精度較高。另外,同化變量的選擇同樣會影響參數(shù)的調(diào)整,最終影響作物產(chǎn)量預(yù)報精度[2]。當(dāng)前常用的同化變量主要是LAI,其次為土壤含水率、生物量等。

當(dāng)然,數(shù)據(jù)同化是一種將外部觀測數(shù)據(jù)納入動態(tài)力學(xué)模型的方法,同化單元尺寸的選擇不僅取決于遙感數(shù)據(jù)反演的同化參數(shù)分辨率,同樣也取決于氣象要素、作物、土壤及田間管理的分辨率。更細(xì)的同化單元則伴隨著巨大的計算力,需要研發(fā)更高效的同化策略及適用于高性能計算的組織架構(gòu)與模式。近年來機(jī)器學(xué)習(xí)方法已與遙感信息進(jìn)行融合,被廣泛應(yīng)用于作物估產(chǎn)研究[8]。機(jī)器學(xué)習(xí)算法具有很強(qiáng)的非線性擬合能力和高計算效率,但其是一種數(shù)據(jù)驅(qū)動的模型,大量的樣本才能使得該模型涵蓋盡可能的分布,從而提升其適用性。受作物育種技術(shù)、品種特性的影響,作物長勢監(jiān)測和估產(chǎn)模型構(gòu)建所需的樣本數(shù)據(jù)一般是近5年的數(shù)據(jù),造成作為訓(xùn)練樣本的數(shù)據(jù)量有限且具有很強(qiáng)的時效性[8]。且使用機(jī)器學(xué)習(xí)實現(xiàn)區(qū)域作物估產(chǎn),樣本采集會耗費大量人力物力。

從提升區(qū)域估產(chǎn)精度與效率的角度出發(fā),本文以河南省鶴壁市石橋村為研究區(qū),基于PROSAIL輻射傳輸模型反演多時相Sentinel-2光學(xué)遙感影像的時序LAI,利用EnKF算法同化時序LAI到作物生長模型中,提供不同長勢的一定數(shù)量單點產(chǎn)量訓(xùn)練數(shù)據(jù),利用建立的機(jī)器學(xué)習(xí)模型反演區(qū)域冬小麥產(chǎn)量,以期實現(xiàn)作物生長模型與機(jī)器學(xué)習(xí)算法的耦合,為準(zhǔn)確、實時、高效地實現(xiàn)區(qū)域作物估產(chǎn)提供新的模式與技術(shù)。

1 材料與方法

1.1 研究區(qū)域

鶴壁市是河南省重要的冬小麥種植區(qū),穩(wěn)定的氣候環(huán)境和少發(fā)的極端氣候為建立和驗證作物估產(chǎn)模型提供良好的氣候條件。為驗證建立的耦合作物生長模型與機(jī)器學(xué)習(xí)算法實現(xiàn)區(qū)域估產(chǎn)的有效性,該研究在鶴壁市中部地區(qū)淇縣橋盟鄉(xiāng)石橋村建立村級區(qū)域級別的試驗與算法驗證基地。該基地是鶴壁市重要的冬小麥種植區(qū)和高標(biāo)準(zhǔn)農(nóng)田建設(shè)示范基地,地勢平整、連片。試驗區(qū)占地309.23 hm2,是河南省農(nóng)業(yè)大數(shù)據(jù)中心的規(guī)?；芾矸N植區(qū)。研究區(qū)具體位置如圖1所示。

圖1 研究區(qū)地理位置Fig.1 Geographical location of study area

試驗區(qū)屬暖溫帶半濕潤型季風(fēng)氣候,四季分明,光照充足,年平均氣溫14.2～15.5℃,年降水量349.2～970.1 mm;年日照時數(shù)1 787.2～2 566.7 h。試驗區(qū)土壤類型以粘土為主, 灌溉條件良好, 主要種植模式為冬小麥和夏玉米輪作。冬小麥播種時間通常情況下在10月上、中旬,出苗期在播種后5～6 d。根據(jù)當(dāng)?shù)貧夂蚯闆r,11月下旬到次年2月下旬,冬小麥處于越冬期;2月底到3月上旬進(jìn)入返青期,隨著氣溫升高,冬小麥隨后開始快速、旺盛生長;3月中、下旬冬小麥開始起身,4月中旬到5月上旬為拔節(jié)孕穗期,隨后中旬進(jìn)入灌漿期,中、下旬為乳熟期,并在6月上旬完成收獲。

1.2 PyWOFOST作物生長模型及其所需數(shù)據(jù)

WOFOST作物生長模型以光攔截和CO2為驅(qū)動因素,應(yīng)用一系列參數(shù)模擬和控制不同類型、不同品種作物受環(huán)境條件影響時的呼吸、蒸騰、物候、干物質(zhì)積累及其分配等生長發(fā)育過程。其對作物生長發(fā)育控制的過程主要分為出苗、開花和成熟3個生育期,可估計特定作物類型每天的LAI、地上生物量和存儲器官總干質(zhì)量(TWSO,即糧食產(chǎn)量)等。另外,該模型可以在潛在模式(無水和養(yǎng)分脅迫引起的限制)、水分限制模式(土壤水分脅迫)或營養(yǎng)限制模式下運行。PyWOFOST模型是WOFOST模型的Python版本,其驅(qū)動數(shù)據(jù)包含氣象數(shù)據(jù)[9]、作物參數(shù)數(shù)據(jù)、土壤數(shù)據(jù)及田間管理數(shù)據(jù)。

1.2.1氣象數(shù)據(jù)

PyWOFOST模型需最低與最高溫度、太陽輻射、降雨量、風(fēng)速等氣象數(shù)據(jù)。研究使用田間建立的小型氣象觀測站數(shù)據(jù),如圖2所示。該氣象站點數(shù)據(jù)包含的逐日氣象要素包含平均氣溫、最高與最低氣溫、風(fēng)速、降雨、日照時數(shù)、水汽壓和輻射等?？紤]到試驗區(qū)范圍、平坦的地形因素及僅有的單個氣象站數(shù)據(jù),估產(chǎn)過程中應(yīng)用氣象站數(shù)據(jù)時未進(jìn)行反距離插值的面域處理。

圖2 氣象站位置及其實物圖Fig.2 Weather station location and weather station

1.2.2作物生長參數(shù)

冬小麥作物生長參數(shù)是PyWOFOST模型中的最重要參數(shù)[10],也是不確定性最大的參數(shù),描述小麥的遺傳特性、發(fā)育性狀和產(chǎn)量性狀等,控制著冬小麥的生長發(fā)育過程,與植株形態(tài)的發(fā)展和籽粒產(chǎn)量的形成直接相關(guān)。該參數(shù)包含不同發(fā)育階段所需要的積溫、光周影響因子、不同生育期的最大光合速率、不同生育期的比葉面積、干物質(zhì)分配系數(shù)、干物質(zhì)和葉片的死亡率等。在應(yīng)用PyWOFOST模型前,需要對冬小麥的作物參數(shù)數(shù)據(jù)進(jìn)行設(shè)置和標(biāo)定,以實現(xiàn)模型的本地化應(yīng)用。作物生長參數(shù)[10-12]見表1。

表1 PyWOFOST模型冬小麥主要作物參數(shù)Tab.1 Main crop parameters in PyWOFOST model

1.2.3土壤數(shù)據(jù)

土壤是作物生長的基礎(chǔ),土壤參數(shù)的準(zhǔn)確性直接影響作物模型的模擬精度。土壤質(zhì)地和結(jié)構(gòu)是土壤的基本物理特性。在冬小麥播種前,對土壤樣品進(jìn)行實地采樣并進(jìn)行理化性質(zhì)分析。測定0～10 cm、10～20 cm、20～50 cm、50～80 cm、80～120 cm、120～160 cm的6個土壤剖面參數(shù)。土壤參數(shù)主要包含土壤質(zhì)量含水率、土壤相對濕度、土壤容重、田間持水量、凋萎濕度等PyWOFOST模型需要輸入的土壤參數(shù)。參數(shù)如表2所示。

表2 PyWOFOST模型土壤參數(shù)Tab.2 Main soil parameters in PyWOFOST model

1.2.4田間管理數(shù)據(jù)

在冬小麥播種期前,實地調(diào)查土壤墑情,在10月22日開始播種,播種品種為鄭麥-119,播種量為210～270 kg/hm2。PyWOFOST模型的田間管理數(shù)據(jù)需要輸入灌溉時間,施肥元素、時間及設(shè)定利用效率。管理數(shù)據(jù)必須輸入作物種植時間、開始生長時間和收獲時間,以便和氣象數(shù)據(jù)進(jìn)行時間上的匹配。表3為田間管理和冬小麥實地考察的生育期。

表3 PyWOFOST模型田間管理和主要生育期參數(shù)Tab.3 PyWOFOST model management and primary reproductive period parameter

1.3 遙感數(shù)據(jù)與觀測數(shù)據(jù)

研究應(yīng)用GEE遙感數(shù)據(jù)平臺獲取Sentinel-2多光譜遙感數(shù)據(jù),該數(shù)據(jù)的具體介紹及獲取方式見文獻(xiàn)[13]。

結(jié)合冬小麥關(guān)鍵生育期、長勢狀況和衛(wèi)星過境時間確定葉面積指數(shù)、土壤含水率、地上生物量、葉片葉綠素含量、冠層光譜等采集工作。地面點數(shù)據(jù)采集時間為2022年3月4日、3月10日、4月2日、4月9日、4月19日、4月29日、5月23日、5月27日。冬小麥葉面積指數(shù)的單點獲取主要應(yīng)用LI-COR公司的LAI-2200型植被冠層分析儀,具體操作參照文獻(xiàn)[14]。采集點坐標(biāo)使用中海達(dá)V30型GPS儀器獲取,與遙感數(shù)據(jù)統(tǒng)一為WGS-84坐標(biāo)系。

對于產(chǎn)量數(shù)據(jù)的獲取,在冬小麥成熟期,依據(jù)不同長勢分布選取74個點作為樣方。在1 m×1 m的樣方范圍內(nèi),人工數(shù)出實際株數(shù)和有效麥穗數(shù),并對樣方范圍內(nèi)的小麥進(jìn)行收割。收割后對樣品進(jìn)行干燥、人工脫粒和稱量,并計算各樣方點產(chǎn)量。在種植過程中將地塊分為10個管理區(qū),為驗證區(qū)域產(chǎn)量估測的精度,在冬小麥?zhǔn)斋@過程中,對每個區(qū)的產(chǎn)量進(jìn)行區(qū)域級別的統(tǒng)計。采樣樣方設(shè)計、采樣點位置和管理區(qū)分布情況如圖3所示。

圖3 采樣樣方設(shè)計、采樣點位置和管理區(qū)分布Fig.3 Winter wheat sampling design, sampling point location and management area distribution

1.4 研究方法

1.4.1葉面積指數(shù)反演算法

不同物候期的LAI能夠反映當(dāng)時作物群體的長勢狀況及特征、光合作用、生物量和產(chǎn)量形成等諸多理化進(jìn)程[15]。研究應(yīng)用PROSAIL輻射傳輸模型反演不同時期的葉面積指數(shù)。具體PROSAIL輻射傳輸模型的參數(shù)敏感性分析、優(yōu)化參數(shù)選擇參照文獻(xiàn)[16],參數(shù)優(yōu)化方法選擇計算效率優(yōu)于查找表法[17]的SUBPLEX算法[18]。

1.4.2EnKF同化算法

EnKF的核心思想是利用集合預(yù)報的方式,基于統(tǒng)計得到的背景誤差協(xié)方差進(jìn)行狀態(tài)變量的最優(yōu)估計,實現(xiàn)對誤差協(xié)方差的更新。EnKF在處理各種數(shù)據(jù)不確定性時比較靈活,算法易于實現(xiàn)和操作,因此被廣泛應(yīng)用于地表數(shù)據(jù)同化研究中。EnKF具體算法及流程參照文獻(xiàn)[19]。

1.4.3RFR算法

應(yīng)用PyWOFOST作物生長模型反演單點產(chǎn)量作為機(jī)器學(xué)習(xí)算法訓(xùn)練數(shù)據(jù)時,該模型反演的產(chǎn)量誤差會對建立的機(jī)器學(xué)習(xí)模型產(chǎn)生影響。為減弱后續(xù)建立的機(jī)器學(xué)習(xí)模型再次產(chǎn)生的誤差傳遞而影響面域估產(chǎn)精度,研究選擇訓(xùn)練速度快、防過度擬合、抗噪性和模型泛化能力強(qiáng)等特點的RFR算法[20]。

1.4.4精度評價

采用決定系數(shù)R2、均方根誤差(Root mean square error,RMSE)、平均絕對誤差(Mean absolute error,MAE)、偏差(Bias)對反演的LAI與單點產(chǎn)量進(jìn)行反演精度評價。對于區(qū)域產(chǎn)量精度評價,計算公式為

P=1-|Rpredict-Rstatics|/Rstatics

(1)

式中P——區(qū)域產(chǎn)量監(jiān)測精度

Rpredict——區(qū)域產(chǎn)量反演結(jié)果

Rstatics——區(qū)域產(chǎn)量統(tǒng)計結(jié)果

2 結(jié)果與分析

2.1 PROSAIL輻射傳輸模型反演葉面積指數(shù)

基于Sentinel-2光學(xué)數(shù)據(jù)的B2、B3、B4、B8波段,應(yīng)用SUBPLEX算法和PROSAIL輻射傳輸模型實現(xiàn)單點尺度LAI的反演。實地LAI單點采樣時間和光學(xué)數(shù)據(jù)的對應(yīng)時間如表4所示。圖4為研究區(qū)不同時期的遙感影像,其時間與實地采樣時間存在差異。實地考察中,冬小麥在4月20日進(jìn)入揚花期,葉片發(fā)育開始減緩,4月28日進(jìn)行全域灌溉,造成4月29日采集樣本較少(僅7個),因此5月2日影像反演的LAI采用4月19日和29日的合并采集數(shù)據(jù)進(jìn)行驗證。研究基于單點LAI的采樣坐標(biāo),應(yīng)用SNAP軟件的“Raster→Export→Extract pixel values”操作步驟提取單點Sentinel-2影像的不同波段光譜數(shù)據(jù)。

表4 Sentinel-2衛(wèi)星影像日期和地面單點LAI采樣日期Tab.4 Sentinel-2 satellite image date and LAI collection samples date

圖4 研究區(qū)不同時期遙感影像Fig.4 Study area remote sensing images at different periods

每期光學(xué)影像數(shù)據(jù)單點尺度的葉面積指數(shù)反演結(jié)果與精度如圖5和表5所示。從圖5與表5可以看出,針對單幅影像LAI的反演情況,4月2日的R2最高,達(dá)到0.744 5。其原因是衛(wèi)星當(dāng)天過境,地面點LAI當(dāng)天采集,兩者數(shù)據(jù)時間匹配良好。5月2日的R2相對3月8日、4月2日和4月7日3期的R2較低,為0.608 1;造成該日R2較低的原因是地面點采集時間與遙感影像時間差距過大。5月13—25日,冬小麥生育期從灌漿期步入成熟期,冬小麥葉片及穗的逐日特征變化迅速、明顯,造成5月17日測量與反演LAI的R2非常小,僅為0.055 1,并存在明顯的低估現(xiàn)象,同時RMSE、MAE、Bias也達(dá)到最大,分別為0.910 3 m2/m2、0.780 0 m2/m2、0.524 6。就整體LAI反演情況而言,該反演方法存在一定的低估現(xiàn)象,但整體反演LAI的R2、RMSE較好,達(dá)到0.863 2、0.719 3 m2/m2,特別是Bias達(dá)到0.090 8,高于蘇偉等[16]在參數(shù)標(biāo)定前應(yīng)用PROSAIL模型反演LAI的精度(R2=0.73,RMSE為0.79 m2/m2)。RMSE和MAE也低于4月7日、5月2日和5月17日,達(dá)到0.719 3 m2/m2和0.574 7 m2/m2。因此,該方法反演LAI的可靠性良好。

表5 單點葉面積指數(shù)反演精度評價Tab.5 Evaluation of retrieval accuracy of single-point LAI

圖5 SUBPLEX+PROSAIL算法單點葉面積指數(shù)反演Fig.5 SUBPLEX+PROSAIL algorithm retrieval single-points LAI

應(yīng)用上述方法進(jìn)行區(qū)域LAI反演,同樣以B2、B3、B4、B8波段作為輸入,逐像元進(jìn)行區(qū)域LAI反演,不同時期反演LAI的值域分布結(jié)果如圖6所示。有研究表明,NDVI與LAI有較強(qiáng)的相關(guān)性[21],為探究區(qū)域LAI反演分布的可靠性,引入?yún)^(qū)域NDVI圖像進(jìn)行分析。區(qū)域NDVI分布情況如圖7所示。從圖7看出,NDVI在4月2日前期能夠反映作物的長勢信息,但在4月2日到5月2日內(nèi)趨于飽和,反映作物長勢信息能力相比LAI較差,且NDVI在達(dá)到飽和后,LAI和NDVI的相關(guān)性較弱。但在5月17日,冬小麥趨于成熟期,NDVI值降低,LAI同樣能夠反映冬小麥長勢分布及情況。

圖6 PROSAIL+SUBPLEX算法反演LAI值域分布Fig.6 LAI value domain distribution map by PROSAIL + SUBPLEX algorithm retrieval

2.2 基于EnKF算法的同化策略與估產(chǎn)實現(xiàn)

2.2.1PyWOFOST作物生長模型參數(shù)敏感性分析

PyWOFOST模型參數(shù)眾多,為對每個參數(shù)或參數(shù)間相互作用影響輸出結(jié)果的程度進(jìn)行定量評價,需對不同參數(shù)及參數(shù)間耦合作用進(jìn)行全局敏感性分析。少數(shù)模型參數(shù)通常是模型輸出的大部分可變性的原因,而大多數(shù)其他參數(shù)可能只有很小的影響。為集中精力校準(zhǔn)一定生物學(xué)意義的少數(shù)敏感性參數(shù)和保證模型模擬效果,又適當(dāng)減少計算量,參照文獻(xiàn)[22],對表1中27個PyWOFOST模型主要參數(shù)及大部分研究中未加入蒸散速率修正因子(CFET)和土壤水分消耗作物群數(shù)量(DEPNR)參數(shù)進(jìn)行敏感性分析[11]。

研究利用Sobol全局敏感性分析算法[3]對目標(biāo)變量TWSO進(jìn)行參數(shù)敏感性分析。Sobol法的參數(shù)樣本數(shù)設(shè)置為2 000×29。研究主要應(yīng)用PyWOFOST作物生長模型的潛在模式,基于2021—2022年冬小麥生長期間的站點氣象數(shù)據(jù)實現(xiàn)參數(shù)的敏感性分析,結(jié)果如圖8所示。潛在模式下對TWSO敏感的前4個主要參數(shù)有AMAXTB1、TDWI、TSUMEM和CVO。

圖8 對TWSO相關(guān)的參數(shù)敏感性分析Fig.8 Parameters sensitivity analysis about TWSO

2.2.2同化變量選取與參數(shù)標(biāo)定方案建立

作物生長模型是一個基于過程的模型,LAI能夠描述作物在不同生長階段的特性,常作為作物模型與遙感數(shù)據(jù)耦合的同化變量[18]。故研究同化變量主要選擇LAI,在校準(zhǔn)LAI時需綜合考慮當(dāng)?shù)氐纳a(chǎn)條件及與LAI相關(guān)參數(shù)的敏感特性。

參考許偉等[19]參數(shù)標(biāo)定方法及為提高單點估產(chǎn)精度,研究同樣選取作物生育期最大葉面積指數(shù)(LAImax)作為分析對象,結(jié)合2021—2022年氣象站提供的氣象數(shù)據(jù),分析PyWOFOST作物模型與LAI相關(guān)的14個重要輸入?yún)?shù)(TBASEM、RGRLAI、SPAN、TBASE、TDWI、Q10、RMR、RMS、RML、RMO、 CVL、CVR、CVS、 CVO)敏感性,分析結(jié)果如圖9所示。本文選取TDWI、TBASE、CVS、CVL作為對LAImax敏感的待優(yōu)化參數(shù)。

圖9 對LAImax相關(guān)的參數(shù)敏感性分析Fig.9 Parameters sensitivity analysis about LAImax

對于參數(shù)標(biāo)定方法,本文利用反演的多時相LAI和粒子群優(yōu)化算法[23](Particle swarm optimization,PSO)優(yōu)化與LAImax相關(guān)的4個敏感參數(shù),并將優(yōu)化后的值輸入到同化LAI的EnKF+PyWOFOST估產(chǎn)系統(tǒng)中,以實現(xiàn)單點產(chǎn)量的估算。同化運行過程中PyWOFOST模型應(yīng)用潛在模式,且集合的成員數(shù)(N)對同化速率和結(jié)果有重要的影響,相關(guān)研究表明N≥30時就能表現(xiàn)出更高的同化精度[12]。同化變量LAI反演精度的標(biāo)準(zhǔn)差設(shè)為本值的10%,以對LAI進(jìn)行擾動設(shè)置[2],生成LAI的集合N=200。基于實測單點產(chǎn)量進(jìn)行驗證分析,結(jié)果如表6和圖10所示。預(yù)測的單點產(chǎn)量與實測產(chǎn)量的R2、RMSE、MAE、Bias分別為0.866 5、468.64 kg/hm2、385.70 kg/hm2和103.08。相對于黃健熙等[12]、謝毅等[7]早前的估產(chǎn)精度,該方法具備較好的單點估測精度。

表6 單點尺度同化估產(chǎn)精度評價Tab.6 Single points yield accuracy estimation

圖10 同化多時相LAI實現(xiàn)單點尺度估產(chǎn)Fig.10 Assimilation multi-temporal LAI to realize single points yield estimation

2.2.3區(qū)域產(chǎn)量估算模型構(gòu)建與分析

研究基于本文的參數(shù)標(biāo)定方法和同化多期LAI到EnKF+PyWOFOST估產(chǎn)系統(tǒng)中,逐像元計算試驗區(qū)每個像元的產(chǎn)量模擬值,其產(chǎn)量分布結(jié)果如表7和圖11a所示。研究區(qū)(309.23 hm2),僅同化LAI過程得到區(qū)域產(chǎn)量結(jié)果花費時間接近144 h,即每公頃算法運行用時約0.47 h。從表7發(fā)現(xiàn),應(yīng)用同化方法的區(qū)域整體估產(chǎn)精度為89.01%,其中管理區(qū)最低、最高區(qū)域估產(chǎn)精度分別為第5管理區(qū)的76.37%和第7管理區(qū)的95.53%。結(jié)合圖3管理區(qū)分布和圖6中反演LAI的分布情況,發(fā)現(xiàn)區(qū)域LAI分布越均勻,估產(chǎn)精度越高。此原因是LAI表達(dá)作物的群體特征,應(yīng)用空間分辨率10 m的遙感數(shù)據(jù),對冬小麥長勢分布的表達(dá)能力有限。若提升長勢分布復(fù)雜的區(qū)域估產(chǎn)精度,建議采用更高分辨率的無人機(jī)數(shù)據(jù)。

表7 區(qū)域統(tǒng)計產(chǎn)量與模擬產(chǎn)量Tab.7 Regional statistical yields and simulated yields

圖11 不同模型冬小麥區(qū)域估產(chǎn)圖Fig.11 Regional yield estimation map of winter wheat based on different models

應(yīng)用機(jī)器學(xué)習(xí)算法實現(xiàn)作物估產(chǎn)的依據(jù)是作物生長因素受多種因素的影響,且其生長過程是一個非常復(fù)雜的生物生理過程,但作物生長狀況、產(chǎn)量可以用一些與其生長過程密切相關(guān)的參數(shù)進(jìn)行表征。

因此機(jī)器學(xué)習(xí)方法是一種與作物生長模型不同的估產(chǎn)方法。由于冬小麥產(chǎn)量形成的非線性和非平穩(wěn)性,在應(yīng)用機(jī)器學(xué)習(xí)進(jìn)行冬小麥產(chǎn)量建模時,需結(jié)合多個生育期的關(guān)鍵影響因子綜合建模,才能有望提升冬小麥產(chǎn)量預(yù)測的精度。劉新杰等[24]利用冬小麥抽穗期和乳熟期的累積NDVI值可實現(xiàn)冬小麥產(chǎn)量的精確估算,黃健熙等[25]建議使用多時序的NDVI值進(jìn)行估產(chǎn)研究。研究參照上述研究結(jié)論,應(yīng)用多時相NDVI和其NDVI累計值作為RFR算法預(yù)測產(chǎn)量的變量。其反演區(qū)域產(chǎn)量結(jié)果如表8和圖11b所示。此過程花費時間為8 min,平均每公頃算法運行用時約1.55 s。從表8發(fā)現(xiàn),區(qū)域整體估產(chǎn)精度雖達(dá)到99.44%,不同管理區(qū)的估產(chǎn)精度仍存在較大差異,如管理區(qū)5和管理區(qū)7的估產(chǎn)精度與物理模型區(qū)域估產(chǎn)精度表現(xiàn)一致。

表8 基于機(jī)器學(xué)習(xí)的冬小麥區(qū)域估產(chǎn)統(tǒng)計Tab.8 Statistics of winter wheat regional yield estimation based on machine learning

為更快地實現(xiàn)區(qū)域估產(chǎn)及避免現(xiàn)場產(chǎn)量采集工作,將上述兩類算法進(jìn)行耦合以實現(xiàn)區(qū)域冬小麥估產(chǎn)研究。研究使用PyWOFOST物理模型反演采集點的單點產(chǎn)量,結(jié)合5期遙感數(shù)據(jù)的NDVI和NDVI累計值作為訓(xùn)練數(shù)據(jù),輸入RFR模型中建立估產(chǎn)模型,得出區(qū)域尺度的產(chǎn)量,結(jié)果如表9和圖11c所示。此過程花費時間44.25 min,平均每公頃用時約8.85 s,其中利用EnKF+PyWOFOST模型花費時長約36 min,占總時長約81.35%。耦合模型的區(qū)域估產(chǎn)精度高于物理模型6.57個百分點,低于機(jī)器學(xué)習(xí)模型3.86個百分點,且在各管理區(qū)域的估產(chǎn)精度分布上,與上述兩者表現(xiàn)出較強(qiáng)的一致性,說明耦合模型具有較強(qiáng)的適用性。

表9 耦合PyWOFOST作物生長模型與機(jī)器學(xué)習(xí)的冬小麥區(qū)域估產(chǎn)統(tǒng)計Tab.9 Regional winter wheat yield estimation coupled PyWOFOST crop growth model with machine learning

3 討論

研究使用PROSAIL輻射傳輸模型結(jié)合Sentinel-2光學(xué)數(shù)據(jù)反演不同時期的LAI,相對于應(yīng)用機(jī)理性弱的機(jī)器學(xué)習(xí)方法反演LAI,該方法具有普適性[26]。PROSAIL模型參數(shù)眾多,特別是采用查找表法計算參數(shù)時,容易造成病態(tài)反演問題,導(dǎo)致反演LAI的可靠性減弱。后續(xù)研究將利用該模型結(jié)合MCMC[27]、VMG[26]等參數(shù)優(yōu)化方法反演LAI,檢驗SUBPLEX+PROSAIL算法組合反演LAI的可靠性和抗噪能力。

同化多期Sentinel-2光學(xué)數(shù)據(jù)反演的變量時,光學(xué)成像時間基本固定,易造成同化時影像數(shù)據(jù)時間過于集中、分散或關(guān)鍵物候期缺乏[28]的情況,對物理模型估算單點產(chǎn)量精度造成一定影響,近而影響耦合模型的區(qū)域估產(chǎn)精度。后期研究中使用不受天氣影響Sentinel-1 SAR數(shù)據(jù)反演同化所需的變量[21]或利用更加豐富的多源遙感數(shù)據(jù)增加同化過程中的數(shù)據(jù)量,以提高估產(chǎn)精度[29]。當(dāng)然,過多的數(shù)據(jù)量會降低計算效率,后續(xù)研究需確定合適同化數(shù)據(jù)量、同化變量的時間與同化步長,在同化數(shù)據(jù)量和估產(chǎn)精度方面做均衡評估。

耦合作物生長模型和機(jī)器學(xué)習(xí)模型存在誤差傳遞問題,即基于作物生長模型估產(chǎn)的單點產(chǎn)量本身存在一定的誤差和機(jī)器學(xué)習(xí)是高度依賴訓(xùn)練樣本的算法,單點產(chǎn)量用于機(jī)器學(xué)習(xí)建模中,會導(dǎo)致誤差的傳遞,而如何解決和描述誤差傳遞的過程和不確定性分析是本文存在的缺陷和難題。研究中僅使用PyWOFOST作物生長模型。作物生長過程是一個非常復(fù)雜的非線性系統(tǒng),具備高度的空間異質(zhì)性,單一的作物模型無法完全表達(dá)和模擬作物生長過程中的所有動態(tài)變量,因此需建立多作物模型估產(chǎn)系統(tǒng)進(jìn)行產(chǎn)量評估,如DSSAT[30]、APSIM-Wheat[31]等模型,依據(jù)一定的權(quán)重實現(xiàn)多模型聯(lián)合產(chǎn)量估算。

當(dāng)前的研究多基于MLR、RFR或SVR方法來探究多時相NDVI、LAI與作物產(chǎn)量的關(guān)系[9]。雖然使用上述幾種方法能夠產(chǎn)生更好的產(chǎn)量估算,但其算法本身不能處理衛(wèi)星的時間序列數(shù)據(jù),致使作物某階段的生長狀態(tài)變量會影響最終產(chǎn)量上體現(xiàn)較弱。相關(guān)研究表明,LSTM模型能夠從時間序列中學(xué)習(xí)歷史特征[32],比從靜態(tài)角度處理輸入數(shù)據(jù)的方法(如RFR方法)具有更大的潛力,若將CNN算法[33]與LSTM模型進(jìn)行融合[34],則模型可從衛(wèi)星的時間序列上提取空間和時間特征,從而提升冬小麥的估產(chǎn)精度。

4 結(jié)論

(1)基于觀測的時序LAI,應(yīng)用PSO優(yōu)化算法逐像元優(yōu)化與LAImax相關(guān)的敏感性參數(shù)(TDWI、TBASE、CVS、CVL),將優(yōu)化參數(shù)輸入到PyWOFOST作物生長模型中,再使用EnKF算法同化時序LAI到PyWOFOST作物生長模型,以調(diào)整對TWSO敏感的4個參數(shù)(AMAXTB1、TDWI、TSUMEM、CVO),得到的單點模擬產(chǎn)量具有較高精度,其R2、RMSE分別達(dá)到0.866 5和468.64 kg/hm2,優(yōu)于相關(guān)研究。

(2)為解決同化遙感信息與作物生長模型在區(qū)域估產(chǎn)中計算低效的問題,應(yīng)用EnKF算法同化PROSAIL模型反演的多期LAI到PyWOTOST作物生長模型中,估測一定數(shù)量冬小麥不同長勢采樣點的產(chǎn)量,將模擬的多個單點產(chǎn)量和該點的多時相NDVI數(shù)據(jù)建立RFR機(jī)器學(xué)習(xí)算法模型,并利用面域多時相NDVI和建立的RFR模型實現(xiàn)區(qū)域冬小麥估產(chǎn)。結(jié)果表明,耦合模型區(qū)域產(chǎn)量模擬精度在試驗區(qū)的小區(qū)域范圍內(nèi)估產(chǎn)精度最低為81.28%,在試驗區(qū)總體估產(chǎn)精度達(dá)到95.58%,有效避免機(jī)器學(xué)習(xí)所需的人工實地采樣和減弱利用同化方法實現(xiàn)區(qū)域估產(chǎn)面臨的巨大算力需求,為快速、準(zhǔn)確實現(xiàn)區(qū)域作物估產(chǎn)提供新的思路。