CN112287534A - 基于nufft的二维磁异常快速正演模拟方法和装置 - Google Patents

Abstract

本申请涉及一种基于NUFFT的二维磁异常快速正演模拟方法、装置、计算机设备和存储介质。本申请通过复杂磁性体模型表示、磁化强度计算、波数域磁位公式推导、一维形函数积分、非均匀波数选取、波数域磁异常场计算、空间域磁异常场计算等主要步骤，利用非均匀快速傅里叶变换和形函数积分实现了二维磁异常场高效、高精度正演模拟。在本申请中，充分利用非均匀快速傅里叶变换的灵活性和基于二次插值形函数积分的高精度优势，解决了波数域磁异常正演模拟局限于规则均匀网格，数值模拟不能同时兼顾计算精度和效率的问题。

Description

基于NUFFT的二维磁异常快速正演模拟方法和装置

技术领域

本申请涉及二维磁异常数值模拟技术领域，特别是涉及一种基于NUFFT的二维磁异常快速正演模拟方法、装置、计算机设备和存储介质。

背景技术

磁法勘探是一种非常经典的物探方法，已广泛应用于固体矿产勘探、石油天然气构造的普查以及水文工程地质等诸多方面，尤其是在矿产勘查中用磁法直接寻找磁铁矿及其共生的磁性矿产是必不可少的。随着勘探深度与难度的增加，实现高效、高精度的磁法勘探成为研究的重点。因此研究复杂条件下的磁异常场高效、高精度正演在磁异常反演成像与人机交互解释建模中有重要作用。

磁异常正演计算是根据磁化率的分布计算磁异常场的变化规律，反演则是根据观测的磁异常数据反推出磁化率的分布。正演是反演的基础，正演计算效率、计算精度直接影响反演计算效率和反演结果质量。对于磁异常正演计算，众多国内外学者进行了研究。对于任意截面形状、任意磁化率分布的二度体磁异常正演计算，一般分采用数值方法计算。基本思路是将复杂二度体模型剖分成许多规则二度体，而后利用规则二度体磁异常的累加来模拟复杂二度体磁异常。文献(刘增群，任意多边形截面二度体磁场正演程序。物探化探计算机技术，1991,13(4)：357-360)计算了一种磁化率为常数的情况下任意多边形截面二度体，将二度体剖分为若干个台阶的组合，利用解析公式计算台阶磁场，并进行累加，这种方法计算精度较高，但计算效率低，且只能计算磁化率为常数的情况；随着快速傅里叶变换扩边法在磁异常正演数值模拟中的广泛应用，波数域方法以其准确性和高效性成为处理大规模复杂模型正演的首选方法。文献(Tontini,F.C.,L.Cocchi,C.Carmisciano.Rapid 3-Dforward model of potential fields with application to the Palinuro Seamountmagnetic anomaly(southern Tyrrhenian Sea,Italy).Journal of GeophysicalResearch,2009.114.)采用快速傅里叶变换，给出了磁化率任意分布情况的磁场波数域表达式，借助快速傅里叶变换扩边法，实现了磁异常快速数值模拟，但在空间域和波数域只能均匀采样，并且由于截断边界效应的影响，采用扩边的方式使得数值模拟精度提高，但随着扩边比的增大计算时间成倍增加。文献(Song J.High-resolution 3-D radar imagingthrough nonuniform fast Fourier transform(NUFFT)[J].Commun Comput Phys,2006,1(1):176-191.)将非均匀快速傅里叶变换(NUFFT)应用到雷达高分辨率成像重建中，结果表明针对波数域非均匀网格的情况下NUFFT在保证计算效率的前提下比FFT波数域插值法成像精度更高。

因此，传统快速傅里叶变换扩边法在应用于波数域磁异常正演模拟时，只能在空间域和波数域均匀采样，局限于规则均匀网格，存在不能兼顾计算精度和效率的问题。

发明内容

基于此，有必要针对上述技术问题，提供一种能够兼顾计算精度与计算速度的基于NUFFT的二维磁异常快速正演模拟方法、装置、计算机设备和存储介质。

一种基于NUFFT的二维磁异常快速正演模拟方法，所述方法包括：

根据二维磁性体的大小，构建包含目标区域的二维磁性体模型；所述观测点坐标包括水平方向和垂向；所述二维磁性体模型是将所述目标区域在水平方向和垂向进行均匀或者非均匀网格剖分，并根据所述二维磁性体的截面形状设定网格剖分的每个节点的磁化率值；

根据所述二维磁性体模型，确定对应的空间域磁位积分表达式；所述空间域磁位积分表达式包含空间域磁化强度；所述空间域磁化强度是根据所述每个节点的磁化率值和已知的目标区域地球主磁场计算得到的；

对所述空间域磁位积分表达式沿水平方向进行一维傅里叶变换，得到波数域磁位一维积分表达式；其中，所述波数域磁位一维积分表达式包含波数域磁化强度；所述波数域磁化强度是对所述空间域磁化强度进行非均匀采样快速傅里叶变换得到的；

采用二次插值的形函数对所述波数域磁位空间域垂向一维积分进行计算，得到波数域磁位离散表达式；

根据水平方向网格剖分的尺寸信息，确定截止频率，根据所述截止频率确定波数选取范围，根据所述波数选取范围和预先设置的选取波数总数，得到波数选取值；

将所述波数选取值代入所述波数域磁位离散表达式，得到波数域磁位值；

根据傅里叶变换的微分性质，由所述波数域磁位值得到波数域磁异常场值；

通过对所述波数域磁异常场值进行一维非均匀快速傅里叶反变换，得到空间域磁异常场值。

在其中一个实施例中，还包括：获取目标区域的地球主磁场分量T_x(x_i,z_j)，T_z(x_i,z_j)。

根据所述每个节点的磁化率值和所述目标区域地球主磁场分量，得到空间域磁化强度为：

M_x(x_i,z_j)＝χ(x_i,z_j)T_x(x_i,z_j)

M_z(x_i,z_j)＝χ(x_i,z_j)T_z(x_i,z_j)

式中，χ(x_i,z_j)表示编号为(x_i,z_j)节点的磁化率值；T_x(x_i,z_j)，T_z(x_i,z_j)分别表示(x_i,z_j)处地球主磁场的x、z分量；M_x(x_i,z_j)、M_z(x_i,z_j)分别表示(x_i,z_j)处磁化强度的x、z分量。

在其中一个实施例中，还包括：空间域磁位积分表达式为：

式中：(x,z)和(x′,z′)分别表示观测点和场源任意一点的坐标；U为磁性体在空间任意一点产生的磁位；π为圆周率常数；μ₀为真空磁导率，μ₀＝4π×10^-7H/m。

对所述空间域磁位积分表达式沿水平方向进行一维傅里叶变换，得到所述波数域磁位一维积分表达式为：

式中：

表示波数域磁位；k表示波数；

分别表示波数域磁化强度的x，z分量；sign表示符号函数：

在其中一个实施例中，还包括：非均匀采样快速傅里叶变换为：

式中：i为虚数单位，α_n为给定离散采样点x_n∈[-N/2,N/2]对应采样点的值，F(α)_j为计算的离散采样点{α_n}傅里叶变换频谱，N表示选取波数总数；

其中，所述非均匀快速傅里叶变换的实现步骤为：

根据x_n临近q个均匀采样点的傅里叶变换基，得到近似非均匀傅里叶变换基

为：

式中：m表示过采样因子，ω_l表示权重因子，s_j为精度因子，[mx_n]表示对mx_n取整；

根据采样点的值和权重因子，计算新傅里叶变换基对应的傅里叶变换系数η_p：

采用均匀的快速傅里叶变换，得到：

式中，F_j表示傅立叶变换之后的频谱。

在其中一个实施例中，还包括：采用二次插值的形函数对所述波数域空间域垂向一维积分进行计算，得到波数域磁位离散表达式，包括：

对所述波数域磁位一维积分表达式进行离散，得到：

式中，I_x(k,z)、I_z(k,z)分别为：

式中，z_j，z_m分别表示单元积分的下限和上限；

将波数域磁场强度表示成单元节点磁化强度的二次插值函数为：

式中，N_j，N_p，N_m表示形函数，将

的表达式带入I_x(k,z)，I_z(k,z)的表达式中，得到：

式中，W′_j、W′_p、W′_m为形函数单元积分，

当波数k≠0，即γ≠0时：

当波数k＝0，即γ＝0时：

在其中一个实施例中，还包括：根据水平方向网格剖分的尺寸信息，确定截止频率为：

式中，k_max表示截止频率；min(Δx)表示水平方向网格剖分矩形尺寸的最小值；

根据所述截止频率确定波数域波数选取范围为(-k_max,k_max)；

根据所述波数选取范围得到波数选取值，包括：

在(0,k_max)上波数依次取值为：

在(-k_max,0)上波数依次取值为：

式中，k_min＝10^-4，N表示选取波数总数。

在其中一个实施例中，还包括：根据傅里叶变换的微分性质，由所述波数域磁位值得到波数域磁异常场值为：

式中，

分别表示波数域磁异常场的x，z分量，

表示波数域磁位值，k表示波数。

一种基于NUFFT的二维磁异常快速正演模拟装置，所述装置包括：

二维磁性体模型构建模块，用于根据二维磁性体的大小，构建包含目标区域的二维磁性体模型；所述观测点坐标包括水平方向和垂向；所述二维磁性体模型是将所述目标区域在水平方向和垂向进行均匀或者非均匀网格剖分，并根据所述二维磁性体的截面形状设定网格剖分的每个节点的磁化率值；

空间域磁位积分表达式确定模块，用于根据所述二维磁性体模型，确定对应的空间域磁位积分表达式；所述空间域磁位积分表达式包含空间域磁化强度；所述空间域磁化强度是根据所述每个节点的磁化率值和已知的目标区域地球主磁场计算得到的；

波数域磁位一维积分表达式确定模块，用于对所述空间域磁位积分表达式沿水平方向进行一维傅里叶变换，得到波数域磁位一维积分表达式；其中，所述波数域磁位一维积分表达式包含波数域磁化强度；所述波数域磁化强度是对所述空间域磁化强度进行非均匀采样快速傅里叶变换得到的；

波数域磁位离散表达式确定模块，用于采用二次插值的形函数对所述波数域空间域垂向一维积分进行计算，得到波数域磁位离散表达式；

波数选取值确定模块，用于根据水平方向网格剖分的尺寸信息，确定截止频率，根据所述截止频率确定波数域波数选取范围，根据所述波数选取范围和预先设置的选取波数总数，得到波数选取值；

波数域磁位值确定模块，用于将所述波数选取值代入所述波数域磁位离散表达式，得到波数域磁位值；

波数域磁异常场值确定模块，用于根据傅里叶变换的微分性质，由所述波数域磁位值得到波数域磁异常场值；

空间域异常场值确定模块，用于通过对所述波数域磁异常场值进行一维非均匀快速傅里叶反变换，得到目空间域磁异常场值。

一种计算机设备，包括存储器和处理器，所述存储器存储有计算机程序，所述处理器执行所述计算机程序时实现以下步骤：

根据二维磁性体的大小，构建包含目标区域的二维磁性体模型；所述观测点坐标包括水平方向和垂向；所述二维磁性体模型是将所述目标区域在水平方向和垂向进行均匀或者非均匀网格剖分，并根据所述二维磁性体的截面形状设定网格剖分的每个节点的磁化率值；

根据所述二维磁性体模型，确定对应的空间域磁位积分表达式；所述空间域磁位积分表达式包含空间域磁化强度；所述空间域磁化强度是根据所述每个节点的磁化率值和已知的目标区域地球主磁场计算得到的；

对所述空间域磁位积分表达式沿水平方向进行一维傅里叶变换，得到波数域磁位一维积分表达式；其中，所述波数域磁位一维积分表达式包含波数域磁化强度；所述波数域磁化强度是对所述空间域磁化强度进行非均匀采样快速傅里叶变换得到的；

采用二次插值的形函数对所述波数域磁位空间域垂向一维积分进行计算，得到波数域磁位离散表达式；

根据水平方向网格剖分的尺寸信息，确定截止频率，根据所述截止频率确定波数选取范围，根据所述波数选取范围和预先设置的选取波数总数，得到波数选取值；

将所述波数选取值代入所述波数域磁位离散表达式，得到波数域磁位值；

根据傅里叶变换的微分性质，由所述波数域磁位值得到波数域磁异常场值；

通过对所述波数域磁异常场值进行一维非均匀快速傅里叶反变换，得到空间域磁异常场值。

一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现以下步骤：

根据二维磁性体的大小，构建包含目标区域的二维磁性体模型；所述观测点坐标包括水平方向和垂向；所述二维磁性体模型是将所述目标区域在水平方向和垂向进行均匀或者非均匀网格剖分，并根据所述二维磁性体的截面形状设定网格剖分的每个节点的磁化率值；

根据所述二维磁性体模型，确定对应的空间域磁位积分表达式；所述空间域磁位积分表达式包含空间域磁化强度；所述空间域磁化强度是根据所述每个节点的磁化率值和已知的目标区域地球主磁场计算得到的；

对所述空间域磁位积分表达式沿水平方向进行一维傅里叶变换，得到波数域磁位一维积分表达式；其中，所述波数域磁位一维积分表达式包含波数域磁化强度；所述波数域磁化强度是对所述空间域磁化强度进行非均匀采样快速傅里叶变换得到的；

采用二次插值的形函数对所述波数域磁位空间域垂向一维积分进行计算，得到波数域磁位离散表达式；

根据水平方向网格剖分的尺寸信息，确定截止频率，根据所述截止频率确定波数选取范围，根据所述波数选取范围和预先设置的选取波数总数，得到波数选取值；

将所述波数选取值代入所述波数域磁位离散表达式，得到波数域磁位值；

根据傅里叶变换的微分性质，由所述波数域磁位值得到波数域磁异常场值；

通过对所述波数域磁异常场值进行一维非均匀快速傅里叶反变换，得到空间域磁异常场值。

上述基于NUFFT的二维磁异常快速正演模拟方法、装置、计算机设备和存储介质，通过复杂磁性体模型表示、磁化强度计算、波数域磁位公式推导、一维形函数积分、非均匀波数选取、波数域磁异常场计算、空间域磁异常场计算等主要步骤，利用非均匀快速傅里叶变换和形函数积分实现了二维磁异常场高效、高精度正演模拟。在本方法中，充分利用非均匀快速傅里叶变换的灵活性和基于二次插值形函数积分的高精度优势，解决了波数域磁异常正演模拟局限于规则均匀网格，数值模拟不能同时兼顾计算精度和效率的问题。

附图说明

图1为一个实施例中基于NUFFT的二维磁异常快速正演模拟方法的流程示意图；

图2为一个实施例中基于NUFFT的二维磁异常快速正演模拟方法的算法流程图；

图3为一个实施例中实现非均匀快速傅里叶变换的流程示意图；

图4为一个实施例中截面为矩形的二维磁性体模型示意图；

图5为一个实施例中空间域均匀剖分，波数域非均匀采样磁异常场水平和垂向分量模拟结果与解析解对比图；

图6为一个实施例中空间域均匀剖分，波数域非均匀采样磁异常场水平和垂向分量模拟结果与解析解的相对误差；

图7为空间域和波数域都非均匀采样磁异常场水平和垂向分量模拟结果与解析解对比图；

图8为空间域和波数域都非均匀采样磁异常场水平和垂向分量模拟结果与解析解的相对误差；

图9为一个实施例中基于NUFFT的二维磁异常快速正演模拟装置的结构框图；

图10为一个实施例中计算机设备的内部结构图。

具体实施方式

为了使本申请的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本申请进行进一步详细说明。应当理解，此处描述的具体实施例仅仅用以解释本申请，并不用于限定本申请。

本申请提供的基于NUFFT的二维磁异常快速正演模拟方法，可以应用于如下应用环境中。根据观测点构建包含目标区域的二维磁性体模型，确定对应的空间域磁位积分表达式，其中空间域磁位积分表达式包含空间域磁化强度，空间域磁化强度是根据所述每个节点的磁化率值和已知的目标区域地球主磁场计算得到的。对空间域磁位积分表达式沿水平方向进行一维傅里叶变换，得到波数域磁位一维积分表达式；其中，波数域磁位一维积分表达式包含波数域磁化强度，波数域磁化强度是对所述空间域磁化强度进行非均匀采样快速傅里叶变换得到的。采用二次插值的形函数对所述波数域磁位空间域垂向一维积分进行计算，得到波数域磁位离散表达式。根据水平方向网格剖分的尺寸信息，确定截止频率，根据截止频率确定波数域波数选取范围，根据波数选取范围和预先设置的选取波数总数，在波数域进行非均匀波数选取，得到波数选取值；将波数选取值代入波数域磁位离散表达式，得到波数域磁位值；根据傅里叶变换的微分性质，由波数域磁位值得到波数域磁异常场值；通过对波数域磁异常场值进行一维非均匀快速傅里叶反变换，得到空间域磁异常场值。

在一个实施例中，如图1所示，提供了一种基于NUFFT的二维磁异常快速正演模拟方法，包括以下步骤：

步骤102，根据二维磁性体的大小，构建包含目标区域的二维磁性体模型。

观测点坐标包括水平方向和垂向；二维磁性体模型是将目标区域在水平方向和垂向进行均匀或者非均匀网格剖分，并根据二维磁性体的截面形状设定网格剖分的每个节点的磁化率值。

构建包含目标区域的二维磁性体模型，将目标区域在水平方向和垂向进行均匀或者非均匀网格剖分，网格剖分灵活，在重点勘探区域水平和垂向网格适当加密，在远离目标区域水平和垂向网格适当稀疏。根据二维磁性体截面形状，对网格剖分的每个节点进行赋值，每个节点的磁化率是一个常数，不同节点的磁化率值不同，以此刻画任意磁化率分布的复杂二维磁性体模型。

步骤104，根据二维磁性体模型，确定对应的空间域磁位积分表达式。

空间域磁位积分表达式包含空间域磁化强度；空间域磁化强度是根据每个节点的磁化率值和已知的目标区域地球主磁场计算得到的。

磁异常即地磁异常，又称磁力异常。地磁场的理论分布是有变化的。而实际上测得的地球磁场强度和理论磁场强度是有区别的，这种区别称为地磁异常。它主要是由地壳内磁性不同的岩石受地磁场磁化而产生的附加磁场。

磁位的定义是把单位强度的磁极从参考点，通常是无穷远，移至所考虑的一点时为反抗磁场而必须做的功。

步骤106，对空间域磁位积分表达式沿水平方向进行一维傅里叶变换，得到波数域磁位一维积分表达式。

其中，波数域磁位一维积分表达式包含波数域磁化强度；波数域磁化强度是对空间域磁化强度进行非均匀采样快速傅里叶变换得到的。

由于构建的维磁性体模型中剖分的尺寸不一定一样，每个网格剖分节点的磁化率值不一定一样，因此，空间域的采样可以是非均匀的。当网格剖分是非均匀的时候，在对空间域磁位积分表达式沿水平方向进行一维傅里叶变换时，波数域磁化强度是对空间域磁化强度进行非均匀采样快速傅里叶变换得到的，可以实现在波数域进行非均匀波数选取，提高计算效率。

步骤108，采用二次插值的形函数对波数域磁位空间域z轴方的一维积分进行计算，得到波数域磁位离散表达式。

保留垂向在空间域，便于在磁异常场变化快的地方加密网格，在磁异常场变化慢的地方稀疏网格，有利于精确模拟复杂地形，并且可以在保证计算精度的前提下尽可能减少计算量。波数域磁位一维积分表达式在垂向是积分形式，通过将垂向剖分成N_z个单元，在每个单元内磁化强度满足二次变化，可以采用基于二次插值的形函数进行一维垂向积分，每个单元内可以推导出近似解，积分精度高。形函数是一种连续函数，满足单元内边界点的给定值和内部连续，可以把单元内任一点的位移用节点位移表示，基于二次插值的形函数可以用三个节点来确定单元内的任一点的位移。

步骤110，根据水平方向网格剖分的尺寸信息，确定截止频率，根据截止频率确定波数域波数选取范围，根据波数选取范围和预先设置的选取波数总数，得到波数选取值。

理论上，采样点选取越多计算精度越高，但同时计算量也相应增加，效率下降。因此，需要在满足计算精度的前提下尽可能减少选取波数总数。通过对不同模型磁异常频谱分析发现，波数选取值绝对值比较小的时候磁异常频谱比较大，随着波数选取值绝对值的增大，磁异常频谱逐渐减小，直到波数选取值绝对值足够大的时候磁异常频谱几乎衰减为零，磁异常频谱几乎衰减为零说明对应频率的磁异常场能量是极低的，因此不需要选取绝对值过大的波数选取值。根据水平方向网格剖分的尺寸信息，确定截止频率，根据截止频率确定波数域波数选取范围，可以在满足精度的前提下尽可能减少选取波数总数，提高方法的计算效率。

步骤112，将波数选取值代入波数域磁位离散表达式，得到波数域磁位值。

相对于空间域磁位表达式，波数域磁位表达式具有计算效率高，准确性高的优点。

步骤114，根据傅里叶变换的微分性质，由波数域磁位值得到波数域磁异常场值。

傅里叶变换的微分性质是一个函数导数的傅里叶变换等于这个函数傅里叶变换乘以因子iw。由于磁场是磁位的一阶导数，可以根据傅里叶变换的微分性质，直接由波数域磁位值得到波数域磁场值。

步骤116，通过对波数域磁异常场值进行一维非均匀快速傅里叶反变换，得到空间域磁异常场值。

通过一维非均匀快速傅里叶反变换，将波数域磁异常场值转换到空间域磁异常场值，可以得到空间域磁异常场值。

上述基于NUFFT的二维磁异常快速正演模拟方法中，通过复杂磁性体模型表示、磁化强度计算、波数域磁位公式推导、一维形函数积分、非均匀波数选取、波数域磁异常场计算、空间域磁异常场计算等主要步骤，利用非均匀快速傅里叶变换和形函数积分实现了二维磁异常场高效、高精度正演模拟。在本方法中，充分利用非均匀快速傅里叶变换的灵活性和基于二次插值形函数积分的高精度优势，解决了波数域磁异常正演模拟局限于规则均匀网格，数值模拟不能同时兼顾计算精度和效率的问题。

在其中一个实施例中，获取目标区域的地球主磁场分量T_x(x_i,z_j)，T_z(x_i,z_j)；

根据每个节点的磁化率值和所述目标区域地球主磁场分量，得到空间域磁化强度为：

M_x(x_i,z_j)＝χ(x_i,z_j)T_x(x_i,z_j)

M_z(x_i,z_j)＝χ(x_i,z_j)T_z(x_i,z_j)

式中，χ(x_i,z_j)表示编号为(x_i,z_j)节点的磁化率值；T_x(x_i,z_j)，T_z(x_i,z_j)分别表示(x_i,z_j)处地球主磁场的x、z分量；M_x(x_i,z_j)、M_z(x_i,z_j)分别表示(x_i,z_j)处磁化强度的x、z分量。

在其中一个实施例中，空间域磁位积分表达式为：

式中：(x,z)和(x′,z′)分别表示观测点和场源任意一点的坐标；U为磁性体在空间任意一点产生的磁位；π为圆周率常数；μ₀为真空磁导率，μ₀＝4π×10^-7H/m。

对本实施例中空间域磁位积分表达式沿水平方向进行一维傅里叶变换，得到波数域磁位一维积分表达式为：

式中：

表示波数域磁位；k表示波数；

分别表示波数域磁化强度的x，z分量；sign表示符号函数：

其中，波数域磁场强度

是对空间域磁场强度M_x(x,z)，M_z(x,z)进行非均匀采样快速傅里叶变换得到的。

由于构建的维磁性体模型中剖分的尺寸不一定一样，每个网格剖分节点的磁化率值不一定一样，因此，空间域的采样可以是非均匀的。通过对空间域磁位积分表达式沿水平方向进行一维傅里叶变换，得到波数域磁位一维积分表达式，当网格剖分是非均匀的时候，在对空间域磁位积分表达式沿水平方向进行一维傅里叶变换时，波数域磁化强度是对空间域磁化强度进行非均匀采样快速傅里叶变换得到的。

在其中一个实施例中，如图2所示，一种基于NUFFT的二维磁异常快速正演模拟方法主要包括七个部分：复杂磁性体模型表示、磁化强度计算、波数域磁位公式推导、一维形函数积分、非均匀波数选取、波数域磁异常场计算和空间域磁异常场计算。复杂磁性体模型表示包括设定目标区域及网格剖分大小、根据二维磁性体截面几何形状，对剖分网格节点磁化率赋值；磁化强度计算包括：根据地球物理实测数据，给定目标区域背景磁场、根据磁化率分布和主磁场，计算磁化强度分量；波数域磁位公式推导包括：采用一维快速傅里叶变换，将空间域磁位转换到波数域、采用一维非均匀快速傅里叶变换，将空间域磁化强度转换到波数域；一维形函数积分是采用二次插值的形函数计算空间域垂向一维积分；非均匀波数选取是根据磁异常频谱变化规律进行非均匀波数选取；波数域磁异常场计算是由傅里叶变换的微分性质，计算波数域磁异常场值；空间域磁异常场计算是对计算的波数域磁异常场进行一维非均匀快速傅里叶反变换，得到空间域磁异常场。

在其中一个实施例中，非均匀采样快速傅里叶变换为：

式中：i为虚数单位，α_n为给定离散采样点x_n∈[-N/2,N/2]对应采样点的值，F(α)_j为计算的离散采样点{α_n}傅里叶变换频谱，N表示选取波数总数；

如图3所示，其中，实现非均匀快速傅里叶变换的步骤为：

步骤302：根据x_n临近q个均匀采样点的傅里叶变换基，得到近似非均匀傅里叶变换基

为：

式中：m表示过采样因子，ω_l表示权重因子，s_j为精度因子，[mx_n]表示对mx_n取整；

步骤304：根据采样点的值和权重因子，计算新傅里叶变换基下对应的傅里叶变换系数η_p：

步骤306：采用均匀的快速傅里叶变换，得到：

式中，F_j表示傅立叶变换之后的频谱。

快速傅里叶变换扩边法在磁异常正演数值模拟中有着广泛应用，但这种方法在磁异常正演数值模拟的应用中只能均匀采样，难以实现计算精度和计算效率的统一。而本实施例中的非均匀快速傅里叶变换是基于将局部插值与FFT相结合来实现非均匀采样的快速傅里叶变换，将本实施例中的非均匀快速傅里叶变换应用于磁异常正演数值模拟中，可以实现波数域的非均匀采样，进而在满足精度的前提下提高计算效率。

在一个实施例中，采用二次插值的形函数对波数域磁位一维积分表达式中空间域垂向的一维积分进行计算，得到波数域磁位离散表达式，包括：

对波数域磁位一维积分表达式进行离散，得到：

式中，I_x(k,z)、I_z(k,z)分别为：

式中，z_j，z_m分别表示单元积分的下限和上限；

将波数域磁场强度表示成单元节点磁化强度的二次插值函数为：

式中，N_j，N_p，N_m表示形函数；

将

的表达式带入I_x(k,z)，I_z(k,z)的表达式中，得到：

式中，W′_j、W′_p、W′_m为形函数单元积分，

当波数k≠0，即γ≠0时：

当波数k＝0，即γ＝0时：

基于二次插值的形函数积分具有高精度的优势。

在一个实施例中，根据水平方向网格剖分的尺寸信息，确定截止频率为：

式中，k_max表示截止频率；min(Δx)表示水平方向网格剖分矩形尺寸的最小值；

根据所述截止频率确定波数域波数选取范围为(-k_max,k_max)；

根据所述波数选取范围得到波数选取值，包括：

在(0,k_max)上波数依次取值为：

在(-k_max,0)上波数依次取值为：

式中，k_min＝10^-4，N表示选取波数总数。

本实施例中的波数选取方式是将波数范围(-k_max,k_max)在对数轴上划分为N个区间，每个区间中波数范围的对数距离为1，采用该方式可以实现在波数绝对值较小的区域密集采样，随着波数绝对值增大，采样逐渐稀疏。通过这种方式可以在满足精度的前提下减少采样点的数量，进而提高计算效率。

在一个实施例中，根据傅里叶变换的微分性质，根据波数域磁位值得到波数域磁异常场值为：

式中，

分别表示波数域磁异常场的x，z分量，

表示波数域磁位值，k表示波数。

在一个具体实施例中，研究区域有一个截面为矩形的模型如图4所示，研究区域范围：水平方向从-600m到600m，垂向从0m到500m。网格数为200×100，水平方向网格间隔为6m，垂向网格间隔为5m。异常体范围沿水平方向从-120m到120m，垂向从250m到400m，异常体的磁化率为0.01，目标区域地球主磁场强度为45000nT，磁倾角为60°，磁偏角为0°。计算水平地面上201个观测点的磁异常。

利用Fortran语言编程实现，运行程序所用的个人电脑配置为：CPU-Inter Corei7-8700，主频为3.2GHz，运行内存为8.00GB。图5为空间域均匀剖分，并采用非均匀波数选取方法计算的截面为矩形异常体产生的磁异常场水平和垂直分量解析解和数值解的结果。从图5中可以看出采用NUFFT分别选取201和401个波数的磁异常场水平和垂直分量模拟结果与解析解吻合的很好，通过图6相对误差图可以看出，磁异常场水平和垂向分量相对误差随着选取波数的增加，计算精度也随之增加。同样对于该模型统计了FFT扩边法与NUFFT方法计算的相对均方根误差及计算时间对比见表1，从表1可以看出在达到相同精度(rrms<0.5％)的情况下，NUFFT相对于FFT扩边法计算效率提高了约5.2倍。

表1截面为矩形模型不同模拟方法统计的相对均方根误差及计算时间对比

其中，4L、8L分别为扩边比，201，401分别为波数域采样点个数，rrms(％)表示相对均方根误差。

B^th表示解析解，B表示数值解。

在另一个具体实施例中，采用图4所示的二维磁异常体模型，但在空间域和波数域都均匀采样，空间域水平方向(0m,150m)网格间隔为3m，在(150m,300m)网格间隔为6m，(300m,600m)网格间隔为12m,在(-600m,0m)关于x＝0对称，垂向网格均匀剖分，网格剖分大小为200×100。图7为空间域和波数域都非均匀采样计算的磁异常水平和垂向分量数值模拟结果，其中，波数采用非均匀波数采样方法分别选201和401个波数，计算结果如图7所示。从图7可以看出解析解与数值模拟结果吻合较好，由图8相对误差可以看出，随着选取波数的增加，计算精度也随之增加，但是与上一个具体实施例的计算结果相比，空间域和波数域都利用非均匀采样傅里叶变换计算精度更高。从而验证了在重点区域适当加密网格，在远离异常体区域网格适当稀疏，可以在保证勘探精度的同时提高计算效率。

应该理解的是，虽然图1的流程图中的各个步骤按照箭头的指示依次显示，但是这些步骤并不是必然按照箭头指示的顺序依次执行。除非本文中有明确的说明，这些步骤的执行并没有严格的顺序限制，这些步骤可以以其它的顺序执行。而且，图1中的至少一部分步骤可以包括多个子步骤或者多个阶段，这些子步骤或者阶段并不必然是在同一时刻执行完成，而是可以在不同的时刻执行，这些子步骤或者阶段的执行顺序也不必然是依次进行，而是可以与其它步骤或者其它步骤的子步骤或者阶段的至少一部分轮流或者交替地执行。

在一个实施例中，如图9所示，提供了一种基于NUFFT的二维磁异常快速正演模拟装置，包括：二维磁性体模型构建模块902、空间域磁位积分表达式确定模块904、波数域磁位一维积分表达式确定模块906、波数域磁位离散表达式确定模块908、波数选取值确定模块910、波数域磁位值确定模块912、波数域磁异常场值确定模块914和空间域异常场值确定模块916，其中：

二维磁性体模型构建模块902，用于根据二维磁性体的大小，构建包含目标区域的二维磁性体模型；观测点坐标包括水平方向和垂向；二维磁性体模型是将所述目标区域在水平方向和垂向进行均匀或者非均匀网格剖分，并根据二维磁性体的截面形状设定网格剖分的每个节点的磁化率值；

空间域磁位积分表达式确定模块904，用于根据二维磁性体模型，确定对应的空间域磁位积分表达式；空间域磁位积分表达式包含空间域磁化强度；空间域磁化强度是根据每个节点的磁化率值和已知的目标区域地球主磁场计算得到的；

波数域磁位一维积分表达式确定模块906，用于对空间域磁位积分表达式沿水平方向进行一维傅里叶变换，得到波数域磁位一维积分表达式；其中，波数域磁位一维积分表达式包含波数域磁化强度；波数域磁化强度是对空间域磁化强度进行非均匀采样快速傅里叶变换得到的；

波数域磁位离散表达式确定模块908，用于采用二次插值的形函数对波数域磁位空间域垂向一维积分进行计算，得到波数域磁位离散表达式；

波数选取值确定模块910，用于根据水平方向网格剖分的尺寸信息，确定截止频率，根据截止频率确定波数域波数选取范围，根据波数选取范围和预先设置的选取波数总数，在波数域进行非均匀波数选取，得到波数选取值；

波数域磁位值确定模块912，用于将波数选取值代入波数域磁位离散表达式，得到波数域磁位值；

波数域磁异常场值确定模块914，用于根据傅里叶变换的微分性质，由波数域磁位值得到波数域磁异常场值；

空间域异常场值确定模块916，用于通过对波数域磁异常场值进行一维非均匀快速傅里叶反变换，得到空间域磁异常场值。

波数选取值确定模块908还包括根据水平方向网格剖分的尺寸信息，确定截止频率为：

式中，k_max表示截止频率；min(Δx)表示水平方向网格剖分矩形尺寸的最小值。

根据截止频率确定波数域波数选取范围为(-k_max,k_max)；

根据波数选取范围得到波数选取值，包括：

在(0,k_max)上波数依次取值为：

在(-k_max,0)上波数依次取值为：

式中，k_min＝10^-4，N表示选取波数总数。

波数域磁异常场值确定模块914还包括根据傅里叶变换的微分性质，由所述波数域磁位值得到波数域磁异常场值为：

式中，

分别表示波数域磁异常场的x，z分量，

表示波数域磁位值，k表示波数。

关于基于NUFFT的二维磁异常快速正演模拟装置的具体限定可以参见上文中对于基于NUFFT的二维磁异常快速正演模拟方法的限定，在此不再赘述。上述基于NUFFT的二维磁异常快速正演模拟装置中的各个模块可全部或部分通过软件、硬件及其组合来实现。上述各模块可以硬件形式内嵌于或独立于计算机设备中的处理器中，也可以以软件形式存储于计算机设备中的存储器中，以便于处理器调用执行以上各个模块对应的操作。

在一个实施例中，提供了一种计算机设备，该计算机设备可以是终端，其内部结构图可以如图10所示。该计算机设备包括通过系统总线连接的处理器、存储器、网络接口、显示屏和输入装置。其中，该计算机设备的处理器用于提供计算和控制能力。该计算机设备的存储器包括非易失性存储介质、内存储器。该非易失性存储介质存储有操作系统和计算机程序。该内存储器为非易失性存储介质中的操作系统和计算机程序的运行提供环境。该计算机设备的网络接口用于与外部的终端通过网络连接通信。该计算机程序被处理器执行时以实现一种基于NUFFT的二维磁异常快速正演模拟方法。该计算机设备的显示屏可以是液晶显示屏或者电子墨水显示屏，该计算机设备的输入装置可以是显示屏上覆盖的触摸层，也可以是计算机设备外壳上设置的按键、轨迹球或触控板，还可以是外接的键盘、触控板或鼠标等。

本领域技术人员可以理解，图10中示出的结构，仅仅是与本申请方案相关的部分结构的框图，并不构成对本申请方案所应用于其上的计算机设备的限定，具体的计算机设备可以包括比图中所示更多或更少的部件，或者组合某些部件，或者具有不同的部件布置。

在一个实施例中，提供了一种计算机设备，包括存储器和处理器，该存储器存储有计算机程序，该处理器执行计算机程序时实现上述方法实施例中的步骤。

在一个实施例中，提供了一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现上述方法实施例中的步骤。

本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的计算机程序可存储于一非易失性计算机可读取存储介质中，该计算机程序在执行时，可包括如上述各方法的实施例的流程。其中，本申请所提供的各实施例中所使用的对存储器、存储、数据库或其它介质的任何引用，均可包括非易失性和/或易失性存储器。非易失性存储器可包括只读存储器(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.一种基于NUFFT的二维磁异常快速正演模拟方法，其特征在于，所述方法包括：

根据二维磁性体的大小，构建包含目标区域的二维磁性体模型；所述观测点坐标包括水平方向和垂向；所述二维磁性体模型是将所述目标区域在水平方向和垂向进行均匀或者非均匀网格剖分，并根据所述二维磁性体的截面形状设定网格剖分的每个节点的磁化率值；

根据所述二维磁性体模型，确定对应的空间域磁位积分表达式；所述空间域磁位积分表达式包含空间域磁化强度；所述空间域磁化强度是根据所述每个节点的磁化率值和已知的目标区域地球主磁场计算得到的；

对所述空间域磁位积分表达式沿水平方向进行一维傅里叶变换，得到波数域磁位一维积分表达式；其中，所述波数域磁位一维积分表达式包含波数域磁化强度；所述波数域磁化强度是对所述空间域磁化强度进行非均匀采样快速傅里叶变换得到的；

采用二次插值的形函数对所述波数域磁位空间域垂向一维积分进行计算，得到波数域磁位离散表达式；

根据水平方向网格剖分的尺寸信息，确定截止频率，根据所述截止频率确定波数选取范围，根据所述波数选取范围和预先设置的选取波数总数，得到波数选取值；

将所述波数选取值代入所述波数域磁位离散表达式，得到波数域磁位值；

根据傅里叶变换的微分性质，由所述波数域磁位值得到波数域磁异常场值；

通过对所述波数域磁异常场值进行一维非均匀快速傅里叶反变换，得到空间域磁异常场值。

2.根据权利要求1所述的方法，其特征在于，根据所述每个节点的磁化率值和已知的目标区域地球主磁场计算得到所述空间域磁化强度，包括：

获取目标区域的地球主磁场分量T_x(x_i,z_j)，T_z(x_i,z_j)；

根据所述每个节点的磁化率值和所述目标区域地球主磁场分量，得到空间域磁化强度为：

M_x(x_i,z_j)＝χ(x_i,z_j)T_x(x_i,z_j)

M_z(x_i,z_j)＝χ(x_i,z_j)T_z(x_i,z_j)

式中，χ(x_i,z_j)表示编号为(x_i,z_j)节点的磁化率值；T_x(x_i,z_j)，T_z(x_i,z_j)分别表示(x_i,z_j)处地球主磁场的x、z分量；M_x(x_i,z_j)、M_z(x_i,z_j)分别表示(x_i,z_j)处磁化强度的x、z分量。

3.根据权利要求2所述的方法，其特征在于，对所述空间域磁位积分表达式沿水平方向进行一维傅里叶变换，得到波数域磁位一维积分表达式，包括：

所述空间域磁位积分表达式为：

式中：(x,z)和(x′,z′)分别表示观测点和场源任意一点的坐标；U为磁性体在空间任意一点产生的磁位；π为圆周率常数；μ₀为真空磁导率，μ₀＝4π×10^-7H/m；

对所述空间域磁位积分表达式沿水平方向进行一维傅里叶变换，得到所述波数域磁位一维积分表达式为：

式中：

表示波数域磁位；k表示波数；

分别表示波数域磁化强度的x，z分量；sign表示符号函数：

4.根据权利要求1所述的方法，其特征在于，所述非均匀采样快速傅里叶变换为：

式中：i为虚数单位，α_n为给定离散采样点x_n∈[-N/2,N/2]对应采样点的值，F(α)_j为计算的离散采样点{α_n}傅里叶变换频谱，N表示选取波数总数；

其中，所述非均匀快速傅里叶变换的实现步骤为：

根据x_n临近q个均匀采样点的傅里叶变换基，得到近似非均匀傅里叶变换基

为：

式中：m表示过采样因子，ω_l表示权重因子，s_j为精度因子，[mx_n]表示对mx_n取整；

根据采样点的值和权重因子，计算新傅里叶变换基对应的傅里叶变换系数η_p：

采用均匀的快速傅里叶变换，得到：

式中，F_j表示傅立叶变换之后的频谱。

5.根据权利要求4所述的方法，其特征在于，所述采用二次插值的形函数对所述波数域空间域垂向一维积分进行计算，得到波数域磁位离散表达式，包括：

对所述波数域磁位一维积分表达式进行离散，得到：

式中，I_x(k,z)、I_z(k,z)分别为：

式中，z_j，z_m分别表示单元积分的下限和上限；

将波数域磁场强度表示成单元节点磁化强度的二次插值函数为：

式中，N_j，N_p，N_m表示形函数；

将

的表达式带入I_x(k,z)，I_z(k,z)的表达式中，得到：

式中，W′_j、W′_p、W′_m为形函数单元积分，

当波数k≠0，即γ≠0时：

当波数k＝0，即γ＝0时：

6.根据权利要求1所述的方法，其特征在于，根据水平方向网格剖分的尺寸信息，确定截止频率，根据所述截止频率确定波数域波数选取范围，根据所述波数选取范围和预先设置的选取波数总数，得到波数选取值，包括：

根据水平方向网格剖分的尺寸信息，确定截止频率为：

式中，k_max表示截止频率；min(Δx)表示水平方向网格剖分矩形尺寸的最小值；

根据所述截止频率确定波数域波数选取范围为(-k_max,k_max)；

根据所述波数选取范围得到波数选取值，包括：

在(0,k_max)上波数依次取值为：

在(-k_max,0)上波数依次取值为：

式中，k_min＝10^-4，N表示选取波数总数。

7.根据权利要求6所述的方法，其特征在于，根据傅里叶变换的微分性质，由所述波数域磁位值得到波数域磁异常场值，包括：

根据傅里叶变换的微分性质，由所述波数域磁位值得到波数域磁异常场值为：

式中，

分别表示波数域磁异常场的x，z分量，

表示波数域磁位值，k表示波数。

8.一种基于NUFFT的二维磁异常快速正演模拟装置，其特征在于，所述装置包括：

二维磁性体模型构建模块，用于根据二维磁性体的大小，构建包含目标区域的二维磁性体模型；所述观测点坐标包括水平方向和垂向；所述二维磁性体模型是将所述目标区域在水平方向和垂向进行均匀或者非均匀网格剖分，并根据所述二维磁性体的截面形状设定网格剖分的每个节点的磁化率值；

空间域磁位积分表达式确定模块，用于根据所述二维磁性体模型，确定对应的空间域磁位积分表达式；所述空间域磁位积分表达式包含空间域磁化强度；所述空间域磁化强度是根据所述每个节点的磁化率值和已知的目标区域地球主磁场计算得到的；

波数域磁位一维积分表达式确定模块，用于对所述空间域磁位积分表达式沿水平方向进行一维傅里叶变换，得到波数域磁位一维积分表达式；其中，所述波数域磁位一维积分表达式包含波数域磁化强度；所述波数域磁化强度是对所述空间域磁化强度进行非均匀采样快速傅里叶变换得到的；

波数域磁位离散表达式确定模块，用于采用二次插值的形函数对所述波数域空间域垂向一维积分进行计算，得到波数域磁位离散表达式；

波数选取值确定模块，用于根据水平方向网格剖分的尺寸信息，确定截止频率，根据所述截止频率确定波数域波数选取范围，根据所述波数选取范围和预先设置的选取波数总数，得到波数选取值；

波数域磁位值确定模块，用于将所述波数选取值代入所述波数域磁位离散表达式，得到波数域磁位值；

波数域磁异常场值确定模块，用于根据傅里叶变换的微分性质，由所述波数域磁位值得到波数域磁异常场值；

空间域异常场值确定模块，用于通过对所述波数域磁异常场值进行一维非均匀快速傅里叶反变换，得到目空间域磁异常场值。

9.一种计算机设备，包括存储器和处理器，所述存储器存储有计算机程序，其特征在于，所述处理器执行所述计算机程序时实现权利要求1至7中任一项所述方法的步骤。

10.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现权利要求1至7中任一项所述的方法的步骤。

Images