基于形函数的非均匀傅里叶变换在地震不规则网格建模中的应用

作者与发表信息 本研究的作者为欧阳芳（Fang Ouyang）、赵建国（Jianguo Zhao）、戴世坤（Shikun Dai）、陈龙伟（Longwei Chen）和王尚旭（Shangxu Wang）。他们分别来自中国石油大学（北京）地球物理学院、中南大学地球科学与信息物理学院以及桂林理工大学地球科学学院。该研究论文发表于国际地球物理学期刊 Geophysics 2021年7-8月刊（第86卷，第4期），具体发表日期为2021年6月1日在线发布。

学术背景与研究目的 本研究属于计算地球物理和地震波数值模拟领域，核心关注点是傅里叶变换算法的改进与应用。在众多地震波正演问题中，例如2.5D频率域波场模拟和全波场模拟，通常需要进行一维（1D）或二维（2D）傅里叶变换，以将频率-波数域（wavenumber-domain）的波场谱转换回空间域。传统的快速傅里叶变换（Fast Fourier Transform, FFT）虽计算高效，但其要求数据在空间或波数域必须均匀采样。这在实际地球物理问题中存在明显限制：一方面，复杂介质和波场分布常常要求灵活的、非均匀的采样策略以捕捉关键特征（如近源区的快速振荡）；另一方面，基于FFT的算法在处理非均匀数据时会遇到困难。

为了解决非均匀采样下的傅里叶变换问题，非均匀快速傅里叶变换（Nonuniform Fast Fourier Transform, NUFFT）被发展出来，它在保留FFT速度优势的同时，允许在非均匀节点上进行计算。然而，现有NUFFT方法在多维情况下（如2D）通常要求网格是多个一维非均匀网格的张量积，这限制了其在描述任意复杂形状波场分布时的灵活性和精度。特别是对于具有快速振荡特性的波场（如近源地震波场），常规NUFFT算法难以保证高精度。因此，开发一种既能适应任意不规则采样，又能为复杂波场（尤其是近源区）提供高精度傅里叶变换的算法，对于推进频率域地震建模及其他地球物理正演问题具有重要意义。

本研究的主要目的，正是基于“形函数”（Shape Function）策略，开发一系列适用于一维和二维场景的非均匀傅里叶变换新算法。这些算法旨在结合有限元方法中常用的三角形网格的灵活性，以更精确地描述复杂波场，并通过对指数项振荡特性的不同处理方式，在计算精度和效率之间取得平衡。

详细工作流程 本研究的工作流程可以分为几个关键步骤：算法设计与推导、数值验证与测试、以及性能综合评估与优化。

第一步：算法设计与理论推导 研究团队基于形函数概念，开发了四类核心算法，分为高精度傅里叶变换和基于NUFFT的快速算法两个系列，每个系列又包含线性和二次形函数两种变体： 1. 高精度傅里叶变换系列：包括线性形函数傅里叶变换（Linear Shape-Function-based Fourier Transform, LSF-FT）和二次形函数傅里叶变换（Quadratic Shape-Function-based Fourier Transform, QSF-FT）。这两种算法的核心思想是，在将整个积分域（波数域）划分为小的线元（1D）或三角形单元（2D）后，仅在每个单元内用线性或二次形函数来近似波场谱本身，而将傅里叶变换核中的复指数项 e^(i k y)（1D）或 e^(i (k_x x + k_y y))（2D）完整保留。然后，在每个单元内进行解析积分。这种处理方式的关键优势在于，它通过Filon积分方法，避免了复指数项快速振荡对数值积分精度的影响，从而理论上可以获得很高的精度。在2D情况下，使用三角形单元作为基本单元，使其能够灵活地离散任何复杂形状的波场支持域。 2. 基于NUFFT的快速算法系列：包括线性形函数NUFFT（LSF-NUFFT）和二次形函数NUFFT（QSF-NUFFT）。与高精度系列相反，这类算法将整个被积函数（即波场谱与复指数项的乘积）作为一个整体，在每个单元内用线性或二次函数来近似。然后，通过计算单元积分得到节点权重系数，最终将变换问题转化为一个标准的NUFFT计算格式。由于复指数项的振荡特性没有被单独解析处理，而是被包含在近似的函数整体中，这类算法继承了标准NUFFT的高计算效率，但在处理快速变化的函数（如近源场）时可能精度受限。

研究详细推导了上述四种算法在一维和二维情况下的局部权重系数计算公式。对于1D情况，直接解析求解积分；对于2D情况，则通过雅可比变换将任意三角形映射到一个标准的等腰直角三角形上，再进行解析积分求解。论文附录提供了这些系数计算的详细公式，包括对奇点（如y=0）的特殊处理。

第二步：数值验证与测试 为了评估新算法的性能，研究设计了系统的数值实验，主要应用于两个典型的地震正演问题：2.5D频率域波场模拟（需要1D逆傅里叶变换）和全波场模拟（需要2D逆傅里叶变换）。

1. 研究模型与数据生成：

   • 模型：使用了三个模型进行测试：(a) 均匀全空间模型，用于与解析解（Ben-Menahem和Singh, 1981）进行对比；(b) 三层水平层状模型，用于生成参考解（通过高精度网格的FEM计算）；© 复杂的BP盐丘模型（Billette和Brandsberg-Dahl, 2005），用于验证算法在复杂介质中的适用性。
   • 波场谱计算：空间域的波场谱 ū(x, k_y, z)（1D情况）和 û(k_x, k_y, z)（2D情况）并非通过新算法生成，而是利用作者团队之前开发的基于有限元（FEM）的频率域正演方法（Zhao et al., 2017）独立计算得到。这确保了测试的焦点完全集中在傅里叶变换算法的性能上。
   • 采样策略：针对波数域数据的特点，研究采用了非均匀采样网格。他们观察到波场能量主要集中在一个以截止波数 k_c（约为剪切波波数 ω/v_s）为中心的有限范围内。因此，将波数域划分为两个区域：有效区域（|k| ≤ k_c 或 k_x^2 + k_y^2 ≤ k_c^2），在此区域采用密集采样；扩展区域（k_c < |k| ≤ 10k_c），在此区域采用稀疏采样，但在截止波数 k_c 附近进行局部加密。这种自适应的非均匀网格正是新算法优势得以发挥的场景。

2. 测试流程与对比：

   • 将计算得到的波场谱，分别用新开发的四种算法（LSF-FT, QSF-FT, LSF-NUFFT, QSF-NUFFT）进行逆傅里叶变换，得到空间域的位移场分量（如 u_z）。
   • 在均匀模型中将结果与解析解对比，在层状模型中将结果与高精度参考解对比。
   • 详细记录了计算时间，以评估算法效率。
   • 在BP盐丘模型中，主要展示了QSF-FT和QSF-NUFFT算法在复杂介质中恢复波场的可行性。

第三步：性能综合评估与算法优化（组合算法开发） 基于前述测试结果，研究团队对四种算法的精度和效率进行了全面分析。他们发现，在2D情况下，高精度算法（QSF-FT）在所有接收点（包括近源点）都表现出极高的精度，但计算速度非常慢；而快速算法（QSF-NUFFT）在远离震源的接收点效率极高且精度可接受，但在近源点由于无法有效处理复指数项的剧烈振荡，误差显著增大。

为了同时获得高精度和高效率，他们提出了一个创新的组合算法（Combined NUFFT）：在近源区域（|z - z_s| 小于一个阈值）使用高精度的QSF-FT算法，而在其他区域则使用高效的LSF-NUFFT或QSF-NUFFT算法。研究还通过分析均匀介质中波场谱的解析表达式，推导出了一个用于判断“近源区”的经验性阈值公式 |z - z_s|_{th} = b / (k_c √(a^2 - 1))，其中 a 和 b 为根据数值实验确定的参数。此外，针对近源区波场变化相对简单的特点，组合算法中对QSF-FT部分采用了更稀疏的波数采样网格，进一步提升了整体计算效率。

主要结果 本研究通过系统的数值实验，得到了以下关键结果：

1. 1D算法结果：无论是对于全空间模型还是三层模型，四种1D新算法（LSF-FT, QSF-FT, LSF-NUFFT, QSF-NUFFT）计算得到的波场与精确解/参考解都吻合得非常好，误差极小。二次形函数算法（QSF-FT和QSF-NUFFT）的精度普遍略高于对应的线性算法。在计算效率上，基于NUFFT的算法（LSF-NUFFT和QSF-NUFFT）显著快于高精度算法（LSF-FT和QSF-FT）。

2. 2D算法结果：

   • 精度对比：二次形函数算法（QSF-FT和QSF-NUFFT）再次展现出比线性算法更高的精度。LSF-FT和QSF-FT在所有接收点（包括近源点）都保持了高精度。然而，LSF-NUFFT和QSF-NUFFT在近源点（在测试中，即靠近震源平面的位置）出现了明显的误差峰值，这与理论预期一致，因为它们没有解析处理近源区复指数项的剧烈振荡。
   • 效率对比：计算时间数据（论文表2）清晰地表明，2D LSF-NUFFT和QSF-NUFFT的计算速度比LSF-FT和QSF-FT快一到两个数量级。例如，对于151x151个空间点、使用23,721个波数采样点的情况，LSF-FT需要约315分钟，而LSF-NUFFT仅需约1.7分钟。
   • 组合算法结果：提出的组合算法（CL-NUFFT/CQ-NUFFT）成功地综合了各类算法的优点。如图14所示，组合算法在近源点的精度显著优于单纯的NUFFT算法，与高精度算法相当；而在远源点，它则保持了NUFFT算法的高效率。同时，由于对近源区QSF-FT计算采用了优化的稀疏网格，组合算法的总计算时间远低于在整个区域使用QSF-FT的时间（论文表3），实现了精度与效率的良好平衡。

3. 复杂模型验证：在BP盐丘模型上的应用表明，QSF-FT和QSF-NUFFT算法能够有效地从波场谱中恢复出空间域的波场，其结果与高精度参考解吻合良好，证明了新算法在复杂非均匀介质中同样具有适用性。

这些结果逻辑连贯：首先验证了算法基本框架的有效性（1D和2D基础测试），然后揭示了不同算法在精度和效率上的固有优缺点（2D详细对比），最后基于这些认识，通过巧妙的组合策略设计出性能更优的新方案（组合算法）。每个步骤的结果都为下一步的分析和优化提供了直接依据和数据支持。

结论与研究价值 本研究成功开发了一套基于形函数的非均匀傅里叶变换新算法，包括LSF-FT、QSF-FT、LSF-NUFFT、QSF-NUFFT及其组合版本CL-NUFFT/CQ-NUFFT。主要结论是：二次形函数算法通常比线性算法更精确；高精度傅里叶变换算法（LSF-FT/QSF-FT）精度高但计算慢；基于NUFFT的算法（LSF-NUFFT/QSF-NUFFT）效率高但在近源点精度不足；而提出的组合算法通过在不同区域智能选择算法，能够为所研究的地震正演问题同时提供令人满意的高精度和高效率。

本研究的科学价值在于，首次将有限元分析中的形函数和三角形网格离散思想系统性地引入到非均匀傅里叶变换中，为解决复杂波场分布的变换问题提供了一个灵活且高精度的理论框架。其应用价值显著：这些算法可以直接用于提高频率域地震波正演模拟的精度和灵活性，特别是在处理近源波场、复杂构造模型以及需要非均匀采样的场景时。尽管测试聚焦于地震模型，但作者指出，这些改进的NUFFT算法在重磁位场正演模拟等其他地球物理问题中同样具有潜在的应用前景。

研究亮点 1. 方法创新：核心亮点是将“形函数”策略与傅里叶变换结合，并创造性地分为“保留指数项积分”（高精度系列）和“整体函数近似”（快速系列）两种技术路径，清晰地揭示了精度与效率取舍的内在机理。 2. 灵活性突破：在2D算法中引入三角形单元作为基本积分单元，这是对传统张量积网格限制的重要突破，使得算法能够与任何不规则采样策略兼容，极大地增强了对复杂计算域的适应能力。 3. 实用性解决方案：提出的“组合算法”是一个极具实用价值的工程解决方案。它不追求单一的“最优”算法，而是基于对问题物理特性（近源/远源场行为不同）的深刻理解，通过智能的区域划分和算法调度，巧妙地达到了全局最优的性能平衡。 4. 系统性的验证：研究通过从简单到复杂（全空间→层状→盐丘）、从1D到2D、从精度到效率的全面对比测试，严谨地评估了所有算法的性能，使结论非常可靠。

其他有价值内容 论文的引言部分对FFT、Gauss-FFT以及各类NUFFT技术的发展脉络进行了清晰的梳理，为读者提供了良好的学术背景。附录部分详细给出了权重系数的推导过程，增强了研究的可重复性和理论深度。此外，作者在文中明确说明了所使用的NUFFT代码来源（Greengard, 2009），并提供了数据可用性声明，符合现代学术研究的规范。