CN113673163B - 一种三维磁异常数快速正演方法、装置和计算机设备 - Google Patents
一种三维磁异常数快速正演方法、装置和计算机设备 Download PDFInfo
- Publication number
- CN113673163B CN113673163B CN202110978731.3A CN202110978731A CN113673163B CN 113673163 B CN113673163 B CN 113673163B CN 202110978731 A CN202110978731 A CN 202110978731A CN 113673163 B CN113673163 B CN 113673163B
- Authority
- CN
- China
- Prior art keywords
- unit
- grid
- magnetic
- dimensional
- 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.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 63
- 230000002159 abnormal effect Effects 0.000 title claims abstract description 38
- 239000011159 matrix material Substances 0.000 claims abstract description 112
- 230000006870 function Effects 0.000 claims abstract description 94
- 230000005415 magnetization Effects 0.000 claims abstract description 68
- 230000008689 nuclear function Effects 0.000 claims abstract description 27
- 230000009466 transformation Effects 0.000 claims abstract description 20
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 18
- 210000004027 cell Anatomy 0.000 claims description 35
- 238000004364 calculation method Methods 0.000 claims description 26
- 238000004590 computer program Methods 0.000 claims description 14
- 230000010354 integration Effects 0.000 claims description 11
- 210000005056 cell body Anatomy 0.000 claims description 8
- 238000000605 extraction Methods 0.000 claims description 6
- 230000035699 permeability Effects 0.000 claims description 5
- 238000010586 diagram Methods 0.000 description 5
- 238000004088 simulation Methods 0.000 description 5
- 238000004458 analytical method Methods 0.000 description 4
- 238000005516 engineering process Methods 0.000 description 3
- 230000003044 adaptive effect Effects 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 230000005405 multipole Effects 0.000 description 2
- 230000001360 synchronised effect Effects 0.000 description 2
- 230000002547 anomalous effect Effects 0.000 description 1
- 239000006185 dispersion Substances 0.000 description 1
- 238000003384 imaging method 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
- 239000003208 petroleum Substances 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
- G06F30/27—Design optimisation, verification or simulation using machine learning, e.g. artificial intelligence, neural networks, support vector machines [SVM] or training a model
-
- 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
- G06F17/141—Discrete Fourier transforms
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/02—Reliability analysis or reliability optimisation; Failure analysis, e.g. worst case scenario performance, failure mode and effects analysis [FMEA]
-
- 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)
- Software Systems (AREA)
- Computational Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Pure & Applied Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- Evolutionary Computation (AREA)
- Mathematical Optimization (AREA)
- Mathematical Analysis (AREA)
- Artificial Intelligence (AREA)
- Medical Informatics (AREA)
- Discrete Mathematics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Geometry (AREA)
- Algebra (AREA)
- Computer Hardware Design (AREA)
- Databases & Information Systems (AREA)
- Measuring Magnetic Variables (AREA)
Abstract
本申请涉及一种三维磁异常数快速正演方法、装置、计算机设备和存储介质。所述方法包括：通过网格剖分形成多个长方体的网格单元；根据地球主磁场三分量和设置的单元体磁化率值得到单元体磁化强度；对预先构建的磁异常核函数三个分量的离散积分算子进行积分，得到核函数单元积分系数的解析表达式及核函数单元系数矩阵，其中核函数单元系数矩阵为Toeplitz矩阵；通过三维离散快速傅里叶变换进行核函数单元系数矩阵和单元体磁化强度三维离散卷积的快速相乘算法，得到待计算区域的磁异常值。本发明借助三维离散快速傅里叶变换实现核函数单元积分系数矩阵与磁化强度的快速相乘，一次可计算整个三维网格单元上的磁异常值。
Description
技术领域
本申请涉及计算机技术领域，特别是涉及一种三维磁异常数快速正演方法、装置、计算机设备和存储介质。
背景技术
随着计算机硬件技术发展及磁法观测数据量的急剧增加，近似的二维反演难以满足符合实际情况的真实三维地质构造解释的需求，因此磁法反演逐步发展到符合实际的三维。三维反演所需的内存和计算量也大幅度增加，因而对计算方法技术提出了更高的要求。正演是反演的基础，并且正演的精度和效率直接影响反演的效果和速度，因此实现高效、高精度的三维磁异常正演是实现三维反演的关键。
目前，针对三维地质体的磁异常的正演计算方法主要分为两种，一种是解析式法，另一种是数值模拟方法。解析式一般针对形状规则或者磁化率按一定规律分布的地质体能够推导出解析表达式，并且计算精度高；数值模拟方法则主要针对地下不均匀地质体。文献(骆遥,姚长利.长方体磁场及其梯度无解析奇点表达式理论研究.石油地球物理勘探,2007,04206):714-719.)推导了长方体磁场无奇异点的解析表达式，对于理论研究具有实际指导意义。文献(Li,K.,Chen,L.W.,Chen,Q.R.,et al.Fast 3D forward modeling ofthe magnetic field and gradient tensor on an undulated surface.AppliedGeophysics,2018,15(3):500-512.)提出了一种基于复杂地形的波数域磁异常及梯度张量快速数值模拟方法，其中采用高斯傅里叶变换克服常规快速傅里叶变换的边界效应的影响，提高了频率域方法的计算精度，但随着高斯点个数的增加计算量也会成倍增加。(肖晓,黄保尚,任政勇,等.基于自适应多层快速多极算法的大规模磁法正演模拟.地球物理学报,2019,62(03):236-246.)利用非结构化四面体网格剖分方式实现了自适应多层快速多极的磁法正演计算，由于采用了非结构化网格，因此对于地形和复杂模型都能很好的拟合，但对于大规模磁异常快速计算，速度较慢。
已有的磁异常数值模拟方法难以满足大规模、精细化三维磁异常快速成像的需求。因此，亟需研究一种高效、高精度的三维磁异常数值模拟方法，以解决上述问题。
发明内容
基于此，有必要针对上述技术问题，提供一种能够提高磁异常正演效率和精度的三维磁异常数快速正演方法、装置、计算机设备和存储介质。
一种三维磁异常数快速正演方法，所述方法包括：
获取待计算区域在三维坐标系中的位置信息，将所述待计算区域沿x、y、z方向分别进行等间隔剖分，形成多个长方体的网格单元；
根据地球主磁场模型，计算每个所述网格单元中心点处的地球主磁场三分量，根据待计算区域中三维地质体的磁化率分布设定所述网格单元的磁化率，得到单元体磁化率值；根据所述地球主磁场三分量和所述单元体磁化率值得到每个所述网格单元中心点处的单元体磁化强度；
对预先构建的磁异常核函数三个分量的离散积分算子进行积分，得到核函数单元积分系数的解析表达式，根据所述核函数单元积分系数的解析表达式得到核函数单元系数矩阵；所述核函数单元系数矩阵为Toeplitz矩阵；
通过三维离散快速傅里叶变换进行所述核函数单元系数矩阵和所述单元体磁化强度三维离散卷积的快速相乘算法，得到所述待计算区域的磁异常值。
在其中一个实施例中，还包括：构建三维坐标系，其中x轴指向正东，y轴指向正北，z轴垂直向下；
获取待计算区域在所述三维坐标系中的位置信息；
将所述待计算区域沿x、y、z方向分别等间隔剖分成Nx,Ny,Nz个小长方体，网格单元间隔分别为Δx、Δy、Δz，形成多个长方体的网格单元。
在其中一个实施例中，还包括：根据地球主磁场模型，计算每个所述网格单元中心点处的地球主磁场三分量为：
其中，表示编号为(l,m,n)单元的中心坐标，其中，l＝0,1,…Nx-1；m＝0,1,…Ny-1；n＝0,1,…Nz-1，α表示当地磁倾角，β表示当地磁偏角，Tx、Ty、Tz分别表示地球主磁场T的x、y、z分量；
根据待计算区域中三维地质体的磁化率分布设定所述网格单元的磁化率，得到单元体磁化率值；
根据所述地球主磁场三分量和所述单元体磁化率值得到每个所述网格单元中心点处的单元体磁化强度为：
其中，Mx、My、Mz分别表示单元体磁化强度的x、y、z分量。
在其中一个实施例中，还包括：以所述网格单元的中心点位置为观测点；
当计算磁异常垂向分量时，获取预先构建的磁异常垂向分量核函数三个分量的离散积分算子为：
其中，μ0表示真空磁导率，(x,y,z)为观测点坐标，表示源点坐标；
(xi,yj,zk)表示编号为(i,j,k)的网格单元的中心坐标，表示编号为(l,m,n)的网格单元的中心坐标,其中，i,l＝0,1,…Nx-1；j,m＝0,1,…Ny-1；k,n＝0,1,…Nz-1，编号为(l,m,n)的单元三个方向积分区间分别为/>至/> 至 至/>
对所述磁异常垂向分量核函数三个分量的离散积分算子进行积分，得到核函数垂向分量三个分量的单元积分系数的解析表达式为：
其中，p,q,r＝1,2为整数，cpqr＝(-1)p(-1)q(-1)r，
在其中一个实施例中，还包括：根据所述核函数单元积分系数的解析表达式得到核函数单元系数矩阵；所述核函数单元系数矩阵为Toeplitz矩阵；所述核函数单元系数矩阵的大小为
在其中一个实施例中，还包括：通过三维离散快速傅里叶变换进行所述核函数单元系数矩阵和所述单元体磁化强度三维离散卷积的快速相乘算法，得到所述待计算区域的磁异常值为：
其中，F和F-1分别表示三维离散傅里叶正反变换算子；表示磁化强度扩展矩阵，M表示网格剖分单元体的磁化强度矩阵，/>为扩展矩阵，/>表示提取矩阵的前Nx×Ny×Nz个元素，即整个三维网格单元上的磁异常值。
在其中一个实施例中，还包括：不同网格单元的单元体磁化率值不同。
一种三维磁异常数快速正演装置，所述装置包括：
网格剖分模块，用于获取待计算区域在三维坐标系中的位置信息，将所述待计算区域沿x、y、z方向分别进行等间隔剖分，形成多个长方体的网格单元；
单元体磁化强度确定模块，用于根据地球主磁场模型，计算每个所述网格单元中心点处的地球主磁场三分量，根据待计算区域中三维地质体的磁化率分布设定所述网格单元的磁化率，得到单元体磁化率值；根据所述地球主磁场三分量和所述单元体磁化率值得到每个所述网格单元中心点处的单元体磁化强度；
核函数单元系数矩阵确定模块，用于对预先构建的磁异常核函数三个分量的离散积分算子进行积分，得到核函数单元积分系数的解析表达式，根据所述核函数单元积分系数的解析表达式得到核函数单元系数矩阵；所述核函数单元系数矩阵为Toeplitz矩阵；
磁异常值计算模块，用于通过三维离散快速傅里叶变换进行所述核函数单元系数矩阵和所述单元体磁化强度三维离散卷积的快速相乘算法，得到所述待计算区域的磁异常值。
一种计算机设备，包括存储器和处理器，所述存储器存储有计算机程序，所述处理器执行所述计算机程序时实现以下步骤：
获取待计算区域在三维坐标系中的位置信息，将所述待计算区域沿x、y、z方向分别进行等间隔剖分，形成多个长方体的网格单元；
根据地球主磁场模型，计算每个所述网格单元中心点处的地球主磁场三分量，根据待计算区域中三维地质体的磁化率分布设定所述网格单元的磁化率，得到单元体磁化率值；根据所述地球主磁场三分量和所述单元体磁化率值得到每个所述网格单元中心点处的单元体磁化强度；
对预先构建的磁异常核函数三个分量的离散积分算子进行积分，得到核函数单元积分系数的解析表达式，根据所述核函数单元积分系数的解析表达式得到核函数单元系数矩阵；所述核函数单元系数矩阵为Toeplitz矩阵；
通过三维离散快速傅里叶变换进行所述核函数单元系数矩阵和所述单元体磁化强度三维离散卷积的快速相乘算法，得到所述待计算区域的磁异常值。
一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现以下步骤：
获取待计算区域在三维坐标系中的位置信息，将所述待计算区域沿x、y、z方向分别进行等间隔剖分，形成多个长方体的网格单元；
根据地球主磁场模型，计算每个所述网格单元中心点处的地球主磁场三分量，根据待计算区域中三维地质体的磁化率分布设定所述网格单元的磁化率，得到单元体磁化率值；根据所述地球主磁场三分量和所述单元体磁化率值得到每个所述网格单元中心点处的单元体磁化强度；
对预先构建的磁异常核函数三个分量的离散积分算子进行积分，得到核函数单元积分系数的解析表达式，根据所述核函数单元积分系数的解析表达式得到核函数单元系数矩阵；所述核函数单元系数矩阵为Toeplitz矩阵；
通过三维离散快速傅里叶变换进行所述核函数单元系数矩阵和所述单元体磁化强度三维离散卷积的快速相乘算法，得到所述待计算区域的磁异常值。
上述三维磁异常数快速正演方法、装置、计算机设备和存储介质，通过获取待计算区域在三维坐标系中的位置信息，将待计算区域沿x、y、z方向分别进行等间隔剖分，形成多个长方体的网格单元；根据地球主磁场模型，计算每个网格单元中心点处的地球主磁场三分量，根据待计算区域中三维地质体的磁化率分布设定网格单元的磁化率，得到单元体磁化率值；根据地球主磁场三分量和单元体磁化率值得到每个网格单元中心点处的单元体磁化强度；对预先构建的磁异常核函数三个分量的离散积分算子进行积分，得到核函数单元积分系数的解析表达式，根据核函数单元积分系数的解析表达式得到核函数单元系数矩阵，其中核函数单元系数矩阵为Toeplitz矩阵；通过三维离散快速傅里叶变换进行核函数单元系数矩阵和单元体磁化强度三维离散卷积的快速相乘算法，得到待计算区域的磁异常值。本发明每个网格单元的核函数积分都可推导出解析解，提高了计算精度，充分利用Toeplitz型矩阵的特点，大大降低核函数矩阵的存储空间和计算时间；同时，借助三维离散快速傅里叶变换实现核函数单元积分系数矩阵与磁化强度的快速相乘，一次可计算整个三维网格单元上的磁异常值，有利于进行井地联合反演研究。
附图说明
图1为一个实施例中三维磁异常数快速正演方法的流程示意图；
图2为一个实施例中计算区域网格剖分示意图；
图3为另一个实施例中三维磁异常数快速正演方法的流程示意图；
图4为一个具体实施例中三维磁异常数快速正演方法的流程示意图；
图5为一个具体实施例中本发明方法计算的地面磁异常结果图；
图6为一个具体实施例中解析解计算的地面磁异常结果图；
图7为一个具体实施例中解析解与本发明方法计算的相对误差图；
图8为一个具体实施例中本发明方法计算的穿过异常体的磁异常结果图；
图9为一个具体实施例中解析解计算的穿过异常体的磁异常图；
图10为一个实施例中三维磁异常数快速正演装置的结构框图；
图11为一个实施例中计算机设备的内部结构图。
具体实施方式
为了使本申请的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本申请进行进一步详细说明。应当理解，此处描述的具体实施例仅仅用以解释本申请，并不用于限定本申请。
本申请提供的三维磁异常数快速正演方法，可以应用于如下应用环境中。其中，终端执行一种三维磁异常数快速正演方法。通过获取待计算区域在三维坐标系中的位置信息，将待计算区域沿x、y、z方向分别进行等间隔剖分，形成多个长方体的网格单元；根据地球主磁场模型，计算每个网格单元中心点处的地球主磁场三分量，根据待计算区域中三维地质体的磁化率分布设定网格单元的磁化率，得到单元体磁化率值；根据地球主磁场三分量和单元体磁化率值得到每个网格单元中心点处的单元体磁化强度；对预先构建的磁异常核函数三个分量的离散积分算子进行积分，得到核函数单元积分系数的解析表达式，根据核函数单元积分系数的解析表达式得到核函数单元系数矩阵，其中核函数单元系数矩阵为Toeplitz矩阵；通过三维离散快速傅里叶变换进行核函数单元系数矩阵和单元体磁化强度三维离散卷积的快速相乘算法，得到待计算区域的磁异常值。其中，终端可以但不限于是各种个人计算机、笔记本电脑、智能手机和便携式设备。
在一个实施例中，如图1所示，提供了一种三维磁异常数快速正演方法，包括以下步骤：
步骤102，获取待计算区域在三维坐标系中的位置信息，将待计算区域沿x、y、z方向分别进行等间隔剖分，形成多个长方体的网格单元。
如图2所示，将该计算区域沿着x、y、z方向分别等间隔剖分成Nx,Ny,Nz个小长方体，网格单元间隔分别为Δx、Δy、Δz。观测点的位置和每个单元的中心点位置重合，观测点的总数为Nx×Ny×Nz个。
步骤104，根据地球主磁场模型，计算每个网格单元中心点处的地球主磁场三分量，根据待计算区域中三维地质体的磁化率分布设定网格单元的磁化率，得到单元体磁化率值；根据地球主磁场三分量和单元体磁化率值得到每个网格单元中心点处的单元体磁化强度。
根据三维地质体的磁化率分布设定网格单元的磁化率值，每个单元的磁化率是一个常数。
步骤106，对预先构建的磁异常核函数三个分量的离散积分算子进行积分，得到核函数单元积分系数的解析表达式，根据核函数单元积分系数的解析表达式得到核函数单元系数矩阵。
核函数单元系数矩阵为Toeplitz矩阵，Toeplitz矩阵的主对角线上的元素相等，平行于主对角线的线上的元素也相等。
步骤108，用于通过三维离散快速傅里叶变换进行核函数单元系数矩阵和单元体磁化强度三维离散卷积的快速相乘算法，得到待计算区域的磁异常值。
地下地质体在任意一点产生的磁异常可以表示为核函数与磁化强度的三维卷积，磁异常垂向分量可表示为：
式中，h表示磁异常的核函数单元积分系数矩阵，表示三维卷积算子，M表示网格单元的磁化强度，Bz表示磁异常垂向分量。
核函数单元系数矩阵和单元体磁化强度三维离散卷积即可得到所求磁异常，这里以一维离散卷积为例详述核函数系数矩阵与磁化强度向量快速相乘算法。
针对一维离散卷积可以表示为
则核函数单元积分系数矩阵h可以表示为矩阵的形式
上式可以进一步表示为
则上述Toeplitz矩阵h只需计算和存储第一行和第一列h＝Toeplitz(t)，进一步将t写成则/>则Toeplitz矩阵与向量的快速相乘可以表示为
式中，F和F-1分别表示一维离散傅里叶正反变换算子；一维磁化强度向量扩展为 表示提取矩阵的前Nx个元素。
三维离散卷积核函数单元积分系数矩阵与磁化强度快速相乘与一维离散卷积类似。通过三维离散快速傅里叶变换进行核函数单元系数矩阵和单元体磁化强度三维离散卷积的快速相乘算法可以一次性提取整个三维网格单元上的磁异常值。
上述三维磁异常数快速正演方法中，通过获取待计算区域在三维坐标系中的位置信息，将待计算区域沿x、y、z方向分别进行等间隔剖分，形成多个长方体的网格单元；根据地球主磁场模型，计算每个网格单元中心点处的地球主磁场三分量，根据待计算区域中三维地质体的磁化率分布设定网格单元的磁化率，得到单元体磁化率值；根据地球主磁场三分量和单元体磁化率值得到每个网格单元中心点处的单元体磁化强度；对预先构建的磁异常核函数三个分量的离散积分算子进行积分，得到核函数单元积分系数的解析表达式，根据核函数单元积分系数的解析表达式得到核函数单元系数矩阵，其中核函数单元系数矩阵为Toeplitz矩阵；通过三维离散快速傅里叶变换进行核函数单元系数矩阵和单元体磁化强度三维离散卷积的快速相乘算法，得到待计算区域的磁异常值。本发明每个网格单元的核函数积分都可推导出解析解，提高了计算精度，充分利用Toeplitz型矩阵的特点，大大降低核函数矩阵的存储空间和计算时间；同时，借助三维离散快速傅里叶变换实现核函数单元积分系数矩阵与磁化强度的快速相乘，一次可计算整个三维网格单元上的磁异常值，有利于进行井地联合反演研究。
在其中一个实施例中，还包括：构建三维坐标系，其中x轴指向正东，y轴指向正北，z轴垂直向下；获取待计算区域在三维坐标系中的位置信息；将待计算区域沿x、y、z方向分别等间隔剖分成Nx,Ny,Nz个小长方体，网格单元间隔分别为Δx、Δy、Δz，形成多个长方体的网格单元。
在其中一个实施例中，还包括：根据地球主磁场模型，计算每个网格单元中心点处的地球主磁场三分量为：
其中，表示编号为(l,m,n)单元的中心坐标，其中，l＝0,1,…Nx-1；m＝0,1,…Ny-1；n＝0,1,…Nz-1，α表示当地磁倾角，β表示当地磁偏角，Tx、Ty、Tz分别表示地球主磁场T的x、y、z分量；
根据待计算区域中三维地质体的磁化率分布设定网格单元的磁化率，得到单元体磁化率值；
根据地球主磁场三分量和单元体磁化率值得到每个网格单元中心点处的单元体磁化强度为：
其中，Mx、My、Mz分别表示单元体磁化强度的x、y、z分量。
在其中一个实施例中，还包括：以网格单元的中心点位置为观测点；
当计算磁异常垂向分量时，获取预先构建的磁异常垂向分量核函数三个分量的离散积分算子为：
其中，μ0表示真空磁导率，(x,y,z)为观测点坐标，表示源点坐标；
(xi,yj,zk)表示编号为(i,j,k)的网格单元的中心坐标，表示编号为(l,m,n)的网格单元的中心坐标,其中，i,l＝0,1,…Nx-1；j,m＝0,1,…Ny-1；k,n＝0,1,…Nz-1，编号为(l,m,n)的单元三个方向积分区间分别为/>至/> 至 至/>
对磁异常垂向分量核函数三个分量的离散积分算子进行积分，得到核函数垂向分量三个分量的单元积分系数的解析表达式为：
其中，p,q,r＝1,2为整数，cpqr＝(-1)p(-1)q(-1)r，
在其中一个实施例中，还包括：根据核函数单元积分系数的解析表达式得到核函数单元系数矩阵；核函数单元系数矩阵为Toeplitz矩阵；核函数单元系数矩阵的大小为
在其中一个实施例中，还包括：通过三维离散快速傅里叶变换进行核函数单元系数矩阵和单元体磁化强度三维离散卷积的快速相乘算法，得到待计算区域的磁异常值为：
其中，F和F-1分别表示三维离散傅里叶正反变换算子；表示磁化强度扩展矩阵，M表示网格剖分单元体的磁化强度矩阵，/>为扩展矩阵，/>表示提取矩阵的前Nx×Ny×Nz个元素，即整个三维网格单元上的磁异常值。
在其中一个实施例中，还包括：不同网格单元的单元体磁化率值不同。
应该理解的是，虽然图1的流程图中的各个步骤按照箭头的指示依次显示，但是这些步骤并不是必然按照箭头指示的顺序依次执行。除非本文中有明确的说明，这些步骤的执行并没有严格的顺序限制，这些步骤可以以其它的顺序执行。而且，图1中的至少一部分步骤可以包括多个子步骤或者多个阶段，这些子步骤或者阶段并不必然是在同一时刻执行完成，而是可以在不同的时刻执行，这些子步骤或者阶段的执行顺序也不必然是依次进行，而是可以与其它步骤或者其它步骤的子步骤或者阶段的至少一部分轮流或者交替地执行。
在另一个实施例中，如图3所示，提供了一种三维磁异常数快速正演方法，包括：
步骤一：计算区域离散：根据异常体的大小对计算区域沿x、y、z方向分别等间隔剖分；
步骤二：设置地球主磁场和模型磁化率：根据公知的地球主磁场模型，计算上述剖分网格单元的主磁场三分量，根据三维地质体磁化率分布设定网格单元的磁化率值；
步骤三：计算磁化强度三分量：根据设定的网格单元主磁场三分量及磁化率，计算磁化强度三分量；
步骤四：给出磁异常满足的三维卷积表达式：地质体产生的磁异常可以表示为核函数与磁化强度的三维卷积；
步骤五：计算核函数积分系数矩阵：对剖分的网格单元积分，得到磁异常核函数单元积分系数矩阵；
步骤六：计算三维网格单元磁异常，采用三维离散傅里叶变换实现核函数矩阵与磁化强度的快速卷积计算。
在其中一个具体实施例中，计算区域内有一个棱柱型异常体，计算区域范围为：x和y方向均从-2000m到2000m，z方向从0m到1000m(z轴垂直向下为正)，x和y方向网格剖分间隔均为40m，z方向网格剖分间隔均10m，整个计算区域剖分为100×100×100个单元，异常体分布范围为：x和y方向均从-400m到400m，z方向从250m到450m，磁化率为0.01(SI)；；研究区域地球主磁场为45000nT，磁倾角为90°，磁偏角为0°，计算整个三维网格上的磁异常。
本发明方法利用Fortran语言编程实现，运行程序所用的个人电脑配置为：CPU-Inter Core i7-8700，主频为3.2GHz，运行内存为8.00GB。计算1×106个观测点的值需要1.18秒。图4为本发明方法计算的地表磁异常垂向分量结果图；图5为采用解析解的方法计算的地表磁异常垂向分量结果图；从两幅图可以看出解析解与本发明方法计算结果形态一致；图6为磁异常的解析解和本发明方法计算的相对误差，可以看出整个平面上相对误差都小于1.1×10-10，进一步可以看出该方法计算精度较高；图7、图8分别为本发明方法和解析解计算的穿过异常体内部(z＝395m)的磁异常，可以看出两幅图形态吻合很好，图9为穿过异常的磁异常异常相对误差等值线图，可以看出在异常体内部，本方法计算精度依然很高。
在一个实施例中，如图10所示，提供了一种三维磁异常数快速正演装置，包括：网格剖分模块1002、单元体磁化强度确定模块1004、核函数单元系数矩阵确定模块1006和磁异常值计算模块1008，其中：
网格剖分模块1002，用于获取待计算区域在三维坐标系中的位置信息，将待计算区域沿x、y、z方向分别进行等间隔剖分，形成多个长方体的网格单元；
单元体磁化强度确定模块1004，用于根据地球主磁场模型，计算每个网格单元中心点处的地球主磁场三分量，根据待计算区域中三维地质体的磁化率分布设定网格单元的磁化率，得到单元体磁化率值；根据地球主磁场三分量和单元体磁化率值得到每个网格单元中心点处的单元体磁化强度；
核函数单元系数矩阵确定模块1006，用于对预先构建的磁异常核函数三个分量的离散积分算子进行积分，得到核函数单元积分系数的解析表达式，根据核函数单元积分系数的解析表达式得到核函数单元系数矩阵；核函数单元系数矩阵为Toeplitz矩阵；
磁异常值计算模块1008，用于通过三维离散快速傅里叶变换进行核函数单元系数矩阵和单元体磁化强度三维离散卷积的快速相乘算法，得到待计算区域的磁异常值。
网格剖分模块1002还用于构建三维坐标系，其中x轴指向正东，y轴指向正北，z轴垂直向下；获取待计算区域在三维坐标系中的位置信息；将待计算区域沿x、y、z方向分别等间隔剖分成Nx,Ny,Nz个小长方体，网格单元间隔分别为Δx、Δy、Δz，形成多个长方体的网格单元。
单元体磁化强度确定模块1004还用于根据地球主磁场模型，计算每个网格单元中心点处的地球主磁场三分量为：
/>
其中，表示编号为(l,m,n)单元的中心坐标，其中，l＝0,1,…Nx-1；m＝0,1,…Ny-1；n＝0,1,…Nz-1，α表示当地磁倾角，β表示当地磁偏角，Tx、Ty、Tz分别表示地球主磁场T的x、y、z分量；
根据待计算区域中三维地质体的磁化率分布设定网格单元的磁化率，得到单元体磁化率值；
根据地球主磁场三分量和单元体磁化率值得到每个网格单元中心点处的单元体磁化强度为：
其中，Mx、My、Mz分别表示单元体磁化强度的x、y、z分量。
核函数单元系数矩阵确定模块1006还用于以网格单元的中心点位置为观测点；
当计算磁异常垂向分量时，获取预先构建的磁异常垂向分量核函数三个分量的离散积分算子为：
其中，μ0表示真空磁导率，(x,y,z)为观测点坐标，表示源点坐标；
(xi,yj,zk)表示编号为(i,j,k)的网格单元的中心坐标，表示编号为(l,m,n)的网格单元的中心坐标,其中，i,l＝0,1,…Nx-1；j,m＝0,1,…Ny-1；k,n＝0,1,…Nz-1，编号为(l,m,n)的单元三个方向积分区间分别为/>至/> 至 至/>
对磁异常垂向分量核函数三个分量的离散积分算子进行积分，得到核函数垂向分量三个分量的单元积分系数的解析表达式为：
/>
其中，p,q,r＝1,2为整数，cpqr＝(-1)p(-1)q(-1)r，
磁异常值计算模块1008还用于通过三维离散快速傅里叶变换进行核函数单元系数矩阵和单元体磁化强度三维离散卷积的快速相乘算法，得到待计算区域的磁异常值为：
其中，F和F-1分别表示三维离散傅里叶正反变换算子；表示磁化强度扩展矩阵，M表示网格剖分单元体的磁化强度矩阵，/>为扩展矩阵，/>表示提取矩阵的前Nx×Ny×Nz个元素，即整个三维网格单元上的磁异常值。
关于三维磁异常数快速正演装置的具体限定可以参见上文中对于三维磁异常数快速正演方法的限定，在此不再赘述。上述三维磁异常数快速正演装置中的各个模块可全部或部分通过软件、硬件及其组合来实现。上述各模块可以硬件形式内嵌于或独立于计算机设备中的处理器中，也可以以软件形式存储于计算机设备中的存储器中，以便于处理器调用执行以上各个模块对应的操作。
在一个实施例中，提供了一种计算机设备，该计算机设备可以是终端，其内部结构图可以如图11所示。该计算机设备包括通过系统总线连接的处理器、存储器、网络接口、显示屏和输入装置。其中，该计算机设备的处理器用于提供计算和控制能力。该计算机设备的存储器包括非易失性存储介质、内存储器。该非易失性存储介质存储有操作系统和计算机程序。该内存储器为非易失性存储介质中的操作系统和计算机程序的运行提供环境。该计算机设备的网络接口用于与外部的终端通过网络连接通信。该计算机程序被处理器执行时以实现一种三维磁异常数快速正演方法。该计算机设备的显示屏可以是液晶显示屏或者电子墨水显示屏，该计算机设备的输入装置可以是显示屏上覆盖的触摸层，也可以是计算机设备外壳上设置的按键、轨迹球或触控板，还可以是外接的键盘、触控板或鼠标等。
本领域技术人员可以理解，图11中示出的结构，仅仅是与本申请方案相关的部分结构的框图，并不构成对本申请方案所应用于其上的计算机设备的限定，具体的计算机设备可以包括比图中所示更多或更少的部件，或者组合某些部件，或者具有不同的部件布置。
在一个实施例中，提供了一种计算机设备，包括存储器和处理器，该存储器存储有计算机程序，该处理器执行计算机程序时实现上述方法实施例中的步骤。
在一个实施例中，提供了一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现上述方法实施例中的步骤。
本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的计算机程序可存储于一非易失性计算机可读取存储介质中，该计算机程序在执行时，可包括如上述各方法的实施例的流程。其中，本申请所提供的各实施例中所使用的对存储器、存储、数据库或其它介质的任何引用，均可包括非易失性和/或易失性存储器。非易失性存储器可包括只读存储器(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 (8)
1.一种三维磁异常数快速正演方法，其特征在于，所述方法包括：
获取待计算区域在三维坐标系中的位置信息，将所述待计算区域沿x、y、z方向分别进行等间隔剖分，形成多个长方体的网格单元；
根据地球主磁场模型，计算每个所述网格单元中心点处的地球主磁场三分量，根据待计算区域中三维地质体的磁化率分布设定所述网格单元的磁化率，得到单元体磁化率值；根据所述地球主磁场三分量和所述单元体磁化率值得到每个所述网格单元中心点处的单元体磁化强度；
以所述网格单元的中心点位置为观测点，当计算磁异常垂向分量时，获取预先构建的磁异常垂向分量核函数三个分量的离散积分算子为：
其中，μ0表示真空磁导率，(x,y,z)为观测点坐标，表示源点坐标；
(xi,yj,zk)表示编号为(i,j,k)的网格单元的中心坐标，表示编号为(l,m,n)的网格单元的中心坐标,其中，i,l＝0,1,…Nx-1；j,m＝0,1,…Ny-1；k,n＝0,1,…Nz-1，编号为(l,m,n)的单元三个方向积分区间分别为/>至/>至至/>
对所述磁异常垂向分量核函数三个分量的离散积分算子进行积分，得到核函数垂向分量三个分量的单元积分系数的解析表达式为：
其中，p,q,r＝1,2为整数，cpqr＝(-1)p(-1)q(-1)r，
根据所述核函数单元积分系数的解析表达式得到核函数单元系数矩阵；所述核函数单元系数矩阵为Toeplitz矩阵；
通过三维离散快速傅里叶变换进行所述核函数单元系数矩阵和所述单元体磁化强度三维离散卷积的快速相乘算法，得到所述待计算区域的磁异常值为：
其中，F和F-1分别表示三维离散傅里叶正反变换算子；表示磁化强度扩展矩阵，M表示网格剖分单元体的磁化强度矩阵，/>为扩展矩阵，/>表示提取矩阵的前Nx×Ny×Nz个元素，即整个三维网格单元上的磁异常值。
2.根据权利要求1所述的方法，其特征在于，获取待计算区域在三维坐标系中的位置信息，将所述待计算区域沿x、y、z方向分别进行等间隔剖分，形成多个长方体的网格单元，包括：
构建三维坐标系，其中x轴指向正东，y轴指向正北，z轴垂直向下；
获取待计算区域在所述三维坐标系中的位置信息；
将所述待计算区域沿x、y、z方向分别等间隔剖分成Nx,Ny,Nz个小长方体，网格单元间隔分别为Δx、Δy、Δz，形成多个长方体的网格单元。
3.根据权利要求2所述的方法，其特征在于，根据地球主磁场模型，计算每个所述网格单元中心点处的地球主磁场三分量，根据待计算区域中三维地质体的磁化率分布设定所述网格单元的磁化率，得到单元体磁化率值；根据所述地球主磁场三分量和所述单元体磁化率值得到每个所述网格单元中心点处的单元体磁化强度，包括：
根据地球主磁场模型，计算每个所述网格单元中心点处的地球主磁场三分量为：
其中，表示编号为(l,m,n)单元的中心坐标，其中，l＝0,1,…Nx-1；m＝0,1,…Ny-1；n＝0,1,…Nz-1，α表示当地磁倾角，β表示当地磁偏角，Tx、Ty、Tz分别表示地球主磁场T的x、y、z分量；
根据待计算区域中三维地质体的磁化率分布设定所述网格单元的磁化率，得到单元体磁化率值；
根据所述地球主磁场三分量和所述单元体磁化率值得到每个所述网格单元中心点处的单元体磁化强度为：
其中，Mx、My、Mz分别表示单元体磁化强度的x、y、z分量。
4.根据权利要求3所述的方法，其特征在于，根据所述核函数单元积分系数的解析表达式得到核函数单元系数矩阵；所述核函数单元系数矩阵为Toeplitz矩阵，包括：
根据所述核函数单元积分系数的解析表达式得到核函数单元系数矩阵；所述核函数单元系数矩阵为Toeplitz矩阵；所述核函数单元系数矩阵的大小为
5.根据权利要求1至4任意一项所述的方法，其特征在于，不同网格单元的单元体磁化率值不同。
6.一种三维磁异常数快速正演装置，其特征在于，所述装置包括：
网格剖分模块，用于获取待计算区域在三维坐标系中的位置信息，将所述待计算区域沿x、y、z方向分别进行等间隔剖分，形成多个长方体的网格单元；
单元体磁化强度确定模块，用于根据地球主磁场模型，计算每个所述网格单元中心点处的地球主磁场三分量，根据待计算区域中三维地质体的磁化率分布设定所述网格单元的磁化率，得到单元体磁化率值；根据所述地球主磁场三分量和所述单元体磁化率值得到每个所述网格单元中心点处的单元体磁化强度；
核函数单元系数矩阵确定模块，用于以所述网格单元的中心点位置为观测点，当计算磁异常垂向分量时，获取预先构建的磁异常垂向分量核函数三个分量的离散积分算子为：
其中，μ0表示真空磁导率，(x,y,z)为观测点坐标，表示源点坐标；
(xi,yj,zk)表示编号为(i,j,k)的网格单元的中心坐标，表示编号为(l,m,n)的网格单元的中心坐标,其中，i,l＝0,1,…Nx-1；j,m＝0,1,…Ny-1；k,n＝0,1,…Nz-1，编号为(l,m,n)的单元三个方向积分区间分别为/>至/>至至/>
对所述磁异常垂向分量核函数三个分量的离散积分算子进行积分，得到核函数垂向分量三个分量的单元积分系数的解析表达式为：
其中，p,q,r＝1,2为整数，cpqr＝(-1)p(-1)q(-1)r， 根据所述核函数单元积分系数的解析表达式得到核函数单元系数矩阵；所述核函数单元系数矩阵为Toeplitz矩阵；
磁异常值计算模块，用于通过三维离散快速傅里叶变换进行所述核函数单元系数矩阵和所述单元体磁化强度三维离散卷积的快速相乘算法，得到所述待计算区域的磁异常值为：
其中，F和F-1分别表示三维离散傅里叶正反变换算子；表示磁化强度扩展矩阵，M表示网格剖分单元体的磁化强度矩阵，/>为扩展矩阵，/>表示提取矩阵的前Nx×Ny×Nz个元素，即整个三维网格单元上的磁异常值。
7.一种计算机设备，包括存储器和处理器，所述存储器存储有计算机程序，其特征在于，所述处理器执行所述计算机程序时实现权利要求1至5中任一项所述方法的步骤。
8.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现权利要求1至5中任一项所述的方法的步骤。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110978731.3A CN113673163B (zh) | 2021-08-25 | 2021-08-25 | 一种三维磁异常数快速正演方法、装置和计算机设备 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110978731.3A CN113673163B (zh) | 2021-08-25 | 2021-08-25 | 一种三维磁异常数快速正演方法、装置和计算机设备 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113673163A CN113673163A (zh) | 2021-11-19 |
| CN113673163B true CN113673163B (zh) | 2023-09-12 |
Family
ID=78545897
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110978731.3A Active CN113673163B (zh) | 2021-08-25 | 2021-08-25 | 一种三维磁异常数快速正演方法、装置和计算机设备 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113673163B (zh) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116244877B (zh) * | 2022-09-05 | 2023-11-14 | 中南大学 | 基于3d傅里叶变换的三维磁场数值模拟方法及系统 |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111400654A (zh) * | 2020-03-13 | 2020-07-10 | 中南大学 | 一种基于Toplitze核矩阵的重力场快速正演方法及反演方法 |
| CN111856598A (zh) * | 2020-06-29 | 2020-10-30 | 中国地质大学（武汉） | 一种磁测数据多层等效源上延拓与下延拓方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| KR101403296B1 (ko) * | 2013-12-09 | 2014-06-03 | 한국지질자원연구원 | 3차원 항공 자력 탐사 시스템 및 이를 이용한 3차원 항공 자력 탐사 방법 |
-
2021
- 2021-08-25 CN CN202110978731.3A patent/CN113673163B/zh active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111400654A (zh) * | 2020-03-13 | 2020-07-10 | 中南大学 | 一种基于Toplitze核矩阵的重力场快速正演方法及反演方法 |
| CN111856598A (zh) * | 2020-06-29 | 2020-10-30 | 中国地质大学（武汉） | 一种磁测数据多层等效源上延拓与下延拓方法 |
Non-Patent Citations (1)
| Title |
|---|
| 磁异常和梯度的频率域三维成像方法;崔亚彤;郭良辉;;物探与化探(第03期);144-152 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113673163A (zh) | 2021-11-19 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN112800657B (zh) | 基于复杂地形的重力场数值模拟方法、装置和计算机设备 | |
| US10338252B1 (en) | Systems and processes for building multiple equiprobable coherent geometrical models of the subsurface | |
| Čuma et al. | Large‐scale 3D inversion of potential field data | |
| CN112287534B (zh) | 基于nufft的二维磁异常快速正演模拟方法和装置 | |
| CN113420487B (zh) | 一种三维重力梯度张量数值模拟方法、装置、设备和介质 | |
| CN111967169B (zh) | 二度体重力异常积分解数值模拟方法和装置 | |
| CN109254327B (zh) | 三维强磁性体的勘探方法及勘探系统 | |
| CN113962077B (zh) | 三维各向异性强磁场数值模拟方法、装置、设备及介质 | |
| CN103278848A (zh) | 基于mpi并行预条件迭代的地震成像正演方法 | |
| Ouyang et al. | Iterative magnetic forward modeling for high susceptibility based on integral equation and Gauss-fast Fourier transform | |
| CN113673163B (zh) | 一种三维磁异常数快速正演方法、装置和计算机设备 | |
| CN114021408A (zh) | 二维强磁场数值模拟方法、装置、设备及介质 | |
| Fregoso et al. | Initializing cross-gradients joint inversion of gravity and magnetic data with a Bayesian surrogate gravity model | |
| Dai et al. | The forward modeling of 3D gravity and magnetic potential fields in space-wavenumber domains based on an integral method | |
| CN110490241A (zh) | 一种水平井参数优化方法及装置 | |
| CN113076678B (zh) | 一种频率域二度体重力异常快速数值模拟方法和装置 | |
| CN113221409B (zh) | 一种有限元和边界元耦合的声波二维数值模拟方法和装置 | |
| CN113779818B (zh) | 三维地质体其电磁场数值模拟方法、装置、设备及介质 | |
| CN113268702B (zh) | 一种频率域磁梯度张量变换方法、装置和计算机设备 | |
| CN115346008A (zh) | 工程勘察三维可视化处理方法、装置、设备及存储介质 | |
| Xu et al. | Solving fractional Laplacian visco-acoustic wave equations on complex-geometry domains using Grünwald-formula based radial basis collocation method | |
| Pletzer et al. | Mimetic interpolation of vector fields on Arakawa C/D grids | |
| CN113656976B (zh) | 一种二维磁梯度张量快速数值模拟方法、装置和设备 | |
| CN113806686B (zh) | 大规模复杂地质体重力梯度快速计算方法、装置和设备 | |
| CN113642189B (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 | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |