Fusion of Remote Sensing, Magnetometric, and Geological Data to Identify Polymetallic Mineral Potential Zones in Chakchak Region, Yazd, Iran

60-9-669-2166 Abstract: Exploration geologists are urged to develop new, robust, and low-cost approaches to identify high potential zones related to underground/unexplored mineral deposits because of increased depletion of ore deposits and high consumption of basic metal production industries. Fusing remote sensing, geophysical and geological data has great capability to provide a complete range of prerequisite data to accomplish this purpose. This investigation fuses remote sensing data, such as Sentinel-2 and Landsat 7, aerial magnetic geophysical data, and geological data for identifying polymetallic mineralization potential zones in the Chakchak region, Yazd province, Iran. Hydrothermal alteration mineral zones and surface and deep intrusive masses, hidden faults and lineaments, and lithological units were detected using remote sensing, aerial magnetic, and geological data, respectively. The exploratory/information layers were fused using fuzzy logic modeling and the multi-class index overlap method. Subsequently, mineral potential maps were generated for the study area. Some high potential zones of polymetallic mineralization were identiﬁed and veriﬁed through a detailed ﬁeld campaign and drilling programs in the Chakchak region. In conclusion, the fusion of remote sensing, geophysical, and geological data using fuzzy logic modeling and the multi-class index overlap method is a robust, reliable, and low-cost approach for mining companies to explore the frontier areas with identical geologic conditions that are alleged to indicate polymetallic mineralization potential.


Introduction
One of the most significant phenomena in geology is the presence of uncertainty. Hence, different modeling methods and fusing information layers are required to be used to achieve appropriate results in determining mineralization potential areas in the metallogenic provinces. A mineral exploration program is required for the collection of information from different sources, which can be derived from superficial to deep-seated phenomena related to ore mineralization. To achieve this goal, it is necessary to interpret and fuse the data collected from the target zone [1][2][3]. Such fusion requires the utilization of Geographic Information System (GIS) techniques. The GIS system is designed to collect, analyze, store, update, and display information in the finest possible way [4][5][6]. Various methods can be used to fuse the exploratory layers derived from different sources to identify the potential zones of mineralization. These methods are usually based on data (data-driven) or based on conceptual models of mineralization and expert knowledge

Geological Setting of the Study Area
The Chakchak region is located in the western region of the Central Iranian Zone, Yazd Block, Yazd province, near the ChakChak Mountain ( Figure 1A,B). The Chakchak Mountain is located approximately 47 km east of Ardakan and 69 km north of Yazd, and it is accessible by an asphalt road. The average height of the region is 1730 m. Figure 1B shows the geology map of the Chakchak area. The consolidation of the Central Iranian basement is attributed to the Late Precambrian Katangan/Pan-African orogenesis. In this zone, rocks of all ages, from Precambrian to Quaternary, and several episodes of orogeny, metamorphism, and magmatism can be distinguished. The Central Iranian Zone in a broad sense comprises the whole area between the North and South Iranian ranges [18,19]. Within the Iranian plate, the Central-East Iran microplate is surrounded by the Great Kavir Fault in the north, by the Nain-Baft Fault in the west and southwest, and by the Harirud Fault in the east. It is bordered by the Upper Cretaceous to Lower Eocene ophiolite and ophiolitic melange. The microplate consists of different structural components, including Kerman-Tabas Block, Yazd Block, and Anarak-Khur Block [19].

Geological Setting of the Study Area
The Chakchak region is located in the western region of the Central Iranian Zone, Yazd Block, Yazd province, near the ChakChak Mountain ( Figure 1A,B). The Chakchak Mountain is located approximately 47 km east of Ardakan and 69 km north of Yazd, and it is accessible by an asphalt road. The average height of the region is 1730 m. Figure 1B shows the geology map of the Chakchak area. The consolidation of the Central Iranian basement is attributed to the Late Precambrian Katangan/Pan-African orogenesis. In this zone, rocks of all ages, from Precambrian to Quaternary, and several episodes of orogeny, metamorphism, and magmatism can be distinguished. The Central Iranian Zone in a broad sense comprises the whole area between the North and South Iranian ranges [18,19]. Within the Iranian plate, the Central-East Iran microplate is surrounded by the Great Kavir Fault in the north, by the Nain-Baft Fault in the west and southwest, and by the Harirud Fault in the east. It is bordered by the Upper Cretaceous to Lower Eocene ophiolite and ophiolitic melange. The microplate consists of different structural components, including Kerman-Tabas Block, Yazd Block, and Anarak-Khur Block [19].
The study area is part of the Zyro Formation and is placed in the Yazd Block. Lithologically, the Zyro Formation consists of pyroclastic, volcanic, and volcanoclastic rocks. Based Remote Sens. 2022, 14, 6018 4 of 24 on stratigraphic observations, the age of the lithological units is the Lower Precambrian. Generally, the study area is composed of metamorphic, volcanic, and volcanoclastic rocks of rhyolite and rhyolite tuff. Alluvial sediments and old alluvial fans cover many parts of the study area ( Figure 2B). The age of alluvial sediments is the Quaternary. The Chakchak area is typically covered with sedimentary units (sandstone, shale, marl, limestone, and conglomerate), igneous units (granite, gabbro, rhyolite, basalt, and diabase), and tuff outcrops ( Figure 2B). The volcanic and intrusive rocks are partially altered in some parts [17]. The texture of the igneous lithological units in the study area is microporphyry. Diabase extrusions with iron oxide impurities and subvolcanic rocks such as quartz monzonite are exposed in some parts of the area. Andesitic tuffs and weathered rhyolites are also observable. Silica veins of low to medium thickness are presented in the western part of the study area. The major geological structure directions of the region are N-S, NW-SE, and NW-SE, which are often cut by younger NE-SW and NW-SE-trending lineaments. Polymetallic mineralizations such as pyrite, chalcopyrite, hematite, iron hydroxide, and titanium oxide were reported in the study area [17]. There are some copper indices, iron deposits, barite veins, kaolinite, and lead, which are quarried locally in the study area [17].
Remote Sens. 2022, 14, x FOR PEER REVIEW 4 of 26 The study area is part of the Zyro Formation and is placed in the Yazd Block. Lithologically, the Zyro Formation consists of pyroclastic, volcanic, and volcanoclastic rocks. Based on stratigraphic observations, the age of the lithological units is the Lower Precambrian. Generally, the study area is composed of metamorphic, volcanic, and volcanoclastic rocks of rhyolite and rhyolite tuff. Alluvial sediments and old alluvial fans cover many parts of the study area ( Figure 2B). The age of alluvial sediments is the Quaternary. The Chakchak area is typically covered with sedimentary units (sandstone, shale, marl, limestone, and conglomerate), igneous units (granite, gabbro, rhyolite, basalt, and diabase), and tuff outcrops ( Figure 2B). The volcanic and intrusive rocks are partially altered in some parts [17]. The texture of the igneous lithological units in the study area is microporphyry. Diabase extrusions with iron oxide impurities and subvolcanic rocks such as quartz monzonite are exposed in some parts of the area. Andesitic tuffs and weathered rhyolites are also observable. Silica veins of low to medium thickness are presented in the western part of the study area. The major geological structure directions of the region are N-S, NW-SE, and NW-SE, which are often cut by younger NE-SW and NW-SE-trending lineaments. Polymetallic mineralizations such as pyrite, chalcopyrite, hematite, iron hydroxide, and titanium oxide were reported in the study area [17]. There are some copper indices, iron deposits, barite veins, kaolinite, and lead, which are quarried locally in the study area [17].

Data Used
Geophysical, remote sensing, and geological data were used in this study. The aerial magnetometry data of the study region were obtained from the Geological Survey of Iran (GSI). Aerial magnetometry surveys were conducted using McFar USA aircraft in 1975. The flight altitude in this collection is 2.5 km, and the distance between the collection profiles is 7.5 km. Sentinel-2 and Landsat-7 remote sensing data were used in this study. Both remote sensing data are cloud-free and acquired in the summer season. Sentinel-2 scenes

Data Used
Geophysical, remote sensing, and geological data were used in this study. The aerial magnetometry data of the study region were obtained from the Geological Survey of Iran (GSI). Aerial magnetometry surveys were conducted using McFar USA aircraft in 1975. The flight altitude in this collection is 2.5 km, and the distance between the collection profiles is 7.5 km. Sentinel-2 and Landsat-7 remote sensing data were used in this study. Both remote sensing data are cloud-free and acquired in the summer season. Sentinel-2 scenes covering the study were acquired (date of acquisition: 28 August 2016) from the European Space Agency (Copernicus Open Access Hub; https://scihub.copernicus. eu/dhus/#/home, accessed on 20 Janaury 2022), which are Level-1C top-of-atmosphere (TOA) reflectance data. They include radiometric and geometric corrections along with orthorectification [20]. Landsat-7 (ETM+) scenes (collection level 2) were obtained (date of acquisition: 22 August 2002) from USGS website (earthexplorer.usgs.gov, accessed on 20 Janaury 2022). Landsat Collection 2 includes scene-based global Level-2 surface reflectance and surface temperature science products. In this study, it was required to combine two scenes to obtain full coverage of the study area. For this purpose, two ETM+ scenes and two Sentinel-2 scenes were used, and the images from each dataset were mosaicked. Sentinel-2 and Landsat-7 data were successfully utilized to identify hydroxylbearing minerals, which are found in the alteration zones associated with ore mineralization in the metallogenic provinces [21,22]. Geological fieldwork was conducted in the study area during August 12th to 15th, the year 2021 and September 20th to 24th, the year 2022. Field photographs and rock samples of hydrothermal alterations and ore mineralizations were collected. A Global Positioning System (GPS) survey was carried out in many locations around the hydrothermal alteration zones and ore mineralizations.

Methodology
In this study, Sentinel-2 and Landsat-7 remote sensing data were used to extract hydrothermal alteration zones. Aeromagnetic geophysical data were utilized to investigate deep intrusive masses. Hidden fault systems in the area were also identified from the magnetic data [22]. The geophysical data used in this analysis provided information for potential mineralization in different levels (superficial, medium, and deep) in the study area. Geological maps of the region were analyzed to reveal the fault systems, geological units, mine locations, and mineral occurrences. Subsequently, the information layers were fused using the fuzzy logic model and the multi-class index overlap method. As a result, two potential maps were produced for the study area. Figure 2 shows the methodological flowchart used in this study.

Preprocessing
The Internal Average Relative Reflectance (IARR) [23] method was implemented in both Sentinel-2 and Landsat-7 data for atmospheric corrections. The IARR uses the relative average of internal reflection: the average reflection of image pixels is calculated, and the radiation intensity values of each pixel are divided by the average of the image pixels. Therefore, the number obtained for each pixel will be the relative and normalized reflection spectrum of that pixel [24]. Bands 1, 2, 3, 4, 11, and 12 of Sentinel-2 and bands 1, 2, 3, 4, 5, and 7 of Landsat-7 were selected and used in this study. These bands contain spectral information for detecting hydrothermal alteration minerals (e.g., iron oxide/hydroxides, clay minerals, and carbonates) [21,25].
For aerial magnetometry data, several corrections were required to be applied to the data for profile positioning and processing networks of total magnetic field strength. The initial processing section includes reviewing and modifying the raw data, determining coordinates of the data, parallax corrections, removing the component attributed to the area of the earth, aligning the data, and deleting all the remaining errors, namely, fine-tuning, networking, and packing. Appropriately, the international geomagnetic reference field (IGRF) was initially corrected in this analysis. The earth's regional field model, or IGRF, was used to remove the earth's magnetic field effect.

Processing Techniques
For processing the remote sensing data (Sentinel-2 and Landsat-7), band ratioing (BR) [25], principal component analysis (PCA) [26] and least-squares fitting (LS-Fit) [9] were implemented to identify hydrothermal alteration zones and minerals. BR is normally used to enhance the spectral variances between bands and to diminish the influences of topography. Allocating one spectral band by another generates an image that provides relative band intensities. The image enhances the spectral discrepancies between bands [25].
Principal component analysis (PCA) produces uncorrelated output bands, to separate out noise components and to diminish the dimensionality of data. Highly correlated remote sensing data transform to uncorrelated output bands by running PCA. This is accomplished by generating a new set of orthogonal axes that have their origin at the data mean and that are rotated so the data variance is maximized [26]. LS-Fit executes a linear band prediction using least-squares fitting, which can identify zones of anomalous spectral response in a dataset. It computes the covariance of the input data and utilizes it to anticipate the selected band as a linear combination of the predictor bands plus an offset. The variance between the actual band and the modeled band is considered and output as an image. Pixels with a large residual (positive and negative) specify the existence of a feature that was not predicted (an absorption band). The modeled band image is also included in the output [9].
For aerial magnetometry data, filters are applied to the corrected data to process the conversions. Data processing with a filter changes amplitude or phase corresponding to the set of sinuous waves that make up the data of each network or profile. These conversions can be returned to the location domain through multiple Fourier transforms of the data in the frequency domain to determine the desired amplitude and then conversion of the modified data spectrum. The filters used in this study include trend removal filter, vertical derivative, analytical signal, tilt derivative, and horizontal derivative. Subsequently, a digital map of the lineaments and deep intrusive masses in the region was obtained from integrating the obtained data.
To separate the aerial magnetometry anomalies, the residual map was generated and examined. This map is obtained by fitting a level (first, second, third, or more degree) as a page or procedure and removing a data flow. The next step is to use a vertical derivative filter. The first vertical derivative (vertical gradient) is the physical equivalent of simultaneously measuring a magnetic field at two vertical points on top of each other, subtracting data, and dividing the results by spatially separating the vertical points of the collection. The second vertical derivative is a vertical gradient of the first vertical derivative. This derivative increases the resolution of high frequencies compared to low frequencies, and this feature is the basis of the application of derivative processing, which eliminates the effects of the long-wavelength region and eliminates the effects of adjacent anomalies. Subsequently, the analytical signal filter was applied. The analytical signal is a function that connects the magnetic field with the derivative [27]. The Equation (1) of this signal is described as follows: (1) In this equation, m is a magnetic anomaly. The dependence of the analytical signal is on the direction of magnetization and the direction of the Earth field. Using analytical signals, masses with the same geometry receive a parity signal. Simultaneously with the symmetry of the maxima, narrow masses form at the top of the edges of the broad masses and just above the center. With the help of an analytical signal, it is also possible to determine the position of a magnetic source, regardless of the remaining magnetism of various sources.
Another technique for determining the edges of anomalies is to apply a tilt angle derivative filter (TDX) [28], which is defined as Equation (2): Concerning the tilt derivative filter, TDX is the tilt angle, and f is the potential field of the collection. This filter has advantages, such as being dimensionless and using tilt angle changes on top of mineral masses to identify them. This analysis used this method to identify deep lines and faults resulting from magnetic data.
In the next step, a horizontal derivative filter was used. The magnitude of the horizontal gradient of the whole potential field, also called the horizontal gradient, is obtained from Equation (3) [29]: In this relation, T is the potential field, ( ∂T ∂x ) and ( ∂T ∂y ) are the horizontal derivatives of the potential field data in the x and y directions, respectively. In this research, the horizontal derivative filter was used to identify faults in the study area. This filter was applied to the data in different directions.

Fusion of the Datasets
Geospatial phenomena and processes are often ambiguously defined and inherently fuzzy. The variables of geospatial processes are scantily confined and regularly defined in terms of linguistic values. Hence, the fuzzy set theory has comprehensive applications in modeling and predicting geospatial results. Additionally, fuzzy sets are applicable for defining indistinct parameters. Fuzzy logic modeling [30] is based on the fuzzy set, which is a form of many-valued logic wherein the truth values of variables might be any real number between 0 and 1 [31]. Fuzzy logic can be considered a practice in which different layers are weighted at different levels, and expert knowledge has a significant role in this process. Various factors such as host rock, structures, geophysical anomalies, remote sensing studies, and documented fieldworks can be considered information layers (fuzzy members) for mineral exploration. Hence, each information layer might contain a specific weight between zero to one. It can be allocated depending on its importance and relation to ore mineralization. The number 1 indicates the degree of complete fuzzy membership, whereas the 0 indicates the absence of fuzzy membership. In some layers, depending on their radius of impact, they may be weighted up to several meters around the desired layer. After determining the mineralization control factors, the data layers are integrated by the fuzzy operators. Data layers are reclassified and networked after classification, data collection, and preparation of the existing data layers.
Typically, fuzzy logic modeling contains three major stages that would be executed on the information data ( Figure 3): (1) fuzzification (encoding); (2) logical integration using inference engine and fuzzy set operations; and (3) defuzzification (decoding) [32,33]. The input data can be categorical or numerical. All the presented data are independently imported, then a fuzzifier converts them into fuzzy values of membership, using a membership function µ A (x), in a continuous range between 0 (also called "false", "non-membership", and "incomplete") and 1 (also called "true", "full membership", and "complete"), respectively. Generally, µ A (x) is specified based on a priori knowledge of the geological system or using training data ( Figure 3) [32,33]. The fuzzy membership values characterize the importance of an input geospatial variable in the geospatial process being modeled. The distinct fuzzy sets are then fused using several parallel and/or serial networks that consecutively fuse fuzzy sets using fuzzy operators. The pattern of the inference engine would imitate the geospatial process being modeled [32].  The multi-class overlap index method is an extended version of the binary index overlap modeling method [35][36][37][38][39]. Each of the j patterns (classes) of the i control map has a Sij weight to indicate the degree of dependence on mineral deposits. The score assigned to the patterns can be a positive integer [35]. There is no restriction to the range of points allocated to the patterns, except that all control maps' range of points assigned to the patterns must be consistent. In other words, they must have the same minimum and maximum values. It is not practical to control the relative advantages of control maps by considering the range of different weight changes for one control map compared to other control maps. The relative importance of a control map compared to other control maps is controlled by the allocation of Wi weights. This score is usually a positive integer. Weighted control maps are then combined using Equation (4) [36,40], which calculates the average weights for each location in Equation (4): where S is the weighted average output value for each location, equal to the sum of Sij multiplied by class j from the control map i and Wi by the weight of the control map [36,41]. Therefore, in addition to flexibility in assigning weight to control patterns, the advantage of the multi-class index overlap method compared to modeling with the binary index overlap method is that it also considers uncertainty and predicts intervention. The weakness of this method is the linear integration of evidence that makes the role of effective processes in mineralization not directly considered.

Remote Sensing Results
For mapping alteration minerals and zones in the study area, BR, PCA, and LS-Fit techniques were implemented in Sentinel-2 and Landsat-7 remote sensing data. To characterize iron oxide/hydroxides (hematite, goethite, and jarosite), the VNIR spectral bands contain the most important information due to electronic transitions of Fe 3+ /Fe 2+ in the VNIR region from 0.45 to 1.2 µm [42,43]. Accordingly, a 4/2 band ratio of Sentinel-2 and a band ratio of 1/3 of Landsat-7 were selected to highlight iron oxide/hydroxides. Hydroxylbearing (Al-OH) alteration and carbonates (muscovite, kaolinite, gypsum, calcite, and dolomite) show spectral absorption features in 2.1-2.5 µm due to overtones and combinations of the fundamental vibrations [44], whereas their spectral reflectance typically occurs in 1.55-1.75 µm in the SWIR regions. These characteristics are matched with band 12 (2.100-2.280 µm) and band 11 (1.565-1.655 µm) of Sentinel-2, as well as band 7 (2.11-2.29 µm) and band 5 (1.57-1.65 µm) of Landsat-7, respectively. Therefore, the 11/12 band ratio of Sentinel-2 and 5/7 band ratio of Landsat-7 were used to map hydroxyl bearing altera- The multi-class overlap index method is an extended version of the binary index overlap modeling method [35][36][37][38][39]. Each of the j patterns (classes) of the i control map has a S ij weight to indicate the degree of dependence on mineral deposits. The score assigned to the patterns can be a positive integer [35]. There is no restriction to the range of points allocated to the patterns, except that all control maps' range of points assigned to the patterns must be consistent. In other words, they must have the same minimum and maximum values. It is not practical to control the relative advantages of control maps by considering the range of different weight changes for one control map compared to other control maps. The relative importance of a control map compared to other control maps is controlled by the allocation of Wi weights. This score is usually a positive integer. Weighted control maps are then combined using Equation (4) [36,40], which calculates the average weights for each location in Equation (4): where S is the weighted average output value for each location, equal to the sum of S ij multiplied by class j from the control map i and W i by the weight of the control map [36,41]. Therefore, in addition to flexibility in assigning weight to control patterns, the advantage of the multi-class index overlap method compared to modeling with the binary index overlap method is that it also considers uncertainty and predicts intervention. The weakness of this method is the linear integration of evidence that makes the role of effective processes in mineralization not directly considered.

Remote Sensing Results
For mapping alteration minerals and zones in the study area, BR, PCA, and LS-Fit techniques were implemented in Sentinel-2 and Landsat-7 remote sensing data. To characterize iron oxide/hydroxides (hematite, goethite, and jarosite), the VNIR spectral bands contain the most important information due to electronic transitions of Fe 3+ /Fe 2+ in the VNIR region from 0.45 to 1.2 µm [42,43]. Accordingly, a 4/2 band ratio of Sentinel-2 and a band ratio of 1/3 of Landsat-7 were selected to highlight iron oxide/hydroxides. Hydroxyl-bearing (Al-OH) alteration and carbonates (muscovite, kaolinite, gypsum, calcite, and dolomite) show spectral absorption features in 2.1-2.5 µm due to overtones and combinations of the fundamental vibrations [44], whereas their spectral reflectance typically occurs in 1.55-1.75 µm in the SWIR regions. These characteristics are matched with band 12 (2.100-2.280 µm) and band 11 (1.565-1.655 µm) of Sentinel-2, as well as band 7 (2.11-2.29 µm) and band 5 (1.57-1.65 µm) of Landsat-7, respectively. Therefore, the 11/12 band ratio of Sentinel-2 and 5/7 band ratio of Landsat-7 were used to map hydroxyl bearing alteration minerals and carbonates in the study area. Figure 4A-D show the spatial distribution of iron oxide/hydroxides and hydroxylbearing (Al-OH) alteration and carbonates using the selected band ratios. Considering Figure 4A,B, high to moderate surface distribution of iron oxide/hydroxides is typically mapped, associated with igneous rocks, the lithological units of sandstone, conglomerate, limestone, and marl. In addition, some parts of the Quaternary background showed high to moderate iron oxide/hydroxides. Regarding Figure 4C,D, high to moderate spatial distribution of hydroxyl-bearing (Al-OH) alteration and carbonates are identified within igneous rocks, limestone, marl, sandstone, and some parts of the Quaternary background.
Remote Sens. 2022, 14, x FOR PEER REVIEW 9 of 26 Figure 4A-D show the spatial distribution of iron oxide/hydroxides and hydroxylbearing (Al-OH) alteration and carbonates using the selected band ratios. Considering Figure 4 A and B, high to moderate surface distribution of iron oxide/hydroxides is typically mapped, associated with igneous rocks, the lithological units of sandstone, conglomerate, limestone, and marl. In addition, some parts of the Quaternary background showed high to moderate iron oxide/hydroxides. Regarding Figure 4C,D, high to moderate spatial distribution of hydroxyl-bearing (Al-OH) alteration and carbonates are identified within igneous rocks, limestone, marl, sandstone, and some parts of the Quaternary background. A PCA image with strong eigenvector loading for indicative reflection and absorption bands (in VNIR and SWIR regions) of a specific mineral or mineral groups (e.g., iron oxide/hydroxide minerals or clay minerals) with opposite signs detects that mineral/mineral groups as bright or dark pixels. If the eigenvector loading is positive in the reflection band of the target mineral the PCA image tone will be bright, and if it is negative, the PCA image tone will be dark for the target mineral [26]. The information typically characterizes in a very small fraction of the total information content for the original bands. After A PCA image with strong eigenvector loading for indicative reflection and absorption bands (in VNIR and SWIR regions) of a specific mineral or mineral groups (e.g., iron oxide/hydroxide minerals or clay minerals) with opposite signs detects that mineral/mineral groups as bright or dark pixels. If the eigenvector loading is positive in the reflection band of the target mineral the PCA image tone will be bright, and if it is negative, the PCA image tone will be dark for the target mineral [26]. The information typically characterizes in a very small fraction of the total information content for the original bands. After implementing the PCA transformation, it is anticipated that the eigenvector loadings specify the spectral signature of the target mineral or mineral groups in the reflection and absorption bands [21,24,26].
Analyzing the eigenvector loading derived from PCA (Table 1, panels A and B) indicates that some PCA images contain useful information for mapping iron oxide/hydroxides and hydroxyl-bearing (Al-OH) alteration and carbonates. The statistical results for Sentinel-2 PCA images show that the PCA2 and PCA4 can map the target alteration minerals. Eigenvector loadings in PCA2 show strong values in band 2 (−0.41) and band 4 (0.44) with opposite signs. Hence, iron oxide/hydroxides can be mapped as bright pixels in the PCA2 image. Strong eigenvector loadings in band 11 (0.78) and band 12 (−0.61) signified that PCA4 can identify hydroxyl-bearing (Al-OH) alteration and carbonates as bright pixels. Figure 5A,B show the resultant images of PCA2 and PCA4 for Sentinel-2 as a pseudo-color ramp. High to moderate surface distribution of iron oxide/hydroxides is typically detected associated with sandstone, conglomerate, marl, and igneous rocks ( Figure 5A). On the other hand, high to moderate values of hydroxyl-bearing (Al-OH) alteration and carbonates are mapped with igneous rocks, limestone, sandstone, and some parts of the Quaternary background ( Figure 5B). Analyzing eigenvector loadings for Landsat-7 PCA images indicates that the PCA3 and PCA4 (Table 1B) contain desired information for detecting iron oxide/hydroxides and hydroxyl-bearing (Al-OH) alteration and carbonates in the study area. The PCA3 has strong loadings of band 1 (0.58) and band 3 (−0.65) with opposite signs; therefore, the PCA3 image maps iron oxide/hydroxides as bright pixels. Eigenvector loadings in PCA4 contain a strong contribution from band 5 (0.64) and band 7 (−0.56). It is evident that the PCA4 image can map hydroxyl-bearing (Al-OH) alteration and carbonates. The pseudo-color ramp for PCA3 and PCA4 images for Landsat-7 are shown in Figure 5C,D. A high to moderate distribution of iron oxide/hydroxides is mapped in the Quaternary background, sandstone, some parts of igneous rocks, and marl ( Figure 5C). A high to moderate spatial distribution of hydroxyl-bearing (Al-OH) alteration and carbonates is identified associated with the limestone, igneous rocks, and some parts of the Quaternary background and sandstone ( Figure 5D). The PCA results are similar to BR outputs. Additionally, PCA mapped iron oxide/hydroxides and hydroxyl-bearing (Al-OH) alteration and carbonates with more details compared to BR (See Figures 4 and 5). Analyzing eigenvector loadings for Landsat-7 PCA images indicates that the PCA3 and PCA4 (Table 1(B)) contain desired information for detecting iron oxide/hydroxides and hydroxyl-bearing (Al-OH) alteration and carbonates in the study area. The PCA3 has strong loadings of band 1 (0.58) and band 3 (−0.65) with opposite signs; therefore, the PCA3 image maps iron oxide/hydroxides as bright pixels. Eigenvector loadings in PCA4 contain a strong contribution from band 5 (0.64) and band 7 (−0.56). It is evident that the PCA4 image can map hydroxyl-bearing (Al-OH) alteration and carbonates. The pseudocolor ramp for PCA3 and PCA4 images for Landsat-7 are shown in Figure 5C,D. A high to moderate distribution of iron oxide/hydroxides is mapped in the Quaternary background, sandstone, some parts of igneous rocks, and marl ( Figure 5C). A high to moderate spatial distribution of hydroxyl-bearing (Al-OH) alteration and carbonates is identified associated with the limestone, igneous rocks, and some parts of the Quaternary background and sandstone ( Figure 5D). The PCA results are similar to BR outputs. Additionally, PCA mapped iron oxide/hydroxides and hydroxyl-bearing (Al-OH) alteration and carbonates with more details compared to BR (See Figures 4 and 5).  The maps were generated using all selected bands of Sentinel-2 (1, 2, 3, 4, 11, and 12) and Landsat-7 (1, 2, 3, 4, 5, and 7) as input bands (predictor bands) and bands 4 and 12 (Sentinel-2) and bands 3 and 7 (Landsat-7) as modeled bands. These maps show the surface distribution of iron oxide/hydroxides and hydroxyl-bearing (Al-OH) alteration and carbonates as dark pixels. As a result, low to moderate values in the pseudo-color ramp show the surface distribution of the target alteration minerals. Figure 6A shows LS-Fit residual band 4 (output) for Sentinel-2. The surface abundance of iron oxide/hydroxides is depicted in igneous rocks, sandstone, limestone, and some parts of the Quaternary background ( Figure 6A). The LS-Fit residual band 12 (output) for Sentinel-2 is shown in Figure 6B. Hydroxyl-bearing (Al-OH) alteration and carbonates are mapped, associated with sand-stone, igneous rocks, marl, and some parts of the Quaternary background ( Figure 6B). Figure 6C shows LS-Fit residual band 3 (output) for Landsat-7. Iron oxide/hydroxides are detected with igneous rocks, sandstone, marl, and some parts of the Quaternary background ( Figure 6C). The LS-Fit residual band 7 (output) for Landsat-7 is displayed in Figure 6D. The spatial distribution of hydroxyl-bearing (Al-OH) alteration and carbonates is typically mapped in sandstone, igneous rocks, and limestone ( Figure 6D). The results of the LS-Fit model are identical to the PCA and BR outputs. Accordingly, the spatial distribution of the alteration minerals with some of the lithological units in the study area is identified and verified using different image processing techniques. The alteration minerals are one of the most important indicators of ore mineralizations, especially where they are associated with igneous lithological units in the study area.

Aerial Magnetometry Results
Figure 7A-C shows the residual map resulting from the application of the trend removal filter using the fit of a plate one-degree, two-degrees, and three-degrees for the study area, respectively. It should be noted that the zones having purple and blue shades in the maps indicate anomalous intrusive masses. In the zones where the dipoles are shaped, the anomalous masses are placed in the center of the dipoles. Figure 7A shows the residual map obtained from the fitting of a first-degree level. This Figure shows some

Aerial Magnetometry Results
Figure 7A-C shows the residual map resulting from the application of the trend removal filter using the fit of a plate one-degree, two-degrees, and three-degrees for the study area, respectively. It should be noted that the zones having purple and blue shades in the maps indicate anomalous intrusive masses. In the zones where the dipoles are shaped, the anomalous masses are placed in the center of the dipoles. Figure 7A shows the residual map obtained from the fitting of a first-degree level. This Figure shows some of the anomalous zones, especially in the northwestern, central, northeastern, and southeastern parts of the study area. Figure 7B shows an almost similar result to the one-degree residual map. The map presented in Figure 7C, which has the highest degree of fit of the curve to a third-degree level, shows the main and real anomalies that are closer to the surface more accurately. Regarding the geology map of the study area (see Figure 1B), igneous rocks (granite, gabbro, and diabase) are mainly placed in the center of the anomalous dipoles. Additionally, the other lithological units such as sandstone and marl, along with some parts of the Quaternary background, showed some anomalous dipoles ( Figure 7C). These sedimentary lithological units might superimpose some intrusive masses.  The map obtained by applying the first vertical derivative filter to the magnetic data is shown in Figure 8A. The vertical derivative maps amplify the high-frequency (lowwavelength) waves; therefore, surface anomalies become clearer, and the effect of deep anomalies will be reduced. As a result, high potential intrusive masses related to mineralization can be revealed. In this map, the analysis of magnetic masses can be observed due to the disappearance effect of large and deep anomalies. The anomalous dipoles are typically associated with igneous rocks (granite, gabbro, and diabase) and in northeastern, northwestern, and southwestern parts of the Quaternary background ( Figure 8A).  The map obtained by applying the first vertical derivative filter to the magnetic data is shown in Figure 8A. The vertical derivative maps amplify the high-frequency (low-wavelength) waves; therefore, surface anomalies become clearer, and the effect of deep anomalies will be reduced. As a result, high potential intrusive masses related to mineralization can be revealed. In this map, the analysis of magnetic masses can be observed due to the disappearance effect of large and deep anomalies. The anomalous dipoles are typically associated with igneous rocks (granite, gabbro, and diabase) and in northeastern, northwestern, and southwestern parts of the Quaternary background ( Figure 8A). Figure 8B shows the resulting map derived from applying the analytical signal filter. This Figure indicates the presence of masses with multiple magnetic anomalies and highlights the approximate position of the edges of the magnetic masses in the study area. In Figure 8B, the anomalous dipoles are mostly found with igneous rocks and the Quaternary background (northeastern, northwestern, and southwestern segments). The TDX was used to identify deep lines and faults resulting from magnetic data. Figure 8C shows the resultant map. Most of the long lineaments trend N-S, NW-SE, NE-SW in the study area. However, some short lineaments strike NE-SW, W-E, N-S are also mapped. is shown in Figure 8A. The vertical derivative maps amplify the high-frequency (lowwavelength) waves; therefore, surface anomalies become clearer, and the effect of deep anomalies will be reduced. As a result, high potential intrusive masses related to mineralization can be revealed. In this map, the analysis of magnetic masses can be observed due to the disappearance effect of large and deep anomalies. The anomalous dipoles are typically associated with igneous rocks (granite, gabbro, and diabase) and in northeastern, northwestern, and southwestern parts of the Quaternary background ( Figure 8A).  Figure 8B shows the resulting map derived from applying the analytical signal filter. This Figure indicates the presence of masses with multiple magnetic anomalies and highlights the approximate position of the edges of the magnetic masses in the study area. In Figure 8B, the anomalous dipoles are mostly found with igneous rocks and the Quaternary background (northeastern, northwestern, and southwestern segments). The TDX Furthermore, the horizontal derivative filter was used to identify fault systems in the study area. This filter was applied to the data in different directions. Figure 9A,B shows the application of horizontal derivative filters in the azimuth direction of 45 and 135 degrees. Deep and hidden faults were identified using the horizontal derivative map and TDX map. It also helped in identifying the orientation of the deep intrusive masses, which are mostly found in the central, northeastern, northwestern, and southwestern parts of the study area ( Figure 9A,B). Concerning the geology map of the study area (see Figure 1), regardless of igneous units, the magnetic anomalous dipoles are mostly placed on sandstone and the Quaternary background. It might show the presence of deep intrusive masses that were covered by sedimentary units. Accordingly, the prevalence of existing structural discontinuities in the deep intrusive masses can build conduit systems for circulation and entrapment of hydrothermal fluids and consequent ore mineralization in the study area. was used to identify deep lines and faults resulting from magnetic data. Figure 8C shows the resultant map. Most of the long lineaments trend N-S, NW-SE, NE-SW in the study area. However, some short lineaments strike NE-SW, W-E, N-S are also mapped. Furthermore, the horizontal derivative filter was used to identify fault systems in the study area. This filter was applied to the data in different directions. Figure 9A,B shows the application of horizontal derivative filters in the azimuth direction of 45 and 135 degrees. Deep and hidden faults were identified using the horizontal derivative map and TDX map. It also helped in identifying the orientation of the deep intrusive masses, which are mostly found in the central, northeastern, northwestern, and southwestern parts of the study area ( Figure 9A,B). Concerning the geology map of the study area (see Figure  1), regardless of igneous units, the magnetic anomalous dipoles are mostly placed on sandstone and the Quaternary background. It might show the presence of deep intrusive masses that were covered by sedimentary units. Accordingly, the prevalence of existing structural discontinuities in the deep intrusive masses can build conduit systems for circulation and entrapment of hydrothermal fluids and consequent ore mineralization in the study area.

Fusion of Exploratory/Information Layers
Fuzzy logic and index overlap methods were used for fusing the exploratory/information layers derived from remote sensing, magnetometry, and geological datasets. The fuzzy logic model was operated to pinpoint the potential zones in the study area. The maps were fused in pairs or in different stages to obtain the final potential maps. To achieve this goal, the importance of ore mineralization factors must be determined. All

Fusion of Exploratory/Information Layers
Fuzzy logic and index overlap methods were used for fusing the exploratory/information layers derived from remote sensing, magnetometry, and geological datasets. The fuzzy logic model was operated to pinpoint the potential zones in the study area. The maps were fused in pairs or in different stages to obtain the final potential maps. To achieve this goal, the importance of ore mineralization factors must be determined. All exploratory/information layers were weighted in this regard [6,45]. The types of layers and created classes, along with the weight of the classes, are presented in Table 2. Subsequently, fuzzy sum and gamma operators were used. In this study, the geological, alteration, and intrusive data layers were fused using the fuzzy summative operator. The fuzzy gamma operator was applied to the fault density layer. Gamma values of 0.8 and 0.9 were utilized to ensure model quality. The geological, remote sensing, and aerial magnetometry were fused using the algebraic addition fuzzy operator, and then, the results were merged with two gamma values of 0.8 and 0.9 with the fault density map. Figure 10 shows the potential map for the study area derived from the fuzzy logic model. Three high potential zones (A-C) were identified, which were typically placed in the central and southwestern parts of the study area ( Figure 10). Most of the documented mineralization precincts in the study area are located in high potential zones.
The multi-class overlap index method was also applied to the exploratory/information layers. To analyze the status of the lineaments and the role of these layers in the mineralization, they were classified into six groups according to the attributes in Table 3. Intrusive masses (detected using aerial magnetometry) were classified into four classes and scored (Table 4). Furthermore, the maps obtained from remote sensing (alteration minerals) were prepared in two classes. The scoring of the layers was implemented based on previous studies [36,40] and expert's opinion. The areas containing iron oxide/hydroxides and hydroxyl-bearing (Al-OH) alteration and carbonates were given a score of 10, and the other points were given a score of zero. The characteristics of rock units and their relationship with ore mineralization in the study area were also assumed. Table 5 shows geological layer classification for integration by multi-class index overlap method.  The multi-class overlap index method was also applied to the exploratory/information layers. To analyze the status of the lineaments and the role of these layers in the mineralization, they were classified into six groups according to the attributes in Table 3. Intrusive masses (detected using aerial magnetometry) were classified into four classes and scored (Table 4). Furthermore, the maps obtained from remote sensing (alteration minerals) were prepared in two classes. The scoring of the layers was implemented based on previous studies [36,40] and expert's opinion. The areas containing iron oxide/hydroxides and hydroxyl-bearing (Al-OH) alteration and carbonates were given a score of 10, and the other points were given a score of zero. The characteristics of rock units and their relationship with ore mineralization in the study area were also assumed. Table 5 shows geological layer classification for integration by multi-class index overlap method.   For implementing the multi-class overlap index method, the iron oxide/hydroxides layer was assigned the highest weight of 8 (i.e., 8 out of 10), the importance of the fault layer was weighted at 6, and the layer of intrusive masses was given the weight of 5, the geological layer and the type of principal rocks (rhyolite granite, etc.) were weighted at 3, respectively. Finally, the hydroxyl-bearing (Al-OH) and carbonates alteration layer was weighted at 1. After weighting the exploratory layers, these layers were fused using the multi-class index overlap method, and a potential map was generated for the study area ( Figure 11). In this map, six high potential zones (A-E) were identified, which were mostly located in the central, southwestern, and northeastern part of the study area ( Figure 11). Several documented mineralization points were delineated in the high potential zones that were identified using the multi-class index overlap method. The analysis of mineral potential maps derived from the fuzzy logic model and the multi-class index overlap method shows that the central and southwestern parts might be high prospective zones for the next stage of the mineral exploration campaign. The mineral potential map derived from the multi-class index overlap method shows more (high to moderate) mineral potential zones in the study area compared to the fuzzy logic model outputs (see Figures 10 and 11). The zones of C and E were only found as high potential zones in the multi-class index overlap method output (see Figure 11). It might show that the mineral potential map derived from the multi-class index overlap method is more sensitive to mineral exploration factors by considering uncertainty and predicting intervention.

Geology and Fieldwork Results
The general trend of mineralization in the region is controlled by geological structural features. Mineralization in the study area consists of copper, iron, barite veins, kaolinite, and lead, which are locally mined. Mineralization zones and hydrothermal alterations in the study area were checked during fieldwork. One of the potential zones contains iron mineralization, which is hosted in a volcano-sedimentary unit (Rizo Formation). Rhyolitic tuffs comprise the alterations zones associated with red iron oxides and feldspars. The altered rhyolites are red in areas with high amount of iron oxides, whereas they are white in areas with high content of feldspars ( Figure 12A). Rhyolitic tuffs are intruded by andesitic dykes with a thickness of 5 to 10 m.
Copper mineralization veins were found with dolomitic units in the central and eastern parts of the study area. These copper veins are formed along the fractures and faults of the host rock. Malachite and chalcopyrite, which were found as the surface expression of copper mineralization, belong to the oxidation zone ( Figure 12B,C). Therefore, the possibility of an enriched mineral zone or supergene zone at the lower level of the veins is feasible. Iron oxide veins and lenses were found in the rhyolites. The thickness of these lenses reached 10 m, and their length reached 30 m. Iron mineralization also replaced dolomites, and it was found as interlayer within rhyolites and acidic tuffs ( Figure 12D). Both the presence of alunite and jarosite in the study area and the secondary iron oxides (goethite and limonite) indicate acid leaching ( Figure 12E). Based on the drilling boreholes up to an average depth of 13 m, there are argillic and propylitic alterations, followed by phyllic and potassic alterations and metal mineralization zones. Pyrite and magnetite mineralization are observable in the cores ( Figure 12F).

Geology and Fieldwork Results
The general trend of mineralization in the region is controlled by geological structural features. Mineralization in the study area consists of copper, iron, barite veins, kaolinite, and lead, which are locally mined. Mineralization zones and hydrothermal alterations in the study area were checked during fieldwork. One of the potential zones contains iron mineralization, which is hosted in a volcano-sedimentary unit (Rizo Formation). Rhyolitic tuffs comprise the alterations zones associated with red iron oxides and feldspars. The altered rhyolites are red in areas with high amount of iron oxides, whereas they are white in areas with high content of feldspars ( Figure 12A). Rhyolitic tuffs are intruded by  Copper mineralization veins were found with dolomitic units in the central and eastern parts of the study area. These copper veins are formed along the fractures and faults of the host rock. Malachite and chalcopyrite, which were found as the surface expression of copper mineralization, belong to the oxidation zone ( Figure 12B,C). Therefore, the possibility of an enriched mineral zone or supergene zone at the lower level of the veins is feasible. Iron oxide veins and lenses were found in the rhyolites. The thickness of these lenses reached 10 m, and their length reached 30 m. Iron mineralization also replaced dolomites, and it was found as interlayer within rhyolites and acidic tuffs ( Figure 12D). Both the presence of alunite and jarosite in the study area and the secondary iron oxides (goethite and limonite) indicate acid leaching ( Figure 12E). Based on the drilling boreholes up Moreover, lenticular bodies, thin and scattered veins of barite, are hosted in the Paleozoic sedimentary rocks. The thickness of the lenticular bodies and veins is up to several centimeters and their length in some cases reaches 30 m ( Figure 13A). Barite veins, filled along the main fractures and faults, are found within carbonate rocks ( Figure 13B). Surface expression of the argillic alteration zone, as industrial soil mine (kaolinite) and sericite-silica veins associated with the mineralization of pyrite and galena, were found in many parts of the study area ( Figure 13C,D).
Moreover, lenticular bodies, thin and scattered veins of barite, are hosted in the Paleozoic sedimentary rocks. The thickness of the lenticular bodies and veins is up to several centimeters and their length in some cases reaches 30 m ( Figure 13A). Barite veins, filled along the main fractures and faults, are found within carbonate rocks ( Figure 13B). Surface expression of the argillic alteration zone, as industrial soil mine (kaolinite) and sericite-silica veins associated with the mineralization of pyrite and galena, were found in many parts of the study area ( Figure 13C,D). A confusion matrix (error matrix) and kappa coefficient [46][47][48] were calculated for the mineral potential maps versus field data. The comparison of the mineral potential maps with field data using a confusion matrix approach and the kappa coefficient shows a very good match, which indicates the overall accuracy of 73% and the kappa coefficient of 0.60, respectively ( Table 6). The producer's accuracy (omission error) indicates the number of correct classified samples of a class (X) that is divided by the total number of reference samples of class (X) (column total). The resulting percentage accuracy indicates the probability that a reference (ground) sample will be correctly classified. The user's accuracy (commission error) shows the number of correctly classified samples of class (X) are divided by the total number of samples that were classified in class (X) (row total). The A confusion matrix (error matrix) and kappa coefficient [46][47][48] were calculated for the mineral potential maps versus field data. The comparison of the mineral potential maps with field data using a confusion matrix approach and the kappa coefficient shows a very good match, which indicates the overall accuracy of 73% and the kappa coefficient of 0.60, respectively ( Table 6). The producer's accuracy (omission error) indicates the number of correct classified samples of a class (X) that is divided by the total number of reference samples of class (X) (column total). The resulting percentage accuracy indicates the probability that a reference (ground) sample will be correctly classified. The user's accuracy (commission error) shows the number of correctly classified samples of class (X) are divided by the total number of samples that were classified in class (X) (row total). The resulting percentage accuracy is indicative of the probability that a sample from the classified image/map represents that class on the ground [47]. The highest producer's accuracy (80%) was obtained for the high potential class and the highest user's accuracy (93%) was achieved for the low potential class. It shows less possibility of confusion between the high potential class and low potential class. For the moderate potential class, the producer's accuracy, and user's accuracy both were 65% (Table 6). It might indicate that the moderate potential class has the capability to be confused with the high or low potential classes. Looking at the user's accuracy for all classes shows that the possibility of confusion between the high potential (66%) and moderate potential (65%) classes is more feasible compared to the low potential class (93%) ( Table 6).

Discussion
The exploration of deep-seated/concealed polymetallic mineral deposits is very challenging due to multiple factors controlling ore mineralization and conditions and mechanisms of hydrothermal alteration systems [49]. It emphasizes the fusion of different sources of geoscience datasets such as geophysical data, geochemical data, remote sensing data, and geological data for polymetallic mineral prospecting. Numerous studies used geoscience information mutually or separately (e.g., geological, geophysical, geochemical, and remote sensing datasets) for polymetallic mineral prospectivity mapping [36,[49][50][51][52][53]. This investigation used remote sensing data to extract hydrothermal alteration zones, aeromagnetic data for identifying deep intrusive masses and hidden fault systems, and geological maps and filed data to show mineral occurrences and mine locations, faults and fractures, and potential host-rock units.
In this analysis, two knowledge-driven spatial mathematical models, namely, the fuzzy logic method and the multi-class index overlap method were implemented to fuse twenty information layers derived from remote sensing, geophysical, and geological datasets for the study area. Iron oxide/hydroxides, hydroxyl-bearing (Al-OH) alteration, and carbonates were mapped using BR, PCA, and LS-Fit techniques applied to Sentinel-2 and Landsat-7 remote sensing data (see . Iron oxide/hydroxides were typically mapped in igneous rocks, sandstone, conglomerate, limestone, and marl. Hydroxyl-bearing (Al-OH) alteration and carbonates were mainly mapped associated with igneous rocks, limestone, marl, and sandstone. Sedimentary rocks such as quaternary unit, limestone, marl, and sandstone can be incorrectly mapped as hydrothermal alteration minerals due to high amounts of detrital clays and carbonate minerals [54,55]. Thus, the detected alteration zones associated with sedimentary units must be omitted and not considered for mineral potential mapping during the analysis.
Deep intrusive masses and fault and fracture systems play important roles in polymetallic ore mineralization [49]. Deep intrusive masses and hidden fault systems in the study area were identified using the trend removal filter, vertical derivative filter, analytical signal filter, TDX filter, and horizontal derivative filters executed to aeromagnetic data (see . Magnetic anomalous zones of deep intrusive masses were detected in the northwestern, central, northeastern, and southeastern parts, which are mostly covered with the Quaternary deposits, sandstone, and marl (see Figures 1B, 7A-C and 8A,B). Deep fault systems in the study area show major directions of N-S, NW-SE, NE-SW, having some intersection with short faults running NE-SW, W-E, N-S (see Figure 8C). The intersection zones might have high potential for polymetallic ore mineralization, which can act significant roles in trapping or conducting (pathways) hydrothermal ore fluid. The spatial distribution of the intersection zones is typically detected in the central, southwestern, and northeastern part of the study area (see Figure 8C).
For fusing the exploratory layers derived from remote sensing, magnetometry, and geological datasets, the fuzzy logic and multi-class index overlap methods were implemented. Potential maps for the study area derived from the methods show some prospective zones, which are typically located in the central, southwestern, and northeastern sectors (see Figures 10 and 11). A field survey of the study area indicated the presence of copper, iron, barite veins, kaolinite, and lead, which are mainly associated with hydrothermal alterations and particularly hosted within igneous units (rhyolite, granite, diabase) and carbonate units. The calculation of the confusion matrix approach and kappa coefficient for the mineral potential maps versus field data confirms a very good match for the results with an overall accuracy of 73% and the kappa coefficient of 0.60. Therefore, the potential zones identified in this study can be considered for consecutive detailed field campaigns and drilling programs.

Conclusions
In this study, multispectral remote sensing data (e.g., Sentinel-2 and Landsat-7), airborne magnetic geophysical data, and geological data were fused using knowledge-driven spatial mathematical models (i.e., fuzzy logic and multi-class index overlap models) for identifying potential zones of polymetallic mineralization in the Chakchak region, Yazd province, Iran. Subsequently, two mineral potential maps were produced for the Chakchak region. Several prospective/anomalous zones for polymetallic mineralization were pinpointed in the central, southwestern, and northeastern parts of the study area. In this analysis, only airborne magnetic geophysical data was available for the study region; consequently, the data were used to identify superficial and deep intrusive masses and deep and hidden faults. However, airborne magnetic data with higher density (centralized measured network) will be valuable for accurate modeling and identification of magnetic anomalies related to iron deposits in the study area. Additionally, more accurate results will be obtained if gravimetric data is available to fuse with the other exploratory/information layers for a regional scale mapping. For the anomalous zones that were identified in this study, it is recommended to investigate using ground-based electromagnetic methods such as self-potential (SP), induced polarization, and resistivity geoelectric methods (IP-Rs) for electrically conductive minerals to explore sulfide mineralization. Furthermore, gravity surveys and ground magnetics could be executed to explore iron deposits and narrow down and confirm airborne magnetic results. The high potential zones can be considered for consecutive detailed field campaigns, geochemical surveys, drilling programs, and laboratory analysis. In conclusion, the fusion of different sources of geoscience datasets using knowledge-driven spatial mathematical models has great capability to identify the potential zones of polymetallic mineralization at a regional scale. This approach is applicable, robust, and affordable for mining companies to explore other regions having similar geologic conditions for polymetallic ore mineralization.