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

姜 鋒，丁志宏，趙 焱

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

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

本文運用CEEMDAN方法對黑河鶯落峽水文站1945-2015年的年徑流量序列在年時間尺度上的多分辨率變化特征[1]進行研究，探討其各個波動分量在各自波動周期上的變化規律，以期為黑河流域水資源的開發、利用與保護工作提供宏觀的科學指導。

1 CEEMDAN基本理論

經驗模態分解[2](Empirical Mode Decomposition，EMD)是一種自適應性數據驅動的時頻分析方法。EMD可以根據一個非線性、非平穩時間序列的局部時變特征而將其分解為一系列平穩的局部分量，具有良好而穩定的分解—重構特性，這些局部分量具有不同的波動頻率的，其波動的頻率和幅值已經過調制，被稱為本征模態函數(Intrinsic Mode Function，IMF)。EMD方法是在數學原理上有別于傅里葉變換和小波變換的一種具有時頻多分辨率特征的數據平穩化方法，已在水文水資源領域得到了廣泛應用[3-6]。

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

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

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

(1)使用EMD計算x(i)=x+β0E1[w(i)](即x的第i次實現)的局部均值以求得第1個殘差：

r1=〈M(x(i))〉

(2)在第1階段(k=1)計算第1階模態：

(3)將r1+β1E2(w(i))實現的局部均值的平均值作為第2個殘差的估計值，將第2階模態定義為：

(4)對于k=3,…,K，計算第k個殘差：

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

(5)計算第k階模態：

(6)返回第(4)步計算下一個k。

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

綜上，經過重構之后，最終殘差滿足：

K為模態的總階數。故，原始序列x可以表示為：

2 實際應用

2.1 基本資料

本文分析選用的水文數據為黑河鶯落峽水文站1945-2015年實測年徑流量序列，見圖1。

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

2.2 年徑流量的CEEMDAN分解

運用CEEMDAN方法(為對比分析，同時采用EEMD方法)對圖1所示序列進行多時間尺度分解，限制標準差 的值取為0.25，結果見圖2～圖6。

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

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

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

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

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

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

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

由圖2～圖6可知：

(1)鶯落峽年徑流量時序具有5階模態，反映了鶯落峽上游的黑河山區流域產匯流系統徑流量在時域中演化的復雜尺度性。

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

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

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

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

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

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

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

2.3 對CEEMDAN分解結果的討論

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

3 結 語

鶯落峽年徑流量分別具有準2～6 a、準4～8 a、準9 a、準13 a、準28 a和準88 a波動周期，71 a來的鶯落峽年徑流量在波動中增加主要是因為大尺度模態的漲幅引起的，未來一段時期內的鶯落峽年徑流量預期將呈現在波動中衰減的總體變化趨勢。在流域社會經濟持續快速發展的現實背景下，這樣的水資源情勢預期對于黑河流域水資源管理工作提出了新的挑戰，更凸顯了實施最嚴格水資源管理制度的必要性和迫切性，以期在維系水生態系統良性循環的基礎上實現黑河流域社會經濟可持續發展。

□

[1] 王文圣，丁 晶，向紅蓮. 水文時間序列多時間尺度分析的小波變換法[J]. 四川大學學報(工程科學版)，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]. 中國農村水利水電，2008，(10)：36-38.

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

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

[6] 李 寧，岳德鵬，于 強，等. 磴口縣地下水埋深時空變化特征[J]. 南水北調與水利科技，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]. 北京大學學報(自然科學版)，2011，47(4)：627-635.

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

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

[11] 吳燕鋒，巴特爾·巴克，李 維，等.基于EEMD的杜尚別市1950-2013年降水多尺度分析[J].干旱區資源與環境，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] 李 鋒，林陽陽，晁蘇全，等.基于CEEMDAN與信息熵的液壓泵故障特征提取方法研究[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的黃河源區年徑流量多時間尺度變化特征研究[J]. 海河水利，2016，(6)：1-6.

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