基于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.