An efficient multivariate random field generator using the fast Fourier transform

It is often convenient to use synthetically generated random fields to study the hydrologic effects of spatial heterogeneity. Although there are many ways to produce such fields, spectral techniques are particularly attractive because they are fast and conceptually straightforward. This paper describes a spectral algorithm for generating sets of random fields which are correlated with one another. The algorithm is based on a discrete version of the Fourier-Stieltjes representation for multidimensional random fields. The Fourier increment used in this representation depends on a random phase angle process and a complex-valued spectral factor matrix which can be readily derived from a specified set of cross-spectral densities (or cross-covariances). The inverse Fourier transform of the Fourier increment is a complex random field with real and imaginary parts which each have the desired coveriance structure. Our complex-valued spectral formulation provides an especially convenient way to generate a set of random fields which all depend on a single underlying (independent) field, provided that the fields in question can be related by space-invariant linear transformations. We illustrate this by generating multi-dimensional mass conservative groundwater velocity fields which can be used to simulate solute transport through heterogeneous anisotropic porous media.     

On a continuous spectral algorithm for simulating non-stationary Gaussian random fields

This paper presents an algorithm for simulating Gaussian random fields with zero mean and non-stationary covariance functions. The simulated field is obtained as a weighted sum of cosine waves with random frequencies and random phases, with weights that depend on the location-specific spectral density associated with the target non-stationary covariance. The applicability and accuracy of the algorithm are illustrated through synthetic examples, in which scalar and vector random fields with non-stationary Gaussian, exponential, Matérn or compactly-supported covariance models are simulated.     

Simulation of intrinsic random fields of order <Emphasis Type="Italic">k</Emphasis> with a continuous spectral algorithm

Intrinsic random fields of order k, defined as random fields whose high-order increments (generalized increments of order k) are second-order stationary, are used in spatial statistics to model regionalized variables exhibiting spatial trends, a feature that is common in earth and environmental sciences applications. A continuous spectral algorithm is proposed to simulate such random fields in a d-dimensional Euclidean space, with given generalized covariance structure and with Gaussian generalized increments of order k. The only condition needed to run the algorithm is to know the spectral measure associated with the generalized covariance function (case of a scalar random field) or with the matrix of generalized direct and cross-covariances (case of a vector random field). The algorithm is applied to synthetic examples to simulate intrinsic random fields with power generalized direct and cross-covariances, as well as an intrinsic random field with power and spline generalized direct covariances and Matérn generalized cross-covariance.     

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.     

Gamma random field simulation by a covariance matrix transformation method

In studies involving environmental risk assessment, Gaussian random field generators are often used to yield realizations of a Gaussian random field, and then realizations of the non-Gaussian target random field are obtained by an inverse-normal transformation. Such simulation process requires a set of observed data for estimation of the empirical cumulative distribution function (ECDF) and covariance function of the random field under investigation. However, if realizations of a non-Gaussian random field with specific probability density and covariance function are needed, such observed-data-based simulation process will not work when no observed data are available. In this paper we present details of a gamma random field simulation approach which does not require a set of observed data. A key element of the approach lies on the theoretical relationship between the covariance functions of a gamma random field and its corresponding standard normal random field. Through a set of devised simulation scenarios, the proposed technique is shown to be capable of generating realizations of the given gamma random fields.     

Covariance functions motivated by spatial random field models with local interactions

Random fields based on energy functionals with local interactions possess flexible covariance functions, lead to computationally efficient algorithms for spatial data processing, and have important applications in Bayesian field theory. In this paper we address the calculation of covariance functions for a family of isotropic local-interaction random fields in two dimensions. We derive explicit expressions for non-differentiable Spartan covariance functions in \({\mathbb{R}}^2\) that are based on the modified Bessel function of the second kind. We also derive a family of infinitely differentiable, Bessel-Lommel covariance functions that exhibit a hole effect and are valid in \({\mathbb{R}}^{d}\), where d > 2. Finally, we define a generalized spectrum of correlation scales that can be applied to both differentiable and non-differentiable random fields in contrast with the smoothness microscale.     

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

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

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

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

A statistical framework of data fusion for spatial prediction of categorical variables

With rapid advances of geospatial technologies, the amount of spatial data has been increasing exponentially over the past few decades. Usually collected by diverse source providers, the available spatial data tend to be fragmented by a large variety of data heterogeneities, which highlights the need of sound methods capable of efficiently fusing the diverse and incompatible spatial information. Within the context of spatial prediction of categorical variables, this paper describes a statistical framework for integrating and drawing inferences from a collection of spatially correlated variables while accounting for data heterogeneities and complex spatial dependencies. In this framework, we discuss the spatial prediction of categorical variables in the paradigm of latent random fields, and represent each spatial variable via spatial covariance functions, which define two-point similarities or dependencies of spatially correlated variables. The representation of spatial covariance functions derived from different spatial variables is independent of heterogeneous characteristics and can be combined in a straightforward fashion. Therefore it provides a unified and flexible representation of heterogeneous spatial variables in spatial analysis while accounting for complex spatial dependencies. We show that in the spatial prediction of categorical variables, the sought-after class occurrence probability at a target location can be formulated as a multinomial logistic function of spatial covariances of spatial variables between the target and sampled locations. Group least absolute shrinkage and selection operator is adopted for parameter estimation, which prevents the model from over-fitting, and simultaneously selects an optimal subset of important information (variables). Synthetic and real case studies are provided to illustrate the introduced concepts, and showcase the advantages of the proposed statistical framework.     

On the use of Empirical Orthogonal Function (EOF) analysis in the simulation of random fields

In several fields of Geophysics, such as Hydrology, Meteorology or Oceanography, it is often useful to generate random fields, displaying the same variabilitity as the observed variables. Usually, these synthetic data are used as forcing fields into numerical models, to test the sensitivity of their outputs to the variability of the inputs. Examples can be found in subsurface or surface Hydrology and in Meteorology with General Circulation Models (GCM). Different techniques have already been proposed, often based on the spectral representation of the random process, with, usually, assumptions of stationarity. This paper suggests that Empirical Orthogonal Function (EOF) analysis, which leads to the decomposition of the covariance kernel on the set of its eigen-functions, is a possible answer to this problem. The convergence and accuracy of the method are shown to depend mainly on the number of EOFs retained in the expansion of the covariance kemel. This result is confirmed by a comparison with the turning band method and a matrix technique. Furthermore, a synthetic example of non-homogencous fields shows the interest of EOF analysis in the direct simulation of such fields.     

Stochastic analysis of groundwater flow in aquifers with leakage

 Stochastic analysis of one- and two-dimensional flow through a shallow semi-confined aquifer with spatially variable hydraulic conductivity K represented by a stationary (statistically homogeneous) random process is carried out by using the spectral technique. The hydraulic head covariance functions for flows in a semi-confined aquifer bounded by a leaky layer above and an impervious stratum below are derived by assuming that the randomness forcing the head variation to originate from the hydraulic conductivity field of the aquifer. The head covariance functions are studied using two convenient forms of the logarithmic hydraulic conductivity process. The results demonstrate the significant reduction in the head variances and covariances due to the presence of a leaky layer. The hydraulic head correlation distance is also reduced greatly due to the presence of the leaky layer.     

A note on smoothness measures for space–time surfaces

The differentiability of a random field has a direct relationship with the differentiability of its covariance function. We review the concept of differentiability of space–time covariance models and random fields, and its implications on predictions. We analyze the change of behavior of the covariance function at the origin and at different space–time lags away from the origin, by using the concept of smoothness which can be considered the geometrical view of the differentiability. We propose a way to measure the smoothness of any covariance function, and apply it to purely spatial and space–time covariance functions.     

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.     

Conditional simulation of spatially variable seismic ground motions based on evolutionary spectra

This paper proposes a computational procedure for the conditional simulation of spatially variable seismic ground motions for long span bridges with multiple supports. The seismic ground motions, with part of their time histories measured at some supports, are regarded as zero‐mean nonstationary random processes characterized by predefined evolutionary power spectral density. To conditionally simulate unknown seismic ground motion time histories at other supports, the Kriging method is first described briefly for the conditional simulation of a random vector comprised of zero‐mean Gaussian variables. The multivariate oscillatory processes characterized by the evolutionary power spectral density matrix are then introduced, and the Fourier coefficients of the oscillatory processes and their covariance matrix are derived. By applying the Kriging method to the random vector of the Fourier coefficients and using the inverse Fourier transform, unknown nonstationary seismic ground motion time histories can be simulated. A numerical example is selected to demonstrate capabilities of the proposed simulation procedure, and the results show that the procedure can ensure unbiased time‐varying correlation functions, especially the cross correlation between known and unknown time histories. The procedure is finally applied to the Tsing Ma suspension bridge in Hong Kong to generate ground accelerations at its multiple supports using limited seismic records. Copyright © 2012 John Wiley & Sons, Ltd.     

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

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

A stochastic framework for downscaling processes of spatial averages based on the property of spectral multiscaling and its statistical diagnosis on spatio-temporal rainfall fields

Spectral multi-scaling postulates a power-law type of scaling of spectral distribution functions of stationary processes of spatial averages, over nested and geometrically similar sub-regions of the spatial parameter space of a given spatio-temporal random field. Presently a new framework is formulated for down-scaling processes of spatial averages, following naturally from the postulate of spectral multi-scaling, and key ingredients required for its implementation are described. Moreover, results from an extensive diagnostic study are presented, seeking statistical evidence supportive of spectral multi-scaling. Such evidence emerges from two sources of data. One is a 13 year long historical record of radar observations of rainfall in southeastern UK (Chenies radar), with high spatial (2 km) and temporal (5 min) resolution. The other is an ensemble of rain rate fields simulated by a spatio-temporal random pulse model fitted to the historical data. The results are consistent between historical and simulated rainfall data, indicating frequency-dependent scaling relationships interpreted as evidence of spectral multi-scaling across a range of spatial scales.     

Permissibility of fractal exponents and models of band-limited two-point functions for fGn and fBm random fields

The fractional Gaussian noise (fGn) and fractional Brownian motion (fBm) random field models have many applications in the environmental sciences. An issue of practical interest is the permissible range and the relations between different fractal exponents used to characterize these processes. Here we derive the bounds of the covariance exponent for fGn and the Hurst exponent for fBm based on the permissibility theorem by Bochner. We exploit the theoretical constraints on the spectral density to construct explicit two-point (covariance and structure) functions that are band-limited fractals with smooth cutoffs. Such functions are useful for modeling a gradual cutoff of power-law correlations. We also point out certain peculiarities of the relations between fractal exponents imposed by the mathematical bounds. Reliable estimation of the correlation and Hurst exponents typically requires measurements over a large range of scales (more than 3 orders of magnitude). For isotropic fractals and partially isotropic self-affine processes the dimensionality curse is partially lifted by estimating the exponent from measurements along fixed directions. We derive relations between the fractal exponents and the one-dimensional spectral density exponents, and we illustrate the relations using measurements of paper roughness.The author would like to acknowledge helpful comments from two anonymous referees.     

An improved spectral turning-bands algorithm for simulating stationary vector Gaussian random fields

We propose a spectral turning-bands approach for the simulation of second-order stationary vector Gaussian random fields. The approach improves existing spectral methods through coupling with importance sampling techniques. A notable insight is that one can simulate any vector random field whose direct and cross-covariance functions are continuous and absolutely integrable, provided that one knows the analytical expression of their spectral densities, without the need for these spectral densities to have a bounded support. The simulation algorithm is computationally faster than circulant-embedding techniques, lends itself to parallel computing and has a low memory storage requirement. Numerical examples with varied spatial correlation structures are presented to demonstrate the accuracy and versatility of the proposal.     

Testing the type of non-separability and some classes of space-time covariance function models

In statistical space-time modeling, the use of non-separable covariance functions is often more realistic than separable models. In the literature, various tests for separability may justify this choice. However, in case of rejection of the separability hypothesis, none of these tests include testing for the type of non-separability of space-time covariance functions. This is an important and further significant step for choosing a class of models. In this paper a method for testing positive and negative non-separability is given; moreover, an approach for testing some well known classes of space-time covariance function models has been proposed. The performance of the tests has been shown using real and simulated data.     

New anisotropic covariance models and estimation of anisotropic parameters based on the covariance tensor identity

 Many heterogeneous media and environmental processes are statistically anisotropic. In this paper we focus on range anisotropy, that is, stochastic processes with variograms that have direction dependent correlation lengths and direction independent sill. We distinguish between two classes of anisotropic covariance models: Class (A) models are reducible to isotropic after rotation and rescaling operations. Class (B) models can be separated into a product of one-dimensional functions oriented along the principal axes. We propose a new Class (A) model with multiscale properties that has applications in subsurface hydrology. We also present a family of Class (B) models based on non-Euclidean distance metrics that are generated by superellipsoidal functions. Next, we propose a new method for determining the orientation of the principal axes and the degree of anisotropy, i.e., the ratio(s) of the correlation lengths. This information reduces the degrees of freedom of anisotropic variograms and thus simplifies the estimation procedure. In particular, Class (A) models are reduced to isotropic and Class (B) models to one-dimensional functions. Our method is based on an explicit relation between the second-rank slope tensor (SRST), which can be estimated from the data, and the covariance tensor. The procedure is conceptually simple and numerically efficient. It is more accurate for regular (on-grid) data distributions, but it can also be used for sparse (off-grid) spatial distributions. In the case of non-differentiable random fields the method can be extended using generalized derivatives. We illustrate its implementation with numerical simulations.