本文介绍一篇发表于《地球物理学报》2024年2月第67卷第2期的原创性研究论文。该论文由戴世坤、朱德祥、张莹等作者共同完成，作者单位包括中南大学有色金属成矿预测与地质环境监测教育部重点实验室、中南大学有色资源与地质灾害探查湖南省重点实验室、中南大学地球科学与信息物理学院，以及西南石油大学地球科学与技术学院。论文题目为《任意起伏地形下重力异常三维正演及并行计算》。研究的核心目标是提升空间-波数域三维重力异常正演算法在复杂地形下的适用性与计算效率，通过引入一种新型的“任意傅里叶变换算法”，并基于CPU-GPU异构平台实现并行加速，最终为大规模复杂条件下的重力勘探精细化反演解释提供高效、高精度的数值模拟工具。

首先，在学术背景方面，该研究属于地球物理勘探，特别是重力勘探领域的数值模拟方法学创新。重力勘探是资源勘查（如油气、金属矿）、地球内部构造研究及环境地球物理调查的重要手段。其核心环节之一是利用数值模拟方法高效、准确地计算地下密度异常体在地表产生的重力异常场。现有的数值模拟方法主要分为空间域和频率域两大类。空间域方法（如解析解、高斯积分、多极子算法等）对简单几何体高效，但对复杂地质体和大规模计算适应性不足。频率域方法（通过傅里叶变换将偏微分方程转化为常微分方程求解）表达式简洁、计算效率高，但其精度和灵活性严重受限于传统傅里叶变换算法（如标准FFT、高斯FFT）。这些传统方法存在采样点必须均匀分布、波数（wavenumber）选取灵活性差、截断效应难以完全消除等固有局限，限制了算法对任意复杂地形和变密度模型的适用性。此外，尽管空间-波数域算法本身具备高度并行性，但在实际应用中，尤其是在处理大规模模型时，其计算效率仍有巨大的提升空间。因此，本研究旨在解决两个关键问题：一是突破传统傅里叶变换在采样和精度上的限制；二是充分利用现代异构计算硬件的优势，将算法计算效率最大化。其研究目标是开发一种采样灵活、精度高、速度快的三维重力异常正演算法，并实现高效的并行计算方案。

其次，研究的工作流程细致而复杂，主要包含两大核心模块的开发与验证：任意傅里叶变换算法的实现 与 CPU-GPU异构并行加速方案的构建，并辅以一系列的数值实验进行验证。

第一个主要流程是理论方法与算法构建。 研究的基础是空间-波数域混合域方法。该方法从描述重力位与密度关系的泊松方程出发，对水平方向进行二维傅里叶变换，将三维偏微分方程转化为一系列不同波数下的一维常微分方程。这些方程结合上下边界条件，构成了边值问题，随后可采用一维有限元法（基于二次插值形函数）进行求解，最后通过傅里叶反变换得到空间域的重力异常场。本研究的关键创新在于取代了传统FFT，引入了“任意傅里叶变换算法”。该算法的核心思想是将二维傅里叶变换分解为两个一维变换，并对每个一维变换积分进行灵活剖分离散。具体步骤为：1）将积分区间离散为多个单元，单元划分可以是均匀或非均匀的；2）在每个离散单元内，使用二次插值形函数对原函数进行高精度拟合；3）基于此拟合，推导出每个单元积分的解析表达式（论文中给出了具体的系数公式，如式13和式15）。这种方法的优势在于：采样点可以任意分布，可根据频谱能量集中的区域进行加密采样，在能量变化平缓处稀疏采样，从而用更少的波数达到所需的精度；同时，由于得到了单元积分的解析解，其积分精度远高于数值积分，并能有效减小傅里叶变换的截断效应。研究还总结了该算法的波数采样规则：波数的最大范围由空间域的最小网格间距决定（基于采样定理），而具体的波数值则根据计算得到的频谱能量分布进行自适应选取，即在能量大、变化剧烈的区域加密采样，反之则稀疏采样。

第二个主要流程是并行计算架构的设计与实现。 研究者深入分析了算法中计算密集且可并行化的部分。主要有两部分：一是对不同波数下常微分方程（最终化为五对角线性方程组）的求解；二是对大量波数进行正/反二维任意傅里叶变换计算。研究采用了CPU-GPU协同的异构并行策略：1）对于求解五对角方程组，其计算过程涉及复杂的循环和迭代，指令控制逻辑较强，因此选用CPU利用OpenMP进行多线程并行计算。主线程分配任务给多个子线程，每个线程独立求解一个波数对应的方程组。2）对于任意傅里叶变换，其计算本质是大量相互独立的乘法和加法运算，非常适合GPU的大规模并行计算。研究者设计了细致的GPU核函数，将计算任务组织成网格（Grid）和线程块（Block）的层次结构。例如，在计算沿x方向的一维积分时，可以将需要计算的波数（kx）和y方向节点数（ny）的二维组合任务映射到二维线程网格中，每个GPU线程负责一个特定波数在某一个y坐标点的积分运算，从而实现了极高的并行度。这种分工充分发挥了CPU擅长复杂逻辑控制和GPU擅长大规模简单数据并行计算的优势。

第三个主要流程是数值算例验证。 为了全面验证算法的正确性、精度和效率，研究者设计了三个模型进行测试。实验使用的硬件平台为配备18核Intel Core i9 CPU和NVIDIA Titan RTX GPU的计算机。程序正确性验证：设计了一个简单棱柱体常密度模型，将新算法得到的重力异常三个分量（gx, gy, gz）的数值解与已知的解析解进行对比。通过计算相对均方根误差（RRMS），无论是在异常体外部（地表）还是内部，误差均低于0.65%，证明了算法的正确性。计算精度与效率对比：设计了一个连续变化的变密度模型，以公认高精度的“高斯快速傅里叶变换”算法的结果作为参考基准。对比了不同波数选取策略下，任意傅里叶变换算法与高斯FFT算法的精度和计算时间。结果显示，在达到相近数值精度（RRMS < 1%）的前提下，本文提出的任意傅里叶变换算法所需选取的波数数量更少，从而计算效率更高。并行加速效果测试：通过计算不同规模模型（从百万到千万级节点）的重力异常，对比了CPU串行算法与CPU-GPU并行算法的耗时。结果表明，并行加速效果显著，对于千万节点级别的模型，正演计算仅需数秒即可完成，相比CPU串行计算，加速比（串行时间/并行时间）达到了40倍以上。实际应用演示：最后，研究者将新算法应用于一个基于真实地形数据的复杂模型进行三维重力异常场数值模拟。模型规模巨大（721×721×421个节点），但利用新算法的并行加速方案，整个计算过程仅耗时17秒，充分证明了该方法处理大规模实际问题的“高效性”与“实用性”。

接下来，研究的主要结果逻辑清晰，环环相扣。1）算法正确性得到验证：棱柱体模型数值解与解析解的高度吻合（RRMS < 1%），为后续所有基于该算法的比较和应用奠定了可信的基础。2）算法精度与效率优势得以量化：在变密度模型中，与高斯FFT这一高精度标杆的对比结果表明，新算法能够以更少的计算量（更少的波数）达到同等精度，这直接体现了其在“灵活性”和“效率”上的理论优势。3）并行加速方案效果显著：千万级节点模型“数秒内”完成计算的测试结果，是整套技术路线（任意变换+异构并行）价值的最终体现，它将理论上的高效潜力转化为了实际应用中可感知的性能飞跃。4）实用性得到最终证明：实际地形数据模拟的成功，是前序所有理论推导、算法实现和性能测试的最终落地，表明该方法不仅是一个“学术原型”，更是能够解决“大规模复杂条件下”实际地球物理问题的有力工具。这些结果层层递进，共同支撑了论文的最终结论。

本研究的结论明确而有力。论文成功开发了一种基于“任意傅里叶变换”和“CPU-GPU异构并行”的空间-波数域三维重力异常快速正演算法。该算法有效克服了传统频率域方法采样不灵活、截断效应影响精度的缺点，实现了在复杂起伏地形和变密度条件下的高精度、高效率数值模拟。其科学价值在于，为频率域重力正演提供了一种更为灵活和精确的傅里叶变换数值实现途径，丰富了地球物理数值模拟的方法体系。其应用价值则更为突出：算法卓越的计算性能（处理千万节点模型仅需数秒）使得对大规模、复杂地质条件下的重力异常场进行精细化三维数值模拟与快速反演成像成为可能，这对于提升重力勘探在深部资源勘查、精细构造解释等领域的作用具有重要意义。

本研究的亮点突出体现在以下几个方面：第一，方法学的创新性：提出的“任意傅里叶变换算法”是其核心理论创新，它通过离散单元内二次插值拟合和解析积分，巧妙地实现了非均匀采样和高精度计算，突破了传统FFT的框架限制。第二，技术集成的先进性：并非仅仅停留在算法理论层面，而是深度融合了高性能计算的最新实践，设计了针对算法特点的CPU-GPU异构并行方案，将理论算法的潜力通过并行编程技术充分释放，体现了显著的工程实现价值。第三，验证体系的完整性：研究从正确性（与解析解对比）、精度（与高精度算法对比）、效率（串行与并行对比）到实用性（实际地形应用）构建了完整的验证链条，结论坚实可信。第四，明确的工程导向：研究始终以解决实际计算瓶颈为导向，无论是追求更少的波数、更高的精度，还是最终实现秒级的千万节点计算，都直指大规模工业级应用的需求。这些亮点共同使得该项研究不仅是一篇优秀的理论方法论文，也是一项具有重要实用价值的技术成果。