于子叶

(中国北京100081中国地震局地球物理研究所)



层状多孔介质的地震电磁异常数值模拟方法研究*

于子叶*

(中国北京100081中国地震局地球物理研究所)

本文针对电激发模型在各向同性分层介质中利用传播矩阵方法对电磁异常产生过程的数值模拟． 层状模型中一部分为饱和多孔介质， 另一部分为普通线弹性介质， 数值模拟过程中将固体场作为主要场单独进行计算， 而将其形变作为等效流体源计算耦合电磁场． 模拟结果显示， 地震波在经过饱和多孔介质边界的时候会产生电磁场， 而其在流体层内部传播的过程中则未产生电磁异常． 模拟结果可以用于解释一些伴随地震产生的电磁异常， 包括地震发生时以及地震波到达时所产生的电磁异常现象， 另一方面对于前震及其它一些无感地震所产生的电磁异常的解释也具有一定的参考价值．

电激发效应 电磁异常 饱和多孔介质 电磁场模拟

引言

大量观测研究表明， 在地震发生过程中会伴随产生电磁异常(Garambois， Dietrich， 2001; Honkuraetal， 2002; Matsushimaetal， 2002; 汤吉等， 2008， 2010； Okuboetal， 2011； Frenkel， 2014)． Honkura等(2002)的研究显示地震波到达之前会产生电磁异常信号， 而Matsushima等(2002)则记录到电磁异常伴随地震波产生; Honkura等(2002)则认为地震发生时检测到的电磁异常产生于震源， 也就是说， 由于地震破裂而产生了一些电磁场扰动， 而与地震波同时到达的电磁异常则与当地场地条件有关． 为了解释这些异常发展了一些电磁场耦合理论， 例如电激发理论(Pride， 1994)、 运动电效应理论(Matsushimaetal， 2002)和压磁效应(Okuboetal， 2011)等． 迄今为止， 地震电磁异常方面的研究多集中于观测层面， 相关理论及数值模拟方法仍不甚成熟， 但是， 地震电磁异常对于地震预警乃至地震前兆异常研究起着非常积极的作用(安张辉等， 2013)， 因此一直是地震学研究的热点问题之一．

介质中由于流体作用产生电磁场的机制是颇有前景的研究方向． 岩石实验研究已表明， 地震波通过含水孔隙介质时可产生电磁信号(Zhuetal， 2000； Thompson， Gist， 2012)． 在理论模拟方面， 考虑到层状介质中格林函数的模拟方法目前已相对成熟， 而且Chen(1999)发展的反透射系数法同样可以高效地计算层状介质中的地震波， 这为层状介质中电磁异常的模拟计算提供了重要的参考． Pride(1994)具体阐述了电激发原理并给出了约束方程． 基于电激发原理产生电磁波的过程大致为： 固体场运动使其中的流体产生相对运动(Kennett， 1984)， 而固体颗粒由于物理或化学作用在其上附着电荷， 同时为保持整体的电中性， 流体中携带有相反电荷， 流体和固体带电层的相对运动在宏观上产生了电磁异常． 随后， Haartsen和Pride(1997)在约束方程的基础上给出了全空间电磁格林函数的表达式； Gao和Hu(2010)以及Ren等(2010)则运用有限差分方法模拟了异常的产生过程； 同时电激发理论被一些实验室观测(Garambois， Dietrich， 2001； Zysermanetal， 2010)所证实． Gao等(2013)根据传播矩阵方法改进了积分方式， 提出了支点积分方法． 上述数值模拟所选取的模型均为孔隙介质模型(张丹等， 2013)， 并且震源对于电磁信号的产生具有一定的影响(Hu， Gao， 2011)， 然而由于受数值模拟方法所限， 对于地震波通过含水层的情况则少有涉及． 鉴于此， 本文拟采用离散波数的方式对电磁异常信号的产生过程进行数值模拟， 以期能够计算不同地层条件下的电磁异常．

本文的模拟实验基于Pride(1994)的电激发理论， 不同于前人所使用的将流体场和电磁场作为总体的传播矩阵进行模拟， 而是在模拟过程中将固体场作为非耦合部分独立进行模拟计算， 并将已得固体场形变作为流体运动的力源． 震源与孔隙介质层分离之后， 震源对于整个电磁异常模拟的作用仅限于产生一个能够影响到固体场扰动的变量， 这样在计算过程中即可对一些更复杂的含水地层所产生的电磁异常信号进行数值模拟．

1 方法原理

1.1 约束方程及固体场计算

图1b给出了加入含水层的孔隙介质模型示意图， 在该图中建立柱坐标系， 原点位于地表． 为了使计算简便， 使z轴通过点源， 向下为正， 第j层边界记为zj， 地面处记为z0． 根据Pride(1994)提出的电磁异常约束方程， 在忽略流体场和电磁场的情况下， 介质条件简化为层状介质(图1a)， 此时约束方程(Haartsen， Pride， 1994， 1997)可简化为

图1 普通层状介质模型(a)和加入含水层的孔隙介质模型(b)Fig.1 Common stratified media model (a) and stratified media model with fluid-saturated porous layer (b)

(1)

式中，H和G分别为体变模量和剪切模量，ρ为介质密度， I为单位矩阵， u为固体场位移， ω为频率， F为震源， 震源存在于s层中， δjs表示仅在s层中存在震源．

引入矢量柱面谐函数作为基函数， 其表达式为

(2)

对固体场向量展开可得(Chen， 1999)

(3)

式中， V为水平和垂直方向上位移和应力的展开系数，M为系数矩阵， 依次为

(4)

式(4)的解为

(5)

式中，j为地层编号，Λ为对角矩阵．

结合层间边界条件和自由边界条件即可求得每层的固体场方程系数aj(Chen， 1999).

1.2 耦合电磁场计算

将流体电磁耦合场部分方程(Pride， 1994)整理成两个关于固体场u， 流体场w和电场E的线性微分方程：

(6)

其中： ρd=iη/(ωκ)， κ为动态渗透率； εd=ε+i/(ωσ)-ρdL2， C为体应变模量， ρf为液体密度， L为低频电激发耦合因子．

由于不存在流体源和电场源， 式(6)中的[C+ωρfI]·u可以作为产生流体运动和电磁场的源． 耦合场的边界条件为

(7)

接下来， 利用基函数式(2)分别对电场和流体场进行展开， 可得

(8)

电磁场计算时已经将u作为已知源， 即固体场波数和式(5)的系数aj是已知的， 将式(8)带入式(6)整理可得

(9)

取

式(9)是τe的一个关于深度z的四次微分方程， 解为

(10)

其中：

Aj是式(14)所对应的四次齐次微分方程通解的系数． 求得τe后， 可以推导出用其表示的电场波数为

(11)

由电场边界条件及流体场边界条件式(7)得式(9)的边界条件为

(12)

利用边界条件求解可得到τe的4个系数， 至此流体层中产生电磁场的计算完成．

流体层内的系数可根据上述方程直接给出， 而流体层外的数据则需根据电磁场边界条件继续计算. 在不存在流体作用和自由电荷影响的情况下， 假设其它层为各项同性均匀介质， 电场方程则可简化为

(13)

其中

式中σ为电导率，μ为磁导率，ε为介质的介电常数， 其边界条件为

(14)

利用式(8)中电场等式展开后， 带入麦克斯韦方程可得无流体时电场系数方程为

(15)

求解得

(16)

式中： C为方程系数， 通过边界条件即可确定； η为式(15)所对应二次方程的根．

至此， 固体场、 流体场、 电磁场波数域展开系数已全部求出， 对其反变换后可得其频率域的空间解为

(17)

图2 电磁异常模拟程序结构框图Fig.2 Layout of EM anomaly simulation program

式中，θ为观测点坐标与x轴的夹角． 固体场与电场矢量反变换形式相同， 均用v表示， 对v进行傅里叶反变换得到其时间域的表达式为

(18)

至此频率域矢量求解完毕．

1.3 程序结构

本文的电磁异常数值计算过程如图2所示． 首先， 利用类传播矩阵方法计算固体场波数域系数， 结合所读取的含水层参数计算电磁场波数域系数； 然后， 根据所得到的电磁场进行反变换求得频率域解； 最后， 进行傅里叶反变换求得其时间域解．

由于Λj存在一个极小的数， 若在系数a的求解过程中将其当成一个整体矩阵求解， 则很容易出现病态矩阵， 故本文采用Chen(1999)的计算方式， 即将整体矩阵分解成小的传播矩阵， 从而降低矩阵的病态性， 但这样使得式(5)中系数a的计算过程相对复杂了一些．

模拟程序采用标准C(C99标准)进行编写， 利用Window平台下的MinGW编译器(GCC版本4.7.0)进行编译， 便于该程序移植到Linux平台．

2 模拟计算与结果比较

2.1 震源模型

本文在研究过程中着重于固体场扰动产生电磁波的过程， 所以震源的作用在于产生地震波， 从而能影响到相关的流体场， 而震源本身的形式并不重要， 因此在模拟中采用相对简单的矩张量点源模型， 震源深度为15 km， 其表达式为

(19)

其中，

(20)

式中：M0=1.5×108N·m， 相当于M6.0地震； 震源时间函数s(t)为地震子波， 表达式为

(21)

式中，Tc=4 s，f0=1 Hz.

2.2 地层模型

模拟的地层为各向同性介质， 共4层． 本文着重研究地震波传播经过层状介质时所产生的电磁波， 为了避免地层间的多次反射波干扰， 将地层参数设为同一数值． 所选取的弹性参数均为典型岩石的弹性参数． 其中2号地层为含水层， 其参数列于表1． 表中孔隙度和盐浓度的选取参考了Gao等(2013)的模拟参数， 电导率σ0和低频电激发耦合因子L0为盐浓度C0与电势ζ的函数， 其中电势ζ为饱和的石英介质25℃时在pH=7的条件下所得的实验值(Haartsen， Pride， 1997)．

表1 含水层参数

2.3 数值结果

将接收点设在含水层界面上， 该点距震中10 km， 方位角为30°， 测定地震波振幅和电场强度． 观测点之所以这样选取， 着重于观测地震波在入射和出射时所产生的电磁异常， 计算结果如图3和4所示．

数值模拟所得图形采用归一化表示， 再乘以相应的系数即所得电场大小． 从图3与图4的对比可以看出: 初始接收到的电场与固体波场的到时相差约1.7 s， 这是由于地震波与电磁波的速度存在差异所致; 地震波进入含水层下界面时产生了电磁场， 二者从下层界面同时出发， 但电磁场提前到达； 电磁场信号与地震波信号呈现出很强的相关性， 而且仅在地震波入射和出射含水层时产生电磁异常； 地震波振幅与所产生的相应的电场振幅呈正相关， 且电磁波较地震波衰减得快．

图3 含水层上界面固体场位移

图4 含水层上界面电场

2.4 观测结果比较

现有很多研究均报道直接观测到了同震电磁信号， 其中比较有代表性的是2013年甘肃岷县漳县MS6.6地震的同震电磁信号(杨皓等， 2015)． 该地震中， 电场异常与地震波同时到达， 且形状相似． 当P波到达时， 电场的水平和垂直方向均产生了明显的电磁异常， 而之后S波到达时所产生的电磁异常幅度较P波所产生的要大， 这与地震波振幅相关． 假设接收点以下介质均为饱和含水介质， 那么观测结果与数值模拟得到的结果较为一致， 即在地震波入射和出射含水层界面时产生电磁异常， 且其异常幅度与地震波振幅相关； 但电场曲线与地震动曲线相比缺少高频特征， 这说明在电磁场与波场的耦合过程中高频部分耦合得不好， 可能与流体性质有关， 在高频运动过程中由于流动性等因素导致流体场与固体场的相对运动减弱， 使得所产生的电磁场减弱， 这也说明Pride(1994)的约束方程可能缺少一个高频阻尼项． 另外， 未观测到与地震破裂过程伴随产生的电磁异常， 这可能是电磁场随深度衰减过快所致， 也可能是在震源破裂过程中并未产生电磁波．

3 讨论与结论

由式(11)可知， 电场波数完全取决于固体场波数域展开系数， 而固体场的波数解则取决于震源， 即电磁场的解是震源的函数， 所以震源函数与电磁场信号存在一定的相关性， 这说明地电磁异常中常用电磁信号与地震信号进行互相关计算是合理的． 式(11)转换成通用形式如下：

(22)

式中，M和N为系数． 可以看出， 电磁场与地震波场存在着一定的相关性． 整理后，Mi中包含了ω0,ω1,ω2和ω3项. 由傅里叶变换的导数性质可知， 电磁场还与固体场的时间域导数相关．

由模拟结果可以看到两个独立的信号， 这说明电磁场仅在震动通过电性变化区域时才会产生， 这也说明如果在实际中由电激发效应产生电磁场， 必要条件是介质孔隙层间流体条件不一致． 在前震、 地震成核阶段以及无感地震中， 均会有一定程度的地下岩石局部形变， 若该形变场通过孔隙条件复杂的流体层， 那么由数值模拟结果可知该过程会产生电磁场， 这就为解释一些同震电磁异常以及无感地震中观测到的电磁场异常提供了一种可能的解释． 鉴于地下条件的复杂性， 电磁异常的产生过程可能并非如Pride(1994)电激发理论中所描述的那样简单， 本研究对地震波场、 流体场和电磁场等3场耦合产生电磁场的过程进行了有现实意义的探索．

本文仅计算了最强的一次场， 并未考虑电磁场在地层中的反射以及不同频率电磁波的衰减， 希望在以后研究中予以改进． 另外， 本文仅用简化的点源进行了解释和模拟， 而未用实际震源进行更具体的模拟， 进一步研究中需要带入真实的震源模型， 并与实际地震进行对比．

安张辉， 杜学彬， 谭大诚， 范莹莹， 刘君， 崔腾发. 2013. 四川芦山MS7.0和汶川MS8.0地震前地电场变化研究[J]. 地球物理学报， 56(11): 3868--3876.

An Z H， Du X B， Tan D C， Fan Y Y， Liu J， Cui T F. 2013. Study on the geo-electric field variation of Sichuan LushanMS7.0 and WenchuanMS8.0 earthquake[J].ChineseJournalofGeophysics， 56(11): 3868--3876 (in Chinese).

汤吉， 詹艳， 王立凤， 徐建郎， 赵国泽， 陈小斌， 董泽义， 肖骑彬， 王继军， 蔡军涛， 徐光晶. 2008. 5月12日汶川8.0级地震强余震观测的电磁同震效应[J]. 地震地质， 30(3): 739--745.

Tang J， Zhan Y， Wang L F， Xu J L， Zhao G Z， Chen X B， Dong Z Y， Xiao Q B， Wang J J， Cai J T， Xu G J. 2008. Coseismic signal associated with aftershock of theMS8.0 Wenchuan earthquake[J].SeismologyandGeology， 30(3): 739--745 (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 WenchuanMS8.0 earthquake[J].ChineseJournalofGeophysics， 53(3): 526--534 (in Chinese).

杨皓， 詹艳， 赵凌强， 陈小斌， 姜峰. 2015. 2013年甘肃岷县漳县MS6.6地震震电磁信号和同震信号观测[J]. 震灾防御技术， 10(增刊): 785--793.

Yang H， Zhan Y， Zhao L Q， Chen X B， Jiang F. 2015. The seismo-electromagnetic signal and coseismic signal of 2013 Minxian-ZhangxianMS6.6 earthquake[J].TechnologyforEarthquakeDisasterPrevention， 10(Suppl): 785--793 (in Chinese).

张丹， 任恒鑫， 黄清华. 2013. 孔隙介质地震电磁信号的数值模拟研究[J]. 地球物理学报， 56(8): 2739--2747.

Zhang D， Ren H X， Huang Q H. 2013. Numerical simulation study of co-seismic electromagnetic signals in porous media[J].ChineseJournalofGeophysics， 56(8): 2739--2747 (in Chinese).

Chen X F. 1999. Seismogram synthesis in multi-layered half-space: PartⅠ.Theoretical formulations[J].EarthquakeResearchinChina， 13(2): 149--174.

Frenkel J. 2014. On the theory of seismic and seismoelectric phenomena in a moist soil[J].JEngMech， 131(9): 879--887.

Gao Y X， Hu H S. 2010. Seismoelectromagnetic waves radiated by a double couple source in a saturated porous medium[J].GeophysJInt， 181(2): 873--896.

Gao Y X， Chen X F， Hu H S， Zhang J. 2013. Early electromagnetic waves from earthquake rupturing:Ⅰ.Theoretical formulations[J].GeophysJInt， 192(3): 1288--1307.

Garambois S， Dietrich M. 2001. Seismoelectric wave conversions in porous media: Field measurements and transfer function analysis[J].Geophysics， 66(5): 1417--1430.

Haartsen M W， Pride R S. 1994. Modeling of coupled electroseismic wave propagation from point sources in layered media[C]∥SEGTechnicalProgramExpandedAbstracts. Los Angeles, California: SEG, 13(1): 1155--1158.

Haartsen M W， Pride S R. 1997. Electroseismic waves from point sources in layered media[J].JGeophysRes， 102(B11): 24745--24769.

Honkura Y， Matshshima M， Oshiman N， Tunçer M K， Bari， Ito A， Iio Y， Iikara A M. 2002. Small electric and magnetic signals observed before the arrival of seismic wave[J].EarthPlanetsSpace， 54(12): e9--e12.

Hu H S， Gao Y X. 2011. Electromagnetic field generated by a finite fault due to electrokinetic effect[J].JGeophysRes， 116(B8): B08302.

Kennett B L N. 1984. Seismic wave propagation in stratified media[J].GeophysJRastrSoc， 86(1): 219--220.

Matsushima M， Honkura Y， Oshiman N， Baris S， Tuncer M K， Tank S B， Celik C， Takahashi F， Nakanishi M， Yoshimura R， Petras R， Komut T， Tolak E， Ito A， Iio Y， Isikara A M. 2002. Seismo-electromagnetic effect associated with the Izmit earthquake and its aftershocks[J].BullSeismolSocAm， 92(1): 350--360.

Okubo K， Takeuchi N， Utsugi M， Yumoto K， Sasai Y. 2011. Direct magnetic signals from earthquake rupturing: Iwate-Miyagi earthquake ofM7.2， Japan[J].EarthPlanetSciLett， 305(1/2): 65--72.

Pride S R. 1994. Governing equations for the coupled electromagnetics and acoustics of porous media[J].PhysRevB， 50(21): 15678--15696.

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[J].EarthquakeScience， 23(2): 167--176.

Thompson A H， Gist G A. 2012. Geophysical applications of electrokinetic conversion[J].LeadEdge， 12(12): 1169--1173.

Zhu Z Y， Haartsen M W， Toksöz M N. 2000. Experimental studies of seismoelectric conversions in fluid-saturated porous media[J].JGeophysRes， 105(B12): 28055--28064.

Zyserman F I， Gauzellino P M， Santos J E. 2010. Finite element modeling of SHTE and PSVTM electroseismics[J].JApplGeophys， 72(2): 79--91.

On numerical simulation method of the seismic electromagnetic anomaly in porous stratified media

Yu Ziye*

(InstituteofGeophysics，ChinaEarthquakeAdministration，Beijing100081，China)

This paper focuses on numerical simulation of the seismic electro-magnetic anomaly generation process in isotropic layered medium by means of propagation-matrix method. For the layered model， some are fluid-saturated porous medium， and the others are linear elastic medium. During the simulation， the solid field is taken as a principal field and calculated individually， and its deformation is taken as the equivalent fluid source. The simulation results show that when the seismic wave goes through the boundary of the fluid-saturated porous medium， the electromagnetic (EM) wave is generated simultaneously， which can be used to illustrate the seismic EM anomaly accompanied by seismic source or the seismic wave. On the other hand， it can elucidate EM signal coupled with the foreshock or the tremor.

electrokinetic effect； electromagnetic anomaly； fluid-saturated porous medium； electromagnetic field simulation

国家重大科学基础设施项目东半球空间环境地基综合监测子午链“子午工程”资助.

2016-01-06收到初稿， 2016-04-13决定采用修改稿.

10.11939/jass.2016.06.006

P319.1+2

A

于子叶. 2016. 层状多孔介质的地震电磁异常数值模拟方法研究. 地震学报, 38(6): 869--877. doi:10.11939/jass.2016.06.006.

Yu Z Y. 2016. On numerical simulation method of the seismic electromagnetic anomaly in porous stratified media.ActaSeismologicaSinica, 38(6): 869--877. doi:10.11939/jass.2016.06.006.

*通讯作者 e-mail: cangye@hotmail.com