基于EEMD的盐城大丰区降水量时间序列多尺度分析

陈 超，李鹏飞，陈尚云，陈 凯，陈则煌

(1.江苏省盐城市大丰区气象局，江苏盐城 224100；2.中国气象局气溶胶与云降水重点开放实验室，江苏南京 210044；3.南京信息工程大学大气物理学院，江苏南京 210044)



基于EEMD的盐城大丰区降水量时间序列多尺度分析

陈 超1，李鹏飞2，3，陈尚云1，陈 凯1，陈则煌2，3

(1.江苏省盐城市大丰区气象局，江苏盐城 224100；2.中国气象局气溶胶与云降水重点开放实验室，江苏南京 210044；3.南京信息工程大学大气物理学院，江苏南京 210044)

摘要[目的] 研究盐城大丰地区年降水量时间序列规律。 [方法] 基于集合经验模态分解(EEMD)的数据处理方法，通过Hilbert变换对降水量数据进行多尺度特征研究，并对研究结果进行显著性和误差检验。 [结果] 近55 a盐城大丰地区年降水量序列可以分解成4个IMF分量和1个趋势分量；主要存在0～0.05和0.45～0.50 Hz 2个频率区间，其间主周期为2.12 a年际周期和27.50 a年代际周期，二者均通过95%显著性检验。EEMD可以有效地重构原始降水量数据变化曲线，与原始数据对比误差百分比始终控制在±0.6%范围内，由白噪声带来的干扰误差几乎可以完全忽略。 [结论] 该研究可为合理开发与充分利用水文资源、防汛抗旱等提供科学依据。

关键词年降水量；时间序列；EEMD；Hilbert变换；显著性检验；盐城大丰地区

近年来，全球温度升高和气候变化对水文气象的影响受到各国越来越多的关注[1]，其中关于降水量多尺度变化的周期性规律研究已成为当今学者们研究的热点之一[2]。随着统计学与计算数学等领域的飞速发展，基于多种研究方法和辅助计算工具对水文历史累积数据深入地再分析有了很大的推进作用。吴亚楠等[3]利用滑动平均、小波分析、M-K突变检验等方法对第二松花江流域65 a间的降水量序列变化特征进行了系统的分析。李光强等[4]利用小波变换方法分析了近41 a来河北地区降水序列的周期变化，得出河北地区中东部、西部以及北部不同地区的周期性和差异性。以上主要采用小波分析法等[3-5]对降水量的周期性进行分析。然而小波分析依然没脱离海森堡测不准原理的束缚，某种尺度下不能在时间和频率上同时具有较高的精度，以及小波是非适应性的导致周期计算不够准确。HHT(Hilbert Huang Transform)技术，作为一种基于Hilbert变换的新的时域分析方法，由于具有较好的自适应性、可以对非线性非平稳的信号进行分析等优点，被许多学者所推崇，已成为取代小波分析方法研究水文气象的周期及变化的主流方法。然而，随着对HHT-Huang的核心经验模式分解算法(Empirical Mode Decomposition,EMD)的深入研究，发现其存在一定的模态混叠问题，导致进行的Hlibert变换所得的瞬时频率与Hlibert谱的准确度降低[6]。近年，Huang等[7]与法国的Flaridrin研究小组尝试利用白噪声来对信号中断所缺失的尺度进行弥补，创造性地提出了基于集合经验模态分解(EEMD)改进HHT的全新方法，这很大程度上克服了模态混叠的问题，保留了原序列的信息。因此，EEMD方法势必将在水文气象领域中得到广泛的运用，为合理开发与充分利用水文资源、防汛抗旱等提供科学依据。笔者基于改进型经验模解-整体经验模态分解法(EEMD)，并以1960～2014年盐城大丰地区年降水量时间序列为例进行应用，探讨近55 a盐城大丰地区年降水量时间序列的变化规律与周期特征，并对所得结果进行置信检验和误差分析。

1资料与方法

1.1资料选取选取盐城大丰地区1960～2014年月降水量时间序列数据作为基础研究资料，数据来源于中国气象科学数据共享服务网。

1.2分析方法经验模态分解EMD是 HHT 变换的主要步骤，作为HHT变换进行分析的前提条件。先通过自适应EMD求取IMF，然后利用Hilbert变换来求取各IMF的瞬时频率以及Hilbert谱。然而由于现实的信号存在诸如脉冲、噪声以及间断信号等其他无关因素的干扰，影响在求取包络过程中对于极值点选取的准确度，一旦所求取的包络不是真实信号的极值点时，就会导致筛选得到的IMF不准确，因而产生模态混叠现象，使IMF失去真实的物理意义[8]。

鉴于受原始信号频率因素影响和EMD自身算法的缺陷等原因产生的模态混叠问题，Huang等[7]曾提出使用中断检测法来处理该问题，但由于此方法需要人为检验而存在效率不高等缺点未被广泛使用。谭善文等[9]提出利用多分辨率思想来实现EMD的分解，但这种方法是以EMD的自适应为代价的。最后由Wu等[10]提出使用添加白噪声的EEMD算法来进行IMF分量筛选，该算法有效地解决了EMD方法中存在的模态混叠等问题，简单有效地改善了HHT变换的准确度问题。

EEMD实质上是一种引用白噪声对信号进行辅助处理的方法，即噪声辅助信号处理。与降噪技术相反，通过将白噪声加入原始信号，然后通过平滑技术去除脉冲干扰。对于时变信号x(t)，可以将x(t)表述为[7]：

X(t)=x(t)+ω(t)

(1)

式中，X(t)是加入白噪声后的信号；ω(t)是符合正态分布的小幅度高斯白噪声，服从下式：

(2)

式中,N为噪声总数；ε为高斯白噪声的幅值；εn为根据IMF重构的信号与原始信号之间的误差值。基于式(1)，利用EMD对X(t)进行分解，得到各阶IMF分量，则X(t)可以表示为：

(3)

向X(t)加入符合(2)式的随机白噪声并重复上述过程，可得：

(4)

由于高斯白噪声具有零均值特性，当进行多次平均后，噪声将相互抵消，并不会影响到原信号的变化情况，则原始信号所对应的IMF分量imfn(t)为：

(5)

因此，x(t)可以分解为：

(6)

式中，imfn(t)为各个IMF分量；rm(t)为趋势项，表征整个事件序列信号的变化趋势。对各个IMF分量分别进行Hilbert变换，可得到其频率、周期以及Hilbert谱H(w，t)。将H(w，t)对时间求积分，即可得Hilbert边际谱h(w)[8]：

(7)

式中，w为瞬时频率；t为时间；T为所研究序列的时间长度。

2结果与分析

图1　1960～2014年盐城大丰地区降水量时间序列变化图Fig.1　Changes of precipitation time series in Dafeng District of Yancheng City from 1960 to 2014

2.1盐城大丰地区降水变化特征由图1可见，盐城大丰地区最大年降水量为1 718.60mm，发生在1991年；最小年降水量为494.00mm，发生在1978年；55a年平均降水量为1 046.15mm。1960～2014年盐城大丰地区相邻年降水量变化波动较大，尤其是在1972、1987、1990年，相邻的两年之间降水量落差情况比较明显。从总体趋势来看，55 a间盐城大丰地区年降水量在波动中增加。对该地区55 a年降水量监测数据运用最小二乘法计算年降水量趋势定量变化一次线性方程为：y= 0.099 8x+847.93，其中气候倾向率为0.998 mm/10 a，表示盐城大丰地区从20世纪60年代开始，年降水量总体上呈递增趋势，增幅为1 mm/10 a。

2.2EEMD在盐城大丰地区降水时间序列值中的应用使用EEMD方法对盐城大丰地区降水量进行IMF分解时，在55 a年降水量值的序列中加入的高斯白噪声序列的参数如下：标准差为0.01，100组。同时，通过平滑技术将导致模态混叠的奇异点进行滤除，而后将得到55 a盐城大丰地区降水量时间序列分解后的真实固有模态函数IMFs。这里的IMFs可以分为4个固有模态函数和1个趋势项，说明55 a盐城大丰地区年降水量序列包含多个时间尺度变化特征(图2)。根据图2对分解得到的IMF各分量的变化情况进行分析，其中IMF1～IMF4的频率值逐渐减小，平均波动幅值逐渐变小，波动周期依次变长。另外，趋势分量r代表了整个时间序列上盐城大丰地区年降水量的变化情况，即减少-增大-减小的勺形变化趋势。而图1展示的是年降水量最小二乘法处理后的结果，代表的是以整个时间序列作为一个整体的降水趋势。图2是对图1更为细致的趋势分解，可以理解为将整个时间序列进行多次划分后展示的降水量趋势。因此，二者并不冲突。

图2　基于EEMD对1960～2014年盐城大丰地区年降水量变化IMFs分解Fig.2　IMFs decomposition of precipitation changes in Dafeng District of Yancheng City based on EEMD from 1960 to 2014

根据IMF平均周期计算公式[6]，利用EEMD法分解后求得IMFs中心频率与平均周期。由表1可知，IMF1～IMF4中心频率分别是0.470、0.160、0.055和0.036 Hz，表明盐城大丰地区年降水量的变化分别存在2.12、6.10 a的年际变化和18.30、27.50 a的年代际变化。

表1　各阶IMF分量统计特征值

根据年降水量数据进行HHT-Huang变换后得到的55 a盐城大丰地区年降水量的Hilbert谱图(图3)。由图3可以看出，色标中红色标记值最大，图中红色标记线主要存在的频率范围为0～0.05和0.45～0.50 Hz。虽然在频率0.05～0.45 Hz也存在标记值较大的标记线，但其范围分散，时间上跨度较大，因此不具有统计显著性。由在HHT-Huang变换的基础上频率对时间积分后得到的结果(图4)可见，信号幅值在整个0～0.50 Hz内有分布，但幅值较大的主要分布在0～0.05和0.45～0.50 Hz，也表明在这2个区间中信号的能量最大，振荡最为剧烈，是整个时间序列中频率显著性最强的2个范围。图4进一步验证了图3中的分析结果，相应地，55 a盐城大丰地区年降水量的变化周期主要存在于2个区间，即20.00 a以上的代际变化和2.00～2.20 a的年际变化。与55 a盐城大丰地区年降水量变化趋势分解图中IMF分量的波动情况(图2)对比，边际谱和IMF分量变化均能给出序列的周期变化，但由于边际谱更能反映各个频率在整个周期上的能量幅值强度，进而更容易辨别整个频率序列上主变化频率所在的范围。因此，综合以上分析，55 a盐城大丰地区年降水量的变化周期主要为2.12 a的年际周期和27.50 a的年代际周期。

图3　1960～2014年盐城大丰地区年降水量时间序列的Hilbert谱图Fig.3　Hilbert spectrogram of precipitation time series in Dafeng District of Yancheng City from 1960 to 2014

图4　1960～2014年盐城大丰地区年降水量时间序列的边际谱图Fig.4　Marginal spectrogram of precipitation time series in Dafeng District of Yancheng City from 1960 to 2014

2.3EEMD降水时间序列值研究结果误差分析EEMD是在EMD的基础上发展起来的，它继承了EMD的自适性，但为避免尺度混合问题，使最终分解的IMFs保持物理上的唯一性，由于从开始之初引入了白噪声，那么对于分解到的IMF判断显得十分重要，即判断信号分解后IMF分量是否真实有效。将基于EEMD变换后得到的IMF分量按照时间序列进行叠加，将叠加后的重构信号进行绘图得到盐城大丰地区55 a年降水量的边际谱图。从图5可看出，重构信号与原始信号几乎重叠。为进一步分析重构信号与原始信号的失真度，将重构信号与原始信号进行数学计算并引入误差分析后进行绘图，得到基于IMF重构波形与原始波形的对比误差变化(图6)，在整个时间序列上，对比误差百分比始终控制在±0.6%范围以内，处于10-3的量级范围之内，对原始信号的研究几乎没有影响。由此可见，基于EEMD的HHT法可以有效地重构原始降水量时间序列值数据变化曲线。

图5　基于IMF重构波形与原始波形的对比Fig.5　Comparison between reconstruction waveform and original waveform based on IMF

2.4EEMD降水时间序列值研究结果显著性检验分析从图7可看出，E点处于95%与99%置信区间，B点较C、D点更接近95%的置信区间，可信水平约为95%。因此，综上B点与E点的可信度最高，即IMF1分量和IMF4分量具有显

图6　基于IMF重构波形与原始波形的对比误差变化Fig.6　Changes of comparative error between reconstruction waveform and original waveform based on IMF

著性意义，这与根据IMF分量中心频率得到的结果一致。因此 IMF1 与 IMF4分量显著性较强，则其对应的盐城市大丰地区年降水量时间序列值的变化周期为主周期，具有统计意义。

注：A点为原始数据；B点为IMF1分量；C点为IMF2分量；D点为IMF3分量；E点为IMF4分量。Note：A was original data；B was MF1 component；C was IMF2 component；D was IMF3 component; E was IMF4 component. 图7　各阶IMF趋势的显著性检验Fig.7　Significance test of IMF trend in each order

3结论

该研究基于EEMD分解的HHT变换分析法对1960～2014年盐城市大丰地区年降水量时间序列值变化进行了系统的研究，得出以下结论：

(1)1960～2014年盐城大丰地区相邻年降水量变化波动较大，从总体趋势来看，55 a间盐城大丰地区年降水量时间序列值在波动中增加，其中气候倾向率为0.998 mm/10 a。

(2)通过EEMD可以将55 a盐城大丰地区年降水量时间序列数据分解成4个IMF分量和1个趋势分量。结合HHT变换边际谱，55 a盐城大丰地区年降水量时间序列值的变化主要有2个周期区间，即2.00～2.20 a的年际变化和20.00 a以上的代际变化，其间主周期主要为2.12 a的年际周期和27.50 a的年代际周期，并通过95%显著性检验。

(3)EEMD可以有效地重构盐城大丰地区原始降水量时间序列数据的变化曲线，与原始数据的对比误差百分比始终控制在±0.6%范围以内，由白噪声带来的干扰误差几乎可以完全忽略。

参考文献

[1] CAI W J，SANTOSO A，WANG G J，et al.Increased frequency of extreme Indian Ocean Dipole events due to greenhouse warming[J].Nature，2014，510(7504)：254-258.

[2] 王兵，李晓东.基于 EEMD 分解的欧洲温度序列的多尺度分析[J].北京大学学报(自然科学版)，2011，47(4)：627-635.

[3] 吴亚楠，李鸿雁，王玉新.第二松花江流域近75年降水量变化特征分析[J].安徽农业科学，2012，40(4)：2257-2259.

[4] 李春强，杜毅光，李保国.1965-2005年河北省降水量变化的小波分析[J].地理科学进展，2010，29(11)：1340-1344.

[5] 辛渝，陈洪武，张广兴，等.新疆年降水量的时空变化特征[J].高原气象，2008，27(5)：993-1003.

[6] 时世晨.EEMD时频分析方法研究和仿真系统设计[D].上海：华东师范大学，2011：36-45.

[7] HUANG N E，SHEN Z，LONG S R.A new view of nonlinear water waves：The Hilbert spectrum[J].Annual review of fluid mechanics，2003，31(1)：417-457.

[8] 陈则煌，张云峰，谢菲，等.EEMD在雷暴日趋势特征分析中的应用[J].热带地理，2015(4)：601-606.

[9] 谭善文，秦树人，汤宝平.Hilbert-Huang变换的滤波特性及其应用[J].重庆大学学报(自然科学版)，2004，27(2)：9-12.

[10] WU Z，HUANG N E.Ensemble empirical mode decomposition：A noise-assisted data analysis method[J].Advances in adaptive data analysis，2011，1(1)：1-41.

Multiscale Analysis of the Precipitation Time Series in Dafeng District of Yancheng City Based on EEMD

CHEN Chao1, LI Peng-fei2,3, CHEN Shang-yun1et a

(1. Dafeng District Meteorological Bureau, Yancheng, Jiangsu 224100; 2. Key Laboratory for Aerosol-Cloud-Precipitation of China Meteorological Administration, Nanjing, Jiangsu 210044; 3. College of Atmospheric Physics, Nanjing University of Information Science & Technology, Nanjing, Jiangsu 210044)

Abstract[Objective] To research the precipitation time series in Dafeng District of Yancheng City. [Method] Based on the data treatment method of ensemble empirical mode decomposition (EEMD), multi-scale characteristics research of precipitation data were researched by Hilbert transform. Significance and error tests on research results were carried out. [Result] The sequence of annual precipitation of Dafeng District in recent 55 years could be broken down into four IMF (Intrinsic Mode Function) components and a trend component. There were mainly two frequency regions, which were 0 - 0.05 Hz and 0.45- 0.50 Hz. The annual cycle was 2.12 years and the decadal cycle was 27.50 years during the main period. Both of them passed the 95% significance test. EEMD precipitation could effectively reconstruct the curve of original data. Compared with the original data, error percentage was always controlled within ± 0.6%. Then the interference error by the white noise could be almost completely ignored. [Conclusion] This research provides scientific basis for the rational development and fully utilization of hydrologic resources, flood control and drought relief.

Key wordsAnnual precipitation; Time series; EEMD; Hilbert transform; Significance test; Dafeng District in Yancheng City

基金项目国家自然科学基金项目(41175003)；江苏高校优势学科建设工程资助项目。

作者简介陈超(1987- )，男，江苏盐城人，工程师，在读硕士，从事雷电过电压防护技术与气象数据统计研究。

收稿日期2016-03-08

中图分类号S 161.6；P 448

文献标识码A

文章编号0517-6611(2016)10-209-03