A Collaborative Change Detection Approach on Multi-Sensor Spatial Imagery for Desert Wetland Monitoring after a Flash Flood in Southern Morocco

Abstract

This study aims to present a technique that combines multi-sensor spatial data to monitor wetland areas after a flash-flood event in a Saharan arid region. To extract the most efficient information, seven satellite images (radar and optical) taken before and after the event were used. To achieve the objectives, this study used Sentinel-1 data to discriminate water body and soil roughness, and optical data to monitor the soil moisture after the event. The proposed method combines two approaches: one based on spectral processing, and the other based on categorical processing. The first step was to extract four spectral indices and utilize change vector analysis on multispectral diachronic images from three MSI Sentinel-2 images and two Landsat-8 OLI images acquired before and after the event. The second step was performed using pattern classification techniques, namely, linear classifiers based on support vector machines (SVM) with Gaussian kernels. The results of these two approaches were fused to generate a collaborative wetland change map. The application of co-registration and supervised classification based on textural and intensity information from Radar Sentinel-1 images taken before and after the event completes this work. The results obtained demonstrate the importance of the complementarity of multi-sensor images and a multi-approach methodology to better monitor changes to a wetland area after a flash-flood disaster.

1. Introduction

Flash floods results in a rapid rise in the water level in a stream or a creek. Flash flooding can cause serious damage within a few hours (e.g., rupture of an embankment, dam failure). The intensity can increase dramatically, especially in cases where intense rainfall results in a rapid surge of rising flood waters. Greater runoff volume can result from various circumstances, such as a higher percentage of impermeable surfaces, compacted soils, soil moisture conditions, soil type and depth, terrain slope, land use, and vegetation, all of which can lead to short time lags between the rainfall occurrence and the peak discharge [1,2,3,4]. Violent rain leads to erosion that can be locally catastrophic. Traces of this hydric erosion in the form of incisions, ravines, and sediment deposits can be measured by imaging with a very high spatial resolution [5]. These extreme events tend to occur at very small spatial and temporal scales (generally affecting areas of up to a few hundred kilometers squared, lasting from minutes to a few hours, evolving very quickly and thus being difficult to forecast) [1,4]. These are some difficulties that are faced during observation and prediction at the regional scale, especially because of the absence of hydrometric networks in the desert environment. Remote sensing techniques can be a great help by providing safe and cost-effective tools for monitoring, mapping, and assessing the transformations and damages caused by flash-flood events [2,4]. During the rainy season, heavy rains cause floods every year in Morocco; these include flash floods, floods of rivers, and mud floods [2,3]. In October 2016, a torrential rain induced an impressive river flood in a large part of southern Morocco, especially in the city of Laayoune and its regions, which are located in the Moroccan Sahara. The event lasted for about 10 h (between 28 and 29 October). Further, the Sakia El Hamra dam, which has a storage capacity that does not exceed 112.6 Mm³, experienced a volume increase from 7 to 203 mm³ during the flash-flood event. The peak of the flood reached more than 3000 m³/s, which is well above the threshold of the dimensioned spillway capacity (410 m³/s), based on the information obtained from the Watershed Agency of Sakia El Hamra and Oued Eddahab. The passage of this exceptional flood caused the overflow of the wadi waters on the crest of the dam of Sakia El Hamra, leading to the degradation of “the downstream slope” and the opening of two breaches of a hundred meters in the body of the dam at the level of the minor bed of the Oued. The height of the two breccia was augmented considerably as the spill gradually approached and almost reached the level of the wadi (Figure 1).

The aim of this study is to map dynamic change and estimate the dynamics of inundation extent with spatially and temporally distributed water and wet/dry soil conditions during a local-scale flash-flood event using a collaboration between two approaches: one based on spectral processing and the other based on categorical processing. The first step was to extract four spectral indices and then detect the changes in multispectral diachronic images from three Sentinel-2 MSI (Mustispectral instrument) images and two Landsat-8 OLI (Operational Land Imager) images acquired prior to and after the event. The second step was performed using pattern classification techniques, namely, linear classifiers based on support vector machines (SVM) with Gaussian kernels [6,7,8]. In the first approach, we determined the extent of the flash flood using temporal and spatial surface changes of inundated areas according to water and wet/dry soil spectral indices. The best discriminating indices were selected from the four spectral indices, and the magnitude and orientation of the change vector were calculated and reclassified to extract the extent of the flood from each time point’s image acquisition and the soil moisture and dryness in these areas. In the literature, monitoring the changes in river response and wetland response to flash floods in a desert environment is very rare work. Our motivation for this study is that this research represents a chance to study the response of this desertic ecosystem to the effect of this very brief and rare event by using spatial data. It is for these reasons that satellite images from different sensors were selected for flash-flood mapping. It must be recognized that the acquired images may have been taken more than one day following the flood peak. This may influence the ground condition, making it entirely different from the flood situation.

2. Material and Methods

In the first part of this section, we present the study area and a detailed explanation of the dataset and the preprocessing step. Finally, we use a flow diagram to illustrate the methodology for the application of optical and synthetic aperture radar (SAR) data.

2.1. Study Area

Laayoune represents the largest city of Southern Morocco. This area is located between latitudes 27.00; 27.25 and longitudes −12.92; −13.38. The town of Laayoune is bounded from the west by the North Atlantic Ocean and is close to the Saquia El Hamra basin, which is a part of the geological unit of Tarfaya–Laayoune–Boujdour (Figure 2). Oued Saquia El Hamra occupies the central part of the Sahara with a wide opening to the Atlantic Ocean. Its area is over 142,865 km². In the rainy season, sedentary farming is based on underground water sources and dry flooding river beds. It is well known that flash floods are unpredictable and rare in the Sahara Desert [9]. Rainfall represents the most random variable among meteorological measurements in arid regions [10]. As a result, the frequency and severity of desert flash floods vary each year and remain undocumented [11]. Desert rainfall varies more spatially and is more random than rainfall in humid regions, so it is often described as irregular. The surface hydrological system’s response to desert rainfall is complex; thus, hydrological modeling is excluded [12]. Warner [11] argued that most desert floods occur because of the unusual nature of the surface rather than the rainfall, since the latter is unlikely to be much more intense than a similar storm type in more humid areas. The reasons for this are multiple [11,12,13]. (i) Very rare organic matter is present in the soil to absorb rainfall. (ii) Most desert areas are geologically predominated by impermeable rocks, which have a high potential for runoff, across all areas of the drainage basin. Consequently, if an abnormal amount of rainfall occurs in a very short time, it accumulates rapidly in flowing water in streams and rivers within flooding areas. (iii) The scarcity of vegetation implies that the raindrop can seal the surface of the soil. (iv) There is a lack of fauna, such as animals, insects, and worms, which make the soil relatively permeable. (v) The formation of runoff is strongly influenced by the structural state of the surface. The surface soil’s structure is degraded by precipitation, resulting in the formation of crusts that promote runoff.

The dunes of Laayoune are considered the fastest in the world, with an average speed of approximately 32 m/year, especially for dunes with a height of over nine meters [14,15]. Also, in the study area, there are two sabkha zones, which represent wet depressions with a saline surface in desert environments. They are generally considered to be floodable wetlands of great environmental interest, and they are strongly linked to numerous environmental issues, such as climate change, water quality, wildlife habitat, and biodiversity [16]. Sabkhas are supratidal, forming along arid coastlines, and are characterized by evaporate–carbonate deposits with some siliciclastic. They also form subaerial, programming, and shoaling-upward sequences that have an average thickness of a meter or less. Also, sabkhas are composed of more pockets of permanent water cavities and salt in the desert, where rainfall is known to be rare [16,17,18].

2.2. Materials

Our methodology is mainly based on determining the size of the area directly affected by floodwater and demonstrating how the combination of optical MSI and OLI sensors and radar data can contribute to deeper knowledge of the area under consideration. The data were processed to reveal and delineate the flooded areas. The Sentinel’s Application Platform (SNAP), the System for Automated Geoscientific Analyses (SAGA), the Geographic Resources Analysis Support System (GRASS), and Quantum GIS (QGIS), which are open-source common architecture for GIS and image processing Toolboxes, were used for the exploitation of earth observation (EO) data. Both optical image types were acquired from the web at https://glovis.usgs.gov/; Radar data was obtained from https://scihub.esa.int/.

2.3. Dataset and Preprocessing

The dataset used in the present research incorporates spaceborne images (

Table 1. Datasets used in this study.

Satellite	Instrument	Mode/Path-Row	Acquisition Date	Use
Sentinel 2 A	MSI	L1C TC SGS	13 October 2016	Two weeks before flash-flood event; used to calculate reference image.
Landsat 8	OLI	LC (105-41)	28 September 2016	One month before flash flood; used to calculate reference image.
Sentinel 2 A	MSI	L1C TL SGS	30 October 2016	One day after the flash-flood event; used for flood extent mapping.
Landsat 8	OLI	LC (205-41)	13 November 2016	Two weeks after flash flood; used for water and moisture/dryness soil information extraction.
Sentinel 2 A	MSI	L1C TL SGS	1 January 2017	Two months after flash-flood event; used for water and moisture/dryness soil information.
Sentinel 1 A	C-SAR	IW GRDH 1SDV(VV,VH)	21 October 2016	One week before flash-flood event.
Sentinel 1 B	C-SAR	IW GRDH 1SDV(VV,VH)	2 November 2016	Three days after flash-flood event.

), including:

• Three Sentinel-2 MSI images, one before and two after the flash-flood event, were used for extracting spectral indices
• Two Landsat-8 OLI images (path 205, row 41), one before and one after the flash-flood event, were used for extracting spectral indices
• Two Sentinel-1 (S1) SAR images, one before and one after the flash-flood event, were utilized for co-registration and classification to extract the flooding area

The bands within the solar reflectance spectral range were used to extract urban, vegetation, water, and wet/dry soil areas. Then, the original digital bands, OLI and MSI, were converted to an at-satellite minimum spectral reflectance value and at-satellite reflectance, and then further converted to surface reflectance [19,20,21,22]. All these data (Figure 3) were georeferenced to the WGS-84 datum and Universal Transverse Mercator zone 28N coordinate system. Each OLI/MSI dataset was clipped to reflect the extent of the study.

The preliminary treatment of the OLI sensor was subject to a significant radiometric and spectral characterization. Calibration represents a conversion to reflectivity and radiance of digital numbers according to Equation (1) for L8 OLI Top-of-Atmosphere (TOA) reflectance and Equation (2) for S2A MSI TOA reflectance for each band based on the metadata of the images [19]:

$${\rho }_{\mathsf{\lambda }}=\frac{M_{\rho }\times Q_{cal}+A_{\rho }}{sin\theta }$$

(1)

where

• ${\rho }_{\mathsf{\lambda }}$ = TOA reflectance,
• M_ρ = Reflectance multiplicative scaling factor for the band (obtained from the metadata file),
• A_ρ = Reflectance additive scaling factor for the band (obtained from the metadata file),
• Q_cal = Level 1-pixel value in digital number (DN),
• θ = Solar Elevation Angle (from the metadata file).

$${\rho }_{\mathsf{\lambda }}=\left(float\right)(\frac{Q_{cal}}{QUANTIFICATION\_VALUE})$$

(2)

• Q_cal = L1C pixel value in DN,
• QUANTIFICATION_VALUE is provided in the metadata file.

The multitemporal images needed to be radiometrically corrected before their further use in this study. This involved normalizing the change in brightness value (DN) because the varying atmospheric conditions and the sun’s position cause a variation in illumination. Radiometric normalization was done based on pseudo-invariant features (PIFs) in the images, which are objects that are spatially well defined and spectrally and radiometrically stable [23]. The mean reference DN values of pre-flood scenarios, i.e., the mean of OLI on 28 September 2016 and MSI on 13 October 2016, were selected as the references or base images for radiometric normalization of the other three post-flood images (MSI on 30 October 2016, OLI on 13 November 2016, and MSI on 1 January 2017). For objects with the same reflectivity on both the dates, the digital counts on each date are simple linear functions of the same variables and, therefore, are linear functions of each other. Equation (3) was used for radiometric normalization [23].

$$D{N}_{1i}=\frac{{\sigma }_{1i}}{{\sigma }_{2i}}D{N}_{2i}+({\mu }_{1i}-\frac{{\sigma }_{1i}}{{\sigma }_{2i}}{\mu }_{2i})$$

(3)

where

• DN_1i = ub>1i = Pixel value from date 1*,
• DN_2i = Pixel value from date 2**,
• σ_1i = Standard deviation (SD) of PIF from date 1*,
• σ_2i = Standard deviation (SD) of PIF from date 2**,
• μ_1i = Mean of PIF from date 1*,
• μ_2i = Mean of PIF from date 2**,
• i = Band number (2,3,4,5,6,7) for OLI and (2,3,4,8A,11,12) for MSI,
• 1* = Mean of OLI and MSI before flood event (reference image),
• 2** = After flood image.

The mean and standard deviation (SD) of the PIFs of each band of images on two dates were retrieved as an example; these values were used with the above formula, resulting in band-to-band radiometric normalization (

Table 2. Mean and standard deviation of Sentinel-2 and Landsat-8 satellite datasets before and after normalization.

	Band before Normalization		Band after Normalization
	Mean	SD	Mean	SD
Reference before flood 1	1956.47	95.21	2050.66	105.14
Reference before flood 2	2364.15	145.78	2390.44	155.55
Reference before flood 3	3291.65	237.01	3335.65	255.66
Reference before flood 4	4050.43	224.13	4066.77	234.66
Reference before flood 5	4676.46	288.21	4705.33	299.77
Reference before flood 6	3558.22	275.94	3588.88	285.98
MSI one day after flood 2	1949.74	92.74	1966.65	100.79
MSI one day after flood 3	2358.06	141.98	2377.09	156.88
MSI one day after flood 4	3286.00	231.30	3289.88	244.54
MSI one day after flood 8A	4034.24	218.96	4055.51	277.77
MSI one day after flood 11	4655.17	266.52	4675.66	288.76
MSI one day after flood 12	3543.82	258.24	3588.74	278.62

).

2.4. Methodology

Optical spatial data enable change detection using several methods, ranging from classic visual-based photointerpretation procedures to semi-automated and automated methods based on spectral characteristics. The success of change detection depends on the image processing and classification procedures employed [24,25,26,27,28,29,30,31,32]. To characterize them, optical sensors use the emission and reflective properties of materials. Optical images are feature-based and enable a fluent interpretation. In addition, a significant number of optical sensors have been used for decades to successfully map floods, as optical satellite sensors can detect floods from space using water’s unique spectral characteristics. [25,29] Although a desert wetland environment presents a multiplicity of spectrally similar features, the near-infrared bands show clear distinctions between these land cover types [26]. In flooded wetlands, open water and flooded vegetation usually absorb near-infrared radiation and are represented in optical images by low DN values, while non-flooded classes emit infrared radiation and are represented by higher DN values [25,26].

Using spectral processing to acquire indices necessitates the analysis of the spectral characteristics of remote sensing imagery to determine their quantitative relationship with wetness/water and dryness in land changes. Albedo, normalized difference water index (NDWI), normalized multi-band drought index (NMDI), and normalized difference moisture index (NDMI) are the principal spectral indices used to characterize moisture in land changes [3,4,33]. Hydrological processing involves acquiring indices to analyze the hydrological characteristics of a digital elevation model (DEM) to detect areas associated with flash floods and erosion caused by rainstorms and sediment transport and accumulation [4,25]. However, even with the increased number of satellite systems and sensors acquiring data with improved spectral, spatial, radiometric, and temporal characteristics, most existing land-cover datasets have been derived from categorical processing based on pixel or object approaches using conventional or advanced pattern recognition techniques, such as random forests (RFs) [34,35], neural networks (NNs) [6,36,37,38,39,40,41,42,43,44], and support vector machines (SVMs) [6,7,8,45,46,47,48]. It is obvious that the collaboration of all these processing methods can enhance mapping by using the multispectral, multitemporal, and spatial information available to derive the best characterization of wet/water and dry land changes from remote sensing images after a flash-flooding event [49,50,51,52,53,54].

To map flooded areas, different remote sensing data and methodologies were used in accordance with the following criteria: (I) multi-sensor data (SAR, multispectral); (II) free (zero-cost) imagery (i.e., we used only free-of-charge data); (III) the availability of data in relation to the flood phase (pre-flood, co-flood, and post-flood); (IV) different spatial resolution of SAR and optical data; and (V) the availability of support data (DEM or other ancillary data). We collected (a) pre-flood multispectral (OLI and MSI) data acquired before the inundation and used their mean as a reference for the change detection analysis; (b) co-flood data, collected after one day of the maximum inundation (MSI) data; (c) multidate post-flood (optical and SAR) data, taken after the flood for mapping dynamic change with temporally distributed water and wet/dry soil conditions [55].

2.4.1. Multispectral Data Methodology

Consider two multispectral images acquired over the same geographical area at two different times. They are expressed as T1 and T2. T1 is the reference image and was obtained by calculating the mean of two images acquired before the flood event from OLI and MSI sensors; T2 corresponds to one of three multitemporal images that were acquired at different times after the flood from OLI and MSI sensors. Then, as shown in Figure 4, the proposed method has six main steps: (1) radiometric correction and normalization; (2) spectral processing; (3) categorical processing; (4) radiometric change detection; (5) categorical change detection; and (6) collaborative change detection.

Spectral Processing and Indices’ Feature Space Choice

Albedo, normalized difference water index (NDWI), normalized multi-band drought index (NMDI), and normalized difference moisture index (NDMI) are the principal indices that are used to characterize moisture/dryness land change [49,50,51,52,53,54,55,56,57,58,59]. The albedo index (Equations (4) and (5)) represents a physical parameter that is affected by land surface conditions, such as soil moisture, vegetation coverage, open water, and dry land. Therefore, albedo (“(4)” and “(5)”) can be used as an important physical parameter reflecting the land’s wet/dry state [49,50,51].

$$Al=\left(0.356{\rho }_b+0.130{\rho }_r+0.373{\rho }_{nir}+0.085{\rho }_{SWIR-1}+0.072{\rho }_{SWIR-2}-0.0018\right)$$

(4)

where ρ_b: blue band (band 2 for OLI or MSI), ρ_nir: near-infrared band (band 5 for OLI or band 8A for MSI), ρr: red band (band 4 for OLI or MSI), ρ_swir-1: swir1 (short wave infrared 1) band (band 6 for OLI or band 11 for MSI), ρ_swir-2: swir2 (short wave infrared 2) band (band 7 for OLI or band 12 for MSI).

For the data normalization process, we must first find the maximum and minimum values of albedo in the study area and then use these values to normalize the data and obtain the albedo normalized index.

$$Albed{o}_n=\left(Al-A{l}_{min}\right)/\left(A{l}_{max}-A{l}_{min}\right)$$

(5)

The NDWI (Equations (6) and (7)) is also used for soil water mapping after periods of flash flood [53,54,55]. Xu [53] similarly concluded that the NDWI could be used to detect the flooding of rice fields in China.

$$NDWI=\left({\rho }_b-{\rho }_{swir-1}\right)/\left({\rho }_b-{\rho }_{swir-1}\right)$$

(6)

The normalization of this index is as follows:

$$NDW{I}_n=\left(NDWI-NDW{I}_{min}\right)/\left(NDW{I}_{max}-NDW{I}_{min}\right)$$

(7)

Furthermore, a new drought index, namely, NMDI (Equations (8) and (9)), was proposed [57] using three bands in the near infrared (NIR) and shortwave infrared (SWIR) for the remote sensing of soil and soil water content (0.86 μm, 1.64 μm, and 2.13 μm):

$$NMDI=700\left({\rho }_{nir}-{\rho }_{swir-1}+{\rho }_{swir-2}\right)/\left({\rho }_{nir}+{\rho }_{swir-1}-{\rho }_{swir-2}\right)$$

(8)

$$NMD{I}_n=\left(NMDI-NMD{I}_{min}\right)/\left(NMD{I}_{max}-NMD{I}_{min}\right)$$

(9)

The NDMI (Equations (10) and (11)) is an index for surface reflectance data [58,59] that uses all six bands, and the coefficients are empirically determined as those that transform a set of training pixels into new dimensions aligned with maximum variability.

$$NDMI=\left({\rho }_{nir}-{\rho }_{swir-1}\right)/\left({\rho }_{nir}+{\rho }_{swir-1}\right)$$

(10)

The normalized index is:

$$NDM{I}_n=\left(NDMI-NDM{I}_{min}\right)/\left(NDM{I}_{max}-NDM{I}_{min}\right)$$

(11)

The images for all water indices’, including NDWI_n and NDMIn, definitely help to improve the separability of the water bodies. Most NDMIn values of water bodies are larger than 0.76, while most NDWI_n values of water bodies are larger than 0.5. As shown in Figure 5a,c, the NDWI_n values of built-up areas are clearly different. Compared with water bodies, vegetation has much smaller NDWI_n values; thus, this index facilitates the distinction between vegetation and water bodies.

Indeed, the built-up characteristics in the NDWI_n image appear in a light-yellow tone with positive values between zero and 0.2. This leads to confusion, particularly in the urban areas, between built-up locations and soil.

After calculating albedon, NDWI_n, NDMIn, and NMDIn (Figure 5), we stacked them in one image file and calculated the correlation matrix to choose the less correlated indices, which were found to be NDMI and albedo, whose relationship has a correlation coefficient of −0.486 (

Table 3. Correlation Matrix.

	NDMIn	Albedon	NDWIn	NMDIn
NDMIn	1	−0.49	0.84	0.88
Albedon	−0.49	1	0.75	0.72
NDWIn	0.84	0.75	1	0.97
NMDIn	0.88	0.72	0.97	1

). This value was confirmed in the feature space (NDMI, albedo). Thus, the NDMI and albedo were chosen as the feature space to better extract the area changes.

As described above, the analysis shows that there is a significant correlation between wet/dry soil, NDMI, and albedo. Figure 6 summarizes the differences in albedo as soil wetness/dryness varies at different phases after flash floods. The figure is based on the Δ (NDMI) and Δ (Albedo) feature space.

Reducing soil wetness leads to a significant increase in albedo, and albedo changes affect the surface radiation balance. In the (Δ (Albedo), Δ (NDMI)) feature space, albedo depends not only on soil moisture but also on soil surface temperature. The line connecting A and B in Figure 6 is considered a high albedo line, and reflects the drought and the high albedo value on the soil surface, which is the upper limit of high albedo, and corresponds to the complete arid land. The line connecting C and D is the low albedo line, representing sufficient soil wetness.

To further study the spatial distribution law of the flash-flood event in the (Δ (NDMI), Δ (Albedo)) feature space, we present a scatterplot for different times after the event (Figure 7).

After statistical regression analysis, we find that Δ (Albedo) and Δ (NDMI) have a significant negative linear correlation, as indicated in Equation (12). The regression equation is:

$$\Delta \left(Albedo\right)=230.426-1.035\Delta \left(NDMI\right)$$

(12)

The correlation coefficient is −0.486. This indicates that Δ(NDMI) and Δ (Albedo) have a moderate negative correlation: i.e., they have a weak degree of association, so a growing statistical trend in soil dryness results in a gradually decreasing statistical trend in NDMI, while the albedo statistical trend increases. In the (Δ (NDMI), Δ (Albedo)) feature space, the wet/dry soil is reflected clearly. According to the study by Verstraete and Pinty [60], different wet/dry soils can be effectively separated by dividing the (Δ (NDMI), Δ (Albedo)) feature space in the principal axis of the cloud in the direction of changing trends (Figure 8).

The method used in this study is based on the analysis of the relationship between the spectral indices of wet soil and dry soil. The idea of change vector analysis (CVA) is that some pixels with different values persist over time at substantially different locations in the feature space. However, the best thematic discriminant feature space is the determinant [61,62,63]. Changes in wet and dry land between pre and post-flash-flood events were determined using CVA with the NDMI and albedo indices. CVA is a multivariate technique [61,62,63,64]. This enables the pixel-by-pixel analysis of spectral bands or products on two axes of a Cartesian plane to evaluate changes between pixel positions (T1) in the reference image and the same pixel (T2) in the post-flood image (Figure 7). The magnitude (M) indicates the intensity of the change and is derived using the Euclidian distance in Equation (9). The result is an image classified into one of four categories: low (0–25), medium (25–50), high (50–75), and extreme (75–100), where the values are the sum of the changes that occurred between the dates.

The aim of our method is to find the best feature space that can best discriminate changes between wet and dry land. To this end, this study aims to find significant negative correlation indices. In this significant correlation feature space, we calculate the magnitude and orientation of the change vector (CV) as follows (Equations (13) and (14)):

$$M=\sqrt{{\left(Albedo1-Albedo2\right)}^2+{\left(NDMI1-NDMI2\right)}^2}$$

(13)

where sub-indices 1 and 2 represent images taken before and after the flash flood, respectively.

The direction of change represents the type of change (θ) detected in a pixel corresponding to T1 and T2; that is, each vector is a function of the combination of positive and negative changes occurring in the spectral band or products analyzed, with θ determined according to Equation (10). The resulting image includes values between 90–180°, representing an increase in drought, or between 270–360°, supporting bright pixels indicating an increase in wetness; values in the ranges 0 to 90° and 180 to 270° represent pixels illustrating no change; therefore, there is no decisive evidence that drought or wetness will increase or decrease.

$$\mathsf{\theta }={\tan }^{-1}\frac{\left(\mathrm{NDMI}1-\mathrm{NDMI}2\right)}{\left(\mathrm{Albedo}1-\mathrm{Albedo}2\right)}$$

(14)

Categorical Processing with SVM

SVM is a categorical processing method that is based on reducing the problem of discrimination to that of searching for the optimal linear hyperplane [6,7,8,45,46,47,48]. Two ideas or tips help achieve this goal:

-

Define the hyperplane as a solution to an optimization problem that separates the different classes of data under constraints whose objective function is expressed using scalar products between vectors and in which the number of “active” constraints or points lying on the boundaries are called support vectors; the middle of the margin is the optimal separating hyperplane that controls the complexity of the model.

-

The selection of the kernel is the most important task in the implementation and the success of the SVM classifier [47,48]. A kernel function is introduced into the scalar product to search for nonlinear separating surfaces implicitly, and the nonlinear transformation of data to a larger feature space is induced. An optimum hyperplane is determined using training data points or the support vector of each class on the hyperplane to be maximized, and its generalization ability is verified using a validation data point. Training vectors, xi, are projected into a higher-dimensional space by the kernel function. The SVM can work with a few training data points; such is the case in this study, with classes having less than 45 samples. The kernel function used in this experiment is a radial basis function (RBF)—defined by Equation (15)—for two input vectors, x_i and x_j, with the maximal margin γ parameter, which is user-defined. The choice of RBF kernel is justified by the frequent use of this method in the literature in the field of image classification and in change detection studies [45,46]. This kernel requires only one parameter to be pre-defined, which makes it easy to implement in contrast to other kernels.

$$K\left(x_i,x_j\right)=e^{-\gamma \Vert x_i-x_j{\Vert }^2},\mathsf{\gamma }>0$$

(15)

2.4.2. SAR Radar Data Methodology

SAR polarization is a determinant factor in flood detection. It has been proven in several previous studies that HH (horizontal/horizontal) -polarized images are more appropriate for flood detection than VV (vertical/vertical)-polarized or cross-polarized images [65,66,67,68,69,70,71,72,73,74,75]. This is mainly because HH polarization gives the best discrimination of backscatter values between dry soil and wet soil and between smooth soil and rough soil, so it enables improved monitoring of the flooded area after an event [65,71]. However, for this study, the available image during the flood has intensity VV and VH polarization. VV polarization was deemed preferable, because the medium incident angle of the data makes this image suitable for flood monitoring and allows for the effective detection of wet areas. Moreover, the same polarization (VV) was chosen for both dates to achieve co-registration for comparing the same co-polarization (VV–VV). This strategy gives better results than using cross-polarization (VH–VV).

Radar-imaging SAR uses the target surface roughness and depicts it in images. Water has a distinctive dielectric property, particularly stagnant water with a specular reflector; therefore, it returns a low radar signal and displays a low backscattering signal in radar images [76,77].

The ability to penetrate cloud cover and detect water makes the satellite SAR system a powerful tool for flash-flood monitoring, as flash floods are often associated with heavy rain, and it is almost impossible for optical remote sensing instruments to image the water in such circumstances. SAR has great potential as a source of relevant and near-real-time information. It represents a unique chance because of the time range of the optical sensor image, which is extremely close to the real-time acquisition of this very brief event. This is the reason for selecting satellite images from different sensors for flash-flood mapping: the images acquired have been taken more than one day after the flood peak. The efficiency of imaging radar SAR in river and wetland flood mapping depends on the wavelength of the sensor and the roughness of the surface of the wetland [78,79,80]. Therefore, long-wavelength SAR better characterizes wetland floods than short-wavelength SAR. However, the finite availability of long-wavelength SARs is paradoxically a real limitation on the use of radar SAR to map and monitor floods in rivers and wetlands [77]. Finally, it may well be that there are limitations for both radar and optical data. However, working with a collaborative approach can help to combine the strong points of different sensor types. For instance, data with a high spatial resolution and low repeat cycle can be combined with data that have a low spatial resolution with high repeat acquisition to create a flood product with high resolution and frequent observations [73,74,75]. It is the association of radar and optical cloud-free data that enables the production of precise maps of both river and wetland flood monitoring. With the results of each sensor, this can be done by combining spectral and categorical processing in a collaborative change detection process. This provides complementary information and increases the accuracy of monitoring. Due to the rapid dynamic nature of this event, temporal and spatial resolutions also influence flash-flood mapping accuracies, so multitemporal, multi-resolution satellite data are needed to monitor floods in a desert river and desert wetland ecosystem.

Specific steps were taken to process the data for radar images (Figure 9). The first step is related to the preprocessing of the data (mosaicking, calibration, co-registration, band stacking). A subset process of images was created to reduce the data size and make it easier to manage. Finally, a speckle filter (Lee filter with a window size of 7 × 7) was applied. Speckle appears in an image as a grainy “salt and pepper” texture. This is caused by random constructive and destructive interference from the multiple returns of scattering occurring within each resolution cell. The second step deals with the main processing of the data. As far as the first method is concerned, for the separation of the inundated areas and the monitoring of wet areas, a supervised classification approach was selected. First, the histogram of the filtered backscatter coefficient of each band in the stacked image was analyzed. The histogram shows a peak of a different magnitude. Low values of the backscatter correspond to the water, and high values correspond to the non-water class. For the water and wet area extraction, we applied the supervised image with a combined textural and intensity classification approach [81,82].

To assess the classification accuracy, it has been suggested that at least 216 randomly selected pixels approximately of each land-cover class should be checked against ground truth verification data to assess the accuracy of the classified images [83,84]. The overall accuracy, user’s accuracy, overall disagreement, quantity disagreement, and allocation disagreement have been used as measures of accuracy [84].

The pair of SAR images obtained from the Sentinel-1C band comprises an image three days following the flood event, which is called the “flood image”, and another one before the flood event, which is called the “reference image”. To distinguish the flooded area, co-registration processing [62,63,64,65,66,84] was applied.

3. Results and Discussion

In the first part of this section, we present the information from monitoring the soil moisture/dryness to detect areas changed by flash floods according to the optical data. Afterward, we present the information obtained from SAR data. Next, we report the results of monitoring the water extent after this flash flood. Finally, we present a discussion about these results.

3.1. Change Mapping Obtained from Optical Data

3.1.1. Characterization of Changed Areas Caused by Flash Floods Using Δ (Albedo) and Δ (NDMI) Change Indices

The change indices Δ (Albedo) and Δ (NDMI) create an output of detected change, where the values vary in [−1;1] and are characterized by the three-mode histograms in Figure 10. The main mode corresponds to the no-change status in the image (value of around 0), and the yellow and cyan areas in the change image and the two secondary modes refer to changes. The first secondary mode in red (Δ (Albedo) ≈ −1 and Δ (NDMI) ≈ −1) corresponds to negative changes (infilling or submersion), while the subsequent secondary mode in blue (Δ (Albedo) ≈ 1 and Δ (NDMI) ≈ 1) corresponds to positive changes (draw-off or emergence). The water submersion (negative changes) appears in red on the change index displayed in the indexed color, while emerging areas (positive changes) appear in blue in Figure 11.

3.1.2. Qualitative Change Analysis by the Color Composite of the (Δ (Albedo), Δ (NDMI)) Feature Space

Figure 12a–c shows the color composite RGB at different times after the flood event, with red indicating Δ (Albedo) and green and blue indicating Δ (NDMI). In these images, the areas in red represent the zones in which the albedo underwent the maximum change between the two dates—before and multiple time points after flooding—and zones in light cyan represent where the NDMI underwent the maximum change between the two dates. It is the riverbed and the sabkha areas whose surfaces changed over time after floods.

3.1.3. Categorical Change Mapping with SVM

The areas in the image can be classified into six categories: water, land, sabkha, building, agricultural area, and sandstorm. Figure 13a–c shows three time points: before the flood event, one day after the event, and two weeks after the event. Visually, the overall performance of the SVM in land-cover classification is good, as it can classify all the pixels effectively. Accuracy analysis was carried out by comparing the classified pixels with real ground truth pixels using a confusion matrix [83,84]. For accuracy assessment purposes, the selection of the ground truth pixels was done by random sampling with equal sample allocation in each stratum (216 randomly selected pixels approximatively of each land-cover class). The overall accuracy of the categorizations using this data subset is 86.45% (see

Table 4. The error matrix showing the overall accuracy (OA) performance expressed in pixel counts.

	Class							Accuracy Assessment %
Reference (Ground Truth)	Building	Sabkha	Sandstorm	Sand	Agriculture	Water	$n_{j+}$	$O_{Er}$	$U_i$	$P_i$	OA
Building	191	3	5	3	2	10	211	10.9	89.1	87.3	86.45
Sabkha	2	186	5	2	10	6	209	12	88	84.7
Sandstorm	7	13	180	4	4	7	211	16.6	83.4	77.8
Sand	6	3	12	188	2	8	219	12.3	87.7	88
Agriculture	3	7	7	6	192	10	219	15	85	88
Water	9	7	19	13	8	180	223	25.1	74.9	81.5
$n_{+j}$	212	216	216	216	216	221	1292
$C_{er}$	12.7	15.3	22.2	12	12	18.5

). This value is weakly acceptable if we consider the difficulties of the studied region and that we used only 221 pixels to train the algorithm. The percentages of correct classifications of pixels as building, sabkha, sandstorm, sand, agriculture, and water are 87.3%, 84.7%, 77.8%, 88%, 88%, and 82.5%, respectively. For these three classifications, the agreement coefficient value is between 88.8–92.5%, Quantity disagreement is between 2.66–3.78, allocation disagreement is between 3.27–4.5, and finally, disagreement is between 5.93–8.28. All these values show that the SVM classification and the real ground-truth pixel data are in relatively good agreement [84].

The map categories (classified pixel by class) are the columns, while the reference categories (reference pixel by class) are the rows; n_i+ and n_+j are the total of the preceding rows and columns, respectively. Abbreviations: Omission Error (O_Er), Commission Error (C_Er), User’s accuracy (U_i), producer’s accuracy (P_i), and overall accuracy (OA).

3.1.4. Monitoring of Moisture/Dryness Wetland Information by Fusion of Spectral Processing and Categorical Processing

The results of the SVM pattern classification before the flash-flood event and each range of time following it were used to extract the limit of the change area. Afterwards, we used the segmentation of the direction (using Equation (10)) and magnitude (using Equation (9)) of the changes that occurred in these change areas given by SVM classification differences between pre-flood and post-flood (at different time points) images. The areas with no change in SVM classification were not used in the segmentation of the CVA process results. Figure 14a,b, and c show the segmentation of the direction (using Equation (10)) and magnitude (using Equation (9)) of the changes that happened in the study area from before the event to the several periods used after the event. There was no change in 93% of the area; 3.5% of the zone presents some degree of wetness, and a change of 3.2% occurred in the drought area, corresponding to the river and sabkha at one day, two weeks, and two months after the event. The magnitude of the soil dryness change in the study area was classified by three categories: soil with low dryness, soil with medium dryness, and soil with high dryness. The magnitude of change in soil wetness was classified according to three categories: soil with high wetness, soil with very high wetness, and open water.

It can be seen from Figure 14 that the main subjects of interest, namely, the river and the sabkha, both show variation in the humidification scale of their surfaces only one day after the floods. Two weeks following the event, only a part of the river remains humid. The observation demonstrates that the southern sabkha presents a drier surface than before the floods occurred, except at the area limits between the sand dune and the edges of this sabkha. Furthermore, the sabkha located in the north remains wet in the majority of its area. Two months following the flooding, we notice that the river returns to its normal state, and that the southern sabkha is noticeably drier, except in the northern part, which is wet over half of its surface. In summary, the collaborative approach combines the benefits of the two process: categorical processing defines the precise edges of the major change classes and cannot discriminate between types of changes; in spectral processing, the pixel values in the original bands are treated as vectors of transformed difference data in addition to a change vector (CV) calculated by subtracting the vectors for all the pixels on different dates. The direction of the CV depicts the type of change, whereas the magnitude of the change corresponds to the length of the CV [63].

3.2. Change Mapping Obtained from SAR Data

The flood extent map obtained from this method (co-registration, classification) was compared and evaluated. The map of changes resulting from the supervised classification, which is based on a combination of the textural and intensity approaches, is demonstrated in Figure 15b. In the figure, the inundated area composed of very wet soil appears in red, which represents areas with increased backscattering due to the considerable level of soil wetness. Areas in blue are the flooded areas that remain covered with water three days after the floods; these are flooded agricultural areas with decreased backscattering. The total water extent was calculated to be 36.2 km² (3620 hectares). The regular water body in the dam (Figure 15a) before the flood event is 12 km² (1200 hectares). Consequently, the flooded area for this time is 24.2 km² (2420 hectares). This classification of the co-registered Sentinel-1 channel results in an overall accuracy of 0.90, and the agreement coefficient value is 89.9%.

3.3. Monitoring River Water Extent after Flash Flood

To be able to extract the water extent, all areas with wet soil were merged and grouped into one water class. The pre-flood extent of a regular water body in the dam (Figure 16a) is 12 km². During the flood peak, one day following the event, the flooded area is 29.44 km² (Figure 16b). Two weeks following the event, the flooded area is 8.48 km² (Figure 16c). Finally, two months after the event, the flooded area is close to zero: 0.35 km² indicating a return to a regular water body in the dam (Figure 16d). The figure represents the variation in the flood extent obtained by the mean of MSI and OLI before the flood event as a reference, MSI one day after the event, OLI two weeks after the event, and MSI two months after the event. In the radar image of the region three days after the flood, the areas in blue, which correspond to a decrease in backscatter, can be designated with certainty as water surfaces. However, the red areas that correspond to an increase in backscatter cannot be designated with certainty as water surfaces, because it could be an increase in surface roughness. Finally, it can be seen that the results obtained from the optical images are in agreement with those obtained with the radar sensor. However, we cannot affirm the presence of water with certainty in the areas with increased backscattering, because it could be due to an increase in the surface roughness after the passage of water in the zone.

3.4. Monitoring Sabkha Water Extent and Soil Wetness after Flash Flood

The extent of a regular water body in the sabkha (Figure 17a) before the flood event is 0.05 km², which represents the permanent water cavity in the southern sabkha. During the flood peak one day following the event (Figure 17b), the extent is 8.35 km² for the southern sabkha and 9.64 km² for the northern sabkha. Three days after the event (Figure 17c), according to the SAR data, the extent is 6.2 km² for the southern sabkha and 5.50 km² for the one in the north. Two weeks following the event (Figure 17d), the extent of the wet area was 1.48 km² for the southern sabkha and 5.75 km² for the one in the north. Finally, two months after the event (Figure 17e), the extent of the wet area is 6.8 km²; it is close to the extent before the event for the southern sabkha, indicating a return to the regular water body in the permanent water cavity of the sabkha. However, the one in the north remains wet with an area of 2.75 km². Thus, we can conclude that the responses of the two sabkhas to this storm are not the same. The southern sabkha has larger transmission losses than the north.

3.5. Discussion

Arid wetlands are important for the conservation of biodiversity, but they are also highly sensitive and vulnerable to climate variability and extreme hydroclimatic events [76,85,86,87]. The amplification of the water cycle, along with the rising frequency and severity of extreme drought and moisture events, are expected to change the spatiotemporal hydrological patterns of wetlands, with potential threats to the functions of ecosystem services [78,85]. Here, by using spatial multi-sensor and multitemporal data, we studied the response to the flash-flood event of two desert wetlands: the river and the sabkha systems. Our observations can be summarized as follows:

In flooded wetlands, open water and flooded vegetation usually absorb near-infrared radiation and are represented by low DN values in optical images, while non-flooded areas generate higher DN values due to the infrared radiation emission. [26,85]. During floods in desert wetland environments and areas with poor vegetation cover, sediments mixed with water are very abundant. These materials have higher reflectance than clear water and can lead to misclassifications. To extract the extent of the water and compensate for these deficiencies, we are interested in the pixel value changes that are caused by this event, instead of focusing on the value of the pixel itself. So, we determined information about the extent of the flash flood using the temporal and spatial surface changes in inundated areas by using water and soil wet/dry spectral indices. The best discriminant indices were selected from four spectral indices to calculate and segment the magnitude and orientation of the CV to extract the extent of the flood in each image acquisition and monitor the soil moisture and dryness in these areas.

The limitations of the presented results are the low availability of optical images in accordance with a high probability of cloud cover during the flash-flooding period. This represents a constraint when using optical sensors rather than SAR sensors for mapping and monitoring river and wetland floods. However, our study is an exception to this rule, since the availability of cloud-free optical data was temporally closer to the event than the SAR data. This may be a consequence of the very brief nature of this event, as well as the good temporal resolution of the Sentinel-2 sensor.

The collaborative approach combines the benefits of the two processes: categorical processing provides the precise edges of the major change classes and cannot discriminate between types of dry/wet changes; with spectral processing, the pixel values in the original bands are treated as vectors of transformed difference data in addition to the change vector (CV), which is calculated by subtracting vectors for all the pixels on different dates. The direction of the CV depicts the type of change, whereas the magnitude of the change corresponds to the length of the CV [63]. So, the (Δ (Albedo), Δ (NDMI)) feature space derived from multi-temporal optical remote sensing imagery allows for the determination of radiometric change, while categorical processing allows for the determination of thematic change.

The proposed method for SAR radar data strongly improves the reliability of flood classifications in arid regions, where sand surfaces generally lead to strong overestimations of the water extent with SAR data [28]. In the co-registered image, areas that correspond to a decrease in backscatter can be designated as water surfaces with certainty. However, areas that correspond to an increase in backscatter cannot be designated as water surfaces with certainty, because it could be caused by an increase in surface roughness due to the transport of materials by river floodwaters. Also, this study shows that the flood extent estimated from radar images is in spatiotemporal continuity with the flood extent estimated from optical images if we consider that the zones of increased backscattering correspond to areas that are flooded.

The southern sabkha has larger transmission losses than the north because of its relatively low depth, as shown in Figure 18, which shows that the northern sabkha has a greater depth than the south, and therefore has a greater capacity to retain rainwater than southern sabkha. There is substantial evaporation, since the sabkha that lies to the north is directly exposed to the cold wind current coming from the north, as seen in the optical satellite images. All these effects make the terminal water remaining after this event decrease rapidly, so this sabkha dries more quickly.

4. Conclusions

This research focused on monitoring the water extent in a wetland after flash floods in the Saharan area. Using collaborative spatial data, we showed the evolution of spatially distributed water in this area after an inundation event from a local-scale flash-flood event. This study used optical data, radar SAR data, and collaborative approaches by combining spectral and categorical processing. Spectral indices derived from multi-temporal optical remote sensing imagery were used to show radiometric changes, while categorical processing was used to show thematic changes. Finally, the fusion of the two processes resulted in the above-presented map of changes. Spectral processing includes spectral indices such as the albedo, the NDWI, the NDMI, and the NMDI, describing the soil reflectance curve.

However, in the desert environment, optical data become as important as radar data and give better results. This work shows that the return to normal water after flooding in an arid environment took approximately two months, and the sabkha’s normal level of water was briefly reduced, resulting in the sabkha becoming drier over time. Also, our methodology determined the size of the area that is directly affected by flash floodwater. It further demonstrated that the evolution of the extent of the floods extracted from the radar image present a similarity with the extent of the floods extracted from the optical images. In addition, this study monitored the moisture condition of the surface following a flash flood in an arid region, with the aim of presenting the underlying causes of developing countries in risks management.

In fact, the results can help improve the management of water regulation structures by developing a methodology for maximizing water storage capacity and reducing the risks of flash floods in Laayoune and its regions. Finally, the achieved results show that the contribution of SAR images to flood mapping is complementary to the optical images, and that merging the spectral and categorical results of multi-sensors increases the precision and the continuity, thus contributing to better desert wetland monitoring.