Introduction

Gypsum (CaSO4⋅2H20), Bassanite (CaSO4⋅0.5H20), and anhydrite (CaSO4) are common Ca-sulfates minerals on the Earth (Bishop et al. 2014). Ca-sulfates (gypsum, bassanite, and anhydrite) typically formed under different temperature conditions in arid evaporation environments (Billo 1987), and they are very important not only for paleo-environmental and paleo-climatic studies but also for industrial applications (components of gypsum boards and cement (Singh and Garg 1995)). Ca-sulfates rocks also play an important role in oil and natural gas exploration in Iran, as they may form the trap where hydrocarbons are accumulated (Gürbüz 2019). In addition, gypsum and bassanite have even been identified in several regions on Mars by both the Observatoire pour la Minéralogie, l´Eau, les Glaces et l´Activité (OMEGA) onboard Mars Express and the Compact Reconnaissance Imaging Spectrometer for Mars (CRISM) onboard Mars Reconnaissance Orbiter (MRO) (Wray et al. 2010, 2011; Sefton-Nash et al. 2012; Weitz et al. 2015). Thus, remote detection of Ca-sulfates is essential to mineral exploration, paleo-climatic studies, and minerals mapping on Mars’ surface.

The spectral properties of Ca-sulfates are the basis of remote detection of gypsum. Previous studies have presented Ca-sulfates’ spectral properties in the visible and near-infrared (VNIR), short-wave infrared (SWIR), mid-infrared (MIR), and thermal infrared (TIR) regions (Hunt 1971; Crowley 1990, 1991; Cloutis et al. 2006; Lane 2007; Harrison 2012; Bishop et al. 2014). Hydrous Ca-sulfates (Gypsum, bassanite) exhibit diagnostic absorption features on VNIR and SWIR bands (a characteristic triplet at 1.446, 1.490, and 1.538 µm; 1.75 µm; 1.9 µm; a triplet at 2.178, 2.217, and 2.268 µm for gypsum) (Howari 2002; Bishop et al. 2014), attributed to vibrations of hydrogen-bonded structural water molecules (Crowley 1991). While, these features are not observed for anhydrous Ca-sulfates (anhydrite) and other evaporite minerals (halite, thenardite, etc.) (Hunt 1971; Bishop et al. 2014). All Ca-sulfates show strong absorption near 4.2–5.0 μm, attributed to overtones and the ν3 SO4 asymmetric stretching vibration, one of the S–O-stretching modes (Bishop and Murad 2005; Bishop et al. 2014). And for TIR, Ca-sulfates exhibit diagnostic spectral emissivity features at 1050–1250, 1000, 500–700, and 400–500 cm–1(Lane 2007). However, the existing satellite sensors do not cover the MIR region, and the rather low spatial resolution of TIR data has a significant impact on mapping accuracy (Khan et al. 2020). Spectral properties on VNIR and SWIR bands are therefore essential to identify hydrous Ca-sulfates’ minerals using multi-spectral and hyper-spectral data.

Different data-driven methods, such as data fusion (Kavak 2005), Decorrelation Stretching (Chapman et al. 1989; Kavak 2005), Linear Mixture Modeling (Bryant 1996), Support Vector Machines (SVM), Spectral Angle Mapper (SAM) (Soltaninejad et al. 2018), Mixture Tuned Match Filtering (MTMF) (Pour et al. 2018), Spectral Feature Fitting (SFF), and Spectral Mixture Analysis (SMA)(Robertson et al. 2016), were applied to multi-spectral and hyper-spectral images for gypsum identification. However, the spectral properties of gypsum may not be adequately considered in these studies. No targeted enhancement based on gypsum spectral properties was tried in these data-driven methods. And the accuracy of these methods relied on the correct selection of training samples, reference spectrum, or spectral endmember, which was always complicated and challenging. In addition, knowledge-based approaches were applied for gypsum mapping, as well. Feature-oriented principle component selection (FPCS) was applied to ASTER data for gypsum mapping in Turkey (Öztan and Süzen 2011; Özyavaş 2016; Gürbüz 2019) and Pakistan (Khan et al. 2020). Same bands (4,6,8,9) were chosen based on the spectral features of gypsum in these studies, but different principle components were utilized for gypsum detection. Least-squares spectral band-fitting algorithm was applied to Airborne Visible/Infrared Imaging Spectrometer (AVIRIS) data for playa evaporite mapping in Death Valley, California, and multiple spectral features of gypsum were chosen for band-fitting analysis (Crowley 1993). Gypsum spectral properties were considered in the FPCS and least-squares spectral band-fitting algorithm; however, the problem of these methods is time-consuming and unstable.

Spectral indices or band ratios are easier to transfer among sensors and can be more robust (Milewski et al. 2019). Various gypsum spectral indices, applying ASTER SWIR bands, (band4 + band8)/band6 (Öztan and Süzen 2011; Özyavaş 2016), band8/band6 (Gürbüz 2019), band4/band9 (Öztan and Süzen 2011), band 4 + band 8/band 6 (Khan et al. 2020), and band 4/ (band 6 + band 9) (Salati et al. 2019), were proposed in recent studies. Most of these spectral indices were designed based on absorption features of gypsum at 2.21 μm (ASTER band6). However, the 2.21 μm absorption feature of gypsum can be confused with the 2.2 μm hydroxyl-bearing minerals absorption (Crowley 1993; Notesco et al. 2016; Milewski et al. 2019). Sulfate index (SI) (Öztan and Süzen 2011) or gypsum index (GI) (Ninomiya and Fu 2019) for ASTER TIR bands (band 10 \(\times \) band 12)/(band 11 \(\times \) band 11) was also proposed, but the ASTER TIR index is highly vulnerable to mixed pixel effect for the low spatial resolution. Besides, a Normalized Differenced Gypsum Ratio (NDGI) was proposed for hyper-spectral data, based on the 1.75 μm gypsum absorption (Milewski et al. 2019). Another hyper-spectral gypsum index was designed for hyper-spectral longwave infrared (LWIR) data, based on gypsum emissivity absorption feature at 8.63 μm (Notesco et al. 2015). However, the application of the hyper-spectral gypsum index is limited, attributing to few hyper-spectral remote-sensing data sources.

All above image processing methods for gypsum mapping have ignored the characteristic 2.26 μm (ASTER band 7) gypsum absorption, not overlapping with most hydroxyl-bearing minerals and other major rock-forming minerals. Although the ankerite also shows a 2.26 μm absorption feature (Mars and Rowan 2011), it is easy to distinguish ankerite and gypsum according to the significant absorption feature difference at 2.21 μm or 2.34 μm between the two minerals.

Previous studies held different opinions about the seasonal variation of image reflectance of rocks. The reflectance of Granite outcrops, covered with lichen and algae, was reported to keep consistent all year across the southwest Australian floristic region (Alibegovic et al. 2015). On the contrary, the image reflectance of rocks was considered variable with seasons of the year, because the physical and chemical properties of rocks may change under heavy precipitation and heat depending on their crystal structure (Yao et al. 2017). Namibian seasonal salt pan crust dynamics were monitored by a dense time-series consisting of 77 Landsat 8 cloud-free scenes, and climate controls on crust formation were discussed (Milewski et al. 2020). Moreover, the performance of ASTER mineral spectral indices was proven to be varying according to the data acquisition season (Hewson and Cudahy 2011). In addition, the mutual transformation between Ca-sulfates (Fig. 1) under laboratory conditions and natural conditions has been presented. Experiments of Harrison (2012) showed that completely dehydrated gypsum (anhydrite) rehydrated to form bassanite at ambient conditions within 1 day. Meanwhile, particulate anhydrite may also be slowly hydrated to form gypsum in a wet environment at low temperatures (Sievert et al. 2005). Under natural conditions, gypsum can be dehydrated to form bassanite and anhydrite under intense solar heating in the arid regions (Warren 2006; Rezaee et al. 2016). Moreover, it was reported that bassanite is present in the dry summer season in Kuwait and hydrates to gypsum in the wet winter months (Gunatilaka et al. 1985). Furthermore, the SWIR absorption feature of Ca-sulfates is closely related to structural water molecules contents (Crowley 1991) (Fig. 2). Consequently, the reflectance of Ca-sulfates may change according to the seasonal change of Ca-sulfates facies. However, to our knowledge, no research has attempted to explore the seasonal variation in image reflectance of Ca-sulfates formations in arid regions.

Fig. 1
figure 1

The mutual transformation between gypsum and bassanite (Gunatilaka et al. 1985)

Full size image
Fig. 2
figure 2

Lab reflectance of gypsum, bassanite, and anhydrite in SWIR region (1.5–2.9 µm) from USGS mineral spectral library. Gypsum sample is HS333.3B, bassanite is GDS145, and anhydrite is GDS42

Full size image

The objectives of this paper are (1) to assess the SWIR absorption feature for gypsum mapping based on published spectral library and ASTER image reflectance, and to design new gypsum indices (the GI1 and GI2), eliminating the interference of hydroxyl-containing minerals for gypsum mapping; (2) to analyze seasonal reflectance variability of gypsum formation pixels, and to access the accuracy of spectral indices for multi-seasonal ASTER data.

Regional setting

The study area is located in the Zagros fold-thrust belt, south of Iran (Salati et al. 2019). Climate appears to be a critical factor affecting the distribution of evaporite mineral facies (Gunatilaka et al. 1985). The average rainfall in the study area is about 190 mm/y, falling mostly between November and April. Summer average temperature is about 37 ℃, the ground surface temperature is up to 60 ℃.

The gypsum formations or gypsum-bearing formations in the study area are Gachsaran Formation (Miocene) and the Hormoz Complex (Precambrian). The Gachsaran Formation contains mainly evaporites, and it also contains marls, limestones, and shales (Bahroudi and Koyi 2004). Traditionally, the Kazerun–Qatar fault divides the Gachsaran basin into two different segments: (1) Fars province in the SE and (2) a Northwestern basin (Stocklin 1968). In Fars Province, Gachsaran Formation consists of three members from bottom to top, namely the Chehel Member (anhydrite or gypsum), the Champeh Member (gypseous marl, chalky-gypseous limestone or dolomite, and nodular to massive gypsum), and the Mol Member (gypseous marl and gypseous limestone) (James and Wynd 1965). The Hormoz Complex appears as Salt plugs on the ground (Bosák et al. 1998). The Bam plug locates in the west of the study area. Petrographic composition of the plug is variable with the occurrence of clastic sedimentary, magmatogenic rock types, and chemogenic, including shales, siltstones, limestones, dolostones, gabbro, diorite, quartz diorite, and propyllitized andesite (Bosák et al. 1998). Gypsum is common as layers and interbeds in sediments, and white gypsum in a tremendous amount covers terrace sediments at the eastern plug margin (Bosák et al. 1998).

Figure 3b shows a simplified geological map of the study area, which was interpreted from a high-resolution Google image with reference to the LAMAZAN and ROYDAR-E-SOFLA geological maps (1/100000) (Moeini et al. 2004; Aghababaei and Motamedi 2006). Based on previous studies (James and Wynd 1965; Stocklin 1968; Bosák et al. 1998; Bahroudi and Koyi 2004; Mortazavi et al. 2017), the strata of the study area were divided into three categories: gypsum and anhydrite formations (GF), gypsum-bearing formations (BGF), and non-gypsum formations (NGF) (Table 1).

Fig. 3
figure 3

modified from Moeini et al. (2004) and Aghababaei and Motamedi (2006), horizontal lines in the polygons represent gypsum and anhydrite formations, black dot pattern represents gypsum-bearing formations

The ASTER image and simplified geological map of the study area. a The ASTER false color composite image of the study area (bands 3, 2, and 1 assigned in red, green, and blue, respectively); b simplified geological map of the study area

Full size image
Table 1 Formations/members and lithologic characters of the study area
Full size table

Methodology

Six ASTER level 1B (L1B) scenes were acquired in dry summer and wet winter over the study area (Table 2). Digital number (DN) values of ASTER L1B were converted to calibrated spectral radiance at-sensor. Atmospheric correction was carried out, respectively, to VNIR bands and SWIR bands by FLAASH model (ENVI 5.3). Thermal Atm Correction Module (ENVI 5.3) was used to approximate and remove the atmospheric contributions from TIR bands radiance data. The datasets and methods used in this study are described in a flowchart (Fig. 4).

Table 2 Acquisition data, solar elevation, and solar azimuth of multi-seasonal ASTER data
Full size table
Fig. 4
figure 4

Flowchart of the datasets and methods applied

Full size image

Acquiring training samples

Spectral characteristics of “pure” pixels were usually used to identify land-cover types that affect target objects extraction accuracy (Feyisa et al. 2014; Fisher et al. 2016) while attempting to design an index to discriminate target objects and other objects accurately. However, different from water or plants, most of the rocks show spectral variability (Murphy et al. 2012; Zhang and Li 2014) (the difference in whole spectral shape, the variations in spectral absorption position, and depth) due to the difference of minerals proportions, chemistry, structure, grain size of same rock types. The spectral variability makes it difficult to obtain “pure” pixels for rocks and minerals. Here, we used approximate “pure” pixels as training samples to evaluate the capacity of ASTER bands to distinguish between gypsum formations and other formations.

The methods used to acquire approximate “pure” pixels (training samples) of geologic formations include Minimum Noise Fraction (MNF) Transform and the Sequential Maximum Angle Convex Cone (SMACC). MNF transform was used to determine the inherent dimensionality of image data, segregate noise in the data, and reduce the computational requirements for subsequent processing (Luo et al. 2016). SMACC is a tool to find spectral endmembers and their abundances. In this study, we applied MNF transform to ASTER NIR and SWIR reflectance (ASTER 1–9 bands, 2004-11-22). Then, the first 6 output channels of the MNF transform were applied to SMACC tool, and 20 spectral endmembers and their abundance images were acquired. A simplified geological map and Google image were used as a reference to identify spectral endmembers and their abundance images for each geologic formation. Thresholds of 0.60–0.75 were used to extract training samples from the abundance images. The thresholds were decided in a process of trial and error with the reference of a simplified geological map and Google image. As shown in Fig. 5, training samples of alluvial terraces and recent deposit are divided into three groups, based on their internal spectral variability and the difference of sediment source. The sediment sources of Qhal1, Qhal2, and Qhal3 are mainly clastic rocks (BK), clastic and evaporites rocks (Aj-Mn and GS), and Hormoz Complex, respectively, according to their surrounding terrain and lithology distribution. Except that a few training samples (about 5%) from different geologic formations are confused for their similar lithology, most training samples are distributed within their geologic formation regions.

Fig. 5
figure 5

Distributions of training samples of geologic formations

Full size image

Spectral analysis

To assess the 2.26 µm absorption for gypsum mapping, and to develop an index that maximizes the differences of image reflectance values between gypsum and other minerals, lab and ASTER SWIR reflectance of gypsum and other minerals were compared.

Lab reflectance comparison between gypsum and other minerals

We compared lab reflectance between gypsum and several minerals with diagnostic absorption features in SWIR region. As shown in Fig. 6a, absorption positions of gypsum in SWIR region (1.4–2.8 μm) are at 1.49 µm, 1.75 µm, 1.9 µm, 2.21 µm, and 2.26 µm. And it can be identified that the 2.21 µm absorption position of gypsum overlaps with hydroxy-bearing minerals (illite, kaolinite, montmorillonite, chlorite, and muscovite) and weathered feldspar minerals, of which 2200 nm absorption features were considered to be a starting alteration to hydroxy-bearing clays (Graham et al. 1970; Hecker et al. 2010). The overlapping of absorption position is the main resource of commission errors for the gypsum indices, designed based on 2.21 μm (ASTER band6) absorption of gypsum. Moreover, the 1.9 µm absorption position of gypsum also overlaps with other hydrated minerals (illite, kaolinite, and montmorillonite), and 1.49 µm absorption position is outside the atmospheric window. To our knowledge, there is no multi-spectral satellites sensor covering 1.75 μm. Thus, it should be noted that the 2.26 µm gypsum absorption does not overlap with other minerals, and probably can be detected by ASTER data (center wavelength of ASTER band7 is exactly at 2.26 µm). However, ASTER band7 has not been used to calculate a gypsum spectral index.

Fig. 6
figure 6

Lab and ASTER reflectance spectra of minerals and formation members. a Mineral spectra from USGS mineral spectral library. Gypsum sample is HS333.3B, illite is GDS4, kaolinite is CM3, montmorillonite is CM20, muscovite is GDS107, oligoclase is HS110.3B, chlorite is HS179.3B, and albite is HS143.3B. b ASTER reflectance of geological formation spectral members

Full size image

ASTER reflectance comparison between GF and other formations

Spectral endmembers reflectance of GF, BGF, and NGF was compared in Fig. 6b. Reflectance curves of GF spectral members (Chm and HS2) show lower reflectance values at 2.21 μm and 2.26 μm in SWIR region, indicating that 2.26 μm reflectance absorption of gypsum can be detected by ASTER image. For BGF, spectral member reflectance of Mln-Cpm shows absorption features at 2.26 μm and 2.33 μm, attributing to Ca-sulfates (gypsum) and carbonates (marl, limestone). The absorption positions of GF and HS1 are similar, but their slopes between 1.56 and 2.21 μm are different. Most NGF (As-Ja, Grm, Aj-Mn) spectral members show absorption features at 2.33 μm for their lithology of carbonates. It can be noted from Fig. 6b that spectral member reflectance of Qhal1 and Qhal3 also has absorption features at 2.21 μm, attributing to clay minerals.

Calculating new gypsum indices

The primary purpose of designing new gypsum indices is to exclude the potential interference of hydroxyl-bearing minerals under the premise of mapping accuracy assurance.

The ASTER SWIR reflectance (4–9 bands) distributions of training samples of GF, BGF, and NGF are shown in Fig. 7. It can be identified that the most significant difference of reflectance distributions between GF (Chm) and other geologic formations is observed in band6 (2.21 μm) and band7 (2.26 μm). However, the reflectance distributions of Chm severely overlap with HS1 and Qhal3 in both two bands. The lithology of HS1 is complicated, except for gypsum, weathered gabbro, diorite, and quartz diorite (illite, montmorillonite, and kaolinite are the main secondary mineral (Pédro 1997; Younis et al. 1997; Deepthy and Balakrishnan 2005)), and propyllitized andesite (altered minerals are mainly chlorite and epidote) also shows absorption feature at 2.21 μm. Alluvium most probably contains a large number of clay minerals (Notesco et al. 2015), showing 2.21 μm absorption. Meanwhile, we found that band4 performed well at discriminating the three geologic formations, and Chm’s reflectance in band4 is greater than Qhal3 and HS1.

Fig. 7
figure 7

Reflectance distributions of training samples of geologic formations and members. Each box plot shows the location of the 10th, 25th, 50th, 75th, and 90th percentiles using horizontal lines, and the circles are 5th and 95th percentiles. The boxes of gypsum and anhydrite formations (Chm) are marked as green in color

Full size image

Four SWIR bands of ASTER were used in developing the new gypsum indices to increase the value contrast between gypsum pixels and other minerals pixels, particularly to distinguish gypsum and hydroxyl-containing minerals. New gypsum indices specified as follows:

$$ {\text{GI}}^{{1}} = \, \left( {\rho \, \left( {{\text{band4}}} \right) \times \rho \, \left( {{\text{band8}}} \right)} \right)/ \, \left( {\rho \, \left( {{\text{band6}}} \right) \times \rho \, \left( {{\text{band7}}} \right)} \right), $$
(1)
$$ {\text{GI}}^{{2}} = \, \left( {\rho \, \left( {{\text{band4}}} \right) \times \rho \, \left( {{\text{band8}}} \right)} \right)/ \, \left( {\rho \, \left( {{\text{band7}}} \right) \times \rho \, \left( {{\text{band7}}} \right)} \right), $$
(2)

where ρ is the reflectance value of ASTER bands: band4 (1.66 μm), band6 (2.21 μm), band7 (2.26 μm), and band8 (2.34 μm). As shown in Fig. 6, band4 (1.66 μm) and band8 (2.34 μm) are reflectance shoulders of gypsum, while gypsum shows absorption features at 2.21 μm (band6) and 2.26 μm (band7). Therefore, multiplying ρ (band4) and ρ (band8), while at the same time dividing by ρ (band6) and ρ (band7) yields relatively large positive values for gypsum pixels compared to other mineral pixels.

Besides, different from 2.21 μm, the absorption position of gypsum at 2.26 μm (band7) does not overlap with clay, weathered feldspar, and muscovite minerals (Fig. 6). Therefore, dividing by ρ (band7, 2.26 μm) enhances the separability between gypsum minerals and hydroxyl-bearing minerals, theoretically, increasing values of gypsum pixels. And, as shown in Fig. 7, the reflectance of Chm (gypsum) in band4 is greater than Qhal3 and HS1. Thus, multiplying ρ (band4, 1.66 μm) may increase the difference of the pixel values between HS1, Qhal3, and Chm, further increasing values of Chm (gypsum) pixels.

Receiver operator characteristic (ROC) curves

The performance of the different gypsum spectral indices was evaluated using receiver operator characteristic (ROC) curves with ground truth information. ROC curves plot the true-positive rate (TPR) against the false-positive rate (FPR) for a range of threshold values (Fawcett 2006; Fisher et al. 2016). When an ROC curve of a spectral index is the closest to the top-left corner, it represents that the spectral index achieves more accurate mapping results than other indices with a limited range of thresholds. Besides, the proportion of area under the curve also indicates the accuracy of spectral indices. A proportion of 100% means a perfect spectral index; a proportion of 50% means a random spectral index (Fisher et al. 2016). When the area under ROC curve of a spectral index is the largest, it represents that the spectral index achieves more accurate mapping results than other indices with conservative or liberal thresholds. Ground truth regions (pixels) of GF, BGF, and NGF were determined from the simplified geological map of the study area (Fig. 1b).

The optimum threshold of each spectral index, which achieved the lowest combined error (100 − TPR + FPR) for GF range (pixels) (Fisher et al. 2016), was also applied to map GF for multi-seasonal ASTER data. And, the results were compared with the reference of the simplified geological map (Fig. 1b).

Results

Analyzing seasonal reflectance variability of GF

To explore whether Ca-sulfates facies of GF change between dry summer and wet winter in the study area, 2.21 μm (ASTER band6) and 2.26 μm (ASTER band7) reflectance of GF spectral member and training samples were compared for winter and summer multi-seasonal ASTER data.

Seasonal reflectance variability of Chm spectral member and training samples are presented in Figs. 8 and 9. As shown in Fig. 8, Chm spectral member reflectance curves shapes, absorption positions, and absorption depth of summer ASTER data (2001-7-25, 2002-7-28, 2003-7-15) are similar. And the same phenomenon is observed for winter ASTER data (2000-11-18, 2004-3-11, and 2004-11-22). Nevertheless, the differences between winter and summer reflectance of Chm spectral member are as follows: (1) the summer ASTER reflectance is higher than winter overall, attributing to larger solar elevation and stronger solar radiation in summer; (2) 2.26 μm (band7) reflectance shows the most significant change, leading to different slopes between 2.21 μm (band6) and 2.26 μm (band7) for summer and winter ASTER data (Fig. 8). These results were further confirmed by comparing reflectance distributions of Chm training samples between summer and winter ASTER data. As shown in Fig. 9, the summer ASTER band6 (2.21 μm) and band7 (2.26 μm) reflectance of Chm training samples are higher than winter overall; furthermore, the reflectance change magnitude of band7 (2.26 μm) is greater than band6 (2.21 μm).

Fig. 8
figure 8

Multi-seasonal ASTER reflectance of Chehel (Chm) spectra member. The study area is dry summer between May and October, and wet winter between November and April

Full size image
Fig. 9
figure 9

Reflectance distributions of training samples of Chehel (Chm) for multi-seasonal ASTER data. Each box plot shows the location of the 10th, 25th, 50th, 75th, and 90th percentiles using horizontal lines, and the circles are 5th and 95th percentiles. The boxes of dry summer ASTER data are marked as red, and blue for wet winter ASTER data

Full size image

Seasonal ASTER reflectance variability of GF pixels was confirmed, attributing to the seasonal mutual transformation between Ca-sulfates, possibly. Furthermore, seasonal variability in absorption depth of band6 (2.21 μm) and band7 (2.26 μm), and the slope of reflectance between them is bound to affect the mapping accuracy of gypsum indices, calculated by band6 and band7. Thus, we evaluated the stability of gypsum indices for winter, summer, and multi-seasonal ASTER data.

Evaluation of gypsum indices for winter ASTER data

The ASTER image, acquired on November 22, 2004, was selected to evaluate the performance of gypsum indices for winter ASTER data. Among the ROC curves, the ROC curve of GI1 was the steepest (closest to the top left in Fig. 10a), and the area under ROC curve of GI1 was the greatest, indicating that the GI1 was the most accurate index for GF mapping using winter ASTER data. The (B4 + B8)/B6 was the second accurate index for GF mapping after GI1 [the ROC curve of (B4 + B8)/B6 was closest to GI1, and the area under the ROC curve is the second greatest after GI1. Area under the ROC curves of the GI1, (B4 + B8)/B6, B4/(B6 + B9), (B10*B12)/(B11*B11) and GI2 were all more than 95% (Table 3), indicating good performance for GF mapping, using winter ASTER data. Using optimum index thresholds, the GI1 achieved the greatest overall accuracy (94.2%), the greatest GF producers' accuracy (93.8%), the lowest combined error (12.1%), second greatest GF users' accuracy (59.2%), and second least FPR (6%). Besides, the overall accuracy of the (B4 + B8)/B6, B4 + B8/B6, (B10*B12)/(B11*B11), GI1 and GI2 were more than 90%, FPR of the (B4 + B8)/B6, B4/(B6 + B9) and GI1 are less than 7%, TPR (identical to producers' accuracy of GF) of the (B4 + B8)/B6, B4 + B8/B6, (B10*B12)/(B11*B11) and GI1 were more than 90%, and the lowest combined errors of (B4 + B8)/B6, (B10*B12)/(B11*B11) and GI1 were less than 20%.

Fig. 10
figure 10

The accuracy of mapping gypsum, anhydrite formations for a range of index thresholds, using Receiver Operator Characteristic (ROC) curves. The circles represent the optimum index thresholds. a ROC curves of GF mapping for winter ASTER data (2004.11.22). b ROC curves of GF mapping for summer ASTER data (2003.7.15)

Full size image

Omission errors were calculated as the number of GF pixels incorrectly classified as a percentage of all GF pixels in each GF category. Chm pixels had fewer omission errors than HS2 pixels across most of the gypsum indices. The GI1 had fewer omission errors for Chm pixels than other indices, while the B8/B6 had fewer omission errors for HS2 pixels. The GI1 and GI2 had 4–6% omission error, the (B10*B12)/(B11*B11), (B4 + B8)/B6 and B4 + B8/B6 had 6–9% omission error, and the (B4 + B8)/B6, B8/B6 and B4/B9 had 10–18% omission error for Chm pixels. While the B8/B6 had an omission error of less than 1%, the B4 + B8/B6 had an omission error of less than 10%, the (B4 + B8)/B6, GI1, and (B10*B12)/(B11*B11) had an omission error of 10–21%, the GI2, B4/(B6 + B9) and B4/B9 had an omission error of more than 50% for HS2 pixels.

Commission errors were calculated as the number of GF pixels incorrectly classified as a percentage of all BGF or NGF pixels in each BGF or NGF category. The B4 + B8/B6 achieved greater commission errors for each BGF or NGF (22–94%). Most of the indices had fewer commission errors of Bk (1–9%), except the B8/B6 (25.21%) and B4 + B8/B6 (71.68%). Due to gypsum layers and interbeds contained in BGF, commission errors were greater for BGF (HS1, Cpm-Mln) across nearly all indices. And, the GI1 and GI2 showed less commission than (B4 + B8)/B6 for HS1, while the (B4 + B8)/B6 had less commission than the GI1 and GI2 for Cpm-Mln. Grm and Aj-Mn commission errors were the fewest for the (B4 + B8)/B6 (5.58% and 2.37%), and the GI1 also had fewer commission errors for Grm (6.00%) and Aj-Mn (2.55%). Commission errors of Qhal were fewer for the GI1 and GI2.

The GI1 and (B4 + B8)/B6 showed better accuracy for GF mapping, using winter ASTER data (2004.11.22). Thus, they were applied to map GF in the study area, using the optimum index thresholds. The results were qualitatively compared, superimposing with the geologic formations boundary. As shown in Fig. 11, both the gypsum indices described the distributions of GF in the study area successfully, and the distribution areas of omission pixels for the two indices were similar, attributing to terrain shadows in ASTER data. We focused on comparing the distribution of commission pixels in the results. In the result of (B4 + B8)/B6, more pixels from alluvail terraces (Qhal) were misclassified into GF, while fewer Qhal commission errors for the GI1, indicating the vital role of band7 (2.26 μm) in excluding the interference of clay minerals for GF mapping. The commission of HS1 was a problem for both GI1 and (B4 + B8)/B6, and less commission for the GI1. Moreover, the HS1 commission pixel distribution of the GI1 matched better with the reported distribution characteristics of gypsum layers and interbeds in HS1. Besides, more pixels of Cpm-Mln were classified into GF for the GI1, but these commission pixels were distributed as east–west stripes, consistent with distribution characteristics of gypsum interlayer in Cpm-Mln (Fig. 11).

Fig. 11
figure 11

The GF mapping performance comparison between the (B4 + B8)/B6 and the GI1 for winter ASTER data (2004.11.22), applying the optimum index thresholds. a Mapping results of the (B4 + B8)/B6. b Mapping results of the GI1

Full size image

Evaluation of gypsum indices for summer ASTER data

A scene of ASTER image, acquired on July 15, 2003, was selected as representative of summer ASTER data to access the accuracy of gypsum indices for GF mapping. Among the indices, the (B4 + B8)/B6 showed the best performance with the steepest ROC curve (Fig. 10b) and the greatest area under ROC curve. The GI1 was second only to the (B4 + B8)/B6 in the accuracy of GF mapping, the ROC curve of GI1 was the closest to (B4 + B8)/B6, and the area under ROC curve of GI1 was second only to (B4 + B8)/B6. The (B4 + B8)/B6, GI1, and GI2 showed better accuracy than other indices with an area under ROC curve of more than 95% (Table 4). While, the B8/B6, B4/B9, B4 + B8/B6 and B4/(B6 + B9) exhibited poor performance in GF mapping with areas under ROC curve of less than 90%. When optimum index thresholds were applied, the B4/B9 showed the greatest overall accuracy (96.8%) and TPR (96.8%) but the least GF users’ accuracy (13.7%), indicating a large number of commission errors. The (B4 + B8)/B6 achieved the second greatest overall accuracy (95.9%), TPR (95.9%), and GF users' accuracy (72.5%), meanwhile, the least combined error (7.7%) and second least FPR (3.5%). Though the GI1 had the greatest GF users’ accuracy (76.3%) and the least FPR (2.7%), but lower overall accuracy (90.6%) and TPR (90.5%). Besides, the overall accuracy and TPR of the GI2 were more than 90%. The B8/B6, B4/(B6 + B9), (B10*B12)/(B11*B11) and GI2 had 30–40% GF users’ accuracy.

Omission and commission errors were also calculated for summer ASTER data. Chm pixels had the least omission error (2.2%) for the GI2. And, the B4/B9, (B4 + B8)/B6, and GI1 had 2–6% omission error, the B4 + B8/B6, B8/B6, and (B10*B12)/(B11*B11) had 9–16% omission error, while B4/(B6 + B9) had omission error of 22.04% for Chm pixels. HS2 pixels had more omission errors across most indices, except B8/B6 (1.16%). The B4/B9, (B4 + B8)/B6 and (B10*B12)/(B11*B11) had 12–37% omission error, the GI1, GI2, B4 + B8/B6, and B4/B9 had 67–88% omission error for HS2.

Commission errors were equally severe with summer ASTER data for BGF (HS1, Cpm-Mln). The (B4 + B8)/B6 showed the best performance in NGF commission errors, with a commission error of less than 9% for each NGF. While the B4/B9 achieved greater commission errors for each BGF or NGF (13–99%). Grm pixels had the least commission (0.41%) for the GI1. And, the B8/B6, (B4 + B8)/B6 had 1–1.6% omission error, the (B10*B12)/(B11*B11), B4/(B6 + B9), and GI2 had 13–39% omission error, while B4 + B8/B6 had omission error of 95.50% for Grm pixels. Qhal pixels had less commission (0.19% and 1.35%) for the GI1 and GI2, proving again that band7 effectively eliminates the interference of clay minerals for gypsum mapping. Most of the indices had fewer commission errors of Aj-Mn (0–8%), except the B4/(B6 + B9) (24.46%) and B4 + B8/B6 (50.76%). In addition, the commission of shadow pixels in As-Ja was a problem for most indices (with 18–91% omission error), except the B4/(B6 + B9) and (B4 + B8)/B6 (with 3–9% omission error). Bk pixels had very few commission errors (0.08%, 0.60% and 1.13%) for the GI1, (B4 + B8)/B6 and GI2, while the (B10*B12)/(B11*B11), B4 + B8/B6, B4/(B6 + B9), and B8/B6 had 5–17% omission error for Bk pixels.

Mapping results of the GI1 and (B4 + B8)/B6, applying to summer ASTER data of July 15, 2003, were also compared qualitatively. Both indices extracted the spatial distribution of Chm member in the study area successfully; however, more omission errors were observed in the marginal zone of Chm regions of the north study area for the GI1. And the distribution areas of other omission pixels in Chm were similar, for the two indices, attributing to terrain shadows or anhydrite distribution. While mapping results of HS2 were quite different (Fig. 12 region A). The (B4 + B8)/B6 showed significantly better performance on recognition of HS2 pixels than GI1, with fewer omission errors. Besides, more commission errors were observed in the results of (B4 + B8)/B6) for Quaternary deposits (Qhal) and HS1, attributing to the overlapping absorption feature on 2.21 μm between gypsum and hydroxyl-containing minerals. While more As-Ja pixels were misclassified into GF for the GI1.

Evaluation of gypsum indices for multi-seasonal ASTER data

Three scenes of summer ASTER images (2001.7.15, 2002.7.28, 2003.7.15) and three scenes of winter ASTER images (2000.11.18, 2004.3.11, 2004.11.22) were selected to evaluate the accuracy of gypsum indices for multi-seasonal ASTER data. ROC curves of (B4 + B8)/B6, GI1 and (B10*B12)/(B11*B11) were plotted (Fig. 13), areas under ROC curves and the lowest combined errors were calculated (Table 5). Overall, the (B4 + B8)/B6, B4/(B6 + B9), (B10*B12)/(B11*B11), GI1 and GI2 showed better accuracy than other indices for multi-seasonal ASTER data in the study area, with an average area under ROC curves of more than 95%. Furthermore, the (B4 + B8)/B6 and GI1 exhibited the most robust mapping accuracy and with areas under ROC curves of more than 98% for six ASTER scenes. For summer ASTER data, the (B4 + B8)/B6 showed the best accuracy with the largest average areas under ROC curves of 99.2%, the lowest combined error of 8.6%. The performance of GI1 and GI2 was secondary to (B4 + B8)/B6, with an average area under ROC curves of 98.5% and 95.6%, while, for winter ASTER data, the GI1 acquired the best accuracy, with the largest average area under ROC curves of 98.7%, the lowest combined error of 10.1%. Besides, the areas under ROC curves of (B4 + B8)/B6, B4/(B6 + B9), (B10*B12)/(B11*B11) and GI2 were all more than 95%, exhibiting good mapping performance for winter ASTER data.

In addition, different seasonal variation characteristics of GF mapping accuracy were observed for different gypsum indices. The (B4 + B8)/B6, B8/B6 and B4 + B8/B6 had greater areas under ROC curves for summer ASTER data. While the B4/(B6 + B9), B4/B9, (B10*B12)/(B11*B11), GI1, and GI2 exhibited greater areas under ROC curves for winter ASTER data. It can be observed in Fig. 13a that ROC curves of (B4 + B8)/B6, applied to summer ASTER data are closer to the left top than winter ASTER data. And the (B4 + B8)/B6 had a greater average area under ROC curves and lower combined error for summer ASTER data (Table 5), indicating better performance for summer ASTER data.

While the GI1 showed better accuracy for winter ASTER data with larger areas under ROC curves and lower combined error. As shown in Fig. 13b, ROC curves of winter ASTER data are steeper. Although ROC curves of 2004.11.22 and 2003.7.15 are close, GI1 achieved better performance on the identification of GF (HS2) pixels in the north of the study area for 2004.11.22 ASTER data (Fig. 11 and Fig. 12). In addition, HS2 pixels had more omission errors for 2003.7.15 ASTER data, and the overall accuracy and TPR of the GI1 was greater for 2004.11.22 ASTER data (Tables 3, 4). For the (B10*B12)/(B11*B11), there was a huge accuracy difference between summer and winter ASTER data. As shown in Fig. 13c, ROC curves for winter ASTER data are much closer to the left top than summer ASTER data. The average area under ROC curve for winter ASTER data was 97.5%, while 93.5% for summer ASTER data. The average lowest combined error for winter ASTER data was 15.4%, while 28.9% for summer ASTER data. Moreover, TPR for winter ASTER data were all more than 90%, while less than 90% for summer ASTER data. Besides, the average area under ROC curve of B4/B9 for winter ASTER data (94.7%) was much greater than summer ASTER data (80.4%), indicating the considerable accuracy difference between summer and winter ASTER data.

Fig. 12
figure 12

The GF mapping performance comparison between the (B4 + B8)/B6 and the GI1 for summer ASTER data (2003.7.15), applying the optimum index thresholds. a Mapping results of the (B4 + B8)/B6. b Mapping results of the GI1

Full size image
Table 3 Gypsum and anhydrite formations classification accuracy statistics for the different index methods, using winter ASTER data (2004.11.22)
Full size table
Table 4 Gypsum, anhydrite formation classification accuracy statistics for the different index methods, using summer ASTER data (2003.7.15)
Full size table
Fig. 13
figure 13

The accuracy of mapping GF for a range of index thresholds, using Receiver Operator Characteristic (ROC) curves. The circles represent the optimum index thresholds. a ROC curves of the (B4 + B8)/B6 for multi-seasonal ASTER data. b ROC curves of the GI1 for multi-seasonal ASTER data. c ROC curves of the (B10*B12)/(B11*B11) for multi-seasonal ASTER data

Full size image

Discussion

Origin of gypsum formation reflectance seasonal variability

ASTER band7 (2.26 µm) appeared obvious seasonal variability for Chm spectral member and GF training samples in this study. It has been proven that the 2.26 µm absorption features of Ca-sulfates were attributed to the proportions of structure molecular water in Ca-sulfates (Hunt 1971; Crowley 1990, 1991; Drake 1995; Cloutis et al. 2006; Harrison 2012; Bishop et al. 2014). Besides, experiments of Bishop (2014) and Harrison (2012) showed that the 2.21 µm and 2.26 µm absorptions decrease in depth, and the absorption positions shift to shorter wavelengths (Harrison 2012), as the gypsum samples dehydrate to bassanite or anhydrite. The climate of the study area is similar to Kuwait, and Ca-sulfate mineral facies were reported to change seasonally in there (Gunatilaka et al. 1985). It can be inferred from these results that seasonal reflectance variability of GF pixels at 2.26 µm may be a response to alternating transformation between bassanite and gypsum with the season. High temperature, extreme drought, and very little vegetation cover provided possible conditions for transformation from gypsum to bassanite in summer (Rezaee and Salari1 2016; Warren 2006). Bassanite was generally considered to be a metastable Ca-sulfates during the conversion between gypsum and anhydrite (Sievert et al. 2005), and only existed steadily in extremely arid climate conditions (Bishop et al. 2014). The summer climate of study area is such arid conditions with few rainfalls, providing probably suitable environment for bassanite. Meanwhile, bassanite was also thought to be a more stable Ca-sulfate phase rather than continuing to rehydrate to form gypsum, without being directly exposed to liquid water (Harrison 2012). Thus, the concentrated precipitation in winter may be the key factor in converting bassanite to gypsum in the study area.

Due to the limitation of the research conditions, we do not verify the mutual transformation of gypsum and bassanite between summer and winter in the field survey. And, it should be noticed that GF reflectance seasonal variability is also under the influence of the different solar elevations and radiation intensity between summer and winter ASTER data. It has been proven that Landsat digital numbers for rocks are proportional to the cosine of the solar zenith angle (Kowalik et al. 1982), and reflectance values are proportional to the slope of the digital numbers versus cosine of the solar zenith angle (Kowalik et al. 1982).

The performance of the proposed new gypsum indices

In this study, two ASTER gypsum indices (the GI1 and GI2) were developed. Accuracy assessment showed that distributions of GF are identified successfully by the GI1 for both summer ASTER data and winter ASTER data. And for winter ASTER data, the GI1 improved the GF mapping accuracy of the study area, achieving the steepest ROC curve (Fig. 10a), the largest average area under ROC curve (98.7%, Table 5), and the lowest average combined error (10.1%, Table 5). While, for summer ASTER data, GF mapping accuracy of the GI1 is only second to (B4 + B8)/B6 with an average area under ROC curve of 98.5% and average lowest combined error of 14.3%. Besides, the (B4 + B8)/B6 also showed better performance for gypsum mapping in the study area, consistent with the results of Özyavaş (2016).

Table 5 Gypsum and anhydrite formation classification accuracy statistics for the different indices, using multi-seasonal ASTER data
Full size table

Similar to the results of previous studies, ASTER band6-based gypsum indices are susceptible to the interference of hydroxyl-bearing minerals. Not only pure gypsum rocks but also other authigenic clays and disseminated gypsum crystals were mapped by (B4 + B8)/B6 in the Kohat Plateau, north Pakistan (Khan et al. 2020). “False kaolinite” identification was reported for SWIR hyper-spectral data in Israel, attributing to the proximity of the absorption at 2.20 µm of gypsum and kaolinite (Notesco et al. 2016). In this study, commission errors were more severe for Qhal and HS1, geologic formations with more hydroxyl-bearing minerals, across the (B4 + B8)/B6, B8/B6, and B4 + B8/B6 (Tables 3 and 4). While, the GI1 and GI2 showed less commission than the three gypsum indices for Qhal and HS1 (Tables 3 and 4). There are usually large amounts of clay minerals in alluvail terraces and recent deposite (Qhal). As shown in Fig. 11, commission pixels of Qhal, distributed on river channels and alluvial fans, for the GI1 are much fewer than the (B4 + B8)/B6. Besides, fewer HS1 pixels were misclassified into GF for the GI1. Previous studies (Bosák et al. 1998) showed that gypsum-containing layers or interlayers are distributed in HS1 in strips, and there is also a large amount of chlorite, epidote, illite, and montmorillonite minerals, formed as secondary weathering minerals on the surface of HS1 regions. These hydroxyl-bearing minerals have an overlapping 2.21 μm absorption position with gypsum. The mapping results of (B4 + B8)/B6 and GI1 in the HS1 region are compared in Fig. 11, it was found that commission pixels of HS1 distributed in several strips for the GI1, consistent with the distribution of gypsum layers and interlayers in HS1. While, in the results of (B4 + B8)/B6, a large number of commission pixels scattered in HS1 disorderly, attributing to distributions of both the gypsum layers, interlayers, and secondary weathering minerals on the surface of rocks. These indicated that the proposed GI1 successfully identifies gypsum layers or interlayers, and effectively eliminates the interference of the hydroxyl-containing minerals in the background.

We also found that the GI2, established using the 2.26 µm absorption feature of gypsum alone, did not achieve the highest accuracy (the mapping accuracy of GI2 is inferior to (B4 + B8)/B6 and GI1 for multi-seasonal ASTER data with an average area under ROC curve of 96.2%). While the GI1, combined using the 2.21 µm and 2.26 µm absorption features of gypsum, showed the best mapping accuracy for winter ASTER data, and the second highest mapping accuracy for summer ASTER data. Furthermore, the GI1 is the most robust gypsum index for multi-seasonal ASTER data with average areas under ROC curve of 98.5% and 98.7% for summer and winter ASTER data. In summary, the 2.21 µm strong absorption feature is critical for gypsum mapping, while the 2.26 µm absorption helps improve the interference of the hydroxyl-containing minerals for gypsum mapping. However, we also found that commission error of carbonates is a problem for the GI1 and GI2, a few pixels of carbonate rocks (As-Ja) were misclassified into GF. Thus, background lithology also should be considered for the selection of gypsum indices.

The performance of gypsum indices for multi-seasonal ASTER data

Previous studies reported that MgOH/carbonate RBD products (B6 + B9)/B8 show better accuracy for summer ASTER data than winter ASTER data, attributing to the difference of terrain shadow and sun altitude (Hewson and Cudahy 2011). Nevertheless, the accuracy of different gypsum indices showed markedly different patterns in terms of the effect of ASTER acquisition season in this study. The B4/(B6 + B9), B4/B9, (B10*B12)/(B11*B11), GI1, and GI2 showed better accuracy for winter ASTER data, and the GI1 achieved the best accuracy. While the (B4 + B8)/B6, B8/B6, B4 + B8/B6 showed better accuracy for summer ASTER data, and the (B4 + B8)/B6 achieved the best accuracy. We also found that gypsum indices, calculated mainly by dividing by the reflectance of band7 (2.26 µm) or band9 (2.40 µm), exhibited better accuracy for winter ASTER data. While gypsum indices, designed based on the band6 (2.21 µm) gypsum absorption, show better accuracy for summer ASTER data. Experimental results showed that seasonal variability of accuracy for gypsum indices is caused by the seasonal change of reflectance curve shape in SWIR region. The weakening of ASTER band7 (2.26 µm) reflectance absorption feature of Chm pixels for summer ASTER data is the reason for the better accuracy of the GI1 and GI2, applied to winter ASTER data.

Conclusions

Three major conclusions are drawn from this study.

  • (1) Seasonal reflectance variability of gypsum pixels was confirmed, and the seasonal variability of 2.26 µm reflectance may be a response to mutual transformation between bassanite and gypsum with the season.

  • (2) The effectiveness of 2.26 µm (ASTER band7) reflectance in gypsum mapping was validated. On the basis of 2.21 µm and 2.26 µm absorption features of gypsum, two new ASTER gypsum indices (GI1 and GI2) were proposed, improving the interference of hydroxyl-containing minerals for gypsum mapping. Furthermore, the GI1 achieved the most robust mapping accuracy for multi-seasonal ASTER data.

  • (3) Gypsum indices showed markedly different accuracy for multi-seasonal ASTER, attributing to the seasonal reflectance variability of GF pixels. The B4/(B6+B9), B4/B9, (B10*B12)/(B11*B11), GI1, and GI2 exhibited better accuracy for winter ASTER data. While, the (B4+B8)/B6, B8/B6, B4+B8/B6 show better accuracy for summer ASTER data.

The next step will be to verify the reason for the seasonal reflectance variability of gypsum pixels in the arid region and explore the possibility of subdividing Ca-sulfates types (anhydrite, bassanite, and gypsum) by ASTER data. The proposed GI1 can be applied in similar arid regions for gypsum mapping.