A batch sliding window method for local singularity mapping and its application for geochemical anomaly identification

Computers & Geosciences 90 (2016) 189–201

Contents lists available at ScienceDirect

Computers & Geosciences journal homepage: www.elsevier.com/locate/cageo

Case study

A batch sliding window method for local singularity mapping and its application for geochemical anomaly identification Fan Xiao a,b,n, Zhijun Chen c,d,nn, Jianguo Chen c,d, Yongzhang Zhou a,b a

School of Earth Science and Geological Engineering, Sun Yat-sen University, Guangzhou 510275, China Guangdong Key Laboratory of Geological Process and Mineral Resources Exploration, Sun Yat-sen University, Guangzhou 510275, China c State Key Laboratory of Geological Processes and Mineral Resources, China University of Geosciences, Wuhan 430074, China d Resources Faculty, China University of Geosciences, Wuhan 430074, China b

art ic l e i nf o

a b s t r a c t

Article history: Received 20 May 2015 Received in revised form 16 October 2015 Accepted 1 November 2015 Available online 2 November 2015

In this study, a novel batch sliding window (BSW) based singularity mapping approach was proposed. Compared to the traditional sliding window (SW) technique with disadvantages of the empirical predetermination of a fixed maximum window size and outliers sensitivity of least-squares (LS) linear regression method, the BSW based singularity mapping approach can automatically determine the optimal size of the largest window for each estimated position, and utilizes robust linear regression (RLR) which is insensitive to outlier values. In the case study, tin geochemical data in Gejiu, Yunnan, have been processed by BSW based singularity mapping approach. The results show that the BSW approach can improve the accuracy of the calculation of singularity exponent values due to the determination of the optimal maximum window size. The utilization of RLR method in the BSW approach can smoothen the distribution of singularity index values with few or even without much high fluctuate values looking like noise points that usually make a singularity map much roughly and discontinuously. Furthermore, the student's t-statistic diagram indicates a strong spatial correlation between high geochemical anomaly and known tin polymetallic deposits. The target areas within high tin geochemical anomaly could probably have much higher potential for the exploration of new tin polymetallic deposits than other areas, particularly for the areas that show strong tin geochemical anomalies whereas no tin polymetallic deposits have been found in them. & 2015 Elsevier Ltd. All rights reserved.

Keywords: Fractal/multifractals Uncertainty of singularity mapping Robust regression Tin polymetallic mineralization Gejiu, Yunnan

1. Introduction Singularity has been considered as one of the most important properties for non-linear natural processes or systems in various branches of earth science such as cloud formation (Schertzer and Lovejoy, 1985), flooding (Malamud et al., 1996), earthquakes (Turcotte, 1997), rainfall (Veneziano and Furcolo, 2002; Veneziano and Yoon, 2013), hurricanes (Sornette, 2004), landslides (Malamud et al., 2004) and hydrothermal mineralization (Cheng, 2007). From a viewpoint of geological application, it can be defined as a special phenomenon with anomalous amount of energy release or material accumulation occurring in narrow spatial-temporal intervals (Cheng, 2007). Taking metallic mineralization process as an n Corresponding author at: School of Earth Science and Geological Engineering, Sun Yat-sen University, Guangzhou 510275, China. nn Corresponding author at: State Key Laboratory of Geological Processes and Mineral Resources, China University of Geosciences, Wuhan 430074, China. E-mail addresses: [email protected] (F. Xiao), [email protected] (Z. Chen).

http://dx.doi.org/10.1016/j.cageo.2015.11.001 0098-3004/& 2015 Elsevier Ltd. All rights reserved.

example, it usually happens with a lot of metal elements precipitation and becomes to be enriched in relatively short geological time and small space scale of ore deposits. Singularity has considered to be a useful property in characterizing local depletion or enrichment of elements from a two dimensional geochemical map for identifying and extracting geochemical anomalies (Cheng, 2007). In practical application, to estimate the singularity exponent, a sliding window (SW) approach was originally developed by Cheng, 2006 on the basis of multifractal relationships between the density of metal and the corresponding area. Since the SW method based singularity mapping technique was proposed, it has been successfully proved to be a powerful tool for mineralization associated geochemical anomaly identification (e.g. Cheng, 2004, 2007, 2012, 2014; Chen et al., 2007; Cheng et al., 2009a, 2009b; Cheng and Agterberg, 2009; Arias et al., 2012; Zuo and Cheng, 2008; Zuo et al., 2009, 2015; Xiao et al., 2012, 2014; Zuo and Wang, 2015) and applied in a wide range of other geo-information anomaly recongnization associated fields such as geophysical exploration (Chen et al., 2015; Wang et al., 2013a), remote sensing (Cheng, 1999a; Neta et al., 2010), and tectonics data modeling

190

F. Xiao et al. / Computers & Geosciences 90 (2016) 189–201

List of notations V

μ(V) C(V) A μ(A) ρ(A)

α ε

Wþ W

Volume The total amount of metal in V The density of metal in V Area The total amount of metal in A The density of metal in A Local singularity index Window size Positive weight in the weights of evidence modeling Negative weight in the weights of evidence modeling

(Wang et al., 2012, 2013b; Zhao et al., 2014). However, on the one hand, the maximum window size in the SW approach, i.e. an important parameter used to define the largest scale where the power-law relationships probably exist between measure and scale for local singularity index estimation, should be empirically predetermined by applicants. On the other hand, the linear regression method used for fitting the power-law relationships on log-log plots is sensitive to outliers. These two issues significantly affect the estimation accuracy of singularity index values in practical applications (Cheng, 1999b; Chen et al., 2007; Zuo et al., 2013; Zhang et al., 2014). In order to overcome the disadvantages in traditional SW based singularity mapping technique, the purpose of this paper is mainly dressed on a procedure to improve SW approach in its accuracy of computation, and a batch SW (BSW) method was therefore developed in this study. It will show that the BSW method was a robust approach for singularity analysis. The concentration values of stream sedimentary geochemical data of tin, i.e. the predominant ore forming element of tin polymetallic mineralization

M Presence of a prospect Absence of a prospect M E Presence of an evidence theme Absence of an evidence theme E C Contrast of the weights of evidence t Student's t-statistic value s2(W +) The variances of W þ s2(W −) The variances of W  N (E ∩ M ) Numbers of unit cells or pixels of E N (E ∩ M ) Numbers of unit cells or pixels of E N (E ∩ M ) Numbers of unit cells or pixels of E N (E ∩ M ) Numbers of unit cells or pixels of E

modeling

∩M ∩M ∩M ∩M

in Gejiu distinct, Yunnan, China, were used to illustrate the implementation of the BSW based singularity analysis approach and a comparison study of BSW and SW in geochemical anomaly identification was conducted.

2. The SW approach The general principle and procedure of the SW approach can be summarized as follows. In the multifractal context, the singularity phenomenon of hydrothermal mineralization can be characterized by power-law models. For convenience but without losing generality, in a three dimensional space, assuming that μ(V) is the total amount of metal in volume V, and C(V), the density of metal in V is therefore equal to μ(V)/V. Obviously, μ(V) and C(V) accordingly change with V. If this relationship is multifractal, the power-law relationship between μ(V) or C(V) and V can be expressed as follows:

Fig. 1. The sketch diagram for illustrating SW method based singularity mapping technique.

F. Xiao et al. / Computers & Geosciences 90 (2016) 189–201

μ(V ) ∝ V α /3

(1)

C (V ) ∝ V α /3 − 1

(2)

where α is the local singularity index that is numerically equal to the exponent of the power-law relationship. Similarly, in a twodimensional space, assuming that μ(A) and ρ(A)¼ μ(A)/A are the total amount of metal and the density of metal within area A, respectively. Replace μ(V) and C(V) with μ(A) and ρ(A), respectively, Eqs. (1) and (2) can be updated to the following:

μ(A) ∝ Aα /2

(3)

ρ(A) ∝ Aα /2 − 1

(4)

To estimate α value in a given position of a two dimensional geochemical map, it defines a set of square windows with center fixed on this position and sizes εi ( εi = (2i − 1)εmin , εmin = ε1 <ε2 < ⋯ < εi⋯ < εn = εmax(i = 1, 2, 3, ⋯n) ), where n is the number of windows, εmin represents the smallest pixel size and εmax is the size of largest window, and calculates average concentration values ρ[A(εi )] within the corresponding areas A(εi ) = εi2. The log–log linear relationship between ρ[A(εi)] and εi derives from Eq. (5):

log ρ[A(εi )] = c + (α − 2)log εi

(5)

where c is a constant. On a log–log plot, the scatters of points (εi, ρ [A(εi)]) can be fitted by means of least-squares (LS), where the coefficient of determination (denoted R2) is used to indicate how well the data points fit a line and the slope is (α-2) (Fig. 1).

191

Similarly, by sliding windows from one position to another, singularity index values for a whole study area are calculated. The distribution of singularity index α in the concerned area can be characterized by the multifractal dimension spectrum function f (α), which implies that the majority values of α are close to 2 for normal area, while relatively minority values of α are larger or smaller than 2 for irregular or unusual areas (Cheng, 2007; Cheng and Agterberg, 2009). Commonly, the location where the singularity index α E2 means no geochemical singularity exists there. α o2 indicates a positive geochemical singularity and indicates anomalous enrichment of the element concentration, whereas α 42 suggests a negative geochemical singularity and may suggest depletion of element concentration (Fig. 1). Many other publications (Cheng, 2004, 2007, 2012; Chen et al., 2007; Cheng and Agterberg, 2009) have summarized more properties of local singularity index α.

3. The BSW method In terms of multifractals, singularity index α also known as Hölder exponent, is essentially the fractal dimension of a measure (e.g. average concentration of geochemical element) within the measure unit (e.g. local range of εmax  εmax in Fig. 1), where the measure is considered to be mono-fractal (Evertsz and Mandelbrot, 1992; Mandelbrot, 1989). Hence, calculation accuracy and uncertainty of the estimated singularity index are strongly dependent on the size of measure unit, i.e. the value of the largest window size εmax. This is because that the fractal relationships in Eqs. (3) and (4) are commonly not applied for any scale of εmax, which has been proved in some practical

raster geochemical map (x, y) singularity estimation for position (η=1, γ=1)(1≤η≤x, 1≤γ≤y) produce m batchs of squre windows with szie ε i (i =1, 2, 3, • • •, n) j) ε min = ε1 < ε 2 < ε 3 < • • • < ε i • • • < ε n = ε(max ( j =1, 2, 3,• • •, m)

j=j+1

log-log plot and RLR fitting of (εi , ρ[A(εi )]) calculate singularity index α(j) and value of RESM(j) yes

j ≤ m? no

singularity index α(η, γ)=α( j) if RESM(j)=min{ RESM } η ≤ x | γ≤ y ?

yes

no singularity index values of geochemical elements Fig. 2. The procedures of BSW method based singularity mapping approach.

η=η+1, γ=γ+1

calculate ρ[A(ε i )]

192

F. Xiao et al. / Computers & Geosciences 90 (2016) 189–201

applications (Cheng, 1999b; Chen et al., 2007; Zuo et al., 2013; Zhang et al., 2014). However, in the practical application of the SW method based singularity mapping technique, εmax is a fixed value that should be empirically predetermined by applicants (Cheng, 2000). Though, theoretically, εmax can be any value that equals to (2n − 1)εmin(n =2, 3, 4, ⋯) , it is difficult or even impossible to determine a reasonable and fixed value for the calculation of singularity index values for all positions in the whole area. This reason is that, to specific measures (e.g. average concentration of geochemical elements), the fractal relationships in Eqs. (3) and (4) probably exist in distinct scales of εmax. In other words, εmax should be flexible to “best” fit the fractal prosperity of a local data set within neighborhood (εmax  εmax) of each location rather than a fixed value for all locations. Thus, in order to improve the calculation accuracy and minimize the uncertainty of singularity index, a BSW method was proposed and provided a robust approach to singularity analysis. In this new method, instead of LS, robust linear regression (RLR) method (Holland and Welsch, 1977; Huber, 1981), which is insensitive to outlier values was employed to fit the scatters of points (εi, ρ[A(εi)]) on log–log plots. The root mean squared error (RMSE) value was used to estimate how well the (j ) data points within various possible largest window sizes εmax ( j = 1, 2, 3, ⋯m) fit straight lines. Therefore, the optimize value of largest window size corresponding to the minimum value of RMSE can be determined for each estimated location. Suppose that a two dimensional rater geochemical map has been created based on element concentration values of geochemical samples in a specific study area. The procedures and practical implementation of BSW based singularity mapping approach can be summarized as follows (Fig. 2): (1) define m batches of square windows with centers located at the estimated position and sizes equal to εi (( j = 1, 2, 3, ⋯m)), (j ) where εmin represents the smallest pixel size, n and εmax are the total number of and the size of largest window in the jth (j ) set of windows, respectively. In theory, the value of εmax equals to 2(n − 1)εmin , where n is usually an integer number larger than 2 to ensure at least three scatters in the following log–log plot so as to make the LS regression there much more statistically significant. (2) estimate the average concentration values ρ[A(εi)] within the corresponding areas A(εi) ¼ εi2 in the jth set of windows and log–log plot (εi, ρ[A(εi)]). (3) using RLR method to fit the scatters of (εi, ρ[A(εi)]) based on linear relationship, i.e. log ρ[A(εi )] = c + (α − 2)log εi , to derive the slope of the line that (α-2) and the value of RMSE. (4) repeat steps (2)–(3) until j4m, which means all the produced m sets of windows have been estimated. (5) find the minimum value of RMSE derived from step (4). The corresponding value of α and maximum window size are regarded as the singularity index and the optimal size of the largest window for singularity analysis of the estimated location, respectively. (6) move to the next position, repeat steps (1)–(5) until the singularity index values for all possible locations within the studied geochemical map have been obtained.

4. Robust linear regression In terms of robust statistics, robust linear regression (RLR) is a form of regression analysis approach designed to circumvent some limitations of classical linear model (i.e. the method of least-squares (LS)), where the regression error is usually assumed to be normally distributed with mean zero and standard deviation (Holland and Welsch, 1977; Huber, 1981; Maronna et al., 2006). Commonly, the LS regression method could have favorable properties if there is no outliers in

analyzed data and its underlying assumptions are true. However, LS estimates for regression models are highly sensitive to outliers (Huber, 1981; Maronna et al., 2006). In particular, if the observed data prone to outliers, the model assumptions of LS are invalidated, and parameter estimates, confidence intervals, and other computed statistics probably become unreliable (Huber, 1981; Maronna et al., 2006). In Fig. 3, we illustrated and compared the estimation results of the LS regression model and RLR method using a modeling data set with outliers (i.e. the second and the last data points that are long distance away from the linear trend of the other data). It shows clearly that the RLR fit is much more less influenced by the outliers than the LS fit. Indeed, RLR has became a cornerstone of robust statistics due to it is not overly affected by violations of regression assumptions by the underlying data-generating process (Maronna et al., 2006). Its principle and computation methods have been widely and detailed documented in many textbooks (e.g. Huber, 1981; Zhou and Huang, 1989; Maronna et al., 2006). In the singularity mapping analysis above mentioned, the linear model of Eq.(5) is usually influenced by outliers due to it derives from the empirical assumptions of fractal relationships in Eq. (1), Eq. (2), Eq. (3) and Eq. (4). This issue has been practically proved to be one of important factors that could affect the estimation accuracy of singularity index values (Cheng, 1999b; Chen et al., 2007; Zuo et al., 2013, 2015). Thus, in order to decrease the influences of outliers in singularity mapping, it is probably more rescannable to use RLR fit rather than LS fit.

5. Students' t-statistic The Student's t-statistic is an useful measure in the weights of evidence (WofE) modeling approach, which not only can be used to estimate the spatial association between point patterns of prospects and evidence themes, but also to reclassify the evidence layers for optimal prediction based on spatial association analysis (Agterberg et al., 1990; Bonham-Carter, 1994; Carranza, 2004; Cheng, 2007; Zuo et al., 2009; Xiao et al., 2012, 2014; Zhao et al., 2012). In the WofE model, there are two important statistics that positive weight (W þ ) and negative weight (W  ) defined to estimate the spatial association between prospects (e.g. known deposits or occurrences) and evidence themes (e.g. geochemical and geophysical anomaly maps) based on conditional probabilities:

W + = ln

P (E/M ) P (E/M )

(6)

W − = ln

P (E /M ) P (E /M )

(7)

where M and M refer to presence and absence of a prospect, respectively; and E and E represent presence and absence of an evidence theme. W þ measures the spatial association between prospects and the presence of evidence theme while W- measures the spatial association between prospects and its absence. The difference C ¼W þ  W  termed the contrast provides a measure of overall spatial association between evidence themes and prospects. Furthermore, in order to evaluate the statistical significance of spatial relationship, Student's t-statistic is defined with:

t=

C = s(C )

C s 2(W +) + s 2(W −)

(8)

where s2(W +) and s2(W −) are the variances of W þ and W  , respectively. These are computed as:

s 2(W +) =

1 1 + N (E ∩ M ) N (E ∩ M )

(9)

F. Xiao et al. / Computers & Geosciences 90 (2016) 189–201

s 2(W −) =

1 1 + N (E ∩ M ) N (E ∩ M )

(10)

where N (E ∩ M ), N (E ∩ M ), N (E ∩ M ), and N (E ∩ M ) represent numbers of unit cells or pixels of E ∩ M , E ∩ M , E ∩ M and E ∩ M on map patterns. Usually, estimates of t-value greater than 1.96 suggests that the spatial association is statistically significant at level of significance α ¼0.05 under some conditions (Agterberg et al., 1990; Bonham-Carter, 1994). Many publications have provided more detailed discussions on this method and derivation of the above equations (e.g. Agterberg et al., 1990; Bonham-Carter, 1994).

6. Case study 6.1. Geological background and geochemical data sets The study area used as example for this study is the Gejiu mineral district, which is located south of and about 200 km far from Kunming, the capital of Yunnan Province, southwest China. It is well known as one of the most important tin mineralization and production areas all over the world because it hosts several world20 Data Ordinary Least Squres Robust Regression

10

0

-10

-20

-30

-40

-50

-60

0

2

4

6

8

10

12

14

16

18

20

Fig. 3. Illustration and comparison for the estimation results of the LS regression model and RLR method.

193

class tin polymetallic deposits, which have been exploited in the past several decades. A simplified geological map for the study area is shown in Fig. 4. The geological units predominately consist of a sequence of Paleozoic to Mesozoic sedimentary strata such as the middle Triassic carbonate rocks of the Gejiu Formation, which is the main country rock hosting most of the known tin polymetallic deposits and igneous rocks including volcanic and intrusive rocks. The Gejiu Batholith, a well exposed granitoid complex rich in Si, K/Na and Al in the center of the study area, has been consensually believed to be one of the most important controlling lithologies for the tin polymetallic mineralization in this area and most of the discovered tin polymetallic deposits are located near it (Zhuang et al., 1996). The N–S trend Gejiu Fault divides the study area into two parts, named eastern and western Gejiu district. The main faults in the eastern district have N–S and E–W trending, whereas the main orientations of faults in the western district are NE–SW or NW–SE trending. These fault systems control the general structure of the mineralization and distribution of tin polymetallic deposits in the study area. Several large tin polymetallic deposits, including Malage, Gaosong, Laochang and Kafang, have been found in the eastern district and are concentrated at the intersections of NNE–SSW and E–W trending faults. In the current study, a total of 1921 stream sediment samples were collected and analyzed by the Chinese National Geochemical Mapping Project as part of the Regional Geochemistry National Reconnaissance (RGNR) Project, which has been carried out since 1979 and covered most of the land surface of China including Gejiu distinct. The sampling density is one site per km2 and four samples per 2 km  2 km cells were mixed to form composite samples. For each simple, the concentration of 39 elements and major oxides including Ag, As, Au, B, Ba, Be, Be, Bi, Cd, Co, Cr, Cu, F, Hg, Mo, Pb, Sb, Sn, Fe2O3, CaO, MnO and SiO2 were analyzed by means of a multi-element analytical system comprising ICP-MS (inductively coupled plasma-mass spectrometry), XRF (X-ray fluorescence) and ICP-AES (inductively coupled plasma-atomic mission spectrometry) as the backbone in combination with other methods (Xie et al., 1997). The data used for this study is the concentration value of tin, i.e. the dominant ore element of tin polymetallic deposits. More details about the sampling and analysis of the stream sediment data have been previously summarized by Xie et al., 1997. 6.2. Tin polymetallic mineralization associated geochemical anomaly identification with the BSW based singularity mapping approach Based on the raster map of tin concentration values obtained by

Fig. 4. Simplified geological map of the study area in Gejiu, Yunnan, China (after the geological map produced by the Yunnan Geological Survey at1:2, 500,000 scale).

194

F. Xiao et al. / Computers & Geosciences 90 (2016) 189–201

Fig. 5. Map showing the spatial distribution of tin concentration values in 1921 stream sediment samples using inverse distance weighted moving average method.

Fig. 6. Map showing the spatial distribution of singularity index (α) values of tin calculated by BSW method (The plus dots are the locations of the examples of estimation of local tin singularity in Fig.7).

inverse distance weighted interpolation method (Fig. 5), the BSW based singularity mapping approach with maximum window sizes (1) (1) (13) (j ) (j ) ( εmax ) ranged from εmax ¼5 km to εmax ¼29 km ( εmax ¼ εmax þ 2  (j1); j¼ 1, 2, …, 13) was used to calculate singularity index values for tin in the studied area. The local singularity map of tin is shown in Fig. 6, where most known tin polymetallic deposits are spatially correlated with the areas of singularity index (α) value less than 2 (α o2). Also,. in order to compare the singularity mapping results calculated by BSW method with SW approach and using RLR instead of LS, the singularity maps shown in Fig. 7 are derived using SW approach and the approach of BSW procedure using LS regression (i.e. the calculation procedure is the same as the proposed BSW method, where only the RLR is replaced with LS in the step (3)). Fig. 7A–C are the singularity mapping results using the traditional BW approach with maximum window size of 7 km, 17 km, and 27 km, respectively. Fig. 7D is the singularity mapping result obtained by the approach of BSW procedure using LS regression. Compared Fig.6 with Fig. 7, it shows that the distribution of singularity index values computed with BSW method (Fig. 6) is smooth with few or even without much high fluctuate values looking like noise points, which typically make a singularity exponent map much roughly and discontinuously since the maximum window size is not reasonable or/and LS regression method is sensitive to outlier values in SW based singularity mapping approach (Fig. 7).

To illustrate the determination of optimal maximum window sizes for singularity value estimation in BSW based singularity mapping method, 13 different locations were selected on the map (j ) as shown on Fig. 6, of which the best maximum window sizes εmax … … (j¼1, 2, , 13) are equal to 5 km, 7 km, 9 km, , 29 km. In Fig. 8, it illustrates in detail how to determine optimal maximum window sizes at the 13 different locations. The corresponding maximum window size and α-value are respectively regarded as the singularity index and the optimal size of the largest window for singularity analysis of the estimated location while minimum value of RMSE is derived. Furthermore, student's t-statistic (Agterberg et al., 1990; BonhamCarter, 1994) was utilized to measure the statistic significance of the strength of spatial correlation between point pattern of the known tin polymetallic deposits and the local tin singularity maps obtained by using BSW method, SW approach and the approach of BSW procedure using LS regression. The student's t-value diagrams were shown in Fig. 9. Fig. 9A is the student's t statistic diagram of singularity map calculated by the proposed BSW method. Fig. 9B–D are the student's t statistic diagrams of singularity maps derived by the traditional BW approach with maximum window size of 7 km, 17 km, and 27 km, respectively. Fig. 9E is the student's t statistic diagrams of singularity map computed by the approach of BSW procedure using LS regression. The maximum t values in Fig. 9A–E are respectively equal to 4.73, 2.48, 3.98, 4.15 and 4.06, which means that the singularity map

F. Xiao et al. / Computers & Geosciences 90 (2016) 189–201

195

Fig. 7. Map showing the spatial distribution of singularity index (α) values of tin calculated by SW approach and by BSW procedure using LS regression (A-maximum window size is 7 km; B-maximum window size is 17 km; C-maximum window size is 27 km; D-BSW procedure using LS regression).

196

F. Xiao et al. / Computers & Geosciences 90 (2016) 189–201

0.1

2.05

0.08

2.00

0.06

1.95

α-value

RMSE

location 1

0.04 0.02

1.90 1.85

location 1 0

1.80 5 7 9 11 13 15 17 19 21 23 25 27 29 Maximum window size

0.35 0.30

5 7 9 11 13 15 17 19 21 23 25 27 29 Maximum window size 3.00

location 2

2.50 2.00 α-value

RMSE

0.25 0.20 0.15 0.10

1.50 1.00 location 2 0.50

0.05

0.00

0.00

5 7 9 11 13 15 17 19 21 23 25 27 29 Maximum window size

5 7 9 11 13 15 17 19 21 23 25 27 29 Maximum window size 0.35

3.00

0.30

2.50 α-value

RMSE

0.25 0.20 0.15 0.10 0.05

location 3

2.00 1.50 1.00 0.50 location 3

0.00

0.00 5 7 9 11 13 15 17 19 21 23 25 27 29 Maximum window size

0.35

3.00

0.30

2.50

0.25

2.00

α-value

RMSE

5 7 9 11 13 15 17 19 21 23 25 27 29 Maximum window size

0.20 0.15 0.10 0.05

location 4

0.00

1.50 1.00 0.50

location 4

0.00 5 7 9 11 13 15 17 19 21 23 25 27 29 Maximum window size

5 7 9 11 13 15 17 19 21 23 25 27 29 Maximum window size

Fig. 8. Examples of determination optimize size of largest windows for local singularity estimation at locations shown in Fig. 5. Plots on left and right show relations between values RMSE and α and maximum window size εmax, respectively. εmax is length of largest side of square cells with widths set equal to 5, 7, 9,…, 29 km. RMSE is root mean squared error value that is used to estimate how well scatters of points (εi, ρ[A(εi)]) fit a straight line on log–log plot using RLR method, where εi is window size that is not larger than εmax and ρ[A(εi)] is mean of concentration values for all 1  1 km cells contained within εi centered on a location. The exponent α is estimated as singularity index value derived from the slope of fitted straight-line between εi and ρ[A(εi)] plus to 2. The corresponding maximum window size and α-value will be regarded as the singularity index and the optimal size of largest window for singularity analysis of the estimated location while minimum value of RMSE is derived.

0.06

2.05

0.05

2.00

0.04

location 9 α-value

RMSE

F. Xiao et al. / Computers & Geosciences 90 (2016) 189–201

0.03 0.02 0.00

1.95 1.90

1.80 5 7 9 11 13 15 17 19 21 23 25 27 29 Maximum window size

0.05

2.20 2.15 α-value

RMSE

5 7 9 11 13 15 17 19 21 23 25 27 29 Maximum window size

location 10

0.04 0.03 0.02

2.10 2.05

0.01

2.00

0.00

1.95

0.050

α-value

RMSE

location 11

0.030 0.020

0.010 0.000

1.905 1.900 1.895 1.890 1.885 1.880 1.875 1.870

2.25 α-value

RMSE

2.30

location 12

2.20 2.15 2.10 2.05

5 7 9 11 13 15 17 19 21 23 25 27 29 Maximum window size

location 13 α-value

RMSE

location 12

2.00 5 7 9 11 13 15 17 19 21 23 25 27 29 Maximum window size

0.80 0.70 0.60 0.50 0.40 0.30 0.20 0.10 0.00

location 11

5 7 9 11 13 15 17 19 21 23 25 27 29 Maximum window size

5 7 9 11 13 15 17 19 21 23 25 27 29 Maximum window size 0.08 0.07 0.06 0.05 0.04 0.03 0.02 0.01 0.00

location 10 5 7 9 11 13 15 17 19 21 23 25 27 29 Maximum window size

5 7 9 11 13 15 17 19 21 23 25 27 29 Maximum window size

0.040

location 9

1.85

0.01

0.06

197

3.50 3.00 2.50 2.00 1.50 1.00 0.50 0.00

5 7 9 11 13 15 17 19 21 23 25 27 29 Maximum window size Fig. 8. (continued)

location 13 5 7 9 11 13 15 17 19 21 23 25 27 29 Maximum window size

198

F. Xiao et al. / Computers & Geosciences 90 (2016) 189–201

0.16

2.15

0.14

2.10 location 5

0.12

α -value

RMSE

0.10

2.05

0.08 0.06

1.85 1.80 5

7

5

α -value

0.15 0.10 0.05 location 6 0.00 5

7

2.50 2.45 2.40 2.35 2.30 2.25 2.20 2.15 2.10 2.05

9 11 13 15 17 19 21 23 25 27 29 Maximum window size

0.35

4.00

0.30

3.50

0.25

3.00 α -value

0.20 0.15 0.10

0.00

5

7

5

7

2.00 1.50

0.00

location 8

0.030

α -value

0.025 0.020 0.015 0.010 0.005 5

7 9 11 13 15 17 19 21 23 25 27 29 Maximum window size

9 11 13 15 17 19 21 23 25 27 29 Maximum window size

2.50

9 11 13 15 17 19 21 23 25 27 29 Maximum window size

0.035

9 11 13 15 17 19 21 23 25 27 29 Maximum window size

location 6

location 7

0.50

0.040

0.000

7

1.00

location 7

0.05

location 5

1.75

9 11 13 15 17 19 21 23 25 27 29 Maximum window size

0.20 RMSE

1.90

0.02

0.25

RMSE

1.95

0.04 0.00

RMSE

2.00

5

7

2.20 2.18 2.16 2.14 2.12 2.10 2.08 2.06 2.04 2.02

9 11 13 15 17 19 21 23 25 27 29 Maximum window size

location 8 5

7

9 11 13 15 17 19 21 23 25 27 29 Maximum window size

Fig. 8. (continued)

calculated by the BSW method could have the larger maximum t value than the other three methods. In other words, the geochemical anomaly identified by the proposed BSW based singularity mapping

technique could have strongest spatial correlation with the known tin polymetallic deposits. It indicates that BSW method can effectively optimize the estimation of singularity exponent index values.

6.00

5.00

5.00

4.00

4.00

3.00

t-value

t-value

F. Xiao et al. / Computers & Geosciences 90 (2016) 189–201

3.00 2.00

2.00

1.00

1.00 0.00 0.60

199

0.00

0.80

1.00

1.20

1.40

1.60

1.80

2.00

2.20

-1.00 0.60 0.80 1.00 1.20 1.40 1.60 1.80 2.00 2.20 2.40 2.60 2.80 3.00 3.20

2.40

singularity index (α)

singularity index (α)

4.50

2.00

3.50

1.00

2.50

t-value

t-value

3.00

0.00

1.50

-1.00

0.50

-2.00 0.60

1.00

1.40

1.80

2.20

2.60

3.00

3.40

3.80

4.20

-0.50 1.00

4.60

1.20

1.40

1.60

singularity index (α)

1.80

2.00

2.20

2.40

2.60

2.80

3.00

singularity index (α)

5.00

4.00

t-value

3.00

2.00

1.00

0.00

-1.00 1.00

1.20

1.40

1.60

1.80

2.00

2.20

2.40

2.60

2.80

3.00

3.20

singularity index (α)

Fig. 9. Student's t-value for singularity index of tin (A-BSW method; B-SW approach with maximum window size is 7 km; C-SW approach with maximum window size is 17 km; D- SW approach with maximum window size is 27 km; E-BSW procedure using LS regression).

Fig. 10. Map showing target areas for prospecting potential of tin polymetallic deposits.

200

F. Xiao et al. / Computers & Geosciences 90 (2016) 189–201

Finally, we use the singularity map derived by BSW method to delineate potential tin polymetallic mineralization areas. In Fig. 9A, it shows that, although there are relatively few known tin polymetallic deposits, the largest t-value 4.73 (where the corresponding α value equals to 1.89 far exceeds to t0.05 ¼2.0 representing 95% confidence level for statistical significance) implies strong positive spatial correlation between the two patterns. In total, 595 local tin singularities have α o1.89 and their combined areas occupy only 7.68% of the whole study area, but 54.2% of the known tin polymetallic deposits occur within this relatively small low-singularity sub-area. The singularity exponent value 1.89 can be regarded as a threshold value to classify the tin geochemical anomaly values into two intervals (i.e. α o 1.89 and α Z1.89), corresponding to high and low anomalies, respectively. The classification results using these two intervals are shown in Fig. 10. It confirms that the target areas within high tin geochemical anomaly could probably have much higher potential to be explored for new tin polymetallic deposits than other areas, particularly for the areas that show strong tin geochemical anomalies whereas no tin polymetallic deposits have been found in them.

7. Conclusions In this paper, a batch sliding window based singularity mapping method has been proposed. The new approach could automatically determine optimal size of largest window and uses robust linear regression that is insensitive to outlier values for local singularity mapping. Therefore, it can overcome the disadvantages in traditional SW based singularity mapping technique: (1) a fixed maximum window size should be empirically predetermined by applicants; and (2) least-squares (LS) linear regression method that is used for fitting power-law relationships on log–log plots is sensitive to outliers. The BSW approach is a useful tool for singularity mapping. It can not only improve calculation accuracy for singularity exponent value due to it determines the optimal maximum window size automatically, but also smoothen the distribution of singularity index values since it utilizes the RLR method for linear regression that is not sensitive to outlier values. In the case study, tin geochemical data have been processed by BSW based singularity mapping approach. The student's t-statistic diagram indicates that there is a strong spatial correlation between high geochemical anomaly and known tin polymetallic deposits. The target areas within high tin geochemical anomaly could probably have much higher potential to explore new tin polymetallic deposits than the other areas, particularly for the areas showing strong tin geochemical anomalies whereas no tin polymetallic deposits have been found in them.

Acknowledgments This research was financially supported by the National Natural Science Foundation of China (Nos. 41502310, 41272361) and the Natural Science Foundation of Guangdong Province, China (No. 2015A030310246).

References Agterberg, F.P., Bonham-Carter, G.F., Wright, D.F., 1990. Statistical pattern integration for mineral exploration. In: Gaál, G., Merriam, D.F. (Eds.), Computer Applications in Resource Estimation. Pergamon Press, Oxford, pp. 1–21. Arias, M., Gumiel, P., Martin-Izard, A., 2012. Multifractal analysis of geochemical anomalies: a tool for assessing prospectivity at the SE border of the Ossa Morena Zone, Variscan Massif (Spain). J. Geochem. Explor. 122, 101–112. Bonham-Carter, G.F., 1994. Geographic Information Systems for Geoscientists:

Modeling with GIS. Pergamon, New York. Carranza, E.J.M., 2004. Weights of evidence modeling of mineral potential: a case study using small number of prospects, Abra, Philippines. Nat. Resource Res. 13, 173–187. Chen, G.X., Cheng, Q.M., Zuo, R.G., Liu, T.Y., Xi, Y.F., 2015. Identifying gravity anomalies caused by granitic intrusions in Nanling mineral district, China: a multifractal perspective. Geophys. Prospect. 63, 256–270. Chen, Z.J., Cheng, Q.M., Chen, J.G., Xie, S.Y., 2007. A novel iterative approach for mapping local singularities from geochemical data. Nonlinear Process. Geophys. 14, 317–324. Cheng, Q.M., 1999a. Multifractality and spatial statistics. Comput. Geosci. 25, 949–961. Cheng, Q.M., 1999b. Spatial and scaling modelling for geochemical anomaly separation. J. Geochem. Explor. 65, 175–194. Cheng, Q.M., 2000. GeoData AnalysisSystem (GeoDAS) for Mineral Exploration: Unpublished User's Guide and Exercise Manual, Material for the Training Workshop on GeoDAS held at York University. Cheng, Q.M., 2004. Quantifying the generalized self-similarity of spatial patterns for mineral resource assessment. Earth Sci. China Univ. Geosci. 29, 733–743 (In Chinese with English abstract). Cheng, Q.M., 2006. GIS-based multifractal anomaly analysis for prediction of mineralization and mineral deposits. In: Harris, J. (Ed.), GIS Applications in Earth Sciences, Geological Association of Canada Special Paper, pp. 289–300. Cheng, Q.M., 2007. Mapping singularities with stream sediment geochemical data for prediction of undiscovered mineral deposits in Gejiu Yunnan Province, China. Ore Geol. Rev. 32, 314–324. Cheng, Q.M., 2012. Singularity theory and methods for mapping geochemical anomalies caused by buried sources and for predicting undiscovered mineral deposits in covered areas. J. Geochem. Explor. 122, 55–70. Cheng, Q.M., 2014. Vertical distribution of elements in regolith over mineral deposits and implications for mapping geochemical weak anomalies in covered areas. Geochem.-Explor. Environ. Anal. 14, 277–289. Cheng, Q.M., Agterberg, F.P., 2009. Singularity analysis of ore-mineral and toxic trace elements in stream sediments. Comput. Geosci. 35, 234–244. Cheng, Q.M., Zhao, P.D., Chen, J.G., Xia, Q.L., Chen, Z.J., Zhang, S.Y., Xu, D.Y., Xie, S.Y., Wang, W.L., 2009a. Application fo sigularity theory in prediction of tin and copper mineral deposits in Gejiu district, Yunnan, China: Weak information extraction and mixing information decomposition. Earth Sci.-J. China Univ. Geosci. 34, 232–242 (In Chinese with English abstract). Cheng, Q.M., Zhao, P.D., Zhang, S.Y., Xia, Q.L., Chen, Z.J., Chen, J.G., Xu, D.Y., Wang, W. L., 2009b. Application of singularity theory in prediction of tin and copper mineral deposits in Gejiu district, Yunnan, China: Information integration and delineation of mineral exploration targets. Earth Science-Journal of China University of Geosciences 34; b, pp. 243–252 (In Chinese with English abstract). Evertsz, C.J.G., Mandelbrot, B.B., 1992. Multifractal measures. In: Peitgen, H.-O., Jürgens, D., Saupe, D. (Eds.), Chaos and Fractals. Springer-Verlag, New York, p. 984, p. Holland, P.W., Welsch, R.E., 1977. Robust regression using iteratively reweighted least-squares. Commun. Stat.: Theory Methods A6, 813–827. Huber, P.J., 1981. Robust Statistics. John Wiley & Sons, Inc., Hoboken, NJ. Malamud, B.D., Turcotte, D.L., Barton, C.C., 1996. The 1993 Mississippi River flood: a one hundred or a one thousand year event? Environ. Eng. Geosci. II, 479–486. Malamud, B.D., Turcotte, D.L., Guzzetti, F., Reichenbach, P., 2004. Landslide inventories and their statistical properties. Earth Surf. Process. Landforms 29, 687–711. Mandelbrot, B.B., 1989. Multifractal measures, especially for the geophysicist. In: Scholz, B., Mandelbrot, B. (Eds.), Fractals in Geophysics. Birkhäuser, Basel, pp. 5–42. Maronna, R.A., Douglas, M.R., Yohai, V.J., 2006. Robust statistics: theory and methods. John Wiley & Sons Ltd, West Sussex. Neta, T., Cheng, Q.M., Bello, R.L., Hu, B., 2010. Upscaling reflectance information of lichens and mosses using a singularity index: a case study of the Hudson Bay Lowlands Cnada. Biogeosciences 7, 2557–2565. Schertzer, D., Lovejoy, S., 1985. The dimension and intermittency of atmospheric dynamics-multifractal cascade dynamics and turbulent intermittency. In: Launder, B. (Ed.), Turbulent Shear Flow. Springer, New York, pp. 7–33. Sornette, D., 2004. Critical Phenomena in Natural Sciences: Chaos, Fractals, Selforgnization and Disorder: Concepts and Tools, second ed. Springer, New York. Turcotte, D.L., 1997. Fractals and Chaos in Geology and Geophysics, second ed. Cambridge Universtity Press, Cambridge. Veneziano, D., Furcolo, P., 2002. Multifractality of rainfall and scaling of intensityduration-frequency curves. Water Resource Res. 38, 1306. Veneziano, D., Yoon, S., 2013. Rainfall extremes, excesses, and intensity-durationfrequency curves: a unified asymptotic framework and new nonasymptotic results based on multifractal measures. Water Resource Res. 49, 4320–4334. Wang, W.L., Zhao, J., Cheng, Q.M., 2013a. Application of singularity index mapping technique to gravity/magnetic data analysis in southeastern Yunnan mineral district, China. J. Appl. Geophys. 92, 39–49. Wang, W.L., Zhao, J., Cheng, Q.M., 2013b. Fault trace-oriented singularity mapping technique to characterize anisotropic geochemical signatures in Gejiu mineral district, China. J. Geochem. Explor. 134, 27–37. Wang, W.L., Zhao, J., Cheng, Q.M., Liu, J.T., 2012. Tectonic-geochemical exploration modeling for characterizing geo-anomalies in southeastern Yunnan district, China. J. Geochem. Explor. 122, 71–80. Xiao, F., Chen, J.G., Agterberg, F., Wang, C.B., 2014. Element behavior analysis and its implications for geochemical anomaly identification: a case study for porphyry

F. Xiao et al. / Computers & Geosciences 90 (2016) 189–201

Cu-Mo deposits in Eastern Tianshan, China. J. Geochem. Explor. 145, 1–11. Xiao, F., Chen, J.G., Zhang, Z.Y., Wang, C.B., Wu, G.M., Agterberg, F.P., 2012. Singularity mapping and spatially weighted principal component analysis to identify geochemical anomalies associated with Ag and Pb–Zn polymetallic mineralization in Northwest Zhejiang, China. J. Geochem. Explor. 122, 90–100. Xie, X.J., Mu, X.Z., Ren, T.X., 1997. Geochemical mapping in China. J. Geochem. Explor. 60, 99–113. Zhang, D.J., Cheng, Q.M., Agterberg, F.P., 2014. An Application of Equal-AreaGrowing Window for Calculating Local Singularity for Mapping Granites in Inner Mongolia. In: Pardo-Igúzquiza, E., Guardiola-Albert, C., Heredia, J., Moreno-Merino, L., Durán, J.J., Vargas-Guzmán, A.J. (Eds.), Mathematics of Planet Earth: Proceedings of the 15th Annual Conference of the International Association for Mathematical Geosciences. Springer Berlin Heidelberg, Berlin, Heidelberg, pp. 73–77. Zhao, J., Wang, W.L., Dong, L.H., Yang, W.Z., Cheng, Q.M., 2012. Application of geochemical anomaly identification methods in mapping of intermediate and felsic igneous rocks in eastern Tianshan, China. J. Geochem. Explor. 122, 81–89. Zhao, J.N., Zuo, R.G., Chen, S.Y., Kreuzer, O.P., 2014. Application of the tectonogeochemistry method to mineral prospectivity mapping: A case study of the Gaosong tin-polymetallic deposit, Gejiu district, SW China. Ore Geol. Rev. 71,

201

719–734. Zhou, F.G., Huang, Y.C., 1989. Applied Linear Regression Analysis. China Renmin University Press, Beijing. Zhuang, Y.Q., Wang, R.Z., Yang, S.P., 1996. Geology of Gejiu Tin-Copper Polymetallic Deposite. Earthquake Publishing House, Beijing, China. Zuo, R.G., Cheng, Q.M., 2008. Mapping singularities-a technique to identify potential Cu mineral deposits using sediment geochemical data, an example for Tibet, west China. Miner. Mag. 72, 531–534. Zuo, R.G., Cheng, Q.M., Agterberg, F.P., Xia, Q.L., 2009. Application of singularity mapping technique to identify local anomalies using stream sediment geochemical data, a case study from Gangdese, Tibet, western China. J. Geochem. Explor. 101, 225–235. Zuo, R.G., Wang, J., 2015. Fractal/multifractal modeling of geochemical data: A review. Journal of Geochemical Exploration, doi.org/10.1016/j.gexplo.2015.04.010. Zuo, R.G., Wang, J., Chen, G.X., Yang, M.G., 2015. Identification of weak anomalies: a multifractal perspective. J. Geochem. Explor. 148, 12–24. Zuo, R.G., Xia, Q.L., Zhang, D.J., 2013. A comparison study of the C-A and S-A models with singularity analysis to identify geochemical anomalies in covered areas. Appl. Geochem. 33, 165–172.