1 引言

地震波在地下含流体孔隙介质中传播时其流相与固相之间会发生相对运动，不仅导致地震衰减和速度频散现象，还会诱发以双电层动电现象为基础的震电效应(Pride，1994；刘洪和李幼铭，1994).震电效应的本质是物质的机械能与电磁能之间的相互耦合作用机制，最初是由Thompson (1936)发现，Pride (1994)建立了针对双相介质电磁场和弹性波场耦合的控制方程.波达波夫(1996)详细论述了震电效应的物理产生机制，胡恒山等(2003)对比研究了波达波夫方程组与Pride震电耦合方程组的传播机理，随后人们针对震电效应的数值模拟方法、近场球面波特征、远场平面波特征以及储层物性影响关系等方面进行了深入研究(Haines和Pride，2006；高永新和胡恒山，2009；崔志文等，2011；丁浩然等，2016；郑海霞等，2016).通过前人大量的实验室测量研究，人们逐渐明确了在不同物理岩样条件、不同频带下震电效应的产生条件和传播规律，并利用其逐步展开实际应用研究(Jouniaux and Pozzi, 1995；石昆法，2001；Zhu and Toksöz，2005；王军等，2010；Li et al., 2014；王军等, 2015, 2016).Schakel等(2012)对震电效应的理论预测和实验测量进行了对比分析研究.由于地下油气储层物性参数(流体类型，孔隙度及渗透率等)与震电效应有着十分密切的联系，一直是人们持续研究的热点课题.另外，地震波在地下含流体孔隙介质中传播时诱发的地震波衰减和频散机理也是人们研究的重要内容.从20世纪50年代以来，许多学者对地震波在双相介质中的衰减机理做了大量研究工作，从研究地震波衰减因素尺度的角度出发，可将地震波衰减机制分为三类，宏观尺度、微观尺度和中观尺度(Müller et al., 2010).由宏观尺度“Biot流”无法解释地震频带(1~1000 Hz)的衰减和频散(Biot, 1956a, b，1962；Dvorkin and Nur, 1993).White等(1975)在考虑中观尺度流体斑块饱和情况时，理论模拟了地震波在地震频带的衰减和频散(White，1975；White et al., 1975).随后人们针对中观尺度的衰减和频散机理研究做了大量的理论研究工作(Vogelaar and Smeulders, 2007；刘炯等，2009；邓继新等，2012；吴建鲁等，2015).值得一提的是，微观尺度喷射流也被认为是地震波频散和衰减的重要原因之一.Song等(2016)基于其建立的喷射流模型，采用实验反演得到的纵横比分布计算了纵波的频散和衰减曲线，结果表明：喷射流机制也能够预测可观的地震频段频散.本文只考虑中观尺度非均匀性导致的流动问题.

在构建特定理想化岩石物理模型的基础上发展特定的地震波衰减理论是研究地震波衰减和频散特征的重要手段和方法，由于实际地下含流体饱和介质的岩性、流体斑块饱和、岩石孔渗性差异等因素的复杂变化，使得刻画的理论模型往往无法与实验测量结果精确匹配，从而验证理论的精确性和适用条件等具有局限性(Batzle et al., 2006；Madonna and Tisato, 2013; Pimienta et al., 2015a, b).Masson等(2006)将低频形式的动态渗透率关系式转换至时间域后，采用准静态时间域交错网格有限差分方法，利用Biot双相介质波动方程首次针对中观尺度非均质模型进行了虚岩石物理理论分析研究，随后将该方法推广至全频带形式(Masson et al., 2006；Masson and Pride, 2007).Quintal等(2011)和Rubino和Holliger(2012)针对中观尺度非均质采用有限元的剖分思想开展了地震波衰减和速度频散的研究工作.Jougnot等(2013)采用准静态有限元方法研究了含流体孔隙介质中充填中观尺度裂隙情况下震电效应.随后，Monachesi等(2015)推导了一维White模型电势与电荷密度的解析表达式，并与数值模拟结果做了对比分析.

本文基于Biot(1941)固结理论的准静态方程，在频率域中采用空间有限差分方法，正演模拟虚岩石物理岩样的地震波衰减和震电效应.与时间域虚岩石物理方法相比，该方法既可以直接求取任一频率下地震波衰减和频散速度，便于预测分析实际岩样测量结果，又消除了针对频率域的动态渗透率公式进行近似化简所带来的误差，也避免了讨论岩样外表面施加的力源函数表达式及时间剖分的稳定性条件等问题.通过数值计算周期性层状介质模型，验证了本文所描述方法的有效性，并进一步求取分析了周期性层状介质两种不同特征单元的渗流电流密度及电势，数值模拟结果表明由中观尺度“局域流”引起的震电效应的电势振幅数值完全在实验室测量范围之内.随后，分析研究了四种不同高渗介质占比值模型的地震衰减、速度频散及震电效应特征，最后，将本文提出的震电效应数值计算方法推广至二维，并求取了二维斑块饱和模型的地震波衰减、速度频散、电势振幅和相位数值结果.

2 双相介质传播方程

由于本文旨在研究中观尺度非均质引起地震频带震电效应，故采用Biot(1941)所推导的低频条件下的线性固结耦合的准静态方程，当所研究的频带范围小于Biot临界频率ω_B时，孔隙流体粘滞边界层效应可忽略不计，体力散度为零，利用该方程可以精确描述地震波在低频条件下的传播特征，线性固结耦合的准静态Biot方程在频率—空间域的一维控制方程为

(1)

其二维控制方程形式为

(2)

其中，τ_xx、τ_zz分别为X、Z方向的正应力，P为流体压力，U_s、U_f为固体位移、渗流位移；λ、μ分别为骨架的拉梅系数和剪切模量，为Biot固结系数，K_d、K_s为骨架和基岩的体积模量，η为流体黏滞系数，k₀为零频极限条件下介质渗透率，ω为圆频率，

(3)

其中，φ为岩石孔隙度，K_f为孔隙流体体积模量.

3 双相介质震电耦合方程

在含流体孔隙介质中，固体表面通常吸附负离子，流体中携带较多的阳离子，流体相可分为吸附层和扩散层，吸附层与流固界面相邻，厚度仅为纳米级别，其外为扩散层.由于扩散层中阳离子受固体表面处阴离子束缚减小，当地震波在双相介质中传播时，产生的流体压力差引发流-固相对运动，且孔隙中流体携带扩散层中的阳离子一起运动，进而使得流体压强相对较低的区域中堆积较多阳离子，流体压强较大区域的阴离子聚集较少，进而产生电势差，即流动电势，形成电流，即为震电效应(石昆法，2001).孔隙流体中渗流位移与电荷密度及电压的定量表达式为(Sill，1983; Jardani et al., 2010; Jougnot et al., 2013)

(4)

(5)

其中，为渗流电流密度，有效电荷密度，φ为电势，σ为电导率.针对砂岩岩样电导率的经验关系表达式为(Archie，1942)

(6)

其中，m=2为胶结指数，σ_w为水的电导率.

Jardani等(2007)给出了针对含水砂岩储层介质有效电荷密度与渗透率之间的经验关系式：

(7)

为分析不同孔隙度条件的震电效应，采用Kozeny-Carman关系式定量计算砂岩岩样孔隙度与渗透率之间关系(Mavko et al., 2009)：

(8)

其中，b为孔隙结构因子，与孔隙介质的挠曲度相关，d为平均岩石颗粒尺寸，本文给定b=0.003、d=8×10^－5m来研究砂岩岩样.

另外，本文利用Krief等(1990)提出的经验关系式来计算不同孔隙度条件下岩石骨架的体积模量和剪切模量：

(9)

(10)

其中，μ_s为基岩的剪切模量.

4 频率域准静态虚岩石物理方法

当外界施加一定的压力P₀时，岩石物理模型受力产生形变，通过对上述两种孔隙介质地震波控制方程微分项进行常规空间离散，即中心空间离散插值格式，如公式(11)，再引入边界条件，即可组装成A·X=b的线性方程组，其中A为包含岩石物理参数信息的刚度矩阵，X=(τ_xx, τ_zz, τ_xz, P, U_s^x, U_s^z, U_f^x, U_f^z)为待求的物理量，b为边界条件.通过求解该方程组即可求取岩石物理模型各处的固体位移值(U_s^x, U_s^z)及正应力(τ_xx, τ_zz)，通过方程(13) 可求取模型的复变模量E(ω)，进而通过方程(12) 获取模型地震波衰减Q^－1(ω)和频散速度V_p(ω).

在计算获取模型渗流位移(U_x^f, U_f^z)值后，进一步通过其与渗流电流和电势之间的定量关系式(4、5) 获取模型的渗流电流和电势分布.

(11)

其中，dz为岩石物理模型的离散空间步长，N为岩石物理模型的总离散点数.

(12)

其中复变模量E(ω)可通过数值模拟获取的应力(τ_xx, τ_zz)和应变(ε_xx, ε_zz)求得(Masson and Pride, 2007)：

(13)

其中，应变通过数值模拟获取的固体位移求得，ε_xx=dU_s^x/dx，ε_zz=dU_s^z/dz.

在上下施加一定的垂直方向外力P₀时，为保证岩石物理模型的封闭条件，使得模型边界的渗流位移U_f为零，则模型的上下边界条件为

(14)

若模拟半空间孔隙介质时，为了保证边界处的流体压力连续，则用边界位置处的流体压力为零来代替封闭条件下的渗流位移为零(Kudarova et al., 2013).

在求解计算所建模型电势时，需满足模型边界的绝缘条件(第二类边界条件)及在模型的一端给定模型的参考电势值(第一类边界条件)：

(15)

5 数值计算结果及分析

为验证频率域一维虚岩石物理方法的有效性，利用周期性层状斑块饱和模型进行数值模拟计算，模型如图 2a所示，饱含水的模型介质中两种不同渗透特征水平地层周期性交互排列，横向无限延伸，仅考虑垂直方向一维应力-应变响应.高、低渗储层的层厚各为20 cm，空间网格步长为0.001 m，施加的外界压力P₀=10³ Pa，采用该模型两种不同特征单元模型(Model 01、Model 02) 进行虚岩石物理实验，Model 01为低渗储层中间嵌入高渗介质情况，Model 02为高渗储层下间嵌入低渗介质情况，如图 3所示，求取了这两种模型的地震波衰减和频散数值结果，并与解析解结果进行了对比(模型解析解公式可参见参考文献(White et al., 1975；吴建鲁等，2015)).与解析解相比，利用本文所描述的准静态频率域有限差分方法求解的两种模型各个频率地震波衰减和频散结果十分吻合，验证了该方法的有效性，也进一步说明了两种模型都可看作表征周期性层状介质模型衰减信息的特征单元.

图 4为针对周期性层状介质模型，正演模拟计算求取的两种特征单元不同频率应力、流体压力、固体位移及渗流位移，由图可见，模型的应力在不同空间位置、不同频率下为一常数，大小与外界施加的压力相等，这与所采用的线性准静态双相介质方程有关，如公式(1) 所示，其物理意义为应力的散度为零，即模型的应力处处相等，并且两种模型的应力分布一致.由图(4c—4h)可知，在不同空间位置和不同频率下流体压力、固体位移及渗流位移数值都是变化的，在高、低渗储层介质分界面处的渗流速度值最大，且两种模型不同频率条件下三者参数数值模拟结果不同.当频率趋于高频极限或低频极限时，在不同空间位置上述三者参数数值随频率的变化逐渐减小，存在高、低频极限特征.

图 5为周期性层状介质模型两种特征单元不同频率渗流电流密度、电势幅值及电势相位角的数值模拟结果.由于模型受压在不同渗透率介质交界处附近会发生流体相对流动(图 4g、4f)，从而形成渗流电流密度，产生电势.由于低渗介质的有效电荷密度相对较高，在不同渗透率储层交界处低渗储层介质一侧产生较强的渗流电流密度(图 5a、5b)，进而求取电势，电势的振幅数量级在10^－3V以上，完全在实验室测量范围之内，且两种模型不同频率和不同空间位置处的电势振幅及相位结果具有明显的差别，说明对于双相介质中渗透率空间变化引起的震电效应差异相对于地震波衰减信息更加敏感，这也为利用双相介质震电效应进行地下储层和岩石岩样的空间渗透率分布特征研究提供了一种重要手段.另外，在高、低频极限处渗流电流密度及电势与频率无关，同样存在频率高、低频极限特征.

由于地下双相介质渗透率复杂多变，含流体介质储层介质的渗透性是评价油气勘探开发的重要指数.为分析研究不同渗透率介质混合情况下的震电效应特征，设计如图 6所示四种不同渗透率介质占比情况一维简单模型，高渗介质占比依次为0.2、0.4、0.6和0.8.所用参数如表 1所示.图 7为利用频率域有限差分准静态虚岩石物理方法求取的四种模型地震波衰减和频散，随着高渗介质的占比值逐渐增加，模型变得更加疏松，地震波衰减逐渐降低且衰减峰值频率向高频移动，地震波频散速度逐渐降低.

图 8为四种模型不同频率及空间位置电势幅值及相位的数值模拟结果.由图可见，随着高渗介质占比值越来越大时，高渗介质中电势振幅极大值向高频移动，且在其介质中不同空间位置的电势(振幅及相位)不变.在低渗介质中电势(振幅及相位)分布特征相对复杂.图 9为四种模型顶部位置不同频率电势振幅及相位角数值计算结果.随着高渗介质占比值逐渐增加，电势振幅及相位角向高频移动，电势振幅幅值小幅增加.

从四种模型模拟获取的地震波衰减信息及震电效应模拟结果可知，地下双相介质储层从低渗储层到高渗储层变化时，相应的地震波衰减、频散、电势振幅及相位变化特征都不相同，当我们从实际地下储层介质或岩石物理实验中获取其更丰富的地震及电磁信息时，对于地下介质的研究也就越深入，预测计算的储层渗透率特征越准确.

进一步将本文提出的震电效应数值模拟方法进行了二维测试.首先给了在低孔、低渗背景中内嵌高孔、高渗孔隙介质的斑块饱和模型(即双重孔隙介质)，如图 10所示.模型在x、z方向的长度为50 cm，空间步长dx、dz为0.005 m，施加的外界压力P₀=10³ Pa，介质参数如表 1所示.通过方程(12)、(13) 计算求得的模型地震波衰减和频散速度如图 11所示.结果显示：模型的地震波衰减表现出“双峰”特征，衰减值较大的峰值频带在地震频带到测井频带，衰减值较小的峰值频带从声波频带到超声波频带，在相应频带内，速度也表现出较强的频散，与Pride等(1994)、巴晶(2010)和Ba等(2011)基于双重孔隙波动理论的解析结果相似.

图 12为50 Hz条件下x、z方向正应力、剪切应力、流体压力及x、z方向固体位移及渗流位移数值模拟结果.结合方程(3)、(4)，可求得二维斑块饱和模型的电势数值模拟结果，图 13为50 Hz条件下计算求得的电势振幅和相位角空间分布结果，同样可以看出，电势的振幅数量级在10^－3V以上，进一步说明了基于中观尺度“局域流”引起的震电效应作用较强，为实验室测量和实际应用提供了一定的借鉴和指导作用.

5 认识与讨论

地震波在地下含流体孔隙介质中传播时，会引起中观尺度的“局域流”，进而产生地震波震电效应.在前人针对中观尺度“局域流”研究的基础上，利用Biot(1941)线性固结双相介质准静态方程，采用频率域有限差分正演模拟方法，利用周期性层状介质模型验证了本文所提出方法的有效性，数值计算了周期性层状介质两种不同特征单元的渗流电流密度及电势，其电势振幅数量级在实验室测量范围之内.进而分析研究了四种不同高渗介质占比值的地震衰减及震电效应特征，随着储层介质渗透性逐渐增加，地震波衰减及震电效应特征变化关系也不相同，这也使得我们可以通过从不同侧面分析研究地下储层介质物性参数与岩石岩样的各种物理参数(地震波衰减、频散、渗流电流密度及电势等)的定量关系，更加准确的获取地下储层介质物性信息.随后，我们又将本文提出的震电效应数值计算方法推广至二维，并求取了二维斑块饱和模型的地震波衰减、速度频散、电势振幅和相位角数值结果.结果显示，求取的地震波衰减曲线呈现“双峰”特征，与前人基于双重孔隙波动理论的解析结果类似.在地震频带，中观尺度“局域流”衰减机制引起的震电效应数值结果较强，这也为我们从实验室测量和实际应用等方面分析含流体双相介质震电效应提供了借鉴和指导.

致谢

十分感谢评审专家对稿件的审阅以及提出了宝贵的修改意见.