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

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

從20世紀(jì)上半葉至今，地球大氣中的CO2濃度不斷升高，變暖趨勢(shì)越來越顯著，全球氣候變化，導(dǎo)致異常氣候事件頻發(fā)[1]。相關(guān)研究表明，盡管農(nóng)業(yè)生產(chǎn)受氣候變化的影響在不同區(qū)域有一定的差異，但整體上弊大于利[2]。我國農(nóng)業(yè)受氣候變化的影響十分嚴(yán)重[3]，在各類農(nóng)業(yè)氣象災(zāi)害中，干旱和洪澇是影響我國農(nóng)作物生產(chǎn)最主要的氣象災(zāi)害。1980年以來，我國洪旱災(zāi)害頻發(fā)。1998年的長江和松花江流域特大洪水事件和2016年的長江中下游及太湖洪澇事件直接帶來上千億元經(jīng)濟(jì)損失[4]。相比洪澇，干旱對(duì)我國農(nóng)業(yè)生產(chǎn)的影響更為嚴(yán)重[5]。有關(guān)統(tǒng)計(jì)數(shù)據(jù)表明，2007至2017年全國農(nóng)作物旱災(zāi)面積年均值為1.53×107hm2，糧食損失年均值為2.14億t，造成直接經(jīng)濟(jì)損失年均值為282.15億元，占GDP均值的0.21%[6]。在全國各種氣象災(zāi)害中，干旱導(dǎo)致的農(nóng)作物受災(zāi)面積占比約56%，洪澇導(dǎo)致的農(nóng)作物受災(zāi)面積占比約24%[7]。此外，低溫冷害、干熱風(fēng)等氣候事件的頻繁發(fā)生也給我國農(nóng)作物安全生產(chǎn)帶來了巨大挑戰(zhàn)。在中國北方麥區(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ì)模型、光能利用率模型、作物生長模擬模型和耦合模型[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ù)測。第三種，將作物產(chǎn)量分解為趨勢(shì)產(chǎn)量和氣象產(chǎn)量兩部分，前者反映由技術(shù)進(jìn)步（如灌溉、施肥、新品種等）導(dǎo)致的產(chǎn)量的長期變化，后者反映由自然氣候要素（光照、降水、輻射等）引起的產(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ù)測精度高，不易產(chǎn)生過擬合等優(yōu)點(diǎn)，不需要顧慮多重共線性的問題，與神經(jīng)網(wǎng)絡(luò)和線性回歸方法相比，其表現(xiàn)更穩(wěn)定，具有較強(qiáng)的抗干擾能力。該算法已經(jīng)廣泛應(yīng)用于干旱監(jiān)測[18]、山體滑坡預(yù)測[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ù)，在中國五大小麥種植區(qū)分別建立估產(chǎn)模型，并驗(yàn)證其在典型災(zāi)害年份的估產(chǎn)精度，為探究氣候變化對(duì)小麥產(chǎn)量的影響，以及建立滿足異常氣候下估產(chǎn)精度的模型提供參考。

1 研究區(qū)概況

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

圖1 中國小麥種植區(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)量為輸出變量，建立樣本庫；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ù)處理

通過查閱資料發(fā)現(xiàn)在我國眾多不同類型氣象災(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í)期長勢(shì)情況的遙感植被指數(shù)指標(biāo)，用于估產(chǎn)模型的建立。具體包括：1）來自全球標(biāo)準(zhǔn)化降水蒸散指數(shù)（Standardized Precipitation Evapotranspiration Index, SPEI）數(shù)據(jù)庫1999—2018年的SPEIbase v2.6數(shù)據(jù)集，時(shí)間尺度為1～48月，分辨率為0.5°；2）來自歐空局哥白尼氣象數(shù)據(jù)中心的1999—2018年的月平均降水量、月平均溫度、14時(shí)風(fēng)速、14時(shí)相對(duì)濕度、日最高溫度，除了日最高溫度分辨率為0.1°外，其他數(shù)據(jù)的空間分辨率均為0.25°；3）來自美國航空航天局的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月、長江中下游冬麥區(qū)10月-次年5月、西北春麥區(qū)3月-8月，因此本文以各麥區(qū)小麥的生育期并集為研究時(shí)段來對(duì)參數(shù)進(jìn)行計(jì)算。

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

2.2.2 干熱風(fēng)指標(biāo) 我國小麥干熱風(fēng)具體可以分為高溫低濕型、雨后青枯型和旱風(fēng)型3種[29]。國家氣象局發(fā)布了2019年干熱風(fēng)國家標(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ā)生在我國北方麥區(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)是無批判的機(jī)械化直接性勞動(dòng),而精英階層也別無二致。無論是資產(chǎn)階級(jí)精英還是普通的工人大眾,都處在這個(gè)時(shí)間空間化的容器里,都沉浸在“意識(shí)的物化”之中,精英的工作也只是在維護(hù)這個(gè)合理性的機(jī)械過程,依然無法跳出物化的命運(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)（下文簡稱為月平均降水距平）。低溫冷害指標(biāo)的計(jì)算方法與洪澇指標(biāo)類似，計(jì)算過程不再贅述，下文簡稱為月平均溫度距平。

2.2.4 遙感植被指數(shù)計(jì)算 首先對(duì)MOD13Q1-NDVI數(shù)據(jù)利用均值濾波從250 m重采樣至1 km，然后用全國土地利用分布圖對(duì)2000—2018年小麥生育期內(nèi)MOD13Q1-NDVI數(shù)據(jù)進(jìn)行耕地掩膜，最后用全國縣級(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ī)森林回歸模型。模型擬合精度越高說明對(duì)氣象產(chǎn)量擬合效果越好，也即對(duì)趨勢(shì)產(chǎn)量擬合效果越好，以此來篩選最適合的趨勢(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ù)測結(jié)果的平均值作為最終預(yù)測結(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ì)的大小來表示各麥區(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ù)測。最后在縣級(jí)尺度和麥區(qū)尺度分別進(jìn)行了精度驗(yàn)證。縣級(jí)尺度驗(yàn)證用各麥區(qū)中所有縣小麥單產(chǎn)真實(shí)值和預(yù)測值的相對(duì)誤差的平均值表示。麥區(qū)尺度驗(yàn)證用各麥區(qū)實(shí)際單產(chǎn)和預(yù)測單產(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ù)測單產(chǎn)是先用各縣預(yù)測單產(chǎn)與各縣實(shí)際種植面積計(jì)算得到麥區(qū)總預(yù)測產(chǎn)量，再用麥區(qū)總預(yù)測產(chǎn)量除以總種植面積得到。

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

用隨機(jī)森林中某個(gè)樹節(jié)點(diǎn)的樣本方差與分裂后兩個(gè)葉子節(jié)點(diǎn)樣本方差進(jìn)行差值運(yùn)算，用該差值來評(píng)價(jià)某個(gè)特征的重要性程度，差值降低越高表示該特征重要性越大。計(jì)算得到模型所有特征重要性后，對(duì)所有特征重要性值作歸一化處理[32]。主要說明的是：對(duì)于某一模型，所有特征重要性值累積為1，特征重要性具體的值本身并無意義，因此分析時(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ū)和長江中下游冬麥區(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ì)誤差最低。綜合來看，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年《中國農(nóng)村統(tǒng)計(jì)年鑒》全國各省農(nóng)作物自然氣候?yàn)?zāi)害損失數(shù)據(jù)為基礎(chǔ)，以各麥區(qū)主要省（直轄市）農(nóng)作物受災(zāi)面積累計(jì)的大小為依據(jù)，對(duì)五大麥區(qū)不同類型災(zāi)害的典型年份進(jìn)行確定。例如，西北春麥區(qū)主要包括甘肅、寧夏以及內(nèi)蒙古西北部地區(qū)，內(nèi)蒙古西北部地區(qū)多為戈壁沙漠地帶，故用寧夏與甘肅農(nóng)作物受災(zāi)面積的累計(jì)來表示該區(qū)整體農(nóng)作物受災(zāi)狀況；西南冬麥區(qū)主要包括重慶與貴州全部、四川中東部、云南中北部以及陜西南部少部分地區(qū)，故用重慶與貴州農(nóng)作物受災(zāi)面積累計(jì)并參考四川、云南農(nóng)作物受災(zāi)面積來表示該區(qū)整體農(nóng)作物受災(zāi)狀況，其他三個(gè)麥區(qū)進(jìn)行類似分析來確定不同類型災(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年，在長江中下游冬麥區(qū)分別為2008年、2001年和2003年，在西南冬麥區(qū)分別為2008年、2001年和2002年。干熱風(fēng)災(zāi)害有較強(qiáng)的區(qū)域性，其主要集中在我國河西地區(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á)到了歷年最長（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)害年份。長江中下游冬麥區(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 長江中下游冬麥區(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ù)測精度整體較高，整體相對(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ū)小麥生長的重要因素。

圖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位、長江中下游冬麥區(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ì)較高。具體來說，西北春麥區(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來說，北部冬麥區(qū)、黃淮海冬麥區(qū)重要性較高的NDVI植被指數(shù)的時(shí)間大致處于小麥生長的拔節(jié)期與抽穗成熟期，而在西北春麥區(qū)與西南冬麥區(qū)，處于小麥種植前期的NDVI植被指數(shù)顯示出較高的重要性。其中值得注意的一點(diǎn)是長江中下游冬麥區(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)能更加綜合的反映氣候條件的變化，特別是考慮作物生長條件的農(nóng)業(yè)氣候指標(biāo)，能更好的衡量農(nóng)業(yè)氣象條件的利弊。本文從影響小麥生長的主要農(nóng)業(yè)自然災(zāi)害（干旱、洪澇、干熱風(fēng)和低溫凍害）出發(fā)，計(jì)算了SPEI干旱指數(shù)、降水距平、干熱風(fēng)指數(shù)和溫度距平來反映研究區(qū)中四個(gè)主要農(nóng)業(yè)自然災(zāi)害的情況，據(jù)此構(gòu)建的估產(chǎn)模型經(jīng)過驗(yàn)證可以有效的進(jìn)行災(zāi)害年份小麥的產(chǎn)量估計(jì)。

2）NDVI指數(shù)（包括NDVI的最大值、生長季的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ù)而帶來的產(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ì)全國小麥空間分布提取也同樣存在一定的誤差，最終選擇用2010年全國土地利用空間分類圖中耕地范圍來代替各縣小麥種植區(qū)進(jìn)行參數(shù)的計(jì)算，這會(huì)在一定程度上帶來誤差。例如對(duì)于NDVI變量，長江中下游冬麥區(qū)的5月25日至6月9日的NDVI重要性顯示很高，其次是4月7日至4月22日的NDVI，前者小麥已經(jīng)進(jìn)入收獲期，后者在抽穗期，這一結(jié)果和理論預(yù)期不符合。究其原因可能是長江中下游的小麥種植面積總體小于水稻的種植面積，小麥?zhǔn)斋@期時(shí)以耕地范圍統(tǒng)計(jì)出來的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)損總量來反映整個(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ì)算，以便分析在小麥不同生長階段發(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ù)測單產(chǎn)平均相對(duì)誤差均低于0.060，各麥區(qū)區(qū)級(jí)實(shí)際單產(chǎn)與預(yù)測單產(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ū)月平均溫度距平重要性與月平均降水距平重要性整體差異比較??；總體來說拔節(jié)期與抽穗期的小麥NDVI相比其他時(shí)間段的NDVI更重要。