徐 寧 ,李發(fā)東 ,張秋英 ,艾治頻 ,冷佩芳 ,舒 旺 ,田 超,李 兆,陳 剛,喬云峰**

(1.中國(guó)科學(xué)院地理科學(xué)與資源研究所生態(tài)系統(tǒng)網(wǎng)絡(luò)觀測(cè)與模擬重點(diǎn)實(shí)驗(yàn)室 北京 100101;2.中國(guó)科學(xué)院大學(xué)中丹學(xué)院北京 100049;3.中國(guó)環(huán)境科學(xué)研究院 北京 100012;4.美國(guó)佛羅里達(dá)州立大學(xué) 塔拉哈西 32306)

盡管全球綠色革命取得了成功,糧食產(chǎn)量提高了100%～200%[1],但全球仍有超過(guò)10 億人面臨饑餓和營(yíng)養(yǎng)不良的困境[2]。在非洲,糧食缺乏尤為嚴(yán)重,21.5%的人口處于高度糧食不安全狀態(tài),高于全球平均水平的9.2%[3-5]。埃塞俄比亞是糧食不安全情況最嚴(yán)重的國(guó)家之一[6],因此農(nóng)業(yè)對(duì)于埃塞俄比亞人民的生活至關(guān)重要,在嚴(yán)酷的氣候條件下糧食產(chǎn)量的變化也嚴(yán)重影響著該國(guó)的糧食安全。

氣候變化影響下的全球糧食安全問題是21 世紀(jì)最重要的挑戰(zhàn)之一[7-11]。通過(guò)研究氣候變化對(duì)玉米(Zea mays)、小麥(Triticum aestivum)和水稻(Oryza sativa)等作物生產(chǎn)的影響[12-21],以及氣候變化對(duì)流域水資源[22-23]、森林[24]、工業(yè)和自然景觀的影響[25],學(xué)界發(fā)現(xiàn)氣候變化是影響作物年產(chǎn)量的最重要因素之一,氣候變化的不確定性將大大增加糧食生產(chǎn)的不確定性。

氣候變化情景是基于一系列氣候關(guān)系和輻射強(qiáng)迫假設(shè)對(duì)未來(lái)氣候的合理描述[26-29]。第六次國(guó)際耦合模式比較計(jì)劃(CMIP6)應(yīng)用了最新的共享社會(huì)經(jīng)濟(jì)路徑-代表性濃度路徑情景(SSP-RCP)模擬不同的未來(lái)發(fā)展情景,并通過(guò)全球氣候模式(GCM)[28]進(jìn)行氣候變化預(yù)測(cè)。GCM 是模擬當(dāng)前和未來(lái)氣候變化的重要工具,可以使研究人員更 好地了解氣候變化在區(qū)域尺度上對(duì)作物生產(chǎn)的影響[30]。

近年來(lái),隨著人工智能的快速發(fā)展,機(jī)器學(xué)習(xí)逐漸被引入各個(gè)應(yīng)用領(lǐng)域的研究中。機(jī)器學(xué)習(xí)應(yīng)用于產(chǎn)量預(yù)測(cè)的主要優(yōu)勢(shì)是能夠基于復(fù)雜的非線性農(nóng)業(yè)數(shù)據(jù)建立作物產(chǎn)量預(yù)測(cè)模型[31]。Gonzalez-Sanchez等[32]選擇了10 種作物,利用5 種機(jī)器學(xué)習(xí)模型構(gòu)建了作物產(chǎn)量預(yù)測(cè)模型并評(píng)估了每種機(jī)器學(xué)習(xí)算法的表現(xiàn)。Filippi 等[33]結(jié)合多屬性數(shù)據(jù)劃分生長(zhǎng)期,并構(gòu)建隨機(jī)森林產(chǎn)量預(yù)測(cè)模型。Khanal 等[34]采集了土壤樣本,結(jié)合多光譜影像數(shù)據(jù),構(gòu)建了玉米產(chǎn)量預(yù)測(cè)模型。綜合來(lái)看,機(jī)器學(xué)習(xí)在作物產(chǎn)量預(yù)測(cè)方面的精度優(yōu)于一般統(tǒng)計(jì)方法,通過(guò)選取適當(dāng)算法、合理添加關(guān)鍵參數(shù)、精準(zhǔn)調(diào)參等優(yōu)化后的機(jī)器學(xué)習(xí)可勝任數(shù)據(jù)稀缺地區(qū)(如埃塞俄比亞)糧食產(chǎn)量預(yù)測(cè)。

埃塞俄比亞約有9100 萬(wàn)人口,位居非洲第二[35],其主要糧食作物有玉米、苔麩(Eragrostis tef)、大麥(Hordeum vulgare)、小麥、高粱(Sorghum bicolor)及各種豆類[35]。該國(guó)農(nóng)業(yè)產(chǎn)業(yè)規(guī)模較小、作物生產(chǎn)力不高,同時(shí)糧食生產(chǎn)結(jié)構(gòu)較為單一,這就直接導(dǎo)致了該地區(qū)農(nóng)業(yè)敏感性較高,抵御氣候波動(dòng)的能力較差[36]。與此同時(shí),嚴(yán)重的水土流失和干旱也對(duì)其農(nóng)業(yè)生產(chǎn)造成不良影響[37-38]。自“一帶一路”倡議中明確把非洲納入“21 世紀(jì)海上絲綢之路”以來(lái),東非國(guó)家在“一帶一路”合作中所處的地位愈發(fā)重要[19],其中埃塞俄比亞的發(fā)展對(duì)“一帶一路”倡議在非洲的推進(jìn)發(fā)揮著重要作用。糧食安全作為保障埃塞俄比亞健康發(fā)展的重點(diǎn)抓手,是我國(guó)國(guó)家戰(zhàn)略推進(jìn)的關(guān)鍵制衡因素,也是實(shí)現(xiàn)聯(lián)合國(guó)可持續(xù)發(fā)展目標(biāo)的國(guó)際研究熱點(diǎn)。

鑒于此,本研究擬以埃塞俄比亞為研究區(qū),基于歷年作物產(chǎn)量數(shù)據(jù)、未來(lái)氣候變化模式模擬數(shù)據(jù)和其他農(nóng)業(yè)相關(guān)的環(huán)境及社會(huì)經(jīng)濟(jì)數(shù)據(jù),針對(duì)埃塞俄比亞產(chǎn)量排名前五的糧食作物(即苔麩、玉米、高粱、小麥和大麥),利用多種機(jī)器學(xué)習(xí)算法訓(xùn)練產(chǎn)量預(yù)測(cè)模型并進(jìn)行預(yù)測(cè),通過(guò)多種模型結(jié)果的分析對(duì)比,確定表現(xiàn)最佳的機(jī)器學(xué)習(xí)算法,進(jìn)一步結(jié)合CMIP6 中提供的SSP-RCP 情景,預(yù)測(cè)研究區(qū)2050 年之前不同未來(lái)情景下糧食產(chǎn)量的變化,并對(duì)研究區(qū)未來(lái)糧食生產(chǎn)情況進(jìn)行分析。

1 材料與方法

1.1 研究區(qū)概況

埃塞俄 比亞位 于非洲 東北部(3°～15°N,33°～48°E),國(guó)土面積112.68 萬(wàn)km2,北部、南部、東北部的沙漠和半沙漠約占領(lǐng)土面積的28%。65%的國(guó)土為可耕地,實(shí)耕13.2 萬(wàn)km2,其中糧食耕地占3/4[36-37]。埃塞俄比亞地區(qū)海拔主要分布在2500～3000 m。國(guó)家劃分3 級(jí)行政區(qū)級(jí)別,分別是州(Region)、區(qū)(Zone)、縣(Woreda),行政區(qū)劃如圖1 所示。埃塞俄比亞分布有7 種氣候類型,兼有熱帶氣候、溫帶氣候以及干旱氣候。高海拔地區(qū)年平均氣溫約為15℃,低海拔地區(qū)年平均氣溫約為25 ℃。此外,埃塞俄比亞平均年降雨量在西部出現(xiàn)極大值(＞2000 mm·a-1),在東北低海拔地區(qū)以及東南部出現(xiàn)極小值(＜200 mm·a-1)[5-6]。

圖1 埃塞俄比亞行政區(qū)劃Fig.1 Administrative division of Ethiopia

1.2 數(shù)據(jù)準(zhǔn)備

1.2.1 糧食產(chǎn)量數(shù)據(jù)

通過(guò)檢索聯(lián)合國(guó)糧食及農(nóng)業(yè)組織統(tǒng)計(jì)數(shù)據(jù)庫(kù)(https://www.fao.org/faostat/en/#home)、埃塞俄比亞中央統(tǒng)計(jì)局?jǐn)?shù)據(jù)庫(kù)(http://www.csa.gov.et/),獲取1995—2020 年埃塞俄比亞市級(jí)5 種主要糧食作物(即苔麩、玉米、小麥、大麥和高粱)單位面積產(chǎn)量數(shù)據(jù),其原始數(shù)據(jù)的產(chǎn)量單位均為t·hm-2。

埃塞俄比亞全年糧食生產(chǎn)活動(dòng)主要?jiǎng)澐譃閮蓚€(gè)主要生產(chǎn)季: 梅赫季(Meher season)和貝爾格季(Belg season)。梅赫季(每年的4—12 月)糧食產(chǎn)量約占全年糧食產(chǎn)量的90%以上,其糧食的主要生長(zhǎng)時(shí)間為7—9月;貝爾格季為每年的2—9 月,糧食的主要生長(zhǎng)時(shí)間為4—5 月。研究中所使用的作物產(chǎn)量數(shù)據(jù)包含該作物當(dāng)年梅赫季產(chǎn)量和貝爾格季產(chǎn)量。由于貝爾格季作物產(chǎn)量較低,且部分地區(qū)不適宜糧食作物種植,該季數(shù)據(jù)在時(shí)間和空間尺度上的缺失均較為嚴(yán)重。

1.2.2 GCM 數(shù)據(jù)

本研究所使用的未來(lái)氣候變化模式數(shù)據(jù)是CMIP6 中參與模式比較的GCM 模擬結(jié)果數(shù)據(jù)。從CMIP6 項(xiàng)目官方網(wǎng)站的下載通道(https://esgf-node.llnl.gov/search/cmip6/)收集了37 個(gè)數(shù)據(jù)較為完整的常用GCM,附表1 中為所選取的GCM 列表及相關(guān)信息(見文后電子版鏈接)。選取GCM 中數(shù)據(jù)相對(duì)較全的9 個(gè)變量: 近地面比濕度(huss)、降水量(pr)、近地面氣壓(ps)、地面下行短波輻射量(rsds)、近地面氣溫(tas)、近地面日最高氣溫(tasmax)、近地面日最低氣溫(tasmin)、近地面東向風(fēng)速(uas)和近地面北向風(fēng)速(vas)。變量信息如表1 所示。

表1 研究初步選用的未來(lái)氣候變化模式氣候參數(shù)列表及信息Table 1 List and information of preliminarily selected Global Climate Model (GCM) variables

本研究所選取的全球氣候模型 GCM 變量均源自于統(tǒng)一的實(shí)驗(yàn)場(chǎng)景設(shè)置,即在模型運(yùn)行次數(shù)(Run)、初始條件(Initialization)、物理方案(Physics)以及強(qiáng)迫數(shù)據(jù)(Forcing) 4 個(gè)方面均采用第1 個(gè)方案(均標(biāo)記為1),即r1i1p1f1。具體而言,本研究包括歷史情景(historical)以及3 個(gè)代表性濃度路徑(Shared Socioeconomic Pathways,SSPs),分別為SSP1-2.6 (SSP-126)、SSP2-4.5 (SSP245)和SSP5-8.5 (SSP585)。這些情景將被分別應(yīng)用于歷史及未來(lái)氣候變化數(shù)據(jù)的計(jì)算與分析中。此方法論的應(yīng)用旨在確保模型輸出的一致性與可比性,為評(píng)估不同氣候情景下的潛在變化提供了堅(jiān)實(shí)的基礎(chǔ)。

1.2.3 格點(diǎn)化再分析數(shù)據(jù)

GCM 輸出結(jié)果是基于大尺度(比如大陸尺度)的氣候數(shù)據(jù),空間分辨率較低,一般在百公里尺度以上。統(tǒng)計(jì)降尺度方法是提高GCM 預(yù)測(cè)精度的常用方法,該方法可以通過(guò)使用真實(shí)或近似真實(shí)的高分辨率歷史氣候數(shù)據(jù)作為參考,將低分辨率的GCM 輸出數(shù)據(jù)轉(zhuǎn)化為較為準(zhǔn)確的高分辨率氣候數(shù)據(jù)[39]。

本研究選擇ERA5 (歐洲中期天氣預(yù)報(bào)中心)數(shù)據(jù)作為高分辨率歷史氣候數(shù)據(jù)參與GCM 降尺度。ERA5 提供了大量大氣、陸地和海洋氣候變量的每小時(shí)估計(jì)值[40-41]。具體使用的數(shù)據(jù)集為ERA5-land的月平均數(shù)據(jù),它與ERA5 的其他部分相比分辨率更高,為0.1°(≈9 km)。選取的變量為1995—2020 年2 m空氣溫度、總降水量、地面下行短波輻射量、10 m東向風(fēng)速、10 m 北向風(fēng)速作為降尺度輸入數(shù)據(jù)。

1.2.4 土壤性質(zhì)數(shù)據(jù)

土壤屬性分布情況在空間尺度上極大程度地影響了糧食種植和產(chǎn)量分布。本研究將土壤屬性參數(shù)作為不隨時(shí)間變化的空間分布屬性參與模型計(jì)算。

選取國(guó)際土壤參考和資料中心(ISRIC)土壤信息數(shù)據(jù)庫(kù)中的SoilGrids250m 2.0[42]數(shù)據(jù)集作為土壤屬性參數(shù)數(shù)據(jù)(https://soilgrids.org/)。SoilGrids 的輸出結(jié)果是6 個(gè)標(biāo)準(zhǔn)深度間隔的全球土壤屬性地圖,空間分辨率為250 m。本研究中首批參數(shù)選擇5～15 cm 深度土壤的容重、陽(yáng)離子交換量、總氮、pH、有機(jī)質(zhì)含量作為土壤參數(shù)。為縮減運(yùn)算量,選擇SoilGrids 中低分辨率(5 km)子數(shù)據(jù)庫(kù)。

1.3 原始數(shù)據(jù)優(yōu)化

1.3.1 GCM 模型評(píng)價(jià)和篩選

使用泰勒?qǐng)D和技能得分進(jìn)行37 個(gè)GCM 表現(xiàn)的比較,泰勒?qǐng)D是氣候模式評(píng)價(jià)中應(yīng)用較為廣泛的一種方法[43-46],它主要通過(guò)相關(guān)系數(shù)、均方根誤差和標(biāo)準(zhǔn)差對(duì)GCM 的表現(xiàn)進(jìn)行綜合評(píng)判。對(duì)于給定的N個(gè)散點(diǎn)數(shù)據(jù)的相關(guān)系數(shù)(R)、模擬場(chǎng)(Xmi)的標(biāo)準(zhǔn)差(dm)、觀測(cè)場(chǎng)(Xoi)的標(biāo)準(zhǔn)差(dg)和均方根誤差(RMSE)的計(jì)算方法分別如下:

此外量化指標(biāo)技能得分(S)[43]可以量化模擬場(chǎng)與觀測(cè)場(chǎng)之間的相關(guān)系數(shù)以及兩者的標(biāo)準(zhǔn)差,以評(píng)估模型模擬能力的強(qiáng)弱。技能得分的計(jì)算方式如下:

式中:R為模擬場(chǎng)與觀測(cè)場(chǎng)的相關(guān)系數(shù);R0為所有模擬場(chǎng)和觀測(cè)場(chǎng)所能達(dá)到的最大相關(guān)系數(shù),即所有模型相關(guān)系數(shù)中的最大值;σf為模擬場(chǎng)標(biāo)準(zhǔn)差dm與實(shí)測(cè)場(chǎng)標(biāo)準(zhǔn)差dg之比,即 σf=dm/dg。

從37 個(gè)GCM 中選取降水量(pr)和近地面氣溫(tas)兩個(gè)變量參與評(píng)價(jià)。計(jì)算各模型技能得分,統(tǒng)計(jì)各個(gè)GCM 的pr 和tas 技能得分,并分別按照兩個(gè)變量對(duì)GCM 進(jìn)行排名并打分。打分原則是技能得分最低(第37 名)的GCM 獲得1 分,技能得分最高(第1 名)的GCM 獲得37 分,以此類推,最終將每個(gè)GCM 兩個(gè)變量的分?jǐn)?shù)平均值作為最終得分。

在氣候模式評(píng)估及其應(yīng)用的相關(guān)研究中,極少直接將表現(xiàn)最優(yōu)的模型作為唯一選擇,往往同時(shí)選取表現(xiàn)較為優(yōu)秀的多個(gè)模型組成最優(yōu)模型集合(MME)[47],并將MME 的加權(quán)平均值輸出為最終變量值。因此本研究選取模型表現(xiàn)平均得分最高且數(shù)據(jù)完整的5 個(gè)GCM 組成MME,并以其得分作為分配權(quán)重的依據(jù),計(jì)算5 個(gè)GCM 中各個(gè)變量的加權(quán)平均值,得到各個(gè)變量的最優(yōu)模型集合平均值(MMME),并作為機(jī)器學(xué)習(xí)模型的輸入數(shù)據(jù)。MMME的計(jì)算方法如下:

式中:N為最優(yōu)模型集中的模型個(gè)數(shù),即為5;Tn為第n個(gè)模型的得分; Modn為第n個(gè)模型的模擬結(jié)果。

1.3.2 變量的相關(guān)性分析和篩選

初步選擇的9 個(gè)GCM 氣候變量以及5 個(gè)土壤變量都是基于其在過(guò)往文獻(xiàn)中的使用頻率,并沒有判斷這些變量是否適用于研究區(qū)。因此,本研究需要先判斷它們是否對(duì)研究區(qū)糧食作物產(chǎn)量具有足夠影響,以及它們之間是否存在共線性過(guò)高的問題。通過(guò)去除高相關(guān)性變量以及數(shù)據(jù)降維處理等方式可以降低變量之間的共線性,去除冗余的特征變量,使機(jī)器學(xué)習(xí)模型具有更高的魯棒性以及更強(qiáng)的可解釋性[48-50]。

使用斯皮爾曼相關(guān)系數(shù)進(jìn)行變量之間的相關(guān)性分析。對(duì)于樣本容量為n的兩個(gè)樣本,其相關(guān)系數(shù)ρ的計(jì)算方法如下:

式中:xi、yi為兩個(gè)變量中的第i項(xiàng)數(shù)據(jù)。

1.4 GCM 模型降尺度

利用機(jī)器學(xué)習(xí)降尺度方法對(duì)GCM 數(shù)據(jù)進(jìn)行空間降尺度。采用隨機(jī)抽樣一致算法、K 近鄰回歸算法(K-neighbors)和線性回歸算法(LR)訓(xùn)練以重采樣GCM 數(shù)據(jù)和ERA5 格點(diǎn)化數(shù)據(jù)為自變量的數(shù)據(jù)集,訓(xùn)練集與驗(yàn)證集的數(shù)據(jù)比為3∶1。進(jìn)一步利用投票回歸算法(VR)將這3 個(gè)回歸器打包為投票機(jī)制下的多模型加權(quán)融合機(jī)器學(xué)習(xí)模型,即集成學(xué)習(xí)模型,將這種集成學(xué)習(xí)算法的輸出數(shù)據(jù)作為GCM 的降尺度數(shù)據(jù)。這種集成學(xué)習(xí)的機(jī)器學(xué)習(xí)方法可以顯著提高GCM 數(shù)據(jù)的偏差修正效果,研究表明該方法可將大部分參數(shù)修正值相對(duì)于觀測(cè)值的R2提升至0.85 以上[51]。

降尺度后的GCM 數(shù)據(jù)與ERA5 數(shù)據(jù)分辨率相同,均高于土壤和作物產(chǎn)量數(shù)據(jù)的分辨率。因此,在輸入模型前使用KD-Tree 算法將土壤和作物產(chǎn)量因子重采樣到GCM 數(shù)據(jù)分辨率上,從而保證輸入模型的數(shù)據(jù)點(diǎn)一一對(duì)應(yīng)。

1.5 模型訓(xùn)練

對(duì)MMME進(jìn)行變量篩選后,利用選出的變量計(jì)算年平均值,并對(duì)4—5 月(貝爾格季作物的主要生長(zhǎng)期)和7—9 月(梅赫季作物的主要生長(zhǎng)期)的上述數(shù)據(jù)取平均值,以保證每種情景下各參數(shù)均具有年平均、貝爾格季平均和梅赫季平均作為氣候自變量,與土壤變量共同組成自變量數(shù)據(jù)集。分別將10 個(gè)糧食作物產(chǎn)量數(shù)據(jù)集(5 種作物梅赫季產(chǎn)量和貝爾格季產(chǎn)量)作為因變量進(jìn)行預(yù)處理,與自變量數(shù)據(jù)形成數(shù)據(jù)庫(kù),用于訓(xùn)練機(jī)器學(xué)習(xí)模型。

初步選取了3 類共6 種機(jī)器學(xué)習(xí)模型進(jìn)行訓(xùn)練,分別為: 1)梯度提升算法(Boosting),包括直方圖梯度提升回歸(HGB)、極端梯度提升隨機(jī)森林(XGBRF)算法以及輕梯度提升樹(LGBM)算法;2)隨機(jī)森林算法,包括隨機(jī)森林(RF)以及極限樹(ET)算法;3) K 近鄰(K-neighbors)算法。

將上述6 種算法分別基于梅赫季苔麩產(chǎn)量數(shù)據(jù)進(jìn)行訓(xùn)練后,使用統(tǒng)計(jì)指標(biāo)評(píng)價(jià)模型表現(xiàn)。然后采用集成學(xué)習(xí)的思想,取表現(xiàn)較好的3 個(gè)算法進(jìn)行打包訓(xùn)練。這里采用的集成學(xué)習(xí)方法是堆疊,即在3個(gè)算法的基礎(chǔ)上再使用一個(gè)線性回歸模型進(jìn)行二次計(jì)算。

參與訓(xùn)練的數(shù)據(jù)被分為訓(xùn)練集(75%)和驗(yàn)證集(25%),其中訓(xùn)練集用于模型訓(xùn)練,驗(yàn)證集用于和模型結(jié)果比對(duì)以判定模型表現(xiàn)。使用模型判決系數(shù)(R2)、平均絕對(duì)誤差(MAE)和解釋方差評(píng)分(EVS)作為評(píng)價(jià)模型效果和判斷模型可信度的依據(jù)。同時(shí)繪制每種模型在驗(yàn)證過(guò)程中對(duì)真實(shí)產(chǎn)量和模型預(yù)測(cè)產(chǎn)量的散點(diǎn)圖,從而為模型精度的評(píng)價(jià)提供直觀判斷。由于數(shù)據(jù)量較大,為方便觀察僅隨機(jī)選取了驗(yàn)證集中的10 000 個(gè)數(shù)據(jù)點(diǎn)進(jìn)行散點(diǎn)圖的繪制。

在6 種機(jī)器學(xué)習(xí)模型的橫向比較中,使用梅赫季苔麩產(chǎn)量數(shù)據(jù)對(duì)6 個(gè)模型進(jìn)行訓(xùn)練和驗(yàn)證,選擇表現(xiàn)較高的3 個(gè)模型進(jìn)行堆疊訓(xùn)練,將訓(xùn)練完成后的最終模型分別應(yīng)用于10 個(gè)因變量糧食數(shù)據(jù)集,進(jìn)行針對(duì)10 個(gè)不同作物數(shù)據(jù)模型精度的驗(yàn)證和評(píng)價(jià),方法與模型橫向比較方法相同。

1.6 產(chǎn)量預(yù)測(cè)

使用訓(xùn)練完成的堆疊模型對(duì)自變量數(shù)據(jù)集進(jìn)行預(yù)測(cè),并輸出預(yù)測(cè)結(jié)果。共輸出30 個(gè)產(chǎn)量預(yù)測(cè)結(jié)果,即10 個(gè)作物數(shù)據(jù)在3 個(gè)未來(lái)情景下的不同產(chǎn)量結(jié)果。進(jìn)一步計(jì)算每個(gè)未來(lái)情景下各作物類別的產(chǎn)量變化,并交叉引用同一作物在不同未來(lái)情景下的產(chǎn)量,分析2021—2050 年埃塞俄比亞主要糧食作物產(chǎn)量趨勢(shì)以及氣候變化對(duì)其的影響。

2 結(jié)果與分析

2.1 GCM 模型評(píng)價(jià)結(jié)果

37 個(gè)GCM 的pr 和tas 變量泰勒?qǐng)D如圖2 所示,附表2 提供了37 個(gè)GCM 兩個(gè)變量下的技能得分、表現(xiàn)評(píng)分以及綜合評(píng)分(見文后電子版鏈接)。

表2 6 種不同機(jī)器學(xué)習(xí)模型的模擬梅赫季苔麩產(chǎn)量的表現(xiàn)Table 2 Prediction performance of teff yield (Meher season) of 6 different machine learning models

從圖2 可以看出在研究區(qū)范圍內(nèi),GCM 之間的表現(xiàn)具有較大差距,這也符合區(qū)域尺度全球氣候模式表現(xiàn)的一般情況。有研究顯示CMCC-CM2-SR5、AWI-ESM-1-1-LR 兩種模式在東非地區(qū)表現(xiàn)優(yōu)良[52],這與本研究結(jié)論一致;EC-Earth系列、MPI-EMS 系列模型在非洲、東亞中亞等地都被證明具有較好的預(yù)測(cè)表現(xiàn)[53],在本研究區(qū)范圍內(nèi)也呈現(xiàn)出相同結(jié)果。FIO-ESM-2-0 模型在非洲以及中亞地區(qū)的模型比較研究中表現(xiàn)并不優(yōu)秀,這與本研究結(jié)論相反,原因可能在于時(shí)間范圍劃定以及用于比對(duì)的數(shù)據(jù)庫(kù)準(zhǔn)確度等不同,ERA5 作為再分析數(shù)據(jù)與眾多研究中使用的實(shí)際觀測(cè)數(shù)據(jù)不同,出現(xiàn)離群值的概率較低,方差分布更平均,但會(huì)與觀測(cè)值有所偏差。

同一GCM 在兩個(gè)變量中的表現(xiàn)也不盡相同。如AWI-CM-1-1-MR、NorESM2-MM、MPI-ESM-1-2-HAM 等模型在兩個(gè)參數(shù)下的排名相差甚至超過(guò)20,在綜合加權(quán)評(píng)分機(jī)制的影響下,這種變量間表現(xiàn)差異較大的GCM 也不在研究考慮選取的范圍內(nèi)。

最終選擇的最優(yōu)模型集合包括以下5 個(gè)GCM:CMCC-CM2-SR5、MPI-ESM1-2-LR、EC-Earth3-Veg-LR、EC-Earth3-Veg 和MPI-ESM1-2-HR,其綜合得分在37 個(gè)GCM 中分別排名2、3、4、8 和9,一些GCM 排名高于上述GCM 但并未入選,原因是在它們?cè)O(shè)定的實(shí)驗(yàn)情景、初始場(chǎng)以及時(shí)間分辨率下不包含全部的9 個(gè)變量,無(wú)法參加后續(xù)變量相關(guān)性分析和模型訓(xùn)練。

2.2 變量篩選結(jié)果

14 個(gè)變量和5 種梅赫季作物數(shù)據(jù)的斯皮爾曼相關(guān)性熱圖如圖3 所示。在氣候變量中,近地面壓強(qiáng)與3 項(xiàng)氣溫指標(biāo)之間呈現(xiàn)強(qiáng)正相關(guān)關(guān)系,這是由于大氣壓強(qiáng)與氣溫線性相關(guān)。在土壤變量中,各變量之間的相關(guān)性明顯高于氣候變量,高相關(guān)性出現(xiàn)在與容重以及土壤有機(jī)碳相交叉的行列中,例如容重與有機(jī)碳相關(guān)系數(shù)為-0.84,呈高度負(fù)相關(guān);容重與總氮和pH 相關(guān)系數(shù)的絕對(duì)值均在0.7 以上,土壤有機(jī)碳與總氮和pH 相關(guān)系數(shù)的絕對(duì)值也均高于0.7,且呈高度負(fù)相關(guān)。

圖3 14 個(gè)變量和5 種梅赫季作物產(chǎn)量的斯皮爾曼相關(guān)性熱圖Fig.3 Spearman correlation heat map with 14 variables and yield of 5 crops in Meher season

從不同變量組之間的關(guān)系來(lái)看,降水量與土壤pH 之間有較強(qiáng)的負(fù)相關(guān)性,原因是土壤中微生物分解以及根系呼吸作用會(huì)產(chǎn)生CO2等酸性物質(zhì),而降水會(huì)將參與中和酸性的鹽離子通過(guò)淋溶作用帶走,使得土壤pH 降低。作物對(duì)于以上全部變量的相關(guān)性反應(yīng)較為穩(wěn)定,相關(guān)系數(shù)范圍在0.13～0.56 之間,未呈現(xiàn)出明顯的趨勢(shì)。這同樣表明作物產(chǎn)量與上述各因子間具有復(fù)雜的內(nèi)在聯(lián)系,但目前理論無(wú)法充分解釋這種關(guān)系。

通過(guò)變量之間的相關(guān)性分析,決定剔除氣候變量中的近地面氣壓、近地面平均日最高與最低氣溫、土壤變量中的容重和土壤有機(jī)碳含量。

2.3 機(jī)器學(xué)習(xí)模型評(píng)價(jià)

2.3.1 不同機(jī)器學(xué)習(xí)模型的比較結(jié)果

對(duì)6 種不同的機(jī)器學(xué)習(xí)模型針對(duì)同一作物的模擬結(jié)果進(jìn)行評(píng)價(jià),糧食產(chǎn)量數(shù)據(jù)選擇梅赫季的苔麩產(chǎn)量。評(píng)價(jià)指標(biāo)的結(jié)果如表2 所示。

從結(jié)果中可以看出綜合表現(xiàn)較好的模型為極端梯度提升隨機(jī)森林、隨機(jī)森林以及極限樹模型。在6 個(gè)模型中表現(xiàn)較差的是直方圖梯度提升和K 近鄰,其中K 近鄰是應(yīng)用最為廣泛且最簡(jiǎn)單的機(jī)器學(xué)習(xí)模型之一,其計(jì)算路徑較短,善于求解線性問題,但在面對(duì)復(fù)雜問題時(shí)解釋性不足;直方圖梯度提升作為梯度提升決策樹模型,雖然其R2和EVS 表現(xiàn)最差,但卻有較低的MAE,說(shuō)明該模型雖然盡量逼近了因變量變化范圍,但是對(duì)因變量變化趨勢(shì)的把握較差,其預(yù)測(cè)結(jié)果幾乎不能解釋因變量方差變化規(guī)律。

整體來(lái)看,在糧食產(chǎn)量長(zhǎng)期預(yù)測(cè)方面,較為復(fù)雜、運(yùn)算量大的機(jī)器學(xué)習(xí)模型的模擬能力占優(yōu)勢(shì)。本研究最終選擇參與堆疊的模型是極端梯度提升隨機(jī)森林、隨機(jī)森林以及極限樹模型。

2.3.2 不同作物模型的比較結(jié)果

將上述3 種模型使用線性回歸模型進(jìn)行堆疊,針對(duì)10 種作物產(chǎn)量數(shù)據(jù)分別進(jìn)行訓(xùn)練,最終得出10 個(gè)作物產(chǎn)量模型,對(duì)它們的模擬能力分別進(jìn)行評(píng)價(jià),評(píng)價(jià)指標(biāo)結(jié)果如表3 所示。

表3 最終模型對(duì)不同作物產(chǎn)量的模擬表現(xiàn)Table 3 Prediction performance of the stacked model for yields of different crops

整體來(lái)看,機(jī)器學(xué)習(xí)模型的預(yù)測(cè)結(jié)果較為準(zhǔn)確,除貝爾格季大麥模型外,所有模型的R2均在可信范圍內(nèi),且與實(shí)際值的誤差都控制在了1 t·hm-2以內(nèi)。其中梅赫季的小麥和玉米產(chǎn)量表現(xiàn)較為突出,R2在0.8 以上。此外,堆疊后模型的表現(xiàn)明顯優(yōu)于單一模型,說(shuō)明集成學(xué)習(xí)的方法將不同模型的優(yōu)點(diǎn)進(jìn)行了整合,顯著提升了模型表現(xiàn)。對(duì)于不同種類作物,梅赫季作物模型表現(xiàn)要明顯優(yōu)于貝爾格季作物模型表現(xiàn),這是由于貝爾格季因變量作物產(chǎn)量數(shù)據(jù)質(zhì)量較差,缺失數(shù)據(jù)較多,經(jīng)預(yù)處理后參與訓(xùn)練的數(shù)據(jù)量?jī)H為梅赫季的60%左右。

貝爾格季大麥模型表現(xiàn)較差,其R2、EVS 均低于0.4～1.0 的可信范圍,同時(shí)MAE 高達(dá)20,甚至超出2020 年產(chǎn)量的兩倍。因此,在后續(xù)分析中不再考慮貝爾格季大麥模型。

2.4 糧食產(chǎn)量預(yù)測(cè)結(jié)果及分析

2.4.1 埃塞俄比亞1995—2020 年糧食產(chǎn)量情況

本研究分析了埃塞俄比亞地區(qū)1995—2020 年糧食產(chǎn)量變化情況(圖4),結(jié)果表明,研究區(qū)中部地區(qū)種植的作物品種最為全面,但單一作物產(chǎn)量低于周邊地區(qū)。這是由于中部地區(qū)氣候條件較好,氣候帶適宜種植,降水條件優(yōu)渥,鮮有干旱發(fā)生,因此多數(shù)作物都可以在此地種植。但是由于中部地區(qū)地勢(shì)較高,多見山地丘陵地形,使適宜種植作物的區(qū)域變得碎片化,因此限制了糧食總產(chǎn)量的增加。研究區(qū)外圍地區(qū)具有各自的作物種植特色,其中苔麩多種植于西北部地區(qū),玉米多種植于東南及東北部地區(qū),大麥及小麥在除了極東部和極西部地區(qū)均有種植,高粱種植則較為分散。

圖4 2020 年埃塞俄比亞5 種主要糧食作物產(chǎn)量(統(tǒng)計(jì)數(shù)據(jù))Fig.4 Yields of 5 main staple crops in 2020 in Ethiopia (statistical data)

玉米是埃塞俄比亞種植最為廣泛的作物,也是全國(guó)總產(chǎn)量最高的糧食作物,這是由于玉米相較于其他糧食作物較為耐旱,在非洲充足的光熱條件下,玉米生長(zhǎng)和種植也相較于其他作物更為容易。反觀小麥,雖然保證了充足的光熱,但干旱導(dǎo)致的水分不足嚴(yán)重影響了小麥的生長(zhǎng)發(fā)育,因此,小麥在埃塞俄比亞主要糧食作物中總產(chǎn)量最低。這種情況也出現(xiàn)在其他的東非國(guó)家,它們都重度依賴小麥進(jìn)口。

1995—2020 年埃塞俄比亞糧食產(chǎn)量變化如圖5所示,整體來(lái)看梅赫季糧食產(chǎn)量整體呈緩慢上升的趨勢(shì),部分地區(qū)有所下降;貝爾格季苔麩、玉米、小麥、大麥平均增長(zhǎng)約0.6 t·hm-2,高粱則有所下降。所有作物在梅赫季產(chǎn)量增減最為不均,中西部地區(qū)增長(zhǎng)近3 t·hm-2,而在北方和南方則出現(xiàn)了近2 t·hm-2左右的下降,這可能是由于自20 世紀(jì)80 年代起東非地區(qū)流行的干旱所致,在干旱較為嚴(yán)重且降水不足的北方和南方,其直接展現(xiàn)出對(duì)作物生長(zhǎng)的抑制[54]。在中西部地區(qū),氣候條件變暖以及高種植密度帶來(lái)的耕作技術(shù)提升緩和了這種狀況。

圖5 1995—2020 年間埃塞俄比亞5 種主要糧食作物產(chǎn)量變化(統(tǒng)計(jì)數(shù)據(jù))Fig.5 Changes of yield of 5 main staple crops in Ethiopia from 1995 to 2020 (statistical data)

2.4.2 埃塞俄比亞2021—2050 年作物產(chǎn)量預(yù)測(cè)

從梅赫季產(chǎn)量變化(圖6)可以看出,苔麩產(chǎn)量增減在0.6 t·hm-2之內(nèi),增長(zhǎng)主要發(fā)生在中西部,大幅減產(chǎn)出現(xiàn)在東南沙漠地區(qū)。玉米的強(qiáng)烈波動(dòng)主要出現(xiàn)在東南沙漠地區(qū),其大部分區(qū)域都出現(xiàn)了2 t·hm-2以內(nèi)的產(chǎn)量增長(zhǎng)。在絕大多數(shù)地區(qū)小麥產(chǎn)量波動(dòng)幅度都在2 t·hm-2以內(nèi),但SSP126 情景下東南地區(qū)產(chǎn)量下降而另外兩種情景下卻在上升。大麥的產(chǎn)量變化情況在空間上類似于苔麩,但其產(chǎn)量變化集中在1.5 t·hm-2以內(nèi),劇烈變化出現(xiàn)在中部地區(qū)。對(duì)于高粱來(lái)說(shuō),SSP126 和SSP585 情景下的變化情況相似,以東南部1.5 t·hm-2左右的增產(chǎn)為主,SSP245 情景下的產(chǎn)量變化則較為平均。

圖6 2021—2050 年SSP126、SSP245、SSP585 情景下梅赫季5 種主要糧食作物產(chǎn)量變化量Fig.6 Yield changes of 5 main staple crops under SSP126,SSP245 and SSP585 scenarios in Meher season from 2021 to 2050

總的來(lái)看,未來(lái)30 年梅赫季5 種主要糧食作物產(chǎn)量增幅以＜2 t·hm-2為主(圖6),不同作物產(chǎn)量變化的空間特征呈現(xiàn)出一定的互補(bǔ)性,例如東南地區(qū)苔麩和大麥產(chǎn)量的下降與小麥和高粱產(chǎn)量的上升。

貝爾格季產(chǎn)量變化可以看出(圖7),苔麩產(chǎn)量變化與梅赫季極為相似,但是其變化量為梅赫季的1/2左右,這與貝爾格季本身產(chǎn)量較低有關(guān)。對(duì)于玉米來(lái) 說(shuō),SSP245 和SSP585情景出現(xiàn)了全國(guó)范圍 內(nèi)1 t·hm-2左右的增產(chǎn),唯有阿法爾地區(qū)出現(xiàn)大幅減產(chǎn),但在SSP126 情景下中南部地區(qū)將出現(xiàn)0.5～0.8 t·hm-2左右的減產(chǎn)。小麥增產(chǎn)集中在中部地區(qū),但中東部和南部地區(qū)在SSP126 情景下出現(xiàn)了0.4 t·hm-2左右的減產(chǎn)。SSP126 情景下的減產(chǎn)情況同樣延續(xù)到了高粱的產(chǎn)量: 全國(guó)大規(guī)模出現(xiàn)＞1 t·hm-2的減產(chǎn),中部地區(qū)的1 t·hm-2的增產(chǎn)發(fā)生在其他兩個(gè)情景中。

圖7 2021—2050 年SSP126、SSP245、SSP585 情景下貝爾格季4 種主要糧食作物產(chǎn)量變化量Fig.7 Yield changes of 4 main staple crops under SSP126,SSP245 and SSP585 scenarios in Belg season from 2021 to 2050

貝爾格季糧食變化情況反映出SSP126 情景與其他兩種情景下的較大差異,即明顯減產(chǎn)的趨勢(shì)。這與部分前人研究結(jié)果類似[2,55],一般是由于環(huán)境條件改變使作物脫離原有適宜生長(zhǎng)范圍。貝爾格季主要播種時(shí)間為夏季,干旱區(qū)夏季種植習(xí)慣和作物生長(zhǎng)規(guī)律與SSP126 情景提供的良好生態(tài)條件不符,可能會(huì)導(dǎo)致部分作物減產(chǎn)。SSP245 情景下作物在中部及南部地區(qū)的產(chǎn)量將普遍增加,這符合過(guò)去26 年作物產(chǎn)量發(fā)展趨勢(shì),進(jìn)一步驗(yàn)證了SSP245 情景與現(xiàn)今社會(huì)發(fā)展情況的一致性。

進(jìn)一步分析往年貝爾格季結(jié)果可以發(fā)現(xiàn)主要糧食作物種植主要集中在中部和南部地區(qū),只有小麥種植范圍向東南方向延伸,而在實(shí)際情況中,相對(duì)于梅赫季存在大量零散種植的情況,貝爾格季糧食種植基本局限于中部和南部地區(qū)。雖然模型給出了貝爾格季其他區(qū)域的預(yù)測(cè)產(chǎn)量,但與實(shí)際情況不符,其參考性較差。因此,具有較高參考性的變化趨勢(shì)應(yīng)該主要在中部和南部地區(qū),這些地區(qū)的產(chǎn)量將隨著未來(lái)情景從SSP585 到SSP126 的轉(zhuǎn)變而減少。

3 討論

總體來(lái)看,2050 年SSP245 情景下作物產(chǎn)量普遍低于SSP126 情景,而SSP585 情景作物產(chǎn)量相較于SSP245 情景則更高。這在梅赫季作物產(chǎn)量變化中較為明顯,即SSP126 情景下作物增產(chǎn)占主導(dǎo),而在另外兩種情景下作物減產(chǎn)規(guī)模有所升高。這說(shuō)明在3個(gè)情景中,SSP245 情景可能是糧食產(chǎn)量增長(zhǎng)最為緩慢的情景。這種現(xiàn)象意味著SSP245 情景對(duì)于糧食作物生產(chǎn)來(lái)說(shuō),表現(xiàn)出情景變化間的拐點(diǎn)特征。本研究認(rèn)為,這種情景間產(chǎn)量先減后增的過(guò)程是由情景所代表的氣候原因造成的。對(duì)于地處干旱區(qū)的埃塞俄比亞來(lái)說(shuō),SSP245 和SSP126 情景下生態(tài)環(huán)境的改善(即干旱的緩解)意味著各種作物的生存環(huán)境都將得到改善,從而帶來(lái)增產(chǎn)的結(jié)果。SSP585 情景中人為活動(dòng)帶來(lái)的溫室氣體排放將加快全球變暖的進(jìn)程,大量研究表明一定程度下的全球變暖可以為作物生長(zhǎng)提供更充足的光合作用原料,這在一定程度上會(huì)帶來(lái)作物的增產(chǎn),這也是SSP245 情景糧食產(chǎn)量相較于另外兩種情景增長(zhǎng)更慢的原因。

有研究曾提出SSP585 情景下的農(nóng)業(yè)發(fā)展會(huì)受到社會(huì)發(fā)展壓力的影響,未來(lái)社會(huì)形勢(shì)從優(yōu)質(zhì)和諧社會(huì)向復(fù)雜矛盾社會(huì)轉(zhuǎn)變時(shí),作物產(chǎn)量增減趨勢(shì)將會(huì)變得更加突出[56],這與本研究結(jié)果相同。這是由于在社會(huì)復(fù)雜化的過(guò)程中,社會(huì)對(duì)糧食供給結(jié)構(gòu)的優(yōu)化程度需求更大,這不僅僅意味著人們對(duì)糧食需求量的增加,同時(shí)代表著人們希望糧食生產(chǎn)結(jié)構(gòu)更加合理化,在最優(yōu)的田間管理?xiàng)l件下,盡可能地利用自然條件的便利最大化生產(chǎn),以對(duì)沖日益加深的社會(huì)矛盾。出于上述原因,在SSP585 情景下,多數(shù)作物呈現(xiàn)出的變化趨勢(shì)與SSP245 情景大致相同,但幅度更大,即在相同地區(qū)更劇烈地增產(chǎn)或減產(chǎn);而相比之下,SSP126 情景中糧食生產(chǎn)結(jié)構(gòu)性改革和社會(huì)生產(chǎn)能力的再分配似乎不那么緊迫,這使得極端值在相應(yīng)圖像中出現(xiàn)的頻率更低。

4 結(jié)論

本文利用機(jī)器學(xué)習(xí)算法將埃塞俄比亞糧食產(chǎn)量統(tǒng)計(jì)數(shù)據(jù)與CMIP6 未來(lái)氣候模式數(shù)據(jù)相結(jié)合,訓(xùn)練模型并預(yù)測(cè)了2021—2050 年不同未來(lái)社會(huì)情景下埃塞俄比亞兩個(gè)主要生產(chǎn)季的5 種糧食作物產(chǎn)量,得出以下結(jié)論:

1)本研究選取的37 個(gè)GCM 中,CMCC-CM2-SR5、MPI-ESM1-2-LR、EC-Earth3-Veg-LR、ECEarth3-Veg、MPI-ESM1-2-HR 在埃塞俄比亞地區(qū)的模擬效果較好,且數(shù)據(jù)完整。采用泰勒?qǐng)D和技能得分結(jié)合排序評(píng)分和加權(quán)平均算法可以方便快捷地評(píng)價(jià)GCM 的表現(xiàn),綜合多個(gè)優(yōu)秀模型可以進(jìn)一步減少誤差。

2)在使用氣候變量和土壤變量預(yù)測(cè)糧食產(chǎn)量的機(jī)器學(xué)習(xí)模型訓(xùn)練過(guò)程中,極端梯度提升隨機(jī)森林、隨機(jī)森林、極限樹3 種模型的表現(xiàn)要優(yōu)于直方圖梯度提升、輕梯度提升以及K 近鄰。

3)未來(lái)30 年間,梅赫季5 種主要糧食作物產(chǎn)量變化以＜2 t·hm-2為主,不同作物空間變化存在一定程度互補(bǔ);在SSP126 情景下貝爾格季出現(xiàn)減產(chǎn)較多,其他兩種情景則以增產(chǎn)為主,這可能與SSP126 情景下的生態(tài)環(huán)境改善與生產(chǎn)季種植習(xí)慣的不匹配有關(guān)。

4)從SSP126 情景到SSP585 情景,隨著社會(huì)矛盾加劇和人為環(huán)境惡化,糧食作物生產(chǎn)結(jié)構(gòu)改革和生產(chǎn)資源再分配的需求增加,導(dǎo)致各種作物的生產(chǎn)力向新出現(xiàn)的適宜地區(qū)轉(zhuǎn)移。研究區(qū)在SSP126 和SSP585 情景下分別會(huì)因?yàn)槊庥诟珊岛蜏厥倚?yīng)加劇獲得更高的糧食作物生產(chǎn)力。