A non-stationary geostatistical approach to multigaussian kriging for local reserve estimation

Multigaussian kriging technique has many applications in mining, soil science, environmental science and other fields. Particularly, in the local reserve estimation of a mineral deposit, multigaussian kriging is employed to derive panel-wise tonnages by predicting conditional probability of block grades. Additionally, integration of a suitable change of support model is also required to estimate the functions of the variables with larger support than that of the samples. However, under the assumption of strict stationarity, the grade distributions and important recovery functions are estimated by multigaussian kriging using samples within a supposedly spatial homogeneous domain. Conventionally, the underlying random function model is required to be stationary in order to carry out the inference on ore grade distribution and relevant statistics. In reality, conventional stationary model often fails to represent complicated geological structure. Traditionally, the simple stationary model neither considers the obvious changes in local means and variances, nor is it able to replicate spatial continuity of the deposit and hence produces unreliable outcomes. This study deals with the theoretical design of a non-stationary multigaussian kriging model allowing change of support and its application in the mineral reserve estimation scenario. Local multivariate distributions are assumed here to be strictly stationary in the neighborhood of the panels. The local cumulative distribution function and related statistics with respect to the panels are estimated using a distance kernel approach. A rigorous investigation through simulation experiments is performed to analyze the relevance of the developed model followed by a case study on a copper deposit.     

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.     

Regional flood frequency analysis using residual kriging in physiographical space

In this article, an approach using residual kriging (RK) in physiographical space is proposed for regional flood frequency analysis. The physiographical space is constructed using physiographical/climatic characteristics of gauging basins by means of canonical correlation analysis (CCA). This approach is a modified version of the original method, based on ordinary kriging (OK). It is intended to handle effectively any possible spatial trends within the hydrological variables over the physiographical space. In this approach, the trend is first quantified and removed from the hydrological variable by a quadratic spatial regression. OK is therefore applied to the regression residual values. The final estimated value of a specific quantile at an ungauged station is the sum of the spatial regression estimate and the kriged residual. To evaluate the performance of the proposed method, a cross‐validation procedure is applied. Results of the proposed method indicate that RK in CCA physiographical space leads to more efficient estimates of regional flood quantiles when compared to the original approach and to a straightforward regression‐based estimator. Copyright © 2010 John Wiley & Sons, Ltd.     

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.     

Broadband Kinematic Stochastic Simulation of an Earthquake Source: a Refined Procedure for Application in Seismic Hazard Studies

To carry out a realistic simulation of earthquake strong ground motion for applied studies, one needs an earthquake fault/source simulator that can integrate most relevant features of observed earthquake ruptures. A procedure of this kind is proposed that creates a broadband kinematic source model. At lower frequencies, the source is described as propagating slip pulse with locally variable velocity. The final slip is assumed to be a two-dimensional (2D) random function. At higher frequencies, radiation from the same running strip is assumed to be random and incoherent in space. The model is discretized in space as a grid of point subsources with certain time histories. At lower frequencies, a realistic shape of source spectrum is generated implicitly by simulated kinematics of slip pulse propagation. At higher frequencies, the original approach is used to generate signals with spectra that plausibly approximate the prescribed smooth far-field source spectrum. This spectrum is set on the basis of the assumedly known regional empirical spectral scaling law, and subsource moment rate time histories are conditioned so as to fit this expected spectrum. For the random function that describes final slip over the fault area, lognormal probability distribution of amplitudes is assumed, on the basis of exploratory analysis of inverted slip distributions. Similarly, random functions that describe local slip rate time histories are assumed to have lognormal distribution of envelope amplitudes. In this way one can effectively emulate expressed ??asperities?? of final slip and occasional occurrence of large spikes on near-source accelerograms. A special procedure is proposed to simulate the spatial coherence of high-frequency fault motion. This approach permits the simulation of fault motion plausibly at high spatial resolution, fulfilling the prerequisite for simulation of strong motion in the vicinity of a fault. A particular realization (sample) of a source created in a simulation run depends on several random seeds, and also on a considerable number of parameters. Their values can be selected so as to take into account expected source features; they can also be perturbed to examine the source-related component of uncertainty of strong motion. The proposed approach to earthquake source specification is well adapted to the study of deterministic seismic hazard: it may be used for simulation of individual scenario events, or suites of such events, as well as for analysis of uncertainty for expected ground motion parameters from a particular class of events. Examples are given of application of the proposed approach to strong motion simulations and related uncertainty estimation.     

Multiscale flow and transport model in three-dimensional fractal porous media

This paper proposes a multiscale flow and transport model which can be used in three-dimensional fractal random fields. The fractal random field effectively describes a field with a high degree of variability to satisfy the one-point statistics of Levy-stable distribution and the two-point statistics of fractional Levy motion (fLm). To overcome the difficulty of using infinite variance of Levy-stable distribution and to provide the physical meaning of a finite domain in real space, truncated power variograms are utilized for the fLm fields. The fLm model is general in the sense that both stationary and commonly used fractional Brownian motion (fBm) models are its special cases. When the upper cutoff of the truncated power variogram is close to the lower cutoff, the stationary model is well approximated. The commonly used fBm model is recovered when the Levy index of fLm is 2. Flow and solute transport were analyzed using the first-order perturbation method. Mean velocity, velocity covariance, and effective hydraulic conductivity in a three-dimensional fractal random field were derived. Analytical results for particle displacement covariance and macrodispersion coefficients are also presented. The results show that the plume in an fLm field moves slower at early time and has more significant long-tailing behavior at late time than in fBm or stationary exponential fields. The proposed fractal transport model has broader applications than those of stationary and fBm models. Flow and solute transport can be simulated for various scenarios by adjusting the Levy index and cutoffs of fLm to yield more accurate modeling results.     

Non-stationary response of structures excited by random seismic processes with time variable frequency content

Seismic random processes are characterized by high non-stationarity and, in most cases, by a marked variability of frequency content. The hypothesis modeling seismic signal as a simple product of a stationary signal and a deterministic modulation function, consequently, is hardly ever applicable. Several mathematical models aimed at expressing the recorded process by means of a system of stationary random processes and deterministic amplitude and frequency modulations are proposed. Models oriented into the frequency domain with subsequent response analysis based on integral spectral resolution and models oriented into the time domain based on the multicomponent resolution are investigated. The resolution into individual components (non-stationary signals) is carried out by three methods. The resolution into intrinsic mode functions seems to possess the best characteristics and yields results almost not differing from the results obtained by stochastic simulation. An example of the seismic response of an existing bridge obtained by two older models and three variants of multicomponent resolution is given.     

Seismic response of a cable‐stayed bridge deck under multi‐component non‐stationary random ground motion

A Markov method of analysis is presented for obtaining the seismic response of cable‐stayed bridges to non‐stationary random ground motion. A uniformly modulated non‐stationary model of the random ground motion is assumed which is specified by the evolutionary r.m.s. ground acceleration. Both vertical and horizontal components of the motion are considered to act simultaneously at the bridge supports. The analysis duly takes into account the angle of incidence of the earthquake, the spatial correlation of ground motion and the quasi‐static excitation. A cable‐stayed bridge is analysed under a set of parametric variations in order to study the non‐stationary response of the bridge. The results of the numerical study indicate that (i) frequency domain spectral analysis with peak r.m.s. acceleration as input could provide more r.m.s. response than the peak r.m.s. response obtained by the non‐stationary analysis; (ii) the longitudinal component of the ground motion significantly influences the vertical vibration of the bridge; and (iii) the angle of incidence of the earthquake has considerable influence on the deck response. Copyright © 2003 John Wiley & Sons, Ltd.     

Simulated precipitation fields with variance-consistent interpolation

Abstract

Gridded meteorological data are available for all of Norway as time series dating from 1961. A new way of interpolating precipitation in space from observed values is proposed. Based on the criteria that interpolated precipitation fields in space should be consistent with observed spatial statistics, such as spatial mean, variance and intermittency, spatial fields of precipitation are simulated from a gamma distribution with parameters determined from observed data, adjusted for intermittency. The simulated data are distributed in space, using the spatial pattern derived from kriging. The proposed method is compared to indicator kriging and to the current methodology used for producing gridded precipitation data. Cross-validation gave similar results for the three methods with respect to RMSE, temporal mean and standard deviation, whereas a comparison on estimated spatial variance showed that the new method has a near perfect agreement with observations. Indicator kriging underestimated the spatial variance by 60–80% and the current method produced a significant scatter in its estimates.

Citation Skaugen, T. & Andersen, J. (2010) Simulated precipitation fields with variance-consistent interpolation. Hydrol. Sci. J. 55(5), 676–686.     

Cokriging for enhanced spatial interpolation of rainfall in two Australian catchments

Rainfall data in continuous space provide an essential input for most hydrological and water resources planning studies. Spatial distribution of rainfall is usually estimated using ground‐based point rainfall data from sparsely positioned rain‐gauge stations in a rain‐gauge network. Kriging has become a widely used interpolation method to estimate the spatial distribution of climate variables including rainfall. The objective of this study is to evaluate three geostatistical (ordinary kriging [OK], ordinary cokriging [OCK], kriging with an external drift [KED]), and two deterministic (inverse distance weighting, radial basis function) interpolation methods for enhanced spatial interpolation of monthly rainfall in the Middle Yarra River catchment and the Ovens River catchment in Victoria, Australia. Historical rainfall records from existing rain‐gauge stations of the catchments during 1980–2012 period are used for the analysis. A digital elevation model of each catchment is used as the supplementary information in addition to rainfall for the OCK and kriging with an external drift methods. The prediction performance of the adopted interpolation methods is assessed through cross‐validation. Results indicate that the geostatistical methods outperform the deterministic methods for spatial interpolation of rainfall. Results also indicate that among the geostatistical methods, the OCK method is found to be the best interpolator for estimating spatial rainfall distribution in both the catchments with the lowest prediction error between the observed and estimated monthly rainfall. Thus, this study demonstrates that the use of elevation as an auxiliary variable in addition to rainfall data in the geostatistical framework can significantly enhance the estimation of rainfall over a catchment.     

Spatial sampling design and covariance-robust minimax prediction based on convex design ideas

This paper presents new ideas on sampling design and minimax prediction in a geostatistical model setting. Both presented methodologies are based on regression design ideas. For this reason the appendix of this paper gives an introduction to optimum Bayesian experimental design theory for linear regression models with uncorrelated errors. The presented methodologies and algorithms are then applied to the spatial setting of correlated random fields. To be specific, in Sect. 1 we will approximate an isotropic random field by means of a regression model with a large number of regression functions with random amplitudes, similarly to Fedorov and Flanagan (J Combat Inf Syst Sci: 23, 1997). These authors make use of the Karhunen Loeve approximation of the isotropic random field. We use the so-called polar spectral approximation instead; i.e. we approximate the isotropic random field by means of a regression model with sine-cosine-Bessel surface harmonics with random amplitudes and then, in accordance with Fedorov and Flanagan (J Combat Inf Syst Sci: 23, 1997), apply standard Bayesian experimental design algorithms to the resulting Bayesian regression model. Section 2 deals with minimax prediction when the covariance function is known to vary in some set of a priori plausible covariance functions. Using a minimax theorem due to Sion (Pac J Math 8:171–176, 1958) we are able to formulate the minimax problem as being equivalent to an optimum experimental design problem, too. This makes the whole experimental design apparatus available for finding minimax kriging predictors. Furthermore some hints are given, how the approach to spatial sampling design with one a priori fixed covariance function may be extended by means of minimax kriging to a whole set of a priori plausible covariance functions such that the resulting designs are robust. The theoretical developments are illustrated with two examples taken from radiological monitoring and soil science.     

Compositional Bayesian indicator estimation

Indicator kriging is widely used for mapping spatial binary variables and for estimating the global and local spatial distributions of variables in geosciences. For continuous random variables, indicator kriging gives an estimate of the cumulative distribution function, for a given threshold, which is then the estimate of a probability. Like any other kriging procedure, indicator kriging provides an estimation variance that, although not often used in applications, should be taken into account as it assesses the uncertainty of the estimate. An alternative approach to indicator estimation is proposed in this paper. In this alternative approach the complete probability density function of the indicator estimate is evaluated. The procedure is described in a Bayesian framework, using a multivariate Gaussian likelihood and an a priori distribution which are both combined according to Bayes theorem in order to obtain a posterior distribution for the indicator estimate. From this posterior distribution, point estimates, interval estimates and uncertainty measures can be obtained. Among the point estimates, the median of the posterior distribution is the maximum entropy estimate because there is a fifty-fifty chance of the unknown value of the estimate being larger or smaller than the median; that is, there is maximum uncertainty in the choice between two alternatives. Thus in some sense, the latter is an indicator estimator, alternative to the kriging estimator, that includes its own uncertainty. On the other hand, the mode of the posterior distribution estimator, assuming a uniform prior, is coincidental with the simple kriging estimator. Additionally, because the indicator estimate can be considered as a two-part composition which domain of definition is the simplex, the method is extended to compositional Bayesian indicator estimation. Bayesian indicator estimation and compositional Bayesian indicator estimation are illustrated with an environmental case study in which the probability of the content of a geochemical element in soil being over a particular threshold is of interest. The computer codes and its user guides are public domain and freely available.     

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.     

Effective dispersion in a chemically heterogeneous medium under temporally fluctuating flow conditions

We investigate effective solute transport in a chemically heterogeneous medium subject to temporal fluctuations of the flow conditions. Focusing on spatial variations in the equilibrium adsorption properties, the corresponding fluctuating retardation factor is modeled as a stationary random space function. The temporal variability of the flow is represented by a stationary temporal random process. Solute spreading is quantified by effective dispersion coefficients, which are derived from the ensemble average of the second centered moments of the normalized solute distribution in a single disorder realization. Using first-order expansions in the variances of the respective random fields, we derive explicit compact expressions for the time behavior of the disorder induced contributions to the effective dispersion coefficients. Focusing on the contributions due to chemical heterogeneity and temporal fluctuations, we find enhanced transverse spreading characterized by a transverse effective dispersion coefficient that, in contrast to transport in steady flow fields, evolves to a disorder-induced macroscopic value (i.e., independent of local dispersion). At the same time, the asymptotic longitudinal dispersion coefficient can decrease. Under certain conditions the contribution to the longitudinal effective dispersion coefficient shows superdiffusive behavior, similar to that observed for transport in s stratified porous medium, before it decreases to its asymptotic value. The presented compact and easy to use expressions for the longitudinal and transverse effective dispersion coefficients can be used for the quantification of effective spreading and mixing in the context of the groundwater remediation based on hydraulic manipulation and for the effective modeling of reactive transport in heterogeneous media in general.     

水文时间序列概率预报方法的通用架构

在准确揭示水文过程变化特性的基础上开展中长期(月尺度及以上)水文预报,是掌握未来水文情势和演变规律,以及研究解决实际水文水资源问题的重要基础.水文时间序列预报方法是揭示未来水文情势和演变规律的重要技术手段.本文首先梳理了目前常用的各类水文序列预报方法,分析讨论了各方法的基本原理和主要缺陷.然后,通过综合分析相关研究成果,总结得到关于水文序列预报方法的4点重要认识:序列预报前应进行序列分解;序列中确定成分和随机成分应分别建模预报;序列预报结果需要估计不确定性;模型集成效果常常优于单个模型效果.最后,提出一个水文时间序列概率预报方法的通用架构.利用该通用架构能够克服常规模型或方法的缺陷,进行物理成因分析的基础上,针对水文序列中不同特性的确定成分和随机成分别进行分析,既可得到准确的确定性预报结果,又可对预报结果的不确定性进行定量评估,并可提高最终预报结果的合理性和可靠性.     

Detecting and quantifying sources of non-stationarity via experimental semivariogram modeling

Conventional geostatistics often relies on the assumption of second order stationarity of the random function (RF). Generally, local means and local variances of the random variables (RVs) are assumed to be constant throughout the domain. Large scale differences in the local means and local variances of the RVs are referred to as trends. Two problems of building geostatistical models in presence of mean trends are: (1) inflation of the conditional variances and (2) the spatial continuity is exaggerated. Variance trends on the other hand cause conditional variances to be over-estimated in certain regions of the domain and under-estimated in other areas. In both cases the uncertainty characterized by the geostatistical model is improperly assessed. This paper proposes a new approach to identify the presence and contribution of mean and variance trends in the domain via calculation of the experimental semivariogram. The traditional experimental semivariogram expression is decomposed into three components: (1) the mean trend, (2) the variance trend and (3) the stationary component. Under stationary conditions, both the mean and the variance trend components should be close to zero. This proposed approach is intended to be used in the early stages of data analysis when domains are being defined or to verify the impact of detrending techniques in the conditioning dataset for validating domains. This approach determines the source of a trend, thereby facilitating the choice of a suitable detrending method for effective resource modeling.     

The generation of spectrum compatible accelerograms for the design of nuclear power plants

Based on the main features of the computer program THGE, the paper reviews the techniques which are available for the construction of an accelerogram whose response spectrum matches a design spectrum. First, a sample accelerogram is generated as the product of the stationary random sequence by a deterministic shape function. It is assumed that the spectral properties of the stationary process have been determined in a previous step, in order to lead to the design spectrum as expected maximum responses of a set of single degree of freedom oscillators. The principal properties of the Fast Fourier Transform which is used for this first step are reviewed. The paper then describes the procedures which are available, both in the frequency domain and in the time domain, to improve the agreement between the response spectrum and the target. It also discusses some related issues, such as the response spectrum calculations, the statistical dependence between the three earthquake components, the duration of the time history, the variability of the secondary response to various samples and the generation of an accelerogram whose response spectra envelop a set of design spectra. The point of view adopted is the one of the structural engineer.     

Conditional bias-penalized kriging (CBPK)

Simple and ordinary kriging, or SK and OK, respectively, represent the best linear unbiased estimator in the unconditional sense in that they minimize the unconditional (on the unknown truth) error variance and are unbiased in the unconditional mean. However, because the above properties hold only in the unconditional sense, kriging estimates are generally subject to conditional biases that, depending on the application, may be unacceptably large. For example, when used for precipitation estimation using rain gauge data, kriging tends to significantly underestimate large precipitation and, albeit less consequentially, overestimate small precipitation. In this work, we describe an extremely simple extension to SK or OK, referred to herein as conditional bias-penalized kriging (CBPK), which minimizes conditional bias in addition to unconditional error variance. For comparative evaluation of CBPK, we carried out numerical experiments in which normal and lognormal random fields of varying spatial correlation scale and rain gauge network density are synthetically generated, and the kriging estimates are cross-validated. For generalization and potential application in other optimal estimation techniques, we also derive CBPK in the framework of classical optimal linear estimation theory.