基于EEMD的南京市降水特征分析投稿.doc

基于EEMD的南京市降水特征分析孙银凤基金项目，自然科学基金项目（NSFC-50979023），水利公益项目（201201026）作者简介：孙银凤( 1987)，女，硕士研究生，研究方向为水资源规划与管理，E-mail： 通讯作者：陆宝宏（1962.10-），男，安徽天长市人，副教授，主要研究方向为水文水资源，E-mail: .,电话:陆宝宏1,2 张巍1王猛1，黄洋1，王盼1，季妤1，穆尼热赛买提1 (1.河海大学 水文水资源学院 江苏 南京 210098;2.河海大学水文水资源与水利工程科学国家重点实验室，江苏 南京 210098)摘要：论文采用EEMD分解方法分析了南京气象站19512009年降水序列的多时间尺度，发现南京市年降水序列存在2年、67年、1415年和20年的准周期变化，降水在70年代之前较小，70年代之后降水较大，整个时间尺度上呈现先减小后升高的平均变化趋势。综合运用Hurst指数及降水距平分析降水的趋势变化结果与EEMD分析结论基本吻合，未来南京市降水可能呈现减少趋势。另外用小波分析对EEMD提取的周期进行验证，小波分析表明南京市降水序列包含2年、79年和1516年的准周期，与EEMD的结果基本一致，说明EEMD的结果具有一定的可信度，可以为该地区的降水预测提供依据。关键词： 降水；多时间尺度；EEMD分解；趋势分析 The characteristic Analysis of precipitation in Nanjing based on EEMD methodSUN Yin-feng1, LU Bao-hong1,2 ZHANG Wei,WANG Meng,HUANG Yang,WANG Pan,JI Yu,MU Ni-re（1.College of Hydrology and Water Resources, Hohai University, Nanjing 210098,China; 2. State key Laboratory of Hydrology Water Resources and Hydraulic Engineering, Hohai University, Nanjing 210098,China）Abstract：Based on the precipitation time series from 1951 to 2009 in Nanjing Weather Observation, Ensemble Empirical Mode Decomposition was used to analyze temporal scales and trend of the annual precipitation. The result shows that the annual precipitation with 2 year, and 6 7year short periods, as well as 20 year of long-term oscillation periods. Hurst index and the precipitation anomaly analysis reveal that the precipitation has a decrease trend before 1970, and then the Precipitation began to increase, lightly. Furthermore, the wavelet method was employed to analyze the precipitation period, it can be found that the annual precipitation has 2 year，79 year and 1516 year periods，which is consistent with the EEMD results, basically.Keywords：Precipitation; multiple time scale; Ensemble Empirical Mode Decomposition; Trend analysis1、引言气候变化以及气候变化对水文气象过程的影响是当今世界的研究热点之一。联合国政府间气候变化专门委员会在第四届评估报告中指出：过去100年（19062005年）全球平均气温升高0.74，预计21世纪末期，全球地表温度平均增温1.16.4（与19801999年相比）1。全球气候变化加之人类活动的影响，近年来极端水文事件频繁发生，降水的变化规律也因此越来越受人们关注。研究降水的周期及变化特征不仅对深入理解气候变化过程与规律，进行气候预测具有重要意义，还可以为区域水资源的开发与合理利用、防汛抗旱提供科学依据。目前，降水特征分析是水文气象上的一个热点问题，很多学者对降水进行了多方面的研究2-4，并取得一定成果。1998年，黄锷提出了一种新的数据处理方法，即经验模态分解（Empirical Mode Decomposition， EMD），它能够在时域和频域上提取数据的信息，并且不需要任何的先验假设，能够保留数据本身的特性，具有良好的自适应性和高效性，非常适合处理非线性、非平稳信号5。EMD方法提出后被广泛应用于水文要素分析上6-8。但是EMD分解存在一个重要缺陷就是模态混淆，为解决这个问题，Wu和Huang在对白噪声进行EMD分解深入研究的基础上，提出了总体平均经验模态分解方法（Ensemble Empirical Mode Decomposition，EEMD）9。与EMD分解相比，EEMD能够提供更好的尺度分解，固有模态函数（Intrinsic Mode Function，IMF）分量可以准确体现时间序列的内在波动特征和趋势变化，在水文气象上应用的效果显著10-11。鉴于此，本文采用EEMD分解对南京市19512009年的年降水序列进行多时间尺度和趋势分析，并结合Hurst指数和降水距平分析、小波分析对结果进行对比验证，为今后南京市降水预报和水资源规划提供参考。2、 方法介绍2.1 EMD的基本原理EMD分解是一种全新的处理非线性、非平稳时间序列的方法，它将整个时间序列分解成一系列的固有模态函数（IMF），每个IMF都满足两个条件：（1）整个时间范围内，每个IMF的局部极值点数和过零点数相等或者最多相差一个；（2）在任意时刻点，局部最大值的包络(上包络线)和局部最小值的包络(下包络线) 平均必须为零。分解的过程基于如下的假设：信号至少有两个极值；数据的局部时域特性是由极值点间的时间尺度唯一确定；如果数据没有极值点但有拐点，则可以通过对数据微分一次或多次求极值，然后再通过积分来获得分解结果。EMD分解的本质是通过数据的特征时间尺度来获得本征波动模式，然后分解数据，可以形象地称之为“筛选(sifting)”过程。具体分解步骤如下：（1）识别出信号x（t）中的所有极值（极大值和极小值）；（2）用三次样条曲线连接局部极大值，形成上外包线emax（t），对极小值进行同样的操作形成下外包线emin（t）；（3）计算上下外包线的均值m1(t)=(emax(t)+emin(t)/2；（4）计算h1（t）=x(t)-m1(t)；（5）将h1（t）作为新的信号序列，重复上述操作直到h1（t）成为一个零均值过程；（6）得出零均值过程h1（t）后，将其作为第一个IMF分量，用c1（t）表示，它表示原始信号最高频的分量；（7）从信号x（t）中减去c1（t），得到去除一个高频分量的新信号过程r1（t），即 （1）将r1（t）作为原始数据，重复（1）（6）筛选过程得到第二个IMF分量c2（t）；如此继续，可以得到c3（t）、c4（t）等，（2）一直到满足迭代条件为止（一般限定趋势项rn（t）为单调函数）。根据式（3）表达的完整性检查，我们可以通过所有的IMF分量和剩余趋势项r来重建数据。（3）式中，rn（t）为趋势项，代表了信号的平均趋势。各个IMF分量c1（t），c2（t）cn(t)分别包含了原始信号包含的不同时间尺度的波动，表示的周期振荡依次由小到大。因此，EMD分解实际上是将原始信号分解为各种不同特征尺度波动，每一个IMF分量可是是线性也可以是非线性的，代表信号的局部特征，每一个IMF都有实际的物理背景相对应。分解的最后一项可以反映数据内在的和自适应的趋势，代表时间序列在整个时间跨度上的变化程度。2.2 EEMD分解原理EMD分解的一个主要挑战是模式混淆（Mode Mixing，MM）现象的频繁出现。模式混叠主要有两种情况，一种是同一个IMF分量表达出两种或两种以上差别较大的时间尺度下的信号成分，另一种情况是两个不同IMF分量表达出相近的时间尺度下的信号成分。Huang指出：MM产生的主要原因是信号的间断性（不连续性）。MM是指EMD分解结果在时频分布上产生严重的锯齿线，EMD分解产生合理IMF分量的能力取决于信号极值点的存在以及极值点的分布间隔，如果信号没有足够的极值点，分解将停止；如果信号中极值点分布间隔不均与，会产生极值点上、下包络线的拟合误差，从而产生模式混淆现象。EEMD是针对EMD方法的不足, 提出的一种噪声辅助数据分析方法。EEMD的原理是利用白噪声具有频率均匀分布的统计特性，在原始信号中加入一定比例的白噪声后进行EMD分解，经过多次计算，最后取平均使加入的白噪声相互抵消，不仅保留了原信号的信息，还可以克服模态混淆问题12。具体步骤如下：（1） 将原始信号加上强度很低的白噪声。（2） 将附加噪声的信号进行EMD分解，得到一组IMF分量。（3） 重复执行上述（1）（2），每次添加不同的白噪声序列，共重复N次。（4） 为了能够消除人工添加噪声的影响，将N次EMD分解得到的IMF分量相加，取其均值作为最终的IMF分量：（4）式中，cj（t）表示对原始信号进行EEMD分解后所得到的第j个IMF分量。本文在EEMD过程中，每次添加白噪声的振幅为合成序列标准差的20%，集合平均次数为200次。采用镜像拓延的方法对数据进行处理，以减小边界效应的影响。3、 实例应用3.1 资料 资料来自国家气象中心整理提供的南京站19512009年的逐日降水资料。将逐日降水资料整理成年降水时间序列，图1给出了南京市59年的降水过程。其中，最大年降水量为1825.8mm，发生在1991年；最小年降水量为534mm，出现在1978年；多年平均降水量为1056mm。图1 南京市19512009年降水序列3.2 南京市降水序列EEMD分解结果分析采用EEMD方法对南京市19512009年的降水序列进行逐步分解，得到4个固有模态函数和一个趋势项（图2），从图中可以看出：（1）南京站19512009年的年降水序列可以分解为4个具有不同波动周期的IMF分量和1个趋势分量r，说明南京市降水59年的降水序列包含多个时间尺度特征。（2）IMF1以及IMF2与原始降水序列具有较高的一致性。IMF1振幅和波动幅度最大。其中，60年代之前振幅波动较大，到7080年代波动相对较小，80年代后期振幅又开始出现较大的波动，整体表现出23年的准周期。IMF2表现出67年左右的准周期，60年代之前变化平稳，波动比较小。60年代之后波动加剧，这种波动一直持续到21世纪。（3）随着阶数的增加，IMF分量的振幅逐渐减小，波长逐渐增加。IMF3具有1415年的准周期，波动幅度在80年代中期之前较小，之后呈现增大趋势。IMF4包含20年的年代际周期波动，波动幅度在80年代初期略有增长，整个时间跨度上波动幅度较稳定。（4）趋势项r表示降水序列整体的变化趋势，图2中趋势项r在整个时间尺度上呈现先减小后增加的趋势变化。以70年代初为分界点，降水由减小趋势变为逐渐增长，没有表现出明显的周期变化，它可能属于更长的周期的组成部分，但是由于时间序列长度的限制，要研究这一时间尺度的变化特征需要更长的降水资料。图2 降水序列EEMD分解结果将EEMD分解的IMF分量和趋势项r相加重构的降水序列，并与原始降水序列进行比较，发现二者相对误差都在1%之内，说明EEMD保留了EMD的特性，能够很好的重构序列，保持原时间序列的变化特征。图3 EEMD分解量重构序列与原降水序列的相对误差分布图利用FFT求平均周期的方法13对各IMF分量进行分析，提取相应的平均周期。图4给出了IMF1的能量-周期关系图，从图中可以看出，能量主要集中在2年的小尺度之内，可以认为IMF1是周期为2年的信号。同样，采用FFT对其他IMF分量进行分析，结果表明IMF2、IMF3和IMF4表示的周期分别为6.6年、14.8年和19.7年，说明南京市降水序列包含多个时间尺度。图4 基于FFT的IMF1周期能量关系图3.3 小波分析结果为验证EEMD分解得到的周期，采用水文气象上常用的morlet连续复小波对南京市的年降水序列进行分析。为了能够提取出降水序列的实际波动，采用降水的距平值进行分析，并用对称外延法将降水序列扩展，以减弱“头部影响”在序列两端产生的误差14-15。图5给出了小波系数实部图和小波方差图。小波系数实部时频图反应了降水序列的时间尺度变化。图中实线表示正值，对应降水量偏多，值越大降水越多；虚线表示负值，对应降水量偏小；值为零则意味着突变点。从图5中可以看出南京市降水存在时间尺度上的复杂嵌套，其中以年际振荡最为强烈，10年内的周期变化最为明显，主要的中心尺度为8年和2年，整个时间跨度上呈现正负相位交替出现，降水量经历丰枯交替变化。代际变化以10年20年尺度为主，三个振荡中心分别位于50年代、70年代以及90年代初期，中心尺度为15年。小波方差图反应了能量随时间尺度的分布，可以确定一个时间序列中各种扰动的相对强度，对应峰值处的尺度为时间序列的主要时间尺度，反应序列的主周期。根据小波方差图可以看出，南京市降水量的小波方差存在3个峰值，依次对应着2年、7-9年以及15-16年的准周期，以79年的年纪周期最为突出。综上，小波分析的结果显示南京市19512009年年降水序列主要存在2年、79年和1516年的周期变化，这与EEMD分析出的结果吻合，说明EEMD在处理降水序列上具有一定的可信度。图5 小波系数实部a和小波方差图3.4 降水趋势分析R/S分析法最初由英国水文学家赫斯特（Hurst，1951年）在研究尼罗河水坝工程时提出，可以用来分析时间序列的分形特征和长期记忆过程，是定量描述时间序列自相似性与长期依赖性的有效方法，其具体原理可参照文献16。R/S分析法得出的Hurst指数H取值为01，当H=0.5时，表明时间序列为相互独立、方差有限的随机序列；当0.5H1时，说明时间序列变化具有持续性，未来的变化与过去一致；当H0.5时，意味着时间序列未来的变化与过去具有反持续性，H越小，反持续性越强。图6给出了南京市降水的Hurst指数分析结果，可以看出，年降水的Hurst指数为0.43891，小于0.5，说明年降水序列具有反持续性，未来降水将呈现与过去相反的趋势。图6 Hurst指数分析结果表1给出了南京市降水序列的年均降水年代距平值。整个年代际内降水距平百分率均在10%之内，说明各年代降水整体变化不大，均在平均值附近波动。5070年代，降水距平值都小于0，说明降水偏少，70年代之后降水距平变为正值，表明降水偏多，与EEMD分析的结果一致。综上，结合hurst指数及EEMD的分析结果，预计南京市未来年降水可能出现减小趋势。表1 19512009年降水年代距平年代506070809021世纪距平-5.78-68.95-1.730.3333.1342.34距平百分比-0.55%-6.53%-0.16%0.03%3.14%4.01%4、 结论（1）论文基于一种新的信号处理方法EEMD分解，将它用到降水序列分析中，即先采用EEMD将降水序列分解成4个IMF分量和一个残留项r，在利用FFT求平均周期的方法分别对IMF分量求周期，并根据趋势项r分析降水序列的整体趋势，结果表明：南京市19512009年的年降水序列具有2年、67年、1415年和20年的准周期变化，降水在整个时间尺度上呈先减小后增加的趋势变化。（2）采用水文气象上常用的morlet复小波对南京市降水序列进行分析验证，发现降水序列中包含2年、79年和1516年的多时间尺度，与EEMD分解的结果吻合，说明EEMD在分析降水周期上具有一定的可靠性。（3）结合Hurst指数及降水距平分析南京市年降水的趋势表明，年降水在70年代之前距平偏小，70年代之后距平偏大，整体表现先减小后增加的趋势，Hurst指数小于0.5，说明降水序列具有反持续性，预测未来降水出现下降趋势。（4）EEMD方法能够将原始的时间序列分解成若干个不同时间尺度的序列，各个组分序列分别体现原始序列的局部化特性，每个组分的影响因子不同。分解得到的IMF分量的数量以及每个IMF分量的频率和振幅由时间序列本身决定，由此研究时间序列的周期更简单、方便、符合实际。（5）水文时间序列的变化规律十分复杂，本文通过对南京市降水序列的研究表明EEMD方法在水文时间序列的多时间尺度和趋势分析上具有很好的效果。但是目前的应用研究仅仅是初步的，怎样确定每个时间尺度的权重并将其应用到水文时间序列的预测中，将是以后研究工作的重点。参考文献： 1 IPCC. Climate change 2007: Impacts, adaptation, and vulnerability: Contribution of Working Group II to the froth assessment report of the intergovernmental panel on climate changeM. Cambridge：Cambridge University Press，2007. 2 高文荣,段春青,陈晓楠,邱林. 基于云推理的年降水预测模型J. 中国农村水利水电,2009,(3). 3 刘德林,刘贤赵,魏兴华,高广路. 基于子波变换的烟台市降水变化特征的多时间尺度研究J. 中国农村水利水电,2007,(11). 4 李帅,刘冀,杜发兴,汤振权. 小波分析在宜昌地区降水变化中的应用J. 中国农村水利水电,2009,(5). 5 Huang N E, Shen Z, Long S R, et al. The Empirical Mode Decomposition and the Hilbert Spectrum for Nonlinear and Non-Stationary Time Series Analysis J. Proceedings: Mathematical, Physical and Engineering Sciences. 1998, 454(1971): 903-995. 6 曹丽青,林振山,纪玲玲. 基于HHT方法的洪家渡水电站年平均流量多尺度演变分析J. 河海大学学报(自然科学版),2008,(4). 7 冯平,丁志宏,韩瑞光,张建伟. 基于EMD的降雨径流神经网络预测模型J. 系统工程理论与实践,2009,(1). 8 邵骏,袁鹏,颜志衡,张文江. 基于HHT的雅鲁藏布江径流变化周期及趋势分析J