一种基于Chikazumi模型的感应-磁化效应三维数值模拟方法与流程

一种基于chikazumi模型的感应
‑
磁化效应三维数值模拟方法
技术领域
1.本发明涉及一种地球物理勘探领域的基于chikazumi模型的感应
‑
磁化效应三维数值模拟方法。


背景技术：

2.航空瞬变电磁法(airbrone electromagnetic method，aem)具有快速、高效、探测范围广等优点，已被广泛应用于矿产资源探勘查，环境监测，基础地质调查等多个领域(lei，2006；zhang，2005)。瞬变电磁法在探测过程中，发射电流的on
‑
time阶段会让地下的磁性粒子改变方向，获得磁化；在电流关断后，磁性粒子会失去磁化强度，接收系统会观测到近似
‑
1次幂率的衰减信号，这一过程被称为感应
‑
磁化效应。国内外实验结果表明，感应
‑
磁化效应不仅存在于多种的岩石及矿物中(buselli，1982；neumann，2006)。在0.1
‑
15um半径范围内的磁铁矿和赤铁矿颗粒也可产生感应
‑
磁化效应(bournas，2017；dabas and skinner，2002)。随着机载发射机偶极矩的增加以及系统噪声水平的数量级的降低，电磁传播过程中的感应
‑
磁化效应已可以被正常检测到(kwan et al.，2015)。目前科学家们在多地均检测到感应
‑
磁化效应的存在(macnae，2017；zhang，2020)。而传统的探测理论，往往会忽略感应
‑
磁化效应的影响，引入较大的误差，无法对探测信息做出正确的解释。因此，对感应
‑
磁化效应进行有效观测和准确解释尤为重要。
3.chikazumi根据观测到的基本衰减规律，推导出了被广泛应用的chikazumi磁化率模型(chikazumi，1978)。lee利用chikazumi模型进行理论计算，验证了基于chikazumi磁化率模型模拟出的电磁响应呈
‑
1次幂率衰减，贴合实际观测到的情况(lee，1984)。kamnev等通过观测磁性岩石样品，确定了响应晚期呈
‑
1.07～
‑
1.13幂率衰减(kamnev，2015)。kozhevnikov等针对地面长导线源激发的地下均匀磁化介质与地面回线源激发的层状磁化介质的数值模拟进行了研究，分析了零频磁化率等性质对于感应
‑
磁化效应的影响(antonov，2017；kozhevnikov，2008，2011，2018)。cowan等推导了地面大回线源下感应
‑
磁化效应垂直和径向分量的解析式，研究有关感应响应与磁化响应的交叉时间的估计方法(cowan，2017)。zhang等研究了长导线源下的感应
‑
磁化效应，进行了层状模型下的一维以及二维数值模拟。并研究了岩、矿石样品的磁化率，验证了磁化效应对于铁磁性矿物探测的有效性(zhang，2020)。然而，由于基于磁化率随时间变化的磁化率模型直接进行时域三维数值模拟存在困难，目前相关正演研究多数集中于一维数值模拟方面。
4.目前，常用的瞬变电磁三维数值模拟方法，主要有有限差分法(wang，1993)、有限元法(assous，1993)以及无网格法(yang，2014)等。其中，有限差分方法是数值模拟最早采用的方法，与其他方法相比，虽然在计算复杂目标体会产生较大误差，但是较为直观，且具有较强的实用性。由于chikazumi磁化率模型在时域的非线性，无法直接在时域进行运算。若采用频率域有限差分方法，由于一般的频率域有限差分方法使用真空磁导率，所以无需考虑磁导率的位置对应问题，在把磁导率变成复数模型直接进行运算时，需要将磁化率展开成矩阵进行除法运算，会导致交错网格位置对应错误以及计算困难，得到错误结果。因
此，如何实现感应
‑
磁化效应的三维数值模拟是待本领域技术人员迫切解决的一个技术问题。


技术实现要素：

5.本发明所要解决的技术问题在于提供一种基于chikazumi模型的感应
‑
磁化效应三维数值模拟方法，针对磁化率模型在时域差分离散困难，无法直接进行时域三维数值模拟的问题，将chikazumi复磁化率模型引入频域maxwell方程组，磁场响含有复磁导率与复频率乘积的形式，传统频域有限差分方法无法求解复磁导率maxwell方程组，提出了复磁导率的矩阵运算形式，优化控制方程的计算方式，根据差分离散后磁场在网格上的位置关系，重新建立了磁导率和磁场响应的科学映射矩阵，定义了磁异常体特别的边界设置方式，实现了利用频域有限差分方法对感应
‑
磁化效应的三维求解。并利用正余弦变换数值滤波算法，将频域结果变换到时域，最终实现了感应
‑
磁化效应时域三维数值模拟。本发明是这样实现的，一种基于chikazumi模型的感应
‑
磁化效应三维数值模拟方法：
6.包括如下步骤：
7.1)、采用频率域加源maxwell旋度方程作为电、磁场分量e
x
、e
y
、e
z
、h
x
、h
y
、h
z
的控制方程；
8.2)、引入chikazumi磁化率模型表征磁化介质磁导率变化过程，建立磁导率分量μ
x
(ω)、μ
y
(ω)、μ
z
(ω)与磁场响应h
x
、h
y
、h
z
的映射矩阵，并带入控制方程；
9.3)、采用非均匀三维yee氏网格对计算区域进行剖分，基于频域有限差分法，推导控制方程在整个计算区域的差分格式，并将其改写为离散项；
10.4)、加载dirichlet边界条件，将磁导率分量与磁场分量进行元素相乘后，保证磁场控制方程的优先运算，将其带入电场控制方程中；
11.5)、将控制方程整理成ax＝b的形式，通过计算x＝a\b得到磁场各分量频率域响应结果；
12.6)、利用正余弦变换数值滤波算法，将磁场响应变换到时域，并对计算结果进行显示；
13.其中步骤2)中，chikazumi磁化率模型表达式为：
[0014][0015]
式(1)中，ω表示角频率，i为虚部，χ0表示零频磁化率，τ1、τ2为时间常数τ的最小值和最大值；
[0016]
对式(1)的复磁化率模型带入控制方程中，可以得出：
[0017][0018][0019]
其中e
x
、e
y
、e
z
为电场强度e在x、y、z三个方向上的分量，h
x
、h
y
、h
z
为磁场强度h在x、y、z三个方向上的分量，j
x
、j
y
、j
z
为导电电流密度j在x、y、z三个方向上的电流源密度，σ
x
、σ
y
、σ
z
代表电导率σ在x，y，z三个方向上的分量；μ
x
(ω)、μ
y
(ω)、μ
z
(ω)代表磁导率μ(ω)在x，y，z三个方向上的分量；
[0020]
将差分矩阵进行离散，并加载dirichlet边界条件，将磁场控制方程带入电场控制方程中，可以得到步骤4)中的运算矩阵为：
[0021][0022]
其中是电场强度e的差分矩阵，代表电场分量e
x
、e
y
、e
z
在不同方向上的微分形式；是磁场h的差分矩阵，代表磁场分量h
x
、h
y
、h
z
在不同方向上的微分形式；对于回线源j，只需要对j
x
、j
y
进行赋值，j
z
分量总为零；将μ
x
(ω)、μ
y
(ω)、μ
z
(ω)根据差分离散后的磁场分量建立映射矩阵，并带入控制方程进行运算；对于磁异常体，网格平面上的磁导率可以用两个相邻立方体的磁导率做面平均进行替代；
[0023]
其中步骤5)中，将式(4)整理成ax＝b的形式，可得到：
[0024][0025][0026][0027]
其中，a
11
、a
22
、a
33
分别为：
[0028][0029]
将计算结果代入步骤6)中的正余弦变换数值滤波算法，可得到感应
‑
磁化响应的时域解。
[0030]
本发明与现有技术相比，有益效果在于：可以克服基于chikazumi磁化率模型感应
‑
磁化效应三维数值模拟困难的问题，实现基于频域有限差分方法的感应
‑
磁化效应衰减过程的三维数值模拟。
附图说明
[0031]
图1是一种基于chikazumi模型的感应
‑
磁化效应三维数值模拟方法示意图；
[0032]
图2是磁化半空间模型下三维数值模拟模拟结果与解析解响应及误差对比；
[0033]
图3复杂磁化模型下，三维数值模拟感应电动势水平切片图；
[0034]
图4复杂磁化模型下，三维数值模拟感应电动势垂直切片图；
具体实施方式
[0035]
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
[0036]
实施例
[0037]
参见图1，一种基于chikazumi模型的感应
‑
磁化效应三维数值模拟方法，包括：
[0038]
1)、采用频域maxwell旋度方程作为磁场分量的控制方程，准静态条件下、结合非均匀yee式网格和dirichlet边界条件，构建网格差分格式，并且网格数目为38
×
38
×
32，其中x、y方向上的网格数目均为38个，z方向上网格数32个，相邻网格的比例系数为2，最小和最大网格步长分别为3m和2560m。
[0039]
2)在整个计算区域内设置电导率、零频磁化率等参数，算例中背景电导率设置为0.01西门子/米，分别在距离地面10m和40m的深度设置一小一大两个磁异常体，分别位于发射线圈的两侧，异常体1边长为12m
×
12m
×
10m，异常体2边长为36m
×
36m
×
40m，设置异常体零频磁化率χ0＝0.01，τ1＝10
‑6，τ1＝106。在计算过程中，网格边缘的电导率可以用相邻的四个立方体的电导率做体平均进行替代，网格平面上的磁导率可以用两个相邻立方体的磁导率做面平均进行替代。
[0040]
3)、根据磁场响应差分离散后网格上的位置关系，建立磁导率和磁场响应的科学映射矩阵。
[0041]
4)、在10
‑1hz
‑
106hz范围内等间隔选取70个频点进行计算，利用控制方程计算磁感应强度b的频域解。
[0042]
5)、判断所有频点是否计算完成，若未完成，重复步骤3)，若完成则输出响应矩阵。
[0043]
6)、选取时间范围为10
‑5‑
10
‑2s，将响应矩阵带入正余弦变换滤波算法中，逐点进行频时变换。
[0044]
7)、判断是否完成全部频点的频时变换，若未完成，重复步骤6)，若完成则输出时域响应结果。
[0045]
8)、根据所选时刻的磁场值，对响应的感应电动势进行结果显示。
[0046]
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。