Multifractal interpolation method for spatial data with singularities

Since the term ‘kriging’ was coined by Georges Matheron in the early 1960s on the basis of Krige’s master’s thesis dealing with interpolation of point samples, geostatistics has been rapidly developed as a branch of science and relevant techniques have been commonly applied in many fields of science for mapping, estimation, simulation, and prediction (Journel and Huijbregts, 1978; Goovaerts, 1997). The International Association for Mathematical Geosciences (IAMG) is proud of the invention and further development of the subject by our members. Kriging and other geostatistical techniques have been widely applied outside of geosciences, where users unaware of its origins and mathematical evolution refer to it simply as a type of spatial analysis. The semivariogram, a function of distance between locations, can measure the spatial autocorrelation between values at locations separated by a distance. Models empirically fitted to semivariograms are used for assigning weights to linear equations whose solutions provide weighted averages for kriging data with stationarity (Goovaerts, 1997; Deutsch and Journel, 2008). Traditionally, kriging is for interpolating data in the neighbourhood and estimating values at locations where no data is available. Interpolation algorithms have been developed for a variety of simple, indicator, and higher-order kriging as well as kriging with transformed and compositional data. Algorithms for interpolation of data with anisotropic spatial association (e.g. Chiles and Delfiner, 1999), mixed categorical and/or continuous data (Journel and Huijbregts, 1978; Goovaerts, 1997), and compositional data (Pawlowsky-Glahn and Olea, 2004), have been created. Case studies comparing these methods are available in the literature (e.g. Park and Jang, 2014). Application of kriging depends heavily on stationarity of the mean and second-order moments involving the variogram and standard deviation of a regionalized random variable. Simple kriging (SK) may be applied if the mean of the data has a known but constant trend, whereas ordinary kriging (OK) may be applied if the trend is constant but unknown. If the trend is unknown but follows some polynomial model, other types of kriging accounting for trends can be used (Hansen et al., 2010). In most cases stationarity of second order moments is also assumed. However, the real data, especially exploratory data involved in characterizing mineralization and hazardous events, often does not meet stationarity requirements because of singularities. Multifractal interpolation method for spatial data with singularities


Introduction
Since the term 'kriging' was coined by Georges Matheron in the early 1960s on the basis of Krige's master's thesis dealing with interpolation of point samples, geostatistics has been rapidly developed as a branch of science and relevant techniques have been commonly applied in many fields of science for mapping, estimation, simulation, and prediction (Journel and Huijbregts, 1978;Goovaerts, 1997). The International Association for Mathematical Geosciences (IAMG) is proud of the invention and further development of the subject by our members. Kriging and other geostatistical techniques have been widely applied outside of geosciences, where users unaware of its origins and mathematical evolution refer to it simply as a type of spatial analysis. The semivariogram, a function of distance between locations, can measure the spatial autocorrelation between values at locations separated by a distance. Models empirically fitted to semivariograms are used for assigning weights to linear equations whose solutions provide weighted averages for kriging data with stationarity (Goovaerts, 1997;Deutsch and Journel, 2008). Traditionally, kriging is for interpolating data in the neighbourhood and estimating values at locations where no data is available. Interpolation algorithms have been developed for a variety of simple, indicator, and higher-order kriging as well as kriging with transformed and compositional data. Algorithms for interpolation of data with anisotropic spatial association (e.g. Chiles and Delfiner, 1999), mixed categorical and/or continuous data (Journel and Huijbregts, 1978;Goovaerts, 1997), and compositional data (Pawlowsky-Glahn and Olea, 2004), have been created.
Case studies comparing these methods are available in the literature (e.g. Park and Jang, 2014). Application of kriging depends heavily on stationarity of the mean and second-order moments involving the variogram and standard deviation of a regionalized random variable. Simple kriging (SK) may be applied if the mean of the data has a known but constant trend, whereas ordinary kriging (OK) may be applied if the trend is constant but unknown. If the trend is unknown but follows some polynomial model, other types of kriging accounting for trends can be used (Hansen et al., 2010). In most cases stationarity of second order moments is also assumed. However, the real data, especially exploratory data involved in characterizing mineralization and hazardous events, often does not meet stationarity requirements because of singularities.
A new approach is the multiple point geostatistics which is a new field of spatial and temporal analysis (Mariethoz and Caers, 2014). Multiple point geostatistics uses training images to quantify and to include structural information into stochastic simulation (Guardiano and Srivastava, 1993;Strebelle, 2002).
The multifractal interpolation method (MIM) based on multifractal theory (Cheng, 1999a(Cheng, , 2000 has been developed for the analysis and interpolation of data with singularities. Multifractal theory integrates both spatial association and local singularities and can enhance and retain the local structure properties (Cheng, 2006a,b). This paper introduces a generalized binomial multiplicative cascade process to demonstrate the generation of one-and two-dimensional data with multi-scale singularities which can be modelled by asymmetrical multifractal distribution. It will then introduce a generalized moving average mathematical model for analyzing and interpolating between data with singularities.
Finally the method will be demonstrated by a onedimensional example.

Multiplicative cascade processes and multifractal distributions
Singular physical, chemical, and biological processes can result in anomalous energy release, mass accumulation, or matter concentration, all of which are generally confined to narrow intervals in space or time (Cheng, 2007a). The end products of these nonlinear processes can be modelled as fractals or multifractals. Singularity is a property of nonlinear natural processes, examples of which include, but are not limited to, cloud formation (Schertzer and Lovejoy 1987), rainfall (Veneziano 2002), hurricanes (Sornette, 2004), flooding (Malamud et al., 1996;Cheng 2008;, landslides (Malamud et al., 2004), forest fires (Malamud et al.,1996), earthquakes (Turcotte 1997), mineral deposits and related geochemical anomalies (Agterberg, 1995;Cheng et al., 1994;Xie and Bao, 2004;Cheng and Agterberg, 2009), solar wind turbulence (Macek, 2007), DNA series (Rosas et al., 2002), heat flow at the mid oceanic ridges (to be published by the author elsewhere) and landscape textures (Plotnick et al., 1993). Multifractal modeling involves quantification of multi-scale singularities and various types of properties associated with spatial distribution of the singularities (Halsey et al., 1986;Feder, 1988;Mandelbrot, 1989;Evertsz and Mandelbrot, 1992). This section introduces an asymmetrical cascade process that generates results with singularities characterized by asymmetrical multifractal models. There are several formalisms for describing multifractal distributions, one of which is the multifractal model based on the partition function (Halsey et al., 1986). This model involves three functions: the mass exponent function or Renyi exponent τ(q), the coarse Hölder exponent α(q), and the fractal spectrum function f(α) (Halsey et al., 1986).
In order to demonstrate the singularities and their distributions in one-or two-dimendional data, I introduce the theories and concepts of multiplicative cascade processes (MCPs), which play a fundamental role in quantifying turbulent intermittency and other nonlinear processes Lovejoy 1985, Schertzer et al. 1997). MCPs have been extensively discussed in the literature (e.g. Gupta and Waymire 1993, Over and Gupta 1996, Menabde and Sivapalan 2000Serinaldi 2010). The model of de Wijs is a simple two-dimensional multiplicative cascade model (de Wijs 1951, Agterberg 2001 described in terms of multiplicative canonical cascade processes (Lovejoy and Schertzer, 2007). Other modifications exist, e.g. a cascade model with functional redistribution rate (Agterberg 2007b); a two-dimensional cascade model with anisotropic partition (Cheng 2005); a generalized two-parameter binomial multiplicative model as proposed by Koscielny-Bunde et al. (2006) and applied for describing multifractal spectra of runoff time series; a three-parameter binomial multifractal model proposed by Macek (2007) and applied to characterize solar wind turbulence data based on a generalized two-scale weighted Cantor set for characterizing asymmetrical multifractal distribution; a two-dimensional cascade model with variable and conditional dependence partition (Cheng, 2012); and a five-parameter binomial multiplicative cascade model has been recently proposed by the author (Cheng, 2014) for describing fundamental multifractal indices characterizing asymmetrical multifractal distribution of real-world data.
The five-parameter generalized multiplicative cascade processes involve the partitioning of with a length L (e.g. in unit of meter) into h equal segments of which m 1 receive d 1 (> 0) and m 2 receive d 2 (> 0) proportion of mass in the previous segment, respectively, where m 1 + m 2 ≤ h. For a closed system with preservation of unit measure (e.g. with total mass M), d 1 + d 2 = 1. Otherwise, d 1 + d 2 < 1 or d 1 + d 2 > 1, corresponding to a loss or gain of mass during the cascade process, respectively. At the nth partition, the segment length will be ε n = L(1/h) n ; the segments are subject to k times segments with measure d 1 /m 1 and n -k times segments with measure d 2 /m 2 , and thus the measures of these segments are ′ κ = M(d 1 /m 1 ) k (d 2 /m 2 ) n-k . Therefore, the numbers of segments with ′ k will be N k = m 1 k m 2 n-k ( k n ). Letting n →∞, we can find series of n and k with k = ξn, 0 ≤ ξ ≤ 1, where ξ is independent of n or k. We then obtain the following relationships where α is the coarse Hölder exponent which quantifies the singularity of the distribution of ′ and the subset of segments with singularity α is an intertwined set which is a fractal with fractal dimension f(α).
The number of segments in each of the intertwined fractals can be expressed as [3] where f(α) is the fractal dimension spectral function, which can be expressed as [4] It can be seen from Equations [1] to [4] that both the measure ′ κ and N k follow power-law relationships at scale ε n . Since the value of ξ falls within the range [0, 1], the value of the singularity index takes any value between α min and α max following the linear relationship between α(ξ) and ξ: where the maximum and minimum values of singularity α are [6] assuming d 1 /m 1 > d 2 /m 2 ; otherwise, the two extremes will be reversed. Accordingly, the corresponding fractal dimensions with special singularities are shown to be The multifractality and symmetry of the multifractal distribution can be characterized by the asymmetry and multifractality indices [8] The asymmetry index corresponds to the ratio of the values m 2 and m 1 ; whereas the multifractality is proportional to the ratio of average measures, (d 1 /m 1 )/(d 2 /m 2 ).

Singularities and nonstationarity
The singularity in the multifractal model introduced in the previous section characterizes how the statistical behaviour varies as the measuring scale changes. For example, at some locations the mean value calculated from the neighbourhood values might be independent of the size of the vicinity within which the values are averaged. In other cases, the mean value might proportionally depend on the size of the vicinity. The former case represents nonsingular location but the latter is for singular location. Singularity property has been commonly observed in geochemical and geophysical quantities (Cheng et al., 1994). Generally speaking, as the size of segment ε n → 0 (n →∞), then the measure defined for each segment (Equation [1]) tends to zero and the number of segments tends to infinity. In order to explain the singularity of geochemical and geophysical quantities according to the notation of the multifractal model shown in Equations [1]-[7], the 'fractal density' of measure with singularity (α) is defined by the author as the ratio of mass ′ (ε n ) over scale ε n α which can be calculated as follows: the new fractal density ρ α has a unit of M over the unit of L α , for example, g/m 0.3 . [10] where ′ α = 1-α represents an index quantifying the local singularity of the measure at locations with singularity α. The ordinary density ρ can be decomposed into two components: fractal density ρ α and ε -′ α , the former is independent of scale ε whereas the latter dependent on the scale. The former component is a non-singular component and the latter is singular component if the singularity index ′ α ≠ 0. The inconvenience property of the measure with following singularity properties implies nonstationarity of the measure or the density (Cheng, 1999a): (1) α = 1 or ′ α = 0, if ρ(ε) = constant, independent of scale size ε (2) α < 1, ′ α > 0, if ρ(ε) is a decreasing function of ε, which normally implies the 'convex' property of ′ (ε) at the location with α (3) α > 1, ′ α < 0, if ρ(ε) is an increasing function of ε, which normally implies the 'concave' property of ′ (ε) at the location with α.
Therefore, the α-values as the fractal dimension of the fractal density (Δα -value as the co-dimension) can be used to characterize the nonlinear structural property of the measure ′ (ε). This approach has already been used for texture analysis of remote sensing Landsat TM images (Cheng, 1997b(Cheng, , 1999c, in multifractal interpolation of geochemical concentration values for mineral exploration (Cheng, 1999a(Cheng, , 2000a(Cheng, , 2000b, in well log curve reconstruction (Li and Cheng, 2001), flood event modelling , in hyperspectral image analysis (Neta et al., 2010), faults and geochemistry (Wang et al., 2013), and in geochemical element concentration mapping (Chen et al., 2007;Zuo et al., 2009).
In order to introduce how singularities can be included in data interpolation, we here introduce the following scaleinvariant property of the measure, ′ (ε) and density, ρ(ε). Due to the properties of power-law type functions we can associate the densities at two different scales (ε n and ε m ) as follows [11] The MIM to be introduced in the next section is developed according to the scale invariance property (Equation [11]).

MIM incorporating spatial association and singularity
Statistical properties derived at one scale may be used to estimate properties at another scale on the basis of the scaling property (Cheng, 1999a(Cheng, , 2000. The main purpose of data interpolation, including kriging, is to estimate values at unknown locations based on the neighbourhood values and their spatial association with the value being estimated. Spatial association represents a type of statistical dependency of values at separate locations. If the value at a location (x) is considered as the realization of a regionalized random variable Z(x), the spatial association or variability can be measured by means of the variogram [12] where the semivariogram γ(x, h) is a function of vector distance h separating locations x and x + h, measuring the symmetrical variability between Z(x) and Z(x + h). Under an assumption of second-order stationary, the semivariogram

(Equation [12]) becomes a function of h that is independent
of location x. This strong assumption on the nature of the regionalized random variable is generally required in kriging. Equation [12] has been commonly used for structural

Multifractal interpolation method for spatial data with singularities 237
The Journal of The Southern African Institute of Mining and Metallurgy VOLUME 115 MARCH 2015 L Multifractal interpolation method for spatial data with singularities analysis and interpolation in geostatistics (Journel and Huijbregts, 1978). It has also been applied for texture analysis in image processing (Atkinson and Lewis, 2000;Herzfeld and Higginson, 1996).
To incorporate both spatial association and singularity in supporting the interpolation model based on Equation [11], the following average density within a small vicinity ⏐ (x 0 , ε) around location x 0 with linear size ε was defined by Cheng (2006) [13] Assume that Equation [11] holds true within a range of window sizes, ε ≤ ε ≤ ε max . Then the average density ρ(ε, x 0 ) within the window ⏐ (x 0 , ε), where it may not contain samples with observed values, can be associated with the average density within the larger window ⏐ (x 0 , ε max ) where it contains samples with observed values and can be estimated by kriging as follows: [14] Equation [14] is a general weighted average model that can be used to estimate the value at the centre of ρ(ε, x 0 ) from the neighbourhood values within ⏐ (x 0 , ε max ) (Cheng, 1999a(Cheng, , b, 2006. Since the above discussions are valid for all dimensions and here we will use E to present the dimension of problem, E=1, 2, 3 stand for 1D, 2D and 3D problems. It has the following properties : (1) If it does not show singularity, α = E or ′ α = 0, then Equation [14] reduces to the ordinary moving average function that has been used commonly in kriging and other data interpolation methods (2) If all locations show the same singularity strength with α = constant or ′ α = constant, then Equation [14] becomes the same as the ordinary moving average function used in kriging and other methods (3) If the singularity varies from location to location, α ≠ constant or ′ α ≠ constant, then Equation [14] is equivalent to the ordinary moving average function multiplied by a scale ratio factor, (ε / ε max ) -Δα , with three possible situations given ε < ε max : [15] These properties indicate that if the data used for interpolation satisfies a multifractal distribution, then Equation [14] must be used as an extended form of the ordinary weighted averaging model. In this case, the scale ratio factor modifies the ordinary average in such a way that if there is positive singularity with ′ α > 0, then the new result is to be increased by a factor (Equation [15]), whereas if ′ α < 0, then the new result is reduced by a factor (Equation [15]). This modification is reasonable because ′ α > 0 and ′ α < 0 correspond to convex and concave properties of the surface ρ(ε, x 0 ) around the location x 0 , respectively.
The new model (Equation [14]) not only describes the spatial association reflected in the calculation of the weight λ, but also incorporates the singularity characterized by the singularity index α. The new model therefore has two obvious advantages: it not only improves the accuracy of the interpolated results but also retains the local structure of the interpolation map. The latter property is essential for geochemical and geophysical data processing and for pattern recognition. This will be demonstrated using the assay values from the Pulacayo sphalerite-quartz vein in Bolivia studied by De Wijs (1951).

Analysis of de Wijs's Bolivia sphalerite data
De Wijs (1951) studied assay values from the Pulacayo sphalerite-bearing quartz vein in Bolivia. Along a drift 118 channel samples had been collected at 2.00 m intervals ( Figure 1A). These values have been analysed by multifractal modeling and spatial analysis (Cheng and Agterberg, 1996;Cheng, 1997b) and they can be approximated by fiveparameter binomial multiplicative cascade models (Cheng, 1999c(Cheng, , 2014. The fractal dimension spectra of the distribution of de Wijs's zinc values are estimated by the gliding window method with order of moment ranging from -20 to 20 and cell size ε ranging from 2 m to 30 m (Cheng, 1999c) ( Figure 2) and fitted by the five-parameter binominal multiplicative cascade model (Cheng, 2014) (seen in  The results indicate that the distribution follows a multifractal model that is nearly symmetrical (Δα maxmin > 0).
In order to show the distribution of singularity and the data interpolation results, Figure 1B shows the resulting distribution of singularity values calculated for the data and the correlation coefficients associated with the linear model fitted after double log-transformation of measure and scale. It shows that the estimated values of α are within a range from 0.6 to 1.4 with correlation coefficients greater than 0.975 ( Figure 1C). Figure 2D illustrates the interpolated and reconstructed results obtained by means of MIM and moving averaging. The yellow line represents the results obtained using MIM with window size 20 m (10 point values) and the thicker red line represents the results obtained using the averaging technique with window size 6 m (two to three point values). The blue dots represent the observed data. Comparing the results obtained using MIM and moving averaging shows that MIM provides better results not only with smaller fitting errors for the observed data, but also that localized multifractality of the data is preserved.

Discussion and conclusions
It has been demonstrated that the multifractal distribution generated by binomial multiplicative cascade processes has multiple singularities that can be quantified by singularity index and fractal dimension spectrum. According to MIM, the singularity of multifractally distributed data can be used in data interpolation for mapping purposes with the localized structural properties (multifractality) preserved. The model used in MIM can be considered as an extended form of the ordinary moving average or weighted average used in various data interpolation methods, including inverse distance weighting and kriging. For most quantities in exploration geochemistry showing singularities, in order to retain the localized structural property, the multifractal interpolation method can be used to extend the ordinary moving average techniques, including ordinary kriging.
Since the singularity can be estimated using various methods, for example, integration of multiple patterns by weights of evidence method (Cheng, 2012) and other anisotropic cascade processes (Cheng, 2004), more general structural property and generalized self-similarity characterized by the singularity can be incorporated in the data interpolation. The multifractal interpolation based on prior knowledge and training images should be further explored