地震的孕育是一个长期、复杂、缓慢的过程，慢滑移或无震蠕滑事件是伴随活动断层地震应力成核的重要部分。地震发生时，断层破裂会释放小部分应变能，而大部分能量在常规地震前后以无震蠕滑或慢滑移的形式进行释放^[1-3]。因此，检测蠕滑时空分布及其演变特征对于地震危险性评估和地震孕震机制研究至关重要。当前，基于GNSS观测网络已涌现出大量针对蠕滑形变信息的检测和反演方法，用于获取震间蠕滑和震后余滑的时空分布特征及演变过程^[4-7]，有力促进了地震震源机制和断层破裂过程的研究。然而，由于GNSS观测资料受海潮及固体潮、地下水、大气荷载、本地效应、共模误差和解算模型残留误差等复杂时空噪声的影响，难以从庞大的GNSS观测网络中快速检测蠕滑形变时空分布信息。通常情况下，GNSS时空滤波、地表形变信息提取和地下断层蠕滑分布反演单独进行。为进一步提高GNSS观测数据信躁比和微形变的检测能力，本文根据断层破裂所产生的地表形变时空相关性高的特点，提出集GNSS网络时空滤波、地表形变信息检测和地下断层蠕滑分布反演于一体的断层蠕滑时空演变分析模型，对覆盖断裂带地表的GNSS位移时空序列进行主成分提取，利用主成分信息快速检测并反演蠕滑断层的时空分布，该方法可显著提高反演信噪比和反演效率。以2005年苏门答腊M_W8.6地震震后余滑和2006年墨西哥Guerrero州慢地震事件为例，对断层的蠕滑过程及演变特征进行分析。

1 GNSS滤波与地表形变提取

对GNSS观测站的时间序列进行预处理，再通过拟合去除长期趋势项、年/半年周期项、阶跃等共模误差的影响；然后根据地壳形变的高空间相关性特点，对残余序列进行主成分时空响应分析^[8]。设存在一个由m个站构成的GNSS观测网络，共观测n个历元，组成时空矩阵为X _m×n，其中每一行代表给定站点N、E或U向所有历元的观测值，每一列代表给定历元在所有站点N、E、U分量的位移值。对矩阵实施中心化：

$ {\mathit{\boldsymbol{X}}^\prime }(i,j) = \mathit{\boldsymbol{X}}(i,j) - \frac{{\sum\limits_{k = 1}^m {{\mathit{\boldsymbol{X}}_0}} (i,k)}}{m} $

(1)

式中，i、j为矩阵行、列号。对中心化后的矩阵X′ _m×n进行奇异值分解，即求正交矩阵U _m×m、V _n×n和由奇异值构成的对角矩阵S _m×n，可表示为：

$ {\mathit{\boldsymbol{X}}^\prime}_{m \times n} = {\mathit{\boldsymbol{U}}_{m \times m}}{\mathit{\boldsymbol{S}}_{m \times n}}{\mathit{\boldsymbol{V}}_{n \times n}}^{\rm{T}} $

(2)

将奇异值按降序排列，前几个较大的奇异值代表位移时空序列的主模式分量。若取前r个奇异值为主模式，则原时空矩阵X可近似表示为：

$ \mathit{\boldsymbol{X}} \approx {\mathit{\boldsymbol{X}}_r} = {\mathit{\boldsymbol{U}}_r}{\mathit{\boldsymbol{S}}_r}{\mathit{\boldsymbol{V}}_r}^{\rm{T}} $

(3)

式中，S _r为前r个较大奇异值组成的对角阵；V _r为r个时间特征向量；U _r为r个空间特征向量。主模式时空响应通常用主模式分量中特征向量除以该特征向量中绝对值的最大值来表示，即对每个时空特征向量分别进行归一化处理，最大时空响应为100%，主模式时空响应可表现高时空的地壳形变信息。

模拟设置已滑动断层，并在断层上方地表覆盖区域按10 km等间隔布设63个测站，根据设定的断层蠕滑类型和滑动参数，基于位错模型^[9-10]，正演地表各站N、E、U向位移时间序列，并利用Fakenet模拟软件^[11]模拟不同信噪比的白噪声和有色噪声。合成较为真实的地表位移时间序列，再利用主成分分析法研究断层不同蠕滑特征与地表形变主成分之间的关系。大量实验结果表明，在信噪比不低于1时，通过地表形变时空分析均能反映实际的断层蠕滑特征，限于篇幅仅列出走滑实验结果(图 1)。设置长50 km、宽10 km、深10 km、走向90°、倾角70°的蠕滑断层片，模拟生成100 d呈指数函数变化的断层走滑时间序列。从提取的第一主成分(PC₁)时间和水平空间响应可以看出，时间响应曲线(图 1(c))同样呈指数形式的递增演变趋势，与事先设定的断层蠕滑特征(图 1(a))一致。N、U向(图 1(d))空间响应程度相对较小，约在±20%之间，且分布无规律。相比之下，E向空间响应分布呈明显规律性，量值主要集中在±80%~±100%。响应梯度高值区与断层位置一致，由此可揭示蠕滑断层的空间分布，断层两侧的空间响应方向相反，且与走向平行，推测断层蠕滑以走向方向为主。在N向空间响应中，断层以南的响应程度明显大于断层以北的响应程度，由此可判断断层倾向为S-N向。由此可见，通过地表GNSS位移主成分时空响应可判断断层的蠕滑特性，由时间响应可确定断层蠕滑的演变特征，由空间响应梯度、方向、大小可判段蠕滑断层的分布与类型^[12]。因此，根据地表GNSS观测网络提取断层形变主成分，应先确定蠕滑断层的发生位置、时间段和蠕滑特征，并将此作为先验信息，再对主成分信息进行断层蠕滑时空分布反演。

2 蠕滑时空反演模型构建

利用地表位移主成分空间分布U _r，反演断层滑移空间主模式分量D _i^[5]：

$ {\mathit{\boldsymbol{U}}_i} = G{\mathit{\boldsymbol{D}}_i}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} (i = 1\backsim r) $

(4)

式中，G为Okada位错模型格林函数。根据滑动断层几何参数(长、宽、深、方位角和倾角)的先验值，由式(4)可知，GNSS观测位移与断层位错呈线性函数关系。设反演所得断层滑移的整个时空分布为L，则整个GNSS位移时空序列X可表示为^[5]：

$ \mathit{\boldsymbol{X}} \approx \sum\limits_{i = 1}^r {{\mathit{\boldsymbol{U}}_i}{\mathit{\boldsymbol{S}}_i}{\mathit{\boldsymbol{V}}_i}^{\rm{T}} = } G\mathit{\boldsymbol{L}} $

(5)

将式(4)代入式(5)得：

$ G\mathit{\boldsymbol{L}} = \sum\limits_{i = 1}^r {G{\mathit{\boldsymbol{D}}_i}{\mathit{\boldsymbol{S}}_i}{\mathit{\boldsymbol{V}}_i}^{\rm{T}}} = G\sum\limits_{i = 1}^r {{\mathit{\boldsymbol{D}}_i}{\mathit{\boldsymbol{S}}_i}{\mathit{\boldsymbol{V}}_i}^{\rm{T}}} $

(6)

由式(6)得：

$ \mathit{\boldsymbol{L}} = \sum\limits_i^r {{\mathit{\boldsymbol{D}}_i}{\mathit{\boldsymbol{S}}_i}{\mathit{\boldsymbol{V}}_i}^{\rm{T}}} $

(7)

通常情况下，为获取断层滑移精细的空间分布，考虑到断层为非均匀滑动，将断层划分为规则的断层片格网，并解算每个断层片的位错量，这会导致解算参数的个数远大于观测方程数，因此附加不同约束条件对保持解的稳定性至关重要。本文考虑的约束条件为：

1)蠕滑空间分布的平滑约束。断层滑动在区域范围内通常为平缓过渡，以避免滑动分布解出现振荡。采用拉普拉斯算子进行平滑约束，使解算结果更加合理。

$ \begin{array}{l} \nabla = \frac{{s(i, j - 1) - 2s(i, j) + s(i, j + 1)}}{{{{(\Delta x)}^2}}} + \\ \ \ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{s(i - 1, j) - 2s(i, j) + s(i + 1, j)}}{{{{(\Delta y)}^2}}} \end{array} $

(8)

式中，s(i, j)为位于第i行、第j列子断层的滑动量，Δx、Δy分别为相邻子断层沿走向和倾向的距离。按式(8)求解所有断层片相应的二阶差分算子，使相邻断层片滑动梯度较小，即

$ d_n^\nabla = H_n^\nabla \cdot s + {\varepsilon _\nabla }, {\varepsilon _\nabla }{\rm{\backsim}}N(0, {\gamma ^{\rm{2}}}\mathit{\boldsymbol{I}}) $

(9)

式(9)等号左边为虚拟观测值，一般假定为0。$H_n^\nabla $为拉普拉斯二阶差分算子；s为各断层片位错量；γ为光滑因子，用来描述断层面上各断层片的平滑程度。

2) 根据前文分析的断层蠕滑特征的先验信息进行约束。设s₀为待估参数s的先验值，ε₊为参数先验误差，转化为观测方程的形式为：

$ {s_0} = \mathit{\boldsymbol{I}} \cdot s + {\varepsilon _ + }, {\varepsilon _ + }{\rm{\backsim}}N(0, \beta \mathit{\boldsymbol{I}}) $

(10)

式中，β用于约束断层滑动量的大小和方向，使其控制在合理范围内。

3 案例分析 3.1 2005年苏门答腊M_W8.6地震震后余滑

苏门答腊俯冲带位于印度洋板块、澳大利亚板块和其他次级板块的交汇部位，历史上苏门答腊地区地壳活动十分活跃，曾发生过多次大地震。2005-03-28印度尼西亚苏门答腊岛附近海域发生M_W8.6级强烈地震，震中位于2.2°N、97.0°E，破坏半径达544 km，有10个GPS站观测到地表数据。该区域板块构造及解算的震后333 d的位移时间序列^[5]如图 2所示。

对由10个测站构成的GNSS观测网络提取主成分，结果见图 3。由图可见，时间响应变化率在震后急速增加，约在第50 d达到最大值，之后开始衰减，并在约第250 d趋于稳定，显示了震后蠕滑随时间指数衰减的过程，也可反映大震之后地壳的粘弹性弛豫形变。空间响应程度最大的区域主要分布在以97°E、1°N为中心，约50 km为半径的区域内。由空间响应梯度可判断断层走向为SE-NW向，N向和E向空间响应量值相同，U向空间响应明显呈对称分布，表明震后蠕滑特征以逆冲为主。以此为反演先验信息，并参考相关研究^{[5, 13]}，确定反演所用断层蠕滑几何参数为：长320 km、宽150 km、深106 km、倾角20°、方位角310°。

将断层沿走向和倾向划分为40 km×30 km大小的断层片，反演每隔40 d的蠕滑时空分布，结果见图 4，图中红色箭头表示滑移方向，从图中可以清晰地看出震后断层滑移的整个时空演变过程。随着震后时间的推移，滑移量和滑移空间的影响范围逐渐增大，直至最后稳定。在时间分布上，约从震后第10 d开始出现明显的断层蠕滑现象，此时滑移量约为20 cm，之后以约0.5 cm/d的速度呈指数形式增大，直到第250 d才趋于稳定，此时达到最大滑移量约150 cm。由此可见，通过GNSS网络滤波反演方法能够快速准确地获取2005年苏门答腊M_W8.6地震震后蠕滑时空分布及演变过程。

3.2 2006墨西哥Guerrero州慢滑移

墨西哥位于拉丁美洲北部，在板块构造上为太平洋板块、美洲板块、加勒比板块和科科斯(COCOS)板块的交汇地区。墨西哥地区地壳构造运动极其复杂和强烈，地震活动频繁，震级较大，地震区分布较广。墨西哥西南部Guerrero州于2006年发生慢滑移事件，为目前世界上观测到的最大慢地震之一，覆盖断层的地表记录中有15个GPS站的观测数据。该区域板块构造背景及GPS站的单日解位移时间序列^[14]见图 5。

提取由15个GPS站构成的GNSS观测网络的主成分信息(图 6)。由图 6(a)可知，时间响应在2006-03开始发生明显的加速增加变化，约在2006-07变化率达到最大值，9月开始出现明显转折，逐步减小并趋于稳定。由图 6(b)可知，U向空间响应显著，且呈对称分布，可判断蠕滑断层走向为NWW-SE向，蠕滑特征以逆冲为主。确定反演所用蠕滑断层几何参数为：长300 km、宽150 km、深106 km、倾角10°、方位角292°^[14]。

图 7为不同时间点滑移时空分布反演结果，图中红色箭头表示蠕滑方向，可以看出，此次慢滑移事件从2006-02开始出现明显的滑移现象，此时滑移量约为10 mm；之后滑移量逐渐增加，空间影响范围也逐渐增大；到2006-09滑移速度开始减慢，并趋于稳定；至2007年初，滑移量达到最大值，约为200 mm。在空间分布上，最大滑移区域主要分布在以99.6°W、17.5°N为中心，半径约100 km的范围内。滑移特征以逆冲滑移为主，呈挤压形变特性，与相关研究结果^[14]吻合。

4 结语

断层破裂引起的地表形变具有高空间相关性特点，提取断层地表的整个GNSS观测网络的位移时空主成分能很好地反映区域地壳形变特征，有效减弱空间不相关噪声(如随机游走和白噪声)的影响。本文集GNSS时空数据处理、断层形变检测、断层蠕滑反演于一体，基于Okada位错模型，顾及拉普拉斯平滑约束和不等式约束，先由地表GNSS观测网络空间主成分反演蠕滑断层空间分布，再与地表GNSS观测网络时间主成分相结合，实现断层蠕滑时空分布的快速检测。以2005年苏门答腊M_W8.6地震震后余滑和2006年墨西哥Guerrero州慢滑移事件为例，该方法能很好地检测并反演震后余滑和慢地震等断层破裂过程和时空演变特征。随着GNSS观测网络的逐步加密和连续、定期观测数据的持续积累，区域GNSS地壳形变监测网可作为一个整体时空观测单元，鉴于地壳形变的高空间相关性，利用整个GNSS网络长时间观测序列主成分时空分析方法能准确得到背景噪声场，最大限度地消除或削弱观测数据中非构造形变信息和非相关噪声，从而提高断层蠕滑检测与反演信噪比及可靠性。