CN110119586A - 轴向电导率各向异性瞬变电磁三分量三维fdtd正演方法 - Google Patents

Abstract

本发明涉及一种轴向电导率各向异性瞬变电磁三分量三维FDTD正演方法。本发明引入轴向电导率张量，构建基于轴向电导率各项异性的电磁场控制方程，将电导率张量代入控制方程推导无源和有源轴向电导率各向异性电场三分量迭代公式，可以实现轴向各向异性三维模型三分量计算。

Description

轴向电导率各向异性瞬变电磁三分量三维FDTD正演方法

技术领域

本发明涉及一种瞬变电磁三分量三维FDTD正演方法，属于地球物理探测技术领域，具体是涉及一种轴向电导率各向异性瞬变电磁三分量三维FDTD正演方法。

背景技术

瞬变电磁法(Transient Electromagnetic Method,TEM)又称时间域电磁法(TimeDomain Electromagnetic Method,TDEM)是一种建立在电磁感应原理基础上的时间域人工源电磁探测方法。利用不接地回线或接地线源向地下发送一次脉冲磁场，在此激发下，地下地质体中激励起的感应涡流将产生随时间变化的感应电磁场，利用线圈或接地电极观测二次场从而达到探测地下地质体的目的。瞬变电磁法以其低成本、高效率、高精度、避免高阻层屏蔽等优点已经被广泛的应用于金属矿资源勘探、煤田采空区探测、环境与水文地质调查等领域。

目前瞬变电磁法主要利用单个(垂直)分量瞬变响应数据，随着技术的发展，瞬变电磁三分量采集与解释已经逐渐被应用到实践中，三分量瞬变电磁响应数据包含更多的地下地质体信息，利用垂直分量数据确定异常体的深度，利用水平分量数据确定异常体的空间位置。目前越来越多的学者认识到电导率各向异性在地下构造中分布的普遍性以及对电磁响应数据解释的影响，尤其是对水平分量响应的影响。

复杂模型的瞬变电磁三维三分量正演主要靠数值方法完成，常用的有有限体积法(Finite Volume Method,FVM)、有限体积分方程法(Volume Integral Equation Method,VIEM)、有限单元法(Finite Element Method,FEM)、时域有限差分法(Finite DifferenceTime Domain,FDTD)。FDTD方法采用Yee结构化网格结构，直接在时间域求解Maxwell方程组获得瞬变电磁响应，能够精确的反映三维地质结构在瞬变电磁场激励下的响应规律。

期刊《石油地球物理学报》刊载文章《浅水域各向异性地层中的瞬变电磁响应分析》提出各项异性地层一维瞬变电磁正演方法，先求解浅水域各向异性地层中频率域响应，然后利用电场的时间域响应和频率域响应的关系求解各向异性地层中的瞬变电磁响应。该方法只能针对一维模型进行正演，计算结果仅为脉冲电场(V.Am-2)。

期刊《地球物理学报》刊载文章《双轴各向异性介质中回线源瞬变电磁三维拟态有限体积正演算法》提出采用模拟离散的有限体积法实现了双轴各向异性地层回线源瞬变电磁三维正演，将计算区域划分为一系列的控制体积单元，采用交错网格对控制方程进行模拟有限体积空间离散，包括旋度算子离散和空间内积离散。基于斯托克斯定理的旋度积分定义公式实现旋度算子离散，中点平均实现电导率双轴各向异性的空间内积离散，从而得到离散化的控制方程，通过均匀全空间中稳定电流回线源的磁场解析表达式得到回线源初始时刻的电磁场分布。该方法求解结果是Bz或者dBz/dt，无法实现轴向电导率各向异性模型水平分量瞬变响应的计算，并且无法计算隧道或者矿井轴向电导率各向异性模型，具有一定的局限性。

发明内容

本发明主要是解决现有技术所存在的不能计算起伏地形、不能适用于隧道和矿井模型、不能计算各向异性模型三分量场值等的技术问题，提供了一种轴向电导率各向异性瞬变电磁三分量三维FDTD正演方法。该方法构建基于轴向电导率各项异性的电磁场控制方程，将电导率张量代入控制方程，采用Yee晶胞格式和直角坐标系进行网格离散，从而实现轴向电导率各项异性三维模型三分量瞬变电磁响应特征的计算。

本发明的上述技术问题主要是通过下述技术方案得以解决的：

1.使用三维时域有限差分方法进行三维计算；

2.引入轴向电导率张量，构建基于轴向电导率各项异性的电磁场控制方程；

3.将电导率张量代入控制方程，推导无源轴向电导率各向异性电场三分量迭代公式；

4.以电流密度的形式加入电导率各向异性麦克斯韦方程组的安培环路定理中实现电流源的加载；

5.对电流源加载过程中考虑关断时间的影响。

因此，本发明具有如下优点：

1.引入轴向电导率张量，构建基于轴向电导率各项异性的电磁场控制方程，将电导率张量代入控制方程推导无源和有源轴向电导率各向异性电场三分量迭代公式，可以实现轴向各向异性三维模型三分量计算；

2.以电流密度的形式加入电导率各向异性麦克斯韦方程组的安培环路定理中实现电流源的加载，可以不受地表起伏或者模型空间形状的影响，实现带地形模型的计算，尤其适用于隧道或者矿井模型正演；

3.计算过程中考虑关断时间的影响，可以有效的模拟实际瞬变电磁探测；

4.对于后期瞬变电磁各向异性的研究具有重要指导意义。

附图说明

为了对本申请进行进一步的说明，提出一部分说明书附图加以说明，本申请的示意性实施例及其说明用于解释本申请，并不构成对本申请的不当限定。

图1为轴向电导率各向异性瞬变电磁三维三分量时域有限差分正演方法具体计算流程图；

图2为轴向电导率各向异性正演程序计算各项同性均匀半空间模型所得的dB_p/dt(p＝x,y,z)与B_p(p＝x,y,z)三分量响应与一维解析解计算结果对比图及精度误差图，其中，a.1为一维解析解和三维FDTD计算结果的dB_x/dt对比图；a.2为一维解析解和三维FDTD计算结果的dB_y/dt对比图；a.3为一维解析解和三维FDTD计算结果的dB_z/dt对比图；a.4为一维解析解和三维FDTD计算结果dB_p/dt(p＝x,y,z)精度误差图；b.1为一维解析解和三维FDTD计算结果的B_x对比图；b.2为一维解析解和三维FDTD计算结果的B_y对比图；b.3为一维解析解和三维FDTD计算结果的B_z对比图；b.4为一维解析解和三维FDTD计算结果B_p(p＝x,y,z)精度误差图；

图3为改变半空间电导率各向异性正演所得的dB_p/dt(p＝x,y,z)与B_p(p＝x,y,z)结果图；a.1)为半空间电导率各向异性的dB_x/dt；a.2为半空间电导率各向异性的dB_y/dt；a.3为半空间电导率各向异性的dB_z/dt；b.1为半空间电导率各向异性的B_x；b.2为半空间电导率各向异性的B_y；b.3为半空间电导率各向异性的B_z；

图4为改变层状地层电导率各向异性正演所得的dB_p/dt(p＝x,y,z)与B_p(p＝x,y,z)结果图；a.1为层状地层电导率各向异性的dB_x/dt；a.2为层状地层电导率各向异性的dB_y/dt；a.3为层状地层电导率各向异性的dB_z/dt；b.1为层状地层电导率各向异性的B_x；b.2为层状地层电导率各向异性的B_y；b.3)为层状地层电导率各向异性的B_z；

图5为改变异常体电导率各向异性正演所得的dB_p/dt(p＝x,y,z)与B_p(p＝x,y,z)结果图；a.1为异常体电导率各向异性的dB_x/dt；a.2为异常体电导率各向异性的dB_y/dt；a.3为异常体电导率各向异性的dB_z/dt；b.1为异常体电导率各向异性的B_x；b.2为异常体电导率各向异性的B_y；b.3为异常体电导率各向异性的B_z；

具体实施方式

下面通过实施例，并结合附图，对本发明的技术方案作进一步具体的说明。

首先对正演方法原理以及公式的推导过程进行一个说明及介绍。

本实施例引入轴向电导率张量，构建基于轴向电导率各项异性的电磁场控制方程，将电导率张量代入控制方程，采用Yee晶胞格式和直角坐标系进行网格离散，推导无源轴向电导率各向异性电场、磁场三分量迭代公式，以电流密度的形式加入电导率各向异性麦克斯韦方程组的安培环路定理中实现电流源的加载，计算过程中考虑关断时间的影响，可以不受地表起伏或者模型空间形状的影响，实现带地形模型的计算，尤其适用于隧道或者矿井模型正演。由于迭代格式中包含了源电流项，并且回线源瞬变电磁的激发源是细导线，实际建模中细导线的尺寸远小于晶胞尺寸，因而不能通过晶胞来模拟细导线。根据法拉第电磁感应定律和安培环路定理，有源所在单元的电磁场值通过积分求解得到。在具体程序中，电流密度等于电流值除以均匀剖分网格的面积。

下面具体描述本实施例的过程。

引入轴向各向异性电导率张量

式中，σ_x，σ_y，σ_z为对应方向上的轴向电导率。

各项异性、有耗、非磁性、无源媒质中的Maxwell方程组表示为：

式中E表示电场强度；H表示磁场强度；μ为磁导率；ε表示介电常数；t表示时间。

在电磁场空间激发长中存在位移电流和传导电流，公式(2)中(b)公式部分表示位移电流，表示传导电流。根据准静态近似条件，瞬变电磁勘探一般忽略位移电流，忽略位移电流后的麦克斯韦方程组为：

由于忽略了位移电流，方程(3)中的(b)式缺少电场对时间的导数，无法构成FDTD计算必需的显式时间步进格式。由于FDTD数值计算的需要，本实施例需要加入一项虚拟位移电流，加入虚拟位移电流后的麦克斯韦方程组为：

其中，γ具有介电常数的量纲，本文称其为虚拟介电常数，包含γ的项具有电流的量纲，称为虚拟位移电流。

结合各向异性介质中的本构关系，在直角坐标系下将公式(4)展开，得到公式(5)和(6)。

和

瞬变电磁勘探采用的是宽频带电磁场，并且低频电磁场是实现测深的主要部分，在进行三维正演时必须考虑低频电磁响应计算结果的可靠性。因此对公式(5)进行低频近似处理(可以参考该领域相关公知技术)得到公式(7)：

公式(6)-(7)为无源区域同轴各向异性电磁场计算的基本控制方程。

公式(6)和(7)均是通过公式(4)经过处理变换过来的，是进行轴向电导率各向异性瞬变电磁三分量三维FDTD正演的主要控制方程。其中公式(6)是在公式(4)的(a)式在直角坐标系下的展开形式，公式(7)是在公式(4)的(b)式在直角坐标系下展开并进行低频近似处理后的形式。

将轴向各向异性电导率张量代入控制方程(6)，采用Yee晶胞格式和直角坐标系对公式(6)-(7)进行网格离散，电场和磁场在空间和时间上的采样约定按照表1设置，在空间上，电场和磁场交替出现，空间离散采用后向差分，在时间轴上，电场和磁场交替采样，采样间隔为半个时间步，时间离散采用中心差分。

表1 Yee元胞中的电场和磁场分量的空间和时间采样约定

其中，i表示X方向上的网格节点，j表示Y方向上的网格节点，k表示Z方向上的网格节点，n表示时间迭代步。

采用Yee晶胞格式和直角坐标系对公式(6)-(7)进行网格离散，空间离散采用后向差分，时间离散采用中心差分。经过推导得到无源轴向各向异性时电场三分量迭代表达式如下：

i表示X方向网格节点数；j表示Y方向网格节点数；k表示Z方向网格节点数；n代表时间迭代步数；Δt表示时间步长；表示对应节点位置、对应坐标方向的轴向电导率；表示对应节点位置、对应坐标方向的电场；表示对应节点位置、对应坐标方向的磁场；Δx_i表示对应节点、对应方向的空间步长；其它分量含义相同，其它公式中字母定义与上述相同。

磁场求解迭代式与电导率无关，各向异性时磁场三分量表达式与各向同性时磁场三分量表达式相同，可以参考本领域的公知技术。

在有源媒质中，Maxwell方程组(4)中的方程(b)必须包含源电流项，修改为：

其中，J_s代表源电流密度

由于激励源电流加载在X轴和Y轴，Z轴上不添加电流，因而方程中仅存在J_sx和J_sy。

其中，J_sx和J_sy分别为电流密度在X方向和Y方向的分量。

因此有源区域电场水平分量迭代公式如下：

下面对本实施例的技术效果进行精度验证，并基于该程序对不同轴向电导率值对三分量响应的影响进行模拟。首先采用各项异性三维FDTD程序计算均匀半空间模型与一维解析解进行三分量对比并验证程序精度；然后分析半空间各向异性、层状地层各项异性、异常体各向异性等情况下瞬变电磁三分量响应特征及规律。

验证提出的轴向电导率各向异性瞬变电磁三分量正演方法及程序的可靠性，计算均匀半空间模型，将结果与一维解析解进行对比并验证精度。

选用发射电流1A、边长30m的回线源，以回线源中心点为坐标原点，以(7，7，0)为接收点，正演瞬变电磁场三分量磁场响应数据，其中发射电流波形采用梯形发射波形，发射电流波的上升沿时间10μs，持续时间10ms，下降沿时间10μs。如图2所示，为三维正演磁场dB_p/dt(p＝x,y,z)和B_p(p＝x,y,z)结果与一维解析解对比图和精度误差图。如图可以发现，各项异性程序计算结果与一维解析解结果基本一致，该方法满足精度要求。

选用发射电流1A、边长30m的回线源，以回线源中心点为坐标原点，以(7，7，0)为接收点，正演瞬变电磁场三分量磁场响应数据，其中发射电流波形采用梯形发射波形，发射电流波的上升沿时间10μs，持续时间10ms，下降沿时间10μs。改变半空间轴向电导率，模型参数如表2：

表2：改变半空间轴向电导率参数

图3分别是改变半空间电导率各向异性正演计算dB_p/dt(p＝x,y,z)和B_p(p＝x,y,z)结果图；图3(a.1-3)分别为dB_x/dt、dB_y/dt、dB_z/dt，图3(b.1-3)分别为B_x、B_y、B_z；通过图3可以发现，半空间轴向电导率对dB_p/dt(p＝x,y,z)和B_p(p＝x,y,z)的影响一致。以dB_p/dt(p＝x,y,z)为例，增大Z方向的轴向电导率，不影响dB_p/dt(p＝x,y,z)瞬变电磁响应；增大X或Y方向的轴向电导率，dB_p/dt(p＝x,y,z)瞬变电磁响应均增大，对于dB_x/dt，Y方向轴向电导率的影响大于X方向，对于dB_y/dt，X方向轴向电导率的影响大于Y方向，对于dB_z/dt，Y方向和X方向轴向电导率的影响一样；dB_p/dt(p＝x,y,z)和dB_p/dt(p＝x,y,z)规律一样。

选用发射电流1A、边长30m的回线源，以回线源中心点为坐标原点，以(7，7，0)为接收点，正演瞬变电磁场三分量磁场响应数据，其中发射电流波形采用梯形发射波形，发射电流波的上升沿时间10μs，持续时间10ms，下降沿时间10μs。改变层状地层轴向电导率，选用K(低高低)型模型，第一层地层厚度为40m，为各项同性介质，电导率为0.01(S.m)；第二层地层厚度为20m，为各向异性介质，模型参数如表3；第三层地层为各项同性介质，电导率为0.01(S.m)：

表3：改变层状地层轴向电导率参数

图4分别是改变层状底层电导率各向异性正演计算dB_p/dt(p＝x,y,z)和B_p(p＝x,y,z)结果图；图4(a.1-3)分别为dB_x/dt、dB_y/dt、dB_z/dt，图4(b.1-3)分别为B_x、B_y、B_z；通过图4可以发现，层状地层轴向电导率对dB_p/dt(p＝x,y,z)和B_p(p＝x,y,z)的影响一致。以dB_p/dt(p＝x,y,z)为例，增大Z方向的轴向电导率，不影响dB_p/dt(p＝x,y,z)瞬变电磁响应；增大X或Y方向的轴向电导率，dB_p/dt(p＝x,y,z)瞬变电磁响应均增大，对于dB_x/dt，Y方向轴向电导率的影响大于X方向，对于dB_y/dt，X方向轴向电导率的影响大于Y方向，对于dB_z/dt，Y方向和X方向轴向电导率的影响一样；dB_p/dt(p＝x,y,z)和B_p(p＝x,y,z)规律一样。

选用发射电流1A、边长30m的回线源，以回线源中心点为坐标原点，以(7，7，0)为接收点，正演瞬变电磁场三分量磁场响应数据，其中发射电流波形采用梯形发射波形，发射电流波的上升沿时间10μs，持续时间10ms，下降沿时间10μs。改变异常体轴向电导率，半空间为各向同性介质，电导率为0.01(S.m)，异常体的规模为(-26m～36m，-26m～36m，50m～100m)，异常体轴向电导率参数如表4：

表4：改变异常体轴向电导率参数

图5分别是改变异常体电导率各向异性正演计算dB_p/dt(p＝x,y,z)和B_p(p＝x,y,z)结果图；图5(a.1-3)分别为dB_x/dt、dB_y/dt、dB_z/dt，图5(b.1-3)分别为B_x、B_y、B_z；通过图5可以发现，异常体轴向电导率对dB_p/dt(p＝x,y,z)和B_p(p＝x,y,z)的影响一致。以dB_p/dt(p＝x,y,z)为例，增大Z方向的轴向电导率，不影响dB_p/dt(p＝x,y,z)瞬变电磁响应；增大X或Y方向的轴向电导率，dB_p/dt(p＝x,y,z)瞬变电磁响应均增大，对于dB_x/dt，改变Y方向轴向电导率使得dB_x/dt出现变号现象，其影响大于改变X方向轴向电导率，对于dB_y/dt，改变X方向轴向电导率使得dB_y/dt出现变号现象，其影响大于改变Y方向轴向电导率，对于dB_z/dt，Y方向和X方向轴向电导率的影响一样；dB_p/dt(p＝x,y,z)和B_p(p＝x,y,z)规律一样。

以上根据实施例对本发明做出了详细描述，然而，所描述是示例性的，本发明不仅仅限于实施例中，本领域技术人员完全能够根据本发明而对其做出各种形式的替换或者变更，在不脱离本发明宗旨和精神的前提下，凡是对本发明做出的各种变更及修饰均为本发明所涵盖的内容，均落入所附权力要求的范围之内。

Claims (8)

1.一种轴向电导率各向异性瞬变电磁三分量三维FDTD正演方法，其特征在于，包括：

控制方程构建步骤，构建基于轴向电导率各项异性的电磁场控制方程；

迭代公式推导步骤，将轴向电导率张量代入控制方程进而推导无源和有源轴向电导率各向异性电场三分量迭代公式；

响应特征计算步骤，利用迭代公式进行响应特征计算。

2.根据权利要求1所述的一种轴向电导率各向异性瞬变电磁三分量三维FDTD正演方法，其特征在于，控制方程构建步骤中，基于下式构建电磁场控制方程：

式中E表示电场强度；H表示磁场强度；ε表示介电常数；t表示时间；各向异性电导率张量

3.根据权利要求1所述的一种轴向电导率各向异性瞬变电磁三分量三维FDTD正演方法，其特征在于，所述迭代公式推导步骤包括：

引入的虚拟位移电流，将所述电磁场控制方程组变换到直角坐标系，将轴向各向异性电导率张量代入控制方程，推导无源轴向各向异性时电场三分量迭代表达式。

4.根据权利要求1所述的一种轴向电导率各向异性瞬变电磁三分量三维FDTD正演方法，其特征在于，所述迭代公式推导步骤中，各向异性电场三分量迭代公式为：

式中，i表示X方向网格节点数；j表示Y方向网格节点数；k表示Z方向网格节点数；n代表时间迭代步数；Δt表示时间步长；表示对应节点位置、对应坐标方向的轴向电导率；表示对应节点位置、对应坐标方向的电场；表示对应节点位置、对应坐标方向的磁场；Δx_i表示对应节点、对应方向的空间步长，其它分量含义与上述相同。

5.一种轴向电导率各向异性瞬变电磁三分量三维FDTD正演系统，其特征在于，包括：

控制方程构建模块，构建基于轴向电导率各项异性的电磁场控制方程；

迭代公式推导模块，将轴向电导率张量代入控制方程进而推导无源和有源轴向电导率各向异性电场三分量迭代公式；

响应特征计算模块，利用迭代公式进行响应特征计算。

6.根据权利要求5所述的一种轴向电导率各向异性瞬变电磁三分量三维FDTD正演系统，其特征在于，控制方程构建模块中，基于下式构建电磁场控制方程：

式中E表示电场强度；H表示磁场强度；ε表示介电常数；t表示时间；各向异性电导率张量

7.根据权利要求5所述的一种轴向电导率各向异性瞬变电磁三分量三维FDTD正演系统，其特征在于，所述迭代公式推导模块引入的虚拟位移电流，将所述电磁场控制方程组变换到直角坐标系，将轴向各向异性电导率张量代入控制方程，推导无源轴向各向异性时电场三分量迭代表达式。

8.根据权利要求5所述的一种轴向电导率各向异性瞬变电磁三分量三维FDTD正演系统，其特征在于，所述迭代公式推导模块中，各向异性电场三分量迭代公式为：

式中，i表示X方向网格节点数；j表示Y方向网格节点数；k表示Z方向网格节点数；n代表时间迭代步数；Δt表示时间步长；表示对应节点位置、对应坐标方向的轴向电导率；表示对应节点位置、对应坐标方向的电场；表示对应节点位置、对应坐标方向的磁场；Δx_i表示对应节点、对应方向的空间步长，其它分量含义与上述相同。