Physics Mathematics
ComputerScience,1970-2008-Computer Science-2006.jsonl.gz.gz
Computer Science
51ba20106d5b88404df05d7bdacfd19a1c183206
Fast Spatial Prediction from Inhomogeneously Sampled Data Based on Generalized Random Fields with Gibbs Energy Functionals

An explicit optimal linear spatial predictor is derived. The spatial correlations are imposed by means of Gibbs energy functionals with explicit coupling coefficients instead of covariance matrices. The model inference process is based on physically identifiable constraints corresponding to distinct terms of the energy functional. The proposed predictor is compared with the geostatistical linear optimal filter (kriging) using simulated data. The agreement between the two methods is excellent. The proposed framework allows a unified approach to the problems of parameter inference, spatial prediction and simulation of spatial random fields.


INTRODUCTION
Spatial prediction of physical variables from samples that are irregularly distributed in space is a task with applications in many fields of science and engineering, including subsurface hydrology [1,2], oil reservoir engineering [3,4], environmental pollutant mapping and risk assessment [5], mining exploration and reserves estimation [6], environmental health studies [7], image analysis [8] and neuroscience [9]. Physical quantities of economical and environmental interest include mineral grades, concentrations of environmental pollutants, soil and rock permeability and flow fields in oil reservoirs. Modeling the variability of such processes is based on the theory of spatial random fields (SRFs) [10]. Knowledge of spatial correlations in SRFs enables (i) generating predictive isolevel maps (ii) estimating prediction uncertainty and (iii) developing simulations that reconstruct probable scenarios conditioned on the data. The classical approach is based on Gaussian SRF's (GSRF's) and various generalizations for non-Gaussian distributions [11,12]. For GSRF's the spatial structure is determined from the covariance matrix, which is estimated from the available sample(s).
Let Ω ∈ R d denote the area of interest, and |Ω| its volume. An SRF state (realization) in Ω can be decomposed into a deterministic trend m x (s), a correlated fluctuation X λ (s), and an independent random noise term, ǫ(s), i.e., X(s) = m x (s) + X λ (s) + ǫ(s). The trend represents large-scale variations obtained in principle by ensemble averaging, i.e. m x (s) = E[X(s)]. In practice, the trend is often determined from a single available realization. The fluctuation represents 'fast variations' describing fine structure above the resolution limit λ. The random noise represents non-resolved inherent variability, purely random additive noise, or non-systematic measurement errors. The fluctuation typically is a second-order stationary SRF, or an intrinsic SRF with second-order station-ary increments [10]. The residual SRF after trend removal is a zero-mean fluctuation: X * (s) = X λ (s) + ǫ(s). The inference process focuses on determining a model for X λ (s) from a set of possibly noisy observations X * (s).
In statistical physics the probability density function (pdf ) of any fluctuation field X(s) governed by an en- In classical geostatistics, the Gaussian joint pdf for a set of fluctuations ij is the inverse covariance matrix, and summation is implied over repeated indices. Instead, Spartan Spatial Random Fields (SSRF's) [13] model spatial correlations in terms of 'interactions'. This change in viewpoint has important consequences for model inference and spatial prediction.

THE FGC MODEL
In [13] the fluctuation-gradient-curvature (FGC) SSRF model is defined and its properties investigated. The continuum FGC model involves the following H[X]: The model is characterized by four parameters: the scale coefficient η 0 , the covariance-shape coefficient η 1 , the characteristic length ξ, and the cutoff wavevector k max . Bochner's permissibility theorem [5] for the positive definiteness of the covariance function requires η 1 > −2 if k max → ∞. In statistical physics terminology, ℓ 1 = η 1 ξ 2−d /(2 η 0 ) and ℓ 2 = ξ 2−d /(2 η 0 ) represent the coupling strengths of the gradient and curvature terms. A coarse-graining kernel is used to cut off the fluctuations at k max [13,14], leading to bandlimited covariance spectral density. If k max is finite, the field's configurations are almost surely differentiable [14]. If k max is infinite, generalized gradient and constraints should be used. The coarse-graining kernel implies that the SRF X λ is a generalized SRF [10]. A moment-based method for parameter estimation was proposed and validated with simulated data [13]. The inference process is based on matching ensemble constraints E[S j (s)] with their sample counterparts, denoted by S j (s), for j = 0, 1, 2. The procedure is extended in [15].
Assume S m = {s 1 , . . . s N } is a set of sampling points on an irregular grid and X * (S m ) = {X * 1 , . . . , X * N } is the respective vector of measurement. On an irregular grid, the translation symmetry of the lattice is lost. The continuum FGC functional is then a more suitable model. For practical purposed, a tractable approximation of the continuum model is needed.
In [15], approximations for the sample averages of S 1 (s) and S 2 (s) are formulated in terms of kernel averages of the data values. In the following we use the notation: c , where the summation is over both indices (i, j = 1, . . . , N ) denotes the average of the quantity A i,j , weighted by the kernel K h (r) with bandwidth parameter h. The pair distance is denoted by s i,j = s i − s j , where r is the Euclidean norm of the distance vector r, and the field increment by Then, the generalized gradient constraint is given by where a 2 1 = s 2 i,j h1 . Sensible estimates of the spacing a 1 should account for the grid topology. E.g., let B 0 be the set of near-neighbor vectors of all the points in S m . If B 0 contains N 0 vectors and ∆ i denote the lengths of the vectors in B 0 , thenâ d where µ 1 and µ 2 are 1 + o(ǫ) constants that depend on the sampling network topology. The µ 1 , µ 2 are defined so as to satisfy asymptotic bias and consistency properties [15]. They introduce explicitly in the problem of model inference four parameters linked to the topology of the sampling network: the spacings a 1 and a 2 that replace the lattice constant, and the bandwidths h 1 and h 2 that determine the range of influence of the averaging kernel. The latter are determined from the consistency principle a 2p The predictions will be denoted by {X λ (z l ), l = 1, . . . , K} and the respective prediction vector byX λ (Z p ). The increments corresponding to V l will be denoted by α p (V l ), p = 1, 2. Typically, single-point prediction is applied sequentially over all points in Z p . Multiple-point prediction is possible in the SSRF framework, but this letter focuses on single-point prediction.

Optimal Linear Prediction
In geostatistics, spatial prediction is based on the Best Linear Unbiased Estimator (BLUE), commonly known as Kriging [1,12]. Different variants of kriging exist, depending on the hypotheses about the normality of the data and the behavior of the mean. These methods are generalizations of the linear minimum mean square error (LMMSE) estimators, also known as Wiener filters [16]. Ordinary kriging (OK) is the most common variety. It is applied to normally distributed data, with an unknown mean that can be considered as locally constant. A single-point prediction is obtained as a superposition X(z l ) = M j=1 λ j X(s j ), where s j are all the points inside a local search neighborhood, B(s l ), around z l . The prediction error is defined as ε(s l ) =X(z l )−X(z l ). The optimal linear coefficients should minimize the mean square error conditional on the zero-bias constraint M j=1 λ j = 1, i.e., the expression E ε 2 (s l ) +µ M j=1 λ j − 1 , where µ is a Lagrange coefficient. This leads to the linear system C X (s i − s j ) λ j + µ = C X (s i − z l ), ∀ i = 1, . . . , M and M j=1 λ j = 1, where C X (r) is the centered covariance function. OK is an exact interpolator, meaning that {X(s i ) = X * (s i ), ∀s i ∈ S m }. Exactitude is not always desirable, since it ignores measurement errors and over-constrains the predictions. For a fixed-size search neighborhood, the numerical complexity of an M -point OK prediction is O(K M 3 ).

Spatial Prediction based on FGC Functional
Once the parameters of the FGC model have been determined from the data, prediction of the SRF at z l is possible by means of at least two different approaches.
First, the corresponding covariance function is determined and spatial predictions are then obtained using the OK predictor. In this case, the only difference introduced by the SSRF functional is the covariance estimator.
Here we propose a different predictor, obtained by maximizing the conditional probability density , the prediction is obtained by maximizing f X [X λ (V l )]. In principle, this requires solving the equation δH fgc /δX λ (z l ) = 0, where δ [.] δX λ (z l ) is the variational derivative of the functional given by Eq. (1) with respect to X λ (z l ). In practice, H fgc is replaced by a discretized estimator,Ĥ fgc . SinceĤ fgc is a bilinear functional, the prediction follows from the solution of the linear equation: H fgc is obtained by means of the estimators {S j (s), j = 0, 1, 2}, using the ergodic hypothesis ds S j (s) ≈ |Ω| S j (s), which leads to: where the spatial averages involve the sampling points and the prediction point as well. In light of (2), (3), and (5), equation (4) leads to the following linear predictor where q 1 = q 2 = 1, q 3 = q 4 = −1, X * l h is the kernel average of the sample, centered at the prediction point and the linear coefficients {c i , i = 1, 2, 3, 4} are given by The summation in g h (V l ) extends over all the N (N −1)/2 non-identical and non-repeating pairs of sampling points. Defining the linear weights the prediction is expressed aŝ

Properties of the FGC Mode Predictor
The present formulation of the FGC mode predictor (FGC-MP) is closer to simple kriging than OK, since the mean is assumed to be known. However, unlike simple kriging the mean does not have to be constant, provided that it changes slowly so that the energy contributions due to the square gradient and curvature of the mean can be ignored compared to the fluctuations. In this respect, the predictor resembles OK, which allows for slow (but unknown) variation of the mean. Predictions of the FGC-MP are independent of η 0 , while the prediction variance is linearly proportional to η 0 .
The FGC-MP is linear and unbiased. Since the joint FGC pdf is Gaussian, the mode estimate is equivalent to the minimum mean square estimate. Thus, the FGC-MP is an optimal linear predictor. The main differences with OK result from the use of the energy functional in the FGC SSRF: (1) The FGC-MP is not an exact interpolator, because it does not use the data at the prediction point. (2) The single-point FGC-MP provides an explicit expression for the prediction, while kriging requires solving a linear system. (3) The FGC-MP does not require specifying a search neighborhood around the prediction point; in kriging definition of a search neighborhood requires an iterative procedure based on cross-validation of the predictions with the data. (4) The FGC-MP incorporates two sets of parameters: the first set determines the spatial dependence of the SRF, while the second set depends on the topology of the sampling network. The influence of the sampling topology is not explicitly accounted for in kriging. (5) The uncertainty estimate involves the SSRF covariance function, for which there are no explicit solutions in d = 2, unlike d = 1, 3 [14]. Obtaining the covariance in d = 2 requires performing numerically a univariate (for isotropic dependence) integration of the spectral density. (6) Regarding multiplepoint estimates the FGC-MP has a numerical complexity O(K 3 ), derived from solving a linear system of K coupled equations at the prediction points, while the numerical complexity of kriging is O(K M 3 ).

PREDICTION USING SIMULATED SAMPLES
At 400 randomly distributed points on a square domain of length L = 100 we generate 100 independent "samples". These represent realizations of a Gaussian random field with m x = 50, and an exponential covariance function C x (r) = σ 2 x exp(− r /b e ), where σ x = 10, and b e = 10. The Cholesky LU decomposition method is used for the simulations. We partition the 400 points into a training set S m of 100 randomly selected points, and a prediction set, Z p , including the remaining points. We use the first set to determine the optimal SSRF parameters, and then predict the values of the field at the locations of the prediction set. The triangular kernel is used in the FGC-predictor mode. Predictions are also generated using the Ordinary Kriging method. The performance of the predictors is evaluated using the bias, the mean absolute error (mae), the root mean square error (rmse), the mean absolute relative error (mare), the root mean square relative error (rmsre) and the linear correlation coefficient (R 2 ). The means are calculated with respect to the values at the 300 prediction points. Statistics of these quantities over the 100 samples are shown in Tables I and II. The kriging predictor is applied with the a priori parameters of the exponential covariance (instead of the inferred covariance model from the data). This choice aims at testing the FGC-Mode Predictor against the "true" model. The results show that the two predictors perform very similarly.

CONCLUSIONS
A fast linear optimal predictor, with applications in the analysis of spatial data, is proposed. The predictor is based on generalized random fields which represent the spatial dependence in terms of interactions. An explicit expression for single-point prediction is obtained. The reduced numerical complexity of the FGC-Mode predictor may promote the use of cross-validation procedures for model parameter inference, instead of the commonly used parametric methods. The SRF representation, which is based by construction on an objective function, provides a unified framework for model parameter estimation, spatial prediction and constrained (respecting the data) simulation. This is in contrast with classical approaches that require the introduction of ad hoc objective functions [4,6] for simulations (e.g., by means of simu-lated annealing.) Expressions for the prediction uncertainty and a linear system for multiple-point prediction have also been derived and will be reported elsewhere [17]. The multiple-point predictor accounts for interactions between the prediction points that may lower the total "energy". Such interactions are missed in singlepoint prediction. Finally, the FGC focuses on short-range interactions, but long-range dependence can be incorporated in the SSRF framework with suitable modifications of the energy functional.