CN113962077A - 三维各向异性强磁场数值模拟方法、装置、设备及介质 - Google Patents

Abstract

三维各向异性强磁场数值模拟方法、装置、设备及介质，对包含目标区域的初始棱柱体模型进行剖分，对剖分得到的小棱柱体的磁化率张量赋值，得到异常棱柱体模型；设定高斯参数，计算离散偏移波数，计算空间域背景场磁场强度，获得磁化强度的计算模型；将空间域异常场磁位和磁化强度满足的三维拉普拉斯方程转为空间波数域一维常微分方程；结合设定的空间波数域异常场需满足的边界条件，将空间波数域异常场满足的边值问题模型转化为等价的变分问题模型并求解，得到空间波数域异常场磁位、空间波数域异常场磁场强度，并通过反傅里叶变换得到空间域异常场磁位以及空间域异常场磁场强度。本发明能够更加准确的对强磁性介质进行磁场数值模拟。

Description

三维各向异性强磁场数值模拟方法、装置、设备及介质

技术领域

本发明属于强磁体数值模拟技术领域，特别涉及一种三维各向异性强磁场数值模拟方法、装置、设备及介质。

背景技术

强磁体的自退磁效应多年来一直是研究的热点。磁介质被磁化后在其表面或内部不均匀处将产生磁荷，这种面磁荷或体磁荷在磁介质内所产生的磁场称为自退磁场。在磁测资料的处理和解释中,在磁性体磁化率小于0.1SI的情况下可忽略退磁效应，此时有效感应磁化强度是磁化率与外部地磁场的乘积,而磁化方向平行于地磁场方向。当磁性体的磁化率逐渐增大时,退磁效应越来越强，在磁性体内部退磁场的作用下,会在磁性体内部产生与外部磁场方向相反的自退磁场，有效感应磁化强度的幅值会减小,而磁化方向也会偏离原来的外部地磁场方向，导致磁异常的幅值减小,形态产生畸变,从而损害磁测数据的处理和解释工作。

岩石多具有磁化率各向异性，尤其是金属磁体，其退磁效应就更加复杂。因此，从各向异性角度去研究地下介质或磁性材料磁场的分布规律，对于寻找石油、天然气、有用矿产资源，研究深部、区域、全球地质构造、地震预报和监测，解决环境、工程、灾害等一系列地质问题具有十分重要的意义。

磁性体为棱柱体、球体、柱体或椭球体等规则形状时，通过正演计算可以得到解析解。文献(Bhattacharyya B K.Magnetic anomalies due to prism-shaped bodies witharbitrarypolarization.Geophysics，1964.29(4):517-531.)文献中提到，绝大部分磁性体为非规则形体，其正演计算时将非规则磁性体剖分成多个规则形体，之后对单个规则形体求解析解，并对多个规则体进行累加，便可得到任意非规则磁性体的磁场分布，但是该方法忽略自退磁场，仅适用于低磁化率小于0.1SI的情况，当磁化率大于0.1SI时，其正演计算复杂且困难。文献(方华竹.任意形状强磁性三度体磁异常的计算.地质学报.1978.(01):65-80.)用加约束的赛德尔迭代法求解关于有效磁荷面密度σ的Fredholm积分方程。文献(Eskola L,Tervo T.1980.Solving the magnetostatic field problem(a case of highsusceptibility)by means of the method of subsections[J].Geoexploration.18(2):79～95.)同样考虑了退磁效应，采用面积分方法进行计算，用物体表面的各种磁场分量来表示物体产生的磁场。文献(Kostrov N P.Calculation of magnetic anomalies causedby 2D bodies of arbitrary shape with consideration of demagnetization[J].Geophysical Prospecting.2007.55(1):p.91-115.)基于体积分(VIE)方法提出采用三角单元的体积分法求解2D问题，能够精确计算相对磁导率在2～20范围内或者离磁性体一个单元大小距离的磁场，但是相对磁导率值更高时，算法的稳定性降低。文献(Ouyang F,ChenL.2019.Iterative magnetic forward modeling for high susceptibility based onintegral equation and Gauss-FFT[J].Geophysics.85(1):1-47.)运用积分方程和高斯傅里叶变换法，对磁性体进行棱柱剖分，实现了强磁体的正演计算，但是在处理复杂强磁性体时需要进行精细剖分，计算效率低，且未考虑磁化率各向异性特征。

综上，目前强磁场数值模拟的研究绝大部分没有考虑自退磁场，即只考虑了弱磁的情况，与强磁场区别较大。部分考虑了自退磁场的强磁场数值模拟方法，仍存在效率低，精度不够，或者只能计算规则异常体的不足。且目前强磁体的磁场数值模拟几乎未曾研究过磁化率各向异性特征，而强磁性矿物多具有磁化率各向异性特征，因此急需提出一种对任意形状各向异性强磁性介质磁场的高效高精度模拟方法。

发明内容

针对目前强磁场数值模拟较少考虑自退磁场和磁化率各向异性特征，以及计算效率较低的问题，本发明旨在提出一种三维各向异性强磁场数值模拟方法、装置、设备及介质。

为实现上述技术目的，本发明提出的技术方案为：

一方面，本发明提供一种三维各向异性强磁场数值模拟方法，包括：

获取包含异常体在内的三维目标区域，建立包含目标区域的初始棱柱体模型；

对包含目标区域的初始棱柱体模型沿x、y、z方向分别进行等间隔剖分，得到多个小棱柱体，对所有小棱柱体的磁化率张量进行赋值，得到异常棱柱体模型；

根据异常棱柱体模型以及给定的x方向和y方向上的高斯参数，计算x和y方向离散偏移波数；

根据地球主磁场模型，计算空间域背景场磁场强度；

根据空间域背景场磁场强度、空间域异常场磁场强度，得到磁化强度的计算模型；

利用二维傅里叶变换将空间域异常场磁位和磁化强度满足的三维拉普拉斯方程转为空间波数域一维常微分方程；

基于空间波数域一维常微分方程，并结合设定的空间波数域异常场磁位需满足的边界条件，将空间波数域异常场磁位满足的边值问题模型转化为等价的变分问题模型；

通过求解变分问题模型，得到空间波数域异常场磁位；

基于空间波数域异常场磁位，求得空间波数域异常场磁场强度；

通过反傅里叶变换将空间波数域异常场磁位以及空间波数域异常场磁场强度转换为空间域异常场磁位以及空间域异常场磁场强度；

判断当前是否满足迭代终止条件，如满足则输出当前计算得到的空间域异常场磁位以及空间域异常场磁场强度。

另一方面，本发明提供一种三维各向异性强磁场数值模拟装置，包括：

第一模块，用于获取包含异常体在内的三维目标区域，建立包含目标区域的初始棱柱体模型；

第二模块，对包含目标区域的初始棱柱体模型沿x、y、z方向分别进行等间隔剖分，得到多个小棱柱体，对所有小棱柱体的磁化率张量进行赋值，得到异常棱柱体模型；

第三模块，用于根据异常棱柱体模型以及给定的x方向和y方向上的高斯参数，计算x和y方向离散偏移波数；

第四模块，用于根据地球主磁场模型，计算空间域背景场磁场强度；

第五模块，用于根据空间域背景场磁场强度、空间域异常场磁场强度，得到磁化强度的计算模型；

第六模块，用于利用二维傅里叶变换将空间域异常场磁位和磁化强度满足的三维拉普拉斯方程转为空间波数域一维常微分方程；

第七模块，用于基于空间波数域一维常微分方程，并结合设定的空间波数域异常场磁位需满足的边界条件，将空间波数域异常场磁位满足的边值问题模型转化为等价的变分问题模型；

第八模块，通过求解变分问题模型，得到空间波数域异常场磁位；

第九模块，用于基于空间波数域异常场磁位，求得空间波数域异常场磁场强度；

第十模块，通过反傅里叶变换将空间波数域异常场磁位以及空间波数域异常场磁场强度转换为空间域异常场磁位以及空间域异常场磁场强度；

第十一模块，用于判断当前是否满足迭代终止条件，如满足则输出当前计算得到的空间域异常场磁位以及空间域异常场磁场强度。

另一方面，本发明提供一种计算机设备，包括存储器和处理器，存储器存储有计算机程序，处理器执行计算机程序时实现以下步骤：

获取包含异常体在内的三维目标区域，建立包含目标区域的初始棱柱体模型；

对包含目标区域的初始棱柱体模型沿x、y、z方向分别进行等间隔剖分，得到多个小棱柱体，对所有小棱柱体的磁化率张量进行赋值，得到异常棱柱体模型；

根据异常棱柱体模型以及给定的x方向和y方向上的高斯参数，计算x和y方向离散偏移波数；

根据地球主磁场模型，计算空间域背景场磁场强度；

根据空间域背景场磁场强度、空间域异常场磁场强度，得到磁化强度的计算模型；

利用二维傅里叶变换将空间域异常场磁位和磁化强度满足的三维拉普拉斯方程转为空间波数域一维常微分方程；

基于空间波数域一维常微分方程，并结合设定的空间波数域异常场磁位需满足的边界条件，将空间波数域异常场磁位满足的边值问题模型转化为等价的变分问题模型；

通过求解变分问题模型，得到空间波数域异常场磁位；

基于空间波数域异常场磁位，求得空间波数域异常场磁场强度；

通过反傅里叶变换将空间波数域异常场磁位以及空间波数域异常场磁场强度转换为空间域异常场磁位以及空间域异常场磁场强度；

判断当前是否满足迭代终止条件，如满足则输出当前计算得到的空间域异常场磁位以及空间域异常场磁场强度。

再一方面，本发明还提供一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现以下步骤：

获取包含异常体在内的三维目标区域，建立包含目标区域的初始棱柱体模型；

对包含目标区域的初始棱柱体模型沿x、y、z方向分别进行等间隔剖分，得到多个小棱柱体，对所有小棱柱体的磁化率张量进行赋值，得到异常棱柱体模型；

根据异常棱柱体模型以及给定的x方向和y方向上的高斯参数，计算x和y方向离散偏移波数；

根据地球主磁场模型，计算空间域背景场磁场强度；

根据空间域背景场磁场强度、空间域异常场磁场强度，得到磁化强度的计算模型；

利用二维傅里叶变换将空间域异常场磁位和磁化强度满足的三维拉普拉斯方程转为空间波数域一维常微分方程；

基于空间波数域一维常微分方程，并结合设定的空间波数域异常场磁位需满足的边界条件，将空间波数域异常场磁位满足的边值问题模型转化为等价的变分问题模型；

通过求解变分问题模型，得到空间波数域异常场磁位；

基于空间波数域异常场磁位，求得空间波数域异常场磁场强度；

通过反傅里叶变换将空间波数域异常场磁位以及空间波数域异常场磁场强度转换为空间域异常场磁位以及空间域异常场磁场强度；

判断当前是否满足迭代终止条件，如满足则输出当前计算得到的空间域异常场磁位以及空间域异常场磁场强度。

与现有技术相比，本发明的优点在于：

1、本发明考虑了很多异常体具有的磁化率各向异性特征，通过对异常棱柱体模型中的所有小棱柱体的磁化率张量进行赋值。具体地小棱柱体的磁化率张量是一个对称张量，有六个独立分量，可以用三个主磁化率χ₁、χ₂和χ₃和三个欧拉角α、β、γ来表示。这样使得本发明提出的方案能够与实际地质情况更相符。

2、本发明考虑到自退磁效应，能够更加准确的对强磁性介质进行磁场数值模拟。

3、进一步地，本发明通过傅里叶变换将三维问题降成一维，对于一维常微分方程采用有限单元法进行求解，每个单元内部采用形函数二次插值，提高了计算精度和计算效率，且并行性好。

附图说明

图1是本发明一实施例中流程图；

图2是本发明一实施例中对包含目标区域的初始棱柱体模型进行剖分的示意图；

图3是坐标旋转示意图，其中(a)表示绕x轴旋转α角度，(b)表示在(a)的基础上绕y轴旋转β角度，(c)表示在(b)的基础上绕z轴旋转γ角度；

图4是本发明一实施例中的目标区域以及异常体的示意图；

图5是本发明一实施例中的模型数值解、解析解及其绝对误差的示意图；其中(a)、(b)、(c)分别代表B_ax的数值解，B_ax的解析解以及B_ax的数值解和解析解的绝对误差；(d)、(e)、(f)分别代表B_ay的数值解，B_ay的解析解以及B_ay的数值解和解析解的绝对误差；(g)、(h)、(i)分别代表B_az的数值解，B_az的解析解以及B_az的数值解和解析解的绝对误差；

图6是本发明一实施例中计算机设备的内部结构图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚明白，下面将以附图及详细叙述清楚说明本发明所揭示内容的精神，任何所属技术领域技术人员在了解本发明内容的实施例后，当可由本发明内容所教示的技术，加以改变及修饰，其并不脱离本发明内容的精神与范围。本发明的示意性实施例及其说明用于解释本发明，但并不作为对本发明的限定。

参照图1，本发明一实施例中，提供一种三维各向异性强磁场数值模拟方法，包括：

(S1)获取包含异常体在内的三维目标区域，建立包含目标区域的初始棱柱体模型；

(S2)对包含目标区域的初始棱柱体模型沿x、y、z方向分别进行等间隔剖分，得到多个小棱柱体，如图2所示。然后对所有小棱柱体的磁化率张量进行赋值，得到异常棱柱体模型；

(S3)根据异常棱柱体模型以及给定的x方向和y方向上的高斯参数，计算x和y方向离散偏移波数；

(S4)根据地球主磁场模型，计算空间域背景场磁场强度；

(S5)根据空间域背景场磁场强度、空间域异常场磁场强度，得到磁化强度的计算模型；

(S6)利用二维傅里叶变换将空间域异常场磁位和磁化强度满足的三维拉普拉斯方程转为空间波数域一维常微分方程；

(S7)基于空间波数域一维常微分方程，并结合设定的空间波数域异常场磁位需满足的边界条件，将空间波数域异常场磁位满足的边值问题模型转化为等价的变分问题模型；

(S8)通过有限单元法求解变分问题模型，得到空间波数域异常场磁位；

(S9)基于空间波数域异常场磁位，求得空间波数域异常场磁场强度；

(S10)通过反傅里叶变换将空间波数域异常场磁位以及空间波数域异常场磁场强度转换为空间域异常场磁位以及空间域异常场磁场强度；

(S11)判断当前是否满足迭代终止条件，如不满足则返回(S8)，如满足则输出当前计算得到的空间域异常场磁位以及空间域异常场磁场强度。

本发明中步骤(S1)中目标区域中的异常体的形状、大小、磁化率分布均不限，可以是任意形状、任意大小、任意磁化率分布的强磁性介质。本发明不光适用于具有磁化率各向同性特征的强磁性介质，也尤其适合具有磁化率各向异性特征的强磁性介质。

考虑异常体多具有磁化率各向异性特征，在本发明的步骤(S2)中，对所有小棱柱体的磁化率进行赋值，磁化率为磁化率张量，单位为SI，如下所示：

磁化率

是一个对称张量，有六个独立分量，可以用三个主磁化率χ₁、χ₂和χ₃和三个欧拉角α、β、γ来表示，其坐标旋转示意图如图3所示。图3中，(a)表示绕x轴旋转α角度，(b)表示在(a)的基础上绕y轴旋转β角度，(c)表示在(b)的基础上绕z轴旋转γ角度。坐标系中点的坐标变化(x,y,z)→(x₁,y₁,z₁)→(x₂,y₂,z₂)→(x_r,y_r,z_r)。

其中三个旋转矩阵分别为:

三个旋转矩阵相乘能够完成坐标系统的转换，能够表征磁化率各向异性特征的磁化率张量

表示如下

其中，

而矩阵D是与三个欧拉角有关的变换矩阵，具体表达式为

每个小棱柱体的磁化率张量都满足式(5)，根据磁化率分布数据对每个小棱柱体顶点的磁化率张量进行赋值，取值可以任意，没有磁化率异常的部分磁化率张量的九个量都为0，以此建立包含任意形状、任意磁化率分布异常体的异常棱柱体模型。

可以理解，在本发明的步骤(S3)中，可以参照本领域中的已有方法中高斯参数的设定方法进行高斯参数的设定，并进行x和y方向离散偏移波数的计算。

在本发明一实施例的步骤(S3)中，给定x方向的高斯点个数N_x，区间[-1,1]上高斯点t_a、高斯系数A_a，其中，a＝1,2,...,N_x；给定y方向的高斯点个数N_y，区间[-1,1]上高斯点t_b、高斯系数A_b，其中，b＝1,2,...,N_y。

在本发明一实施例的步骤(S3)中，x和y方向离散偏移波数，通过以下方法计算：

式中，

其中：k_x表示x方向的偏移波数，Δk_x表示x方向基波数，NN_x表示异常棱柱体模型x方向上小棱柱体的剖分个数，Δx表示异常棱柱体模型x方向上小棱柱体的单元长度；k_y表示y方向的偏移波数，Δk_y表示y方向基波数，NN_y表示异常棱柱体模型y方向上小棱柱体的剖分个数，Δy表示异常棱柱体模型y方向上小棱柱体的单元长度。

在本发明一实施例的步骤(S4)中，包括：

根据地球主磁场模型IGRF，计算异常棱柱体模型中每个棱柱体顶点处的地球主磁场强度，将其作为空间域背景场磁场强度H₀，为数值模拟中的背景场，即无异常时的磁场，单位为A/m。空间域背景场磁场强度H₀其三个方向的分量分别表示为H_0x、H_0y、H_0z：

H_0x＝||H₀||·cos(af)·cos(bta) (10)

H_0y＝||H₀||·cos(af)·sin(bta) (11)

H_0z＝||H₀||·sin(af) (12)

其中：||H₀||表示背景场H₀的L2范数，af为目标区域磁倾角，bta为目标区域磁偏角。

在本发明一实施例的步骤(S5)中，磁化强度的计算模型为：

其中H_a表示异常棱柱体模型中每个棱柱体顶点处由异常体产生的磁场强度，即其空间域异常场磁场强度，为数值模拟中的异常场，即由磁化率异常所产生的磁场，单位为A/m，其三个分量分别为H_ax、H_ay、H_az。空间域总磁场强度H为背景场与异常场之和。

在本发明一实施例的步骤(S6)中，空间域异常场磁位U_a和磁化强度M满足的三维拉普拉斯方程：

▽²U_a＝▽·M (14)

式中，▽²为拉普拉斯算子，展开为

▽·M表示对磁化强度M求散度。式(14)展开即为：

对式(14)进行二维傅里叶变换，得到空间波数域一维常微分方程：

其中

表示空间波数域异常场磁位，

为波数域磁化强度，k_x、k_y分别为x、y方向的偏移波数。公式(15)即空间波数域异常场磁位

所满足的一维常微分方程。

步骤(S7)中，为了得到控制方程的定解，需给出合适的边界条件。在笛卡尔坐标系下，取Z轴垂直向下为正向，取水平地面为上边界Z_min，取地下离异常体足够远处为下边界Z_max，其上、下边界条件满足：

上边界：

下边界：

其中，

联立公式(15)、(16)和(17)，得到空间波数域异常场磁位满足的边值问题模型：

运用变分法，将空间波数域异常场磁位满足的边值问题模型转化为等价的变分问题模型：

在本发明中，通过有限单元法求解变分问题模型，得到空间波数域异常场磁位

空间波数域异常场磁位

所满足的常微分方程的右端项中为包含有背景场和异常场，而异常场未知，因此采用迭代求解。空间域总磁场强度H为背景场与异常场之和。本发明在初始时刻，结合方程(13)和方程(18)，可知磁化强度M由背景场H₀和异常场H_a之和与磁化率的乘积所决定，且H＝H₀+H_a,而异常场H_a未知，因此初次迭代假设H_a为0，将初始的空间域总磁场强度H用背景场H₀替代，从而将一维偏微分方程变为一维常微分方程求解，得到第一个异常场，之后再将得到的异常场和背景场之和作为新的空间域总磁场强度进行下一次求解，以此迭代获得空间波数域异常场磁位

对于空间波数域异常场磁位满足的变分问题模型，采用基于二次插值的一维有限单元法进行求解，能够充分考虑计算精度和计算效率，能够实现对复杂模型的高精度模拟，且可以利用追赶法实现对角线性方程组的快速求解。

本发明一实施例的步骤(S9)中，通过下式求得空间波数域异常场磁场强度：

其中，i为虚数，d/dz表示对z求微分，

表示空间波数域异常场磁位对z求微分。

最后由异常场磁感应强度B_a(单位为T)与异常场磁场强度H_a的关系(21),可求得磁感应强度B_a，进而得到B_a的三个分量B_ax，B_ay，B_az。

B_a＝μH_a (21)

其中，μ为介质的绝对磁导率，单位为H/m。

绝对磁导率μ与

之间的关系如方程(22)所示，μ₀为真空中磁导率，μ₀＝4π×10^{- 7}H/m。

本发明的步骤(S10)中，通过反傅里叶变换将空间波数域异常场磁位

以及空间波数域异常场磁场强度

转换为空间域异常场磁位U_a以及空间域异常场磁场强度H_a。

可以理解，预设的迭代终止条件是指预先设置的模型计算约束条件，用于约束整个模型进行性能计算的过程趋向收敛，以使模型能够输出满足条件的结果。本发明中可以将(S11)中的迭代终止条件设置为：每次迭代后的空间域总磁场强度H为空间域背景场磁场强度H₀和当次计算得到的空间域异常场磁场强度H_a之和。用j表示上一次迭代，j+1表示本次迭代，设置以下迭代收敛条件：

|H^j+1-H^j|/H^j+1＜10^-4 (21)

当满足以上迭代收敛条件时，迭代停止，其中H^j上一次迭代计算得到的空间域总磁场强度，H^j+1为本次计算得到的磁场总强度。

当然，实际应用中，本领域技术人员也可基于现有技术、本领域的惯用技术手段或者公知常识，设定其他的迭代终止条件，不局限于本申请上述优选实施例中所述设置的迭代终止条件。

下面对本发明提供的三维各向异性强磁场数值模拟方法的精度和效率进行检验。

测试电脑配置为i5-4590，主频3.30GHz，内存12GB。

目标区域为三维棱柱体结构，其大小为1000m×1000m×1000m。目标区域中的异常体为各向异性三轴椭球体结构，其长轴、中轴、短轴分别为180m、170m、160m，异常体中心位于目标区域的中心(500m，500m，500m)处，其示意图如图4所示。

异常体其磁化率张量主轴为(2SI，2SI，2SI)，欧拉角(60°，60°，30°)。目标区域背景磁场强度为50000nT，磁倾角45°，磁偏角9°。x，y，z方向节点数取101,101,101。在单线程的情况下，采取Gauss-FFT，达到收敛条件需迭代7次，用时70.49s，占用内存1.98GB，占用内存较低且效率较高。其B_ax，B_ay，B_az与解析解的绝对误差均与场值相差两个数量级以上，满足数值计算要求。图5为模型计算B_ax，B_ay，B_az分量的数值解、解析解及其绝对误差，图5中(a)、(b)和(c)分别代表B_ax的数值解，B_ax的解析解以及B_ax的数值解和解析解的绝对误差；图5中(d)、(e)和(f)分别代表B_ay的数值解，B_ay的解析解以及B_ay的数值解和解析解的绝对误差；图5中(g)、(h)和(i)分别代表B_az的数值解，B_az的解析解以及B_az的数值解和解析解的绝对误差。

本发明在考虑自退磁效应的基础上，将磁化率由各向同性发展到各向异性，使得磁异常刻画与实际地质体更加接近，且可以计算VTI、HTI、TTI等典型各向异性介质模型的磁场分布；运用傅里叶变换将三维问题降至一维，只保留z方向，运用一维有限单元法，单元内部采用形函数二次插值，从而对微分方程进行迭代求解，极大的提高了计算精度和计算效率，算法并行度好，占用内存少。因此本发明提供了一种可以计算任意形状、任意磁化率分布的各向异性强磁性介质的磁场数值模拟方法，为磁测资料解释和磁异常反演奠定了基础。

本发明一实施例中提供发明提供一种三维各向异性强磁场数值模拟装置，包括：

第一模块，用于获取包含异常体在内的三维目标区域，建立包含目标区域的初始棱柱体模型；

第二模块，对包含目标区域的初始棱柱体模型沿x、y、z方向分别进行等间隔剖分，得到多个小棱柱体，对所有小棱柱体的磁化率张量进行赋值，得到异常棱柱体模型；

第三模块，用于根据异常棱柱体模型以及给定的x方向和y方向上的高斯参数，计算x和y方向离散偏移波数；

第四模块，用于根据地球主磁场模型，计算空间域背景场磁场强度；

第五模块，用于根据空间域背景场磁场强度、空间域异常场磁场强度，得到磁化强度的计算模型；

第六模块，用于利用二维傅里叶变换将空间域异常场磁位和磁化强度满足的三维拉普拉斯方程转为空间波数域一维常微分方程；

第七模块，用于基于空间波数域一维常微分方程，并结合设定的空间波数域异常场磁位需满足的边界条件，将空间波数域异常场磁位异常场满足的边值问题模型转化为等价的变分问题模型；

第八模块，通过求解变分问题模型，得到空间波数域异常场磁位；

第九模块，用于基于空间波数域异常场磁位，求得空间波数域异常场磁场强度；

第十模块，通过反傅里叶变换将空间波数域异常场磁位以及空间波数域异常场磁场强度转换为空间域异常场磁位以及空间域异常场磁场强度；

第十一模块，用于判断当前是否满足迭代终止条件，如满足则输出当前计算得到的空间域异常场磁位以及空间域异常场磁场强度。

上述各模块功能的实现方法，可以采用前述各实施例中相同的方法实现，在此不再赘述。

在本一个实施例中，提供了一种计算机设备，该计算机设备可以是服务器，其内部结构图可以如图6所示。该计算机设备包括通过系统总线连接的处理器、存储器、网络接口和数据库。其中，该计算机设备的处理器用于提供计算和控制能力。该计算机设备的存储器包括非易失性存储介质、内存储器。该非易失性存储介质存储有操作系统、计算机程序和数据库。该内存储器为非易失性存储介质中的操作系统和计算机程序的运行提供环境。该计算机设备的数据库用于存储样本数据。该计算机设备的网络接口用于与外部的终端通过网络连接通信。该计算机程序被处理器执行时以实现上述三维各向异性强磁场数值模拟方法。

本领域技术人员可以理解，图6中示出的结构，仅仅是与本申请方案相关的部分结构的框图，并不构成对本申请方案所应用于其上的计算机设备的限定，具体的计算机设备可以包括比图中所示更多或更少的部件，或者组合某些部件，或者具有不同的部件布置。

在一个实施例中，提供了一种计算机设备，包括存储器和处理器，该存储器存储有计算机程序，该处理器执行计算机程序时实现上述实施例中三维各向异性强磁场数值模拟方法的步骤。

在一个实施例中，提供了一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现上述实施例中三维各向异性强磁场数值模拟方法的步骤。

本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的计算机程序可存储于一非易失性计算机可读取存储介质中，该计算机程序在执行时，可包括如上述各方法的实施例的流程。其中，本申请所提供的各实施例中所使用的对存储器、存储、数据库或其它介质的任何引用，均可包括非易失性和/或易失性存储器。非易失性存储器可包括只读存储器(ROM)、可编程ROM(PROM)、电可编程ROM(EPROM)、电可擦除可编程ROM(EEPROM)或闪存。易失性存储器可包括随机存取存储器(RAM)或者外部高速缓冲存储器。作为说明而非局限，RAM以多种形式可得，诸如静态RAM(SRAM)、动态RAM(DRAM)、同步DRAM(SDRAM)、双数据率SDRAM(DDRSDRAM)、增强型SDRAM(ESDRAM)、同步链路(Synchlink)DRAM(SLDRAM)、存储器总线(Rambus)直接RAM(RDRAM)、直接存储器总线动态RAM(DRDRAM)、以及存储器总线动态RAM(RDRAM)等。

以上实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。

以上所述实施例仅表达了本申请的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本申请构思的前提下，还可以做出若干变形和改进，这些都属于本申请的保护范围。因此，本申请专利的保护范围应以所附权利要求为准。

Claims (14)

1.三维各向异性强磁场数值模拟方法，其特征在于，包括：

获取包含异常体在内的三维目标区域，建立包含目标区域的初始棱柱体模型；

对包含目标区域的初始棱柱体模型沿x、y、z方向分别进行等间隔剖分，得到多个小棱柱体，对所有小棱柱体的磁化率张量进行赋值，得到异常棱柱体模型；

根据异常棱柱体模型以及给定的x方向和y方向上的高斯参数，计算x和y方向离散偏移波数；

根据地球主磁场模型，计算空间域背景场磁场强度；

根据空间域背景场磁场强度、空间域异常场磁场强度，得到磁化强度的计算模型；

利用二维傅里叶变换将空间域异常场磁位和磁化强度满足的三维拉普拉斯方程转为空间波数域一维常微分方程；

基于空间波数域一维常微分方程，并结合设定的空间波数域异常场磁位需满足的边界条件，将空间波数域异常场磁位满足的边值问题模型转化为等价的变分问题模型；

通过求解变分问题模型，得到空间波数域异常场磁位；

基于空间波数域异常场磁位，求得空间波数域异常场磁场强度；

通过反傅里叶变换将空间波数域异常场磁位以及空间波数域异常场磁场强度转换为空间域异常场磁位以及空间域异常场磁场强度；

判断当前是否满足迭代终止条件，如满足则输出当前计算得到的空间域异常场磁位以及空间域异常场磁场强度。

2.根据权利要求1所述的三维各向异性强磁场数值模拟方法，其特征在于，根据磁化率分布数据对每个小棱柱体顶点的磁化率张量进行赋值，磁化率张量单位为SI，如下所示：

是一个对称张量，有六个独立分量，用三个主磁化率χ₁、χ₂和χ₃和三个欧拉角α、β、γ来表示，其中，

磁化率张量

是一个对称张量，有六个独立分量，用三个主磁化率χ₁、χ₂和χ₃和三个欧拉角α、β、γ来表示，而矩阵D是与三个欧拉角有关的变换矩阵，具体表达式为：

3.根据权利要求2所述的三维各向异性强磁场数值模拟方法，其特征在于，给定高斯参数包括：给定x方向的高斯点个数N_x，区间[-1,1]上高斯点t_a、高斯系数A_a，其中，a＝1,2,...,N_x；给定y方向的高斯点个数N_y，区间[-1,1]上高斯点t_b、高斯系数A_b，其中，b＝1,2,...,N_y。

4.根据权利要求3所述的三维各向异性强磁场数值模拟方法，其特征在于，x和y方向离散偏移波数，通过以下方法计算：

式中，

其中：k_x表示x方向的偏移波数，Δk_x表示x方向基波数，NN_x表示异常棱柱体模型x方向上小棱柱体的剖分个数，Δx表示异常棱柱体模型x方向上小棱柱体的单元长度；k_y表示y方向的偏移波数，Δk_y表示y方向基波数，NN_y表示异常棱柱体模型y方向上小棱柱体的剖分个数，Δy表示异常棱柱体模型y方向上小棱柱体的单元长度。

5.根据权利要求1至4中任一项所述的三维各向异性强磁场数值模拟方法，其特征在于，根据地球主磁场模型IGRF，计算异常棱柱体模型中每个棱柱体顶点处的地球主磁场强度，将其作为空间域背景场磁场强度H₀。

6.根据权利要求5所述的三维各向异性强磁场数值模拟方法，其特征在于，磁化强度的计算模型为：

其中H_a表示异常棱柱体模型中每个棱柱体顶点处由异常体产生的磁场强度，即其空间域异常场磁场强度，H表示空间域总磁场强度。

7.根据权利要求6所述的三维各向异性强磁场数值模拟方法，其特征在于，对空间域异常场磁位U_a和磁化强度M满足的三维拉普拉斯方程

进行二维傅里叶变换，得到空间波数域一维常微分方程：

其中

为拉普拉斯算子，

表示对磁化强度M求散度；

表示空间波数域异常场磁位，

为波数域磁化强度，k_x、k_y分别为x、y方向的偏移波数。

8.根据权利要求7所述的三维各向异性强磁场数值模拟方法，其特征在于，空间波数域异常场磁位需满足的边界条件为：

上边界：

下边界：

其中，

得到空间波数域异常场磁位满足的边值问题模型，如下：

运用变分法，将空间波数域异常场磁位满足的边值问题模型转化为等价的变分问题模型。

9.根据权利要求7所述的三维各向异性强磁场数值模拟方法，其特征在于，通过有限单元法求解变分问题模型，得到空间波数域异常场磁位

10.根据权利要求7所述的三维各向异性强磁场数值模拟方法，其特征在于，通过下式求得空间波数域异常场磁场强度：

其中，i为虚数，d/dz表示对z求微分。

11.根据权利要求1、2、3、4、6、7、8、9或10所述的三维各向异性强磁场数值模拟方法，其特征在于，迭代终止条件设置为：

|H^j+1-H^j|/H^j+1＜10^-4

其中H^j表示第j次迭代计算得到的空间域总磁场强度，H^j+1表示第j+1次计算得到的空间域总磁场强度，空间域总磁场强度为空间域背景场磁场强度和当次计算得到的空间域异常场磁场强度之和。

12.一种三维各向异性强磁场数值模拟装置，其特征在于，包括：

第一模块，用于获取包含异常体在内的三维目标区域，建立包含目标区域的初始棱柱体模型；

第二模块，对包含目标区域的初始棱柱体模型沿x、y、z方向分别进行等间隔剖分，得到多个小棱柱体，对所有小棱柱体的磁化率张量进行赋值，得到异常棱柱体模型；

第三模块，用于根据异常棱柱体模型以及给定的x方向和y方向上的高斯参数，计算x和y方向离散偏移波数；

第四模块，用于根据地球主磁场模型，计算空间域背景场磁场强度；

第五模块，用于根据空间域背景场磁场强度、空间域异常场磁场强度，得到磁化强度的计算模型；

第六模块，用于利用二维傅里叶变换将空间域异常场磁位和磁化强度满足的三维拉普拉斯方程转为空间波数域一维常微分方程；

第七模块，用于基于空间波数域一维常微分方程，并结合设定的空间波数域异常场磁位需满足的边界条件，将空间波数域异常场磁位满足的边值问题模型转化为等价的变分问题模型；

第八模块，通过求解变分问题模型，得到空间波数域异常场磁位；

第九模块，用于基于空间波数域异常场磁位，求得空间波数域异常场磁场强度；

第十模块，通过反傅里叶变换将空间波数域异常场磁位以及空间波数域异常场磁场强度转换为空间域异常场磁位以及空间域异常场磁场强度；

第十一模块，用于判断当前是否满足迭代终止条件，如满足则输出当前计算得到的空间域异常场磁位以及空间域异常场磁场强度。

13.一种计算机设备，包括存储器和处理器，存储器存储有计算机程序，其特征在于，处理器执行计算机程序时实现权利要求1、2、3、4、6、7、8、9或10所述三维各向异性强磁场数值模拟方法的步骤。

14.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，计算机程序被处理器执行时实现权利要求1、2、3、4、6、7、8、9或10所述三维各向异性强磁场数值模拟方法的步骤。

Images