鲁铁定 谢建雄

1 东华理工大学测绘工程学院，南昌市广兰大道418号，330013 2 东华理工大学江西省数字国土重点实验室，南昌市广兰大道418号，330013 3 漳州市测绘设计研究院，福建省漳州市龙溪南路5-58号，363000

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)向待分解信号中加入随机高斯白噪声：

xi(t)=x(t)+nωi(t)

(1)

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

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

(2)

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

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

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

(3)

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

1.2 多尺度排列熵(MPE)

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

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

(4)

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

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

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

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

2 仿真算例分析

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

Y=9sin(0.6πt)+6cos(2πt)+

11sin(3πt)+noise

(5)

图1 仿真信号Fig.1 Simulation signal

图2 各阶IMF与仿真信号的相关系数Fig.2 Correlation coefficients of each IMF and simulated signal

图3 各阶IMF的值

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

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

图4 3种方案降噪后的信号与原参考信号的残差结果Fig.4 Residual results of the signal after thenoise reduction of the three schemes and the original reference signal

表1 不同方法的降噪评价指标

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

3 实例分析

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

图5 BJFS和NVSK站高程时间序列Fig.5 Elevation time series of BJFS and NVSK stations

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

图6 各阶IMF与原始序列的相关系数Fig.6 Correlation coefficients of each IMFand original sequence

图7 各阶IMF的值

图8 基于不同方法获取的高程时间序列对比Fig.8 Comparison of elevation time series obtained by different methods

图9 基于不同方法滤除的噪声序列对比Fig.9 Comparison of noise sequences filtered by different methods

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

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

dnSNR=10 lg(Ps/Pg)

(6)

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

表2 各种方法的降噪效果评价指标

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

4 结 语

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