近日，一项关于分层介质中三维点源声波场高效数值模拟的研究取得了重要进展。这篇题为《Algorithm for Acoustic Wavefield in Space-Wavenumber Domain of Vertically Heterogeneous Media Using NUFFT》的研究论文，由来自中南大学资源与安全工程学院/地球科学与信息物理学院的张颖与戴思坤（通讯作者）共同完成，并于2025年2月9日在线发表于开源学术期刊《Mathematics》（卷13，期4，文号571）上。该研究针对地球物理勘探与地震学中的核心计算难题，提出了一种结合非均匀快速傅里叶变换的高效求解算法，为三维声波方程的精确、快速求解提供了创新性方案。
本研究隶属于计算地球物理学和科学计算领域，核心目标是解决三维分层介质中点源声波方程的数值求解问题。地震波在分层介质中的传播特性研究，对于地震学、石油勘探、地下成像等领域至关重要。传统解析方法虽精度高，但仅适用于简单模型，对复杂的分层结构束手无策。因此，基于波动方程的数值模拟成为主流工具。
然而，在三维声波场的数值模拟中，效率与精度往往难以兼顾。传统的数值方法（如有限差分法、有限元法等）在面对大规模、复杂模型时，面临着计算成本高昂、内存需求巨大的挑战。尤其是在需要考虑诸如空气层这类低波速、需要密集网格剖分的场景时，计算负担更为沉重。为突破这一瓶颈，研究者旨在开发一种既能保持高精度，又能大幅降低内存占用和计算时间的新型算法。本研究正是在此背景下展开，其核心目标是将复杂的三维问题降维简化，并引入灵活的频谱变换技术，从而在计算效率和数值精度之间找到最佳平衡点。
本研究设计并实现了一套完整、创新的数值求解流程，主要包含以下五个关键步骤：
问题域的转换：从三维空间到空间-波数域 研究的第一步是对三维频域声波方程（控制方程）进行数学变换，实现问题的降维。研究者首先对三维空间域中的声波方程在水平方向（x，y方向）进行二维傅里叶变换。由于假定介质波速仅在垂直方向（z方向）变化（即垂向非均匀介质），这一变换成功地将原本的三维偏微分方程，转化为一个仅关于深度z的一维常微分方程。该方程的解被称为“空间-波数域”的波场。这一步骤是本研究的核心思想之一，因为它将庞大的三维问题瞬间简化为一系列独立的一维问题，每个一维问题对应于一个特定的水平波数对（kx， ky），从而在理论上极大地降低了对计算内存的需求。
边界值问题的建立与求解：一维有限元法 转换后得到的一维方程，连同其上下边界条件（分别对应只有上行波和下行波的辐射边界条件），构成了一个完整的边值问题。研究者采用伽辽金有限元法来求解这个一维问题。通过离散化垂直方向，并利用二次插值形函数，最终将边值问题转化为一个大型线性方程组。由于形函数是二次的，每个单元涉及三个节点，最终组装得到的全局刚度矩阵是一个五对角矩阵。为了高效求解这个方程组，研究者采用了经典的托马斯算法。该算法专门用于求解三对角或块三对角系统，对于这里的五对角结构同样高效。通过求解这个系统，研究者获得了在给定频率和水平波数下的、整个垂直剖面上的波场分布（即在“空间-波数域”中的解）。
波数选择的策略与优化 为了从空间-波数域最终恢复出三维空间波场，需要在水平波数域进行足够精确的采样。波数选择的合理与否直接关系到最终结果的精度和计算量。本研究没有盲目地采用均匀密集采样，而是基于对点源声波场物理特性的深刻理解，提出了一个高效的采样策略。研究者分析了全空间点源格林函数在波数域的表达式，发现波场能量主要集中在临界波数kc（kc = ω/c， ω为角频率，c为波速）附近，当波数绝对值大于2kc后，波场谱值已衰减至近乎为零。因此，研究者将整个模型的波数采样范围统一设定为[-2kcmax, 2kcmax]，其中kcmax取整个模型中各层的最小波速所对应的临界波数。这样，既能捕获波场的核心频谱信息，又避免了在无效波数区域进行不必要的计算，为后续的高效逆变换奠定了基础。
核心创新：应用非均匀快速傅里叶变换 获得空间-波数域的解后，最后一步是通过水平方向的二维逆傅里叶变换，将其重构回三维空间域波场。这是本研究最具创新性的环节。传统的快速傅里叶变换存在一个固有局限：一旦空间域采样确定，波数域的采样（间隔和范围）也就随之固定，缺乏灵活性。为了在保证精度的前提下进一步提高计算效率，并允许非均匀采样（这在实际应用中非常有用），研究者引入了第三类非均匀快速傅里叶变换。NUFFT是一种能够处理非均匀网格点数据与变换结果的广义傅里叶变换算法。 研究采用的方法基于高斯核卷积与去卷积的“网格化”思想。其具体工作流程是：首先，利用高斯函数将空间-波数域中可能非均匀的波数采样点，卷积平滑并插值到一个均匀的中间网格上；然后，对这个均匀网格上的数据进行标准的快速傅里叶变换，得到空间域均匀网格上的中间结果；最后，再通过高斯核的去卷积和插值，将结果映射回最终期望的空间域网格点上（该网格同样可以是任意非均匀的）。这一方法的优势在于，它允许研究者独立、灵活地选择空间域和波数域的采样策略，无需受制于传统FFT的互锁关系，从而可以用更少的采样点达到所需的精度，显著节省了计算资源。
流程整合与算法验证 以上步骤被整合为一个完整的算法流程：输入模型参数和源信息后，算法首先为每个频率和每个选定的水平波数对，求解一维有限元问题，得到空间-波数域波场；然后，利用NUFFT对所有波数进行二维逆变换，合成最终的三维空间波场。为了验证算法的正确性、精度、适应性和效率，研究者设计并执行了系统的数值实验。
实验在一台配置为Intel i7-11800H处理器、64GB内存的计算机上进行，主要包含以下四个部分：
正确性验证：与全空间解析解对比 研究者首先设计了一个简单的全空间均匀介质模型，将本文算法计算的波场与已知的解析解进行对比。结果显示，无论是波场的实部还是虚部，数值解与解析解的波形都高度吻合。定量分析采用相对均方根误差作为精度指标，其实部和虚部的RRMS误差均小于1%。这一结果充分证明了算法从理论推导到程序实现的正确性，为后续更复杂模型的测试提供了信心。
半空间模型测试：验证含空气层计算的准确性 为检验算法处理波速剧烈变化和需要密集网格剖分场景的能力，研究者设计了一个半空间模型：上层为250米厚的空气层（波速340 m/s），下层为250米厚的固体介质（波速3000 m/s）。点源位于空气层中。由于空气波速极低，波长很短，通常需要非常密集的垂直网格才能准确模拟，这对算法的内存和效率是巨大考验。研究者在不同频率（5 Hz， 10 Hz， 20 Hz）下计算了波场。结果显示，算法能够清晰、准确地刻画出波场在半空间中的传播、在界面处的反射等物理现象，并且在不同频率下，波长的变化规律（频率越高，波长越短，波场振荡越剧烈）也得到正确体现。这表明，该算法能够有效处理包含空气层这类“苛刻”条件的模型，且保持了高精度。
三层介质模型测试：验证对复杂分层结构的适应性 为进一步验证算法对典型分层地质结构的适应性，研究者构建了一个三层介质模型，各层波速从浅到深递增（1000， 2000， 3000 m/s）。计算结果表明，波场在频域表现出明显的分层特征。随着深度增加、波速增大，相同频率下的波长变长，分辨率降低，这与实际地震波传播规律完全一致。该实验成功证明了算法在计算分层介质声波场分布方面的准确性和鲁棒性，为将其应用于更实际的地质模型奠定了基础。
敏感性分析与效率评估 这是本研究的亮点之一，通过系统的参数测试，量化了算法性能。
本研究成功提出并验证了一种用于求解垂向非均匀介质中三维点源声波方程的高效数值算法。其核心贡献在于：通过水平方向的傅里叶变换实现降维，将三维问题转化为一系列一维问题；利用非均匀快速傅里叶变换实现灵活、高效的谱域与空间域转换。这两者的结合，使得算法在保证高数值精度的同时，实现了极低的内存消耗和可观的计算速度。
该研究的科学价值在于，为求解一类重要的偏微分方程（具有垂向变化系数的亥姆霍兹方程）提供了全新的思路和方法框架。其应用价值尤为显著：该算法为散射波地震数值模拟中的背景场计算提供了高效、精确的工具，这是全波形反演等关键技术的前置步骤。同时，算法本身的高效性使其能够在普通个人计算机上处理大规模模型，降低了高性能计算的门槛，具有强大的工程应用潜力。
论文还对波数选择的物理依据进行了深入探讨，基于点源格林函数的性质确定了高效的采样范围（-2kc ~ 2kc），这并非简单的经验设定，而是有坚实的物理和数学基础，体现了研究的严谨性。此外，研究在引言部分对现有数值方法（如有限差分、边界元、谱元法等）进行了简要回顾，指出了它们在处理大规模分层介质问题时面临的共同挑战，从而自然而然地引出了本工作的必要性和创新点。