A comprehensive scheme for lithological mapping using Sentinel-2A and ASTER GDEM in weathered and vegetated coastal zone, Southern China

3.1 Remote Sensing Dataset and Image Pre-processing

In this contribution, both the Copernicus Sentinel-2A (S2-A) level-1C image on 11 February 2017 and ASTER GDEM v2 data (30-m resolution) were selected for use in Table 1. As a new mission, the sentinel satellites can provide routine observations and data continuity of ERS, SPOT, ASTER, etc. [25]. The S2-A was launched on 23 June 2015 for the operational needs of the Copernicus programme, carrying a single Multi-Spectral Instrument (MSI) in 13 spectral bands from the Visible and the Near Infrared (VNIR) to the Shortwave Infrared (SWIR) at three spatial resolutions (10, 20, 60 m) [26]. S2-A and another Sentinel-2B launched on 7 March 2017 comprise a system of twin-polar orbiting satellites, phased at 180^∘ to each other, and the cooperation of the twin satellites provides a temporal resolution of 5 days. Results provided by Reference [27] also show that Sentinel-2 products are of high quality and consistency with the good radiometric calibration.

Here, we used S2-A Level-1C image because it currently provides freely available multispectral bands at relatively high spatial resolution comparing to other free data, which can satisfy the fundamental demands of outcrops extraction. Besides, previous study [28] has demonstrated that S2-A data turn out to be adequate for lithological discrimination, which yielded higher classification accuracy than Landsat 8 OLI and ASTER imagery. Therefore, two adjacent scenes were freely downloaded from the Sentinels Scientific Data Hub (https://scihub.copernicus.eu/) in the JPEG 2000 format based on several criteria including acquisition date, cloud (~2%) and vegetation cover. The images were orthorectified in the UTM/WGS84 projection and geometrically corrected.

Comparing to ASTER GDEM v1, ASTER GDEM v2 data are characterized by increased horizontal and vertical accuracy, larger coverage and less artifacts. The image was obtained from Geospatial Data Cloud site, Computer Network Information Center, Chinese Academy of Sciences (http://www.gscloud.cn). Apart from the spectral information of S2-A, topographic features derived from DEM contribute to exclude pseudo extracted outcrops further. In addition, as depicted in Table 1, the complementary references such as the high-resolution Google Earth^TM images, camera photos (~2000) and coordinates (~3 accuracy) of field observation points (tuff: 783; granite: 59) were also utilized for accuracy assessment.

In this study, those data-sets were (pre-) processed using the Sentinel Application Platform (SNAP) v6.0 and ENVI v5.3 software packages. We conducted the atmospheric correction to S2-A images using sen2cor plugin v2.5.5 [29] on SNAP, transforming the reflectance from Top-Of-Atmosphere (TOA) Level 1C to Bottom-Of-Atmosphere (BOA) Level 2A. This plugin is developed based on the Atmospheric and Topographic Correction for Satellite Imagery (ATCOR) [30] and it can reduce atmospheric and topographic effects in satellite imagery and extract physical surface properties based on the libRadtran radiative transfer model [31]. After atmospheric correction, only the 10m spatial resolution VNIR bands and the 20m spatial resolution Vegetation red edge and SWIR bands of Sentinel-2 data were considered. All these ten image bands were cubically resampled to 10 m spatial resolution and layer stacked to one file. Then, the two scenes were mosaicked based on the seamless mosaic workflow built in ENVI 5.3; for another, ASTER GDEM v2 image was also cubically resampled from 30m to 10m spatial resolution in order to match with the S2-A image. At last, the sentinel-2 data and ASTER GDEM v2 image were clipped over to the spatial extent of the study area.

3.2 Spectral Features of Rocky Outcrops

As mentioned, tuff strata and granite plutons dominate the area. In order to reveal the spectral features of different outcrops, several representative image-derived observation points were selected (T1-T4 and G1-G4 in Figure 1). Regarding to Figure 3a, reflectance curves of tuff and granite outcrops appear to be similar which exhibit a profound reflection peak at Band 11 (SWIR, central wavelength: 1610nm), while higher reflection values of granite outcrops might be attribute to the enrichment of quartz mineral in high reflectivity during weathering [32]. All the reflection valleys appear at Band 2 (Blue, central wavelength: 490nm).

Figure 3

(a) Spectral Z-profile of representative bedrock outcrops (T1-T4 and G1-G4 in Figure 1) derived from the S2-A image mentioned later, using the Z-profile (spectrum) module built in ENVI 5.3. (b) Spectral Z-profile of background features in the S2-A image.

Figure 3b illustrates the spectral characteristics of different objects derived from the study area. It further reveals that the vegetation and water bodies exhibit distinct spectral responses from rocky outcrops. However, there are only subtle differences in the reflectivity level between two lithologies and soil, building as well as road. These similarities will hinder spectrum-based outcrops recognition and introduce a lot of false detected results. Herein, tuff outcrops show the more similar spectral characteristics to the background features than granite outcrops, thus probably resulting in more pseudo detected outcrop patches.

3.3 Matched Filtering

As a powerful and widely used tool in the communication and signal processing domain, MF has been further developed to the application of remote-sensing target detection [33] recently. It is a partial unmixing approach with a normal distribution of the target and background classes, which can maximize the response of the endmembers derived from the image, field measurement or laboratory, and minimize the response of the composite unknown background [34]. Finally, a single-band image with matching degree of each pixel value ranging from 0 to 1 (where 1 represents a perfect match), can be obtained

In this project, we used the image-derived reference spectra, that is, eight typical points (or ROIs - regions of interest) in the image corresponding to the known rocky exposures (T1-T4 and G1-G4 in Figure 1), as the matching templates for the MF-based preliminary extraction of bedrock outcrops. Those carefully selected outcrops are all well-exposed and have a planar expression on the image, or in other words, can be seen remotely. Chemical weathering process is the key reason for the selection of the image-derived reference spectra instead of laboratory measured spectra. Another reason to choose as many as eight outcrop points is that one type of rock may have several different spectra. These preparations are prerequisite for forming a more effective template.

For the purpose, running the MF technique, using reference spectra of the representative points, can enhance the outcrop-related information with bright pixels. However, as expected in the spectral analysis section, soil, building and road pixels in the image were also highlighted in Figure 4a and 4b. Herein, the MF-based image of granite obviously outperforms the results of tuff, i.e., the background in granite enhanced image (Figure 4b) were well suppressed. And despite using the fractal method mentioned later (Figure 7a and 7b), pixel patches related to the rocky outcrops as well as soil, building and road were obtained after MF. In Figure 4c, abundant building, linear rural highway and patches in granite plutons were extracted due to similar spectral features. In Figure 4d, upon the visual interpretation, most resulting patches fell into the granite plutons, some patches were identified into tuff outcrops and some patches also fell into the background features. Hence, in order to obtain accuracy lithological mapping in the sparsely outcropping areas, more processings independent of spectral nature must be considered further.

Figure 4

(a) and (b) Matched filter result of tuff and granite outcrops using the image-derived spectra. (c) and (d) Resulting patches of tuff and granite following MF analysis based on the fractal algorithm mentioned later.

3.4 Spectral Indices and Terrain Analysis

Spectral indices are a wide and common approach with orthogonal or mathematical transformations of multispectral data, which can suppress unnecessary effects and produce more sensitive spectral reflectance [35]. Additionally, NDXI (X: vegetation, soil and water) also helps compensate for differences both in illumination within an image due to slope and aspect in the rugged terrains. Here, the spectral indices used contain NDVI, NDSI and MNDWI representing the existence of vegetation, soil and water, respectively. The equation of the three indices are calculated as:

Where ρGreen, ρRed, ρNIR and ρMIR are the reflectance of S2-A data bands 3, 4, 8 and 11, respectively. Note: NDVI is the world’s well-known normalized difference vegetation index; MNDWI is the modified NDWI proposed by Xu [36] to reveal subtler features of water body; and NDSI is the normalized difference soil index [37].

The false color composite image of NDXI is shown in Figure 5a. In theory, NDXI values range from −1.0 to +1.0. Higher values of NDXI they are, more likely the area is to be covered by the corresponding enhanced objects. In this regard, geological exposures with removal of vegetation and soil and non-water coverage, resulted in decreased values of NDXI. Based on this principle, the objective of these indices was to mask or filter out the pseudo detected “outcrops” caused by the complex mixture of ground objects during MF processing. It is worth noting that there is no order of priority between these two steps, and if performing the spectral indices firstly, the entire study area will be restricted to a limited interpreted region mainly focusing on exposed rocks for the subsequent outcrop detection.

In this case, the relevant cutting thresholds for the difference spectral indices were over 0.322 for NDVI, over 0.302 for NDSI, and over 0.202 for MNDWI, and they were determined according to both empirical values and the internal characteristics of the study area. Next, the extracted three vector shapefiles were combined via spatial union analysis into a mask (Figure 5b) to eliminate the pseudo detected “outcrops”.

Figure 5

(a) False color composite image made with the spectral index products. (b) Mask consisting of NDXI vector shapefiles including NDVI (green), MNDWI (blue) and NDSI (brown).

Except the spectral information, terrain features are proved to be another important indicator of geology [38, 39] such as the identification of outcrop occurrence [14, 16]. On the one hand, most common geological exposures in the field tend to occur in the rugged or steep terrain. In this region, the occurrence of outcrops is also the case. Obviously, the spectral reflectance patterns of bedrock will help us to identify the outcrops with similar spectral characteristics and exclude the different ones. On the other hand, crushed stones could be removed away from bedrock exposures to other flat areas by transporting action of running water, weathering, human excavations, etc. Inasmuch, the terrain or micro-geomorphology also has the potential to distinguish in situ from ex situ outcrops pertinent to spectral similarities.

Furthermore, there are subtle difference in topographic features between the tuff strata and granite plutons, that is, granite plutons are characterized in relatively gentle landforms comparing with tuff strata, because the latter is prone to suffer from weathering, particularly spheroidal weathering [40]. Notably, here only the overall topographic features of tuff and granite were considered, and inevitably, some granite outcrops can occur in steep region and some tuff outcrops are found in flat area. In order to remove most of false detected patches, this analysis will come at the price of the loss of true patches.

Summarily, slope attribute was selected as the interpreted sign to limit the outcropping areas of tuff strata and granite plutons. Slope product extracted from ASTER GDEM (30m) was incorporated for improving the outcrops extraction results. In Figure 6, the slope layer (T₁₋₂-T₁₋₃, granite plutons) in red color with 50% transparency where granite outcrops are most likely to occur, was overlaid on the hillshade map consistently. Detected granite outcrop patches that fall into the red layer were further remained. Then, according to the different topographic features, threshold of (T₁-max) was utilized to limit the potential outcropping areas of tuff outcrops. The following results or accuracy also prove the feasibility of terrain analysis. Note that: the thresholds were determined by fractal algorithm mentioned later (Figure 7c).

Figure 6

The matching slope layers (red: granite; blue: tuff) overlying on the hillshade map, derived from ASTER GDEM.

Figure 7

Fractal ln (N) versus ln (DN) schema of MF: (a) tuff; (b) granite, and (c) Terrain analysis. Note that: the wireframe represents the range of selected threshold.

3.5 DN-A Fractal Scheme

Setting appropriate thresholds is essential to yield accurate detection results or reduce the interpretation target area in the enhanced image. Here, the fractal C-A (concentration-area) algorithm, which is one of the most advanced unsupervised classification methods [41], was introduced to determine the cutting thresholds of MF and terrain analysis, respectively. With respect to the fractal method, Shen [42] suggested that it can be used to constitute the power-law relationship magnitude and frequency on a log–log paper, and a number of straight-line segments (scale-less intervals), which usually correspond to the real-world features (classes) on the ground, can also be manually or automatically obtained to yield a group of cut-off values for image segmentation or target detection. Inasmuch, for pixel values and pixel value frequency distribution as well as spatial and geometrical properties of image, the DN-A (equivalent to C-A) method [43] can be expressed as follows:

where A (DN) denotes the area occupied by pixels with DN values equal to and greater than a certain threshold (s), and the exponent D is fractal dimension.

Figure 7 shows the fractal schema for the MF and terrain analysis in this case. Note that: 1) The scaleless line segments reflecting the features of interest (real pixel patches) are determined by visual interpretation and marked by black wireframes; 2) R² indicates the goodness of fit, and R² of each linear segment is kept ≥ 0.9 (except the segment D₂ of Figure 7b); 3) N₁, N₂... are the proportion that pixels falling within the segment D₁, D₂..., respectively, take up of the total pixels involved in fractal calculation; 4) T_n−1 (n=2, 3, 4...) is the threshold value between two adjacent segments; 5) The parameters for barren D₁ and D₂ are also displayed in this diagram; 6) The “oscillation” of D₃ segments might be attributed to the mixture pixels or random noise (Figure 7a); 7) The computer programme in MATLAB language is available from Zhao et al. [44] on request.

According to these criteria, T₁-max was selected as the MF thresholds of tuff and granite. With respect to terrain analysis, thresholds of T₁₋₂-T₁₋₃ and T₁-max were used to limit the outcropping areas of granite and tuff, respectively.

3.6 Accuracy Assessment

In generally, the Kappa coefficient is often adopted as a standard measure of classification accuracy [45], which measures the agreement between a classification image and its ground truth. However, actually there is no standard or reference basemap to validate the accuracy of the outcrops detection results, so the Kappa coefficient was inapplicable for accuracy assessment in this study. In addition, it should be admitted that in these cases, it is impossible to validate the whole detected results using a single criterion, just like the way using multiple methods to extract geological outcrops. Thus, the semi-automatic multilevel validation method is expressed as follows:

• A detected pixel patch falling within the known lithostratigraphic region (100-m buffered) may be real (see Figure 4, 8 and 9). The procedure was to overlay the vector shapefiles of the lithostratigraphic regions and the acquired pixel patches;

Figure 8

(a) Obtained tuff patches on the S2-A true color image; the number annotations (1~4) and transparent wireframes in the image are consistent with the subplots number of 1~4. Frequency (b) and cumulative frequency (c) distribution histogram of the patch areas. I-III are the verification photos and Google earth image of tuff.

Figure 9

(a) Obtained granite patches on the S2-A true color image. Frequency (b) and cumulative frequency (c) distribution histogram of the patch areas from the minimum to an area of 20 km² with a bin size of 0.5 km². I-III are the verification photos and Google earth image of granite.

• A detected pixel patch that is near to a field validation point (< 0.1km) may be real (Figure 8, 9), especially when it appears in the same gully/valley. The procedure was to create a 100m buffer for each of the validation points (842 point-vectors) firstly, and overlay the vector shapefiles of the validation points and the obtained pixel patches;

• A detected patch corresponding to a non-vegetated “outcrop” in the Google Earth™image may correspond to real-world outcrops;

• Other visual cues. For example, rocky outcrops exposed by water erosion or human engineering.

Here, we assumed that any detected pixel patch satisfying case (2), or simultaneously satisfying case (1), (3) and (4) can be classified as “real outcrop”, and the interpretation ratio is simply “% of the correctly classified patches over the total number of obtained patches”.

4.1 Tuff Outcrops

After the previous series of operations, the resulting patches of tuff outcrops are displayed in Figure 8, therein lies 3 insert maps zoomed in detailedly. A total of 916 patches (3139 pixels) were extracted to represent the real-world rocky exposures in this image. At first glance, most pixel patches fell within the known stratigraphic terrane (with a 100m buffer). To be specific, the tuff patches are generally small and tend to occur in cluster comparing to granite outcrops (Figure 8-1~3). Those pseudo patches in purple color were recognized falling onto the granite plutons and buildings. Those uncertain patches in green color mainly occurred in the bare land, which cannot be well constrained by the four validation criteria. Figure 8b and 8c illustrate the frequency distribution and cumulative frequency of the tuff patch areas from the minimum to an area of 4 km² with a bin size of 0.1 km², respectively. The peak value appears in the range 0–0.1 km² containing ~51% of patches, which is equal to the area of a pixel in this S2-A image. 69% of patches are < 0.2 km² and approximately 90% of patches are < 0.6 km².

Subsequently, based on verification criteria, Table 2 shows the two correct interpretation ratios (interpretation accuracy) of preliminary extraction (MF), elimination of false patches (NDXI + Terrain) of tuff and granite, respectively. Firstly, the accuracy after MF was 41~43.1% (4589~4827/11200) and it is the foundation of following NDXI + terrain analysis. The truth is that despite without diagnostic spectral features, it is also characterized in topographical and spatial aspect on the DEM image. Secondly, the accuracy after NDXI and terrain reached 82.4~83.2% (755~762/916), which had improved by ~40%, and correspondingly, 91.8% of the raw pixel patches in Figure 8 were removed in this process. The improvement of accuracy also demonstrates that elimination of the false patches can be sometimes more important than extraction of the real ones [16].

4.2 Granite Outcrops

In Figure 9, a total of 575 patches (12122 pixels) were extracted from the S2-A image. Comparatively, most tuff patches were consisted of one or two pixels (~70%), however in addition to showing the same small-scale feature with tuff patches, the granite patches also contained few large pixel patches with over 160 or even 300 pixels (16-30 km²). As presented in Figure 9b and c, the peak value appears in the range 0–0.5 km² containing 60% of patches. 74.1% of patches are < 1 km² and approximately 90% of patches are < 3.5 km². Those large-scale patches, for example 12 patches over 16 km², mainly resulted from human excavations to ore occurring in granite bodies, such as quarrying and open-pits.

As presented in Table 2, the accuracy after MF was 62.6~67.9% (431~467/688), and after NDXI and terrain analysis, the accuracy reached 91.3~92.3% (525~531/575) which had improved by ~15%. Meanwhile, 16.42% of the raw pixel patches in Figure 9 were removed in this process. Notably, the interpretation accuracy of granite is higher than the tuff from MF to NDXI + terrain scenario. Among the reasons for the higher interpretation accuracy of granite outcrops was that they show more diagnostic spectral features comparing to tuff outcrops in Figure 2.

5.1 Spatial and Spectral Resolution

The goal of this project was to extract the pixel patches of geological outcrops from the S2-A data. In the heavily vegetated area, the originally bedrocks are generally covered by vegetation or loess, leaving sparse and limited outcropping areas. Meanwhile, outcrops of various lithologies may have completely different spectral patterns due to different rock compositions (note that the spectral patterns of tuff and granite appear to be similar). Consequent physical or chemical weathering will further alter the structure and chemical component of outcrops, which exacerbates the complexity of lithological mapping. On the contrary, Human excavations to ore such as quarrying and open-pits are among the reasons for making it have greater chances of being detected with remote sensing. But it can be seen from the results that artificial constructions could not always produce large-scale outcrops, i.e. in both sides of the highway, and the corresponding numbers are notably beyond our expectation. Hence, although natural outcrops are a scarce resource in the sparsely outcropping areas in which exactly hinder the development of geological mapping, ecological investigation and so on, abundant human-made outcrops may contribute to reveal more geological information as well as the interactions of humans and nature on the Earth’s surface.

Inasmuch, free Sentinel-2A image with relatively high resolution was used to balance against the accuracy required in this research. However, this data’s spatial or spectral resolution still limit its application in such areas resulting in the cases of few false pixel patches and many undetected outcrops. Besides, despite with a high accuracy of over 90% using multiple methods, there is still room for improvement. Upon visual inspection of the entire pixel patches, the remaining false patches were usually caused by footpath (Figure 8a) and cement-related grave with the similar spectral pattern and topographic characteristics to real outcrops. There are also many outcrops remaining to be undetected. So the availability of finer spatial or spectral resolution sensors such as Worldview-3 [46, 47] may offer a possible partial resolution to the extraction of geological outcrops in such a challenging terrain.

5.2 More Interpreted Signs

Spectral patterns and topographic features serve as the two major RS-oriented interpretation signs in the text, playing a vital role in limiting and locating the outcropping areas. Here, some non-outcrop pixel patches were identified because of similar spectral and topographic characteristics to those real outcrops. So, more interpretation signs are demanded to perform more accurate bedrock outcrops interpretation. As a future development of this work, the proposed approach will be integrated with more interpreted signs including vegetation, thermal inertial property and so on:

Vegetation has generally been a noise or a hindrance of masking geological information, and most of the world’s mineral production are obtained from low to moderately vegetated land surface instead of heavily vegetated areas. Whereas in geobotanical remote sensing [48], the spectral behavior of vegetation that can be responsive to differences in soil-lithology characteristics, is thought to have potential in deciphering geological information. For instance, Wang et al. [16] suggested that the chestnut trees in Hebei province of China can be regarded as an empirical discriminant sign indicating the presence of gneiss strata; Grebby et al. [49] have explicitly tested the impact of vegetation on remote lithological mapping capability and highlighted that useful lithological information can be obtained even with the limited spectral resolution and ubiquitous vegetation cover.

Thermal inertial can act as a universal interpreted sign for distinguishing rocks from other ground objects, because it is an inherent physical property of the rocks [6]. Specifically, rocks (higher thermal inertia), particularly denser and non-porous, tend to retain the heat obtained during the sunlight exposure more effectively than unconsolidated regolith such as soils (lower thermal inertia). In this regard, thermal inertial property may offer more important reference values for outcrops extraction than vegetation.

In summary, the formation of RS-oriented geological interpreted signs is the process of transforming subjective experiences from geologists to executable technical operations in computer. In order to perform accurate RS interpretation, more interpreted signs should be introduced, despite coming at the price of the decreasing number of detected pixel patches, but the proportion of true outcrops in the whole detected patches will increase gradually.

5.3 Accuracy Assessment Criteria

Since the complexity of bedrock outcrops mentioned above, it is not possible to use either single interpreted sign or single method to provide a satisfactory interpretation result. Similarly, a single criterion is also insufficient to validate the whole detected pixel patches in the study area. For this purpose, the authors mainly started from the three aspects containing lithostratigraphic regions, field observation points and finer spatial resolution images, utilizing three corresponding data including geological maps, field surveys and finer resolution RS images to assist verification, respectively. As a direct verification standard, field observation points obviously outperform the other two indirect standards because of their higher survey precision, while the latter are also necessary in areas without access to field validation, but should be used with more consideration.

5.4 Pros and Cons of Weathering Process and the Identification of Protoliths

In coastal zone, chemical weathering is a common process for all rocks, which can change the rock component and produce regolith on the surface. The spectral features of weathered outcrops changed and became unpredictable without field measurement. However, the same kind of lithology produces similar weathering products, while different lithologies shows different weathering characteristics, i.e., differential weathering. The process can form different topographic features of tuff and granite for terrain analysis. Therefore, differential weathering is regarded as a necessary condition for coastal lithological mapping. In addition, when the weathered outcrops were extracted, their protoliths must be defined by the corresponding fresh exposures. In this study, several prospecting trenches were excavated to reveal the protoliths of weathered outcrops. Trenches distributed in the gentle landforms show that the extracted patches characterize a kind of weathering products from fresh granites (Figure 9-II), and those distributed in the steep landforms demonstrate the extracted patches describe weathering products derived from fresh tuffs (Figure 8-I).