近年来，大陆岩石圈应力状态、强度分布及层间耦合作用逐渐成为地球物理学界关注的热点^[1]。理解地幔物质流动(下地壳流动)在岩石圈底部(上地壳底部)拖曳作用如何传递应力，能增强对强震孕育机制的理解^[2]。上地壳发震层集中分布了大陆地壳震源，其应力传递和集中是地震发生的直接动力来源，发震层应力状态不仅受岩石圈底部(上地壳底部)拖曳力的控制，也受边界构造和物质分布等因素的影响^[3]。

青藏高原东南缘的西北部受高原物质溢出作用，西南部受印度板块持续推挤影响，东部则有强硬的四川盆地阻挡，3面力源汇聚造成该区复杂的构造背景和频繁的地震活动。本文利用EGM2008重力模型的11~36阶球谐系数计算研究区岩石圈底面地幔对流应力场，结合震源机制解数据，采用区域应力张量阻尼反演法得到芦山地震前和现今的发震层应力场，探讨不同区域的力学耦合情况与强震发生的关系。

1 青藏高原东南缘地幔对流应力场

Runcorn^[4-5]从牛顿层状粘滞流体力学出发，假设地幔为均匀牛顿粘滞体，对流局限于2个刚性的球形边界之中，其下边界重力效应可忽略，上边界处于静水平衡状态，由此提出利用卫星重力数据研究地幔对流应力场的基本理论。

1.1 计算方法

基于上述假设，地幔对流动力学方程(即Navier-Stokes方程)可简化为:

$ - \mu {\nabla ^2}\vec v = - \nabla P + \vec g\rho $

(1)

式中，μ为地幔粘滞系数，$\vec v$为地幔流体速度矢量，P为流体压强，$\vec g$为重力加速度矢量，ρ为地幔差异密度场，$\nabla $为梯度算符。式(1)左边还应有惯性力、加速度力及科里奥利力，但3者相对于重力小8~20个数量级，因此可忽略不计^[6]。

通过式(1)得到地幔对流对岩石圈底面产生的北向应力分量为^[7]:

$ {\sigma _N}\left( {\theta ,\lambda } \right) = \sum\limits_{n = 2}^\infty {\frac{{Mg}}{{4{\rm{ \mathsf{ π} }}{a^2}}}{{\left( {\frac{R}{a}} \right)}^{n + 1}}\frac{{2n + 1}}{{n + 1}} \cdot \frac{{{\rm{d}}{S_n}\left( {\theta ,\lambda } \right)}}{{{\rm{d}}\theta }}} $

(2)

东向应力分量为:

$ \begin{array}{l} {\sigma _E}(\theta , \lambda ) = \sum\limits_{n = 2}^\infty {\frac{{Mg}}{{4{\rm{ \mathsf{ π} }}{a^2}}}} {\left( {\frac{R}{a}} \right)^{n + 1}} \cdot \\ \frac{{2n + 1}}{{n + 1}} \cdot \frac{1}{{\sin \theta }} \cdot \frac{{{\rm{d}}{S_n}(\theta , \lambda )}}{{{\rm{d}}\lambda }} \end{array} $

(3)

其中，

$ \begin{array}{l} {S_n}\left( {\theta , \lambda } \right) = \sum\limits_{m = 0}^n {\bar P_n^m\left( {\cos \theta } \right) \cdot } \\ \left[ {\bar C_n^m\cos \left( {m\lambda } \right) + \bar S_n^m\sin \left( {m\lambda } \right)} \right] \end{array} $

(4)

式中，θ、λ分别为计算点的余纬和经度，M为地球质量，R为地球平均半径，a为地幔对流层上边界球面半径，g为计算点的正常重力值，${\bar C_n^m}$和${\bar S_n^m}$为完全归一化的n阶m级球谐函数系数，${\bar P_n^m\left( {\cos \theta } \right)}$为完全归一化的n阶m级缔合勒让德函数。使用WGS84提供的正常重力计算公式^[8]:

$ \begin{array}{l} g = 9.780\;326\;8\;\left[ {1 + 0.005\;302\;4{{\cos }^2}\left( \theta \right) - } \right.\\ \left. {\;\;\;\;\;\;\;0.000\;005\;8{{\sin }^2}\left( {2\theta } \right)} \right]\;\left( {{\rm{m}}/{{\rm{s}}^2}} \right) \end{array} $

(5)

完整的递推公式如下^[7]:

$ \begin{array}{l} \bar P_0^0\left( {\cos \theta } \right) = 1\;, \;\;\bar P_1^0\left( {\cos \theta } \right) = \sqrt 3 \cos \theta \;, \;\;\\ \;\;\;\;\;\;\;\;\;\;\;\;\bar P_1^1\left( {\cos \theta } \right) = \sqrt 3 \sin \theta \, \end{array} $

(6)

$ \begin{array}{l} \bar P_n^m\left( {\cos \theta } \right) = \sqrt {\frac{{(2n - 1)(2n + 1)}}{{(n - m)(n + m)}}} \cos \theta \bar P_{n - 1}^m\left( {\cos \theta } \right) - \\ \sqrt {\frac{{(2n + 1)(n + m - 1)(n - m - 1)}}{{(2n - 3)(n + m)(n - m)}}} \bar P_{n - 2}^m\left( {\cos \theta } \right) \end{array} $

(7)

$ \bar P_n^n\left( {\cos \theta } \right) = \sqrt {\frac{{2n + 1}}{{2n}}} \sin \theta \bar P_{n - 1}^{n - 1}\left( {\cos \theta } \right) $

(8)

实践证明，上述递推公式运算稳定、精度较高，其微商由式(6)~(8)容易推得^[9-10]。

计算应力分量时，还需选择阶数n的上下限，恰当划分低阶与高阶球谐系数，从而合理区分地幔对流尺度问题。综合分析认为，2~3、4~10、11~36阶球谐函数系数分别对应核幔边界、中下层地幔和上地幔(包括软流圈)的物质对流情况，本文研究的是青藏高原东南缘岩石圈底面软流圈内发生的地幔对流引起的拖曳力场，故阶数n选取11~36。

1.2 计算结果

使用EGM2008地球重力场模型的11~36阶完全归一化球谐系数，地球质量M取5.977×10²⁴kg，地球平均半径R取6 371.393 km，对流层上边界球面半径a取6 271.393 km。利用式(2)~(8)计算青藏高原及其邻区(20°~45°N，75°~115°E)岩石圈底面地幔对流应力场，网格间距为1°×1°，结果如图 1所示。图中箭头指示应力方向，长短代表应力大小。

Runcorn^[4-5]指出，高重力区的汇聚地幔流产生壳下挤压应力场，而低重力区对应上升的发散地幔流产生拉张应力场。由图 1可知，处于重力低值中心的青藏高原主体岩石圈底部存在发散地幔流，以32°N左右为界，分别向北和南产生拖曳应力; 而准噶尔、塔里木、柴达木等盆地处于汇聚地幔流区域，这些均与Runcorn提出的观点一致。应力大小则以105°E左右为界，西部较大，量级约为10⁶~10⁷ Pa; 东部明显减小，量级约为10³~10⁵ Pa。上述结果与黄培华等^[11]及熊熊等^[12]的结果大致相符，细微之处略有不同，与时变重力模型、阶数和计算方法的选取有关。

图 2是采用相同方法计算的青藏高原东南缘(20°~35°N，95°~107°E)岩石圈底部地幔对流引起的拖曳应力场，网格间距取0.5°×0.5°。

从应力分布来看，青藏高原东南缘主要受到3个方面的地幔对流冲击:西北部有来自青藏高原的发散地幔流，东部特别是四川盆地附近则有强烈的西向地幔流，这2股地幔流于98.5°E左右产生碰撞导致其各自转为南北向流动，南向地幔流与西南部来自印度板块推挤作用的北东向物质流在川滇菱形块体北端汇聚，构成青藏高原东南缘独特的地幔对流格局。

2 青藏高原东南缘发震层应力场

地震是岩层在区域应力场作用下，累积应变能达一定程度时突发破裂和错动而产生的。震源机制解是用地球物理学方法判别断层类型和发震机制的一种方法，它不仅能描述断层产状，还可以反映地震发生时断层的具体滑动情况。因此，利用震源机制解可以有效地反演应力场。

本文以芦山地震为时间分割点，将芦山地震发生前后的发震层应力场分别与地幔对流应力场进行相关性分析，观察芦山地震前与现今的应力状态及耦合性质的变化，提取强震的异常特征并对未来趋势作出分析。

2.1 研究方法

统计青藏高原东南缘1970年以来613次5级以上地震的震源深度发现，大多数震源位于10~25 km的发震层内，而使用大量震源机制解资料推断得到的区域应力场是发震层应力场的真实反映。

区域应力场阻尼反演法通过构建可调阻尼参数模型，引入平滑约束，最大限度地抑制邻区内的应力差异，并用最小二乘法求得稳定解^[13]。该方法反演速度快，能够较大程度上反映应力场的真实变化^[14]。Lund等^[15]将由该方法编制成的时空应力场反演(spacial and temporal stress invertion, SATSI)程序编写成可自动运算的MATLAB程序包，称之为MSATSI。

本文将收集到的震源机制解数据进行网格化处理，利用MSATSI程序运算，得到青藏高原东南缘部分网格的应力方位角及倾伏向，并进行薄板样条插值，最终得到网格间距为0.5°×0.5°的青藏高原东南缘发震层应力场。

2.2 数据资料

震源机制解数据主要来自美国哈佛大学CMT目录及前人文献，如表 1所示。编号1为CMT目录提供的1976~2015-05-01的强震数据，编号2时间截止到芦山地震前，编号3~6来自文献[16-19]。不同时间段的震源机制解可反映不同时期的应力状态，利用编号2~6共1 091个震源机制解反演得到芦山地震震前的发震层应力场; 利用编号1、3~6共1 131个震源机制解得到青藏高原东南缘现今的构造应力场。

2.3 反演结果

对研究区进行1°×1°的网格划分，并将每个网格点附近1°范围内的震源机制解归到该网格点上共同计算其应力张量。所有数据均选取节面Ⅰ，并将走向换算为与其垂直的倾向，使输入的文件格式符合软件要求。反演时软件自动选取最佳阻尼参数。

首先利用震前1 091个震源机制解反演得到图 3(a)所示的区域应力场，图中某些网格点因地震数目无法达到软件要求而缺失。对最大主压应力进行样条插值后得到与地幔对流应力点对点对应的0.5°×0.5°发震层应力场(图 3(b))。图 4(a)为加入了芦山地震及之后的震源机制解的反演结果，同样进行插值处理得到图 4(b)。

反演结果表明，青藏高原东南缘西北部受到来自青藏高原NE向的挤压作用; 东部则有华南块体阻挡产生的NW、NNW向应力，该阻挡作用导致来自高原的E向应力在菱形块体及其西侧发生顺时针旋转，菱形块体内部呈现北部SE向、南部SSE向应力分布情况; 西南部由于印度板块持续推挤产生NE向应力。这一结果与郭祥云等^[14]、曹颖等^[20]的研究结果相吻合。

比较图 3、4发现，芦山地震发生前后某些应力矢量存在细微偏转，特别是龙门山断裂带附近，地震发生前四川盆地一侧应力呈现NWW向，另一侧的松潘-甘孜区域则持续往E方向挤压，应变能恰在龙门山断裂处集中; 地震发生后，集中现象有所缓解，并向四川盆地内部推移。

3 地幔对流应力场与发震层应力场耦合关系分析

对于岩石圈型块体，其底部拖曳力在块体中传递较为顺利，对上层的发震层应力控制作用较强; 而地壳型块体多数存在壳内软弱带，应力传递过程中通过软弱带时易发生不可预知的偏转，使得上层应力脱离底部控制而独立存在。据此，可以研究岩石圈底部地幔对流应力方向与其上层发震层应力方向间的相关性，探讨青藏高原东南缘不同区域的力学耦合情况，分析其强震孕震机制和未来地震危险性。

点对点计算0.5°×0.5°地幔对流应力和发震层应力网格数据间相关系数，仅比较应力方位角，不考虑倾伏向，将对应网格点上2个应力方位角的差值进行归一化，即可得到应力方位角的相关性。若方向完全相同，相关系数为1;若完全相反，相关系数为0，由此绘制相关性图像。图 5采用芦山地震震前发震层应力场计算相关性，得到震前的耦合情况，其中图 5(a)为等值线图，图 5(b)为画有断层的彩虹图，不同颜色代表不同相关性:相关系数位于0.75~1之间的红黄色部分耦合程度高，为力学耦合区域; 0.45~0.75的绿色部分耦合程度居中; 0~0.45之间的蓝紫色部分则耦合程度低，为力学解耦区域。

由图 5可知，青藏高原东南缘大部分位于耦合与解耦的中间地带。耦合区域基本上按块体分布，东侧的华南块体耦合程度高，与青藏高原相连的西北部藏北块体和巴颜喀拉块体存在部分高度耦合带，西南部的滇西南块体则有小部分耦合。力学解耦主要发生在松潘-甘孜块体附近，毗邻耦合强烈的华南块体，中间为地震多发的龙门山断裂带。同时解耦的还有围绕川滇菱形块体的安宁河、则木河以及红河断裂带附近，东南部滇西南块体和拉萨块体之间的南汀河断裂带西南段也存在解耦现象。

图 6是现今应力耦合图。比较图 5、6可看出，红黄色区域逐渐减小、蓝紫色区域逐渐增大，这说明整体耦合程度有所降低，某些耦合区域在逐渐解耦。

从图 5、6中发现，强震大多发生在耦合程度变化较大的区域。这些区域通常是断裂带一侧相对另一侧耦合程度较高，而地震往往发生在耦合一侧。尤其在龙门山断裂带附近，完全解耦的松潘-甘孜块体与强烈耦合的华南块体碰撞，地震频发。此外，强震发生后耦合程度明显降低。

4 孕震机制及危险性分析

结合块体类型分析孕震原因，如图 7所示。若在地幔对流方向一致区域同时存在两类块体，岩石圈型块体在岩石圈范围内力学耦合，其壳内应力方向与地幔对流拖曳力方向较为一致; 地壳型块体下部存在壳内软弱带，软弱带下伏应力在带中传递时方向易发生大幅度变化，使岩石圈上层发震层应力方向与底部地幔对流应力方向有较大出入，耦合程度低。此时，两类应力耦合强弱交界区域中处于发震层深度的部分，应力不断集中，应变能积累，直至岩层破裂产生地震。

综上所述，青藏高原东南缘地幔对流应力场与发震层应力场耦合程度强弱交界带为地震危险性较高的区域。从现今耦合关系来看，龙门山断裂带仍处于耦合程度变化剧烈的区域，具备孕育强震的储能条件，地震危险性较高，极有可能再次发生地震。其他危险区域大致有:岷江断裂带附近、鲜水河-安宁河-则木河断裂带附近、红河断裂带及南汀河断裂带附近。

5 结语

位于青藏高原东南缘南部的滇中块体和滇西南块体覆盖了云南大部分区域，是研究区域内岩石圈最薄的地方，软流圈上涌构造明显。但与松潘-甘孜块体不同，软流圈上涌减薄了岩石圈，却没有使该区的地壳增厚，这也从侧面说明这一区域可能并没有发育较为成熟的壳内软弱带。水平方向上，GPS显示的上地壳形变仍然延续了北侧的挤压作用，围绕东构造结的旋转使得方向变为NS或NNE向，快波表现出来的却是近EW向的形变。由于岩石圈减薄，南侧区域的SKS剪切波更多反映的是软流圈流动方向，这种EW向流动可能拖动着岩石圈下部一起运动，从而上下层解耦。但该区剖面的速度结构没有显示大面积的解耦构造(即软弱带分布)，可能是软流圈上涌造成的EW向拖曳并不强烈，而块体的侧向挤压和断裂带的相互作用占据了主导地位，地壳内还没有形成解耦构造，但未来解耦的可能性较大。岩石圈形变依旧耦合，但与软流圈完全解耦。