轴向各向异性地层瞬变电磁三分量响应特征

郭 建 磊

(中煤科工集团 西安研究院有限公司，陕西 西安 710077)

0 引言

在沉积岩地区由于层理发育导致地下介质电阻率随电流方向发生变化，地层表现出很强的各向异性导电特性。地层电导率各向异性的研究最早可以追溯到1920年，Schlumberger在进行电法勘探时就发现了沉积地层的电性各向异性现象[1]。在近几十年的历程中，较多的地球物理学者对地层各向异性进行了研究[2-3]，目前学者们对电导率各向异性与矢量电磁场的研究取得了初步进展[4]。

对于地层电导率而言，各向同性与各向异性的区别在于地层电导率与方向的关系。地层各向同性指地层电导率不随方向的改变而变化，地层各向异性指地层电导率随方向的改变而变化。地层多发育为连续沉积地层[5]，当沉积层具有明显层理时，导电大地不同方向的导电能力发生变化，地下岩层表现出宏观电导率各向异性；另外地层在形成过程中有大量的孔隙和裂隙，孔隙和裂隙会被其他固体物质或者流体物质充填，最终导致地层具有强烈的各向异性特性[6-8]；地层开采会产生一定的扰动，改变原岩的构造应力，引起上覆岩层的破坏和采空区裂隙分布的改变，产生不同形态、不同高度和密度的采空裂隙和裂缝[9]，经过一定时间后裂隙和裂缝会充水形成裂隙含水体，由于构造和含水充填的原因，裂隙含水体会产生电导率各向异性特性[10-11]。

瞬变电磁法目前是探测地层电性分布、地层富水性以及采空区积水的主要地球物理探测方法，在煤田、金属矿等领域应用广泛。随着计算机的发展以及瞬变电磁理论的进步，众多学者在瞬变电磁法各向异性三维正演、瞬变电磁三分量解释及应用等方面做了大量工作，比如：周建美等基于有限体积法实现双轴各向异性瞬变电磁三维正演，通过正演模拟发现水平方向电导率变化严重影响瞬变电磁磁场垂直分量响应，垂直方向电导率变化对磁场垂直分量响应几乎没有影响[12]；刘亚军等基于有限体积法实现任意各向异性瞬变电磁三维正演，发现主轴水平方向电导率各向异性严重影响瞬变电磁垂直分量的磁场响应，且倾斜各向异性比水平各向异性的影响更为严重，通过水平各向异性信号能较清晰地判断出各向异性主轴方向[13]；程久龙等基于时域有限差分法研究地层各向异性对矿井瞬变电磁的影响特征，发现相比HTI-X介质，HTI-Y介质和VTI介质各向异性系数的变化对瞬变电磁场的影响显著，并结合实例进行了说明[14]。

目前，瞬变电磁法无论从硬件到软件均已经实现了三分量数据采集与解释[15-16]，并且模型试验和实际资料处理结果表明三分量解释方法能够利用更多的瞬变电磁数据，获得更高分辨率的处理结果[17]。瞬变电磁法的解释方法主要有视电阻率解释法[18]、浮动薄板解释法[19]、烟圈理论解释法[20]、一维反演解释法[17]和拟地震成像法[21]，但上述解释方法均是基于将地层假设为各向同性介质为前提进行。前人研究结果表明介质各向异性特性会严重影响瞬变电磁响应，尤其是三维异常体的各向异性影响会更为突出，因此，基于各向同性模型的传统成像及反演成像技术不再适用于各向异性数据的处理解释，这导致传统的基于各向同性模型的成像方法和一维反演技术不再适用于各向异性数据的处理与反演解释，并且随着瞬变电磁技术的发展以及计算机计算性能的快速提升，三分量采集及三分量联合反演解释等成为目前实现瞬变电磁精细探测的前沿技术[15-17]。前人的研究工作更多的是探讨电导率各向异性对瞬变电磁磁场垂直分量响应的影响，并未涉及对磁场水平分量响应的影响。研究地层各向异性和含水体各向异性对回线源瞬变电磁三分量磁场响应的影响特征，有利于地质精细化勘探。该项研究最便捷的途径就是三维正演计算。

时域有限差分算法经过近些年的快速发展，其理论不断完善并趋于成熟。孙怀凤等[22]通过安培环路定理方程将矩形回线源电流密度加入麦克斯韦方程组，实现了回线源的加入，实现TEM全波场正演模拟；卢绪山[23]在孙怀凤的基础上利用Open-ACC Fortran语言进行了以异构并行算法开发，加快了计算速度；李展辉等[24]将复频率完全匹配层吸收边界应用到三维瞬变电磁时域有限差分正演中，给出了CFS-PML在TEM正演中参数设置准则；孙怀凤等[25]后续又在常规FDTD的迭代基础上开发了Crank-Nicolson FDTD正演算法，在保证求解精度的同时计算效率大幅提高。

本文基于三维时域有限差分算法，引入轴向电导率张量构建各向异性电磁场控制方程，进而推导无源和有源条件下轴向各向异性电场三分量迭代公式，其中磁场分量的迭代计算公式与各向同性磁场迭代计算公式一致，以电流密度的形式加入各向异性麦克斯韦方程组安培环路定理实现电流源的加载，实现各向异性三维瞬变电磁三分量正演。通过与各向同性及各向异性半空间模型解析结果对比验证算法的准确性，分别建立各向异性半空间模型、层状模型、含水体模型并进行正演计算，研究x、y、z轴各向异性分别对瞬变电磁∂Bx/∂t、∂By/∂t和∂Bz/∂t响应的影响。

1 轴向各向异性有限差分瞬变电磁正演算法

1.1 正演算法原理

各向异性、有耗、无源媒质中的Maxwell方程组表示为[14]：

式中：E为电场强度，B为磁感应强度，H为磁场强度，ε为虚拟介电常数，t为时间，σ为张量电导率。

地球物理勘探中一般忽略位移电的作用。为了构成显式的时间迭代格式，满足三维时域有限差分计算要求，加入人工虚拟介电常数项，方程(1b)变为：

(2)

采用经典Yee晶胞格式(图1)对计算区域进行网格剖分，电场分量位于六面体的棱边中心点，磁场分量位于六面体的面中心点，每一个电场(磁场)分量均被相应的4个磁场(电场)分量所包围，电场分量和磁场分量在时间轴上交替采样(图2)。

图1 Yee晶胞网格剖分示意Fig.1 Schematic diagram of Yee unit cell meshing

图2 电磁场时间采样分布示意Fig.2 Schematic diagram of time sampling distribution of electromagnetic field

空间上采用后向差分，时间上采用中心差分，将式(2)在tn+1/2时刻进行空间和时间离散，整理后得到电场强度随时间的推进公式：

(3)

为分析轴向各向异性瞬变电磁响应特征，将任意各向异性简化为轴向各向异性，σ为各个节点上等效轴向各向异性电导率张量，通常表示为：

(4)

将σ分别代入矩阵P、Q，求得:

(5)

以差分代替微分，采用非均匀网格形式[22]进行模型剖分，将式(3)在各分量所在空间节点处进行空间和时间离散，得到电导率各向异性下电场强度迭代表达式：

式中：i、j、k分别表示x、y、z方向网格剖分的节点号，Δx、Δy、Δz分别表示剖分网格的长度，Δt表示时间间隔。

采用回线源中加入阶跃电流的方式进行电流源加载，采用梯形波发射波形，考虑发射电流的上升沿、持续时间和下降沿，有源区域式(2)可修改为：

(9)

式中Js为源电流密度；因电流在水平方向，因此方程中仅存在Jsx和Jsy。在直角坐标系中加入电流密度展开后，将发射回线所在网格按照差分格式正常离散,得到有源区域网格电场迭代形式[22]：

磁场迭代公式不包含电导率项，各向异性磁场迭代公式与各向同性迭代公式相同，具体推导过程可以参看参考文献[22,23]。

瞬变电磁法是宽频带勘探方法，加载梯形波作为激发源。选用第一类边界条件，即在6个边界面上将电场和磁场分量全部强制赋零，计算过程中就需要满足时间域和空间域的稳定性条件，即Courant稳定性条件:

(12)

式中v表示电磁波的相速度，保证电磁场传播的因果关系。

1.2 正演算法验证

采用各向同性和各向异性半空间模型，通过与一维数值解的比较验证本文算法的精度。采用笛卡尔正交坐标系，模型中发射回线采用190 m×190 m的方形回线，回线源放置于xoy平面且回线源中心位于坐标系中心，发射电流为1 A，电流方向呈顺时针，上升沿和下降沿均为1 μs，脉冲持续时间10 ms，二次场采样时间为10 ms，模型剖分均匀网格尺寸为10 m，模型剖分网格数为221×221×200。

首先，建立各向同性半空间模型。模型电导率为0.01 S/m，考虑到水平分量磁场响应在回线源中心点位置结果为0，故将接收点选择在(-60 m,-60 m,0 m) 处，分别接收∂Bx/∂t、∂By/∂t、∂Bz/∂t响应，与一维解析结果[16]进行对比。对比结果(图3)显示：∂Bx/∂t、∂By/∂t、∂Bz/∂t响应的解析解与数值解的误差均在10%以下，在0.02～10 ms时间范围内误差均在5%以下，∂Bz/∂t的响应误差稍大于∂Bx/∂t和∂By/∂t的响应误差，接收点位于发射回线的对角线上，∂Bx/∂t与∂By/∂t响应曲线重合且响应误差一致。

图3 各向同性半空间模型垂直分量磁场响应(a)与相对误差(b)Fig.3 Isotropic half-space model vertical component magnetic field response (a) and relative error (b)

然后，建立如图4所示的各向异性半空间模型，模型轴向电导率参数如表1所示。在回线源中心点处接收∂Bz/∂t响应，图5为各向异性半空间模型有限差分数值解与一维解析解[30]的对比结果。通过图5可以发现：3个模型的∂Bz/∂t响应一维解析结果与数值解的误差均在5%以下，满足计算精度要求。

图4 各向异性半空间模型Fig.4 Anisotropic half-space models

图5 各向异性半空间模型垂直分量磁场响应(a)与相对误差(b)Fig.5 Anisotropic half-space model vertical component magnetic field response (a) and relative error (b)

表1 图4模型的轴向电导率Table 1 Axial conductivity of the Fig.4 model

通过各向同性和各向异性半空间模型对比，验证了本文算法的准确性与精度。本文主要分析地层的轴向各向异性对磁场三分量的影响特征，因此为考虑成图方便，下文中将接收点坐标点选择在(-60 m,-60 m,0 m) 处。

2 各向异性瞬变电磁响应特征

2.1 各向异性半空间模型

基于上述各向异性半空间模型进行正演，研究回线源瞬变电磁三分量响应(图6)。图6中，由于在双对数坐标系中无法表示负值，因此将感应电动势为负值的部分取绝对值后绘制于双对数坐标系中，后续绘图采用如上约定。通过∂Bx/∂t响应曲线

图6 各向异性半空间模型瞬变电磁磁场三分量响应Fig.6 Three-component TEM response of anisotropic half-space models

可以发现：增大x轴电导率后响应幅值整体增大；增大y轴电导率后曲线出现“变号”现象，“变号”时间约为0.1 ms，响应幅值增大且增大幅度大于增大x轴电导率；增大z轴电导率对响应几乎没有影响，与背景响应重合。通过∂By/∂t响应曲线可以发现：增大x轴电导率后曲线出现“变号”现象，“变号”时间为0.1 ms左右，响应幅值增大且增大幅度大于增大y轴电导率；增大y轴电导率后响应幅值整体增大；增大z轴电导率对响应几乎没有影响，与背景响应重合。通过∂Bz/∂t响应曲线可以发现：增大x、y轴电导率后响应幅值整体增大且响应曲线重合，此时接收点位于对角线上，x、y轴电导率各向异性对∂Bz/∂t响应的影响程度一致；增大z轴电导率对∂Bz/∂t响应几乎没有影响，与背景响应重合。

为研究轴向电导率各向异性对三分量磁场等值线整体分布的影响，基于上述模型，选取关断后0.1 ms时刻的三分量磁场值绘制地面三分量磁场强度等值线(图7)，并以各向同性半空间模型三分量磁场响应分布图为参考背景。通过∂Bx/∂t等值线分布图可知，增大x轴电导率后等值线图左、右两侧幅值增强，等值线渐变为近似椭圆形；增大y轴电导率后等值线图左、右两侧幅值增强且增大幅度大于增大x轴电导率，等值线图左、右两侧出现“变号”现象；增大z轴电导率后等值线图几乎没有变化。通过∂By/∂t等值线图可知，增大x轴电导率后等值线图左、右两侧幅值增强且增大幅度大于增大y轴电导率，等值线图左、右两侧出现“变号”现象；增大y轴电导率后等值线图左、右两侧幅值增强，等值线渐变为椭圆形；z轴电导率增大后等值线图几乎没有变化。通过∂Bz/∂t等值线图可知，增大x、y轴电导率后等值线图幅值增强，等值线由同心圆向外扩散模式变为2个极值椭圆中心向外扩散模式，增大z轴电导率后等值线图几乎没有变化。

图7 地面磁场三分量0.1 ms时刻分布情况Fig.7 Ground three-component magnetic field distribution at 0.1 ms

图6和图7均证明水平方向电导率对磁场三分量的影响较大，垂直方向电导率对磁场三分量基本没有影响，这主要是因为回线源产生的电流主要是水平方向的，因此三分量响应主要受水平方向电导率的影响，垂直方向电导率的影响很小；y方向的电导率对∂Bx/∂t响应的影响大于x方向的，x方向的电导率对∂By/∂t响应的影响大于y方向的。

2.2 各向异性层状地层

沉积岩具有层理结构。当沉积层具有明显层理时，导电大地不同方向的导电能力发生变化，地下岩层表现出宏观电导率各向异性；另外地层在形成的过程中有大量的孔隙和裂隙，孔隙和裂隙会被其他固体物质或者流体物质充填，最终导致地层具有强烈的各向异性特性。因此，研究层状地层各向异性对回线源瞬变电磁三分量磁场响应的影响特征具有现实意义。

以各向同性半空间模型为背景模型，设地层电导率为0.01 S/m，基于背景模型设置不同深度的地层为各向异性，分别改变上部地层(深度在0～100 m)、中间地层(深度在100～200 m)、下部地层(深度在100 m以下的地层)以及全部地层(地表以下地层)为各向异性地层，改变各向异性地层的3个方向电导率，电导率张量参数如表2所示。

表2 层状模型的轴向电导率Table 2 Anisotropy parameters of axial conductivity

采用与前文1.2章节中正演算法相同的发射源和发射电流参数，经过三维计算得到层状各向异性地层瞬变电磁三分量响应结果(图8)。图8表明：增大x轴电导率，改变上部地层与改变全部地层的∂By/∂t响应出现“变号”现象；增大y轴电导率，改变上部地层与改变全部地层的∂Bx/∂t响应出现“变号”现象；增大水平方向电导率，改变上部地层与改变全部地层的三分量响应曲线重合，改变中间地层和改变下部地层的三分量响应曲线重合，改变上部地层与改变全部地层的三分量响应幅值大于改变中间地层和改变下部地层的三分量响应幅值；增大z轴电导率，所有模型的瞬变电磁三分量响应曲线全部重合。

通过图8再次验证了:水平方向电导率对磁场三分量的影响较大，垂直方向电导率对磁场三分量基本没有影响;y方向电导率对∂Bx/∂t响应的影响大于x方向电导率，x方向电导率对∂By/∂t响应的影响大于y方向电导率;浅部地层的电导率各向异性对三分量磁场响应的影响占主导成分。

图8 各向异性层状模型瞬变电磁三分量响应Fig.8 Three-component TEM response of anisotropic layered models

2.3 各向异性含水体模型

设置各向异性含水体模型，研究含水体各向异性对瞬变电磁磁场三分量响应的影响。各向异性含水体模型如下：在电导率为0.01 S/m的各向同性半空间中放置一个顶部埋深50 m、规模为110 m×110 m×50 m的含水体。当含水体为各向同性时，轴向电导率均为1 S/m；分别增大含水体的x、y、z轴电导率,设置模型为各向异性，电导率张量如表3所示。

表3 含水体模型的轴向电导率Table 3 Anisotropy parameters of axial conductivity

采用与前文相同的发射源和发射电流参数，将含水体中心分别放于发射源正下方、x轴正向55 m、y轴正向55 m处(图9)。图10为含水体模型瞬变电磁磁场三分量响应结果，图11为其x、y、z轴各向异性瞬变电磁水平分量响应对比。

图9 含水体相对于回线源位置俯视图Fig.9 Top view of the position of the hydrate relative to the loop source

由图10a和图11可知，当含水体位于源中心正下方时，增大x、y轴电导率后，∂Bx/∂t、∂By/∂t、∂Bz/∂t响应幅值大于增大z轴电导率和各向同性模型，且x轴各向异性模型与y轴各向异性模型的∂Bx/∂t、∂By/∂t、∂Bz/∂t响应曲线基本重合，z轴各向异性模型与各向同性模型的∂Bx/∂t、∂By/∂t、∂Bz/∂t响应曲线基本重合。此时，接收点x方向距含水体中心点距离与y方向相等，x轴各向异性模型的∂By/∂t响应幅值大于∂Bx/∂t响应幅值，y轴各向异性模型的∂Bx/∂t响应幅值大于∂By/∂t响应幅值，z轴各向异性模型的∂Bx/∂t响应曲线与∂By/∂t响应曲线重合;这再次说明了x轴各向异性对∂By/∂t分量响应的影响大于∂Bx/∂t，y轴各向异性对∂Bx/∂t分量响应的影响大于∂By/∂t。

将含水体沿x轴正向移动55 m并保持深度不变，此时接收点x方向距含水体中心点距离大于y方向。通过图10b可知，增大y轴电导率后∂Bx/∂t、∂By/∂t、∂Bz/∂t响应幅值增大且幅度大于增大x轴电导率。这主要是因为异常体位于x轴，根据“烟圈理论”，通过异常体位置的磁力线主要与y轴电导率方向平行，受其影响较大。将含水体沿y轴正向移动55 m并保持深度不变，此时接收点x方向距含水体中心点距离小于y方向。通过图10c可知，增大x轴电导率后∂Bx/∂t、∂By/∂t、∂Bz/∂t响应幅值增大且幅度大于增大y轴电导率，这主要是因为异常体位于y轴，根据“烟圈理论”，通过异常体位置的磁力线主要与x轴电导率方向平行，受其影响较大。通过图11可知，模型位于回线源正下方时响应反映均大于模型中心位于偏x轴正向和偏y轴正向，根据“烟圈理论”，位于回线源正下方的异常体受回线源所激发的二次场能量最强。

图10 各向异性含水体瞬变电磁三分量响应Fig.10 Transient electromagnetic three-component response of anisotropic aquifer

图11 x、y、z轴各向异性瞬变电磁水平分量响应对比Fig.11 x,y,z axis anisotropic transient electromagnetic horizontal component response comparison diagram

3 结论

本文基于三维时域有限差分算法，引入轴向电导率张量构建各向异性的电磁场控制方程，实现电导率轴向各向异性瞬变电磁三维正演。通过与各向同性及各向异性半空间模型一维解析结果对比,验证算法的准确性。

建立各向异性半空间模型、层状模型、含水体模型并进行正演计算，研究x、y、z轴各向异性分别对瞬变电磁∂Bx/∂t、∂By/∂t和∂Bz/∂t响应的影响,得到如下结论：

1)通过各向异性半空间模型发现：水平方向电导率对磁场三分量的影响较大，垂直方向电导率对磁场三分量基本没有影响，三分量响应主要受水平方向电导率的影响，垂直方向电导率的影响很小；y方向电导率的∂Bx/∂t响应的影响大于x方向电导率，同理,x方向电导率的∂By/∂t响应的影响大于y方向电导率；

2)通过层状各向异性模型发现：各向异性地层所处深度不同对瞬变电磁三分量磁场响应影响特征不同，认为浅部地层的电导率各向异性对三分量磁场响应的影响占主导成分；

3)通过各向异性含水体模型发现：接收点与各向异性含水体的相对位置不同会影响瞬变电磁三分量磁场响应，当接收点x方向距离含水体中心小于y方向时，x轴各向异性对三分量响应影响程度大于y轴各向异性，反过来，当接收点x方向距离三维地质体中心距离大于y方向时，x轴各向异性对三分量响应影响程度小于y轴各向异性，并且模型位于回线源正下方时响应反映均大于模型中心位于偏x轴正向和偏y轴正向。

研究发现地层电导率轴向各向异性对瞬变电磁磁场三分量响应具有较大影响，为提高瞬变电磁精细勘探，需要进一步研究不同地电类型、不同探测方式的轴向各向异性瞬变电磁三分量磁场响应特征，并进一步研究轴向电导率对三分量磁场的物理影响机理，为三分量联合探测及精细反演提供指导依据。

致谢：本研究在孙怀凤等开发的TEM3DFDTD程序基础上完成，在此表示感谢！