A spectral approach to estimation and smoothing of continuous spatial processes

This paper is concerned with developing computational methods and approximations for maximum likelihood estimation and minimum mean square error smoothing of irregularly observed two-dimensional stationary spatial processes. The approximations are based on various Fourier expansions of the covariance function of the spatial process, expressed in terms of the inverse discrete Fourier transform of the spectral density function of the underlying spatial process. We assume that the underlying spatial process is governed by elliptic stochastic partial differential equations (SPDE's) driven by a Gaussian white noise process. SPDE's have often been used to model the underlying physical phenomenon and the elliptic SPDE's are generally associated with steady-state problems.A central problem in estimation of underlying model parameters is to identify the covariance function of the process. The cumbersome exact analytical calculation of the covariance function by inverting the spectral density function of the process, has commonly been used in the literature. The present work develops various Fourier approximations for the covariance function of the underlying process which are in easily computable form and allow easy application of Newton-type algorithms for maximum likelihood estimation of the model parameters. This work also develops an iterative search algorithm which combines the Gauss-Newton algorithm and a type of generalized expectation-maximization (EM) algorithm, namely expectation-conditional maximization (ECM) algorithm, for maximum likelihood estimation of the parameters.We analyze the accuracy of the covariance function approximations for the spatial autoregressive-moving average (ARMA) models analyzed in Vecchia (1988) and illustrate the performance of our iterative search algorithm in obtaining the maximum likelihood estimation of the model parameters on simulated and actual data.     

A class of non-stationary covariance functions with compact support

This article describes the use of non-stationary covariance functions with compact support to estimate and simulate a random function. Based on the kernel convolution theory, the functions are derived by convolving hyperspheres in \(\mathbb{R}^n\) followed by a Radon transform. The order of the Radon transform controls the differentiability of the covariance functions. By varying spatially the hyperspheres radius one defines non-stationary isotropic versions of the spherical, the cubic and the penta-spherical models. Closed-form expressions for the non-stationary covariances are derived for the isotropic spherical, cubic, and penta-spherical models. Simulation of the different non-stationary models is easily obtained by weighted average of independent standard Gaussian variates in both the isotropic and the anisotropic case. The non-stationary spherical covariance model is applied to estimate the overburden thickness over an area composed of two different geological domains. The results are compared to the estimation with a single stationary model and the estimation with two stationary models, one for each geological domain. It is shown that the non-stationary model enables a reduction of the mean square error and a more realistic transition between the two geological domains.     

Precise estimation of covariance parameters in least-squares collocation by restricted maximum likelihood

Precise estimates of the covariance parameters are essential in least-squares collocation (LSC) in the case of increased accuracy requirements. This paper implements restricted maximum likelihood (REML) method for the estimation of three covariance parameters in LSC with the Gauss-Markov second-order function (GM2), which is often used in interpolation of gravity anomalies. The estimates are then validated with the use of an independent technique, which has been often omitted in the previous works that are confined to covariance parameters errors based on the information matrix. The crossvalidation of REML estimates with the use of hold-out method (HO) helps in understanding of REML estimation errors. We analyzed in detail the global minimum of negative log-likelihood function (NLLF) in the estimation of covariance parameters, as well, as the accuracy of the estimates. We found that the correlation between covariance parameters may critically contribute to the errors of their estimation. It was also found that knowing some intrinsic properties of the covariance function may help in the scoring process.     

A comparison of approaches to include outcrop information in overburden thickness estimation

The estimation of overburden sediment thickness is important in hydrogeology, geotechnics and geophysics. Usually, thickness is known precisely at a few sparse borehole data. To improve precision of estimation, one useful complementary information is the known position of outcrops. One intuitive approach is to code the outcrops as zero thickness data. A problem with this approach is that the outcrops are preferentially observed compared to other thickness information. This introduces a strong bias in the thickness estimation that kriging is not able to remove. We consider a new approach to incorporate point or surface outcrop information based on the use of a non-stationary covariance model in kriging. The non-stationary model is defined so as to restrict the distance of influence of the outcrops. Within this distance of influence, covariance parameters are assumed simple regular functions of the distance to the nearest outcrop. Outside the distance of influence of the outcrops, the thickness covariance is assumed stationary. The distance of influence is obtained thru a cross-validation. Compared to kriging based on a stationary model with or without zero thickness at outcrop locations, the non-stationary model provides more precise estimation, especially at points close to an outcrop. Moreover, the thickness map obtained with the non-stationary covariance model is more realistic since it forces the estimates to zero close to outcrops without the bias incurred when outcrops are simply treated as zero thickness in a stationary model.     

THE COLLOCATION APPROACH TO THE INVERSION OF GRAVITY DATA1

Gravity data inversion can provide valuable information on the structure of the underlying distribution of mass. The solution of the inversion of gravity data is an ill-posed problem, and many methods have been proposed for solving it using various systematic techniques. The method proposed here is a new approach based on the collocation principle, derived from the Wiener filtering and prediction theory. The natural multiplicity of the solution of the inverse gravimetric problem can be overcome only by assuming a substantially simplified model, in this case a two-layer model, i.e. with one separation surface and one density contrast only. The presence of gravity disturbance and/or outliers in the upper layer is also taken into account. The basic idea of the method is to propagate the covariance structure of the depth function of the separation surface to the covariance structure of the gravity field measured on a reference plane. This can be done since the gravity field produced by the layers is a functional (linearized) of the depth. Furthermore, in this approach, it is possible to obtain the variance of the estimation error which indicates the precision of the computed solution. The method has proved to be effective on simulated data, fulfilling the a priori hypotheses. In real cases which display the required statistical homogeneity, good preliminary solutions, useful for a further quantitative interpretation, have also been derived. A case study is discussed.     

Linear spatial interpolation: Analysis with an application to San Joaquin Valley

The properties of linear spatial interpolators of single realizations and trend components of regionalized variables are examined in this work. In the case of the single realization estimator explicit and exact expressions for the weighting vector and the variances of estimator and estimation error were obtained from a closed-form expression for the inverse of the Lagrangian matrix. The properties of the trend estimator followed directly from the Gauss-Markoff theorem. It was shown that the single realization estimator can be decomposed into two mutually orthogonal random functions of the data, one of which is the trend estimator. The implementation of liear spatial estimation was illustrated with three different methods, i.e., full information maximum likelihood (FIML), restricted maximum likelihood (RML), and Rao's minimum norm invariant quadratic unbiased estimation (MINQUE) for the single realization case and via generalized least squares (GLS) for the trend. The case study involved large correlation length-scale in the covariance of specific yield producing a nested covariance structure that was nearly positive semidefinite. The sensitivity of model parameters, i.e., drift and variance components (local and structured) to the correlation length-scale, choice of covariance model (i.e., exponential and spherical), and estimation method was examined. the same type of sensitivity analysis was conducted for the spatial interpolators. It is interesting that for this case study, characterized by a large correlation length-scale of about 50 mi (80 km), both parameter estimates and linear spatial interpolators were rather insensitive to the choice of covariance model and estimation method within the range of credible values obtained for the correlation length-scale, i.e., 40–60 mi (64–96 km), with alternative estimates falling within ±5% of each other.     

Improving the assessment of uncertainty and risk in the spatial prediction of environmental impacts: a new approach for fitting geostatistical model parameters based on dual Kriging in the presence of a trend

Modern methods of geostatistics deliver an essential contribution to Environmental Impact Assessment (EIA). These methods allow for spatial interpolation, forecast and risk assessment of expected impact during and after mining projects by integrating different sources of data and information. Geostatistical estimation and simulation algorithms are designed to provide both, a most likely forecast as well as information about the accuracy of the prediction. The representativeness of these measures depends strongly on the quality of the inferred model parameters, which are mainly defined by the parameters of the variogram or the covariance function. Available data may be sparse, trend affected and of different data type making the inference of representative geostatistical model parameters difficult. This contribution introduces a new method for best fitting of the geostatistical model parameters in the presence of a trend, which utilizes the empirical and theoretical differences between Universal Kriging and trend-predictions. The method extends well known approaches of cross validation in two aspects. Firstly, the model evaluation is not only limited to sample data locations but is performed on any prediction locations of the attribute in the domain. Secondly, it extends the measure used in cross validation, based on a single point replacement by using error curves. These allow defining rings of influence representing errors resulting from separate variogram lags. By analyzing the different variogram lags the fit of the complete covariance can be assessed and the influence of the several model parameters separated. The use of the proposed method in an EIA context is illustrated in a case study related on the prediction of mining-induced ground movements.     

Error covariance calculation for forecast bias estimation in hydrologic data assimilation

To date, an outstanding issue in hydrologic data assimilation is a proper way of dealing with forecast bias. A frequently used method to bypass this problem is to rescale the observations to the model climatology. While this approach improves the variability in the modeled soil wetness and discharge, it is not designed to correct the results for any bias. Alternatively, attempts have been made towards incorporating dynamic bias estimates into the assimilation algorithm. Persistent bias models are most often used to propagate the bias estimate, where the a priori forecast bias error covariance is calculated as a constant fraction of the unbiased a priori state error covariance. The latter approach is a simplification to the explicit propagation of the bias error covariance. The objective of this paper is to examine to which extent the choice for the propagation of the bias estimate and its error covariance influence the filter performance. An Observation System Simulation Experiment (OSSE) has been performed, in which ground water storage observations are assimilated into a biased conceptual hydrologic model. The magnitudes of the forecast bias and state error covariances are calibrated by optimizing the innovation statistics of groundwater storage. The obtained bias propagation models are found to be identical to persistent bias models. After calibration, both approaches for the estimation of the forecast bias error covariance lead to similar results, with a realistic attribution of error variances to the bias and state estimate, and significant reductions of the bias in both the estimates of groundwater storage and discharge. Overall, the results in this paper justify the use of the traditional approach for online bias estimation with a persistent bias model and a simplified forecast bias error covariance estimation.     

Calibration framework for a Kalman filter applied to a groundwater model

The paper presents a novel approach to the setup of a Kalman filter by using an automatic calibration framework for estimation of the covariance matrices. The calibration consists of two sequential steps: (1) Automatic calibration of a set of covariance parameters to optimize the performance of the system and (2) adjustment of the model and observation variance to provide an uncertainty analysis relying on the data instead of ad-hoc covariance values. The method is applied to a twin-test experiment with a groundwater model and a colored noise Kalman filter. The filter is implemented in an ensemble framework. It is demonstrated that lattice sampling is preferable to the usual Monte Carlo simulation because its ability to preserve the theoretical mean reduces the size of the ensemble needed. The resulting Kalman filter proves to be efficient in correcting dynamic error and bias over the whole domain studied. The uncertainty analysis provides a reliable estimate of the error in the neighborhood of assimilation points but the simplicity of the covariance models leads to underestimation of the errors far from assimilation points.     

Statistical Travel-time Tomography in Terms of Stacking Velocity

— Velocity evaluation is a key step in seismic analysis. The covariance of the true velocity field must be known when interpolating or simulating velocities from well measurements using geostatistical methods. In addition, inversion procedures often require information pertaining to this covariance. Traditionally it has been taken to be the covariance of stacking velocities. We present a simple example to show that this approximation can lead to significant errors. Better methods, such as those of Touati (1996) and Iooss (1998), use the variance of prestack picked travel times as a function of offset to infer that of the velocities. In this paper we extend their results on the estimation of the covariance of the reflected traveltimes, and obtain an explicit expression for the covariance of the square of the stacking slowness as a function of the covariance of the velocities. Although we are not able to invert the formula analytically to yield an explicit estimator for these parameters, the results obtained using it furnish a good and quick estimation of the velocity's covariance. This is illustrated with synthetic examples.     

A novel geostatistical approach combining Euclidean and gradual-flow covariance models to estimate fecal coliform along the Haw and Deep rivers in North Carolina

Incorporating flow in the covariance function is important for geostatistical water quality estimation that accounts for hydrological transport. Very few studies have successfully incorporated flow due to various reasons including implementation difficulties. To address this critical issue, we introduce here the first implementation of a flow weighted covariance model that uses gradual flow, and we use this model in a novel hybrid Euclidean/Gradual-flow covariance model to estimate fecal coliform along the Haw and Deep rivers from 2006 to 2010. The hybrid Euclidean/Gradual-flow model results in a 12.4% reduction in estimation mean square error compared to the Euclidean model, indicating that this is the first study to successfully incorporate gradual flow and demonstrate an improvement in estimation accuracy over the purely Euclidean approach. Furthermore, results show that the Euclidean/Gradual-flow model is more accurate and easier to implement than the Euclidean/Pipe-flow model. Our assessment found that 96 river miles were detected as being impaired according to the Euclidean/Gradual-flow method, which is more than twice the 39 river miles found according to the Euclidean estimate. These results demonstrate that the Euclidean/Gradual-flow model substantially increase the sensitivity in detecting fecal impairment, which provide critical new information for watershed management and public health protection measures.     

Geomagnetic field effect on the ionospheric contribution to the error of navigation satellite systems

When using global navigation satellite systems (GNSSs) for high-precision measurements, one should consider high-order errors. The ionospheric second-order error caused by the geomagnetic field is approximately proportional to the total electron content. Therefore, this error can be taken into account by modifying the coefficients in an “ionosphere-free” combination of GNSS measurements at two frequencies. This study checks the approximations underlying this modification. We reveal that these approximations are valid and the results depend weakly on the accuracy of ionospheric parameters used a priori for calculating the coefficients of the modified two-frequency formula. In addition, we investigate how the choice of a model of the Earth’s magnetic field affects the second-order ionospheric error.     

Unscented transformation with scaled symmetric sampling strategy for precision estimation of total least squares

The errors-in-variables (EIV) model is a nonlinear model, the parameters of which can be solved by singular value decomposition (SVD) method or the general iterative algorithm. The existing formulae for covariance matrix of total least squares (TLS) parameter estimates don’t fully consider the randomness of quantities in iterative algorithm and the biases of parameter estimates and residuals. In order to reflect more reasonable precision information for TLS adjustment, the derivative-free unscented transformation with scaled symmetric sampling strategy, i.e. scaled unscented transformation (SUT), is introduced and implemented. In this contribution, we firstly discuss the existing various solutions of TLS adjustment and covariance matrices of TLS parameter estimates and derive the general first-order approximate cofactor matrices of random quantities in TLS adjustment. Secondly, based on the combination of TLS iterative algorithm and calculation process of SUT, we design the two SUT algorithms to calculate the biases and the second-order approximate covariance matrices. Finally, the straight line fitting model and plane coordinate transformation model are used to demonstrate that applying SUT for precision estimation of TLS adjustment is feasible and effective.     

A Tailored Computation of the Mean Dynamic Topography for a Consistent Integration into Ocean Circulation Models

Geostrophic surface velocities can be derived from the gradients of the mean dynamic topography—the difference between the mean sea surface and the geoid. Therefore, independently observed mean dynamic topography data are valuable input parameters and constraints for ocean circulation models. For a successful fit to observational dynamic topography data, not only the mean dynamic topography on the particular ocean model grid is required, but also information about its inverse covariance matrix. The calculation of the mean dynamic topography from satellite-based gravity field models and altimetric sea surface height measurements, however, is not straightforward. For this purpose, we previously developed an integrated approach to combining these two different observation groups in a consistent way without using the common filter approaches (Becker et al. in J Geodyn 59(60):99–110, 2012; Becker in Konsistente Kombination von Schwerefeld, Altimetrie und hydrographischen Daten zur Modellierung der dynamischen Ozeantopographie 2012). Within this combination method, the full spectral range of the observations is considered. Further, it allows the direct determination of the normal equations (i.e., the inverse of the error covariance matrix) of the mean dynamic topography on arbitrary grids, which is one of the requirements for ocean data assimilation. In this paper, we report progress through selection and improved processing of altimetric data sets. We focus on the preprocessing steps of along-track altimetry data from Jason-1 and Envisat to obtain a mean sea surface profile. During this procedure, a rigorous variance propagation is accomplished, so that, for the first time, the full covariance matrix of the mean sea surface is available. The combination of the mean profile and a combined GRACE/GOCE gravity field model yields a mean dynamic topography model for the North Atlantic Ocean that is characterized by a defined set of assumptions. We show that including the geodetically derived mean dynamic topography with the full error structure in a 3D stationary inverse ocean model improves modeled oceanographic features over previous estimates.     

Parameter estimation in ensemble based snow data assimilation: A synthetic study

Estimating erroneous parameters in ensemble based snow data assimilation system has been given little attention in the literature. Little is known about the related methods’ effectiveness, performance, and sensitivity to other error sources such as model structural error. This research tackles these questions by running synthetic one-dimensional snow data assimilation with the ensemble Kalman filter (EnKF), in which both state and parameter are simultaneously updated. The first part of the paper investigates the effectiveness of this parameter estimation approach in a perfect-model-structure scenario, and the second part focuses on its dependence on model structure error. The results from first part research demonstrate the advantages of this parameter estimation approach in reducing the systematic error of snow water equivalent (SWE) estimates, and retrieving the correct parameter value. The second part results indicate that, at least in our experiment, there is an evident dependence of parameter search convergence on model structural error. In the imperfect-model-structure run, the parameter search diverges, although it can simulate the state variable well. This result suggest that, good data assimilation performance in estimating state variables is not a sufficient indicator of reliable parameter retrieval in the presence of model structural error. The generality of this conclusion needs to be tested by data assimilation experiments with more complex structural error configurations.     

Evaluating the utility of the ensemble transform Kalman filter for adaptive sampling when updating a hydrodynamic model

This paper compares two Monte Carlo sequential data assimilation methods based on the Kalman filter, for estimating the effect of measurements on simulations of state error variance made by a one-dimensional hydrodynamic model. The first method used an ensemble Kalman filter (EnKF) to update state estimates, which were then used as initial conditions for further simulations. The second method used an ensemble transform Kalman filter (ETKF) to quickly estimate the effect of measurement error covariance on forecast error covariance without the need to re-run the simulation model. The ETKF gave an unbiased estimate of EnKF analysed error variance, although differences in the treatment of measurement errors meant the results were not identical. Estimates of forecast error variance could also be made, but their accuracy deteriorated as the time from measurements increased due in part to model non-linearity and the decreasing signal variance. The motivation behind the study was to assess the ability of the ETKF to target possible measurements, as part of an adaptive sampling framework, before they are assimilated by an EnKF-based forecasting model on the River Crouch, Essex, UK. The ETKF was found to be a useful tool for quickly estimating the error covariance expected after assimilating measurements into the hydrodynamic model. It, thus, provided a means of quantifying the ‘usefulness’ (in terms of error variance) of possible sampling schemes.     

Norm-dependent covariance permissibility of weakly homogeneous spatial random fields and its consequences in spatial statistics

 Permissibility of a covariance function (in the sense of Bochner) depends on the norm (or metric) that determines spatial distance in several dimensions. A covariance function that is permissible for one norm may not be so for another. We prove that for a certain class of covariances of weakly homogeneous random fields, the spatial distance can be defined only in terms of the Euclidean norm. This class includes commonly used covariance functions. Functions that do not belong to this class may be permissible covariances for some non-Euclidean metric. Thus, a different class of covariances, for which non-Euclidean norms are valid spatial distances, is also discussed. The choice of a coordinate system and associated norm to describe a physical phenomenon depends on the nature of the properties being described. Norm-dependent permissibility analysis has important consequences in spatial statistics applications (e.g., spatial estimation or mapping), in which one is concerned about the validity of covariance functions associated with a physically meaningful norm (Euclidean or non-Euclidean).     

A mathematical model to predict the inelastic response of a steel frame: Formulation of the model

The purpose of this research is to use data from experiments to formulate a mathematical model that will predict the non-linear response of a single-storey steel frame to an earthquake input. The process used in this formulation is system identification. The form of the model is a second-order non-linear differential equation with linear viscous damping and Ramberg—Osgood type hysteresis. The damping coefficient and the three parameters in the hysteretic model are to be established. An integral weighted mean squared error function is used to evaluate the [goodness of fit] between the model's response and the structure's response when both are subjected to the same excitation. The function includes errors in displacement and acceleration and is integrated from zero to a time T, which may be the full duration of the recorded response or only a portion of it. The parameters are adjusted using a modified Gauss-Newton method until the error function is minimized. The computer program incorporating these steps in the system identification process is verified with simulated data. Results given in the paper show that in every case the program converges in few iterations to the assigned set of parameters.     

Simulation of normal distributed smooth fields by Karhunen-Loéve expansion in combination with kriging

Simulation of multigaussian stochastic fields can be made after a Karhunen-Loéve expansion of a given covariance function. This method is also called simulation by Empirical Orthogonal Functions. The simulations are made by drawing stochastic coefficients from a random generator. These numbers are multiplied with eigenfunctions and eigenvalues derived from the predefined covariance model. The number of eigenfunctions necessary to reproduce the stochastic process within a predefined variance error, turns out to be a cardinal question. Some ordinary analytical covariance functions are used to evaluate how quickly the series of eigenfunctions can be truncated. This analysis demonstrates extremely quick convergence to 99.5% of total variance for the 2nd order exponential (‘gaussian’) covariance function, while the opposite is true for the 1st order exponential covariance function. Due to these convergence characteristics, the Karhunen-Loéve method is most suitable for simulating smooth fields with ‘gaussian’ shaped covariance functions. Practical applications of Karhunen-Loéve simulations can be improved by spatial interpolation of the eigenfunctions. In this paper, we suggest interpolation by kriging and limits for reproduction of the predefined covariance functions are evaluated.     

Simulation of normal distributed smooth fields by Karhunen-Loéve expansion in combination with kriging

Simulation of multigaussian stochastic fields can be made after a Karhunen-Loéve expansion of a given covariance function. This method is also called simulation by Empirical Orthogonal Functions. The simulations are made by drawing stochastic coefficients from a random generator. These numbers are multiplied with eigenfunctions and eigenvalues derived from the predefined covariance model. The number of eigenfunctions necessary to reproduce the stochastic process within a predefined variance error, turns out to be a cardinal question. Some ordinary analytical covariance functions are used to evaluate how quickly the series of eigenfunctions can be truncated. This analysis demonstrates extremely quick convergence to 99.5% of total variance for the 2nd order exponential (‘gaussian’) covariance function, while the opposite is true for the 1st order exponential covariance function. Due to these convergence characteristics, the Karhunen-Loéve method is most suitable for simulating smooth fields with ‘gaussian’ shaped covariance functions. Practical applications of Karhunen-Loéve simulations can be improved by spatial interpolation of the eigenfunctions. In this paper, we suggest interpolation by kriging and limits for reproduction of the predefined covariance functions are evaluated.