Why do we need and how should we implement Bayesian kriging methods

The spatial prediction methodology that has become known under the heading of kriging is largely based on the assumptions that the underlying random field is Gaussian and the covariance function is exactly known. In practical applications, however, these assumptions will not hold. Beyond Gaussianity of the random field, lognormal kriging, disjunctive kriging, (generalized linear) model-based kriging and trans-Gaussian kriging have been proposed in the literature. The latter approach makes use of the Box–Cox-transform of the data. Still, all the alternatives mentioned do not take into account the uncertainty with respect to the distribution (or transformation) and the estimated covariance function of the data. The Bayesian trans-Gaussian kriging methodology proposed in the present paper is in the spirit of the “Bayesian bootstrap” idea advocated by Rubin (Ann Stat 9:130–134, 1981) and avoids the unusual specification of noninformative priors often made in the literature and is entirely based on the sample distribution of the estimators of the covariance function and of the Box–Cox parameter. After some notes on Bayesian spatial prediction, noninformative priors and developing our new methodology finally we will present an example illustrating our pragmatic approach to Bayesian prediction by means of a simulated data set.     

Wavelet analysis residual kriging vs. neural network residual kriging

 This paper deals with the problem of spatial data mapping. A new method based on wavelet interpolation and geostatistical prediction (kriging) is proposed. The method – wavelet analysis residual kriging (WARK) – is developed in order to assess the problems rising for highly variable data in presence of spatial trends. In these cases stationary prediction models have very limited application. Wavelet analysis is used to model large-scale structures and kriging of the remaining residuals focuses on small-scale peculiarities. WARK is able to model spatial pattern which features multiscale structure. In the present work WARK is applied to the rainfall data and the results of validation are compared with the ones obtained from neural network residual kriging (NNRK). NNRK is also a residual-based method, which uses artificial neural network to model large-scale non-linear trends. The comparison of the results demonstrates the high quality performance of WARK in predicting hot spots, reproducing global statistical characteristics of the distribution and spatial correlation structure.     

Geostatistical mapping of geomorphic variables in the presence of trend

Mapping geomorphic variables geostatistically, specifically by kriging, runs into difficulties when there is trend. The reason is that the variogram required for the kriging must be of residuals from any trend, which in turn cannot be estimated optimally by the usual method of trend surface analysis because the residuals are correlated. The difficulties can be overcome by the use of residual maximum likelihood (REML) to estimate both the trend and the variogram of the residuals simultaneously. We summarize the theory of REML as it applies to kriging in the presence of trend. We present the equations to show how estimates of the trend are combined with kriging of residuals to give empirical best linear unbiased predictions (E‐BLUPs). We then apply the method to estimate the height of the sub‐Upper‐Chalk surface beneath the Chiltern Hills of southeast England from 238 borehole data. The variogram of the REML residuals is substantially different from that computed by ordinary least squares (OLS) analysis. The map of the predicted surface is similar to that made from kriging with the OLS variogram. The variances, however, are substantially larger because (a) they derive from a variogram with a much larger sill and (b) they include the uncertainty of the estimate of the trend. Copyright © 2006 John Wiley & Sons, Ltd.     

Comparison of ordinary kriging and artificial neural network for spatial mapping of arsenic contamination of groundwater

In this technical note, we investigate the hypothesis that ‘non-linearity matters in the spatial mapping of complex patterns of groundwater arsenic contamination’. The spatial mapping pertained to data-driven techniques of spatial interpolation based on sampling data at finite locations. Using the well known example of extensive groundwater contamination by arsenic in Bangladesh, we find that the use of a highly non-linear pattern learning technique in the form of an artificial neural network (ANN) can yield more accurate results under the same set of constraints when compared to the ordinary kriging method. One ANN and a variogram model were used to represent the spatial structure of arsenic contamination for the whole country. The probability for successful detection of a well as safe or unsafe was found to be atleast 15% larger than that by kriging under the country-wide scenario. The probability of false hopes, which is a serious issue in public health monitoring was found to be significantly lower (by more than 10%) than that by kriging.     

Simultaneous estimation of the parameters of the Hurst�CKolmogorov stochastic process

Various methods for estimating the self-similarity parameter (Hurst parameter, H) of a Hurst–Kolmogorov stochastic process (HKp) from a time series are available. Most of them rely on some asymptotic properties of processes with Hurst–Kolmogorov behaviour and only estimate the self-similarity parameter. Here we show that the estimation of the Hurst parameter affects the estimation of the standard deviation, a fact that was not given appropriate attention in the literature. We propose the least squares based on variance estimator, and we investigate numerically its performance, which we compare to the least squares based on standard deviation estimator, as well as the maximum likelihood estimator after appropriate streamlining of the latter. These three estimators rely on the structure of the HKp and estimate simultaneously its Hurst parameter and standard deviation. In addition, we test the performance of the three methods for a range of sample sizes and H values, through a simulation study and we compare it with other estimators of the literature.     

Multigaussian kriging for point-support estimation: incorporating constraints on the sum of the kriging weights

In the geostatistical analysis of regionalized data, the practitioner may not be interested in mapping the unsampled values of the variable that has been monitored, but in assessing the risk that these values exceed or fall short of a regulatory threshold. This kind of concern is part of the more general problem of estimating a transfer function of the variable under study. In this paper, we focus on the multigaussian model, for which the regionalized variable can be represented (up to a nonlinear transformation) by a Gaussian random field. Two cases are analyzed, depending on whether the mean of this Gaussian field is considered known or not, which lead to the simple and ordinary multigaussian kriging estimators respectively. Although both of these estimators are theoretically unbiased, the latter may be preferred to the former for practical applications since it is robust to a misspecification of the mean value over the domain of interest and also to local fluctuations around this mean value. An advantage of multigaussian kriging over other nonlinear geostatistical methods such as indicator and disjunctive kriging is that it makes use of the multivariate distribution of the available data and does not produce order relation violations. The use of expansions into Hermite polynomials provides three additional results: first, an expression of the multigaussian kriging estimators in terms of series that can be calculated without numerical integration; second, an expression of the associated estimation variances; third, the derivation of a disjunctive-type estimator that minimizes the variance of the error when the mean is unknown.     

Links, comparisons and extensions of the geographically weighted regression model when used as a spatial predictor

In this study, we link and compare the geographically weighted regression (GWR) model with the kriging with an external drift (KED) model of geostatistics. This includes empirical work where models are performance tested with respect to prediction and prediction uncertainty accuracy. In basic forms, GWR and KED (specified with local neighbourhoods) both cater for nonstationary correlations (i.e. the process is heteroskedastic with respect to relationships between the variable of interest and its covariates) and as such, can predict more accurately than models that do not. Furthermore, on specification of an additional heteroskedastic term to the same models (now with respect to a process variance), locally-accurate measures of prediction uncertainty can result. These heteroskedastic extensions of GWR and KED can be preferred to basic constructions, whose measures of prediction uncertainty are only ever likely to be globally-accurate. We evaluate both basic and heteroskedastic GWR and KED models using a case study data set, where data relationships are known to vary across space. Here GWR performs well with respect to the more involved KED model and as such, GWR is considered a viable alternative to the more established model in this particular comparison. Our study adds to a growing body of empirical evidence that GWR can be a worthy predictor; complementing its more usual guise as an exploratory technique for investigating relationships in multivariate spatial data sets.     

Gaussian process classification: singly versus doubly stochastic models, and new computational schemes

The aim of this paper is to compare four different methods for binary classification with an underlying Gaussian process with respect to theoretical consistency and practical performance. Two of the inference schemes, namely classical indicator kriging and simplicial indicator kriging, are analytically tractable and fast. However, these methods rely on simplifying assumptions which are inappropriate for categorical class labels. A consistent and previously described model extension involves a doubly stochastic process. There, the unknown posterior class probability f(·) is considered a realization of a spatially correlated Gaussian process that has been squashed to the unit interval, and a label at position x is considered an independent Bernoulli realization with success parameter f(x). Unfortunately, inference for this model is not known to be analytically tractable. In this paper, we propose two new computational schemes for the inference in this doubly stochastic model, namely the “Aitchison Maximum Posterior” and the “Doubly Stochastic Gaussian Quadrature”. Both methods are analytical up to a final step where optimization or integration must be carried out numerically. For the comparison of practical performance, the methods are applied to storm forecasts for the Spanish coast based on wave heights in the Mediterranean Sea. While the error rate of the doubly stochastic models is slightly lower, their computational cost is much higher.     

Declustering experimental variograms by global estimation with fourth order moments

The variogram is a key parameter for geostatistical estimation and simulation. Preferential sampling may bias the spatial structure and often leads to noisy and unreliable variograms. A novel technique is proposed to weight variogram pairs in order to compensate for preferential or clustered sampling . Weighting the variogram pairs by global kriging of the quadratic differences between the tail and head values gives each pair the appropriate weight, removes noise and minimizes artifacts in the experimental variogram. Moreover, variogram uncertainty could be computed by this technique. The required covariance between the pairs going into variogram calculation, is a fourth order covariance that must be calculated by second order moments. This introduces some circularity in the calculation whereby an initial variogram must be assumed before calculating how the pairs should be weighted for the experimental variogram. The methodology is assessed by synthetic and realistic examples. For synthetic example, a comparison between the traditional and declustered variograms shows that the declustered variograms are better estimates of the true underlying variograms. The realistic example also shows that the declustered sample variogram is closer to the true variogram.     

Spatial prediction on river networks: comparison of top‐kriging with regional regression

Top‐kriging is a method for estimating stream flow‐related variables on a river network. Top‐kriging treats these variables as emerging from a two‐dimensional spatially continuous process in the landscape. The top‐kriging weights are estimated by regularising the point variogram over the catchment area (kriging support), which accounts for the nested nature of the catchments. We test the top‐kriging method for a comprehensive Austrian data set of low stream flows. We compare it with the regional regression approach where linear regression models between low stream flow and catchment characteristics are fitted independently for sub‐regions of the study area that are deemed to be homogeneous in terms of flow processes. Leave‐one‐out cross‐validation results indicate that top‐kriging outperforms the regional regression on average over the entire study domain. The coefficients of determination (cross‐validation) of specific low stream flows are 0.75 and 0.68 for the top‐kriging and regional regression methods, respectively. For locations without upstream data points, the performances of the two methods are similar. For locations with upstream data points, top‐kriging performs much better than regional regression as it exploits the low flow information of the neighbouring locations. Copyright © 2012 John Wiley & Sons, Ltd.     

Spatial prediction of river channel topography by kriging

Topographic information is fundamental to geomorphic inquiry, and spatial prediction of bed elevation from irregular survey data is an important component of many reach‐scale studies. Kriging is a geostatistical technique for obtaining these predictions along with measures of their reliability, and this paper outlines a specialized framework intended for application to river channels. Our modular approach includes an algorithm for transforming the coordinates of data and prediction locations to a channel‐centered coordinate system, several different methods of representing the trend component of topographic variation and search strategies that incorporate geomorphic information to determine which survey data are used to make a prediction at a specific location. For example, a relationship between curvature and the lateral position of maximum depth can be used to include cross‐sectional asymmetry in a two‐dimensional trend surface model, and topographic breaklines can be used to restrict which data are retained in a local neighborhood around each prediction location. Using survey data from a restored gravel‐bed river, we demonstrate how transformation to the channel‐centered coordinate system facilitates interpretation of the variogram, a statistical model of reach‐scale spatial structure used in kriging, and how the choice of a trend model affects the variogram of the residuals from that trend. Similarly, we show how decomposing kriging predictions into their trend and residual components can yield useful information on channel morphology. Cross‐validation analyses involving different data configurations and kriging variants indicate that kriging is quite robust and that survey density is the primary control on the accuracy of bed elevation predictions. The root mean‐square error of these predictions is directly proportional to the spacing between surveyed cross‐sections, even in a reconfigured channel with a relatively simple morphology; sophisticated methods of spatial prediction are no substitute for field data. Copyright © 2007 John Wiley & Sons, Ltd.     

Probability density estimation in stochastic environmental models using reverse representations

The estimation of probability densities of variables described by stochastic differential equations has long been done using forward time estimators, which rely on the generation of forward in time realizations of the model. Recently, an estimator based on the combination of forward and reverse time estimators has been developed. This estimator has a higher order of convergence than the classical one. In this article, we explore the new estimator and compare the forward and forward–reverse estimators by applying them to a biochemical oxygen demand model. Finally, we show that the computational efficiency of the forward–reverse estimator is superior to the classical one, and discuss the algorithmic aspects of the estimator.     

Testing the correctness of the sequential algorithm for simulating Gaussian random fields

The sequential algorithm is widely used to simulate Gaussian random fields. However, a rigorous application of this algorithm is impractical and some simplifications are required, in particular a moving neighborhood has to be defined. To examine the effect of such restriction on the quality of the realizations, a reference case is presented and several parameters are reviewed, mainly the histogram, variogram, indicator variograms, as well as the ergodic fluctuations in the first and second-order statistics. The study concludes that, even in a favorable case where the simulated domain is large with respect to the range of the model, the realizations may poorly reproduce the second-order statistics and be inconsistent with the stationarity and ergodicity assumptions. Practical tips such as the multiple-grid strategy do not overcome these impediments. Finally, extending the original algorithm by using an ordinary kriging should be avoided, unless an intrinsic random function model is sought after.     

A geostatistical re-interpretation of gravity surveys in the Yagoua,Cameroon region

Since 1960, many gravity studies have been carried out in the Yagoua region of northern Cameroon. Gravity data was collected over a wide area of approximately 11628 km². These data are insufficient, irregular, scattered and do not efficiently permit gravity field downward and upward continuations, derivatives and other operations that might require regular gridded data. Some anomalies on the Collignon map (1968), may correlate with known geological structure but do not appear on maps by Louis (1970) and Poudjom et al. (1996). To produce regular gridded gravity data and better control anomalies due to geological structures, the kriging method was applied to a 188-data baseline. Several variogram models were tested for this purpose. It was found that a spherical variogram model is the best; it has produced a new kriging dataset of about 10,100 data and a new map of kriged Bouguer data. This map contains positive anomalies in the Maroua-Mindif and Maga areas on the Collignon (1968) map, which were not present on Louis (1970) and Poudjom et al. (1996) maps. The positive anomalies of Guibi-Doukoula and Yagoua, not separated on the Louis (1970) and Poudjom et al. (1996) maps, show up as clearly distinct as previewed by Collignon (1968). The new results can be used for subsequent gravimetric studies.     

Space–time forecasting of groundwater level using a hybrid soft computing model

Forecasting of space–time groundwater level is important for sparsely monitored regions. Time series analysis using soft computing tools is powerful in temporal data analysis. Classical geostatistical methods provide the best estimates of spatial data. In the present work a hybrid framework for space–time groundwater level forecasting is proposed by combining a soft computing tool and a geostatistical model. Three time series forecasting models: artificial neural network, least square support vector machine and genetic programming (GP), are individually combined with the geostatistical ordinary kriging model. The experimental variogram thus obtained fits a linear combination of a nugget effect model and a power model. The efficacy of the space–time models was decided on both visual interpretation (spatial maps) and calculated error statistics. It was found that the GP–kriging space–time model gave the most satisfactory results in terms of average absolute relative error, root mean square error, normalized mean bias error and normalized root mean square error.     

Exploring a valid model for the variogram of an isotropic spatial process

The variogram is one of the most important tools in the assessment of spatial variability and a crucial parameter for kriging. It is widely known that an estimator for the variogram cannot be used as its representator in some contexts because of its lack of conditional semi negative definiteness. Consequently, once the variogram is estimated, a valid family must be chosen to fit an appropriate model. Under isotropy, this selection is carried out by eye from the observation of the variogram estimated curve. In this paper, a statistical methodology is proposed to explore a valid model for the variogram. The statistic for this approach is based on quadratic forms depending on smoothed random variables which gather the underlying spatial variation. The distribution of the test statistic is approximated by a shifted chi-square distribution. A simulation study is also carried out to check the power and size of the test. Reference bands, as a complementary graphical tool, are calculated. An example from the literature is used to illustrate the methodologies presented.     

Application of modified Patankar schemes to stiff biogeochemical models for the water column

In this paper, we apply recently developed positivity preserving and conservative Modified Patankar-type solvers for ordinary differential equations to a simple stiff biogeochemical model for the water column. The performance of this scheme is compared to schemes which are not unconditionally positivity preserving (the first-order Euler and the second- and fourth-order Runge–Kutta schemes) and to schemes which are not conservative (the first- and second-order Patankar schemes). The biogeochemical model chosen as a test ground is a standard nutrient–phytoplankton–zooplankton–detritus (NPZD) model, which has been made stiff by substantially decreasing the half saturation concentration for nutrients. For evaluating the stiffness of the biogeochemical model, so-called numerical time scales are defined which are obtained empirically by applying high-resolution numerical schemes. For all ODE solvers under investigation, the temporal error is analysed for a simple exponential decay law. The performance of all schemes is compared to a high-resolution high-order reference solution. As a result, the second-order modified Patankar–Runge–Kutta scheme gives a good agreement with the reference solution even for time steps 10 times longer than the shortest numerical time scale of the problem. Other schemes do either compute negative values for non-negative state variables (fully explicit schemes), violate conservation (the Patankar schemes) or show low accuracy (all first-order schemes).     

Geostatistical radar–raingauge merging: A novel method for the quantification of rain estimation accuracy

Compared to other estimation techniques, one advantage of geostatistical techniques is that they provide an index of the estimation accuracy of the variable of interest with the kriging estimation standard deviation (ESD). In the context of radar–raingauge quantitative precipitation estimation (QPE), we address in this article the question of how the kriging ESD can be transformed into a local spread of error by using the dependency of radar errors to the rain amount analyzed in previous work. The proposed approach is implemented for the most significant rain events observed in 2008 in the Cévennes-Vivarais region, France, by considering both the kriging with external drift (KED) and the ordinary kriging (OK) methods. A two-step procedure is implemented for estimating the rain estimation accuracy: (i) first kriging normalized ESDs are computed by using normalized variograms (sill equal to 1) to account for the observation system configuration and the spatial structure of the variable of interest (rainfall amount, residuals to the drift); (ii) based on the assumption of a linear relationship between the standard deviation and the mean of the variable of interest, a denormalization of the kriging ESDs is performed globally for a given rain event by using a cross-validation procedure. Despite the fact that the KED normalized ESDs are usually greater than the OK ones (due to an additional constraint in the kriging system and a weaker spatial structure of the residuals to the drift), the KED denormalized ESDs are generally smaller the OK ones, a result consistent with the better performance observed for the KED technique. The evolution of the mean and the standard deviation of the rainfall-scaled ESDs over a range of spatial (5–300 km²) and temporal (1–6 h) scales demonstrates that there is clear added value of the radar with respect to the raingauge network for the shortest scales, which are those of interest for flash-flood prediction in the considered region.     

Comparison between the indirect approach and kriging with samples of different support for estimation using samples of different length

In this paper we compare two estimation methods to deal with samples of different support: (1) the indirect approach using accumulation and (2) kriging with samples of different support. These two methods were tested in a simple example. The estimates of the two methods were compared against a benchmark scenario. The benchmark consisted of kriging using a complete set of samples on the same support. The effects of the nugget effect, variogram range and type on the weight of long samples, the estimate, and the error variance were assessed. Kriging with samples of different support led to lower error variance and to estimates closer to the estimates of the benchmark scenario. Furthermore, in the case of spatially continuous attributes (low nugget effect), the indirect approach assigns greater weight to long samples than kriging with samples of different support. A cross validation study comparing the two methods with a database from a bauxite deposit was performed. The results of the cross validation study showed that kriging with samples of different support resulted in more precise estimates.     

Transition probability/Markov chain analyses of DNAPL source zones and plumes

At sites where a dense nonaqueous phase liquid (DNAPL) was spilled or released into the subsurface, estimates of the mass of DNAPL contained in the subsurface from core or monitoring well data, either in the nonaqueous or aqueous phase, can be highly uncertain because of the erratic distribution of the DNAPL due to geologic heterogeneity. In this paper, a multiphase compositional model is applied to simulate, in detail, the DNAPL saturations and aqueous-phase plume migration in a highly characterized, heterogeneous glaciofluvial aquifer, the permeability and porosity data of which were collected by researchers at the University of Tübingen, Germany. The DNAPL saturation distribution and the aqueous-phase contaminant mole fractions are then reconstructed by sampling the data from the forward simulation results using two alternate approaches, each with different degrees of sampling conditioning. To reconstruct the DNAPL source zone architecture, the aqueous-phase plume configuration, and the contaminant mass in each phase, one method employs the novel transition probability/Markov chain approach (TP/MC), while the other involves a traditional variogram analysis of the sampled data followed by ordinary kriging. The TP/MC method is typically used for facies and/or hydraulic conductivity reconstruction, but here we explore the applicability of the TP/MC method for the reconstruction of DNAPL source zones and aqueous-phase plumes. The reconstructed geometry of the DNAPL source zone, the dissolved contaminant plume, and the estimated mass in each phase are compared using the two different geostatistical modeling approaches and for various degrees of data sampling from the results of the forward simulation. It is demonstrated that the TP/MC modeling technique is robust and accurate and is a preferable alternative compared to ordinary kriging for the reconstruction of DNAPL saturation patterns and dissolved-phase contaminant plumes.