基于機(jī)器學(xué)習(xí)與氣象災(zāi)害指標(biāo)的蘋(píng)果相對(duì)氣象產(chǎn)量預(yù)測(cè)

羅 琦 茹曉雅 姜 元 馮 浩 于 強(qiáng) 何建強(qiáng)

(1.西北農(nóng)林科技大學(xué)旱區(qū)農(nóng)業(yè)水土工程教育部重點(diǎn)實(shí)驗(yàn)室,陜西楊凌 712100;2.西北農(nóng)林科技大學(xué)中國(guó)旱區(qū)節(jié)水農(nóng)業(yè)研究院,陜西楊凌 712100;3.中國(guó)科學(xué)院水利部水土保持研究所黃土高原土壤侵蝕與旱地農(nóng)業(yè)國(guó)家重點(diǎn)實(shí)驗(yàn)室,陜西楊凌 712100;4.陜西省氣象局秦嶺和黃土高原生態(tài)環(huán)境氣象重點(diǎn)實(shí)驗(yàn)室,西安 710016)

0 引言

目前,我國(guó)黃土高原是全世界最大的蘋(píng)果(MaluspumilaMill.)種植區(qū)之一,其蘋(píng)果種植面積和產(chǎn)量分別占全球的25.2%和26.3%[1-2]。由于蘋(píng)果種植業(yè)的經(jīng)濟(jì)效益較高,黃土高原大部分傳統(tǒng)農(nóng)田都已改建為蘋(píng)果園,蘋(píng)果產(chǎn)業(yè)已成為當(dāng)?shù)剞r(nóng)村經(jīng)濟(jì)的支柱產(chǎn)業(yè)之一,是當(dāng)?shù)剞r(nóng)民的主要收入來(lái)源。因此,及時(shí)準(zhǔn)確地預(yù)測(cè)蘋(píng)果產(chǎn)量對(duì)當(dāng)?shù)靥O(píng)果產(chǎn)業(yè)的發(fā)展具有重要意義。

黃土高原大部分地區(qū)光熱資源豐富、晝夜溫差大,對(duì)蘋(píng)果樹(shù)的生長(zhǎng)發(fā)育和結(jié)果極為有利[3]。降水、氣溫和太陽(yáng)輻射是影響蘋(píng)果生長(zhǎng)發(fā)育的重要?dú)庀笠?但氣象災(zāi)害也會(huì)對(duì)蘋(píng)果產(chǎn)量形成過(guò)程造成極大的影響。黃土高原地區(qū)干旱少雨,屬于典型生態(tài)脆弱區(qū),容易發(fā)生各種氣象災(zāi)害。已有相關(guān)研究分析了各類(lèi)氣象災(zāi)害對(duì)蘋(píng)果產(chǎn)量的影響。例如,RODRIGO[4]發(fā)現(xiàn)春季花期凍害很少導(dǎo)致果樹(shù)死亡,但會(huì)嚴(yán)重?fù)p害蘋(píng)果花芽,是限制蘋(píng)果產(chǎn)量的關(guān)鍵因素;韓文靜等[5]采用連陰雨指標(biāo)分析發(fā)現(xiàn)蘋(píng)果著色成熟期和采收期受到連陰雨危害的程度高于其他生育期;戴安然等[6]研究結(jié)果表明,旱災(zāi)是限制中國(guó)蘋(píng)果生產(chǎn)發(fā)展的主要因素之一,會(huì)使果樹(shù)正常生長(zhǎng)受到抑制,導(dǎo)致葉片萎蔫、落果、枝條抽干等。此外,干旱指標(biāo)中標(biāo)準(zhǔn)化降水蒸散發(fā)指數(shù)(Standardized precipitation evapotranspiration index,SPEI)對(duì)降水和潛在蒸散發(fā)同樣敏感,1—6個(gè)月尺度的SPEI可用于評(píng)估農(nóng)業(yè)干旱,且該指標(biāo)還被廣泛應(yīng)用于不同作物的產(chǎn)量預(yù)測(cè)模型中[7]。

為了量化氣象因子和氣象災(zāi)害因子對(duì)蘋(píng)果產(chǎn)量的影響,有必要建立上述因子和蘋(píng)果產(chǎn)量之間的數(shù)學(xué)映射關(guān)系。但是不同氣象因子和氣象災(zāi)害因子之間存在較強(qiáng)的共線性,并且這些因子與蘋(píng)果產(chǎn)量之間通常為非線性關(guān)系,將傳統(tǒng)的線性模型應(yīng)用于作物產(chǎn)量預(yù)測(cè)往往不能獲得理想的預(yù)測(cè)結(jié)果。近年來(lái),機(jī)器學(xué)習(xí)算法由于其眾多優(yōu)勢(shì)而逐漸被廣泛應(yīng)用于作物產(chǎn)量預(yù)測(cè)研究[8]。這類(lèi)方法從訓(xùn)練數(shù)據(jù)集中提取信息建立統(tǒng)計(jì)模型,隨后根據(jù)測(cè)試數(shù)據(jù)集評(píng)估模型模擬精度,從而在自變量和非自變量之間建立非線性關(guān)系[9-10]。相關(guān)研究表明機(jī)器學(xué)習(xí)算法可在蘋(píng)果產(chǎn)量預(yù)測(cè)中取得較好的效果。例如,LI等[11]利用氣象數(shù)據(jù)和支持向量機(jī)(Support vector machine,SVM)算法對(duì)陜西省28個(gè)基地縣的蘋(píng)果進(jìn)行了估產(chǎn),但研究區(qū)域只限制于陜西省;KHAN等[12]基于農(nóng)業(yè)生產(chǎn)總值建立了基于列文伯格馬夸爾特(Levenberg-Marquardt optimization,LM)、尺度共軛梯度(Scale conjugate gradient back propagation,SCG)、貝葉斯正則化反向傳播神經(jīng)網(wǎng)絡(luò)(Bayesian regularization back propagation,BRBP)算法的蘋(píng)果產(chǎn)量預(yù)報(bào)模型,但其選擇的模型輸入變量并不包括影響蘋(píng)果產(chǎn)量形成的環(huán)境因素;景輝等[13]采用氣象數(shù)據(jù)和多元線性回歸(Multiple linear regression,MLR)、反向傳播神經(jīng)網(wǎng)絡(luò)(Back propagation neural network,BPNN)算法建立了蘋(píng)果產(chǎn)量早期預(yù)測(cè)模型,但并未考慮氣象災(zāi)害因子對(duì)蘋(píng)果產(chǎn)量的影響,且只考慮了各生育期內(nèi)氣象因子對(duì)蘋(píng)果產(chǎn)量的影響。總體而言,上述研究的不足在于研究區(qū)域往往較小,所選擇的蘋(píng)果產(chǎn)量預(yù)測(cè)模型輸入特征變量較少考慮蘋(píng)果生長(zhǎng)季內(nèi)不同月份氣象因子的影響,以及氣象災(zāi)害因子對(duì)蘋(píng)果產(chǎn)量的影響。

本文以我國(guó)黃土高原蘋(píng)果產(chǎn)區(qū)為研究區(qū)域,基于氣象、空間和氣象災(zāi)害3類(lèi)特征變量與產(chǎn)區(qū)蘋(píng)果單產(chǎn)統(tǒng)計(jì)數(shù)據(jù),選擇梯度提升樹(shù)(Gradient boosting decision tree,GBDT)、SVM和BRBP神經(jīng)網(wǎng)絡(luò)3種機(jī)器學(xué)習(xí)算法和MLR算法,構(gòu)建黃土高原蘋(píng)果相對(duì)氣象產(chǎn)量預(yù)測(cè)模型。首先,采用斯皮爾曼相關(guān)性分析確定影響黃土高原蘋(píng)果相對(duì)氣象產(chǎn)量的最重要?dú)庀筇卣髯兞?然后以最重要?dú)庀筇卣髯兞?、氣象?zāi)害特征變量和空間特征變量為輸入變量,選擇蘋(píng)果相對(duì)氣象產(chǎn)量預(yù)測(cè)模型的最佳模型輸入特征變量組合,最后基于不同生育期和生長(zhǎng)季內(nèi)各月份的最佳模型輸入特征變量組合,分析不同模型預(yù)測(cè)蘋(píng)果產(chǎn)量的提前期,以期為黃土高原蘋(píng)果產(chǎn)量早期預(yù)測(cè)提供科學(xué)依據(jù)和技術(shù)參考。

1 材料與方法

1.1 研究區(qū)域

研究區(qū)域?yàn)槲覈?guó)黃土高原蘋(píng)果產(chǎn)區(qū),主要包括河南、山西、陜西和甘肅4省在內(nèi)的共86個(gè)蘋(píng)果基地縣(34°22′～37°42′N(xiāo),104°53′～112°46′E,海拔84～4 216 m,圖1)。黃土高原屬于半干旱大陸性季風(fēng)氣候,是我國(guó)最大的蘋(píng)果優(yōu)生區(qū),其多項(xiàng)氣象因素能滿足優(yōu)質(zhì)蘋(píng)果生長(zhǎng)需求。研究區(qū)內(nèi)年均溫為6～16℃,年降水量為201～1 010 mm,年日照時(shí)數(shù)為1 294～2 900 h。

圖1 黃土高原86個(gè)蘋(píng)果生產(chǎn)基地縣和氣象站點(diǎn)分布圖

1.2 數(shù)據(jù)來(lái)源與處理

1981—2017年共37年的逐日氣象觀測(cè)數(shù)據(jù)(包括降水量、平均氣溫、最高氣溫、最低氣溫、日照時(shí)數(shù)及空氣相對(duì)濕度等)來(lái)源于中國(guó)氣象科學(xué)數(shù)據(jù)共享網(wǎng)(CMDSSS,http:∥cdc.cma.gov.cn/)。氣象觀測(cè)站空間位置信息(包括經(jīng)度、緯度、海拔)來(lái)自各縣市內(nèi)的氣象觀測(cè)站。河南、山西、陜西和甘肅共86個(gè)縣級(jí)行政區(qū)的蘋(píng)果單產(chǎn)數(shù)據(jù)來(lái)自《河南統(tǒng)計(jì)年鑒》[14]、《山西統(tǒng)計(jì)年鑒》[15]、《陜西統(tǒng)計(jì)年鑒》[16]和《甘肅發(fā)展年鑒》[17],以及各市人民政府依申請(qǐng)公開(kāi)數(shù)據(jù)(表1)。本文剔除了各縣市單產(chǎn)缺失年份,共獲得1 852個(gè)有效單產(chǎn)數(shù)據(jù)。

表1 黃土高原蘋(píng)果產(chǎn)區(qū)86個(gè)生產(chǎn)基地縣產(chǎn)量記錄年份

數(shù)據(jù)歸一化處理可消除不同數(shù)據(jù)的量綱影響,提升模型各項(xiàng)性能[18]。因此,模型建立前統(tǒng)一對(duì)輸入數(shù)據(jù)進(jìn)行歸一化處理,最后將模型輸出值進(jìn)行反歸一化處理。計(jì)算公式為

(1)

式中xn——?dú)w一化后的輸出數(shù)據(jù)

x——待歸一化的輸入數(shù)據(jù)

xmax——待歸一化輸入數(shù)據(jù)最大值

xmin——待歸一化輸入數(shù)據(jù)最小值

1.3 模型輸入特征變量和目標(biāo)變量

1.3.1模型輸入特征變量?jī)?yōu)選

依據(jù)文獻(xiàn)[19-20],黃土高原產(chǎn)區(qū)的蘋(píng)果全生育期為本年11月至次年10月,可主要?jiǎng)澐譃槁淙~期(本年11月)、休眠期(本年12月—次年2月)、萌芽幼果期(次年3—5月)、果實(shí)膨大期(次年6—8月)和著色成熟期(次年9—10月)。本文選取不同蘋(píng)果生育期內(nèi)的9個(gè)氣象特征變量、3個(gè)空間特征變量和3個(gè)氣象災(zāi)害特征變量作為構(gòu)建蘋(píng)果相對(duì)氣象產(chǎn)量預(yù)測(cè)模型的輸入變量,具體特征變量如表2所示。

表2 蘋(píng)果相對(duì)氣象產(chǎn)量模擬模型輸入特征變量

太陽(yáng)輻射數(shù)據(jù)根據(jù)各氣象站點(diǎn)觀測(cè)的日照時(shí)數(shù)和Angstrom-Prescott經(jīng)驗(yàn)公式[21]進(jìn)行估算,公式為

(2)

其中

Ra=(24×60/π)Gscdr(ωssinφsinδ+

cosφcosδ)

(3)

dr=1+0.033cos(2πJ/365)

(4)

δ=0.409sin(2πJ/365-1.39)

(5)

ωs=arccos(-tanφtanδ)

(6)

N=24ωs/π

(7)

式中N——最大可能日照時(shí)數(shù),h

ωs——太陽(yáng)時(shí)角,rad

φ——緯度,rad

δ——太陽(yáng)磁偏角,rad

J——年內(nèi)某天的日序數(shù)

dr——日地間相對(duì)距離的倒數(shù)

Gsc——太陽(yáng)常數(shù),取0.082 MJ/(m2·min)

Ra——大氣外總輻射,MJ/(m2·min)

as、bs——經(jīng)驗(yàn)系數(shù),取0.19和0.54[22]

n——實(shí)際日照時(shí)數(shù),h

Rs——地表總輻射,MJ/(m2·min)

1個(gè)月尺度標(biāo)準(zhǔn)化降水指數(shù)(SPEI)根據(jù)Thornthwaite法[23]通過(guò)R語(yǔ)言中的“SPEI”包計(jì)算得到,具體計(jì)算過(guò)程主要分為3個(gè)步驟。

(1)計(jì)算潛在蒸散量PET,計(jì)算式為

(8)

(9)

A=6.75×10-7H3-7.71×10-5H2+

1.792×10-2H+0.49

(10)

式中A——常數(shù)

H——年熱量指數(shù)

PETi——月潛在蒸散量,mm/d

(2)計(jì)算氣候水平衡,計(jì)算式為

Di=Pi-PETi

(11)

式中Di——月降水量與蒸散量差值,mm/d

Pi——月降水量,mm/d

(3)對(duì)Di數(shù)據(jù)序列標(biāo)準(zhǔn)化,采用三參數(shù)的log-logistic概率分布F(x)對(duì)其進(jìn)行擬合,計(jì)算出每個(gè)Di對(duì)應(yīng)的SPEI值,計(jì)算式為

(12)

(13)

(14)

式中SPEI——1個(gè)月尺度標(biāo)準(zhǔn)化降水指數(shù)

其中C0、C1、C2、d1、d2、d3為常數(shù)項(xiàng),取2.515 517、0.802 853、0.010 328、1.432 788、0.189 269、0.001 308。

為消除氣象變量多年周期性循環(huán)的影響,本文將蘋(píng)果多年生長(zhǎng)季內(nèi)的氣象特征變量進(jìn)行平均,再將其與目標(biāo)變量(蘋(píng)果相對(duì)氣象產(chǎn)量)進(jìn)行斯皮爾曼相關(guān)性分析(Spearman correlation analysis),選擇對(duì)蘋(píng)果產(chǎn)量影響最大的重要?dú)庀筇卣髯兞?以減少輸入數(shù)據(jù)的復(fù)雜程度。在每類(lèi)氣象特征變量中,選擇與目標(biāo)變量的相關(guān)系數(shù)絕對(duì)值最大且顯著的特征變量。如果同類(lèi)特征變量之間的相關(guān)系數(shù)小于0.5且其與目標(biāo)變量顯著相關(guān)時(shí),則可選擇多個(gè)特征變量。

1.3.2模型目標(biāo)變量構(gòu)建

蘋(píng)果產(chǎn)量的形成受多種因素的影響,主要包括氣象要素、技術(shù)措施和其他因素,相應(yīng)地蘋(píng)果產(chǎn)量可以分解為氣象產(chǎn)量、趨勢(shì)產(chǎn)量和隨機(jī)產(chǎn)量。為了能更好地反映氣象要素造成的蘋(píng)果產(chǎn)量波動(dòng),降低時(shí)間和地域?qū)Ξa(chǎn)量影響的限制,采用相對(duì)氣象產(chǎn)量作為目標(biāo)變量定量評(píng)估氣象要素對(duì)產(chǎn)量形成的影響。計(jì)算公式為

(15)

其中

(16)

yc=y-yt+yr

(17)

式中y——實(shí)際產(chǎn)量,kg/hm2

yw——相對(duì)氣象產(chǎn)量,正值表示氣象要素有利于作物生長(zhǎng)發(fā)育即產(chǎn)量增加,負(fù)值表示產(chǎn)量減少

yt——趨勢(shì)產(chǎn)量,kg/hm2,可采用一階傅里葉擬合來(lái)估算[24]

yc——?dú)庀螽a(chǎn)量,kg/hm2

yr——隨機(jī)產(chǎn)量,kg/hm2,本文忽略不計(jì)

a0——常系數(shù),取-2.94×1010～2.17×1011

w——常系數(shù),取0.30～0.90

ai——各級(jí)系數(shù),取-2.17×1011～2.94×1010

bi——各級(jí)系數(shù),取-2.09×104～1.22

M——傅里葉展開(kāi)級(jí)數(shù),取1

t——年份

1.4 蘋(píng)果相對(duì)氣象產(chǎn)量模擬模型

1.4.1機(jī)器學(xué)習(xí)算法

(1)梯度提升樹(shù)

梯度提升樹(shù)(GBDT)是基于Booting方法的集成模型。該類(lèi)方法采用加法模型和前向分步算法實(shí)現(xiàn)學(xué)習(xí)優(yōu)化[25]。GBDT回歸算法中以回歸樹(shù)作為弱學(xué)習(xí)器,學(xué)習(xí)任務(wù)時(shí),利用損失函數(shù)的負(fù)梯度作為當(dāng)前學(xué)習(xí)器的偽殘差,根據(jù)偽殘差擬合回歸樹(shù),直到滿足要求停止重復(fù)此過(guò)程。GBDT算法一般比隨機(jī)森林算法的訓(xùn)練結(jié)果更準(zhǔn)確[9],收斂速度快且不易出現(xiàn)過(guò)擬合。R語(yǔ)言中提供了GBDT算法包“gbm”,所需要的參數(shù)包括:分布函數(shù)(distribution,選擇“gaussian”)、樹(shù)終節(jié)點(diǎn)的最小個(gè)數(shù)(n.minobsinnode,取10)、學(xué)習(xí)率(shrinkage,取0.01)、回歸樹(shù)數(shù)量(n.trees,取200)以及單棵回歸樹(shù)最大深度(interaction.depth,以1為步長(zhǎng)在[0,10]內(nèi)取值)。

(2)支持向量機(jī)

支持向量機(jī)(SVM)通過(guò)非線性映射將自變量映射到高維的特征空間,在高維特征空間中尋找一個(gè)最優(yōu)超平面,使得所有訓(xùn)練樣本距離該最優(yōu)分類(lèi)面誤差最小[26]。對(duì)于回歸應(yīng)用,SVM算法主要利用核函數(shù)構(gòu)造線性回歸方程求得最優(yōu)超平面[27]。該算法可利用R語(yǔ)言提供的“e1071”包來(lái)實(shí)現(xiàn),所需參數(shù)包括:核函數(shù)類(lèi)型(Kernel,定為高斯徑向基核函數(shù))、核函數(shù)系數(shù)(gamma,以0.01為步長(zhǎng)在[0.01,1.00]內(nèi)選擇)以及懲罰系數(shù)(cost,以0.1為步長(zhǎng)在[0.1,100]內(nèi)確定)。

(3)貝葉斯正則化反向傳播神經(jīng)網(wǎng)絡(luò)

貝葉斯正則化反向傳播神經(jīng)網(wǎng)絡(luò)(BRBP)融合傳統(tǒng)反向傳播神經(jīng)網(wǎng)絡(luò)算法和貝葉斯理論自動(dòng)選擇最優(yōu)的正則化參數(shù),能夠避免過(guò)擬合問(wèn)題,極大提高模型的泛化能力[28]。采用包含輸入層、隱含層和輸出層的3層神經(jīng)網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu),隱含層的節(jié)點(diǎn)數(shù)根據(jù)試湊法,以預(yù)測(cè)誤差最小為原則來(lái)確定。Matlab軟件中的“newff”函數(shù)可以實(shí)現(xiàn)該類(lèi)算法,所需參數(shù)包括:輸入節(jié)點(diǎn)數(shù)a、輸出節(jié)點(diǎn)數(shù)b、隱含層節(jié)點(diǎn)數(shù)l、節(jié)點(diǎn)傳遞函數(shù)(隱含層為tanh,輸出層為purelin)、訓(xùn)練函數(shù)(trainbr)、迭代次數(shù)k、學(xué)習(xí)率η以及目標(biāo)函數(shù)誤差(goal)。其中,參數(shù)k、η、goal經(jīng)反復(fù)試算后分別取1 000、0.1、10-5。

(4)多元線性回歸

多元線性回歸(MLR)是一種廣泛應(yīng)用于單個(gè)目標(biāo)變量與2個(gè)及以上輸入特征變量的線性回歸方法[29]。目標(biāo)變量通常受多種特征變量的影響,因此多元線性回歸相對(duì)于一元線性回歸更具適用性。在多元線性回歸中,方程的參數(shù)通過(guò)最小二乘法獲得[30]。R語(yǔ)言提供的“l(fā)m”函數(shù)可實(shí)現(xiàn)MLR算法。

1.4.2蘋(píng)果相對(duì)氣象產(chǎn)量模擬模型構(gòu)建

如圖2(圖中M、G、D分別表示最重要?dú)庀鬄?zāi)害特征變量、空間特征變量、氣象災(zāi)害特征變量)所示,以GBDT算法為例,簡(jiǎn)述蘋(píng)果相對(duì)氣象產(chǎn)量yw預(yù)測(cè)模型的建立過(guò)程。首先,基于蘋(píng)果全生育期最重要?dú)庀筇卣髯兞?逐次添加空間特征變量和氣象災(zāi)害特征變量,形成不同的模型輸入特征變量組合。其次,應(yīng)用GBDT算法建立模型估算蘋(píng)果相對(duì)氣象產(chǎn)量,以yw估算值和實(shí)測(cè)值之間的誤差最小化為優(yōu)化目標(biāo),獲得最佳模型輸入特征變量組合。然后,逐次添加蘋(píng)果落葉休眠期、萌芽幼果期、果實(shí)膨大期、著色成熟期內(nèi)的最佳模型輸入變量組合,構(gòu)建不同生育期的蘋(píng)果相對(duì)氣象產(chǎn)量預(yù)測(cè)模型,以模型測(cè)試集均方根誤差(RMSE)最小化為依據(jù),確定蘋(píng)果相對(duì)氣象產(chǎn)量預(yù)測(cè)的提前生育期。最后,逐次添加蘋(píng)果生長(zhǎng)季內(nèi)各月份的最佳模型輸入特征變量組合,構(gòu)建蘋(píng)果生長(zhǎng)季內(nèi)不同月份的蘋(píng)果相對(duì)氣象產(chǎn)量預(yù)測(cè)模型,以模型測(cè)試集歸一化RMSE最小化為依據(jù),確定蘋(píng)果相對(duì)氣象產(chǎn)量預(yù)測(cè)的提前月份。并且在不同生育期和生長(zhǎng)季內(nèi)不同月份的模型預(yù)測(cè)精度分析過(guò)程中,將相對(duì)氣象產(chǎn)量轉(zhuǎn)換為實(shí)際產(chǎn)量,比較其與以相對(duì)氣象產(chǎn)量為目標(biāo)變量的機(jī)器學(xué)習(xí)模型精度,從而判斷蘋(píng)果相對(duì)氣象產(chǎn)量模擬模型構(gòu)建方法精度。

圖2 蘋(píng)果相對(duì)氣象產(chǎn)量模擬模型流程圖

模型模擬前,隨機(jī)選取70%觀測(cè)數(shù)據(jù)作為訓(xùn)練數(shù)據(jù)集,30%數(shù)據(jù)作為測(cè)試數(shù)據(jù)集,以10折交叉驗(yàn)證誤差為原則對(duì)模型參數(shù)進(jìn)行尋優(yōu)。為降低模型精度極大值和極小值的影響,模型的最終精度為100次運(yùn)行誤差的50%分位數(shù)。

1.5 模型精度評(píng)價(jià)

采用皮爾遜相關(guān)系數(shù)(Pearson correlation coefficient,r)和均方根誤差(Root mean square error,RMSE)2個(gè)統(tǒng)計(jì)指標(biāo)來(lái)評(píng)價(jià)模型模擬精度。

2 結(jié)果與分析

2.1 影響蘋(píng)果產(chǎn)量的氣象特征變量

本文分析了目標(biāo)變量(蘋(píng)果相對(duì)氣象產(chǎn)量)與氣象特征變量之間的相關(guān)關(guān)系(圖3,P<0.05),發(fā)現(xiàn)氣溫特征變量Tmean、Tmin、Tmax、Td、AT、NT與目標(biāo)變量之間存在正相關(guān)關(guān)系,其中Tmax與目標(biāo)變量的相關(guān)性最高(r=0.26)。為了避免輸入因子間共線性對(duì)模型精度的影響,可同時(shí)選擇與Tmax相關(guān)系數(shù)小于0.5的其他氣溫特征變量,其中符合要求的兩個(gè)氣溫特征變量為T(mén)min與NT,但是它們之間存在較強(qiáng)的相關(guān)性,且Tmin與目標(biāo)變量之間的相關(guān)性強(qiáng)于NT,因此最終只選擇Tmax和Tmin作為模型輸入的氣溫特征變量。需水特征變量RH與目標(biāo)變量之間呈負(fù)相關(guān)關(guān)系,而供水特征變量P、輻射特征變量R與目標(biāo)變量之間的相關(guān)性系數(shù)較小。因此,本文最終選取Tmax、Tmin、RH、P、R作為蘋(píng)果相對(duì)氣象產(chǎn)量模擬模型輸入的最重要?dú)庀筇卣髯兞俊?/p>

圖3 目標(biāo)變量(蘋(píng)果相對(duì)氣象產(chǎn)量)與氣象特征變量關(guān)聯(lián)熱圖

2.2 基于多種特征變量的不同蘋(píng)果產(chǎn)量模型模擬精度分析

分別使用蘋(píng)果整個(gè)生長(zhǎng)季內(nèi)的最重要?dú)庀筇卣髯兞拷M合(M),最重要?dú)庀筇卣髯兞亢涂臻g特征變量組合(M+G),以及最重要?dú)庀筇卣髯兞俊⒖臻g特征變量和氣象災(zāi)害特征變量組合(M+G+D)作為不同的模型輸入特征變量組合,分別驅(qū)動(dòng)基于GBDT、SVM、BRBP和MLR算法的蘋(píng)果相對(duì)氣象產(chǎn)量模擬模型,并比較不同模型測(cè)試數(shù)據(jù)集的模擬精度(圖4,圖中虛線表示1∶1線;**表示P<0.01)。當(dāng)僅以最重要?dú)庀筇卣髯兞縈作為模型輸入變量特征組合時(shí),SVM模型具有最大的r值(0.70),而B(niǎo)RBP模型的預(yù)測(cè)誤差最小(RMSE為0.46)且數(shù)據(jù)點(diǎn)均勻分布在1∶1線的兩側(cè)(圖4a、4d、4g、4j)。當(dāng)以M+G作為模型輸入特征變量集時(shí),與僅以M作為模型輸入特征變量的情形相比,GBDT和BRBP的RMSE分別降低2%、2.17%,同時(shí)GBDT、BRBP和MLR模型的r分別提升4.76%、2.99%和1.72%,但SVM模型的精度無(wú)變化(圖4b、4e、4h、4k)。可見(jiàn)在模型輸入氣象特征變量中再加入空間特征變量,可以提高蘋(píng)果相對(duì)氣象產(chǎn)量模擬模型的預(yù)測(cè)精度。進(jìn)一步在上述模型輸入變量中再添加氣象災(zāi)害特征變量(M+G+D),結(jié)果表明增加氣象災(zāi)害特征變量可使GBDT、BRBP和MLR模型預(yù)測(cè)精度進(jìn)一步提升,其中GBDT模型的模擬精度提升效果最為明顯,其r增加16.67%,RMSE減少10.20%,而SVM模型的模擬精度則有所下降(圖4c、4f、4i、4l)。

總體而言,基于GBDT和BRBP這2種機(jī)器學(xué)習(xí)算法和M+G+D輸入特征變量組合建立的蘋(píng)果相對(duì)氣象產(chǎn)量預(yù)測(cè)模型,其模擬精度普遍優(yōu)于MLR算法模擬。這可能是因?yàn)樘O(píng)果相對(duì)氣象產(chǎn)量與最重要?dú)庀筇卣髯兞?、空間特征變量和氣象災(zāi)害特征變量之間存在較強(qiáng)的非線性關(guān)系,而機(jī)器學(xué)習(xí)算法比多元線性回歸方法更能準(zhǔn)確描述這種非線性關(guān)系。因此,本文采用M+G+D組合作為最佳的模型輸入特征變量組合。

2.3 模型模擬誤差空間分析

以蘋(píng)果整個(gè)生長(zhǎng)季內(nèi)的M+G+D特征變量組合作為模型輸入變量,驅(qū)動(dòng)基于不同算法建立的蘋(píng)果相對(duì)氣象產(chǎn)量模擬模型,可進(jìn)一步分析不同模型在黃土高原各省份的總體模擬精度(圖5,圖中每個(gè)箱體包含各省蘋(píng)果生產(chǎn)基地縣的蘋(píng)果相對(duì)氣象產(chǎn)量模擬結(jié)果,箱體上、下限為數(shù)據(jù)的上四分位數(shù)和下四分位數(shù),中間實(shí)線為數(shù)據(jù)中位數(shù),箱體外的兩條線為上邊緣和下邊緣,超出上下邊緣的值為異常值)。由圖5可知,不同模型在不同省份的表現(xiàn)存在較大差異,其中BRBP模型在河南省和陜西省表現(xiàn)最佳,SVM模型在山西省和甘肅省表現(xiàn)較佳。河南、山西、陜西和甘肅各省最佳模型的 RMSE均值分別為0.47、0.11、0.39和0.20。在這4個(gè)省份中,山西省的模型總體估計(jì)誤差相對(duì)較小,4個(gè)模型的RMSE平均值僅為0.24,且誤差分布相對(duì)較為集中。但河南省RMSE箱型圖的箱體最長(zhǎng),這可能由于本文中河南省僅有3個(gè)蘋(píng)果基地縣,數(shù)據(jù)量較少使得箱型圖的分布范圍較大。

圖5 基于最重要?dú)庀?、空間和氣象災(zāi)害特征變量組合和不同算法的模擬模型在4個(gè)不同省份模擬蘋(píng)果相對(duì)氣象產(chǎn)量的RMSE箱型圖

2.4 蘋(píng)果不同生育期輸入特征變量組合對(duì)模型模擬精度的影響

進(jìn)一步探究在M+G+D輸入變量組合下,不同模型在不同生育期(落葉休眠期、萌芽幼果期、果實(shí)膨大期、著色成熟期)的蘋(píng)果產(chǎn)量預(yù)測(cè)精度(圖6)。當(dāng)以相對(duì)氣象產(chǎn)量作為目標(biāo)變量時(shí),對(duì)于GBDT和BRBP模型,各生育期模型的模擬精度無(wú)明顯變化,其RMSE分別保持在0.44～0.46和0.44～0.45之間,這表明這2種模型在蘋(píng)果各個(gè)生育期內(nèi)均能獲得相對(duì)較高的相對(duì)氣象產(chǎn)量預(yù)測(cè)精度(圖6a)。在不同蘋(píng)果生育期,SVM模型的RMSE先減小后增加,次年8月產(chǎn)量預(yù)測(cè)的RMSE具有最低值(0.51),可見(jiàn)SVM模型可在果實(shí)膨大期獲取最好的預(yù)測(cè)結(jié)果。MLR模型模擬的RMSE隨著蘋(píng)果生長(zhǎng)發(fā)育的推進(jìn)而逐漸減小,并在次年8月獲得最佳的預(yù)測(cè)結(jié)果,與次年10月該模型的蘋(píng)果相對(duì)氣象產(chǎn)量預(yù)測(cè)精度僅相差0.01,表明該模型可在蘋(píng)果果實(shí)膨大期就獲取較為理想的蘋(píng)果相對(duì)氣象產(chǎn)量預(yù)測(cè)結(jié)果。當(dāng)將相對(duì)氣象產(chǎn)量轉(zhuǎn)換為實(shí)際產(chǎn)量時(shí),GBDT、BRBP、MLR模型在各生育期的模擬精度均呈現(xiàn)減小趨勢(shì)(圖6b)。具體而言,GBDT和BRBP模型在蘋(píng)果各個(gè)生育期內(nèi)均能獲得相對(duì)較高的實(shí)際產(chǎn)量預(yù)測(cè)精度,其RMSE分別保持在4.42～4.90 kg/hm2和4.62～5.20 kg/hm2之間。MLR模型在果實(shí)膨大期的預(yù)測(cè)精度與著色成熟期僅相差0.03 kg/hm2,表明該模型在果實(shí)膨大期也能獲取較為理想的蘋(píng)果產(chǎn)量預(yù)測(cè)精度。SVM模型的RMSE先減小后增加,并且在果實(shí)膨大期獲取最好的預(yù)測(cè)結(jié)果(5.04 kg/hm2)。

圖6 基于GBDT、SVM、BRBP和MLR算法的不同模擬模型在蘋(píng)果不同生育期模擬黃土高原蘋(píng)果產(chǎn)量的精度對(duì)比

2.5 蘋(píng)果生長(zhǎng)季內(nèi)不同月份的模型預(yù)測(cè)精度

在生長(zhǎng)季內(nèi)不同月份對(duì)蘋(píng)果相對(duì)氣象產(chǎn)量進(jìn)行預(yù)測(cè)模擬并分析其精度,可確定出蘋(píng)果產(chǎn)量預(yù)測(cè)的最佳提前期。根據(jù)生長(zhǎng)季內(nèi)不同月份各模型模擬結(jié)果歸一化RMSE,可知4種模型的模擬精度都呈現(xiàn)出類(lèi)似的變化趨勢(shì)(圖7a)。即隨著蘋(píng)果生育期的推進(jìn),模型輸入變量不斷增多,各模型的歸一化RMSE逐漸減少。當(dāng)以相對(duì)氣象產(chǎn)量作為目標(biāo)變量時(shí),對(duì)于GBDT和BRBP算法,基于生長(zhǎng)季內(nèi)所有月份輸入特征變量的結(jié)果最佳。與基于生長(zhǎng)季所有月份輸入變量的模型模擬結(jié)果歸一化RMSE相比,在次年9月這2個(gè)模型的歸一化RMSE分別與之僅相差0.62%和1.04%,精度非常相近。SVM算法可在次年8月獲取最佳的預(yù)測(cè)精度?；贛LR算法的蘋(píng)果相對(duì)氣象產(chǎn)量預(yù)測(cè)模型可在次年9月獲得最小的歸一化RMSE。當(dāng)將相對(duì)氣象產(chǎn)量轉(zhuǎn)換為實(shí)際產(chǎn)量時(shí),GBDT和BRBP算法在次年9月的產(chǎn)量預(yù)測(cè)精度分別與次年10月僅相差0.38%和1.34%,而SVM和MLR算法分別在次年8月和次年9月模型預(yù)測(cè)精度最佳(圖7b)。綜上所述,本文建立的4種不同模型均可在蘋(píng)果成熟期前1～2個(gè)月實(shí)現(xiàn)蘋(píng)果相對(duì)氣象產(chǎn)量的相對(duì)準(zhǔn)確預(yù)測(cè)。

圖7 基于GBDT、SVM、BRBP和MLR算法的不同模擬模型在蘋(píng)果生長(zhǎng)季內(nèi)各月份模擬黃土高原蘋(píng)果產(chǎn)量的歸一化RMSE

3 討論

3.1 影響蘋(píng)果產(chǎn)量的特征變量和產(chǎn)量去趨勢(shì)方法

氣象要素和氣象災(zāi)害對(duì)蘋(píng)果產(chǎn)量具有重要的影響。本文選擇了蘋(píng)果產(chǎn)量形成過(guò)程中具有代表性的氣象變量和氣象災(zāi)害變量,包括Tmean、Tmin、Tmax、Td、AT、NT、RH、P、R、N4≤0℃、CR和SPEI。其次,氣象特征變量篩選結(jié)果表明Tmin、Tmax、RH、P、R是影響黃土高原蘋(píng)果相對(duì)氣象產(chǎn)量預(yù)測(cè)模型的精度的重要?dú)庀筇卣髯兞?這與前人的研究結(jié)果相一致。白秀廣等[31]和張強(qiáng)[32]研究表明,氣溫對(duì)蘋(píng)果產(chǎn)量有正向影響,降水量、空氣相對(duì)濕度和日照時(shí)數(shù)對(duì)蘋(píng)果產(chǎn)量有負(fù)面影響。本文中降水量和太陽(yáng)輻射與蘋(píng)果相對(duì)氣象產(chǎn)量之間的相關(guān)系數(shù)較小,這可能是因?yàn)辄S土高原降水不均勻特征和灌溉等因素影響了降水量與蘋(píng)果產(chǎn)量間的相關(guān)性[33]。此外,本文中蘋(píng)果產(chǎn)量的樣本量可能也影響了太陽(yáng)輻射與蘋(píng)果產(chǎn)量之間的相關(guān)性。隨后,在最重要?dú)庀筇卣髯兞恐幸来渭尤肟臻g特征變量和氣象災(zāi)害特征變量,均可進(jìn)一步提升基于GBDT、BRBP和MLR算法的蘋(píng)果相對(duì)氣象產(chǎn)量模擬模型的預(yù)測(cè)精度。當(dāng)模型輸入特征組合為“氣象特征變量M+空間特征變量G+氣象災(zāi)害特征變量D”時(shí)模型精度最高。劉峻明等[34]研究表明,氣象站點(diǎn)空間位置會(huì)影響作物的積溫和物候,但由于氣候條件往往存在一定的年際變化,導(dǎo)致空間位置和氣象要素的相關(guān)性較弱。所以,作物產(chǎn)量預(yù)測(cè)模型中需要同時(shí)考慮氣象特征變量和空間特征變量。FENG等[35]研究發(fā)現(xiàn),需要在作物產(chǎn)量預(yù)測(cè)模型中加入氣象災(zāi)害,以此減少氣象災(zāi)害給作物產(chǎn)量預(yù)測(cè)結(jié)果帶來(lái)的不確定性。然而,本文中加入氣象災(zāi)害特征變量會(huì)導(dǎo)致SVM模型的預(yù)測(cè)精度降低,這可能是因?yàn)樵撍惴ㄒ资軘?shù)據(jù)維度的影響,而氣象災(zāi)害因子的加入增加了樣本維度,從而導(dǎo)致模型性能受損[36]。

蘋(píng)果相對(duì)氣象產(chǎn)量?jī)H受氣象條件的影響,但是該產(chǎn)量不能通過(guò)測(cè)量直接獲取,因此需對(duì)實(shí)際蘋(píng)果產(chǎn)量數(shù)據(jù)進(jìn)行預(yù)處理,而數(shù)據(jù)處理方法的選擇就甚為重要。本文采用了一階傅里葉對(duì)所有基地縣的蘋(píng)果單產(chǎn)去趨勢(shì)得到蘋(píng)果相對(duì)氣象產(chǎn)量,這可以在一定程度上消除其他因素對(duì)產(chǎn)量的影響,從而能準(zhǔn)確地反映蘋(píng)果產(chǎn)量隨氣象因子的非線性變化。然而,產(chǎn)量去趨勢(shì)其他方法還包括:線性回歸、多項(xiàng)式回歸、指數(shù)回歸、對(duì)數(shù)回歸、滑動(dòng)平均、Logistic回歸、HP濾波、局部加權(quán)回歸等[37-40]。此外,由于社會(huì)、經(jīng)濟(jì)發(fā)展不平衡,不同區(qū)域的作物產(chǎn)量變化趨勢(shì)也會(huì)存在差異性[41]。因此,在今后的研究中需要在各蘋(píng)果生產(chǎn)縣采用不同的去趨勢(shì)方法獲得更加準(zhǔn)確、客觀的蘋(píng)果相對(duì)氣象產(chǎn)量,從而提升模型的模擬精度。

3.2 基于機(jī)器學(xué)習(xí)算法的蘋(píng)果相對(duì)氣象產(chǎn)量模擬模型

機(jī)器學(xué)習(xí)模型能夠較好地建立氣象災(zāi)害因子與產(chǎn)量之間的關(guān)系,而這是作物模型普遍無(wú)法實(shí)現(xiàn)的[42]。在模擬黃土高原蘋(píng)果相對(duì)氣象產(chǎn)量的過(guò)程中,本文發(fā)現(xiàn)非線性算法GBDT、BRBP優(yōu)于線性算法MLR。當(dāng)以M+G+D組合作為模型輸入特征變量組合來(lái)驅(qū)動(dòng)本文所建立的模型模擬蘋(píng)果相對(duì)氣象產(chǎn)量時(shí),機(jī)器學(xué)習(xí)模型(GBDT、BRBP模型)的r和RMSE分別比線性模型(MLR模型)高11.11%～22.22%和低10.20%,這表明模型輸入特征變量與蘋(píng)果相對(duì)氣象產(chǎn)量之間存在較強(qiáng)的非線性關(guān)系,機(jī)器學(xué)習(xí)算法比線性算法能夠更為準(zhǔn)確地描述這種非線性關(guān)系。前人研究也發(fā)現(xiàn)氣溫、太陽(yáng)輻射、降水量、空氣相對(duì)濕度等氣象要素、空間要素和氣象災(zāi)害要素與蘋(píng)果產(chǎn)量之間存在非線性關(guān)系[11,24,34,43-44],本文結(jié)果與這些結(jié)論基本相同。然而,基于SVM算法建立的蘋(píng)果相對(duì)氣象產(chǎn)量預(yù)測(cè)模型的RMSE均高于線性模型,這可能是因?yàn)檩斎胱兞恐g存在較強(qiáng)的線性關(guān)系,破壞了SVM模型的預(yù)測(cè)性能,從而導(dǎo)致SVM模型預(yù)測(cè)結(jié)果與線性模型相近甚至更差[45]。此外,本文還探究了在蘋(píng)果生長(zhǎng)季內(nèi)不同月份各模型預(yù)測(cè)蘋(píng)果相對(duì)氣象產(chǎn)量的精度。結(jié)果表明,當(dāng)以M+G+D特征變量組合作為模型輸入變量時(shí),本文建立的4種模型均可在果實(shí)成熟前的1～2個(gè)月獲得較高精度的蘋(píng)果相對(duì)氣象產(chǎn)量預(yù)測(cè)結(jié)果。本文所構(gòu)建的蘋(píng)果相對(duì)氣象產(chǎn)量模擬模型在蘋(píng)果生長(zhǎng)季內(nèi)不同月份的預(yù)測(cè)中模擬精度存在差異,這一方面可能因?yàn)椴煌路輾庀筇卣髯兞繉?duì)蘋(píng)果產(chǎn)量的影響并非同等重要[13]。另一方面可能是由于特征變量的輸入個(gè)數(shù)隨著月份的遞增而增加,機(jī)器學(xué)習(xí)模型對(duì)更多數(shù)量的輸入特征通常模擬精度可能更佳。已有研究表明小麥、玉米、水稻等糧食作物均可在收獲前的1～2個(gè)月實(shí)現(xiàn)較為準(zhǔn)確的早期產(chǎn)量預(yù)測(cè),這些結(jié)果也與本文的結(jié)果較為一致[21,46]。并且將相對(duì)氣象產(chǎn)量轉(zhuǎn)換為實(shí)際產(chǎn)量后,4種模型在蘋(píng)果生長(zhǎng)季內(nèi)不同月份的模擬精度變化趨勢(shì)基本一致,僅在數(shù)值上存在部分差異,可能是蘋(píng)果實(shí)際產(chǎn)量受到技術(shù)措施等其他因素影響所導(dǎo)致的。這表明本文所構(gòu)建的蘋(píng)果產(chǎn)量模型具有較好的預(yù)測(cè)精度。

3.3 氣象特征變量和氣象災(zāi)害特征變量對(duì)蘋(píng)果相對(duì)氣象產(chǎn)量模擬精度的影響

本文通過(guò)篩選,最終確定了預(yù)測(cè)蘋(píng)果相對(duì)氣象產(chǎn)量的輸入特征變量為:氣象特征變量Tmax、Tmin、RH、P、R,氣象災(zāi)害特征變量N4≤0℃、CR、SPEI,空間特征變量Lat、Lon、Ele。以機(jī)器學(xué)習(xí)模型(GBDT、SVM、BRBP、MLR模型)在蘋(píng)果生長(zhǎng)季內(nèi)不同月份的表現(xiàn)來(lái)看,空間特征變量在模型輸入中是固定不變的,但是隨著蘋(píng)果生育期從冬季休眠期(本年11月)發(fā)展到收獲期(次年10月),輸入模型的最重要?dú)庀筇卣髯兞亢蜌庀鬄?zāi)害特征變量的時(shí)間序列逐漸遞增,而此過(guò)程中所有模型的RMSE均逐漸減小(圖7)。這表明氣象特征變量和氣象災(zāi)害特征變量在模型性能的后續(xù)提升過(guò)程中發(fā)揮著關(guān)鍵作用。因此,相比于空間特征變量,氣象特征變量和氣象災(zāi)害特征變量對(duì)蘋(píng)果相對(duì)氣象產(chǎn)量模擬模型精度的影響更大。

3.4 蘋(píng)果相對(duì)氣象產(chǎn)量模擬模型的不足

本文建立的黃土高原蘋(píng)果相對(duì)氣象產(chǎn)量預(yù)測(cè)模型也存在一定的不足。首先,本文所建立模型的輸入數(shù)據(jù)較為單一,然而蘋(píng)果的生長(zhǎng)發(fā)育過(guò)程非常復(fù)雜,產(chǎn)量形成過(guò)程中受到眾多因素的影響,例如土壤[40,47]、灌溉、施肥以及園藝管理措施[48]等。因此,在未來(lái)的研究中應(yīng)考慮將上述因子納入模型的輸入特征變量組合,這將會(huì)進(jìn)一步提高蘋(píng)果相對(duì)氣象產(chǎn)量模型的模擬精度。其次,本文選用的3種機(jī)器學(xué)習(xí)算法均為黑箱模型,并不能定量化解釋輸入變量對(duì)目標(biāo)變量的影響[49-50]。最后,本文還未能將所選的最佳預(yù)測(cè)模型在其他不同蘋(píng)果種植區(qū)域和不同的氣候環(huán)境下進(jìn)行充分驗(yàn)證。例如,LI等[11]建立了陜西省縣域產(chǎn)量預(yù)測(cè)模型,并預(yù)測(cè)了未來(lái)氣候條件下28個(gè)蘋(píng)果基地縣的產(chǎn)量。因此,在隨后的研究中擬將本文建立的蘋(píng)果相對(duì)氣象產(chǎn)量預(yù)測(cè)模型推廣至整個(gè)黃土高原的蘋(píng)果種植適宜區(qū),并利用其分析未來(lái)氣候變化條件下黃土高原產(chǎn)區(qū)的蘋(píng)果相對(duì)氣象產(chǎn)量變化。

4 結(jié)論

(1)基于GBDT和BRBP機(jī)器學(xué)習(xí)算法的蘋(píng)果相對(duì)氣象產(chǎn)量模擬模型的精度均優(yōu)于基于傳統(tǒng)線性回歸算法的模型,這表明最重要?dú)庀筇卣髯兞?、空間特征變量和氣象災(zāi)害特征變量與蘋(píng)果產(chǎn)量之間存在較強(qiáng)的非線性關(guān)系,而機(jī)器學(xué)習(xí)算法能夠更為準(zhǔn)確地描述這種非線性關(guān)系。

(2)以最重要?dú)庀筇卣髯兞俊⒖臻g特征變量和氣象災(zāi)害特征變量作為模型輸入變量組合,GBDT和BRBP算法模型在各個(gè)生育期內(nèi)均能獲得相對(duì)較高的蘋(píng)果相對(duì)氣象產(chǎn)量預(yù)測(cè)精度,SVM和MLR算法模型可在果實(shí)膨大期獲取較為理想的蘋(píng)果相對(duì)氣象產(chǎn)量模擬結(jié)果。

(3)基于最重要?dú)庀筇卣髯兞?、空間特征變量和氣象災(zāi)害特征變量組合,機(jī)器學(xué)習(xí)算法(GBDT、SVM、BRBP)和線性算法(MLR)的產(chǎn)量模擬模型均可在蘋(píng)果成熟期前1～2個(gè)月實(shí)現(xiàn)對(duì)黃土高原蘋(píng)果相對(duì)氣象產(chǎn)量較為準(zhǔn)確的預(yù)測(cè)。