GPS约束下唐山地区断裂带现今运动速度特征非连续接触模拟*

孟庆筱

1)中国地震局第一监测中心，天津 300180

2)中国地震局地震研究所(地震大地测量重点实验室)，武汉 430071

3)中国地震局地壳应力研究所武汉科技创新基地，武汉 430071

GPS约束下唐山地区断裂带现今运动速度特征非连续接触模拟*

孟庆筱1，2，3)

1)中国地震局第一监测中心，天津 300180

2)中国地震局地震研究所(地震大地测量重点实验室)，武汉 430071

3)中国地震局地壳应力研究所武汉科技创新基地，武汉 430071

以1999～2007年GPS观测结果为基础，使用二维接触有限元方法计算唐山地区主要断裂带的运动速度状态。模拟结果表明，现今地壳形变运动状态下唐山地区主要断裂带运动状态具有差异性。由于山海关块隆对唐山块陷的推挤作用，唐山断裂带以右旋拉张运动为主，同时在蓟运河断裂带位置处保持着快速的挤压应力积累，具有较高的地震危险性。

唐山;GPS;有限元;接触单元;断裂带

研究认为，1976年唐山7.8级强烈地震与深大断裂没有表现出直接联系，而是发生在唐山块陷内部北东向的唐山断裂带上。边界断裂的存在使得被围陷的唐山块体得到暂时的平衡，形成相对锁闭的区域。区域构造应力场发生急剧变化时，锁闭平衡状态被打破，从而发生唐山大地震［1］。2012-05-20唐山地区又发生M4.7地震，说明唐山断裂构造的活动性明显增强。研究唐山地区断裂带的现今形变特征，对于唐山地区地震危险性评价具有重要意义。本文通过构建唐山地区二维有限元接触模型，利用中国大陆网络工程复测资料1999～2007年GPS结果进行边界条件约束，使用经过信度检验的有限元模拟结果分析唐山断裂、丰台-野鸡坨断裂、蓟运河断裂、宁河-昌黎断裂、大李庄-古马断裂等唐山地区主要断裂带的现今形变特征，为唐山地区地震危险性分析提供依据。

1 唐山地区构造概况

唐山块陷位于阴山-燕山断褶带和冀鲁断块坳陷之间，四周被深大断裂所包围，块体中央有唐山断裂带，唐山地震就发生在这个块体内部的断裂带上。唐山块陷在形态上为一个走向北东东的菱形块体，四周被具有新活动的4条大断裂和深断裂所包围，断块的南界为宁河-昌黎深断裂，北界为丰台-野鸡坨大断裂，二者走向均为北东东;东界是滦县-乐亭大断裂，西界为蓟运河深断裂，二者走向均成北西向。其活动特征分述如下［1］:

1)宁河-昌黎深断裂:为唐山块陷和乐亭块陷的分界线，向北东延入渤海，向西与沧东断裂相接。断裂产生于吕梁运动，新活动强烈，断裂两盘第四纪沉积物落差500 m以上。在断裂带的西北盘上，有一平行于该断裂的北东东向河流纵比降转折带，与断裂的北西盘上升、南东盘下降有关。

2)丰台-野鸡坨大断裂:为蓟县块隆和唐山块陷的交界线，控制古生代地层的发育。断裂以北为巨厚的震旦系地层，以南主要为古生代地层，西段断面倾向西北，东段倾向南东。在榛子镇附近为断面转折的枢纽部位，新构造运动明显。东段是山区和平原的分界线，西段断裂两侧第四系厚度之差可达300～400 m，丰润西南沿断裂出现第四系沉积物厚度等值线的密集带。

3)滦县-乐亭大断裂:断裂控制山海关块隆和唐山块陷的活动，长期前者上升、后者下沉。断面倾向北东，为高角度的逆冲断层。在马城-野鸡坨之间有4条平行断层，形成于前古生代，新构造运动强烈。

4)蓟运河深断裂:为唐山块陷和沧县块隆的分界线，向西北延伸，与宝坻-香河东西向断裂相连。香河宝坻段为重力梯度带，东南段梯度缓断裂控制古生代地层的发育和中生代盆地的形成。断裂形成于前古生代，新构造运动强烈，不但控制蓟运河的流向，而且以断裂为界，两侧水系流向截然不同。

5)唐山断裂带:该断裂由唐山-古冶断裂、陡河断裂和唐山-巍山-长山南坡断裂3条平行断裂组成，其中唐山-古冶断裂南西段走向N30°E，东北段走向N50°E;陡河断裂东北段为一向北西倾的正断层，西南段由平行的4条小断层组成，断面倾向北西。最西一条为正断层，东边3条为逆断层，全长约50 km;唐山-巍山-长山南坡断裂由一些断断续续的北东向断层组成，多为向北西倾的逆断层。为便于建模，将以上3条断裂合并为唐山断裂带。

根据以上论述，多期地质构造运动导致唐山地区呈现为被深大断裂围陷、分割的复杂块体结构。由于地壳深部各圈层之间存在的相互作用在构造线位置最为突出，因此断裂带的展布对于地壳运动的水平速度场具有明确的控制作用，需要建立足够真实的有限元模型来反映研究区内断裂构造的复杂性以及相邻块体之间沿断裂可能发生的相对滑动。

2 模型和方法

2.1 唐山地区二维有限元接触模型

研究范围为 117.5 ～ 119.0°E、38.9 ～ 40.4°N内的矩形区域。根据邓启东［2］给出的唐山地区构造情况和断裂带展布，本文构建了1个二维矩形弹性薄片模型(图2)，主要涵盖Ⅰ-冀中块陷、Ⅱ-乐亭块陷、Ⅲ-山海关块隆、Ⅳ-蓟县块隆、Ⅴ-唐山块陷(图1(a))。本研究使用GPS观测结果作为边界条件进行约束，其时间区段为1999～2007年。由于模拟时间尺度较短且运动速率较小，地壳变形过程基本处于稳定状态。因此，在本研究中忽略动力学效应，将模型简化为静力学问题。本文采用二维弹性有限元模型模拟地壳水平运动，可以大大降低计算成本，将其作为精确三维地质体模型的合理近似，从而实现浅层岩石圈模拟［3］。

为了提高模型的可分析性，在保证模型客观反映基本问题的前提下，对研究区内的次级断裂带进行简化。纳入模型的断层共8条:(1)丰台-野鸡坨断裂、(2)蓟运河断裂、(3)唐山断裂、(4)宁河-昌黎断裂、(5)滦县-乐亭断裂、(6)沧西断裂、(7)滦县-卢龙断裂、(8)大李庄-古马断裂(图1(a))。因此，该模型可以满足唐山地区断裂带现今变形模拟分析的需求。

采用ANSYS有限元模拟环境，模型中构造块体和易变形带采用弹性连续本构关系，采用面单元对研究区进行单元剖分，见表1［4］。断裂带作为本研究中的重要因素，为了更加真实地反映断裂带附近破碎岩体较围岩强度更低的实际情况，本文通过介质物性参数弱化的方法在模型中构建了宽度约为断层长度1%的断裂过程带(图3)［5］。

图1 唐山地区地震构造、GPS站速度和有限元模拟结果Fig.1 Seismic tectonics，velocity by GPS observation and FEM simulation results in Tangshan area

图2 有限元模型Fig.2 FEM model

图3 断裂过程带示意图［5］Fig.3 Schech of the process zone

为了充分考虑相邻块体可能沿着断裂发生的相对滑动，将断裂带作为不连续面进行处理。使用conta175、targe169并以接触对的形式进行建模，即接触面之间用假想的法向和切向弹簧链接，其力学性质通过法向接触刚度、切向接触刚度和滑动摩擦系数来控制(表 2)［2，6］。

表1 Plane182单元的物理力学参数［4］Tab.1 Mechanical parameters of the unit no.182

表2 非接触单元的物理力学参数［2］Tab.2 Mechanical parameters of non-contact unit

2.2 GPS数据及边界条件

采用ITRF框架下中国大陆1999～2007年站速度结果，通过坐标旋转，得到区域无旋转坐标框架下GPS速度场(图1(b))。该结果反映了唐山地区地壳的水平运动状态，并可以清晰地观察断层两侧以及不同块体上的GPS速度特征。将靠近模型矩形区域边界处的GPS测点的实际观测速度作为距离最近处节点的位移约束，通过样条插值方法得到关键控制点之间其他节点的边界条件。

3 计算结果与讨论

3.1 地壳运动模拟结果及其信度检验

使用ANSYS计算得到的地壳运动模拟速度与GPS实际观测结果的对比如图1(b)。可以发现，GPS观测值和有限元模拟值的一致性较好，说明模拟结果可以较好地反映研究区地壳运动的实际。图4(a)、(b)分别为模拟结果以及与GPS观测结果的东西向分量、南北向分量的误差比对结果。可见，大部分测站的模拟速度都在GPS观测误差范围之内。东西方向分量和南北方向分量的总体标准偏差分别为1.52和1.34 mm/a，模拟结果信度较高，可以直接用于唐山地区断裂带现今形变特征研究。

图1(b)的地壳运动模拟结果显示，唐山地区地壳形变场具有较为明显的不均匀性和分界性，总体呈现出向西运动的速度场状态，研究区自北向南，除个别点外，速度矢量方向由正西向逐渐转变为南西向，不均匀性的同时，呈现较为规律的扇形分布。唐山块陷内部GPS观测点运动速率明显小于乐亭块陷、山海关块隆和蓟县块隆等外围地块，说明唐山块陷在四周深大断裂的围陷下，在1999～2007年间呈现为相对锁闭平衡状态。

3.2 断裂带运动速度分析

依据单元接触对的位移可以得到断层运动速率，包括平行于断层的水平走滑分量(右旋/左旋，右旋为正)和垂直断层面的水平拉张分量(拉张/挤压，拉张为正)，见表3和图5。

本文选择5条断裂带进行重点分析。运动速率计算结果显示:

1)丰台-野鸡坨断裂:该断裂主要以左旋挤压运动为主，其中左旋走滑速率为2.61 mm/a，挤压速率为1.39 mm/a。这与张四昌［7］使用地质构造、地震活动图像和震源机制解等资料综合分析得到的结论一致。该部位相对较小的挤压速率说明，丰台-野鸡坨断裂作为唐山块陷与蓟县块隆的主要边界，其挤压效应不甚明显，其变形情况与区域形变场具有相对较好的一致性，该断裂带活动性相对较弱。

2)蓟运河断裂:该断裂主要以左旋挤压运动为主，其中左旋走滑速率为6.31 mm/a，挤压速率为5.27 mm/a。该位置挤压速率较大，说明唐山块体在受到北东方向上的山海关块体的推挤下，在蓟运河断裂处受到阻碍，并进行快速的应变积累。该断裂带在现今构造环境下具有相对较大的活动性，地震危险性相对较大。

图4 地壳运动模拟结果和GPS观测结果的误差比对(误差棒代表GPS观测结果的误差范围)Fig.4 Comparision of crustal movement by GPS observation with FEM simulation

图5 唐山块体内断裂构造运动速度计算结果Fig.5 Movement rates of the faults in Tangshan block

3)唐山断裂:主要呈现右旋走滑和拉张状态，其右旋走滑速率为2.88 mm/a，拉张速率为0.21 mm/a。这与郭慧、江娃利等人［8］在唐山地区一系列钻探和孙家楼槽探得到的结果具有一致性。其微弱的拉张速率说明唐山块陷在锁闭状态下，块体内部的唐山断裂带所处的构造环境相对稳定，暂时不具备发生大震的可能性。

4)宁河-昌黎断裂:该断裂主要以右旋拉张状态为主，其走滑速率为1.25 mm/a，拉张速率为2.26 mm/a，说明乐亭块陷对唐山块陷具有明显的拉张牵引作用。滦县-乐亭断裂主要呈现为挤压逆冲状态，其挤压速率为1.4 mm/a，右旋运动速率为0.8 mm/a，说明山海关块隆对唐山块陷具有明显的挤压逆冲作用，这与虢顺民［1］的认识一致。

表3 唐山地区主要断裂带运动速率模拟结果及与地质结论比较Tab.3 Comparison of movement rates from FEM simulation with geological analysis

4 结论

本文利用1999～2007年GPS数据和二维接触有限元方法，对唐山地区断裂带运动进行数值模拟，计算了唐山块陷内部及其周边块体内相关断裂带的运动状态，得到以下主要结论:

1)总体而言，唐山地区地壳形变场具有较为明显的不均匀性和分界性，总体呈现出向西运动的速度场状态。研究区自北向南，除个别点以外，速度矢量方向由正西方向逐渐转变为南西方向，不均匀性的同时，呈现较为规律的扇形分布。唐山块陷内部GPS观测点运动速率明显小于乐亭块陷、山海关块隆和蓟县块隆等外围地块，说明唐山块陷在四周深大断裂的围陷下，在1999～2007年间呈现为相对锁闭平衡状态。

2)唐山断裂带在现今构造环境下的活动性相对较弱，不具备发生大震的危险性。但是山海关南西西向的推挤运动导致唐山块体与冀中块陷结合处的蓟县断裂呈现出快速的应变积累状态，其断层活动速率较高，具有发生较大规模地震事件的可能性。

1 虢顺民，李义志.唐山地震区域构造背景和发展模式的讨论［J］.地质科学，1977(10):305 －321.(Guo Shunmin，Li Yizhi.Discussion on the regional structural background and the seismogenic model of the Tangshan earthquake［J］.Scientia Geologica sinica，1977(10):305 －321)

2 邓起东，张培震，冉勇康，等.中国活动构造基本特征［J］.中国科学:D辑:地球科学，2002(12)1 020－1 030.(Deng Qidong，Zhang Peizhen，Ran Yongkang，et al.The basic characteristics of active tectonic structure［J］.Science in China:Series D，2002(12):1 020 －1 030)

3 雷显权，陈运平，赵炯洋.天山现今地壳变形的非连续接触模型模拟［J］.中南大学学报:自然科学版，2011(9):2 754 －2 762.(Lei Xianquan，Chen Yunping，Zhao Jiongyang.Discontinuous contact model to simulate current crustal deformation in Tianshan Mountain［J］.Journal of Central South University:Science and Technology，2011，42(9):2 754－2 762)

4 陈顒，黄庭芳.岩石物理学［M］.合肥:中国科学技术大学出版社，2009.(Chen Yong，Huang Tingfang.Rock physics［M］.Hefei:Press of University of Science and Technology of China，2009)

5 Vermilye J M，Scholz C H.The process zone:a microstructural view of faultgrowth［J］.Journal of Geophysical Research，1998，103(B6):12 223 –12 237.

6 安关峰，高大钊.接触面弹粘塑性本构关系研究［J］.土木工程学报，2001，34(1):88 － 91.(An Guanfeng，Gao Dazhao.Research on elastic-visco-plastic constitution of interfaces［J］.Civil Engineering Journal，2001，34(1):88 －91)

7 张四昌，刁桂琴.唐山地震序列的构造过程［J］.中国地震，1992(2):75 － 82.(Zhang Sichang，Diao Guiqin.The tectonic process of the Tangshan earthquake sequence［J］.Earthquake Research in China，1992，8(2):73 －80)

8 郭慧，江娃利，谢新生.钻孔与探槽揭示1976年河北唐山Ms7.8地震发震构造晚第四纪强震活动［J］.中国科学:地球科学，2011(7):1 009 －1 028.(Guo Hui，Jiang Wali，Xie Xinsheng.Late-Quaternary strong earthquakes on the seismogenic fault of the 1976 Ms7.8 Tangshan earthquake，Hebei，as revealed by drilling and trenching［J］.Sci China Earth Sci，2011，41(7):1 009 －1 028)

SIMULATING CURRENT CRUSTAL DEFORMATION IN TANGSHAN AREA WITH DISCONTINUOUS CONTACT MODEL CONSTRAINED BY GPS

Meng Qingxiao1，2，3)
1)First Crust Monitoring and Application Center，CEA，Tianjin 300180
2)Key Leboratory of Earthquake Geodesy，Institute of Seismology，CEA，Wuhan 430071
3)Wuhan Base of Institute of Crustal Dynamics，CEA，Wuhan430071

The current crustal deformation of the Tangshan area was analyzed with FEA numerical method under constraint of GPS observation，using physical contact model to obtain the discontinuous motion of faults.The results of simulation show that the movement of the faults in Tangshan area is different from each others.The faults in Tangshan area are right lateral slipping for Tangshan Block is pushed by Shanhaiguan plate.Jiyunhe fault zone has accumulated higher stain，therefore seismic risk there is higher.

Tangshan area;GPS;FEM;contact model;fault zone.

P227

A

1671-5942(2014)05-0040-05

2013-09-12

地震科技星火基金青年项目(XH14067Y);中国地震局地震研究所所长基金项目(IS201266113);中国地震局监测预报司震情跟踪专项(2014020212)。

孟庆筱，男，1986年生，助理工程师，主要从事工程地震、水库诱发地震、数值模拟和构造形变分析等研究。E－mail:mqxgucas@163.com。