2. 中国地震局地球物理研究所, 北京 100081
2. Institute of Geophysics, China Earthquake Administration, Beijing 100081, China
对区域地震危险性进行评估时需要使用该区域的地震目录,尤其是强震的地震目录.地震危险性分析使用概率模型,需要大量的样本以保证其可信性.但是由于现代地震观测时间短,强震数量很少,有些地区甚至没有发生强震的记录,但是具有发生强震的危险性.使用统计地震学模型分析区域危险性时,其准确性难以验证,因此需要开展理论地震活动性模拟方面的研究.
目前国际上使用的地震活动模拟模型很多,比较有代表性的有:(1) 利用速度与状态参数关系准则求解速度量与状态量组成的偏微分方程组的数值解 (Dieterich, 1994),并且以速度大小作为破裂是否发生的标度,从而产生理论地震目录 (Dieterich, 1995, Dieterich et al., 2015),此模拟方法可以得到滑动速率,从而可以求得地震的破裂长度;然而这种方法计算量较大,计算耗时较长,模拟得到单一地震的发生过程,较实际观测存在较大偏离,还没有真正应用到地震危险性分析中;(2) 利用断层物理参数模拟特定的实际断层的地震活动性以及应力场变化规律 (Ben-Zion et al., 2003, Ben-Zion, 2012),此方法将断层的蠕滑区域和地震区域分开赋予强度值,从而这种方法可以细致模拟断层不同区域的地震特性,但是对于断层的物性的空间分布的认识程度要求较高,目前该方法用于圣安德烈斯断层的地震活动性研究中 (Zigone et al., 2015);(3) 基于细胞自动机模型的地震活动性模拟方法,通过空间与时间的离散化进行地震模拟,但是并没有考虑应力传递机制 (李锰和杨峰, 2011).
多断层模型地震活动性模拟方法引入断层间传递机制 (Robinson and Benites, 1995),在此基础上考虑破裂传递和自动愈合效应,使之更接近实际破裂过程 (Robinson et al., 2005),并将改进的模型用于研究新西兰惠灵顿地区的地震活动性 (Zhou et al., 2006, Robinson et al., 2011).此方法被应用于模拟中国川西地区的地震活动性,分析该区域的强震分布特性,并得出龙门山断裂带具有超强地震 (7.5级) 发生的危险的结论 (周仕勇, 2008).
本研究发展并改进了Robinson和Zhou等 (2005)建立的模型,使新模型可以使用GPS反演的断层滑动速率 (Shen et al., 2005, Tong et al., 2014) 作为震间剪应力增长源.模拟方法中将断层模型由单层模型转化为两层模型,上半层为地震区域,下半层为滑动区域.除此之外改进了求解算法,使计算效率大大提高.改进的模拟方法用于模拟中国太原地区的地震活动性.首先,我们将模拟结果与实际观测资料进行对比,验证此模拟方法用于太原地区的合理性;而后,通过模拟结果分析得到太原地区的区域地震危险性特征;最后,本研究分析了整个区域和每条断层上的地震随时间分布的统计分布特征,并分析了用Poisson分布来近似地震序列置信度.
2 地震活动性模拟方法本研究采用基于库伦破裂准则的准静态地震活动性模拟的方法进行地震活动性模拟,通过比较促使断层滑动的剪应力和阻碍断层滑动的摩擦阻力的大小来判断断层是否发生滑动,即地震是否发生.本研究中剪应力的变化来源于断层的滑动速率,断层的滑动速率是由GPS数据反演得到的.根据GPS反演模型,本研究中断层模型分为两层,如图 1所示,上层为锁住部分,即地震区域;下层为滑动部分.锁住区域震间无滑移,积累剪应力至其达到强度,地震发生时该区域发生同震滑移并释放应力;滑动区域震间以匀速相对滑移,没有应力积累,地震发生时该区域不发生同震滑移.区域地震危险性分析着重研究地震时产生的同震位移和破裂面积的大小,因此对地震区域进行网格划分,并研究其应力变化规律.
具体模拟过程如下:
(1) 首先,对断层的地震区域进行网格划分.破裂单元基本尺度的大小根据想要模拟的最小震级来选取.网格尺度越小,可以模拟的地震的最小震级就越小,但是由于网格数量增加,计算耗时也越长.
(2) 计算每个子单元的剪应力.震间网格区域的剪应力和强度改变主要来源于滑动区域的匀速滑动,地震发生时网格区域的剪力和强度改变主要来源于地震区域网格的破裂.其中各单元剪切应力的走向分量与倾向分量的计算如下式:
(1) |
式中τs(i, m, t), τd(i, m, t) 分别表示的断层面上 (i, m) 号网格上的剪切应力在t时刻沿着走向和倾向的分量,其中i表示研究区域断层的编号,m是i号断层上划分的子单元的编号;τs0和τd0分别表示 (i, m) 号网格在地震活动性模型中设置在子断层单元上的初始剪应力在走向和倾向方向的分量;(n, l) 网格表示对 (i, m) 号网格的应力状态产生影响的网格,其中n是研究区域断层的编号,l是n号断层上划分的子单元的编号;Nf是研究区断层的数量 (本研究为8);Ln是n号断层上分割的断层单元的总数量;us(n, l, t) 和ud(n, l, t) 是t时刻传递到 (i, m) 断层单元的 (n, l) 断层单元发生的脆性破裂产生的位错;Kss(i, m; n, l),Ksd(i, m; n, l),Kds(i, m; n, l),Kdd(i, m; n, l) 分别是 (n, l) 单元发生的走向和倾向方向的位错对 (i, m) 单元走向和倾向方向的剪切应力影响的传递系数 (triggering coefficient),可由弹性位错理论求得 (Okada, 1992);Vsn和Vdn是n号断层锁定深度以下的滑动部分的滑动速率分量,由GPS观测资料反演得到;K′ ss(i, m; n),K′ sd(i, m; n),K′ ds(i, m; n),K′ dd(i, m; n) 分别表示n号断层的下层滑动部分走向和倾向方向的滑动对 (i, m) 单元走向和倾向方向的剪切应力影响的传递系数.每一个单元的剪应力的计算公式如下:
(2) |
(3) 抑制断层滑动的摩擦力,又可称之为断层的静态强度.每一个网格上的静态强度由下式确定:
(3) |
式中 (i, m) 仍为单元的编号,Sstatic(i, m, t) 为断层单元的摩擦阻力即静态强度;μstatic为静态摩擦系数;L为断层单元处的静岩压 (压缩方向为负,拉张方向为正);δσn是构造运动加载以及研究区域t时刻前所有断层单元的破裂产生的位移造成的 (i, m) 单元上正应力的变化,计算方法与断层上的剪应力的计算方法相同;Pbase为空隙压的恒定部分;δP(i, m, t) 是所有断层单元的破裂的位移造成 (i, m) 单元的空隙压力的变化.当断层发生破裂时其强度变化为动态破裂强度,动态破裂强度为
(4) |
式中μdyn=dyn·μstatic为动态摩擦系数.
(4) 由上两步可以计算每一个子单元每一时刻的剪应力和强度,根据库伦破裂准则确定每一条断层最先破裂的时间,从中找出最先破裂的子单元,该单元的位置就是地震的初始破裂点.计算每个网格的破裂时间时我们采用牛顿法与二分法相结合的计算方法,首先使用二分法找到包含破裂时间的单调区间,在单调区间内使用牛顿法找到破点时间点.牛顿法与二分法相结合的计算方法相比于之前模型中使用二分法的计算方法,计算步数减少了一个量级,大大提高了计算效率.
破裂单元的强度降为动态强度,剪应力也相应降低,剪应力降为
(5) |
这里arr为捕获因子,取值在0.0~dyn之间,这样才能保证破裂发生后剪应力低于强度.
(5) 本研究考虑了破裂自发愈合过程,这样使模拟过程更接近实际的动态过程.子单元发生破裂以后,该子单元的强度立即由静态破裂强度降低为动态破裂强度,然而破裂单元在环境应力作用下会愈合,经过一定的愈合时间强度恢复到静态破裂强度 (Somerville et al., 1999).本研究中破裂愈合时间设定为3.0 s.
(6) 地震区域网格发生破裂时,每次破裂都将产生扰动应力场,影响所有单元的应力值.扰动应力场可能引起其他单元发生破裂,因此要重复步骤 (2)-(5),搜索其他破裂单元.若扰动应力场引起同断层上的其他单元发生破裂,宏观上表现为断层破裂的传播或地震的增长;若扰动应力场引起其他断层上的单元,使其发生破裂,宏观上表现为研究区多地震事件同时发生.本研究考虑到应力场的传播过程,单元破裂产生的应力场以一定速度传递到其他子单元,本研究考虑到这一传递过程造成的时间延迟,并选取传递速度为模型介质中S波的传播速度的0.9倍.直至每一次破裂传递到所有单元,而无新的破裂发生,一次地震事件视为结束.
(7) 一次地震事件中初始破裂点与最终破裂单元的破裂时间差即为地震过程的持续时间.地震事件的地震矩计算公式如下:
(6) |
G为介质的剪切模量,A为子单元的面积,un为破裂断层上第n个子单元在破裂过程中发生的累积位错.地震事件的矩震级MW可由经验公式来确定:
(7) |
本研究以上述步骤完成地震活动性的模拟,由断层的构造模型和滑动速率计算得出区域理论合成地震目录.
3 太原地区地震活动性模拟应用我们改进的模型可以将GPS反演结果作为构造加载进行输入,断层上的剪应力积累来自下层滑动区域的构造运动.我们选取经度35°-38°,纬度110°-113°的太原盆地地区作为研究区域.山西省的省会太原位于该区域内,是山西省政治经济的中心.该地区历史上曾发生过破坏严重的地震,造成了重大人员财产损失,因此对该区域的研究具有很重要的意义.该区域内部的地质构造主要有北部的太原盆地和南部的临沂盆地,地处于鄂尔多斯周缘断裂区域 (杨国华等, 2002).本研究以这两个盆地周缘断层为研究基础,主要考虑太原盆地东西两侧的交城断裂和太古断裂,太原盆地和临汾盆地之间的霍山山前断裂、和贾村断裂,以及临汾盆地周缘的浮山山前断裂、罗云山山前断裂和娥眉台地北缘断裂.
利用高彬等人根据地震活动性建立起断层的几何模型 (高彬等, 2016)(图 2),根据地表GPS资料和水准观测资料反演断层之上的滑动速率,表 1中的结果是由高彬等人利用GPS反演方法计算得到的结果 (高彬, 2016).表 1中滑动区域下界的确定主要根据GPS反演使用的滑动区域下界,设定为90 km,由于1号断层和2号断层在50 km深处相交,避免两断层相互截断1号断层和2号断层的滑动下界设定为50 km;3号断层在70 km深处插到5号断层下方,所以3号断层滑动下界设定为70 km.
模拟不仅需要断层的物理信息,还需要一些介质物理参数,表 2中列出了我们根据地震学观测结果估计的介质物理参数.利用以上参数按照我们改进过的模拟方法对断层模型进行活动性的模拟,进而分析太原地区的地震活动性.
本研究利用上述断层模型和介质参数模拟了20000年的地震目录.我们采用随机初始应力设置,为了避免初始应力对模拟结果的影响,我们截断了前3000年的地震目录,使用后17000年的地震目录进行统计分析.图 3展示了模拟的17000年的地震目录的时间震级序列和太原地区1970年以来的实际观测记录的震级时间图像.
图 4展示了不同滑动速率下的模拟结果的震级频度曲线与1970年以来的时间观测记录的震级频度曲线的对比图.由于缺乏可靠的跨断层形变观测资料,论文中所应用的远场GPS观测资料所估计的断层滑动速率具有较大的不确定性 (高彬, 2016).图 4(a-c)分别为直接使用GPS反演的滑动速率,反演结果×50%和反演结果×25%作为构造加载得到的模拟结果.从对比模拟结果可以看出,如果直接使用GPS反演得到的滑动速率进行模拟,5.1~5.3级的模拟地震发生频率比实际观测记录的发生频率大5倍.随着滑动速率的降低,模拟结果的地震发生频率也随之降低,当我们使用反演滑动速率×25%作为构造加载进行模拟时5.1~5.3级的模拟地震发生频率与实际观测较为接近.如果地震活动性在时间上是稳定的,我们是可以使用这一模拟方法对GPS反演得到的滑动速率进行校正.从图像可以看出模拟结果在强震区域与古登堡-里克特分布有较大差异,而服从特征地震分布 (Schwartz and Coppersmith, 1984, Yamashita, 1993, Parsons and Geist, 2009, Kagan et al., 2012).该区域1970年以来实际观测记录的b值为1.0950,3次模拟结果的b值分别为0.8964,1.1514和0.9515,都与实际观测的b值相近.
由模拟结果可以得到太原地区区域地震危险性空间分布,图 5为模拟地震的空间分布图.图 5显示,交城断裂与太古断裂之间的太原盆地中6级以上地震发生概率较高,并且7级以上地震发生概率也相对较高;除此之外,模拟结果显示研究区域南部的峨眉台断裂带与罗云山断裂带6级以上地震的发生概率也较高;这些模拟结果与实际观测资料比较吻合.由物理模拟方法得到的模拟地震只能分布在模拟使用的断层面上,因此断层模型的平面假设和模型中断层数量的局限性导致了模拟地震分布的集中性,这与实际的观测资料有一定的差异,也是纯物理模型的一点局限.
为了更直观地分析模拟结果的空间分布特征,图 6显示了模拟地震的空间概率密度分布,其中图 6a显示了6级以上模拟结果的概率密度分布.由图 6a可以看出,娥眉台地北缘断裂和罗云山山前断裂6级以上地震危险性较高,太原盆地两侧的交城断裂和太古断裂的危险性次之,霍山山前断裂、浮山山前断裂和贾村断裂的危险性较低.这与实际地震的分布比较相似,但也存在一定的差异.实际观测的地震目录显示交城断裂、太古断裂和罗云山山前断裂6级以上地震活动性较强,与模拟结果相同;然而,模拟得到的峨眉台地北缘断裂的地震危险性比实际观测的地震发生率高了很多.存在这一差异,主要是由于峨眉台地北缘断裂的断层长度较长,而且GPS模拟结果显示走向滑动速率相对较大.
图 6b显示了7级以上模拟结果的概率密度分布.由图 6b可以看出太原盆地地区的交城断裂上7级以上地震危险性相对较高,罗云山山前断裂、峨眉台地北缘断裂和浮山山前断裂具有一定的7级以上地震危险性.实际观测资料显示研究区域内,只有罗云山断裂北段上发生过7级以上地震,与我们的模拟有一定的差距.由于模拟地震序列时间为20000年,而实际资料只有1000年;而且我们的模拟结果显示交城断裂的7级以上地震复现周期为4000年,峨眉台地北缘断裂的7级以上复现周期为10000年,罗云山山前断裂南段和浮山山前断裂的7级以上地震复现周期为20000年;因此7级以上地震的模拟结果的区域分布比实际地震的区域分布要广一些,这一结果是合理的.
图 7为6级以上模拟地震的时间间隔分布以及与其λ值相同的Poisson分布的时间间隔的对比图,模拟得到的地震序列的年平均发生率λ=0.0129.由图 7可以看出整个地区的地震序列近似服从Poisson分布,对于模拟结果使用λ=0.0129的Poisson分布进行χ2检验,检验结果其符合Poisson分布的置信水平为99.0%.实际观测资料为1038年以来的地震目录,实际观测资料的6级以上的地震的年平均发生率为λ=0.0082,两者结果相差约1/3.
图 8为八条断层的6级以上地震时间间隔分布图像.由图 8中a和h表示的第1条断层和第8条断层的时间间隔分布与Poisson分布有一定的相似性,而第2~7条断层的时间间隔分布与Poisson分布有很大的差别.表 3中展示了每条断层上6级以上地震年平均发生率λ与该断层用Poisson分布拟合的置信度.虽然整个地区的6级以上地震的地震活动性可以用Poisson分布近似,但是单断层的6级以上地震活动性并不都能用Poisson模型很好地近似,其中6级以上地震年平均发生率相对较高的断层的地震间隔与Poisson分布有较高的相似性,但是年平均发生率相对较低的断层的地震间隔与Poisson有很大的差异,主要是因为这些断层6级以上地震的时间间隔样本也比较少,并不能形成很有规律的统计分布.
本文进一步发展了Robinson等 (2011)和Zhou等 (2006)的地震活动性模拟模型,使得模型可以使用GPS反演得到的断层滑动速率进行地震活动性模拟,并且用双断层模型对模型进行了检验.最后我们把模型应用到太原地区,对太原地区进行了地震活动性模拟.
我们将改进的模拟方法运用到中国太原地区,结果显示我们的模拟结果与1897年以来的实际地震记录的震级频度曲线吻合得很好,并且模拟的b值与实际观测记录的b值相差在合理范围内.由于观测记录时间有限,观测记录的最大震级只有5.5级,因此震级时间曲线只有4.7级至5.5级可以进行对比.对于这个问题我们可以再缩小网格的尺度使得模拟震级减小,由于我们使用的弹性位错理论在距离缩小时可能会产生一定的歧义结果,所以这一解决方法还有待进一步的研究;另一个解决方法是扩大模拟区域,使得可以得到的实际观测记录的样本量更多,这需要我们建立大区域的断层模型,并且需要更大区域内GPS反演的数据结果.
模拟结果在空间分布上与实际观测资料相吻合,只是在南部的峨眉台断裂带区域模拟结果高于实际观测资料.我们分析此结果产生的原因是峨眉台断裂带的长度较长,这个断层的GPS反演的走向滑动速率也相对较大,导致该断层的模拟地震活动性很高.因此我们有理由认为我们得到的该断层的GPS滑动速率结果相对于实际滑动速率是偏大的.
整个区域的6级地震活动性模拟结果的时间间隔分布近似服从Poisson分布,采用Poisson分布拟合置信度高达99.0%.但是单断层结果并不都能够服从Poisson分布,6级以上地震高发断层的时间间隔分布近似服从Poisson分布,而地震发生率低的断层则不服从Poisson分布.此结果表明,地震样本量比较大时,区域地震的时间分布可以用Poisson分布估计.
最后我们讨论了模拟方法的稳定性,实验证明,初始应力场的选取对模拟结果的统计特征影响不大.为了分析随机初始应力场对模拟结果的影响,我们选取了几组随机初始应力场,并且对模拟结果的空间概率密度图像进行了分析.图 9和图 10分别展示了4组完全不同的初始应力场的6级以上模拟地震和7级以上模拟地震的空间概率密度分布.从该图像我们可以得出,初始应力场的设置对于模拟结果的空间分布影响不大.区域地震活动性的空间分布主要取决于断层分布、断层滑动速率大小以及断层尺度等因素.我们同时计算了不同初始应力场下的每年平均发生概率.计算结果表明,不同初始应力场下的6级以上地震的年平均发生率在0.011次/a和0.013次/a之间波动.该结果证明,初始应力场模拟结果的时间分布影响不大.
我们感谢中国地震局地球物理研究所吴建平研究员提出分断层讨论地震活动性的宝贵建议,感谢两位审稿人对这篇文章提出建设性的意见.
Ben-Zion Y, Eneva M, Liu Y F. 2003. Large earthquake cycles and intermittent criticality on heterogeneous faults due to evolving stress and seismicity. Journal of Geophysical Research:Solid Earth, 108(B6): 2307. DOI:10.1029/2002JB002121 | |
Ben-Zion Y. 2012. Episodic tremor and slip on a frictional interface with critical zero weakening in elastic solid. Geophysical Journal International, 189(2): 1159-1168. DOI:10.1111/j.1365-246X.2012.05422.x | |
Dieterich J. 1994. A constitutive law for rate of earthquake production and its application to earthquake clustering. J. Geophys. Res., 99(B2): 2601-2618. DOI:10.1029/93JB02581 | |
Dieterich J H. 1995. Earthquake simulations with time-dependent nucleation and long-range interactions. Nonlinear Processes in Geophysics, 2(3-4): 109-120. DOI:10.5194/npg-2-109-1995 | |
Dieterich J H, Richards-Dinger K B, Kroll K A. 2015. Modeling injection-induced seismicity with the physics-based earthquake simulator RSQSim. Seismological Research Letters, 86(4): 1102-1109. DOI:10.1785/0220150057 | |
Gao B, Zhou S Y, Jiang C S. 2016. Estimate of dip angles of faults around Ordos block based on earthquakes. Chinese Journal of Geophysics, 59(7): 2444-2452. DOI:10.6038/cjg20160711 | |
Kagan Y Y, Jackson D D, Geller R J. 2012. Characteristic earthquake model, 1884-2011, R.I.P. Seismological Research Letters, 83(6): 951-953. DOI:10.1785/0220120107 | |
Li M, Yang F. 2011. Simulation study on heterogeneous cellular automation influencing seismicity factors. Inland Earthquake, 25(3): 205-214. | |
Okada Y. 1992. Internal deformation due to shear and tensile fault in a half space. Bulletin of the Seismological Society of America, 92(2): 1018-1040. | |
Parsons T, Geist E L. 2009. Is there a basis for preferring characteristic earthquakes over a Gutenberg-Richter distribution in probabilistic earthquake forecasting?. Bulletin of the Seismological Society of America, 99(3): 2012-2019. DOI:10.1785/0120080069 | |
Robinson R, Benites R. 1995. Synthetic seismicity models of multiple interacting faults. Journal of Geophysical Research, 1001(B9): 18229-18238. DOI:10.1029/95JB01569 | |
Robinson R, Zhou S Y, Johnston S, et al. 2005. Precursory accelerating seismic moment release (AMR) in a synthetic seismicity catalog:a preliminary study. Geophysical Research Letters, 32(7): L07309. DOI:10.1029/2005GL022576 | |
Robinson R, van Dissen R, Litchfield N. 2011. Using synthetic seismicity to evaluate seismic hazard in the wellington region, New Zealand. Geophysical Journal International, 187(1): 510-528. DOI:10.1111/j.1365-246X.2011.05161.X | |
Schwartz D P, Coppersmith K J. 1984. Fault behavior and characteristic earthquakes:examples from the Wasatch and San Andreas Fault Zones. Journal of Geophysical Research:Solid Earth, 89(B7): 5681-5698. DOI:10.1029/JB089iB07p05681 | |
Shen Z K, Lü J N, Wang M, et al. 2005. Contemporary crustal deformation around the southeast borderland of the Tibetan Plateau. J. Geophys. Res, 110(B11): B11409. DOI:10.1029/2004JB003421 | |
Somerville P, Irikura K, Graves R, et al. 1999. Characterizing crustal earthquake slip models for the prediction of strong ground motion. Seismological Research Letters, 70(1): 59-80. DOI:10.1785/gssrl.70.1.59 | |
Tong X P, Smith-Konter B, Sandwell D T. 2014. Is there a discrepancy between geological and geodetic slip rates along the San Andreas Fault System?. Journal of Geophysical Research:Solid Earth, 119(3): 2518-2538. DOI:10.1002/2013JB010765 | |
Yamashita T. 1993. Application of fracture mechanics to the simulation of seismicity and recurrence of characteristic earthquakes on a fault. Journal of Geophysical Research, 98(B7): 12019-12032. DOI:10.1029/93JB00713 | |
Yang G H, Wang M, Han Y P, et al. 2002. The movement of trend and dynamic character in Shanxi fault zone. Earthquake Research in China, 18(2): 148-156. | |
Zhou S Y, Johnston S, Robinson R, et al. 2006. Tests of the precursory accelerating moment release model using a synthetic seismicity model for wellington, New Zealand. Journal of Geophysical Research, 111(B5): B05308. DOI:10.1029/2005JB003720 | |
Zhou S Y. 2008. Seismicity simulation in western sichuan of China based on the fault interactions and its implication on the estimation of the regional earthquake risk. Chinese Journal of Geophysics, 51(1): 165-174. | |
Zigone D, Ben-zion Y, Campillo M. 2015. Modelling non-volcanic tremor, slow slip events and large earthquakes in the Guerrero subduction zone (Mexico) with space-variable frictional weakening and creep. Geophysical Journal International, 202(1): 653-669. DOI:10.1093/gji/ggv174 | |
高彬, 周仕勇, 蒋长胜. 2016. 基于地震活动性资料估计鄂尔多斯块体周缘构造断层面倾角. 地球物理学报, 59(7): 2444–2452. DOI:10.6038/cjg20160711 | |
李锰, 杨峰. 2011. 影响地震活动性因素的非均匀细胞自动机模拟研究. 内陆地震, 25(3): 205–214. | |
杨国华, 王敏, 韩月萍, 等. 2002. 山西断裂带活动趋势与动态特征. 中国地震, 18(2): 148–156. | |
周仕勇. 2008. 川西及邻近地区地震活动性模拟和断层间相互作用研究. 地球物理学报, 51(1): 165–174. | |