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

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

Info

Publication number
CN116911146A
CN116911146A CN202311184544.3A CN202311184544A CN116911146A CN 116911146 A CN116911146 A CN 116911146A CN 202311184544 A CN202311184544 A CN 202311184544A CN 116911146 A CN116911146 A CN 116911146A
Authority
CN
China
Prior art keywords
dimensional
domain
wave number
thread
gravity
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.)
Granted
Application number
CN202311184544.3A
Other languages
English (en)
Other versions
CN116911146B (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 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方向的剖分计算波数kx和ky。
4）根据计算模型参数设置离散节点上的剩余密度。
5）计算基于形函数插值的二维全息傅里叶正变换需要的积分系数。
6）基于步骤5）的结果将目标区域的剩余密度进行水平方向全息傅里叶正变换，并基于GPU多线程加速，得到二维波数域剩余密度。
7）将二维波数域剩余密度代入二维波数域下的一维控制方程，采用基于二次插值的一维有限单元法，并加载上、下边界条件，得到五对角方程组。再采用CPU上的OpenMP加速以追赶法求解方程组，得到二维波数域下的重力位。最后利用空间-波数域重力位与重力场、梯度张量之间的关系，得到空间-波数域重力场以及梯度张量。
8）计算基于形函数插值的二维全息傅里叶反变换需要的积分系数。
9）将空间-波数域重力场进行二维全息傅里叶反变换，并基于GPU多线程加速，得到空间域重力场以及梯度张量。
所述的方法，所述的步骤1）中，重力位与空间域密度满足泊松方程为：
其中，为拉普拉斯算子，U为空间域重力位，(x,y,z)为空间观测点坐标，π为圆周 率，G为万有引力常量，ρ为空间域密度。
然后沿水平方向做二维傅里叶正变换，得到二维波数域下的一维控制方程：
其中，表示偏导，是空间-波数域重力位，是空间-波数混合域 密度，和分别是x和y方向傅里叶变换后的波数，k为中间变量，。
且在对基于一维控制方程表示的一系列平面波进行求解时，应结合上下边界条件：
其中代表上边界，代表下边界。
所述的方法，所述的步骤2）中，基于x，y和z方向对目标区域进行离散以实现剖分时，是通过均匀剖分或不均匀剖分来将目标区域剖分为多个空间单元，从而形成多个离散节点，离散节点的数量即空间单元的数量+1。
所述的方法，所述的步骤3）中，采样规则为：在频谱能量大且变化剧烈的区间加密 采样，在频谱能量小且变化缓慢的区间稀疏采样，其中区间为波数（、）的范围区间。其 中波数的基本单位即基频的表达式为：
其中tmax为目标区域x或y方向计算范围的最大值，tmin为空间域x或y方向计算范围 的最小值。当空间域x、y方向上的采样点数为Nx和Ny、采样间隔为Δx和Δy时，根据采样规则 求得x和y方向对应的波数kx和ky，波数节点数量分别为、，π为圆周率。
所述的方法，所述的步骤4）中，给定每个节点的剩余密度为ρ，其中异常体剩余密度与周围剩余密度不同。
所述的方法，所述的步骤5）中，通过下式计算积分系数：
其中，，分别为沿x方向第j个空间单元上三个节点的傅里 叶正变换积分系数，，，分别为沿y方向上第j个空间单元的傅里 叶正变换节点积分系数，e表示自然对数的底数，i为虚数单位。
所述的方法，所述的步骤6）中，是根据步骤5得到的积分系数，对目标区域的剩余 密度沿水平方向进行二维全息傅里叶正变换：
首先根据下式对剩余密度沿x方向进行一维积分计算得到：
其中表示在不同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方向上的波数个数；
然后根据下式对沿y方向进行一维积分计算得到波数域下的：
其中表示在不同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（k x ，k y ）表示波数域的重力异常场、梯度张量，沿kx方向进行一维积分得到F x （x，k y ）:
计算过程中利用GPU进行多线程加速，以GPU的各个线程单元来对每个空间节点x下的沿kx方向傅里叶变换积分进行并行计算。
然后对F（k x ，k y ）沿ky方向进行一维积分得到，F x （x，k y ）表示空间域重力异常 场以及梯度张量：
计算过程中利用GPU进行多线程加速，以GPU的各个线程单元来对每个空间节点y下的沿ky方向傅里叶变换积分进行并行计算。
本发明的技术效果在于，本发明实现了基于二次插值形函数任意采样傅里叶变换的空间-波数域重力异常数值模拟算法，并采用CPU-GPU进行加速，实现了三维重力异常高效、高精度数值模拟。利用水平方向二维傅里叶变换把空间域引力位满足的三维偏微分方程转化为不同波数满足的一维常微分方程，将三维问题转化为一维常微分问题求解，避免了大型系数矩阵方程组的求解，减少了计算量及存储需求，而且节点规模越大，效果越明显，不同波数之间常微分方程相互独立，具有高度并行性。常微分方程上下边界问题具有严格的边界条件，与物理边界吻合，解决上下边界问题。水平方向二维傅里叶变换采用任意采样傅里叶变换，与重力场水平方向左右、前后边界吻合，解决重力场水平方向左右、前后边界问题，数值模拟结果与重力场解析解高度吻合。基于空间波数域算法以及任意傅里叶方法的并行性，在CPU上并行求解五对角方程组，GPU上计算正反二维傅里叶变换，实现CPU-GPU并行加速方案，兼顾计算效率与计算精度。该方法适用于任意复杂模型三维重力异常高效、高精度数值模拟，为大规模重力实测数据的高效、精细反演提供重要技术支持。
附图说明
图1为本发明实施例中棱柱体模型的平面图。
图2为本发明实施例中重力异常场在地面z = 0 m数值解的平面等值线图。
图3为本发明实施例中重力异常场在地面z = 0 m解析解的平面等值线图。
具体实施方式
本实施例所提供的基于任意采样傅里叶变换的空间-波数域重力异常数值模拟及CPU-GPU加速方法，包括以下步骤：
1）引入重力位满足的泊松方程，对该方程进行水平方向二维傅里叶变换，得到二维波数域下的一维控制方程组。
这里的重力位与空间域密度满足泊松方程(Blakely R J. 1996. Potential theory in gravity and magnetic applications: London, Cambridge University Press.)：
，
其中，为拉普拉斯算子，U为空间域重力位，(x,y,z)为空间观测点坐标，π为圆周 率，G为万有引力常量，值为，ρ为空间域密度。再对式(1)沿水平方向做二维 傅里叶正变换，得到式（2）空间-波数域常微分方程即二维波数域下的一维控制方程：
.
其中，表示偏导，是空间-波数域重力位，是空间-波数混 合域密度，和分别是x和y方向傅里叶变换后的波数，k为中间变量，。
这里在对基于一维控制方程表示的一系列平面波进行求解时，应结合上下边界条件：
其中代表上边界，代表下边界。从而使边界条件与重力位物理边界条件 一致，空间-波数域重力异常上下边界问题得到准确表示。
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，波数节点数量分别为、。
4）根据计算模型参数设置离散节点上的剩余密度。本实施例中给定每个节点的剩余密度为ρ，其中异常体剩余密度与周围剩余密度不同。
5）计算基于形函数插值的二维任意傅里叶正变换需要的积分系数。
对于基于形函数的任意傅里叶正变换，其核心思想是将二维傅里叶正变换转化为两个一维傅里叶正变换，并把一维傅里叶正变换离散为多个单元积分累加和，而离散单元中的原函数可以用二次插值形函数拟合，然后求出单元积分的解析表达式，最终表示为：
其中，，为沿x方向第j个单元上三个节点的傅里叶变换积 分系数，，，分别为沿y方向上第j个单元的傅里叶变换节点积分系 数，该积分系数具体为：
其中e表示自然对数的底数，i为虚数单位。根据式(6)和(7)，x、y、kx以及ky确定以后，即可求出任意采样傅里叶变换沿两个方向的积分系数，并且当x、y方向剖分固定时，单元形函数不变，积分系数也不会改变，因此可以将积分系数提前计算出来并存储，以重复利用，从而减少冗余计算，提高计算效率。
6）将目标区域的剩余密度进行水平方向任意傅里叶正变换，采用GPU多线程加速，得到二维波数域剩余密度。
本步骤中是利用上一步骤中计算出来的积分系数，对目标区域的剩余密度 沿水平方向进行二维任意傅里叶正变换，第一步根据式(4)对剩余密度沿x方向进行 一维积分计算得到，具体的计算公式为：
其中表示在不同z值下的。
由于每个波数下的积分计算相互独立，因此沿x方向傅里叶变换积分计算可以 并行执行。本实施例中利用GPU众多的线程单元来进行并行计算，这里线程单元采用二维线 程块和二维线程网格来组织，GPU的线程单元间通过线程索引( threadIdx%x, threadIdx% y )和线程网格索引( blockIdx%x, blockIdx%y )定位自己在线程网格中的索引，这样来 使每个GPU的线程单元准确找到线程数据对应的逻辑位置。一个线程块中有a×b个线程，线 程网格中有m×n个线程块，线程块中的一个线程单元负责一个波数计算：即Nx次乘法和加 法运算。因此有：，，Nkx表示kx方向上下的波数个数。
第二步根据式(5)对沿着y方向进行一维积分计算得到波数域下的，计算过程类似沿x方向，同样采用GPU多线程并行加速，即根据下式对 沿y方向进行一维积分计算得到波数域下的：
其中表示在不同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方向第j个空间单元上三个节点的 傅里叶反变换积分系数，，，分别为沿ky方向上第j个空间单元的 傅里叶反变换节点积分系数，e表示自然对数的底数，i为虚数单位。
9）将空间-波数域重力场进行二维任意傅里叶反变换，并采用GPU加速，得到空间域重力场以及梯度张量。
利用步骤8中计算出来的积分系数，对波数域的重力异常场以及梯度张量沿着水平方向进行二维任意傅里叶反变换，计算加速过程类似正变换过程。
F（k x ，k y ）表示波数域的重力异常场、梯度张量，沿kx方向进行一维积分得到F x （x，k y ）:
计算过程中利用GPU进行多线程加速，以GPU的各个线程单元来对每个空间节点x下的沿kx方向傅里叶变换积分进行并行计算；
然后对F x （x，k y ）沿ky方向进行一维积分得到，表示空间域重力异常场以 及梯度张量：
计算过程中利用GPU进行多线程加速，以GPU的各个线程单元来对每个空间节点y下的沿ky方向傅里叶变换积分进行并行计算。
基于以上步骤，进行数值模拟验证过程如下：
设定如图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×2 km×2 km，异常体剩余密度=2000，=2000。根据步骤3确定波数选取范围 为~，波数选取个数为253×253。用解析解的计算结果为参照，验证方法的正 确性。
图2、图3分别为重力异常场在地面z = 0 m数值解以及解析解的平面等值线图。由 图可知，数值解和解析解等值线图形态相似，场值以及梯度张量吻合得很好，且在z = 0 m 平面、、、、、、、以及9个分量的相对均方根误差分别为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，128 GBRAM. GPU显卡型号NVIDIA TITAN RTX，CUDA版本号11.6，运行于Linux编译环境。本实施例的计算节点数为201×201×101=15813251，计算总时间为0.7s，占用CPU+GPU内存为819+363MB，由此可见本发明计算速度快，占用内存小，具有很高的计算效率。

Claims (10)

1.一种三维重力场全息数值模拟及CPU-GPU加速方法，其特征在于，包括以下步骤：
1）对重力位与空间域密度满足泊松方程进行水平方向二维傅里叶变换，得到二维波数域下的一维控制方程；
2）在直角坐标系中建立能够包围三维剩余密度异常体的最小长方体目标区域，并基于x，y和z方向对目标区域进行离散以实现剖分并形成离散节点；其中三维剩余密度异常体是根据计算需要建立的；
3）利用基于形函数插值的积分区间为无限的全息傅里叶变换的采样规则，根据x和y方向的剖分计算波数kx和ky；
4）根据计算模型参数设置离散节点上的剩余密度；
5）计算基于形函数插值的二维全息傅里叶正变换需要的积分系数；
6）基于步骤5）的结果将目标区域的剩余密度进行水平方向全息傅里叶正变换，并基于GPU多线程加速，得到二维波数域剩余密度；
7）将二维波数域剩余密度代入二维波数域下的一维控制方程，采用基于二次插值的一维有限单元法，并加载上、下边界条件，得到五对角方程组；再采用CPU上的OpenMP加速以追赶法求解方程组，得到二维波数域下的重力位；最后利用空间-波数域重力位与重力场、梯度张量之间的关系，得到空间-波数域重力场以及梯度张量；
8）计算基于形函数插值的二维全息傅里叶反变换需要的积分系数；
9）将空间-波数域重力场进行二维全息傅里叶反变换，并基于GPU多线程加速，得到空间域重力场以及梯度张量。
2.根据权利要求1所述的方法，其特征在于，所述的步骤1）中，重力位与空间域密度满足泊松方程为：
；
其中，为拉普拉斯算子，U为空间域重力位，(x,y,z)为空间观测点坐标，π为圆周率，G为万有引力常量，ρ为空间域密度；
然后沿水平方向做二维傅里叶正变换，得到二维波数域下的一维控制方程：
；
其中，表示偏导，/>是空间-波数域重力位，/>是空间-波数混合域密度，kx和ky分别是x和y方向傅里叶变换后的波数，k为中间变量，/>；
且在对基于一维控制方程表示的一系列平面波进行求解时，应结合上下边界条件：
；
；
其中代表上边界，/>代表下边界。
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，波数节点数量分别为N kx 、N ky ，π为圆周率。
5.根据权利要求1所述的方法，其特征在于，所述的步骤4）中，给定每个节点的剩余密度为ρ，其中异常体剩余密度与周围剩余密度不同。
6.根据权利要求1所述的方法，其特征在于，所述的步骤5）中，通过下式计算积分系数：
；
；
其中，/>，/>分别为沿x方向第j个空间单元上三个节点的傅里叶正变换积分系数，/>，/>，/>分别为沿y方向上第j个空间单元的傅里叶正变换节点积分系数，e表示自然对数的底数，i为虚数单位。
7.根据权利要求1所述的方法，其特征在于，所述的步骤6）中，是根据步骤5得到的积分系数，对目标区域的剩余密度沿水平方向进行二维全息傅里叶正变换：
首先根据下式对剩余密度沿x方向进行一维积分计算得到/>：
；
其中表示在不同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方向上的波数个数；
然后根据下式对沿y方向进行一维积分计算得到波数域下的/>：
；
其中表示在不同z值下的/>；
计算过程中同样利用GPU进行多线程加速，以GPU的各个线程单元来对每个波数ky中的沿y方向傅里叶变换积分进行并行计算。
8.根据权利要求1所述的方法，其特征在于，所述的步骤7）中，基于二次插值函数的一维有限单元法构成对角方程组Ku=P，其中K为五对角阵的系数矩阵；u为求解空间-波数域重力位未知量，P为右端项；从而形成Nkx×Nky个五对角方程组；然后利用二维波数域下重力位与重力场、梯度张量之间的关系，采用OpenMP加速追赶法求解方程组得到二维波数域下的重力位：
；
；
其中指的是二维波数域下的重力位，/>为二维波数域下的重力场，/>、/>和/>分别为x、y和z三个方向上的分量，/>为二维波数域下的重力梯度张量，/> 为9个梯度张量分量，i为虚数单位。
9.根据权利要求1所述的方法，其特征在于，所述的步骤8）中，通过下式计算积分系数，
；
；
其中，/>，/>分别为沿kx方向第j个空间单元上三个节点的傅里叶反变换积分系数，/>，/>，/>分别为沿ky方向上第j个空间单元的傅里叶反变换节点积分系数，e表示自然对数的底数，i为虚数单位。
10.根据权利要求1所述的方法，其特征在于，所述步骤9）中，利用步骤8中计算出来的积分系数，对波数域的重力异常场以及梯度张量沿着水平方向进行二维全息傅里叶反变换：
F（k x ，k y ）表示波数域的重力异常场、梯度张量，沿kx方向进行一维积分得到F x （x，k y ）:
；
计算过程中利用GPU进行多线程加速，以GPU的各个线程单元来对每个空间节点x下的沿kx方向傅里叶变换积分进行并行计算；
然后对F x （x，k y ）沿ky方向进行一维积分得到f（x，y），f（x，y）表示空间域重力异常场以及梯度张量：
；
计算过程中利用GPU进行多线程加速，以GPU的各个线程单元来对每个空间节点y下的沿ky方向傅里叶变换积分进行并行计算。
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 true CN116911146A (zh) 2023-10-20
CN116911146B 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 (14)

* 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. 歩様バランス定量化方法および歩様バランス定量化装置
US20080177497A1 (en) * 2007-01-19 2008-07-24 Nintendo Co., Ltd. Storage medium having acceleration data processing program stored thereon, storage medium having game program stored thereon, and acceleration data processing apparatus
US20120232871A1 (en) * 2011-03-10 2012-09-13 Ivan Priezzhev Method for 3-d gravity forward modeling and inversion in the wavenumber domain
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的大规模传热异构并行仿真方法

Patent Citations (14)

* 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. 歩様バランス定量化方法および歩様バランス定量化装置
US20080177497A1 (en) * 2007-01-19 2008-07-24 Nintendo Co., Ltd. Storage medium having acceleration data processing program stored thereon, storage medium having game program stored thereon, and acceleration data processing apparatus
US20120232871A1 (en) * 2011-03-10 2012-09-13 Ivan Priezzhev Method for 3-d gravity forward modeling and inversion in the wavenumber domain
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
SHIKUN DAI: "The forward modeling of 3D gravity and magnetic potential fields in space-wavenumber domains base on integral method", GOOGLE *
戴世坤 等: "空间—波数混合域二度体重力异常正演", 《石油地球物理勘探》, vol. 1383, no. 6 *
戴世坤;陈轻蕊;李昆;凌嘉宣;: "重力异常场空间-波数混合域三维数值模拟", 地球物理学报, no. 05 *
赵明;赵亮;YANN CAPDEVILLE;: "基于谱元-简正振型耦合方法的核幔边界D″区地震波波形模拟方法研究", 地球物理学报, no. 04 *

Also Published As

Publication number Publication date
CN116911146B (zh) 2024-01-19

Similar Documents

Publication Publication Date Title
Biancolini Fast radial basis functions for engineering applications
Ren et al. Fast 3‐D large‐scale gravity and magnetic modeling using unstructured grids and an adaptive multilevel fast multipole method
CN109375280B (zh) 一种球坐标系下重力场快速高精度正演方法
CN106777598B (zh) 任意磁化率分布复杂磁性体磁场梯度张量数值模拟方法
CN106646645B (zh) 一种重力正演加速方法
Martyshko et al. Studying the structural features of the lithospheric magnetic and gravity fields with the use of parallel algorithms
Ren et al. Gravity gradient tensor of arbitrary 3D polyhedral bodies with up to third-order polynomial horizontal and vertical mass contrasts
Taube et al. A high‐order discontinuous Galerkin method with time‐accurate local time stepping for the Maxwell equations
Zhao et al. High-accuracy 3D Fourier forward modeling of gravity field based on the Gauss-FFT technique
CN103278848A (zh) 基于mpi并行预条件迭代的地震成像正演方法
CN108984939A (zh) 基于3D Gauss-FFT的三维重力场正演方法
Wang et al. Improved preconditioned conjugate gradient algorithm and application in 3D inversion of gravity-gradiometry data
Dai et al. Three-dimensional numerical modeling of gravity anomalies based on Poisson equation in space-wavenumber mixed domain
CN115292973A (zh) 一种任意采样的空间波数域三维磁场数值模拟方法及系统
Raevskiy et al. On the solution of inverse problems of gravimetry by the modified method of S-approximations
CN107942399A (zh) 一种大距离位场向上延拓计算方法
Dai et al. The forward modeling of 3D gravity and magnetic potential fields in space-wavenumber domains based on an integral method
CN116911146B (zh) 三维重力场全息数值模拟及cpu-gpu加速方法
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
Qiang et al. A fast forward algorithm for three-dimensional magnetic anomaly on undulating terrain
Wang et al. Efficient 2D modeling of magnetic anomalies using NUFFT in the Fourier domain
Yin et al. A fast 3D gravity forward algorithm based on circular convolution
CN115640720A (zh) 一种基于距离控制网格加密的自引力仿真方法
Bucha et al. Cap integration in spectral gravity forward modelling up to the full gravity tensor

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