3维线程格 gpu_基于CUDA的GPU并行优化重力三维反演

重力勘探由于其成本较低、施工方法方便等, 被广泛应用于大尺度的地质异常体勘查、大范围找矿普查、以及小比例尺密度三维地质建模等工作中。目前常用的反演方法有两种, 2.5维联合3维界面反演[和三维物性反演。通过对地下介质进行三维网格剖分, 反演每个网格单元的物性, 从而得到地质异常体的物性分布与形态特征[。

对于复杂结构目标体的三维反演, 由于受到数据计算与存储量的限制, 主要采用的是三维物性反演。即便如此, 重力场计算的三维正演模拟和反演中依旧面临着很大的数据处理量, 使用传统的线性反演, 往往数百万个网格点的计算就会超过个人计算机的内存处理范围。由于计算机内存与CPU计算效率的瓶颈, 很大程度上限制了数据处理解释工作的品质。为了解决这个问题, 近十几年来, 并行计算成为了主流的研究方向。Chen等[采用了基于C语言的CUDA代码实现了重力的三维正演成像; 林巍等[利用GPU遗传算法实现了海洋重力三维反演; Cuma等[实现了基于CPU的OpenMP算法加速与基于GPU的OpenACC算法加速, 提供了在服务器等高性能计算机集群上实现重力反演计算的方法。从硬件结构上来看, 相较于目前仅有数个、数十个核的CPU, 拥有数百甚至数千个核的GPU在处理大数据量的简单运算, 如重力反演大矩阵运算上, 有着天然的优势。

笔者在PGI Fortran平台上采用传统Fortran语言与CUDA语言混编, 将重加权正则化共轭梯度(Re-weighting Regularized Conjugate Gradient, RRCG)三维重力场反演算法进行重组, 使其可以在个人台式计算机与笔记本计算机上使用GPU进行并行计算, 正则化的使用使反演结果稳定可信、具有物理意义; 放弃存储敏感度矩阵A, 从而使得计算效率与数据处理量同时得到提高; 引入了线程多次调用与敏感度范围(footprint)搜索, 实现了存储扩容与反演加速。为了验证程序的正确性, 我们采用相同算法的MATLAB传统代码进行了比对测试。该方法可以有效的在个人计算机上解决上亿网格点的反演计算问题, 为解决大数据量的三维重力反演提供了新的方法思路。

1 正演与反演

1.1 正演算法

虽然对于重力分量体积分的解析解是可以直接求取的, 但为考虑并行计算的效率, 笔者仍采用常用的离散化单点积分:将地下介质划分成均匀的立方体网格, 并假设网格内部密度均匀统一, 经过网格离散化以后, 在重力测量中常用的垂直分量计算公式可以表示为

离散形式。其中:γ表示万有引力常数; Nm表示立方体网格总数; k表示网格点编号; 矢量r'= (x', y', z'), r=(x, y, z)分别表示观测点和源网格中心点的位置; z和z'分别表示观测点和源中心点在垂直方向上的位置; Δ x, Δ y, Δ z分别表示划分立方体网格沿相对坐标轴三个方向的长度; ρ表示立方体网格的密度。<