引言

四川和云南地区是中国大陆强震多发地区．近年来，在中国地震局的年度地震趋势会商会上，有关专家多次将川滇菱形地块及其部分邻区划定为地震重点危险区(闻学泽，2000)．鲜水河—安宁河—则木河—小江断裂带作为川滇活动地块的东边界((张培震等，2003)，是一条大型的左旋走滑活动断裂带，其强震活动最为频繁(图 1)，在过去的300多年中曾发生过多次强震活动，其中多数发生在鲜水河断裂上，而安宁河断裂、则木河断裂和小江断裂自1833年后均未发生过M≥7强震．因此，针对川滇活动地块东边界的中长期大地震潜势的研究，对地震预测和防震减灾都具有非常重大的意义．

图  1  川滇地区主要活动断裂和地震分布

Figure  1.  Distribution of main active faults and earthquakes in the Sichuan-Yunnan block

下载: 全尺寸图片 幻灯片

闻学泽等(2008)通过小震静定位的时空分布特征，结合已有的地震活动性参数和形变测量分析结果，探讨了川西安宁河断裂、则木河断裂地区最近30年来地震平静的含义，并估计了安宁河断裂的两个闭锁段潜在地震的最大可能震级．易桂喜等(2004)和朱艾斓等(2009)基于地震重定位结果，沿鲜水河—安宁河—则木河—小江断裂带进行b值和局部复发时间的空间扫描，对该断裂带上现今应力积累的闭锁段位置进行了推断．王辉等(2012)利用40多年的小震资料分析了川滇地区1981年以来的19次M≥6.5地震与b值的时空分布关系，探讨了该地区强震活动与b值变化趋势之间的关系．

上述研究均为通过地震精确定位、地质调查以及地震活动参数来估计地震危险性．目前在对川滇地区应力应变状态的数值模拟方面也开展了很多相关研究. 例如，李玉江等(2009)模拟分析了云南地区构造应力应变场的年变特征，王辉等(2010)分析了断层分布和几何形态对川西及其邻区应变分配的影响，陈连旺等(2011)讨论了流变特性对青藏高原构造变形的影响．由于受模拟计算时间步长限制或巨大重力作用可能造成有限元网格强烈坍塌变形等影响，多数研究仅简单假设地壳在长期重力作用下岩体已基本处于重力均衡状态，计算时只考虑静岩压力之外的构造应力而忽略了重力作用(张东宁等，2007； 陈连旺等，2008； 王辉等，2008； 张怀等，2009)．但是青藏高原高海拔地形蓄积的重力势能，是高原物质侧向挤出的重要动力学因素，同时也是控制川滇地区地壳水平向偏应力(构造应力)和水平向偏应变的重要地球动力学因素之一(张健，石耀霖，2002； Yang，Liu，2009)．因此，充分考虑青藏高原在重力作用影响下的地球动力学特征是模拟地壳和上地幔构造应力状态的重要前提(祝爱玉，张东宁，2014)．

鉴于此，本文采用三维有限元数值模拟方法，以最新的川滇地区深部结构、岩石圈流变特性、构造应力场以及GPS观测给出的地块运动为约束，建立了包括川滇地块东边界中南段(即安宁河—则木河—小江断裂带)的三维有限元数值模型，研究在速度边界条件和重力势能等动力因素共同作用下该断裂带的应力分布特性，探索该断裂带现今应力状态及其与地震活动的关系，同时给出强震活动的主要控制因素以及未来可能的强震危险区．

1.   三维黏弹性模型

1.1   三维有限元模型的建立

利用已有的探测结果以及最新的南北地震带地区地壳上地幔三维高精度速度结构模型和物性参数，根据地质、地应力、地形变资料确定的区域内及周边构造地块间的相互作用形式，本文建立了包括安宁河—则木河—小江断裂带及其邻区的三维精细黏弹性有限元实体模型(图 2)，其范围为101°E—104°E、24°N—30°N．该模型考虑了研究范围内的主要活动断裂： 小江断裂、则木河断裂、安宁河断裂、大凉山断裂和丽江—小金河断裂(邓起东等，2002)．由于本文重点研究安宁河—则木河—小江断裂带上的应力状态，所以在该断裂带上采用多种尺寸的摩擦接触单元来反映断层，而在丽江—小金河断裂和大凉山断裂则采用一定宽度的特定材料参数单元模拟活动断裂(软弱单元)．与软弱单元相比，摩擦接触单元所模拟的活动断裂更能接近真实，但其收敛比较困难．模型中包括了安宁河—则木河—小江断裂带东西两边的华南地块和川滇地块．

图  2  川滇地块东边界中南段断裂带地区的三维有限元模型

Figure  2.  The three-dimensional finite element model of the Anninghe- Zemuhe-Xiaojiang fault zone

下载: 全尺寸图片 幻灯片

如图 2所示，该模型深度为200 km，分为地表、上地壳、中地壳、下地壳、地幔顶层(岩石圈深度)和上地幔(软流圈部分)等6个分层．采用美国地球物理中心发布的ETOPO1地形高程数据，其分辨率为1 s(National Oceanic and Atmospheric Administration，2014); 上地壳、中地壳和下地壳分层数据引自Stolk等(2013)最新结果; 岩石圈深度结果安美建和石耀霖(2007)的结果．用六面体单元对本文模型进行网格划分，横向和地壳深度的单元网格大小均为5 km，地幔深度的单元网格尺寸为10 km，共划分约30万个单元．

1.2   材料参数

人工地震测深和天然地震观测研究表明，川滇地区的地壳厚度变化剧烈，速度结构横向不均匀性较强．本文参考吴建平等(2006)和王椿镛等(2008)的波速结果，得到了研究区的弹性模量，如表 1所示．其具体计算公式为

式中，ν为泊松比，v_P和v_S分别为地壳某一层的P波和S波速度，E为弹性模量.

表  1  本文有限元模型的材料参数

Table  1.  Material parameters of the finite element model built up in this paper

地层	华南地块			川滇地块
	弹性模量 /（109Pa）	泊松比	黏滞系数 /(1021Pa	弹性模量 /（109 Pa）	泊松比	黏滞系数/(1021Pa
地表	90	0.25	2	80	0.25	1
上地壳	90	0.25	2	80	0.25	1
中地壳	100	0.30	500	90	0.30	100
下地壳	120	0.35	10	100	0.35	10
上地幔顶部	100	0.35	10	100	0.35	10
断层	10	0.49	1	10	0.49	1
软流圈	100	0.35	0.1	100	0.35	0.1
注：各地层的密度均为2 700 kg/m3.

下载: 导出CSV 

| 显示表格

利用软弱夹层断裂材料表示大凉山断裂、丽江—小金河断裂．软弱夹层断裂材料的弹性模量设置得明显低于周围地块介质，其泊松比大于周围地块介质．该泊松比的设定，除通过地震波进行计算外，还参考了以往的数值模拟工作，最终设置的泊松比较地震波计算结果偏大，主要是为了体现长期地质演化过程中岩石圈介质更倾向于流变性．根据该区域相对偏高的平均大地热流值，可推测川滇菱形地块的下地壳介质强度较低、相对较软，而菱形地块东边的华南地块和四川盆地的地壳介质则相对较硬(石耀霖，曹建玲，2008)．有限元模型介质的黏滞系数也列于表 1．

1.3   初始条件和边界条件

初始条件是地球动力学数值模拟中最为困难的问题．就现有的观测手段而言，初始地应力场不能直接获得，只能借助于反演等手段来得到，但是这种地应力状态不具有唯一性，因此本文在模型中加入重力作用. 由于岩石圈的长期流变性质，其应力状态趋向于“静岩压力”应力状态，即重力均衡状态，以此作为初始条件，然后施加边界条件．本文对青藏高原东南缘现今地壳形变的GPS观测结果(王敏等，2008)进行插值，给出了研究区各段边界的速度值．模型的水平位移约束条件为该速度值与计算时间步长(100年)的乘积，垂直方向位移保持自由; 模型的上表面为自由表面．为简化模型，将模型底部所有节点垂直方向(z方向)上的位移约束为0，但可以在水平方向上运动．

2.   模型检验

本文在考虑重力势能作用的条件下，将模型趋于重力均衡的状态作为模型的初始应力状态，然后施加边界条件．经过100个时间步(时间步长为100年)，即1万年的加载作用，以最终得到的稳定地壳构造应变场和应力场作为研究地区的背景力学环境，来计算此后模型的应力应变以及主要断层的应力状态．模拟实验过程中，依据安宁河—则木河—小江断裂带地区的现今地壳运动场、地壳构造应力场、断层滑动速率以及在该边界条件和重力作用下该地区的地壳隆升和沉降特征等方面对计算结果的可靠性进行检验．

2.1   水平方向地壳形变特征

近年来，GPS观测结果提供了高精度的中国大陆地壳运动定量数据，并且取得了一批基于GPS观测结果的研究成果(王琪等，2002； 王敏等，2008)．GPS观测结果显示，川滇活动地块的主要运动方向为SE150°—160°，北部的鲜水河附近区域的运动方向约为SE120°，而南部的安宁河—则木河—小江断裂带运动方向变为SE165°．华南地块相对川滇地块比较稳定，其内部没有明显的褶皱和活动断裂，地震活动性相对较弱．GPS观测结果表明，华南地区和东南沿海整体向E和ESE方向运动，其内部的差异运动并不明显．图 3给出了安宁河—则木河—小江断裂带地壳形变速率的计算结果和观测结果. 可以看出： 其计算值与测量值在方向和大小上都比较接近； 安宁河—则木河—小江断裂带附近区域的水平地壳形变方向为SE向．

图  3  地壳形变速率的计算值与测量值的对比

Figure  3.  Comparison of calculated crustal deformation rates with observed ones

下载: 全尺寸图片 幻灯片

2.2   地壳应力场特征

构造应力场的研究结果(许忠淮，2001； 崔效锋等，2006)一致认为青藏高原北部地区的最大主压应力方向为NNE，从青藏高原的东北部到东南部，其最大主压应力方向由NE，ENE，近E-W及SE转为SES，呈现出以青藏高原为中心的辐射状图象．为了进一步说明本文模型的可靠性，图 4给出了研究区域的最大主压应力场特征，并将其与研究区域的地壳构造应力状态进行了对比．其中测量值引自世界应力分布图 2008年发布的最新结果(Heidbach et al，2010)．对比分析表明，本文以GPS观测资料插值产生边界条件所得到的中地壳水平方向最大主压应力方向数值模拟结果，与实际观测和研究资料基本一致．

图  4  最大主压应力方向模拟结果与观测结果的对比

Figure  4.  Comparison of calculated maximum principal compressional stress directions with observed ones

下载: 全尺寸图片 幻灯片

3.   数值模拟结果

3.1   川滇地块东边界中南段应力状态

本文从应力角度分析了川滇地块东边界中南段(即安宁河—则木河—小江断裂带)的剪应力场分布特征，其模拟结果如图 5a所示． 可以看出，安宁河—则木河—小江断裂带地区存在多个应力集中的地点，如东川、巧家、西昌及石棉附近的应力场明显偏大．闻学泽(2000)综合断裂几何结构活动习性、历史地震及地球形变资料，对川滇地块东边界上特性部位的终止破裂扩展条件进行了分析，得到了5处较大的障碍体，其中有3处在本文的研究范围内： ① 新民附近区域，位于鲜水河断裂最南端与安宁河断裂的转折区； ② 西昌附近区域，从安宁河断裂—则木河断裂的转折区； ③ 巧家附近区域，位于则木河断裂与小江断裂之间．本文从应力分析的角度进一步验证了该结论； 同时还发现，在东川附近也有非常明显的应力集中现象．上述模拟结果与现场资料相互验证，进一步加深了对该区域的认识，即在现今区域应力场条件下，川滇地块东边界中南段表现出明显的力学分段特性．这对理解该区域的强震发生条件具有十分重要的意义．

图  5  (a) 安宁河—则木河—小江断裂带地区剪应力分布特征； (b) 断层面上2.0≤M≤5.9地震震源分布(引自易桂喜等， 2004)

Figure  5.  (a) Distribution characteristics of tangential traction (in unit of MPa) on Anninghe-Zemuhe- Xiaojiang fault zone; (b) Epicentral distribution of 2.0≤M≤5.9 on the fault zone (after Yi et al， 2004)

下载: 全尺寸图片 幻灯片

图 5a中的红色圆点(仅标识震中位置)表示安宁河—则木河—小江断裂带自1327年以来所发生的M≥6.5强震活动，从鲜水河断裂南端—安宁河断裂—则木河断裂—小江断裂方向分析，可以看出： 在鲜水河断裂南端，有一个剪应力集中区，1327年发生过M7 ½ 强震；在新民过渡区即石棉附近，1480年发生过M7 ½ 强震； 在西昌转折区发生过3次强震，即1536年M7 ½ 、1489年M6 ¾和1732年M6 ¾ 地震；在东川应力集中区，发生过两次强震事件，即1733年M6 ¾和1966年M6 ½ 地震．另外，还有其它一些强震如1850年M7 ½ 、1713年M6 ¾ 、1500年M7、1789年M7和1909年M6 ½ 地震均发生在剪应力相对集中区域．从上述历史强震活动与剪应力的对应关系可以看出，剪应力的分布与强弱地震的活动和断层的活动性特征有很好的对应关系，并且多数强震活动均发生在特定的构造应力环境条件下．沿安宁河断裂的中段，强震都发生在剪应力集中的部位； 沿则木河断裂，强震主要发生在与其它方向断裂相交汇的部位(如西昌附近); 而在巧家附近，虽然有明显的剪应力集中区域，但迄今为止尚未发生强震．由此推断，巧家附近将是未来强震发生的可能区域，为强震活动危险区域．另外在石棉附近也有明显的剪应力集中区域，1480年发生过的强震不是在其剪应力最集中的地点，所以石棉附近也可能是未来强震活动的危险区域．

图 5给出了安宁河—则木河—小江断裂带地区剪应力分布与小震精定位结果的对比分析，图 5b给出了重新定位的小震震源分布(易桂喜等，2004)．由本文计算的断裂带剪应力(图 5a)可以看出，在安宁河—则木河段出现两处剪应力较小的断层段，即冕宁—西昌段和西昌—普格段，其中冕宁—西昌段正好与图 5b的闭锁断层带相对应，西昌—普格段则正好与1850年四川西昌地震的主破裂面相对应．同时可以看出，小震多发区(如石棉，西昌，巧家附近)与剪应力集中区对应得很好．

3.2   库仑应力触发

由上述断裂带应力状态与历史地震活动的对应关系可以看到，大多数地震特别是强震均发生在应力相对集中区域，但也有小部分强震发生在应力偏小区域，尤为明显的是1850年发生在则木河断裂上的M7 ½ 地震和1952年发生在安宁河断裂的 M6 ¾ 地震．下面从库仑应力触发的角度探讨这一问题．

模拟地震活动引起的破裂段将作为影响研究区域应力应变的介质变化参数的输入条件，即修改强震破裂段的摩擦系数，模拟地震活动导致的断裂带应力/应变状态和应变能的变化情况，以研究强震活动对研究区域内的几条主要活动断裂带今后强震活动危险性的影响． 本文的输入地震主要依据Wen等(2008)整理的安宁河—则木河—小江断裂带1327年以来所发生强震的时间、 地点及震级等，具体如表 2所示．

表  2  安宁河—则木河—小江断裂带地区1327年以来发生的M≥6 ½地震

Table  2.  TheM≥6 ½ earthquakes on Anninghe Zemuhe Xiaojiang fault zone since the year 1327

序号	年份	东经/°	北纬/°	M	地点
1	1327	102.08	29.73	7 ½	四川雅安
2	1480	102.21	28.86	7 ½	四川岳西
3	1489	102.3	27.89	6 ¾	四川西昌
4	1500	103.16	24.87	7	云南宜良
5	1536	102.19	28.23	7 ½	四川西昌北
6	1713	103.24	25.47	6 ¾	云南寻甸南
7	1725	103.04	25.13	6 ¾	云南嵩明—宜良
8	1732	102.4	27.7	6 ¾	四川西昌西南
9	1733	103.09	26.37	7 ¾	云南东川
10	1786	102.04	29.87	7	四川泸定南
11	1789	102.96	24.29	7	云南华宁
12	1833	103	25.0	8	云南嵩明南
13	1850	102.53	27.37	7 ½	四川西昌
14	1909	103.15	24.35	6 ½	云南弥勒
15	1952	102.18	28.41	6 ¾	四川冕宁南
16	1966	103.15	26.1	6 ½	云南东川

下载: 导出CSV 

| 显示表格

库仑破裂应力变化ΔCFS可表示为

式中：Δτ是由第一次地震引起的邻近断层的滑动方向上的剪应力变化；Δσ_n是由第一次地震引起的邻近断层的断层面法向上的正应力变化，Δσ_n>0意味着断层拉张加强；μ为断层摩擦系数；ΔP为孔隙压力的变化． ΔCFS可以用来估计一次地震能否导致另一次地震孕育过程更接近或更远离破裂失稳． 库仑破裂应力变化理论已经被成功地应用于确定那些库仑应力变化为负值(ΔCFS < 0)，即处于松弛状态断层的情况． 正的库仑破裂应力变化能使断层破裂，从而触发地震；反之，负的库仑应力变化抑制断层破裂，使地震发生的可能性降低，安宁河—则木河—小江断裂带区域因此成为“应力影区”(Reasenberg，Simpson，1992；Stein et al，1997；石耀霖，曹建玲，2010)．

表 3给出了表 2中的16次强震在后续地震地点15 km深度处的库仑应力变化值，最后一行给出了地震引起的库仑应力变化的累积量． 因为一般主震能够触发地震活动的阈值为0.01 MPa(Harris，1998)，所以表中库仑应力变化值大于0.01 MPa的以黑体标记． 可以看出，1500年M7、 1713年M63/4、1732年M63/4、 1733年M73/4、 1786年M7和1789年M7这些地震没有任何应力触发的可能，因为其库仑应力变化均小于阈值，除了1732年M63/4和1786年M7地震(因为这两次地震均发生靠近数值模型边界区域，受边界条件影响剧烈，所以予以排除)．从图 5a不难看出，这些地震均发生在应力集中非常明显的区域，进一步说明了地震的发生与断裂带的应力状态有密切的关系．

从表 3最后一行可以看出，库仑应力变化的累积量大于阈值0.01的地震分别有1480年M7 ½ 、1536年M7 ½ 、1725年M6 ¾ 、1833年M8、1850年M7 ½ 、1909年M6 ½和1952年M6 ¾ 等7次地震．其中1850年发生在则木河断裂的M7 ½ 地震和1952年发生在安宁河断裂的M6 ¾ 地震均发生在断裂带上剪应力最小的段落上，其累积阈值都大于0.01，而其它阈值大于0.01的地震则发生在相对那几个大的应力集中区以外的较小应力集中区．这一现象说明在剪应力较小的区域，强震的发生主要是应力触发起主要作用．

表  3  安宁河—则木河—小江断裂带地区的库仑应力变化

Table  3.  Variation of coulomb stress of Anninghe-Zemuhe-Xiaojiang fault zone

地震	第一列地震分别在各地震发生地点15 km深度处产生的库仑应力变化/MPa
	1480年 M7 ½	1489年 M6 ¾	1500年M7	1536年M7 ½	1713年M6 ¾	1725年M6 ¾	1732年M6 ¾	1733年M7 ¾	1786年M7	1789年M7	1833年M8	1850年M7 ½	1909年M6 ½	1952年M6 ¾	1966年M6 ½
1327年M7 ½	0.0426	-0.0019	0	0.0077	-0.0006	0.0001	-0.0006	-0.0003	-0.7669	0	-0.0001	0.0001	0	0.0088	0.0002
1480年M7 ½		0.0110	0	0.0513	-0.0003	-0.0001	0.0076	0	0.0043	0	0	0.0036	0	0.1019	0
1489年M6 ¾			0	0.0201	-0.0006	0	-0.3705	-0.0009	0	0	-0.0003	0.0101	0	0.0061	0.0002
1500年M7				-0.0001	-0.0156	0.3632	0.0004	-0.0016	0	-0.0151	0.2323	0.0004	0.0231	0	0.0021
1536年M7 ½					-0.0007	-0.0009	-0.0491	-0.0015	0.0027	-0.0002	-0.0008	0.0894	0.0003	-0.0489	0.0003
1713年M6 ¾						0.0293	0	-0.0057	0	-0.0003	-0.0059	0.0005	0.0018	-0.0001	0.0044
1725年M6 ¾							0	-0.0029	0	-0.0088	1.2354	0.0006	0.0150	-0.0001	0.0039
1732年M6 ¾								-0.0009	0	0	-0.0002	0.0447	0	0.0044	0.0002
1733年M7 ¾									0.0002	-0.0017	-0.0019	0.0066	0.0042	-0.0006	-0.2007
1786年M7										-0.0001	-0.0007	-0.0003	0.0002	0.0355	0.0004
1789年M7											-0.2677	0.0012	0.2404	0	0.0040
1833年M8												0.0049	0.2155	-0.0001	0.0354
1850年M7 ½													0.0001	0.0149	0.0028
1909年M6 ½														0	0.0001
1952年M6 ¾															0
1966年M6 ½
总计	0.0426	0.0091	0	0.0790	-0.0178	0.3916	-0.4122	-0.0138	-0.7597	-0.0262	1.1901	0.1618	0.5006	0.1218	-0.1467

下载: 导出CSV 

| 显示表格

由上述分析可以看出，本文研究区域发生的强震，一方面与区域应力环境有密切的联系，另一方面与应力触发有明显的关系，或者由于两者的共同控制，因此导致了现有的安宁河—则木河—小江断裂带复杂的分段特性．图 6给出了研究区域历史强震活动引起的累积库仑应力变化. 可以看出，安宁河断裂、则木河断裂中段、小江断裂北段及大凉山断裂的南段，库仑应力变化值比较大，其所受库仑应力触发的影响比较大，可能是未来地震活动的危险区域．

图  6  安宁河—则木河—小江断裂带地区的累积库仑应力变化

Figure  6.  The accumulative Coulomb stress (in uint of MPa) of Anninghe-Zemuhe- Xiaojiang fault zone

下载: 全尺寸图片 幻灯片

4.   讨论与结论

本文通过三维有限元数值模拟方法探讨了安宁河—则木河—小江断裂带复杂的分段特性及其与地震活动的关系．模拟过程中考虑了重力作用，并采用最新的位移消除算法解决了由于重力导致的模型巨型坍塌问题，这是本文的一个重要创新点．

在重力势能和速度边界条件的共同作用下，讨论了安宁河—则木河—小江断裂带上的剪应力分布，得到了4个重要的剪应力集中区，即石棉附近、西昌附近、巧家附近和东川附近．其中前3个应力集中区与闻学泽(2000)得到的障碍体位置相同．本文将历史强震和小震精定位的结果与断裂带上应力状态进行了对比，发现无论是强震还是小震，大多数均发生在剪应力相对集中区域．但也有个别强震，如1850年发生在则木河断裂的M7 ½ 地震和1952年发生在安宁河断裂的M6 ¾ 地震，发生在剪应力非常小的区域．为了解释这一现象，本文模拟了16次地震引起的库仑应力变化. 结果表明： 在应力明显集中的区域，地震活动主要由区域应力分布控制； 在应力非常小的区域，地震活动由库仑应力触发作用控制； 在应力相对比较明显的区域，地震活动由区域应力状态和应力触发共同控制．

综上所述，从断层面的应力分布角度来看，石棉、巧家、西昌及其附近区域是未来可能的地震高危险区； 从库仑应力触发的角度来看，安宁河段和巧家—东川段的未来地震活动性也相对偏高．

在本文撰写过程中，得到中国地震局地球物理研究所李永华研究员的指导和建议，作者在此表示诚挚的感谢．