0 引言

强震孕育、发生与震后调整过程始终伴随着地壳形变发生，因此地壳形变是地震过程中最直接的现象之一(武艳强等，2016)。近30年来，以GNSS(Global Navigation Satellite System)为代表的空间对地观测技术迅猛发展，随着中国地壳运动观测网络和中国大陆构造环境监测网络(CMONOC)的实施，产出了丰富的GNSS观测结果(Wang et al，2001、2020；Liang et al，2013；李强等，2012；Zheng et al，2017)，为区域地壳变形与强震过程研究提供了高精度、高时空分辨率、高可靠性的数据结果。自1999年中国大陆开展大规模GNSS观测以来，西部地区发生了7次M≥7.0地震，分别为2001年昆仑山口西8.1级地震、2008年于田7.3级地震、2008年汶川8.0级地震、2010年玉树7.1级地震、2013年芦山7.0级地震、2014年于田7.3级地震和2017年九寨沟7.0级地震。上述GNSS观测资料和强震震例的积累为研究中国大陆西部地区地壳形变场与M≥7.0强震的孕育关系提供了可靠的数据约束和研究对象。

构造地震的孕育发生是地壳相对运动造成应力、应变在局部区域积累和集中，最终导致断层错动释放应变能的过程。在整个地震周期中，同震应变释放和震后调整过程均比较显著，且可被多种观测手段识别，相关研究极大地推进了对地震现象和物理过程的认识(Miller et al，1993；Sathyaseelan et al，2016；邵志刚等，2011；李强等，2012；Wang et al，2001；Shen et al，2009；Wu et al，2013；郑智江等，2017)。另一方面，从观测资料中提取大地震发生前的孕震晚期变形特征对于认识强震孕育过程、开展强震物理预测探索具有重要意义，是值得研究的重要科学问题。本文利用1999年以来CMONOC产出的GNSS流动观测资料，解算得到了中国大陆西部2组震前GNSS速度场，在误差分析的基础上构建最小二乘配置模型(江在森等，2010；Wu et al，2011)，通过GNSS变形场分析，研究了走滑型和倾滑型强震的孕震形变场共性特征，初步分析了强震前区域地壳变形动态特征与大地震孕育晚期的可能关联，探索GNSS形变场用于强震危险区预测的技术途径。

1 GNSS速度场结果

本文GNSS速度场结果包括2组，第一组为1999~2001期速度场结果，观测时间为1999年和2001年，数据解算由王敏等(2003)采用GAMIT/GLOBK软件和QOCA软件完成(Herring et al，2010a、2010b，Dong et al，1998)，具体解算策略可参阅相关文献(王敏等，2003；Wang et al，2020)。第二组为1999~2007期速度场结果，观测时间为1999年、2001年、2004年和2007年，采用GAMIT/GLOBK软件解算得到。为了削弱2001年昆仑山口西8.1级地震的影响，在进行GLOBK融合解算时，对地震显著影响区的测站施加了等速度约束，即约束地震前后的GNSS速度保持一致，该处理方法可以有效去除同震及震后3年的影响(Wu et al，2015)。为了直观反映变形信息、增强2组结果的可比性，基于1999~2001期和1999~2007期2组ITRF(International Terrestrial Reference Frame)参考框架下的GNSS速度场，采用Shen等(2005)和江在森等(2009)的方法，通过基准变换扣除整体运动(V=R×ω，V为刚性运动，R为椭球半径，ω为欧拉矢量)。图 1给出了区域参考基准的速度场，图中二级活动地块边界带来自文献(张培震等，2003)，其中图 1(a)中包括692个测站、图 1(b)包括708个测站。

图 1给出的GNSS速度场结果显示，青藏高原内部存在较大尺度的NS-NNE向挤压变形、青藏高原中部地区存在明显的近EW向拉张变形，并且沿羌塘和巴颜喀拉地块边界带存在明显的左旋剪切变形分布。青藏高原东北缘地区的变形特征以近NE向挤压变形为主，南北地震带中段的龙门山断裂带受青藏高原东向运动的加载效应，表现出较大空间尺度的挤压变形分布。川滇菱形地块的地壳运动速率明显快于与之相邻的羌塘地块、巴颜喀拉地块和华南地块，其东边界表现为显著的左旋剪切变形特征。新疆天山地区的变形方式以近NS向挤压变形为主，天山西段的变形量值明显高于中东段。具体而言，图 1(a)结果显示2001年昆仑山口西8.1级地震的震源区是中国大陆西部近EW向水平差异运动最显著的区域(陈兵等，2003)，在区域基准下36°N以北的地壳运动以西向为主、36° N以南的地壳运动以东向为主，表现为明显的EW向左旋剪切变形特征；图 1(b)结果显示巴颜喀拉地块东部自西向东表现为缓慢衰减的东向运动特征，量值由5~8mm/a逐渐衰减，至龙门山断裂带趋近于0，表明该区存在大尺度挤压应变积累现象(江在森等，2009)。

2 GNSS应变率计算方法与残差分析

相对于GNSS速度场结果，应变率能够直接给出变形量大小及空间差异信息，并且不同的应变矢量可以反映不同方向的变形分布。为了客观描述GNSS变形场，本文采用稳定性强、抗差性好的最小二乘配置球面应变率计算方法(武艳强等，2009；江在森等，2010；Wu et al，2011)。该方法首先通过统计方式获得GNSS观测速度信号的协方差函数，进而构建GNSS速度场连续分布函数模型，然后基于应变率与运动速率的偏导关系(Melchior，1978)解算出待定点的应变率，并据协方差传播定律获得不同应变率张量的误差分布，具体的计算流程可参阅文献(Wu et al，2011)。

球面最小二乘配置应变率计算方法的优势在于尽量减少GNSS应变率计算中的人为因素影响，其中对结果影响最大的信号协方差函数采用统计方法得到。GNSS数据(图 1)的协方差参数统计结果见图 2，由图可见，拟合函数为f(d)=Ae^-k2d2(其中，f(d)为协方差值，A、k为待定参数，d为测点间的距离)，显示GNSS速度观测信号的协方差统计值与函数有较高的拟合度。

由于应变率结果是否可靠取决于输入GNSS速度的拟合程度，对图 1中GNSS速度场的拟合残差进行统计分析，得到1999~2001期的拟合残差绝对值的均值为± 1.3mm/a、1999~2007期则为± 0.8mm/a，与GNSS速度的误差量值相当；同时，表现出后者的精度优于前者，其原因在于前者仅采用2期观测数据，而后者采用了4期观测数据。图 3给出的GNSS速度拟合残差的直方图结果显示，残差分布满足正态性和无偏性，具有随机特性。为了考察残差分布的空间特征，图 3(c)、3(d)给出了残差的空间分布，结果显示绝大多数测点的速度残差满足随机特性，少数测点的速度残差较大，例如图 3(d)中(103°E、23°N)处的测点。通过与图 1对比，发现其原因在于这些测点与周围测点的速度明显不协调，残差较大表明其模型值较小，即对周边测点速度的影响较小，该特点从另一方面表明最小二乘配置方法的噪声滤波功能对奇异测点的影响具有较好的抑制作用。

3 GNSS应变率场结果与M≥7.0强震孕育关系分析

由于地震的破裂错动发生在局部，震源区与邻近区域会存在明显的变形差异；同震释放的应变能是由震前积累起来的，因此在地震孕育过程中也会存在一定程度的时空变形差异，这种时空形成差异有可能被高精度的GNSS观测识别到。下面针对中国大陆西部的7次M≥7.0强震震前的地壳变形特征进行具体分析。

3.1 1999~2001期结果分析

本文重点研究的M≥7.0震例均发生在1999~2001期数据的观测区间之后，该时段的GNSS结果完全不受上述地震破裂的影响，主要反映了震间期构造变形和几次地震的震前变形信息。在对应变率场进行具体分析前，首先对其误差分布进行简要说明：EW向正应变率的平均误差为± 0.24×10^-8/a、NE向剪应变率的平均误差为± 0.16×10^-8/a、SN向正应变率的平均误差为± 0.21×10^-8/a。由于图 4(a)、4(b)给出的结果均基于上述应变率张量计算得到，误差量值不足图中相邻2条等值线的间隔，因此图 4中的变形特征可以被有效识别。GNSS主应变率结果显示，中国大陆西部地区的变形模式以NS-NNE向挤压变形、EW- NWW向拉张变形为主，变形基本呈连续分布特征，但存在明显的局部差异性。快速变形区主要分布在喜马拉雅地区、川滇地区、青藏高原中西部地区、青藏高原东北缘地区和天山西部地区。

考虑到7次地震中有4次强震(2001年昆仑山口西8.1级、2010年玉树7.1级、2014年于田7.3级和2017年九寨沟7.0级地震)属于走滑机制，为了研究上述走滑型地震前的区域变形特征，图 4(a)给出了1999~2001期最大剪应变率分布。结果显示，青藏高原的最大剪应变率整体偏高，其中喜马拉雅地区、青藏高原中西部地区、川滇地区、祁连地块中部等地区的剪切应变高值尤为突出。另外，新疆天山西段地区也处于最大剪应变率高值区。具体而言，4次走滑型强震均发生在最大剪应变率高值区边缘，其中，2001年昆仑山口西8.1级、2014年于田7.3级地震发生在青藏高原中西部地区的最大剪应变率高值区边缘，2010年玉树7.1级、2017年九寨沟7.0级地震发生在川滇地区的最大剪应变率高值区边缘。

7次地震中有3次强震(2008年于田7.3级、2008年汶川8.0级和2013年芦山7.0级地震)属于倾滑机制，其发震断层为近NS向和NE向。由于上述3个地震的动力成因主要为青藏高原东向运动，因此图 4(b)给出了EW向应变率分布。结果显示，青藏块体西部(92.5°E以西)以近EW向拉张变形为主，反映了青藏高原地壳物质东向运动而造成高原西部的EW向拉张特征。2008年于田7.3级张性地震发生在该EW向拉张区边缘。另一方面，青藏块体东部(92.5°E~100°E)以近EW向挤压变形为主，反映了青藏高原物质东向运动受到华南地块的阻挡，形成了空间跨度约900km的压缩区；汶川8.0级、芦山7.0级地震均发生在该大尺度挤压区的东边缘。

3.2 1999~2007期结果分析

上述基于1999~2001期观测给出的GNSS应变率结果显示，中国大陆西部4次走滑型M≥7.0地震与最大剪应变率高值区边缘、2次逆冲型强震与挤压应变高值区边缘、1次拉张型强震与拉张应变高值区边缘存在一定对应关系。然而，由于1999~2001期观测间隔较短、测站数量较少，致使应变率的精度和空间分辨率存在一定程度的不足。另外，2001年昆仑山口西8.1级地震的发生形成了长达426km的地表破裂，最大地表同震左旋水平位移达6.4m，最大垂直位移达4m(任金卫等，2005)，该事件必然会对研究区域的变形分布产生重要影响，特别是断层解耦导致的巴颜喀拉地块东向运动加速对区域变形场的影响深远。因此，基于1999~2007期数据(在数据解算阶段扣除了昆仑山口西8.1级地震的同震影响和震后3年的快速变形影响)，对观测结束后10年内的6次M≥7.0地震进行分析。

1999~2007期GNSS应变率误差分布如下：EW向正应变率的平均误差为±0.17×10^-8/a、NE向剪应变率的平均误差为± 0.12×10^-8/a、SN向正应变率的平均误差为± 0.17×10^-8/a，误差量值明显小于1999~2001期结果。由于图 5给出的结果均基于上述应变率张量计算得到，误差量值不足图中相邻2条等值线的间隔，因此图中的变形特征可以被有效识别。

图 5结果显示，2008~2018年的3次走滑型地震同样与最大剪应变率高值区边缘有较好的对应关系，同时3次倾滑型地震位于EW向应变率高值区边缘。对比图 4(a)和图 5(a)的结果显示，最大剪应变率动态调整最显著的区域为青藏高原地区，表现为青藏高原中西部最大剪应变率高值区范围的扩大和极值区量值的减小；同时，青藏高原北部地区的最大剪应变率有所减小，青藏高原东部地区的最大剪应变率量值有所增加。另外，青藏高原南部地区由于观测数据的变化导致的最大剪应变率差异不在本文的讨论范围。针对2007年之后的3次走滑型M≥7.0强震震中的最大剪应变量值，2010年玉树7.1级和2014年于田7.3级地震的震中位置最大剪应变率变化不显著，2017年九寨沟7.0级地震震中的最大剪应变率略有增加，由1999~2001期结果的(1.10±0.06)×10^-8/a增加到1999~2007期的(1.23±0.04)×10^-8/a。

对比图 4(b)、5(b)给出的EW向应变率，结果显示昆仑山口西8.1级地震对EW向应变率最显著的影响为，青藏高原东部地区EW向挤压变形范围和量值的增加，同时青藏高原西部地区EW向拉张区域存在收缩调整。针对2007年之后的3次倾滑型M≥7.0强震震中的EW向应变率量值，2008年于田7.3级地震震中的EW向应变率由1999~2001期结果的(0.98±0.30)×10^-8/a减小为1999~2007期的(0.34±0.21)×10^-8/a，2008年汶川8.0级地震震中的EW向应变率结果由(-1.54±0.14)×10^-8/a变为(-1.62±0.08)×10^-8/a，2013年芦山7.0级地震震中的EW向应变率由(-1.75±0.14)×10^-8/a变为(-2.24±0.09)×10^-8/a。

3.3 GNSS旋转率场的动态特征

邹镇宇(2015)经过公式推导显示，微元的旋转与剪应变不同，不会随着所取方向发生变化，因此GNSS旋转率在反映走滑断层变形特征方面具有独特优势。下面分别对2期GNSS资料的旋转率动态结果进行分析。

GNSS旋转率场结果(图 6)显示，中国西部可以划分为多个旋转变形单元，其中塔里木地块以顺时针旋转为主、柴达木地块和巴颜喀拉地块中西部以逆时针旋转为主、羌塘地块和拉萨地块以绕喜马拉雅东构造结顺时针旋转为特征、川滇菱形地块及东部地区以逆时针旋转为主。其中，最显著的为喜马拉雅东构造结附近的顺时针旋转区，反映了青藏东南部绕喜马拉雅东构造结的顺时针扭转变形特征，最大顺时针旋转率达到-4.5×10^-8rad/a。上述有规律的旋转变形单元，并非都以一级或二级活动地块边界为分界，相反多数单元均跨越地块边界带，其原因在于中国大陆西部的块体边界多以走滑运动(其中左旋走滑与逆时针旋转对应、右旋走滑与顺时针旋转对应)为主，而走滑断层两侧的差异运动刚好可以形成旋转变形单元(武艳强，2012)。

对比图 6(a)和图 6(b)结果，显示2001年昆仑山口西8.1级地震对柴达木地块和巴颜喀拉地块影响最为显著，同时对祁连地块也有一图 6结果还显示昆仑山口西地震对逆时针旋转率的增强大致截止于龙日坝断裂带，由于旋转量主要反映了走滑变形的影响范围，因此可以初步认为截至2007年底，昆仑山口西地震对东昆仑断裂带东段、龙门山断裂带和甘孜-玉树断裂的走滑运动加载影响有限。

4 讨论与结论

区域地壳变形过程与强震孕育发生密切相关，其中，震间期的断层闭锁、震前期的应力扰动、同震期断层错动、震后期调整过程均会对地壳变形场产生影响，因此地壳应变场中包含了强震孕育信息。本文研究的中国大陆西部7次M≥7.0地震与青藏高原的运动与变形密切相关，反映了印度板块NE向推挤青藏高原致使高原内部积累了大量的应变能并通过地震释放的过程。

(1) 本文采用的GNSS应变率计算方法具有稳定性高、抗差性强的特点(Wu et al，2011)，特别是通过统计方法获取的GNSS速度协方差函数可以有效降低应变率计算中的人为影响，保证计算结果的客观性。同时，对应变率场进行分析前，均进行了误差分析，可以保证识别到的震前变形信息有效。

(2) GNSS主应变率结果显示青藏高原主要表现为NS-NNE向挤压、EW-NWW向拉张变形特性。青藏高原中东部地区呈现由南到北、自西向东主压应变率方向逐步向东偏转的有序分布特征。青藏高原EW向应变率场结果表明，青藏块体西部(92.5°E以西)EW向以拉张为主、青藏块体东部(92.5°E~100°E)EW向以挤压为主，反映了青藏高原物质东向运动受到了华南地块的阻挡的动力学特征。

(3) GNSS旋转率场结果显示，中国大陆西部呈现交替旋转特性，在东构造结附近有一个显著的顺时针旋转区(极值为-4.5×10^-8rad/a)，向外过渡到逆时针旋转区。这一整体特征反映了青藏东南部大体上以东构造结顺时针扭转变形为主，以羌塘地块-川滇菱形地块向东-东南运动最快。柴达木地块表现为逆时针旋转变形特征(极值为1.0×10^-8rad/a)，塔里木地块表现为顺时针旋转变形特征(极值为-1.0×10^-8rad/a)，并且有由南向北交替衰减的趋势。

(4) 中国大陆西部7次M≥7.0强震均发生在与其构造背景一致的应变率高值区边缘，其中4次走滑型强震发生在剪应变率高值区边缘与走滑断层的交汇区，2次逆冲型强震发生在EW向挤压应变率高值区边缘与逆断层的交汇区，1次拉张型强震发生在EW向拉张应变高值区边缘与正断层的交汇区。关于上述现象，其背后的机制有2种解释，第一种解释是处于孕震晚期的发震断层存在强闭锁，该闭锁效应导致断层两侧变形分布在较宽范围、且近断层变形较平缓；而相邻断层段由于处于孕震中期或早期，变形主要分布在断层两侧较近的区域内、且近断层变形较剧烈；第二种解释是处于孕震晚期的断层由于距离上次强震的离逝时间较长，因而上一个强震事件的震后粘弹性加载效应较弱，而相邻断层段由于距离上次事件的离逝时间较短而粘弹性加载效应较强。总体而言，上述2种机制均可能导致强震更容易发生于应变快速积累高值区边缘，这种现象在意大利地区也被GNSS观测到(Riguzzi et al，2012)。

(5) 本文关注的7次M≥7.0地震均发生在巴颜喀拉地块及边界地区，实际上该强震过程开始于1997年玛尼7.5级地震，这些地震的孕育发生与巴颜喀拉地块的运动与变形密切相关，前人研究表明，玛尼地震对2001年昆仑山口西地震的触发作用不明显(陈兵等，2003；汪建军等，2009)。GNSS变形动态结果显示，昆仑山口西地震导致了青藏高原东部地区较大尺度EW向挤压应变增强现象，在一定程度上反映了巴颜喀拉地块加速东向运动引起的变形调整过程，以震源区东部挤压变形范围和量值的增大尤为明显，该效应有利于2008年汶川地震和2013年芦山地震的孕育发生。另外，GNSS旋转率场结果显示昆仑山口西地震的主要影响约达1200km范围内，表现为逆时针旋转的增强；该调整增强现象大致截止到龙日坝断裂带。因此，截至2007年的观测结果显示，昆仑山口西地震对东昆仑断裂带东段、龙门山断裂带和甘孜-玉树断裂的走滑运动加载影响有限。

综上所述，本文基于GNSS观测给出了中国大陆西部地区地壳变形场的定量结果，并在此基础上分析了研究区域7次M≥7.0地震震前的形变场特征，识别出了一些可能与强震孕育晚期相关的形变现象，这些现象可作为提取中长期强震危险区时地壳形变场异常判据的参考。从物理机制角度来看，上述现象在一定程度上反映了强震孕育晚期的断层强闭锁效应，据此通过系统分析中国大陆西部应变率场结果，结合断层构造变形背景，认为需关注安宁河-则木河-小江断裂带、阿尔金断裂带、东昆仑断裂带东段、嘉黎断裂、西秦岭北缘断裂、六盘山断裂带、天山断裂中段等的M≥7.0强震危险性。同时，由于震例数较为有限，本研究还有待后续更多震例的检验。

致谢: 本文的观测数据来源于“中国地壳运动观测网络”，1999~2001期GNSS速度场结果由王敏研究员提供，在文稿撰写过程中得到了中国地震局第一监测中心庞亚瑾副研究员的帮助，在此一并感谢。