2. 中国地震局第二监测中心, 西安 710054;
3. 河北省地震局流动测量队, 河北保定 071000
2. The Second Monitoring Center, China Earthquake Administration, Xi'an 710054, China;
3. Mobile Survey Team of Hebei Earthquake Agency, Baoding 071000, Hebei, China
邢台地震后, 中国地震局开始在京津唐等地区尝试开展流动重力监测工作, 对流动重力的发展起到了积极作用(赵炜等, 1981)。1976年7月28日唐山地震发生前, 观测数据记录到的重力变化为98μGal(贾民育等, 1985), 业内人士对这一结果的肯定, 为重力在地震监测预报工作中的作用奠定了良好基础(项爱民等, 2007)。近年来, 随着地震监测预报研究工作的不断深入及观测技术的不断发展(项爱民等, 2007;刘冬至等, 2007), 以及对地震动态重力监测网监测能力要求的不断提升, 2009年10月中国地震局完成了对大华北地区主要活动构造带上的区域流动重力复测网的调整、优化和改造, 新增绝对重力观测点的同时, 对测网间有效连接进行了改进, 改进后的大华北地震动态重力监测网对MS≥5.0地震有了较好的监测能力(李辉等, 2010)。
重力观测数据受地球深部、浅部及其外部空间所有物质源的密度变化及位移等综合因素的影响, 近年来,利用小波多尺度分解方法特有的低阶小波细节尺度不变属性, 可以很好地将位场数据进行区域异常和局部异常分离, 该方法作为滤波工具被广泛使用在地震数据处理中(宋治平等, 2003;万永革等, 2003;张燕等, 2003)。Fedi等(1998)利用离散小波进行位场分离, 并利用Minimum Entropy Compactness准则来选取区域场;Ucan等(2000)利用小波多尺度分解研究了1个理论模型的异常场分解, 实验结果证明该方法不会对观测数据造成明显的影响, 同时也不会受到区域场与残差场功率谱叠加的显著影响, 取得了较好效果;侯遵泽等(1997)利用小波多尺度分解方法对中国大陆地区布格重力异常数据进行处理, 成功地对中国大陆地壳密度差异进行了反演;高德章等(2000)利用二维小波多尺度分解方法对中国东海及临区自由空间重力异常进行分解, 得到了沉积基底面和莫霍面产生的重力异常;刘少明等(2004)利用小波多尺度分解方法对三峡流动重力资料进行处理, 并分析了网格间距对分解结果的影响以及各阶小波与重力变化波长的关系;刘芳等(2013)利用小波多尺度分解方法对大华北地区流动重力资料进行位场分离;李大虎等(2014)对四川地区流动重力资料进行位场分离与异常特征提取;郭树松等(2014)利用小波多尺度分解方法对青藏高原东北缘地区的相对重力场进行动态变化分析, 并结合震前异常分布情况研究小波多尺度分解技术在地震预测中的应用。
鉴于河北省内及其周边邻近地区地震活动水平较高, 1966年以来发生了1966年宁晋-隆尧7.2级震群、1967年河间6.3级地震、1976年唐山7.8级地震、1998年张北6.2级地震等多次6级以上地震。为更好地加强对河北省震情的跟踪分析, 本文基于小波多尺度分解方法,对研究区范围内2010~2014年期间8期流动重力观测数据进行重力场有效场源分离, 得到不同阶次的小波变换逼近场与小波变换细节场, 并结合期间河北省及其邻区地震活动性进行震例回溯分析, 在总结经验的同时, 能够为后续震情研判提供较好的参考。
1 资料选取及处理方法中国地震局于2009年10月完成了对大华北地区重力区域网的调整、优化和改造, 因此, 本文选取监测网整改完成后的2010~2014年的重力观测数据进行处理、分析。研究区内重力观测点位分布及1976年以来研究区内M≥6地震的震中分布见图 1。
小波分析方法由Mallat于20世纪80年代提出, 近些年小波变换方法发展迅速, 并被广泛应用于地球物理信号处理领域(Mallat, 1989)。通过选取合适的小波基和小波尺度因子参数, 可将位场数据分解成不同尺度空间下的细节场和逼近场。重力小波细节场可抑制深部场响应, 突出局部场信息, 反应区域较为精细的地壳结构特征;小波逼近场则反应地球深部近莫霍面的重力变化情况。本文将采用双正交小波基函数“bior3.5”对各期重力位场数据进行4阶小波多尺度分解。
二维重力场的小波多尺度分解表达式可简写为(刘天佑, 2007、2009)
$ \Delta g(x) = {A_4}G + {D_4}G + {D_3}G + {D_2}G + {D_1}G $ |
式中, A4G为4阶小波变换低频逼近部分, DjG为第j阶小波变换高频细节部分, 其中, j=1, …, 4。
功率谱分析法作为小波多尺度分解的辅助手段, 可以定量地估算不同尺度下场源的近似埋深, 从而定量地说明各不同深度处的重力变化情况。
2 资料处理结果及分析通过对预处理数据进行4阶小波多尺度分解, 以1408-1403期(即2014年8月与2014年3月2期观测数据的相对重力变化, 命名规则下同)观测数据处理结果为例(图 2), 可以看出:受地表等多种因素影响, 1阶小波细节结果变化比较复杂, 相比其他阶次, 其结果比较零散, 等值圈闭合范围较小, 反应地表重力变化情况;而2阶、3阶小波细节结果随阶次的增加, 等值圈闭合范围随之变大, 说明随深度增加, 重力变化范围相对集中;4阶小波细节结果反应近莫霍面距离重力场变化情况, 出现较大的成片区域连在一起, 说明中下地壳深部变化范围更为集中。
为了更好地将小波多尺度分解结果应用到震情跟踪分析的研判中, 且考虑到整合后的大华北重力监测网的监测能力, 本文将结合2010~2014年110°~120°E, 36°~42°N范围内MS4.5以上地震, 总结震前各阶小波细节变化特点, 为后续震情研判提供一定参考。其中, 参与统计的地震基本信息及震中分布, 见表 1、图 3。
图 4为1008-1003期小波细节多尺度分解结果, 通过图 4可以看出,2010年6月5日山西阳曲5.0级地震前2阶、3阶小波尺度分解结果在距震中一定范围内出现正负重力变化异常区, 重力差异变化达50μGal以上, 地震发生在重力差异变化大的带上。
图 5为1103-1008期小波细节多尺度分解结果, 通过图 5可以看出,2011年3月7日五寨4.5级地震前2阶小波尺度分解结果在距震中一定范围内出现50μGal以上的重力差异变化。
图 6、7分别为1403-1308期、1408-1403期4阶小波细节多尺度分解结果, 通过图 6、7可以看出,2014年9月6日涿鹿M4.8地震前, 1403-1308期、1408-1403期2阶小波细节多尺度分解结果在距震中一定范围内出现50μGal以上的重力差异变化。结合表 1可以看出, 2014年9月6日涿鹿M4.8地震震源深度为20km, 小波细节重力场异常变化持续时间较长, 这可能与震源较深有关。
图 8为1103-1003期小波细节多尺度分解结果, 通过图 8可以看出,2010年6月5日阳曲5.0级地震2阶、3阶小波尺度分解结果在距震中一定范围内出现50μGal以上的重力差异变化, 2011年3月7日五寨M4.5地震前2阶、3阶小波尺度分解结果在距震中一定范围内出现50μGal以上的重力差异变化。
图 9为1408-1308期小波细节多尺度分解结果, 通过图 9可以看出,2014年9月6日涿鹿4.8级地震前3阶、4阶小波尺度分解结果在距震中一定范围内出现50μGal以上的重力差异变化。此外, 结合2014年9月6日涿鹿M4.8地震震源深度20km及4阶小波细节重力场异常显著变化这一现象,可以认为小波多尺度分解结果中某一区域出现高阶异常变化时, 未来发震的震源深度可能越大。
小波变换多尺度分析结果在统计意义上等效于不同深度处重力变化情况。根据位场频谱理论, 功率谱斜率与场源埋深成正相关, 因此通过功率谱直线段斜率或横轴切点计算各阶小波细节所对应的平均场源深度(杨文采等, 1978;Cianciara et al,1976;Syberg, 1972;Bhimasankaram, 1977)。本文通过对研究区范围内预处理数据及小波细节多尺度分解结果进行功率谱计算, 得出以下结论:对预处理数据进行功率谱计算, 其结果可以比较清楚地分成4个不同斜率的阶段(图 10中红色线段);对各阶小波变换多尺度分解结果进行功率谱计算, 其结果中均有明显的直线段(图 11中红色线段), 因此对预处理数据进行4阶小波多尺度分解的过程合理。通过计算可以得知,1阶小波细节计算的近似深度约为1km, 反映地表重力变化情况;2阶小波细节计算的近似深度约为11km;3阶小波细节计算的近似深度约为23km, 反映中上地壳重力场变化情况;4阶小波细节计算的近似深度约为36km, 反映地壳深部重力场变化情况。
通过对2010~2014年大华北范围内的重力观测数据进行数据预处理、小波多尺度分解计算及将其与地震活动性相结合进行分析、功率谱计算, 得出如下结论:
(1) 小波多尺度分解方法可以很好地将重力位场数据分解成不同尺度空间下的深部场和浅部场、区域场和局部场;
(2) 功率谱计算作为小波多尺度分解的辅助手段, 可以定量地估算不同尺度下场源的近似埋深;
(3) 小波多尺度分解中各阶细节结果反应不同深度下的重力变化情况, 结合震例进行回溯性分析, 可以发现:当2阶、3阶、4阶小波细节结果出现大于50μGal的重力差异变化时, 在异常变化区域一定范围内(基本以正负异常变化中心为圆心, 半径200km范围内), 半年内会发生4.5级以上地震的可能性较大;
(4) 小波细节变化结果中出现异常的阶次越高、阶数越多时, 往往发生的地震震级越大、震源深度越深;
(5) 在数据处理结果分析过程中, 应将不同时间尺度下的处理结果结合起来进行分析, 其中半年尺度、1年尺度实验结果对1年时间内发生地震的预测效果较好。
综上所述,利用小波多尺度分解方法对重力位场数据进行处理, 可以较好地将处理结果应用于震情研判工作。但由于本文所用数据时间尺度及该时间段内震例有限, 多年时间尺度的重力位场小波多尺度分解结果可能受自身变化幅度及所选震级较小等因素影响, 在震例回溯方面效果不佳, 后续将考虑从更长时间尺度及更大震级方面进行实验、分析。
致谢: 感谢“重力台网中心”提供相关数据以及匿名专家对本论文提出宝贵意见。
高德章、侯遵泽、唐建, 2000, 东海及邻区重力异常多尺度分解, 地球物理学报, 43(6): 842-849. DOI:10.3321/j.issn:0001-5733.2000.06.013 |
郭树松、刘芳、祝意青, 2014, 小波多尺度分解在地震预测中的应用, 大地测量与地球动力学, 34(4): 34-38. |
侯遵泽、杨文采, 1997, 中国重力异常的小波变换与多尺度分析, 地球物理学报, 40(1): 85-95. DOI:10.3321/j.issn:0001-5733.1997.01.010 |
贾民育、刘敬宽、吳兵等, 1985, 论唐山地震前重力变化的可靠性, 地壳形变与地震, 5(3): 249-257. |
李大虎、丁志峰、梁明剑等, 2014, 四川地区流动重力资料的位场分离与异常特征提取, 地震学报, 36(2): 261-274. DOI:10.3969/j.issn.0253-3782.2014.02.011 |
李辉、徐如刚、申重阳, 2010, 大华北地震动态重力监测网分形特征研究, 大地测量与地球动力学, 30(5): 15-18. |
刘冬至、邢乐林、徐如刚等, 2007, FG5/232绝对重力仪的实验观测结果, 大地测量与地球动力学, 27(2): 114-118. |
刘芳、祝意青、陈石, 2013, 华北时变重力场离散小波多尺度分解, 中国地震, 29(1): 124-131. DOI:10.3969/j.issn.1001-4683.2013.01.014 |
刘少明、申重阳、孙少安等, 2004, 小波多尺度分解特征分析, 大地测量与地球动力学, 24(2): 34-41. |
刘天佑, 2007, 位场勘探数据处理新方法, 北京: 科学出版社.
|
刘天佑, 2009, 重磁勘探软件手册(GMS3.0)方法原理, 武汉: 中国地质大学.
|
宋治平、武安绪、王梅等, 2003, 小波分析方法在形变数字化资料处理中的应用, 大地测量与地球动力学, 23(4): 21-27. |
万永革、齐福荣、孟晓春等, 2003, 中国大陆及华北地区地震资料的小波分析, 大地测量与地球动力学, 23(4): 28-33. |
项爱民、孙少安、李辉, 2007, 流动重力运行状态及质量评价, 大地测量与地球动力学, 27(6): 109-114. |
杨文采、郭爱缨、谢玉清等, 1978, 重磁异常频率域解释的理论与方法(下), 物化探研究报导, 78(2): 134-178. |
张燕、吴云、刘永启等, 2003, 潮汐形变资料中地震前兆信息的识别与提取, 大地测量与地球动力学, 23(4): 34-39. |
赵炜、马丽, 1981, 重力重复测量精度讨论, 地壳形变与地震, (3): 23-34. |
Bhimasankaram V L S, 1977, Interpretation of gravity anomalies due to finite inclined dikes using Fourier transformation, Geophysics, 42(1): 51-59. DOI:10.1190/1.1440713 |
Cianciara B, Marcak H, 1976, Interpretation of gravitv anomalies by means of local power spectra, Geophys Prospect, 24(2): 273-286. DOI:10.1111/j.1365-2478.1976.tb00925.x |
Fedi M, Quarta T, 1998, Wavelet analysis for the regional-residual and local separation of potential field anomalies, Geophys Prospect, 46(5): 507-525. DOI:10.1046/j.1365-2478.1998.00105.x |
Mallat S G, 1989, A theory for multiresolution signal decomposition:the wavelet representation, IEEE Trans Pattern Anal Mach Intell, 11(7): 674-693. DOI:10.1109/34.192463 |
Syberg F J R, 1972, A Fourier method for the regional-residual problem of potential fields, Geophys Prospect, 20(1): 47-75. DOI:10.1111/j.1365-2478.1972.tb00619.x |
Ucan O N, Seker S, Albora A M, et al, 2000, Separation of magnetic fields in geophysical studies using a 2-D multi-resolution Wavelet analysis approach, J Balk Geophys Soc, 3(3): 53-58. |