二维强磁场数值模拟方法、装置、设备及介质与流程

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


背景技术：

2.磁法勘探是通过观测和分析由岩石、矿石或其他探测对象磁性差异引起的磁异常，进而研究地下地质构造的一种地球物理勘探方法，主要应用于矿产资源勘探中。而自然界中，铁矿大多都是高磁化率和强剩磁的，且退磁效应的存在会导致异常幅值变小，当不考虑退磁效应影响时所估算出来的资源量会比实际资源量要少很多。因此强磁场的勘探对于对磁测资料的正确处理、准确的地质解释有重要意义。
3.理论上，复杂形体重磁异常正演可以分为三度体和二度体，解释人员可以使用二维模型、2.5维模型或者三维模型来模拟不同复杂程度的场源。实际情况下，所有的地质体都是三维的，磁异常正演应该都以三维模型为基础，但是三维模型的计算量往往都非常大。在实际的重磁勘探过程中，有很多线性地质体，比如断层、接触带等，这类地质体走向方向的尺度远比垂直其走向方向的尺度大，那么它们的实际场源分布就可以用走向方向无限延伸的二度体代替，不仅可以大大减少计算时间，相应的反演算法也比较容易实现。
4.然而目前而强磁勘探技术领域中，对于强磁场的数值模拟大部分仍停留在空间域中，空间域数值模拟方法应用于大规模复杂介质条件下强磁场的数据处理和精细反演成像时效率低，因此，研究一种高效、高精度的强磁场数值模拟有重要意义。


技术实现要素：

5.本发明着眼于二维强磁勘探，旨在提出一种二维强磁场数值模拟方法、装置、设备及介质，其适用于二维强磁性介质的高效、高精度的数值模拟计算。正演是反演的基础，该发明能有效地提高二维强磁场反演的计算效率，对于扩大强磁场勘探应用范围和地质解释准确度都有重要意义。
6.为实现上述技术目的，本发明提出的技术方案为：
7.一方面，本发明提供一种二维强磁场数值模拟方法，包括：
8.对于沿y轴方向无限延伸的异常体，在xoz平面确定目标区域，所述异常体xoz方向的整个截面包含在所述目标区域中，对目标区域进行网格剖分，得到多个长方形单元，对各长方形单元顶点处的磁化率χ进行赋值，磁化率χ为常数；
9.根据磁化率张量、空间域背景场磁场强度、空间域异常场磁场强度，得到空间域磁化强度；
10.利用一维傅里叶变换将空间域异常场磁位和空间域磁化强度满足的二维偏微分方程转为空间波数混合域一维常微分方程；
11.基于空间波数混合域一维常微分方程，并结合设定的空间波数混合域异常场磁位需满足的边界条件，将空间波数混合域异常场磁位满足的边值问题模型转化为等价的变分
问题模型；
12.求解变分问题模型，得到空间波数混合域异常场磁位；
13.基于空间波数混合域异常场磁位，求得空间波数混合域异常场磁场强度；
14.通过反傅里叶变换将空间波数混合域异常场磁位以及空间波数混合域异常场磁场强度转换为空间域异常场磁位以及空间域异常场磁场强度；
15.判断当前是否满足迭代终止条件，如满足则输出当前计算得到的空间域异常场磁位以及空间域异常场磁场强度。
16.另一方面，本发明提供一种二维强磁场数值模拟装置，包括：
17.第一模块，用于对沿y轴方向无限延伸的异常体，在xoz平面确定目标区域，所述异常体xoz方向的整个截面包含在所述目标区域中，对目标区域进行网格剖分，得到多个长方形单元，对各长方形单元顶点处的磁化率χ进行赋值，磁化率χ为常数；
18.第二模块，用于根据磁化率张量、空间域背景场磁场强度、空间域异常场磁场强度，得到空间域磁化强度；
19.第三模块，利用一维傅里叶变换将空间域异常场磁位和空间域磁化强度满足的二维偏微分方程转为空间波数混合域一维常微分方程；
20.第四模块，基于空间波数混合域一维常微分方程，并结合设定的空间波数混合域异常场磁位需满足的边界条件，将空间波数混合域异常场磁位满足的边值问题模型转化为等价的变分问题模型；
21.第五模块，用于求解变分问题模型，得到空间波数混合域异常场磁位；
22.第六模块，基于空间波数混合域异常场磁位，求得空间波数混合域异常场磁场强度；
23.第七模块，通过反傅里叶变换将空间波数混合域异常场磁位以及空间波数混合域异常场磁场强度转换为空间域异常场磁位以及空间域异常场磁场强度；
24.第八模块，判断当前是否满足迭代终止条件，如满足则输出当前计算得到的空间域异常场磁位以及空间域异常场磁场强度。
25.另一方面，本发明提供一种计算机设备，包括存储器和处理器，存储器存储有计算机程序，处理器执行计算机程序时实现以下步骤：
26.对于沿y轴方向无限延伸的异常体，在xoz平面确定目标区域，所述异常体xoz方向的整个截面包含在所述目标区域中，对目标区域进行网格剖分，得到多个长方形单元，对各长方形单元顶点处的磁化率χ进行赋值，磁化率χ为常数；
27.根据磁化率张量、空间域背景场磁场强度、空间域异常场磁场强度，得到空间域磁化强度；
28.利用一维傅里叶变换将空间域异常场磁位和空间域磁化强度满足的二维偏微分方程转为空间波数混合域一维常微分方程；
29.基于空间波数混合域一维常微分方程，并结合设定的空间波数混合域异常场磁位需满足的边界条件，将空间波数混合域异常场磁位满足的边值问题模型转化为等价的变分问题模型；
30.求解变分问题模型，得到空间波数混合域异常场磁位；
31.基于空间波数混合域异常场磁位，求得空间波数混合域异常场磁场强度；
32.通过反傅里叶变换将空间波数混合域异常场磁位以及空间波数混合域异常场磁场强度转换为空间域异常场磁位以及空间域异常场磁场强度；
33.判断当前是否满足迭代终止条件，如满足则输出当前计算得到的空间域异常场磁位以及空间域异常场磁场强度。
34.再一方面，本发明还提供一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现以下步骤：
35.对于沿y轴方向无限延伸的异常体，在xoz平面确定目标区域，所述异常体xoz方向的整个截面包含在所述目标区域中，对目标区域进行网格剖分，得到多个长方形单元，对各长方形单元顶点处的磁化率χ进行赋值，磁化率χ为常数；
36.根据磁化率张量、空间域背景场磁场强度、空间域异常场磁场强度，得到空间域磁化强度；
37.利用一维傅里叶变换将空间域异常场磁位和空间域磁化强度满足的二维偏微分方程转为空间波数混合域一维常微分方程；
38.基于空间波数混合域一维常微分方程，并结合设定的空间波数混合域异常场磁位需满足的边界条件，将空间波数混合域异常场磁位满足的边值问题模型转化为等价的变分问题模型；
39.求解变分问题模型，得到空间波数混合域异常场磁位；
40.基于空间波数混合域异常场磁位，求得空间波数混合域异常场磁场强度；
41.通过反傅里叶变换将空间波数混合域异常场磁位以及空间波数混合域异常场磁场强度转换为空间域异常场磁位以及空间域异常场磁场强度；
42.判断当前是否满足迭代终止条件，如满足则输出当前计算得到的空间域异常场磁位以及空间域异常场磁场强度。
43.与现有技术相比，本发明的优点在于：
44.1、适用于二维强磁勘探数值模拟；
45.2、考虑到自退磁效应，能够更加准确地对强磁性介质进行数值模拟；
46.3、通过傅里叶变换将空间域磁位满足的二维偏微分方程转化为不同波数相互独立的一维常微分方程，采用有限单元法求解不同波数满足的常微分方程，提高了计算效率和计算精度。
附图说明
47.图1是本发明一实施例中流程图；
48.图2是本发明一实施例中的目标区域以及异常体的示意图；
49.图3是本发明一实施例中空间域异常场磁场强度的x分量、z分量的解析解和数值解以及相对误差结果图，(1)为空间域异常场磁场强度ha其x分量解析解和数值解结果图，(2)为空间域异常场磁场强度ha其x分量解析解和数值解的相对误差，(3)为空间域异常场磁场强度ha其z分量解析解和数值解结果图，(4)为空间域异常场磁场强度ha其z分量解析解和数值解的相对误差；
50.图4是本发明一实施例中计算机设备的内部结构图。
具体实施方式
51.为使本发明实施例的目的、技术方案和优点更加清楚明白，下面将以附图及详细叙述清楚说明本发明所揭示内容的精神，任何所属技术领域技术人员在了解本发明内容的实施例后，当可由本发明内容所教示的技术，加以改变及修饰，其并不脱离本发明内容的精神与范围。本发明的示意性实施例及其说明用于解释本发明，但并不作为对本发明的限定。
52.参照图1，本发明一实施例中，提供一种二维强磁场数值模拟方法，包括：
53.(s1)对于沿y轴方向无限延伸的异常体，在xoz平面确定目标区域，所述异常体xoz方向的整个截面包含在所述目标区域中，对目标区域进行网格剖分，得到多个长方形单元，对各长方形单元顶点处的磁化率χ进行赋值，磁化率χ为常数；
54.(s2)根据磁化率张量、空间域背景场磁场强度、空间域异常场磁场强度，得到空间域磁化强度；
55.(s3)利用一维傅里叶变换将空间域异常场磁位和空间域磁化强度满足的二维偏微分方程转为空间波数混合域一维常微分方程；
56.(s4)基于空间波数混合域一维常微分方程，并结合设定的空间波数混合域异常场磁位需满足的边界条件，将空间波数混合域异常场磁位满足的边值问题模型转化为等价的变分问题模型；
57.(s5)求解变分问题模型，得到空间波数混合域异常场磁位；
58.(s6)基于空间波数混合域异常场磁位，求得空间波数混合域异常场磁场强度；
59.(s7)通过反傅里叶变换将空间波数混合域异常场磁位以及空间波数混合域异常场磁场强度转换为空间域异常场磁位以及空间域异常场磁场强度；
60.(s8)判断当前是否满足迭代终止条件，如满足则输出当前计算得到的空间域异常场磁位以及空间域异常场磁场强度，如果不满足则返回(s5)。
61.本发明中步骤(s1)中沿y轴方向无限延伸的异常体为强磁性介质，其xoz方向的整个截面的形状、大小均不限，可以是圆柱形、长方形等其他规则或者不规则的形状。异常体磁化率分布不限。根据目标区域已有的地质解释资料和成果，各长方形单元顶点处的磁化率χ进行赋值，一个长方形单元顶点相当于一个节点，不同节点处的磁化率χ的取值可以不同。
62.本发明另一实施例中，步骤(s2)根据磁化率张量、空间域背景场磁场强度、空间域异常场磁场强度，得到空间域磁化强度，包括：
63.(s2.1)根据地球主磁场模型igrf，计算各长方形单元顶点处的地球主磁场强度，将其作为空间域背景场磁场强度h0；即为数值模拟中的背景场，即无异常时的磁场，单位为a/m。
64.(s2.2)空间域磁化强度m表示为：
[0065][0066]
其中ha表示空间域异常场磁场强度，h表示空间域总磁场强度，h＝h0 ha。
[0067]
本发明另一实施例中，步骤(s3)中，空间域异常场磁位和空间域磁化强度满足的二维偏微分方程为：
[0068]
[0069]
对所述二维偏微分方程进行一维傅里叶变换，得到空间波数混合域一维常微分方程：
[0070][0071]
其中为梯度算子，e
x
和ez分别为x和z方向上的单位向量，u表示空间域异常场磁位，波数k
x
为x方向上的偏移波数，表示空间波数混合域异常场磁位，分别为x和z方向上的空间波数混合域磁化强度，i为虚数单位。
[0072]
可以理解，在本发明的步骤(s3)中，可以参照本领域中的已有方法对x方向上的偏移波数k
x
进行设定。
[0073]
在本发明另一实施例中，x方向上的偏移波数k
x
通过以下步骤获得：
[0074]
给定x方向的高斯点个数n
x
，区间[-1,1]上高斯点ta以及高斯系数aa，其中，a＝1,2,...,n
x
；
[0075]
x方向上的偏移波数k
x
，如下：
[0076][0077]
其中，δk
x
表示x方向上的基波数，nn
x
表示x方向上小长方形单元的剖分个数，δx为x方向上小长方形单元的单元长度。
[0078]
本发明的步骤(s4)中，空间波数混合域异常场磁位需满足的边界条件，包括上边界条件和下边界条件，在笛卡尔坐标系下，取z轴垂直向下为正向，取水平地面为上边界z
min
，取地下离异常体足够远处为下边界z
max
。为了得到控制方程的定解，需要给出合适的边界条件，本领域技术人员可以根据实际清楚给出合适的边界条件。
[0079]
在本发明另一实施例中，空间波数混合域异常场磁位需满足的边界条件，如下：
[0080]
上边界条件为：
[0081][0082]
下边界条件为：
[0083][0084]
联立可得空间波数混合域异常场磁位满足的边值问题模型，如下：
[0085][0086]
运用变分原理，将空间波数混合域异常场磁位满足的边值问题模型转化为等价的变分问题模型：
[0087][0088]
空间域异常场磁场强度，h表示空间域总磁场强度
[0089]
本发明的步骤(s5)中，通过有限单元法求解变分问题模型，得到空间波数混合域异常场磁位具体地，采用基于二次插值的一维有限单元法求解，能够保了计算效率和计算精度，且可以利用追赶法实现对角线性方程组的快速求解。
[0090]
根据空间波数混合域异常场磁位求解空间波数混合域异常场磁场强度
[0091]
通过反傅里叶变换将空间波数混合域异常场磁位以及空间波数混合域异常场磁场强度转换为空间域异常场磁位u以及空间域异常场磁场强度ha。
[0092]
空间域总磁场强度h为空间域背景场磁场强度与空间域异常场磁场强度之和，即h＝h0 ha。本发明在初次迭代时，设初次迭代时的空间域异常场磁场强度ha为0，即初次迭代时的空间域总磁场强度h即为空间域背景场磁场强度h0，从而将一维偏微分方程变为一维常微分方程求解，得到第一个空间域异常场磁场强度(上标1代表第1次迭代，即第1次迭代得到的空间域异常场磁场强度)，和h0之和即为第一次迭代后得到的空间域总磁场强度h1。第一次迭代后得到的空间域总磁场强度h1作为下一次迭代时的空间域总磁场强度进行下一次求解，以此规律不断迭代，直至收敛。
[0093]
可以理解，预设的迭代终止条件是指预先设置的模型计算约束条件，用于约束整个模型进行性能计算的过程趋向收敛，以使模型能够输出满足条件的结果。本发明中可以将(s8)中的迭代终止条件设置为：
[0094][0095]
其中hj表示第j次迭代计算得到的空间域总磁场强度，h
j 1
表示第j 1次计算得到的空间域总磁场强度，初次迭代时，空间域背景场磁场强度h0作为初次迭代的空间域总磁场强度h0。
[0096]
当满足以上迭代收敛条件时，迭代停止。
[0097]
当然，实际应用中，本领域技术人员也可基于现有技术、本领域的惯用技术手段或者公知常识，设定其他的迭代终止条件，不局限于本技术上述优选实施例中所述设置的迭代终止条件。
[0098]
下面对本发明提供的二维强磁场数值模拟方法的正确性、精度和效率进行验证。
[0099]
以下操作均利用fortran95语言编程计算实现，电脑配置为：intel core i3-4150 cpu，主频为3.50ghz，内存为12gb。
[0100]
设计沿y轴无限延伸的二维圆柱体异常体模型。背景磁场为h0＝50000nt，异常体磁化率χ＝5si，磁倾角α＝45
°
，磁偏角β＝5
°
。研究区域为：x方向[-1000m，1000m]，z方向[0，1000m]；网格数为400
×
400，水平采样间隔为5m，垂向采样间隔为2.5m；异常圆柱体圆心坐标为(0，500m)，半径为r＝200m，异常体沿y方向上无限延伸r，如图2所示。在单线程情况下，达到收敛条件需迭代9次，用时15.146s。空间域异常场磁场强度的x分量、z分量的解析解和数值解及各观测点的相对误差如图3所示，图3中：(1)为空间域异常场磁场强度ha其x分量解析解和数值解结果图，(2)为空间域异常场磁场强度ha其x分量解析解和数值解的相对误差，(3)为空间域异常场磁场强度ha其z分量解析解和数值解结果图，(4)为空间域异常场磁场强度ha其z分量解析解和数值解的相对误差。
[0101]
本发明考虑了自退磁效应的影响，对于一些二维强磁性体磁异常的刻画更加准确；而且兼顾数值模拟的精度和效率，充分利用了空间域和波数域方法的优势，将空间域磁位满足的二维偏微分方程沿水平方向进行一维傅里叶变换，将其转换成不同波数相互独立的一维常微分方程，运用有限单元法，单元内部采用形函数二次插值，从而对微分方程进行迭代求解，大大提高了计算精度和计算效率，本发明为磁测数据的处理、解释以及磁场正反演奠定了基础。
[0102]
本发明一实施例中提供发明提供一种二维强磁场数值模拟方法，包括：
[0103]
第一模块，用于对沿y轴方向无限延伸的异常体，在xoz平面确定目标区域，所述异常体xoz方向的整个截面包含在所述目标区域中，对目标区域进行网格剖分，得到多个长方形单元，对各长方形单元顶点处的磁化率χ进行赋值，磁化率χ为常数；
[0104]
第二模块，用于根据磁化率张量、空间域背景场磁场强度、空间域异常场磁场强度，得到空间域磁化强度；
[0105]
第三模块，利用一维傅里叶变换将空间域异常场磁位和空间域磁化强度满足的二维偏微分方程转为空间波数混合域一维常微分方程；
[0106]
第四模块，基于空间波数混合域一维常微分方程，并结合设定的空间波数混合域
异常场磁位需满足的边界条件，将空间波数混合域异常场磁位满足的边值问题模型转化为等价的变分问题模型；
[0107]
第五模块，用于求解变分问题模型，得到空间波数混合域异常场磁位；
[0108]
第六模块，基于空间波数混合域异常场磁位，求得空间波数混合域异常场磁场强度；
[0109]
第七模块，通过反傅里叶变换将空间波数混合域异常场磁位以及空间波数混合域异常场磁场强度转换为空间域异常场磁位以及空间域异常场磁场强度；
[0110]
第八模块，判断当前是否满足迭代终止条件，如满足则输出当前计算得到的空间域异常场磁位以及空间域异常场磁场强度。
[0111]
上述各模块功能的实现方法，可以采用前述各实施例中相同的方法实现，在此不再赘述。
[0112]
在本一个实施例中，提供了一种计算机设备，该计算机设备可以是服务器，其内部结构图可以如图4所示。该计算机设备包括通过系统总线连接的处理器、存储器、网络接口和数据库。其中，该计算机设备的处理器用于提供计算和控制能力。该计算机设备的存储器包括非易失性存储介质、内存储器。该非易失性存储介质存储有操作系统、计算机程序和数据库。该内存储器为非易失性存储介质中的操作系统和计算机程序的运行提供环境。该计算机设备的数据库用于存储样本数据。该计算机设备的网络接口用于与外部的终端通过网络连接通信。该计算机程序被处理器执行时以实现上述二维强磁场数值模拟方法。
[0113]
本领域技术人员可以理解，图4中示出的结构，仅仅是与本技术方案相关的部分结构的框图，并不构成对本技术方案所应用于其上的计算机设备的限定，具体的计算机设备可以包括比图中所示更多或更少的部件，或者组合某些部件，或者具有不同的部件布置。
[0114]
在一个实施例中，提供了一种计算机设备，包括存储器和处理器，该存储器存储有计算机程序，该处理器执行计算机程序时实现上述实施例中二维强磁场数值模拟方法的步骤。
[0115]
在一个实施例中，提供了一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现上述实施例中二维强磁场数值模拟方法的步骤。
[0116]
本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的计算机程序可存储于一非易失性计算机可读取存储介质中，该计算机程序在执行时，可包括如上述各方法的实施例的流程。其中，本技术所提供的各实施例中所使用的对存储器、存储、数据库或其它介质的任何引用，均可包括非易失性和/或易失性存储器。非易失性存储器可包括只读存储器(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)等。
[0117]
以上实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。
[0118]
以上所述实施例仅表达了本技术的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本技术构思的前提下，还可以做出若干变形和改进，这些都属于本技术的保护范围。因此，本技术专利的保护范围应以所附权利要求为准。