尼泊尔地区GNSS站点是大陆俯冲带站点最为密集的区域之一，国内外学者利用连续GNSS测量研究结果揭示喜马拉雅造山带在不同空间尺度和不同孕震阶段的变形特征，对尼泊尔地震的同震滑移、震间闭锁和区域流变结构进行深入研究^[1]。然而，对于尼泊尔地区而言，以往的GNSS数据分析大多局限于孕震周期的某一阶段(震间、同震和震后)，缺乏对多孕震阶段时间序列的统一联合分析。事实上，尼泊尔地区的GNSS时间序列包含了丰富的信号特征，既有构造因素引起的变形信号(地震周期变形)，也有环境因素引起的变形信号(如降水变化、冰后回弹等机制)，这些变形信号相互耦合，难以单独进行观测约束。因此，联合2015-04-25尼泊尔M_W7.8地震前后的GNSS观测数据，采用统一的误差分析和噪声分析手段，对于揭示喜马拉雅中段地区地表的长期变形特征尤为必要。

本文收集了尼泊尔近场密集GNSS站点数据，通过研究完整的时间序列数据处理流程，精确获得地表形变位移，并模拟GNSS测站位置的时间变化。数据处理过程考虑测站震间速度、同震位移和地震后指数松弛，研究结果能广泛用于推断尼泊尔M_W7.8地震主要的地震源机制和断层几何尺度。

1 数据收集与解算

本文数据主要由连续观测站和流动观测站两部分组成，其中连续站包含中国境内的中国大陆构造环境监测网络(陆态网络)基准站(CMONOC)、尼泊尔境内的连续观测站(NEPAL)及中国地震局震后加密的连续观测站(CEA)，流动站主要是由中国地质大学(武汉)(CUG)和中国地震局地震研究所科考队野外采集，测站统计情况见表 1。

图 1展示了地震震中、GNSS站点分布情况及几种测站建站类型，其中连续观测站以深锚式天线墩和水泥天线墩为主，流动观测站观测方式有强制对中(岩石地钉、水泥天线墩)和架设脚架(地表埋石)两种方式。本文收集震后流动站5期观测数据(2015-05、2015-08、2015-11、2016-04、2016-12)，每期观测36 h，保证至少有一个整日观测时段。观测设备包括Trimble 5700、R7、NETR8、NETR9及Topcon NET-G3A等多种型号，外业观测完成后需对数据作预处理、格式整理和数据质量检测。天线高的量取方式有斜高和直高两种，天线外观尺寸不同量取的位置也有差异; 数据解算中天线代码由5个字符组成，前2个字符表示量取方式，后3个字符表示量取位置，如Trimble 5700、R7的天线垂高和斜高量取代码分别为DHARP和SLBGN。

本文采用GNSS软件GAMIT/GLOBK处理收集的数据，为减少解算误差，在GAMIT基线处理中加入日/月极潮、固体潮、海潮等改正，同时对卫星轨道、大气压模型、水汽分布、多路径效应及天线相位中心进行改正; GLOBK平差处理中合并IGS全球解文件(igs1、igs2、igs3、igs4、igs5、igs6、igs7)，以获得这些站点在ITRF2008框架下的坐标和中误差。

2 原理与方法 2.1 时间序列模型拟合

GNSS测站观测的时间序列可用如下线性多项式来表示^[2]：

$ \begin{gathered} y\left(t_i\right)=a+b t_i+c \sin \left(2 \pi t_i\right)+d \cos \left(2 \pi t_i\right)+ \\ e \sin \left(4 \pi t_i\right)+f \cos \left(4 \pi t_i\right)+ \\ \sum\limits_{j=1}^{n_h} g_j H\left(t_i-T_{\mathrm{eq}}\right)+\sum\limits_{j=1}^{n_j} h_j H\left(t_i-T_{\mathrm{post}}\right) t_i+ \\ \sum\limits_{j=1}^{n_g} k_j \exp \left(-\left(t-T_{\mathrm{post}}\right) / \tau_j\right) H\left(t_i-T_{\mathrm{post}}\right)+v_i \end{gathered} $

(1)

式中，t_i(i=1, …, N)为观测时间，系数a、b、c、d、e、f分别为测站时间序列起始位置、运动速率、年周期运动和半年周期运动，H(*)为阶跃函数，系数g_j为地震发生时刻T_eq发生的突变值，n_g表示有n个地震发生(n_h、n_k含义相同)，h_j为地震发生之后测站运动速率的变化，k_j为震后T_post时刻松弛或滑移位移的大小，τ_j为测站震后松弛时间常数，v_i为测量误差(是一个与时间无关的量，即期望E(v_i)=0)。

对于尼泊尔地震，发生时刻(2015-04-25 M_W7.8，2015-05-12 M_W7.3)和震后时间是已知的，也就是式中T_eq、T_post已知，n_g=n_h=n_k=1，τ_j的最佳估值采用模型不拟合度最小方法确定^[2]。

2.2 周期性分析

GNSS时间序列可以展现地表在时间变化下的特性，如因板块运动的线性效应及周期性效应等。在周期性效应中，最为常见的为年周期与半年周期，该效应一般被认为由3种原因造成^[3]：1)由日、月潮引力造成的影响; 2)由太阳为热源造成的一系列周期变化，如大气压力、非潮汐所引起的海平面变化、地下水效应及热膨胀效应、季节性的降雨等; 3)由GNSS解算模型改正误差引起的周期性变化，如多路径效应、天线相位中心改正等。Dong等^[3]分析全球429个GPS连续站周期性效应显示，水平方向的年周期振幅范围在1~3 mm，高程方向则在4~10 mm。

在GNSS数据解算的时候，已经充分考虑第1)和第3)两种原因造成的扰动，因此时间序列中的周期性效应更多来源于第2)种因素的影响。本文采用GARCE时变卫星重力场来计算水文荷载变化引起的地表位移，并用以改正GNSS时间序列中的周期项。

利用GRACE反演地表三分量位移变化(Δe，Δn，Δh)的计算公式为^[4-6]：

$\left\{\begin{aligned} \Delta e= & \frac{R}{\sin \theta} \sum\limits_{n=1}^{\infty} \sum\limits_{m=0}^n \bar{P}_{n m}(\cos \theta) m\left[-\Delta \bar{C}_{n m} \cdot\right. \\ & \left.\sin (m \varphi)+\Delta \bar{S}_{n m} \cos (m \varphi)\right] \frac{l^{\prime}}{1+k^{\prime}}{ }_n \\ \Delta n= & -R \sum\limits_{n=1}^{\infty} \sum\limits_{m=0}^n \frac{\partial}{\partial \theta} \bar{P}_{n m}(\cos \theta)\left[\Delta \bar{C}_{n m} \cdot\right. \\ & \left.\cos (m \varphi)+\Delta \bar{S}_{n m} \sin (m \varphi)\right] \frac{l_n^{\prime}}{1+k_n^{\prime}} \\ \Delta h= & R \sum\limits_{n=1}^{\infty} \sum\limits_{m=0}^n \bar{P}_{n m}(\cos \theta) \cdot\left[\Delta \bar{C}_{n m} \cdot \cos (m \varphi)+\right. \\ & \left.\Delta \bar{S}_{n m} \cdot \sin (m \varphi)\right] \cdot \frac{h^{\prime}{ }_n}{1+k^{\prime}} \end{aligned}\right. $

(2)

式中，l′_n、h′_n和k′_n为负荷勒夫数，$\Delta \bar{C}_{n m} $、$\Delta \bar{S}_{n m} $为GRACE的时变球谐系数，θ为地心余纬，$\bar{P}_{n m} $为完全规格化的勒让德函数，R为地球平均半径。

本文使用CSR提供的RL05(Level-2 Release-05)月重力场产品GSM，时间跨度为2011~2016年。由于GNSS时间序列数据中没有删除大气及非潮汐海洋负荷影响，所以在GRACE的GSM数据中加回了GAC产品，以保持与GNSS数据一致。具体数据处理步骤为：1)使用SLR计算的C₂₀替换原C₂₀项^[4]，修正GRACE轨道系统误差; 2)GRACE的所有产品都基于CM框架，不提供一阶项，Swenson等^[4]利用重力场球谐系数之间的关系及海底压强模型约束得到一阶项并验证了其有效性，本文采用该一阶项结果; 3)采用高斯平滑滤波方法去除GRACE数据截断误差，选取平滑半径为400 km^[5]; 4)GRACE数据同次不同阶之间存在相关性，需要进行去相关性处理^[6]。

2.3 共模误差分析

GNSS时间序列共模误差采用PCA空域滤波方法^{[3, 7]}。由于存在各种未知原因，GNSS实际数据中不可避免地存在大量缺失数据^[8]，通常采用插值补齐缺失数据，但这可能造成数据失真，因此本文采用存在缺失数据的GNSS时间序列主成分分析法。

对于区域网GPS站点日解时间序列，每个坐标分量(E、N和U三个方向)的残差坐标序列(去趋势和去均值)都可构造一个m×n维矩阵X，其中n表示测站数，m表示历元数，且满足关系式m≥n。矩阵X的列表示某一测站单分量的所有历元数据，行表示某一历元有观测的所有测站集合。采用存在缺失数据的PCA方法，将基于Dong等^[7]的经验正交函数改造如下：

$ \begin{gathered} a_k\left(t_i\right)=\sum\limits_{j=1}^{j \in S_i} X\left(t_i, j\right) v_k(j)+ \\ \sum\limits_{j=1}^{j \notin S_i} X\left(t_i, j\right) v_k(j), k=1, 2, \cdots, n \end{gathered} $

(3)

式中，时间序列在t_i历元拥有观测数据的测站构成子集S_i，则S_i集合包含的最大测站数为n，v_k是协方差矩阵B的第k个特征向量，协方差矩阵B的特征值呈降序排列(λ₁≥λ₂≥…≥λ_n)。第一主成分对应着最大的特征值λ₁，具有整个区域最多的信息，往往反映整个网的共同变化趋势。

2.4 噪声分析

对GNSS时间序列进行噪声分析主要有2个目的^[9]：1)获得合理的参数估计误差。时间序列不是单纯的高斯白噪声，而是包含与时间相关的噪声，直接采用最小二乘求解会掩盖部分误差; 2)本文收集的GNSS测站有4种不同的布点方式，噪声分析可以探究GNSS建站方式对坐标时间序列噪声模型的影响。

GNSS时间序列的功率谱是经过傅立叶变换的平方除以采样点数N获得的，定义为单位频带内的信号功率，反映信号功率随频率的变化情况，即信号功率在频域的分布状况^[9]：

$ P_x(f)=P_0\left(\frac{f}{f_0}\right)^k $

(4)

式中，f为频率; P₀和f₀为常数; k为待求的谱指数，取值范围－3~1，当－3＜k＜－1时称为典型的布朗模型(classical Brownian motion)，当－1＜k＜1时称为分式白噪声(fractional white noise)。

3 时间序列分析结果 3.1 周期性变形特征

本文利用GRACE改正GNSS时间序列中由水文造成的部分周期影响。经过GRACE改正后，GNSS坐标时间序列3个分量加权均方根分别平均降低26.4%、15.5%、45.0%。

图 2(单位mm)为GNSS测站三分量的周年项、半周年项振幅和相位，图中线段长度代表振幅大小，方向表示相位。相位代表了周期项波峰出现的时间，尼泊尔近场GNSS测站E向周年项波峰和半周年项波峰分别出现在夏季和冬季附近; N向周年项波峰和半周年项波峰出现在春夏之交和秋冬之交; U向周年项波峰和半周年项波峰大多出现在初春和冬、夏季。对比GRACE水文荷载时间序列中的周年项和半周年项的相位和振幅发现，GRACE比GNSS振幅要小且存在相位差，说明除了水文因素影响，GNSS时间序列还存在其他物理因素作用。图 3以KKN4和SNDL测站为例对比了GNSS时间序列中的周期信号和GRACE模型的季节性水文荷载时间序列，亦可以得出相同结论。

3.2 共模误差特征

图 4为利用存在缺失数据的PCA方法求得的GNSS站点前三主成分空间向量和各特征值贡献率。结果显示，第一主成分对坐标序列各分量的贡献率分别为61.91%、53.91%和53.88%; 第二和第三主成分空间向量随空间变化波动较大，不能反映整个网的共同特征。第一主成分空间向量很平稳，有非常一致的空间分布模式，具有整个区域最多的信息，反映整个网的空间共同变化趋势，故采用第一主成分计算共模误差：

$ \varepsilon\left(t_i, j\right)=a_1\left(t_i\right) v_1(j) $

(5)

式中，a₁为第一主成分, v₁为协方差矩阵的第一特征向量, t_i为历元，j为测站数。

3.3 噪声特性

图 5为KKN4和SDNL两个测站的噪声时间序列及其功率谱，可以看出，当频率为1 Hz和2 Hz时功率谱存在明显波动，表示为周年和半周年振幅。

图 6为研究区GNSS站点的噪声类型统计结果，可以看出，谱指数多数位于－1~0之间，表明多数站点既有闪烁噪声特性也有白噪声特性。选定最佳的噪声类型WN+FN，构建随机模型如下：

$ \boldsymbol{D}=\boldsymbol{\sigma}_{\mathrm{w}}^2 \boldsymbol Q_{\mathrm{w}}+\boldsymbol{\sigma}_{\mathrm{F}}^2 \boldsymbol Q_{\mathrm{F}} $

(6)

式中，$ \boldsymbol{\sigma}_{\mathrm{w}}^2$、$\boldsymbol{\sigma}_{\mathrm{F}}^2 $分别表示白噪声、闪烁噪声的方差，Q_w、Q_F分别表示白噪声、闪烁噪声协因数矩阵。

采用方差分量估计算法定量计算闪烁噪声和白噪声大小，同时估计震间速度、同震和震后位移。

表 2为噪声模型计算测站的N分量与E分量震间速度值、尼泊尔地震同震及震后位移值，上行为经过噪声模型WN+FN处理分析的结果，下行为WN噪声模型分析结果，后面的常数是两者误差的倍数。统计发现，WN+FN模型得出的震间速度、同震和震后位移中误差分别约为WN模型结果的6~9倍、4~7倍和2~4倍，WN模型低估了误差值。图 7为两类观测站震后时间序列。

3.4 同震和震后位移

利用时间序列分析结果获得同震和震后位移，由图 8可知，地震造成中国藏南地区在水平方向产生永久形变^[10]，同震垂直位移场尼泊尔境内KKN4站和NAST站处于隆升区，分别抬升125 cm和61 cm，而尼泊尔境内CHLM站和西藏聂拉木县WT11站处于下降区，分别下降57 cm和12.7 cm，此外J041站下降7.8 cm，TNG3站下降7.5 cm。利用指数模型对尼泊尔地震震后20个月的GNSS数据进行分析^[11]可知，震后松弛时间常数τ_j=80，计算每个GNSS测站的震后松弛位移(图 9)发现，尼泊尔地震震后形变主要发生在北部区域，且最大震后位移形变可达约60 mm。

4 结语

本文系统研究GNSS时间序列后处理分析流程，利用GAMIT/GLOBK软件包解算2011~2017年尼泊尔近场58个GPS连续观测站和41个震后流动观测站，获得ITRF2008高精度坐标序列，并以此为基础，研究GRACE重力数据改正GNSS时间序列的周期项、共模误差和噪声特性。得出以下结论：

1) GRACE改正GNSS数据，使得时间序列3个分量的WRMS平均降低15.52%、26.41%和45.06%;

2) 采用存在缺失数据的PCA方法获得第一主成分并计算共模误差，3个分量有效率达到61.91%、53.91%和53.88%;

3) 频谱分析确定WN+FN为最佳噪声模型，新的噪声模型得出的震间速度、同震和位移中误差是常用的WN模型结果的约6~9倍、4~7倍和2~4倍，即采用WN模型得出的中误差低估了真实值。

4) 通过对尼泊尔地震近场GNSS站点时间序列进行分析可知，结合Mencin等^[11]提供的震后余滑模型，尼泊尔地震震后形变主要发生在北部区域，且最大震后位移形变可达约60 mm。