A Remote Sensing ‐ Based Application of Bayesian Networks for Epithermal Gold Potential Mapping in Ahar ‐ Arasbaran Area, NW Iran

: Mapping hydrothermal alteration minerals using multispectral remote sensing satellite imagery provides vital information for the exploration of porphyry and epithermal ore mineralizations. The Ahar ‐ Arasbaran region, NW Iran, contains a variety of porphyry, skarn and epithermal ore deposits. Gold mineralization occurs in the form of epithermal veins and veinlets, which is associated with hydrothermal alteration zones. Thus, the identification of hydrothermal alteration zones is one of the key indicators for targeting new prospective zones of epithermal gold mineralization in the Ahar ‐ Arasbaran region. In this study, Landsat Enhanced Thematic Mapper+ (Landsat ‐ 7 ETM+), Landsat ‐ 8 and Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) multispectral remote sensing datasets were processed to detect hydrothermal alteration zones associated with epithermal gold mineralization in the Ahar ‐ Arasbaran region. Band ratio techniques and principal component analysis (PCA) were applied on Landsat ‐ 7 ETM+ and Landsat ‐ 8 data to map hydrothermal alteration zones. Advanced argillic, argillic ‐ phyllic, propylitic and hydrous silica alteration zones were detected and discriminated by implementing band ratio, relative absorption band depth (RBD) and selective PCA to ASTER data. Subsequently, the Bayesian network classifier was used to synthesize the thematic layers of hydrothermal alteration zones. A mineral potential map was generated by the Bayesian network classifier, which shows several new prospective zones of epithermal gold mineralization in the Ahar ‐ Arasbaran region. Besides, comprehensive field surveying and laboratory analysis were conducted to verify the remote sensing results and mineral potential map produced by the Bayesian network classifier. A good rate of agreement with field and laboratory data is achieved for remote sensing results and consequential mineral potential map. It is recommended that the Bayesian network classifier can be broadly used as a valuable model for fusing multi ‐ sensor remote sensing results to generate mineral potential map for reconnaissance stages of epithermal gold exploration in the Ahar ‐ Arasbaran region and other analogous metallogenic provinces around the world. phyllic zone. 3, successfully indicates the known epithermal gold mineralizations and several new high prospective zones in the study area. The approach developed in this study is a cost ‐ effective technique that can be used for epithermal gold exploration in metallogenic provinces before costly geophysical and geochemical studies. Briefly, this study suggests that geostatistical techniques (e.g., Bayesian network model, Fuzzy model, Artificial Neural Network Model etc.) are valuable approaches to fuse thematic layers of the multi ‐ sensor imagery for generating the remote sensing ‐ based mineral potential map for metallogenic provinces. The mineral exploration community and mining companies can consider the remote sensing ‐ based mineral potential map as an economical and cost ‐ effective tool for mineral prospecting before pricey geophysical and geochemical surveys in the metallogenic provinces.

Obtaining information from multi-sensor remote sensing satellite data can produce relevant results for detailed mapping of hydrothermal alteration zones [12]. The integration of the multi-sensor remote sensing results using geostatistical techniques can quickly produce a mineral potential map, which indicates the high potential zones of hydrothermal ore mineralizations [40]. Mineral potential map of a region is generally realized as the predictive classification of each spatial unit contains a particular combination of spatially coincident predictor patterns as mineralized or barren zones [41,42]. A Bayesian network is a type of statistical model (probabilistic graphical model), which represents a set of variables and their conditional dependencies through a Directed Acyclic Graph (DAG) [41,43,44]. It predicts the likelihood that anyone of several possible known causes was the contributing factor [45]. Therefore, the Bayesian network is a suitable model for fusing thematic layers derived from multi-sensor remote sensing satellite data to generate a mineral potential map.
In this study, Landsat-7 ETM+, Landsat-8 and ASTER multispectral remote sensing datasets were used to identify hydrothermal alteration zones associated with epithermal gold mineralization and producing thematic layers, which were afterward synthesized in the Bayesian networks for mineral potential mapping in the Ahar-Arasbaran region, NW Iran ( Figure 1). This region is a well-endowed terrain hosting numerous known epithermal gold deposits, several porphyry and skarn Cu-Mo deposits, Fe skarn deposits, Cu-Au porphyry deposits and many other Cu-Mo-Au vein mineralizations [46][47][48][49][50]. The deposits are associated with extensive hydrothermal alteration mineral zones such as iron oxide/hydroxides, advanced argillic, argillic, phyllic and propylitic [48,51,52]. The Ahar-Arasbaran region has a high potential for exploring new prospective zones of epithermal gold and many other ore mineralizations. Pazand et al. [52] used ASTER satellite data for hydrothermal alteration mapping in the Ahar area, NW Iran. Some geo-referenced hydrothermal alteration maps were produced using RBD (relative absorption band depth), principal component analysis (PCA), minimum noise fraction (MNF) and matched filtering (MF) image processing techniques for reconnaissance stages of porphyry copper exploration in the Ahar area. Furthermore, Pazand and Hezarkhani [48] generated a favorability map for Cu porphyry mineralization using fuzzy modeling in the Ahar-Arasbaran zone, NW Iran. There is no comprehensive remote sensing research available for mapping hydrothermal alteration zones in the Ahar-Arasbaran region using multi-sensor satellite imagery at a regional scale. This study characterizes an extensive remote sensing analysis using Landsat-7 ETM+, Landsat-8 and ASTER datasets, detailed fieldwork and laboratory analysis for mineral potential mapping. Therefore, the primary purposes of the research are: (1) to map hydrothermal alteration mineral zones using Landsat-7 ETM+, Landsat-8 and ASTER datasets by implementing the band ratio, PCA, RBD and selective PCA image processing techniques; (2) to generate mineral potential map by fusing the alteration thematic layers using the Bayesian networks; and (3) to verify the high potential zones by checking the detailed global positioning system (GPS) surveying in the field and analyzing several microphotographs of hydrothermal alteration minerals and gold mineralization and X-ray diffraction (XRD) analysis of collected rock samples from alteration zones.

Geology of the Ahar-Arasbaran Region
The Ahar-Arasbaran region covers an area (approximately 5000 km 2 ), which is located between latitudes 38°07′N and 38°52′N and longitudes 46°15′E and 47°30′E ( Figure 1). This zone is a part of Lesser Caucasus metallogenic zone and corresponding to tectono-magmatism activity from Jurassic to Quaternary [46,47,53,54]. The volcano-plutonic belt of Arasbaran-Lesser Caucasus is a mountainous and uplifted region that trending NW-SE from Georgia (Republic of Azerbaijan) to the Talesh region (Iran) [50]. Magmatic rocks in the Ahar-Arasbaran region containing tholeiitic, calc-alkaline, high calcium calc-alkaline, shoshonitic, adakitic, alkaline sodic and potassic rocks, which are formed in a continental margin of a subduction zone (subduction to post-collision stages) [50]. Cretaceous units (limestone and shale), flysch deposits, Paleocene and Eocene volcanic rocks are also exposed in the study area ( Figure 1). Several intrusive bodies having different sizes are penetrated in the Eocene and Cretaceous volcanic-sedimentary rocks and caused folding, alteration and mineralization [49,51]. Structural trends of folds, faults, dykes and veins are mostly NW-SE, E-W and NE-SW, which show the main stresses that affected the study area [50,51].
The intrusion of the Oligo-Miocene batholiths into the Cretaceous to Eocene sedimentary and volcano-sedimentary deposits along with hydrothermal fluids is formed intensive alteration halos in the Eocene volcanic rocks [55]. The alteration zones such as argillic, silica and alunite are associated with Cu, Au, Mo, Ag, Pb and Zn mineralizations [49]. Moreover, several skarn zones are formed in the contact zone of intrusive masses with Cretaceous limestone [51]. A variety of ore mineralization zones were identified in the Ahar-Arasbaran region, including Fe, Cu, Pb-Zn, Cu-Au, Cu-Mo, Au-Ag, Fe-Au, which occurred in the form of sprains, veins, stokes and in relation to the skarn zones [49,51]. The gold mineralization in the study area is observed in the form of epithermal veins [55]. The Masjed Daghi (Siahrood) and AliJavad valley (Anjerd) are considered to be Au-Cu porphyry deposits. The Sharaf-abad, Hize-jan, Nabi-jan, Zailig, Miveh-roud, Safi-Khanloo, Noqdouz, Anniqh and Khoyneh-roud are known as epithermal gold deposits in the study area [55].

Remote Sensing Data and Pre-Processing
The Landsat-7 ETM+, Landsat-8 and ASTER satellite remote sensing datasets were used in this study. Technical characteristics of Landsat-7 ETM+, Landsat-8 and ASTER remote sensing sensors are summarized in Table 1. A Landsat-7 ETM+ scene (Path/Raw: 168/33) covering the Ahar-Arasbaran region was acquired on 15 June 2001. A level 1T (terrain corrected) Landsat 8 scene (Path/Raw: 168/33) was also acquired on 10 June 2016 for the study area. Seven level 1B ASTER scenes covering the study area were acquired from 8 to 29 June 2002-2004. The data were obtained from the U.S. Geological Survey's Earth Resources Observation System (EROS) Data Center (EDC) (https://earthexploere.usgs.gov/ and https://glovis.usgs.gov). The scenes were cloud-free and have been already georeferenced to the UTM zone 38 North projection using the WGS-84 datum. For converting Landsat-7 ETM+ digital numbers to spectral radiance or exoatmospheric reflectance (reflectance above the atmosphere), the Landsat Calibration technique was adopted from Chander et al. [60]. This technique uses the published post-launch gain and offset values [61,62]. The mathematical details of the technical performance can be found in Chander et al. [60]. For Landsat 8 and ASTER datasets, Internal Average Relative Reflectance (IARR) was utilized. The IARR calibration method normalizes images to a scene average spectrum [61,63]. This is particularly effective for reducing imaging spectrometer data to relative reflectance in an area where no ground measurements exist and little is known about the scene [61,63]. It works best for arid areas with no vegetation. The IARR calibration is performed by calculating an average spectrum for the entire scene and using this as the reference spectrum. Apparent reflectance is calculated for each pixel of the image by dividing the reference spectrum into the spectrum for each pixel. The atmospheric correction was implemented to ASTER data after Crosstalk correction [64].
Moreover, the 15 m VNIR bands of ASTER were resampled to the 30 m SWIR bands using the cubic convolution technique. A masking procedure was applied to the remote sensing datasets for removing the effects of vegetation and Quaternary deposits. Normalized Difference Vegetation Index (NDVI) was calculated for the remote sensing datasets. As a result, a masking procedure was executed to the remote sensing datasets for eliminating the influences of sparse vegetation in the study area. For Quaternary deposits, we used geological map of the study area to identify the location of the Quaternary units, then a masking procedure was implemented to the remote sensing datasets. The ENVI (Environment for Visualizing Images, http://www.exelisvis.com) version 5.2 and ArcGIS version 10.3 (Esri, Redlands, CA, USA) software packages were employed for processing Landsat-7 ETM+, Landsat-8 and ASTER data.

Image Processing Techniques
The main objective of image processing techniques implemented in this analysis is to map hydrothermal alteration zones for generating thematic layers from multi-sensor remote sensing satellite datasets. Then, the thematic layers are fused using a Bayesian network model for producing a mineral potential map of the Ahar-Arasbaran region. Fieldwork and laboratory analysis are used to verify the results. A view of the methodological flowchart applied in this study is shown in Figure 2.

Band Ratio
The band ratio technique is one of the most applicable image processing techniques for mapping hydrothermal alteration minerals and zones such muscovite, jarosite, gossan, advanced argillic, argillic-phyllic, propylitic and hydrous silica-affected zones [23,39,65,66]. The digital number (DN) value of a band is partitioned by the DN value of other band, which highlights particular spectral features related to minerals or materials that planned to map [23]. Relative Absorption Band Depth (RBD) uses three-point ratio formulation for detecting typical absorption features related to a specific mineral or alteration zone [67]. For a specific absorption characteristic, the numerator is the sum of the bands demonstrating the shoulders and the denominator is the band positioned adjoining the absorption feature minimum [67]. Therefore, the absorption intensities attributed to Al-OH, Fe,Mg-OH, Si-OH and CO3 can be formulated for mapping advanced argillic, argillic-phyllic, propylitic and hydrous-silica alteration zones [35].

Principal Component Analysis
Principal Component Analysis (PCA) is a statistical approach that broadly and successfully used for decorrelation and enhancing the spectral contrast in remote sensing imagery [71]. This method transforms a number of correlated variables into several uncorrelated variables that termed PCs [72]. The eigenvector loadings (uncorrelated linear combinations) of variables were selected in a consistent way that each PC contains a smaller variance of extracted linear combination, sequentially [71,73]. The eigenvector loadings include key information linked to spectral features, which are anticipated from spectral bands of a remote sensing sensor [74]. For instance, a PC contains strong eigenvector loadings for indicative bands (reflection and absorption bands) of an alteration mineral with opposite signs enhances that mineral as bright pixels (if loading is positive in reflection band) or dark pixels (if loading is negative in reflection band) in the PC image [74,75].
In this study, the PCA method was implemented to some selected bands of Landsat-7 ETM+, Landsat-8 and ASTER using a covariance matrix for mapping hydrothermal alteration minerals. For identifying iron oxide-affected zones (gossan), bands 1, 3, 4 and 5 of Landsat-7 ETM+, bands 2, 4, 5 and 6 of Landsat-8 and bands 1, 2, 3 and 4 of ASTER were selected. The selected bands cover the iron oxide/hydroxide spectral properties in the VNIR region [3][4][5]. The eigenvector matrix for the selected bands and satellite sensors for mapping iron oxide/hydroxides are shown in Table 3A-C. Bands 1, 4, 5 and 7 of Landsat-7 ETM+, bands 2, 5, 6 and 7 of Landsat-8 and bands 1, 3, 4 and 6 of ASTER were used for detecting hydroxyl-bearing minerals. These bands cover the reflectance and absorption features of OH-minerals in the VNIR and SWIR regions [1,2]. Table 4A-C shows the eigenvector matrix for the selected bands and satellite sensors for mapping hydroxyl-bearing minerals. The reflectance properties and absorption intensities related to Al-OH, Fe, Mg-OH and CO3 can be mapped by ASTER VNIR+SWIR bands [23,35,38]. Bands 1, 4, 6 and 7 of ASTER were utilized for mapping advanced argillic zone. Bands 1, 3, 5 and 6 of ASTER were executed to detect argillic-phyllic zone. Bands 1, 3, 5 and 8 of ASTER were performed for discriminating propylitic alteration zone. Table 5A-C shows eigenvector matrix for the selected bands of ASTER for mapping advanced argillic, argillic-phyllic and propylitic alteration zones. After implementing the algorithms for all band ratios and PCAs, firstly the obtained DN values were normalized, then the X+3S was used to obtain definite anomaly. It means all the DN values showing the number more than the X+3S have been considered as target alteration minerals and zones.

Bayesian Networks Model
A Bayesian network is an interpreted directed acyclic graph (DAG), which is able to model uncertain relationships between variables in a complex system [76][77][78][79]. The mathematical concepts of the Bayesian networks model can be summarized as follows [43,77]. The subclass x belongs to a class of a set of classes ω1, ω2, …, ωn, if a class is defined by the highest conditional probability. The conditional probability is calculated using Equation (1): where P(x) is the non-conditional probability and P(ωi) is the prior probability of each class. The prior probability is calculated by dividing the number of samples in each class by the total number of samples [43]. In this method, a probability distribution function (PDF) is assigned for each class. Then, the training data is exploited to estimate the parameters involved in the PDF. The covariance matrix and the mean vector are calculated as the parameters of a Gaussian probability function provided that the data is normally distributed [76]. In other words, it is mathematically formulated as follows: In this equation (Equation (2)), m is the number of variables, which is added to μi and Σi of the mean vector and an m*m covariance matrix of the ith class that calculated using Equations (3) and (4): Bayesian networks model uses a structural graph known as a DAG to represent the knowledge about different domains or random variables [41]. The DAG is defined by the nodes and the directed edges. The former and the latter represent random variables and the relationship among variables, respectively, as it is shown in Figure 3. As can be seen from the direction of the arrow in Figure 3, there is a direct relationship between xi and xj. The xi (known as the parent node) is a dependent variable of the xj (known as an offspring node) [43]. There are different forms of Bayesian networks (See Reference [41] and references therein). One of the most popular forms of Bayesian networks is Naive Bayes (NB) classifier [80,81]. It is a simple structured algorithm with a single parent node and a number of offspring nodes [76,79,80]. A typical NB classifier diagram is shown in Figure 4. It is not only straightforward and easy to construct but also, no training procedure is required in the NB classifier [81]. The NB classifier undertakes comprehensive conditional independence between characteristics, which is impracticable for several predictor patterns utilized in mineral potential mapping [41]. In this study, the NB classifier was used for fusing the thematic layers derived from Landsat-7 ETM+, Landsat-8 and ASTER satellite sensors for generating a mineral potential map for the Ahar-Arasbaran region.

Fieldwork Data and Laboratory Analysis
The locations of hydrothermal alteration zones and their spatial relation with epithermal gold mineralization were systematically investigated using Global positioning system (GPS) survey in the study area (several field campaigns from June to August 2018). A handheld GPS (Garmin, Etrex Vista Hcx), with an average accuracy of 7 m, was used to record the hydrothermal alteration locations. Numerous field photographs and rock samples (120 samples) were collected from the alteration zones and ore mineralization. Rock samples were utilized for laboratory analysis to prepare thin and polished sections of altered rocks and ore mineralization as well as X-ray diffraction (XRD) analysis. Mineralogical compositions were analyzed using an Asenware AW-XDM 300 X-ray diffractometer (voltage: 40 Kv, current: 30 mA, step time: 1s and step size: 0.05° 2θ) at the Zarazma Mineral Studies Company, Tehran, Iran. Besides, the confusion matrix (error matrix) and Kappa Coefficient were calculated for hydrothermal alteration mineral mapping derived from remote sensing analysis versus field data.

Generating Thematic Layers Using Multi-Sensor Remote Sensing Data
Figure 5A-C shows iron oxide/hydroxide zones (gossan) derived from 3/1 band ratio of Landsat-7 ETM+, 4/2 band ratio of Landsat-8 and 2/1 band ratio of ASTER, respectively. Figure 5A shows the spatial distribution of iron oxide/hydroxide minerals derived from the Landsat-7 ETM+ band ratio as red pixels. Most of the documented gold mineralizations are associated with iron oxide/hydroxide zones (gossan), especially in the northern and northeastern parts of the study area. The spatial distribution of iron oxide/hydroxide minerals in the Landsat-8 band ratio image ( Figure 5B) is almost similar to Landsat-7 ETM+ resultant image. But, it is extensive in some locations in the northwestern and southeastern parts of the selected subset scene. Figure 5C shows the ASTER band ratio resultant image. The surface abundance of iron oxide/hydroxides in this image is lower compared to the Landsat-7 ETM+ and Landsat-8 results. However, the high concentration of iron oxide/hydroxides was mapped in the northwestern part of the study area using the ASTER band ratio ( Figure 5C). Regarding the geological map of the Ahar-Arasbaran region (see Figure 1), iron oxide/hydroxide minerals were mapped along with geological lineament features and igneous rocks (granite, granodiorite, biotite granite, andesite, dasite and basalt), volcano sedimentary units and massive and bedded limestone.
Typically, hydroxyl-bearing (Al-OH and Fe,Mg-OH) minerals and carbonates zones were mapped in Figure 6A-C using the 5/7 band ratio of Landsat-7 ETM+ (A), 6/7 band ratio of Landsat-8 (B) and 4/9 band ratio of ASTER (C). The green pixels depict OH-alteration and carbonates, which are normally concentrated in igneous rock (granite, granodiorite, biotite granite, andesite and dasite), volcano sedimentary units and limestone. The OH-alteration minerals are more strongly mapped in the Landsat-8 and ASTER resultant images compared to Landsat-7 ETM+ image ( Figure 6A-C). Almost all of the documented gold occurrences have an adjoining spatial relationship with hydroxyl-bearing alteration minerals; it is particularly observable in the Landsat-8 resultant image ( Figure 6B). It may be due to the high signal to noise radiometer performance of Landsat-8 data, which allows detecting subtle variation in surface conditions [58].
ASTER band ratios were used to specifically map the surface distribution of hydrothermal alteration zones in the study area. Figure 7A-C shows the advanced argillic alteration zone derived from 4/6 (A), the argillic-phyllic alteration zone derived from 5/6 (B) and the propylitic alteration zone derived from 5/8 (C), respectively. Concerning the geology map of the study area (see Figure 1), the advanced argillic alteration zone is corresponded to igneous, volcano sedimentary units and limestone; the argillic-phyllic alteration zone is associated with granite, granodiorite, andesite, dasite, rhyolite, trachyte, limestone units and sedimentary rocks; the propylitic alteration zone is typically concentrated with andesite, dasite, volcano sedimentary units and limestone (Figure 7A-C). The high surface abundance of argillic-phyllic and propylitic alteration zones was mainly mapped in the northwestern part of the study area. The spatial distribution of the advanced argillic alteration zone ( Figure 7A) is intensely matched with hydroxyl-bearing mineral zones that mapped by Landsat-7 ETM+ and Landsat-8 band ratio images ( Figure 6A,B). The documented gold mineralizations have closer spatial relationship with the advanced argillic alteration zone compared to the argillic-phyllic and propylitic alteration zones in the study area.   Detailed mapping of advanced argillic, argillic-phyllic, propylitic and hydrous silica-affected alteration zones was obtained using the RDB1 (4 + 6/5), RDB2 (5 + 7/6), RDB3 (6 + 9/7 + 8) and RDB4 (5 + 8/6 + 7) of ASTER ( Figure 8). Red pixels show advanced argillic zones, which are mostly distributed in the eastern and southeastern parts of the selected subset scene. Comparison to the geological map of the study area (see Figure 1), suggests that the advanced argillic zones are typically associated with granite and granodiorite rocks. Some of the documented gold mineralizations show close spatial relationship with the advanced argillic zones, especially in the eastern part of the study area ( Figure 8). Argillic-phyllic alteration zone depict as green pixels. This alteration zone is distributed in many parts of the study area, which are normally associated with andesite, dasite, volcano sedimentary units and sedimentary rocks (e.g., sandstone, siltstone, marl and conglomerates). Due to high content of detrital clays (montmorillonite, illite and kaolinite) in the sedimentary units, argillic-phyllic alteration zone could also be mapped with exposures of sedimentary rocks [35]. The surface abundance of hydrous silica-affected alteration zone (blue pixels) is low and mostly detected in the southwestern and northwestern parts of the study zone ( Figure 8). The hydrous silica zone was commonly identified with sedimentary units (conglomerates and sandstone), although this alteration zone is correspondingly adjacent to some of the gold mineralization zones in the northwestern part of the study area. Propylitic zone (yellow pixels) was strongly mapped in the selected subset scene (Figure 8). With regard to the geology map of the study area (see Figure 1), the spatial distribution of the propylitic zone typically corresponds with massive and bedded limestone, volcano sedimentary units and intermediate to mafic igneous rocks. It is because carbonates and alteration products of mafic minerals contain the strong contribution of CO3 and Fe,Mg-OH mineral groups, which produce similar spectral features to propylitic alteration zone. However, this alteration zone is one of the dominant mineral assemblages that mapped near to the gold mineralization zones, especially in the northwestern and northern parts of the study area ( Figure 8). Figure 8. The RDB1(4 + 6/5), RDB2 (5 + 7/6), RDB3 (6 + 9/7 + 8) and RDB4 (5 + 8/6 + 7) of ASTER shows advanced argillic, argillic-phyllic, propylitic and hydrous silica-affected alteration zones in the study area overlaid on hill shade. Table 3 shows the eigenvector loadings derived from PCA for mapping iron oxide/hydroxides (gossan) using bands 1, 3, 4 and 5 of Landsat-7 ETM+, bands 2, 4, 5 and 6 of Landsat-8 and bands 1, 2, 3 and 4 of ASTER. Analyzing the eigenvector loadings for Landsat-7 ETM+ selected bands (1, 3, 4 and 5) indicates that the PCA3 contains unique contribution (magnitude and sign of eigenvector loadings) of iron oxide/hydroxide minerals. The PCA3 has moderate loadings of band 1 (0.420) and strong loadings of band 5 (−0.780) with opposite signs (Table 3A). Band 1 (0.45-0.52 μm) of Landsat-7 ETM+ is positioned at absorption features of iron oxide/hydroxides (band 1 is considered an absorption band herein), while band 5 (1.55-1.75 μm) of Landsat-7 ETM+ is positioned at reflectance properties of iron oxide/hydroxides (band 5 is considered a reflection band herein). Thus, iron oxide/hydroxide minerals appear as dark pixels in the PCA3 due to negative sing in the reflection band (band 5), which were subsequently converted to bright pixels by negation. Figure 9A shows the resultant PCA3 image. Iron oxide/hydroxide minerals (red pixels) are mainly represented in the northern and northwestern parts of the study area, which are associated with some of the gold occurrences. However, a number of epithermal gold mineralizations do not show the spatial relationship with high abundance of iron oxide/hydroxide minerals, which are located in the southern and western parts of the study area.
Analysis of the eigenvector loadings of Landsat-8 selected bands (2, 4, 5 and 6) shows that the PCA2 can be used for mapping oxide/hydroxide minerals. The PCA2 contains moderate to strong contribution of bands 2 (−0.374) and 5 (−0.555) as absorption bands and strong contribution of band 6 (0.741) with a positive sign as reflection band (Table 3B). As a result, iron oxide/hydroxide minerals manifest as bright pixels in the PCA2 image ( Figure 9B). The spatial distribution of iron oxide/hydroxide minerals (red pixels) in Landsat-8 results is identical with Landsat-7 ETM+ PCA3 image but it is stronger in some parts, mainly in the southern and western sectors. Iron oxide/hydroxide minerals can be detected using the PCA3 derived from ASTER selected bands (1, 2, 3 and 4). The PCA3 shows moderate to strong loadings in absorption bands, including band 1 (−0.396), band 2 (−0.539) and band 3 (−0.259) with a negative sign and strong and positive loading in band 4 (0.695) as reflection band (Table 3C). Hence, iron oxide/hydroxide minerals represent bright pixels ( Figure 9C). A higher abundance of iron oxide/hydroxide minerals was mapped in the ASTER PCA2 image compared to Landsat-7 ETM+ and Landsat-8 PCA images, which is typically matched with most of the gold mineralizations. The pixels contain iron oxide/hydroxide minerals mapped by PCA images show a better spatial relationship with the gold mineralization zones compared to band ratio images (see Figure 5A-C). It indicates that the selective PCA can specially detect the alteration pixels in the spatial domain. Table  4 shows the eigenvector loadings derived from PCA for mapping hydroxyl-bearing minerals using bands 1, 4, 5 and 7 of Landsat-7 ETM+, bands 2, 5, 6 and 7 of Landsat-8 and bands 1, 3, 4 and 6 of ASTER. Considering the eigenvector loadings contain unique contribution of hydroxyl-bearing minerals in the absorption and reflection bands, it is discernible that the PCA4 includes the unique contribution of OH-minerals for all selected datasets (Table 4A-C). The PCA4 derived from Landsat-7 ETM+ selected bands (1, 4, 5 and 7) shows a strong negative loading in band 5 (−0.599) and a strong positive loading in band 7 (0.692) ( Table 4A). Because of negative loading in the reflection band (band 5), the hydroxyl-bearing minerals are represented as dark pixels in the PCA4, which are inverted to bright pixels by multiplication to −1, subsequently ( Figure 10A). Surface distribution of hydroxyl-bearing minerals (green pixels) depicts in the PCA4 image of Landsat-7 ETM+. The PCA4 derived from Landsat-8 selected bands (2, 5, 6 and 7) contains a strong negative loading of band 6 (−0.643) (the reflection band) and a strong positive loading of band 7 (0.694) (the absorption band) (Table 4B). Therefore, the PCA4 image was negated to depict the OH-minerals as bright pixels. Figure 10B shows the resultant image. For ASTER selected bands (1, 3, 4 and 6), the PCA4 has a strong positive contribution of band 4 (0.567) and a strong negative contribution of band 6 (−0.649) (Table 4C). Hence, the hydroxyl-bearing minerals appear as bright pixels in the PCA4 image. Figure 10C manifests the spatial distribution of the OH-minerals as green pixels in the PCA4 image of ASTER. Comparison of the PCA images to the band ratio images (see Figure 6A-C) suggests that the pixels detected in the selective PCA method show a closer spatial relationship to the gold mineralization zones and have a stronger manifestation in the image-maps.  Table 5 shows the eigenvector loadings for mapping advanced argillic, argillic-phyllic and propylitic alteration zones using ASTER bands such as bands 1, 4, 6 and 7 for the advanced argillic zone, bands 1, 3, 5 and 6 for the argillic-phyllic zone and bands 1, 3, 5 and 8 for the propylitic zone. Considering the magnitude and sign of eigenvector loadings for mapping advanced argillic zone (Table 5A), it is evident that the PCA3 contains spectral information to map advanced argillic zone due to a strong negative loading in band 4 (−0.549) and a strong positive loading in band 6 (0.509). Dark pixels depict the alteration zone due to a negative sign in the reflection band (band 4), which are afterward converted to bright pixels. The analysis of eigenvector loadings for mapping argillic-phyllic zone indicates that the PCA4 can mainly detect argillic-phyllic zone because of the strong contribution of bands 5 (−0.733) and 6 (0.679) with inverse signs (Table 5B). Muscovite (as a typical and dominant mineral in the phyllic zone) shows strong absorption in band 6 of ASTER, while lower absorption in band 5 of ASTER [35,37]. Thus, band 5 is assumed to be a reflection band and band 6 is considered as a strong absorption band herein. As a result, argillic-phyllic zone manifests as dark pixels due to negative sign in the band 5 (reflection band). Then, dark pixels were inverted to bright pixels by negation. Propylitic alteration zone can be mapped in the PCA4 image because of strong eigenvector loadings in band 5 (−0.712) and band 8 (0.696) with opposed signs (Table 5C). Herein, band 5 is pondered as reflection band and band 8 is deliberated as absorption band. Fe,Mg-OH and CO3 mineral groups (propylitic zone: chlorite, epidote and calcite) have high absorption properties in band 8 (2.295-2.365 μm) and reflection (very low absorption) features in band 5 (2.145-2.185 μm) of ASTER [35,39]. Accordingly, propylitic alteration zone appears as dark pixels that were negated to bright pixels in the PCA4 image. Figure 11 shows PCA image-map derived from the PCA3 image for advanced argillic mapping, the PCA4 image for argillic-phyllic zone mapping and the PCA4 image for propylitic zone mapping.
The spatial distribution of advanced argillic zones is stronger in the northeastern parts and weaker in southeastern part of the study area compared to RDBs image-map (see Figure 8). The advanced argillic zone resulting from PCA shows remarkable vicinity to the gold mineralization ( Figure 11). The argillic-phyllic zone shows nearly similar surface distribution to RDBs image-map. However, the high concentration of propylitic zone was mapped in the northwestern part of the study area in the PCA image-map compared to RDBs image-map (see Figure 8). On the other hand, the propylitic zone shows lower spatial distribution in the northeastern and southeastern parts of the PCA image-map ( Figure 11). Figure 11. The PCA image-map derived from ASTER selected PCAs for mapping advanced argillic, argillic-phyllic and propylitic alteration zones in the study area that overlaid on hill shade.

Fusing Thematic Layers Using Naive Bayes (NB) Classifier
The thematic layers of hydrothermal alteration zones derived from Landsat-7 ETM+, Landsat-8 and ASTER datasets were fused using the NB classifier to generate a mineral potential map for the Ahar-Arasbaran region. A DAG was designed for the thematic layers produced by image processing techniques in this study ( Figure 12). Eight distinct layers were employed as independent predictive layers, including iron oxide minerals derived from Landsat-7 ETM+, hydroxyl-bearing minerals derived from Landsat-7 ETM+, iron oxide minerals derived from Landsat-8, hydroxyl-bearing minerals derived from Landsat-8, iron oxide minerals derived from ASTER, advanced argillic alteration derived from ASTER, argillic-phyllic alteration derived from ASTER and propylitic alteration derived from ASTER. The DAG was used to integrate the predictor variables. It yields a posterior probability map showing the probability of gold mineralization occurrences. Subsequently, the following steps were taken to generate the posterior probability map. To train the DAG, 25 known gold mineralizations in the study area were selected as positive sites and 26 non-mineralized locations were selected as non-deposit (negative) sites, which have already been verified by field survey. In the next stage, the thematic layers (alteration image-maps) were resampled to a cell size of 150 * 150 m and a buffer zone of 300 m was considered around the positive and negative sites. The training data, the pixels superimposed by the positive and negative sites, containing a total of 468 pixels. Each pixel was considered as a vector of 8 arrays, including the values of 8 thematic (alteration) layers. To train the model, 70% of these pixels were used, while 30% of the pixels were used to validate the model generated. The calculation of the confusion matrix shows a total accuracy of 85.1%, which indicates that the model has been hypothesized and established. Having the trained NB model, all the data were used as the input of the model to generate the posterior probability map. However, the map is also required subsequent classification; thus the natural breaks algorithm [81] was used for classification of the posterior probability map. Three threshold values (produced by the foregoing algorithm) were used to generate a four-class map showing the probability of epithermal gold occurrences. The classes are highly probable (red), probable (green), moderately probable (yellow) and improbable (gray). As a result, a mineral potential map for the Ahar-Arasbaran region was produced ( Figure 13). Most of the known gold mineralizations are located in the highly probable (red) zone, although a small number of the gold occurrences can be seen in the probable (green) and moderately probable (yellow) zones. Many high probable zones in the northwestern, northern, northeastern, southeastern and southwestern parts of the study area contain high potential for undiscovered epithermal gold mineralizations ( Figure 13).

Verifying the Results Using Field Data and Laboratory Analysis
Several GPS surveys were carried out in different parts of the Ahar-Arasbaran region for verifying the mineral potential map and discovering new prospective zones of epithermal gold mineralizations, especially in highly probable zones. Numerous field photographs and rock samples were collected from different types of alteration zones related to gold mineralization such as advanced argillic, argillic-phyllic, propylitic and hydrous silica. In this investigation, some of the gold mineralization areas (highly probable zone), such as Zailig, Noghdouz, Javan-Sheikh, Nabi-Jan and Sonajil, were selected for a detailed field excursion, petrographic study and XRD analysis. The advanced argillic alteration, argillic-silica alteration, silica alteration and propylitic alteration were identified in the Zailig area ( Figure 14A-D). The advanced argillic alteration is the most extensive alteration zone in the vicinity of gold mineralizations in the Zailig area ( Figure 14A,B). The silica alteration is identified in the form of silica major clasts along with iron oxides ( Figure 14C). The other type of alteration zones is argillic-silica alteration, which is placed around the silica veins associated with gold mineralization ( Figure 14D). Figure 15A,B) shows microphotographs of argillic-silica alteration. Primary plagioclase replaced by sericite, clay minerals and jarosite ( Figure 15A). Recrystallized quartz and relics of plagioclase are surrounded by clay minerals ( Figure 15B). Propylitic alteration zone were also found as distal alteration zone in the Zailig area ( Figure 16A-D). Secondary minerals for instance, chlorite, epidote and calcite replaced original mineralogy (feldspars) as vesicular and amygdaloidal textures in the propylitic zone ( Figure 16B). Microphotographs of the propylitic zone show that the phenocrysts of plagioclase are replaced by chlorite, epidote and calcite ( Figure 16C,D). Quartz is phenocrystalline and anhedral in the background, while plagioclase is euhedral and partially replaced by epidote ( Figure 16C). The amygdaloidal texture is observable in Figure 16D, which amygdales are filled with calcite and quartz.   Some typical silicified and breccia (quartz veins) zones occur in altered granitic and andesitic rocks in the Noghdouz gold mineralization area (Figures 17 and 18). The specimens of silicified zone show breccia and clastic textures. The cement and major clasts of the breccia textures are composed of silicate minerals ( Figure 17A-C). Epithermal gold mineralization occurs in the breccia zone (quartz veins) in the altered granitic host rocks. This mineralization is also associated with advanced argillic alteration ( Figure 18A,B).  Iron oxide alteration zone (limonitic-hematite rocks) and oxidized breccia with banded chalcedonic quartz are identified in Javan-Sheikh gold mineralization area ( Figure 19A,B). The well-developed gossan covers (limonitic-hematite-silicic rocks) show rough and geologic relief features compared to surrounding altered rocks ( Figure 19A). The size of gossan covers are around 200 to 300 m that are surrounded by more extensive zones of propylitic and phyllic-argillic alteration zones. Although, silicified zone is also associated with the gossan covers, partially. The epithermal gold mineralization of the Nabi-Jan area is located in quartz-sulfide veins and developed at the top of an intrusive body of granodiorite ( Figure 20A,B). Gold mineralization is typically in the zones where intensely silicified and located in the advanced argillic alteration. In the Nabi-Jan area, the distal alteration zone is also propylitic alteration. The Sonajil gold mineralization occurs as a stock-work of thin quartz veins in granitoid rocks. The development of the argillic-phyllic alteration zone along with the siliceous zones and iron oxides were identified in the Kalijan area. Sphalerite, galena, chalcopyrite and pyrite are main sulfide mineralization associated with native gold mineralization ( Figure 21A-E). Quartz, iron oxide/hydroxide and minor calcite are gangue minerals.
Mineralogical compositions of hydrothermal alteration zones were investigated by XRD analysis. Thirty samples from different hydrothermal alteration zones were analyzed for this study. Representative XRD analysis of samples collected from the iron oxide/hydroxide alteration (gossan covers), advanced argillic alteration, argillic-phyllic alteration, propylitic alteration and hydrous silica alteration (silicified zone) are shown in this paper ( Figure 22A-E). Goethite, jarosite, gypsum and quartz are mineral phases that detected in the gossan cover ( Figure 22A). In the advanced argillic alteration ( Figure 22B), muscovite, illite, kaolinite, gypsum, orthoclase, albite and quartz are main mineralogical phases. The predominant minerals detected in the argillic-phyllic alteration are kaolinite, muscovite, illite, jarosite, albite and quartz ( Figure 22C). Epidote, chlorite, calcite, albite and quartz are identified in the propylitic alteration ( Figure 22D). Quartz, albite, jarosite, goethite, calcite, chlorite, gypsum and dolomite are observed in the silicified alteration zone ( Figure 22E).    In this analysis, confusion matrix and Kappa Coefficient [82][83][84][85][86] were used for assessing the accuracy of alteration mineral mapping derived from remote sensing analysis versus systematic GPS surveys collected from different alteration zones during fieldwork in the study area. Thirty representative GPS points were used for calculating the confusion matrix and Kappa Coefficient in this paper. Table 6 shows the locations of hydrothermal alteration zones recorded by a systematic GPS survey. Table 7 shows the confusion matrix for alteration mineral mapping versus field data. The results show the overall accuracy of 76.66% and Kappa Coefficient of 0.71. The advanced argillic alteration, argillic-phyllic and propylitic classes show the producer's accuracy of 83%, while the producer's accuracy for the iron oxide/hydroxide and hydrous silica classes is 67%. The highest user's accuracy is achieved for the argillic-phyllic and propylitic classes (100%), whereas the lowest user's accuracy is recorded for the iron oxide/hydroxide class (50%). The advanced argillic has the user's accuracy of 83% and hydrous silica class shows the user's accuracy of 67% (Table 7). Accordingly, the accuracy assessment results indicate that the alteration mineral mapping has appropriate match (overall accuracy: 76.66%) and very good degree of agreement (Kappa Coefficient: 0.71) with field data. However, some spectral mixing and confusion between alteration classes are also distinguishable. The iron oxide/hydroxide and hydrous silica classes show the highest feasibility for spectral mixing and confusion compared to other classes. The propylitic and argillic-phyllic classes contain the lowest spectral mixing and confusion. The advanced argillic class has some spectral mixing and confusion with the argillic-phyllic class.
Iron oxide/hydroxide zones (gossan cover), hydroxyl-bearing (Al-OH and Fe,Mg-OH) minerals and carbonates zones, advanced argillic, argillic-phyllic, propylitic and hydrous silica (silicified zone) alteration zones were mapped using band ratio, RBD and selective PCA image processing techniques. Using band ratios of 3/1 (Landsat-7 ETM+), 4/2 (Landsat-8) and 2/1 (ASTER) identify the spatial distribution of iron oxide/hydroxide zones, which are mainly associated with lineament features and igneous rocks, volcano sedimentary units and massive and bedded limestone (See Figure 5A-C). The documented epithermal gold occurrences mostly show close spatial locations with detected iron oxide/hydroxide zones. The PCA3 image of Landsat-7 ETM+ selected bands (1, 3, 4 and 5), the PCA2 image of Landsat-8 selected bands (2, 4, 5 and 6) and the PCA3 image of ASTER selected bands (1, 2, 3 and 4) were also represented iron oxide/hydroxide spatial distribution in the study area (see Figure  9A-C). The identified iron oxide/hydroxide zones using PCA are characteristically better matched with most of the gold mineralizations compared to band ratio images. Using band ratios of 5/7 (Landsat-7 ETM+), 6/7 (Landsat-8) and 4/9 (ASTER) detect the hydroxyl-bearing minerals and carbonates zones (see Figure 6A-C), which are generally matched with igneous rock, volcano sedimentary units and limestone. The gold mineralizations are typically located in the high abundance zones of hydroxyl-bearing/carbonate minerals. The advanced argillic, argillic-phyllic and propylitic alteration zones are mapped using ASTER band ratios of 4/6 (advanced argillic), 5/6 (argillic-phyllic), 5/8 (propylitic), respectively (see Figure 7A-C). The advanced argillic alteration shows closer spatial location with the gold mineralizations in comparison with the argillic-phyllic and propylitic alteration zones. The PCA4 image of Landsat-7 ETM+ selected bands (1, 4, 5 and 7), Landsat-8 selected bands (2, 5, 6 and 7) and ASTER selected bands (1, 3, 4 and 6) detects the surface distribution of hydroxyl-bearing minerals in the study area (see Figure 10A-C). The pixels detected in the PCA images show a stronger manifestation of OH-minerals compared to band ratio images and closer spatial relationship to the documented gold mineralization zones.
Implementing the RDB1 (4 + 6/5), RDB2 (5 + 7/6), RDB3 (6 + 9/7 + 8) and RDB4 (5 + 8/6 + 7) of ASTER reveal the advanced argillic, argillic-phyllic, propylitic and hydrous silica-affected alteration zones in the study area, comprehensively (see Figure 8). Some of the gold mineralizations in the eastern part of the study area are mainly situated in the advanced argillic zones. The hydrous silica zone was also mapped near some of the gold mineralization zones in the northwestern part of the study area. The propylitic zone is one of the main mineral assemblages associated with gold mineralization zones in the northern and northwestern parts of the study area. Only few gold occurrences were identified in the argillic-phyllic alteration zone. The PCA3 image derived from ASTER bands 1, 4, 6 and 7 for advanced argillic mapping, the PCA4 image derived from ASTER bands 1, 3, 5 and 6 for argillic-phyllic zone mapping and the PCA4 image derived from ASTER bands 1, 3, 5 and 8 for propylitic zone mapping show surface abundance of advanced argillic, argillic-phyllic and propylitic zone with some spatial discrepancies (see Figure 11) compared to RDBs image-map (see Figure 8). Notable vicinity to the documented gold mineralizations was mapped in the advanced argillic zone, which is detected with the PCA technique.
The produced thematic layers (see the DAG diagram in Figure 12) derived from band ratio and PCA image processing techniques are fused using the NB classifier. Consequently, a mineral potential map for the Ahar-Arasbaran region is produced (see Figure 13), which includes four classes such as highly probable, probable, moderately probable and improbable. Maximum numbers of the known gold occurrences are situated in highly probable class, while some of the gold mineralizations are located in the probable and moderately probable classes. Accordingly, several parts of the study area, such as the northwestern, northern, northeastern, southeastern and southwestern sectors, can be considered to be highly prospective zones for epithermal gold mineralizations and may contain undiscovered Au deposits (see Figure 13).
Detailed field expedition, petrographic study and XRD analysis in some of the prospective areas located in the highly probable zone show the presence of hydrothermal alteration zones associated with gold mineralizations. Extensive alteration mineral assemblages of the advanced argillic and argillic-silica alteration zones are found in the vicinity of gold mineralizations in the Zailig area (see Figure 14A-D). Microphotographs of argillic-silica alteration show that primary plagioclase replaced by sericite, clay minerals and jarosite and relics of plagioclase are surrounded by clay minerals. The distal alteration zone in the Zailig area is propylitic alteration zone, which contains chlorite, epidote and calcite that replaced original mineralogy (feldspars) as vesicular and amygdaloidal textures (see Figure 16A-D). In the Noghdouz area, gold mineralization is occurred in the breccia zone (quartz veins) in the altered granitic host rocks, which is associated with advanced argillic alteration (see Figures  17A-C and 18A,B). Limonitic-hematite rocks and oxidized breccia with banded chalcedonic quartz are identified in Javan-Sheikh gold mineralization area (see Figure 19A,B), which are surrounded by propylitic and phyllic-argillic alteration zones. In the Nabi-Jan area, gold mineralization is associated with quartz-sulfide veins hosted by granodiorite (see Figure 20A,B), which strongly silicified and placed in advanced argillic alteration. Development of the argillic-phyllic alteration zone associated with the siliceous zones and iron oxides in granitoid rocks were identified with gold mineralization in the Sonajil area. Native gold mineralization is associated with sphalerite, galena, chalcopyrite and pyrite (see Figure 21A-E).
The XRD analysis of rock samples collected from different alteration zones is verified the presence of hydrothermal alteration minerals, including (i) goethite, jarosite, gypsum and quartz in the gossan cover; (ii) muscovite, illite, kaolinite, gypsum, orthoclase, albite and quartz in the advanced argillic alteration, (iii) kaolinite, muscovite, illite, jarosite, albite and quartz in the argillic-phyllic alteration; (iv) epidote, chlorite, calcite, albite and quartz in the propylitic alteration (see Figure 22A-E). The accuracy assessment results show the overall accuracy of 76.66% and Kappa Coefficient of 0.71 for hydrothermal alteration mapping using remote sensing datasets. It indicates that the alteration mineral mapping contains a suitable match and a very good degree of agreement with field data. Analyzing the producer's accuracy and user's accuracy shows that some spectral mixing and confusion between alteration classes are also feasible, especially for iron oxide/hydroxide and hydrous silica alteration groups and the advanced argillic and the argillic-phyllic alteration groups. Accordingly, the mineral potential map produced in this study using multi-sensor remote sensing imagery and Bayesian network model is viable and can be broadly applicable for epithermal gold exploration in the Ahar-Arasbaran region.

Conclusion
This investigation was accomplished to produce a mineral potential map for prospecting epithermal gold mineralization in the Ahar-Arasbaran region, NW Iran using multi-sensor remote sensing satellite imagery (e.g., Landsat-7 ETM+, Landsat-8 and ASTER) and the Bayesian network model. Iron oxide/hydroxide zones, hydroxyl-bearing minerals and carbonates zones, advanced argillic, argillic-phyllic, propylitic and silicified alteration zones were mapped in the Ahar-Arasbaran region using band ratio, RBD and selective PCA image processing techniques. The NB classifier was successfully implemented to fuse the thematic layers of hydrothermal alteration zones derived from the multi-sensor satellite imagery. As a result, a mineral potential map for the Ahar-Arasbaran region was produced, which highlighted the prospective zones as highly probable, probable and moderately probable zones. The northwestern, northern, northeastern, southeastern and southwestern parts of the study area were considered high potential zones for epithermal gold mineralizations, which might have undiscovered epithermal gold deposits. The high potential zones were verified by field and laboratory analysis such as systematic GPS surveying, analyzing several microphotographs of hydrothermal alteration minerals and ore mineralization and XRD analysis of collected rock samples from alteration zones. The advanced argillic and argillic-silica alteration zones were typically found in the vicinity of gold mineralizations. However, limonitic-hematite rocks, oxidized breccia and propylitic alteration zones were also documented as high potential zones in the study area. The field and laboratory results verified that the mineral potential map of the Ahar-Arasbaran region successfully indicates the known epithermal gold mineralizations and several new high prospective zones in the study area. The approach developed in this study is a cost-effective technique that can be used for epithermal gold exploration in metallogenic provinces before costly geophysical and geochemical studies. Briefly, this study suggests that geostatistical techniques (e.g., Bayesian network model, Fuzzy model, Artificial Neural Network Model etc.) are valuable approaches to fuse thematic layers of the multi-sensor imagery for generating the remote sensing-based mineral potential map for metallogenic provinces. The mineral exploration community and mining companies can consider the remote sensing-based mineral potential map as an economical and cost-effective tool for mineral prospecting before pricey geophysical and geochemical surveys in the metallogenic provinces.