GPS连续站观测作为主要的形变监测手段之一，可以有效捕获地震周期运动，涵盖了地壳应力积累引起的震间形变，断层突然破裂产生的同震形变，以及地壳和上地幔逐步恢复到稳态的震后形变过程，对于揭示地震孕育发生和调整起到重要作用(Bock et al., 1993; Bock and Melgar, 2016; Larson et al., 1997;Gan et al., 2007; Liang et al., 2013; Zhao et al., 2017).研究GPS时间序列的震后形变具有非常重要的意义，一方面，建立合适的震后形变模型有助于获得更准确的震间形变和同震形变，准确的震间和同震形变是研究板块运动、应力应变积累、地震孕育和破裂过程的基础(Bürgmann and Thatcher, 2013).另一方面，由于地表捕获的震后形变的时空变化和地壳深部物理状态直接相关，研究GPS观测的震后形变有助于研究震后形变机制、区域岩石圈流变学结构和断裂带的力学性质(Pollitz et al., 2000; Bürgmann and Dresen, 2008; 谭凯等, 2007; 朱守彪和蔡永恩, 2009; 张晁军等, 2009; Panet et al., 2010; 孟国杰等, 2016).此外，GPS跟踪站的震后形变参数也是维持动态参考框架的重要组成部分，如ITRF2014动态参考框架已提供IGS站的震后形变参数供用户使用(Altamimi et al., 2016).

在GPS时间序列震后形变估计方法上，许多研究(Gonzalez-Ortega et al., 2013;Wang and Fialko, 2014；Ding et al., 2015；Tobita, 2016)采用分步法计算各模型参数，首先利用震前数据估计震前速度，震前和震后测站的坐标差估计同震变形，然后从GPS时间序列中扣除震前速度和同震形变得到震后形变.这种方法计算的速度忽视了大量震后的观测数据，而同震形变受随机误差影响且包含了部分震后形变，进而影响震后形变参数.因此，对GPS时间序列进行整体建模估算各模型参数更严谨(Altamimi et al., 2016).另一方面，GPS测站观测到的震后形变由相同的驱动机制引起，必然存在内在联系，采用逐站点逐分量估计的策略则忽略了它们的空间相关性.Savage和Svarc (2009)和Barbot等(2009)基于震后形变的空间相关性，采用主成分分析法(Principal Component Analysis, PCA)提取共性的震后形变时间演化过程，但Savage和Barbot的研究仍采用分步法，没有避免震后形变和速度、同震形变参数的相互干扰.

基于现有的参数估计方法，本文提出迭代PCA参数估计方法，该方法兼顾震后形变衰减特征的空间相关性和整体建模的思想，并通过蒙特卡罗方法合成模拟数据验证了迭代PCA方法在参数估计上的优势，最后基于该方法研究了2016年M_W7.8新西兰Kaikoura地震的震后形变特点，定量分析了震后形变对地表速度的影响.

原始GPS时间序列包含了站点速度，天线移动或设备更换引起的阶跃，同震和震后形变，年周期和半年周期项，还有一些未模型化的瞬态形变和噪声.考虑到历元间分量的相关性非常小(Zhang, 1996)，t_i历元的坐标时间序列可以模型化为(Nikolaidis, 2002)：

(1)

其中，t_i表示以t₀起的历元，a是初始历元t₀的位移值，b表示震间速度，c, d, e和f分别表示年周期和半年周期项系数，g为T_g历元的阶跃，h为T_h历元的同震形变，基于地壳线性黏弹性流变机制的指数函数表示T_k历元振幅为k衰减常数为τ的震后形变，ε为观测噪声，H代表Heaviside阶梯函数.

常用的震后形变模型还有基于余滑的对数模型 (Marone et al., 1991)，它们被广泛应用于1992年M_W7.2 Landers，1999年M_W7.1 Hector Mine，2004年M_W6.0 Parkfield，2010年M_W7.2 El Mayor-Cucapah等震后形变研究中(Shen et al., 1994; Savage and Svarc, 1997, 2008; Savage and Langbein, 2009; Gonzalez-Ortega et al., 2013).

对有n个测站，时间尺度为m天的区域GPS观测网的时间序列，假定m>n, 构建一个m×n的观测矩阵X，这里每一行代表同一历元所有站的三个方向的坐标分量，每一列代表一个站点所有历元的坐标分量，每个站点N、E、U方向共占据三列，每一列均需要去均值化.协方差矩阵定义为B=X^TX，对协方差矩阵B进行特征值分解得到B=VΛV^TX, V是相互正交的特征向量，Λ有k个非零的特征值，特征值λ_k和特征向量e(k)根据特征值降序排列.

第k主成分为

(2)

利用PCA分析震后形变信息，可以提取区域内共性震后形变，抑制噪声和其他因素影响，越大的主成分对观测值贡献越大(Guttman，1950)，第一主成分对应整个观测网的共性震后形变(Savage and Svarc, 2009; Barbot et al., 2009)，贡献小的主成分可能与小尺度或局部性的非构造运动或噪声有关.

迭代PCA方法在PCA的基础上融合了整体建模的方法，并通过迭代避免各参数的耦合, 它的核心思想分为三个部分：一是顾及区域内GPS站点的震后形变的内在联系，对所有站的震后形变信息进行PCA分析，搜索整网最优衰减常数τ；二是对各站点坐标时间序列NEU方向分别整体建模，避免PCA方法分步估计参数的缺陷；三是由于震后形变衰减常数和震后形变、同震形变和速度参数存在一定耦合，通过迭代逐步分离参数间的相互影响.迭代过程以衰减常数τ为突破点，衰减常数τ不再变化时，其他参数相应稳定.

迭代PCA方法处理区域GPS网时间序列的流程(见图 1)为：首先，对数据预处理，包括粗差探测和空缺值内插，采用滑动窗口法探测并剔除孤立粗差，窗口尺寸设为180天，将残差超过3倍的四分位距和中值之差作为粗差(Nikolaidis, 2002)；PCA分析时，如果某一历元某一个站存在空缺需要将该历元所有站数据删除，这将浪费大量有价值数据，为尽可能保留多有效观测数据需要预先对缺失数据插值，具体方法参照Dong等(2006).其次，逐站NEU方向分别整体建模计算初始模型参数，根据模型参数扣除速度、阶跃、同震形变、年周期半年周期项，仅保留震后形变和未模型化的残差.然后，对所有站各方向的震后形变进行PCA分析，得到各主成分.最后，对第一主成分进行模型拟合，估算最佳衰减常数τ，并将新的τ代入各站各分量，重新估算新模型参数.多次迭代，直到各参数不再变化时，退出循环.各站共性的震后形变时间演化过程归因于共同的震后形变驱动机制, 而第一主成分认为是震后形变在时间域上共同的衰减过程，对应的系数表示不同的空间响应(Savage and Svarc, 2009; Barbot et al., 2009).由于震后形变模型是非线性的，本文采用试错搜索法(Dong et al., 2006)试算一定范围内所有τ值，产生最小卡方值(观测值和模型值残差平方和)的τ作为最优衰减常数，再基于最优τ值进行线性最小二乘估计其他模型参数.

图 1

Fig. 1

图 1 迭代PCA方法流程图 Fig. 1 Flow chart of the iterative PCA method

为了验证迭代PCA方法的优势，本文利用模拟数据对比该方法与分步PCA方法和单站整体建模方法在参数估计上的表现.根据GPS时间序列模型(式(1))采用蒙特卡罗方法合成模拟数据，共合成1000组数据，每组5个GPS时间序列，其中，时间范围为2010—2017年，时间序列参数分别在给定范围内随机采样，初始值a的采样范围0~1000 mm，速度b为1~50 mm·a^-1，周期项c, d, e, f为0.1~5 mm，同震形变h为1~150 mm，震后形变振幅k为1~100 mm，衰减常数τ为1~200天，观测噪声ε设计为方差为5 mm的随机高斯白噪声.为了简化问题，模拟数据仅设计一次地震事件，发震时刻2011.0027，不含其他阶跃.将随机采样的衰减常数、速度值、同震和震后形变参数作为已知真值，统计三种方法估算的模型参数与真值的差值RMS(表 1)，可以看出分步PCA方法得出的同震形变偏差较大，单站整体建模的衰减常数偏差较大，迭代PCA方法可以扬长避短，整体表现最稳健可靠.其中，估算的衰减常数比整体建模方法提升约83%，估算的同震形变比分步PCA方法提升约94%，估算的震后形变和速度比现有方法分别提升42%和47%以上.

表 1 Table 1

表 1 三种方法估算的参数误差 Table 1 RMS of parameters estimated by three methods 迭代PCA 分步PCA 单站整体建模 衰减常数(天) 3.0715 3.0272 18.5387 震后形变振幅(mm) 0.4572 0.7894 0.8675 同震形变(mm) 0.5208 8.5195 1.0469 速度(mm·a-1) 0.0874 0.2623 0.1644

表 1 三种方法估算的参数误差 Table 1 RMS of parameters estimated by three methods

为了进一步分析衰减常数、同震和震后形变，速度参数之间的耦合性，本文统计了两两参数之间的关系，图 2展示了存在相关性的参数误差分布.从图中可以看出当衰减常数τ增大，三种方法估算τ的误差随之增大，单站单分量整体建模方法的τ值的偏差和增幅最大(图 2a)；衰减常数τ越小，震后初期衰减越快，三种方法估算的同震形变的误差越大，分步PCA方法采用坐标差计算同震形变吸收了部分震后形变而整体偏大，尤其当衰减常数较小时，偏差非常大(图 2b)；衰减常数τ越大，三种方法计算的震后形变振幅估算偏差越大，迭代PCA方法最优，分步PCA方法在衰减常数较小时表现略差，而整体建模方法在衰减常数较大时表现略差(图 2c)；分步PCA方法利用震前数据估算速度，因此与衰减常数无关，而迭代PCA与整体建模方法估算的速度随着衰减常数减小而减小(图 2d); 速度偏差与震后形变偏差具有非常明显的负相关关系，迭代PCA方法通过迭代的方式一定程度上降低了速度和震后形变的耦合性(图 2e).参数偏差的分布从侧面揭示了估算的参数误差主要与震后形变衰减形态有关，震后初期衰减越快，估算的同震形变偏差越大，但震后形变与速度越容易分离，速度与震后形变偏差越小；反之，震后初期衰减越慢，同震形变偏差越小，衰减常数、震后形变振幅和速度偏差越大.

图 2

Fig. 2

图 2 模拟实验的参数关系图 (a)衰减常数与衰减常数估算偏差的关系；(b)衰减常数与同震形变估算偏差的关系；(c)衰减常数与震后形变偏差的关系；(d)衰减常数与速度偏差的关系；(e)震后形变估算偏差与速度估算偏差的关系.红色实点是迭代PCA方法对应值，绿色实点是分步PCA方法对应值，蓝色实点是单站单分量建模方法对应值. Fig. 2 Relationships between parameters estimated in the simulation data experiment (a) Decay constant and decay constant bias; (b) Decay constant and coseismic deformation bias; (c) Decay constant and postseismic amplitude bias; (d) Decay constant and velocity bias; (e) Postseismic amplitude bias and velocity bias. Red dots denote the parameters estimated by iterative PCA method. Green dots are two-step PCA method. Blue dots are single mode modeling method.

新西兰地区位于太平洋和澳大利亚板块边界，太平洋板块相对于澳大利亚板块西西南向移动，速率大约40 mm·a^-1(USGS)，太平洋和新西兰南岛的澳大利亚板块浅部逆冲倾滑造成了2016年11月13日新西兰Kaikoura M_W7.8地震.该地震震中位于Hope断层以南20 km处，破裂始发于北坎特伯雷，随后破裂向北扩展约150 km直到马尔堡, 进而引发多条断裂带破裂(Litchfield et al., 2016; Bradley et al., 2017).破裂区的地质构造复杂，至少有12条断层发生破裂，这些断裂有左旋的，右旋的，也有逆冲的，是记录到的最复杂的地震之一(Shi et al., 2017).该次地震是自1855年M_W8.2-8.3威拉拉帕湖地震以来记录到的最大地震(GeoNet，2016)，整个新西兰地区均有震感，南岛的北部地区大面积受到破坏，并在南岛和北岛的东海岸引发海啸(Hamling et al., 2017).

本文选取了SOPAC处理的新西兰境内37个GPS连续站原始坐标时间序列(http://sopac-ftp.ucsd.edu/pub/timeseries/measures/ats/Global/)，数据时间跨度为各站起始观测至2017年6月10日.由于震后早期以余滑为主(Segall，2010；孟国杰等，2016)，因此，M_W7.8新西兰Kaikoura地震的震后形变模型采用对数模型.利用迭代PCA方法分析了37个GPS站坐标时间序列，迭代了7次后退出循环，τ依次为11、9、7、6、5、4、4天.图 3a显示了PCA分析的前三个主成分，第一主成分占所有信号的82.57%，可以看作共性的震后形变信号，其他主成分主要是误差和未模型化的残留信号，图 3b为第一主成分和模型拟合值，衰减常数τ为4天的模型很好的拟合了第一主成分.震后形变方向和振幅大小(图 3c和3d)显示：水平方向上，南岛北东部和北岛最南部震后形变较大，南岛北西部-南岛北东部-北岛最南部的震后形变有逆时针旋转的形态；垂直方向上，南岛大部分地区和北岛最南部发生抬升，离震中最近的cmbl站抬升幅度最大，与震中北东方向的地表破裂区发生隆起(Hamling et al., 2017; Shi et al., 2017)基本一致，北岛北部大部分地区发生小幅下沉.

图 3

Fig. 3

图 3 GPS震后形变的前三主成分和第一主成分系数 (a)前三主成分(蓝线是第一主成分，橙线是第二主成分，黄线是第三主成分)；(b)第一主成分以及模型值(黑点是第一主成分，红线是模型值，蓝线是发震时刻)；(c)各站NE方向的第一主成分系数(黄星代表震中，红色和蓝色箭头代表不同比例尺的系数大小)；(d)各站U方向的第一主成分系数(黄星代表震中，红色箭头代表系数为正，蓝色箭头代表系数为负). Fig. 3 First three principal modes and the coefficient of first principal mode of GPS postseismic deformation (a) First three principal modes. The blue line is first mode, orange line is second mode and yellow line is third mode. (b) First principal mode and its model value. Dark dots denote first principal mode. Red line denotes model value. Blue line marks the time of earthquake occurred. (c) The coefficients of first principal mode in NE direction. Yellow star is the epicenter. Red and blue arrows are the coefficients on different scales. (d) The coefficients of first principal mode in vertical direction. Yellow star is the epicenter. Red and blue arrows are the positive and negative coefficients, respectively.

图 4为GPS连续站震后形变时间序列(篇幅有限，仅列出震后形变最大的4个站)，为清晰显示震后形变特征，已扣除时间序列中的速度、周期项、阶跃和同震形变.黑色实点为观测值，表示残差和震后形变之和，红色实线表示震后形变的模型值，RMS值显示各站震后形变与模型很好的拟合.此外，还可以看出震后形变初期衰减非常快，衰减速度逐渐变缓，震后形变还在持续.cmbl站是观测到的震后形变最大点，截至6月10日，NEU方向累计的震后形变分别达到107 mm, 135 mm和187 mm.

图 4

Fig. 4

图 4 GPS震后形变时间序列观测值和模型值 黑点表示震后形变观测值，红线为震后形变模型值，蓝线为地震时刻. Fig. 4 Observed and modeled postseismic deformation from GPS time series Black dots are observed postseismic deformation. Red lines are modeled postseismic deformation. Blue lines are time of earthquake occurred.

有几点需要说明，第一，新西兰地区历史地震多发，部分站点受到2004年11月22日M_W7.1(cncl, hoki, quar, waka, mtjo)，2009年7月15日M_W7.8(cncl, dund, dunt, hoki, kara, quar, waka, lytt, mqzg, mtjo, ous2)，2010年9月3日M_W7.0(mqzg, lytt)，2013年7月21日M_W6.5(cmbl, wgtn, wgtt), 2013年8月16日M_W6.5(cmbl)地震的影响.PCA分析M_W7.8 Kaikoura地震的震后形变时，需要预先扣除这4次地震的震后形变.第二，各站观测时间差异较大，为保留尽可能多有效信号又防止过度插值，PCA分析取所有站起始和终止时间的第75个百分位数作为起始和终止时间，时间窗为2013.3986—2017.4397，其中震后共209天.第三，PCA分析时间窗内所有站的震后形变信号，但各站独立估计参数时，分别对各站所有观测历元数据进行建模.

GPS观测的地表形变是各类形变的综合效应(孟国杰等，2016)，而震后松弛常常需要几年甚至几十年才能完全释放(Segall，2010)，在这期间，若未分离震后形变的影响而直接将观测的GPS速度替代震间速度，并以此研究稳定的板块运动是非常不严谨的.本文定量地研究震后形变对地表速度的影响，这里将问题简化，考虑单次地震的影响，周期性运动和噪声对速度影响较小，在此忽略，将式(1)简化为：

(3)

对其求导，那么t_i时刻地表GPS站的速度为：

(4)

同理，对数震后形变模型对应t_i时刻地表GPS站的速度为：

(5)

式(4)和(5)表明地表观测的GPS速度为稳定的震间速度和震后形变效应贡献部分的叠加.以cmbl为例，震后形变采用对数模型，从图 5可以看出，NEU方向的震间速度分别为28.83 mm·a^-1, -20.42 mm·a^-1, -7.35 mm·a^-1.地震发生时，快速的震后松弛响应导致地表形变速度NEU方向分别达到2958.17 mm·a^-1，3684.07 mm·a^-1，5108.12 mm·a^-1, 远高于震间速度，随着时间推移，站点速度逐渐收敛于稳定的震间速度.截至2017年6月10日NEU方向的速度仍为133.58 mm·a^-1，112.05 mm·a^-1，175.58 mm·a^-1.

图 5

Fig. 5

图 5 cmbl站的震间速度和地表观测速度对比 蓝线为震间速度，红线为地表速度. Fig. 5 Comparison between the interseismic velocity and the observed surface velocity of cmbl site Blue lines denote interseismic velocity. Red lines denote observed surface velocity.

(1) 本文通过蒙托卡罗方法合成模拟数据的过程有一些粗糙之处有待进一步完善，比如误差设置为高斯白噪声而忽略了有色噪声，以及同震和震后形变、衰减常数参数之间没有关联而完全随机采样.

(2) 迭代PCA方法将震后形变映射到时间域上共同的衰变特点和空间域上不同的响应，时间域的共同衰变过程表现为具有相同的衰减常数τ.衰减常数τ与深部物理性质和状态相关，从物理意义来看，基于余滑的对数模型的衰减常数τ与摩擦特性(a-b)和刚度κ有关，基于黏弹性松弛的指数模型的衰减常数τ与黏滞系数η和剪切模量μ有关，因此在一定范围内具有相同的衰减常数，而适用空间范围应结合地区的深部物理性质来判断.

(3) 迭代PCA方法中的PCA分析将信号映射到特征向量组成的空间，它根据信号强度来提取信号的主要成分，但不能区分信号成分或类别.因此，在PCA分析时程序会预先扣除模型里较强的信号，如同震形变、速度、周期项等，仅保留震后形变和残余信号；而个别站可能由于水文和环境变化等影响还存在较强的未模型化的信号，为避免干扰，PCA分析时需要排除这些站.信号的盲源问题仍是国内外研究的难点和热点，下一步研究将尝试其他盲源分离方法.

(4) 本文分析的新西兰地区37个GPS站分布较广，有助于把握整个新西兰地区Kaikoura地震震后形变的空间分布和时间衰减特点，但近震区站点稀疏，地质构造十分复杂，研究震区更细节的震后形变特点和震后形变机制则需要更多关键点的数据.

本文在现有参数估计方法的基础上扬长避短，提出了针对震后形变的迭代PCA参数估计方法，并利用模拟数据测试新方法在参数估计上的表现.参数误差的分析显示速度、同震和震后形变等参数与震后形变的衰减特征即衰减常数有关，速度与震后形变存在耦合.迭代PCA方法通过PCA分析区域内震后形变，提高震后形变信号的信噪比，从而估算出更准确的衰减常数，在此基础上融合了整体建模方法估算速度、同震形变和震后形变振幅等参数，并通过迭代降低速度和震后形变的耦合.因此，迭代PCA方法表现更稳定可靠，可以有效提高参数估计的精度，尤其当时间序列中的震后形变信噪比较低时，迭代PCA方法的优势更为明显.迭代PCA方法在速度、同震和震后形变参数估算精度的提高对GPS时间序列分析，以及地壳形变、同震破裂以及震后形变机制等相关研究具有重要意义.

此外，本文将迭代PCA方法应用到新西兰地区的GPS时间序列参数估计中，有效的分离Kaikoura地震的震后形变，为研究该地区的震后形变机制以及岩石圈流变学提供了更可靠的数据源.震后形变的时空分布显示：在时间域上，衰减常数为4天，震后形变衰减较快；空间域上，在南岛北东部和北岛最南部水平方向上震后形变较大；垂直方向上，南岛大部分地区和北岛最南部发生抬升，北岛北部大部分地区发生小幅下沉.cmbl站发生最大震后形变，且仍在持续松弛中，截至2017年6月10日震后松弛仅为同震形变的9%(cmbl站同震形变NEU方向分别为2403.54 mm，1406.32 mm和1031.84 mm).在此基础上，本文还推导了震后形变对地表速度的贡献理论公式，并以cmbl站为例定量分析了Kaikoura地震的震后形变对速度的贡献，可以看出震后形变对地表速度的影响非常大，因此，震间形变的研究必须剥离历史地震的震后效应影响.