融合多源時空數據的冬小麥產量預測模型研究

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

(1.河南省農業科學院農業經濟與信息研究所， 鄭州 450002； 2.農作物種植監測與預警河南省工程實驗室， 鄭州 450002)

0 引言

及時準確掌握區域農作物產量信息，能夠為糧食生產宏觀調控、經濟政策制定和農作物保險提供決策支持，對服務國家糧食安全戰略具有重要意義。小麥是我國糧食作物之一，也是我國重要的口糧作物，對糧食安全和居民生活都具有重要意義。

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

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

1 材料與方法

1.1 研究區域概況

本文將河南省作為研究區域。河南省位于我國中部的黃河中下游地區，地處我國第二階梯向第三階梯的過渡地帶，北、西、南三面的太行山、伏牛山、桐柏山、大別山沿省界呈半環形分布；中、東部為黃淮海沖積平原；西南部為南陽盆地。全省耕地面積8.110 3×106hm2，占中國總耕地面積的6.24%，是我國小麥的最適生態區，小麥是河南省最主要的糧食作物，全省小麥播種面積、總產量均居全國首位。

1.2 數據獲取與處理

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

1.2.1遙感數據

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

1.2.2氣象數據

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

1.2.3土壤含水率數據

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

1.2.4數字高程數據

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

1.2.5冬小麥單產數據

河南省縣級冬小麥單產數據來源于《河南省統計年鑒》，2005—2019年間具有有效小麥單產數據的縣市105個。

1.3 研究方法

1.3.1隨機森林回歸模型

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

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

1.3.2特征變量重要性評價

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

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

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

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

1.3.3模型精度評價

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

2 結果與分析

2.1 特征變量與冬小麥產量時空相關性分析

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

圖2 特征變量與冬小麥產量的時空相關性Fig.2 Spatiotemporal correlations between characteristic variable and wheat yield

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

2.2 特征變量重要性分析

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

表1 動態特征變量重要性指標Tab.1 Importance of dynamic feature variables

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

圖3 影響因子重要性統計圖Fig.3 Importance of impact factors

2.3 不同時間段模型預測精度比較

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

圖4 不同時間段模型預測結果對比Fig.4 Comparisons of prediction results of different periods

2.4 單產預測誤差空間分布特征

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

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

3 討論

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

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

4 結論

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

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

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

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