廖力 李平恩 刘盼 杨建思

摘要：从弹性动力学方程出发模拟了鲁甸地震在包谷垴—小河断裂地震自发破裂过程，探讨影响鲁甸地震破裂的因素。研究结果表明：鲁甸地震的左旋走滑的震源机制以及震级主要是受背景应力场的影响，断层滑移分布受到断层几何结构、水平应力场方向及相对大小的影响，非平面复杂断层几何结构是导致鲁甸地震复杂的滑动位移分布的原因。

关键词：鲁甸地震;自发破裂;曲线有限差分方法;包谷垴—小河断裂

中图分类号：P315.332文献标识码：A文章编号：1000-0666（2019）01-0024-09

0引言

2014年8月3日云南魯甸发生6.5级地震，震源深度约12km。此次地震震中位于包谷垴—小河断裂及昭通—鲁甸断裂附近（27.1°N，103.3°E），是青藏高原巴颜喀拉块体及其邻近地区最新地震活跃期内发生的一次中等强度地震（邓起东等，2014）。张勇等（2014）利用地震波形数据反演了鲁甸地震破裂过程，结果显示是一次共轭破裂地震事件，约70%的地震能量释放在近SN向的包谷垴—小河断裂上;主震震源机制表明此次地震为一次高倾角左旋走滑地震。张振国等（2014）依据张勇等（2014）反演的破裂过程模拟了鲁甸地震的强地面运动并进行了烈度预测。现场地质调查发现鲁甸地震的地表破裂均在包谷垴—小河断裂带（李西等，2014），显示此次地震的发震断裂为包谷垴—小河断裂，余震精定位（房立华等，2014）及无人机航测技术分析（于江等，2018）的结果也支持这一论断。因此，本研究拟从弹性动力学方程出发模拟鲁甸地震在包谷垴—小河断裂的地震自发破裂过程，探讨分析影响鲁甸地震滑动位移分布的主要因素。

1研究方法与模型构建

1.1数值方法

本文利用曲线有限差分方法模拟包谷垴—小河断裂（图1）上的鲁甸地震的动力学自发破裂过程。该方法采用分裂节点描述断层间断面，网

格沿着断层以及自由表面等不规则界面划分。通过重新定义曲面断层分裂节点物理量（Dayetal，2005;Dalguer，Day，2007），曲线坐标系中的弹性动力学方程可以与摩擦准则约束的物理变量牵引力结合起来。该方法既保持了传统有限差分方法的特点，如计算效率高、容易施加等，又在模拟复杂断层方面有较高灵活性，可以用来模拟现实中常见的复杂几何断层，如倾斜、弯曲、跳跃等断层（张伟，2006;张振国，2014）。

1.2断层结构

包谷垴—小河断裂东起昭通—鲁甸断裂龙头山段，向西延伸至莲峰断裂，平均走向N30°W（李克昌等，1981）。震后地质调查表明，该断裂断层产状为60°∠60°（李西等，2014），断层长度约40km，图1给出了包谷垴—小河断裂的不规则地表迹线示意图。

现场地质调查只能给出断层在地表的出露，对于其深部构造则不能给出更多约束信息。断层的三维复杂几何结构一般由观测资料反演获得，特别是近场强震记录、空间观测技术等资料。综合这些因素，本文构造基础断层三维结构，即断层在地表的走向参考现场地质调查给出的地表迹线进行约束，深部则按照倾角为90°的断层向下延拓，整个断层模型沿走向方向长约60km，沿倾向方向宽20km。

现有的观测资料表明，断层两侧岩体的接触面并非简单平面，而是存在不同空间尺度的粗糙表面，数值计算表明粗糙表面会对自发破裂以及地震波辐射产生重要影响。因此笔者在基础三维断层面结构基础之上附加自相似扰动，扰动尺度约为1km。

地震，尤其是破裂面积较大的地震，其发生一般是由于一定区域达到破裂条件开始破裂并向外传播。数值模拟地震的破裂过程也需要一个成核区，在成核区内人为设置一个高于断层强度的剪切应力值，使之处于起始破裂状态，之后引起断层面上其他区域的自发破裂。根据地震精定位结果（房立华等，2014），鲁甸地震的起始破裂区域发生在断层模型沿走向约20km处，参考Xie等（2015）的震源机制解及张勇等（2014）的震源破裂过程及精定位结果，震源断层参数设定为断层走向165°，倾角90°，震源矩心深度3km，地震自发破裂深度设定为6km，震源成核区在龙头山附近。

1.3应力状态与动力学参数

地质观测资料以及研究区域的中小地震震源机制反演结果显示研究区域的断层以走滑断层为主。基于云南地区地震震源机制分析，孙业君等（2017）给出了该区域内三轴构造应力分布与相对大小。根据其研究结果，本文中包谷垴—小河断裂带的应力比设定为：

式中：σH，σh分别是最大、最小水平主压应力;σv为垂直主压应力，各应力基本关系表达式为：

式中：ρ和ρw分别表示上覆岩体和水的密度;h表示质点到地表的深度;g是重力常数。断层及周围岩体在构造应力的驱动下发生形变及移动，当形变积累的应力能力超过最大强度时，断层发生相对错动，释放应变能，产生地震波辐射。本文选取大尺度范围内的空间观测结果作为数值模拟中的应力场方向，并设定最大水平主压应力σH的方位角为330°，其他两个主压应力轴与σH相互垂直，三者构成右手系。图1给出了2个水平应力的方向和相对大小，其中箭头对分别表示最大水平主压应力σH和最小水平主压应力σh，σv垂直于水平面。

本文用简单速度弱化准则（Slip-weakeninglaw）作为摩擦准则（Ida，1972;Andrews，1976），临界滑移距离Dc=0.4m，在小于临界滑移距离时，静摩擦系数逐渐减小直至临界滑移距离转变为动摩擦系数，也反映了断层面滑动过程中应力释放过程（图2）。刘博研和史保平（2011）利用原地应力测量数据（郭啟良等，2009）和简单倾斜断层模型分析汶川地震动力学参数，得到静摩擦系数为0.5，动摩擦系数为0.2，由于钻井深度在400m左右，该研究仅限于浅部地震有一定参考意义。本文设定静摩擦系数μs为0.4，动摩擦系数μd为0.24。

1.4速度结构

包谷垴—小河断裂地处青藏高原东南缘（邓起东等，2002），位于印度板块与欧亚板块碰撞带的东部，速度结构复杂，对地震破裂过程以及强地面运动产生重要调制作用。简单的速度模型不能反映真实地震发生时对近场造成的地面运动和地震灾害。本文采用的是区域三维速度结果（郑晨等，2016），对研究区域介质可以提供较好的约束。

2研究结果

包谷垴—小河断裂上鲁甸地震的模拟主要是非平面上断层的动力学模拟，即地震触发之后的传播、扩展到愈合等一系列的自发破裂过程。为了精细地反映出地震的破裂过程，动力学模拟采用60m网格。

图3模拟了鲁甸地震后15s内的自发破裂过程，本文截取的是从破裂开始到扩展，直至破裂接触地表6个比较典型的位移分布结果。从图中可以看出，破裂在断层上的传播具有较强的非对称性，破裂从成核区（五角星）开始向地表传播，在地表形成了破裂。图4显示地震矩能量释放主要集中在第1～10s，矩震级MW=6.56。自发破裂模拟的震源时间函数和滑动位移分布与张勇等（2014）利用近震、远震数据联合反演的结果也较为接近。

从图5可以看出，鲁甸地震自发破裂向深部传播过程中，由于受到趋于静水围压的三轴应力的作用，断层深部的有效剪切应力趋于零，滑动在传播到断层深部边界之前就逐渐停止而愈合。破裂在断层走向两侧是正常传播，地震的最终滑移量分布（图5）显示较强的非均匀性。最主要的是滑动位移从成核区向NW方向传播直至地表，显示出较强的左旋走滑性质，这与区域背景应力场的性质关系较为密切。

3讨论

3.1断层几何弯曲的影响

本文采用的断層信息是从房立华等（2014）研究图中取点，以断层N30°W（李克昌等，1981）插值得到，断层倾角以60°从地表向下拓展。因此，模型中的断层面几乎是一个平面，而实际断层面常为不规则几何弯曲结构，从而导致断层面上更为复杂的应力状态。为研究断层几何弯曲对自发破裂的影响，将模型中断层向ES方向延拓50km，按照断层的平均走向在包谷垴—小河断裂段设置一些几何弯曲，自发破裂的成核区设置在整个断层的中部（龙头山段），保持震源深度、模型的背景应力场等参数设置不变模拟鲁甸地震自发破裂过程。从图6结果可以看到，鲁甸地震自发破裂的方式和设置模型较为一致，都是从成核区向NW方向扩展，但破裂释放的能量相对小一点，震级达到约MW6.35。相比于前述模型能量释放比较均匀，在此模型中破裂能量释放更快更集中，破裂的时间主要集中于0～5s，这种自发破裂的特点可能与成核区的断层几何弯曲结构有关。图7则显示了鲁甸地震自发破裂在断层走向和倾向上的最终滑移量，可以看到，断层的最终滑移都是发生在走向上，这也显示了鲁甸地震左旋走滑的特点。

3.2背景应力场的影响

3.2.1最大主压应力方向的影响

孙业君等（2017）给出了该区域内最大主压应力的方向，但因其采用近于水平主应力中量值较大的主应力方向近似，因此存在最大可能数十度的误差。因此，本文模拟了将水平最大主压应力顺时针和逆时针分别偏转20°，即最大主压应力方向分别为310°，350°的情形，模拟结果如图8所示。从模拟结果可以发现在最大主压应力向逆时针旋转（从310°到350°）的过程中，自发破裂强度逐渐变大，地震震级逐渐变大，在地表的滑动位移也更大。这表明本文模型中逆时针旋转的构造应力更有利于断层上的整体破裂。

3.2.2三轴应力场大小相对值的影响

鲁甸地震震中区域位于以巧家为中心的三轴应力场大小相对值R从高向低变化的过渡带（孙业君等，2017），为研究三轴应力场大小相对值对鲁甸地震自发破裂的影响，本文还模拟了R值设定为0.8，0.6和0.4的滑动分布。R值越大，在走滑断层背景应力场表示水平向2个主应力大小越接近。从图9可以看出：当R=0.8时，无论在地表还是断层深部，滑动位移相比其他两种情况的模拟结果都更大，R=0.6时的滑动位移比R=0.4时稍大，震级也稍大。以上说明水平向2个主应力大小越接近时，地震破裂释放的能量越大。

3.2.3背景应力场的影响

将模型设定的背景应力场按以下方式进行设置，应力比为：

各应力基本关系表达式为：

在此应力场设置下，断层性质表现为正断层，其他参数设置保持不变，模拟鲁甸地震后0～8s的自发破裂过程，断层倾向滑动速率结果如图10所示。从图中可以看出，鲁甸地震自发破裂从震源成核区出发，向震源两侧扩展传播，直到遇到人工强边界而停止。因为断层几何结构比较均匀，因此破裂显示出比较明显的均匀分布特性。断层破裂在此应力场设置下与前述2个模型的破裂方式完全不同，破裂主要是在断层倾向上发生，显示出比较明显的正断层性质，断层的最终滑移分布见图11。由此可见背景应力场是影响自发破裂方式的重要因素。

3.3滑动摩擦系数的影响

将模型的背景应力场按3.2.3节的方式进行设置，但将动摩擦系数在断层走向80km以南区域设为0.24，以北区域设为0.27，研究动摩擦系数对滑动位移分布的影响。由模型采用的滑动弱化准则可知，高的动摩擦系数使断层面上的剪应力更高，使应力降更小，从而使最终断层滑移更小。将成核区设在沿断层走向90km处，模拟结果（图12）与3.2节成核区2侧均匀的位移分布不同，成核区以南的断层滑移更大，而成核区以北因为较大的动摩擦系数导致较小的应力降，使最终滑动位移较小。

4结论

本文以鲁甸地震为研究对象，考虑区域三维速度结构差异、区域背景应力场、地震自发破裂断裂带几何结构、建立能反映地表起伏的断层自发破裂有限差分模型。再以现场地质调查的滑动破裂结果为约束条件模拟鲁甸地震自发破裂过程，得到以下结论：

（1）本文模拟的鲁甸地震自发破裂过程是基于云南区域的背景应力场，鲁甸断层的几何资料以及最新三维速度结构建立模型进行模拟，以最新震源机制反演结果的矩心深度以及震源参数设定自发破裂的成核区，模拟的震源时间函数与张勇等（2014）的反演结果较为接近，自发破裂模拟的地表破裂约6km，与鲁甸地震后地质调查发现的地表破裂也较为吻合。

（2）研究區域的背景应力场三轴主应力相对大小决定了鲁甸地震发震震级、走滑型的震源机制，主应力的方向变化会影响滑动位移分布。

（3）背景应力场对地震自发破裂起决定性作用，而不同的应力场参数包括水平主应力方向的改变以及相对大小的改变，会影响地震破裂的传播、滑动位错的分布以及地震波的辐射。

（4）断层的几何结构会改变断层面上应力场的大小，从而控制断层自发破裂，特别是发震断层的几何结构、参数等存在较强的非均匀性时，不同的起始破裂区域往往能导致相差较大的地震。地震起始破裂区域的位置对地震破裂形态有重要作用，可能对相应强地面运动规律与地震灾害分布也有重要影响。

（5）在滑动弱化准则控制下的动摩擦系数会影响断层破裂引起的滑动位移分布，动摩擦系数越小的断层单元最终滑动位错越大。

感谢南方科技大学张振国助理教授、中国科技大学张文强博士对自发破裂代码修改的建议和意见。

参考文献：

邓起东，程绍平，马冀，等.2014.青藏高原地震活动特征及当前地震活动形势[J].地球物理学报，57（7）：2025-2042.

邓起东，张培震，冉勇康，等.2002.中国活动构造基本特征[J].中国科学，32（12）：1020-1030.

房立华，吴建平，王未来，等.2014.云南鲁甸MS6.5地震余震重定位及其发震构造[J].地震地质，36（4）：1173-1185.

郭啟良，王成虎，马洪生，等.2009.汶川MS8.0大震前后的水压致裂原地应力测量[J].地球物理学报，52（5）：1395-1401.

李克昌，候学英，赵维城，等.1981.滇东北地区地震地质特征[J].地震研究，4（1）：57-63.

李西，张建国，谢英情，等.2014.鲁甸MS6.5地震地表破坏及其与构造的关系[J].地震地质，36（4）：1280-1291.

刘博研，史保平.2011.MS8.0汶川地震断层的应力状态以及对余震危险性的影响[J].地球物理学报，54（4）：1002-1009.

孙业君，赵小艳，黄耘，等.2017.云南地区震源机制及应力场特征[J].地震地质，39（2）：390-407.

于江，张彦琪，李西，等.2018.无人机航测技术在2014年鲁甸MS6.5地震震区活动构造调查中的应用[J].地震研究，41（2）：166-172.

张伟.2006.含起伏地形的三维非均匀介质中地震波传播的有限差分算法以及在强地面运动模拟中的应用[D].北京：北京大学.

张勇，许力生，陈运泰，等.2014.2014年8月3日云南鲁甸MW6.1（MS6.5）地震破裂过程[J].地球物理学报，57（9）：3052-3059.

张振国，孙耀充，徐建宽，等.2014.2014年8月3日云南鲁甸地震强地面运动初步模拟及烈度预测[J].地球物理学报，57（9）：3038-3041.

张振国.2014.三维非平面断层破裂动力学研究[D].合肥：中国科学技术大学.

郑晨，丁志峰，宋晓东.2016.利用面波频散与接收函数联合反演青藏高原东南缘地壳上地幔速度结构[J].地球物理学报，59（9）：3223-3236.

AbstractInthisstudy，wesimulatetheprocessofspontaneousruptureoftheLudianearthquakeatBaogunaoXiaohefault，anddiscusstheinfluencingfactorsoftheLudianearthquakerupture.TheresultsshowthatthefocalmechanismoftheleftlateralstrikeslipandthemagnitudeoftheLudianearthquakearemainlyaffectedbythebackgroundstressfield.Thedistributionoffaultslipisaffectedbythechangeoffaultgeometry，directionandrelativesizeofhorizontalstressfield.NonplanarcomplexfaultgeometrystructurecausesthecomplexslidingdisplacementdistributionoftheLudianearthquake.

Keywords：Ludianearthquake;spontaneousrupture;curvefinitedifferencemethod;BaogunaoXiaohefault