中国地震  2021, Vol. 37 Issue (2): 400-414
山东乳山震群遗漏地震检测及震源区活动特征分析
王鹏1,2, 郑建常1     
1. 山东省地震局, 济南 250102;
2. 东营市地震监测中心, 山东东营 257000
摘要:对震群活动特征的深入研究可以为区域地震危险性分析和地震活动趋势判断提供有效的依据,但受台网布局和震群位置的影响,地震目录中往往会遗漏一些地震,而地震目录的完整性将会影响震群活动特征分析的可靠性。因此,本文利用基于GPU加速的模板匹配方法对山东乳山震群2014年5月至2015年6月期间固定台站和流动台阵记录的连续波形进行扫描,检测出7500个地震事件,为台网目录的4倍,使完备震级从1.0级降至0.2级。基于完备地震目录,利用传染型余震序列(ETAS)模型分析了乳山震群的活动特征。ETAS模型中的随机成分可以作为外力作用触发地震的指标。研究结果表明,乳山震群外力触发的地震所占比例较大,结合前人研究结果,初步推测流体作用可能是外力因素之一,其在乳山震群的发展过程中起一定的作用,但不同时期外力作用对微震活动的触发强度不同。
关键词ETAS模型    乳山震群    模板匹配    完备震级    外力作用    
Analysis of Seismic Characteristics and Missing Earthquake Detection of the Rushan Earthquake Swarm in Shandong Province
Wang Peng1,2, Zheng Jianchang1     
1. Shandong Earthquake Agency, Jinan 250102, China;
2. Dongying Earthquake Monitoring Center, Dongying 257000, Shandong, China
Abstract: The study of swarm activity characteristics can provide effective basis for regional seismic risk analysis and seismic activity trend judgment. However,due to the influence of network layout and swarm location,some earthquakes are often missed,and the integrity of earthquake catalogue will affect the reliability of swarm activity characteristics analysis. Earthquake swarms have occurred frequently in North China since 2013, especially in Rushan,Shandong Province. More than 13000 small earthquakes were recorded in Rushan area on October 1,2013, and the maximum was ML5.0 on May 22,2015. During this period,the seismicity increased and decreased for many times. In this paper,the GPU accelerated template matching method was used to scan the continuous waveforms recorded by the fixed stations and mobile arrays during the period from May 2014 to June 2015 of the Rushan earthquake swarm in Shandong Province. A number of 7500 earthquake events were detected,which was 4 times of the catalogue of the network,and the complete magnitude reduced from 1.0 to 0.2. Based on the completed earthquake catalogue,the activity characteristics of Rushan earthquake group were analyzed by using the model of infectious aftershock series(ETAS). The random components in ETAS model is used as the index of external force triggering earthquake. The Rushan earthquake swarm began to become intensively active in October 2013, while the Rushan array began to operate on May 5,2014. As a result,we added the detection catalogue to the earthquake catalogue of the whole Rushan earthquake swarm,divided the whole development process of the earthquake swarm into five active periods,and analyzed the seismicity of the Rushan earthquake swarm and its variation characteristics in each stage with the help of ETAS model. The results showed that the proportion of the earthquakes triggered by the external force of the Rushan earthquake swarm was large. Combined with the previous research results,it was preliminarily speculated that the fluid effect may be one of the factors of the external force and play a certain role in the development process of the Rushan earthquake swarm,but the triggering intensity of the external force on the microseismic activity was different in different periods.
Key words: ETAS model     Rushan earthquake swarm     Template matching     Magnitude of completeness     External force    
0 引言

2013年以来华北地区震群活动频发,其中山东乳山震群尤为突出。2013年10月1日乳山地区地震丛集活动持续至2019年,共记录到小震1.3万余次,其中ML≥3.0地震31次,最大震级为2015年5月22日ML5.0,期间伴随着多次地震活动增强—减弱的过程。

乳山震群所在的胶东半岛地区,位于华北地块东缘与苏鲁造山带山东段的交汇位置,在大地构造划分上属于华北克拉通与扬子克拉通的结合部位(晁洪太等,2001周瑶琪等,2015)。受苏鲁造山带与古太平洋板块俯冲的双重影响,该区域岩浆活动十分频繁,自太古代到新生代均有发现岩浆活动,其中以燕山晚期岩浆活动最为剧烈,同期还发生了大规模的岩浆侵入作用。在火山活动区域,高热流区会导致较小的低粘度累积弹性应变,相对于主余型序列而言,往往会产生更多的震群(Ben-Zion et al,2006Yang et al,2009)。根据历史地震记录,胶东半岛地区历史上大震较少,却经常发生震群活动(李冬梅等,2011),但这些小震序列(震群)持续时间通常很短,一般持续几天,最多不会超过1个月。而乳山震群持续6年,且频度和强度均较大,这在大陆东部地区甚为罕见。对乳山震群震源区活动特征的深入研究可为区域地震危险性分析和地震趋势判断提供有效的依据。

火山或地热地区的震群发震机制通常可以用流体运移模型来解释,在地壳中运移的流体包括水和挥发性气体(CO2)等,其通常被认为是由中地壳到上地幔这一深度的岩浆源提供的(Takahata et al,2003Ohba et al,2011)。在前人对乳山震群发震机制的研究中,也猜测流体可能是震群的触发因素(Zheng et al,2017郑建常等,2019),但始终缺乏直接的观测证据,而在震群发震构造(曲均浩等,201420152019郑建常等,2015a2015b李铂等,2017张斌等,2017)、震源参数(苗庆杰等,2016王鹏等,2016郑建常等,2016刘方斌等,2018)等问题的讨论上也存在不同见解。究其原因为,乳山地区断裂错综复杂,次级断裂(近SN向)、隐伏断裂(NE、NW向)较为发育,断裂格局总体呈现羽状或雁列式形态(刘善宝,2005),震中区并无活动断层,距离最近的乳山断裂15km。

统计学方法可以从复杂的现象中探究本质规律,自Ogata(1988)在“大森-宇津”公式(Omori,1894Utsu,1961)的基础上提出了“传染型余震序列”(Epidemic-type Aftershock Sequence,ETAS)模型后,该模型已受到广泛的关注,其在研究地震活动变化规律(Ogata,19922001Zhuang,2000蒋海昆等,2007蒋长胜等,2013a2013b2013c20152017)、探测应力变化(雷兴林等,2008Lei et al,2017)、识别与流体作用有关的诱发地震(Hainzl et al,2005Lei et al,2008龙锋等,2010蒋海昆等,20112012)等方面得到了长足发展。

统计地震学方法很大程度上依赖于地震目录,地震目录的不完备势必会影响到基于地震目录的相关统计学结果的可靠性(李智超等,2014)。地震目录的完备性主要依赖于研究区域地震台网的监测能力,而乳山震群周围台站分布格局不够理想,地震序列短时间的丛集发生通常会造成地震的遗漏。通过遗漏地震检测来使地震目录更加完整(Schaff,2008),可以提高分析结论的可靠性。用于微震检测的模板匹配滤波技术(Gibbons et al,2006)提出后,该技术发展迅速并得到了广泛的应用,在检测地脉动(Shelly et al,2007)、研究余震序列变化(Peng et al,2009侯金欣等,2017王鹏等,2016)、确定发震构造(Yang et al,2009谭毅培等,2016)、分析震源区应力变化(Meng et al,2013)、识别重复地震(马腾飞,2015)、监视核爆(Zhang et al,2015)等诸多方面发挥着重要的作用。

本文利用基于GPU并行计算的模板匹配方法检测乳山震群的遗漏地震,基于完备地震目录对该震群ETAS模型的统计参数进行分析,并初步推测乳山震群的发震机制。

1 理论方法

本文所用模板匹配方法的处理流程参照Peng等(2009),并与侯金欣等(2017)的处理步骤类似,具体方法和处理流程可见Wang等(2019),在此不再赘述,仅对ETAS模型做简单介绍。

ETAS模型条件强度函数λ(t)(Ogata,1988蒋海昆等,2012)为

$ \lambda (t) = \mu (t) + K\sum\limits_{}^{{t_i} < t} {\frac{{{{\rm{e}}^{\alpha \left({{M_i} - {M_z}} \right)}}}}{{{{\left({t - {t_i} + c} \right)}^p}}}} $ (1)

式中,右侧第二项为第i次地震对地震发生率的贡献,i遍历地震序列中的所有地震,Mz为参考震级,其必须大于等于最小完整性震级Mcp代表地震序列衰减的快慢,p大则衰减快,反之则衰减慢;α表示触发次级余震的能力,α越大表示越强的高阶余震的激发能力;μ(t)最初的物理含义为剔除地震丛集后单位时间内的理论地震数目,也表示由于外力作用触发微震活动的强弱,与地震间的自激发过程无关(Hainzl等,2005);c为数值极小的正值,其保证了右侧第二项的分母不为0,与主震后短时间内数据的完备程度(蒋海昆等,2007)和震源区的应力状态有关(Narteau et al,2009);K为与主震震级及bpαcμ参数有关的地震活动水平,其中b为G-R关系中的比例系数(宋金等,2009)。

对于同一组观测数据,数学上可采用多个统计模型对其进行拟合,并可得到多组符合条件的参数。Ogata(1988)利用最大似然法估计ETAS模型参数,并利用最小化模型的信息量准则(Akaike information criterion,AIC)进行最佳模型参数的判别(Akaike,1974)。AIC准则既考虑模型对观测数据拟合的优劣,也会惩罚为了提高拟合程度而无限制地增加模型参数的行为。ETAS模型的对数似然值为

$ {\lg L = \sum\limits_{i = 1}^{{N_{{\rm{AIC}}}}} {\lg } \lambda (t) - \int\limits_{{t_b}}^{{t_e}} \lambda (t){\rm{d}}t} $ (2)
$ {{\rm{AIC}} = - \lg L + 2k} $ (3)

其中,tbte为用于AIC计算的起止时间;NAIC为计算时段内的地震数目,为了保证计算结果的可靠性,计算时段内的地震目录必须保证基本完备;k为待拟合参数个数,本文中k=5,取AIC最小时对应的模型参数作为最佳模型。

由式(1)可知,ETAS模型所表征的地震活动率包括两部分,其中右侧第二项是由于“余震激发余震”所导致的地震丛集,右侧第一项是与余震激发无关、因外力作用而产生的地震活动。因而,通过对ETAS模型两部分地震活动进行分离,将外力触发地震比例Rb用式(4)表示,而自激发地震所占比例Rt=1-Rb,可以根据各自所占的比重来讨论外因触发与地震活动的关系(蒋海昆等,2011)

$ {R_b} = \frac{{\int_{{t_b}}^{{t_e}} \mu (t){\rm{d}}t}}{{\int_{{t_b}}^{{t_e}} \lambda (t){\rm{d}}t}} $ (4)
2 数据处理

乳山震群周围100km内共有9个固定地震台站,均位于震群的北部、西部和东部,而南侧无固定台站。乳山流动台阵于2014年5月架设完成,共计18套宽频带地震仪,采样率为100Hz(图 1(b)),平均台间距约为0.5km,距离震群最大距离约10km。乳山台阵对震群的包络性较好,在震群南侧架设的宫家岛台(GJDT)弥补了震群南侧海岛无固定台站的缺陷。乳山台阵自2014年5月投入使用,乳山震群持续时间较长,震群最大地震发生在2015年5月,地震活动主要集中在2015年6月,后续地震活动开始逐渐衰减。因此,选择2014年5月至2015年6月期间固定地震台站和流动台阵的连续波形进行遗漏地震检测。

图 1 乳山震群和地震台站空间分布 (a)中黑色三角形代表固定台站,蓝色三角形代表乳山台阵,灰色圆圈为台网定位结果,红色圆圈为重定位结果,红色方框为研究区域,红色五角星为乳山震群;(b)为图(a)中黑色方框区域放大

首先,对连续波形数据去均值、去线性趋势,并采用零相移4阶Butterworth带通滤波器对2~10Hz进行带通滤波。选取台网目录(中国地震台网中心)中ML≥2.0的地震作为模板地震,根据台网目录中发震时刻截取乳山台阵中距离乳山震群最近且分布较好的RSH、GJDT、DFS、NHNC、WHUT和XJZT等6个台站的三分量记录作为模板地震波形,再依据震相报告读取模板地震的P和S波到时,计算并选取信噪比大于5的109个地震作为模板地震,最大震级为ML4.5。然后,按照模板匹配方法的数据处理流程(Wang et al,2019),利用模板地震与连续波形互相关计算后做叠加得到互相关系数波形,为了减少误检测事件,本文选取12倍的绝对离差中位数(Median Absolute Deviation,MAD)作为判定新检测地震的阈值,并设定5s时间范围内仅检测1次地震事件,选取互相关系数最大的遗漏事件,最后确定检测事件的震级和发震时刻。

图 2为检测到的一次ML2.7遗漏地震的示例。其中,图 2(a)为互相关系数和MAD的倍数,可见在1s的时间内互相关系数出现多次高值,且均高于12倍MAD的阈值,但根据假定只取互相关系数最大的点,即图中红色圆圈所对应的时间;图 2(b)表示了遗漏地震所在的波形序列的位置,可见在遗漏地震之前连续发生了2次较大地震,分别是2014年9月16日14时42分59秒的ML4.1和14时43分33秒的ML3.9地震,而该遗漏地震几乎淹没在前面较大地震的尾波里,幅值是前面较大地震的1/10,不易分辨,这也是地震漏检测的原因;图 2(c)2(d)为RSH台和DFS台记录的模板地震在连续波形上的扫描结果,可以看到模板地震与检测出的遗漏事件相关性较好,相关系数为0.61。同时,模板地震的S波的幅值与遗漏事件相当,因此根据S波振幅比计算的遗漏地震的震级也为ML2.7。依据检测事件与所对应的模板事件具有相同走时的假设,确定检测地震事件的发震时刻为2014年9月16日14时44分03秒。

图 2 模板匹配方法检测遗漏地震示例(以ML2.7地震为例) (a)为经过偏移、叠加和取平均值后的互相关系数和MAD的倍数图,其中相关系数大小用竖线表示,MAD的倍数用橙色点表示;(b)中红色点线框为遗漏地震的位置,即图(a)所示的时间范围;(c)、(d)为DFS台和RSH台记录的模板地震和遗漏地震事件的波形对比,其中红色为2014年11月19日3时24分30秒的ML2.7模板地震,蓝色为该时段内的连续波形,检测地震事件为2014年9月16日14时44分3秒的ML2.7地震
3 结果分析 3.1 遗漏地震检测结果

2014年5月5日至2015年6月30日期间,共检测到7500个地震事件,其中109个为模板事件的自检测,7391个为检测地震事件,检测地震事件也包括了台网已记录但未用作模板的地震事件。山东地震台网记录到地震1937个,检测地震事件是台网目录的近4倍,由模板匹配方法获得的地震目录简称为检测目录。根据检测目录和台网目录的频度对比(图 3)可见,与台网目录相比,检测目录中的遗漏地震主要集中在ML≤1.0的震级范围,1.0<ML≤2.5的震级范围内遗漏地震事件较少,ML>2.5时几乎没有遗漏的地震,这与震群所在区域台网的监测能力有一定关系。

图 3 检测目录和台网目录频度和累积频度对比

利用最大曲率法对比分析了台网目录和检测目录的完备震级,该方法选取震级-频率曲线(Frequency Magnitude Curve)中斜率最大值所对应的震级作为Mc(Wiemer et al,2000)。在实际应用中,Mc震级通常对应非累积频度-震级分布(Non Cumulative Frequency Magnitude Distribution)中具有最大地震频度的震级(黄亦磊等,2016)。

图 3可以看出,在非累积频度-震级分布中,检测目录和台网目录的最大值对应的完备震级分别为0.1级和0.9级,由于最大曲率方法得到的完备震级通常偏小,利用最大似然法进行拟合,发现检测目录和台网目录分别在震级0.2和1.0级以上b值的拟合较为稳定,最小完整性震级分别取0.2级和1.0级,计算的b值分别为0.95±0.01和0.89±0.03。基于模板匹配方法的遗漏地震检测将研究时段内乳山震群的完备震级从1.0级降至0.2级,弥补了台网目录中1.0级以下地震的缺失,使地震目录更加完整。

3.2 利用ETAS模型分析地震活动特征

乳山震群自2013年10月1日开始,期间经历了多个增强—减弱的阶段,根据日频度和累积频度随时间的演化特征(图 4),将整个震群发展过程分为5个活动时段:第一阶段为2013年10月1日至2014年1月6日,2013年10月1日12时7分55秒震群初始地震ML3.8地震发生后,地震活动短时间内集中发生,日频次达到高值,并在接下来的2个月内逐渐减弱至2014年1月6日;第二阶段开始于2014年1月7日,乳山震群发生了ML4.7地震,随后地震活动再一次增强,日活动频度逐渐减弱后又出现了短时间的增强,至2014年4月3日又趋于减弱;第三阶段则开始于2014年4月4日的ML4.6地震,与前一阶段类似,地震活动在短时间的增强后迅速减弱;第四阶段为震群最大地震2015年5月22日ML4.6地震发生后,日频度达到最高值,随后地震活动经过半年多时间的调整,期间也伴随着增强—减弱的过程,至2015年底基本趋于稳定;第五阶段为2016年以后,地震活动较弱,总体较为稳定,期间仅有个别地震活动小幅度增强的现象。

图 4 乳山震群日频度和累积频度随时间的演化 左侧纵坐标表示日频度,用黑色点线表示;右侧纵坐标表示累积频度,用红色线表示;蓝色点线将整个发展过程分为5个阶段,黑色点线又将每个阶段分为若干个子阶段

乳山台阵于2014年5月5日开始运行,故遗漏地震检测只针对2014年5月至2015年6月期间,相对于整个震群的发展过程,该时间位于第三活动阶段和第四阶段(图 5 (a)),从该时段内检测目录和台网目录随时间的变化对比(图 5 (b))来看,尽管检测目录和台网目录在低震级段的地震数量有所差别,但其反映的地震活动时段一致,例如在2015年5月22日前后地震活动有明显的差异。同时,在模板匹配方法检测时,台网记录的连续波形中有几天数据缺失,造成了检测目录中有几天地震缺失,但其并不是地震丛集时段,不会对地震活动性的分析造成影响。

图 5 乳山震群震级、频度-时间图 (a)为检测目录(黑点)和台网目录(灰点)的震级随时间变化图,其中红色点线框为本文利用模板匹配方法检测的时间段;(b)为图(a)红色框的放大,并加入累积频度的统计,其中蓝点为检测目录,黑点为台网目录,红色线为累积频度曲线

根据震群发展过程,每个阶段又分为若干个小的阶段(图 4)。将检测得到的完备震级加入到台网目录中,利用ETAS模型对这些不同的阶段分别进行反演,得到各模型参数的变化特征,从而得出外力作用和自激发作用对震群发展的影响。

不同的截止震级对ETAS模型参数的反演影响较大(蒋海昆等,2012蒋长胜等,2013c),要保证每个阶段的ETAS模型参数反演中截止震级必须大于完备震级。因此首先利用震级-序号法考察完备震级在震群发展过程中的变化(图 6),震级-序号法是按地震发生的次序排列,通过分析不同震级的地震概率密度分布来定性分析完备震级Mc的方法,其中地震概率密度最大的位置所对应的震级即为完备震级Mc。之所以使用地震发生次序,而不是地震时间,主要是为了避免余震或震群中“丛集事件”的影响(蒋长胜等,2013a)。这里以500个地震作为窗长,以0.1级为间隔,计算窗内各个震级地震的概率密度,可以看到在震群发生初期完备震级在0.9级左右,在检测目录加入后使得完备震级明显降低,在检测的2014年5月至2015年6月期间,完备震级基本维持在0.1级左右,在2015年6月之后完备震级又上升到0.9级左右。为了分析乳山震群整个发展过程中ETAS模型参数的趋势变化,并能更好地与其他阶段的参数变化进行对比,因此统一将截止震级取为0.9级,按照图 4中的阶段划分,反演了各个阶段的ETAS模型参数,见表 1,其中第一阶段(Ⅰ)至第四阶段(Ⅳ)以及每个阶段的子阶段(Ⅲ1、Ⅲ2,Ⅳ1、Ⅳ2)的拟合图见图 7

图 6 乳山震群震级-序号图 色标表示计算窗内地震在某个震级段的概率密度分布

表 1 乳山震群不同阶段的ETAS模型参数

图 7 乳山震群各阶段ETAS模型拟合结果

对于第一阶段,背景地震或外力作用的比例为35.5%,说明有35.5%的地震是由外力作用触发的,剩余64.5%的地震由余震的自激发产生。α值为0.947,低于典型的构造地震(2.0),震群型地震的α值通常小于1(Ogata,2001)。日本地区震群的α值一般介于0.35~0.85之间,而一般的主余序列则在1.2~3.1之间(Ogata,1992)。α值表征余震的激发能力,这一阶段的α值接近1,表明有较强的余震的激发能力。p值为余震衰减率,这里p值(0.99)与典型构造地震(1.0左右)相当,表明序列属于正常衰减。同时,在这一阶段模型理论值与观测值拟合较好,表明参数估计较为准确,外力作用在震群发生初期占一定比重,但主要还是受余震的自激发作用影响。

第二阶段的地震活动伴随着较大地震(2014年1月7日ML4.7地震)短时间内集中发生,随后逐渐衰减,期间发生过1次日频度突然增强的现象。这一阶段外力作用所占比例有所下降,从第一阶段的35.5%降至24.3%,表明外力作用较弱,p值无变化,仍为0.99,表明序列仍处于正常衰减,但α值较高(1.258),且高于第一阶段,表明余震的激发能力增强。

第三阶段与第二阶段类似,伴随着较大地震(2014年4月4日ML4.6地震)的发生,地震活动出现丛集,随后趋于衰减,由于地震活动速率变化较大,将该阶段分为2个子阶段(Ⅲ1、Ⅲ2图 4),其中Ⅲ1阶段外力作用所占比例有所上升(39.9%),表明外力作用有增强的趋势,这种短期丛集的地震受外力作用触发的较多。α值明显下降,由1.258降至0.362,表明余震的激发能力明显减弱,与外力作用增强一致。同时p值也明显增大,变为1.5,表明序列的衰减作用增强。在这一阶段,模型理论值与观测值拟合较差,实际地震发生率(图 7(c)蓝色实线)明显高于拟合的理论曲线(图 7(c)灰色点线),表明大森定律型自触发只有很小的促进作用,外力作用引起的地震活动的比例可能被低估(Lei et al,2008)。在Ⅲ2阶段,外力作用所占比例继续增大(45.9),说明将近一半的地震是由于外力作用触发的,同时α值明显升高(1.444),表明这一活动时期地震活动的自激发作用也相应增强,p值也再次回到1.01,序列呈现天然地震的正常衰减速率。

第四阶段地震活动速率明显增加,也伴随着震群最大地震(2015年5月22日ML5.0地震)的发生,日频次达到最高值,随后快速下降。将这一活动时段也分为2个子阶段(Ⅳ1、Ⅳ2图 4),其中Ⅳ1阶段外力作用较上一阶段有所降低(35.2%),约有1/3的地震是由外力因素触发的,而α值依然保持较高值(1.002),表明这一阶段更多的是受最大地震(ML5.0)的影响,激发了短期内的余震发生,p值为1.06,表明仍在正常衰减。Ⅳ2阶段在衰减(p值)变化不大的情况下,α值略微下降(0.969),说明受最大地震(ML5.0)的影响减弱,同时受外力作用明显增强(43.6%),与Ⅲ1阶段类似,该阶段的实际地震发生率(图 7(f)蓝色实线)也明显高于拟合的理论曲线(图 7(f)灰色点线),表明外力作用也可能被低估(Lei et al,2008)。外力作用在这一阶段可能起主要作用,而地震的自激发作用次之。

第五阶段随着地震活动的减弱,外力作用比例减弱,α值和p值也有所下降,说明这一阶段外力作用和自激发作用均在减弱,同时衰减也相应减弱,与这一时期持续较长时间的弱地震活动水平有关。

我们还考察了乳山震群的ETAS模型参数Rbαp随时间的演化规律(图 8),可以看到除个别点外,Rbα的变化趋势较为一致,总体呈现先上升后下降的趋势,而p值除在个别点外,其他时段变化较小,基本维持在1.0左右,表明乳山震群一直保持着天然地震正常的平均衰减速率。Rb值在整个过程中基本维持在20%~50%之间,表明外力作用在震群的发展过程起到一定的作用,其在某些阶段甚至接近50%,但总体来看仍是震群前期地震的自激发起主导作用。

图 8 ETAS模型参数随时间的变化
4 讨论 4.1 不同截止震级对ETAS模型反演结果的影响

通过对2014年5月至2015年6月期间的遗漏地震进行检测,使得完备震级大幅减小,在完备震级随时间的演化中也能明显看到这种特征。2014年5月至2015年6月这一时间段对应Ⅲ2、Ⅳ1阶段(图 4),因此,测试这2个阶段在截止震级分别取0.1级和0.9级时ETAS模型的拟合结果,比较不同参数的结果(图 9)。从图 9可以看到,外力触发地震的比例Rb差别较大,截止震级取0.1级时2个阶段的Rb值明显低于截止震级取0.9级的情况。由于任何一个地震均能以一定概率触发后续地震,低于Mc的大量微震事件将可能显著改变背景地震概率(Rb值)(蒋长胜等,2013b)。在III2阶段,随着截止震级的减小,α值和μ值明显减小,表明序列触发微震的能力和外力作用触发微震活动的能力均有所减弱。在Ⅳ1阶段,随着截止震级的减小,α值和μ值则有所增加,表明自激发作用和外力作用均有所增强。这种现象可能与2个阶段的地震活动特征不同有关,III2阶段地震活动趋势较为平稳,而Ⅳ1阶段由于震群最大地震的发生,地震活动频度迅速增加,日频度达到最大值(图 5),震级越大,激发高阶余震的能力越强,对后续地震的影响也就越大。尽管α值和μ值均有所增加,受震群最大地震的影响,外力触发地震的比例却很低(3.7%),地震活动中仍是自激发作用的贡献更大。另一方面由于震群中最大地震发生后,短期内ETAS模型参数随时间变化较为明显且波动较大,ETAS模型参数与地震序列持续时间也有一定关系(蒋长胜等,2013c),2个阶段的持续时间相差较大,可能也是一方面因素。此外,震源区的破裂特征、断层的愈合速度、应力加载、周边构造场的调整等诸多因素均会使地震序列的衰减特征和激发余震能力出现明显差异。

图 9 不同截止震级下ETAS模型拟合结果对比 (a)、(b)截止震级为0.1级时2个阶段的拟合结果;(c)、(d)截止震级为0.9级时2个阶段的拟合结果

由于ETAS模型假定每一次地震均可触发其后发生的任何震级的地震,因此,截止震级在确保地震目录完整性的同时,也可能切断了截止震级以下地震与以上地震的触发关系,从而造成“链接缺失”现象(Wang et al,2010)。通过遗漏地震检测,使得完备震级减小,而截止震级的减小尽可能地减少了这种现象的出现,因此地震目录的完整性对ETAS模型参数反演至关重要。

经过遗漏地震检测,在保证地震目录完整性的基础上,降低了截止震级,计算得到的ETAS模型参数也更加可靠。从理论曲线和实际地震活动的拟合来看,截止震级取0.1级时拟合较好,同时AIC值也更低(图 9)。但由于乳山台阵于2014年5月才架设运行,而乳山震群2013年10月开始集中活动,因此基于台阵波形数据的遗漏地震检测的时间段仅是震群发展过程中的一部分(图 5)。为了考察乳山震群整个发展过程中ETAS模型参数的趋势变化(图 8),同时能更好地与其他阶段进行对比(图 7),对于整个震群发展过程取截止震级0.9级以上的完备地震目录进行反演计算。将截止震级设为高于完备震级(0.9级)的某个震级,利用ETAS模型进行反演,尽管各参数的数值不完全相同,但反映出各参数随时间的变化规律是一致的。

4.2 发震机制讨论

本文利用外力作用引起的地震与自激发地震的比例系数Rb来考察外力因素对震群活动的影响。除个别阶段外,Rbα的总体变化趋势较为一致,呈现出先升高后降低的过程,而p值除个别阶段外,总体变化较小,维持在天然地震正常的衰减速率。Rb值在整个过程中基本维持在20%~50%之间,表明外力作用在震群的发生、发展过程中均占不小的比例,尤其在震群初期或较大地震发生之前,地震可能是受外力因素触发或间接引起的应力调整引发的,但从整个过程的比例来看,余震的自激发作用所占比重更大,发挥着更大的作用。震群发展后期(V阶段)外力作用和自激发作用均在减弱,地震活动呈现较长时间的弱活动水平。

前人在对乳山震群空间分布、震中迁移速率、ETAS模型、高频衰减系数等方面的研究中推测乳山震群可能是由流体入侵触发的(郑建常等,2015b2016Zheng et al,2017),但目前仍没有确切的观测证据。通过分析ETAS模型参数特征和其随时间演化的规律,也可以看到外力因素在乳山震群的发生、发展过程中起着一定的作用,而流体入侵或运移也可能是外力因素之一,今后还需结合其他依据做详细的论证。

5 结论

本文利用基于GPU并行运算的模板匹配滤波技术,对于乳山震群2014年5月5日至2015年6月30日期间的连续地震波形进行检测,识别出7500个地震事件,其中109个为模板事件的自检测,7391个为检测地震事件。检测事件是台网目录的4倍,利用较为准确的b值拟合方法得到的震群完备震级从1.0级降至0.2级。

乳山震群2013年10月开始集中活动,而乳山台阵自2014年5月5日才开始运行,因此,将检测目录补充到整个乳山震群的地震目录中,并划分整个震群发展过程为5个活动时段,借助ETAS模型分析了乳山震群的地震活动性以及在各个阶段的变化特征。结果表明,乳山震群外力触发的地震所占比例较大,结合前人研究,初步推测流体作用可能是外力因素之一,且不同时期外力对微震活动的触发强度不同。

参考文献
晁洪太、于慎谔、李家灵等, 2001, 山东半岛地区活断层研究, 东北地震研究, (4): 1-8. DOI:10.3969/j.issn.1674-8565.2001.04.001
侯金欣、王宝善, 2017, 2014年鲁甸MS6.5地震前后地震活动性, 地球物理学报, 60(4): 1446-1456.
黄亦磊、周仕勇、庄建仓, 2016, 基于地震目录估计完备震级方法的数值实验, 地球物理学报, 59(4): 1350-1358.
蒋长胜、吴忠良、韩立波等, 2013a, 地震序列早期参数估计和余震概率预测中截止震级Mc的影响: 以2013年甘肃岷县-漳县6.6级地震为例, 地球物理学报, 56(12): 4048-4057.
蒋长胜、吴忠良、尹凤玲等, 2015, 余震的序列参数稳定性和余震短期发生率预测效能的连续评估——以2014年云南鲁甸6.5级地震为例, 地球物理学报, 58(11): 4163-4173.
蒋长胜、吴忠良、庄建仓, 2013b, 地震的"序列归属"问题与ETAS模型——以唐山序列为例, 地球物理学报, 56(9): 2971-2981.
蒋长胜、庄建仓、龙锋等, 2013c, 2013年芦山MS7.0地震序列参数的早期特征: 传染型余震序列模型计算结果, 地震学报, 35(5): 661-669.
蒋长胜、庄建仓、吴忠良等, 2017, 两种短期概率预测模型在2017年九寨沟7.0级地震中的应用和比较研究, 地球物理学报, 60(10): 4132-4144. DOI:10.6038/cjg20171038
蒋海昆、宋金、吴琼等, 2012, 基于ETAS模型对三峡库区流体触发微震活动的定量检测, 地球物理学报, 55(7): 2341-2352.
蒋海昆、杨马陵、孙学军等, 2011, 暴雨触发局部地震活动的一个典型例子: 2010年6月广西凌云-凤山交界3级震群活动, 地球物理学报, 54(10): 2606-2619. DOI:10.3969/j.issn.0001-5733.2011.10.018
蒋海昆、郑建常、吴琼等, 2007, 传染型余震序列模型震后早期参数特征及其地震学意义, 地球物理学报, 50(6): 1778-1786. DOI:10.3321/j.issn:0001-5733.2007.06.018
雷兴林、马胜利、闻学泽等, 2008, 地表水体对断层应力与地震时空分布影响的综合分析——以紫坪铺水库为例, 地震地质, 30(4): 1046-1064. DOI:10.3969/j.issn.0253-4967.2008.04.021
李铂、崔鑫、娄世平等, 2017, 山东乳山地区震群特征及发震背景再研究, 震灾防御技术, 12(3): 491-500.
李冬梅、周翠英、董翔等, 2011, 山东地区震群活动与周围地区中强地震的关系, 华北地震科学, 29(4): 21-26. DOI:10.3969/j.issn.1003-1375.2011.04.004
李智超、黄清华, 2014, 基于概率完备震级评估首都圈地震台网检测能力, 地球物理学报, 57(8): 2584-2593.
刘方斌、曲均浩、李亚军等, 2018, 山东乳山地震序列震源机制解一致性参数特征, 地震地质, 40(5): 1086-1099.
刘善宝, 2005. 山东乳山金青顶金矿田成矿规律及其成矿远景研究. 硕士学位论文. 西安: 长安大学.
龙锋、杜方、阮祥等, 2010, 用ETAS模型分析自贡矿井注水触发地震, 中国地震, 26(2): 164-171. DOI:10.3969/j.issn.1001-4683.2010.02.004
马腾飞, 2015. 连续波形中"重复地震"的识别——与WFSD相联系的汶川地震序列分析. 博士学位论文. 北京: 中国地震局地球物理研究所.
苗庆杰、刘希强, 2016, 山东乳山震群剪切波分裂参数的时间演化特征分析, 地震学报, 38(2): 220-231.
曲均浩、蒋海昆、李金等, 2015, 2013—2014年山东乳山地震序列发震构造初探, 地球物理学报, 58(6): 1954-1962.
曲均浩、刘瑞峰、李金等, 2014, CAP方法反演2014年山东乳山M4.2、M4.0地震震源机制解, 地震工程学报, 36(4): 1076-1080. DOI:10.3969/j.issn.1000-0844.2014.04.1076
曲均浩、王长在、刘方斌等, 2019, 乳山序列地震分布与震源区速度结构的关系, 地震地质, 41(1): 99-118. DOI:10.3969/j.issn.0253-4967.2019.01.007
宋金、蒋海昆, 2009, 序列衰减与余震激发研究进展及应用成果, 地震地质, 31(3): 559-571. DOI:10.3969/j.issn.0253-4967.2009.03.017
谭毅培、邓莉、曹井泉等, 2016, 2015年河北滦县震群发震机理分析, 地球物理学报, 59(11): 4113-4125. DOI:10.6038/cjg20161115
王鹏、郑建常、谭毅培, 2016, 利用重复地震研究山东乳山地区地壳介质波速变化, 地震学报, 38(5): 728-738.
张斌、苏道磊、范建柯等, 2017, 基于自适应量子遗传算法对胶东半岛地区乳山震群重定位及构造特征分析, 地球物理学进展, 32(3): 1080-1088.
郑建常、李冬梅, 2019, 基于误差分布的震源区波速比反演及其应用: 乳山震群源区介质性质变化研究, 地球物理学报, 62(5): 1693-1703.
郑建常、林眉、王鹏等, 2015a, CAP方法反演震源机制的误差分析: 以胶东半岛两次显著中等地震为例, 地球物理学报, 58(2): 453-462.
郑建常、曲利、曲均浩等, 2015b, 乳山震群大小台网地震定位的对比研究, 地震地质, 37(4): 953-965.
郑建常、王鹏、徐长朋等, 2016, 乳山震群震源谱参数的稳健反演, 地球物理学报, 59(11): 4100-4112. DOI:10.6038/cjg20161114
周瑶琪、张振凯、梁文栋等, 2015, 山东东部晚中生代构造-岩浆活动及原型盆地恢复, 地学前缘, 111(1): 137-156.
Akaike H, 1974. A new look at the statistical model identification. In: Parzen E, Tanabe K, Kitagawa G. Selected Papers of Hirotugu Akaike. New York: Springer, 215~222.
Ben-Zion Y, Lyakhovsky V, 2006, Analysis of aftershocks in a lithospheric model with seismogenic zone governed by damage rheology, Geophys J Int, 165(1): 197-210. DOI:10.1111/j.1365-246X.2006.02878.x
Gibbons S J, Ringdal F, 2006, The detection of low magnitude seismic events using array-based waveform correlation, Geophys J Int, 165(1): 149-166. DOI:10.1111/j.1365-246X.2006.02865.x
Hainzl S, Ogata Y, 2005, Detecting fluid signals in seismicity data through statistical earthquake modeling, J Geophys Res Solid Earth, 110(B5): B05S07.
Lei X L, Huang D J, Su J R, et al, 2017, Fault reactivation and earthquakes with magnitudes of up to MW4.7 induced by shale-gas hydraulic fracturing in Sichuan Basin, China, Sci Rep, 7(1): 7971. DOI:10.1038/s41598-017-08557-y
Lei X L, Yu G Z, Ma S L, et al, 2008, Earthquakes induced by water injection at-3km depth within the Rongchang gas field, Chongqing, China, J Geophys Res Solid Earth, 113(B10): B10310. DOI:10.1029/2008JB005604
Meng X F, Peng Z G, Hardebeck J L, 2013, Seismicity around parkfield correlates with static shear stress changes following the 2003 MW6.5 san simeon earthquake, J Geophys Res Solid Earth, 118(7): 3576-3591. DOI:10.1002/jgrb.50271
Narteau C, Byrdina S, Shebalin P, et al, 2009, Common dependence on stress for the two fundamental laws of statistical seismology, Nature, 462(7273): 642-645. DOI:10.1038/nature08553
Ohba T, Daita Y, Sawa T, et al, 2011, Coseismic changes in the chemical composition of volcanic gases from the owakudani geothermal area on hakone volcano, Japan, Bulletin of Volcanology, 73(4): 457-469. DOI:10.1007/s00445-010-0445-9
Ogata Y, 1992, Detection of precursory relative quiescence before great earthquakes through a statistical model, J Geophys Res Solid Earth, 97(B13): 19845-19871. DOI:10.1029/92JB00708
Ogata Y, 2001, Increased probability of large earthquakes near aftershock regions with relative quiescence, J Geophys Res Solid Earth, 106(B5): 8729-8744. DOI:10.1029/2000JB900400
Ogata Y, 1988, Statistical models for earthquake occurrences and residual analysis for point processes, J Am Stat Assoc, 83(401): 9-27. DOI:10.1080/01621459.1988.10478560
Omori F, 1894, On the after-shocks of earthquakes, J Coll Sci Imp Univ Tokyo, 7: 111-200.
Peng Z G, Zhao P, 2009, Migration of early aftershocks following the 2004 parkfield earthquake, Nat Geosci, 2(12): 877-881. DOI:10.1038/ngeo697
Schaff D P, 2008, Semiempirical statistics of correlation-detector performance, Bull Seismol Soc Am, 98(3): 1495-1507. DOI:10.1785/0120060263
Shelly D R, Beroza G C, Ide S, 2007, Non-volcanic tremor and low-frequency earthquake swarms, Nature, 446(7133): 305-307. DOI:10.1038/nature05666
Takahata N, Yokochi R, Nishio Y, et al, 2003, Volatile element isotope systematics at ontake volcano, Japan, Geochemical Journal, 37(3): 299-310. DOI:10.2343/geochemj.37.299
Utsu T, 1961, A statistical study on the occurrence of aftershocks, Geophys Mag, 30: 521-605.
Wang P, Wang B S, 2019, ETAS model analysis on the chang island earthquake swarm in Shandong Province, China, Earthq Res China, 33(4): 617-631.
Wang Q, Jackson D D, Zhuang J C, 2010, Missing links in earthquake clustering models, Geophys Res Lett, 37(21): L21307.
Wiemer S, Wyss M, 2000, Minimum magnitude of completeness in earthquake catalogs: examples from Alaska, the western United States, and Japan, Bull Seismol Soc Amer, 90(4): 859-869. DOI:10.1785/0119990114
Yang W, Ben-Zion Y, 2009, Observational analysis of correlations between aftershock productivities and regional conditions in the context of a damage rheology model, Geophys J Int, 177(2): 481-490. DOI:10.1111/j.1365-246X.2009.04145.x
Yang H F, Zhu L P, Chu R S, 2009, Fault-plane determination of the 18 April 2008 Mount Carmel, Illinois, earthquake by detecting and relocating aftershocks, Bull Seismol Soc Amer, 99(6): 3413-3420. DOI:10.1785/0120090038
Zhang M, Wen L X, 2015, Seismological evidence for a low-yield nuclear test on 12 May 2010 in North Korea, Seismol Res Lett, 86(1): 138-145. DOI:10.1785/02201401170
Zheng J C, Li D M, Li C Q, et al, 2017, Rushan earthquake swarm in eastern China and its indications of fluid-triggered rupture, Earthq Sci, 30(5-6): 239-250. DOI:10.1007/s11589-017-0193-4
Zhuang J C, 2000, Statistical modelling of seismicity patterns before and after the 1990 Oct 5 Cape Palliser earthquake, New Zealand, New Zeal J Geol Geophys, 43(3): 447-460. DOI:10.1080/00288306.2000.9514901