融合多源時空數(shù)據(jù)的冬小麥產(chǎn)量預(yù)測模型研究

王來剛 鄭國清 郭 燕 賀 佳 程永政

(1.河南省農(nóng)業(yè)科學(xué)院農(nóng)業(yè)經(jīng)濟與信息研究所， 鄭州 450002； 2.農(nóng)作物種植監(jiān)測與預(yù)警河南省工程實驗室， 鄭州 450002)

0 引言

及時準(zhǔn)確掌握區(qū)域農(nóng)作物產(chǎn)量信息，能夠為糧食生產(chǎn)宏觀調(diào)控、經(jīng)濟政策制定和農(nóng)作物保險提供決策支持，對服務(wù)國家糧食安全戰(zhàn)略具有重要意義。小麥?zhǔn)俏覈Z食作物之一，也是我國重要的口糧作物，對糧食安全和居民生活都具有重要意義。

冬小麥產(chǎn)量形成是一個復(fù)雜的過程，中間涉及非常多的生理生化過程，氣象條件、土壤條件、地理環(huán)境、物候信息等都是其產(chǎn)量估算需要考慮的因素，產(chǎn)量可看成是一段時期內(nèi)多個影響因子相互疊加的結(jié)果。遙感技術(shù)在農(nóng)業(yè)中已廣泛應(yīng)用，可以從遙感數(shù)據(jù)中提取出各種相關(guān)信息用于產(chǎn)量預(yù)測。特別是植被指數(shù)，如歸一化差異植被指數(shù)(Normalized difference vegetation index， NDVI)，已被廣泛使用[1-3]。其他指數(shù)如增強植被指數(shù)(Enhanced vegetation index，EVI)[4]、作物水分脅迫指數(shù)[5]、綠色植被指數(shù)、土壤調(diào)節(jié)植被指數(shù)[6]等，也被用于作物產(chǎn)量預(yù)測。氣象數(shù)據(jù)(如降水量、空氣溫度等)[7-9]和土壤條件數(shù)據(jù)(如土壤含水率和溫度)[10]，在產(chǎn)量預(yù)測中經(jīng)常被用作作物生長環(huán)境指標(biāo)。日光誘導(dǎo)葉綠素?zé)晒膺b感數(shù)據(jù)(SIF)為植物進行光吸收后葉綠素a 釋放出的微弱電磁輻射信號[11]，與光合速率緊密耦合。近年來，SIF數(shù)據(jù)日益增多，為作物估產(chǎn)提供了一條新的有效途徑[12-14]?；谶b感數(shù)據(jù)的作物產(chǎn)量預(yù)測模型主要有2種：作物模擬模型和經(jīng)驗統(tǒng)計模型。盡管作物模擬模型精確地模擬了作物生長的物理過程，但需要精細的土壤和氣象數(shù)據(jù)，這阻礙了它們的大規(guī)模應(yīng)用。相比之下，經(jīng)驗統(tǒng)計模型簡單，需要較少的輸入數(shù)據(jù)，因此被廣泛用作基于過程的模型的常見替代方案。上述相關(guān)研究多集中在利用植被指數(shù)和氣象數(shù)據(jù)預(yù)測作物產(chǎn)量，將SIF數(shù)據(jù)作為產(chǎn)量預(yù)測因子的研究相對較少。

隨機森林回歸模型是一種基于機器學(xué)習(xí)的統(tǒng)計方法，它能處理高維度的數(shù)據(jù)并且估算結(jié)果具有較高的準(zhǔn)確率。楊北萍等[15]將氣象數(shù)據(jù)和遙感數(shù)據(jù)作為輸入變量，基于隨機森林回歸方法預(yù)測了水稻產(chǎn)量。王鵬新等[16]通過隨機森林回歸算法獲取玉米主要生育時期各個特征變量的權(quán)重，構(gòu)建了加權(quán)特征變量與玉米單產(chǎn)間的回歸模型。程千等[17]采用偏最小二乘回歸、支持向量機回歸和隨機森林回歸3種機器學(xué)習(xí)算法進行冬小麥產(chǎn)量估測，其中隨機森林回歸模型估測精度最高。本文以河南省為研究區(qū)域，擬通過探索分析多年小麥生長季的遙感數(shù)據(jù)、氣象數(shù)據(jù)、土壤含水率數(shù)據(jù)、地形數(shù)據(jù)與小麥單產(chǎn)的相關(guān)性，利用隨機森林回歸模型融合多源時空數(shù)據(jù)構(gòu)建冬小麥產(chǎn)量預(yù)測模型，以期提高大尺度冬小麥產(chǎn)量預(yù)測精度。

1 材料與方法

1.1 研究區(qū)域概況

本文將河南省作為研究區(qū)域。河南省位于我國中部的黃河中下游地區(qū)，地處我國第二階梯向第三階梯的過渡地帶，北、西、南三面的太行山、伏牛山、桐柏山、大別山沿省界呈半環(huán)形分布；中、東部為黃淮海沖積平原；西南部為南陽盆地。全省耕地面積8.110 3×106hm2，占中國總耕地面積的6.24%，是我國小麥的最適生態(tài)區(qū)，小麥?zhǔn)呛幽鲜∽钪饕募Z食作物，全省小麥播種面積、總產(chǎn)量均居全國首位。

1.2 數(shù)據(jù)獲取與處理

研究需要的數(shù)據(jù)包括2005—2019年河南省小麥生長季的遙感數(shù)據(jù)、氣象數(shù)據(jù)、土壤含水率數(shù)據(jù)、數(shù)字高程和小麥單產(chǎn)數(shù)據(jù)，所有數(shù)據(jù)利用小麥種植空間分布數(shù)據(jù)進行掩膜，并疊加河南省縣級行政邊界，經(jīng)柵格均值計算，得各縣區(qū)的變量均值。

1.2.1遙感數(shù)據(jù)

衛(wèi)星遙感參數(shù)分別選擇增強型植被指數(shù)(EVI)和日光誘導(dǎo)葉綠素?zé)晒?SIF)來反映小麥長勢和光合作用。EVI是一個廣泛用于監(jiān)測作物長勢和預(yù)測產(chǎn)量的植被指數(shù)[18]，與NDVI相比，EVI對較高的冠層葉面積指數(shù)敏感，受大氣氣溶膠影響較小。本文使用了美國國家航空航天局(NASA)MODIS EVI 產(chǎn)品(MOD13A1, Collection 6)，該產(chǎn)品的空間分辨率為500 m，時間分辨率為16 d。近年來，隨著技術(shù)的發(fā)展，SIF遙感數(shù)據(jù)不斷增多、空間分辨率不斷提高，使得利用 SIF 遙感數(shù)據(jù)直接預(yù)測作物產(chǎn)量成為可能[19-20]。SIF數(shù)據(jù)是從全球空間連續(xù)日光誘導(dǎo)葉綠素?zé)晒鈹?shù)據(jù)集產(chǎn)品(CSIF)獲得(DOI:10.6084/m9.figshare.6387494)，其時間分辨率為4 d和全球空間分辨率為0.05°。應(yīng)用最大值合成技術(shù)生成冬小麥生長季(10月至次年5月)月時間尺度的EVI和SIF。

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

從Terra Climate數(shù)據(jù)集獲得小麥整個生育期每月的降水量、最高溫度、最低溫度數(shù)據(jù)作為影響產(chǎn)量形成的特征變量。這些數(shù)據(jù)集對1958—2019年間的全球地表具有高空間分辨率(約4 km)，可用于區(qū)域作物產(chǎn)量預(yù)測[21-22]。

1.2.3土壤含水率數(shù)據(jù)

土壤含水率數(shù)據(jù)來自國家青藏高原科學(xué)數(shù)據(jù)中心的中國土壤水分?jǐn)?shù)據(jù)集(DOI: 10.5281/zenodo.4738556)，時間分辨率為月，空間分辨率為0.05°。

1.2.4數(shù)字高程數(shù)據(jù)

來源于美國國家航空航天局(NASA)ASTERG-DEM，柵格尺寸為30 m×30 m。

1.2.5冬小麥單產(chǎn)數(shù)據(jù)

河南省縣級冬小麥單產(chǎn)數(shù)據(jù)來源于《河南省統(tǒng)計年鑒》，2005—2019年間具有有效小麥單產(chǎn)數(shù)據(jù)的縣市105個。

1.3 研究方法

1.3.1隨機森林回歸模型

隨機森林是由多棵分類回歸樹(Classification and regression tree，CART)構(gòu)成的組合分類模型[23]。該研究將各年份各縣市的小麥產(chǎn)量因子特征變量數(shù)據(jù)和產(chǎn)量數(shù)據(jù)進行集成，共同構(gòu)成隨機森林的樣本數(shù)據(jù)集，通過自助法(Bootstrap)從原始樣本集采樣得到構(gòu)建N棵樹所需的N個子集，每次未被抽到的數(shù)據(jù)稱為袋外數(shù)據(jù)(Out-of-bag，OOB)，用來進行內(nèi)部誤差估計和變量重要性評價；生成每棵樹時，從規(guī)模為M的特征變量集中隨機選擇m個變量(m

圖1 基于隨機森林算法的冬小麥產(chǎn)量預(yù)測流程圖Fig.1 Flow chart of winter wheat yield prediction based on random forest algorithm

1.3.2特征變量重要性評價

隨機森林算法是眾多決策樹并行式集成學(xué)習(xí)的方法，除了可對數(shù)據(jù)集進行分類和回歸外，還可以進行變量重要性分析、奇異值檢測等[24]，能夠解釋若干自變量對因變量的作用。通過模型內(nèi)部重要性評價結(jié)果，分析不同影響因子對小麥產(chǎn)量的影響程度。這種方法是直接分析評價每種特征變量對模型預(yù)測準(zhǔn)確率的影響，基本思想是重新排列特征變量的順序，觀測模型準(zhǔn)確率的降低程度。對于不重要的特征變量，這種方法對模型準(zhǔn)確率的影響很小，但是對于重要特征變量卻會極大降低模型的準(zhǔn)確率。具體過程為：

(1)對于隨機森林中的每一棵決策樹，使用相應(yīng)的袋外數(shù)據(jù)(OOB)來計算它的袋外數(shù)據(jù)誤差，記為eOOB1。

(2)隨機地對OOB數(shù)據(jù)所有樣本的特征變量加入噪聲干擾，隨機改變樣本在特征變量處的值，再次計算它的袋外數(shù)據(jù)誤差，記為eOOB2。

(3)假設(shè)隨機森林中有Ntree棵樹，則特征變量的重要性指標(biāo)為∑(eOOB2-eOOB1)/Ntree。之所以可以用這個表達式作為相應(yīng)特征的重要性度量值，是因為：若給某個特征隨機加入噪聲之后，袋外的準(zhǔn)確率大幅度降低，則說明這個特征對于樣本的分類結(jié)果影響很大，即重要程度比較高。

1.3.3模型精度評價

采用決定系數(shù)(Coefficient of determination，R2)、均方根誤差(Root mean square error，RMSE)和相對誤差(Relative error，RE)3個指標(biāo)評價模型精度。其中R2評價預(yù)測模型擬合能力，RMSE和RE評價預(yù)測值和實測值離散程度。

2 結(jié)果與分析

2.1 特征變量與冬小麥產(chǎn)量時空相關(guān)性分析

為研究特征變量與小麥產(chǎn)量之間的關(guān)系，分別從時間和空間上進行了各特征變量數(shù)據(jù)與產(chǎn)量之間相關(guān)性分析，見圖2，箱線圖是特征變量與產(chǎn)量相關(guān)性的時間模式，相關(guān)性的空間模式基于具有最高相關(guān)系數(shù)的月份，即方框圖中的紅點。在時間相關(guān)性模式上，遙感參數(shù)EVI和SIF與產(chǎn)量都呈正相關(guān)，其中4月相關(guān)性達到最高峰，相關(guān)系數(shù)可達0.7左右，隨后由于小麥冠層衰老和籽粒形成而開始下降，使得遙感參數(shù)與產(chǎn)量之間的相關(guān)性降低。與遙感參數(shù)不同，土壤含水率與產(chǎn)量之間的最大相關(guān)系數(shù)出現(xiàn)時期并不突出，主要在10月和次年3、4月相關(guān)性比較大，但整體相關(guān)性低于遙感參數(shù)。氣象變量中，各生長階段降水量和最低溫度與小麥氣象波動產(chǎn)量為正相關(guān)，降水量與產(chǎn)量之間的相關(guān)性最大的月份與土壤含水率基本相同，最低溫度與產(chǎn)量相關(guān)性最大的月份為3月。3月小麥主要處于拔節(jié)-孕穗階段，最低溫度是影響小麥產(chǎn)量的重要因素。最高溫度與產(chǎn)量主要呈負相關(guān)性，但相關(guān)系數(shù)較小。

圖2 特征變量與冬小麥產(chǎn)量的時空相關(guān)性Fig.2 Spatiotemporal correlations between characteristic variable and wheat yield

選擇特征變量與小麥產(chǎn)量相關(guān)性最高的月份，分析特征變量與小麥產(chǎn)量相關(guān)性的空間分布特征。遙感參數(shù)EVI和SIF在豫北平原和豫東平原與產(chǎn)量相關(guān)性較高，相關(guān)系數(shù)最高達0.8。豫西豫南山地丘陵地區(qū)和部分種植模式混雜的縣區(qū)相關(guān)性較低，相關(guān)系數(shù)大部分在0～0.2之間，這些差異可能是由于豫西、豫南山地丘陵地區(qū)地塊破碎，EVI和SIF混合像元較多，遙感參數(shù)難以真實反映小麥長勢情況。土壤含水率與產(chǎn)量相關(guān)性在空間分布上，無明顯特征。降水量在灌溉條件較差的丘陵地區(qū)與產(chǎn)量相關(guān)性較高，豫北地區(qū)相關(guān)性略低，可能由于豫北地區(qū)灌溉條件較好，小麥生長對降水量依賴性不高。最高溫度在全省與小麥產(chǎn)量主要呈負相關(guān)，在山地丘陵地區(qū)相關(guān)性較高。最低溫度在豫北地區(qū)相關(guān)性較高，豫北地區(qū)冬小麥易遭受晚霜凍害，最低溫度是影響產(chǎn)量的一項重要制約因素。

2.2 特征變量重要性分析

冬小麥特征變量可以分為動態(tài)特征變量(遙感和氣象變量等)和靜態(tài)特征變量(高程)。利用2005—2019年的冬小麥特征動態(tài)變量數(shù)據(jù)與產(chǎn)量分別進行隨機森林的OOB重要性分析，按照模型準(zhǔn)確率降低的程度對特征變量由大到小排序(表1)。4月EVI和SIF、2月降水量排名位居前3位。開花期是冬小麥體內(nèi)新陳代謝最旺盛的生長時期, 正是小麥產(chǎn)量形成的關(guān)鍵時期。河南省冬小麥的開花期主要集中在4月，EVI和SIF與冬小麥葉面積指數(shù)、生物量等苗情指標(biāo)具有較強的相關(guān)性，因此在各特征變量中重要性更高，重要性指標(biāo)在0.4左右。2月是河南省冬小麥返青期，是促進苗情轉(zhuǎn)化升級的關(guān)鍵時期，降水量是一個重要的影響因素。重要性指標(biāo)第4～6位分別為10月土壤含水率、4月最低溫度和2月EVI。10月是冬小麥播種出苗時期，土壤墑情是出苗品質(zhì)的關(guān)鍵因子。4月初小麥處于孕穗期，很容易受到低溫凍害，特別是豫北地區(qū)。2月EVI代表了小麥生長前期的總體苗情，也是影響產(chǎn)量的重要因子。

表1 動態(tài)特征變量重要性指標(biāo)Tab.1 Importance of dynamic feature variables

對EVI、SIF、土壤含水率、降水量、最高溫度、最低溫度和高程7大類特征變量進行基于隨機森林的袋外OOB重要性分析(圖3)。結(jié)果表明，EVI、SIF和高程對小麥產(chǎn)量的重要性遠大于其他因素，重要性指標(biāo)均超過0.45，說明除了遙感參數(shù)外，地形對產(chǎn)量的影響也非常大。其次是土壤含水率和降水量，說明含水率是影響小麥產(chǎn)量的重要環(huán)境因子。最后是最高溫度和最低溫度，重要性指標(biāo)只有0.10左右。綜上，描述作物生長狀況EVI和SIF對產(chǎn)量預(yù)測的貢獻大于描述作物生長環(huán)境的氣象和土壤因子。

圖3 影響因子重要性統(tǒng)計圖Fig.3 Importance of impact factors

2.3 不同時間段模型預(yù)測精度比較

為了深入研究不同時間段的冬小麥單產(chǎn)預(yù)測精度，并結(jié)合單產(chǎn)預(yù)測實際工作需求，分別輸入10月—次年2月、10月—次年3月、10月—次年4月、10月—次年5月的特征變量，對小麥產(chǎn)量進行隨機森林建模，結(jié)果見圖4。由圖4可知，基于小麥整個生長階段10月—次年5月和10月—次年4月特征變量的估產(chǎn)精度較高，R2分別為0.85和0.84，RMSE均最小，2個時間段模型預(yù)測精度基本相同。河南省冬小麥抽穗期和開花期主要集中在4月，開花期是冬小麥體內(nèi)新陳代謝最旺盛的生長時期，正是小麥產(chǎn)量形成的關(guān)鍵時期，溫度和降水量與小麥產(chǎn)量密切相關(guān)[25]。僅有小麥生長前期的10月—次年2月和10月—次年3月特征變量的估產(chǎn)精度相對偏低，RMSE均大于900 kg/hm2。因此，4月底或5月初是冬小麥產(chǎn)量預(yù)測的最佳時間。

圖4 不同時間段模型預(yù)測結(jié)果對比Fig.4 Comparisons of prediction results of different periods

2.4 單產(chǎn)預(yù)測誤差空間分布特征

輸入冬小麥各生長階段的EVI、SIF、土壤含水率、氣象要素和高程等全部變量，以小麥產(chǎn)量為目標(biāo)變量構(gòu)建隨機森林模型，預(yù)測2018、2019年河南省各縣冬小麥產(chǎn)量，圖5為冬小麥產(chǎn)量預(yù)測相對誤差與數(shù)字高程疊加空間分布情況。整體上，2年模型預(yù)測相對誤差空間分布基本相似，產(chǎn)量預(yù)測相對誤差均在20%以內(nèi)，平原地區(qū)模型預(yù)測相對誤差大部分在10%以內(nèi)，而豫西和豫南丘陵山地模型預(yù)測相對誤差高于平原地區(qū)，這可能由于丘陵山地區(qū)域地塊破碎，而模型輸入數(shù)據(jù)EVI、SIF等數(shù)據(jù)空間分辨率較低，混合像元現(xiàn)象比較嚴(yán)重，很難準(zhǔn)確反映小麥實際長勢情況。另外，信陽地區(qū)和開封部分縣小麥種植地塊零星破碎，預(yù)測精度也相對偏低。2年數(shù)據(jù)對比可以看出，2019年預(yù)測相對誤差低于2018年，特別在豫北地區(qū)，2018年預(yù)測相對誤差明顯較大。這可能由于2018年4月初河南豫北地區(qū)發(fā)生一次較重的晚霜凍害，不同小麥品種抗凍害能力差別較大，產(chǎn)量損失程度更是不一致，而該預(yù)測模型無品種特征變量。

圖5 河南省冬小麥產(chǎn)量預(yù)測相對誤差空間分布Fig.5 Relative error spatial distributions of winter wheat yield prediction in Henan Province

3 討論

研究表明，SIF與小麥產(chǎn)量呈正相關(guān)，這一結(jié)果得到了已有研究的支持[26-27]。研究區(qū)內(nèi)，SIF和EVI與小麥產(chǎn)量的相關(guān)性相近，EVI對綠葉結(jié)構(gòu)、葉綠素含量和生物量的變化比較敏感，而SIF有著直接指示植被光合作用的巨大潛力，2種數(shù)據(jù)可能共同反映一些與地上作物生物量相關(guān)的信息。本文研究的小麥產(chǎn)量預(yù)測模型輸入變量重點考慮了氣象條件、植被生長狀態(tài)、地形地貌等，但在實際生產(chǎn)過程中小麥產(chǎn)量會受到多種環(huán)境因素的影響，以及各種人為因素造成的產(chǎn)量波動，這些都會影響最終的預(yù)測精度。氣象條件中雖然用到了最高溫度和最低溫度等特征變量，但顯然預(yù)測模型對災(zāi)害信息反應(yīng)不敏感，在對2018年冬小麥產(chǎn)量預(yù)測中，受到晚霜凍害的區(qū)域小麥預(yù)測精度偏低。根據(jù)實際調(diào)查，小麥凍害除了與最低溫度有關(guān)之外，與品種、種植密度和土壤特性都有明顯的相關(guān)性。因此，未來的研究中可以增加災(zāi)害分布圖層作為輸入變量，從而提高模型預(yù)測精度。

隨機森林算法在預(yù)測產(chǎn)量上具有很大的潛力，但并不意味著該算法一定是預(yù)測產(chǎn)量的最優(yōu)算法，后期研究將對比深度神經(jīng)網(wǎng)絡(luò)、卷積神經(jīng)網(wǎng)絡(luò)等深度學(xué)習(xí)算法在作物預(yù)測產(chǎn)量的優(yōu)缺點，解決算法對于極端天氣等突發(fā)事件的適應(yīng)性[28-29]。這項研究受到了一些不確定性的困擾，其中一個問題是SIF和EVI的低信噪比和粗糙的時空分辨率[30-31]，可能無法真實反映小麥長勢空間特征的細節(jié)。特別是在山地丘陵地塊破碎的地區(qū)，產(chǎn)量預(yù)測精度都有不同程度下降。因此，如果要提高遙感數(shù)據(jù)在作物產(chǎn)量預(yù)測中的作用，時空連續(xù)高分辨率數(shù)據(jù)是非常重要的。隨著衛(wèi)星遙感技術(shù)發(fā)展，將實現(xiàn)具有更高空間和時間分辨率的EVI和SIF。例如，歐空局設(shè)計的FLEX(Fluorescence explorer)，提供全球尺度上的植被葉綠素?zé)晒鉁y量衛(wèi)星數(shù)據(jù)，目前其空間分辨率設(shè)定為300 m×300 m，軌道幅寬為150 km，非常有利于小尺度精準(zhǔn)化的作物估產(chǎn)研究[32]。

4 結(jié)論

(1)以遙感數(shù)據(jù)、氣象數(shù)據(jù)、土壤含水率數(shù)據(jù)和地形數(shù)據(jù)為特征變量，分析了2005—2019年小麥產(chǎn)量與特征變量的相關(guān)性時空特征，基于隨機森林算法對特征變量進行了重要性分析，并建立了不同時間段的小麥產(chǎn)量預(yù)測模型。

(2)SIF是指示植被光合變化的有效探針，與小麥產(chǎn)量呈高度正相關(guān)，在小麥產(chǎn)量預(yù)測特征變量重要性分析中，與EVI起著同等重要作用，是小麥產(chǎn)量預(yù)測的一項重要因子。

(3)在7大類特征變量中，EVI、SIF和高程對小麥產(chǎn)量的重要性遠大于其他因素，重要性指標(biāo)都超過0.45，其次是土壤含水率和降水量，最后是最高溫度和最低溫度。高程為靜態(tài)特征變量，在空間有差異，年度間無變化。

(4)基于隨機森林算法以小麥生長階段10月—次年5月和10月—次年4月為特征變量構(gòu)建的產(chǎn)量預(yù)測模型精度較高，R2分別為0.85和0.84，RMSE分別為821.55、832.01 kg/hm2。在產(chǎn)量預(yù)測誤差空間分布特征上，全省相對誤差均在20%以內(nèi)，平原地區(qū)模型預(yù)測相對誤差大部分在10%以內(nèi)，而豫西和豫南丘陵山地模型預(yù)測相對誤差高于平原地區(qū)。