三维重力场全息数值模拟及CPU-GPU加速方法

摘要

本发明实现了一种三维重力场全息数值模拟及CPU‑GPU加速方法。利用水平方向二维傅里叶变换把空间域引力位满足的三维偏微分方程转化为不同波数满足的一维常微分方程，将三维问题转化为一维常微分问题求解，减少了计算量及存储需求，不同波数之间常微分方程相互独立，具有高度并行性；一维常微分方程上下边界问题具有严格的边界条件，与物理边界吻合；水平方向二维傅里叶变换采用全息傅里叶变换，与重力场水平方向边界吻合，因此三维偏微分方程描述的重力位物理信息得到完整、精确的数值模拟；基于算法的并行性，实现CPU‑GPU并行加速方案，兼顾计算效率与计算精度。该方法适用于任意复杂模型三维重力异常高效、高精度数值模拟。

说明书

技术领域

本发明涉及重力勘探领域，特别涉及一种三维重力场全息数值模拟及CPU-GPU加速方法。

背景技术

重力勘探广泛应用于油气勘探、金属矿勘查、地球内部构造、区域地质构造以及环境地球物理等研究领域，研究针对大规模复杂条件下重力异常场高效、高精度三维数值模拟对重力异常反演以及人机交互解释有着重要作用。

重力异常场数值模拟分为空间域和频率域，相比于空间域，频率域表达式简洁，计算效率高效，近年来，基于微分法（戴世坤，赵东东，张钱江等. 2018. 起伏地形条件下空间波数混合域重力异常三维数值模拟(英文). 应用地球物理，15(3-4): 513-523）和基于积分法（Dai S K, Chen Q R, Li K, et al. 2022. The forward modeling of 3D gravityand magnetic potential fields in space-wavenumber domains based on anintegral method. Geophysics, 2022, 87(3): G83-G96）的空间-波数域重力异常三维数值模拟等研究一定程度上提高了三维重力异常数值模拟的效率，但其效率和精度受限于傅里叶变换，仍有很大的提高空间。

发明内容

为了解决现有空间-波数域三维重力数值模拟中傅里叶变换存在的效率问题以及精度问题，本发明提供一种基于二次插值形函数的全息傅里叶变换方法，并通过GPU实现其并行加速。本发明中全息傅里叶变换方法的积分区间与重力场水平方向左右、前后边界吻合，减少截断效应。利用全息傅里叶算法以及空间-波数域算法的高度并行性，CPU多核并行执行五对角方程组求解，GPU多线程并行执行全息采样傅里叶变换，极大地提升了计算效率。

为了实现上述技术目的，本发明的技术方案是，

一种三维重力场全息数值模拟及CPU-GPU加速方法，包括以下步骤：

1）对重力位与空间域密度满足泊松方程进行水平方向二维傅里叶变换，得到二维波数域下的一维控制方程。

2）在直角坐标系中建立能够包围三维剩余密度异常体的最小长方体目标区域，并基于x，y和z方向对目标区域进行离散以实现剖分并形成离散节点。其中三维剩余密度异常体是根据计算需要建立的。

3）利用基于形函数插值的积分区间为无限的全息傅里叶变换的采样规则，根据x和y方向的剖分计算波数k

4）根据计算模型参数设置离散节点上的剩余密度。

5）计算基于形函数插值的二维全息傅里叶正变换需要的积分系数。

6）基于步骤5）的结果将目标区域的剩余密度进行水平方向全息傅里叶正变换，并基于GPU多线程加速，得到二维波数域剩余密度。

7）将二维波数域剩余密度代入二维波数域下的一维控制方程，采用基于二次插值的一维有限单元法，并加载上、下边界条件，得到五对角方程组。再采用CPU上的OpenMP加速以追赶法求解方程组，得到二维波数域下的重力位。最后利用空间-波数域重力位与重力场、梯度张量之间的关系，得到空间-波数域重力场以及梯度张量。

8）计算基于形函数插值的二维全息傅里叶反变换需要的积分系数。

9）将空间-波数域重力场进行二维全息傅里叶反变换，并基于GPU多线程加速，得到空间域重力场以及梯度张量。

所述的方法，所述的步骤1）中，重力位与空间域密度满足泊松方程为：

其中，

然后沿水平方向做二维傅里叶正变换，得到二维波数域下的一维控制方程：

其中，

且在对基于一维控制方程表示的一系列平面波进行求解时，应结合上下边界条件：

其中

所述的方法，所述的步骤2）中，基于x，y和z方向对目标区域进行离散以实现剖分时，是通过均匀剖分或不均匀剖分来将目标区域剖分为多个空间单元，从而形成多个离散节点，离散节点的数量即空间单元的数量+1。

所述的方法，所述的步骤3）中，采样规则为：在频谱能量大且变化剧烈的区间加密采样，在频谱能量小且变化缓慢的区间稀疏采样，其中区间为波数（

其中t

所述的方法，所述的步骤4）中，给定每个节点的剩余密度为ρ，其中异常体剩余密度与周围剩余密度不同。

所述的方法，所述的步骤5）中，通过下式计算积分系数：

其中

所述的方法，所述的步骤6）中，是根据步骤5得到的积分系数，对目标区域的剩余密度

首先根据下式对剩余密度

Fx(kx,y)=∑j=1Mx[f(xj1,y)⋅I(xj1,kx)+f(xj2,y)⋅I(xj2,kx)+f(xj3,y)⋅I(xj3,kx)]

其中

计算过程中利用GPU进行多线程加速，以GPU的各个线程单元来对每个波数k

然后根据下式对

F(kx,ky)=∑j=1My[Fx(kx,yj1)⋅J(yj1,ky)+Fx(kx,yj2)⋅J(yj2,ky)+Fx(kx,yj3)⋅J(yj3,ky)]

其中

计算过程中同样利用GPU进行多线程加速，以GPU的各个线程单元来对每个波数k

所述的方法，所述的步骤7）中，基于二次插值函数的一维有限单元法构成对角方程组Ku=P，其中K为五对角阵的系数矩阵。u为求解空间-波数域重力位未知量，P为右端项。从而形成N

g˜=[g˜xg˜yg˜z]=[ikxU˜ikyU˜∂U˜∂z]

T̃=[T˜xxT˜yxT˜zxT˜xyT˜yyT˜zyT˜xzT˜yzT˜zz]=[−kx2U˜−kxkyU˜−kx∂U˜∂z−kxkyU˜−ky2U˜−ky∂U˜∂z−kx∂U˜∂z−ky∂U˜∂z∂2U˜∂z2]

其中

所述的方法，所述的步骤8）中，通过下式计算积分系数，

其中

所述的方法，所述步骤9）中，利用步骤8中计算出来的积分系数，对波数域的重力异常场以及梯度张量沿着水平方向进行二维全息傅里叶反变换：

Fx(x,ky)=∑j=1Mkx[F(kxj1,ky)⋅I−(kxj1,x)+F(kxj2,ky)⋅I−(kxj2,x)+F(kxj2,ky)⋅I−(kxj2,x)]

计算过程中利用GPU进行多线程加速，以GPU的各个线程单元来对每个空间节点x下的沿k

然后对

f(x,y)=∑j=1Mky[Fx(x,kyj1)⋅J−(kyj1,y)+Fx(x,kyj2)⋅J−(kyj2,y)+Fx(x,kyj3)⋅J−(kyj3,y)]

计算过程中利用GPU进行多线程加速，以GPU的各个线程单元来对每个空间节点y下的沿k

本发明的技术效果在于，本发明实现了基于二次插值形函数任意采样傅里叶变换的空间-波数域重力异常数值模拟算法，并采用CPU-GPU进行加速，实现了三维重力异常高效、高精度数值模拟。利用水平方向二维傅里叶变换把空间域引力位满足的三维偏微分方程转化为不同波数满足的一维常微分方程，将三维问题转化为一维常微分问题求解，避免了大型系数矩阵方程组的求解，减少了计算量及存储需求，而且节点规模越大，效果越明显，不同波数之间常微分方程相互独立，具有高度并行性。常微分方程上下边界问题具有严格的边界条件，与物理边界吻合，解决上下边界问题。水平方向二维傅里叶变换采用任意采样傅里叶变换，与重力场水平方向左右、前后边界吻合，解决重力场水平方向左右、前后边界问题，数值模拟结果与重力场解析解高度吻合。基于空间波数域算法以及任意傅里叶方法的并行性，在CPU上并行求解五对角方程组，GPU上计算正反二维傅里叶变换，实现CPU-GPU并行加速方案，兼顾计算效率与计算精度。该方法适用于任意复杂模型三维重力异常高效、高精度数值模拟，为大规模重力实测数据的高效、精细反演提供重要技术支持。

附图说明

图1为本发明实施例中棱柱体模型的平面图。

图2为本发明实施例中重力异常场在地面z = 0 m数值解的平面等值线图。

图3为本发明实施例中重力异常场在地面z = 0 m解析解的平面等值线图。

具体实施方式

本实施例所提供的基于任意采样傅里叶变换的空间-波数域重力异常数值模拟及CPU-GPU加速方法，包括以下步骤：

1）引入重力位满足的泊松方程，对该方程进行水平方向二维傅里叶变换，得到二维波数域下的一维控制方程组。

这里的重力位

其中，

其中，

这里在对基于一维控制方程表示的一系列平面波进行求解时，应结合上下边界条件：

其中

2）在直角坐标系中建立一个长方体目标区域，其中三维剩余密度异常体在目标区内，将目标区域x，y和z方向进行离散。其中长方体目标区域是能够将整个三维剩余密度异常体全部包裹在内的最小的长方体形状的区域。而三维剩余密度异常体的具体参数是根据计算的需要来建立的。其中对长方体目标区域的x、y、z方向可以均匀剖分，也可以不均匀剖分，从而将该区域分为若干空间单元，x、y、z方向的剖分节点数量分别设为N

3）利用基于形函数插值的任意采样傅里叶变换的采样规则，根据x和y方向的剖分计算波数k

本实施例中利用任意采样傅里叶变换的采样规则：波数选取范围延用采样定理，在频谱能量大且变化剧烈的区间加密采样，在频谱能量小且变化缓慢的区间稀疏采样，区间是指的波数（k

其中t

4）根据计算模型参数设置离散节点上的剩余密度。本实施例中给定每个节点的剩余密度为ρ，其中异常体剩余密度与周围剩余密度不同。

5）计算基于形函数插值的二维任意傅里叶正变换需要的积分系数。

对于基于形函数的任意傅里叶正变换，其核心思想是将二维傅里叶正变换转化为两个一维傅里叶正变换，并把一维傅里叶正变换离散为多个单元积分累加和，而离散单元中的原函数可以用二次插值形函数拟合，然后求出单元积分的解析表达式，最终表示为：

Fx(kx,y)=∑j=1Mx[f(xj1,y)⋅I(xj1,kx)+f(xj2,y)⋅I(xj2,kx)+f(xj3,y)⋅I(xj3,kx)]

F(kx,ky)=∑j=1My[Fx(kx,yj1)⋅J(yj1,ky)+Fx(kx,yj2)⋅J(yj2,ky)+Fx(kx,yj3)⋅J(yj3,ky)]

其中

其中e表示自然对数的底数，i为虚数单位。根据式(6)和(7)，x、y、k

6）将目标区域的剩余密度进行水平方向任意傅里叶正变换，采用GPU多线程加速，得到二维波数域剩余密度。

本步骤中是利用上一步骤中计算出来的积分系数，对目标区域的剩余密度

Fx(kx,y)=∑j=1Mx[f(xj1,y)⋅I(xj1,kx)+f(xj2,y)⋅I(xj2,kx)+f(xj3,y)⋅I(xj3,kx)]

其中

由于每个波数

第二步根据式(5)对

F(kx,ky)=∑j=1My[Fx(kx,yj1)⋅J(yj1,ky)+Fx(kx,yj2)⋅J(yj2,ky)+Fx(kx,yj3)⋅J(yj3,ky)]

其中

另外对于z方向来说，其每一层都要做一次x和y方向上的二维任意采样傅里叶变换。

7）将二维波数域剩余密度带入二维波数域一维空间域控制方程，采用基于二次插值的一维有限单元法，并加载上、下边界条件，得到五对角方程组，OpenMP加速实现追赶法求解方程组，得到一维空间域重力位，再利用空间-波数域重力位与重力场、梯度张量之间的关系，得到空间-波数域重力场以及梯度张量。

本步骤中基于二次插值函数的一维有限单元法构成对角方程组Ku=P，其中K为系数矩阵，该矩阵为五对角阵。u为求解空间-波数域重力位未知量，P为右端项。这种五对角方程组的个数有N

g˜=[g˜xg˜yg˜z]=[ikxU˜ikyU˜∂U˜∂z]

T̃=[T˜xxT˜yxT˜zxT˜xyT˜yyT˜zyT˜xzT˜yzT˜zz]=[−kx2U˜−kxkyU˜−kx∂U˜∂z−kxkyU˜−ky2U˜−ky∂U˜∂z−kx∂U˜∂z−ky∂U˜∂z∂2U˜∂z2]

可求得空间-波数域的重力场以及重力梯度张量。其中

8）计算基于形函数插值的二维任意傅里叶反变换需要的积分系数。

对于基于形函数的任意傅里叶反变换计算过程类似二维正变换，正变换分别沿x,y方向进行一维积分，而反变换分别对

通过下式计算积分系数，

其中

9）将空间-波数域重力场进行二维任意傅里叶反变换，并采用GPU加速，得到空间域重力场以及梯度张量。

利用步骤8中计算出来的积分系数，对波数域的重力异常场以及梯度张量沿着水平方向进行二维任意傅里叶反变换，计算加速过程类似正变换过程。

Fx(x,ky)=∑j=1Mkx[F(kxj1,ky)⋅I−(kxj1,x)+F(kxj2,ky)⋅I−(kxj2,x)+F(kxj2,ky)⋅I−(kxj2,x)]

计算过程中利用GPU进行多线程加速，以GPU的各个线程单元来对每个空间节点x下的沿k

然后对

f(x,y)=∑j=1Mky[Fx(x,kyj1)⋅J−(kyj1,y)+Fx(x,kyj2)⋅J−(kyj2,y)+Fx(x,kyj3)⋅J−(kyj3,y)]

计算过程中利用GPU进行多线程加速，以GPU的各个线程单元来对每个空间节点y下的沿k

基于以上步骤，进行数值模拟验证过程如下：

设定如图1的棱柱体模型，计算范围x方向-10 km~10 km，y方向-10 km~10 km，z方向计算范围设为0~7.5 km，剖分网格节点个数201×201×201，三个方向均匀剖分，水平方向网格间距均为100 m，垂直方向50 m，异常体平面位置如图1所示，4个角上得异常体位置对称，距离边界200 m，异常体处于同一深度，顶部埋深为500 m，异常体大小均为2 km×2km×2 km，异常体剩余密度

图2、图3分别为重力异常场在地面z = 0 m数值解以及解析解的平面等值线图。由图可知，数值解和解析解等值线图形态相似，场值以及梯度张量吻合得很好，且在z = 0 m平面

本实施例是基于以下配置计算机进行计算：CPU Intel(R) Core(TM) i9-9980XE，主频3.0GHZ，128 GBRAM. GPU显卡型号NVIDIA TITAN RTX，CUDA版本号11.6，运行于Linux编译环境。本实施例的计算节点数为201×201×101=15813251，计算总时间为0.7s，占用CPU+GPU内存为819+363MB，由此可见本发明计算速度快，占用内存小，具有很高的计算效率。