基于极点对称模态分解的北京市降水特征分析

段志鹏，李继清

(华北电力大学可再生能源学院，北京 102206)

0　引　言

北京市是我国水资源紧缺最严重的城市之一。在全球气候变暖的背景下，北京市在21世纪以来，全年降水明显减少，更是加剧了北京市的缺水危机，因此研究北京市降水的变化规律以及未来的演变趋势，成为人们普遍关注的中心问题。目前已经有许多学者对北京市降水的变化规律及演变趋势等方面做了许多研究，朱龙腾[1]等采用经验模态分解(EMD)方法分析了北京市1951～2009年降水时间序列的多尺度振荡变化，分析结果表示短期内北京市的年降水量将在波动中不断减少；孙振华[2]等通过对北京市16个雨量站55年的降水资料进行统计计算和频率分析，探讨了北京市降水的时空分布特征及变化趋势，发现20世纪90年代之后北京市连续干旱，2008年之后年降水可能会进入一个相对较丰的时期；王文[3]等利用1961-2008年北京11个台站的气候观测资料，分析了过去48 a北京城区和郊区的气温和降水变化趋势，发现无论郊区和城区的年降水量都是呈减少趋势。上述研究成果从不同角度分析了北京市气温、降水的周期和变化特征，不仅有利于深入理解气候变化过程与规律，进行气候预测，更对北京市社会经济发展具有重要意义。

气候系统是非线性复杂系统，相关研究表明，气候系统除了具有非线性和非平稳性，还有层次性，许多大小不一的时空尺度构成了多层次结构[4]，传统趋势分析中，大尺度循环和趋势变化很可能混合在一起，不能分辨出趋势变化还是周期震荡。极点对称模态分解(ESMD)方法是近年来发展的局部自适应时间序列分析技术，是在经验模态分解方法基础上改进的一种新的非线性、非平稳时间序列分析方法。较之常用的小波分析方法，ESMD可以根据时间序列局部时变特征进行自适应地时频分解，完全摆脱了傅里叶变换的束缚，并可得到极高的时频分辨率，非常适合对非平稳、非线性时间序列进行分析。与广泛应用于气候水文要素分析的EMD方法相比，该方法解决了EMD方法中的“模态混叠”(或称“频率交叉”)问题，能够提供更好的尺度分解，模态分量可以准确体现时间序列的波动特征和趋势变化，在气候数据分析方面具有较大优势。

鉴于此本文采用ESMD分解对北京市1951-2015年的年降水序列进行多时间尺度和突变、趋势分析，并结合小波分析、Mann-Kendall 法和重标极差(R/S)分析法对ESMD方法的结果进行对比验证，为今后北京市降水预报和水资源规划管理提供参考。

1　研究方法

1.1　经验模态分解的基本原理

1998年，黄锷提出了经验模态分解(Empirical Mode Decomposition，EMD)方法，这种新的数据处理方法能够在时域和频域上提取数据的信息，很好的保留数据本身的特性，具有良好的自适应性，非常适合处理非线性、非平稳信号[5]。EMD方法是把不同周期的波动从原信号中分离出来，并且该波动是平稳的，最后得到的是趋势分量。不同尺度的波动被定义成为本征模态函数(Intrinsic Mode Function，IMF)。经验模态分解的计算过程如下：

(2)计算h1(t)=x(t)-m1(t)；

(3)将h1(t)作为新的信号序列，重复上述操作直到h1(t)成为一个零均值过程；

(4)得出零均值过程h1(t)后，将其作为第一个IMF分量c1(t)，它表示原始信号最高频的分量；

(5)从原始信号x(t)中减去c1(t)，得到去除一个高频分量的新信号过程r1(t)，即：

r1(t)=x(t)-c1(t)

(6)将r1(t)作为原始数据，重复(1)～(4)步骤，依次可得到c2(t)、c3(t)、c4(t)等，一直到满足迭代条件为止，则原始信号重构为：

(1)

式中：rn(t)为趋势项，代表信号的平均趋势。

1.2　极点对称模态分解原理

极点对称模态分解(Extreme-point Symmetric ModeDecomposition，ESMD)方法是经验模态分解方法的新发展，该方法解决了EMD方法中筛选终止和“模态混叠”(或称“频率交叉”)问题，并在气候数据分析方面具有较大优势[6]：①擅长寻找变化趋势，能够观测序列中分离出年际变化趋势和气候变化总趋势，有助于探究全球气候变暖问题；②擅长异常诊断，能够从分解模态中发现异常时段与频段，有利于气候异常研究；③擅长时-频分析，先进的直接插值法能够直观地分析各时间尺度上的频率变化。

ESMD方法是由两部分组成：第一部分是模态分解，可产生数个模态与一条最佳自适应全局均线；第二部分是时-频分析，利用直接插值法计算固有模态的瞬时频率，分析各时间尺度上的频率变化，判断各时间尺度下突变发生的时间。ESMD方法模态分解过程如下：

(1)找出序列X的所有极大值和极小值点，记作Ei(1≤i≤n)。

(2)用线段连接所有相邻的极点并将线段中点依次记为Fi(1≤i≤n-1)，并通过一定方式在左右两端添补边界中点F0和Fn。

(3)利用n+1个中点构建p条内插曲线L1，…，Lp(p≥1)，计算其均值曲线L*=(L1+…+Lp)/p。

(4)对X-L*序列重复上述步骤，直到|L*|≤ε(ε是允许误差)或者筛选次数达到预设的最大值K，得到第一个经验模M1。

(5)对剩余序列X-M1重复上述4个步骤，直到剩余序列R只剩于一定数量的极点，便可分别得到经验模M2，M3,…。

(6)在限定区间[Kmin,Kmax]内改变最大筛选次数K值，重复上述5个步骤。然后计算方差比率σ/σ0，并画出它随K的变化图，在图中找出σ/σ0最小值对应的K0，以K0作为限制条件再次重复上述5个步骤，最后剩余项R就是序列X的自适应全局均线。

经过分解，原始的时间序列X可表示为X=∑Mi+R，即利用ESMD将时间序列X分解成一系列经验模态和一个剩余变量。

在第二部分时-频分析中，考虑到原有的希尔伯特变换把原本离散的信号转换成解析函数处理，这样的处理方式不可避免的要受到各种数学理论的限制。面对离散数据，在分析过程中应当尊重其离散性特征，故此提出了针对数据的“直接插值法”，借此不但可以直观地体现各模态的振幅与频率的时变性，还可明确地获知总能量变化[7]。直接插值法基本思路如下：

(1)寻找极值点，计算两个相邻极大值点和相邻极小值点之间的时间差；

(2)将这些时间段视为局部周期赋给其中点，画出时间-周期对应点图；

(3)将这些局部周期值取倒数得到局部频率，再做三次样条插值得到光滑的时间-频率变化曲线。

2　应用实例

2.1　研究区概况

北京市地处华北平原西北端，海河流域中部，土地面积16 410 km2，山区面积约占总面积的61%。地理环境特殊，海拔高度变差大，山地和平原之间过渡急剧，界限清晰。北京市属于典型的北温带半湿润大陆性季风气候，夏季高温多雨，冬季寒冷干燥，降水季节分配不均。地形的复杂多变造成了北京地区天气状况的多样性。暴雨和强对流天气是北京地区夏季(6-8月)主要的灾害性天气[8,9]，导致北京地区旱涝灾害频繁。如图1所示为北京市1951-2015近65 a的降水过程。其中，年最大降水量为1 407 mm，发生在1957年；年最小降水量为267 mm，出现在1965年；多年平均降水量为589.75 mm。对北京市的年降水量进行线性拟合，可以看出从1951-2015年年平均降水量整体呈下降趋势，线性降水减少速率为3.47 mm/a。

图1　北京市1951-2015年降水时间序列Fig.1 Annual precipitation time series from 1951 to 2015 at Beijing station

本文研究所用年降雨资料来源于中国气象局国家气候中心提供的北京市1951-2015年月降雨资料，是经过该机构整编审查后发布的，故资料的可靠性基本可以得到保证。对资料的代表性进行分析，如图2所示，绘制北京市年降雨量模比系数差积曲线，可以看出随时间的增长，曲线变幅越来越小，年降雨系列在50年以上时，模比系数差积曲线稳定趋近于“1”，说明此系列资料具有一定的稳定性。

图2　北京市年降雨量模比系数差积曲线Fig.2 Beijing annual precipitation model ratio coefficient curve

2.2　降水变化的多尺度分析

ESMD方法具有自适应性和基于信号的局部变化特性，适用于非线性、非平稳序列的时频分析，与其他方法相比可以更好地提取非线性序列的变化趋势。因此利用ESMD方法对北京市1951-2015年年平均降水量的时间序列进行了分解，计算出当方差比率最小时最佳筛选次数为26次，此时得到的分解结果最好，为4个模态分量(模态1～4)以及一个趋势余量R，如图3所示。

图3　降水时间序列ESMD分解的模态分量及趋势项Fig.3 Modal components and trend of precipitation time series by ESMD

为验证分解结果的可靠性，用分解得到的模态1～4和趋势余量R进行合成，可得到一个重构序列，发现重构序列与原始的降水序列完全吻合，这说明ESMD方法分解结果是可信的。图2中ESMD分解的4个模态分量各自都具有独立的代表性，不会在同一时刻出现时间尺度相同的振荡，且每个模态分量一般都具有物理意义，它们各自反映了原序列中固有的不同特征尺度的振荡。为解析降水序列固有的不同特征尺度的振荡，利用周期图法[10]估算各个分解模态的平均周期，进而将每种尺度信号波动频率和振幅对原数据总体特征影响用方差贡献率表示出来，如表1所示。

表1　降水序列各分量的方差贡献率Tab.1 Contribution rates of ESMD decomposition for mean precipitation

表1给出了不同尺度模态分量所表征的不同时间尺度波动的平均周期和各模态分量的方差贡献率。由降水序列的各模态分量可以看出，在年际尺度上，北京市年平均降水以2.6 a(模态1)和4.3 a(模态2)的周期振荡为主；在年代际尺度上，则以14 a(模态3)和21.7 a(模态4)的周期振荡为主。

结合表1和图2可以看出：北京市65 a的降水序列包含多个时间尺度特征。其中模态1[图2(a)]表示的2.6 a周期方差贡献率最大，达到了32.47%，信号振荡非常明显，表明降水量变化呈现减少—增加—减少的循环变化，基本反映了20世纪70年代以前北京市降水丰枯交替的变化情况；模态2[图2(b)]表示的4.3 a周期方差贡献率约为16.74%，20世纪70年代以前振幅波动较大，而70年代以后振幅逐渐减小，趋于平稳，与模态1表现出相似特征；模态3[图2(c)]表示的14 a周期方差贡献率约为27.5%，仅次于模态1，振幅在20世纪60年代末期至20世纪90年代中期变化平稳，波动相对较小，而在20世纪50年代初期至60年代末期和21世纪之后，振幅明显增大，说明相比于其他时期，这两个时期的降水量变化幅度较大；模态4[图2(d)]表示的21.7年周期方差贡献率仅为1.37%，整个时间跨度上波动幅度比较稳定。对比模态3和模态4可以发现，模态3的振幅变化大而模态4的振幅变化显著变小，这表明降水异常主要发生在周期为准14 a的时间尺度上。趋势项R[图2(e)]为最佳自适应全局均线对应于年际降水变化，其方差贡献率为21.92%，反映了过去65 a降水序列的整体变化趋势，从图2(e)中可以看出，趋势项在整个时间尺度上呈现非线性的逐渐下降趋势。1951-1970年降水量逐年减少，1970-1990年降水量总体变化不大，1990-2015年降水量又逐年下降，且下降趋势较为明显。

小波分析是近年来十分热门的一个前沿领域，因其对信号处理具有特殊优势，已被广泛应用于气象和气候序列的时频结构分析中[11]。为了对比ESMD方法和小波分析在信号处理中的性能，采用水文气象上常用的Morlet小波变换对北京市的年降水距平序列进行分析，并且为了减弱序列开始点和结束点附近的边界效应，采用对称外延法将降水序列扩展[12]。如图4所示为北京市近65 a降水距平序列的小波系数实部时频分布图。该图反应了降水序列各种周期的强弱和分布，小波系数为正时，表示降水量相对偏多，值越大降水越多；小波系数为负时表示降水量偏少。从图3中可以看出北京市降水存在多时间尺度特征，降水量变化存在准2、5～8、10～15和准25 a的周期。整个时间跨度上正负相位交替出现，说明降水量经历丰枯交替变化。综上所述，表2给出了小波与ESMD分析周期结果对比，二者结果有一定差异，而且呈现出周期越长差异越大的规律。究其原因，小波分析是以傅里叶变换作为理论依据，存在一定的局限性，如小波基的选取、固定的基函数和恒定的多分辨率等问题[13]。与小波分析相比，ESMD方法摆脱了傅里叶变换理论的束缚，具有较强的灵活性和自适应性，依据数据自身特点进行分解，更有助于探索事物的内在规律。

图4　小波系数实部时频分布图Fig.4 The real part distribution of wavelet coefficients

表2　ESMD与小波分析周期结果对比　aTab.2 Comparison of cycle results by ESMD and Wavelet Analysis

2.3　不同时间尺度下降水突变分析

将分解得到的4个固有模态利用直接插值法绘出其瞬时频率与瞬时振幅的时变图，如图5所示，图5中F代表频率，A代表振幅。这样的分布图可以明确表达出模态的振荡强度与变化快慢程度。观察F1和A1可知，在2.6 a尺度下，突变点较多，降雨偏多和偏少交替相对比较频繁，如1957年时发生了突然的低频、大振幅振荡，对比原始数据可以发现，1957年北京市降水量仅487 mm，而1956年降水量达到了1 116 mm，下降了约56.4%；对比F2和A2，4.3 a尺度下，突变点在1961年和1995年；而且F1和F2存在一个频率恒定的时期，这是其他方法不可见的结果。F3和A3对比可知，在14 a尺度下，在1988年附近发生突变；观察F4和A4，发现在21.7 a的年代际尺度上，20世纪70年代以后的模态频率持续增高，而振幅却一直减小，说明在21.7年时间尺度上北京地区年降雨量持续减少。

图5　降水数据各模态频率(F)与振幅(A)的时变图Fig.5 Time-varying figure for each modal frequency and amplitude of precipitation data

Mann-Kendall 法(简称M-K) 是一种非参数统计检验方法，M-K法能很好地揭示时间序列的突变及其趋势变化特征[14]。为了验证北京市年降水量突变的显著性，作M-K检验曲线，可以确定突变发生在1990年附近，突变发生后北京市降水量呈现减少的趋势。对比两种方法的突变分析结果，不难看出ESMD方法具有更大的优势。M-K方法仅能找到整个时间序列的突变年份，而ESMD方法时频分析通过不同模态瞬时频率与瞬时振幅的对应关系图，可以找出不同时间尺度下的突变年份，分析结果更加详细全面，特别对局部时间序列的降水变化有较好分析，更好地为水文预报提供参考。

2.4　降水变化的趋势分析

图6为运用ESMD方法对降水序列进行不同时间尺度的分解后趋势项对原始数据拟合情况，表示降水序列整体的变化趋势。从图6中可以看出，在整个时间尺度上呈现非线性的逐渐下降趋势，预测到北京市降水量短期内呈现波动中下降的趋势。

图6　最佳自适应全局均线对降雨数据的拟合情况图Fig.6 The fitting of the best adaptive global mean curve to the rainfall data

R/S分析法最初由英国水文学家赫斯特(Hurst,1951年)在研究尼罗河水坝工程时提出，该方法认为降水量的变化是具有自相似性和长程依赖性的时间序列，Hurst指数H是定量描述时间序列自相似性与长程依赖性的有效方法[15,16]。利用R/S分析法侧面验证ESMD方法对于趋势判断的正确性。北京市年降水的Hurst指数H=0.649 38，大于0.5，说明年降水序列具有持续性，预测未来降水将持续减少。因此运用ESMD方法对降水序列的趋势判断具有可靠性，这种趋势变化特征可以为北京市水资源的合理利用和有效管理提供科学的依据。

3　结　论

ESMD是一种适用于非线性、非平稳序列的信号分析方法，将其应用于气候要素时间序列，可提取气候变化的固有时间尺度，不但能够从数年的观测序列中分离出年际和年代际变化趋势，还能够从其中分离出气候变化的总趋势。本文利用ESMD方法对北京市1951-2015年降水序列进行分解，可以得到以下结论：

(1)北京市65 a的降水序列包含多个时间尺度特征。在年际尺度上具有2.6 a和4.3 a的周期变化，在年代际尺度上具有14 a和21.7 a的周期变化，其中年际周期变化在北京市整体降水量变化中占据主导地位。

(2)不同时间尺度下突变发生的时间不同。在2.6 a尺度下，突变点较多，降雨偏多和偏少交替相对比较频繁；4.3 a尺度下，突变发生在1961年和1995年附近；14 a尺度下，在1988年附近发生突变。

(3)北京市降水量在研究时期内呈现“下降—平稳—下降”的非线性下降趋势，且后期下降较为明显，因此预测未来几年内北京市年降水量仍将持续减少。

在全球变暖背景下，自然因素和人类活动影响相互叠加，水文时间序列的变化规律更加复杂，本文通过对北京市降水序列的研究表明ESMD方法在水文序列的多时间尺度和突变、趋势分析上具有较好的效果。但是目前的研究仅处于初步阶段，如何确定各个模态分量的影响因子，计算其权重并最终应用到水文时间序列预测中，将是以后研究的重点工作。

□

参考文献：

[1]朱龙腾，陈远生，李璐，等．1951-2009年北京市降水变化情势分析[J]．水资源保护，2012，28(3)：42-46．

[2]孙振华，冯绍元，杨忠山，等．1950-2005年北京市降水特征初步分析[J]．灌溉排水学报，2007，(2)：12-16．

[3]王文，张薇，蔡晓军.近50 a来北京市气温和降水的变化[J]．干旱气象，2009，27(4)：350-353．

[4]薛春芳，侯威，赵俊虎，等．集合经验模态分解在区域降水变化多尺度分析及气候变化响应研究中的应用[J]．物理学报，2013，62(10)：109-203．

[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:903-995．

[6]王金良，李宗军．可用于气候数据分析的ESMD方法[J]．气候变化研究快报，2014，(3):1-5．

[7]王金良，李宗军．极点对称模态分解方法——数据分析与科学探索的新途径[M]．北京：高等教育出版社，2015：7-16．

[8]北京志编委．北京志地质矿产水利气象卷水利志[M]．北京：北京出版社，2000．

[9]北京市水利局．北京水旱灾害[M]．北京：中国水利水电出版社，1999．

[10]李强，吴健，许正文，等．利用EMD方法提取太阳活动周期成分[J]．空间科学学报，2007，27(1):1-6 ．

[11]魏凤英．现代气候统计诊断与预测技术[M]．2版.北京：气象出版社，2007:36-104．

[12]王文圣，丁晶，向红莲．小波分析在水文学中的应用研究及展望[J]．水科学进展，2002，13(4)：515-520．

[13]杨永峰，吴亚峰．经验模态分解在振动分析中的应用[M]．北京：国防工业出版社，2013:52-53．

[14]白晶，延军平，苏坤慧．1958- 2007 年秦岭南北气候变化的差异性分析[J]．陕西师范大学学报， 2010，38(6): 98-105．

[15]常远勇，侯西勇，毋亭，等．1998-2010年全球低纬度降水时空特征分析[J]．水科学进展，2012，(4)．

[16]江田汉，邓莲堂．Hurst指数估计中存在的若干问题：以在气候变化研究中的应用为例[J]．地理科学，2004，24(2)：177-182．