Fast and accurate calculation of gravitational effects on a regional or global scale with complex density environment is a critical issue in gravitational forward modelling. Most existing significant developments with tessroid-based modelling are limited to homogeneous density models or polynomial ones of a limited order. Moreover, the total gravitational effects of tesseroids are often calculated by pure summation in these methods, which makes the calculation extremely time-consuming. A new efficient and accurate method based on tesseroids with a polynomial density up to an arbitrary order in depth is developed for 3D large-scale gravitational forward modelling. The method divides the source region into a number of tesseroids, and the density in each tesseroid is assumed to be a polynomial function of arbitrary degree. To guarantee the computational accuracy and efficiency, two key points are involved: (1) the volume Newton’s integral is decomposed into a one-dimensional integral with a polynomial density in the radial direction, for which a simple analytical recursive formula is derived for efficient calculation, and a surface integral over the horizontal directions evaluated by the Gauss–Legendre quadrature (GLQ) combined with a 2D adaptive discretization strategy; (2) a fast and flexible discrete convolution algorithm based on 1D fast Fourier transform (FFT) and a general Toepritz form of weight coefficient matrices is adopted in the longitudinal dimension to speed up the computation of the cumulative contributions from all tesseroids. Numerical examples show that the gravitational fields predicted by the new method have a good agreement with the corresponding analytical solutions for spherical shell models with both polynomial and non-polynomial density variations in depth. Compared with the 3D GLQ methods, the new algorithm is computationally more accurate and efficient. The calculation time is significantly reduced by 3 orders of magnitude as compared with the traditional 3D GLQ methods. Application of the new algorithm in the global crustal CRUST1.0 model further verifies its reliability and practicability in real cases. The proposed method will provide a powerful numerical tool for large-scale gravity modelling and also an efficient forward engine for inversion and continuation problems.

在具有复杂密度环境的区域或全球范围内快速准确地计算引力效应是引力正演建模中的关键问题。基于 tessroid 的建模的大多数现有重大发展仅限于均匀密度模型或有限阶的多项式模型。而且，在这些方法中，立方体的总引力效应往往是通过纯求和来计算的，这使得计算极其耗时。为 3D 大规模引力正演建模开发了一种基于多项式密度高达任意深度的 tesseroids 的新的有效和准确的方法。该方法将源区域划分为多个tesseroids，每个tesseroid中的密度被假定为任意阶的多项式函数。为了保证计算的准确性和效率，涉及到两个关键点：（1）将体积牛顿积分分解为径向多项式密度的一维积分，推导出简单的解析递推公式以实现高效计算, 以及由高斯-勒让德积分 (GLQ) 结合 2D 自适应离散化策略评估的水平方向上的表面积分；(2) 在纵向维度上采用基于一维快速傅里叶变换 (FFT) 和一般托普利茨形式的权重系数矩阵的快速灵活的离散卷积算法，以加速所有曲面的累积贡献的计算。算例表明，新方法预测的引力场与具有多项式和非多项式密度深度变化的球壳模型的相应解析解具有良好的一致性。与 3D GLQ 方法相比，新算法在计算上更加准确和高效。与传统的 3D GLQ 方法相比，计算时间显着减少了 3 个数量级。新算法在全球地壳CRUST1.0模型中的应用进一步验证了其在实际案例中的可靠性和实用性。所提出的方法将为大规模重力建模提供强大的数值工具，也为反演和延拓问题提供高效的正演引擎。