CN116911146B - 三维重力场全息数值模拟及cpu-gpu加速方法 - Google Patents

三维重力场全息数值模拟及cpu-gpu加速方法 Download PDF

Info

Publication number
CN116911146B
CN116911146B CN202311184544.3A CN202311184544A CN116911146B CN 116911146 B CN116911146 B CN 116911146B CN 202311184544 A CN202311184544 A CN 202311184544A CN 116911146 B CN116911146 B CN 116911146B
Authority
CN
China
Prior art keywords
dimensional
domain
wave number
thread
space
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
Application number
CN202311184544.3A
Other languages
English (en)
Other versions
CN116911146A (zh
Inventor
戴世坤
朱德祥
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Central South University
Original Assignee
Central South University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Central South University filed Critical Central South University
Priority to CN202311184544.3A priority Critical patent/CN116911146B/zh
Publication of CN116911146A publication Critical patent/CN116911146A/zh
Application granted granted Critical
Publication of CN116911146B publication Critical patent/CN116911146B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • 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
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F9/00—Arrangements for program control, e.g. control units
    • G06F9/06—Arrangements for program control, e.g. control units using stored programs, i.e. using an internal store of processing equipment to receive or retain programs
    • G06F9/46—Multiprogramming arrangements
    • G06F9/50—Allocation of resources, e.g. of the central processing unit [CPU]
    • G06F9/5005—Allocation of resources, e.g. of the central processing unit [CPU] to service a request
    • G06F9/5027—Allocation of resources, e.g. of the central processing unit [CPU] to service a request the resource being a machine, e.g. CPUs, Servers, Terminals
    • 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/14—Force analysis or force optimisation, e.g. static or dynamic forces
    • 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
    • Y02D—CLIMATE CHANGE MITIGATION TECHNOLOGIES IN INFORMATION AND COMMUNICATION TECHNOLOGIES [ICT], I.E. INFORMATION AND COMMUNICATION TECHNOLOGIES AIMING AT THE REDUCTION OF THEIR OWN ENERGY USE
    • Y02D10/00—Energy efficient computing, e.g. low power processors, power management or thermal management

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • Algebra (AREA)
  • Computational Mathematics (AREA)
  • Geometry (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Evolutionary Computation (AREA)
  • Pure & Applied Mathematics (AREA)
  • Databases & Information Systems (AREA)
  • Computer Hardware Design (AREA)
  • Complex Calculations (AREA)

Abstract

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

Description

三维重力场全息数值模拟及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 gravity andmagnetic potential fields in space-wavenumber domains based on an integralmethod.Geophysics,2022,87(3):G83-G96)的空间-波数域重力异常三维数值模拟等研究一定程度上提高了三维重力异常数值模拟的效率，但其效率和精度受限于傅里叶变换，仍有很大的提高空间。
发明内容
为了解决现有空间-波数域三维重力数值模拟中傅里叶变换存在的效率问题以及精度问题，本发明提供一种基于二次插值形函数的全息傅里叶变换方法，并通过GPU实现其并行加速。本发明中全息傅里叶变换方法的积分区间与重力场水平方向左右、前后边界吻合，减少截断效应。利用全息傅里叶算法以及空间-波数域算法的高度并行性，CPU多核并行执行五对角方程组求解，GPU多线程并行执行全息采样傅里叶变换，极大地提升了计算效率。
为了实现上述技术目的，本发明的技术方案是，
一种三维重力场全息数值模拟及CPU-GPU加速方法，包括以下步骤：
1)对重力位与空间域密度满足泊松方程进行水平方向二维傅里叶变换，得到二维波数域下的一维控制方程。
2)在直角坐标系中建立能够包围三维剩余密度异常体的最小长方体目标区域，并基于x，y和z方向对目标区域进行离散以实现剖分并形成离散节点。其中三维剩余密度异常体是根据计算需要建立的。
3)利用基于形函数插值的积分区间为无限的全息傅里叶变换的采样规则，根据x和y方向的剖分计算波数kx和ky。
4)根据计算模型参数设置离散节点上的剩余密度。
5)计算基于形函数插值的二维全息傅里叶正变换需要的积分系数。
6)基于步骤5的结果将目标区域的剩余密度进行水平方向全息傅里叶正变换，并基于GPU多线程加速，得到二维波数域剩余密度。
7)将二维波数域剩余密度代入二维波数域下的一维控制方程，采用基于二次插值的一维有限单元法，并加载上、下边界条件，得到五对角方程组。再采用CPU上的OpenMP加速以追赶法求解方程组，得到二维波数域下的重力位。最后利用空间-波数域重力位与重力场、梯度张量之间的关系，得到空间-波数域重力场以及梯度张量。
8)计算基于形函数插值的二维全息傅里叶反变换需要的积分系数。
9)将空间-波数域重力场进行二维全息傅里叶反变换，并基于GPU多线程加速，得到空间域重力场以及梯度张量。
所述的方法，所述的步骤1中，重力位与空间域密度满足泊松方程为：
其中，为拉普拉斯算子，U为空间域重力位，(x,y,z)为空间观测点坐标，π为圆周率，G为万有引力常量，ρ为空间域密度。
然后沿水平方向做二维傅里叶正变换，得到二维波数域下的一维控制方程：
其中，表示偏导，/>是空间-波数域重力位，/>是空间-波数混合域密度，kx和ky分别是x和y方向傅里叶变换后的波数，k为中间变量，/>
且在对基于一维控制方程表示的一系列平面波进行求解时，应结合上下边界条件：
其中zmin代表上边界，zmax代表下边界。
所述的方法，所述的步骤2中，基于x，y和z方向对目标区域进行离散以实现剖分时，是通过均匀剖分或不均匀剖分来将目标区域剖分为多个空间单元，从而形成多个离散节点，离散节点的数量即空间单元的数量+1。
所述的方法，所述的步骤3中，采样规则为：在频谱能量大且变化剧烈的区间加密采样，在频谱能量小且变化缓慢的区间稀疏采样，其中区间为波数(kx、ky)的范围区间。其中波数的基本单位即基频的表达式为：
其中tmax为目标区域x或y方向计算范围的最大值，tmin为空间域x或y方向计算范围的最小值。当空间域x、y方向上的采样点数为Nx和Ny、采样间隔为Δx和Δy时，根据采样规则求得x和y方向对应的波数kx和ky，波数节点数量分别为Nkx、Nky，π为圆周率。
所述的方法，所述的步骤4中，给定每个节点的剩余密度为ρ，其中异常体剩余密度与周围剩余密度不同。
所述的方法，所述的步骤5中，通过下式计算积分系数：
其中I(xj1,kx)，I(xj2,kx)，I(xj3,kx)分别为沿x方向第j个空间单元上三个节点的傅里叶正变换积分系数，J(yj1,ky)，J(yj2,ky)，J(yj3,ky)分别为沿y方向上第j个空间单元的傅里叶正变换节点积分系数，e表示自然对数的底数，i为虚数单位。
所述的方法，所述的步骤6中，是根据步骤5得到的积分系数，对目标区域的剩余密度ρ(x,y,z)沿水平方向进行二维全息傅里叶正变换：
首先根据下式对剩余密度ρ(x,y,z)沿x方向进行一维积分计算得到ρx(kx,y,z)：
其中Fx(kx,y)表示在不同z值下的ρx(kx,y,z)；
计算过程中利用GPU进行多线程加速，以GPU的各个线程单元来对每个波数kx中的沿x方向傅里叶变换积分进行并行计算；其中线程单元采用二维线程块和由二维线程块所组成的二维线程网格来组织，线程单元之间通过线程索引(threadIdx％x,threadIdx％y)和线程网格索引(blockIdx％x,blockIdx％y)定位自己在线程网格中的索引，以使每个线程单元对应到每个线程数据的逻辑位置，其中线程数据包括要进行傅里叶变换的数据f(xj,y)以及积分系数I(xj,kx)；一个线程块中有a×b个线程单元，线程网格中有m×n个线程块，线程块中的一个线程单元负责一个波数计算：即Nx次乘法和加法运算，因此有：m＝(a+Nkx-1)/a，n＝(b+Ny-1)/b，Nkx表示kx方向上的波数个数；
然后根据下式对ρx(kx,y,z)沿y方向进行一维积分计算得到波数域下的
其中F(kx,ky)表示在不同z值下的
计算过程中同样利用GPU进行多线程加速，以GPU的各个线程单元来对每个波数ky中的沿y方向傅里叶变换积分进行并行计算。
所述的方法，所述的步骤7中，基于二次插值函数的一维有限单元法构成对角方程组Ku＝P，其中K为五对角阵的系数矩阵。u为求解空间-波数域重力位未知量，P为右端项。从而形成Nkx×Nky个五对角方程组。然后利用二维波数域下重力位与重力场、梯度张量之间的关系，采用OpenMP加速以追赶法求解方程组得到二维波数域下的重力位：
其中指的是二维波数域下的重力位，/>为二维波数域下的重力场，/> 和/>分别为x、y和z三个方向上的分量，/>为二维波数域下的重力梯度张量，为9个梯度张量分量，i为虚数单位。
所述的方法，所述的步骤8中，通过下式计算积分系数，
其中分别为沿kx方向第j个空间单元上三个节点的傅里叶反变换积分系数，/> 分别为沿ky方向上第j个空间单元的傅里叶反变换节点积分系数，e表示自然对数的底数，i为虚数单位。
所述的方法，所述步骤9中，利用步骤8中计算出来的积分系数，对波数域的重力异常场以及梯度张量沿着水平方向进行二维全息傅里叶反变换：
F(kx，ky)表示波数域的重力异常场、梯度张量，沿kx方向进行一维积分得到Fx(x，kky)：
计算过程中利用GPU进行多线程加速，以GPU的各个线程单元来对每个空间节点x下的沿kx方向傅里叶变换积分进行并行计算。
然后对Fx(x，ky)沿ky方向进行一维积分得到f(x，y)，f(x，y)表示空间域重力异常场以及梯度张量：
计算过程中利用GPU进行多线程加速，以GPU的各个线程单元来对每个空间节点y下的沿ky方向傅里叶变换积分进行并行计算。
本发明的技术效果在于，本发明实现了基于二次插值形函数任意采样傅里叶变换的空间-波数域重力异常数值模拟算法，并采用CPU-GPU进行加速，实现了三维重力异常高效、高精度数值模拟。利用水平方向二维傅里叶变换把空间域引力位满足的三维偏微分方程转化为不同波数满足的一维常微分方程，将三维问题转化为一维常微分问题求解，避免了大型系数矩阵方程组的求解，减少了计算量及存储需求，而且节点规模越大，效果越明显，不同波数之间常微分方程相互独立，具有高度并行性。常微分方程上下边界问题具有严格的边界条件，与物理边界吻合，解决上下边界问题。水平方向二维傅里叶变换采用任意采样傅里叶变换，与重力场水平方向左右、前后边界吻合，解决重力场水平方向左右、前后边界问题，数值模拟结果与重力场解析解高度吻合。基于空间波数域算法以及任意傅里叶方法的并行性，在CPU上并行求解五对角方程组，GPU上计算正反二维傅里叶变换，实现CPU-GPU并行加速方案，兼顾计算效率与计算精度。该方法适用于任意复杂模型三维重力异常高效、高精度数值模拟，为大规模重力实测数据的高效、精细反演提供重要技术支持。
附图说明
图1为本发明实施例中棱柱体模型的平面图。
图2为本发明实施例中重力异常场在地面z＝0m数值解的平面等值线图。
图3为本发明实施例中重力异常场在地面z＝0m解析解的平面等值线图。
具体实施方式
本实施例所提供的基于任意采样傅里叶变换的空间-波数域重力异常数值模拟及CPU-GPU加速方法，包括以下步骤：
1)引入重力位满足的泊松方程，对该方程进行水平方向二维傅里叶变换，得到二维波数域下的一维控制方程组。
这里的重力位U与空间域密度满足泊松方程(Blakely R J.1996.Potentialtheory in gravity and magnetic applications:London,Cambridge UniversityPress.)：
其中，为拉普拉斯算子，U为空间域重力位，(x,y,z)为空间观测点坐标，π为圆周率，G为万有引力常量，值为6.67×10-11N·m3/kg2，ρ为空间域密度。再对式(1)沿水平方向做二维傅里叶正变换，得到式(2)空间-波数域常微分方程即二维波数域下的一维控制方程：
其中，表示偏导，/>是空间-波数域重力位，/>是空间-波数混合域密度，kx和ky分别是x和y方向傅里叶变换后的波数，k为中间变量，/>
这里在对基于一维控制方程表示的一系列平面波进行求解时，应结合上下边界条件：
其中zmin代表上边界，zmax代表下边界。从而使边界条件与重力位物理边界条件一致，空间-波数域重力异常上下边界问题得到准确表示。
2)在直角坐标系中建立一个长方体目标区域，其中三维剩余密度异常体在目标区内，将目标区域x，y和z方向进行离散。其中长方体目标区域是能够将整个三维剩余密度异常体全部包裹在内的最小的长方体形状的区域。而三维剩余密度异常体的具体参数是根据计算的需要来建立的。其中对长方体目标区域的x、y、z方向可以均匀剖分，也可以不均匀剖分，从而将该区域分为若干空间单元，x、y、z方向的剖分节点数量分别设为Nx、Ny、Nz。这样也就形成了多个离散节点，离散节点的数量即空间单元的数量+1。
3)利用基于形函数插值的任意采样傅里叶变换的采样规则，根据x和y方向的剖分计算波数kx和ky。其中傅里叶变换是本实施例中所采用的积分区间为无限的全息傅里叶变换，其推导过程在后续的步骤5中进行具体说明。
本实施例中利用任意采样傅里叶变换的采样规则：波数选取范围延用采样定理，在频谱能量大且变化剧烈的区间加密采样，在频谱能量小且变化缓慢的区间稀疏采样，区间是指的波数(kx、kxy)的范围区间。基频表达式为(陈龙伟,戴世坤,吴美平.应用任意采样点数FFT算法时离散频率计算.地球物理学进展,2016,31(01):164-169)：
其中tmax为空间域x或y方向计算范围的最大值，tmin为空间域x或y方向计算范围的最小值。当空间域水平方向x、y上采样点数Nx和Ny、采样间隔为Δx和Δy时，根据上述采样规则求得x和y方向对应的波数kx和ky，波数节点数量分别为Nkx、Nky。
4)根据计算模型参数设置离散节点上的剩余密度。本实施例中给定每个节点的剩余密度为ρ，其中异常体剩余密度与周围剩余密度不同。
5)计算基于形函数插值的二维任意傅里叶正变换需要的积分系数。
对于基于形函数的任意傅里叶正变换，其核心思想是将二维傅里叶正变换转化为两个一维傅里叶正变换，并把一维傅里叶正变换离散为多个单元积分累加和，而离散单元中的原函数可以用二次插值形函数拟合，然后求出单元积分的解析表达式，最终表示为：
其中I(xj1,kx)，I(xj2,kx)，I(xj3,kx)为沿x方向第j个单元上三个节点的傅里叶变换积分系数，J(yj1,ky)，J(yj2,ky)，J(yj3,ky)分别为沿y方向上第j个单元的傅里叶变换节点积分系数，该积分系数具体为：
其中e表示自然对数的底数，i为虚数单位。根据式(6)和(7)，x、y、kx以及ky确定以后，即可求出任意采样傅里叶变换沿两个方向的积分系数，并且当x、y方向剖分固定时，单元形函数不变，积分系数也不会改变，因此可以将积分系数提前计算出来并存储，以重复利用，从而减少冗余计算，提高计算效率。
6)将目标区域的剩余密度进行水平方向任意傅里叶正变换，采用GPU多线程加速，得到二维波数域剩余密度。
本步骤中是利用上一步骤中计算出来的积分系数，对目标区域的剩余密度ρ(x,y,z)沿水平方向进行二维任意傅里叶正变换，第一步根据式(4)对剩余密度ρ(x,y,z)沿x方向进行一维积分计算得到ρx(kx,y,z)，具体的计算公式为：
其中Fx(kx,y)表示在不同z值下的ρx(kx,y,z)。
由于每个波数kx下的积分计算相互独立，因此沿x方向傅里叶变换积分计算可以并行执行。本实施例中利用GPU众多的线程单元来进行并行计算，这里线程单元采用二维线程块和二维线程网格来组织，GPU的线程单元间通过线程索引(threadIdx％x,threadIdx％y)和线程网格索引(blockIdx％x,blockIdx％y)定位自己在线程网格中的索引，这样来使每个GPU的线程单元准确找到线程数据对应的逻辑位置。一个线程块中有a×b个线程，线程网格中有m×n个线程块，线程块中的一个线程单元负责一个波数计算：即Nx次乘法和加法运算。因此有：m＝(a+Nkx-1)/a，n＝(b+Ny-1)/b，Nkx表示kx方向上下的波数个数。
第二步根据式(5)对ρx(kx,y,z)沿着y方向进行一维积分计算得到波数域下的计算过程类似沿x方向，同样采用GPU多线程并行加速，即根据下式对ρx(kx,y,z)沿y方向进行一维积分计算得到波数域下的/>
其中F(kx,ky)表示在不同z值下的
另外对于z方向来说，其每一层都要做一次x和y方向上的二维任意采样傅里叶变换。
7)将二维波数域剩余密度带入二维波数域一维空间域控制方程，采用基于二次插值的一维有限单元法，并加载上、下边界条件，得到五对角方程组，OpenMP加速实现追赶法求解方程组，得到一维空间域重力位，再利用空间-波数域重力位与重力场、梯度张量之间的关系，得到空间-波数域重力场以及梯度张量。
本步骤中基于二次插值函数的一维有限单元法构成对角方程组Ku＝P，其中K为系数矩阵，该矩阵为五对角阵。u为求解空间-波数域重力位未知量，P为右端项。这种五对角方程组的个数有Nkx×Nky个。本实施例通过将一个三维的大规模方程组，转换为有Nkx×Nky个独立的一维有限元方程组，以达到计算量小、存储需求少以及并行性高的效果。利用二维波数域下重力位与重力场、梯度张量之间的关系，采用OpenMP加速以追赶法求解方程组得到二维波数域下的重力位。
可求得空间-波数域的重力场以及重力梯度张量。其中指的是二维波数域下的重力位，/>为二维波数域下的重力场，/>和/>分别为x、y和z三个方向上的分量，/>为二维波数域下的重力梯度张量，/> 为9个梯度张量分量，i为虚数单位。
8)计算基于形函数插值的二维任意傅里叶反变换需要的积分系数。
对于基于形函数的任意傅里叶反变换计算过程类似二维正变换，正变换分别沿x,y方向进行一维积分，而反变换分别对kx，ky进行一维积分运算，同样也需要计算反变换需要的积分系数，计算过程类似正变换积分系数过程。
通过下式计算积分系数，
其中分别为沿kx方向第j个空间单元上三个节点的傅里叶反变换积分系数，/> 分别为沿ky方向上第j个空间单元的傅里叶反变换节点积分系数，e表示自然对数的底数，i为虚数单位。
9)将空间-波数域重力场进行二维任意傅里叶反变换，并采用GPU加速，得到空间域重力场以及梯度张量。
利用步骤8中计算出来的积分系数，对波数域的重力异常场以及梯度张量沿着水平方向进行二维任意傅里叶反变换，计算加速过程类似正变换过程。
F(kx，ky)表示波数域的重力异常场、梯度张量，沿kx方向进行一维积分得到Fx(x，ky)：
计算过程中利用GPU进行多线程加速，以GPU的各个线程单元来对每个空间节点x下的沿kx方向傅里叶变换积分进行并行计算；
然后对Fx(x，ky)沿ky方向进行一维积分得到f(x，y)，f(x，y)表示空间域重力异常场以及梯度张量：
计算过程中利用GPU进行多线程加速，以GPU的各个线程单元来对每个空间节点y下的沿ky方向傅里叶变换积分进行并行计算。
基于以上步骤，进行数值模拟验证过程如下：
设定如图1的棱柱体模型，计算范围x方向-10km～10km，y方向-10km～10km，z方向计算范围设为0～7.5km，剖分网格节点个数201×201×201，三个方向均匀剖分，水平方向网格间距均为100m，垂直方向50m，异常体平面位置如图1所示，4个角上得异常体位置对称，距离边界200m，异常体处于同一深度，顶部埋深为500m，异常体大小均为2km×2km×2km，异常体剩余密度ρ1＝2000kg·m-3，ρ1＝2000kg·m-3。根据步骤3确定波数选取范围为-3.9×10-2～3.9×10-2，波数选取个数为253×253。用解析解的计算结果为参照，验证方法的正确性。
图2、图3分别为重力异常场在地面z＝0m数值解以及解析解的平面等值线图。由图可知，数值解和解析解等值线图形态相似，场值以及梯度张量吻合得很好，且在z＝0m平面gx、gy、gz、Txx、Tyy、Tzz、Txy、Txz以及Tyz9个分量的相对均方根误差分别为0.057％，0.057％，0.054％，0.20％，0.20％，0.19％，0.12％，0.20％、0.20％，异常体内部的重力异常场的相对均方根误差均小于1％，表明基于任意采样傅里叶变换的空间-波数域三维重力异常算法正确。
本实施例是基于以下配置计算机进行计算：CPU Intel(R)Core(TM)i9-9980XE，主频3.0GHZ，128GBRAM.GPU显卡型号NVIDIA TITAN RTX，CUDA版本号11.6，运行于Linux编译环境。本实施例的计算节点数为201×201×101＝15813251，计算总时间为0.7s，占用CPU+GPU内存为819+363MB，由此可见本发明计算速度快，占用内存小，具有很高的计算效率。

Claims (7)

1.一种三维重力场全息数值模拟及CPU-GPU加速方法，其特征在于，包括以下步骤：
1)对重力位与空间域密度满足泊松方程进行水平方向二维傅里叶变换，得到二维波数域下的一维控制方程；
2)在直角坐标系中建立能够包围三维剩余密度异常体的最小长方体目标区域，并基于x，y和z方向对目标区域进行离散以实现剖分并形成离散节点；其中三维剩余密度异常体是根据计算需要建立的；
3)利用基于形函数插值的积分区间为无限的全息傅里叶变换的采样规则，根据x和y方向的剖分计算波数kx和ky；
4)根据计算模型参数设置离散节点上的剩余密度；
5)计算基于形函数插值的二维全息傅里叶正变换需要的积分系数；
6)基于步骤5的结果将目标区域的剩余密度进行水平方向全息傅里叶正变换，并基于GPU多线程加速，得到二维波数域剩余密度；
7)将二维波数域剩余密度代入二维波数域下的一维控制方程，采用基于二次插值的一维有限单元法，并加载上、下边界条件，得到五对角方程组；再采用CPU上的OpenMP加速以追赶法求解方程组，得到二维波数域下的重力位；最后利用空间-波数域重力位与重力场、梯度张量之间的关系，得到空间-波数域重力场以及梯度张量；
8)计算基于形函数插值的二维全息傅里叶反变换需要的积分系数；
9)将空间-波数域重力场进行二维全息傅里叶反变换，并基于GPU多线程加速，得到空间域重力场以及梯度张量；
所述的步骤6中，是根据步骤5得到的积分系数，对目标区域的剩余密度ρ(x,y,z)沿水平方向进行二维全息傅里叶正变换：
首先根据下式对剩余密度ρ(x,y,z)沿x方向进行一维积分计算得到ρx(kx,y,z)：
其中Fx(kx,y)表示在不同z值下的ρx(kx,y,z)；
计算过程中利用GPU进行多线程加速，以GPU的各个线程单元来对每个波数kx中的沿x方向傅里叶变换积分进行并行计算；其中线程单元采用二维线程块和由二维线程块所组成的二维线程网格来组织，线程单元之间通过线程索引(threadIdx％x,threadIdx％y)和线程网格索引(blockIdx％x,blockIdx％y)定位自己在线程网格中的索引，以使每个线程单元对应到每个线程数据的逻辑位置，其中线程数据包括要进行傅里叶变换的数据f(xj,y)以及积分系数I(xj,kx)；一个线程块中有a×b个线程单元，线程网格中有m×n个线程块，线程块中的一个线程单元负责一个波数计算：即Nx次乘法和加法运算，因此有：m＝(a+Nkx-1)/a，n＝(b+Ny-1)/b，Nkx表示kx方向上的波数个数；
然后根据下式对ρx(kx,y,z)沿y方向进行一维积分计算得到波数域下的
其中F(kx,ky)表示在不同z值下的
计算过程中同样利用GPU进行多线程加速，以GPU的各个线程单元来对每个波数ky中的沿y方向傅里叶变换积分进行并行计算；
所述的步骤7中，基于二次插值函数的一维有限单元法构成对角方程组Ku＝P，其中K为五对角阵的系数矩阵；u为求解空间-波数域重力位未知量，P为右端项；从而形成Nkx×Nky个五对角方程组；然后利用二维波数域下重力位与重力场、梯度张量之间的关系，采用OpenMP加速以追赶法求解方程组得到二维波数域下的重力位：
其中指的是二维波数域下的重力位，/>为二维波数域下的重力场，/>和/>分别为x、y和z三个方向上的分量，/>为二维波数域下的重力梯度张量，/> 为9个梯度张量分量，i为虚数单位；
所述步骤9中，利用步骤8中计算出来的积分系数，对波数域的重力异常场以及梯度张量沿着水平方向进行二维全息傅里叶反变换：
F(kx,ky)表示波数域的重力异常场、梯度张量，沿kx方向进行一维积分得到Fx(x,ky):
计算过程中利用GPU进行多线程加速，以GPU的各个线程单元来对每个空间节点x下的沿kx方向傅里叶变换积分进行并行计算；
然后对Fx(x,ky)沿ky方向进行一维积分得到f(x,y)，f(x,y)表示空间域重力异常场以及梯度张量：
计算过程中利用GPU进行多线程加速，以GPU的各个线程单元来对每个空间节点y下的沿ky方向傅里叶变换积分进行并行计算。
2.根据权利要求1所述的方法，其特征在于，所述的步骤1中，重力位与空间域密度满足泊松方程为：
其中，为拉普拉斯算子，U为空间域重力位，(x,y,z)为空间观测点坐标，π为圆周率，G为万有引力常量，ρ为空间域密度；
然后沿水平方向做二维傅里叶正变换，得到二维波数域下的一维控制方程：
其中，表示偏导，/>是空间-波数域重力位，/>是空间-波数混合域密度，kx和ky分别是x和y方向傅里叶变换后的波数，k为中间变量，/>
且在对基于一维控制方程表示的一系列平面波进行求解时，应结合上下边界条件：
其中zmin代表上边界，zmax代表下边界。
3.根据权利要求1所述的方法，其特征在于，所述的步骤2中，基于x，y和z方向对目标区域进行离散以实现剖分时，是通过均匀剖分或不均匀剖分来将目标区域剖分为多个空间单元，从而形成多个离散节点，离散节点的数量即空间单元的数量+1。
4.根据权利要求1所述的方法，其特征在于，所述的步骤3中，采样规则为：在频谱能量大且变化剧烈的区间加密采样，在频谱能量小且变化缓慢的区间稀疏采样，其中区间为波数(kx、ky)的范围区间；其中波数的基本单位即基频的表达式为：
其中tmax为目标区域x或y方向计算范围的最大值，tmin为空间域x或y方向计算范围的最小值；当空间域x、y方向上的采样点数为Nx和Ny、采样间隔为Δx和Δy时，根据采样规则求得x和y方向对应的波数kx和ky，波数节点数量分别为Nkx、Nky，π为圆周率。
5.根据权利要求1所述的方法，其特征在于，所述的步骤4中，给定每个节点的剩余密度为ρ，其中异常体剩余密度与周围剩余密度不同。
6.根据权利要求1所述的方法，其特征在于，所述的步骤5中，通过下式计算积分系数：
其中I(xj1,kx)，I(xj2,kx)，I(xj3,kx)分别为沿x方向第j个空间单元上三个节点的傅里叶正变换积分系数，J(yj1,ky)，J(yj2,ky)，J(yj3,ky)分别为沿y方向上第j个空间单元的傅里叶正变换节点积分系数，e表示自然对数的底数，i为虚数单位。
7.根据权利要求1所述的方法，其特征在于，所述的步骤8中，通过下式计算积分系数，
其中分别为沿kx方向第j个空间单元上三个节点的傅里叶反变换积分系数，/>分别为沿ky方向上第j个空间单元的傅里叶反变换节点积分系数，e表示自然对数的底数，i为虚数单位。
CN202311184544.3A 2023-09-14 2023-09-14 三维重力场全息数值模拟及cpu-gpu加速方法 Active CN116911146B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202311184544.3A CN116911146B (zh) 2023-09-14 2023-09-14 三维重力场全息数值模拟及cpu-gpu加速方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202311184544.3A CN116911146B (zh) 2023-09-14 2023-09-14 三维重力场全息数值模拟及cpu-gpu加速方法

Publications (2)

Publication Number Publication Date
CN116911146A CN116911146A (zh) 2023-10-20
CN116911146B true CN116911146B (zh) 2024-01-19

Family

ID=88360728

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202311184544.3A Active CN116911146B (zh) 2023-09-14 2023-09-14 三维重力场全息数值模拟及cpu-gpu加速方法

Country Status (1)

Country Link
CN (1) CN116911146B (zh)

Citations (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US3952186A (en) * 1975-02-10 1976-04-20 The United States Of America As Represented By The Secretary Of The Navy Apparatus for the generation of a two-dimensional discrete fourier transform
JP2000200261A (ja) * 1998-12-29 2000-07-18 Hitachi Ltd フ―リエ変換方法、シミュレ―ション方法およびプログラム記録媒体
WO2007052631A1 (ja) * 2005-10-31 2007-05-10 Bycen Inc. 歩様バランス定量化方法および歩様バランス定量化装置
CN105334542A (zh) * 2015-10-23 2016-02-17 中南大学 任意密度分布复杂地质体重力场快速、高精度正演方法
KR20160107702A (ko) * 2015-03-05 2016-09-19 인하대학교 산학협력단 파수-공간-시간 영역 공통중간점 모음자료 합성을 위한 효율적인 음향 파동 방정식 유한차분법 모델링 방법
CN107577641A (zh) * 2017-08-21 2018-01-12 吉林大学 一种基于gpu并行的重力梯度张量数据快速密度反演方法
CN112800657A (zh) * 2021-04-15 2021-05-14 中南大学 基于复杂地形的重力场数值模拟方法、装置和计算机设备
CN113420487A (zh) * 2021-08-25 2021-09-21 中南大学 一种三维重力梯度张量数值模拟方法、装置、设备和介质
CN113591030A (zh) * 2021-08-17 2021-11-02 东北大学 基于多gpu的重力梯度数据灵敏度矩阵压缩及调用方法
CN114490011A (zh) * 2020-11-12 2022-05-13 上海交通大学 N体模拟在异构架构的并行加速实现方法
CN114611062A (zh) * 2022-02-23 2022-06-10 中国地质科学院地球物理地球化学勘查研究所 三维重力快速反演优化方法、系统、存储介质和电子设备
CN116258042A (zh) * 2023-01-31 2023-06-13 重庆励颐拓软件有限公司 一种基于ddm的大规模传热异构并行仿真方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP5127242B2 (ja) * 2007-01-19 2013-01-23 任天堂株式会社 加速度データ処理プログラムおよびゲームプログラム
US8700372B2 (en) * 2011-03-10 2014-04-15 Schlumberger Technology Corporation Method for 3-D gravity forward modeling and inversion in the wavenumber domain

Patent Citations (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US3952186A (en) * 1975-02-10 1976-04-20 The United States Of America As Represented By The Secretary Of The Navy Apparatus for the generation of a two-dimensional discrete fourier transform
JP2000200261A (ja) * 1998-12-29 2000-07-18 Hitachi Ltd フ―リエ変換方法、シミュレ―ション方法およびプログラム記録媒体
WO2007052631A1 (ja) * 2005-10-31 2007-05-10 Bycen Inc. 歩様バランス定量化方法および歩様バランス定量化装置
KR20160107702A (ko) * 2015-03-05 2016-09-19 인하대학교 산학협력단 파수-공간-시간 영역 공통중간점 모음자료 합성을 위한 효율적인 음향 파동 방정식 유한차분법 모델링 방법
CN105334542A (zh) * 2015-10-23 2016-02-17 中南大学 任意密度分布复杂地质体重力场快速、高精度正演方法
CN107577641A (zh) * 2017-08-21 2018-01-12 吉林大学 一种基于gpu并行的重力梯度张量数据快速密度反演方法
CN114490011A (zh) * 2020-11-12 2022-05-13 上海交通大学 N体模拟在异构架构的并行加速实现方法
CN112800657A (zh) * 2021-04-15 2021-05-14 中南大学 基于复杂地形的重力场数值模拟方法、装置和计算机设备
CN113591030A (zh) * 2021-08-17 2021-11-02 东北大学 基于多gpu的重力梯度数据灵敏度矩阵压缩及调用方法
CN113420487A (zh) * 2021-08-25 2021-09-21 中南大学 一种三维重力梯度张量数值模拟方法、装置、设备和介质
CN114611062A (zh) * 2022-02-23 2022-06-10 中国地质科学院地球物理地球化学勘查研究所 三维重力快速反演优化方法、系统、存储介质和电子设备
CN116258042A (zh) * 2023-01-31 2023-06-13 重庆励颐拓软件有限公司 一种基于ddm的大规模传热异构并行仿真方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
The forward modeling of 3D gravity and magnetic potential fields in space-wavenumber domains base on integral method;Shikun Dai;Google;全文 *
基于谱元-简正振型耦合方法的核幔边界D″区地震波波形模拟方法研究;赵明;赵亮;Yann Capdeville;;地球物理学报(04);全文 *
空间—波数混合域二度体重力异常正演;戴世坤 等;《石油地球物理勘探》;第1383-1389卷(第6期);全文 *
重力异常场空间-波数混合域三维数值模拟;戴世坤;陈轻蕊;李昆;凌嘉宣;;地球物理学报(05);全文 *

Also Published As

Publication number Publication date
CN116911146A (zh) 2023-10-20

Similar Documents

Publication Publication Date Title
Biancolini Fast radial basis functions for engineering applications
CN106777598B (zh) 任意磁化率分布复杂磁性体磁场梯度张量数值模拟方法
CN109375280B (zh) 一种球坐标系下重力场快速高精度正演方法
CN106646645B (zh) 一种重力正演加速方法
CN110133644B (zh) 基于插值尺度函数法的探地雷达三维正演方法
Zhao et al. High-accuracy 3D Fourier forward modeling of gravity field based on the Gauss-FFT technique
CN108197389A (zh) 二维强磁性体磁场的快速、高精度数值模拟方法
CN108984939A (zh) 基于3D Gauss-FFT的三维重力场正演方法
Henneböhl et al. Spatial interpolation in massively parallel computing environments
Dai et al. Three-dimensional numerical modeling of gravity anomalies based on Poisson equation in space-wavenumber mixed domain
CN115292973A (zh) 一种任意采样的空间波数域三维磁场数值模拟方法及系统
Wang et al. Improved preconditioned conjugate gradient algorithm and application in 3D inversion of gravity-gradiometry data
Ouyang et al. Iterative magnetic forward modeling for high susceptibility based on integral equation and Gauss-fast Fourier transform
Hou et al. Fast density inversion solution for full tensor gravity gradiometry data
CN107942399A (zh) 一种大距离位场向上延拓计算方法
CN116911146B (zh) 三维重力场全息数值模拟及cpu-gpu加速方法
CN107748834B (zh) 一种计算起伏观测面磁场的快速、高精度数值模拟方法
Li et al. Fast 3D forward modeling of the magnetic field and gradient tensor on an undulated surface
Nemitz et al. Topological sensitivity and FMM-accelerated BEM applied to 3D acoustic inverse scattering
Wang et al. Efficient 2D modeling of magnetic anomalies using NUFFT in the Fourier domain
Qiang et al. A fast forward algorithm for three-dimensional magnetic anomaly on undulating terrain
Saucedo-Zendejo A novel meshfree approach based on the finite pointset method for linear elasticity problems
Shutyaev et al. Numerical solution of the problem of variational data assimilation to restore heat fluxes and initial state for the ocean thermodynamics model
CN115640720A (zh) 一种基于距离控制网格加密的自引力仿真方法
Lan et al. High-precision 3D reconstruction of multiple magnetic targets based on center weighting method

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