引言

随着数字地震观测技术的不断提高， 近场观测的发展和积累， 使得人们对震源破裂过程的细节有了更多地了解， 其复杂性已成为当前震源研究中普遍重视的课题． 对震源破裂过程的研究可以通过波形反演拟合等方法获得对震源尺度、 震级大小、 断层面上滑动分布的时空不均匀性、 位错过程等的认识， 即从运动学角度展现震源过程的复杂性． 通过对大量大震震例的总结分析， 震源运动学反演得到了长足的发展， 建立了越来越细化的运动学模型． 尽管运动学模型能很好地拟合记录波形， 却没有考虑应力变化因素， 不能从力学角度解释复杂的震源破裂过程的物理本质(张海明， 陈晓非， 2003)， 因此我们需要进行震源动力学破裂过程的研究． 震源动力学的研究主要采用数值模拟方法， 基于弹性动力学原理并结合描述地震断层破裂面相互作用的本构关系， 研究地震断层破裂的成核、 扩展传播和停止等过程， 这对于预测强地面运动和理解地震破裂发生、 传播和停止的物理机制起到了关键作用(Kostrov， Das， 1989； Scholz， 1989)． 对于一些大地震破裂过程的研究已经证实了我们对地震破裂传播和应力释放的理解， 也证明了摩擦作用在地震过程中起到的重要作用(Wald， Heaton， 1994； Beroza， Mikumo， 1996)． 大地震断层上的破裂传播是由断层上的摩擦性质所控制的， 摩擦控制着断裂的起始、 发展和愈合(Madariaga, Olsen， 2000)． 同时， 实验技术的进步使得岩石破裂和摩擦行为的实验研究有了重要进展，其中摩擦本构关系和失稳机制的研究成为研究的热点(Dieterich， 1979； Ruina，1983)． 进行震源动力学破裂过程研究， 需要精确有效的震源模拟数值方法， 进而有限差分方法(Madariaga， 1976； Andrews， 1976a)、 离散元方法(Dalguer et al， 2001a， b， 2002)、 有限元方法(Oglesby et al， 2000a, b； Aagaard et al， 2001)和边界积分方程方法(Fukuyama， Madariaga， 1995； Aochi et al， 2000a, b)等得到了长足发展和广泛应用．

1.   复杂的震源动力学破裂传播过程

震源动力学破裂传播过程是非常复杂的， 与断层的几何形态、 断层与自由表面的相互作用等均有着密切联系．

作为零阶近似， 地震断层可以用简单的平面断层来模拟． 但是， 地震断层的表面并不是平面形态， 而是有很多的分支、 连接、 跨跃等． 一般大的破坏性地震也并不是发生在简单的单个平面断层上， 而是发生在由若干个子断层组成的复杂非平面断层体系中， 如1992年Landers地震, 1999年集集地震, 2008年汶川地震等(Aochi， Fukuyama， 2002； Zhang et al， 2003， 2010； Shen et al， 2009)． 那么复杂的非平面几何断层形态对破裂传播过程有什么影响呢？ 20世纪90年代， Yamashita和Umeda 1994)通过研究平行断层间的相互作用得出， 动力学裂纹的相互作用在某些情况下可以形成凹凸体和障碍体．之后， Aochi2000(a，b)研究了分支断层对动力破裂过程的重要影响， 结果表明， 断层的复杂几何形态 造成了断层面上的应力非均匀性， 进而导致了动力学破裂传播的多样化． Aochi和Fukuyama(2002)成功模拟了1992年Landers地震的动力学破裂过程． 由于小的Kickapoo 分支断层的存在， 破裂产生了“跳跃”， 使得其传播方向发生了改变． Oglesby和Day(2001a，b)虽然没有进行实际的地震动模拟， 但他们用简单的动力学模型也说明了断层几何形态对于解释滑动及逆冲断层近场强震动中的重要作用． 汶川M_S8.0地震后， 研究者们对其运动学破裂过程进行了反演． 沈正康等研究结果表明， 地震破裂断层存在两个滑动剧烈的地区： 映秀、 北川的滑动量近10m; 南坝滑动量为5 m (Shen et al， 2009)． 经深入观察发现， 这些滑动剧烈的局部地区与断层的几何非均匀性(弯折、 跨跳等)有一致的对应关系． 陈晓非(2010)根据地震断层的几何非均匀性与等效平面地震断层的剪切破裂强度等动力学参数之间的关系， 建立平面断层模型， 模拟计算了汶川地震的自发破裂动力学过程， 较好地解释了断层滑动复杂非均匀分布与断层的几何非均匀性之间的对应关系， 并指出这种几何非均匀性是震源运动学反演中复杂滑动分布的重要影响因素之一． 可见， 复杂的非平面断层几何形态对断层破裂传播过程有重要影响． 如果我们要得到震源破裂过程更详细的信息， 就需要将断层几何形态的复杂性考虑在内．

一般而言， 大地震都会在地表形成破裂形迹， 且多为浅源地震． 断层与地表倾斜相交是断层破裂复杂性的重要来源之一， 是震源破裂过程研究中的一个不可忽视的因素． 地震动的很多重要特征都是由于断层与地表有一定的倾角所造成的， 如地震动的上盘效应等(Oglesby et al， 1998， 2000a，b; Mikumo et al, 2003)． 从震源辐射的能量表象定理出发， 当断层接近地表且倾角很小时， 断层与自由表面相互耦合作用的地震波总辐射能远远大于由断层本身滑动位移所释放出的辐射能量， 随着断层与自由表面之间的距离进一步增加， 这种耦合作用随之减小(李彦恒等, 2009)． 地震波辐射能量与断层动力学破裂过程直接相关， 进而浅源断层与自由表面的强烈耦合作用对于破裂过程有着重要影响． 此外， 从地表反射回来的地震波可能会改变断层面上的应力场， 进而改变断层面上的破裂运动和近场强地震波辐射(Oglesby et al， 2000a，b)． 因此， 考虑自由表面的影响， 研究半空间甚至更为接近实际情况的分层半空间介质模型的破裂行为具有重要意义(Archuleta， Frazier， 1978； Nielsen， 1998)． 当断层平行或垂直于地表时， 讨论自由表面的影响比较简单， 可充分利用断层的几何对称性． 例如垂直走滑断层(图1a)， 可采用镜像技术， 将断层的反射影像置于自由表面的另外一侧来考虑自由表面的影响(Archuleta， Day， 1980； Miyatake， 1980a， b； Tse， Rice， 1986； Mikumo et al， 1987； Quin， 1990)． 但这一技术并不适用于断层与地表斜交等几何不对称性的倾斜断层(图1b)．

图  1  自由表面对垂直走滑断层(a)及倾斜断层(b)的影响

Figure  1.  Effect of free surface on vertical strike-slip fault (a) and oblique fault (b)

下载: 全尺寸图片 幻灯片

2.   数值模拟方法

由于震源动力破裂过程的复杂性， 只有极少数简单模型下的特殊问题能够求得解析解(Kostrov， 1964， 1966)． 因此， 一般地震断层的动力学破裂问题需采用数值方法求解． 目前， 常用的数值方法包括： 有限差分方法(Andrews， 1976a， b)、 有限元方法(Day, 1977)、 边界积分方程方法(Das， Aki， 1977; Andrews， 1985; Fukuyama， Madariaga， 1995; Tada， Yamashita， 1997； Tada et al, 2000)和离散元方法(Dalguer et al, 2001a, b,2002)等．

有限元方法对解决连续介质力学问题非常有效， 这种方法很容易实现模拟介质的非均匀性、 几何非线性和物理非线性， 但对于非连续介质， 则需要构造特殊单元来反映这种不连续性． 有限元法(Day, 1977； Oglesby et al， 1998， 2000a， b； Aagaard， 1999； Aagaard et al， 2001)网格划分灵活， 可以处理任意几何形态断层和起伏地表问题． 利用这种方法， Oglesby等(a, b)研究了倾斜断层动力学； Duan和Oglesby(2006, 2007)研究了平行走滑断层和弯折断层． 刘启方等(2003)采用显式有限元方法和多次透射边界技术建立了一种完全显式解耦模拟断层动力破裂过程的有限元方法， 该方法既可以用于均匀成层的介质， 也可以模拟复杂场地中的断层破裂问题． 但是有限元法运算量巨大， 计算时间过长， 成本高， 难以精确模拟三维实际问题．

离散元方法是一种处理不连续介质问题的数值模拟方法， 其理论基础是结合不同本构关系的牛顿第二定律， 采用动态松弛法求解方程， 最初主要应用在岩土、 矿冶等领域． Dalguer等(2001a， b)将二维离散元方法引入到地震动力学破裂模拟中， 成功模拟了1999年台湾集集地震的动力学剪切破裂过程， 解释了自由表面的存在及断层的几何不对称性是造成上下盘近场强震动存在较大差别的原因． 之后， Dalguer等(2002)又利用三维离散元方法研究了2000年Tottori地震的动力学破裂过程， 从应力变化角度分析了此次地震的前震、 主震及余震的发震机制及空间分布特征．

Andrews(1976a)和Madariaga(1976)最先独立地将有限差分法引入到震源动力破裂传播问题中， 目前该方法已成为研究断层破裂传播问题的重要手段． 未改进的有限差分方法由于单元模型的形状主要局限于正方形(体)或长方形(体)， 因此仅能用于研究水平或垂直断层模型． 但为了拓展有限差分方法在倾斜断层中的应用， Zhang等(2006利用Pitarka (1999)提出的非均匀格点距的四阶交错网格有限差分方法对倾斜断层进行了研究， 通过改变矩形网格的长宽间距， 使得给定倾角的倾斜断层面恰好通过矩形的斜边来模拟计算倾斜断层的破裂问题． 该方法不再要求断层表面必须平行于有限差分法的网格， 而可以任意倾斜， 但主要应用于平面断层． Zhang和Chen(2006c)提出了基于任意曲线坐标网格的有限差分方法， 用于研究地壳内具有任意倾角的曲面断层的破裂传播动力学问题． 罗扬(2007)将该方法用于模拟二维倾斜断层的自相似破裂和自发破裂； 祝贺君(2008)做了进一步研究， 尝试了二维非平面断层的自发破裂模拟； 朱德瀚等(2010)则基于该方法提出了以断层破裂面为约束面的断层面贴体网格有限差分算法， 模拟计算具有任意倾角的三维地震断层破裂动力学问题．

边界积分方法最早由Das和Aki(1977)用于模拟断层破裂动力学. 该方法是一种半数值半解析方法， 它通过离散化， 将描述地震断层破裂面上物理量变化的积分方程转化为易于求解的线性方程组， 再结合摩擦准则求解出地震断层破裂扩展的动力学过程． 不同于有限差分和有限元方法， 边界积分方程方法将域问题转化为边界问题， 使得求解问题的空间维数降低了一维， 计算效率提高(所谓的域方法， 需要对所考虑的整个区域进行离散化求解， 对于三维问题， 需求解大型方程组， 计算量较大). 并且该方法在处理复杂非平面断层问题上有优势(Aochi et al， 2000a， b； Aochi， Fukuyama， 2002； Aochi， Madariaga， 2003； Fukuyama， Mikumo， 2006)． 但由于边界积分方程的建立强烈依赖于格林函数， 而弹性动力学方程的格林函数对于全空间介质模型才有简单的解析解， 故早期该方法的研究仅适用于全空间介质模型． 为了拓展边界积分方程方法的应用， 张海明等首次建立了基于半空间格林函数的边界积分方程求解断层上的自发破裂问题， 能够模拟计算均匀半空间中任意倾斜平面断层的自发破裂动力学问题(Zhang, Chen, 2006a,b; Chen, Zhang, 2006)． 但由于该方法需要知道格林函数的解析表达式， 所以难以求解非均匀介质和具有起伏地表的断层破裂动力学问题．

地震震源的模拟涉及解弹性动力学波传播的偏微分方程的数值解． 进行数值模拟的方法很多， 但在这些方法中， 没有所谓的“最好”， 因为每种方法都有它的优缺点， 而方法的选择则视具体的问题以及计算机的计算能力而定(Day et al， 2005)． 有限差分方法已经被广泛地应用到数值模拟中， 但缺点是该方法不适宜处理复杂的边界问题． 有限元方法可以自由描述复杂的断层几何形态， 然而该方法网格划分非常复杂， 由于模拟断层的动力破裂必须将断层的可能破裂区域都包含在计算区域内， 因此计算区域较大， 计算时间较长. 而边界积分方程方法的优点就是仅考虑边界， 使问题简化， 计算量相对减小， 但该方法是建立在弹性互易定理的基础上的， 所以原则上不宜处理非线性问题， 不能用于非均匀介质中． 不管哪种方法， 对于一个三维的动力学破裂问题研究， 最好采用高效合适的算法， 且一般需要采用并行计算技术进行数值计算来提高计算效率．

3.   滑动摩擦准则

震源动力学的数值模拟研究其实是建立弹性动力学方程， 并在设定的合理初始条件和边界条件下， 求解动力学方程， 获取断层滑动量、 滑动速率及应力分布等动力学参数的过程． 地震过程在力学上可以被视为断层由静摩擦转为动摩擦的过程. 一般而言， 断层面上的应力被认为是动力学模型中最基本的参量． 对于震源破裂过程的动力学模拟， 摩擦准则起着重要作用． 大地震断层上的破裂传播是由断层上的摩擦性质控制的， 摩擦控制着断裂的起始、 发展和愈合． 对于不同的地震， 震源的破裂过程也不相同， 理解摩擦准则对于破裂过程的控制作用能为我们更好地理解震源过程的多样性和非均匀性提供物理基础． 当断层破裂后， 产生两个新的破裂面， 在很高的围压下， 两个新的断层面之间开始相对摩擦滑动， 这是一个典型的动力接触问题， 需要建立合理的本构关系模型来描述断层两盘之间的相互作用力(刘启方等， 2003)．

目前， 较为常用的摩擦准则是滑动弱化准则(Ida， 1972； Andrews， 1976a， b； Day， 1982a， b)、 速率弱化准则(Carlson， Langer， 1989； Cochard， Madariaga， 1994； Aagaard et al， 2001)和速率-状态依赖摩擦准则(Dieterich, 1979; Scholz, 1998)．

滑动弱化模型(图2a)最早由Ida (1972)提出， 他假定内聚力可以用位移间断的函数表示. 在这个模型中， 内聚区的应力与破裂面之间的滑动距离有关． 岩石摩擦实验表明， 随着滑动量的增大， 物质性质变弱， 最终进入稳态滑动模式． 在滑动弱化准则中， 破裂尖端是应力集中的部位(Ohnaka, 1992)． 当断层上一点的剪应力超过屈服应力强度τ_p时， 该点开始破裂并发生滑动. 随着位错Δu的增加， 剪应力以线性方式下降， 直至位错距离达到临界滑动弱化距离D_c后， 剪应力下降到稳定的剩余应力水平τ_r． 滑动弱化模型有很好的物理基础. Dieterich (1979)岩石破裂的实验结果表明， 破裂点的应力降不是瞬时产生的， 而是要经过一定的滑动距离， 应力才能下降到动摩擦力. 这与滑动弱化模型中应力与滑动距离成反比的假定一致(刘启方等, 2003)． 该准则数学上形式非常简单， 便于数值操作， 因此在震源动力学的模拟中被广泛使用(Andrews, 1976a, b; Madariaga, 1976; Harris et al， 1991； Harris， Day， 1993)． 但若要描述断层行为的一个完整循环过程， 则还需要修改(Okubo, 1989)．

图  2  滑动摩擦准则. (a) 滑动弱化准则； (b) 速率弱化准则； (c) 速率-状态依赖摩擦准则.D_c为临界弱化距离； τ_p为屈服应力； τ₀为初始应力； τ_r为剩余应力； Δu为位错

Figure  2.  Slip friction law. (a) Slip-weakening law; (b) Slip rate-weakening law; (c) Rate and state dependent friction law. D_c is critical slip-weakening distance; τ_p is yield strength; τ₀ is initial stress; τ_r is residual stress; Δu is dislocation

下载: 全尺寸图片 幻灯片

在利用滑动弱化准则进行数值模拟时， 需要预先设定的一个非常重要的参数就是临界弱化距离D_c． Ide和Takeo (1997)提出用波形反演的时空滑动速率函数来评价滑动弱化距离． Guatteri和Spudich(2000)展示了由滑动弱化距离与破裂应力降之间的消长关系而造成的分辨率的限制． 最近， Mikumo等(2003提出了一个利用近断层地震波形来评价滑动弱化距离的简单方法． 他们假设应力破裂时间非常接近于滑动速率达到峰值的时间． 如果破裂速度不发生剧烈变化， 滑动弱化曲线在应力破裂时间没有梯度的突变的话， 那么这个假设是有效的. 这种方法使得我们可以将滑动弱化距离与破裂区分开来．

通过大量的实验室试验以及野外观测资料分析， 目前已经有很多关于脆性岩层中D_c的研究． Ide和Takeo (1997)通过地震波形反演的方法评价了1995年Kobe地震中断层上的应力变化， 结果表明， D_c在脆性岩层中的量级为十至几十厘米， 而在0—4 km范围内D_c大于1 m． Ohnaka (1992)提出在脆性岩层(4—12 km)范围内， D_c为一常数； 而在0—4 km和韧性岩层(12—20 km)中逐渐变大(图3)． Marone和Kilgore (1993)的研究发现， 当断层泥厚度较大时， D_c较大． 根据对集集地震的动力学模拟结果表明， D_c与最终的滑动量成正比( Zhang et al， 2003， 2010)． 大的D_c会阻碍地震成核过程， 那么在这样的区域， 余震活动也会减弱． 对集集地震的研究也证实了这一点， 具有大的D_c的区域， 很少有余震发生(Zhang et al， 200)． 但具体到实际的地震动力学模拟， 还应根据具体情况来设定合理的D_c．

图  3  临界弱化距离D_c与深度 的依赖关系(引自Scholz, 1998)

Figure  3.  Depth dependence of critical slip weakening distance D_c

下载: 全尺寸图片 幻灯片

实验结果还发现(Dieterich， 1979； Ruina， 1983)， 摩擦与滑动速率和状态也有关系， 这就涉及到另外两个摩擦准则： 滑动速率弱化准则和速率-状态依赖摩擦准则(Okubo, 1989)． 滑动速率弱化准则(图2b)是指当断层上一点的剪切应力超过屈服应力强度时， 该点开始破裂并滑动. 随着滑动速率的增加， 剪切应力以某种方式下降. 当滑动速率达到某一值时， 剪应力下降到一个稳定的值． 当速率降低时， 摩擦增加， 控制着滑动愈合的方式． 而速率-状态依赖模型(图2c)是指当断层上一点的剪切应力超过应力强度时， 该点开始破裂并滑动． 剪切应力是位错和滑动速率的函数. 一般讲， 随着位错的增加， 剪切应力下降. 当位错到达临界滑动弱化距离时， 剪应力下降到动摩擦力. 然后， 随着滑动速率的下降剪应力有一个强化的过程， 由动摩擦力上升到静摩擦力(Zhang et al， 2003)．

如上所述， 滑动弱化准则在应用到整个地震循环模拟时受限， 因为其本身并未包含断层失稳后再次闭锁、 强度回升的机制， 无法解释粘滑的重复发生． 而速率-状态依赖摩擦准则把摩擦系数的变化与滑动速率变化联系起来， 并允许断层面通过一个或多个状态变量来演变， 这符合Dieterich (1979)和Ruina (1983)的实验结果. 在稳定滑动阶段， 滑动速率为主导因素． 在模拟地震循环的不同阶段， 从震前蠕滑成核、 动力学不稳定扩展到断裂愈合、 破裂停止以及震后余滑等(Tse， Rice, 1986; Dieterich, 1992)有重要应用． Dieterich (1992)基于该准则用数值模拟方法研究了地震成核过程， 揭示了不均匀性情况下断层错动及摩擦滑动加速现象， 并对速率-状态依赖的摩擦本构关系进行了改进． David和Peter (1999)采用颗粒流模型， 结合速率-状态依赖摩擦本构关系模拟了地震的孕育、 发震现象． Kato和Hirasawa (1996)基于该准则研究了远场加载速率对成核尺寸及断层强度的影响． 何昌荣(2000)以速率-状态依赖摩擦定律为基础， 模拟了均匀和非均匀断层的地震成核过程与前兆现象， 并讨论了正应力的影响． 刘桂萍等(2004)采用Gomberg建立的受速率-状态依赖摩擦准则控制的数学物理模型， 模拟了大地震后由区域静应力场变化引起的区域触发地震现象． 王威等(2008)则采用速率-状态依赖的摩擦本构关系， 结合断层弹性介质模型， 讨论了远场应变速率对非均匀断层成核过程的影响， 系统研究了各物理量在断层失稳过程中的变化． 此外， 滑动速率-状态依赖摩擦定律还被用于描述地震速率和相关的地震模式的变化(Parsons et al, 2000; Toda et al, 2002)． 早期利用速率-状态依赖摩擦准则进行断层滑动的时空演化过程研究多是考虑均匀断层面上的摩擦性质， 而对于非均匀介质的研究则主要基于实验室中对参数a和b随深度的变化研究(马胜利等, 2003)． 基于岩石摩擦的速率-状态依赖摩擦定律， 速度弱化(a-b＜0)是断层失稳的必要条件， 速度弱化向速度强化(a-b>0)的转换控制着断层从不稳定滑动向稳定滑动的转换． 因此， 断层摩擦滑动的速度转换还涉及断层带上地震活动特征、 地震成核深度及慢地震的机制等重要问题． Dunham等(2011a， b)研究了弹性和强速率弱化摩擦准则控制下的平面和非平面断层的破裂动力学问题， 讨论了近断层非弹性变形对于阻碍大尺度I型破裂的影响．

Fukuyama和Madariaga (1998)研究了滑动弱化和速率弱化两种摩擦准则下的破裂行为， 发现破裂的起始主要由滑动弱化摩擦准则控制， 而愈合主要由速率弱化控制． 二者并非相互对立， 而是互为补充的． 需要强调的是， 如果单纯地研究一次地震的动力学破裂传播过程， 滑动弱化准则比较有效Okubo, 1989; Cocco, Bizzarri， 2002)； 然而， 如果我们强调的是断层在地震循环中的力学条件， 则速率-状态依赖准则能对破裂面上的应力状态提供更好地描述．

4.   讨论与结论

震源动力学研究不同于运动学反演， 主要从力学角度探讨地震破裂成核、 传播和停止等过程的物理机制． 宽频带数字地震观测技术的发展以及反演复杂震源特征能力的提高， 使得震源动力学破裂的复杂性逐渐被揭示． 非平面断层几何形态、 自由表面与埋深较浅断层间的耦合作用等均为造成复杂动力学破裂过程的原因． 当然还有其它因素， 如不均匀初始应力场分布和物质强度等． 为了讨论这些复杂来源对断层动力学破裂传播过程的影响， 大量的数值模拟方法被广泛应用， 比较常用的如有限元方法、 离散元方法、 有限差分方法及边界积分方程方法等． 这些方法各有其优缺点， 在实际应用中可根据具体条件选择． 动力学数值模拟往往计算量较大， 采用并行计算技术可节省时间． 断层摩擦的本构关系是了解震源过程和机制的重要基础， 滑动弱化和速率-状态依赖摩擦准则是震源动力学实验研究的焦点． 但由于实验室条件等的限制， 对岩石摩擦本构关系仍没有很好地理解， 特别是摩擦滑动微观过程研究的缺乏， 造成了实验结果向实际震例应用外推的困难， 进而对于滑动弱化距离的给定， 目前仍在进一步的探索中． 而利用动力学模型模拟断层的破裂过程时， 主要的困难还是初始参数的合理选择．

在震源动力学理论不断发展完善的基础上， 拓展其在地震预测和地震工程中的实际应用是关键． 由于对工程结构造成破坏的地震动频段远远超过1 Hz， 因此， 在模拟和预测地震动时需要考虑高频地震动的预测方法． 探讨震源高频辐射地震动的机制， 是震源力学研究的一个主要方向． 通过对震源动力学的深入研究， 可以为运动学模型提供更合理的物理基础， 为地震动特别是近断层强地面运动预测提供帮助．

非常感谢审稿专家对本文修改提出的指导性建议．