二度体重力异常积分解数值模拟方法和装置与流程

本申请涉及二度体重力异常计算技术领域，特别是涉及一种二度体重力异常积分解数值模拟方法、装置、计算机设备和存储介质。

背景技术：

重力勘探是一种以地下介质密度差异作为物性基础的一种地球物理勘探方法。该方法施工方便、成本低、效率高、勘探深度大，已广泛应用于地球深部构造研究、大地及区域地质构造研究、构造单元划分、隐伏岩体或岩层及断裂的探测、成矿远景区的划分、石油、天然气或煤等远景盆地圈定等方面。在实际生产中，常把走向方向尺度远比垂直走向方向尺度大的地质体用走向方向无限延伸的二度体代替。针对走向方向无限延伸的二度体模型，高效、高精度数值模拟作为反演成像与定量解释的基础，一直是研究的焦点。

目前，有关二度体的数值模拟方法主要分为空间域方法和波数域方法。空间域方法是通过异常场解析式直接求取空间任意点的精确异常场，其优点是原理简单、结果精确，缺点是解析式比较复杂，推导过程繁琐，且当计算大量位置点异常场数据时，速度较慢。波数域方法是根据场源产生异常场的傅里叶变换解析表达式，通过数值方法计算该异常频谱的反傅里叶变换得到空间域异常场。相对于空间域方法而言，其优点是表达式简洁，计算效率较高。随着快速傅里叶变换扩边法和高斯快速傅里叶变换法gauss-fft在重力异常正演数值模拟中的广泛应用，频率域方法以其简洁性、准确性和高效性成为处理大规模复杂模型正演的首选方法。文献(tontini,f.c.,l.cocchi,c.carmisciano.rapid3-dforwardmodelofpotentialfieldswithapplicationtothepalinuroseamountmagneticanomaly(southerntyrrheniansea,italy).journalofgeophysicalresearch,2009.114.)利用傅里叶变换，详细推导了全空间的波数域表达式，借助快速傅里叶变换算法，实现了重力异常快速数值模拟，但在空间域和波数域只能均匀采样，并且由于截断边界效应的影响，采用扩边的方式使得数值模拟精度相对较低)；文献(wu,l.,tian,g.high-precisionfourierforwardmodelingofpotentialfields.geophysics,2014,79(5):g59-g68.)提出了一种重力异常正演模拟的gauss-fft方法，该方法有效克服了快速傅里叶变换方法的边界截断效应问题，提高了数值模拟精度，但在空间域和频率域只能均匀采样，并且增加了成倍的计算量，效率相对较低。文献(leejy,greengardl.thetype3nonuniformfftanditsapplications.journalofcomputationalphysics,2005,206(1):1-5.)将非均匀快速傅里叶变换(nufft)应用于核磁共振模拟成像中，取得了很好的效果。

目前，不管是传统快速傅里叶变换扩边法还是gauss-fft法，在空间域和波数域只能均匀采样，只适用于规则测线，且难以兼顾计算精度与计算速度。

技术实现要素：

基于此，有必要针对上述技术问题，提供一种能够兼顾计算精度与计算速度的二度体重力异常积分解数值模拟方法、装置、计算机设备和存储介质。

一种二度体重力异常积分解数值模拟方法，所述方法包括：

根据观测点在笛卡尔坐标系下的位置信息，构建包含目标区域的二度体模型；所述观测点坐标包括x轴方向和z轴方向；所述二度体模型是将所述目标区域分成若干个矩形，所述二度体的截面形状用于确定每个矩形的空间域异常体剩余密度值；

根据所述二度体模型，确定对应的空间域重力异常表达式；所述空间域重力异常表达式包含空间域异常体剩余密度值；

对所述空间域重力异常表达式沿x轴方向进行一维傅里叶变换，得到波数域重力异常表达式；其中，所述波数域重力异常表达式包含波数域异常体剩余密度值；所述波数域异常体剩余密度值是对所述空间域异常体剩余密度值进行非均匀采样快速傅里叶变换得到的；

根据所述矩形沿x轴方向的尺寸信息，确定截止频率，根据所述截止频率确定波数域波数采样值范围，根据所述波数采样值范围得到波数采样值；

将所述波数采样值代入所述波数域重力异常表达式，得到波数域重力异常场值；

通过对所述波数域重力异常场值进行一维非均匀快速傅里叶反变换，得到目标区域内任一点的空间域重力异常场值。

在其中一个实施例中，还包括：空间域重力异常表达式为：

式中：表示空间域重力异常场；表示万有引力常数；表示x方向剖分的矩形个数；表示z方向剖分的矩形个数；表示所述观测点坐标；表示编号为的矩形中心坐标；表示编号为的矩形的空间域异常体剩余密度值，表示矩形的x轴方向尺寸，矩形的z轴方向尺寸。

在其中一个实施例中，还包括：波数域重力异常表达式为：

式中：表示波数域重力异常场；k表示波数；表示波数域异常体剩余密度值；表示符号函数：

。

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

式中：i为虚数单位，为给定离散采样点对应采样点的值，为计算的离散采样点傅里叶变换频谱，n表示采样点总数；

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

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

式中：m表示过采样因子，表示权重因子，为精度因子，表示对取整；

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

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

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

在其中一个实施例中，还包括：将所述目标区域分成若干个矩形，每个矩形的x轴方向矩形尺寸、z轴方向矩形尺寸可以相同，也可以不同；在重力异常场变化快的地方加密网格，在重力异常场变化慢的地方稀疏网格。

在其中一个实施例中，还包括：根据所述矩形沿x轴方向的尺寸信息，确定截止频率为：

式中，表示截止频率；表示x轴方向矩形尺寸的最小值。

在其中一个实施例中，还包括：

根据所述截止频率确定波数域波数采样值范围为；

根据所述波数采样值范围和预先设置的采样点总数，在波数域进行非均匀采样，得到波数采样值，包括：

在上波数依次取值为：

在上波数依次取值为：

式中，，n表示预先设置的采样点总数。

一种二度体重力异常积分解数值模拟装置，所述装置包括：

二度体模型构建模块，用于根据观测点在观测点坐标中的位置信息，构建包含目标区域的二度体模型；所述观测点坐标包括x轴方向和z轴方向；所述二度体模型是将所述目标区域分成若干个矩形，所述二度体的截面形状用于确定每个矩形的空间域异常体剩余密度值；

空间域重力异常表达式确定模块，用于根据所述二度体模型，确定对应的空间域重力异常表达式；所述空间域重力异常表达式包含空间域异常体剩余密度值；

波数域重力异常表达式确定模块，用于对所述空间域重力异常表达式沿x轴方向进行一维傅里叶变换，得到波数域重力异常表达式；其中，所述波数域重力异常表达式包含波数域异常体剩余密度值；所述波数域异常体剩余密度值是对所述空间域异常体剩余密度值进行非均匀采样快速傅里叶变换得到的；

波数采样值确定模块，用于根据所述矩形沿x轴方向的尺寸信息，确定截止频率，根据所述截止频率确定波数域波数采样值范围，根据所述波数采样值范围得到波数采样值；

波数域重力异常场值确定模块，用于将所述波数采样值代入所述波数域重力异常表达式，得到波数域重力异常场值；

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

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

根据观测点在笛卡尔坐标系下的位置信息，构建包含目标区域的二度体模型；所述观测点坐标包括x轴方向和z轴方向；所述二度体模型是将所述目标区域分成若干个矩形，所述二度体的截面形状用于确定每个矩形的空间域异常体剩余密度值；

根据所述二度体模型，确定对应的空间域重力异常表达式；所述空间域重力异常表达式包含空间域异常体剩余密度值；

对所述空间域重力异常表达式沿x轴方向进行一维傅里叶变换，得到波数域重力异常表达式；其中，所述波数域重力异常表达式包含波数域异常体剩余密度值；所述波数域异常体剩余密度值是对所述空间域异常体剩余密度值进行非均匀采样快速傅里叶变换得到的；

根据所述矩形沿x轴方向的尺寸信息，确定截止频率，根据所述截止频率确定波数域波数采样值范围，根据所述波数采样值范围和预先设置的采样点总数，在波数域进行非均匀采样，得到波数采样值；

将所述波数采样值代入所述波数域重力异常表达式，得到波数域重力异常场值；

通过对所述波数域重力异常场值进行一维非均匀快速傅里叶反变换，得到目标区域内任一点的空间域重力异常场值。

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

根据观测点在笛卡尔坐标系下的位置信息，构建包含目标区域的二度体模型；所述观测点坐标包括x轴方向和z轴方向；所述二度体模型是将所述目标区域分成若干个矩形，所述二度体的截面形状用于确定每个矩形的空间域异常体剩余密度值；

根据所述二度体模型，确定对应的空间域重力异常表达式；所述空间域重力异常表达式包含空间域异常体剩余密度值；

对所述空间域重力异常表达式沿x轴方向进行一维傅里叶变换，得到波数域重力异常表达式；其中，所述波数域重力异常表达式包含波数域异常体剩余密度值；所述波数域异常体剩余密度值是对所述空间域异常体剩余密度值进行非均匀采样快速傅里叶变换得到的；

根据所述矩形沿x轴方向的尺寸信息，确定截止频率，根据所述截止频率确定波数域波数采样值范围，根据所述波数采样值范围和预先设置的采样点总数，在波数域进行非均匀采样，得到波数采样值；

将所述波数采样值代入所述波数域重力异常表达式，得到波数域重力异常场值；

通过对所述波数域重力异常场值进行一维非均匀快速傅里叶反变换，得到目标区域内任一点的空间域重力异常场值。

上述二度体重力异常积分解数值模拟方法、装置、计算机设备和存储介质，通过构建包含目标区域的二度体模型，确定对应的空间域重力异常表达式，其中空间域重力异常表达式包含空间域异常体剩余密度值，通过对空间域重力异常表达式沿x轴方向进行一维傅里叶变换，得到波数域重力异常表达式，其中波数域异常体剩余密度值是对空间域异常体剩余密度值进行非均匀采样快速傅里叶变换得到的，根据矩形沿x轴方向的尺寸信息，确定截止频率，根据截止频率确定波数域波数采样值范围，根据波数采样值范围和预先设置的采样点总数进行波数域非均匀采样，得到波数采样值，根据波数采样值，代入波数域重力异常表达式中，得到波数域重力异常场值，最后通过对波数域重力异常场值进行一维非均匀快速傅里叶反变换，得到目标区域内任一点的空间域重力异常场值。采用本方法进行二度体重力异常积分解数值模拟，通过借助非均匀采样快速傅里叶变换，能够实现空间域和波数域非均匀采样，可以有效克服快速傅里叶变换边界效应导致计算精度低，又能保证精度的前提下解决gauss-fft方法随着高斯点数增加计算效率低的问题。有效兼顾了重力异常数值模拟的计算精度与计算效率，解决了波数域重力异常数值模拟只适合规则测线、数值模拟方法不能同时兼顾计算效率和计算精度的问题。

附图说明

图1为一个实施例中二度体重力异常积分解数值模拟方法的流程示意图；

图2为一个实施例中二度体重力异常积分解数值模拟方法的算法流程图；

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

图4为一个实施例中截面为圆形二度体模型示意图；

图5为一个实施例中波数域非均匀采样模拟结果与理论值对比图；

图6为一个实施例中波数域非均匀采样模拟结果与理论值的相对误差；

图7为另一个实施例中空间域和波数域都非均匀采样模拟结果与理论值对比图；

图8为另一个实施例中空间域和波数域都非均匀采样模拟结果与理论值的相对误差；

图9为另一个实施例中二度体重力异常积分解数值模拟装置的结构框图；

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

具体实施方式

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

本申请提供的二度体重力异常积分解数值模拟方法，可以应用于如下应用环境中。根据观测点构建包含目标区域的二度体模型，确定对应的空间域重力异常表达式，其中空间域重力异常表达式包含空间域异常体剩余密度值，通过对空间域重力异常表达式沿x轴方向进行一维傅里叶变换，得到波数域重力异常表达式，其中波数域异常体剩余密度值是对空间域异常体剩余密度值进行非均匀采样快速傅里叶变换得到的，根据矩形沿x轴方向的尺寸信息，确定截止频率，根据截止频率确定波数域波数采样值范围，根据波数采样值范围和预先设置的采样点总数进行波数域非均匀采样，得到波数采样值，根据波数采样值，代入波数域重力异常表达式中，得到波数域重力异常场值，最后通过对波数域重力异常场值进行一维非均匀快速傅里叶反变换，得到目标区域内任一点的空间域重力异常场值。

在一个实施例中，如图1所示，提供了一种二度体重力异常积分解数值模拟方法，包括以下步骤：

步骤102，根据观测点在观测点坐标中的位置信息，构建包含目标区域的二度体模型。

观测点坐标包括x轴方向和z轴方向；二度体模型是将目标区域分成若干个矩形，矩形的截面形状用于确定每个矩形的空间域异常体剩余密度值。

具有一定走向，且沿走向方向异常值变化不明显的异常，称为二度异常，它是二维空间坐标的函数，引起二度异常的地质体称为二度体，如层状矿体、脉状矿体、向斜、背斜等。

构建包含目标区域的二度体模型，将目标区域分成若干个矩形，每个小矩形的几何尺寸,可以相同，也可以不同，剖分矩形的网格剖分灵活，可以在重力异常场变化快的地方加密网格，在重力异常场变化慢的地方适当稀疏网格。根据二度体的截面形状，对矩形密度进行赋值，每个矩形的密度是一个常数，不同矩形密度值不同，以此表示任意截面形状、任意密度分布的二度体模型。

步骤104，根据二度体模型，确定对应的空间域重力异常表达式；空间域重力异常表达式包含空间域异常体剩余密度值。

由于实际地球内部的物质密度分布非常不均匀，因而实际观测重力值与理论上的正常重力值总是存在着偏差，这种在排除各种干扰因素影响之后，仅仅是由于物质密度分布不匀而引起的重力的变化，就称为重力异常。重力场正演是指根据密度分布计算重力场，反演是指根据观测重力值计算密度分布。正演是反演的基础，正演计算的效率直接影响反演计算的效率，而正演计算精度直接影响反演结果的质量。在二度体重力异常积分解数值模拟方法中，可以根据构建的二度体模型得到对应的空间域重力异常表达式。

步骤106，对空间域重力异常表达式沿x轴方向进行一维傅里叶变换，得到波数域重力异常表达式；其中，波数域重力异常表达式包含波数域异常体剩余密度值；波数域异常体剩余密度值是对空间域异常体剩余密度值进行非均匀采样快速傅里叶变换得到的。

由于构建的二度体模型中剖分的矩形尺寸不一定一样，每个矩形对应的密度值不一定一样，因此，空间域的采样可以是非均匀的。通过对空间域重力异常表达式沿x轴方向进行一维傅里叶变换，得到波数域重力异常表达式，当剖分的矩形尺寸不一的时候，波数域异常体剩余密度值是对空间域异常体剩余密度值进行非均匀采样快速傅里叶变换得到的。

步骤108，根据矩形沿x轴方向的尺寸信息，确定截止频率，根据截止频率确定波数域波数采样值范围，根据波数采样值范围和预先设置的采样点总数，在波数域进行非均匀采样，得到波数采样值。

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

步骤110，将波数采样值代入所述波数域重力异常表达式，得到波数域重力异常场值。

相对于空间域重力异常表达式，波数域重力异常表达式的优点是表达式简洁，计算效率较高。

步骤112，通过对波数域重力异常场值进行一维非均匀快速傅里叶反变换，得到目标区域内任一点的空间域重力异常场值。

通过一维非均匀快速傅里叶反变换，将波数域重力异常场值转换到空间域重力异常场值，可以得到目标区域内任一点的空间域重力异常场值。

上述二度体重力异常积分解数值模拟方法中，通过构建包含目标区域的二度体模型，确定对应的空间域重力异常表达式，其中空间域重力异常表达式包含空间域异常体剩余密度值，通过对空间域重力异常表达式沿x轴方向进行一维傅里叶变换，得到波数域重力异常表达式，其中波数域异常体剩余密度值是对空间域异常体剩余密度值进行非均匀采样快速傅里叶变换得到的，根据矩形沿x轴方向的尺寸信息，确定截止频率，根据截止频率确定波数域波数采样值范围，根据波数采样值范围和预先设置的采样点总数进行波数域非均匀采样，得到波数采样值，根据波数采样值，代入波数域重力异常表达式中，得到波数域重力异常场值，最后通过对波数域重力异常场值进行一维非均匀快速傅里叶反变换，得到目标区域内任一点的空间域重力异常场值。采用本方法进行二度体重力异常积分解数值模拟，波数域非均匀采样可以根据重力异常场的频谱，对应波数域采样的间隔，实现根据波数域异常场变化规律进行采样，在波数域异常场变化较快的地方多采样，在场变化缓慢的地方少采样，以最少的采样个数达到数值模拟精度要求，进一步提高计算效率，有效兼顾了重力异常数值模拟的计算精度与计算效率，解决了波数域重力异常数值模拟只适合规则测线、数值模拟方法不能同时兼顾计算效率和计算精度的问题。

在其中一个实施例中，空间域重力异常表达式为：

式中：表示空间域重力异常场；表示万有引力常数；表示x方向剖分的矩形个数；表示z方向剖分的矩形个数；表示所述观测点坐标；表示编号为的矩形中心坐标；表示编号为的矩形的空间域异常体剩余密度值。

对本实施例中空间域重力异常表达式沿x轴方向进行一维傅里叶变换，得到波数域重力异常表达式为：

式中：表示波数域重力异常场；k表示波数；表示波数域异常体剩余密度值；表示符号函数：

其中，波数域异常体剩余密度值是对空间域异常体剩余密度值进行非均匀采样快速傅里叶变换得到的。

由于构建的二度体模型中，矩形尺寸是根据重力场的分布剖分的，不一定均匀剖分，本具体实施例中的波数域重力异常表达式可以实现在波数域进行非均匀采样，提高计算效率。

在其中一个实施例中，如图2所示，一种二度体重力异常积分解数值模拟方法主要包括四个部分：二度体模型表示、波数域重力异常公式推导、波数域非均匀采样和空间域重力异常计算。二度体模型构建包括设定目标区域及矩形几何尺寸，设定矩形个数和根据二度体形状对每个矩形进行赋值，以此刻画不同截面二度体；波数域重力异常公式推导是对空间域重力异常表达式进行一维傅里叶变换；波数域非均匀采样是根据重力异常频谱变换规律进行非均匀采样；空间域重力异常计算是对计算的非均匀采样频谱进行一维nufft反变换，得到空间域重力异常。

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

式中：i为虚数单位，为给定离散采样点对应采样点的值，为计算的离散采样点傅里叶变换频谱，n表示采样点总数；

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

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

式中：m表示过采样因子，表示权重因子，为精度因子，表示对取整；

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

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

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

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

在其中一个实施例中，将目标区域分成若干个矩形，每个矩形的x轴方向尺寸、z轴方向尺寸可以相同，也可以不同；在重力异常场变化快的地方加密网格，在重力异常场变化慢的地方稀疏网格。采用这种方式可以在空间域中根据重力异常场的分布情况决定采样间隔，在重力异常场变化快的地方多采样，在重力异常场变化慢的地方少采样，从而提高计算效率。

在其中一个实施例中，根据矩形沿x轴方向的尺寸信息，确定截止频率为：

式中，表示截止频率；表示x轴方向矩形尺寸的最小值。

根据截止频率确定波数域波数采样值范围为，再根据波数采样值范围得到波数采样值，包括：

在上波数依次取值为：

在上波数依次取值为：

式中，，n表示采样点总数。

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

在一个具体实施例中，目标区域有一个截面为圆形的二度体如图4所示，目标区域范围为：x轴方向从-128km到128km，z轴方向从0km到100km，和均为1km，可知x轴方向剖分矩形个数为256个，z轴方向剖分矩形个数为100个，圆心坐标为(0km,80km)，圆的半径为5km，剩余密度为10kg/m3，计算测线为水平地面上的重力异常。

利用fortran语言编程实现，运行程序所用的个人电脑配置为：cpu-intercorei7-8700，主频为3.2ghz，运行内存为8.00gb。

图5为空间域均匀剖分无限长圆柱体产生的重力异常波数域非均匀采样模拟结果与理论值对比图，其中，nufft表示非均匀快速傅里叶变换，n表示采样点总数。从图5中可以看出采用非均匀快速傅里叶变换选取150和256个波数域采样点总数，模拟结果与理论值吻合的很好，精度较高。图6表示本实施例中波数域非均匀采样模拟结果与理论值的相对误差，可以看出，波数域采样点总数为256时比波数域采样点总数为100时，模拟结果与理论值的相对误差更小，计算精度更高。

对于同样的二度体模型，统计了快速傅里叶变换扩边法（fft扩边法）、高斯快速傅里叶变换法（gauss-fft方法）与非均匀快速傅里叶变换法（nufft方法）计算的相对均方根误差及计算时间，对比如表1，可以看出在达到相同精度的情况下，非均匀快速傅里叶变换法计算所用时间更短，计算效率更高。

其中，4l、8l为快速傅里叶变换扩边法的扩边比，2、3分别为gauss-fft中采用高斯点的个数，rrms(%)表示相对均匀方根误差。

表示解析解，表示数值解。

在另一个具体实施例中，目标区域范围为：x轴方向从-128km到128km，z轴方向从0km到100km，在空间域进行非均匀采样，在(0km,15km)范围内，剖分的矩形为0.5km，在(15km,98km)范围内，矩形为1km，在(98km,128km)范围内，矩形为2km，-128km到0km关于（0,0）对称。采用非均匀快速傅里叶变换选取150和256个波数域采样点总数，计算结果如图7所示，可以看出，模拟结果与理论值吻合的很好，精度较高。图8表示本实施例中波数域非均匀采样模拟结果与理论值的相对误差，可以看出，波数域采样点总数为256时比波数域采样点总数为150时，模拟结果与理论值的相对误差更小，计算精度更高。

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

在一个实施例中，如图9所示，提供了一种二度体重力异常积分解数值模拟装置，包括：二度体模型构建模块902、空间域重力异常表达式确定模块904、波数域重力异常表达式确定模块906、波数采样值确定模块908、波数域重力异常场值确定模块910和空间域重力异常场值确定模块912，其中：

二度体模型构建模块902，用于根据观测点在观测点坐标中的位置信息，构建包含目标区域的二度体模型；观测点坐标包括x轴方向和z轴方向；二度体模型是将目标区域分成若干个矩形，矩形的截面形状用于确定每个矩形的空间域异常体剩余密度值；

空间域重力异常表达式确定模块904，用于根据二度体模型，确定对应的空间域重力异常表达式；空间域重力异常表达式包含空间域异常体剩余密度值；

波数域重力异常表达式确定模块906，用于对空间域重力异常表达式沿x轴方向进行一维傅里叶变换，得到波数域重力异常表达式；其中，波数域重力异常表达式包含波数域异常体剩余密度值；波数域异常体剩余密度值是对空间域异常体剩余密度值进行非均匀采样快速傅里叶变换得到的；

波数采样值确定模块908，用于根据矩形沿x轴方向的尺寸信息，确定截止频率，根据截止频率确定波数域波数采样值范围，根据波数采样值范围和预先设置的采样点总数，在波数域进行非均匀采样，得到波数采样值；

波数域重力异常场值确定模块910，用于将波数采样值代入波数域重力异常表达式，得到波数域重力异常场值；

空间域重力异常场值确定模块912，用于通过对波数域重力异常场值进行一维非均匀快速傅里叶反变换，得到目标区域内任一点的空间域重力异常场值。

二度体模型构建模块902还包括将目标区域分成若干个矩形，每个矩形的x轴方向尺寸、z轴方向尺寸可以相同，也可以不同；在重力异常场变化快的地方加密网格，在重力异常场变化慢的地方稀疏网格。

波数采样值确定模块908还包括根据所述矩形沿x轴方向的尺寸信息，确定截止频率为：

式中，表示截止频率；表示x轴方向矩形尺寸的最小值。

波数采样值确定模块908还包括根据截止频率确定波数域波数采样值范围为；

根据波数采样值范围得到波数采样值，包括：

在上波数依次取值为：

在上波数依次取值为：

式中，，n表示采样点总数。

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

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

本领域技术人员可以理解，图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）等。

以上实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。

以上所述实施例仅表达了本申请的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本申请构思的前提下，还可以做出若干变形和改进，这些都属于本申请的保护范围。因此，本申请专利的保护范围应以所附权利要求为准。