震级上限是指一个地区(潜在震源区或断层)可能发生的最大地震的震级,其在地震危险性分析中具有非常重要的作用[1]. 目前,由于历史地震记录时间不够长,记录周期不能跨越最大地震的复发周期,因此震级上限不能根据地震目录来确定,而往往是根据构造或者地质学原则确定的[2]. 在我国地震危险性分析中,震级上限的确定主要有以下几种途径:对于有记载发生过破坏性地震的地区,通常是根据历史地震来确定震级上限;在未记载到发生过破坏性地震的地区,往往是根据该区域的构造特征,与相邻区域进行比较研究后确定;再者就是根据区域内活动构造的长度、位错量、滑移速率等参数来确定[3]. 上述确定震级上限的方法都需要很好地了解研究区内的地震构造情况. 对于那些地质构造特征不明确的地区,该法则不适用. 因此美国在编制2008美国中东部地震区划图时,统计了全球克拉通地区和板块边缘地区的最大地震的震级,以此来确定美国中东部地区的震级上限[2]. 事实上,这一方法是在地震活动性较弱和地震构造不明确地区的一种替代方法.
我们可将地震的发生看成是一个随机过程,假设地震记录的时间无限长,则记录到的地震事件样本容量可看成是无限大的,这时样本中最大地震的震级趋近于震级上限. 然而在实际情况中我们记录的地震事件频度是有限的,因此需要采用统计模型来估计震级上限. 本文中我们试图采用含有震级上限的截断的震级-频度关系(truncated Gutenberg-Richter Frequency-Nagnitude,简称截断的G-R 关系)[4]来估计震级上限.
我国东北地震区的地震活动性较弱,记录到的最大地震震级为6 3/4 . 区内地震地质工作基础较差. 在先前的研究中人们往往是将地震区内最大地震的震级加上一定的数值并取整得到东北地震区的震级上限,因此震级上限的确定具有较大的不确定性,可能会低估震级上限值. Nueller[5](2010) 研究表明,震级上限在弱地震活动性地区对地震危险性的影响比其在强地震活动性地区的影响要大. 为此,我们尝试使用东北地震区的地震目录,根据截断的G-R 关系式,来估计东北地震区的可能会发生的最大震级. 该法可作为确定研究区震级上限的途径之一.
2 方法一般来说,小震的发生频度要比大震高的多,大小地震的比例关系可用Gutenberg and Richter[6]提出的关系式来表征:
(1) |
上式即为Gutenberg-Richter Frequency-Nagnitude (简称G-R 关系),其中N(m)为震级大于等于m的地震频数,a、b为参数,可用实际地震目录回归得到.
在(1) 式中并未考虑地震震级的上限,即认为震级是可以无限大的,这显然不符合物理规律. Cornell 等提出了考虑震级上限的震级-频度关系,即截断的G-R 关系[4]:
(2) |
则震级概率密度函数为
(3) |
其中N(m)为震级大于等于m的地震频数,N(m0) 为所有震级大于等于震级下限m0 的地震频数,b为表示大地震频数和小地震频数比例关系的一个值.
Mu 为震级上限,通俗地说震级上限即一个地区可能发生的最大地震的震级,从概率上讲,震级上限即表示一个地区发生超过该震级地震的概率趋于零.
可使用非线性最小二乘法或者极大似然法求取N(m0) ,b,Mu. 最优Mu 使截断的G-R 关系标准差最小或似然函数最大. 本文中将详细讲解使用最大似然法来求取Mu.
本文中使用的是不同观测时间段的地震目录.
由于在不同的时间段地震目录完整的起始震级不同,因此,在写出似然函数之前需要将概率密度函数进行简单的归一,归一后的概率密度函数可写成:
(4) |
则似然函数为
(5) |
其中N为总的地震频数,ni为mi+dm的地震频数,dm为震级间隔. 对数似然函数可写为
(6) |
为了求取Mu、b,则应求取对数似然函数的极值,对数似然函数分别对Mu、b求偏导得
(7) |
求上述方程组即可求Mu、b.
3 资料选取及完整性分析本文研究区域为东北地震区(图 1 中多边形),东北地震区主要包括黑龙江、吉林大部、辽宁西北部、内蒙古东北部及河北北部部分地区,这些地区的地震活动性特征较为一致,因此统一划分为东北地震区1) . 需要特别说明的是东北地震区不包含1975 年发生海城地震的区域.
1)中国地震局地质研究所.中国地震区、带划分修改方案简要说明.2008.
选取了东北地震区的地震目录(图 1 中圆圈),震级标度为面波震级MS. 研究所用的历史强震资料来源于《中国历史强震目录(公元前23 世纪—公元1911年)》[7],1912—1990 年的资料来自《中国近代地震目录(公元1912—公元1990 年)》[8],1990 年以后的资料来自《1990年后中国及邻区地震目录》2) .
2)吕悦军.1990年后中国及邻区地震目录.中国地震动参数区划图技术报告.2010.
对于地震活动性来说,地震目录的完整性和可靠性直接关系到最终研究结果的准确性. 为此我们对仪器记录地震目录和历史地震目录做了完整性分析. 对于仪器记录地震,我们采用Frankel[9]在美国地震区划图中使用的方法,即画出地震的G-R 关系曲线和随时间变化的累积频度曲线,观察这两条曲线线性关系是否优良. 图 2即为东北地震1970年以来M≥3.0级地震的G-R 关系曲线和M3.0~3.5 地震的时间累积频度曲线,可以看出线性关系良好,因此可以认为东北地区仪器记录地震从1970 年来是完整的. 对于历史强震目录(M≥4.7) ,据研究表明从1920年来是完整的[10-11].
为了充分利用东北地震记录的历史地震资料,我们希望获得M≥6.0地震的完整时间段,然而由于东北地震区符合M≥6.0的地震比较少,无法采用G-R 关系曲线和时间累积频度曲线来判断其完整性. 为此,我们计算了多个观测时间段的M≥6.0地震的年发生率,并与编制中国地震动参数区划图[12](胡聿贤等,2001) 中东北地区的M≥6.0的地震年发生率做了比较(图 3) . 图 3 中黑色实线为中国地震动参数区划图中东北地震区的G-R关系曲线. 我们将公元1900—2009 年(图 3 中实心三角号)、公元1850—2009 年(图 3 中实心圆)、公元1820—2009年(图 3 中实心方框)三个观测时间段的M≥6.0的地震年发生率与中国地震动参数区划图中的地震年发生率做了比较,可以发现公元1850—2009年观测时间段的M≥6.0 的地震年发生率较为合适,未过高或过低估计6 级以上地震的发生率,因此本文使用的6 级以上地震的观测时间段为1850—2009年(表 1) . 这也避免了因高震级地震的年发生率估计过高而导致计算得到的震级上限偏高的嫌疑.
根据上述描述的方法及选取的地震目录,计算了东北地震区的震级上限Mu 及b值. 式(5) 是关于Mu 和b函数,我们计算了很多格点(b,Mu)上的对数似然函数值lnL,设最大的似然函数值为lnLmax. 我们将所有格点的似然函数值均减去最大的似然函数值,即lnL-lnLmax,称之为相对似然函数值,则似然函数最大处的值为零. 图 4a即是相对似然函数值的等值线图,可以看出在极值附近等值线为闭合椭圆形,这说明采用截断的G-R 模型是能够较好地计算得到震级上限Mu 和b的最优值. 极值点(图中+)对应的纵、横坐标分别为震级上限Mu 和b值的最优值. 图 4b为震级-频度曲线,黑色实线是采用本研究中得到的震级上限和b值所作的东北地震区震级-频度曲线,灰色实线为根据中国地震动参数区划图中的参数所作的曲线,可以看出采用本研究得到的震级上限,计算得到的高震级地震发生率比《中国地震动参数区划图》中的要高. 这说明在中国地震动参数区划图中,可能由于震级上限估计地偏小而低估了高震级地震的年发生率.
一般认为,误差符合正态分布,本研究中Mu 和b相互独立,因此参数Mu 和b的联合误差分布可认为符合自由度为2的关于独立变量lnL-lnLmax 的2 分布. 通过查表可得与需要的置信区间对应的变量值(文中即相对似然函数值),根据该相对似然函数值的等值线可计算出Mu 和b在该置信区间内的范围. 然而在本研究中,由于样本容量相对较小,地震目录中震级相对较大的地震不够多,采用上述方法很难计算出Mu 的置信区间. 为此,我们以东北地震区发生的最大地震的震级M=6.8和中国大陆发生地震的最大震级M=8.5分别作为Mu 的下界和上界,以此作为物理约束边界,然后再根据误差的χ2 分布来计算Mu95%的置信区间(表 2) . 表 2中第三列即为Mu 的95%置信区间.
上述采用地震的面波震级,是根据郭履灿3)的经验公式将地方震级ML 换算而成的. 然而近年来有学者研究认为ML 与MS 差别不大,在实际应用中无需换算[13-15]. 为此,我们直接使用未经震级换算的地震目录计算了Mu 和b(表 3) . 图 5a为关于Mu 和b的相对似然函数值等值线,可以看出在极值附近等值线为闭合椭圆形,这说明采用截断的G-R 模型能够较好地计算得到震级上限Mu 和b的最优值. 极值点(图中+)对应的纵、横坐标分别为震级上限Mu 和b值的最优值. 图 5b为震级-频度曲线,黑色实线是采用本研究中得到的震级上限和b值所作的东北地震区震级-频度曲线,灰色实线为根据中国地震动参数区划中的参数所作的曲线,可以看出在本研究中M≤5.5和M≥6.5地震的发生率均高于中国地震区划图中的值. 本研究中低震级的地震年发生率高于《中国地震动参数区划图》是由于未进行震级的转换造成的,而高震级地震的发生率高于《中国地震动参数区划图》则可能是由于《中国地震动参数区划图》中震级上限的估计偏小造成的.
3)郭履灿.华北地区的地方性震级ML 和面波震级MS 经验关系. 1971.全国地震工作会议资料.1-10.
为了探明本研究中使用的截断的G-R模型的可靠性,我们将截断的G-R 模型与广泛使用的G-R 模型做了比较. 采用Akaike 信息准则[16](Akaike information criterion,简称AIC)来判定模型的合理性:
其中k为模型中参数的个数,L为最大似然函数值. 一般认为AIC 值越小,模型越合理. 表 2 和表 3 中最后一列分别为截断的G-R 模型和G-R 模型的AIC 值,可以看出前者的AIC 值大于后者的AIC 值,因此可以认为在本研究中截断的G-R 模型可能并不优于G-R 模型. 然而我们发现截断的G-R 模型的最大似然值均大于G-R 模型的极大似然值,这说明前者的拟合优度是好于后者的,同时也表明截断的G-R 模型AIC 值偏大是由于其多了一个参数Mu,从而使其复杂度增加造成的. 形成上述现象的主要原因我们认为是由于地震目录不够长,目录中震级较大的地震偏少造成的.
目前人们在计算b值时广泛采用G-R 模型,认为采用该模型计算得到的b值是可信的. 为此,我们比较了分别采用截断的G-R 模型和G-R 模型计算的b值. 从表 2和表 3中的第三列可以看出,采用这两个模型计算出的b值差别很小,相对差值小于5%. 这说明采用截断的G-R 关系不仅能够计算出Mu,也能计算出合理可信的b值. 这也表明截断的G-R 模型在功能上是优于G-R 模型的.
为了进一步验证本文中所用方法的稳定性,在这一部分我们将考虑震级的不确定性,即震级误差. 由于震级误差的存在,我们往往高估了地震的活动率[17-18],原因是:震级误差是符合高斯分布的,高斯分布是一种对称分布,即将一个M=4.8 级地震记录成5.0级地震的概率和将一个M=5.2级地震记录成M=5.0级地震的概率是相同的. 然而地震的震级分布并不是对称分布的,有1个M=5.2 级的地震记录会对应有3个M=4.8级的地震记录(假设b=1) ,这就会导致地震目录中地震震级整体偏高[19]. Felzer在计算美国加州地区的地震活动性参数时主张将每个地震的震级减去其对应的误差ΔM,然后使用修正过的地震目录来计算地震活动性参数[19]. 在我国,有研究表明我国测算的震级总体趋势偏高[20]. 然而各台网在计算地震震级时并不给出震级误差,因此我们无法将每个地震的震级减去其对应的误差. 为此我们把地震震级误差看成是一致的,即将地震震级统一减去一个误差值ΔM. 参考Felzer[19]给出的ΔM值,文中对于1970 年以来记录的地震,将ΔM取0.2,对1970 年以前的地震统一取ΔM为0.33.
表 4为采用震级误差修正后的地震目录计算得到的Mu 和b值. 对于不经震级转换的地震目录,在进行了震级修正后能够较好地得到Mu 和b值. 对于经过震级转换的地震目录,进行震级转换后,由于在Mu 轴方向似然函数没有与Mu 唯一对应的极大值,因此不能求取Mu的最优值,但仍能得到b的最优值.
为了避免因计算方法的因素而造成结果的偶然性,我们又采用最小二乘法进行了计算. 表 5为采用上述四套地震目录计算得到的Mu 和b值. 可以看出Mu 的大小保持在7.5 左右. 这说明了使用东北地震区的地震目录,根据截断的G-R 关系式来求取Mu 的值是可行的,得出的结果是合理可信的.
上述仅分析了在东北地震区采用截断的G-R 关系式来求取Mu 的可靠性. 事实上我们更加关心文中方法的普适性. 高孟潭等4) 采用文中的方法计算了日本东北及邻近海域地区、琉球海沟、马尼拉海沟及美国南加州地区的震级上限,通过与其他研究结果的比较,发现采用文中方法计算得到的震级上限是合理可靠的. 因此可以认为文中采用截断的G-R关系式来求取震级上限的方法具有普适性.
4)高孟潭,徐伟进,周本刚.采用截断的G-R 关系式来求取震级上限的普适性研究.2011
5 结论与讨论本文根据截断的G-R 关系模型,使用东北地震区的地震目录,计算了东北地震区最大可能发生的地震的震级,结果表明东北地震区的震级上限Mu 大于7.0 级,且应为Mu=7.5 左右. 计算中我们考虑了不同震级的转换、震级误差的修正以及计算方法的影响. 最终结果表明,不论采用何种方案进行计算,东北地震区的震级上限值均始终保持在7.5左右,这说明我们采用本文中方法计算得到的东北地震区的震级上限值是合理可信的.
从目前的地震活动性水平来看,东北地震区地震活动性较弱. 有研究表明在地震活动较弱的地区,震级上限对地震危险性的影响是显著的[5]. 在编制《中国地震动参数区划图》[12]中,东北地震区的震级上限被定为7.0级,比较本文的研究结果,该值可能估计的偏小,从而将直接导致高震级地震的年发生率偏低,最终会导致对东北地震区的地震危险性水平估计不足. 因此本文的研究结果可为将来东北地震区的地震危险性研究、地震区划等提供参考. 最新的地质研究表明东北地震区的依兰-伊通断裂带部分段落为晚更新世-全新世活动断裂,并且历史上发生过7.0级以上强震[21]. 该研究结果意味着东北地震区的震级上限大于7.0级,与本文研究结果一致,从地质学上印证了本文研究方法和研究结果的科学性与合理性.
[1] | Ward S N. More on Mmax. Bull. Seism. Soc. Am. , 1997, 87(5): 119-1208. |
[2] | Petersen M D, Frankel A D, Harmsen S C, et al. Documentation for the 2008 Update of the United States National Seismic Hazard Maps. U. S. Geol. Surv. Open-File Rept , 2008, 2008. |
[3] | 冉洪流. 潜在震源区震级上限不确定性研究. 地震学报 , 2009, 31(4): 396–402. Ran H L. Research on uncertainty of upper limit earthquake magnitude in potential seismic source zone. Acta Seismologica Sinica (in Chinese) (in Chinese) , 2009, 31(4): 396-402. |
[4] | Cornell C A, Vanmarcke E H. The major influences on seismic risk. in Proc. of the Fourth World Conference on Earthquake Engineering, January 1969, Santiago, Chile, A-1, 69-93. |
[5] | Mueller C S. The influence of maximum magnitude on seismic-hazard estimates in the central and eastern united states. Bull. Seism. Soc. Am. , 2010, 100(2): 699-711. DOI:10.1785/0120090114 |
[6] | Gutenberg B, Richter C F. Frequency of earthquakes in California. Bull. Seism. Soc. Am. , 1944, 34(4): 185-188. |
[7] | 国家地震局震害防御司. 中国历史强震目录(公元前23世纪-公元1911年). 北京: 地震出版社, 1995. Department of Earthquake Disaster Prevention, National Earthquake Administration. The Catalogue of Chinese Historical Strong Earthquakes (23th Century BC-1911 AD) (in Chinese). Beijing: Seismological Press, 1995. |
[8] | 中国地震局震害防御局. 中国近代地震目录(公元1912-公元1990年). 北京: 科学技术出版社, 1999. Department of Earthquake Disaster Prevention, China Earthquake Administration. The Catalogue of Modern Earthquakes in China (A. D. 1912-A. D. 1990, MS≥4.7) (in Chinese). Beijing: China Science and Technology Publishing House, 1999. |
[9] | Frankel A C, Mueller C, Barnhard T, et al. National Seismic-Hazard Maps: Documentation June 1996, U. S. Geol. Surv. Open-File Rept, 1996, 96-532. |
[10] | 黄玮琼, 李文香, 曹学锋. 中国大陆地震资料完整性研究之二-分区地震资料基本完整的起始年分布图像. 地震学报 , 1994, 16(4): 423–432. Huang W Q, Li W X, Cao X F. Completeness analysis for earthquake data in China Main Land Ⅱ. Acta Seismologica Sinica (in Chinese) (in Chinese) , 1994, 16(4): 423-432. |
[11] | 吴兆营, 薄景山, 刘志平, 等. 东北地震区b值和地震年平均发生率的统计分析. 东北地震研究 , 2005, 21(3): 27–32. Wu Z Y, Bo J S, Liu Z P, et al. Statistical analysis on b value and earthquake yearly average occurrence ratio in north-east seismic zone. Seismological Research of Northeast China (in Chinese) (in Chinese) , 2005, 21(3): 27-32. |
[12] | 胡聿贤, 高孟潭. GB18306-2001《中国地震动参数区划图》宣贯教材. 北京: 中国标准出版社, 2001 . Hu Y X, Gao M T. A Tutorial of the GB 18306-2001, Seismic Ground Motion Parameter Zoning Map of China (in Chinese). Beijing: China Standard Press, 2001 . |
[13] | 张宏志, 刁桂苓, 赵明淳, 等. 不同标度震级关系和台基影响问题探讨. 中国地震 , 2007, 23(2): 141–146. Zhang H Z, Diao G L, Zhao M C, et al. Discussion on relation between different earthquake magnitude scales and effect of seismic station site on magnitude estimation. Earthquake Research in China (in Chinese) (in Chinese) , 2007, 23(2): 141-146. |
[14] | 刘瑞丰, 陈运泰, 任枭, 等. 中国地震台网震级的对比. 地震学报 , 2007, 29(5): 467–476. Liu R F, Chen Y T, Ren X, et al. Comparison between different earthquake magnitudes determined by China seismograph network. Acta Seismologica Sinica (in Chinese) (in Chinese) , 2007, 29(5): 467-476. |
[15] | 汪素云, 王健, 俞言祥, 等. 基于中国地震台网观测报告的ML 与MS经验关系. 中国地震 , 2010, 26(1): 14–22. Wang S Y, Wang J, Yu Y X, et al. The empirical relation between ML and MS based on bulletin of seismological observations of Chinese stations. Earthquake Research in China (in Chinese) (in Chinese) , 2010, 26(1): 14-22. |
[16] | Akaike H. A new look at the statistical model identification. IEEE Transactions on Automatic Control , 1974, 19(6): 716-723. DOI:10.1109/TAC.1974.1100705 |
[17] | Tinti S, Mulargia F. Effects of magnitude uncertainties on estimating the parameters in the gutenberg-richter frequency-magnitude law. Bull. Seis. Soc. Am. , 1985, 75(6): 1681-1697. |
[18] | Rhoades D A. Estimation of the Gutenberg-Richter relation allowing for individual earthquake magnitude uncertainties. Tectonophysics , 1996, 258(1-4s): 71-83. |
[19] | Felzer Karen R. Appendix I-Calculating California seismicity rates, in Earthquake Rate Model 2.2 of the Working Group on California Earthquake Probabilities. U. S. Geol. Surv. Open-File Rept, 2007, 2007-1437.I. |
[20] | 迟振才, 迟天峰. 两种震级标度讨论. 东北地震研究 , 2000, 16(4): 9–15. Chi Z C, Chi T F. Discussion on two magnitude scales. Seismological Research of Northeast China (in Chinese) (in Chinese) , 2000, 16(4): 9-15. |
[21] | 闵伟, 焦德成, 周本刚, 等. 依兰—伊通断裂全新世活动的新发现及其意义. 地震地质 , 2011, 33(1): 141–150. Min W, Jiao D C, Zhou B G, et al. The significance of discovery on Holocene activity on the Yilan-Yitong fault in Northeast China. Seismology and Geology (in Chinese) (in Chinese) , 2011, 33(1): 141-150. |