0 引言

大量观测研究表明地震孕育和发生过程中伴随着地磁场的变化(徐文耀, 2003)。但是, 地球内部产生的地球基本磁场对地磁场的贡献远远大于震磁效应产生的变化, 而与地震孕育密切相关的磁场变化很小, 以至于被湮没在仪器记录到的数据中, 所以传统的时间序列曲线很难辨别地磁异常。因此, 在较强的空间正常变化磁场的背景下识别和提取与震源破裂相关的磁场微弱变化特征是努力的方向(冯志生等, 2010)。在地磁台站记录的资料中含有大量的非稳态成分, 包括数据的偏移、趋势上升或下降、突变(Dunlop et al，1994)。这些信号成分反映了地磁信号的重要特征, 并且它们的频域特性又都是随时间而变化的。而地磁场是综合叠加构成的场, 其中包含了丰富的频谱成分(詹志佳等, 1992), 用数字信号处理的方法对地磁数据进行时频分析, 有可能提取震磁前兆异常(陈伯舫等, 2006；范莹莹等, 2010)。

数字信号处理中最基本和最常用的方法是傅里叶分析, 但是它只能提供全局频率特性而缺乏局部信息, 而这对于非平稳信号至关重要(马世伟, 2000)。针对此类问题, 发展出了时频分析方法, 该方法能够在时频联合域的同时分析信号的时间和频率特性。平滑伪魏格纳-维勒分布(SPWVD)是其中的一种, 该方法能够在抑制交叉项的同时, 又能在信噪比较低的情况下准确分辨出微弱信号。何康等(2013、2014)率先将魏格纳-维勒分布(WVD)时频分析方法引入到地磁数据的研究与震例分析中，其使用平滑伪魏格纳-维勒分布(SPWVD)对2008~2013年川滇地区M_S≥5.5地震(共10个)前3个月内各地磁台站4.4h分量的振幅进行研究和统计, 发现在6次地震前多台4.4h分量出现同步高值变化。张明东等(2018)使用相同方法对中国大陆的一些显著地震进行分析, 发现一些地震在震前3个月内, SPWVD时频分析地磁结果存在高值异常且可能为震前一年的极值, 另外, 还发现高值频次在震前显著增大的异常。

以上学者主要研究了2016年以前的显著地震, 对近几年川滇地区的显著地震和长时间尺度振幅结果并无详细研究。由此, 笔者使用该方法对四川省内近年来显著的2次强震进行回溯性研究, 主要分析了康定地震、九寨沟地震震前半年及震后3个月内4.4h周期分量振幅的变化特征, 并根据中强地震目录所得到的川滇地区中强震活跃期和中强震平静期, 研究台站4.4h周期长时间尺度的振幅变化情况。

1 计算方法及处理过程

对WVD分析的定义为

$ W\left({t, f} \right) = \int_{ - \infty }^{ + \infty } {s * \left({t - \frac{1}{2}\tau } \right)} s\left({t + \frac{1}{2}\tau } \right){{\rm{e}}^{ - j2{\rm{ \mathsf{ π} }}f\tau }}{\rm{d}}\tau $

(1)

式中, s*(t)为s(t)的共轭；τ为信号的时间延迟；f为信号的频率。

对多分量信号而言, 使用WVD会出现交叉项, 以致产生虚假信号, 因此需要添加核函数来抑制交叉项, 常用对WVD进行加核函数的有伪魏格纳-维勒分布(PWVD)和平滑伪魏格纳-维勒分布(SPWVD)。

其中, SPWVD分布为

$ W_Z^S\left({t, f} \right) = \int_{ - \infty }^{ + \infty } {h\left(\tau \right)s * \left({t - u - \frac{1}{2}\tau } \right)} s\left({t - u + \frac{1}{2}\tau } \right)g\left(u \right)h\left(\tau \right){{\rm{e}}^{ - j2{\rm{ \mathsf{ π} }}f\tau }}{\rm{d}}\tau $

(2)

在对台站数据进行时频分析前, 首先对Z分量小时值数据进行预处理。采用K-L最佳直线拟合法去除Z分量的长期变化并实现标准化；对去除线性趋势后的Z分量数据使用三点线性平滑滤波进行处理, 将信号分离为主要周期为6~24h和2~6h的两部分；然后对2~ 6h部分的数据进行Hilbert变换, 构建解析信号以期消除由负频率成分产生的交叉干扰项, 最后再进行SPWVD时频分析。

根据前人的研究及震例成果, 本文只选用2~6h中的4.4h对四川及其邻区的地磁数据进行时频分析。同时, 考虑到外空磁场活动会对4.4h振幅变化产生较大的影响, 在进行震例研究之前, 根据同时段∑K_p指数数值, 将每天∑K_p指数大于20的4.4h振幅高值干扰去掉, 设定成背景值, 从而消除外空磁场的影响。本文使用地震及台站如表 1、图 1所示。

2 结果 2.1 康定地震时频分析

2014年11月22日在鲜水河断裂带中南段四川康定县发生了M_S6.3地震, 此次地震为带有微小正断分量的左旋走滑型地震。震中附近500km范围内能够提供分析数据的地磁台站共8个, 如图 1所示。由于巴塘数据质量不太稳定, 马边台在此期间无数据, 故未使用这2个台站。

笔者统计了康定地震前6个月和震后3个月震中附近地磁台站4.4h周期的计算结果。从结果看, 9月21日8个台站的幅值很明显地同步出现了高值的异常现象(图 2), 其中5个台站的当日幅度为震前6个月的最高值, 当天距离发震的时间为52天。当日磁指数ΣK_p=6, 为磁静日, 外空场对其影响忽略不计, 异常可信。具体统计异常信息如表 2所示。

2.2 九寨沟地震时频分析

2017年8月8日在青藏高原东缘四川省九寨沟县发生M_L7.1强烈地震, 此次地震呈左旋走滑类型, 发震构造为位于岷江断裂与塔藏断裂之间的树正断裂(易桂喜等, 2017)。震中附近400km范围内能够提供分析数据的地磁台站共9个, 由于巴塘数据质量不太稳定, 玛曲和临夏台在此期间缺数较多, 故未使用这3个台站。

笔者统计了九寨沟地震震中附近地磁台站地磁数据的前6个月和震后3个月4.4h周期的时序曲线结果。从结果来看, 5月14日6个台站(成都、崇州、江油、合作、舟曲、天水)同步出现了高值异常现象, 其中, 成都、崇州、江油3个台站的幅值为震前几个月的最高值, 其余台站当日幅值较高, 但高值并不突出(图 3)。同步异常时间距离发震时刻为53天。当日ΣK_p=17, 为磁静日, 外空场对其影响忽略不计, 异常可信。具体统计异常信息如表 3所示。

2.3 中强地震活跃时段和平静时段时频结果分析

笔者从中国地震台网中心的地震目录中统计了2014年1月~2019年2月四川省及其邻区5.5级及以上地震信息(图 4), 从图 4可以看出, 2014年四川及邻区一年时间内连续发生了5次M_S6以上强震(且都发生在川滇地区)。而2015年~2017年8月川滇地区仅有1次5.5级地震(理塘M_L5.5)发生。因此, 可以明显分辨出2014年川滇地区的强震活跃(表 4、图 4)和2015、2016年相对的强震平静。笔者对2014年强震活跃期和2015年强震平静期、2017年12月中强地震平静期至2018年12月开始的中强地震活跃期的川滇地区地磁台站的4.4h振幅时序曲线进行了研究, 取得了一些新认识。

2.3.1 2014~2015年4.4h时间序列

笔者选择了川滇地区8个质量较好、数据较连续的地磁台站进行分析, 台站分布范围较广, 能够满足区域性研究的样本要求。这里, 笔者尝试给出定量描述该方法地震前兆的异常指标。从各台4.4h周期时间序列结果来看, 各台基于背景振幅的高值相差较大, 如果统一定义异常幅度会造成同步异常次数过多, 虚报和误报率也高。因此, 提出了基于每个台站的实际高值和背景值设定的相对合适的异常阈值, 这一点与前人的判定标准不同(何康等, 2014)。本文各台使用阈值如表 5所示。

从各台4.4h周期的结果(图 5)来看, 在川滇地区发生地震频次高、震级大的2014年, 4.4h周期的幅值较高。其中, 西昌、道孚、平地、永胜的幅值是该时间段内的最高值。

2015年部分台站也出现高值, 但这种高值的次数和台站的数量低。此外, 从多台同步性来讲, 2014年出现了多次同步高值异常, 这里定义的多台同步异常为区域内至少4个地磁台站出现超阈值异常变化。2014年7月6日、2014年8月26日、2014年9月21日均出现了同步高值异常, 当日均为磁静日, 分析认为是强震前的异常反应；2014年12月~2015年3月, 楚雄等5个台站(图 5淡蓝色区域)在1月份的多天同步出现较高数值的幅值异常, 分析认为可能是川滇菱块内2015年1月14日四川金口河M_L5.2地震的前兆异常(表 6)。2015年3月~2015年12月期间, 几个台站4.4h周期的幅值在此阶段内的最高值较小, 出现高值的次数也少。从表 6还可以看出, 2014年川滇地区6级以上强震前均有同步超阈值异常反应, 且异常幅值最大的台站距离震中较近。

2.3.2 2017年12月至2018年12月的4.4h时序曲线分析

从地震活动性(图 4)可以看出, 2017年底~2018年9月初是川滇地区的地震平静期, 期间该区域无5级以上地震发生。从所选台站的4.4h周期结果来看, 2017年12月~2018年7月各台的振幅有明显增大的趋势, 至2018年7月4号同步达到所选时间段的最高值, 异常出现后的65天(3个月内), 即2018年9月8日发生了云南墨江M_L6.2地震。而后4.4h振幅变化较弱, 2018年10月1日又出现了1次较高幅值的同步异常, 异常出现后的30天(3个月内), 即10月31日发生了西昌M_L5.4中等地震。之后幅值变小、变化也较小并相对平稳发展, 后续无5.5级以上中强震发生(图 6)。

从2017年12月地震平静期到2018年12月地震活跃期的震例统计可以看出, 平静时段内该区域还是出现了4次虚报现象(无M_S≥5.0地震对应), 但同步异常数量较低, 异常最大振幅也较低(表 7)。

3 结论与讨论

本文使用平滑伪魏格纳-维勒时频分析方法, 对康定地震、九寨沟地震进行震例回溯性分析, 以期为川滇地区中强地震预测和区域地震活动性分析提供一定依据。研究发现：康定地震之前, 附近台站4.4h周期的幅值存在8个台站同步高值异常, 异常时间距离地震发震时刻为52天, 异常出现时间在地震前3个月内, 显示出短临特征。九寨沟地震前附近台站4.4h周期存在6个台站同步高值异常, 时间上也显示短临特征。但是九寨沟震例高值的普遍性特征不太突出, 只有成都、崇州、江油等3个台站在震前有大幅度的高值现象。结合以往统计的该区域强震震例，认为在川滇地区强震震前3个月内周围台站可能出现多台同步高值异常现象(异常指标初步定为4个台站), 高值可能是该时间前后3个月乃至1年尺度内的极大值, 强震震前同步异常台站分布范围较广。有些强震震前高值异常频次会增多(芦山地震、岷县地震), 但统计显示大多数强震前多台同步高值异常是最普遍的现象, 且时间上呈现短临特征。中等地震前3个月内周围台站可能也会出现多台(至少4台)同步幅值较高的异常现象, 短时间内幅值较高的频次可能也会变大。

此外, 笔者还从中强地震活跃期和中强地震平静期2个时段出发, 探讨了2014年、2018年4.4h周期的幅值在活跃期和平静期的变化特征。研究发现：中强地震活跃期间, 4.4h周期的幅值会增大, 多出现4台以上同步的高值变化特征, 其中强震震级越大, 同步变化的台站数量越多。在地震平静期间, 4.4h周期幅值较小且变化较平稳。因此, 从已有地震活动性强弱的时段划分来研究川滇地区地磁台站4.4h幅值变化特征, 可以为预测区域地震活动的强弱提供参考。

何康等(2013、2014)指出, 地磁Z分量滤波后4.4h周期在磁静日的时频分布一致性较差, 在扰日期间相距较近的地磁台站时频分布一致性较好, 说明4.4h周期的磁场既体现了区域性特征, 同时也受外空场的影响。因此, 为了说明各地磁台4.4h周期的变化特征是与地震孕育和发生密切相关的前兆异常信息, 使用了同时段K_p指数来排除外空场的影响。由此可见, 地磁4.4h周期分量信号主要体现了地壳磁场和感应磁场的变化, 前者是由地壳和上地幔磁性岩石产生的, 其大小基本不随时间而变化, 但是与发震直接相关的块体构造运动能够引起区域性的磁性变化；后者是外源场在地球内部生成的感应电流的磁场, 直接产生感应场的源是分布在地壳和地幔中的感应电流。总体来看, 感应场由于包含了更多的地壳和地幔中的信息, 在孕震过程中可能导致观测的地磁数据出现异常变化, 从而包含了地磁场4.4h周期出现变化, 这种变化突出表现在震前区域性地磁台站时频分析后该周期的幅度、频次以及台站出现异常的时间同步性上。

地震的孕育和发生过程中会伴随着地磁场的变化, 但是对震磁现象物理机制方面的研究尚处于探索阶段。目前形成了多种假说来解释上述变化, 如感应磁效应、压磁效应、动电效应、热磁效应及微破裂等机制。其中, 动电效应和微破裂机制被认为是对震磁关系最有可能解释得通的机制。微破裂机制认为在地震孕育过程中, 震源区的应力缓慢积累, 岩石受到不同程度的挤压或拉伸, 以至于岩石的形态和赋存空间发生变化。Pride(1994)进一步研究认为岩石空隙中流体是产生电磁变化的重要原因。微裂隙扩展时, 间隙的流体会侵入新生及扩展裂隙中, 导致震源区岩石的电性或磁性发生改变, 进而导致地球内部的感应磁场发生变化。

虽然SPWVD时频方法在对相关震例的地磁数据分析取得了较明显的效果, 但是需要说明的是, 这些特征都是在已发地震后进行回溯总结时所统计的；该方法在地震平静期也存在一定程度漏报地震的情况。因此, 以该方法来准确预报地震(地震的时间、地点、震级)还有较远的距离。另一方面, 有些台站在选取震例前后的时间段内存在缺数甚至长时间缺数的情况；还有些台站得地磁数据受到较严重的环境破坏和人为干扰, 数据质量存在问题。这些情况给震例分析和研究带来一定困难。总的来说, 通过使用SPWVD时频分析方法对川滇地区地磁台站进行时频分析, 得到4.4h周期分量时序曲线变化特征, 多台站同步异常的幅值、频次、数量和异常出现的时间可以对整个川滇地区(川滇菱形块体及其附近区域)未来3个月至1年的地震活动性变化趋势提供统计预测, 对地震形势做宏观判定, 为地震预测预报服务。

致谢: 在论文成稿之前得到了四川省地震局诸多专家的指导, 感谢匿名专家提出的修改意见。数据来源于国家地震前兆台网中心。