Question

板内地震的非线性模型研究

Answer 1

李丽　张国民

（中国地震局分析预报中心，北京　100036）

石耀霖

（中国科技大学研究生院，北京　100039）

摘要　板内地震对人类造成的损失十分巨大。本文应用非线性动力学模型，以大陆构造块体为参照对象，对大陆地震的成组孕育和发生过程进行了模拟研究。非线性动力学模型的理论计算得到了一系列与中国大陆地震较类似的人工地震图像，如模型中的地震活动在时间上有轮回性、在地点上有条带性，并在不同轮回中有条带迁移的现象。模型中应力变化十分复杂，但仍能为预测地震提供概略性的信息。

关键词　板内地震　地球动力学　非线性模型

1　引言

地震分为板内地震和板间地震两种。大多数板内地震位于大陆地区，由于板内地震有时发生在人口密集地区，因而对人类社会造成了极大的破坏。板内地震从数目上仅占全球地震的15%，但造成的破坏占全球地震灾害的85%。大约三分之一的板内地震发生在中国。本世纪以来，中国承受的地震灾害占全球地震损失的50%，因此板内地震的研究对中国减轻地震灾害有十分重要的意义。

多年来，不少学者对地震孕育的动力学过程进行了多方面的探讨。其中，由Burridge和Knopoff提出的弹簧滑块模型最具代表性^[1]。继后，Byerlee^[2]、Dieterich^[4]和Rice^[9]等将诸如滑动弱化等复杂本构关系引入到一个自由度的单滑块系统中。1991年，朱元清和石耀霖用多个耦合非线性单元系统模拟了多个平行断层带的地震活动^[11]。近年来，石耀霖（1994）、张国民（1993）、耿鲁明（1993）等又进一步开展了应用非线性动力学模型模拟地震活动的研究^[6,10]。

本文采用一种非线性动力学模型对中国大陆板内地震的孕育和发生过程作了进一步的模拟研究。

2　动力学模型的建立

中国大陆位于欧亚板块的东南部，属欧亚板块内的中国－东南亚亚板块，是典型的板内区。这个具有整体性的大陆地块内部被活动断裂分割成各种不同级别的构造块体，构成了中国大陆地壳的层次块体结构图像。在研究地壳结构、地震构造和强震活动分布等的基础上，高维明等^[5]给出了中国大陆地震构造带的划分，他们在中国大陆地震区内划分出31个地震带和70多个孕震区。根据中国大陆地质构造的上述框架和大震带状分布图像，作为对强震图像的最初级的抽象，把我国大陆地震区抽象为含有若干地震带的大孕震系统。且每个地震带又含有若干个强震的孕震体。其结构如图1a所示。在将图1a转化为物理力学模型时需做如下假设：

图1　大陆地震构造抽象框图（a）及6×8模型图（b）

（1）由于大震的孕育一般可能需要数百至数千年乃至更长的时间，因此在整个孕震过程中，岩石的应力、应变积累、在一定程度上表现为流变体的积累过程。而地震发生则是一个瞬间的快速失稳破裂过程，其间发生震源破裂、岩体错动、应力下降等效应，并对周围岩石介质产生瞬间作用，故可将其看成是一个不连续、非线性的过程。为描述震源的这种性质，本文将震源简化成Maxwell体和刚性滑块串联的简单元件。文中称这些元件为基本元件。

（2）将整个大陆构造块体简化为含若干地震活动带的统一的孕震系统，每一个地震活动带上又可有若干孕震区。在物理模型中将地震活动带简化成若干基本元件的并列排列，而且若干地震带组合来模拟孕震系统。

（3）在相邻两排的基本元件之间连接一个耦合元件，来描述地震活动带之间和带内各段之间的相互耦合作用。耦合元件由弹簧和阻尼器构成。

其物理力学模型如图1b。考虑到现有运算设备的条件和运算的能力，本文将用6×8个弹簧－滑块－阻尼器组成的模型模拟孕震大系统。亦即在孕震系统中，含有6个地震带，每个地震带又含有8个孕震体。

须要说明的是，图1所给出的模型只是对中国大陆板块构造的一个极其初步的抽象和简化。中国大陆作为一个相对完整的构造块体，内含一系列次级块体及边界带（地震带），在印度板块和太平洋板块的动力边界作用下，形成一个统一的孕震大系统。系统内含有多个地震带，每带又含有多个震源体。这就是用图1模拟的基础。至于中国大陆构造的不均一性及其孕震系统的三维特性等问题，尚待今后更深入的研究。

3　动力学方程组及其求解简述

对图1b所示的力学系统的求解已由石耀霖等给出。震源区长期孕震阶段与发震时快速破裂阶段的力学特性截然不同，因此对每个震源体均从孕震阶段（即应力应变的积累阶段）和发震阶段（即滑块瞬间滑动阶段）分别进行力学模拟和计算。

3.1　应力积累阶段（孕震阶段）

设模型由m×n（m排n列）个基本元件和（m-1）×（n-1）个耦合元件组成。在孕震阶段每个基本元件和耦合元件均满足Maxwell体应力应变关系。假设：①拉应力为正，压应力为负；顺时针剪切力为正，逆时针剪切力为负；②串联元件的总应变率等于各个元件应变率之和；③在定常速率作用下，模型中各节点（连接周围元件的点）保持平衡，即与第i个节点相连的各元件应力之和应为零。

根据约束条件，当给定初始条件后，就可得到各单元上的应力值。

3.2　滑块错动阶段（地震发生）

构造地震是一种力学失稳现象，可将其看成是由任意微小的渐进变化的区域应力或位移而引起的断面上有限的突然的动力滑动。Burridge（1973）等曾提出描述这种现象的有限应力强度极限及其失稳准则。

模拟中根据该准则作一定的简化^[7]。设定串联元件的总应变调整量等于各元件应变调整量之和，并且各元件应力增量应满足节点平衡，即与第i个节点相连的各元件应力增量之和应为零，且系统总应变在滑动瞬间不产生变化，即：

第30届国际地质大会论文集　第5卷　现代岩石圈运动　地震地质

根据上述方法及设定条件，可列出一个线性方程组，对其求解，即可得到各元件应力调整值。实际计算中，给定应力初始值，可求解各元件的应力随时间的变化。当元件上应力值满足滑块滑动条件时，求解线性方程组得到各元件上应力调整后的应力值，并以其为新的应力初始值重新计算。整个模拟过程即为以上两步骤的不断重复。

4　理论结果及其与实际地震活动的初步对比

应用上述方程组的求解结果，并在设定各种参数的取值之后，即可计算出6×8模型的理论结果。在参数选取时，考虑到地壳岩石介质的差异不大，同时为计算简便，因此本文对各基本元件的弹性模量和粘性系数取相同的值。重要的是各基本元件的静摩擦强度取值不同，代表各震源区介质破裂强度的差异。而动摩擦强度与静摩擦强度之比为1∶1.25。边界条件取为定常应变速率。

图2为邱竞男^[8]等给出的我国近百年来强震活动的4个轮回图。从图中可见，每个轮回包含地震活动平静和活跃两种状态。其能量释放和应力变化曲线都与理论结果吻合。

张国民等指出^[10]，中国大陆的每个活跃期中的强震频度都随时间呈指数变化（图3）。对本世纪以来中国大陆地区地震时空分布的研究表明，不同的活跃期内，主体活动地区各不相同。本模型正是针对中国大陆地震分布的这种时空不均匀性特征而设计和计算的。

图2　我国近百年来强震活动的4个轮回图^[8]

图3　本世纪我国大陆强震活动的指数型分析（据张国民等^[10]）

4.1　模型地震特征

图4给出模型的计算曲线。时间取值范围为：t=980至t=1160。这时段中的地震活动大体上可以分为4个周期，每个周期包括一个应力积累、地震孕育的平静阶段和一个地震连发的活跃阶段。这4个周期为：t_Ⅰ∈[980,1040];t_Ⅱ∈[1040,1090];t_Ⅲ∈[1090,1120];t_Ⅵ∈[1120,1160]。图4与图2和图3给出的中国大陆地震活动情况十分类似。

图4　[980,1160]时段模型系统的平均应力图（a）、平均应变能量图（b）、累积释放能量图（c）、累积频度图（d）和M——t图（e）a、b、c图纵坐标为相对数值

图5是上述四个轮回的每一轮回中地震活动总图像，共分四组。6×8个网格代表48个元件，图框左边和上边的数字标定了各方块代表的元件行与列号。每组左图表示地震发生的顺序，右图表示地震的分布与震级大小。由图可见，每个轮回中地震的主体活动带都不同，而且随着时间的进程，每个轮回中的主体活动带显示一定的迁移趋向。这种迁移有一定的规律，大致是由南向北，再由北向南。在I轮回中，1带→6带；在Ⅱ轮回中，5带→1带；在Ⅲ轮回中，3带→2带→5带；在Ⅳ轮回中，6带→2带。从上述模型计算结果可看到，每个活跃期有一个地震活动主体地区。而不同活跃期中，主体活动地区呈迁移特点。这种模型结果可用以理解中国大陆地区地震轮回活动中主体活动地区的存在及其特性转移。

图5　由模型计算得到的各轮回地震空间分布

4.2　模型的应力图像

图6给出了模型在t∈[980,1050]时段内的应力随时间的变化图像。本文计算结果图中的时间是相对的，即以计算开始的时刻为零时刻，时间迭代步长为5，也就是当没有地震发生时每隔5个时间点进行一次迭代，当有地震发生时，从地震发生的时刻起，每5个时间点作一次迭代。图6中的[980,1050]时段就是计算了980个相对时间点以后的图像。图中的数字是相对应力值，代表每个孕震区的应力水平，每个元件在某时刻的相对应力值由以下公式给出：

第30届国际地质大会论文集　第5卷　现代岩石圈运动　地震地质

式中：σ_i,j指第i,j元件上的应力值；

为元件的动摩擦强度；

为元件的静摩擦强度。当且仅当V_i,j=10时，该元件破裂。图6方框中“－”表示该元件上的应力值低于动摩擦强度；“9＋”表示元件上的应力比9大接近于10，但还未破裂；“＊”表示在此时间点上元件破裂。从图6中可以看到，当t=982时，整体应力水平较低，但由于元件本身及系统的不均匀性，在经历一定的时间以后，一系列的破裂活动给出了整个系统的应力不均匀图像。应力经过15个时间点以后，由于在定常边界速率下系统整体水平逐渐增长，到t=997时，E3,8元件的应力值超过9,E6,6的应力值为9。当t=1001时E3,8首先破裂，并导致E3,7应力升高而发生破裂，这两次破裂导致E6,6应力下降。t∈[1001,1011]时段中，系统经过震后应力调整和在边界动力作用下应力的连续积累，使E1,5、E2,2、E6,6应力超过9。当t=1014时，E1,5、E2,2破裂，E6,6仍保持“9＋”状态。t=1019,E6,6破裂。由于系统在定常应变速率的边界动力作用下，其应力、应变和应变能的积累是持续不断、始终不停地进行的。所以到t=1035前，整个系统积累应力和应变能已达到很高的水平，因而在t=1035时，连续发生了5次破裂，这是这一时段中地震活动的高潮。由于t=1035前后连续发生了5次地震的高潮活动，从而使整个系统的应力、应变和应变能大量释放。因此，在t=1035之后，系统应力水平迅速下降到较低水平状态，使得t=1040后系统进入到长时间的应力积累过程，然后又开始下一轮的地震活动。模型的计算结果显示伴随着某元件的破裂，系统整体的应力水平都有所下降，但与破裂元件在同一地震带上的元件应力都有所升高，离破裂元件较近的元件应力升高较大。其他带上元件的应力都有所下降，离破裂元件越近的元件上应力值下降越多。另外，并不是每当元件的应力值升到9接近10时，元件都会破裂，例如t=1035时，E5,8和E6,8都达到“9＋”，但由于E5,8先破裂，导致E5,7应力升高而随之破裂（t=1039），这样使与E5,8同一排的E6,8应力大幅度减少，因而E6,8在经历了一个临近破裂的过程以后，应力下降，在此时段中没有发生破裂。由此看来，6×8系统的应力图像很复杂。

每当系统有元件破裂，系统的总平均应力值就有所下降，然后系统应力进一步积累，直到一定程度以后，又有元件破裂。模拟结果表明，孕震系统在边界动力作用下经历应力积累、大释放、再积累、再大释放的过程，相应的地震活动呈现为平静、高潮，又平静、又高潮的似周期或轮回活动图像。

5　讨论

在模型计算中，可以看到，在每个轮回中，各地震带的总体应力水平是相同的，只是有些带上应力分布比较均匀，因而不易发生大规模的破裂。而在另外一些条带上，应力集中地分布于几个孕震区中，从而使得这个条带在这一轮回中比较活跃。也就是说，比较活跃的地震带上总应力并不比不活跃的地震带上应力高，之所以某些条带地震比较活跃是因为在这些条带上应力分布不均匀。

考察图6发现，在一个轮回的地震活动开始时的应力图像对这一轮回中破裂发生的地点有较好的预示作用。如t=1011中应力超过9的三个元件在这一轮回中都发生了破裂。某一时刻的应力图像也可以预示下一时刻的地震发生地点，如t=997时E3,8的应力超过9，在t=1001时这个元件就破裂了。但应力图像又是十分复杂的，某个元件应力达到或超过9，并不一定就会马上破裂，甚至在这一轮回中都不会破裂。也就是说，虽然发震前应力场有较明确的显示，但这种显示与必然发震不是一一对应的。

图6　系统在980～1050时段中的应力图像

致谢　本文是应用非线性动力学模型对大陆地震孕育发生过程的初步研究。工作中得到耿鲁明、刘杰的热情帮助，特此感谢。

参考文献

[1]R.Burridge,L.Knopoff.Modeland Theoretical Seismicity.Bull.Seismol.Soc.Am.,1967,57:341～371.

[2]D.J.Byerlee.Friction of Rocks.Pure Appl.Geophy.,1978,116:615～626.

[3]Chen Huikai.Linear network and system.Translated by Wang Zhaoming,et al..Beijing:Electricity Industry Press,1988.

[4]J.Dieterich.Constitutive properties of faults with simulated gouge.Mechanical behavior of crustal rocks.A.G.U.Monograph,1981,24:103～120.

[5]Gao Weiming.et al..The boundary of seismologically risk regions and active tectonics.Earthquake 1996,（2）.

[6]Geng Luming,et al..Study on the relationship between the field and the source during the preparation and the occurrence of earthquakes.Earthquake Research in China,1994,9（4）:310～319.

[7]Li Li,Zhang Guomin and Shi Yaolin.A dynamic model for continental earthquakes.Seismology and Geology,1996,18（Suppl.Jun.）:9～20.

[8]Qiu Jingnan,et al..Discussion on the stages of seismicity in this century in China.Earthquake1986,（6）,41～47.

[9]J.R.Rice,et al..Dynamic motion of a single degree of freedom system following a rate and state dependent friction law.J.G.R.,1986,91:521～530.

[10]Zhang Guomin,et al..Computer models for the cycle activity of earthquakes in China mainland.Earthquake Research in China,1993,9（1）:20～30.

[11]Zhu Yuanqing.et al..Dynamically nonlinear models in the research on seismicity.A.G.S.1991,34（1）:20～31.