GNSS垂向坐标时间序列周年项振幅的时变特征及成因分析

贾彦锋， 朱新慧， 孙付平， 肖凯， 柯能

1 信息工程大学地理空间信息学院， 郑州 450000 2 西安卫星测控中心航天器在轨故障诊断与维修重点实验室， 西安 710043

0 引言

GNSS坐标时间序列中包含着明显的周期性信号，对这种周期信号进行精确地建模和分析不仅有助于获取精确的基准站位置和速度，进而实现高精度地球参考框架的建立和维持，而且对于研究板块构造运动、地表质量迁移等地球动力学过程也具有重要的意义(Chen et al., 2013; 姜卫平等, 2013).处理坐标时间序列中周期信号的传统做法是利用具有常振幅和相位的谐波函数进行拟合，但事实上这些季节性信号的振幅和相位是随时间变化的(Blewitt and Lavallée, 2002; Murray and Segall, 2005; 张鹏等, 2007; Bennett, 2008; 明锋等, 2016; Klos et al., 2018).这种振幅和相位随时间变化的周期性信号被称作“准周期信号(Bennett, 2008)”或“时变季节性信号”.有研究表明，忽略季节性信号的时变特性会对基准站速度及其不确定度的估计精度产生影响(Blewitt and Lavallée, 2002; Bennett, 2008; Bogusz and Klos, 2016; Klos et al., 2018)，同时在残差序列中遗留大量的随机性季节信号(Davis et al., 2012)，导致残差时间序列出现相关性误差，这些相关性误差会使其频谱呈现出幂律噪声的特性(Blewitt and Lavallée, 2002; Williams et al., 2004; Langbein, 2008)，进而影响基于残差时间序列的误差频谱分析结果的正确性.除此之外这些时变季节性信号还会掩盖一些短期的非周期性信号(Bennett, 2008).因此，利用坐标时间序列进行精细化研究时，需要考虑振幅和相位的时变特性带来的影响.

针对这种时变季节性信号相关学者提出了不同的提取和分析方法.Davis等(2006)使用了分段连续线性多项式来表示季节性信号的振幅变化.Bennett(2008)提出了一种灵活的半参数方法来提取GPS坐标时间序列中的准周期信号.Davis等(2012)将时变季节性信号看作常振幅周期信号和零均值随机性分量的和，并利用卡尔曼滤波进行了信号的提取.Chen等(2013)将奇异谱分析引入到坐标时间序列时变季节性信号的提取中，并取得了较好的效果.伍浩琛等(2021)提出了一种基于Huber函数M估计的GPS坐标时间序列时变振幅周期信号的估计方法.此外有关学者还提出利用滑动最小二乘、小波分解和契比雪夫多项式函数等方法来提取时变季节性信号(Janusz, 2015; Klos et al., 2018).但是，已有研究针对的主要是时变季节性信号的提取，对于其产生的机制仅仅做了少量推测性的分析，认为其与降水、降雪等地球物理学因素有关(Murray and Segall, 2005; Bennett, 2008; Freymueller, 2009)，并未进行深入探究.此外，虽然环境负载和热膨胀效应已被证实是GNSS坐标时间序列周年变化的主要因素(姜卫平等，2018)，但从严格意义上讲，无法直接推测出两者也是周年项振幅变化的主要因素.而对振幅变化的量级大小和特征的研究更是很少涉及.

本文的主要目的是基于全球分布的GNSS测站垂向坐标时间序列，从坐标非线性变化的机制(环境负载、热膨胀效应以及GNSS数据处理策略)入手，深入分析季节性信号振幅变化的量级大小和时空分布特征，探究引起振幅时变特性的因素并建立相应的定量关系，为更好地理解和提取时变季节性信号提供参考.由于全球大部分测站的坐标时间序列中均存在周年项，且周年项的振幅往往大于半周年项等其他周期项，振幅的变化对于周年项的影响更大，因此本文将周年项作为研究对象.

1 数据选取与处理

1.1 数据选取

本文采用的数据主要有GNSS坐标时间序列数据、环境负载数据以及热膨胀模型数据.

GNSS垂向坐标时间序列数据是来自于斯克里普斯轨道和永久阵列中心(Scripps Orbit and Permanent Array Center, SOPAC)提供的全球IGS14框架下的“Clean_Detrend”单天解坐标时间序列(http:∥sopac.ucsd.edu/).SOPAC提了三个版本的坐标时间序列，分别是基于GIPSY-OASIS软件得到的JPL时间序列、基于GAMIT/GLOBK软件计算得到的SOPAC时间序列以及经过两者联合处理得到的COMB组合时间序列.本文在分析周年振幅变化的特征以及探究其与环境负载和热膨胀效应之间的关系时采用的是组合时间序列，因为组合解具有更高的精度和可靠性(Ferland and Piraszewski, 2009; Bock et al., 2019)；在研究GNNS数据处理策略对周年项振幅变化的影响时采用的是JPL和SOPAC的单一解.为了更加科学可靠地分析周年项振幅的变化特征和成因，需要对测站进行筛选，选取的原则是：坐标时间序列的时间跨度大于10年，数据质量较好且全球分布.根据以上原则，分别从组合解和单一解中选取了468和368个测站的数据用于研究，这些测站包含部分国际GNSS服务(International GNSS Service，IGS)的基准站和美国连续运行参考站(Continuously Operating Reference Stations，CORS)网的测站.

环境负载数据采用的是法国地球科学实验室与学院(École & observatoiredes sciences de la Terre，EOST)负载服务提供的地面测站环境负载在地球形心(Center of Figure, CF)框架下的三维位移时间序列(http:∥loading.u-strasbg.fr/ITRF/CF/)，该产品涵盖了ITRF2014基准站在内的约7500个全球分布的测站.其中大气负载(Atmospheric Loading, ATML)选取的是基于海洋对大气压力变化存在逆气压响应的假设，采用欧洲中期天气预报中心(European Centre for Medium-Range Weather Forecasts, ECMWF)提供的3 h分辨率的地表气压数据计算的ATMIB位移时间序列；水文负载(Continental Water Storage Loading, CWSL)来自于MERRA2模型计算的MERRA2时间序列；非潮汐海洋负载(Non-tidal ocean Loading, NTOL)采用的是根据ECCO2海洋底压模型计算的ECCO2位移时间序列.

热膨胀模型数据包含测站的温度数据和观测墩信息.其中温度数据由ECMWF ERA5 0.25°×0.25°全球温度格网的双线性内插得到；测站观测墩的信息(材质以及高度)来自于SOPAC提供的测站log文件(ftp:∥garner.ucsd.edu/archive/garner/docs/station_logs/).这些数据被用来计算测站的热膨胀位移.

1.2 时变周年项振幅的提取方法

为了研究周年项振幅变化的特征和成因，需要对其时变振幅进行提取.在振幅的提取过程中，本文利用了分段最小二乘的方法，其主要思想是：第一，利用巴特沃斯带通滤波器对待提取时间序列数据进行滤波得到只含有周年项作为主信号的时间序列(截止频率设置为0.75 cpy和1.25 cpy)，以排除半周年项等其他周期项的影响；第二，对滤波后的坐标时间序列整体做最小二乘解算得到周年项的平均振幅和初始相位；第三，将初始相位代入分段后的坐标时间序列利用最小二乘计算出该段的振幅，即假设周年项的相位不随时间变化.由于周年项的振幅变化不只局限于不同年份之间，在一个周期内波峰和波谷对应的幅值也具有很大的差异性，因此在对坐标时间序列进行分段时，不是以一个周期作为分段的标准，而是采用半个周期的分段原则，即以周年项的波峰和波谷为参照进行分段.

2 周年项振幅变化的时空分布特征

为了量化周年项振幅变化的量级，更好地分析高程方向周年项振幅的变化特点及成因，本文采用振幅的绝对偏差(Absolute Deviation, AD)和平均绝对偏差(Mean Absolute Deviation, MAD)来衡量周年项振幅的波动性：

(1)

2.1 全球分布特征

图1反映了周年项振幅平均绝对偏差的空间和大小分布情况.

从全球来看，测站周年振幅MAD的均值为0.88 mm，最大能够达到3.02 mm.在选取的测站中，周年项振幅的MAD主要分布在0.5～1 mm之间，占比60.47%，大于1 mm的比例为29.7%.经分析计算，以常振幅的方法进行拟合，最大会造成2～3 mm的周年项残留，因此在利用坐标时间序列进行精细化研究时周年项振幅的变化不可忽略.

从图1a中可以看出，南半球测站周年振幅的波动幅度整体要小于北半球.而研究表明，作为周年项主要来源的环境负载在南半球引起的测站垂向位移也要小于北半球(Davis et al., 2004; Jiang et al., 2013; Li et al., 2014, 2020; 龚国栋等, 2017; 常航, 2019)，这两个现象之间具有很好的一致性，说明环境负载很可能是周年项振幅变化的潜在因素之一.此外，北半球中高纬度地区的测站周年振幅波动性普遍较大.

2.2 时间分布特征

图2描述了周年项振幅绝对偏差随时间的变化.其中，蓝色的离散点代表的是测站在不同历元下周年振幅的绝对偏差，红色的曲线代表周年振幅在不同历元下绝对偏差的平均值(此处以周年项波峰和波谷处的历元表示不同波段对应的时间).从图中可以看出，全球测站周年振幅波动的平均水平是随时间变化的.从1995—2002年，周年振幅绝对偏差的全球平均水平呈下降趋势，在2002年下半年达到最低点，然后又开始上升，并在2014年附近达到了最高点，之后又开始在波动中减小.整体来看，全球测站周年项振幅的偏差呈现出近似周期性的变化，但是受限于坐标时间序列的时间跨度，并不能完全确定这一现象，仍需要数据量的积累和进一步的研究.同时可以看出振幅的变化不仅仅局限于不同年度之间，在年内，波峰和波谷对应的振幅也存在一定的差异.

图2 全球测站周年振幅绝对偏差随时间的变化

3 周年项振幅变化的成因分析

本文从引起坐标季节性变化的物理机制入手，采用坐标时间序列与不同因素引起的测站位移之间周年振幅变化的一致性，以及坐标时间序列经不同因素改正前后周年振幅平均绝对偏差的变化来分析它们对周年项振幅变化的影响和贡献.

振幅变化的一致性(Consistency of Amplitude Variation, CAV)定义为

(2)

其中，c和t分别表示不同因素引起的基准站位移与坐标时间序列之间周年振幅变化一致(相对于各自平均周年振幅同时变大或同时变小)的历元数和总历元数.

周年项振幅MAD的变化dMAD定义为

dMAD=MAD-MADcorrected，

(3)

其中，MAD和MADcorrected分别为改正前后坐标时间序列的周年振幅变化的平均绝对偏差.一致性越高，dMAD越大表明该因素对周年项振幅变化的影响和贡献越大.

3.1 环境负载对周年项振幅变化的贡献

在对EOST负载服务提供的环境负载时间序列和SOPAC提供的坐标时间序列进行历元的统一和对齐之后，一方面分别采用带通滤波器和分段最小二乘对两者进行周年项振幅的提取并计算两者变化的一致性，另一方面分别利用大气负载、水文负载、非潮汐海洋负载以及总的环境负载(ATML+CWSL+NTOL)对坐标时间序列进行改正，并计算改正前后周年项振幅MAD的变化，得到的结果如图3、图4、表1和表2所示.

图3 环境负载位移与坐标时间序列之间周年项振幅变化的一致性(左)和改正前后坐标时间序列周年项振幅MAD的变化(右)的空间分布，从上至下依次为大气负载、水文负载、非潮汐海洋负载和总的环境负载

图4 环境负载位移与坐标时间序列之间周年项振幅变化的一致性(a)和改正前后坐标时间序列周年项振幅MAD的变化(b)的大小分布

表1 环境负载位移与坐标时间序列之间周年项振幅变化的一致性的统计数据

表2 环境负载改正前后坐标时间序列周年项振幅MAD变化的统计数据

从振幅变化的一致性水平来看，各地球物理因素引起的负载位移和坐标之间的振幅变化有着很高的同步性，它们之间的平均一致性均能够达到50%以上.其中，平均一致性最高的是水文负载，可以达到64.28%.大部分测站的一致性水平在40%～80%之间，10%～20%的测站负载位移与坐标时间序列之间振幅变化的一致性能够达到80%以上，而位于安大略湖畔的KNGS和位于格陵兰岛的MSVG测站的水文负载与坐标时间序列之间振幅变化的一致性为100%，即完全同步.另外，在全球分布上，水文负载在北美和欧洲南部地区与坐标时间序列之间振幅变化的一致性明显高于其他地区，尤其是在格陵兰岛的东部沿海地区、地中海附近；非潮汐海洋负载在新西兰和美国中部具有较高的一致性.

从周年项振幅MAD的变化来看，经过改正，60%以上的测站周年项振幅的MAD出现了减小，减小的幅度大都在0～1 mm，少数测站能够达到1～2 mm.整体来讲，水文负载改正的效果最好，经过改正有70.95%的测站周年振幅变化的平均绝对偏差减小，平均减小程度为0.17 mm.在全球分布上，MAD的变化呈现出一定的区域特性.其中水文负载在北美的五大湖附近的改正效果十分可观的，该地区的测站经过水文改正周年振幅的MAD值出现了较大的减小；大气负载和非潮汐海洋负载的改正效果分布则比较均匀.

另外，我们发现，无论是在一致性水平上，还是改正前后MAD的变化上，经过叠加得到的总环境负载相比于单个因素并没有明显的提升，初步推测可能与由三种负载获取总的环境负载的方式有关.Clarke(2005)认为通过累加获取总环境负载的方法没有明确的理论支撑，盲目的累加可能会损失负载模型的性能.

总的来看，环境负载对于坐标时间序列振幅变化具有很大的贡献，是周年项振幅变化的原因之一.就不同的因素而言，水文负载整体对周年项振幅变化的影响和贡献最大，其次是大气负载和非潮汐海洋负载.至于水文负载在三种负载中贡献最大的原因，本文推测与引起水文负载的因素有关.水文负载与降水、降雪等气象因素具有十分密切的关系，而这些因素在不同的年份甚至在年内都存在较大的差异，这些差异会最终通过水文负载反映在坐标时间序列的振幅变化上.就不同的地区而言，北美洲和欧洲地区测站周年项振幅的变化受水文负载的影响较大.

3.2 热膨胀效应对周年项振幅变化的贡献

热膨胀效应引起的测站垂向位移包含两部分，一部分是温度变化引起的观测墩热弹性形变(Thermal Expansion of the Monument，TEM)，该部分采用的是王锴华(2019)提出的一种改进TEM模型；另一部分是所在基岩的热弹性形变(Thermal Expansion of the Bedrock，TEB)，该部分利用的是Fang等(2014)提出的统一球体热传导模型.对以上两部分相加得到热膨胀效应引起的测站垂向总位移，接着计算其与U方向坐标时间序列之间振幅变化的一致性和经过热膨胀效应改正后周年项振幅MAD的变化，得到结果如图5和图6所示.

图5 热膨胀效应位移与坐标时间序列之间周年项振幅变化一致性的空间(a)和大小(b)分布

根据计算结果，热膨胀效应与坐标时间序列之间振幅变化的平均一致性为60.43%，一致性最高可以达到96.67%，最低为11.11%.其中一致性在50%以上的测站占总测站数的80.96%，有24.51%的测站一致性水平在70%以上.从一致性的全球分布来看，北美东南部和欧洲中部的测站整体较高，尤其是在东欧地区，两者之间的一致性普遍在80%以上.

另外，经过热膨胀效应的改正有75.88%的测站周年项振幅的MAD出现了减小，最大减小幅度为1.3 mm，进一步说明了热膨胀效应是全球范围内测站周年振幅变化的重要原因之一.其中在北美的北部和东部振幅MAD的减小幅度整体较大.除此之外，在格陵兰岛西部和欧洲东部地区振幅的波动性也出现了比较明显的下降，说明在这些地区振幅的变化受热膨胀效应的影响比较大.

图7反映了通过叠加环境负载和热膨胀效应得到的综合改正模型的改正效果(dMAD/MAD×100%).从图中可以看出，经过环境负载和热膨胀效应的综合改正虽然有65%以上的测站周年项振幅的波动幅度出现了下降，但是下降幅度最大也只占到了改正前波动幅度的73%，仍有27%以上的残留，并且有接近35%的测站振幅的波动性不降反增.这一现象一方面可能是由环境负载模型和热膨胀效应模型的误差以及改正方法的缺陷(见前文分析)造成的，另一方面也说明了周年项振幅的变化并非完全由环境负载和热膨胀效应引起.而GNSS数据处理模型及策略很可能也是造成振幅波动的原因之一.

图7 环境负载和热膨胀效应的综合改正效果

3.3 GNSS数据处理模型及策略对周年项振幅变化的影响

为了探究坐标时间序列周年项振幅的变化是否与GNSS数据处理模型及策略有关，本文选取JPL和SOPAC两个不同的坐标时间序列数据，分析了两者之间周年振幅变化的一致性和平均绝对偏差的差异性.JPL和SOPAC时间序列采用的都是SOPAC数据库统一提供的GNSS原始观测数据和测站信息，两者的不同之处在于采用了不同的数据处理策略和模型，主要表现在(Bock et al., 2019)：第一，在处理策略上，JPL采用了精密单点定位的方法，而SOPAC采用的是分布式网平差的方法；第二，在海洋潮汐模型上，JPL采用了IERS 2010 convention提供的模型，而SOPAC采用的是FES04模型；第三，SOPAC相比于JPL在处理过程中进行了二阶和三阶电离层改正，同时采用了BERN 15参数光压模型.因此，通过计算两者周年振幅变化的一致性和平均绝对偏差的差异，即可反映出GNSS数据处理模型及策略对坐标时间序列周年振幅变化的影响.图8和图9分别给出了JPL与SOPAC坐标时间序列之间周年振幅变化的一致性和平均绝对偏差差值(MADJPL-MADSOPAC)的分布情况.

图8 JPL和SOPAC坐标时间序列周年项振幅变化一致性的空间(a)和大小(b)分布

图9 JPL和SOPAC坐标时间序列周年项振幅MAD差异性的空间(a)和大小(b)分布

根据计算结果，我们发现JPL和SOPAC坐标时间序列之间周年项振幅的变化并非完全一致，两者的平均一致性为70.70%，一致性在80%以上的测站仅有28.75%，在90%以上的更是不足10%，而且仍有部分测站一致性在50%以下.而从两者平均绝对偏差的差异来看，如果数据处理模型及策略对周年项振幅的变化没有影响，那么两者的差值应当满足均值为零的正态分布，并且其方差不应过大.但是从计算结果来看，两者差值的平均值为-0.13 mm，并不为0，而且在对其进行Shapiro-Wilk正态检验时，发现在0.05的显著性水平下，检验结果为两者差值不符合正态分布，即在95%的置信度下认为数据处理模型及策略会引起坐标时间序列周年项振幅的变化.

已有大量的研究表明，GNSS数据处理模型及策略的不完善会使坐标时间序列产生虚假的周期信号(Penna and Stewart, 2003; Deng et al., 2017; Allahverdi-Zadeh et al., 2016; 姜卫平等, 2018)，其中与周年项相近的信号的振幅很可能也是随时间变化的，这种变化会反映在坐标时间序列中，进而造成周年项振幅的变化.

3.4 讨论

结合周年项振幅波动幅度的空间分布，我们可以看出：

(1)在北美洲和欧洲地区，水文负载和热膨胀效应不论是在一致性上还是MAD的减小幅度上都要优于其他两种因素，说明在以上两个地区测站周年项振幅的变化受水文负载和热膨胀效应影响较大.引起水文负载和热膨胀效应的因素主要是降雨、降雪、温度等气象因素.以北美洲为例，北美洲地跨热带、温带、寒带，气候复杂多样，再加上内部山脉大都为南北或近似南北走向，其气候很不稳定，冬季时而寒冷，时而解冻(代祖玲, 1987).因此，由水文负载和热膨胀效应引起的测站位移的周年振幅也会发生较大的波动，进而造成坐标时间序列周年项振幅随时间发生变化.

(2)在亚洲及其他地区没有一种因素的一致性和改正效果明显优于其他因素，说明在这些地区周年项振幅的波动性是受多种因素综合影响的结果.

4 结论

本文通过对GNSS垂向坐标时间序列中周年项的时变振幅进行提取，计算了其变化的量级和大小，深入分析了全球测站周年项振幅的时变特征，然后在此基础上，从环境负载、热膨胀效应以及GNSS数据处理模型和策略入手，探究了引起周年项振幅变化的原因，确定了不同因素对周年项振幅变化的影响和贡献，得出了以下结论：

(1)在全球范围内，高程方向上的周年项振幅普遍呈现出一定的时变特性，其平均波动幅度在1 mm左右，其中有近30%的测站波动性大于1 mm，最大可以超过3 mm，忽略这种时变特性最大会造成2～3 mm的周年项残留.

(2)周年项振幅的波动性在空间上呈现出差异性的分布.其中北半球测站振幅的波动性整体要大于南半球，而且北半球中高纬度地区和整个亚洲地区的测站周年振幅波动性普遍较大.在时间分布上，全球测站周年振幅绝对偏差的平均水平是随时间变化的，并且呈现出近似的周期特性.

(3)环境负载和热膨胀效应是引起全球测站垂向坐标时间序列中周年项振幅时变特性的两个重要因素.其中由两者引起的测站垂向位移与GNSS垂向坐标之间周年项振幅变化的平均一致性能够达到60%左右，而且经过两者改正分别有68%和76%的测站周年项振幅波动性出现了下降.而在三种环境负载中，对周年项振幅变化贡献最大的是水文负载，其次依次是大气负载和非潮汐海洋负载.除此之外，GNSS数据处理模型及策略的不完善也是引起周年项振幅变化的原因之一.

致谢感谢许雪晴博士在基岩热膨胀效应计算方法上提供的帮助和建议，感谢SOPAC提供的GNSS测站坐标时间序列和测站信息，感谢EOST负载服务提供的环境负载数据和ECMWF提供的全球温度格网数据.