2020, Vol. 63
2. 中国科学院声学研究所声场声信息国家重点实验室, 北京 100190;
3. 哈尔滨工业大学航天科学与力学系, 哈尔滨 150001
2. State Key Laboratory of Acoustics, Institute of Acoustics, Chinese Academy of Sciences, Bejing 100190, China;
3. Department of Astronautics and Mechanics, Harbin Institute of Technology, Harbin 150001, China
地震诱导电磁现象是当前地球科学领域十分关注的前沿问题.观测资料表明了地震发生过程中可以观测到电磁变化，包括提前于地震波到达的早期电磁波信号(Fujinawa and Takahashi, 1998；Okubo et al., 2011)和同步于地震波到达的伴随电磁场(Nagao et al., 2000；汤吉等，2013；Gao et al., 2016).其中早期电磁波携带着震源以及传播路径的有用信息，对地震预警有重要意义.伴随电磁场则携带着地震波以及信号接收区的介质信息(高永新和胡恒山，2009；Huang et al., 2015；Ren et al., 2016)，有助于研究浅层地表的介质特性(Dzieran et al., 2019).研究地震诱导的电磁场，可以使我们从电磁学角度认识对地震孕育的物理过程，对地震减灾有重要意义(黄清华，2004；赵国泽等，2012).
然而目前对地震电磁场产生的物理机制并不清楚，而且对地震电磁耦合特性的认识也不足，从而限制了地震电磁场在地震减灾中的应用.已有研究表明，动电效应是产生地震电磁场的一种可能机制，即流体饱和孔隙介质在外力(如断层滑动和地震波)的作用下，孔隙流体携带净剩离子运动，并产生电磁场的物理现象(Pride，1994；Revil et al., 1999a, b；Jouniaux et al., 2000；Frenkel，2005；Glover et al., 2012).为了定量描述孔隙介质中由动电效应导致的震电耦合波场，Pride(1994)通过严格的数学推导，得到了一套震电耦合方程组，简要地说，该方程组是描述多孔介质弹性波传播的Biot(1956a, b)方程组和描述电磁波传播的Maxwell方程组的耦合.基于Pride方程，前人针对不同模型进行了理论模拟研究(Haartsen and Pride, 1997；Garambois and Dietrich, 2002；Ren et al., 2010；Zheng et al., 2015；Grobbe and Slob, 2016；Guan et al., 2018)，并预测动电效应能够激发两种震电转换信号，一种是同步于地震波到达的伴随电磁场，其只存在于地震波的扰动区域；另一种是地震波到达物性分界面时产生的界面电磁波.
为解释地震诱导电磁场现象，Gao和Hu(2010)给出了全空间模型双力偶点震源激发的地震波场的解析表达式，并对双力偶源激发的电磁场的辐射图案作出了描述.Slob和Muller(2016)对格林函数作了进一步的补充，可用于任意点源所激发的地震波场和电磁场.然而仅有全空间模型的解是不够的，实际地层总是存在边界，例如自由地表.针对水平分层模型，高永新和胡恒山(2009)研究了地震双力偶点源激发的震电波场.在此基础上，Hu和Gao(2011)采用点源叠加的方法实现了对有限断层模型激发震电响应的模拟，并发现震源破裂终止后断层附近存在残余电场这一现象，很好地解释了Nagao等(2000)观测到的现象.Ren等(2010, 2012)也研究了水平分层模型中地震点源和有限断层激发的电磁场，并分析了电磁场的特征.Gao等(2013a, b)从理论上研究了地震早期电磁波的构成和空间分布特征.这些理论模拟表明，地震震源可以激发出两种电磁信号，一种是震源激发的电磁波，另一种是伴随地震波的电磁场.
然而以上研究主要是针对比较规则的模型，而实际地层是比较复杂的，例如存在起伏的地表和不规则的内部界面.这些不规则界面对地震电磁信号的生成和传播会产生影响.前人提出了一些计算震电耦合波场的数值方法，如时域有限元法(Han and Wang, 2001)、频域有限元法(Zyserman等，2010)以及有限差分法(Guan and Hu, 2008；Gao等，2019)，这些研究主要是面向勘探问题.在天然地震尺度(大于100 km)和频率(几赫兹以下)情况下，含不规则界面特别是起伏地表的地层中震电耦合具有何种特征，目前还缺乏充分的认识.因此，本文利用商业有限元软件COMSOL Multiphysics来模拟不规则地层中震源激发的电磁场，通过对不同介质模型进行模拟，研究起伏地表和地下界面对震电波场的激发和传播的影响，为今后地震电磁观测资料的解释提供参考.
1 理论方法COMSOL Multiphysics软件基于有限元理论开发，可用于多物理场耦合的计算.鉴于有限元方法的普遍性，本文将不具体介绍其方法和原理.本节仅阐述本研究所涉及的理论公式的推导以及边界条件的设置.由于地下介质的复杂性，本文只针对一些简单、但是典型的模型开展理论模拟.
1.1 震电耦合控制方程我们基于Pride震电耦合方程组模拟震电耦合波场.遵循认识物理现象从简单到复杂的规律，本文从二维SHTE模式(Haartsen and Pride, 1997)震电波场入手，研究起伏地形对震电波场激发、传播的影响(该模式仅考虑SH模式地震波和TE模式电磁波的耦合).
如图 1，考虑波在x－z平面内传播，其中x轴向右为正，z轴向下为正.当仅考虑SHTE模式时，Pride方程组可以简化为如下两个方程：
|
图 1 空气/多层孔隙介质模型剖面示意图 Fig. 1 Schematic drawing of air/multiple-porous-layer model |
|
(1) |
|
(2) |
其中，uy和Ey为待求量，分别代表固相位移和电场，G为剪切模量，Fy和fy分别为作用在整个孔隙介质和孔隙流体上的体积力密度，Cy为外加电流源，ε为孔隙介质的介电常数，μ为磁导率，ρ=(1－ф)ρs+фρf为孔隙介质的等效密度，ρs为固相基质密度，ρf为孔隙流体密度，η为孔隙流体的黏滞系数.σ、κ和L分别为孔隙介质的动态电导率、动态渗透率和动电耦合系数，它们都是频率的函数，具体的表达式可参见文献(Pride，1994).在低频范围，这三个变量随频率变化不大，可以取其静态值σ0，κ0和
本文采用COMSOL软件中PDE(Partial Differential Equation)接口进行求解，需要将方程(1)和(2)改写为如下亥姆霍兹方程形式：
|
(3) |
对照方程(1)、(2)和(3)，可以得到
|
(4) |
|
(5) |
|
(6) |
|
(7) |
在采用数值方法模拟波动问题时，需要消除人工截断边界带来的反射.与单纯的地震波和电磁波模拟不同，在模拟震电耦合波场时需要在截断边界同时吸收地震波和电磁波.在众多方法中，完全匹配层(PML)广泛应用于弹性波和电磁波的数值模拟，它能使传入吸收层的波随传播距离快速衰减，从而达到对波的吸收.本文使用Zhang和Shen(2010)给出的复频移完全匹配层(CFS-PML)方法.在PML中，将波动方程的空间偏导项∂/∂x和∂/∂z做以下替换：
|
(8) |
|
(9) |
其中sx和sz为复扩展函数，这里取
|
(10) |
|
(11) |
其中dx(x)和dz(z)为衰减因子，βx(x)和βz(z)为比例因子，αx(x)和αz(z)为频移因子.它们的具体形式如下：
|
(12) |
|
(13) |
|
(14) |
d0和β0是dx(x)和βx(x)在PML外边界处的值.α0是αx(x)在PML内边界(即PML与计算区域的界面)处的值.dz(z)、βz(z)和αz(z)有相同的表达形式.L是PML的宽度.pd, pβ, pα在1~4范围内取值，通常取2.这里取pd=pβ=2，pα=1.
|
(15) |
Vs是横波波速，R是平面波垂直入射到PML外边界时的理论反射系数.本文取PML厚度(网格点数)N为40，R=1×10－5.另外取α0=πfc，其中fc是震源时间函数的中心频率，β0=2.
以上PML设置对于地震波能够有较好的吸收.值得说明的是，PML达到一定厚度时才能充分吸收同时又不会产生较大的虚假反射.但是由于同样频率下电磁波波长比地震波的波长大很多，当PML的厚度能够刚好满足地震波吸收时，并不能满足对电磁波的吸收.本文采用Gao等(2019)的方案：(1)只允许PML吸收地震波，即只对方程(1)中的空间偏导项做式(8)和(9)所示的替换；(2)在PML外部增加了一个附加层，使得电磁波自然扩散，传播至附加层外边界时振幅衰减得足够小，反射至计算区域的电磁波可以忽略不计.附加层的网格尺寸可以比计算区域大很多，因此无需过多网格数量就能很好地吸收电磁波.本文中附加层网格尺寸取最小电磁波波长的1/10，设置10层网格.计算结果表明，PML和附加层的组合能够很好地吸收地震波和电磁波.在获得震电波场的频率域响应后，采用快速傅里叶变换得到时间域波形.
2 数值模拟 2.1 正确性验证 2.1.1 与水平分层模型对比我们首先在水平分层模型下验证COMSOL软件求解震电耦合波场的正确性.对于该模型，我们可以采用解析的方法求解点源激发的SHTE模式震电波场(Gao等，2019).
考虑如图 2所示半空间模型，计算区域(Interior modelled area)尺寸为100 km×30 km.其中PMLs是用来吸收地震波的完全匹配层，Additional layers是用来吸收电磁波的附加层.采用的震源为位于地表下10 km强度为1 N的力源Fy，其距坐标原点和右计算边界水平距离都为50 km.接收点R距离源水平距离为30 km，深度为地表下1 m，详细介质参数见表 1中孔隙介质1.力源的时间函数为雷克子波，其频率域函数表达式为
|
图 2 水平分层半空间模型示意图 Fig. 2 Schematic drawing of the half-space model |
|
表 1 孔隙介质的介质参数 Table 1 Properties of the porous layers |
|
(16) |
其中，中心频率fc=0.5 Hz，时间延迟t0=2 s.
图 3分别给出了采用解析方法(实线)与有限元方法(虚线)计算得到的位移场和电场信号.在位移场Uy中可以看到S波信号，其包括震源激发的直达S波及其地表产生的反射S波.电场Ey中可以看到两个信号，其中EM是电磁波信号，包括震源激发的直达电磁波及其在地表的反射电磁波，由于电磁波速度比地震波快很多，在震源激发后几乎瞬间传到观测点；另一个与S波到时相同的是伴随S波电场信号，这些电场信号在后面会有详细的介绍.可以看到，两种方法的结果一致，说明采用COMSOL软件计算震电耦合波场是正确有效的.
|
图 3 半空间模型中解析解(实线)和有限元数值解(虚线)震电响应对比 (a)位移Uy；(b)电场Ey.“S”和“EM”分别代表由源产生的S波和电磁波. Fig. 3 Comparison of analytical (solid lines) and finite element (dashed lines) seismoelectric signals for the half-space model (a) Displacement Uy; (b) Electric field Ey. The symbols "S" and "EM" denote the S and EM waves generated by the source, respectively. |
互易性也是检验数值模拟正确性的另一种有效方法，根据Pride和Haartsen(1996)，格林函数张量之间有如下互易关系：
|
(17) |
式中GuC代表电流源激发的位移场的格林函数，GEF代表力源激发的电场的格林函数，其中上标代表源的种类，下标代表场量.由于频率域中的－iω对应于时间域中的时间一阶导数∂/∂t，故方程(17)表示电流源在r2处产生的位于r1处的速度场与点力源在r1处产生的位于r2处的电场相等.
我们取2.2.2节中图 8所示的山丘模型进行互易性验证.在r2(x=80 km, z=10 km)处添加一个点电流源Cy=1 A·m-2，在r1 (x=20 km, z=10 km)处接收速度场信号Vy；然后在r1处添加一个点力源Fy=1 N，并在r2处接收电场信号Ey.图 4中实线为点力源激发的电场，虚线为电流源激发的速度场.可以发现两者波形和幅度完全一致，再次验证了采用COMSOL软件计算震电耦合波场的有效性.
|
图 8 山丘模型示意图 Fig. 8 Schematic drawing of the hill model |
|
图 4 山丘模型震电波场的互易性验证 Fig. 4 Reciprocity test of the seismoelectric wavefield for the hill model |
本节中，我们先从具有水平地表的半空间模型入手，研究震电波场的特点，然后在此基础上考虑存在起伏地表和断层情况，考察震电波场的特征.
2.2.1 水平地表模型我们先考虑图 2所示的半空间模型，介质参数与2.1.1节所用相同(表 1中孔隙介质1).图 5给出了三个不同时刻位移和电场的波场快照，其中不同颜色色标代表位移和电场的幅度大小，后面图形色标有相同的意义.t=3 s时刻直达S波尚未传播到地表，从位移(图 5a)中可以看到震源激发的直达S波，其波阵面是一个圆环.在电场(图 5d)中可以看到与S波相同形状的圆环波阵面，这是伴随S波的电场信号.除此以外，电场信号中还观察到了不同于位移场中的信号，在图 5d中圆环之外区域还充满了蓝色亮斑，这是由震源激发的电磁波所产生的电场扰动.因为电磁波速度比地震波快得多，其刚被震源激发就“瞬间”传到空间各个位置，并引起了电场扰动.t=6.5 s直达S波刚好入射到地表，从电场(图 5e)中可以发现空气层中存在覆盖范围较大的红色亮斑，这是直达S波在地表处激发的透射电磁波.地层中存在大范围的蓝色亮斑则是S波在地表激发的反射电磁波，由于电磁波速度远比地震波快，因此透射和反射电磁波在激发后几乎瞬间传到空间中的各个点.t=10 s时刻，在位移场(图 5c)中可以看到直达S波在地表产生了反射S波，在电场(图 5f)中则可以看到地层中伴随反射S波的同震电场信号.
|
图 5 半空间模型中位移和电场在t=3 s、t=6.5 s和t=10 s时刻的波场快照 Fig. 5 Snapshots of the solid displacement and electric field for the horizontal half-space model at the times t=3 s, t=6.5 s and t=10 s |
在电场(图 5f)中还可以看到，S波沿着地表传播时，在空气层中同样产生了电场信号.该电场信号与图 5e中直达S波在地表产生的且以“光速”传播的电磁波不同，其是S波沿地表传播时所激发的倏逝电磁波.倏逝电磁波是由于地震波超过临界入射角入射所产生的非均匀电磁波，又称为消逝波或隐失波.Ren等(2016)对固体-孔隙介质分界面产生的倏逝电磁波做了详细的诠释，并指出该电磁波对于流体-孔隙介质分界面(Gao et al., 2017)和空气-孔隙介质分界面也存在.本文的模拟验证了这个结论，可以看到该信号的幅度随着离地表距离的增加迅速减小，关于倏逝电磁波衰减特征的描述可参见Ren等(2018).
图 6显示了一个水平接收阵列上的位移和电场记录.该接收阵列由201个接收点组成，x坐标从0 km至100 km，间隔500 m，z向深度均为1 m.在位移场(图 6a)中可以看到S波信号，在电场记录(图 6b)中可以看到与S波相同到时的同震电场.此外，还可以观察到两个近似为直线的事件.最先到达的是由源激发的直达电磁波EMd，对应于图 5d中圆环之外区域充满的蓝色亮斑，对于距离震源最近的观测点，该电磁波比直达S波早约5 s；第二个事件EM1是直达S波在源正上方地表处产生的界面电磁波信号，对应于图 5e中源正上方地表处的红色亮斑，其到时约为5 s，与直达S波从震源传播到地表花费的时间10 km/2 km·s-1=5 s一致.
|
图 6 半空间模型中201个接收点上的震电波场 (a)位移Uy；(b)地表下1 m处观测的电场Ey；(c)地表以上1 m空气中的电场Ey.接收点x坐标为0~100 km, 间隔500 m. Fig. 6 Seismoelectric fields at an array of 201 receivers for the horizontal half-space model (a) Solid displacement Uy; (b) Electric field Ey received 1 m below the surface; (c) Electric field Ey received 1 m above the surface. The horizontal coordinates of the receivers are from 0 km to 100 km, with interval 500 m. |
由图 5可知，地表上方空气中也存在电磁扰动.为研究空气中电磁信号的特征，考虑位于地表上方空气中布置一排接收点，x坐标仍为0~100 km，间隔500 m，z均为地表上方1 m.得到空气中的电场信号如图 6c所示，对比图 6c和图 6b，可以发现地表上方观测到的电场信号与地下1 m深度处的信号几乎完全相同.空气中不存在震电耦合，但却出现了类似伴随场的信号，这主要是因为地震波经过时产生了倏逝电磁波.图 7为水平源距30 km地表上方和下方1 m两个观测点上的电场信号对比，可以看出两个观测点上的信号完全相同.
|
图 7 地表下方和上方(空气中)1 m处观测的电场对比 Fig. 7 Comparison of the electric fields observed 1 m below (subsurface) and above (air) the surface |
本节我们研究存在起伏地表的模型中震电信号的特点.为了突出地表起伏对震电信号的影响，我们假设地下为均匀半空间介质，其参数取表 1中的孔隙介质1.考虑两种起伏地表情况，一是含有凸起地表的山丘模型，二是含有凹陷地表的凹陷模型，山丘的高度和凹陷的深度相同，它们相对震源的距离也相同.
(a) 山丘模型
图 8为山丘模型图，考虑位于(x=20 km, z=10 km)处的点力源激发.山丘位于模型中间位置，形貌为半个正弦波.其顶点坐标为(x=50 km, z=－4 km)，即顶点高度为4 km.山丘底部宽度为10 km.计算时采用自由三角形网格，山丘内部和附近网格做加密处理.模型中内部计算区域中最大的网格尺寸为500 m，最小100 m.
图 9给出了不同时刻下位移和电场的波场快照.t=3 s以及t=6.5 s时刻位移和电场波场快照与水平地表模型的结果(见图 5a、5b、5d、5e)相似；在t=15 s时刻S波刚到达山丘左侧与地表交点(图 9g)，并产生了散射电磁波(蓝色的亮斑), 散射电磁波是地震波在不规则界面(如拐点)处产生的电磁波，相对于规则界面产生的电场信号更强(即在该时刻山丘左侧与地表交点处的蓝色亮斑比该点左侧的红色亮斑区域更大)；t=19.5 s时刻S波到达山丘右侧面(图 9h)并激发出电磁波，该电磁波在山丘右侧面的左边表现为红色亮斑，右边则为蓝色亮斑，表明界面内外两侧电场具有相反的极性.
|
图 9 山丘模型中位移和电场在t=3 s、t=6.5 s、t=15 s和t=19.5 s时刻下的波场快照 Fig. 9 Snapshots of the solid displacement and electric field for the hill model at the time t=3 s, t=6.5 s, t=15 s and t=19.5 s |
图 10显示了沿山丘左侧一个水平接收阵列上的位移和电场记录.该接收阵列由91个接收点组成，x坐标从0 km至45 km，间隔500 m，z向深度均为1 m.在位移场(图 10a)中可以看到三个事件，其中Sd是直达S波信号，S1、S2是在山丘左右两侧与地表交点处产生的散射地震波.在电场记录(图 10b)中可以看到与这三个S波相同到时的信号，它们是伴随这些S波的同震电场.此外，还可以观察到几个近似为直线的事件.其中，首先出现的是由源激发的直达电磁波EMd，对应于图 9e中圆环之外区域充满的蓝色亮斑；第二个直线事件EM1是直达S波在源正上方地表处产生的电磁波信号，对应于图 9f中源正上方地表处的红色亮斑；第三个直线事件EM2是直达S波在山体左侧地表交点产生的散射电磁波，对应于图 9g中山体左侧地表交点处的蓝色亮斑，其到时约为14 s，与直达S波从震源传播到该交点处的时间一致，即
|
图 10 山丘左侧91个接收点上的震电波场 (a)位移Uy；(b)电场Ey.接收点坐标为(x=[0, 45 km], z=1 m). Fig. 10 Seismoelectric fields at an array of 91 receivers on the left side of the hill (a) Solid displacement Uy; (b) Electric field Ey. The coordinates of the receivers are (x=[0, 45 km], z=1 m). |
图 11显示了沿山丘右侧一个水平接收阵列上观测到的位移和电场记录.该阵列由91个接收点组成，x坐标从55 km至100 km，间隔500 m，z向深度均为1 m.可以发现，在山丘右侧也可以观测到图 10b中所示的电磁波信号，即震源直接激发的电磁波(EMd)以及S波分别在源正上方地表处(EM1)、山丘左侧地表交点处(EM2)和山丘右侧界面(EM3)产生的电磁波.另外值得说明的是，这些电磁波信号比地震波更早到达.例如，对于x=100 km处的观测点，直达电磁波比直达S波早40 s到达，即使是产生于山丘右侧界面的电磁波(EM3)也比S波早20 s到达.
|
图 11 山丘右侧91个接收点上的震电波场 (c)位移Uy；(d)电场Ey.接收点坐标为(x=[55 km, 100 km], z=1 m). Fig. 11 Seismoelectric fields at an array of 91 receivers on the right side of the hill (c) Solid displacement Uy; (d) Electric field Ey.The coordinates of the receivers are (x=[55 km, 100 km], z=1 m). |
(b) 凹陷模型
再来考察如图 12所示的凹陷地表模型，与山丘模型相似，只是地形由凸起变为凹陷.
|
图 12 凹陷模型示意图 Fig. 12 Schematic drawing of the hollow model |
图 13为凹陷地表模型不同时刻下位移和电场的波场快照，可以看到与山丘模型相似的现象：(1)直达S波在源正上方水平地表处激发出电磁波(图 13f)；(2)S波在凹陷左侧界面也产生了电磁波，且该电磁波的电场信号在界面两侧极性相反(图 13g)；(3)在凹陷与地表的右交点处还产生了散射S波，在电场快照中可以观察到伴随这些散射S波的同震电场信号(图 13h).
|
图 13 凹陷模型中位移和电场在t=3 s、t=6.5 s、t=16 s和t=21 s时刻下的波场快照 Fig. 13 Snapshots of the solid displacement and the electric field for the hollow model at the time t=3 s, t=6.5 s, t=16 s and t=21 s |
图 14显示了沿凹陷左侧一个水平接收阵列上观测到的位移和电场记录.在位移场(图 14a)中可以看到三个事件，分别是直达S波信号Sd和凹陷左右两侧与地表交点处产生的散射地震波S1、S2.在电场记录(图 14b)中除了伴随这些S波的同震电场外，还可以看到由源处激发的直达电磁波EMd(对应于图 13e中圆环之外区域充满的蓝色亮斑)和S波在源正上方地表处产生的界面电磁波信号EM1(对应于图 13f中源正上方地表处的红色亮斑).在t=15 s时刻出现了一个幅度较强的电场信号EM2，其是S波入射凹陷左侧界面时产生的界面电磁波，对应于图 13g中凹陷左侧面的红色亮斑和蓝色亮斑.从图 13 h中可以看到S波在凹陷区右侧地表交点处也产生了散射电磁波，但是该电磁波信号比EM1和EM2幅度弱很多，在图 14b中难以分辨出来.
|
图 14 凹陷左侧91个接收点上的震电波场 (a)位移Uy；(b)电场Ey.接收点坐标为(x=[0, 45 km], z=1 m). Fig. 14 Seismoelectric fields at an array of 91 receivers on the left side of the hollow (a) Solid displacement Uy; (b) Electric field Ey. The coordinates of the receivers are (x=[0, 45 km], z=1 m). |
图 15分别给出了山丘模型(实线)和凹陷模型(虚线)在x=40 km观测点处的位移和电场信号.从图 15a位移uy中可以看到三个信号，其中t=10 s到达的Sd信号是直达S波；t=16 s左右到达的S1信号是直达S波到达山丘或凹陷与地表的左交点处激发出的散射S波，凹陷模型的散射S波振幅比山丘模型的大，且两者反向；t=25 s左右到达的S2信号是直达S波到达山丘或凹陷与地表的右交点处激发出的散射S波，且山丘模型中的信号振幅是明显大于凹陷模型的.从图 15b电场Ey中可以看到由源处产生的直达电磁波EMd、直达S波在源正上方地表处产生的界面电磁波EM1以及三个明显的伴随电场信号Sd、S1、S2.此外，还可以观察到t=15 s时刻直达S波入射山丘和凹陷地区左侧界面时产生了电磁波信号EM2，其中凹陷模型产生的电磁波幅度高于山丘模型.对于山丘模型，还可以观测到电磁波信号EM3，其是直达S波到达山丘模型右侧界面时产生的电磁波，这个电磁信号没有出现在凹陷模型中.
|
图 15 山丘与凹陷模型同一接收点处的震电波场比较，接收点坐标为(x=40 km, z=1 m) Fig. 15 Comparison of seismoelectric signals at the same receiver for the hill and hollow models, The coordinates of the receiver are (x=40 km, z=1 m) |
通过对比水平地表模型和起伏地表模型的模拟结果，我们发现，当地表平坦的时候，震源激发的地震波传播至震源上方地表处可以产生电磁波，此外(不考虑地下界面的存在)将不会再产生电磁波(图 5和6).但若地表存在起伏，地震波在遇到不规则地表(如山丘和凹陷)时则会产生电磁波.
2.2.3 断层模型一直以来断层都备受地球物理工作者的关注，研究断层区域震电效应对地震预警有着重要意义.考虑如图 16所示的模型，断层右侧介质相对于左侧有所下沉，地层参数取表 1中的三种介质.源类型和时间函数与水平地表模型中的一致.在地表处等距布置81个接收点用于接收地震波和电磁波信号, x坐标从60 km至100 km，间隔500 m，z向深度均为1 m.
|
图 16 断层模型示意图 Fig. 16 Schematic drawing of the fault model |
图 17为不同时刻的震电波场快照，在位移场图中可以看到震源激发的S波在不同界面处发生了反射和透射，在电场图中可以观察到伴随直达、反射和透射S波的同震信号，此外也可以发现这些地震波在源正上方地表处、源正下方分层处、断层处以及内部拐点处都产生了电磁波信号.
|
图 17 断层模型中位移和电场在t=2 s、t=4.5 s、t=18 s和t=35.5 s时刻下的波场快照 Fig. 17 Snapshots of the solid displacement and the electric field in fault model at the time t=2 s, t=4.5 s, t=18 s and t=35.5 s |
图 18画出了81个接收点上接收的位移和电场记录.在位移记录中可以看到由于地下发生了多次反射、透射，位移记录中存在多个地震波信号.在电场图中除了有伴随这些地震波的电场信号之外，还可以看到在地震波到达之前出现了多个表现为“直线”的电磁波信号，这是由地震波在地下界面多次激发产生的电磁波.图 19给出了(x=80 km，z=1 m)处接收点上的位移场和电场.S波在t=32 s左右到达，其后为多次反射波(图 19a).在电场记录(图 19b)中可以观测到伴随这些地震波的同震信号.另外，在t=0 s出现的是震源激发的电磁波，其后至S波到达前出现的是S波在地下不同界面处激发的电磁波，其幅度相比震源激发的电磁波和伴随地震波的同震信号弱.
|
图 18 断层模型中81个接收点上的震电波场 (a)位移Uy；(b)电场Ey.接收点坐标为(x=[60 km, 100 km], z=1 m). Fig. 18 Seismoelectric fields at an array of 81 receivers for the fault model (a) Solid displacement Uy; (b) Electric field Ey, The coordinates of the receivers are (x=[60 km, 100 km], z=1 m). |
|
图 19 断层模型中接收点(x=80 km, z=1 m)上的震电波场 (a)位移Uy；(b)电场Ey. Fig. 19 Seismoelectric signals at the receiver (x=80 km, z=1 m) for the fault model (a) Solid displacement Uy; (b) Electric field Ey. |
本文采用有限元软件COMSOL模拟了点力源激发的SHTE模式震电波场，通过在水平分层模型下与解析方法获得的结果对比以及互易性测试，我们验证了采用该软件计算震电耦合波场的正确性和有效性，并分别模拟了含起伏地表和断层模型的震电响应，主要结论如下：
(1) 存在三种类型的电磁信号，第一种是震源直接激发的电磁波；第二种是地震波在分界面处激发的电磁波，分界面包括自由表面(即空气-孔隙介质界面，如图 5e)和地下不同介质界面(即孔隙介质-孔隙介质分界面)；第三种是伴随地震波的同震电场信号.震源激发的电磁波比地震波速度快很多，比地震波更早到达观测台站.这预示着实际地震发生时，存在比地震波早到的电磁信号，如果能识别此种电磁信号并应用于地震预警，那么可以大大提高预警时间，对减少地震灾害具有重要意义.
(2) 当地表为不存在起伏的水平界面且不考虑地下界面时，震源激发的地震波仅在震源正上方处地表产生电磁波.当震源和观测点之间存在起伏地表时，地震波在起伏表面上(如图 9h)和拐点处(如图 9g山丘与地表交点处和图 13h凹陷区与地表处交点处)也可以激发出电磁波, 这些电磁波也比地震波更早到达观测点.这预示着，实际地震中可能存在出现于震源发动之后且地震波到达之前的电磁波信号，此种电磁波信号对于地震预警同样具有价值.
(3) 地震波传播至地表并沿着地表传播时，在空气层中同样产生了电场信号，将观测点放置在空气层中接收到的电磁信号与地下接收点接收的结果相同，这是因为地震波以超过临界折射角的角度入射到地表，激发出了倏逝电磁波(Ren et al., 2016)，从图 5f可以观测到该倏逝电磁波幅度远离地表时迅速减小，因此仅在靠近地表的地方具有较强的幅度.Ujihara等(2004)在日本2003年5月26日宫城县M7.1地震发生后，在震中区附近进行电磁监测，并记录到了多个余震中的同震电磁信号.为了验证记录到的同震电磁信号不是由仪器晃动引起的，他们将磁传感器悬挂于空中，结果发现悬挂于空气中的磁场传感器同样记录到了伴随地震波的同震磁场信号.本文的模拟结果为解释Ujihara等(2004)的观测提供了理论证据.
本文模拟时仅考虑二维SHTE模式震电耦合波，所考虑的源也是简单的点力源，但值得说明的是，所得到的结论同样适用于二维PSVTM模式以及三维模型下的震电耦合波场.可以预期的是，地震震源激发的P波和SV波在地表处以及地下分界面处同样会激发出比地震波速度更快的电磁波，对于远离震源的台站而言，这些电磁波信号将比地震波更早到达.后续将针对二维PSVTM模式以及三维模型下的震电耦合波场开展理论模拟，并研究震电波场的特征，为今后地震电磁信号的解释和利用电磁信号服务于地震减灾打下基础.
Biot M A. 1956a. Theory of propagation of elastic waves in a fluid-saturated porous solid. Ⅰ. Low-frequency range. The Journal of the Acoustical Society of America, 28(2): 168-178. DOI:10.1121/1.1908239 |
Biot M A. 1956b. Theory of propagation of elastic waves in a fluid-saturated porous solid. Ⅱ. Higher frequency range. The Journal of the acoustical Society of America, 28(2): 179-191. DOI:10.1121/1.1908241 |
Dzieran L, Thorwart M, Rabbel W, et al. 2019. Quantifying interface responses with seismoelectric spectral ratios. Geophysical Journal International, 217(1): 108-121. DOI:10.1093/gji/ggz010 |
Frenkel J. 2005. On the theory of seismic and seismoelectric phenomena in a moist soil. Journal of Engineering Mechanics, 131(9): 879-887. DOI:10.1061/(ASCE)0733-9399(2005)131:9(879) |
Fujinawa Y, Takahashi K. 1998. Electromagnetic radiations associated with major earthquakes. Physics of the Earth and Planetary Interiors, 105(3-4): 249-259. DOI:10.1016/S0031-9201(97)00117-9 |
Gao Y X, Hu H S. 2009. Numerical simulation and analysis of seismoelectromagnetic wave fields excited by a point source in layered porous media. Chinese Journal of Geophysics (in Chinese), 52(8): 2093-2104. DOI:10.3969/j.issn.0001-5733.2009.08.018 |
Gao Y X, Hu H S. 2010. Seismoelectromagnetic waves radiated by a double couple source in a saturated porous medium. Geophysical Journal International, 181(2): 873-896. |
Gao Y X, Ch en, X F, Hu H S, et al. 2013a. Early electromagnetic waves from earthquake rupturing: Ⅰ.theoretical formulations. Geophysical Journal International, 192(3): 1288-1307. DOI:10.1093/gji/ggs096 |
Gao Y X, Chen X F, Hu H S, et al. 2013b. Early electromagnetic waves from earthquake rupturing:Ⅱ.validation and numerical experiments. Geophysical Journal International, 192(3): 1308-1323. DOI:10.1093/gji/ggs097 |
Gao Y X, Harris J M, Wen J, et al. 2016. Modeling of the coseismic electromagnetic fields observed during the 2004 MW6.0 Parkfield earthquake. Geophysical Research Letters, 43(2): 620-627. DOI:10.1002/2015GL067183 |
Gao Y X, Wang M Q, Hu H S, et al. 2017. Seismoelectric responses to an explosive source in a fluid above a fluid-saturated porous medium. Journal of Geophysical Research:Solid Earth, 122(9): 7190-7218. DOI:10.1002/2016JB013703 |
Gao Y X, Wang D D, Yao C, et al. 2019. Simulation of seismoelectric waves using finite-difference frequency-domain method:2-D SHTE mode. Geophysical Journal International, 216(1): 414-438. |
Garambois S, Dietrich M. 2002. Full waveform numerical simulations of seismoelectromagnetic wave conversions in fluid-saturated stratified porous media. Journal of Geophysical Research:Solid Earth, 107(B7): ESE 5-1-ESE 5-18. DOI:10.1029/2001JB000316 |
Glover P W J, Walker E, Jackson M D. 2012. Streaming-potential coefficient of reservoir rock:A theoretical model. Geophysics, 77(2): D17-D43. |
Grobbe N, Slob E C. 2016. Seismo-electromagnetic thin-bed responses:Natural signal enhancements?. Journal of Geophysical Research:Solid Earth, 121(4): 2460-2479. DOI:10.1002/2015JB012381 |
Guan W, Hu H S. 2008. Finite-difference modeling of the electroseismic logging in a fluid-saturated porous formation. Journal of Computational Physics, 227(11): 5633-5648. DOI:10.1016/j.jcp.2008.02.001 |
Guan W, Shi P, Hu H S. 2018. Contributions of poroelastic-wave potentials to seismoelectromagnetic wavefields and validity of the quasi-static calculation:a view from a borehole model. Geophysical Journal International, 212(1): 458-475. DOI:10.1093/gji/ggx417 |
Haartsen M W, Pride S R. 1997. Electroseismic waves from point sources in layered media. Journal of Geophysical Research:Solid Earth, 102(B11): 24745-24769. DOI:10.1029/97JB02936 |
Han Q Y, Wang Z J. 2001. Time-domain simulation of SH-wave-induced electromagnetic field in heterogeneous porous media:A fast finite-element algorithm. Geophysics, 66(2): 448-461. |
Hu H S, Gao Y X. 2011. Electromagnetic field generated by a finite fault due to electrokinetic effect. Journal of Geophysical Research:Solid Earth, 116(B8): B08302. DOI:10.1029/2010JB007958 |
Huang Q H, Ren H X, Zhang D, et al. 2015. Medium effect on the characteristics of the coupled seismic and electromagnetic signals. Proceedings of the Japan Academy, Series B, 91(1): 17-24. DOI:10.2183/pjab.91.17 |
Jouniaux L, Bernard M L, Zamora M, et al. 2000. Streaming potential in volcanic rocks from Mount Pelée. Journal of Geophysical Research:Solid Earth, 105(B4): 8391-8401. DOI:10.1029/1999JB900435 |
Nagao T, Orihara Y, Yamaguchi T, et al. 2000. Co-seismic geoelectric potential changes observed in Japan. Geophysical Research Letters, 27(10): 1535-1538. DOI:10.1029/1999GL005440 |
Okubo K, Takeuchi N, Utsugi M, et al. 2011. Direct magnetic signals from earthquake rupturing:Iwate-Miyagi earthquake of M7.2, Japan. Earth and Planetary Science Letters, 305(1-2): 65-72. DOI:10.1016/j.epsl.2011.02.042 |
Pride S. 1994. Governing equations for the coupled electromagnetics and acoustics of porous media. Physical Review B, 50(21): 15678-15696. DOI:10.1103/PhysRevB.50.15678 |
Pride S R, Haartsen M W. 1996. Electroseismic wave properties. The Journal of the Acoustical Society of America, 100(3): 1301-1315. DOI:10.1121/1.416018 |
Ren H X, Huang Q H, Chen X F. 2010. A new numerical technique for simulating the coupled seismic and electromagnetic waves in layered porous media. Earthquake Science, 23(2): 167-176. DOI:10.1007/s11589-009-0071-9 |
Ren H X, Chen X F, Huang Q H. 2012. Numerical simulation of coseismic electromagnetic fields associated with seismic waves due to finite faulting in porous media. Geophysical Journal International, 188(3): 925-944. DOI:10.1111/j.1365-246X.2011.05309.x |
Ren H X, Huang Q H, Chen X F. 2016. Numerical simulation of seismo-electromagnetic fields associated with a fault in a porous medium. Geophysical Journal International, 206(1): 205-220. DOI:10.1093/gji/ggw144 |
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. Pure and Applied Geophysics, 175(8): 2853-2879. DOI:10.1007/s00024-018-1823-z |
Revil A, Pezard P A, Glover P W J. 1999a. Streaming potential in porous media:1.Theory of the zeta potential. Journal of Geophysical Research:Solid Earth, 104(B9): 20021-20031. DOI:10.1029/1999JB900089 |
Revil A, Schwaeger, H, Cathles L M, et al. 1999b. Streaming potential in porous media:2.Theory and application to geothermal systems. Journal of Geophysical Research:Solid Earth, 104(B9): 20033-20048. DOI:10.1029/1999JB900090 |
Slob E, Mulder M. 2016. Seismoelectromagnetic homogeneous space Green's functions. Geophysics, 81(4): F27-F40. DOI:10.1190/geo2015-0337.1 |
Tang J, Zhan Y, Wang L F, et al. 2013. Electromagnetic coseismic effect associated with aftershock of Wenchuan MS8.0 earthquake. Chinese Journal of Geophysics (in Chinese), 53(3): 526-534. DOI:10.3969/j.issn.0001-5733.2010.03.006 |
Ujihara N, Honkura Y, Ogawa Y. 2004. Electric and magnetic field variations arising from the seismic dynamo effect for aftershocks of the M7.1 earthquake of 26 May 2003 off Miyagi Prefecture, NE Japan. Earth, Planets and Space, 56(2): 115-123. |
Zhang W, Shen Y. 2010. Unsplit complex frequency-shifted PML implementation using auxiliary differential equations for seismic wave modeling. Geophysics, 75(4): T141-T154. DOI:10.1190/1.3463431 |
Zhao G Z, Wang L F, Zhan Y, et al. 2012. A new electromagnetic technique for earthquake monitoring-CSELF and the first observational network. Seismology and Geology (in Chinese), 34(4): 576-585. |
Zheng X B, Hu H S, Guan W, et al. 2015. Simulation of the borehole quasistatic electric field excited by the acoustic wave during logging while drilling due to electrokinetic effect. Geophysics, 80(5): D417-D427. DOI:10.1190/geo2014-0506.1 |
Zyserman F I, Gauzellino P M, Santos J E. 2010. Finite element modeling of SHTE and PSVTM electroseismics. Journal of Applied Geophysics, 72(2): 79-91. |
高永新, 胡恒山. 2009. 水平分层孔隙介质中点源激发的震电波场数值模拟及分析. 地球物理学报, 52(8): 2093-2104. DOI:10.3969/j.issn.0001-5733.2009.08.018 |
黄清华. 2004. 电磁学手法在地震研究中的应用. 石油地球物理勘探, 39(S1): 75-79, 84. |
汤吉, 詹艳, 王立凤, 等. 2013. 汶川地震强余震的电磁同震效应. 地球物理学报, 53(3): 526-534. DOI:10.3969/j.issn.0001-5733.2010.03.006 |
赵国泽, 王立凤, 詹艳, 等. 2012. 地震预测人工源极低频电磁新技术(CSELF)和第一个观测台网. 地震地质, 34(4): 576-585. DOI:10.3969/j.issn.0253-4967.2012.04.004 |