鄭芳芳，王文圣，張嵐婷

(四川大學(xué)水利水電學(xué)院，四川 成都 610065)

中長(zhǎng)期水文預(yù)測(cè)對(duì)防汛抗旱、水資源管理與保護(hù)等具有重要意義[1]。水文學(xué)者對(duì)中長(zhǎng)期水文預(yù)測(cè)開(kāi)展了大量研究，常用途徑有模糊分析法、灰色系統(tǒng)分析方法、小波分析方法[2-3]、人工神經(jīng)網(wǎng)絡(luò)模型(ANN)、支持向量機(jī)(SVM)、最近鄰抽樣回歸模型(NNBR)及其相關(guān)耦合模型等[4-5]?；谙嗨圃斫⒌腘NBR在水文中長(zhǎng)期預(yù)測(cè)中進(jìn)行了應(yīng)用[6-7]，隨后提出了基于小波分析的NNBR模型、NNBR與SVM的耦合模型[8]、自適應(yīng)NNBR和ANN的耦合模型[9]，研究表明NNBR模型及其組合預(yù)測(cè)精度較高，并且通常操作簡(jiǎn)便，無(wú)需識(shí)別參數(shù)。集合經(jīng)驗(yàn)?zāi)B(tài)分解方法(Ensemble Empirical Mode Decomposition，EEMD)能將具有非線性、非平穩(wěn)特點(diǎn)的徑流序列從不同時(shí)間尺度逐級(jí)分解出來(lái)[10-11]，分解序列能較好地反映原始序列真實(shí)的多時(shí)間演變特征和周期成分，具有良好的局部自適應(yīng)性和直觀性。將EEMD與NNBR模型耦合，可充分利用2種方法的優(yōu)勢(shì)，不失為一種新的組合方式。文獻(xiàn)[10]將該組合方式用于降水預(yù)測(cè)中并表明組合模型是可行的。但EEMD也存在不足，比如采用三次樣條插值對(duì)信號(hào)上、下極值點(diǎn)進(jìn)行擬合時(shí)，由于極值分布不均勻?qū)е聰M合的包絡(luò)線不能正確反映信號(hào)的包絡(luò)形狀， 且端點(diǎn)處無(wú)法確定極大值和極小值，同時(shí)對(duì)極值點(diǎn)進(jìn)行插值擬合時(shí)會(huì)產(chǎn)生較大的誤差[12]。針對(duì)這2個(gè)問(wèn)題，本文采用極值中心三次樣條插值法對(duì)EEMD進(jìn)行改善并與NNBR結(jié)合，建立了改進(jìn)的EEMD-NNBR耦合模型。將改進(jìn)的EEMD-NNBR耦合模型應(yīng)用于金沙江屏山站年徑流預(yù)測(cè)中，研究表明建議模型是有效的。

1 改進(jìn)的EEMD-NNBR耦合模型

1.1 改進(jìn)的集合經(jīng)驗(yàn)?zāi)B(tài)分解方法

集合經(jīng)驗(yàn)?zāi)B(tài)分解方法(EEMD)非常適用于處理非線性、非平穩(wěn)的徑流序列，它可以將復(fù)雜的徑流序列分解成不同時(shí)間尺度的IMF分量，揭示其內(nèi)在的徑流周期。加入新的徑流信息后，分解過(guò)程可以自適應(yīng)調(diào)整[13]。但EEMD分解過(guò)程中，需找出序列所有的局部極大和極小值點(diǎn)，而序列端點(diǎn)處不可能正好是極值點(diǎn)，即使是極值點(diǎn)，也不可能既是極大值點(diǎn)又是極小值點(diǎn)，導(dǎo)致在進(jìn)行插值擬合時(shí)，端點(diǎn)處會(huì)產(chǎn)生較大擬合誤差。因此，為了充分利用端點(diǎn)值，防止得到的IMF分量在端點(diǎn)處產(chǎn)生畸變，在年徑流序列前后端各延拓一個(gè)序列均值。序列均值計(jì)算如下：

(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ù)信號(hào)來(lái)設(shè)置，沒(méi)有確定的計(jì)算公式。

白噪聲具有頻譜均勻分布的特性，在添加白噪聲后，序列中原本不同時(shí)間尺度的信號(hào)會(huì)自動(dòng)分離到適當(dāng)?shù)膮⒄粘叨壬先?，避免隱含尺度的出現(xiàn)。

步驟2對(duì)xt(i)進(jìn)行EEMD分解，得到IMFct(ij)(j=1,2,…M；M為IMF分量個(gè)數(shù))和趨勢(shì)rt(i)。ct(ij)為EEMD對(duì)序列xt(i)進(jìn)行分解得到的第j個(gè)IMF分量，rt(i)為EEMD對(duì)序列xt(i)進(jìn)行分解得到的趨勢(shì)項(xiàng)分量。由于三次樣條曲線良好的穩(wěn)定性和收斂性，EEMD分解過(guò)程中常采用三次樣條插值方法對(duì)上、下極值點(diǎn)進(jìn)行插值擬合，但通常年徑流序列的極值點(diǎn)分布并不均勻，導(dǎo)致擬合的包絡(luò)線出現(xiàn)過(guò)包絡(luò)和欠包絡(luò)的現(xiàn)象，不能反映正確的包絡(luò)形狀，從而造成較大的誤差，影響預(yù)測(cè)精度。為了改善這種現(xiàn)象，應(yīng)用極值中心三次樣條插值方法到EEMD分解年徑流序列過(guò)程中，其具體分解過(guò)程如下。

找出xt(i)序列所有的局部極大和極小值點(diǎn)。端點(diǎn)和極大值點(diǎn)之間、端點(diǎn)和極小值點(diǎn)之間分別用直線段連接，形成上、下包絡(luò)折線xmax和xmin。對(duì)所有極值點(diǎn)時(shí)刻的xmax和xmin取平均，得到極值中心點(diǎn)。再對(duì)極值中心進(jìn)行三次樣條插值，擬合出平均包絡(luò)線再用減去平均包絡(luò)線，獲得去掉低頻序列的新序列。該過(guò)程即為極值中心三次樣條插值方法。對(duì)新序列重復(fù)上述操作，整個(gè)過(guò)程循環(huán)至平均包絡(luò)線趨近于零時(shí)停止，然后截?cái)鄡啥搜油夭糠?，得到表示原始序列中最高頻的第一個(gè)IMF分量ct(i1)。一般IMF分量必須滿足：①過(guò)零點(diǎn)數(shù)和極值點(diǎn)數(shù)差值為0或1；②xmax和xmin每個(gè)時(shí)刻的平均值為零。用序列減去未截?cái)嗟腸t(i1)得到剩余值序列r1，然后把r1序列作為一個(gè)新的原始序列，重復(fù)上述步驟得到ct(i2)，ct(i3)，…，ct(iM)，當(dāng)剩余值序列為常數(shù)、單調(diào)函數(shù)或僅有一個(gè)極值的函數(shù)時(shí)，分解過(guò)程結(jié)束，此時(shí)的剩余值序列即為趨勢(shì)項(xiàng)rt(i)。

步驟3步驟1和步驟2重復(fù)NE次，對(duì)IMF分量進(jìn)行平均運(yùn)算，得到EEMD的IMF分量為：

( 3 )

式中ct(j)——第j個(gè)IMF分量，由EEMD對(duì)原始序列進(jìn)行分解得到。

此時(shí)，原始序列被EEMD方法分解成M個(gè)IMF和一個(gè)趨勢(shì)項(xiàng)：

( 4 )

式中rt——最后殘余分量，代表序列整體趨勢(shì)。

1.2 最近鄰抽樣回歸模型(NNBR)

NNBR模型無(wú)需預(yù)先假設(shè)隨機(jī)過(guò)程的分布規(guī)律，是依靠數(shù)據(jù)驅(qū)動(dòng)的非參數(shù)模型。NNBR模型根據(jù)研究對(duì)象的不同分為單因子模型和多因子模型[14]，本文采用單因子模型對(duì)(j=1, 2,…,M)和趨勢(shì)項(xiàng)建模并進(jìn)行預(yù)測(cè)。

設(shè)某一個(gè)分解序列yt(t=1, 2,…,n)，yt依賴于前p個(gè)歷史值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個(gè)歷史模式矢量，然后將其按照相似程度由高到低進(jìn)行排列：D1(i),D2(i),… ,DK(i)，與之對(duì)應(yīng)的后續(xù)值記為yj(i)(j=1, 2 ,… ,K)。矢量維數(shù)一般由yt的自相關(guān)圖和偏自相關(guān)圖綜合確定[15]，同時(shí)要求。最近鄰模式矢量數(shù)K的選取一般采用K[16]。

當(dāng)前模式矢量的后續(xù)值預(yù)測(cè)如下：

( 5 )

( 6 )

判別當(dāng)前模式Di與歷史模式Dt相似程度的指標(biāo)有歐氏距離、余弦相似度、耦合序位值、DTW距離等，文獻(xiàn)[14]研究表明基于歐式距離的NNBR模型預(yù)測(cè)精度較高，且歐式距離計(jì)算簡(jiǎn)便，故本文采用歐式距離指標(biāo)來(lái)選擇最近鄰模式矢量。歐式距離越小，則表示二者越近鄰(相似)，歐式距離計(jì)算公式為：

(7)

式中l(wèi)t(i)——Dt與Di間的歐氏距離；yij、ytj——Di、Dt的第j個(gè)元素。

1.3 EEMD-NNBR耦合模型

將EEMD與NNBR結(jié)合構(gòu)建EEMD-NNBR耦合模型，其基本步驟如下。

步驟1利用改進(jìn)的EEMD方法對(duì)年徑流時(shí)間序列進(jìn)行分解，得到年時(shí)間尺度由小到大的IMF分量序列(j=1, 2, …,M)和一個(gè)趨勢(shì)項(xiàng)rt。

步驟2分別對(duì)各IMF分量(j=1, 2, …,M)和趨勢(shì)項(xiàng)rt建立最近鄰抽樣回歸模型并進(jìn)行預(yù)測(cè)，將各預(yù)測(cè)值相加得到年徑流量的預(yù)測(cè)值，即：

(8)

2 改進(jìn)的EEMD-NNBR耦合模型在年徑流預(yù)測(cè)中的應(yīng)用

采用金沙江屏山站1940—2016年的年徑流序列(圖1)來(lái)驗(yàn)證EEMD-NNBR耦合模型的適用性。金沙江位于長(zhǎng)江上游，長(zhǎng)江河源有北源楚瑪爾河，正源沱沱河和南源當(dāng)曲，沱沱河與當(dāng)曲匯合處為通天河起點(diǎn)，終點(diǎn)為玉樹(shù)市巴塘河口；巴塘河口至宜賓市岷江河口稱為金沙江，以下始稱長(zhǎng)江。從青海江源到四川宜賓干流河長(zhǎng)3 481 km，流域面積約47.32萬(wàn)km2，多年平均降水量在600～1 500 mm，部分地區(qū)超過(guò)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年的年徑流序列

2.1 模型建立

為充分利用信息，采用逐步向下滑動(dòng)方式進(jìn)行建模，即根據(jù)屏山站1940—2006年的年徑流序列預(yù)測(cè)2007年年徑流量，然后在序列中加入2007年的實(shí)測(cè)年徑流量，預(yù)測(cè)2008年年徑流量。以此類推，滑動(dòng)預(yù)測(cè)2009—2016年的年徑流量。

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

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

a) IMF1

b) IMF2

c) IMF3

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

e) IMF5

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

a) IMF1

b) IMF2

c) IMF3

d) IMF4

e) IMF5

f) 趨勢(shì)項(xiàng)圖4 EEMD方法分解得到的IMF分量

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

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

將獲得的5個(gè)IMF分量和趨勢(shì)項(xiàng)數(shù)據(jù)代入NNBR模型進(jìn)行預(yù)測(cè)，通過(guò)參數(shù)試錯(cuò)分析，確定P=3，K=8。同理，利用1940—2007年的年徑流序列建立EEMD-NNBR模型并預(yù)測(cè)2008年年徑流量。以此類推，分別預(yù)測(cè)2009—2016年的年徑流量。

2.2 預(yù)測(cè)結(jié)果

屏山站年徑流預(yù)測(cè)結(jié)果見(jiàn)表3。采用平均相對(duì)誤差對(duì)模型改進(jìn)前后的預(yù)測(cè)精度進(jìn)行評(píng)定比較，分別為10.08%、8.59%。

表3 屏山站年徑流預(yù)測(cè)精度分析結(jié)果

2.3 對(duì)比分析與討論

由表4可知，總體上EEMD-NNBR模型不如改進(jìn)的EEMD-NNBR模型預(yù)測(cè)精度高。另外，由表3可知，對(duì)于2008年徑流極大值和2011年徑流極小值而言，改進(jìn)的EEMD-NNBR耦合模型預(yù)測(cè)的相對(duì)誤差分別為11.14%、29.56%，較改進(jìn)前的誤差13.10%和36.75%有明顯的降低，說(shuō)明改進(jìn)的EEMD-NNBR耦合模型對(duì)過(guò)程線的極值點(diǎn)預(yù)測(cè)效果更好。

EEMD-NNBR耦合模型在進(jìn)行端點(diǎn)延拓和采用極值中心三次樣條插值的方法后，提高了對(duì)年徑流序列的預(yù)測(cè)精度。

3 結(jié)論

本文提出了改進(jìn)的EEMD技術(shù)，改進(jìn)技術(shù)采用極值中心三次樣條插值方法，減緩了EEMD分解過(guò)程出現(xiàn)的過(guò)包絡(luò)、欠包絡(luò)現(xiàn)象和端點(diǎn)效應(yīng)。在此基礎(chǔ)上建立改進(jìn)的EEMD-NNBR耦合模型，模型充分利用了EEMD分解非線性、非平穩(wěn)序列的優(yōu)勢(shì)和NNBR預(yù)測(cè)過(guò)程簡(jiǎn)便、預(yù)測(cè)精度高的特點(diǎn)。以屏山站年徑流量預(yù)測(cè)為例，研究結(jié)果表明改進(jìn)的EEMD-NNBR耦合模型在年徑流序列預(yù)測(cè)中模型適用性好，預(yù)測(cè)精度高，同時(shí)能良好地預(yù)測(cè)過(guò)程線的極值點(diǎn)。EEMD-NNBR耦合模型是在基于各種模型優(yōu)點(diǎn)基礎(chǔ)上建立的一種模型，值得進(jìn)一步加以推廣和應(yīng)用。