CN113591030B - 基于多gpu的重力梯度数据灵敏度矩阵压缩及调用方法 - Google Patents
基于多gpu的重力梯度数据灵敏度矩阵压缩及调用方法 Download PDFInfo
- Publication number
- CN113591030B CN113591030B CN202110941656.3A CN202110941656A CN113591030B CN 113591030 B CN113591030 B CN 113591030B CN 202110941656 A CN202110941656 A CN 202110941656A CN 113591030 B CN113591030 B CN 113591030B
- Authority
- CN
- China
- Prior art keywords
- gpu
- index
- sensitivity matrix
- layer
- gradient data
- 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
Links
- 239000011159 matrix material Substances 0.000 title claims abstract description 103
- 230000035945 sensitivity Effects 0.000 title claims abstract description 77
- 230000005484 gravity Effects 0.000 title claims abstract description 43
- 238000000034 method Methods 0.000 title claims abstract description 43
- 230000006835 compression Effects 0.000 title claims abstract description 31
- 238000007906 compression Methods 0.000 title claims abstract description 31
- 238000004364 calculation method Methods 0.000 claims abstract description 26
- 238000012795 verification Methods 0.000 abstract 1
- 238000012545 processing Methods 0.000 description 4
- 238000011160 research Methods 0.000 description 4
- 238000012360 testing method Methods 0.000 description 4
- 230000001133 acceleration Effects 0.000 description 3
- 230000008451 emotion Effects 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 3
- 238000004891 communication Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 230000004888 barrier function Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000007667 floating Methods 0.000 description 1
- 230000006870 function Effects 0.000 description 1
- 230000002068 genetic effect Effects 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000005065 mining Methods 0.000 description 1
- 230000010287 polarization Effects 0.000 description 1
- 230000003595 spectral effect Effects 0.000 description 1
- 238000010998 test method Methods 0.000 description 1
Classifications
-
- 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
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
-
- 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
- G06F9/5044—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 considering hardware capabilities
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Software Systems (AREA)
- Mathematical Physics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Data Mining & Analysis (AREA)
- General Engineering & Computer Science (AREA)
- Computational Mathematics (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Computing Systems (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明提供一种基于多GPU的重力梯度数据灵敏度矩阵压缩及调用方法，涉及地球物理正反演技术领域；首先建立重力梯度数据各分量的灵敏度矩阵的等价计算公式，然后计算各GPU中分得灵敏度矩阵元素的范围和数量等参数，利用GPU并行计算压缩后的灵敏度矩阵，最后获取待参加计算的灵敏度矩阵元素索引，换算其在压缩矩阵中的索引并读取数值，从而实现大规模重力梯度数据不同梯度分量或全张量梯度数据的快速正反演；经验证，本发明能够有效提高重力梯度数据正反演的计算规模和大规模重力梯度数据联合反演的计算效率；还适用于重、磁异常等其他类型的位场数据正反演领域，也适用于单GPU的情况。
Description
技术领域
本发明涉及地球物理正反演技术领域，具体涉及一种基于多GPU的重力梯度数据灵敏度矩阵压缩及调用方法。
背景技术
重力梯度数据是一种高精度的地球物理数据，具有更高的信噪比。和应用重力异常相比，全张量梯度数据包含更多的地质信息，利用多分量梯度数据实现联合反演可获得更高分辨率的地下空间三维密度分布。位场数据正反演都需要灵敏度矩阵计算，然而少量的观测数据和地下划分网格层数也能产生大尺寸的灵敏度矩阵，严重降低了正反演计算效率，同时占用了大量的内存，致使普通计算设备无法完成正反演。
并行计算作为当今科技创新的重要手段已应用于不同的科研领域中。随着近20年的技术发展，目前基于集群机的MPI(Message Passing Interface)、基于CPU(CentralProcessing Unit)多线程的OpenMP(Open Multi-Processing)和基于GPU(GraphicsProcessing Unit)的CUDA(Computed Unified Device Architecture)等并行方法已较为广泛地应用于以文献[1-5]为代表的非震地球物理数据正反演研究中(其中，文献[1]为：陈召曦,孟小红,郭良辉等.基于GPU并行的重力,重力梯度三维正演快速计算及反演策略.地球物理学报,2012,55(12):4069-4077.；文献[2]为：胡祥云,李焱,杨文采等.大地电磁三维数据空间反演并行算法研究.地球物理学报,2012,55(12):3969-3978.；文献[3]为：M,Zhdanov M S.Massively parallel regularized 3D inversion of potential fieldson CPUs and GPUs.Computers&Geosciences,2014,62:80-87.；文献[4]为：Zhang Z Y,TanH D,Wang K P,et al.Two-dimensional inversion of spectral induced polarizationdata using MPI parallel algorithmin data space.Applied Geophysics,2016,13(1):13-24.；文献[5]为：Hou Z,Huang D,Wang E,et al.3D density inversion of gravitygradiometry data with a multilevel hybrid parallel algorithm.AppliedGeophysics,2019,16(2):141-152.)。现有的研究表明，并行计算能够有效地提高运算速度。其中，使用GPU的方法能够减少或避免在CPU集群节点间出现的大量通信开销等问题，对属于数据并行模型的算法(例如：位场数据正反演等)具有更好的加速效果。但受带宽等硬件属性限制，目前利用单GPU的方法在处理复杂迭代反演时会遇到通信开销增大、计算效率降低等情况；而采用多GPU则可以减少每个GPU的数据吞吐量、提高算法并行度。因此，通过多GPU开展大规模数据正反演的运算加速研究是十分必要的。
为了处理更大规模数据的反演问题，除了开展加速算法研究，还需要解决灵敏度矩阵内存占用的问题。硬件技术和矩阵压缩算法是目前解决该问题的两个重要手段。其中，硬件技术主要指借助并行计算中集群节点内存和GPU显存等硬件，获得更大的存储空间；矩阵压缩算法主要包括以文献[6-8]为代表的小波变换法(其中，文献[6]为：Li Y,OldenburgD W.Fast inversion of large-scale magnetic data using wavelet transforms anda logarithmic barrier method.Geophysical Journal International,2003,152(2):251-265.；文献[7]为：Martin,R.,Monteiller,V.,Komatitsch,D.,et al.Gravityinversion using wavelet-based compression on parallel hybrid CPU/GPU systems:application to southwest Ghana.Geophysical Journal International,2013,195,1594-1619.；文献[8]为：Sun,S.,Yin,C.,Gao,X.,Liu,Y.,Ren,X.Compressed SensingForward Modeling and Multiscale Gravity Inversion Based on WaveletTransform.Applied Geophysics,2018,15(2),342-352.)、以文献[9-10]为代表的几何格架等效压缩存储技术(其中，文献[9]为：姚长利,郝天珧,管志宁等.重磁遗传算法三维反演中高速计算及有效存储方法技术.地球物理学报,2003,46(2):252-258.；文献[10]为：姚长利,郑元满,张聿文.重磁异常三维物性反演随机子域法方法技术.地球物理学报,2007,50(5):1576-1583.)和文献[11-13]等研究中使用的其他减少灵敏度矩阵内存占用的算法(其中，文献[11]为：Portniaguine O,Zhdanov M S.3-D magnetic inversion with datacompression and image focusing.Geophysics,2002,67(5):1532-1541.；文献[12]为：秦朋波,黄大年.重力和重力梯度数据联合聚焦反演方法.地球物理学报,2016,59(6):2203-2224.；文献[13]为：Hou Z,Wei X,Huang D.Fast density inversion solution for fulltensor gravity gradiometry data.Pure and Applied Geophysics,2016,173,509-523.)。在并行计算框架内提供矩阵压缩算法的软件接口将有利于进一步提高反演的大数据计算能力。
因此，开展重力梯度数据灵敏度矩阵压缩及调用的并行方法研究是非常重要的，对于深部地质找矿具有实际意义。
发明内容
针对现有技术的不足，本发明提出一种基于多GPU的重力梯度数据灵敏度矩阵压缩及调用方法，包括：
步骤1：建立重力梯度数据各分量的灵敏度矩阵的等价计算公式；
步骤2：针对多GPU的并行运算，计算各GPU中分得灵敏度矩阵元素的范围、矩阵元素值的真实索引；
步骤3：计算压缩后的重力梯度数据各分量的灵敏度矩阵；
步骤4：当需要重力梯度数据各分量的灵敏度矩阵参与计算时，获取矩阵元素值的真实索引并换算成压缩后的灵敏度矩阵中的索引；
步骤5：根据换算后的索引读取压缩后的灵敏度矩阵中的数值，实现灵敏度矩阵的调用。
所述步骤1中等价计算公式表示为：
Gxx(i,j,m,n)＝Gxx(|i-m|+1,|j-n|+1,1,1) (1)
Gyy(i,j,m,n)＝Gyy(|i-m|+1,|j-n|+1,1,1) (2)
Gzz(i,j,m,n)＝Gzz(|i-m|+1,|j-n|+1,1,1) (3)
式中，Gxx(i,j,m,n)、Gyy(i,j,m,n)、Gzz(i,j,m,n)、Gxy(i,j,m,n)、Gxz(i,j,m,n)、Gyz(i,j,m,n)分别表示等价后的灵敏度矩阵，索引i、j、m、n表示测网中第i行第j列的测点索引与地下某一层中第m行第n列对应的网格索引。
所述步骤2包括：
步骤2.1：计算每个GPU分得的最后一个网格的索引last_m：
last_m＝base+local_m (7)
式中，base表示排序在当前GPU之前的各GPU中网格数量的和，用于将当前GPU中网格索引转换为真实索引；local_m表示各个GPU中分得的网格单元数量；
步骤2.2：计算当前GPU中分得的网格所在层数的范围：
first_layer＝base/N+1 (8)
last_layer＝(last_m-1)/N+1 (9)
式中，first_layer和last_layer分别表示分得网格所在的第一层和最后一层；N表示每层网格单元的数量；
步骤2.3：计算每个GPU分得的压缩矩阵大小：
size_mat＝(last_layer-first_layer+1)×N (10)
式中，size_mat表示每个GPU分得的压缩矩阵尺寸，即矩阵元素数量；
步骤2.4：计算用于将当前GPU中压缩灵敏度矩阵元素索引转换为真实压缩矩阵索引的参数base_mat：
base_mat＝(first_layer-1)×N (11)
其中，base_mat表示排序在当前GPU之前的各GPU中压缩矩阵元素数量的和。
所述步骤2还包括计算各GPU中分得的网格单元数量local_m：
local_m＝M/num_gpu (12)
式中，M为划分的网格单元总数量；num_gpu表示使用的GPU的数量；/表示求商运算，结果向下取整。
进一步地，当M不能被num_gpu整除时，则网格单元数量local_m表示为：
local_m＝(M-rem_m)/num_gpu+1 (13)
式中，rem_m表示余数。
所述步骤3中重力梯度数据各分量的灵敏度矩阵的计算公式为：
式中，Gxx、Gyy、Gzz、Gxy、Gxz、Gyz分别表示重力梯度数据Vxx、Vyy、Vzz、Vxy、Vxz、Vyz的灵敏度矩阵；xi＝x-ξi，yj＝y-ηj，zk＝z-ζk，μijk＝(-1)i+j+k；γ表示万有引力常数；x、y和z为任意观测点坐标；ξ、η和ζ分别为网格单元在X、Y、Z轴方向上的对顶点坐标。
所述步骤4包括：
步骤4.1：利用公式(20)换算每个GPU中分得的网格单元的列索引i_col：
i_col＝i_local+base (20)
式中，i_local表示各GPU中的索引；
步骤4.2：计算第i_col个网格单元所在层数以及在该层的行、列数：
i_x＝((i_col-1)％N)％nx+1 (21)
i_y＝((i_col-1)％N)/nx+1 (22)
i_z＝(i_col-1)/N+1 (23)
式中，i_x、i_y分别表示网格单元所处的列数、行数；％表示求余数的运算；nx为每条测线上的测点数；i_z表示网格单元所在层数；
步骤4.3：利用灵敏度矩阵元素的行索引计算测点的行、列索引：
obs_x＝(j_row-1)％nx+1 (24)
obs_y＝(j_row-1)/nx+1 (25)
式中，obs_x、obs_y分别为测点的列、行索引；j_row为灵敏度矩阵元素的行索引；
步骤4.4：计算当前GPU中压缩后的灵敏度矩阵元素索引：
p_i＝|obs_x-i_x|+1 (26)
p_j＝|obs_y-i_y|+1 (27)
index＝(i_z-first_layer)×N+(p_j-1)×nx+p_i (28)
式中，index为当前GPU中压缩后的灵敏度矩阵元素索引值；通过读取当前GPU中第index个压缩后的灵敏度矩阵元素值，完成对元素值的调用。
本发明的有益效果是：
本发明提出了一种基于多GPU的重力梯度数据灵敏度矩阵压缩及调用方法，能够有效解决重力梯度数据正反演中的灵敏度矩阵内存占用问题。和其他方法相比，具有压缩程度高、计算速度快等优点，适用于百万量级网格、甚至千万量级以上网格的大规模反演计算。本发明对于深部找矿中获取精细化三维密度分布、确定岩性与划分矿体、构造具有重要作用。
附图说明
图1为本发明中基于多GPU的重力梯度数据灵敏度矩阵压缩及调用方法流程图。
具体实施方式
下面结合附图和具体实施实例对发明做进一步说明。
如图1所示，一种基于多GPU的重力梯度数据灵敏度矩阵压缩及调用方法，包括：
步骤1：根据几何格架等效压缩存储技术，每一层地下网格中任意单元的几何格架都与该层中排序第1的单元的几何格架存在如下关系；即重力梯度数据各分量的灵敏度矩阵的等价计算公式为:
Gxx(i,j,m,n)＝Gxx(|i-m|+1,|j-n|+1,1,1) (1)
Gyy(i,j,m,n)＝Gyy(|i-m|+1,|j-n|+1,1,1) (2)
Gzz(i,j,m,n)＝Gzz(|i-m|+1,|j-n|+1,1,1) (3)
式中，Gxx(i,j,m,n)、Gyy(i,j,m,n)、Gzz(i,j,m,n)、Gxy(i,j,m,n)、Gxz(i,j,m,n)、Gyz(i,j,m,n)分别表示等价后的灵敏度矩阵，索引i、j、m、n表示测网中第i行第j列的测点索引与地下某一层中第m行第n列对应的网格索引。
对获取的重力梯度数据Vxx、Vyy、Vzz、Vxy、Vxz、Vyz首先进行网格化等常规处理。
步骤2：针对多GPU的并行运算，计算各GPU中分得灵敏度矩阵元素的范围、矩阵元素值的真实索引；包括：
步骤2.1：根据负载平衡和最大化并行度的原则，在正反演算法中，特别是在迭代反演中，将出现大量与网格单元总数量有关的计算，因此需要利用GPU对这些计算并行化，则计算各GPU中分得的网格单元数量local_m：
local_m＝M/num_gpu (12)
式中，M为划分的网格单元总数量；num_gpu表示使用的GPU的数量；/表示求商运算，结果向下取整；
当M不能被num_gpu整除时，则网格单元数量local_m表示为：
local_m＝(M-rem_m)/num_gpu+1 (13)
式中，rem_m表示余数；即前rem_m个GPU分得公式(13)结果的数量，其余GPU的数量比上式少1；
步骤2.2：假设划分的网格单元总数为M，由于与网格单元有关的、总个数为M的计算被分配给各GPU，相当于将M个网格进行分配，则每个GPU分得的最后一个网格的索引last_m表示为：
last_m＝base+local_m (7)
式中，base表示排序在当前GPU之前的各GPU中网格数量的和，用于将当前GPU中网格索引转换为真实索引；local_m表示各个GPU中分得的网格单元数量；
步骤2.3：计算当前GPU中分得的网格所在层数的范围：
first_layer＝base/N+1 (8)
last_layer＝(last_m-1)/N+1 (9)
式中，first_layer和last_layer分别表示分得网格所在的第一层和最后一层(层的排序是由浅及深的)；N表示测点数，亦为每层网格单元的数量；
步骤2.4：计算和分得网格相关的灵敏度矩阵的尺寸大小，即每个GPU分得的压缩矩阵大小：
size_mat＝(last_layer-first_layer+1)×N (10)
式中，size_mat表示每个GPU分得的压缩矩阵尺寸，即矩阵元素数量；
步骤2.5：计算用于将当前GPU中压缩灵敏度矩阵元素索引转换为真实压缩矩阵索引的参数base_mat：
base_mat＝(first_layer-1)×N (11)
其中，base_mat表示排序在当前GPU之前的各GPU中压缩矩阵元素数量的和；
步骤3：由步骤1可知，仅计算每一层地下网格中排序第1的单元的几何格架即可，计算压缩后的重力梯度数据各分量的灵敏度矩阵为：
式中，Gxx、Gyy、Gzz、Gxy、Gxz、Gyz分别表示重力梯度数据Vxx、Vyy、Vzz、Vxy、Vxz、Vyz的灵敏度矩阵；xi＝x-ξi，yj＝y-ηj，zk＝z-ζk，μijk＝(-1)i+j+k；γ表示万有引力常数；x、y和z为任意观测点坐标；ξ、η和ζ分别为网格单元在X、Y、Z轴方向上的对顶点坐标，即每层第1个网格单元的空间分布范围；
步骤4：当需要重力梯度数据各分量的灵敏度矩阵参与计算时，获取矩阵元素值的真实索引并换算成压缩后的灵敏度矩阵中的索引；包括：
步骤4.1：根据步骤2可知，调用灵敏度矩阵元素需要知道其所在的行、列数，但由于并行化的是与网格单元总数量有关的计算，利用公式(20)换算每个GPU中分得的网格单元的列索引i_col：
i_col＝i_local+base (20)
式中，i_local表示各GPU中的索引，i_col即第i_col个网格单元的序号；
步骤4.2：计算第i_col个网格单元所在层数以及在该层的行、列数：
i_x＝((i_col-1)％N)％nx+1 (21)
i_y＝((i_col-1)％N)/nx+1 (22)
i_z＝(i_col-1)/N+1 (23)
式中，i_x、i_y分别表示网格单元所处的列数、行数；即分别对应的测点、测线序号；％表示求余数的运算；nx为每条测线上的测点数；i_z表示网格单元所在层数；
步骤4.3：利用灵敏度矩阵元素的行索引计算测点的行、列索引：
obs_x＝(j_row-1)％nx+1 (24)
obs_y＝(j_row-1)/nx+1 (25)
式中，obs_x、obs_y分别为测点的列、行索引，即第obs_y条测线上第obs_x个测点；j_row为灵敏度矩阵元素的行索引；obs_y、obs_x、i_y和i_x即步骤1中的索引i、j、m、n；
步骤4.4：计算当前GPU中压缩后的灵敏度矩阵元素索引：
p_i＝|obs_x-i_x|+1 (26)
p_j＝|obs_y-i_y|+1 (27)
index＝(i_z-first_layer)×N+(p_j-1)×nx+p_i (28)
式中，index为当前GPU中压缩后的灵敏度矩阵元素索引值；通过读取当前GPU中第index个压缩后的灵敏度矩阵元素值，完成对元素值的调用；
步骤5：根据换算后的索引读取压缩后的灵敏度矩阵中的数值，实现灵敏度矩阵的调用。
采用一组双立方体模型检验方法的有效性：设在笛卡尔坐标系下存在两个完全相同的立方体，它们的剩余密度均为0.5g/cm3，在x轴方向上分布范围分别为-700～-300m、+300～+700m，y轴方向上分布范围均为-200～+200m，z轴方向上分布范围均为200～600m；令观测数据均匀分布于z＝0m的平面上，x轴、y轴方向上分布范围均为-980～+980m、间隔均为40m，为了便于测试，观测数据量仅选为50×50；其中，z轴垂直向下为正方向，并设地下空间最大反演深度为1000m，等分为20层，则此时灵敏度矩阵规模为2500×50000；对正演的数据加入5％的高斯白噪声以模拟真实情况。
本组模型试验是在如下配置的工作站上运行测试的：HP Z8 G4，CPU：2*IntelXeon 4110，GPU：4*NVIDIA Quadro P2000，内存：64GB。测试方式如下：首先，根据上述建模方式正演计算理论模型的重力全张量梯度数据；其次，利用该全张量梯度数据中的Vxx、Vyy、Vzz、Vxy、Vxz、Vyz开展三维非线性联合反演，在反演过程中记录或估算灵敏度矩阵占用内存；最后，分析矩阵占用内存情况，评价本方法的应用效果。
本试验中将灵敏度矩阵存储为双精度浮点型数据，约占5.59GB内存，由本方法的发明内容可知，灵敏度矩阵被有效压缩为原大小的1/2500，仅占约2.29MB的内存，而采用了4个GPU开展并行化后，单个GPU中分得的灵敏度矩阵最终仅占不足1MB的显存。这相当于在单个GPU中，原灵敏度矩阵内存被压缩了1万倍，大幅降低了灵敏度矩阵元素的存储数量和内存占用空间，也间接地提升了计算速度。因此，本方法被证明具有实用性与可行性。
Claims (4)
1.一种基于多GPU的重力梯度数据灵敏度矩阵压缩及调用方法，其特征在于，包括：
步骤1：建立重力梯度数据各分量的灵敏度矩阵的等价计算公式；
步骤2：针对多GPU的并行运算，计算各GPU中分得灵敏度矩阵元素的范围、矩阵元素值的真实索引；
步骤2.1：计算每个GPU分得的最后一个网格的索引last_m：
last_m＝base+local_m (7)
式中，base表示排序在当前GPU之前的各GPU中网格数量的和，用于将当前GPU中网格索引转换为真实索引；local_m表示各GPU中分得的网格单元数量；
所述各GPU中分得的网格单元数量local_m：
local_m＝M/num_gpu (12)
式中，M为划分的网格单元总数量；num_gpu表示使用的GPU的数量；/表示求商运算，结果向下取整；
当M不能被num_gpu整除时，则网格单元数量local_m表示为：
local_m＝(M-rem_m)/num_gpu+1 (13)
式中，rem_m表示余数；
步骤2.2：计算当前GPU中分得的网格所在层数的范围：
first_layer＝base/N+1 (8)
last_layer＝(last_m-1)/N+1 (9)
式中，first_layer和last_layer分别表示分得网格所在的第一层和最后一层；N表示每层网格单元的数量；
步骤2.3：计算每个GPU分得的压缩矩阵大小：
size_mat＝(last_layer-first_layer+1)×N (10)
式中，size_mat表示每个GPU分得的压缩矩阵尺寸，即矩阵元素数量；
步骤2.4：计算用于将当前GPU中压缩灵敏度矩阵元素索引转换为真实压缩矩阵索引的参数base_mat：
base_mat＝(first_layer-1)×N (11)
其中，base_mat表示排序在当前GPU之前的各GPU中压缩矩阵元素数量的和；
步骤3：计算压缩后的重力梯度数据各分量的灵敏度矩阵；
步骤4：当需要重力梯度数据各分量的灵敏度矩阵参与计算时，获取矩阵元素值的真实索引并换算成压缩后的灵敏度矩阵中的索引；
步骤5：根据换算后的索引读取压缩后的灵敏度矩阵中的数值，实现灵敏度矩阵的调用。
2.根据权利要求1所述的一种基于多GPU的重力梯度数据灵敏度矩阵压缩及调用方法，其特征在于，所述步骤1中等价计算公式表示为：
Gxx(i,j,m,n)＝Gxx(|i-m|+1,|j-n|+1,1,1) (1)
Gyy(i,j,m,n)＝Gyy(|i-m|+1,|j-n|+1,1,1) (2)
Gzz(i,j,m,n)＝Gzz(|i-m|+1,|j-n|+1,1,1) (3)
式中，Gxx(i,j,m,n)、Gyy(i,j,m,n)、Gzz(i,j,m,n)、Gxy(i,j,m,n)、Gxz(i,j,m,n)、Gyz(i,j,m,n)分别表示等价后的灵敏度矩阵，索引i、j、m、n表示测网中第i行第j列的测点索引与地下某一层中第m行第n列对应的网格索引。
3.根据权利要求1所述的一种基于多GPU的重力梯度数据灵敏度矩阵压缩及调用方法，其特征在于，所述步骤3中重力梯度数据各分量的灵敏度矩阵的计算公式为：
式中，Gxx、Gyy、Gzz、Gxy、Gxz、Gyz分别表示重力梯度数据Vxx、Vyy、Vzz、Vxy、Vxz、Vyz的灵敏度矩阵；xi＝x-ξi，yj＝y-ηj，zk＝z-ζk，μijk＝(-1)i+j+k；γ表示万有引力常数；x、y和z为任意观测点坐标；ξ、η和ζ分别为网格单元在X、Y、Z轴方向上的对顶点坐标。
4.根据权利要求1所述的一种基于多GPU的重力梯度数据灵敏度矩阵压缩及调用方法，其特征在于，所述步骤4包括：
步骤4.1：利用公式(20)换算每个GPU中分得的网格单元的列索引i_col：
i_col＝i_local+base (20)
式中，i_local表示各GPU中的索引；
步骤4.2：计算第i_col个网格单元所在层数以及在该层的行、列数：
i_x＝((i_col-1)％N)％nx+1 (21)
i_y＝((i_col-1)％N)/nx+1 (22)
i_z＝(i_col-1)/N+1 (23)
式中，i_x、i_y分别表示网格单元所处的列数、行数；％表示求余数的运算；nx为每条测线上的测点数；i_z表示网格单元所在层数；
步骤4.3：利用灵敏度矩阵元素的行索引计算测点的行、列索引：
obs_x＝(j_row-1)％nx+1 (24)
obs_y＝(j_row-1)/nx+1 (25)
式中，obs_x、obs_y分别为测点的列、行索引；j_row为灵敏度矩阵元素的行索引；
步骤4.4：计算当前GPU中压缩后的灵敏度矩阵元素索引：
p_i＝|obs_x-i_x|+1 (26)
p_j＝|obs_y-i_y|+1 (27)
index＝(i_z-first_layer)×N+(p_j-1)×nx+p_i (28)
式中，index为当前GPU中压缩后的灵敏度矩阵元素索引值。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110941656.3A CN113591030B (zh) | 2021-08-17 | 2021-08-17 | 基于多gpu的重力梯度数据灵敏度矩阵压缩及调用方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110941656.3A CN113591030B (zh) | 2021-08-17 | 2021-08-17 | 基于多gpu的重力梯度数据灵敏度矩阵压缩及调用方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113591030A CN113591030A (zh) | 2021-11-02 |
| CN113591030B true CN113591030B (zh) | 2024-01-30 |
Family
ID=78258239
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110941656.3A Active CN113591030B (zh) | 2021-08-17 | 2021-08-17 | 基于多gpu的重力梯度数据灵敏度矩阵压缩及调用方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113591030B (zh) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116911146B (zh) * | 2023-09-14 | 2024-01-19 | 中南大学 | 三维重力场全息数值模拟及cpu-gpu加速方法 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107577641A (zh) * | 2017-08-21 | 2018-01-12 | 吉林大学 | 一种基于gpu并行的重力梯度张量数据快速密度反演方法 |
| CN110398782A (zh) * | 2019-07-17 | 2019-11-01 | 广州海洋地质调查局 | 一种重力数据和重力梯度数据联合正则化反演方法 |
| CN112147709A (zh) * | 2020-08-03 | 2020-12-29 | 中国海洋石油集团有限公司 | 一种基于部分光滑约束的重力梯度数据三维反演方法 |
| CN112199859A (zh) * | 2020-10-26 | 2021-01-08 | 东北大学 | 一种重力梯度数据联合反演的方法 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| KR102330829B1 (ko) * | 2017-03-27 | 2021-11-24 | 삼성전자주식회사 | 전자 장치에서 증강현실 기능 제공 방법 및 장치 |
| US11175146B2 (en) * | 2017-05-11 | 2021-11-16 | Anantak Robotics Inc. | Autonomously moving machine and method for operating an autonomously moving machine |
-
2021
- 2021-08-17 CN CN202110941656.3A patent/CN113591030B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107577641A (zh) * | 2017-08-21 | 2018-01-12 | 吉林大学 | 一种基于gpu并行的重力梯度张量数据快速密度反演方法 |
| CN110398782A (zh) * | 2019-07-17 | 2019-11-01 | 广州海洋地质调查局 | 一种重力数据和重力梯度数据联合正则化反演方法 |
| CN112147709A (zh) * | 2020-08-03 | 2020-12-29 | 中国海洋石油集团有限公司 | 一种基于部分光滑约束的重力梯度数据三维反演方法 |
| CN112199859A (zh) * | 2020-10-26 | 2021-01-08 | 东北大学 | 一种重力梯度数据联合反演的方法 |
Non-Patent Citations (8)
| Title |
|---|
| 3D density inversion of gravity gradiometry data with a multilevel hybrid parallel algorithm;Hou Zhen-Long etc;APPLIED GEOPHYSICS;第16卷(第2期);141-152 * |
| Focusing geophysical inversion images;Portniaguine O, Zhdanov M S;Geophysics;第64卷(第2期);874-887 * |
| Multi-GPU parallel algorithm design and analysis for improved inversion of probability tomography with gravity gradiometry data;Zhenlong Hou etc;Applied Geophysics;18-27 * |
| Three-dimensional gravity modelling and focusing inversion using rectangular meshes;Commer, M;Geophysical Prospecting;第59卷(第5期);966-979 * |
| 基于深度加权的重力梯度数据联合相关成像反演;侯振隆 等;东 北 大 学 学 报 ( 自 然 科 学 版 );第41卷(第11期);1628-1632 * |
| 基于灵敏度矩阵压缩的重力梯度数据非线性反演方法;候振隆 等;中国地球科学联合学术年会 2019;1528-1529 * |
| 重力梯度欧拉反褶积及其在文顿盐丘的应用;候振隆 等;石油地球物理勘探;第54卷(第2期);472-480 * |
| 重磁遗传算法三维反演中高速计算及有效存储方法技术;地球物理学报(第2期);252-258 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113591030A (zh) | 2021-11-02 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Komatitsch | Fluid–solid coupling on a cluster of GPU graphics cards for seismic wave propagation | |
| CN110045432B (zh) | 基于3d-glq的球坐标系下重力场正演方法及三维反演方法 | |
| Hou et al. | 3D density inversion of gravity gradiometry data with a multilevel hybrid parallel algorithm | |
| CN113591030B (zh) | 基于多gpu的重力梯度数据灵敏度矩阵压缩及调用方法 | |
| Hou et al. | Full tensor gravity gradiometry data inversion: Performance analysis of parallel computing algorithms | |
| CN111856598A (zh) | 一种磁测数据多层等效源上延拓与下延拓方法 | |
| Cuomo et al. | A GPU parallel implementation of the local principal component analysis overcomplete method for DW image denoising | |
| Hou et al. | Fast density inversion solution for full tensor gravity gradiometry data | |
| Hou et al. | Joint nonlinear inversion of full tensor gravity gradiometry data and its parallel algorithm | |
| CN117725354B (zh) | 一种大数据量重力与重力梯度联合的快速正反演方法及系统 | |
| Hou et al. | 3D inversion of vertical gravity gradient with multiple graphics processing units based on matrix compression | |
| Hou et al. | Multi-GPU parallel algorithm design and analysis for improved inversion of probability tomography with gravity gradiometry data | |
| CN117761789A (zh) | 一种基于视觉自注意力机制的大地电磁二维成像方法 | |
| CN112346139B (zh) | 一种重力数据多层等效源延拓与数据转换方法 | |
| Tu et al. | Comparative investigation of parallel spatial interpolation algorithms for building large-scale digital elevation models | |
| Sheng et al. | The improved residual node density based gravity forward method and its application | |
| Yin et al. | A fast 3D gravity forward algorithm based on circular convolution | |
| Peng et al. | Efficient 3D inversion of potential field data using fast proximal objective function optimization algorithm | |
| CN113238284B (zh) | 一种重磁快速正演方法 | |
| CN112904407B (zh) | 复杂地形与干扰条件下的微动勘探方法 | |
| Hariri et al. | Graph variational autoencoder for detector reconstruction and fast simulation in high-energy physics | |
| CN113204057B (zh) | 一种基于多层级方法重磁快速反演方法 | |
| Lu et al. | A parallel numerical algorithm by combining MPI and OpenMP programming models with applications in gravity field recovery | |
| Protasov | Parallel implementation of 3D seismic beam imaging | |
| CN116794735B (zh) | 航空磁矢量梯度数据等效源多分量联合去噪方法及装置 |
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 |