CN113656976B - 一种二维磁梯度张量快速数值模拟方法、装置和设备 - Google Patents
一种二维磁梯度张量快速数值模拟方法、装置和设备 Download PDFInfo
- Publication number
- CN113656976B CN113656976B CN202110978629.3A CN202110978629A CN113656976B CN 113656976 B CN113656976 B CN 113656976B CN 202110978629 A CN202110978629 A CN 202110978629A CN 113656976 B CN113656976 B CN 113656976B
- Authority
- CN
- China
- Prior art keywords
- dimensional
- grid
- magnetic gradient
- matrix
- gradient tensor
- 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.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 63
- 238000004088 simulation Methods 0.000 title claims abstract description 33
- 239000011159 matrix material Substances 0.000 claims abstract description 100
- 230000005415 magnetization Effects 0.000 claims abstract description 91
- 238000004364 calculation method Methods 0.000 claims abstract description 90
- 230000006870 function Effects 0.000 claims description 74
- 230000014509 gene expression Effects 0.000 claims description 33
- 230000002159 abnormal effect Effects 0.000 claims description 27
- 230000009466 transformation Effects 0.000 claims description 27
- 238000004590 computer program Methods 0.000 claims description 14
- 230000035699 permeability Effects 0.000 claims description 6
- 238000000605 extraction Methods 0.000 claims description 4
- 238000010276 construction Methods 0.000 abstract description 3
- 239000006185 dispersion Substances 0.000 abstract description 3
- 238000004458 analytical method Methods 0.000 description 7
- 238000010586 diagram Methods 0.000 description 4
- 230000005856 abnormality Effects 0.000 description 2
- 230000003044 adaptive effect Effects 0.000 description 2
- 230000001360 synchronised effect Effects 0.000 description 2
- 238000004422 calculation algorithm Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000010354 integration Effects 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 239000004973 liquid crystal related substance Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- 238000012876 topography Methods 0.000 description 1
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)
- Data Mining & Analysis (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- General Engineering & Computer Science (AREA)
- Databases & Information Systems (AREA)
- Algebra (AREA)
- Software Systems (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- Computing Systems (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分量。
在其中一个实施例中，还包括：根据所述二维模型对预先构建的二维磁梯度张量公式进行离散，得到二维磁梯度的解析表达式为：
其中，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)单元的核函数积分系数。
在其中一个实施例中，还包括：根据所述二维磁梯度的解析表达式，得到网格单元的核函数积分系数为：
其中，μ0表示真空磁导率，Δx、Δz分别表示x方向和z方向单元网格间隔。
在其中一个实施例中，还包括：将所述网格单元的所述核函数积分系数表示为矩阵形式为：
将扩展为得到核函数积分系数扩展矩阵为：
将所述网格单元的磁化强度扩展为其中扩展的其他元素用零补齐，得到磁化强度扩展矩阵为：
在其中一个实施例中，还包括：将磁梯度张量表示为所述核函数积分系数和所述磁化强度的卷积：
将卷积扩展为所述核函数积分系数扩展矩阵和所述磁化强度扩展矩阵的二维卷积：
通过二维离散傅里叶变换将所述核函数积分系数扩展矩阵和所述磁化强度扩展矩阵的二维卷积计算转换为乘积计算：
其中，F表示二维离散傅里叶正变换算子；
再进行傅里叶反变换，得到所述计算区域的磁梯度张量为：
其中，F-1表示二维离散傅里叶反变换算子；表示提取矩阵的前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分量。
在其中一个实施例中，还包括：对二维梯度张量公式进行离散，得
式中，Bxz表示磁梯度张量水平分量，Bzz表示磁梯度张量垂直分量，μ0表示真空磁导率，(x,z)表示观测点坐标，(x′,z′)表示源点坐标。
进一步地，得到二维磁梯度的解析表达式为：
其中，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)，根据二维磁梯度的解析表达式，得到网格单元的核函数积分系数为：
其中，μ0表示真空磁导率，Δx、Δz分别表示x方向和z方向单元网格间隔。
在其中一个实施例中，还包括：将网格单元的核函数积分系数表示为矩阵形式为：
将扩展为得到核函数积分系数扩展矩阵为：
将网格单元的磁化强度扩展为其中扩展的其他元素用零补齐，得到磁化强度扩展矩阵为：
在其中一个实施例中，还包括：将磁梯度张量表示为核函数积分系数和磁化强度的卷积：
将卷积扩展为所述核函数积分系数扩展矩阵和所述磁化强度扩展矩阵的二维卷积：
通过二维离散傅里叶变换将核函数积分系数扩展矩阵和磁化强度扩展矩阵的二维卷积计算转换为乘积计算：
其中，F表示二维离散傅里叶正变换算子；
再进行傅里叶反变换，得到计算区域的磁梯度张量为：
其中，F-1表示二维离散傅里叶反变换算子；表示提取矩阵的前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 本发明方法与解析式法相对误差统计
在一个实施例中，如图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还用于根据二维模型对预先构建的二维磁梯度张量公式进行离散，得到二维磁梯度的解析表达式为：
其中，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还用于根据二维磁梯度的解析表达式，得到网格单元的核函数积分系数为：
其中，μ0表示真空磁导率，Δx、Δz分别表示x方向和z方向单元网格间隔。
矩阵扩展模块1208还用于将网格单元的核函数积分系数表示为矩阵形式为：
将扩展为得到核函数积分系数扩展矩阵为：
将网格单元的磁化强度扩展为其中扩展的其他元素用零补齐，得到磁化强度扩展矩阵为：
磁梯度张量确定模块1210还用于将磁梯度张量表示为核函数积分系数和磁化强度的卷积：
将卷积扩展为所述核函数积分系数扩展矩阵和所述磁化强度扩展矩阵的二维卷积：
通过二维离散傅里叶变换将核函数积分系数扩展矩阵和磁化强度扩展矩阵的二维卷积计算转换为乘积计算：
其中，F表示二维离散傅里叶正变换算子；
再进行傅里叶反变换，得到计算区域的磁梯度张量为：
其中，F-1表示二维离散傅里叶反变换算子；表示提取矩阵的前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 (7)
1.一种二维磁梯度张量快速数值模拟方法，其特征在于，所述方法包括：
根据二维异常体的大小信息，构建计算区域的二维模型，将所述计算区域沿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分量；
根据所述二维模型对预先构建的二维磁梯度张量公式进行离散，得到二维磁梯度的解析表达式为：
其中，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)单元的核函数积分系数；
根据所述二维磁梯度的解析表达式，得到网格单元的核函数积分系数为：
其中，μ0表示真空磁导率，Δx、Δz分别表示x方向和z方向单元网格间隔；
将所述网格单元的磁化强度和所述核函数积分系数表示为矩阵形式，并进行扩展，得到磁化强度扩展矩阵和核函数积分系数扩展矩阵；
通过二维离散傅里叶变换将所述核函数积分系数扩展矩阵和所述磁化强度扩展矩阵的二维卷积计算转换为乘积计算，再进行傅里叶反变换，得到所述计算区域的磁梯度张量。
2.根据权利要求1所述的方法，其特征在于，将所述网格单元的磁化强度和所述核函数积分系数表示为矩阵形式，并进行扩展，得到磁化强度扩展矩阵和核函数积分系数扩展矩阵，包括：
将所述网格单元的所述核函数积分系数表示为矩阵形式为：
将扩展为得到核函数积分系数扩展矩阵为：
将所述网格单元的磁化强度扩展为其中扩展的其他元素用零补齐，得到磁化强度扩展矩阵为：
3.根据权利要求2所述的方法，其特征在于，通过二维离散傅里叶变换将所述核函数积分系数扩展矩阵和所述磁化强度扩展矩阵的二维卷积计算转换为乘积计算，再进行傅里叶反变换，得到所述计算区域的磁梯度张量，包括：
将磁梯度张量表示为所述核函数积分系数和所述磁化强度的卷积：
将卷积扩展为所述核函数积分系数扩展矩阵和所述磁化强度扩展矩阵的二维卷积：
通过二维离散傅里叶变换将所述核函数积分系数扩展矩阵和所述磁化强度扩展矩阵的二维卷积计算转换为乘积计算：
其中，F表示二维离散傅里叶正变换算子；
再进行傅里叶反变换，得到所述计算区域的磁梯度张量为：
其中，F-1表示二维离散傅里叶反变换算子；表示提取矩阵的前Nz行，Nx列。
4.根据权利要求1至3任意一项所述的方法，其特征在于，不同网格单元的磁化率值不同。
5.一种二维磁梯度张量快速数值模拟装置，其特征在于，所述装置包括：
网格剖分模块，用于根据二维异常体的大小信息，构建计算区域的二维模型，将所述计算区域沿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分量；
核函数积分系数确定模块，用于根据所述二维模型对预先构建的二维磁梯度张量公式进行离散，得到二维磁梯度的解析表达式为：
其中，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)单元的核函数积分系数；
根据所述二维磁梯度的解析表达式，得到网格单元的核函数积分系数为：
其中，μ0表示真空磁导率，Δx、Δz分别表示x方向和z方向单元网格间隔；
矩阵扩展模块，用于将所述网格单元的磁化强度和所述核函数积分系数表示为矩阵形式，并进行扩展，得到磁化强度扩展矩阵和核函数积分系数扩展矩阵；
磁梯度张量确定模块，用于通过二维离散傅里叶变换将所述核函数积分系数扩展矩阵和所述磁化强度扩展矩阵的二维卷积计算转换为乘积计算，再进行傅里叶反变换，得到所述计算区域的磁梯度张量。
6.一种计算机设备，包括存储器和处理器，所述存储器存储有计算机程序，其特征在于，所述处理器执行所述计算机程序时实现权利要求1至4中任一项所述方法的步骤。
7.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现权利要求1至4中任一项所述的方法的步骤。
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 CN113656976A (zh) | 2021-11-16 |
| CN113656976B true 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 (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| 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 | 中南大学 | 一种频率域磁梯度张量变换方法、装置和计算机设备 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6969991B2 (en) * | 2002-12-11 | 2005-11-29 | The Board Of Trustees Of The Leland Stanford Junior University | Correction of the effect of spatial gradient field distortions in diffusion-weighted imaging |
| JP6543639B2 (ja) * | 2014-02-10 | 2019-07-10 | シーアール ディベロップメント アーベー | サンプルにおける等方性拡散及び／又は異方性拡散を定量化するための方法 |
-
2021
- 2021-08-25 CN CN202110978629.3A patent/CN113656976B/zh active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| 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 (1)
| Title |
|---|
| 全张量磁梯度数据解释的均衡边界识别及深度成像技术;舒晴等;《地球物理学报》;第61卷(第4期);第1539-1548页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113656976A (zh) | 2021-11-16 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN112800657B (zh) | 基于复杂地形的重力场数值模拟方法、装置和计算机设备 | |
| Komatitsch et al. | Porting a high-order finite-element earthquake modeling application to NVIDIA graphics cards using CUDA | |
| 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 | |
| CN103278848B (zh) | 基于mpi并行预条件迭代的地震成像正演方法 | |
| CN113420487A (zh) | 一种三维重力梯度张量数值模拟方法、装置、设备和介质 | |
| CN116629079B (zh) | 混合有限元空间构造及求解线弹性力学问题的方法及装置 | |
| CN108021533A (zh) | 一种基于广义坐标系求解任意色散材料电磁特性的方法 | |
| CN113076678B (zh) | 一种频率域二度体重力异常快速数值模拟方法和装置 | |
| CN114021408A (zh) | 二维强磁场数值模拟方法、装置、设备及介质 | |
| Chen et al. | 3D gravity interface inversion constrained by a few points and its GPU acceleration | |
| CN113642189B (zh) | 一种基于积分解的重力梯度张量快速数值模拟方法和装置 | |
| CN113656976B (zh) | 一种二维磁梯度张量快速数值模拟方法、装置和设备 | |
| CN113673163B (zh) | 一种三维磁异常数快速正演方法、装置和计算机设备 | |
| Yamaguchi et al. | Fast crustal deformation computing method for multiple computations accelerated by a graphics processing unit cluster | |
| Xu et al. | Solving fractional Laplacian visco-acoustic wave equations on complex-geometry domains using Grünwald-formula based radial basis collocation method | |
| CN113806686B (zh) | 大规模复杂地质体重力梯度快速计算方法、装置和设备 | |
| Li et al. | A subspace expanding technique for global zero finding of multi-degree-of-freedom nonlinear systems | |
| Filippone et al. | CUDA dynamic active thread list strategy to accelerate debris flow simulations | |
| Zapata et al. | A GPU parallel finite volume method for a 3D Poisson equation on arbitrary geometries | |
| Rubio et al. | Optimizing isotropic and fully-anisotropic elastic modelling on multi-GPU platforms | |
| DAmbrosio et al. | Lava flow modeling by the sciara-fv3 parallel numerical code | |
| CN111274665A (zh) | 一种模拟裂纹尖端应力场的数值方法和系统 | |
| Farooq et al. | A pruning algorithm preserving modeling capabilities for polycrystalline data | |
| Agata et al. | Several hundred finite element analyses of an inversion of earthquake fault slip distribution using a high-fidelity model of the crustal structure |
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 |