SUPERVISED CLASSIFICATION METHODS FOR AUTOMATIC DAMAGE DETECTION CAUSED BY HEAVY RAINFALL USING MULTITEMPORAL HIGH RESOLUTION OPTICAL IMAGERY AND AUXILIARY DATA

In the context of climate change and rising frequency of extreme hydro-meteorological events around the world, flood risk management and mapping of heavy rainfall-related damages represent an ongoing critical challenge. For decades now, remote sensing has been largely used to investigate spatial and temporal changes in land use and water resources. Today, different satellite products provide fast and crucial knowledge for the study of hydrological disasters over large areas, possibly in remote regions, with high spatial resolution and high revisit frequency. Yet, until now, few works have sought to detect the full range of extreme rainfall-related damages with optical imagery, especially those caused by intense rainwater runoff beyond the direct vicinity of major waterways. The work presented in this paper focuses on the Aude severe weather event of October 15th, 2018, in the South of France, for which more than a thousand claims for agricultural disaster were registered, both related to river overflowing and rainwater runoff. The full resources of ground truths, contextual information, land use as well as digital elevation model (DEM) combined to high resolution and high frequency optical imagery (Sentinel-2, Pléiades) are used to develop an automatic damage detection method based on supervised classification algorithms. Through the combination of several indicators characterizing heterogeneous spectral variations among agricultural plots following the event, a Gaussian process classifier achieved various classification accuracies up to 90% on a large comparable and independent photo-interpreted validation sample. This work builds great expectations for applications in other areas with contrasted climate, topography and land cover.


INTRODUCTION
Floods arise in the aftermath of extreme (either in magnitude, rate or duration) rainfall events. Depending on the territories, they appear and persist over different time scales from hours to weeks, mainly triggered by surface and subsurface runoff. During such events, overland flow occurs under various configurations such as in high damp slopes, lowlands and floodplains, small channels or dry thalwegs (Cerdà et al., 2021). Rainwater-related damages can thus take many forms and emerge in multiple places. Because of their intensity and often dramatic social and economic consequences, the vast majority of research and operational activities on flood damage detection mainly focuses on overflowing of major waterways (e.g. Copernicus EMS Rapid Mapping activity in particular through SERTIT). In this context, assessing flood extent from satellite images has been a pressing topic in remote sensing of natural disasters for decades (Rahman, Di, 2017). Many studies have established methods to identify flooded and flood-prone areas from satellite imagery through direct identification of water bodies, either with optical instruments like Landsat (Yamagata, Akiyama, 1988;Swain et al., 2020) and Sentinel-2 (Pulvirenti et al., 2020;Goffi et al., 2020), or with SAR products (Hostache et al., 2007;Matgen et al., 2011, Rambour et al., 2020). Yet, very few studies have undertaken the detection of the full range of heavy rainfallrelated damages, especially those caused by intense rainwater runoff. These disturbances often take place during short, hardly observable time periods, potentially anywhere with little topography and outside the direct vicinity of major waterways. * Corresponding author Some authors have developed change detection techniques to detect areas affected by flood and erosion from a heavy rainfall, but over long time intervals and with medium spatial resolution (30 m from Landsat, Dhakal et al., 2002), or with higher resolution but strictly limited to overflowing (4 m from Kompsat-2, Byun et al., 2015). In addition, these works have no real aim at replicability (only one study area and over a single event). Numerous related topics and methodologies mostly relying on multitemporal analysis have been investigated for years in the remote sensing community, including the study of gullies (Fadul et al., 1999) and soil erosion (Dwivedi et al., 1997;Begueria, 2006;Sepuru, Dube, 2018), landslides (Danneels and Havenith, 2007;Mwaniki et al., 2015), or agricultural losses and crop yields (Pantaleoni et al., 2007). Beyond the traditional pixel-oriented methods, object-oriented change detection approaches have recently become more popular to detect land cover changes (Huang et al., 2018). In France, around 85% of natural disaster claims between 1982 and 2010 were flood-related. However, only half of them were likely to be associated with river overflowing (Breil et al., 2016). Cerbelaud et al. (2020) laid out the potential of a combined use of multitemporal multispectral high resolution images (Sentinel-2), very high resolution (VHR) optical imagery (Pléiades) and contextual information (plot-based approach) to detect different types of damages on a large area, including far from waterways, resulting from extreme weather events. Using two Sentinel-2 images acquired 10 days before and after the storm and focusing on a small-sized sample of flooded agricultural lands, they found that most damaged plots present distinctive intra-plot variability for the relative difference of specific spectral indices in comparison to unaffected areas. Because of its transient nature and potential to appear anywhere in a given territory, intense rainwater runoff is very difficult to observe and detect per se. Its various consequences in the form of erosion, mudslides, deposits and uprooting are on the other hand usually sizeable and durable, although damages caused to human infrastructures are generally quickly cleaned up. The main assumption at the heart of this study is that these scars (deterioration of vegetation, exposed bare soil, sediment deposit) can be traced back to specific signatures through spectral variations in time and in space. However, detecting these signatures requires to be able to control for the unique and distinctive reactions of each land cover to intense precipitations and seasonal changes. The originality of this work is thus built upon (i) a contextual plot-based approach so as to work with piecewise constant land cover areas and (ii) a change detection method relying on the combined search for spectral (use of indices involving several bands), temporal (before and after the event) and spatial (intraplot heterogeneities) variations. A supervised classification method for automatic damage detection caused by heavy rainfall is thus developed so as to efficiently identify recently flooded areas, not only nearby but also far away from waterways, where intense rainwater runoff can be accountable for. The strength of the methodology also resides in deriving classification capability from a large sample of certified ground truths, change images only twenty days apart and with a 10 m spatial resolution, along with validation data based on very high resolution post-event images from Pléiades satellites. After presenting the materials and detailing the approach, results are interpreted and discussed in the context of future studies on the robustness and replicability of the methodology to different times of the year and contrasted regions in terms of climate, land use and topography.

Study area, territorial subdivision and land use
This study focuses on the Aude flooding event of October 15 th , 2018, in the south-west part of France. Within a 24h-period, close to 300 mm of rainfall was measured around the city of Carcassonne, and up to 250 mm fell locally in only 6 hours (Lebouc et al., 2019), leading to both flooding by river overflowing and intense rainwater runoff. For this work, the region of interest extends over roughly 1 150 km 2 in the Fresquel, Orbiel and Aude watersheds ( Figure 1). As this study's design relies on operating at a plot scale, the department official land cadastre dataset was retrieved under QGIS 3.6 to obtain the most natural territorial subdivision with piecewise constant land use. The area of interest thus comprised more than 200 000 plots under the 2018 version of the land cadastre, with an average area of 0.5 ha and a high standard deviation of 1.7 ha due to the disparities between urban and rural territories. The OSO French land cover product available at a 10 m spatial resolution was used to determine the main land use (LU) category within each plot. According to the OSO 2018 raster, most of the study zone is covered up with grasslands (30% of pixels), vineyards (24%), forests (20%), built-up (14%) and other types of cultures (cereals, protein crops, sunflower etc., 12%).

Multispectral satellite optical imagery products
Two complementary types of satellite products were acquired for this work (Table 1). First, the earliest post-event very high resolution (VHR) Pléiades image (ortho product) from November 3 rd was used to confirm and/or to identify traces of overflowing and rainwater runoff on agricultural lands, grasslands, roads and other diverse works. With a 0.5 m multispectral resampled spatial resolution, this product was a valuable tool to differentiate damaged from non-affected areas by photo-interpretation. Different types of damages were thus identifiable: landslides and mudslides, gully erosion, sediment deposit and vegetation uprooting. Secondly, multiple images from one Sentinel-2 tile (T31TDH from Sentinel-2A satellite, level 2A treatment, Flat Reflectance data, corrected for atmospheric effects and with cloud and shadow masks) were acquired at different dates and crosscorrected for spatial offsetting. Along with a medium to high spatial resolution of 10 m for visible and near infrared bands, Sentinel-2 (S2) data was the key asset to this work through its great revisit frequency (around 5 days) allowing to closely monitor spectral variations before and after the event (Table 1) Table 1. Imaging capabilities of Pléiades and Sentinel-2 instruments in the visible and near infrared bands (NIR).

Selection of ground truths and photo-interpreted plots
Following this hydrological catastrophe, around 900 claims for agricultural disaster were registered and certified by the local authorities (DDTM 11). From these declarations, more than 1 000 damaged lands were geo-referenced and the corresponding plots were tagged on the department land cadastre. After visual inspection of the complete sample thanks to the Pléiades image, 310 plots with clearly identifiable damages were selected to form the ground truths training sample (hereafter referenced as DGT class). In addition, 362 damaged but yet unclaimed plots were hand-picked independently by photo-interpretation in the surrounding areas for validation purposes (hereafter referenced as DPI class). Lastly, 480 undamaged lands were chosen closeby as well as further away for the control group (hereafter referenced as UDPI class). The resulting tagged plots can be seen in Figure 2. The visual approval of only 310 out of 1 000 claimed damaged plots for this study can be explained through multiple rationale. First, without a very high resolution prior image close enough in time, a fair amount of damages could not be unquestionably certified at the given spatial resolution of 0.5 m.
Second, when the damage was mainly related to ponded water and did not involve modifications in the structure of the field (e.g. vineyard), the absence of noticeable marks on the ground after water had receded between the time of the event and acquisition of the image made it hard to confirm the damage. Moreover, during this same time interval (3 weeks), small damages were likely to have been repaired by land owners. Finally, various areas either contained or were surrounded with built-up, high trees and shadows at their edge that made it hard to discriminate the spectral signature of the environment from a possible damage. In order to make sure that the resulting dataset is representative of the study area and that damaged classes are somewhat comparable, Table 2 presents the main characteristics of each class. In addition, Figure 3 shows the distribution of plots' major land use, overall and by class, excluding the built-up category. Logically, vineyards are over-represented within the DGT class because of their high economic value, with land owners usually awaiting consequent financial compensations for harvesting losses.

Selection of spectral indicators, auxiliary data and supervised classification algorithms for damage detection
Following the work of Cerbelaud et al. (2020), this study focuses on information derived from soil and vegetation-specific spectral indices based on the visible and near infrared 10 m bands of Sentinel-2 images: NDVI for Normalized Difference Vegetation Index, NDWI for Normalized Difference Water Index (Gao, 1996), BI for Brightness Index and SAVI for Soil Adjusted Vegetation Index (Qiu et al., 2017; see Table 2 from Cerbelaud et al., 2020 for details). Two formulations of the BI were tested here, using either the red and NIR bands, or including the green band as well. Physics-wise, this study's core assumption is that river overflowing and intense rainwater runoff can be traced back in the form of erosion, mudslides, various deposits and uprooting, to specific signatures (presence of vegetation, sediments, barren soils) through spectral variations in time and space. Therefore, new images (called change images) were produced for each spectral index k by computing the pixel-by-pixel relative difference (RD) between the closest cloud-free Sentinel-2 images before and after the event (October 5 th and 25 th ), written . Particular attention was thus paid to the correction applied to the Sentinel-2 tiles for spatial offsetting (see Figure 4 for the final NDVI change image). In order to work at a plot scale, intra-plot pixel statistics θ -mean, maximum, minimum and standard deviation -were estimated on each change image to produce plot-specific vectors of characteristics ( ) ( Figure 5). The underlying goal was to highlight singular spectral variations within damaged plots with high values of , implying a decrease in the index and thus a deterioration in vegetation cover, as well as the presence of spatial heterogeneities. Because some land covers react in distinctive ways to intense precipitations, having the most natural territorial subdivision with piecewise constant land use was crucial in this approach. These vectors were then completed with other characteristics such as mean and maximum slope (in °) and major land use category. Summary statistics were first performed on the damaged ground truths DGT class and the undamaged photo-interpreted UDPI class to determine the potential of each ( ) variable to discriminate the two classes. In order to make the best of and determine the optimal combination of variables to detect damages, a Pearson correlation matrix was then computed on the same dataset. Indeed, the association of non collinear indicators is likely to lead to the highest explanatory power and finest classification capacity. In order to develop an automatic damage detection method, two samples, one for training and the other for validation, were thus put together. First, the UDPI class of 480 plots was split in half with a random seed under Python 3.7 both for training and validation (later referenced as UDPI*) purposes. The ground truths DGT class then completed the learning sample (550 plots total, 11.3 km 2 ) while the damaged photo-interpreted DPI class filled in the validation sample (602 plots, 14.6 km 2 ). Lastly, plot-wise supervised classifications were achieved based on the combination of LU, mean slope as well as ( ) with three different classification methods: (i) the k-nearest neighbors algorithm (k-NN); (ii) the Multi-Layer Perceptron neural network (two hidden layers and 100 neurons per layer); (iii) the Gaussian process classifier (squared-exponential kernel).

RESULTS
Descriptive statistics (Figure 6) performed on the DGT and the UDPI classes, involving almost 800 plots and covering a 17.1 km 2 area, confirm while slightly revise the findings of Cerbelaud et al. (2020). First of all, substantial overlaps between the two classes for all statistics of both formulations of the brightness index (BI 1 using the red and NIR bands, BI 2 including the green band as well) indicate their lack of relevance to identify traces of intense rainwater runoff. On the contrary, NDVI, NDWI and to a lesser extent SAVI seem to be suitable candidates for efficient detection. For all indices, the large overlap observed for the minimum statistics between undamaged and affected plots reveals the metric's inadequacy. On the other hand, the mean, the maximum and the standard deviation turn out to be promising candidates for discrimination and confirm the physical explanation and methodology chosen in this study. Indeed, damaged plots show greater spatial heterogeneity and higher average values of the relative difference for every index, indicating noteworthy reactions inside affected and likely deteriorated lands.  Table 3. Pearson correlation coefficients between ( ) over DGT class (310 plots) and UDPI class (480).
The Pearson correlation matrix between the remaining promising indicators was computed on the same dataset in Table 3. All coefficients are statistically significant at the 0.1% significance level (from Student table with degree of freedom around 800). Greener boxes point out lower correlation coefficients (chosen arbitrarily below 0.3) between indicators, implying that rather independent information can be extracted from them. First of all, each type of statistics is strongly correlated within all spectral indices. Combination of distinct types is thus considered. Because of larger overlap between classes for ( ) ( Figure 6) and greater correlation coefficients with other variables for ( ) , ( ) seems to be the most appropriate candidate in combination with the mean of any of the three indices. Maximum statistics appear to be self-explanatory for all indices. Figure 7 exhibits 2D scatter plots of the two classes in order to visualize classification potential of some indicators. Good candidates would induce almost separate data clouds for DGT and UDPI classes. Several sets of variables were thus tested. Comparable results were obtained with the three types of classifiers (not shown here). Therefore, the Gaussian process classifier was retained for its simplicity and because of the physical meaning it conveys. Beyond binary classification, probabilities can be derived from it to estimate confidence in each class (Figures 8 and 9). Table 4 shows results from 2-variable classifications on the validation sample overall and by class. With potentially one pixel to induce detection, using indicators involving the intra-plot maximum of the RD pixels of NDVI or NDWI (lines 3 and 4) unsurprisingly leads to overestimation of damaged lands with smaller producer accuracy on UDPI* plots. In order to avoid false positives as much as possible, more stable scores across the two classes are obtained with the combined use of standard deviation and mean statistics (lines 1 and 2). As expected, the use of ( ) brings lower results (lines 5 and 6). Finally, inclusion of land use in the classification process proved beneficial with accuracies risen by around 2% and less false positives with better predictions of the UDPI* class (Table 5). The most satisfying classification obtained with a Gaussian process classifier, scoring 91.0 % accuracy, was thus with three variables: ( ) , ( ) and LU. Processing of the entire 1 150 km 2 area then allowed to produce maps displaying probability of damage classification given by the Gaussian process classifier. Different close-ups can be seen in Figures 8 and 9. Both overflowing and intense rainwater runoff are thus easily identifiable in flagged plots (Figure 8).

DISCUSSION
This work mainly relied on the use of the closest cloud-free preevent and post-event Sentinel-2 images around October 15 th , 2018. In order to take advantage of the durable effects of intense rainwater runoff on soils' spectral signatures, multiple post event images could be involved in the classification process. Similarly, because they do not make up for evidence of intense water runoff,    the various impacts of seasonal changes and local anthropogenic activities on spectral indexes between the pre and post event images need to be accounted for through several sensitivity analyses. To quantify this unwanted effect, the change detection approach could be tested using two images (more or less distant in time) between which no particular intense rainfall activity occurred. A plot-wise demeaning could be applied to all change images at the pixel level, allowing to get rid of the average evolution of all pixels associated with a given land use type. Once trained on the Aude learning sample from October 2018, applying the classification method to other extreme weather events in the same area at different times of the year (e.g. spring, summer) will also be enlightening. Further tests on other regions with contrasted climate, topography and land cover will also provide key elements as to the replicability of the approach. Inclusion of Sentinel-2 SWIR bands (B11 and B12) at 20 m resolution, which would allow to derive NDMI (moisture) change images, could also improve detection of rainwater-related damages, albeit downgrading spatial resolution. Several other improvements will also be investigated using hydrological tools involving topographical considerations as well as patch-based classification methods with post event VHR Pléiades images. Since it is based on optical imagery, it has to be acknowledged that this method would be arduously applicable to cloud-prone regions (e.g. small mountainous islands, tropics). Eventually, the output probability maps will further be used to evaluate and validate conceptual or physics-based intense rainwater runoff mapping models such as IRIP © (Indicateur de Ruissellement Intense Pluvial).

CONCLUSION
Considering that only half of flood-related damages are likely to be associated with river overflowing in France, this work focuses on the Aude 2018 heavy rainfall event to develop an automatic damage detection method not only nearby but also far away from waterways where intense rainwater runoff can be accountable for. A plot-based supervised classification method was thus achieved relying on the combined search for spectral, temporal and spatial variations in Sentinel-2 change images acquired only twenty days apart with a 10 m spatial resolution. For this, a large learning sample of certified ground truths along with validation data based on photo-interpretation from very high resolution (0.5 m) post-event images from Pléiades satellites were used. Following optimal selection of input indicators and training of a Gaussian process classifier, classification accuracy up to 90% was reached. These results emphasize the need for combined analysis of contextual, spectral, spatial and temporal information and thus confirm the great potential in multispectral multitemporal high resolution optical imagery to detect the full range of heavy rainfall-related damages. A replicable method would be very precious in order to enhance local hydrological knowledge and flood risk management, as well as to quickly identify post-event damages in the context of increasing extreme hydro-meteorological events around the world.