CN113656976A - 一种二维磁梯度张量快速数值模拟方法、装置和设备 - Google Patents

一种二维磁梯度张量快速数值模拟方法、装置和设备 Download PDF

Info

Publication number
CN113656976A
CN113656976A CN202110978629.3A CN202110978629A CN113656976A CN 113656976 A CN113656976 A CN 113656976A CN 202110978629 A CN202110978629 A CN 202110978629A CN 113656976 A CN113656976 A CN 113656976A
Authority
CN
China
Prior art keywords
dimensional
grid
magnetic gradient
magnetic
kernel function
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.)
Granted
Application number
CN202110978629.3A
Other languages
English (en)
Other versions
CN113656976B (zh
Inventor
王旭龙
柳建新
范平阳
郭友军
王永斐
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Central South University
Original Assignee
Central South University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Central South University filed Critical Central South University
Priority to CN202110978629.3A priority Critical patent/CN113656976B/zh
Publication of CN113656976A publication Critical patent/CN113656976A/zh
Application granted granted Critical
Publication of CN113656976B publication Critical patent/CN113656976B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

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
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/30—Assessment of water resources

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Software Systems (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • Computing Systems (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • Measuring Magnetic Variables (AREA)

Abstract

本申请涉及一种二维磁梯度张量快速数值模拟方法、装置、计算机设备和存储介质。所述方法包括：通过模型构建、磁化强度计算、磁梯度张量计算公式离散、核函数单元积分系数计算、核函数单元积分系数矩阵和磁化强度的快速相乘、二维网格单元的磁梯度张量计算等步骤实现了二维磁梯度张量快速数值模拟模拟。本发明解决了现有二维磁梯度张量正演方法计算精度和计算效率低、无法实现复杂地质体的高效、精细化磁法勘探的问题。

Description

一种二维磁梯度张量快速数值模拟方法、装置和设备
技术领域
本申请涉及计算机技术领域，特别是涉及一种二维磁梯度张量快速数值模拟方法、装置、计算机设备和存储介质。
背景技术
磁法勘探的定量解释依赖于正演的计算精度。根据场源的产状，地质解释人员可以采用二维模型或者三维模型来模拟不同场源的分布。其中，对于沿着走向方向变化不大的线性地质体引起的磁异常，可以利用二度体模型进行解释。尽管在客观上不存在真正的二度体，但三维模型通常需要较大的内存需求以及庞大的计算量，因此现实中常用走向方向尺度远大于垂直走向尺度的地质体用二度体代替，大大减少了计算时间，相应的反演也容易实现。
目前，针对二维磁梯度张量数值模拟一般采用解析方法和数值方法。文献(Jia,Z,Meng,L.S.Magnetic gradients produced by a 2D homogeneous polygonalsourceMagnetic gradients produced by a 2D polygon.Geophysics,2009.74(1)：L1-L6.)推导了多边形截面二度体磁梯度张量的解析表达式，并对解析式中的奇异性进行了分析，有效改进了解析式法的稳定性。解析式法计算精度高，但是针对任意密度分布、任意截面形状的地质体难以推导其解析表达式，且当进行大规模计算时其计算效率较低。文献(Jeshvaghani,M.S.,Darijani,M.Two-dimensional geomagnetic forward modelingusing adaptive finite element method and investigation of the topographiceffect.Journal of Applied Geophysics,2014,105,169-179.)采用自适应有限单元法实现了基于复杂模型及起伏地形的二度体磁异常数值模拟。该方法适用于任意截面形状、任意磁化率分布的复杂地质体的数值模拟，适应性强，并且该方法在复杂模型边界可以网格适当加密，以提高模拟精度，但随着网格节点的增加，形成的大型稀疏矩阵较大，求解效率低。异常体模拟的现有技术存在效率不高的问题。
发明内容
基于此，有必要针对上述技术问题，提供一种能够提高异常体模拟效率的二维磁梯度张量快速数值模拟方法、装置、计算机设备和存储介质。
一种二维磁梯度张量快速数值模拟方法，所述方法包括：
根据二维异常体的大小信息，构建计算区域的二维模型，将所述计算区域沿x、z方向进行均匀网格剖分，得到多个网格单元，根据所述二维异常体的截面形状和磁化率分布，对所述网格单元的磁化率进行赋值；
根据已知的地球正常场，计算所述网格单元地球主磁场的x，z分量，根据所述网格单元地球主磁场的x，z分量和所述网格单元的磁化率，得到所述网格单元的磁化强度；
根据所述二维模型对预先构建的二维磁梯度张量公式进行离散，得到二维磁梯度的解析表达式，根据所述二维磁梯度的解析表达式，得到网格单元的核函数积分系数；
将所述网格单元的磁化强度和所述核函数积分系数表示为矩阵形式，并进行扩展，得到磁化强度扩展矩阵和核函数积分系数扩展矩阵；
通过二维离散傅里叶变换将所述核函数积分系数扩展矩阵和所述磁化强度扩展矩阵的二维卷积计算转换为乘积计算，再进行傅里叶反变换，得到所述计算区域的磁梯度张量。
在其中一个实施例中，还包括：根据已知的地球正常场，计算所述网格单元地球主磁场的x，z分量，根据所述网格单元地球主磁场的x，z分量和所述网格单元的磁化率，得到所述网格单元的磁化强度为：
Mx(xi,zj)＝χ(xi,zj)Tx(xi,zj)
Mz(xi,zj)＝χ(xi,zj)Tz(xi,zj)
其中，Tx(xi,zj)为所述网格单元地球主磁场的x分量，Tz(xi,zj)为所述网格单元地球主磁场的z分量，χ(xi,zj)表示编号为(xi,zj)单元的磁化率值，Mx(xi,zj)、Mz(xi,zj)分别表示(xi,zj)处磁化强度的x、z分量。
在其中一个实施例中，还包括：根据所述二维模型对预先构建的二维磁梯度张量公式进行离散，得到二维磁梯度的解析表达式为：
Figure BDA0003228312330000031
Figure BDA0003228312330000032
其中，Bxz表示磁梯度张量水平分量，Bzz表示磁梯度张量垂直分量，(x,z)表示观测点坐标，(x′,z′)表示源点坐标(xm,zn)表示编号为(m,n)单元的中心坐标，(x′i,z′j)表示编号为(i,j)单元的中心坐标；Mx(x′i,z′j)、Mz(x′i,z′j)表示编号为(i,j)单元的水平和垂向磁化强度；fx(xm-x′i,zn-z′j)，fz(xm-x′i,zn-z′j)表示观测点在(xm,zn)处针对编号为(i,j)单元的核函数积分系数。
在其中一个实施例中，还包括：根据所述二维磁梯度的解析表达式，得到网格单元的核函数积分系数为：
Figure BDA0003228312330000033
Figure BDA0003228312330000041
其中，μ0表示真空磁导率，Δx、Δz分别表示x方向和z方向单元网格间隔。
在其中一个实施例中，还包括：将所述网格单元的所述核函数积分系数表示为矩阵形式为：
Figure BDA0003228312330000042
将
Figure BDA0003228312330000043
扩展为
Figure BDA0003228312330000044
得到核函数积分系数扩展矩阵为：
Figure BDA0003228312330000045
将所述网格单元的磁化强度
Figure BDA0003228312330000046
扩展为
Figure BDA0003228312330000047
其中扩展的其他元素用零补齐，得到磁化强度扩展矩阵为：
Figure BDA0003228312330000048
在其中一个实施例中，还包括：将磁梯度张量表示为所述核函数积分系数和所述磁化强度的卷积：
Figure BDA0003228312330000051
将卷积
Figure BDA0003228312330000052
扩展为所述核函数积分系数扩展矩阵和所述磁化强度扩展矩阵的二维卷积：
Figure BDA0003228312330000053
通过二维离散傅里叶变换将所述核函数积分系数扩展矩阵和所述磁化强度扩展矩阵的二维卷积计算转换为乘积计算：
Figure BDA0003228312330000054
其中，F表示二维离散傅里叶正变换算子；
再进行傅里叶反变换，得到所述计算区域的磁梯度张量为：
Figure BDA0003228312330000055
其中，F-1表示二维离散傅里叶反变换算子；
Figure BDA0003228312330000056
表示提取矩阵的前Nz行，Nx列。
在其中一个实施例中，还包括：不同网格单元的磁化率值不同。
一种二维磁梯度张量快速数值模拟装置，所述装置包括：
网格剖分模块，用于根据二维异常体的大小信息，构建计算区域的二维模型，将所述计算区域沿x、z方向进行均匀网格剖分，得到多个网格单元，根据所述二维异常体的截面形状和磁化率分布，对所述网格单元的磁化率进行赋值；
磁化强度确定模块，用于根据已知的地球正常场，计算所述网格单元地球主磁场的x，z分量，根据所述网格单元地球主磁场的x，z分量和所述网格单元的磁化率，得到所述网格单元的磁化强度；
核函数积分系数确定模块，用于根据所述二维模型对预先构建的二维磁梯度张量公式进行离散，得到二维磁梯度的解析表达式，根据所述二维磁梯度的解析表达式，得到网格单元的核函数积分系数；
矩阵扩展模块，用于将所述网格单元的磁化强度和所述核函数积分系数表示为矩阵形式，并进行扩展，得到磁化强度扩展矩阵和核函数积分系数扩展矩阵；
磁梯度张量确定模块，用于通过二维离散傅里叶变换将所述核函数积分系数扩展矩阵和所述磁化强度扩展矩阵的二维卷积计算转换为乘积计算，再进行傅里叶反变换，得到所述计算区域的磁梯度张量。
一种计算机设备，包括存储器和处理器，所述存储器存储有计算机程序，所述处理器执行所述计算机程序时实现以下步骤：
根据二维异常体的大小信息，构建计算区域的二维模型，将所述计算区域沿x、z方向进行均匀网格剖分，得到多个网格单元，根据所述二维异常体的截面形状和磁化率分布，对所述网格单元的磁化率进行赋值；
根据已知的地球正常场，计算所述网格单元地球主磁场的x，z分量，根据所述网格单元地球主磁场的x，z分量和所述网格单元的磁化率，得到所述网格单元的磁化强度；
根据所述二维模型对预先构建的二维磁梯度张量公式进行离散，得到二维磁梯度的解析表达式，根据所述二维磁梯度的解析表达式，得到网格单元的核函数积分系数；
将所述网格单元的磁化强度和所述核函数积分系数表示为矩阵形式，并进行扩展，得到磁化强度扩展矩阵和核函数积分系数扩展矩阵；
通过二维离散傅里叶变换将所述核函数积分系数扩展矩阵和所述磁化强度扩展矩阵的二维卷积计算转换为乘积计算，再进行傅里叶反变换，得到所述计算区域的磁梯度张量。
一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现以下步骤：
根据二维异常体的大小信息，构建计算区域的二维模型，将所述计算区域沿x、z方向进行均匀网格剖分，得到多个网格单元，根据所述二维异常体的截面形状和磁化率分布，对所述网格单元的磁化率进行赋值；
根据已知的地球正常场，计算所述网格单元地球主磁场的x，z分量，根据所述网格单元地球主磁场的x，z分量和所述网格单元的磁化率，得到所述网格单元的磁化强度；
根据所述二维模型对预先构建的二维磁梯度张量公式进行离散，得到二维磁梯度的解析表达式，根据所述二维磁梯度的解析表达式，得到网格单元的核函数积分系数；
将所述网格单元的磁化强度和所述核函数积分系数表示为矩阵形式，并进行扩展，得到磁化强度扩展矩阵和核函数积分系数扩展矩阵；
通过二维离散傅里叶变换将所述核函数积分系数扩展矩阵和所述磁化强度扩展矩阵的二维卷积计算转换为乘积计算，再进行傅里叶反变换，得到所述计算区域的磁梯度张量。
上述二维磁梯度张量快速数值模拟方法、装置、计算机设备和存储介质，通过根据二维异常体的大小信息，构建计算区域的二维模型，将计算区域沿x、z方向进行均匀网格剖分，得到多个网格单元，根据已知的地球正常场，计算网格单元地球主磁场的x，z分量，根据网格单元地球主磁场的x，z分量和设置的网格单元的磁化率，得到网格单元的磁化强度；根据二维模型对预先构建的二维磁梯度张量公式进行离散，得到二维磁梯度的解析表达式以及网格单元的核函数积分系数；通过二维离散傅里叶变换将核函数积分系数扩展矩阵和磁化强度扩展矩阵的二维卷积计算转换为乘积计算，再进行傅里叶反变换，得到计算区域的磁梯度张量。本申请通过模型构建、磁化强度计算、磁梯度张量计算公式离散、核函数单元积分系数计算、核函数单元积分系数矩阵和磁化强度的快速相乘、二维网格单元的磁梯度张量计算等步骤实现了二维磁梯度张量快速数值模拟模拟。
附图说明
图1为一个实施例中二维磁梯度张量快速数值模拟方法的流程示意图；
图2为另一个实施例中二维磁梯度张量快速数值模拟方法的流程示意图；
图3为一个具体实施例中构建的二维模型示意图；
图4为一个具体实施例中解析解计算的磁梯度张量水平分量平面等值线图；
图5为一个具体实施例中本发明方法计算的磁梯度张量水平分量平面等值线图；
图6为一个具体实施例中磁梯度张量水平分量相对误差平面等值线图；
图7为一个具体实施例中磁梯度张量水平分量地面测线相对误差曲线图；
图8为一个具体实施例中解析解计算的磁梯度张量垂向分量平面等值线图；
图9为一个具体实施例中本发明方法计算的磁梯度张量垂向分量平面等值线图；
图10为一个具体实施例中磁梯度张量垂向分量相对误差平面等值线图；
图11为一个具体实施例中磁梯度张量垂向分量地面测线相对误差曲线图；
图12为一个实施例中二维磁梯度张量快速数值模拟装置的结构框图；
图13为一个实施例中计算机设备的内部结构图。
具体实施方式
为了使本申请的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本申请进行进一步详细说明。应当理解，此处描述的具体实施例仅仅用以解释本申请，并不用于限定本申请。
本申请提供的二维磁梯度张量快速数值模拟方法，可以应用于如下应用环境中。其中，终端执行一种二维磁梯度张量快速数值模拟方法。通过根据二维异常体的大小信息，构建计算区域的二维模型，将计算区域沿x、z方向进行均匀网格剖分，得到多个网格单元，根据已知的地球正常场，计算网格单元地球主磁场的x，z分量，根据网格单元地球主磁场的x，z分量和设置的网格单元的磁化率，得到网格单元的磁化强度；根据二维模型对预先构建的二维磁梯度张量公式进行离散，得到二维磁梯度的解析表达式以及网格单元的核函数积分系数；通过二维离散傅里叶变换将核函数积分系数扩展矩阵和磁化强度扩展矩阵的二维卷积计算转换为乘积计算，再进行傅里叶反变换，得到计算区域的磁梯度张量。其中，终端可以但不限于是各种个人计算机、笔记本电脑、智能手机和便携式可穿戴设备。
在一个实施例中，如图1所示，提供了一种二维磁梯度张量快速数值模拟方法，包括以下步骤：
步骤102，根据二维异常体的大小信息，构建计算区域的二维模型，将计算区域沿x、z方向进行均匀网格剖分，得到多个网格单元，根据二维异常体的截面形状和磁化率分布，对网格单元的磁化率进行赋值。
根据二维异常体的大小，构建包含所有目标区域的二维模型，确定二维模型在x，z方向的起始位置；其次，将整个研究区域沿着平行于x，z坐标轴进行均匀网格剖分；最后，根据二维异常体截面形状和磁化率分布，对网格剖分的每个单元进行赋值，每个单元的磁化率是一个常数，不同网格单元的磁化率值不同，以此刻画任意磁化率分布的二维磁异常模型。
步骤104，根据已知的地球正常场，计算网格单元地球主磁场的x，z分量，根据网格单元地球主磁场的x，z分量和网格单元的磁化率，得到网格单元的磁化强度。
步骤106，根据二维模型对预先构建的二维磁梯度张量公式进行离散，得到二维磁梯度的解析表达式，根据二维磁梯度的解析表达式，得到网格单元的核函数积分系数。
步骤108，将网格单元的磁化强度和核函数积分系数表示为矩阵形式，并进行扩展，得到磁化强度扩展矩阵和核函数积分系数扩展矩阵。
步骤110，通过二维离散傅里叶变换将核函数积分系数扩展矩阵和磁化强度扩展矩阵的二维卷积计算转换为乘积计算，再进行傅里叶反变换，得到计算区域的磁梯度张量。
上述二维磁梯度张量快速数值模拟方法中，通过根据二维异常体的大小信息，构建计算区域的二维模型，将计算区域沿x、z方向进行均匀网格剖分，得到多个网格单元，根据已知的地球正常场，计算网格单元地球主磁场的x，z分量，根据网格单元地球主磁场的x，z分量和设置的网格单元的磁化率，得到网格单元的磁化强度；根据二维模型对预先构建的二维磁梯度张量公式进行离散，得到二维磁梯度的解析表达式以及网格单元的核函数积分系数；通过二维离散傅里叶变换将核函数积分系数扩展矩阵和磁化强度扩展矩阵的二维卷积计算转换为乘积计算，再进行傅里叶反变换，得到计算区域的磁梯度张量。本申请通过模型构建、磁化强度计算、磁梯度张量计算公式离散、核函数单元积分系数计算、核函数单元积分系数矩阵和磁化强度的快速相乘、二维网格单元的磁梯度张量计算等步骤实现了二维磁梯度张量快速数值模拟模拟。
在其中一个实施例中，还包括：根据已知的地球正常场，计算网格单元地球主磁场的x，z分量，根据网格单元地球主磁场的x，z分量和网格单元的磁化率，得到网格单元的磁化强度为：
Mx(xi,zj)＝χ(xi,zj)Tx(xi,zj)
Mz(xi,zj)＝χ(xi,zj)Tz(xi,zj)
其中，Tx(xi,zj)为网格单元地球主磁场的x分量，Tz(xi,zj)为网格单元地球主磁场的z分量，χ(xi,zj)表示编号为(xi,zj)单元的磁化率值，Mx(xi,zj)、Mz(xi,zj)分别表示(xi,zj)处磁化强度的x、z分量。
在其中一个实施例中，还包括：对二维梯度张量公式进行离散，得
Figure BDA0003228312330000101
Figure BDA0003228312330000102
式中，Bxz表示磁梯度张量水平分量，Bzz表示磁梯度张量垂直分量，μ0表示真空磁导率，(x,z)表示观测点坐标，(x′,z′)表示源点坐标。
进一步地，得到二维磁梯度的解析表达式为：
Figure BDA0003228312330000103
Figure BDA0003228312330000104
其中，Bxz表示磁梯度张量水平分量，Bzz表示磁梯度张量垂直分量，(x,z)表示观测点坐标，(x′,z′)表示源点坐标(xm,zn)表示编号为(m,n)单元的中心坐标，(x′i,z′j)表示编号为(i,j)单元的中心坐标；Mx(x′i,z′j)、Mz(x′i,z′j)表示编号为(i,j)单元的水平和垂向磁化强度；fx(xm-x′i,zn-z′j)，fz(xm-x′i,zn-z′j)表示观测点在(xm,zn)处针对编号为(i,j)单元的核函数积分系数。
给定观测点坐标(xm,zn)和编号为(i,j)的单元坐标(x′m,z′n)，根据二维磁梯度的解析表达式，得到网格单元的核函数积分系数为：
Figure BDA0003228312330000111
Figure BDA0003228312330000112
其中，μ0表示真空磁导率，Δx、Δz分别表示x方向和z方向单元网格间隔。
在其中一个实施例中，还包括：将网格单元的核函数积分系数表示为矩阵形式为：
Figure BDA0003228312330000121
将
Figure BDA0003228312330000122
扩展为
Figure BDA0003228312330000123
得到核函数积分系数扩展矩阵为：
Figure BDA0003228312330000124
将网格单元的磁化强度
Figure BDA0003228312330000125
扩展为
Figure BDA0003228312330000126
其中扩展的其他元素用零补齐，得到磁化强度扩展矩阵为：
Figure BDA0003228312330000127
在其中一个实施例中，还包括：将磁梯度张量表示为核函数积分系数和磁化强度的卷积：
Figure BDA0003228312330000128
将卷积
Figure BDA0003228312330000129
扩展为所述核函数积分系数扩展矩阵和所述磁化强度扩展矩阵的二维卷积：
Figure BDA00032283123300001210
通过二维离散傅里叶变换将核函数积分系数扩展矩阵和磁化强度扩展矩阵的二维卷积计算转换为乘积计算：
Figure BDA00032283123300001211
其中，F表示二维离散傅里叶正变换算子；
再进行傅里叶反变换，得到计算区域的磁梯度张量为：
Figure BDA0003228312330000131
其中，F-1表示二维离散傅里叶反变换算子；
Figure BDA0003228312330000132
表示提取矩阵的前Nz行，Nx列。
二维磁梯度张量满足的二维积分，离散化后得到的系数矩阵维数非常高，这就造成两方面困难：一是存储系数矩阵需要耗费大量计算机内存；二是对系数矩阵的运算(主要指其与向量相乘运算)需要花费很长时间。本发明方法只需对系数矩阵计算一次，并采用二维离散傅里叶变换一次可以得到整个二维磁梯度张量异常。以一个具体实施例为例：本发明方法网格剖分为300×300，常规空间域计算整个二维剖面所占用内存(假设存储一个数需要8字节)为，即3002×3002×8B＝60.3GB，而本发明方法所需内存仅需要599×599×8B＝0.0027GB，并采用快速傅里叶变换实现矩阵与向量的快速相乘。占用计算机内存小，计算时间大大缩短。
在其中一个实施例中，还包括：不同网格单元的磁化率值不同。
应该理解的是，虽然图1的流程图中的各个步骤按照箭头的指示依次显示，但是这些步骤并不是必然按照箭头指示的顺序依次执行。除非本文中有明确的说明，这些步骤的执行并没有严格的顺序限制，这些步骤可以以其它的顺序执行。而且，图1中的至少一部分步骤可以包括多个子步骤或者多个阶段，这些子步骤或者阶段并不必然是在同一时刻执行完成，而是可以在不同的时刻执行，这些子步骤或者阶段的执行顺序也不必然是依次进行，而是可以与其它步骤或者其它步骤的子步骤或者阶段的至少一部分轮流或者交替地执行。
在另一个实施例中，如图2所示，提供了一种二维磁梯度张量快速数值模拟方法，包括：
步骤一：二维模型构建：根据探测目标大小设定网格尺寸，根据目标体磁化率分布，设置剖分单元磁化率值；
步骤二：磁化强度计算：根据地球正常场数据，给定目标区域背景磁场，根据磁化率分布和主磁场，计算磁化强度各分量值；
步骤三：磁梯度张量计算公式离散：根据上述模型参数对二维磁梯度张量计算公式进行离散；
步骤四：核函数单元积分系数计算：对每个网格单元进行积分，计算核函数单元积分系数；
步骤五：核函数系数矩阵与磁化强度快速相乘：采用二维离散傅里叶变换实现核函数与磁化强度快速乘积；
步骤六：二维剖面磁梯度张量计算：对核函数系数矩阵与磁化强度快速相乘结果进行二维傅里叶反变换，计算得到整个剖面磁梯度张量异常。
在一个具体实施例中，目标区域有一个截面为规则矩形的二度体模型如图3所示，计算区域范围为：x方向从-500m到500m，z方向从0m到500m(z轴垂直向下为正)，水平网格间隔均为5m，垂向网格间隔为2.5m，水平和垂向网格个数为200×200，截面为矩形的异常体分布范围为：x方向从-100m到100m，z方向从200m到300m，异常区域磁化率为0.03，背景磁化率为0；给定地球正常磁场强度为56000nT，磁倾角为60°，磁偏角为0°，通过本发明方法计算了整个二维剖面上的磁梯度张量水平和垂向分量。
本发明方法利用Fortran语言编程实现，运行程序所用的个人电脑配置为：CPU-Inter Core i7-8700，主频为3.2GHz，运行内存为8.00GB。计算整个二维剖面40000个点所需时间约为27毫秒，可以看出本发明方法计算效率很高。图4、图5分别为解析式法和本发明方法计算的二维剖面的磁梯度梯度张量水平分量，从形态上看，两者是一致的，图6为水平梯度的解析解和本发明方法计算的相对误差，可以看出整个平面相对误差都小于1.92×10-4；图7为水平梯度地面测线的相对误差，可以看出地面相对误差总体小于4.73×10-6；图8、图9分别为解析式法和本发明计算的磁梯度张量垂向分量平面等值线图，从两幅图可以看出形态吻合很好，图10为垂向梯度的解析解和本发明方法计算的相对误差，可以看出整个平面相对误差都小于4.47×10-5；图11为垂向梯度地面测线的相对误差，可以看出地面相对误差总体小于2.66×10-6；表1给出了磁梯度张量各分量整个二维剖面的相对误差进行统计，由表1可知，本发明算法精度很高。
表1 本发明方法与解析式法相对误差统计
Figure BDA0003228312330000151
在一个实施例中，如图12所示，提供了一种二维磁梯度张量快速数值模拟装置，包括：网格剖分模块1202、磁化强度确定模块1204、核函数积分系数确定模块1206、矩阵扩展模块1208和磁梯度张量确定模块1210，其中：
网格剖分模块1202，用于根据二维异常体的大小信息，构建计算区域的二维模型，将计算区域沿x、z方向进行均匀网格剖分，得到多个网格单元，根据二维异常体的截面形状和磁化率分布，对网格单元的磁化率进行赋值；
磁化强度确定模块1204，用于根据已知的地球正常场，计算网格单元地球主磁场的x，z分量，根据网格单元地球主磁场的x，z分量和网格单元的磁化率，得到网格单元的磁化强度；
核函数积分系数确定模块1206，用于根据二维模型对预先构建的二维磁梯度张量公式进行离散，得到二维磁梯度的解析表达式，根据二维磁梯度的解析表达式，得到网格单元的核函数积分系数；
矩阵扩展模块1208，用于将网格单元的磁化强度和核函数积分系数表示为矩阵形式，并进行扩展，得到磁化强度扩展矩阵和核函数积分系数扩展矩阵；
磁梯度张量确定模块1210，用于通过二维离散傅里叶变换将核函数积分系数扩展矩阵和磁化强度扩展矩阵的二维卷积计算转换为乘积计算，再进行傅里叶反变换，得到计算区域的磁梯度张量。
磁化强度确定模块1204还用于根据已知的地球正常场，计算网格单元地球主磁场的x，z分量，根据网格单元地球主磁场的x，z分量和网格单元的磁化率，得到网格单元的磁化强度为：
Mx(xi,zj)＝χ(xi,zj)Tx(xi,zj)
Mz(xi,zj)＝χ(xi,zj)Tz(xi,zj)
其中，Tx(xi,zj)为网格单元地球主磁场的x分量，Tz(xi,zj)为网格单元地球主磁场的z分量，χ(xi,zj)表示编号为(xi,zj)单元的磁化率值，Mx(xi,zj)、Mz(xi,zj)分别表示(xi,zj)处磁化强度的x、z分量。
核函数积分系数确定模块1206还用于根据二维模型对预先构建的二维磁梯度张量公式进行离散，得到二维磁梯度的解析表达式为：
Figure BDA0003228312330000161
Figure BDA0003228312330000162
其中，Bxz表示磁梯度张量水平分量，Bzz表示磁梯度张量垂直分量，(x,z)表示观测点坐标，(x′,z′)表示源点坐标(xm,zn)表示编号为(m,n)单元的中心坐标，(x′i,z′j)表示编号为(i,j)单元的中心坐标；Mx(x′i,z′j)、Mz(x′i,z′j)表示编号为(i,j)单元的水平和垂向磁化强度；fx(xm-x′i,zn-z′j)，fz(xm-x′i,zn-z′j)表示观测点在(xm,zn)处针对编号为(i,j)单元的核函数积分系数。
核函数积分系数确定模块1206还用于根据二维磁梯度的解析表达式，得到网格单元的核函数积分系数为：
Figure BDA0003228312330000163
Figure BDA0003228312330000171
其中，μ0表示真空磁导率，Δx、Δz分别表示x方向和z方向单元网格间隔。
矩阵扩展模块1208还用于将网格单元的核函数积分系数表示为矩阵形式为：
Figure BDA0003228312330000172
将
Figure BDA0003228312330000173
扩展为
Figure BDA0003228312330000174
得到核函数积分系数扩展矩阵为：
Figure BDA0003228312330000175
将网格单元的磁化强度
Figure BDA0003228312330000176
扩展为
Figure BDA0003228312330000177
其中扩展的其他元素用零补齐，得到磁化强度扩展矩阵为：
Figure BDA0003228312330000178
磁梯度张量确定模块1210还用于将磁梯度张量表示为核函数积分系数和磁化强度的卷积：
Figure BDA0003228312330000181
将卷积
Figure BDA0003228312330000182
扩展为所述核函数积分系数扩展矩阵和所述磁化强度扩展矩阵的二维卷积：
Figure BDA0003228312330000183
通过二维离散傅里叶变换将核函数积分系数扩展矩阵和磁化强度扩展矩阵的二维卷积计算转换为乘积计算：
Figure BDA0003228312330000184
其中，F表示二维离散傅里叶正变换算子；
再进行傅里叶反变换，得到计算区域的磁梯度张量为：
Figure BDA0003228312330000185
其中，F-1表示二维离散傅里叶反变换算子；
Figure BDA0003228312330000186
表示提取矩阵的前Nz行，Nx列。
关于二维磁梯度张量快速数值模拟装置的具体限定可以参见上文中对于二维磁梯度张量快速数值模拟方法的限定，在此不再赘述。上述二维磁梯度张量快速数值模拟装置中的各个模块可全部或部分通过软件、硬件及其组合来实现。上述各模块可以硬件形式内嵌于或独立于计算机设备中的处理器中，也可以以软件形式存储于计算机设备中的存储器中，以便于处理器调用执行以上各个模块对应的操作。
在一个实施例中，提供了一种计算机设备，该计算机设备可以是终端，其内部结构图可以如图13所示。该计算机设备包括通过系统总线连接的处理器、存储器、网络接口、显示屏和输入装置。其中，该计算机设备的处理器用于提供计算和控制能力。该计算机设备的存储器包括非易失性存储介质、内存储器。该非易失性存储介质存储有操作系统和计算机程序。该内存储器为非易失性存储介质中的操作系统和计算机程序的运行提供环境。该计算机设备的网络接口用于与外部的终端通过网络连接通信。该计算机程序被处理器执行时以实现一种二维磁梯度张量快速数值模拟方法。该计算机设备的显示屏可以是液晶显示屏或者电子墨水显示屏，该计算机设备的输入装置可以是显示屏上覆盖的触摸层，也可以是计算机设备外壳上设置的按键、轨迹球或触控板，还可以是外接的键盘、触控板或鼠标等。
本领域技术人员可以理解，图13中示出的结构，仅仅是与本申请方案相关的部分结构的框图，并不构成对本申请方案所应用于其上的计算机设备的限定，具体的计算机设备可以包括比图中所示更多或更少的部件，或者组合某些部件，或者具有不同的部件布置。
在一个实施例中，提供了一种计算机设备，包括存储器和处理器，该存储器存储有计算机程序，该处理器执行计算机程序时实现上述方法实施例中的步骤。
在一个实施例中，提供了一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现上述方法实施例中的步骤。
本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的计算机程序可存储于一非易失性计算机可读取存储介质中，该计算机程序在执行时，可包括如上述各方法的实施例的流程。其中，本申请所提供的各实施例中所使用的对存储器、存储、数据库或其它介质的任何引用，均可包括非易失性和/或易失性存储器。非易失性存储器可包括只读存储器(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 (10)

1.一种二维磁梯度张量快速数值模拟方法，其特征在于，所述方法包括：
根据二维异常体的大小信息，构建计算区域的二维模型，将所述计算区域沿x、z方向进行均匀网格剖分，得到多个网格单元，根据所述二维异常体的截面形状和磁化率分布，对所述网格单元的磁化率进行赋值；
根据已知的地球正常场，计算所述网格单元地球主磁场的x，z分量，根据所述网格单元地球主磁场的x，z分量和所述网格单元的磁化率，得到所述网格单元的磁化强度；
根据所述二维模型对预先构建的二维磁梯度张量公式进行离散，得到二维磁梯度的解析表达式，根据所述二维磁梯度的解析表达式，得到网格单元的核函数积分系数；
将所述网格单元的磁化强度和所述核函数积分系数表示为矩阵形式，并进行扩展，得到磁化强度扩展矩阵和核函数积分系数扩展矩阵；
通过二维离散傅里叶变换将所述核函数积分系数扩展矩阵和所述磁化强度扩展矩阵的二维卷积计算转换为乘积计算，再进行傅里叶反变换，得到所述计算区域的磁梯度张量。
2.根据权利要求1所述的方法，其特征在于，根据已知的地球正常场，计算所述网格单元地球主磁场的x，z分量，根据所述网格单元地球主磁场的x，z分量和所述网格单元的磁化率，得到所述网格单元的磁化强度，包括：
根据已知的地球正常场，计算所述网格单元地球主磁场的x，z分量，根据所述网格单元地球主磁场的x，z分量和所述网格单元的磁化率，得到所述网格单元的磁化强度为：
Mx(xi,zj)＝χ(xi,zj)Tx(xi,zj)
Mz(xi,zj)＝χ(xi,zj)Tz(xi,zj)
其中，Tx(xi,zj)为所述网格单元地球主磁场的x分量，Tz(xi,zj)为所述网格单元地球主磁场的z分量，χ(xi,zj)表示编号为(xi,zj)单元的磁化率值，Mx(xi,zj)、Mz(xi,zj)分别表示(xi,zj)处磁化强度的x、z分量。
3.根据权利要求2所述的方法，其特征在于，根据所述二维模型对预先构建的二维磁梯度张量公式进行离散，得到二维磁梯度的解析表达式，包括：
根据所述二维模型对预先构建的二维磁梯度张量公式进行离散，得到二维磁梯度的解析表达式为：
Figure FDA0003228312320000021
Figure FDA0003228312320000022
其中，Bxz表示磁梯度张量水平分量，Bzz表示磁梯度张量垂直分量，(x,z)表示观测点坐标，(x′,z′)表示源点坐标(xm,zn)表示编号为(m,n)单元的中心坐标，(x′i,z′j)表示编号为(i,j)单元的中心坐标；Mx(x′i,z′j)、Mz(x′i,z′j)表示编号为(i,j)单元的水平和垂向磁化强度；fx(xm-x′i,zn-z′j)，fz(xm-x′i,zn-z′j)表示观测点在(xm,zn)处针对编号为(i,j)单元的核函数积分系数。
4.根据权利要求3所述的方法，其特征在于，根据所述二维磁梯度的解析表达式，得到网格单元的核函数积分系数，包括：
根据所述二维磁梯度的解析表达式，得到网格单元的核函数积分系数为：
Figure FDA0003228312320000023
Figure FDA0003228312320000031
其中，μ0表示真空磁导率，Δx、Δz分别表示x方向和z方向单元网格间隔。
5.根据权利要求4所述的方法，其特征在于，将所述网格单元的磁化强度和所述核函数积分系数表示为矩阵形式，并进行扩展，得到磁化强度扩展矩阵和核函数积分系数扩展矩阵，包括：
将所述网格单元的所述核函数积分系数表示为矩阵形式为：
Figure FDA0003228312320000032
将
Figure FDA0003228312320000033
扩展为
Figure FDA0003228312320000034
得到核函数积分系数扩展矩阵为：
Figure FDA0003228312320000035
将所述网格单元的磁化强度
Figure FDA0003228312320000036
扩展为
Figure FDA0003228312320000037
其中扩展的其他元素用零补齐，得到磁化强度扩展矩阵为：
Figure FDA0003228312320000041
6.根据权利要求5所述的方法，其特征在于，通过二维离散傅里叶变换将所述核函数积分系数扩展矩阵和所述磁化强度扩展矩阵的二维卷积计算转换为乘积计算，再进行傅里叶反变换，得到所述计算区域的磁梯度张量，包括：
将磁梯度张量表示为所述核函数积分系数和所述磁化强度的卷积：
Figure FDA0003228312320000042
将卷积
Figure FDA0003228312320000043
扩展为所述核函数积分系数扩展矩阵和所述磁化强度扩展矩阵的二维卷积：
Figure FDA0003228312320000044
通过二维离散傅里叶变换将所述核函数积分系数扩展矩阵和所述磁化强度扩展矩阵的二维卷积计算转换为乘积计算：
Figure FDA0003228312320000045
其中，F表示二维离散傅里叶正变换算子；
再进行傅里叶反变换，得到所述计算区域的磁梯度张量为：
Figure FDA0003228312320000046
其中，F-1表示二维离散傅里叶反变换算子；
Figure FDA0003228312320000047
表示提取矩阵的前Nz行，Nx列。
7.根据权利要求1至6任意一项所述的方法，其特征在于，不同网格单元的磁化率值不同。
8.一种二维磁梯度张量快速数值模拟装置，其特征在于，所述装置包括：
网格剖分模块，用于根据二维异常体的大小信息，构建计算区域的二维模型，将所述计算区域沿x、z方向进行均匀网格剖分，得到多个网格单元，根据所述二维异常体的截面形状和磁化率分布，对所述网格单元的磁化率进行赋值；
磁化强度确定模块，用于根据已知的地球正常场，计算所述网格单元地球主磁场的x，z分量，根据所述网格单元地球主磁场的x，z分量和所述网格单元的磁化率，得到所述网格单元的磁化强度；
核函数积分系数确定模块，用于根据所述二维模型对预先构建的二维磁梯度张量公式进行离散，得到二维磁梯度的解析表达式，根据所述二维磁梯度的解析表达式，得到网格单元的核函数积分系数；
矩阵扩展模块，用于将所述网格单元的磁化强度和所述核函数积分系数表示为矩阵形式，并进行扩展，得到磁化强度扩展矩阵和核函数积分系数扩展矩阵；
磁梯度张量确定模块，用于通过二维离散傅里叶变换将所述核函数积分系数扩展矩阵和所述磁化强度扩展矩阵的二维卷积计算转换为乘积计算，再进行傅里叶反变换，得到所述计算区域的磁梯度张量。
9.一种计算机设备，包括存储器和处理器，所述存储器存储有计算机程序，其特征在于，所述处理器执行所述计算机程序时实现权利要求1至7中任一项所述方法的步骤。
10.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现权利要求1至7中任一项所述的方法的步骤。
CN202110978629.3A 2021-08-25 2021-08-25 一种二维磁梯度张量快速数值模拟方法、装置和设备 Active CN113656976B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110978629.3A CN113656976B (zh) 2021-08-25 2021-08-25 一种二维磁梯度张量快速数值模拟方法、装置和设备

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110978629.3A CN113656976B (zh) 2021-08-25 2021-08-25 一种二维磁梯度张量快速数值模拟方法、装置和设备

Publications (2)

Publication Number Publication Date
CN113656976A true CN113656976A (zh) 2021-11-16
CN113656976B CN113656976B (zh) 2023-09-01

Family

ID=78481879

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110978629.3A Active CN113656976B (zh) 2021-08-25 2021-08-25 一种二维磁梯度张量快速数值模拟方法、装置和设备

Country Status (1)

Country Link
CN (1) CN113656976B (zh)

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20040113615A1 (en) * 2002-12-11 2004-06-17 The Board Of Trustees Of The Leland Stanford Junior University Characterization of the effect of spatial gradient field distortions in diffusion-weighted imaging
US20160356873A1 (en) * 2014-02-10 2016-12-08 Cr Development Ab Method for quantifying isotropic diffusion and/or anisotropic diffusion in a sample
CN112287534A (zh) * 2020-10-21 2021-01-29 中南大学 基于nufft的二维磁异常快速正演模拟方法和装置
CN112748471A (zh) * 2020-12-29 2021-05-04 中国地质大学（武汉） 一种非结构化等效源的重磁数据延拓与转换方法
CN113268702A (zh) * 2021-05-20 2021-08-17 中南大学 一种频率域磁梯度张量变换方法、装置和计算机设备

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20040113615A1 (en) * 2002-12-11 2004-06-17 The Board Of Trustees Of The Leland Stanford Junior University Characterization of the effect of spatial gradient field distortions in diffusion-weighted imaging
US20160356873A1 (en) * 2014-02-10 2016-12-08 Cr Development Ab Method for quantifying isotropic diffusion and/or anisotropic diffusion in a sample
CN112287534A (zh) * 2020-10-21 2021-01-29 中南大学 基于nufft的二维磁异常快速正演模拟方法和装置
CN112748471A (zh) * 2020-12-29 2021-05-04 中国地质大学（武汉） 一种非结构化等效源的重磁数据延拓与转换方法
CN113268702A (zh) * 2021-05-20 2021-08-17 中南大学 一种频率域磁梯度张量变换方法、装置和计算机设备

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
QINGZHU LI等: "Magnetic Object Positioning Based on Second-Order Magnetic Gradient Tensor System", pages 1 - 12, Retrieved from the Internet <URL:《网页在线公开：https://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=9469737》> *
舒晴等: "全张量磁梯度数据解释的均衡边界识别及深度成像技术", 《地球物理学报》, vol. 61, no. 4, pages 1539 - 1548 *

Also Published As

Publication number Publication date
CN113656976B (zh) 2023-09-01

Similar Documents

Publication Publication Date Title
Börner Numerical modelling in geo-electromagnetics: advances and challenges
CN112800657B (zh) 基于复杂地形的重力场数值模拟方法、装置和计算机设备
Spitzer A 3-D finite-difference algorithm for DC resistivity modelling using conjugate gradient methods
CN112287534B (zh) 基于nufft的二维磁异常快速正演模拟方法和装置
Fischer et al. On the C1 continuous discretization of non‐linear gradient elasticity: a comparison of NEM and FEM based on Bernstein–Bézier patches
CN113420487A (zh) 一种三维重力梯度张量数值模拟方法、装置、设备和介质
CN110346835B (zh) 大地电磁的正演方法、正演系统、存储介质及电子设备
CN103278848B (zh) 基于mpi并行预条件迭代的地震成像正演方法
Pardo et al. A survey on direct solvers for Galerkin methods
CN116842813B (zh) 一种三维大地电磁正演数值模拟方法
CN113076678B (zh) 一种频率域二度体重力异常快速数值模拟方法和装置
Mazzia et al. A comparison of numerical integration rules for the meshless local Petrov–Galerkin method
CN111967169A (zh) 二度体重力异常积分解数值模拟方法和装置
Zang et al. Isogeometric boundary element method for steady-state heat transfer with concentrated/surface heat sources
CN113673163B (zh) 一种三维磁异常数快速正演方法、装置和计算机设备
CN113642189B (zh) 一种基于积分解的重力梯度张量快速数值模拟方法和装置
Boroomand et al. Dynamic solution of unbounded domains using finite element method: discrete Green's functions in frequency domain
CN113656976B (zh) 一种二维磁梯度张量快速数值模拟方法、装置和设备
Xu et al. Solving fractional Laplacian visco-acoustic wave equations on complex-geometry domains using Grünwald-formula based radial basis collocation method
CN113221409A (zh) 一种有限元和边界元耦合的声波二维数值模拟方法和装置
CN113806686B (zh) 大规模复杂地质体重力梯度快速计算方法、装置和设备
Wan et al. Numerical solutions of incompressible Euler and Navier-Stokes equations by efficient discrete singular convolution method
CN114036806A (zh) 基于热导率各向异性介质的三维地温场数值模拟方法
Woodward et al. On the application of the method of difference potentials to linear elastic fracture mechanics
Zapata et al. A GPU parallel finite volume method for a 3D Poisson equation on arbitrary geometries

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
GR01 Patent grant
GR01 Patent grant