大华北地壳动力学参数及三维有限元模型的建立*

詹自敏 陈连旺 陆远忠

(中国地震局地壳应力研究所，北京 100085)

大华北地壳动力学参数及三维有限元模型的建立*

詹自敏 陈连旺 陆远忠

(中国地震局地壳应力研究所，北京 100085)

利用华北地区的断层数据、高程数据、波速结构数据、弹性常数以及多种流变机制参数，建立大华北地区三维有限元数值模型，并通过给每个单元赋属性的方法，实现物性参数的连续变化。以模型边界附近实测GPS数据插值作为边界位移条件，进行数值模拟。模拟结果表明，大华北区域位移场与实测位移场基本一致。

华北;有限元模型;数值模拟;GPS;地壳动力学

1 引言

大华北地区多次发生7级以上强震，给人民的生命财产安全造成了极大的威胁。

有限元的方法是在有限差分的基础上发展起来的数值模拟方法，运用有限元方法建立的地块实体模型，可以模拟地块深部应力场的变化和强震的孕育的过程。建立大华北地区三位实体模型，需要比较详尽的参数，而华北地区大量的监测数据，对于建立华北地区有限元模型提供了良好的条件。

早在1980年王仁等［1］便对华北地区进行了二维有限元数值模拟，采用弹性材料，并用初始应力场和应力释放的方法模拟了邢台、河间、渤海、海城、唐山5次地震的发生;2003年，文献［2，3］用非连续变形和有限元相结合的方法对华北地区进行了二维数值模拟，同样采用了弹性的本构关系;1998年，刘昌铨等［4］建立了华北的三维有限元弹性模型;1999年及其后，陈连旺等［5，6］分别利用弹性、黏弹性本构关系模拟了华北地区三维构造应力场以及强震应力场随时间的变化特征;2005年，刘勉等［7］建立了华北地区三维模型，分析了华北地区应变及应变能密度和震源机制;2007年，刘峡［8］在其博士论文中对华北地区进行了较详细的三维有限元数值模拟，建立了均匀分层和实际分层两个三维模型以及一个二维模型，使用了接触单元模拟断层，分析了华北地区的形变场和应力场，并探讨了其动力学机制;2010年，胡勐乾［9］也在其硕士学位论文中进行了华北地区的三维数值模拟，建立整个华北地区的分块模型，选取了多条弱化断裂，使用位移场作为边界条件，模拟华北地区速度场和应力场并进行了断层性质的分析。

然而以往的模型大多未能全面涵盖本区的地形、物性、实际断层几何特性，有的模型边界为简单矩形，未考虑块体构造特性。本文将在建立较完善的华北地区地壳动力学参数库的基础上，建立华北地区三维黏弹性数值模型，并用此模型模拟华北地区的形变场，探索地震孕育规律。

2 大华北地区地壳动力学参数

2.1 华北构造分区及主要的活动断裂

根据中国区域构造可以将华北地区划分为鄂尔多斯活动地块、华北平原活动地块和鲁东黄海活动地块，其中华北平原活动地块又可以分为3个三级活动地块：太行山次级活动地块、冀鲁次级活动地块和豫淮次级活动地块［10］。

华北地块的周围被深大断裂所围限，其内部也发育多条断裂，主要走向有北北东、北东、北西和近东西向4组［11］，其主要断裂如表1所示。

2.2 地形数据与深部介质物性参数

地形数据来源于 http：//srtm.csi.cgiar.org/。地形数据包括经度、纬度和高程值，高程的基准点为海平面。

地壳上地幔S波速度成果为文献［12，13］通过面波层析成像的方法得到的结果。依据理想各向同性介质中弹性常数剪切模量、泊松比和P波、S波速度之间的定量关系，计算得到的杨氏模量在三维立体空间内画散点图如图1所示。

岩石圈随深度不同具有不同的温压条件，岩石在高温高压下会出现一定的流变性质。2002年，臧绍先等［14］建立了一种华北岩石圈流变结构的初步模型，给出了3种不同岩石强度的计算公式(表2)。 s-1［15，16］。对于大陆岩石圈应变率一般取10-14～10-15s-1，然而深部岩石应变率的变化比较复杂，想要确定深部岩石的应变率比较困难，所以本文模型暂不考虑应变率随深度的变化，计算时使用全球平均结果10-15

图1 华北地区杨氏模量分布散点图Fig.1 Scattered diagram of Young modulus in North China

本文根据华北岩石组成的一般模型［17，18］，结合地下岩石的流变特性等性质［19］，再综合利用华北岩石圈三维流变结构的初步模型［14］，给出华北地下岩石分层以及岩石特性(表2)。

3 华北地区三维有限元模型的建立

3.1 模型边界与断层的选取

本研究区约为北纬30～43°、东经105～125°，北边界依燕山构造带和阴山一带而行，直至辽东;南边界以秦岭-大别山为界，直至黄海中部;西边界以贺兰山、六盘山为界;东边界位于124.5°E附近的海域中。模型边界全长5 320 km，东西向长度最大为1 660 km，南北向长度最大为1 350 km，深度为200 km。

断层在模型中的位置如图2所示。根据表1中各条断层的走向、倾向和倾角等数据，实现断层的倾斜，每一条断层用一个小体表示。

3.2 网格划分

网格为上层细网格、下层粗网格，模型主要采用六面体网格划分的方法，这种划分有两个难题：一是每一层网格要求的尺寸大小不一致，如何实现层与层之间六面体网格的过渡和转换;二是断层的倾斜，使得断层附近不可能建立完全直立的六面体单元，如何即保证网格是六面体，又实现断层的倾斜。

对于第一个问题，不同层之间不同尺寸的单元网格的连接，采用金字塔过渡单元来解决(图3 (a))。模型中，结合表3的分层，0～30 km为第一大层，此层又可分为4个小层，每层约为8 km;40～80 km为第二大层，100～200 km为第三大层，30～40 km和80～100 km为过渡层。

表1 华北地区主要断裂一览表Tab.1 Major faults in North China

表2 模型选用岩石分层及岩石性质参数列表Tab.2 Parameters of rocks in different layers of this model

关于第二个问题，在断层附近采取扫略方法实现六面体网格的划分，而在距离断层较远的地方，使用正常的六面体网格划分(图3(b)、(c))。

网格划分的整体效果如图3(a)所示。模型共有254 156个单元，234 056个节点。

图2 模型中断层的位置及名称示意图Fig.2 Location and name of selected faults in the model

图3 华北模型网格划分及倾斜断层示意图Fig.3 Meshing and the inclined fault of the model

3.3 地形起伏

为了在模型顶面实现地形的起伏，并考虑到地表起伏的高度相对模型的范围来说较小，一般不会出现夹角过大的单元，所以采用直接改变地表节点Z坐标值的方法，即将模型地表的每一个节点的坐标值提取出来，在地形参数库中运用插值的方法得到相应的高程值，然后改变地表节点的Z坐标，实现模型的地表高程起伏(图4，将模型地表的高程值放大50倍的效果图)。

3.4 物性参数

使用PRONY模型模拟黏弹性，PRNOY模型的核方程为：

图4 华北模型地表起伏示意图Fig.4 Surface topographic undulation of the model

其中G为剪切模量，K为体变模量，G∞和K∞为时间无穷大时剪切模量和体变模量的取值，而Gi和Ki则分别为第i个PRONY元件的剪切和体变模量。nG和nK为PRONY元件的个数(简单情况下取1)。t为时间和分别为剪切模量和体变模量的松弛时间。

模拟中用弱化带模拟断层，将断层的物性参数赋为周边块体的十分之一。

本模型与前人模型的重要区别在于将物性参数赋值到每个单元之上，实现了物性参数较为连续的变化，相对来说更接近实际。因此，本模型介质属性的组数等于单元个数，共有254 156组介质属性。由于目前波速结构结果分辨率的限制，导致一些单元的参数是相同的。

4 华北地区GPS位移场模拟与分析

为了验证模型内部材料属性是否能够反应真实情况，本文将中国地震局第一形变监测中心获得的GPS年平均位移数据［20］插值得到三维模型边界处的位移数值作为三维模型6个侧面的边界条件。由于深部运动速率与GPS速率之间的差异尚无明确结论，作为一种近似，本文所施加的侧面边界条件未考虑沿深度的变化。模型的底面在垂直方向固定不动，这近似于一种重力均衡的表现。在上述边界条件作用及考虑重力载荷下，对模型进行静态计算，并将计算得到的地表位移场与GPS实测位移场进行比较(图5，红色箭头表示华北地区GPS年位移场，每一个数据值均为一个观测点的数值，没有进行插值;蓝色箭头为数值模拟结果)。

为进一步分析华北块体内部的模拟情况，将华北地区分为西北、中北、东北、西南、中南、东南6块，分别对每个块区的模拟值与GPS观测值的区域平均结果进行比较(表3)。

图5 华北地区GPS实测位移场与模型模拟结果对照图Fig.5 Comparison between GPS observed displacement values and the simulated values

表3 GPS观测值与模型模拟值分区比较Tab.3 Comparsion between GPS observed values and the simulated values in different districts

从图5及表3可以看出，除河北中东部分地区差别较大以及山西和河南境内部分观测点方向不一致外，大部分地区位移场的大小和方向的分布与观测GPS位移场基本一致，其总体走势也基本相同，模型模拟出来的位移场与真实的GPS位移场分布具有很高的相似性，除东北区的角度相差较大外，绝大部分区域的数值大小和角度方向均相差不大。

5 结论

本文建立的大华北地区三维精细结构模型所模拟出来的华北地区位移场与GPS观测位移场符合较好，说明本文所建立的大华北地区参数数据来源可靠，具有真实的地理信息依据，能够真实地反应华北地区深部构造状况。

致谢 感谢李延兴、黄忠贤、马保起研究员在GPS观测数据、深部波速结构数据及地形高程数据等方面提供的帮助与指导。

1 王仁，等.华北地区地震迁移规律的数学模拟［J］.地震学报，1980，2(1)：32-42.

2 白武明，林邦慧，陈祖安.1976年唐山大震发生对华北地区各地块运动与变形影响的数值模拟研究［J］.中国科学(D辑)，2003，33(增刊)：99-107.

3 陈祖安，等.1966年以来华北地区一系列七级大震破裂过程的数值模拟［J］.地球物理学报，2003，46(3)：373-381.

4 刘昌铨，刘明军，嘉世旭.利用华北北部深部地球物理资料数值模拟地壳应力场［J］.地震学报，1998，20(3)：240-249.

5 陈连旺，等.华北地区断层运动与三维构造应力场的演化［J］.地震学报，2001，23(4)：349-361.

6 陈连旺，等.1966年邢台地震引起的华北地区应力场动态演化过程的三维黏弹性模拟［J］.地震学报，2001，23 (5)：480-491.

7 Liu Mian and Yang Youqing.Contrasting seismicity between the north China and south China blocks：Kinematics and geodynamics［J］.Geophysics Research Letters，2005，32，L12310，doi：1029/2005GL023048.

8 刘峡.华北地区现今地壳运动及形变动力学数值模拟［D］.中国科学技术大学，2007.

9 胡勐乾.并行计算三维数值模拟在华北地区现今构造变形分析中的应用研究［D］.中国地震局地质研究所，2010.

10 韩竹军，等.华北地区活动地块与强震活动［J］.中国科学(D辑)，2003，33(增刊)：108-118.

11 李铁朋，等.华北地区Ms≥6.5级地震震源断层参数的研究［J］.地球物理学进展，2007，22(1)：95-103.

12 黄忠贤，等.中国东部海域岩石圈结构面波层析成像［J］.地球物理学报，2009，52(3)：653-662.

13 Zhongxian Huang，et al.The lithosphere of North China Craton from surface wave tomography［J］.Earth and Planetary Science Letters，2009，288：164-173)

14 臧绍先，等.华北岩石圈三维流变结构的一种初步模型［J］.中国科学(D辑)，2002，32(7)：588-597.

15 魏荣强，臧绍先.大陆岩石圈流变结构研究进展及存在的问题［J］.地球物理学进展，2007，22(2)：359-364.

16 魏荣强，臧绍先.岩石圈流变结构的一种新的应变率约束［J］.地球物理学报，2004，47(6)：1 029-1 034.

17 吴宗絮，等.华北大陆地壳——上地幔岩石学结构与演化［J］.岩石矿物学杂志，1994，13(2)：106-115.

18 邢集善，刘建华，赵晋泉.华北板内深部构造［J］.山西地震，2002，4：3-12.

19 周永胜，何昌荣.地壳主要岩石流变参数及华北地壳流变性质研究［J］.地震地质，2003，25(1)：109-122.

20 张静华，等.用GPS测量结果研究华北现今构造形变场［J］.大地测量与地球动力学，2004，(3)：40-46.

CRUSTAL DYNAMIC PARAMETERS AND 3D FINITE ELEMENT MODEL OF NORTH CHINA

Zhan Zimin，Chen Lianwang and Lu Yuanzhong
(Institute of Crustal Dynamics，CEA，Beijing 100085)

Using the data of fault and terrain，seismic velocity structure and rarious rheological parameters of North China，to build a crustal dynamic parameters library and build a 3D finite element numerical simulation model of this region based on the library.Different from other models，an attribute is set for every element，so that the variation of attribute is continuous.By using the interpolated GPS data as displacement boundary conditions，the calculated displacement field is similar to the observational displacement field.This model has great significance for crustal dynamics study，for example，the research on the stress and strain fields of spatio-temporal evolution in North China，the characteristics of strong earthquakes activities，and the analysis of stress field after strong earthquakes and so on.

North China;finite element model;numerical simulation;GPS;crustal dynamics

1671-5942(2011)Supp.-0028-05

2011-01-14

中央级公益性科研院所基本科研业务专项重点项目(ZDJ2009-06);中央级公益性科研院所基本科研业务专项重大项目(ZDJ2007-01)

詹自敏，女，1987年生，硕士研究生，主要研究方向为地壳动力学数值模拟.E-mail：zhanzimin@gmail.com

P315.8

A