基于PCA方法的GPS台站时间序列分析

基于PCA方法的GPS台站时间序列分析
张晶;王小亚;胡小工
【摘要】基于10 a以上的全球GPS台站数据,利用主成分分析法及其他数据处理方法,对台站时间序列进行预处理和结果分析,研究其中的非线性周期规律,探讨时间序列的主要影响机制.结果表明,主成分分析法可以将台站残差时空矩阵分解成若干正交成分,GPS台站时间序列的东西方向具有线性漂移趋势,全球大部分GPS台站都存在非线性周期规律,周年项和半周年周期占据主导地位.
【期刊名称】《大地测量与地球动力学》
【年(卷),期】2019(039)006
【总页数】7页(P613-619)
【关键词】时间序列;主成分分析;小波变换;GPS
【作者】张晶;王小亚;胡小工
【作者单位】中国科学院上海天文台,上海市南丹路80号,200030;中国科学院大学,北京市玉泉路19号甲,100049;中国科学院上海天文台,上海市南丹路80
号,200030;中国科学院上海天文台,上海市南丹路80号,200030
【正文语种】中文
【中图分类】P228
GPS台站时间序列相对于其他空间大地测量技术呈现出明显的非线性变化，这些特征有助于研究台站非线性运动的物理机制，进而改正各种误差模型[1]。

为此，
本文基于主成分分析法(PCA)，将全球GPS台站坐标残差时空矩阵分解成若干正
交成分，分析和评估时间序列中各种非线性变化特征和统计规律，提取各分量坐标时间序列主分量的噪声特性、周年特性和季节特性，为进一步分析台站非线性变化机制奠定基础。

1 基于主成分分析法的台站时间序列分析
1.1 台站时间序列特征
理想情况下，GPS台站时间序列的噪声特性应表现为纯白噪声，然而由于地球物
理效应及与GPS技术相关的系统误差，会产生潜在的有色噪声，如站点不稳定性
导致闪烁噪声和随机游走噪声[2]。

连续GPS台站测量可以检测地壳的水平运动和垂直运动[3]。

图1给出了GPS台站在水平和竖直方向的坐标速度场，图中显示了形式精度优于0.3 mm/a的台站，该速度场及其所反映出的板块运动整体情况与
其他学者结果一致[4]。

大气压负荷、非潮汐海洋负荷和水文负荷与GPS台站垂直位移强相关[5]，GPS数据处理模型的不完善及未模型化的GPS系统误差也可能导致虚假的非线性位移，包括电离层延迟高阶项、大气潮汐、海洋潮汐、对流层延迟处理方式等[6]。

因此，在台站坐标时间序列的模型化拟合过程中，应最大限度地
减少参数个数，并选择合适跨度的数据序列进行拟合，以提高拟合参数的精度[7]。

图1 形式精度优于0.3 mm/a的GPS台站水平和垂直速度场Fig.1 Horizontal and vertical velocity fields of GPS stations with a form accuracy better than
0.3 mm/a
1.2 主成分分析法机理
PCA方法的数学原理和具体步骤如下。

1)降维思想。

PCA采用降维的思想，即将原有的多个属性转化为少数综合属性，
已有属性通过适当的线性组合变换成相互独立的新属性，让数据本身体现出“空间响应均匀分布”特性。

设X是m×n矩阵，xi为其元素，将X视为n个m维点的集合(对坐标残差序列来说，n代表台站数，m代表时间序列)，寻找一个k维空间(k<n)，使得原始数据点与X内所有数据点在k维空间中的投影点最接近，从而可以用得到的这些投影点来描述m个原始数据点，保证信息损失最小。

当k=1时，挑选一条过原点的直线，使原始数据点与到该直线上的投影最接近，记该直线为L1，其单位向量为u1，则xi到L1的投影为所谓接近的准则是误差平方和达到最小。

由勾股定理，有：
(1)
由于固定，上式的最小化问题即是使最大化，利用矩阵表示则是求u使得uTXTXu最大。

由二次型极值的性质可知，最大值就是XTX的最大特征值λ1，达到这个最大值的u就是λ1所对应的特征向量u1，由此决定的L1就是对数据具有最优拟合的一维空间。

然后进行二维空间拟合，求一个与u1正交的单位向量u2，由u1和u2产生二维空间，作m个点到这个二维空间的投影，使得这些投影点与原始数据点最接近，满足的单位向量u2使得达到最小。

同样由勾股定理得：
(2)
上式达到最小即是达到最大，u2即为XTX第二大特征值所对应的特征向量。

依此类推，若设XTX的特征值λ1>λ2>…>λk>0，则其对应的标准正交特征向量
分别为u1，u2，…，uk。

可以得到结论：对任意k(0<k<n)，在所有可能的k维
子空间，以XTX的前k个标准正交特征向量u1，u2，…，uk张成的子空间，使
x1，x2，…，xn与其在子空间上的投影具有最小误差平方和。

2)主成分对原始数据的恢复。

记yi=Xuj, j=1,2,…,k,表示X的k个主成分，用矩阵
表示为：
Y=[y1,y2,…，yk]=XU
(3)
其中，U=(u1,u2,…，uk)。

由于
YTY=UTXTXU=UT(UΛUT)U=Λ
(4)
主成分之间相互正交，且第j个主成分的模为由于特征向量U的正交性有X=YUT，表示用主成分可恢复原始数据。

若选择部分主成分对原始数据逼近，可得到最优化近似为：
(5)
对任意p(1<p<k)，逼近误差平方和为：
dp≜
(6)
式中，tr(XTX)表示X的总变差，λj表示主成分yj对数据X的变差贡献，表示前p 个主成分对数据X的总变差贡献率。

从数据恢复损失信息最小来说，用主成分逼
近是最优的。

3)主成分计算。

采用经典正交分解的方法计算主成分。

首先要对数据矩阵X进行
中心化处理，将X转化成一个列均值为0的中心化矩阵。

计算X的协方差阵B的特征值λ1>λ2>…>λk>0及其对应的标准正交特征向量u1,u2,…,uk，B矩阵及其
谱分解表示为：
B=XTX=UΛUT
(7)
其中，U=(u1,u2，…，uk)，Λ=diag(λ1,λ2,…,λk)。

计算主成分对总变差的累积贡献率：
αp≜
(8)
与已选定的累积贡献率c进行对比，确定使αp>c最小的p值。

计算p个主成分：yj=Xuj，j=1,2,…,p
(9)
对GPS台站残差时间序列组成的时空矩阵X(m个历元n个站点组成的m×n阵，不失一般性，设m>n)进行经典正交分解或奇异值分解，获取协方差矩阵B的特
征向量阵U、特征值对角阵Λ，其中U和Λ均为n×n阵。

主成分A表示为：
(i=1,2,…,m;j=1,2,…，n)
(10)
对特征值对角阵Λ进行降序排列并重排主成分以及特征向量，通过探查每个主成
分特征向量的空间响应，选择特征向量具有空间统一响应的主成分来进行分析。

2 实验验证和精度评估
2.1 数据准备
根据He等[8]利用4种技术得到的多年SINEX周(日)解以及全球并置站信息，引
入方差分量估计加权方法，建立大型综合法方程系统，解算综合地球参考框架和地
球定向参数时间序列。

基于自主研发的综合软件STRF进行调优和自动化处理，生成J2005.0时刻的台站坐标时间序列。

PCA方法的核心是对数据进行正交分解，而正交分解要求时空数据矩阵连续无间断。

因此，本文选择数据缺失低于15%、数据跨度10 a以上的GPS台站进行坐
标时间序列分析，代表性GPS台站共299个，时间序列采样间隔为7 d。

选择的GPS台站中，北美和欧洲西部最为密集，赤道附近台站数量较少，南北半球分布
不均[9]。

2.2 台站时间序列特征提取
2.2.1 时间序列数据预处理
大部分GPS台站的时间序列都存在明显跳变和间断点[10]，会影响模型解算的精
度和可靠性。

在分析GPS台站坐标序列时，需要同时在时域和频域来研究其全貌
和局部性质，从而更加有效地探测GPS台站坐标时间序列所含的非线性周期信息。

因此，首先要对选择的GPS台站残差坐标时间序列进行数据预处理[11]。

广义离群检测算法GESD (generalized extreme studentized deviate)可以探测
跳变的时间点，检验服从近似正态分布的一个单变量数据集中的一个或多个粗差。

本文采用小波变换结合广义离群检测算法，将时间序列中含有跳变的GPS台站坐
标周解从原始坐标时间序列中剔除。

同时，时间序列趋势项通常称为台站运动速度，其中异常的水平速度可能与台站的稳定有关，可对比推测可能存在的漂移运动[12]。

本文对有明显趋势项的台站进行去趋势处理并进行对比分析。

以AZCN站为例，设其未发生跳变，经过综合后得到的台站残差序列如图2所示，可以看出，影响了2006～2010年的结果。

如图2(a)中，3个方向均存在显著的
线性项特征，为此，将该站在约化儒略日 55 425.0前后的站坐标进行分段处理，去除趋势项和跳变，解算后的时间序列如图2(b)。

对比发现，时间序列变得更平稳，进而周期性分析将更加可靠。

图2 AZCN台站残差序列图(去除跳变和趋势项前后的分段线性项)Fig.2 Residual sequence diagram of the AZCN station (segment linear term before and after the hop and trend terms are removed)
图3为CHPI站去除趋势项前后的时间序列，可以看出，其在东西方向表现出明显的线性趋势项，推测可能存在水平漂移运动；垂向坐标时间序列也存在较小的整体趋势，这可能与该地区的地下水变化、台站稳定性和其他未知效应有关。

扣除线性趋势项后(图3(b))，整体变化和周期项信号更加突出。

BJFS站是研究中国大陆及其邻近地区水平位移场的核心站之一(图4)。

可以看出，BJFS站水平方向的波动比垂直方向大，说明该台站垂直方向比水平方向稳定。

其中南北向振幅约为-10～10 mm，东西向约为-20～40 mm；数据预处理之后，南北向振幅约为-5～5 mm，东西向约为-20～20 mm，垂直方向变化不大。

2015-12南北、东西向存在微小跳变，参考相关资料可以推测，跳变是由于接收机和天线更换及数据处理中GPS卫星截止高度角改变所致。

为完成PCA对时空数据矩阵的正交分解，需要对时间序列中缺失的数据进行最小二乘内插补齐，获取3组完整的初始时空数据矩阵。

通过对所有10 a以上较连续GPS台站时间序列的分析和探测，消去其中存在的跳变、线性趋势项和野值，弥补缺失历元数据，为周期项特征探测和分析准备好“干净”的数据。

图3 CHPI站去除趋势项前后的时间序列Fig.3 CHPI station time series before and after the trend item is removed
图4 BJFS站去除跳变前后的时间序列Fig.4 BJFS station time series before and after the transition is removed
2.2.2 用主成分分析方法分析时间序列
对预处理后的时空矩阵进行主成分分析，每组数据分别计算299个主成分、空间特征向量及主成分特征值。

由于PCA方法的特性，在GPS台站坐标时间序列中提
取主成分分量是独立的，因此可以对每个方向的分量进行特征提取。

3个方向的贡献率及主分量统计见表1，主分量对应的空间向量表示的是时空数据的空间相应特征——东西、南北和垂向的坐标时空矩阵第1主成分贡献率分别达
到10%、12%和9%，第2主
表1 3个方向主要分量的累计贡献率Tab.1 Cumulative contribution rate of main components in three directions主分量ENU10.096 4080.117 4260.085 73120.155 9750.182 9740.152 30930.197 0330.232 2590.207 09940.234 4010.273 1380.249 14850.269 5120.308 9920.289 539
成分的贡献率达到6%、7%和7%，其他主成分的贡献率更低。

由此可见，越靠前的主成分分量越占据主导地位。

分别对各主成分分量进行傅里叶分析，提取周期项，得到图5～7。

表2统计了台站3个方向主成分的主要周期项。

表2 GPS台站周期信号及幅度Tab.2 Cycle signal and amplitude of GPS stations分量N/mmE/mmU/mm11 a(2.0)、0.5 a(1.5)1 a(41.7)1 a(9.8)、0.5
a(2.9)、0.25 a(2.8)21 a(8.0)1 a(25.9)1 a(8.1)31 a(7.7)1 a(11.2)1 a(4.3)461.4周(2.0)、1 a(1.9)1 a(6.2)、0.5 a(4.7)1 a(3.1)51 a(5.6)1 a(11)1 a(3.9)
图5 台站时间序列PCA分析后E方向的主要分量Fig.5 Main component of E direction after PCA analysis of station’s time s eries
图6 台站时间序列PCA分析后N方向的主要分量Fig.6 Main components of the N direction after PCA analysis of station’s time series
可以看出，大部分台站在3个方向上的主要周期项均为周年项，同时垂直方向的
第1主分量存在半年和季节性变化周期；另外南北方向的第1主分量、东西方向
的第4分量也存在半年周期项。

幅度方面，东西方向的幅值远大于另外2个方向，前3个主分量的周年项振幅依次为41.7 mm、25.9 mm和11.2 mm，说明GPS
台站主要表现出东西方向的位移。

根据幅度值可以看出，3个方向各主分量的幅值
基本符合从大到小的排列规律，与贡献率值基本吻合。

GPS台站时间序列的垂直方向表现出明显的季节性。

由于目前的地球参考框架模
型仅考虑台站的线性运动情况，因此台站的周期信号均被残差序列所吸收。

造成台站周期信号的可能原因较多，如地表温度变化、大气潮汐效应、陆地水和海洋环流等地理现象或者观测技术系统误差。

从图6看出，东西方向的周期性运动振幅小
于其长期项运动振幅，水平方向的主要分量中1 a周期项最为突出。

从上面的分析可以看出，目前线性STRF地球参考框架没有考虑非线性的周期运动和震后形变，这些信息就包含在台站坐标的残差序列中。

GPS台站坐标残差序列
有明显的周年项和季节项，这些周期项可能是由地球物理现象如地球表面物质随时间变化产生的负荷效应(如海潮、大气潮汐效应、积雪、土壤水、海洋非潮汐影响、地表温度变化、陆地水和海洋环流等)引起的，也可能是因潮汐频率和GPS卫星轨道频率之间差异引起的长周期影响导致的海潮改正误差放大或GPS观测技术系统
误差。

从表2看出，东西向的周期性运动振幅小于其长期项运动振幅，水平方向
的主要分量中1 a周期项最为突出。

图7 台站时间序列PCA分析后U方向的主要分量Fig.7 Main components of the U direction after PCA analysis of station’s time series
3 结语
本文利用全球299个具有10 a以上、观测较为连续的GPS台站数据，通过主成
分分析法进行台站时间序列变化分析。

结果表明，主成分分析法可以将时空矩阵分解成若干正交成分，有效提取时间序列中的周年项、半年项等周期项。

在GPS台
站东西向、南北向和垂向3个方向的周期性特征中，1 a周期具有一定的普遍性。

另外大部分台站受板块运动影响，在东西向有较大的线性漂移速度。

台站坐标残差序列信息研究及其结果分析对沉降或变形监测及趋势预报、地球物理相关方向研究和地球参考框架精度及GPS数据处理模型误差精度提高等具有重要的参考意义，
为下一步寻找其机制、探讨其与地球物理现象的相关性奠定了基础。

参考文献
【相关文献】
[1] 程宗颐, 张飞鹏, 董大南,等. 利用GPS监测中国地壳的垂向季节性变化[C].中国地球物理学会年会, 北海，2002(Cheng Zongyi, Zhang Feipeng, D ong Da’nan, et al. Monitoring the Vertical Seasonal Variation of the Chinese Crust Using GPS[C].Annual Meeting of the Chinese Geophysical Society,Beihai,2002)
[2] 王小亚. GPS在地球物理方面的应用[D]. 上海：中国科学院上海天文台, 2002(Wang Xiaoya. Application of GPS in Geophysics[D]. Shanghai: Shanghai Astronomical Observatory,CAS, 2002)
[3] 田亮. GPS测站坐标非线性变化规律分析与机制研究[D]. 郑州：信息工程大学, 2011(Tian Liang. Analysis and Mechanism Research of Nonlinear Variation Law of GPS Station Coordinates[D]. Zhengzhou: Information Engineering University, 2011)
[4] 吴伟伟. 华北地区GPS连续站坐标序列特征研究[D]. 北京：中国地震局地震预测研究所, 2014(Wu Weiwei. Study on the Coordinates of GPS Continuous Stations in North
China[D].Beijing: Institute of Earthquake Forecasting, CEA, 2014)
[5] 田亮, 刘武凤, 张建东,等. 基于GPS坐标残差序列的全球测站非线性变化规律统计[J]. 地理空间信息, 2013(4):70-71(Tian Liang, Liu Wufeng, Zhang Jiandong, et al. Statistics of Nonlinear Variation of Global Stations Based on GPS Coordinate Residual Sequences[J]. Geospatial Information, 2013(4): 70-71)
[6] 符养. 中国大陆现今地壳形变与GPS坐标时间序列分析[D]. 上海：中国科学院上海天文台, 2002(Fu Yang. Current Crustal Deformation and GPS Coordinate Time Series Analysis in Mainland China[D]. Shanghai: Shanghai Astronomical Observatory,CAS, 2002)
[7] Dong D, Dickey J O, Chao Y, et al. Geocenter Variations Caused by Atmosphere, Ocean and Surface Ground Water[J]. Geophysical Research Letters, 1997, 24(15):1 867-1 870 [8] He B, Wang X Y. Combination of Terrestrial Reference Frames Based on Space Geodetic Techniques in SHAO: Methodology and Main Issues[J]. Research in Astronomy and Astrophysics, 2017, 17(9):1-14
[9] 叶叔华, 黄珹. 天文地球动力学[M].济南：山东科学技术出版社, 2000(Ye Shuhua, Huang Cheng. Astronomical Geodynamics[M]. Ji’nan:Shandong Science and Technology Press, 2000)
[10] Altamimi Z, Métivier L, Collilieux X. ITRF2008 Plate Motion Model[J]. Journal of Geophysical Research: Solid Earth, 2012, 117(B7)
[11] Altamimi Z, Boucher C, Willis P. Terrestrial Reference Frame Requirements within GGOS Perspective[J]. Journal of Geodynamics, 2005, 40(4):363-374
[12] 张鹏, 蒋志浩, 秘金钟,等. 我国GPS跟踪站数据处理与时间序列特征分析[J]. 武汉大学学报:信息科学版, 2007, 32(3):251-254(Zhang Peng, Jiang Zhihao, Bei Jinzhong, et al. Data Processing and Time Series Feature Analysis of GPS Tracking Stations in China[J]. Geomatics and Information Science of Wuhan University, 2007, 32(3):251-254)。