0 引言

地震的孕育和发生是一个复杂过程, 这一过程引起的电磁场变化, 早在100多年前便已引起关注(丁鉴海等, 2006; 钱家栋, 2010)。随着观测密度和精度的大幅提升, 与地震相关的电磁场变化报道增多(Fraser-Smith et al, 1990; 关华平等, 1995; Park, 1996; Nagao et al, 2002; 郝锦绮等, 2003; 杜爱民等, 2004; Hattori, 2004; 丁鉴海等, 2006; 顾左文等, 2006; 李琪等, 2006, 2016; 杜学彬等, 2008; 陈化然等, 2009; 汤吉等, 2010; 钱家栋, 2010; 赵国泽等, 2010; Uyeda, 2013; 张建国, 2017; Enomoto et al, 2017; Zhang et al, 2017; 解滔等, 2018; 姚休义等, 2018; Lei et al, 2018; Liu et al, 2019), 地震电磁研究逐渐发展成相对独立领域。地震电磁异常通常出现在发震前不久, 对于地震短临预测及防震减灾具有潜在的应用价值(Hattori et al, 2018)。因此, 20世纪末至今, 震前电磁信号一直是研究热点, 各国投入大量的人力物力, 经过数十年的观测和研究, 已初步证实孕震过程中电磁现象的客观存在。

在实际观测中, 电磁数据包含各种场源信号：电离层磁层扰动产生的电磁信号(空间场源信号); 直流电车、汽车等产生的电磁信号(人文干扰信号); 地震活动产生的电磁信号(地震电磁信号)。其中：空间场源信号最强, 在较大空间范围内(可达上千千米)具有一致性, 是一种广域性信号; 人文干扰信号是一种局部信号, 可在数十千米范围内被观测到, 且信号在人文活动较少的夜间相对较弱; 地震电磁信号由于在地下传播过程中的衰减, 多被认为是一种局部性信号。已有研究表明, 地震电磁信号通常比前两者弱(Hattori et al, 2004, 2013a)。为了有效识别、提取地震电磁信号, 研究者提出一系列数据处理方法, 如：小波分析(李琪等, 2006; Han et al, 2011)、主成分分析(Uyeda et al, 2002; 韩鹏等, 2009)、分形分析(荣扬名等, 2012)、盲源信号分离法(周媛媛等, 2019)等。Han等(2014)提出一种基于远场参考的地震地磁数据处理方法, 该方法基于目标台站与观测台站信号能量的相关性, 利用最小二乘法构造预测模型, 利用参考台站数据计算目标台站信号能量的理论值, 并与实际观测结果对比, 从而实现目标台站磁场异常检测及识别。本研究以2003年1月日本茨城县5.1级地震为例, 探讨远场参考方法在超低频(ULF)地震地磁数据处理中的应用, 进而检测此次地震发生前局部磁场异常。

1 观测数据

日本茨城县位于太平洋板块、菲律宾板块、北美板块交汇处, 地震频发。尤其是茨城县东部海域, 由于太平洋板块的俯冲下插, 地壳活动剧烈, 几乎每年有5级以上地震发生, 是一个天然的俯冲带地震研究试验场。因此, 日本气象厅1912年在茨城县建立Kakioka观象台(下文简称KAK台), 陆续开展地磁、地电、地震、大气电场等地球物理观测(http://www.kakioka-jma.go.jp/intro/enkaku.html#taishou), 迄今已积累上百年观测资料。

2003年初, 在茨城县东部海域发生一组中强地震, 最大地震震级5.1, 震源深度46 km。在此次地震发生前约1周, 同一深度处发生4.7级地震。2次地震与KAK台震中距分别为77 km和86 km。据调查, 2002年底距KAK台约18 km处发生一次4.2级地震, 震源深度54 km。3次地震详细信息见表 1, 信息均来自日本气象厅https://www.data.jma.go.jp/svd/eqev/data/ bulletin/index_e.html)。3次地震及KAK台空间分布见图 1。

为了考察上述地震前可能存在的磁场变化, 对KAK台站磁场观测数据展开分析。同时选用距KAK台站一千多千米的Kanoya台站(下文简称KNY台)作为远场参考, 以排除空间活动干扰。KAK台与KNY台均隶属于日本气象厅, 观测数据质量较高, 尤其是KAK台站, 环境噪音较小, 是国际基准台。2个台站均使用磁通门磁力仪(采样率为1 Hz, 分辨率为0.01 nT), 对地磁场绝对值进行连续测量。

选取2002年11月1日至2003年3月1日KAK台地磁观测数据(http://www.kakioka-jma.go.jp/cgi-bin/plot/plotSetNN.pl), 涵盖3次地震发生前后各1个月的数据, 以KNY台作为远场参考, 分析超低频(ULF)地震地磁数据处理方法的应用。

2 远场参考法

地表观测的电磁信号包含多种场源, 主要为空间活动及人文干扰产生的信号。对于空间场源信号, Hattori等(2013b)选用多个数据质量高的国际基准台, 考察了超低频段空间场源信号在不同区域的特征, 发现在0.001—0.100 Hz频段, 空间场源信号在上千千米范围内仍具有较好一致性。对于人文干扰信号, Han等(2014)分析了KAK台站记录的1天磁场数据频谱, 发现在凌晨电车停运时段, 能量背景明显降低。因此, 选取每日2:30—4:00的磁场数据进行分析, 以减小人文噪音, 同时选用KNY台站作为远场参考。

为了高效、准确地检测并识别地震电磁信号, 对子夜时段磁场数据进行离散小波分析。在以往研究中, 地震电磁信号优势频段在0.01 Hz上下(Hattori et al, 2018)。经小波变换, 提取KAK台站磁场数据(0.010±0.005)Hz频带的信号, 并计算每日该段信号能量, 与参考台站KNY所得结果进行对比。Han等(2014)分析发现, KAK台站磁场能量变化与远参考台KNY台站高度一致, 相关系数高达0.94。基于此, 利用最小二乘法建立KAK台站磁场能量变化模型(Han et al, 2014), 公式如下

$ \mathit{Z}_{{\rm{KAK}}}^{\rm{r}} = \beta \cdot {Z_{{\rm{KNY}}}} + c $

(1)

其中, Z_KNY为参考台站KNY磁场能量变化的观测值, Z_KAK^r为目标台站KAK磁场能量变化的预测值, β和c为基于2个地磁台实际观测数据计算的线性系数。通常, 若只存在空间场源信号, 观测值与预测值具有高度一致性。当孕震过程中出现局部异常磁场信号时, 可以发现观测值相对于预测值出现偏离。通过计算二者之比, 可较好检测出异常磁场信号, 公式如下

$ P = \frac{{{Z_{{\rm{KAK}}}}}}{{Z_{{\rm{KAK}}}^{\rm{r}}}} $

(2)

式中, Z_KAK为KAK台站磁场能量变化的观测值。

3 结果

对KAK台站及参考台KNY台站2002年11月1日至2003年3月1日磁场Z(垂直)分量数据进行分析, 获得子夜(2:30—4:00 LT)磁场信号在(0.010±0.05)Hz频段的每日能量数值。基于式(1), 利用参考台KNY所得结果, 计算KAK台磁场能量变化的预测值Z_KAK^r, 并在此基础上, 利用式(2), 计算KAK台磁场能量观测值与预测值的比值P, 结果见图 2。

图 2(a)中红色曲线为KAK能量变化预测值, 蓝色曲线为其观测值。由图 2(a)可见, 在大多数时间内, 观测值与预测值具有较好一致性, 这是因为, 绝大多数磁场信号为空间场源信号, 其能量在KAK台和远参考台KNY具有均一性。但在2002年12月中旬及2003年1月上中旬, 2条曲线出现明显不一致, 观测值大幅高于预测值, 说明此时段KAK台磁场信号除空间场源信号, 还混有其他局部场源信号。图 2(b)展示了观测值与预测值之比P的变化, 可见在大多数时间内P值在1处上下波动, 说明观测值和预测值大多数时间较吻合, 对比图 2(a)可知, 观测值高于预测值时间段, P值明显高于1。因此, P值可以较好地反映观测值的异常升高, 进而识别、检测出局部磁场异常信号。图 2(c)展示了KAK台磁场能量预测值及观测值的10天滑动均值变化。图 2(d)展示了P值的10天滑动均值变化, 其中紫色虚线为P = 1参考值。显然, 10天滑动均值可以抑制随机扰动, 使得异常更加明显、聚焦。

图 2(e)展示了分析时段内KAK台100 km范围内2003年1月10、21日分别发生的4.7级、5.1级地震。图 2(f)展示了KAK台100 km范围内每日地震能量E_S, 计算公式如下(Hattori et al, 2006)

$ {E_{\rm{S}}} = \sum\limits_{{\rm{1day}}} {E_{\rm{S}}^\prime } $

(3)

$ E_{\rm{S}}^\prime = \frac{{{{10}^{4.8 + 1.5M}}}}{{{r^2}}} $

(4)

式中, M为地震的震级, r为震源到台站的距离。已有研究表明, 地震电磁信号能否被检测到不仅与地震大小相关, 还与地磁台站到震源的距离相关(Hattori et al, 2004)。由于E_S指数不仅考虑了震级影响, 也考虑了台站到震中的距离, 应比单纯震级参数存在优势。由图 2(e)可知, 2个4.5级以上地震发生前几天至1周, 地磁场出现明显的异常变化。但2002年12月中旬出现的磁场异常, 似乎无地震对应, 但由图 2(f)可见, 在2002年12月下旬, 异常出现大约1周后, 在距KAK台仅约18 km处发生一个4.2级地震。

为了进一步查清KAK台磁场异常出现原因, 对该台原始数据波形进行分析, 将震前异常时间段(2003年1月8号)的磁场数据与震后无异常时间段(2003年2月24号)数据进行对比, 结果见图 3, 其中红色曲线为KAK台数据, 蓝色曲线为远参考台KNY数据, 磁场X为南北分量, Y为东西分量, Z为垂直分量, 为了便于比较, 均值已剔除。由图 3可见：在无异常时间段, 参考台和KAK台磁场三分量变化具有较好一致性; 在异常时间段, KAK台磁场Z分量出现明显的脉冲式异常信号, 参考台KNY无类似变化。由此表明, 图 2检测到的磁场异常确为KAK台出现的局部磁场异常。

4 讨论

上述分析结果表明, 远场参考法可有效识别、检测磁场局部异常。其关键是寻找合适的参考台站。选取原则如下：不能太远, 以便目标台站与观测台站空间场源信号具有一致性, 磁场能量变化具有较好相关性; 不能太近, 使得目标台站出现的磁异常在参考台能衰减到可以忽略, 参考台站数据不含目标台局部异常信号。在实际中, 需要根据能量相关性选取合适的参考台站, 既要求相关系数足够高, 背景变化一致, 又要求相关系数不能等于1, 且异常在2个台站表现不同, 这样才能在观测背景中识别局部异常变化。

由式(1)可知, 计算目标台站磁场能量预测值的关键是估算正确的β和c值。当数据量较少时, 直接用最小二乘回归可能存在较大偏差, 因为局部磁场异常会对回归结果产生影响。此时, 应选取不含局部磁场异常的时间段数据, 即可准确估算2个台站空间场源信号的能量相关系数, 并计算β和c值。此为, 可使用加权最小二乘, 降低局部磁场异常数据在回归系数计算时的权重。文中分别使用以上2种方法计算β和c值, 所得结果几乎相等, 最大相对误差不超过2%。

在本研究中, 考察远参考台与目标台站的磁场信号的能量, 由于能量不包含相位信息, 计算简便, 可以快速识别磁场异常。但正因为其不含相位信息, 难以实现异常信号的精确提取。Harada等(2004)提出基于参考台的转换函数法, 通过构建目标台站和参考台之间的转换函数, 实现局部磁场异常信号的提取。周媛媛等(2019)提出基于参考台的盲源信号分离法, 可实现局部磁场异常信号的提取。将本研究方法与上述2种方法相结合, 开发更有效的地震地磁信号识别、检测、提取方法, 将是一项值得深入开展的工作。

由图 2可见, 若只考察地震震级, 将会出现一些对应不佳的异常。若考察震级, 同时考虑震中距因素, 则异常与地震对应效果变好。这是因为, 距台站较近的某些地震, 虽然震级不大, 仍可能产生较弱的局部磁场异常。因此, 在分析磁场异常时, 不能只关注大地震, 还应该对台站周边震中距较小的地震有所考量。此外, 2003年1月茨城县发生2次中等地震, 难以准确判断异常与哪个地震对应。为了获得异常出现的一般规律, 需要基于大量数据开展统计分析。Han等(2014)分析了KAK台时长10年的磁场数据, 发现异常多出现在地震发生前1—2周。本研究检测的异常多出现在震前约1周, 与之前的统计分析具有较好的一致性。

通过考察磁场异常的原始波形数据, 发现Z分量异常更明显, 见图 3, 说明垂直分量与水平分量磁场幅度相比的谐波振幅比法(郭明瑞等, 2019)或极化法(Hayakawa et al, 1996)存在一定合理性, 且异常信号多在相对高频出现, 表现为脉冲式类噪音信号, 与Uyeda等(2002)对伊豆群岛群发地震相关磁场信号波形的分析结果一致性较好。Dunson等(2011)曾在美国加州“QuakeFinder”地震电磁监测台网观测到类似信号。

5 结论

本研究以2003年1月日本茨城县5.1级地震为例, 探讨远场参考在超低频(ULF)地震地磁数据处理中的应用。通过对茨城县KAK台及远参考台KNY的磁场数据进行分析, 可以得到以下结论：①KAK台在地震发生前约1周存在明显的磁场异常; ②KAK台2002年12月中旬出现的异常可能与相距仅18 km的4.2级地震相关; ③分析异常出现前后的原始波形数据发现, 磁场Z分量异常显著; ④异常出现时间与已有统计结果具有较好的一致性。

上述结果表明, 2003年1月日本茨城县5.1级地震发生前存在局部磁场异常, 验证了远场参考方法的有效性, 相关成果可为地震地磁数据处理提供参考。

研究使用的参考台地磁数据由日本气象厅提供, 在此表示感谢。