A New Non-stationary High-order Spatial Sequential Simulation Method

A new non-stationary, high-order sequential simulation method is presented herein, aiming to accommodate complex curvilinear patterns when modelling non-Gaussian, spatially distributed and variant attributes of natural phenomena. The proposed approach employs spatial templates, training images and a set of sample data. At each step of a multi-grid approach, a template consisting of several data points and a simulation node located in the center of the grid is selected. To account for the non-stationarity exhibited in the samples, the data events decided by the conditioning data are utilized to calibrate the importance of the related replicates. Sliding the template over the training image generates a set of training patterns, and for each pattern a weight is calculated. The weight value of each training pattern is determined by a similarity measure defined herein, which is calculated between the data event of the training pattern and that of the simulation pattern. This results in a non-stationary spatial distribution of the weight values for the training patterns. The proposed new similarity measure is constructed from the high-order statistics of data events from the available data set, when compared to their corresponding training patterns. In addition, this new high-order statistics measure allows for the effective detection of similar patterns in different orientations, as these high-order statistics conform to the commutativity property. The proposed method is robust against the addition of more training images due to its non-stationary aspect; it only uses replicates from the pattern database with the most similar local high-order statistics to simulate each node. Examples demonstrate the key aspects of the method.

similar patterns in different orientations, as these high-order statistics conform to the

Introduction
The uncertainty of a spatially distributed and varying geological attribute can be quantified by analyzing the variation in this attribute within a set of geostatistical or stochastic simulations. Second-order spatial simulation methods (Journel and Huijbregts 1978;David 1988;Goovaerts 1997;Chiles and Delfiner 2012) have been used to quantify spatial uncertainty. Since the early 1990s, a new multiple-point statistics (MPS) spatial simulation framework has been developed to overcome the limitations of previous approaches (Guardiano and Srivastava 1993;Gómez-Hernández and Srivastava 2021;Strebelle 2002Strebelle , 2021Journel 2005;Liu et al. 2006;Arpat and Caers 2007;Hu and Chugunova 2008;de Vries et al., 2008;Mariethoz et al. 2010;Honarkhah and Caers 2010;De Iaco and Maggio 2011;Stien and Kolbjørnsen 2011;Chatterjee et al. 2012;Lochbühler et al. 2013;Rezaee et al. 2013;Toftaker and Tjelmeland 2013;Mustapha et al. 2014;Strebelle and Cavelius 2014;Zhang et al. 2017). MPS methods are able to capture and produce more complex spatial patterns than second-order spatial simulation methods by using multiple-point spatial templates that replace the well-established two-point spatial statistics. In the MPS framework, sample data are used as the conditioning data, while the multiple-point spatial relations of the simulations are captured from patterns extracted from a training image (TI) or analogue. As a result, MPS simulations exhibit the statistics of the TI instead of the data, which raises the issue of statistical conflicts between the samples and the TI in circumstances in which sample data are relatively abundant, such as in mining applications (Osterholt and Dimitrakopoulos 2007;Goodfellow et al. 2012). High-order spatial simulation methods (HOSIM) have been developed to address this issue and to provide a new geostatistical simulation framework that deals with complex spatial patterns. HOSIM does not require distributional assumptions and is mathematically consistent (Dimitrakopoulos et al. 2010;Dimitrakopoulos 2010a,b, 2011a;Mustapha et al. 2011Mustapha et al. , 2014Tamayo-Mas et al. 2016;Minniakhmetov and Dimitrakopoulos 2016;Minniakhmetov et al. 2018;Yao et al. 2018). It should be noted that in HOSIM methods, the conditional probability distribution function (CPDF) is modelled by a set of orthogonal polynomials, Legendre polynomials in particular (Abramowitz 1974). Under the stationarity assumption, high-order spatial cumulants (Dimitrakopoulos et al. 2010) are inferred by averaging the same-order cumulants of the samples extracted from either the data set alone, if sufficient data are provided, or both data and TI in the case of sparse data. The high-order statistics enable the CPDF model to capture the complex spatial structures and connectivity of the high values within the data and to later reproduce them via the simulation process (Mustapha and Dimitrakopoulos 2011a, b;Minniakhmetov et al. 2018;de Carvalho et al. 2019). Stationarity is a major assumption made by the commonly employed spatial simulation methods and may not always be easily accommodated in applications. The traditional approach assumes that a mean function (drift) exists as a linear or polynomial trend. The so-called universal kriging (David 1977) considering the drift, however, may be difficult to apply to some cases (Cressie 1986). More recent research adopts the local probability distribution to model the non-stationary distribution, but only second-order statistics are considered (Machuca-Mory and Deutsch 2013). With respect to high-order simulations, the same assumption also results in the use of all patterns extracted from data and the TI for the estimation of spatial cumulants, leading to the slow convergence of the orthogonal polynomials used for modelling CPDFs. In practice, this calculation may result in numerical instabilities (Boyd and Ong 2011). Thus, a non-stationary approach is used herein to address these issues. For each simulation, a set of patterns are chosen from the data and TI based on their similarities to the data event of the simulation pattern. Then, only the chosen patterns are used for inferring the high-order statistics of the simulation pattern. Similar patterns could be generated from a simple CPDF model with faster convergence using a lower number of polynomial terms.
The high-order sequential simulation method introduced here first finds a data event of a fixed size n in the immediate neighbourhood of each simulation node, and a set of n +1 high-order spatial statistics is then calculated from the data event. A novel similarity measure is then introduced and utilized to compare the extracted statistics from the simulation grid to those extracted from all of the patterns in the TI; the statistics with the highest similarity score are then used for the inference of the high-order statistics at that node. This is the non-stationary aspect of the method, which is the key difference between this method and previous high-order simulation methods Dimitrakopoulos 2010b, 2011b;Minniakhmetov et al. 2018, Minniakhmetov andYao et al. 2018Yao et al. , 2020Yao et al. , 2021. The proposed method avoids an excessive number of patterns extracted from the TI due to the selectivity of the method in terms of choosing the number of related patterns. In addition, the method is able to efficiently identify patterns with similar statistics in different orientations without reducing the computational efficiency of the algorithm.
In the following sections, the proposed method and corresponding algorithm are first detailed. Then, examples and comparisons are presented, and conclusions follow.

Spatial Random Field
Given a probability space ( , F, P), a function Z (X ) : D ⊂ R N → R is a realvalued spatial random field given that for any integer number n ≥ 0 and ∀r ∈ R n+1 , a subset A r {Z (X ) ≤ r } ∈ F exists and its probability is defined by P, where D is the spatial domain in a two-or three-dimensional space R N (N 2, 3) and X ⊂ D is a set of n + 1 points in this domain. Note that X {x, x 1 , . . . , x n } and Z (X ) {Z (x), Z (x 1 ), . . . , Z (x n )} and thus the binary operation ≤ in A r {Z (X ) ≤ r } is true only if all elements on the left-hand side are less than or equal to the elements on the right hand side of the equation.

Template, Pattern, Data Event and Neighbourhood
Consider a random field Z (X ) with a fixed structure in space and a spatial template T that is characterized by a set of lag vectors (Mustapha and Dimitrakopoulos 2011b), that is, T {h 1 , . . . , h n }. The template connects a specific location x to each of the positions in the neighbourhood N x x + T {x 1 , . . . , x n }. Without loss of generality, hereafter the real-valued spatial random field associated to the location x and its neighbourhood N x is denoted by Z (x, N x ), and its realization is denoted by T represents the realization of the field at the position of the neighbourhood nodes.

Invariant High-Order Statistics
Vieta's formula (Funkhouser 1930) presents a relation between the coefficients and the roots of a polynomial. For a polynomial of degree n, p(X ) X n + c 1 u 1 X n−1 + · · · + c n u n , with m ∈ {1, . . . , n} and the normalization factor c m ( n 0 ptm )σ m , where σ is the standard deviation of the x in the TI. Furthermore, a recurrence relation could present the coefficients of a polynomial q(X ) (X − x 0 ) p(X ) X n+1 + c 1 u 1 X n + · · · + c n+1 u n+1 , which has one extra root x 0 , that is, {x 0 , x 1 , . . . , x n }. Then, it follows that Vieta's formula transforms the input vector of the roots of the polynomial [x 1 , . . . , x n ] T into the vector of the coefficients [u 1 , . . . , u n ] T , with the advantage of being invariant to the ordering of the domain, given that the polynomial is invariant under a re-ordering of its roots. Herein, Eq. (2) is used to transform either the patterns or the data events of the simulation and TI, such that (4)

High-Order Statistics Distance Vector
Given two sets of data events d N xs and d N xs of order n, in order to develop an L2-norm distance measure, one must compare all possible ordering of these two data events, which results in n × (n!) number of operations. Hence, for a pattern of order n 5, the number of operations becomes 600 per pattern, which is computationally expensive and can only operate on small TIs. In a new approach (Abolhassani et al. 2017), the L2-norm distance is calculated for the high-order statistics vectors of the data events where the high-order statistics of the data event U s are calculated from Eq. (2). It is worth mentioning that the number of operations for this new distance measure is dramatically reduced to n × 2 n−1 per simulation node versus n × (n)! in L2-norm. For n 5, #op(L2 − norm) 600 versus #op(high − ord) 80.

Modelling the Spatial Non-Stationary Joint Cumulative Distribution Function (CDF) as a Finite Sum of Disjoint CDFs
Given a probability space, the real-valued spatial random field Z (x, N x ) is fully characterized by its non-stationary joint cumulative probability distribution function F(x, Z (x, N x )). In this work, this function is modelled by the statistical ensemble (Landau and Lifshitz 1980) given by which is a sum over m mutually exclusive stationary CDFs, F i s, with a set of nonstationary mixing coefficients, known as a state set with binary variables φ i (x 0 ) ∈ {0, 1}. At any position x 0 , the value of only one of the coefficients is 1, and the rest are zero, that is, φ k 1 and φ i 0 for i k.

Mutually Exclusive CDFs
On the grid of an input image with a fixed template T , Eq. (6) is considered to have a mutually exclusive set of stationary CDFs only if it satisfies the condition ⇒ The high -order statistics of two patterns match where D(.) is the high-order statistics distance vector calculated from Eq. (5) to compare two sets of high-order statistics vectors of two data events d N x 1 and d N x 2 . In addition, U Z (x 1 , N x 1 ) and U Z (x 2 , N x 2 ) are the high-order statistics of two patterns, introduced in Sect. 2.3, with elements calculated from relation (2). According to Eq. (7), if the statistics of the data event of two patterns are close, then their CDFs are equal to an identical stationary CDF, defined as F k (Z (x, N x )), and, and further, that their high-order statistics are the same. This also implies that φ k (x 1 ) φ k (x 2 ) 1. A set of mutually exclusive CDFs are used as the basis functions to span a more complex non-stationary CDF using Eq. (6). At each position on the grid, only one basis function is responsible for the simulation of the node; that is, each basis function F k (Z (x, N x )) simulates a set of nodes D k on the grid. The D k s are mutually disjoint, and their union is the entire grid D ∪ m k 1 D k . As a result, for each simulation node, the entire TI is searched for nodes with similar statistics. Considering Eq. (7), the resulting set of nodes are generated from one of the basis CDFs and used for estimating the high-order statistics of that basis CDF, and the simulation node value is then sampled from the CDF.

High-Order Transformation-Invariant Sequential Simulation
The goal is to simulate the non-stationary, real-valued random field {Z (x s , N x s )|x s ∈ s } at the position of all simulation nodes s with a fixed template T. A sequential multi-grid process is used when provided with a sparse set of data In the hierarchy of the sequential multi-grid simulation, the coarse templates are firstly adopted to complete the first-stage simulation. As more nodes are simulated, finer-scale templates are used to search for the conditioning data. The order of the sequential simulation is shown in Fig. 1.
At each sequence, the path is chosen randomly and saved into a vector containing the indices of the visiting nodes I s . Each successive random variable Z (x s , N x s ) at node {x s |s ∈ I s } is conditioned to a data event d N xs and the neighbourhood nodes that are selected from the set of previously simulated nodes and the sample data (Goovaerts 1997). A template T {h 1 , . . . , h n } is formed spatially by connecting each data event node to the simulation node. The first goal is to simulate the high-order statistics of the pattern, that is, U Z (x s , N xs ) , and then to sample the simulation node Z x s , N x s , Fig. 1 The hierarchy of the sequential multi-grid simulation. The order of the simulation starting from the blue nodes (sample data) and the red nodes (simulated first), followed by black nodes, orange nodes and eventually the green nodes. The simulation of each stage is conditioned on the previous simulated nodes given the CPDF P(U Z (x s , N xs ) |U d Nx s , U Z (x t , N x t ) ). In this case, the U (.) function refers to the high-order statistics vector with the elements calculated from Eq. (2). This probability is complex and cannot be simulated directly unless it is represented by a model with a set of parameters θ ∈ that are independent from the data event and TI. The parameter θ is optimized to express the data event and TI. Thus, this probability can be decomposed using the total probability rule into An estimation of the argument of the integral further simplifies Eq. (8). Figure 2 represents the term P(θ |U d Nx s , U Z (x t , N x t ) ) as a function of θ . The contribution of this function is negligible except for a narrow band near an optimal value for the parameter's maximum a posteriori estimation (MAP) θ MAP . Consequently, Eq. (8) can be estimated as In Eq. (9), C MAP is the normalization factor to ensure the validity of the probability. Equation (9) implies that the probability distribution function of the node x s can be calculated if θ MAP is known. To estimate θ MAP based on (Fig. 2), θ θ MAP when P(θ |U d Nx s , U Z (x t , N x t ) ) is maximized for θ ∈ ; that is, Using Bayes' rule, the right-hand side of the equation is extended into In this case, one assumes a uniform prior in the parameter space , with the marginal probability in TI P(U Z (x t , N x t ) ) remaining independent from θ . Hence, the θ MAP is equivalent to the maximum likelihood estimation θ MLE , and Eq. (11) simplifies into Each node in TI is only conditioned on its neighbours N x t and the parameter set θ . Hence, the joint distribution in Eq. (12) can be further decomposed into A probability is maximized if the logarithm of that probability is maximized.
Furthermore, at the maximum, the derivative with respect to θ should be zero; that is, Equation (15) is solved for an optimal solution θ MLE by modelling the P(U Z (x t , N x t ) |θ , U d Nx s , U d Nx t ) with an exponential family in Sect. 2.8.

High-Order Simulation Model
The exponential family is used here to model the likelihood function in Eq. (15) with the parameter set θ [θ 1 , . . . , θ n+1 ] T as In Eq. (16), W is introduced as a weight function based on the similarity measure of the data event d N xs and d N x t . It ensures that the TI patterns with more similar data events contribute more to the likelihood function. W has a probability distribution over [0, 1]; W ML is the maximum likely similarity that is estimated by where D(., .) is given by Eq. (5) as the high-order statistics distance vector, a distance measure between two data events, and d N x t is the mean square distance vector of the data event over the TI. Note that the weight W is defined to be a positive function that decreases as the distance between the two data events increases. σ 0 2 in Eq. (16) is to ensure that a proper probability is established for Eq. (16).
Substituting Eq. (16) into Eq. (15), the solution of the θ MLE is Applying the recurrence relation (3), Eq. (20) is rewritten in terms of the simulation data event U d Nx s and the simulation node Z (x s ) (the extra root in Vieta's formula (3)); that is, where u m can be calculated with Eq.
(2) and c m ( n 0 ptm )σ m , in which σ is the standard deviation of the x in the TI. Eventually, the simulation node Z (x s ) can be sampled from the distribution in Eq. (21).

Algorithm
1. Define a random path s to visit all positions on the simulation grid s .

Examples and Comparisons
Applied aspects of the simulation method presented in the previous section are examined in this section using two-dimensional horizontal layers from the Stanford V fluvial reservoir (Mao and Journel 1999). Selected layers represent exhaustive images that are then sampled to generate data sets subsequently used to generate simulations with the proposed new non-stationary high-order spatial transformation-invariant simulation method, HOSTSIM. The statistics of the simulated realizations, namely, histogram, variogram map and third-order L-shaped cumulant map, are then compared to the ones calculated from the initial data. In addition, the simulations generated by the HOSTSIM method are compared to those generated by the well-known FILTERSIM multiple-point simulation method ) through its implementation available in the public domain (Remy et al. 2009). Figure 3 shows two exhaustive images used for the simulations E 1 and E 2 ; each one is a horizontal section of the complete data with size of 100 × 100m 2 . Two sets of data are then selected from each section E 1 and E 2 , one sampled at 625 and the other at 156 locations, respectively, resulting in four different data sets, namely DS 1 (625 points from E 1 ), DS 2 (156 points from E 1 ), DS 3 (625 points from E 2 ) and DS 4 (156 points from E 2 ), as shown in Fig. 4. The TI used for all simulations is a different section selected from the Stanford V fluvial reservoir and is shown in Fig. 5.
In each example, one of the data sets in Fig. 4 and the TI in Fig. 5 are used as input for simulations, and two realizations are generated by each method noted above. The exhaustive images, the input data and the realizations generated by HOSTSIM and FILTERSIM are illustrated in Figs. 6, 7, 8 and 9 for data samples DS 1 through DS 4 , along with their respective validation graphs. The validation graphs are generated for both the input data and the HOSTSIM simulation and include the histogram, variogram map and third-order L-shaped spatial cumulant map. Note that the secondorder variogram map of an image represents the variogram for each lag direction, that is, r ∈ [0, 65] and θ ∈ 0, π 30 , . . . , π , in a polar coordinate system. The third-order L-shaped spatial cumulant of an image is generated by using an L-shaped template, which represents two orthogonal lag vectors from a common point. For all possible pairs of horizontal and vertical lags, the template is shifted over the image to extract the third-order cumulant by averaging the values calculated for each pattern within the image Dimitrakopoulos 2010b, 2011a). Figures 6 and 7 show the results generated from DS 1 and DS 3 with 625 data samples. There are low-contrast regions in the HOSTSIM realizations (light blue/green area seen toward the bottom of E 1 in Fig. 6, and in the bottom-right and top-center of the E 2 realization in Fig. 7) in both realizations. This behaviour is expected, as a denser set of data is available in these regions, and HOSTSIM can match similar patterns within the TI and better simulate the corresponding nodes. This is also the case in the high-value areas on the simulated realization (red-coloured regions, or channels). In general, both HOSTSIM and FILTERSIM perform well in reproducing the horizontal channels, which are along the preferential direction in both the exhaustive image and the denser data samples. However, HOSTSIM shows better performance than FILTERSIM in reproducing the connectivity along the vertical channels, as can be seen from the major vertical channels in the left part of the exhaustive image and  The E 1 exhaustive image is in top row (left), and the input image DS 2 is in the second row. The two HOSTSIM and FILTERSIM simulations are in the second and third columns (above), followed by the histogram, variogram map and third-order spatial cumulant map (below) Fig. 9 Example 2: The E 2 exhaustive image is shown in the top row (left), and the input image DS 4 is provided in the second row. The two HOSTSIM and FILTERSIM simulations are in the second and third columns (above), followed by the histogram, variogram map and third-order spatial cumulant map (below) the related realizations. In this situation, the concept of weighted high-order spatial statistics plays an important role; accordingly, HOSTSIM demonstrates an advantage because the samples along the vertical channels are given more weight according to their similarity to the sample data events even when they are less frequent. Although there are some cases where the connectivity of the channels is not well generated, for example in the top-right part of both realizations in Fig. 7, the narrow channel is not well represented in the data, and the HOSTSIM cannot accurately estimate the highorder statistic of the data event to match appropriate patterns from either the data or the TI. The histogram and variogram graphs of the HOSTSIM simulations are similar to those calculated for DS 1 and DS 3 in both Figs. 6 and 7. Similarity of the third-order spatial cumulants can also be observed from Figs. 6 and 7.

Example 2
The HOSTSIM realizations are presented in Figs. 8 and 9. As the number of input data are reduced to 156 samples, the quality of the results is reduced (relative to Example 1). The low-contrast regions are, however, simulated in these examples (light blue/green areas). Some parts of the channels are not generated in both realizations (the center-left part of the channel in Fig. 8 and the center-left and bottom-left part of the channels in Fig. 9) due to a lack of high-value samples in those regions, which are necessary in order to represent narrow channels. On the other hand, some channels are merged in the HOSTSIM realizations (center-right part of both Figs. 8 and 9), which is due to a lack of low-value samples between two channels in the input data samples. The validation graphs of the HOSTSIM simulations are similar to the ones generated for the input data samples in Figs. 8 and 9, although they are somewhat degraded due to the sparsity of the input data. While the HOSTSIM realizations generated for this example are somewhat degraded compared with the previous example, there is still better connectivity in the data than in the results from FILTERSIM; for example, the channel in the center-right part of Fig. 8 is partially generated in HOSTSIM realizations but not well presented in the FILTERSIM realizations. Low-contrast regions are better represented in the HOSTSIM realizations as well (see the bottom-left part of Fig. 8).
Regarding the reproduction of high-order spatial statistics, the third-order L-shaped spatial cumulants are well generated for the lags < 50 m. Considering that replicates with lags larger than 50 m are relatively few because the samples are taken from an area of 100 × 100 m 2 , given the sparsity of the samples, the inference of spatial cumulants within a short-scale range is more reliable in terms of comparisons.
To complete our analysis, two commonly used similarity measures in computer vision are applied here to compare the realizations generated by HOSTSIM and FIL-TERSIM for all cases. These measures are the peak signal-to-noise ratio (PSNR) and the structural similarity index (SSIM) (Salomon 2007). The PSNR is a point-to-point difference measure between two images in the logarithmic scale; in particular, the realization image is compared with the exhaustive image from point to point. Realizations that are more similar to the exhaustive image result in greater values, and where PSNR values approach infinity, the images are identical. Hence, this measure is useful for qualitatively comparing two realizations generated for an exhaustive image, as the

Conclusions
A new high-order, non-stationary stochastic simulation method was introduced in this paper. Given a spatially sparse set of data and a TI, this method simulates realizations of attributes of interest on a finer spatial grid. The method is designed to use only training patterns that respect the local high-order statistics of the simulation pattern by comparing the spatial high-order statistics of the data events of both the training and the target simulation patterns. The high-order statistics of the data event of each simulation pattern is first calculated and compared with high-order statistics of the data event calculated for the training patterns with a similar spatial template. The contribution of each training pattern in simulating a pattern is determined by a weight calculated based on the similarity of the high-order statistics of the training and simulation patterns. The calculated weights are designed to be invariant against ordering the data event nodes. The high-order statistics of the simulation pattern is then estimated by calculating the non-stationary maximum likelihood estimate over the training patterns by giving more weight to the most similar patterns, after which a sample is drawn from the distribution function fitted to the estimated high-order statistics. The use of the high-order statistics in this simulation method ensures the ability to capture and reproduce complex patterns in the generated simulations. Meanwhile, using the non-stationarity similarity measure ensures that only training patterns that respect the local statistics of the simulation are used in the estimation of the high-order spatial statistics of the pattern. Notably, (i) non-stationary weight reduces the numerical instabilities that may be encountered in stationary high-order simulation methods, and (ii) the transformation-invariant quality of the measure ensures that different orientations of each TI are considered for the simulation of the pattern. The latter allows for data-driven simulations even without the use of a TI. Two set of examples have been used for testing the proposed method. The first contains 625 data points and the second contains a very sparse set of 156 data points. The simulation results are calculated on a 100 × 100 spatial grid. The training patterns are extracted from a 100 × 100 grid image. Despite the sparsity of the data in both examples, the method is able to simulate the channels in the initial exhaustive image. This is due to the utilization of a non-stationary rotation invariant weighting for the contribution of the patterns based on the similarity of the high-order statistics of their data events. I (reference) and K (reconstructed) of size m × n are considered. Multiple smaller windows are generated from both images. Consider x and y representing two windows from I and K . SSIM is calculated by SSIM(x, y) 2μ x μ y + c 1 2σ xy + c 2 μ 2 x + μ 2 y + c 1 σ 2 x + σ 2 y + c 2 , ( 2 4 ) where μ x and μ y are the mean values over x and y windows. The parameters σ x , σ y and σ xy are the variance of x and y and the covariance between x and y, respectively. The values c 1 and c 2 are constants.