本文介绍了一项题为《任意起伏地形下重力异常三维正演及并行计算》的原创性研究工作。该项研究由戴世坤、朱德祥、张莹、李昆、陈轻蕊、凌嘉宣、田红军共同完成，研究者主要来自中南大学有色金属成矿预测与地质环境监测教育部重点实验室、有色资源与地质灾害探查湖南省重点实验室以及地球科学与信息物理学院，其中李昆来自西南石油大学地球科学与技术学院。该研究成果于2024年2月发表在《地球物理学报》（Chinese Journal of Geophysics）第67卷第2期。
一、 研究的学术背景
本研究属于地球物理勘探，特别是重力勘探领域的数值模拟方法研究。重力勘探广泛应用于油气勘探、金属矿勘查、地球内部构造研究等领域。随着勘探精度要求的提高和对复杂地质条件（如剧烈起伏地形）研究的深入，发展高效、高精度的三维重力异常场数值模拟方法，对于实现大规模复杂地质条件下的重力异常精细化反演成像与综合解释具有至关重要的作用。
数值模拟方法主要分为空间域和频率域两大类。空间域方法（如长方体、多面体解析解、数值积分等）对于形状规则、密度均匀的异常体有优势，但对于大规模复杂地形和物性分布的模拟，计算量巨大。频率域方法（利用傅里叶变换）因其表达式简洁、计算高效而展现出优势，但其精度和适用范围受到傅里叶变换算法的制约。传统的快速傅里叶变换（Fast Fourier Transform, FFT）算法需要均匀采样，且存在截断效应（Truncation Effect），即使通过扩边或高斯FFT（Gauss-FFT）等方法改进，仍存在采样不灵活、在某些情况下需选取大量波数才能保证精度的问题，影响了计算效率。
此外，无论是空间域还是频率域方法，在面对大规模三维模型（节点数可达千万级别）时，串行计算耗时过长。高性能计算（尤其是CPU-GPU异构并行）为解决此问题提供了途径。然而，之前基于微分法的空间-波数域重力异常正演算法虽具有良好的并行潜力，但并未实现完整的并行加速方案。
因此，本研究旨在解决两个核心问题：一是发展一种采样更灵活、精度更高、截断效应更小的新型傅里叶变换算法以改进空间-波数域正演的核心环节；二是充分利用算法的并行性，设计并实现一套高效的CPU-GPU异构并行计算方案，从而大幅提升任意起伏地形下大规模三维重力异常正演的计算效率和适用范围。
二、 详细工作流程
本研究的工作流程可以概括为三个核心部分：新算法的理论构建、算法正确性与性能验证、并行计算方案实现与测试。具体步骤如下：
第一， 基于“任意傅里叶变换”的空间-波数域正演算法理论构建。 本研究的核心创新是提出并应用了一种名为“任意傅里叶变换”（Arbitrary Fourier Transform）的新算法。该算法是空间-波数域重力正演的关键步骤，用于实现空间域与波数域之间的高效、高精度转换。 1. 算法原理：首先，将二维傅里叶变换分解为两个连续的一维傅里叶变换。然后，对每个一维积分进行离散化处理，将整个积分区间划分为多个单元。 2. 核心创新点：在每个离散单元内，不再像传统方法那样依赖固定的采样点，而是采用二次插值形函数来拟合原函数。这使得算法可以对空间域和波数域进行灵活的非均匀剖分与采样。 3. 解析求解：基于二次插值形函数，推导出每个单元积分的解析表达式。这样，整个傅里叶变换就转化为对所有单元解析结果求和的过程。由于单元积分系数在剖分确定后可以预先计算并存储，且不同波数下的计算相互独立，该算法天然具备高度并行性。 4. 采样规则：研究者总结了波数选取的准则：在频谱能量大且变化剧烈的区间（通常是低频部分）进行加密采样，在频谱能量小且变化平缓的区间（高频部分）进行稀疏采样。这确保了在保证精度的前提下，尽可能减少所需计算的波数数量，从而提高效率。波数的最大取值范围由空间域的最小网格间距决定（满足采样定理），基频由计算区域的总范围决定。
第二， 算法验证与性能对比研究。 为了验证新算法的正确性、精度和效率，研究者设计并进行了多组数值实验。 1. 研究模型与对象： * 棱柱体模型（验证正确性）：设计了一个密度均匀（2000 kg·m⁻³）的规则长方体模型，置于一定范围的模拟区域中心。模型参数明确，可以获得重力异常场的解析解。计算区域被离散为251(x)×251(y)×151(z)个节点。 * 连续密度模型（对比性能）：设计了一个密度在水平方向连续变化的复杂模型。由于没有解析解，研究采用公认精度高的高斯FFT（n_g=4）算法的计算结果作为“参考解”，用于对比评估本文新算法的精度和效率。 * 实际地形数据模型（验证实用性）：使用真实地形数据进行大规模三维重力异常场数值模拟，模型规模达到721×721×421个节点，以测试算法处理实际复杂问题的能力。 2. 实验方法与处理过程： * 正确性验证：使用新算法计算棱柱体模型在地表和异常体内部两个平面的重力异常三个分量（g_x, g_y, g_z）。将数值解与解析解进行对比，绘制等值线图，并计算相对均方根误差（Relative Root Mean Square Error, RRMS）。 * 精度与效率对比：在连续密度模型上，分别应用基于任意傅里叶变换的算法和基于高斯FFT的算法进行正演计算。比较两种算法在不同波数选择下的计算结果与参考解的误差（RRMS），同时记录各自的耗时，从而在相同精度水平下对比计算效率。 * 并行效果测试：设计不同规模（节点数从十万到千万级）的模型，分别使用CPU串行版本、CPU并行（OpenMP）版本和CPU-GPU异构并行版本的算法进行计算，记录并对比计算时间，计算并行加速比（串行时间/并行时间），以量化并行方案的性能提升。
第三， CPU-GPU异构并行计算方案的设计与实现。 本研究详细解剖了正演算法的计算流程，识别出两个可并行的核心部分：求解五对角线性方程组和进行任意傅里叶变换。 1. 并行策略设计： * 求解方程组（CPU并行）：在空间-波数域方法中，每个波数对应一个独立的一维常微分方程边值问题，通过有限元法离散后形成一个五对角线性方程组。这些方程组彼此独立，可以并行求解。由于追赶法求解过程存在循环和逻辑控制，适合在多核CPU上使用OpenMP指令进行线程级并行。 * 任意傅里叶变换（GPU并行）：任意傅里叶变换算法包含大量规则且相互独立的乘加运算。研究者将计算任务映射到NVIDIA CUDA平台上，由GPU上成千上万个线程并行执行。具体而言，将不同波数、不同y坐标或不同深度z的计算任务分配给不同的GPU线程块（Block）和线程（Thread），通过设计一维、二维或三维的线程层次结构，高效地组织大规模并行计算。 2. 工作流程整合：整个正演程序由CPU主线程控制。主程序首先在CPU上分配任务，将不同波数对应的方程组求解任务分发给多个CPU线程并行计算（OpenMP并行）。同时，将正、反傅里叶变换的巨大计算量提交给GPU，启动成千上万个CUDA线程并行执行。最终，CPU整合所有波数的计算结果，完成正演过程。
三、 主要研究结果
第一， 算法正确性得到验证。 对棱柱体模型的计算结果表明，无论是地表（z=0 m）还是异常体内部（z=500 m），新算法得到的重力异常场数值解与解析解的等值线图几乎完全重合。定量分析显示，所有三个重力分量在内外两个平面上的相对均方根误差（RRMS）均远小于1%（在0.48%至0.65%之间）。这充分证明了基于任意傅里叶变换的空间-波数域正演算法的数学正确性和高精度。
第二， 新算法在效率与精度上展现优势。 在连续密度模型的对比实验中，与高斯FFT算法相比，在达到相近数值精度（RRMS相当）的前提下，本文提出的任意傅里叶变换算法所需的波数更少。这意味着在每次正演计算中，需要求解的方程组数量和傅里叶变换的计算量都相应减少，从而显著提高了计算效率。这验证了新算法“采样灵活、积分精度高、计算速度快”的理论优势。
第三， CPU-GPU并行方案取得显著加速效果。 并行性能测试结果明确显示，相比传统的CPU串行算法，本研究实现的CPU-GPU异构并行算法将计算效率提升了40倍以上。对于节点数达到千万级别的超大规模模型，正演计算时间仅需数秒。文中具体提到，对一个包含721×721×421个节点（约2.19亿个网格单元）的实际地形模型进行正演，耗时仅17秒。这一结果充分证明了并行方案对于处理大规模实际问题的巨大价值。
第四， 实际应用验证了方法的实用性与高效性。 将新算法应用于实际地形数据的重力异常场数值模拟，成功完成了大规模、复杂地形条件下的快速计算。这不仅是一个算例，更直接证明了该方法能够高效处理实际地球物理勘探中的数据，具备重要的工程应用价值。
四、 结论
本研究成功发展了一套适用于任意起伏地形的高效、高精度三维重力异常快速正演及并行计算方案。其科学价值和应用意义主要体现在： 1. 方法学创新：提出的“任意傅里叶变换”算法，突破了传统FFT必须均匀采样的限制，实现了根据频谱特征灵活、自适应选取波数的能力，在保证精度的同时减少了必要的计算量。 2. 计算性能突破：通过精细设计的CPU-GPU异构并行架构，将算法中计算密集且可并行的两部分分别交由CPU和GPU高效处理，实现了对超大规模模型（千万级节点）的秒级正演，极大地提升了计算能力。 3. 实践指导意义：该研究为大规模复杂地质条件下的重力勘探数据精细解释与三维反演提供了强有力的工具。高效的正演是高效反演的基础，本方法有助于推动重力勘探在复杂地区的深度应用。
五、 研究亮点
六、 其他有价值的内容
文中还简要回顾了重力数值模拟方法的发展历程，梳理了空间域方法（如解析解、体积分转面积分、线积分、多极子算法、有限元/有限差分法）和频率域方法（如各种改进的FFT算法）的主要进展与优缺点，为本研究工作的定位和创新点提供了清晰的学术背景。此外，研究对GPU并行中的线程层次结构（一维、二维、三维网格与块组织）进行了简要说明，体现了并行实现的具体技术细节。