Search

SUMMARY Although slip inversion has been almost always reported to be inherently non-unique, we prove that Fredholm integral equations of the first kind for slip inversion are mathematically of a unique solution, theoretically assuring that earthquake rupture can be properly reconstructed if there exist a sufficiently large number of measurements. The statement about non-uniqueness of slip inversion can be misleading and misunderstanding, since it is theoretically not true but simply caused by a practical issue of lack of measurements. We propose an inequality-constrained regularized inversion of slip distributions on multiple faults. The method implements a physically more general inequality constraint to accommodate more complex dislocation models and allows reconstructing a complex earthquake mechanism on multiple faults from measurements. The corresponding inequality constraints are a natural extension of positivity constraints proposed by Olson & Apsel and Hartzell & Heaton, or the equivalent inequality constraints that only allow slips to take place between −45° and 45° around the main rupture direction. The regularization parameter is chosen by minimizing the mean squared errors of inverted slip solutions. The proposed method is applied to the 2016 Kumamoto Mw 7.0 earthquake with GNSS measurements. Our slip inversion results show that the 2016 Kumamoto earthquake is only of magnitude Mw 6.7 and ruptures shallowly up to 10 km under the surface. More precisely, about 94 and 88 per cent of the energy released from Hinagu and Futagawa faults take place up to this depth, with the maximum slips of 4.81 and 7.89 m on each fault, respectively.

Accelerometers are supposed to be physical instruments for measuring the acceleration of a moving object. Although huge technological advance is made in hardware of accelerometers over more than one century, accelerometers have been persistently designed and fabricated mechanically under the framework of forward problems of damped mass-spring systems. However, accelerometers are essentially inverse ill-posed source problems from the mathematical point of view, implying that small measurement errors of equivalent displacements are inherently amplified significantly such that accelerations output from accelerometers can become extremely noisy, numerically incorrect and physically meaningless in the case of high sampling rates. The ill-posedness of accelerometers cannot be rigorously solved mechanically but is always implicitly circumvented approximately. As a result, accelerometers theoretically can only produce approximate outputs of acceleration. Here we present the concept of computerized accelerometers in distinguishment of mechanical ones by applying regularization to compute accelerations from equivalent displacements. The acceleration is rigorously reconstructed mathematically as a solution to the inverse ill-posed source problem of acceleration. The concept of computerized accelerometers theoretically warrants precise measurement of acceleration without approximation, is valid for nonlinear damping as well and provides a turning point for accelerometers from a mechanical approximation to a fully rigorous computerized instrument.

Measurements in remote sensing geodesy have been well known to be of speckle noise nature. Although a number of despeckling algorithms have been proposed mainly based on the local weighted statistics in the engineering literature, there are relatively few studies on the statistical adjustment methods for processing the measurements contaminated with the speckle or multiplicative errors. We develop the least squares (LS)-based adjustment methods for the remote sensing measurements with multiplicative errors and trends, evaluate the accuracy of the parameter estimates, and derive the corresponding formulas to estimate the variance of the unit weight. Simulation examples are used to illustrate the developed theory and methods.

Akaike’s Bayesian information criterion (ABIC) has been widely used in inverse ill-posed problems. Little has been done to investigate its statistical aspects. We present an alternative derivation of the marginal distribution of measurements for ABIC under the assumption of normal distributions and show that the principle of ABIC is to statistically estimate the variances of measurements and prior data by maximizing the marginal distribution of measurements. The determination of the regularization parameter with ABIC is essentially equivalent to estimating the relative weighting between measurements and prior data. We prove that ABIC theoretically would produce a biased estimate of the variance of measurements. Since the prior mean is generally unknown but arbitrarily treated as zero in inverse ill-posed problems, ABIC is shown to fail to produce any reasonable estimate for the prior variance. Although ABIC is constructed under the Bayesian framework, it essentially plays more or less the same role as biased regularization from the frequentist’s point of view. ABIC error evaluation cannot be performed under the Bayesian framework but should be more appropriately done with the frequentist’s standpoint in terms of mean squared errors. ABIC is sensitive to prior distributions. In the case of non-informative prior distribution, ABIC leads to the conventional weighted least squares (LS) estimate of parameters and cannot be used to solve inverse ill-posed problems. It is not linked to the regularization parameter but only straightforwardly produces an unbiased estimator for the noise level of measurements, which is only applicable numerically for well-posed problems but not for inverse ill-posed problems. Numerical simulated examples are used to demonstrate the statistical performances of ABIC.

Although the weighted least squares (LS) method is straightforward to deal with errors-in-variables (EIV) models, it results in the biased estimates of parameters and the variance of unit weight. The total least squares (TLS) method is statistically rigorous and optimal but is computationally much more demanding. This paper aims at constructing an improved weighted LS estimate of parameters from the perspective of multiplicative error models, which is expected to mainly maintain the advantages of computational simplicity of the weighted LS method and the optimality of the weighted TLS method to some extent but almost remove the bias of the weighted LS estimate and free the computational burden of the TLS method. The statistical aspects of the weighted LS estimate developed from the perspective of multiplicative error models have been analyzed and an almost unbiased estimate of the variance of unit weight proposed. Finally, an N-calibrated weighted LS estimate has been constructed from the perspective of multiplicative error models.

Difference methods have been routinely used to compute velocity and acceleration from precise positioning with global navigation satellite systems (GNSS). A low sampling rate (say a rate not greater than 1 Hz, for example) has been always implicitly assumed for applicability of the methods, because random measurement errors are significantly amplified, either proportional to the sampling rate in the case of velocity or square-proportional to the sampling rate in the case of acceleration. Direct consequences of a low sampling rate are the distortion of the computed velocity and acceleration waveforms and the failure to obtain almost instantaneous values of velocity and acceleration. We reformulate the reconstruction of velocity and acceleration from very high-rate (50 Hz) precise GNSS as an inverse ill-posed problem and propose the criterion of minimum mean squared errors (MSE) to regularize solutions of velocity and acceleration. We successfully apply the MSE-based regularized method to reconstruct the very high-rate velocity and acceleration waveforms, the peak ground velocity (PGV) and the peak ground acceleration (PGA) from 50 Hz precise point positioning (PPP) position waveforms for the 2013 Lushan Mw6.6 earthquake. The reconstructed results of velocity and acceleration are shown to be in good agreement with the motion patterns in the PPP position waveforms and correctly recover the earthquake signal. The reconstructed GNSS-based PGA values are a few hundred times smaller than those from the strong motion seismometers.

Although maximum likelihood has been widely used to estimate unknown parameters in stochastic models of random errors, we show that the method cannot be used to estimate variance components...

High-rate GNSS precise point positioning (PPP) has been playing a more and more important role in providing precise positioning information in fast time-varying environments. Although kinematic PPP is commonly known to have a precision of a few centimeters, the precision of high-rate PPP within a short period of time has been reported recently with experiments to reach a few millimeters in the horizontal components and sub-centimeters in the vertical component to measure seismic motion, which is several times better than the conventional kinematic PPP practice. To fully understand the mechanism of mystified excellent performance of high-rate PPP within a short period of time, we have carried out a theoretical error analysis of PPP and conducted the corresponding simulations within a short period of time. The theoretical analysis has clearly indicated that the high-rate PPP errors consist of two types: the residual systematic errors at the starting epoch, which affect high-rate PPP through the change of satellite geometry, and the time-varying systematic errors between the starting epoch and the current epoch. Both the theoretical error analysis and simulated results are fully consistent with and thus have unambiguously confirmed the reported high precision of high-rate PPP, which has been further affirmed here by the real data experiments, indicating that high-rate PPP can indeed achieve the millimeter level of precision in the horizontal components and the sub-centimeter level of precision in the vertical component to measure motion within a short period of time. The simulation results have clearly shown that the random noise of carrier phases and higher order ionospheric errors are two major factors to affect the precision of high-rate PPP within a short period of time. The experiments with real data have also indicated that the precision of PPP solutions can degrade to the cm level in both the horizontal and vertical components, if the geometry of satellites is rather poor with a large DOP value.

Tienstra’s method was developed to solve parameter adjustment with linear equality constraints, which has been otherwise often carried out by directly applying the conventional (or standard) method of Lagrange multipliers to quadratic optimization problems with a positive definite matrix. We analyze the computational complexity of the celebrated Tienstra’s method and compare it with that of the method of Lagrange multipliers. We show that Tienstra’s method is not only statistically elegant but also interestingly of significant computational advantage over the method of Lagrange multipliers to solve weighted least squares adjustment with linear equality constraints, with the saving of computational cost by a minimum of about 38% to a maximum of 100%. Tienstra’s method can be very important for large scale problems of adjustment and inversion with linear equality constraints.

A moment magnitude Mw9.0 earthquake hit northeastern Japan at 14:46:18 (Japan Standard Time), March 11, 2011. We have obtained 1 s precise point positioning solutions for 1198 GEONET stations. Although GPS position time series have been routinely investigated and used as waveforms for dynamic inversion of earthquakes, we focus on exploring the spatial displacement features of GEONET stations for this earthquake. A movie inspection of high-rate GPS waveforms leads us to find that 76.21% of the GEONET stations in the Japanese islands subsided suddenly within 1 s between 14:59:45 and 14:59:46, Japan local time, with an average displacement of $$-\,2.43$$ , 2.83 and $$-\,4.75$$ mm in the east, north and vertical components, respectively, about 15 min after the 2011 Tohoku earthquake. We have performed different types of independent tests, namely measurement error analysis, processing the GEONET data with a different software system, a statistical hypothesis testing under a simple assumption of sign distributions, the test computation of the displacement field outside of the Japanese islands and an independent test with the Japanese strong motion borehole network KiK-net, to see whether these sudden movements actually occurred. The first four independent tests are passed almost without any doubt, and the direction of the average sudden displacements is roughly consistent tectonically with the direction of subduction of the Pacific plate. Because there are only 78 KiK-net borehole stations available for an independent seismic test, the KiK-net results are marginally consistent with those of GEONET. In the daily seismological and geophysical practice, one may then conclude that the sudden movement within the second is real after passing these five independent tests. However, a further epoch-by-epoch check pinpoints a few more seconds with even a higher probability of sudden displacement from the 20-min three-component high-rate GPS waveforms after the main shock, or more precisely, the seconds between 14:59:04 and 14:59:05, 15:01:04 and 15:01:05, and 15:03:39 and 15:03:40 with 80.80, 84.14 and 85.89% of the GEONET stations simultaneously moving upward, southward and westward, respectively. Although these probabilities are very high, it may hardly be imagined that a large scale of sudden movements could occur repeatedly between 14:59:04 and 15:03:40. The high-rate GPS results imply that some detected sudden movements after the earthquake could be unidentified artifacts of GPS data processing, though we cannot rule out the possibility that the detected sudden movements in Japan after the 2011 Tohoku Mw9.0 earthquake are real physical signals.

Although total least squares (TLS) has been widely applied, variance components in an errors-in-variables (EIV) model can be inestimable under certain conditions and unstable in the sense that small random errors can result in very large errors in the estimated variance components. We investigate the effect of the random design matrix on variance component (VC) estimation of MINQUE type by treating the design matrix as if it were errors-free, derive the first-order bias of the VC estimate, and construct bias-corrected VC estimators. As a special case, we obtain a bias-corrected estimate for the variance of unit weight. Although TLS methods are statistically rigorous, they can be computationally too expensive. We directly Taylor-expand the nonlinear weighted LS estimate of parameters up to the second-order approximation in terms of the random errors of the design matrix, derive the bias of the estimate, and use it to construct a bias-corrected weighted LS estimate. Bearing in mind that the random errors of the design matrix will create a bias in the normal matrix of the weighted LS estimate, we propose to calibrate the normal matrix by computing and then removing the bias from the normal matrix. As a result, we can obtain a new parameter estimate, which is called the N-calibrated weighted LS estimate. The simulations have shown that (i) errors-in-variables have a significant effect on VC estimation, if they are large/significant but treated as non-random. The variance components can be incorrectly estimated by more than one order of magnitude, depending on the nature of problems and the sizes of EIV; (ii) the bias-corrected VC estimate can effectively remove the bias of the VC estimate. If the signal-to-noise is small, higher order terms may be necessary. Nevertheless, since we construct the bias-corrected VC estimate by directly removing the estimated bias from the estimate itself, the simulation results have clearly indicated that there is a great risk to obtain negative values for variance components. VC estimation in EIV models remains difficult and challenging; and (iii) both the bias-corrected weighted LS estimate and the N-calibrated weighted LS estimate obviously outperform the weighted LS estimate. The intuitively N-calibrated weighted LS estimate is computationally less expensive and shown to statistically perform even better than the bias-corrected weighted LS estimate in producing an almost unbiased estimate of parameters.

This chapter discusses statistical and numerical aspects of constrained and unconstrained errors-in-variables (EIV) models. The parameters in an EIV model can often be estimated by using three categories of methods: the conventional weighted least squares (LS) method, normed orthogonal regression, and the weighted total least squares (TLS) method. The conventional weighted LS method is of significantly computational advantage but not rigorous statistically. We systematically investigate the effects of random errors in the design matrix on the weighted LS estimates of parameters and variance components, construct the N-calibrated almost unbiased weighted LS estimator of parameters and derive almost unbiased estimates for the variance of unit weight. Although orthogonal regression can be used to estimate the parameters in an EIV model, it is not statistically optimal either. The weighted TLS method is most rigorous and optimal to statistically estimate the parameters in an EIV model at the cost of substantially increasing computation. We reformulate an EIV model as a nonlinear adjustment model without constraints and investigate the statistical effects of nonlinearity on the nonlinear TLS estimate, including the first order approximation of accuracy, nonlinear confidence region and bias of the nonlinear TLS estimate. Closed form solutions to coordinate transformation have been presented as well. Finally, we prove that variance components in an EIV model with the simplest stochastic structure are not estimable.

Multiplicative error models should become more and more important in geodesy, since modern measurement technology on the basis of electromagnetic wave has clearly demonstrated that measurements of this type contain two types of random errors: fixed random errors and baseline-length dependent random errors. Although a number of the estimators of the variance of unit weight are derived from the least-squares-based adjustment methods for multiplicative error models recently, we know very little about their statistical performances. We first derive the variances of the estimates of the variance of unit weight in multiplicative error models. We find that the second order term of random errors will not affect the unbiasedness of an estimate of the variance of unit weight, if such a term is generated from the nonlinearity of models and/or least-squares-based nonlinear objective functions. The result is surprising, since the second order term of random errors has been well known to create the biases in both the estimate of parameters and the measurement corrections in the literature of nonlinear adjustment and nonlinear regression. Simulations are carried out to confirm the statistical analysis and to numerically compare the performances of different estimates of the variance of unit weight in multiplicative error models. From the simulation results, we recommend the estimate of the variance of unit weight with the bias-corrected weighted LS solution, followed by the two estimates with the ordinary LS solutions and the first estimate with the weighted LS solution.

The purpose of this note is to derive alternative formulae for parameter estimation in partial errors-in-variables models to those given in Xu et al. (J Geod 86:661–675, 2012). We show that these new formulae are more compact and more straightforward to interpret the corrections of the random elements of the design matrix. Nevertheless, each formula has its own computational advantage, depending on the relationship between the number of measurements y and that of the independent random elements of the design matrix. The original formula is computationally more efficient, if the number of measurements y is significantly larger than that of the independent random elements of the design matrix; otherwise, the new alternative formula is much more efficient.

Although total least squares (TLS) is more rigorous than the weighted least squares (LS) method to estimate the parameters in an errors-in-variables (EIV) model, it is computationally much more complicated than the weighted LS method. For some EIV problems, the TLS and weighted LS methods have been shown to produce practically negligible differences in the estimated parameters. To understand under what conditions we can safely use the usual weighted LS method, we systematically investigate the effects of the random errors of the design matrix on weighted LS adjustment. We derive the effects of EIV on the estimated quantities of geodetic interest, in particular, the model parameters, the variance–covariance matrix of the estimated parameters and the variance of unit weight. By simplifying our bias formulae, we can readily show that the corresponding statistical results obtained by Hodges and Moore (Appl Stat 21:185–195, 1972) and Davies and Hutton (Biometrika 62:383–391, 1975) are actually the special cases of our study. The theoretical analysis of bias has shown that the effect of random matrix on adjustment depends on the design matrix itself, the variance–covariance matrix of its elements and the model parameters. Using the derived formulae of bias, we can remove the effect of the random matrix from the weighted LS estimate and accordingly obtain the bias-corrected weighted LS estimate for the EIV model. We derive the bias of the weighted LS estimate of the variance of unit weight. The random errors of the design matrix can significantly affect the weighted LS estimate of the variance of unit weight. The theoretical analysis successfully explains all the anomalously large estimates of the variance of unit weight reported in the geodetic literature. We propose bias-corrected estimates for the variance of unit weight. Finally, we analyze two examples of coordinate transformation and climate change, which have shown that the bias-corrected weighted LS method can perform numerically as well as the weighted TLS method.

Although total least squares has been substantially investigated theoretically and widely applied in practical applications, almost nothing has been done to simultaneously address the estimation of parameters and the errors-in-variables (EIV) stochastic model. We prove that the variance components of the EIV stochastic model are not estimable, if the elements of the random coefficient matrix can be classified into two or more groups of data of the same accuracy. This result of inestimability is surprising as it indicates that we have no way of gaining any knowledge on such an EIV stochastic model. We demonstrate that the linear equations for the estimation of variance components could be ill-conditioned, if the variance components are theoretically estimable. Finally, if the variance components are estimable, we derive the biases of their estimates, which could be significantly amplified due to a large condition number.

Modern observation technology has verified that measurement errors can be proportional to the true values of measurements such as GPS, VLBI baselines and LiDAR. Observational models of this type are called multiplicative error models. This paper is to extend the work of Xu and Shimada published in 2000 on multiplicative error models to analytical error analysis of quantities of practical interest and estimates of the variance of unit weight. We analytically derive the variance-covariance matrices of the three least squares (LS) adjustments, the adjusted measurements and the corrections of measurements in multiplicative error models. For quality evaluation, we construct five estimators for the variance of unit weight in association of the three LS adjustment methods. Although LiDAR measurements are contaminated with multiplicative random errors, LiDAR-based digital elevation models (DEM) have been constructed as if they were of additive random errors. We will simulate a model landslide, which is assumed to be surveyed with LiDAR, and investigate the effect of LiDAR-type multiplicative error measurements on DEM construction and its effect on the estimate of landslide mass volume from the constructed DEM.

Ionospheric delay effect is a critical issue that limits the accuracy of precise Global Navigation Satellite System (GNSS) positioning and navigation for single-frequency users, especially in mid- and low-latitude regions where variations in the ionosphere are larger. Kriging spatial interpolation techniques have been recently introduced to model the spatial correlation and variability of ionosphere, which intrinsically assume that the ionosphere field is stochastically stationary but does not take the random observational errors into account. In this paper, by treating the spatial statistical information on ionosphere as prior knowledge and based on Total Electron Content (TEC) semivariogram analysis, we use Kriging techniques to spatially interpolate TEC values. By assuming that the stochastic models of both the ionospheric signals and measurement errors are only known up to some unknown factors, we propose a new Kriging spatial interpolation method with unknown variance components for both the signals of ionosphere and TEC measurements. Variance component estimation has been integrated with Kriging to reconstruct regional ionospheric delays. The method has been applied to data from the Crustal Movement Observation Network of China (CMONOC) and compared with the ordinary Kriging and polynomial interpolations with spherical cap harmonic functions, polynomial functions and low-degree spherical harmonic functions. The statistics of results indicate that the daily ionospheric variations during the experimental period characterized by the proposed approach have good agreement with the other methods, ranging from 10 to 80 TEC Unit (TECU, 1 TECU = 1 × 1016 electrons/m2) with an overall mean of 28.2 TECU. The proposed method can produce more appropriate estimations whose general TEC level is as smooth as the ordinary Kriging but with a smaller standard deviation around 3 TECU than others. The residual results show that the interpolation precision of the new proposed method is better than the ordinary Kriging and polynomial interpolation by about 1.2 TECU and 0.7 TECU, respectively. The root mean squared error of the proposed new Kriging with variance components is within 1.5 TECU and is smaller than those from other methods under comparison by about 1 TECU. When compared with ionospheric grid points, the mean squared error of the proposed method is within 6 TECU and smaller than Kriging, indicating that the proposed method can produce more accurate ionospheric delays and better estimation accuracy over China regional area.

The effects of incorrect weights on the weighted least squares (WLS) estimate, its accuracy, relative goodness and reliability have been thoroughly investigated in statistical and geodetic literature. Although the variance of unit weight has been one of the most important quantities in statistics and geodesy, little has been done to understand the effect of incorrect weights of observations and incorrect prior information on this quantity, except for an earlier study of Koch who showed that incorrect weights of observations would create a bias in the estimation of the variance of unit weight, though they do not affect the unbiasedness of the WLS estimates of parameters and the corrections to the observations. Since Koch did not directly give the formula to compute the bias, we will derive the bias of the estimated variance of unit weight due to incorrect weights of observations in this paper. In the case of incorrect prior information, both incorrect prior mean and incorrect prior weights can independently contribute a bias to the estimate of the variance of unit weight. Two simulated examples clearly show that the bias of the estimated variance of unit weight due to incorrect weights can be much larger than the variance of unit weight itself.

Adjustment has been based on the assumption that random errors of measurements are added to functional models. In geodetic practice, we know that accuracy formulae of modern geodetic measurements often consist of two parts: one proportional to the measured quantity and the other constant. From the statistical point of view, such measurements are of mixed multiplicative and additive random errors. However, almost no adjustment has been developed to strictly address geodetic data contaminated by mixed multiplicative and additive random errors from the statistical point of view. We systematically develop adjustment methods for geodetic data contaminated with multiplicative and additive errors. More precisely, we discuss the ordinary least squares (LS) and weighted LS methods and extend the bias-corrected weighted LS method of Xu and Shimada (Commun Stat B29:83–96, 2000) to the case of mixed multiplicative and additive random errors. The first order approximation of accuracy for all these three methods is derived. We derive the biases of weighted LS estimates. The three methods are then demonstrated and compared with a synthetic example of surface interpolation. The bias-corrected weighted LS estimate is unbiased up to the second order approximation and is of the best accuracy. Although the LS method can warrant an unbiased estimate for a linear model with multiplicative and additive errors, it is less accurate and always produces a very poor estimate of the variance of unit weight.

High-rate GPS has been widely used to construct displacement waveforms and to invert for source parameters of earthquakes. Almost all works on internal and external evaluation of high-rate GPS accuracy are based on GPS relative positioning. We build an experimental platform to externally evaluate the accuracy of 50-Hz PPP displacement waveforms. Since the shake table allows motion in any of six degrees of freedom, we install an inertial measurement unit (IMU) to measure the attitude of the platform and transform the IMU displacements into the GPS coordinate system. The experimental results have shown that high-rate PPP can produce absolute horizontal displacement waveforms at the accuracy of 2–4 mm and absolute vertical displacement waveforms at the sub-centimeter level of accuracy within a short period of time. The significance of the experiments indicates that high-rate PPP is capable of detecting absolute seismic displacement waveforms at the same high accuracy as GPS relative positioning techniques, but requires no fixed datum station. We have also found a small scaling error of IMU and a small time offset of misalignment between high-rate PPP and IMU displacement waveforms by comparing the amplitudes of and cross-correlating both the displacement waveforms.

The numerical integration method has been routinely used to produce global standard gravitational models from satellite tracking measurements of CHAMP/GRACE types. It is implemented by solving the differential equations of the partial derivatives of a satellite orbit with respect to the unknown harmonic coefficients under the conditions of zero initial values. From the mathematical point of view, satellite gravimetry from satellite tracking is the problem of estimating unknown parameters in the Newton's nonlinear differential equations from satellite tracking measurements. We prove that zero initial values for the partial derivatives are incorrect mathematically and not permitted physically. The numerical integration method, as currently implemented and used in satellite gravimetry and statistics, is groundless. We use three different methods to derive new local solutions to the Newton's nonlinear governing differential equations of motion with a nominal reference orbit. Bearing in mind that satellite orbits can now be tracked almost continuously at unprecedented high accuracy, we propose the measurement-based perturbation theory and derive global uniformly convergent solutions to the Newton's nonlinear governing differential equations of motion. Since the solutions are global uniformly convergent, they are able to extract smallest possible gravitational signals from modern and future satellite tracking measurements for global high-precision, high-resolution gravitational models. By directly turning the nonlinear differential equations of satellite motion into the nonlinear integral equations, we reformulate the links between satellite tracking measurements and the global uniformly convergent solutions to the Newton's governing differential equations as a condition adjustment model with equality constraints., 35 pages, 1 figure

The ionospheric delay is the predominant error source limiting the required positioning accuracy of GNSS whose estimation accuracy affects the effectiveness of ionospheric delay correction directly, especially in mid- and low-latitude regions where the irregularities of ionosphere are prominent. And a certain spatiotemporal variation exists in ionosphere which is caused by the fact that the intensity of ionospheric delay depends on the solar activities. Considering the spatiotemporal distribution characteristics of the ionosphere, data from the Crustal Movement Observation Network Of China (CMONOC) are processed to achieve the ionospheric TEC grid modeling in real time based on spatial variability and correlation theory, and its performance is analyzed under different solar activities. By processing the measured data, the accuracy of Kriging modeling under solar maximum (2014.02) and solar minimum (2015.03) since 2014 are analyzed and verified by the data from external rover stations located at different latitudes and longitudes. The statistics show that the inner precision can still reach about 2.63 TECU, and the Mean Square Error (MSE) of estimated grids VTEC is about 4.52 TECU under intense solar activity condition. The effect of correction at each rover station and the overall average of correction is generally 80 and 90 % or more, respectively, which is nearly equal to than that under low solar activity. This work reveals that the Kriging approach presented here which takes into account the spatial correlation and variability structure of the ionosphere TEC is appropriate for the areas with complex ionospheric variations such as China, and is also suitable for the situation under different solar activities.

A regularized solution is well-known to be biased. Although the biases of the estimated parameters can only be computed with the true values of parameters, we attempt to improve the regularized solution by using the regularized solution itself to replace the true (unknown) parameters for estimating the biases and then removing the computed biases from the regularized solution. We first analyze the theoretical relationship between the regularized solutions with and without the bias correction, derive the analytical conditions under which a bias-corrected regularized solution performs better than the ordinary regularized solution in terms of mean squared errors (MSE) and design the corresponding method to partially correct the biases. We then present two numerical examples to demonstrate the performance of our partially bias-corrected regularization method. The first example is mathematical with a Fredholm integral equation of the first kind. The simulated results show that the partially bias-corrected regularized solution can improve the MSE of the ordinary regularized function by 11%. In the second example, we recover gravity anomalies from simulated gravity gradient observations. In this case, our method produces the mean MSE of 3.71 mGal for the resolved mean gravity anomalies, which is better than that from the regularized solution without bias correction by 5%. The method is also shown to successfully reduce the absolute maximum bias from 13.6 to 6.8 mGal.

SUMMARY Strong motion accelerographs have been deployed worldwide to monitor the ground shaking of the Earth and the recorded accelerograms have been used to recover the velocities and displacements by integration. In spite of their fundamental importance in seismology and earthquake engineering, few works address the error estimates of the derived velocities and displacements. We show that the error estimates of the velocities and displacements obtained from accelerograms in the earthquake literature approach to zero as the sampling interval of accelerographs tends to zero; these are erroneous from the statistical point of view. As a result, we present a set of formulae to correctly estimate the errors (or variances) of the integrated velocities and displacements from accelerograms. In addition, we also derive the covariances between the velocities and displacements.

The integer least squares (ILS) problem, also known as the weighted closest point problem, is highly interdisciplinary, but no algorithm can find its global optimal integer solution in polynomial time. We first outline two suboptimal integer solutions, which can be important either in real time communication systems or to solve high dimensional GPS integer ambiguity unknowns. We then focus on the most efficient algorithm to search for the exact integer solution, which is shown to be faster than LAMBDA in the sense that the ratio of integer candidates to be checked by the efficient algorithm to those by LAMBDA can be theoretically expressed by rm, where r[m] and m is the number of integer unknowns. Finally, we further improve the searching efficiency of the most powerful combined algorithm by implementing two sorting strategies, which can either be used for finding the exact integer solution or for constructing a suboptimal integer solution. Test examples clearly demonstrate that the improved methods can perform significantly better than the most powerful combined algorithm to simultaneously find the optimal and second optimal integer solutions, if the ILS problem cannot be well reduced.

The LLL reduction of lattice vectors and its variants have been widely used to solve the weighted integer least squares (ILS) problem, or equivalently, the weighted closest point problem. Instead of reducing lattice vectors, we propose a parallel Cholesky-based reduction method for positive definite quadratic forms. The new reduction method directly works on the positive definite matrix associated with the weighted ILS problem and is shown to satisfy part of the inequalities required by Minkowski’s reduction of positive definite quadratic forms. The complexity of the algorithm can be fixed a priori by limiting the number of iterations. The simulations have clearly shown that the parallel Cholesky-based reduction method is significantly better than the LLL algorithm to reduce the condition number of the positive definite matrix, and as a result, can significantly reduce the searching space for the global optimal, weighted ILS or maximum likelihood estimate.

SUMMARY The method of generalized cross-validation (GCV) has been widely used to determine the regularization parameter, because the criterion minimizes the average predicted residuals of measured data and depends solely on data. The data-driven advantage is valid only if the variance–covariance matrix of the data can be represented as the product of a given positive definite matrix and a scalar unknown noise variance. In practice, important geophysical inverse ill-posed problems have often been solved by combining different types of data. The stochastic model of measurements in this case contains a number of different unknown variance components. Although the weighting factors, or equivalently the variance components, have been shown to significantly affect joint inversion results of geophysical ill-posed problems, they have been either assumed to be known or empirically chosen. No solid statistical foundation is available yet to correctly determine the weighting factors of different types of data in joint geophysical inversion. We extend the GCV method to accommodate both the regularization parameter and the variance components. The extended version of GCV essentially consists of two steps, one to estimate the variance components by fixing the regularization parameter and the other to determine the regularization parameter by using the GCV method and by fixing the variance components. We simulate two examples: a purely mathematical integral equation of the first kind modified from the first example of Phillips (1962) and a typical geophysical example of downward continuation to recover the gravity anomalies on the surface of the Earth from satellite measurements. Based on the two simulated examples, we extensively compare the iterative GCV method with existing methods, which have shown that the method works well to correctly recover the unknown variance components and determine the regularization parameter. In other words, our method lets data speak for themselves, decide the correct weighting factors of different types of geophysical data, and determine the regularization parameter. In addition, we derive an unbiased estimator of the noise variance by correcting the biases of the regularized residuals. A simplified formula to save the time of computation is also given. The two new estimators of the noise variance are compared with six existing methods through numerical simulations. The simulation results have shown that the two new estimators perform as well as Wahba's estimator for highly ill-posed problems and outperform any existing methods for moderately ill-posed problems.

Satellite orbits have been routinely used to produce models of the Earth’s gravity field. In connection with such productions, the partial derivatives of a satellite orbit with respect to the force parameters to be determined, namely, the unknown harmonic coefficients of the gravitational model, have been first computed by setting the initial values of partial derivatives to zero. In this note, we first design some simple mathematical examples to show that setting the initial values of partial derivatives to zero is generally erroneous mathematically. We then prove that it is prohibited physically. In other words, setting the initial values of partial derivatives to zero violates the physics of motion of celestial bodies.

Geodetic adjustment theory has been developed on the basis of a linear or nonlinear Gauss-Markov model, in which the random errors of measurements are always assumed to be independent of the true values of measurements themselves and naturally added to the functional model. However, modern geodetic instruments and geodetic imaging systems have clearly shown that the random errors of such measurements consist of two parts: one is of local nature and has nothing to do with the quantity under observation, and the other is proportional to the true value of measurement. From the statistical point of view, these two types of errors are called additive and multiplicative errors, respectively. Obviously, the conventional geodetic adjustment theory and methods for Gauss-Markov models with additive errors cannot theoretically meet the need of processing measurements contaminated by mixed additive and multiplicative random errors. This paper presents an overview of parameter estimation methods for processing mixed additive and multiplicative random errors. More specifically, we discuss two types of methods to estimate parameters in a mixed additive and multiplicative error model, namely, quasi-likelihood and least-squares-based methods. From this point of view, we extend the conventional adjustment theory and methods and give a solid theoretical foundation to process geodetic measurements contaminated by mixed additive and multiplicative random errors. Finally, we further discuss parameter estimation with prior information.

Space geodesy has brought profound convenience to our daily life with positioning-based products on one hand and provided a challenging opportunity to contribute fundamentally to mathematics and statistics on the other hand. The use of carrier phase observables has given rise to new observation models we have never encountered in any course/lecture of statistics and/or adjustment theory. This chapter is to provide a tutorial on mixed integer linear models. We first classify real-valued and mixed integer linear models and, then, accordingly define the corresponding conventional and mixed integer least squares (ILS) problems. Since the weighted ILS problem or the closest point problem is NP-hard, we start with suboptimal integer solutions. Reduction and/or decorrelation methods are briefly reviewed for practical applications. Integer unknown parameters are then exactly solved under a general framework of integer programming (aided with decorrelation and/or reduction methods) and represented/estimated from the statistical point of view. As an indispensable and fundamental element of integer least squares estimator, the Voronoi cell is shown to be extremely complex, both computationally and in shape, and has to be bounded with figures of simple shape. As a direct application, we obtain lower and upper probabilistic bounds for the probability with which the integers are correctly estimated. Finally, we briefly discuss an integer hypothesis testing problem.

Although variance and covariance components have been extensively investigated and a number of elegant formulae to compute them have been derived, nothing is known, without any ambiguity, about their estimability in the case of a fully unknown variance–covariance matrix. We prove that variance and covariance components in this case are not estimable, thus clarifying the ambiguity of the literature on the topic and correcting some erroneous statements in the literature. We also give a new theorem on the estimability of a linear function of variance and covariance components. Then we propose a new method to estimate the variance–covariance matrix with special structure, which can presumably be represented by, at most, r(r + 1)/2 independent parameters to guarantee its estimability in such a subspace, by directly implementing the positive definiteness of the matrix as constraint to the restricted maximum likelihood method, where r is the number of redundant measurements. Therefore, our estimates of the variance and covariance components always reconstruct a positive definite matrix and are always physically meaningful.

Although real-valued linear models, whether or not of full rank, have been thoroughly investigated and are well documented, very little is known about statistical and probabilistic aspects of a mixed integer linear model, which arose from space geodesy and serves as the standard starting model for precise positioning using the global positioning system (GPS). Voronoi cells play a fundamental role in the least squares estimation of the integer unknowns of the model. In this paper, we first develop a method to construct Voronoi cells and study how to fit figures of simple shape to a Voronoi cell, both from inside and outside. We then derive a number of new lower and upper bounds on the probability that the integers of the model are correctly estimated. Finally, we discuss the tests of two hypotheses on the integer mean

Regularization has been applied by implicitly assuming that the weight matrix of measurements is known. If measurements are assumed to be heteroscedastic with different unknown variance components, all regularization techniques may not be proper to apply, unless techniques of variance component estimation are directly implemented. Although variance component estimation techniques have been proposed to simultaneously estimate the variance components and provide a means of regularization, the regularization parameter is treated as if it were also an extra variance component. In this paper, we assume no prior information on the model parameters and do not treat the regularization parameter as an extra variance component. Instead, we first analyze the biases of estimated variance components due to the regularization parameter and then propose bias-corrected variance component estimators. The results have shown that they work very well. Finally, we propose and investigate through simulations an iterative scheme to simultaneously estimate the variance components and the regularization parameter, in order to eliminate the effect of regularization parameter on variance components and the effect of incorrect prior weights or initial variance components on the regularization parameter.

Kaula’s rule of thumb has been used in producing geopotential models from space geodetic measurements, including the most recent models from satellite gravity missions CHAMP. Although Xu and Rummel (Manuscr Geod 20 8–20, 1994b) suggested an alternative regularization method by introducing a number of regularization parameters, no numerical tests have ever been conducted. We have compared four methods of regularization for the determination of geopotential from precise orbits of COSMIC satellites through simulations, which include Kaula’s rule of thumb, one parameter regularization and its iterative version, and multiple parameter regularization. The simulation results show that the four methods can indeed produce good gravitational models from the precise orbits of centimetre level. The three regularization methods perform much better than Kaula’s rule of thumb by a factor of 6.4 on average beyond spherical harmonic degree 5 and by a factor of 10.2 for the spherical harmonic degrees from 8 to 14 in terms of degree variations of root mean squared errors. The maximum componentwise improvement in the root mean squared error can be up to a factor of 60. The simplest version of regularization by multiplying a positive scalar with a unit matrix is sufficient to better determine the geopotential model. Although multiple parameter regularization is theoretically attractive and can indeed eliminate unnecessary regularization for some of the harmonic coefficients, we found that it only improved its one parameter version marginally in this COSMIC example in terms of the mean squared error.

SUMMARY Stress tensors are regularly determined from fault-slip and/or earthquake focal mechanism data in structural geology and seismotectonics. The inverse problem is non-linear, unless empirical rules of rupture and friction are employed. All the methods used to date to determine the stress tensor to this non-linear inverse problem are of a local nature, and thus cannot guarantee that the global optimal stress tensor has been obtained. We apply a hybrid global optimization method to find the global optimal stress tensor. Although the inverse problem is non-linear, the effect of non-linearity on the biases of stress tensors has not been investigated. We will examine the biases of inverted stress tensors and their effect on the principal orientations of stress and the shape parameter of the stress ellipsoid. The biases of stress parameters have been shown to be comparable with the estimated stress parameters numerically. We compare the accuracy of the four stress parameters, with and without taking the errors in fault planes into account. If the errors in fault planes are not taken into account, the stress parameters are too optimistically estimated by a factor of 4 to 9 in the example. We will also mathematically reformulate the assumption that the directions of maximum shear stress represent those of slips on fault planes as two functionally independent but equivalent constraints of equality. The new formulation is computationally more effective and provides a correct method for calculating the accuracy of stress parameters.

The previous work of Xu on discrete nonlinear filtering is extended to continuous systems. The new results are summarized as follows: (1) a second-order unbiased prediction of the true state governed by a vector stochastic differential equation is worked out; (2) a set of coupled differential equations for a new truncated second-order nonlinear filter and its variance–covariance matrix are derived from the frequentist point of view. The new filter is proved to be unbiased to the second-order approximation; and, most importantly, (3) comparison of the new filtering and accuracy results with the literature on nonlinear filtering has indicated that more than 40 years of nonlinear filtering of continuous systems may have foundational problems.

Although simple probabilistic models for directions in two–dimensional (2D) and three–dimensional (3D) spaces are well known (respectively, the von Mises and Fisher distributions), probabilistic models for directions, planes and referential systems in an n –dimensional ( n D) space are either not popular or not available. The results for directions in 2D and 3D can be generalized for (i) arbitrary directions; (ii) planes; (iii) referential systems; and (iv) symmetric random tensors, and in turn, all these results can be generalized to n D spaces. In this paper we will focus on and build up isotropic probabilistic models for directions, planes and referential systems in n D space.