三种GNSS 高程时序降噪方法的效果对比分析

范小猛，胡川，张重阳，李成洪

(重庆交通大学 土木工程学院,重庆 400074)

0 引言

在全球卫星导航系统(GNSS)坐标时间序列中，不仅包含真实信号，还含有白噪声(WN)、闪烁噪声(FN)、随机漫步噪声(RWN)等多种有色噪声.合理的分离出信号和噪声，有助于获取更准确的测站运动参数估值，对于精化速度场、建立地球框架以及揭示构造变形运动等有着重要的意义[1-2].

经验模态分解(EMD)、整体经验模态分解(EEMD)和小波降噪是三种较常用的噪声分离方法.HUANG等[3]在1998 年讨论了EMD 在时间序列分析中的应用问题.张双成等[4]将基于EMD 的时间序列分析法引入GNSS 坐标序列分析，对中国境内及周边9 个国际GNSS 服务(IGS)站三个坐标分量上的序列数据进行处理，有效地分离出线性趋势项、周期项和噪声项，证明了EMD 法降噪的有效性.随着研究的不断深入，研究人员发现EMD 存在模式混叠的现象[5-6]，为了改善这个问题，WU 等[7]在EMD 的基础上提出了EEMD.张恒璟等[8]采用EMD 和EEMD分别对BJFS 站和URUM 站10 a 左右的高程数据进行处理，分离序列趋势项和周期项运动，发现EEMD可以明显减弱模式混叠现象，但不能完全消除模式混叠，同时EEMD分解时所加入的WN 会导致其分解结果不唯一.小波分析具有多分辨率、低熵性、去相关性等特点，得到了GNSS 时间序列分析研究人员的重视[9-11].魏冠军等[12]采用小波变换阈值降噪对高程分量的坐标时序进行处理，有效减弱了噪声信息，证明了该方法降噪的可行性，但样本序列数据时间长度在3 a 左右，没有涉及到更长时间跨度时序的降噪分析.马俊等[13]提出了小波包系数信息熵法，对KMIN 站近5 a 的三个坐标分量上的时序进行了降噪处理，得到了仅受WN 影响的时间序列，消除了有色噪声的影响，但同样没有对更长时间跨度的坐标时间序列进行分析.

目前，已经有众多学者对EMD、EEMD 和小波降噪有了深入研究，但少有文献考虑到数据长度可能会对降噪效果产生的影响.随着GNSS基准站的不间断观测，相关机构已积累了20 多年的坐标时间序列数据，为研究上述方法在不同时序长度下的降噪效果提供了可能.针对上述问题，本文通过对比分析EMD、EEMD 和小波降噪在不同GNSS 高程坐标时序长度下的降噪效果，采用信噪比(SNR)、相关系数和均方根误差(RMSE)定量评价它们的降噪性能，为相关研究提供一定参考.

1 降噪方法

1.1 经验模态分解

EMD 是一种基于非线性、非平稳数据的自适应时频信号分解方法，在不同尺度上分析数据，提取有用信息[3].其基本思想是通过筛选，把原始信号X(t)分解成m个不同时间尺度的本征模态分量(IMF)和残差项rn，分解流程如下：

1)寻找原始信号中所有的局部极大、极小值点，利用三次样条函数分别拟合上述点，得到关于时间的上、下包络线，计算两条包络线的平均包络线m1(t)，用原始信号X(t)减去平均包络得到新的待分解信号h1(t)

2)判断h1(t)是否满足以下条件：一是待分解的信号中的极值点的个数与过零点的个数应相等或最多相差1；二是在任意点上，上下包络线的均值为0，即上下包络线应对称.如果满足上述条件，将h1(t)作为分解得到的第一个分量 I MF1(t)，然后进行第3)步操作，否则将h1(t)作为新的原始信号重复第1)步操作.

3)计算剩余信号r1=X(t)-IMF1(t)，然后将r1作为新的待分解信号.

4)判断r1是否满足以下条件：r1为单调变化的残差序列，且不能继续分解得到IMF 分量.如果满足，EMD 分解过程结束，残差项rm=r1，否则重复流程1)～3)操作.原始信号经过m次分解可以表示为

为了获取不含噪声的GNSS 坐标时间序列，需要对上式中分解出的IMF 分量进行重构.首先，计算各分量与原始序列间的相关系数；然后识别第一个取局部极小值的相关系数对应的分量IMFk，将其作为信号与噪声之间的分界层；最后，将分界层与之前的IMF 分量层相加作为高频噪声，然后用原始信号减去重构高频噪声得到降噪后信号[3].

1.2 整体经验模态分解

EMD 的缺点之一是模式混叠：信号被EMD 分解为m个不同的IMF 分量后，不同频率的信号会存在于同一个IMF 分量中，即一个IMF 分量包含着不同尺度的信号.EEMD 通过在分解前向原始信号中添加白噪声w(t)，生成待分解信号，改善模式混叠问题.

式中，下标i表示第i次分解.其余分解步骤与EMD 相同，此处不再赘述.将多次分解得到的尺度相同的IMF 分量整体取平均，作为当前尺度的IMF 分量，从而得到最终的一系列不同尺度的IMF 分量数据，按照相关系数准则重构信号即可得到降噪之后的GNSS坐标时间序列.本文根据文献[14]实验结果，将加入的WN 与原始信号的标准差之比设为0.1，整体平均次数设置为100 次.

1.3 小波降噪

多分辨率分析是小波分析的重要概念，它是利用特定的小波基函数的伸缩变换，将GNSS 坐标时间序列在不同分辨率上进行分解，得到能够反映时间序列概貌的低频部分以及能够反映时间序列细节的高频部分，进而对不同分辨率上的信号进行分析[12].小波阈值降噪是在多分辨率分析的基础上，通过设置一定的阈值对分解之后的各个尺度上的信号进行处理，达到分离噪声的目的.小波阈值降噪的具体步骤如下：

1)小波分解.选取合适的小波基函数，对GNSS坐标时间序列进行小波分解，得到不同尺度上的小波系数.

2)阈值降噪.对分解得到的各个尺度上的信号进行非线性阈值降噪处理，剔除小波系数小于阈值的小波项.

3)重构信号.将剔除小波项之后的各尺度上的信号重构得到降噪后的GNSS 坐标时间序列.

经过实验对比，本文采用Coif3 小波对原始数据进行分解，根据软阈值规则对数据进行降噪处理.

2 实验数据和方案

2.1 实验数据

由于我国境内的IGS 站时间序列在水平方向上主要以线性为主，而在高程方向上表现出复杂的非线性运动[15-16]，因此本文以高程方向上的坐标时序为例进行分析.为保证实验所需时序长度，选取我国境内BJFS、KMIN、LHAS、SHAO、URUM、WUHN 共6 个IGS 站点所对应的高程方向单日解坐标时序，各站点完整数据的起止历元均为1999.17—2020.39，截取其中2000—2019 年的观测数据作为本次实验的原始数据.

采用线性拟合去除原始数据中的趋势项，并且根据3σ 准则剔除孤立值.另外，为了避免由于数据缺失造成模式混叠现象，以及满足小波分析要求数据具有均匀取样和零均值的特性，利用线性插值对数据缺失的部分进行插补，得到适用于降噪分析的连续数据.

2.2 实验方案

实验方案共分为三种：

方案一：将每个站点共20 a 的坐标时间序列数据依次分为4 段，分别为2000—2004、2005—2009、2010—2014、2015—2019 时间段，得到24 个时序长度为5 a 的样本序列；

方案二：将每个站点共20 a 的坐标时间序列数据平均分为2 段，分别为2000—2009、2010—2019，得到12 个时序长度为10 a 的样本序列；

方案三：将每个站点2000—2019 共20 a 的坐标时间序列数据分别整体作为一段，得到6 个时序长度为20 a 的样本序列.

降噪效果通过降噪后信号与原始信号之间的SNR、相关系数以及RMSE 来判断.其计算公式分别为：

式中：X(t)为原始序列；为降噪后的真实信号；N为采样点个数.SNR 体现的是噪声信号在信号中所占的比重，其值越大，表示噪声含量越低；相关系数越接近1，表示两者的相似度越高，降噪效果越好；RMSE 表示的是两者之间的偏差程度，其值越小代表降噪效果越好.

3 实验结果与分析

3.1 降噪结果

经过计算，三种方法在不同样本序列长度下所分解出的分量个数相同，从5 a 到20 a 分别为10、11、12 个.对于EMD 和EEMD 而言，根据各IMF 分量与原样本序列的相关系数得到的信号与噪声间分界层不同：EMD 法得到的不同样本序列分界层数分别为第5、4、5 层，而EEMD 法得到的不同样本序列分界层数分别为第2、5、5 层，根据上述分界层对样本序列进行降噪处理.

图1 和图2 描述了不同IGS 站点高程方向上经过降噪的坐标时间序列.从图中可以发现三种方法都可以有效的对坐标时间序列数据进行降噪，降噪后信号表现出明显的周期性变化，曲线变化也较为平滑.其中，图1 中BJFS 站的5 a-2 样本序列以及KMIN 站的5 a-4 样本序列的光滑度明显优于其余相同时间跨度的样本序列，其原因可能是在这两个样本数据中缺失历元数较少，数据较为完整.例如BJFS 站的5 a-2样本序列的缺失历元数为54，仅占样本总数的3%，远小于其它曲线光滑度较低的样本序列缺失历元数，说明时间序列中缺失数据的存在会对降噪效果产生一定影响.

图1 BJFS 站和KMIN 站降噪结果对比

图2 SHAO 站和URUM 站降噪结果对比

图2 中SHAO 和URUM 两个站点，EMD 降噪后的曲线在个别时间段呈现出明显的异常波动.其中SHAO 站的5 a-2 样本序列数据在2007 年前后波动最大，最大达到了80 mm，5 a-4 样本序列数据中异常波动较多，但波动幅度较小，其余站点也存在相似的异常波动情况.这种异常波动现象可能是原始时间序列数据中存在长时间的连续中断，如SHAO 站连续中断的历元数为411，缺值率达到了22%，即便是后续对缺失处进行了插补，依然会对EMD 降噪结果造成偏差.

此外，在上述出现异常波动的时间点，EEMD 和小波降噪得到的序列趋势相似，振幅为3～5 mm，趋势较平缓，说明即使数据存在长时间的连续缺失，这两种方法也可以有效减弱插值的影响，削弱信号中存在的噪声.

3.2 降噪结果对比分析

为了进一步说明EMD、EEMD 和小波降噪的效果，分别计算不同样本长度下各方法降噪结果的SNR、相关系数和RMSE，其中样本长度为5 a 和10 a 时，将各指标的平均值作为最终结果，如表1 所示.

由表1 可知：

表1 各站点评价参数统计 m

1) EMD 受初始数据质量的影响较大.在数据较完整的站点，EMD 降噪序列各评价指标略差或持平于EEMD 和小波降噪结果；而在缺失比例大的站点，各评价指标值明显较差.

2)相比于EMD，EEMD 不仅改善了模式混叠问题，还在一定程度上提高了EMD 分解的精度，即便是在缺值数较多的SHAO 站，EEMD 降噪序列的相关系数值大部分在0.7 以上，降噪效果稳定.

3)在5 a 和10 a 的样本时间序列数据中，小波降噪序列的评价参数值明显优于EMD 和EEMD；当时间跨度为20 a 时，EEMD 降噪序列的评价参数与小波相近.这说明在这三种方法中，小波降噪的适用性更好，而EEMD 更适用于对时间跨度大的坐标时间序列进行降噪分析.

4)大部分站点的EEMD 降噪序列和小波降噪序列随着时间跨度的增加，三个评价参数值均有所增加，难以根据某一个参数准确说明其降噪性能的强弱.

为了进一步探究EMD、EEMD 和小波降噪对有色噪声的剔除能力，分别计算三种方法降噪序列的谱指数，谱指数越接近零，在一定程度上说明有色噪声的振幅越小，对真实信号的影响也越小[13]，结果如图3 所示.

图3 三种方法得到的降噪序列的谱指数

由图3 可以得出以下结论：1)EEMD 降噪后时序的谱指数大于EMD 降噪序列的谱指数，剔除有色噪声的能力更好，但是否因EEMD 方法所加入WN 产生的影响还有待进一步研究；2)小波分析降噪在不同时间跨度的时序数据下都可以更好的抑制有色噪声，使得降噪后的信号更加可靠.

4 结论

本文通过设置不同样本长度的降噪实验，对EMD、EEMD 以及小波分析三种方法的降噪效果进行对比，并利用SNR、RMSE 和相关系数定量分析了其降噪性能.实验表明：

1)EMD 易受坐标时间序列自身质量的影响，降噪效果相对较差；EEMD 和小波降噪在时间序列存在长时间中断，质量较差时，依然可以得到较好的降噪结果.

2)随着时序样本长度的增加，三种方法的降噪性能均有提升.其中，小波降噪的适用性最广泛，在不同的样本数据中，降噪性能均优于EMD 和EEMD；EEMD 在时间跨度较长的情况下降噪性能更好，与小波降噪性能接近；EMD 降噪效果相对最差.

3)小波降噪相比于EMD 和EEMD，所得到的降噪序列中含有的有色噪声更少，更有利于对数据进行后续的分析和使用.