Quantitative 3D Association of Geological Factors and Geophysical Fields with Mineralization and Its Significance for Ore Prediction: An Example from Anqing Orefield, China

Abstract

Quantitative 3D spatial association of geological factors and geophysical fields with orebodies is critical for ore prediction. The Anqing orefield, a principal Cu–Fe orefield in China, is closely associated with the Yueshan intrusion. By compiling the data from drilling and tunnelling exploration, Controlled Source Audiofrequency Magnetotelluric (CSAMT) surveying and the computational modelling of magmatic intrusion’s cooling process, we constructed models of the Yueshan intrusion, ore-favourable carbonate formation, orebodies, resistivity field and volume strain field. These models are used as evidential factors to analyse their spatial association with mineralization by the weights-of-evidence (WofE) method. The location of orebodies is closely related to the shape of the contact zone of the Yueshan intrusion. The spaces with the distance ≤200 m to the concaves that were selected by minimum principal curvature (|K_min| ≥ 0.0025) from contact zones, are very favourable for localization of orebodies. Most orebodies are not located in the spaces of the lowest resistivity, suggesting that the lowest resistivity cannot be used as an indicator for mineralization. The spaces with higher positive volumetric strain have higher positive weights with orebodies, implying that the mineralization is positively related to the positive volumetric strain. The spaces of all evidential factors that had positive correlation with mineralization were integrated to create a 3D prospectivity map by calculating posterior probability. Five areas with higher posterior probability, indicating higher prospectivity potential, are selected as targets for future exploration.

1. Introduction

Prediction of ore body location is critical for mineral exploration [1,2,3,4], which is generally achieved by analysing the associations of geological factors, geophysical and geochemical fields with mineralization. Such associations are commonly represented as 2D maps, and quantitatively analysed in the platform of 2D Geographic Information System (GIS), which have definitely facilitated the understanding of mineralization systems and predictive exploration of mineral resources [5,6,7]. However, such 2D studies are insufficient to present and analyse complicated mineralization systems. Unfortunately, most ore deposits, especially of hydrothermal metallic ore deposits, are complexly formed and distributed in 3D space. Therefore, the 3D models and 3D GIS have become more and more important for researching complex geological systems and ore deposits [8,9,10,11].

The 3D geological modelling, which is a class of computer techniques for geological characterization [12,13], has been developed as a robust and useful tool with the advancement in computer technologies, geomathematics and computational 3D graphics. It is widely used for presenting complex geological systems and managing geological information. Generally, these various classes of information have played important roles in understanding geological systems and exploring mineral resources [14,15,16,17,18,19,20,21]. The techniques for 3D geological modelling mainly include two groups, surface-based and volume-based modelling. The surface-based modelling is the common method for representing the surface of the complex geological factors in 3D. Several algorithms have been developed for interpolation, surface construction and rendering in order that the models can approximate the shape of real geological factors as closely as possible by using limited sampling data [22]. These commonly used algorithms include Delaunay triangulation for constructing Triangulated Irregular Network (TIN) model [23,24,25], Discrete Smooth Interpolation (DSI) for surface smoothing [26,27] and Marching Cubes for isosurface extraction [28]. The specialized commercial geological modelling software have been developed by adopting the above algorithms. Benefiting from these software tools, it has become convenient to model the complex and irregular geological factors in 3D virtual space. Voxel-based modelling is commonly used to represent the geological objects with continuous attribute values, such as geophysical and geochemical fields. Inverse Distance Weighted (IDW) [29,30] and Kriging [31,32] have been approved methods in interpolation which is important for constructing field models. Voxelization is a key approach to convert the surface model into a block model with some methods, such as Flood-fill and Octree-based Divisive Algorithm [8,9].

Spatial analysis that is used for inquiring spatial information and revealing the characteristic of geological objects, is widely applied in evaluating and predicting mineral resources by combining with 3D GIS. Nowadays more and more quantitative mathematical methods have been developed to appraise the association of geological factors, geochemistry field and geophysics field with mineralization, and carry out 3D prospectivity mapping [33,34,35,36,37,38,39,40]. These methods are divided into two classes, knowledge-driven methods and data-driven methods [41,42]. The common knowledge-driven methods include Fuzzy logic [43], Analytical hierarchy process [6], Index Overlay [44] and Inference Networks and Decision tress in expert systems [45]. The common data-driven methods include Weights of Evidence (WofE) [46,47], Logic Regression [48] and Neural Networks [49]. The WofE is a statistical method based on Bayes theorem of conditional probability, and it has been extensively and successfully applied in prediction of mineral potential in 2D by integrating relevant and reliable information that are usually from geological, geophysical and geochemical surveys [47,50,51,52,53,54]. It has been successfully applied in 3D quantitative spatial analysis and prospectivity mapping in recent years [35,39,40].

The Anqing orefield is so attractive because it hosts the Anqing Cu–Fe skarn deposit that is the largest Cu–Fe deposit in the Yangtze River metallogenic belt. The skarn ore deposits are an important type of metallic deposits. The direct contact zone of the intermediate-acid intrusion with the carbonate rocks is definitely the key factor for controlling skarn orebodies [55,56,57]. Since the first drill intersected the high-grade Cu–Fe skarn in 1960 in this ore field, it has been explored intensively, resulting in the Cu–Fe deposits discovered extensively and irregularly in the field. These deposits and their related felsic intrusion, the Yueshan intrusion, have been studied by many researchers [58,59,60,61,62,63,64,65,66,67]. Most of these studies focus on the geological and geochemical characteristics, genesis and dating of ore-formation and magmatism for understanding mineralization processes. Only a few studies focus on the architecture and dynamics of the mineralization system for directly facilitating ore prediction and exploration. Most of orebodies in this ore field are located in the contact zones of the Yueshan intrusion which is complicated in shape and highly various in occurrence. These ore-favourable zones were formed when the crust was subjected to dilation deformation, which is a direct result of the coupled mechano–thermo–hydrological processes during syn-tectonic cooling of the intrusion [3,22,68,69]. However, the credible prediction of orebodies at depth is the toughest challenge for mineral exploration in this field. For increasing ore reserves through deep exploration, the Tongling Nonferrous Metals Group Co., Ltd. had conducted a Controlled Source Audiofrequency Magnetotelluric (CSAMT) survey of 26.6 km² on a grid spacing of 200 m × 40 m and a high precision magnetic survey of 34.2 km² on grid spacing of 100 m × 40 m. Targeting at the anomalies from these geophysical surveying, 12 drills with total accumulative depth of 16,725.39 m have been completed. Unfortunately, no orebody has been discovered by these works. It suggests that the CSAMT and high precise magnetic surveys are not effective enough to discover orebodies at depth by traditional study, although the Anqing deposit was initially discovered by checking magnetic anomalies. Nowadays it is necessary to comprehensively appraise such a puzzling situation for facilitating predictive discovery of orebodies at depth. In this paper, we apply the methods of 3D modelling and spatial statistics to analyse the spatial association of the Yueshan intrusion, resistivity field and volume strain field with orebodies in 3D and, finally, to achieve an innovative prediction of orebodies at depth.

2. Geological Setting and Ore Deposits

2.1. Geological Setting

The Anqing orefield, located in the famous Tongling-Anqing Cu skarn district in eastern China, is a principal Cu–Fe ore field in the central segment of the Yangtze River metallogenic belt [3,22,46,61]. There are five groups of major structures in this orefield (Figure 1), including NW-, approximately N–S- and NE-trending folds and faults, and approximately E–W-, NNW- and NNE-trending faults. The NW- and NS-trending folds might have been formed immediately after the early Triassic, while the NE-trending folds might have been formed immediately after the middle Jurassic. The approximately E–W-trending faults consist mainly of normal faults and cut through the NE-trending folds and the Yueshan intrusion; the NNW- and NNE-trending faults are certainly the latest, because they cut through all other structures and orebodies. The approximately E–W-trending normal fault is parallel to and immediately nearby the orebodies [3,22].

The key ore-controlling factor in this ore field is the Yueshan intrusion (shown in Figure 1), the largest copper mineralized intrusion with an outcropped area of about 11 km² in the district. It is composed mostly of diorite with minor quartz diorite. The U–Pb zircon isotopic age of the Yueshan intrusion is 138.7 Ma, or Mid Early Cretaceous, when the crust in this region was turned from compression into extension and tectonic regime [59,62,64].

The sedimentary rocks of Paleozoic through to Mesozoic are distributed around the Yueshan intrusion in the ore field. On the surface, the direct walls of the intrusion are of the Tertiary. The most favourable wall rocks for hosting Cu- and Cu–Fe-skarn deposits, are littoral to neritic carbonates interbedded with bathyal facies beds that alternated with marine-continental clastic rocks, especially where they occur in contact with the diorite intrusion (Figure 1). These carbonate rocks are distributed in the Yueshan Formation of Middle Triassic (T₂y, dolomite and limestone breccia) and the Nanlinghu Formation of Lower Triassic (T₁n, limestone) [22].

2.2. Ore Deposits

There are hundreds of orebodies discovered in the Anqing orefield. These orebodies are mainly of Cu and Cu–Fe skarns, minor of Fe skarns. Their total reserves are more than 48.8 Mt copper ores at average grade of 1.3% Cu and 100 Mt iron ores at average grade of 48% Fe. Most of these orebodies occur in the contact zone between the Yueshan intrusion and its wall rocks (Figure 1), marble and dolomite marble of the Low to Mid Triassic. Only minor orebodies are within the intrusion, but close to the contact zone and related to the carbonate xenoliths in the diorite intrusion. The orebodies around the Yueshan intrusion are locally distributed. More than 99% of total Cu and Fe reserves proven in the ore field are distributed in two major deposits, the Anqing Cu–Fe deposit and the Zhuchong Fe–Cu deposit.

The Anqing Cu–Fe deposit with copper-iron ore reserves of 41 Mt Cu @ 1.3% and Fe @ 45% is located in the south contact zone between the diorite of the Yueshan intrusion and the marble of the T₁n and dolomite marble of T₂y. The orebodies are mainly composed of Cu sulfide-bearing massive magnetite and unmineralized skarns, with minor disseminated-copper sulphides in diorite. The lines C–C′, D–D′ and E–E′ in Figure 2 show that the orebodies are irregularly shaped, E–W-trending and steeply dipping toward the south. The south contact zones are not completely mineralized, which is related to its topographic variation. From east to west, the “step-shaped” contact zone becomes a “tongue shaped” contact zone. In the former, the contact zone of the intrusion with the carbonate hanging wall changes from a steep-dip into a gentle-dip; while in the latter, the marble is surrounded by the diorite. The orebodies are only located in the steep segments of the “step-shaped” contact zone and the tips of the “tongue-shaped” contact zone [3,22].

The Zhuchong Fe-Cu deposit, with iron ores of 51.62 Mt Fe @ 50.1% and copper ores of 4.4 Mt Cu @ 1.15%, is located in the north contact zone between the diorite of the Yueshan intrusion and the dolomite marble of the Yueshan formation. The major orebodies composed of Cu sulphide-bearing massive magnetite occur as irregular lenses scattered in concaves of the intrusion’s contact zone (shown in Line A–A′ and B-B′ of Figure 2). Minor orebodies composed of disseminated-copper sulphides in diorite and skarn are distributed as irregular veins within the intrusion close to the contact zone. The orebodies are E–W-trending and gently northward dipping.

3. Methods and Algorithms

3.1. TIN and DSI Methods for Modeling Geological Factors

The surface-based geometric modelling is the most commonly used for representing the shape of geological factors, such as topography, geological formations, intrusions, alteration zones, faults and orebodies [70,71,72]. Due to its efficiency in data storage and ability to accommodate irregularly spaced boundary points, the TIN method is increasingly popular for surface-based modelling. TIN modelling is a representation of a continuous surface that consisted entirely of triangular facets. The non-overlapping irregular triangular meshes are created from the scattered data set without repeating points according to certain rules. The Delaunay triangulation is the most common tessellation algorithm [25,73]. The empty circumcircle and Max–Min Angle criterion is the basic rule of the Delaunay algorithm. A TIN surface is defined as a set of triangles, which are continuous and not overlapped. It maximizes the minimum angle of all the angles of the triangles in the triangulation, and they tend to avoid sliver triangles and hold uniqueness of triangular mesh [73,74].

For optimizing the TIN surface models, Mallet [26,27] developed the DSI algorithm. Through node interpolation, the optimized triangular meshes are continuous, non-overlapping, irregular, and have acute triangles and approximate sides. The DSI optimization is dependent on the topotaxy of grid nodes, not limited by space dimensions. If a set of known grid nodes satisfied a certain constraint conditions, mutual contact nets can be constructed between discrete data points, meanwhile the value of an unknown node will be obtained by solving a linear equation.

3.2. Interpolation Algorithms for 3D Block Modeling

The 3D block modelling has been widely applied for presenting field models, such as temperature, resistivity and strain. In a 3D block model, every block corresponds to the actual position and attribute information, which can be easily queried. The data sets for block modelling are generated by interpolation from the discrete data. The 3D interpolation methods, such as Kriging and IDW, are the mathematic core for 3D block modelling.

The Kriging algorithm is a geostatistical interpolation technique that considers both the distance and the degree of variation between known data points when estimating values in unknown areas [31]. A kriged estimate is a weighted linear combination of the known sample values around the point to be estimated. The cornerstone of the Kriging algorithm is the semivariogram function, or the variogram as it is generally referred to by geostatisticians. In the interpolation process, the weights not only depend on the distance between sampling and unknown points, but also on whole distribution of all the sampling points.

The IDW is to estimate the assigned values of unknown points by a weighted average of the values available at the known points and it resorts to the inverse of the distance to each known point when assigning weights [29]. Assuming that in the neighbourhood of the point, P(x, y, z), which needs an interpolated value A, there are discrete datasets Q_i (x_i, y_i, z_i), i = 1, 2, …, n. The value A_p is obtained by interpolation of a weighted average for Q_i. The weight is related with the distance between interpolated point P and interpolating point Q_i. Generally, A_p is defined as:

$$A_P=\{\begin{array}{ll}\frac{{{\displaystyle \sum }}_{i=1}^n\frac{Q_i}{d_i^k}}{{{\displaystyle \sum }}_{i=1}^n\frac{1}{d_i^k}}& \mathrm{if} \left(x,y,z\right)\ne \left(x_i,y_i,z_i\right) \mathrm{for} \mathrm{all} i\\ Q_i& \mathrm{if} \left(x,y,z\right)=\left(x_i,y_i,z_i\right) \mathrm{for} \mathrm{all} i\end{array}$$

(1)

where $d_i=\sqrt{{\left(x-x_i\right)}^2+{\left(y-y_i\right)}^2+{\left(z-z_i\right)}^2}$, k, power exponent for estimation, is a positive real number.

3.3. The WofE Method for 3D Spatial Analysis

The WofE which is based on a loglinear form of Bayesian rule is a discrete multivariate geostatistical method. It is widely used in prediction and evaluation of mineral resources by associating the mineral occurrence with a series of geological data [46,47]. In the WofE model, the mineralization occurrences (orebodies) are used for training, and the geological (or exploration) data are served as evidential factors by calculating the weights and contrasts of all the evidential factors, which are used as binary maps with mineralization occurrences, the appropriate exploration criterion that are used as binary predictive maps, are combined to calculate the posterior probability (P_post). For each binary map, “1” means the evidence is present, “0” for absent, and “null” for unknown or missing. These weights and contrasts are helpful to express the spatial association of evidential factors with mineralization occurrences, and the P_post is useful for plotting predictive maps and locating orebodies.

In the process of 3D association analysis using WofE, the study area which contains the mineralization occurrences (orebodies) and all evidential factors is subdivided into N(T) identical blocks in 3D space. The occurrences and evidential factors are used in voxel form. The training blocks that contain occurrences are generated from the surface-based 3D models of orebodies by voxelization. The voxelization is the approach that generating volume datasets by discretizing continuous objects on a regular grid of voxels in 3D Euclidean space. It is usually realized through the method of Flood-fill Algorithm, Octree-based Divisive Algorithm and distance transformation [5,9,74,75,76,77]. If there are D training blocks, the number of blocks that don’t contain occurrences are $\mathrm{N}\left(\overline{D}\right)=\mathrm{N}\left(T\right)-\mathrm{N}\left(D\right)$. The prior probability that any given block will contain training blocks is $\mathrm{P}\left\{D\right\}=\mathrm{N}\left(D\right)/\mathrm{N}\left(T\right)$, and expressed as odds by:

$$\mathrm{O}\left\{D\right\}=\frac{\mathrm{P}\left\{D\right\}}{1-\mathrm{P}\left\{D\right\}}=\frac{\mathrm{N}\left(D\right)}{\mathrm{N}\left(T\right)-\mathrm{N}\left(D\right)}$$

(2)

For the evidential factor (binary map) B_j, the number of blocks where the evidence occurred is N(B_j), and the number of blocks where the evidence did not occur is $\mathrm{N}\left(\overline{B_j}\right)=\mathrm{N}\left(T\right)-\mathrm{N}\left(B_j\right)$. According to Bayesian rule, the conditional probability of choosing a block with mineralization occurring, where the evidence B_j occurred is:

$$\mathrm{P}\left\{D|B_j\right\}=\frac{\mathrm{P}\left\{D\cap B_j\right\}}{\mathrm{P}\left\{B_j\right\}}=\frac{\mathrm{N}\left(D\right)\cap \mathrm{N}\left(B_j\right)}{\mathrm{N}\left(B_j\right)}$$

(3)

The weights for the evidential factor B_j can be defined as:

$$\begin{array}{ll}{W}_j^{+}=\ln \frac{\mathrm{P}\left\{B_j|D\right\}}{\mathrm{P}\left\{B_j|\overline{D}\right\}}& \mathrm{for}\mathrm{evidence} \mathrm{of} B_j\mathrm{present}\\ W_j^{-}=\ln \frac{\mathrm{P}\left\{\overline{B_j}|D\right\}}{\mathrm{P}\left\{\overline{B_j}|\overline{D}\right\}}& \mathrm{for}\mathrm{evidence} \mathrm{of} B_j\mathrm{absent}\end{array}$$

(4)

where $W_j^{+}$ is the positive weight and $W_j^{-}$ is the negative weight, and they have opposite meanings in terms of identification of spatial association. It indicates a positive spatial association between the evidence B_j and mineralization when $W_j^{+}>0$ or $W_j^{-}<0$, conversely it indicates a negative association when $W_j^{+}<0$ or $W_j^{-}>0$.

More details on the derivation of these formulae are given in Bonham-Carter [46].

The contrast, C_j, is used to express the strength of the spatial association between evidential factors and mineralization. It is defined as:

$$C_j=W_j^{+}-W_j^{-}$$

(5)

In this equation, when $C_j>0$, the spatial association is positive; and $C_j<0$ for the negative spatial association. In the special case of $C_j=0$, it means a lacking relationship. The higher C_j values indicate more significant spatial association between the mineralization and the given evidential factor B_j.

The P_post is calculated by combining n binary predictive maps which are conditionally independent with respect to the mineralization occurrences [46]. The P_post is expressed by log likelihood ratio:

$$\ln \mathrm{O}\left\{D|\left(B_1^k\cap B_2^k\cap \dots \cap B_n^k\right)\right\}=\ln \mathrm{O}\left\{D\right\}+{\displaystyle \sum }_1^n{W}_j^k$$

(6)

where k is the status of evidence factor B_j in given block, and

$$W_j^k=\{\begin{array}{ll}{W}_j^{+}& \mathrm{for} \mathrm{evidence} \mathrm{of} B_j\mathrm{present}\\ W_j^{-}& \mathrm{for} \mathrm{evidence} \mathrm{of} B_j \mathrm{absent}\\ 0& \mathrm{for} \mathrm{evidence} \mathrm{of} B_j \mathrm{missing}\end{array}$$

(7)

Supposing that

$$f=\ln \mathrm{O}\left\{D\right\}+{\displaystyle \sum }_1^n{W}_j^k$$

(8)

the P_post can be calculated by using

$$P_{\mathrm{post}}=\mathrm{P}\left\{D|\left(B_1^k\cap B_2^k\cap \dots \cap B_n^k\right)\right\}=\frac{e^f}{1+e^f}$$

(9)

The threshold values obtained from capture-efficiency curves [78], are used to effectively separate potential areas. The P_post map, which created by the combination of binary predictive maps and P_post value, is used to demarcate the exploration targets. The higher P_post values indicate the high probability that the ore body is distributed.

4. Spatial Association of Intrusion, Resistivity and Volume Strain Fields with Orebodies

4.1. 3D Spatial Association of the Yueshan Intrusion with Orebodies

Liu et al. [22] had constructed 3D models of the Yueshan intrusion and its related orebodies. Since then, 263 drill holes with a total length of 176,058 m and about 5000 m of underground tunnels have been finished. These works have provided abundant information of the underground geology and led to the discovery of many orebodies. By integrating all the newly-acquired information from these works on the Micromine 11.0 platform (MICROMINE Head Office, Perth, Australia), we reconstructed the surface-based models of the Yueshan intrusion, ore-favourable carbonate formations, T₁n and T₂y, as well as all orebodies (Figure 3a).

Compared to the previous models of Liu et al. [22], these new 3D models reveal the geometric features and spatial association of the intrusion and orebodies:

(1)

The 3D Yueshan intrusion exhibits extreme variation in attitude and topography of its contact zone, and such variations have strong constraints on uneven localization of orebodies around the intrusion. Both the south contact zone and north contact zone host more than 99% ores in the ore field. Their common feature that is nearly E–W-trending is completely different from both the west and the east contact zone that are nearly S–N-trending (Figure 3a). It suggests that the E–W-trending contact zone is favourable for skarn mineralization.

(2)

The south contact zone hosts much more Cu reserves than the north one. It also displays distinct differences from the north one in occurrence and topography. The south contact zone has an extremely irregular surface with a wide range of dips, from northward (inward) 40°–60° to southward (outward) 20°–45°, while the north one has gentle waved surface with a stable northward (outward) tip about 25°–40° (Figure 3b). These suggest that the rapid change of occurrence might have made a difference for Cu mineralization in contact zone of the Yueshan intrusion. Particularly, the south contact zoon is more favourable for formation of orebodies that the north.

(3)

Almost all of orebodies are covered by the strata, T₁n and T₂y. No orebody is distributed in the east branch of the Yueshan intrusion, where no carbonate rock is present. It illustrates that the mineralization is controlled by the presence of carbonate rocks in direct contact with diorite rocks.

(4)

The orebodies are unevenly localized in the vicinity of the intrusion’s contact zone, especially around the concaves (Figure 4a,b). Remarkably, such uneven localization of orebodies is related to the curvatures [79] of the contact surface. By using the minimum principal curvature (K_min) to describe the topography of the intrusion’s contact zone (Figure 5a), it is evident that most orebodies are localized around the contact zone with |K_min| > 0.0025 (Figure 5b).

The study area for Yueshan intrusion and its surrounding areas are divided into a total of 278,0526 blocks by using a 40 m × 40 m × 40 m grid size, and the surface model of orebodies is voxelized into 3186 blocks with the same dimensions in the study area. We selected the concaves as abnormal zones (Figure 5b) from the contact zone of the Yueshan intrusion with |K_min| > 0.0025. The 3D buffers that used to identify the influence range and degree of some point, line and surface geological objects [80,81,82], are constructed around these zones by calculating Euclidean distances between individual blocks in the study area. To analyze the spatial association of the Yueshan intrusion with orebodies, a series of WofE calculations were processed for different buffers with block dimensions. The calculation results (

Table 1. Weights for 3D buffers around the chosen abnormal zones of Yueshan intrusion contact.

Buffer Distance (m)	Training Blocks	3D Buffers Blocks	Positive Weight	Negative Weight	Contrast
0–40	1670	13,257	4.833	−0.739	5.572
40–80	635	21,356	3.285	−0.215	3.5
80–120	419	39,987	2.223	−0.127	2.349
120–160	142	35,011	1.267	−0.033	1.3
160–200	115	48,498	0.729	−0.019	0.748
200–240	53	52,795	−0.132	0.002	−0.135
240–280	24	42,926	−0.718	0.008	−0.726
280–320	23	50,185	−0.917	0.011	−0.928
320–360	22	53,040	−1.017	0.012	−1.029
360–400	23	47,611	−0.864	0.01	−0.874
400–440	23	48,613	−0.885	0.01	−0.896
440–480	17	45,318	−1.117	0.011	−1.128
480–520	16	46,371	−1.201	0.012	−1.213
520–560	2	38,500	−3.095	0.013	−3.108
560–600	2	34,698	−2.991	0.012	−3.003

, Figure 6a) show that, all the training blocks are located in the buffers within the distance of 600 m which occupies 22.23% of the total blocks. From 40 m to 600 m, the positive weights decline approximately to −3.095 from 4.833, and the contrasts decline to −3.108 from 5.572. For the buffers within 200 m from the contact zone, they occupy 5.69% of the total blocks, but include 93.57% of training blocks; their positive contrasts suggest that the concaves zone we selected from the Yueshan intrusion contact must have positive constraints on the localization of orebodies (Figure 7). This indicates that the enrichment of the orebody is closely related to distance from contact zone of the Yueshan intrusion, much closer to the intrusion’s contact, and contributing more to mineralization.

4.2. 3D Resistivity Field and Its Association with Orebodies

Geophysical prospecting had played an important role in the discovery of the Anqing deposit. Although the CSAMT surveying have been carried out in an area of 26.4 km², few orebodies have been discovered by this surveying because of a lack of quantitative analysis between the resistivity and mineralization in 3D space. Through using the IDW method to interpolate the resistivity data (ρ) from the CSAMT surveying, we constructed the 3D resistivity field with the block dimensions of 40 m × 40 m × 40 m (Figure 8a), and the space in the study area without surveying data were excluded from the spatial analysis. The Cu–Fe skarn orebodies that are especially composed of Cu-sulphide bearing massive magnetite ores, are generally recognized as the low resistivity geological factors. However, the 3D modelling results show that almost no orebodies occur in the space of ρ <500 Ω·m (Figure 8b). There are few orebodies in the space of ρ >7000 Ω·m (Figure 8c). A majority of orebodies are localized in the spaces with ρ between 500 and 2500 Ω·m (Figure 8d), where they are much larger than the training blocks within them.

We use the WofE method to analyse spatial association between the space of different resistivity and orebodies. The calculation results (

Table 2. Weights for different spaces of resistivity field at 500 Ω·m intervals.

Resistivity (Ω·m)	Training Blocks	Resistivity Field Blocks	Positive Weight	Negative Weight	Contrast
0–500	7	108,326	−3.757	0.160	−3.917
500–1000	261	131,683	−0.332	0.061	−0.393
1000–1500	569	88,472	0.850	−0.207	1.056
1500–2000	428	62,106	0.919	−0.153	1.072
2000–2500	114	46,792	−0.125	0.008	−0.133
2500–3000	96	37,656	−0.080	0.004	−0.084
3000–3500	102	31,350	0.165	−0.008	0.173
3500–4000	84	26,459	0.140	−0.006	0.146
4000–4500	72	24,442	0.065	−0.002	0.068
4500–5000	59	20,238	0.055	−0.002	0.056
5000–5500	56	16,482	0.208	−0.005	0.214
5500–6000	31	13,868	−0.212	0.004	−0.215
6000–6500	32	12,117	−0.044	0.001	−0.045
6500–7000	19	10,945	−0.465	0.006	−0.471
7000–7500	14	10,260	−0.706	0.007	−0.713
7500–8000	8	9005	−1.136	0.009	−1.144
greater than 8000	34	69,258	−1.729	0.084	−1.813

) shows that, for the spaces of ρ <500 Ω·m, there are only seven training blocks, and the positive weight and contrast are all negative, indicating that it is unlikely to host orebodies there. In the spaces of 1000 Ω·m ≤ ρ ≤ 2000 Ω·m, the contrasts are positive and much greater than other spaces (Table 2, Figure 6b), the blocks occupy 20.93% of the resistivity field and contain 50.20% of training blocks. It indicates that these spaces are more closely associated with mineralization in the Anqing orefield than other spaces, however, no tendency can be discovered easily in such spatial association. These WofE study results demonstrate that the orebodies cannot be detected as lowest resistivity geological factors by the CSAMT surveying and the spatial association between orebodies and the resistivity field is weak, which accounts for the reason why the targets delineated mainly by the CSAMT have no-identified orebodies.

4.3. 3D Volume Strain Field and Its Association with Orebodies

By using computation geodynamics modelling, Liu et al. [22] studied the coupled mechano–thermo–hydrological dynamics during the Yueshan intrusion’s cooling process-modelling. The computational modelling results show that the deposition of ores or formation of ore-hosting space is related to the volumetric strain increment (vsi) of which positive value present with dilation deformation [22]. In order to analyse qualitatively the spatial association of volumetric strain with mineralization, we constructed a 3D volume strain field (Figure 9) by using vsi and calculate weights of spaces for different interval vsi with training blocks.

The calculation results (

Table 3. Weights for different spaces of volume strain field at 0.2% intervals.

Vsi (%)	Training Blocks	Volume Strain Field Blocks	Positive Weight	Negative Weight	Contrast
less than 0	1027	1,764,862	−0.678	0.618	−1.296
0–0.2	619	490,349	0.097	−0.022	0.119
0.2–0.4	540	267,098	0.569	−0.085	0.654
0.4–0.6	435	156,733	0.886	−0.089	0.975
0.6–0.8	279	58,388	1.432	−0.071	1.502
greater than 0.8	286	43,096	1.762	−0.078	1.838

, Figure 6c) show that when vsi < 0, the positive weight is less than zero and the negative weight is greater than zero. The blocks in this space occupy almost the entire study area, but only a small part of orebodies is located there. When vsi > 0, the positive weight is greater than zero and the negative weight is less than zero. The contrasts show a tendency to rise along with the increase of vsi. These illustrate quantitatively that the mineralization is closely related to the dilation zones, and the spaces with higher vsi facilitate the formation of orebodies. When vsi ≥ 0.6%, the strain field blocks occupy 3.65% of total blocks containing 19.74% of training blocks, and the contrasts are both greater than 1, indicating that there is a better correlation between volume strain and mineralization in these intervals. By comparing the blocks of vsi ≥ 0.6% with the Yueshan intrusion, we find that the shapes of higher dilation zones are generally consistent with the contact zone of the intrusion, whereas the spatial association of the volume strain field with orebodies is not as strong as the contact of intrusion.

5. Quantitative Integration of Ore-Related Information and Ore Prediction

5.1. Binary Predictive Maps of Exploration Criteria

Based on the analysis above, we found that these evidential factors have positive correlation with mineralization: (1) the 3D buffers around the chosen areas that the distance within 200 m to the Yueshan intrusion; (2) the carbonate rocks that are distributed in the strata, T₁n and T₂y; (3) resistivity field with resistivity between 1000 and 2000 Ω·m; and (4) volume strain field of which vsi is greater than 0.6%. Although the correlation of carbonate rocks, resistivity and volumetric strain with orebodies is not as strong as the contact of the Yueshan intrusion, the exploration criteria from the evidential factors, B, C, V and R (

Table 4. Exploration criteria determined for prospectivity mapping from evidential factors.

Evidential Factors	Exploration Criterion	Favourable Range for BPM 1	Conditional Probability (%)	Positive Weight	Negative Weight	Contrast
Contact zone of the Yueshan intrusion	B: 3D buffers (around chosen contact zones)	0–200 (m)	1.884	2.818	−2.676	5.494
Stratigraphy	C: carbonate	T1n and T2y	0.896	2.064	−0.567	2.631
3D volume strain field	V: volumetric strain increment	greater than 0.6%	0.655	1.584	−0.158	1.742
3D resistivity field	R: resistivity	1000–2000 (Ω·m)	0.556	0.869	−0.453	1.322

¹ BPM: binary predictive map.

, Figure 10), are very necessary for calculating posterior probability and selecting exploration targets in this ore field. The assumption of conditional independence is very important for every possible pair of the binary predictive maps [45], it can be examined by applying all the binary maps in mineral occurrence through chi-square testing [81]. Generally, the conditional independence assumption for 3D weights-of-evidence modelling is violated [35,82]. These exploration criteria are converted into binary predictive maps, for which the presence and absence (or missing) of an evidence feature are respectively denoted by B⁺, B⁻, C⁺, C⁻, V⁺, V⁻, R⁺, R⁻ and R⁰. The maximum of conditional probability for co-occurrences of mineral occurrence and every evidential factor is 1.884%, and the prior probability, P{D}, is 0.115%.

5.2. Prospectivity Mapping and Ore Prediction

The binary predictive maps mentioned above are combined to calculate P_post by the Formulae (6)–(9) and create the prediction map, the calculation results are shown in

Table 5. Combination of binary predictive maps with different posterior probability.

CofBPMs 1	Posterior Probability (%)	Training Blocks	Total Blocks
B+C+V+R+	64.003	60	355
B+C+V+R0	42.476	71	3604
B+C+V+R−	31.720	190	4925
B+C+V−R+	23.722	180	2215
B+C+V−R0	11.438	350	18,077
B+C−V+R+	11.346	41	296
B+C+V−R−	7.515	626	11,334
B+C−V+R0	5.047	125	10,568
B+C−V+R−	3.236	72	5034
B+C−V−R+	2.190	210	8168
B+C−V−R0	0.921	511	65,414
B−C+V+R+	0.726	0	2456
B+C−V−R−	0.581	543	28,119
B−C+V+R0	0.303	6	7490
B−C+V+R−	0.191	0	5551
B−C+V−R+	0.128	0	9270
B−C+V−R0	0.053	3	70,145
B−C−V+R+	0.053	0	7794
B−C+V−R−	0.033	0	30,461
B−C−V+R0	0.022	0	27,843
B−C−V+R−	0.014	0	25,708
B−C−V−R+	0.009	0	120,053
B−C−V−R0	0.004	134	1,857,926
B−C−V−R−	0.002	64	457,720

¹ CofBPMs: combination of binary predictive maps.

. We plot capture-efficiency curves by using the cumulative proportion of total blocks against its corresponding training blocks from high to low posterior probability (Figure 11a), as well as posterior probability curves against cumulative proportion of total blocks (Figure 11b). Two thresholds are defined through comparing these two curves. The threshold 1 corresponds with the P_post of 2.190%, which is greater than the maximum of conditional probability, 1.884%, and the threshold 2 yields a P_post of 0.582%, which is greater than the minimum of conditional probability, 0.556%. In addition, they are all greater than the prior probability, 0.115%, suggest that the study area can be demarcated into three classes of potential blocks for ore prospecting according to the P_post values, (1) for high potential, P_post ≥2.190%; (2) for medial potential, 2.190% > P_post ≥0.582%; and (3) for low potential, P_post <0.582% (Table 5).

The Figure 12 shows the spatial association of the Yueshan intrusion and orebodies with combination of binary predictive maps for high and medial potential blocks. The high potential blocks, B⁺C⁺V⁺R⁺, B⁺C⁺V⁺R⁰, B⁺C⁺V⁺R⁻, B⁺C⁺V⁻R⁺, B⁺C⁺V⁻R⁰, B⁺C⁻V⁺R⁺, B⁺C⁺V⁻R⁻, B⁺C⁻V⁺R⁰, B⁺C⁻V⁺R⁻ and B⁺C⁻V⁻R⁺, occupy 2.32% of total blocks and contain 60.42% of known orebodies (Figure 12a), and for medial potential blocks, B⁺C⁻V⁻R⁰, B⁻C⁺V⁺R⁺ and B⁺C⁻V⁻R⁻, occupy 3.45% of total blocks and contain 33.08% of the training blocks (Figure 12b). The P_post maps (Figure 13a) and its east, north and vertical direction sections (Figure 13b–d) obviously show that the spaces with high P_post are well consistent with the contact of the Yueshan intrusion, demonstrating that the contact of the Yueshan intrusion has strong constraints on localization of orebodies in the Anqing orefield, especially where concaves occur.

We can select exploration targets with different priority for future mineral exploration by P_post values. The favourable targets can be identified in high and medial potential blocks where no orebody was discovered and had never been drilled. Figure 14 shows five (A–E) targets selected from high potential blocks. The accurate location of any target blocks can be easily queried in the 3D block model.

6. Conclusions

This paper analysed the spatial association of the Yueshan intrusion, carbonate rocks, resistivity and volume strain field with orebodies by using 3D WofE modelling. The Yueshan intrusion and its shape of contact zone, as well as the carbonate rocks in T₁n and T₂y, are closely related to mineralization. The spaces with the distance ≤200 m to the contact surface, where |K_min| is greater than 0.0025, have strong positive constraints on localization of orebodies. The orebodies cannot be detected as lowest resistivity geological factors as generally recognized in this ore field. There is no good correlation between resistivity and mineralization, and it is hard to target orebodies individually by resistivity with a certain interval from CSAMT surveying. For the volume strain field, a high value of volume strain indicates the stronger dilation deformation, and the spaces with high positive vsi are favourable for the localization of orebodies.

The prospectivity map was plotted by integrating the binary predictive maps that were created from the contact zone of the Yueshan intrusion, strata in T₁n and T₂y, resistivity and volume strain field. The study area was demarcated into high, medial and low potential blocks by P_post value. Five targets have been preferentially determined from high potential blocks for future mineral exploration in the Anqing orefield.