長時間序列氣象數(shù)據(jù)結(jié)合隨機森林法早期預(yù)測冬小麥產(chǎn)量

劉峻明，和曉彤，王鵬新，黃健熙

?

長時間序列氣象數(shù)據(jù)結(jié)合隨機森林法早期預(yù)測冬小麥產(chǎn)量

劉峻明1,3，和曉彤1,3，王鵬新2,3，黃健熙1

（1. 中國農(nóng)業(yè)大學(xué)土地科學(xué)與技術(shù)學(xué)院，北京 100083；2. 中國農(nóng)業(yè)大學(xué)信息與電氣工程學(xué)院，北京 100083；3. 農(nóng)業(yè)部農(nóng)業(yè)災(zāi)害遙感重點實驗室，北京 100083）

冬小麥生育早期的產(chǎn)量預(yù)測對于制定冬小麥整個生長期的精準(zhǔn)管理決策具有重要參考意義。該文基于隨機森林算法，采用1990—2015年河南省小麥平均拔節(jié)期至平均抽穗期地面觀測氣象數(shù)據(jù)與統(tǒng)計產(chǎn)量數(shù)據(jù)，分別提取不同穗分化期的溫濕度、降水等47個氣象要素和小麥種植區(qū)經(jīng)緯度、高程3個空間要素，共計50個參數(shù)作為特征變量集，以實際單產(chǎn)、氣象產(chǎn)量和相對氣象產(chǎn)量分別作為目標(biāo)變量，構(gòu)建多種變量組合模型對冬小麥產(chǎn)量進行回歸預(yù)測，并結(jié)合袋外數(shù)據(jù)重要性結(jié)果對產(chǎn)量影響因子進行分析。研究結(jié)果表明：1）使用氣象產(chǎn)量和相對氣象產(chǎn)量作為目標(biāo)變量建模的預(yù)測效果優(yōu)于單產(chǎn)模型，決定系數(shù)2均達到0.8以上，氣象產(chǎn)量的平均絕對誤差（mean absolute error，MAE）和均方根誤差（root mean square error，RMSE）分別為415和558 kg/hm2，相對氣象產(chǎn)量的MAE和RMSE分別為0.07和0.09；2）相較于氣象特征，空間特征在產(chǎn)量預(yù)測中起決定性作用，且小花分化期以及抽穗開花期的氣象特征產(chǎn)量預(yù)測精度高于其他穗分化期；3）在氣象特征中，利用袋外數(shù)據(jù)變量重要性得出平均溫度、最低溫度、負(fù)積溫、最高溫度在不同生育階段對產(chǎn)量的影響程度。該研究結(jié)果為冬小麥生育早期產(chǎn)量預(yù)測提供了新的思路和方法。

模型；溫度；隨機森林；產(chǎn)量預(yù)測；冬小麥

0 引 言

冬小麥的產(chǎn)量受生產(chǎn)技術(shù)水平、氣象條件等多重影響，其生長環(huán)境是非常復(fù)雜的非線性系統(tǒng)，產(chǎn)量可看成是一段時期內(nèi)溫度、降水量、生長發(fā)育狀況等多個影響因子相互疊加的結(jié)果。目前關(guān)于農(nóng)作物產(chǎn)量預(yù)測的研究方法多以作物生長機理模型或經(jīng)驗?zāi)Ｐ蜑橹鱗1]。作物生長機理模型通過輸入研究區(qū)作物、氣象、土壤和栽培措施等參數(shù)，對作物生長及產(chǎn)量形成過程進行動態(tài)模擬，適用于點尺度的模擬和預(yù)測，基于數(shù)據(jù)同化方法可以將作物生長機理模型與遙感數(shù)據(jù)進行同化，實現(xiàn)面尺度的作物產(chǎn)量預(yù)測，但由于參數(shù)眾多，誤差累積往往導(dǎo)致精度不夠高[2-4]；經(jīng)驗?zāi)Ｐ屯ㄟ^構(gòu)建作物影響因子與單產(chǎn)之間的統(tǒng)計關(guān)系進行估產(chǎn)，如利用衛(wèi)星數(shù)據(jù)反演的氣象參數(shù)[5]、遙感植被指數(shù)[6]等，建立與其農(nóng)作物產(chǎn)量之間的線性回歸模型，但通常存在經(jīng)驗特征強、容易過擬合的缺點。隨著計算機技術(shù)的發(fā)展，近些年來機器學(xué)習(xí)方法也被應(yīng)用到作物產(chǎn)量預(yù)測研究中，取得較好的效果。黎銳等[7]利用多時相遙感數(shù)據(jù)和支持向量機（support vector machine，SVM）對冬小麥進行估產(chǎn)，該方法核函數(shù)的選取以及參數(shù)的確定具有經(jīng)驗性，對精度有所影響[8]；姜新等[9-10]分別建立了基于葉面積指數(shù)、生物量等作物參數(shù)以及農(nóng)機水平的人工神經(jīng)網(wǎng)絡(luò)（artificial neural network，ANN）產(chǎn)量預(yù)測模型，但其訓(xùn)練過程存在收斂速度慢和易陷入局部極值的問題，且模型參數(shù)較為復(fù)雜[11-12]。隨機森林（random forest，RF）是一種基于分類與回歸樹的機器學(xué)習(xí)算法，由Breiman于2001年提出[13]，相較于傳統(tǒng)決策樹構(gòu)建方法，其優(yōu)越性體現(xiàn)在同等運算率下的高預(yù)測精度，以及在處理多維特征上對多重共線性并不敏感的特性[14]，目前在農(nóng)作物產(chǎn)量預(yù)測方面，也取得了優(yōu)于SVM、BP神經(jīng)網(wǎng)絡(luò)等算法的性能，且模型構(gòu)建過程相對簡單[15-16]；Everingham等[17]基于高光譜數(shù)據(jù)，認(rèn)為RF能較好地進行作物產(chǎn)量預(yù)測；Tulbure等[18]利用RF識別了影響柳枝稷產(chǎn)量的包括降水、土壤肥料等多種關(guān)鍵因子。上述研究雖達到較好預(yù)測效果，但在特征選擇方面欠缺深入研究以及較少考慮多時期氣象特征對作物的影響。

因此本文基于隨機森林算法，以河南省為研究區(qū)域，結(jié)合地面觀測氣象數(shù)據(jù)與冬小麥實際單產(chǎn)數(shù)據(jù)，選擇不同特征及目標(biāo)變量構(gòu)建冬小麥產(chǎn)量預(yù)測模型，并對模型預(yù)測結(jié)果及氣象特征影響進行分析，以期為區(qū)域冬小麥產(chǎn)量預(yù)測提供服務(wù)。

1 材料與方法

1.1 研究區(qū)域

河南省位于31°23′～36°22′N，110°22′～116°38′E，屬北亞熱帶濕潤氣象與暖溫帶半濕潤季風(fēng)氣象的過渡氣象，年均降水量為500～1 000 mm，降水季節(jié)分布不均，全年50%的降水集中在夏季。河南省冬小麥多種植冬性、弱冬性品種，一般在9月中下旬至10月上旬播種，12月中旬后進入越冬期，翌年2月下旬開始進入返青期，3月中下旬進入拔節(jié)期，4月中下旬進入抽穗期，5月底至6月初成熟。

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

河南省冬小麥單產(chǎn)資料來自《河南省統(tǒng)計年鑒》[19]中以縣級行政區(qū)為單元的冬小麥單產(chǎn)數(shù)據(jù)，在1990—2015年間連續(xù)種植冬小麥的縣市共106個，空間位置信息來自縣市內(nèi)氣象觀測站，分布情況如圖1，剔除各縣市單產(chǎn)缺失年份，共獲得2 740個有效單產(chǎn)數(shù)據(jù)，將這些數(shù)據(jù)作為本研究的單產(chǎn)樣本數(shù)據(jù)。氣象資料來自中國氣象科學(xué)數(shù)據(jù)共享服務(wù)網(wǎng)的《中國地面氣象資料數(shù)據(jù)集（V3.0）》，提取1990—2015共26 a的逐日氣象觀測資料，對于沒有氣象資料的縣市，采用地理信息系統(tǒng)的插值方法生成相應(yīng)的氣象數(shù)據(jù)。

圖1 河南省冬小麥種植區(qū)分布Fig.1 Distribution of winter wheat planting area in Henan Province

1.3 研究方法

1.3.1 隨機森林回歸模型的構(gòu)建

1）特征集及目標(biāo)變量構(gòu)造

由于拔節(jié)至抽穗期間是冬小麥發(fā)育最為敏感的階段，是影響產(chǎn)量最重要的時期[20]。將該時段內(nèi)氣象、物候、地理位置作為產(chǎn)量影響因子。隨機森林特征變量集如下：

首先，考慮總研究時段內(nèi)的氣象要素，提取河南省每年平均拔節(jié)至平均抽穗期內(nèi)（3月14日—5月8日）的最高氣溫（max）、最低氣溫（min）、累計降水量（P）、負(fù)積溫（AT）、有效積溫（AT），總計5個特征。

其次，為體現(xiàn)不同時間階段氣象條件的影響，根據(jù)小麥穗分化進程[21]，將拔節(jié)至抽穗期按8 d為單位進行分段，得到7個時間段，與穗分化期大致對應(yīng)關(guān)系為：小花分化前期（3月14日—3月21日）、小花分化后期（3月22日—3月29日)、雌雄蕊分化期（3月30日—4月6日)、藥隔前期（4月7日—4月14日)、藥隔后期（4月15日—4月22日)、四分體時期（4月23日—4月30日)、抽穗開花期（5月1日—5月8日)。生成各穗分化期內(nèi)的max、平均氣溫（avg）、min、最大氣溫日較差（T）、平均相對濕度（hu）、P，總計42個特征。

最后，采用冬小麥種植區(qū)縣市內(nèi)氣象觀測站的經(jīng)緯度（onat）和高程數(shù)據(jù)（le）共計3個變量作為空間特征。

農(nóng)作物的產(chǎn)量可以劃分為趨勢產(chǎn)量、氣象產(chǎn)量和誤差部分。趨勢產(chǎn)量用于表達因技術(shù)革新或社會進步因素而形成的產(chǎn)量，氣象產(chǎn)量用于表示氣象因子變化所帶來的短期波動的產(chǎn)量，其他因素導(dǎo)致的表示為誤差部分[22]。計算公式如下：

式中為實際產(chǎn)量，kg/hm2，y為趨勢產(chǎn)量，kg/hm2，y為氣象產(chǎn)量，kg/hm2，為誤差部分。

趨勢產(chǎn)量對氣象產(chǎn)量存在一種平滑的作用，可以看作是以時間為自變量而進行的線性或非線性模擬。本文假設(shè)農(nóng)業(yè)技術(shù)的提高對作物產(chǎn)量的影響呈平穩(wěn)變化，忽略誤差因素，對單產(chǎn)時間序列采用5a滑動平均法計算得到趨勢產(chǎn)量，對于起始年份趨勢值，依次倒推往年相對應(yīng)的產(chǎn)量數(shù)據(jù)，河南省各年單產(chǎn)均值和趨勢產(chǎn)量曲線如圖2。為消除年代間的農(nóng)業(yè)生產(chǎn)水平差異，利用各年各縣市實際單產(chǎn)減去趨勢產(chǎn)量得到氣象產(chǎn)量，同時利用氣象產(chǎn)量除以趨勢產(chǎn)量獲得相對氣象產(chǎn)量以反映年際間氣象差異造成的產(chǎn)量波動。分別使用單產(chǎn)、氣象產(chǎn)量以及相對氣象產(chǎn)量作為隨機森林模型的目標(biāo)變量，構(gòu)建3種產(chǎn)量預(yù)測模型。公式如下：

式中y為相對氣象產(chǎn)量。

圖2 河南省1990—2015年冬小麥趨勢產(chǎn)量

2）特征變量相關(guān)性分析

圖3為1990—2015年間不同穗分化期氣象要素以及空間特征的關(guān)聯(lián)熱圖，黃色表示強相關(guān)，藍色表示弱相關(guān)。圖3a為拔節(jié)至抽穗期總時段內(nèi)空間特征（onatle），氣象特征（maxminPATAT）之間的相關(guān)性，可以看出，on和le之間相關(guān)性較強，這是由河南省西高東低的地勢造成的；由于氣象因子年際差異較大，經(jīng)緯度、高程與各氣象要素之間的相關(guān)性均較弱；圖3b～圖3h分別為小花分化前期、小花分化后期、雌雄蕊分化期、藥隔前期、藥隔后期、四分體時期和抽穗開花期各氣象要素之間的相關(guān)性，可以看出，從拔節(jié)期開始max和avg以及avg和min之間在各時段內(nèi)存在較強相關(guān)性，但max和min的相關(guān)性較弱；T、hu和P在各時段內(nèi)的相關(guān)性均較弱；相比其他時期，藥隔后期的氣象特征之間相關(guān)性較高。因此，最終選取onat、le作為空間特征變量，不同時期的maxavgminT、PATAT作為氣象特征變量。

注：Lon、Lat、Ele、Tmax、Tavg、Tmin、Td、Ps、Rhu、NAT、AAT、yw分別表示經(jīng)度、緯度、高程、最高氣溫、平均氣溫、最低氣溫、最大氣溫日較差、累計降水量、平均相對濕度、負(fù)積溫、有效積溫、相對氣象產(chǎn)量。

3）構(gòu)建過程

隨機森林是由多棵分類回歸樹（classification and regression tree，CART）構(gòu)成的組合分類模型[13]，各年份各縣市的特征數(shù)據(jù)和產(chǎn)量數(shù)據(jù)進行集成共同構(gòu)成隨機森林的樣本數(shù)據(jù)集，通過自助法（bootstrap）從原始樣本集采樣得到構(gòu)建棵樹所需的個子集，每次未被抽到的數(shù)據(jù)稱為袋外數(shù)據(jù)（out-of-bag，OOB），用來進行內(nèi)部誤差估計和變量重要性評價；生成每棵樹時，從規(guī)模為的特征變量集中隨機選擇個變量（），對于回歸，采用均方差作為節(jié)點分裂標(biāo)準(zhǔn)，遞歸執(zhí)行選取最優(yōu)分枝的操作。由于隨機森林采用樣本和特征的雙重隨機抽樣構(gòu)建決策樹，因此即使不對決策樹進行剪枝操作也不會出現(xiàn)傳統(tǒng)CART決策樹過擬合的現(xiàn)象[14]。最后將這些樹的結(jié)果取平均值即為目標(biāo)變量的預(yù)測值。在RF算法中，根據(jù)試驗結(jié)果實時優(yōu)化決策樹數(shù)目和創(chuàng)建分枝所需特征變量個數(shù)這2個參數(shù)。

1.3.2 變量重要性分析

隨機森林可以解釋若干自變量對因變量的作用。通過模型內(nèi)部重要性結(jié)果，分析不同特征變量對產(chǎn)量的影響程度。其基本思想是：對于變量V，首先計算每棵樹相應(yīng)的袋外數(shù)據(jù)OOB的誤差率EROOB；然后，對袋外數(shù)據(jù)中的變量V值進行隨機的序列改變，而其他所有變量在保持不變的情況下，重新計算改變順序后的袋外數(shù)據(jù)OOB的誤差率EROOB，通過分析袋外數(shù)據(jù)序列改變時誤差的增加情況來估計某一特征變量的重要程度[13]。變量V引起袋外誤差增加的越大，精度減少的越多，說明該變量越重要。變量V的重要性表示為

式中為隨機森林算法中樹的數(shù)量，為個特征中變量的位置。

1.4 模型精度評價與驗證

采用以下3個指標(biāo)作為評價模型擬合程度優(yōu)劣，即決定系數(shù)（coefficient of determination，2）、均方根誤差（root mean square error，RMSE）和平均絕對誤差（mean absolute error，MAE）[23-24]。

2 結(jié)果與分析

2.1 不同目標(biāo)變量驗證結(jié)果對比

以1990—2009年2 107組數(shù)據(jù)作為訓(xùn)練樣本，使用全部特征變量，分別以單產(chǎn)、氣象產(chǎn)量、相對氣象產(chǎn)量為目標(biāo)變量構(gòu)建隨機森林產(chǎn)量預(yù)測模型，根據(jù)經(jīng)驗及多次試驗，將RF算法的參數(shù)和分別設(shè)為500和15。以2010—2015年633組數(shù)據(jù)作為驗證樣本輸入各個預(yù)測模型，結(jié)果如圖4。從圖中可以看出，單產(chǎn)模型的2達到0.71，MAE和RMSE分別達到1 213、1 387 kg/hm2，樣本分布較為松散，預(yù)測結(jié)果整體上比實測值要低；而氣象產(chǎn)量和相對氣象產(chǎn)量的預(yù)測效果顯著優(yōu)于單產(chǎn)模型，趨勢線與1:1線的交叉點均在0刻度附近，大部分樣本聚集在1:1線周圍，有較高的擬合度，決定系數(shù)2均達到0.8以上，氣象產(chǎn)量的MAE和RMSE分別達到415、558 kg/hm2，相對氣象產(chǎn)量的MAE和RMSE分別達到0.07和0.09，相較于單產(chǎn)模型表現(xiàn)出更小的偏差。

產(chǎn)量總體趨勢是增長的，這與品種改良、技術(shù)進步以及田間管理等措施有關(guān)。在沒有消除趨勢產(chǎn)量的情況下，基于1990—2009年數(shù)據(jù)構(gòu)建的隨機森林單產(chǎn)模型中的最大值會較大概率低于后面的年份，導(dǎo)致預(yù)測結(jié)果整體偏低；氣象產(chǎn)量和相對氣象產(chǎn)量的擬合程度較高則說明二者均能較好地去除產(chǎn)量年際間生產(chǎn)技術(shù)水平的影響，且相對氣象產(chǎn)量的趨勢線與1:1線最為貼合，說明在氣象差異的影響下，相對氣象產(chǎn)量相較于氣象產(chǎn)量更能突出空間差異造成的影響，利用隨機森林能達到較好的預(yù)測效果。

注：*，P<0.05；**，P<0.01；下同。

2.2 不同特征變量類型對模型的影響分析

分別使用47個氣象特征、3個空間特征以及使用全部特征變量作為特征集，相對氣象產(chǎn)量作為目標(biāo)變量，RF算法的參數(shù)保持不變設(shè)為500，分別設(shè)為13、1和15。利用2010—2015年河南省產(chǎn)量數(shù)據(jù)進行驗證，結(jié)果見圖5。由圖可知，僅使用氣象特征預(yù)測結(jié)果較差，趨勢線較為平緩；僅使用空間特征預(yù)測精度雖然有所提高，但由于預(yù)測結(jié)果僅取決于地理位置，樣本點沿1:1線呈水平分布，利用RF算法在相同地區(qū)不同的產(chǎn)量數(shù)據(jù)會被預(yù)測為相同的值，遠離1:1線的樣本點可以推斷為受到了氣象要素的影響；而在同時使用空間特征和氣象特征對產(chǎn)量進行預(yù)測時，由圖5c所示，大部分樣本點聚集程度較高，相比前者預(yù)測結(jié)果表現(xiàn)出更小的偏差，相同地區(qū)不同的相對氣象產(chǎn)量因為氣象要素的加入更貼近于真值?？梢?，在構(gòu)建隨機森林對冬小麥產(chǎn)量進行預(yù)測時，需要考慮空間要素對產(chǎn)量造成的影響。

圖5 不同特征變量RF預(yù)測的相對氣象產(chǎn)量的結(jié)果對比

2.3 不同穗分化期對模型的影響分析

為探究不同時段的冬小麥預(yù)測精度，分別將小花分化前期等的氣象要素和空間特征作為特征集，對相對氣象產(chǎn)量進行建模，結(jié)果如圖。由圖可知，基于冬小麥小花分化后期以及抽穗開花期的估產(chǎn)精度比其他時期高，2達到0.8以上。小花分化期是決定穗數(shù)和粒數(shù)的關(guān)鍵時期，作物對氣象變化敏感，而抽穗期是冬小麥將有機物從營養(yǎng)器官轉(zhuǎn)移到籽粒的階段，氣象要素與小麥千粒質(zhì)量密切相關(guān)，這2個穗分化期的氣象要素對產(chǎn)量影響最大，故估產(chǎn)精度最高。藥隔后期的產(chǎn)量預(yù)測值與實際值的偏差最大，這是由于相比其他時期，該時期內(nèi)氣象要素之間存在較強相關(guān)性，存在多重共線性，削弱了空間特征的影響所造成的。綜上，利用小花分化后期和抽穗開花期的氣象條件通過RF算法可達到較好的產(chǎn)量預(yù)測效果。

圖6 不同穗分化期RF預(yù)測結(jié)果對比

2.4 模型誤差空間分布

將模型參數(shù)和分別設(shè)為500和15，使用全部特征變量，以相對氣象產(chǎn)量為目標(biāo)變量構(gòu)建隨機森林模型，然后對2014和2015年河南省小麥種植區(qū)分別進行驗證，結(jié)果如圖7。根據(jù)河南省各縣市統(tǒng)計數(shù)據(jù)[19]，冬小麥單產(chǎn)常年東南部偏高，西北部偏低。究其原因，主要是因為在緯度高的地區(qū)，太陽高度角比較小，存在斜射現(xiàn)象，單位面積的地表獲得的太陽輻射少，作物物候期較晚，實際產(chǎn)量較低；緯度較低的區(qū)域物候期較早，實際產(chǎn)量較高?？紤]到河南省地勢基本上是西高東低，西部山區(qū)溫度相對較低，不能滿足冬小麥正常生長發(fā)育所需要的活動積溫，對小麥高產(chǎn)穩(wěn)產(chǎn)有一定影響。從圖7中可以看出，相對氣象產(chǎn)量的相對誤差均在±0.2內(nèi)，預(yù)測結(jié)果整體上呈現(xiàn)東南部偏低，西北部偏高的分布趨勢。2014—2015年，西部地區(qū)的三門峽市、洛陽市和北部的安陽市冬小麥單產(chǎn)分別低于河南省平均產(chǎn)量26%～40%左右，而模型預(yù)測相對誤差平均在0.15左右，東南部地區(qū)的商丘市、周口市和駐馬店市的小麥單產(chǎn)分別高于全省平均產(chǎn)量的18%～22%左右，常年高于全省均值的20%左右，模型預(yù)測相對誤差平均在-0.18左右。顯然，模型預(yù)測結(jié)果在實際產(chǎn)量較低時易被高估，在實際產(chǎn)量較高時易被低估，這是由于RF算法預(yù)測結(jié)果為多棵決策樹投票得到的結(jié)果所致，算法本身傾向于數(shù)據(jù)的平均狀態(tài)[25]。

2.5 特征變量重要性分析

利用1990—2015年的冬小麥特征數(shù)據(jù)與相對氣象產(chǎn)量分別進行基于RF的OOB重要性和相關(guān)性分析，按照OOB誤差對特征變量進行由大到小排序，結(jié)果見表1。重要性排名前3為均為空間特征，依次為緯度（at）、高程（le）和經(jīng)度（on），且他們的||也很顯著，說明空間位置對相對氣象產(chǎn)量的影響起著主導(dǎo)地位，決定了其本底數(shù)值。第4～7位分別為小花分化后期的平均溫度（avg）小花分化前期的最低溫（min）、拔節(jié)至抽穗期總時段的負(fù)積溫（AT）和抽穗開花期的最高溫度（max），代表4個不同的氣象指標(biāo)，平均溫度代表了小麥基本的熱量需求，最低溫度、負(fù)積溫和最高溫度則分別代表了極端低溫條件、持續(xù)低溫累積、極端高溫條件對小麥的脅迫影響，因為負(fù)積溫主要是在初級累積，這個結(jié)果也反映了該地區(qū)初期易遭受晚霜凍害、后期易受干熱風(fēng)影響的氣候特點，這與以往研究結(jié)果[26-27]相一致。

從||看，除了空間位置的重要性與||基本一致外，氣象要素的重要性與||的關(guān)系沒呈現(xiàn)出明顯規(guī)律性。如各階段最大氣溫日較差（T）的||在0.22～0.29之間，但其重要性并未體現(xiàn)出來，而各階段平均氣溫（avg）的||更只有0.01～0.09，但小花分化后期的平均溫度（avg）的重要性相對較為顯著，表明單因素氣象特征與相對氣象產(chǎn)量的相關(guān)性不太具有意義。

圖7 2014—2015年河南省相對氣象產(chǎn)量相對誤差分布圖

河南省是晚霜凍害發(fā)生較為頻繁的地區(qū)[28]，國內(nèi)學(xué)者多以最低氣溫作為該研究區(qū)的氣象指標(biāo)，張雪芬等[29]發(fā)現(xiàn)低溫發(fā)生頻率與產(chǎn)量之間相關(guān)性顯著，但產(chǎn)量受諸多因素影響[30]，低溫并不必然導(dǎo)致災(zāi)害發(fā)生，隨著品種和耕作措施的改進，即使發(fā)生低溫，后期仍可能獲得豐產(chǎn)。這可解釋了表1中氣象特征重要性和相關(guān)性不一致的問題，利用單一氣象指標(biāo)建立回歸分析，預(yù)測結(jié)果可能會存在較大偏差。隨機森林重要性的概念反映的是多因素疊加作用下，該因素的影響程度，并將這種疊加效果在宏觀層面體現(xiàn)出來，其重要優(yōu)勢之一是能處理具有多維特征的數(shù)據(jù)，并且不用做特征篩選，這有助于從更綜合的層面上來分析多個因素的疊加影響[31]。

表1 特征變量重要性及變量與相對氣象產(chǎn)量之間相關(guān)性分析

3 討 論

受地勢影響，河南省冬小麥種植面積常年東多西少，單產(chǎn)東南高西北低[19]，呈現(xiàn)顯著的空間分布特征?？臻g位置主要影響作物所獲得的積溫和物候，但由于氣象條件年際差異大以及樣本時間周期不是很長，導(dǎo)致空間位置與氣象環(huán)境因素之間的相關(guān)性較弱，因此在構(gòu)建隨機森林的過程中，不能僅考慮氣象要素特征，空間差異性是需要考慮的重要因素。

理論上氣象產(chǎn)量和相對氣象產(chǎn)量僅受氣象條件影響，但因它們的值不能直接通過測量獲取，所以高度依賴于去趨勢方法的選擇。本文基于5a滑動平均方法對單產(chǎn)進行了去趨勢而得到氣象產(chǎn)量和相對氣象產(chǎn)量，雖然可在一定程度上消除單產(chǎn)的年際影響，能突出氣象要素對產(chǎn)量的影響，該方法具有一定的局限性，可以考慮其他去趨勢方法，盡可能僅保留氣象特征對產(chǎn)量的影響，減小影響因素的干擾。

在時間段劃分方面，本文根據(jù)河南省平均拔節(jié)時間和平均抽穗時間，將該時段均勻劃分，每年采用相同的時間節(jié)點作為小麥穗分化期的劃分依據(jù)，但小麥穗分化進程受到溫度、水分以及品種特性在內(nèi)的諸多因素的影響[32-33]，因此，為了進一步提高模型的預(yù)測精度，可考慮氣象變化的年際差異，更為準(zhǔn)確地劃分時間段。

由于總體樣本數(shù)量偏少，災(zāi)害樣本更少，本文未將災(zāi)害年和非災(zāi)害年分開建模，基于較少樣本訓(xùn)練建立的模型可能會產(chǎn)生較大偏離。本文選用1990—2009年用于訓(xùn)練，2010—2015年用于驗證，帶有一定的主觀性，但考慮到2010—2015年間，既有比較正常的樣本，也有災(zāi)害樣本，既有比較正常的年份，也有災(zāi)害年份，所以這段時間作為驗證年份具有合理性。

4 結(jié) 論

基于地面觀測氣象數(shù)據(jù)、空間特征與冬小麥實際單產(chǎn)數(shù)據(jù)，利用隨機森林算法對冬小麥產(chǎn)量進行了回歸試驗，對預(yù)測結(jié)果的空間分布狀況進行了分析，并結(jié)合袋外數(shù)據(jù)重要性對產(chǎn)量影響因子進行了探討。結(jié)果如下：

1）隨機森林算法在預(yù)測產(chǎn)量上具有很大潛力，使用氣象產(chǎn)量和相對氣象產(chǎn)量作為目標(biāo)變量建模的預(yù)測效果優(yōu)于單產(chǎn)模型，相對氣象產(chǎn)量的預(yù)測效果最優(yōu)，其決定系數(shù)2達到0.84，平均絕對誤差（mean absolute error，MAE）和均方根誤差（root mean square error，RMSE）分別達到0.07和0.09左右。

2）空間特征在所構(gòu)建的隨機森林產(chǎn)量預(yù)測模型中起到了重要的作用，在此基礎(chǔ)上增加氣象特征可以使預(yù)測結(jié)果表現(xiàn)出更小的偏差，2達到0.88，MAE和RMSE分別達到0.06和0.08，且利用冬小麥小花分化期以及抽穗開花期的氣象特征進行估產(chǎn)，精度要高于其他穗分化期，說明該時段環(huán)境的變化對最終產(chǎn)量造成的影響更大。

3）利用袋外數(shù)據(jù)變量重要性得出除了空間特征之外，冬小麥小花分化期的平均溫度和最低溫度、拔節(jié)至抽穗期間的負(fù)積溫、抽穗開花期的最高溫度4個氣象指標(biāo)對產(chǎn)量影響較大。

[1] Capa-Morocho M, Rodríguez-Fonseca Belén, Ruiz-Ramos M . Crop yield as a bioclimatic index of El Ni?o impact in Europe: Crop forecast implications[J]. Agricultural and Forest Meteorology, 2014, 198/199: 42－52.

[2] 王靜，李新. 基于作物生長模型和多源數(shù)據(jù)的融合技術(shù)研究進展[J]. 遙感技術(shù)與應(yīng)用，2015，30(2)：209－219. Wang Jing, Li Xin.Research progress of fusion technology based on crop growth model and multi-source data[J].Remote Sensing Technology and Application, 2015, 30(2): 209－219. (in Chinese with English abstract)

[3] 黃健熙，賈世靈，馬鴻元，等. 基于WOFOST模型的中國主產(chǎn)區(qū)冬小麥生長過程動態(tài)模擬[J]. 農(nóng)業(yè)工程學(xué)報，2017，33(10)：222－228. Huang Jianxi, Jia Shiling, Ma Hongyuan, et al. Dynamic simulation of growth process of winter wheat in main production areas of China based on WOFOST model[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2017, 33(10): 222－228. (in Chinese with English abstract)

[4] Brown J N, Zvi H, Dean H, et al. Seasonal climate forecasts provide more definitive and accurate crop yield predictions[J]. Agricultural and Forest Meteorology, 2018, 260/261: 247－254.

[5] Basist A, Dinar A, Blankespoor B, et al. Use of satellite information on wetness and temperature for crop yield prediction and river resource planning[J]. Climate Smart Agriculture, 2018, 52: 77－104.

[6] Zhang S, Liu L. The potential of the MERIS Terrestrial Chlorophyll Index for crop yield prediction[J]. Remote Sensing Letters, 2014, 5(8):10.

[7] 黎銳，李存軍，徐新剛，等. 基于支持向量回歸（SVR）和多時相遙感數(shù)據(jù)的冬小麥估產(chǎn)[J]. 農(nóng)業(yè)工程學(xué)報，2009，25(7)：114－117. Li Rui, Li Cunjun, Xu Xingang, et al. Winter wheat yieldestimation based on suport vector machine regression andmulti-temporal remote sensing data[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2009, 25(7): 114－117. (in Chinese with English abstract)

[8] Zai Songmei, Jia Yanhui, Wen Ji, et al. Study on grain yield prediction of irrigation district based on least squares support vector machine[J]. Agricultural Science and Technology, 2009, 10(6): 1－3, 6.

[9] 姜新. 河南省糧食產(chǎn)量影響因素和預(yù)測方法研究[J]. 中國農(nóng)學(xué)通報，2019，35(1)：154－158. Jiang Xin. Study on the factors affecting grain yield and forecasting methods in henan province[J]. Chinese Agricultural Science Bulletin, 2019, 35(1): 154－158. (in Chinese with English abstract)

[10] Pandey A, Mishra A. Application of artificial neural networks in yield prediction of potato crop[J]. Russian Agricultural Sciences, 2017, 43(3): 266－272.

[11] Kaul M, Hill R L, Walthall C . Artificial neural networks for corn and soybean yield prediction[J]. Agricultural Systems, 2005, 85(1): 1－18.

[12] Uno Y, Prasher S O, Lacroix R, et al. Artificial neural networks to predict corn yield from compact airborne spectrographic imager data[J]. Computer and Electronics in Agriculture, 2005, 47(2): 149－161.

[13] Breiman L. Random forests[J]. Machine Learning, 2001, 45(1): 5－32.

[14] 方匡南，吳見彬. 隨機森林方法研究綜述[J]. 統(tǒng)計與信息論壇，2011，26(3)：32－38. Fang Kuangnan, Wu Jianbin. A review of random forest method research[J]. Statistics and Information Forum, 2011, 26(3): 32－38.

[15] Anna Chlingaryan, Salah Sukkarieh, Brett Whelan. Machine learning approaches for crop yield prediction and nitrogen status estimation in precision agriculture: A review[J]. Computers and Electronics in Agriculture, 2018, 151: 61－69.

[16] Dhivya Elavarasan, Durai Raj Vincent, Vishal Sharma, et al. Forecasting yield by integrating agrarian factors and machine learning models: A survey[J]. Computers and Electronics in Agriculture, 2018, 155: 257－282.

[17] Everingham Y, Sexton J, Skocaj D, et al. Accurate prediction of sugarcane yield using a random forest algorithm[J]. Agronomy for Sustainable Development, 2016, 36(2): 1－9.

[18] Tulbure M G, Wimberly M C, Boe A, et al. Climatic and genetic controls of yields of switchgrass, a model bioenergy species[J]. Agriculture Ecosystems and Environment, 2012, 146(1): 121－129.

[19] 河南省統(tǒng)計局. 河南農(nóng)村統(tǒng)計年鑒1990－2015[M]. 北京：中國統(tǒng)計出版社，2016.

[20] 王紹中，田云峰，郭天財，等. 河南小麥栽培學(xué)[M]. 北京：中國統(tǒng)計出版社，2011.

[21] 崔金梅，郭天財. 小麥的穗[M]. 北京：中國農(nóng)業(yè)出版社，2006.

[22] 房世波. 分離趨勢產(chǎn)量和氣候產(chǎn)量的方法探討[J].自然災(zāi)害學(xué)報，2011，20(6)：13－18. Fang Shibo. Discussion on methods of separating trend yield and climate yield[J]. Journal of Natural Disasters, 2011, 20(6): 13－18. (in Chinese with English abstract)

[23] Jamieson P D, Porter J R, Wilson D R. A test of the computer simulation model ARCWHEAT1 on wheat crops grown in New Zealand[J]. Field Crops Research, 1991, 27(4): 337－350.

[24] Zhu Y, Li Y, Feng W, et al. Monitoring leaf nitrogen in wheat using canopy reflectance spectra[J]. Canadian Journal of Plant Science, 2006, 86(4): 1037－1046.

[25] Zhou Zhihua, Ji Feng. Deep forest: Towards an alternative to deep neural networks[J]. Machine Learning, 2017, 497: 3553－3559.

[26] 張榮榮，寧曉菊，秦耀辰，等. 1980年以來河南省主要糧食作物產(chǎn)量對氣象變化的敏感性分析[J]. 資源科學(xué)，2018，40(1)：137－149. Zhang Rongrong, Ning Xiaoju, Qin Yaochen, et al. Sensitivity analysis of main grain crop yields to climate change in Henan Province since 1980[J]. Resources Science, 2018, 40(1): 137－149. (in Chinese with English abstract)

[27] 成林，李彤霄，劉榮花. 主要生育期氣象變化對河南省冬小麥生長及產(chǎn)量的影響[J]. 中國生態(tài)農(nóng)業(yè)學(xué)報，2017，25(6)：931－940. Cheng Lin, Li Tongxiao, Liu Ronghua. Effects of climate change during main growth period on winter wheat growth and yield in Henan Province[J].Chinese Journal of Eco-Agriculture, 2017, 25(6): 931－940. (in Chinese with English abstract)

[28] 馮玉香，何維勛，孫忠富，等. 我國冬小麥霜凍害的氣候分析[J]. 作物學(xué)報，1999，25(3)：335－340. Feng Yuxiang, He Weixun, Sun Zhongfu, et al. Climatological study onfrost damage of winter wheat in China[J]. Acta Agronomy Sinica, 1999, 25(3): 335－340. (in Chinese with English abstract)

[29] 張雪芬，鄭有飛，王春乙，等. 冬小麥晚霜凍害時空分布與多時間尺度變化規(guī)律分析[J]. 氣象學(xué)報，2009，67(2)：321－330. Zhang Xuefeng, Zheng Youfei, Wang Chunyi, et al. Spatial-temporal distribution and multiple-temporal scale variationanalyses of winter wheat late freezing injury[J].Acta Meteorologica Sinica, 2009, 67(2): 321－330. (in Chinese with English abstract)

[30] 朱虹暉，武永峰，宋吉青，等. 基于多因子關(guān)聯(lián)的冬小麥晚霜凍害分析：以河南省為例[J]. 中國農(nóng)業(yè)氣象，2018，39(1)：59－68. Zhu Honghui, Wu Yongfeng, Song Jiqing, et al. Analysis of winter frost damage of winter wheat based on multi-factor correlation: A case study of Henan Province[J]. Chinese Journal of Agricultural Meteorology, 2018, 39(1): 59－68. (in Chinese with English abstract)

[31] Gislason P, Benediktsson J. Random forests for land cover classification[J]. Pattern Recognition Letters, 2006, 27(4): 294－300.

[32] Mo X , Liu S , Lin Z , et al. Prediction of crop yield, water consumption and water use efficiency with a SVAT-crop growth model using remotely sensed data on the North China Plain[J]. Ecological Modelling, 2005, 183(2/3):301-322.

[33] Krupnik T J, Ahmed Z U, Timsina J, et al. Untangling crop management and environmental influences on wheat yield variability in Bangladesh: An application of non-parametric approaches[J].Agricultural Systems, 2015, 139: 166－179.

Early prediction of winter wheat yield with long time series meteorological data and random forest method

Liu Junming1,3, He Xiaotong1,3, Wang Pengxin2,3, Huang Jianxi1

(1.,,100083,2.,,100083,3.,,100083,)

Early prediction of winter wheat yield is of great significance for the formulation of precise management decisions for the whole growth period of winter wheat. The yield of winter wheat is affected by production technology level and climatic conditions. This study analyzed the feasibility of early prediction of winter wheat yield with long time series meteorological data and random forest method in Henan Province. Winter wheat was planted in a total of 106 counties (cities) in Henan province. Based on the ground observation meteorological data and the winter wheat statistical yield data from the year of 1990 to 2015, we extracted 47 climatic factors such as temperature, humidity and precipitation in different growth stages from wheat jointing to heading stage, and 3 spatial factors of latitude, longitude and elevation. A total of 50 parameters were used as a set of feature variables. The actual yield, meteorological yield and relative meteorological yield were used as the target variables respectively, and a random forest yield prediction model with multiple variables was constructed. The data from the year of 1990 to 2009 were used as training samples to construct the model and the forests constructed were validated with data from the year of 2010 to 2015. The yield impact factors were analyzed by combining the data importance results outside the bag. The results showed that: 1) The prediction results by using meteorological yield and relative meteorological yield as the target variables were better than the yield model. For the meteorological yield and relative meteorological yield models, the values of determination coefficient2were both above 0.8, the values of mean absolute error (MAE) and root mean square error (RMSE) of meteorological yield were 415 and 558 kg/hm2, respectively, and the values of MAE and RMSE of relative meteorological yield were 0.07 and 0.09, respectively; 2) The spatial characteristics played an important role in the improving the random forest yield model. However, if the model included only spatial parameters, the predicted values were horizontally distributed along 1:1 line and the different yields in the same region by using random forest algorithm were predicted as the same values. The values far from 1:1 line might be affected by meteorological factors. Therefore, on this basis, adding meteorological features improved the prediction accuracy with smaller deviations, higher2(0.88), and smaller MAE and RMSE (0.06 and 0.08). 3) The model prediction was also affected by crop growing stages. The accuracy based on the meteorological features of winter wheat florets differentiation and heading and flowering stage was higher than the other spike differentiation periods, indicating that the environmental changes during this period have a greater impact on the final yield; The predicted results at the late drug interval had the larger deviation from the actual yield. It was because the meteorological factors had strong correlation and it weakened the impacts of spatial characteristics. 4) Based on the importance of outside the bag data, In the meteorological features, the average temperature and minimum temperature of winter wheat floret differentiation period, the spatial characteristics parameters were important. In addition, the negative accumulated temperature from the jointing to heading stage, and the maximum temperature at heading and flowering stage had great influence on yield. During the model establishment, we didn’t differentiate disaster from non-disaster year because the sample sizes were small. However, during the model validation, the data were from both normal and disaster years, which could ensure the reliability of the prediction model. Thus, the winter wheat yield prediction based on random forest should consider both spatial and meteorological characteristics parameters. The results of this study provide new ideas and methods for early prediction of winter wheat yield.

models; temperature; random forest;yield prediction; winter wheat

2018-09-20

2019-02-28

國家自然科學(xué)基金項目（41471342）

劉峻明，副教授，博士，主要從事地理信息系統(tǒng)和定量遙感研究。Email：liujunming2000@163.com

10.11975/j.issn.1002-6819.2019.06.019

S127；S512.1+1

A

1002-6819(2019)-06-0158-09

劉峻明，和曉彤，王鵬新，黃健熙. 長時間序列氣象數(shù)據(jù)結(jié)合隨機森林法早期預(yù)測冬小麥產(chǎn)量[J]. 農(nóng)業(yè)工程學(xué)報，2019，35(6)：158－166. doi：10.11975/j.issn.1002-6819.2019.06.019 http://www.tcsae.org

Liu Junming He Xiaotong Wang Pengxin, Huang Jianxi. Early prediction of winter wheat yield with long time series meteorological data and random forest method[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2019, 35(6): 158－166. (in Chinese with English abstract) doi：10.11975/j.issn.1002-6819.2019.06.019 http://www.tcsae.org