Correlation characteristics of the geoelectric field with multi-azimuth and multi-pole distance for the same field
-
摘要: 针对5个典型地震震中周边的50余个场地记录到的地电场数据展开分析,计算了各场地不同方位之间的地电场相关系数δxy或同方位不同极距之间的地电场相关系数δxx,以此表征场地地电场变化的复杂性。结果显示:地震活动平静期间这些场地的地电场相关系数δxy和δxx 具有相对稳定性,中强地震孕育发生期间震中附近多数场地的δxy和δxx会出现明显下降,且这种现象具有方向性、时间准同步性,并随着震中距的增加而逐步消失。从机理上推断,δxy和δxx的下降反映了震中附近局部场地地电场变化的非均匀性。Abstract: At a given site, the correlation coefficient δxy between geoelectric field with different azimuths, or δxx between geoelectric field with differernt polar distances for the same azimuth can express the complexity of site geoelectric field variation. This paper analyzed the geoelectric field correlation coefficient δxy and δxx based on the geoelectric field data from more than fifty sites around five typical earthquake epicenters. The results show that δxy and δxx of a site are relatively stable during a quiet period of seismic activity. During a moderate earthquake, the δxy and δxx for most of the sites near the epicenter decrease significantly. This phenomenon is directional and quasi-synchronous in time, and gradually disappears with epicentral distance increasing. Furthmore, in viewpoint of mechanism, it is deduced that the decrease in δxy and δxx reflects the non-uniformity of geoelectric field variation in the local sites near the epicenter.
-
Keywords:
- geoelectric field /
- multi-azimuth /
- multi-pole distance /
- correlation
-
引言
地震发生时,地震预警系统(earthquake early warning system,缩写为EEWs)可以在破坏性地震波到达前向预警目标区发布预警信息,这对于减轻人员伤亡和减少次生灾害至关重要(Allen,Kanomori,2003)。目前,世界上地震较活跃的国家和地区都已经建立或者正在建立地震预警系统,如日本(Kamigaichi et al,2009)、美国(Cua et al,2009)、意大利(Zollo et al,2009)、中国福建(Zhang et al,2016)和中国台湾(Wu,2002)地区等。震级作为地震预警系统中的关键预警信息,其准确性对地震预警系统极为重要。地震预警系统主要依据震级的大小来评估地震对目标区域的影响程度和范围,因此,如何从一个正在发生的地震事件中较为准确地估算出震级,是地震预警系统研发中亟需解决的关键难题之一。
初至地震波(通常为3 s以上)包含了与震级相关的信息(Kanamori,2005;Zollo et al,2006),因此,地震预警中震级估算普遍利用初至地震波的特征参数来建立预测震级的经验公式。这些特征参数在频域和时域内代表初至地震波中与震级相关的信息。在频域内,一般认为震级越大地震波周期越长,据此研究人员提出了卓越周期τ p (Nakamura,1988)、最大卓越周期${\tau ^{\max }_{\mathrm{p}}} $ (Allen,Kanamori,2003)和平均卓越周期τc (Kanamori,2005)等参数;在时域内,相同震源距下震级越大地表变形越大,据此研究人员提出了位移幅值P d (Wu,Zhao,2006)、速度平方积分(integral of the squared velocity,缩写为IV2)(Festa et al,2008)和累积速度绝对值(cumulative absolute velocity,缩写为 CAV)(Alcik et al,2009)等参数。除了前述早期提出的经典特征参数外,近年来也不断提出新的特征参数,如平均对数周期τlog (Ziv,2014)和累积位移绝对值(cumulative absolute displacement,缩写为CAD)(马美帅等,2022)等,这些新的特征参数改善了震级预测效果,并促进了震级预测理论的发展。然而,震级估算十分复杂,需要综合考虑震源、传播路径、场地条件、监测仪器等多方面因素的影响(林彬华等,2021),尽管单一特征参数与震级有一定的相关性,但其仅能代表初至地震波中与震级相关的某方面信息,较难反映出地震的全部特征,从而导致震级估算的效果不稳定(Shieh et al,2008;Ziv,2014)。为了改善震级估算效果,一些学者提出采用两个或两个以上特征参数来估算震级,例如将τc与P d相乘来判断地震是否具有破坏性(Wu, Kannmori,2005),或将${\tau ^{\max }_{\mathrm{p}}} $和τc各自估算震级的平均值作为最终震级(Hsiao et al,2009),或采用${\tau ^{\max }_{\mathrm{p}}} $和P d各自估算后取平均震级(Chung et al,2019),以及将频域参数${\tau ^{\max }_{\mathrm{p}}} $和时域参数CAD耦合为新的特征参数Sdτ等(Wang et al,2022)。这些研究均表明了利用更多的特征参数比采用单一特征参数更有利于震级估算。
近年来,随着人工智能的兴起,机器学习在地震预警领域的应用进入了新阶段,在震相捡拾(Perol et al,2018;Wang et al,2021)、地震定位(Perol et al,2018)、事件判别(Li et al,2018)、地震动峰值预测(Hsu et al,2013)等方面都取得了显著的成果。相较人为定义的经验公式,机器学习的优势在于可以自动地建立更为复杂的输入与输出关系,实现更为精准的预测。机器学习的震级估算方法,根据输入的不同可以分为两类:一类是直接将初至地震波作为输入,自动从初至地震波中提取特征,实现端到端的震级估算,以深层卷积神经网络最为典型(Mousavi,Beroza,2020;位栋梁等,2022;Wang et al,2023);另一类是将五个以上的特征参数(时域参数和频域参数)作为输入,融合初至地震波的多方面信息从而实现震级的估算,主要有浅层人工神经网络(Böse et al,2012;杨黎薇等,2018)、支持向量机(朱景宝等,2021)等方法。尽管当前的机器学习方法显著地改善了震级估算效果,但忽略了地震预警系统估算震级过程中的信息变化。在地震预警系统的实际应用过程中,随着触发台站数量的不断增加,震源距估算是从无到有,从不准确到准确的过程。目前,机器学习方法并不能考虑震源距在地震预警系统的前述变化过程。
本文基于机器学习领域的高斯过程回归(Gaussian processes regression,缩写为GPR)提出了一种多输入的震级估算方法。该方法包括两个版本:无需震源距参数的GPR-M和利用震源距参数以提高估算准确性的GPR-M-R方法。拟利用日本的强震记录对该方法进行训练和测试,并采用智利的强震记录测试GPR方法的泛化能力,之后再通过中国的典型地震事件检验GPR方法的实际应用效果。此外,将GPR方法与广泛使用的${\tau ^{\max }_{\mathrm{p}}} $方法和P d方法进行对比,以进一步评估其性能,以期为地震预警提供一种更准确、更灵活的震级估算方法。
1. 地震数据
本文选取日本Kiban-Kyoshin Network (KiK-net)和智利SIBER-RISK (Simulation Based Earthquake Risk and Resilience of Interdependent Systems and Networks)数据库的地面加速度记录开展算法的研究和测试,在选取加速度记录时遵循以下原则:① MW≥4.0;② 因为地震预警系统存在“预警盲区”的问题(郭凯等,2016),为保证初至地震波中至少具有3 s P波用于估算震级(Kanamori,2005;Wu,Zhao,2006;Hsiao et al,2011;金星等,2012;Zhang et al,2016;Cuéllar et al,2018;Kohler et al,2020),且尽可能地包含日本近海地震事件(王延伟等,2020),所以震源距为25—200 km;③ 地震记录的三分向合成加速度的峰值需大于2 cm/s2 (Cheng et al,2014),且信噪比不小于10 dB (Perol et al,2018)。对已筛选出的记录进行以下常规处理:① 检查基线并校正;② 采样频率统一为100 Hz;③ 自动捡拾P波到时并人工校验。经以上数据选取和处理后,从日本KiK-net数据库总共筛选出3万3 698条竖向加速度记录,包含了1997年10月21日至2021年12月18日的3 271次MW4.0—9.0地震事件。将筛选出的KiK-net强震记录数据按照7 ∶ 3的比例随机划分为训练数据集和测试数据集,其中训练数据集含2万3 588条记录,测试数据集含1万零110条记录。从智利SIBER-RISK数据库中总共筛选出5 353条竖向加速度记录,包含了1985年3月3日至2021年7月21日的1 033次MW4.0—8.8地震事件,用于测试GPR方法的泛化能力。训练、测试和泛化数据集的记录数量的关系如图1所示。
2. 高斯过程回归震级估算方法
2.1 高斯过程回归
GPR是一个非参数模型,具有严格的统计学习理论基础,对处理高维、非线性等复杂问题具有很好的适应性,同时GPR还具备避免过拟合、泛化能力强、超参数自适应获取以及输出具有概率意义等优点,这使得GPR成为机器学习领域的重要方法(Seeger,2004;Rasmussen,Williams,2006;Park et al,2011),并在时间序列分析(Rohani et al,2018)、图像处理(Hua et al,2023)和自动控制(Nguyen-Tuong et al,2008)等领域中取得了显著成果。GPR的基本原理是通过协方差核函数的运算,将样本数据映射到高维空间,基于训练数据的训练过程就是对超参数的选择过程。关于GPR的详情请参阅Seeger (2004)与Rasmussen和Williams (2006),此处不再赘述。GPR的网络结构如图2所示。
2.2 基于GPR的震级估算
在地震预警系统中,随着触发台站的增加,震源距参数也被不断更新,因此震级估算方法需适应震源距的变化(Wu,Kanomori,2005)。据此,本文在有、无震源距两种情况下估算震级,将特征参数作为输入的震级估算定义为GPR-M方法,将特征参数和震源距作为输入的震级估算定义为GPR-M-R方法,两种方法都将预测的震级作为输出。本文选取频域、时域和时频域三类共十个特征参数作为GPR-M和GPR-M-R方法的输入。
频域类参数估算震级不需要震源距,估算震级的依据是:震级越大,地震波的长周期成分越丰富(Nakamura,1988)。将以下三个典型频域参数作为特征参数输入GPR-M方法:
1) 最大卓越周期$ {\tau ^{\max }_{\mathrm{p}}}$。Allen和Kanamori (2003)提出利用卓越周期τp (Nakamura,1988)的最大值估算震级,即${\tau ^{\max }_{\mathrm{p}}} $方法,计算公式为:
$$ \left\{\begin{array}{l} {\tau }_{{\mathrm{p}}}=2 \text{π} \sqrt{\dfrac{{X}_{\mathrm{i}}}{{D}_{i}}}\text{,}\\ {\tau }_{\mathrm{p}}^{\mathrm{m}\mathrm{a}\mathrm{x}}=\mathrm{max}\left\{{\tau }_{{\mathrm{p}}}\right\} \text{,} \end{array}\right. $$ (1) 式中,$ {X_{{i}}} = \alpha {X_{{{i }- 1}}} + x_{{i}}^{\text{2}} $,$ {D_{{i}}} = \alpha {D_{{{i }- 1}}} + ( {\text{d}}x/{\text{d}}t ) _{{i}}^{\text{2}} $,τp为i时刻测定的卓越周期,xi和Di分别为平滑后的竖向速度值及其导数,α为平滑系数,一般取0.99 (Allen,Kanamori,2003)。
2) 平均卓越周期τc。Kanamori (2005)提出利用卓越周期τp的平均值来估算震级,即τc方法,计算公式如下:
$$ {\tau }_{\mathrm{c}}=2 \text{π} \sqrt{\frac{{\displaystyle\int }_{0}^{\mathrm{\Delta }t}{u}^{2} ( t ) \mathrm{d}t}{{\displaystyle\int }_{0}^{\mathrm{\Delta }t}{v}^{2} ( t ) \mathrm{d}t}} \text{,} $$ (2) 式中:u为地震波位移;v为地震波速度;Δt为地震波时长,单位为s。
3) 平均对数周期τlog。Ziv (2014)提出使用功率谱的平均对数周期τlog估算震级,其计算公式为
$$ \lg {\tau }_{\mathrm{l}\mathrm{o}\mathrm{g}}=\frac{\displaystyle\sum\limits _{{i}}\left({P}_{{i}}^{\mathrm{*}} ( {w}_{{i}} ) {{\mathrm{lg}}}\dfrac{1}{{w}_{{i}}} \right)}{\displaystyle\sum\limits _{{i}}{P}_{{i}}^{\mathrm{*}} ( {w}_{{i}} ) } \text{,} $$ (3) 式中,w为功率谱的频率,${P}_{{i}}^{\mathrm{*}} $为0.1—10 Hz频率范围内每0.1个log (频率)的功率谱纵坐标值。
时域类参数估算震级需要已知震源距,估算震级的依据是:在相同震源距情况下,一般震级越大,地表变形越大。以下六个典型时域参数作为特征参数输入GPR-M方法:
1) 位移、速度和加速度的幅值。Wu和Zhao (2006)首先利用初至3 s P波的位移幅值P d预测震级,并且给出了P d、震源距和震级的经验公式;速度幅值P v和加速度幅值P a与地表变形的最大速度相关(Wu,Kanamori,2005),因此,P v和P a与震级存在相关性。
2) 速度平方积分IV2。Festa等(2008)基于地表变形的能量变化提出IV2估算震级,计算公式为:
$$ \mathrm{I}{\mathrm{V}}^{2}={\int }_{0}^{\mathrm{\Delta }t}{v}^{2} ( t ) \mathrm{d}t \text{,} $$ (4) 式中,v为地震波速度,Δt为地震波时长,单位为s。
3) 累积绝对速度值CAV。该参数包含了地表加速度的过程信息,最早被伊斯坦布尔地震预警系统用于判别地震大小(Alcik et al,2009),其计算公式为:
$$ \mathrm{C}\mathrm{A}\mathrm{V}={\int }_{0}^{\mathrm{\Delta }t}\left|a ( t ) \right|\mathrm{d}t \text{,} $$ (5) 式中a (t)为t时刻的地震波加速度。
4) 累积绝对位移值CAD。该参数包含了地表位移的过程信息,马美帅等(2022)首先提出利用该参数估算震级,其计算公式为:
$$ \mathrm{C}\mathrm{A}\mathrm{D}=\sum_{i}\left|{\int }_{0}^{\mathrm{\Delta }t}{\int }_{0}^{\mathrm{\Delta }t}a ( t ) \mathrm{d}t\mathrm{d}t\right| \text{.} $$ (6) 时频类特征参数在有无震源距时都可以估算震级,估算震级的依据是地震波的时域和频域都包含了与震级相关的信息,综合利于两方面信息更有利于震级估算。Wang等(2022)将CAD和${\tau ^{\max }_{\mathrm{p}}} $相乘得到的时频类特征参数Sdτ用于震级估算,Sdτ的计算公式为
$$ {S}_{\mathrm{d}\tau }=\sum _{i}\left|{D}_{i}\right|{\tau }_{\mathrm{p}}^{\mathrm{m}\mathrm{a}\mathrm{x}} \text{.} $$ (7) 3. 结果及分析
3.1 GPR的训练
由于GPR方法中协方差函数在有限输入点集上要求是正定的,且是一个满足Mercer条件的对称函数,故协方差函数等价于核函数(何志昆等,2013)。常见的协方差函数有径向基核函数、马顿(Matérn)核函数、指数核函数、二次有理核函数等(Seeger,2004;孙斌等,2012)。为选择合适的协方差核函数,将均方误差(mean squared error,缩写为MSE)$ [ $式(8)$ ] $作为损失函数。在初至地震波 3 s 时对样本进行训练分析,径向基核函数、马顿核函数、指数核函数和二次有理核函数的MSE分别为0.114 5,0.126 7,0.118 4和0.110 3,其中二次有理核函数的MSE小于其它协方差函数,故本文采用二次有理协方差函数$ [ $式(9)$ ] $训练GPR方法。
$$ \mathrm{M}\mathrm{S}\mathrm{E}=\frac{1}{n}\sum _{i=1}^{n}{ ( {y}_{{l}}-{y}_{i} ) }^{2} \text{,} $$ (8) 式中,n为样本数,yl为估算震级,yi为实际震级。
$$ {{\boldsymbol{K}}}_{\mathrm{R}\mathrm{Q}} ( {{\boldsymbol{x}}}^{i}\text{,} {{\boldsymbol{x}}}^{j} ) ={{\boldsymbol{\sigma }}}_{{\mathrm{f}}}^{2}{\left[1 + \frac{{ ( {{\boldsymbol{x}}}^{i}-{{\boldsymbol{x}}}^{j} ) }^{{\mathrm{T}}}{\boldsymbol{M}} ( {{\boldsymbol{x}}}^{i}-{{\boldsymbol{x}}}^{j} ) }{2\alpha }\right]}^{-\mathrm{\alpha }} \text{,} $$ (9) 式中:M为l −2的对角矩阵,l为关联性测度超参数;$ {\boldsymbol{\sigma}} _{\text{f}}^2 $为核函数的信号方差;α为核函数的形状参数。参数集合$ {\boldsymbol{\theta }}= \left\{ {M\text{,} {\boldsymbol{\sigma}} _{\text{f}}^2 \text{,} {\boldsymbol{\sigma}} _{\text{n}}^2} \right\} $为包含所有超参数的向量,通过极大似然法求得。首先建立训练样本条件概率的负对数似然函数$L ( {\boldsymbol{\theta}} ) = - {\log {{p}}} ( {{\boldsymbol{y}} \,|{\boldsymbol{X}} \text{,} {\boldsymbol{\theta}} } ) $,并令其对超参数θ求偏导,然后采用共轭梯度法对偏导数进行最小化以得到超参数的最优解,即:
$$ L ( {\boldsymbol{\theta}} ) =\frac{1}{2}{{\boldsymbol{y}}}^{\mathrm{T}}{{\boldsymbol{K}}}_{{{y}}}^{-1}{\boldsymbol{y}} + \frac{1}{2}\lg\left|{{\boldsymbol{K}}}_{y}\right|+\frac{n}{2}\lg2 \text{π} \text{,} $$ (10) $$ \frac{{\partial } L ( {\boldsymbol{\theta}} ) }{{\partial } {{\boldsymbol{\theta}} }_{i}}=-\frac{1}{2}\mathrm{t}\mathrm{r}\left({{\boldsymbol{K}}}^{-1}\frac{{\partial } {\boldsymbol{K}}}{{\partial } {{\boldsymbol{\theta}} }_{i}}\right)-\frac{1}{2}{{\boldsymbol{y}}}^{{\mathrm{T}}}{{\boldsymbol{K}}}_{y}^{-1}\frac{{\partial } {\boldsymbol{K}}}{{\partial } {{\boldsymbol{\theta}} }_{i}}{{\boldsymbol{K}}}_{y}^{-1}{\boldsymbol{y}} \text{,} $$ (11) 式中,tr ()表示矩阵的迹。
3.2 确定P d及$ {{{{\boldsymbol{\tau}} ^{{\mathbf{max}} }_{\mathbf{p}}}}}$的经验公式
利用训练数据集拟合震级与P d和${\tau ^{\max }_{\mathrm{p}}} $的经验公式,如式(12)所示。选取P d和${\tau ^{\max }_{\mathrm{p}}} $作对比的原因是:${\tau ^{\max }_{\mathrm{p}}} $方法与震级的经验公式中未使用震源距参数,与GPR-M方法的参数选择具有共同因子,且频域参数震级估算方法的对比结果表明${\tau ^{\max }_{\mathrm{p}}} $方法估算震级的准确性高于其它频域参数方法(Nazeri et al,2017;王延伟等,2020),该方法已经应用于地震预警系统 (Hsiao et al,2009;Chung et al,2019);P d方法与震级的经验公式中使用了震源距参数,与GPR-M-R方法的参数选择具有共同因子,且特征参数震级估算方法的对比研究表明P d方法估算震级的准确性明显优于其它参数方法(Leyton et al,2018;王延伟等,2020),该方法也被广泛应用于地震预警系统 (Hsiao et al,2009)。据此本文选择了以上两种特征参数方法进行对比。
$$ \mathrm{lg} ( \mathit{{\mathrm{Para}}} ) {=} \mathit{a} { } \mathit{M} \mathrm + \mathit{b} {\lg } \mathit{R} { +{{ c}}} \text{,} $$ (12) 式中:Para为特征参数(P d或${\tau ^{\max }_{\mathrm{p}}} $);M为震级;R为震源距;a,b,c为待定的回归系数;${\tau ^{\max }_{\mathrm{p}}} $不需要震源距R项时,c=0;P d需要震源距R项时,c ≠ 0。
地震预警系统中的震级估算是随时间持续更新的过程,利用更长的初至地震波可以提高震级估算的准确性(Kanamori,2005;Wu,Zhao,2006;Ziv,2014;Peng et al,2017),初至P波到达3—10 s时,特征参数的经验公式可由训练数据集得到,利用最小二乘法回归得到经验公式(12)的回归系数,如表1所示。
表 1 初至3—10 s地震波的Pd和$ {\tau ^{\max }_{\mathrm{p}}} $方法的拟合系数表Table 1. Fitting coefficients for the Pd and $ {\tau ^{\max }_{\mathrm{p}}} $ methods based on the first 3−10 s seismic waves初至地震波
窗长/s$P_{{\mathrm{d}}} $方法 $ {\tau ^{\max }_{\mathrm{p}}}$方法 a b c a b 3 0.639 4 −3.987 2 −0.840 7 0.309 6 −2.049 4 4 0.681 4 −3.910 9 −0.966 4 0.320 6 −2.048 3 5 0.708 0 −3.740 4 −1.104 4 0.329 8 −2.041 8 6 0.723 6 −3.518 2 −1.248 1 0.339 8 −2.016 5 7 0.738 2 −3.292 8 −1.378 0 0.348 6 −1.975 3 8 0.740 9 −3.070 7 −1.478 3 0.352 6 −1.955 5 9 0.737 8 −2.897 5 −1.538 8 0.358 3 −1.911 7 10 0.733 4 −2.791 4 −1.561 7 0.364 7 −1.875 0 3.3 利用初至3 s P波测试数据集的震级估算结果
在地震预警系统的研究和应用中,普遍认为初至3 s P波估算的震级较好地兼顾了地震预警的时效性和准确性(Kanamori,2005;Wu,Zhao,2006;Peng et al,2017;Wang et al,2022)。因此,改善初至3 s P波估算震级的效果至关重要。图3为GPR-M方法与${\tau ^{\max }_{\mathrm{p}}} $方法估算震级的误差分布及其直方图。从误差的定性分布来看:GPR-M方法估算震级的误差离散性明显小于${\tau ^{\max }_{\mathrm{p}}} $方法且误差更多地分布在$ [ $−0.5, 0.5$ ] $之内,GPR-M方法比${\tau ^{\max }_{\mathrm{p}}} $方法估算震级的误差更小、误差更接近于零。随着震级的增加,两者对MW≥6.5地震事件的震级估算偏小(误差分布下倾),使得MW4.0—6.4范围内的准确率明显高于MW6.5—9.0范围内的准确率。这与众多研究中的震级饱和问题是一致的(Wu,Zhao,2006;王延伟等,2020;Wang et al,2022)。从误差的定量结果看:GPR-M方法的误差标准差、误差均值和误差绝对值均值明显低于${\tau ^{\max }_{\mathrm{p}}} $方法(GPR-M方法分别为0.353 9,−0.091 0和0.274 7,${\tau ^{\max }_{\mathrm{p}}} $方法分别为0.745 5,−0.142 8和0.687 4);GPR-M方法的误差标准差比${\tau ^{\max }_{\mathrm{p}}} $方法降低了约52.53%;在MW4.0—6.4范围内,GPR-M方法震级估算的准确率约是${\tau ^{\max }_{\mathrm{p}}} $方法的1.59倍(GPR-M方法为88.75%,${\tau ^{\max }_{\mathrm{p}}} $方法为55.94%);在MW6.5—9.0范围内,即使两者表现出震级饱和现象,GPR-M方法震级估算的准确率也约是${\tau ^{\max }_{\mathrm{p}}} $方法的1.41倍(GPR-M方法为47.89%,${\tau ^{\max }_{\mathrm{p}}} $方法为33.98%)。
图 3 初至3 s P波时GPR-M (a)和$ {\tau ^{\max }_{\mathrm{p}}}$(b)估算震级的误差分布及其直方图图中的散点表示每条记录估算震级的误差,蓝色实线和红色虚线表示±0.5,σ为误差标准差,μ为误差均值,ω为误差绝对值均值,准确率定义为误差在$ [ $−0.5, 0.5$ ] $的记录数与记录总数的比值,下同Figure 3. Distribution of estimated magnitude errors by GPR-M (a) and $ {\tau ^{\max }_{\mathrm{p}}}$ (b) at initial 3 s P-waveScatter points show magnitude estimation errors based on per record. Blue solid lines and red dashed lines indicate ±0.5 . σ is error standard deviation,μ is mean error,ω is mean absolute error,accuracy is the ratio of records with errors in $ [ $−0.5,0.5$ ] $ to total records,the same below图4为GPR-M-R方法与P d方法估算震级的误差分布及其直方图。从误差的定性分布来看:GPR-M-R方法估算震级的误差更多地分布在$ [ $−0.5, 0.5$ ] $之内,随着震级的增加,对于MW≥6.5的地震事件同样表现出了震级低估现象。从误差的定量结果看:GPR-M-R方法的误差标准差、误差均值和误差绝对值均值都优于P d方法(GPR-M-R方法分别为0.338 4,0.078 1和0.240 7;P d方法分别为0.570 2,−0.094 6和0.444 3);GPR-M-R方法的误差标准差比P d方法降低了约40.65%;在MW4.0—6.4范围内,GPR-M-R方法震级估算的准确率是P d方法的1.45倍(GPR-M-R方法为90.54%,P d方法为62.28%);在MW6.5—9.0范围内,GPR-M-R方法震级估算的准确率是P d方法的1.06倍(GPR-M-R方法为48.67%,P d方法为45.89%)。
3.4 不同初至地震波时长的测试数据集的震级估算结果
将初至地震波时长从3 s增加至10 s,测试GPR-M方法与$ {\tau ^{\max }_{\mathrm{p}}}$方法、GPR-M-R方法与P d方法持续估算震级的效果。图5a为GPR-M方法与$ {\tau ^{\max }_{\mathrm{p}}}$方法所估算震级的误差标准差随初至地震波时长的变化,可见:两种方法的误差标准差随着初至地震波时长的增加逐步减小,但GPR-M方法的误差标准差始终小于$ {\tau ^{\max }_{\mathrm{p}}}$方法,GPR-M方法的误差标准差比$ {\tau ^{\max }_{\mathrm{p}}}$方法降低了52.53%—57.21%。值得注意的是,GPR-M方法在3 s时的误差标准差(0.353 9)仍小于$ {\tau ^{\max }_{\mathrm{p}}}$方法在10 s时的误差标准差(0.668 6)。图5b为GPR-M-R方法与P d方法所估算震级的误差标准差随初至地震波时长的变化,可见:两种方法的误差标准差同样随着初至地震波时长的增加逐步减小,GPR-M-R方法的误差标准差始终小于P d方法,GPR-M-R方法的误差标准差比P d方法降低了约37.72%—40.65%,并且GPR-M-R方法在3 s时的误差标准差(0.338 4)小于P d方法在10 s时的误差标准差(0.404 8)。
图 5 初至3—10 s地震波的不同方法估算震级的误差标准差对比(a) GPR-M与$ {\tau ^{\max }_{\mathrm{p}}}$方法对比;(b) GPR-M与P d方法对比Figure 5. Comparison of the standard deviation of errors in estimating magnitudes for initial 3−10 s seismic waves(a) Comparison of GPR-M method with $ {\tau ^{\max }_{\mathrm{p}}}$ method ;(b) Comparison of GPR-M method with P d method图6和图7给出了GPR-M方法与$ {\tau ^{\max }_{\mathrm{p}}}$方法。GPR-M-R方法与P d方法估算不同范围震级的准确率随初至地震波时长的变化,可见:在MW4.0—6.4范围内,GPR-M方法的准确率由88.75%提高至91.68%,$ {\tau ^{\max }_{\mathrm{p}}}$方法的准确率由55.94%提高至61.23%,GPR-M方法的准确率始终约是$ {\tau ^{\max }_{\mathrm{p}}}$的1.50—1.59倍(图6a);在MW6.5—9.0范围内,GPR-M方法的准确率由47.89%提高至71.45%,$ {\tau ^{\max }_{\mathrm{p}}}$方法的准确率由33.98%提高至46.81%,GPR-M方法的准确率始终约为$ {\tau ^{\max }_{\mathrm{p}}}$方法的1.41—1.53倍(图6b);在MW4.0—6.4范围内,GPR-M-R方法的准确率在由90.54%提高至94.80%,P d方法的准确率由62.28%提高到79.02%,GPR-M-R方法的准确率始终是P d方法的1.20—1.45倍 (图7a);在MW6.5—9.0范围内,GPR-M-R方法的准确率由48.67%提高至78.12%,P d方法的准确率由45.89%提高至68.17%,GPR-M-R方法的准确率约为P d方法的1.06—1.15倍(图7b)。综上表明,GPR-M方法和GPR-M-R方法估算震级准确率的稳定性高于$ {\tau ^{\max }_{\mathrm{p}}}$方法和P d方法。
3.5 泛化数据集的震级估算结果
上述测试中,训练数据集和测试数据集均采用了日本的地震记录,两数据集都包含一些相似的区域信息,例如相似的震源、传播路径、场地条件、监测仪器等,这些区域信息可能会使得GPR-M和GPR-M-R方法仅适用于日本地区,而无法用于其它地区。因此,为了检验日本地震记录训练的GPR-M方法和GPR-M-R在其它地区估算震级的效果,需要测试其泛化能力。泛化数据集选取智利的地震记录,因为智利的地震记录没有用于训练GPR-M方法和GPR-M-R方法,也未被用于拟合$ {\tau ^{\max }_{\mathrm{p}}}$和P d方法的经验公式。考虑到误差随震级的分布可以最为详细地展示震级估算结果,故给出了泛化数据集的震级估算误差分布图(图8,图9),以P波到达后3 s,5 s,8 s,10 s的结果为例。
由图8和图9可以看出GPR-M方法和GPR-M-R方法估算震级的误差分布更集中在$ [ $−0.5, 0.5$ ] $范围内,两者估算震级的误差标准差、误差均值和误差绝对值均值分别优于$ {\tau ^{\max }_{\mathrm{p}}}$方法和P d方法,GPR-M方法的误差标准差比$ {\tau ^{\max }_{\mathrm{p}}}$方法降低了53.08%—55.13%,GPR-M-R方法的误差标准差比P d方法降低了35.88%—36.59%,这表明GPR-M方法和GPR-M-R方法的泛化能力和估算震级的稳定性分别优于$ {\tau ^{\max }_{\mathrm{p}}}$方法和P d方法。
3.6 震例测试
为了检验GPR震级估算方法对于实际震例的可靠性与合理性,选取了我国的三次典型震例进行验证分析。需要注意的是,我国常用震级为MS,对于MW<8.5的地震,MW与MS之间的差异可以忽略不计(Peng et al,2017),因此,本文给出的GPR,$ {\tau ^{\max }_{\mathrm{p}}}$和P d方法都可以直接用于我国地震的震级估算。选取的三次震例分别为:2019年6月22日四川宜宾MS5.4地震,震中位置为(28.4°N,104.8°E),震源深度为10 km;2021年5月21日云南漾濞MS6.4地震,震中位置为(25.7°N,99.9°E),震源深度为8 km;2008年5月12日四川汶川MS8.0地震,震中位置为(31.0°N,103.4°E),震源深度为14 km。三次地震的震中位置及所用的最先触发的四个台站分布如图10所示。在未知震源距和已知震源距两种情况下,将首个监测台站的P波到时记为0时刻,对每个触发台站P波到达后3 s 至 10 s,以1 s间隔持续估算震级。
图11为震源距未知、三次地震最先触发四个台站时,利用$ {\tau ^{\max }_{\mathrm{p}}}$方法和GPR-M方法持续估算震级的结果。可以看出,$ {\tau ^{\max }_{\mathrm{p}}}$方法明显低估了三次地震的震级,并且估算震级偏差明显大于GPR-M方法:对于宜宾MS5.4地震,采用$ {\tau ^{\max }_{\mathrm{p}}}$方法的最大震级偏差为−1.59,采用GPR-M为−0.54;对于漾濞MS6.4地震,采用$ {\tau ^{\max }_{\mathrm{p}}}$方法的最大偏差为−2.34,采用GPR-M为−0.39;对于汶川MS8.0地震,采用$ {\tau ^{\max }_{\mathrm{p}}}$方法的最大偏差为−3.08,采用GPR-M为−1.40。此外,GPR-M估算震级的离散性比$ {\tau ^{\max }_{\mathrm{p}}}$方法更小。
图 11 使用GPR-M方法(左)和$ {\tau ^{\max }_{\mathrm{p}}}$方法(右)持续估算我国2019年宜宾MS5.4 (a)、2021年漾濞MS6.4 (b)和2008年汶川MS8.0 (c)地震震级的结果Figure 11. Results of continuous magnitude estimation for the 2019 Yibin MS5.4 (a),2021 Yangbi MS6.4 (b),and 2008 Wenchuan MS8.0 (c) earthquakes using the GPR-M-R (left panels) and $ {\tau ^{\max }_{\mathrm{p}}}$ (right panels) methods图12为震源距已知、三次地震最先触发四个台站时,用P d方法和GPR-M-R持续估算震级的结果。可以看出,GPR-M-R估算的震级更接近实际震级,且偏差小于P d方法:对于宜宾MS5.4地震,GPR-M-R的最大偏差为+ 0.47,P d方法为−0.56;对于漾濞MS6.4 地震,GPR-M-R的最大偏差为−0.27,P d方法为−0.57;对于汶川MS8.0 地震,GPR-M-R的最大偏差与P d方法相同,均为−1.59,但GPR-M-R估算的震级更接近实际震级,且离散性更小。
4. 讨论与结论
为改善地震预警系统中震级估算的准确性,本文提出了能够同时适用于有震源距和无震源距两种情况下进行震级估算的GPR震级估算方法。该方法以时域、频域以及时频域多个特征参数作为输入(GPR-M),或以特征参数和震源距作为输入(GPR-M-R),以震级作为输出,通过融合初至地震波中多方面信息实现震级估算。GPR-M方法在未使用震源距的情况下提高了震级估算的准确性,GPR-M-R方法在使用震源距的情况下进一步提高震级估算的准确性。在这两种情况下,利用日本的强震记录进行训练和测试,利用智利的强震记录进行泛化能力测试,利用我国的三次震例进行可靠性和合理性检验,并与当前广泛采用的$ {\tau ^{\max }_{\mathrm{p}}}$方法和P d方法进行了对比,结论如下:
1) 利用日本测试数据集初至3—10 s地震波估算震级时,GPR-M方法和GPR-M-R方法估算震级效果优于$ {\tau ^{\max }_{\mathrm{p}}}$方法和P d方法,GPR-M方法估算震级的误差标准差比$ {\tau ^{\max }_{\mathrm{p}}}$方法降低了52.53%—61.20%,GPR-M-R方法比P d方法降低了37.72%—41.21%。
2) 对于震级较大的地震(MW≥6.5),GPR-M方法和GPR-M-R方法的震级饱和现象无$ {\tau ^{\max }_{\mathrm{p}}}$方法和Pd方法明显,GPR-M方法估算震级的准确率比$ {\tau ^{\max }_{\mathrm{p}}}$方法提高了1.4—1.5倍,GPR-M-R方法估算震级的准确率比P d方法提高了1.2—1.45倍。
3) 基于智利强震记录的泛化能力测试结果表明,GPR-M方法和GPR-M-R方法震级估算的准确性相较于$ {\tau ^{\max }_{\mathrm{p}}}$方法和P d方法更高,泛化能力更强,受区域特性影响更小。
4) GPR-M方法和GPR-M-R方法较好地估算了我国三次地震的震级,且准确度高于$ {\tau ^{\max }_{\mathrm{p}}}$方法和P d方法,表明基于GPR的震级估算方法应用于我国地震预警系统具有较好的可靠性与合理性。
本文所提出的GPR震级估算方法在测试数据集和泛化数据集中表现出了较好的震级估算效果,有效地提高了震级估算的准确性。但对于震级较大的地震,震级饱和现象依然存在。震级饱和问题已经成为了震级估算研究中的学术难题,因为大地震的断层破裂过程有几十秒至上百秒,且破裂过程十分复杂,是否可以使用初始几秒的断层破裂信息来确定最终破裂的规模在理论上仍无统一结论(Olson,Allen,2005;Ide,2019)。此外,本文在进行训练GPR方法时,训练数据集里中小地震记录(4.0≤MW≤6.4,2万2 207条)的占比远大于大地震(6.5≤MW≤ 9.0,1 381条),数据的不均匀性是否影响了GPR-M方法的震级估算效果,也需要进一步研究。尽管GPR方法在智利地震和我国震例的测试中表现出良好的震级估算效果,但区域性的差异(如传播介质和场地条件等)对GPR方法的影响程度尚待明确。
日本防灾科学技术研究所、智利SIBER-RISK数据库和中国地震局工程力学研究所强震动观测中心为本研究提供数据支持,审稿专家为本文提出了诸多宝贵的意见,作者在此一并表示感谢。
-
图 1 2019年1—3月蒙城台(a)、温泉台(b)和密山台(c)同场地多方位之间地电场相关系数δxy及其稳定性Δδxy、最大地磁日K指数的变化以及地电场电极布设示意图
Figure 1. Correlation coefficient δxy between multiple azimuths and their stabilities Δδxy,maximum daily K-index of geomagnetic field as well as elelctrodes layout for the sites Mengcheng (a),Wenquan (b) and Mishan (c) during the period from January to March,2019
图 2 2019年1—3月蒙城台(a)、温泉台(b)和密山台(c)同场地长、短极距之间的相关系数δxx及其稳定性Δδxx和最大地磁日K指数的变化
Figure 2. Correlation coefficient δxx between the long and short pole distances of the same site and their stabilities Δδxx as well as maximum daily K-index of geomagnetic field for the sites Mengcheng (a),Wenquan (b) and Mishan (c) in the period from January to March,2019
图 3 地电场非均匀变化对相关系数的影响
(a) 原始序列y1x、y2x和序列y1x’、y2x及其相关系数;(b) 序列y1x和y2x以任意比例整体变化;(c) 原始波形数据、双向地铁干扰数据及原始相关系数;(d) 双向、NS向、EW向叠加地铁干扰后的相关系数
Figure 3. Effect of correlation coefficient non-uniform variation of geoelectric field
(a) Original series y1x,y2x and series y1x’,y2x and their correlation coefficient;(b) Series y1x and y2x overall change at any ratio;(c) Raw and two-direction subway interference waveform data as well as the original correlation coefficient; (d) Correlation coefficient after stacking subway interference in two-direction,NS and EW
图 4 青海门源MS6.4地震周边台站分布 (a)及松山台和黄羊川台附近场地不同方向间地电场相关系数δxy (b)和长短极距间相关系数δxx (c)的变化
Figure 4. Distribution of stations around the Menyuan MS6.4 earthquake in Qinghai Province and the correlation coefficient δxy between different azimuths (b) and δxx between long and short pole distance (c) for the sites Songshan and Huangyangchuan
图 5 呼图壁MS6.2和精河MS6.6地震周围台站分布(a)及红浅台和温泉台附近场地不同方向间地电场相关系数(b)和长短极距间相关系数(c)
Figure 5. Distribution of stations around the Hutubi MS6.2 and Jinghe MS6.6 earthquakes (a) and the correlation coefficient between different azimuths (b) and those between long and short pole distance (c) for the sites Hongqian and Wenquan
图 7 九寨沟MS7.0地震周边台站分布与相关系数及方位角变化
(a) 震中及台站分布;(b) 方位角α变化;(c) δxy和δxx变化;(d) ΔESP变化
Figure 7. Distribution of stations around the Jiuzhaigou MS7.0 earthquake and variation of correlation coefficient and azimuth
(a) Epicenter and station distribution;(b) Variation of azimuth α;(c) Variation of δxy and δxx;(d) Variation of ΔESP
表 1 2019年1月1日至3月31日同场地多方位之间地电场相关系数δxy及其稳定性Δδxy统计
Table 1 Multi-azimuth and multi-pole distance correlation coefficient δxy and their stabilities Δδxy from 1 January 2019 to 31 March 2019
台站 δ12 δ13 δ23 Δδ12 Δδ13 Δδ23 台站 δ12 δ13 δ23 Δδ12 Δδ13 Δδ23 德都 0.95 0.95 0.95 0.05 0.05 0.05 绥化 0.10 0.90 0.60 0.50 0.10 0.40 密山 −0.90 0.95 −0.90 0.05 0.05 0.05 通河 −0.90 −0.85 0.95 0.10 0.15 0.05 望奎 0.45 0.35 −0.40 0.35 0.35 0.20 肇东 0.80 −0.20 −0.70 0.10 0.20 0.10 林甸 0.80 0.55 0.10 0.20 0.20 0.30 嘉山 −0.95 0.95 −0.95 0.05 0.05 0.05 蒙城 −0.20 0.65 −0.85 0.20 0.15 0.15 高邮 0.45 0.70 0.90 0.15 0.10 0.10 海安 −0.10 0.60 0.70 0.20 0.20 0.20 兰州 0.55 0.80 0.10 0.25 0.20 0.50 山丹 0.90 0.95 0.90 0.10 0.05 0.10 平凉 0.95 −0.80 −0.95 0.05 0.10 0.05 瓜州 0.80 0.95 0.95 0.10 0.05 0.05 寺滩 −0.40 0 0.85 0.30 0.60 0.15 乌什 −0.80 0.60 −0.05 0.20 0.20 0.15 和田 0 0.25 −0.25 0.60 0.15 0.15 温泉 0.90 0.95 0.90 0.10 0.02 0.10 红浅 0.30 0.80 0.80 0.30 0.10 0.10 注:δ变动范围取Δδxy正常变化的最大绝对值。 表 2 2019年1月1日至3月31日同场地长、短极距之间的相关系数δxx及其稳定性Δδxx统计
Table 2 Correlation coefficients δxx between the long and short pole distances and their stabilities Δδxx for the same sites from 1 January 2019 to 31 March 2019
台站 δ11 δ22 δ33 Δδ11 Δδ22 Δδ33 台站 δ11 δ22 δ33 Δδ11 Δδ22 Δδ33 德都 0.90 0.90 0.95 0.10 0.10 0.05 绥化 0.30 0.70 0.50 0.70 0.30 0.50 密山 0.95 0.95 0.95 0.05 0.05 0.05 通河 0.95 0.95 0.95 0.05 0.05 0.05 望奎 0.50 0.95 0.80 0.50 0.05 0.20 肇东 0.90 0.90 0.90 0.10 0.10 0.10 林甸 0.95 0.95 0.95 0.05 0.05 0.05 嘉山 0.90 0.95 0.95 0.10 0.05 0.05 蒙城 0.90 0.90 0.95 0.10 0.10 0.05 高邮 0.90 0.95 0.90 0.10 0.05 0.10 海安 0.90 0.90 0.95 0.10 0.10 0.05 兰州 0.90 0.40 0.60 0.10 0.40 0.40 山丹 0.95 0.90 0.95 0.05 0.10 0.05 平凉 0.90 0.95 0.90 0.10 0.05 0.10 瓜州 0.95 0.95 0.95 0.05 0.05 0.05 寺滩 −0.20 0.90 0.70 0.20 0.10 0.30 乌什 0.95 0.80 0.80 0.05 0.10 0.20 和田 0.80 0.40 0.60 0.20 0.40 0.40 温泉 0.95 0.95 0.95 0.05 0.05 0.05 红浅 0.70 0.55 0.80 0.20 0.35 0.20 注:δ变动范围取Δδxx正常变化的最大绝对值。 表 3 距门源MS6.4地震震中900 km范围内台站的地电场相关系数统计
Table 3 Statistics on correlation coefficient of geoelectrical field for the stations 900 km away from the epicenter of the Menyuan MS6.4 earthquake
台站 震中距
/km场地异常变化
距发震时间/d台站 震中距
/km场地异常变化
距发震时间/d台站 震中距
/km场地异常变化
距发震时间/dδxy δxx δxy δxx δxy δxx 武威 78 − − 中卫 319 16 20 汉王 568 25 − 金银滩 98 − − 白水河 374 23 21 凤翔 585 − − 拦隆口 103 25 − 都兰 348 25 − 瓜州 592 25 30 古丰 110 35 − 大武 377 23 23 宝鸡 642 − − 黄羊 132 40 40 嘉峪关 378 − − 乌加河 675 − − 山丹 132 40 32 银川地电 405 25 − 乾陵 687 − − 红沙 147 35 − 银川 417 30 − 周至 705 − − 松山 177 15 15 固原彭 433 − − 成都 778 − − 寺滩 206 − − 天水 455 − − 合阳 814 − − 高台 246 51 51 石嘴山 473 35 − 临汾 885 − − 兰州 266 − − 平凉 501 25 34 注:“−”表示震前未显示出异常变化;黄羊川台2015年6月中旬δxy和δxx的显著下降对应于11月23日青海祁连MS5.2地震(震中距为241 km)。 表 4 距典型震例震中600 km以内台站的相关系数及大地电场优势方位角α异常台站占总台站比例统计
Table 4 Statistics on correlation coefficient and the azimuth α variation of telluric field for the stations within 600 km from typical earthquake epicenters
地震 时段/mo 异常台站所占比例 地震 时段/mo 异常台站所占比例 δxy δxx α δxy δxx α 门源MS6.4 24 68% 36% 40% 松原MS5.7 20 100% 50% 41.7% 呼图壁MS6.2 24 100% 100% 66.7% 九寨沟MS7.0 12 100% 60% 25% 精河MS6.6 24 100% 100% 100% -
陈小斌,赵国泽. 2009. 自动构建大地电磁二维反演的测点中心网格[J]. 地球物理学报,52(6):1564–1572. doi: 10.3969/j.issn.0001-5733.2009.06.018 Chen X B,Zhao G Z. 2009. Automatic construction of a site-centered grid (SCG) for 2D magnetotelluric inversion[J]. Chinese Journal of Geophysics,52(6):1564–1572 (in Chinese).
范莹莹,杜学彬,Zlotnicki J,谭大成,刘君,安张辉,陈军营,郑国磊,解滔. 2010. 汶川MS8.0大震前的电磁现象[J]. 地球物理学报,53(12):2887–2898. doi: 10.3969/j.issn.0001-5733.2010.12.012 Fan Y Y,Du X B,Zlotnicki J,Tan D C,Liu J,An Z H,Chen J Y,Zheng G L,Xie T. 2010. The electromagnetic phenomena before the MS8.0 Wenchuan earthquake[J]. Chinese Journal of Geophysics,53(12):2887–2898 (in Chinese).
黄清华,刘涛. 2006. 新岛台地电场的潮汐响应与地震[J]. 地球物理学报,49(6):1745–1754. doi: 10.3321/j.issn:0001-5733.2006.06.022 Huang Q H,Liu T. 2006. Earthquakes and tide response of geoelectric potential field at the Niijima station[J]. Chinese Journal of Geophysics,49(6):1745–1754 (in Chinese).
黄清华,林玉峰. 2010. 地震电信号选择性数值模拟及可能影响因素[J]. 地球物理学报,53(3):535–543. Huang Q H,Lin Y F. 2010. Numerical simulation of selectivity of seismic electric signal and its possible influences[J]. Chinese Journal of Geophysics,53(3):535–543 (in Chinese).
马钦忠,冯志生,宋治平,赵卫国. 2004. 崇明与南京台震前地电场变化异常分析[J]. 地震学报,26(3):304–312. doi: 10.3321/j.issn:0253-3782.2004.03.009 Ma Q Z,Feng Z S,Song Z P,Zhao W G. 2004. Study on the variation characteristics of the geoelectric field preceding earthquakes[J]. Acta Seismologica Sinica,26(3):304–312 (in Chinese).
毛桐恩,席继楼,王燕琼,尹淑芝. 1999. 地震过程中的大地电场变化特征[J]. 地球物理学报,42(4):520–528. doi: 10.3321/j.issn:0001-5733.1999.04.010 Mao T E,Xi J L,Wang Y Q,Yin S Z. 1999. The variation characteristics of the telluric field in the process of earthquake[J]. Chinese Journal of Geophysics,42(4):520–528 (in Chinese).
钱复业,赵玉林. 2005. 地电场短临预报方法研究[J]. 地震,25(2):33–40. doi: 10.3969/j.issn.1000-3274.2005.02.005 Qian F Y,Zhao Y L. 2005. Study on geoelectric field method for short-term and impending earthquake prediction[J]. Earthquake,25(2):33–40 (in Chinese).
谭大诚,赵家骝,席继楼,杜学彬,徐建明. 2010. 潮汐地电场特征及机理研究[J]. 地球物理学报,53(3):544–555. Tan D C,Zhao J L,Xi J L,Du X B,Xu J M. 2010. A study on feature and mechanism of the tidal geoelectrical field[J]. Chinese Journal of Geophysics,53(3):544–555 (in Chinese).
谭大诚,王兰炜,赵家骝,席继楼,刘大鹏,于华,陈军营. 2011. 潮汐地电场谐波和各向波形的影响要素[J]. 地球物理学报,54(7):1842–1853. doi: 10.3969/j.issn.0001-5733.2011.07.018 Tan D C,Wang L W,Zhao J L,Xi J L,Liu D P,Yu H,Chen J Y. 2011. Influence factors of harmonic waves and directional waveforms of tidal geoelectrical field[J]. Chinese Journal of Geophysics,54(7):1842–1853 (in Chinese).
谭大诚,赵家骝,刘小凤,范莹莹,刘君,陈军营. 2014. 自然电场的区域性变化特征[J]. 地球物理学报,57(5):1588–1598. doi: 10.6038/cjg20140522 Tan D C,Zhao J L,Liu X F,Fan Y Y,Liu J,Chen J Y. 2014. Features of regional variations of the spontaneous field[J]. Chinese Journal of Geophysics,57(5):1588–1598 (in Chinese).
谭大诚,辛建村,王建军,范莹莹,王玮铭. 2019. 大地电场岩体裂隙模型的应用基础与震例解析[J]. 地球物理学报,62(2):558–571. doi: 10.6038/cjg2019L0584 Tan D C,Xin J C,Wang J J,Fan Y Y,Wang W M. 2019. Application foundation and earthquake case analysis of the telluric field rock crack model[J]. Chinese Journal of Geophysics,62(2):558–571 (in Chinese).
汤吉,詹艳,王立凤,董泽义,赵国泽,徐建郎. 2010. 汶川地震强余震的电磁同震效应[J]. 地球物理学报,53(3):526–534. Tang J,Zhan Y,Wang L F,Dong Z Y,Zhao G Z,Xu J L. 2010. Electromagnetic coseismic effect associated with aftershock of Wenchuan MS8.0 earthquake[J]. Chinese Journal of Geophysics,53(3):526–534 (in Chinese).
王书明,底青云,王若,苏晓璐,Mohamed A. 2018. 三维MCSEM利用电磁场分解消除空气波效应[J]. 地球物理学报,61(2):742–749. doi: 10.6038/cjg2018L0033 Wang S M,Di Q Y,Wang R,Su X L,Mohamed A. 2018. Removing airwave effects using electromagnetic fields decomposition for 3D MCSEM[J]. Chinese Journal of Geophysics,61(2):742–749 (in Chinese).
辛建村,谭大诚. 2017. 地电场多测向日变波形相位关联特征[J]. 地震学报,39(4):604–614. doi: 10.11939/jass.2017.04.014 Xin J C,Tan D C. 2017. Phase correlation features of geoelectric field diurnal waveforms in multi-orientation[J]. Acta Seismologica Sinica,39(4):604–614 (in Chinese).
赵国泽,陆建勋. 2003. 利用人工源超低频电磁波监测地震的试验与分析[J]. 中国工程科学,5(10):27–33. doi: 10.3969/j.issn.1009-1742.2003.10.005 Zhao G Z,Lu J X. 2003. Monitoring & analysis of earthquake phenomena by artificial SLF waves[J]. Engineering Science,5(10):27–33 (in Chinese).
Gershenzon N I,Gokhberg M B,Yunga S L. 1993. On the electromagnetic field of an earthquake focus[J]. Phys Earth Planet Inter,77(1/2):13–19.
Ren H X,Wen J,Huang Q H,Chen X F. 2015. Electrokinetic effect combined with surface-charge assumption:A possible generation mechanism of coseismic EM signals[J]. Geophys J Int,200(2):835–848. doi: 10.1093/gji/ggu435
Ren H X,Huang Q H,Chen X F. 2018. Quantitative understanding on the amplitude decay characteristic of the evanescent electromagnetic waves generated by seismoelectric conversion[J]. Pure Appl Geophys,175(8):2853–2879. doi: 10.1007/s00024-018-1823-z
Varotsos P,Alexopoulos K. 1984a. Physical properties of the variations of the electric field of the earth preceding earthquakes,I[J]. Tectonophysics,110(1/2):73–98.
Varotsos P,Alexopoulos K. 1984b. Physical properties of the variations of the electric field of the earth preceding earthquakes. Ⅱ . Determination of epicenter and magnitude[J]. Tectonophysics,110(1/2):99–125.
Zhang P,Zhao G Z,Tang J,Liu F,Sun W H,Han B,Wang L F,Zhang L. 2018. Establishment of the ELF network in Yunnan and electromagnetic precursory monitoring results of the Yangbi MS5.1 earthquake on March 27,2017[J]. Earthquake Research in China,32(2):254–264.