全球GPS测站垂向周年变化统计改正模型的建立

傅彦博，孙付平，朱新慧，刘 婧

信息工程大学，河南 郑州 450001

随着各类空间大地测量技术的迅猛发展，毫米级地球参考框架正逐步地建立并应用于大地测量领域中。温度变化被发现为GPS测站垂向周年变化作出了一定比例的毫米量级贡献[1-3]；季节性温度变化引发的测站水泥墩与基岩的热胀冷缩被发现与测站周年变化存在很强的正相关性，会造成测站高程方向的季节性变化[2-8]；大气负荷、水文负载、非潮汐海洋负载等负载变化是引起测站垂向周年变化的主要因素[9-17]；高阶电离层延迟是造成中纬度测站垂向周年，东西方向半周年以及低纬测站垂向半周年运动的主要原因之一[18-20]；文献[21]统计了全球492个GPS测站坐标时间序列的周期项，发现周年项在全球具有普遍性；文献[22]发现了南北半球周年项初始相位存在约180°的系统差。由于GPS测站的非线性变化成因复杂，目前国际上还未建立包含多种机制影响的非线性变化理论改正模型。若基于现有全球GPS实测数据，分析测站垂向坐标时间序列中周期项的全球分布规律，并根据此规律构建全球测站垂向非线性变化统计改正模型，将会从整体上削弱全球测站垂向非线性运动，并为从全球角度研究测站垂向非线性变化提供参考。

本文在全球范围内选取了461个连续观测时段在2年以上的GPS测站(其中北半球357个，南半球104个)，提取了测站垂向周年项的振幅与初始相位，分析得出了周年项振幅与初始相位的全球分布规律函数，基于此函数构建了全球GPS测站垂向周年变化统计改正模型，并对模型的改正效果进行了测试与评估。

1 数据源与提取周期项

1.1 数据源

本文采用的测站坐标残差数据可通过ITRF2008官方网站(http:∥itrf.ensg.ign.fr/ITRF_solutions/2008/)下载。为了便于后续数据分析，本文对原始数据进行了粗差剔除与零均值化处理。

1.2 周期项提取

根据谐波分析的基本思想，含有若干周期项的时间序列可以看成由一系列不同频率的周期波叠加而成，因此只需要确定某一周期波的准确频率以及参考历元，就能计算出其对应的振幅和初始相位[3]。本文采用周期图法对测站垂向坐标残差时间序列进行了频谱分析，探测并统计了周期项。鉴于篇幅，只列出了BJFS站垂向坐标残差时间序列的频谱图，如图1所示。

图1 BJFS站垂直方向坐标时间序列的频率谱Fig.1 Frequency spectrum of vertical coordinates time series of BJFS station

图1中能量最高频点的频率为0.002 79 Hz，可以得出所对应的周期为358.4 d，接近于一周年。确定了周期项所对应的频率后，可以利用最小二乘拟合的方法得出周期项的振幅和初始相位[23]。

1.3 全球GPS测站垂向时间序列周期项统计

本文对全球461个GPS测站垂向时间序列所包含的各类周期项进行了统计，如表1所示。

表1 全球GPS测站各类周期项统计

在实际探测中，发现测站时间序列中很少存在严格的周年项，主要为类周年项，因此本文在提取周年振幅与初始相位时，将周期介于350～380 d 之间的周期项都视为周年项[21-24]。

从表1可以看出，周年项在全球具有普遍性，因此研究周年项振幅与初始相位的全球分布规律，是建立全球测站垂向非线性变化统计改正模型的重要手段。

2 全球测站周年项分布规律

建立基于周年项的统计改正模型，需要分析周年项振幅与初始相位的全球分布规律。考虑到南北半球存在以下差异：①南北半球存在非对称构造变化[25]。②南北半球的海洋与陆地面积比例不同，导致南北半球的气候变化、气象条件、GPS测站部署有所差异[26]。因此，分别对南北半球测站垂向周年项的全球分布规律进行研究。

2.1 周年项振幅全球分布规律

图2和图3分别展示了北半球、南半球测站垂向周年项振幅随纬度的分布规律。从中可看出，南北半球垂向周年项振幅随纬度的总体变化趋势都呈现了单一的凹凸性，且总体趋势较为平缓。由于二次项函数的凹凸性单一，不存在拐点，且函数值随着横坐标的变化不会出现剧烈的震荡，结合振幅随纬度的分布特征，对垂向周年项振幅全球分布规律进行建模时，用二次项拟合的方法最合适；考虑到函数在赤道处的连续性，使用分段拟合的方法，最终得出了南北半球振幅的分布规律函数。北半球垂向周年项振幅分布规律函数为

AN(B)=-0.002 417B2+0.206B+3.2

(1)

式中，B为纬度，下同。

南半球垂向周年项振幅分布规律函数为

AS(B)=-0.001 378B2+0.091B+3.2

(2)

对南北半球的振幅分布规律函数进行了t检验(显著性水平为0.05，零假设为拟合函数值与振幅值来自同一分布)，检验结果如表2所示。

表2 南北半球振幅分布规律函数显著性检验结果

经过显著性检验，在显著性水平为0.05的条件下，南北半球振幅拟合函数的H值都为0，则接受零假设，即模型拟合值与真实振幅值为同一分布，说明二次项拟合函数能准确刻画振幅的分布特征。

从图2和图3中可以看出两个明显的现象：①随着纬度的增加，南北半球测站垂向周年项振幅呈现先增大后减小的规律；②北半球的垂向周年项振幅普遍高于南半球。

针对第1个现象，由于环境负载与测站垂向坐标季节性温度变化存在强相关性[15]，对于季节性变化剧烈的中纬度地区的GPS测站，其环境负载导致的垂向位移普遍高于其他地区；同时，环境负载是引起测站垂向周年变化的主要因素[9]。因此，中纬度地区的周年振幅高于低、高纬度测站。

针对第2个现象，由于地表环境负载是引起测站垂向周年变化的主要因素，同时北半球测站环境负载的垂向绝对平均坐标值普遍高于南半球测站[9]，因此导致北半球测站的周年项振幅普遍高于南半球测站。

2.2 周年项初始相位全球分布规律研究

提取周年项的初始相位时，需要统一各测站垂向时间序列的参考历元[3]。在兼顾时间序列长度与测站数目的基础上，本文选用2002年1月1日为参考历元。由于部分测站时间序列的观测时段不包含参考历元，因此参与初始相位全球分布规律分析的北半球、南半球测站分别为298、89个，测站数目总和占全球测站数目的84%。

图4和图5分别展示了北半球、南半球测站垂向周年项初始相位随纬度的分布规律。

根据南北半球垂向周年项初始相位的分布特征，同2.1节，使用分段二项式拟合的方法得出了垂向周年项初始相位的全球分布规律函数。北半球垂向周年项初始相位分布规律函数为

φN(B)=-0.000 624B2+0.051 2B-0.25

(3)

南半球垂向周年项初始相位分布规律函数为

φS(B)=0.000 615 8B2-0.045 57B-0.25

(4)

对南北半球的初始相位分布规律函数进行了t检验(显著性水平为0.05，零假设为拟合函数值与振幅值来自同一分布)，检验结果如表3所示。

表3 南北半球初始相位分布规律函数显著性检验结果

根据表3可知，模型拟合值与真实初始相位值为同一分布，说明二次项拟合函数能准确刻画初始相位的分布特征。

图5、图6之间最大的差异在于拟合曲线的凹凸。在中纬度地区，南北半球测站垂向周年项初始相位之间存在约180°的系统差，即垂向周年位移趋势相反。由于南北半球测站的季节性温度变化相反，测站环境负载与测站季节性温度变化存在强相关性[15]，且环境负载是导致测站垂向周年变化的主要因素[9]，因此在同一时间段内，季节性温度变化最剧烈的中纬度地区会呈现南北半球测站垂向周年位移趋势相反的现象，低、高纬度地区测站的季节性变化相对较弱，因此其垂向周年位移不会呈现明显的相反趋势。

3 建模及改正效果分析

3.1 建立统计改正模型

在只考虑周年项的条件下，统计改正模型时间序列M(t)为

M(t)=A(B)·cos[ω1at+φ(B)]

(5)

GPS测站时间序列采样周期为7 d/次，可以得出周年角频率ω1a

(6)

根据全球GPS测站垂向周年项振幅与初始相位的全球分布规律函数，建立了全球GPS测站垂向周年变化统计改正模型。北半球的GPS测站垂向周年变化统计改正模型为

(7)

式中，B为测站纬度；t0为测站时间序列的起始历元；t为拟合值所对应的历元，下同。

南半球的GPS测站垂向周年变化统计改正模型为

(8)

3.2 模型改正效果分析

利用建立的统计改正模型对全球461个测站的垂向时间序列进行了改正，并且根据时间序列的RMS对改正效果进行了评价。

设原始时间序列为O(t)，模型序列为M(t)，改正后的剩余残差序列为R(t)，原始时间序列所包含历元数为n。则原始时间序列O(t)的RMS为

(9)

改正模型序列M(t)的RMS为

(10)

剩余残差序列R(t)的RMS为

(11)

本文利用了模型改正比例km及剩余残差比例kr，评估模型改正效果

(12)

本文在全球低、中、高纬度各选取了具代表性的4个测站，分别统计了模型改正比例及剩余残差比例，如表4所示。

表4 所选各测站信息及改正效果统计

图6是低纬度4个测站的改正效果(其中蓝色线表示原始时间序列，黑色线表示模型改正序列， 红色线表示剩余残差序列，下同)。图7是中

纬度4个测站的改正效果。图8是高纬度4个测站的改正效果。

图2 北半球垂向周年项振幅随纬度分布规律Fig.2 Distribution law of vertical annual terms’ amplitudes of northern hemisphere with latitude

图3 南半球垂向周年项振幅随纬度分布规律Fig.3 Distribution law of vertical annual terms’ amplitudes of southern hemisphere with latitude

图4 北半球垂向周年初始相位随纬度分布规律Fig.4 Distribution law of vertical annual terms’ initial phases of northern hemisphere with latitude

图5 南半球垂向周年项初始相位随纬度分布规律Fig.5 Distribution law of vertical annual terms’ initial phases of southern hemisphere with latitude

图6 低纬度测站改正效果Fig.6 Correction effect of stations located in low latitudes

图7 中纬度测站改正效果Fig.7 Correction effect of stations located in middle latitudes

图8 高纬度测站改正效果Fig.8 Correction effect of stations located in high latitudes

结合表4与图6—图8，可以看出中纬度测站的平均模型改正比例最大，平均剩余残差比例最小，总体改正效果最好；低、高纬度测站的改正效果均不如中纬度测站。可能有以下两点原因：①中纬度地区由于季节性温度效应较为明显[6-7]，环境负载、热膨胀效应等因素与季节性温度变化存在强相关性[7-9]，且环境负载是导致测站垂向周年运动的主要因素[9]，因此垂向周年变化是中纬度地区测站垂向位移的主要周期性变化；低、高纬度测站由于受季节性温度效应影响较小，其他周期性变化在测站垂向位移中的占比高，导致中纬度测站垂向周年振幅随纬度的分布比低、高纬度测站更为集中，改正效果更好。②与原因①同理，中纬度测站垂向周年初始相位值分布比低、高纬度测站分布更为集中，因此与低、高纬度测站相比，中纬度测站的改正效果更好。

从图6—图8可看出，测站原始时间序列的主周期变化与模型序列中的周年变化并不完全契合，因此原始序列中非线性变化的不是严格的周年变化，而是类周年变化，印证了文献[27]的结论。

本文统计了全球461个GPS测站的改正效果，如图9和图10所示。

图9 全球GPS测站模型改正比例统计Fig.9 Statistics of model correction ratio of global GPS stations

图10 全球GPS测站剩余残差比例统计Fig.10 Statistics of residuals ratio of global GPS stations

根据统计结果，在所有测站中，剩余残差比例达50%～70%的测站占多数，可知统计改正模型能削弱多数测站30%～50%的垂向坐标残差。在模型改正比例达30%～40%的测站中，存在一部分模型改正比例达到130%～140%的测站，是因为此类测站的垂向周年项振幅小于模型拟合值。综合分析统计结果，改正模型能大幅度削弱GPS测站的垂向坐标残差并提升大规模测站数据处理效率。

3.3 统计改正模型评价

3.3.1 模型优点

(1) 本文建立的统计改正模型能削弱全球多数测站30%～50%的垂向坐标残差，精确分离测站的线性与非线性运动。

(2) 模型能统一对全球测站垂向坐标进行改正，提高了数据处理效率。

(3) 模型能直观地显示测站垂向周年项振幅以及初始相位的全球分布规律，为研究全球范围内周年性地球物理运动提供参考。

3.3.2 模型缺点

(1) 由于全球测站所处地区的地球物理条件各异，引起周年项的因素复杂，部分测站的周年项振幅与模型拟合值相差较大，导致其改正效果不理想。

(2) 二次曲线受局部异常值影响较大，异常值的存在会影响二项式拟合效果。

(3) 由于数据观测时段限制，有16%的测站不包含用于提取初始相位的参考历元，导致建模时数据利用率降低，影响此类测站的改正效果。

3.3.3 模型改进展望

(1) 半周年项也较为普遍地存在于全球GPS测站垂向坐标残差时间序列中。总结半周年项的全球分布规律，优化现有统计改正模型，对一些含有半周年项的测站将有更好的改正效果。

(2) 研究测站水平方向非线性变化项的全球分布规律，建立GPS测站三维非线性变化统计改正模型。

(3) 由于二项式拟合对异常值较为敏感，今后改进模型时，可以考虑使用抗差性更好的方法进行更精确的拟合。

(4) 对全球GPS测站数据使用最新的精密数据处理策略进行重新处理，使所有测站的坐标时间序列中包含相同的参考历元，得出更精确的周年项初始相位分布规律函数，精化全球GPS测站垂向周年变化统计改正模型。

4 结 论

本文基于实测数据，研究发现了全球GPS测站垂向坐标残差时间序列中周年项的全球分布规律，建立了全球GPS测站垂向周年变化统计改正模型，分析了模型的改正效果，得出如下结论。

(1) 周年项是全球GPS测站中普遍存在的周期项，南北半球周年项振幅、北半球周年项初始相位随纬度的增加呈现先增后减的规律，南半球周年项初始相位的变化规律与北半球相反，南北半球中纬度地区测站周年项初始相位之间存在约180°的系统差，低、高纬度地区测站的周年项初始相位之间则不存在明显系统差。

(2) 本文建立的统计改正模型能削弱全球多数测站30%～50%的垂向坐标残差，可以精确分离测站的线性与非线性运动，说明建立统计改正模型是必要且可行的。

(3) 测站周年项的周期并不严格等于一周年，各个测站周年项的真实周期存在微小偏差。