Digital Mapping of Land Cover Changes Using the Fusion of SAR and MSI Satellite Data

Abstract

The aim of this article is to choose the most appropriate method for identifying and managing land cover changes over time. These processes intensify due to human activities such as agriculture, urbanisation and deforestation. The study is based in the remote sensing field. The authors used four different methods of satellite image segmentation with different data: Synthetic Aperture Radar (SAR) Sentinel-1 data, Multispectral Imagery (MSI) Sentinel-2 images and a fusion of these data. The images were preprocessed under segmentation by special algorithms and the European Space Agency Sentinel Application Platform (ESA SNAP) toolbox. The analysis was performed in the western part of Lithuania, which is characterised by diverse land use. The techniques applied during the study were: the coherence of two SAR images; the method when SAR and MSI images are segmented separately and the results of segmentation are fused; the method when SAR and MSI data are fused before land cover segmentation; and an upgraded method of SAR and MSI data fusion by adding additional formulas and index images. The 2018 and 2019 results obtained for SAR image segmentation differ from the MSI segmentation results. Urban areas are poorly identified because of the similarity of spectre signatures, where urban areas overlap with classes such as nonvegetation and/or sandy territories. Therefore, it is necessary to include the field surveys in the calculations in order to improve the reliability and accuracy of the results. The authors are of the opinion that the calculation of the additional indexes may help to enhance the visibility of vegetation and urban area classes. These indexes, calculated based on two or more different bands of multispectral images, would help to improve the accuracy of the segmentation results.

1. Introduction

Land and soil are limited nonrenewable resources that disappear over time. This process is particularly accelerated by the intensity of various human activities that have transformed and fragmented ecosystems of the land, as well as the functionality and performance of shared land resources. The anthropogenic influences on climate are the emission of greenhouse gases and changes in land use, such as urbanization and agriculture. Of course, soil erosion, soil organic matter decline, soil contamination and compaction and surface temperature also contribute to this decline, but human activities, such as land development, intensive land use and land abandonment, only accelerate these processes [1,2,3]. The scientist excreted three highest-emitting regions (Latin America, Southeast Asia and sub-Saharan Africa) that dominated global emissions growth from 1961 to 2017, driven by the rapid and extensive growth of agricultural production and related land-use change. From the point of view of environmental protection, landscape morphology and functioning, economically and socially, land cover is considered to be a very important object of research [4]. Firstly, stabilizing the local surface temperature level and achieving a balance between anthropogenic emissions by sources and removals by sinks of greenhouse gases is very important, now and in the future globally [5].

The term ‘land cover’ itself is similarly described by various authors and international organisations: it is the physical (bio) cover of the Earth’s surface, including meadows, forests, artificial surfaces, waters, etc. [6,7,8]. Land cover is constantly changing and is affected by economic, political, cultural and other societal factors [9,10,11,12,13,14,15,16,17,18,19,20]. Long-term land cover observations can be used to assess the factors that affect the landscape to assess their extent, to predict trends in landscape change, to identify patterns of development, to prevent destructive processes [21] and to assess natural or man-made damage to the landscape. Monitoring land cover change is a relevant topic but, at the same time, a challenge to conduct with as few resources as possible and to obtain the most efficient process. From a historical–temporal point of view, it is easy to see that changes in land cover have been largely influenced by human activity. It is often important to know how land cover changes over time, and to assess what changes we can expect in the near future, as well as the impact this will have on people’s lives. Today, the observation of changes at the Earth’s surface is becoming easier with the advent of remote sensing and, more specifically, with various satellite images. Data from various remote sensing applications, such as panchromatic, multispectral instruments (MSI) and hyperspectral, synthetic aperture radar (SAR) images, covering different parts of the electromagnetic spectrum, are obtained from different artificial Earth observation satellites. These data can be processed and used for a variety of tasks but, in many cases, the use of only one type of image may be insufficient for land cover observations. The fusion of different types of data becomes an excellent solution for observing and analysing the surface of the Earth to obtain more information about an object. Information collected from different image sensors and properly processed by synthesising satellite images allows data gaps to be filled and more accurate results to be obtained [22]. Various studies on the fusion of different data can be found in the scientific literature [23,24,25,26,27,28,29,30,31,32,33]. Such works began in the 1980s, but some of the processes have not been automated.

Some authors used optical images with different resolution (10–30 m). Data fusion results were images with 31 spectral bands, and a 30 m spatial resolution was performed and used to identify soil information. The modelling process produced an increase of approximately 10% in prediction accuracy compared to that obtained from single-sensor images [26]. The data fusion technics apply when scientists have optic image limitations because of cloudy times. The authors identify the process of combining data from multiple sources to produce more accurate, consistent and concise information than that provided by any individual data source [27]. The literature review by Jayme Garcia and Arnal Barbedo highlights the fusion techniques and accuracy. The Convolutional Neural Networks (CNN), Backpropagation Neural Network (BPNN), Random Forest (RF), Partial Least Squares Regression (PLSR), Dempster–Shafer and the Kalman filter and least squares support vectors machine gave the best mean precision of the result [28]. The authors explain that in the case of agriculture, image data fusion has been particularly effective at the orbital level, both for artificially increasing the spatial resolution of sources with revisit frequencies and for compensating cloud cover using the information present in SAR images. In these cases, the improvement can exceed 50%. There is a study by N. Kussul and coauthors [29] in which vegetation was classified using fusion of satellite images EO-1 and RADARSAT-2. In their study, they distinguished four classes: corn, soybeans, sunflower, and winter/spring wheat. They show that the result obtained with the data fusion is 2.5% more accurate compared to the result where only SAR images were used and 6.6% more accurate than with optical images only. Nicola Clerici and coauthors [30] used Sentinel-1 and Sentinel-2 in their study. The aim of the study was to investigate the precision of land cover classification using different algorithms, classifying SAR and MSI images separately and their fusion. Their result shows that regardless of the algorithm used, image fusion gives a more accurate result. They obtained a 58% more accurate result using the random forest algorithm than using SAR images alone for land cover classification, and it was 16% more accurate than using only optical images. The same Sentinel data were used by other authors, such as Julien Denize et al. [31], who show that SAR and MSI data fusion is 21% more accurate than SAR-only images and 7% more accurate than MSI-only images. Audrey Mercier and coauthors [32] also performed better in image fusion, which is 20% more accurate than using SAR data alone and 5% more accurate than the results obtained using optical images alone. A very similar result was obtained by Armugha Khan, Himanshu Govil, Gaurav Kumar and Rucha Dave [33], who distinguished five land cover classes using S-1 and S-2 satellite imagery—waters, urban areas, ravines, barren land and crops. Their result is 25% more accurate with image fusion than SAR data alone and 5% more accurate with optical images alone.

There is no consensus on the most appropriate method for monitoring and assessing land cover change, so this study focuses on these processes and their improvement. The different MSI and SAR data fusion techniques have been used for the segmentation of land cover in this study. The purpose of this article is to choose the most appropriate data fusion technique for identifying and managing land cover changes over time.

2. Study Area

An area of 2000 km² in western Lithuania was selected for the study (Figure 1). This territory was chosen due to the diversity of land cover, which includes one of Lithuania’s largest cities (Klaipėda), less urbanised areas (Palanga, Kretinga, Kretingalė, Gargždai, etc.), main roads, side roads, wooded areas, cultivated fields, meadows, sandy areas, coastal dunes, open waters (Baltic Sea), lagoons (Curonian Lagoon), closed waters (lakes and ponds), rivers and swamps.

Six classes of land cover (hydrography, forests, vegetation areas, sand dunes, built-up areas and nonvegetation) were selected for the study because they were recommended by other authors for land cover/land use classification [7,34,35,36,37,38,39,40,41,42,43].

In the study area, the hydrography class includes the sea, lagoon, water bodies and rivers. Areas belonging to the forest class are covered with forests (including young stands). The class of vegetation areas are those covered with green vegetation. Sand dunes include the coastal zone, quarries and other similar areas. The class of built-up areas includes settlements, artificial surfaces and industrial zones, while no vegetation areas include felled areas, bare soil and areas where vegetation does not grow. The methods for preparing fusion of the data from the different satellite images described in this article were used to monitor and visualise changes in the Earth’s surface.

3. Materials and Methods

3.1. Data Acquisition and Pre-Processing

The Sentinel-1 SAR and Sentinel-2 MSI satellite images distributed by the European Space Agency (ESA) are available for free download from the Copernicus Open Access Hub website. These images were selected based on their resolution, spectral bands, recording frequency and coverage area [44]. The SAR and MSI data provided by the ESA have different characteristics, both in terms of their physical properties and the principle of obtaining the images themselves; therefore, pre-processing of the data is required for data preparation and fusion (Figure 2).

Section 3.1.1 and Section 3.1.2 provide detailed information on the pre-processing of data for the study.

3.1.1. Pre-Processing of Sentinel-1 SAR Satellite Data

Sentinel-1 is the first satellite of the Copernicus Global Monitoring Program for Environmental Monitoring. The satellite was developed to observe the subsystems of the Earth, the atmosphere, the oceans and the continents [45,46]. For this study, 13 S1A images from the Sentinel-1 satellite were chosen. It was Interferometric Wide swath (IW), Single Look Complex (SLC) products with dual VV+VH polarisation. The processing of images was performed using Graph Builder from the ESA SNAP toolbox.

In the first pre-processing step, the satellite position and velocity information was refined (Apply Orbit File). High-accuracy image positioning is important in the process of fusion of different data. The exact position of the satellite is not known when the SAR image is taken in real time, so only abstract information is captured in the product metadata. The ESA SNAP toolbox uses the Apply Orbit File to update the satellite position and speed information. This uses orbital files that are available a few days or weeks after the image is captured. In ESA SNAP, these files are automatically downloaded [45]. The data are then calibrated and the satellite images acquire the true values of the reflected rays (sigma, beta and gamma); the bursts and subswaths in the original image were merged into a single image, and the image was converted to form uniformly sized grids. In addition, a Speckle Filter was applied, which reduced speckles in the image, making it difficult to interpret. For this study, a Lee Sigma filter was chosen, which uses a standard deviation to reduce speckles in the image. Pixel contrast is important in identifying different ground cover objects, so the polarisation bands of images need to be transformed to a logarithmic scale to highlight them. This feature shifts bright values towards the middle, and dark values are stretched over a wider range of colours, resulting in greater contrast between some of the different land objects. Original SAR images are often distorted and their position appears mirrored due to the orbit of the satellite; therefore, after the initial processing of the SAR images, a geographic orientation was performed as the last step to ensure that the images are geographically correct. The SRTM 1 sec HGT terrain model, 10 m pixel size, and WGS84/UTM zone 34N coordinate system were used for this operation. An example of an image processed and prepared for analysis is shown in Figure 3.

3.1.2. Pre-Processing of Sentinel-2 MSI Satellite Data

Mission Sentinel-2 is dedicated to monitoring variability in surface conditions and changes in vegetation during the growing season. These satellites systematically capture 13 different multispectral channels. The study used MSI level-2A (L2A)-type satellite images, which were selected according to an important criterion-cloud percentage. Images in June, July and August of 2018 and 2019, with cloud cover levels below 15%, were used. The L2A MSI images provided by ESA were already partially processed, i.e., the data were corrected due to atmospheric and terrain errors. Therefore, the downloaded data were immediately suitable for data fusion and did not require any additional processing by software. However, to avoid data redundancy and speed up further processing, image cropping was performed, which is a subset of data selection. Solar position, azimuth, or quality indicators are not required for data fusion, so only MSI image bands were selected, eliminating aerosol and cloud-sensitive tapes (B1, B9 and B10 tapes). Subset tools of the ESA SNAP toolbox software were used. The final data product consisted of bands B2, B3, B4, B5, B6, B7, B8, B8a, B11 and B12. An example of an image prepared for analysis is shown in Figure 4.

3.2. Segmentation and Land Cover Classification

In the study, we used four methods of image processing and fusion:

• The coherence of two SAR images;
• The method when SAR and MSI images are segmented separately and the results of segmentation are fused;
• The method when SAR and MSI data are fused before land cover segmentation;
• The upgraded method of SAR and MSI data fusion by adding additional formulas and index images.

3.2.1. The Coherence of Two SAR Images

In the study, we estimated land cover classification using only SAR images and MSI images as a support to fill in the gaps. The methodology of Magdalena Fitrzyk (2019) was used when the change in the land cover was identified using the coherence of electromagnetic waves, and the result is presented in RGB composition [47].

Two images of SAR with equal territorial coverage in early June 2018 and early 2019 were selected to test this method. The following SNAP tools were used for pre-processing both images: Apply Orbit File, back geocoding, enhanced spectral diversity, coherence estimation, debursting and terrain correction. The back geocoding operator co-registers two S-1 SLC split products (master and slave) of the same subswath, using the orbits of the two products and a Digital Elevation Model (DEM). This created one image, for which we could estimate coherence. After coherence estimation for the same raw SAR images, the second pre-processing has to be carried out. The scheme of the second pre-processing is given in Section 3.1.1. These two preprocessing stages result in three images for the fusion: one coherence SAR image and two preprocessed SAR images (Sigma0_vv_2018 and Sigma0_vv_2019). For the fusion, the SNAP Stack tool was used. After fusion, the images were converted to the logarithmic scale, creating two new bands: average (averange_sigma) and difference (diff_sigma).

$$diff\_sigma=Sigma0\_vv\_2018-Sigma0\_vv\_2019$$

(1)

$$average\_sigma=\frac{Sigma0\_vv\_2018+Sigma0\_vv\_2019}{2}$$

(2)

All three bands (coherence, diff_sigma and average_sigma) were then used for the RGB composite and further analysed. The coherence band was used to open this RGB image correctly for the Red band; for the Green band: average_sigma; and for the Blue band: diff_sigma. The processing schema of the coherence method is presented in Figure 5. The results are presented in Section 4.1.

3.2.2. Separate SAR and MSI Segmentation Method

Land cover segmentation was performed separately for the summer SAR and MSI images of 2018, using the supervised segmentation method and the regression-type Random Forest (RF) algorithm. The pixel and coordinate results were then combined into an overall segmentation result. The data pre-processing is described in Section 3.1.1 and Section 3.1.2. Land cover classification was performed separately for SAR and MSI images, according to Section 3.3. The processing scheme of the separate SAR and MSI segmentation method is shown in Figure 6.

The analysis of the obtained results is presented in Section 4.1.

3.2.3. SAR and MSI Fusion Technique

In the study, the SAR and MSI images were combined into a common image after pre-processing. The merging of the images was performed using the ESA SNAP software, using the Collocate tool. The SAR image was the Master image. This means that all other image views are resampled according to their geographical position and grid size. The SAR image was selected as the primary one because the images are extremely accurate after applying the orbit file during image preprocessing. Resampling by bilinear interpolation was used for data aggregation. Bilinear interpolation (bilinear filtering) is a technique for calculating values of a grid location based on nearby grid cells. The key difference is that it uses the average digital number (DN) of the cell centers closest to the four pixels (Figure 7) [48].

The four nearest-neighbour image pixel centre coordinates $f\left(n_{10}, n_{20}\right)$, $f\left(n_{11}, n_{21}\right)$, $f\left(n_{12}, n_{22}\right)$ and $f\left(n_{13}, n_{23}\right)$ geometrically transformed the image g(n1, n2), as computed by [49]:

$$g\left(n_1, n_2\right)=A_0+A_1{n}_1+A_2{n}_2+A_3{n}_1{n}_2$$

(3)

where $A_0$, $A_1$ and $A_2$ (the bilinear weights) were found by solving the matrix:

$$\left[\begin{array}{c}{A}_0\\ A_1\\ A_2\\ A_3\end{array}\right]={\left[\begin{array}{cccc}1& n_{10}& n_{20}& n_{10}{n}_{20}\\ 1& n_{11}& n_{21}& n_{11}{n}_{21}\\ 1& n_{12}& n_{22}& n_{12}{n}_{22}\\ 1& n_{13}& n_{23}& n_{13}{n}_{23}\end{array}\right]}^{-1}\left[\begin{array}{c}f\left(n_{10},n_{20}\right)\\ f\left(n_{11},n_{21}\right)\\ f\left(n_{12},n_{22}\right)\\ f\left(n_{13},n_{23}\right)\end{array}\right]$$

(4)

where $g\left(n_1, n_2\right)$ is defined as a linear combination of the grey levels of its four nearest neighbours. The linear combination defined by Equation (4) is, in fact, the value assigned to $g\left(n_1, n_2\right)$ when the best (least squares) planar fit is made to these four neighbours. This process of optimal averaging produces a visually smoother result.

Regardless of the interpolation approach that is used, it is possible that the mapping coordinates $A_1\left(n_1, n_2\right),$ $A_2\left(n_1, n_2\right)$ do not fall within the pixel ranges.

We obtained one GEOTIFF format image after SAR and MSI data fusion and the data used in the segmentation process are described in Section 3.3. The processing schema of SAR and MSI data fusion are given in Figure 8.

The results are given in Section 4.1.

3.2.4. Image Fusion Technique with Additional Indexes

In this study, we decided to update the fusion technique of the images (Section 3.2.3) with derivative images, i.e., indexes. The accuracy of the identification of all land cover classes can be improved by using derivative products of the images. Different land cover objects reflect multispectral waves differently so, in order to obtain more accurate land cover segmentation results, indices that are typically calculated from two or more different bands of the multispectral image can also be used.

After the segmentation process (Section 3.3) of the least-identified and overlapping urbanised, nonvegetation and sandy areas, three additional indices were chosen to highlight the vegetation class: Normalised Vegetation Difference Index (NDVI), Sentinel-2 Red Edge Position Index Sentinel-2 Red-Edge Position Index (S2REP), Green Normalised Difference Vegetation Index (GNDVI) and one to highlight the urban class—Normalised Difference Built-Up Index (NDBI).

NDVI is the most commonly used index in the remote sensing analysis of vegetation, which is also widely used in land cover mapping [30].

Clerici et al., (2017) argued that it may also be useful to include information from the red-edge bands of the electromagnetic spectrum for land use classification. The Sentinel-2 S2REP index is based on linear interpolation using Sentinel-2 bands B5 and B6, both of which are for red edge and suitable vegetation identification [30].

GNDVI was calculated using data from the visible green spectrum instead of the visible red spectrum (as in the NDVI index) and the near-infrared spectrum. It is useful for measuring the intensity of photosynthesis and monitoring plant stress; due to the peculiarities of the spectral bands used in the indicator, it is more sensitive compared with NDVI. This index not only helps to identify vegetation but also highlights differences between hydrographic and vegetation classes (hydrographic objects acquire negative values) [50].

All three of these indices are designed to better identify vegetation in satellite imagery. Analysing the intermediate results of land cover classification, it was observed that urban areas are difficult to distinguish and mix with other classes, so it was decided to use an index that would highlight urban areas. One of these is the Normalised Difference Built-Up Index (NDBI). This index highlights urbanised areas with higher shortwave reflections in the infrared spectrum [51,52].

In addition to the above indices for multispectral images, calculations were also performed for the processed SAR images. To reduce band imbalance, the ratio between the two polarisations is calculated using the VV and VH bands, according to the formula below:

$$Sigma\_ratio\_db=Sigma0\_VV\_db-Sigma0\_VH\_db$$

(5)

These indices were calculated for each image separately, regardless of how many of them are intended to be used in the fusion process. All of the calculated indices in the process were used as additional bands, so the final multispectral photo product consisted of B2, B3, B4, B5, B6, B7, B8, B8a, B11, B12, NDVI, S2REP, GNDVI and NDBI bands, and the SAR product- $Sigma0\_VV\_db$, $Sigma0\_VH\_db$ and $Sigma\_ratio\_db$ bands. The results and statistical analysis are presented in Section 4.1 and Section 4.2.

3.3. Image Segmentation and Identification of Land Cover Change

The accuracy of a map of land cover classes depends not only on the data used to compile it but also on the method used for segmentation and the algorithm used. Three methods are used to segment digital images: supervised, unsupervised and object-based segmentation. In this study, the supervised segmentation method was chosen. When classifying land cover, the set of samples are used to define the classes. The operator collects sample areas by studying the maps of the selected territory, as well as aerial photographs and orthophotographic maps, checking selected sites using field surveys. The goal is to identify a set of areas that accurately reflect the spectral values in each land cover class. The classification algorithm or classifier is trained to select classes using a group of spectral features and a training sample. When signatures for each class are identified, each pixel in the image is compared to those signatures and marked as belonging to the class that it most numerically matches [17,53].

Various authors focusing on land cover classification studies using SAR and MSI images mainly use two main classification algorithms: Support Vector Machine (SVM) and Random Forest (RF) [29,30,31,32,54,55,56]. RF is one of the most popular and effective supervised training algorithms. The basic idea of RF is to form an accurate classifier by combining solutions from multiple binary solution trees, grown using different subsets of data from the original data set, and randomly selected subsets of features from the feature set. Such a set of binary trees is characterised by resilience to learning and, as the number of trees grows, the convergence of the common error to a stable value [57]. Each decision tree is trained by randomly separating a portion of the data from the training sample, and the remaining data is used for testing. As the number of decision trees increases, the error for the test portion of the data decreases. The significance of the characteristics is measured using the remaining data. The selected characteristics depend on the number of data included in the remaining data, so different amounts of data are taken for the selection of the most important characteristics. By eliminating the features in this way, the less important features are removed, as long as the most important spectral features remain. Analysing the work by Kussul, Clerici, Valbuena Calderón and Posada, Haas and Ban, Gerrells, Sun, Denize and Mercier, which compares the results obtained using the SVM and RF algorithms, shows that the RF algorithm often yields more accurate ground cover classification results; this algorithm was chosen for this study.

After land cover segmentation, the raster was converted to GeoTIFF-BigTIFF format, so that the image does not lose its geolocation information and is suitable for further analysis with GIS applications. Often, there are single misclassified grids in the segmentated raster that look like noise. In this case, a reclassification of the grids should be performed: the value of another class with which it shares four boundaries was transferred to the individual grids.

Rasters processed and prepared for two different periods to identify land cover changes were converted to vector data. This step was performed to remove cells that did not undergo a change in land cover, leaving only those areas where the change was recorded.

The change in land cover classes was estimated by intersect (Union tool in GIS app); this vector classifies land cover layers for different periods. A new data layer was obtained, from which it is possible to see which two classes were intersected.

The first filtering was performed according to the overlapping grid identifier, which means that those grids that have the same land cover class identifier in the images of different periods were removed; in other words, the grids in which the land cover change did not occur. Grids in which there was a change between vegetation and uncovered vegetation classes were also removed because it is a seasonal change in land use rather than a change in land cover. When grids with changes are selected, they are merged according to the land cover class identifier, but without creating multipart geometry objects, which means that objects of the same class that have a common boundary will be merged.

After the merging, a second filtration was performed, during which areas smaller than 5 ha were removed. This number was chosen based on the scale of the map being created, which in this case is M 1:50,000. On a map of this scale, an area of 5 ha is a rectangle of 4 × 5 mm², so it is visible on both the digital and printed maps. At the same time, small changes are distinguished from significant, larger changes. After filtration, changes in land cover over the analysed period were obtained.

4. Results

To review the process of all the research presented in the article, a general methodological scheme for the identification of land cover changes has been developed (Figure 9).

The results section is divided into two parts: the results of data segmentation and comparisons, in visual terms, and the results of land use change and accuracy assessment in 2018 and 2019 in statistical terms.

4.1. Data Segmentation and Comparison Results

After the pre-processing and coherence assessment (Section 3.2.1) of two SAR images (2018, 2019), three data bands were obtained: SAR images after pre-processing in 2018; SAR images after pre-processing in 2019; and the coherence band. According to the methodology described in Section 3.3, the coherence of the data for 2018 and 2019, and their aggregation to assess the changes, was performed. The results are shown in Figure 10.

The results obtained from the assessment of land cover changes in 2018 and 2019 were visually analysed and compared with RGB satellite images from the same period. Figure 10 shows that the forests and urban areas were sufficiently clear and precise and that the riverbed did not change throughout the year. However, when analysing arable land, meadows, pastures, and other vegetation areas, annual differences have been identified that occupy a considerable area. Comparing these areas with the MSI RGB image shows that their actual value has not changed, meaning that if that area was used for agriculture, it was used for the same purpose a year later, but there was a change between the stage of vegetation growth. In 2018, it was observed that the vegetation season started later than in 2019, or other crops were grown. Coherently prepared SAR data are suitable for visual analysis and comparison but are completely unsuitable for mathematical analysis. It is not possible to calculate the changes between classes. This is because it is possible to see which grids belong to which class (at a small scale), but analysing them separately shows that they have neither a uniform colour value specific to that class nor an identifier to distinguish object classes. Consequently, this method is suitable only for initial inspection of the site and identification of the most altered site.

SAR and MSI preprocessed data were segmented separately in the next study method. The results of segmentation are shown in Figure 11.

Analyzing the images in Figure 11, it was found that the results obtained from the MSI images are much more in line with the real situation in the land cover than the images classified by SAR. In the SAR-segmented photograph, forest areas are highlighted, but all other object classes show significant discrepancies with the actual situation. Meanwhile, the analysis of MSI-classified images shows a fairly well-defined urban area but, in places, they overlap with areas of similar colour spectrum, i.e., sandy areas, dry ploughed land, clear-cuts, and so on. MSI images also show a very good distinction between smaller and larger wooded areas, larger water bodies, and rivers. Smaller water bodies remain unidentified (as would be in SAR images). This methodology highlights the problems of SAR and MSI images quite well, i.e., the ability to identify different ground cover classes using one or another data source. Classifying SAR and MSI images separately also results in too large a difference in the results to be later fused into a single map of the land cover classes.

In the next study, the prepared SAR and MSI images were fused into one image before the segmentation process (Section 3.2.3). The results of June 2018 and 2019 were obtained after calculating the fused image (Figure 12).

Compared with the previously obtained results (Figure 11), this result is better, but a closer analysis of the obtained images shows that the urban areas were relatively poorly separated. They often overlap with the no vegetated land cover class or with sandy areas. To improve accuracy, it was decided to include the indices described in Section 3.2.4 in the segmentation process to highlight the vegetation and urban classes. The segmentation view, including the indices NDVI, S2REP, GNDVI, NDBI and $Sigma\_ratio\_db$, is shown in Figure 13.

The statistical results from this method are presented in Section 4.2 and the recorded changes in land cover are presented for the period 2018–2019.

4.2. Results of Change Assessment

The survey identified 511 changes in land cover class in the selected area during the summer periods of 2018 and 2019. The statistics are given in Figure 14. The abscissa shows the change in the area of land cover class, from one class to another. Nonvegetated areas turned into forests (0–1), sand dunes (0–2) and urbanised areas (0–3); forests turned into uncovered areas (1–0) and vegetation areas (1–4); urban areas became nonvegetated areas (3–0), forests (3–1), vegetation areas (3–4) and even hydrography (3–5); and vegetation areas turned into forests (4–1), urban areas (4–3) and the hydrography class (4–5).

In the period 2018–2019, most land cover changes were recorded in two class pairs: 0–1 and 4–1. Nonvegetated and vegetation areas of the land cover changed to forests. The least changes were recorded in the sand and hydrography classes. To estimate the actual change in all class pairs, 20 random objects (or all if fewer changes were identified) were selected and the result was checked against the 2018 and 2019 RGB images. Figure 15 is a diagram summarising the results.

This shows that as many as 5 of the 12 pairs of land cover classes listed have a 100% negative result, meaning that no actual change has occurred at those locations and classification inaccuracies have been identified. This analysis proves, once again, that despite the use of SAR data in addition to MSI data, the identification and segregation of vegetation areas, urban areas and nonvegetated areas still remains difficult. This was also found in the work of other authors [30,58].

5. Discussion and Conclusions

To summarise the findings of this study:

• After preprocessing and coherence extraction of the SAR and the two-period images, the result is only relevant for the initial inspection of area changes and the detection of the most altered areas.
• The results from using the second method show that the individual segmentation of SAR and MSI images differs drastically, and the field studies and accuracy estimation required the evaluation of reliability.
• The result of segmentation of the combined SAR and MSI data showed that the classes of urban areas, nonvegetated areas and sandy areas are poorly separated due to a similar spectral signature. Therefore, it was decided to include NDVI, S2REP and GNDVI indices to improve accuracy and highlight the class of vegetation areas. NDBI was included to highlight the class of urban areas.
• Additional indices improved the result of segmentation, but there are still errors in identifying urban areas.
• Changes that were falsely identified during the qualitative accuracy check of the identified changes (92.08% of all changes checked) were False Positive results and no False Negative results were observed in the analysis of the images. Although changes are incorrectly identified in some identified cases, visual inspection (especially when potential locations for potential inaccuracies are known) and manual correction would still use less time than not automating all the process.

In our study, the method that provides the most accurate results for land cover segmentation confirms Clerici’s [30] finding that the inclusion of additional indices enriches the image and yields better segmentation results. As many as three indices were used in this work to highlight the vegetation class (NDVI, S2REP and GNDVI), but questions arise as to whether this is not a surplus. In order to ascertain which vegetation indices are most beneficial for the segmentation result, a broader analysis of the literature and, if necessary, research is provided. Future work will aim to refine the method so that it is optimal in terms of time and quality compared with similar manual work.

The applied image merging techniques also confirm the conclusions of other authors [26,27,28,29,30,31,32,33] that the merging of two different sensors enriches the image information. Moreover, SAR images may compensate for the low cloud content of MSI images, especially if multiple SAR images are used. However, it is interesting that in all the analysed authors, the percentage accuracy of the fused and separately segmented images varies quite strongly in some cases. This may be influenced by the satellite imagery used, the different ground cover classes used and the different classification algorithms used. It is also very important to note that in the articles analysed, the research is carried out on different areas that are in different latitudes, which may also be the reason why the results obtained by different researchers are so different.

Our plan for future work is to carry out studies and projections of the impact of land cover/land use classes on future climate change and expand the study spatially and temporally.