Lithostratigraphic Classification Method Combining Optimal Texture Window Size Selection and Test Sample Purification Using Landsat 8 OLI Data

2.1 Workflow of Lithostratigraphic Classification

The workflow for this method comprises several stages shown in Figure 1. The main stages are: semivariogram based optimal texture window size selection, GLCM based textural features extraction, SVM based lithostratigraphic classification, test sample set size calculation and selection, sample purification rules creation based on textural information of training samples and purified test samples, accuracy assessment based on test samples with sample purification, and results analysis.

Figure 1

Flowchart of the lithostratigraphic classification method combining optimal texture window size selection and test sample purification technique

2.2 Semivariogram-based Window Size Selection for Texture Extraction

The semivariogram proposed by Matheron [47] is a spatial statistical method. The semivariogram can show how the variance of a variable varies with different sampling intervals on one image, which can represent the degree or scale of spatial variability. Semivariance γ(h) is determined by Equation (1) [21].

where h is the lag (or spatial sampling interval) and is a vector in both distance and direction because of the existence of anisotropy. X( i) is a sample value at location i, and X(i+h) is another sample value separated from X(i) with the sample interval h. N(h) is the total amount of pairs of X(i) and X(i+h).

The increasing and second derivative of semivariance are calculated from Equation (2) [21] and (3) [21] respectively.

where [Δγ(h)_i] is the increasing of semivariance and [γ''(h)_i] is the second derivative of semivariance, i represents the calculating times of the semivariogram with different sampling intervals.

The range of semivariogram (threshold of autocorrelation) associates the feature scale of the landscape. Within the range, the autocorrelation exists among variables, while beyond the range, autocorrelation may be untenable. So the range can be deemed as the measurement of similarity between variables, and it can indicate the size of a spatial object, a spatial phenomenon or a spatial pattern, as well as the size of texture primitive. So the range of semivariogram can be used to estimate the optimal texture window size. More specifically, the range is set as the sampling interval h_i corresponding to [Δγ(h)_i] whose value is first less than or equal to zero with the increasing of h and [γ"(h)_i] whose value is first greater than or equal to zero with the increasing of h. More details, please refer to Section 4.1.

2.3 Principles of Textural Feature Expression Based on GLCM

GLCMis a common method for describing texture by studying the spatial correlation of pixel gray. As a statistical method, it has been widely used to extract the features of texture statistics from images. Based on GLCM, Haralick [13] proposes 14 textural features to extract from the texture statistics of images. Commonly used features for lithostratigraphic classification include Angular Second Moment (ASM), Homogeneity (Homo), Entropy (Ent), Correlation (Cor), and Contrast (Con) [48, 49]. In this study, taking into account the time consumption, four important features, ASM, Ent, Cor, and Con, are selected for extracting textural features of the study area.

ASM, known as Uniformity or Energy, is the sum of squares of the elements value in the GLCM [48]. It reflects gray scale distribution uniformity and texture roughness of the image [50]. ASM is high when the image is a uniform and regularly changing texture pattern. ASM is extracted by Equation (4) [51].

Ent, determined by Equation (5) [48], shows the complexity of the image [50]. The higher the Ent, the more complex image [48].

Cor, calculated by Equation (6) [48], measures the dependency of grey levels in the row or column direction [15]. Therefore, it reflects the correlation between local grey levels [52]. Cor is high when grey levels of an image are even [50].

Con, extracted by Equation (7) [51], measures the distribution of values in the GLCM and local changes in the image [15]. It reflects the clarity of the image and the depth of the textural grooves [53]. Con is high when the textural groove is deep, and the image is clear; the Con is low when the textural groove is shallow, and the image is fuzzy [50].

Where P(i,j|d,θ) is a GLCM with step size d and angle θ, (i, j) means position. μ_x, μ_y, σ_x, and σ_y are the means and standard deviations of rows and columns.

For this paper, the step sizes of one and direction of 45∘ were set for computing the P(i,j|d,θ). The window size is selected by semivarogram. Then, according to Equations (4)–(7), ASM, Ent, Cor, and Con can be extracted.

2.4 Support Vector Machine Classifier

The SVM studied by Vapnik et al. [54, 55, 56, 57] is a supervised machine learning algorithm which is based on statistical learning theory. The SVM seeks a hyperplane that satisfies the classification requirements, and ensures the distance between the training data and hyperplane is as far as possible. The learning strategy of SVM is through interval maximization, which is applied to convex quadratic programming problem solving.

The rationale of SVM can be summarized as follows [58]:

1. Map the points in low-dimensional space into high-dimensional space based on the kernel function, so that they become linearly separable.

2. Use the principle of linear division to determine the optimal hyperplane which can separate the two classes correctly, and the classification margin is at the maximum. The classification margin is the gap between the hyperplane and the closest training samples. The closest training samples are support vectors.

The optimal classification function is [58]:

where {(X₁, y₁),(X₂, y₂),...,(X_n, y_n)} are training samples, X_i∈ Rⁿ, y_i∈ (-1,1). X = (X₁, X₂,..., X_n) means the input vectors, X_i(i = 1,2,..., N) are support vectors. a_i and b are undetermined coefficients of the hyperplane. K(X, X_i) is the kernel function including Linear, Polynomial, Radial Basis Function (RBF), and Sigmoid Kernel. The Radial Basis Function is chosen for this study.

2.5 Accuracy Assessment Combining Test Sample Purification

3.1 Study Area

The study area (Figure 2) is located in the northwest margin of the Tarim Basin, which lies in the Xinjiang province of Western China, with an eastern longitude range from 77∘35' to 78∘10' and a northern latitude range from 40∘0' to 40∘25'. The average elevation of the study area is approximately 1885 m. The region has a temperate-continental climate, resulting in low precipitation, poor vegetation development, and good rock outcrops, which is the ideal area for using satellite remote sensing technology for geological research. The main stratigraphic units in this area include Є-O, S, D, C-P, N₂, Q₁, Q₃₋₄ [60]. The 1:200,000 geological map (Xinjiang Geological Survey, 1965a, b) of the area was used as the reference map [60], shown in Figure 3. Main rock units in the study area and their lithostratigraphic description are shown in Table 1.

Figure 2

Location of study area.

Figure 3

The lithostratigraphic map of the study area.

3.2 Data and Data Pre-processing

The image used in the present study was acquired by Landsat 8 OLI on July 30, 2013. The data (LC81480322013211LGN00), was provided by USGS EROS data center [46]. The image has 0.75% of cloud cover, which mainly concentrated in the upper corner. Therefore, no cloud was presented in the image of study area, as shown in Figure 4. Figure 4 shows the false color composition of the study area. The study uses the multispectral and panchromatic data of the image, which are detailed in Table 2. The images we obtained are Level-1 digital number data products that have been radiometrically and geometrically corrected by USGS. The images were processed as follows:

Figure 4

Landsat 8 false color composite image of the study area (R = band 5 (NIR, 0.845–0.885), G = band 4 (Red, 0.630–0.680), and B = band 3 (Green, 0.525–0.600)).

1. Radiometric calibration was conducted for VNIR–SWIR images to acquire the top of atmosphere reflectance data.

2. FLAASH atmospheric correction with the Mid-Latitude Summer atmospheric model and the rural aerosol model were conducted for the top of atmosphere reflectance data to acquire the surface reflectance data.

3. For ease of visually interpretation and enhancing the separability between different category, the surface reflectance data were fused with the panchromatic band by Gram-Schmidt Pan Sharpening.

4. C-correction, a topographic correction method based on the linear relationship between the image pixels and the cosine of the solar relative incidence angle [61], was conducted for the images to remove the shadows effect. ASTGTM DEM (Figure 5) with 30 m spatial resolution was used in Topographic Correction process.

Figure 5

ASTGTM DEM of the study area.

After data pre-processing, the fused VNIR–SWIR bands were used to extract texture features. All the procedures of data processing in this section were conducted in ENVI 5.1.

3.3 Vegetation Elimination in Remote Sensing Images

A small amount of vegetation was present in the used image, which needs to be considered, as shown in Figure 4. In this figure, the vegetation is shown in magenta. The texture features are described based on the spatial variability associated with the digital number (DN) of the local pixels. Therefore, the vegetation in the image will hide or damage the texture information. In order to remove vegetation, Equation (14) [62, 63] was used to calculate the normalized difference vegetation index (NDVI). In this paper, the NDVI threshold is set to 0.3, marking off the vegetation pixels, which in Figure 6 is green.

Figure 6

Vegetation pixels in green with 0.3 < normalized difference vegetation index (NDVI) < 1.0.

4.1 Window Size Selection for Texture Extraction

As shown in Figure 7(a), the representative sampling area with a size of 600 × 600 pixels is used for calculating semivariance in this experiment. In order to accurately reflect the variation of semivariance along the spatial sampling distance, the lags are set from 3 to 45, and the step of the lag is set as 2. Band 6 with the largest standard deviation, contains the maximum amount of information. To reduce the amount of computation, the semivariograms in the horizontal and vertical directions are calculated by band 6 of the sampling area images. The synthetic semivariances are the average of the horizontal and vertical semivariograms, displayed in Figure 7(b).

Figure 7

(a) the sampling area with a size of 600 × 600 pixels; and (b) synthetic semivariograms.

From Figure 7(b), it can be seen that the range of the synthetic semivariogram is not obvious. To achieve pre-estimation of the reasonable window size for texture extraction, the increasing of the synthetic semivariance [Δγ(h)_i] and second derivative of the synthetic semivariance [γ''(h)_i] are displayed in Figure 8(a) and (b). According to Figure 8, the sampling distance for the image shown in Figure 8(a) is 11, corresponding to a 23 × 23 spatial window, and the sampling distance for image shown in Figure 8(b) is 15, corresponding to a 31 × 31 spatial window. By employing the texture window size pre-estimation method proposed in Section 2.2, the theoretical optimal texture window size is determined as 31 × 31.

Figure 8

Computed results of: (a) the increase of the synthetic semivariance; (b) the second derivative of the synthetic semivariance.

4.2 Spatial Textural Feature Analysis

4.2.1 Spatial Textural Feature Extraction

According to Figure 8, the texture extraction window size of 31 × 31 is used for extracting Ent, Cor, and Con. According to seven bands of Landsat 8 data, four textural features of GLCM, ASM, Ent, Cor, and Con, were calculated to obtain 28 texture images. In order to study the effectiveness of the different texture images for lithostratigraphic identification and create sample purification rules, the texture values of the training samples were extracted from the texture images. The median texture value of each category is shown in Figure 9.

Figure 9

Median texture value for the seven rock types: (a) median ASM value; (b) median Ent value; (c) median Cor value; and (d) median Con value. TFi means the texture images were derived from spectral band i.

Figure 9 provides the following information:

1. The ASM feature distinguishes Q₃₋₄ from other rock types well, and the value calculated from band 5 is largest. Q₃₋₄ is generally larger on ASM values, which may be because its texture is more uniform and regular-changing than other rock types.

2. N₂ has the highest Ent value, which indicates its surface is most complex. Q₃₋₄ has the lowest Ent, TF5, and TF6 values, which means it is simpler than other rocks.

3. The Cor feature can separate Q₁from other rocks. The Cor values of Q₁ are relatively small, which means Q₁ has a small local correlation on the image. Cor values curves of S and Q₃₋₄ fluctuate greatly, which may due to the fact that texture images derived from different bands have a variable effect on rock identification [19].

4. N₂ shows a significant difference in the Con feature compared to other rock types, which may result from its deep textural trench. From TF1 to TF7, between the Con values, there is an obvious difference between Є-O and C-P, and the between S, D, and Q₁.

In short, the application of GLCM textural features to the classification process can help improve the separation of different rock types. For the same texture, the texture images obtained from different bands have different sensitivities for rock identification. The different textures calculated in the same band have different sensitivities for lithostratigraphic identification as well. For example, an ASM image derived from band 5 distinguishes Q₃₋₄ from other rock types, however, Ent and Con images calculated from band 5 separate N₂ well. Overall, different texture images have abilities to identify different rocks.

4.3 Lithostratigraphic Classification

There are seven rock types in the study area including Є-O, S, D, C-P, N₂, Q₁, and Q₃₋₄, detailed in Table 1. The training sample amounts of each rock type are displayed in Table 6.

The training samples of these rocks were selected based mainly on the reference map. Based on GLCM, four features, including ASM, Ent, Cor, and Con, were extracted. In this paper, a SVM classifier was used to classify the combinations of texture and spectral images. The combinations of texture and spectral images classified by SVM are: spectral features (SF), SF+ASM, SF+Ent, SF+Cor, SF+Con, SF+ Ent+Con, SF+Ent+Cor+Con and SF+ASM+Ent+Cor+Con, as shown in Table 5.

4.4.1 Calculation of Sample Set Size

Before application, classification maps constructed from images should be assessed for accuracy [34]. An accuracy assessment process begins with sampling design. The first step in sample design is to calculate sample set size for the test samples. With z = 1.96, O = 0.74, and d = ±5 %, the resulting sample set size from Equation (9) is n = 295, marked as sample A1. Dependent on sample A1, the first accuracy assessment is carried out, and the mean user’s accuracies (U_i) are shown in Table 7.

Because N is adequately large (exceeding six million pixels in this study), for Equation (10), the second item in the denominator can be neglected. Once n is defined, proportional allocation [65] can be used to determine the sample size per strata in this paper. In addition, the stratified sampling [42, 66, 67, 68] strategy is employed and systematic sampling [36, 59] is then applied for every class. In this research, test samples select from the entire original image with step sizes of 80 pixels row by row and column by column. Table 7 provides information that is needed to calculate sample set size from Equation (10). Among them, the mapped area proportions (W_i) is the mean of all classification results. Based on Table 7, specified S(Ô) for 0.01, the sample set size calculated from Equation (10) is n = 1246, marked as sample B1. n is the overall sample set size. The allocation of the number of samples per class was determined by proportion, detailed in Table 8. When we applied sample purification rules (Table 9) to samples A1 and B1, we derived sample A2 and B2, detailed in Table 8.

4.4.2 Creation of Sample Purification Rules

According to problems that were mentioned in Section 1, test sample purification was implemented to improve the reliability of classification accuracy. Based on Section 2.5.2, Section 4.2, and ranges of ASM, Ent, Cor, and Con values, sample purification rules were created, detailed in Table 9.

5.1 Verification of Texture Window Size Selection Method Based on Semivariogram

To verify the effectiveness of the semivariogram-based texture extraction window selection method, different ASMs with different window sizes combined with spectral features are classified by SVM. To reduce computation of verification, ASM extractions are based on window size from 5 × 5 to 41 × 41 with a step of 6. To maintain the comparability between different classification results, the training and test sample sets were the same. The test sample set used is sample B2. Table 10 demonstrates the overall accuracies (OA) and kappa with different ASM window sizes in SVM based supervised classifications.

Table 10 shows the change in OA and Kappa with the change of ASM extraction window in SVM classifications. The highest OA is 88.09% corresponding to the texture extraction window of 23 × 23, which is just one of the pre-estimated window sizes. The second highest OA is 88.01% corresponding to the texture extraction window of 29 × 29, which is very close to the pre-estimated window of 31 × 31.

The experiment results verified that the spatial statistics-based texture window size pre-estimation is effective.

5.2 Benefit of Sample Purification

Figure 10 shows the OA and kappa of different combinations derived from test samples A1, A2, B1, and B2. As shown in the figure, the OA and Kappa obtained from the four test samples are different for the same classification result. The differences between samples A1 and B1, or samples A2 and B2 may result from the difference of sample set size and sampling method. For accuracies derived from samples A1 and A2, or samples B1 and B2, the difference may due to sample purification. The OA and Kappa derived from test samples with sample purification are higher than original test samples, because the error in the test samples that have been purified is smaller than those that have not. The error of test sample will reduce classification accuracy. After sample purification, the error of samples A2 and B2 is smaller than original test samples. Accuracies based on samples A2 and B2, which reduced the error of test samples, is more precise than that based on original test samples.

Figure 10

The overall accuracies (OA) and kappa of different combinations derived from test samples A1, A2, B1, and B2: (a) OA curves; (b) kappa curves. Sample A1 and B1 are the original sample set selected by simple random sampling and stratified sampling, sample A2 and B2 are the sample set with sample purification selected by simple random sampling and stratified sampling. Ci is the different combinations classified by SVM, detailed in Table 5.

Table 11 demonstrates the PA and UA of different classifications obtained from samples B1. Table 12 shows the PA and UA of different classifications obtained from samples B2. By analyzing Table 11 and 12, it can be explicitly seen that UA and PA for lithostratigraphic classification have been improved when the ASM feature is added. Regardless of whether the texture features were added or not, the recognition accuracies of S and D is relatively high. This may be due to their prominent color features. When combined with texture features, the recognition accuracy of Є-O and Q₃₋₄ is significantly improved. It may result from their different s urface roughness, which would assist rock. For C-P and N₂, their texture features are similar to Q₁, and the spectral characteristics of C-P are similar to Є-O, resulted in their poor classification accuracy. Furthermore, besides UA for N₂ and Q₁, the rest of the features were also improved after the test samples were purified. In summary, lithostratigraphic classification combined with test sample purification i s more r eliable for lithostratigraphic identification.

5.3 Effect of Texture features for Lithostratigraphic Classification

According to Figures 11 and 12, it also can be seen that lithostratigraphic classification accuracy has improved when spectral features are combined with textural features, except when the Cor feature is used. We concluded that when the selected suitable textural features are fused with spectral features classification accuracy is improved. In this study area, ASM, Ent, and Con are suitable textural features and Cor is an unsuitable textural feature for lithostratigraphic classification. Accuracy of SF + Ent+ Con is higher than SF+ Ent+ Cor+ Con and SF+ ASM+ Ent+ Con+ Cor, and it may cause by possible conflicts between ASM, Ent, Con, and Cor.

Figure 11

Classification results of different combinations classified by SVM; Ci means different combinations classified by SVM, detailed in Table 5.

The classification results of only spectral features and spectral features fused with different textures are shown in Figure 11. The error of classification mainly occurs in the boundary of each rock, which is not random. Figure 11C2 is the classification result of SF+ ASM, which is closest to the actual lithostratigraphic distribution. Furthermore, Figure 11C3 and C6 which was respectively the classification results of SF+Ent and SF+ Ent+ Con, is also very close to the reference map. Figure 11C4, C7, and C8, the classification result of SF+ Cor, SF+ Ent+ Cor, and SF+ Ent+ Cor+ Con, have serious salt and pepper effect. It may result from that Cor feature is not conducive to the separation of rocks.

Figure 12 shows the local comparison with reference map, the classification results of spectral features, and of SF+ ASM. In classification results of spectral features and SF+ ASM, recognition of N₂ and C-P are poor. Identification of N₂ and C-P have been improved, and C-P misclassified as Q₁ is also slightly reduced when ASM feature is added. Moreover, recognition of Q₃₋₄ and Є-O are significantly improved, and Q₃₋₄misclassified into Q₁, N₂, and D is notably reduced, and the internal salt and pepper effect is also obviously reduced.

Figure 12

Zoom-in and comparison of reference map and classification results of: (a) reference map; (b) SVM with only spectral features (SF); and (c) SVM with spectral features fused with ASM (SF+ ASM).

Improved accuracy is obtained when the spectral features and suitable texture feature are integrated in SVM for lithological classification. The accuracy is highest when using SF+ ASM for lithological classification. This can be attributed to two factors. The first factor is that the texture features keeps the separability of objects, as they reveal structure information of objects and their relationship with the surrounding environment in image. The second factor is that there is a large difference in the uniformity of gray distribution or texture thickness between different rocks in the study area, so that SF+ ASM yields the highest classification accuracy.