Lithology classification of igneous rocks using C-band and L-band dual-polarization SAR data

4.1 Data preprocessing

The extraction of backscatter coefficients, decomposition parameters, and texture information for both Sentinel-1 SLC data and PALSAR SLC data were performed in the Sentinel’s Application Platform (SNAP) toolbox from ESA [36]. All features are shown in Table 1.

The digital number (DN) values of the SAR images were calculated to linear-scale values and then we converted them into decibels (dB). The images were multi-look processing in range and azimuth direction to generate ground ranged square pixels. Filtering was then performed using a refined Lee filter with a 7 × 7 window, as it better preserves the polarization features before filtering with reduced noise [37]. Finally, Doppler terrain correction was performed. SRTM 1 Arc-Second Global elevation data was used as Digital Elevation Model (DEM) for terrain correction, and bilinear interpolation was used as an image interpolation method to resample the resolution of Sentinel-1 data to 10 m and PALSAR data to 15 m. For Sentinel-1, we obtained the backscatter coefficients of VV and VH. Then the ratio VH–VV was calculated. For PALSAR, the ratio HV–HH was similarly calculated, except for the backscatter coefficients of HH and HV.

The decomposition feature extraction requires a complex matrix in radiometric calibration. We used the polarization matrix generation tool to obtain the covariance matrix (C2) of the dual-pol data. Entropy (H), scattering angle (α), and anisotropy (A) extracted based on the covariance matrix were used to describe the scattering mechanism of the rocks. The entropy–angle (H–α) decomposition technique was proposed by Cloude and Pottier [38] for full-polarization data. The parameter α is directly related to the mean physical scattering mechanism. The parameter H defines the randomness of the scattering. H is plotted on the X-axis and α is plotted on the Y-axis. When H is low (H < 0.3), the system is weakly depolarized, and when H is high, the ensemble-averaged scatterer exhibits a depolarized state. The parameter A is used to complement the ratio relationship between the eigenvalues of the covariance matrix and can improve the discrimination of different scatterings when H is large. We also performed multi-look process, filtering, and Doppler terrain correction for decomposition features.

Texture information quantitatively describes relationships of DN values of neighboring pixels [21]. It is a combination of macroscopic and microscopic structures of different topographic features. GLCM was used to extract texture information from SAR data [22]. Haralick et al. [39] proposed various features based on GLCM, and these features are widely used for texture extraction from geological units [24]. GLCM was then applied using a window of 9 × 9 pixels. We obtained the mean, variance, homogeneity, contrast, dissimilarity, entropy, and correlation of GLCM for each polarized image. They represent the contrast features, orderliness features, and statistics features of the VV and VH (HH and HV) images. The quantization levels parameter was set to 64.

4.2 Lithology classification

We used the features extracted from two SAR data to classify igneous rocks, respectively, and explore the applicability of the AdaBoost algorithm in igneous rocks classification using SAR data. The AdaBoost algorithm is based on mixing several classifiers with weak classification ability to become a strong classifier with strong classification ability [40]. The package adabag developed by Alfaro et al. was implemented through R (V.4.0.5) package to perform AdaBoost classification [41]. The number of iterations in the model was set to 100.

To compare the classification performance of AdaBoost, random forest (RF) and XGBoost algorithm were also used to classify igneous rocks. They are based on ensemble methods [42,43]. RF algorithm predicts classes based on inputs by constructing multiple decision trees and is not prone to overfitting [44,45]. RF algorithm used the package randomForest through R (V.4.0.5), developed by Liaw and Wiener [46]. Two important parameters of the RF algorithm are the number of trees (ntree), which indicates the number of decision trees contained in the model, and the number of input variables (mtry), which indicates the number of variables used for the binomial tree in a given node. Ntree was set as 500, and several articles have shown that the model will have stable accuracy with this value while a larger ntree will influence the run rate of the model [47,48]. Mtry is often the square root or one-third of the number of features. In this study, the number of features is 23 and the value of mtry was set to 5. Unlike the bagging ensemble learning of RF, the XGBoost algorithm uses boosting ensemble learning with multiple decision trees with correlation for decision making [49]. The package xgboost developed by Chen et al. was implemented through R (V.4.0.5) package to perform the classification of igneous rocks. We chose to use the default parameters in the xgboost package for classification [50].

Meanwhile, mature machine learning algorithms naïve Bayes, K-Nearest Neighbors (KNN), and Support Vector Machine (SVM) were equally used to explore the feasibility of using dual-pol SAR data for igneous rock classification [51]. Among them, naïve Bayes and SVM were implemented by the e1071 package through R (V.4.0.5) [52], and KNN was implemented by the kknn package through R (V.4.0.5) [53]. For the above three classification algorithms, we tried the common kernel functions provided by the R package for each algorithm separately, and the rest of the parameters were the default.

Both training data and verification data are from field survey and 5-fold cross-validation tests were used to evaluate the classification performance. Cross-validation is a statistical method to assess the accuracy of modeling and prediction procedures [54]. K-fold cross-validation is widely used in machine learning [55,56]. Four-fold data had been used for training and 1-fold data for validation each time. It is noteworthy that each kind of rock is divided into different folds in equal proportion.

4.3 Feature reduction

To remove noise from the classification process and avoid overfitting, feature reduction is important. In reducing the features, we considered both the importance metrics of the classifiers and the pairwise separability.

For two kinds of dual-pol SAR data with 26 features, we chose the importance index provided by AdaBoost and RFs, because these two classifiers had better OA. The randomForest package provides two important metrics: Mean Decrease in Accuracy (MDA) and Mean Decrease in Gini. MDA was chosen for this study because it is widely considered to be the most effective variable importance measure for RFs [57]. Another variable selection method is the importance provided by the adabag package, which considers the gain of the Gini index given by the variables in the tree and the weights of the tree.

The above two variable selection methods were applied to rank the features provided by the two SAR data. Meanwhile, the J-M distance was used to express the separability between categories [58]. Higher values indicate stronger separability between samples. The results of J-M distance were also used as a reference for feature reduction [59]. Based on the importance metrics ranking results, we removed the bottom three covariates from the respective data. The most suitable variable selection method was determined by the OA of the classification before and after the feature reduction.

4.4 Accuracy evaluation

For accuracy assessment, we focused on analyzing the classification performance of the AdaBoost classifier with the highest OA. We also tested the difference in the OA of igneous rocks before and after using the texture component.

The importance of texture features was confirmed again in addition to the feature importance rank. The confusion matrix was used for the evaluation of the classification effectiveness of igneous rocks. It not only accurately describes the confusion status of each rock type, but also provides a measure of the OA and kappa coefficient. Precision (user’s accuracy), recall (producer’s accuracy), and F1-score were calculated simultaneously.

4.5 Oxide content test

Since the chemical composition of rocks indirectly affects the backscattering coefficient by influencing the dielectric constant and weathering characteristics of rocks, we measured the chemical composition (oxide content) of five granitoid samples as well as three syenite porphyry samples.

5.1 Characterization of backscattering coefficients and decomposition parameters

The VV and VH backscatter coefficients of igneous rocks were obtained from dual-pol C-band SAR data. Figure 3a represents the boxplots of VV and VH. The horizontal line in each box indicates the median of the data. The black circle is the mean value. Granitoid and tuff samples have similar mean VV values of −14.6 and −14.2 dB, respectively. Syenite porphyry has a slightly smaller mean VV around −15 dB. Compared to VV, the mean VH of the three types of samples are more different. Syenite porphyry has the smallest mean VH value (−28.8 dB) with a large standard deviation, and tuff has a relatively large mean (−26.6 dB).

Figure 3

Boxplots of backscatter coefficients for igneous rocks from (a) Sentinel-1 (VV and VH) and (b) PALSAR (HH and HV).

The HH and HV backscatter coefficients were obtained from dual-pol L-band SAR data. Figure 3b represents the boxplots of HH and HV. The mean HV values are similar for granitoid and syenite porphyry (−30.8 dB) and slightly larger for tuff at −30.5 dB. The mean HH value (−20.3 dB) and standard deviation for tuff samples are the largest for the three igneous rock types. Granitoid has the smallest mean HH value (−22.9 dB).

Figure 4 shows the false-color composite images of backscatter coefficients of two types of dual-pol SAR data, which produce different responses due to the different wavelengths and polarization methods. Taking the areas A, B, and C in the figure as examples, Sentinel-1 and PALSAR images have obvious differences in multiple outcrops. For example, some intrusive rocks are more visible in the PALSAR false-color composite image, and some water systems can be seen that are not visible in Sentinel-1, even though PALSAR has a relatively low spatial resolution (15 m). The water systems and rock orientation are more hierarchical in the Sentinel-1 false-color composite image.

Figure 4

RGB false-color composite images from (a) Sentinel-1 and (b) PALSAR backscatter coefficients.

For the parameters H and α extracted from Sentinel-1, granitoid and tuff samples have similar mean values, while syenite porphyry differs from them. The mean H of the granitoid is about 0.05 smaller than that of the other two types. The mean A value of granitoid and tuff samples is 0.8, while syenite porphyry is slightly larger at 0.9. The median α of the three types of igneous rock samples is almost the same, and the mean value of the syenite porphyry is 1.3° smaller than that of granitoid and tuff samples (Figure 5a).

Figure 5

Boxplots of H–α–A distributions for igneous rocks from (a) Sentinel-1 and (b) PALSAR. The black circle is the mean value.

Among the decomposition parameters extracted by PALSAR, the mean H values of granitoid and tuff are similar, and syenite porphyry is 0.03 larger than them. The mean α values of all three types of igneous rock samples are around 13° and syenite porphyry is slightly larger at 13.1°. The mean A values of the three types of samples are between 0.76 and 0.78. Similar to the boxplots of backscatter coefficients, H, A, and α are not sufficient to distinguish the three categories of igneous rocks, due to high intra-class variance (Figure 5b).

Figure 6a and b shows the distribution of igneous rock samples in the H–α planes of two dual-pol SAR data. In the H–α plane of Sentinel-1, the samples are concentrated in the space of 0.2 < H < 0.7 and 5° < α < 25°. As for the H–α plane of PALSAR data, the samples are concentrated in the space of 0.4 < H < 0.7 and 5° < α < 22.5°. The covariance matrix (C2) of the dual-pol SAR data contains only two eigenvalues, so the boundary of the H–α plane consists of two symmetrical curves about the α equals to 45°. Among them, the samples are more compactly distributed in the H–α plane of PALSAR data.

Figure 6

Distribution of igneous rock samples in the H–α plane of (a) Sentinel-1 and (b) PALSAR.

5.2 Feature reduction and evaluation of classification results

Figure 7 shows the top six important features from both data using MDA and Adaboost-based Gini index variables selection methods. It can be seen that texture features and backscatter coefficients are very important for igneous rock classification. As the values of the bottom-ranked importance indicators are close to zero, they are not shown.

Figure 7

Feature importance: (a) top six scores of C-band data features calculated by MDA, (b) top six scores of L-band data features calculated by MDA, (c) top six scores of C-band data features calculated by Gini index based on AdaBoost, and (d) top six scores of L-band data features calculated by Gini index based on AdaBoost.

Tables 2 and 3 show the parameters corresponding to the highest five J-M distances for each rock combination. J-M distance, MDA, and Gini index provide different results, but they have a certain pattern. For instance, we do not find any polarization decomposition parameters in Figure 7, Tables 2 and 3, either H, A, or α. However, texture parameters and backscatter coefficients appear frequently.

We removed the three features that were ranked low based on MDA and Gini index. According to MDA results, C_VV_Contrast, C_VH_Entropy, and C_VH_Correlation were removed from the C-band features, and L_Alpha, L_HV_Correlation, and L_HH_Homogeneity were removed from the L-band features. According to the Gini-based selection results, C_Anisotropy, C_VH_Entropy, and C_Entropy were removed from the C-band features, and L_HV_Dissimilarity, L_HH_Variance, and L_Alpha were removed from the L-band features.

We compared the OA before and after feature reduction to select the most promising classifier (Figure 8). Gini index feature selection method improved the classification accuracy of some classifiers, but MDA was better adapted for most classifiers. Ultimately, the RF-based MDA metric was used for variable selection, and for most classifiers, the MDA-based feature selection method improved the OA.

Figure 8

The overall accuracies of different classifiers of (a) Sentinel-1 and (b) PALSAR before and after feature reduction.

Ensemble learning algorithm maintained higher OA than other classifiers, whether feature selection or not. For both C-band and L-band SAR data, the AdaBoost algorithm has the highest classification accuracy of 0.5 and 0.42, respectively. In the following, we only analyzed the classification results of the AdaBoost classifier in detail. Using the AdaBoost classifier, the OA with and without texture features is shown in Figure 9. The OA of igneous rocks is improved by 2% for both C-band and L-band dual-pol data after adding texture features. The error bar highlights the advantage of 5-fold cross-validation.

Figure 9

Overall accuracies of cross-validation with and without texture features.

Regarding the confusion matrices (Tables 4 and 5), the classification is better using C-band SAR data (OA = 50%, kappa = 0.24) and slightly worse for L-band SAR data (OA = 42%, kappa = 0.13). The F1-score is defined as the summed average of the precision and recall (Table 6). For Sentinel-1, tuff has the highest F1-score of 0.55 and syenite porphyry has the lowest of 0.42. For PALSAR, the highest F1-score is also from tuff at 0.45 and the lowest F1-score is 0.38 for granitoid.

5.3 Oxide content test

Table 7 shows the mean values of oxide content for these two types of samples. The average SiO₂ content of the granitoid samples is 68.73%, which is consistent with the SiO₂ range of acidic igneous rocks. Granitoid samples have higher values of SiO₂, Fe₂O₃, CaO, and Na₂O. Syenite porphyry samples have higher values of TiO₂, Al₂O₃, FeO, MgO, K₂O, and P₂O₅.

6.1 Backscattering characteristics

Backscatter coefficients are sensitive to the moisture content, dielectric constant, and roughness of the feature. They are widely used in agriculture and environmental remote sensing [60,61]. In the domain of geology, backscatter coefficients are used to discriminate large intrusive rocks as well as baked edges because they have significantly different roughness or dielectric constants from the surrounding sedimentary or extrusive rocks [19]. However, for areas without significant intrusive bodies, backscatter coefficients are difficult to correlate directly with lithology [35]. For different types of igneous rocks, the differences in mineral composition are subtle and the surface morphology may be similar after long-term weathering.

The surface roughness is considered to be an important factor in the control of radar backscatter by the geological body [62,63], and the dielectric constant can also affect the backscatter intensity. According to Zhao [15], for neutral and acidic rocks, the content of CaO, MgO, Fe₂O₃, and Na₂O correlate most strongly and positively with the dielectric constant. Table 7 shows that granitoid samples have relatively high levels of CaO, Fe₂O₃, Na₂O, which may contribute to the high dielectric constants of granitoid samples. The black mica in tuff is oxidized during weathering, in which divalent iron is oxidized to trivalent iron, which may change the surface roughness and thus affect backscatter coefficients.

Since it is difficult to discriminate lithology directly by visual interpretation on SAR images, we compared other large-scale features, such as rock formation. They show that different wavelengths of SAR data contain different information. The combination of multi-band backscatter coefficients may provide greater support for lithology classification.

6.2 Polarization decomposition

The polarization decomposition parameters do not seem to play a useful role in this study. The three types of igneous rock samples are compactly distributed in the H–α plane, especially PALSAR data. The H–α decomposition can provide support for lithological classification [26], but the small contribution in this study is caused by the limitation of the polarization mode of the dual-pol SAR data. The full-polarization H–α plane is divided into nine regions based on the scattering mode, which can better distinguish the scattering characteristics of ground objects [38]. Compared with HH–VV polarized SAR data, VV–VH and HH–HV polarized SAR data are more difficult to distinguish various scattering mechanisms [64]. Ji and Wu [64] found that a large amount of low-entropy dipole scatter and low-entropy multiple scattering is transferred to low-entropy surface scatter in the VV–VH polarized H–α plane, and low entropy multiple scatter is transferred to low entropy surface scatter in the HH–HV plane. Some of the low entropy surface scatter diffuses into medium entropy surface scatter. Different types of igneous rocks tend to have different scattering characteristics due to different formation causes and formation times, but in this study, almost all samples are in the region of low entropy surface scattering. Moreover, α extracted in the L-band ranks low in the feature importance ranking. The above results indicate that the VV–VH or HH–HV dual-pol H–α planes are not good classification criteria in the rock exposure region. Due to the limitation of the satellite sensor, we did not acquire full-polarization SAR data, and the applicability of the full-polarization H–α plane in lithology classification may be explored in the subsequent study.

6.3 Texture features

Even in high-resolution radar images, lithology does not correlate directly with all texture features. Feature selection helps us select appropriate texture features in the classification process [65]. Intrusive rocks are formed deep inside the surface of the earth whereas extrusive rocks are formed at the surface of the earth when magma finds a way to eject or pour out of the surface. Different formation processes may generate different textural features, which may be the reason why textural features play an important role in classification.

6.4 Lithology classification and accuracy evaluation

The accuracy of igneous rock classification obtained based on dual-pol SAR parameters varies depending on (i) the classification algorithm (i.e., various types of machine learning algorithms) and (ii) the SAR dataset (i.e., Sentinel-1 or PALSAR). Feature reduction methods may be useful only for a certain classifier, but after experimentation, appropriate feature filtering methods can improve classification accuracy and reduce the risk of overfitting. Based on the importance ranking provided by the classifier, we refer to the J-M distance. Because of the similarity of the results, we finally followed the importance ranking of the classifier. It is worth noting that the AdaBoost algorithm has the best classification results for both types of dual-pol SAR data, although the training time of the model is slightly longer. It might be an interesting direction to explore the role of different sub-classifiers in the AdaBoost algorithm if they can be constructed for lithology classification using SAR data.

The spatial resolutions of the two dual-pol SAR data we chose are 10 m for Sentinel-1 and 15 m for PALSAR. Since the samples are chosen with an outcrop area as big as possible, the influence of resolution on the lithological discrimination can be ignored. Higher resolution L-band SAR data needs to be acquired in future research and this is an important factor to improve the lithology classification accuracy.

Granitoid tend to have better accuracy in lithology classification using optical remote sensing data because minerals such as quartz, which is widely distributed in granitoid, have special spectral characteristics [35,66]. The reason for the low accuracy of granitoid samples in lithology discrimination using SAR data may be that granitoid bodies are affected by weathering. The dielectric constant or surface roughness of weathered layers is similar to the surrounding rocks or sediments. Another possibility is that the radar echoes are affected by the special rock stratification phenomenon in the study area. The borehole data provided by the Natural Geological Archives Data Center provided us with a new idea [67]. The borehole data obtained near the study area in 2009 showed that granitoid with a depth of about 8 m was developed at the surface, and gabbro was developed at the bottom of the granitoid. Several boreholes suggest that the surface was covered with shallow Quaternary sediments. SAR data in L-band could potentially penetrate to a depth of about 1 m, and another rock at a depth of 8 m would not be detected [68,69]. Where granitoid layers are thin, SAR may detect other types of rocks. The penetrating nature of SAR may cause lithological confusion, but it has a penetrating effect on surface gravel and rock weathering surfaces, which is the advantage of SAR over optical remote sensing.

SAR images have advantages in the case of similar spectral features. Lu et al. [35] achieved acceptable results in the identification of dolomite and limestone using Sentinel-1 extracted features combined with SRTM data. They confirmed that Sentinel-1 is effective in the discrimination of sedimentary rocks and that the backscatter coefficients, decomposition parameters, texture, and elevation each play an important role. In this study, we focus our research target on igneous rocks. The relationship between lithology and elevation was found through field investigations, and the reason for not using topographic data was the desire to classify by the parameters provided by the dual-pol SAR data. In our study, the contribution of decomposition parameters is very small. This may indicate that decomposition parameters are effective in the discrimination of sedimentary rocks, but ineffective in the classification of igneous rocks.

More approaches will be applied in the future to improve classification accuracy. The roughness may vary within a pixel due to the limitation of SAR image resolution. Since different bands have different sensitivity to surface roughness, combining multi-frequency SAR data can provide more information for lithology classification [25,62]. Meanwhile, applying the idea of time-series to lithology discrimination might be a good idea to improve the quality of classification features.

Overall, this article proves the effectiveness of classifying igneous rocks using only dual-pol SAR data and providing new ideas for prospecting and large-scale geological mapping.