基于震源应力与GPS应变分析喜马拉雅东构造结及其邻区地壳变形特征*

张 晨，季灵运，朱良玉，徐 晶

(中国地震局第二监测中心，陕西 西安 710054)

0 引言

印度板块与欧亚板块之间的相互挤压碰撞形成了青藏高原，这是地质历史上最重要的造山事件。喜马拉雅东构造结地区属于该构造碰撞的前缘，位于青藏高原东南缘、印度大陆东北以及缅甸板块西北角。该地区地貌复杂，构造活动剧烈，且地震频发。1900年以来，该地区发生5.5级以上地震近百次，且大多数集中于区域内部的主要地震带上。因此，该地区被称为地震研究天然实验场(杨建亚等，2017；王凯悦等，2021)。

地壳应力是地震孕育和发生的直接驱动力(Hardebeck，Michael，2006)，地壳应力场的方向也是评估断层潜在滑移的一个重要参数，地壳应力的时空变化能够较好地反映地震与应力之间的关系特征，因而能在地震危险性和风险评估中发挥至关重要的作用(Heidbach，Ben-Avraham，2007；Martinez-Garzon，2014；Walsh，Zoback，2016；吴微微，2020)。现有推算构造应力场的方法主要包括：测井扩张破裂、实地应力测量中的水压破裂和应力解除(Zoback，1992；Reynolds，2006)、震源机制解(徐纪人等，2008；王晓山等，2015；黄骥超等，2016；崔华伟等，2017)、断层滑移测量(Angelier，1979)等。其中，利用震源机制解求解构造应力场是研究应力构造的主要方法(谢富仁等，1993；崔效锋，谢富仁，1999；Hardebeck，Michael，2006；Wan，2010；Martinez-Garzon，2014；杨帆等，2019；王晓山等，2020)。在所有利用震源机制解反演应力场的方法中，阻尼应力张量反演(SATSI)是较为成熟的计算方法(Hardebeck，Michael，2006；Zhao，2013；Luo，2015；王晓山等，2015)，该方法被广泛应用于青藏高原地区的应力场研究(徐纪人等，2008；侯强等，2018；李君等，2019)。

随着现代大地测量学的发展，GPS大地测量技术已成为地震监测以及形变测量的主要手段(王敏，沈正康，2020)，国内外学者利用GPS数据获取了许多构造活跃区域的应变率场(Allmendinger，2007)，尤其是喜马拉雅造山带(Bilham，1997；Vernant，2015)、青藏高原(张培震等，2002；张清志等，2005)以及阿尔金断裂带(郑文俊等，2009)。前人利用GPS观测手段和震源机制解，获取全球范围的应力方向分布图(Zoback，1992)，该结果为研究全球地壳应力应变特征提供了指导作用，但针对喜马拉雅东构造结及其周边地区构造应力场的精细分析尚显不足。

因此，本文利用最新发表的GPS速度场数据，计算喜马拉雅东构造结及其邻区的地壳应变率场，然后收集该地区≥2.0地震的震源机制解，利用MSATAI程序包反演研究区的应力场，最后对比GPS主应变方向与震源应力场方向，分析研究区地壳深浅部构造变形特征。

1 GPS数据与应变率场

为了获取研究区的现今GPS应变率场，本文使用Wang和Shen(2020)收集的1991—2016年的长期GPS观测数据(图1)，参考框架为欧亚大陆，采用Shen等(2015)提出的一种通过速度场计算应变的方法。该方法较为稳健，不依赖于对数据的某种假设，且平衡了稳定性与分辨率两方面，计算结果更为全面可靠。计算结果如图2所示。

图1 喜马拉雅东构造结及其邻区速度场

由图2a可知，最大主应变分布自西向东、由北至南呈现顺时针偏转。区域内部的大型走滑断裂，其最大主压应变也较高，表明其具有较为强烈的构造变形活动。在96°E以西、30°N以北的青藏高原地区，应变以拉张为主，方向由近EW向逐渐转为NW向；而30°N以南的喜马拉雅区域，主要以SN-NE向挤压为主，主压应变方向近似与板块边界相垂直，表明该地区受到来自印度板块的NNE挤压作用。在96°E以东地区，主要包括地震活动最剧烈的青藏高原东南缘区域——川滇菱形块体，该区域受华南地块的阻挡，主应变复杂多变，整体表现为绕喜马拉雅构造结呈SN-NNW-NW向的顺时针偏转。巴颜喀拉块体内部主压应变方向为近EW向，川滇菱形块体东边界主压应变方向从北向南逐渐由近EW向转为近NW向。最大剪应变率高值区主要位于川滇菱形块体东边界，甘孜—玉树断裂带、鲜水河断裂带、安宁河断裂带、则木河断裂带、小江断裂带及其周边次级断裂(图2b)。本文研究区应变率场结果与前人研究得出的空间分布特性一致。

图2 喜马拉雅东构造结及其邻区面应变率(a)和最大剪应变率(b)分布

2 震源机制解数据与应力反演

本文使用的数据为1976年1月1日至2020年9月25日喜马拉雅东构造结及周边地区(23°～34°N，90°～105°E)的7 092个震源机制解，对应震级范围为2～7.9级。其中，420个震源机制解数据来源于GCMT(Global Centroid-Moment-Tenso)，对应震级范围为4.7～7.9级，其余数据来源于四川地震台，震级范围为2～6.2级。按照Aki和Richards(1980)的分类方法对震源机制解数据进行分类(图3)，分别是：正断层，逆断层，走滑断层以及斜滑断层，具体分类依据见表1。

图3 喜马拉雅东构造结及其邻区震源机制解分布

表1 本文研究震源机制解划分依据

通常情况下，由于研究区空间应力场分布的非均匀性，我们常将其划分成更为细小的多个子区域，并独立拟合出每个子区域的应力张量，但这种模型下的应力场空间变化存在难以解释的情况。因此本文采用阻尼应力张量反演法尽可能消除应力变化产生的虚假信息，该方法能使相邻子区域的应力张量分量之间的差异最小化，也令震源机制解与应力张量的匹配程度达到最高。鉴于该方法在空间应力场分布反演上具有极大潜力，Material-Garzon等(2014)在SATSI算法的基础上开发了MATLAB环境中运行计算的MSATAI程序包。在进行应力张量反演时，将震源机制解数据按照0.5°×0.5°的网格进行分区，共计660个子区域，如图4a所示。在输入参数设置时，为了保证更多的子区域能够反演得到应力结果，将每个网格地震最少个数设置为1，但由于每个子区域至少需要包含4个地震数据才能计算，此时，程序会根据该子区域周围的震源机制解对该区域的应力张量反演进行约束。最后，对于该区域，本文选取最佳阻尼系数为1.2(图4b)。

图4 0.5°×0.5°网格地震数目统计(a)及模型拟合曲线图(b)

本文的反演结果提供了3个主应力轴的方向以及反映3个主应力相对大小的参数：应力形因子(Pesicek，2010)，可表示为：

(1)

式中：、和分别是最大主压应力、中间主压应力以及最小主压应力；主要反映中间应力更接近于最大主压应力或是最小主压应力，其取值为0～1。当<05时，求得的最大主张应力轴(最小主压应力轴)方向可靠性较最大主压应力轴的方向相对较高，此时，中间应力呈现为压应力，与的大小表现为负相关，即越小，中间应力轴压应力越明显；当>05时，最大主压应力轴的方向可靠性更大，中间应力呈现张应力，与的大小表现为正相关。极端情况下，=1时，最大主压轴确定、最大主张轴不确定；=05时，两轴均确定；=0时，最大主张轴确定、最大主压轴不确定。

利用上述震源机制解数据，通过MSATAI程序包反演得到研究区的应力结果，如图5a所示，蓝色线段代表最大主压应力轴，红色线段代表最小主压应力轴。线段的长度表示应力轴倾角的大小，线段越短倾角越大，越长则倾角越小。

图5 震源应力场反演结果(a)、应力形因子R与最大水平应力SHmax分布结果(b)

为了更加清晰直观地展示该区域的应力状态，本文基于Lund和 Townend(2007)提出的最大水平应力计算方法，获取了该区域的最大水平应力分布情况，如图5b所示，图中红色直线表示最大水平应力max，最大水平应力的背景值为应力形因子。

3 研究结果

3.1 应力场反演结果

从图5a可以看出，在99°E以西地区，为NE向，呈NWW或NW走向；在羌塘块体内部，垂直于断裂走向，则与断裂走向平行；自西向东、由北向南存在NNE-NE向偏转，轴倾角也逐渐增大，其中怒江、澜沧江断裂处轴倾角近垂直(Qiu，Qiao，2017)；自西向东呈现出NWW-NW向偏转，且该区域的轴倾角较小，接近于水平。在缅甸弧俯冲带区域指向NNE向，指向SEE向。

在99°E以东地区，由北至南呈现EW-SEE-NS的偏转，轴倾角也存在较大的变化。巴颜喀拉块体内部、川滇菱形块体东边界轴倾角较小，其西边界轴倾角较大；表现为NS-NNE-EW的顺时针偏转，轴倾角较小。该区域主要是地震活动最强烈的川滇菱形块体，应力方向在区域内存在很大差异，川滇菱形块体东边界主要为主压应力方向是SEE向的走滑断层，其西边界则主要是主压方向为SSE向的正断结构，其内部西北区域为正断型应力状态，主要集中于金沙江断裂带以及理塘断裂，龙门山断裂则主要为逆断型应力状态。

从图5b可以看出，区域内最大水平应力max方向呈现出较为明显的分区特征。96°E以西地区，max主要为NNE向，且在羌塘地块内部<0.5，表明该区域处于双轴压缩状态；96°～99°E范围内，max方向呈现出NNE-NE-WE的偏转，该范围属于应力轴方向偏转的过渡区域，值大部分都大于0.5，表明与双轴拉张；99°E以东地区，max表现为围绕喜马拉雅构造结的顺时针偏转特征。本文所获应力场分布与已有研究结果基本一致(徐志刚，2017；Zhao，2013)。

本文所反演的应力场结果与该区域内的地震类型较为吻合，可见，使用该方法能较为准确地反演出区域震源应力场状态，这为我们分析地壳深部构造应力方向提供了一种可靠约束。

3.2 GPS主应变率与最大水平应力SHmax

GPS观测可提供地壳浅部的应变结果，而震源机制解所反演的应力场结果则从地壳深部揭示了地球动力学过程(Pan，2020)。将两者相结合分析，有助于认识研究区的应力状态以及构造活动。因此，将GPS主压应变率方向与利用震源机制解所获取的最大水平应力max进行比较分析，如图6所示，从而更好地认识研究区的地壳构造变形特征。由图6可知，整体上，最大水平应力max方向与GPS主压应变的方向几乎一致，呈现出顺时针偏转的现象。为了精细分析图6反映的各区域主要特征，将研究区分为A、B、C、D、E、F 六个区域。其中，A区主要包括羌塘地块的东南部分，B区主要是印度板块东北部的阿萨姆地区，E区则主要是滇西北区域，F区域主要是滇南块体范围，C区域范围最小，位于羌塘块体中部。A、E、F三个区域将川滇菱形块体的西边界包含于其中，D区域包含川滇菱形块体的东边界以及块体内部。

图6 GPS主压应变率方向与最大水平应力SHmax方向对比

按照上述分区，本文统计了该区域两种方向的夹角，如图7所示。其中，角度小于10°的占40%，10°～20°占比25.1%；夹角均值在A区为39.33°，在B区为37.35°，在C区为38.8°，在E、F区分别为42.63°、28.71°。GPS主压应变方向与最大水平应力轴方向分布方向一致的D区夹角均值为5.3°。

图7 GPS主压应变方向与最大水平应力SHmax方向夹角统计(a)、分区夹角均值(b)

由上述应力场反演与统计结果可知，在川滇菱形块体东边界以及块体内部，最大水平应力max优势方向与GPS主压应变几乎相同，且由北向南呈现NS-NW的偏转，表明该处地表浅部与上地壳深部受到了相同机制的驱动力；青藏高原向东南部挤出的作用力在川滇菱形块体东边界一直延续至小江断裂。而在川滇菱形块体西边界，最大水平应力max优势方向与GPS主压应变展现出不一致性，位于滇西南块体的红河断裂以及澜沧江断裂附近区域也存在此现象，说明在上地壳深、浅部可能具有不相同的驱动机制。

4 讨论

本文所获取的GPS主压应变方向与max方向分别从地壳浅部与深部提供了认识该区域构造活动的两种方式。但应力分布可能受到局部构造环境或者其它外部因素影响，因此在对结果的分析中需要综合考虑其它因素。

4.1 川滇菱形块体西边界SHmax与GPS主压应变方向差异分析

杨晓松和马瑾(2003)研究表明，青藏高原及其周边地区表现为地壳内部解耦现象，且在其东南缘该现象更为明显。印度板块向欧亚大陆的俯冲作用带来NNE向的压力，该作用力与青藏高原东南向的逃逸挤出相互作用(图8)，可能是造成区域A、B、E区差异较大的主要原因。

图8 喜马拉雅东构造结及其邻区地壳变形模式的示意图

常利军等(2015)指出，在青藏高原东构造结内部及其周边地区，测量的快波方向分别与最大左旋剪切方向和最大拉伸方向一致；通过GPS应变计算所得的地幔各向异性与测量的快波方向的一致性进一步说明该地区岩石圈变形属于垂直连贯变形模式。然而，本文反演结果中A区GPS主压应变与max方向存在较大差异，产生这种差异的原因可能有2个：一是震源机制解获取的地壳应力场是地壳弹性层以上的应力方向，而常利军等(2015)获取的快波方向是上地幔与地壳累积效应，观测应力方向的深部不一致；二是震源机制解的应力场方向是近50 a的结果，而快波方向通常被认为是岩石晶格长达百万年尺度的构造作用下产生的方向，两者在时间尺度上存在较大差异。因此，本文结果与通过快波方向得到的结果可能并不是同一种物理量。此外，该分区的震源机制解数据较少，也可能对结果有一定影响。

4.2 滇西南块体SHmax与GPS主压应变差异分析

F区主要包含红河断裂以及澜沧江断裂两条主要断裂带。GPS观测显示该区地壳主要围绕喜马拉雅东构造结旋转，方向为近SN向或NEE向；地球物理等多方面观测资料显示，缅甸微板块下方呈现高度倾斜的俯冲状态，层析成像结果显示在印缅山脉下方有一个高速异常体向东俯冲(Pesicek，2010)，而同时缅甸微板块整体随印度板块向北移动。快波结果表现出近EW向的形变，该处软流圈上涌，致使岩石圈减薄，SKS剪切波更多反映的是软流圈的方向，拖拽岩石圈下部一起运动，导致上下解耦，但速度场剖面结果显示该区域没有成熟的壳内软弱带，可能的原因是岩石圈的下沉引起如软流层拖拽力不强烈。这些使得地壳内耦合，岩石圈与软流圈解耦(侯强等，2018)，可能是引起F区中GPS主压应变与max方向分布杂乱的原因。

4.3 喜马拉雅东构造结SHmax与GPS主压应变差异分析

地震各向异性特性可以反映区域应力场、变形特征以及断裂活动特征，有助于研究地壳深部构造以及动力学机制。王凯悦等(2021)研究表明，在活动断裂或其附近快波偏振方向与断裂走向相一致。在拉萨地块东部地区，郭铁龙和高原(2020)关于剪切波分裂研究结果表明快波偏振方向为NS或NNE，这与本文B区右上角拉萨地块内部水平最大主压应力方向吻合；在察隅断裂南部，近震直达S波分裂结果表明，快波偏振优势方向与断裂走向一致，为NW向(黄臣宇等，2021)，与本文B区左下角最大水平应力方向一致；在印度板块东北角的阿萨姆地区，应力反演结果与印度板块向欧亚板块NE或NNE向挤压作用方向一致；而该区域GPS观测数据较少，这可能是导致B区max与GPS主压应变方向存在差异的原因。

5 结论

本文主要利用喜马拉雅东构造结及其邻区≥2.0震源机制解计算得到max方向与GPS主压应变方向，并对二者进行比较，以研究该区的变形模式。研究结果表明：

(1)在川滇菱形块体东边界以及块体内部，最大水平应力max方向与GPS主压应变方向基本一致，且由北向南呈现NS-NW向偏转，方向变化具有一定的连续性，说明该处地表浅部与上地壳深部受到了相同机制的驱动力；青藏高原向南部挤出的作用力在川滇菱形块体东边界从鲜水河断裂、安宁河断裂、则木河断裂一直延续至小江断裂。

(2)在川滇菱形块体西边界，max方向与GPS主压应变呈现出了显著的差异性，可能由于其构造环境复杂，深、浅部存在差异运动，在上地壳深、浅部具有不相同的驱动机制；在红河断裂附近区域，max方向与GPS主压应变方向也存在差异，可能是由于缅甸板块的俯冲作用力与岩石圈的下沉引起的软流层拖拽力等的共同结果。

感谢四川地震台龙锋高级工程师为本文研究提供的震源机制解数据。