刘博研 史保平 雷建设

1）中国北京 100085 中国地震局地壳应力研究所地壳动力学重点实验室

2）中国北京 100049 中国科学院大学地球科学学院

引言

前人对于动态和静态应力变化的研究表明，同震过程导致了断层的相互作用和地震的触发．大震后区域静态库仑应力变化常常被用于解释区域地震活动性速率的变化（Reasenberg，Simpson，1992；Simpson，Reasenberg，1994；Harris，Simpson，1998；Todaetal，1998）、主震断层外余震的发生（Dieterich，1972，1994；Rybicki，1973；Das，Scholz，1981；Steinetal，1992，1994；Dieterich，Kilgore，1996；Gross，Kisslinger，1997；Gross，Bürgmann，1998）以及即将失稳断层的地震发生概率的变化（Todaetal，1998）；而动态模型则被用于研究由动态应力变化引起的断层间相互作用和地震的触发（Harris，Day，1993；Gombergetal，1997，1998）．1992年美国Landers地震为我们提供了一些小范围内断层相互作用和远距离地震活动性触发的实例．King等（1994）计算得到Landers地震对后续Big BearMS6．5地震产生的库仑破裂应力增量约为0．2MPa，认为是Landers地震加速了Big Bear地震的发生．Spudich等（1995）还讨论了关于Big Bear余震是由1992年Landers地震触发的问题，并指出无论用动态触发还是静态触发，都很难解释这3．5个小时的时间延迟．Gomberg等（1998）解释了摩擦失稳模型中近距离（静态和动态）触发和远距离（动态）触发的主要特征，并讨论了静态和动态加载是怎样改变了断层的失稳时间．Harris（1998）则详细地总结了应力触发、应力影区及其对地震危险性的影响．

关于汶川地震应力触发的相关研究受到了国内外专家学者的重视．汶川地震震中位于龙门山断裂的中段，震源深度14km（中国地震台网中心，2013a），破裂沿着北东走滑、西部倾斜的龙门山逆冲断层带向前扩展．该地震引起地壳深部的岩石破裂长达300多千米（陈运泰，2008；王卫民等，2008）．Toda等（2008）利用几条主要活动断裂插值得到了汶川地震产生的应力变化分布．万永革等（2009）和单斌等（2009）计算了汶川地震导致的周围断层应力变化，解朝娣等（2010）在其基础上探讨了汶川地震引起的周边断裂应力变化后地震活动性随时间变化的情况．Parsons等（2008）还初步探讨了汶川地震产生的应力变化和对周围断层的影响．

汶川地震发生近5年后，2013年4月20日8点2分46秒，我国四川省雅安市芦山县（30．3°N，103．0°E）发生了MS7．0地震，震区发生大面积的滑坡和建筑物损毁．地震定位和震源机制结果显示该地震是发生在龙门山断裂带上的又一次逆冲型地震，也是汶川地震后发生在龙门山断裂带上最强的一次地震．芦山MS7．0地震发生在龙门山断裂带的南端，在汶川地震发生时这一地段没有发生破裂滑动．“芦山地震是否为汶川地震的余震”在学界引起了一定的争论．陈运泰等（2013）认为芦山地震是汶川地震的最强余震；王卫民等（2013）指出“芦山地震在宏观上可视为汶川地震一次‘迟到’的强余震”；而刘杰等（2013）认为两次地震的余震区存在约45km的间隔，芦山MS7．0地震不是汶川地震的余震．

然而，无论持有哪种观点，各位专家都不否认芦山MS7．0与汶川MS8．0地震有密切的联系．汶川地震发生后，雅安芦山地区的库仑应力增加的事实是肯定的（Parsonsetal，2008）．那么汶川地震究竟是怎样触发芦山地震？是否还会有类似芦山地震的强余震在周围地区发生？这才是研究汶川地震后区域库仑应力变化的最主要目的和最核心的问题．

本文利用含（滑移）速率和状态的摩擦定律（Dieterich，1994），结合汶川地震前后的地震活动性水平，深入阐述了芦山MS7．0地震的静态应力触发过程，从地震发生概率的变化上理解芦山地震的发震根源，并计算周边的鲜水河断层和熊坡断层的发震概率变化，从而为区域地震活动性以及危险性评价提供重要的参考．

1 含（滑移）速率和状态的摩擦定律

图1 弹簧－滑块模型示意图，由此得到的滑块运动方程为：τ＝－k（δ－δ0）．其中τ为剪切摩擦应力，k为弹簧的有效弹性系数，δ和δ0分别为滑动位移和初始滑动位移Fig．1 Spring-slide block model．The equation of motion for the slider is τ＝－k（δ－δ0），where k，δ，δ0，τare the spring constant，the displacement of the block，the displacement of the right-hand end of the spring，and the shear friction，respectively

含（滑移）速率和状态的摩擦本构关系为我们提供了一个对断层属性进行复杂的定量实验观测的基本框架（Dieterich，1979，1981，1994；Ruina，1983）．最简单地表述含（滑移）速率和状态的断层应力是由Ruina（1983）基于Dieterich（1994）的理论提出的．如果将其推广为多态变量（Rice，Gu，1983），则可以写作式中，τ和σn分别为剪切应力和正应力，˙δ为滑移速率，θi为状态变量．参数μ0，A和Bi为通过实验得到的系数．带有星号的项为标准化常量．Aln（˙δ／˙δ＊）表征了类似于小形变的不规则体或障碍体产生于滑动面的粘性阻力．Biln（θi／θ＊i）表征了与接触时间成正比的两表面间的化学附着状态（Kanamori，Brodsky，2004）．假设断层的摩擦过程由一维弹簧－滑块模型描述（图1），则剪切摩擦应力τ＝－k（δ－δ0）．其中k为弹簧的有效弹性系数，δ和δ0分别为滑动位移和初始滑动位移．若剪应力变化率˙τ为常数且不等于0，则其加载过程由τ（t）＝τ0＋˙τt描述，那么可得到滑动速率为（Kanamori，Brodsky，2004）

式中， 为初始滑动速率，H＝－k／σn＋（Bi／Dci）为模型参数（Dieterich，1994）．当地震过程中，瞬态滑移速率从每秒几厘米变化到每秒几米，或者长期滑动速率由每年几个毫米变到几个厘米时，可以看作瞬态滑移速率无穷大（1／˙δi＝0）．此时断层突然失稳，由此可得到失稳时间为

进一步，假设第n次子事件的失稳时间tf＝nΔt，其中Δt是第n－1次事件到第n次事件的时间间隔，那么可以得到第n次子事件的滑动速率

在初始时刻t0，由于主震使得周围断层突然有一个Δτ的应力加载，那么基于弹簧－滑块模型（图1），有kδ0＝τ0（t＜t0）和kδ0＝τ0＋Δτ（t＞t0）．为保证在t＝t0时滑动位移和滑动速率是连续的，有δ（t＝t0－ε）＝δ（t＝t0＋ε）和˙δ（t＝t0＋ε）＝˙δ（t＝t0－ε）exp［Δτ／（Aσn）］，即滑动速率增加了exp［Δτ／（Aσn）］倍．此时失稳时间变为

将式（4）代入式（5），可得地震（事件）次数

此时，n和tf（n）都是离散变量．假设很多类型的地震是通过背景场地震活动性的扰动所引起，而这个扰动又来自之前的地震造成的应力场状态的改变．将瞬态地震活动性速率R定义为R＝dn／dtf（n），并把n和tf（n）都看做连续函数．当参考剪切应力变化率˙τr与实际剪切应力变化率˙τ相等时，有

式中，r是在参考场或背景场下区域内剪切应力变化率˙τr下的常数稳态背景地震活动性；ta为余震的持续时间，即余震地震活动性恢复到背景场的特征时间，与剪应力的加载速率相关．根据Dieterich（1994）的理论，余震持续时间可定义为ta＝Aσn／˙τ．当t＜ta时，式（7）符合大森余震衰减定律；而当t＞ta时，地震活动性恢复到背景场r的水平．需要说明的是，式（7）中，R和r都为震级M的函数．

由于地震活动性与区域应力场有关，那么区域应力场的改变必将导致相应的地震危险性概率的变化．假定地震序列在时间t内，对震级M以上的发震概率服从泊松模型（Dieterich，Kilgore，1996）

那么，利用式（8）的概率分布，结合式（7）可以计算在应力Δτ的扰动下造成大于M级地震的发震概率．

2 发震概率计算

截至2013年2月28日，汶川地震已有M＞5．0的余震128次，M＞4．0的余震841次，M＞2．0的余震2 174次（中国地震台网中心，2013b），而且余震还在持续．图2显示了汶川地震前（1970-01-01—2008-05-11）后（2008-05-12—2010-05-11）的地震分布及M＞2．0地震的年发生率．震后汶川周边地区的应力状态发生了明显地变化（Parsonsetal，2008；Todaetal，2008；万永革等，2009；单斌等，2009；解朝娣等，2010），从而导致了地震活动性分布的显著不同．

图2 汶川地震前（a）、后（b）的M＞2．0地震分布及地震年发生率网格大小为0．05°×0．05°，右侧色标显示的年发生率为网格内M＞2．0的年发生率．红色圆圈代表M＞5．0地震，黑色圆圈代表2．0＜M＜5．0地震，灰色实线代表断层位置Fig．2 Distribution of M＞2．0earthquakes and seismicity rate of Longmenshan faultbefore（a）and after（b）Wenchuan earthquake Red and black circles represent M＞5．0and 2．0＜M＜5．0earthquakes，respectively．The seismicity on colorbar on the right is in the grid of 0．05°×0．05°．The gray line is the fault

假设汶川地区的余震持续时间为85年（申文豪等，2013），大于设定震级M的地震在该地区的年发生率为0．01次，A＝0．01，σn＝10MPa（Dieterich，Kilgore，1996）．在不同应力扰动的情况下，根据式（6）和（7）计算所得发生大于设定震级M的概率如图3所示．从图3可见，大应力扰动可导致地震发生率的快速上升，因此在受到大应力扰动的地区地震发生的时间会大幅提前．

图3 不同应力扰动下发生大于设定震级M的概率变化Fig．3 The probability of earthquake occurrence of magnitude＞M on different shear stress disturbance

对于雅安地区（30．0°—30．7°N，102．5°—103．5°E），汶川地震后库仑应力增加了0．1MPa（Parsonsetal，2008）．从1970年1月1日—2008年1月1日，共发生了23次M＞4．0地震，其中仅发生一次M＞6．0的地震．那么根据背景场地震活动性速率公式r≥M＝N／Δt（其中N为发生大于震级M的地震次数，Δt为时间间隔）可以求出r≥6．0＝0．026 3次／年．依然假设A＝0．01，σn＝10MPa（Dieterich，Kilgore，1996），那么在汶川地震后雅安地区的地震活动性概率变化如图4所示．从图4可以看出，发震断层的余震持续时间对于发震概率的影响并不大，而受到扰动之后的几年内发震概率的增幅较快．芦山地震发生在汶川地震后近5年，而由于受汶川地震的影响，雅安地区此时发生M6．0地震的概率已经增加到25%—30%．

图4 雅安地区发生M＞6．0地震的概率变化不同颜色曲线代表了不同的余震持续时间，黑色箭头所指位置为现今雅安地区对应的发震概率Fig．4 The probability of earthquake occurrence of M＞6．0in Ya＇an area Different colors give difference aftershock durations，the black arrow shows the earthquake probability in Ya＇an area now

事实上，汶川地震之后不仅雅安地区的库仑应力增加，其周边的鲜水河断层和熊坡断层也存在库仑应力的明显增加（Parsonsetal，2008），这些断层仍然是今后需要防范的重点地区．图5显示了汶川地震后鲜水河断层不同震级的发震概率变化．根据中国地震台网中心给出的地震目录，从汶川地震后到2013年2月28日，鲜水河断层已发生M＞4．0地震10次，其中2010年4月27日发生一次M5．0地震，此时M＞5．0的发震概率已超过80%，而至今未发生M＞6．0的地震，而M＞6．0的发震概率已约为60%，所以鲜水河断层附近将成为M＞6．0地震的重点防范地区．鲜水河断层未来发生M＞6．0地震概率由表1给出．

图5 汶川地震后鲜水河断层发生地震的概率变化实箭头所指位置为现今鲜水河断层对应的M＞6．0的发震概率，虚箭头所指的位置为2010年4月27日发生的M5．0地震对应的发震概率Fig．5 The probability of earthquake occurrence of Xianshuihe fault after Wenchuan earthquake The black arrow shows the probabilities of earthquake occurrence of M＞6．0on Xianshuihe fault at present，the dashed arrow points out the probabilities of earthquake occurrence when the M5．0earthquake happened on April 27，2010

表1 鲜水河断层未来发生M＞6．0的地震概率Table 1 The probabilities of earthquake occurrence of M＞6．0on Xianshuihe fault

图6显示了汶川地震后熊坡断层不同震级的发震概率变化．根据中国地震台网给出的地震目录，从汶川地震后到2013年2月28日，尚未发生M＞4．0地震，此时M＞4．0的发震概率已接近90%，虽然M＞5．0的发震概率仅为26%，但是该地区将来仍旧有发生中强地震的危险性．熊坡断层未来发生M＞4．0地震概率由表2给出．

图6 汶川地震后熊坡断层发生地震的概率变化黑色虚线所指的位置为现今对应的发震概率Fig．6 The probabilities of earthquake occurrence of Xiongpo fault after Wenchuan earthquake The dashed line shows the probabilities of earthquake occurrence at present

表2 熊坡断层未来发生M＞4．0地震概率Table 2 The probabilities of earthquake occurrence of M＞4．0on Xiongpo fault

3 讨论与结论

芦山MS7．0地震的发生重新引起了对汶川地震的热议．本文利用含（滑移）速率和状态的摩擦定律（Dieterich，1994），结合汶川地震前后的地震活动性水平，从地震发生概率的变化上解释了芦山地震的发震根源；计算得到了芦山地震发生时，雅安地区发生M＞6．0地震的概率已经增加了25%—30%．此外，本文对于库仑应力增加的鲜水河断层和熊坡断层也进行了发震概率的计算．根据中国地震台网中心给出的地震目录，从M＞4．0地震的发生情况来看，计算结果与实际情况基本符合．鲜水河断层从汶川地震后至今未发生M＞6．0地震，而M＞6．0的发震概率已约为60%；熊坡断层自汶川地震以来尚未发生M＞4．0地震，此时M＞4．0的发震概率已接近90%．所以，我们认为鲜水河断层附近将成为M＞6．0地震的重点防范地区，熊坡地区将来仍旧有发生中强地震的危险性．

然而，在计算过程中存在很多不确定性．首先，参数的选取可能带来一定误差，比如A取0．01，σn取10MPa（Dieterich，Kilgore，1996）等，与实际不完全一致；其次，库仑应力状态的计算过程依赖模型的选取，不同研究者可能给出不同的应力变化值（Parsonset al，2008；Todaetal，2008；单斌等，2009；万永革等，2009；解朝娣等，2010）．本文仅选用Parsons等（2008）结果进行计算也会带来一定的误差；最后，地震目录特别是小震目录的缺失使得对背景场地震活动性的估计不够准确．所以，应用条件概率模型会更具可靠性．但是由于缺乏发震断层的滑移历史、发震的准周期及上次地震至今的时间，条件概率模型在目前的研究中，尤其对我国西南地区来讲，仍然难以应用．

陈运泰．2008．汶川特大地震的震级和断层长度［J］．科技导报，26（10）：26-27．

陈运泰，杨智娴，张勇，刘超．2013．浅谈芦山地震［J］．地震学报，35（3）：285-295．

刘杰，易桂喜，张志伟，官致君，阮祥，龙锋，杜方．2013．2013年4月20日四川芦山M7．0级地震介绍［J］．地球物理学报，56（4）：1404-1407．

单斌，熊熊，郑勇，刁法启．2009．2008年5月12日MW7．9汶川地震导致的周边断层应力变化［J］．中国科学：D辑，39（5）：537-545．

申文豪，刘博研，史保平．2013．MW7．9汶川地震余震序列触发机制研究［J］．地震学报，35（4）：461-476．

万永革，沈正康，盛书中，徐晓枫．2009．2008年汶川大地震对周围断层的影响［J］．地震学报，31（2）：128-139．

王卫民，郝金来，姚振兴．2013．2013年4月20日四川芦山地震震源破裂过程反演初步结果［J］．地球物理学报，56（4）：1412-1417．

王卫民，赵连锋，李娟，姚振兴．2008．四川汶川8．0级地震震源过程［J］．地球物理学报，51（5）：1403-1410．

解朝娣，朱元清，Lei Xinglin，于海英，虎雄林．2010．MS8．0汶川地震产生的应力变化空间分布及其对地震活动性的影响［J］．中国科学：D辑，40（6）：688-698．

中国地震台网中心．2013a．中国地震台网（CSN）地震目录［EB／OL］．［2013-05-05］．http：∥www．csndmc．ac．cn／newweb／cgi-bin／csndmc／csn＿catalog＿p003．pl？mode＝catalog＆ot＝20080512＿0627595．

中国地震台网中心．2013b．中国地震台网（CSN）地震目录［EB／OL］．［2013-05-05］．http：∥www．csndmc．ac．cn／newweb／data／csn＿catalog＿p002．jsp．

Das S，Scholz C H．1981．Theory of time dependent rupture in the Earth［J］．JGeophysRes，86（B7）：6039-6051．

Dieterich J H．1972．Time dependence friction as a possible mechanism for aftershocks［J］．JGeophysRes，77（20）：3771-3781．

Dieterich J H．1979．Modeling of rock friction：1．Experimental results and constitutive equations［J］．JGeophysRes，84（B5）：2161-2168．

Dieterich J H．1981．Constitutive properties of faults with simulated gouge［C］∥MechanicalBehaviorofCrustalRocks：TheHandinVolume，GeophysicalMonograrphSeries．Washington D C：AGU：103-120．

Dieterich J H．1994．A constitutive law for rate of earthquake production and its application to earthquake clustering［J］．JGeophysRes，99（B2）：2601-2618．

Dieterich J H，Kilgore B D．1996．Implications of fault constitutive properties for earthquake prediction［J］．PNatlAcad SciUSA，93（9）：3787-3794．

Gomberg J，Beeler N M，Blampied L M，Bodin P．1998．Earthquake triggering by transient and static deformations［J］．JGeophysRes，103（B10）：24347-24358．

Gomberg J，Blampied L M，Beeler N M．1997．Transient triggering of near and distant earthquakes［J］．BullSeismolSoc Am，87（2）：294-309．

Gross S，Bürgmann R．1998．Rate and state of background stress estimated from the aftershocks of the 1989Loma Prieta，California，earthquake［J］．JGeophysRes，103（B3）：4915-4927．

Gross S，Kisslinger C．1997．Estimating tectonic stress rate and state with Landers aftershocks［J］．JGeophysRes，102（B4）：7603-7612．

Harris R A．1998．Introduction to special section：Stress triggers，stress shadows，and implications for seismic hazard［J］．JGeophysRes，103（B10）：24347-24358．

Harris R A，Day S M．1993．Dynamics of fault interaction：Parallel strike-slip fault［J］．JGeophysRes，98（B3）：4461-4472．

Harris R A，Simpson R．1998．Suppression of large earthquakes by stress shadows：A comparison of Coulomb and rateand-state failure［J］．JGeophysRes，103（B10）：24439-24451．

Kanamori H，Brodsky E E．2004．The physics of earthquakes［J］．RepProgPhys，67（8）：1429-1496．

King G C P，Stein R S，Lin J．1994．Static stress changes and the triggering of earthquakes［J］．BullSeismolSocAm，84（3）：935-953．

Parsons T，Ji C，Kirby E．2008．Stress changes from the 2008Wenchuan earthquake and increased hazard in the Sichuan basin［J］．Nature，454：509-510．

Reasenberg P A，Simpson R W．1992．Response of regional seismicity to the static stress change produced by the Loma Prieta earthquake［J］．Science，255（5052）：1687-1690．

Rice J R，Gu J C．1983．Earthquake aftereffects and triggered seismic phenomena［J］．PureApplGeophys，121（2）：187-219．Ruina A L．1983．Slip instability and state variable friction laws［J］．JGeophysRes，88（B12）：10359-10370．

Rybicki K．1973．Analysis of aftershocks on the bases of dislocation theory［J］．PhysEarthPlanetInter，7（4）：409-422．

Simpson R W，Reasenberg P A．1994．Earthquake-induced static stress changes on central California faults［J］．USGeol SurvProfPap，1550-F：55-89．

Spudich P，Steck L K，Hellweg M，Fletcher J B，Baker L M．1995．Transient stresses at Parkfield，California，produced by theM7．4Landers earthquake of June 28，1992：Observations from UPSAR dense seismograph array［J］．JGeophysRes，100（B1）：675-690．

Stein R S，King G C P，Lin J．1992．Change in failure stress on the southern San Andreas fault system caused by the 1992magnitude＝7．4Landers earthquake［J］．Science，258（5086）：1328-1332．

Stein R S，King G C P，Lin J．1994．Stress triggering of the 1994M＝6．7Northridge，California，earthquake by its predecessors［J］．Science，265（5177）：1432-1435．

Toda S，Lin J，Meghraoui M，Stein R S．2008．12May 2008M＝7．9Wenchuan，China，earthquake calculated to increase failure stress and seismicity rate on three major fault systems［J］．GeophysResLett，35（17）：L17305．

Toda S，Stein R S，Reasenberg P A，Dieterich J H，Yoshida A．1998．Stress transferred by the 1995MW＝6．9Kobe，Japan，shock：Effect on aftershocks and future earthquakes probabilities［J］．JGeophysRes，103（B10）：24543-24565．