The Sequential Generation of Gaussian Random Fields for Applications in the Geospatial Sciences †

This paper presents practical methods for the sequential generation or simulation of a Gaussian two-dimensional random field. The specific realizations typically correspond to geospatial errors or perturbations over a horizontal plane or grid. The errors are either scalar, such as vertical errors, or multivariate, such as x, y, and z errors. These realizations enable simulation-based performance assessment and tuning of various geospatial applications. Both homogeneous and non-homogeneous random fields are addressed. The sequential generation is very fast and compared to methods based on Cholesky decomposition of an a priori covariance matrix and Sequential Gaussian Simulation. The multi-grid point covariance matrix is also developed for all the above random fields, essential for the optimal performance of many geospatial applications ingesting data with these types of errors.


Introduction and Motivation
This paper identifies a specific and practical subclass of homogeneous Gaussian two-dimensional (2D) random fields and presents a simple, fast, sequential method to generate discrete realizations over a (horizontal) grid for the purpose of Monte Carlo simulation-based analyses.Let us term this method Fast Sequential Simulation (FSS) for brevity of further description.FSS can be considered an extension of the sequential generation of a first order Gauss-Markov process from a 1D function of time to a 2D function of horizontal space.Although FSS was derived independently, it is also demonstrated a special case of Sequential Gaussian Simulation which is commonly used in the Geostatistics community.In particular, FSS is an unconditional simulation with simplicity and speed due to both exponential correlation in the spatial directions and an ordered generation over an evenly spaced grid of horizontal locations.Although other applications of Sequential Gaussian Simulation are more general (conditional or unconditional, irregularly spaced points, random generation paths, arbitrary valid correlation functions, etc.), many applications do not require these generalities but do require speed, preferably with a simple and direct implementation.
The paper first addresses scalar random fields, i.e., where typically represents a scalar error or perturbation at grid location .The desired variance and spatial correlations for the are specifiable with FSS, and the multi-grid point covariance matrix derived.The paper then generalizes FSS to the generation of multivariate Gaussian two-dimensional random fields, i.e., , where is a vector of arbitrary dimension .Finally, the paper generalizes FSS results even further to non-homogeneous Gaussian two-dimensional random fields, where the variance and spatial correlations are a function of grid location .Some of the practical techniques presented for the sequential generation of both multivariate and non-homogeneous random fields are believed to be new and somewhat innovative.
Example realizations of scalar, multivariate, homogenous, and non-homogenous random fields are presented throughout the paper, as well as various theoretical properties, insights, and proofs, the latter contained in appendices.FSS is also compared to equivalent methods based on (1) Cholesky decomposition of a pre-computed a priori covariance matrix; and (2) Sequential Gaussian Simulation as implemented in various statistical packages.FSS is demonstrated to be many orders of magnitude faster than all of these other generation methods, as well as being a simpler implementation.
The ability to simulate errors across a horizontal grid with specifiable expected magnitudes (variance) and interrelationships (correlations) is an important capability in support of the Geospatial Sciences and supported by the FSS method presented in this paper.For example, the errors can represent elevation errors across a Digital Elevation database, horizontal errors in the location of vertices across a GIS database, horizontal and vertical errors in the locations of control points across a control point database, etc.All of these errors are essentially a function of horizontal location, i.e., representable as a two-dimensional random field.
These simulated errors can be used to modify corresponding "truth" data in a simulation environment.Subsequent performance of various down-stream applications can then be meaningfully assessed, including modification (tuning) of their algorithms for optimal and reliable performance.Alternatively, in an operational environment, the applications themselves can have an embedded simulation capability in order to represent the effects of errors in their input data of known (specifiable) a priori characteristics.The effect is relative to the application's output product and usually represented graphically.The simulation of tens of millions of errors within a few seconds and hundreds of millions within 30 s on a laptop computer is desired.
Previously, relevant errors have sometimes been simulated as homogeneous errors solely as an assumption for reduced complexity and/or increased speed.However, many realistic applications correspond to data with non-homogeneous error characteristics; for example, data sets previously fused from other data sets with differing error (accuracy) characteristics.This paper addresses both types of errors.The non-homogeneous techniques presented in this paper essentially preserve the speed associated with the technique presented for the homogeneous case, typically reducing the speed by only a factor of two or three.The corresponding non-homogeneous characteristics are not totally general, but still adequate for many applications.
Finally, a common theme throughout this paper is the practical computation and need for a multi-grid point covariance matrix corresponding to or at multiple grid point locations.It is used by various applications to predict the accuracy of their input data and properly weight it within their various algorithms.
The authors of [1][2][3][4][5] discuss random fields in general, including their generation or simulation.Generation techniques include those based on Cholesky decomposition and Sequential Gaussian Simulation.In addition, these references discuss interpolation of a random field's realization based on Kriging.These references are relatively standard in the geostatistics community.They, along with other references from this community, are referenced per specific topic throughout the remainder of this paper, including appendices.
As detailed later, FSS was derived independently of Gaussian Sequential Simulation but is equivalent in specified circumstances.FSS is also directly related to both generalized multi-grid point covariance matrices [6] and strictly positive definite correlation functions [7] that have applications in the Geospatial community.Recent applications of FSS include evaluating conflation methods [8] and various geospatial algorithms [9].

Roadmap
Sections 1-4 of this paper define the scalar homogeneous Gaussian 2D random field, the fast sequential generation algorithm FSS, and related practical aspects.Section 5 compares FSS to more typical generation methods such as those based on Cholesky decomposition or Sequential Gaussian Simulation, particularly with respect to timing or throughput.Section 6 then extends the FSS technique to a multivariate homogeneous Gaussian 2D random field, and Section 7 to a non-homogeneous Gaussian 2D random field.

A Scalar Gaussian 2D Random Field and Its Sequential Generation
In this section of the paper, we define a scalar, homogeneous, Gaussian two-dimensional (2D) random field, typically corresponding to perturbations or errors.We also present FSS, the algorithm for the fast and sequential generation of a discrete and specific realization of this random field that can be used for various Monte Carlo simulation-based analyses.
We assume that the random field corresponds to -error for specificity.Indices and correspond to a two-dimensional grid in the -horizontal plane: is in the " " direction and is the first grid index, is in the " " direction and is the second grid index.Specifically, corresponds to -error in meters at grid point .

Statistical Characteristics
is normally (Gaussian) distributed, and has a mean value of zero and a specifiable one-sigma , i.e., is normally distributed for all grid locations .Its spatial correlation across the grid is separable, i.e., has the (normalized) correlation function , where and are the absolute values of the component-wise differences in the location of two arbitrary grid points.This function is represented as: (1) where and are specifiable spatial correlation distance constants (meters) and and specifiable grid spacing (meters) in the and the directions in the horizontal plane, respectively.Note that and .
Figure 1 presents an example of , with m, m, and m.
The spatial correlation is applicable to any pair of grid points within the entire grid that are separated by meters in the -direction and meters in the -direction.The use of two different spatial correlation distance constants allows for specification of different correlation characteristics in each of the horizontal directions.Regarding the a priori statistics of in a more formal manner: (meters) at an arbitrary location within the grid.
and when .Section 3.0 of this paper also presents the covariance matrix associated with two or more , each associated with a different grid point .

Core Grid-Generation Equation
Equation ( 3) is the core grid-generation equation for FSS: ( The integers and correspond to points in the grid, , and . is a random sample of Gaussian white noise, and is normally distributed , where .That is, given a desired , , and , a corresponding value of is computed per the above. is the spatial correlation of the scalar error between adjacent grid points in the (or ) direction ( , unit-less), and is the spatial correlation of the scalar error between adjacent grid points in the (or ) direction ( , unit-less) Therefore, we can also write the spatial correlation function as a function of grid units only: Figure 2 illustrates the grid of errors generated based on Equation (3) and corresponding to a specific realization of the scalar, homogeneous, two-dimensional random field.All errors in the light orange area affect the error and those in the light blue area do not.Also, based on the specific sequential grid generation algorithm presented in Section 2.3, an error in the light blue area (e.g., ) may be generated prior to even though it does not affect it.

Sequential Generation Algorithm for Realization over a pxq Grid
The following presents a specific FSS algorithm for a discrete realization of over a grid based on Equation ( 3 " corresponds to a random number (realization) from a probability distribution; for example, in matlab this is implemented as " ".Appendix A of this paper presents a direct proof that the above FSS algorithm generates a realization of a two-dimensional random field with the specifiable statistical properties presented in Section 2.1.Appendix B also demonstrates its mathematical equivalence to a corresponding Sequential Gaussian Simulation approach for completeness.The latter must specify separable exponential correlation in the spatial directions and a fixed grid with a specific ordered path across it for generation of the realization.Also, depending on how it is implemented, it may or may not take advantage of the need to use the realization at only three previously generated grid points for generation at the current grid point, as does FSS.If it did in an efficient manner using pre-computed Kriging weights and little overhead due to flexibility and complexity, its speed could approach that of FSS.

Grid Spacing
Equation (3) and the above algorithm should typically incorporate grid spacing and equivalent to approximately one-ninth or less their respective spatial correlation distance constant, insuring at least 0.9 correlation with an adjacent grid point, i.e., and , or equivalently, m and m.Of course, this "rule-of-thumb" is application dependent.For example, if very high spatial correlation between adjacent grid points is of interest, spacing should be closer.

Grid Buffer
As shown in Appendix A of this paper, the statistical properties of the are based on steady-state properties over an assumed infinite horizontal grid.Thus, for an actual application (realization) necessarily using a finite grid, the "final" grid should have a "buffer" on two edges of the computed grid to ensure that steady-state has essentially been reached.This is illustrated in Figure 3, where the buffer is yellow and the final grid is green.Placement of the buffer corresponds to the specific sequential grid generation algorithm presented earlier that starts at the top of the grid and always proceeds from right to left.The width of the top buffer should correspond to the equivalent of approximately two times the spatial correlation distance constant in the -direction, and the width of the side buffer should correspond to the equivalent of approximately two times the spatial correlation distance constant in the -direction.
More specifically, width of top buffer (# grid rows) = and width of side buffer (# grid columns) = 2 .Or equivalently, if and equal the value 0.9, 19 grid rows and 19 grid columns.This will ensure generation of errors throughout the final grid with the desired statistical properties.

Example Realizations: Surface Plots
This section presents surface plots of the -error over a subset of a 2D final grid generated using the sequential algorithm of Section 2.3.Example 1 corresponds to specified m, and specified (thus m).Assuming a grid spacing in both the y and x directions of 1 m ( ), this corresponds to spatial distance constants equal to m.In this particular case, the spatial distance constants and were derived from the specified and , given assumed grid spacing and , not vice versa.The spatial distance constants were computed for information only.That is, there are two basic but equivalent approaches for the specification, application, and interpretation of spatial correlation, the particular approach selected based on convenience: Approach 1-specify spatial correlation by the values of and (unitless) directly, implement Equation (3), and then interpret location-dependent results in terms of grid units.Given assumed grid spacing and (meters), the spatial distance constants and (meters) can be derived for information purposes only.Approach 2-specify spatial correlation by the values of and (meters) and grid spacing and (meters), compute and , implement Equation (3), and then interpret location-dependent results in terms of y-x horizontal space in meters.The approach works well when the generated random field is to correspond to the a priori statistics and spatial resolution of a specific application of interest in the Geospatial Sciences.m).As expected, the above realizations over portions of the grid do not have a mean -error of zero nor a standard deviation about the mean of 10 m.However, when sample statistics are computed over numerous realizations, the corresponding mean and standard deviation approach 0 m and 10 m, respectively, matching the common a priori statistics used to generate the realizations.

Example Realizations: Sample Statistics
Sample correlation functions or correlograms were computed for two independent realizations across a final z-grid 1000 × 1000 in size.A priori statistics were specified with a fixed mean value of 0 and a standard deviation about the mean of m.The first realization corresponded to a priori correlations represented by and , and is presented in were plotted as a function of horizontal distance in the x-direction and horizontal distance in the y-direction, and are different as expected per the values of r and s.The second realization corresponded to , and is presented in Figure 9. Three different horizontal directions were evaluated: x, y, and 45 degrees between.Note that results for the latter are different even though because the FSS correlation model is not isotropic.(Note that 45 degrees yields the direction with maximum difference.)In general, both plots demonstrate nearly identical results between the true and sample correlation functions-not unexpected because FSS has virtually no approximations and because the random field is ergodic and the number of samples within a given realization large.200x200 grid; a priori mean = 0 m, sigma = 10 m, r = 0.9, s = 0.9; correlation spatial direction: y-direction (blue); multiple realizations typically of interest in the geostatistics community, computed in the y-direction across the horizontal grid.(See [5] for example, definitions of the correlogram and semivariogram.)Note the reasonable variability of the sample semivariograms corresponding to each of the five realizations about the common theoretical semivariogram.

Multi-Grid Point Covariance Matrix
If there are specific scalar errors of interest associated with arbitrary and different grid locations , their corresponding (joint) covariance matrix is symmetric and positive definite (valid) since all the grid point errors have the same variance and have inter-grid point correlation between pairs corresponding to a normalized strictly positive definite correlation function (spdcf) , i.e., the multi-grid point covariance matrix equals: (5) where is an vector such that and the , , correspond to an ordered list of the grid point locations.Also, and are the y and x distances in meters in the horizontal plane between the ordered points ; directly multiplies each element of the matrix.(Alternatively, the spatial correlation function and distances could have been written based on grid units.)Note that the above is an a priori covariance matrix corresponding to the various considered as random variables, not specific realizations.See Reference [7] regarding the properties of a spdcf such that the above is guaranteed a valid (symmetric and positive definite) covariance matrix regardless the size of .In general, just because the absolute value of correlation between two arbitrary grid point locations is less than 1.0, this by itself does not ensure a valid multi-grid point covariance matrix for The FSS generation of a realization of over a 2D grid as presented in Section 2.3 did not require use of the explicit multi-grid point covariance matrix in the generation process.So, why is the calculation of this covariance matrix of interest in terms of specific applications of generated errors or perturbations?A major reason is as follows: An "analysis module" may generate the simulated grid of errors and apply a subset to "truth data" and then pass the composite data to a "down-stream" application such that its performance can be assessed in the presence of errors.Many such applications also require knowledge of the multi-grid point covariance matrix corresponding to the composite data for purposes of proper weighting of the composite data in various estimation procedures (Kalman Filter and Weighted Least Squares estimators) for the parameters (state vectors) of interest to the application [6].Of course, these applications can simply be passed, along with the composite data, the corresponding and the parameters that define , such that the applications can then build the appropriate multi-grid point covariance matrix themselves.

Homogeneity and Gaussian Joint Probability Density
The scalar errors generated using the FSS sequential generation technique are Gaussian distributed as they are a linear combination of the , which are Gaussian distributed by definition.(The linear combination is demonstrated explicitly in Appendix A.) Also, the corresponding grid of -errors corresponds to a wide-sense homogenous random field since the variance and correlation of errors are invariant of specific absolute grid location(s).In addition, since the errors are Gaussian distributed, a wide-sense homogeneous random field is equivalent to a homogeneous random field [2].
Any finite collection of at different grid locations contained in the vector has a corresponding Gaussian joint probability density function defined as follows: (6) where is the multi-grid point covariance matrix, and is the matrix determinant.Thus, probabilities can be assigned in a straightforward manner to any absolute or relative confidence interval of interest.
Finally, it is worth noting that all of the multi-grid point covariance matrices computed per this paper are valid, regardless the specific underlying probability distribution of errors.This is true for both the scalar error of Sections 1-5, multivariate errors of Section 6, and non-homogenous errors of Section 7.This is discussed in References [6,7], which also allow for the use of any valid family of spdcf.In addition, the authors of [6] discuss the importance of such covariance matrices, other generation methods for such covariance matrices, and how to generate corresponding probability error ellipsoids.Note that in Reference [6], these covariance matrices are more generally termed "multi-state vector error covariance matrices".

Interpolation into the Grid
The FSS technique as described in Section 1 and 2 of this paper generates a realization of a random field at grid point locations only.This is perfectly adequate for many applications since the grid can be very large and dense.However, if scalar errors are desired between grid point locations, interpolation of the at the four enclosing grid point locations may be performed.Also, the multi-grid point covariance (Equation ( 5)) can be easily modified for a corresponding set of interpolated points.Simply modify the distances between grid points to corresponding distance between the locations of the interpolated points.These distances may be represented in either non-integer units of grid spacing or corresponding and distances in meters in the horizontal plane.

Related Effects
Dependent on the interpolation method, the actual a priori statistics (one-sigma value, spatial correlations) for the interpolated may be somewhat different than the assumed values as represented by the modified multi-grid point error covariance discussed above.The latter is consistent with the a priori statistics of the random field assuming no interpolation.(For a priori statistics, the interpolated value of z is treated as a random variable, not an estimate of a realization.)For nearest neighbor interpolation, there are no explicit differences between the actual and assumed statistics as the location of the point is assumed a (nearest) grid point location.Of course, there can be implicit differences: for example, if two points for interpolation are close but not identical in location, they can be assigned to the same grid point with corresponding 100% spatial correlation between their errors whereas the actual spatial correlation is less.
For bi-linear interpolation, there are differences due to the "averaging" of uncorrelated components of error in the surrounding used during bi-linear interpolation.The higher the correlation between adjacent grid points (the larger and ), the less the effect (differences).If the recommended a priori correlation between adjacent grid points (0.9 or greater per Section 2.3.1) is used during grid generation, the effect is minimal.For example, if , the actual a priori one-sigma value for the interpolated points is 0%-5% less than the a priori one-sigma value for the grid points, and the actual a priori correlation between interpolated points 0%-3% less than for grid points at corresponding distances.The actual value for the percentage difference is dependent on how close the interpolated point(s) is to a grid point.Figure 11 illustrates bilinear interpolation.

Comparison to Alternate Generation Methods
This paper presents FSS, a fast and efficient sequential method for the generation of a 2D grid of errors or perturbations.There is also an implied, associated multi-grid point covariance matrix (e.g., Equation ( 5)), but this covariance matrix is not needed in the generation process.On the other hand, the spatial correlation of errors with this generation method is limited to a specific spdcf family of spatial correlation functions, i.e., . (Albeit, reasonably general in that the distance constants and are specifiable).
There are two other general approaches to the generation (simulation) over a 2D grid: (1) matrix square roots; and (2) Sequential Gaussian Simulation.The latter is sequential and, as mentioned previously, more general than FSS.The former is also more general, but not sequential.They are described in more detail in the next subsection.

Timing Comparisons among Simulation Techniques
Figure 12 shows CPU time comparisons among five different techniques for simulating perturbations for a square grid by varying the number of points n; where is the number of points along one side of the grid.Computation times were measured with a PC laptop with Intel i5 dual core 2.3 GHz CPUs and 8 GB of memory.
The objective was to measure the CPU performance of the main computation (sans overhead setup) for each method.Efforts were made to match the modeling parameters among all methods as closely as possible.Testing assumed unconditional, homogeneous, and isotropic models only (actually, for FSS the model was approximately isotropic, as ).The following describes what main computations were timed in each of the five methods, and are listed in ascending order of computational speed gain according to Figure 12.
(1) Principal matrix square root (using Matlab function SQRTM) (7) where is the nxn principal matrix square root of ; is a n × 1 vector realization of n independent N(0,1) distributed random variables, and is the n × 1 vector of perturbations corresponding to the random variable z or over the grid.was assumed to be a full, and positive definite matrix, i.e., the a priori covariance matrix corresponding to the random field at the n different points in the grid.Matlab uses the Schur decomposition technique to compute SQRTM for a general square matrix, which can be further sped up for symmetric and real matrices.over a 2D square grid (Note that FSS is cutoff at ~15 s due to reaching the system memory limit).
(2) Cholesky decomposition (using Matlab function CHOL) (8) where L is the lower triangular nxn matrix from the Cholesky decomposition, , where is the conjugate transpose of L; is a n × 1 vector realization of n independent N(0,1) distributed random variables, and is the n × 1 vector of perturbations corresponding to the random variable z or over the grid.
was assumed to be a full and positive definite matrix.
(3) PREDICT.GSTAT (version 1.0, 19 April 2014) [10] is the algorithm based upon Pebesma [11] as implemented and tested in the "R" (version 3.0.2,64 bit) statistical package.The following R script is an example of parameters used to time unconditional simulation on a 100 × 100 grid.Note that only the execution of the PREDICT function was timed.
(5) Fast Sequential Simulation (FSS) is the technique described in this paper, and was coded and timed as a Matlab function.The main required parameters used were grid spacing ( ), standard deviation for the random variable ( ), and the spatial distance correlation constants ( ), as all described in Sections 2.3-2.4 of this paper.

Discussion of Timing Results
The principal matrix square root (SQRTM) and Cholesky decomposition (CHOL) methods were provided to serve as a starting benchmark.While they are the least practical for large n, their main benefit is providing an exact solution for any spatial distribution of points and any a priori spatial statistics (valid covariance matrix).Figure 12 shows Cholesky providing roughly half an order of magnitude speed gain over SQRTM.
PREDICT.GSTAT and VISIM provide implementations of standard geostatistical techniques for Sequential Gaussian Simulation. Figure 12 shows that they provide comparable speed performance, and are 1-2 orders of magnitude faster than SQRTM and Cholesky.Their main advantage is providing broad flexibility for general purpose modeling among conditional and unconditional simulations.Moreover, additional speed efficiencies can be achieved when simulating multiple realizations with a fixed parameter set, which is not captured in Figure 12.E.g., following a single random path through the locations, PREDICT.GSTAT reuses results for each of the subsequent simulations [10].
FSS is the technique proposed in this paper.The two main advantages of FSS are (1) speed gain, e.g., three orders of magnitude faster than the next fastest technique as shown in Figure 12; and (2) simplicity of operation, e.g., requiring only three main parameters.Note that the FSS curve in Figure 12 is cutoff at ~15 s, which was due to reaching the memory limit for the grid size (n > 2 × 10 8 points).However, this constraint can be easily overcome by performing the computation with a local moving window versus storing the entire grid into system memory.The speed gain of FSS makes simulation of considerably denser grids more practical compared to the other methods.With this capability, our conjecture is that for those applications requiring interpolation, less expensive bilinear or nearest neighbor interpolation could be adequate for very dense grids versus more expensive Kriging in coarser grids.The variogram (correlation) model is constrained to an exponential function with FSS, which makes it less flexible than the GSTAT (Sequential Gaussian Simulation) methods.However, the tradeoff in speed gain and simplicity of implementation offers practical and useful advantages to motivate a potentially broader community of users.

Extension of FSS to a Multivariate Gaussian 2D Random Field
The FSS core grid generation equation, Equation (3), can be extended from a scalar error to a (multivariate) error vector over a 2D grid in a straightforward manner.The more general case is presented directly below, with special but practical subcases presented in following subsections that include simpler notation., the Hadamard product (term by term product) of two matrices, where and , correspond to matrix row , column .
Note that the above constraint that is a positive definite matrix is not satisfied for all possible combinations of , , and desired (valid) steady state error covariance , in which case Equation (9) and its statistics are no longer valid.consist of each element of multiplied by the scalar value , .

Diagonal Covariance Subcase
Another special, but practical, subcase of Equation ( 9) is when the specified is a diagonal matrix.This allows for any values of and , i.e., different specifiable spatial correlations for each of the two directions for each of the error components.Additionally, of course, this allows for different variances specified along the diagonal elements of ; also, the constraint on is always satisfied.Note that this special case is simply equivalent to the scalar case for each of the components applied independently.
The resultant system of equations are identical to Equation ( 9) except that we have the following diagonal matrices: (12) where .
Note that each component of has its own spatial correlation function with specifiable distance constants.

Example Realizations
This section presents quiver plots of multivariate error over a subset of a 2D final grid generated using the sequential algorithm discussed in Section 6.The multivariate error corresponds to a two-dimensional vector ( ).The corresponding covariance matrix is diagonal with a common variance for the two components of error for convenience, i.e., .Similarly, spatial correlations corresponding to a specific spatial direction are common for convenience, i.e., and .Figure 13 (automatically scaled) corresponds to m, and and .Figure 14 (automatically scaled) corresponds to m, and and .

Multi-Grid Point Covariance Matrix
For the special case of a diagonal covariance matrix , the corresponding covariance matrix for a collection of at arbitrary grid points has a convenient and valid representation: (13) where is an vector such that and the , , correspond to an ordered list of the grid point locations.Also, corresponds to the matrix Hadamard (element by element) product, the diagonal matrix , and corresponds to the spatial correlation function associated with component of .

General Case with Constraint Enforced
Referring back to the general case of Section 6, the following presents two examples for .Assume that the two components of error correspond to -error and -error for specificity.From Equation ( 9), can have any combination of values such that each is within the positive interval (0,1) and , a function of the desired and the , is a symmetric and positive definite matrix.

Figure 13.
Realization of two-dimensional multivariate errors over a 2D grid: high spatial correlation in the grid's or y-direction and high spatial correlation in the grid's or -direction.
Subcase 1: Assume that and , to narrow down the possible combinations; therefore, (14) which is always positive definite for any and .
Subcase 2: Assume that and , therefore, (15) which is positive definite if , where The left portion of Figure 15 plots the upper and lower bounds for given the desired value of and assuming that ; the right side assuming that .As can be seen from the above, the larger the absolute value of the correlation between the error components x and y, the closer must be to .Note that any multi-grid point covariance matrix for this general case must be assembled "term-by-term" using the a priori statistics presented in Equation ( 9), i.e., there is no convenient functional form for the cross-covariances in the multi-grid point covariance matrix similar to those presented for the special case of common spatial correlation among the components of and the special case of a diagonal covariance matrix presented earlier.

Extension of FSS to a Non-Homogeneous 2D Random Field
This section of the paper extends the FSS of a scalar homogeneous Gaussian 2D random field to a scalar non-homogeneous Gaussian 2D random field.In particular, the specified values for , , and (variance and spatial correlation parameters) corresponding to are either explicitly or implicitly a function of grid location .There is no one "right way" to do the extension.Two general methods are presented below, each practical but with different characteristics regarding the form of non-homogeneity represented.Each method can also compute a corresponding multi-grid point covariance matrix, necessary for many applications as discussed earlier in Section 3.For one method, this covariance matrix is exact, for the other, an approximation.The best technique, when both non-homogeneity characteristics and possible multi-grid point covariance matrix approximations are taken into account, is application dependent.(The methods presented below can also be extended in a straightforward manner to multivariate non-homogeneous Gaussian 2D random fields.)

Method 1: Convex Combination
The core grid generation equation and corresponding sequential algorithm for scalar errors (Sections 2.2 and 2.3) is simply exercised different times, either sequentially with the results saved temporally, or in parallel in order to save storage.(Of course, the grid size and spacing remains constant each time.)The number of times is typically two, i.e., . Each uses a different set of specified , , and .Thus, after the above is performed there are grids, each homogenous and in accordance with the , , and specified for use with that particular grid.Each grid is uncorrelated with the others.
The grids of , designated for , are then combined based on a convex combination into a final grid of .That is, at each location in the grid: , where and ( The specification of the values, also symbolized as for convenience, can be as simple or as complicated as appropriate over the locations across the grid.However, their recommended values are in accordance with the following: Let us assume that is to be exclusively the value across the various in Region i of the grid; hence, in this region, all .In addition, let us define Region i-j as a "buffer region" from Region i into Region j.In this buffer region, varies linearly from 1-0 corresponding to the at the start to the end of the buffer region, respectively.Furthermore, of course, throughout Region i-j.Finally, for all locations in .See Figure 16 as an example for .Note that the width of the buffer region Region i-j should be at least twice the maximum of the corresponding spatial distance constants associated with and , expressed in grid unit. This ensures reasonable spatial correlation across the buffer region.If there were no buffer region, the spatial correlation between two points, one anywhere in Region i and the other anywhere in Region j, would be 0, i.e., there would be an unwanted abrupt change across the boundary of the two regions.Assume that different scalar in the (final) grid are of interest regarding a corresponding multi-grid point covariance matrix.Each of these corresponds to their own unique location in the grid, and are ordered in a known fashion sequentially for , and placed into an vector , where .Such a vector can also be defined for the same ordered locations for each of the different realizations as , .Therefore, based on Equation (17) we have: , where is an diagonal matrix for with the appropriate values of down its diagonal.For example, if the first component of corresponds to grid location , the first diagonal component of equals .
Let us represent the corresponding multi grid point covariance matrix for the as .(See Section 3 for how this matrix is computed given the corresponding , , and .)The multi grid point covariance matrix for is computed as follows: (19) A nice feature of Method 1 is that the above representation for the multi-grid point covariance matrix is exact.also corresponds to the following: (20) where is an vector such that and the , , correspond to the ordered list of the grid point locations; corresponds to the explicit correlation between two such points; the matrix entries "_" indicate symmetry., etc.

Example Realizations
The following non-homogeneous realization combines two homogeneous realizations with and , respectively.Region 1 of the displayed portion of the final grid consists of =1-10, Region 1-2 = 11-30, and Region 2 = 31-60.(For each region, the corresponding = 1-60.)Use of the typical assignment scheme for the values of (and was employed.Figure 17 below presents the results.The same experiment was performed (but different realization) except that there was no Region 1-2, i.e., a non-typical scheme.Figure 18 below presents the results.

Method 2: Functional Variation of a Priori Statistics
With the second method, the core grid generation equation and corresponding sequential algorithm for scalar errors (Sections 2.2 and 2.3) is implemented only once, but modified as follows: where is a random sample of Gaussian white noise distributed , and where .Also, , and , that is, the spatial distance constants can be considered a function of as well.
Figure 18.An abrupt transition between Region 1 and Region 2, each with their own specified a priori statistics.
As indicated above, the values for , , and , and hence , are a function of the grid location .In addition, for Method 2, they are determined by the bilinear interpolation of such specified values over a less-dense grid overlaying the grid of errors to be generated.For example, if the 2D grid of errors to be generated is 900 × 1000, the grid for interpolation might be an evenly spaced 4 × 3 parameter grid overlying the denser grid.Each of the corresponding 12 parameter grid points contains the specified values for , , and for the corresponding local region around the parameter grid point.Note that is a function of the interpolated values of , , and ; hence, is also recalculated in Equation (21) for every grid location .See Figure 19 for a graphical representation of the interpolation parameter grid.Each interpolation parameter grid location contains a unique set of values for , , and .Also, the spacing between interpolation parameter grid points should be at least twice the maximum of the corresponding spatial distance constants associated with that grid point and the other interpolation parameter grid points immediately surrounding it, expressed in grid unit.This ensures that both the desired and the computed approximation of the a priori statistics corresponding to the across the dense grid are approximately met and reasonably reliable (see Section 7.2.2),respectively.(This also assumes that the appropriate buffer relative to the "final" grid is included as well-see Section 2.3.2) Once defined appropriately, Equation ( 21) is then implemented via a direct counterpart to the algorithm described in Section 2.3.The latter simply utilizes the appropriate , , and values that vary with location.

Example Realizations
The following examples correspond to a parameter interpolation grid overlaying a 2D grid.The results corresponding to a displayed portion of the final grid are presented.All 49 sets of parameters were identical except for four sets corresponding to an interior rectangle near the center of the final grid.Let us term the 45 common sets as Group 1 and the other four common sets as Group 2. In Figure 20 below, the Group 1 set contain values m, .Group 2 sets contain values m, .In Figure 21 below, the Group 1 sets contain m, .Group 2 sets contain m, .

Statistics and Multi-Grid Point Covariance Matrix
Corresponding a priori statistics are no longer straightforward for this method, but can be approximated.In particular, corresponding to a specific location is the corresponding bilinear interpolated value.The spatial correlation function corresponding to different locations is the average of spatial correlation functions, each corresponding to the bilinear interpolated values for and for that location.Of course, these statistics reflect non-homogeneity, i.e., are a function of the specific locations of interest.The corresponding approximation for the multi-grid point covariance matrix corresponding to scalar errors at different grid locations is represented as follows:    Because the average of a collection of strictly positive definite correlation functions (spdcfs) is an spdcf itself, the above is guaranteed a valid covariance matrix regardless of the fact that the various can vary in value (Reference [7]).
Note that if the grid points consist of widely spaced subgroups of points such that the scalar error at any grid point in one subgroup has (approximated) low correlation (e.g., less than 0.1) with the scalar error at any grid point in any other subgroup, a higher fidelity representation for the multi-grid point covariance matrix can be achieved as follows: Use the representation in Equation ( 22) to compute a "sub-multi-grid point covariance matrix" for each subgroup of grid points, and then combine them into the (final) multi-grid point covariance matrix by placing (in order) each sub-multi-grid point covariance matrix down the block diagonals with values of zero for all off-diagonal (cross-covariance) blocks.
Finally, to generate a multi-grid point covariance matrix for a non-homogeneous multivariate instead of that for a non-homogeneous scalar , the same general procedure presented in this section can be extended in a straightforward manner using methods described in References [6,7].

Summary and Conclusions
Practical methods for the sequential generation of two-dimensional random fields were presented, and their corresponding multivariate covariance matrices derived.The corresponding methods were based on FSS, which was also compared to Sequential Gaussian Simulation and other approaches.Although less general, FSS was shown to be clearly superior in terms of speed and simplicity, primarily due to assumed separable exponential spatial correlation and simple ordered generation over an evenly spaced grid.FSS methods presented in the paper are applicable to performance assessment and tuning of geospatial applications in a simulation environment, as well as near-real time display of the effect of errors on applications by the applications themselves.

A.2.1. Detailed Derivations
The following derives Equation (A3) by the statistical properties of Equation (A1): ; ; The above utilized the following: by definition, , , and ; by the properties of a geometric series, , which is applicable since and in the above.Define and for convenience.
The latter representation is due to the fact that .Thus, .Similarly, .By the same procedure, it is also follows that .
is the best estimate of the realization at the appropriate horizontal location.is the corresponding simulated value, where corresponds to a mean-zero Gaussian random number with variance , i.e., the variance of the Kriging solution relative to the value of the true realization.
(More details can be found in the literature for simple Kriging [1,3], conditional and unconditional simulation [1], sequential simulation [14], Sequential Gaussian Simulation [4], and Gaussian probability distributions, the conditional mean, and the conditional variance [15].) Let us now assume specific grid point locations of interest and a Gaussian two-dimensional random field with a priori standard deviation and separable exponential spatial correlation with correlation between adjacent grid point locations and .A representative set of 14 grid locations ( ) are represented graphically as follows:  3), corresponds to and to .)We now implement the Kriging equations (B1) for Sequential Gaussian Simulation using the value of generated by FSS and detailed in the previous paragraph.In addition, the same additive Gaussian random number will be added to the Kriging solution for per the Sequential Gaussian Simulation procedure since, as will be demonstrated below, .In support of the Kriging solution, the a priori cross-covariance matrix between and and the a priori covariance matrix for are computed in accordance with the assumed statistics ( ) of the random field presented previously.Correspondingly, both and are full (no zero elements), but the product is only non-zero for the three elements of which correspond to the nearest three grid locations 8, 9, and 13 to the point to be simulated.(It also follows that a pre-computed "compressed" 1 × 3 version of can actually be used as common Kriging weights for the realizations at the three nearest grid points.), identical to that generated using the FSS method of this paper.In addition, , which is equal to .These equivalences and the use of only the nearest three grid points were enabled due to both the separable exponential spatial correlation and the regular grid of realizations generated in a simple, ordered fashion.(Thus, for example, if grid point #8 were moved +0.5 grid units in the y-direction, there would be 7 instead of 3 non-zero scalar weights.) The above was an arbitrary, but specific, example.A formal analytic proof that the Kriging weight row vector has only non-zero weights -rs, s, and r corresponding to the three nearest grid points is relatively easy and done by direct verification that , where is a row vector consisting of all zeroes except for the non-zero scalar weights -rs, s, and r at the appropriate locations corresponding to the three nearest grid points (see Figure B1).Similarly, in order that , must equal ), which is easily verified by direct evaluation of .Thus, this appendix has both demonstrated and proven that FSS is equivalent to Gaussian Sequential Simulation under appropriate circumstances.

Figure 1 .
Figure 1.An example of the separable spatial correlation function; in the plot of , and have signed values.

Figure 4
presents the realization results of Example 1 based on Approach 1.Note that the remaining realization examples in this paper use Approach 1 as well, as it is most convenient.

Figure 4 .
Figure 4. Example 1-Realization of -error with high spatial correlation between adjacent grid points.

Figure 5
Figure 5 corresponds to Example 2, a new realization with the same m, but with (thus m).Figure 6 corresponds to Example 3, a new realization with the same m, but with (thusm).As expected, the above realizations over portions of the grid do not have a mean -error of zero nor a standard deviation about the mean of 10 m.However, when sample statistics are computed over numerous realizations, the corresponding mean and standard deviation approach 0 m and 10 m, respectively, matching the common a priori statistics used to generate the realizations.Finally, Figure7below presents Example 4, a new realization with the same m, but with and (thus, ), i.e., different spatial correlations in the two directions.

Figure 6
Figure 5 corresponds to Example 2, a new realization with the same m, but with (thus m).Figure 6 corresponds to Example 3, a new realization with the same m, but with (thusm).As expected, the above realizations over portions of the grid do not have a mean -error of zero nor a standard deviation about the mean of 10 m.However, when sample statistics are computed over numerous realizations, the corresponding mean and standard deviation approach 0 m and 10 m, respectively, matching the common a priori statistics used to generate the realizations.Finally, Figure7below presents Example 4, a new realization with the same m, but with and (thus, ), i.e., different spatial correlations in the two directions.

Figure 5 .
Figure 5. Example 2-Realization of -error with low spatial correlation between adjacent grid points.

Figure 6 .
Figure 6.Example 3-Realization of -error with very high spatial correlation between adjacent grid points.

Figure 8 .
Figure 8. Sample statistics corresponding to 1000 × 1000 grid with different a priori correlations.

Figure 9 .
Figure 9. Sample statistics corresponding to 1000 × 1000 grid with the same a priori correlations for the x and y directions, evaluated across three different directions.

Figure 10 .
Figure 10.Semivariograms corresponding to 200 × 200 grid with the same a priori correlations for the x and y directions, evaluated across the y-direction for five different realizations.

Figure 10
Figure10corresponds to a set of five independent realizations, this time over a much smaller 200 × 200 grid.A priori statistics were identical to those corresponding to Figure9except that .Also, the sample and theoretical statistics computed correspond to the semivariogram,

Figure 11 .
Figure 11.Bilinear interpolation of four surrounding grid points.

Figure 12 .
Figure 12.Time comparison among methods for unconditional simulation of a scalar random fieldover a 2D square grid (Note that FSS is cutoff at ~15 s due to reaching the system memory limit).
valid (symmetric and positive definite) covariance matrix which satisfies the following:

Figure 14 .
Figure 14.Realization of two-dimensional multivariate errors over a 2D grid: high spatial correlation in the grid's or y-direction and lower spatial correlation in the grid's or -direction.

Figure 15 .
Figure 15.Flow-down of constraints to spatial correlation bounds.

Figure 16 .
Figure 16.Example of region layout over a grid.

Figure 17 .
Figure 17.A smooth transition from Region 1 to Region 2, each with their own specified a priori statistics.

Figure 19 .
Figure 19.An example of an Interpolation Parameter Grid.
ordered list of the grid point locations.