区域GNSS基准网解算及其时间序列分析

尹 彤，鲁欣宇，聂建亮，蒲永峰，武慧琳

（1.自然资源部测绘标准化研究所，西安 710054；2.河南省遥感测绘院，郑州 450003；3.陕西测绘地理信息局，西安 710054；4.自然资源部第二地形测量队，西安 710054）

0 引言

随着高精度全球卫星导航系统（global navigation satellite system，GNSS）网络的不断建设和完善，至今，在GNSS基准站上已经积累了二十多年的观测数据，这些数据成果将形成长期GNSS坐标时间序列，其中蕴含着很多地表变化及地球框架变化的周期和非周期的运动信息。如何实现区域 GNSS网的联合解算，从而获取长期连续的坐标时间序列，是进行基准站稳定性分析的前提。在连续的GNSS台站的时间序列中，包含地壳、地表运动信息，主要分为构造信息和非构造信息。对非构造运动信息需要进行分类剔除。如海潮运动、天线相位中心、地表大气负荷、地表水迁移等，才能精确获取纯净的时间序列，从而分析基准站的稳定性。

目前国内 GNSS网解算都停留在综合平差阶段，没有考虑长期时间序列噪声和误差分析过程，故成果精度不高。本文探讨使用区域GNSS数据，采用加米特（GAMIT）/格洛布克（GLOBK）软件进行区域联合解算，获取基准站的坐标时间序列，然后利用主成分分析的方法提取 GNSS时序中的共模噪声，并在时序分析中将其剔除，继而基于纯净的 GNSS时间序列，进一步提高基准坐标和运行速度的解算精度。

1 基准网解算策略及时间序列分析方法

1.1 区域GNSS基准网解算策略

采用GAMIT/GLOBK软件进行区域GNSS基准网解算。GAMIT/GLOBK软件是一款精密的GNSS定位和定轨软件。它由美国麻省理工学院（Massachusetts Institute of Technology, MIT）和加州大学圣地亚哥分校的美国斯克里普斯海洋研究所（Scripps Institution of Oceanography, SIO）共同开发和维护。

GAMIT处理双差观测量，利用双差直接消除卫星及测站的时间误差，采用最小二乘原理进行参数估计进行平差计算。软件适合进行上千千米的超长基线解算，点位精度可以达到毫米级别，是全球公认的高精度 GNSS长基线处理软件。GAMIT软件可以实现超长基线的GNSS多星座数据混合解算，它由许多不同功能的模块组成, 这些模块可以独立地运行。

GLOBK是一个卡尔曼滤波器，其主要目的是综合空间大地测量和经典大地测量的初步处理结果完成数据的后处理。它的输入数据一般是一些准观测量，如测站坐标、地球自转参数、卫星轨道及它们的方差和协方差。GLOBK可以融合多种基线成果格式，实现多源基线成果的联合平差。GLOBK软件的主要功能有：1）实现多源基线成果和多天综合基线成果的联合平差；2）获取 GNSS跟踪站的坐标分量时间序列、计算站点运动速度场。GLOBK平差计算主要分为：1）输入文件的格式转换；2）平差和输出命令文件准备；3）平差命令globk或glred、成果输出提取。通过GLOBK软件可以获取多时段的GNSS站点的坐标时间序列。

1.2 GNSS时间序列分析方法

GNSS时间序列分析目的就是进一步改正和过滤 GNSS数据解算过程中未考虑的噪声，获取更精确的基站位置和运动状态（位置和速度）。共模误差（common mode error，CME）是GNSS 坐标时间序列呈现周期性的主要原因之一，需对CME进行改正，从而减少台站坐标时间序列坐标分量的上下波动，提高时间序列的信噪比，减弱坐标序列周期信号。主成分分析方法（principal components analysis，PCA）的应用，能够识别出规模较大的区域网络的 CME，作为第一主成分或主要成分的线性组合，并对单个站点进行“区域调整”。主成分方法可以提高CME提取的可靠性，并顾及形式坐标误差。

主成分分析是一种通过线性变换来简化数据集的方法，它将数据变换到一个新的坐标系统中，使得数据投影的第一大方差在第一个坐标（第一主成分）上，第二大方差在第二个坐标（第二主成分）上，依次类推。主成分分析经常用减少数据集的维数，同时保持数据集的对方差贡献最大的特征。旨在利用降维思想，把多指标转化为少数几个综合指标。在GNSS基准站的长期时间序列中，运用PCA主成分分析可以把随时间变化的序列分解成时间域的主分量和空间域的特征向量。主成分分析模型由 GNSS坐标向量残差时间序列构造观测矩阵R(t,r) （= 1, 2, ···,；= 1, 2, ···,）：t表示时间域中个时间元素的第个；r表示空间域里分布的个研究对象中的第个，其中坐标序列需首先经过内插实现连续化（用个站组成的 GNSS站网构造观测矩阵）,同时分别为每个地心分量（北、东和高）创建分量矩阵。基于观测矩阵的元素协方差的每个元素的对称矩阵B(t,r)定义为

主成分分析方法通过对协方差矩阵(,r)进行特征值分解，分解为和（为正交矩阵，为具有个非零对角特征值{}(≥)的对角矩阵），其计算方法为

那么，采用特征值和特征向量矩阵中对应列的降阶来定义各阶主成分为

式中：(t)为(t,r)的次主成分；v为第个特征值对应的特征向量。

PCA是去除CME的有效算法，所得到的正交分量（主分量）按变化幅度从大到小排序。PCA方法的基础是选择主分量，主分量超过一半的对应归一化响应的值大于25%。CME的计算方法为

式中：t为时间域的个时间元素；r为空间域分布的第个研究对象；(t)为式（3）中(,r)的次主成分；v为第个特征值对应的特征向量；为前个主成分。

这种方法通过正交变换重建滤波残差，从而获取 GNSS残差序列的至少 6个主成分，其比重依次降低。进而获取 GNSS坐标序列残差中的共模误差。PCA方法通过减去占数据方差中的主要分量来提取共模误差，然后利用最大似然估计对其进行分析。这样不仅数据的均方根值会减小，随机部分也更接近白噪声。

2 区域GNSS网综合平差及其时间序列分析

2.1 数据收集及整理

本文采用中国境内及周边的国际 GNSS服务（International GNSS Service, IGS）站的原始观测数据作为解算数据，收集了14个IGS站从2013—2019年至少 6 a的数据。基线解算软件采用GAMIT/GLOBK10.6，时间序列分析软件采用准观测组合分析（quasi-observation combination analysis,QOCA）软件的姆洛德（Mload） 模块和 PCA模块。采用 IGS站的名称及时间跨度如表1所示，点位分布如图1所示。

表1 IGS站点名称及时间跨度表

图1 14个IGS站点位分布图

本文收集的其他改正模型资料包括全球大气模型（Atmospheric-M）、全球积雪模型（Snow-M）、全球土壤湿度模型（Soil-M）、非潮汐改正模型（Ecco-M）等的数据，时间跨度为2013—2019年，如表2所示。

表2 改正模型名称及时间跨度度统计表

从表2的时间可以判断，时间序列中如果需要加入 Mload的各种改正的话，必须选用 2 013.0~2 019.0 a的这个时段。所以本次数据处理在最终时序分析中使用的数据时段为2 013.0~2 019.0 a。

2.2 基线解算及综合平差

本文基线解算采用6 a的连续数据，考虑到半月周期，数据间隔采用7天一个时段，共计312个时段。采用GMAIT软件对14个IGS站进行长基线解算。

对312个时段的基线解算进行精度统计。单天解的标准化均方根（normalized root mean square,NRMS）代表了基线解算的精度情况，通常情况下，NRMS小于0.3。除了11个时段的NRMS值超限之外，其他301个时段的NRMS值统计表3所示。

表3 单天解NRMS统计表

2.3 平差获取坐标时间序列

本文将采用 GLOBK软件对所有符合要求的单时段基线解进行整网约束平差。平差选择 6个IGS站作为固定站，这6个分别为YSSK、KHAJ、PIMO、KIT3、SUWM、PIMO，对固定站的平差约束量采用其坐标和速度的一倍中误差，平差基准采用国际地球参考框架（The International Terrestrial Reference Frame, ITRF）2014，参考历元为2010.0。

通过格尔雷德（GLRED）解算整体综合解，获得站点在各个时段的坐标值，并形成坐标分量的长期坐标时间序列（N、E、U）。此处仅呈现出我国境内的LHAZ（拉萨站）站的单天剔除粗差之后的时间序列，如图2所示。

从图2可以看出，在LHAZ站上，、及三个方向上均有显著的周期波动，尤其在方向上波动较大：方向波动振幅约6 mm、方向振幅约7 mm、方向振幅约20 mm。

图2 LHAZ站网解时间序列（剔除粗差之后）

2.4 主成分分析（PCA）

采用主成分分析（PCA）模型时：1）去除时间序列中的阶跃项、并对序列的趋势项进行改正，改正后的序列只剩下年周期项和半年周期项；2）对序列残差进行去除周期项分析，序列周期信息包括了由于地壳质量负荷引起的周期项，地壳质量负荷改正主要包括了如：大气改正、积雪改正、土壤湿度改正、非潮汐改正，本文采用全球大气资料Atmospheric-M、全球积雪资料Snow-M、全球土壤湿度资料Soil-M、非潮汐改正资料Ecco-M等模型，根据坐标序列中的坐标值进行差值计算，获取所有站点的坐标时间序列值对应的质量负荷改正量，并实现对时间序列的改正，改正后获得的所有站点坐标序列将作为一个整体序列群输到主成分分析模型中去。

对改正后的相对纯净的时间序列数据进行PCA主成分分析。从PCA分析来看，CME的存在是造成残差序列中存在明显的周期变化的主要原因，尤其在竖直方向，周期性表现尤为明显。为了对 CME进行分离，在滤波过程中，首先由 PCA法分别计算得到14个IGS站集群、在同一时间域的前12个主成分对3个坐标分量（、、）的空间响应，其中：主成分 1对坐标分量（、、）的贡献率分别为 70.83%、69.32%、68.23%，主成分2对坐标分量（、、）的贡献率分别为12.89%、14.43%、10.03%，主成分 3对坐标分量（、、）的贡献率为 5.43%、9.56%、3.32%。

14个IGS站集群在同一时间域的前3个主成分，累计贡献率累计为 89.166%、85.453%、83.95%，所有第一、第二、第三主成分的累计贡献率综合了大部分的空间响应（>87%）。

以上统计表明，采用PCA方法可得到前3个主成分作为该数据群的CME，然后在各个站点的坐标序列中，需要针对CME进行滤波处理。这样滤波之后得到纯净的GNSS站点坐标时间序列。最后，基于纯净的时间序列计算精确的站点坐标和速度。

2.5 PCA滤波结果分析

通过地表负荷改正和 PCA滤波处理剔除了CME之后，各个站点时间序列曲线的周年运动的振幅有所减小：平均年振幅变化减少最小为3.22%；最大减少了65.23%；平均减少了32.52%。

图3为LHAZ站的结果。从图3可以看出，曲线振幅减小了，收敛程度更高了。

图3 LHAZ站PCA滤波前后时序

在 ITRF2014框架下计算站点的框架坐标和速度场，并将结果与IGS网站公布的ITRF2014框架的坐标和速度进行比较，差异情况如图4、图5所示。

图4 滤波后坐标值与IGS成果对比

图5 滤波后速度值与IGS结果对比

从图4、图5可以看出，14个 IGS站点经过滤波处理后的解算的坐标与 IGS成果差异小于9 mm，其速度值与IGS成果差异小于3 mm/a。

3 结束语

本文利用 GAMIT/GLOBK软件实现了区域GNSS网络的综合解算和平差，并采用 PCA方法对长期坐标时间序列进行了稳定性分析。重点介绍了如何采用滤波的方法来探测和剔除共模误差，从而获得纯净的坐标序列，提高了坐标精度。得到结论如下：

1）在全球GNSS网络平差中，需对约束站给予合适的坐标和速度约束才能保证平差结果的精确。

2）共模误差是GNSS坐标时间序列程序周期性信号的主要因素，采用主成分滤波可以有效去除共模误差的影响。

3）实施负荷改正和主成分分析等区域滤波后，基站坐标时间序列趋于平稳，波动振幅有所减小，大大提高了时间序列的信噪比。

4）时间序列分析目的就是分类各种信息，并加以改正、剔除、过滤，从而获取纯净的时间序列。基于纯净的时间序列数据平差得到更高精度的站点位坐标和速度场。