海原-六盘山断裂带现今地震危险性的数值模拟分析

蒋锋云，季灵运,2，赵 强

1.中国地震局第二监测中心，陕西 西安 710054；2.防灾科技学院，河北 三河 065201

0 引言

地震是在区域构造应力作用下，应力在活动断裂带上不断积累并达到极限状态，而后突发失稳破裂的结果(Scholz，1990，1998; 张培震等, 2003，2013)，弹性回跳理论表明，当断层处于震间期时，随着构造作用力的持续加载，断层面逐渐闭锁，断层两侧地壳运动在断层面及附近不断积累应变能，直至发生破裂。因此，对活动断裂震间期断层面闭锁特征，相关区域构造变形与应变积累的研究，是地震研究不可或缺的内容。基于此，许多学者利用GPS观测得到的地壳水平运动速度场，对海原-六盘山断裂带构造变形、应变积累、断层面滑动分布及闭锁特征进行了深入的研究。这些研究对认识海原-六盘山断裂带现今构造变形特征以及强震发生的危险性，理解青藏高原东北部扩展变形机理具有重要的意义。然而这些研究大多从运动学的角度出发，很少涉及地壳介质属性和断层力学性质参数等动力学方面的研究和讨论。实际上，区域地壳形变本身受到这些参数的影响和制约，断裂带介质力学参数与断层面上的滑动速率、应力状态、区域地壳运动速度场等密切相关。如果不考虑这些力学参数对断裂带运动的影响，模拟结果可能和真实情况会存在出入。如Mildon et al.(2019)研究表明，库仑预应力场的存在和断层复杂结构对地震触发和强震危险性研究具有不可忽视的作用。也有学者借助有限元技术，顾及地壳介质和断裂带力学参数对包括海原-六盘山断裂带在内的区域进行形变场的模拟，并讨论了断裂带力学参数和地壳介质属性对区域地壳形变场的影响，给出了合理的参数解释。如He et al.(2013)利用GPS资料，借助有限元技术探讨了断裂带摩擦系数和地壳介质流变学参数对区域地壳速度场地的影响。研究表明断层低摩擦系数和恰当的地壳流变学参数是海原-六盘山区域现今地壳运动速度图像控制的主要因素。石富强等(2018)同样以GPS资料作为约束，采用三维有限元技术，模拟研究了青藏高原东北缘断层剪切力学性能(剪切模量)对区域地壳运动速度场图像的控制作用，进而在最优模型的基础上分析了当前青藏高原东北缘不同断裂的应力状态。结果表明，不同的断裂其剪切力学性能差异较大，对区域速度场的控制作用也不尽相同，解释了不同学者给出断裂摩擦性能存在差异的原因。这些力学模拟，对认识区域地壳介质属性和断裂带力学性质参数具有重要的意义。需要指出的是，由于各自研究侧重点的不同，上述研究将某一断裂带或断裂带的某一段的力学参数(摩擦系数或者剪切模量)由浅至深仅仅作为一个参数，并没有考虑到沿着深度方向上断层力学参数的变化。实际上断层深浅部力学参数可能存在较大的差异，如按照负位错理论，离开断层的远场变形主要由断层深部无震蠕滑(剪切滑动)来控制，而断层近场形变主要由断层浅部断层震间期运动来控制(部分闭锁；Reinoza et al., 2015)。

因此，文中将以海原-六盘山断裂带附近GPS观测地壳水平运动速度场作为约束，通过构建三维弹性有限元模型，借助非线性反演方法遗传算法，参考石富强等(2018a)的文章，将断层看做正交各向异性介质构成弱化带，来反演海原-六盘山断裂带沿着断层面平行方向上剪切模量深浅部的差异。进一步结合区域构造特征及地震活动分析海原-六盘山断裂现今地震危险性。

1 模型构建

海原-六盘山断裂带位于青藏高原东北部，该区域处于青藏高原北东向扩展的前沿地带，新生代构造变形和地震活动十分强烈。海原-六盘山断裂为该区域内主干断裂带，是青藏高原东北部的地貌边界与构造边界，也是青藏高原块体与阿拉善块体、鄂尔多斯块体的一级块体活动边界带(郑文俊等, 2019；李锦轶等，2019)。该断裂带自西向东由近东西走向的冷龙岭断裂、金强河断裂、毛毛山-老虎山断裂、狭义海原断裂和近南北走向的六盘山断裂组成。从地震活动性来看，广义海原断裂及周边地区历史强震非常活跃，曾发生过1920年海原8.5级和1927年古浪8.0级两次8级以上地震以及多次7级左右地震。两次8级以上地震之间长达220 km的天祝破裂空段(Gaudemer et al., 1995)有历史地震记载以来，没有发生过6级以上地震，存在强震发生的危险。六盘山断裂带中段也存在地震空区或者破裂空段，断裂带北段及其附近曾发生1219年固原7级地震、1306年固原7级地震和1622年固原北7级地震等；断裂带南段是公元600年天水—陇县间地震的发震断裂(王师迪等, 2018)；断裂带中段的隆德—陇县之间，历史上无强震与大地震记载，属于地震空区或者破裂空段，长约70 km，也具有强震发生的危险性(图1)。

图1 研究区构造特征及历史强震

文中主要研究海原-六盘山弧形断裂带断层深浅部平行断层面方向剪切模量的差异，因此在模型构建时，为了减少周边次级断裂的影响，将模型范围限制于沿断裂两侧100 km的范围内，断层的展布主要参考徐锡伟等(2017)最新活断层研究结果，将海原-六盘山断裂带近似为连续折线展布。

通过对所调研的地质资料进行假设，构建海原-六盘山断裂带三维有限元模型：考虑海原-六盘山地区地壳速度结构(陈九辉等, 2005)，以及地壳孕震层20 km厚深度(段星北, 1997; 张国民等, 2002)，将模型从地面延伸至30 km深。弹性力学参数参考Crust1.0(Laske et al., 2013)以及嘉世旭和张先康(2008)与张洪双等(2015)给出地壳速度结构，地壳弹性模量采用统一的值为85 GPa，泊松比为0.25。

1.1 有限元断层模型及网格划分

断层的模拟是文中研究的重点。相关研究表明，可以借助断层两盘之间的摩擦系数(Hergert and Heidbach, 2010)或者断裂带的剪切模量(Fialko et al., 2002)来模拟断层的力学性能。对应于有限元模拟中，通常采用接触摩擦模型，或者断裂带介质弱化模型来模拟断层的运动及变形。从长期构造演化来看(万年甚至百万年)，断层作为弱化带来处理是没有问题的。同样，接触摩擦模型如果来研究长期构造演化最后其实也是连续的滑动变形(时间步长大于地震复发周期)。但对于研究相对短期的(几十年)震间断层非连续变形而言，这二者之间就会存在差异。采用介质弱化模型来模拟断层短期变形行为，会使得断层内垂直断层走向力的平衡被破坏，而实际上在长期的构造作用下，断层与周边介质已经处于力平衡状态。因此介质弱化模型使得模拟结果会存在误差，而相对震间断层短期变形而言，忽略这种差异可能会使得模拟结果的合理性和可靠性存在问题。而接触摩擦模型，能够很好的反映断层两盘变形的非连续性，从接触摩擦理论的角度考虑，也具有明确的物理意义，更接近于实际情况。但摩擦接触计算是个非线性问题，不仅涉及到计算量大的问题，还涉及到解的收敛问题(李玉江等, 2009)。特别是文中如果考虑断层面深浅部摩擦系数的不同，不仅使得收敛变得困难，同时将其作为目标参数进行非线性反演，需要的计算量是非常巨大的，使得反演几乎不可能完成。

石富强等(2018b)在介质弱化带模型基础上引入正交各向异性理论，用正交各向异性本构方程替代传统弱化带中的线弹性本构方程，实现了弹性模量与剪切模量相互独立。即将断层带看作夹在块体之间的类似层合板的特殊介质，其在垂直于断层面方向剪切模量与周边介质一致，使之不易于发生大的挤压或拉张变形；而在平行断层面方向剪切模量相比周边介质要小，从而能够沿着断层面产生较大的滑动变形。这种正交各向异性模型，相比传统的介质弱化模型，一定程度上克服了由于降低弹性模量使得断层内垂直断层走向的力的平衡被破坏的情况，使得模拟结果和摩擦接触模型更为接近。同时，由于其是连续模型，又能避免摩擦接触非线性引起的计算结果的收敛和计算量过大的问题。因此，文中采用基于正交各向异性理论的介质弱化模型来模拟断断层的变形行为。根据复合材料力学基础知识，材料正交各向异性断层本构关系如下，断层内应力应变关系在其三个主方向上满足关系：

(1)

其中ε=[ε11ε22ε33ε12ε23ε13]T；σ=[σ11σ22σ33σ12σ23σ13]T分别为三个主方向的应变和应力分量；K柔度矩阵。断层介质的变形主要表现为走向和倾向的剪切行为，其他方向上的力学性质与周围介质一致，因此正交各项异性介质本构关系的柔度矩阵可以表述为：

(2)

其中E，v分别为周边介质弹性模量和泊松比；Gslip为断层在剪切方向上(走向和倾向)的剪切模量，独立于周边介质参数E，v。为了后面反演的需要，引入一无量纲参数λ=Gslip/G0，其中G0=E/2(1+v)为周边地壳介质剪切模量。断层的剪切方向上的剪切模量通常小于周边地壳介质剪切模量，所以λ通常介于0到1之间。

参照利用负位错进行断层面滑动分布反演时断层面的设置情况(李强等, 2014；Li et al., 2016；郝明等, 2017)，将断层沿走向和倾向，划分为不同的区，不同的区具有独立的平行断层面方向的剪切模量。海原-六盘山断裂断层浅部倾角较大，为了建模方便，将其均看作直立断层。垂直断层方向剪切模量和周围介质一致，而断层泊松比相比周围介质较大，取为0.28。考虑到文中以断层为中心构建模型，且模型尺度较小，GPS观测结果可能更多的反映了地壳浅部的运动，因此对模型进行格网划分的时候对于20 km以上，深度方向上采取较为密集的网格划分，每5 km划分一个网格，20 km以下采用较为稀疏的网格网划分每10 km划分一个网格。考虑到计算量和研究目的，在网格划分的时候，在水平面内也是在断层附近网格较密，而离开断层网格较稀疏。整个有限元数值模型采用SOLID186单元进行网格化，共划分3850个单元，18690个节点(图2)。

图2 有限元模型及网格划分

1.2 模型边界约束

在数值模拟研究中，模型的边界条件对最后的结果起着至关重要的作用，合理的边界条件是模拟能够取得成功的关键因素之一。从平面图(图3)上来看，文中模型边界区域均有一定数量的GPS观测地表水平运动速度场分布，将其附近GPS速度采用克里金差值方法内插到边界上作为模型水平向边界约束。GPS观测得到的地壳运动速度，是深部地壳运动在地表的直接反应，因此认为在文中给定的模型深度上，地壳运动速度沿深度方向不发生变化。不考虑深部复杂的地球动力学作用，约束模型底部垂向位移为0，水平向不作约束，可以自由移动。参照传统介质弱化模型断层模量的取值方法，在正交各向异性断层模型中，将平行断层面方向剪切模量和周围介质的比值的系数作为待确定参数，以模型内部断层两侧100 km以内GPS观测速度场作为约束，通过不断调整比例系数，使得观测值与拟合值的方差最小，从而得到最优的平行断层运动方向上断层的剪切模量的分布。GPS速度场采用Wang and Shen(2020)最新公开发表的1999—2016年长期速度场成果。目标函数如下：

图3 研究区GPS观测速度场(1999—2016年)

(3)

uix，vix，uiy，viy分别为x方向的模拟值和观测值以及y方向的模拟值与观测值。

1.3 反演方法

目前，基于限元软件的APDL语言与MATLAB进行交互式样访问的技术已经非常成熟(夸克工作室, 2002)，因此，对于文中的有限元非线性反演，可以方便的利用MATLAB提供的非线性反演工具遗传算法来实现。遗传算法GA(Genetic Algorithm)是根据生物进化思想而启发得出的一种全局优化算法，在本质上是一种不依赖具体问题的直接搜索方法。直接以目标函数值作为搜索信息，它仅仅使用适应度函数值来度量个体的优良程度，不涉及目标函数值求导求微分的过程，具有广泛的适应性。同时，遗传算法又具有群体搜索的特性，且基于概率规则，而不是确定性规则，具有获得全局最优解的特性。因此，其在各学科和领域中得到了广泛的应用。当然，它也有自身的缺点，比如编码不规范、搜索效率低，计算量大，过早收敛等问题。文中将通过调整遗传算法多个控制参数，扩大种群数和增加进化代数来获得全局最优解。通过反演得到的每一代的目标函数的最佳值和平均值图像(图4)可以看出经过多达300多次进化，目标函数最小值趋于稳定，且和平均值得差异很小，说明找到了全局最优解。

图4 目标函数收敛过程

2 模拟结果可靠性评价分析

对于模拟结果的检验一般追求的是模拟结果与观测结果一致，但实际上模拟结果和实际观测结果不可能完全一致，为了定量的给出模拟结果的可靠程度，石富强等(2018a)对He et al.(2013)的检验方法做了一定的补充之后，认为只要当模拟值与观测值之差位于观测结果误差椭圆之内时，模拟结果较为可靠，若超出误差椭圆，则可信度降低。将第i个观测点的观测结果与模拟结果差矢量的模与其对应方位误差椭圆极半径之比δi，作为第i个观测点的评价参数，当δi≤1时，i点的模拟结果落在该点观测值误差椭圆之类。具体计算公式如下：

(4)

(5)

(6)

从模拟结果和实际观测结果的对比来看(图5)，除了靠近断裂个别点之外，模拟结果和实际观测结果在速度大小和方向上基本一致。反映了在周缘稳定阿拉善地块、鄂尔多斯地块阻挡下，地壳运动沿着海原六盘山弧形断裂带呈顺时针旋转特征，构造上则表现为从海原断裂带的左旋走滑向六盘山构造带的逆冲推覆转换。从对模拟结果的残差统计分析上来看(图6)，接近90%站点的残差在1 mm以内，且绝大多数对于位于观测值误差椭圆之内，整个区域测点δi整体较小，算术平均值小于1，表明模拟结果与观测结果整体上具有较好的一致性。

图5 模拟速度场和GPS观测结果的对比

图中柱状灰度深浅代表文中公式(4)中δi的大小

应变作为微小量，不受参考框架的影响，能更好的反映出模拟结果观测结果之间是否具有一致性。因此，分别利用GPS站点观测值和模拟值，借助Shen et al.(2015)计算应变的程序，分别使用GPS站点上的模拟值和观测值，计算了研究区最大剪应变率和主应变率场(图7)。结果显示，模拟结果计算的应变率场空间分布特征和观测结果计算的应变率场基本上一致，只是在量值上略小，可能与模型没有考虑中下地壳介质粘弹性效应有关。从长期变形来看，介质的粘弹性效应对区域形变场的影响不可忽视，但由于文中模拟的是较短时间震间期地壳运动，且模型尺度较小，模拟残差在可以接受的范围内，因此可以近似认为模拟结果是可靠的。从实际观测值计算的最大剪应变率和主应变率场来看，研究区西段，应变高值区主要集中在毛毛山断裂、老虎山断裂和香山-天景山断裂之间，说明香山-天景山断裂和金强河断裂、毛毛山断裂、老虎山断裂一起承担鄂尔多斯地块和青藏地块之间在该区域的差异变形。模拟结果计算虽然也显示了这一特征，但不论是主压应变率的强度，还是最大剪应变率强度，相比实际观测计算结果都要小，可能与模型没有考虑香山-天景山等次级断裂有关。而考虑多条断裂更复杂的模型可能会对模拟结果得到改进，但同时会增加反演参数和计算量。狭义海原断裂带用模拟结果和实际观测结果分别计算应变率的差异要比金强河断裂、毛毛山断裂、老虎山断裂小得多，其高值区主要集中于断裂带附近，反映狭义海原断裂带两侧地壳差异运动明显。而六盘山断裂带利用模拟结果和观测结果分别计算的应变率场，在空间分布特征、量值上都有较好的一致性。从应变率图上来看，六盘山断裂带在其北段以东平凉附近有一明显的局部应变高值区，且主张应变和断裂近乎垂直，经过仔细分析认为，这一现象可能是由于位于平凉市内单个GPS站速度(GSPL)相对周围测站速度较大引起，与断裂的整体运动相关性不大。总体来看，六盘山断裂附近一定范围内应变不明显，可能反映六盘山断裂和汶川地震之前的龙门山断裂所处的情况类似，由于断裂高度闭锁，近场几乎没有相对变形。从主应变率来看，整个区域明显受到近北东向主压应变率和与垂直方向主张应变率的控制，且自西向东，主应变率方向总体上存在顺时针旋转特征。从主应变率和断裂的夹角来看，毛毛山断裂、老虎山断裂和狭义海原断裂带整体上以左旋运动为主，而六盘断裂带整体上，特别是相对断裂中远场，垂直断裂显示轻微的逆冲特征，和地质学上的认识较为一致。

图中色卡对应颜色表示最大剪应变率大小，十字叉表示主压、主张应变率大小和方向

3 讨论

在一个地震复发周期中，地震发生时，断层在短时间内发生快速滑动破裂，断层内阻碍断层运动的强度快速降低，即对应平行断层剪切模量快速减小，和周围介质剪切模量的比值也快速减小。随着地震离逝时间的增加，断层面逐渐愈合，断层力学属性逐渐增高，即对应平行断层走向剪切模量逐渐增大，和周围介质剪切模量的比值也逐渐增大，直到下次地震来临之前，断层力学属性又增加到较强的状态，即对应平行断层走向剪切模量几乎和周围介质相同，比值接近于1。因此，沿着断层面方向剪切模量和周边介质剪切模量的比值，一定程度上能够反映断层的闭锁特征。即比值越小表明断层越处于蠕滑状态，越容易滑动，越不容易积累应变能，比值越大表明断层更可能处于粘结状态，越不容易滑动，更容易积累应变能。因此可以用沿着断层面方向剪切模量和周边介质剪切模量的比值，即上文中所述的λ大小来表示断层面的闭锁程度。

从反演结果来看(图8)，六盘山断裂中南段闭锁程度较高，比例系数接近于1(s8和s9段)，深度最深达到30 km(s9段)，而其北段(s7段)和海原断裂交汇段，仅仅在浅部(0～5 km)存在一定程度的闭锁。这和利用负位错方法反演得到的六盘山断裂带断层面闭锁结果基本一致(李强等，2014; Li et al., 2016；郝明等，2017)，也和地震地质资料研究得出的六盘山断裂带中南段处于强震破裂空段具有一致性(M7专项工作组，2012)。整个狭义海原断裂带(s4—s6段)闭锁程度均较低，仅在浅部存在一定程度的闭锁，比例系数均在0.4以下，可能反映了1925年海原8.5级地震之后，整个断裂仍然处于震后调整状态。西段金强河断裂、毛毛山断裂、老虎山断裂浅部闭锁程度较低，而在5～20 km存在一定程度的闭锁，表明该断裂可能存在浅部蠕滑，而较深的位置存在能量积累，具有强震发生的危险性，这和相关研究的结论也基本吻合。

色块代表每个分块平行断层面的剪切模量和周围介质剪切模量的比值，白色圆圈代表距断层两侧各10 km范围内小震(3.0级以上)精定位结果在断层面上的投影

为了将数值模拟结果和小震精定位结果进行对比，利用中国地震局地球物理研究所房立华提供的小震重新定位目录，将1980以来发生在研究区内海原-六盘山断裂两侧各10 km范围内3.0级以上地震投影到断层面并叠加到数值模拟结果之上。结果显示，反演结果和地震活动沿着断层面空间分布特征具有较好的对应关系：六盘山断裂带中南段和北段浅部地震活动较为稀疏，和反演得到的断层高度闭锁相对应；海原断裂带整体地震活动稀疏，虽然反演结果显示断层闭锁程度较低，但可能与海原断裂带1925年8.5级地震使得整条断裂发生贯通性破裂有关(Gaudemer et al., 1995)；金强河断裂、毛毛山断裂(s1，s2段)5～20 km地震活动也较为稀疏，反演显示该区域具有一定闭锁程度，而老虎山断裂(s3段)由地表至深度25 km处，均显示地震活动较为密集，和其5～20 km存在一定程度的闭锁区相矛盾，根据老虎山断裂近场InSAR最新成果(Jolivet et al., 2012)，认为可能是反演所使用的GPS观测站点较为稀疏，空间分辨率不够所致。

虽然从模拟效果来看，文中的模拟结果在可接受的范围内，但从有限元模型构建的合理性、可靠性来看，仍然存在一些需要改进和深入的问题。①考虑更复杂的断层模型、地壳介质的粘弹性效应和初始应力场等方面的问题，会使模拟结果更好，也会提高模拟结果的科学性和合理性；②对于断层的模拟，选择什么样的模型与所研究的问题有关，如对于地质长时间尺度的模拟而言，传统的介质弱化模型是合理的，而对于短期的同震、震间或者震后形变而言，可能接触摩擦模型或基于正交各项异性的介质弱化模型更为合理；③仅以GPS观测的地壳水平运动速度场作为约束来反演断层剪切模量，缺乏InSAR、水准等垂向资料的约束，如果加入这些约束资料可能会增加反演可信度；④将断层统一看做直立断层，对于垂直于断层方向上倾滑运动(拉张或挤压分量)的模拟可能不够，特别是六盘山断裂，和龙门山断裂一样具有铲形特征，深部倾角较小。

此外，需要说明的是，相比基于块体运动的负位错模型而言，文中不需要考虑块划分对模拟结果的影响，但受限于反演方法带来的计算量和计算效率的问题，使得断层网格划分不宜过密，反演结果无法像负位错反演一样做到更高的分辨率，因此，文中反演得到的断层面网格划分较粗，反演结果不够精细。

4 结论

文中基于正交各向异性理论的介质弱化模型来模拟断层的变形行为，将平行断层面方向的剪切模量和周围介质剪切模量的比值作为反演参数，以海原-六盘山断裂附近现今GPS观测地壳水平运动速度场作为约束，通过构建三维有限元模型，采用遗传算法，反演了海原-六盘山断裂平行断层面的剪切模量分布。进一步结合研区域构造特征和断层附近地震精定位结果分析了海原-六盘山断裂带现今强震危险性。结果表明：六盘山断裂带中南段闭锁程度较高，沿断层面3.0级以上地震动活动稀疏，未来发生强震的危险性较高，而其北段仅在浅部有一定程度闭锁，深部3.0级以上地震活动相比中南段更为密集，未来发生强震的可能性较低；狭义海原断裂带整体闭锁程度不高，3.0级以上地震活动较为稀疏，反映在1920年海原8.5级地震之后，该断裂仍然处于震后调整状态，推测距下一次强震发生仍较为遥远；西段金强河断裂、毛毛山断裂和老虎山断裂在浅部存在一定程度的蠕滑，3.0级以上地震活动较为密集，但在5～20 km处具有一定的闭锁，并可与金强河断裂带、毛毛山断裂带的小震活动稀疏区域相对应，推测未来具有发生强震的背景。而老虎山断裂，虽然反演结果显示在5～20 km同样具有一定程度的闭锁，但沿断层面地震由地表至深部地震活动较为密集，可能处于贯通性蠕滑状态，未来发生强震的可能性不大。

致谢：中国地震局房立华研究员提供了研究区小震精定位结果，在此表示感谢。