GNSS coordinate time series denoising analysis combined with weighted wavelet and EEMD
-
摘要: 针对GNSS坐标时间序列中有用信号与噪声难以准确分离这一问题，本文提出加权小波Z变换(weighted wavelet Z-transform, WWZ)和集合经验模态分解(ensemble empirical mode decomposition，EEMD)的降噪方法. 通过对西北地区70个陆态网络连续站垂向坐标时间序列的降噪处理，分别采用均方根误差(root mean squared error，RMSE)、信噪比(signal to noise ratio，SNR)、闪烁噪声(flicker noise，FN)振幅及速度不确定度为评价指标，验证了本文方法的降噪效果在一定程度上优于小波降噪和EEMD降噪. 结果显示：WWZ-EEMD相比小波降噪和EEMD降噪，降噪后信号序列RMSE分别降低了0.331 mm、 0.757 mm，SNR分别提高了1.911 dB、3.635 dB；FN振幅及速度不确定度均有明显改善，验证了本文降噪方法的有效性.
-
关键词:
- GNSS坐标时间序列 /
- 降噪 /
- 小波 /
- 集合经验模态分解(EEMD) /
- 速度不确定度
Abstract: Aiming at the problem that it is difficult to accurately separate the useful signal and noise in the GNSS coordinate time series, this paper proposes a noise reduction method based on combined weighted wavelet Z-transform (WWZ) and set empirical mode decomposition (EEMD). Through the noise reduction processing of the vertical coordinate time series of 70 continuous stations in the northwest region, the root mean square error (RMSE), signal-to-noise ratio (SNR), flicker noise (FN) amplitude and velocity uncertainty are used as the evaluation indicators respectively, which verifies that the noise reduction effect of the method in this paper is superior to wavelet noise reduction and EEMD noise reduction to a certain extent. The results show that compared with wavelet denoising and EEMD denoising, the RMSE of signal sequence after denoising is reduced by 0.331 mm and 0.757 mm respectively, and the SNR is increased by 1.911 dB and 3.635 dB respectively; The uncertainty of FN amplitude and velocity has been significantly improved, which verifies the effectiveness of the noise reduction method in this paper. -
0. 引 言
近年来，我国GPS基准网的建设取得了长足的发展，为研究各种不同的地球物理现象积累了大量可靠的观测数据. 由于受外部环境因素的影响，GNSS坐标时间序列中不可避免地存在噪声，极大地影响了测站的真实运动特征[1-2]. 因此，如何更加精确有效地去除GNSS坐标时间序列中的噪声，成为获取更加真实信号的关键[3].
目前已有许多学者对小波分析、经验模态分解(empirical mode decomposition，EMD)、集合经验模态分解(ensemble empirical mode decomposition，EEMD)等降噪方法[4-6]有了深入研究，以上方法均能够对坐标时间序列中噪声起到抑制作用，但均有一定局限性. 小波分析具有多分辨率、去相关性等特点，马俊等[7]利用小波包系数信息熵法对5 a的坐标时间序列进行降噪处理，消除了有色噪声的影响；Ji等[8]提出了一种加权小波分析，对27个测站进行信号提取，结果表明该方法比传统小波分析所提取的信号要更接近真实信号. 而针对EMD会产生模态混叠这一现象，Wu等[9]在其基础上提出了EEMD，极大地削弱了模态混叠现象；张恒璟等[10]采用EEMD对BJFS站和URUM站10 a的高程数据进行处理，结果显示相比EMD，EEMD更好地提取坐标时间序列中的周期项和趋势项.
为能够更准确地提取坐标时间序列中的有用信息，本文基于以上理论方法提出了一种加权小波Z变换(weighted wavelet Z-transform，WWZ)与EEMD相结合的降噪方法，即WWZ-EEMD降噪方法. 通过对西北地区70个测站的降噪处理，验证了本文方法能够有效地提取噪声分量中的有用信息，从而改善了降噪效果.
1. 原理及方法
1.1 加权小波多尺度分析
小波分析具有极好的时频特性和多尺度分解的能力，传统的小波处理方法没有对日解误差进行考虑，然而GPS坐标时间序列的精度是随时间而变化的，因此利用日解误差来构造权重因子，通过不断地迭代来更加有效地提取时间序列中的信号[8].
首先对原始坐标时间序列
$ x $ 进行小波多尺度分解，采用相关系数法进行信噪分离，则在时间序列的第$ i $ 个历元时，信号$ {s_{{i}}} $ 和噪声$ {e_{{i}}} $ 可表示为[8]$$ x_i=s_i+e_i,\left(i=0,1,2\cdots,N-1\right) $$ (1) 由于GPS坐标时间序列中没有考虑时间相关性，日解观测误差只与噪声有关，则将噪声与权重因子相乘构造出新的时间序列为[8]
$$ x'_{{i}} = {s_{{i}}} + {p_{{i}}}{e_{{i}}} $$ (2) 权重因子
$ {p_{{i}}} $ 可通过日解观测误差来得到[8]$$ {p_{{i}}} = \sqrt {\frac{N}{{\displaystyle\sum\limits_{j = 0}^{N - 1} {\sigma _{{i}}^{{2}}/\sigma _{{j}}^{{2}}} }}} $$ (3) 式中：
$ N $ 为观测数；$ \sigma _{{i}}^{{2}} $ 为原始序列在第$ i $ 个历元时观测误差的平方. 经过一次小波分解与重构之后，权重因子也不断更新[8]$$ {\sigma _{{{x'}_{{i}}}}} = {p_{{i}}}{\sigma _{{i}}} $$ (4) 利用小波连续的分解和重构，提取到信号和噪声序列，当相邻两次信号差小于阈值0.001时，则迭代完成，具体迭代过程可参考文献[8].
1.2 EEMD
考虑到EMD模态混叠的问题，EEMD通过在分解前在原始信号中添加白噪声(white noise，WN)构成新信号的方法来改善模态混叠的问题，其分解步骤与EMD一致，再将多次分解的相同尺度下的本征模态函数(intrinsic mode function，IMF)分量取平均，作为该次分解得到的IMF分量，最终得到由EEMD得到的不同尺度下的IMF分量[11]. 按照相关准则进行信号重构，可以得到处理后的坐标时间序列.
$$ {X_{{i}}}\left( t \right) = X\left( t \right) + {w_{{i}}}\left( t \right)$$ (5) 式中，
$ i $ 表示第$ i $ 次分解，其余步骤与EMD相同 ，可参考文献[12].1.3 联合算法原理
对原始坐标序列进行加权小波多尺度分解，加权迭代处理后得到噪声序列和信号序列，具体过程如下，对原始坐标时间序列
$ X=\left\{{X}_{1},{X}_{2},\cdots, {X}_{n}\right\} $ 进行加权小波多尺度分解和重构，则得到$$ X = {S_1} + {n_1}$$ (6) 式中：
$ X $ 为原始坐标时间序列；$ {S_1} $ 为重构信号序列；$ {n_1} $ 为重构噪声序列.对得到的噪声序列再进行EEMD，得到不同尺度下的IMF分量，利用相关系数法重构信号
$ {S_2} $ 和噪声$ {n_2} $ 序列，将两次算法所拟合的信号序列累加则得到降噪后的信号$ \hat X $ 为$$ \hat X = {S_1} + {S_2} $$ (7) 2. 数据获取与与预处理
本文数据来自中国地震局GNSS数据产品服务中心平台，其解算策略为采用GAMIT/GLOBK软件，结合众多参数进行平差来获得原始坐标时间序列. 选取西北地区70个站点的垂向坐标时间序列，时间跨度为2011—2021共10 a的原始数据进行研究，站点分布如图1所示. 由于GNSS测站数据解算策略、观测条件等影响，导致坐标时间序列中存在粗差和数据缺失的情况[13]. 为了更好分析测站的运动特征，应该对原始坐标时间序列进行预处理：首先利用四分位距法探测并剔除粗差，即当坐标时间序列中观测值大于四分位距统计量的3倍时，应当予以剔除；采用最小二乘法对原始坐标时间序列进行线性趋势项的拟合来去除趋势项；为保证数据均匀且连续，采用基于克里金卡尔曼滤波模型的动态时空插值方法[14]对去趋势序列进行插值处理.
![]()
图 1 西北地区陆态网络连续站站点分布图Figure 1. Distribution map of continuous stations of land network in northwest China2.1 基于主成分分析的时空滤波
目前，共模误差是一种局部尺度的误差已无争议，但其来源与本质尚未定论，研究显示，滤波后坐标时间序列中仍含有除WN外的其他噪声[15]. 本文采用最常用的主成分分析(principal component analysis，PCA)滤波方法[16]进行共模误差的剔除. PCA滤波作为一种线性降维的方法，可以极大程度的来保证滤波后数据的特征. 经统计西北地区70个测站滤波后U分量的残差时间序列均方根(root mean square，RMS)有所减少，平均降低了40.7%，说明较好地去除了站点共模误差，提高了数据精度.
经过以上数据处理，极大程度地去除了粗差和共模误差等影响，为后面降噪处理提供了更加有利的数据基础.
3. 实验结果与分析
3.1 降噪结果及对比分析
GNSS坐标时间序列中包含多种类型的噪声及信号，具有较强的时频特征. 降噪的主要目的是能够更加准确地分析季节性运动信号及构造运动，小波分析与EEMD均具有良好的多尺度时频分析特征，因此组合两者优势对西北地区GNSS坐标时间序列进行降噪分析.
利用加权小波分析对处理后的坐标时间序列进行多尺度分解，但在在分解过程中，小波基函数的选取和分解层数的确定也尤为重要[17]. 为此对常用的处理时间序列较有优势的小波基函数及其特点进行了统计，结果如表1所示. 由表1可知，CoifN小波函数具有良好的时间序列分析特性，因此，本文一般选择Coif5小波进行小波分解. 根据小波分解的特性，随着分解层数的增加，坐标时间序列各层信号被划分越来越细，但计算效率也随之增加. 以GSJN站垂向坐标时间序列为例，分别进行不同层数的分解(5~8层)，如表2所示，通过对比5~8不同分解层数的RMSE和MAE值，本实验选择分解层数为7层.
表 1 常用小波基函数及其特点Table 1. Common wavelet basis functions and their characteristics小波基函数 支撑长度 消失矩阶数 对称性 特点 Haar小波 1 1 对称 时域不连续，
频域局部性差dbN小波 2N N 近似对称 随着N的增加
光滑性也提高CoifN小波 6N 2N 近似对称 具有比dbN更好
的对称性SymN小波 2N N 近似对称 更适合图像
处理领域表 2 不同分解层数的拟合精度对比Table 2. Comparison of fitting accuracy of different decomposition levelsmm 拟合精度 分解层数 5 6 7 8 RMSE 0.657 0.615 0.614 0.614 MAE 0.252 0.240 0.239 0.240 图2为小波7层分解后重构的各层系数，由图2可知，在进行7层小波分解后高频系数愈发平稳，在细节系数d1、d2、d3、d4中表现出的随即项成分较多，周期项较少，其相应的拟合效果要小于其他层.
![]()
图 2 GSJN站各层重构小波系数Figure 2. Reconstructed wavelet coefficients of each layer of GSJN station采用相关系数法进行信噪分离，将局部相关系数最小的一层定为边界层. 由表3可知，对于垂向分量，信号与噪声之间的边界层为D6. 将边界层之前的(含边界层)重构为噪声序列，之后的重构为信号序列，图3为GSJN站U分量加权小波重构后的信号及噪声序列. 经过加权小波所拟合的时间序列更能反应原始时间序列的走势.
表 3 GSJN站各层重构分量相关系数Table 3. Correlation coefficient of reconstruction components of each layer of GSJN station分解层数 相关系数 D1 0.3819 D2 0.3134 D3 0.2755 D4 0.2679 D5 0.2465 D6 0.2278 D7 0.2857 ![]()
图 3 GSJN站U分量加权小波重构后的信号及噪声序列Figure 3. Signal and noise sequence after GSJN station U component weighted wavelet reconstruction尽管小波分解对原始信号进行了有效的信号提取，但噪声序列中难免仍会存在一些有用的周期项或季节项信号，因此对重构后的噪声序列采用EEMD，添加噪声参数设置为0.2，重复次数设为300，大部分站点序列分解的本征模态分量IMF为12个，SNHY站分解的分量为11个，其分解的IMF均包含一列原始信号、其他高低频模态分量及趋势项，根据相关系数准则重构信号与噪声序列. 图3为GSJN站噪声序列中提取的信号和噪声，结果表现出所提取的信号序列中仍有较强的规律性，主要表现为周期性.
通过以上两次提取，将所拟合的信号序列重构叠加后得到最终降噪后的信号，如图4(b)所示，对于长度为10 a的坐标时间序列，经过本文方法降噪后的时间序列不仅能反映出原始时间序列的走势，而且两者之间差异明显减少. 经统计，本文方法、小波降噪及EEMD降噪后GSJN站所获得的残差平均振幅分别为0.5996 mm、0.7836 mm、0.9574 mm，较单一方法均有明显降低，结果表明本文方法可以更加有效地提取到坐标时间序列中的有用信息，进而极大地提高了噪声中有用信息的利用率.
![]()
图 4 GSJN站噪声序列中提取的信号和噪声Figure 4. Signal and noise extracted from GSJN station noise sequence为定量说明本文方法和小波分析及EEMD的降噪效果，分别计算不同方法下降噪结果的信噪比(signal to noise ratio，SNR)和RMSE. SNR值越大，表示信号中噪声含量少，降噪效果越好. RMSE则反映降噪后信号与真实信号的偏差，对于非线性的GNSS坐标时间序列，RMSE越小，降噪效果越好. 统计结果如图5所示.
![]()
图 5 各测站不同方法降噪后RMSE和SNRFigure 5. RMSE and SNR after noise reduction by different methods at each station由图5可知，本文方法对各测站U分量降噪后的RMSE和SNR均有所改善，其RMSE明显降低，SNR有较大提高，其原因在于对小波分解重构后的噪声序列进行了二次提取，极大地避免了噪声中有用信息的丢失. 经统计在对西北地区70个测站的降噪处理中，本文方法降噪效果均优于EEMD，相比小波降噪，有56个站点为本文方法降噪效果最佳，GSDX、GSMX、GSWD、QHDL、QHTT、QHYS、SNHY、SNXY、XJBE、XJDS、XJML、XJTC、XJWL、XJWQ这14个站小波降噪效果略好，但与本文方法降噪效果相差不大. 综合来看，本文方法降噪后信号序列的RMSE均值为1.469 mm，对比小波降噪和EEMD降噪，分别降低了0.331 mm、 0.757 mm；本文方法SNR均值为7.758 dB，分别提高了1.911 dB、3.635 dB，结果表明：本文方法的降噪效果在一定程度上优于上述两种方法.
3.2 闪烁噪声与速度不确定度
为进一步验证本文降噪效果，将闪烁噪声(flicker noise, FN)振幅和测站速度不确定度作为精度指标[18]，采用统一的WN+FN组合模型分别对降噪后的数据进行测站速度估计. FN振幅以及速度不确定度平均值统计结果如表所示.
由表4可知，在WN+FN组合模型下，站点速度不确定度由降噪前0.808 mm/a降低到0.100 mm/a以下，表明小波降噪、EEMD降噪以及本文方法降噪对坐标时间序列噪声都有不同程度的剔除效果. 通过对比分析，不同降噪方法对FN振幅及速度不确定都的改善效果具有较强的一致性. 对于FN，改善效果由高到低依次为：本文联合降噪方法、小波降噪、EEMD降噪. 对于站点速度不确定度，本文方法改善效果较好，相对于小波降噪和EEMD降噪分别平均降低了31.6%、43.8%. 结果显示本文联合加权小波及EEMD的组合方法降噪效果最佳.
表 4 不同方法降噪后FN及速度不确定度平均值Table 4. Average value of flicker noise and speed uncertainty after noise reduction by different methods处理方法 原始数据 预处理数据 小波 EEMD 本文方法 FN/(mm.a−0.25) 18.242 9.286 1.355 1.803 0.894 速度不确定/(mm.a−1) 0.808 0.419 0.060 0.073 0.041 4. 结 语
本文利用EEMD分解对噪声序列中有用信号进行再次重构提取，提出了一种WWZ-EEMD的组合降噪方法. 通过对研究区域70个测站垂向坐标时间序列的降噪处理，结果表明：本文降噪方法相比传统方法，能够更有效地提取高频噪声信息，进而提高去噪精度，具有更低的FN和较高的SNR；将FN振幅和速度不确定度作为评价指标，结果显示本文方法可以更加有效地降低FN振幅以及速度不确定度，进一步验证了本文方法降噪效果.
-
![]()
图 1 西北地区陆态网络连续站站点分布图
Figure 1. Distribution map of continuous stations of land network in northwest China
![]()
图 2 GSJN站各层重构小波系数
Figure 2. Reconstructed wavelet coefficients of each layer of GSJN station
![]()
图 3 GSJN站U分量加权小波重构后的信号及噪声序列
Figure 3. Signal and noise sequence after GSJN station U component weighted wavelet reconstruction
![]()
图 4 GSJN站噪声序列中提取的信号和噪声
Figure 4. Signal and noise extracted from GSJN station noise sequence
![]()
图 5 各测站不同方法降噪后RMSE和SNR
Figure 5. RMSE and SNR after noise reduction by different methods at each station
表 1 常用小波基函数及其特点
Table 1 Common wavelet basis functions and their characteristics
小波基函数 支撑长度 消失矩阶数 对称性 特点 Haar小波 1 1 对称 时域不连续，
频域局部性差dbN小波 2N N 近似对称 随着N的增加
光滑性也提高CoifN小波 6N 2N 近似对称 具有比dbN更好
的对称性SymN小波 2N N 近似对称 更适合图像
处理领域
下载: 导出CSV
表 2 不同分解层数的拟合精度对比
Table 2 Comparison of fitting accuracy of different decomposition levels
mm 拟合精度 分解层数 5 6 7 8 RMSE 0.657 0.615 0.614 0.614 MAE 0.252 0.240 0.239 0.240
下载: 导出CSV
表 3 GSJN站各层重构分量相关系数
Table 3 Correlation coefficient of reconstruction components of each layer of GSJN station
分解层数 相关系数 D1 0.3819 D2 0.3134 D3 0.2755 D4 0.2679 D5 0.2465 D6 0.2278 D7 0.2857
下载: 导出CSV
表 4 不同方法降噪后FN及速度不确定度平均值
Table 4 Average value of flicker noise and speed uncertainty after noise reduction by different methods
处理方法 原始数据 预处理数据 小波 EEMD 本文方法 FN/(mm.a−0.25) 18.242 9.286 1.355 1.803 0.894 速度不确定/(mm.a−1) 0.808 0.419 0.060 0.073 0.041
下载: 导出CSV
-
[1] 陈祥, 杨志强, 田镇, 等. GA-VMD与多尺度排列熵结合的GNSS坐标时序降噪方法[J]. 武汉大学学报(信息科学版), 2023, 48(9): 1425-1434. [2] 李昭, 姜卫平, 刘鸿飞, 等. 中国区域IGS基准站坐标时间序列噪声模型建立与分析[J]. 测绘学报, 2012, 41(4): 496-503. [3] 张恒璟, 龙安森, 文汉江. EEMDAN的CORS站高程时间序列分析方法[J]. 测绘科学, 2020, 45(2): 29-34. [4] 戴海亮, 孙付平, 姜卫平, 等. 小波多尺度分解和奇异谱分析在GNSS站坐标时间序列分析中的应用[J]. 武汉大学学报(信息科学版), 2021, 46(3): 371-380. [5] 姜卫平, 王锴华, 李昭, 等. GNSS坐标时间序列分析理论与方法及展望[J]. 武汉大学学报(信息科学版), 2018, 43(12): 2112-2123. [6] 张双成, 何月帆, 李振宇, 等. EMD用于GPS时间序列降噪分析[J]. 大地测量与地球动力学, 2017, 37(12): 1248-1252. [7] 马俊, 曹成度, 姜卫平, 等. 利用小波包系数信息熵去除GNSS站坐标时间序列有色噪声[J]. 武汉大学学报(信息科学版), 2021, 46(9): 1309-1317. [8] JI K P, SHEN Y Z, WANG F W. Signal extraction from GNSS position time series using weighted wavelet analysis[J]. Remote sensing, 2020, 12(6): 992. DOI: 10.3390/rs12060992
[9] WU Z H, HUANG N. Ensemble empirical mode decomposition: a noise-assisted data analysis method[J]. Advances in adaptive data analysis, 2009, 1(1): 1-41. DOI: 10.1142/S1793536909000047
[10] 张恒璟, 程鹏飞. 基于经验模式分解的CORS站高程时间序列分析[J]. 大地测量与地球动力学, 2012, 32(3): 129-134. [11] 范小猛, 胡川, 张重阳, 等. 三种GNSS高程时序降噪方法的效果对比分析[J]. 全球定位系统, 2022, 47(1): 68-73. [12] HUANG N, SHEN Z, LONG S R, et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[J]. Proceedings of the royal society A, 1998, 454(1971): 903-995. DOI: 10.1098/rspa.1998.0193
[13] 黄立人. GPS基准站坐标分量时间序列的噪声特性分析[J]. 大地测量与地球动力学, 2006, 26(2): 31-33,38. [14] SHUMWAY R, STOFFER D S. An apprpach to time series smoothing and forecasting using the EM algorithm[J]. Journal of time series analysis, 1982, 3(4): 253-264. DOI: 10.1111/J.1467-9892.1982.TB00349.X
[15] 王健, 许安安, 周伯烨. 顾及共模误差的大区域GPS网坐标时间序列噪声分析[J]. 测绘通报, 2018(4): 6-9,56. [16] 殷海涛, 甘卫军, 熊永良, 等. PCA空间滤波在高频GPS定位中的应用研究[J]. 武汉大学学报(信息科学版), 2011, 36(7): 825-829. [17] 邱小梦, 陶国强, 王奉伟, 等. LMD和小波阈值的GNSS坐标时间序列降噪应用[J]. 测绘科学, 2021, 46(8): 28-32,48. [18] 杨兵, 杨志强, 田镇, 等. 联合EMD-HD和小波分解的GNSS坐标时间序列降噪分析[J]. 测绘学报, 2022, 51(9): 1881-1889. -
期刊类型引用(2)
1. 王辰晨，徐工. 时变周期信号特征提取及精度分析. 信息记录材料. 2024(06): 242-245 . 
百度学术
2. 田斌，孙州. 优化ICEEMDAN的工频磁异常信号去噪方法. 自动化与仪表. 2024(11): 55-59+68 . 百度学术
其他类型引用(0)
计量
- 文章访问数: 195
- HTML全文浏览量: 65
- PDF下载量: 30
- 被引次数: 2
