基于径向多项式密度Tesseroid体的大尺度重力场快速正演计算

Fast gravitational forward modelling at large scales using tesseroids with polynomial density in depth

在大尺度重力模拟中,为了考虑地球曲率的影响,通常将质量体沿经度、纬度及径向方向剖分为多个Tesseroid体,然后根据叠加原理将所有Tesseroid体重力效应之和近似为质量体的总重力场.在现有大多数方法中,总重力场的计算一般通过将各Tesseroid体的重力效应直接累加获得,尽管这种做法十分简单,但对于大尺度的复杂密度模型,其计算量十分巨大.为此,本文充分利用重力位、重力加速度及重力梯度张量核函数积分在经度方向的卷积特性,提出一种基于快速离散卷积算法和优化重力公式的高精度大尺度重力场模拟方法.该方法采用具有径向任意阶多项式密度的Tesseroid体来描述复杂密度分布,并以变密度Grombein优化公式为基础进行重力场正演模拟,从而有效克服经典球坐标系重力公式引起的极轴观测点奇异问题;同时联合Gauss-Legendre求积(GLQ)和水平方向的自适应网格剖分,实现重力积分的高精度计算,最后结合快速离散卷积算法获得质量体的重力效应.本文方法的加速策略通过如下途径实现:①假定经度方向的Tesseroid体具有相同尺寸,观测点位于同一高度且其采样间隔与Tesseroid体尺寸相等,由此建立与重力积分核对应的Toeplitz权系数矩阵,并利用Toeplitz对角线元素相等的特性,计算其第一列和第一行元素来复原完整矩阵,以此大大降低计算内存及计算时间;②利用基于FFT的离散卷积算法实现权系数矩阵与密度矩阵乘积的快速计算,从而进一步提高计算效率.数值算例表明,本文方法可实现近地表至卫星高度重力位、重力加速度及重力梯度张量的高精度计算,且其计算效率相较于直接累加法高出约3~4个数量级.该方法在WINTERC-G全球模型重力效应计算中的应用进一步证实了算法的可靠性和灵活性.本文快速算法是经典3D-GLQ方法的重要补充,对研究大尺度重力正反演问题、地形效应及地球内部结构等有重要意义.

In large-scale gravitational forward modelling,mass bodies are usually divided into a number of tesseroids along the longitudinal,latitudinal and radial directions to take the effect of the Earth's curvature into account.The total effect is then obtained as a cumulative contribution of all tesseroids in terms of superposition.In most existing methods,the calculation of the total gravitational effects is generally achieved by directly adding the gravity response of each tesseroid together.Although this approach is very simple,the computational cost becomes extremely large when the density model is complex.In this study,we make full use of the convolution relation of the Newton's kernels with respect to the longitude,and propose an accurate and efficient method for large-scale gravitational forward modelling based on the fast discrete convolution algorithm and the well-established optimized formulas for the gravitational field.In the proposed method,the density in each tesseroid is assumed to be polynomial in depth to adapt to a complex density environment,and the forward modeling is carried out based on the Grombein's optimized formulas with a vertically variable density,so as to effectively overcome the polar singularity induced by the spherical coordinate system.In addition,the Gauss Legendre quadrature(GLQ)and the 2D adaptive discretization strategy are combined to ensure the accuracy of the numerical integration of the Newton's integral.The fast discrete convolution algorithm is then utilized to speed up the forward calculation.The key points of the acceleration technique involves:①In the longitudinal direction,the tesseroids should have the same size,the observation points are equally spaced and the sampling interval is equal to the size of the tesseroid,so that the weight coefficient matrix can be established into a general Toeplitz form.Considering the fact that the diagonal elements of a Toeplitz are constant,the whole weight coefficient matrix can then be restored by only calculating the elements in its first column and first row.As a result,the computation time and memory cost are significantly reduced;②The product of the weight coefficient matrix and density matrix is then efficiently calculated by using the discrete convolution algorithm based on FFT.This further improves the computation speed.Numerical tests show that the proposed method can obtain accurate results for the gravitational potential,the gravitational acceleration,and the gravitational gradient tensor from near surface to satellite height,and its computational efficiency is about 3~4 orders of magnitude higher than the traditional 3D-GLQ method based on pure summation.The application to the WINTERC-G global model further confirms the reliability and flexibility of the proposed method.Our algorithm is an important supplement to the classical 3D-GLQ method,and offers a significant tool for solving large-scale gravity forward and inverse problems,analyzing topographic effects,and understanding the internal structure of the Earth.