Urban Flood Detection Using TerraSAR-X and SAR Simulated Reflectivity Maps

Abstract

Synthetic Aperture Radar (SAR) imagery is a vital tool for flood mapping due to its capability to acquire images day and night in almost any weather and to penetrate through cloud cover. In rural areas, SAR backscatter intensity can be used to detect flooded areas accurately; however, the complexity of urban structures makes flood mapping in urban areas a challenging task. In this study, we examine the synergistic use of SAR simulated reflectivity maps and Polarimetric and Interferometric SAR (PolInSAR) features in the improvement of flood mapping in urban environments. We propose a machine learning model employing simulated and PolInSAR features derived from TerraSAR-X images along with five auxiliary features, namely elevation, slope, aspect, distance from the river, and land-use/land-cover that are well-known to contribute to flood mapping. A total of 2450 data points have been used to build and evaluate the model over four different areas with different vegetation and urban density. The results indicated that by using PolInSAR and SAR simulated reflectivity maps together with five auxiliary features, a classification overall accuracy of 93.1% in urban areas was obtained, representing a 9.6% improvement over using the five auxiliary features alone.

1. Introduction

1.1. General

Flooding is among the most devastating natural disasters on earth and is the most frequent hazard in Canada, mainly due to heavy rainfall or rapid snow melt [1]. With climate change, flooding is expected to increase worldwide, particularly in high latitude regions like Canada [2,3,4]. Thus, the level of destruction caused by floods is high, especially in urban areas where damage to infrastructure leads to irreversible losses.

Remote sensing has shown a great potential to detect flooded areas thanks to optical and Synthetic Aperture Radar (SAR) satellite images [2,5,6,7,8,9,10]. However, flood situations are very likely to be cloudy, a condition in which passive optical remote sensing cannot be used. On the other hand, SAR remote sensing can penetrate cloud cover in nearly all-weather conditions, regardless of the time of day or night, which makes them useful in flood situations.

In recent years, SAR data has been recognized as a powerful tool for flood mapping [1,6,11,12,13]. SAR satellite images such as Sentinel-1A [5,6,14], TerraSAR-X [7,11,12,15], COSMO-SkyMed [13,16], RADARSAT-2 [17], Radarsat Constellation Mission [1], and ALOS-2/PALSAR-2 [18] have been employed to improve flood detection. Furthermore, sensors with a higher spatial resolution (3 m or more) such as TerraSAR-X, RADARSAT-2, and COSMO-Sky Med are well suited for urban areas [7,13,17,19].

Several studies have indicated that SAR imagery can perform better than optical imagery in flood mapping [20,21]. In rural flood mapping, the SAR backscatter intensity is a key factor in distinguishing flooded and non-flooded areas. Thus, SAR-based flood detection in rural areas has been extensively explored [12,15,22,23]. On the contrary, SAR-based flood detection in urban areas is challenging due to various complex urban backscatter patterns, including double bounce, shadow, and layover, which cause misclassification and increase false alarms [5,7,11,13,16,17,18]. Because of this issue, the majority of SAR-based studies in the literature have focused on rural flood detection and left urban flood detection largely unexplored.

A limited number of studies in the literature have focused on using SAR images for urban flood detection. These studies demonstrated that the use of multi-temporal Interferometric SAR (InSAR) coherence together with SAR intensity works efficiently in urban flood mapping [6,7,14,16,24,25,26]. Multi-temporal InSAR coherences have been found to be effective for flood detection even when using medium-resolution SAR images like Sentinel-1A [27]. InSAR coherence is the normalized cross-correlation between two interferometric SAR images, and it is sensitive to any variation in the geometric properties of the objects and thus to flood-induced changes. InSAR coherence is however sensitive to variations of the geometric properties of the scatterers, and it has been suggested to be considered as complementary information along with other data for flood detection [7]. Additionally, the authors of this paper showed that incorporating the multi-temporal InSAR phase improves flood detection in urban areas [28]. InSAR phase is assumed to remain unchanged as long as no changes are made to the spatial arrangement of objects [29]. We, therefore, believe that the InSAR phase can be used as a feature in flood detection.

Polarimetric SAR (PolSAR) decompositions are also a powerful tool to detect complicated urban backscatter patterns as they present the average scattering mechanism which is associated with the physical properties of an object [30]. However, it is not feasible to construct polarimetric decompositions for single polarized images as polarimetric decompositions require fully polarized or at least dual polarized data. Additionally, PolSAR filtering techniques can preserve the polarimetric properties, which have been demonstrated to be affected by any changes to the geophysical properties of the objects (e.g., surface roughness) [31]. Thus, specific scattering mechanisms caused by complex urban structures such as double bounce can be detected by single polarized PolSAR data [30,32,33].

SAR simulated reflectivity maps have been used for identifying complex backscatter patterns in urban areas [34]. SAR simulation produces a reliable approximation of real SAR data in the azimuth-range plane based on object geometry, surface characteristics, or certain weather conditions. It has been proven that SAR simulators, in conjunction with Light Detection and Ranging (LiDAR) data, can resolve layover and shadow, which are sources of geometric distortion in SAR images [6]. Shadowed areas can be misclassified as flooded due to their low backscatter, while layover results in the misclassification of flooded areas as non-flooded due to the layover strong return. However, they have not been thoroughly investigated in SAR-based flood detection [35].

The purpose behind using reflectivity maps is to help reducing the effect of geometric distortions such as shadows and layovers in the SAR image that can be confused with the backscatter patterns in flooded areas. Therefore, we believe that incorporating simulated reflectivity maps can reduce the false alarms in a predictive model for flood mapping. RaySAR simulator is used along with LiDAR data to extract simulated SAR reflectivity maps including all-reflections, double bounce, shadow, and layover. Considering this, urban flood mapping using SAR images is challenging mainly due to: (1) the existence of complex infra-structure, which causes different backscatter patterns, such as double-bounce and triple-bounce; and (2) the presence of large shadow and layover areas in the SAR images, especially in urban areas with steep terrain, which makes double-bounce detection more challenging and may also affect the flooded area’s classification accuracy.

The contribution of this study is to examine the synergistic use of multi-temporal Polarimetric and Interferometric SAR (PolInSAR) and SAR simulated reflectivity maps for flood detection in urban areas. High-resolution TerraSAR-X images are used to generate multi-temporal PolInSAR features for detecting floods in urban areas, including SAR intensity, InSAR coherence, InSAR phase, and PolSAR Boxcar filter images. Moreover, the RaySAR simulator is used along with LiDAR data to extract SAR simulated reflectivity maps, including all-reflections, double bounce, shadow, and layover. A set of auxiliary features namely elevation, slope, aspect, distance from the river, and land use/land cover (LULC) are used in the flood detection model. These features have been shown to effectively contribute to flood mapping [36,37]. We used a Random Forest (RF) algorithm to distinguish the flooded and non-flooded areas. The proposed method is tested on the 2017 Ottawa River flood in four different areas consisting of urban and sub-urban textures.

One of the main limitations of this study is the number of control points for model training. Based on the available medium-resolution optical and drone images, a limited number of labeled flooded data points (i.e., Ground Control Points (GCPs)) have been produced. Accordingly, an equal number of labeled non-flooded data points were generated to ensure a balanced data set, resulting in limited data points for model training. Moreover, GCP distribution is not uniform which may affect the universality of the results; thus, a data point shuffling is used in this study. Esfandiari et al. (2020) showed that the Random Sample Consensus (RANSAC) method can be used in conjunction with a RF classifier to use pseudo training data points and address the data point limitation [33]. Since this issue was already addressed in the literature, we did not incorporate it in this paper.

1.2. Study Area

This study uses the 2017 Ottawa River flood across Ottawa city, the capital city of Canada, and Gatineau, Quebec, as shown in Figure 1A, as a case study. The Area of Interests (AOIs) include four different areas with different LULC, namely Area1, Area2, Area3, and Area4 as shown in Figure 1B. The 2017 spring flood of Ottawa was primarily caused by exceptionally heavy precipitation and melting snow [38].

1.3. Datasets

This study uses five single polarized SSC StripMap TerraSAR-X images, covering the pre-flood to post-flood periods together with a high-resolution Digital Terrain Model (DTM) for SAR simulation and PolInSAR processing. Historical National Aerial Surveillance Program (NASP) images, captured by a non-metric camera mounted on a drone, are then used to generate the GCPs and evaluate the results. The dataset is shown in

Table 1. Dataset.

Data	Ground Sampling Distance	Date
		Pre-Flood	Co-Flood	Post-Flood
TerraSAR-X (The data was granted based on the proposal approval from the European Space Agency. Accessed on 7 August 2020)	3 m	31 March 2017 22 April 2017 (nearest)	3 May 2017	25 May 2017 (nearest) 27 June 2017
DEM (https://open.canada.ca/data/en/dataset/0fe65119-e96e-4a57-8bfe-9d9245fba06b)	1 m	-	-	2020
NASP UAV images (https://pscanada.maps.arcgis.com/apps/MapSeries/index.html?appid=fd5c6a7e5e5f4fb7909f67e40e781e06)	-	-	-	12 May 2017
River Network (The data was downloaded from https://open.ottawa.ca/)	-	-	-	-
Water level (https://wateroffice.ec.gc.ca/report/historical_e.html?stn=02KF005&dataType=Daily&parameterType=Level&year=2017&mode=Graph)	-	-	-	-
LULC (https://www.openstreetmap.org/#map=2/71.3/-96.8)	-	-	-	-

.

PolInSAR data used in this study include multi-temporal phase, multi-temporal coherence, and multi-temporal Boxcar filter images as complementary information to multi-temporal intensities. For SAR simulation, we utilized all-reflections, double bounce, shadow, and layover features that we anticipate are effective for flood mapping in urban areas; and the five auxiliary features include elevation, slope, aspect, distance from the river, and LULC. To prepare the elevation layer, we used a high resolution 1 m LiDAR DTM and resampled it to 3 m to correspond with the TerraSAR-X images. The resampled elevation layer was then used to produce the slope and aspect features in ArcGIS Pro. The river network map was employed to generate the distance from river layer using the Euclidean distance tool within the ArcGIS Pro.

GCPs were generated by comparing historical NASP images with 2017 historical Google images. A total number of 2450 GCPs were generated and dispersed as follows: 750 points in Area1, 700 points in Area2, 300 points in Area3, and 700 points in Area4. The number of flooded and non-flooded points is even in all areas. Area3 has fewer points compared to other areas due to the limited number of flooded points that are visible from comparing NASP and historical Google images. The flood peak occurred on 7 May 2017, for which no TerraSAR-X images were available. The closest date to the flood peak on which we have a TerraSAR-X image was 3 May 2017. The amount of water level on 3 May 2017 is the closest to that of 7 May 2017. Thus, the closest water level in the TerraSAR-X image captured on 3 May 2017 makes it the best available candidate for the co-flood image. Moreover, NASP images were not available on 3 May 2017, so images taken on 12 May 2017 were chosen since the water level on 12 May 2017 was similar to that of 3 May 2017 as shown in Figure 2.

2. Methodology

The whole flowchart of the work is shown in Figure 3. To simplify the explanation of the methodology, we divide the work into the following sections: InSAR, PolSAR, SAR simulation, and Random Forest Model, as highlighted in the flowchart. We tested different scenarios combining features from multi-temporal InSAR, PolSAR and SAR Simulation along with auxiliary flood mapping features to find the optimum set that can result in maximum accuracy for flood mapping.

2.1. Interferometric SAR

SAR interferometry utilizes the phase difference between two complex SAR images obtained from slightly different positions and/or at different times. Figure 4 shows the geometry of an along-track satellite InSAR system, where $S_p$ and $S_s$ indicate the primary and secondary locations of SAR satellites; H is the altitude of satellite while h is the elevation of an arbitrary point. The distance between the two satellites in the plane perpendicular to the orbit is called the interferometer baseline ($B$); its projection to the slant range is the parallel baseline ($B_{\parallel }$) and its projection perpendicular to the slant range is the perpendicular baseline $\left(B_{\perp }\right)$. The incidence angle (θ) is the angle between the SAR beam and the perpendicular axis to the local topography. Tilt angle (α) refers to the angle between the horizontal plane and interferometer baseline.

SAR interferograms are generated by cross-multiplying the primary SAR image with the complex conjugate of the secondary images [13]. As a result, the interferometric phase is the phase difference between the two images, while the interferometric amplitude is the multiplication of the amplitude of the primary image and the secondary image. Assuming $R_p and R_{s }$ are range distances from the primary image and the secondary image, the primary (${\phi }_p)$, secondary (${\phi }_s$), and interferogram phase $(\Delta \phi ={\phi }_{int})$ are as follows, where $\lambda$ is the wavelength:

$${\phi }_{p }=\frac{4\pi }{\lambda }{R}_p$$

(1)

$${\phi }_s=\frac{4\pi }{\lambda }{R}_S$$

(2)

$$\Delta \phi ={\phi }_{p }-{\phi }_s=\frac{4\pi }{\lambda }\Delta R$$

(3)

The residual interferogram phase then can be written as follows, where ${\phi }_{def}$, deformation phase, is the phase change due to the displacement in line-of-sight direction, ${\phi }_{topo}$ is the phase change due to topography, ${\phi }_{atm}$ is the phase change due to the atmospheric retardation, $\Delta {\phi }_{orb}$ is the residual phase due to orbit errors, and ${\phi }_{noise}$ is the remaining phase noise due to other variables such as thermal noise and coregistration errors [39]:

$$\Delta \phi ={\phi }_{int }={\phi }_{def}+{\phi }_{topo}+{\phi }_{atm}+\Delta {\phi }_{orb}+{\phi }_{noise }$$

(4)

In order to remove the topography phase, interferograms are terrain corrected using a DEM [39]. As a result, Equation (4) can be written as follows, where $\Delta {\phi }_{topo}$ is the residual phase due to DEM errors.

$$\Delta \phi ={\phi }_{int }={\phi }_{def}+{\phi }_{atm}+\Delta {\phi }_{topo}+\Delta {\phi }_{orb}+{\phi }_{noise }$$

(5)

This absolute interferogram phase change is ambiguous and needs to be unwrapped [40]. Unwrapping is the process of removing ambiguity from the wrapped phase. After phase unwrapping, there are still four terms to be corrected in Equation (5) to get the deformation phase. These four terms, i.e., ${\phi }_{atm}, \Delta {\phi }_{topo,} \Delta {\phi }_{orb},$ and ${\phi }_{noise},$ include spatially correlated and spatially uncorrelated error terms. The spatially correlated parts are assumed to be temporally uncorrelated; therefore, they can be estimated and removed by first using high-pass filtering in time and then low-pass filtering in space [41]. This results in the deformation phase based on which InSAR phase features were generated, and spatially uncorrelated error can be modeled as noise [39]. InSAR coherence represents the similarity of radar reflection between two SAR images on a pixel-by-pixel basis and provides a quantitative measure of the interferogram noise. Using N independent images, an estimate of coherence can be calculated as shown in Equation (6), where $u$ represents the amplitude of SAR images, and u* indicates its complex conjugate; this equation was used to generate InSAR coherence ($\widehat{\gamma }$) features used in this study:

$$\widehat{\gamma }=\frac{{{\displaystyle \sum }}_{i=1}^N{u}_{1i}{u}_{2i}^{*}}{\sqrt{{{\displaystyle \sum }}_{i=1}^N{\left|u_{1i}\right|}^2}\sqrt{{{\displaystyle \sum }}_{i=1}^N{\left|u_{2i}\right|}^2}}$$

(6)

2.2. Polarimetric SAR

The SAR Polarimetry concept involves analyzing the polarization state of an electromagnetic (EM) wave. Different SAR systems transmit and receive EM waves in different polarizations, including horizontal and vertical polarization. Different objects generate different combinations of coherent speckle noise and random scattering effects. As a result, for the classification of SAR data, it is important to calculate the average or dominant scattering mechanism. Consequently, speckle filtering is an indispensable initial step. The speckle noise model of a SAR image can be expressed as follows, where $x$ is the pixel without noise, and $s$ represents the speckle:

$$y\left(i,j\right)=x\left(i,j\right) s\left(i,j\right)$$

(7)

Assuming an infinite homogeneous sample, the speckle-free pixels can be extracted from Equation (7) as below, where $\overline{y}$ is the mean of the sample:

$$x\left(i,j\right)=\overline{y}\left(i,j\right)$$

(8)

Although the number of homogenous pixels is not infinite in reality, the optimal scenario in classical speckle filters is to average the N finite pixels. PolSAR boxcar filtering is an effective way that preserves the phase information while reducing speckle [42]. The Boxcar filter reduces the speckle inherent to SAR images through local averaging, similar to Equation (8). The reason Boxcar filtering was chosen over coherent decompositions is that coherent decompositions cannot be calculated with the single polarized TerraSAR-X data in this study, as such decompositions need fully polarized data [31]. Therefore, Equation (8) was used to produce Boxcar PolSAR features used in this study.

2.3. SAR Simulation

There are three steps to SAR simulation: (1) modeling, (2) sampling by a simulator, and (3) scatterer analysis (S. J. Auer 2011). Modeling provides relevant input information about the objects being simulated. After modeling, each detected signal in the scene is sampled using a ray tracer, which involves deriving its 3D position and amplitude, bounce level information, flags to indicate its specular direction, and intersection points to identify the reflecting surface [43,44,45,46]. In the end, a scatterer analysis is performed to provide SAR simulated reflectivity maps representing multiple bounce returns including single bounce, double bounce, and all-reflections.

One of the main problems in SAR simulation, especially in dense urban areas, is how to address the multipath effect due to the summation of signal rays overlaid in one single pixel as the received radar signal. A common solution to resolve multipath is to recover the signal propagation path. For the simulation of radar signal propagation, information such as sensor parameters, the geometry of the objects, and physical scene characteristics need to be specified [45,46]. The scene geometry parameters are used to define the position of the virtual SAR sensor in the modeled scene.

RaySAR simulator which is used in this study is developed for analyzing local urban scenes under the assumption that the local incidence angle of the radar signal is constant [13]. Therefore, a parallel light source and a virtual orthographic camera are considered to represent a radar transmitter and receiver, respectively. Thus, the far-field coordinates of signals can be directly simulated based on the position of the scene center, light source, and virtual orthographic camera as well as the size of the simulated image [34,44]. The scene center corresponds to the center of the 3D object model. The light source and the virtual orthographic camera (i.e., virtual SAR sensor) are therefore the same.

In this study, we used simplified diffuse and specular reflection models to simulate the radiometry of SAR images with a focus on geometrical correctness. Consequently, we followed [31] and developed a processing chain of simulation that uses digital elevation models as input data. The processing chain is illustrated in Figure 5. High-resolution 1 m Digital Terrain Model (DTM) was downloaded from the Natural Resources Canada (NRCan) website. The modified Digital Surface Model (DSM) and normalized Digital Surface Model (nDSM) were generated using the input DTM. The modified DSM represents the geodetic height with UTM horizontal coordinates, with only building height information embedded; height information related to vegetation has been removed from it using a building footprint provided by Open Street Map (OSM) (https://www.openstreetmap.org). The nDSM is derived by subtracting DTM from the modified DSM.

In this context, three sets of 3D object models including DTM, modified DSM, and nDSM were employed to produce reflectivity maps as indicated in Figure 5a. These simulated reflectivity maps were then geocoded using the geoinformation of the DSM, along with the orbit and projection parameters of the real SAR image. Geocoded simulated reflectivity maps include all-reflections and double bounce features derived from the modified DSM and one all-reflections layer from each DTM and nDSM, as shown in Figure 5b. These simulated reflectivity maps are then post-processed as described in [31] to produce four SAR simulation features representing all-reflections, double bounce, shadow, and layover as shown in Figure 5c.

2.4. Random Forest Model

RF is a supervised machine learning model and a powerful tree-based classifier commonly used for general-purpose classification and regression, especially when the number of variables is large relative to the number of observations [47]. Several studies have shown that RF has outperformed other algorithms such as Maximum Likelihood, Artificial Neural Networks, and Support Vector Machines for flood detection and other applications [48,49]. Thus, we used RF to identify flooded areas from non-flooded ones in this study.

The RF classifier was trained and evaluated with different combinations of features, including the SAR simulation and PolInSAR features, along with other auxiliary features that were found to be significant for flood mapping in [36], including elevation, slope, aspect, distance from the river, and LULC. Bayesian Optimization was used to tune RF hyperparameters. The importance of features is computed in the feature selection process based on their weighted mean square error in the node splitting that is based on a random subset of features for each tree [47]. Moreover, to prepare the data points, an equal number of flooded and non-flooded GCPs are distributed across the floodplain to minimize data imbalance.

3. Results

In this section, PolInSAR and SAR simulation analyses are explained. For each combination of features derived from these analyses, a single model is trained using the data points from all four AOIs, i.e., Area1–4. These RF models were then evaluated separately using the dataset sample from each specific area. Consequently, 75% of the data points (2450 ground truth points) were used to train the model over the whole AOIs, and then 25% of data in each specific area was used for validation. The training dataset was randomly selected to avoid autocorrelation.

For PolInSAR analysis, a consecutive set of SSC single polarized TerraSAR-X images spanning the pre-flood, co-flood, and post-flood periods are used, as listed in Table 1, to create both PolSAR and InSAR features. PolSAR and InSAR analyses were conducted separately using Geomatica Banff. For PolSAR analysis, we first radiometrically calibrated the TerraSAR-X images using sigma naught. Sigma naught or backscatter coefficient refers to the normalized return of radar signal from an object, which is measured as per unit area on the ground. Next, the calibrated complex images were converted into detected data. Then, a 5*5 Boxcar filter was applied to reduce speckle noise by local averaging while preserving polarimetric information. Therefore, we produced five locally averaged TerraSAR-X images for each predefined area as PolSAR features, as listed in

Table 2. Generated features from PolInSAR processing over each area.

Analysis (Total Number of Features)	Features	Number of Features
InSAR (8)	InSAR coherence InSAR phase	4 4
PolSAR (5)	Boxcar filter	5
PolInSAR (18)	InSAR coherence InSAR phase Boxcar filter SAR intensities	4 4 5 5

. The PolSAR Boxcar filter images derived from co-flood TerraSAR-X in all four areas are shown in Figure 6.

For InSAR analysis, we used the TerraSAR-X image captured on 3 May 2017 as the primary image, and the other four images as secondary, which resulted in four interferometry pairs. The minimum and maximum temporal baseline of the interferometry pairs used in this study are 12 days and 56 days, respectively. During InSAR processing, the topography-related phase is estimated and removed from the interferograms using a DEM. InSAR coherence and InSAR phase associated with each interferogram are then extracted. As a result, we had four InSAR coherence and four InSAR phase features for each study area as listed in Table 2. Figure 7 illustrates sample InSAR coherence images over four areas.

For SAR simulation analysis, we used the open-source ray tracer RaySAR which is a modified version of POV-Ray implemented in MATLAB [50]. We used a TerraSAR-X StripMap SSC image captured on 3 May 2017 (flood peak), presented in a geographic coordinate system with ellipsoidal correction, as our real SAR image. Ellipsoidal correction is a process in which the SAR image is being resampled and projected from radar coordinate to map coordinate [44,50]. Next, the following parameters were derived from the metadata of the TerraSAR-X image to define the geometry of the scene being simulated: (1) the azimuth angle; (2) the frame mean height; (3) incidence angles of the scene center and its four corners: used for interpolating the incidence angle at the scene center, assumed to be locally constant over the entire scene; (4) East and North pixel spacing which is equal to the azimuth and range pixel spacing in SAR coordinate system.

As the second input data, three sets of 3D object models including DTMs, modified DSMs, and nDSMs based on LiDAR data were used. Using DTMs, modified DSMs that only represent terrain and constructed objects, not vegetation, were generated by using a DTM, buildings footprints, and buildings elevation. nDSM is derived then by subtracting DTM from the modified DSM. The sampling of these object models should be the same as SAR image sampling. Therefore, DEM object models were resampled to 3 m using a bilinear interpolation method. We used the same coordinate system for both data sources, i.e., WGS 84 of UTM, and ellipsoidal heights. To perform ray tracing, object models were converted to the POV-Ray format using AccuTrans. Consequently, the discrete position of the simulated signals, as well as their strength and reflection levels are generated as explained in [43]. These features are then used to create the SAR simulated reflectivity maps. These simulated images that are in the Radar coordinate system, i.e., azimuth-range geometry, are then geocoded to create SAR simulation features as listed in

Table 3. SAR simulation outputs.

Input Object Models	SAR Simulated Reflectivity Maps from RaySAR	Geocoded Simulation Features Based on Reflectivity Maps
DTM	All-reflections	All-reflections, Double Bounce, Layover, Shadow
Modified DSM	All-reflections + Double Bounce
nDSM	All-reflections

. Details of the geocoding process can be found in [34].

The geocoded simulation outputs are post-processed to generate four features, namely all-reflections, double bounce, layover, and shadow, as described in [34] for each AOIs. Figure 8 illustrates these four features generated over AOIs. Moreover, the generated flood maps over Area1–4 are shown in Figure 9.

For each AOI, 1731 scenarios were examined as listed in

Table 4. Features used in different scenarios of A Categories.

A Categories (Number of Features)	Feature Types	Name of Features	All Possible K-Combinations
A1- SAR Intensity (5)	SAR intensities	Pre-, Nearest pre-, co-, nearest post-, and post-flood intensities	${\displaystyle {\displaystyle \sum }_{k=1}^5}\left(\begin{array}{c}5\\ k\end{array}\right)=31$
A2- InSAR (8)	InSAR phases, InSAR coherences	Pre-, nearest pre-, nearest post-, and post-flood InSAR coherences/phases	${\displaystyle {\displaystyle \sum }_{k=1}^8}\left(\begin{array}{c}8\\ k\end{array}\right)=255$
A3- PolSAR (5)	Boxcar filtered images	Pre-, Nearest pre-, co-, nearest post-, and post-flood Boxcar filtering	${\displaystyle {\displaystyle \sum }_{k=1}^5}\left(\begin{array}{c}5\\ k\end{array}\right)=31$
A4- Simulation (4)	Reflectivity maps	All-reflections, double bounce, layover, shadow	${\displaystyle {\displaystyle \sum }_{k=1}^4}\left(\begin{array}{c}4\\ k\end{array}\right)=15$
A5- Baseline (5)	Auxiliary features	Elevation, slope, aspect, distance from the river, LULC	${\displaystyle {\displaystyle \sum }_{k=1}^5}\left(\begin{array}{c}5\\ k\end{array}\right)=31$
			Total of 363 scenarios

and

Table 5. Features used in different scenarios of B Categories.

B Categories (Number of Features (k))	Included Selected Features from A Category	Total Examined Combinations = (2k − 1) × k
B1- PolInSAR (9)	A1 + A2 + A3	(2 × 9 − 1 = 17) × 9 =153
B2- All SAR (12)	A1 + A2 + A3 + A4	(2 × 12 − 1 = 23) × 12 = 276
B3- Without Simulation (14)	A1 + A2 + A3 + A5	(2 × 14 − 1 = 27) × 14 = 378
B4- All Categories (17)	A1 + A2 + A3 + A4 + A5	(2 × 17 − 1 = 33) × 17 = 561
		Total of 1368 scenarios

and each scenario had a different number of features. To identify the most important features for flood mapping, different combinations of features were tested as different scenarios in a two-step procedure: A category scenarios, and B category scenarios. In A category, features are divided into: (A1) SAR intensity, (A2) InSAR, (A3) PolSAR, (A4) SAR simulation, and (A5) five auxiliary features as the baseline, shown in Table 4. In B categories, selected features from A category were integrated to build up the B category: (B1) PolInSAR, (B2) All SAR, (B3) without SAR simulation, and (B4) All categories. Table 5 shows included features in each B category.

For A categories, we explored a total of 363 scenarios covering all possible k-combinations of features in each category. First, an RF algorithm was trained using all possible k-combinations of features within each category A1–5, as shown in Table 4. Then, the feature importance was estimated by permuting out-of-bag observations among the trees; Consequently, the resulting best features contributing to the top three highest overall accuracies within each category A1–5 were marked as “selected features”, shown in

Table 6. Selected features from A Categories.

A Categories	Selected Features (Number of Features)
A1- SAR Intensity	Co-flood and nearest pre-flood and post-flood intensity (3)
A2- InSAR	All four InSAR coherences, nearest post-flood InSAR phase (5)
A3- PolSAR	Co-flood PolSAR Boxcar filtered image (1)
A4- Simulation	Double bounce, Shadow, All-reflections (3)
A5- Auxiliary	Elevation, Slope, Aspect, Distance from the River, LULC (5)

. These “selected features” in A1–5 categories are then combined to be tested as different B1–4 categories, as listed in Table 5.

For B categories, we combined the selected features from A categories that led to the highest overall accuracies. Since including correlated features can affect the RF model and reduce the overall accuracy, we started training RF with a single feature within each B1–4 category and then added the rest of the features one by one until we used all features within that certain category. The features were then removed one by one until the very first feature we started with are left. Then, the feature importance was estimated by permuting out-of-bag observations among the trees, similar to what we did in the A category. As a result, features are prioritized based on their importance according to the RF algorithm. The most important features associated with the top three highest accuracies achieved in B categories are integrated and listed in

Table 7. Selected features in B category.

B Categories	Selected Features in B Categories
	Intensities	InSAR	PolSAR	Simulation	Auxiliary
B1- PolInSAR	Co-, nearest pre-, and post-flood intensities	All four InSAR coherences, nearest post-flood InSAR phase	Co-flood Boxcar filtered image	-	-
B2- All SAR	Co-, nearest pre-, and post-flood intensities	All four InSAR coherences, nearest post-flood InSAR phase	Co-flood Boxcar filtered image	Shadow, Double Bounce, All-reflections	-
B3- All without Simulation	Co-, nearest pre-, and post-flood intensities	All four InSAR coherences, nearest post-flood InSAR phase	Co-flood Boxcar filtered image,	-	All Auxiliary Features
B4- All Categories	Co-, nearest pre-, and post-flood intensities	All four InSAR coherences, nearest post-flood InSAR phase	Co-flood Boxcar filtered image,	Shadow, Double bounce, All-reflections	All Auxiliary Features

as selected features in B categories.

4. Discussion

As stated in the previous section, feature selection was performed in a two-step procedure where we tested different scenarios in two categories (i.e., A category and B category) to find the most effective features. Since RF is a binary classifier, the accuracy is evaluated by comparing classification overall accuracy. Consequently, the highest classification overall accuracies achieved in the A category and B category are reported in

Table 8. Acquired classification accuracies of A Categories: “All” columns indicate the accuracies where all features were included while “selected features” columns show the best accuracies achieved in each category using the selected features.

A Categories	Classification Overall Accuracies
	Area1		Area2		Area3		Area4
	All Features	Selected Features	All Features	Selected Features	All Features	Selected Features	All Features	Selected Features
A1- Intensity	84.8	85.3	85.85	86.1	60.3	61	86.1	86.5
A2- InSAR	82.2	82.5	82.7	82.9	60.3	60.6	83.6	83.7
A3- PolSAR	86.9	86.9	86.2	86.2	63.2	63.2	87.1	87.1
A4- Simulation	53.4	53.4	53.8	53.8	52.4	52.4	54.5	54.5
A5- Baseline (auxiliary)	83.1	83.6	83.5	83.6	82.5	82.5	82.7	82.8

and

Table 9. Highest classification accuracies achieved in B category.

B Categories	Classification Overall Accuracies
	Area1	Area2	Area3	Area4
B1- PolInSAR	88.4	89.1	61.7	90.1
B2- All SAR	88.6	89.3	61.7	90.3
B3- All without Simulation	91.3	92.2	82.1	92.8
B4- All Features	92.6	93.5	83.8	93.2

. Overall accuracies are computed by dividing the number of data points that were correctly labeled as flooded and non-flooded points to the number of all labeled data points as follows:

$$Overall Accuracy=\frac{TP+TN}{TP+TN+FP+FN}$$

(9)

Comparing B category scenarios over different AOIs, as listed in Table 9, shows that the results of Area1, Area2, and Area4 are exhibiting a similar trend, while Area3 has a lower accuracy than the other areas, which is due to the lower data points available in Area3, and its different and more vegetated LULC compared to other areas. Overall, among various features, the most important ones for flood mapping, which produced the highest accuracies in each area, were derived from the incorporation of all A categories, i.e., the integration of SAR intensities, InSAR, PolSAR, SAR simulation, and auxiliary features, as listed in B4 category in Table 9. The confusion matrix associated with the highest accuracy achieved in A categories and B categories over all four AOIs is listed in

Table 10. The performance measurements True Positive (TP), True Negative (TN), False Positive (FP), and False Negative (FN) associated with the highest accuracies achieved using selected features in A category for all four AOIs as highlighted in Table 8.

Category	TP (%)	TN	FP	FN	TP	TN	FP	FN
	Area1				Area2
A1	92.4	78.2	21.8	7.6	92.8	79.4	20.6	7.2
A2	93.4	71.6	28.4	6.6	93.5	72.3	27.7	6.5
A3	97.8	76	24	2.2	99.4	73	27	0.6
A4	9.8	97	3	90.2	9.8	97.8	2.2	90.2
A5	92.8	74.4	25.6	7.2	87.9	79.3	20.7	12.1
Category	Area3				Area4
A1	43.1	78.9	21.1	56.9	92.2	80.8	19.2	7.8
A2	61.6	59.6	40.4	38.4	93.9	73.5	26.5	6.1
A3	75.3	51.1	48.9	24.7	89.9	84.3	15.7	10.1
A4	4.8	100	0	95.2	10.8	98.2	1.8	89.2
A5	91.6	73.4	26.6	8.4	87.7	77.9	22.1	12.3

and

Table 11. The associated performance measurements True Positive (TP), True Negative (TN), False Positive (FP), and False Negative (FN) associated with the highest accuracies achieved using selected features in B category for all four AOIs as listed in Table 9.

Category	TP (%)	TN	FP	FN	TP	TN	FP	FN
	Area1				Area2
B1	95.4	81.4	18.6	4.6	96.9	81.3	18.7	3.1
B2	95	82.2	17.8	5	97.1	81.5	18.5	2.9
B3	93.9	88.7	11.3	6.1	92.9	91.5	8.5	7.1
B4	92.9	92.3	7.7	7.1	95.4	91.6	8.4	4.6
Category	Area3				Area4
B1	64.8	58.5	41.5	35.2	97.1	83.1	16.9	2.9
B2	63.8	59.6	40.4	36.2	97.1	83.5	16.5	2.9
B3	92.6	71.6	28.4	7.4	94.7	90.9	9.1	5.3
B4	87.9	79.7	20.3	12.1	93.3	93.1	6.9	6.7

. Investigating all scenarios in different categories reveals that the False Positive (FP) rate stands above the False Negative (FN) as shown in Figure 10, which means the proposed model overestimated the flood areas. Although the best practice is to remove all false results, we believe that overestimating flooded areas is better than underestimating as a high FN rate can result in missing damaged areas in post-disaster recovery efforts.

In the A5 category, we used five auxiliary features, i.e., elevation, slope, aspect, distance from the river, and LULC, that are reported in the literature to be effective in flood mapping as our baseline. A5 category reported accuracies ranging from 82.5% to 83.6% in different AOIs. Comparing A5 to B1, where PolInSAR features were involved, showed that using PolInSAR features improved the model by 4.8%, 5.5%, and 7.3% in Area1, Area2, and Area4. Comparing B4 to B3, where PolInSAR features are considered together with five auxiliary features, i.e., all features except for simulation ones, revealed that adding SAR simulated features improved accuracy by 1.3%, 1.3%, and 0.4% in Area1, Area2, and Area4.

Table 12. Average classification overall accuracy improvement over Area1, Area2, and Area4. This table compares the proposed models accuracy to the existing models accuracy: e.g., it shows that the proposed B4 model improved flood detection by an average of 7.1%, 10%, and 6.3% compared to the SAR based models used in existing studies, and by an average of 9.7% compared to the non-SAR based model, i.e., by using auxiliary features.

Proposed\Existing	SAR-Based Methods			Non-SAR-Based Method
	A1: Intensity	A2: InSAR	A3: PolSAR	A5: Auxiliary
B4: Incorporated PolInSAR, SAR Simulation, and Auxiliary Features	7.1%	10%	6.3%	9.7%
B2: PolInSAR and SAR Simulation	3.4%	6.5%	2.7%	6.1%

provides a comparison of classification overall accuracy improvement of the proposed models and four existing models commonly used for flood detection in rural and urban areas. Comparing A5 to B2, where all SAR features are used, showed that synergistic use of simulated images and PolInSAR features without including auxiliary features outperformed the baseline by increasing the accuracy by 5%, 5.7%, and 7.5% in Area1, Area2, and Area4, respectively; the average improvement over those areas are shown in Table 12 as 6.1%. When we used PolInSAR and simulation features in combination with the auxiliary features as in the B4 category, we reached even higher accuracy with an increase of 9%, 9.9%, and 10% compared to the baseline in Area1, Area2, and Area4, respectively; the average improvement over those areas are shown in Table 12 as 9.7%. Additionally, comparing B4 to A1, where we used multi-temporal SAR intensities alone, revealed that the proposed method achieved an accuracy improvement of 7.3%, 7.4%, and 6.7%, respectively in Area1, Area2, and Area4, which is a common practice for flood mapping in the literature; the average improvement over those areas are shown in Table 12 as 7.1%.

Exploring all scenarios associated with simulation features, i.e., A4, B2, and B4 categories, we noticed that including shadow features always reduced false negatives and helped the classifier to improve flood detection more effectively. In this regard, the B4 category showed a higher accuracy compared to B3, which demonstrates that simulation features, especially shadow, layover, and double bounce features that were the most important, improved flood detection by 1.3%, 1.3%, and 0.4% in Area1, Area2, and Area4, respectively; comparing B3 to B4 shows that the synergistic use of SAR simulation and PolInSAR outperformed employing PolInSAR features alone by an average of 1% improvement.

The results of this study are consistent with the PolInSAR outcome of our previous study [28] and support the assumption that investigating the trend in the time series of PolInSAR features concerning their LULC class can improve identifying the complex urban backscatter patterns.

5. Conclusions

In this work, we investigated the effect of employing InSAR, PolSAR, and SAR-simulated reflectivity maps in the improvement of urban flood detection using a RF classifier. Examining the synergistic use of multi-temporal PolInSAR features and SAR-simulated reflectivity maps for urban flood detection is the main contribution of this study. These SAR features are employed along with a set of auxiliary features that have been reported in the literature to contribute to flood mapping, including elevation, slope, aspect, distance from the river, and LULC. Multi-temporal SAR intensity, InSAR coherence, and InSAR phase features were generated using InSAR analysis. Additionally, PolSAR analysis was employed to generate PolSAR Boxcar filter features. Simulated reflectivity maps features used in this study included all-reflections maps, double-bounces, shadows, and layovers. Different combinations of these features were examined to select the most effective features in flood detection.

In general, the results showed that employing PolInSAR and SAR simulated features along with auxiliary features in an RF model can improve flood detection in urban areas. The highest accuracy was achieved from a subset of features derived from the performed SAR analyses, i.e., InSAR, PolSAR, and SAR simulation. The results showed that the incorporation of single polarized high-resolution TerraSAR-X images with auxiliary features outperforms the models that merely use the auxiliary features.

The 2017 Ottawa River flood was explored using a set of high-resolution TerraSAR-X images. Four study areas with different LULCs were examined, for which GCPs were sampled by comparing historical NASP images with 2017 historical Google images. We employed multi-temporal SAR intensities, multi-temporal PolInSAR features including InSAR phase, InSAR coherence, and PolSAR Boxcar filter, as well as SAR simulated features including all-reflections, double bounce, shadow, and layover.

Feature selection was performed in a two-step procedure where we tested different scenarios in two categories to find the most effective features. Among all investigated features, the best accuracy of 93.1% was achieved by considering Co-flood and nearest pre-flood and post-flood intensity, all four InSAR coherence, nearest post-flood InSAR phase, co-flood PolSAR Boxcar filtered image, double bounce, shadow, all-reflections, elevation, slope, distance from the river, and LULC. The results revealed that the synergistic use of SAR simulated reflectivity maps and PolInSAR features together with auxiliary features can improve urban flood detection by an average of 9.6% when compared to the baseline. Moreover, without considering any auxiliary features, the exclusive use of PolInSAR and SAR simulated features outperformed the baseline by 7%, which indicates the effectiveness of using SAR features in flood mapping. This study indicates that SAR simulated features are effective in decreasing false positives, which helped the proposed method exceed flood detection accuracy by merely use of PolInSAR features. Furthermore, the proposed method outperformed the exclusive use of SAR intensity by an average of 7.1% accuracy improvement, while SAR intensity is the basis for most existing SAR-based flood mapping methods. Given the promising results of using TerraSAR-X in flood detection, our future work with TerraSAR-X will focus on using fully polarized images as well as on examining other study areas with different terrains and LULC.

It needs to be acknowledged that the proposed model uses a greater number of input features compared to some other flood mapping techniques; however, this is compensated for the improved accuracy. Since accurate flood mapping can contribute to saving lives, reducing damage, and providing accurate damage estimates, the promisingly improved proposed flood mapping method is a good candidate for flood detection in urban areas. Especially knowing that the proposed model is independent of optical imagery and fully operates based on SAR and some static features (e.g., elevation models and slope), the proposed flood mapping method can be used in all weather conditions and at any time of the day which speeds of the flood mapping and consequently reduces the response time. However, it should be noted that although the proposed model worked well in urban and semi-urban areas (i.e., Area1, Area2, and Area4 of the second study area over Ottawa), it did not work as well in a densely vegetated area (i.e., Area3 of the second study area); thus, we recommend that the proposed model be extrapolated to areas with similar LULC.