0 引言

地震活动伴随着断层错动的同时, 使地下介质积累的部分能量以应变能的形式传递到周围, 继而改变了一定区域的介质应力状态。因此, 通过对地震导致区域应力场变化的研究, 可以分析地震对周围环境产生的动力学影响。通常, 一次地震的发生会伴随余震或后续破裂事件, 静态应力触发即是研究地震在近场区产生的静态应力变化(不随时间变化), 并根据变化的大小来分析主震对其后续余震的触发作用(King et al, 1994; Toda et al, 1998; Harris et al, 2000; 万永革, 2001; Freed, 2005; 王莹, 2010; 盛书中等, 2015)。

北京时间2017年8月9日7时27分, 新疆博尔塔拉州精河县发生M_W6.3地震, 震中位于44.27°N、82.89°E, 震源深度11km, 震中距精河县城约37km, 该地震是一次逆冲型为主的事件。据中国地震台网中心(http://www.ceic.ac.cn/ )统计, 在距震中20km以内区域人口总数不超过800人, 截至8月9日12时, 暂无人员死亡。但此次地震强度高, 破坏性大, 波及博乐市、伊犁州、伊宁市等地区, 造成了较大的经济损失。

新疆地区地震活动频繁, 发育有多条活动断层, 在距震中200km的范围内分布着不少于20条第四纪活动断层(沈军等, 2014; 徐锡伟, 2015)。在这样复杂的地质背景下, 大的地震活动会使周围区域或构造的地震活动性发生明显改变(Toda et al, 1998)。截至9月12日共记录到340次精河地震余震事件 ^① 。本文对精河地震产生的静态应力变化进行研究, 并结合余震的空间分布进一步考察地震静态应力触发的效果, 以期为大震后危险区的判定提供理论依据(万永革, 2001; 盛书中等, 2015)。

① 房立华, 2017, 私人通讯

1 研究方法与数据

静态库仑应力变化综合地反映正应力、剪应力变化的大小, 更直观地体现了静态应力场的变化情况, 因此, 本文计算了由精河地震造成的周围区域静态库仑应力变化。King等(1994)对库仑破裂应力变化ΔCFS定义如下

$ \Delta {\rm{CFS}} = \Delta \tau + \mu '\Delta {\sigma _{\rm{n}}} $

(1)

式中, Δτ为剪切应力变化(定义与滑动方向一致为正); μ'为等效摩擦系数, 与孔隙流体和断层面介质特性有关, 一般取值为0.2~0.8(Harris, 1998; Cotton et al, 1997), 本文简化地以花岗岩的等效摩擦系数0.4为例; Δσ_n为正应力变化(拉张为正, 压缩为负)。综合以上参数计算ΔCFS, 并对其空间分布进行初步的探讨。

本文中使用的震源模型为张勇^②的有限断层位错模型, 该模型以2km×2km的网格对主震断层进行划分, 共划分出180个子断层, 并分别给出了各子断层面的经度、纬度、深度、走向、倾角、滑动角及滑动量等, 其中, 走向、倾角、滑动角是定量, 分别为92°、60°、92°。利用Okada(1992)总结的根据有限断层位错量计算位移场和应变场的解析表达式, 由胡克定律计算得到应力张量, 结合库仑破裂应力的定义(式(1)), 最终求得某断层面上的静态库仑破裂应力变化ΔCFS。如果ΔCFS为正值, 表明此次地震对该断层滑动有触发作用; 相反, 若ΔCFS为负值, 则该断层滑动受到抑制。大量研究(Harris, 1998、2000;Freed, 2005)表明, 库仑应力变化达到0.01MPa时便会有效地促使后续地震的发生。因此, 对于本文计算的库仑应力变化达到触发阈值0.01MPa的区域, 则被认为是地震危险区。

② http://www.cea-igp.ac.cntpxw275885.html

由于库仑应力反映了接收平面上的正应力和剪应力的相互作用, 因此, 平面的空间位置(经度、纬度、深度)和几何参数(平面外法向)是计算过程中所需要的。本文在震中附近主要断层上的库仑应力变化计算中, 参考了沈军等(2014)和徐锡伟(2015)的新疆地区地质构造图, 选取晚更新世-全新世以来的19条活动断层, 并将这些断层作为接收面计算其上的静态库仑应力变化, 分析此次地震对周围断层的应力触发效果。表 1给出了接收断层的几何参数。

3 研究结果

根据库仑应力的计算步骤, 本文逐一讨论了精河地震产生的位移场、应力场、周围主要断层上的静态库仑应力变化以及主震对余震的触发作用。

3.1 精河地震产生的位移场

本文的震源模型为有限矩形源模型, 利用Okada(1992)总结的解析表达式计算空间位移场及其对空间的偏导数。图 1为地下11km处(震源深度)垂直向和水平向的位移分布图。由图 1可见, 距震中40km范围内垂直位移变化明显。震中南侧呈现隆升趋势, 北侧为大面积沉降, 且越靠近震中, 位移变化量越大, 这与此次地震的逆断性质相吻合。研究区域中垂直向位移最大值约为6cm(隆升), 最小值约为-5.7cm(沉降)。水平位移方向(箭头的指向)呈现震中南北两侧向中心汇聚、东西两侧向外扩散的特点。这是由于断层两盘相向运动带动南北两侧物质向震中汇聚, 继而震中过剩的物质沿东西向外排开。为了进一步研究此次地震对地表位移产生的影响, 将计算深度取为0km, 结果如图 2所示。

对比图 1、2发现, 垂直向与水平向的位移场分布大体一致, 地表隆升最大值约为6.6cm, 沉降最大值约为-1.8cm。但仍有2个局部区域存在明显差异:①地表垂向位移为正的区域向北扩张, 负值区则被压缩。这是由于发震断层为向南倾斜, 随着深度逐渐减小, 断层线越向北移动, 因此, 断层上盘呈现的隆升区域向北扩张。②震中附近南北两侧水平位移的向外扩散。由图 1、2中子断层的位置可见, 扩散区刚好位于断层上盘, 因此, 该现象可解释为, 断层南侧物质向震中聚集导致物质过剩, 而地表几乎没有垂直向下力的约束, 故过剩的物质会向周围散开, 呈现出图 2的位移分布情况。

3.2 精河地震产生的应力场

由3.1节得到的位移场对空间的偏导数获得应变场, 结合应力、应变本构关系得到应力场。由于各点的应力状态为1个二阶张量, 很难直观地展示其变化情况, 所以, 本节仅考虑水平方向的应力, 简要分析其变化情况, 这有助于理解精河地震的发生对周围产生的动力学影响(万永革等, 2015)。

得到的水平应力场为2×2的矩阵, 形式如下

$ \mathit{\boldsymbol{\sigma }}{\rm{ }} = \left[ {\begin{array}{*{20}{c}} {{\sigma _{xx}}}&{{\sigma _{xy}}}\\ {{\sigma _{yx}}}&{{\sigma _{yy}}} \end{array}} \right]{\rm{ }} $

(2)

通过以下2个公式计算 σ特征值——最大主应力σ_Hmax和最小主应力σ_Hmin。正值表示为拉应力, 负值表示为压应力。

$ {\sigma _{{H_{{\rm{max}}}}}} = \frac{1}{2}\left[ {{\sigma _{xx}} + {\sigma _{yy}} + \sqrt {{{({\sigma _{xx}} - {\sigma _{yy}})}^2} + 4\sigma _{xy}^2} } \right] $

(3)

$ {\sigma _{{H_{{\rm{min}}}}}} = \frac{1}{2}\left[ {{\sigma _{xx}} + {\sigma _{yy}} - \sqrt {{{({\sigma _{xx}} - {\sigma _{yy}})}^2} + 4\sigma _{xy}^2} } \right] $

(4)

由于应力主轴相互垂直, 求出σ_Hmax轴与x轴(北向)的夹角θ, 即可确定主轴方向。θ可由下式求得

$ \theta = \frac{1}{2}{\rm{ta}}{{\rm{n}}^{ - 1}}\left({\frac{{2{\sigma _{xy}}}}{{{\sigma _{xx}} - {\sigma _{yy}}}}} \right) $

(5)

面应力可用来判断在该平面应力场条件下物质的状态。当面应力大于0时, 即拉应力为优势应力, 表现为物质膨胀; 小于0时, 压应力占优势, 表现为物质被压缩。面应力σ_Hc计算公式如下

$ {\sigma _{{H_{\rm{c}}}}} = \sigma {H_{{\rm{max}}}} + {\sigma _{{H_{{\rm{min}}}}}} $

(6)

本节计算的平面应力场深度为11km, 主应力和面应力分布如图 3所示。计算结果为:研究区内最大主应力的最大值约为0.57MPa(拉应力), 最小主应力的最小值约为-0.26MPa(压应力); 面应力最大值约为0.36MPa, 最小值为-0.09MPa。由图 3可见, 面应力膨胀区和压缩区近似南北、东西对称分布, 结合3.1节水平位移方向展布可知, 断层错动对南北两侧产生“拖拽”作用, 继而物质受拉膨胀, 过剩物质沿EW向排出, 继而东西两侧物质受压收缩。膨胀区拉应力占主导, 故在图 3中表现为红色箭头长度大于黑色箭头, 且红色箭头方向更接近真实受力情况, 整体表现为以震中为中心向外辐射状拉张。压缩区压应力占主导, 即黑色箭头更长, 其方向接近该点的受力情况, 压应力的东西分量表示两侧物质受到压缩的程度, SN向分量则表示断层错动带动周围的运动趋势。

3.3 精河地震对周围主要断层静态库仑应力的影响

将3.2节中得到的应力场投影到某断层面上, 就可以计算出精河地震在该断层面上产生的库仑应力变化ΔCFS。本节研究的主要断层位于43°~46°N、81°~85°E区域范围内, 结合表 1中各接收断层的走向、倾角、滑动角及其空间位置(经度、纬度、深度), 便可得到该断层面上的正应力Δσ_n和剪应力Δτ, 再利用式(1)求解各断层面上的静态库仑应力变化ΔCFS。表 1给出了各接收断层上的ΔCFS最大值和最小值, 图 4显示了各接收断层面的库仑应力变化情况。

图 4中红色框标注的2处断裂为ΔCFS大于0.01MPa的区域, 其一为库松木契克山前断裂中段, 该处静态库仑应力变化达0.03MPa; 其二为四棵树-古尔图南断裂西段, 其静态库仑应力变化达0.02MPa。我们认为2个区域的断层达到了静态应力触发的阈值(0.01MPa), 即为地震危险区, 该处地震活动时间可能被提前。除此之外, 距震中最近的库松木契克山前断裂东段的静态库仑应力变化达到了负峰值-0.1MPa, 表示精河地震在该段产生了抑制破裂的效果, 降低了该处地震的发生概率。

在距震中200km范围内的这些主要断层中, 库仑应力变化大体呈现3种趋势:①由西向东递减, 包括托里断裂、四棵树-古尔图南断裂、亚马特断裂西北段、博罗科努断裂北-中段、依连哈尔比尕断裂、蒙马拉尔断裂、博尔博松断裂; ②由西向东递增, 包括巴尔雷克断裂东北段、达尔布特断裂、亚马特断裂东段、喀什河断裂中段、库松木契克山前断裂西段、库松木契克断裂、科古琴断裂; ③巴尔雷克断裂、巩留南断裂库仑应力变化则分别在-600、-200Pa的量级附近波动。

3.4 精河地震对余震触发的研究

近年来, 许多学者结合静态库仑应力变化来研究主震对余震的触发作用。大量研究结果表明, 大震后库仑应力加载区的余震在数量上占有绝对优势(盛书中等, 2015), Deng等(1997)发现1932~1992年间发生在南加州的5.0级以上余震中约有85%发生在库仑应力增加区(ΔCFS大于0)。同样地, 本节结合精河地震产生的库仑应力与余震^②分布, 进一步考察地震静态应力触发的效果。

首先, 对余震的深度分布进行统计(图 5), 得到51%的余震震源深度为8~12km, 23%的余震为4~8km。因此, 本节分别计算10、6km深度的库仑应力变化, 将深度为8~12km的余震投影在10km的平面内, 4~8km的余震投影在6km的平面内, 结果如图 6所示。

由图 6可见, 绝大多数余震分布在库仑应力加载区, 为了精确地计算出余震被触发的比例, 将余震破裂面作为接收面计算该面上的库仑应力变化。目前普遍认为, 余震与主震具有相同的震源机制, 因此, 将接收面的几何参数(走向、倾角、滑动角)设为与主震的一致。计算结果显示, 在10km深度处, 148个余震的库仑应力变化为正, 25个余震的为负, 即85.5%的余震处在库仑应力加载区; 在6km深度处, 82个余震的库仑应力变化为正, 11个余震的为负, 即87%的余震处在库仑应力加载区。这进而说明, 此次精河地震静态触发模型对余震分布有较好的解释。

4 结论与讨论

本文采用张勇^②的有限断层破裂模型, 基于Okada(1992)给出的计算断层滑动产生的位移变化量及其空间偏导数的解析表达式, 分别计算此次精河地震产生的位移场、应力场、周围主要断层上的静态库仑应力变化以及主震对余震的触发作用。得到以下认识。

(1) 精河地震产生的垂直位移场在不同深度上存在差异。由11km深度至地表, 隆升区域向北扩张, 反映了此次地震断层为向南倾的逆断型; 地表隆升最大值约为6.6cm, 沉降最大值约为1.8cm。水平位移场呈现南北侧向震中汇聚, 且越靠近震中, 位移量越大, 而震中东西侧呈现向外“流出”; 地表水平位移在距震中较近的上盘呈现向四周扩散的趋势。

(2) 精河地震产生的水平面应力场展布南北侧物质主要受到指向震中的拉张力作用; 东西两侧物质主要受到因震中过剩物质EW向排出产生的东西向挤压力作用。

(3) 精河地震的发生导致了周围地区库仑应力的明显改变。位于震中西侧约20km的库松木契克山前断裂中段和震中东北部约50km的四棵树-古尔图南断裂西段的库仑应力加载尤为明显, 均大于0.01MPa, 即这2处为地震危险区。

(4) 余震主要分布在震中西侧的库仑应力加载区。在深度为8~12km的余震事件中, 约有85.5%受到主震的触发作用; 在深度为4~8km的余震事件中, 约有87%受到主震的触发作用。

通过计算库仑应力变化来研究地震触发已得到广泛的关注, 但目前大多数计算库仑应力变化的方法并未考虑构造应力场的影响。事实上, 如果区域的初始应力状态远低于破裂水平, 那么即使很大的库仑应力增加, 也很难达到应力触发的效果; 如果区域已经接近临界破裂状态, 一个很小的库仑应力增加就可能会触发地震。然而, 区域的初始应力状态, 特别是应力量值(万永革等, 2006)很难精确确定。如果这些参数能得到精确估计, 就会使地震静态应力触发的研究有很大的进展。

研究中对接收断层面作了相关简化近似处理, 比如, 主要断层面中同一断层的倾角和滑动角默认为相同, 但实际上各分段间会存在一些几何差异; 余震震源机制的不确定性很大, 尤其是震级较小的余震, 将其与主震震源机制近似处理来研究静态应力触发会使计算结果产生一定误差。王莹(2011)等认为在断层几何参数中, 滑动角对库仑应力变化影响最大。今后的研究中, 应在计算库仑应力时更多考虑断层几何参数的不确定性。

致谢: 衷心感谢两位审稿专家提出了十分宝贵的意见。本文采用北京大学张勇教授的精河地震有限断层模型以及中国地震局地球物理研究所房立华研究员提供的余震精定位数据, 在此一并表示感谢。