林 瀅，邵懷勇

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

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

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

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

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

1.1 研究區(qū)概況

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

1.2 數(shù)據(jù)來(lái)源

試驗(yàn)獲取的數(shù)據(jù)包括2001-2015年的遙感、氣候、土壤和產(chǎn)量數(shù)據(jù)。所有數(shù)據(jù)空間分辨率統(tǒng)一為1 km×1 km，時(shí)間分辨率統(tǒng)一為一個(gè)月，并且所有變量將基于小麥生長(zhǎ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ù)來(lái)自于MOD13Q1，周期為16 d，空間分辨率為250 m×250 m。

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

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

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

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

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

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

1.3 研究方法

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

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

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

(1)

式中，F(xiàn)(x)是一個(gè)組合模型，hi是單一決策樹，Y表示輸出變量，I表示特征函數(shù)。試驗(yàn)使用的樹數(shù)量為100。數(shù)據(jù)不被提取的概率為1-1/N，收斂為1/e≈0.368，即有約37%的訓(xùn)練數(shù)據(jù)不會(huì)被運(yùn)用到單個(gè)模型構(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 精度評(píng)估指標(biāo)

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

(5)

(6)

(7)

(8)

1.3.3 重要性評(píng)價(jià)

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

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

2 結(jié)果與分析

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

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

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

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

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

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

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

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

2.3 影響因子重要性評(píng)價(jià)

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

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

3 討 論

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

4 結(jié) 論

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