CN108197389A - 二维强磁性体磁场的快速、高精度数值模拟方法 - Google Patents
二维强磁性体磁场的快速、高精度数值模拟方法 Download PDFInfo
- Publication number
- CN108197389A CN108197389A CN201810006227.5A CN201810006227A CN108197389A CN 108197389 A CN108197389 A CN 108197389A CN 201810006227 A CN201810006227 A CN 201810006227A CN 108197389 A CN108197389 A CN 108197389A
- Authority
- CN
- China
- Prior art keywords
- magnetic field
- dimentional
- ferromagnetic
- directions
- rectangle
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Measuring Magnetic Variables (AREA)
Abstract
本发明提供一种二维强磁性体磁场的快速、高精度数值模拟计算方法，其通过二维强磁性体模型表示、高斯参数设计、离散偏移波数计算、波数域加权系数计算、磁化强度计算、一维离散傅里叶反变换、空间域磁异常计算、迭代收敛判断等步骤，实现了二维强磁性体磁场数值模拟在效率和精度上的统一。本发明解决了目前强磁性体磁场数值模拟方法计算精度低，无法满足大规模强磁测数据精细反演成像的问题，有助于开展大规模强磁测数据二维磁化率精细反演成像、人机交互建模和解释的研究。
Description
技术领域
本发明涉及一种面向航空磁法勘探的数值模拟方法，特别是一种面向二维强磁性体磁法勘探的高效、高精度数值模拟方法。
背景技术
磁法是矿产资源勘探的重要手段，是寻找一些金属矿类的主要方法之一。当今的高精度磁力仪如质子磁力仪、光泵磁力仪等，其灵敏度已分别达到0.01nT和0.001nT。因此，相应的高精度处理转换以及反演解释逐步受到重视。但是矿产大多都是高磁化率的，这给在磁测资料的处理和解释中带来很大困难。不考虑强磁效应的影响，必然会引起计算结果的误差，并严重影响地质解释的效果。
针对复杂磁化率分布条件下强磁性体磁场数值模拟，众多国内外学者进行了研究。文献(Eskola L,Tervo T.Solving the magnetostatic field problem(a case ofhigh susceptibility)by means of the method of subsections.Geoexploration,1980,18(2):79-95.)将磁性空间表示为线性均匀的磁导率和剩磁的子区间，在考虑退磁效应的前提下，采用面积分方法进行了计算，数值计算量要小于体积分方法，但是解的准确性与计算机能够求解的代数方程组的大小有关。
文献(Purss M B J.A new iterative method for computing the magneticfield at high magnetic susceptibilities.Geophysics,2005,70(5):53.)用一组具有不同半径和磁矩的球体来表示非椭球形状地质体，并通过迭代方法来计算各个球体之间的相互作用，从而减小近似估算退磁因子的误差。该方法计算速度快，并能精确地给出复杂形体退磁系数的估算值，但是对于磁化率很高的情况下仍存在较大误差。
文献(Kostrov N P.Calculation of magnetic anomalies caused by 2Dbodies of arbitrary shape with consideration of demagnetization.GeophysicalProspecting,2007,55(1):91-115.)其基于体积分(VIE)方法提出了一种采用三角单元的体积分方法，该方法极大地克服了VIE方法的缺陷，能够十分精确计算相对磁导率在2～20范围内，以及离磁性体不到一个单元大小距离的磁场，但是该方法的计算效率较低。
已有的强磁性体磁场数值模拟方法，对计算机的计算资源要求普遍较高，对于大规模的磁测数据模拟来说，计算速度较慢。因此寻找一种计算效率与计算精度高的强磁性体磁场数值模拟方法，对于实现强磁性体磁测数据二维反演成像具有重要的现实意义。
发明内容
本发明的目的在于提供一种二维强磁性体磁场的快速、高精度数值模拟计算方法，该发明解决了目前强磁性体磁场数值模拟方法计算精度低，计算时间长，无法满足大规模强磁测数据精细反演成像的问题，有助于开展大规模强磁测数据二维磁化率精细反演成像、人机交互建模和解释的研究。
为了实现上述技术目的，本发明的技术方案是：
一种二维强磁性体磁场的快速、高精度数值模拟计算方法，包括以下步骤：
步骤一：二维强磁性体模型表示：
确定目标区域以及目标区域内的二维强磁性体，建立包含所有目标区域的长方形模型，使得目标区域内的二维强磁性体完全包含在该长方形模型中；在长方形模型中，任取一点作为原点，建立二维直角坐标系x0z，利用由多根x方向以及z方向的平行线组成的网格线将该长方形模型剖分为若干个小长方形，长方形模型的x、z方向小长方形的剖分个数分别为Nx、Nz；各个小长方形x方向边长均相同；在长方形模型的z方向上同一层的各小长方形其z方向边长相同，ζ表示长方形模型中的目标区域中任一点的z坐标，ζ的取值范围是[zmin，zmax]，zmin表示目标区域中所有点中的最小z坐标，zmax表示目标区域中所有点中的最大z坐标。在目标区域中的平行于x轴且分别对应不同z坐标的一组平行线为观测线，所有观测线对应的z坐标的取值范围是[zmin，zmax]。
根据每个小长方形内嵌入的目标区域所对应的磁化率的分布，对每个小长方形的磁化率进行相应的赋值，每个小长方体磁化率为常值，不同小长方体的磁化率取值不同，以此刻画任意磁化率分布复杂的二维强磁性体模型；将位于空气部分的小长方体的磁化率值设为零，以此刻画起伏地形。
步骤二：高斯参数设计：
给定x方向的高斯点个数Lx、区间[-1,1]上高斯点ta及高斯系数Aa；其中，高斯点个数Lx＝4，a＝1,2,…,Lx。ta及Aa在高斯点个数确定后可通过查表得到。
步骤三：根据空间域剖分参数和高斯参数，计算离散偏移波数，具体过程如下：
式中，
a＝1,2,…,Lx
kx表示x方向的偏移波数，Δkx表示x方向基波数，Nx表示长方形模型其x方向小长方形的剖分个数，Δx表示小长方形其x方向的边长。
步骤四：波数域加权系数计算：
式中，μ表示真空磁导率，Δzk表示长方形模型其第k层小长方形的z方向边长，i表示虚数单位。zl表示当前计算的观测线所对应的z坐标；二维强磁性体磁场包括x，z两个方向的磁场，分别为Bx(xi,zl)和Bz(xi,zl)；分别表示二维强磁性体磁场x方向的磁场Bx(xi,zl)在x，z两个方向的波数域加权系数；分别表示二维强磁性体磁场z方向的磁场Bz(xi,zl)在x,z两个方向的波数域加权系数，sgn()为符号函数
波数域加权系数计算结果如下：
(1)当zl-ζ＞0时
(2)当zl-ζ＜0时
步骤五：地球主磁场设置：
首先，根据公知地球主磁场模型IGRF，计算各小长方形中心位置(xi,zk)处的地球主磁场二分量Tx(xi,zk)，Tz(xi,zk)；其中：Tx(xi,zk)、Tz(xi,zk)分别表示(xi,zk)处地球主磁场的x、z分量；(xi,zk)表示编号为(i,k)的小长方形的几何中心坐标，i＝1,2,…,Nx，k＝1,2,…,Nz；Nz表示长方形模型其z方向小长方形的剖分个数。
步骤六：设置二维强磁性体磁场初始值，对二维强磁性体磁场进行迭代循环计算，直到计算出的二维强磁性体磁场满足设置的迭代收敛条件为止。
(1)二维强磁性体磁场初始值设置：
给出初始二维强磁性体磁场：Bx(xi,zl)＝0，Bz(xi,zl)＝0；Bx(xi,zl)，Bz(xi,zl)分别表示二维强磁性体磁场的x，z两个方向的磁场分量；
(2)磁化强度计算：
Mx(xi,zk)＝χ(xi,zk)(Tx(xi,zk)+Bx(xi,zl)) (15)
Mz(xi,zk)＝χ(xi,zk)(Tz(xi,zk)+Bz(xi,zl)) (16)
其中：χ(xi,zk)表示编号为(i,k)的小长方形的磁化率值；
(3)采用一维快速傅里叶变换算法(本领域的常规算法)，快速计算波数域磁化强度二分量：
其中，为z方向第k层小长方形波数域磁化强度的x、z分量，k＝1,2,…,Nz，i表示虚数单位，χ(xi,zk)表示编号为(i,k)的小长方形的磁化率值。
(4)根据波数域磁化强度和波数域加权系数，计算目标区域不同观测线上各自对应的波数域磁场
式中，表示z坐标为zl的观测线上x、z方向波数域磁场。
(5)采用一维快速傅里叶反变换算法，快速计算目标区域不同观测线上各自对应的x、z方向的磁场分量
其中，Aa为高斯加权系数。
(6)将(5)中得到的计算结果更新为下一次二维强磁性体磁场迭代循环计算的初始值，按照步骤(2)至(5)中的方法，迭代循环计算二维强磁性体磁场，直到计算出的二维强磁性体磁场满足设置的迭代收敛条件为止。
在二维强磁性体磁场迭代循环计算过程中，每次更新磁化强度公式如下：
(Mx(xi,zk))n＝χ(xi,zk)(Tx(xi,zk)+(Bx(xi,zl))n-1),n＝2,3,4...Nd (23)
(Mz(xi,zk))n＝χ(xi,zk)(Tz(xi,zk)+(Bz(xi,zl))n-1),n＝2,3,4...Nd (24)
其中Nd表示迭代的次数。磁场的迭代循环计算过程中，每循环一次得到一个计算结果；
(Bx(xi,zl))n表示第n次循环计算中得到的z坐标为zl的观测线上对应的x方向的磁场分量；
(Bx(xi,zl))n-1表示第n-1次循环计算中得到的z坐标为zl的观测线上对应的x方向的磁场分量；
(Bz(xi,zl))n表示第n次循环计算中得到的z坐标为zl的观测线上对应的z方向的磁场分量；
(Bz(xi,zl))n-1表示第n-1次循环计算中得到的z坐标为zl的观测线上对应的z方向的磁场分量。
迭代收敛条件如下：
其中rrms(Bx(xi,zl))、rrms(Bz(xi,zl))分别表示二维强磁性体磁场其x，z两个方向的磁场分量Bx(xi,zl)与Bz(xi,zl)的相对均方根误差。
本发明是一个有机整体，包含两个关键环节：第一个是在特殊的模型剖分方式下，给出了精确的长方形组合模型磁场波数域计算公式；第二个是二维强磁性体磁场快速、高精度迭代正演算法。基于上述基础，本发明实现了任意磁化率分布、二维强磁性体磁场数值模拟在计算效率和计算精度上的统一。
与现有技术相比，本发明具有以下优点：
(1)模型剖分方法简单、灵活，很容易刻画任意磁化率分布复杂磁性体以及起伏观测线；
(2)能够实现强磁性体磁场的快速、高精度数值模拟，可以满足大规模强磁性体磁测数据精细反演成像、人机交互建模和解释的需求；
(3)大规模数值模拟时，占用计算机内存少，并行性好，计算效率和计算精度高。
附图说明
图1为本发明的流程图；
图2为二维强磁性体模型示意图；
图3为强磁性体不迭代计算值；
图4为强磁性体迭代计算值
图5为Bx分量迭代收敛趋势图
图6为Bz分量迭代收敛趋势图；
图中符号说明如下：
N：表示迭代收敛的次数；
rrms：表示相对均方根误差。
本发明目的的实现、功能特点及优点将结合实施例，参照附图做进一步说明。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施方式作进一步地详细描述。
如图1所示，图1为本发明一种二维强磁性体磁场的快速、高精度数值模拟计算方法的流程图。本发明包括以下步骤：
步骤一：二维强磁性体模型表示：
如图2所示，图2为二维强磁性体模型示意图。二维强磁性体磁场就是二维强磁性体(即图2中的异常体)产生的，异常体所在的区域叫异常区域。
首先，确定目标区域以及目标区域内的二维强磁性体，建立包含所有目标区域的长方形模型，使得目标区域内的二维强磁性体所在的区域完全包含在该长方形模型中。
在长方形模型中，任取一点作为原点，建立二维直角坐标系x0z，利用由多根x方向以及z方向的平行线组成的网格线将该长方形模型均匀剖分为若干个小长方形，长方形模型的x、z方向小长方形的剖分个数分别为Nx、Nz；各个小长方形x方向边长均相同；在长方形模型的z方向上同一层的各小长方形其z方向边长相同。ζ表示长方形模型中的目标区域中任一点的z坐标，ζ的取值范围是[zmin，zmax]，zmin表示目标区域中所有点中的最小z坐标，zmax表示目标区域中所有点中的最大z坐标；在目标区域中的平行于x轴且分别对应不同z坐标的一组平行线为观测线。
根据每个小长方形内嵌入的目标区域所对应的磁化率的分布，对每个小长方形的磁化率进行相应的赋值，每个小长方体磁化率为常值，不同小长方体的磁化率取值不同，以此刻画任意磁化率分布复杂的二维强磁性体模型；将位于空气部分的小长方体的磁化率值设为零，以此刻画起伏地形。
步骤二：高斯参数设计：
给定x方向的高斯点个数Lx、区间[-1,1]上高斯点ta及高斯系数Aa；本实施例中：高斯点个数Lx＝4，a＝1,2,…,Lx。ta及Aa在高斯点个数确定后可通过查表得到。
步骤三：根据空间域剖分参数和高斯参数，计算离散偏移波数，具体过程如下：
式中，
a＝1,2,…,Lx
kx表示x方向的偏移波数，Δkx表示x方向基波数，Nx表示长方形模型其x方向小长方形的剖分个数，Δx表示小长方形其x方向的边长。
步骤四：波数域加权系数计算：
式中，μ表示真空磁导率，Δzk表示长方形模型其第k层小长方形的z方向边长，i表示虚数单位；zl表示当前计算的观测线所对应的z坐标；二维强磁性体磁场包括x，z两个方向的磁场，分别为Bx(xi,zl)和Bz(xi,zl)；分别表示二维强磁性体磁场x方向的磁场Bx(xi,zl)在x，z两个方向的波数域加权系数；分别表示二维强磁性体磁场z方向的磁场Bz(xi,zl)在x,z两个方向的波数域加权系数，sgn()为符号函数
波数域加权系数计算结果如下：
(1)当zl-ζ＞0时
(2)当zl-ζ＜0时
步骤五：目标区域的地球主磁场设置：
根据公知地球主磁场模型IGRF，计算各小长方形中心位置(xi,zk)处的地球主磁场二分量Tx(xi,zk)，Tz(xi,zk)；其中：Tx(xi,zk)、Tz(xi,zk)分别表示(xi,zk)处地球主磁场的x、z分量；(xi,zk)表示编号为(i,k)的小长方形的几何中心坐标，i＝1,2,…,Nx，k＝1,2,…,Nz；Nz表示长方形模型其z方向小长方形的剖分个数。
步骤六：设置二维强磁性体磁场初始值，对二维强磁性体磁场进行迭代循环计算，直到计算出的二维强磁性体磁场满足设置的迭代收敛条件为止。
(1)二维强磁性体磁场初始值设置：
给出初始二维强磁性体磁场：Bx(xi,zl)＝0，Bz(xi,zl)＝0；Bx(xi,zl)，Bz(xi,zl)分别表示二维强磁性体磁场的x，z两个方向的磁场分量；
(2)磁化强度计算：
Mx(xi,zk)＝χ(xi,zk)(Tx(xi,zk)+Bx(xi,zl)) (40)
Mz(xi,zk)＝χ(xi,zk)(Tz(xi,zk)+Bz(xi,zl)) (41)
(3)采用一维快速傅里叶变换算法(本领域的常规算法)，快速计算波数域磁化强度二分量：
其中，为z方向第k层小长方形波数域磁化强度的x、z分量，k＝1,2,…,Nz，i表示虚数单位，χ(xi,zk)表示编号为(i,k)的小长方形的磁化率值。
(4)根据波数域磁化强度和波数域加权系数，计算目标区域不同观测线上各自对应的波数域磁场
式中，表示z坐标为zl的观测线上的x、z方向波数域磁场分量。
(5)采用一维快速傅里叶反变换算法，快速计算目标区域不同观测线上各自对应的x、z方向的磁场分量
其中，Aa为高斯加权系数。
(6)将(5)中得到的计算结果更新为下一次二维强磁性体磁场迭代循环计算的初始值，按照步骤(2)至(5)中的方法，迭代循环计算二维强磁性体磁场，直到计算出的二维强磁性体磁场满足设置的迭代收敛条件为止。
在二维强磁性体磁场迭代循环计算过程中，每次更新磁化强度公式如下：
(Mx(xi,zk))n＝χ(xi,zk)(Tx(xi,zk)+(Bx(xi,zl))n-1),n＝2,3,4...Nd (48)
(Mz(xi,zk))n＝χ(xi,zk)(Tz(xi,zk)+(Bz(xi,zl))n-1),n＝2,3,4...Nd (49)
其中Nd表示迭代的次数。磁场的迭代循环计算过程中，每循环一次得到一个计算结果；
(Bx(xi,zl))n表示第n次循环计算中得到的z坐标为zl的观测线上对应的x方向的磁场分量；
(Bx(xi,zl))n-1表示第n-1次循环计算中得到的z坐标为zl的观测线上对应的x方向的磁场分量；
(Bz(xi,zl))n表示第n次循环计算中得到的z坐标为zl的观测线上对应的z方向的磁场分量；
(Bz(xi,zl))n-1表示第n-1次循环计算中得到的z坐标为zl的观测线上对应的z方向的磁场分量；
迭代收敛条件设置如下：
其中rrms(Bx(xi,zl))、rrms(Bz(xi,zl))分别表示二维强磁性体磁场其x，z两个方向的磁场分量Bx(xi,zl)与Bz(xi,zl)的相对均方根误差。
下面给出一具体实施例，验证本发明所提出的方法用于计算任意磁化率分布的二维强磁性体磁场的效率和精度。
目标区域有一个长方形强磁性体(即目标区域内的二维强磁性体)，目标区域范围为：x方向从-1000m到1000m，z方向从0m到1000m(z轴向下为正)；长方形强磁性体的展布范围为：x方向从-500m到500m，z方向从200m到700m；磁化率为1SI；目标区域的地球主磁场为45000nT，磁偏角为5度，磁倾角为45度。
将目标区域利用网格线剖分成200×200个大小相同的小长方形，在目标区域中确定观测线，观测线是平行于x轴且分别对应不同z坐标的一组平行线。迭代计算目标区域不同深度(即不同z坐标)观测线磁场。
本发明方法利用Fortran语言编程实现，运行程序所用的个人台式机配置为：CPU-Inter Core i5-4590,主频为3.3GHz,内存为16GB。图3与图4分别为不采用迭代算法与采用迭代算法的等值线图，从图中对比可以看出两者存在明显的差异。图5、图6分别为磁场Bx与磁场Bz迭代收敛趋势图，从图中可以看出，当迭代次数达到15次左右时，两个分量的迭代误差都达到0.1％以下，收敛速度较快，计算时间约为3.37s。
以上所述仅是本发明的优选实施方式，本发明的保护范围并不局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。
Claims (4)
1.一种二维强磁性体磁场的快速、高精度数值模拟计算方法，其特征在于，包括以下步骤：步骤一：二维强磁性体模型表示：
确定目标区域以及目标区域内的二维强磁性体，建立包含所有目标区域的长方形模型，使得目标区域内的二维强磁性体完全包含在该长方形模型中；在长方形模型中，任取一点作为原点，建立二维直角坐标系x0z，利用由多根x方向以及z方向的平行线组成的网格线将该长方形模型剖分为若干个小长方形，长方形模型的x、z方向小长方形的剖分个数分别为Nx、Nz；各个小长方形x方向边长均相同；在长方形模型的z方向上同一层的各小长方形其z方向边长相同，ζ表示长方形模型中的目标区域中任一点的z坐标，ζ的取值范围是[zmin，zmax]，zmin表示目标区域中所有点中的最小z坐标，zmax表示目标区域中所有点中的最大z坐标；在目标区域中的平行于x轴且分别对应不同z坐标的一组平行线为观测线；
根据每个小长方形内嵌入的目标区域所对应的磁化率的分布，对每个小长方形的磁化率进行相应的赋值，每个小长方体磁化率为常值，不同小长方体的磁化率取值不同，以此刻画任意磁化率分布复杂的二维强磁性体模型；将位于空气部分的小长方体的磁化率值设为零，以此刻画起伏地形。
步骤二：高斯参数确定：
给定x方向的高斯点个数Lx、区间[-1,1]上高斯点ta及高斯系数Aa，其中a＝1,2,…,Lx；
步骤三：根据空间域剖分参数和高斯参数，计算离散偏移波数：
式中，
a＝1,2,…,Lx
kx表示x方向的偏移波数，Δkx表示x方向基波数，Nx表示长方形模型其x方向小长方形的剖分个数，Δx表示小长方形其x方向的边长；
步骤四：波数域加权系数计算：
式中，μ表示真空磁导率，Δzk表示长方形模型其第k层小长方形的z方向边长，i表示虚数单位；zl表示当前计算的观测线所对应的z坐标；二维强磁性体磁场包括x，z两个方向的磁场，分别为Bx(xi,zl)和Bz(xi,zl)；分别表示二维强磁性体磁场x方向的磁场Bx(xi,zl)在x，z两个方向的波数域加权系数；分别表示二维强磁性体磁场z方向的磁场Bz(xi,zl)在x,z两个方向的波数域加权系数，sgn()为符号函数
波数域加权系数计算结果如下：
(1)当zl-ζ＞0时
(2)当zl-ζ＜0时
步骤五：地球主磁场设置：
根据公知地球主磁场模型IGRF，计算各小长方形中心位置(xi,zk)处的地球主磁场二分量Tx(xi,zk)，Tz(xi,zk)；其中：Tx(xi,zk)、Tz(xi,zk)分别表示(xi,zk)处地球主磁场的x、z分量；(xi,zk)表示编号为(i,k)的小长方形的几何中心坐标，i＝1,2,…,Nx，k＝1,2,…,Nz；Nz表示长方形模型其z方向小长方形的剖分个数；
步骤六：设置二维强磁性体磁场初始值，对二维强磁性体磁场进行迭代循环计算，直到计算出的二维强磁性体磁场满足设置的迭代收敛条件为止。
2.根据权利要求1所述的二维强磁性体磁场的快速、高精度数值模拟计算方法，其特征在于，步骤二中，高斯点个数Lx＝4。
3.根据权利要求1所述的二维强磁性体磁场的快速、高精度数值模拟计算方法，其特征在于，步骤六的实现方法如下：
(1)二维强磁性体磁场初始值设置：
给出初始二维强磁性体磁场：Bx(xi,zl)＝0，Bz(xi,zl)＝0；Bx(xi,zl)，Bz(xi,zl)分别表示二维强磁性体磁场的x，z两个方向的磁场；
(2)磁化强度计算：
Mx(xi,zk)＝χ(xi,zk)(Tx(xi,zk)+Bx(xi,zl)) (15)
Mz(xi,zk)＝χ(xi,zk)(Tz(xi,zk)+Bz(xi,zl)) (16)
其中：χ(xi,zk)表示编号为(i,k)的小长方形的磁化率值；
(3)采用一维快速傅里叶变换算法，快速计算波数域磁化强度二分量：
其中，为z方向第k层小长方形波数域磁化强度的x、z分量，k＝1,2,…,Nz，i表示虚数单位，χ(xi,zk)表示编号为(i,k)的小长方形的磁化率值；
(4)根据波数域磁化强度和波数域加权系数，计算目标区域不同观测线上各自对应的波数域磁场
式中，表示z坐标为zl的观测线上的x、z方向波数域磁场分量；
(5)采用一维快速傅里叶反变换算法，快速计算目标区域不同观测线上各自对应的x、z方向的磁场分量
其中，Aa为高斯加权系数；
(6)将(5)中得到的计算结果更新为下一次二维强磁性体磁场迭代循环计算的初始值，按照步骤(2)至(5)中的方法，迭代循环计算二维强磁性体磁场，直到计算出的二维强磁性体磁场满足设置的迭代收敛条件为止。
4.根据权利要求3所述的二维强磁性体磁场的快速、高精度数值模拟计算方法，其特征在于，步骤六的(6)中，
在二维强磁性体磁场迭代循环计算过程中，每次更新磁化强度公式如下：
(Mx(xi,zk))n＝χ(xi,zk)(Tx(xi,zk)+(Bx(xi,zl))n-1),n＝2,3,4...Nd (23)
(Mz(xi,zk))n＝χ(xi,zk)(Tz(xi,zk)+(Bz(xi,zl))n-1),n＝2,3,4...Nd (24)
其中Nd表示迭代的次数；磁场的迭代循环计算过程中，每循环一次得到一个计算结果；
(Bx(xi,zl))n表示第n次循环计算中得到的z坐标为zl的观测线上对应的x方向的磁场分量；
(Bx(xi,zl))n-1表示第n-1次循环计算中得到的z坐标为zl的观测线上对应的x方向的磁场分量；
(Bz(xi,zl))n表示第n次循环计算中得到的z坐标为zl的观测线上对应的z方向的磁场分量；
(Bz(xi,zl))n-1表示第n-1次循环计算中得到的z坐标为zl的观测线上对应的z方向的磁场分量；
迭代收敛条件设置如下：
其中rrms(Bx(xi,zl))、rrms(Bz(xi,zl))分别表示二维强磁性体磁场其x，z两个方向的磁场分量Bx(xi,zl)与Bz(xi,zl)的相对均方根误差。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810006227.5A CN108197389A (zh) | 2018-01-04 | 2018-01-04 | 二维强磁性体磁场的快速、高精度数值模拟方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810006227.5A CN108197389A (zh) | 2018-01-04 | 2018-01-04 | 二维强磁性体磁场的快速、高精度数值模拟方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN108197389A true CN108197389A (zh) | 2018-06-22 |
Family
ID=62587771
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810006227.5A Pending CN108197389A (zh) | 2018-01-04 | 2018-01-04 | 二维强磁性体磁场的快速、高精度数值模拟方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108197389A (zh) |
Cited By (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109254327A (zh) * | 2018-10-30 | 2019-01-22 | 桂林理工大学 | 三维强磁性体的勘探方法及勘探系统 |
| CN110543611A (zh) * | 2019-08-15 | 2019-12-06 | 桂林理工大学 | 低纬度磁异常数据磁化极计算方法和装置 |
| CN111967169A (zh) * | 2020-10-21 | 2020-11-20 | 中南大学 | 二度体重力异常积分解数值模拟方法和装置 |
| CN112287534A (zh) * | 2020-10-21 | 2021-01-29 | 中南大学 | 基于nufft的二维磁异常快速正演模拟方法和装置 |
| CN113640886A (zh) * | 2021-08-12 | 2021-11-12 | 桂林理工大学 | 强磁性二度体的勘探方法和勘探系统 |
| CN113640887A (zh) * | 2021-08-12 | 2021-11-12 | 桂林理工大学 | 复杂强磁性体的航空勘探方法和勘探系统 |
| CN113656750A (zh) * | 2021-10-20 | 2021-11-16 | 中南大学 | 基于空间波数域的强磁介质的磁感应强度计算方法 |
| CN113962077A (zh) * | 2021-10-20 | 2022-01-21 | 中南大学 | 三维各向异性强磁场数值模拟方法、装置、设备及介质 |
-
2018
- 2018-01-04 CN CN201810006227.5A patent/CN108197389A/zh active Pending
Cited By (14)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109254327A (zh) * | 2018-10-30 | 2019-01-22 | 桂林理工大学 | 三维强磁性体的勘探方法及勘探系统 |
| CN109254327B (zh) * | 2018-10-30 | 2020-11-20 | 桂林理工大学 | 三维强磁性体的勘探方法及勘探系统 |
| CN110543611A (zh) * | 2019-08-15 | 2019-12-06 | 桂林理工大学 | 低纬度磁异常数据磁化极计算方法和装置 |
| CN110543611B (zh) * | 2019-08-15 | 2022-11-25 | 桂林理工大学 | 低纬度磁异常数据磁化极计算方法和装置 |
| CN112287534B (zh) * | 2020-10-21 | 2022-05-13 | 中南大学 | 基于nufft的二维磁异常快速正演模拟方法和装置 |
| CN112287534A (zh) * | 2020-10-21 | 2021-01-29 | 中南大学 | 基于nufft的二维磁异常快速正演模拟方法和装置 |
| CN111967169A (zh) * | 2020-10-21 | 2020-11-20 | 中南大学 | 二度体重力异常积分解数值模拟方法和装置 |
| CN113640886A (zh) * | 2021-08-12 | 2021-11-12 | 桂林理工大学 | 强磁性二度体的勘探方法和勘探系统 |
| CN113640887A (zh) * | 2021-08-12 | 2021-11-12 | 桂林理工大学 | 复杂强磁性体的航空勘探方法和勘探系统 |
| CN113640886B (zh) * | 2021-08-12 | 2023-08-29 | 桂林理工大学 | 强磁性二度体的勘探方法和勘探系统 |
| CN113640887B (zh) * | 2021-08-12 | 2023-09-12 | 桂林理工大学 | 复杂强磁性体的航空勘探方法和勘探系统 |
| CN113656750A (zh) * | 2021-10-20 | 2021-11-16 | 中南大学 | 基于空间波数域的强磁介质的磁感应强度计算方法 |
| CN113962077A (zh) * | 2021-10-20 | 2022-01-21 | 中南大学 | 三维各向异性强磁场数值模拟方法、装置、设备及介质 |
| CN113962077B (zh) * | 2021-10-20 | 2022-08-09 | 中南大学 | 三维各向异性强磁场数值模拟方法、装置、设备及介质 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108197389A (zh) | 二维强磁性体磁场的快速、高精度数值模拟方法 | |
| CN106777598B (zh) | 任意磁化率分布复杂磁性体磁场梯度张量数值模拟方法 | |
| CN102798898B (zh) | 大地电磁场非线性共轭梯度三维反演方法 | |
| CN107024723B (zh) | 一种二度体磁场数值计算方法 | |
| CN109254327B (zh) | 三维强磁性体的勘探方法及勘探系统 | |
| CN111856599B (zh) | 一种基于pde的磁测数据等效源化极与类型转换方法 | |
| Jahandari et al. | Comparison between staggered grid finite–volume and edge–based finite–element modelling of geophysical electromagnetic data on unstructured grids | |
| CN109725361B (zh) | 基于磁梯度张量不变量的一种磁性目标定位方法 | |
| CN104750917A (zh) | 分层介质粗糙面电磁散射系数的确定方法 | |
| CN110346835A (zh) | 大地电磁的正演方法、正演系统、存储介质及电子设备 | |
| Liu et al. | Three-dimensional magnetotellurics modeling using edgebased finite-element unstructured meshes | |
| Ren et al. | Recursive analytical formulae of gravitational fields and gradient tensors for polyhedral bodies with polynomial density contrasts of arbitrary non-negative integer orders | |
| CN107748834B (zh) | 一种计算起伏观测面磁场的快速、高精度数值模拟方法 | |
| CN109101463A (zh) | 一种广域层状大地格林函数的多精度求解方法 | |
| Jeshvaghani et al. | Two-dimensional geomagnetic forward modeling using adaptive finite element method and investigation of the topographic effect | |
| CN114004127A (zh) | 二维主轴各向异性强磁场数值模拟方法、装置、设备及介质 | |
| CN116165722A (zh) | 一种采用高斯牛顿法的回线源瞬变电磁三维快速反演方法 | |
| CN113673163B (zh) | 一种三维磁异常数快速正演方法、装置和计算机设备 | |
| Sheng et al. | The improved residual node density based gravity forward method and its application | |
| CN103559376B (zh) | 基于等效磁矩的潜航平台磁异特征场表征方法 | |
| CN113640886B (zh) | 强磁性二度体的勘探方法和勘探系统 | |
| Zhao et al. | Sparse matrix canonical grid method for three-dimension rough surface | |
| Zhao et al. | An Accelerated 3D-Magnetostatic Field Computation Method With Tetrahedral Surface Integrals | |
| WEI et al. | Forward modeling and inversion of 3‐D cross‐hole electromagnetic fields | |
| CN113139289B (zh) | 一种基于积分方程的退磁效应下磁测数据的正反演方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| RJ01 | Rejection of invention patent application after publication | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20180622 |