2023, Vol. 43
2. 中国科学院精密测量科学与技术创新研究院大地测量与地球动力学国家重点实验室，武汉市徐东大街340号，430077;
3. 天津市地质工程勘测设计院有限公司，天津市红旗南路261号，300191;
4. 中国地震局第一监测中心，天津市耐火路7号，300180
高精度GNSS坐标时序可为板块运动、地壳变形等研究提供基础数据。GNSS坐标时序具有时域相关性和空间域相关性[1]，时域相关性的分析方法一般包括极大似然估计法(maximum likelihood estimate, MLE)和最小二乘估计法[2]。在区域连续GNSS坐标时序中存在一种空间域误差，即共模误差CME。共模误差处理方法较多[3]，最早使用堆栈滤波法来减弱共模误差，该方法假设CME在区域内均匀分布，但该假设在空间跨度较大的观测网内不成立[4]。PCA可消除CME，弥补堆栈滤波法的不足：刘晓祥等[5]使用PCA对陆态网224个GNSS基准站坐标时序进行分析，结果表明，华北地区、西北地区、云南地区在N、E、U方向上提取的主分量贡献率变化明显，其中U方向变化最显著；明锋等[1]使用PCA和ICA分别对区域网CME进行分析，结果表明，ICA能更有效地提取观测网CME；李斐等[6]分别使用PCA和ICA对南极半岛地区15个GNSS站进行滤波分析，结果表明，ICA的滤波效果优于PCA。
由于PCA可能会剔除测站的原始信息，导致某一主成分包含不同的物理模式，进而出现虚假特征。因此为减少PCA对测站原始信息造成的损失，有效利用CME中的高频信息，本文以华北地区GNSS坐标时序为例，提出将PCA与改进的完全集合经验模态分解法ICEEMDAN相结合的IC-PCA方法，消除共模误差。ICEEMDAN可精确分离出时间序列中的噪声和真实信号[7]，可对PCA滤波剔除的共模误差进行二次分解，提取出其中有用的IMF，有效减少PCA对原始信息造成的损失。利用CME中的高频信息，可使滤波后的结果更加符合原始序列的变化趋势。
1 GNSS数据预处理华北地区具有较为密集的GNSS基准站，且部分站点建站时间较早，坐标时序跨度较长。本文收集华北地区2011-01-04~2021-05-04连续11 a的陆态网络GNSS基准站时间序列，数据来源于中国地震局第一监测中心(https://www.fmac.ac.cn)。由于GNSS基准站原始坐标时序会受到设备变更、天线更换、地震等外界因素干扰，也会因通信链路终端损坏或电源故障等因素随机丢失部分数据[8]，因此本文首先剔除华北地区37个基准站中缺失率大于25%的10个站点，利用剩余27个站点进行分析；然后按照3倍标准差(3σ)准则对GNSS坐标时序进行去粗差处理并将其剔除；最后进行插值处理。本文采用KKF(kriged Kalman filter)[9]方法，借助GMIS(GNSS missing data interpolation)软件对数据进行修复，得到不包含粗差且完整的时间序列。
在提取共模误差前需要检验区域内各测站之间的相关性。本文采用KMO因子分析法[6]对时间序列进行检验，经检验选取的27个GNSS站点在N、E、U三个方向上的KMO值分别为0.977、0.963、0.976，说明原始数据适合进行因子分析。
2 结果分析 2.1 PCA和IC-PCA结果分析在进行PCA滤波前，需要对检验后的27个站点在N、E、U三个方向上的时间序列进行因子分析，因子分析得分较低的测站结果见表 1。
|
表 1 因子分析结果 Tab. 1 The results of factor analysis |
由表 1可知，N方向上PC1因子得分最低的2个测站为NMDW和TJBD，可作为异常站点剔除。同理，在E、U方向上分别剔除SXGX、BJGB和SXXX、SXGX四个站点。可以看到，N、E、U三个方向上PC1分量的因子得分高于PC2和PC3，说明PC1分量更适合提取CME。使用筛选后的测站再次进行PCA滤波处理，得到N、E、U方向上的空间响应(图 1)。
|
图中箭头方向表示主分量的空间响应方向; 箭头长短表示主分量的空间响应度, 越接近于1则空间响应度越高; 箭头向上表示空间响应为正, 向下为负。 图 1 N、E、U方向上的PC1分量对应的空间响应 Fig. 1 SR of PC1 components in N, E, U directions |
由图 1可见，N、E、U三个方向上各测站的空间响应SR基本一致，SR1的空间响应方向相同且空间响应度超过85%，因此本文使用N、E、U三个方向上的PC1计算CME。由于PCA是提取残差序列方差最大的PC计算CME，因此PCA滤波可能会剔除测站的部分原始信息，或残差时序中存在某些未建模误差PCA未被剔除。为减少PCA滤波对原始信息造成的损失，使用ICEEMDAN法对PCA滤波剔除的CME进行二次分解。ICEEMDAN法分解得到的IMF分层能够更好地展现原有坐标时序，避免有效信息的损失。使用排列熵[10](permutation entropy, PE)分辨时序分解后IMF的频率，将选取出的干净信号和第一次PCA滤波后的结果相加得到最终结果。随机选取NMER测站和HEYY测站作为实验对象，结果如图 2所示。
|
图 2 IC-PCA滤波前后结果 Fig. 2 The results before and after regional filter using IC-PCA method |
由图 2可见，经过PCA滤波后，N、E方向上的时间序列与原始序列相比过于平滑，与原始变化趋势不符；使用ICEEMDAN法优化后的曲线更加符合原始序列的变化趋势。NMER测站N、E方向上原始时序在2012~2013年、2015~2016年和2019~2020年间的下降趋势、HEYY测站N方向上原始时序在2012-06~2015-06的下降趋势及E方向上原始序列在2011-04~2012-08的明显上升趋势均得到显著响应。U方向上NMER测站在2020-10~2021-05以及HEYY测站在2011-05~2019-08、2018-07~2019-02出现的PCA滤波时序和原始序列振幅差距较大的问题也得到改善，优化后的结果与原始序列的变化趋势基本相同。对PCA和IC-PCA优化方法进行对比，结果见表 2。
|
|
表 2 部分测站N、E、U方向上PCA及IC-PCA滤波前后RMS变化 Tab. 2 RMS in N, E, U directions before and after PCA and IC-PCA filtering |
由表 2可见，经过ICEEMDAN法优化后，PCA滤波结果残差的RMS得到有效降低，最终的滤波结果也更加符合原始序列的变化趋势。同时可以看出，N方向上BJSH、JIXN两个测站PCA滤波方法的RMS残差分别为56.6%和54.4%，远大于其他测站。造成上述差异的原因主要是BJSH、JIXN测站的原始序列和其他测站相比数值较小，因此与其他测站一起计算的CME不符合其原始序列，出现虚假特征。但经过ICEEMDAN法优化得到的结果更加接近原始序列，图 3为BJSH、JIXN测站N、E方向上的滤波前后结果对比。
|
图 3 N、E方向上的IC-PCA滤波结果 Fig. 3 IC-PCA filtering results in N, E directions |
由图 3可见，PCA滤波结果和原始序列相差较大，而IC-PCA滤波结果和原始序列相差较小，说明IC-PCA滤波方法可以在一定程度上减少PCA造成的虚假特征。由滤波后测站N、E、U三个方向上残差占比的RMS变化可知，IC-PCA方法能够更加精确地提取测站时序中的CME，且原始时序的信息损失更小。
2.2 ICA结果分析为更好地验证IC-PCA方法的有效性，使用ICA滤波方法对原始序列进行处理，并对滤波结果进行分析。第一次进行ICA滤波时未剔除BJSH和JIXN测站，图 4为随机选取的HELQ测站ICA滤波结果。
|
图 4 HELQ测站第一次ICA滤波结果 Fig. 4 The results of the first ICA filter of HELQ station |
由图 4可见，N、E方向上的ICA滤波结果与原始序列的差异较大，推测ICA滤波可能受到BJSH和JIXN测站较小的原始序列数据影响，拟合效果较差。因此去除BJSH和JIXN测站后再次进行ICA滤波处理，根据N、E、U三个方向上所有独立分量对应的SR进行分类，挑选相同SR方向的独立成分IC。在N、E方向上使用14个独立分量计算CME，在U方向上使用7个独立分量计算CME。图 5为N、E、U方向上挑选的IC分量对应的SR图。
|
图 5 N、E、U方向上IC分量对应的空间响应 Fig. 5 SR of IC component in N, E, U directions |
由图 5可见，各方向挑选的独立分量所对应的空间响应方向均一致，且空间响应度较高，可由此计算各方向的CME。为方便比较，以NMER、HEYY测站为例进行ICA滤波处理，结果见图 6。
|
图 6 第二次ICA滤波结果 Fig. 6 The second ICA filtering results |
由图 6可见，对测站进行ICA滤波后的结果与原始序列基本相同，仅在个别时间段内存在偏差。同样使用RMS比较ICA滤波前后的变化，表 3为各测站ICA滤波前后RMS变化情况。
|
|
表 3 部分测站N、E、U方向上ICA滤波前后的RMS变化 Tab. 3 RMS in N, E, U directions before and after ICA filtering |
由表 3可见，ICA滤波在N、E方向上的残差出现0和负值，说明N、E方向上的滤波结果相较于原始序列无变化或者滤波后CME增加，因此N、E方向上的ICA滤波结果无法采用。U方向上ICA滤波结果的平均残差相较于PCA滤波结果减小1.887%，相较于IC-PCA滤波结果增加1.423%，说明IC-PCA滤波具有较好的优化效果。ICA滤波方法容易受到原始序列数值异常的影响，滤波后的结果与原始序列变化趋势相差较大。由此可知，在华北地区，ICA去除区域性共模误差的效果不如IC-PCA。
3 结语1) PCA滤波后BJSH、JIXN两个测站在N、E方向上均出现虚假特征，滤波结果与原始序列相差较大；IC-PCA滤波方法可有效减少PCA滤波造成的虚假特征，滤波结果与原始序列接近。
2) IC-PCA滤波后N、E、U三个方向上的残差RMS分别降低1.93%、1.92%和7.60%，说明IC-PCA方法可以减少PCA方法对原始序列造成的信息损失，且相较于ICA也更加符合原始序列的变化趋势。
| [1] |
明锋, 杨元喜, 曾安敏. 共模误差PCA与ICA提取方法的比较[J]. 大地测量与地球动力学, 2017, 37(4): 385-389 (Ming Feng, Yang Yuanxi, Zeng Anmin. Analysis and Comparison of Common Mode Error Extraction Using Principal Component Analysis and Independent Component Analysis[J]. Journal of Geodesy and Geodynamics, 2017, 37(4): 385-389)
(  0) |
| [2] |
姜卫平, 王锴华, 李昭, 等. GNSS坐标时间序列分析理论与方法及展望[J]. 武汉大学学报: 信息科学版, 2018, 43(12): 2 112-2 123 (Jiang Weiping, Wang Kaihua, Li Zhao, et al. Prospect and Theory of GNSS Coordinate Time Series Analysis[J]. Geomatics and Information Science of Wuhan University, 2018, 43(12): 2 112-2 123)
( 0) |
| [3] |
Freymueller J T. Seasonal Position Variations and Regional Reference Frame Realization[M]. Berlin: Springer, 2009
( 0) |
| [4] |
Chen Q, Dam T, Sneeuw N, et al. Singular Spectrum Analysis for Modeling Seasonal Signals from GPS Time Series[J]. Journal of Geodynamics, 2013, 72: 25-35 DOI:10.1016/j.jog.2013.05.005
( 0) |
| [5] |
刘晓祥, 高二涛, 罗益, 等. 利用主成分分析法分析GNSS坐标时间序列[J]. 大地测量与地球动力学, 2021, 41(1): 43-48 (Liu Xiaoxiang, Gao Ertao, Luo Yi, et al. Analysis of Coordinate Time Series of CMONOC GNSS Fiducial Stations Using Principal Component Analysis[J]. Journal of Geodesy and Geodynamics, 2021, 41(1): 43-48)
( 0) |
| [6] |
李斐, 李文浩, 张胜凯, 等. 基于独立分量分析的南极半岛GNSS网区域滤波[J]. 地球物理学报, 2019, 62(9): 3 279-3 295 (Li Fei, Li Wenhao, Zhang Shengkai, et al. Spatiotemporal Filtering for Regional GNSS Network in Antarctic Peninsula Using Independent Component Analysis[J]. Chinese Journal of Geophysics, 2019, 62(9): 3 279-3 295)
( 0) |
| [7] |
Colominas M A, Schlotthauer G, Torres M E. Improved Complete Ensemble EMD: A Suitable Tool for Biomedical Signal Processing[J]. Biomedical Signal Processing and Control, 2014, 14: 19-29 DOI:10.1016/j.bspc.2014.06.009
( 0) |
| [8] |
Gazeaux J, Williams S, King M, et al. Detecting Offsets in GPS Time Series: First Results from the Detection of Offsets in GPS Experiment[J]. Journal of Geophysical Research: Solid Earth, 2013, 118(5): 2 397-2 407 DOI:10.1002/jgrb.50152
( 0) |
| [9] |
Liu N, Dai W J, Santerre R, et al. A MATLAB-Based Kriged Kalman Filter Software for Interpolating Missing Data in GNSS Coordinate Time Series[J]. GPS Solutions, 2018, 22(1): 25 DOI:10.1007/s10291-017-0689-3
( 0) |
| [10] |
Dong D N, Fang P, Bock Y, et al. Spatiotemporal Filtering Using Principal Component Analysis and Karhunen-Loeve Expansion Approaches for Regional GPS Network Analysis[J]. Journal of Geophysical Research: Solid Earth, 2006, 111(B3)
( 0) |
2. State Key Laboratory of Geodesy and Earth's Dynamics, Innovation Academy for Precision Measurement Science and Technology, CAS, 340 Xudong Street, Wuhan 430077, China;
3. Tianjin Geological Engineering Survey and Design Institute Co Ltd, 261 South-Hongqi Road, Tianjin 300191, China;
4. The First Monitoring and Application Center, CEA, 7 Naihuo Road, Tianjin 300180, China