地球物理学报  2012, Vol. 55 Issue (7): 2341-2352   PDF    
基于ETAS模型对三峡库区流体触发微震活动的定量检测
蒋海昆1 , 宋金1 , 吴琼1 , 李金2 , 曲均浩3     
1. 中国地震台网中心, 北京 100045;
2. 中国地震局地震预测研究所, 北京 100036;
3. 中国地震局地质研究所, 北京 100029
摘要: 针对ETAS模型参数估计方法(MLE)中的初值敏感性问题,提出GA+MLE算法,以GA结果作为MLE算法的初始输入,对结果进行精细计算.通过ETAS模型研究三峡库区微震活动在快速加载及缓慢卸载两种状态下的流体触发、地震自激发及微震活动衰减特征,讨论库水渗透及加卸载过程的可能影响.结果显示:(1)库水快速加载阶段ETAS模型参数μ、α、p及流体触发地震所占比例Rb均显示由小变大、又由大变小的变化过程,但p值的统计差异不显著;在库水缓慢卸载阶段,μRb持续减小;(2)平均来看,库水快速加载阶段流体对微震活动显示较强的外因触发作用,同一条件下序列地震自激发明显增强、衰减相对较慢;水位缓慢卸载阶段,流体对地震活动的触发影响相对较弱,地震自激发不强、衰减相对较快;(3)分阶段来看,蓄水初期库水作用对微震活动的外因触发影响较弱,随库水位的升高及作用时间的增长,流体渗透逐渐发挥作用,孔隙压逐渐增大,流体外因触发作用明显增强,大多数微震活动缘于流体的直接触发(Rb≥95%);足够长的时间之后,由于地下数公里范围在新的载荷及渗透条件下趋于新的平衡,流体渗透影响趋于稳定,孔隙压趋于常数,孔隙压变化趋于0,流体对微震活动的触发作用逐渐减弱.
关键词: 三峡水库      微震活动      ETAS模型      流体触发      渗透      加卸载     
Quantitative investigation of fluid triggering on seismicity in the Three Gorge Reservoir area based on ETAS model
JIANG Hai-Kun1, SONG Jin1, WU Qiong1, LI Jin2, QU Jun-Hao3     
1. China Earthquake Networks Center, Beijing 100045, China;
2. Institute of Earthquake Science, China Earthquake Administration, Beijing 100036, China;
3. Institute of Geology, China Earthquake Administration, Beijing 100029, China
Abstract: To solve the sensitivity problem of the initial value in the MLE algorithm for ETAS model parameter estimation, the GA+MLE method has been proposed in this paper. We take the results of GA (genetic algorithm) as the initial input of MLE, and then estimate the ETAS model parameters more accurately. Under two conditions of quick load and slow unload, the features of fluid triggering, self-generation of earthquakes and sequence decay in the Three Gorge Reservoir area has been studied by the ETAS model. The potential influence of fluid intrusion and load-unload process on seismicity has been discussed roughly. The results show that for quick load phase, the ETAS model parameters, μ, α and p, as well as the Rb, the ratio of fluid-induced earthquakes to total earthquakes, have a positive correlation, all of which showing a similar changing tendency: from small to large, and then to small. But the statistic difference of p value is not very obvious. For slow unload phase, μ and Rb values decrease unceasingly. Averagely, the fluid triggering on seismicity is something strong during the quick load phase. Meanwhile, the self-generation of earthquakes enhances obviously and decay tends to slow down. During the slow unload stage, the fluid triggering on microearthquake activity is relatively weak. In this situation the self-generation is weak and the decay is somewhat quick. For different phases, the influence of fluid triggering on seismicity is weak during the initial phase of water storage. With the time going and with the increase of water level, the pore pressure increases gradually due to the fluid intrusion, and the fluid triggering on seismicity increases obviously. Most of microearthquakes during this phase are caused by the direct fluid triggering (Rb≥95%). After the time long enough, the medium several kilometers below the reservoir becomes saturated under the new condition and the influence of fluid intrusion tends to be stable. In this case, the pore pressure approaches constant and the variation of pore pressure approaches zero, therefore the fluid triggering on microearthquakes decreases gradually..
Key words: Three Gorge Reservoir      Microearthquake activity      ETAS model      Fluid triggering      Fluid intrusion      Load-unload process     
1 引言

一般认为流体渗透而非水体载荷是水库诱发地震活动的主要因素[1],但实际上,微震活动亦显示出与周期性载荷变化相关联的特性[2].由于三峡水库水位变化具有大体相似的时间变化周期(年尺度),因而为对比研究流体渗透及加卸载过程对微震活动的影响提供了可能.

三峡水库位于中国湖北宜昌三斗坪,混凝土重力坝,坝顶高程185 m,最高蓄水位175 m,对应库容393亿立方米.数十年来,对三峡库区地质构造、 地震活动、变形特征等已有非常多的研究.从构造上粗略来看,三峡库区碳酸盐类地层分布广,不同程度地发育有各种岩溶,库区基岩主要为前震旦纪斜长花岗岩,透水性弱[3],有多条NNW、NNE 向断裂穿过库坝区,断层规模不大,较为陡立,倾角在60°~80°之间[4],大多数(约90%)断层胶结良好,其余断层的单位吸水量介于5×10-5 ~l×10-4 m3 · (min· m2)-1 [5].从微震活动来看,自2003 年5 月三峡开始蓄水以来,随库水位的逐渐升高,库区微震活动显著增加.

水体载荷变化可以从库水位记录粗略得到,但渗透及所导致的地下介质孔隙压变化一般无法直接测量.近期的研究显示,传染型余震序列模型(Epidemic Type Aftershock Sequencemodel,简称ETAS模型,详见第2节)参数μ 值的时间变化,与地下孔隙压力增加率相一致[6],这为间接研究地下介质孔隙压变化趋势提供了一种可能.ETAS 模型最初为描述“余震激发余震"现象而提出.ETAS 模型引入无限嵌套的自相似思想,假定任何一次地震均能产生自己的高阶余震,这些高阶余震均遵循修改的Omori公式[7-8].随研究的不断深入,ETAS 模型各参数的物理含义及所表达的物理过程被不断拓展.ETAS模型所表达的地震活动包括两部分:一部分是由于地震自激发(self-generation)所导致的地震群集,另一部分是与地震自激发无关、完全由于外因作用而产生的地震活动.因而,有可能通过ETAS 模型对两部分地震活动进行分离,进而探讨外因触发或自激发对地震活动的贡献[9],该项工作近期的一个重要进展即是流体触发地震活动的定量检测[6, 10].流体引发的地震活动是由于内部应力触发和外部应力变化所导致的地震活动的总和,前者属于地震自激发范畴,后者则可能缘于流体引起的孔隙压力变化.Vogtland震群的研究结果显示,ETAS 模型参数μ 值的时间变化与地下孔隙压力增加率相一致[6],这意味着,利用μ值变化可间接了解地下介质内部孔隙压力的变化.Lei等利用类似的方法对由注水引起的荣昌地震序列进行研究,计算了不同阶段外力触发(流体作用)和内部应力触发(地震间自激发)两种机制所导致地震活动的比例[10];近期基于ETAS模型对一次大暴雨过程诱发(触发) 的微震活动进行了研究,发现由于岩溶构造有利于流体的快速下渗,因而前期干旱前提下的短时间大暴雨过程,导致地下介质孔隙压的快速增加,其对微震活动的触发作用远强于所研究过的水库蓄水作用的影响[11].

据此,本文力图利用三峡水库地震台网记录的微震资料,结合库水位变化,通过ETAS 模型研究加、卸载过程中微震活动的外因触发、地震自激发及微震活动衰减特征,初步探讨水库诱发地震活动与库水渗透及库水加卸载过程的关系.

2 ETAS模型及模型参数求解算法的优化

ETAS模型条件强度函数λ(t)为[7-8]

(1)

右侧求和式为第i次地震对地震发生率的贡献,i遍取所有地震,常数Kpαc适用于所有的iMz为参考震级(必须大于等于最小完备震级Mc).p表征地震活动衰减的快慢,其含义与修改的大森公式中p值相同,p大则衰减快,p小则衰减慢.α 用于度量一次地震产生不同级别高阶“余震"的效率,当Mz取震级下限Mc时,愈大的α 表示愈强的高阶余震激发能力[12].μ(t)最直接的物理含义是剔除地震群集之后单位时间内的理论地震数,因而μ(t)表征由于外力作用而产生的“背景"地震活动率,与地震间自激发过程无关.本文中,μ(t)数值大小表征流体外因触发微震活动的强弱[7].c为数值极小的正常数,与修改的大森公式中c值类似,以往一般认为c值与主震后短时间内数据的完备程度有关,其最主要的作用在于保证(1) 式分母不为0[13],近期的研究显示,c可能还与震源地方的应力状态有关[14].K为与主震震级及bpαcμ 等参数有关的地震活动水平,其中bG-R关系比例系数.

对同一组观测数据,数学上可采用多个统计模型对其进行描述或拟合,即可得到多组符合条件的模型参数.Ogata利用最大似然法估计ETAS 模型参数,并利用AIC 进行最佳模型参数判别[7-8].AIC 既考虑模型对观测数据拟合的优劣,也对无限制地增加参数个数以提高拟合程度的行为进行惩罚. ETAS模型的对数似然值lgL[8]:

(2)

(3)

k为待定参数个数,本文中k=5,取具有最小AIC 的模型参数作为最佳模型.

tbte 为用于AIC 计算的起、止时间,tb 取决于序列资料的记录质量,前提是保证tbtte 时段内地震基本完备,NAICtbtte 时段内的地震数.

上述MLE 算法的计算软件要求预先给定式(1) 式中5个待定参数较为“适当"的初始值,不“适当"的初始值可能导致计算过程的溢出(overflow error)1).与“真值"的接近程度,是初始值“适当"与否的标准.遗憾的是,实际问题中“真值"是待求参数.算例数据的试算结果显示,MLE 算法似乎具有强的初值敏感性,当初始值愈发接近“真值"时,求解结果愈发趋近“真值"且较为稳定;当初始值与“真值"有一定差异时,求解结果与“真值"之间差异会变得较为显著,某个参数初始值的微小改变,会导致求解结果的较大差异2).为合理确定MLE 算法的初始值,本文首先通过遗传算法(Genetic Algorithm,简称GA,下同)在足够大的范围内求取GA 意义下的最优ETAS模型参数,再以之作为MLE 算法的初始输入,对结果进行进一步的精细计算(即GA+ MLE 算法).GA 算法中选择哪一组搜索结果作为进一步产生下一组解的依据,取决于观测值和模型计算值之间的比较,即目标函数F最小.目标函数F定义为:

(4)

式中N(i)为[0,ti]时段内实际发生的地震数,NETAS(i) 是[0,ti]时段内基于ETAS模型计算得到的理论地震数,n为总时段数.式(1) 自tbte 积分可得[tb,te]时段内的NETAS:

(5)

式(5) 右侧两项分别为时间[tb,te]内流体触发的地震数Nb 及地震自激发产生的地震数Nt.据此,外因触发地震所占比例Pb=Nb/NETAS,自激发地震所占比例Pt=1-Pb.式(5) 第二项直接采用Simpson 数值积分方法[15]进行计算.

1) Ogata Y.Statistical Analysis of Seismicity-Updated Version (SASeis2006).2006:12-19

2) 蒋海昆等.2010年7月,水库地震序列统计特征及与水库加卸载关系研究.国家“十一五"科技支撑计划项目子专题中期检查材料.

GA 作为数值估计的一种优化算法,具有全域搜索、无需求导、效率较高而又普遍适用的特点,其最终结果并不是搜索目标的真值,但却是规定条件下对真值的最优化逼近.GA 参数的选择对收敛速度有一定影响,仿石耀霖(1994)推荐的做法3),规定种群数目保持为64,交配及变异概率分别为0.90 及0.02,随机种子数为5.目标函数的尺度变换取指数形式,并规定保留每代的最优模型.目标函数尺度变换的目的在于提高搜索效率,使算法既不在早期限于局部极值、又能在后期诸多较优模型中放大它们之间的差别,以搜索差异很小的许多模型中的最优模型.指数形式的尺度变换即起到将目标函数接近最小值的部分拉伸、接近最大值的部分压缩这样一种作用.依据试算结果,GA 迭代300次之后结果已基本稳定,因而规定GA 最高迭代300次,GA 搜索范围及迭代步长列于表 1.

3) 石耀霖.遗传算法及其在地球科学中的应用,中国科学院中国科学技术大学研究生院讲义.1994.

表 1 余震序列E1AS模型参数确定的GA搜索范围及迭代步长 Table 1 Searching range and iterative step for estimation of ETAS model parameters by GA

GA 和GA+ MLE 分别给出各自“最优"的ETAS模型参数,其中后者以GA 结果作为MLE 算法的初始输入.以模型输出与实际观测数据的吻合程度对GA 及GA+MLE 结果进行评价,吻合程度高者作为最终的ETAS模型参数.即令

(6)

式中N(i)为实际累积地震频次,若等时间间隔统计,则i为时间间隔序数(i=1,2,…,n),n为总时段数;若依地震顺序依次计算,则N(i)=ii=1,2,3,……,nn为地震次数.NGA(i)和NMLE(i)分别为利用GA 及GA+MLE 算法估计ETAS 模型参数、再经式(5) 计算得到的累积地震频次.若δMLE< δGA则取GA+ MLE 结果为最优模型参数,反之则取GA 结果为最优模型参数.

以[N(i),NGA(i)]及[N(i),NMLE(i)](i=1,2,…,n)在二维线性坐标上作图,理论上若NGA(i) 或NMLE(i)与N(i)非常接近,则会得到一条过原点、斜率为1 的直线.据此,分别以数据集[N(i),NGA(i)]及[N(i),NMLE(i)](i=1,2,…,n)进行最小二乘拟合y=a+bx,其中y={N(i)},x={NGA(i)} 或x={NMLE(i)},理论上好的模型参数应有a→0 及b→1.图 1给出一个确定最优ETAS模型参数的实例,该实例中GA+ MLE 结果明显优于GA 结果.换言之,以GA 搜索300次迭代得到的结果,作为GA+MLE 算法的输入,两者相结合可以在一定程度提高ETAS模型参数估计的效率及准确性.

图 1 确定最优ETAS模型参数的一个实例 (三峡研究区,2008-09-25-10-31,ML≥0.7,MZ = 0.7,地震样本数171) (a)累积地震频次随时间的变化,依地震顺序逐次计算.曲线为实际地震累积频次,“〇"和“□”为基于ETAS模型计算得到的累积地震频次,ETAS模型参数分别由GA及GA+MLE方法得到;(b)基于ETAS模型计算的累积频次与实际地震累积频次之间的关系,余同(a). Fig. 1 A case for determination of the optimal ETAS model parameters. Sanxia region,Sep.25,2008 - Oct.31,2008.ML≥0.7,Mz = 0.7,the sample number is 171 (a) Cumulated earthquake frequency vs time,calculated one by one according to the orders of the earthquakes.Curve is the real cumulated earthquake frequency,“〇" and “□” are cumulated earthquake frequency giving by ETAS model with GA and GA+MLE respectively; (b) Cumulated earthquake frequency by ETAS model vs real cumulated earthquake frequency,the meanings of symbols are same as (a).
3 资料及加、卸载时段划分

三峡工程早期地震监测采用有人值守地震台的组网方式,最早于1958年开始监测.至三峡工程下闸蓄水的2003 年,已有45 年观测天然地震的时间[16].为保证一致的观测精度,本文使用三峡数字地震台网(简称三峡台网,下同)的观测资料.三峡台网包括24个高增益子台,其中库坝区4 个子台为24位数采、1s周期的反馈地震计,其余为16 位数采、1s周期的开环地震计.空间上库首区三斗坪至庙河段布设8个子台,平均台距约10km;香溪至碚石段布设14个子台,平均台距15~20km;台网西部网沿,在碚石南、北两侧布设2 个子台(参见图 2).三峡台网于2000年8月开始试运行,2001年8 月正式运行,库坝区定位精度1km,其余区域小于2km[17].

图 2 三峡地区ML≥0.0 地震分布(2003-05-19-2009-12-31) (★大顼位置,台站,▲2008年11月26日胡家坪Ml4.6级地震) Fig. 2 Earthquake distribution with ML≥0.0 in Sanxia reservoir and surrounding area (from May 19,2003 to 31 Dec.,2009)

2003年5 月19 日三峡水库蓄水以来,库区及附近小震活动明显增加(参见图 4a),空间上主要分布于巴东东瀑口、官渡口及楠木园一带(图 2).关于三峡台网的微震监测能力,以往研究认为在坝址至楠木园段(90×30km2)为ML≥0.5级,在楠木园至碚石段监测能力优于ML(0.7 级[16-17].针对图 2 所示的地震密集区域,取窗长12个月、滑动步长1个月,由G-R关系依次确定每个时间窗内的最小完备震级Mc,结果如图 3 所示,图 3 还给出相应的地震频次及G-R关系b值,可见地震频次随蓄水过程(时间)而逐渐增加,每个窗口内参与计算的地震数均大于300次.b值总体上较为稳定,数值变化范围非常小,介于0.91~1.18 之间,但似有随时间略为降低的趋势.最小完备震级Mc 在不同时段略有差异:2006年上半年之前为ML0.6,2006 年下半年之后大部分时间段为ML0.7,少数时段为ML0.8.对台网的进一步调查发现,自蓄水以来台网运行状况基本稳定,尤其2006年前后监测台网及相应台站并无改变,因而2006 年下半年之后G-R关系线性段转折点Mc(最小完备震级据此确定)由ML0.6改变为ML0.7、个别时段还达到ML0.8(图 3图 4),可能缘于该时段研究区地震震级结构的真实变化,而不是检测能力降低、地震漏记等引起.换言之,最小完备震级Mc 可统一确定为ML0.6.

图 3 三峡地区最小完备震级Mc(○)、G-R关系b值(☆)及地震频次(-) (统计时间窗及时间滑动步长分别为12个月及1个月) Fig. 3 Minimum completeness magnitude Mc(○), b value of G-R relationship (☆) and earthquake frequency (-) in Sanxia reservoir and surrounding area
图 4 确定最小完备震级Mc的两个实例(上图为频次直方图,下图为G-R关系) Fig. 4 Two cases for determination of the minimum completeness magnitude Mc. (upper plot is histogram of frequency,lower plot is G-R relationship) (a) 2004 年1-12 月,b=1.07,Mc = 0.6 (b) 2008 年1-12 月,b=1.03,Mc = 0.8.

关于三峡库区微震震源深度,到目前为止未有系统的精确结果.由三峡台网测定结果来看(图 5),绝大多数ML≥0.6地震震源深度介于3~10km 之间,其中又以5~7km 深度最为密集.能够确定深度的5872次地震中,80%的深度小于9km,90%的深度小于11km,95%的深度小于13km,99%的深度小于18km.2138 次地震的精确测定结果显示,不同区域震源深度存在明显的差异,例如仙女山断裂附近平均震源深度为5.6km,而巴东地区平均震源深度则仅为2.6km[18].

图 5 三峡地区ML≥0.6地震深度分布(左测坐标) 及震源深度分布的累计概率分布(右侧坐标) (图 2圈定范围内2003-06-01-2009-12-31期间共发生ML≥ 0.6地震6311次,其中439次无法确定深度,图中表示为0 km) Fig. 5 Focal depth distribution with ML≥0.6 in Sanxia reservoir and surrounding area (scale is in left axis) and its accumulative probability curve (scale is in right axis)

图 6a图 2所示区域ML≥0.0地震M-t图及水位曲线,图 6bML≥0.6级地震的月频次.可见自2003年5月下旬蓄水以来,微震活动显著增强,2006年9 月完成156 m 蓄水后,频次持续增加,2008年10月达到173m 水位前后,频次达到最高,11月22日在库首区胡家坪发生ML4.6级地震(震源深度8km),这是三峡水库自2003年蓄水以来发生的最大地震.震源机制解揭示的胡家坪发生ML4.6级地震的发震断层为NNW 向,逆冲左旋错动4),是在库水荷载与库水下渗的共同作用下,沿仙女山断裂发生的构造型水库诱发地震[19].

4) 湖北省地震局,2008年11月22日秭归县屈原镇Ms4.1 地震宏观烈度调查报告,2008年12月

图 6 三峡水库蓄水前后库区及附近ML≥0.0微震M-t图、 库水位变化(a)及ML≥0.6月频次统计(b) Fig. 6 M-t plot with ML≥0.0,water level curve (a) and monthly histogram with ML≥0.6 (b) in Sanxia reservoir and surrounding area

三峡水库蓄水过程分为几个阶段(参见图 6a),2003年5月19日至6月15日,完成135m 蓄水目标,库容达140×108 m3;之后直至2006年9月蓄水季之前,水位介于135~138 m 之间,变化幅度较小;2006年9月20日至10月27日,完成156m 蓄水目标,水位提高约20m;2008年9月开始的蓄水季,预定目标是最高水位175m,但由于上游来水不足,11月11日达到173m 后开始缓慢下降;2009年夏季过后于9 月15 日重新开始蓄水,目标仍然是175m,同样由于上游来水不足、下游干旱严重等原因,水位于11月24日达到171m 后开始缓慢下降. 由于2007 年之后每年的最低水位一般维持在145m,因而2008、2009年的两次蓄水过程,水位上升幅度约为30m.可见自2006年以来,在水位总体持续增高的前提下,三峡库水位呈现“快速上升(1 个月左右)-缓慢下降(10 个月左右)"的循环变化过程,这对库水渗透及对库区浅层地壳的载荷作用将产生一定影响,也为研究库区微震活动与库水加卸载和渗透过程之间的可能关系提供了条件.水位“快速上升"、“缓慢下降"具体时段划分如表 2所列.

表 2 与蓄水有关的三峡微震活动统计时段划分 Table 2 Statistical stages for microseismic activity in Shanxia region concerning with water level variation in Sanxia reservoir and surrounding area
4 基于ETAS 模型参数变化对流体触发微震活动过程的研究 4.1 流体对不同震级下限地震活动的触发影响

几次水位快速上升阶段(R1、R2、R3 及R4 阶段)的μ 值变化可分为大体一致的三组(图 7),一组震级下限M0 介于0.6~0.8之间,随蓄水过程的发展,不同加载阶段μ 值变化显著;另一组M0 介于1.2~1.4 之间,μ 值较小且变化相对平稳;第三组M0 介于0.9~1.1 之间,μ 值及其变化幅度介于前二者之间.有两个显著特点,一是不同震级下限μ 值具有类似的变化趋势,但小地震μ 值的变化幅度大,较大地震μ 值的变化幅度小;二是同一时期较大地震的μ 值相对较小、较小地震的μ 值相对较大.

图 7 水位快速上升阶段ETAS模型参数μ值的变化(横坐标R1、R2、R3、R4分别与表 2中水位快速上升阶段相对应,右侧数字为震级下限M0) Fig. 7 μ values in different quick raising stages of the water level, μ is one of the ETAS model parameters

由于μ 值大小表征流体对地震触发作用的强弱,因而前一个特点意味着,尽管同是快速加载状态,但不同时期流体对微震活动的外因触发强度有差异,与渗透过程的持续发展过程有关;后一个特点则给出一个非常重要的认识,即库水渗透及加卸载过程对小地震有较强的外因诱发(触发)作用,对较大地震而言则流体诱发(触发)作用较弱.这从不同阶段平均μ 值随震级下限的变化亦可清晰看出,无论加载还是卸载过程,平均μ 值随震级下限M0 的提高而逐渐减小(图 8).

图 8 库水加、卸载阶段ETAS模型参数μ值平均值随震级下限M0风的变化 (●水位快速上升阶段,★水位缓慢下降阶段) Fig. 8 ● value variation with the lower magnitude M0 in loading and unloading stages,★ is one of the ETAS model parameters
4.2 不同加、卸载阶段的流体触发特征

为保证有足够的统计样本及较为显著的变化特征,参考图 7,选择M0≥0.6、M0≥0.7、M0≥0.8三组地震数据,计算不同加、卸载阶段的最优ETAS 模型参数及流体触发地震所占比例,研究库水位快速上升(加载)和缓慢下降(卸载)阶段,ETAS 模型参数的数值特征,粗略讨论其物理含义.

(1) 库水快速上升(加载)阶段

定性来看,库水快速加载阶段ETAS模型参数μαp及流体触发地震所占比例Rb 之间的变化趋势定性一致,均有一个由小变大、又由大变小的过程(图 9),其中p值变化的统计差异不显著.

图 9 库水位快速上升阶段ETAS模型参数特征 横坐标R1、R2、R3、R4分别与表 2中水位快速上升阶段相对应;右侧数字为震级下限Mc 〇MC>0. 7,AMC>0.8;★及误差棒:平均值及相应的均方差范围 Fig. 9 Characteristics of ETAS model parameters in different quick raising stages of the water level

三峡水库2003年5-6月完成135m 蓄水目标之后、直至2006年9月之前,水位介于135~138m 之间,变幅不大,库水载荷作用相对稳定.2006年9 月开始的第一次水位快速上升(R1阶段),μαp值较小(图 9(a-c)),均值分别约为0.54、0.66 及0.83,流体触发地震所占比例平均约为0.47(图 8d),表明水位首次快速上升阶段,库水作用对微震活动的外因触发影响较弱,系统自激发不强,微震活动亦衰减较慢.

2007年9月开始的R2阶段,水位从145 m 上升到156m.这一时期尽管未超过R1阶段曾达到过的最高水位,但由μ 值表征的流体外因触发作用明显增强,平均μ 值升高至2.59(图 9a),95% 以上的微震活动缘于流体的直接触发(Rb>0.95;图 9d),地震自激发能力略有增强,平均α 值约为0.80(图 9b),地震活动衰减较R1 阶段有所加快,平均p值约为1.22(图 7c).这表明,随库水位的升高,流体渗透逐渐发挥作用,导致孔隙压逐渐增大、介质裂隙强度降低,因而R2阶段流体作用较前期R1阶段明显增强.

2008年9 月开始的第三次水位快速上升(R3 阶段),库水位从145 m 上升到173 m,水位快速升高约28m,较此前曾达到过的最高水位156m 也升高了17m.三峡库区流体渗透导致孔隙压变化的一维模拟研究结果显示5),这一阶段库区地下数公里范围内孔隙压趋于常数,孔隙压变化趋于0.因而流体渗透影响趋于稳定,平均μ 值与R2 阶段基本一致,约为2.57(图 9a),仍保持较强的流体触发作用但变化已经不大.这一时期微震活动的显著增加(参见图 6b),可能缘于地震自激发能力的明显增强,因为该时段α 值明显升高至约2.11(图 9b),而流体触发地震所占比例Rb 则由此前R2阶段的0.95降低至约0.54(图 9d);由于序列自激发地震数量的明显增多,即流体触发无关的“天然"地震增多,使得这一时期微震活动衰减继续加快,平均p值增大至1.73 (图 9c).

5)蒋海昆等.《水库地震序列统计特征与水库加卸载关系研究》,国家科技支撑计划项目专题(2008BAC38B03-01)研究报告第5 章:水库诱发地震活动的定量检测及可能机理讨论.2011年8月.

2009 年9 月开始的R4 阶段,库水位再次从145m 上升到171 m,升高约26 m,但未超过2008 年173m的最高水位.这一时期地下固定深度处孔隙压变化微小,流体渗透影响减弱.因而这一时段μ 值下降明显,平均μ 值约为0.70(图 9a);与此相一致,流体触发地震所占比例降低至约0.20(图 9d); 地震自激发能力也明显降低,衰减有所减慢,平均αp值分别约为0.51和0.76(图 9(bc)).

(2) 库水位缓慢下降(卸载)阶段

图 10为水位缓慢下降阶段ETAS模型参数μαp及流体触发地震所占比例Rb的数值分布.D1、 D2和D3 分别是前述各快速加载阶段(R1、R2 及R3)之后、持续较长时期的水位缓慢下降阶段(表 2).定性来看,在库水缓慢卸载阶段,μRb 正相关,持续减小;考虑参数计算的离散程度之后,αp的变化趋势不明显,但粗略来看,较大的μ 值对应较小的α 及较大的p,较小的μ 值对应较大的α 及较小的p.

图 10 库水位缓慢下降阶段ETAS模型参数特征横坐标D1、D2、D3分别与表 2中水位缓慢下降段相对应;右侧数字为震级下限McMc≥0.6,□ Mc≥0.7,△ Mc≥0.8;★及误差棒:平均值及相应的均方差范围 Fig. 10 Characteristics of ETAS model parameters in different slow decreasing stages of the water level

D1阶段库水位从156m 缓慢下降至145m,水位降低约16m.尽管载荷降低,但由于渗透过程的持续发展,这一时期流体对地震活动的外因触发作用较此前的快速加载R1阶段反而明显增强,平均μ 值约等于1.98(图 10a);与此相对应,绝大多数微震活动缘于流体的直接触发(Rb≈0.95,图 10d).由于“天然"地震所占比例不高,这一阶段地震自激发仍然不强,平均α 约等于0.43(图 10b),但可能由于力学卸载的影响,微震活动衰减却非常快,平均p值约为2.24(图 10c).

D2阶段库水位仍然从156m 缓慢下降至145m,此阶段平均μ 值降至约1.29(图 10a),表明流体对微震活动的直接触发较D1 阶段有所减弱,比此前的快速上升阶段(R2)更是明显降低;与此相对应,流体触发地震所占比例也降低至0.33(图 10d).这一阶段地震自激发仍然不高,与此前D1 阶段大体类似,平均α 约等于0.38(图 10b).序列衰减较D1 阶段明显减缓,平均p值约为1.08(图 10c),已接近天然地震的平均衰减速率[20].

D3阶段库水位从173m 缓慢下降至145m,水位降低约28m.这一缓慢卸载过程流体对地震活动直接的外因触发作用继续减弱,平均μ 值降至约0.81(图 10a),流体触发地震所占比例Rb 降低至0.19(图 10d).地震自激发尽管比此前水位快速上升阶段(R3)明显降低,但较之D2阶段则有所提高,平均α 约等于1.01(图 10b).由于约80%的微震活动缘于地震自激发(Rb≈0.2),因而微震活动衰减仍显示与“天然"地震类似的衰减特征,平均p值约为1.22(图 10c).

4.3 加、卸载阶段流体触发平均统计特征的对比

依据前述各加、卸载阶段ETAS模型参数及流体触发地震所占比例的计算结果,表 3 给出所有加载阶段及所有卸载阶段ETAS 模型参数的平均结果.由表 3可见,快速加载阶段,μα 相对较高,p相对较小;缓慢卸载阶段则正好相反.这表明,平均来看,库水快速加载阶段流体对微震活动具有较强的外因触发作用,同一条件下地震自激发也明显增强,衰减相对较慢;水位缓慢卸载阶段,流体对地震活动的触发影响相对较弱,序列地震自激发强度不高,衰减相对较快.由表 3 还可看出,无论快速加载还是缓慢卸载,流体外因触发地震所占比例Rb 差异不大.这可能意味着,由于地震自激发与流体触发定性正相关(参见图 10、11),因而尽管不同加、卸载阶段Rb 数值有差异,但就所有时段平均结果来看,流体直接触发与地震自激发所产生地震大体上各占一半.需要指出的是,由于不同加、卸载阶段ETAS 模型参数变化较大(参见图 910),因而表 3所列参数波动范围(均方差SD)较大.

表 3 加、卸载阶段ETAS模型参数统计特征(ML≥0.6) Table 3 Average characteristics of ETAS model parameters in loading and unloading stages (ML≥0.6)
5 结论

(1) 针对已有ETAS 模型参数估计MLE 算法中的初值敏感性问题,提出GA+ MLE 算法.首先利用GA 算法,在足够大的数据空间中求取GA 意义下的最优ETAS 模型参数,以之作为MLE 算法的初始输入,对结果进行精细计算.以模型输出与实际观测数据的吻合程度对模型优劣进行评判.结果显示,GA+MLE 算法得到的最佳ETAS模型与实际资料具有较高的吻合程度.

(2) 三峡库区自2003年5月下旬蓄水以来微震活动显著增强,2006 年9 月完成156m 蓄水后微震频次持续增加,2008年10月达到173m 水位前后频次达到最高,2008年11月22日发生蓄水以来的最大地震(ML4.6).三峡库区及附近微震震源深度介于3~10km 之间,以5~7km 深度最为密集.三峡水库蓄水过程具有阶段性特征,在总体持续增高的前提下,自2006年以来库水位变化呈现“快速上升(1个月左右)-缓慢下降(10 个月)"的循环加卸载过程.

需要指出的是,2003 年5 月水位快速上升之后,微震活动增加明显(参见图 6).但之后直至2006 年9月之前,库水位变化平稳,这一时期尽管微震活动持续不断,但频次并未有明显增加.2006 年9 月之后,水库蓄、排水过程出现明显的年变,库水加、卸载作用明显增强,微震频次亦明显增加,2008年10 月达到173m 水位前后,频次达到最高,11月22日发生蓄水以来的最大地震---胡家坪ML4.6 级地震.这表明,从现象上看,水库蓄水导致库区地震活动显著增加,在此基础上,库水位的明显起伏及进一步提高会导致库区地震活动在前期基础上进一步增强.关于这一现象的机理6),限于篇幅将另文讨论.

6) 蒋海昆等.2011 年8 月.《水库地震序列统计特征与水库加卸载关系研究》,国家科技支撑计划项目专题(2008BAC38B03-01) 研究报告第5 章:水库诱发地震活动的定量检测及可能机理讨论.

(3) 库水快速加载阶段,ETAS模型参数μαp及流体触发地震所占比例Rb 的变化趋势大体一致,随时间的推移,先由小变大、再由大变小的过程,但p值变化的统计差异不显著.这意味着,外因触发地震活动强度随渗透过程的发展而变化,在库水快速加载阶段,地震自激发及微震活动衰减与外因触发强弱正相关.在库水缓慢卸载阶段,μRb 持续减小,αp的变化趋势不明显,但粗略来看,较大的μ 值对应较小的α 及较大的p、较小的μ 值对应较大的α 及较小的p.即在缓慢卸载阶段,流体触发作用较强时,序列自激发及衰减相对较弱,流体触发作用较弱时,序列自激发及序列衰减相对较强,这与快速加载阶段正好相反.

(4) 由于库水载荷增加相对于地体压力来说非常小,而流体渗透又是一个缓慢、滞后的过程,因而在加载初期,库水作用(渗透及加卸载)对微震活动的外因触发影响较弱、系统自激发不强、序列衰减较慢,μαpRb 值较小.随库水位升高及蓄水时间的持续增长,流体渗透逐渐发挥作用,导致孔隙压逐渐增大,裂隙强度降低,由μ 值表征的流体外因触发作用明显增强,这一阶段流体触发地震所占比例Rb 明显增大,大多数微震活动缘于流体的直接触发(Rb≥95%),地震自激发能力也随之增强,序列衰减略有加快(αp增大).足够长的时间之后,由于库区地下数公里范围内介质在新的载荷条件下逐渐趋于新的平衡,孔隙压趋于常数,孔隙压变化趋于0,流体外因触发作用逐渐减弱6)(具体模拟计算结果将另文叙述),μ 值明显降低,流体触发地震所占比例Rb 亦降至约20%左右.

(5) 就快速加载及缓慢卸载阶段平均来看,流体在快速加载阶段对微震活动具有较强的外因触发作用,同一条件下序列地震自激发也明显增强、衰减相对较慢;缓慢卸载阶段则流体的触发影响相对较弱,这一条件下地震间自激发能力不高,微震活动衰减相对较快.

需要指出的是,较浅部位孔隙压变化对库水位变化的时间响应快、变化幅度大,较深处则响应滞后时间长、变化幅度小[11, 21],因而即使在水位缓慢下降阶段,深部孔隙压仍可能处于增加状态(孔隙压变化滞后),但总体上孔隙压变化幅度随时间逐渐变小[21].若仅考虑水库蓄水之后库水位的总体抬升,而不顾及其间的水位变化的话,蓄水导致的孔隙压变化是单调增加的,但增幅随时间及深度仍然逐渐变小、直至趋于0[11],6).因而,蓄水初期水库地震的显著增加主要与流体渗透导致的孔隙压增加有关; 同时,不同阶段ETAS 模型参数的变化特征,可能也受库水位快速增加或缓慢降低这一准年变化过程所导致的孔隙压波动影响.

(6) 由于μ 表征库水外因触发地震的强弱,因而从机理上不难理解Rbμ 值定性正相关这一事实.依据ETAS 模型的原始定义,任何一次地震均激发自身的高阶余震,因而在库水加载阶段,若流体直接触发的地震活动增强,这些地震激发的高阶余震数量亦会增加,即地震自激发增强,反之亦然,由此亦可理解库水加载阶段μα 之间定性正相关这一观测事实.在足够长的时间之后,渗透影响趋于稳定,此时力学影响将有所现象,而卸载是与加载相反的力学作用过程,因而卸载条件下μαp之间的关系与加载条件下相反.但在蓄水早期阶段,由于渗透导致孔隙压变化的影响,这一认识并不成立.

参考文献
[1] Gough D I, Gough W I. Time dependence and trigger mechanisms for the Kariba (Rhodesia) earthquakes. Eng. Geol. , 1976, 10(2-4): 211-217. DOI:10.1016/0013-7952(76)90021-1
[2] El Hariri M, Abercrombie R E, Rowe C A, et al. The role of fluids in triggering earthquakes: observations from reservoir induced seismicity in Brazil. Geophys. J. Int. , 2010, 181: 1566-1574.
[3] 于品清, 周明礼, 杨淑贤. 未来三峡水库区诱发岩溶型水库地震的可能震级. 华南地震 , 1992, 12(2): 81–87. Yu P Q, Zhou M L, Yang S X. Possible magnitude of induced Karst earthquake in the reservoir which will be built in the coming years. South China Journal of Seismology (in Chinese) (in Chinese) , 1992, 12(2): 81-87.
[4] 高士钧, 王青云, 宋慧珍, 等. 长江三峡地区地壳应力场与地震. 北京: 地震出版社, 1992 : 4 -9. Gao S J, Wang Q Y, Song H Z, et al. Crustal Stress Field and Earthquakes in Changjiang Sanxia Region (in Chinese). Beijing: Earthquake Press, 1992 : 4 -9.
[5] 于品清. 从水文地质条件探讨未来三峡水库发生构造型水库地震的可能性. 华南地震 , 1993, 13(1): 76–83. Yu P Q. Discussion on the possible tectonic earthquake in Changjiang three-gorge reservoir area in relation to hydrogeologic conditions. South China Journal of Seismology (in Chinese) (in Chinese) , 1993, 13(1): 76-83.
[6] Hainzl S, Ogata Y. Detecting fluid signals in seismicity data through statistical earthquake modeling. J. Geophys. Res. , 2005, 110: B05S07. DOI:10.1029/2004JB003247
[7] Ogata Y. Statistical models for earthquake occurrences and residual analysis for point processes. J. Am. Stat. Assoc. , 1988, 83(401): 9-27. DOI:10.1080/01621459.1988.10478560
[8] Ogata Y. Statistical model for standard seismicity and detection of anomalies by residual analysis. Tectonophysics , 1989, 169(1-3): 159-174. DOI:10.1016/0040-1951(89)90191-1
[9] Ogata Y. Space-time point-process models for earthquake occurrences. Ann. Inst. Stat. Math. , 1998, 50(2): 379-402. DOI:10.1023/A:1003403601725
[10] Lei X L, Yu G Z, Ma S L, et al. Earthquakes induced by water injection at -3km depth within the Rongchang gas field, Chongqing. China. J Geophys. Res. , 2008, 113: B10310. DOI:10.1029/2008JB005604
[11] 蒋海昆, 杨马陵, 孙学军, 等. 暴雨触发局部地震活动的一个典型例子: 2010年6月广西凌云-凤山交界3级震群活动. 地球物理学报 , 2007, 54(10): 1778–1786. Jiang H K, Yang M L, Sun X J, et al. A typical example of locally triggered seismicity in the boundary area of Lingyun and Fengshan following the large rainfall event of June 2010. Chinese J. Geophys. (in Chinese) , 2007, 54(10): 1778-1786.
[12] Ogata Y. Increased probability of large earthquakes near aftershock regions with relative quiescence. J. Geophys. Res. , 2001, 106(B5): 8729-8744. DOI:10.1029/2000JB900400
[13] 蒋海昆, 郑建常, 吴琼, 等. 传染型余震序列模型震后早期参数特征及其地震学意义. 地球物理学报 , 2007, 50(6): 1778–1786. Jiang H K, Deng J C, Wu J, et al. Earlier statistical features of ETAS model parameters and their seismological meanings. Chinese J. Geophys. (in Chinese) , 2007, 50(6): 1778-1786.
[14] Clement Narteau, Svetlana Byrdina, Peter Shebalin, et al. Common dependence on stress for the two fundamental laws of statistical seismology. Nature , 2009, 462(7273): 642-646. DOI:10.1038/nature08553
[15] 易大义, 蒋叔豪, 李有法. 数值方法. 杭州: 浙江科学技术出版社, 1984 : 58 -113. Yi D Y, Jiang S H, Li Y F. Numerical Method (in Chinese). Hangzhou: Zhejiang Science and Technology Press, 1984 : 58 -113.
[16] 杨晓源, 张绿茵. 三峡水库诱发地震监测简介. 四川地震 , 2001(2): 1–5. Yang X Y, Zhang L Y. Introduction on monitoring system of earthquakes induced by Sanxia Reservoir. Earthquake Research In Sichuan (in Chinese) (in Chinese) , 2001(2): 1-5.
[17] 谢蓉华, 宋澄, 邵玉坪. 我国水库地震监测数字化进展. 水电站设计 , 2008, 24(24): 79–82. Xie R H, Song C, Shao Y P. Progress of digital monitoring on induced earthquakes in Chinese mainland. Design of Power Dam (in Chinese) (in Chinese) , 2008, 24(24): 79-82.
[18] 廖武林, 张丽芬, 姚运生. 三峡水库地震活动特征研究. 地震地质 , 2009, 31(4): 707–714. Liao W L, Zhang L F, Yao Y S. Characteristics of seismicity in the three gorges reservoir area. Seismology and Geology (in Chinese) (in Chinese) , 2009, 31(4): 707-714.
[19] 车用太, 陈俊华, 张莉芬, 等. 长江三峡工程库首区胡家坪MS4.1水库诱发地震研究. 地震 , 2009, 29(4): 1–13. Che Y T, Chen J H, Zhang L F, et al. Study of the reservoir induced Hujiaping MS4.1 Earthquake in the three gorges dam area. Earthquake (in Chinese) (in Chinese) , 2009, 29(4): 1-13.
[20] Utsu T, Ogata Y, Matsuura R S. The Centenary of the Omori formula for a decay law of aftershock activity. J. Phys. Earth. , 1995, 43(1): 1-33. DOI:10.4294/jpe1952.43.1
[21] Xinglin L. Possible roles of the Zipingpu reservoir in triggering the 2008 Wenchuan earthquake. Journal of Asian Earth Sciences , 2011, 40(4): 844-854. DOI:10.1016/j.jseaes.2010.05.004