林 瀅,邵懷勇
(成都理工大學地球科學學院,四川成都 610059)
我國是人口大國,糧食安全至關重要。小麥是我國主要糧作物之一。據統計,2018年我國小麥播種面積2 427萬hm2,總產量13 143萬t,占全球小麥產量的17.62%[1]。及時準確掌握小麥生長信息,并進行估產分析對農田管理、農業政策制定等具有重要意義[2-3]。目前作物產量預測主要采用經驗統計模型和面向過程的作物生長模型。傳統方法中,許多國家的作物產量統計數據通常是將地面觀測與產量報告相結合計算得出。如Reynolds[4]在2010年將實時衛星圖像引入到地理信息系統(GIS)和糧食及農業組織(FAO)的作物特定水平衡(CSWB)模型中,開發了一種可操作的作物單產模型,并取得較好效果。但是基于地面實地考察的數據收集成本高昂、耗時且效率低下,同時無法確定過程中的可靠性,因此可能會導致作物產量評估效果不佳[5-6]。且經驗模型往往因地理位置、作物品種和生長季節而異,模型的空間泛化能力低。遙感技術的出現為農業研究提供了新的方法[7]。充分利用遙感數據,可實現農作物長勢監測與產量估算,研究作物遙感估產的基本機理與方法[8-9]。如吳炳方[10]利用每旬的最大NDVI圖像對全國范圍的作物進行遙感長勢監測,估算作物種植的面積。作物估產模型包括農業技術轉移決策支持系統(DSSAT)、農業生產系統模擬器(APSIM)、捕捉大面積作物-天氣關系的模式(MCWLA)和世界糧食研究模型(WOFOST)[11]。雖然模型可以更高的精度模擬作物產量,但運行模型需要大量的參數輸入(例如氣候變量、肥料、灌溉、土壤和水文特征),費時費力,當實驗區較大時會因無法獲取某些數據而受到限制[12]。隨著人工智能快速興起,機器學習作為新技術,在數據挖掘與分析方面展現了強大的功能,從而為農業應用提供了新的技術與途徑(包括作物類型分類和產量預測),并推動了農業的發展[13]。因其在處理多維數據集方面的強大能力,機器學習技術將為改進產量預測模型提供強有力的支持。近年來,已有機器學習方法運用于作物產量估測,如人工神經網絡[14]、高斯過程[15]等。然而,許多研究基于整個生長季節選擇變量,在收獲日期之前實際上難以估計最終產量。此外,很少有研究確定冬小麥產量預測的最佳訓練時間段,而找到能夠很好地反映冬小麥估產的最佳時間窗,將有助于提高作物估產模型的應用效果。
隨機森林(random forest,RF)是一種需要模擬和迭代的基于分類樹的算法[16]。它可以在不增加運算量的情況下保持良好的準確性[17],并且可以評價自變量的重要性,避免回歸分析中多重共線性現象[18]。目前,隨機森林算法已應用于農業研究,例如Yvette Everingham使用隨機森林算法根據不同大小預測范圍對甘蔗產量進行預估并取得良好效果[19];王娣運用隨機森林算法建立水稻估產模型,并進行模型精度評價[20]。
已有學者對比不同機器學習算法在作物估產中的效果[21],但就不同時間段訓練樣本對預測模型精度影響的討論不多。本研究嘗試以河南省113個縣的冬小麥為例,運用隨機森林算法,探討河南省訓練樣本選擇的最佳時間段,并分析不同影響因素對產量預測的相對重要性,以期提高該算法在作物估產的應用效果。
河南省地處北緯31°23′-36°22′、東經110°21′-116°39′之間,位于黃河中下游,地勢西高東低,大部分區域屬暖溫帶氣候。河南省小麥種植區屬于黃淮海冬麥區,耕作制度為一年兩熟。土壤肥沃,生產環境佳,適宜小麥生長,是我國冬小麥的核心生產區之一。2017年,河南省冬小麥播種總面積達547.5萬hm2,總產量355億kg,占全國總產的四分之一[1],因此其小麥高產穩產對全國小麥生產與供求平衡具有重要影響[22]。
試驗獲取的數據包括2001-2015年的遙感、氣候、土壤和產量數據。所有數據空間分辨率統一為1 km×1 km,時間分辨率統一為一個月,并且所有變量將基于小麥生長像素進行掩膜,并按縣求平均值。數據處理主要在ArcGIS進行。
1.2.1 遙感數據
歸一化植被指數(normalized vegetation index,NDVI)和增強型植被指數(enhanced vegetation index,EVI)與作物產量有較好的相關性[23-25],因而常被用于作物估產研究,將NDVI與EVI結合可為作物估產提供更多的信息[26]。河南省2014和2015年的兩種植被指數來自于MOD13Q1,周期為16 d,空間分辨率為250 m×250 m。
1.2.2 氣候數據
氣候對農作物的產量、生長發育、種植制度均有重要影響[27-28]。參照前人的研究[29],本研究選取每月最高溫度(Tmax)、最低溫度(Tmin)、干旱指數(Psdi)和降水量(Pr)作為氣溫變化要素參與作物產量預測。采用Terra Climate數據集提取出研究時間段內所需氣候指標[30],在GEE平臺處理數據并計算每個縣的氣候變量。
1.2.3 土壤數據
土壤理化性質是作物產量的關鍵影響因素[31-32]。本研究選取土壤深度、土壤質地、有機碳含量、pH、陽離子交換容量、地表土層容重和地下土壤層容重作為土壤影響因子。土壤理化性質數據來自世界土壤數據庫(HWSD)[33]。
1.2.4 冬小麥產量數據
冬小麥產量數據收集自《中國農業年鑒》和縣級統計數據,個別區域數據缺失。
冬小麥不同生育階段的形態、生理特征不同,因而估產選擇不同生長階段的數據進行產量預測的精度不同。河南省冬小麥的一般9月份播種,來年6月收獲。本研究從冬小麥生長季中抽取8個時長不同的時間段(9-5月、9-6月、10-3月、10-4月、10-5月、11-3月、11-4月、12-3月)數據,其中以2001-2013年的數據作為測試數據訓練模型,進行十折交叉驗證。用訓練好的模型分別預測2014和2015年河南的冬小麥產量,與實際產量做精度對比選擇出最佳的樣本選擇時間段。
1.3.1 隨機森林算法預測作物產量
隨機森林算法由一系列不同的回歸樹(CART)組成基于多個分類樹的算法,它對數據集的適應能力較強,能有效地運行大數據集。由于隨機性的引入,隨機森林法不容易陷入過擬合并且具有很好的抗噪聲能力,提高了學習的穩定性。目前,該算法已在生態學[35]、水利水電[36]、地災評估[37]等領域有所應用。主要公式如下:
(1)
式中,F(x)是一個組合模型,hi是單一決策樹,Y表示輸出變量,I表示特征函數。試驗使用的樹數量為100。數據不被提取的概率為1-1/N,收斂為1/e≈0.368,即有約37%的訓練數據不會被運用到單個模型構建中。邊緣函數如下:

(2)
該式表達模型的可靠性,即函數值越大,分類越可靠。分類器的通用表達如下:
PE*=Pxy[mg(X,Y)<0]
(3)
其中,(X,Y)為概率空間。隨著決策樹數目的增加,PE*序列將變為:
Pxy{Pθ[h(X,θ)]=Y-maxPθ[h(X,θ)=J]}<0
(4)
當樹的數量增加,泛化誤差總是收斂的。具體模型構建原理參照Breiman等[16]方法。
1.3.2 精度評估指標
用決定系數(R2)、均方根誤差(root mean square error,RMSE)[38]、平均絕對誤差(mean absolute error,MAE)[39]和誤差系數評估模型預測的精度。
(5)
(6)
(7)
(8)

1.3.3 重要性評價
為了探索模型中不同預測變量的重要性,可計算基于均方誤差(MSE)的預測精度下降的平均值。準確性平均值的下降表明,如果排除該特定變量,則隨機森林模型預測精度也將下降。預測變量準確度的平均值下降幅度越大,認為該變量就越重要[40]。
本研究以產量為因變量,NDVI、EVI、Tmax、Tmin、Psdi、Pr和土壤水分七個每月變化的因素為自變量,運用R語言做隨機森林回歸,可以得到自變量的相對重要程度。
由表1可知,整體上利用不同時段數據訓練出來的模型對小麥產量的預測精度沒有太大差異。在2014年,11-3月和12-3月時段的模型預測精度較好,R2分別是0.81和0.80,MAE和RMSE值均最小。2015年,只有12-3月時段的模型預測精度最好,R2為0.81,MAE和RMSE也都最小。兩年相比,2014年的預測精度整體上優于2015年。2015年有6個時間段RMSE大于1 000 kg·hm-2,預測精度不夠好。

表1 2014、2015年隨機森林算法在不同時間段小麥產量預測精度對比Table 1 Comparison of accuracy in different time period by random forest algorithm
綜合分析表明,12-3月時段的R2在2014和2015年均大于0.8,因而將該時段作為河南省冬小麥最佳訓練樣本選擇時間段。該時段的預測精度較高可能是因為該時間段冬小麥處于返青期,植株生長及氣候特征相關性較高。從擬合曲線上看,在低產區,預測值低于實際值;在高產區,實際值略高于預測值(圖1)。

圖1 2014、2015年小麥預測產量與實際產量散點圖Fig.1 Scatter plot of predicted and actual yield of wheat in 2014 and 2015
將隨機森林算法通過12-3月的數據預測的2014和2015年冬小麥產量及誤差系數進行空間可視化,結果如圖2所示。整體上,河南省冬小麥2014年和2015年的實際產量空間分布狀態相當,均呈現西低東高的狀態。對比發現,兩年預測產量分布與實際產量特點大體相似,但在高產區存在整體預測結果較低的情況。結合誤差系數可以發現,預測結果整體效果較好,大部分區縣的誤差系數介于-0.1~0.1之間,存在個別預測值過高或過低的情況。2014年研究區西部存在局部估產過高的情況,2015年東部有個別過低估計產量的區縣。

圖2 2014、2015年小麥預測產量及誤差系數空間分布圖Fig.2 Distribution of wheat yield and error coefficient
用R語言對冬小麥整個生長期數據進行隨機森林建模回歸分析,結果(圖3)表明,月降水對小麥產量的影響遠大于其他因素,重要性達到了29.79,即當降水發生變化時,對模型精度的影響最大,與水分是影響作物產量的重要環境因子的事實相符。其次是月最低溫度和增強植被指數,重要性分別是23.76和22.64。NDVI、干旱指數、土壤水分的重要性相當,分別是21.05、21.05和21.04;最后是月最高溫度,重要性只有20.75。

圖3 影響因子重要性統計圖Fig.3 Importance of impact factors
作物產量快速預測可為糧食銷售決策制定提供參考依據,同時可以指導整地、除草、施肥等田間管理。本研究利用隨機森林算法實現了冬小麥產量的預測。利用作物生長模型如農業生產系統模擬器(APSIM)進行估產時,通常需要大量的輸入變量,不僅有許多假設,而且輸入變量的數據大多需要田間測量,在大范圍的產量預測方面具有一定的局限性,且大范圍的田間調查需要消耗大量的人力物力。隨機森林等機器學習算法的一個顯著優點是在較少的假設下,可以組合能免費獲取的遙感等數據,通過信息挖掘較好地實現大范圍的作物產量預測,過程相對簡單,且具有通用性的潛力,但相比于作物模型,該方法無法表達各因素對產量影響的具體機理。同時,本研究發現,利用不同生長時段的樣本建模,模型的預測精度不同,表明變量的時段是模型非常重要的影響因素之一。劉峻明等[41]探討了基于隨機森林算法的河南省冬小麥平均拔節期到平均抽穗期(3-5月)氣象要素對產量預測的作用,主要引入氣象特征。本研究主要探討基于冬小麥整個生長期(9月-來年6月)的最佳預測時間窗的選擇,同時考慮到氣候因子對產量預測的影響,引入植被指數等與作物生長相關的要素,使訓練模型的因子選擇更豐富。由于氣候、地形、緯度等差異,作物產量預測的最佳時間窗口可能存在著空間差異性。此外,農戶的管理措施(如作物品種、播種量等)對作物產量也具有較大影響,這些信息應該包括在未來的建模方法中,即將管理投入與模型結合是未來研究的一個有前景的途徑。
選擇2001-2013年河南省冬小麥八個時間段數據作為訓練樣本,應用隨機森林算法建立基于訓練樣本的估產模型,進而預測2014和2015年產量,其中12-3月的估算精度在兩年都最佳,預測產量與實際產量在空間分布上也基本一致,該時段可作為隨機森林算法估算河南省冬小麥產量選擇訓練樣本的最佳時間;在影響因子中,月降水對冬小麥產量的重要性最大,其次是月最低溫度和EVI,NDVI、干旱指數、土壤水分和月最高溫度重要性相當。