联合EMD-HD和小波分解的GNSS坐标时间序列降噪分析

杨 兵，杨志强，田 镇，陈 祥

长安大学地质工程与测绘学院，陕西 西安 710054

长期积累的GNSS基准站坐标时间序列广泛应用于大地测量学、地球物理科学和环境学等领域，为高精度全球参考框架、地震等自然灾害的监测与预警，以及地球深部动力学研究提供了重要的基础数据[1-4]。然而，由于受到数据处理策略、模型误差和接收机噪声等GNSS技术误差，以及包括环境负荷、测站热膨胀效应在内的地球物理因素的影响，GNSS基准站坐标时间序列中除反映构造运动的趋势项外，还包括了各种非线性信号(主要包括周年项、半周年项)和噪声，严重制约了其应用[5-8]。因此，GNSS测站坐标时间序列的降噪分析是目前的研究热点。

经验模态分解(empirical mode decomposition,EMD)作为一种非平稳、非线性信号的处理方法[9]，能够较好地识别并剔除地震信号[10]、GNSS坐标时间序列[11-13]和高精度测量仪器的数据[14-16]中的噪声。EMD方法将GNSS坐标时间序列自适应地分解成一系列从高频到低频的本征模态函数(intrinsic mode function，IMF)，再基于一定的筛选准则识别噪声和信号分量，从而达到降噪效果[9,16]。为了避免过度降噪和降噪不足，选取合适的筛选准则至关重要。此外，EMD方法分解得到的IMF分量中可能存在尺度互异或相似的信号，即所谓的模态混叠效应，使得信号分量中混有部分噪声，严重影响降噪效果[9,17-18]。

国内外学者对EMD方法的筛选准则进行了大量研究。文献[10]利用基于连续均方误差(consecutive mean square error，CMSE)的EMD滤波方法[19]有效压制了微震信号中的随机噪声。文献[11,13]选择能量密度E和平均周期T的乘积ET[20]作为分界指标，验证了该指标识别GNSS坐标时间序列噪声的可行性。文献[12]在各IMF分量与原始信号的相关系数[11,21]发生明显变化处开始重构信号，显著提高了GNSS坐标时间序列的信噪比。但是，文献[12]的方法需要将各IMF分量与原始信号的相关系数和预先设置的阈值进行比较，某些情况下信号的相关性过强或过弱，会造成噪声分量的错误识别。文献[10—13]均是根据筛选准则第一次显著变化判断分界分量，具有一定的主观性，使降噪过程不稳定。文献[11—13]选用的测站仅包括我国区域内的IGS站，不能较全面地反映EMD方法对我国GNSS测站的降噪效果。

针对EMD的模态混叠效应，文献[17—18]研究了EMD分解过程的特点，表明连续数据的EMD分解过程中也存在模态混叠效应，并分别提出了集合经验模态分解和互补集合经验模态分解算法以改善该现象。虽然这两种改进算法在一定程度上抑制了模态混叠效应，但需要预先加入噪声，造成计算量增加，影响计算效率的同时也可能导致降噪不足[22]。

事实上，数据分布形状的概率密度函数在一定程度上可以反映两种信号的差异，那么度量概率相似性的Hausdorff距离(Hausdorff distance,HD)便可以界定噪声分量与信号分量[16]。文献[16]通过空气动力学数据证明了HD方法在降噪方面的有效性和优势。因此，本文将该筛选准则引入EMD方法对GNSS坐标时间序列的降噪处理中。而另一方面，小波分解(wavelet decomposition，WD)作为一种时频分析方法，能够基于给定的小波基函数和分解层数对时域信号进行多尺度分解，将得到的高频分量和低频分量的小波系数进行阈值滤波、重构，从而实现含噪信号的降噪处理[23-25]。

基于以上分析，本文利用小波分解对产生模态混叠效应的IMF分量进行处理，探索其对EMD方法降噪和模态混叠效应的改善效果，提出将Hausdorff距离与EMD相结合并联合小波分解的降噪算法EMD-HD&WD，对我国大陆构造环境监测网络(crustal movement observation network of China,CMONOC)的GNSS垂向坐标时间序列进行降噪处理。本文在引入新筛选准则的同时改善模态混叠效应,以拓展EMD方法在GNSS数据处理中的应用，并提出将GNSS测站的运动模型的精度作为降噪效果对比的指标，以研究降噪方法对GNSS运动模型的优化效果。

1 原理与方法

1.1 Hausdorff距离(HD)原理

HD通过计算两个点集之间的概率密度函数来度量相似性。其他度量距离方法多以点集中两个点的距离表述空间目标的距离，而HD则兼顾目标的整体波形和点集中的异常值[16]。其原理如下：

有限集合A={a1,a2,…,am}和B={b1,b2,…,bn}分别是GNSS测站的坐标时间序列和EMD算法分解得到的IMF分量所对应的概率密度函数值(probability density function，PDF)。集合A中元素a到有限集合B之间的距离定义为

(1)

那么，定义A和B之间的单向HD为

(2)

式中，max为最大值操作。

同理，B和A之间的单向HD可以定义为

(3)

单向HD具有非对称性，这是极大极小函数的特性。绝对HD为2个单向HD的最大值，即

HD(A,B)=max{h(A,B),h(B,A)}

(4)

为了利用HD准则筛选IMF分量，首先利用EMD算法将GNSS测站的坐标时间序列x(t)分解为n个IMF分量，然后估计x(t)及各IMF分量的PDF，分别记为PDF(x)与PDF(IMF)={PDF(IMF(1))，PDF(IMF(2))，…，PDF(IMF(i))，…,PDF(IMF(n))}，最后根据式(2)—式(4)计算x(t)与各IMF分量之间PDF的HD值，记为HD(i)

HD(i)=HD{PDF(x),PDF(IMF(i))}

(5)

式中，i=1，2,…,n。

HD的第1个全局极大值点所对应的阶数记作m，则第一阶至第m阶的IMF分量为噪声分量，其余为信号分量，即

(6)

式中,arg max{·}为全局极大值操作。

1.2 EMD-HD&WD算法

小波分解具有良好的时频特性和多尺度分解能力，可以对产生模态混叠效应的IMF分量进行多分辨率分析，剔除其中的噪声，而HD根据GNSS坐标时间序列与EMD分解得到的各IMF分量的概率密度函数进行相似性度量，判断噪声与信号分量。本文结合小波分解和HD的特点和优势，提出了EMD-HD&WD算法。该算法一方面引入了HD筛选准则，可以避免出现其他指标的主观性或陷入局部最小值的情况，减少某些情况下的误判，另一方面联合了小波分解，能够在不需要先验信息的情况下改善模态混叠效应，提高EMD的降噪能力。

EMD-HD&WD算法的具体步骤如下：

(1) 依据EMD算法对GNSS坐标时间序列进行分解，得到n阶IMF分量

(7)

式中，IMFi(t)代表从原始信号中分离出来的频率从高到低的IMF分量；R(t)代表趋势项分量。

(2) 基于核密度函数估计原始信号及各IMF分量的PDF。

(3) 利用HD比较IMF分量与原始信号的PDF的相似性，确定m。

(4) 将m阶之前的各IMF分量重构,得到初始噪声N0(t)

(8)

(5) 选取coif5小波基和8阶分解层数[25]，对m之后的IMF分量(不包括趋势项分量)进行小波分解降噪，得到第j阶IMF分量的降噪后信号S1j(t)和噪声N1j(t)

IMFj(t)=S1j(t)+N1j(t)
(j=m+1,m+2,…,n)

(9)

(6) 将初始噪声与小波分解得到的噪声重构得到最终噪声N1(t)，最终信号S(t)则为S1(t)与趋势分量之和

(10)

1.3 评价指标

1.3.1 复合指标T值

传统降噪方法的质量评价指标主要包括均方根误差(RMSE)、信噪比(SNR)、相关系数(ρ)及平滑度(r)。试验表明，单纯依靠某一种或两种指标，可能会出现不一致的现象[26]。文献[26]提出通过变异系数定权将归一化后的均方根误差与平滑度线性组合的复合指标，以评价小波降噪的效果。本文在文献[26]的基础上，采用变异系数定权方法将上述4种传统评价指标进行融合。4种传统评价指标的计算方法可以参考文献[27],本文重点介绍复合指标构建过程中各指标的归一化和定权方法。

各指标的归一化处理如下

(11)

式中，RMSE、SNR、ρ和r分别表示均方根误差、信噪比、相关系数及平滑度；PRMSE(i)、PSNR(i)、Pρ(i)和Pr(i)分别为测站i的各指标经归一化处理的结果；max(·)和min(·)分别为取最大值和最小值操作。

各指标的变异系数定权公式为

(12)

式中，CV=σ/u表示各传统指标的变异系数，σ和u分别表示各归一化指标的方差和均值；W表示各归一化指标在复合指标中的权重。则复合指标T可表示为

T=PRMSE×WPRMSE+PSNR×WPSNR+Pρ×WPρ+Pr×WPr

(13)

根据4个传统指标在数据降噪过程中的意义与特征以及归一化和变异系数定权的计算过程可知，在对比不同降噪方法时，复合指标T越小，表示降噪效果越好。

1.3.2 GNSS坐标时间序列函数模型精度

GNSS坐标时间序列的函数模型[28-29]可以表示为

(14)

式中，y(t)是测站在时刻t下的位移；tR是参考时间；a为截距；v是测站运动速度；nF是频率个数，通常用来拟合季节性变化；sk和φk分别是频率ωk的振幅和相位，本文只考虑周年项(k=1)和半周年项(k=2)；nJ是出现阶跃的次数；bj是阶跃时刻tj时测站在地心笛卡儿坐标中的瞬时位移向量；H是Heaviside阶跃函数；ε(t)为随机过程。

研究表明，白噪声加闪烁噪声的噪声组合模型[6,24]是中国区域测站垂向坐标时间序列的最优噪声模型。此外，为了直观地比较降噪前后GNSS时间序列的噪声振幅变化和降噪算法的效果，降噪前后的噪声模型保持一致[6,24]。因此，本文在白噪声加闪烁噪声的模型下利用最大似然估计法[29]解算降噪前后GNSS垂向坐标时间序列的模型参数，得到测站的速度不确定度和闪烁噪声振幅。通过对比降噪前后速度不确定度和闪烁噪声振幅判断模型精度，属于内符合精度。速度不确定度和闪烁噪声振幅越小，表明降噪效果越好[7,28-29]。GNSS垂向坐标时间序列函数模型的计算过程及相关误差分析的公式参见文献[29]的式(7)、式(15—31)及附录。

2 试验与分析

2.1 GNSS数据及处理

以GNSS为主的CMONOC为我国大陆地壳运动监测提供了高精度、高时空分辨率的基础数据。考虑到观测时长和数据缺失情况，本文选取了CMONOC中1999.16—2020.22年的149个GNSS测站坐标时间序列，其中1999.16—2020.22年共10个测站；2010.32—2020.22年共104个测站；2012.19—2019.75年共16个测站；其余测站的观测时间跨度详见表1。本文采用测站的空间分布和观测时长的详细信息如图1和表1所示。本文所选取的GNSS测站的垂向坐标时间序列是由中国地震局GNSS数据产品服务平台(http:∥www.cgps.ac.cn)所提供。GNSS观测数据经GAMIT基线解算处理、GLOBK约束平差，得到GNSS坐标时间序列[4,30]。

表1 部分GNSS基准站的观测时间Tab.1 Observation times of GNSS reference stations

图1 选取的GNSS基准站的分布和观测时长Fig.1 Locations and observation durations of GNSS continuous stations used in this study

由于GNSS测站设备更换、观测条件和数据处理策略等因素的影响，致使GNSS坐标时间序列存在粗差和数据缺失情况。为了后续测站运动特征分析，在进行降噪前，需要对原始GNSS垂向坐标时间序列进行预处理[29,31]:首先基于四分位距统计量识别和剔除粗差，即当GNSS垂向坐标时间序列中的数值大于3倍的GNSS坐标时间序列的四分位距统计量时，可认为是粗差，予以剔除；然后采用分段三次Hermite插值法对GNSS坐标时间序列中缺失的数据进行补全。

2.2 结果分析

本文对CMONOC的149个测站垂向坐标时间序列采用EMD-HD&WD算法和文献[12—14,16]的方法进行降噪处理，计算复合指标T值、降噪前后GNSS垂向坐标时间序列的运动参数及闪烁噪声振幅。为方便表述，文献[12—14,16]中以CMSE、ET、相关系数和HD为筛选准则的EMD降噪方法及处理结果分别记为EMD-CMSE、EMD-ET、EMD-CORR、EMD-HD，降噪前和本文算法的结果记为raw和EMD-HD&WD。

GNSS坐标时间序列降噪的主要目的是获取更真实的构造与季节性运动信号(趋势项和周期项)，因此本文在对比各降噪方法的效果之前，首先分析了降噪后的构造及季节性信号是否受到影响。从信号降噪的原理出发，数据处理过程并未涉及测站的趋势项，因而本文的方法不会影响稳态的构造运动信号。周年项和半周年项的振幅是衡量季节性运动的重要指标[4,28]。为进一步分析降噪过程是否会对季节性运动信号产生影响，本文计算了降噪前后149个GNSS时间序列的周年项和半周年项的平均振幅(表2)。EMD-CMSE、EMD-CORR、EMD-ET、EMD-HD和EMD-HD&WD所得到的周年、半周年项的振幅均略小于降噪前的结果，但相差不大，因此本文所采用的降噪方法并没有平滑掉测站的地球物理信号。图2为不同降噪方法对HEZJ测站的降噪结果。

图2 不同降噪方法得到的噪声与信号(以HEZJ站为例)Fig.2 The time series of noise and signal by different methods for HEZJ station

表2 GNSS坐标时间序列降噪前后GNSS时间序列的周年项和半周年项的平均振幅Tab.2 The average amplitudes of annual and semi-annual signal on GNSS coordinate time series with different strategies mm

2.2.1 最优筛选准则的确定

图3是不同降噪方法所确定的m和T值。本文根据图3(a)计算了不同筛选准则得到的m所占总阶数的比重。结果表明，EMD-HD得到的比重平均值为0.52，大于0.48(EMD-CMSE)、0.44(EMD-CORR)、0.43(EMD-ET)，说明文献[13]的降噪结果与文献[14]相近，EMD-HD所确定的m较其他准则整体偏大。表3为图3(b)中各种降噪方法所得到的T值的平均值。从表3可知，EMD-CMSE、EMD-CORR、EMD-HD和EMD-ET的149个T值的平均值分别是0.61、0.85、0.52、0.59。结合图3(a)和表3可以发现，EMD-CORR与EMD-ET的降噪结果相近，但是T值却相差0.24，原因是因为前者需要预先设置阈值，但该阈值仅对部分测站比较有效[12]；EMD-ET确定为噪声的IMFs分量整体小于EMD-CMSE，但两者的T值却接近，这是由于比较的T值是149个测站数据的平均值，而EMD-ET对于部分测站具有较好的作用[11,13]，改善了该方法的T值平均值；从m的比重来看，EMD-CMSE、EMD-CORR和EMD-ET所识别的噪声的IMFs分量小于EMD-HD的结果，与T值的结果一致，说明前3种方法降噪不足。以上分析表明：EMD-CORR和EMD-ET仅适用于部分测站；HD优于其他筛选准则。

表3 基于不同准则得到的T值平均值Tab.3 Averages of T obtained by different criteria

图3 不同降噪方法所确定的m和T值Fig.3 m and T obtained by different methods

文献[23,25]表明，白噪声可以被各类方法很好地识别并剔除，但对闪烁噪声的剔除效果却不一致。图4是降噪前后GNSS垂向坐标时间序列的闪烁噪声及相对于降噪前的增益率，并给出149个GNSS测站时序在降噪前后的闪烁噪声的平均值，见表4。结果表明，与降噪前的GNSS坐标时间序列相比，经不同方法降噪处理后的GNSS坐标时间序列的闪烁噪声振幅得到显著改善，均小于5 mm·a-0.25。其中，EMD-HD的平均改正率为79.7%，增益率大于85%的测站数占比为54.4%。表4的结果与图3(a)相一致，表明HD的增益效果均高于其他准则的结果。

图4 降噪前后GNSS坐标时间序列的闪烁噪声振幅及相对于降噪前的增益率Fig.4 Amplitudes of flicker noise of GNSS coordinate time series before and after noise reduction and its improvement rates gain rate relative to before noise reduction

图5给出了降噪前后GNSS垂向坐标时间序列的速度不确定度及相对于降噪前的增益率，以进一步分析本文算法对GNSS测站的运动模型的优化效果。表4是降噪前后149个GNSS测站的速度不确定度的平均值。由表4可知，速度不确定度从降噪前的0.54 mm·a-1减少到0.200 mm·a-1以下。由图5可以得出，EMD-HD的平均改正率为79.7%，增益效果均高于其他准则。通过对比图5和图4的结果可知，不同降噪方法对闪烁噪声振幅与速度不确定度的改善效果表现出较强的一致性，说明降噪处理不仅可以剔除GNSS坐标时间序列中的噪声，而且可以提高GNSS测站运动速率的估值精度，两者具有统一性。

表4 降噪前后的GNSS坐标时间序列的闪烁噪声振幅和速度不确定度的平均值Tab.4 Averages of flicker noise amplitude and velocity uncertainty of GNSS coordinate time series before and after noise reduction

图5 降噪前后GNSS垂向坐标时间序列的速度不确定度及相对于降噪前的增益率Fig.5 Amplitudes of velocity uncertainty of GNSS coordinate time series before and after noise reduction and its improvement rates gain rate relative to before noise reduction

从以上分析可以看出，T值与噪声振幅和速度不确定度的结果虽然并非完全一致，但均表明，HD得到的结果好于其他准则的结果，说明HD能够更准确地确定噪声分量。

2.2.2 模态混叠效应的改善

图6是EMD-HD&WD和EMD-HD的闪烁噪声振幅和速度不确定度估值及EMD-HD&WD相对于EMD-HD的增益率。从表4和图6可知，EMD-HD&WD方法降噪后的速度不确定度和闪烁噪声振幅的平均值分别是0.065 mm·a-1和1.710 mm·a-0.25，相对于降噪前平均改正率均为88.4%；相对于EMD-HD方法，速度不确定度方面，96%的测站得到增益，平均增益率为40.4%；仅有6个测站有平均12.8%的负增益。闪烁噪声振幅方面，96.6%的测站得到增益，平均增益率为40.1%；仅有5个测站有平均14.8%的负增益。结果表明，EMD-HD&WD算法改善模态混叠效应的同时进一步提高了GNSS坐标时间序列的精度。

图6 EMD-HD&WD和EMD-HD的运动模型参数及EMD-HD&WD相对于EMD-HD的增益率Fig.6 The motion model parameters of EMD-HD&WD and EMD-HD and the gain rate of EMD-HD&WD compared to EMD-HD

3 结 语

本文针对基于EMD的GNSS基准站坐标时间序列的噪声识别与剔除过程中存在的问题，提出了EMD-HD&WD算法，并通过对CMONOC的149个GNSS垂向坐标时间序列的降噪处理验证了该算法的有效性及效果，得出以下结论：

(1) HD优于现有的其他筛选准则。相对于CMSE等准则，EMD-HD和EMD-HD&WD更好地识别和剔除了GNSS垂向坐标时间序列中的噪声,并且降噪处理后测站的速度不确定度和噪声振幅显著小于其他方法，证实了HD作为筛选准则的可靠性和普适性。

(2) EMD-HD&WD方法能较好地改善EMD的模态混叠效应。EMD-HD&WD方法有效地将低频IMF分量中混入的高频信号分离出来，进一步降低了测站的速度不确定度和闪烁噪声，在保留坐标时间序列中季节变化特征的同时具有较好的光滑性，达到更好的降噪效果，一定程度改善了EMD方法的模态混叠效应。

(3) 结合GNSS坐标时间序列分析的理论，本文提出将速度不确定度和噪声振幅作为降噪评价指标，效果良好。经EMD-HD&WD算法对GNSS坐标时间序列进行处理后，噪声振幅和速度不确定度的平均值从降噪前的14.35 mm·a-0.25和0.54 mm·a-1降低至1.71 mm·a-0.25和0.065 mm·a-1,噪声振幅和速度不确定度的平均改正率均为88.4%。本文将传统的降噪方法的评价指标和GNSS测站运动参数进行综合评价降噪效果，为降噪方法的可靠性评价提供新的依据。

(4) 本文所提的降噪算法EMD-HD&WD能够有效识别并剔除GNSS坐标时间序列的噪声，提高GNSS坐标时间序列和垂向速度场的精度，不但有助于提取GNSS坐标时间序列中季节性变化并解释其物理机制，而且还为基于GNSS技术的地球参考框架的建立与精化，以及板块构造运动、地球质量变迁等地球动力学过程等研究提供了高精度数据。