*Article* **Mapping Impact of Tidal Flooding on Solar Salt Farming in Northern Java using a Hydrodynamic Model**

#### **Anang Widhi Nirwansyah 1,2,\* and Boris Braun 2**

Received: 19 August 2019; Accepted: 10 October 2019; Published: 12 October 2019

**Abstract:** The number of tidal flood events has been increasing in Indonesia in the last decade, especially along the north coast of Java. Hydrodynamic models in combination with Geographic Information System applications are used to assess the impact of high tide events upon the salt production in Cirebon, West Java. Two major flood events in June 2016 and May 2018 were selected for the simulation within inputs of tidal height records, national seamless digital elevation dataset of Indonesia (DEMNAS), Indonesian gridded national bathymetry (BATNAS), and wind data from OGIMET. We used a finite method on MIKE 21 to determine peak water levels, and validation for the velocity component using TPXO9 and Tidal Model Driver (TMD). The benchmark of the inundation is taken from the maximum water level of the simulation. This study utilized ArcGIS for the spatial analysis of tidal flood distribution upon solar salt production area, particularly where the tides are dominated by local factors. The results indicated that during the peak events in June 2016 and May 2018, about 83% to 84% of salt ponds were being inundated, respectively. The accurate identification of flooded areas also provided valuable information for tidal flood assessment of marginal agriculture in data-scarce region.

**Keywords:** mapping impact; tidal flood; hydrodynamic model; solar salt farming

#### **1. Introduction**

Globally, coastal flooding have been devastating events causing cost for human environment, increase property damage, and around 20 million people are exposed to present high tide levels and 200 million to storm tide levels [1,2]. Currently, the Intergovernmental Panel for Climate Change (IPCC) report suggests that the global mean sea levels will increase 36–71 cm by 2100 based on Representative Concentration Pathway (RCP) 4.5 mid emissions scenario [3]. This situation may increase the vulnerability of coastal regions, especially of cities, due to demographic trends and economic expansion [4,5]. Meanwhile, in developing countries where various types of agricultural activities dominated the local economies, the impact is either ignored or simplified using rough estimates because of low expected losses [6–8]. Moreover, local types of agriculture such as solar salt production in tropical countries are also facing the impact of tidal flooding in particular location, which is overlooked in the global discussion. This type of agriculture, however, has the potential to generate revenue from salt in various aspects, not only in terms of salt product quality, but also for tourism, or even partly for coastal research centers [9].

Currently, solar salt production is acknowledged to be a marginal economic sector, especially in Indonesia [10]. It manually operates through traditional technology by using solar evaporation [11], locally referred to as '*maduranese*' method. The process starts in the saltpan. Seawater is let into the first and largest concentrating pond, or concentrator, through an inlet [12]. Most of the salt farmers

are producing sea salt only during the dry season (April–October). The timing of the process highly depends on weather conditions, as rain reduces salinity and clouds decelerate evaporation [13]. Tidal flooding upon solar salt farming areas frequently occurs during high tide in the production period and thus threatens the production and distribution processes (see Figure 1). Between the 5th and 8th of June 2016, high tide flooding inundated hundreds of hectares of salt ponds in Cirebon, West Java, which is one of the major producer sites for salt in Indonesia [14,15]. Concurrently, the similar astronomic phases during the inundation events, showed the increase of high water level and low water level [16]. Temporarily, there was another nuisance flood between 23rd and 25th May 2018 along the north coast of Java (locally referred to '*Pantura*') during high tide [17]. Both events were narrated widely in both electronic and printed media.

**Figure 1.** Setting of traditional solar salt production area including: (**a**) salt evaporation pan and channel; (**b**) nearest pond to the sea; and (**c**) inundated pond due to high tide during fieldwork on 7 January 2018, 11:00 UTC (17:00 local time).

Tidal hydrodynamics in the Java Sea are complicated, due to their rough shallow bottom topography, diverse types of coastlines, and the interference of tidal waves propagating from the Pacific Ocean, Indian Ocean, and South China Sea [18]. Koropitan and Ikeda [19] previously investigated the implication of the barotropic tides into four tidal harmonic constituents using a three-dimensional (3D) hydrodynamic model combined with observation data and have suggested that the semi-diurnal *M*2 component dominates over Java Sea. Several studies show that wind factor has a minimum contribution to tidal propagation in Java Sea [18–20]. Tidal flooding in the northern part of Java periodically rises in July and August during the East Monsoon period [21]. In recent years, the tidal inundation comes not only at high tide but even at the regular tide in some areas along *Pantura* [22]. Furthermore, the local economy, such as salt production which is dependent on coastal conditions, is eventually disrupted during these events.

Tides caused by the gravitational effects of sun and moon are periodic and very predictable [23–25]. Tidal floods (also defined as "nuisance" flooding) are occurring more often during seasonal high tides or minor wind events, and the frequency is likely to escalate intensely in the forthcoming decades [26]. Currently, the impact of tidal flood is usually modeled using planar approach in geographic information

analysis. This approach assumes areas lower than a particular elevation to be inundated utilizing digital elevation model (DEM) and geographic information system (GIS) [27,28]. Geophysical processes including bottom friction or motion transfer are not considered in these particular models [29]. Ultimately, the uncertain behavior of the coastal system during a coastal flooding event is still a challenge in this model [30]. Previous work on GIS modeling presents various resolutions of DEM, which suggest using high resolution of elevation data to increase accuracy [24,31,32]. However, it is the hydrodynamics of the tide that is responsible for the size of the tide range, high and low waters momentum, and tidal characteristic, as well as the speed and timing of the tidal current [33]. An integrated approach considering both aspects is recommended, particularly for smaller areas and in cases where details are essential [29,34,35].

In this study, a simulation of the tidal flooding on salt production area is presented. Furthermore, investigations of high tide flooding in the salt production area of Cirebon are not available. The method, implemented in past events using hydrodynamic model with additional inputs on behavior of the coastal system during a tidal flooding event. This approach retrieves the values of the flood impact on the parcel of salt pond from two-dimensional (2D) (floodplain flow) tidal inundation simulation. Against the above background, the objectives of this study are to: (i) validate the tidal flooding that happened in June 2016 and May 2018 using a hydrodynamic model; (ii) analyze the highest tidal elevation and factors associated with flooding using tidal constituent and wind; (iii) plots tidally inundated area upon solar salt production area by considering the spatial distribution and the depths. Finally, this research offers better accuracy of analysis on the distribution of tidal flooding in salt production areas within limitation of tidal flood data.

#### **2. Location of the Study Area**

Cirebon is located 6◦30 —7◦00 S and 108◦40 —108◦48 E. It covers an area of 990.36 km2. Administratively, Cirebon is a part of the West Java Province and is bordered by the Java Sea, and by Indramayu in the north, Kuningan in the south, Central Java Province in the east, and Majalengka in the west (Figure 2). It is a typical lowland with an average elevation of 0–25 m and covers 64,500 hectares. Cirebon connects the capital city of Jakarta with major cities in central and east Java. Cirebon has 40 districts, 424 villages, and 12 sub-districts. Based on BPS [36], the port city has an approximate population of 2.1 million, with 2205 inhabitants per km2 and a population growth averaging of 1.28% per year.

**Figure 2.** Geographical situation in Cirebon within typical coastal lowland adjacent to Java Sea and simulation coverage area.

Along with agriculture, salt production is shaping the local economy in the coastal region of Cirebon. The salt ponds cover 7819.32 hectares and provide jobs for 3707 people, such as pond owners, salt workers, and intermediaries [37,38]. Salt production predominantly takes place in low-lying areas dominated by alluvial deposits alongside with mangrove ecosystems. The salt production period in Cirebon begins during southeast monsoon. Most of the farmers start to store seawater in April, May, or June depending on the weather. They start to collect the brine daily and generate yields 0.5-1 ton/hectare/day. The dry season begins in March and ends in September, with a mean temperature of 32.8 ◦C, while the rainy season usually lasts from October to February, with an average rainfall of 1300-1500 mm/year and an average temperature of 24.2 ◦C [36]. The tidal regime is dominated by a mixed semidiurnal type and experiences two high and two low tides of different scales each lunar day. This tidal characteristic dominates the tidal cycle along Java sea [39].

The Java Sea is mainly identified as shallow water within roughly rectangular morphology, a mean depth of 50 m, a length of 950 km, and a width of 440 km [19,40]. The tidal range in the Java Sea is approximately 1.2–2 m, with peak values around Surabaya, Madura, and Bali [41,42]. The Java Sea is strongly governed by the monsoon climate. The northwest monsoon (NWM) reaches its peak between December and February (DJF) and it is usually characterized by frequent rainfalls and windy periods, while the Southeast monsoon (SEM) extends from June to August (JJA) and is usually characterized by much lower rainfalls [43].

#### **3. Materials and Model Description**

#### *3.1. Data Acquisition*

This study has used several data to simulate post-events of tidal floods. Firstly, the bathymetry and land topography information for the domain areas were handled as two main inputs for the model. Bathymetry of the Java Sea has been generated using gridded national bathymetry of Indonesia (BATNAS) provided from the Geospatial Information Agency (we referred to BIG: '*Badan Informasi Geospasial*' in Indonesian) (http://tides.big.go.id/DEMNAS/) within a 6 arc-second resolution. This data has been produced through the inversion of gravity anomaly of altimetry by adding sounding data carried with single and multi-beam surveys, which has better resolution in coastal areas than GEBCO (30 arc-second) [44,45]. Land topography data was resolved using the DEMNAS (0.27 arc-second resolution) also from BIG. DEMNAS is national seamless digital elevation data which already constructed within assimilated data of IFSAR (5-m resolution), TerraSAR-X (5-m resolution) and ALOS PALSAR (resolution 11.25 m), by adding stereo-plotting mass-point data [45]. This research draws on a previous approach by Tehrany [46] and Zalite [47] to utilize detailed topographic data for flood models. Here, five subsets (1309-12, 1309-14, 1309-21, 1309-23, and 1309-24) of DEM within the 0.27 arc-second spatial resolution were employed, and merged into a single raster data using GIS.

Secondly, tidal level data for Cirebon waters were captured from local tidal station in Cirebon port operated by BIG. In this case, we used hourly data of both selected period of simulations. As additional input, the simulation included wind data in the form of wind velocity, and wind direction of the Jatiwangi station. This data was extracted from OGIMET online meteorological database (https://www.ogimet.com/). Lastly, this model involved the latest (updated 2015) of the salt parcel dataset taken from the ministry of marine and fisheries affair during the salt inventory mapping project together with the national geospatial agency and PT. Garam [38]. At this point, salt parcel can clearly separate each pond by scale of 1: 15,000 in polygon format and was referenced into the World Geodetic System 1984. Details of data needed on this research are mentioned in Table 1.

**Table 1.** Data requirement for salt farming impact due to tidal flood.

#### *3.2. Model Setup*

This research applied a numerical hydrodynamic model (HDM) to forecast run-up and tidal inundation in the salt production area in Cirebon. HDMs originated through resolving Laplace Tidal Equations and using bathymetry data as boundary conditions [48]. The module of MIKE 21 package was used, as one of the most widely used hydrodynamic model in computation by Danish Hydraulic Institute (DHI), including the assessment of hydrographical sequences in non-stratified waters, coastal flooding and storm surge, inland flooding, and overflow [49–51]. The MIKE package also represents user-friendly GIS interfaces and provides better possibilities to simulate the flooding using elevation data and bathymetry [52].

The model employed MIKE 21 Flexible Mesh (FM) to simulate water levels and tidal floodings in selected events. These tools were utilized using the input of tidal gauge records to indicate the spatial variability of tidal flood characteristics of two events. The model was applied for two separated months, June 2016 (event A) and May 2018 (event B), which covered the occurrence of the selected tide flood events. The unstructured triangular mesh with 87,103 nodes and 170,501 elements was generated in the simulation and covered 11,515.20 km2 (see Figure 3). The mesh file in ASCII format included information of the coordinates and bathymetry for each node point in the mesh [50]. The grid dimension differs by 2800 m in the northeast ocean boundary, the smaller grid size around 450 m in the outland, and 120 m in the inland area.

**Figure 3.** Mesh generation of Cirebon waters and part of Java Sea, bathymetry (in meters), comprising of 87,103 nodes and 170,501 triangular elements.

The tidal flood events, which were triggered by high tides, were simulated through forcing tidal elevations at open borders, winds, and temperatures. The tidal height was calculated by hourly local station measurement using a harmonic approach. This classical harmonic analysis represents the tidal forcing as a set of spectral lines, demonstrating the predetermined set of sinusoids at specified frequencies [53,54]. This stage resulted in nine tidal components (*M*2, *S*2, *N*2, *K*2, *K*1, *O*1, *P*1, *M*4, and *MS*4) that correspond to specific physical phenomena such as the period of the moon around the earth or friction against the seabed in shallow seas.

Following step of this model was to include the wind energy and set tilt facility to declare water level correction along the boundary of the waters using Navier-Stokes equations. These equations deliver an appropriate model for wave overtopping and overcome sophisticated hydrodynamics, including wave breaking and its theoretical limitation [55,56]. Furthermore, the flood and drying (FAD) ability of this model assisted the water run-up simulation and executed the inundation process of high tides. This scheme has been alternated to describe the coastal situation, where it can be flooded at one time but dry at other times. This study use recommended value in Thambas [50], *hdry* = 0.005 m, flooding depth *hflood* = 0.05 m and wetting depth *hwet* = 0.1 m.

It should be noted that the elevation data is used in this simulation without considering water surface evaporation. Due to the high complexity of the site study and the limitations of the model, the following steps were considered for the tidal flood simulation: bottom friction is based on Manning's approach, with the ranges of friction coefficients from 40 for water to 32 for land [57]. A Manning number in the range 20-40 m 1/3/s is typically applied with an advised value of 32 m 1/3/s if other data is unavailable [58]. The Manning number relates to the flow path and peak time of flooding and does not have a significant effect on flood distribution and depth [59]. Furthermore, the simulation included horizontal eddy viscosity using the Smagorinsky type within a value of 0.28 [49].

After the work with MIKE had been completed, the result was exported using "MIKE2Grid" for further spatial processes. This produced an ASCII file that is readable in ArcGIS. Importing this file and reorganizing the classes produced the inundation map of MIKE 21. Finally, we superimposed the grid data to the salt parcel dataset and used the tidal simulation as a basis in inundation analysis process. This plot dataset has a shapefile format, which is suitable for further handling in ArcGIS. The inundation map was created by using grid data of both simulations in ArcGIS. The two-top water levels of the selected simulations, which were identified with the maximum value in the data series that was used as a benchmark for inundation analysis. In this research, certain assumptions were made, such as that no precipitation data inputs were used during the period of the incident as it may raise the inundation level, no sea-level rise and land subsidence were considered in the simulation, as there is still no strong local evidence of both factors in the research location. The overall steps of this research are summarized in Figure 4.

**Figure 4.** Diagram of simulation of tidal flood and impact mapping procedure including validation process.

#### **4. Results**

#### *4.1. Validation of Tidal Simulation*

To get a sound validation of our model for tidal simulation, this study used tidal data of the open boundary model that were obtained from the global tidal model. Following the step from Ningsih et al. [60], the model compared the results of simulated sea level from MIKE 21 with tidal station in Cirebon from BIG and also the global tidal model TPXO9 (this model can be accessed on http://volkov.oce.orst.edu/tides/global.html). TPXO9 is the latest version of TPXO-series [61,62], which includes global tidal solutions with 1/6◦ resolution that fit, in least-square, both the Laplace's equation and also the long track averaged data from TOPEX/Poseidon and Jason (on T/P tracks since 2002) [61,63]. Tidal constituents of the tide record from BIG tidal station perform comparable values with tidal constituents from the global tide modeling TPXO in Indonesian waters [64]. At this point, the wind factor was excluded and focused on gravitational force only. Moreover, the tidal current velocity was verified with Tidal Model Driver (TMD). This free MATLAB package offers harmonic constituents for tide models, making predictions of tide height and also currents [65]. Verification points of the selected simulations (event A and B) are located in Tawangsari, Pangenan, and Bungko (as pointed P1-P3 in Figure 5). The model exposed the statistical correlation using the Pearson value (r) of the three locations with general tidal model of TPXO9, and presented the value of the Root Mean Square (RMS) error. The RMS error was calculated with:

$$\chi\_{\rm RMS} = \sqrt{(\sum\_{i=1}^{n} x\_i^2)/n} \tag{1}$$

where *xi* is the *i* th point of the chosen area, were calculated in the region 6◦S–7◦S, 108◦E–109◦E (the Java Sea). The overall locations verified excellent correlations between simulated outcomes and those of TPXO9 and TMD.

**Figure 5.** Salt production area for tidal simulation and validation points (P1-P3).

The Pearson correlation of the tidal height simulation is in the range of 0.903–0.908 for event A and of 0.848–0.903 for event B with the RMS Error within approximately 0.069–0.100 m. For the *u*-velocity component, the correlation shows a coefficient around 0.833–0.965 with RMS Error of about 0.023–0.0196 m/s for event A. For event B, the correlation coefficients range between 0.570 and 0.877 with a RMS Error around 0.019–0.190 m/s. Furthermore, the υ-velocity component shows a good agreement of the correlation coefficient, whereas the number of the correlation is about 0.683–0.824 with RMS Error about 0.040–0.061 m/s. Although there were some inconsistencies between *u*-velocity components in MIKE 21 simulation on event B, overall, the simulation results managed well with the TMD data. Lower values of RMS Error suggested the appropriate model to the data points; likewise, values of Pearson close to the maximum point of one (value of 1) indicated that the model has a strong correlation to the water level data [66,67]. Detailed values of Pearson Correlation and RMS Error between simulation and global tide model of TPXO9 for water elevation and those TMD for tidal velocity elements are shown in Table 2. The illustration for tidal height (ζ), *u* and υ-velocity at verification points can be seen in Figures 6–8.

**Table 2.** Pearson Coefficient (r) and Root Mean Square (RMS) Error between simulation rates in MIKE and TPXO9 for water elevation and velocity components from TMD.

**Figure 6.** Comparison of water level between simulation and TPXO on June 2016 (left) and May 2018 (right) at (**a**) P1—Tawangsari; (**b**) P2—Pangenan; and (**c**) P3—Bungko.

**Figure 7.** Comparison of *u*-velocity component between MIKE and TMD on June 2016 (left) and May 2018 (right) at (**a**) P1—Tawangsari; (**b**) P2—Pangenan; and (**c**) P3—Bungko.

**Figure 8.** Comparison of υ-velocity at verification points and TMD on June 2016 (left) and May 2018 (right) at (**a**) P1—Tawangsari; (**b**) P2—Pangenan; and (**c**) P3—Bungko.

Spring tides (at the new moon phase) appeared in between the flooding events of event A, which was recorded from June 1st–5th, 2016. In event B, the high tides were tracked from the 23rd to the 28th of May during the end of the moon phase. In both of these conditions, the rise of waves is a natural phenomenon due to moon force (*M*2) [20]. However, on the three tide validation points, the increase of high water levels (HWL) and low water levels (LWL) during the inundation events can clearly be compared to the same astronomic phases tides data before and after the events (see black boxes in Figure 6 above). Moreover, the horizontal (*u*) velocity on three sample locations show typical performances within range of 0-10 m/s and vertical (*v*) velocity in range of 0-28 m/s (see Figures 7 and 8). Here, it can be seen that there are differences velocity level between the simulation and the TPXO model in three sample locations, which are explained in the previous section.

In the next step, the wind air pressure data were engaged in the model as an additional factor in tidal propagation. Wind data from OGIMET has been entered into the simulation for both selected periods. Adding hourly wind data into the model shows minor differences of amplitudes and phases of tidal constituents. The simulations show that there is an insignificant difference in the tidal pattern due to the relatively small effect of wind, as the velocity of wind dominantly emerged from the north with a maximum speed of 6.8 m/s during event A and an average velocity around 0.81 m/s. For event B, an extreme increasing velocity of 60.48 m/s is recorded in the simulation. Ultimately, the average wind speed is approximately 0–1.49 m/s. Here, the assumption has been made that typically calm wind in both selected periods has a minor impact on tidal floods along the coast.

Wind velocity confirmed a minimum correlation to the water level (Pearson correlation 0.1902-0.1905 for event A and –0.021 to –0.031 for event B). The negative correlation probably relates to the minimum velocity of the wind, as OGIMET provides hourly datasets in both periods of simulation. The typically calm wind during these periods provides a better situation for evaporation in the salt production process. Nevertheless, high tide continually increases the potential of inundation in coastal areas where salt production takes place. At the same time, the water elevation of the simulation also shows significant correlation to the observation data from tide gauges (Figure 9a,b). Figure 9c,d present the peak tide levels for both events. Here, the series of surface water (*z*) data from MIKE 21 simulation are also used to determine the nine tidal components. Here, Table 3 presents the variability of tidal constituents in both periods of the model.

**Table 3.** Tidal amplitudes constituents in Cirebon resulted from simulation.

Based on the calculation, the wind involved in the simulation showed correlation values to tidal gauge observation 0.761 with RMSE of 0.134 m for event A and 0.79 with RMSE 0.120 m for event B. As a result, surface elevation at the peaks of both tidal events in Cirebon reached 0.38 m (event A) on the 2nd of June 2016 12:00 UTC and 0.40 m (event B) on the 25th of May 2018 11:00 UTC. Gurumoorthi and Venkatachalapathy [68] and Pugh [69] mentioned that the relative importance of diurnal and semi-diurnal components differ with geographical position and can be calculated by the formulation factor:

$$F = (O\_1 + K\_1) / (M\_2 + S\_2) \tag{2}$$

In Equation (2), the average constituents confirm the typical mixed, predominantly semi-diurnal tide within 0.73 and 0.68 for both events A and B. Considerably, the amplitudes of semidiurnal tidal constituents were higher than the diurnal tides.

**Figure 9.** Scatterplot and RMS error of simulated surface water elevation with tide gauge observation on: (**a**) event A; (**b**) event B within peak level of water during simulation on; (**c**) 2 June 2016 12:00 UTC (steps 36); and (**d**) 25 May 2018 11:00 UTC (steps 587).

#### *4.2. Maximum Tidal Height and Exposed Salt Production Area*

Understanding the impact of tidal flood dispersal on the coastal area demands a model of inundated area caused by tide water level [24]. Equilibrium flood mapping or the "bathtub" approach compare the maximum total water level and ground height. At those places, where the land is lower than the expected maximum water level, it will be flooded [70]. The expected water depth for each salt pond can have major implications for tidal management, especially for vulnerability measurements as damage is often associated with the depth of inundation and its duration. The simulation showed that both tidal floods were forecasted to be generated by meteorological factors. Here, *M*2 tidal response provides the dominant influence (which both events record 0.143) in amplitudes of Cirebon waters. As a result, the surface elevation during the maximums of both tidal events have been exported in ArcGIS using Mike2Grid tools and visualized water level and the spatial distribution of the inundation upon salt production area (see Figure 10).

**Figure 10.** *Cont.*

(b)

**Figure 10.** Simulation of water level during peak period of: (**a**) event A and (**b**) event B.

The simulation map shows that tidal inundation occurs along the coastline of Cirebon during peak water levels. The grid data was superimposed with detailed DEMNAS to investigate the impact of tidal flooding upon solar salt production land. A reclassification process in GIS elaborates the tidal dynamics and flood depth upon salt pond in the study area. Here, each salt parcel has a single value of depth level through spatial joint between both vector types of water level and parcel of salt pond datasets (Figure 11). Thus, this results in the appropriate value for inundation for each pond that has been impacted by the tidal occurrence for both events.

(a)

**Figure 11.** *Cont.*

**Figure 11.** Estimated inundation level for each pond of during the highest tide of (**a**) event A and (**b**) event B (each parcel contains single value of inundation depth).

As mentioned in the previous section, the results regarding the submerged area of around 0–0.38 m for (event A) and 0–0.40 m for (event B), have significantly affected the salt production areas in Cirebon. Based on previous maps (Figure 11), it can be seen that around 1990.55 ha of salt production pond in Losari were inundated during event A and 1992.07 ha during event B. This district is also recorded as the most impacted area due to both tidal flood events (99.92% and 99.99% of total cultivated area in Losari). The salt production area of Gebang, which is located in the west part of the study area, has also been flooded up to 816.32 ha (100%) during the events A and B. At the same time, a slight increase of flood coverage has occurred in Kapetakan due to both tidal events. During the peak level of event A, almost 56.15% or 1,538.96 ha were exposed to tidal floodings, and 57.22% or 1,568.34 ha suffered inundation according to our simulations. In the middle part, tidal heights of both selected simulations in events A and B have submerged Suranenggala, Gunungjati, Mundu, and Astanajapura to a lesser degree in terms of total area, but with more significant percentages of inundated salt pond (approximately 49–99% in both events). The model presents the areas of inundation on A and B events as it is drawn in Figure 12.

**Figure 12.** Tidal flood accumulative distribution area of inundation in Cirebon due to high tide at event A and event B.

This model shows relatively small differences in terms of the affected area during both simulation periods, as the wind factor has less impact and relatively similar water level. The simulated inundated salt production areas for A and B peak events are estimated to be 6489 ha and 6570 ha, respectively, which equals to 83 % and 84.2 % of the total salt production area in Cirebon. Overall, the peak depth of > 0.35 m dominates the tidal flood sequence, with 41.9% and 45.5% of the area being inundated to such a degree. This depth level has a significantly larger effect in destructing dikes compared to a lower flood level. Meanwhile, around 16–17 % of ponds are relatively safe from these events, as they are located further inland on higher elevations. The less impacted area is estimated within >0-5 cm depth, which covers 1.7% (event A) and 4.1% (event B) of the total inundated area. During the flood, salt production was postponed and stopped until the water receded. It has to be noticed that higher flood levels also take a longer time to recede, which prolongs the preparation and the pre-production process, thus worsening the impact of the floods upon the salt production. Estimations of the salt production area that have been inundated based on simulation is presented in Table 4.

**Table 4.** Estimated areas and percentage of salt production in Cirebon affected by tidal flood in selected events (area in hectare).

#### **5. Discussion**

This paper presents simulation of the inundated areas upon salt farming due to tidal flood events. Tidal floods that occurred upon salt production area were triggered by the high tide events in Java Sea. That two events were similarly situated to tidal incidents located along the southern part of Java adjacent with Indian Ocean [71,72]. Both periods studied present similar tidal elevations. Inundation dominantly occurred in the western and eastern parts of the region. The total impacted salt production area in both events were about 6489.4 ha (event A) (83%) and 6570 ha (84%) (event B). As illustrated by Châu [73], the inundated area may be overrated based on data characteristics and the methods employed. Different resolutions of DEM may result in different total area of inundation. Furthermore, bathymetry, wind velocity, and Manning coefficient also correlate with the hydrodynamic process of tidal forcing.

This method, which relies on tidal characteristics and hydrodynamic parameters, leads to a usefulness tidal flood mapping for salt production areas. This idea improves the marine environment evaluation through cost-effective technique and limited data collection in particular coastal regions [74], and improved flood forecasting [56]. Based on performance, the hydrodynamic simulation's high degree of confidence with the global tide model can be placed as input to identify inundated area during tide events. The results express the significance of gaining reliable datasets into calibration and validation processes [75]. The availability of spatial data for the study area, including DEM, bathymetry, meteorological data, and salt parcel area, also give beneficial support for the models. Although wind data do not confirm significantly in the simulation performance, there were secure connections between our models with tidal records from local station (see Figure 9).

In this study, DEMNAS (0.27 arc-second) was the highest resolution elevation data available in Cirebon and representative for the simulation and tidal inundation mapping; however, more accurate results would be achievable through higher resolutions [34,76,77] such as LiDAR-derived DTMs [56,78], and more extended tidal gauges data. Previous work by Seenath et al. [34] acquired 10-m DEM in the flood modeling component and delivered relatively higher RMS error. Additionally, smaller frequencies of simulation, i.e., 5-sec [79,80] and 6-min [81] will improve the model's stability. While there are advantages to use a hydrodynamic model for tidal flood mapping, the required computation time and the resolution of input data may limit its application in practice, especially for larger areas (compare Seenath et al. [34]). In this case, DEMNAS performs better resolution on water depth visually, but within a more extended handling time. This data was also previously exported into polyline format along with bathymetry in the preparation step. Processing simulations take a much longer time than calculating a general bathtub model. The hourly data during the 30-day period of simulation took almost 48 hours (using a standard PC with Intel I5 and 8GB of RAM).

As Cirebon salt production operates in a very traditional manner, the local salt farmers rely on daily harvests and the tidal cycle. Meanwhile, the ability to recover from disasters is far from sufficient. A comprehensive risk analysis of the salt production area is urgently needed to be better prepared to deal with the more prominent impacts of tidal floods on coastal areas. Tidal flood simulations have the potential ability to lead for better evaluations, including the potential damage loss on case-based analysis. This data will enable farmers and stakeholders to better respond to future hazards and to build capacity to improve the quality of livelihoods in the tidally flooded areas.

#### **6. Conclusions and Future Works**

This study has developed a method to identify the tidal flood impact in different types of agriculture areas in the coastline where the tide is generally forced by local factors, using hydrodynamic models. The model simulates typical aspects of tidal flood in Cirebon coastal regions, where lunar force dominates the domestic tidal properties during salt production periods. The method allows for critical identification points of flooding in the simulation, rather than using sets of scenarios. This study also underlines the main interest in the possible analysis of marginal agriculture along the coast, where the tidal hazard may continue in the future and would make a more significant impact. Meanwhile, there is still a limited number of studies that emphasizes the exposure of tidal hazard in local traditions economic activities. Tidal flood impact mapping can be beneficial to increase awareness of salt farmers to the flood occurrences. Additionally, the uncertainty and volatility level of this type of flooding is driving the local government to put more attention, especially for countermeasure planning and efficient mitigation strategies. Using higher-resolution DEM, such as LiDAR, and echo-sounder survey data for detailed bathymetry can lead to an improvement for tidal flood assessments accuracy. Although this particular model is initiated for local-case study, it is believed that this technique can be developed at a regional scale with data limitation. Finally, this research will be more substantial to include the benefit-cost (B/C) analysis of the post tidal flood events.

**Author Contributions:** Conceptualization, Anang Widhi Nirwansyah and Boris Braun; methodology, software, validation, formal analysis and writing-original draft preparation, Anang Widhi Nirwansyah; writing-review & editing, Anang Widhi Nirwansyah and Boris Braun; supervision, Boris Braun.

**Funding:** The authors would like to acknowledge financial support provided by Indonesia Endowment Fund for Education (LPDP) within code of LPDP: PRJ-115 /LPDP.3/2017 provided to Anang Widhi Nirwansyah, Institute of Geography (University of Cologne) for publication funding and Universitas Muhammadiyah Purwokerto, Indonesia.

**Acknowledgments:** This project has been part of doctoral research at the Institute of Geography, University of Cologne, Germany. The authors express their gratefulness to the anonymous reviewers for their valuable advice. Part of the data used for this paper, such as bathymetry, DEM, salt parcel on scale 1: 15,000, administrative boundary was made possible through the support of Indonesian National Geospatial Agency (BIG), and Indonesia Ministry of Marine and Fisheries Affairs (KKP). Finally, we would like to thank Junika A. Fathonah (UNDIP) for helping us during MIKE simulation, and salt farmers in Cirebon for valuable information during fieldwork.

**Conflicts of Interest:** The authors declare no conflicts of interest and the founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

#### **References**

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Review* **UAV-Based Structural Damage Mapping: A Review**

#### **Norman Kerle 1,\*, Francesco Nex 1, Markus Gerke 2, Diogo Duarte 3,4 and Anand Vetrivel 5**

Received: 22 November 2019; Accepted: 23 December 2019; Published: 26 December 2019

**Abstract:** Structural disaster damage detection and characterization is one of the oldest remote sensing challenges, and the utility of virtually every type of active and passive sensor deployed on various air- and spaceborne platforms has been assessed. The proliferation and growing sophistication of unmanned aerial vehicles (UAVs) in recent years has opened up many new opportunities for damage mapping, due to the high spatial resolution, the resulting stereo images and derivatives, and the flexibility of the platform. This study provides a comprehensive review of how UAV-based damage mapping has evolved from providing simple descriptive overviews of a disaster science, to more sophisticated texture and segmentation-based approaches, and finally to studies using advanced deep learning approaches, as well as multi-temporal and multi-perspective imagery to provide comprehensive damage descriptions. The paper further reviews studies on the utility of the developed mapping strategies and image processing pipelines for first responders, focusing especially on outcomes of two recent European research projects, RECONASS (Reconstruction and Recovery Planning: Rapid and Continuously Updated Construction Damage, and Related Needs Assessment) and INACHUS (Technological and Methodological Solutions for Integrated Wide Area Situation Awareness and Survivor Localization to Support Search and Rescue Teams). Finally, recent and emerging developments are reviewed, such as recent improvements in machine learning, increasing mapping autonomy, damage mapping in interior, GPS-denied environments, the utility of UAVs for infrastructure mapping and maintenance, as well as the emergence of UAVs with robotic abilities.

**Keywords:** drone; computer vision; point clouds; machine learning; CNN; GAN; first responder; RECONASS; INACHUS

#### **1. Introduction**

#### *1.1. Structural Damage Mapping with Remote Sensing*

The first documented systematic post-disaster damage assessment attempt with remote sensing technology dates back to 1906, when parts of earthquake-affected San Francisco were mapped with a 20 kg camera that was raised on a series of kites some 800 m above the disaster scene [1]. This makes damage mapping one of the oldest applications in the remote sensing domain, but also one of the few that continues to elude robust operational solutions, and which remains a subject of active research. Since the early pioneering days, nearly every type of active and passive sensor has been mounted on

airborne platforms that range from tethered to autonomous or piloted, as well as satellites operating in different orbital or network configurations, to attempt increasingly automated damage detection [2,3]. However, despite more than a century of research and tremendous technological developments both on the hardware and the computing side, operational image-based damage mapping, such as through the International Charter "Space and Major Disasters" or the Copernicus Emergency Management Service (EMS), continues to be a largely manual exercise (e.g., [4,5]).

Charter and EMS activations center on a particularly challenging type of damage mapping. Both need to respond to a wide range of natural and anthropogenic disaster types, and the first maps are expected to be available within hours of image acquisition, while the particular damage patterns and their recognition are subject to a number of variables. Building typologies, spatial configurations, and construction materials differ, and recognizable damage indicators are strongly dependent on the type of hazard and its magnitude. Image type, in terms of spatial and spectral characteristics, as well as incident angle, but also environmental conditions such as haze or cloud cover, differ enormously, further challenging the development of generic and widely applicable damage detection algorithms. Satellite-based damage mapping has the additional disadvantage that damage that may be quite variably expressed on each of the building's facades, its roof, as well as its interior, is largely reduced to a single dimension, the quasi-vertical perspective that centers on the roof. Damage detection in reality is then supported by the use of proxies, such as evidence of nearby debris or damage clues associated with particular shadow signatures [6,7]. There have been some notable successes in satellite-based damage mapping, especially related to cases where radar data have an advantage, in particular interferometric [8] and polarimetric synthetic aperture radar [9]. Where damage patterns are structurally characteristic, such as foundation walls remaining after the 2011 Tohoku (Japan) tsunami, simple backscatter intensity has also been used to detect damage [10]. Increasingly advanced machine learning algorithms, including convolutional neural networks (CNN), are used to detect different forms of building damage with radar data [11].

Efforts to process optical satellite data for rapid damage mapping are also moving in the machine learning direction. This includes methods based on artificial neural networks [12], and increasingly also CNN [13–16]. Studies vary in terms of mapping ambition, with many only aiming at a binary classification (damage/no damage; [12]), and there is no evidence yet of emerging methodologies being used operationally. However, the recently released xBD satellite dataset containing more than 700,000 building damage labels and corresponding to 8 different disaster types [17] will help in developing and benchmarking novel methodologies.

#### *1.2. Scope of the Review*

Automated satellite-based damage mapping has thus shown limited progress, at least in terms of versatile methodologies that can readily map structural damage caused by different event types in diverse environments. At the same time, the proliferation and rapidly growing maturity of unmanned aerial vehicles (UAVs/drones) in recent years has created vast new prospects for rapid and detailed structural damage assessment, which are the focus of this review. We do not consider historical, mainly military systems, such as unpiloted reconnaissance aircraft that date back to World War II. Rather, we focus on the suite of platforms that evolved from remote-controlled (mainly hobbyist) planes and helicopters, with the first documented scientific studies on UAV-based disaster response dating back to about 2005 [18]. The review also does not include non-structural damage assessment, such as studies on crop or forest damage. It also does not cover issues of UAV communication (e.g., use of UAVs to create ad hoc communication networks over disaster areas), nor studies on drone network or scheduling optimization. For both good reviews already exist (e.g., [19,20]).

The review includes peer-reviewed publications indexed in Scopus and Web of Science, focusing on research on automated damage detection rather than provision of data for visual assessment, and is not meant to be exhaustive. While the topic is a niche within the remote sensing domain and the amount of studies remains relatively small, a number of application papers without significant novelty exist, which are excluded here. The article is built on a recent conference contribution [21], though the focus of that paper on the results of two European research projects is expanded here to a comprehensive review study. In addition to tracing relevant technical and methodological developments in damage detection, we synthesize the current state of the art and evaluate current and emerging research directions. In addition, we assess the actual usability and practical value of emerging methods for operational damage mapping, including for local mapping by first responders. In the following section relevant publications on the use of UAVs for structural damage mapping are reviewed, sorted by increasing technical sophistication, and a summary is provided in Table 1.

#### **2. UAV-Based Damage Mapping**

#### *2.1. Scene Reconnaissance and Simple Imaging*

The principal advantage of a UAV in a disaster situation is its vantage point, a flexible position that can provide both synoptic and detailed views of a potentially complex scene, as well as overcome access limitations. Early studies thus focused on scene imaging, aiding disaster responders by supplying a relatively low-cost aerial perspective [22]. Taking advantage of increasingly efficient structure from motion (SfM) and 3D reconstruction concepts emerging at the time (e.g., [23]), in some early studies data were already processed to derive georeferenced images [24], terrain information/digital elevation models (DEM) [18], or orthophotos [25]/orthomosaics [26]. In cases without a full processing pipeline and where no suitable DEM data existed, pseudo-orthorectified images (assuming constant terrain height) were created. Instead of still images also video data were transmitted in real time to allow visual damage inspection [27].

In the years following the initial studies, little methodological progress in damage mapping was made, despite advances in off-the-shelf UAV systems, or the emergence of ArduPilot in 2007 for improved UAV flight stability, or Pix4D(Pix4D, Switzerland) in 2011 for easier photogrammetric image processing. A range of studies appeared that essentially still focused on image provision or simple photogrammetric processing, using remote-controlled helicopter systems [28], multi-copters [29–32], or fixed-wing UAVs [33–35].

#### *2.2. Texture- and Segmentation-Based Methods*

Initial attempts to extract damage information automatically from UAV data were based on segmentation- and texture-based approaches, using mono-temporal imagery. Fernandez Galarreta et al. [36] processed UAV imagery of an 2012 Emilia Romagna (Italy) earthquake site into detailed 3D models. The work adapted and expanded earlier approaches developed for the airborne (piloted) Pictometry system that yields similar oblique, overlapping, and multi-perspective imagery. Also, those images had been photogrammetrically processed [37], and used for structural damage assessment [38,39]. The analysis of [36] focused on geometric damage indicators such as slanted walls or deformed roofs, as well as presence of debris piles (Figure 1). In addition, object-based image analysis (OBIA) was carried out on the images to extract damage features such as cracks or holes, but also identification of those damage features intersecting with apparent load-carrying structural elements. A similar OBIA strategy was used by [40] to identify damage in Mianzhu city, affected by the 2008 Wenchuan earthquake.

**Figure 1.** Damages identified from unmanned aerial vehicle (UAV)-derived point clouds and from object-based image analysis (OBIA) processing. (**a**) Inclination in walls, (**b**) openings (turquois), cracks (magenta), and damage crossing beams, (**c**,**d**) detailed point cloud and segment orientation angles [adapted from 36].

#### *2.3. Conventional Classifiers*

The work of Fernandez Galarreta et al. demonstrated the significance of geometric information in damage detection, in particular of openings in roofs and façades. Vetrivel et al. [41] advanced the work by developing a method to isolate individual buildings from a detailed image-derived point cloud covering a neighborhood of Mirabello (Italy) comprising nearly 100 buildings. Each of those was then subjected to a search for openings attributable to seismic damage, such as partial roof collapses or holes in the façades, a focus similar to [42]. The gaps were identified based on Gabor wavelets as well as histogram of gradient (HoG) orientation features. Two basic machine learning algorithms, Support Vector Machine (SVM) and Random Forest (RF), were used to identify damaged regions based on the radiometric descriptors, with a success rate of approximately 95%. However, the work also illustrated how the segmentation of point clouds is frequently hindered by artefacts and data gaps. In [43], an approach was developed to overcome this problem: after projecting the initial point cloud-derived 3D segments into image space, a subsequent segmentation using both geometric and radiometric features yielded more accurate and complete building segments.

Manned system with 5 cameras comparable to Pictometry. Though not a UAV category, some relevant Pictometry-based studies are included in the review. Usehigh-resolution airborne data, with no specific platform being indicated.

1

The work in [41] work also showed the limitations of HoG and Gabor filters in the classification of complex scenes, and of global feature representations on general. The latter cause problems when scene and image characteristics vary, which is typically the case between different disaster areas or in multi-temporal assessments. The work described in [68] moved towards descriptors that are more generalizable and invariant to image characteristics. The method was built on the Visual Bag of Words approach and focused on the detection of rubble, debris piles, and severe spalling. The method performed well on individual UAVs and also Pictometry data sets of Mirabello (Italy) and Port-au-Prince (Haiti), respectively, but also on a dataset that combined the two airborne datasets with transverse street-level images. The limitation of the method is that it is grid-based and can only identify general damage patches, i.e., grid cells affected by one or more of the damage types considered, a limitation also evident in the study of [50], who used RF on superpixels. A detailed localization and characterization (size, shape, etc.) of damages of a specific type would be preferable, though this will come at the cost of increased processing time.

#### *2.4. Advanced Machine Learning and the Emergence of CNN*

Image classification used for damage mapping increasingly made use of machine learning, in particular SVM and RF [41,54,69] or different boosting algorithms, such as AdaBoost [38] or XGBoost [58], and moving towards more advanced scene understanding and semantic processing. However, the features used were typically hand-crafted (such as HoG or Gabor, or other point feature descriptors related to spectral, textural, and geometrical properties [54]), and emerging work had shown that in deep learning approaches CNN could actually learn features and their representation directly from the image pixel values [70]. Thus, the damage detection work proceeded in this direction, hypothesizing that image classification would benefit from the micropropagation of 3D point cloud features. The work described in [55] applied a multiple-kernel-learning framework on several sets of diverse aerial images, and showed that combining the radiometric and geometric information yields higher classification accuracies. The processing was based on Simple Linear Iterative Clustering (SLIC) superpixels, meaning that damage was again only identified in patches, though those were labelled with specific prediction scores. Song et al. [66] also worked with SLIC superpixels, though unlike in [55] where they had formed the basis for the ML analysis, here first a CNN-based damage detection was carried out directly on the image, and the SLIC segments were then used in combination with mathematical morphology to refine the results. In [67] a similar approach was taken, except that instead of SLIC a multi-resolution segmentation was carried out, to allow features naturally occurring at different spatial scales to be used effectively. The CNN approach developed in [55] was also used by Cusicanqui et al. [71], who reasoned that video data are often available before suitable still photographs (e.g., acquired by police or the media). In the study it was thus tested whether 3D reconstructions based on video data could offer similar support, and it was indeed shown that a binary damage classification based on deep learning applied to SLIC superpixels and the 3D models led to results comparable to those based on still photographs.

The particular significance of the work in [55] for disaster response and search and rescue was that the method demonstrated significant transferability, which has become a frequent focus in recent literature. A model trained with a sufficient number of samples (e.g., trained before an actual event) performed well when then applied to a new disaster scene, supporting a rapid analysis without the need for extensive retraining. This approach can help to overcome the traditional limitation of CNN, i.e., their need for a large amount of labelled training data. A different approach was taken by Li et al. [62], who used a convolutional autoencoder (CAE) that was trained using unlabelled post-disaster imagery based on SLIC superpixels, with results being finetuned by a CNN classifier. In follow-up work [63] the authors in addition employed a range of data augmentation methods, such as data blurring or rotating, to enlarge the number of samples. The resulting pre-training improved the overall damage detection accuracy by 10%.

Disaster scenarios are frequently characterized by imperfect image data availability, and a rapid response effort has to make do with what exists. In this respect it is valuable to be able to incorporate images of different types and scales into the training model. Duarte et al. [56] trained a CNN with different types of aerial imagery to classify post-disaster satellite data of Port-au-Prince. Although information coming from the different image resolutions evidently improved the model and classification accuracy, the approach still failed to capture smaller damage features. The work also focused on determining the effect of multi-scale information on the CNN activation layers as a proxy for improved damage recognition, while not allowing a detailed assessment of where the classification improvement originated in terms of false positives and negatives, or specific damage types. Later work focused on multi-resolution feature fusion and its effect on building damage classification [57]. It showed that such a fusion is useful and can improve the overall accuracy, though it still failed to show which specific damage types are identified, and how well they are captured.

Earlier work had shown how highly variable the expression of structural damage is in vertical and oblique data [6]. The former essentially only considers the damage expressed in the roof, and in addition makes use of proxies such as debris piles for specific shadow configurations [7]. Significant additional information is also encoded in the façade information, as already explained in Section 2.1. However, the OBIA-based approach used for example in [36] tends towards overfitting and lacks the efficiency and transferability of deep learning. While a focus on façades is appealing, their actual delineation in imagery poses its own challenges, especially when considering aspects such as occlusion or environmental effects such as shadows (Figure 2). The work described in [51] thus focused on developing an efficient method to extract façades that were subsequently assessed for damage using CNN. The approach made use of a point cloud calculated from vertical imagery acquired in an initial UAV survey. From the sparse point cloud, the building roofs were segmented and the building façades hypothesized, which in turn was used to extract the actual façades from oblique UAV images. The patch-based damage classification had an overall accuracy of approximately 80%, though the work also demonstrated the significant challenge of damage identification on façades, due to architectural complexities and associated diverse shadow patterns, but also occlusion (by external features such as vegetation, or internal ones such as balconies).

**Figure 2.** Typical problems for image processing posed by shadow and occlusion [51].

It stands to reason that some ambiguities can be resolved by analysing multi-perspective data (views of a given façade from different angles that go beyond regular stereoscopic overlap), but also by incorporating multi-temporal data where available. The majority of the studies described above only used post-disaster imagery. However, in the last few years the availability of high spatial resolution pre-event reference imagery has been growing rapidly. This has led to additional methodological developments that built on the segmentation- and texture-based damage detection described above, extending them into a multi-temporal framework. Vetrivel et al. [47] used pre- and post-earthquake data of L'Aquila (Italy) and focused on the identification of 3D segments missing in the post-disaster data as an indicator of damage. Both voxel- and segment-based approaches were tested, and finally a composite segmentation method that subjects an integrated pre- and post-event point cloud to plane-based segmentation was chosen. Although working with conventional airborne data, in [61] those assumptions were also tested in a CNN framework, where 6 different multi-temporal approaches were compared against 3 mono-temporal ones. It was concluded that a multi-temporal approach with 3 views at each the pre- and post-event epoch performed best. Also, here smaller damage features eluded detection. However, the authors expect better results with UAV data, given that the problem of occlusion can be reduced through more flexible image acquisition.

#### *2.5. Levels of Disaster Damage Mapping*

Early efforts in disaster response with satellite imagery identified damaged areas more generally, while airborne data were used to detect specific damage proxies, usually debris piles (e.g., [72]). Especially in more recent years, overall classification accuracy and f-scores have been the most commonly used metrics to assess the efficacy of a given damage mapping method, and to judge progress within the discipline. However, this focus neglects an inherent incomparability of many of the studies produced to date, and the absence of a generally agreed upon damage scale. The introduction of the European Macroseismic Scale 1998 (EMS-98) led to a broad homogenization and alignment of efforts, by grouping structural building damage in 5 categories, D1 (negligible/slight damage)–D5 (destruction) [73]. Building on its common use in satellite-based damage detection (e.g., [74–76]), later its utility for UAV-based damage mapping was explored. For example, [38] classified building damage according to EMS-98, though recognizing the diversity and ambiguity of the observed damage patterns the study did not aim at automatic damage classification, except in cases where the 3D model clearly showed complete collapse (D5). Also, studies [36,48,77] used this scale as a basis, with [31] even adding a 6th damage level.

One consequence of the continuing challenge of image-based damage mapping is that, while D1 and D5 are comparatively easy to determine but intermediate damage stages are not, many studies have departed from the 5-level classification scheme. The work in [50] opted instead for a 4-class approach (intact, light, medium, and heavy damage), while several studies grouped damage into 3 classes. However, even within one such category damage levels/class names vary, limiting comparability. For example, Zeng et al. [40] mapped intact, damaged, and destroyed buildings, while Vetrivel et al. [47] termed the classes undamaged, lower levels of damage, and highly damaged/collapsed, and Song et al. [66] distinguished intact, semi-collapsed, and collapsed buildings, with differences in class definition going beyond semantics. However, the majority of recent studies opted for a simple binary classification, either explicitly mapping both damaged and undamaged structures (e.g., [12,49,55,65,67]), or only mapping damage in general in a single class [51,61]. In addition, there are studies that focused on the identification of specific damage types, such as holes in the roof [41,42], or dislocated roof tiles and cracks along walls [36]. Others mixed damage and proxy classes, such as [62], who mapped damaged and undamaged structures, but also debris as a separate class. Creative choice of class names is further hindering a comparison between different studies. Li et al. [63] used the classes mildly damaged and ruins, while Xu et al. [54] mapped categories including roof, ground, debris, and small objects. The difficulty of image-based damage mapping has led to a focus on severe damage classes (D4-5), making studies such as [42] that expressly focus on lesser damage (D2-3) an exception. Approaches

based on deep learning are particularly suited for binary classification, which is another reason why in the interest of automation only a single damage class is now frequently considered.

#### *2.6. The Special Case of Infrastructure Damage Mapping*

The focus of this review is on structural building damage. However, one of the fastest growing UAV application areas in recent years is infrastructure monitoring and detection of damage indicators related to wear and degradation, such as of roads, bridges, or tunnels. The lines between disciplines have blurred, with studies such as by Dominici et al. [32] addressing both regular structures and infrastructure. Furthermore, from a methodological perspective studies focusing on crack or spalling assessment along bridges or tunnels are also relevant for the disaster damage mapping community, and damage to infrastructure caused by disaster events naturally also falls under the scope of this review. For this reason, papers marking key developments in infrastructure monitoring and damage mapping are briefly reviewed here.

A recent review by Dorafshan and Maguire [52] provides an overview of the specific challenges of bridge inspection and maintenance, and how UAVs, both with active and passive sensors, are starting to become a commonly used tool. In an early study by Whang et al. [22], a UAV with two coaxial rotors was developed to perform somewhat autonomous bridge inspection, within limits even in GPS-denied areas beneath the bridge. In addition, the system was able to place a small autonomous rover on the bridge using ultrasonic localization, and which provided images for damage inspection. However, few details about the actual methods and system performance are provided in the paper. The authors of [44] focused on the detection of small fatigue cracks on bridges, assessing the value of active illumination, and carrying out controlled laboratory experiments to determine detection limits and optimal mapping approaches.

Increasingly, the focus has been on image- or laser-based 3D reconstruction of the bridge or tunnel in question, as a basis for visual or automated damage identification. In [78] the accuracy and thus utility of such 3D models was assessed, and [45] also assessed how well complex bridge structures can be reconstructed with SfM methods, in addition attempting 3D volume calculations or major spalling instances. The work of [79] expressly focused on seismic damage detection on bridges, also using UAV-based 3D reconstructions, though here starting with pre-event Building Information Modelling (BIM) data that were updated with the detected damage. Akbar et al. [46] addressed structural health monitoring (SHM) of tall structures, focusing on comprehensive 3D model creation through speeded up robust features (SURF), and on the detection of simulated damage features on large concrete slabs, though providing little detail on the actual damage detection algorithm.

Deep learning with CNN is also being used in SHM. In [53] an AlexNet network was trained to detect small cracks in concrete walls, reporting accuracies of nearly 95%, and also testing network transferability. Comparable accuracies were reported by Liang [64], who in addition also tested GoogleNet and VGG-16 networks to detect earthquake damage on a bridge.

#### **3. Damage Product and System Usability**

Post-disaster damage mapping serves a specific purpose, i.e., providing timely, accurate, and actionable information to a range of stakeholders. Those include civil protection agencies planning emergency response actions, but also incident commanders and first responders operating at the actual disaster site. One of the consequences of the growing availability of UAV technology is a declining need to rely on formal protocols such as the Charter or EMS, and instead allowing actual site-based damage mapping. It is thus surprising that the usability of data acquisition pipelines (including planning tools, hardware components, and data processing routines), but also of resulting damage mapping products, has scarcely been considered in the literature reviewed in this paper. This section briefly introduces two recent research projects with a strong focus on UAV-based structural damage assessment, and from which a number of publications reviewed in this paper emerged. In

these projects also a range of different end users participated, and their evaluation of the developed damage mapping procedures is also summarized.

#### *3.1. Damage Detection in Two European Research Projects*

RECONASS (Reconstruction and Recovery Planning: Rapid and Continuously Updated Construction Damage, and Related Needs Assessment; www.reconass.eu) and INACHUS (Technological and Methodological Solutions for Integrated Wide Area Situation Awareness and Survivor Localization to Support Search and Rescue Teams; www.inachus.eu) were research projects funded through the 7th Framework of the European Union, and which ran with some overlap from 2013 until the end of 2018. The focus of RECONASS was to create a system for monitoring and damage assessment for individual high-value buildings, based on a range of internally installed sensors that included accelerometers, inclinometers, and position tags, with data getting processed in a finite element structural stability model to determine damages caused by seismic activity or by either interior or exterior explosions. UAV-based 3D reconstruction of the building exterior and detailed damage mapping were carried out to patch data gaps caused by failed sensor nodes, as well as to validate model outputs. The progressively developed methods were tested in a series of experiments, culminating in a pilot where a 3-story reinforced concrete building was first subjected to an explosion of 400 kg TNT placed 13 m away, and later by a 15 kg charge detonated within the structure itself. End users, including the German Federal Agency for Technical Relief (THW), were present to assess the utility of the system.

The purpose of INACHUS was to assist disaster response and urban search and rescue forces by providing early and increasingly detailed information on damage hotspots and the likely location of survivors. Different UAV platforms, but also ground-based and portable laser scanning instruments, were used to map a damaged structure. One research focus was on scene reconstruction and damage mapping based on optical imagery from a low-cost UAV. The French remote sensing lab ONERA also deployed various larger UAVs that carried different laser scanners, in part with proprietary data processing solutions. The major pilots were also assessed by a group of end users.

#### *3.2. Tests with End Users in Two European Research Projects*

Both RECONASS and INACHUS included a number of pilot experiments, where first individual components or sets thereof, and later the entire systems were tested under relatively realistic conditions. For the explosion experiments in Sweden data were acquired using an Aibot X6 Hexacopter carrying a Canon D600 camera with a Voigtländer 20 mm lens. In addition to reference data, images were acquired after both the exterior and the interior blasts, with a ground sampling distance (GSD) of approximately 1.5 cm. From those images, detailed 3D point clouds were calculated and analyzed. The data proved suitable to identify damage-related openings, such as infill walls damaged or blown-out by the blasts, as well as cracks and debris. Additionally, subtle façade deformations could be detected and quantified (Figure 3), both using only the post-detonation point cloud, as well as in a comparison with pre-event reference data. It was also shown how a BIM model of the structure could be automatically updated, both to visualize and catalogue detailed damage information. THW deployed a LEICA TM30 total station to survey the structure from 4 reference points, using 16 prisms mounted on the structure. While the total station has the advantage that a structure can be continuously monitored for minute deformations—critical when rescue personnel operates near or within weakened structured—the UAV-derived data provided damage data of comparable quality, with greater flexibility and lower cost, including the roof that ground-based surveys cannot see, and potentially operated from a safer distance. The building was further surveyed by a Riegl VZ400 terrestrial laser scanner (TLS), which also confirmed the high quality of the UAV-derived 3D models.

**Figure 3.** UAV-derived point clouds of reinforced concrete structure with brick in-fill walls subjected to exterior and interior detonations. Openings, cracks, and debris piles, as well as subtle deformation in the façades were automatically detected.

Four INACHUS pilot experiments were conducted at 4 different sites in France and Germany, and included buildings in the process of being demolished, as well as an urban search and rescue training site (Training Base Weeze in Germany). In response to criticism by end users in RECONASS as to the high cost of the Aibot UAV (ca. 40,000 Euro), in INACHUS low-cost DJI drones (Phantom 4 and Mavic Pro) were used. Following the research directions described in Sections 2.3 and 2.4, the work focused less on simple scene reconstruction, but on integration with other spatial data, as well as advanced data analysis, including with CNN. For each of the pilots, the building in question was also surveyed by ONERA using different UAV-borne laser instruments, as well as with a TLS, to detect the respective strengths of the individual systems. The initial experiments with UAV-based laser scanners failed. First a Riegl VZ-1000 instrument (weight of about 10 kg) was deployed on a Yamaha RMAX helicopter (weight > 60 kg), though the acquired data suffered from artefacts and were not useful. Also, data acquired with a Velodyne HDL32 (weight of only 1.3 kg) deployed on a VARIO BENZIN helicopter (weight just under 10 kg) proved unusable for damage detection, owing to the very unstable platform. For the final pilot, a high quality Riegl VUX-1 was mounted on a stable DJI Matrice 600 hexacopter platform. The data were excellent, though the combined system is also very costly (>80,000 EUR) and requires expert knowledge for flight planning and execution, as well as data processing. The mapping with optical data focused on using data acquired with the built-in cameras of the Phantom 4 and Mavic Pro (costs of < 2000 Euro), and advanced along the computer vision and machine learning trajectory described earlier. The 3D data obtained from the optical imagery were of comparable quality to the VUX-1 data while also providing native color information, better spatial detail, and full coverage also of façades (Figure 4). The expectation that the airborne laser data would patch the one principal weakness of photogrammetry, the inability to map dark interior spaces through openings (as a means of possibly locating trapped survivors), was also not met. The data on openings and connected interior spaces were primarily delivered from the tripod-mounted ground-based laser scanner, though here the limited flexibility and occlusion by the building's structural elements also prevented a complete mapping of openings.

While commercial UAVs by DJI and other makers have clearly reached high levels of cost-benefit, stability, and reliability, most are also not designed to be survey-grade instruments working in real time. For rapid search and rescue support it is vital to provide usable information quickly. For that reason, in INACHUS a procedure was developed to process the data with minimal delay. Working with the ability of the Mavic Pro to stream images during flight, a procedure was built that (i) downloads images right after acquisition, (ii) builds a progressively extended sparse 3D model of the scene using established SfM methods, (iii) applies CNN to detect damage, and (iv) orthorectifies the images using

the 3D model. By the time the UAV lands after a maximum flight duration of about 25 min, all processing is done and the damage map available. A smart phone app was also built that allows this procedure to be executed together with a standard laptop (Figure 5). Details about the app and data processing workflow can be found in [59], while more information about the optimized CNN that was made available on GitHub can be found in.

**Figure 4.** Point cloud representation of an INACHUS pilot structure in Lyon, France, calculated from optical imagery acquired with a low-cost commercial drone (Phantom 4, DJI), showing damage detected through machine learning (red).

**Figure 5.** Workflow of the app developed for near real-time damage mapping. Images are streamed to a laptop computer and processed immediately after acquisition. A convolutional neural network (CNN)-based damage detection algorithm is applied, and a progressively built sparse 3D model is used to orthorectify them. By the time the UAV lands, an orthomosaic displaying the damage is finished (adapted from [59]).

#### *3.3. Validation*

At every pilot, different end users were present and undertook a detailed assessment of every tool produced and tested. The RECONASS system was evaluated by THW at the pilot site, and more extensively in a dedicated workshop at ISCRAM 2017 by a total of 11 specialist end users, representing both governmental and non-governmental emergency response organizations, as well as organizations involved in the creation of damage maps. It was concluded that the UAV-based element met all previously established user requirements, principally the detection of all externally expressed damage types and their annotation both on imagery but also a 3D model and a BIM, as well as the provision of 3D volume calculations, all in GIS-ready format. The final system received a maximum score of 10/10.

At the final INACHUS pilot that took place in Roquebillière, France, in November 2018 a total of 25 end users from 8 countries participated, representing USAR teams and other civil protection organizations. They followed individual demonstrations of all technical tools developed and graded them. Of all hard- and software or procedure solutions developed in INACHUS, the 3D mapping and damage detection with a light-weight commercial UAV scored highest (overall 4.5 out of 5). The high score does not so much represent a high level of technical sophistication, but rather the simplicity, both in terms of off-the-shelf hardware and an automated flight planning and damage mapping routine. The end users especially appreciated the simple, low-cost approach that provided accurate and useful information in near-real time, without the need for a highly specialized operator.

#### *3.4. Limitations*

Despite the positive evaluations, the end user assessment also revealed limitations of the developed damage mapping solution. Legal restrictions of drone deployment continue to pose challenges, though problems are less severe for lighter platforms, and in addition first responder and civil protection organizations tend to operate under different legal frameworks. A clear disadvantage of small multi-copter UAV platforms is their comparatively small operating range and flight duration. The limited spatial scope of RECONASS and INACHUS matched their abilities well, but damage assessment over larger affected areas requires different solutions. Off-the-shelf UAVs come equipped with high quality optical cameras, though the computer vision processing to generate 3D point clouds fails for dark image patches such as shadow or smaller building openings. For this reason, openings and possible survival spaces in the pilot structures could not be mapped, and here active sensors have a clear advantage. Commercial UAVs also tend to be closed and largely proprietary systems, meaning that it is not easily possible, if at all, to exchange or add sensors, or to install processing units such as a DJI Manifold (China) or NVIDIA Jetson TX2 (USA) to push more autonomy in onboard image processing or dynamic flight path adjustment onto the drone. Several of these limitations are the focus of ongoing research, as explained in the following section.

#### **4. Outlook and New Developments**

The literature reviewed in this paper mirrors a rapidly developing discipline that in only a few years moved from largely descriptive imaging of disaster scenes to fully automated analysis procedures that build on state-of-the-art methods originating, in particular, in the computer science domain. At the same time, limits persist in hard- and software, in operational damage mapping procedures, but also in the conceptual basis of how images can be related to the actual meaning and significance of damage, which are addressed in this section.

#### *4.1. Improvements in Machine Learning*

For all the sophistication of machine leaning approaches to recognize patterns and features, some open questions persist. The black box nature of deep learning approaches means that the specific effect of certain training labels remains unclear, challenging efforts to optimize the training efficiency for specific damage features. Training to map only specific indicators such as cracks or object dislocations is thus challenging, compounded by the scarcity of large training samples for individual damage features. Also, solutions developed to date still tend to be patch/grid-based, highlighting damage in general, but not specific features. This, however, is highly scale dependent, with high resolution image data, for example, also yielding small superpixels that allow precise damage identification [50].

Work such as in [56,57] tends to focus on activation layers that indicate the presence and approximate position of damage (Figure 6), rather than the creation of actual damage maps. From a user perspective, more clarity on the specific damage type, but also more precise location, shape, and size, would be preferable. In addition, the nature of CNN-based studies prevents insights into how specifically a network with superior overall accuracy performs in terms of reducing false positives or negatives.

**Figure 6.** Examples of building damage detection via CNN activation layers, using aerial pre- and post-earthquake façade images. Bright activation colors show damage hotspots (adapted from [61]).

To overcome the problem of the large number of training samples needed in CNN analysis, recent work has shown how Generative Adversarial Networks (GAN) can effectively enlarge sample databases, which has already been shown to benefit the identification of damage of road furniture [60]. GAN seem to be particularly useful in anomaly detection [80], where training does not focus on a potentially large number of specific damage features or indicators, but rather where a comprehensive understanding of normal, undamaged scenes is created, based on which anomalies such as damage are identified. GAN have been mainly used in applications with smaller variabilities than are typical for urban scenes (i.e., indoor environments with fixed cameras). Their use in urban scenes is, therefore, an additional challenge that could be compensated by only using very large and comprehensive datasets of undamaged scenes, to prevent the generation of many false positives.

#### *4.2. Mapping Autonomy*

Traditional UAV surveys were based on pre-defined flight plans or manual piloting supported by video streams from the instrument, with data getting processed after landing of the aircraft, or through pipelines such as described in [59]. A more ideal scenario would be for the UAV to carry out an initial, for example vertical, survey over a pre-defined area, identify hotspot and damage candidate areas based on limited real-time processing, followed by a more detailed and multi-perspective survey of those marked areas. The work in [51] showed how data from an initial coarse vertical survey can be used to guide a more local assessment. Such a procedure can be implemented based on streamed data that are processed in near-real time, and adjusted flight path instructions uploaded. Alternatively, data can be processed on the UAV itself. Work described in [81,82] showed how even microdrones can perform analysis based on deep neural networks to facilitate autonomous navigation. UAVs with greater payloads have been fitted with more powerful computing units, such as NVIDIA Jetson TX2, which are capable of facilitating advanced real-time object tracking [83] or image segmentation [84].

In follow-up work to INACHUS, H2020 project PANOPTIS (Development of a Decision Support System for increasing the Resilience of Transportation Infrastructure; www.panoptis.eu) focuses on road surface and road corridor damage assessment to detect signs of gradual wear and decay, as well as the ability to respond rapidly to a disaster situation. This is done with a hybrid UAV platform (DeltaQuad from Vertical Solutions) that allows both corridor mapping of a fixed-wing platform and hovering for detailed mapping. Also, here a Jetson TX2 will be used to advance data processing on the drone itself, for both navigation and damage detection.

#### *4.3. Indoor Mapping*

UAVs have brought structural damage mapping within touching distance of buildings. Nevertheless, critical damage evidence is frequently hidden from sight, e.g., where internal load-carrying structures are compromised. In addition, damage assessment, such as defined in INACHUS, also includes support for first responders in the search for victims or survivors trapped internally, though different cavity mapping strategies only had limited success. Even with a TLS, interior cavities with connections to the outside could only be detected to a limited extent (Figure 7).

**Figure 7.** Voids within the photogrammetric model shown in Figure 4, obtained with a terrestrial laser scanning system. (**a**) Estimated size of open spaces observed through openings, (**b**) distance of voids to the edge of the building.

Recent work has demonstrated UAVs operating increasingly autonomously and effectively in interior, largely-GPS-denied spaces [85]. There has been a surge in research on UAV-based indoor mapping, both with single platforms and swarms. Most make use of visual SLAM to map their GPS-denied environment (e.g., [86,87]), or focusing on continuity mapping when transiting between outdoor and indoor places [88]. Others have experimented with localizing via sensors such as ultrasound [89], and the works cited in 4.2 on autonomous navigation and mapping are also relevant here. One element of improved indoor 3D reconstruction and damage mapping will be a more effective use of artificial lighting that, for example, improved the detection of small cracks in [44]. Another line of research has focused on the engineering of UAV platforms that can change shape to facilitate their entering and operating in tight spaces [90].

The damage detection work of INACHUS will also be advanced in indoor environments with H2020 project INGENIOUS (The First Responder of the Future: A Next Generation Integrated Toolkit for Collaborative Response, increasing protection and augmenting operational capacity; www.ingenious-firstresponders.eu). The focus will be on the use of drone swarms for indoor mapping to support first responders in unknown and potentially dark, smoke-filled, and hazardous indoor settings, using UAV platforms of different sizes and with different sensor load and ability, with focus on collaboration and optimization.

#### *4.4. The Age of Drones with Robotic Abilities*

UAVs tend to be fragile, susceptible to wind, and, thanks to inexpensive GPS and IMU components subject to positional inaccuracies, best operated away from structures. However, better platform control, use of collision avoidance through use of sensors or depth sensing, as well as progress in robotics and mechatronics have resulted in novel research directions. For example, in the age of aging infrastructure efforts have been spreading towards UAV-supported maintenance. This implies a number of challenges. Infrastructure is diverse and includes complicated indoor spaces such as chimneys [91], but also roads, tunnels, and bridges. Solutions are emerging to carry out day-to-day monitoring to detect defects or signs of decay, but also damage after a disaster event or accident (e.g., [92]). Such works increasingly extend into another emerging line of development, blending UAV-based abilities with robotics and mechatronics solutions. Here, UAVs are not only used to map and model infrastructure spaces, but also to carry actuator arms to place sensors for in-situ measurements [93,94], interact with objects [95,96], perform physical tests [97], or to carry out limited repairs.

#### **5. Conclusions**

Structural damage mapping with remote sensing has been a continuous research problem for decades, and for rapid operational disaster response, such as through the Charter or Copernicus EMS, reliable automated methods continue to be lacking. However, substantial progress has been made in the last decade that resulted primarily in rapid developments in UAV technology, computer vision, and in advanced image data processing with machine learning, in particular deep learning with CNN, all of which was assessed in this review. This includes a detailed analysis of the progress in image-based damage mapping that has moved from providing largely descriptive overview imagery to automated scene mapping with advanced machine learning.

The paper has shown how image-derived 3D point clouds allow a highly detailed and accurate scene reconstruction, and how the coupling of the geometric information with the original image information allows very advanced feature recognition. Classifier training is also starting to overcome the challenge of, in particular, CNN-based methods, requiring millions of training samples. The development of unsupervised CNN approaches (such as Auto-encoders) or Generative Adversarial Networks (GAN) could represent a step forward in this direction. Newer approaches are improving the efficiency, but also the transferability of classifiers, critical to be able to respond quickly to a disaster event. Comprehensive tests with first responders and urban search and rescue personnel showed that, in particular, solutions with light-weight off-the-shelf drones strike a very good compromise of high information quality and ready usability.

Developments continue at a rapid pace, with significant research efforts now being focused on UAV-based mapping in indoor settings, on UAVs also being equipped with mechatronic abilities to allow the deployment of additional sensors or to carry out repairs, though newer networks also allow more sophisticated and robust deep learning solutions. Nevertheless, more effort is needed to understand better the actual meaning and significance of specific damage evidence. In addition, UAVs need to become more autonomous to increase the efficiency of damage mapping operations. Finally, progress in the processing of UAV-based imagery, in particular through advanced machine learning, must eventually lead to fully automated and accurate damage mapping with optical satellite imagery.

**Author Contributions:** Conceptualization, N.K..; Methodology, N.K.; Formal Analysis, N.K., F.N., M.G., D.D., A.V.; Investigation, N.K., F.N., M.G., D.D., A.V.; Writing—Original Draft Preparation, N.K.; Writing—Review & Editing, N.K., F.N., M.G., D.D., A.V.; Visualization, N.K., F.N., M.G., D.D., A.V.; Supervision, N.K.; Funding Acquisition, N.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was funded by the EU-FP7 projects RECONASS (grant no. 312718) and INACHUS (grant no. 607522), as well as H2020 project PANOPTIS (grant no. 769129).

**Acknowledgments:** We thank Pictometry, Inc. for providing the Haiti and Italy imagery used in this study, and the DigitalGlobe Foundation (www.digitalglobefoundation.com) for providing satellite images on Italy and Ecuador. We also appreciate the comments made by three anonymous reviewers.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*ISPRS International Journal of Geo-Information* Editorial Office E-mail<EMAIL_ADDRESS>www.mdpi.com/journal/ijgi

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18