该文档属于类型a，是一篇报告单一原创研究的学术论文。
关于复杂条件二维重力场及重力张量场空间波数域正演方法的研究报告
一、 研究作者、机构与发表信息
本研究由中南大学地球科学与信息物理学院的李颖梅、戴世坤、李昆、陈轻蕊和凌嘉宣共同完成。该研究成果以论文形式发表在《物探化探计算技术》（computing techniques for geophysical and geochemical exploration）期刊2019年3月出版的第41卷第2期上。
二、 研究的学术背景与目标
本研究属于应用地球物理学领域，具体聚焦于重力和重力梯度（重力张量）勘探中的正演模拟核心问题。随着传感器技术的飞速发展，重力场与重力张量场的测量技术日益精进，为反演揭示地下精细的密度结构提供了更丰富、更高质量的数据基础。然而，正演计算作为反演解释的基石，其效率和精度直接决定了后续反演和人机交互解释的可行性与效果。
传统上，针对二度体的重力及重力张量正演计算方法主要分为解析法和数值法两大类。解析法对于形状规则、密度分布简单的模型有效，但面对起伏地形、任意复杂截面形状及任意密度分布等“复杂条件”时，其解析解的推导变得极其困难甚至不可能，且存在奇异点计算和数值不稳定问题。数值法，如有限元法或有限差分法等微分法，以及基于规则体叠加或直接数值积分的积分法，虽然适应性更强，但往往计算效率较低，难以满足大规模、高分辨率反演对正演计算速度的迫切需求。谱方法（如快速傅里叶变换法）利用卷积定理在频率域进行快速计算，但传统的FFT方法在处理有限区域问题时存在截断效应（边界效应）和零波数奇异问题，影响计算精度。
因此，本研究旨在解决一个关键的科学与工程问题：如何实现在任意密度分布和起伏地形观测条件下，对复杂二维地质体进行重力场及重力张量场的高效、高精度正演计算。其具体目标包括：1) 提出一种新的、适用于复杂条件的空间-波数混合域正演算法；2) 该算法需兼具高计算效率与高数值精度；3) 为后续实现高效的精细化重力反演及人机交互解释提供可靠的技术支撑。
三、 研究的详细工作流程
本研究的工作流程是一个从理论推导、算法设计到数值验证的完整链条，主要包括以下几个核心步骤：
步骤一：复杂地质体模型的离散化与问题转化。 研究首先将任意截面形状、任意密度分布的复杂二度体剖分为一系列小的矩形单元（网格）。每个矩形单元内的密度被视为常数值，通过调整网格大小和密度赋值，可以灵活且精细地逼近任意复杂的地质构造。这种剖分方式简单直观，特别适合于反演和人机交互建模。基于重力异常的叠加原理，整个复杂二度体产生的重力场及重力张量场，即可转化为所有这些小矩形二度体（即“矩形二度体组合模型”）产生的场的总和。研究给出了这一组合模型在空间域计算的精确积分表达式（论文中的公式1至4），这构成了后续算法推导的出发点。
步骤二：核心算法的理论推导——空间-波数混合域表达式的建立。 这是本研究最核心的创新环节。为了将二维空间域的卷积计算转化为更高效的形式，研究对步骤一中得到的空间域积分公式（公式1-4）的两边，在水平方向（x方向）施加了一维傅里叶变换。通过一系列严谨的数学推导，研究团队得到了重力垂直分量g_z、水平分量g_x以及重力张量分量g_zx和g_xx在空间-波数混合域的全新解析表达式（论文中的公式5至8）。这些公式的特点是：在波数（k）域，场值表示为关于深度（z）的函数，且表达式中包含了对矩形单元在垂直方向（z方向）积分的解析形式。特别地，公式中出现了密度函数的一维傅里叶变换形式（公式9），这为后续利用快速傅里叶变换（FFT）加速计算埋下了伏笔。
步骤三：高效高精度数值算法的实现——Gauss-FFT与样条插值。 获得混合域表达式后，需要将其转换回最终需要的空间域场值，即进行一维傅里叶反变换。本研究在此关键环节引入了Gauss-FFT方法。这是一种高精度的数值积分方法，用于计算傅里叶反变换积分。与传统的标准FFT扩边法相比，Gauss-FFT方法避免了零波数的奇异计算问题，并且能有效克服传统FFT法因数据截断产生的边界效应，从而在保证计算效率（比标准FFT扩边法更快）的同时，显著提高了计算精度。这是本研究方法实现“高精度”目标的关键技术之一。
对于起伏地形上的观测，研究提出了一种巧妙的策略：首先，根据地形的起伏范围，确定观测区域内从最低点到最高点之间的一系列不同高度的水平网格线。然后，利用上述基于Gauss-FFT的新算法，高效计算出所有这些水平网格线上的重力场及重力张量场值。最后，对于实际起伏地形上的任意观测点，采用三次样条插值方法，根据其水平坐标和高度，从上下相邻的两条水平网格线的计算结果中进行插值，从而获得该起伏观测点上的场值。这一策略的优势在于，它将复杂的、非规则观测面的正演问题，转化为多个规则（水平）观测面的正演与一维插值问题，既利用了新算法在水平观测下的高效高精度特性，又通过成熟的插值技术解决了地形起伏问题，避免了直接计算非规则网格的复杂性。
步骤四：算法验证与性能测试。 为了验证所提算法的正确性、精度和对复杂条件的适应性，研究设计了两个模型算例。第一个算例使用一个简单的矩形均匀密度体模型，在水平观测线上计算其重力场和重力张量场，并将新算法的数值解与已知的解析解进行对比，通过绘制曲线图和计算相对均方根误差来定量评估精度。第二个算例则模拟了更为现实的场景：一个具有正弦起伏的观测地形。研究通过改变用于插值的水平网格线的数量（从3条到更多条），系统分析了起伏地形上计算精度随插值密度变化的规律，并统计了起伏线上数值解相对于解析解的相对均方根误差。所有计算均在指定配置的计算机上完成，并记录了计算时间，以评估算法的效率。
四、 研究的主要结果
本研究通过理论推导和数值实验，获得了以下关键结果：
在水平观测条件下的精度验证结果： 对于网格剖分为100×100的矩形模型，在水平测线上，新算法计算的重力场（g_x, g_z）数值解与解析解高度吻合。相对均方根误差非常小，其中重力场的最大相对误差仅为0.004%，重力张量场（g_xx, g_zx）的最大相对误差约为0.4%。误差分布显示，由于截断效应，场值在计算区域边界处误差稍大；同时，重力张量由于信号变化更剧烈，其计算精度总体上略低于重力场，这与理论预期一致。这些数据强有力地证明了本文提出的空间-波数混合域算法结合Gauss-FFT反变换，在水平观测情况下具有极高的计算精度。
在起伏地形条件下的算法性能结果： 针对正弦起伏地形，研究发现起伏地形上的计算精度（通过相对均方根误差衡量）随着用于插值的水平网格线数量的增加而迅速提高并趋于稳定。当仅使用3条水平线进行插值时，起伏地形上重力场的误差已小于0.5%，重力张量场误差约为2.5%。当水平网格线数量增加到11条以上时，误差几乎不再下降，表明计算精度已达到最优。在采用15条水平线进行插值计算时，整个算法（包括计算所有水平线场值和插值）的耗时仅为0.020秒。这一结果验证了“水平计算+插值”策略对于处理起伏地形的有效性和高效性。图6展示了当插值线为11条时，起伏地形上数值解与解析解的曲线对比，两者吻合度极高，直观地证明了方法的可靠性。
算法的计算效率与可扩展性： 研究明确指出，该算法的计算量主要与模型的网格剖分数量相关，且由于核心计算在波数域分解为多个独立的一维积分，算法具有良好的并行潜力。测试表明，在网格数为100×100时计算迅速；当网格数增大到1000×1000时，计算时间约为1.3秒，显示出算法耗时随网格数量近似线性增长的良好可扩展性，这为处理大规模、高分辨率的实际地质问题奠定了基础。
实际模型应用展示： 为进一步展示算法的实用性，研究应用该算法计算了文献中给出的“格尔木—花海子剖面”岩石圈二维密度结构模型（图7）在地表水平观测线上的重力场和重力梯度（图8）。虽然此部分未与解析解对比，但成功计算出了符合模型特征的理论异常曲线，展示了该方法处理复杂密度分布模型的能力。
五、 研究的结论与价值
本研究成功提出并验证了一种基于空间-波数混合域表达式、Gauss-FFT反变换和三次样条插值的复杂条件二维重力场及重力张量场高效高精度正演算法。主要结论和价值体现在： 1. 科学方法论价值： 提出了一种新的正演计算框架。通过将二维空间卷积拆分为波数域多个一维积分的累加，并利用Gauss-FFT进行高精度反变换，从根本上改进了传统谱方法，实现了效率与精度的平衡。 2. 技术应用价值： 该算法能够有效应对起伏地形和任意密度分布这两大复杂条件，计算结果精度高、速度快。为利用高精度重力及重力梯度数据进行精细化的地下密度结构反演提供了高效可靠的正演工具，是推动重力勘探走向定量化、精细化解释的关键技术环节。 3. 扩展性价值： 论文指出，该方法的核心思想（矩形剖分、波数域转化、Gauss-FFT）还可推广应用于三维重磁问题的求解，展示了其方法论的普适潜力。
六、 研究的亮点
七、 其他有价值的内容
研究在引言部分对重力及重力张量正演方法的发展脉络进行了清晰的梳理，将现有方法分为解析法和数值法（微分法与积分法），并指出了各类方法的优缺点，为本文工作的定位和创新点的阐述提供了坚实的学术背景。此外，文献综述部分引用了大量国内外相关研究，显示了研究工作的继承性与前沿性。模型算例的设计从简单到复杂（水平线->起伏线），从精度验证到性能测试，逻辑严谨，论证充分。最后，将算法应用于一个接近实际的岩石圈密度结构模型，增强了研究的应用说服力。