Means and Covariance Functions for Geostatistical Compositional Data: an Axiomatic Approach

This work focuses on the characterization of the central tendency of a sample of compositional data. It provides new results about theoretical properties of means and covariance functions for compositional data, with an axiomatic perspective. Original results that shed new light on geostatistical modeling of compositional data are presented. As a first result, it is shown that the weighted arithmetic mean is the only central tendency characteristic satisfying a small set of axioms, namely continuity, reflexivity, and marginal stability. Moreover, this set of axioms also implies that the weights must be identical for all parts of the composition. This result has deep consequences for spatial multivariate covariance modeling of compositional data. In a geostatistical setting, it is shown as a second result that the proportional model of covariance functions (i.e., the product of a covariance matrix and a single correlation function) is the only model that provides identical kriging weights for all components of the compositional data. As a consequence of these two results, the proportional model of covariance function is the only covariance model compatible with reflexivity and marginal stability.     

Calibration of imperfect models to biased observations

The problem of assimilating biased and inaccurate observations into inadequate models of the physical systems from which the observations were taken is common in the petroleum and groundwater fields. When large amounts of data are assimilated without accounting for model error and observation bias, predictions tend to be both overconfident and incorrect. In this paper, we propose a workflow for calibration of imperfect models to biased observations that involves model construction, model calibration, model criticism and model improvement. Model criticism is based on computation of model diagnostics which provide an indication of the validity of assumptions. During the model improvement step, we advocate identification of additional physically motivated parameters based on examination of data mismatch after calibration and addition of bias correction terms. If model diagnostics indicates the presence of residual model error after parameters have been added, then we advocate estimation of a “total” observation error covariance matrix, whose purpose is to reduce weighting of observations that cannot be matched because of deficiency of the model. Although the target applications of this methodology are in the subsurface, we illustrate the approach with two simplified examples involving prediction of the future velocity of fall of a sphere from models calibrated to a short-time series of biased measurements with independent additive random noise. The models into which the data are assimilated contain model errors due to neglect of physical processes and neglect of uncertainty in parameters. In every case, the estimated total error covariance is larger than the true observation covariance implying that the observations need not be matched to the accuracy of the measuring instrument. Predictions are much improved when all model improvement steps were taken.     

BLP approach for the estimation of variogram parameters

Variograms of hydrologic characteristics are usually obtained by estimating the experimental variogram for distinct lag classes by commonly used estimators and fitting a suitable function to these estimates. However, these estimators may fail the conditionally positive-definite property and the better results for the statistics of cross-validation, which are two essential conditions for choosing a valid variogram model. To satisfy these two conditions, a multi-objective bilevel programming estimator (MOBLP) which is based on the process of cross-validation has been developed for better estimate of variogram parameters. This model is illustrated with some rainfall data from Luan River Basin in China. The case study demonstrated that MOBLP is an effective way to achieve a valid variogram model.     

Analysis of human risks due to dam-break floods—part 1: a new model based on Bayesian networks

Dam breaks have catastrophic consequences for human lives. This paper presents a new human risk analysis model (HURAM) using Bayesian networks for estimating human risks due to dam-break floods. A Bayesian network is constructed according to a logic structure of loss-of-life mechanisms. The nodes (parameters) and the arcs (inter-relationships) of the network are quantified with historical data, existing models and physical analyses. A dataset of 343 dam-failure cases with records of fatality is compiled for this purpose. Comparison between two existing models and the new model is made to test the new model. Finally, sensitivity analysis is conducted to identify the important parameters that lead to loss of life. The new model is able to take into account a large number of important parameters and their inter-relationships in a systematic structure; include the uncertainties of these parameters and their inter-relationships; incorporate information derived from physical analysis, empirical models and historical data; and update the predictions when information in specific cases is available. The application of this model to the study of human risks in a specific dam-break case is presented in a companion paper.     

High-resolution regional gravity field model of Pakistan based on best residual terrain model (RTM) of topography and interpolation techniques

This study focuses on the development of absolute gravity model for Pakistan based on best possible residual terrain model of gravity using residual terrain modeling technique. The datasets used for the development of model are observed gravity, global gravity models, and Shuttle Radar Topographic Mission (SRTM30) elevation data. The residual terrain modeling technique has been used in the remove-restore procedure for smoothing the observed gravity field. Different topographic elevation models were tested in the model selection and one best possible model with minimum mean and standard deviation was selected for residual terrain effects. Least square collocation technique has been used for quality control and error estimates. The best possible covariance model was established from residual gravity for onward prediction of gravity anomalies at the earth surface for error and prediction analysis. The residual terrain effect of gravity, value of free air anomaly from EGM96, and observed free air anomaly are added to normal gravity to compute the absolute gravity at earth surface. The prediction of these parameters is made by employing Lagrange interpolation with least square adjustment. The results are compared with ～5% randomly selected data points not utilized for the development of covariance function and/or model development. Spline interpolation technique has also been used for the prediction of gravity field-related parameters. Lagrange interpolation exhibits relatively superior results over spline-based interpolation. This is as per expectation due to the reason that additional gridding for spline interpolation filters the signal part as well. This fact is evident from the results of spline interpolation of Grid-I and Grid-II with relatively better prediction results in Grid-I. This version of the model is capable of prediction having limiting error of 30 mGal. The predicted results show that 96.16% of prediction data falls within above-mentioned limit with Lagrange interpolation technique with least square adjustment for whole Pakistan area. The adverse effect of gridding is absent in case of Grid-I due to relatively flat areas and predicted data matches totally with control values for both spline as well as Lagrange interpolations. However, in case of Grid-II which includes high mountains of Himalaya, gridding effect is present and the accuracy of the predicted results falls to ～92%. The computed results have been compared with absolute values predicted using EGM96 and EGM2008 models as well. The gravity field recovered with PAKGM model is much better, i.e., ～ 96.16%, than both with EGM96 and EGM2008 which is about 85% only.     

Polynomial expansions in likelihoods for spatial data: A case study

This paper illustrates the computational benefits of polynomial representations for quantities in the likelihood function for the spatial linear model based on the power covariance scheme. These benefits include a comprehensive study of likelihoods and maximum likelihood estimators for data. For simplicity, we focus on a relatively simple covariance scheme and data observed at equal intervals along a transect; we briefly indicate how generalizations to more complicated covariance functions and higher dimensions will operate.     

Orthonormal residuals in geostatistics: Model criticism and parameter estimation

In linear geostatistics, models for the mean function (drift) and the variogram or generalized covariance function are selected on the basis of the modeler's understanding of the phenomenon studied as well as data. One can seldom be assured that the most appropriate model has been selected; however, analysis of residuals is helpful in diagnosing whether some important characteristic of the data has been neglected and, ultimately, in providing a reasonable degree of assurance that the selected model is consistent with the available information. The orthonormal residuals presented in this work are kriging errors constructed so that, when the correct model is used, they are uncorrelated and have zero mean and unit variance. It is suggested that testing of orthonormal residuals is a practical way for evaluating the agreement of the model with the data and for diagnosing model deficiencies. Their advantages over the usually employed standardized residuals are discussed. A set of tests are presented. Orthonormal residuals can also be useful in the estimation of the covariance (or variogram) parameters for a model that is considered correct.     

Site response analysis using one-dimensional equivalent-linear method and Bayesian filtering

Site response analysis is crucial to define the seismic hazard and distribution of damage during earthquakes. The equivalent-linear (EQL) is a numerical method widely investigated and used for site response analysis. Because several sources of uncertainty are involved in this type of analysis, parameters defining the numerical models need to be identified from in-situ measurements. In this paper, a Bayesian inference method to estimate the expected values and covariance matrix of the model parameters is presented. The methodology uses data from downhole arrays recorded during earthquakes. Two numerical applications show the good performance and prediction capabilities of the proposed approach.     

Quadratic covariance estimation and equivalence of predictions

This article illustrates the use of linear and nonlinear regression models to obtain quadratic estimates of covariance parameters. These models lead to new insights into the motivation behind estimation methods, the relationships between different methods, and the relationship of covariance estimation to prediction. In particular, we derive the standard estimating equations for minimum norm quadratic unbiased translation invariant estimates (MINQUEs) from an appropriate linear model. Connections between the linear model, minimum variance quadratic unbiased translation invariant estimates (MIVQUEs), and MINQUEs are examined and we provide a minimum norm justification for the use of one-step normal theory maximum likelihood estimates. A nonlinear regression model is used to define MINQUEs for nonlinear covariance structures and obtain REML estimates. Finally, the equivalence of predictions under various models is examined when covariance parameters are estimated. In particular, we establish that when using MINQUE, iterative MINQUE, or restricted maximum likelihood (REML) estimates, the choice between a stationary covariance function and an intrinsically stationary semivariogram is irrelevant to predictions and estimated prediction variances.     

Euclidean Distance Matrix Analysis (EDMA): Estimation of mean form and mean form difference

Euclidean Distance Matrix Analysis (EDMA) of form is a coordinate free approach to the analysis of form using landmark data. In this paper, the problem of estimation of mean form, variance-covariance matrix, and mean form difference under the Gaussian perturbation model is considered using EDMA. The suggested estimators are based on the method of moments. They are shown to be consistent, that is as the sample size increases these estimators approach the true parameters. They are also shown to be computationally very simple. A method to improve their efficiency is suggested. Estimation in the presence of missing data is studied. In addition, it is shown that the superimposition method of estimation leads to incorrect mean form and variance-covariance structure.     

Time-dependent Inversion of Surface Subsidence due to Dynamic Reservoir Compaction

We introduce a novel, time-dependent inversion scheme for resolving temporal reservoir pressure drop from surface subsidence observations (from leveling or GPS data, InSAR, tiltmeter monitoring) in a single procedure. The theory is able to accommodate both the absence of surface subsidence estimates at sites at one or more epochs as well as the introduction of new sites at any arbitrary epoch. Thus, all observation sites with measurements from at least two epochs are utilized. The method uses both the prior model covariance matrix and the data covariance matrix, which incorporates the spatial and temporal correlations between model parameters and data, respectively. The incorporation of the model covariance implicitly guarantees smoothness of the model estimate, while maintaining specific geological features like sharp boundaries. Taking these relations into account through the model covariance matrix enhances the influence of the data on the inverted model estimate. This leads to a better defined and interpretable model estimate. The time-dependent aspect of the method yields a better constrained model estimate and makes it possible to identify non-linear acceleration or delay in reservoir compaction. The method is validated by a synthetic case study based on an existing gas reservoir with a highly variable transmissibility at the free water level. The prior model covariance matrix is based on a Monte Carlo simulation of the geological uncertainty in the transmissibility.     

Markov Models for Cross-Covariances

Markov models based on various data screening hypotheses are often used because they reduce the statistical inference burden. In the case of co-located cokriging, the commonly used Markov model results in the cross-covariance being proportional to the primary covariance. Such model is inappropriate in the presence of a smoothly varying secondary variable defined on a much larger volume support than the primary variable. For such cases, an alternative Markov screening hypothesis is proposed that results in a more continuous cross-covariance proportional to the secondary covariance model. A parallel development of both Markov models is presented. A companion paper provides a comparative application to a real data set.     

The Second-Order Stationary Universal Kriging Model Revisited

Universal kriging originally was developed for problems of spatial interpolation if a drift seemed to be justified to model the experimental data. But its use has been questioned in relation to the bias of the estimated underlying variogram (variogram of the residuals), and furthermore universal kriging came to be considered an old-fashioned method after the theory of intrinsic random functions was developed. In this paper the model is reexamined together with methods for handling problems in the inference of parameters. The efficiency of the inference of covariance parameters is shown in terms of bias, variance, and mean square error of the sampling distribution obtained by Monte Carlo simulation for three different estimators (maximum likelihood, bias corrected maximum likelihood, and restricted maximum likelihood). It is shown that unbiased estimates for the covariance parameters may be obtained but if the number of samples is small there can be no guarantee of good estimates (estimates close to the true value) because the sampling variance usually is large. This problem is not specific to the universal kriging model but rather arises in any model where parameters are inferred from experimental data. The validity of the estimates may be evaluated statistically as a risk function as is shown in this paper.     

生态系统碳通量估算中耦合涡度协方差与遥感技术研究进展

陆地生态系统CO^₂和水热通量的长期观测研究一直是国际上关注的热点问题。截止目前,利用微气象学原理的涡度协方差技术是唯一能直接测定生物圈与大气间物质与能量通量的标准方法,成为国际通量观测网络的主要技术。但是涡度协方差技术的测定仍然是一种小尺度观测方法,其观测结果难于直接外推到更大尺度。同时,缺乏区域、跨尺度生态系统及其时空动态观测数据一直是限制碳循环研究的主要障碍,而遥感技术的发展可望在不远的将来使大尺度、高分辨生态系统变化的长期定量观测成为可能。这些问题在当今集中体现在如何建立通量—遥感的跨尺度观测体系,并有效地将有限的通量站点测量数据与大尺度遥感资料以及生态模型有机地结合。总结过去耦合涡度协方差技术与遥感技术的工作,主要在以下3个层面展开：①涡度协方差技术与遥感技术对碳通量估算的相互验证;②涡度协方差技术为遥感反演提供地面参数;③遥感观测解译辅助分析通量贡献区（footprint）。集中在这3个方面进行探讨,通过总结各方面的研究特点与进展,可望为未来在这个领域开展工作理顺思路。     

Dependence of Bayesian Model Selection Criteria and Fisher Information Matrix on Sample Size

Geostatistical analyses require an estimation of the covariance structure of a random field and its parameters jointly from noisy data. Whereas in some cases (as in that of a Matérn variogram) a range of structural models can be captured with one or a few parameters, in many other cases it is necessary to consider a discrete set of structural model alternatives, such as drifts and variograms. Ranking these alternatives and identifying the best among them has traditionally been done with the aid of information theoretic or Bayesian model selection criteria. There is an ongoing debate in the literature about the relative merits of these various criteria. We contribute to this discussion by using synthetic data to compare the abilities of two common Bayesian criteria, BIC and KIC, to discriminate between alternative models of drift as a function of sample size when drift and variogram parameters are unknown. Adopting the results of Markov Chain Monte Carlo simulations as reference we confirm that KIC reduces asymptotically to BIC and provides consistently more reliable indications of model quality than does BIC for samples of all sizes. Practical considerations often cause analysts to replace the observed Fisher information matrix entering into KIC with its expected value. Our results show that this causes the performance of KIC to deteriorate with diminishing sample size. These results are equally valid for one and multiple realizations of uncertain data entering into our analysis. Bayesian theory indicates that, in the case of statistically independent and identically distributed data, posterior model probabilities become asymptotically insensitive to prior probabilities as sample size increases. We do not find this to be the case when working with samples taken from an autocorrelated random field.     

Characterization and Quantification of Uncertainty in?Coregionalization Analysis

Coregionalization analysis has been presented as a method of multi-scale analysis for multivariate spatial data. Despite an increasing use of this method in environmental and earth sciences, the uncertainty associated with the estimation of parameters in coregionalization analysis (e.g., sills and functions of sills) is potentially high and has not yet been characterized. This article aims to discuss the theory underlying coregionalization analysis and assess the robustness and limits of the method. A theoretical framework is developed to calculate the ergodic and fluctuation variance-covariance matrices of least-squares estimators of sills in the linear model of coregionalization. To adjust for the positive semidefiniteness constraint on estimated coregionalization matrices, a confidence interval estimation procedure for sills and functions of sills is presented. Thereafter, the relative importance of uncertainty measures (bias and variance) for sills and structural coefficients of correlation and determination is assessed under different scenarios to identify factors controlling their uncertainty. Our results show that the sampling grid density, the choice of the least-squares estimator of sills, the positive semidefiniteness constraint, the presence of scale dependence in the correlations, and the number and range of variogram models, all affect the level of uncertainty, sometimes through multiple interactions. The asymptotic properties of variogram model parameter estimators in a bounded sampling domain impose a theoretical limit to their accuracy and precision. Because of this limit, the uncertainty was found to be high for several scenarios, especially with three variogram models, and was often more dependent on the ratio of variogram range to domain extent than on the sampling grid density. In practice, in the coregionalization analysis of a real dataset, the circular requirement for sill estimates in the calculation of uncertainty measures makes the quantification of uncertainty very problematic, if not impossible. The use of coregionalization analysis must be made with due knowledge of the uncertainty levels and limits of the method.     

Maximum Likelihood Estimation of Spatial Covariance Parameters

In this paper, the maximum likelihood method for inferring the parameters of spatial covariances is examined. The advantages of the maximum likelihood estimation are discussed and it is shown that this method, derived assuming a multivariate Gaussian distribution for the data, gives a sound criterion of fitting covariance models irrespective of the multivariate distribution of the data. However, this distribution is impossible to verify in practice when only one realization of the random function is available. Then, the maximum entropy method is the only sound criterion of assigning probabilities in absence of information. Because the multivariate Gaussian distribution has the maximum entropy property for a fixed vector of means and covariance matrix, the multinormal distribution is the most logical choice as a default distribution for the experimental data. Nevertheless, it should be clear that the assumption of a multivariate Gaussian distribution is maintained only for the inference of spatial covariance parameters and not necessarily for other operations such as spatial interpolation, simulation or estimation of spatial distributions. Various results from simulations are presented to support the claim that the simultaneous use of maximum likelihood method and the classical nonparametric method of moments can considerably improve results in the estimation of geostatistical parameters.     

Outlier Detection for Compositional Data Using Robust Methods

Outlier detection based on the Mahalanobis distance (MD) requires an appropriate transformation in case of compositional data. For the family of logratio transformations (additive, centered and isometric logratio transformation) it is shown that the MDs based on classical estimates are invariant to these transformations, and that the MDs based on affine equivariant estimators of location and covariance are the same for additive and isometric logratio transformation. Moreover, for 3-dimensional compositions the data structure can be visualized by contour lines. In higher dimension the MDs of closed and opened data give an impression of the multivariate data behavior.     

An application of the deterministic variogram to design-based variance estimation

When estimating the mean value of a variable, or the total amount of a resource, within a specified region it is desirable to report an estimated standard error for the resulting estimate. If the sample sites are selected according to a probability sampling design, it usually is possible to construct an appropriate design-based standard error estimate. One exception is systematic sampling for which no such standard error estimator exists. However, a slight modification of systematic sampling, termed 2-step tessellation stratified (2TS) sampling, does permit the estimation of design-based standard errors. This paper develops a design-based standard error estimator for 2TS sampling. It is shown that the Taylor series approximation to the variance of the sample mean under 2TS sampling may be expressed in terms of either a deterministic variogram or a deterministic covariance function. Variance estimation then can be approached through the estimation of a variogram or a covariance function. The resulting standard error estimators are compared to some more traditional variance estimators through a simulation study. The simulation results show that estimators based on the new approach may perform better than traditional variance estimators.     

Multivariate Spatial Modeling for Geostatistical Data Using Convolved Covariance Functions

Soil pollution data collection typically studies multivariate measurements at sampling locations, e.g., lead, zinc, copper or cadmium levels. With increased collection of such multivariate geostatistical spatial data, there arises the need for flexible explanatory stochastic models. Here, we propose a general constructive approach for building suitable models based upon convolution of covariance functions. We begin with a general theorem which asserts that, under weak conditions, cross convolution of covariance functions provides a valid cross covariance function. We also obtain a result on dependence induced by such convolution. Since, in general, convolution does not provide closed-form integration, we discuss efficient computation. We then suggest introducing such specification through a Gaussian process to model multivariate spatial random effects within a hierarchical model. We note that modeling spatial random effects in this way is parsimonious relative to say, the linear model of coregionalization. Through a limited simulation, we informally demonstrate that performance for these two specifications appears to be indistinguishable, encouraging the parsimonious choice. Finally, we use the convolved covariance model to analyze a trivariate pollution dataset from California.