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

?

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

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

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

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

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

0 引 言

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

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

1 材料與方法

1.1 研究區(qū)域

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

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

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

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

1.3 研究方法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3）構(gòu)建過程

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

1.3.2 變量重要性分析

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

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

1.4 模型精度評(píng)價(jià)與驗(yàn)證

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

2 結(jié)果與分析

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

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

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

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

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

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

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

2.3 不同穗分化期對(duì)模型的影響分析

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

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

2.4 模型誤差空間分布

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

2.5 特征變量重要性分析

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

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

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

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

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

3 討 論

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

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

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

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

4 結(jié) 論

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

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

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

3）利用袋外數(shù)據(jù)變量重要性得出除了空間特征之外，冬小麥小花分化期的平均溫度和最低溫度、拔節(jié)至抽穗期間的負(fù)積溫、抽穗開花期的最高溫度4個(gè)氣象指標(biāo)對(duì)產(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] 王靜，李新. 基于作物生長(zhǎng)模型和多源數(shù)據(jù)的融合技術(shù)研究進(jìn)展[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模型的中國(guó)主產(chǎn)區(qū)冬小麥生長(zhǎng)過程動(dòng)態(tài)模擬[J]. 農(nóng)業(yè)工程學(xué)報(bào)，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í)相遙感數(shù)據(jù)的冬小麥估產(chǎn)[J]. 農(nóng)業(yè)工程學(xué)報(bào)，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ù)測(cè)方法研究[J]. 中國(guó)農(nóng)學(xué)通報(bào)，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ī)森林方法研究綜述[J]. 統(tǒng)計(jì)與信息論壇，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)計(jì)局. 河南農(nóng)村統(tǒng)計(jì)年鑒1990－2015[M]. 北京：中國(guó)統(tǒng)計(jì)出版社，2016.

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

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

[22] 房世波. 分離趨勢(shì)產(chǎn)量和氣候產(chǎn)量的方法探討[J].自然災(zāi)害學(xué)報(bào)，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)量對(duì)氣象變化的敏感性分析[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] 成林，李彤霄，劉榮花. 主要生育期氣象變化對(duì)河南省冬小麥生長(zhǎng)及產(chǎn)量的影響[J]. 中國(guó)生態(tài)農(nóng)業(yè)學(xué)報(bào)，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] 馮玉香，何維勛，孫忠富，等. 我國(guó)冬小麥霜凍害的氣候分析[J]. 作物學(xué)報(bào)，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] 張雪芬，鄭有飛，王春乙，等. 冬小麥晚霜凍害時(shí)空分布與多時(shí)間尺度變化規(guī)律分析[J]. 氣象學(xué)報(bào)，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]. 中國(guó)農(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

國(guó)家自然科學(xué)基金項(xiàng)目（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

劉峻明，和曉彤，王鵬新，黃健熙. 長(zhǎng)時(shí)間序列氣象數(shù)據(jù)結(jié)合隨機(jī)森林法早期預(yù)測(cè)冬小麥產(chǎn)量[J]. 農(nóng)業(yè)工程學(xué)報(bào)，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