CN108984939A - 基于3D Gauss-FFT的三维重力场正演方法 - Google Patents

Abstract

本发明提供一种基于3D Gauss‑FFT的三维重力场正演方法，包括以下步骤：步骤S100：三维场源体建模；步骤S200:3D Gauss‑FFT正变换计算三维重力场偏移频谱；步骤S300：3D Gauss‑FFT逆变换计算单个高斯节点的三维重力场；步骤S400:累加求和得到三维重力场分布。该方法首先提出3D Gauss‑FFT算法，即在每一个积分区域用高斯型求积公式代替传统的矩形求积公式，大大提高三维重力场计算精度，同时由于FFT算法内核的使用，计算效率得到保证，模型试验表明本方法正演精度与空间域解析解相当。

Description

基于3D Gauss-FFT的三维重力场正演方法

技术领域

本发明涉及勘探地球物理技术领域，具体的涉及一种基于3D Gauss-FFT的三维重力场正演方法。

背景技术

重力勘探作为一项最基本的地球物理方法，已被广泛运用于固体矿产资源勘探、岩石圈壳幔结构研究、地形校正、大地水准面测量、水文地质勘察、考古勘察等。测量得到的重力场数据与地下密度异常体有着密切关系，可被直接用于反演得到地下三维密度分布，从而推断地下构造运动、岩石圈结构构造等科学问题。正演的速度和精度决定了反演的可行性及可信度，因此在一个高效的反演程序中，快速高精度正演算法尤为重要。

在直角坐标系下，地下场源场被剖分为规则的直立棱柱体，传统的重力场正演方法有空间域解析解法，频率域FFT法。然而空间域解析解法计算零误差，但计算异常耗时；频率域FFT法计算效率高而计算精度难以保证。这一矛盾在三维重力场反演中尤为显著，不精确的正演结果以及低效率的正演方法将制约着反演的进行。在当前计算机硬件难以突破的背景下，想要大幅度提高计算效率，只能寻求新的正演方法。

发明内容

本发明的目的在于提供一种基于Gauss-FFT的三维重力场正演方法，该发明解决了现有直角坐标系下三维重力场正演算法正演效率低，计算精度差的技术问题。

本发明提供一种基于Gauss-FFT的三维重力场正演方法，包括以下步骤：

步骤S100：根据地下三维异常体的形态和尺寸，设置正演模型并剖分所述正演模型得到等体积的多个直立棱柱体小单元体，以任一所述小单元体的几何中心处为高斯节点，设定所述正演模型中高斯节点个数，每一个方向上高斯节点数为2，可以更多，查询[-1，1]上高斯节点系数表，将其转换为[0，1]上相应的节点值和系数值；

步骤S200：对第i个高斯节点采用3D Gauss-FFT正变换计算三维重力场偏移频谱：

步骤S210：分别设置X、Y、Z方向上第i个高斯节点和第i个高斯节点处的高斯系数数值对(λ_ix,α_ix),(λ_iy,α_iy),(λ_iz,α_iz)，其中，λ_ix为X方向上第i个高斯节点处的高斯系数数，λ_iy为Y方向上第i个高斯节点处的高斯系数数，λ_iz为Z方向上第i个高斯节点处的高斯系数数，α_ix为X方向上第i个高斯节点的值，α_iy为Y方向上第i个高斯节点的值，α_iz为Z方向上第i个高斯节点的值，离散密度ρ(x_m,y_n,z_l)乘以高斯偏移因子得到第i个高斯点上的偏移密度其中，e为自然底数；x_m为X方向上第m段的中心点坐标，y_n为Y方向上第n段的中心点坐标，z_l为Z方向上第l段的中心点坐标，α_ix为X方向上第i个高斯节点的值，α_iy为Y方向上第i个高斯节点的值，α_iz为Z方向上第i个高斯节点的值；Δk_x为k_x方向上频率的间隔，Δk_y为k_y方向上的频率间隔，Δk_z为k_z方向上频率的间隔。

步骤S220：对所述偏移密度进行三维离散傅里叶变换(3D DFT)，得到第i个高斯点上的密度分布的偏移频谱

步骤S230：对于所述密度分布的偏移频谱乘以大地滤波因子，得到第i个高斯点上的三维重力异常偏移频谱其中，k_xp为k_x方向上第p个点处的频率值，k_yq为k_y方向上第q个点处的频率值，k_zw为k_z方向上第w个点处的频率值；

步骤S300：对第i个高斯节点采用3D Gauss-FFT逆变换计算第i个高斯节点上的重力场响应；

步骤S400：取i＝i+1重复所述步骤S200～S300，直至i＝n时停止，其中n为总的高斯节点数，累加所得各高斯节点上的重力场响应得到所述地下三维异常体的三维重力场分布。

进一步地，所述步骤S100包括以下步骤：

步骤S110：根据地下三维异常体的形态和尺寸，设置用于包含所述三维异常体的棱柱体目标区域，所设置的棱柱体目标区域能将整个三维异常体完全嵌入其中；

步骤S120：设置所述三维异常体的剖分段数，确定所述三维异常体X、Y、Z方向上的剖分间隔，根据所述剖分段数和所述剖分间隔将所述目标区域剖分成体积相等的直立棱柱体小单元体，每一个小单元体的几何中心处为高斯节点，剖分段数的设置可根据问题的实际需要以及计算机的实际性能进行设置；

步骤S130：根据所述三维异常体的密度分布，对每一个剖分后的小直立棱柱体单元密度进行赋值，得到ρ(x_m,y_n,z_l)。

进一步地，所述步骤S300包括以下步骤：

步骤S310：对第i个高斯节点上的三维重力场偏移频谱做三维离散傅里叶逆变换(3D IDFT)，得到第i个高斯点处三维重力场偏移频谱空间域值；

步骤S320：将所述第i个高斯点处三维重力场偏移频谱空间域值乘以高斯逆偏移因子和高斯系数，得到第i个高斯节点上的重力场响应。

进一步地，所述步骤S220中所述三维离散傅里叶变换(3D DFT)包括以下步骤：

步骤S221：将排列好的偏移密度数据输入到Fortran库函数fftw3，程序将自动得到偏移密度的频谱。

进一步地，所述步骤S300中3D IDFT逆变换包括以下步骤：

步骤S310：将排列好的第i个点处所对应的三维重力异常频谱输入到Fortran库函数fftw3中，程序将自动计算其逆变换值并输出；

步骤S320：将程序输出的结果再乘以第i个高斯点所对应的高斯系数以及逆偏移因子，得到第i个高斯节点上的重力异常分布。

相对现有技术，本发明的技术效果：

本发明提供的基于Gauss-FFT的三维重力场正演方法，首先提出3D Gauss-FFT算法，即在每一个积分区域用高斯型求积公式代替传统的矩形求积公式，大大提高三维重力场计算精度，同时由于FFT算法内核的使用，计算效率得到保证，模型试验表明本方法正演精度与空间域解析解相当。

本发明提供的基于Gauss-FFT的三维重力场正演方法，核心内容为：针对传统基于FFT的三维重力场正演方法，提出用3D Gauss-FFT算法代替FFT算法，在每一个单元体上用更高精度的高斯型积分代替传统FFT算法中的矩形积分，有效减小了传统离散傅里叶变换中由于边界效应、强加周期、截断效应等所引起的诡源效应，大幅度提高了正演计算精度，同时由于FFT内核的使用，正演计算效率也得到了保证。模型试验表明，本发明所提出的方法使得三维重力场正演精度相比于传统FFT算法提高了近两个数量级，有效压制了诡源效应。本发明所提出的方法可直接用于三维重力场正演，井地重力场联合反演。

本发明提供的基于Gauss-FFT的三维重力场正演方法，在传统重力场频域正演中，使用快速傅里叶变换来提高计算效率。在现有技术中采用传统离散傅里叶变换算法时，在每一个积分区域中使用矩形求积公式，因此其正演精度较低，尤其是在截断边界处。本发明提出了3D Gauss-FFT方法，该方法在每一个积分区域中使用高斯型求积公式代替传统的矩形求积公式，使得计算精度大大提高。

本发明提供的基于Gauss-FFT的三维重力场正演方法，由于传统FFT内核的调用，本发明所提出的方法同样具有传统频域正演算法的高效性。模型算例表明，本发明所提出的方法能够有效减小诡源效应的影响，使得重力异常分布能够正确反映地下场源分布。诡源效应的消除对于三维重力异常反演至关重要，正确的异常分布才有可能得到正确的密度异常体分布。此外，与传统重力场频域算法相比，本发明所提出的方法在正演精度方面提高了两个数量级，在计算效率方法与传统方法相当，算例证明本发明所提出的方法是一种快速高精度的三维重力场正演方法，可直接用于三维重力场反演以及矿产水文勘查。

本发明提供的基于Gauss-FFT的三维重力场正演方法，在每一个小的积分区域中，运用高精度高斯型积分代替传统FFT算法中的矩形积分，大幅度提高直角坐标系下三维重力场正演精度，同时由于使用FFT内核，计算效率也得到了保证。

具体请参考根据本发明的基于Gauss-FFT的三维重力场正演方法提出的各种实施例的如下描述，将使得本发明的上述和其他方面显而易见。

附图说明

图1是本发明提供基于Gauss-FFT的三维重力场正演方法的流程示意图；

图2是本发明优选实施例地下三维场源网格剖分及观测点排列示意图；

图3是本发明优选实施例沿着y＝995m断面上的密度分布和各种方法计算的重力异常及其误差示意图；(a)为沿着该断面的密度分布；(b)为空间域解析解法计算结果；(c)为标准3D FFT算法计算结果；(d)为使用2点的3D Gauss-FFT算法正演结果；(e)为使用标准3D FFT算法计算结果相对于解析解法的绝对误差；(f)为使用3D Gauss-FFT算法正演结果相对于解析解法的绝对误差；

图4是本发明优选实施例中合成复杂密度模型及其理论垂直重力异常示意图；(a)沿着y＝2420m断面密度分布；(b)沿z＝-20m平面上理论垂直重力异常；

图5是本发明优选实施例合成复杂密度模型沿着z＝-20m平面所计算的重力异常及其绝对误差示意图；其中，(a)为采用标准2D FFT算法所得重力异常结果；(b)为采用标准3D FFT算法所得重力异常结果；(c)为采用4点的2D Gauss-FFT算法所得重力异常结果；(d)为采用4点的3D Gauss-FFT算法所得重力异常结果；(e)为采用标准2D FFT算法所得绝对误差结果；(f)为采用标准3D FFT算法所得绝对误差结果；(g)为采用4点的2D Gauss-FFT算法所得绝对误差结果；(h)为采用4点的3D Gauss-FFT算法所得绝对误差结果；

图例说明如下：

1、P(x,y,z)：表示观测点坐标；

2、(ξ_m,ζ_n,η_l)：表示编号为(m,n,l)的单元体的坐标位置；

3、Δx,Δy,Δz：表示场源在三个方向上的剖分间隔。x轴指北为正，y轴指东为正，z轴向下为正。

具体实施方式

构成本申请的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。

参见图1，本发明提供的基于Gauss-FFT的三维重力场正演方法，包括以下步骤：

步骤S100：根据地下三维异常体的形态和尺寸，设置正演模型并剖分所述正演模型得到等体积的多个直立棱柱体小单元体，以任一所述小单元体的几何中心处为高斯节点，设定所述正演模型中高斯节点个数，每一个方向上高斯节点数为2，可以更多，查询[-1，1]上高斯节点系数表，将其转换为[0，1]上相应的节点值和系数值；

步骤S200：对第i个高斯节点采用3D Gauss-FFT正变换计算三维重力场偏移频谱：

步骤S210：分别设置X、Y、Z方向上第i个高斯节点和第i个高斯节点处的高斯系数数值对(λ_ix,α_ix),(λ_iy,α_iy),(λ_iz,α_iz)，其中，λ_ix为X方向上第i个高斯节点处的高斯系数数，λ_iy为Y方向上第i个高斯节点处的高斯系数数，λ_iz为Z方向上第i个高斯节点处的高斯系数数，α_ix为X方向上第i个高斯节点的值，α_iy为Y方向上第i个高斯节点的值，α_iz为Z方向上第i个高斯节点的值，离散密度ρ(x_m,y_n,z_l)乘以高斯偏移因子得到第i个高斯点上的偏移密度

步骤S220：对所述偏移密度进行三维离散傅里叶变换(3D DFT)，得到第i个高斯点上的密度分布的偏移频谱

步骤S230：对于所述密度分布的偏移频谱乘以大地滤波因子，得到第i个高斯点上的三维重力异常偏移频谱

步骤S300：对第i个高斯节点采用3D Gauss-FFT逆变换计算第i个高斯节点上的重力场响应；

步骤S400：取i＝i+1重复所述步骤S200～S300，直至i＝n时停止，其中n为总的高斯节点数，累加所得各高斯节点上的重力场响应得到所述地下三维异常体的三维重力场分布。

上述步骤S300中得到的仅为单个高斯节点处的重力异常分布，如果每个方向上都运用2点高斯积分，那么需要重复上述步骤S300～S400 8次(2×2×2)，将每次计算结果乘以相应的系数，并累加求和，得到最终的三维重力异常分布。

优选的，所述步骤S100包括以下步骤：

步骤S110：根据地下三维异常体的形态和尺寸，设置用于包含所述三维异常体的棱柱体目标区域。所设置的棱柱体目标区域能将整个三维异常体完全嵌入其中。

步骤S120：设置所述三维异常体的剖分段数，确定所述三维异常体X、Y、Z方向上的剖分间隔，根据所述剖分段数和所述剖分间隔将所述目标区域剖分成体积相等的直立棱柱体小单元体，每一个小单元体的几何中心处为高斯节点。剖分段数的设置可根据问题的实际需要以及计算机的实际性能进行设置。

步骤S130：根据所述三维异常体的密度分布，对每一个剖分后的小直立棱柱体单元密度进行赋值，得到ρ(x_m,y_n,z_l)，其中，x_m为X方向上第m段的中心点坐标，y_n为Y方向上第n段的中心点坐标，z_l为Z方向上第l段的中心点坐标。

优选的，所述步骤S300包括以下步骤：

步骤S310：对第i个高斯节点上的三维重力场偏移频谱做三维离散傅里叶逆变换(3D IDFT)，得到第i个高斯点处三维重力场偏移频谱空间域值；

步骤S320：将所述第i个高斯点处三维重力场偏移频谱空间域值乘以高斯逆偏移因子和高斯系数，得到第i个高斯节点上的重力场响应。

优选的，所述步骤S220中所述三维离散傅里叶变换(3D DFT)包括以下步骤：

步骤S221：将排列好的偏移密度数据输入到Fortran库函数fftw3，并给定每一个方向上的数据点个数、输出的三维数组名称、以及正变换标志变量值。

程序将自动调用快速傅里叶变换蛇形算法得到偏移密度的频谱，而用户无需参与具体内部操作。

优选的，所述步骤S300中3D IDFT逆变换包括以下步骤：

步骤S310：将排列好的第i个点处所对应的三维重力异常频谱输入到Fortran库函数fftw3中，程序将自动计算其逆变换值并输出。

步骤S320：将程序输出的结果再乘以第i个高斯点所对应的高斯系数以及逆偏移因子，得到第i个高斯节点上的重力异常分布。

具体的，该方法中三维重力异常频谱表达式的推导如下：

地下任意分布的异常体都可以剖分成如图2所示的规则的直立棱柱体，将整个场源剖分成M×N×L个规则的单元体，对于编号为(m,n,l)的直立棱柱体而言，其中m＝1,2…,M；n＝1,2…,N；l＝1,2…,L，其几何中心位置为(ξ_m,ζ_n,η_l)，该直立棱柱体的剩余密度ρ(ξ_m,ζ_n,η_l)看成常数。对于三维频域重力异常正演来说，观测点P(x,y,z)的坐标与网格剖分一致。而对于二维频域重力场正演来说，观测点位于一个固定高度的水平面上(如图2所示)。

在空间域，由一个直立棱柱体在任意观测点P(x,y,z)所产生的重力位U可以写为一个三维积分公式：

其中，G＝6.67×10^-11N·m²/kg²是重力常数，是观测点到场源(ξ,ζ,η)之间的距离。

对上述三维重力位公式做三维傅里叶变换，得到包含于场源体内的任意观测点处的重力位频谱响应为：

其中，k_x,k_y,k_z是坐标系直角坐标系下x、y、z方向上所对应的频率域的波数。

上面(1)和(2)式中的格林函数的三维频谱可以写为

其中为三个方向上的总波数,，π为圆周率

将(3)式带入(2)式，可得：

重力异常的频谱可以通过对上式(4)求垂向导数得到：

使用重力异常的叠加性质，整个三维场源所产生的重力异常的频谱可以通过累加M×N×L个独立直立棱柱体的频谱得到：

其中：

可以看出，在公式(6)右侧的左成因子类似于大地滤波因子，这就表明重力异常可以看做是三维密度分布在大地滤波之后的版本。

那么整个三维场源内部的观测点上的重力异常Δg可以通过三维傅里叶逆变换得到：

至此，便可以利用三维傅里叶变换的方法求取三维重力异常的频谱，然后再实施一次逆变换，就可以得到三维重力异常分布。

第二点：3DGauss-FFT算法的应用。

在实际中处理的都是离散形式的数据，这就意味着上述公式(7)和(8)中的连续表达式要被相应的离散形式所取代。公式中的连续傅里叶变换将被离散傅里叶变换所取代，在这个过程中，不可避免地带来了截断误差，强加周期，边界效应等一系列问题，使得重力场的正演精度大大降低，尤其是在边界处。因此本发明提出3D Gauss-FFT算法，用高斯型求积公式代替传统的矩形积分，大大提高正演精度。

为了方便计算，以Δx,Δy,Δz为剖分间隔均匀采样，采样点数为M×N×L。每一个点的位置可以表示为

当x₀＝0.5Δx,y₀＝0.5Δy,z₀＝0.5Δz时，采样点便和剖分单元几何中心点重合。假设采样点M,N,L都是偶数，那么常见的离散波数的取法为：

那么离散密度数据ρ(x_m,y_n,z_l)的三维离散傅里叶变换以及其逆变换分别可以写为：

公式(11)和(12)是三维连续傅里叶变换的矩形积分近似形式。注意到公式(7)和公式(11)具有相同的形式，因此在实际中，公式(11)被用来计算三维密度分布的频谱。将所得到的密度分布频谱再乘以大地滤波因子之后便得到重力异常的频谱最后对重力异常频谱做三维傅里叶逆变换就得到三维重力异常Δg(x,y,z)。

为了提高逆傅里叶变换的精度，在每一个积分区域使用高斯型数值积分代替传统的矩形积分。常用的一维高斯积分公式形如：

其中K是高斯节点数；A_k和t_k是[-1,1]上的高斯节点和高斯系数。

这里三维高斯积分的积分点个数设为I_x,I_y,I_z，[0,1]上高斯系数和高斯节点分别为(λ_ix,α_ix),(λ_iy,α_iy)和(λ_iz,α_iz),ix＝1,2…I_x；iy＝1,2…I_y；iz＝1,2…I_z.首先将(11)式中频谱在k_x,k_y,k_z方向上偏移α_ix,α_iy,α_iz，得到偏移的频谱

重新整理得：

由公式(15)可知，偏移频谱可以通过对ρ(x_m,y_n,z_l)乘以高斯偏移因子然后做三维离散傅里叶变换得到。

将上述偏移频谱乘以相应的偏移大地滤波因子这里得到第i个高斯点处的重力异常分布

基于高斯积分公式，3D Gauss-IDFT可以通过对偏移频谱的加权求和得到：

其中，加权系数为λ_ixλ_iyλ_iz。重新整理得：

当把偏移频谱看做是另一组新的离散数据，那么公式(18)中的A_ix,iy,iz与公式(12)中的3DIDFT有相同的形式。因此，可以直接运用快速傅里叶变换(FFT)算法来提高计算效率。此外，在实际中由于很难得到高斯节点处的空间域数据，因此在正变换中仍然使用矩形积分公式。

将上述的3D Gauss-IDFT方法运用到三维重力异常的计算中，得到本发明提供的方法。该方法基本包括以下步骤：

1)离散密度数据ρ(x_m,y_n,z_l),选择高斯节点和高斯系数(λ_ix,α_ix),(λ_iy,α_iy),(λ_iz,α_iz)；

2)乘以高斯偏移因子:

3)通过传统3D DFT计算密度分布的偏移频谱

4)乘以大地滤波因子(公式6)得到:

5)通过3D IDFT(公式18)计算A_ix,iy,iz；

6)乘以逆高斯偏移因子和高斯系数:

7)累加求和，遍历I_x,I_y,I_z(公式17),得到:Δg(x_m,y_n,z_l)。

上述步骤基于3DGauss-FFT算法的三维重力异常正演方法流程同样适用于二维的情况。

注意到上述方法计算的是三维重力异常，作为选择，二维傅里叶变换方法可以用来计算某一z₀平面上的二维重力异常。一个由M×N个直立棱柱体组成的质量层(l)所产生的重力异常频谱可以表示为：

其中，

那么地下所有层的三维密度分布的重力异常频谱可以写成累加求和的形式

其中，在深度z方向上的剖分可以是任意的。

对上式(21)做二维傅里叶逆变换可得空间域重力异常分布Δg(x,y,z₀)

其中水平观测平面z₀必须小于场源，也就是说水平观测面必须在场源正上方。

与三维时相似，二维逆Gauss-FFT可以写为：

其中

以下结合具体实施例，对本发明提供方法进行详细说明。

下面将给出两个模型算例以佐证本发明所提出的方法的有效性。首先设计一个简单三维密度模型，用于测试本发明方法对诡源效应消除方面的作用。该三维密度模型在x,y,z三个方向上的起始范围为：0-2000m,0-2000m,0-1000m，分别剖分为200×200×100个单元体。场源是一个立方体，密度差为1000kg/m³，在x方向上范围为800-1200m，在y方向上范围为600-1400m，在z方向上范围为700-900m；分别剖分为40×80×20个小单元体。图3(a)是y＝995m断面上的密度分布。

为了展示本发明所提出的方法的计算精度，计算整个场源区域的三维重力异常。图3(b)-(d)展示了y＝995m断面上的理论重力异常、用传统3D FFT算法计算的重力异常以及适用3D Gauss-FFT算法所计算的重力异常。图3(e)和图3(f)分别为相应的绝对误差。

如图3(b)所示，垂直断面上的理论重力异常是关于异常体中心对称分布的，也就是说零重力异常等值线位于z＝800m，改线刚好穿过异常体的中心，正的重力异常分布在该剖面的上方，而负的重力异常分布在其下方。使用传统3D FFT方法正演结果也具有相似的异常分布特征，但是在研究区域的上边界处出现了明显的假的负异常，这种加的异常称之为诡源，这种效应称之为诡源效应。使用传统3D FFT算法正演结果的绝对误差高达4mGal(图3e)，诡源的出现严重影响真实场源，在反演计算中将会导致虚假异常的产生，给资料解释带来麻烦。相反，使用本发明所提出的3D Gauss-FFT算法正演结果和理论结果几乎完全一致，其相应的最大绝对误差仅为0.06mGal。更为重要的是，本发明所提出的方法不会产生诡源。

为了更进一步去对比本发明所提出的方法的性能，并对比现有方法与本发明提供方法进行对比，接下来设计一个复杂密度异常体，后续实施例1～3和对比例1～6均对该模型进行处理。

该模型包括一个倾斜台阶和一个浅层立方体。研究区域在x,y,z三个方向上的范围都是0-4800m，剖分成120×120×120个等间隔的直立棱柱体单元。倾斜台阶模型和浅层立方体模型的密度差分别设为：2000kg/m³和-2000kg/m³。观测点位于z＝-20m的一个水平面上，测点个数为120×120。图4(a)展示了y＝2420m断面上的密度异常分布；图4(b)为沿z＝-20m平面上理论垂直重力异常。由图4(a)和(b)可见，重力异常的分布形态与异常体模型密度分布有着直接对应关系，也即图4(b)中的负重力异常值对应于图4(a)中的负密度异常，而4(b)中的正的重力异常值是倾斜台阶模型的直接反应。

作为对比，本算例同样包含了二维频域正演方法。如前所述，这种方法仅仅适用于计算高于场源体的某一水平面上的重力场，因此在这里选取观测面高度为z＝-20m。图5展示了使用不同方法(2D FFT法，2D Gauss-FFT法，3D FFT法，3D Gauss-FFT法)在该观测平面上的重力异常值。

对比例1

采用空间域方法处理上述模型。

对比例2

采用标准2D FFT方法处理上述模型。

对比例3

采用标准3D FFT方法处理上述模型。

对比例4

采用标准2点2D Gauss-FFT方法处理上述模型。

对比例5

采用标准4点2D Gauss-FFT方法处理上述模型。

对比例6

采用标准8点2D Gauss-FFT方法处理上述模型。

实施例1

采用本发明提供的方法进行2点3D Gauss-FFT方法处理上述模型。

实施例2

采用本发明提供的方法进行4点3D Gauss-FFT方法处理上述模型。

实施例3

采用本发明提供的方法进行8点3D Gauss-FFT方法处理上述模型。

实施例1～3和对比例1～6的计算时间、最大绝对误差和所需内存结果，列于表1中。所得正演结果如图5中的(a)～(h)所示。

图5(a)和图5(b)可以看出，使用传统2D和3D FFT方法所得到的正演结果明显被低估。与空间域解析解对比，传统频域正演方法所得结果的重力异常分布有较大畸变，尤其是传统3D FFT法。从图5(e)和图5(f)可以看出，其最大绝对误差分别为4.5mGal和9mGal。

图5(c)和图5(d)展示了使用4节点的2D Gauss-FFT方法和4节点的3D Gauss-FFT方法的正演结果，其正演的重力异常分布特征与空间域解析解法完全相同。图5(g)和图5(h)分别展示了其相对于空间域解析解的绝对误差分布，可以看出，最大绝对误差一般分布于密度场源变换最剧烈的边界处。相比于标准FFT方法(图5(e)和图5(f))，发现使用Gauss-FFT方法的正演精度提高了近两个数量级。

表1合成复杂密度模型的重力场正演计算性能统计对比。

方法	计算时间(s)	最大绝对误差(mGal)	所需内存(Mb)
				空间域方法	11857.8	-	13.82
标准2D FFT	0.3	4.52	27.64
				标准3D FFT	0.3	8.77	27.64
2点2D Gauss-FFT	3.7	1.396	27.64
				4点2D Gauss-FFT	14.7	0.014	27.64
8点2D Gauss-FFT	57.8	0.009	27.64
				2点3D Gauss-FFT	10.4	1.171	82.92
4点3D Gauss-FFT	81.2	0.031	82.92
				8点3D Gauss-FFT	652.6	0.030	82.92

表1列出了不同正演方法在计算时间、最大绝对误差和内存需求方面的综合比较。结果表明空间域正演方法精度最高，但同时也是最耗时的。相反标准2D和3D FFT正演方法需要更少的时间但是经度更低，在计算时间和计算精度之间最佳的折中便是使用Gauss-FFT方法。

本算例还对比了不同的高斯节点数对正演综合性能的影响，如表1所示。在Gauss-FFT方法中，随着高斯节点数的增加，正演精度随着增加，但是更多的节点数意味着需要更多的计算时间。对于使用8个高斯节点的计算时间是使用4个高斯节点计算时间的数倍，而在经度上的提高却是微乎其微的。因此，可以得出结论，使用4点的Gauss-FFT方法足以得到勘探尺度中重力异常正反演的要求。其次，使用3D Gauss-FFT方法的计算时间比使用相同节点数的2D Gauss-FFT方法的计算时间更长，这是由于前者得到的是整个场源中的三维重力异常，而后者的2D Gauss-FFT方法仅仅得到的是某一平面上的重力异常。因此可以得出结论：如果只想得到某一平面上的二维重力异常，4点的2D Gauss-FFT方法具有不可比拟的优势，但是如果想计算整个场源的三维重力异常分布，那么3D Gauss-FFT法将会更为有效。

对比2D和3D Gauss-FFT方法，发现2D Gauss-FFT方法具有更高的正演精度，这主要是因为在2D Gauss-FFT(或者FFT)法中，在深度方向上的积分是精确的，而在3D Gauss-FFT(或者FFT)中只能通过数值积分(矩形积分或者高斯型积分)得到，这一步将不可避免地引入误差。

关于在内存消耗方面，空间域方法所需内存最小，因为在空间域正演中，仅仅需要存储三维密度分布矩阵和一个二维观测点矩阵。使用2D FFT(或者Gauss-FFT)方法所占用的内存是空间域方法的二倍，这是因为在频域，所有的矩阵都定义为复数型。由于需要存储三维密度分布以及相应的三维观测点数据，3D Gauss-FFT(或者FFT)法正演需要最大的内存消耗。

本领域技术人员将清楚本发明的范围不限制于以上讨论的示例，有可能对其进行若干改变和修改，而不脱离所附权利要求书限定的本发明的范围。尽管己经在附图和说明书中详细图示和描述了本发明，但这样的说明和描述仅是说明或示意性的，而非限制性的。本发明并不限于所公开的实施例。

通过对附图，说明书和权利要求书的研究，在实施本发明时本领域技术人员可以理解和实现所公开的实施例的变形。在权利要求书中，术语“包括”不排除其他步骤或元素，而不定冠词“一个”或“一种”不排除多个。在彼此不同的从属权利要求中引用的某些措施的事实不意味着这些措施的组合不能被有利地使用。权利要求书中的任何参考标记不构成对本发明的范围的限制。

Claims (5)

1.一种基于3D Gauss-FFT的三维重力场正演方法，其特征在于，包括以下步骤：

步骤S100：根据地下三维异常体的形态和尺寸，设置正演模型并剖分所述正演模型得到等体积的多个直立棱柱体小单元体，以任一所述小单元体的几何中心处为高斯节点，设定所述正演模型中高斯节点个数，每一个方向上高斯节点数为2，可以更多，查询[-1，1]上高斯节点系数表，将其转换为[0，1]上相应的节点值和系数值；

步骤S200：对第i个高斯节点采用3D Gauss-FFT正变换计算三维重力场偏移频谱：

步骤S210：分别设置X、Y、Z方向上第i个高斯节点和第i个高斯节点处的高斯系数数值对(λ_ix,α_ix),(λ_iy,α_iy),(λ_iz,α_iz)，其中，λ_ix为X方向上第i个高斯节点处的高斯系数数，λ_iy为Y方向上第i个高斯节点处的高斯系数数，λ_iz为Z方向上第i个高斯节点处的高斯系数数，α_ix为X方向上第i个高斯节点的值，α_iy为Y方向上第i个高斯节点的值，α_iz为Z方向上第i个高斯节点的值，离散密度ρ(x_m,y_n,z_l)乘以高斯偏移因子得到第i个高斯点上的偏移密度其中j为虚数单位；

步骤S220：对所述偏移密度进行三维离散傅里叶变换(3D DFT)，得到第i个高斯点上的密度分布的偏移频谱

步骤S230：对于所述密度分布的偏移频谱乘以大地滤波因子，得到第i个高斯点上的三维重力异常偏移频谱

步骤S300：对第i个高斯节点采用3D Gauss-FFT逆变换计算第i个高斯节点上的重力场响应；

步骤S400：取i＝i+1重复所述步骤S200～S300，直至i＝n时停止，其中n为总的高斯节点数，累加所得各高斯节点上的重力场响应得到所述地下三维异常体的三维重力场分布。

2.根据权利要求1所述的基于3D Gauss-FFT的三维重力场正演方法，其特征在于，所述步骤S100包括以下步骤：

步骤S110：根据地下三维异常体的形态和尺寸，设置用于包含所述三维异常体的棱柱体目标区域，所设置的棱柱体目标区域能将整个三维异常体完全嵌入其中；

步骤S120：设置所述三维异常体的剖分段数，确定所述三维异常体X、Y、Z方向上的剖分间隔，根据所述剖分段数和所述剖分间隔将所述目标区域剖分成体积相等的直立棱柱体小单元体，每一个小单元体的几何中心处为高斯节点，剖分段数的设置可根据问题的实际需要以及计算机的实际性能进行设置；

步骤S130：根据所述三维异常体的密度分布，对每一个剖分后的小直立棱柱体单元密度进行赋值，得到ρ(x_m,y_n,z_l)，其中，x_m为X方向上第m段的中心点坐标，y_n为Y方向上第n段的中心点坐标，z_l为Z方向上第l段的中心点坐标。

3.根据权利要求2所述的基于3DGauss-FFT的三维重力场正演方法，其特征在于，所述步骤S300包括以下步骤：

步骤S310：对第i个高斯节点上的三维重力场偏移频谱做三维离散傅里叶逆变换，得到第i个高斯点处三维重力场偏移频谱空间域值；

步骤S320：将所述第i个高斯点处三维重力场偏移频谱空间域值乘以高斯逆偏移因子和高斯系数，得到第i个高斯节点上的重力场响应。

4.根据权利要求3所述的基于Gauss-FFT的三维重力场正演方法，其特征在于，所述步骤S220中所述三维离散傅里叶变换(3D DFT)包括以下步骤：

步骤S221：将排列好的偏移密度数据输入到Fortran库函数fftw3，程序将自动得到偏移密度的频谱。

5.根据权利要求4所述的基于3D Gauss-FFT的三维重力场正演方法，其特征在于，所述步骤S300中3D IDFT逆变换包括以下步骤：

步骤S310：将排列好的第i个点处所对应的三维重力异常频谱输入到Fortran库函数fftw3中，程序将自动计算其逆变换值并输出；

步骤S320：将程序输出的结果再乘以第i个高斯点所对应的高斯系数以及逆偏移因子，得到第i个高斯节点上的重力异常分布。