阳江新湖水库ML3.0级以上地震应力触发效应的统计检验

王力伟 黄柳芳 吴国瑞

摘 要：为研究新湖水库有波形记录以来的10次ML3.0级以上地震之间是否有静态库仑应力触发现象。文章首先使用双差定位方法对这10次地震进行了精定位，然后使用CAP方法反演了2次ML4.0级左右地震的震源机制解和最佳震源深度，使用初动和振幅比方法反演了其他8个震级为ML3.0～ML3.7的小地震的震源机制解。最后，使用PSGRN/PSCMP软件包，计算并统计了不同计算参数组合下，10次地震之间的静态库仑应力触发结果。结果发现：新湖库区9次后发生的地震中，6次地震被已发生地震所触发，占66.7%，1次地震被已发生地震所抑制，占11.1%，2次地震与已发生地震无关系，占比22.2%。新湖库区10次ML3.0级地震之间，存在微弱的静态库仑应力触发关系。

关键词：新湖震群；震源机制；双差定位；库仑应力；统计检验

中图分类号：P315 文献标志码：A 文章编号：2095-2945（2018）25-0001-06

Abstract： In order to study whether there are static Coulomb stress triggering phenomena among 10 earthquakes with ML≥3.0 since the waveform records of Xinhu Reservoir. In this paper， the 10 earthquakes were first precisely located using the double difference location method， and then the focal mechanism solution and the optimum focal depth of the two earthquakes with ML=4.0 or so were inversed by the CAP method. The focal mechanism solutions of eight other small earthquakes with magnitude ML3.0～3.7 are inversed using the method of initial motion and amplitude ratio. Finally， using PSGRN/PSCMP software package， the static Coulomb stress triggering results of 10 earthquakes under different combinations of calculated parameters are calculated and calculated. The results show that 6 of the 9 earthquakes occurring in Xinhu reservoir area were triggered by the occurrence of earthquakes， occupying 66.7%； 1 earthquake was suppressed by the occurrence of earthquakes， taking up 11.1%； 2 earthquakes had nothing to do with the occurrence of earthquakes， accounting for 22.2% of the total. There is a weak static Coulomb stress triggering relationship among the 10 earthquakes with ML=3.0 in Xinhu reservoir area.

Keywords： Xinhu earthquake swarm； focal mechanism； double difference location； Coulomb stress； statistical test

引言

新湖水庫位于阳西县程村镇北部，是以防洪、灌溉为主，兼有养殖、供水等效益的中型水利工程。工程兴建于1958年10月，1963年12月竣工，总库容3768万m3，设计洪水水位21.22m，正常水位20.30m，死水位15m。库区属亚热带湿润季风气候，流域内降雨分配不均，4～9月丰水期雨量占全年总雨量的86%以上，枯水期为10～3月，占全年的14%，水库多年平均降雨量2350mm，年均径流量1575mm。

新湖水库所在的阳江地区构造活动强烈，1969年7月26日曾发生Ms6.4级地震[1]，但新湖水库背景地震活动却较弱，1970年以来，仅发生ML3.0级以上地震10次（图1），最大为2018年3月20日ML4.2级地震，该地震之前1个月内，库区曾密集发生2次ML3.0级以上地震。我们初步研究表明：2018年3月17日ML3.5级地震被2月18日ML3.2级地震触发，静态库仑应力增加0.0016Mpa，2018年3月20日ML4.2级地震被前两次地震联合触发，静态库仑应力增加0.0119Mpa，接近应力触发典型值0.01-0.1Mpa[2]，3次地震之间的存在应力触发关系，符合地震破裂断层已处于不稳定的边缘，地震可以看成围绕这一临界状态的涨落，小应力扰动也可以触发地震的假设[3]。但是，新湖水库全部ML3.0级以上的10次地震之间有无应力触发关系，尚没有文献进行研究，本文利用静态库仑应力，通过计算和统计大量参数组合的结果，对新湖库区ML3.0级以上地震间是否存在应力触发关系进行进一步研究。

越来越多的研究认为，当构造应力水平达到岩石的破裂临界极限时，应力触发现象普遍存在[4-16]，根据产生应力扰动源的不同，可进一步分为固体潮触发（日月运动引起）[8，9，10，11，12，13]，水库注水触发[14]，地震引起的静态库仑应力触发[4，6，7，15，16]，地震引起的动态库仑应力触发[5]等。Ide S et al.通过统计全球1000次大中型地震与固体潮的关系后发现，当潮汐应力增强时，全球大地震发生的比例也会增加，且新月和满月时，发生地震的震级更容易往大的方向扩展[10]。Lei研究发现紫坪铺水库蓄水在汶川地震破裂断层上造成的静态库仑应力增加大于0.1Mpa，二者间可能存在触发关系[14]。Deng研究南加州中强地震之间应力触发关系发现，95%的6级以上地震发生在库仑应力增加区域[15]。Stein研究土耳其北安那托利亚断层1939-1992年之间10个6.7级以上地震应力触发关系后，发现9个地震均被其他地震触发[16]。King等在研究landers地震与其余震的触发关系

后，认为0.1～0.6Mpa的静态库仑应力升高可触发地震[7]。关于应力触发阈值问题，Cochran等[9]研究发现全球5.5级以上的浅层逆冲型地震与固体潮有相关性，应力扰动幅值仅0.005Mpa。静态库仑应力触发的典型值为0.01Mpa，0.01Mpa的应力扰动能影响到余震的空间位置，且引起大地震的提前或延后几十年发生[17，18]。

为研究新湖库区10次ML3.0级以上地震之间是否存在静态库仑应力触发关系，本文首先对新湖库区有数字记录以来的10次ML3.0级地震进行了双差定位[19]，以减少地震间位置的不确定性对计算结果的影响[20，21]，然后，利用CAP方法[22]反演了2次ML4.0级左右地震的震源机制解和最佳拟合震源深度，利用Snoke[23]初动和振幅比方法反演了其他8个事件的震源机制解。最后，利用PSGRN＼PSCMP软件包[24]，考虑震源（源断层）和被触发的地震（接收断层）的震源深度、摩擦系数、走向、倾向、滑动角以及震源机制解的两个节平面等参数均有一定取值范围，利用不同计算参数组合，计算了各次地震被前面发生的地震的静态库仑应力触发关系。

图中阴影区域为水库水域，黑色沙滩球为采用CAP方法和Snoke方法反演的震源机制解，黑色五角星为精定位过后的震源位置，黑色圆点为台网地震目录中初始位置。

1 库仑应力变化

地震的发生可认为是地下断层面上所积累的剪切应力大于与正应力和摩擦系数有关的静摩擦力，造成闭锁的断层面开始活动，断层发生破裂进而引起地震，这一物理过程可以用岩石破裂的库仑应力准则来描述[6，7，15-18]，库仑应力用加载到地震破裂面的剪切应力和正应力来综合描述外来应力扰动对地震发震断层的影响，库仑应力变化为正，则促进断层发生破裂，库仑应力变化为负，这抑制断层发生破裂[6，7，15-18]。实际观测中，因加载到地震断层破裂面的构造应力的绝对值难以测量，但地震发生后，由于地震位错产生的位移和应变，可通过弹性位错理论计算，将这种应力应变扰动，投影到接收断层面上形成断层面上的正应力和剪切应力。按库仑应力计算公式组合：

2 数据和统计方法

2.1 地震精定位

因新湖水库地震监测台网较稀疏，为减少地震发生位置对静态库仑应力计算的影响，我们首先采用双差定位方法[19]对新湖水库10次ML3.0级以上地震进行了精定位。双差定位方法是一种相对定位方法，其基于如果两个地震事件之间的距离远小于其到地震台站之间的距离和介质速度变化的不均匀尺度，则这两个地震到台站的路径效应是相同的。因此，两地震之间同一震相的走时差便是两事件之间空间位置的差异。由于减小了速度模型对定位结果的影响，可以获得较高精度的地震间的相对位置，因此可减少静态库仑应力计算时地震之间相对位置的不确定性对计算结果的影响。我们利用广东省数字地震台网提供的10次地震的观测报告，仅使用震相到时差数据，利用2015年珠江口人工结构探测实验L1测线得到的一维速度模型[26]，对这10次事件进行了精定位。精定位结果如图1所示，可图1可见，除2007年11月26日ML3.7，精定位前后距离修正0.782km，2018年3月17日ML3.5，精定位前后距离修正0.656km，2008年2月4日ML3.2，精定位前后距离修正1.55km，其他7次地震精定位前后位置差别不大。

2.2 震源机制

本文首先采用CAP方法[22]反演了2007年12月4日 ML3.9和2018年3月20日ML4.2两次ML4.0级左右地震的震源机制解和最佳拟合震源深度。CAP方法将整个波形按照P波部分Pnl和面波部分Sur分别裁剪，将三分量记录中的5个窗口内的观测波形与理论计算波形进行互相关比较（图2），且可以在波形拟合时，给不同分量以不同的权重进行反演。充分考虑波形各个部分对震源机制的贡献，防止反演中面波权重过大的影响。之后，使用格点搜索方法，在适当的范围内循环搜索地震的深度、方位角、倾角、滑动角，得到相对误差最小时地震的震源机制和震源深度。

反演使用广东省数字地震台网记录到的震中距在300km范围内的宽频带波形数据，使用频率-波速域的FK方法[27]计算理论地震图，理论地震图计算所采用的一维速度模型，采用2015年珠江口海陆联测项目L1测线的人工深地震探测结果（图2a）[26]，计算了该一维速度模型下，不同震中距台站的理论格林函数，并按相同的窗长将实际和理论数据裁剪为Pnl和Sur两部分。P波窗长35s，滤波范围0.02～0.15Hz，Sur波窗长70s，滤波范围0.02～0.1Hz。然后利用格点搜索法即互相关方法，计算二者的互相关系数、振幅和到时差。图2给出了以2007年12月4日新湖ML3.9级地震为例的计算实例，经过筛选后的10个台站记录的38个震相，理论波形在不同方位角台站上都能较好得与实际记录波形匹配，且绝对振幅大小也相似，94.7%的互相关系数在0.6以上，各震相均符合较好，结果可信。

（a）理论地震图计算采用的一维速度模型；（b）CAP反演過程所使用的台站分布；（c）波形拟合误差随深度分布图，在深度空间搜索误差最小的CMT解，当震源深度为11.7km时，误差最小，为最优解；（d）最优11.7km深度理论计算波形和实际观测波形拟合效果如图示。其中，红线代表理论波形，黑线代表实际观测波形，波形右上角标注了理论波形和观测波形的相关系数，右下角标注了相应震相的时间移动，每行最左边标注了台站的台网标识和台站名。

对于其他8个震级较小的地震，本文采用Snoke方法[23]利用初至震相的初动符号和P波S波的振幅比计算了震源机制解。该方法利用双力偶点源模型，采用格点尝试法求取所有满足小于所设定的矛盾数上限的机制解。对求得的这些震源机制解，通过比较理论计算与实际观测得到的P波、SV波、SH波初动符号和振幅比，选择其中拥有矛盾符号数最少和振幅比残差最小的机制解，作为最佳震源机制解。我们选用信噪比较高的波形数据，手动标注信噪比较高的Pn、Pg初动极性和部分SH和SV的极性，并直接在速度记录上，在震相到时前3个周期内量取P波和S波的最大振幅。全部震源机制解的反演结果见表1和图1。从表1和图1可见，该区震源机制解以走滑和正断类型为主，节面取向与阳江地区NNW30°的主压应力场方向相符[28]。

2.3 统计方法

利用准备好的各次地震精定位后的震源位置、震源深度、震源机制解和地震震级等数据，利用PSGRN/PSCMP软件包[24]，采用均匀弹性半空间的位错理论[29]，计算了各次地震被之前发生的地震静态库仑应力触发结果，即源断层选择已发生的地震，接收断层选择即将发生的地震，断层破裂面分别选择震源机制解的两个节面，源断层的破裂面尺度和位错量根据Wells & Coppersmith 经验关系根据地震震级计算得到[29]。因小地震的震源机制解的精度没有大地震高，且无法根据余震空间分布或活动断裂展布等信息判断真实破裂节面，因此，我们分别对源地震断层破裂面的节面、节面的走向、倾向和滑动角、节面深度以及接收地震断层破裂面上的摩擦系数、节面、节面的走向、倾角和滑动角以及节面深度取一定的参数变化范围，对走向、倾角和滑动角分别变化±10°，步长2°，深度变化±3km，步长2km，等效摩擦系数取值0.4～0.8，步长0.2，并考虑不同的投影节面。分别计算了不同参数组合下，各次地震被之前发生的地震的静态库仑应力触发影響，最后，统计两次地震之间不同计算参数组合下的静态库仑应力变化的均值，图3以2018年3月20日ML4.2级地震为例，展示了不同计算参数组合下，该地震被前面已经发生的9次地震的静态库仑应力触发结果。从图3以及表2中给出的各次地震静态应力触发的均值来看，按地震发生的时间顺序，已发生的地震对该次ML4.2级地震静态应力变化的均值分别为113.5Pa、-45.5Pa、2256.8Pa、14.3Pa、0.1Pa、5.9Pa、20.1Pa、112.9Pa和428.4Pa，2018年3月20日ML4.2级地震被前面已发生的9次地震累计静态库仑变化总和为2906.5Pa，约0.0029Mpa，虽小于应力触发阈值0.01Mpa。但静态库仑变化为正，说明该次地震被前面9次地震所触发，地震发生时间被提前。

图3中，横坐标代表该次地震计算得到的静态库仑应力结果，纵坐标代表相应静态库仑应力变化数值所对应的百分比，各幅图的标题代表该次地震的发生日期和震级。

3 结果分析

表2详细给出了新湖库区10次ML3.0级以上地震之间的静态库仑应力计算结果，从各次地震间的静态库仑应力变化的均值来看，2018年3月20日ML4.2级地震被前面9次地震静态库仑应力触发值分别为113.5Pa、-45.5Pa、2256.8Pa、14.3Pa、0.1Pa、5.9Pa、20.1Pa、112.9Pa和428.4Pa，2018年3月20日ML4.2级地震被前面已发生的9次地震累计静态库仑变化总和为2906.5Pa，约0.0029Mpa，静态库仑变化为正，说明该次地震被前面9次地震所触发，地震发生时间被提前；2018年3月17日ML3.5级地震被前面8次地震静态库仑应力触发值分别为-4.2Pa、-2.4Pa、-201.7Pa、-1.3Pa、-0.5Pa、0.4Pa、246.8Pa和4573.1Pa，2018年3月17日ML3.5级地震被前面已发生的8次地震累计静态库仑变化总和为4610.2Pa，约0.0046Mpa，静态库仑变化为正，说明该次地震被前面8次地震所触发，地震发生时间被提前；2018年2月18日ML3.2级地震被前面7次地震静态库仑应力触发值分别为-12.4Pa、20.5Pa、-3686.5Pa、3.9Pa、0.5Pa、1.8Pa和63.7Pa，2018年2月18日ML3.2级地震被前面已发生的7次地震累计静态库仑变化总和为-3583.7Pa，约-0.0036Mpa，静态库仑变化为负，说明该次地震被前面7次地震所抑制，地震发生时间被延后；同样，2010年7月10日ML3.1级地震被前面已发生的6次地震累计静态库仑应力变化总和为-12.9Pa，静态库仑应力变化为负，说明该次地震被前面6次地震所抑制，地震发生时间被延后；2008年2月4日ML3.2级地震被前面已发生的5次地震累计静态库仑变化总和为7337.9Pa，约0.0073Mpa，静态库仑变化为正，说明该次地震被前面5次地震所触发，地震发生时间被提前；2008年2月1日ML3.1级地震被前面已发生的4次地震累计静态库仑变化总和为32046.0Pa，约0.032Mpa，大于应力触发阈值0.01Mpa，静态库仑变化为正，说明该次地震被前面4次地震所触发，地震发生时间被提前；2007年12月9日ML3.2地震被前面已发生的3次地震累计静态库仑变化总和为3434.6Pa，约0.0034Mpa，静态库仑变化为正，说明该次地震被前面3次地震所触发，地震发生时间被提前；2007年12月4日ML3.9地震被前面已发生的2次地震累计静态库仑变化总和为366.0Pa，静态库仑变化为正，说明该次地震被前面2次地震所触发，地震发生时间被提前；2007年11月27日ML3.3地震被2007年11月26日ML3.7级地震所触发，静态库仑应力增加1611.4Pa。从表2和以上分析可见，除了2008年2月1日ML3.1级地震被已发生的4次地震应力触发值0.032Mpa，大于应力触发阈值0.01Mpa外，其它地震之间的累计静态库仑应力触发值均小于应力触发的常见阈值0.01Mpa，6次地震被已发生的地震累计应力触发大小量级为0.001Mpa，与固体潮引起的应力扰动幅值相当[9]，另外还有2次地震被已发生的地震累计应力触发值较小，量级小于0.0001Mpa，远小于固体潮引起的应力扰动值，可视为无触发关系[32-35]。因此，9次地震中，6次地震被已发生地震所触发，占比66.7%，1次地震被已发生地震所抑制，占比11.1%，2次地震与已发生地震静态应力触发无关系，占比22.2%。

4 结论及讨论

本文针对阳江新湖水库有波形记录以来的10次ML3.0级地震，利用双差定位得到的精确的震源相对位置，CAP方法和Snoke方法得到的震源机制解数据，利用PSGRN/PSCMP软件计算并讨论了均匀弹性半空间介质模型下，新湖水库10次ML3.0级以上地震之间的静态库仑应力触发现象。通过统计不同计算参数组合下，地震之间的静态库仑变化结果，我们发现，2008年2月1日ML3.1级地震被已发生的4次地震应力触发值为0.032Mpa，大于应力触发阈值0.01Mpa，该次地震被已发生地震显著触发，其他6次地震之间的累计静态库仑应力触发值大小量级为0.001Mpa，与固体潮引起的应力扰动幅值相当[9]，另外2次地震累计应力触发值较小，量级小于0.0001Mpa，远小于固体潮引起的应力扰动值，可视为无触发关系。因此，9次地震中，6次地震被已发生地震所触发，占比66.7%，1次地震被已发生地震所抑制，占比11.1%，2次地震与已发生地震静态应力触发无关系，占比22.2%。根据小应力扰动也能触发地震的假设，新湖库区10次ML3.0级以上地震之间，存在微弱的静态库仑应力触发关系[36-38]。

致謝

本文事件波形和观测报告来至广东省数字地震台网，CAP方法采用Geotaos软件包[14]，SNOKE方法采用FOCMEC交互式程序包[31]，精定位采用HypoDD[19]，库仑应力采用PSGR

N/PSCMP软件包[24]，在此一并表示感谢。

参考文献：

[1]钟贻军，任镇寰.1969年阳江6.4级地震发震构造研究[J].大地测量与地球动力学，2003，23（4）：92-98.

[2]Scholz C H. The mechanics of earthquakes and faulting. Second edition[M]. The mechanics of earthquakes and faulting /. Cambridge University Press， 2002：496.

[3]吴忠良.地震震源物理中的临界现象[M].地震出版社，2000.

[4]万永革，吴忠良，周公威，等.地震静态应力触发模型的全球检验[J].地震学报，2002，24（3）：302-316.

[5]龚萱，陈棋福，彭志刚，等.2010年智利8.8级地震在北京房山岩体附近的动态触发活动[J].地球物理学报，2014，57（1）：115-128.

[6]Toda S， Stein R S， Richards-Dinger K， et al. Forecasting the evolution of seismicity in southern California： Animations built on earthquake stress transfer[J]. Journal of Geophysical Research Solid Earth， 2005，110（B5）.

[7]King G C P， Stein R S， Lin J. Static stress changes and the triggering of earthquakes[J]. Bull.seism.soc.am， 1994，78（3）：935-953.

[8]左玉玲.月球与太阳对地震的触发作用[J].国际地震动态，1990（12）：32-33.

[9]Cochran E S， Vidale J E， Tanaka S. Earth tides can trigger shallow thrust fault earthquakes.[J]. Science， 2004， 306（5699）：1164-6.

[10]Ide S， Yabe S， Tanaka Y. Earthquake potential revealed by tidal influence on earthquake size-frequency statistics[J]. Nature Geoscience， 2016，9（11）.

[11]解朝娣，吴小平，雷兴林，等.长周期潮汐与全球地震能量释放[J].地球物理学报，2013，56（10）：3425-3433.

[12]冯向东，魏东平.地震活动性与日月引潮力相关性统计分析[J].国际地震动态，2007，2007（5）：9-15.

[13]尹祥础，刘月.加卸载响应比——地震预测与力学的交叉[J].力学进展，2013，43（6）：555-580.

[14]Lei X. Possible roles of the Zipingpu Reservoir in triggering the 2008 Wenchuan earthquake[J]. Journal of Asian Earth Sciences， 2011，40（4）：844-854.

[15]Deng J， Sykes L R. Evolution of the stress field in southern California and triggering of moderate-size earthquakes： A 200-year perspective[J]. Journal of Geophysical Research Solid Earth， 1997，102（B5）：9859-9886.

[16]Stein R S， Barka A A， Dieterich J H. Progressive failure on the North Anatolian fault since 1939 by earthquake stress triggering[J]. Geophysical Journal International， 1997，128（3）：594-604.

[17]Reasenberg P A， Simpson R W. Response of regional seismicity to the static stress change produced by the loma prieta earthquake[J]. Science，1992，255（5052）：1687-1690.

[18]Rydelek P A， Sacks I S. Large earthquake occurrence affected by small stress changes[J]. Translated World Seismology， 2000，89（3）：822-828.

[19]Waldhauser F. HypoDD： A computer program to compute double-difference earthquake locations[J]. 2001.

[20]石耀霖，曹建玲.庫仑应力计算及应用过程中若干问题的讨

论——以汶川地震为例[J].地球物理学报，2010，53（1）：102-110.

[21]王力伟，陈棋福.以2008年汶川8.0级地震为例分析地震触发的静态库仑应力计算的不确定性[J].中国地震， 2010，26（3）：3-16.

[22]Zhao L S， Helmberger D V. Source estimation from broadband regional seismograms[J]. Bull.seismol.soc.am， 1994，84（1）：91-104.

[23]Snoke J A， Munsey J W， Teague A G， et al. A program for focal mechanism determination by combined use of polarity and SV-P amplitude ratio data[J]. Earthquake Notes， 1984，55（3）：15-40.

[24]Wang R J. A Simple Orthonormalization method for stable and efficient computation of green's functions[J]. Bull. Seism. Soc. Am.， 1999，89（3）：733-741.

[25]Parsons T， Ji C， Kirby E. Stress changes from the 2008 Wenchuan earthquake and increased hazard in the Sichuan basin[J]. Nature， 2009，454（7203）：509-510.

[26]Zhang X， Ye X， Lv J， et al. Crustal structure revealed by a deep seismic sounding profile of Baijing-Gaoming-Jinwan in the Pearl River Delta[J]. Journal of Ocean University of China， 2018，17（1）：186-194.

[27]Zhu L， Rivera L A. A note on the dynamic and static displacements from a point source in multilayered media[J]. Geophysical Journal of the Royal Astronomical Society， 2002，148（3）：619-627.

[28]康英，杨选，吕金水，等.广东及邻区地震的震源机制特征[J].中国地震，2005，21（3）：320-331.

[29]Okada Y. Internal deformation due to shear and tensile faults in a half space[J]. Bull. Seismol. Soc. Am.， 1992，82：1018-1040.

[30]Wells D L， Coppersmith K J. New empirical relationships among magnitude， rupture length， rupture width， rupture area， and surface displacement[J]. Bull.seism.soc.am， 1994，84（4）：974-1002.

[31]刘泽民，倪红玉，张炳，等.基于FOCMEC方法反演震源机制解的交互式程序研制与使用[J].华北地震科学，2015，33（1）：19-24.

[32]彭利媚，魏娅玲.由ML≥4.0地震震源机制解推断四川理县-黑水断裂构造特征[J].华南地震，2016，36（3）：9-16.

[33]郭增建，郭安宁，李健梅，等.基于汶川地震震例用震兆共迁方法对地震预测的讨论[J].华南地震，2016，36（4）：8-13.

[34]黄元敏，杜龙，邵叶，等.信宜-廉江断裂地震活动状态分析[J].华南地震，2016，36（4）：14-19.

[35]徐晓枫，王惠琳，胡久常，等.铺前-清澜断裂带附近地震的重定位及其构造意义的初步分析[J].华南地震，2017，37（2）：10-16.

[36]王小娜，邓志辉，叶东华，等.2016年10月广东乳源ML2.8级震群活动特征[J].华南地震，2017，37（3）：1-8.

[37]刘川琴，李发，刘东旺，等.金寨地震序列震源参数研究[J].华南地震，2017，37（3）：69-73.

[38]朱亮，缪鹏.基于应变能分析南北地震带中段地震活动对华东地区的影响[J].华南地震，2017，37（2）：53-57.