根据临界点理论,地震的孕育过程是一种临界现象(Vere-Jones,1977),这一临界现象的典型特征是大地震之前往往存在地震活动的相关长度增长及主震周围矩释放加速(AMR)现象。Bowman等(1998)研究发现,几次大震前震源区地震能量会出现明显幂率加速释放过程,符合临界理论的标度率,进而使用震前累积应变释放曲线的幂率拟合和线性拟合的均方根残差比最小来定量描述AMR的“临界”时、空尺度。随后人们对世界上其他地区的震例进行了类似的研究,Sykes和Jaume(1990)研究表明美国旧金山几次强震前都有5.0级左右中等地震活动的矩释放加速过程; Ellsworth 等(1981)指出,强震前小震发生率并没有明显变化,而中等地震的累积应变能呈指数增长。我国地震学家在长期预报研究实践中也将AMR现象作为强震预测的主要前兆之一(刘蒲雄,1982)。杨文政和马丽(1999)用类似方法对中国震例开展研究,发现三分之二以上的中强震在震前较长时间均出现了地震矩释放加速的现象; Jiang和Wu(2005)的研究指出,中国大陆有60%的M_S≥6.0地震前存在稳定的地震矩释放加速特征,但在多时、空尺度的二维空间,一些地震前的矩释放表现出加、减速交替出现的复杂现象。更有研究表明,主震前临界区范围内地震活动呈现加速释放,而在较小的孕震区范围则出现地震活动的减弱态势(苗青壮,蒋海昆,2008; Papazachos et al.,2000,2006)。

华北地区是一个内部结构复杂和断裂活动强烈的地区。按其主要地震构造带和区域地震构造格局,可划分为两个一级活动地块:鄂尔多斯活动块体和太行冀鲁西活动块体,其分界带为山西地堑系,其中太行冀鲁西地块又可分为太行山块体、华北平原块体和胶辽块体,其分界带为太行山山前断裂带和郯庐断裂带,还有华北北部的阴山—燕山块体,其分界带为燕山—渤海断裂带(Xu,Gao,2000; Huang,Guo,1998)。华北地区的强烈地震几乎全部发生在这些活动块体的边界带上,其中山西地震带发生M_S>7地震8次,华北平原带3次,郯庐地震带3次以及燕山—渤海地震带8次,给我们的国计民生带来了巨大的损失。因此,华北地区一直是我国地震监测防预的重点监视区和研究地震前兆机理、检验地震预测理论和探索地震数值预测方法的一个重要的“实验场地”。

本文在以往研究的基础上,以1989年以来华北地区12次M_S≥5.0中强地震为例,利用地震矩释放程度的时—空扫描方法,系统验证华北地区中强震前地震能量的释放特征,进而讨论AMR现象在华北地区是否具有普遍性的前兆意义。

“地震破裂时间分析”法是利用给定区域的前兆性地震活动来预测未来地震的发震时间和强度的一种中长期预测方法,其函数定义一般可表示为(Bufe,Varnes,1993)

S(t)=A+B(t_f-t)^m.(1)

式中,S(t)是[0,t]时段内地震活动的累积应变能,t_f为预期地震的破裂时间,t为前兆性地震事件的发震时间,t<t_f。A、B、m为模型参数,其中A在物理上表示[0,t_f]时段内地震活动的累积Benioff应变能。若A已知,则由A减去实际已经释放的Benioff应变,剩余则是预期在t_f时刻要发生目标地震的应变能(蒋海昆等,2009)。幂指数m是加速释放模型中的标度参数,直接反映Benioff应变释放的状态和释放速率的大小。m<1表示Benioff应变加速释放,m越小则加速释放特征越显著; m>1表征 Benioff应变为减速释放状态,相应的m越大减速释放特征越明显; m=1表示Benioff应变匀速释放。

对临界区域的估算依据是使区域内的样本数据能被幂率关系拟合得最好,“最优搜索办法”(Bowan et al.,1998)是以强震震中为圆心,用不同半径范围内的观测数据连续计算曲率参数C值(震前Benioff应变释放的幂律拟和均方差与线性拟和均方差之比),参数C为最小且幂指数m<0.8的圆型区域被认为是主震的临界区域。基于上述理论,杨文政和马丽(1999)利用空间临界圆重叠率最高的方法来寻找未来强震的临界区域; 也有研究认为这种“最优搜索”仅仅是对区域数据本身的拟合,显著存在的问题是相同区域构造和相近震级的中强震之间的临界时、空尺度差异较大 Kossobokovv.2007.Earthquake forecast/prediction:problem of verification,accuracy and limitations.ICTP.。因此蒋长胜和吴忠良(2009)仅利用m值作为描述震前矩加速释放“程度”的参量,在时—空—地震序列的完整性震级组成的三维空间内,使用固定时间尺度和多空间尺度系统研究了2008年3月21日新疆于田M_S7.3强震前AMR现象存在的客观性。研究表明矩加速释放特征在时间上有逐渐逼近临界状态的趋势,空间上与实际震中位置也有较好的对应的关系。这种利用m值的时—空演化来研究地震孕育的物理过程和寻找未来中强震的临界区域的方法称为地震矩释放程度的时—空扫描。这种方法实质上是假定实际震中和发震时刻未知,使用固定的时—空尺度,在时—空坐标上扫描震前的矩释放程度,考察研究区中强震前AMR现象与实际震中和发震时刻之间的相关性。

笔者以1989年以来华北地区12次M_S≥5.0地震作为研究对象,选取(34°～42°N,108°～123°E)的空间范围为研究区域。研究区中强地震的活动表现出明显的分区、交替、成组发生的特征,在同一地震活跃期中发生的中强地震,往往时间相连,地点相近,构成一个典型的地震成组集中活动区,而与此同时其他地区几乎没有6级以上地震发生。1900年以来,华北地区中强地震的活动经历了5次周期性轮回过程(孙加林,章瑞,2007),目前正处于第五活跃期的末期或第五个平静期(图1)。第五活跃期以1989年山西大同和阳高间M_S5.9地震为起始,至今已发生了12次M_S≥5.0的中强地震。如图2所示,这12次地震主要分布在燕山—渤海地震带和山西地震带,而在活动地块内部及郯庐地震带、华北平原地震带和黄海地震带仅有个别地震发生,在空间上分区成组分布的特征明显。

图1 华北地区1900年以来MS≥5.0地震M-t图(图中标出5次中强震跃时段)
Fig.1 M-t plot of MS≥5.0 earthquakes in north China since 1900(five periods of moderate strong earthquake activity are marked in the Figure)

图2 1989年以来华北地区MS≥5.0地震分布(底图中黑色线段为构造线(据邓启东等,2003),左下角给出研究区示意图)
Fig.2 Distribution of MS≥5.0 earthquakes in north China since 1989(the black line in the base map show the main tectonics line(by Deng et al., 2003),schematic diagram of the research area in the low left corner)

图3 1989年以来华北地区完整性震级MC时序图
Fig.3 Time sequence of completeness magnitude MC in north China since 1989

本文在矩释放时—空扫描模型参数的选取中,结合前人给出的有关AMR现象中震级与“临界”时、空尺度的关系(Ben-Zion,Lyakhovsky,2002; 彭克银,尹祥础,2003)分别计算了12次震例前的前兆性地震序列在时—空(t-r)二维空间内m值的分布特征,然后用归一化的方法给出华北地区中强震前具有普遍意义的AMR现象的“临界”时—空尺度。近年来在许多地震预测模型的研究中对所使用地震目录是否删除余震都采取了折衷的做法,同时给出删除余震前、后的两种预测结果。在以往对AMR现象的研究中,蒋长胜(2004)研究指出余震删除与否,对m值计算结果的影响并不明显; 李宇彤和蒋长胜(2010)对东北地区中强震的研究表明,使用删除余震目录得到的时—空尺度与前人研究结果更接近。因此本文使用G-K法(Gardner,Knopoff,1974)删除余震后的地震目录,将(1)式中的t_f固定为研究震例的实际发震时刻; t为震前2～14 a,以1 a为

表1 1989年以来华北地区MS≥5.0地震震源参数
Tab.1 Source parameters of MS≥5.0 earthquakes in the north China since 1989

注:-表示Benioff累积应变减速释放; +表示Benioff累积应变加速释放; 0表示Benioff累积应变匀速释放.

图4 删除余震后的12次中强震前m值在时—空二维空间的分布
(a)1991年山西忻州MS5.1地震;(b)1991年河北唐山MS5.1地震;(c)1992年黄海MS5.3地震;(d)1995年山东
苍山MS5.2地震;(e)1995年河北唐山MS5.0地震;(f)2006年河北文安MS5.1地震;(g)1989年山西
大同MS5.7地震;(h)1991年山西大同MS5.8地震;(i)1996年内蒙包头MS6.4地震;(j)1998年
河北张北MS6.2地震;(k)1999年河北张北MS5.6地震;(l)1999年山西浑源MS5.6地震
Fig.4 Two dimensional spatial and temporal distribution of scaling coefficient m before 12 moderate

(a)Shanxi Xizhou M_S5.1 earthquake in 1991;(b)Hebei Tangshan M_S5.1 earthquake in 1991;(c)Huanghai M_S5.3 earthquake

in 1992;(d)Shandong Cangshan M_S5.2 earthquake in 1995;(e)Hebei Tangshan M_S5.0 earthquake in 1995;(f)Hebei Wen'an

M_S5.1 earthquake in 2006;(g)Shanxi Datong M_S5.7 earthquake in 1989;(h)Shanxi Datong M_S5.8 earthquake in 1991;

(i)Inner Mongolia Baotou M_S6.4 earthquake in 1996;(j)Hebei Zhangbei M_S6.2 earthquake in 1998;(k)Hebei Zhangbei

M_S5.6 earthquake in 1999;(l)Shanxi Hunyuan M_S5.6 earthquake in 1999

滑动窗长; r则是以实际震中为圆心,10 km为滑动步长,在30～250 km的范围内多尺度半径选取地震样本进行Benioff应变释放曲线的幂律拟和,由此得到m值及m值的拟合标准差Δm。在对累积Benioff应变释放曲线进行非线性拟合时,给定参数的初值,进行迭代计算函数极小点参数值。拟合参数的结果通常与初值的选取有关,因此选取多初值拟合,用拟合的最小残差值确定“拟合最佳解”的办法来提高拟合的精确度(蒋长胜,2004)。在计算结果中为确保m值的可靠性,删除了样本数少于5个、非线性拟合与线性拟合均方差残差比C>1、拟合标准差Δm>0.5的计算结果。因本文仅考虑地震矩加速释放模型(即m<1的情况),图4给出研究时段内华北地区12次中强震前m值在t-r二维空间内克里金插值法绘出的等值线分布图,表1中t-r列给出各次震例前m最低值出现的时—空二维尺度。

由图4可见,12次中强震前均存在或强或弱的AMR现象,但m值空间分布较为复杂,低值出现的时—空尺度差异显著。如1998年河北张北M_S6.2、1999年M_S5.6和2006年河北文安M_S5.1地震前m低值现象明显,均出现m<0.3的时—空分布; 1992年黄海M_S5.3和1991年山西大同M_S5.8地震前出现大范围矩释放的加速与减速交替发生的复杂图像; 1995年山东苍山M_S5.2地震与1995年河北唐山M_S5.0地震前累积应变能几乎为匀速释放状态; 1991年山西忻州M_S5.1、1989年山西大同M_S5.7和1999年山西浑源M_S5.6地震前,距离震中较远的区域出现应变能释放加速,而离震中较近的孕震区则表现出地震矩的匀速释放特征。为提取华北地区具有普遍应用价值的AMR时—空特征尺度,本文将12次中强震按震级大小分为两类,即M_S5.0～5.5和M_S5.6～6.4,同时定义了地震矩加速特征q值为

q=1/n∑ⁿ_i=1[1-m_i(t,r)·(RMS_i-nonlinear(t,r))/(RMS_i-linear(t,r))].(2)

式中,n为震例个数,t和r是计算Benioff应变曲线的时、空参数,m_i(t,r)为i第个震例的幂指数m值在二维空间的分布; RMS_i-linear和RMS_i-nonlinear分别为第i个震例在t-r二维空间中计算样本的线性拟合和幂律拟合的均方根残差,其中RMS_i-linear/RMS_i-nonlinear即为Bowman等(1998)定义的曲率参数,用于判断临界圆“最优搜索半径”的C值。这里的矩加速特征q值实际上是每个震例的特征值q_i归一化后,多震例所反映的地震矩加速释放的平均时、空特征。图5为提取出的华北地区中强震前AMR“优势”时—空特征尺度,即M_S5.0～5.5地震优势时间窗t=4 a、扫描半径r=110 km; M_S5.6～6.4地震优势时间窗t=6 a、扫描半径r=180 km,此参数与前人研究结果较为接近,基本符合最佳搜索半径r与震级的线性相关。

图5 6次MS5.0～5.5地震(a)和6次MS5.6～6.4地震(b)矩加速特征值在t-r二维空间的分布(图中空白色块为地震数目不足或m值拟合标准差过大的位置)
Fig.5 Distribution of AMR characteristic of six MS5.0～5.5(a)and six MS5.6～6.4(b)earthquakes in the t-r two dimensional space(the blank white block indicate the location of earthquake number is not enough or the fitting standard deviation of m-value is too large)

按照地震的临界性概念,一次大地震的孕育初期,孕震系统并非一直处于临界状态,而是随着周围应力状态的不断演化逐渐接近并进入临界

状态。在这一过程中,通常情况下孕震区中小地震活动的频度和强度不断增加,应变能释放加速,直到强震发生,临界状态被破坏(Ben-Zio,Lyakhovsky,2002)。如果主震的震中在临界点附近,主震发生前就可以通过对周围地震的频度或强度在一定时—空范围的扫描寻找未来震中的位置。本文利用上述提取出的华北地区矩加速释放时—

图6 2006年7月4日河北文安MS5.1地震前的矩释放程度空间扫描
(a)震前1年m值分布;(b)震前9个月m值分布;(c)震前6个月m值分布;(d)震前3个月m值分布;
(e)震前1个月m值分布;(f)震时m值分布
Fig.6 Spatial scanning map of seismic moment release before Hebei Wen'an MS5.1 earthquake on Jul.4,2006(a)m-value distribution in 1 year before the earthquake;(b)m-value distribution in 9 months before the earthquake;(c)m-value distribution in 6 months before the earthquake;(d)m-value distribution in 3 months before the earthquake;(e)m-value distribution in 1 month before the earthquake;(f)m-value distribution just before the earthquake

空扫描模型参数对研究时段内12次M_S≥5.0中强震进行时空扫描。首先将研究区划分为0.15°×0.15°的网格空间,网格节点作为未来可能的震中位置,地震可能的“破裂时间”t_f假定为震前1年、9个月、6个月、3个月、1个月和实际发震时刻,分别计算震前矩释放程度m值的空间分布。扫描结果显示,6次M_S5.0～5.5地震前m值空间

分布均杂乱无序,没有明显的应力优势集中区,并且与实际震中和发震时刻之间没有关联性可言。图6为2006年7月4日河北文安M_S5.1地震前AMR现象的扫描图像,震中附近没有出现明显的矩释放加速区域,而是随着实际发震时刻的逐渐逼近,AMR现象逐渐远离实际震中,外围远场的AMR演化图像也没有任何与实际震例时间相依的

图7 1998年1月10日河北张北MS6.2地震前的矩释放程度空间扫描
(a)震前1年m值分布;(b)震前9个月m值分布;(c)震前6个月m值分布;(d)震前3个月m值分布;
(e)震前1个月m值分布;(f)震时m值分布
Fig.7 Spatial scanning map of seismic moment release before Hebei Zhangbei MS6.2 earthquake on Jun.10,1998 (a)m-value distribution in 1 year before the earthquake;(b)m-value distribution in 9 months before the earthquake;(c)m-value distribution in 6 months before the earthquake;(d)m-value distribution in 3 months before the earthquake;(e)m-value distribution in 1 month before the earthquake;(f)m-value distribution just before the earthquake

图8 震前孕震区的Benioff应变释放曲线
(a)1991年山西忻州MS5.1地震;(b)1991年河北唐山MS5.1地震;(c)1992年黄海MS5.3地震;(d)1995年山东
苍山MS5.2地震;(e)1995年河北唐山MS5.0地震;(f)2006年河北文安MS5.1地震;(g)1989年山西大同
MS5.7地震;(h)1991年山西大同MS5.8地震;(i)1996年内蒙古包头MS6.4地震;(j)1998年河北
张北MS6.2地震;(k)1999年河北张北MS5.6地震;(l)1999年山西浑源MS5.6地震
Fig.8 Benioff strain release curves in seismogenic zone before the earthquake (a)Shanxi Xizhou MS5.1 earthquake in 1991;(b)Hebei Tangshan MS5.1 earthquake in 1991;(c)Huanghai MS5.3 earthquake in 1992;(d)Shandong Cangshan MS5.2 earthquake in 1995;(e)Hebei Tangshan MS5.0 earthquake in 1995;(f)Hebei Wen'an MS5.1 earthquake in 2006;(g)Shanxi Datong MS5.7 earthquake in 1989;(h)Shanxi Datong MS5.8 earthquake in 1991; (i)Inner Mongolia Baotou MS6.4 earthquake in 1996;(j)Hebei Zhangbei MS6.2 earthquake in 1998;(k)Hebei Zhangbei MS5.6 earthquake in 1999;(l)Shanxi Hunyuan MS5.6 earthquake in 1999

规律可言。对于6次M_S5.6～6.4震例的AMR扫描结果则理想很多,6次中强震前都出现较长时间的AMR现象,且空间分布较为集中与实际震中有很好的对应关系。以1998年1月10日河北张北M_S6.2地震为例(图7),震前1年至6个月,在实际震中附近出现明显、稳定的矩加速释放区,随着周围应力场的调整,震前3个月至发震,震中附近的异常区域逐渐收缩,m<0.5的低值区域从不同方位向震中区逐渐逼近,时间上也表现出逐渐趋近于临界状态的趋势,这种时—空上的一致性说明地震的孕育过程与AMR现象在物理上有一定的相关性。与此同时在沂沭带南段也逐渐形成一个M<1的矩加速释放区,但在之后的时间里,沂沭带南段及周围地区并未发生M_S≥5地震,也说明了AMR现象及“破裂时间分析”等方法在实际地震预测中的不确定性和复杂性。为进一步考察华北地区中强震前震中区附近的矩加速释放情况,分别以实际震中为圆心,采用与时—空扫描相同的模型参数(即M_S5.0～5.5的地震半径r=110 km,时间窗t=4 a; M_S5.6～6.4的地震半径r=180 km,时间窗t=6 a),计算了12次M_S≥5.0 地震震中区附近的Benioff累积应变能释放情况,相应的Benioff应变释放线性和非线性拟合情况(图8)。从6次M_S5.0～5.5地震的Benioff应变释放曲线来看,除1992年黄海M_S5.3地震外均表现为近匀速释放状态,与对m值的时—空扫描图像 结果基本一致,即震前在孕震区没有明显的应变能加速释放或地震频次增加的现象。而6次

M_S5.6～6.4中强震虽然在时—空扫描中存在大面积、稳定的AMR现象,但对于破裂成核点附近的Benioff应变能释曲线仅4例地震存在矩加速释放现象,这说明AMR现象仅是中强震孕育发生的模式之一。从表1中可以清晰的看出,12次震例中仅5例震前在破裂成核点附近存在加速释放特征,3例震前表现出应变的减速释放特征,另外4例为匀速释放。值得一提的是,在M_S5.0～5.5震例中仅1例在震前表现出应变能加速,占总数的17%; 而M_S5.6～6.4的震例中有4例在震前存在矩加速释放特征,占总数的67%。因此认为华北地区M_S≥5.5中强震前的AMR现象更具有一定的普遍性,但却不是充要条件。

中小地震丛集现象往往被认为是强震发生的前兆或构造活动的外释信号,成为了地震预测的重要依据之一。2012年以来,华北地区最值得关注的震情活动是环渤海地区的中等地震及震群事件的频发,先后发生了辽宁灯塔M5.1、内蒙通辽M5.3、吉林前郭M5.8和河北滦县M4.8、山东莱州M4.6、乳山M4.3等中等地震,以及山东长岛、莱州、乳山、辽宁海城、盖州和河北赞皇、蔚县等一系列围绕渤海及郯庐断裂带北段的典型事件,这足以说明环渤海地区的应力背景处于高值水平。本文应用提取出的华北地区矩加速释放模型参数,采用删除余震并剔除干扰事件(M≥5.0)后M_C2.5的地震目录,假定发震时刻为2014年1月 26日,来检验目前环渤海地区的AMR现象及Benioff应变释放曲线情况,结果如图9所示。

图9 环渤海地区AMR空间扫描图(a)和Benioff应变释放曲线(b)
Fig.9 AMR spatial scanning map(a)and Benioff strain release curves(b)in the Bohai rim region

图9a中显示,在环渤海及黄海地区集中出现了北北西向展布的大面积的低值区域,矩加速释放特征显著。Benioff应变释放曲线表明自2012年开始该区域处于应变能加速累积的阶段,曲率参数C值为0.62,低于华北地区1989年以来12次中强震前的任何一次(图9b)。

本文利用地震矩释放程度m值对1989年以来华北地区12次中强震前的AMR现象进行了回顾性震例检验。(1)通过对值在t-r二维空间的图像分析,提取出了研究区与AMR现象有关的特征时—空尺度,即M_S5.0～5.5地震在半径r=110 km,时间窗长t=4 a,M_S5.6～6.4 地震在半径r=180 km,时间窗长t=6 a的情况下,AMR现象更具有“优势”和普遍性。(2)利用得到的特征时—空尺度对每个震例进行震前1 a至震时,依次递减3个月在空间范围的扫描。结果显示,除1992年黄海M_S5.3地震外的5次M_S5.0～5.5地震前AMR现象没有形成明显的优势集中区,与实际震中和发震时刻没有关联性; 6次M_S5.6～6.4地震前在震中区附近均出现显著、稳定的距加速释放区域,值随着实际发震时刻的临近而逐渐减小或向震中位置逐渐收缩,AMR现象出现的空间区域与实际震中也有较好的对应关系。(3)考察了震前震中区附近的Benioff应变释放情况,12次震例中有5例在破裂成核点附近存在AMR现象,与空间扫描结果基本一致,说明华北地区M_S≥5.5中强震前AMR现象存在一定的普遍性,与区域地震孕育的物理过程有一定的相关性。(4)对研究区目前地震矩释放程度时—空扫描的结果显示,环渤海地区存在显著的矩加速释放特征,Benioff应变释放曲线表明背景应力水平较高。(5)通过上述研究,不难发现地震能量的累积释放与时间的关系比简单的用AMR模型来描述要复杂很多,一些地震前没有AMR现象出现,出现AMR现象也不一定会有中强地震发生,因此用矩释放程度的时—空扫描方法对未来中强震发震地点的预测需要与其它地震学预测方法相互映证,综合考虑。

本研究在成文过程中得到山东省地震局刘希强研究员、山长仑高级工程师的帮助,审稿专家提出了诸多有益建议,在此深表感谢。