Three-dimensional forward modeling and response analysis of semi-airborne transient electromagnetic with arbitrary anisotropy based on finite volume method

XiaoTeng JIAGN; HuaiFeng SUN; Dong LIU; ShangBin LIU; Zhen WANG; XiaoSheng ZHOU; Jing Yang

doi:10.6038/pg2024HH0168

2024 , Vol. 39 >Issue 5: 2046 - 2058

DOI: https://doi.org/10.6038/pg2024HH0168

Three-dimensional forward modeling and response analysis of semi-airborne transient electromagnetic with arbitrary anisotropy based on finite volume method

XiaoTeng JIAGN ^,¹^,³ ,
HuaiFeng SUN ^,¹^,³^,⁴^,* ,
Dong LIU ¹^,²^,³ ,
ShangBin LIU ¹^,³^,⁴ ,
Zhen WANG ¹^,³^,⁵ ,
XiaoSheng ZHOU ¹^,² ,
Jing Yang ¹^,³

Expand

¹ Geotechnical and Structural Engineering Research Center, Shandong University, Jinan 250061, China
² Guangxi Communications Investment Group Co., Ltd., Nanning 530022, China
³ Laboratory of Earth Electromagnetic Exploration, Shandong University, Jinan 250061, China
⁴ Advanced Exploration and Transparent City Innovation Centre, Shandong Research Institute of Industrial Technology, Jinan 250061, China
⁵ Beijing Research Institute of Shandong University, Beijing 100086, China

Received date: 2023-04-24

Online published: 2024-12-19

Copyright

Fold

Abstract

Most research on transient electromagnetic phenomena in the field of airborne electromagnetics focuses on three-dimensional isotropic media, with limited studies on anisotropy. Interpreting data from regions with significant electrical anisotropy using isotropic models can lead to significant deviations from reality. To analyze the impact and extent of electrical anisotropy on transient electromagnetic signals from long conductive sources, this study employs a finite volume-based three-dimensional forward modeling algorithm for arbitrary anisotropic media. By solving the discretized full tensor conductivity time-domain electromagnetic field equations, a methodology is developed. First, based on Maxwell's equations, and employing the background and anomaly separation theory, the computational domain is divided into a homogeneous half-space and an anomalous body. The electromagnetic response of the homogeneous half-space serves as the background field, while the electromagnetic response of the anomalous body forms the anomaly field. The curl equation of the electric field anomaly is further curled to obtain the curl equation for the secondary electric field. Integrating the curl equation for the secondary electric field within a control volume yields an integral equation. Applying Gauss's theorem converts the volume integral equation into a surface integral equation, which is then discretized for independent computation on the surface. Employing backward Euler algorithm in the time domain leads to the final discrete equation.To enhance accuracy and efficiency of time-domain electromagnetic field solution, the algorithm combines time segmentation with direct equation solving. The algorithm's validity is verified through one-dimensional anisotropic layered model calculations. Subsequently, the impact of principal axis anisotropy on attenuation voltage is simulated and analyzed. Arbitrary anisotropic models are then simulated and the effects of anisotropic dip and azimuth angles on the response are investigated. Finally, the algorithm's simulation capabilities for complex geological models are tested. Experimental results demonstrate the algorithm's effectiveness in three-dimensional arbitrary anisotropic forward modeling of transient electromagnetic phenomena, offering high precision in simulation outcomes.

Key words： Semi-airborne transient electromagnetic; Limited volume; Anisotropy; 3D forward modeling

Cite this article

XiaoTeng JIAGN , HuaiFeng SUN , Dong LIU , ShangBin LIU , Zhen WANG , XiaoSheng ZHOU , Jing Yang . Three-dimensional forward modeling and response analysis of semi-airborne transient electromagnetic with arbitrary anisotropy based on finite volume method[J]. Progress in Geophysics, 2024 , 39(5) : 2046 -2058 . DOI: 10.6038/pg2024HH0168

0 引言

半航空瞬变电磁法是通过地面长导线源发射信号，在飞行器上使用接收线圈接收感应电磁场的地球物理勘探方法.相对于地面采集系统，半航空瞬变电磁法采用空中接收方式具有作业效率高的优点；相对于航空系统，该方法由于采用地面发射，具有数据信噪比高、勘探深度大的优势.半航空瞬变电磁法现已广泛应用于矿产勘查、地下水探勘等诸多地球物理探测领域(Yu et al., 1997；Fountain, 1998；Weidelt，1999；闫述等，2002；Commer and Newman, 2004；王昌学等，2006；吕国印，2007；薛国强等，2007；沈金松和郭乃川，2008；Streich，2009；孙怀凤等，2013；Jahandari and Farquharson, 2014；彭荣华等，2016；任秀艳，2018；刘亚军等，2019；王振水等，2021；杨静等，2022；李雪峰等，2023；Sun et al., 2023).

在早期的研究中，大多将地下介质认为是各向同性介质，然而，随着理论和实践技术的进步，发现电各向异性现象广泛存在于地下介质中.因此，一些学者对电磁法的各向异性开展了研究.刘云鹤等(2018)系统阐述了电磁法各向异性的国内外研究现状.徐世浙等(1983)基于有限单元法完成了二维大地电磁主轴各向异性正演计算；Osella和Martinelli(1993)研究了各向异性对大地电磁正演模拟的影响；陈桂波等(2009)利用积分方程法计算了三轴各向异性地层中大地电磁响应；林长佑等(2004)对采集的大地电磁实测数据进行了一维各向异性反演与分析；杨长福等(2005)、Pek和Santos(2006)都对大地电磁数据进行了各向异性反演计算；喻国等(2021)完成对目标函数及其梯度的计算，借助非线性共轭梯度优化技术，实现了对大地电磁各向异性介质的反演.Um等(2012)开发并验证了一种具有自适应时间步长的迭代有限元时域方法，用于模拟三维各向异性地球介质中的瞬态电磁场；Yin等(2016)在Um的研究方法基础上，基于有限单元法实现了时间域航空电磁法的三维各向异性研究，并分析了各向异性对响应的影响；周建美等(2014)基于有限体积法研究了瞬变电磁双轴各向异性正演并分析了电各向异性对响应的影响，发现水平方向上的电导率比垂直方向电导率对响应影响更大.

然而，针对半航空瞬变电磁法的研究仍主要集中在各向同性介质，对于地下不良地质体的宏观电各向异性的半航空瞬变电磁的研究还未见报道.本文从有限体积理论出发，对时域Maxwll方程组进行空间离散，利用隐式差分格式的后退欧拉法进行时间域离散，实现了半航空瞬变电磁三维任意各向异性正演.为提高计算效率，本文采用非均匀化网格与时间分段等步长相结合的方式，减少了网格数量与求解次数.最后，通过对主轴各向异性、倾斜各向异性以及任意各向异性模型进行模拟计算，对比分析了不同电各向异性模型对瞬变电磁数据的影响.

1 三维正演

1.1 理论推导

基于Maxwell方程组，采用背景场、异常场分离理论，将计算区域分解为均匀半空间和异常体，将均匀半空间的电磁响应作为背景场，而异常体的电磁响应为异常场.因此，电、磁场由以下两部分组成：

其中，E_b、H_b分别为电场背景场和磁场背景场；E_s和H_s分别为电场异常场和磁场异常场.背景场的计算首先需要求解均匀半空间模型的频率域响应，之后通过正、余弦变换将频率域的响应转换到时间域的响应.将 E_b、H_b带入麦克斯韦方程组，并与总场的麦克斯韦方程组做差，可以得到 E_s和H_s旋度方程，即二次电磁场方程，对 E_s旋度方程再取旋度，并把 H_s旋度方程代入，就得到二次电场的双旋度方程：

对于电导率任意各向异性介质，

为3×3对称张量：

可以通过对主轴各向异性电导率张量做图 1所示的欧拉旋转得到.首先将原始坐标系绕z轴逆时针旋转角度α_s得到一个新坐标系，然后将其绕新坐标系的x轴逆时针旋转角度α_d得到中间坐标系，最后绕中间坐标系的z轴逆时针旋转角度α_l得到最终坐标系.其中，α_s、α_d、a_l三个角度分别为各向异性走向角、各向异性倾角、各向异性偏角，每次旋转都可以得到对应角度的旋转矩阵，将旋转矩阵与主轴各向异性电导率张量相乘就可以得到

显示原图|下载原图ZIP|生成PPT

图1 电导率张量形式欧拉旋转示意图

Fig 1 Euler rotation in conductivity tensor form

1.2 方程离散

有限体积法基本思想是将计算区域划分为一系列不重复的控制体积，将待解的偏微分方程在控制体积内积分得到积分方程，然后离散积分方程得到线性方程组进行求解.本文采用Yee(1966)网格对计算域进行剖分，并基于Yee网格划分控制体积对二次电场的方程进行空间离散，以任意一点E_y为例对控制体积进行说明：如图 2所示，图中深色区域表示E_y的控制体积，由4个相邻的Yee网格组成，控制体积的顶点是Yee网格表面的形心，控制体积的棱边长度、表面面积以及体积可由对应的Yee的尺寸计算得到.

显示原图|下载原图ZIP|生成PPT

图2 y方向控制体积

Fig 2 Control volume of y axis

将二次电场的双旋度方程(3)在控制体积V内进行积分，得到积分方程为：

根据高斯公式，公式(5)中的双旋度体积分可以转化为二次电场旋度与表面外法向量叉乘的面积分，然后将面积分项离散到控制体积的六个表面独立运算.对于积分方程中的体积分项，电场的变化率与空间无关，可以将体积分项分解为电导率的体积分独立计算.但是，在各向异性介质中，电流密度的计算需要考虑介质在不同方向上的电导率，因此电流密度在每个方向上的大小不仅取决于该方向上的电场分量，还受到垂直于该方向的另外两个方向上电场分量的影响，每个方向上的电流密度计算公式为：

因此，首先将电流密度离散到三个方向，再将每个方向的体积分项分解为电导率的体积分独立计算.以y方向的电导率体积分为例：

完成空间离散后，采用后退欧拉法对控制方程在时间域进行离散，获得的离散方程以E_y为例：

所有的离散方程组成一个大型的方程组，将其写为矩阵形式为：

其中C和D为不同时刻空间任意一点二次电场系数矩阵，R为由背景电导率等相关参数构成的系数矩阵.

1.3 初始条件与边界条件

初始和边界条件是求解微分方程的关键，不同的初始条件和边界条件决定了微分方程不同解.本文中，我们采用了阶跃波形作为电磁正演模拟的发射源激励信号，并沿其上升沿激励.在t=0时刻，不存在发射电流，电磁场的初始条件为：

电磁波会随着距离的增加而逐渐减弱，并且在远离发射源的地方会变得非常微弱，甚至不可测量.因此，在计算电磁场时，我们需要在计算区域的外边界上设置一个适当的边界条件，以模拟这种衰减.在本文中，我们选择了Dirichlet边界条件模拟电磁场的随距离的衰减，电场切向分量在计算区域外边界处为零，即：

1.4 网格剖分与时间步长

随时间迭代求解方程组(9)是算法最耗时部分，因此选择合理的网格剖分方式以及迭代步长对于提高计算速度至关重要.方程组(9)可简化为：

其中A为系数矩阵，与模型网格剖分、介质电导率性质以及时间步长有关.为提高计算速度，改进数值稳定性，本文在网格划分时采用下述方式：发射源附近采用均匀网格，然后按固定比例因子向边界逐渐扩大，同时限定最大网格尺寸，保障数值稳定，这种剖分方式大大减少了网格数量，提高了计算效率.一般来讲，时间步长越小计算精度越高，但这导致正演时间越长，因此在满足计算精度的条件下应选择较大的时间步长.为使得数值模拟更加准确和稳定，本文采用时间分段等步长，具体而言，将整个模拟时间区间分成若干个等长的时间段，并在每个时间段内采用相同的时间步长来进行模拟，这样不需在同时间段内重复计算系数矩阵A，只需更新背景场的信息，提高了正演效率，也有效地保证模拟的精度和稳定性.另外，采用时间分段等步长的方法还可以很方便地控制模拟的时间范围和总步数，因为每个时间段的长度相等，所以只需要控制时间段的个数就可以控制模拟的总时间范围，而总步数也可以直接计算得到.

2 主轴各向异性三维正演精度验证与响应规律分析

2.1 精度验证

为验证本文半航空瞬变电磁三维各向异性正演程序对各向异性模型正演的准确性及精度，考虑具有准解析解的一维横向各向同性(VTI)模型.电阻率参数及层厚如图 3所示，其中第一层与第三层为各向同性，第二层为各向异性层，接收点坐标为(0, 100, 110).此模型的准解析解采用正余弦变换算法，由频率域转换而来，一维频率域各向异性响应采用软件EMmod计算.为验证本文开发算法的正确性，采用有限元数值模拟方法的对上述模型正演计算，计算结果如图 4所示.

显示原图|下载原图ZIP|生成PPT

图3 一维VTI模型示意图

Fig 3 Schematic diagram of one-dimensional VTI model

显示原图|下载原图ZIP|生成PPT

图4 一维层状VTI模型数值解与解析解对比

Fig 4 Comparison between numerical and analytical solutions of one-dimensional layered VTI model

图 4给出了三个不同测点处模型的电动势衰减曲线与准解析解以及有限元计算结果的对比，从图中可以看出，在整个观测时间范围内，数值解与准解析解、有限元结果吻合较好，满足计算精度要求.

2.2 层状主轴各向异性模型

利用本文开发的半航空瞬变电磁三维各向异性正演程序进一步研究主轴各向异性对半航空瞬变电磁响应的影响.一维主轴各向异性模型如图 5所示.其中第二层电阻率值参考表 1.

显示原图|下载原图ZIP|生成PPT

图5 层状主轴各向异性模型

Fig 5 Layered principal axis anisotropic model

表1 不同模型电阻率

Table 1 Resistivity of different models

	model 0	model 1	model 2	model 3
ρ_x(Ω·m)	50	10	50	50
ρ_y(Ω·m)	50	50	10	50
ρ_z(Ω·m)	50	50	50	10

图 6给出了层状主轴各向异性模型在不同接收点(0, 30, 110)、(0, 100, 110)、(0, 300, 110)、(300, 300, 110)处的计算结果.从图中可得：在早期时间道内，层状主轴各向异性模型与各向同性模型的脉冲响应值差别很小，且衰减趋势一致，说明主轴各向异性在早期对脉冲响应的影响较小.在偏移距较小时，主轴各向异性的脉冲响应值略大于各向同性模型，在偏移距最大时，呈相反规律.但无论偏移距如何变化，主轴各向异性模型的脉冲响应值在10^-5~10^-4 s时间道的计算结果存在明显拐点，这是由于模型第一层和第三层为各向同性而第二层为各向异性，导致电磁波在不同方向上传播的速度和路径长度不同，从而电磁波到达接收点的时间存在差异，这就引起了脉冲响应存在不同的幅度和时间延迟，从而导致拐点的出现.在晚期时间道，层状主轴各向异性模型的脉冲响应逐渐高于层状各向同性模型，且响应差距随着偏移距增加而逐渐变大.在偏移距较小时，z主轴各向异性的层状模型在晚期的衰减电压明显小于x轴跟y轴，说明z方向电导率各向异性影响程度大于水平方向电导率.

显示原图|下载原图ZIP|生成PPT

图6 层状主轴各向异性模型计算结果

Fig 6 Calculation results of low resistance anisotropic layered model

2.3 三维各向异性模型

在一维模型的基础上，分析三维主轴各向异性对半航空瞬变电磁响应的影响.三维主轴各向异性倾斜低阻异常体模型如图 7所示，其异常体电阻率值电阻率值参考表 1.

显示原图|下载原图ZIP|生成PPT

图7 三维各向异性模型

Fig 7 3D anisotropic model

图 8给出了三维各向异性模型在不同接收点(0, 30, 110)、(0, 100, 110)、(0, 300, 110)、(300, 300, 110)处的计算结果.从图中可得：在早期时间道内，三维各向异性模型与各向同性模型的脉冲响应值几乎相同，且衰减趋势一致，说明主轴各向异性在早期对脉冲响应的影响较小.；且不存在层状主轴各向异性模型在10^-5~10^-4 s时间道内的衰减电压拐点，进一步验证了拐点是由于首层为各向同性第二层为各向异性导致的.在偏移距最小时，model 3在晚期时间道与各向同性模型的计算结果差异明显，说明z方向电导率各向异性在偏移距较小时对脉冲响应的影响大于水平方向，但随着偏移距的增大，与各向同性模型的计算结果差异变小，是由于异常体体积较小所致.

显示原图|下载原图ZIP|生成PPT

图8 三维各向异性模型各向异性计算结果

Fig 8 Calculation results of anisotropy for 3D anisotropic model

3 任意各向异性三维正演响应规律分析

3.1 任意各向异性层状模型响应分析

3.1.1 各向异性倾角

在一维层状主轴各向异性模型model 1的基础上进一步分析任意各向异性对半航空瞬变电磁响应的影响.在之前模型的基础上，取4个不同数值的各向异性倾角，角度数值分别为0°、30°、60°、90°, 设定各向异性偏角与各向异性走向角均为零，分析各向异性倾角对半航空瞬变电磁响应的影响.

图 9给出了各向异性偏角与各向异性走向角均为零，各向异性倾角为4个不同的数值时的层状模型在不同接收点处的计算结果.从图中可得：各向异性倾角为0°和90°的模型的脉冲响应曲线几乎重合，且与各向异性倾角为30°和60°的模型的脉冲响应曲线差异较大，这是由于各向异性倾角为90°时同样为主轴各向异性模型，与model 0是一致的，而另外两个模型为任意各向异性模型.在偏移距较小时，各向异性倾角的不同值对早期脉冲响应无明显影响，曲线基本重合；在晚期，任意各向异性模型的脉冲响应值大于主轴各向异性模型，且脉冲响应的值上升，这是因为任意各向异性介质中，电磁波在传播过程中会受到路径的影响，从而影响衰减电压的分布，当电磁波沿着电导率变化较大的路径传播时，衰减电压会增加，而当电磁波沿着电导率变化较小的路径传播时，衰减电压会减小.四个模型在10^-5~10^-4 s时间道内，脉冲响应同样存在明显拐点，是由于第一层和第三层为各向同性介质、第二层为各向异性介质造成的.

显示原图|下载原图ZIP|生成PPT

图9 层状任意各向异性模型计算结果

Fig 9 Calculation results of arbitrary anisotropic layered model

3.1.2 各向异性走向角

本节分析各向异性走向角对半航空瞬变电磁信号的影响，取4个不同数值的各向异性走向角，角度数值分别为0°、30°、60°、90°.

图 10给出了各向异性偏角与各向异性倾角均为零，各向异性走向角为4个不同的数值时层状模型在不同接收点(0, 30, 110)、(0, 100, 110)、(0, 300, 110)、(300, 300, 110)处的计算结果.从图中可得：各向异性走向角为0°和90°的模型均为主轴各向异性模型，脉冲响应曲线几乎重合，30°和60°的模型为任意各向异性模型，30°和60°的模型的脉冲响应值与0°和90°的模型的脉冲响应差异明显比各向异性倾角的差异小，说明不同各向异性走向角对脉冲响应的影响小于不同各向异性倾角.不同偏移距下的各向异性走向角的不同值对早期脉冲响应无明显影响，曲线基本重合；在晚期，任意各向异性模型的脉冲响应值小于主轴各向异性模型.四个模型在10^-5~10^-4 s时间道内，脉冲响应同样存在明显拐点，是由于第一层和第三层为各向同性介质、第二层为各向异性介质造成的.

显示原图|下载原图ZIP|生成PPT

图10 任意各向异性层状模型计算结果

Fig 10 Calculation results of arbitrary anisotropic layered model

3.2 各向同性背景中的三维各向异性模型

考虑到实际异常体可能为明显电各向异性的矿体，但周围介质可能为各向同性介质，因此建立如图 11所示模型，层状介质均为各向同性，异常体为任意各向异性，其电阻率值参考表 1，将各向异性倾角和走向角取四组不同数值，分别为0°、15°、30°和45°.

显示原图|下载原图ZIP|生成PPT

图11 异常体模型

Fig 11 Abnormal model

图 12给出了各向异性走向角与各向异性倾角分别为4个不同的数值、各向异性偏角为零的地质模型在不同接收点处的计算结果.从图中可得：在早期时间道内，不同旋转角度的地质模型的脉冲响应值差别很小，且衰减趋势一致，说明不同旋转角度在早期对脉冲响应的影响较小.在偏移距较小时，model 1与model 2在晚期10^-3~10^-2 s时间道内与主轴各向异性模型的计算结果有明显差异，说明旋转角度在后期对脉冲响应影响较大；在偏移距较大时，不同角度的地质模型的脉冲响应无明显差异，是由于异常体体积较小所致.

显示原图|下载原图ZIP|生成PPT

图12 任意各向异性模型计算结果

Fig 12 Calculation results of arbitrary anisotropy model

3.3 各向异性背景中的三维各向同性模型

实际工程中低阻异常体是我们关心区域，为检验本程序对电各向异性介质中的低阻异常的响应情况，建立如图 13所示模型.其中第二层介质为各向异性，第一层、第三层以及含水体为各向同性，含水体电阻率值参考表 1，将各向异性倾角和走向角取四组不同数值，分别为0°、15°、30°和45°.

显示原图|下载原图ZIP|生成PPT

图13 异常体模型

Fig 13 Abnormal model

图 14给出了各向异性走向角与各向异性倾角分别为4个不同的数值、各向异性偏角为零的地质模型在不同接收点(0, 30, 110)、(0, 100, 110)、(0, 300, 110)、(300, 300, 110)处的计算结果.从图中可得：在10^-5~10^-4 s时间道内计算结果存在明显拐点，整体响应与层状主轴各向异性响应存在明显差异，但不同旋转角度的模型的脉冲响应值差别很小，且衰减趋势一致，推测是由于异常体体积较小或旋转角度较小所致.

显示原图|下载原图ZIP|生成PPT

图14 任意各向异性模型计算结果

Fig 14 Calculation results of arbitrary anisotropy model

4 结论

本文基于有限体积法对电磁场扩散方程进行空间离散，结合后退欧拉法进行时间域离散，实现了接地源半航空瞬变电磁有限体积三维任意各向异性正演，并结合pardiso求解器直接求解，通过一维VTI模型验证了该算法的精度和有效性.由主轴各向异性层状模型的模拟结果可知，响应会由于层状介质的电导率差异而产生拐点，而且z方向电导率影响程度大于水平电导率；三维主轴各向异性异常体由于设定尺寸较小，无明显规律.通过对任意各向异性的模拟计算结果表明，各向异性倾角对响应的影响较大，而各向异性走向角影响较小.进一步通过对复杂模型的计算，验证了开发的接地源半航空瞬变电磁有限体积三维任意各向异性正演程序的有效性，为半航空瞬变电磁三维正演算法提供了参考，为进一步计算带地形的复杂模型以及实测数据的反演提供正演算法基础.

感谢审稿专家提出的修改意见和编辑部的大力支持！

References

Publishing order | Descend order by publishing year | Descend order by cited within

Chen

G B

, Wang

H N

, Yao

J J

. Modeling of electromagnetic responses of 3-D electrical anomalous body in a layered anisotropic earth using integral equations. Chinese Journal of Geophysics, 2009, 52(8 2174 2181.

模态框（Modal）标题

Abstract

Cite this article

0 引言

1 三维正演

1.1 理论推导

图1 电导率张量形式欧拉旋转示意图

1.2 方程离散

图2 y方向控制体积

1.3 初始条件与边界条件

1.4 网格剖分与时间步长

2 主轴各向异性三维正演精度验证与响应规律分析

2.1 精度验证

图3 一维VTI模型示意图

图4 一维层状VTI模型数值解与解析解对比

2.2 层状主轴各向异性模型

图5 层状主轴各向异性模型

表1 不同模型电阻率

图6 层状主轴各向异性模型计算结果

2.3 三维各向异性模型

图7 三维各向异性模型

图8 三维各向异性模型各向异性计算结果

3 任意各向异性三维正演响应规律分析

3.1 任意各向异性层状模型响应分析

3.1.1 各向异性倾角

图9 层状任意各向异性模型计算结果

3.1.2 各向异性走向角

图10 任意各向异性层状模型计算结果

3.2 各向同性背景中的三维各向异性模型

图11 异常体模型

图12 任意各向异性模型计算结果

3.3 各向异性背景中的三维各向同性模型

图13 异常体模型

图14 任意各向异性模型计算结果

4 结论

References