汪 浩，岳建平，向云飞，李乐乐

(河海大学地球科学与工程学院，江苏 南京 211100)

随着GPS研究的不断深入，已证实连续长期GPS坐标时间序列中不仅存在白噪声(white noise，WN)，还包含时空相关的有色噪声，如谱指数分别为-1、-2的闪烁噪声(flicker noise，FN)、随机漫步噪声(random walk noise，RW)及非整数谱指数幂律噪声(power law noise，PL)。文献[1—2]分析指出全球GPS站坐标时序的最优噪声模型为WN+FN。文献[3]认为通过空间滤波去除共模误差后可明显提高香港地区GPS基准站坐标时序的信噪比，最优噪声模型为可变白噪声(variable white noise，VW)+FN。实际上，GPS站点的噪声特性较复杂，各个站点上N、E、U分量的噪声特性可能各有不同[4- 6]。

除上述随机噪声外，GPS站坐标时序中还存在诸多非构造形变信号，其中固体潮、海洋潮与极潮改正在GPS数据后处理中多已被采纳，而地表质量负载包括大气压、非潮汐海洋及水文负载等引起的测站形变尚未修正。众多研究表明，地表质量负载形变与GPS周年运动呈强相关性。文献[7]定量计算地表质量负载引起的地表位移可解释全球IGS站高程方向约40%的周年项振幅。文献[8]利用海潮与地表质量负载造成的地壳非构造形变修正中国陆态网络GPS基准站坐标时序，可降低其垂向周年项振幅的37%，并使RMS值降低约1 mm。文献[9]发现地表质量负载虽能改变中国区域11个IGS站的噪声特性并降低U方向周年运动，但周日、半周日大气潮汐才是造成其周年运动的主要因素之一。本文将着重分析地表质量负载改正前后格陵兰岛区域GPS测站坐标时序的噪声特性、速度场及周期项振幅的变化情况，有利于深入探究格陵兰岛的地壳构造形变及运动规律。

1 数 据

1.1 GPS数据

为了更加精确地估计噪声特性、速度及周期等参数，一般要求GPS坐标时序的时间跨度不少于2.5 a。本文选取了ITRF2008框架下格陵兰岛区域12个连续GPS观测站2013年1月—2016年12月期间的单天坐标时间序列，数据来源于Nevada Geodetic Laboratory，采用GAMIT/GLOBK10.5对原始数据进行后处理。因格陵兰岛冰原覆盖占总面积的83.7%，气候严寒，交通不便，中南部点位密集且数据质量较好，故所选GPS站点分布于格陵兰岛中南部沿岸地区(如图1所示，图幅范围为20°W—60°W，59°N—75°N)。GPS坐标时序中普遍存在异常值及数据缺失等问题，首先采取三倍四分位距法(interquartile range，IQR)剔除粗差[10]，再利用三次多项式插值法对缺失值进行插补，直至获得连续高精度GPS站点坐标时序数据。

1.2 地表质量负载形变

采用的地表质量负载形变由法国斯特拉斯堡大学下属机构EOST Loading Service提供。大气压、非潮汐海洋与水文负载形变分别由欧洲中期天气预报中心(ECMWF)、海洋环流及气候估值(ECCO)模型与全球路面数据同化系统(GLDAS)等提供的数据计算获得[11- 12]，三者时间分辨率有所差异，大气压与水文负载为3 h，而非潮汐海洋负载则是24 h，空间分辨率相同，均为0.25°×0.25°。图2显示了MIK2站地表质量负载形变时间序列，从图2及其他站点计算结果可知，U方向地表质量负载形变最大，E、N方向次之，与文献[7—9]研究成果一致。针对所选站点而言，地表质量负载形变U方向振幅在10 mm以上，E、N方向振幅约为2～3 mm，在高精度GPS地壳形变应用中其影响不可忽略。

2 数据分析方法

2.1 极大似然估计

GPS坐标时序噪声分析可选用极大似然估计(maximum likelihood estimation，MLE)、最小范数二次无偏估计法(minimum norm quadratic unbiased estimation，MINQUE)等进行[13]。本文采用的CATS软件是通过MLE对GPS坐标时序的噪声特性、速度及周期项振幅等进行估计的[14]。基本原理为：针对给定的GPS日解坐标时间序列x建立线性函数模型，在估计噪声分量和线性函数待求参数时必须使其发生概率(l，即似然)达到最大，即

(1)

(2)

2.2 数据处理策略

GPS站坐标时序的噪声模型可为白噪声+有色噪声，但白噪声模型中VW只能反映GPS站坐标分量质量的好坏，不可作为最优噪声模型[5]，在此不予考虑，仅考虑WN与其他有色噪声的组合情况，利用CATS软件估计WN、WN+FN、WN+RW、WN+FN+RW及WN+PL 5种噪声模型的MLE值[9,15]，并以文献[4]提出的保守估计准则判断最优噪声模型，即当噪声模型的拟合自由度(具体见表1)相同时，MLE值较大则较优；不同时，则显著区分的阈值为：ΔMLE=拟合自由度差×2.6。

表1 不同噪声模型拟合自由度

3 地表质量负载改正前后GPS坐标时序特征变化分析

3.1 谱指数变化分析

采用WN+PL模型估计各个站点地表质量负载改正前后的谱指数，结果见表2。地表质量负载改正后各个站点N、E、U方向谱指数普遍减小，分别占GPS站总数的75.0%、91.7%、91.7%，且U方向谱指数范围发生明显变化，由-1.280 2～-0.749 5变为-1.273 5～-0.952 3，即更加接近-1，意味着闪烁噪声的成分随之增加，N、E方向的谱指数虽有一定程度减小，但其范围仍主要集中在-0.8～-0.5之间，表明幂律噪声可能是有色噪声的主要成分。

表2 地表质量负载改正前后谱指数

3.2 确定最优噪声模型

仅谱指数分析无法明确各站点的噪声特性，进一步利用MLE值判断最优噪声模型，所选站点的最优噪声模型主要为WN+PL，少量为WN+FN，具体表现为：地表质量负载改正前KSNB(N)、VFDG(E)、LYNS/MIK2/SRMP/WTHG(U)6个站点方向分量，以及地表质量负载改正后KSNB/MIK2(E)、DGJG/KELY/KSNB/MSVG/RINK/VFDG/WTHG(U)9个站点方向分量的最优噪声模型为WN+FN，其余均为WN+PL。统计可知：地表质量负载改正使得8.3%、25.0%、75.0%的N、E、U方向最优噪声模型发生改变，说明地表质量负载对U方向噪声特性影响最大，E、N方向次之，与前文结果一致。由图3可知，地表质量负载改正后，U方向WN+FN比例提升25.0%，甚至超过WN+PL的16.7%，证明其可增加U方向闪烁噪声的成分，对其他方向影响较小，与谱指数分析结果相同。

3.3 测站速度场变化

根据最优噪声模型估计测站速度场，所选12个GPS站点N、E、U方向的速度范围分别为-21.63～-11.83 mm/a、9.96～20.02 mm/a、0.82～17.45 mm/a，均值为-15.92、16.15、6.32 mm/a；经地表质量负载改正后，N、E、U方向速度范围变为-21.67～-11.87 mm/a、9.99～20.08 mm/a、0.81～17.44 mm/a，均值为-15.93、16.18、5.94 mm/a。由此可知：若不考虑非构造因素的影响，格陵兰岛地区每年约向东南方向偏移22.71 mm，向上抬升5.94 mm，进一步揭示了该地区的地壳运动规律。图4显示了地表质量负载改正前后站点三维速度差值，水平方向速度差值在0.1 mm/a以内；垂向差值为-0.8～0 mm/a，均为负值，表明地表质量负载对测站速度场影响较小(特别是水平方向)，平均能够降低高程方向速度约0.36 mm/a。

3.4 周期项振幅变化

GPS站点各个分量的周期特性可由功率谱分析法得出，如图5所示，KELY站存在明显的1 a项和0.5 a项，且U方向1 a项功率谱能量在地表质量负载改正后明显降低，分析其他站点时结果大致相同，鉴于SOPAC在进行GPS站坐标时序线性拟合时均考虑了1 a项与0.5 a项，本文功率谱分析结果与其一致。计算最优噪声模型下各方向分量的1 a项与0.5 a项振幅(如图6所示)，U方向周期项振幅(6 mm以内)远大于N、E方向(2 mm以内)。地表质量负载改正后，N、E方向1 a项振幅普遍增加，平均增幅约为18.5%、46.4%，0.5 a项振幅增加不明显，平均增幅约为2.9%、6.6%；U方向1 a项、0.5 a项振幅降低明显，平均降幅44.1%、14.2%。综上，地表质量负载引起的非构造形变是造成格陵兰岛GPS站垂向季节性信号(特别是1 a项与0.5 a项)的重要原因之一，却增强水平方向的季节性运动，但不排除受水平方向季节性运动较弱和其他系统误差等影响。

4 结 论

本文利用地表质量负载形变修正格陵兰岛区域12个GPS站4年的坐标时间序列，主要探讨了地表质量负载对GPS坐标时序噪声特性、速度场及周期项振幅的影响，得出了以下结论：

(1) 地表质量负载改正后，GPS站坐标时序噪声特性发生一定的改变，谱指数普遍降低，且U方向谱指数更接近-1，即闪烁噪声成分有所增加；最优噪声模型仍以WN+PL和WN+FN为主，但N、E、U方向8.3%、25.0%、75.0%的测站噪声模型发生改变，特别是U方向WN+FN模型比例提高25.0%，比WN+PL高16.7%，与谱指数分析结果一致。

(2) 地表质量负载对测站N、E方向速度影响较小(±0.1 mm/a以内)，但可平均降低U方向速度约0.36 mm/a，所选格陵兰岛区域地壳运动规律为：东南方向偏移22.71 mm/a，向上抬升5.94 mm/a。

(3) 地表质量负载形变修正后可降低U方向1 a项、0.5 a项的44.1%、14.2%振幅，相反，加大了N、E方向周期项振幅，1 a项振幅提升18.5%、46.4%，0.5 a项振幅增加较少，但水平方向周期项振幅因量级(2 mm以内)较小，可能受系统误差影响较大，有待进一步研究及验证。