Bayesian updating revisited

Bayesian updating methods provide an alternate philosophy to the characterization of the input variables of a stochastic mathematical model. Here, a priori values of statistical parameters are assumed on subjective grounds or by analysis of a data base from a geologically similar area. As measurements become available during site investigations, updated estimates of parameters characterizing spatial variability are generated. However, in solving the traditional updating equations, an updated covariance matrix may be generated that is not positive-definite, particularly when observed data errors are small. In addition, measurements may indicate that initial estimates of the statistical parameters are poor. The traditional procedure does not have a facility to revise the parameter estimates before the update is carried out. alternatively, Bayesian updating can be viewed as a linear inverse problem that minimizes a weighted combination of solution simplicity and data misfit. Depending on the weight given to the a priori information, a different solution is generated. A Bayesian updating procedure for log-conductivity interpolation that uses a singular value decomposition (SVD) is presented. An efficient and stable algorithm is outlined that computes the updated log-conductivity field and the a posteriori covariance of the estimated values (estimation errors). In addition, an information density matrix is constructed that indicates how well predicted data match observations. Analysis of this matrix indicates the relative importance of the observed data. The SVD updating procedure is used to interpolate the log-conductivity fields of a series of hypothetical aquifers to demonstrate pitfalls and possibilities of the method.     

Bayesian updating revisited

Bayesian updating methods provide an alternate philosophy to the characterization of the input variables of a stochastic mathematical model. Here, a priori values of statistical parameters are assumed on subjective grounds or by analysis of a data base from a geologically similar area. As measurements become available during site investigations, updated estimates of parameters characterizing spatial variability are generated. However, in solving the traditional updating equations, an updated covariance matrix may be generated that is not positive-definite, particularly when observed data errors are small. In addition, measurements may indicate that initial estimates of the statistical parameters are poor. The traditional procedure does not have a facility to revise the parameter estimates before the update is carried out. alternatively, Bayesian updating can be viewed as a linear inverse problem that minimizes a weighted combination of solution simplicity and data misfit. Depending on the weight given to the a priori information, a different solution is generated. A Bayesian updating procedure for log-conductivity interpolation that uses a singular value decomposition (SVD) is presented. An efficient and stable algorithm is outlined that computes the updated log-conductivity field and the a posteriori covariance of the estimated values (estimation errors). In addition, an information density matrix is constructed that indicates how well predicted data match observations. Analysis of this matrix indicates the relative importance of the observed data. The SVD updating procedure is used to interpolate the log-conductivity fields of a series of hypothetical aquifers to demonstrate pitfalls and possibilities of the method.     

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.     

Singularity and Nonnormality in the Classification of Compositional Data

Geologists may want to classify compositional data and express the classification as a map. Regionalized classification is a tool that can be used for this purpose, but it incorporates discriminant analysis, which requires the computation and inversion of a covariance matrix. Covariance matrices of compositional data always will be singular (noninvertible) because of the unit-sum constraint. Fortunately, discriminant analyses can be calculated using a pseudo-inverse of the singular covariance matrix; this is done automatically by some statistical packages such as SAS. Granulometric data from the Darss Sill region of the Baltic Sea is used to explore how the pseudo-inversion procedure influences discriminant analysis results, comparing the algorithm used by SAS to the more conventional Moore–Penrose algorithm. Logratio transforms have been recommended to overcome problems associated with analysis of compositional data, including singularity. A regionalized classification of the Darss Sill data after logratio transformation is different only slightly from one based on raw granulometric data, suggesting that closure problems do not influence severely regionalized classification of compositional data.     

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.     

A nonparametric view of generalized covariances for kriging

Fitting trend and error covariance structure iteratively leads to bias in the estimated error variogram. Use of generalized increments overcomes this bias. Certain generalized increments yield difference equations in the variogram which permit graphical checking of the model. These equations extend to the case where errors are intrinsic random functions of order k, k=1, 2, ..., and an unbiased nonparametric graphical approach for investigating the generalized covariance function is developed. Hence, parametric models for the generalized covariance produced by BLUEPACK-3D or other methods may be assessed. Methods are illustrated on a set of coal ash data and a set of soil pH data.     

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.     

Bayesian estimation of reservoir properties—effects of uncertainty quantification of 4D seismic data

This paper shows a history matching workflow with both production and 4D seismic data where the uncertainty of seismic data for history matching comes from Bayesian seismic waveform inversion. We use a synthetic model and perform two seismic surveys, one before start of production and the second after 1 year of production. From the first seismic survey, we estimate the contrast in slowness squared (with uncertainty) and use this estimate to generate an initial estimate of porosity and permeability fields. This ensemble is then updated using the second seismic survey (after inversion to contrasts) and production data with an iterative ensemble smoother. The impact on history matching results from using different uncertainty estimates for the seismic data is investigated. From the Bayesian seismic inversion, we get a covariance matrix for the uncertainty and we compare using the full covariance matrix with using only the diagonal. We also compare with using a simplified uncertainty estimate that does not come from the seismic inversion. The results indicate that it is important not to underestimate the noise in seismic data and that having information about the correlation in the error in seismic data can in some cases improve the results.     

Variance–Covariance Matrix of the Experimental Variogram: Assessing Variogram Uncertainty

Assessment of the sampling variance of the experimental variogram is an important topic in geostatistics as it gives the uncertainty of the variogram estimates. This assessment, however, is repeatedly overlooked in most applications mainly, perhaps, because a general approach has not been implemented in the most commonly used software packages for variogram analysis. In this paper the authors propose a solution that can be implemented easily in a computer program, and which, subject to certain assumptions, is exact. These assumptions are not very restrictive: second-order stationarity (the process has a finite variance and the variogram has a sill) and, solely for the purpose of evaluating fourth-order moments, a Gaussian distribution for the random function. The approach described here gives the variance–covariance matrix of the experimental variogram, which takes into account not only the correlation among the experiemental values but also the multiple use of data in the variogram computation. Among other applications, standard errors may be attached to the variogram estimates and the variance–covariance matrix may be used for fitting a theoretical model by weighted, or by generalized, least squares. Confidence regions that hold a given confidence level for all the variogram lag estimates simultaneously have been calculated using the Bonferroni method for rectangular intervals, and using the multivariate Gaussian assumption for K-dimensional elliptical intervals (where K is the number of experimental variogram estimates). A general approach for incorporating the uncertainty of the experimental variogram into the uncertainty of the variogram model parameters is also shown. A case study with rainfall data is used to illustrate the proposed approach.     

基于不确定性分析框架的动态环状河网水质模型——以温州市温瑞塘河为例

以温州市温瑞塘河为背景,基于不确定性分析的框架,开发了动态环状河网水质模型。用HSY(Hornberger,Spear and Young)算法作水质参数的不确定性分析,求得模型参数的空间分布,从而提高模型使用的可靠性、降低决策的风险度。模型由水量子模型及水质子模型两部分组成:水量子模型采用圣维南方程并用四级解法求解;水质子模型采用CSTR(Continuously Stirred Tank Reactor)模型的机理,并结合环状河网特征作了修正,由于只需解常微分方程,因而避免了矩阵求解,显著提高了计算效率,使得用HSY算法进行不确定性分析成为可能。模型选取温州市温瑞塘河流域的鹿城区河网进行了首次应用,用HSY算法率定水质参数并讨论了参数的不确定性。对参数率定结果进行验证,结果比较理想。最后简要讨论了参数不确定性传递。     

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.     

A geostatistical basis for spatial weighting in multivariate classification

Earth scientists and land managers often wish to group sampling sites that are both similar with respect to their properties and near to one another on the ground. This paper outlines the geostatistical rationale for such spatial grouping and describes a multivariate procedure to implement it. Sample variograms are calculated from the original data or their leading principal components and then the parameters of the underlying functions are estimated. A dissimilarity matrix is computed for all sampling sites, preferably using Gower's general similarity coefficient. Dissimilarities are then modified using the variogram to incorporate the form and extent of spatial variation. A nonhierarchical classification of sampling sites is performed on the leading latent vectors of the modified dissimilarity matrix by dynamic clustering to an optimum. The technique is illustrated with results of its application to soil survey data from two small areas in Britain and from a transect. In the case of the latter results of spatially weighted classifications are compared with those of strict segmentation. An appendix lists a Genstat program for a spatially constrained classification using a spherical variogram as an example.     

Improving the Ensemble Estimate of the Kalman Gain by Bootstrap Sampling

Using a small ensemble size in the ensemble Kalman filter methodology is efficient for updating numerical reservoir models but can result in poor updates following spurious correlations between observations and model variables. The most common approach for reducing the effect of spurious correlations on model updates is multiplication of the estimated covariance by a tapering function that eliminates all correlations beyond a prespecified distance. Distance-dependent tapering is not always appropriate, however. In this paper, we describe efficient methods for discriminating between the real and the spurious correlations in the Kalman gain matrix by using the bootstrap method to assess the confidence level of each element from the Kalman gain matrix. The new method is tested on a small linear problem, and on a water flooding reservoir history matching problem. For the water flooding example, a small ensemble size of 30 was used to compute the Kalman gain in both the screened EnKF and standard EnKF methods. The new method resulted in significantly smaller root mean squared errors of the estimated model parameters and greater variability in the final updated ensemble.     

Stationary and Isotropic Vector Random Fields on Spheres

This paper presents the characterization of the covariance matrix function of a Gaussian or second-order elliptically contoured vector random field on the sphere which is stationary, isotropic, and mean square continuous. This characterization involves an infinite sum of the products of positive definite matrices and Gegenbauer??s polynomials, and may not be available for other non-Gaussian vector random fields on spheres such as a ?? ² or log-Gaussian vector random field. We also offer two simple but efficient constructing approaches, and derive some parametric covariance matrix structures on spheres.     

A modification of minimum norm quadratic estimation of a generalized covariance function for use with large data sets

Marshall and Mardia (1985) and Kitanidis (1985) have suggested using minimum norm quadratic estimation as a method to estimate parameters of a generalized covariance function. Unfortunately, this method is difficult to use with large data sets as it requires inversion of an n × n matrix, where n is number of observations. These authors suggest replacing the matrix to be inverted by the identity matrix, which eliminates the computational burden, although with a considerable loss of efficiency. As an alternative, the data set can be broken into subsets, and minimum norm quadratic estimates of parameters of the generalized covariance function can be obtained within each subset. These local estimates can be averaged to obtain global estimates. This procedure also avoids large matrix inversions, but with less loss in efficiency.     

Kriging in a global neighborhood

The kriging estimator is usually computed in a moving neighborhood; only the data near the point to be estimated are used. This moving neighborhood approach creates discontinuities in mapping applications. An alternative approach is presented here, whereby all points are estimated using all the available data. To solve the resulting large linear system the kriging estimator is expressed in terms of the inverse of the covariance matrix. The covariance matrix has the advantage of being positive definite and the size of system which can be solved without encountering numerical instability is substantially increased. Because the kriging matrix does not change, the estimator can be written in terms of scalar products, thus avoiding the more time-consuming matrix multiplications of the standard approach. In the particular case of a covariance which is zero for distances greater than a fixed value (the range), the resulting banded structure of the covariance matrix is shown to lead to substantial computational savings in both run time and storage space. In this case the calculation time for the kriging variance is also substantially reduced. The present method is extended to the nonstationary case.     

Dating Quaternary fault scarps formed in alluvium using morphologic parameters

Ages of fault scarps, as well as those other types of transport-limited slopes, can be estimated by comparing their morphology with the morphology of scarps of known age. Age estimates are derived by fitting the scarp profiles to synthetic profiles generated using a diffusion equation or, alternatively, by classification using a linear discriminant function. The usefulness of morphology-derived age estimates depends on the relative importance of non-age-related morphologic variation. Data from more than 200 scarp profiles demonstrate that morphologic variation not related to scarp age can introduce significant uncertainties into morphology-derived age estimates.     

Estimation of a multidimensional covariance function in case of anisotropy

A model of a multivariate covariance function with an ellipsoidal directional correlation scale has been developed. The axes of the ellipsoidal scale are related to the eigenvalues and eigenvectors of a matrix B which characterizes the ellipsoid of the range of influence. The matrix B is found to be related to a matrix T which can be estimated directly from sparse sampling data and can be used to determine estimates of the matrix B. The method has been applied to both two-dimensional and three-dimensional cases. The numerical results show that the satisfactory accuracy is obtained with sparse sampling data from an anisotropic random function.     

Some properties of principal component scores

Commonly used methods for calculating component scores are reviewed. Means, variances, and the covariance structures of the resulting sets of scores are examined both by calculations based on a large set of electron microprobe analyses of melilite (supplied by D. Velde)and by a survey of recent geological applications of principal component analysis. Most of the procedures used to project raw data into the new vector space yield uncorrelated scores. In exceptions so far encountered, correlations between scores seem to have been occasioned by the use of unstandardized variables with components calculated from a correlation matrix. In a number of cases substantive interpretations of such correlations have been proposed. A different set of correlations results for the same data if scores are computed from standardized variables and components based on the covariance matrix. If unscaled components are rotated by the varimax procedure, the result is a return to the original space. In the work reported here, nevertheless, scores calculated from varimax-rotated scaled vectors are uncorrelated.     

Nonstationarity of the Mean and Unbiased Variogram Estimation: Extension of the Weighted Least-Squares Method

When concerned with spatial data, it is not unusual to observe a nonstationarity of the mean. This nonstationarity may be modeled through linear models and the fitting of variograms or covariance functions performed on residuals. Although it usually is accepted by authors that a bias is present if residuals are used, its importance is rarely assessed. In this paper, an expression of the variogram and the covariance function is developed to determine the expected bias. It is shown that the magnitude of the bias depends on the sampling configuration, the importance of the dependence between observations, the number of parameters used to model the mean, and the number of data. The applications of the expression are twofold. The first one is to evaluate a priori the importance of the bias which is expected when a residuals-based variogram model is used for a given configuration and a hypothetical data dependence. The second one is to extend the weighted least-squares method to fit the variogram and to obtain an unbiased estimate of the variogram. Two case studies show that the bias can be negligible or larger than 20%. The residual-based sample variogram underestimates the total variance of the process but the nugget variance may be overestimated.