基于多源時空數據的冬小麥產量預測模型

劉靖宇, 李遠斌, 董 昊, 馮俊辰, 丁云鴻

(哈爾濱師范大學計算機科學與信息工程學院,黑龍江哈爾濱 150025)

由于人口的增長和人們生活水平的提高,糧食需求正在迅速增長。Baj?elj等認為,糧食需求將在21世紀中葉以后超過糧食產量[1-3]。小麥是世界上最重要的3種糧食作物(小麥、水稻和玉米)之一,全球種植面積超過2.2億hm2/年[4]。中國是全球第一大小麥生產國,其產量約占全球小麥產量的18%[5]。隨著糧食需求的增加,中國需要將糧食產量提高36%才能養活全國[6]。因此,及時準確掌握各地區農作物的產量信息,能夠為糧食生產宏觀調控、經濟政策制定和農作物保險提供信息支持,對服務國家糧食安全戰略具有重要意義。一般而言,基于過程的作物生長模擬模型和統計模型是用于產量預測的常用模型。作物生長模擬模型(如WOFOST[7]、DSSAT[8]、APSIM[9]、STICS[10]和MONICA[11])可以模擬農作物的生長發育和產量形成過程,能夠研究作物產量和環境條件之間的相互作用[12-14]。然而,大規模運行作物模型非常耗時[15-16],還需要大量來自實地統計的數據(如品種特征、管理信息、土壤特性等)以校準模型[2,17]。因此,難以將該方法擴展到多種作物和地區[18-20]。與基于過程的作物模型相比,統計模型(如傳統的統計模型和機器學習模型)在大規模產量預測中更加有效,因此應用更廣泛[20]。已有研究使用傳統的回歸模型進行產量預測。如Lobell等使用多元線性回歸(MLR)模型預測1980—2003年美國加利福尼亞州的作物產量,且僅使用2～3個氣候變量的簡單方程解釋超過2/3觀測到的產量變化[21]。但與非線性回歸模型相比,傳統的統計模型(如線性回歸模型)準確性相對較低。實際上作物產量與大多環境因素之間的關系通常是非線性的[22-24]。機器學習是一種先進的統計技術,可以分析變量與產量之間的非線性關系[25-28]。最近,有研究開發出基于統計的作物產量預測模型,如人工神經網絡(ANN)、最小絕對值收斂和選擇算子算法(LASSO)、支持向量機(SVM)和隨機森林 (RF)[29-32]。Khanal等使用線性回歸和6種機器學習算法預測玉米產量并比較它們的性能[33]。Leng等使用傳統的線性回歸模型和RF模型預測1980—2010年美國玉米產量的變化[34]。眾多研究結果表明,機器學習模型的性能優于線性回歸模型。此外,部分機器學習模型能夠識別預測變量對產量的相對重要性。作物產量受多種因素影響,如作物品種特征、土壤特征、病蟲害等[35],但是先前的大多數研究都基于單一因素預測作物產量,如氣候數據。Anderson等認為,衛星遙感可以連續監測不同光譜波段的作物生長,并為作物產量預測提供有用的補充信息[36-38]。其中,可見光和近紅外數據可以用來計算各種植被指數以監測作物生長,故其在預測作物產量方面最具優勢[39-42]。自從Tucker提出第1個歸一化植被指數(NDVI)以來[43],幾種流行的植被指數[NDVI、增強植被指數(EVI)]已被廣泛應用于農業領域[44-45]。EVI是NDVI的改進版,對較高的冠層葉面積指數敏感,受大氣氣溶膠的影響較小,是最常用的作物產量預測指數。但植被指數以綠化率為基礎,對溫度、水汽壓、吸收輻射等氣象因素引起的植被生理變化不敏感。近年來,日光誘導葉綠素熒光(SIF)這一衛星指數進入眾多研究者的視線。SIF是植物葉綠素吸收光能后產生的熒光,反映光合作用效率和響應。在農業、生態學和氣候研究中,SIF被廣泛應用,幫助了解植物對環境的適應性和生長狀態,為農作物生產和生態系統保護提供有價值的信息[46-49]。大量研究結果表明,日光誘導葉綠素熒光(SIF)與生物量直接相關[50-53],可以直接反映作物的呼吸作用,對環境變化作出及時、準確的響應。Somkuti等直接使用SIF來預測作物產量,并取得比植被指數更好的結果[54-56]。同時,由EVI再度改進出一個新的植被指數,即無藍光增強植被指數(EVI2),它不僅繼承了EVI的全部優點,還可以消除亞像素和云等造成的藍光波段的影響。Franch等認為,基于氣候數據和衛星數據的作物模型的產量預測精度普遍優于僅基于氣候數據的模型[57-59]。Guan等認為,各種衛星產品具有共同的重疊和互補信息有利于產量預測[60]。但如何更好地結合衛星數據和其他環境因素來進行作物產量預測還有待進一步研究。冬小麥的生長周期可細分為11個階段,分別是發芽期、出苗期、分蘗期、越冬期、返青期、拔節期、孕穗期、抽穗期、開花期、灌漿期、成熟期。Li等將這11個時期分為4個階段[基于生長期(Bogp)]：T1,發芽期—分蘗期(10—11 月);T2,分蘗期—拔節期(12月至次年3月);T3,拔節期—抽穗期(3—4月);T4,抽穗期—成熟期(5—6 月),并按照這4個階段整理數據集[61]。但這樣的時間劃分方式存在很大的干擾因素。如10月為播種階段,5月下旬至6月上旬為收獲階段,人為因素影響巨大。本研究將多源時空數據(包括時空數據、氣候數據、土壤數據、衛星數據)按照月份進行拆分,去除受人為因素影響最大的10月以及次年的5、6月[基于時空(Boam)],采用LASSO、嶺回歸(RIDGE)、支持向量機回歸(SVR)、隨機森林、XGBoost、輕型梯度提升算法(LightGBM)等6個機器學習算法,對山東省2009—2020年13個縣(市區)的冬小麥產量進行預測,并進行相關研究。本研究旨在解決以下幾個問題：(1)采用本研究的時間劃分方式整理數據集訓練模型,其性能與以往的研究相比是否有改進?這6種機器學習模型誰的性能最強?最佳預測時間段是哪幾個月?(2)多源時空數據對預測模型是否有積極影響?時空變量對預測模型是否有影響? EVI2與SIF對提高預測精度的作用誰更大?

1 材料與方法

1.1 研究區域

山東省是我國冬小麥主產區之一,屬溫帶大陸性季風氣候。該地區主要盛行的輪作方式是冬小麥和夏玉米,本研究僅針對冬小麥進行分析。同時,山東省內各地所種植冬小麥的品種相近,所以忽略基因型的影響(圖1)。

1.2 數據和預處理

本研究共使用了作物產量和種植面積、氣候數據、衛星數據、土壤數據、時空數據等數據,數據集描述見表1。本研究收集的多源數據具有不同的時空分辨率,因此首先將所有變量聚合為縣級空間分辨率和月份時間分辨率,再匯總每個縣(市、區)的月平均變量。因為各變量每個月對冬小麥的生長具有不同的影響,所以將每個變量按月份劃分為6個因子,作為預測模型的輸入。

表1 數據集的詳細描述

1.2.1 作物產量和種植面積 基于山東省各縣(市、區)的農業統計年鑒(http://tjj.shandong.gov.cn/col/col6279/index.html),本研究收集13個縣(市、區)2009—2020年冬小麥種植面積(hm2)和總產量數據(kg)。有一些縣域缺失某些年份的數據,再剔除氣候數據、衛星數據和土壤數據等不全的異常數據,共得到110條樣本數據。通過冬小麥總產量和種植面積計算得到冬小麥產量(kg/hm2)。

1.2.2 氣候數據 氣候數據來源于MERRA-2數據集[62]。MERRA-2是一套由美國國家航空航天局(National Aeronautics and Space Administration,NASA)提供的大氣再分析資料。本研究使用其中的8個氣象因子,包括地表溫度(℃)、降水量(cm)、相對濕度(%)、蒸發量(cm)、平均風速(m/s)、長波輻射(W/m2)、短波輻射(W/m2)和日照時數(峰值,h)。

小麥屬于喜涼性作物,溫度、降水、太陽輻射和光照是影響其產量的關鍵因素。本研究所選的8個氣象因子中,地表溫度、降水量、相對濕度和蒸發量是反映溫度和降水的變量,長波輻射、短波輻射和日照時數是反映太陽輻射和光照的變量,這7個變量已被廣泛應用于作物產量研究[63]。另外,風速是一個很少被考慮的氣候因素,它可以通過影響葉片氣體和熱交換來影響植物生長,或通過改變葉片邊界層的厚度來改變水分狀態[64-65]。

1.2.3 衛星數據 本研究使用SIF和EVI2共2個衛星植被指數。SIF是近年來研究光合活動的一種新指數,已廣泛應用于作物監測和產量預測。SIF數據下載自OCO-2衛星(http://data.globalecology.unh.edu)的SIF再分析數據集(GOSIF)。GOSIF是一個再分析數據集,基于來自OCO-2、MODIS數據和氣候數據的SIF數據,與直接從OCO-2探測中聚合的粗分辨率SIF相比,GOSIF具有更精細的空間分辨率,全球連續覆蓋和更長的記錄。EVI2是EVI的改進版本,它僅由紅光波段和近紅外波段2個波段計算得出。與EVI相比,它不僅繼承了EVI的全部優點,還可以消除亞像素和云等造成的藍光波段的影響。EVI2數據下載自美國國家航空航天局的MOD13A3產品(https://search.earthdata.nasa.gov/search),該產品的數據不僅包含NDVI和EVI等2個現成的植被指數,還包括紅光、近紅外、藍光、中紅外等4個波段的原始數據。EVI2由公式(1)計算得出。

(1)

式中：NIR表示近紅外波段;RED表示紅光波段。

1.2.4 土壤數據 土壤理化性質是作物產量的關鍵影響因素,本研究共使用6個土壤變量,其中土壤容重、土壤深度、有機碳含量、酸堿度pH值和黏土的陽離子交換能力等5個變量[63]來源于世界土壤數據庫(HWSD,http://www.tpdc.ac.cn/zh-hans/data/844010ba-d359-4020-bf76-2b58806f9205/)[66],HWSD里中國地區的數據源為中國科學院南京土壤研究所提供的1995年全國第2次土地調查的土壤數據。使用的最后1個土壤變量為根區濕度(kg/m2),來源于美國哥達德空間飛行中心(GSFC)和美國國家環境預報中心(NCEP)聯合開發的全球高分辨率的陸面模擬系統(GLDAS),空間分辨率為0.25°×0.25°,時間分辨率為逐月。

1.2.5 地理基礎數據 作物生長狀況和生長環境具有空間異質性,相鄰縣在某一年的作物產量通常是相似的。空間自相關可以通過地理坐標(經緯度)來解釋[67]。

1.2.6 數據預處理 EVI2、SIF和HWSD土壤數據均是柵格數據,其中EVI2和SIF均是用ENVI軟件從原始數據中導出表格數據,再根據經緯度求縣級月平均數據。HWSD土壤數據則是先用ArcGIS軟件從原始數據庫中導出所需要的表格數據,再根據經緯度求縣級平均數據。為了消除指標之間的量綱影響,對所有數據作最大值、最小值歸一化處理。

共整理出3個數據集：第一,按照Boam劃分方式,先將10月以及次年5月、6月的數據剔除掉,再把剩下的數據按照月分辨率整理成Boam數據集。第二,將數據按照Bogp劃分方式整理出來的Bogp數據集。第三,先通過收集到的全部數據計算出每個變量的日均值,再刪除經緯度坐標,以此得到Nost(無時空變量)數據集。

1.3 研究方法

采用最小絕對收縮和選擇算子算法(LASSO)、嶺回歸(RIDGE)、支持向量機回歸(SVR)、隨機森林(RF)、極端梯度提升回歸(XGBoost)、輕量級梯度提升算法(LightGBM)等6種機器學習算法構建冬小麥產量預測模型,并比較它們的性能。

1.3.1 最小絕對值收斂和選擇算子算法和嶺回歸 LASSO由Robert Tibshirani于1996年首次提出,是一種正則化回歸分析方法,通過變量選擇和正則化可以提高模型的預測精度和可解釋性。LASSO通過使用特征選擇生成簡約模型,并通過懲罰系數絕對值之和來避免過擬合。RIDGE是1970年由Hoerl和 Kennard提出的一種算法,也是一種正則化回歸方法,但RIDGE使用L2正則化,而LASSO使用L1正則化。LASSO和RIDGE的模型使用Python 3.9中的scikit-learn包進行訓練。

1.3.2 支持向量機回歸 SVR是SVM的一種回歸方法,是一種基于核的算法。與支持向量機類似,SVR找到一個回歸平面來最小化所有輸入到這個超平面的距離。一般而言,SVR需要1個核函數將原始空間的所有輸入映射到高維空間,然后在該特征空間中構造一個線性函數來平衡誤差最小化和過擬合。最常用的核函數有線性核、多項式核和高斯徑向基核。此外,需要調優的超參數是懲罰系數C和核系數gamma,該模型使用Python 3.9中的scikit-learn包進行訓練。

1.3.3 隨機森林 RF最早由Leo Breiman和Adele Cutler在1995年提出。RF模型是一種用于回歸和分類的集成學習方法,它由許多決策樹組成。在訓練階段將訓練集分成多個不同的子訓練集,再用每個子訓練集生成單個決策樹,每個決策樹訓練完成后都會給出一個預測結果。最后對每個決策樹的預測結果取平均值,得到最終的預測結果。此外,RF對噪音不敏感,具有很強的魯棒性。決策樹的數量和最大深度是研究中需要調優的超參數。該模型使用Python 3.9中的scikit-learn包進行訓練。

1.3.4 極端梯度提升算法 XGBoost是2014年華盛頓大學的陳天奇博士提出的一種梯度增強框架下的算法。XGBoost算法的基本原理是在完整數據上構建多個弱學習器,并將所有弱學習器的建模結果進行聚合,以獲得更好的回歸或分類性能。XGBoost結合了正則化模型來防止過擬合,弱學習器可以是回歸樹或線性模型。本研究的 XGBoost 模型是基于決策樹建立的。該模型使用Python 3.9中的XGBoost包進行訓練。

1.3.5 輕型梯度提升算法 LightGBM是微軟于2017年提出的boosting框架,其基本原理與XGBoost一樣,使用基于學習算法的決策樹,只是在框架上做了優化(重點在模型訓練速度上的優化)。LightGBM已經被應用于許多不同類型的數據挖掘任務(如分類、回歸和排序),并顯示出極好的準確性。LightGBM算法包含基于梯度的單邊采樣(GOSS)和互斥特征捆綁(EFB)2種新技術,這2種技術分別便于處理大量數據實例和大量特征。因此,與其他類似算法(如XGBoost)相比,LightGBM在計算速度和內存占用方面的表現明顯更好。每棵樹的葉子數量、迭代速度、樹的最大深度、每張葉子記錄的最小數量、每次迭代隨機選擇的特征的比例、每次迭代使用的數據的比例是LightGBM算法中需要調優的主要參數。該模型使用Python 3.9中的lightgbm包進行訓練。

1.3.6 模型評估 本研究使用留一法交叉驗證來評估模型的實用性,即使用除某縣某年的其他所有縣所有年份的數據來訓練模型,再對目標年份進行預測。該方法是一種被廣泛使用的交叉驗證方法,其具有實現簡單、通用性強、能避免過擬合問題等優點。為了評估模型的性能,本研究選擇確定系數(R2)和均方根誤差(RMSE)作為評估指標。R2是最常用于評價回歸模型優劣程度的指標,其取值范圍為0～1,無單位。R2越大(接近于1),所擬合的回歸方程越優。RMSE衡量觀測值與真實值之間的偏差,常用來作為機器學習模型預測結果衡量的標準。RMSE越小,預測結果越接近真實值。RMSE的計算公式見公式(2)。R2和RMSE由Python 3.9中的scikit-learn包里的函數計算得出。

(2)

1.3.7 試驗設計 本研究共設計5個試驗來探究問題。試驗1：分別用Bogp數據集和Boam數據集訓練6個機器學習模型,比較哪個數據集訓練出的冬小麥預測模型性能更強,再找出最強模型。試驗2：Boam數據集中,氣候數據、衛星數據及土壤數據中的根區濕度按照月份劃分,其他為固定變量。本試驗抽取10個不同時間段(11月至次年1月、11月至次年2月、11月至次年3月、11月至次年4月、12月至次年2月、12月至次年3月、12月至次年4月、1—3月、1—4月、2—4月)的數據對模型進行訓練,用得到的R2和RMSE選出預測效果最佳的時間段。為了避免機器學習模型的過擬合現象,每個時間段最低包含3個月的數據。試驗3：先用包含年份、經緯度和氣候變量的數據對模型進行訓練并預測,依次向數據集增加土壤變量和衛星變量,并分別對模型進行訓練測試,最后對比3次訓練得到的評估指數,以此探究多源數據對產量預測結果的影響。試驗4：用Nost數據集訓練模型,將預測結果與用Boam數據集得到的結果對比,以此探究時空變量對預測結果的影響。試驗5：對衛星數據進行拆分,先用EVI2和其他變量作為模型的輸入,再用SIF和其他變量作為模型的輸入,最后用EVI2+SIF和其他變量作為模型的輸入,以此測試2個衛星因子的性能。

2 結果與分析

2.1 數據集的時間性能研究

由圖2可知,除嶺回歸模型外,其余5種模型均顯示,Boam劃分方式略優于Bogp劃分方式。對比Boam數據集訓練出來的6個預測模型可知,XGBoost模型的預測性能最強(R2=0.843 6,RMSE=0.085 8),其次是LightGBM模型(R2=0.680 7,RMSE=0.102 3)和RF模型(R2=0.518 7,RMSE=0.116 4),而嶺回歸模型(R2=0.096 5,RMSE=0.162 4)的性能最差。可能是因為大多數變量與結果之間的關系是非線性的。

試驗2用Boam數據集來完成,并僅用基于XGBoost算法的冬小麥產量預測模型來分析。10個時間段數據訓練出的模型的R2和RMSE見圖3。可見10個模型的性能差距不大,其中使用11月至次年1月數據訓練出來的模型性能最好,R2、RMSE分別為0.850 0、0.085 1。故11月至次年1月是最佳預測時間段。

2.2 對預測變量的研究

試驗3的R2和RMSE見圖4。可見添加土壤變量后,R2由0.664 3增加到0.794 9,RMSE由0.103 8減小到0.091 0。最后加上衛星變量,R2又增加到0.843 6,RMSE減小到0.085 8。可見隨著數據的多源化,預測精度在逐步增加。圖5是試驗4的R2和RMSE。可見添加時空變量后,R2由0.713 5增加到0.843 6,RMSE由0.106 4減小到0.085 8,說明添加時空變量對預測結果有一定的積極作用。

試驗5中3次訓練得到的R2和RMSE見圖6。與SIF(R2=0.732 2,RMSE=0.097 4)相比,EVI2(R2=0.797 5,RMSE=0.089 0)對預測結果的增幅更大。但用同時包含兩者的數據訓練出來的模型(R2=0.843 6,RMSE=0.085 8)性能更加強大。SIF和EVI2對預測結果的增益方式不同,故累加兩者之后預測效果更強。

3 結論與討論

3.1 數據集的時間性能研究

與Bogp劃分方式相比,按照Boam劃分方式整理數據集訓練的模型性能更強。不可否認的是,以前研究的時間劃分有一定的道理,但是沒有排除人為因素的影響,就像10月不僅是發芽期還是播種期,5月和6月不僅是成熟期還是收獲期,這必然對預測結果產生一定的影響。排除人為干擾,僅用11月至次年4月的數據進行預測,同時把變量按照月份進行劃分,時間劃分更細致。

測試的6種機器學習模型中,基于XGBoost算法的冬小麥產量預測模型性能穩居第一,但是此模型也有缺點,如占用內存較大、運行時間較長。與之相比,LightGBM算法作為XGBoost算法的改進版,占用內存較小,運行速度較快。而LightGBM算法是犧牲了部分精度才獲取到最優速度。因此,是使用基于XGBoost算法的模型還是使用基于LightGBM算法的模型,又或是其他機器學習模型,要視具體情況而定。

試驗2旨在找到最佳預測時間段,因此抽取10個不同時間段的數據進行預測,發現使用11月至次年1月數據訓練模型時,得到R2最大值為0.850 0,RMSE最小值為0.085 1。可見,11月至次年1月是最佳預測時間段。11月至次年1月的預測結果優于11月至次年4月的預測結果。據此推斷,數據集中的74個變量存在無關變量或消極變量。進而計算每個變量的貢獻度,貢獻度的取值范圍為0～1,全部變量的貢獻度總和為1。

由表2可知,上年11月日照時數、上年12月日照時數、1月日照時數、3月日照時數、4月日照時數、3月EVI2等6個變量貢獻度為0,即無關變量。其中,3月SIF、2月EVI2、上年12月SIF等3個變量為貢獻度最高的3個變量,其貢獻度均在0.1之上。值得注意的是,貢獻度最高的3個變量全都是衛星因子,由此可預見衛星數據在農學領域的潛力。此外,其他變量貢獻度均在0～0.1之間。

表2 變量貢獻度

3.2 對預測變量的研究

泱泱華夏地大物博,不同時間不同地區的氣候、人文、土壤等均不相同。試驗3探究了多源數據對冬小麥產量預測模型的影響,數據集由多源數據整理而成,結果表明,每增加1組不同來源的數據,預測結果精度便增加1分。可見,數據越多源,預測結果精度越高。本研究共使用時空數據、氣候數據、土壤數據、衛星數據等4個來源的數據進行預測,像社會因子(灌溉面積、化肥農藥等)、高程數據等均未使用。故本研究還可以繼續收集更多源的數據整理進數據集,以使模型性能更強,預測結果更準確。

試驗4探究時空變量對預測結果的影響,發現與用不包含時空變量的Nost數據集訓練出的預測模型相比,用加上時空變量之后的Boam數據集訓練出的預測模型精度更高。精度提升較低,這是因為僅選擇山東省作為研究區域,如果把華北平原甚至全中國、全世界作為研究區域,那時空變量一定會給預測結果帶來更大的提升。

試驗5探究EVI2和SIF對預測結果的影響。在農業領域,已有學者使用衛星數據來源預測作物產量,使用最多的是NDVI和EVI這2個植被指數。本研究使用的EVI2指數是EVI的改良版,至今未曾被用于農業領域。它僅由紅光波段和近紅外波段計算得出,不受藍光波段的影響,故性能強于EVI,本研究探究其和近期興起的植被指數SIF的性能誰更強。結果表明,EVI2的性能強于SIF,但兩者的疊加效果又強于兩者單獨的作用。可見,2個植被指數的作用域不同,預測時可疊加兩者使模型性能變得更強。

3.3 結論

本研究開發基于機器學習算法的冬小麥產量模型,通過整合多源時空數據(時空數據、氣候數據、土壤數據、衛星數據)來預測冬小麥產量,并得出以下結論。

(1)與以往的研究相比,Boam時間劃分方式更合理,因而預測精度更高。測試的6個冬小麥產量預測模型中,XGBoost模型在預測冬小麥產量方面具有最優性能。6個機器學習模型的性能排序依次是XGBoost模型(R2=0.843 6,RMSE=0.085 8)、LightGBM模型(R2=0.680 7,RMSE=0.102 3)、RF模型(R2=0.518 8,RMSE=0.116 4)、SVR模型(R2=0.495 6,RMSE=0.118 3)、LASSO模型(R2=0.273 4,RMSE=0.135 0)、RIDGE(R2=0.096 5,RMSE=0.162 4)。結果表明,11月至次年1月是進行冬小麥產量預測的最佳時間段。由6個月縮減到3個月,不僅可以大大縮減整理數據集的工作量,還可以增強模型的性能。數據集中存在無關變量。3月SIF、2月EVI2、上年12月SIF等3個變量為貢獻度最高的3個變量,均為衛星因子,由此可預見衛星數據在農學領域的潛力。

(2)更多源的數據能給預測結果帶來更大的增益。首先,在數據集中依次增加土壤數據和衛星數據后,模型的預測性能逐漸增強。其次,時空變量對預測結果具有積極的作用。用Boam數據集得到的預測結果比用沒加時空變量(年份和經緯度)的數據得到的預測結果精度更高。最后,與SIF相比,EVI2的性能更勝一籌,而同時使用SIF和EVI2共2個因子的模型性能強于僅適用SIF或EVI2的模型。因此,2個衛星植被指數的作用域不同,在進行冬小麥產量預測時完全可以疊加兩者使用。

本研究證明,以XGBoost為主的幾個機器學習算法在冬小麥產量預測方面的能力很強,且建模方法可以通過使用公開的數據擴充本研究的數據集,進而廣泛應用于全世界其他地區。