鄭芳芳,王文圣,張嵐婷
(四川大學(xué)水利水電學(xué)院,四川 成都 610065)
中長期水文預(yù)測對防汛抗旱、水資源管理與保護(hù)等具有重要意義[1]。水文學(xué)者對中長期水文預(yù)測開展了大量研究,常用途徑有模糊分析法、灰色系統(tǒng)分析方法、小波分析方法[2-3]、人工神經(jīng)網(wǎng)絡(luò)模型(ANN)、支持向量機(jī)(SVM)、最近鄰抽樣回歸模型(NNBR)及其相關(guān)耦合模型等[4-5]。基于相似原理建立的NNBR在水文中長期預(yù)測中進(jìn)行了應(yīng)用[6-7],隨后提出了基于小波分析的NNBR模型、NNBR與SVM的耦合模型[8]、自適應(yīng)NNBR和ANN的耦合模型[9],研究表明NNBR模型及其組合預(yù)測精度較高,并且通常操作簡便,無需識別參數(shù)。集合經(jīng)驗?zāi)B(tài)分解方法(Ensemble Empirical Mode Decomposition,EEMD)能將具有非線性、非平穩(wěn)特點的徑流序列從不同時間尺度逐級分解出來[10-11],分解序列能較好地反映原始序列真實的多時間演變特征和周期成分,具有良好的局部自適應(yīng)性和直觀性。將EEMD與NNBR模型耦合,可充分利用2種方法的優(yōu)勢,不失為一種新的組合方式。文獻(xiàn)[10]將該組合方式用于降水預(yù)測中并表明組合模型是可行的。但EEMD也存在不足,比如采用三次樣條插值對信號上、下極值點進(jìn)行擬合時,由于極值分布不均勻?qū)е聰M合的包絡(luò)線不能正確反映信號的包絡(luò)形狀, 且端點處無法確定極大值和極小值,同時對極值點進(jìn)行插值擬合時會產(chǎn)生較大的誤差[12]。針對這2個問題,本文采用極值中心三次樣條插值法對EEMD進(jìn)行改善并與NNBR結(jié)合,建立了改進(jìn)的EEMD-NNBR耦合模型。將改進(jìn)的EEMD-NNBR耦合模型應(yīng)用于金沙江屏山站年徑流預(yù)測中,研究表明建議模型是有效的。
集合經(jīng)驗?zāi)B(tài)分解方法(EEMD)非常適用于處理非線性、非平穩(wěn)的徑流序列,它可以將復(fù)雜的徑流序列分解成不同時間尺度的IMF分量,揭示其內(nèi)在的徑流周期。加入新的徑流信息后,分解過程可以自適應(yīng)調(diào)整[13]。但EEMD分解過程中,需找出序列所有的局部極大和極小值點,而序列端點處不可能正好是極值點,即使是極值點,也不可能既是極大值點又是極小值點,導(dǎo)致在進(jìn)行插值擬合時,端點處會產(chǎn)生較大擬合誤差。因此,為了充分利用端點值,防止得到的IMF分量在端點處產(chǎn)生畸變,在年徑流序列前后端各延拓一個序列均值。序列均值計算如下:
(1)
式中x0——序列前端延拓值;xn+1——序列后端延拓值。
改進(jìn)的EEMD方法的步驟如下。
步驟1在延拓后的徑流序列xt(t=0,1, 2,…,n,n+1)中加入高斯白噪聲序列,即:
xt(i)=xt+nt(i)
( 2 )
式中xt(i)——第i次添加高斯白噪聲得到的序列。
高斯白噪聲添加次數(shù)NE通常取50或100,其標(biāo)準(zhǔn)差要根據(jù)原始序列標(biāo)準(zhǔn)差確定,一般取原始序列標(biāo)準(zhǔn)差的Nstd倍,Nstd根據(jù)信號來設(shè)置,沒有確定的計算公式。
白噪聲具有頻譜均勻分布的特性,在添加白噪聲后,序列中原本不同時間尺度的信號會自動分離到適當(dāng)?shù)膮⒄粘叨壬先ィ苊怆[含尺度的出現(xiàn)。
步驟2對xt(i)進(jìn)行EEMD分解,得到IMFct(ij)(j=1,2,…M;M為IMF分量個數(shù))和趨勢rt(i)。ct(ij)為EEMD對序列xt(i)進(jìn)行分解得到的第j個IMF分量,rt(i)為EEMD對序列xt(i)進(jìn)行分解得到的趨勢項分量。由于三次樣條曲線良好的穩(wěn)定性和收斂性,EEMD分解過程中常采用三次樣條插值方法對上、下極值點進(jìn)行插值擬合,但通常年徑流序列的極值點分布并不均勻,導(dǎo)致擬合的包絡(luò)線出現(xiàn)過包絡(luò)和欠包絡(luò)的現(xiàn)象,不能反映正確的包絡(luò)形狀,從而造成較大的誤差,影響預(yù)測精度。為了改善這種現(xiàn)象,應(yīng)用極值中心三次樣條插值方法到EEMD分解年徑流序列過程中,其具體分解過程如下。
找出xt(i)序列所有的局部極大和極小值點。端點和極大值點之間、端點和極小值點之間分別用直線段連接,形成上、下包絡(luò)折線xmax和xmin。對所有極值點時刻的xmax和xmin取平均,得到極值中心點。再對極值中心進(jìn)行三次樣條插值,擬合出平均包絡(luò)線再用減去平均包絡(luò)線,獲得去掉低頻序列的新序列。該過程即為極值中心三次樣條插值方法。對新序列重復(fù)上述操作,整個過程循環(huán)至平均包絡(luò)線趨近于零時停止,然后截斷兩端延拓部分,得到表示原始序列中最高頻的第一個IMF分量ct(i1)。一般IMF分量必須滿足:①過零點數(shù)和極值點數(shù)差值為0或1;②xmax和xmin每個時刻的平均值為零。用序列減去未截斷的ct(i1)得到剩余值序列r1,然后把r1序列作為一個新的原始序列,重復(fù)上述步驟得到ct(i2),ct(i3),…,ct(iM),當(dāng)剩余值序列為常數(shù)、單調(diào)函數(shù)或僅有一個極值的函數(shù)時,分解過程結(jié)束,此時的剩余值序列即為趨勢項rt(i)。
步驟3步驟1和步驟2重復(fù)NE次,對IMF分量進(jìn)行平均運算,得到EEMD的IMF分量為:
( 3 )
式中ct(j)——第j個IMF分量,由EEMD對原始序列進(jìn)行分解得到。
此時,原始序列被EEMD方法分解成M個IMF和一個趨勢項:
( 4 )
式中rt——最后殘余分量,代表序列整體趨勢。
NNBR模型無需預(yù)先假設(shè)隨機(jī)過程的分布規(guī)律,是依靠數(shù)據(jù)驅(qū)動的非參數(shù)模型。NNBR模型根據(jù)研究對象的不同分為單因子模型和多因子模型[14],本文采用單因子模型對(j=1, 2,…,M)和趨勢項建模并進(jìn)行預(yù)測。
設(shè)某一個分解序列yt(t=1, 2,…,n),yt依賴于前p個歷史值yt-1,yt-2,…,yt-p。定義歷史模式矢量Dt=(yt-1,yt-2,…,yt-p),其后續(xù)值為yt(t=p+1,p+2 ,… ,n),p為矢量Dt的維數(shù)。已知當(dāng)前模式矢量Di=(yi-1,yi-2,… ,yi-p),按照水文相似性原理,在歷史模式矢量中找出與當(dāng)前模式矢量最相似(最近鄰)的K個歷史模式矢量,然后將其按照相似程度由高到低進(jìn)行排列:D1(i),D2(i),… ,DK(i),與之對應(yīng)的后續(xù)值記為yj(i)(j=1, 2 ,… ,K)。矢量維數(shù)一般由yt的自相關(guān)圖和偏自相關(guān)圖綜合確定[15],同時要求。最近鄰模式矢量數(shù)K的選取一般采用K[16]。
當(dāng)前模式矢量的后續(xù)值預(yù)測如下:
( 5 )

( 6 )
判別當(dāng)前模式Di與歷史模式Dt相似程度的指標(biāo)有歐氏距離、余弦相似度、耦合序位值、DTW距離等,文獻(xiàn)[14]研究表明基于歐式距離的NNBR模型預(yù)測精度較高,且歐式距離計算簡便,故本文采用歐式距離指標(biāo)來選擇最近鄰模式矢量。歐式距離越小,則表示二者越近鄰(相似),歐式距離計算公式為:
(7)
式中l(wèi)t(i)——Dt與Di間的歐氏距離;yij、ytj——Di、Dt的第j個元素。
將EEMD與NNBR結(jié)合構(gòu)建EEMD-NNBR耦合模型,其基本步驟如下。
步驟1利用改進(jìn)的EEMD方法對年徑流時間序列進(jìn)行分解,得到年時間尺度由小到大的IMF分量序列(j=1, 2, …,M)和一個趨勢項rt。
步驟2分別對各IMF分量(j=1, 2, …,M)和趨勢項rt建立最近鄰抽樣回歸模型并進(jìn)行預(yù)測,將各預(yù)測值相加得到年徑流量的預(yù)測值,即:
(8)

采用金沙江屏山站1940—2016年的年徑流序列(圖1)來驗證EEMD-NNBR耦合模型的適用性。金沙江位于長江上游,長江河源有北源楚瑪爾河,正源沱沱河和南源當(dāng)曲,沱沱河與當(dāng)曲匯合處為通天河起點,終點為玉樹市巴塘河口;巴塘河口至宜賓市岷江河口稱為金沙江,以下始稱長江。從青海江源到四川宜賓干流河長3 481 km,流域面積約47.32萬km2,多年平均降水量在600~1 500 mm,部分地區(qū)超過1 800 mm。由屏山站1940—2016年徑流數(shù)據(jù)可知,金沙江流域徑流的年內(nèi)、年際分配極不均勻。連續(xù)最小三月(一般2—4月)占年徑流量6%~14%,連續(xù)最大三月(一般7—9月)占年徑流量41%~67%,多年的最大月徑流量與最小月徑流量比值在4~16之間。年徑流的變差系數(shù)CV為0.16,最大、最小年徑流的比值為1.95。

圖1 屏山站1940—2016年的年徑流序列
為充分利用信息,采用逐步向下滑動方式進(jìn)行建模,即根據(jù)屏山站1940—2006年的年徑流序列預(yù)測2007年年徑流量,然后在序列中加入2007年的實測年徑流量,預(yù)測2008年年徑流量。以此類推,滑動預(yù)測2009—2016年的年徑流量。
以2007年年徑流量的預(yù)測過程為例,分別運用改進(jìn)的EEMD方法和EEMD方法對屏山站1940—2006年的年徑流序列進(jìn)行多時間尺度分解。通過參數(shù)試錯分析,白噪聲標(biāo)準(zhǔn)差取原始序列標(biāo)準(zhǔn)差的4倍,即取Nstd=4;添加100次高斯白噪聲,即NE=100。圖2所示,EEMD方法采用三次樣條插值方法得到的上、下包絡(luò)線存在多處過包絡(luò)和欠包絡(luò),這樣導(dǎo)致2個同樣性質(zhì)的極值點間的包絡(luò)線不單調(diào),不能很好地反映序列的頻率特征。改進(jìn)的EEMD方法分解年徑流序列時,包絡(luò)線為各極值點的折線連接,相比三次樣條插值方法得到的包絡(luò)線,過包絡(luò)和欠包絡(luò)現(xiàn)象明顯降低,且在進(jìn)行端點延拓后,減緩了端點效應(yīng)。年徑流序列分解后最終得到5個IMF分量以及一個趨勢項rt,見圖3、4。調(diào)用Matlab中hhspectrum函數(shù)求得的各IMF分量頻率統(tǒng)計特征值見表1、2。由表1可以看出,基于改進(jìn)的EEMD方法,屏山站年徑流序列存在約為3、5、8、11、15年的周期。5個IMF分量中,IMF1的振幅相對較大,頻率相對較高;IMF2、IMF3、IMF4和IMF5的振幅相對較小,頻率較低。由表2可知,EEMD方法分解得到的年徑流序列存在約為3、6、14、38、65年的周期,由于資料長度有限,不能確定屏山站年徑流序列是否存在65年周期變化。5個IMF分量中,IMF1、IMF2和IMF3的振幅相對較大,頻率相對較高;IMF4、IMF5的振幅相對較小,頻率較低。由于2種方法分解得到的周期存在一定的差別,故本文運用小波分析來識別屏山站1940—2006年的年徑流序列周期,結(jié)果表明存在3、9、12、17年左右的周期,可見改進(jìn)EEMD方法識別周期成分是有效的,進(jìn)一步說明其改進(jìn)是有意義的。

圖2 折線包絡(luò)線與三次樣條包絡(luò)線比較

a) IMF1

b) IMF2

c) IMF3

d) IMF4圖3 改進(jìn)的EEMD方法分解得到的IMF分量

e) IMF5

f) 趨勢項續(xù)圖3 改進(jìn)的EEMD方法分解得到的IMF分量

a) IMF1

b) IMF2

c) IMF3

d) IMF4

e) IMF5

f) 趨勢項圖4 EEMD方法分解得到的IMF分量

表1 改進(jìn)的EEMD分解得到的IMF分量頻率統(tǒng)計特征值

表2 EEMD分解得到的IMF分量頻率統(tǒng)計特征值
將獲得的5個IMF分量和趨勢項數(shù)據(jù)代入NNBR模型進(jìn)行預(yù)測,通過參數(shù)試錯分析,確定P=3,K=8。同理,利用1940—2007年的年徑流序列建立EEMD-NNBR模型并預(yù)測2008年年徑流量。以此類推,分別預(yù)測2009—2016年的年徑流量。
屏山站年徑流預(yù)測結(jié)果見表3。采用平均相對誤差對模型改進(jìn)前后的預(yù)測精度進(jìn)行評定比較,分別為10.08%、8.59%。

表3 屏山站年徑流預(yù)測精度分析結(jié)果
由表4可知,總體上EEMD-NNBR模型不如改進(jìn)的EEMD-NNBR模型預(yù)測精度高。另外,由表3可知,對于2008年徑流極大值和2011年徑流極小值而言,改進(jìn)的EEMD-NNBR耦合模型預(yù)測的相對誤差分別為11.14%、29.56%,較改進(jìn)前的誤差13.10%和36.75%有明顯的降低,說明改進(jìn)的EEMD-NNBR耦合模型對過程線的極值點預(yù)測效果更好。
EEMD-NNBR耦合模型在進(jìn)行端點延拓和采用極值中心三次樣條插值的方法后,提高了對年徑流序列的預(yù)測精度。
本文提出了改進(jìn)的EEMD技術(shù),改進(jìn)技術(shù)采用極值中心三次樣條插值方法,減緩了EEMD分解過程出現(xiàn)的過包絡(luò)、欠包絡(luò)現(xiàn)象和端點效應(yīng)。在此基礎(chǔ)上建立改進(jìn)的EEMD-NNBR耦合模型,模型充分利用了EEMD分解非線性、非平穩(wěn)序列的優(yōu)勢和NNBR預(yù)測過程簡便、預(yù)測精度高的特點。以屏山站年徑流量預(yù)測為例,研究結(jié)果表明改進(jìn)的EEMD-NNBR耦合模型在年徑流序列預(yù)測中模型適用性好,預(yù)測精度高,同時能良好地預(yù)測過程線的極值點。EEMD-NNBR耦合模型是在基于各種模型優(yōu)點基礎(chǔ)上建立的一種模型,值得進(jìn)一步加以推廣和應(yīng)用。