张恒璟 ，程鹏飞

(1. 辽宁工程技术大学 测绘学院，辽宁 阜新 123000； 2. 中国测绘科学研究院，北京 100830；3. 国家测绘产品质量检验测试中心，北京 100830)

一、引 言

1998年黄鄂教授提出了研究非线性非平稳信号分析技术的希尔伯特-黄变换方法[1](Hilbert-Huang transformation，HHT)，它包括经验模式分解(empirical mode decomposition，EMD)和希尔伯特变换两个步骤，已经被广泛应用于爆破振动信号、机械振动信号、桥梁振动信号、阵列声波信号、地震波信号、电力谐波信号、语音信号、金融时间序列信号等社会生产建设的诸多领域。

EMD具有类似于小波分解的滤波性质，国际上已有不少关于EMD滤波去噪算法的研究结果，但EMD方法本身还不能用完整的数学公式加以描述或定义，相比较小波完善的理论体系，EMD技术发展并完善的路还很长，且需要更多的人去关注和研究。本文以EMD滤波去噪问题为切入点，总结了近年来专家学者提出的基于EMD的滤波去噪算法，分析了各自算法的适用条件及可能的改进，最后以GPS高程时间序列噪声提取的具体问题，分析了EMD去噪算法相对于传统周期拟合模型获取噪声的优势。

二、经验模式分解及其滤波性质

经验模式分解的基本方法是：通过不断的“筛选”，将信号分解为一系列频率由高到低的本征模态函数分量(intrinsic mode function，IMF)和残差项，见式(1)，具体的过程可参考文献[1]。

(1)

式中，N为分解得到IMF的个数；rN(t)为分解的残差项；ci(t)表示第i个IMF分量。

分解产生的IMF函数需要满足两个条件：一是待分解信号中极值点个数与过零点个数相等或至多相差1；二是在任意点上，由局部极大值定义的上包络线与局部极小值定义的下包络线的均值为零。这种分解基于数据自身的局部特征，与傅里叶频谱分析分解成恒定振幅与频率的一系列正弦函数和小波分解需预先给定基函数相比，EMD是一种自适应的信号分解方法。

由于受信号中断的影响，EMD会产生模式混叠现象[2]，不仅混淆了信号的时频分布特性，更重要的是使信号IMF分量的物理意义变得不清晰。为了减弱模式混叠现象的影响，Wu Zhaohua等人在EMD方法的基础上提出了整体经验模式分解方法[3](enable EME，EEMD)。这是一种噪声辅助的数据分析方法，其基本思路是在分解之前，将原始序列信号加入一定信噪比的白噪声，生成待分解信号，经过EMD分解后生成一系列的IMF分量；由于每次随机生成的白噪声都是不同的，这样反复EMD分解过程，得到多次分解的IMF分量，将对应相同分解尺度的IMF分量取平均，即得到了最终信号分解的IMF分量。

Flandrin研究了分形布朗运动的增量过程，即分形高斯噪声的EMD分解，获得了一系列高低频IMF分量，试验发现EMD具有类似小波分解的二进滤波性质[4-6]；Wu Zhaohua采用EMD分解高斯白噪声，通过数值试验发现EMD等效于一个二进滤波器，且分解的IMF分量都服从正态分布，在以周期对数为横轴的坐标系中，IMF分量的傅里叶频谱曲线几乎相同[7]；Wu Zhaohua又研究了当EMD分解的筛选次数改变时的滤波性质[8]，以白噪声和delta函数为例，发现当筛选次数从有限次到无限次时，滤波的比率从二进滤波的2逐渐变化为1，但实际上筛选次数不可能无限多，因此比率总是大于1。经验模式分解类似于小波的滤波性质，为其在信号去噪方面的应用提供了有利的条件。

三、基于EMD的滤波去噪算法

1. 软阈值法消噪

Donoho等设计了基于小波的软阈值去噪算法，分离了信号中混合的加性高斯白噪声[9-10]。基于对EMD滤波性质的试验结果，类似于小波分解的滤波去噪功能，Boudraa于2005年提出了基于EMD的软阈值去噪算法[11]，基本原理是：将原始信号x(t)经过EMD分解后，计算每一个IMF分量消噪的软阈值，分别进行去噪。

首先计算每一个IMF序列的绝对中位差(median absolute deviation，MAD)

MADj=median{|cj(t)-median(cj(t))|}

(2)

式中，median()表示取序列元素从小到大排列后的中位值，如果序列长度为奇数，则取中间的元素；如果序列长度为偶数，则取中间两个元素的算术平均值。

然后计算IMF所含的噪声水平(即噪声中误差)

(3)

再计算每一个IMF消噪的软阈值

(4)

式中，L是信号长度。

最后消除IMF分量的噪声

(5)

软阈值法消噪的公式中，下标j表示第j个IMF分量，j=1，2，…，N，N是EMD分解得到的IMF分量的个数。

软阈值法消噪的滤波器设计是基于分解信号中包含有高斯白噪声的前提，去噪性能比较优越；如果信号中包含有色噪声的干扰，则设计的滤波器能否有效去噪，还需要进一步试验研究。

2. 部分IMF重构算法

(6)

此处k如果取1，则去噪信号与原始信号相同，表明原始信号中不包含噪声。基于EMD分解的连续均方误差CMSE准则如下

(7)

由式(7)可知，CMSE值相当于第k个IMF分量的能量密度，从数学公式上称为该IMF分量的均方误差(mean squared error，MSE)。当k取2～N时可依次计算出CMSE值，将CMSE值第一次显著发生变化处的k值(记为临界点js)作为干净信号的起点，见式(8)，重构后的信号表示见式(9)，临界点以前的IMF为信号中包含的噪声，即从原始信号中扣除重构信号得到了部分IMF重构算法的噪声。

(8)

(9)

部分IMF重构算法在判断临界点时，依据的是CMSE值第一次显著发生变化的位置，这一显著性是通过试验观察得到的，缺少数学规则和证明。EMD分解得到一系列频率由高到低的IMF分量，每一个IMF的CMSE值如果比较接近，就带来了判断临界点的困难。单纯依靠CMSE准则的判断就需要改进，此时可以考虑信号的实际物理意义，根据信号EMD分解后每一个IMF分量的时频特性，并结合其他的准则(如平均周期)对部分IMF重构算法作进一步完善。

3. 最优信号重构算法

Weng在2008年提出基于EMD的最优信号重构算法[14](optimal signal reconstruction，OSR)。具体做法是：将原始不含噪声的信号d(t)加入一定量的高斯白噪声w(t)后，得到含噪信号x(t)；将含噪信号采用EMD分解成一系列IMF分量和残差项。假设对IMF分量线性加权，得到去噪后的重构信号如下

(10)

式中，ai表示第i个IMF分量的权重系数。

(11)

J1取得最小值时，令上式对ai求导结果等于0，计算出每一个IMF对应的权重系数。

后半部分原始信号用来测试去噪效果，将后半部分信号进行EMD分解，将权重系数回代到式(10)中，通过计算后半部分重构信号与后半部分原始信号的MSE来衡量去噪性能。

Weng还提出了双向最优信号重构算法(bidirectional OSR，BOSE)，对不同层次IMF分量加权的同时考虑每个IMF分量中信号之间的水平相关性(即IMF分量各个采样值之间的相关性)。从Weng采用仿真信号的试验结果来看，高频IMF分量的权重系数较小，低频权重系数较大，可知噪声主要包含在信号的高频部分。

OSR方法的局限性是必须知道原始干净信号，这样才能基于MSE最小准则估计每个IMF分量的权重系数。但一般情况下试验得到的原始信号已经包含了噪声，原始干净信号无法获得，因此基于EMD方法的OSR去噪算法在应用中就受到了限制。

4. GPS高程时间序列的滤波去噪

GPS高程时间序列噪声分析时，首先需要建立合适的序列拟合模型获取残差序列，将拟合残差作为GPS高程序列包含的噪声。统一的GPS高程时间序列周期拟合模型如下

y(ti)=a+bti+csin(2πti)+dcos(2πti)+

Tkj)ti/τj)H(ti-Tkj)+vi

(12)

式中，ti(i=1，2，…，N)为以年为单位的时间，在此以每日解构成坐标时间序列；待求系数a为序列的平均值；b为线性速率；c、d和e、f分别为年周期和半年周期项的系数；gj则表示由于各种原因引起的阶跃式坐标突变(如远场大地震引起的同震位移、由于仪器或天线变更引起的位移、甚至由于某种有清楚原因引起的点位变化等)；hi和kj表示更加复杂的阶跃变化改正项；最后一项vi为GPS高程时间序列周期函数拟合模型的残差，即下一步分析的序列观测噪声。

国内外学者采用式(12)的截断模型或完整模型，获得了GPS高程时间序列的噪声，并以此为基础，采用极大似然估计方法定量分析噪声类型和含量[15-22]。笔者也曾分别采用抗差功率谱估计和EMD分解的方法，初步分析了GPS高程时间序列的周期特性[23-24]，功率谱估计结果表明序列中并不存在精确的半年、年或两年周期项；EMD分解后各个IMF分量的周期也不是一个恒定的数值，而是一个明显的时变序列，说明了GPS高程时间序列信号的非平稳性特征，这与传统的研究结果一致。

笔者在文献[24]中研究了EMD分解方法在GPS高程时间序列分析中的初步应用。GPS高程时间序列信号经过EMD分解为一系列频率由高到低的本征函数模态分量和残差项，通过残差项与邻近的IMF分量合并，并结合各个IMF方差贡献率指标，识别GPS高程时间序列信号的非线性趋势项并分离。中低频的IMF分量是序列明显的季节性变化项，包括明显的近似月、双月、三个月、半年、年、双周年信号，但并不存在严格精确的周期。因此采用周期拟合模型获得的序列残差中仍包含着周期成分，以此作为噪声的分析结果就存在着偏差。

基于EMD的软阈值去噪算法，是假设信号中包含加性高斯白噪声条件下提出的，GPS高程序列信号具有非平稳性特征，序列中既包含白噪声也包含有色噪声的影响，这已经被大量研究所证实，因此软阈值去噪算法应用于GPS高程序列信号去噪，理论上具有局限性，去噪效果需要进一步的试验验证。

基于EMD的部分IMF重构算法，通过构建一个低通滤波器，基于CMSE准则，识别出信号中包含的噪声。通过笔者在文献[24]中的试验可知，GPS高程序列信号的低频部分为序列的非线性趋势项，中低频的IMF分量是序列明显的季节性周期信号，高频的IMF分量既包含有序列噪声，还有在目前研究的序列长度区间上暂时无法识别的短期周期项，也可以一并作为噪声处理。这样构建的低通滤波器就为滤波GPS序列噪声提供了方法上的可行性。

四、结束语

采用周期函数模型拟合GPS高程时间序列时采用固定的半年和年周期的做法，与序列实际的周期运动情形不符，传统功率谱技术的谱分辨率低，小波分解的基底固定的也违背了序列本身的特征，因此传统的法拟合序列周期项进而获得序列噪声的做法存在着缺陷，似乎没有一个最佳有效的解决方案。

本文将基于EMD分解的信号去噪算法，引入到GPS高程时间序列分析领域，是一个新的尝试，为GPS高程时间序列信号的周期特性和噪声研究开辟了一个新的途径。

参考文献：

[1] HUANG N E，WU M L，LONG S R，et al. The Empirical Mode Decomposition and the Hilbert Spectrum for Nonlinear and Non-stationary Time Series Analysis[J].Proceedings of the Royal Society，1998，454(1971)：903-995.

[2] HUANG N E，SHEN Z，LONG S R. A New View of Nonlinear Water Waves：The Hilbert Spectrum[J]. Annual Review of Fluid Mechanics，1999，31：417-457.

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

[4] FLANDRIN P，RILLING G，GONALVES P. Empirical Mode Decomposition as a Filterbank[J]. IEEE Signal Processing Letters，2004，11(2)：112-114.

[5] FLANDRIN P，GONCALVES P，RILLING G. Detrending and Denoising with Empirical Mode Decomposition[C]∥European Signal Processing Conference(EUSIPCO2004).[S.l.]:Cite SeerX.2004.

[6] FLANDRIN P. Empirical Mode Decompositions as Data-driven Wavelet-like Expansions[J].International Journal of Wavelets: Multiresolution and Information Processing，2004，2(4)：1-20.

[7] WU Z，HUANG N E. A Study of the Characteristics of White Noise Using the Empirical Mode Decomposition Method[J]. Proceedings of the Royal Society，2004，460(2046)：1597-1611.

[8] WU Z，HUANG N E. On the Filtering Properties of the Empirical Mode Decomposition[J]. Advances in Adaptive Data Analysis，2010，2(4)：397-414.

[9] DONOHO D L，JOHNSON I M. Ideal Spatial Adaptation by Wavelet Shrinkage[J].Biometrika，1994，81(3)：425-455.

[10] DONOHO D L. De-noising by Soft-thresholding[J].IEEE Transactions on Information Theory,1995，l41(3)：613-627.

[11] BOUDRAA A O，CEXUS J C，SAIDI Z. EMD-based Signal Noise Reduction[Z].World Academy of Science，Engineering and Technology，2005，2:33-36.

[12] 张义平.爆破振动信号的HHT分析与应用[M].北京：冶金工业出版社，2008：83-84.

[13] BOUDRAA A O，CEXUS J C.EMD-based Signal Filtering[J].IEEE Transactions on Instrumentation and Measurement，2007，56(6)：2196-2202.

[14] WENG B，BARNER K E. Optimal Signal Reconstruction Using the Empirical Mode Decomposition[J].EURASIP Journal on Advances in Signal Processing，2008：1-12.

[15] LANGBEIN J，JOHNSON H. Correlated Errors in Geodetic Time Series: Implications for Time-dependent Deformation[J].Journal of Geophysical Research，1997，102(B1)：591-603.

[16] ZHANG J，BOCK Y，JOHNSON H，et al. Southern California Permanent GPS Geodetic Array: Error Analysis of Daily Position Estimates and Site Velocities[J].Journal of Geophysical Research，1997，102(B8)：18035-18055.

[17] NIKOLAIDIS R. Observation of Geodetic and Seismic Deformation with the Global Positioning System[D].Califormia：University of Califormia，2002.

[18] 符养.中国大陆现今地壳形变与GPS坐标时间序列分析[D].上海：中国科学院上海天文台，2002.

[19] DONG D，FANG P，BOCK Y，et al. Spatiotemporal Filtering Using Principal Component Analysis and Karhunen-Loeve Expansion Approaches for Regional GPS Network Analysis[J].Journal of Geophysical Research，2006，111(B3)：1-16.

[20] 黄立人，符养.GPS连续观测站的噪声分析[J].地震学报，2007，29(2)：197-202.

[21] 袁林果，丁晓利，陈武，等.香港GPS基准站坐标序列特征分析[J].地球物理学报，2008，51(5)：1372-1384.

[22] 蒋志浩，张鹏，秘金钟，等.顾及有色噪声影响的CGCS2000下我国CORS站速度估计[J].测绘学报，2010，39(4)：355-363.

[23] 张恒璟，程鹏飞，郭英.我国IGS基准站高程时间序列抗差功率谱估计[J].测绘科学，2012，37(3)：52-53.

[24] 张恒璟，程鹏飞.基于经验模式分解的CORS站高程时间序列分析[J].大地测量与地球动力学，2012，32(3)：129-134.