基于CEEMDAN的黑河鶯落峽年徑流量多時(shí)間尺度變化特征研究

姜 鋒，丁志宏，趙 焱

(1.甘肅省水文水資源局, 蘭州 730000；2.中水北方勘測(cè)設(shè)計(jì)研究有限責(zé)任公司規(guī)劃設(shè)計(jì)處，天津 300222；3.黃河水利科學(xué)研究院水資源研究所，鄭州 450003)

黑河是我國(guó)第2大內(nèi)陸河，發(fā)源于青海省祁連縣，介于東經(jīng)96°42′～102°04′、北緯39°45′～42°40′，流域面積12.8 萬(wàn)km2，干流總長(zhǎng)928 km，其上游分東西2支，東支俄博河發(fā)源于俄博灘東的金瑤嶺，自東向西流，河長(zhǎng)80 km；西支野牛溝發(fā)源于鐵里干山，自西向東流，河長(zhǎng)190 km。東西2大支流在黃藏寺匯合后折向北流始稱黑河。出山口鶯落峽以上的祁連山區(qū)為黑河上游，海拔高程多在3 000～5 000 m，年降水量在350 mm以上，氣候陰濕，地貌景觀垂直地帶性鮮明，是黑河的產(chǎn)流區(qū)，鶯落峽水文站控制流域面積10 009 km2。

本文運(yùn)用CEEMDAN方法對(duì)黑河鶯落峽水文站1945-2015年的年徑流量序列在年時(shí)間尺度上的多分辨率變化特征[1]進(jìn)行研究，探討其各個(gè)波動(dòng)分量在各自波動(dòng)周期上的變化規(guī)律，以期為黑河流域水資源的開(kāi)發(fā)、利用與保護(hù)工作提供宏觀的科學(xué)指導(dǎo)。

1 CEEMDAN基本理論

經(jīng)驗(yàn)?zāi)B(tài)分解[2](Empirical Mode Decomposition，EMD)是一種自適應(yīng)性數(shù)據(jù)驅(qū)動(dòng)的時(shí)頻分析方法。EMD可以根據(jù)一個(gè)非線性、非平穩(wěn)時(shí)間序列的局部時(shí)變特征而將其分解為一系列平穩(wěn)的局部分量，具有良好而穩(wěn)定的分解—重構(gòu)特性，這些局部分量具有不同的波動(dòng)頻率的，其波動(dòng)的頻率和幅值已經(jīng)過(guò)調(diào)制，被稱為本征模態(tài)函數(shù)(Intrinsic Mode Function，IMF)。EMD方法是在數(shù)學(xué)原理上有別于傅里葉變換和小波變換的一種具有時(shí)頻多分辨率特征的數(shù)據(jù)平穩(wěn)化方法，已在水文水資源領(lǐng)域得到了廣泛應(yīng)用[3-6]。

但是， EMD在分解過(guò)程中可能會(huì)產(chǎn)生模態(tài)混淆現(xiàn)象(也稱混頻現(xiàn)象)，即相似的尺度存在于不同的IMF中或者在一個(gè)IMF中存在完全不同的尺度，而理想的分解結(jié)果應(yīng)該是每個(gè)IMF的變量尺度均是相似的，這主要是由于其算法所特有的局部篩選性質(zhì)制約而造成的。為了改善EMD的混頻現(xiàn)象，Huang等人提出了集合經(jīng)驗(yàn)?zāi)B(tài)分解[7](Ensemble Empirical Mode Decomposition，EEMD)，該方法是運(yùn)用總體平均思想來(lái)對(duì)添加了高斯噪聲的原始序列集合進(jìn)行EMD分解之后再求其集合平均值，增強(qiáng)了EMD的適應(yīng)性，以此作為EMD的最終模態(tài)分解結(jié)果。目前，EEMD在水文水資源領(lǐng)域也得到了廣泛應(yīng)用[8-11]，但是，EEMD也尚存在一些不足，諸如IMF分量存在殘留噪聲、每次EMD需要添加不同幅值的高斯噪聲、每次EMD可能產(chǎn)生不同的IMF個(gè)數(shù)等問(wèn)題，這使得最后的集合平均值求解存在困難。為此，Huang等人又通過(guò)使用增加互補(bǔ)高斯噪聲對(duì)的方法提出了互補(bǔ)集合經(jīng)驗(yàn)?zāi)B(tài)分解[12](Complementary Ensemble Empirical Mode Decomposition，CEEMD)，以此減輕IMF存在噪聲殘留的問(wèn)題。但是，CEEMD的數(shù)學(xué)完整性有待證明，計(jì)算工作量相應(yīng)增大，而且EMD分解可能產(chǎn)生不同個(gè)數(shù)IMF的問(wèn)題仍然未能解決。

為解決EEMD存在的上述問(wèn)題，Torres等人于2011年提出了具適應(yīng)性噪聲的完全集合經(jīng)驗(yàn)?zāi)B(tài)分解[13](Complete Ensemble Empirical Mode Decomposition with Adaptive Noise，CEEMDAN)，對(duì)EEMD進(jìn)行了顯著改進(jìn)并已在多個(gè)領(lǐng)域得到了應(yīng)用[14-16]，但是該方法仍存在分解早期可能存在虛假模態(tài)以及個(gè)別模態(tài)包含殘留噪聲等2個(gè)問(wèn)題[18]。2014年，Torres等人又對(duì)CEEMDAN進(jìn)行了改進(jìn)[17]，實(shí)例結(jié)果顯示改進(jìn)后的CEEMDAN方法可以精確地分解用來(lái)進(jìn)行測(cè)試的波形復(fù)雜的正弦疊加信號(hào)和電聲門信號(hào)，與EMD、EEMD、CEEMDAN初始算法相比，改進(jìn)后的CEEMDAN方法所分解得到的IMF分量不存在虛假模態(tài)和模態(tài)混疊，完美解決了CEEMDAN初始算法所存在的問(wèn)題。

CEEMDAN的具體算法為[17]：設(shè)x為待分解序列，令Ek(·)為通過(guò)EMD產(chǎn)生第k階模態(tài)的算子，令M(·)是產(chǎn)生待分解序列局部均值的算子，令w(i)是均值為0、方差為1的高斯噪聲，βk=ε0std (rk)，k≥1，βk=ε0std (x)/std {E1[w(i)]}>0，ε0為初始噪聲和原始序列的理想信噪比SNR的倒數(shù)，x(i)=x+w(i)，〈·〉是在實(shí)現(xiàn)中求取平均值的算子，可見(jiàn)E1(x)=x-M(x)，則：

(1)使用EMD計(jì)算x(i)=x+β0E1[w(i)](即x的第i次實(shí)現(xiàn))的局部均值以求得第1個(gè)殘差：

r1=〈M(x(i))〉

(2)在第1階段(k=1)計(jì)算第1階模態(tài)：

(3)將r1+β1E2(w(i))實(shí)現(xiàn)的局部均值的平均值作為第2個(gè)殘差的估計(jì)值，將第2階模態(tài)定義為：

(4)對(duì)于k=3,…,K，計(jì)算第k個(gè)殘差：

rk=〈M(rk-1+βk-1Ek(w(i)))〉

(5)計(jì)算第k階模態(tài)：

(6)返回第(4)步計(jì)算下一個(gè)k。

重復(fù)進(jìn)行第(4)步至第(6)步，直至所求得的殘差滿足以下條件之一時(shí)可停止計(jì)算：①滿足IMF條件；②局部極值點(diǎn)的個(gè)數(shù)小于3個(gè)；③不能被EMD進(jìn)一步分解為止。

綜上，經(jīng)過(guò)重構(gòu)之后，最終殘差滿足：

K為模態(tài)的總階數(shù)。故，原始序列x可以表示為：

2 實(shí)際應(yīng)用

2.1 基本資料

本文分析選用的水文數(shù)據(jù)為黑河鶯落峽水文站1945-2015年實(shí)測(cè)年徑流量序列，見(jiàn)圖1。

圖1 鶯落峽1945-2015年的年徑流量序列

2.2 年徑流量的CEEMDAN分解

運(yùn)用CEEMDAN方法(為對(duì)比分析，同時(shí)采用EEMD方法)對(duì)圖1所示序列進(jìn)行多時(shí)間尺度分解，限制標(biāo)準(zhǔn)差 的值取為0.25，結(jié)果見(jiàn)圖2～圖6。

由圖2～圖6可知，CEEMDAN和EEMD均將鶯落峽1945-2015年的年徑流量序列分解為5階模態(tài)，其中包括4個(gè)IMF分量(圖2～圖5)和1個(gè)趨勢(shì)項(xiàng)Res分量(圖6)。

以CEEMDAN分解結(jié)果為分析基準(zhǔn)(以下同)，采用納什效率系數(shù)(Nash-Sutcliffe efficiency coefficient，NSE)對(duì)2種方法的分解結(jié)果進(jìn)行對(duì)比分析，相應(yīng)的2個(gè)IMF1、IMF2、IMF3、IMF4、Res之間的納什效率系數(shù)分別為0.987 1、0.915 4、0.811 6、0.545 7、0.824 5，可見(jiàn)CEEMDAN較EEMD在低頻率(長(zhǎng)周期)尺度上的分解精度高，2者在第1階模態(tài)上的相似程度最高。

圖2 鶯落峽年徑流量序列的IMF1分量

圖3 鶯落峽年徑流量序列的IMF2分量

圖4 鶯落峽年徑流量序列的IMF3分量

圖5 鶯落峽年徑流量序列的IMF4分量

圖6 鶯落峽年徑流量序列的Res分量

由圖2～圖6可知：

(1)鶯落峽年徑流量時(shí)序具有5階模態(tài)，反映了鶯落峽上游的黑河山區(qū)流域產(chǎn)匯流系統(tǒng)徑流量在時(shí)域中演化的復(fù)雜尺度性。

(2)第1階模態(tài)IMF1是振幅最大、頻率最高、周期最短的一個(gè)波動(dòng)，依次下去的第2～5階模態(tài)的振幅逐漸減小、頻率逐漸降低、周期逐漸變長(zhǎng)。

(3)第1階模態(tài)IMF1具有準(zhǔn)2～6 a波動(dòng)周期，在71 a間，其平均振幅1.81 億m3，最大振幅3.46 億m3，最小振幅0.42 億m3。

(4)第2階模態(tài)IMF2具有準(zhǔn)4～8 a波動(dòng)周期，在71 a間，其平均振幅1.18 億m3，最大振幅2.27 億m3，最小振幅0.002 億m3；自1969-1980年為低幅振蕩時(shí)段，平均振幅0.59 億m3，最大振幅1.14 億m3，最小振幅0.20 億m3；自1992-2008年為低幅振蕩時(shí)段，平均振幅0.58 億m3，最大振0.90 億m3，最小振幅0.002 億m3。

(5)第3階模態(tài) IMF3具有準(zhǔn)9 a和準(zhǔn)13 a波動(dòng)周期，在71 a間，其平均振幅0.95 億m3，最大振1.62 億m3，最小振幅0.32 億m3。

(6)第4階模態(tài) IMF4具有準(zhǔn)28 a波動(dòng)周期，在71 a間，其平均振幅0.90 億m3，最大振1.14 億m3，最小振幅0.76 億m3。

(7)第5階模態(tài)Res分量顯示的是鶯落峽年徑流量的宏觀變化趨勢(shì)，1969年和2013年分別是Res曲線的谷值點(diǎn)和峰值點(diǎn)。1945-1969年，鶯落峽年徑流量總體呈衰減趨勢(shì)，減幅為1.65%；1970-2013年，鶯落峽年徑流量在總體上呈增加趨勢(shì)，增幅為2.11%；自2014年開(kāi)始，鶯落峽年徑流量總體上又開(kāi)始新一輪循環(huán)。以上述峰谷值年份所形成的半波動(dòng)周期推測(cè)，Res分量具有準(zhǔn)88 a波動(dòng)周期，振幅3.31 億m3。

(8)圖1中所示的鶯落峽年徑流量的增加主要是由第4模態(tài)IMF4和趨勢(shì)分量Res的增幅所形成的。

2.3 對(duì)CEEMDAN分解結(jié)果的討論

由圖2～圖6可見(jiàn)，準(zhǔn)28 a和準(zhǔn)88 a周期的大尺度模態(tài)在整體和全局上控制著整個(gè)鶯落峽年徑流量序列的變化過(guò)程，與之相比，準(zhǔn)2～6 a、準(zhǔn)4～8 a、準(zhǔn)9 a、準(zhǔn)13 a周期的中短尺度模態(tài)則以其更大的振蕩幅度而刻畫(huà)著整個(gè)鶯落峽年徑流量序列的變化細(xì)節(jié)。徑流量是流域產(chǎn)匯流系統(tǒng)的主要輸出變量，受降水量的影響最為顯著，而降水量的周期波動(dòng)又是全球和區(qū)域物理、氣候復(fù)雜巨系統(tǒng)演化的結(jié)果。據(jù)研究，鶯落峽年徑流量波動(dòng)的準(zhǔn)2～6 a、準(zhǔn)4～8 a、準(zhǔn)9 a、準(zhǔn)13 a、準(zhǔn)28 a周期與太陽(yáng)黑子和海-氣相互作用等有密切關(guān)系[19]。

3 結(jié) 語(yǔ)

鶯落峽年徑流量分別具有準(zhǔn)2～6 a、準(zhǔn)4～8 a、準(zhǔn)9 a、準(zhǔn)13 a、準(zhǔn)28 a和準(zhǔn)88 a波動(dòng)周期，71 a來(lái)的鶯落峽年徑流量在波動(dòng)中增加主要是因?yàn)榇蟪叨饶B(tài)的漲幅引起的，未來(lái)一段時(shí)期內(nèi)的鶯落峽年徑流量預(yù)期將呈現(xiàn)在波動(dòng)中衰減的總體變化趨勢(shì)。在流域社會(huì)經(jīng)濟(jì)持續(xù)快速發(fā)展的現(xiàn)實(shí)背景下，這樣的水資源情勢(shì)預(yù)期對(duì)于黑河流域水資源管理工作提出了新的挑戰(zhàn)，更凸顯了實(shí)施最嚴(yán)格水資源管理制度的必要性和迫切性，以期在維系水生態(tài)系統(tǒng)良性循環(huán)的基礎(chǔ)上實(shí)現(xiàn)黑河流域社會(huì)經(jīng)濟(jì)可持續(xù)發(fā)展。

□

[1] 王文圣，丁 晶，向紅蓮. 水文時(shí)間序列多時(shí)間尺度分析的小波變換法[J]. 四川大學(xué)學(xué)報(bào)(工程科學(xué)版)，2002，34(6)：15-17.

[2] N E Huang，Z Shen，S R Long，et al．The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[J]．Proceedings of the Royal Society A：Mathematical，Physical and Engineering Sciences，1998，454：903- 995.

[3] 劉奎建，丁志宏. 基于EMD的北洛河天然年徑流量變化特征分析[J]. 中國(guó)農(nóng)村水利水電，2008，(10)：36-38.

[4] 邵 駿，袁 鵬，顏志衡，等. 基于HHT的雅魯藏布江徑流變化周期及趨勢(shì)分析[J]. 中山大學(xué)學(xué)報(bào)(自然科學(xué)版)，2010，49(1)：125-130.

[5] 丁志宏，張金良，馮 平. 河流系統(tǒng)水沙變量的聯(lián)合分布模型研究[J]. 人民黃河，2012，34(8)：21-23，26.

[6] 李 寧，岳德鵬，于 強(qiáng)，等. 磴口縣地下水埋深時(shí)空變化特征[J]. 南水北調(diào)與水利科技，2017，15(3)：49-54,79.

[7] Z Wu, N E Huang. Ensemble empirical mode decomposition: a noise-assisted data analysis method[J]. Advances in Adaptive Data Analysis，2009，1(1):1-41.

[8] 王 兵，李曉東. 基于EEMD分解的歐洲溫度序列的多尺度分析[J]. 北京大學(xué)學(xué)報(bào)(自然科學(xué)版)，2011，47(4)：627-635.

[9] 姚欣明，陳元芳，顧圣華，等. EEMD-NNBR模型在降水預(yù)測(cè)中的應(yīng)用[J]. 水電能源科學(xué)，2014，32(12)：11-13，16.

[10] 張余慶，陳昌春，姚 鑫，等.江西省信江流域極端降水時(shí)空變化特征[J]. 水土保持研究，2015，22(4)：189-194，200.

[11] 吳燕鋒，巴特爾·巴克，李 維，等.基于EEMD的杜尚別市1950-2013年降水多尺度分析[J].干旱區(qū)資源與環(huán)境，2015，29(6)：152-157.

[12] J R Yeh, J S Shieh, N E Huang. Complementary ensemble empirical mode decomposition: a novel noise enhanced data analysis method[J]. Advances in Adaptive Data Analysis，2010，2(2):135-156.

[13] M E Torres, M A Colominas, G Schlotthauer, et al. A complete ensemble empirical mode decomposition with adaptive noise[M]∥ Proc. 36th IEEE Int. Conf. on Acoust., Speech and Signal Process, ICASSP 2011. Prague: Czech Republic, 2011:4 144-4 147.

[14] X Navarro, F Poree, G Carrault. EGG removal in preterm EGG combining empirical mode decomposition and adaptive filtering[C]∥ Proc. 37th, IEEE Int. Conf. on Acoust., Speech and Signal Process, ICASSP 2012. IEEE，2012:661-664.

[15] J Han, M Van der Bann. Empirical mode decomposition for seismic time-frequency analysis[J]. Geophysics，2013，78(2):9-19.

[16] 李 鋒，林陽(yáng)陽(yáng)，晁蘇全，等.基于CEEMDAN與信息熵的液壓泵故障特征提取方法研究[J]. 機(jī)床與液壓，2016，44(19)：192-195.

[17] M A Colominas, G Schlotthauer, M E Torres. Improved complete ensemble EMD: a suitable for biomedical signal processing[J]. Biomedical Signal Processing and Control，2014，14(11):19-29.

[18] 丁志宏，張金萍，趙 焱.基于CEEMDAN的黃河源區(qū)年徑流量多時(shí)間尺度變化特征研究[J]. 海河水利，2016，(6)：1-6.

[19] 藍(lán)永超，沈永平，林 紆，等.黃河上游徑流豐枯變化特征及其環(huán)流背景[J]. 冰川凍土，2006，28(6)：951-955.