Use of aerial multispectral images for spatial analysis of flooded riverbed-alluvial plain systems: the case study of the Paglia River (central Italy)

ABSTRACT Image processing and classification techniques are widely used for land use definition. They can also provide interesting applications in fluvial geomorphology, for outlining morpho-sedimentary features (bars, channels, banks and floodplain) at various temporal stages, in order to monitor the evolution of river systems. Frequent monitoring is especially important for streams, in terms of flood risk in urban areas. This study shows how techniques of supervised analysis can be applied to river systems, also under particular conditions, like after flood events (when large portions of riverbed and alluvial plain are covered with mud). The procedure starts from the classical photogrammetric techniques, based on multispectral classification, and goes on with post processing operations of pixel aggregation and shadow treatment. The classification also uses the elevation information provided by Digital Surface Model produced by photogrammetry. This paper introduces a new technique of remote sensing in fluvial areas that allows for both the identification and classification of the fluvial features in a post flooding condition. Application of the procedure over time permits the evolution of the fluvial dynamics to be monitored in an accurate and inexpensive way, particularly for flood event conditions which lead to major changes in the dynamics of riverbeds.


Introduction
Remote sensing is an important tool for risk mitigation and management of natural disasters (Van Westen 2000; Baiocchi et al. 2010Baiocchi et al. , 2014; Furtun a & Holobȃc a 2013; . Multispectral images are used in various fields for rapid feature identification over large areas; furthermore, the short time interval in image capturing allows for the temporal analysis and monitoring of feature variations (Coppin & Marvin 1996;Lu et al 2004;Sonka et al. 2014).
Monitoring is essential for understanding the dynamics of river systems, particularly in Italian streams that experience rapid channel changes and human disturbance leading to increase in the risk of flooding in anthropic areas.
There are several examples in the literature of the application of remote sensing techniques to Italian streams (Surian & Rinaldi 2003;Scorpio et al. 2015); nevertheless, the use of automatic procedures of image analysis is not so usual in the scientific literature, especially for Italian rivers. Several authors (Smith 1997;Mertes 2002;Hudson & Colditz 2003;Alsdorf et al. 2007) have used multitemporal satellite data for remote sensing of river stage and discharge in order to produce flood inundation maps, in near real time.
Some examples of automated remote sensing applied to post-event images were used by Jiang and Friedland (2015) to accurately differentiate urban debris areas from non-debris areas. Moreover, the use of remote sensing techniques appears to be fundamental even for flood risk management (Franci et al. 2014;Amarnath & Rajah 2015). Bryant and Gilvear (1999) used remote sensing to quantify geomorphic changes in the Tay River (Scotland), coupling the identification of land-use changes in riparian areas and the identification of bathymetric changes (derived from wavelength rangebandwidth and reflectance) for the riverbed. Moreover, this study highlighted that the detection of bathymetric change is possible only in shallow gravel-bed rivers in the absence of high turbidity. In the present case, multispectral images were captured as soon as the flood event took place in order to take into account the morphological changes, but a different approach had to be used due to the high water turbidity and the presence of mud covering large areas.
This paper shows the results of a study aimed at evaluating the planar view changes which affected the riverbed of the Paglia River (Tiber River basin, central Italy) as a consequence of the flooding event of November 2012, in the Umbria Region. The goal of the study was to develop automatic and/or semi-automatic procedures for the morphological and sedimentary analysis of riverbeds using multispectral aerial images. These procedures would allow for the identification of morphological changes of planar river features and their comparison, by means of the Geographic Information System, with the previous state of the river in order to rapidly map the changes that occurred.
In this paper, supervised classification techniques were applied to a fluvial system. In particular, the complex morphological configuration of the riverbed due to the extreme flood event was analysed. It caused the deposition of sediments, both fine (mud) and very coarse (gravel pebbles), over large portions of the riverbed and floodplain and the extensive removal of riparian vegetation.
Moreover, in this work the elevation data derived from a Digital Surface Model (DSM) generated by a photogrammetric process was applied to improve the results of the automatic supervised classification.
The procedure using these automated tools could be useful for future monitoring of river system changes.

Data
Geological, geomorphological, climatic and hydrological characteristics of the study area The Paglia River (Figure 1) rises in south-eastern Tuscany from Mt. Amiata (1738 m a.s.l.) and is one of the main tributaries on the right-hand side valley of the Tiber River. Its confluence is located in southern Umbria, after about 60 km. The area of the Paglia River Basin is 1320 km 2 .
The Paglia River, especially in the reaches of its lower valley, shows critical characteristics from the morpho-sedimentary dynamic point of view. The Paglia riverbed is in a state of sediment-limited non-equilibrium and is characterized by an intensely active vertical incision. Recently, the thalweg has cut through unconsolidated alluvial sediments and is now actively incising into the bedrock; this is formed by overconsolidated marine clays (formation of 'Fabro's Clays', of Pliocene), which represent the basis of the regressive series of Citt a della Pieve (Ambrosetti et al. 1987;Cattuto et al. 1988). Until to the middle of the last century, the riverbed was very large (up to 600 m wide), moderately braided, with a strongly moving sinuous main channel, set at heights very close to the floodplain (wandering type -Church 1983;Carson 1984). Over time, but especially in the last 50 years, the morphology of the river has changed drastically: the main channel has entrenched and its course has shrunk considerably, leaving most of the bars, which were nonvegetated and active in the past, but now are inactive and covered with dense vegetation of tall trees (abandoned riverbed). Thus, the fluvial course has gradually become a 'single-channel with low-sinuosity' channel type.
Previous studies based on the comparison of historical cartographic documents (maps of the Gregorian Cadastre, aerial photos and orthophotos coming from various periods) attempted to quantify the changes in the last 200 years (Cencetti et al. 2004(Cencetti et al. , 2014. During a period of about 130 years, from 1821 to 1954, the riverbed was affected by a modest degree of narrowing (the average width decreased from 202 m to 155 m); later, from 1954 to 1999, the riverbed shrunk by about 70% (the average width decreased to 46 m in 1999).
In the last decade (until surveying of the situation in the summer of 2012), due to some major flood events, the trend reversed and the riverbed regained part of its original area: the average width increased from 46 m in 1999 to 60 m in September 2012. Presently, large parts of the natural channel migration zones have been destined for agricultural use or for infrastructures and urban areas (Orvieto Scalo -Ciconia).
In the Paglia River Basin (Figure 2), the maximum rainfall occurs in autumn (in November, about 140 mm on average) and the minimum in summer, when the maximum water deficit occurs as a result of relatively high temperatures (average temperature is about 25 C in July-August).
The mean annual discharge is extremely variable, as highlighted in Figure 2; in the period from 1994 to 2015, the values go from a minimum of 0.04 m 3 /s (recorded in 1995) to a maximum of 8.55 m 3 /s (recorded in 2004).
In November 2012, an extreme flood event occurred that produced extensive damage to hydraulic engineering works, agricultural, commercial and manufacturing infrastructures, and movable and immovable property of private owners. The low-pressure covering the area brought exceptional amounts of heavy rain over much of the region: 307 mm in 72 hours in the surroundings of Orvieto, and 230 mm in the area of Marsciano, just north of the area studied. The meteorological situation on 11 November 2012 was described in detail by the 'Centro Funzionale Decentrato -Regione Umbria' in the event report (Regione Umbria 2012): it is represented by a wide and deep 'trough' with a N-S axis, extending from the Norwegian Sea to North Africa (Atlas Mountains). In that situation, the ascending branch of the polar jet stream was located on the Western Mediterranean Sea (Figure 3). In the following hours, a stream of hot air, followed on its western edge by a cold front, produced convective activity that developed over the Mediterranean Sea and headed quickly towards the mainland. The first storms affected Western Umbria, passing rapidly eastward; in the afternoon of the same day a series of rain systems, continuously fed by warm, moist air on the Tyrrhenian Sea, affected Umbria with very heavy rainfall. In the evening, the trough was cut-off and a closed minimum was isolated, thanks to the presence of an Atlantic anticyclone in Northern Europe. In the Mediterranean Sea, in the lower atmospheric layers, the warm moist air stream continued, animated by the minimum slowly transiting in North Algeria. During the following day (12 November) and up to the late evening, this allowed the formation of a series of convective systems producing widespread, intense and persistent rainfall in Tuscany, northern Lazio and Umbria, accompanied by intense electrical activity. Only during the night of 12 November, did the shift of the low-pressure core towards the East interrupt the feeding of the convective systems that exhausted their energy with scattered thunderstorms during the night until the first hours of Tuesday 13 November.
The hydraulic flow of the Paglia River reached values estimated at about 2300 m 3 /s at the hydrometric station of Orvieto Scalo; according to the estimate of the Civil Protection authority, this discharge corresponds to a return period between 100 and 200 years. In Figure 2(B) this value is represented by a red dashed line.

Photogrammetric data
The morphological-sedimentological 'post-event' configuration of the riverbed and alluvial plain, both in terms of planimetry and elevation, was obtained by the photogrammetric processing of 40 aerial multispectral frames (four bands: red, green and blue + near infrared regions (NIR)) acquired in December 2012 by means of an UltraCam Xp (Vexcel) digital metric camera. They are part of a single strip covering an area of about 27 £ 3 km, with a longitudinal overlap of about 60%, and were provided with internal and external orientation parameters.
The digital images were acquired at an average flying height Hm = 2500 m using a camera with a focal length c = 100.5 mm. Therefore, the equivalent scale on the sensor is about 1:25,000 and the pixel size (equal to 6 microns) generates a GSD (Ground Sample Distancepixel size on the ground) of about 15 cm.
The ground control points (GCPs) required for the photogrammetric orientation of the images were measured through a Global Navigation Satellite System (GNSS) survey in Network Real Time Kinematic mode (accuracy 1-2 cm) using a TOPCON GR-5 GPS/GLONASS double frequency geodetic receiver and a Topcon FC-236 controller, using the corrections sent by the GNSS permanent stations network GPSUMBRIA (Fastellini et al. 2008). Figure 4 shows the photogrammetric strip acquired and the GCP distribution (yellow points). The internal and external orientation, the photogrammetric restitution and the Digital Terrain Model (DTM) extraction were performed by means of a digital photogrammetric workstation using the software Socet Set 5.6.0. by BAE Systems. The project was georeferenced in the ETRF2000 Datum with orthometric heights. For mapping purposes, ETRF2000 can be considered very close to WGS84.
The relative and absolute orientation of the frames was performed by aerial triangulation with a bundle adjustment algorithm. The tie points of the frames were automatically identified using an image-matching algorithm. The GCPs, whose coordinates were known from the GNSS survey, were manually identified on the frames and collimated one by one.
The results of aerial triangulation, expressed in root mean square errors, can be summarized as follows: 0.532 for the image (pixels); 0.066 for the X or long (m); 0.068 for the Y or lat. (m); 0.226 for the Z or elevation (m). The DSM extracted in a grid format with 1 m grid steps and with a vertical resolution of 0.5 m was used to generate an orthoimage (mosaic) of the entire study area (Figure 5), and was used to improve the results of the multispectral classification.

Method
Classification A multispectral classification of the orthoimage was performed in order to obtain a land-use map, regarding some classes closely connected to the stream corridor. Some indexes that exploit different surface reflectivity at different spectral bands were used: the NDVI -Normalized Difference Vegetation Index - (Govaerts & Verhulst 2010) that provides information on the health and vigor of vegetation and the NDWI -Normalized Difference Water Index (Gao 1996).
The NDVI index is expressed as follows:  The NDVI computation generates an image characterized by various vegetation classes, typically no vegetation, sparse, moderate and dense vegetation ( Figure 6).
These vegetation classes can be used for multispectral classification with the other Regions Of Interest (ROI) subsequently defined to obtain a land-use map.
The NDWI has an expression similar to NDVI, but a green band is used instead of a red one: In this case, water takes positive values close to 1, while land and vegetation assume lower values (Figure 7).
The orthoimage was classified using a Supervised Pixel Based methodology using the maximum likelihood algorithm (Ahmad 2012) in a test area chosen as representative of the whole study region.
The first step of the processing was the identification of specific land-use classes (ROI) through the selection of some pixels representing the different types of land cover. These classes were chosen by photographic interpretation of the orthoimages, except for those related to vegetation, which were automatically identified by the NDVI calculation. After the ROI creation, their separability was evaluated through functions that calculate the spectral distance (Jeffrie Matusita, Transformed Divergence - Swain & Davis 1978;Mausel et al. 1990;Jensen 1996) because there may be classes  characterized by very similar radiometry. Thus, it is possible to identify and correct any ambiguity caused by different classes with a similar radiometry by adding or distributing differently more representative pixels in the scene.
Using the maximum likelihood algorithm, the pixels were assigned to a specific class according to the maximum probability of belonging criterion.
Buildings and streets were very well defined, and vegetation was automatically classified, obtaining a good separation between meadows and dense vegetation. Instead, some fields without vegetation were wrongly classified as roads or mud because of their very similar radiometry. In addition, the vegetation delineation algorithm, using NDVI for separation between classes, correctly identified the vegetated areas, but had some problems in correctly identifying the leafless trees located within the stream corridor.

Post classification
The results of the supervised classification as described in the above section, although showing good qualitative results, need some improvement in order to be used in the post-processing step.
In particular, it is necessary to improve the following two aspects.
Excessive fragmentation: the reclassified images have an unrealistic fragmentation; the element classes, the object of classification, by their nature are not so fragmented. Moreover, the fragmentation can complicate later analyses relative to the automated measurement of some geomorphological indexes. This analysis is one of the future goals of the research project. Shadows: the orthophoto shows some areas with shadows due to the time of day (afternoon) when the images were acquired. Such shadows continue to be present in the final result.
Therefore, it is necessary to identify which rules with their respective thresholds are required to eliminate or minimize the issues listed above.
The proposed procedure can be explained by the following steps (see Figure 8).
(1) Conversion of cells classified as shadow (by the multispectral supervised analysis) into NODATA.
(2) Cell aggregation using a moving window of 9 £ 9 cells with mode function. At this step, the r.neighborhood module of GRASS GIS was used. With this tool to each cell, a value was assigned equal to the statistic mode (the most frequent value in a data set) of the cell values in the 9 £ 9 moving window. (3) Filling the NODATA pixels from shadow using a moving 33 £ 33 window always through r.neighborhood. The moving window dimension of 33 £ 33 was applied as this is the smallest area able to fill all the NODATA present in the classified image. Figure 8 shows how the procedure is able to solve the issues related to excessive fragmentation and the presence of shadow. Figure 9 shows the image classification, final output of the proposed procedure, compared to the original orthophoto.
A qualitative pixel-based control shows that moving water, shallow water, gravel and fallen vegetation were correctly identified. For these classes the classification errors were very scattered and confined to areas where the NODATA were filled with the moving 33 £ 33 window.
Some aspects that are highlighted by arrows in Figure 9 still remain unsolved: some pixels adjacent to the channel were classified as 'fields', whereas they were actually covered by mud deposits (blue arrow); the classification did not distinguish tall tree riparian vegetation from low tree vegetation or crops (orchards, olive groves and vineyards) located outside the stream corridor (red arrow).

Reclassification improvement using DSM images
The DSM extracted by means of photogrammetric techniques was used to improve the image classification affected by the problems described above. DSM was detrended using the R (R Core Team 2013) library gstat (Pebesma 2004). Such procedure defines, through a trend surface analysis, the best-fit plane for the DSM. Detrending the DSM consists in removing the longitudinal valley slope with the effect of 'levelling the DSM to the horizontal' (Figure 10).  This procedure allowed two elevation levels to be selected over the detrended DSM. The elevation of 7.2 m identifies the upper limit of the stream corridor, composed by the recent riverbed (channel and gravel/nonvegetated bars) and inactive bars, and the elevation of 11.5 m identifies the lower limit above the tall vegetation present. Figure 10 shows the performed detrending procedure; part B shows the detrended DSM. In Figure 10(B), the dashed black line highlights the area where the elevation is below the level of 7.2 m, attributed to the stream corridor.
The result of supervised classification was corrected according to the following criteria: stream corridor areas (elevation 7.2 m), previously classified as 'fields', were reclassified as 'mud'; areas with an elevation above 11.5 m, previously classified as 'generic vegetation' or 'demolished vegetation', were reclassified as 'tall trees'.

Discussion
The results of the supervised classification consist in a grid raster map, having a cell dimension of 1 m, where the original orthophoto is classified in 8 different classes, as follows: moving water, shallow water, gravel, mud, field, vegetation, fallen vegetation and high trees. Figure 11 shows the results with and without the correction procedure based on the elevation from detrended DSM, as described above. The original classification is shown on the left, while the classification corrected by DSM is shown on the right: a new class 'tall trees' has been added (arrow in Figure 9) and some areas which had been erroneously classified as 'fields', are now classified as 'mud' (see the oval line in Figure 11).
The results analysis was carried out both qualitatively and quantitatively. Quality control was carried out at the pixel level and showed how moving water, shallow and cloudy water, gravel and fallen vegetation were correctly identified. However, some errors still remain in limited areas where the NODATA were filled with the r.neighborood algorithm.
The most difficult class to be recognized was the mud class, often confused with the field, and vegetation classes, not distinguishable from tall trees, shrubs and farm. These issues were resolved by using the DSM correction.
Quantitative analysis of the results was carried out comparing the supervised classification and a classification independently performed using a classical stereoscopic photointerpretation technique. It should be noted that such a photo interpretation is the only method available for a validation procedure.
Photointerpretation was performed with the aim of classifying the main features of the river, distinguishing 5 classes: channel, gravel-cobble bars, abandoned channel and vegetated bars, floodplain inundated by the November 2012 event and floodplain not inundated by the November 2012 event. Figure 12 shows the results of the photointerpretation used for the validation procedure. The classes obtained using the supervised classification and those obtained from the photointerpretation are listed in Table 1. Figure 11. Classification results using the detrending DSM procedure. The oval line shows some areas incorrectly classified as 'fields' and now classified as 'mud.' Figure 13 shows the results of the two independent techniques superimposed on the same test area. Since the two analyses are based on different assumptions and modalities, they define different classes, and, therefore, it is not possible to make an accurate comparison.
The performance evaluation was carried out matching the frequency of classes of the supervised method (left column in Table 1) with the classes from the photointerpretation (right column in Table 1) and considering the compatibility among the different combinations ( Figure 14).
Inside the channel it is expected to find moving water and shallow water, but a high percentage of gravel, mud and fallen vegetation is also admissible. Conversely the presence of vegetation and fields is attributable as errors.
Histogram A in Figure 14 shows that the supervised method used in the channel already provides quite accurate results even without the DSM correction: 66% moving water and 17% shallow water  (correct classification), 11% gravel and fallen vegetation (acceptable classification) and 4% for fields (wrong classification).
The DSM correction slightly improved the final results converting the class 'field' (not acceptable) to 'mud' (acceptable), because in this specific case the major classification errors occur between them.
Areas assigned from photointerpretation to gravel-cobble bars were considered correctly classified when assigned to the following classes: gravel, mud, fallen vegetation, generic vegetation and shallow water; they were considered a misclassification when falling into field and moving water classes.
Also in this case, the already good results were improved applying the DSM correction which allowed the errors to be eliminated, transforming the class 'field' into 'mud'.
In the abandoned channel areas, the classes considered correct were: all classes of vegetation, gravel, mud and shallow water that may be present in the secondary channels. Such an assumption is even more correct in this case concerning a riverbed just after an extreme flood event. On the other hand, the classes moving water and fields are considered as misclassifications.
The results of the supervised classification, without the DSM correction, are the following (histogram C in Figure 14): vegetation (generic and fallen)mudgravel and shallow water (correct classification) 61%; moving water (misclassification) 0% and field (misclassification) 39%. In this case, using the algorithm, it was somewhat difficult to discriminate between 'field' and 'mud' classes, because they have very similar radiometric values; for this reason, an important error is introduced, which, using the DSM correction, was significantly reduced, from 39% to 27%.
Furthermore, the DSM improved the classification results via identification of tall trees class (13.6%), which is distinguished from the general vegetation class: the vegetation decreased from 21.4% to 11%, and fallen vegetation from 15% to 11.8%.
In the classes 'inundated floodplain by the November 2012 event' and the 'not inundated floodplain by the November 2012 event', the supervised classification provided satisfactory results. In fact, the frequency of the only incorrect class (moving water) was 0%.
The supervised algorithm worked very well in the floodplain where there was no critical situation.
Indeed, it can be stated that in the floodplain the supervised technique supplied a more accurate classification, with respect to the photointerpretation, in identifying areas characterized by deposition of fine sediments (mud), coarse sediments (gravel) and water stagnation (shallow water). Moreover, the identification of such elements using classical techniques requires a significant effort on the part of the operator, as this is a time consuming and highly detailed job.
On the other hand, the DSM correction was showed to be essential for identifying tall trees class, typical of riparian areas since without the DSM correction tall trees were classified as generic vegetation and grouped with shrubs and crops (mainly vineyards).

Conclusions
The processing of multispectral aerial images can provide useful information in fluvial geomorphology and for evaluating the morphological and sedimentary features of floodplains.
The analysis covered a reach of the Paglia River and was carried out using digital aerial photographs taken in the days just after the extreme flood event that occurred in November 2012.
The particular configuration of the riverbed made the analysis more difficult for two reasons: wide sections of riverbed and floodplain were covered by both very fine (mud) and very coarse sediments (gravel and pebbles); riparian vegetation, characterized by the autumnal slow growing rate, had been largely removed and/or covered by mud.
The supervised classification of multispectral images provided satisfactory results in terms of ROI separability, with a limited percentage of classification errors due to the spectral similarity between mud and field classes.
The classification was refined in post-processing using the r.neighborhood module of GRASS GIS, allowing for the filling in of the missing value pixels in areas covered by 'shadows' and reducing the aspect of excessive fragmentation.
However, some critical issues, linked to the effects of the flooding that occurred in November 2012, still remain open: the supervised classification algorithm is not able to accurately distinguish between fine sediments (mud) and fields classes, because they have very similar radiometric values; the classification does not distinguish the tall trees class in the riparian vegetation from the other vegetation or crops (orchards, olive groves and vineyards), which are present outside of the stream corridor.
In order to reduce such problems, a correction criterion was introduced; it is based on height values from a DSM produced in the photogrammetric restitution.
A validation procedure was carried out comparing the output of supervised classification (without and with DSM correction) to a classical stereoscopic photointerpretation technique.
The validation output ( Figure 14) showed that the good results obtained with the supervised automatic classification after post-processing were refined and improved by the DSM corrections.
Examining the results in detail, it is possible to deduce the follows: in the channel (stream corridor), the percentage of erroneous classifications was reduced from 4% (without DSM) to 0% (with DSM); in areas of gravel and cobble bars (stream corridor), the percentage of erroneous classifications was reduced from 15% (without DSM) to 2% (with DSM); in the floodplain, where the automatic supervised classification provided the best results, the errors were negligible and the multispectral classification technique provided a better identification of classes, like deposits of fine sediment (mud), coarse sediment (gravel) and shallow water, than the classical photointerpretation method; in more complex areas such as the abandoned channel areas, the DSM corrections reduced the misclassification from 39% to 27%.
Moreover, the application of DSM after the supervised classification was essential for identifying tall trees class, typical of riparian areas, that otherwise would have been grouped with a generic vegetation type.
The automatic or semi-automatic procedure presented in this paper can be easily replicated and provides a useful tool for monitoring riverbeds and related alluvial plains.
It is important to note that if the aerial photos are taken long after the extreme flood events, the supervised classification would be facilitated as the issues related to the extended presence of mud and gravel in not common areas (such as floodplain and abandoned channel) and the green-brown colour of the foliage would be less relevant.
As a future development, a possibility could be to use the results provided from the supervised analysis after the post-processing procedure and the DSM correction as input for further processing, aimed at evaluating parameters and geomorphological indexes (e.g. number of channels, sinuosity, etc.), useful for the characterization of the overall state of the river.