Fast gravitational forward modelling at large scales using tesseroids with polynomial density in depth
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.
Regional and global scalesGravitational modellingFFT-based accelerationTesseroid with polynomial density in depthSpherical coordinates