三维黏弹性数值模拟华北盆地地震空间分布与构造应力积累关系

柳 畅，石耀霖，郑 亮，朱伯靖

1 中国科学院大学，中国科学院计算地球动力学重点实验室，北京 100039

2 Laboratoire de Géologie，CNRS－UMR8538，Ecole Normale Supérieure，Paris 75231

1 引 言

华北盆地位于中国大陆北部，地处燕山以南，太行山以东，南以广饶—聊城断裂为界与鲁西隆起相邻，东与渤海湾相接［1］.华北克拉通在早前寒武纪克拉通化后的演化中遭受到强烈的破坏和再造，使得200km左右厚而冷的岩石圈转变为80～120km热的岩石圈，虽然华北克拉通的破坏时间尚存争议，但是克拉通破坏事件已是共识［2-5］.华北盆地裂陷区新生代构造活动强烈，在早第三纪经过强烈的拉张断陷，形成盆-岭构造；晚第三纪起，由于区域沉降而成为大型坳陷盆地.在裂陷区沉降时，周围地区相对抬升为山地.

尽管华北盆地远离大陆活动边界带及板块俯冲带，但历史以来华北盆地板内地震活动频繁（图1），破坏性地震时常发生［6-8］.自1966年至1976年的10年间，华北地区相继发生了1966年邢台MS7.2级，以及1967年河间MS6.3级地震，1969年渤海MS7.4级，和1976年唐山MS7.8级大地震，地震历史记录1697年三河—平谷曾经发生M8.0级大地震.邢台、唐山地震是我国大陆地区为数不多的两次强震群型地震.邢台地震序列包含两次不在一处的强主震（3月8日马兰MS6.8级，3月22日东汪MS7.2级地震）和30余次MS≥5.0地震；唐山地震序列包含3次强震（唐山市MS7.8级主震、同日滦县MS7.1级地震、同年11月15日宁河MS6.9级地震）和30多次 MS≥5.0强余震［8］.

大量地震勘探成果所揭示的华北盆地地壳结构与余震震源定位工作，表明华北地震空间分布特征为：

其一是，在地震区域分布上，华北盆地板内大地震基本发生在地壳的薄弱地带（Moho面上隆地带，图2），或者地壳厚度的急剧变化带，如唐山地震、邢台地震、河间地震［9-13］.

其二是，华北板内地震并不是遍及整个地壳范围内，而是多集中在地壳的一定深度范围内，主震一般在上地壳底部9～15km深度范围，余震多发生在上地壳与中地壳的5～25km深度范围内，在沉积层与下地壳中鲜见有余震发生［6，10，14-17］.

图1 研究区域 ）红色三角表示城市地名，白色圆点表示华北地震主震，蓝色实线表示下文（图7提到的纵剖面.Fig.1 Map of the North China basinRed triangles indicate city names.White dots indicate the epicenters of（ main sho）cks.Blue line indicates the cross section along several earthquake zonesin Fig.7.

邢台地震主震发生在37°6′N，114°18′E，主震震源深度为9km.余震优势层在10～22km深度范围，主要分布在中地壳的深度范围［6］.地震断层面走向N35°E，倾向N125°E，断层总长度50km.余震空间分布总体呈走向上的椭圆状分布，MS7.2级主震位于震区中部，其余MS6.0级以上地震以近等间距对称地排列在主震的南北两侧［18-19］.

唐山MS7.8级地震主震发生在39°24′N，118°30′E，震源深度为11km，最大余震栾县MS7.1地震深度为10km，宁河地震深度为17km.余震分布在5～22km 范围内，为上地壳与中地壳 深度范围［6，14-17］.陈祖安等［20］和张之立等［21］利用P波资料测定震源参数，其结果表明唐山地震初始破裂的走向为NE30°，但总体走向 NE49°，断层面倾向SE，倾角76°.断层破裂方式为不对称的双侧破裂，向NE传播70km，向SW传播45km，同震平均位错136cm，最大位错达173cm，在地表造成长8km、宽30m的地表裂缝带.张宏志等［22］利用现今唐山余震定位的结果认为唐山震区内存在明显的差异，并将其划分出3个区域：唐山MS7.8级主震所在的中区、滦县MS7.1级余震所在的东区和宁河MS6.9级余震所在的西区.中区内的地震条带走向从南向北以NNE向延伸20km，在主震震中处出现一个明显的转折，转向NE延伸约50km.东区内的地震呈“丁”字形的2个条带，SE向长约20km，NE向长约30km，此处发生最大的MS7.1级栾县余震.西区内的地震没有密集成带，没有明显优势方向的区域，长度30km.地震区自西南到东北总长140km，现今唐山余震与1976—1979年的余震震中分布情况大体一致.

图2 华北盆地Moho面深度图白色圆圈表示华北地震主震投影.Fig.2 Depth of Moho surface beneath the North China basin White dots indicate the epicenters of main shocks in the North China basin.

河间地震发生在38°29′N，116°28′E，震源深度为30km.地震活动在38°10′N—39°05′N，116°00′E—117°00′E的范围内，绝大部分余震多密集分布在主震 东 南 侧 的 38°20′N—38°30′N，116°25′E—116°44′E的近似一个椭圆形小区内，能量集中在该区域内释放，地震影响分布广，破坏小（http：／／www.csi.ac.cn／manage／html／4028861611c5c2b a0111c5c558b00001／zhenli／zhwei／html／zhenli084.htm）.

三河 M8.0级古地震发生在39°30′N，117°9′E.其主震震源深度通过地表等震线分布研究确定在15km［23］.该地震的地表破裂带展布在京东夏垫一带.该地表破裂带呈NE方向、以雁列方式展布，长度约l0km，活动方式为兼具右旋走滑的正倾滑活动.通过地表陡坎调查、探槽及浅钻研究，认为该次地震的同震最大垂直位移为3.16m［24］.

渤海MS7.4级地震主震发生在38°20′N，119°27′E，魏光兴等［25］通过研究渤海地震震源分布，认为主震震源深度为30km.而束沛镒等［26］利用远震P波波形反演渤海地震的震源参数认为震源深度在25km.渤海地震是典型的主-余震型地震，发生在郯-庐大断裂的西侧，震中区位于一系列北西向构造与郯-庐大断裂的交汇点.余震震中一般分布在37°55′N—38°30′N，119°15′E—119°55′E 范围内的椭圆，其长轴约为73km，短轴约为51km.

为解释华北盆地地震的孕震成因以及如上地震空间分布特征，历来已有诸多研究从各方面进行了定性或者定量的探讨.

近30年来，在华北盆地开展了大量的石油地震勘探、地震波折射、地震波广角反射剖面、层析成像和接收函数的研究工作，以探讨华北盆地地壳结构与地震发生成因之间的关系.结论认为华北地区地壳构造环境是强震的孕育基础，具体结论可归纳如下三点：

（1）震区高、低速体交错的包体现象是强震孕育的重要基础.如邢台、唐山震区.

（2）震区Moho界面的隆起和地壳力学性质不同的多层结构是强震发生的深部条件.如邢台、唐山、河间震区.

（3）震区上地壳地幔的上涌变形作用使能量在地壳大规模积累与释放.如邢台、唐山震区.

与此同时，不同学者对华北构造应力场的特征与演化进行了大量的计算，针对以上总结的三种推论，分别对地壳内高、低速体的存在，Moho面的隆起，和上地幔的上涌三种因素在华北地震系列孕震过程中的作用，进行了定量计算研究.

宋惠珍等［27］利用有限元方法，研究了唐山地震三维震源断层的应力场，认为唐山地震可能是水平和垂直力同时作用下产生的结果.郑治真等［28］曾用二维有限元方法研究了块体在水平挤压作用下唐山大地震的孕震过程.张东宁等［29］讨论了冀中坳陷西侧太行山稳定地壳的阻挡作用与滑脱构造发育的关系，模拟结果显示地幔热物质的侵人在冀中坳陷上地壳造成张应力区，随时间推移上地壳的拉张应力状态基本不变.王继存等［30］利用二维数值方法研究了唐山地震从孕震到发震的力学过程.梅世蓉等［31］基于二维Maxwell体黏弹性模型，进行了唐山地震孕震过程的数值模拟，并首次提出地震成因的“坚固体模型”，其结果认为唐山大地震的前兆现象与坚固体的软化有密切关系.冯德益等［32］利用华北地区岩石圈12层三维弹性非均匀模型，在模型四周施加力的边界条件，并在模型底部施加0～100MPa的上拱力进行了应力场计算，讨论了唐山地震孕育背景.其结果表明：（1）唐山地震震源位于华北地壳中的坚固体之内，该坚固体能在地壳内形成局部的应力集中；（2）华北岩石圈上部块体基底承受的地幔上拱力能在地壳内形成高度应力集中区，因而可能在唐山MS7.8级地震的孕育中起较大作用.尹京苑和梅世蓉［33］在弹性模型中同时考虑了地壳浅层沉积盆地，地壳内高、低速夹层结构，超壳深断裂及Moho面的局部上隆与上地幔热物质的侵人地壳这几种因素.其结果表明：（1）高速体单独存在可导致应力在高速体内的局部集中，并在其边缘和4个角区形成应力的高度集中，有利于强震的孕育发生；（2）Moho面局部隆起对壳内应力场的影响较深浅断层的影响规模更大.但是尹京苑和梅世蓉的结果又强调Moho面隆起造成的应力增强的密集带不在隆起区的正上方，而是在其斜坡上.孙武城等［34］曾设计了光弹模拟实验，测定在Moho面呈波浪形起伏的情况下，构造应力引起的地壳内应力场的分布状态.其结果发现，地壳中最大剪应力既不在波峰处，也不在波谷处，而是位于 Moho面起伏的“拐点”上方.肖兰喜等［35］将地壳介质视为Maxwell体，运用三维有限元方法，在模型四周施加力的挤压的边界条件下，探讨分层地壳结构中存在高速体、低速体和深部断裂时，地壳应力的集中及应力集中随时间演化的特征.并同时考虑了地幔以1cm／a速率上隆时，50年内地壳模型平均应力扰动值随时间的变化.从施加的应力边界条件上来看，肖兰喜等考虑的是Maxwell体的应力松驰的问题，并不一定与实际的应力积累过程相符.宋治平在弹性模型中考虑了包体模型应力集中问题.朱守彪等［1］利用有限元方法对华北盆地发震断层下方可能存在的速度异常体进行了力学分析，计算结果显示：无论速度异常区是低速异常还是高速异常，在水平构造力的作用下（无论是挤压还是引张力作用），应力及应变能密度都在空间某些区域集中积累.

以上工作从数值模拟角度对地震的孕震成因进行解释，基本支持从地壳结构所得到的三点认识，认为地壳内的高低速包裹体、Moho面的隆起和壳幔的相互作用在华北地震孕育过程起到重要的作用.并认为地震孕育的动力源存在于两方面.

（1）一部分研究认为华北地震的成因首先考虑板块之间互相挤压造成的水平压应力，如邢台、唐山强震这样的大陆地震，是由不同板块碰撞时的水平压力所引起的，高速体的存在，使高速体中产生应力集中现象，断层的存在也会使应力在断裂带产生应力集中.

（2）另一部分研究却对板块的水平运动是引起大陆强震的唯一力源表示怀疑，认为地震的动力来源并非主要来自地壳的上部，而应该来自上地幔的作用.在远离板块边界的唐山地区，仅用地壳块体之间的水平向相互作用无法解释唐山大地震的成因.同时GPS及大地水准测量等观测结果并不能支持唐山地区具有大的水平向运动.只有上地幔的变形运动才能提供在地壳内积累巨大地震能量的动力来源.

需要指出的是，根据现今对华北岩石圈各层位的流变结构的认识，结合其考虑的问题与研究对象，在以上一些计算之中模型所涉及的物质性质或者边界条件尚有些不足之处.比如在考虑长期地质作用下（水平板块构造作用与垂直方向的壳幔相互作用），地壳应力积累问题的时候，模型选择使用完全弹性体的本构关系，以及施加载力的边界条件或者固体位移的边界条件并不十分恰当.因此，比如，孙武城等［34］、尹京苑和梅世蓉［33］分别在光弹实验与理论模型计算的结果中得到起伏的分层介质界面波峰没有应力集中的现象，进行由此推论地壳的上隆处也不会出现应力集中现象的结论也还有待商榷.

本文结合实际地震勘探结果所揭示的华北盆地三维岩石圈结构，与温度计算所揭示的地壳分层流变特征，建立华北盆地岩石圈三维黏弹性有限元模型.综合华北盆地现今地表应变率，与余震震源机制与GPS所揭示的最大主压力方向，确定位移挤压边界条件，数值模拟华北盆地在长时期的水平构造挤压作用下地壳各层位应力积累过程.在此基础上探讨地壳结构尤其是地壳薄弱带（Moho面上隆处）在应力积累与地震孕育过程中的控制作用，并分析华北盆地地震空间分布与构造应力积累的关系.

2 研究区域地质构造特征

2.1 华北岩石圈分层结构

在过去30年中，大量的石油地震勘探、地震深反射剖面、层析成像以及接收函数的研究工作，揭示了华北盆地地壳结构［9，11-13，36-40］.地震测深资料分析表明，华北地壳可分为上、中、下三部分，Moho界面明显波浪起伏（图2），起伏的波峰在不同地区又各自不同，局部Moho面波浪起伏剧烈，梯度较大.地壳厚度总体趋势为西厚东薄，在西北张家口—怀来地区地壳厚度较大可达43km左右，而往东至唐山，地壳则减薄到28～32km.华北盆地震区均有Moho面上隆现象，如：邢台震区Moho面在震区局部隆起，最薄地壳厚度为31km［9-10，13，41］.唐山震区地壳上隆，厚度减薄为28～32km，但周边地壳厚度却急剧增加到40km，地壳厚度的最大差异达到10 km 左右［10，13，40，42-43］.河间震区 Moho面局部隆起，地壳厚度减薄为30km左右，周边地区地壳厚度在横向约60km范围内增厚到36km左右［13］.三河震区较周边地区 Moho面有较弱的2～3km上隆幅度［12-13］.

P波与面波层析成像与接收函数工作揭示了现今华北岩石圈西厚东薄的特征，厚度在80～120km之间［44-46］.

2.2 壳、幔物质流变结构

岩石圈热状态决定其力学性质，大陆岩石圈一般存在脆性的上地壳、柔性的下地壳和较强上地幔这种三明治式的分层流变结构.臧绍先等［47］利用华北地区的地震P波速度资料和大地热流资料建立了华北岩石圈的三维波速分布及温度分布.考虑了摩擦滑动、脆性破裂及蠕变三种主要的流变机制在岩石圈中的作用.计算了华北岩石圈流变强度及黏度的三维分布，认为上地壳上部为脆性区，中地壳可以大部分是以脆性破裂为主的脆性区，部分是以蠕变为主的延性区，而下地壳几乎均是以蠕变为主的延性区.安美建和石耀霖［48］用地震波速资料计算了中国大陆上地幔温度场，为地壳底部温度提供了约束，并对中国大陆岩石圈温度场进行了重新计算.石耀霖和曹建玲［49］根据以上的岩石圈温度结果，结合中国大陆地壳成分模型的研究成果，及基于GPS观测的地表应变速率，对华北地区及中国大陆其它地区岩石圈三维流变结构进行了估计［49］.同时，Li等［50-52］应用不同的岩石圈流变结构和会聚速率对苏鲁—大别以及喜马拉雅造山带的动力学进行了系统的研究.

3 三维黏弹性岩石圈模型及应力积累

3.1 有限元模型的建立

本文建立如图1所示区域范围内三维黏弹性Maxwell体岩石圈模型（见图3）.计算中的坐标系为X轴朝东、Y轴朝北、Z轴朝下.模型深度为从地面0至岩石圈底部100km深（岩石圈厚度取平均值而并未考虑横向厚度变化）.模型结构具体考虑所示区域地形高差（ETOPO30），并结合上述地震勘探成果中显示的地壳分层结构，Moho面的起伏，以及不同层位岩石介质的黏滞系数的差异.

计算所涉及的黏弹性介质的3个物质参数分别为：杨氏模量E，泊松比ν和黏滞系数η.岩石杨氏模量E 与泊松比ν依据波速计算［9，11-13，36-40］.同一层位的物质参数相同，一般更深的层位地震波速较大，岩石的杨氏模量亦增加，依据所选择物质参数，华北盆地岩石圈100km深度范围内介质杨氏模量随深度的变化曲线见图4.

结合臧绍先等［47］与石耀霖等［49］对华北岩石圈的流变参数的计算结果，对所建立的三维黏弹性Maxwell体岩石圈模型各层位赋予相应黏滞系数.相应的有限元模型（图3）介质的黏滞系数与介质松驰时间分别为：沉积层（编号1）3.2×1022Pa·s，30000年左右；上地壳（2）1.3×1023Pa·s，50000年左右；中地壳上层（3）2.0×1022Pa·s，8000年左右；中地壳下层（4）9.5×1021Pa·s，2800年；下地壳（5）3.0×1020Pa·s，60年；岩石圈上地幔（6）2.0×1021Pa·s，300年［49］.

针对以上黏弹性Maxwell体模型，在万年的时间尺度内考察岩石圈的应力积累过程，上地壳与中地壳上层将表现为脆性，中地壳下层、下地壳与岩石圈上地幔介质的流变性质则可以充分得到体现（取决于黏弹性介质的松驰时间）.由于Maxwell体具有这种可以自然处理脆-韧性转变的优点，本文有限元模型实际可以有效地描述弹性层与黏弹性层耦合条件下的应力积累过程.

100km的岩石圈分为6层，1层为沉积层，2层为上地壳，3层和4层分别为中地壳上层与下层，5层为下地壳，6层为岩石圈上地幔（图5）.用六面体单元对模型进行规则网格划分，横向上，单元网格大小统一为5km；深度向，地壳内单元网格大小为1km，而岩石圈上地幔内单元网格大小为5km.单元总数为896000个.

3.2 边界条件和初始条件

大量的华北地震系列震源机制与主压应变研究表明，华北地区构造应力场在各地区极为统一（包括震前与震后调整），最大主压应力主轴接近水平，在NEE—SWW方向，最小主压应力在NNW—SSE方向［53-55］.华北地区的应变测量认为华北应变率量级在0.5×10-9／a左右［54］.因此本文计算的边界条件为：在最大主压应力方向上施加挤压的位移边界条件.在侧面边界上，认为模型受到了来自于NEE75°方向上3mm／a的位移挤压，并将该位移边界条件线性分解推广到模型北侧面与东侧面边界作为水平速度约束条件；而南面与西面侧边界水平速度约束为0；且假定该条件从地表到100km深度保持一致，这未必与实际情况相符，但不妨作为一级近似；垂直方向位移保持自由.上表面为自由边界，即法向应力和剪应力均为零.对于底部边界，鉴于在本文研究的目标只考虑板块构造加载在华北岩石圈应力积累的作用，并不考虑该区域软流圈地幔-岩石圈相互作用，因此暂且将底面垂直方向速度约束为0，而水平方向自由.在黏弹性问题中，边界条件随时间的变化也是重要的问题，在目前的模拟中，我们初步假定边界位移速率不随时间变化.

图3 华北盆地岩石圈模型［9，11-13，36-40］，尺寸为888km×777km×100km（X×Y×Z）岩石圈中不同层位用数字标出：编号1为沉积层，编号2为上地壳，编号3为中地壳上层，编号4为中地壳下层，编号5为下地壳，编号6为岩石圈上地幔.Fig.3 Model of the lithosphere of North China basin［9，11-13，36-40］with the size of 888km×777km×100km（X×Y×Z）Different layers are marked by numbers.

图4 华北岩石圈物质杨氏模量随深度变化图［9，11-13，36-40］Fig.4 Diagram of Young′s modulus in the lithosphere of the North China basin［9，11-13，36-40］

图5 华北岩石圈各层位岩石介质黏滞系数在邢台、河间与唐山震区随深度变化图［49］Fig.5 Diagram of viscosity in the lithosphere of the North China basin［49］

初始条件是构造应力场模拟中最困难的问题，我们对深部三维应力分布和应力演变历史几乎没有定量的资料.在这种情况下，只能先假定初始应力为0，然后计算在定长的边界位移速率条件下应力的演变.虽然我们不可能模拟真实的应力状态，但可以了解在定长边界位移速率下应力增长率.应力增长率高的地方，未必一定是目前应力最高的地方，但如果初始应力类似，则较高应力增长率的地区则更有可能现今应力绝对值较大，本文主要就应力增长率与地震活动性的关系进行讨论.

4 应力积累计算结果

计算所使用的程序是利用ANSYS有限元软件黏弹性处理模块，程序计算的可靠性已经通过大量实例验证.本文采用1年为一个时间步长，根据时间步长逐步加载边界位移约束.经过了20000个时间步，即20000年的加载作用，计算华北盆地岩石圈模型各层位的应力积累过程与应力积累分布特征.本文的主要目标是计算华北岩石圈各层位对构造加载的响应，在计算过程中只考虑在边界条件作用下构造应力的积累变化，而忽略模型本身的重力因素.

黏弹性Maxwell体在外部载荷下的变形，不仅与边界条件及其随时间的变化有关，而且与初始条件和以前的应力演变历史有关.鉴于初始条件缺乏测量资料，只能假定为零应力、初始应变速率为0.边界条件也假定了位移速率为常量、不随时间变化.在这种条件下，开始的数百年内，下地壳、地幔等黏滞系数较小、弛豫时间较短（数十到数百年）的柔性可以占主导层位，在压缩位移下的应力增长与柔性介质内的应力松弛将达到平衡，应力维持在一个较低水平而不再增加.相反上地壳黏滞系数高（弛豫时间达数万年）的层位，在几百到两万年的期间内，弹性仍然占主导，应力随压缩几乎可以接近线性的速率增长.该结论在我们的另一研究中已经得到证明［56］.本研究中重点讨论在数百年以上的应力积累过程及其对地震活动性的影响.

我们将模型计算结果中最大剪切应力分布，在华北盆地强震震源深度处的横切面上，以及穿过若干震区的纵深剖面上显示出来.并在邢台震区进行理论钻孔，得到岩石圈各层位在不同时刻的应力增长率，以分析华北盆地岩石圈各层位对板块构造作用的时间响应.

震源定位工作表明邢台地震带（MS7.2级）、河间地震带（MS6.3级）和唐山地震带（MS7.8级）的主震震源深度在上地壳9～15km 之间［6，10，14］.因此，我们取上地壳内10km深度的横向切面上应力积累10000年后的最大剪应力分布图（见图6），以揭示应力积累与华北地震的横向空间分布的关系.

图6显示整个研究区域的应力积累水平与Moho面的起伏等值线（图2）呈现良好的对应关系.Moho面隆起处均有应力集中现象，如邢台、河间、宁河、唐山和三河地区.这些震区的最大剪应力积累水平明显地或者一定程度上高于周边地区的应力积累.太行山以西地区的最大剪应力积累水平整体低于太行山以东的华北盆地的最大剪应力积累.这种现象可以解释为：太行山以西的Moho面整体西倾（华北地壳西厚东薄）但局部较少起伏，应力积累分布均匀，而太行山以东华北盆地内地壳起伏不平，多处的Moho面隆起所造成的应力集中在整体上提高了华北盆地内部的构造应力的积累水平.比如，在邢台、河间、聊城（图6与图7）三处Moho面的隆起造成了一个三角形地带的明显的应力集中现象.边界条件作用10000年后，在整个研究区域最大剪应力积累的最大值在河间Moho面隆起处为3.2MPa，最小值在太行山以西仅为1.6MPa，两者量值近乎2倍的关系.可见Moho面上隆引起的应力积累变化作用之显著.

图6 模型10km深度处横切面上10000年后最大剪应力积累分布黑点表示华北地震主震的投影，并标出相应震级.XT：邢台地震，HJ：河间地震，NH：宁河地震，TS：唐山地震，LX：滦县地震，SH：三河地震，BH：渤海地震.Fig.6 Maximum shear stress accumulation distribution after 10000years at 10km depthBlack dots indicate projection of the hypocenters of main shocks designated with name and magnitude.

我们给出穿过邢台地震带、河间地震带与唐山地震带的纵深剖面（位置见图1中蓝实线），10000年后的最大剪应力积累水平，并将震中50km范围内地震投影到该剖面上，以分析各层应力积累水平与地震成层［6-7，10，22，57-58］分布的关系.

图7a显示了该剖面在10000年后的最大剪应力积累情况.邢台震区Moho面隆起处，宁河震区、唐山区Moho面隆起处，以及河间震区Moho面隆起处，均有明显的应力集中现象.各Moho面隆起处的上地壳底部应力集中较强，整个中地壳上层应力积累水平次之，而中地壳下层再次，下地壳没有出现应力集中，在岩石圈地幔内部的应力整体积累均匀.

图7b详细地显示了邢台震区地壳范围40km深度以上的最大剪应力积累细节.在脆性上地壳与中地壳上层内地壳分层界面的平坦处应力积累分布均匀，而在分层界面的隆起处产生应力集中现象，但是在分界面下凹处并未有应力集中出现.导致这一现象的原因是，在匀速挤压的速度边界条件下，低黏性的中地壳下层与下地壳的应力积累与松驰达到动态平衡而导致其应力不会继续增长［56］，水平载荷主要由上覆弹性层（上地壳与中地壳上层）承担，在该层较薄的Moho面上隆处会产生应力集中.10000年后的最大应力积累居于邢台震区上地壳的底部，大小为2.5MPa.中地壳上层的应力积累集中水平次之，中地壳下层再次，柔性下地壳应力集中水平最弱，其岩石黏性松驰效应使应力积累至一定水平后不再增长.

将邢台主震与余震投影到该剖面，邢台主震（MS7.2）震源深度9km处于近上地壳底部，在最大的应力增长率的范围之内；而大部分余震分布优势层在9～22km之间［6］，与应力积累量次之的脆性中地壳上层深度范围能较好地对应；少部分余震发生在较小应力积累的中地壳下层和下地壳.因此，可以推断邢台地震余震的成层分布现象与邢台地壳各层位的应力积累水平有关，而在给定的应变率的挤压边界条件下，决定该层应力积累水平的因素在于各层岩石的黏滞系数.

图7 （a）显示纵剖面（位置见图1）上10000年后的最大剪应力积累分布，（b）和（c）分别显示了（a）中所示黑色虚线方框即地壳深度范围内的最大剪应力积累分布的细节（a）中灰色的竖直实线表示在邢台震区的理论钻孔以测量整个深度范围最大剪应力随时间的变化.黑色三角表示震区地名，Xingtai：邢台，Hejian：河间，Ninghe：宁河，Tangshan：唐山，Luanxian：滦县.数字1、2、3、4、5、6分别表示岩石圈各层位，同图2.白色星星表示各次地震主震，灰色点表示余震在该剖面的投影［6，10，15-17，22］.Fig.7 （a）Maximum shear stress accumulation distribution after 10000years in the section indicated in Fig.1.（b）and（c）Stress accumulation in the area of black dash line rectangle in Fig.（a）The gray solid line in Fig.（a）indicates the theoretical drilling near the epicenter of Xingtai earthquake.Numbers mark the layers，which are as the same as those in the model of Fig.2.White stars indicate the hypocenters of the main shocks.Gray dots are the projection of hypocenters of aftershocks［6，10，15-17，22］.

唐山震区的两个Moho面隆起处唐山与宁河的上地壳与中地壳上层也有明显的应力集中现象，并且同样在各自区域的上地壳底部均有相应的最大的应力积累，应力在该区的优先积累增长，从而可能引发1978年的唐山MS7.8级主震（11km），滦县 MS7.1级余震（10km）与宁河 MS6.9级余震（17km）［6，10，15］.并且由于唐山与宁河两地上地壳底部的应力集中导致了两地之间上地壳的应力积累水平远远高于其它部分，为唐山—宁河地震带的地震发生提供应力基础.10000年后唐山地震震源处最大剪应力积累量最大值为2.9MPa.同时地震震源定位的结果［6，10，22］表明，唐山地震的余震发震优势层在10～22km，该深度在应力积累有较大水平的中地壳范围.

与此类似，河间震区Moho面隆起处，上地壳底部与整个中地壳内也有应力集中现象，最大应力积累在上地壳底部，中地壳上层应力集中程度比下层较强，而下地壳应力分布均匀，岩石圈上地幔内部的应力整体积累均匀.但是在界面下凹处并未产生应力集中现象.河间的地震记录资料较少，地震影响分布广，破坏小，国家地震数据中心将河间地震震源定为 30km（http：／／www.csi.ac.cn／manage／html／4028861611c5c2ba0111c5c558b00001／zhenli／yxy／ html／zhenli009.htm）.韦吉生等［59］通过震源定位认为2003年8月16日赤峰MS5.1级地震可能发生在25km左右的下地壳内，推论赤峰地区下地壳温度较低，可表现为脆性.如果河间地震震源深度为30km，那么，如何解释河间地震主震发生在华北盆地柔性下地壳内，将是一个比较困难的问题.

在我们建立的华北岩石圈结构模型中并未考虑渤海地区的地壳精细结构.赖晓玲等［60］结合张渤地区人工地震探测揭示的地壳结构与当地大地热流值，认为渤海地区为一上地幔隆起中心，最薄处地壳厚度仅有28km；同时该地区也是高热流值中心与地壳垂直形变中心.渤海地震位于上隆中心的边缘.根据本文以上的研究结果可以认为，渤海地震震区在长时间的应力积累过程中也应该是应力集中区.而考虑魏光兴等［25］与束沛镒等［26］将渤海地震的震源深度分别定位在30km与25km的情况，则渤海地震震源处于渤海震区下地壳范围.是否可以就此推论该地区下地壳岩石介质可能有较高的力学强度？如何解释高强度的下地壳岩石介质与高地表热流值这二者的并存关系，或需更多细致的研究工作加以探讨.将来在数值模拟渤海地区构造应力积累的研究中，亦可以考虑该地区地壳岩石流变性质与华北盆地中心区域的地壳岩石流变性质之间可能存在的差异.

为更直观地显示华北岩石圈不同深度的应力随时间增长趋势，我们以邢台震区为例，进行从地表至地下100km的理论钻孔（位置见图7），取最大剪应力在不同时刻的应力增长率，如图8所示.

从黏弹性 Maxwell体的本构关系可知，对Maxwell体介质突然施加恒定位移速率边界作用，在瞬间加载时介质表现为虎克弹性固体；在保持恒定应变速率加载作用下，介质中的应力随时间以σ0（1-exp（t／τ））形式增长，τ为黏弹性介质松驰时间，t为时间.取小时间段落（显著小于介质的松驰时间），增长速率可近似于线性；当时间显著超出介质的松驰时间范围时，介质表现基本为黏性主导，应力趋于稳定值σ0，其大小与介质黏滞系数及加载的应变速率成正比［56］，应力增长速率则趋近于0.因此邢台理论钻孔不同层位的最大剪应力的增长速率取决于各层岩石介质的松驰时间.我们模型中各层介质的松驰时间：上地壳为50000年左右，中地壳上层为8000年左右，中地壳下层为2800年，下地壳为60年，岩石圈上地幔为300年.

图8 邢台震区从地表到地下100km的理论钻孔中最大剪应力增长率随时间变化图（a—f）分别表示第1年，第500年，第2000年，第5000年，第10000年与第20000年时的应力增长率.Fig.8 Maximum shear stress rate distribution with depth in the theoretical drilling near the epicenter of Xingtai earthquake at different times

在边界条件加载的第1年（图8a），增长速率表现为弹性杨氏模量控制的线性增长，不同层位的应力增长率随相应层位的杨氏模量的增大而增加.边界条件作用500年时（图8b），加载时间超过下地壳（25～32km）的松驰时间，其应力已经趋于稳定不再增长，应力增长率为0.在应力积累5000年时（图8d），时间接近中地壳下部（19～25km）的松驰时间，该韧性层的应力松驰效应已经开始体现，应力增长率较之前有所减小；此时下地壳与岩石圈上地幔的应力增长速率已为0.当时间持续至10000年，至20000年时，中地壳上层的应力增长速率已经略有所下降，上地壳的应力能保持稳定的增长率，最大剪应力增长率峰值在上地壳底部，大小为0.25kPa／a.

5 讨 论

华北盆地构造活动复杂多样.华北盆地经历了前寒武纪克拉通化的强烈破坏和再造过程，太平洋板块向华北克拉通下的俯冲和弧后扩张作用，以及印度板块与欧亚板块碰撞的力的传递作用［10，61-63］.华北盆地板内地震可能由多种因素共同所致.现今普遍流行的观点认为华北盆地构造活动的力源来自于水平的板块构造作用或者垂直向的地幔-岩石圈的相互作用，或者两者的共同作用.华北盆地的地壳结构在地震孕育中提供了基础条件，比如：Moho面的隆起和震源下方的低速包裹体的存在.

我们的计算主要考虑了在板块构造挤压的边界条件下，华北盆地的Moho面的隆起对应力积累的影响.尹京苑和梅世蓉［33］曾就此问题进行了讨论.其结果认为Moho面局部隆起对壳内应力场产生了非常大的影响（超过深、浅断层对应力的影响规模）.在这一点上，我们的结论与尹京苑和梅世蓉的观点一致.但是，他们的结论中进一步强调Moho面隆起造成的应力增强的密集带不在隆起区的正上方，而是在其斜坡上.并且，孙武城等［34］设计的光弹模拟实验测定在Moho面呈波浪形起伏的情况下构造应力场的分布状态，其结果也认为地壳中最大剪应力既不在波峰处，也不在波谷处，而是位于Moho面起伏的“拐点”上方，支持尹京苑和梅世蓉的结论.不可否认基于弹性理论的理论计算与光弹实验结果的正确性，但是依据现今不论由野外观测还是室内岩石实验所揭示的华北盆地岩石圈流变分层结构的认识而言，考虑长时期的应力积累过程的时候，对有不同流变结构的地壳统一使用弹性性质建立模型似乎并不妥当.地震的孕育和发生不仅仅与地壳非均匀性引起的弹性应力场的扰动值有关，更重要的是与该扰动值随时间的变化有关.在实际的长期地质构造作用的地壳应力积累过程中，上、中、下地壳将由各自的岩石流变性质决定其不同的力学表现，或脆性、或韧性、或者脆韧性转换.一般来说，在较短的时间尺度内考虑地壳的变形，可以用弹性变化来描述，如瞬间的地震应力触发.但是长期的构造应力作用下强震的孕育过程中，就必须考虑岩石的蠕变性质.

我们的计算结果显示应力积累集中在Moho面起伏的上隆处，而不是在拐点与波谷处.本文的模型中包含有脆性的上地壳与中地壳上层、以及韧性的中地壳下层、下地壳与岩石圈上地幔.在长时间的挤压的边界条件下，介质分层界面的上隆处会因挤压而产生压应变集中，同时由于黏性中地壳下层与下地壳的应力松驰效应，两者的共同作用使上伏弹性层（上地壳与中地壳上层）在界面上隆处应力得到不断的积累产生较大的应力集中现象，快速的应力增长有利于在该处应力首先达到岩石的破裂强度而发生优势破裂，引发主震.

在中地壳应力的增长速率次之，一般认为深部岩石的破裂强度比浅部更大，构造应力的积累尚未达到该层岩石的破裂强度.但是，在不论是邢台、河间还是唐山的地震带中，大部分的余震都集中在该深度范围，我们认为其中的可能性为主震的应力触发作用［64-65］.

下地壳应力由于黏性松驰效应导致其增长速率随时间不断减小，且下地壳岩石的强度与浅部岩石相比更大，构造应力的积累与主震的应力触发的共同作用也很难使其产生破裂.因此该层有较少的地震记录出现.

对于震源下方速度异常体对应力积累的控制作用，朱守彪等［1］在最新的研究中认为：无论异常区是低速异常还是高速异常，在水平构造力的作用下（挤压或是引张力作用），都能造成低速体上方的应力集中现象，并且挤压的边界条件下，弹性低速体上方最大剪应力的增长率为0.41kPa／a.在我们的计算结果中，Moho上隆的邢台震源处的最大剪应力增长速率为0.25kPa／a.震源下方速度异常体与Moho面隆起这两种因素导致的应力积累量值相当.但是需要指出的是，在朱守彪计算所使用的介质中低速度体包裹的形状为横向长轴为18km，纵向短轴为10km的椭圆.该椭圆的曲率与华北盆地部分地区Moho面的起伏曲率相当.因此该计算结果中的最大剪应力积累计算结果0.41kPa／a同时包含了界面起伏与低速包裹体这两者共同作用造成的应力集中效应.如果简单将本文得到的由Moho面起伏所造成的应力集中速率0.25kP／a从朱守彪的研究结果0.41kPa／a中减去，那么单纯低速体所造成的应力集中速率为0.16kPa／a.

袁金荣等［62］认为华北盆地地震可能由地幔对流作用在Moho面隆起的下方造成的上拱力所致.曾融生等［66］提出华北盆地强震以及凹陷形成的模型，认为在水平板块构造应力场的背景中，上地幔热物质向地壳下部入侵，它所产生的扰动应力场不仅在横向是不均匀的，而且在垂向也不均匀.地幔的入侵能够在地壳上部产生足够大的伸张应力场，同时在地壳中部或下部产生水平切应力场.冯德益等［32］通过数值模拟认为地幔对流的上拱力引起了地壳的应力场的变化.因其计算结果的应力与本文不同，因此不方便统一比较.

总之，Moho面的隆起（地壳薄弱带），如震源下方低速度体的存在与壳幔相互作用一样，都在一定程度上对华北盆地的地震孕育起着重要的作用.而Moho面的起伏在地壳应力集中的控制作用的现象并不为华北盆地地震现象中所独有.在2008年汶川MS8.0级大地震的龙门山断裂带处，地壳厚度在龙门山断裂带从青藏高原向四川盆地的急剧减小近20km的因素，为龙门山断裂带在地壳缩短的构造挤压情况下的应力集中积累提供了重要控制作用［56］.

曾融生等［66］认为华北盆地地壳的力学性质随深度而改变，强震可能是由中部地壳的塑性形变以及上部地壳的脆性断层所组成的，即所谓两层破裂的震源模型.臧绍先等［47］计算了华北岩石圈流变强度及黏度的三维空间分布.其结果认为上地壳上部为脆性区，中地壳可以大部分是以脆性破裂为主的脆性区，部分是以蠕变为主的延性区，而下地壳几乎均是以蠕变为主的延性区；壳下岩石圈上部是以脆性破裂为主或以蠕变为主的高强度区.通过本文模型分层应力积累结果表明，最大剪应力增长率最大值在相应震区的脆性上地壳底部，在脆性中地壳上部应力能保持持久增长，而在韧性的中地壳下部应力增长率次之，在更加韧性的下地壳应力增长速率为0.这种应力积累特征与华北盆地的主震深度、余震成层分布结果相对应.我们的计算结果支持华北岩石圈地壳流变层状结构，即脆性（上地壳）－脆性较弱（中地壳上层）－较弱韧性（中地壳下层）－较强韧性（下地壳）－韧性（岩石圈上地幔）的分层结构.因此华北盆地软弱下地壳的存在，使得华北盆地岩石圈的流变结构很难用统一的强地壳-弱地幔的形式来概括描述.

与此同时，华北盆地软弱的下地壳的存在也正好解释了，从Airy重力均衡的角度而言的华北盆地平坦地表与Moho面起伏之间的矛盾［67］，而且这一现象在世界各地均有体现.石耀霖等［68］曾论述一些山没有山根［69］（阿巴拉契亚山脉的Scottish部分），美国盆岭区引张系数达1.4～4.0，但 Moho面平坦而没有上涌［70］，类似现象在英国大陆架深反射探测中也被观测到［71］，这些都可以用下地壳岩石弱于Moho面下的地幔的模型来解释；较强的Moho面没有沉陷或上涌，而是软弱的下地壳岩石横向流动补偿了地表载荷.

目前，对华北地区构造应力场现有研究表明看法尚不一致［62］.一是依据地质构造体系的发育，认为该区为NWW—SEE向的拉张应力场；二是依据地震的震源机制解、地应力解除测量与GPS观测资料，认为该区为NEE—SWW的水平挤压应力场.

但是绝大多数研究者均支持华北地区有稳定的NEE—SWW 方面的最大主压应力，比如：李钦祖［53］、汪素云［72］、刁桂苓［7］、白武明［55］、张宏志［22］和陈连旺［63］由华北地震余震震源机制与华北构造应力场数值模拟给出的结果，认为华北构造应力场存在较好的一致性，其最大压应力轴为NEE，最小压应力轴为NNW向，并且都近于水平向.李延兴［54］根据近十多年来华北地区GPS网的观测资料，分别计算了邢台、渤海、海城和唐山4次MS7.0级以上大地震震中区现今应变场的主应变参数.结果表明，由GPS观测得到的邢台、渤海、唐山地震震中区的主压应变轴方向与震前震中及周围地区的主压应力轴方向是一致的或基本上是一致在NEE方向.因此本文所采用在NEE75°方向的挤压边界条件是合理的.

华北盆地断层发育丰富，但是在本文的模型中，没有考虑断层.在应力数值模拟中，一般将断层考虑为低速体，低速体的存在势必会使断层周围产生应力的集中现象［1］，因此在本研究中忽略断层，以单独考虑Moho面隆起这一因素对构造应力积累集中的贡献.

6 结 论

本文结合地震探测所揭示的华北盆地岩石圈结构，使用通过基于大地热流测量的温度计算与实验室高温高压岩石流变实验研究结果，建立了华北盆地岩石圈三维黏弹性有限元模型，从零应力下开始，在NEE方向上施加按匀速挤压的位移边界条件，数值模拟华北岩石圈各层位在长期板块构造挤压作用下的应力积累过程.分析了华北地震空间分布与构造应力积累的关系，探讨了地壳结构与地壳流变分层性质对应力积累的影响.得到以下结论：

（1）Moho面隆起与地壳岩石各层位的黏滞系数是华北盆地地震孕育的重要的控制因素.Moho面隆起区域往往是地幔上涌、温度较高.高黏度脆性上地壳较薄的部位，在长期均恒挤压变形速率下，出现应力集中现象.华北盆地若干震区如邢台、河间、唐山、渤海均处于在Moho面隆起上方.

（2）华北震区的应力在长时期的积累过程中，在脆性的上地壳与中地壳上层，应力表现近于线性增长趋势，上地壳底部较其它深度有最大值的应力增长率（如邢台震源深度处的最大剪应力增长率为0.25kPa／a），主震可以在应力积累至岩石破裂强度时被促发；在脆、韧性转换的中地壳，应力增长速率次之，但是由于其弹性（弹性模量）较高，因此在载荷速率突然变化（边界位移或区域内部发生大地震而造成应力瞬间变化时），其应力变化率高，华北地震的大部分余震可能在该层位为主震所触发；而在柔性的下地壳应力增长近于指数形式，稳定状态之后其应力增长速率为零.

（3）模型中应力在不同层位的积累速率的差异与华北盆地地震的成层分布现象支持华北岩石圈流变结构的分层特征：脆性（上地壳）－较弱脆性（中地壳上层）－较弱韧性（中地壳下层）－较强韧性（下地壳）－韧性（岩石圈上地幔）的分层结构.

（References）

［1］朱守彪，张培震，石耀霖.华北盆地强震孕育的动力学机制研究.地球物理学报，2010，53（6）：1409－1417.Zhu S B，Zhang P Z，Shi Y L.A study on the mechanisms of strong earthquake occurrence in the North China basin.Chinese J.Geophys.（in Chinese），2010，53（6）：1409－1417.

［2］Gao S，Roberta R L，Yuan H L，et al.Recycling lower continental crust in the North China craton.Nature，2004，432（7019）：892－897.

［3］邓晋福，莫宣学，赵海玲等.中国东部岩石圈根／去根作用与大陆“活化”──东亚型大陆动力学模式研究计划.现代地质，1994，8（3）：349－356.Deng J F，Mo X X，Zhao H L，et al.Lithosphere root／de－rooting and activation of the east China continent.Geoscience（in Chinese），1994，8（3）：349－356.

［4］Menzies A M，Fan W，Zhang M.Palaeozoic and Cenozoic lithoprobes and the loss of＞120km of Archaean lithosphere，Sino－Korean craton，China.The Geological Society，1993，76（1）：71－81.

［5］Zhu R X，Zheng T Y.Destruction geodynamics of the North China Craton and its Paleoproterozoic plate tectonics.Chinese Science Bulletin，2009，54（19）：3354－3366.

［6］毛桐恩，陈锦标，姚家榴.邢台、唐山震区震源层及其顶、底部介质条件的综合分析.地震学报，1996，18（1）：89－96.Mao T E，Chen J B，Yao J L.Analysis of the upper and bottom boundary of Xingtai and Tangshan earthquake zone.Acta Seismologica Sinica （in Chinese），1996，18（1）：89－96.

［7］刁桂苓，于利民，李钦祖等.唐山和澜沧地震序列震源区应力场的对比分析.地震学报，1995，17（3）：305－311.Diao G L，Yu L M，Li Q Z，et al.Comparison of the stress field between Tangshan and Lancang earthquake series.Acta Seismologica Sinica （in Chinese），1995，17（3）：305－311.

［8］梅世蓉.邢台地震孕育发生模型及其前兆机理探讨.地震，1999，19（1）：1－10.Mei S R.Model of preparation and generation of xingtai earthquake and mechanism of its precursors.Earthquake（in Chinese），1999，19（1）：1－10.

［9］邵学钟，张家茹，章思亚等.邢台地震区深部构造背景的地震转换波探测和研究.地球物理学报，1993，36（5）：610－620.Shao X Z，Zhang J R，Zhang S Y，et al.Study of deep structures in xingtai earthquake area by method of converted waves of earthquakes.Chinese J.Geophys.（in Chinese），1993，36（5）：610－620.

［10］曾融生，张少泉，周海南等.唐山地震区的地壳结构及大陆地震成因的探讨.地震学报，1985，7（2）：126－140.Zeng R S，Zhang S Q，Zhou H N，et al.Crustal structure of Tangshan epicentral region and its relation to the seismogenic process of a continental earthquake.Acta Seismologica Sinica （in Chinese），1985，7（2）：126－140.

［11］邵学钟，张家茹，章思亚等.唐山大震区深部构造的初步研究.地球物理学报，1986，29（1）：29－40.Shao X Z，Zhang J R，Zhang S Y，et al.Investigation of deep structures in Tangshan earthquake area.Chinese J.Geophys.（in Chinese），1986，29（1）：29－40.

［12］刘启元，王峻，陈九辉等.1976年唐山大地震的孕震环境：密集地震台阵观测得到的结果.地学前缘，2007，14（6）：206－212.Liu Q Y，Wang J，Chen J H，et al.Seismogenic tectonic environment of 1976great Tangshan earthquake：results given by dense seismic array observations.Earth Science Frontiers（in Chinese），2007，14（60）：205－213.

［13］王峻，刘启元，陈九辉等.首都圈地区的地壳厚度及泊松比.地球物理学报，2008，52（1）：58－65.Wang J，Liu Q Y，Chen J H，et al.The crustal thickness and poisson′s ratio beneath the capital circle region.Chinese J.Geophys.（in Chinese），2008，52（1）：58－65.

［14］Shedlock K M，Baranowsk J，Xiao W W，et al.The Tangshan aftershock sequence.J.Geophys.Res.，1987，92（B3）：2792－2803.

［15］藏绍先，杨军亮.我国华北等地区板内地震的深度分布及其物理背景.地震地质，1984，6（5）：67－75.Zang S X，Yang J L.Distribution of focal depths for intraplate earthquakes in North China and its physical explanation.Seismology and Geology （in Chinese），1984，6（5）：67－75.

［16］Nábělek J，Chen W P，Ye H.The Tangshan earthquake sequence and its implications for the evolution of the North China basin.J.Geophys.Res.，1987，92：12615－12628.

［17］叶洪，张文郁.初论华北板内地震断层特征及其与地壳、上地幔动力学过程的关系.地震地质，1980，2（1）：23－37.Ye H，Zhang W Y.The characteristics of intraplate earthquake faults in North China and their relationship to the dynamical processes in earth′s crust and uppermost mantle.Seismology and Geology （in Chinese），1980，2（1）：23－37.

［18］陈运泰，林邦慧，林中洋等.根据地面形变的观测研究1966年邢台地震的震源过程.地球物理学报，1975，18（3）：164－182.Chen Y T，Lin B H，Lin Z Y，et al.The focal mechanism of the 1966Hsingta earthquake as inferred from the ground deformation observations.Chinese J.Geophys.（in Chinese），1975，18（3）：164－182.

［19］徐锡伟，于贵华，王峰等.1966年邢台地震群的发震构造模型——新生断层形成？先存活断层摩擦粘滑？中国地震，2000，16（4）：364－378.Xu X W，Yu G H，Wang F，et al.Seismogenic model for the 1966Xingtai Earthquakes－nucleation of new－born fault or strick－slip of pre－existing fault？Earthquake Research in China （in Chinese），2000，16（4）：364－378.

［20］陈祖安，白武明，林邦慧等.1966年以来华北地区一系列七级大震破裂过程的数值模拟.地球物理学报，2003，46（3）：373－381.Chen Z A，Bai W M，Lin B H，et al.Numerical simulation for rupture processes of a series of strong earthquakes（MS＞7）in North China since 1966.Chinese J.Geophys.（in Chinese），2003，46（3）：373－381.

［21］张之立，王成宝，方兴等.唐山地震破裂过程的雁行断裂模式及理论和试验的模拟.地震学报，1989，（3）：291－301.Zhang Z L，Wang C B，Fang X，et al.An echelon fracture model for the Tangshan earthquake sequence as well as simulation in theory and experiment.Acta Seismologica Sinica （in Chinese），1989，（3）：291－301.

［22］张宏志，刁桂苓，赵英萍等.邢台地区近年的震源机制.大地测量与地球动力学，2007，27（6）：91－95.Zhang H Z，Diao G L，Zhao Y P，et al.Recent focal mechanism in Xingtai region.Journal of Geodesy and Geodynamics（in Chinee），2007，27（6）：91－95.

［23］潘波，许建东，刘启方.1679年三河-平谷8级地震近断层强地震动的有限元模拟.地震地质，2009，31（1）：69－83.Pan B，Xu J D，Liu Q F.Simulations of the near－fault strong ground motion of the 1679Sanhe－Pinggu M8.0earthquake.Seismology and Geology （in Chinese），2009，31（1）：69－83.

［24］冉勇康，邓起东，杨晓平等.1679年三河—平谷8级地震发震断层的古地震及其重复间隔.地震地质，1997，19（3）：194－201.Ruan Y K，Deng Q D，Yang X P，et al.Paleo－earthquakes and recurrence interval on the seismogenic fault of 1697 Sanhe－Pinggu M8.0earthquake，Hebei and Beijing.Seismology and Geology （in Chinese），1997，19（3）：194－201.

［25］魏光兴，季同仁，李秉峰.渤海地震序列及其特征.地震地质，1984，6（1）：22－29.Wei G X，Ji T R，Li B F.Characteristics of magnitude 7.4 earthquake equence in the Bohai area.Seismology and Geology （in Chinese），1984，6（1）：22－29.

［26］束沛镒，李幼铭，铁安等.利用远震P波波形反演渤海地震的震源参数.地球物理学报，1983，26（1）：31－38.Shu P Y，Li Y M，Tie A，et al.An inversion of the focus parameters of Bohai earthquake from the teleseismic P wave shape.Chinese J.Geophys.（in Chinese），1983，26（1）：31－38.

［27］宋惠珍，高维安，孙君秀等.唐山地震震源应力场的数值模拟研究——三维有限单元法在计算震源应力场中的应用.西北地震学报，1982，（3）：32－38.Song H Z，Gao W A，Sun J X，et al.Numerical modelling the stress field in Tangshan earthquake—application of 3D modelling stress field by FEM.Northwest China Seismological Journal（in Chinese），1982，（3）：32－38.

［28］郑治真，刘元壮，胡祚春.唐山地震孕育过程的数学模拟——（二）地震孕育过程的数学模拟.地震研究，1984，（3）：51－58.Zhen Z Z，Liu Y Z，Hu Z C.The numerical simulation of preparatory process of Tangshan earthquake.（2）A numerical simulation of preparatory process of earthquake.Journal of Seismological Research （in Chinese），1984，（3）：51－58.

［29］张东宁，曾融生.冀中坳陷滑脱构造动力的数值模拟.地震学报，1995，17（4）：414－421.Zhang D N，Zeng R S.Numerical modelling the detachment of Jizhong depression area.Acta Seismologica Sinica （in Chinese），1995，17（4）：414－421.

［30］王继存，续春荣.唐山地震断层破坏及其力学过程的数值模拟.地震地质，1989，11（4）：71－76.Wang J C，Xu C R.Fault failure by Tangshan earthquake and numerical analogue of its mechanical process.Seismology and Geology （in Chinese），1989，11（4）：71－76.

［31］梅世蓉，梁北援.唐山地震孕震过程的数值模拟.中国地震，1989，（3）：42－47.Mei S R，Liang B Y.A mathematical simulation for seismogenetic process of the Tangshan earthquake.Earthquake Research in China （in Chinese），1989，（3）：42－47.

［32］冯德益，刘喜兰，舒立德等.唐山7.8级强震孕育条件及前兆场时空演化的三维数值模拟计算与分析.华北地震科学，1996，（4）：43－52.Feng D Y，Liu X L，Shu L D，et al.Three－dimensional numerical simulation and analysis of seismogenic condition and temporo－spatial evolution of precursor field of Tangshan M＝7.8great earthquake.North China Earthquake Sciences（in Chinese），1996，（4）：43－52.

［33］尹京苑，梅世蓉.邢台强震区的深部构造对强震孕育影响的三维数值模拟.地球物理学报，1999，42（增刊）：131－140.Yin J Y，Mei S Y.3－D simulation for the effect of deep structures to the preparation of strong earthquake in Xingtai area.Chinese J.Geophys.（in Chinese），1999，42（Suppl.）：131－140.

［34］孙武城，李松林，张晓普.华北大陆地壳-上地幔波浪断裂型构造的光弹模拟实验.华北地震科学，1983，1（1）：21－27.Sun W C，Li S L，Zhang X P.Photoelastic experiment on the failure structure of the fluctuating Moho surface in North China continent. North China Earthquake Sciences （in Chinese），1983，1（1）：21－27.

［35］肖兰喜，朱元清，张少泉等.邢台震区深部构造与强震孕育关系的探讨.地震学报，1999，21（6）：597－607.Xiao L X，Zhu Y Q，Zhang S Q，et al.The relationship between the deep－level structure in crust and brewing of strong earthquakes in Xingtai area.Acta Seismologica Sinica（in Chinese），1999，21（6）：597－607.

［36］宋治平，尹祥础，梅世蓉.地震孕育体源流变模型（二）——应变场及其应用.地震学报，2004，26（2）：121－131.Song Z P，Yin X C，Mei S R.3－D rheologic model of earthquake preparation（2）.Strain field and its applications.Acta Seismologica Sinica （in Chinese），2004，26（2）：121－131.

［37］金安蜀，刘福田，孙永智.北京地区地壳和上地幔的三维P波速度结构.地球物理学报，1980，23（2）：172－182.Jin A S，Liu F T，Sun Y Z.Three－dimensional P velocity structure of the crust and upper mantle under Beijing region.Chinese J.Geophys.（in Chinese），1980，23（2）：172－182.

［38］刘福田，曲克信，吴华等.华北地区的地震层面成象.地球物理学报，1986，29（5）：442－449.Liu F T，Qu K X，Wu H，et al.Seismic tomography of North China region.Chinese J.Geophys.（in Chinese），1986，29（5）：442－449.

［39］林真明，邵学钟，陈学波.1966年邢台地震区地震测深资料再解释.华北地震科学，1990，（3）：34－37.Lin Z M，Shao X Z，Chen X B.Reanalyzing for the seismic sounding data in Xingtai seismic area in 1966.North China Earthquake Sciences（in Chinese），1990，（3）：34－37.

［40］孙若昧，刘福田.京津唐地区地壳结构与强震的发生——1.P波速度结构.地球物理学报，1995，38（5）：599－607.Sun R M，Liu F T.Crust structure and strong earthquake in Beijing，Tianjing，Tangshan area：I.P wave velocity structure.Chinese J.Geophys.（in Chinese），1995，38（5）：599－607.

［41］于湘伟，陈运泰，王培德.京津唐地区中上地壳三维P波速度结构.地震学报，2003，25（1）：1－14.Yu X W，Chen Y T，Wang P D.Three－dimensional P wave velocity structure in Beijing－Tianjing－Tangshan area.Acta Seismologica Sinica （in Chinese），2003，25（1）：1－14.

［42］梅世蓉，梁北援.唐山地震孕震过程的数值模拟.中国地震，1989，（3）：42－47.Mei S R，Liang B Y.A mathematical simulation for seismogenetic process of the Tangshan earthquake.Earthquake Research in China （in Chinese），1989，（3）：42－47.

［43］张家茹，邵学钟，章思亚等.邢台强震的深部三维构造背景和孕震条件.地震地质，1994，（4）：21－27.Zhang J R，Shao X Z，Zhang S Y，et al.Background of 3－D deep structure and dynamic conditions for occurrence of strong Xingtai earthquake.Seismology and Geology （in Chinese），1994，（4）：21－27.

［44］赖晓玲，张先康，郑需要.利用三维弯曲界面的反演方法重建唐山震区莫霍界面三维构造形态.地震学报，1997，19（5）：506－516.Lai X L，Zhang X K，Zheng X Y.Rebuilding the Moho surface of Tangshan earthquake zone using bending surface inversion.Acta Seismologica Sinica （in Chinese），1997，19（5）：506－516.

［45］An M J，Feng M，Zhao Y.Destruction of lithosphere within the north China craton inferred from surface wave tomography.Geochemistry Geophysics Geosystems，2009，10：Q08016.

［46］Huang J L，Zhao D P，Zheng S H.Lithospheric structure and its relationship to seismic and volcanic activity in southwest China.J.Geophys.Res.，2002，111：2255.

［47］臧绍先，李昶，宁杰远等.华北岩石圈三维流变结构的一种初步模型.中国科学（D），2002，32（7）：588－597.Zang S X，Li C，Ning J Y，et al.3Drheology model of the lithosphere of North China basin.Science China（D）（in Chinese），2002，32（7）：588－597.

［48］安美建，石耀霖.中国大陆地壳和上地幔三维温度场.中国科学（D），2007，37（6）：736－745.An M J，Shi Y L.3Dtemperature of crust and lithosphere mantle in Chinese continent.Science China（D）（in Chinese），2007，37（6）：736－745.

［49］石耀霖，曹建玲.中国大陆岩石圈等效粘滞系数的计算和讨论.地学前缘，2008，15（3）：83－95.Shi Y L，Cao J L.Effective viscosity of China continental lithosphere.Earth Sciences Frontier（in Chinese），2008，15（3）：83－95.

［50］Li Z H，Gerya T V.Polyphase formation and exhumation of high－to ultrahigh－pressure rocks in continental subduction zone：Numerical modeling and application to the Sulu ultrahigh－pressure terrane in eastern China.J.Geophys.Res.，2009，114：B09406.

［51］Li Z H，Gerya T V，Burg J P.Influence of tectonic overpressure on P－T paths of HP－UHP rocks in continental collision zones：Thermomechanical modeling.Journal of Metamorphic Geology，2010，28（3）：227－247.

［52］Li Z H，Xu Z Q，Gerya T V.Flat versus steep subduction：contrasting modes for the formation and exhumation of high－to ultrahigh－pressure rocks in continental collision zones.Earth and Planetary Science Letters，2011，301（1－2）：65－77.

［53］李钦祖.华北地壳应力场的基本特征.地球物理学报，1980，23（4）：376－388.Li Q Z.General features of the stress field in the crust of North China.Chinese J.Geophys.（in Chinese），1980，23（4）：376－388.

［54］李延兴，徐杰，陈聚忠等.邢台、渤海、海城和唐山大地震震中区现今应变场的基本特征.华北地震科学，2006，24（2）：36－39.Li Y X，Xu J，Chen J Z，et al.Basic characteristics of current strain field in the epicentral area of Xingtai，Bohai，Haicheng and Tangshan strong earthquake.North China Earthquake Sciences（in Chinese），2006，24（2）：36－39.

［55］白武明，林邦慧，陈祖安.1976年唐山大震发生对华北地区各地块运动与变形影响的数值模拟研究.中国科学，2003，33（增刊）：99－107.Bai W M，Lin B H，Chen Z A.Numerical modelling the block motion and deformation induced by Tangshan earthquake.Science in China （D）（in Chinese），2003，33（Suppl.）：99－107.

［56］柳畅，朱伯靖，石耀霖.粘弹性数值模拟龙门山断裂带应力积累及大震复发周期.地质学报，2012，86（1）：157－169.Liu C，Zhu B J，Shi Y L.Stress accumulation of the Longmenshan fault and recurrence interval of Wenchuan earthquake based on viscoelasticity simulation.Acta Geologica Sinica （in Chinese），2012，86（1）：157－169.

［57］李轶群，王健.唐山余震区中小地震震源机制解分区特征的初步研究.中国地震，2008，24（2）：150－157.Li Y Q，Wang J.The preliminary study on focal mechanism of small and moderate earthquakes in divided zones around Tangshan.Earthquake Research in China （in Chinese），2008，24（2）：150－157.

［58］兰从欣，邢成起，苗春兰等.近年首都圈地区中小地震震源机制解及其特征分析.华北地震科学，2005，24（3）：22－25.Lan C X，Xing C Q，Miao C L，et al.The focal mechanism solutions and their characteristics for moderate and small earthquakes in Capital Area in recent years.North China Earthquake Sciences（in Chinese），2005，24（3）：22－25.

［59］韦生吉，倪四道，崇加军等.2003年8月16日赤峰地震：一个可能发生在下地壳的地震？地球物理学报，2009，52（1）：111－119.Wei S J，Ni S D，Chong J J，et al.The 16August 2003 Chifeng earthquake：Is it a lower crust earthquake？Chinese J.Geophys.（in Chinese），2009，52（1）：111－119.

［60］赖晓玲，李松林，孙译.渤海及邻区3次7级以上地震的深部构造背景.大地测量与地球动力学，2007，27（1）：31－33，54.Lai X L，Li S L，Sun Y.Deep tectonic background of three MS＞7.0strong earthquakes in Bohai and its adjacent region.Journal of Geodesy and Geodynamics （in Chinese），2007，27（1）：31－33，54.

［61］徐菊生，袁金荣，高士钧.利用GPS观测结果研究华北地区现今构造应力场.地壳形变与地震，1999，19（2）：81－89.Xu J S，Yuan J R，Gao S J.Research on the precent tectonic stress field in North China with GPS observations.Crustal Deformation and Earthquake （in Chinese），1999，19（2）：81－89.

［62］袁金荣，徐菊生，高士钧.利用GPS观测资料反演华北地区现今构造应力场.地球学报，1999，20（3）：232－238.Yuan J R，Xu J S，Gao S J.Back analysis of present tectonic stress field in North China region using GPS data.Acta Geoscientia Sinica （in Chinese），1999，20（3）：232－238.

［63］陈连旺，陆远忠，张杰等.华北地区三维构造应力场.地震学报，1999，21（2）：140－149.Chen L W，Lu Y Z，Zhang J，et al.There dimensional tectonic stress field in North China.Acta Seismologica Sinica（in Chinese），1999，21（2）：140－149.

［64］Harris R A.Introduction to special section：Stress triggers，stress shadows，and implications for seismic hazard.J.Geophys.Res.，1998，103（B10）：24347－24358.

［65］石耀霖.关于应力触发和应力影概念在地震预报中应用的一些思考.地震，2001，21（3）：1－7.Shi Y L.Stress triggers and stress shadows：How to apply these concepts to earthquake prediction.Earthquake，2001，21（3）：1－7.

［66］曾融生，朱露培，何正勤等.华北盆地强震的震源模型兼论强震和盆地的成因.地球物理学报，1991，34（3）：288－301.Zeng R S，Zhu L P，He Z Q，et al.Model of the earthquake source of North China basin strong earthquake and discussion the mechanism of extensional basin.Chinese J.Geophys.（Acta Geophysica Sinica）（in Chinese），1991，34（3）：288－301.

［67］熊熊，许厚泽.地幔对流诱发的莫霍面起伏与下地壳塑性流动.测绘学报，2000，29（增刊）：21－25.Xiong X，Xu H Z.Undulation of Moho and ductile flow of lower crust induced by mantle flow.Acta Geodaetica et Cartographica Sinica （in Chinese），2000，29（Suppl.）：21－25.

［68］石耀霖，朱守彪.中国大陆震源机制深度变化反映的地壳、地幔流变特征.地球物理学报，2003，46（3）：359－365.Shi Y L，Zhu S B.Contrast of rheology in the crust and mantle near Moho revealed by depth variation of earthquake mechanism in continental China.Chinese J.Geophys.（in Chinese），2003，46（3）：359－365.

［69］McKenzie D P，Nimmo F，Jackson A J，et al.Characteristics and consequences of flow in the lower crust.J.Geophys.Res.，2000，105（B5）：11029－11046.

［70］Gans P B.An open－system two layer crustal stretching model for the eastern Great basin.Tectonics，1987，6（1）：1－12.

［71］Kusznir H N，Matehews H D.Deep seismic reflections and the deformational mechanics of the continental lithosphere.J.Petrol.，1988，Special Lithosphere Issue：63－87.

［72］汪素云，陈培善.中国及邻区现代构造应力场的数值模拟.地球物理学报，1980，23（1）：35－45.Wang S Y，Chen P S.A numerical simulation of the present tectonic stress field of china and its vicinity.Chinese J.Geophys.（Acta Geophysica Sinica）（in Chinese），1980，23（1）：35－45.