关于代码无法运行情况，请下载完整程序包尝试运行。
链接：esmd程序包
提取码：zzl8
引言
上世纪六十年代，著名气象学家洛伦兹提出了提出了著名的“蝴蝶效应”，这种混沌学概念一经发表，便引起广泛关注，成为当今最富盛名的科学概念之一。
混沌——是指发生在确定性系统中的貌似随机的不规则运动，一个确定性理论描述的系统，其行为却表现为不确定性一不可重复、不可预测，这被证明为非线性动力系统中固有的属性。
气候系统，作为一个典型的混沌、非线性系统，其研究与预测有着许多的不确定性。对于混沌系统的研究，我们试图在无序中找出有序，对应的，在气候中，我们需要找出气候变化的趋势与跃变。
气候系统的特点，除了非线性、非平稳性外，还具有层次性。这意味着，一个地区的气候变化，往往是由多个层次的变化影响叠加的。换句话说，一个地区气候的变化，是由气候周期、大尺度环流、小尺度人为影响等多种因素决定。
同样的气象数据，也是一组非线性、非平稳的时间序列，且导致变化的物理机制并不明确。对于这类随机数据，我们需要将其分解成不同频率模态，从中寻找可能的变化规律，分别探寻不同模态的物理机制。
在模态分解算法，常用波谱分析方法，并衍生出三代信号处理技术。
第一代：Fourier谱分析：事先给定基函数（三角函数）
第二代：Wavelet分析：事先给定基函数（Morlet小波、墨西哥帽小波）
第三代：HHT方法：基函数来自数据本身（let data speak for themselves
我们即将要介绍的EMD与ESMD方法，就是基于第三代HHT所发展的。
EMD与其局限性
Hibert-Huang变换(EMD)，是一种数据自适应处理方法，不需要预先取定基函数或窗口长度，其分解模态不但频率可变振幅也可变，适用于非线性非平稳信号。
所谓数据自适应是指：依据数据自身的时间尺度特征来进行信号分解，无须预先设定任何基函数。这一点与建立在先验性的谐波基函数和小波基函数上的傅里叶分解与小波分解方法具有本质性的差别。因此，理论上它可以适用于任何任何类型的信号分解。
即：任何复杂数据都可用EMD方法分解为有限个周期振荡分量/本征函数（IMF: Intrinsic Mode Function）;它经过Hilbert变换后可以产生瞬时频（instantaneous
frequencies)。
在物理上，如果瞬时频率有意义，那么函数必须是对称的，局部均值为零，并且具有相同的过零点和极值点数目。在此基础上，NordneE.Huang等人提出了本征模函数（Intrinsic Mode Function，简称IMF）的概念。本征模函数任意一点的瞬时频率都是有意义的。Huang等人认为任何信号都是由若干本征莫函数组成，任何时候，一个信号都可以包含若干个本征模函数，如果本征模函数之间相互重叠，便形成复合信号。EMD分解的目的就是为了获取本征模函数，然后再对各本征模函数进行Hilbert变换，得到Hilbert谱。（EMD方法基本知识）
本征函数（IMF）满足两个条件：
1、data中的极值个数和0点个数必须相等或最多相差一个；
2、在任何一点，局部极大值和极小值组成的两
条包络线的均值为0。(Huang et al. 1998)
听起来非常抽象，简单解释一下：
首先，包络线（Envelope）是几何学里的概念，代表一条曲线与某个曲线族中的每条线都有至少一点相切 （曲线族即一些曲线的无穷集 ，它们有一些特定的关系。）
在EMD算法中，上下包络线即为将各极值点用三次样条插值(cubic spline line)连接形成，连接极大值点构成上包络线(upper envelope)，连接极小值点构成下包络线(lower envelope)。
首先，绘制出数据的上下包络线
根据上下包络线，绘制出均值包络线
原始信号减均值包络线，可以得到中间信号
判断得到的信号是否满足IMF的两个条件，满足，则接受，将其作为一个IMF，不满足，则重复上述步骤迭代，使用上述方法得到第一个IMF后，用原始信号减IMF1，作为新的原始信号，再通过1）~4）的分析，可以得到IMF2，以此类推，完成EMD分解。
上述过程称为筛选(shifting)
（参考EMD基础理论）
EMD的停止准则：
停止准则决定了产生一个本征模函数(IMF)筛选过程执行的数目，请参考浅谈EMD
同样的EMD也有着一些局限性，包括：
- 筛选次数难以确定
- 相邻模态混叠现象
- 结果不稳定
为了解决EMD中的模态混叠(mode mixing)混叠问题，EMD方法的提出者又建立集成经验模态分解(Ensemble empirical mode decomposition, EEMD)方法，其基本想法是：将原信号多次添加白噪声进行 EMD 分解再取平均。模态由于是多次加噪平均的结果，对扰动不敏感，被认为是分解稳定的并以此来理解物理现象。然而，直接添加噪声往往会造成不利结果，甚至会污染信号。对此，对EMD又有了新的发展——ESMD，由青岛理工大学理学院王金良与李宗军两位老师合作研发。
ESMD原理与流程
ESMD借鉴了 EMD 的思想，将外部包络线插值改为内部极点对称插值。可分为两部分
- 模态分解（ESMD)，得到模态
- 时频分析，得到各模态随时间变化特征
ESMD的算法流程如下（本人随手画的）：
结果与EMD类似，将原序列分解为经验模态和一个剩余变量（X=ΣMi+R）
ESMD中的停止准则比之EMD更加简单，借用“最小二
乘”的想法来优化最后剩余模态使其成为整个数据的
“自适应全局均线”，并由此来确定最佳筛选次数。当筛选次数达到最佳值时，停止迭代筛选。
时-频分析的流程如下：
抛弃传统积分变换概念，将Hibert变化改为使用针对数据的“直接插值法”，直观体现各模态振幅和频率的时变性。步骤如下：
- 寻找极值点，计算两个相邻极大值点和相邻极小值点之间的时间差;
- 将这些时间段视为局部周期赋给其中点，画出时间－周期对应点图;
- 将这些局部周期值取倒数得到局部频率，再做三次样条插值得到光滑的时间－频率变化曲线。
包括 Hilbert 变换在内的所有积分变换在分析时–频变化方面都存在固有缺陷，而ESMD抛弃了频谱分析依靠积分变换的传统观念，创造性地提出了针对数据的“直接插值(DI)法”。
关于Hibert变换的局限性，理论比较复杂，这里选用ESMD方法提出者书里的一句话说明：
In fact, no matter how the integral transform is defined, it is actually a uniform running-mean processing. Now that the processing is done on thedata, why not calculate the instantaneous frequency from the data in a directway?
简单来说，目前对于频率变换的方法基本为计算平均周期(频率)的基本过零方法（ zero-crossing method），但无法很好的瞬时变化，即，时间分辨率较低。
为此，ESMD提出了direct interpolating（直接插值法）方法来计算瞬时频率。
综上，从理论分析而言，ESMD方法在寻找变化趋势、异常诊断、时-频分析与能量变化方面都有着较好的优势。
ESMD应用
ESMD方法比较与EEMD法而言应用较少，但应用前景依然广阔，可运用于气候、海洋等领域中，下面简单举几个例子：
海-气通量的研究
海气通量的观测是气象观测中的重要部分，在边界层中，湍流可雷诺分解为：
当稳定时，可表示为：
然而，大多数时间观测到的为非平稳信号，此时波动会包含非湍流成分
此时：
红圈部分即为我们要提取的湍流通量，尝试用ESMD方法分解。
前五个模态即为我们需要的湍流成分，判断需要基于经验与物理机制（−5/3 power law）。
模态1-5合成后与原信号傅里叶谱图对比：
降水特征分析
数据：1951-2015年北京月降水资料，变化如下：
不同模态对应方差与周期：
频谱分析可用来观察各模态的突变点。
径流预测
径流预测通常用coupled decomposition-prediction的组合模型来预测，即，使用模态分解＋预测的组合方法，流程如下：
使用各方法对准确性进行了对比，可以发现ESMD的误差最小：
代码
disp('It is running. Please wait ...');
clear all;
javaaddpath(pwd);
esmd= esmd4j.Esmd();
% load data
data=load('winddata.txt');
Y=data(:,1);
delt_t = 0.05; % sampling period
t=esmd.init_t(length(Y),delt_t);
%parameters setting:
minLoop=1;
maxLoop=40;
extremeNumR=4; % >=4
jianGeNum = 1; %
rList = esmd.getVarianceRatio(t, Y, minLoop, maxLoop, extremeNumR);
[minVar,idx]=min(rList);
optLoop=idx+minLoop-1;
disp('optimal loop is: ');
fprintf('optLoop=%d\n',optLoop);
%-----------------------------------
x=minLoop:maxLoop;
figure(1)
plot(x,rList)
%parameters setting:
esmd.getSift(t, Y, optLoop, extremeNumR);
%-------------------------------------------
len = esmd.yImfsR.size();
figure(2)
for i=1:len
subplot(len,1,i);
plot(t, esmd.yImfsR.get(i-1));
%xlabel('t');
if i==1
ylabel('Y');
elseif i==len
ylabel('R');
else
ylabel(strcat('Imf',num2str(i-1)));
end
end
%-------------------------------------------
len=esmd.interfs.size();
figure(3)
%frequency distribution figure
for i=1:len
plot(t,esmd.interfs.get(i-1))
hold on
end
%-------------------------------------------
figure(4)
%frequency Amplitude figure
for i=1:len
subplot(2*len,1,2*i-1)
plot(t,esmd.interfs.get(i-1))
ylabel(strcat('F',num2str(i)))
subplot(2*len,1,2*i)
plot(t,esmd.upperEvelops.get(i-1))
ylabel(strcat('A',num2str(i)))
end
%-------------------------------------------
figure(5)
%AdaptGlobalMeanCurve on Y
len = esmd.yImfsR.size();
plot(t,Y);
hold on;
R=esmd.yImfsR.get(len-1);
plot(t,R);
%-------------------------------------------
figure(6)
%plot Y-R
plot(t,Y-R)
%-------------------------------------------
figure(7)
plot(t,esmd.energy)
%clear all;
%javarmpath(pwd);
先运行代码1，根据输出图设置optLoop，或采取默认值，再运行代码2。
总结
总体而言，ESMD方法是EMD方法的一种新发展，在理论山解决了一些EMD方法的局限性，比之更为常用的EEMD方法具有不同的优势。
然而此方法较为冷门，目前应用较少，可作为一种新的方法来了解尝试。
笔者认为，统计与算法，永远只是工具，沉迷于理论的研究并不是气象科研人员的工作。以应用为主，在应用过程中，探索物理机制与原理，才是气象工作者的目的。
参考文献
- J. L. Wang and Z. J. Li., “Extreme-Point Symmetric Mode Decomposition Method for Data Processing,” Ad_x0002_vances in Adaptive DataAnalysis, Vol.5,2013.
- Li, H.-F., Wang, J.-L. and Li, Z.-J. (2013) Application of ESMD method to air-sea flux investigation. International Journal of Geosciences, 4, 8-11.
- Yi-zhen Li, Chun-fang Yue; Prediction and analysis of non-stationary runoff extreme sequence based on ESMD combination prediction model. Water Supply 1 June 2020; 20 (4): 1439–1452.
- 段志鹏,李继清.基于极点对称模态分解的北京市降水特征分析[J].中国农村水利水电,2018(03):17-21.
2371
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
