三维数值模拟在华北地区现今构造变形分析中的应用

胡勐乾 邓志辉 陆远忠 宋 键路 雨 朱秀云 孙 锋

1)环境保护部核与辐射安全中心，北京 100082

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

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

4)青岛市地震局，青岛 266034

0 引言

华北地区地质构造演化过程复杂、构造变形独特、地震活动性强，历来是众多学者研究的焦点。不论在地质资料、活动断层的研究，或者是地球物理资料、GPS观测的密度和时间上，华北地区都已有众多研究成果(郭慧等，2011)。前人对华北地区构造变形的数值模拟研究取得了丰富的成果(邓志辉等，2011a，b)，但也存在一些科学问题需要进一步深入研究。

在数值模拟实验中，模拟结果的对比是工作的重要环节之一，其用来检验数值模型的合理性及数据的可靠性，这在很大程度上决定了数值模拟结果的可信度。遗憾的是，这个环节常常被人们所忽视，或者仅进行某一方面的结果对比就认定数值模型的可靠性。

本研究在综合分析华北地区地质构造、地球物理资料的基础上，利用高性能并行计算技术，建立华北地区的三维有限元模型，以高精度的GPS观测研究结果为运动边界，对华北地区的运动速度场、断层运动性质、应变场三方面进行数值模拟结果对比，从而验证华北地区数值模型的可靠性与合理性。

1 华北地区三维有限元模型及边界条件

华北地区(106°～124°E，31°～42°N)，位于燕山、阴山以南，秦岭、大别山以北，贺兰－六盘断裂带以东一直到黄海。华北地区是中国大陆地震活动最强的地区之一，地震强度和频度非常高，因此也是研究程度最高的地区之一。

本文数值模拟工作是基于并行版ANSYS数值模拟软件平台实施。

1.1 三维有限元模型

1.1.1 边界范围及分层

根据中国活动地块划分及活动断裂分布，结合地壳、上地幔的波速结构、地震活动和GPS资料等，确定模型的几何边界范围是 99.8°～121.4°E，27.9°～42.3°N，区域总面积约为3840000km2(图1)。模型的边界范围包括了华北活动地块的绝大部分以及周边青藏、西域、南华和东北亚地块的部分地区。模型在纵向上分为均匀、水平的4层，即上地壳、中地壳、下地壳以及部分上地幔，总厚度为100km。

图1 华北地区有限元模型地表边界及活动断裂分布Fig.1 Finite element model boundary and active faults in North China.

1.1.2 活动断裂

华北地区的断裂分布非常密集(图1)，为了更加准确地模拟活动断裂的运动特征，模型中加入已知的所有晚更新世以来活动的断裂。由于断层空间分布复杂，在建模过程中充分考虑了断层的走向及延伸的几何形态，对其局部的细节进行一定的简化(表1)。

表1 华北地区有限元模型的活动断裂Table1 Active faults of the finite element model in North China

续表1

1.1.3 活动地块及介质参数

中国大陆晚新生代和现代构造变形以活动地块运动为主要特征。活动地块边界构造活动强烈，绝大多数强震都发生在其边界的活动构造带上。中国大陆及邻区可以划分出6个Ⅰ级活动地块区，它们是:青藏、西域、南华、滇缅、华北和东北亚(张培震等，2003)。其中华北活动地块区由鄂尔多斯、华北平原和鲁东－黄海3个Ⅱ级活动地块构成。

根据张培震等(2003)对中国大陆活动地块的划分，模型的介质参数在每一层上划分为8个区域:鄂尔多斯地块、华北平原地块、鲁东－黄海地块、阿拉善地块、青藏地块、燕山地块、华南地块以及地块边界带部分(图2)。利用Huang等(2003)计算得到的1°×1°中国大陆地区地壳上地幔波速结构，对数据进行筛选、求平均值，计算得到上述8个区域分别在上地壳、中地壳、下地壳及上地幔4层上的弹性模量(表2)。在模型中将断层处理为软弱带，是地块介质强度的10%。

1.1.4 网格划分

利用ANSYS软件中的solid187三维实体单元对华北地区有限元模型进行网格划分，为保证模型网格密度和计算质量，单元最大边长为25km，模型共划分为416582个单元，582392个节点。图2和图3分别表示了华北地区有限元模型的表面网格和整体网格分布，从图中可看出有限元模型的网格密度、介质分区及断层分布。

1.2 边界条件

1.2.1 模型边界条件使用的GPS数据

图2 华北地区有限元表面网格分布图Fig.2 Outside mesh of the finite element model in North China.

表2 华北地区各区域弹性模量Table 2 Block modulus in North China

目前，GPS技术已广泛应用于大地测量、资源勘查、地壳运动等领域。利用GPS测量数据作为数值模型的边界条件已有众多研究者进行了尝试，并取得了一定的成果(刘孙君，2003;郑勇，2005;刘峡等，2006;曹建玲等，2009)。本研究使用的中国速度场为1999—2004年和2004—2007年研究区域及周边GPS站点相对于欧亚板块的运动速率。其中，中国速度场1999—2004年数据是由1999，2001，2004年3期的观测结果得到的运动速度，而中国速度场2004—2007年数据是由2004与2007年2期的观测结果得到的运动速度。为了得到较好的GPS插值结果，选取了研究区及其邻区共979个台站的数据进行插值计算，得到华北有限元模型的边界约束条件。

1.2.2 确定模型边界条件的施加方式

图3 华北地区有限元模型整体网格图Fig.3 Overall mesh of the finite element model in North China.

GPS台站的数据并不能直接作为边界条件施加在模型上，经模拟实验发现，这样做将会在施加位移节点处引起较大的应力集中，模拟结果与实际差别较大。因此，首先需要对华北及邻区GPS观测数据进行插值，并提取模型边界处的GPS插值结果作为模型的边界条件，模型的底面在垂向上固定约束。图4和图5分别显示了华北地区1999—2004年和2004—2007年2个时期GPS观测的地表水平速度和模型边界约束。此为地壳表层的GPS约束、而在深度上GPS约束如何变化，将通过确定一个评价标准，并在这个评价标准下判断各种约束施加方式优劣来确定深部边界条件。

1.2.2.1 模拟结果的评价标准

对于在深度上如何施加GPS约束的问题，将通过比较所有GPS站点数值模拟结果与观测结果的符合度来确定，并且涉及断层带弹性模量的取值。根据模拟结果与观测结果，求出其离散度P，公式如下:

其中:n为区域内GPS点的数量;VEi为第i个GPS点的EW向模拟速度;VNi为第i个GPS点的NS向模拟速度;UEi为第i个GPS点的EW向观测速度;UNi为第i个GPS点的NS向观测速度。如果求得的离散度P越小，则说明GPS符合度越大;反之亦然。

1.2.2.2 几种不同载荷施加方式的比较

在深度上GPS约束的施加方式，模拟实验假定了3种方式，分别为:1)仅在地表层边界施加GPS约束;2)沿深度GPS约束大小相同;3)沿深度GPS约束逐渐增大。根据这3种约束施加方式所求得的离散度P，然后在合理的GPS施加方式中，分别将断层带的弹性模量设为地块介质的1/10、1/5和1/2，比较结果见表3。

图4 1999—2004年GPS观测的地表水平速度和模型边界约束Fig.4 1999－2004 surface horizontal velocity and model boundary constraints from GPS.

图5 2004—2007年GPS观测的地表水平速度和模型边界约束Fig.5 2004－2007 surface horizontal velocity and model boundary constraints from GPS.

从表3可以看出，GPS施加方式对模拟结果的影响较大，而断层弹性模量变化的影响较小。结果显示，GPS约束沿深度大小相同计算得到的模拟结果符合度更高，这可能与模型为弹性本构模型有关。

综上所述，采用大地直角坐标系，在华北地区数值模拟计算中使用弹性本构关系，模型表面为自由面，模型的底面在法向方向上固定，切向上自由。模型外边界施加的GPS约束沿深度大小不变，并且断层弹性模量取周围介质的1/10。

表3 华北地区模型各计算方案及离散度Table3 Calculation plans and dispersion in the North China model

图6 1999—2004年数值模拟地表速度场与观测速度场对比图Fig.6 Contrast of the 1999－2004 velocity field between the numerical simulation and observation.

2 数值模拟结果对比研究

数值模拟结果与实际资料对比分析是数值模拟不可缺少的一个过程，这些资料的比对在很大程度上可以验证模型是否合理。在确定了华北地区数值模型以及边界条件的基础上，将从运动速度场、断层运动性质、应变场来进行数值模拟结果对比，从而验证数值模型。

2.1 运动速度场模拟结果分析对比

在GPS边界约束条件下，将模拟结果所对应的GPS站点运动速度与实际观测的GPS运动速度进行对比，可以评价物理模型和数值模拟计算的可靠性。图6和图7分别表示了1999—2004年和2004—2007年2个时期以GPS插值结果作为数值模型边界约束条件模拟得到的地表运动速度场和GPS站点的观测结果之间的对比。

图7 2004—2007年数值模拟地表速度场与观测速度场对比图Fig.7 Contrast of 2004－2007 velocity field between the numerical simulation and observation.

从图5，6中可以看出，数值模拟速度与GPS观测运动速度大部分比较吻合，反映了比较好的一致性。模拟结果显示，地壳运动速率的整体分布呈现由西向东、自北向南逐渐增大的特点。华北地区北部、燕山地块的地壳运动速率总体较小，＜2mm/a，而在华北地区中部、南部以及华南地块，运动速率逐渐增大。在燕山地块地壳运动速度方向总体为SW向，而在华北地区及华南地块地壳运动速度方向逐渐转为SE向。

另外，2004—2007年华北地区地壳运动速率大小普遍高于1999—2004年，反映了2004—2007年华北地区地壳运动速度加快。而地壳运动速率的方向整体并无太大变化，只是在华北地区中、南部以及华南地块，1999—2004年速率方向近于SE向，而2004—2007年速率方向转为SSE向。

在个别区域，模拟结果与GPS观测结果也存在一定差异，其中，差异较大的点主要分布在华北平原块体的东北部，其原因可能与2006年7月4日文安5.1级地震的孕育和发生有关，也可能是GPS观测误差所致。

根据各地块内部的GPS数据，分别统计了1999—2004年和2004—2007年各地块运动速度的平均GPS观测值和平均模拟值，并进行了比较。表4、图8，9分别显示了比较结果。

从统计数据及图表可以看出，大部分地块运动速度的平均观测值与平均模拟值相差不大，少数有较大差别，如1999—2004年的阿拉善地块，2004—2007年的鄂尔多斯地块，可能是由于块体内部的GPS站点稀疏且个别台站误差较大引起平均值变化巨大。而2004—2007年与1999—2004年相比，各地块内部无论是平均GPS观测值还是平均模拟值都提高了1mm/a左右，反映了2004—2007年华北地区地壳运动速度加快。

表4 各地块运动速度的平均观测值和平均模拟值Table4 Averaged observation and averaged simulation of velocity of blocks

图8 1999—2004年各地块运动速度的平均GPS观测值与平均模拟值对比图Fig.8 Contrast of 1999－2004 block motion velocity between averaged GPS observation and averaged simulation.

2.2 断层运动性质模拟结果分析对比

华北地区有限元模型在建模过程中考虑了区域内所有晚更新世以来的活动断裂，根据数值模拟的计算结果可以得出每条断层两侧的运动速率，进一步确定断层活动性质。将模拟所得的断层活动性质与地质调查得到的活动断裂运动性质进行比较，可在一定程度上验证数值模型的合理性。由于断层建模时都按垂直处理，并且模型边界条件中GPS速度的施加也是在水平方向上的，所以只考虑走滑性质的断层活动。

图9 2004—2007年各地块运动速度的平均GPS观测值与平均模拟值对比图Fig.9 Contrast of 2004－2007 block motion velocity between averaged GPS observation and averaged simulation.

表5 断层走滑运动性质模拟结果与地质调查结果比较Table5 Contrast of faults slip between simulation results and geological survey results

续表5

表5列出了断层运动性质模拟结果与实际地质调查结果的比较，表中只统计了实际地质调查结果中具有走滑性质的断层。根据数值模拟断层走滑运动性质的结果，华北地区表现为右旋走滑的断裂带有银川－吉兰泰断陷带、山西断陷带、郯庐断裂带和唐山－磁县断裂带，表现为左旋走滑的断裂带有河套断陷带、张家口－渤海断裂带、渭河断裂带和安阳－徐州断裂带。

模拟结果与地质调查结果不符合的数据在表格中以灰色底纹显示。统计中共有48条具有走滑性质的活动断裂，1999—2004年模拟结果共有43条断裂活动性质符合，符合率为89.6%;2004—2007年模拟结果共有45条断裂活动性质符合，符合率为93.8%。

从统计结果的符合率上来看，数值模拟的结果能较好地反映断层的运动性质，说明华北地区数值模型具有一定的合理性。

图10 华北地区剔除整体运动后各点的运动(杨国华等，2003)Fig.10 GPS points movement excluding whole movement in North China(after YANG Guo-hua et al.，2003).

图11 华北及邻区主应变图(张静华等，2004)Fig.11 Principal strain in North China and adjacent areas(after ZHANG Jing-hua et al.，2004).

杨国华等(2003)根据高精度GPS复测资料，给出了华北1999—2001年的水平位移场(图10)，燕山块体本身的运动性质并不完全一致，若以中间为界，东部以西向运动为主(偏北)，西部以北向运动为主;燕山构造带东段(张渤带)以左旋走滑运动为主，燕山构造带西段(阴山带)以拉张裂陷运动为主，银川带以右旋走滑为主，鄂尔多斯块体南界(秦岭－渭河带)以左旋走滑运动为主，山西断陷带以拉张型右旋走滑为主，其他块(带)的相对运动则不明显。

表6给出了利用1999—2004年数据进行数值模拟得出的断层运动性质与根据GPS水平位移场得出的断层运动性质之间的比较，不符合的在表格中以灰色底纹显示。可以看出，数值模拟得到的断层运动性质与实测GPS资料得出的结果吻合得较好。其中，山西断陷带的恒山断裂和五台山北麓断裂与山西断陷带整体右旋的运动性质不符，其原因可能是由于这2条断裂的走向与山西断陷带的总体走向有较大的夹角所致。值得注意的是，华山山前断裂模拟结果为左旋，与地质调查结果(右旋)不符，但杨国华等(2003)的结果认为，该断裂所属的秦岭－渭河断裂带以左旋走滑为主，其原因可能是由于华山山前断裂的运动性质发生了改变。

2.3 应变场模拟结果分析对比

华北地区数值模拟的结果显示，华北地区的水平应变场总体上比较一致:第1主应变的方向为NNW－SSE，第3主应变(最大压应变)的方向为NEE－SWW。这与前人利用震源机制解结果、GPS测量结果等资料反演得到的华北地区主应变场比较一致(张静华等，2004;张跃刚等，2006;杨国华等，2007)。

图11为张静华等(2004)得出的华北及邻区主应变图，图12和13分别为利用1999—2004年和2004—2007年GPS数值模拟得出的华北及邻区主应变方向。对比图12和13可以看出，2张图具有比较好的一致性。华北平原地块、鄂尔多斯地块和燕山地块第1主应变方向为NNW－SSE;鲁东－黄海地块第1主应变方向大致为近NS向。图11和12整体比较一致，说明这个时期华北地区主应变方向没有太大变化。

表6 杨国华等(2003)结果与数值模拟结果比较Table6 Contrast of YANG Guo-hua et al.，2003 results with simulation results

3 结论与展望

3.1 主要结论

本研究基于并行版ANSYS数值模拟软件平台，根据华北地区活动地块划分及活动断裂分布，综合地质、地球物理研究成果，建立了三维有限元模型，模拟了华北地区现今构造变形特征，主要认识如下:

图12 1999—2004年数值模拟华北及邻区主应变方向Fig.12 1999－2004 simulated principal strain direction in North China and adjacent areas.

图13 2004—2007年数值模拟华北及邻区主应变方向Fig.13 2004－2007 simulated principal strain direction in North China and adjacent areas.

(1)将模拟得到的断层运动性质与地质调查得到的断层运动性质进行比较显示，数值模拟结果能较好地反映断层的运动性质，说明华北地区地壳运动速率整体分布呈现由西向东、自北向南逐渐增大的特点。华北地区北部、燕山地块的地壳运动速率总体较小，而在华北地区中部、南部以及华南地块，运动速率逐渐增大。

(2)模型具有一定的合理性。另外，利用1999—2004年数据进行数值模拟得出的断层运动性质与根据高精度GPS复测资料反演水平位移场得出的断层运动性质进行比较，结果显示，二者也吻合得较好。

(3)数值模拟结果显示，华北地区的水平应变场总体上比较一致:第1主应变的方向为NNW－SEE，第3主应变(主压应变)的方向为NEE－SWW。这与前人利用震源机制结果、GPS测量结果等资料反演得到的华北地区主应变场比较一致。将数值模拟得出的华北及邻区主应变场与张静华等(2004)得出的华北及邻区主应变图进行比较，两者也具有比较好的一致性。

3.2 存在的问题及进一步研究方向

活动断裂在数值模拟的工作中起到非常重要的作用，虽然本文考虑了模拟区域内绝大部分的活动断裂，但活动断裂都作为垂直考虑。这直接导致了对于模拟活动断裂运动性质而言，只能判断其走滑运动性质并与实际地质调查结果对比，而不能根据模拟结果判断其运动性质为正断层或逆断层。相信建立更接近实际断层形态的三维模型，对于进一步的数值模拟工作有很大帮助。

另外，与时间相关的动力学过程也是非常值得研究的问题。如汶川地震对华北地区的影响，由于岩石的蠕变机制，对于较远区域的影响作用可能需要较长的时间，所以考虑岩石蠕变性的数值模型进行数值模拟，其结果更有意义。

曹建玲，石耀霖，张怀，等.2009.青藏高原GPS位移绕喜马拉雅东构造结顺时针旋转成因的数值模拟［J］.科学通报，54(2):224—234.

CAO Jian-ling，SHI Yao-lin，ZHANG Huai，et al.2009.Numerical simulation of GPS observed clockwise rotation around the eastern Himalayan Syntax in the Tibetan plateau［J］.Chinese Science Bulletin，54(2):224—234(in Chinese).

邓志辉，宋键，孙君秀，等.2011a.数值模拟方法在地震预测研究中应用的初步探讨(Ⅰ)［J］.地震地质，33(3):660—669.doi:10.3969/j.issn.0253－4967.2011.03.015.

DENG Zhi-hui，SONG Jian，SUN Jun-xiu，et al.2011a.Preliminary study on application of numerical simulation methods to earthquake prediction research(Ⅰ)［J］.Seismology and Geology，33(3):660—669(in Chinese).

邓志辉，胡勐乾，周斌，等.2011b.数值模拟方法在地震预测研究中应用的初步探讨(Ⅱ)［J］.地震地质，33(3):670—683.doi:10.3969/j.issn.0253－4967.2011.03.016.

DENG Zhi-hui，HU Meng-qian，ZHOU Bin，et al.2011b.Preliminary study on application of numerical simulation methods to earthquake prediction research(Ⅱ)［J］.Seismology and Geology，33(3):670—683(in Chinese).

郭慧，江娃利，谢新生.2011.对1976年河北唐山MS7.8地震地表破裂带展布及位移特征的新认识［J］.地震地质，33(3):506—524.doi:10.3969/j.issn.0253－4967.2011.03.002.

GUO Hui，JIANG Wa-li，XIE Xin-sheng.2011.New evidence for the distribution of surface rupture zone of the 1976 MS7.8 Tangshan earthquake［J］.Seismology and Geology，33(3):506—524(in Chinese).

刘孙君.2003.GPS约束下青藏高原地壳运动位移场模拟及应力应变分析［D］:［学位论文］.武汉:中国科学院测量与地球物理研究所.

LIU Sun-jun.2003.Numerical simulation for displacement of crustal movement in Tibetan plateau and stress，strainanalysis constrained by GPS data［D］.Master thesis.Institute of Geodesy and Geophysics，Chinese Academy of Sciences，Wuhan(in Chinese).

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

LIU Xia，FU Rong-shan，YANG Guo-hua，et al.2006.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).

杨国华，张风霜，韩月萍，等.2007.华北地区现今水平运动场的动态特征［J］.地震，27(2):1—8.

YANG Guo-hua，ZHANG Feng-shuang，HAN Yue-ping，et al.2007.Dynamic characteristics of nowadays horizontal movement field in North China region［J］.Earthquake，27(2):1—8(in Chinese).

张静华，李延兴，郭良迁，等.2004.用GPS测量结果研究华北现今构造形变场［J］.大地测量与地球动力学，24(3):40—46.

ZHANG Jing-hua，LI Yan-xing，GUO Liang-qian，et al.2004.Study on present-day deformation and strain field in North China by use of GPS data［J］.Crustal Deformation and Earthquake，24(3):40—46(in Chinese).

张培震，邓起东，张国民，等.2003.中国大陆的强震活动与活动地块［J］.中国科学(D辑)，33(增刊):12—20.

ZHANG Pei-zhen，DENG Qi-dong，ZHANG Guo-min，et al.2003.Active tectonics blocks and strong earthquakes in the continent of China［J］.Science in China(Ser D)，46(Suppl):13—24.

张跃刚，帅平，胡新康，等.2006.从GPS观测看华北地区的形变场演化［J］.大地测量与地球动力学，26(1):36—41.

ZHANG Yue-gang，SHUAI Ping，HU Xin-kang，et al.2006.Study on evolution of deformation field in North China area with GPS data［J］.Journal of Geodesy and Geodynamics，26(1):36—41(in Chinese).

郑勇.2005.青藏高原及中国地区大陆地壳岩石层变形演化动力学数值模拟［D］:［学位论文］.合肥:中国科学技术大学.

ZHENG Yong.2005.Dynamic simulation of lithospheric evolution from the Tibetan plateau and China mainland［D］.Doctoral dissertation.University of Science and Technology of China，Hefei(in Chinese).

Huang Zhongxian，Su Wei，Peng Yanju，et al.2003.Rayleigh wave tomography of China and adjacent regions［J］.Journal of Geophysical Research，108(B2):2073.