收稿日期:2015-04-16
基金项目:国家科技支撑项目专题(2012BAK19B04-01)资助.
基金项目:国家科技支撑项目专题(2012BAK19B04-01)资助.
郯庐断裂带中南段前兆综合信息量地震预测方法研究*
(1.山东省地震局,山东 济南 250014; 2.河北省地震局,河北 石家庄 050021)
Study on Prediction Method of Earthquake Precursor Synthetic Information along the Mid-south Segment of Tanlu Fault Zone
GENG Jie1,WANG Mei1,SHAN Lian-jun2
(1. Earthquake Administration of Shandong Province,Jinan 250014,Shandong,China)(2. Earthquake Administration of Hebei Province,Shijiazhuang 050021,Hebei,China)
Tanlu Fault Zone; earthquake prediction; precursor synthetic information; temporal and spatial analysis; wavelet transform; EMD decomposition technique
备注
摘要
全文
图/表
参考文献
结合地震前兆信息量综合数学表达式和小波变换、EMD分解技术,遴选郯庐断裂带中南段及其附近地区67个台点91个测项的连续性、可靠性较高前兆观测资料,多频域剖析前兆观测资料的变化特征、多层次挖掘地震前兆观测异常、最大程度地展示前兆群体映震效能。结果 显示,地震综合前兆信息量时序曲线,在ML≥5.0地震前会出现不同程度的高值异常(S时序≥0.1),一般高值异常出现1～6个月内,研究区有发生ML≥5.0地震的可能; 地震前兆信息量综合空间扫描的结果显示,未来的震中区多位于最早出现S窗≥0.2异常区的边缘区域,但异常区域的迁移规律性特征不明显,这可能与资料分析时段研究区内发生地震的震级较小及前兆台点的空间分布不匀有一定的关系。
Combined with the comprehensive mathematical expressions of earthquake precursory information,the wavelet and EMD decomposition technique,we selected the precursory observation data with good continuity and high reliability of 91 items at 67 stations in the mid-south segment of Tanlu Fault Zone and its adjacent area in order to analyzing the variation characteristics of precursory observation data in multi frequency domain,multi level digging into objective earthquake precursory anomaly and showed reflecting earthquake ability of group precursor. The results showed that there appeared the high value anomaly(Stime≥0.10)in different level before the ML≥5.0 earthquakes on the time series curve of earthquake synthetic information,and might occur ML≥5.0 earthquakes after high value anomaly appeared within 1～6 months in the study area. The spatial scanning results of synthesis precursory information showed that the future epicenter area located in the edge of the Swindow≥0.2 abnormal area where appeared earliest. But the migration regularity of anomaly region was not obvious,which may have a certain relationship with the earthquake magnitude was smaller and spatial distribution of precursory stations was unevenness in the study area in the process of analyzing data.
引言
郯庐断裂带作为中国大陆东部最显著的深大断裂,其未来的震情趋势一直引起我们的高度关注。根据历史地震活动与现阶段中小震活动特点,该区存在发生中强地震的背景。本文研究在对该带发生过的中强以上地震前出现的前兆异常进行综合分析的基础上,较完整地收集了2007年以来与研究区相关的省份的数字化前兆观测资料及相关台站和观测信息(干扰识别)资料,采用时频分析方法对研究区前兆观测资料进行了异常提取,“边研究边应用”地开展了研究区前兆异常综合信息量提取方法及异常时间和空间综合演化特征分析。
1 技术思路和技术途径
地震预报的研究与实践表明(平建军等,1999,2003a; 王炜等,1996),单项前兆手段的异常时空分布是极其复杂的。没有一种单项前兆手段的异常在所有地震前都出现,亦没有在任何一种单项前兆手段的异常变化出现后都有相应的地震发生,即单项前兆手段与地震的关系具有不确定性。尽管每一个前兆异常都不能反映地震孕育发生的全过程,但它却可能反映地震孕育过程的某些信息,这些信息组合起来在一定的时间、空间内构成一个综合的地震前兆信息场,这个信息场的时、空演变将有可能反映地壳介质破裂的震源物理过程。因此,深入开展地震前兆场的综合信息量及其时、空动态演变研究,对于探寻地震孕育发生的普遍规律,提取地震预报的判据和指标,提高地震预报水平具有至关重要的意义。
目前我国的地震监测有地形变、地磁、地电、水位、水温、水化学、地应力等多种学科,各学科甚至同一学科不同观测方法所测的物理量都不相同,要想将这些量纲各异的各种地震监测资料有机地组合起来,进行地震预报综合分析,首先要设法将每一种地震观测异常统一转化为同一信息量值,罗兰格(2002)对此进行了研究,提出了地震异常前兆信息量综合数学表达式法。经对实际地震观测资料进行计算,认为该表达式不仅可计算提取各种地震观测资科异常持续变化的前兆信息量,还可定量地描述异常结束后其前兆信息的延续性,从而使其地震前兆信息更加突出合理。这就为我们进一步对比研究地震异常群体前兆信息特征、开展地震前兆综合信息量的地震预报奠定了基础。
前兆观测资料中,往往存在着趋势上升、下降、周期等非平稳变化,这就需要选择平滑滤波、傅里叶变换、变化率等数学方法,事先进行消周期、消趋势项的数据预处理工作。地震前兆信息可能存在于不同的频带中,既可能存在于高频段,也可能存在于低频段,还可能兼而有之。在现行的数据预处理工作中,对前兆观测数据或去低频,或消高频,通常会导致遗漏一些有用的异常信息,使异常信息不能完全表征出来。近几年越来越被重视且得到广泛应用的小波变换和经验模态(EMD)分解技术(程正兴,武铁敦,1994; 武安绪等,2006),为解决这一问题带来了便利。通过对原始前兆观测资料进行小波变换或EMD分解,可有效地把各种频率成份从中分离出来,且分离出来的时频序列资料数据变化平稳,非常便于应用地震前兆信息量综合数学表达式对其进行无量纲化震兆信息提取。
本文结合地震前兆信息量综合数学表达式和小波变换、EMD分解技术,分析处理了研究区及其附近地区67个台点、91个测项的前兆观测资料,以期为多频域剖析前兆观测资料的变化特征、客观自然地多层次挖掘地震异常、最大程度地展示前兆群体异常在前兆综合信息量时序和空间上的演化特征,为有效的进行震情预测提供可信度较高的依据。
2 地震综合前兆信息量提取方法
3 研究区地震前兆综合信息量的计算及其映震效能分析
4 分析结论
(1)小波和EMD分解技术的应用,可有效地把各种频率成份从原始前兆观测资料中分离出来,多频率域的充分展示前兆观测资料的异常信息,增加了对地震异常的辨识和挖掘力度。将地震前兆信息量综合数学表达式和小波、EMD分解技术结合在一起,分析处理前兆观测资料,为多频域剖析前兆观测资料的变化特征,客观自然地多层次挖掘地震异常、最大程度地展示前兆群体映震效能提供了新的途径。
(2)分析研究区地震综合前兆信息量时序曲线,在ML≥5.0地震前会出现不同程度的高值异常(S时序≥0.10),一般高值异常出现后1～6个月内研究区有发生ML≥5.0地震的可能。
(3)对研究区地震前兆信息量综合空间扫描的结果显示,未来的震中区多位于最早出现S窗≥0.2异常区的边缘区域,但异常区域的迁移规律性特征不明显,这可能与资料分析时段研究区内发生地震的震级较小及前兆台点的空间分布不匀有一定的关系。前兆台点的空间分布不匀某种程度上导致前兆异常分布不匀,也会导致等值线连线在一定程度上存在失真现象。
3.1 研究区前兆观测资料基本情况本文研究区域涉及辽宁、河北、北京、天津、山东、河南、江苏、安徽8个省(市)(图1),在对上述区域各台站测项观测资料进行充分收集、分析和调研基础上,结合2012年度全国前兆观测资料预效能评估结果,系统收集处理了各省(市、区)预报效能较好、干扰因素相对较少、观测数据可信度较高且资料完整性相对较好的67个观测台(点)合计91个测项2007-01-01～2014-06-10资料进行分析计算。
分析计算所使用原始数据来源于上述省份及国家台网中心数字化前兆台网Oracle原始和预处
理数据库,在此基础上生成整点值或日值数据,作为计算的基础数据。原始观测数据中,由于人为、仪器原因(标定、校测和故障)等因素,经常会产生观测值的非正常变化,直接影响数据的分析,对于这些非正常变化、非客观的资料动态,不同形态特征要进行具体分析,针对经常出现的非正常变化情况,在计算前分别对其进行了预处理,本次研究的主要预处理方法有:
图1 研究区前兆信息量分析台(点)分布图
Fig.1 Distribution map of the stations(points)for precursory information analysis in the study area缺值:由于仪器原因,在某些时间段内缺值,数据中以“999999”表示。对于缺值这种情况,依据一定的方法进行了插值。常用的插值法有线性插值法、抛物线插值法、拉格朗日插值法和多元线性回归拟合外推插值法等。
突降:数据中某些点由于仪器等干扰原因造成突降。这种情况采取了根据突升前后值,将突降后数据接回突降以前的处理。
突跳:由于已知的干扰方面原因,数据出现几个点高于或低于周围点的突跳变化情况,根据突跳点两端数据的平均值,插入突跳处,以消除突跳点。
3.2 单项前兆信息量的计算与分析对于研究区域每一前兆观测资料(形变、电磁、流体等),首先判别其是否存在非平稳变化。若存在非平稳变化,则使用小波变换或 EMD分解技术,将其分解为一组10阶的时频序列,分析各时频序列映震能力,并选择出映震效能最优的时频序列,应用地震前兆信息量综合数学表达式,提取其无量纲前兆信息量。
对研究地区前兆观测资料进行小波或经验模态(EMD)分解,结果显示该方法对非平稳观测数据的异常信息具有较强的挖掘和识别能力,部分原始数据曲线上异常特征不明显的测项,经识别后,显示有明显的异常信息存在。
3.3 地震前兆信息量综合时序和空间分析对研究区域前兆观测资料,统一按上述方法提取前兆信息量后,即可开展时、空综合分析。
3.3.1 地震前兆信息量综合时序2007-01-01～2014-06-10期间,研究区及其附近地区发生ML≥5.0地震4次,显著震群事件1次,我们对上述地震进行地震前信息量综合时序回溯性检验。对经过小波变换或EMD分解后的单项前兆异常信息量提取结果按(5)式进行综合信息量综合时序分析,结果见表1、图2。
除安徽安庆地震前信息量综合时序曲线未出现异常外,其余4次地震前均显示出现明显的高值异常出现,异常出现在震前1～105 d,异常幅度最大为0.158。上述研究区内发生的5次地震均发生在前兆综合信息量时序曲线的高值异常出现后,由于研究时段内所发生的地震震级差异不大,在前兆综合信息量时序曲线上显示,异常幅度和异常持续时间差异不大。
表1 研究区及其附近地区ML≥5.0地震前兆综合信息量时序异常特征(S时序≥0.1)(2007-01-01～2014-06-10)
Tab.1 The time series anomaly characteristics(Stime≥0.1)of earthquake precursor information with ML≥5.0 earthquakes in the study area and its vicinity from Jan.1,2007 to June.10,20141995年9月20日山东苍山MS5.2地震是发生在郯庐断裂带上的显著地震事件之一,苍山地震前兆综合信息量时序值计算结果显示(图3),1995年8月10日和8月20日分别出现S时序为0.28和0.23的高值异常,震后时序曲线迅速恢复到正常背景值波动。
图2 研究区前兆综合信息量时序曲线(2012-01～2014-09)
Fig.2 The time sequence curve of earthquake precursor comprehensive information in the study area from Jan.2012 to Sep.20142015年度山东省地震趋势会商分析,2012-01-01～2014-09-30对山东地区前兆测项观测资料进行时序前兆综合信息量计算结果显示,2013年11月23日莱州ML5.0地震前,于2013-08-10～2013-09-10出现S时序≥0.2的高值异常时段,并于2013年9月10日出现最高值达0.301; 2014-01-01～2014-09-30 S时序计算结果,仅在
2014年7月31日出现0.216的异常值,8月恢复至正常波动,再无高值异常显示,至2015年3月,山东地区无ML≥5.0地震发生。
图3 1995年苍山MS5.2地震前后前兆综合信息量时序曲线
Fig.3 The time sequence curve of earthquake precursor comprehensive information Before and after the Cangshan MS5.2 earthquake in 19953.3.2 地震前兆信息量综合空间分析对2007-01-01～2014-06-10时段研究区前兆测项资料经过小波或EMD分解后的单项前兆异常信息量提取结果按式(6)进行综合信息量空间分析,结果见表2。
对前述5次地震进行地震前信息量综合空间回溯性检验表明,除安徽安庆地震前信息量综合空间扫描未出现S窗≥0.2的异常区外,其余4次地震前均显示有明显的异常区出现(图4)。
对2007-01-01～2014-06-10时段研究区综合信息量空间扫描结果,2007-01-01～2014-05-31无S窗≥0.2的异常区出现。研究区在2012-06-10出现S窗≥0.2的异常区后,2012-06-30在该异常区西侧又出现一个S窗≥0.2的异常区,之后相邻的两个异常区逐渐扩展,至2012-07-10达最大,之后逐渐收缩。2012年7月20日高邮ML5.3地震发生在最初出现S窗≥0.2异常区边缘,S窗最大值0.98。
表2 研究区及其附近地区ML≥5.0地震前综合信息量(S窗≥0.2)空间演化特征(2007-01-01～2014-06-10)
Tab.2 Spatial evolution characteristics of earthquake precursor comprehensive information(Stime≥0.20) before ML≥5.0 earthquakes in the study area and its vicinity from Jan 1. 2007 to Jun.10,2014
图4 2013年11月23日山东莱州地震前地震前兆综合信息量空间演化特征(a)2013-06-30;(b)2013-07-10;(c)2013-07-31;(c)2013-08-10;(d)2013-09-10
Fig.4 The spatial evolution map of earthquake precursor comprehensive information before Shandong Laizhou ML5.0 earthquake on Jan. 1,2013 in Shandong Province2013年1月23日灯塔ML5.1地震的回溯性检验显示:2012-11-10研究区开始出现S窗≥0.2的异常区,异常区逐渐扩大,至2012-12-10达到最大后逐渐收缩,至2013-01-20后异常区域消失。2013年1月23日灯塔ML5.1地震发生在最初出现S窗≥0.2异常区边缘,S窗最大值0.999 6。
2013-01-20～2013-06-20期间,研究区内无S窗≥0.2异常区出现。2013-06-30在沂沭断裂带南段西侧首先出现S窗≥0.2的异常区,该异常区逐渐扩展,至2013-07-10达最大,并且在莱州ML5.0地震区域附近出现异常区,2013-07-31后莱州地震附近区域的异常区逐渐收缩,至2013-08-10完全消失,与此同时,前述沂沭断裂带南段西侧的异常区也经历了逐渐缩小的过程。2013年11月23日莱州ML5.0地震发生在最初出现S窗≥0.2异常区边缘,S窗最大值0.997(图4)。
2014年4月1日发生乳山ML4.6震群,该震群是2014年度山东地区最为显著的地震事件,在华北地区也十分突出。研究区空间综合信息量扫描结果显示,2013-12-10开始出现S窗≥0.2的异常区,之后异常区逐渐扩展,至2014-01-20达最大后逐渐缩小(2014-02-10)。2014-02-20在异常区西侧出现异常区,2014-03-10乳山震群附近的异常区消失。乳山震群发生在最初出现S窗≥0.2异常区边缘,S窗最大值0.997。
- 程正兴,武铁敦.1994.小波的发展与应用[J].微机发展,(5):8-10.
- 贾炯,平建军,刁桂苓,等.2010.华东地区地震综合前兆信息量及其映震能力分析[J].华北地震科学,28(2):48-52.
- 李广鑫,程式,李振兴.1992.综合积分预报方法在四川地区的应用研究[J].四川地震,(3):41-46.
- 罗兰格.2002.我国地震综合预报方法研究的回顾与展望[J].华北地震科学,20(4):1-l8.
- 平建军,王秀英.1999.地震前兆信息量的研究(二)地震前兆信息量的综合数学表达式[J].华北地震科学,17(3):8-l4.
- 平建军,张永仙,单连君,等.2003a.华北地区地震短期异常分析及其前兆信息量的提取与计算[J].华北地震科学,21(3):1-7.
- 平建军,张永仙,单连君,等.2013b.震灾防御技术,8(4):397-407.
- 平建军,张永仙,单连君.2013c.华北地区地震短期异常分析及其前兆信息量的提取与计算[J].华北地震科学,21(3):1-7.
- 王炜,吴耿锋,黄冰树,等.1996.基于模糊神经网络和符号的地震预报专家系统NGESEP[J].中国地震,12(4):339-346.
- 武安绪,林向东,穆会泳,等.2006.EMD新技术在数字波形预处理中的初步应用[J].华南地震,26(1):l33-138.
- Morlet J. G.,Arens G.,Fourgeau E..1982. Wave propagation and sampling theory:complex signal and scattering in muti 2 layered media[J].Geophysics,47(2):203-211.
对工作区域每一前兆观测资料(地壳形变、地电、地磁、水化、水位、水温、应力、应变等)首先判别其是否存在非平稳变化。若存在,则使用小波变换或EMD分解技术,将其分解为一组多阶时频序列,分析各时频序列映震能力,并挑出映震效能最优的时频序列,应用地震前兆信息量综合数学表达式,提取其无量纲前兆信息量; 如不存在,则可根据地震前兆信息量综合数学表达式,直接提取计算其无量纲前兆信息量。地震前兆信息量综合数学表达式及其式中各参数的含义,请参见有关文献(平建军等,1999,2013b,c; 贾炯等2010; 李广鑫等,1992),这里不再赘述。
2.1 单项前兆信息量的提取数字化前兆资料多为典型的观测信号,为非平稳过程,具有不稳定性、变化快等时频特点。在多种天然与人为因素(不包括观测因素)影响下,实际前兆观测资料往往存在趋势和多种周期成份非平稳的动态变化。其中某些周期性变化,可与具有一定周期性的构造活动以及地震孕育发生过程产生力学耦合作用,这类周期性变化,即为可能的地震异常。这就是说,前兆数据序列中仅是某些周期成分变化才与地震的孕育和发生有关。小波变换和EMD分解技术的特点,正好适合从多种周期变化中提取出与地震有关的周期性异常信息。
2.1.1 小波变换基础对于离散序列信号f(x), 在小波函数ψ(t)∈L2(R)中,尺度因子(伸缩因子)a和平移因子b(a,b∈R),也需要离散化,应用离散小波变换(DWT)作为不同频率的信息识别基础,即
DWTa,b=∫Rf(t)ψm,n(t)dt.(1)
在计算中,采用a=2k。随着k的增加,信号从最高频向低频分解。当k=0时,信号为采样频率,k=1时,将频率二等分,依此类推。
对于数字信号可以近似地表示为
f(x)Aj f(x)=ajk f(x)+djk f(x).(2)
其中,ajk f(x)与djk f(x)分别为信号f(x)在分辨率为j时的近似部分与分解(或细节)部分,正交展开系数ajk与djk分别为离散近似与离散细节。Mallat算法就是ajk与djk分别满足
{ajk=∑∞i=-∞h^-n-iaj-i,
djk=∑∞i=-∞g^-n-iaj-i.(3)
而Mallat的重构算法为
aj-1k=∑∞hn-2iaji+∑∞gn-2iaji.(4)
其中,hn与gn分别是尺度函数和小波函数的滤波器,并且h^-n=h-n,g^-n=g-n。
2.1.2 小波选择与频率分析在地球物理领域,Morlet小波应用比较广泛。因此,我们在分处理地震活动与前兆信号时,均采用Morlet小波(Morlet et al.,1882)。常用的是复值Morlet小波。
小波变换系数aik与djk实际上是窗口小波变换。当j固定时,随着k的变化,aik与dik都占满了整个时间域,而且没有重叠。随着j的变化,占满了整个频率域,也没有重叠。设信号f(x)在时间域[0,T]中采样,共采了N个点,采样间隔为Δt=T/N。{fi}(i=0,1,…,N-1)的频带为[0,Ω],求出近似部分的系数{a0k}(k=0,1,…,N-1),其频带为[0,Ω],而ajk与djk(j=1,2,…,l=0,1,…,2j-1)的频带如何呢?a0k的频带为[0,Ω],a1k的频带为[0,Ω/2],d1k的频率为[Ω/2,Ω]。以此类推,ajk与djk的频带分别为[0,Ω/2j]与[Ω/2j,Ω/2j-1]。
2.1.3 不同频率的信息识别思路由于数字化前兆的观测精度提高,在干扰因素排除后,对趋势异常与短期异常的识别与排除也是一大问题,而通过小波变换方法,对不同频率范围内的信息(或高频与低频信息)进行识别与分解。根据式(2)以及小波分解的近似部分aik f(x)、细节部分djk f(x)与频率的关系,对观测资料进行近似部分(低频)与细节部分(高频)信息分离。根据观测量的物理含义(如水动态观测资料中反映固体潮汐的日波、半日波、半月波、月波等信息),在分析干扰因素的基础上,从近似部分中确定出趋势变化信息; 从细节部分识别出短期异常。
2.1.4 小波基的选取前兆模拟观测资料的观测周期主要为日观测和整点值观测,数字化观测采样周期主要为分钟值观测和整点值观测。观测资料中除含有地壳应力应变的信息外,还含有诸如降雨、温度、开采和仪器故障等不同频率的干扰成分,但大多数资料都呈现出规律性较好的年变化形态,因此,在小波分析时,选用正则性较好的离散正交小波(Daubechies小波)进行分析。采用db4小波对数据进行处理,首先考虑到db4小波在地球物理学研究范围被广泛采用; 其次考虑到db4小波是紧支撑正交小波,可以使Mallat算法更快捷。它的光滑性也可以更高精度地模拟和分析信号。db4小波随着分解层级的增加其正则性也增加,它抑制了该多项式信号在零阶和一阶的部分信号,而仅对该信号的二阶部分及噪声进行分解。因此,分解的细节信号部分db4中包含了噪声信号的不规则性,其余各层细节中的信号周期性随层级的增加而增大。另外也考虑到db4小波在时域和频域局部化方面的强劲性。
本研究对于db4小波分析,将经预处理后的前兆观测值分解成10个细节(10阶)时频序列,多频域分析前兆观测资料的变化特征和各时频序列映震能力,客观自然地多层次挖掘地震异常,在此基础上选择最优时频序列,应用地震前兆信息量综合数学表达式,提取其无量纲前兆信息量,以最大程度地展示前兆群体异常映震效能,开展震情预测。
对研究区域前兆观测资料,统一按上述方法,提取单项前兆信息量后,即可开展如下时间和空间综合分析。
2.2 地震前兆信息量综合分析2.2.1 地震前兆信息量综合时序分析由于实际收集到的研究区域数字前兆观测资料的起止时间有长有短,为抑制不同时段因存在的数字前兆观测资料数目的多少而对综合结果所产生的影响,在对研究区域进行地震综合前兆信息量时序分析时,根据每一数字前兆观测资料地震信息量的计算结果,逐时刻地截取并计算研究区域地震前兆综合信息量时序值:
S时序=1/N∑Ni=1Si.(5)
式中,S时序为某时刻工作区域地震综合前兆信息量时序值; Si为该时刻某前兆观测资料的前兆信息量值; N为该时刻所有前兆观测资料总数。
2.2.2 地震前兆信息量综合空间分析根据研究区域每一数字前兆观测资料地震信息量的计算结果及其台站的空间经、纬度,再考虑同一观测台站不同观测项目地震信息量间的差异性,以及前兆观测台站在空间分布上的不均匀性对地震信息场所产生的影响,按式(6)逐时刻地截取并绘制研究区域地震综合前兆信息量空间分布图,通过研究中强震前区域震兆信息空间演化特征,从而建立其中强震空间发震地点判据指标。
S窗=1/N∑Ni=1Si.(6)
式中,S窗为某时刻工作区域地震综合前兆信息量时序值; Si为该时刻某前兆观测资料的前兆信息量值; N为该时刻该经纬度节点扫描范围内前兆观测资料项目总数。
