地球物理学报  2012, Vol. 55 Issue (1): 105-116   PDF    
用单元降刚法探索中国大陆强震远距离跳迁及主体活动区域转移
杨树新1,2, 陆远忠1, 陈连旺1, 叶际阳1, 米琦3     
1. 中国地震局地壳应力研究所, 北京 100085;
2. 北京交通大学土木建筑工程学院, 北京 100044;
3. 中国科学院研究生院, 北京 100049
摘要: 中国大陆强震有成组活动、远距离跳迁和不同时段形成主体活动地区的特征.本文利用强震发生位置处单元降刚法,对中国大陆地区强震的远距离跳迁和主体活动地区转移机理进行了数值模拟研究.得出的新认识是:(1)在地壳中存在初始应力场的环境中,已发生强震区部分丧失承载能力(模拟中作为显著降低单元组的弹性模量来处理),可以引起大范围兆帕量级的应力场调整,它是后续强震可远距离跳迁的主要因素;(2)一个活动期中,中国大陆强震主体活动地区及其迁移,受主要活动断层分布、初始应力场和边界载荷的配置方式的综合影响,但在十年左右的活动幕中,边界载荷的配置方式可能是控制主体活动地区及主体活动地区转移的重要因素.
关键词: 单元刚度      应力场调整      数值模拟      强震迁移      初始应力场      中国大陆     
The mechanism of long-distance jumping and the migration of main active areas for strong earthquakes occurred in the Chinese continent
YANG Shu-Xin1,2, LU Yuan-Zhong1, CHEN Lian-Wang1, YE Ji-Yang1, MI Qi3     
1. Institute of Crustal Dynamics, China Earthquake Administration, Beijing 100085, China;
2. School of Civil Engineering, Beijing Jiaotong University, Beijing 100044, China;
3. Graduate University of Chinese Academy of Sciences, Beijing 100049, China
Abstract: Strong earthquakes occurred in the Chinese continent are usually characterized by grouped activity, long-distance jumping migration, and different main activity areas formed in different times. In the present study a three-dimensional (3-D) finite element model was set up for the Chinese continent involving surface topography, major active fault zones as well as initial stress field to study the mechanism of the long-distance jumping migration and the migration of main active areas for the strong earthquakes occurred in the Chinese continent by using numerical simulations, after considering the birth and death of element groups under different boundary loading actions. Our results show that (1) in an environment where always exists an initial stress field in the Earth's crust, the area where a strong earthquake has occurred has no longer the ability to support high stress (we suppose that the moduli of elasticity of elements within the area be reduced), which can lead to stress adjustment to the order of MPa in a large area, and may be one of the main factors affecting the long-distance jumping migration of the follow-up strong earthquakes, and (2) it is hard to accurately predict where the follow-up strong-earthquakes will occur because it could be affected by many factors such as the loading manner, geological structure, active fault zones, initial stress field, stress field adjustment induced by strong earthquakes, but in an active period of earthquakes, the major factor to affect the activity area and migration of strong earthquakes is the boundary loading configurations. These results suggested that it is helpful to predict the trend of strong-earthquake migration by investigating various kinds of boundary loading manners and the relationship between stress field adjustment induced by strong earthquakes and the regions where strong-earthquakes occurred.
Key words: Element stiffness      Stress field adjustment      Numerical simulation      Strong earthquake migration      Initial stress field      Chinese continent     
1 引 言

中国大陆地区强震活动有一些基本特征:强震沿主要活动断裂带或活动地块边界分布;强震活动在时间轴上具有活跃与平静相交替的特征.在活跃期内,强震呈成组的群体活动状态发生.强地震活动的各个活跃期都有其强震的主体活动地区[1-2];相继发生的强震震中,除强余震外,迁移的距离可达2000km 以上.在变化缓慢的边界驱动或上地幔拖拽作用下何以能发生短时间的(数天-数月)强震震中千米量级的迁移?成组活动强震的主体活动区的转移的原因是什么?研究这些问题对于指导年度尺度的强震预测和防震减灾是重要的,它也是推进从经验性、统计性地震预报向动力的、物理的地震预报转变的必要探索.

近十年来,用临界地壳状态下强震本身产生的库仑应力扰动来解释后续强震活动进行了大量的研究.从1997年开始,提出了多种模型计算地震产生的库仑破裂应力,研究了库仑破裂应力对后续地震的影响,震例研究结果分布在世界的不同地区[3-5].2008年5月汶川8.0 级大震后,Toda等[6]计算了汶川地震产生的库仑应力,认为库仑应力可能导致鲜水河断层、东昆仑断层大部以及部分岷江断层更加接近临界状态,其地震发生率增加.SteinRS.[7]乐观地认为:“通过计算断层活动所发生的应力,科学家发现,大地震之间是会互相影响的.这个令人兴奋的发现,让预测地震成为可能".但是,也有研究结果不支持这些观点,发现一些后续地震发生在库仑应力减小的区域.万永革等[8]研究得出中国及邻区的地震(以走滑断层为主)不容易观测到地震静态应力触发现象;周仕勇[9]定量计算了川西地区一条断层的破裂产生的库仑应力在其他断层面上的投影,得出5条断层中发生的强震的相互影响,加载、减载和影响很小的情况都存在.为了改进库仑应力对后续地震触发作用的研究,一些研究者利用临界状态下地震矩加速释放地区来约束库仑应力增加的区域,以提高对应效果[10-11],用断层滑动的速率、状态依从模型或下地壳的黏滞性等来说明库仑应力延迟触发后续地震[12-13],但是强震产生的静态应力扰动(或触发效应)随距离的衰减很快,难以解释强震远距离迁移和主体活动区域的形成和转移的特征.

将数值模拟引入应力场扰动、地震活动迁移的研究,可以在研究中涉及复杂的介质结构来模拟地震动力环境.Cianetti等[14]研究1992年Landers和1999Hector矿山地震的同震、震后应力变化,其目的是探讨重力载荷、库仑摩擦以及区域应力场的方向、大小和随深度的变化是如何影响断层的位错和应力变化的.Zller等[15]建立了在弹性半空间中嵌入一个由脆性和蠕变特性网格组成的非均匀矩形断层的模型,得到了地震成丛、强震的准周期性等各种观察到的地震活动性特征,还得到了与近年来国外广泛研究讨论的地震矩加速释放模型相一致的结果,并支持地震活动的临界点动力学的观点.陈连旺等[16]通过建立较精细的川滇地区三维有限元模型,模拟了川滇地区主要活动断裂的强震活动对于其他活动断裂潜在强震孕育进程的库仑破裂应力加卸载效应.陶玮等[17]对中国强震活动的主体地区的形成机制进行了探讨,他们建立了二维黏弹性有限元模型,模拟韧性层的变形及应力场分布.认为岩石圈韧性层的变形控制了强震活动主体地区的主要位置和形态,而周围板块边界作用是形成主体地区的主要驱动力.郑勇等[18]用有限元方法探讨了欧亚大陆的碰撞对中国大陆岩石圈形变和应力场的影响以及它们与强震活动性的关系.结果认为,印度板块的碰撞对中国大陆的强震活动性有重要影响.Parsons[19]用三维有限元模型研究了1906 年美国旧金山地震的应力恢复过程.他指出,该地震产生的应力影区(库仑应力减小区)推迟了预估的本地强震发生的时间.陈连旺等[20]根据板块构造理论,结合华北地区应力场演化的历史进程,利用华北地区地壳构造应力场3D 有限元模型,数值模拟华北地区动力边界作用力发生改变时其内部构造应力场所发生的相应变化,定量研究了印度板块和太平洋板块作用力变化所引起的应力场和能量场的分区加卸载效应.

上述众多的数值模拟研究得出的共同认识是,对强震活动有重大影响的因素,除边界驱动作用外,强震的发生(发震断层的滑动)及其引起的应力场扰动是十分重要的.强震产生的应力场的动态和静态扰动,对其后的余震的时空分布、对于处于临界状态地区的远程小震活动(瞬时触发的小震和改变小震发生率)确有重大影响.但是,无论是简单的解析计算或较复杂的数值模拟,用位错产生的应力场扰动得到的结果都是局部的,都难以解释中国大陆地区强震空间分布的大的格局变化.Lu 等[21]曾经用单元生死的方法数值模拟了强震区完全卸载引起的应力场大范围调整,初步解释了上述现象.但是,在实际地壳介质中很少出现强震后震源区完全没有承载能力的状况,并且这种丧失的承载能力还会逐渐恢复.在最近20年中,很多研究人员依据震中附近区域的微震观测,成功地证实由强震引起的震中区表层和地震断层数千米深处,地壳介质的物性变化(波速变化,先减小然后逐渐缓慢回升)[22-24].徐煜坚等[25]在研究邢台、渤海、唐山等强震迁移时,利用降低单元弹性模量,用有限元数值模拟了强震震源区应力场的变化,解释在一个地震序列中应力的转移.本文应用文献[21]建立的中国大陆有限元模型,通过逐次降低最高应力单元组的弹性模量的方法,在模拟强震的发生并使初始应力场中的震源区部分丧失承载能力导致强震迁移方面进行了更深入的研究.

2 3D 数值模型

本文按文献[21]的做法建立了中国大陆三维有限元数值模拟模型.以中国大陆的实际边界为参考向外扩展数百千米,以消除边界效应对大陆内部的影响.模型分7层,共有28640个六面体单元,33165个节点(图 1).为了显示清楚,在垂直方向显示放大了15倍.图中包括20条主要活断层,其位置参考了马杏垣[26]主编的中国岩石圈动力学图集.图中序号为加外载荷的边界的序号.

图 1 有限元模型示意图 Fig. 1 Sketch map of finite element model

表 1给出了各层介质的力学参数初始弹性模量E,泊松比υ,表中注明了赋予20 条断层的弹性模量.表 2给出了计算初始应力场采用的边界条件.

表 1 模型参数设置 Table 1 Parameter values used in the model
表 2 边界载荷(边界序号见图 1) Table 2 Boundary load (the boundary H) referred to Fig. 1)

图 2a绘出了用商业软件ansys计算得到的中国大陆第3层(10~15km)单元地壳等效应力场的分布.图中等值线值是等效应力,单位为Pa其定义为:

(1)

图 2 (a)初始等效应力及中国7级以上强震震中分布;(b)地表位移向量方向与GPS实测向量方向的比较.在360°范围内,55%的方向差值小于40° Fig. 2 (a) Initial equivalent stress and the distribution of earthquakes(M≥7) epicenter in China; (b) The comparison of ground displacements between calculated values and observed values

式中,σxσyσzσxyσyzσzx为直角坐标中的各应力张量的分量.图 2a中用蓝色实心圈绘出了中国有记录以来至2008年9月的7级以上强震的震中.粗略比较可以看出,这些强震大都发生在高应力区和主要断层带上.下面的计算和分析均用这些强震作为初始应力场,来探讨强震的发生和现今不同边界载荷对中国大陆强震空间分布和迁移的影响.由于加入了包括地形影响的重力作用和最高达100 MPa的周边边界力,因此计算得到的地壳内15km 深部的等效应力水平约为70~1000 MPa左右,与实测原地应力测量和由自重推算的应力量级一致[27].

图 2b绘制了在上述初始应力状态下,给模型边界施加由实测GPS值插值于表面边界节点的位移,计算得到的表面位移向量方向场(蓝色箭头).实测GPS值为李延兴等[28]提供的由2004年与2007年2期的观测结果扣除欧亚板块的刚体运动得到的速度方向场(红色箭头).为了便于比较,图中仅绘出了与GPS实测点最近距离的节点(994个)的位移向量方向.定量统计表明,在全部994个比较点位中,两种向量差值小于40°的占55%.

比较图 2a图 2b强震空间分布和GPS 位移场与由模型计算的等效应力场和位移场,可以认为本模型有一定的合理性,可以作为进一步研究的基础.

3 震源区部分丧失承载能力引起的强震远距离跳迁的数值模拟

一个强震的发生会在近距离范围内产生附加的应力应变扰动.邱泽华等[29]计算了不同震级地震引起10-8应变变化(应力约为0.0001 MPa)的最大传递距离.一个7~8 级地震最大传递距离为410~1300km.强震后的大量余震与这些附加的应力扰动(库仑应力)有密切关系.但是,伴随强震及其大量余震发生后,另一个十分重要而常被忽略的事实是强震区丧失承载能力并引起大范围的应力调整和转移.震后发震断层的余滑和强震带上的应力松弛[30]以及中国大陆地区强震在发震断层段重复时间为数十年到数百年,乃至几千年[31-32]也佐证了强震区在相当长的时间内丧失承载能力、不能重复发生强震.实际上,一个强震的发生,由地震波辐射的能量仅占地震释放能量的5% 左右[33],绝大部分能量损失在震源区的其他非弹性的和非力学的能量耗损中.这就导致震源区及其附近区域断层在一段时间部分丧失了承载能力,并引起非弹性非动态应力场大范围调整和应力集中区转移.为此,我们在前述模型的基础上,把已经产生的应力场和应变场作为初始背景场,重新在各边界施以年速率为厘米级的位移作为外载荷,研究显著降低模型中当前包含最大应力单元的一组单元(从其表层到底层的单元)的弹性模量后应力应变场的变化,用它来模拟强震发生引起的应力场调整.用被显著降低弹性模量的一组单元模拟已发生强震的震源.在调整后的应力场中再寻找包含最高应力单元的单元组,作为后续强震的发震位置,在不变的位移载荷和原有的初始背景应力状态下,再显著降低这组单元的弹性模量,又引起新的应力应变场的调整和应力集中区转移.与此同时,上一次被降低了弹性模量的单元组的弹性模量恢复性增加其弹模数值.我们把这些逐次被降低弹模,继后又逐次恢复性增加弹模的单元组称为“模拟地震",它们同时也引起应力场的不断大范围调整.与文献[21]把模拟地震的弹性模量降为10-6(“杀死相应单元组")不同,我们试验了多种降低“模拟地震"的单元组的弹性模量并采用随后每步逐渐恢复其弹性模量的方法,在单元弹性模量降低或恢复中,单元的泊松比也作相应变化,以使单元的剪切模量、体积模量合理.表 4举例列出了三种降刚方式下20步运行中“模拟地震"单元组的弹性模量的取值.第1 行为加载步序号,第2行“模拟地震"发生时,该单元组弹性模量降为0.1GPa,随后,它们逐渐恢复,直至20步结束时,恢复到9.6GPa(所有“模拟地震"作同样降刚值处理).第3 行举例列出了24986 号单元为“模拟地震"发生时,其弹性模量降为其原值的1/4.47,随后,它们逐渐恢复,直至20 步结束时,恢复到原值.这后一种降刚和逐渐恢复至原刚度的取值方法,是考虑强震的发生,非弹性的和非力学的能量耗损约为地震波辐射能的20 倍;我们用“模拟地震"单元降刚至其原值的1/4.47来模拟这种效应.

根据计算结果,图 3a绘出了按表 3中方式8的位移加载下,经过20次发生模拟地震后等效应力分布.图 3b在模型轮廓(图中示意显示了断层带位置)平面图上,标出了20次模拟地震的序号,边界位置编号见图 1,位移加载于该边界从地表到深40km的所有单元的外表面节点上.图 4 绘制了我国大陆第4活动期(1966~1976 年)7 级以上16 次强震的震中分布,图中也标注了地震迁移序号.

表 3 边界荷载配置方式(单位:m) Table 3 Boundary loads of the model (unit: m)
图 3 等效应力分布与模拟地震发生顺序 (a)等效应力分布(表 3中第8种加载方式,按表 4第2行方式降刚);(b)模拟地震发生顺序(表 3中第8种加载方式,按表 4第2行方式降刚); (c)等效应力分布(表 3中第8种加载方式,按表 4第3行方式降刚);(d)模拟地震发生顺序(表 3中第8种加载方式,按表 4第3行方式降刚) Fig. 3 Equivalent stress distribution and the sequence of simulation earthquake occurrence (a) Equivalent stress distribution(Tab 3-No. 8 boundary load) ; (b) The sequence of simulation earthquake occurrence (Tab 3-No. 8 boundary load) ; (c) Equivalent stress distribution (Tab 3-No. 8 boundary load) ; (d) The sequence of simulation earthquake occurrence (Tab 3-No. 8 boundary load).
图 4 1966〜1976年中国大陆7级 以上地震震中分布(按发震时序) Fig. 4 The earthquake(M^7) epicenter distribution map between 1966〜1976 in Chinamainland (litt by the earthquake time sequence)

图 3可以看出,第8 种加载方式下,20 次模拟地震的空间分布与我国大陆第4活动期实际发生的强震有相似的特征:一是位置分布在东、中、西部广大范围;二是强震可能相继成丛发生,三是强震可远距离跳迁.图 5a图 5b分别绘制了这20次模拟地震和实际发生的16次地震震中距的直方图,它们有较好的相似性.距离在100km 以内的,模拟地震和实际地震分别占总数的16% 和20%;距离在1000km 以上分别占总数的66%和79%.最远的跳迁分别是3800km 和3901km.为了比较,图 3c图 3d也绘出了按表 4第3行的降刚模式,计算得到的经过20次发生模拟地震后等效应力分布和20次模拟地震的序号,与图 3a图 3b对比,两者结果相似.多种降刚模式试算比较(包括文献[21]中用“杀死单元来模拟地震")可知,尽管模拟地震的位置和发生顺序以及模拟地震引起的应力扰动状态有变化,但模拟地震的发生均引起大范围兆帕量级的应力场调整,是后续强震可远距离跳迁的主要影响因素.本文后面,我们仅讨论了表 4第二行降刚方式的结果.

图 5 相继震中距直方图 (a)第8加载方式下的模拟地震(总数20次)(b)1966〜1976年中国7级以上地震(总数16次) Fig. 5 Histogramof the distance between the earthquakes (a) The simulation earthquakes (total 20) in No. 8 boundary load; (b) Earthquakes occurred in China mainland from 1966 to 1976 (M>7, total 16)
表 4 各步最高应力单元组的弹性模量 Table 4 Elastic modulus in highett stress element set of each sub step
4 震源区部分丧失承载能力引起的应力场的调整变化

为了研究模拟地震的空间分布形成机理,我们重点探讨它们引起的应力场的变化.表 5 各列顺序列出了在图 3所示各次模拟地震的序号、最大应力单元号、相继最大应力单元的转移距离、该模拟地震未发生时与其前一步单元等效应力增量、以及模拟地震引起震源处的等效应力变化及其百分比.图 6举例绘出了第3号模拟地震(22565号单元组)引起的等效应力调整变化(模拟地震前后两步的等效应力差)的等值线图.

表 5 模拟地震发生前后地震所在单元的等效应力变化 Table 5 The changes of equivalent stressin the elements which the simulation earthquake occurred
图 6 第3号模拟地震(22565号单元组)引起的等效应力调整变化(Pa) Fig. 6 The change of equivalent stress (Pa) triggered by No. 3 simulation earthquake (No. 22565 element set)

可以看出,在本模型的初始应力场中模拟地震失去部分承载能力引起的应力调整很大,第3 号模拟地震在它的震源处产生了72%的应力降,而在相距2315km 的第4号模拟地震处,引起达5.82 MPa的等效应力升高,这是位错引起的库仑应力变化不可比拟的(兆帕级的应力变化可能是强震及强震主体活动区可以远距离迁移的重要原因).模拟地震部分失去承载能力引起的应力调整的空间分布,总的趋势仍然是距离越近,引起的应力调整量越大.模拟地震在其附近引起很大的应力变化,从而相继发生后续强震,是直观、可以理解的(见表 4 中2 号,6 号,10号等模拟地震,引起后续近处地震处数百兆帕的变化).但这种应力调整的空间分布不规则,不是简单的距离的函数.为什么在距模拟地震2000km 以上的地方,兆帕量级的应力调整(见表 4 中,1,4,9号等共10个)会使之达到最高等效应力,从而变为后继的模拟地震(远距离强迁移的强震)?甚至,在6,8,10 号模拟地震位置,各自前一个模拟地震(即4,7,9号地震)引起的等效应力变化减小(兆帕量级),但这些单元仍然变为该步模型中最高等效应力单元.实际上,最高应力单元位置的变化受边界作用大小、方式、模型结构、断层结构和初始应力场等多因素组合的制约,20次逐次模拟地震丧失承载能力的影响是累加的,其先后次序也有关系.因此,要由前一个强震推算后续强震的准确位置是复杂困难的.仔细研究每一个模拟地震前震源处应力变化,大致有三种情况:第一种是模拟地震发生前,该处等效应力一直变化不大,突然上升模型最高等效应力(3,7,11 号地震),全部近距离迁移地震均属于此类;第二种是多数模拟地震(如,2,4,5,15,17,18,19,20号),震前等效应力呈折线上升或起伏上升至模型最高等效应力;第三种是模拟地震(包括6,8,10,12,13,14,16号模拟地震)发生前几步等效应力变化不大,它们上升至模型最高应力单元是由于其他单元应力相对减小所致.由此推断,强震前震源处应力逐步起伏增高或突然增高是较普遍的现象.远距离跳迁的强震,均属于后两种.图 7举例表示这3种模拟地震等效应力顺序变化过程.

图 7 模拟地震单元等效应力随演化顺序变化过程 (a)第11号模拟地震单元;(b)第18号模拟地震单元;(c)第12号模拟地震单元. Fig. 7 The change of the simulation earthquake elements' equivalent stress with the evolution sequence (a) No. 11 simulation earthquake elements; (b) No. 18 simulation earthquake elements; (c) No. 12 simulation earthquake elements.

由此可见,强震远距离跳迁主要是因为在某区域已经积累了高应力的条件下,已发生强震区域部分丧失承载能力引起的兆帕量级的应力变化,使高应力区域中应力容易集中处应力更升至最高所致.表 5中6,7号模拟地震单元的降刚没有引起等效应力的降低,反而引起等效应力升高,可能因为它们处于大断层的最端部、应力集中程度很高所致;6号地震后,相邻仅6km的7号模拟地震发生.

由于在模拟中我们设计了模拟地震发生后其所在单元组的弹性模量会随后逐渐恢复性增长,计算结果也显示,这些单元组虽然其应力大小受其后别处模量地震发生的影响出现起伏升降,但总的趋势是应力逐渐上升.例如第1号地震应力降是16.3%,所在24988单元在第20步计算后,其等效应力恢复至未发生模拟地震前的84%.由于本次研究集中于模拟强震活动的空间转移,未考虑其时间因素,模型中的黏滞性影响和应力恢复的时间效应是今后进一步研究的问题.

5 边界位移加载方式对强震活动主体地区转移影响的数值模拟

我国大陆应力应变场的分布主要受到印度洋板块和太平洋板块的驱动作用的控制,那么,大陆的强震主体活动地区的分布和转移是否主要受到从边界传递的不同的加载组合的影响?文献[3]建立了二维黏弹性有限元模型来模拟韧性层的变形及应力场分布.他们的研究结果说明了边界作用对大陆内部的应力场有重要影响.但是,仅靠缓慢的边界作用和下地壳、上地幔的黏滞性难以解释十年作用时间尺度的大陆强震活动的时空变化.我们研究发现,虽然模拟地震远距离迁移受多种因素的影响,难以事先预测,但是边界加载的配置方式与模拟地震引起的应力场格局的变化共同作用,可以用来对强震活动主体地区的转移提供解释.我们用9 种不同位移边界配置研究了本文模型中模拟地震引起的应力场演变和模拟地震空间分布的转移.每一种边界配置,仍然进行20步模拟地震的演算.表 3列出了9种加载方式的设置.图 8给出了8种位移加载配置下20次模拟地震的空间分布图像(第8 种加载方式已在图 3绘出).图中标出了各模拟地震的顺序号,每幅图的左下角表示了加载方式设置号(见表 3)并以图边的箭头示意出加载方式.

图 8 不同加载方式下模拟地震位置随演化顺序的分布 Fig. 8 Distribution of the simulation earthquakes location in different boundary loads

图 8中同样可以看出,不管边界加载方式如何,除了图中北边界由于xy方向固定或加载边界拐角处应力集中引起的模拟地震外(它们表现为这9种加载方式下,同样位置均有同样模拟地震发生),模拟地震都集中在断层带上、特别是其端部,且远距离跳迁的现象仍普遍存在.此外,可以看到模拟地震的总体空间分布与加载方式关系密切:① 总体看,在单边界加载方式下(第1~6 种加载),随着加载的南移、东移,模拟地震主体地区也由西向南、向东迁移,每种加载的第1、2 号模拟地震所在位置都离相应加载边界近.②第1和第6种加载作用下,模拟地震都发生在北纬35°~45°之间且分别发生在东经95°的西东部,没有交叉.③第5边界东南-西北方向的加载促使华北模拟地震的发生.④ 第3 边界西南-东北方向的加载促使西南地区模拟地震的发生.⑤在多边加载方式下(第7~9 种加载),模拟地震东、中、西部都有.合成单方向加载的结果可以大致圈出相应多边界加载模拟地震大致发生的主要区域.例如第1、2和第3 种加载就大致指出第7 种加载方式下模拟地震发生的位置.加入第5种加载,就大致可以勾画出第8种加载方式模拟地震发生的位置.第2、4种加载方式模拟地震分布较广.第9种加载方式的结果说明加入东北地区向西的位移约束,对模拟地震的发生位置影响不大.

6 结论和讨论

利用数字等高线的资料建立了包含地面高程、主要活动断层以及初始应力场的中国大陆三维有限元模型,在不同的边界加载作用下,对中国大陆地区强震的远距离跳迁和主体活动地区转移机理进行了数值模拟研究.得出了两个认识:(1)在地壳中始终存在初始应力场的环境中,已发生强震区部分丧失承载能力,可以引起大范围兆帕量级的应力场调整,它可能是后续强震远距离跳迁的主要因素;(2)尽管后续强震的位置受加载方式、地质构造、活断层、初始应力场以及发生强震的顺序等多种因素影响,难以准确预测,但是一个活跃幕中,中国大陆强震(模拟地震组)主体活动地区及其迁移,受边界载荷的配置方式影响很大.研究多种组合的边界加载与强震发生区域的关系,有助于对其做出趋势预测.这个结论与用杀死单元来模拟强震的研究结果一致.本文的不同在于,用降低强震震源区的弹性模量使之部分丧失承载能力,更合乎地震发生的实际,而且可以进一步研究地震的应力降、震后的应力恢复和地震的重发.

本文对强震空间分布及迁移的机理研究提出了一种新的探索性思路,得到了有参考价值的结果.但有两个问题值得讨论,一是选择最高应力单元而不是达到其破坏强度的单元为模拟地震,不够合理.一方面没有科学合理的方法给定各单元强度,另一方面人为确定一个阈值应力使数个单元同时被杀死,也不符合中国大陆强震没有同时发生的实际.本文选择最高应力单元降低其弹性模量模拟地震是当前一种可行的选择.二是虽然我们设计了模拟地震两种弹性模量降低恢复的方式(见表 4),采用中、高、低3种初始应力场,也在模型底部加入了20cm 的位移约束进行数值模拟试验,得出的基本结论相同.但要应用于强震预报的实践,具体预测强震的迁移过程和活动期幕的主体活动地区、几何模型、物性设置、边界加载方式以及地幔的拖拽作用都是值得深入研究的课题.例如,修改模型,加入发生在印度板块东边界的南亚的一些特大地震对中国大陆强震活动的影响[34],对于汶川8.0级大震后强震活动的转移都是重要的,值得我们进一步研究.此外,探索强震活动(模拟地震组空间分布和转移)的时间特征,而不是本文中应力场随载荷步的变化,同时利用模拟地震组复活研究强震再孕育的时间过程;细分单元尺度(水平尺度平均为10km 左右),进行大规模并行计算,探讨制约中强地震活动规律的机理;用实测资料(如边界GPS 位移资料,应力测量资料),对初始应力场和边界加载作更实际的约束等是我们进一步研究的目标.

致谢

感谢李延兴研究员为本文提供了2004 与2007年2期的GPS观测资料.最后,对审稿者的宝贵意见致以诚挚谢意.

参考文献
[1] 马宏生, 张国民, 刘杰, 等. 中国大陆及其邻区强震活动与活动地块关系研究. 地学前缘 , 2003, 10(Suppl.): 75–80. Ma H S, Zhang G M, Liu J, et al. Correlation between strong earthquake activity and active crustal-block in China mainland and its adjacent regions. Earth Science Frontiers (in Chinese) , 2003, 10(Suppl.): 75-80.
[2] 洪汉净, 刘培洵, 于泳, 等. 百年来中国大陆强震活动的微动态分期及其空间分布. 地震地质 , 2003, 25(3): 394–402. Hong H J, Liu P X, Yu Y, et al. Temporal and spatial distributions of strong earthquakes in China's continent in the past century. Seismology and Geology (in Chinese) , 2003, 25(3): 394-402.
[3] Stein R S. The role of stress transfer in earthquake occurrence. Nature , 1999, 402(6762): 605-609. DOI:10.1038/45144
[4] Toda S, Stein R S. Did stress triggering cause the large off-fault aftershocks of the 25 March 1998 Mw=8.1 Antarctic plate earthquake?. Geophys. Res. Lett. , 2000, 27(15): 2301-2304. DOI:10.1029/1999GL011129
[5] Nalbant S S, Steacy S, Sieh K, et al. Earthquake risk on the Sunda trench. Nature , 2005, 435(7043): 756-757. DOI:10.1038/nature435756a
[6] Toda S, Lin J, Meghraoui M, et al. 12 May 2008 M=7.9 Wenchuan, China, earthquake calculated to increase failure stress and seismicity rate on three major fault systems. Geophys. Res. Lett. , 2008, 35: L17305. DOI:10.1029/2008GL034903.
[7] Stein R S. Earthquake conversations. Scientific American , 2003, 288(1): 72-79. DOI:10.1038/scientificamerican0103-72
[8] 万永革, 吴忠良, 周公威, 等. 地震静态应力触发模型的全球检验. 地震学报 , 2002, 24(3): 302–316. Wan Y G, Wu Z L, Zhou G W, et al. Global test of seismic static stress triggering model. Acta Seismologica Sinica (in Chinese) , 2002, 24(3): 302-316.
[9] 周仕勇. 川西及邻近地区地震活动性模拟和断层间相互作用研究. 地球物理学报 , 2008, 50(1): 165–174. Zhou S Y. Seismicity simulation in Western Sichuan of China based on the fault interactions and its implication on the estimation of the regional earthquake risk. Chinese J. Geophys. (in Chinese) , 2008, 50(1): 165-174.
[10] Bowman D D, G King C P. Accelerating seismicity and stress accumulation before large earthquakes. Geophys. Res. Lett. , 2001, 28(21): 4039-4042. DOI:10.1029/2001GL013022.
[11] Mignan A, King G C P, Bowman D. Accelerating seismicity and stress accumulation model. Geophysical Research Abstracts , 2007, 9: 26-44.
[12] 刘桂萍, 傅征祥, 李钢, 等. 摩擦状态——速率依从的区域地震触发模型研究. 地震 , 2004, 24(1): 177–183. Liu G P, Fu Z X, Li G, et al. Study on the frictional state-velocity dependent model of regional triggering seismicity. Earthquake (in Chinese) , 2004, 24(1): 177-183.
[13] 万永革, 沈正康, 曾跃华, 等. 青藏高原东北部的库仑应力积累演化对大地震发生的影响. 地震学报 , 2007, 29(2): 115–129. Wan Y G, Shen Z K, Zeng Y H, et al. Evolution of cumulative coulomb failure stress in Northeastern Qinghai-Xizang (Tibetan) Plateau and its effect on large earthquake occurrence. Acta Seismologica Sinica (in Chinese) , 2007, 29(2): 115-129.
[14] Cianetti S, Giunchi C, Cocco M. Three-dimensional finite element modeling of stress interaction: An application to Landers and Hector Mine fault systems. J. Geophys. Res. , 2005, 110: B05S17. DOI:10.1029/2004JB003384.
[15] Z?ller G, Hainzl S, Ben-Zion Y, et al. Earthquake activity related to seismic cycles in a model for a heterogeneous strike-slip fault. Tectonophysics , 2006, 423(1-4): 137-145. DOI:10.1016/j.tecto.2006.03.007
[16] 陈连旺, 张培震, 陆远忠, 等. 川滇地区强震序列库仑破裂应力加卸载效应的数值模拟. 地球物理学报 , 2008, 51(5): 1411–1421. Chen L W, Zhang P Z, Lu Y Z, et al. Numerical simulation of loading/unloading effect on Coulomb failure stress among strong earthquakes in Sichuan-Yunnan area. Chinese J. Geophys. (in Chinese) , 2008, 51(5): 1411-1421.
[17] 陶玮, 洪汉净, 刘培洵. 中国大陆及邻区强震活动主体地区形成的数值模拟. 地震学报 , 2000, 22(3): 271–277. Tao W, Hong H J, Liu P X. The finite element modeling of strong earthquake major regions in Chinese mainland and its vicinity. Acta Seismologica Sinica (in Chinese) , 2000, 22(3): 271-277.
[18] 郑勇, 傅容珊, 熊熊. 中国大陆及周边地区现代岩石圈演化动力学模型. 地球物理学报 , 2006, 49(2): 415–427. Zheng Y, Fu R S, Xiong X. Dynamic simulation of lithospheric evolution from the modern China mainland and its surrounding areas. Chinese J. Geophys. (in Chinese) , 2006, 49(2): 415-427.
[19] Parsons T. Post-1906 stress recovery of the San Andreas fault system calculated from three-dimensional finite element analysis. J. Geophys. Res. , 2002, 107(8): 1-13.
[20] 陈连旺, 陆远忠, 郭若眉, 等. 边界作用力变化引起的华北地区应力场分区加卸载效应. 大地测量与地球动力学 , 2007, 27(6): 10–16. Chen L W, Lu Y Z, Guo R M, et al. Regional loading-unloading effects induced by boundary force changes on stress field in North China. Journal of Geodesy and Geodynamics (in Chinese) , 2007, 27(6): 10-16.
[21] Lu Y Z, Yang S X, Chen L W, et al. Mechanism of the spatial distribution and migration of the strong earthquakes in China inferred from numerical simulation. Journal of Asian Earth Sciences , 2011, 40(3): 990-1001.
[22] Chao K, Peng Z G. Temporal changes of seismic velocity and anisotropy in the shallow crust induced by the 1999 October 22 M6.4 Chia-Yi, Taiwan earthquake. Geophys. J. Int. , 2009, 179(3): 1800-1816. DOI:10.1111/gji.2009.179.issue-3
[23] Zhao P, Peng Z G. Depth extent of damage zones around the central Calaveras fault from waveform analysis of repeating earthquakes. Geophys. J. Int. , 2009, 179(3): 1817-1830. DOI:10.1111/j.1365-246X.2009.04385.x.
[24] Rubinstein J L, Uchida N, Beroza G C. Seismic velocity reductions caused by the 2003 Tokachi-Oki earthquake. J. Geophys. Res. , 2007, 112: B05315. DOI:10.1029/2006JB004440.
[25] 徐煜坚, 罗焕炎, 虢顺民, 等. 华北北部地区地质模型与强震迁移. 北京: 地震出版社, 1985 . Xu Y J, Luo H Y, Guo S M, et al. The Geological Model and Strong Earthquake Migration in Northern Huabei (in Chinese). Beijing: Seismological Press, 1985 .
[26] 马杏垣. 中国岩石圈动力学地图集. 北京: 中国地图出版社, 1989 . Ma X Y. Lithospheric Dynamics Atlas of China (in Chinese). Beijing: Sinomaps Press, 1989 .
[27] 吴满路, 张春山, 廖椿庭, 等. 青藏高原腹地现今地应力测量与应力状态研究. 地球物理学报 , 2005, 48(2): 327–332. Wu M L, Zhang C S, Liao C T, et al. The recent state of stress in the central Qinghai-Tibet Plateau according to in-situ stress measurements. Chinese J. Geophys. (in Chinese) , 2005, 48(2): 327-332.
[28] 李延兴, 杨国华, 李智, 等. 中国大陆活动地块的运动与应变状态. 中国科学(D辑) , 2003, 33(Suppl.): 65–81. Li Y X, Yang G H, Li Z, et al. Movement and strain conditions of active blocks in the Chinese mainland. Science in China(Series D) (in Chinese) , 2003, 33(Suppl.): 65-81.
[29] 邱泽华, 石耀霖. 地震应变变化的理论量级与最大传递距离关系. 大地测量与地球动力学 , 2003, 22(4): 99–106. Qiu Z H, Shi Y L. Relationship between theoretic magnitude and longest propagation distance of seismic strain change. Journal of Geodesy and Geodynamics (in Chinese) , 2003, 22(4): 99-106.
[30] Avage J C, Langbein J. Postearthquake relaxation after the 2004 M6 Parkfield, California, earthquake and rate-and-state friction. J. Geophys. Res. , 2008, 113: B10407. DOI:10.1029/2008JB005723.
[31] 闻学泽. 中国大陆活动断裂段破裂地震复发间隔的经验分布. 地震学报 , 1999, 21(6): 616–622. Wen X Z. Distribution of empirical recurrence intervals of segment-rupturing earthquakes on active faults of the Chinese mainland. Acta Seismologica Sinica (in Chinese) , 1999, 21(6): 616-622.
[32] 梁海华, 候建军, 刘树文, 等. 中国构造应力场与大震复发周期关系的数值模拟. 地震地质 , 1999, 21(1): 52–56. Liang H H, Hou J J, Liu S W, et al. Tectonic stress field and large earthquake recurrence period in China. Seismology and Geology (in Chinese) , 1999, 21(1): 52-56.
[33] 陈运泰. 地震参数-数字地震学在地震预测中的应用. 北京: 地震出版社, 2003 . Chen Y T. Seismic Parameters-Applications of Numerical Seismology in Earthquake Prediction (in Chinese). Beijing: Seismological Press, 2003 .
[34] 汪一鹏, 马谨, 李传友. 南北地震带强震迁移特征及其与南亚地震带的联系. 地震地质 , 2007, 29(1): 1–14. Wang Y P, Ma J, Li C Y. The migration characteristics of strong earthquakes on the north-south seismic belt and its relation with the south Asia seismic belt. Seismology and Geology (in Chinese) , 2007, 29(1): 1-14.