数值模拟方法在地震预测研究中应用的初步探讨(Ⅰ)

邓志辉 宋 键 孙君秀 陶京玲 胡勐乾马晓静 姜 辉 李 红

(中国地震局地质研究所，北京 100029)

数值模拟方法在地震预测研究中应用的初步探讨(Ⅰ)

邓志辉 宋 键 孙君秀 陶京玲 胡勐乾马晓静 姜 辉 李 红

(中国地震局地质研究所，北京 100029)

在地震危险性分析和预测研究中，自20世纪80年代初开始引入数值模拟方法以来，随着科学技术的不断进步，经历了由二维到三维，由线性到非线性，由弹性到粘弹性，由单场到多场耦合分析的改进和发展。地震的孕育和发生是复杂的物理过程，地震前的异常表现更是各种各样，但地震前的应变能量积累是地震发生的必要条件。地震预测分析必须首先考虑应变能量积累的状态。由于地球内部的难入性，直接测量震源深处的应力应变是很困难的事情，利用数值分析方法，建立地壳上地幔三维动力学模型，模拟岩层变形过程，是当前研究地壳能量转移、积累最有效的方法之一。

地震预测 数值模拟 等效应力 应变能密度

0 引言

科学的发展常常经历由定性的现象描述到定量的数学求解。数学问题的求解包括解析分析和数值模拟两种方法。解析分析是用数学分析的方法(比如微分、积分、特殊方程等)，对实际情况列出方程，用解析的方式求出函数解。数值模拟则是用计算机来做实验，通过数值计算和图像显示的方法，达到对工程问题和物理问题乃至自然界各类问题研究的目的。计算机数值模拟是一项综合应用技术，它已经与理论分析、试验研究成为科学技术探索研究的3个相互依存、不可缺少的手段。正如美国著名数学家拉克斯所说“科学计算是关系到国家安全、经济发展和科技进步的关键性环节，是事关国家命脉的大事”。

随着计算机性能和软件功能的不断提高，特别是高性能并行计算技术的进步，数值模拟已经成为解决复杂问题的重要手段。地质构造变形是复杂的地球科学问题，虽然经过一代一代的科学家的努力已经取得了重大的进展，但仍然无法得到各个地点不同时间的解析解，越来越多的学者认识到数值模拟是分析其时空分布规律的有效方法。地质构造变形数值模拟技术是从工程模拟技术移植并发展起来的，它们有很多相通之处，又有很大的区别。工程上的研究对象多为非常具体、规则、均匀和精确的模型，模拟结果要反映模型的细节，而地质构造变形数值模拟的对象是几何模型非常不规则，物理性质非常不均匀，多场耦合非常不清楚，边界条件非常不确定的地质体，所以对计算机性能和软件功能要求更高，对模拟研究人员要求知识面更广。

在地震危险性分析和预测研究中，自20世纪80年代初开始引入数值模拟方法(汪素云等，1980;王仁等，1980;罗焕炎等，1982;England et al.，1982;梅世蓉，1989)以来，随着科学技术的不断进步，经历了由二维到三维，由线型到非线性，由弹性到粘弹性，由单场到多场耦合分析的改进和发展(张东宁等，1999;刘洁等，1999;Parsons，2002;刘峡等，2006;陈连旺等，2008;胡勐乾等，2010)。迄今这一方法，尤其是有限元方法，已成为该领域中基本的研究手段之一。

总体来说，以往的工作已经取得了许多有意义的成果，特别是在模拟地壳应力应变场的时空分布规律上，成果更为丰富。本文将在前人工作的基础上，结合当前的地震形势和作者近期完成的一些研究结果，试图对数值模拟方法在地震预测中一些可能的应用进行举例分析。

1 地震能量积累分析

地震的孕育和发生是复杂的物理过程，地震前的异常表现更是各种各样，但地震前的应变能量积累是地震发生的必要条件。地震预测分析必须事先考虑应变能量积累的状态。由于地球内部的难入性，直接测量震源深处的应力应变是很困难的事情。利用数值分析方法，建立地壳上地幔三维动力学模型，模拟岩层变形过程，是当前研究地壳能量转移、积累最有效的方法之一。这里以青藏高原东构造结周边地区现今构造变形数值模拟分析为例，研究其应变能量分布的特征。

研究区范围取90°～100°E，25°～32°N，区内主要发育11条断裂构造。以这些断裂为边界，可把研究区划分为11个块体 (图1)。

在深度方向，将模型分为3层:上地壳、下地壳和岩石圈上地幔，总共120km。根据前人研究成果，各层深度略有不同，综合分析后选取了各层深度参数(表1)。单元类型采用的是Ansys软件Solid 186类型，地块内网格尺寸≤10km，断层带网格尺寸≤3km。模型中一共划分了130,574个单元，579,383个节点。

为了更好地模拟青藏高原脆性上地壳和柔性下地壳的变形，模型的上地壳采用弹性模型，下地壳及其以下选用的是粘弹性模型。地球物理观测和实验研究表明，青藏高原地壳厚度相当于正常地壳厚度的2倍，而且壳幔结构非常复杂，流变结构呈现分层性，具有脆性的上地壳和非常柔软的下地壳，多处都存在不均匀分布熔融层。很多学者对柔性下地壳进行了研究，结果表明青藏高原的下地壳黏度较周边的印度板块、塔里木盆地和华南地块都低。模型中弹性模量、泊松比的选取就是根据前人(赵文津等，2001;王椿镛等，2008;Denghai，2010)地震波接收函数研究中提供的P波、S波进行反演计算得到的(表1)。

对断层带采取弱化的方式，杨氏模量取值为两侧地块平均值的1/10。

本研究的模型边界条件采用位移约束，主要采用GPS观测块体位移速率数据，其来源包括以下3部分:一是“中国地壳运动观测网络工程”2007年的观测数据;二是国家自然科学基金项目“喜马拉雅东构造结周边地区主要断裂现今运动的GPS观测研究”新的观测数据，三是国内外有关文献。

由于南部边界主要在印度和缅甸境内，GPS观测点较少，目前仅有印度西隆的shill点和“中国地壳运动观测网络工程”缅甸密支那基准站可用，结合Paul等(2001)、Malaimani等(2008)、Larson等(1999)的研究，考虑到模型中主边界断裂简化为直立断层，扣除掉一些边界断裂俯冲作用吸收的1/3位移，作为南边界位移的约束速度。模型的顶面自由，底面垂向固定，水平向自由。由于目前对地表以下不同深度层位的运动速度变化不清楚，模型的侧面垂向自由，水平向将地表到模型底部取相同的位移约束条件。

由于地壳在长期的重力作用下，岩体已基本处于重力均衡状态，所使用的模型参数和边界位移数据也是在目前条件下的地球物理探测和测量结果。因此，本文在研究中没有考虑重力的作用。

本研究的计算步长取1a，共进行了1 000a的模拟分析，得到各节点的位移大小、方向、位移场和等效应力分布。模拟结果总体上与已有的GPS观测结果基本一致，青藏高原地壳水平位移场围绕东构造结顺时针转动。图2为研究区等效应力分布图，从图2中可见，应力集中区主要分布在东构造结周边的块体边界断裂带上，特别是嘉黎断裂东南段、墨脱断裂、阿帕龙断裂、印度-缅甸俯冲带以及川滇交界地区，应该注意这些地区未来发生强震的危险性。

图2 东构造结周边地区等效应力分布图Fig.2 Equivalent stress distribution around the Eastern Himalayan Syntaxis.

2 地震能量转移分析

在相互关联的构造区域内，块体与块体之间，断层与断层之间是相互作用的，一个区域发生地震，岩体快速位移和变形，必然引起相关联块体和断裂的位移和变形，应变能也将发生相应的转移。应用数值模拟分析的方法将可以帮助人们预测应变能可能的转移过程，为地震跟踪和强化监测提供参考靶区。这里以1989年大同5.8级地震和张北6.2级地震为例说明地震能量转移分析方法。

研究区包括首都圈地区，其空间范围为38.7°～41°N，112°～118°E，垂向深度50km。模型由上地壳浅层地壳、上地壳多震层、中地壳软弱层和下地壳及其与壳幔过渡带组成(嘉世旭等，2005)。其中，浅层地壳，基本上按照沉积层底面为界设置;中地壳软弱层即非连续分布的低速高导层。模型中各层的厚度及物理力学参数见表2。

分析系统采用ANSYS软件，模型采用20节点等参立体单元，遵循Maxwell粘弹性本构关系，并按线性粘弹性进行处理。将有强烈活动的断裂作为非线性接触单元处理，其他断裂作为连续介质处理。

模型的最大压应力主轴方向为NEE-SWW，设为9MPa，最小压应力主轴为NNW-SSE，设为3MPa，只考虑水平方向作用，未考虑垂直方向的作用。模型底面受垂向约束，可将其理解为岩石圈重力均衡的中性面，即相当于此面上、下呈镜像对称状，作为对岩石圈实际状况的粗略近似。

模拟从发生强地震的周期入手，模拟各个时期地壳运动和应力场的演化。围绕首都圈地区潜在地震危险性预测，计算时间从大同地震前开始，共计算了3个时间段的应力场和应变场，即:大同地震前、大同地震后和张北地震前。本研究采用降低岩石部分刚度的方法模拟地震的发生。模拟发现2次地震前都有多震层应变能密度显著增强的异常特征。

图3展示了大同地震前多震层应变能密度的分布等值线，从图3中可见，大同地区出现了高应变能密度异常区，应变能大的地区有利于育孕强震。1989年就在高应变能密度异常区发生了5.8级地震。

图3 大同地震前多震层应变能密度的分布(单位:J/m3)Fig.3 Strain energy density distribution before 1989 Datong earthquake(MS 5.8).

图4 是大同地震后多震层应变能密度分布等值线，图4中显示出大同地震后，由于能量释放，应力场进行了重新调整。除大同及其周围地区应变能下降外，其他地区应变能密度变化不大。

此后，经过9年的应力调整，应变能密度分布出现了新的变化。图5是张北地震前多震层应变能密度分布等值线，从图5中可见，在张北地震前，震中附近地区应变能密度显著增加，说明该地区在失稳前有能量的积累和强化。

图4 大同地震后多震层应变能密度的分布等值线(单位:J/m3)Fig.4 Strain energy density distribution after 1989 Datong earthquake(MS 5.8).

图5 张北地震前多震层应变能密度的分布等值线(单位:J/m3)Fig.5 Strain energy density distribution before 1998 Zhangbei earthquake(MS 6.2).

3 结论与讨论

地震的孕育和发生是复杂的物理过程，地震前的异常表现更是各种各样，但地震前的应变能量积累是地震发生的必要条件。地震预测分析必须首先考虑应变能量积累的状态。由于地球内部的难入性，直接测量震源深处的应力应变是很困难的事情，利用数值分析方法，建立地壳上地幔三维动力学模型，模拟岩层变形过程，是当前研究地壳能量转移、积累最有效的方法之一。

青藏高原东构造结附近地区现今构造变形数值模拟结果总体上与已有的GPS观测结果基本一致，青藏高原地壳水平位移场围绕东构造结顺时针转动。现今有效应力集中区主要分布在东构造结周边的块体边界断裂带上，特别是嘉黎断裂东南段、墨脱断裂、阿帕龙断裂、印度-缅甸俯冲带以及川滇交界地区，应该注意这些地区未来发生强震的危险性。

在相互关联的构造区域内，块体与块体之间，断层与断层之间是相互作用的，一个区域发生地震，岩体快速位移和变形，必然引起相关联块体和断裂的位移和变形，应变能也将发生相应的转移。1989—1998年首都圈构造变形数值模拟结果清晰地显示应变能密度高值区从1989年大同5.8级地震震中区转移到1998年张北6.2级地震的震中区。说明应用数值模拟分析的方法将可以帮助人们预测应变能可能的转移过程，为地震跟踪和强化监测提供参考靶区。

陈连旺，张培震，陆远忠，等.2008.川滇地区强震序列库仑破裂应力加卸载效应的数值模拟［J］.地球物理学报，51(5):1411—1421.

CHEN Lian-wang，ZHANG Pei-zhen，LU Yuan-zhong，et al.2008.Numerical simulation of loading/unloading effect on Coulomb failure stress among strong earthquakes in Sichuan-Yunnan area ［J］.Chinese J Geophys，51(5):1411—1421(in Chinese).

胡勐乾，邓志辉，陆远忠.2010.有限元数值模拟方法在华北地区地震地质研究中的应用进展［J］.地震地质，32(1):162—173.

HU Meng-qian，DENG Zhi-hui，LU Yuan-zhong.2010.Advances in application of finite elementnumerical simulation to seismogeology in North China［J］.Seismology and Geology，32(1):162—173(in Chinese).

嘉世旭，张先康.2005.华北不同构造块体地壳结构及其对比研究［J］.地球物理学报，48(3):611—620.

JIA Shi-xu，ZHANG Xian-kang.2005.Crustal structure and comparison of different tectonic blocks in North China［J］.Chinese JGeophys，48(3):611—620(in Chinese).

刘洁，宋惠珍.1999.用数值模拟方法评估华北北部地震危险性［J］.地震地质，21(3):221—228.

LIU Jie，SONG Hui-zhen.1999.Estimation of earthquake risk for northern North China by numericalmodeling ［J］.Seismology and Geology，21(3):221—228(in Chinese).

刘峡，傅容珊，杨国华，等.2006.用GPS资料研究华北地区形变场和构造应力场［J］.大地测量与地球动力学，26(3):33—39.

LIU Xia，FU Rong-shan，YANG Guo-hua et al.2006.2D deformation field and tectonic stress field constrained by GPS observations in North China［J］.Journal of Geodesy and Geodynamics，26(3):33—39(in Chinese).

罗焕炎，徐煜坚，宋惠珍，等.1982.青藏高原近代隆起原因及其与地震关系的有限元分析［J］.地震地质，4(1):31—37.

LUO Huan-yan，XU Yu-jian，SONG Hui-zhen，etal.1982.Finite elementanalysis for recentQinghai-Xizang Plateau uplifting and its relation to seismicity［J］.Seismology and Geology，4(1):31—37(in Chinese).

梅世蓉，梁北援.1989.唐山地震孕震过程的数值模拟［J］.中国地震，5(3):9—17.

MEIShi-rong LIANG Bei-yuan.1989.A mathematic simulation for seismogenic process of the Tangshan earthquake［J］.Earthquake Research in China，5(3):9—17(in Chinese).

王椿镛，楼海，吕智勇，等.2008.青藏高原东部地壳上地幔S波速度结构:下地壳流的深部环境［J］.中国科学(D 辑)，38(1):22—32.

WANG Chun-yong，LOU Hai，LÜ Zhi-yong，et al.2008.S-wave crustal and upper mantle's velocity structure in the eastern Tibetan Plateau-Deep environment of lower crustal flow ［J］.Science in China(Ser D)，38(1):22—32(in Chinese).

王仁，何国琦，殷有泉.1980.华北地区地震迁移规律的数学模拟［J］.地震学报，2(1):32—42.

WANG Ren，HE Guo-qi，YIN You-quan.1980.A mathematical simulation for the pattern of seismic transference in North China［J］.Acta Seismologica Sinica，2(1):32—42(in Chinese).

汪素云，陈培善.1980.中国及其邻区现代构造应力场的数值模拟［J］.地球物理学报，23(1):35—45.

WANG Su-yun，CHEN Pei-shan.1980.The numerical simulation of contemporary tectonic stress field in China and its adjacent area［J］.Acta Geophysica Sinica，23(1):35—45(in Chinese).

张东宁，许忠淮.1999.中国大陆岩石层动力学数值模型的边界条件［J］.地震学报，21(2):133—139.

ZHANG Dong-ning，XU Zhong-huai.1999.Boundary conditions of the dynamic numericalmodel for the Chinesemainland lithosphere［J］.Acta Seismologica Sinica，21(2):133—139(in Chinese).

赵文津及INDEPTH项目组.2001.喜马拉雅山及雅鲁藏布江缝合带深部结构与构造研究［M］.北京:地质出版社.

ZHAOWen-jin，INDEPTH Team.2001.Deep structure and tectonics of Himalaya Mountains and Yaluzangbu suture［M］.Geological Publishing House，Beijing(in Chinese).

Denghai Bai，Martyn JUnsworth，Max A Meju，et al.2010.Crustal deformation of the eastern Tibetan plateau revealed bymagnetotelluric imaging［J］.Nature Geoscience，1—5.doi:10.1038.

England P，McKenzie D.1982.A thin viscous sheetmodel for continental deformation［J］.Geophys JR Astr Soc，70:295—321.

Larson K M，Burgmann R，Bilham R，etal.1999.Kinematics of the India-Eurasia collision zone from GPSmeasurements［J］.JGeophys Res，104:1077—1093.

Malaimani E C，Ravi Kumar N，Akilan A，et al.2008.GPS ～ geodesy with GNSS receivers for Indian plate kinematics'studieswith the recent plate velocities estimated from GNSS data［J］.J Ind Geophys，12(3):109—114.

Parsons T.2002.Post-1906 stress recovery of the San Andreas Fault system calculated from three-dimensional finite element analysis［J］.Journal of geophysical research，107(8):1—13.

Paul J，Burgmann R，Gaur V K，et al.2001.Themotion and active deformation of India［J］.Geophysical Research letters，28(4):647—650.

PRELIM INARY STUDY ON APPLICATION OF NUMERICAL SIMULATION METHODS TO EARTHQUAKE PREDICTION RESEARCH(Ⅰ)

DENG Zhi-hui SONG Jian SUN Jun-xiu TAO Jing-ling HU Meng-qian MA Xiao-jing JIANG Hui LIHong
(Institute of Geology，China Earthquake Administration，Beijing 100029，China)

Earthquake preparation and occurrence is a complex physical process.Although the earthquake abnormalities are varied，the strain energy accumulation is requisite before an earthquake.Earthquake prediction analysismust consider the strain energy accumulation process.As hard to go into the Earth's interior，directmeasurement of stress and strain in deep focus is very difficulty.The use of numerical analysis，which constructs three-dimensional dynamicmodels of the crustand uppermantle to simulate the rock deformation process，is currently one of themost effectivemethods to study the crustal energy transfer and accumulation.

The simulation result of current crustal deformation is consistent with the existing GPS data around the Eastern Himalayan Syntaxis and its surrounding areas，in that the crustal horizontal displacement field of the eastern Tibetan Plateau rotates clockwise around the Eastern Himalayan Syntaxis.Current effective stress concentration areasmainly distribute along the block boundary fault belts around the Eastern Himalayan Syntaxis，especially along the southeast section of Jiali Fault，Moto Fault，Apalong Fault，India-Myanmar subduction zone and the Sichuan-Yunnan border region.It should be noted the risk of future strong earthquakes in these areas.

In the adjacent interconnected tectonic areas，the blocks and faults are interrelated and interacted each other.When an earthquake occurs in a region，the rapid displacement and deformation of rock will inevitably lead to displacement and deformation of the associated blocks and faults;strain energy will transfer from one region to others.The numerical simulation results of deformation process in the Capital area from 1989 to 1998 clearly show that the high strain energy concentration region shifted from Datong area where 1989 earthquake(MS5.8)occurred to Zhangbei area where 1998 earthquake happened.It illustrates that the application of numerical simulation analysis method may help us predict the possible strain energy transfer process，thus，providing the reference target regions for earthquakemonitoring.

earthquake prediction，numerical simulation，effective stress，strain energy density

P315.61

A

0253-4967(2011)03-0660-10

10.3969/j.issn.0253-4967.2011.03.015

2011-06-01收稿，2011-08-22改回。

中国地震局地质研究所基本科研业务专项(IGCEA1001)和国家自然科学基金(40841016、40372131、40702056)共同资助。

邓志辉，男，1962年出生，1992年在国家地震局地质研究所获大地构造物理专业博士学位，研究员，现主要研究方向为构造物理、地震预测方法研究，电话010-62009089，E-mail:deng 6789@163.com。