二维强磁性体磁场的快速、高精度数值模拟方法技术

本发明专利技术提供一种二维强磁性体磁场的快速、高精度数值模拟计算方法，其通过二维强磁性体模型表示、高斯参数设计、离散偏移波数计算、波数域加权系数计算、磁化强度计算、一维离散傅里叶反变换、空间域磁异常计算、迭代收敛判断等步骤，实现了二维强磁性体磁场数值模拟在效率和精度上的统一。本发明专利技术解决了目前强磁性体磁场数值模拟方法计算精度低，无法满足大规模强磁测数据精细反演成像的问题，有助于开展大规模强磁测数据二维磁化率精细反演成像、人机交互建模和解释的研究。

【技术实现步骤摘要】
二维强磁性体磁场的快速、高精度数值模拟方法
本专利技术涉及一种面向航空磁法勘探的数值模拟方法，特别是一种面向二维强磁性体磁法勘探的高效、高精度数值模拟方法。
技术介绍
磁法是矿产资源勘探的重要手段，是寻找一些金属矿类的主要方法之一。当今的高精度磁力仪如质子磁力仪、光泵磁力仪等，其灵敏度已分别达到0.01nT和0.001nT。因此，相应的高精度处理转换以及反演解释逐步受到重视。但是矿产大多都是高磁化率的，这给在磁测资料的处理和解释中带来很大困难。不考虑强磁效应的影响，必然会引起计算结果的误差，并严重影响地质解释的效果。针对复杂磁化率分布条件下强磁性体磁场数值模拟，众多国内外学者进行了研究。文献(EskolaL,TervoT.Solvingthemagnetostaticfieldproblem(acaseofhighsusceptibility)bymeansofthemethodofsubsections.Geoexploration,1980,18(2):79-95.)将磁性空间表示为线性均匀的磁导率和剩磁的子区间，在考虑退磁效应的前提下，采用面积分方法进行了计算，数值计算量要小于体积分方法，但是解的准确性与计算机能够求解的代数方程组的大小有关。文献(PurssMBJ.Anewiterativemethodforcomputingthemagneticfieldathighmagneticsusceptibilities.Geophysics,2005,70(5):53.)用一组具有不同半径和磁矩的球体来表示非椭球形状地质体，并通过迭代方法来计算各个球体之间的相互作用，从而减小近似估算退磁因子的误差。该方法计算速度快，并能精确地给出复杂形体退磁系数的估算值，但是对于磁化率很高的情况下仍存在较大误差。文献(KostrovNP.Calculationofmagneticanomaliescausedby2Dbodiesofarbitraryshapewithconsiderationofdemagnetization.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,…,Lxkx表示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方向的本文档来自技高网...

【技术保护点】
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；步骤三：根据空间域剖分参数和高斯参数，计算离散偏移波数：...

【技术特征摘要】
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,…,Lxkx表示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所述的二维...

【专利技术属性】
技术研发人员：李昆，陈龙伟，陈轻蕊，强建科，戴世坤，张钱江，赵东东，
申请(专利权)人：中南大学，
类型：发明
国别省市：湖南,43