GPS高程时间序列中包含各类信号和噪声，因受测站外部环境、观测技术误差等因素的影响，表现出明显的周期性变化。利用传统周期函数模型估计的周期项振幅为常数，但实际上，时间序列周期性信号的振幅是随时间变化的^[1-2]，因此采用传统方法获得的分析结果不能很好地反映高程时间序列的真实运动。

将随机噪声进行有效剔除是获取合理可靠的周期性变化信号的关键。刘焕玲等^[3]利用经验模态分解(EMD)方法分析IGS站高程时间序列，但该方法通过Hilbert频谱图识别信息IMF时存在较大的主观性，缺乏统一的标准；张双成等^[4]结合相关系数原理，应用EMD方法获得的重构信号与原始时间序列较为一致，但不同测站各分量的噪声特性存在差异^[5]，并且当EMD处理的时间序列存在中断时会产生较为明显的模式混叠现象^[6]，导致部分站点依据相关系数首个局部极小值识别噪声IMF分界时出现区别效果不明显的情况，从而影响降噪结果的准确性。

针对上述问题，综合考虑整体经验模态分解(EEMD)能有效减弱模式混叠及多尺度排列熵(MPE)可以反映时间序列复杂程度的优势，提出一种EEMD-MPE阈值降噪方法，并通过仿真算例和部分IGS站高程时间序列对该方法进行验证。

1 算法原理 1.1 EEMD算法原理

EEMD^[7]的本质是基于高斯白噪声在所有频率上具有相等能量的特性，通过在待分解信号中加入高斯白噪声后再进行多次EMD处理，从而实现减弱或消除模式混叠的方法。EEMD算法的基本步骤如下：

1) 向待分解信号中加入随机高斯白噪声：

$ {x_i}\left( t \right) = x\left( t \right) + n{\omega _i}\left( t \right) $

(1)

式中，x(t)为待分解信号；n为加入白噪声的幅值系数；ω_i(t)为加入的白噪声，i=1, 2, …, N。

2) 采用EMD处理每个x_i(t)，获得k个IMF和1个余项r_i(t)：

$ {x_i}\left( t \right) = \mathop \sum \limits_{j = 1}^k {c_{ij}}\left( t \right) + {r_i}\left( t \right) $

(2)

式中，c_ij(t)为第i次分解得到的第j个IMF。

3) 循环步骤1)和2)，但每次均加入不同的随机高斯白噪声。

4) 计算分解所得的IMF的总体均值，得到EEMD的最终结果：

$ x\left( t \right) = \frac{1}{N}\left( {\mathop \sum \limits_{i = }^N \mathop \sum \limits_{j = 1}^k {c_{ij}}\left( t \right) + \mathop \sum \limits_{i = 1}^N {r_i}\left( t \right)} \right) $

(3)

关于EEMD的2个重要参数，采用Wu等^[7]的建议，设置整体平均次数N=100，白噪声的幅值系数为0.2。

1.2 多尺度排列熵(MPE)

多尺度排列熵是一种能够反映信号不确定性和不规则性的非线性分析方法，与排列熵^[8]相比具有更好的稳定性和更强的抗噪声能力，其基本思想是对时间序列多尺度化后，计算对应尺度上的排列熵。多尺度排列熵的计算步骤参考文献[9]。

为了使多尺度排列熵值位于[0, 1]，通常对其进行归一化处理：

$ \overline {{\rm{MPE}}\left( X \right)} = \frac{{{\rm{MPE}}\left( X \right)}}{{{\rm{ln}}\left( {m!} \right)}} $

(4)

对于关键参数的设置，根据郑近德等^[10]的建议，取m=6，t=1 s，计算各IMF分量在尺度因子s位于[1,15]的排列熵均值作为最终MPE值。MPE值越接近1，反映序列的随机波动性越大，不规则程度越高。排列熵阈值一般取0.55~0.6，即MPE值大于设定阈值的分量随机程度较大，视为噪声分量，直接舍去。

1.3 基于EEMD-多尺度排列熵(MPE)的降噪新方法

通过MPE值判断出噪声分量的层数，阈值的选取尤为重要，但目前关于排列熵阈值的设定缺乏明确的标准。选取某一固定阈值判断出噪声分量层数，该方法难免会导致有用信息被当成噪声剔除，造成信号失真。针对这一问题，本文结合IMF分量MPE值的分布特点，提出一种新的降噪方法，主要步骤如下：

1) 利用EEMD处理GPS高程时间序列，获得一系列IMF和1个余项；

2) 分别计算出各阶IMF和余项的MPE值；

3) 根据MPE值将各IMF细分为噪声IMF、混合IMF和信息IMF三类，即MPE值大于0.65的分量归类为噪声IMF，MPE值位于[0.55, 0.65]之间的分量归类为混合IMF，MPE值小于0.55的分量归类为信息IMF；

4) GPS高程时间序列的有用信息主要集中在混合IMF的"干净"部分和信息IMF中，因此对高频噪声IMF直接剔除，然后采用改进阈值函数^[11]处理混合IMF；

5) 重构利用阈值降噪后所得的结果与信息IMF即可获得最终的输出序列。

2 仿真算例分析

由于GPS高程时间序列包含的噪声种类较多，频率分布范围较广，为保证与实际情况相符，需构造复合仿真信号来模拟GPS高程时间序列的实测信号。仿真信号的采样频率为365 Hz，采样点个数为4 000，信号包括3个周期项和随机噪声项，其中随机噪声项(noise)包含高斯白噪声和有色噪声。仿真信号如图 1所示，其表达式为：

$ Y = 9{\rm{sin}}(0.6{\rm{ \mathsf{ π} }}t) + 6{\rm{cos}}(2{\rm{ \mathsf{ π} }}t) + {\rm{ }}11{\rm{sin}}(3{\rm{ \mathsf{ π} }}t) + {\rm{noise}} $

(5)

从图 1可以看出，仿真信号设置合理，与实际情况较为接近，采用EEMD算法自适应分解仿真信号，获得若干个IMF及余项。为了更好地展现本文方法的优势，设置如下实验方案进行降噪效果对比：方案1，采用相关系数法；方案2，采用MPE法(MPE阈值取0.6^[10])；方案3，采用本文MPE-阈值新方法。计算各阶IMF与仿真信号的相关系数和MPE值，结果分别如图 2和3所示。

由图 2可知，相关系数的局部极小值首次出现在第4阶IMF，方案1识别前4阶IMF为噪声IMF，因此重构IMF₅及剩下的IMF，获得降噪后的信号。图 3为EEMD处理的各阶IMF的MPE值分布情况，由图可见，IMF₁~IMF₅的MPE值均位于[0.6, 0.9]，超过设定的阈值，因此方案2识别前5个IMF为噪声IMF，重构IMF₆之后的分量获取降噪后的信号。根据方案3的方法，由图 3可知，IMF₁~IMF₄的MPE值均大于0.65，属于噪声IMF，表现出较强的随机性，因此直接从仿真信号中滤除，而IMF₅的MPE值位于[0.55, 0.65]之间，表明IMF₅属于混合IMF，采用阈值降噪方法处理IMF₅，重构阈值降噪后所得结果与余下的IMF，即可获得最终的降噪结果。

为了分析不同方案的处理效果，选用均方根误差(RMSE)和信噪比(SNR)^[12]指标进行降噪质量评价。一般认为，RMSE值越小，降噪后的信号与原始参考信号越接近，降噪质量越好；而SNR则相反，其值越高，降噪效果越显著。

3种方案降噪后的残差结果如图 4所示，表 1统计了3种方案降噪后的RMSE和SNR值。

由图 4和表 1可知，方案3获得的降噪后信号与原始参考信号的残差结果最小，即RMSE值最小，SNR值最大，表明该方案降噪性能最优；而方案1的降噪效果最差，这是由于方案1未能准确识别出噪声层数；尽管方案2的降噪性能略优于方案1，但该方案将混合IMF当成噪声IMF进行了剔除，在降噪的同时将分量信号中的有效成分也一并舍去，造成信号失真。

3 实例分析

为进一步检验本文方法的有效性，以BJFS和NVSK站高程时间序列(图 5)为例进行降噪分析，数据来源于SOPAC (Scrips Orbit and Permanent Array Center)提供的去除了均值、趋势项、粗差、同震跳变和非地震跳变影响后的GPS单天高程时间序列。

从图 5可以看出，两个IGS站高程时间序列中仍存在部分粗差，因此采用"3σ"准则进行粗差识别和剔除，以提高数据的质量。为能有效并准确地提取序列中的周期项，采用3种基于EEMD的降噪方法对高程时间序列进行处理和分析，以验证本文方法的优越性。

利用EEMD算法分别对两个IGS站高程时间序列进行分解，然后分别计算出各阶IMF与原始时序的相关系数和MPE值，结果见图 6和7。

采用相关系数法判定噪声分量分界时，由图 6可知，BJFS站的相关系数局部极小值首次出现在IMF₃，所以从IMF₄开始进行重构获得降噪后的序列；而NVSK站的相关系数局部极小值首次出现在IMF₈，识别出前8阶IMF为高频噪声分量，显然该方法失效。由图 7可见，两个IGS站IMF₁~IMF₄的MPE值均位于[0.7, 0.9]之间，利用MPE法判定噪声分量时发现，前4阶IMF超过了设定的阈值(0.6)，重构IMF₅之后的分量便可获取降噪后的序列；采用本文MPE-阈值法进行降噪处理时发现，IMF₁~IMF₄的MPE值均大于0.65，属于噪声IMF，因此首先将IMF₁~IMF₄从原始序列中直接剔除，而IMF₅的MPE值位于[0.55, 0.65]之间，表明IMF₅属于混合IMF，所以将IMF₅作阈值降噪处理，重构阈值降噪后的序列和余下的分量即可得到最终降噪后的高程时间序列。图 8为基于不同方法的降噪结果对比，图 9为基于不同方法滤除的噪声序列对比。

图 8和9表明，对于BJFS站，MPE法降噪效果优于相关系数法，而本文MPE-阈值法降噪后的序列与原始时间序列的变化趋势更吻合，在最大限度滤除噪声的基础上较好地保留了时间序列的有用成分，与另外两种方法相比，降噪效果最佳。对于NVSK站，由于相关系数法无法准确判断出噪声分界层数，该方法失效；而MPE法尽管能够准确判断出噪声分界层数，但不能识别出混合IMF，导致降噪效果欠佳；本文MPE-阈值法能够将混合IMF中残留的噪声去除，同时保留信号的有用成分，降噪的结果优于MPE法。

由于实测信号的有效信号和噪声功率均未知，仍采用信噪比评价降噪效果是不准确的，因此本文引入降噪误差比(d_nSNR)^[13]指标评价降噪质量，d_nSNR值越小降噪效果越显著。

$ {d_{{\rm{nSNR}}}} = 10\;\lg ({P_s}/{P_g}) $

(6)

式中，P_s为含噪信号的功率，P_g为滤除的噪声功率。同时，选用平滑度(r)^[12]作为降噪质量评价指标，r值越小信号序列越光滑，降噪效果越好。表 2统计了各种方法的降噪效果评价指标。

由表 2可知，仅从平滑度(r)来看，采用本文方法处理BJFS站数据结果的r值最小，仅为0.001 6，表明本文方法的降噪性能最好，同时对比3种方法的d_nSNR值也验证了这一观点；但对于NVSK站，MPE法和本文方法的r值均为0.003 9，无法体现降噪质量的优劣，而本文方法的d_nSNR值小于MPE法，说明本文方法的降噪性能优于MPE法，同时在评价降噪性能方面也反映了d_nSNR指标比平滑度(r)指标更具有适用性。

4 结语

本文综合EEMD算法和MPE的优势，提出一种基于改进阈值函数的EEMD-MPE联合降噪方法。该方法对EEMD处理后的混合IMF进行降噪，并重构降噪后的序列与余下分量，获得最终降噪后的序列。通过对仿真信号进行详细的降噪分析，验证本文方法的降噪性能优于相关系数法和MPE法。将该方法应用于IGS站GPS高程时间序列的降噪处理中，结果表明，本文方法在降噪性能上具有明显的优势，获得的降噪后时间序列更符合基准站的实际运动情况，可为进一步研究地壳运动及分析形变模式提供更加稳定、可靠的数据基础。