任亞飛， 田 帥， 邵馨葉,2， 邵建龍*

(1.昆明理工大學(xué) 信息工程與自動(dòng)化學(xué)院，云南 昆明 650500；2.云南電網(wǎng)有限責(zé)任公司 信息中心，云南 昆明 650200)

參考作物騰發(fā)量(ET0)是計(jì)算蒸散量和估算作物實(shí)際需水量的重要參數(shù)，在灌溉制度的選擇、作物種植結(jié)構(gòu)優(yōu)化、田間水分平衡理論等方面有著廣泛的應(yīng)用。ET0的計(jì)算受到多個(gè)氣象參數(shù)的影響，包括溫度、相對濕度、風(fēng)速、日照時(shí)長、氣壓和太陽輻射量等，因此很難建立一個(gè)簡單的方程來獲得不同氣候條件下的準(zhǔn)確估計(jì)值[1-3]。目前基于氣候變量估算ET0的方法大致可以分為三類：輻射法、溫度法和合成法[4]。但是上述方法無法完全體現(xiàn)出ET0計(jì)算過程中非線性、動(dòng)態(tài)性和高度復(fù)雜性的特點(diǎn)，導(dǎo)致其計(jì)算結(jié)果在不同區(qū)域存在截然不同的精度，適用性較差。1998年由聯(lián)合國糧農(nóng)組織提出的經(jīng)典Peman-Monteith公式被普遍認(rèn)為是目前計(jì)算ET0的唯一精確方法[5]。然而，PM模型需要大量的氣象數(shù)據(jù)作為輸入，計(jì)算過程涉及到這些因素之間復(fù)雜的非線性計(jì)算。因此，為了使其更加方便的應(yīng)用于農(nóng)業(yè)生產(chǎn)指導(dǎo)中，有必要需要尋找一個(gè)更加簡單、精確的模型來計(jì)算ET0。

近年來，隨著神經(jīng)網(wǎng)絡(luò)理論的迅速發(fā)展，基于此的BP神經(jīng)網(wǎng)絡(luò)模型在ET0預(yù)測中得到了廣泛的應(yīng)用。柳燁等[6]建立了基于溫濕度的ET0估算模型，通過對陜西省6個(gè)站點(diǎn)的數(shù)據(jù)進(jìn)行分析，證明該模型不僅能得到較為精確的預(yù)測結(jié)果且具有更強(qiáng)的魯棒性；余江斌等[7]利用天氣預(yù)報(bào)的可測因子建立了BP-ET0預(yù)測模型，發(fā)現(xiàn)其預(yù)測值與PM公式的計(jì)算值相接近，且操作簡單。這些都證明了BP神經(jīng)網(wǎng)絡(luò)在灌溉量預(yù)測上的應(yīng)用是可行和有效的[7]。但BP神經(jīng)網(wǎng)絡(luò)也存在一些缺點(diǎn)，如網(wǎng)絡(luò)結(jié)構(gòu)難以確定、易陷入局部極小值、收斂速度慢等[8]。針對這些局限性問題，已有學(xué)者提出了一些改進(jìn)措施以改善BP神經(jīng)網(wǎng)絡(luò)的不足，如人工蜂群算法(ABC)、遺傳算法(GA)、混合PSO-BP等群智能優(yōu)化算法被用于優(yōu)化BP神經(jīng)網(wǎng)絡(luò)以提高其性能[9-11]。但是這些算法仍無法克服局部極小值的問題。為此，本文提出了一種新的結(jié)合平均影響值算法(MIV)和改進(jìn)的非線性隨機(jī)遞減蝙蝠算法(Nonlinear Random Decreasing Bat Algorithm，NRDBA)的反向傳播算法(NRDBA-BP)，首先采用MIV算法和SPSS軟件對輸入變量進(jìn)行篩選，以優(yōu)化BP神經(jīng)網(wǎng)絡(luò)的結(jié)構(gòu)，降低計(jì)算量。同時(shí)利用改進(jìn)的蝙蝠算法對BP神經(jīng)網(wǎng)絡(luò)的初始閾值和權(quán)值進(jìn)行優(yōu)化，以提高BP神經(jīng)網(wǎng)絡(luò)的收斂速度和預(yù)測精度。

1 NRDBA-BP混合算法

1.1 BP神經(jīng)網(wǎng)絡(luò)原理

BP神經(jīng)網(wǎng)絡(luò)是采用誤差反向傳播算法訓(xùn)練的多層前饋網(wǎng)絡(luò)。該方法的訓(xùn)練速度快，非線性擬合能力強(qiáng)，是目前最流行和應(yīng)用最為廣泛的一種方法。大量研究表明，單隱含層的BP神經(jīng)網(wǎng)絡(luò)模型即可逼近任意精度的連續(xù)函數(shù)。BP神經(jīng)網(wǎng)絡(luò)的經(jīng)典結(jié)構(gòu)如圖1所示，由3層組成：輸入層、隱含層和輸出層。在BP神經(jīng)網(wǎng)絡(luò)模型中，信號的前向傳播和誤差的反向傳播交替執(zhí)行，通過梯度下降算法來使誤差函數(shù)逐漸減小[12]。一旦確立了目標(biāo)誤差和訓(xùn)練參數(shù)，網(wǎng)絡(luò)模型就可以通過連續(xù)校正各層之間的權(quán)重和閾值來逼近任意精度的非線性連續(xù)函數(shù)。誤差函數(shù)為

圖1 神經(jīng)網(wǎng)絡(luò)經(jīng)典結(jié)構(gòu)模型

(1)

其中Pn是該神經(jīng)網(wǎng)絡(luò)模型的計(jì)算輸出，On是期望的輸出，N是訓(xùn)練數(shù)據(jù)集的個(gè)數(shù)。

1.2 蝙蝠算法

fi=fmin+(fmax-fmin)·β，

(2)

(3)

(4)

式中β是隨機(jī)數(shù)且β∈[0，1]，x*是全局最佳位置，fmax、fmin是頻率的最大、最小值，t為迭代次數(shù)。

此外蝙蝠還采用另一種隨機(jī)移動(dòng)策略以更好地對搜索空間進(jìn)行探索。算法每次迭代的過程中隨機(jī)生成一個(gè)數(shù)λ，若λ>ri，則根據(jù)式(5)更新位置，否則隨機(jī)產(chǎn)生一個(gè)新的位置[14]：

(5)

(6)

(7)

1.3 改進(jìn)的蝙蝠算法(NRDBA)

文獻(xiàn)[13]已經(jīng)證明BA算法優(yōu)于其他基于進(jìn)化的群智能算法，如粒子群算法(PSO)、遺傳算法(GA)等。本文為了進(jìn)一步提高該算法的全局搜索和局部搜索能力，提出了一種新的改進(jìn)算法，將粒子群算法中慣性權(quán)重的思想引入到速度更新公式中。

通過對蝙蝠算法基本原理的分析，發(fā)現(xiàn)算法中蝙蝠的位置受到蝙蝠移動(dòng)速度的影響，同時(shí)蝙蝠當(dāng)前速度又受上一代移動(dòng)速度的影響。根據(jù)蝙蝠算法的基本公式可知，上一代蝙蝠的移動(dòng)速度前的系數(shù)始終為1，這將不利于算法前期對搜索空間的快速搜索和算法后期的局部精細(xì)搜索[15]。為此，本文借鑒粒子群算法的思想，將非慣性權(quán)重引入到速度更新公式中，同時(shí)為了避免算法后期陷入局部極值，在非線性遞減的基礎(chǔ)上引入隨機(jī)因子，使得每一代蝙蝠的個(gè)體速度具有一定的隨機(jī)性，在算法前期使種群快速地分布在整個(gè)搜索空間，在后期陷入局部極值時(shí)，隨機(jī)的速度增加了其跳出局部極值的機(jī)會。具體表達(dá)式如下：

w=wmax-(0.7+0.3·rand())·(wmax-wmin)·

(8)

式中wmax=0.9和wmin=0.4為慣性權(quán)重w的最大和最小值，rand()為在[0，1]上服從均勻分布的隨機(jī)數(shù)，t和tmax分別為算法的當(dāng)前運(yùn)行次數(shù)和最大運(yùn)行次數(shù)，betarnd()為服從beta分布的隨機(jī)數(shù)[16]，k∈[0.1，0.9]為慣性偏離因子，本文取0.8。

此時(shí)，慣性權(quán)重w在算法的前期數(shù)值較大，蝙蝠的飛行速度也大，此時(shí)蝙蝠具有出色的全局探測能力，可以快速地找到局部最優(yōu)解；相反當(dāng)后期w較小時(shí)，蝙蝠的局部搜索能力大大增強(qiáng)，可以在局部最優(yōu)解的附近進(jìn)行精細(xì)搜索[17]。同時(shí)，隨機(jī)因子的引入使得算法后期w也可能得到一個(gè)相對大的值，從而增加算法跳出局部最小值的概率，而且使得蝙蝠種群的速度在每代和每一代的個(gè)體之間都具有隨機(jī)性，增加了種群的多樣性，可以加快蝙蝠算法的收斂速度，進(jìn)一步避免算法后期陷入局部極值[18]。此時(shí)，速度更新公式為

(9)

1.4 NRDBA-BP神經(jīng)網(wǎng)絡(luò)模型

為了得到更加精確的ET0預(yù)測值，本文采用NRDBA算法對BP神經(jīng)網(wǎng)絡(luò)的初始參數(shù)進(jìn)行優(yōu)化，整個(gè)模型的算法實(shí)現(xiàn)流程如圖2所示。具體實(shí)現(xiàn)步驟如下：

圖2 NRDBA-BP算法流程圖

(1)確定BP神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)，若模型具有m個(gè)輸入層和n個(gè)隱含層，則模型的結(jié)構(gòu)為m-n-1，蝙蝠算法的搜索空間維數(shù)為[19]D=(m+1)·n+(n+1)；

(2)初始化蝙蝠算法的參數(shù)，對蝙蝠種群的數(shù)量N、個(gè)體位置x、速度v、脈沖頻率f、脈沖響度A和速率r隨機(jī)進(jìn)行初始化；

(3)確定適應(yīng)度函數(shù)：

(10)

(4)更新蝙蝠搜索脈沖、位置和速度，獲得局部最佳位置和全局最佳位置；

(5)生成隨機(jī)數(shù)λ，若λ>ri則根據(jù)式(5)對當(dāng)前的最優(yōu)解進(jìn)行隨機(jī)擾動(dòng)，否則通過隨機(jī)飛行產(chǎn)生一個(gè)新解；

(6)若[λ

(7)根據(jù)式(10)計(jì)算所有蝙蝠的適應(yīng)度值并對其排序，得到全局最佳位置的適應(yīng)度值fg，若fg達(dá)到網(wǎng)絡(luò)訓(xùn)練精度(ε)或當(dāng)前迭代次數(shù)達(dá)到最大值，則算法結(jié)束，否則轉(zhuǎn)到步驟4繼續(xù)更新蝙蝠的速度和位置；

(8)蝙蝠算法求出fg對應(yīng)的的全局最佳位置即是BP神經(jīng)網(wǎng)絡(luò)最佳的權(quán)值和閾值。

2 實(shí)驗(yàn)數(shù)據(jù)分析

2.1 實(shí)驗(yàn)數(shù)據(jù)來源

本研究使用的數(shù)據(jù)來源于中國氣象數(shù)據(jù)網(wǎng)提供的中國地面氣候資料日值數(shù)據(jù)集(V3.0)中商丘站2014-01-01至2016-06-17共900天的逐日風(fēng)速、日照時(shí)長、平均氣溫、最高氣溫、最低氣溫和平均相對濕度6個(gè)氣象資料，并以此作為原始數(shù)據(jù)，采用修正后的Peman-Monteith公式計(jì)算ET0，具體數(shù)據(jù)見表1。

表1 2014-01-01至2016-06-17商丘地區(qū)逐日氣象數(shù)據(jù)與ET0

本文中ET0參考值采用修正后的Peman-Monteith公式計(jì)算，公式為

(11)

其中ET0為參考作物騰發(fā)量(mm/d)，Rn是作物表面的凈輻射量(MJ/(m2·d))，G是地面熱通量(MJ/(m2·d))，T是2 m高度處的日平均氣溫(℃)，u2為2 m高度處的風(fēng)速(m/s)，es為飽和蒸汽壓(kPa)，ea為實(shí)際蒸汽壓力(kPa)，Δ為飽和水汽壓與溫度曲線的斜率(kPa/℃)，γ是濕度常數(shù)(kPa/℃)。

2.2 MIV特征值篩選變量

由表1數(shù)據(jù)可看出，利用PM公式計(jì)算ET0時(shí)涉及到的氣象因素較多，若將其直接用于BP神經(jīng)網(wǎng)絡(luò)建模，不僅會增加網(wǎng)絡(luò)結(jié)構(gòu)設(shè)計(jì)難度和算法的運(yùn)行時(shí)間，而且會對其精確度產(chǎn)生影響。因此選擇盡可能少的影響ET0的主要?dú)庀笠蛩赜脕斫Ｊ怯肂P神經(jīng)網(wǎng)絡(luò)預(yù)測ET0時(shí)另一迫切需要解決的問題。

在從所有因素中選擇關(guān)鍵決定因素時(shí)主成分分析和因子分析被認(rèn)為是有效的。但是使用這些分析時(shí)，需保證主成分之間相互獨(dú)立且各因素之間無相互作用。然而，本研究中各氣象因子之間具有明顯的相關(guān)性[20]，這與前提條件相沖突。而簡單的相關(guān)分析沒有考慮變量之間的相互作用，很可能導(dǎo)致錯(cuò)誤的結(jié)果。MIV算法被認(rèn)為是評估變量在神經(jīng)網(wǎng)絡(luò)中對結(jié)果影響的最佳指標(biāo)之一[21]。故本文將平均影響值(MIV)與BP神經(jīng)網(wǎng)絡(luò)相結(jié)合，篩選出對ET0影響較大的參數(shù)進(jìn)行建模，同時(shí)利用SPSS軟件進(jìn)行進(jìn)一步的相關(guān)性分析，避免出現(xiàn)具有強(qiáng)相互作用的變量參與建模。

MIV的主要原理是在網(wǎng)絡(luò)模型中訓(xùn)練自變量，通過觀察不同變量對結(jié)果的影響值來確定參數(shù)對預(yù)測結(jié)果的影響。其具體計(jì)算過程是：首先利用訓(xùn)練集訓(xùn)練BP神經(jīng)網(wǎng)絡(luò)，然后將訓(xùn)練集中每個(gè)變量分別增加10%和減少10%，得到兩個(gè)新的樣本集A1和A2。最后將A1和A2分別使用已經(jīng)訓(xùn)練好的神經(jīng)網(wǎng)絡(luò)進(jìn)行預(yù)測得出兩組預(yù)測值B1和B2，其差值即為影響變化值(IV)并對其求平均，其值即為平均影響值(MIV)。值得一提的是，MIV的絕對值大小代表各類參數(shù)對因變量產(chǎn)生的實(shí)際影響大小，正負(fù)值則代表相關(guān)的方向。此時(shí)可將MIV絕對值較大的參數(shù)挑選出來進(jìn)行建模，從而實(shí)現(xiàn)對自變量的篩選[22]。

求出各自變量的MIV值后，計(jì)算第i個(gè)自變量對因變量的相對貢獻(xiàn)率：

(12)

對上文中提到的影響ET0的6個(gè)氣象因素進(jìn)行MIV變量篩選，其結(jié)果如表2所示。

表2 各影響因素MIV特征值

根據(jù)MIV絕對值數(shù)值越大時(shí)其相對于網(wǎng)絡(luò)輸出結(jié)果所產(chǎn)生的影響也越大的原理，由表2可知，排名前4的因素為風(fēng)速、平均氣溫、最低氣溫和平均相對濕度。根據(jù)以往資料可知輸入變量的獨(dú)立性也會影響B(tài)P神經(jīng)網(wǎng)絡(luò)的擬合精度，因此本文利用SPSS 25軟件對其進(jìn)行進(jìn)一步的獨(dú)立性分析。各自變量之間的皮爾遜相關(guān)系數(shù)如表3所示[24]。

表3 各氣象因素與ET0相關(guān)系數(shù)

由表3可以看出風(fēng)速、平均氣溫、平均相對濕度之間彼此獨(dú)立，而最低氣溫和平均氣溫之間具有強(qiáng)相關(guān)性，且平均氣溫對ET0的影響程度大，故本文選擇風(fēng)速、平均氣溫和平均相對濕度作為輸入變量，ET0為輸出變量，且輸入變量三者的CCR為86.27%，滿足BP神經(jīng)網(wǎng)絡(luò)的建模條件。

2.3 輸入數(shù)據(jù)預(yù)處理

由于BP神經(jīng)網(wǎng)絡(luò)參與建模的各因素之間的數(shù)量級不同，若不對其進(jìn)行歸一化處理，將會導(dǎo)致較大的誤差。因此在BP神經(jīng)網(wǎng)絡(luò)的訓(xùn)練之前，將訓(xùn)練和預(yù)測的所有影響因素標(biāo)準(zhǔn)化到0～1的區(qū)間內(nèi)，即[25]

(13)

其中xnorm、x0、xmin、xmax分別是歸一化值、原始數(shù)據(jù)、最小值和最大值。

2.4 確定BP神經(jīng)網(wǎng)絡(luò)的結(jié)構(gòu)

由上一節(jié)可知BP神經(jīng)網(wǎng)絡(luò)的輸入和輸出變量分別為3和1，然而，至今尚未有一個(gè)明確的公式可以用來確定隱含層節(jié)點(diǎn)數(shù)。根據(jù)式(13)可以確定隱含層節(jié)點(diǎn)數(shù)的范圍，然后用試錯(cuò)法確定最優(yōu)的節(jié)點(diǎn)數(shù)[26]。

通常情況下，實(shí)驗(yàn)者會先基于經(jīng)驗(yàn)公式劃分粗略的范圍，然后將其代入神經(jīng)網(wǎng)絡(luò)進(jìn)行逐一的驗(yàn)證，從中選擇預(yù)測誤差最小的隱含層節(jié)點(diǎn)數(shù)[27]：

(14)

式中m為輸入層節(jié)點(diǎn)數(shù)，n為輸出層節(jié)點(diǎn)數(shù)，a∈[0，10]。

由式(14)可以確定隱含層節(jié)點(diǎn)數(shù)l的范圍為[2，12]，取該區(qū)間的整數(shù)值依次進(jìn)行實(shí)驗(yàn)，其結(jié)果如圖3所示。

圖3 不同隱含層節(jié)點(diǎn)數(shù)的預(yù)測誤差對比

由圖3可知，使BP神經(jīng)網(wǎng)絡(luò)誤差最小的隱含層節(jié)點(diǎn)數(shù)為6，故BP神經(jīng)網(wǎng)絡(luò)的結(jié)構(gòu)為3-6-1。

3 基于NRDBA-BP的ET0預(yù)測模型

根據(jù)以上分析，本文分別建立BP模型、BA-BP模型、NRDBA-BP模型和PSO-BP模型對ET0進(jìn)行預(yù)測。其中BA-BP和NRDBA-BP算法使用相同的參數(shù)設(shè)置：N=50，x∈[-2，2]，v∈[-5，5]，f∈[0，2]，A0=0.5ones(N，1)，α=0.95，r0=0.9ones(N，1)，γ=0.95。PSO-BP算法的參數(shù)為：學(xué)習(xí)因子c1=c2=2，慣性權(quán)重wmax=0.7，wmin=0.4，迭代次數(shù)tmax=100，粒子種群個(gè)數(shù)N=20。

本研究共收集了900組數(shù)據(jù)，使用其中的675組數(shù)據(jù)作為訓(xùn)練集，其余225組數(shù)據(jù)作為測試集。圖4為4種模型的預(yù)測結(jié)果與實(shí)際結(jié)果的對比圖，圖5為4種模型的預(yù)測誤差對比圖。同時(shí)為了更加明顯的說明4種模型的預(yù)測效果，本文選取以下統(tǒng)計(jì)指標(biāo)作為評價(jià)模型的依據(jù)，其計(jì)算公式如下：

(1)均方誤差：

(15)

(2)平均絕對誤差：

(16)

(3)決定系數(shù)：

(17)

從圖4中可以看出NRDBA-BP的預(yù)測結(jié)果與BP模型、BA-BP模型、PSO-BP模型相比更加接近真實(shí)值，說明NRDBA-BP模型比BP、BA-BP、PSO-BP模型具有更好的預(yù)測性能。圖5中不同模型預(yù)測的誤差對比進(jìn)一步說明了這個(gè)情況，從中可以看出在4種預(yù)測模型中，從整體上看NRDBA-BP模型的誤差與其他3種模型相比較小，表明其預(yù)測精度最高，與真實(shí)值最接近。

圖4 不同模型預(yù)測曲線比較

圖5 不同模型預(yù)測誤差對比

為了定量地確定最優(yōu)預(yù)測模型，采用上述3種統(tǒng)計(jì)指標(biāo)對模型進(jìn)行評價(jià)和比較。表4記錄了BP、BA-BP、PSO-BP和NRDBA-BP 4種模型預(yù)測ET0時(shí)的相關(guān)指標(biāo)。其中NRDBA-BP模型的MSE、MAE與BP模型相比降低了10.11%和8.41%，與BA-BP模型相比降低了5.9%和3.2%，與BA-BP模型相比降低了4.58%和3.16%，且4種預(yù)測模型中NRDBA-BP的R2最大，這表明NRDBA-BP模型在ET0預(yù)測方面優(yōu)于BP、BA-BP和PSO-BP算法。因此綜合分析，本文所提出的NRDBA-BP預(yù)測模型具有良好的預(yù)測性能和可信賴性，其預(yù)測值更加接近真實(shí)的ET0，說明NRDBA-BP模型在ET0預(yù)測方面具有很大的應(yīng)用潛力。

表4 模型指標(biāo)對比

4 結(jié)語

本文提出了一種基于非線性隨機(jī)遞減權(quán)重改進(jìn)的NRDBA算法，并將其應(yīng)用于BP神經(jīng)網(wǎng)絡(luò)預(yù)測ET0中以優(yōu)化其初始閾值和權(quán)值。該方法能夠有效解決傳統(tǒng)BP神經(jīng)網(wǎng)絡(luò)由于初始閾值和權(quán)值隨機(jī)選取所帶來訓(xùn)練速度慢、易陷于局部極值的問題。同時(shí)本文的另一個(gè)意義在于為ET0神經(jīng)網(wǎng)絡(luò)建模篩選關(guān)鍵因素提供了一種新的思路，通過MIV算法和SPSS軟件實(shí)現(xiàn)對變量的篩選，從而將互不相干的關(guān)鍵因素挑選出來參與建模。在此基礎(chǔ)上構(gòu)建了基于風(fēng)速、平均氣溫和平均相對濕度的三因素NRDBA-BP預(yù)測模型。實(shí)驗(yàn)結(jié)果表明，所提出的模型顯示出了較高的精度，對發(fā)展節(jié)能高效的新型灌溉技術(shù)具有一定的指導(dǎo)作用，同時(shí)所建立的模型對于減少工作量和利用較少氣象因素下建立預(yù)測模型具有重要意義。