A graph-theoretical algorithm for detecting configuration defects in triangular geodetic networks

Summary The problem to detect configurational defects in geodetic networks is solved by a graph-theoretical algorithm, here applied to triangular geodetic networks and being presented as a computer program in the Appendix. Based on an analysis of the incidence matrix the algorithm detects, for instance, missing vertical directions which cause two types of deficiencies. In case of only vertical direction measurements from one point to another, but no counter vertical direction measurements backwards, the rank deficiency of the first type is identified. Furtheron if there are “bare” points with no vertical direction measurements at all, the rank deficiency of the second type is found. The algorithm has proved a rank deficiency of 4+13=17 in theSW Finland triangular network which before has been found as the surprizing rank defect ofTAGNET 3d-operational adjustment.     

Use of potential coefficient models for geoid undulation determinations using a spherical harmonic representation of the height anomaly/geoid undulation difference

 This paper suggests that potential coefficient models of the Earth's gravitational potential be used to calculate height anomalies which are then reduced to geoid undulations where such quantities are needed for orthometric height determination and vertical datum definition through a potential coefficient realization of the geoid. The process of the conversion of the height anomaly into a geoid undulation is represented by a height anomaly gradient term and the usual N–ζ term that is dependent on elevation and the Bouguer anomaly. Using a degree 360 expansion of 30′ elevations and the OSU91A potential coefficient model, a degree 360 representation of the correction terms was computed. The magnitude of N–ζ reached –3.4 m in the Himalaya Mountains with smaller, but still significant, magnitudes in other mountainous regions. Received: 6 May 1996; Accepted: 30 October 1996     

Precise alpine geoid determination without digital terrain models

The short wavelength geoid undulations, caused by topography, amount to several decimeters in mountainous areas. Up to now these effects are computed by means of digital terrain models in a grid of 100–500m. However, for many countries these data are not yet available or their collection is too expensive. This problem can be overcome by considering the special behaviour of the gravity potential along mountain slopes. It is shown that 90 per cent of the topographic effects are represented by a simple summation formula, based on the average height differences and distances between valleys and ridges along the geoid profiles, δN=[30.H.D.+16.(H−H′).D] in mm/km, (error<10%), whereH, H′, D are estimated in a map to the nearest 0.2km. The formula is valid for asymmetric sides of valleys (H, H′) and can easily be corrected for special shapes. It can be used for topographic refinement of low resolution geoids and for astrogeodetic projects. The “slope method” was tested in two alpine areas (heights up to 3800m, astrogeodetic deflection points every 170km ²) and resulted in a geoid accuracy of ±3cm. In first order triangulation networks (astro points every 1000km ²) or for gravimetric deflections the accuracy is about 10cm per 30km. Since a map scale of 1∶500.000 is sufficient, the method is suitable for developing countries, too.     

A new method for computing the ellipsoidal correction for Stokes's formula

 This paper generalizes the Stokes formula from the spherical boundary surface to the ellipsoidal boundary surface. The resulting solution (ellipsoidal geoidal height), consisting of two parts, i.e. the spherical geoidal height N ₀ evaluated from Stokes's formula and the ellipsoidal correction N ₁, makes the relative geoidal height error decrease from O(e ²) to O(e ⁴), which can be neglected for most practical purposes. The ellipsoidal correction N ₁ is expressed as a sum of an integral about the spherical geoidal height N ₀ and a simple analytical function of N ₀ and the first three geopotential coefficients. The kernel function in the integral has the same degree of singularity at the origin as the original Stokes function. A brief comparison among this and other solutions shows that this solution is more effective than the solutions of Molodensky et al. and Moritz and, when the evaluation of the ellipsoidal correction N ₁ is done in an area where the spherical geoidal height N ₀ has already been evaluated, it is also more effective than the solution of Martinec and Grafarend. Received: 27 January 1999 / Accepted: 4 October 1999     

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     

Accurate formula expressing the difference between the normal gravity and its radial component

A simple formula is presented giving the value of γ−γ _r to better than 0.001 mgal associated with an arbitrary reference ellipsoid, where γ is the normal gravity and γ _r is its radial component. Further simplifications of this formula are possible, depending on the desired accuracy. Since in the actual field g−g_r equals γ−γ _r to a good approximation, this formula makes it possible to work in terms of g_r rather than in terms of the measured quantity g. Such a choice is attractive mainly because the spherical harmonic expansion of g_r is very simple.     

Somigliana–Pizzetti gravity: the international gravity formula accurate to the sub-nanoGal level

 The Somigliana–Pizzetti gravity field (the International gravity formula), namely the gravity field of the level ellipsoid (the International Reference Ellipsoid), is derived to the sub-nanoGal accuracy level in order to fulfil the demands of modern gravimetry (absolute gravimeters, super conducting gravimeters, atomic gravimeters). Equations (53), (54) and (59) summarise Somigliana–Pizzetti gravity Γ(φ,u) as a function of Jacobi spheroidal latitude φ and height u to the order ?(10⁻¹⁰ Gal), and Γ(B,H) as a function of Gauss (surface normal) ellipsoidal latitude B and height H to the order ?(10⁻¹⁰ Gal) as determined by GPS (`global problem solver'). Within the test area of the state of Baden-Württemberg, Somigliana–Pizzetti gravity disturbances of an average of 25.452 mGal were produced. Computer programs for an operational application of the new international gravity formula with (L,B,H) or (λ,φ,u) coordinate inputs to a sub-nanoGal level of accuracy are available on the Internet. Received: 23 June 2000 / Accepted: 2 January 2001     

Polyhedral approximations in physical geodesy

Summary A procedure is derived for the upward continuation of unevenly spaced gravity data. The topographic relief is approximated by a polyhedron with triangular faces and vertices placed at small distances around the surface of a sphere. The usual Fredholm integral equation of the second kind is modified considering the discontinuity of the normal vector. It is solved by successive approximations assuming the unknown function is linear inside each face at every step of the iteration process. An approximate formula to obtain the anomalous potential from the Bouguer anomaly is discussed. The potential of a homogeneous polyhedron is derived and used to compute relief corrections to the geoid undulations. Numerical applications are presented with respect to the Romanian territory.Partially presented at theJoint Symposium of IGC and IGC, Graz, Austria, 11–17 September 1994     

The absolute determination of the gravitational acceleration at Sydney, Australia

An absolute measurement of the gravitational acceleration “g” has been made at the National Standards Laboratory, Chippendale, N.S.W., Australia. The determination was made by studying the free motion of a body projected vertically upwards in a vacuum and the time between its initial and final passages through two horizontal planes of known vertical separation was measured. The measured value ofg at a point 12 metres above the floor in room B. 37 of the National Standards Laboratory is 9.7967134 m/s² The corresponding value at floor level at the BMR gravity station is 9.796717 m/s² Paper presented at the meeting of the International Gravimetric Commission, Paris 7–11 September 1970.     

The optimal Mercator projection and the optimal polycylindric projection of conformal type – case-study Indonesia

As a conformal mapping of the sphere S ² _R or of the ellipsoid of revolution E ² _A , _B the Mercator projection maps the equator equidistantly while the transverse Mercator projection maps the transverse metaequator, the meridian of reference, with equidistance. Accordingly, the Mercator projection is very well suited to geographic regions which extend east-west along the equator; in contrast, the transverse Mercator projection is appropriate for those regions which have a south-north extension. Like the optimal transverse Mercator projection known as the Universal Transverse Mercator Projection (UTM), which maps the meridian of reference Λ₀ with an optimal dilatation factor &ρcirc;=0.999 578 with respect to the World Geodetic Reference System WGS 84 and a strip [Λ₀−Λ _W ,Λ₀ + Λ _E ]×[Φ _S ,Φ _N ]= [−3.5^∘,+3.5^∘]×[−80^∘,+84^∘], we construct an optimal dilatation factor ρ for the optimal Mercator projection, summarized as the Universal Mercator Projection (UM), and an optimal dilatation factor ρ₀ for the optimal polycylindric projection for various strip widths which maps parallel circles Φ₀ equidistantly except for a dilatation factor ρ₀, summarized as the Universal Polycylindric Projection (UPC). It turns out that the optimal dilatation factors are independent of the longitudinal extension of the strip and depend only on the latitude Φ₀ of the parallel circle of reference and the southern and northern extension, namely the latitudes Φ _S and Φ _N , of the strip. For instance, for a strip [Φ _S ,Φ _N ]= [−1.5^∘,+1.5^∘] along the equator Φ₀=0, the optimal Mercator projection with respect to WGS 84 is characterized by an optimal dilatation factor &ρcirc;=0.999 887 (strip width 3^∘). For other strip widths and different choices of the parallel circle of reference Φ₀, precise optimal dilatation factors are given. Finally the UPC for the geographic region of Indonesia is presented as an example. Received: 17 December 1997 / Accepted: 15 August 1997     

Unbiased vs biased estimation of GPS phase ambiguities from dual-frequency code and phase observables

With access to dual-frequency pseudorange and phase Global Positioning System (GPS) data, the wide-lane ambiguity can easily be fixed. Advantage is taken of this information in the linear combination of the above four observables for base ambiguity estimation (i.e. of N ₁ and N ₂). Starting points for our analysis are the Best Linear Unbiased Estimators BLUE₁ and BLUE₂. BLUE₁ is the best one (with minimum mean square error, MSE) if the ionosphere effect is negligible. If this is not the case, BLUE₂ has the smallest variance, but not necessarily the least mean square error. Hence, both estimators may suffer from a non-optimal treatment of the ionosphere bias. BLUE₁ ignores possible ionosphere bias, while BLUE₂ compensates for this bias in a less favourable way by eliminating it at the price of increased noise. As an alternative, linear estimators are derived, which make a compromise between the ionosphere bias and the random observation errors. This leads to the derivation of the Best Linear Estimator (BLE) and the Restricted Best Linear Estimator (RBLE) with minimum MSE. The former is generally not very useful, while the RBLE is recommended for practical use. It is shown that the MSE of the RBLE is limited by the variances of BLUE₁ and BLUE₂, i.e.

The exterior Airy/Heiskanen topographic–isostatic gravity potential, anomaly and the effect of analytical continuation in Stokes' formula

This study deals with the external type of topographic–isostatic potential and gravity anomaly and its vertical derivatives, derived from the Airy/Heiskanen model for isostatic compensation. From the first and the second radial derivatives of the gravity anomaly the effect on the geoid is estimated for the downward continuation of gravity to sea level in the application of Stokes' formula. The major and regional effect is shown to be of order H ³ of the topography, and it is estimated to be negligible at sea level and modest for most mountains, but of the order of several metres for the highest and most extended mountain belts. Another, global, effect is of order H but much less significant Received: 3 October 1997 / Accepted: 30 June 1998     

GPS signal diffraction modelling: the stochastic SIGMA-δ model

The SIGMA-Δ model has been developed for stochastic modelling of global positioning system (GPS) signal diffraction errors in high precision GPS surveys. The basic information used in the SIGMA-Δ model is the measured carrier-to-noise power-density ratio (C/N0). Using the C/N0 data and a template technique, the proper variances are derived for all phase observations. Thus the quality of the measured phase is automatically assessed and if phase observations are suspected to be contaminated by diffraction effects they are weighted down in the least-squares adjustment. The ability of the SIGMA-Δ model to reduce signal diffraction effects is demonstrated on two static GPS surveys as well as on a kinematic high-precision GPS railway survey. In cases of severe signal diffraction the accuracy of the GPS positions is improved by more than 50% compared to standard GPS processing techniques. Received: 27 July 1998 / Accepted: 1 March 1999     

Computer algebra solution of the GPS N-points problem

A computer algebra solution is applied here to develop and evaluate algorithms for solving the basic GPS navigation problem: finding a point position using four or more pseudoranges at one epoch (the GPS N-points problem). Using Mathematica 5.2 software, the GPS N-points problem is solved numerically, symbolically, semi-symbolically, and with Gauss–Jacobi, on a work station. For the case of N > 4, two minimization approaches based on residuals and distance norms are evaluated for the direct numerical solution and their computational duration is compared. For N = 4, it is demonstrated that the symbolic computation is twice as fast as the iterative direct numerical method. For N = 6, the direct numerical solution is twice as fast as the semi-symbolic, with the residual minimization requiring less computation time compared to the minimization of the distance norm. Gauss–Jacobi requires eight times more computation time than the direct numerical solution. It does, however, have the advantage of diagnosing poor satellite geometry and outliers. Besides offering a complete evaluation of these algorithms, we have developed Mathematica 5.2 code (a notebook file) for these algorithms (i.e., Sturmfel’s resultant, Dixon’s resultants, Groebner basis, reduced Groebner basis and Gauss–Jacobi). These are accessible to any geodesist, geophysicist, or geoinformation scientist via the GPS Toolbox () website or the Wolfram Information Center ().

Construction of anisotropic covariance functions using Riesz-representers

A reproducing-kernel Hilbert space (RKHS) of functions harmonic in the set outside a sphere with radius R ₀, having a reproducing kernel K ₀(P,Q) is considered (P, Q, and later P _n being points in the set of harmonicity). The degree variances of this kernel will be denoted σ_{0 n} . The set of Riesz representers associated with the evaluation functionals (or gravity functionals) related to distinct points P _n ,n = 1,…,N, on a two-dimensional surface surrounding the bounding sphere, will be linearly independent. These functions are used to define a new N-dimensional RKHS with kernel (a _n >0)

The CARIB97 high-resolution geoid height model for the Caribbean Sea

A 2×2 arc-minute resolution geoid model, CARIB97, has been computed covering the Caribbean Sea. The geoid undulations refer to the GRS-80 ellipsoid, centered at the ITRF94 (1996.0) origin. The geoid level is defined by adopting the gravity potential on the geoid as W ₀=62 636 856.88 m²/s² and a gravity-mass constant of GM=3.986 004 418×10¹⁴ m³/s². The geoid model was computed by applying high-frequency corrections to the Earth Gravity Model 1996 global geopotential model in a remove-compute-restore procedure. The permanent tide system of CARIB97 is non-tidal. Comparison of CARIB97 geoid heights to 31 GPS/tidal (ITRF94/local) benchmarks shows an average offset (h–H–N) of 51 cm, with an Root Mean Square (RMS) of 62 cm about the average. This represents an improvement over the use of a global geoid model for the region. However, because the measured orthometric heights (H) refer to many differing tidal datums, these comparisons are biased by localized permanent ocean dynamic topography (PODT). Therefore, we interpret the 51 cm as partially an estimate of the average PODT in the vicinity of the 31 island benchmarks. On an island-by-island basis, CARIB97 now offers the ability to analyze local datum problems which were previously unrecognized due to a lack of high-resolution geoid information in the area. Received: 2 January 1998 / Accepted: 18 August 1998     

Evaluation of deflections of the vertical on the sphere and the plane: a comparison of FFT techniques

This paper presents a set of efficient formulas to evaluate the deflections of the vertical on the sphere using gridded data. The Vening-Meinesz formula, the topographic indirect effect on the deflections of the vertical as well as the terrain corrections are expressed as both 2D and 1D convolutions on the sphere, and consequently can be evaluated by the 2D and the 1D fast Fourier transform (FFT). When compared with the results obtained from pointwise integration, the use of the 1D FFT gives identical results, and therefore these results were used as control values in this paper. The use of the spherical 2D FFT improves significantly the computational efficiency with little sacrifice of accuracy (0.6^″ rms difference from the 1D FFT results). The planar 2D FFT, which is as efficient as the spherical 2D FFT, gives worse results (1.2^″ rms difference from the 1D FFT results) because of the extra approximations. Received: 27 February 1996; Accepted: 24 January 1997     

Green's function solution to spherical gradiometric boundary-value problems

 Three independent gradiometric boundary-value problems (BVPs) with three types of gradiometric data, {Γ _rr }, {Γ _{r θ},Γ _{r λ}} and {Γ_θθ−Γ_λλ,Γ_θλ}, prescribed on a sphere are solved to determine the gravitational potential on and outside the sphere. The existence and uniqueness conditions on the solutions are formulated showing that the zero- and the first-degree spherical harmonics are to be removed from {Γ _{r θ},Γ _{r λ}} and {Γ_θθ−Γ_λλ,Γ_θλ}, respectively. The solutions to the gradiometric BVPs are presented in terms of Green's functions, which are expressed in both spectral and closed spatial forms. The logarithmic singularity of the Green's function at the point ψ=0 is investigated for the component Γ _rr . The other two Green's functions are finite at this point. Comparisons to the paper by van Gelderen and Rummel [Journal of Geodesy (2001) 75: 1–11] show that the presented solution refines the former solution. Received: 3 October 2001 / Accepted: 4 October 2002     

On non-linear low–low SST observation equations for the determination of the geopotential based on an analytical solution

In satellite data analysis, one big advantage of analytical orbit integration, which cannot be overestimated, is missed in the numerical integration approach: spectral analysis or the lumped coefficient concept may be used not only to design efficient algorithms but overall for much better insight into the force-field determination problem. The lumped coefficient concept, considered from a practical point of view, consists of the separation of the observation equation matrix A=BT into the product of two matrices. The matrix T is a very sparse matrix separating into small block-diagonal matrices connecting the harmonic coefficients with the lumped coefficients. The lumped coefficients are nothing other than the amplitudes of trigonometric functions depending on three angular orbital variables; therefore, the matrix N=B ^T B will become for a sufficient length of a data set a diagonal dominant matrix, in the case of an unlimited data string length a strictly diagonal one. Using an analytical solution of high order, the non-linear observation equations for low–low SST range data can be transformed into a form to allow the application of the lumped concept. They are presented here for a second-order solution together with an outline of how to proceed with data analysis in the spectral domain in such a case. The dynamic model presented here provides not only a practical algorithm for the parameter determination but also a simple method for an investigation of some fundamental questions, such as the determination of the range of the subset of geopotential coefficients which can be properly determined by means of SST techniques or the definition of an optimal orbital configuration for particular SST missions. Numerical results have already been obtained and will be published elsewhere. Received: 15 January 1999 / Accepted: 30 November 1999     

Applications of Lodrigues Matrix in 3D Coordinate Transformation

Three transformation models (Bursa-Wolf, Molodensky, and WTUSM) are generally used between two data systems transformation. The linear models are used when the rotation angles are small; however, when the rotation angles get bigger, model errors will be produced. In this paper, we present a method with three main terms: ① the traditional rotation angles θ , φ ,ψ are substituted with a , b, c which are three re-spective values in the anti-symmetrical or Lodrigues matrix; ② directly and accurately calculating the formula of seven parameters in any value of rotation angles; and ③ a corresponding adjustment model is established. This method does not use the triangle function. Instead it uses addition, subtraction, multiplication and division, and the complexity of the equation is reduced, making the calculation easy and quick.