2024, Vol. 44
2. 长沙理工大学交通运输工程学院，长沙市万家丽南路二段960号，410114
区域GNSS坐标时间序列通常存在具有空间相关性的误差，即共模误差(common mode error，CME)[1]。一般认为CME由未准确建模的负荷形变残余信号和GNSS技术性及解算误差组成。准确估计和消除CME对提高GNSS数据质量和坐标速率估计精度至关重要[2-3]。为此，学者们提出一些时空滤波方法来探测CME，如区域叠加滤波法(stacking)[1]、区域加权叠加滤波法[4]、主成分分析法(PCA)、Karhunen-Loeve扩展法(KLE)[5]和独立成分分析法(ICA)[6]等。然而，在大尺度区域的时空滤波中，时空滤波结果容易受到个别异常单站、不同地理因素和测站之间距离的影响，增加了CME估计的复杂性[2]。研究表明，CME在大尺度区域上是非空间均匀分布的，尤其在间隔距离超过500 km后，GNSS参考站间的相关性大幅度减弱[7]。在大尺度区域GNSS坐标序列PCA/ICA时空滤波中，很多分量很难表现出所有站点空间响应的一致分布，导致共模分量判断困难[8]。一些研究将大尺度区域划分为若干个子区域进行滤波分析，但子区域的划分通常依赖主观判断或第三方数据产品，并非依据原始观测数据本身，普适性不强[9]。
针对CME估计的复杂性和非均匀分布，本文提出一种自适应时空滤波方法，自适应地将大尺度区域细分为多个子区域，以确保滤波结果的可靠性。通过对陆态网GNSS站垂向序列数据进行滤波分析，验证本文方法的应用效果。
1 数据来源本文实验数据来源于中国地震局提供的陆态网261个参考站的原始GNSS坐标时间序列，时间跨度为2011-01-01~2017-12-31。按照数据完整率大于90%的要求进行筛选，并剔除数据质量较差的时间序列，最终使用数据质量较好的184个GNSS观测站作为本次实验对象。
2 研究方法 2.1 PCA时空滤波方法假设研究区域中有n个GNSS观测站，每个观测站的观测历元为m，构成m×n维的原始GNSS坐标时间序列矩阵A，且m＞n。
对GNSS坐标残差时间序列进行预处理，然后采用PCA时空滤波方法[2]提取共模误差，k个主分量计算得到的共模误差εPCA为：
| $ \varepsilon^{\mathrm{PCA}}=\sum\limits_{k=1}^r v_k p_k, 1 \leqslant r <n $ | (1) |
式中，
本文自适应时空滤波方法的步骤如下：
1) 数据的准备与预处理。首先，利用3倍中误差法对原时间序列进行粗差探测及去除。利用式(2)拟合GNSS坐标时间序列[10]：
| $ \begin{gathered} y\left(t_i\right)=a+b t_i+c \sin \left(2 \pi t_i\right)+d \cos \left(2 \pi t_i\right)+ \\ e \sin \left(4 \pi t_i\right)+f \cos \left(4 \pi t_i\right)+ \\ \sum\limits_{i=1}^{n_g} g_i H\left(t_i-T_{g_i}\right)+\varepsilon_i \end{gathered} $ | (2) |
式中，
2) 拟合后，将常数项、趋势项和周期项从模型中剔除，对部分由于天线移动、板块构造运动、人为干扰等原因造成的阶跃现象进行改正。
3) 自适应划分子区域。对经过预处理的大尺度区域GNSS坐标残差时间序列，随机选取大尺度区域内不同的参考站作为初始站点，通过迭代逐个计算各站点间的空间距离和Pearson相关系数。通过设定阈值来约束测站间距离和Pearson相关系数，将其自适应地划分为不同的子区域。参考以往学者在GNSS坐标时空滤波研究中的结果，设置约束距离为450 km，Pearson相关系数值为0.35 [11-12]。
4) 使用PCA时空滤波方法依次提取各子区域的GNSS坐标残差时间序列的CME[13]。
3 结果与讨论 3.1 PCA时空滤波结果为直观表示时间分量在不同空间分布下的空间响应幅度，将空间响应的幅值缩放区间设置为-100%~100%，相应缩放量乘入对应的时间分量里。图 1为使用PCA时空滤波方法对整块大尺度区域经过数据预处理的GNSS残差时间序列提取的前3个时间分量及其对应的空间分量结果。实验结果显示，PC1对应的特征值累积贡献率达到54%，且其对应的空间响应表现出良好的空间一致性，其对应的平均归一化空间响应值(绝对值)为74.8%。PC2和PC3没有显示出明显的空间一致性。
|
图 1 整体PCA时空滤波时间与空间分量结果 Fig. 1 Time and space component results of overall PCA spatiotemporal filtering |
在整体PCA时空滤波实验中，最终选取第1个时间分量代表整块区域的CME。将CME从原始残差序列中剔除后，184个GNSS参考站坐标残差时间序列平均RMS值由原来的7.3 mm减少为4.8 mm，滤波前后残差时间序列的平均RMS值减少34.2%。图 2为使用整体PCA时空滤波时，各站点序列的RMS减少率。可以看出，传统的整体PCA时空滤波对非均匀误差的削弱效果有限，只对中国中部地区效果较好。
|
图 2 整体PCA时空滤波RMS减少率 Fig. 2 RMS reduction ratio of the overall PCA spatiotemporal filtering |
采用本文自适应PCA滤波方法对相同站点坐标残差时间序列数据进行滤波分析。图 3为随机选取的3个初始站点所得的3组分区结果。3组实验中，分别将整块区域自适应地划分为13、12、14个子区域。分区过程中，个别站点与距离范围内的其他站点的Pearson相关系数较低，没有划分到任一子区域中，在后续滤波中相关站点不作处理。最终在每组随机实验中分别留下了满足约束条件的177、175、178个GNSS站，均保留了90%以上的原始站点数。
|
图 3 自适应分区的3组子区域结果 Fig. 3 Three sets of subregion results for adaptive partitioning |
在自适应划分子区域之后，分组对每一个子区域进行PCA时空滤波。结果显示，每个子区域的PC1与其对应的空间响应表现出良好的空间一致性，因此选择自适应划分得到的子区域的第1个主分量作为其对应的共模误差。分别计算3组随机实验的残差时间序列在去除CME前后的单站RMS值，RMS减少百分比见图 4。
|
图 4 自适应PCA时空滤波RMS减少率 Fig. 4 RMS reduction ratio of the adaptive PCA spatiotemporal filtering |
图 5为3组实验中滤波前后站点平均RMS值变化。结果表明，预处理后的GNSS残差时间序列经过传统的PCA时空滤波后的平均RMS值为4.8 mm，相较于滤波前平均RMS值减少约34.2%，而使用自适应PCA时空滤波后的3组平均RMS值分别约为4.4 mm、4.5 mm、4.4 mm，相较于滤波前平均RMS值分别减少约39.7%、38.4%、39.7%。综上，自适应PCA时空滤波方法滤波结果更优，且3组随机实验结果相近，表明初始站点的选取对滤波结果影响较小。
|
图 5 滤波前后平均RMS值 Fig. 5 Average RMS values before and after filtering |
GNSS坐标时间序列中通常存在与时间相关的噪声，选取幂律噪声(power-law noise, PL) 加白噪声(white noise, WN)的组合来描述GNSS坐标时间序列随机特性[10]。以整体PCA时空滤波与第1组自适应时空滤波结果为例，使用Hector软件的PL+WN模型对滤波前后的残差时间序列进行噪声估计[14]。图 6为滤波前后177个GNSS参考站残差时间序列的PL幅值分布。
|
图 6 滤波前后幂律噪声振幅 Fig. 6 Amplitude of the power-law noise before and after adaptive filtering |
滤波前后PL的变化见表 1。可以看出，相较于滤波前，PCA时空滤波方法使PL的幅值平均下降约40.0%；自适应PCA时空滤波方法使PL的幅值平均下降约57.8%。另外，使用PCA时空滤波方法可使WN的幅值平均下降约7.6%，使用自适应PCA时空滤波方法可使WN的幅值平均下降11.3%。由此可以看出，针对大尺度区域，相较于整体PCA时空滤波方法，自适应PCA时空滤波方法能够进一步减弱噪声的影响。
|
表 1 不同滤波方法对幂律噪声振幅的削弱程度 Tab. 1 Weakening degree of power-law noise amplitude by different filtering methods |
本文将大区域自适应地划分为若干子区域，并在每个子区域内应用PCA时空滤波，以更好地适应大尺度区域GNSS坐标序列非均匀共模误差分布和异常点影响。
陆态网垂向GNSS坐标序列的时空滤波分析实验表明，自适应方法可避免针对大尺度区域的人为分区对结果的影响，有更好的适用性；同时，自适应时空滤波可以剔除异常站点，削弱异常站点对时空滤波结果的影响。相对于整体PCA时空滤波方法，自适应时空滤波方法能更有效地减小GNSS残差时间序列的RMS值，并有效削弱坐标残差序列中的幂律噪声，提高数据质量与可靠性。
| [1] |
Wdowinski S, Bock Y, Zhang J, et al. Southern California Permanent GPS Geodetic Array: Spatial Filtering of Daily Positions for Estimating Coseismic and Postseismic Displacements Induced by the 1992 Landers Earthquake[J]. Journal of Geophysical Research: Solid Earth, 1997, 102(B8)
(  0) |
| [2] |
Dong D, 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) |
| [3] |
姜卫平, 王锴华, 邓连生, 等. 热膨胀效应对GNSS基准站垂向位移非线性变化的影响[J]. 测绘学报, 2015, 44(5): 473-480 (Jiang Weiping, Wang Kaihua, Deng Liansheng, et al. Impact on Nonlinear Vertical Variation of GNSS Reference Stations Caused by Thermal Expansion[J]. Acta Geodaetica et Cartographica Sinica, 2015, 44(5): 473-480)
( 0) |
| [4] |
Nikolaidis R. Observation of Geodetic and Seismic Deformation with the Global Positioning System[D]. San Diego: University of California, 2002
( 0) |
| [5] |
Shen Y Z, Li W W, Xu G C, et al. Spatiotemporal Filtering of Regional GNSS Network's Position Time Series with Missing Data Using Principle Component Analysis[J]. Journal of Geodesy, 2014, 88(1): 1-12 DOI:10.1007/s00190-013-0663-y
( 0) |
| [6] |
Yan J, Dong D N, Bürgmann R, et al. Separation of Sources of Seasonal Uplift in China Using Independent Component Analysis of GNSS Time Series[J]. Journal of Geophysical Research: Solid Earth, 2019, 124(11): 11 951-11 971 DOI:10.1029/2019JB018139
( 0) |
| [7] |
田云锋, 沈正康. GPS观测网络中共模分量的相关加权叠加滤波[J]. 地震学报, 2011, 33(2): 198-208 (Tian Yunfeng, Shen Zhengkang. Correlation Weighted Stacking Filtering of Common-Mode Component in GPS Observation Network[J]. Acta Seismologica Sinica, 2011, 33(2): 198-208)
( 0) |
| [8] |
Ming F, Yang Y X, Zeng A M, et al. Spatiotemporal Filtering for Regional GPS Network in China Using Independent Component Analysis[J]. Journal of Geodesy, 2017, 91(4): 419-440 DOI:10.1007/s00190-016-0973-y
( 0) |
| [9] |
He X X, Hua X H, Yu K G, et al. Accuracy Enhancement of GPS Time Series Using Principal Component Analysis and Block Spatial Filtering[J]. Advances in Space Research, 2015, 55(5): 1 316-1 327 DOI:10.1016/j.asr.2014.12.016
( 0) |
| [10] |
李昭, 姜卫平, 刘鸿飞, 等. 中国区域IGS基准站坐标时间序列噪声模型建立与分析[J]. 测绘学报, 2012, 41(4): 496-503 (Li Zhao, Jiang Weiping, Liu Hongfei, et al. Noise Model Establishment and Analysis of IGS Reference Station Coordinate Time Series Inside China[J]. Acta Geodaetica et Cartographica Sinica, 2012, 41(4): 496-503)
( 0) |
| [11] |
田云锋. GPS位置时间序列中的中长期误差研究[D]. 北京: 中国地震局地质研究所, 2011 (Tian Yunfeng. Study on Intermediate and Long-Term Errors in GPS Position Time Series[D]. Beijing: Institute of Geology, CEA, 2011)
( 0) |
| [12] |
贺小星. GPS坐标序列噪声模型估计方法研究[D]. 武汉: 武汉大学, 2016 (He Xiaoxing. Study on the Noise Model of GPS Coordinates Time Series[D]. Wuhan: Wuhan University, 2016)
( 0) |
| [13] |
Li W W, Shen Y Z. The Consideration of Formal Errors in Spatiotemporal Filtering Using Principal Component Analysis for Regional GNSS Position Time Series[J]. Remote Sensing, 2018, 10(4): 534 DOI:10.3390/rs10040534
( 0) |
| [14] |
Bos M S, Fernandes R M S, Williams S D P, et al. Fast Error Analysis of Continuous GNSS Observations with Missing Data[J]. Journal of Geodesy, 2013, 87(4): 351-36 DOI:10.1007/s00190-012-0605-0
( 0) |
2. School of Traffic and Transportation Engineering, Changsha University of Science and Technology, 960 Second Section of South-Wanjiali Road, Changsha 410114, China