采用混合DST-FD技术对引力势及其导数进行的大规模三维数值模拟

《Journal of Applied Geophysics》：Large-scale three-dimensional numerical simulation of gravitational potential and its derivatives using a hybrid DST-FD technique

　　胡一遥|童晓忠|苏杰宇|魏康贵•该方法可计算三维大尺度重力势、重力场矢量以及重力梯度张量。•所提出的加速方法有效避免了传统有限差分法中所需的大规模矩阵求逆或迭代求解过程。•该算法在保证高精度的同时显著减少了计算时间。•该方法已被应用于三维SEG/EAGE盐模型重力异常的计算中。

引言
重力勘探作为一种基础地球物理技术，其原理在于测量地球重力场的微小变化。通过分析重力异常信号，可以重建地下密度分布，这一技术广泛应用于油气资源评估（Kadirov, 2000; Choi et al., 2017; Chen et al., 2021; Xu et al., 2024）、深部矿产勘探（Mosher and Farquharson, 2013; Kamm et al., 2015; Upadhyay and Sharma, 2020; Vatankhah et al., 2023）、深部基底结构识别（Tschirhart et al., 2016; Sampietro and Capponi, 2022; Kobayashi et al., 2025）以及环境岩土工程调查（Pedersen and Bastani, 2016; Zhang et al., 2026）等领域。近年来，随着中国“深空、深海、深地”战略的推进以及“一带一路”国际合作的拓展，人们对重力勘探的精度提出了更高要求，同时也需要更强的处理大规模数据集的能力以及更好的应对复杂问题的能力（Ren et al., 2023）。在此背景下，作为重力数据解释和反演成像核心的3D正演建模已成为一个前沿研究方向，提升其精度和效率对于突破深地探测中的瓶颈至关重要。重力正演建模的理论基础是控制重力势的泊松方程，其数学本质是在狄利克雷或罗宾边界条件下求解偏微分方程（PDE）（Blakely, 1996）。当前主流的3D正演建模方法大致可分为空间域方法和傅里叶域（或频域）方法，这种分类是根据求解控制方程所在的域来确定的。这些方法在计算效率、边界条件处理以及复杂模型适应性方面各有优势，共同推动了重力勘探向深地及深海领域的深入发展。

空间域解法通常分为基于积分公式的前向计算（解析方法）和基于微分公式的前向计算（数值方法）。解析方法是指通过累加各种基本体（如棱柱体）的贡献来计算重力效应，其中的积分要么以封闭形式求解，要么通过数值方法进行最终积分。然后利用解析解的线性叠加原理，即可快速且高精度地计算出三维复杂密度分布下的重力势、重力场矢量以及重力梯度张量。空间域解析方法的典型策略包括：将计算域分解为多个矩形棱柱体（Li and Chouteau, 1998; Nagy et al., 2000），将计算域分解为多个多面体元素（Petrovic, 1996; Ren et al., 2017; Chen et al., 2018），以及直接采用三维数值积分法（Asgharzadeh et al., 2007; Roussel et al., 2015）。基于积分的方法存在一个显著的计算缺陷，即被积函数取决于观测点的位置。因此，每个观测点都需要进行一次体积积分，这就导致计算时间与观测点数量成线性关系。对于具有大量观测点且密度分布范围较大的情况，这将带来极高的计算成本（Chen et al., 2012）。相比之下，基于微分表达式的正演建模原理将问题转化为由重力势控制的边值问题（BVP），进而计算其他衍生量（如重力场矢量或重力梯度张量）。这一类别的主要方法包括有限差分法（FD）（Farquharson and Mosher, 2009; Pan et al., 2021; Hu et al., 2023）、有限元法（Cai and Wang, 2005; Tong et al., 2024）、谱元法（Gharti et al., 2018; Martin et al., 2017）以及无网格法（Long and Farquharson, 2019）。基于微分方程的正演方法能更有效地适应源物理特性的变化。然而，对于网格离散化程度很高的大规模模型，会出现一个关键限制：计算成本会急剧上升，同时内存消耗也会大幅增加。

傅里叶域正演建模是通过对空间域的解析或近似解析表达式应用积分变换（如傅里叶变换、拉普拉斯变换和小波变换）来获得它们的频谱，随后通过这些频谱的逆变换得出正演响应值（Bhattacharyya and Navolio, 1976）。与空间域方法相比，傅里叶域重力正演建模在数学简洁性和计算效率方面具有明显优势。不过，其精度会受到吉布斯现象以及离散傅里叶变换中的波数截断误差的影响。快速傅里叶变换算法的出现是一个重要的进步，它将传统DFT的计算复杂度从O(N2)降低到了O(N log N)，极大地加快了波数域与空间域之间的数据转换速度。因此，基于频域的3D重力响应正演建模能够在保持计算精度的同时，高效地模拟复杂的地质模型（Tontini et al., 2009; Chen and Liu, 2019; Yin et al., 2023; Wan et al., 2024; Wang et al., 2024）。经典FFT算法的一个众所周知的局限性是，在处理被有限窗口截断的数据时会产生吉布斯振荡伪影，这对于非周期信号而言表现为边界能量泄漏，这就需要采用边界扩展方法来减轻频谱失真（Wu and Tian, 2014）。不过，这些缓解策略会在一定程度上降低正演建模的计算效率。为了解决这个问题，Wu和Tian（2014）提出了高斯-FFT技术，并证明了其在傅里叶域正演建模中的优异性能（Wu, 2016; Zhao et al., 2018; Roy and Wu, 2023），该技术能有效避免零波数处的除零奇点问题。通过将物理场问题转化为频域进行求解，可以利用频域方法大幅提高计算效率，但其缺点在于属于近似计算而非精确计算，而且仅适用于规则网格上的正演建模。

最近，Dai等人（2020, 2022, 2024）提出了一种创新的混合傅里叶谱方法。该方法的理论框架利用了2D高斯-FFT在水平方向上的波数解耦特性，能够将原本强耦合的三维边值问题分解为多个针对不同波数的独立一维常微分方程边值问题，从而提升计算效率。与直接在空间域求解重力势的边值问题相比，这种方法能实现更快的计算速度。不过，在数值求解这些一维常微分方程边值问题时，仍然需要求解一组线性方程（Dai et al., 2022）。此外，水平方向上应用2D傅里叶变换本质上也是基于傅里叶谱方法的概念，这就要求必须遵循相应的周期性边界条件，这对获取大规模问题的高精度正演响应构成了巨大挑战。离散正弦变换是一种傅里叶变换的变体，广泛用于信号处理、图像压缩和数值分析领域。在偏微分方程的数值求解中，经常运用离散正弦变换来简化边界条件的处理并加快计算速度，进而提升计算效率（Schumann and Sweet, 1988; Feng and Zhao, 2020; Morin and Paxus, 2024）。通过巧妙选择并将离散正弦变换技术与有限差分法相结合，可以显著提高求解偏微分方程的效率（Wang et al., 2016; Gatiso et al., 2021; Paux et al., 2025; Tong et al., 2025）。这种方法为大规模地球物理正演建模问题的快速计算带来了巨大希望。

重力势及其导数的三维数值模拟面临着双重计算瓶颈：由于全空间离散化，计算复杂度呈数量级上升，而内存存储需求则呈指数级增长。这些挑战使得高效、精确的大尺度重力正演建模和反演工作变得极为困难。为了解决这一计算问题，本文提出了一种创新的混合离散正弦变换-有限差分（DST-FD）耦合算法用于数值模拟。该算法的核心机制是在水平方向上应用2D离散正弦变换，这种变换将重力势的3D有限差分公式转化为空间域中一系列相互独立的1D垂直方向常差分方程，每一个方程对应一个不同的波数。这种方法不仅显著提升了正演计算速度，大幅降低了计算机内存存储需求，同时还避免了傅里叶谱方法中固有的周期性边界问题。通过设计一个球形密度模型，并将重力势及其导数的数值解与其解析结果进行比较，可以验证该方法的准确性和适用性。通过SEG/EAGE盐穹模型的数值模拟，检验了混合DST-FD算法在大型正演建模中的有效性。最后，利用长白山的实际地形数据对重力场进行高效的正演计算，结果表明我们的正演建模算法为极端地形区域的深度勘探数值模拟提供了一种新的方法。

章节节选
重力势的三维边值问题
根据牛顿万有引力定律和高斯定理，空间域中密度分布所产生的重力势遵循以下三维泊松方程（Blakely, 1996）：
?2Uxyz/?x2 + ?2Uxyz/?y2 + ?2Uxyz/?z2 = ?4πGρxyz
其中，(x, y, z)表示三维标准笛卡尔坐标系，Uxyz表示重力势，ρxyz表示残余密度，G=6.672×10?11N?m2/kg2则是经验物理常数。

方程（1）的解

数值实验
我们利用非常适合数值计算的软件环境MATLAB，开发了一套用于三维重力正演建模的混合DST–FD程序。为了验证所提出的正演算法的正确性和性能，我们在球形模型以及三维SEG/EAGE盐穹模型上进行了数值实验。本研究中的所有数值模拟都是在一台配备64GB内存的台式电脑上完成的。

重力场的三个分量是通过对……进行近似得到的

实际地形的应用
为了验证所提出的正演建模算法在大型复杂模型中的工程适用性，我们在中国长白山的一个代表性区域开展了实证研究，研究中使用了数字高程模型。如图12所示，研究区域的范围为50公里×50公里×1.769公里，该区域的密度设定为2580千克/立方米，而周边区域的密度则设为0千克/立方米。该区域的高程从……变化至……

结论
利用空间重力势满足非齐次狄利克雷边界条件的特性，我们提出了一种利用DST加速有限差分法的数值模拟方法，即混合DST-FD方法。该方法的3D计算域被离散化为直角矩形棱柱单元，其中x方向和y方向的单元尺寸均匀，而z方向的单元尺寸则不均匀。控制重力势的泊松方程……

作者贡献声明
胡一遥：方法论、形式分析、研究工作、初稿撰写、审稿与编辑。童晓忠：概念构建、方法论、项目管理、资金筹集、监督工作、审稿与编辑。苏杰宇：研究工作、形式分析、资源协调、审稿与编辑。魏康贵：数据整理、验证、软件应用、形式分析、可视化处理。

利益冲突声明
作者声明他们没有已知的可能会影响本文所述工作的财务利益或个人关系。

致谢
本项工作得到了中国国家深地重大科技项目（编号2025ZD1009704）、湖南省自然科学基金（编号2023JJ30659）以及大学生创新创业训练计划项目（编号S202510533440）的支持。作者们还要感谢中国海洋大学刘颖先生在讨论中所给予的宝贵帮助。同时，我们也感谢主编V. Paoletti以及两位匿名审稿人的辛勤工作。