改進的EEMD-NNBR耦合模型在年徑流預測中的應用

鄭芳芳，王文圣，張嵐婷

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

中長期水文預測對防汛抗旱、水資源管理與保護等具有重要意義[1]。水文學者對中長期水文預測開展了大量研究，常用途徑有模糊分析法、灰色系統分析方法、小波分析方法[2-3]、人工神經網絡模型(ANN)、支持向量機(SVM)、最近鄰抽樣回歸模型(NNBR)及其相關耦合模型等[4-5]。基于相似原理建立的NNBR在水文中長期預測中進行了應用[6-7]，隨后提出了基于小波分析的NNBR模型、NNBR與SVM的耦合模型[8]、自適應NNBR和ANN的耦合模型[9]，研究表明NNBR模型及其組合預測精度較高，并且通常操作簡便，無需識別參數。集合經驗模態分解方法(Ensemble Empirical Mode Decomposition，EEMD)能將具有非線性、非平穩特點的徑流序列從不同時間尺度逐級分解出來[10-11]，分解序列能較好地反映原始序列真實的多時間演變特征和周期成分，具有良好的局部自適應性和直觀性。將EEMD與NNBR模型耦合，可充分利用2種方法的優勢，不失為一種新的組合方式。文獻[10]將該組合方式用于降水預測中并表明組合模型是可行的。但EEMD也存在不足，比如采用三次樣條插值對信號上、下極值點進行擬合時，由于極值分布不均勻導致擬合的包絡線不能正確反映信號的包絡形狀， 且端點處無法確定極大值和極小值，同時對極值點進行插值擬合時會產生較大的誤差[12]。針對這2個問題，本文采用極值中心三次樣條插值法對EEMD進行改善并與NNBR結合，建立了改進的EEMD-NNBR耦合模型。將改進的EEMD-NNBR耦合模型應用于金沙江屏山站年徑流預測中，研究表明建議模型是有效的。

1 改進的EEMD-NNBR耦合模型

1.1 改進的集合經驗模態分解方法

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

(1)

式中x0——序列前端延拓值；xn+1——序列后端延拓值。

改進的EEMD方法的步驟如下。

步驟1在延拓后的徑流序列xt(t=0,1, 2,…,n,n+1)中加入高斯白噪聲序列，即：

xt(i)=xt+nt(i)

( 2 )

式中xt(i)——第i次添加高斯白噪聲得到的序列。

高斯白噪聲添加次數NE通常取50或100，其標準差要根據原始序列標準差確定，一般取原始序列標準差的Nstd倍，Nstd根據信號來設置，沒有確定的計算公式。

白噪聲具有頻譜均勻分布的特性，在添加白噪聲后，序列中原本不同時間尺度的信號會自動分離到適當的參照尺度上去，避免隱含尺度的出現。

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

找出xt(i)序列所有的局部極大和極小值點。端點和極大值點之間、端點和極小值點之間分別用直線段連接，形成上、下包絡折線xmax和xmin。對所有極值點時刻的xmax和xmin取平均，得到極值中心點。再對極值中心進行三次樣條插值，擬合出平均包絡線再用減去平均包絡線，獲得去掉低頻序列的新序列。該過程即為極值中心三次樣條插值方法。對新序列重復上述操作，整個過程循環至平均包絡線趨近于零時停止，然后截斷兩端延拓部分，得到表示原始序列中最高頻的第一個IMF分量ct(i1)。一般IMF分量必須滿足：①過零點數和極值點數差值為0或1；②xmax和xmin每個時刻的平均值為零。用序列減去未截斷的ct(i1)得到剩余值序列r1，然后把r1序列作為一個新的原始序列，重復上述步驟得到ct(i2)，ct(i3)，…，ct(iM)，當剩余值序列為常數、單調函數或僅有一個極值的函數時，分解過程結束，此時的剩余值序列即為趨勢項rt(i)。

步驟3步驟1和步驟2重復NE次，對IMF分量進行平均運算，得到EEMD的IMF分量為：

( 3 )

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

此時，原始序列被EEMD方法分解成M個IMF和一個趨勢項：

( 4 )

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

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

NNBR模型無需預先假設隨機過程的分布規律，是依靠數據驅動的非參數模型。NNBR模型根據研究對象的不同分為單因子模型和多因子模型[14]，本文采用單因子模型對(j=1, 2,…,M)和趨勢項建模并進行預測。

設某一個分解序列yt(t=1, 2,…,n)，yt依賴于前p個歷史值yt-1,yt-2,…,yt-p。定義歷史模式矢量Dt=(yt-1,yt-2,…,yt-p)，其后續值為yt(t=p+1,p+2 ,… ,n)，p為矢量Dt的維數。已知當前模式矢量Di=(yi-1,yi-2,… ,yi-p)，按照水文相似性原理，在歷史模式矢量中找出與當前模式矢量最相似(最近鄰)的K個歷史模式矢量，然后將其按照相似程度由高到低進行排列：D1(i),D2(i),… ,DK(i)，與之對應的后續值記為yj(i)(j=1, 2 ,… ,K)。矢量維數一般由yt的自相關圖和偏自相關圖綜合確定[15]，同時要求。最近鄰模式矢量數K的選取一般采用K[16]。

當前模式矢量的后續值預測如下：

( 5 )

( 6 )

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

(7)

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

1.3 EEMD-NNBR耦合模型

將EEMD與NNBR結合構建EEMD-NNBR耦合模型，其基本步驟如下。

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

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

(8)

2 改進的EEMD-NNBR耦合模型在年徑流預測中的應用

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

圖1 屏山站1940—2016年的年徑流序列

2.1 模型建立

為充分利用信息，采用逐步向下滑動方式進行建模，即根據屏山站1940—2006年的年徑流序列預測2007年年徑流量，然后在序列中加入2007年的實測年徑流量，預測2008年年徑流量。以此類推，滑動預測2009—2016年的年徑流量。

以2007年年徑流量的預測過程為例，分別運用改進的EEMD方法和EEMD方法對屏山站1940—2006年的年徑流序列進行多時間尺度分解。通過參數試錯分析，白噪聲標準差取原始序列標準差的4倍，即取Nstd=4；添加100次高斯白噪聲，即NE=100。圖2所示，EEMD方法采用三次樣條插值方法得到的上、下包絡線存在多處過包絡和欠包絡，這樣導致2個同樣性質的極值點間的包絡線不單調，不能很好地反映序列的頻率特征。改進的EEMD方法分解年徑流序列時，包絡線為各極值點的折線連接，相比三次樣條插值方法得到的包絡線，過包絡和欠包絡現象明顯降低，且在進行端點延拓后，減緩了端點效應。年徑流序列分解后最終得到5個IMF分量以及一個趨勢項rt，見圖3、4。調用Matlab中hhspectrum函數求得的各IMF分量頻率統計特征值見表1、2。由表1可以看出，基于改進的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年的年徑流序列周期，結果表明存在3、9、12、17年左右的周期，可見改進EEMD方法識別周期成分是有效的，進一步說明其改進是有意義的。

圖2 折線包絡線與三次樣條包絡線比較

a) IMF1

b) IMF2

c) IMF3

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

e) IMF5

f) 趨勢項續圖3 改進的EEMD方法分解得到的IMF分量

a) IMF1

b) IMF2

c) IMF3

d) IMF4

e) IMF5

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

表1 改進的EEMD分解得到的IMF分量頻率統計特征值

表2 EEMD分解得到的IMF分量頻率統計特征值

將獲得的5個IMF分量和趨勢項數據代入NNBR模型進行預測，通過參數試錯分析，確定P=3，K=8。同理，利用1940—2007年的年徑流序列建立EEMD-NNBR模型并預測2008年年徑流量。以此類推，分別預測2009—2016年的年徑流量。

2.2 預測結果

屏山站年徑流預測結果見表3。采用平均相對誤差對模型改進前后的預測精度進行評定比較，分別為10.08%、8.59%。

表3 屏山站年徑流預測精度分析結果

2.3 對比分析與討論

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

EEMD-NNBR耦合模型在進行端點延拓和采用極值中心三次樣條插值的方法后，提高了對年徑流序列的預測精度。

3 結論

本文提出了改進的EEMD技術，改進技術采用極值中心三次樣條插值方法，減緩了EEMD分解過程出現的過包絡、欠包絡現象和端點效應。在此基礎上建立改進的EEMD-NNBR耦合模型，模型充分利用了EEMD分解非線性、非平穩序列的優勢和NNBR預測過程簡便、預測精度高的特點。以屏山站年徑流量預測為例，研究結果表明改進的EEMD-NNBR耦合模型在年徑流序列預測中模型適用性好，預測精度高，同時能良好地預測過程線的極值點。EEMD-NNBR耦合模型是在基于各種模型優點基礎上建立的一種模型，值得進一步加以推廣和應用。