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.     

Intuitive derivation of loop inverses and array algebra

Array algebra forms the general base of fast transforms and multilinear algebra making rigorous solutions of a large number (millions) of parameters computationally feasible. Loop inverses are operators solving the problem of general matrix inverses. Their derivation starts from the inconsistent linear equations by a parameter exchangeX→L ₀, where X is a set of unknown observables,A ₀ forming a basis of the so called “problem space”. The resulting full rank design matrix of parameters L₀ and its ℓ-inverse reveal properties speeding the computational least squares solution expressed in observed values . The loop inverses are found by the back substitution expressing ∧X in terms ofL through . Ifp=rank (A) ≤n, this chain operator creates the pseudoinverseA ⁺. The idea of loop inverses and array algebra started in the late60's from the further specialized case,p=n=rank (A), where the loop inverse A ₀ ⁻¹ (AA ₀ ⁻¹ )^ℓ reduces into the ℓ-inverse A^ℓ=(A^TA)⁻¹A^T. The physical interpretation of the design matrixA A ₀ ⁻¹ as an interpolator, associated with the parametersL ₀, and the consideration of its multidimensional version has resulted in extended rules of matrix and tensor calculus and mathematical statistics called array algebra.     

A test of significance for the Helmert-Kubik-Problem of Weight-determination

As it has been shown by Kubik it is possible to get an estimate, , of the reciprocal of the weight-matrix in an adjustment problem. If we want to see whether this new estimate differssignificantly from our a priori valueQ ₀ it is necessary to know the distribution function of the elements , the ’s being the elements of . This distribution is found in the present article and it is shown that it is not identical with any of the distributions well known from statistical textbooks. Furthermore a way of computing this new distribution is presented. Finally the connection with the chi-square distribution is explored and it is proved that the chi-square-distribution may be used as an approximation for a large number of over-determinations.     

The ionospheric eclipse factor method (IEFM) and its application to determining the ionospheric delay for GPS

A new method for modeling the ionospheric delay using global positioning system (GPS) data is proposed, called the ionospheric eclipse factor method (IEFM). It is based on establishing a concept referred to as the ionospheric eclipse factor (IEF) λ of the ionospheric pierce point (IPP) and the IEF’s influence factor (IFF) . The IEF can be used to make a relatively precise distinction between ionospheric daytime and nighttime, whereas the IFF is advantageous for describing the IEF’s variations with day, month, season and year, associated with seasonal variations of total electron content (TEC) of the ionosphere. By combining λ and with the local time t of IPP, the IEFM has the ability to precisely distinguish between ionospheric daytime and nighttime, as well as efficiently combine them during different seasons or months over a year at the IPP. The IEFM-based ionospheric delay estimates are validated by combining an absolute positioning mode with several ionospheric delay correction models or algorithms, using GPS data at an international Global Navigation Satellite System (GNSS) service (IGS) station (WTZR). Our results indicate that the IEFM may further improve ionospheric delay modeling using GPS data.     

An efficient short-ARC orbit computation

Short-arc orbit computations by numerical or analytical integration of equations of motion traditionally utilized in geodetic and geodynamic satellite positioning are relatively involved and computationally expensive. However, short-arc orbits can be evaluated more efficiently by means of least squares polynomial approximations. Such orbit computations do not significantly increase the computation time when compared to widely used semi-short-arc techniques which utilize externally generated orbits. The sufficiently high-degree polynomial approximation of the second time derivatives , and evaluated from a gravitational potential model at regular (two-minute) intervals and everaged initial conditions (position and velocity vectors at the beginning, the middle and the end of a pass) reproduces the U.S. Defense Mapping Agency precise ephemeris of the Navy Navigation Satellites (NNSS) to about 5 cm RMS in each coordinate. To achieve this level of orbit shape resolution for NNSS satellites, the gravitational potential model should not be truncated at less than degree and order 10. Contribution of the Earth Physics Branch No. 1034.     

Molodensky’s problem in gravity space: A review of the first results

In the last year a new formulation of Molodensky's problem has been given, in which the gravity vector has been considered as the independent variable of the problem, while the position vector is the dependent. This new approach has the great advantage to transform the problem of Molodensky which is of free boundary type, into a fixed boundary problem for a non linear differential equations. In this paper the first results of the study of the new approach are summarized, without going into many mathematical details. The problem of Molodensky for the rotating earth is also discussed.     

Elliptical harmonic series and the original Stokes problem with the boundary of the reference ellipsoid

Since the earth is closer to a revolving ellipsoid than a sphere, it is very important to study directly the original model of the Stokes' BVP on the reference ellipsoid, where denotes the reference ellipsoid, is the Somigliana normal gravity, andh is the outer normal direction of. This paper deals with: 1) simplification of the above BVP under preserving accuracy to , 2) derivation of computational formula of the elliptical harmonic series, 3) solving the BVP by the elliptical harmonic series, and 4) providing a principle for finding the elliptical harmonic model of the earth's gravity field from the spherical harmonic coefficients ofg. All results given in the paper have the same accuracy as the original BVP, that is, the accuracy of the BVP is theoretically preserved in each derivation step.     

A model comparison in least squares collocation

The well known least squares collocation model (I) $$\ell = Ax + \left[ {\begin{array}{*{20}c} O \\ I \\ \end{array} } \right]^T \left[ {\begin{array}{*{20}c} s \\ {s' + n} \\ \end{array} } \right]$$ is compared with the model (II) $$\ell = Ax + \left[ {\begin{array}{*{20}c} R \\ I \\ \end{array} } \right]^T \left[ {\begin{array}{*{20}c} s \\ n \\ \end{array} } \right]$$ The basic differences of these two models in the framework of physical geodesy are pointed out by analyzing the validity of the equation $$s' = Rs$$ that transforms one model into the other, for different cases. For clarification purposes least squares filtering, prediction and collocation are discussed separately. In filtering problems the coefficient matrix R becomes the unit matrix and by this the two models become identical. For prediction and collocation problems the relation s′=Rs is only fulfilled in the global limit where s becomes either a continuous function on the earth or an intinite set of spherical harmonic coefficients. Applying Model (II), we see that for any finite dimension of s the operator equations of physical geodesy are approximated by a finite matrix relation whereas in Model (I) the operator equations are applied in their correct form on a continuous, approximate function \(\tilde s\) .     

Rectangular rotation of spherical harmonic expansion of arbitrary high degree and order

In order to move the polar singularity of arbitrary spherical harmonic expansion to a point on the equator, we rotate the expansion around the y-axis by \(90^{\circ }\) such that the x-axis becomes a new pole. The expansion coefficients are transformed by multiplying a special value of Wigner D-matrix and a normalization factor. The transformation matrix is unchanged whether the coefficients are \(4 \pi \) fully normalized or Schmidt quasi-normalized. The matrix is recursively computed by the so-called X-number formulation (Fukushima in J Geodesy 86: 271–285, 2012a). As an example, we obtained \(2190\times 2190\) coefficients of the rectangular rotated spherical harmonic expansion of EGM2008. A proper combination of the original and the rotated expansions will be useful in (i) integrating the polar orbits of artificial satellites precisely and (ii) synthesizing/analyzing the gravitational/geomagnetic potentials and their derivatives accurately in the high latitude regions including the arctic and antarctic area.     

Horizontal strain and precision of EDM from repeated surveys of the Precise Geodetic Net of Japan

A simple statistical approach has been applied to the repeated electro-optical distance measurements (EDM) of 1,358 lines in the Tohoku district of Japan to obtain knowledge about the precision of EDM and the possible accumulation of strain. The average time interval between measurements is about seven or eight years. It is shown that the whole data of the difference between distance measurements repeated over a given lineD are interpreted in terms of EDM errors comprising distance proportional systematic errors and standard errors expressed by the usual form . The rate of horizontal deformation must therefore be much smaller than the strain rates of about 0.7 0.8 ppm over 7 to 8 years which have been hitherto expected.     

On the computation and approximation of ultra-high-degree spherical harmonic series

Spherical harmonic series, commonly used to represent the Earth’s gravitational field, are now routinely expanded to ultra-high degree (> 2,000), where the computations of the associated Legendre functions exhibit extremely large ranges (thousands of orders) of magnitudes with varying latitude. We show that in the degree-and-order domain, (ℓ,m), of these functions (with full ortho-normalization), their rather stable oscillatory behavior is distinctly separated from a region of very strong attenuation by a simple linear relationship: , where θ is the polar angle. Derivatives and integrals of associated Legendre functions have these same characteristics. This leads to an operational approach to the computation of spherical harmonic series, including derivatives and integrals of such series, that neglects the numerically insignificant functions on the basis of the above empirical relationship and obviates any concern about their broad range of magnitudes in the recursion formulas that are used to compute them. Tests with a simulated gravitational field show that the errors in so doing can be made less than the data noise at all latitudes and up to expansion degree of at least 10,800. Neglecting numerically insignificant terms in the spherical harmonic series also offers a computational savings of at least one third.     

Graph theory for analyzing pair-wise data: application to geophysical model parameters estimated from interferometric synthetic aperture radar data at Okmok volcano,Alaska

Graph theory is useful for analyzing time-dependent model parameters estimated from interferometric synthetic aperture radar (InSAR) data in the temporal domain. Plotting acquisition dates (epochs) as vertices and pair-wise interferometric combinations as edges defines an incidence graph. The edge-vertex incidence matrix and the normalized edge Laplacian matrix are factors in the covariance matrix for the pair-wise data. Using empirical measures of residual scatter in the pair-wise observations, we estimate the relative variance at each epoch by inverting the covariance of the pair-wise data. We evaluate the rank deficiency of the corresponding least-squares problem via the edge-vertex incidence matrix. We implement our method in a MATLAB software package called GraphTreeTA available on GitHub (https://github.com/feigl/gipht). We apply temporal adjustment to the data set described in Lu et al. (Geophys Res Solid Earth 110, 2005) at Okmok volcano, Alaska, which erupted most recently in 1997 and 2008. The data set contains 44 differential volumetric changes and uncertainties estimated from interferograms between 1997 and 2004. Estimates show that approximately half of the magma volume lost during the 1997 eruption was recovered by the summer of 2003. Between June 2002 and September 2003, the estimated rate of volumetric increase is \((6.2 \, \pm \, 0.6) \times 10^6~\mathrm{m}^3/\mathrm{year} \). Our preferred model provides a reasonable fit that is compatible with viscoelastic relaxation in the five years following the 1997 eruption. Although we demonstrate the approach using volumetric rates of change, our formulation in terms of incidence graphs applies to any quantity derived from pair-wise differences, such as range change, range gradient, or atmospheric delay.     

Fast error analysis of continuous GPS observations

It has been generally accepted that the noise in continuous GPS observations can be well described by a power-law plus white noise model. Using maximum likelihood estimation (MLE) the numerical values of the noise model can be estimated. Current methods require calculating the data covariance matrix and inverting it, which is a significant computational burden. Analysing 10 years of daily GPS solutions of a single station can take around 2 h on a regular computer such as a PC with an AMD Athlon^TM 64 X2 dual core processor. When one analyses large networks with hundreds of stations or when one analyses hourly instead of daily solutions, the long computation times becomes a problem. In case the signal only contains power-law noise, the MLE computations can be simplified to a process where N is the number of observations. For the general case of power-law plus white noise, we present a modification of the MLE equations that allows us to reduce the number of computations within the algorithm from a cubic to a quadratic function of the number of observations when there are no data gaps. For time-series of three and eight years, this means in practise a reduction factor of around 35 and 84 in computation time without loss of accuracy. In addition, this modification removes the implicit assumption that there is no environment noise before the first observation. Finally, we present an analytical expression for the uncertainty of the estimated trend if the data only contains power-law noise. Electronic supplementary material The online version of this article (doi: ) contains supplementary material, which is available to authorized users.     

Ambiguity resolution strategies using the results of the International GPS Geodynamics Service (IGS)

Resolving the initial phase ambiguities of GPS carrier phase observations was always considered an important aspect of GPS processing techniques. Resolution of the so-called wide-lane ambiguities using a special linear combination of theL ₁ andL ₂ carrier and code observations has become standard. New aspects have to be considered today: (1) Soon AS, the so-called Anti-Spoofing, will be turned on for all Block II spacecrafts. This means that precise code observations will be no longer available, which in turn means that the mentioned approach to resolve the wide-lane ambiguities will fail. (2) Most encouraging is the establishment of the new International GPS Geodynamics Service (IGS), from where high quality orbits, earth rotation parameters, and eventually also ionospheric models will be available. We are reviewing the ambiguity resolution problem under these new aspects: We look for methods to resolve the initial phase ambiguities without using code observations but using high quality orbits and ionospheric models from IGS, and we study the resolution of the narrow-lane ambiguities (after wide-lane ambiguity resolution) using IGS orbits.     

Estimating crustal deformations from a combination of baseline measurements and geophysical models

The estimation of crustal deformations from repeated baseline measurements is a singular problem in the absence of prior information. One often applied solution is a free adjustment in which the singular normal matrix is augmented with a set of inner constraints. These constraints impose no net translation nor rotation for the estimated deformations X which may not be physically meaningful for a particular problem. The introduction of an available geophysical model from which an expected deformation vector \(\bar X\) and its covariance matrix \(\sum _{\bar X} \) can be computed will direct X to a physically more meaningful solution. Three possible estimators are investigated for estimating deformations from a combination of baseline measurements and geophysical models.     

Fast error analysis of continuous GNSS observations with missing data

One of the most widely used method for the time-series analysis of continuous Global Navigation Satellite System (GNSS) observations is Maximum Likelihood Estimation (MLE) which in most implementations requires $\mathcal{O }(n^3)$ operations for $n$ observations. Previous research by the authors has shown that this amount of operations can be reduced to $\mathcal{O }(n^2)$ for observations without missing data. In the current research we present a reformulation of the equations that preserves this low amount of operations, even in the common situation of having some missing data.Our reformulation assumes that the noise is stationary to ensure a Toeplitz covariance matrix. However, most GNSS time-series exhibit power-law noise which is weakly non-stationary. To overcome this problem, we present a Toeplitz covariance matrix that provides an approximation for power-law noise that is accurate for most GNSS time-series.Numerical results are given for a set of synthetic data and a set of International GNSS Service (IGS) stations, demonstrating a reduction in computation time of a factor of 10–100 compared to the standard MLE method, depending on the length of the time-series and the amount of missing data.     

Approximating the Bayesian estimate of the standard deviation in a linear model

The Bayesian estimates _b of the standard deviation σ in a linear model—as needed for the evaluation of reliability—is well known to be proportional to the square root of the Bayesian estimate (s ²) _b of the variance component σ² by a proportionality factor involving the ratio of Gamma functions. However, in analogy to the case of the respective unbiased estimates, the troublesome exact computation ofa _b may be avoided by a simple approximation which turns out to be good enough for most applications even if the degree of freedom ν is rather small. Paper presented to the Int. Conf. on “Practical Bayesian Statistics”, Cambridge (U.K.), 8.–11. July 1986.