Optimized formulas for the gravitational field of a tesseroid

Various tasks in geodesy, geophysics, and related geosciences require precise information on the impact of mass distributions on gravity field-related quantities, such as the gravitational potential and its partial derivatives. Using forward modeling based on Newton’s integral, mass distributions are generally decomposed into regular elementary bodies. In classical approaches, prisms or point mass approximations are mostly utilized. Considering the effect of the sphericity of the Earth, alternative mass modeling methods based on tesseroid bodies (spherical prisms) should be taken into account, particularly in regional and global applications. Expressions for the gravitational field of a point mass are relatively simple when formulated in Cartesian coordinates. In the case of integrating over a tesseroid volume bounded by geocentric spherical coordinates, it will be shown that it is also beneficial to represent the integral kernel in terms of Cartesian coordinates. This considerably simplifies the determination of the tesseroid’s potential derivatives in comparison with previously published methodologies that make use of integral kernels expressed in spherical coordinates. Based on this idea, optimized formulas for the gravitational potential of a homogeneous tesseroid and its derivatives up to second-order are elaborated in this paper. These new formulas do not suffer from the polar singularity of the spherical coordinate system and can, therefore, be evaluated for any position on the globe. Since integrals over tesseroid volumes cannot be solved analytically, the numerical evaluation is achieved by means of expanding the integral kernel in a Taylor series with fourth-order error in the spatial coordinates of the integration point. As the structure of the Cartesian integral kernel is substantially simplified, Taylor coefficients can be represented in a compact and computationally attractive form. Thus, the use of the optimized tesseroid formulas particularly benefits from a significant decrease in computation time by about 45 % compared to previously used algorithms. In order to show the computational efficiency and to validate the mathematical derivations, the new tesseroid formulas are applied to two realistic numerical experiments and are compared to previously published tesseroid methods and the conventional prism approach.     

Calculation of strongly singular and hypersingular surface integrals

Efficient numerical computation of integrals defined on closed surfaces in ℝ³ with non-integrable point singularities that arise in physical geodesy is discussed. The method is based on the use of polar coordinates and the definition of integrals with non-integrable point singularities as Hadamard finite part integrals. First the behavior of singular integrals under smooth parameter transformations is studied, and then it is shown how they can be reduced to absolutely integrable functions over domains in ℝ². The correction terms that usually arise if the substitution rule is formally applied, in contrast to absolutely integrable functions, are calculated. It is shown how to compute the regularized integrals efficiently, and, numerical efforts for various orders of singularity are compared. Finally, efficient numerical integration methods are discussed for integrals of functions that are defined as singular integrals, a task that typically arises in Galerkin boundary element methods. Received: 15 April 1997 / Accepted: 7 May 1998     

Representation of planar integral-transformations by 4-D wavelet decomposition

In one way or the other, numerical methods for the evaluation of integral operators can often be related to the solution of the so-called Galerkin equations. For convolution operators and exponentials with purely imaginary exponents as base functions the Galerkin matrix becomes diagonal and this fact is the core of the FFT techniques, used in Physical Geodesy. For non-convolution operators the FFT technique is not applicable. This paper aims at the development of a technique, which can also be applied for non-convolution operators. This technique is based on the use of wavelets as base functions. In this case the Galerkin matrix is not diagonal but (after thresholding) very sparse and this leads to methods, which are similarly efficient as FFT in the convolution case. The paper starts with the theoretical background for n-dimensional wavelet analysis and the representation of integral operators with respect to those wavelet bases. The resulting algorithm is tested for convolution and non-convolution operators.     

关于重心基准的坐标转换方法研究

基准统一与转换是测量工程中的重要问题,文中通过数值计算表明,为了克服坐标转换模型中参数之间的强相关问题,采用坐标重心化策略,得到的结果与四参数方法基本相同,重心化方法减小了参与计算的坐标数值,增强数值计算的可靠性,随着现代计算技术(计算设备)的广泛应用,两者的差异并不明显。各种坐标转换算法的可靠性依赖于公共点的精度与公共点的分布形式,在极坐标方法中尤其要避免边长悬殊的短边参与参数解算。      

近区地形直接与间接影响的棱柱模型算法

基于Helmert第二压缩法进行边值解算时需要计算地形压缩对重力的直接影响和对（似）大地水准面的间接影响。计算近区直接、间接影响的传统积分算法仍是二重积分形式。该算法以网格中心点处的积分核作为网格积分核的平均值的计算模式在一定程度上引入了近似误差。另外,直接、间接影响的传统积分算法在中央区存在奇异性,需单独计算中央网格地形影响,因而增加了计算的复杂性。为此,本文推导了近区地形直接、间接影响的棱柱模型公式,一方面提高了地形影响的计算精度;另一方面中央区不存在奇异性,从而简化了计算过程。为避免棱柱模型存在的平面近似误差,可使用顾及地球曲率的棱柱模型算法计算地形影响。最后通过试验得出结论,在（似）大地水准面精度要求较高的应用中,应尽量使用顾及地球曲率的棱柱模型算法计算地形影响。     

On the accurate numerical evaluation of geodetic convolution integrals

In the numerical evaluation of geodetic convolution integrals, whether by quadrature or discrete/fast Fourier transform (D/FFT) techniques, the integration kernel is sometimes computed at the centre of the discretised grid cells. For singular kernels—a common case in physical geodesy—this approximation produces significant errors near the computation point, where the kernel changes rapidly across the cell. Rigorously, mean kernels across each whole cell are required. We present one numerical and one analytical method capable of providing estimates of mean kernels for convolution integrals. The numerical method is based on Gauss-Legendre quadrature (GLQ) as efficient integration technique. The analytical approach is based on kernel weighting factors, computed in planar approximation close to the computation point, and used to convert non-planar kernels from point to mean representation. A numerical study exemplifies the benefits of using mean kernels in Stokes’s integral. The method is validated using closed-loop tests based on the EGM2008 global gravity model, revealing that using mean kernels instead of point kernels reduces numerical integration errors by a factor of ~5 (at a grid-resolution of 10 arc min). Analytical mean kernel solutions are then derived for 14 other commonly used geodetic convolution integrals: Hotine, Eötvös, Green-Molodensky, tidal displacement, ocean tide loading, deflection-geoid, Vening-Meinesz, inverse Vening-Meinesz, inverse Stokes, inverse Hotine, terrain correction, primary indirect effect, Molodensky’s G1 term and the Poisson integral. We recommend that mean kernels be used to accurately evaluate geodetic convolution integrals, and the two methods presented here are effective and easy to implement.     

重力异常向上延拓严密改化模型及向下延拓应用

重力异常向上延拓全球积分模型在航空重力测量数据质量评估和向下延拓迭代计算等领域具有广泛的应用。为了消除积分核函数奇异性影响,需要对该模型进行基于积分恒等式的移去-恢复转换及全球积分域的分区改化处理。在此过程中,传统改化处理方法往往忽略了全球积分过渡到局域积分引起的积分恒等式偏差影响,从而导致不必要的计算模型误差,最终影响向上延拓计算结果的可靠性,甚至影响向下延拓迭代解算结果的稳定性。针对此问题,本文开展了重力异常向上延拓积分模型改化及向下延拓应用分析研究,依据实测数据保障条件和积分恒等式适用条件要求,导出了重力异常向上延拓积分模型的分步改化公式,提出了补偿传统改化模型缺陷的修正公式,并将最终的严密改化模型应用于重力异常向下延拓迭代解算。使用超高阶地球位模型EGM2008作为标准位场开展数值计算检验,分别对重力异常向上延拓分步改化模型的计算精度及在向下延拓迭代解算中的应用效果进行了检核评估,验证了采用严密改化模型的必要性和有效性。     

四种改进积分法的低空扰动引力计算

针对Stokes积分方法计算扰动引力中计算点从空中趋近地面时存在积分奇异和不连续的问题,该文提出了去中央奇异点法、奇异点积分值修正法、中央格网加密算法和改进积分式法4种改进Stokes积分的计算公式,并进行了实验计算。计算结果表明:近地空间范围内,4种改进算法都能在一定程度上改进原始积分的奇异性问题;相同条件下,奇异点积分值修正法和改进积分式法计算精度最高,适宜于低空计算;改进积分式法通过理论推导,得到了从球外部到球面统一、连续且无奇异的改进Stokes积分公式,理论严谨。     

Is Newton's iteration faster than simple iteration for transformation between geocentric and geodetic coordinates?

Summary Two iterative algorithms for transformation from geocentric to geodetic coordinates are compared for numerical efficiency: the well known Bowring's algorithm of 1976, which employs the method of simple iteration, and the recent (1989) algorithm by Borkowski, which employs the Newton-Raphson method. The results of numerical tests suggest that the simple iteration method implemented in Bowring's algorithm executes approximately 30% faster than the Newton-Raphson method implemented in Borkowski's algorithm. Only two iterations of each algorithm are considered. Two iterations are sufficient to produce coordinates accurate to the comparable level of 1E-9 m, which exceeds the requirements of any practical application. Therefore, in the class of iterative methods, the classical Bowring's algorithm should be the method of choice.     

Truncated geoid and gravity inversion for one point-mass anomaly

The truncated geoid, defined by the truncated Stokes' integral transform, an integral convolution of gravity anomalies with the Stokes' function on a spherical cap, is often used as a mathematical tool in geoid computations via Stokes' integral to overcome computational difficulties, particularly the need to integrate over the entire boundary spheroid. The objective of this paper is to demonstrate that the truncated geoid does, besides having mathematical applications, have physical interpretation, and thus may be used in gravity inversion. A very simple model of one point-mass anomaly is chosen and a method for inverting its synthetic gravity field with the use of the truncated geoid is presented. The method of inverting the synthetic field generated by one point-mass anomaly has become fundamental for the authors' inversion studies for sets of point-mass anomalies, which are published in a separate paper. More general applications are currently under investigation. Since an inversion technique for physically meaningful mass distributions based on the truncated geoid has not yet been developed, this work is not related to any of the existing gravity inversion techniques. The inversion for one point mass is based on the onset of the so-called dimple event, which occurs in the sequence of surfaces (or profiles) of the first derivative of the truncated geoid with respect to the truncation parameter (radius of the integration cap), its only free parameter. Computing the truncated geoid at various values of the truncation parameter may be understood as spatial filtering of surface gravity data, a type of weighted spherical windowing method. Studying the change of the truncated geoid represented by its first derivative may be understood as a data enhancement method. The instant of the dimple onset is practically independent of the mass of the point anomaly and linearly dependent on its depth. Received: 26 September 1996 /Accepted: 28 September 1998     

不规则形状地下空间开挖条件下地表沉陷预计方法研究

地下工程施工或地下采矿导致地表沉陷，甚至造成突然坍塌或沿层面滑移等灾害事故。文中根据概率积分法的基本原理，基于褶曲构造地层，任意形状空间开挖条件下地表点在任意方向的移动与变形值预计方法，并顾及复杂地质和开挖等因素，推导出不规则形状地下空间开挖条件下地表移动与变形预计公式。由于地表沉陷预计公式都是不可积重积分式，因此，采用数值积分方法解算这些不能求出原函数的方程式和变积分限问题。     

复合Simpson公式在线路中边桩坐标计算中的应用

基于线路线任一点切线方位角计算通式，推出了线路中边桩坐标计算的积分通用公式；借助于复合Simpson公式，给出了线路中边桩坐标计算的通用数值模型，并详细讨论了积分区间等分数n的取值问题。通用模型的建立，为线路中边桩的极坐标设法奠定了坚实的基础。     

A comparison of methods for the inversion of airborne gravity data

Four integral-based methods for the inversion of gravity disturbances, derived from airborne gravity measurements, into the disturbing potential on the Bjerhammar sphere and the Earths surface are investigated and compared with least-squares (LS) collocation. The performance of the methods is numerically investigated using noise-free and noisy observations, which have been generated using a synthetic gravity field model. It is found that advanced interpolation of gravity disturbances at the nodes of higher-order numerical integration formulas significantly improves the performance of the integral-based methods. This is preferable to the commonly used one-point composed Newton–Cotes integration formulas, which intrinsically imply a piecewise constant interpolation over a patch centered at the observation point. It is shown that the investigated methods behave similarly for noise-free observations, but differently for noisy observations. The best results in terms of root-mean-square (RMS) height-anomaly errors are obtained when the gravity disturbances are first downward continued (inverse Poisson integral) and then transformed into potential values (Hotine integral). The latter has a strong smoothing effect, which damps high-frequency errors inherent in the downward-continued gravity disturbances. An integral method based on the single-layer representation of the disturbing potential shows a similar performance. This representation has the advantage that it can be used directly on surfaces with non-spherical geometry, whereas classical integral-based methods require an additional step if gravity field functionals have to be computed on non-spherical geometries. It is shown that defining the single-layer density on the Bjerhammar sphere gives results with the same quality as obtained when using the Earths topography as support for the single-layer density. A comparison of the four integral-based methods with LS collocation shows that the latter method performs slightly better in terms of RMS height-anomaly errors.     

Ellipsoidal terrain correction based on multi-cylindrical equal-area map projection of the reference ellipsoid

An operational algorithm for computation of terrain correction (or local gravity field modeling) based on application of closed-form solution of the Newton integral in terms of Cartesian coordinates in multi-cylindrical equal-area map projection of the reference ellipsoid is presented. Multi-cylindrical equal-area map projection of the reference ellipsoid has been derived and is described in detail for the first time. Ellipsoidal mass elements with various sizes on the surface of the reference ellipsoid are selected and the gravitational potential and vector of gravitational intensity (i.e. gravitational acceleration) of the mass elements are computed via numerical solution of the Newton integral in terms of geodetic coordinates {,,h}. Four base- edge points of the ellipsoidal mass elements are transformed into a multi-cylindrical equal-area map projection surface to build Cartesian mass elements by associating the height of the corresponding ellipsoidal mass elements to the transformed area elements. Using the closed-form solution of the Newton integral in terms of Cartesian coordinates, the gravitational potential and vector of gravitational intensity of the transformed Cartesian mass elements are computed and compared with those of the numerical solution of the Newton integral for the ellipsoidal mass elements in terms of geodetic coordinates. Numerical tests indicate that the difference between the two computations, i.e. numerical solution of the Newton integral for ellipsoidal mass elements in terms of geodetic coordinates and closed-form solution of the Newton integral in terms of Cartesian coordinates, in a multi-cylindrical equal-area map projection, is less than 1.6×10^–8 m²/s² for a mass element with a cross section area of 10×10 m and a height of 10,000 m. For a mass element with a cross section area of 1×1 km and a height of 10,000 m the difference is less than 1.5×10^–4m²/s². Since 1.5× 10^–4 m²/s² is equivalent to 1.5×10^–5m in the vertical direction, it can be concluded that a method for terrain correction (or local gravity field modeling) based on closed-form solution of the Newton integral in terms of Cartesian coordinates of a multi-cylindrical equal-area map projection of the reference ellipsoid has been developed which has the accuracy of terrain correction (or local gravity field modeling) based on the Newton integral in terms of ellipsoidal coordinates.Acknowledgments. This research has been financially supported by the University of Tehran based on grant number 621/4/859. This support is gratefully acknowledged. The authors are also grateful for the comments and corrections made to the initial version of the paper by Dr. S. Petrovic from GFZ Potsdam and the other two anonymous reviewers. Their comments helped to improve the structure of the paper significantly.     

轨道交通中线测设统一数学模型的建立

根据对常用的两种中线测设数学模型的分析，从线路中线曲率半径的特点出发，提出了中线整体积分数学模型，并导出其计算公式。为了实现整体积分模型的计算，对龙贝格积分算法进行了改进，并用于武汉市轨道交通一号线一期工程的中线测设。比较和分析表明，从实用性和精度等方面来说，整体积分模型是目前轨道交通最优的中线测设模型。     

An Alternative Orbit Integration Algorithm for GPS-Based Precise LEO Autonomous Navigation

Single-epoch point positioning with the global positioning system (GPS) is as accurate in low orbit as it is on the ground: typically a three-dimensional rms accuracy of 20 to 30 m as the selective availability turns to zero. This is achieved at any observation epoch without orbit dynamic information. With sophisticated models and filtering techniques onboard the spacecraft, the orbit accuracy of a Low Earth Orbiter (LEO) can be improved to a few meters using the civilian broadcast GPS signals. To achieve this accuracy autonomously in real time, an efficient onboard computing processor is required to carry out the sophisticated orbit integration and filtering process. In this paper, a new orbit integrator is presented that computes the nominal orbit states (the position and velocity) and the state transition equations with numerical methods of integral equation, instead of differential equation usually used for orbit computation. The algorithm is simple, and can be easily embedded in an onboard processor. The numerical results demonstrate that the proposed method of the integral equation provides precise orbit predictions over several orbits. The sequential filter based on the above integrator allows the use of simple orbit state equations to efficiently correct dynamical model errors with precise GPS measurements or improve the orbits using GPS navigaion solutions from the 3D rms accuracy of 26 m to 3.7 m within a few hours of tracking. ? 2001 John Wiley & Sons, Inc.     

Integral formulas for computing the disturbing potential, gravity anomaly and the deflection of the vertical from the second-order radial gradient of the disturbing potential

Integral formulas are derived which can be used to convert the second-order radial gradient of the disturbing potential, as boundary values, into the disturbing potential, gravity anomaly and the deflection of the vertical. The derivations are based on the fundamental differential equation as the boundary condition in Stokes’s boundary-value problem and the modified Poisson integral formula in which the zero and first-degree spherical harmonics are excluded. The rigorous kernel functions, corresponding to the integral operators, are developed by the methods of integration.     

计算子午线弧长的三类算法及其分析比较"

多项式展开算法是计算子午线弧长的传统方法,为了研究利用数值积分算法和常微分方程数值解法进行子午线弧长计算的可行性与可靠性,本文选取大地纬度自0°至90°的3组样本数据(间隔距离分别为1°、1'、1″),分别基于多项式展开数值积分算法和常微分方程数值解法,计算得到各组样本数据的子午线弧长,并通过算法计算结果精度和运算速度两个方面对数值算法的质量进行了评价。计算结果表明:数值积分算法和常微分方程数值解法均可以得到与多项式展开算法精度相同的结果;数值积分算法可通过减小步长以提高计算结果精度,但运算速度急剧降低;3阶、4阶的Runge-Kutta算法不仅运算结果精度高,而且运算速度也比传统算法快3倍多,表明了常微分方程数值解法更适用于子午线弧长的大数据计算。     

A geodetic approach for the recovery of global kinematic plate parameters

Modern techniques of precise geodetic positioning are capable of monitoring global tectonic movements. We can avoid the tremendous effort of observing those point motions at every place on the earth, if we accept the model of rigid tectonic plates, which allows us to extrapolate from discrete point observations to the appertaining plates. The target of describing plate kinematics is the determination of its kinematic parameters, which are the coordinates of the rotation pole and the rotational velocity of each tectonic plate. A mathematical model is presented, which is capable of including geodetic observations (point coordinate shifts, distance changes) as well as geophysical quantities (sea floor spreading rates, earthquake slip vectors). The parameter estimation procedure is derived and demonstrated in simulated examples. Finally a global geodetic network for space techniques is designed, which provides an optimum parameter estimation.     

GPS/GLONASS组合单点定位及精度分析

介绍了基于广播星历的GPS/GLONASS组合导航单点定位的数学模型,分析了组合导航的技术难点。在GPS伪距法单点定位的基础上进行组合导航定位,其中GLONASS卫星坐标运用四阶龙格—库塔(Runge-Kutta)数值积分方法求得,利用一种新的不需要进行轨道拟合的编程方法来进行计算。以IGS跟踪站提供的观测数据为例,分别采用GPS、GLO-NASS和GPS/GLONASS三种方式组合进行伪距法单点定位,同时比较分析了不同权重选择对组合定位精度的影响。