朱秀芳，李石波，石振君，任笑，程昌秀

（1. 北京師范大學(xué)遙感科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100875；2. 北京師范大學(xué)環(huán)境演變與自然災(zāi)害教育部重點(diǎn)實(shí)驗(yàn)室，北京 100875；3. 中國(guó)地質(zhì)大學(xué)土地科學(xué)技術(shù)學(xué)院，北京 100083；4. 北京師范大學(xué)地理科學(xué)學(xué)部，北京 100875）

從20世紀(jì)上半葉至今，地球大氣中的CO2濃度不斷升高，變暖趨勢(shì)越來(lái)越顯著，全球氣候變化，導(dǎo)致異常氣候事件頻發(fā)[1]。相關(guān)研究表明，盡管農(nóng)業(yè)生產(chǎn)受氣候變化的影響在不同區(qū)域有一定的差異，但整體上弊大于利[2]。我國(guó)農(nóng)業(yè)受氣候變化的影響十分嚴(yán)重[3]，在各類農(nóng)業(yè)氣象災(zāi)害中，干旱和洪澇是影響我國(guó)農(nóng)作物生產(chǎn)最主要的氣象災(zāi)害。1980年以來(lái)，我國(guó)洪旱災(zāi)害頻發(fā)。1998年的長(zhǎng)江和松花江流域特大洪水事件和2016年的長(zhǎng)江中下游及太湖洪澇事件直接帶來(lái)上千億元經(jīng)濟(jì)損失[4]。相比洪澇，干旱對(duì)我國(guó)農(nóng)業(yè)生產(chǎn)的影響更為嚴(yán)重[5]。有關(guān)統(tǒng)計(jì)數(shù)據(jù)表明，2007至2017年全國(guó)農(nóng)作物旱災(zāi)面積年均值為1.53×107hm2，糧食損失年均值為2.14億t，造成直接經(jīng)濟(jì)損失年均值為282.15億元，占GDP均值的0.21%[6]。在全國(guó)各種氣象災(zāi)害中，干旱導(dǎo)致的農(nóng)作物受災(zāi)面積占比約56%，洪澇導(dǎo)致的農(nóng)作物受災(zāi)面積占比約24%[7]。此外，低溫冷害、干熱風(fēng)等氣候事件的頻繁發(fā)生也給我國(guó)農(nóng)作物安全生產(chǎn)帶來(lái)了巨大挑戰(zhàn)。在中國(guó)北方麥區(qū)，干熱風(fēng)災(zāi)害的發(fā)生一般會(huì)導(dǎo)致小麥減產(chǎn)5%～10%，在受災(zāi)嚴(yán)重的年份減產(chǎn)率能達(dá)到20%～30%[8]。

隨著遙感技術(shù)的不斷發(fā)展與創(chuàng)新應(yīng)用，許多專家學(xué)者將遙感技術(shù)應(yīng)用到農(nóng)作物估產(chǎn)中[9]，為作物估產(chǎn)研究開辟了新途徑[10]。從模型建立的理論角度出發(fā)，當(dāng)前主流的作物估產(chǎn)模型可以分為四類：經(jīng)驗(yàn)統(tǒng)計(jì)模型、光能利用率模型、作物生長(zhǎng)模擬模型和耦合模型[11]。探究氣候變化對(duì)作物產(chǎn)量的影響，大多采用經(jīng)驗(yàn)統(tǒng)計(jì)模型中的氣象估產(chǎn)模型。氣象估產(chǎn)模型發(fā)展較早，始于20世紀(jì)70年代末，基本的氣象估產(chǎn)模型有三種。第一種，直接建立氣象因子和作物產(chǎn)量之間的回歸模型；第二種，首先計(jì)算相鄰兩年作物產(chǎn)量差和氣象因子差，然后建立作物產(chǎn)量差和氣象因子差之間的回歸模型，進(jìn)而進(jìn)行作物產(chǎn)量預(yù)測(cè)。第三種，將作物產(chǎn)量分解為趨勢(shì)產(chǎn)量和氣象產(chǎn)量?jī)刹糠?，前者反映由技術(shù)進(jìn)步（如灌溉、施肥、新品種等）導(dǎo)致的產(chǎn)量的長(zhǎng)期變化，后者反映由自然氣候要素（光照、降水、輻射等）引起的產(chǎn)量的短期波動(dòng)。

常用的氣象估產(chǎn)模型構(gòu)建方法包括線性和非線性回歸兩大類[12-13]。隨著數(shù)據(jù)挖掘技術(shù)的發(fā)展，神經(jīng)網(wǎng)絡(luò)、支持向量機(jī)、隨機(jī)森林等新的數(shù)據(jù)挖掘手段也應(yīng)用于估產(chǎn)模型的構(gòu)建[14-16]。Breiman[17]在2001年提出了隨機(jī)森林算法。它是一種集成學(xué)習(xí)算法，具有穩(wěn)定性好，預(yù)測(cè)精度高，不易產(chǎn)生過(guò)擬合等優(yōu)點(diǎn)，不需要顧慮多重共線性的問題，與神經(jīng)網(wǎng)絡(luò)和線性回歸方法相比，其表現(xiàn)更穩(wěn)定，具有較強(qiáng)的抗干擾能力。該算法已經(jīng)廣泛應(yīng)用于干旱監(jiān)測(cè)[18]、山體滑坡預(yù)測(cè)[19]、洪澇[20]和火災(zāi)風(fēng)險(xiǎn)[21]分析等方面。但目前應(yīng)用隨機(jī)森林算法進(jìn)行異常氣候條件下的作物估產(chǎn)的研究甚少。目前大多研究也并沒有驗(yàn)證所建立的估產(chǎn)模型在災(zāi)害年份的估產(chǎn)效果。為此，本文基于隨機(jī)森林算法，以氣象估產(chǎn)模型為基礎(chǔ)，結(jié)合氣象災(zāi)害指數(shù)和遙感植被指數(shù)，在中國(guó)五大小麥種植區(qū)分別建立估產(chǎn)模型，并驗(yàn)證其在典型災(zāi)害年份的估產(chǎn)精度，為探究氣候變化對(duì)小麥產(chǎn)量的影響，以及建立滿足異常氣候下估產(chǎn)精度的模型提供參考。

1 研究區(qū)概況

本文研究范圍（圖1）包括西北春麥區(qū)、北部冬麥區(qū)、黃淮海冬麥區(qū)、長(zhǎng)江中下游冬麥區(qū)與西南冬麥區(qū)[22-23]。西北春麥區(qū)位于黃淮上游三大高原的交匯地帶，該區(qū)由寧夏全區(qū)、甘肅省大部分地區(qū)，內(nèi)蒙古西北部地區(qū)構(gòu)成，冬天寒冷、夏季炎熱，春季與秋季多風(fēng)，天氣較為干燥，白天和晚上溫差很大。降水不足是影響該區(qū)小麥作物生長(zhǎng)的重要因素，部分地區(qū)春小麥生長(zhǎng)后期有干熱風(fēng)危害。北部冬麥區(qū)位于我國(guó)中緯度地帶，由北京全市、天津全市、河北、山西以及陜西省部分地區(qū)組成，屬于暖溫帶季風(fēng)區(qū)域。該區(qū)冬寒春旱，降水不足且分布不均勻。黃淮海冬麥區(qū)位于黃河中下游，由山東省全部、河南省大部分區(qū)域以及陜西省與山西省南部、安徽省與江蘇省北部、河北南部的少部分地區(qū)組成。該區(qū)北部小麥在低溫年份有遭受低溫冷害的風(fēng)險(xiǎn)，南部地區(qū)氣溫較高，冬小麥返青期不明顯。降水量南多北少，東多西少，小麥生育期降雨可基本滿足，但北部仍偶爾有旱災(zāi)發(fā)生，需灌溉。長(zhǎng)江中下游冬麥區(qū)由湖北全省、湖南全省、江西全省、浙江全省、上海全市，安徽省與浙江省大部分地區(qū)組成。該區(qū)大部分地區(qū)的自然條件適宜小麥生長(zhǎng)，但該區(qū)降水極不平衡，南部降水過(guò)多不適宜小麥生長(zhǎng)，北部降水較少時(shí)有干旱發(fā)生。西南冬麥區(qū)由重慶全市、貴州全省、四川省東部與云南省北部組成。該麥區(qū)地形復(fù)雜，包括山地、盆地、平原等多種地形，其中山地為該區(qū)主要地形，復(fù)雜的地形也導(dǎo)致該區(qū)域氣候環(huán)境差異較大，農(nóng)作物種植、成熟時(shí)間差異較大，除平原地區(qū)外，農(nóng)田地塊較為破碎。

圖1 中國(guó)小麥種植區(qū)劃Fig. 1 Wheat growing region of China

2 技術(shù)路線與研究方法

本文技術(shù)路線如圖2所示，主要包括以下步驟：1）利用原始?xì)庀髷?shù)據(jù)與遙感數(shù)據(jù)建立相應(yīng)氣象災(zāi)害指標(biāo)與遙感植被指數(shù)；2）利用直線滑動(dòng)平均法和HP濾波法擬合趨勢(shì)產(chǎn)量，選擇最佳擬合方法；3）以氣象災(zāi)害指標(biāo)、遙感植被指數(shù)為輸入變量，實(shí)際產(chǎn)量為輸出變量，建立樣本庫(kù)；4）利用隨機(jī)森林算法在不同麥區(qū)建立估產(chǎn)模型并進(jìn)行精度驗(yàn)證；5）篩選出各麥區(qū)不同類型災(zāi)害年份，并進(jìn)行估產(chǎn)驗(yàn)證；6）對(duì)不同麥區(qū)輸入變量進(jìn)行重要性評(píng)價(jià)。

圖2 技術(shù)路線Fig. 2 Technical flowchart

2.1 數(shù)據(jù)獲取與預(yù)處理

通過(guò)查閱資料發(fā)現(xiàn)在我國(guó)眾多不同類型氣象災(zāi)害中，對(duì)農(nóng)作物生產(chǎn)影響較大的災(zāi)害類型主要包括旱災(zāi)、洪澇、低溫凍害以及干熱風(fēng)四種[24-28]。為此，本文收集了氣象數(shù)據(jù)與遙感數(shù)據(jù)，用于構(gòu)建表示干旱、洪澇、低溫冷害、干熱風(fēng)四種災(zāi)害的指標(biāo)，以及表示小麥不同時(shí)期長(zhǎng)勢(shì)情況的遙感植被指數(shù)指標(biāo)，用于估產(chǎn)模型的建立。具體包括：1）來(lái)自全球標(biāo)準(zhǔn)化降水蒸散指數(shù)（Standardized Precipitation Evapotranspiration Index, SPEI）數(shù)據(jù)庫(kù)1999—2018年的SPEIbase v2.6數(shù)據(jù)集，時(shí)間尺度為1～48月，分辨率為0.5°；2）來(lái)自歐空局哥白尼氣象數(shù)據(jù)中心的1999—2018年的月平均降水量、月平均溫度、14時(shí)風(fēng)速、14時(shí)相對(duì)濕度、日最高溫度，除了日最高溫度分辨率為0.1°外，其他數(shù)據(jù)的空間分辨率均為0.25°；3）來(lái)自美國(guó)航空航天局的2000—2018年的16 d合成的MOD13Q1-NDVI數(shù)據(jù)，空間分辨率為250 m。

2.2 輸入?yún)?shù)計(jì)算

五大麥區(qū)范圍廣，小麥種植時(shí)間差異較大，各麥區(qū)的小麥生育期間分別為黃淮海冬麥區(qū)9月-次年6月、北部冬麥區(qū)9月-次年6月、西南冬麥區(qū)8月-次年7月、長(zhǎng)江中下游冬麥區(qū)10月-次年5月、西北春麥區(qū)3月-8月，因此本文以各麥區(qū)小麥的生育期并集為研究時(shí)段來(lái)對(duì)參數(shù)進(jìn)行計(jì)算。

2.2.1 干旱指標(biāo) 本文選用SPEI8-6、SPEI6-10、SPEI6-10、SPEI5-8、SPEI7-12分別表示西北春麥區(qū)、北部冬麥區(qū)、黃淮海冬麥區(qū)、長(zhǎng)江中下游冬麥區(qū)、西南冬麥區(qū)中小麥整個(gè)生育期內(nèi)的干旱程度。其中，SPEIm-n為m月份時(shí)間尺度為n的SPEI干旱指數(shù)，用來(lái)描述m月份過(guò)去n個(gè)月整體的干旱情況。具體計(jì)算流程是首先將SPEI干旱指數(shù)重采樣至1 km空間尺度，用全國(guó)耕地范圍數(shù)據(jù)進(jìn)行耕地掩膜，然后以縣域范圍內(nèi)SPEI干旱指數(shù)均值作為描述某縣小麥生育期內(nèi)整體干旱程度的指標(biāo)。

2.2.2 干熱風(fēng)指標(biāo) 我國(guó)小麥干熱風(fēng)具體可以分為高溫低濕型、雨后青枯型和旱風(fēng)型3種[29]。國(guó)家氣象局發(fā)布了2019年干熱風(fēng)國(guó)家標(biāo)準(zhǔn)，采用日最高氣溫、14時(shí)空氣相對(duì)濕度和14時(shí)風(fēng)速組合，結(jié)合20 cm土壤相對(duì)濕度確定了3種類型小麥干熱風(fēng)指標(biāo)。不同地區(qū)干熱風(fēng)類型不同，評(píng)定標(biāo)準(zhǔn)也不同，同時(shí)干熱風(fēng)災(zāi)害主要發(fā)生在我國(guó)北方麥區(qū)，具有局地性特征。為了統(tǒng)一各麥區(qū)干熱風(fēng)指標(biāo)，本文采用了“三三三”標(biāo)準(zhǔn)對(duì)各麥區(qū)小麥生育期并集內(nèi)干熱風(fēng)有效天數(shù)進(jìn)行統(tǒng)計(jì)，也即當(dāng)小麥生育期并集內(nèi)某天同時(shí)滿足日最高氣溫>30℃、14時(shí)空氣相對(duì)濕度<30%、14時(shí)風(fēng)速>3 m/s三個(gè)條件時(shí)即為一個(gè)干熱風(fēng)有效天。干熱風(fēng)有效天數(shù)指小麥生育期并集內(nèi)干熱風(fēng)有效天的累積。

其次是普遍性特征,這表現(xiàn)在人的意識(shí)普遍被物化。在現(xiàn)代性社會(huì)中,工人的生產(chǎn)活動(dòng)是無(wú)批判的機(jī)械化直接性勞動(dòng),而精英階層也別無(wú)二致。無(wú)論是資產(chǎn)階級(jí)精英還是普通的工人大眾,都處在這個(gè)時(shí)間空間化的容器里,都沉浸在“意識(shí)的物化”之中,精英的工作也只是在維護(hù)這個(gè)合理性的機(jī)械過(guò)程,依然無(wú)法跳出物化的命運(yùn)。由此可見,人的命運(yùn)已深陷在物化意識(shí)的結(jié)構(gòu)之中。

2.2.3 洪澇與低溫冷害指標(biāo) 洪澇指標(biāo)的計(jì)算流程是：首先在像元尺度上計(jì)算小麥生育期內(nèi)某月份平均降水的歷年均值，然后用當(dāng)月平均降水實(shí)際值與歷年均值作差得到當(dāng)月平均降水距平值，最后取縣域內(nèi)耕地像元距平值的均值作為描述該縣當(dāng)月洪澇情況的指標(biāo)（下文簡(jiǎn)稱為月平均降水距平）。低溫冷害指標(biāo)的計(jì)算方法與洪澇指標(biāo)類似，計(jì)算過(guò)程不再贅述，下文簡(jiǎn)稱為月平均溫度距平。

2.2.4 遙感植被指數(shù)計(jì)算 首先對(duì)MOD13Q1-NDVI數(shù)據(jù)利用均值濾波從250 m重采樣至1 km，然后用全國(guó)土地利用分布圖對(duì)2000—2018年小麥生育期內(nèi)MOD13Q1-NDVI數(shù)據(jù)進(jìn)行耕地掩膜，最后用全國(guó)縣級(jí)矢量對(duì)各縣耕地范圍內(nèi)NDVI取均值，作為各縣的遙感植被指數(shù)。

2.3 趨勢(shì)產(chǎn)量擬合

在氣象估產(chǎn)模型中，通常將小麥實(shí)際產(chǎn)量分解為趨勢(shì)產(chǎn)量、氣象產(chǎn)量和隨機(jī)產(chǎn)量三部分，其中隨機(jī)產(chǎn)量往往忽略不計(jì)。因此，實(shí)際單產(chǎn)可以分解為趨勢(shì)單產(chǎn)和氣象單產(chǎn)之和。

式中：yt為小麥實(shí)際單產(chǎn)，Yt為趨勢(shì)單產(chǎn)，Yc為氣象單產(chǎn)，單位為kg/hm2。

本文利用直線滑動(dòng)平均法（包括3 a和5 a直線滑動(dòng)平均法）和HP濾波法[30]對(duì)五大麥區(qū)所有縣小麥實(shí)際單產(chǎn)進(jìn)行趨勢(shì)產(chǎn)量擬合，用實(shí)際單產(chǎn)減去趨勢(shì)單產(chǎn)得到氣象單產(chǎn)。以氣象單產(chǎn)為輸出變量，以四類氣象災(zāi)害指標(biāo)為輸入變量，分別在五大麥區(qū)構(gòu)建隨機(jī)森林回歸模型。模型擬合精度越高說(shuō)明對(duì)氣象產(chǎn)量擬合效果越好，也即對(duì)趨勢(shì)產(chǎn)量擬合效果越好，以此來(lái)篩選最適合的趨勢(shì)產(chǎn)量擬合方法。模型擬合精度用驗(yàn)證樣本的決定系數(shù)R2，均方根誤差（RMSE）和平均絕對(duì)誤差（MAE）進(jìn)行評(píng)價(jià)。

2.4 隨機(jī)森林回歸與驗(yàn)證

隨機(jī)森林算法的主要思想是從原始訓(xùn)練集中有放回地抽取k個(gè)樣本，且每個(gè)樣本的樣本容量均與原始訓(xùn)練集的大小一致；然后對(duì)每個(gè)樣本分別進(jìn)行回歸樹建模，得到k個(gè)建模結(jié)果，最后以每一棵回歸樹預(yù)測(cè)結(jié)果的平均值作為最終預(yù)測(cè)結(jié)果[31]。隨機(jī)森林估產(chǎn)模型的輸入變量包括趨勢(shì)單產(chǎn)1個(gè)、干熱風(fēng)有效天數(shù)1個(gè)、SPEI干旱指數(shù)1個(gè)、生育期內(nèi)月平均降水距平N個(gè)、生育期內(nèi)月平均溫度距平N個(gè)、生育期內(nèi)每16 d的NDVI植被指數(shù)M個(gè)，合計(jì)2N+M+3個(gè)。需要注意的是，不同麥區(qū)小麥的生育期不同，因此生育期內(nèi)月平均溫度距平和月平均降水距平的個(gè)數(shù)N在不同麥區(qū)將發(fā)生改變。以樣本集中3/4樣本子集作為訓(xùn)練樣本建立回歸模型，1/4樣本子集作為驗(yàn)證樣本進(jìn)行估產(chǎn)模型擬合精度的驗(yàn)證。隨機(jī)森林回歸模型擬合精度用模型決定系數(shù)R2、驗(yàn)證樣本的均方根誤差、平均絕對(duì)誤差以及平均相對(duì)誤差進(jìn)行評(píng)價(jià)。

2.5 災(zāi)害年份估產(chǎn)精度驗(yàn)證

采用麥區(qū)內(nèi)主要?。ㄖ陛犑校┺r(nóng)作物受災(zāi)面積累計(jì)的大小來(lái)表示各麥區(qū)整體受干旱、洪澇、低溫冷害三類氣象災(zāi)害導(dǎo)致的農(nóng)作物受災(zāi)情況。另外，由于年鑒統(tǒng)計(jì)數(shù)據(jù)中不包含干熱風(fēng)災(zāi)害統(tǒng)計(jì)數(shù)據(jù)，因此對(duì)于不同麥區(qū)干熱風(fēng)災(zāi)害年份的確定，主要參考了其他專家學(xué)者的相關(guān)研究。災(zāi)害年份確定之后，將五大麥區(qū)災(zāi)害年份對(duì)應(yīng)的樣本進(jìn)行剔除，對(duì)各麥區(qū)重新建立隨機(jī)森林回歸模型，然后將各麥區(qū)旱災(zāi)、洪澇災(zāi)害、低溫冷害、干熱風(fēng)災(zāi)害年份對(duì)應(yīng)的樣本的自變量代入重新訓(xùn)練好的模型進(jìn)行災(zāi)害年份產(chǎn)量預(yù)測(cè)。最后在縣級(jí)尺度和麥區(qū)尺度分別進(jìn)行了精度驗(yàn)證?？h級(jí)尺度驗(yàn)證用各麥區(qū)中所有縣小麥單產(chǎn)真實(shí)值和預(yù)測(cè)值的相對(duì)誤差的平均值表示。麥區(qū)尺度驗(yàn)證用各麥區(qū)實(shí)際單產(chǎn)和預(yù)測(cè)單產(chǎn)的相對(duì)誤差進(jìn)行評(píng)價(jià)。其中，各麥區(qū)實(shí)際單產(chǎn)是先用各縣實(shí)際單產(chǎn)與各縣實(shí)際種植面積計(jì)算得到麥區(qū)實(shí)際總產(chǎn)量，再用實(shí)際總產(chǎn)量除以總種植面積得到。同理，各麥區(qū)預(yù)測(cè)單產(chǎn)是先用各縣預(yù)測(cè)單產(chǎn)與各縣實(shí)際種植面積計(jì)算得到麥區(qū)總預(yù)測(cè)產(chǎn)量，再用麥區(qū)總預(yù)測(cè)產(chǎn)量除以總種植面積得到。

2.6 特征重要性評(píng)價(jià)方法

用隨機(jī)森林中某個(gè)樹節(jié)點(diǎn)的樣本方差與分裂后兩個(gè)葉子節(jié)點(diǎn)樣本方差進(jìn)行差值運(yùn)算，用該差值來(lái)評(píng)價(jià)某個(gè)特征的重要性程度，差值降低越高表示該特征重要性越大。計(jì)算得到模型所有特征重要性后，對(duì)所有特征重要性值作歸一化處理[32]。主要說(shuō)明的是：對(duì)于某一模型，所有特征重要性值累積為1，特征重要性具體的值本身并無(wú)意義，因此分析時(shí)關(guān)注的是特征重要性值之間的相對(duì)差異。

3 研究結(jié)果

3.1 趨勢(shì)產(chǎn)量擬合方法

分別采用3 a直線滑動(dòng)平均法、5 a直線滑動(dòng)平均法和HP濾波法對(duì)五大麥區(qū)所有縣小麥實(shí)際單產(chǎn)進(jìn)行趨勢(shì)產(chǎn)量擬合，進(jìn)而計(jì)算氣象單產(chǎn)，建立氣象單產(chǎn)和四類氣象災(zāi)害指標(biāo)之間的隨機(jī)森林回歸模型，對(duì)其進(jìn)行驗(yàn)證，結(jié)果如圖3所示。北部冬麥區(qū)和長(zhǎng)江中下游冬麥區(qū)中HP方法R2最高，但3 a直線滑動(dòng)平均法的均方根誤差和平均絕對(duì)誤差最低；黃淮海冬麥區(qū)和西北春麥區(qū)中3 a直線滑動(dòng)平均法R2最高，均方根誤差、平均絕對(duì)誤差最低；西南冬麥區(qū)5 a直線滑動(dòng)平均法R2最高，但3 a方法的均方根誤差、平均絕對(duì)誤差最低。綜合來(lái)看，3 a直線滑動(dòng)平均法擬合趨勢(shì)產(chǎn)量最好。因此，本文采用3 a直線滑動(dòng)平均法對(duì)趨勢(shì)產(chǎn)量進(jìn)行擬合，并進(jìn)行下一步估產(chǎn)模型建立。

圖3 氣象產(chǎn)量與氣象因子擬合精度Fig. 3 Fitting accuracy of meteorological yield and meteorological factors

3.2 隨機(jī)森林估產(chǎn)模型構(gòu)建與精度驗(yàn)證結(jié)果

表1為各麥區(qū)隨機(jī)森林估產(chǎn)模型擬合精度統(tǒng)計(jì)。各麥區(qū)估產(chǎn)模型擬合R2均達(dá)到了0.98以上，平均相對(duì)誤差均低于0.073。其中，黃淮海冬麥區(qū)估產(chǎn)模型擬合精度最高，決定系數(shù)R2為0.992、均方根誤差為125.183 kg/hm2，平均絕對(duì)誤差為76.703 kg/hm2、平均相對(duì)誤差為0.019；西北春麥區(qū)估產(chǎn)模型擬合精度相對(duì)最低，決定系數(shù)R2為0.990、均方根誤差為200.835 kg/hm2、平均絕對(duì)誤差為131.655 kg/hm2、平均相對(duì)誤差為0.073。綜上所述，用本文提出的方法對(duì)各麥區(qū)建立的隨機(jī)森林估產(chǎn)模型擬合精度整體很高，能夠滿足大范圍尺度估產(chǎn)精度。

表1 各麥區(qū)估產(chǎn)模型擬合精度Table 1 Fitting accuracy of yield estimation models in each wheat growing region

3.3 災(zāi)害年份估產(chǎn)精度驗(yàn)證

3.3.1 災(zāi)害年份確定 圖4至圖8為各麥區(qū)農(nóng)作物災(zāi)損統(tǒng)計(jì)情況，以2000—2019年《中國(guó)農(nóng)村統(tǒng)計(jì)年鑒》全國(guó)各省農(nóng)作物自然氣候?yàn)?zāi)害損失數(shù)據(jù)為基礎(chǔ)，以各麥區(qū)主要?。ㄖ陛犑校┺r(nóng)作物受災(zāi)面積累計(jì)的大小為依據(jù)，對(duì)五大麥區(qū)不同類型災(zāi)害的典型年份進(jìn)行確定。例如，西北春麥區(qū)主要包括甘肅、寧夏以及內(nèi)蒙古西北部地區(qū)，內(nèi)蒙古西北部地區(qū)多為戈壁沙漠地帶，故用寧夏與甘肅農(nóng)作物受災(zāi)面積的累計(jì)來(lái)表示該區(qū)整體農(nóng)作物受災(zāi)狀況；西南冬麥區(qū)主要包括重慶與貴州全部、四川中東部、云南中北部以及陜西南部少部分地區(qū)，故用重慶與貴州農(nóng)作物受災(zāi)面積累計(jì)并參考四川、云南農(nóng)作物受災(zāi)面積來(lái)表示該區(qū)整體農(nóng)作物受災(zāi)狀況，其他三個(gè)麥區(qū)進(jìn)行類似分析來(lái)確定不同類型災(zāi)害的典型年份。

圖4 西北春麥區(qū)農(nóng)作物受災(zāi)面積統(tǒng)計(jì)Fig. 4 Area of crop affected by disaster in northwest spring wheat area

圖5 北部冬麥區(qū)農(nóng)作物受災(zāi)面積統(tǒng)計(jì)Fig. 5 Area of crop affected by disaster in northern winter wheat area

圖8 西南冬麥區(qū)農(nóng)作物受災(zāi)面積統(tǒng)計(jì)Fig. 8 Area of crop affected by disaster in southwest winter wheat area

由圖可以看出，低溫冷害、干旱和洪澇災(zāi)害嚴(yán)重年份在西北春麥區(qū)分別為2004年、2000年和2013年，在北部冬麥區(qū)分別為2006年、2001年和2007年，在黃淮海冬麥區(qū)分別為2005年、2001年和2003年，在長(zhǎng)江中下游冬麥區(qū)分別為2008年、2001年和2003年，在西南冬麥區(qū)分別為2008年、2001年和2002年。干熱風(fēng)災(zāi)害有較強(qiáng)的區(qū)域性，其主要集中在我國(guó)河西地區(qū)、黃淮海平原等地區(qū)[33]，即西北春麥區(qū)、北方冬麥區(qū)以及黃淮海冬麥區(qū)。已有研究表明[26,29,34-35]，1960年至2017年河西地區(qū)干熱風(fēng)發(fā)生總次數(shù)和日數(shù)先緩慢減少后迅速增多，其中2004年的干熱風(fēng)發(fā)生總天數(shù)最高，持續(xù)天數(shù)也達(dá)到了歷年最長(zhǎng)（12 d）；1961年至2015年黃淮海平原地區(qū)干熱風(fēng)日數(shù)總體呈下降趨勢(shì)，但是在2001年黃淮海平原地區(qū)仍發(fā)生了一次嚴(yán)重的干熱風(fēng)災(zāi)害，區(qū)域平均干熱風(fēng)日數(shù)達(dá)8.1 d。因此，本文將 2004年作為西北春麥區(qū)干熱風(fēng)災(zāi)害年份，2001年作為北部冬麥區(qū)與黃淮海冬麥區(qū)的干熱風(fēng)災(zāi)害年份。長(zhǎng)江中下游冬麥區(qū)與西南冬麥區(qū)未進(jìn)行干熱風(fēng)災(zāi)害年份的確定。

圖6 黃淮海冬麥區(qū)農(nóng)作物受災(zāi)面積統(tǒng)計(jì)Fig. 6 Area of crop affected by disaster in Huang-Huai-Hai winter wheat area

圖7 長(zhǎng)江中下游冬麥區(qū)農(nóng)作物受災(zāi)面積統(tǒng)計(jì)Fig. 7 Area of crop affected by disaster in middle and lower reaches of Yangtze River winter wheat area

3.3.2 災(zāi)害年份估產(chǎn)精度驗(yàn)證結(jié)果 表2為縣級(jí)和麥區(qū)級(jí)兩個(gè)尺度的驗(yàn)證結(jié)果。各麥區(qū)所有縣單產(chǎn)平均相對(duì)誤差的最大值為西北春麥區(qū)洪澇災(zāi)害年份的估產(chǎn)精度（0.060），最小值為黃淮海冬麥區(qū)洪澇災(zāi)害年份的估產(chǎn)精度（0.017），平均相對(duì)誤差的均值為0.032。黃淮海冬麥區(qū)災(zāi)害年份的估產(chǎn)精度高于其他麥區(qū)災(zāi)害年份的估產(chǎn)精度，西北春麥區(qū)災(zāi)害年份估產(chǎn)精度相對(duì)最低。在不同類型災(zāi)害年份，五大麥區(qū)單產(chǎn)預(yù)測(cè)精度整體較高，整體相對(duì)誤差最大為0.049，最小為0.001，均值為0.018。

表2 精度驗(yàn)證結(jié)果Table 2 Accuracy verification results

3.4 輸入特征重要性排序結(jié)果

不同小麥種植區(qū)有各自的生態(tài)環(huán)境特點(diǎn)，小麥在不同麥區(qū)受到的異常氣候?yàn)?zāi)害類型差異較大。由圖9可以看出，西北春麥區(qū)和北部冬麥區(qū)變量重要性排序中，干熱風(fēng)有效天數(shù)重要性在所有變量中排名分別達(dá)到了第3和第13位，在所有氣象類因子重要性排名中則分別達(dá)到了第1位和第6位，而在其他麥區(qū)，干熱風(fēng)災(zāi)害因子重要性都很低，這表明干熱風(fēng)災(zāi)害是影響西北春麥區(qū)和北部冬麥區(qū)小麥生長(zhǎng)的重要因素。

圖9 各麥區(qū)輸入變量重要性排序Fig. 9 Order of importance of input variables in five wheat growing regions

各麥區(qū)SPEI干旱指數(shù)在所有氣象類因子重要性排序中均達(dá)到較高排名，分別為西北春麥區(qū)第2位、北部冬麥區(qū)第2位、黃淮海冬麥區(qū)第1位、長(zhǎng)江中下游冬麥區(qū)第1位、西南冬麥區(qū)第3位，表明干旱仍是五大麥區(qū)影響小麥生產(chǎn)的重要因素。

不同生育階段溫度變化對(duì)旱地春小麥產(chǎn)量的影響大小依次為灌漿成熟期、開花灌漿期、播種出苗期、孕穗開花期、拔節(jié)孕穗期、出苗分蘗期、分蘗拔節(jié)期[36]，需水關(guān)鍵期為出苗期、拔節(jié)期、開花期[37]，本文西北春麥區(qū)不同月份降水與溫度參數(shù)重要性計(jì)算統(tǒng)計(jì)結(jié)果與其大致相符。冬小麥在拔節(jié)期對(duì)水分脅迫敏感性最強(qiáng)，其次是抽穗灌漿期、返青期、灌漿成熟期[38]，除西南冬麥區(qū)外，其他冬麥區(qū)降水參數(shù)重要性統(tǒng)計(jì)結(jié)果與其大致相符，不同月份的溫度參數(shù)重要性分布較為均勻，整體差異不大。

各麥區(qū)不同時(shí)間NDVI植被指數(shù)重要性差異較大，每個(gè)麥區(qū)都有一個(gè)時(shí)間段的NDVI植被指數(shù)重要性相對(duì)較高。具體來(lái)說(shuō)，西北春麥區(qū)為16號(hào)變量（3月22日至4月6日）、北部冬麥區(qū)為16號(hào)變量（4月7日至4月22日）、黃淮海冬麥區(qū)為16號(hào)變量（4月7日至4月22日）、西南冬麥區(qū)為4號(hào)變量（9月14日至9月29日）?？傮w來(lái)說(shuō)，北部冬麥區(qū)、黃淮海冬麥區(qū)重要性較高的NDVI植被指數(shù)的時(shí)間大致處于小麥生長(zhǎng)的拔節(jié)期與抽穗成熟期，而在西北春麥區(qū)與西南冬麥區(qū)，處于小麥種植前期的NDVI植被指數(shù)顯示出較高的重要性。其中值得注意的一點(diǎn)是長(zhǎng)江中下游冬麥區(qū)的18號(hào)變量（5月25日至6月9日的NDVI）顯示重要性很高，其次是15號(hào)變量（4月7日至4月22日），前者小麥已經(jīng)進(jìn)入收獲期，后者在抽穗期，理論上后者才是更重要，出現(xiàn)這個(gè)問題的原因詳見后文的討論。

4 討論

4.1 本研究的主要科學(xué)貢獻(xiàn)

1）傳統(tǒng)的氣象估產(chǎn)模型的輸入變量往往是原始的氣象要素，如降水量、溫度等，而非氣候指標(biāo)。氣候指標(biāo)是一定氣候條件下的單項(xiàng)氣候要素或多種氣候要素綜合的特征量，如干旱指數(shù)，熱浪指數(shù)，高溫指數(shù)等。氣候指標(biāo)能更加綜合的反映氣候條件的變化，特別是考慮作物生長(zhǎng)條件的農(nóng)業(yè)氣候指標(biāo)，能更好的衡量農(nóng)業(yè)氣象條件的利弊。本文從影響小麥生長(zhǎng)的主要農(nóng)業(yè)自然災(zāi)害（干旱、洪澇、干熱風(fēng)和低溫凍害）出發(fā)，計(jì)算了SPEI干旱指數(shù)、降水距平、干熱風(fēng)指數(shù)和溫度距平來(lái)反映研究區(qū)中四個(gè)主要農(nóng)業(yè)自然災(zāi)害的情況，據(jù)此構(gòu)建的估產(chǎn)模型經(jīng)過(guò)驗(yàn)證可以有效的進(jìn)行災(zāi)害年份小麥的產(chǎn)量估計(jì)。

2）NDVI指數(shù)（包括NDVI的最大值、生長(zhǎng)季的NDVI累計(jì)值等）是遙感估產(chǎn)中最常用的植被指數(shù)。傳統(tǒng)的遙感統(tǒng)計(jì)估產(chǎn)模型往往只依托于遙感指數(shù)，而不考慮非遙感數(shù)據(jù)。本研究綜合使用遙感數(shù)據(jù)和氣象數(shù)據(jù)進(jìn)行估產(chǎn)模型的建立。相比單純的遙感數(shù)據(jù)，加入氣候數(shù)據(jù)能夠增加模型估產(chǎn)結(jié)果的穩(wěn)定性，減少由傳感器衰變，大氣條件差異等造成的植被指數(shù)本身年際間的變化而產(chǎn)生的估產(chǎn)誤差。相比單純使用氣候指數(shù)，遙感數(shù)據(jù)能夠在一定程度上反映人類活動(dòng)的影響，例如人工灌溉對(duì)干旱的抵御作用，從而避免僅僅考慮氣候?yàn)?zāi)害指數(shù)而帶來(lái)的產(chǎn)量低估。

3）本文利用隨機(jī)森林算法對(duì)各個(gè)麥區(qū)估產(chǎn)模型的輸入變量的重要性進(jìn)行了分析，可以幫助我們更好的理解影響各個(gè)麥區(qū)小麥產(chǎn)量的主要因素，為不同區(qū)域的小麥防災(zāi)減災(zāi)措施的制定提供參考，也可以為其他區(qū)域的小麥估產(chǎn)研究和其他作物估產(chǎn)模型的建立提供輸入變量篩選的方法借鑒。

4.2 本研究的不足之處

1）計(jì)算2000—2018年小麥種植區(qū)各縣四類災(zāi)害指標(biāo)與遙感植被指數(shù)時(shí)，考慮到工作量和試驗(yàn)可行性，以及每年對(duì)全國(guó)小麥空間分布提取也同樣存在一定的誤差，最終選擇用2010年全國(guó)土地利用空間分類圖中耕地范圍來(lái)代替各縣小麥種植區(qū)進(jìn)行參數(shù)的計(jì)算，這會(huì)在一定程度上帶來(lái)誤差。例如對(duì)于NDVI變量，長(zhǎng)江中下游冬麥區(qū)的5月25日至6月9日的NDVI重要性顯示很高，其次是4月7日至4月22日的NDVI，前者小麥已經(jīng)進(jìn)入收獲期，后者在抽穗期，這一結(jié)果和理論預(yù)期不符合。究其原因可能是長(zhǎng)江中下游的小麥種植面積總體小于水稻的種植面積，小麥?zhǔn)斋@期時(shí)以耕地范圍統(tǒng)計(jì)出來(lái)的NDVI指數(shù)很大程度上是水稻的信息，5月25日至6月9日大致在早稻的拔節(jié)孕穗期和中稻的移栽期，此時(shí)發(fā)生干旱對(duì)水稻的產(chǎn)量影響很大。

2）本文根據(jù)各省（直轄市）農(nóng)作物自然災(zāi)害損失數(shù)據(jù)進(jìn)行各麥區(qū)災(zāi)害年份的確定，一方面災(zāi)損數(shù)據(jù)并不只針對(duì)小麥，另一方面災(zāi)損數(shù)據(jù)是按照行政單元統(tǒng)計(jì)的，而麥區(qū)是跨行政單元的。我們以麥區(qū)包括的主要省（直轄市）的災(zāi)損總量來(lái)反映整個(gè)麥區(qū)的災(zāi)損情況，由此確定各麥區(qū)不同災(zāi)害類型的典型災(zāi)害年份，可能會(huì)產(chǎn)生一定誤差。在后續(xù)的研究中可以考慮根據(jù)時(shí)間序列氣象數(shù)據(jù)進(jìn)行不同類型災(zāi)害年份的確定。

3）本文選擇的SPEI干旱指數(shù)和干熱風(fēng)有效天數(shù)是對(duì)小麥生育期整體災(zāi)害情況的描述，然而發(fā)生在不同生育期的相同災(zāi)害對(duì)作物產(chǎn)量的影響是不同的，在后續(xù)的研究中需要在小麥不同生育期對(duì)災(zāi)害指標(biāo)進(jìn)行計(jì)算，以便分析在小麥不同生長(zhǎng)階段發(fā)生的氣象災(zāi)害對(duì)產(chǎn)量的影響差異。

5 結(jié)論

本文基于隨機(jī)森林算法，以氣象災(zāi)害指標(biāo)、遙感植被指數(shù)、趨勢(shì)產(chǎn)量為輸入變量，構(gòu)建了五大麥區(qū)的估產(chǎn)模型，并在典型災(zāi)害年份進(jìn)行了估產(chǎn)精度的驗(yàn)證，得到如下主要結(jié)論：

1) 精度驗(yàn)證結(jié)果顯示，各麥區(qū)回歸模型的決定系數(shù)R2均達(dá)到了0.95以上，驗(yàn)證樣本的平均相對(duì)誤差均低于0.073，均方根誤差均低于200.835 kg/hm2，平均絕對(duì)誤差均低于131.655 kg/hm2。各麥區(qū)所有縣災(zāi)害年份實(shí)際單產(chǎn)與預(yù)測(cè)單產(chǎn)平均相對(duì)誤差均低于0.060，各麥區(qū)區(qū)級(jí)實(shí)際單產(chǎn)與預(yù)測(cè)單產(chǎn)的相對(duì)誤差均低于0.049。

2）輸入特征的重要性在不同麥區(qū)具有一定的差異性，與地區(qū)氣候差異情況有關(guān)。SPEI干旱指數(shù)在五大麥區(qū)重要性均較高，表明干旱仍然是五大麥區(qū)影響小麥生產(chǎn)的重要因素；干熱風(fēng)指標(biāo)重要性在西北春麥區(qū)與北部春麥區(qū)較高，符合干熱風(fēng)主要發(fā)生在河西與黃淮海地區(qū)的局地性特點(diǎn)；各麥區(qū)月平均溫度距平重要性與月平均降水距平重要性整體差異比較小；總體來(lái)說(shuō)拔節(jié)期與抽穗期的小麥NDVI相比其他時(shí)間段的NDVI更重要。