CN113779853B - 一种时域电磁感应-磁化效应分数阶三维数值模拟方法 - Google Patents

Abstract

本发明涉及一种时域电磁感应‑磁化效应分数阶三维数值模拟方法，通过将科尔‑科尔分数阶磁化率模型引入Maxwell方程组，采用频率域有理函数逼近算法和时间域卷积递归运算的方法，实现时间域磁感应强度的离散运算。改进了以Maxwell方程组为控制方程的电、磁场迭代过程，最后基于有限差分算法对控制方程各偏导项进行近似，并推导出电场和磁场各分量迭代方程。最终实现了时域电磁感应‑磁化效应分数阶三维数值模拟。本发明目的在于，可以克服目前研究方法仅能进行时域电磁感应‑磁化效应﹣1次幂率衰减的一维数值模拟，实现感应‑磁化效应分数次幂率衰减过程的三维数值模拟。

Description

一种时域电磁感应-磁化效应分数阶三维数值模拟方法

技术领域

本发明涉及一种地球物理勘探领域的时域电磁感应-磁化效应分数阶三维数值模拟方法。

背景技术

时间域电磁法(Time-Domain Electromagnetic Method，TEM)具有快速、经济、勘探深度大等优点，广泛应用于资源勘探和环境监测等领域。铁磁性矿物、磁性基岩，以及我国广泛分布的第四纪黄土、风化层及红土层具备磁化特征，时域电磁法在磁性地质条件下探测的过程中，在发射电流的on-time阶段，介质中的磁性粒子改变方向，获得磁化；发射电流关断后，磁性粒子逐渐失去磁化强度，电磁场观测系统将会检测到近似-1次幂律的衰减信号，这一过程被称为感应-磁化效应。时域电磁法在磁性地质环境的勘探中，若沿用传统的探测理论及方法，忽视磁化效应，将会导致电阻率解释过程中引入较大误差、无法对磁性矿产资源准确识别。

俄罗斯科学院石油地质与地球物理研究所的Kozhevnikov团队在地面时域电磁感应-磁化效应方面开展了系列研究，研究了长导线源激发下的均匀磁化介质感应-磁化效应数值模拟方法，和地面回线源激发下层状磁化介质的感应-磁化效应数值模拟方法,分析了层状结构对感应-磁化效应的影响(Antonov，2017；Kozhevnikov，2008，2011，2018)。Clark采用COMSOL软件对感应-磁化效应进行模拟，证明了土壤的磁化效应能够降低电磁传感器对地质目标体的检测能力(Clark，2014)。Sattel等基于深部硫化矿物和浅层SPM磁异常进行建模，研究了磁化效应和深部电性异常的电磁传播特征及识别方法(Sattel，2014；Sattel，2015)。Cowan等推导了地面大回线源的感应-磁化效应垂直和径向分量的解析式，并进行一维数值模拟，研究了感应响应与磁化响应交叉时间的估计方法(Cowan，2017)。近年来，皇家墨尔本理工大学Macnae团队基于Chikazumi磁化率模型，研究了航空电磁数据中磁化效应的评估方法，并应用于实测感应-磁化数据处理中，得到了更符合实际地质资料的解释结果。

以上证实了感应-磁化效应的研究能够提高时域电磁法在磁性地质中的勘探精度，现有模拟方法以基于Chikazumi磁化率模型的整数阶、一维模拟为主。然而，实际地质探测结果及实验室结果均表明，感应-磁化效应晚期斜率并不严格遵照Chikazumi模型的-1次幂律衰减，而是呈分数阶扩散规律，磁化率模型有待改进，引入分数阶微分算子。随着计算机性能提高及分数阶算子求解方法的不断拓展，需要实现感应-磁化效应数值模拟从一维、整数阶，向三维、分数阶跨越。因此，如何实现感应-磁化效应的分数阶三维数值模拟是待本领域技术人员迫切解决的一个技术问题。

发明内容

本发明所要解决的技术问题在于提供一种时域电磁感应-磁化效应分数阶三维数值模拟方法，通过将科尔-科尔分数阶磁化率模型引入Maxwell方程组，采用频率域有理函数逼近算法和时间域卷积递归运算的方法，实现时间域磁感应强度的离散运算。改进了以Maxwell方程组为控制方程的电、磁场迭代过程，最后基于有限差分算法对控制方程各偏导项进行近似，并推导出电场和磁场各分量迭代方程。最终实现了时域电磁感应-磁化效应分数阶三维数值模拟。本发明目的在于，可以克服目前研究方法仅能进行时域电磁感应-磁化效应﹣1次幂率衰减的一维数值模拟，实现感应-磁化效应分数次幂率衰减过程的三维数值模拟。

本发明是这样实现的，一种时域电磁感应-磁化效应分数阶三维数值模拟方法：

包括如下步骤：

1)、采用准静态条件下，无源Maxwell旋度方程作为电、磁场分量E_x、E_y、E_z、H_x、H_y的控制方程；

2)、引入科尔-科尔磁化率模型表征磁化介质磁导率变化过程，采用频率域有理函数逼近算法实现科尔-科尔磁化率模型的有理逼近；

3)、对磁化强度M求解公式进行频率域到时间域的转换，采用卷积递归运算的方法，完成时间域磁感应强度的离散运算；

4)、采用非均匀三维Yee氏网格对计算区域进行剖分，基于有限差分法，推导控制方程在整个计算区域的差分迭代格式；

5)、加载C-PML边界条件，实现电、磁场迭代运算；

6)、迭代计算结束后，提取磁场各分量响应，并对计算结果进行显示；

进一步地，步骤2)中，科尔-科尔磁化率模型表达式为：

式(1)中，ω表示角频率，i为虚部，χ_∞表示无穷频率下磁化率，χ₀表示零频磁化率，τ_c为时间常数，α为频散系数，取值在0与1之间，控制弛豫时间的分布；

对式(1)进行频率域有理函数逼近，可以得出：

其中c_j，d_k(j＝1，2，…n，k＝1，2，…m)为未知系数，m和n为常数，设置为n-m＝1；

通过对公式(2)进行误差极小化得到c_j，d_k(j＝1，2，…n，k＝1，2，…m)的最佳取值，利用近似有理函数χ(ω)具有不同的极点，通过对公式(2)转换到时间域，可以得出一系列指数函数的和：

进一步地，步骤3)中，磁化强度M时间域求解公式为：

其中H为磁场强度；

以公式(4)中求和公式第一项的计算为例，展开可得：

其中H的上角标表示磁场的取值时刻；

将公式(5)表达为前一时刻的卷积项，可得：

通过观察公式(5)及公式(6)可以得到卷积项的递归计算公式：

进一步地，步骤4)中，电场、磁场各分量的迭代顺序为(以x方向为例)：

其中E_x、B_x、H_x、M_x分别表示x方向的电场、磁感应强度、磁场及磁化强度，其上角标表示对应的时刻，其中ε为介电常数，μ₀为真空中导磁率，σ为电导率，Δt_j及At_j+1表示j及产1时刻的时间步长。

本发明与现有技术相比，有益效果在于：可以克服目前研究方法仅能进行时域电磁感应-磁化效应-1次幂率衰减的一维数值模拟，实现感应-磁化效应分数次幂率衰减过程的三维数值模拟。

附图说明

图1是一种时域电磁感应-磁化效应分数阶三维数值模拟方法示意图；

图2是磁化半空间模型下三维数值模拟模拟结果与解析解响应及误差对比；

图3不同磁化率参数下的感应电动势不同幂律衰减响应对比；

图4含有表面磁化层及低阻异常的感应电动势响应对比；

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

实施例

参见图1，一种时域电磁感应-磁化效应分数阶三维数值模拟方法，包括：

1)、采用散度方程作为磁场垂直分量的控制方程，准静态条件下、无源Maxwell旋度方程作为其他电磁场分量的控制方程结合非均匀Yee式网格和DuFort-Frankel方法构建显式差分格式，网格数目为107×107×53，其中x、y方向上的网格数目均为107个，z方向上网格数53个，最小和最大网格步长分别为10m和120m，除地面外计算区域其余5面最外层8个网格加载C-PML层。

2)在整个计算区域内设置电导率、磁化率参数、C-PML系数等参数，算例中背景电导率设置为0.01西门子/米，磁化率参数设置为χ_∞＝0,χ₀＝0.01,α＝0.5,τ_c＝1。

3)、将计算区域的电性参数代入到发射在空中的地下电场响应表达式中，计算初始时刻的x、y方向的电场响应。根据Maxwell旋度方程计算下半个时刻磁场响应。

4)、将参与迭代的CPU序列转化为GPU序列。

5)、将电场初始场值代入迭代方程中在整个计算区域更新电场值E_x、E_y、E_z，保存前一时刻与当前时刻电场值。

6)、将磁场初始场值代入迭代方程中在整个计算区域更新磁场值H_x、H_y、H_z，保存前一时刻与当前时刻磁场值。

7)、根据当前时刻的磁场值，计算磁化强度M_x、M_y、M_z，保存前一时刻与当前时刻磁化强度值。

8)根据当前时刻的磁场值及磁化强度M_x、M_y、M_z，计算当前时刻的感应电动势，即磁感应强度偏导项的值。

9)、是否完成全部迭代，如未完成，重复5-8步骤，如果完成全部迭代，将计算结果从GPU释放出来，并对计算结果进行显示，完成全部计算。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (1)

1.一种时域电磁感应-磁化效应分数阶三维数值模拟方法其特征在于，包括如下步骤：

1)、采用准静态条件下，无源Maxwell旋度方程作为电、磁场分量E_x、E_y、E_z、H_x、H_y的控制方程；

2)、引入科尔-科尔磁化率模型表征磁化介质磁导率变化过程，采用频率域有理函数逼近算法实现科尔-科尔磁化率模型的有理逼近；

3)、对磁化强度M求解公式进行频率域到时间域的转换，采用卷积递归运算的方法，完成时间域磁感应强度的离散运算；

4)、采用非均匀三维Yee氏网格对计算区域进行剖分，基于有限差分法，推导控制方程在整个计算区域的差分迭代格式；

5)、加载C-PML边界条件，实现电、磁场迭代运算；

6)、迭代计算结束后，提取磁场各分量响应，并对计算结果进行显示；

其中步骤2)中，科尔-科尔磁化率模型表达式为：

式(1)中，ω表示角频率，i为虚部，χ_∞表示无穷频率下磁化率，χ₀表示零频磁化率，τ_c为时间常数，α为频散系数，取值在0与1之间，控制弛豫时间的分布；

对式(1)进行频率域有理函数逼近，可以得出：

其中c_j,d_k(j＝1,2,…n,k＝1,2,…m)为未知系数,m和n为常数，设置为n-m＝1；

通过对公式(2)进行误差极小化得到c_j,d_k(j＝1,2,…n,k＝1,2,…m)的最佳取值，利用近似有理函数χ(ω)具有不同的极点，通过对公式(2)转换到时间域，可以得出一系列指数函数的和：

其中步骤3)中，磁化强度M时间域求解公式为:

其中H为磁场强度；

以公式(4)中求和公式第一项的计算为例，展开可得：

其中H的上角标表示磁场的取值时刻；

将公式(5)表达为前一时刻的卷积项，可得：

通过观察公式(5)及公式(6)可以得到卷积项的递归计算公式：

其中步骤4)中，电场、磁场各分量的迭代顺序为(以x方向为例)：

其中E_x、B_x、H_x、M_x分别表示x方向的电场、磁感应强度、磁场及磁化强度，其上角标表示对应的时刻，其中ε为介电常数，μ₀为真空中导磁率，σ为电导率，Δt_j及Δt_j+1表示j及j+1时刻的时间步长。