On the foundation of collocation in physical geodesy

The boundary value problem in physical geodesy is nowadays mostly presented with the use of an advanced stochastic model by Krarup-Moritz. This model includes a primary Gauss-Markov model and an adjoining Wiener-Hopf model. Degenerations of the Wiener-Hopf section are found in thesingular auto-covariance matrix of the residuals. The non-singular inverse of the auto-covariance matrix of the signal is proved to be a generalized inverse of the singular auto-covariance matrix of the residuals. The joint model is given a non-stochastic evaluation for a case with spherical external surface (using a non-singular inverse). These findings will not prevent a successful application of the model, which has important merits, specially when using suitablea priori values for the stochastic parameters in the covariance functions. A method for quadratic unbiased estimation ofa priori variances is presented in an introductory section. It is meant to be of value when using a solution of the boundary value problem with the collocation technique based on the classical Gauss-Markov solution. (Bjerhammar (1963).)     

不同扰动位泛函间的积分变换广义核函数

基于物理大地测量边值问题的解,利用一阶边界算子定义,推导重力异常Δg、单层密度μ、大地水准面高N,垂线偏差ε、扰动重力δg等扰动场元的解。利用球谐函数的正交特性,通过对核函数的算子运算,可以得到上述扰动场元的有关逆变换公式。相对经典物理大地测量公式应用的边界面条件,笔者将含有因子r的对应扰动场元反演关系的公式称为广义积分公式。针对常用的重力异常Δg、大地水准面高N,垂线偏差ε、扰动重力δg计算,重点分析它们之间的变换关系,给出利用某个选定扰动场元计算其他扰动场元的广义积分公式。同时,通过对积分边界面的讨论,分析经典公式与广义积分公式的差异和联系。最后,给出所有外部扰动场元与核函数映射的关系表。     

Helmert扰动位及其积分核函数的椭球实用公式

借助以地心参考椭球面为边界面的第二大地边值问题的理论,基于Helmert空间的Neumann边值条件,给定Helmert扰动位的椭球解表达式,并详细推导第二类勒让德函数及其导数的递推关系、Helmert扰动位函数的椭球积分解以及类椭球Hotine积分核函数的实用计算公式,便于后续椭球域第二大地边值问题的实际研究。     

On a four-dimensional integrated geodesy

In contrast to continuous global considerations of time dependent boundary value problems an attempt is made to define4D-linear observation equations in the framework of integrated geodesy for discrete, more or less regional and local applications (deformation analysis) where time variations in position and in the gravity field have to be considered. The derivation is a strict analogue and extension of the3D integrated approach. In addition the construction of time dependent covariance functions is discussed, which are necessary to solve for unknown displacements and changes in the gravity potential in the generalized least squares collocation model.     

An algorithm for the inverse of a multivariate homogeneous polynomial of degreen

Summary Various geodetic problems (the free nonlinear geodetic boundary value problem, the computation of Gauß-Krüger coordinates or UTM coordinates, the problem of nonlinear regression) demand theinversion of an univariate, bivariate, trivariate, in generalmultivariate homogeneous polynomial of degree n. The new algorithm which is oriented towardsSymbolic Computer Manipulation is based upon the algebraic power base computation with respect toKronecker-Zehfu product structure leading to the solution of a system oftriangular matrix equations: Only the first row of the inverse triangular matrix has to be computed. TheSymbolic Computer Manipulation program of the GKS algorithm is available from the authors.     

GPS在物理大地测量中的应用及GPS边值问题

针对GPS技术逐步用于物理大地测量的现状，从物理大地测量的基本原理及方法入手，分析了GPS对传统物理大地测量理论及方法所产生的影响、GPS在解决物理大地测量中的一些难题所发挥的作用以及对物理大地测量的一些功能所造成的变化。重点讨论了物理大地测量边值问题的基本属性的改变、第二边值问题的应用、高程系统的确定、GPS边值问题的定义、特征及求解。最后，对GPS与物理大地测量进一步结合所面临的有关问题进行了简要论述。     

The free versus fixed geodetic boundary value problem for different combinations of geodetic observables

Various formulations of the geodetic fixed and free boundary value problem are presented, depending upon the type of boundary data. For the free problem, boundary data of type astronomical latitude, astronomical longitude and a pair of the triplet potential, zero and first-order vertical gradient of gravity are presupposed. For the fixed problem, either the potential or gravity or the vertical gradient of gravity is assumed to be given on the boundary. The potential and its derivatives on the boundary surface are linearized with respect to a reference potential and a reference surface by Taylor expansion. The Eulerian and Lagrangean concepts of a perturbation theory of the nonlinear geodetic boundary value problem are reviewed. Finally the boundary value problems are solved by Hilbert space techniques leading to new generalized Stokes and Hotine functions. Reduced Stokes and Hotine functions are recommended for numerical reasons. For the case of a boundary surface representing the topography a base representation of the solution is achieved by solving an infinite dimensional system of equations. This system of equations is obtained by means of the product-sum-formula for scalar surface spherical harmonics with Wigner 3j-coefficients.     

The free versus fixed geodetic boundary value problem for different combinations of geodetic observables

Various formulations of the geodetic fixed and free boundary value problem are presented, depending upon the type of boundary data. For the free problem, boundary data of type astronomical latitude, astronomical longitude and a pair of the triplet potential, zero and first-order vertical gradient of gravity are presupposed. For the fixed problem, either the potential or gravity or the vertical gradient of gravity is assumed to be given on the boundary. The potential and its derivatives on the boundary surface are linearized with respect to a reference potential and a reference surface by Taylor expansion. The Eulerian and Lagrangean concepts of a perturbation theory of the nonlinear geodetic boundary value problem are reviewed. Finally the boundary value problems are solved by Hilbert space techniques leading to new generalized Stokes and Hotine functions. Reduced Stokes and Hotine functions are recommended for numerical reasons. For the case of a boundary surface representing the topography a base representation of the solution is achieved by solving an infinite dimensional system of equations. This system of equations is obtained by means of the product-sum-formula for scalar surface spherical harmonics with Wigner 3j-coefficients.     

The nullfield method for the ellipsoidal Stokes problem

The ellipsoidal Stokes problem is one of the basic boundary-value problems for the Laplace equation which arises in physical geodesy. Up to now, geodecists have treated this and related problems with high-order series expansions of spherical and spheroidal (ellipsoidal) harmonics. In view of increasing computational power and modern numerical techniques, boundary element methods have become more and more popular in the last decade. This article demonstrates and investigates the nullfield method for a class of Robin boundary-value problems. The ellipsoidal Stokes problem belongs to this class. An integral equation formulation is achieved, and existence and uniqueness conditions are attained in view of the Fredholm alternative. Explicit expressions for the eigenvalues and eigenfunctions for the boundary integral operator are provided. Received: 22 October 1996 / Accepted: 4 August 1997     

Linear homotopy solution of nonlinear systems of equations in geodesy

A fundamental task in geodesy is solving systems of equations. Many geodetic problems are represented as systems of multivariate polynomials. A common problem in solving such systems is improper initial starting values for iterative methods, leading to convergence to solutions with no physical meaning, or to convergence that requires global methods. Though symbolic methods such as Groebner bases or resultants have been shown to be very efficient, i.e., providing solutions for determined systems such as 3-point problem of 3D affine transformation, the symbolic algebra can be very time consuming, even with special Computer Algebra Systems (CAS). This study proposes the Linear Homotopy method that can be implemented easily in high-level computer languages like C++ and Fortran that are faster than CAS by at least two orders of magnitude. Using Mathematica, the power of Homotopy is demonstrated in solving three nonlinear geodetic problems: resection, GPS positioning, and affine transformation. The method enlarging the domain of convergence is found to be efficient, less sensitive to rounding of numbers, and has lower complexity compared to other local methods like Newton–Raphson.     

Singularity free formulations of the geodetic boundary value problem in gravity-space

The idea of transforming the geodetic boundary value problem into a boundary value problem with a fixed boundary dates back to the 1970s of the last century. This transformation was found by F. Sanso and was named as gravity-space transformation. Unfortunately, the advantage of having a fixed boundary for the transformed problem was counterbalanced by the theoretical as well as practical disadvantage of a singularity at the origin. In the present paper two more versions of a gravity-space transformation are investigated, where none of them has a singularity. In both cases the transformed differential equations are nonlinear. Therefore, a special emphasis is laid on the linearized problems and their relationships to the simple Hotine-problem and to the symmetries between both formulations. Finally, in numerical simulation study the accuracy of the solutions of both linearized problems is studied and factors limiting this accuracy are identified.     

超定边值问题的差分法

物理大地测量面临着越来越多的数据:高程异常、垂线偏差、重力异常、重力梯度等,因此出现了超定边值问题,本文采用求解偏微分方程最简单而又最常用的差分法,对这一问题进行了初步的研究。     

Numerical solution of the linearized fixed gravimetric boundary-value problem

The fixed gravimetric boundary-value problem (FGBVP) represents an exterior oblique derivative problem for the Laplace equation. Terrestrial gravimetric measurements located by precise satellite positioning yield oblique derivative boundary conditions in the form of surface gravity disturbances. In this paper, we discuss the boundary element method (BEM) applied to the linearized FGBVP. In spite of previous BEM approaches in geodesy, we use the so-called direct BEM formulation, where a weak formulation is derived through the method of weighted residuals. The collocation technique with linear basis functions is applied for deriving the linear system of equations from the arising boundary integral equations. The nonstationary iterative biconjugate gradient stabilized method is used to solve the large-scale linear system of equations. The standard MPI (message passing interface) subroutines are implemented in order to perform parallel computations. The proposed approach gives a numerical solution at collocation points directly on the Earth’s surface (on a fixed boundary). Numerical experiments deal with (i) global gravity field modelling using synthetic data (surface gravity disturbances generated from a global geopotential model (GGM)) (ii) local gravity field modelling in Slovakia using observed gravity data. In order to extend computations, the memory requirements are reduced using elimination of the far-zone effects by incorporating GGM or a coarse global numerical solution obtained by BEM. Statistical characteristics of residuals between numerical solutions and GGM confirm the reliability of the approach and indicate accuracy of numerical solutions for the global models. A local refinement in Slovakia results in a local (national) quasigeoid model, which when compared with GPS-levelling data, does not make a large improvement on existing remove-restore-based models.     

Fast numerical solution of the linearized Molodensky problem

 When standard boundary element methods (BEM) are used in order to solve the linearized vector Molodensky problem we are confronted with two problems: (1) the absence of O(|x|⁻²) terms in the decay condition is not taken into account, since the single-layer ansatz, which is commonly used as representation of the disturbing potential, is of the order O(|x|⁻¹) as x→∞. This implies that the standard theory of Galerkin BEM is not applicable since the injectivity of the integral operator fails; (2) the N×N stiffness matrix is dense, with N typically of the order 10⁵. Without fast algorithms, which provide suitable approximations to the stiffness matrix by a sparse one with O(N(logN) ^s ), s≥0, non-zero elements, high-resolution global gravity field recovery is not feasible. Solutions to both problems are proposed. (1) A proper variational formulation taking the decay condition into account is based on some closed subspace of co-dimension 3 of the space of square integrable functions on the boundary surface. Instead of imposing the constraints directly on the boundary element trial space, they are incorporated into a variational formulation by penalization with a Lagrange multiplier. The conforming discretization yields an augmented linear system of equations of dimension N+3×N+3. The penalty term guarantees the well-posedness of the problem, and gives precise information about the incompatibility of the data. (2) Since the upper left submatrix of dimension N×N of the augmented system is the stiffness matrix of the standard BEM, the approach allows all techniques to be used to generate sparse approximations to the stiffness matrix, such as wavelets, fast multipole methods, panel clustering etc., without any modification. A combination of panel clustering and fast multipole method is used in order to solve the augmented linear system of equations in O(N) operations. The method is based on an approximation of the kernel function of the integral operator by a degenerate kernel in the far field, which is provided by a multipole expansion of the kernel function. Numerical experiments show that the fast algorithm is superior to the standard BEM algorithm in terms of CPU time by about three orders of magnitude for N=65 538 unknowns. Similar holds for the storage requirements. About 30 iterations are necessary in order to solve the linear system of equations using the generalized minimum residual method (GMRES). The number of iterations is almost independent of the number of unknowns, which indicates good conditioning of the system matrix. Received: 16 October 1999 / Accepted: 28 February 2001     

The application of commutative algebra to geodesy: two examples

Algebra, in particular commutative algebra, is applied here to provide a general unified solution to nonlinear systems of equations encountered in geodesy. Starting with the “Abelian group”, the “polynomial ring” is defined and used to form generators of ideals. By applying Buchberger or polynomial resultant algorithms, these generators are reduced to simple structures often comprising a univariate polynomial in one of the unknowns. The advantage of the proposed unified approach is that it provides exact solutions to geodetic nonlinear systems of equations without the traditional requirements of linearization, iterations or approximate starting values. The commutative algebraic approach therefore alleviates the need for isolated exact solutions to various geodetic nonlinear systems of equations. The procedure is applied to GPS meteorology to compute refraction angles, and Helmert’s one-to-one mapping of topographical points onto the reference ellipsoid.     

FFT实现物理大地测量边值问题类型转换

本文讨论了用快速付里叶变换(FFT)由大地水准面高度数据反算重力异常数据,由此提出了用FFT实现物理大地测量过值问题类型转换这一思想,并与正交基函数法进行了比较,从而得出运用FFT这一方法能达到简化边值问题本身和运算步骤以及节省计算时间的目的。     

The convergence of harmonic reduction to an internal sphere

The boundary value problem of physical geodesy has been solved with the use of a harmonic reduction down to an internal sphere using a discrete procedure. (For gravity cf. Bjerhammar 1964 and for the potential cf. Bjerhammar 1968). This was a finite-dimensional approach mostly with one-to-one correspondence between observations and unknowns on the sphere. Earlier studies were made with the use of surface elements (on the sphere) with constantgravity. Integration over the surface elements was replaced by a discrete approach with the use of the distance to a point in the centre of the surface element. See Bjerhammar (1968) and (1969). This approach was later presented as a “reflexive prediction” technique for a weakly stationary stochastic process. Bjerhammar (1974, 1976). Krarup (1969) minimized the L²-norm of the potential on the internal sphere. It will here be proved that the two solutions are identical for a proper choice of the radii of the internal spheres. The proof is given for a spherical earth with selected choice of “carrier points”. The convergence problem is discussed. The L²-norm solution is found convergent for the fully harmonic case. Uniform convergence is obtained in the non-harmonic case with the use of the original procedure applied in accordance with the theorems of Keldych-Lavrentieff and Yamabe.     

The spherical horizontal and spherical vertical boundary value problem – vertical deflections and geoidal undulations – the completed Meissl diagram

 In a comparison of the solution of the spherical horizontal and vertical boundary value problems of physical geodesy it is aimed to construct downward continuation operators for vertical deflections (surface gradient of the incremental gravitational potential) and for gravity disturbances (vertical derivative of the incremental gravitational potential) from points on the Earth's topographic surface or of the three-dimensional (3-D) Euclidean space nearby down to the international reference sphere (IRS). First the horizontal and vertical components of the gravity vector, namely spherical vertical deflections and spherical gravity disturbances, are set up. Second, the horizontal and vertical boundary value problem in spherical gravity and geometry space is considered. The incremental gravity vector is represented in terms of vector spherical harmonics. The solution of horizontal spherical boundary problem in terms of the horizontal vector-valued Green function converts vertical deflections given on the IRS to the incremental gravitational potential external in the 3-D Euclidean space. The horizontal Green functions specialized to evaluation and source points on the IRS coincide with the Stokes kernel for vertical deflections. Third, the vertical spherical boundary value problem is solved in terms of the vertical scalar-valued Green function. Fourth, the operators for upward continuation of vertical deflections given on the IRS to vertical deflections in its external 3-D Euclidean space are constructed. Fifth, the operators for upward continuation of incremental gravity given on the IRS to incremental gravity to the external 3-D Euclidean space are generated. Finally, Meissl-type diagrams for upward continuation and regularized downward continuation of horizontal and vertical gravity data, namely vertical deflection and incremental gravity, are produced. Received: 10 May 2000 / Accepted: 26 February 2001     

Nonlinear least-squares method via an isomorphic geometrical setup

The resolution of a nonlinear parametric adjustment model is addressed through an isomorphic geometrical setup with tensor structure and notation, represented by a u-dimensional “model surface” embedded in a flat n-dimensional “observational space”. Then observations correspond to the observational-space coordinates of the pointQ, theu initial parameters correspond to the model-surface coordinates of the “initial” pointP, and theu adjusted parameters correspond to the model-surface coordinates of the “least-squares” point . The least-squares criterion results in a minimum-distance property implying that the vector Q must be orthogonal to the model surface. The geometrical setup leads to the solution of modified normal equations, characterized by a positive-definite matrix. The latter contains second-order and, optionally, thirdorder partial derivatives of the observables with respect to the parameters. This approach significantly shortens the convergence process as compared to the standard (linearized) method.