Introduction

Land cover classification is a critical process in remote sensing, providing spatially explicit information at different scales for numerous environmental applications [1]. Such information is widely applied to issues that require practical geospatial solutions like land cover change detection [2], environmental monitoring [3,4], fossil fuel exploration [5], and landmine detection [6]. Land classification techniques, likewise, play a crucial role in the integration of Earth observation data into comprehensive interdisciplinary issues on sustainable development goals achieving [7,8], in particular, combat climate change and its impacts [9], reverse land degradation [10] and halt biodiversity loss [11], protect water-related ecosystem for safety water supply [12], and provide support for food security and sustainable agriculture [13,14].

Today, most classification methods are mainly divided into supervised and unsupervised [15]. However, in remote sensing, the supervised classification methods are the most appropriate for the majority of the thematic tasks because applying these methods can establish the characteristics of the output classes, unlike unsupervised ones. A training sample set is used to set the characteristics of classes in supervised classification methods. Such a set contains the signatures of features of each class.

The input data for classification is heterogeneous geospatial data, which can be represented in the form of raster layers. To combine such layers into a single array, it is customary to form a geospatial data cube [16]. From ordinary datasets, data cubes differ by integrating different data types into a coherent and interoperable structure [17,18]. After the cube's formation, the training sample's signatures must be determined in each layer. Hence, each layer is a feature of the training sample.

A geospatial data cube for specific classification tasks can comprise numerous layers. Including multitemporal data, a data cube aims to distinguish dynamic objects that change significantly during a specific period or vary much from each other at different stages of their development, e.g. vegetation cover, wetlands in particular [12]. Different physical aspects of multisource data, e.g. optical and radar, highlight diverse object traits and variations in land cover types [19]. Multiple ancillary data, e.g. data on geomorphology, hydrology, or phenology, help to differentiate land cover types due to their context [20].

However, the redundancy of the geospatial data cube causes two significant problems [21]. Firstly, the processing of such a data cube has high computational complexity. Secondly, since the signatures of the training sample are defined in each layer of the geospatial data cube, the separability of the training sample can be reduced if the layers are either incorrectly created or irrelevant to the selected thematic task. In turn, the low separability of the training sample leads to a decrease in classification accuracy [22].

In light of the above, optimization of the geospatial data cube is seen as a solution to the mentioned problems [23,24]. Among approaches to reduce input data, the Principal Component Analysis (PCA) [25] and the Minimum Noise Fraction (MNF) [26] are the most widely used. There are also similar methods, for example, Noise-Adjusted Principal Components (NAPC) [27], Independent Component Analysis (ICA) [28], Non-Negative Matrix Factorization (NMF) [29] and Spatio-Spectral Decomposition (SSD) [30]. However, a common disadvantage of the considered approaches is that they do not consider the training sample's structure (in particular separability) and the selected classifier's specificity.

The presented study aims to enhance land cover classification by selecting the cube layers, the training sample separability of which will be the highest among other options. For this purpose, the optimization technique of the geospatial data cube was developed. It has two goals: enhancement of land cover classification and reduction of geospatial data cube size.

Hence, in the relevant sections of this article, the separability assessment of the training sample, the optimization technique of the geospatial data cube, and the experiment conducted to demonstrate the effectiveness of the developed technique are described.

Methods

This section presents the optimization technique of a geospatial data cube. Since this optimization is based on training sample separability, the objection function is the developed separability index of the training sample (SITS). Thus, the training sample separability assessment is also presented below as an algorithm for SITS calculation.

Assessment of the training sample separability

Separability is one of the training sample characteristics that affect classification accuracy. This characteristic shows the extent to which signatures representing different classes do not overlap. A low degree of separability is inherent in a high level of training sample mixing. In turn, this leads to a significant number of misclassified objects in the classification. Thus, the training sample separability is directly proportional to the classification accuracy.

The algorithm depicted in the flowchart(Figure 1) describes the separability assessment of the training sample. The first step implies classifier training by the training sample. Importantly, the supervised classification method must be the same as the one selected to classify the geospatial data cube further. Moreover, due to the proposed, separability depends on its structure (i.e., the set of layers) and the selected supervised classification method.

In the second step, the classifier is used to classify each signature from the training sample set. The third step is the formation of the confusion matrix [31] for the classification obtained in the previous step.

The fourth and final step is calculating the SITS. This index quantifies the separability of training samples by measuring the ratio of correctly classified training samples to the total number of training samples. In other words, SITS equals the overall accuracy [31] based on the confusion matrix obtained in the previous step. The calculation of the SITS is shown in the following formula:

𝑆𝐼𝑇𝑆 = ∑ 𝑥 𝑖𝑖 𝐾 𝑖=1 𝑁 , (1)

where 𝐾 is the number of classes, N is the total number of training sample signatures, 𝑥 𝑖𝑖 is the number of class i signatures classified as class i (i.e. diagonal elements of the obtained confusion matrix that correspond to correctly classified signatures).

The values of the considered index range from 0 to 1. In this case, the value 0 shows that the training sample is entirely mixed (minimum separability), and the value 1 corresponds to the training sample, which is entirely separable (maximum separability).

Optimization technique

This technique is an optimization procedure that aims to reduce the number of layers of the geospatial data cube and increase the separability of the training sample, the signatures of which are defined in each layer of this cube. The objective function implies using the SITS. Thus, geospatial data cube optimization can be described as a search of the minor number of cube layers for which the training sample has the highest SITS value among all other sets of cube layers.

The flowchart of the technique algorithm is shown in Figure 2.

Figure 2: Algorithm of the optimization technique

The initial data and their characteristics will be introduced below for a detailed description of the technique algorithm.

Let the initial geospatial data cube have the following form: 𝐺𝐶𝐷 𝑖𝑛𝑖𝑡𝑖𝑎𝑙 = {𝐿 1 , 𝐿 2, … , 𝐿 𝑁 }, where 𝐿 𝑖 is the layer i of the geospatial data cube, and N is the total number of layers included in the initial geospatial data cube.

Then, as shown in Figure 2, the SITS value of the initial training sample, which has the signatures defined in each layer of the initial geospatial data cube, is first calculated. Let this value be 𝑆𝐼𝑇𝑆 𝑖𝑛𝑖𝑡𝑖𝑎𝑙 .

Let us assign the values obtained above to the variables corresponding to the optimal set of layers of the geospatial data cube and the corresponding value of the SITS. Thus, we have:

𝐺𝐶𝐷 𝑜𝑝𝑡𝑖𝑚𝑎𝑙 ≔ 𝐺𝐶𝐷 𝑖𝑛𝑖𝑡𝑖𝑎𝑙 , 𝑆𝐼𝑇𝑆 𝑜𝑝𝑡𝑖𝑚𝑎𝑙 ≔ 𝑆𝐼𝑇𝑆 𝑖𝑛𝑖𝑡𝑖𝑎𝑙 . Next, an iterative procedure follows, in which the following steps are performed at each iteration.

Step 1. At the iteration i, the current geospatial data cube 𝐺𝐶𝐷 𝑜𝑝𝑡𝑖𝑚𝑎𝑙 , consisting of N-(i-1) layers, is decomposed into N-(i-1) cubes. Each newly formed cube is obtained by discarding one of the layers from the current cube. Then, each of the newly created geospatial data cubes will have the following form:

𝐺𝐶𝐷 1 ′ = 𝐺𝐶𝐷 − {𝐿 1 } = {𝐿 2 , 𝐿 3, … , 𝐿 𝑁−(𝑖−1) }, 𝐺𝐶𝐷 2 ′ = 𝐺𝐶𝐷 − {𝐿 2 } = {𝐿 1 , 𝐿 3 … , 𝐿 𝑁−(𝑖−1) }, … 𝐺𝐶𝐷 𝑁−(𝑖−1) ′ = 𝐺𝐶𝐷 − {𝐿 𝑁−(𝑖−1) } = {𝐿 1 , … , 𝐿 𝑁−(𝑖−1)−1 }.

Since the number of layers decreases by one at each iteration, the obtained cubes will contain N-i layers. Therefore, the following is valid:

|𝐺𝐶𝐷 1 ′ | = |𝐺𝐶𝐷 2 ′ | = ⋯ = |𝐺𝐶𝐷 𝑁−(𝑖−1) ′ | = 𝑁 − 𝑖.

Thus, the generated cubes can be written in the form of the following set:

𝐶 = {𝐺𝐶𝐷 1 ′ , 𝐺𝐶𝐷 2 ′ , … , 𝐺𝐶𝐷 𝑁−(𝑖−1) ′ }.

Step 2. For each newly formed cube, the SITS value is calculated for the training sample, the signatures of which are defined in each cube layer. Then, the value of SITS for a particular cube 𝐺𝐶𝐷 𝑡 ′ will be denoted as 𝑆𝐼𝑇𝑆 𝑡 ′ . Thus, a set containing the value of the SITS for each newly formed cube will be obtained:

𝑆 = {𝑆𝐼𝑇𝑆 1, ′ 𝑆𝐼𝑇𝑆 2 ′ , … , 𝑆𝐼𝑇𝑆 𝑁−(𝑖−1) ′ }.

Step 3. Among the obtained cubes, the one with the highest value of the SITS is selected. Such a cube will be denoted as 𝐺𝐶𝐷 𝑖 . The selected cube can be expressed as follows:

𝐺𝐶𝐷 𝑖 = {𝐺𝐶𝐷 𝑡 ′ ∈ 𝐶|𝑆𝐼𝑇𝑆 𝑡 ′ = max{𝑆}}. 𝑆𝐼𝑇𝑆 𝑖 denotes the value of the SITS of the geospatial data cube 𝐺𝐶𝐷 𝑖 .

Step 4. The values of the variables 𝑆𝐼𝑇𝑆 𝑖 and 𝑆𝐼𝑇𝑆 𝑜𝑝𝑡𝑖𝑚𝑎𝑙 are compared, and two options are considered:

1) if 𝑆𝐼𝑇𝑆 𝑖 < 𝑆𝐼𝑇𝑆 𝑜𝑝𝑡𝑖𝑚𝑎𝑙 , then the execution of the optimization algorithm will be interrupted, and further steps will be ignored. The optimal geospatial data cube will be the one obtained in the previous iteration, namely 𝐺𝐶𝐷 𝑜𝑝𝑡𝑖𝑚𝑎𝑙 . Accordingly, the SITS of the training sample of the corresponding cube has the value 𝑆𝐼𝑇𝑆 𝑜𝑝𝑡𝑖𝑚𝑎𝑙 .

2) if 𝑆𝐼𝑇𝑆 𝑖 ≥ 𝑆𝐼𝑇𝑆 𝑜𝑝𝑡𝑖𝑚𝑎𝑙 , then the variable 𝑆𝐼𝑇𝑆 𝑜𝑝𝑡𝑖𝑚𝑎𝑙 is assigned the value of the variable 𝑆𝐼𝑇𝑆 𝑖 , i.e.:

𝑆𝐼𝑇𝑆 𝑖 ≔ 𝑆𝐼𝑇𝑆 𝑜𝑝𝑡𝑖𝑚𝑎𝑙 .

Step 5. At this step, as at the previous one, two options are considered: 1) if the number of layers of the obtained cube 𝐺𝐶𝐷 𝑜𝑝𝑡𝑖𝑚𝑎𝑙 is 1, i.e.:

|𝐺𝐶𝐷 𝑜𝑝𝑡𝑖𝑚𝑎𝑙 | = 1, then the execution of the optimization algorithm will be interrupted. The optimal geospatial data cube will be the one obtained at the current iteration, namely -𝐺𝐶𝐷 𝑜𝑝𝑡𝑖𝑚𝑎𝑙 . Accordingly, the SITS value of the training sample of the corresponding cube is -𝑆𝐼𝑇𝑆 𝑜𝑝𝑡𝑖𝑚𝑎𝑙 .

2) if the number of layers of the obtained geospatial data cube 𝑮𝑪𝑫 𝒐𝒑𝒕𝒊𝒎𝒂𝒍 is greater than 1, i.e.:

|𝐺𝐶𝐷 𝑜𝑝𝑡𝑖𝑚𝑎𝑙 | > 1, then a new iteration will start, and the actions specified in step 1 will be performed again.

The result of the described optimization procedure is the geospatial data cube with the set of layers that achieves the highest value of the SITS among all other considered sets. This geospatial data cube will be used for further land cover classification.

Experiment

The experiment was conducted to test the effectiveness of the developed technique. It consisted of carrying out a classification of the selected study area. For this purpose, an initial geospatial data cube with an excessive number of layers was formed. Then, the optimization technique was applied, resulting in an optimized cube. Finally, two classifications were obtained -before and after the optimization.

Study area

Since 2007, the Ukrainian network of test sites has provided validation and calibration of various remote sensing techniques and satellite-based products, including land cover classification [32]. The proposed technique was tested at the site within the Shatsk National Natural Park (SNNP). It is situated in the northwest of Ukraine, within Volyn' oblast, between 51º 28'25"N and 23º 49'29"E. The SNNP encompasses highly heterogeneous natural landscapes, like forests, peat bogs, transitional mires, meadows, and lakes.

The site comprises more than 100 georeferenced sample plots and gives comprehensive ground truth information about the representative landscapes of the West Polissia region (Figure 3).

Training sample

Six broad land cover classes that characterize the study area were defined: artificial surfaces, tree-covered areas, grassland, agricultural areas, water bodies, and wetlands.

The given classes varied considerably both in spatial extent and heterogeneity. The smallest class of artificial surfaces included diverse features of built-up areas and transport units. While the biggest ones, like tree-covered areas and wetlands, included various sub-types that could still be quite homogenous due to the big extent. Water bodies represent the most homogeneous class. Therefore, the number of training pixels of each class also varied disproportionally. The overall number of all training pixels accounted for 6474. Table 1 shows labels, descriptions, and training pixel amounts for the land cover classes assigned for the experiment.

Initial geospatial data cube

The experimental classification focuses on wetlands, the prevailing landscape of the test site and one of the most important for conservation within the SNNP. Evident differences in the seasonal development of wetlands and other vegetative land cover classes help distinguish them and require the application of multitemporal data [12].

The primary data source for forming the geospatial data cube was Sentinel-2 satellite imagery [33]. The images were selected for 4 dates (04.07.2018, 05.12.2018, 06.01.2018 and 10.14.2018) with minimal or no cloudiness.

Each Sentinel-2 image contains 13 spectral bands. At the preprocessing stage, atmospheric correction was performed for each image to eliminate the influence of the atmosphere and calculate the pixel values corresponding to the surface reflectance (bottom of atmosphere). During this procedure, 3 bands (B1, B9 and B10) that consider the effects of aerosols and water vapour on reflectance were removed. The spectral bands of the Sentinel-2 image have different spatial resolutions, namely 10, 20 and 60 meters.

Next, a complete set of normalized difference indices was calculated for each image. The following combnatorial formula is used to calculate such an index:

𝐼𝑛𝑑𝑒𝑥 = 𝑏 𝑖 − 𝑏 𝑗 𝑏 𝑖 + 𝑏 𝑗 , 𝑖 ≠ 𝑗,

where 𝑏 𝑖 is the spectral band i. This set could be presented in the following form:

𝑁𝐷𝐼 = { 𝑏 𝑖 − 𝑏 𝑗 𝑏 𝑖 + 𝑏 𝑗 |𝑖, 𝑗 ∈ {1,2 … ,10}, 𝑖 ≠ 𝑗}.

The cardinality of this set (i.e. number of normalized difference indices of one image) is calculated by the formula below:

|𝑁𝐷𝐼| = 𝐶 𝑛 𝑚 = 𝑛! 𝑚! (𝑛 − 𝑚)! = 10! 2! (10 − 2)! = 8! * 9 * 10 1 * 2 * 8! = 90 2 = 45,

where variable n corresponds to the number of image bands, and m is the number of arguments in the index calculation formula.

As a result, 45 spectral indices were obtained for each image. Another component of the input cube was the geomorphological data obtained from the ALOS PALSAR DEM [34]. In particular, this data contains the height above sea level and the slope. Their spatial resolution is 12.5 m.

All the above-described data must be spatially regularized to form the input geospatial data cube. It involves bringing all layers to the same spatial resolution, map projection, and size. With this in mind, all layers were scaled to a spatial resolution of 10 m, transformed to Universal Transverse Mercator Projection, Zone 34N (EPSG:32634) and resized so that all layers lie within the study area.

So, the input geospatial data cube contained 222 raster layers, namely 40 spectral bands of 4 different time Sentinel-2 satellite images, 180 corresponding spectral indices, and 2 raster layers of geomorphological parameters. The layers of the input cube are described in Table 2. Hence, signatures of the initial training sample were assigned in each layer of the initial cube. To assess the separability of the training sample, the SITS was calculated using Formula 1:

𝑆𝐼𝑇𝑆 𝐼𝑛𝑖𝑡𝑖𝑎𝑙 = ∑ 𝑥 𝑖𝑖 𝐾 𝑖=1 𝑁 = 6469 6474 ≈ 0.9992.

As seen above, this training sample had 5 misclassified signatures.

Optimized geospatial data cube

After applying the developed optimization technique, the geospatial data cube size was reduced from 222 to 40 layers. The selected layers are listed below in Table 3. Along with the cube optimization, the training sample signatures were reassigned according to the selected cube layers. The separability of that training sample was assessed by the SITS calculated below:

𝑆𝐼𝑇𝑆 𝑂𝑝𝑡𝑖𝑚𝑎𝑙 = ∑ 𝑥 𝑖𝑖 𝐾 𝑖=1 𝑁 = 6474 6474 = 1.

This training sample had no misclassified signatures, so the optimized training sample was entirely separated.

Land cover classifications

In order to test the effectiveness of the developed technique, the classifications using the initial cube and the optimized one were compared.

Firstly, the classification was obtained using the initial cube and the appropriate training sample. Then, the classification was obtained using the optimized cube and the appropriate training sample. These classifications are depicted in Figure 4. The classifications above were obtained using Mahalanobis distance [35] as a supervised classification method. Exactly for this method, the separability assessment of the training sample was carried out for both initial and final geospatial data cubes.

Accuracy assessment

Classification accuracy assessment involved independent verification of initial and final land cover maps using proportionate stratified random samplings. This sampling technique produces sample set sizes directly related to the size of the classes and is widely used in assessing the classification accuracy of classes disproportionate in their extent. As the required sample size for a class, 0.01% of the total classified pixels of this class were analyzed. Thus, test sample sets were equal to 300 pixels for each land cover map. Satellite images (QuickBird) of high spatial resolution, available for 2018 in the Google Earth Pro app, were used for verification as reference data.

Table 4 shows the confusion matrix of the initial land cover map. In addition to the two dimensions ("Reference" and "Prediction"), this matrix shows metrics such as producer accuracy (PA) and user accuracy (UA) for each class [31]. The accuracy of the obtained land cover classification was assessed by indicators of overall accuracy and the kappa coefficient [31]. The overall accuracy value was calculated by the following formula: Table 5 shows the confusion matrix of the final land cover map.

The same indicators were selected for the final land cover classification as for the initial one. Thus, the value of overall accuracy value is stated below:

Discussion

The aims of the developed optimization technique are stated as layers reduction of the geospatial data cube and enhancement of the classification. Thus, the experiment result should be considered in terms of these two aspects. Firstly, the size of the optimized cube was 40 layers, whereas the initial one contained 222 layers. Therefore, the number of layers was reduced by 5.55 times.

Secondly, classification enhancement was evidenced by increasing indicators such as overall accuracy and the kappa coefficient. Namely, the overall accuracy increased by 0.04 from 0.9 to 0.94, and the kappa coefficient increased by 0.06 from 0.86 to 0.92. Since the classification focused on the wetlands, the accuracy of this class should be considered individually. Thereby, both user and producer accuracy of wetlands class were significantly increased, viz. by 0.11 from 0.71 to 0.82 and by 0.4 from 0.96 to 1, respectively.

Conclusion

This article presents an optimization technique to reduce geospatial data cube size and enhance land cover classification. The technique is based on the separability of the training sample, which is defined in each layer of the geospatial data cube. To assess the separability, the appropriate index (i.e. SITS) was developed and used as an object function in the technique frame. The algorithm of the optimization technique implies stepwise band discarding to define the optimal set of the geospatial data cube layers. Such a set has the highest value of SITS among other options.

The conducted experiment implied techniques application to the land cover classification of the highly heterogeneous natural landscapes in the case of the Shatsky National Natural Park. This classification covered six land cover classes where wetlands are prevailing. The technique's effectiveness was approved by geospatial data cube reduction and classification accuracy enhancement, evidenced by the increase in such indicators as overall accuracy and kappa coefficient.

Further research should be aimed at technique application in other study areas and thematic tasks. Also, the separability assessment of the training sample could be extended by additional criteria. For example, the kappa coefficient could substitute overall accuracy as the basis of the developed separability index.