基于CEEMDAN和小波包多阈值的GNSS高程时间序列去噪方法

于宏旭 文汉江 刘焕玲 董 杰 蔺文奇

1 辽宁工程技术大学测绘与地理科学学院，辽宁省阜新市玉龙路88号，123000 2 中国测绘科学研究院，北京市莲花池西路28号，100036

常用的GNSS高程方向去噪方法有Kalman滤波[1]、小波去噪[2-3]、经验模态分解(EMD)及其改进算法[4-7]等。但这几种方法都存在不足：Kalman滤波在处理非线性时间序列时效果不佳；小波去噪在设置阈值时只考虑了高频噪声的影响，对于低频噪声的去噪效果并不理想[8]；CEEMDAN算法舍弃高频分量，保留低频分量，容易丢失有用信息且忽略了低频分量中包含的噪声。小波包多阈值分解基于频率高低分段选取合适的阈值，能够高效剔除各频段的噪声，同时保留高频段中的有用信息[9]。

基于以上研究，本文结合CEEMDAN和小波包多阈值去噪的优势，引入一种CEEMDAN+小波包多阈值方法对GNSS高程时间序列进行去噪。

1 算法原理

1.1 CEEMDAN算法原理

CEEMDAN算法原理如下[10]：

(1)

(2)

残差为：

(3)

(4)

残差为：

(5)

3)按照上述2个步骤进行分解，直到待分解信号只有2个极值点，即为单调函数时停止计算，此时CEEMDAN中原始待分解信号可表示为：

(6)

1.2 小波包多阈值原理

小波包多阈值变换包括小波包分解、阈值处理、小波重构等3个步骤[11]。最优小波基、最佳分解层数和小波包去噪阈值的选取是小波包分解的关键因素。本文选择较常用的dbN和symN小波基函数。分解层数的选择对小波包去噪效果尤为重要，任何信号都存在一个去噪效果最好的分解层数，通常分解层数越大，噪声和信号表现的不同特性越明显，越有利于二者的分离；但分解层数越大，重构信号失真越严重，会在一定程度上影响去噪效果。常用信噪比(SNR)和均方根误差(RMSE)综合确定分解层数：

(7)

(8)

式中，x(i)为原始信号，x′(i)为降噪后的信号，m为信号长度。RMSE反映去噪后信号与原始信号之间的区别，RMSE越小则信号去噪效果越好；SNR反映信号与噪声的比例，重构信号的信噪比越高则信号的重构效果越好。

小波包多阈值去噪方法可根据频率的高低选择阈值，避免了单阈值去噪不能充分考虑信号与噪声分布、极易去掉中高频信号中有用信息、产生过度去噪的缺点。目前小波包分解常用的阈值准则有固定阈值(sqtwolog)准则、自适应阈值(rigrsure)准则与极大极小阈值(minimaxi)准则。sqtwolog准则是将小波包系数全部置0，能够较强地去除噪声，适合处理高频信号；rigrsure准则是基于无偏似然估计原理的自适应阈值准则，仅将部分系数置0，适合处理中低频信号；minimaxi准则是固定形式的阈值准则，是一种较为保守的处理方法，适合处理中低频信号。因此，本文在GNSS高程时间序列去噪时，采用rigrsure准则处理低频段，采用minimaxi准则处理中频段，采用sqtwolog准则处理高频段[12]。对于阈值函数类型，选取最常用的软阈值函数。本文去噪流程见图1。

图1 本文去噪流程Fig.1 Denoising flow chart in this paper

2 仿真实验

采用MATLAB开展仿真实验。因高程时间序列季节性变化近似于正弦曲线，且包含周期项、半周期项和趋势项等信息，因此将y=6sin(0.01t)+2cos(0.03t)-0.005t+f(t)+w(t)作为原始信号，其中，f为有色噪声，w为白噪声，信号长度为3 650。根据高程时间序列的噪声特点，将有色噪声振幅设置为白噪声振幅的2倍。利用CEEMDAN对信号进行分解，在添加高斯白噪声时，为削弱添加的白噪声对降噪结果的影响，设置整体平均次数为100，添加的白噪声信噪比为0.2[13]。信号、添加的白噪声和有色噪声分别如图2(a)～2(c)所示，加入噪声后的混合信号如图2(d)所示。

图2 仿真实验的模拟信号Fig.2 The analog signal used in the simulationexperiment

利用CEEMDAN分解得到10个频率逐渐降低的IMF分量和1个残余项(图3)。从图3看出，随着信号分解层数的增加，噪声逐渐减少。

单独使用CEEMDAN去噪时，利用分界点将IMF分为噪声IMF和信号IMF。然而，在噪声IMF中仍有部分有用信息，在信号IMF中也包含部分噪声。为克服该缺点，本文不进行信号与噪声分界点的判断，而是采用小波包多阈值法对每个IMF分量进行小波包多阈值去噪。首先计算每个IMF分量中的根节点能量占IMF总能量的百分比，根据能量不同将节点分为低、中、高频3个部分，然后分别采用rigrsure、minimaxi和sqtwolog准则进行数据处理。

在利用小波包多阈值去噪时，选择对称性与正则性较好的sym小波(sym4、sym5、sym6)和db小波(db3、db4、db5)作为小波基函数，分解层数定为3、4、5层。利用RMSE、SNR评价最终降噪结果。经过多次对比实验(图4)，发现小波基选择sym5、分解层数为4层时得到的RMSE最小，SNR最大，说明此时小波包多阈值去噪效果最优。

图4 利用RMSE和SNR判断最佳小波基与分解层数Fig.4 Using RMSE and SNR to determine the best wavelet basis and the number of decomposition layers

利用小波包多阈值去噪对各IMF分量去噪，将去噪后的IMF分量与残余项进行重构，得到去噪后的信号，结果见图5。

图5 CEEMDAN+小波包多阈值去噪结果和噪声Fig.5 CEEMDAN and wavelet packet multi-thresholddenoising results and noise

为验证本文方法的去噪效果，分别采用本文方法(CEEMDAN+小波包多阈值)、EMD、CEEMDAN、小波去噪和小波包多阈值去噪分析仿真信号，结果见图6。选择SNR、RMSE评价5种方法的降噪效果，结果见表1。

图6 5种去噪方法对比Fig.6 Comparison of five denoising methods

表1 去噪效果比较

从图6可以看出，5种方法均能对模拟信号去噪，但EMD去噪结果过于平缓，丢失了部分信息；与小波包多阈值去噪、小波去噪、CEEMDAN去噪相比，CEEMDAN+小波包多阈值方法去噪后RMSE最小、SNR最大，说明该方法去噪效果最优。

3 实例分析

选取JPL发布的拉萨站(LHAZ)2000～2020年单天解高程时间序列进行降噪分析，原始时间序列中已剔除序列初值、地震和非地震跳变的影响，高程时序缺值或者阶跃不影响高程序列降噪效果。利用3倍中误差准则去除时间序列中剩余的粗差，提高数据质量。LHAZ站原始高程时间序列见图7。

图7 LHAZ站预处理后的高程时间序列Fig.7 Height time series after preprocessing at LHAZ station

利用CEEMDAN算法将时间序列分解成12个IMF分量和1个残余项，将各IMF分量进行小波包多阈值去噪，得到降噪后的各IMF分量，并与残余项重构得到去噪后的信号，结果见图8。

图8 LHAZ站CEEMDAN+小波包多阈值去噪信号和所剔除的噪声Fig.8 CEEMDAN and wavelet packet multi-threshold denoising signal and eliminated noise at LHAZ station

使用上节中5种方法对LHAZ站高程时间序列进行去噪处理，结果见图9，去噪效果见表2。结合图9和表2可以看出，5种方法均能很好地去除信号中的噪声，其中，CEEMDAN+小波包多阈值去噪法的SNR最大、RMSE最小，表明其去噪效果最好。小波去噪过程中固定了小波阈值，只对低频信号进行分解，没有对高频信号进行分解，导致无法保留高频信号中的有用信息。EMD去噪时在峰值附近会产生过度去噪的现象。

图9 5种去噪方法对比Fig.9 Comparison of five denoising methods

表2 5种去噪方法效果比较

4 结 语

GNSS高程时间序列中含有大量噪声，本文针对EMD、CEEMDAN、小波、小波包多阈值等传统去噪方法去噪不彻底或去噪过度的缺点，尝试利用CEEMDAN +小波包多阈值的方法对信号进行去噪，并分别利用仿真信号、LHAZ站高程时间序列进行实验。结果表明，5种去噪方法均可以达到去噪效果，但本文方法去噪后各项指标均最优，效果最好。