High-Order Spatial Simulation Using Legendre-Like Orthogonal Splines

High-order sequential simulation techniques for complex non-Gaussian spatially distributed variables have been developed over the last few years. The high-order simulation approach does not require any transformation of initial data and makes no assumptions about any probability distribution function, while it introduces complex spatial relations to the simulated realizations via high-order spatial statistics. This paper presents a new extension where a conditional probability density function (cpdf) is approximated using Legendre-like orthogonal splines. The coefficients of spline approximation are estimated using high-order spatial statistics inferred from the available sample data, additionally complemented by a training image. The advantages of using orthogonal splines with respect to the previously used Legendre polynomials include their ability to better approximate a multidimensional probability density function, reproduce the high-order spatial statistics, and provide a generalization of high-order simulations using Legendre polynomials. The performance of the new method is first tested with a completely known image and compared to both the high-order simulation approach using Legendre polynomials and the conventional sequential Gaussian simulation method. Then, an application in a gold deposit demonstrates the advantages of the proposed method in terms of the reproduction of histograms, variograms, and high-order spatial statistics, including connectivity measures. The C++ course code of the high-order simulation implementation presented herein, along with an example demonstrating its utilization, are provided online as supplementary material. Electronic supplementary material The online version of this article (10.1007/s11004-018-9741-2) contains supplementary material, which is available to authorized users.

herein, along with an example demonstrating its utilization, are provided online as supplementary material.
Keywords Stochastic simulation · Orthogonal splines · High-order spatial statistics · Non-Gaussian distribution · Spatial complexity

Introduction
Geostatistical simulations are used to quantify the uncertainty of spatially distributed attributes of interest describing mineral deposits, petroleum reservoirs, hydrogeological horizons, environmental contaminants, and other spatially variant natural phenomena. Since the 1990s, multiple-point spatial simulation (MPS) methods and variations (Guardiano and Srivastava 1993;Strebelle 2002;Journel 2005Journel , 2018Zhang et al. 2006;Arpat and Caers 2007;Chugunova and Hu 2008;de Vries et al. 2009;Straubhaar et al. 2011;De Iaco and Maggio 2011;Honarkhah 2011;Strebelle and Cavelius 2014;Chatterjee et al. 2012;Lochbühler et al. 2014;Mustapha et al. 2014;Rezaee et al. 2013;Toftaker and Tjelmeland 2013;Zhang et al. 2017;others) have been developed to advance the simulation methods beyond the past generation of second-order spatial statistics, which were typically based on Gaussian processes (e.g., Journel and Huijbregts 1978;David 1988;Goovaerts 1998;Chilès and Delfiner 1999). A limitation of MPS approaches is that they are largely algorithmic and do not consistently account for the high-order spatial relations in the available sample data. Patterns and complex spatial relations are derived from the so-termed training images (TIs), or geological analogues, rather than from sample data; this is a critical topic for relatively data-rich type applications, where data statistics have been shown to not be reproduced by simulated realizations based on MPS methods (Osterholt and Dimitrakopoulos 2018;Goodfellow et al. 2012). To address some of these limits, high-order simulation techniques for complex and non-Gaussian spatially distributed variables have also been developed Dimitrakopoulos 2010, 2011;Tamayo-Mas et al. 2016;Minniakhmetov and Dimitrakopoulos 2017a), based on generating conditional distributions through Legendre polynomials (Lebedev 1965) and high-order spatial cumulants. Yao et al. (2018) developed a new computational model to significantly reduce the computation cost of the method. The high-order simulation approach does not require any transformation of initial data and makes no assumptions about the related probability distribution function. The approach reproduces high-order spatial statistics of sample data and a training image. The high-order spatial statistics are shown to capture directional multiple-point periodicity, connectivity of extreme values, and complex spatial architecture . However, polynomial approximations do not always converge to analytic functions (Runge 1901;Boyd and Ong 2009;Fornberg and Zuev 2007). In addition, the high-order polynomials are very sensitive to rounding errors for values near the endpoints of an approximation domain; therefore, even if the interpolants converge in theory, they will diverge rapidly when computed (Platte et al. 2011). This is critical for the simulation of extreme values, as they are located at the endpoints of an approximation domain. In an effort to improve upon the limitations of polynomial approximation, a spline approximation of complex multidimensional functions is considered here (Piegl 1989;Hughes et al. 2005;Ruiu et al. 2016).
Splines are piecewise-defined polynomial functions in which pieces are connected by some condition of smoothness. The places where these pieces meet are called knots, and two adjacent knots form a knot interval, hereafter referenced simply as an interval. Knot locations have a significant impact on the quality and flexibility of approximation, particularly in the approximation of functions with discontinuities (López de Silanes et al. 2001;Sinha and Schunck 1992) and functions with locally high gradients (Malagù et al. 2014). Furthermore, through the proper choosing of the knot sequence, splines can accurately approximate very complex functions, such as the shapes of three-dimensional objects in a computer-aided geometric design (Hoschek and Lasser 1993;Park and Lee 2007). Therefore, splines are chosen herein to approximate complex multidimensional joint distributions. The most commonly used mathematical formulation in different applications of splines are B-splines (from basis spline). The construction of B-splines is straightforward and simple to implement; however, the high-order simulation framework proposed by Mustapha and Dimitrakopoulos (2010) assumes the orthogonality of basis functions, and, therefore, splines in the form of B-splines are not suitable for a high-order spatial simulation approach.
In this paper, Legendre-like splines (Wei et al. 2013) are used, which are shown to be orthogonal and can be easily integrated in the high-order simulation framework. There are two user-defined parameters used for the constructing of Legendre-like splines: the order of splines and the maximum number of knots. In practice, cubic splines (order 3) are commonly used (Hughes et al. 2005;Piegl 1989), as they provide efficient smooth approximation. For the cubic splines, the first four Legendre-like splines are defined at two endpoint knots of the approximation domain and Legendre polynomials up to order 3. Next, Legendre-like splines are constructed by adding an additional knot per Legendre-like spline until the user-defined maximum number of knots is reached. Increasing the number of knots improves the approximation and describes more complex relations in the available data in the same way that the high-order polynomials capture the complex behavior of the function to be approximated. Thus, the maximum number of knots reflects the maximum order of high-order spatial relations that can be calculated from the available data. This spline approach aims to improve the estimation of the conditional probability density function (cpdf) and overcome the limitations of polynomial approximations. In addition, the proposed approach provides a general framework for high-order simulation techniques. For example, by using only one interval for spline construction, the technique becomes the one proposed by Dimitrakopoulos (2010, 2011).
The paper is organized as follows. First, the high-order simulation framework is outlined. Then, two systems of basis functions are outlined: Legendre polynomials (Lebedev 1965) and Legendre-like orthogonal splines. In the following section, the capabilities of both systems are compared using a fully known dataset to demonstrate the advantages of orthogonal splines in simulating connected high values. Next, the proposed approach is applied to a gold deposit and compared with the sequential Gaussian simulation approach in terms of the reproduction of histograms, variograms, high-order spatial statistics, and the connectivity of high values. Discussion and conclusions follow. Supplementary material available online provides the C++ course code of the high-order sequential simulation implementation detailed in Sect. 2.

Sequential High-Order Simulation
Let Z(u i ) be a stationary ergodic random field indexed in R n , where u i ∈ D ⊆ R n (n 1, 2, 3), i 1 . . . N and where N is the number of points in a discrete grid D ⊆ R n . Random variables indexed on the grid D ⊆ R n are denoted by Z i ≡ Z (u i ), whereas their outcomes are denoted by z i z(u i ). The focus of high-order simulation techniques is to simulate the realization of the random field Z (u i ) for all nodes of a grid D with a given set of conditioning data d n {z(u α ), α 1 . . . n}.
The joint probability density function f (u 0 , u 1 , . . . u N ; z 0 , z 1 , . . . z N |d n ) of the random field Z (u i ) can be decomposed into the product of conditional univariate distributions using the basic concept of sequential simulation (Journel and Alabert 1989;Journel 1994;Dimitrakopoulos and Luo 2004) Accordingly, the random path of visiting all grid nodes is defined first. Then, starting from the first node in the random path, the value z i is simulated based on the estimated cpdf f (u i ; z i |z 1 , . . . , z i−1 , d n ). Finally, the simulated value is added to the set of conditional data, and the process is repeated until all grid nodes in the random path are visited. Eventually, any resulting simulation represents a realization of the complex joint distribution f (u 0 , u 1 , . . . u N ; z 0 , z 1 , . . . z N |d n ). Without loss of generality, let u 0 be the first node in the random path. According to Bayes' rule (Lee 2012) where f (u 0 , u 1 , . . . , u n ; z 0 , d n ) is a joint probability density function and f (u 1 , . . . , u n ; d n ) can be calculated as f (u 1 , . . . , u n ; d n ) f (u 0 , u 1 , . . . , u n ; ξ 0 , d n )dξ 0 . (3) In this paper, the joint probability density function f (u 0 , u 1 , . . . , u n ; z 0 , d n ) is approximated using Legendre polynomials and Legendre-like orthogonal splines.

Approximation of a Joint Probability Density Using Orthogonal Functions
Let f (z) be a probability density function of a random variable Z defined on [a, b] and let ϕ 1 (z), ϕ 2 (z), . . . be a complete system of orthogonal functions in [a, b], then f (z) can be approximated by the finite number ω of functions ϕ 1 (z), ϕ 2 (z), . . . ϕ ω (z) where L m are coefficients of approximation. The system of functions where δ mk 1, m k 0, m k is the Kronecker delta, and, therefore, ∀k 0 .
By definition where E stands for mathematical expectation. The coefficients L m can be estimated from the available data. Similarly, a joint probability density function can be approximated as where |dz| dz 0 dz 1 · · · dz n . By definition When considering the spatial locations u {u 0 , u 1 , . . . , u n } of random variables Z 0 , Z 1 , . . . Z n , the coefficients L k 0 ,k 1 ,...,k n can be estimated from available data using Eqs. (9) and (10) by calculating (11) where values z k i , i 0 . . . n are taken from the available data z k i ∈ d n and the given training image, and separated by lags Finally, high-order sequential simulations are generated using the following algorithm: Algorithm A.1 1. Define a random path for visiting all unsampled nodes on the simulation grid. 2. For each node u 0 in the path: a. Find the closest sampled grid nodes u 1 , Calculate the coefficients L k 0 ,k 1 ,...,k n using Eq. (11). e. Build the cpdf f (u 0 ; z 0 |z 1 , . . . z n ) for the random variable Z 0 at the unsampled location u 0 given the conditioning data z 1 , . . . z n at the corresponding neighbors u 1 , u 2 , . . . u n using Eqs. (2) and (8). f. Draw a uniform random value in [0, 1] to generate a simulated value z 0 from the conditional distribution f (u 0 ; z 0 |z 1 , . . . z n ).
g. Add z 0 to the set of sample data and the previously simulated values. 3. Repeat Steps 2a-g for the next points in the random path defined in Step 1. Mustapha and Dimitrakopoulos (2010) proposed using a Legendre series as a set of basis functions ϕ 1 (z), ϕ 2 (z), . . .. The Legendre polynomial P k of order k is defined as in Lebedev (1965)

Legendre Polynomials
The set of Legendre polynomials {P k (z)} k forms a complete basis set on the interval [− 1, 1], and, accordingly, the function f (u 0 ; z 0 |z 1 , . . . z n ) can be approximated using Eqs. (2) and (8). By their construction, the order of Legendre polynomials corresponds to the order of high-order spatial statistics of the probability function f (u 0 ; z 0 |z 1 , . . . z n ). However, there are practical limitations when using Legendre polynomials for the approximation of functions in multidimensional space. This is discussed in Sect. 3.

Legendre-Like Orthogonal Splines Approximation
In the present work, Legendre-like splines (Wei et al. 2013) are used as a set of basis functions. These splines are constructed using Legendre polynomials and the linear combination of B-splines. B-splines of order r in a variable t ∈ [a, b] are piecewise polynomials defined over the domain The points t i , i 0 . . . m max are called knots. Each piece, a B-spline of order r, can be derived using de Boor's formula (de Boor 1978) B-splines do not form an orthogonal basis; however, Wei et al. (2013) introduced orthogonal splines based on the combination of B-splines and a set of knot sequences. The first r + 1 splines are defined as the Legendre polynomials up to order r The subsequent splines are constructed on subsets T m For example, the first and second subsets are then, the remaining Legendre-like splines S k (t), k r +2 . . . r +m max are determined by where f m (t) is the determinant of the matrix: The examples of orthogonal splines of order r 3 and the knot sequence T [−1, −1, −1, −1, −0.6, −0.2, 0.2, 0.6, 1, 1, 1, 1] are presented in Figs. 1 and 2. The first r + 1 splines are defined on the knot sequence with only one interval [− 1, 1] and are, thus, simply Legendre polynomials up to the order r (Fig. 1). For each subsequent spline, the knot sequence is updated by adding a knot from the initial knot sequence T , e.g., the fifth spline ( In this work, the initial values are linearly transformed into the [− 1, 1] values range and divided into m max intervals. There are two parameters that have a significant impact on the quality of approximation: the maximum number of intervals m max Fig. 1 The first four Legendre-like splines over the knot sequence T Fig. 2 The last four Legendre-like splines over the knot sequence T and the order of splines r. These parameters reflect the maximum order of high-order spatial statistics that can be captured from the available data. The first parameter is the order of splines r. High values of the order r lead to a non-stable approximation (Runge 1901;Platte et al. 2011) and high computational costs, whereas low values affect the continuity and smoothness of approximation. For example, zero-order splines are good for the approximation of the stepwise function because each spline is a constant function. Splines with r 1 are used for the approximation of continuous, but not smooth, functions, i.e., a polygonal line. In practice, cubic splines, i.e., r 3, are commonly used (Hughes et al. 2005;Piegl 1989). The second parameter is the maximum number of intervals m max . Low values of m max lead to an approximation that is close to a polynomial case, for which limitations are discussed in Sect. 4. For example, an approximation with m max 1 corresponds to the Legendre polynomial approximation presented by Mustapha and Dimitrakopoulos (2010). High values of m max result in overfitting or poor predictive performance, as it overreacts to minor fluctuations in the data. In addition, an approximation with a high value of m max affects the variability of the simulations because it directly samples values from the initial available data and pastes them into simulations. To choose values of r and m max , different measures are tested. The widely known Kolmogorov-Smirnov statistics test (Stephens 1974) that indicates whether two data samples come from the same distribution is not utilized herein because the related quantile-quantile plot is reproduced well for a very wide range of r and m max values; thus, such a statistical test does not provide guidance on selecting suitable parameters. Other approaches, including comparing high-order spatial statistics maps or connectivity properties, are hard to quantify by a single number and complex to implement. In this work, a simple and fast measure of the quality of approximation is used. This quality of approximation is expressed in terms of the number of grid nodes where splines fail to approximate the conditional distribution. At these nodes, the high-order simulation method produces numerical artefacts, such as outliers or noise values. Outliers can be easily detected by comparing the value at nodes with their local neighborhood average value. The average number of outliers is calculated for different values of r and m max (Fig. 3). According to Fig. 3, as the number of intervals m max is increased, the quality of approximation improves. At the same time, increasing the order of splines r decreases the quality of approximation. For cubic splines (order r 3), the reasonable number of intervals is 30, as it provides the same quality of approximation as 50, demonstrates better predictive performance, and is computationally less expensive. The corresponding order of high-order spatial statistics is m max + r 33.

Testing the Simulation Approach with a Fully Known Dataset
The high-order simulation method presented in the previous section is tested with a fully known dataset obtained from an image of a fracture network downloaded from a texture synthesis website (http://br.depositphotos.com/5211338/stock-photodry-terrain.html) and displayed in Fig. 4. Gray-scale values of the image are transformed to the [0, 1] domain. The reference image (Fig. 5a) and the TI (Fig. 5b) are taken from different parts of the image and have sizes 150 × 150 and 200 × 400 grid nodes, respectively. The dataset (Fig. 5c) is generated from the reference image with a uniform random sampling of 225 values (1% of the image points). Three different systems of functions for the high-order simulation approach (hosim) and the sequential Gaussian simulation (sgsim) method are compared next: (a) Legendre-like splines (r 3, m max 30, the corresponding order of high-order spatial statistics is m max + r 33), (b) sgsim method, (c) Legendre polynomials of order 10 (the corresponding order of high-order spatial statistics is 10), and (d) Legendre polynomials of order 20 (the corresponding order of high-order spatial statistics is 20). Hereafter, the system of functions based on splines and Legendre polynomials are correspondingly called hosim-splines and hosim-polynomials. The simulation using hosim-splines (Fig. 6a) shows a stable reproduction of spatially connected structures. Simulations using hosim-polynomials (Fig. 6c, d) have less connected features than the simulation using splines. sgsim (Fig. 6b) fails to reproduce the spatial continuity of high values. Table 1 shows the average value, median, and variance for the sample data, reference image, TI, and simulated realizations; note that only hosim-polynomials of order 20 are included in the comparisons that follow. All methods reproduce well the low-order statistics of the sample data and the TI. Figure 7 shows the quantile-quantile (QQ) plots of ten simulated realizations of each simulation approach with the sample data. The QQ plots for the simulations using hosim-splines are represented by red lines. The 45°black line represent QQ plots of the sample data with the sample data. The blue line represents the QQ plot of the reference image with the sample data. The green line represents the QQ plot of the training image with the sample data. The QQ plots for the simulations using sgsim are depicted by gray dashed lines and the QQ plots for the simulations using hosimpolynomials are shown by gray solid lines. Overall, the QQ plots of simulations using hosim-splines are consistent with the QQ plot of the sample data and the reference image, whereas QQ plots of simulations using hosim-polynomials and sgsim slightly deviate from the QQ plot of the sample data. Figure 8 shows variograms along the north-east and north-west directions; that is, directions of the main continuity of high values, calculated for the simulations using hosim-splines (the red lines), the sample data (dots), the reference image (the blue line), the TI (the green line), the simulations using hosim-polynomials (the gray solid lines), and the simulations using sgsim (the gray dashed lines). All techniques demonstrate reasonable reproduction of the second-order statistics of the sample data.
For the calculation of the third-order and fourth-order spatial statistics, the estimations of the high-order moment are used ) Fig. 6 The simulation results are as follows: a the simulation using hosim-splines, b the simulation using sgsim, and c, d the simulations using hosim-polynomials of orders 10 and 20, respectively  where N h1h2 is the number of elements of replicates found on lags h 1 and h 2 , and N h1h2h3 is the number of elements of replicates found on distances h 1 , h 2 , and h 3 . To highlight the connectivity property along the north-east (NE) and north-west (NW) directions, the third-order moments are calculated for binary images with a cut-off value of 0.82 (95th percentile). An example of a binary image is shown in Fig. 9a. The third-order spatial statistics are estimated based on a template with directional vectors along the NE and NW directions (Fig. 9b), i.e., h 1 (idx, idy) and h 2 (− jdx, jdy), respectively, where i, j 1 . . . 30 and the lag discretization along x and y is dx dy 1 pixel. The physical meaning of the third-order moment of the binary image is straightforward-it is the probability of having high values at the three points separated by lags h 1 and h 2 (Minniakhmetov and Dimitrakopoulos 2017b). The red-orange values represent the average sizes of connected high values along the NE and NW directions. In the third-order indicator moment map of the reference image, the average sizes of the interconnected high values are 10 and 20 pixels along the NE and NW directions, respectively.
The third-order moments are calculated for simulations and averaged to account for differences between the realizations. The high-order simulation technique using hosim-splines and hosim-polynomials (Fig. 10b, e) reproduce the third-order moment map of the sample data (Fig. 10a), the reference image (Fig. 10c), and the TI (Fig. 10d), as can be seen from the similar size of the red-orange value areas in the corresponding figures. The moment map of the simulation using sgsim (Fig. 10f) does not reproduce connectivity along the NE and NW directions; the size of the red-shaded area is 8 × 10 pixels, compared to a 10 × 20-pixel area in the reference image's moment map (Fig. 10c).
The fourth-order spatial statistics are estimated for binary images with a cut-off value of 0.82 (95th percentile) based on a template with directional vectors NE h 1 (idx, idy), NW h 2 (− jdx, jdy), and south-west (SW) h 3 (−kdx, −kdy), where i, j, k 1 . . . 30 and lag discretization along x and y is dx dy 1 pixels. Similarly to  top subfigures a, b) and north-west direction (bottom subfigures c, d) of the sample data (dots), reference image (blue line), training image (green line), simulations using hosim-splines (red lines), simulations using hosim-polynomials (gray lines), and simulations using sgsim (gray dashed lines) Fig. 9 The third-order indicator moment calculation for the reference image: a the binary values of the reference image for a cut-off value of 0.82 (95th percentile), the black lines represent L-type template for calculation of the third-order spatial statistics; b the third-order spatial moments of the binary image Fig. 10 The third-order indicator moments: a sample data samples, b the simulation using hosim-splines, c the reference image, d the training image, e a simulation using hosim-polynomials, and f a simulation using sgsim the third order, the fourth-order moments are calculated for simulations and averaged to account for differences in the various realizations. The red-orange areas along the axes of the fourth-order spatial statistics (Fig. 11) represent the high values along the NE, NW, and SW directions. According to Fig. 11, the fourth-order moment map for the simulation using hosim-splines (Fig. 11b) reproduce the sizes of fractures along the NE, NW, and SW directions in the fourth-order moment of the sample data (Fig. 11a), the reference image (Fig. 11c), and the TI (Fig. 11d). The fourth-order moment map of the simulation using hosim-polynomials (Fig. 11e) shows a smaller connectivity of fractures along the NE and SW directions. The spatial statistics map of the simulation Fig. 11 The fourth-order indicator moments: a sample data samples, b the simulation using hosim-splines, c the reference image, d the training image, e a simulation using hosim-polynomials, and f a simulation using sgsim using sgsim (Fig. 11f) does not reproduce the connectivity of fractures along the NE and SW directions.
The connectivity of high values is measured using the function presented by Journel and Alabert (1989). As in the above examples, the cut-off value is equal to 0.82 (95th percentile). Figure 12 shows P50 statistics of the connectivity measure along the NE (top subfigures) and NW directions (bottom subfigures). The P50 statistics of connectivity are calculated for the simulations using the proposed techniques (red solid line), hosim-polynomials (gray solid line), and sgsim method (gray dashed line). The connectivity measures of the reference image (blue line) and the TI (green line) falls within the P10 and P90 statistics of the connectivity measure in the simulations Fig. 12 The connectivity functions along the north-east (top subfigures) and north-west directions (bottom subfigures) 95th percentile: reference image (blue line), training image (green line), P50 statistics for simulations using hosim-splines (red solid line), P10 and P90 statistics for simulations using hosim-splines (red dash-dot line), P50 statistics for simulations using hosim-polynomials (gray solid line), and P50 statistics for simulations using sgsim (gray dashed line) using hosim-splines (red dash-dot lines), whereas the connectivity of the simulations using hosim-polynomials and sgsim is lower, on average, than the connectivity of the TI and the reference image.

Application at a Gold Deposit
Data from a gold deposit are used as a case study to demonstrate the intricacies and advantages of the high-order spatial simulation method described above. In addition, the method is compared with the sgsim approach for the reproduction of histograms, variograms, high-order spatial statistics, and the connectivity of high and extreme values.
The deposit is 2 km by 2 km wide and extends to a depth of 500 m. Sample data are available from 288 exploration drillholes. Blast-hole data are also available for the deposit and used in the construction of a training image. The three-dimensional TI is defined on 405 × 445 × 86 grid blocks of size 5 × 5×5 m 3 . The simulation grid The two-dimensional sections in Fig. 14 show that the simulation using the proposed method reproduces the spatial distribution of grades and the continuity of high grades. The areas with high values in Fig. 14 are in good agreement with the drillhole data ( Fig. 15) and the TI (Fig. 13). The simulation using sgsim (Fig. 16) exhibits a greater number of disconnected structures with high values and sparsely distributed outliers.
These observations are confirmed by a quantitative analysis in a subsequent validation by (1) mean and variance comparison, (2) QQ plots between drillhole data and simulated values, (3) variogram validation, (4) high-order spatial cumulant validation, and (5) connectivity measure. Table 2 shows the average value, median, and variance for the drillhole data, the TI, and the simulations. Both methods reproduce well the low-order statistics of the drillhole data and the TI.  Figure 17 shows the QQ plots of the simulated realizations and the drillhole data. Quantiles of simulations using hosim-splines are shown by red lines. Quantiles of simulations using the sgsim method are shown by gray lines. In addition, the QQ plots of the training image and the drillhole data are depicted by the blue line. The closer these curves are to the 45°black line in the graph, the better they reproduce the distribution of the drillhole data. Both methods provide simulations consistent with the drillhole data in terms of distributions. Figure 18 presents variograms for the northsouth and east-west directions. Simulations using hosim-splines (red lines) share the second-order statistics of drillhole data and the TI (blue lines). The simulations using the sgsim method (gray lines) preserve the second-order statistics of the drillhole data (black dots).
Applying a zero-mean transformation, the third-order cumulants can be calculated using Eq. (22) with lags h 1 (idx, 0) h 2 (0, jdy) indexed by i 1 . . . 7, j 1 . . . 7, where dx and dy are distances between drillholes, that is, 100 m × 100 m. Figure 19 shows the comparison of cumulant maps for sample data, the TI, and the simulations. The values along axes reflect variograms along their corresponding directions because the third-order moment E(Z 2 (x)Z (x + h)) has similar spatial relations as the second-order moment E(Z (x)Z (x + h)). However, the square term Z 2 (x) in  (Fig. 19c) reproduce red areas along the x-y axes and yellow-green areas in the cumulant map of the drillhole data (Fig. 19a) and the TI (Fig. 19b). These areas reflect the size of connected high grades and are equal to approximately 400 m along the x-axis and 300 m along the y-axis. The cumulant map for the simulation using sgsim (Fig. 19d)  The fourth-order cumulants are calculated using the following equation   where N h1h2h3 is the number of elements of replicates found on distances h 1 and h 2 , and m 2 (h) is the second-order moment along direction h, which is equal to the covariance for a zero-mean random field. The lags h 1 (idx, 0), h 2 (0, jdy), and h 3 (0, kdz) are indexed by i 1…7, j 1…7, and k 1…7, where, dx, dy, and dz are distances between data samples, that is, 100 m × 100 m × 5 m. The high-order Fig. 20 The fourth-order spatial cumulant maps of the a drillhole data, b training image (TI), c a simulation using hosim-splines, and d a simulation using sgsim Fig. 21 The connectivity along x (left subfigure) and y (right subfigure) for the 99th percentile: the training image (TI, blue lines), P50 statistics for hosim simulations (red solid lines), P10 and P90 statistics for hosim simulations (red dashed lines), and P50 statistics for simulations using sgsim (gray lines) cumulants calculated reflect the complex structures of orebodies ). According to Figs. 19 and 20, the size of connected structures is reproduced in simulations using the proposed method (Figs. 19c,20c). This can also be traced in the vertical profiles (Figs. 13c, 14c, 15c). The fourth-order cumulant map of the simulation using sgsim (Fig. 20d) has a rather small red area in comparison with structures in the cumulant maps of the drillhole data (Fig. 20a) and the TI (Fig. 20b).
The connectivity along the x and y axes is analyzed using the connectivity measure presented by Journel and Alabert (1989). The cut-off value is equal to 5 ppm (99th percentile). The P10, P50, and P90 statistics of connectivity measures are calculated for simulations using the hosim-splines method and depicted by red lines in Fig. 21. Solid lines represent the P50 of connectivity measures, whereas dashed lines show the P10 and P90 statistics. The connectivity of the simulations using the proposed method (red lines) remains close to the connectivity measure of the TI (blue lines). The P50 statistics of connectivity measure calculated using sgsim simulations (gray lines) is quite far from the connectivity of the TI. Thus, despite reproducing the histograms and variograms, Gaussian simulation methods fail to reproduce an important property of the connectivity of high values.

Conclusions
This paper presents a novel approach for the high-order simulation of continuous variables based on Legendre-like orthogonal splines. Splines are flexible tools for the approximation of complex probability density functions. Using different knot sequences, orders of splines, and smoothness of piecewise polynomials, it is possible to obtain a stable approximation that reproduces the spatial connectivity of the extreme values. The simulations are consistent with the spatial statistics of the sample data and share the high-order spatial statistics of the available data and the training image.
The proposed approach is also compared with the conventional second-order approach sequential Gaussian simulation and the high-order simulation method using Legendre polynomials. The approach using splines exhibits a more stable approximation of the conditional probability density function (cpdf) and a better representation of the spatial connectivity of extreme values. The applied connectivity measure confirms the results obtained by analyzing the high-order statistics and demonstrates the limitations of Gaussian simulation methods in the characterization of a mineral deposit. In addition, the proposed approach provides a general framework for high-order simulation techniques. For example, by using just one interval for spline construction, the technique reproduces the method proposed by Dimitrakopoulos (2010, 2011).
Further research will address the simulation of categorical variables using splines of order zero and the simulation of multiple correlated continuous and discrete variables within a general framework. In addition, the adaptive knot sequence will be investigated for a better approximation of the cpdf.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.