陈 超，李鹏飞，肖稳安，张春龙，吕东波

(1.江苏省盐城市大丰区气象局，江苏 盐城 224100；2.中国气象局 气溶胶与云降水重点开放实验室，江苏 南京 210044；3.黑龙江省气象灾害防御技术中心，黑龙江 哈尔滨 150030；4.中国气象科学研究院 高纬度雷电研究实验室，黑龙江 哈尔滨 150030)

基于EEMD与Markov Chain的雷暴日动态与趋势预测
——以盐城地区为例

陈 超1,2，李鹏飞2,3,4*，肖稳安2,4，张春龙2,3,4，吕东波3,4

(1.江苏省盐城市大丰区气象局，江苏 盐城 224100；2.中国气象局 气溶胶与云降水重点开放实验室，江苏 南京 210044；3.黑龙江省气象灾害防御技术中心，黑龙江 哈尔滨 150030；4.中国气象科学研究院 高纬度雷电研究实验室，黑龙江 哈尔滨 150030)

以盐城地区为例，通过滑动t突变检验、EEMD(集合经验模态分解，Ensemble Empirical Mode Decompositon)以及Markov Chain(马尔科夫链)统计模型，对盐城地区54 a年雷暴日数据进行动态与趋势预测研究，并对周期研究结果进行显著性和误差检验，结果表明：盐城地区54 a年雷暴日呈现减少的趋势，突变年为1966年，其54 a年雷暴日序列可以分解成4个IMF(Intrinsic Mode Function)分量和1个趋势分量。主要存在0～0.05、0.45～0.50 Hz两个频率区间，相应的其周期为2.22 a年际周期和27 a年代际周期(可信度均超过95%)。EEMD重构雷暴日数据与原始数据的误差百分比处于±0.8%范围内，其马尔科夫链长周期预测显示未来年雷暴日处于25～35 d之间的概率为45%，大于35 d的概率为32%。

盐城地区；年雷暴日；EEMD；Hilbert变换；马尔科夫链

近年来，全球温度升高以及气候变化，导致极端强对流天气时常发生，雷电灾害也愈发频繁。由于敏感性元件在电子电器等智能电气领域广泛应用，导致设备的防雷电过电压、雷电电磁脉冲等的能力降低，由雷电灾害带来的经济损失呈现指数增长。因而探究雷暴日发生动态周期等规律受到越来越多专家学者们的关注[1-5]。

虽然近10年国内闪电定位系统发展迅速，但是由于其监测系统装备时间较短，可供研究基础数据的时间尺度较小等原因，并不能深度挖掘雷电较长年代际演变的规律，而年雷暴日作为衡量一个地区一年中雷电放电的天数，从建国初期就有完整的记录，数据时间跨度较大以及累积量大。因此，以雷暴日作为研究对象进行动态与趋势预测的统计分析，具有非常高的实用价值[6]。在对气象数据的年际变化趋势、周期特征、极端事件等方面的研究过程中[7-8]，小波分析法作为探究周期特征的主要手段被广泛使用。但是，最新研究表明：小波分析的自适应性具有相当大的局限性[9]。而由Huang等专家提出的基于EEMD分解的HHT(Hilbert Huang Transform)模型，可以对非线性非平稳信号进行分析，而且其自身具有非常好的自适应性等优点，势必逐渐替代小波分析法，成为多个领域内进行多尺度周期分析统计的主流方法。同样，对合理分析雷暴日的多尺度周期规律，具有目前其他算法无法比拟的准确度。盐城地区位于江苏省东部、黄海之滨，是长江三角洲重要的区域性中心城市，也是著名的“中国麋鹿之乡”、“中国优秀旅游城市”。本文基于滑动t突变检验，集合经验模态分解-希尔伯特黄变换法以及马尔科夫链模型，对盐城地区1960～2013年共54年雷暴日序列进行动态与趋势预测研究，为盐城地区相关政府部门进一步制定防雷减灾决策提供有力的科学参考依据。

1 方法介绍

1.1 滑动t突变检验法

气候突变的主要表现是气候数据在时空上从一个统计特性到另一统计特性的急剧变化。即从一种相对稳定的时态(或者是一种稳定趋势持续变化)跳跃式地变化成另外一种相对稳定的时态的现象[10]。

(1)

1.2HHT-EEMD法

由于HHT-Huang的核心经验模式分解算法(EMD，EmpiricalModeDecomposition)存在一定的模态混叠问题[11-14]。即一个IMF(IntrinsicModeFunction)中包含差异极大的特征时间尺度，或者相近的特征时间尺度分布在不同的IMF中，导致相邻的IMF波形相互混叠而难以辨认。由此可能降低Hlibert变换后结果的准确度，间接影响了多尺度分析的结果，造成多尺度周期的误判。然而，在EMD改进的EEMD法中，核心算法是进行加入白噪声信号[15]，验证发现，经过添加白噪声后的分解算法可以有效解决了EMD算法中的模态混叠问题，并且白噪声在计算过程中将相互抵消，基本不会影响原始数据的分解。

可以将时变信号x(t)表述为：

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

(2)

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

(3)

式中：ε是白噪声幅值；N是噪声的总数；εn是重构的信号与原始信号对比后的计算误差值。根据式(2)，将EMD对X(t)进行分解，得到各阶IMF分量，则X(t)可以表示为：

(4)

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

(5)

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

(6)

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

(7)

式中：imfn(t)为各个IMF分量；rm(t)为趋势项，最后对各个IMF分量分别进行Hilbert变换[13]。

1.3 马尔科夫链法

马尔科夫链是一类特殊的马尔科夫过程，其状态与时间都是离散形式，其原理是基于变量的现在状态与动向，从而对将来的状态与动向进行预测分析。由于具有的马尔科夫性对历史数据的需求不多，因而可通过计算状态转移概率预测内部状态的变化，适用性较好[16]。

若随机过程{x(t),t∈T}对于任意整数和状态a1,a2,a3,…,an∈A条件概率满足：

(8)

其一步转移矩阵为：

(9)

k步转移概率可由方程p(k)=pk得到。

2 盐城地区雷暴日数据分析

2.1 资料选取

盐城地区1960～2013年逐月雷暴日数据来源于江苏省气象局资料中心。将逐月雷暴日资料整理成年雷暴日时间序列，其中最大年雷暴日为54d，发生在1964年；最小年雷暴日为13d，出现在2013年，多年平均雷暴日为33.3d。最早的初雷日发生在1996年，最晚的终雷日发生在1997年。初雷日与终雷日间天数最长为320d。

从图1可得，盐城地区1960～2013年的相邻年雷暴日变化波动较强烈，相邻的年份之间年雷暴日变化情况明显的主要有1966、1987、1998年。对54a年雷暴日数据运用(最小二乘法)计算得到一次线性方程为：y=-0.3248x+678.39，其中气候倾向率为-3.248d/10a，表示盐城地区从1960年开始，年雷暴日总体上呈现出递减趋势，每10a减少3.248d。

图1 盐城地区1960～2013年年雷暴日序列变化图

在图2滑动t突变检验图中，发现在1966年发生了一次较为显著的突变，显著水平为0.01。可以得出盐城地区54a年雷暴日在1966年发生一次显著水平超过99%的突变。

2.2 盐城地区雷暴日序列EEMD分析结果

对盐城地区54a年雷暴日数据进行EEMD法分解，所使用的白噪声的选用为100组标准差为0.01的序列。分析得到1个趋势项r和4个不同尺度的固有模态函数IMFs，可见原始信息包含了不同尺度的变化特征，如图3所示。

根据图3中基于EEMD对盐城地区54a年雷暴日变化趋势分解图可得，IMFn(n=1,2,3,4)波动频率与幅值逐渐减小。r代表的趋势分量如图1显示的一次函数趋势近似均为逐渐减少的总趋势。

HHT变换中边际谱表示的是信号的某一频率在各个时刻的幅值之和，这意味着可以比较准确地反映信号的实际频率成分[13]。因此为探究年雷暴日序列值的周期性变化规律，对年雷暴日序列值的变化求出其边际谱具有重要实用价值。由于IMF是一种极值点数、过零点数均相等对称包络线，为确定每个IMFn(n=1,2,3,4)的周期，对其极大值点数和数据总量进行统计，计算公式如下[7]：

(10)

根据公式(10)，求得的IMFn(n=1,2,3,4)平均周期见表1。

图2 盐城地区1960～2013年雷暴日滑动t突变图

图3 基于EEMD对盐城地区54 a年雷暴日变化趋势分解图

由表1所列，IMFn(n=1,2,3,4)中心频率分别是0.45、0.19、0.09、0.04 Hz，相应求得盐城地区年雷暴日存在2种尺度周期变化，分别是年际变化(5.4 a、2.22 a)和年代际变化(27 a、10.8 a)。

图4为经过EEMD分解后，进行希尔伯特黄变换法(HHT-Huang)变换后得到的54年雷暴日的Hilbert谱图(图4-a)、边际谱图(图4-b)。由图4-a可得，在频率范围为0～0.05 Hz和0.45～0.50 Hz之间出现较为均匀集中的红色标记点。在0.05～0.45 Hz之间也出现了部分标记点，但是标记点的分布显得较为分散。图4-b是边际谱，其可以显著展示各个频率段在整个时间尺度上的能量变化[17]，可见能量较高的幅值部分的分布同样是在0～0.05 Hz和0.45～0.50 Hz之间。因此，推断主要的周期范围存在于2.0～2.2 a的年际尺度和20年以上的代际尺度2个区间。综上分析，盐城地区54 a年雷暴日的变化周期主要为：2.22 a的年际周期和27 a的年代际周期。

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

图4 盐城地区54 a年雷暴日的Hilbert谱、边际谱图

3 研究结果误差分析

EEMD是EMD算法的拓展和延续，它具有EMD的自适性，同时也由于加入高斯白噪声解决的尺度混叠的问题。对分解后的IMFn(n=1,2,3,4)数据进行序列叠加重构，来判断本文分析结果的准确度以及与原始数据的误差百分比。由图5-a可知为重构信号与原始信号2条曲线几乎重叠，难以分辨。将重构波形与原始波形进行对比误差百分比分析，结果显示，对比误差的百分比在整个时间序列上始终保持在±0.8%范围内，即千分之一的量级。因此，对原始信号的变换失真度几乎没有影响。

图5 基于IMFn重构波形与原始波形的对比与误差图

4 研究结果显著性检验分析

图6为各阶IMF趋势的显著性检验图，它可以判别通过EEMD分解后得到的IMF分量是具有实际意义的分量还是仅为单纯的噪音分量。从图6可知，C、D点处于B点以下，可信水平低于95%。B点较C、D点几乎处于95%的置信区间，可信水平为95%。E点则处于95%～99%置信区间之间。因此B点与E点的可信度最高(95%)，即IMF1分量和IMF4分量更具有显著性统计意义，这与根据IMFs分量中心频率得到的结果一致。

5 年雷暴日长期预测

根据1960～2013年盐城地区年雷暴日数据，计算出其转移矩阵[18]，运用渐近性计算出未来若干年当地雷暴日分布概率，对盐城地区年雷暴日进行长期预测。

图6 各阶IMF趋势的显著性检验

首先将雷暴日划分为以下4个等级：第1级，雷暴日天数≤15 d，其编码为1；第2级，雷暴日天数15～25 d，编码为2；第3级，雷暴日天数25～35 d编码为3；第4级，雷暴日天数≥35 d，编码为4。通过对盐城地区1960～2013年的雷暴日数据等级划分，计算出盐城地区年雷暴日各等级转移规律，具体见表2。

表2 盐城地区年雷暴日各等级转移规律

由表2可得一步转移概率p：

令n=0代表1960年，n=1代表1961年，以此类推，令{X(n),n=0,1,2,3,…}表示盐城地区年雷暴日1960年以后各年雷暴日，显然为齐次马尔科夫链，由此计算出各年雷暴日预测等级：

因为转移概率具有渐近性质，即极限概率分布，所以当n趋向无穷大时，p(n)趋近于固定的一组概率，及其预测的概率分布趋向于一组固定的概率分布。

通过matlab编程求解出当n趋向于无穷时，p(n)概率分布(图7)。从图7中可以看出未来盐城地区雷暴日数主要分别在25～35d之间，其概率达到45%。其次是雷暴日数大于35d的天数段，分布概率达到32%。换而言之，盐城地区未来的雷暴日基本分布在25d以上。

图7 雷暴日数等级概率分别

6 结论

本文利用多种统计分析法对盐城地区54a年雷暴日时间序列值变化进行了系统的研究，并概括为以下3个结论：

(1)盐城地区54a年雷暴日变化波动较大，年雷暴日总体趋势在波动中减少，气候倾向率为-3.248d/10a。滑动t突变检验发现在1966年发生一次显著水平超过99%的突变。

(2)基于EEMD的盐城地区54年雷暴日数据分解成4个IMF分量和1个趋势分量。结合HHT谱变换，54a年雷暴日的变化主要有2.22a的年际周期和27a的年代际周期，并通过95%显著性检验。重构变化曲线与原始数据对比误差百分比数量级仅为0.1%，高斯白噪声带来的干扰误差几乎可以完全忽略。

(3)基于马尔科夫链预测未来盐城地区年雷暴日数主要在25～35d内，其分布概率达到45%。其次是年雷暴日数大于35d的天数段，分布概率为32%。

由于闪电资料数据时间尺度较短，以及近3年雷暴日的资料数据还有待进一步整理。因此，文中在对年雷暴日的预测验证分析和与闪电资料数据进行对比分析方面会以此立足，进一步深入研究。

[1] 央美,旦增格列,达瓦泽仁,等.那曲地区雷暴天气时空变化特征及影响因素[J].应用气象学报,2014(6):751-760.

[2] 刘正源,姜苏,曹洪亮.呼和浩特市雷暴气候特征及其变化分析[J].江西农业学报,2013，25(2):79-82.

[3] 王纪军,郭红晨,卢广建.河南省雷暴日数时空分布的非均一性特征[J].自然灾害学报,2009,18(4):115-119.

[4] 左迎芝，孙健，韩通，等.基于两种观测资料的雁江区闪电活动规律分析[J].江西农业学报,2016,28(7):116-120.

[5] 王学良,张科杰,张义军,等.雷电定位系统与人工观测雷暴日数统计比较[J].应用气象学报,2014(6):741-750.

[6] 高燚,周方聪,劳小青.1999～2011年海南岛雷电灾害特征分析[J].自然灾害学报,2014，23(5):253-262.

[7] 王兵,胡娅敏,杜尧东,等.小波分析和EEMD方法在广州气温及降水的多尺度分析中的差异分析[J].热带气象学报,2014,30(4):769-776.

[8] 冯民学,周宇,虞敏,等.2013年江苏省ADTD闪电定位系统资料的评估分析[J].科学技术与工程,2015,15(7):79-84.

[9] 朱宁辉,白晓民,董伟杰.基于EEMD的谐波检测方法[J].中国电机工程学报,2013,33(7):92-98.

[10] 刘宇峰,原志华,吴林,等.基于湿润指数的近55年安康地表干湿变化趋势及周期研究[J].自然灾害学报,2016,25(3):11-21.

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

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

[13] 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.

[14] 朱正.Hilbert-Huang变换及其在目标方位估计和水声通信中的应用研究[D].哈尔滨:哈尔滨工程大学,2013:18-40.

[15] Wu Z, Huang N E. Ensemble empirical mode decomposition: a noise assisted data analysis method[J]. Advances in Adaptive Data Analysis, 2009, 1(1): 1-41.

[16] Putze A, Derome L, Donato F, et al. A markov chain monte carlo technique to sample transport and source parameters of Galactic cosmic rays[J]. Memorie Della Societa Astronomica Italiana, 2015, 462(1): 863.

[17] 谭善文.多分辨希尔伯特-黄(Hilbert-Huang)变换方法的研究[D].重庆:重庆大学,2001:15-25.

[18] Yan Y Q, Zhang D F. Markov chain model of navigation based on frequence[J]. Application Research of Computers, 2007, 24(3): 41-43.

(责任编辑：曾小军)

Dynamic and Trend Prediction of Thunderstorm Days in Yancheng Area Based on EEMD and Markov Chain

CHEN Chao1,2, LI Peng-fei2,3,4*, XIAO Wen-an2,4, ZHANG Chun-long2,3,4, LV Dong-bo3,4

(1. Meteorological Bureau of Dafeng District, Yancheng City, Jiangsu Province, Yancheng 224100, China; 2. Key Laboratory of Aerosol and Cloud Precipitation, Chinese Meteorological Administration, Nanjing 210044, China; 3. Meteorological Disaster Prevention Technology Center of Heilongjiang Province, Harbin 150030, China; 4. Laboratory of High-latitude Lightning, Chinese Academy of Meteorological Sciences, Harbin 150030, China)

Through adopting the movingt-mutation test, EEMD (Ensemble Empirical Mode Decomposition) and Markov chain statistic model, the author studied the dynamic and trend prediction of annual thunderstorm day number in the past 54 years in Yancheng area, and conducted the significance and error test for the periodicity research results. The results indicated that the annual thunderstorm day number in the past 54 years in Yancheng area revealed a decreasing trend, and the mutation year was 1966. The sequence of annual thunderstorm day number in 54 years could be decomposed into fourIMF(Intrinsic Mode Function) components and a trend component. There existed two main frequency scopes: 0～0.05 Hz and 0.45～0.50 Hz, and their corresponding period was 2.22 years and 27 years (their reliability all exceeded 95%). The error percentage between the reconstructed thunderstorm day data by EEMD and the original thunderstorm day data ranged from -0.8% to 0.8%. The long-period prediction by Markov chain shows that the probability which the number of annual thunderstorm days in the future is 25～35 d will be about 45%, and the probability which the number of annual thunderstorm days in the future is over 35 d will be about 32%.

Yancheng area; Annual thunderstorm day; EEMD; Hilbert transform; Markov chain

2017-02-28

国家自然科学基金资助项目(41175003)；江苏高校优势学科建设工程资助项目(PAPD)；黑龙江省气象局青年英才项目。

陈超(1987─)，男，江苏盐城人，硕士研究生,研究方向：气象数据挖掘统计与雷电防护技术。*通讯作者：李鹏飞。

P44

A

1001-8581(2017)06-0094-06