林 瀅，邵懷勇

(成都理工大學(xué)地球科學(xué)學(xué)院，四川成都 610059)

我國是人口大國，糧食安全至關(guān)重要。小麥?zhǔn)俏覈饕Z作物之一。據(jù)統(tǒng)計，2018年我國小麥播種面積2 427萬hm2，總產(chǎn)量13 143萬t，占全球小麥產(chǎn)量的17.62%[1]。及時準(zhǔn)確掌握小麥生長信息，并進(jìn)行估產(chǎn)分析對農(nóng)田管理、農(nóng)業(yè)政策制定等具有重要意義[2-3]。目前作物產(chǎn)量預(yù)測主要采用經(jīng)驗統(tǒng)計模型和面向過程的作物生長模型。傳統(tǒng)方法中，許多國家的作物產(chǎn)量統(tǒng)計數(shù)據(jù)通常是將地面觀測與產(chǎn)量報告相結(jié)合計算得出。如Reynolds[4]在2010年將實時衛(wèi)星圖像引入到地理信息系統(tǒng)(GIS)和糧食及農(nóng)業(yè)組織(FAO)的作物特定水平衡(CSWB)模型中，開發(fā)了一種可操作的作物單產(chǎn)模型，并取得較好效果。但是基于地面實地考察的數(shù)據(jù)收集成本高昂、耗時且效率低下，同時無法確定過程中的可靠性，因此可能會導(dǎo)致作物產(chǎn)量評估效果不佳[5-6]。且經(jīng)驗?zāi)Ｐ屯虻乩砦恢?、作物品種和生長季節(jié)而異，模型的空間泛化能力低。遙感技術(shù)的出現(xiàn)為農(nóng)業(yè)研究提供了新的方法[7]。充分利用遙感數(shù)據(jù)，可實現(xiàn)農(nóng)作物長勢監(jiān)測與產(chǎn)量估算，研究作物遙感估產(chǎn)的基本機(jī)理與方法[8-9]。如吳炳方[10]利用每旬的最大NDVI圖像對全國范圍的作物進(jìn)行遙感長勢監(jiān)測，估算作物種植的面積。作物估產(chǎn)模型包括農(nóng)業(yè)技術(shù)轉(zhuǎn)移決策支持系統(tǒng)(DSSAT)、農(nóng)業(yè)生產(chǎn)系統(tǒng)模擬器(APSIM)、捕捉大面積作物-天氣關(guān)系的模式(MCWLA)和世界糧食研究模型(WOFOST)[11]。雖然模型可以更高的精度模擬作物產(chǎn)量，但運行模型需要大量的參數(shù)輸入(例如氣候變量、肥料、灌溉、土壤和水文特征)，費時費力，當(dāng)實驗區(qū)較大時會因無法獲取某些數(shù)據(jù)而受到限制[12]。隨著人工智能快速興起，機(jī)器學(xué)習(xí)作為新技術(shù)，在數(shù)據(jù)挖掘與分析方面展現(xiàn)了強(qiáng)大的功能，從而為農(nóng)業(yè)應(yīng)用提供了新的技術(shù)與途徑(包括作物類型分類和產(chǎn)量預(yù)測)，并推動了農(nóng)業(yè)的發(fā)展[13]。因其在處理多維數(shù)據(jù)集方面的強(qiáng)大能力，機(jī)器學(xué)習(xí)技術(shù)將為改進(jìn)產(chǎn)量預(yù)測模型提供強(qiáng)有力的支持。近年來，已有機(jī)器學(xué)習(xí)方法運用于作物產(chǎn)量估測，如人工神經(jīng)網(wǎng)絡(luò)[14]、高斯過程[15]等。然而，許多研究基于整個生長季節(jié)選擇變量，在收獲日期之前實際上難以估計最終產(chǎn)量。此外，很少有研究確定冬小麥產(chǎn)量預(yù)測的最佳訓(xùn)練時間段，而找到能夠很好地反映冬小麥估產(chǎn)的最佳時間窗，將有助于提高作物估產(chǎn)模型的應(yīng)用效果。

隨機(jī)森林(random forest，RF)是一種需要模擬和迭代的基于分類樹的算法[16]。它可以在不增加運算量的情況下保持良好的準(zhǔn)確性[17]，并且可以評價自變量的重要性，避免回歸分析中多重共線性現(xiàn)象[18]。目前，隨機(jī)森林算法已應(yīng)用于農(nóng)業(yè)研究，例如Yvette Everingham使用隨機(jī)森林算法根據(jù)不同大小預(yù)測范圍對甘蔗產(chǎn)量進(jìn)行預(yù)估并取得良好效果[19]；王娣運用隨機(jī)森林算法建立水稻估產(chǎn)模型，并進(jìn)行模型精度評價[20]。

已有學(xué)者對比不同機(jī)器學(xué)習(xí)算法在作物估產(chǎn)中的效果[21]，但就不同時間段訓(xùn)練樣本對預(yù)測模型精度影響的討論不多。本研究嘗試以河南省113個縣的冬小麥為例，運用隨機(jī)森林算法，探討河南省訓(xùn)練樣本選擇的最佳時間段，并分析不同影響因素對產(chǎn)量預(yù)測的相對重要性，以期提高該算法在作物估產(chǎn)的應(yīng)用效果。

1 數(shù)據(jù)與方法

1.1 研究區(qū)概況

河南省地處北緯31°23′-36°22′、東經(jīng)110°21′-116°39′之間，位于黃河中下游，地勢西高東低，大部分區(qū)域?qū)倥瘻貛夂?。河南省小麥種植區(qū)屬于黃淮海冬麥區(qū)，耕作制度為一年兩熟。土壤肥沃，生產(chǎn)環(huán)境佳，適宜小麥生長，是我國冬小麥的核心生產(chǎn)區(qū)之一。2017年，河南省冬小麥播種總面積達(dá)547.5萬hm2，總產(chǎn)量355億kg，占全國總產(chǎn)的四分之一[1]，因此其小麥高產(chǎn)穩(wěn)產(chǎn)對全國小麥生產(chǎn)與供求平衡具有重要影響[22]。

1.2 數(shù)據(jù)來源

試驗獲取的數(shù)據(jù)包括2001-2015年的遙感、氣候、土壤和產(chǎn)量數(shù)據(jù)。所有數(shù)據(jù)空間分辨率統(tǒng)一為1 km×1 km，時間分辨率統(tǒng)一為一個月，并且所有變量將基于小麥生長像素進(jìn)行掩膜，并按縣求平均值。數(shù)據(jù)處理主要在ArcGIS進(jìn)行。

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

歸一化植被指數(shù)(normalized vegetation index，NDVI)和增強(qiáng)型植被指數(shù)(enhanced vegetation index，EVI)與作物產(chǎn)量有較好的相關(guān)性[23-25]，因而常被用于作物估產(chǎn)研究，將NDVI與EVI結(jié)合可為作物估產(chǎn)提供更多的信息[26]。河南省2014和2015年的兩種植被指數(shù)來自于MOD13Q1，周期為16 d，空間分辨率為250 m×250 m。

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

氣候?qū)r(nóng)作物的產(chǎn)量、生長發(fā)育、種植制度均有重要影響[27-28]。參照前人的研究[29]，本研究選取每月最高溫度(Tmax)、最低溫度(Tmin)、干旱指數(shù)(Psdi)和降水量(Pr)作為氣溫變化要素參與作物產(chǎn)量預(yù)測。采用Terra Climate數(shù)據(jù)集提取出研究時間段內(nèi)所需氣候指標(biāo)[30]，在GEE平臺處理數(shù)據(jù)并計算每個縣的氣候變量。

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

土壤理化性質(zhì)是作物產(chǎn)量的關(guān)鍵影響因素[31-32]。本研究選取土壤深度、土壤質(zhì)地、有機(jī)碳含量、pH、陽離子交換容量、地表土層容重和地下土壤層容重作為土壤影響因子。土壤理化性質(zhì)數(shù)據(jù)來自世界土壤數(shù)據(jù)庫(HWSD)[33]。

1.2.4 冬小麥產(chǎn)量數(shù)據(jù)

冬小麥產(chǎn)量數(shù)據(jù)收集自《中國農(nóng)業(yè)年鑒》和縣級統(tǒng)計數(shù)據(jù)，個別區(qū)域數(shù)據(jù)缺失。

1.3 研究方法

冬小麥不同生育階段的形態(tài)、生理特征不同，因而估產(chǎn)選擇不同生長階段的數(shù)據(jù)進(jìn)行產(chǎn)量預(yù)測的精度不同。河南省冬小麥的一般9月份播種，來年6月收獲。本研究從冬小麥生長季中抽取8個時長不同的時間段(9-5月、9-6月、10-3月、10-4月、10-5月、11-3月、11-4月、12-3月)數(shù)據(jù)，其中以2001-2013年的數(shù)據(jù)作為測試數(shù)據(jù)訓(xùn)練模型，進(jìn)行十折交叉驗證。用訓(xùn)練好的模型分別預(yù)測2014和2015年河南的冬小麥產(chǎn)量，與實際產(chǎn)量做精度對比選擇出最佳的樣本選擇時間段。

1.3.1 隨機(jī)森林算法預(yù)測作物產(chǎn)量

隨機(jī)森林算法由一系列不同的回歸樹(CART)組成基于多個分類樹的算法，它對數(shù)據(jù)集的適應(yīng)能力較強(qiáng)，能有效地運行大數(shù)據(jù)集。由于隨機(jī)性的引入，隨機(jī)森林法不容易陷入過擬合并且具有很好的抗噪聲能力，提高了學(xué)習(xí)的穩(wěn)定性。目前，該算法已在生態(tài)學(xué)[35]、水利水電[36]、地災(zāi)評估[37]等領(lǐng)域有所應(yīng)用。主要公式如下：

(1)

式中，F(xiàn)(x)是一個組合模型，hi是單一決策樹，Y表示輸出變量，I表示特征函數(shù)。試驗使用的樹數(shù)量為100。數(shù)據(jù)不被提取的概率為1-1/N，收斂為1/e≈0.368，即有約37%的訓(xùn)練數(shù)據(jù)不會被運用到單個模型構(gòu)建中。邊緣函數(shù)如下：

(2)

該式表達(dá)模型的可靠性，即函數(shù)值越大，分類越可靠。分類器的通用表達(dá)如下：

PE*=Pxy[mg(X,Y)<0]

(3)

其中，(X,Y)為概率空間。隨著決策樹數(shù)目的增加，PE*序列將變?yōu)椋?/p>

Pxy{Pθ[h(X,θ)]=Y-maxPθ[h(X,θ)=J]}<0

(4)

當(dāng)樹的數(shù)量增加，泛化誤差總是收斂的。具體模型構(gòu)建原理參照Breiman等[16]方法。

1.3.2 精度評估指標(biāo)

用決定系數(shù)(R2)、均方根誤差(root mean square error，RMSE)[38]、平均絕對誤差(mean absolute error，MAE)[39]和誤差系數(shù)評估模型預(yù)測的精度。

(5)

(6)

(7)

(8)

1.3.3 重要性評價

為了探索模型中不同預(yù)測變量的重要性，可計算基于均方誤差(MSE)的預(yù)測精度下降的平均值。準(zhǔn)確性平均值的下降表明，如果排除該特定變量，則隨機(jī)森林模型預(yù)測精度也將下降。預(yù)測變量準(zhǔn)確度的平均值下降幅度越大，認(rèn)為該變量就越重要[40]。

本研究以產(chǎn)量為因變量，NDVI、EVI、Tmax、Tmin、Psdi、Pr和土壤水分七個每月變化的因素為自變量，運用R語言做隨機(jī)森林回歸，可以得到自變量的相對重要程度。

2 結(jié)果與分析

2.1 不同時間段預(yù)測精度比較

由表1可知，整體上利用不同時段數(shù)據(jù)訓(xùn)練出來的模型對小麥產(chǎn)量的預(yù)測精度沒有太大差異。在2014年，11-3月和12-3月時段的模型預(yù)測精度較好，R2分別是0.81和0.80，MAE和RMSE值均最小。2015年，只有12-3月時段的模型預(yù)測精度最好，R2為0.81，MAE和RMSE也都最小。兩年相比，2014年的預(yù)測精度整體上優(yōu)于2015年。2015年有6個時間段RMSE大于1 000 kg·hm-2，預(yù)測精度不夠好。

表1 2014、2015年隨機(jī)森林算法在不同時間段小麥產(chǎn)量預(yù)測精度對比Table 1 Comparison of accuracy in different time period by random forest algorithm

綜合分析表明，12-3月時段的R2在2014和2015年均大于0.8，因而將該時段作為河南省冬小麥最佳訓(xùn)練樣本選擇時間段。該時段的預(yù)測精度較高可能是因為該時間段冬小麥處于返青期，植株生長及氣候特征相關(guān)性較高。從擬合曲線上看，在低產(chǎn)區(qū)，預(yù)測值低于實際值；在高產(chǎn)區(qū)，實際值略高于預(yù)測值(圖1)。

圖1 2014、2015年小麥預(yù)測產(chǎn)量與實際產(chǎn)量散點圖Fig.1 Scatter plot of predicted and actual yield of wheat in 2014 and 2015

2.2 產(chǎn)量預(yù)測空間分布特征

將隨機(jī)森林算法通過12-3月的數(shù)據(jù)預(yù)測的2014和2015年冬小麥產(chǎn)量及誤差系數(shù)進(jìn)行空間可視化，結(jié)果如圖2所示。整體上，河南省冬小麥2014年和2015年的實際產(chǎn)量空間分布狀態(tài)相當(dāng)，均呈現(xiàn)西低東高的狀態(tài)。對比發(fā)現(xiàn)，兩年預(yù)測產(chǎn)量分布與實際產(chǎn)量特點大體相似，但在高產(chǎn)區(qū)存在整體預(yù)測結(jié)果較低的情況。結(jié)合誤差系數(shù)可以發(fā)現(xiàn)，預(yù)測結(jié)果整體效果較好，大部分區(qū)縣的誤差系數(shù)介于-0.1～0.1之間，存在個別預(yù)測值過高或過低的情況。2014年研究區(qū)西部存在局部估產(chǎn)過高的情況，2015年東部有個別過低估計產(chǎn)量的區(qū)縣。

圖2 2014、2015年小麥預(yù)測產(chǎn)量及誤差系數(shù)空間分布圖Fig.2 Distribution of wheat yield and error coefficient

2.3 影響因子重要性評價

用R語言對冬小麥整個生長期數(shù)據(jù)進(jìn)行隨機(jī)森林建?；貧w分析，結(jié)果(圖3)表明，月降水對小麥產(chǎn)量的影響遠(yuǎn)大于其他因素，重要性達(dá)到了29.79，即當(dāng)降水發(fā)生變化時，對模型精度的影響最大，與水分是影響作物產(chǎn)量的重要環(huán)境因子的事實相符。其次是月最低溫度和增強(qiáng)植被指數(shù)，重要性分別是23.76和22.64。NDVI、干旱指數(shù)、土壤水分的重要性相當(dāng)，分別是21.05、21.05和21.04；最后是月最高溫度，重要性只有20.75。

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

3 討 論

作物產(chǎn)量快速預(yù)測可為糧食銷售決策制定提供參考依據(jù)，同時可以指導(dǎo)整地、除草、施肥等田間管理。本研究利用隨機(jī)森林算法實現(xiàn)了冬小麥產(chǎn)量的預(yù)測。利用作物生長模型如農(nóng)業(yè)生產(chǎn)系統(tǒng)模擬器(APSIM)進(jìn)行估產(chǎn)時，通常需要大量的輸入變量，不僅有許多假設(shè)，而且輸入變量的數(shù)據(jù)大多需要田間測量，在大范圍的產(chǎn)量預(yù)測方面具有一定的局限性，且大范圍的田間調(diào)查需要消耗大量的人力物力。隨機(jī)森林等機(jī)器學(xué)習(xí)算法的一個顯著優(yōu)點是在較少的假設(shè)下，可以組合能免費獲取的遙感等數(shù)據(jù)，通過信息挖掘較好地實現(xiàn)大范圍的作物產(chǎn)量預(yù)測，過程相對簡單，且具有通用性的潛力，但相比于作物模型，該方法無法表達(dá)各因素對產(chǎn)量影響的具體機(jī)理。同時，本研究發(fā)現(xiàn)，利用不同生長時段的樣本建模，模型的預(yù)測精度不同，表明變量的時段是模型非常重要的影響因素之一。劉峻明等[41]探討了基于隨機(jī)森林算法的河南省冬小麥平均拔節(jié)期到平均抽穗期(3-5月)氣象要素對產(chǎn)量預(yù)測的作用，主要引入氣象特征。本研究主要探討基于冬小麥整個生長期(9月-來年6月)的最佳預(yù)測時間窗的選擇，同時考慮到氣候因子對產(chǎn)量預(yù)測的影響，引入植被指數(shù)等與作物生長相關(guān)的要素，使訓(xùn)練模型的因子選擇更豐富。由于氣候、地形、緯度等差異，作物產(chǎn)量預(yù)測的最佳時間窗口可能存在著空間差異性。此外，農(nóng)戶的管理措施(如作物品種、播種量等)對作物產(chǎn)量也具有較大影響，這些信息應(yīng)該包括在未來的建模方法中，即將管理投入與模型結(jié)合是未來研究的一個有前景的途徑。

4 結(jié) 論

選擇2001-2013年河南省冬小麥八個時間段數(shù)據(jù)作為訓(xùn)練樣本，應(yīng)用隨機(jī)森林算法建立基于訓(xùn)練樣本的估產(chǎn)模型，進(jìn)而預(yù)測2014和2015年產(chǎn)量，其中12-3月的估算精度在兩年都最佳，預(yù)測產(chǎn)量與實際產(chǎn)量在空間分布上也基本一致，該時段可作為隨機(jī)森林算法估算河南省冬小麥產(chǎn)量選擇訓(xùn)練樣本的最佳時間；在影響因子中，月降水對冬小麥產(chǎn)量的重要性最大，其次是月最低溫度和EVI，NDVI、干旱指數(shù)、土壤水分和月最高溫度重要性相當(dāng)。