Unmanned Aerial Vehicle Systems for Remote Estimation of Flooded Areas Based on Complex Image Processing

Floods are natural disasters which cause the most economic damage at the global level. Therefore, flood monitoring and damage estimation are very important for the population, authorities and insurance companies. The paper proposes an original solution, based on a hybrid network and complex image processing, to this problem. As first novelty, a multilevel system, with two components, terrestrial and aerial, was proposed and designed by the authors as support for image acquisition from a delimited region. The terrestrial component contains a Ground Control Station, as a coordinator at distance, which communicates via the internet with more Ground Data Terminals, as a fixed nodes network for data acquisition and communication. The aerial component contains mobile nodes—fixed wing type UAVs. In order to evaluate flood damage, two tasks must be accomplished by the network: area coverage and image processing. The second novelty of the paper consists of texture analysis in a deep neural network, taking into account new criteria for feature selection and patch classification. Color and spatial information extracted from chromatic co-occurrence matrix and mass fractal dimension were used as well. Finally, the experimental results in a real mission demonstrate the validity of the proposed methodologies and the performances of the algorithms.


Introduction
In the repertory of natural disasters, floods often cause the greatest material damage [1]. For example, floods in 2013 constituted 31% of the economic losses resulting from natural disasters [2]. Therefore, the forecasting, prevention, detection, monitoring, and flood damage assessment are of paramount importance for public authorities and people. Because early warning is essential for limiting the loss of life and property damage, many works examine real time flood detection systems [1,3,4]. For example, the integration of several existing technologies was used in Taiwan to develop a coastal flooding warning system [3].
The problem that we are solving in this paper is the evaluation of small flooded areas in rural zones with the aid of a cheap solution based on processing of images taken by the unmanned aerial system designed by the authors. The result is necessary to evaluate the post-flood damage by the local authorities (to determine relief funds) and assurance companies (to determine payments). Even small flooded areas are investigated and classified. For this purpose two sub-problems must be solved. First is the optimal coverage of the area to be monitored, from the point of view of energy consumption (a limiting factor, especially if the UAV is electrically powered through a battery) and trajectory length. The second sub-problem is the detection and estimation in terms of flooded areas of the covered regions.
For the purpose of flood detection and monitoring, different sensors have been used: optical sensors, multi-spectral sensors and synthetic aperture radars (SARs). Satellite remote sensing imagery offers less spatial and temporal resolution than aircraft and UAVs, but a larger field of view. It was successfully used on large-scale geographic analysis on the flood overflow area. For example, images from Landsat Thematic Mapper/Enhanced Thematic Mapper Plus (TM/ETM+) sensors were used to monitor the floods near Lena River (Siberia) with a spatial resolution of 30 m and a temporal resolution of 2.6 days [5]. On the other hand satellite images in the visible and near infrared spectrum are highly dependent on cloud conditions whereas radar transmitters and receivers can be used independently of weather conditions [6]. Based on the surface water extent, measured with a microwave remote sensor (Radiometer for Earth Observation System AMSR-E and AMSR2), the authors in [7] implemented a method for detecting floods at large scale. In [8] the authors propose to combine aerial thermal data with high resolution RGB images in order to quickly and accurately identify the sub fluvial springs of a stream. Both cameras, thermal and action, are installed on a rotor platform which is able to support more weight, but has a smaller surveillance area.
Combining information from space, aerial and ground [6], an integrated system using different technologies (remote sensing, the Global Navigation Satellite System (GNSS), data transmission, and image processing) was implemented in China for monitoring and evaluating flood disasters. Because of the high cost, aircrafts use SAR only for serious and urgent flood cases. The spatial resolution is of 3 m at 9 km swath and the monitoring is in real time, independent of weather. Among the methods for flood detection, the interpretation of optical remote images is widely used and also gives the best results concerning price and accuracy.
In order to detect the flood by image analysis, three solutions usually appear in the literature: (a) use of images from satellites [9][10][11]; (b) use of images from fixed cameras on the ground [4,12,13]; and (c) use of images from aircrafts or UAVs [14].
The Advanced Land Observing Satellite Phased Array type L-band Synthetic Aperture Radar (ALOS PALSAR) satellite [11] provides multi-temporal data which maps large zones of flood via a classification into flood and non-flood areas. Based on this classification and on images taken at pre-flood and post-flood time instants, information about the flooding hazard is provided. In general the satellite applications for flood detection, like those presented in [10] from Spot-5 imagery, or in [15] by WorldView-2 satellite imagery, are based on high spatial resolution images and have the disadvantage of being high-cost solutions, hence less approachable for public use. In addition, these solutions have the disadvantage of being sensitive to weather patterns (clouds or other obscuring weather features will render them useless).
An alternative approach for monitoring flood disasters is the system of fixed cameras proposed in [12] which is based on the dynamic detection of floods via intrusions of objects in the video frames. These objects are separated by segmentation from the rest in the image.
To monitor and evaluate the flood disasters, concatenated images, created by photomosaic generation, can be useful. Thus, the gaps or duplications of flooded regions, in different analyzed images, are avoided. In this case, the UAV solution is a cheaper and more flexible one which can ensure superior image resolution even under adverse weather conditions. In this direction, the authors in [16] developed a solution for detection and evaluation of the dynamic evolution of the flood based on a collaborative team of UAVs. More recently a multicopter-based photogrammetry procedure was used to evaluate the effect of an earthquake on complex architectural landscapes [17]. Also, Feng et al. [18] used a UAV for urban flood monitoring and showed that such platforms can provide accurate flood maps. In their proposed method, the authors show how the acquired images are ortho-rectified and combined into a single image. Subsequently, the flood detection is realized through feature extraction from gray co-occurrence matrix and forest tree classifier methods. Boccardo et al. [19] compare the main advantages/disadvantages of fixed-wing UAV versus rotor platforms for area surveillance.
So, rotor platforms can be used only for very small areas or isolated buildings, while for small and medium areas fixed-wing UAVs are recommended. For large areas, UAV teams with pre-positioned stand-by can successfully perform the aerial surveillance of the disaster affected areas. Systems using UAVs are able to operate at lower altitude and to acquire multi-view, repetitive images with high resolutions [20]. These systems (fixed-wing type) are used to provide large image blocks to perform an image-based registration on multi-temporal datasets.
Control of a team of collaborative agents (UAVs in our case) is challenging, especially so under external perturbations, loss or delay of communication, etc. Therefore, the usual approach is to have a hierarchical control structure: the lower-level controller (the "autopilot" implemented on-board tracks a given reference) and the higher-level controller ("mission flight management") provides a reference trajectory [21].
Any mission has specific design requirements for the trajectory generation procedure [22,23]. Foremost in observation missions (surveillance, photogrammetry, target tracking, etc.) is to maintain a constant velocity or to allow variation within a narrow band (such that the photos taken cover uniformly the area under observation [24]). Whenever a team of UAVs is considered, additional issues appear (e.g., collision and avoidance constraints [25]). Not in the least, the trajectory has to be feasible (in the sense that it respects the UAV dynamics) Additional limitations on trajectory length, obstacle and collision avoidance are also encountered.
A promising framework is the flat representation of dynamical systems [26]. This approach expresses the states and inputs independently though a so-called flat output (which "hides" the underlying link between the states and inputs). Relatively recent work, has concentrated on providing flat characterizations which handle well numerical issues and have a manageable complexity [27,28]. In this sense, B-spline functions are an interesting choice: their geometrical properties lead to a good flat output parametrization and allow construct optimization problems which integrate costs and constraints in order to obtain the desired results [29]. Assuming that all low-level control loops are already designed such that a predefined trajectory is followed accurately and the payload is stabilized, we can reduce the path generation problem to an optimization problem where various constraints, parameters and costs are taken into account. To conclude, a flat representation which accounts for the low-level dynamics of the autopilot (approximated by first and second-order dynamics) and uses the properties of B-spline basis functions will provide a comprehensive and flexible framework [30]. In particular, it is possible to penalize trajectory length and energy costs, guarantee obstacle avoidance and pass through or within a pre-specified distance from a priori given way-points.
In order to detect and segment the flooded regions from the acquired images, texture and color analysis can be employed. Texture analysis techniques are a subject extensively studied in literature [31][32][33][34], but all suggested solutions are tailored to the specific context of the application considered. Moreover, there is a substantial interest in studying methods using the grey level co-occurrence matrix for texture characterization, but extremely little work is done when multi-spectral (color) co-occurrence is concerned [33,34]. All the above image acquisition strategies impose strict constraints on the photographs' capture during the UAV mission, i.e., photographs have to be captured: at a constant height (low/medium/high-the classification is relative, depending on context and application); such that there is a predefined overlap between neighboring photographs and there are no gaps in the area of interest (such that a photo-mosaiced image covering the entire area is computed). While there are many specialized software applications which can merge photographs with partial overlap to generate a continuous mapping and detect features of interest, there are still open issues in the generation of the flight path to be followed by a UAV [22]. This apparently simple problem has a number of intricacies: turn maneuvers of the UAV should not "cut" into the shape of the area under observation, maximal distance between consecutive path lines has to be respected and, not in the least, the UAV operational costs (energy, time of travel) should be minimized [23,24].
In this paper, we implemented new solutions concerning the system concept, path generation and image segmentation. As first novelty, we propose a multilevel system, with two components, terrestrial and aerial, as support for image acquisition and transmission from a specified region. The terrestrial component comprises a Ground Control Station (GCS, as coordinator or master node at distance), more Ground Data Terminals (GDTs, as a fixed nodes network for data acquisition and transmission), and a launcher. The aerial component contains mobile nodes (UAV-fixed wing type). The communication is established via internet (GDTs-GCS) or direct radio (in rest). This hybrid network has the advantage of extending operational area. The fixed-wing type UAV for image acquisition was developed by the authors in the Multisensory Robotic System for Aerial Monitoring of Critical Infrastructures (MUROS) project [35] funded by the Romanian National Research Program Space Technology and Advanced Research (STAR) from the Romanian Space Agency (ROSA) [36]. The proposed system is completely autonomous, apart from the take-off stage where a human operator is needed, and can be monitored and controlled at distance from the operational field. The area to be monitored is covered with the aid of a trajectory designed by a suitable optimization problem while the acquired images are analyzed in order to detect and assess the extent of floods. The second novelty refers to previous trajectory design implementations. So, the main contributions are: (a) the full dynamics (GCS + autopilot levels)-described in the flat representation and (b) the area under surveillance-partitioned between UAVs, such that the workload is balanced and collision with another UAV is impossible. The third novelty is a new solution for detection and quantitative evaluation of flooded small areas, based on the gliding box algorithm and local image processing. The advantages of this solution are the following: lower cost compared to a manned aircraft or a satellite solution, better resolution than a satellite solution, and the possibility of operating on cloudy conditions. The proposed method simultaneously uses pixel distribution and color information taking into account the chromatic co-occurrence matrix and mass fractal dimension on color components. The features used are not fixed as in [18] but rather they are being adapted to each application and environment condition. Results of the feature selection (especially associated with color channels) eliminate the temporal (colors of vegetation) and geographical influences (soil and vegetation colors, buildings and infrastructures).
The rest of the paper is organized as follows: in Section 2, first, the model of the UAV system based on hybrid wireless network is presented and second, the methodology and algorithms for image processing with the aim of flooded area detection, segmentation and estimation are described and implemented. The results obtained from images acquired with a fixed-wing UAV, designed by a team including the authors, are reported and analyzed in Section 3. For image acquisition, a path generated by the method introduced in Section 2 is used. Finally, the conclusions and discussions are reported in Section 4 and, respectively, Section 5.

Instruments and Methods
It is difficult and expensive to obtain precise data of the flood size within a certain small area from aerial photographs. As it was stated in Section 1, in this paper we propose a cheap and accurate solution to estimate the size of the dispersed small flooded areas. The solution is based on image segmentation obtained by a hybrid aerial-ground network integrated in internet. Three important sides are investigated: (a) the configuration of the network (which was partially implemented in the MUROS project [35] and will be finalized in SIMUL project [36]); (b) the trajectory control; and (c) the image processing for flooded area detection and estimation. The entire system is monitored and controlled remotely by GCS, via the Internet.

UAV System
The proposed system is configured as a hybrid network both with fixed nodes (terrestrial part) and mobile nodes (aerial part). The terrestrial part consists of the following components, which are considered at fixed locations during the mission ( Figure 1): Ground Data Terminals (GDTs), Launchers (Ls), Ground Control Station (GCS) and Image Processing Unit (IPU). The aerial part contains mobile nodes (UAVs, fixed wing type) which fly over a specified flooded zone. GCS is located at distance from the operational field and the communication is made via GSM + Internet. Four wireless communication channels were used: GCS-GDT (GSM + internet), UAV-GDT (radio) and L-GDT (radio), and UAV-UAV (radio). GDT-GCS connection uses a modem GPRS/LTE as router via Ethernet interface. It is a Virtual Private Network. The block diagram of the system consists of several modules, wired to a common control bus. Each module contains a Central Processor Unit (CPU), a Power Supply Unit (PSU), and a Controller Area Network (CAN) adaptor. The wireless module is characterized by the following: (a) radio modem; (b) frequency: telemetry-[3.3 GHz-3.5 GHz], video-2.4 GHz; (c) Data rate: telemetry-230 kbps, video-analog PAL; (d) range: telemetry-20 km, video-15 km. The significance of the module abbreviations in Figure 1 and their functions are presented in Table 1.   Figure 1 and their functions are presented in Table 1. Figures 2 and 3 present the principal components of UAV system, used for flood detection: UAV MUROS, GCS, GDT with ID box, Payload with camera, and Launcher.

FMCU Flight and Mission
Control Unit -Coordinates the flight mission; -Provides the platform's stability and quick response in case of disturbances that may deflect the drone from its pre-defined route or its removal from the flight envelope; -Allows for manual piloting by an operator on the ground; -Implements the automated low-level control loops which assure path tracking.

Trajectory Control
For image acquisition, the UAVs must follow specific trajectories such as simultaneously cover the monitored area ( Figure 4). As stated earlier, we propose to use flat output characterizations to describe the dynamics of the UAVs and further use B-spline parameterizations of the flat output in order to enforce various constraints and penalize some desired cost in the resulting constrained optimization problem.

Trajectory Control
For image acquisition, the UAVs must follow specific trajectories such as simultaneously cover the monitored area ( Figure 4). As stated earlier, we propose to use flat output characterizations to describe the dynamics of the UAVs and further use B-spline parameterizations of the flat output in order to enforce various constraints and penalize some desired cost in the resulting constrained optimization problem. Let us consider the nonlinear dynamics in standard notations [37]: such that the states and inputs can be expressed in terms of ) (t z and its higher-order derivatives: Further, let us consider the simplified UAV dynamics with north, east, down directions ( n p , e p and h) and yaw angle  as states: The autopilot is assumed to control directly the fligth-path angle γ , airspeed a V and roll angle accordingly chosen. Note that the closed-loop dynamics of the autopilot are simplified to first-order dynamics (a reasonable assumption in many circumstances).
Using the flat output we may express the dynamics in their flat representation as follows: Let us consider the nonlinear dynamics in standard notations [37]: where x(t) ∈ R n is the state vector and u(t) ∈ R m is the input vector. The system (1) is called differentially flat if there exists z(t) ∈ R m such that the states and inputs can be expressed in terms of z(t) and its higher-order derivatives: Further, let us consider the simplified UAV dynamics with north, east, down directions (p n , p e and h) and yaw angle ψ as states: The autopilot is assumed to control directly the fligth-path angle γ, airspeed V a and roll angle ϕ through input elements γ c , V a c and ϕ c , respectively: with parameters b γ , b V a and b ϕ accordingly chosen. Note that the closed-loop dynamics of the autopilot are simplified to first-order dynamics (a reasonable assumption in many circumstances).
Using the flat output z = [z 1 z 2 z 3 ] T = [p n p e h] T we may express the dynamics in their flat representation as follows: The major difficulty lies in the fact the constraints and costs are expressed as functions of state and input which do not necessarely translate well in formulations involving the flat output z. The usual solution is to parametrize the flat output after some basis functions (B d,i (t)): and to find the parameters P i which are, in some sense, feasible and optimal. Here, the parameter d denotes the degree of the B-spline functions. That is, each B-spline function can be seen as a piecewise combination of polynomial terms of degree d. Due to the particularities of the construction, a B-spline function of order d is continuous at least up to its (d-1) derivative. B-splines, due to their properties [30], permit to express the constrained optimization problem in terms of their control points P i (grouped here in column form in matrix P): , P) are short-hand notations which denote the cost, equality and inequality constraints, respectively. The cost can impose any penalization we deem necessary (length of the trajectory, input variation or magnitude, etc.) and constraints cover way-point validation, input magnitude constraints, etc. In general, a problem like Equation (8) is nonlinear and hence difficult to solve (particular solutions exploit the geometrical properties of the B-spline functions and/or heuristic methods). Considering multiple UAVs further increases the difficulty of the problem. In particular, we need to decide how the way-points are partitioned between the UAVs. One, rather cumbersome, solution is to attach to each way-point a binary variable and force that at least one of the UAVs passes through it. In practice, this can be relaxed, without any loss of generality to a condition which assumes that each UAV covers a contiguous part of the surveilled region. Moreover, it makes sense to partition the regions into areas of equal length parallel with the direction of travel. Then, each UAV has to cover its own independent region with additional collision avoidance constraints which may become active around the edges (since the UAV make turns which get out from under their surveillance area). To cover this possibility we may consider collision avoidance at the autopilot level (proximity sensors) or, more robustly, at the GCS level by either introducing additional constraints in the trajectory design procedure or, preferably, by changing the start and end points for each of the agent (such that neighboring points are reached at different moments in time).

Image-Based Flood Detection System
In order to fulfill the mission of detection, segmentation and estimation of the flooded areas, successive images are taken with constant rate on the pre-determined trajectory, like in the above section. For flooded area estimation, a patch-based segmentation was used. So, each image is partitioned in small boxes (e.g., in our application, patches of dimension 50 × 50 pixels), using a partitioning algorithm of images [38]. Note that the patch dimension is chosen depending on the image resolution and the texture of the segmented RoI (in our case, the flood). From a cluster of such patches (boxes), manually selected, a group of efficient features for flood detection is established based on a performance indicator. The features are used to create two classes: flood class (F) and non-flood class (NF). The propose method for image processing and interpretation has two phases: the learning phase and the mission phase. Both the learning images and test images were captured by the same camera device. Because the characteristics of the flood images can differ for each application, the learning phase is necessary to establish the class representatives and the signature patch structure. In the mission phase, a trajectory covering the investigated area is established. The acquired images are concatenated and processed to create an orthophotoplan without overlapping and without creating gaps. To this end, an overlap of 60% between two adjacent images is necessary to create an orthophotoplan. Then, they are indexed with an ID number in chronological order and are partitioned in the same way as in the learning phase. Based on the features selected in the learning phase, a similarity criterion is used to assign each patch to the class F or NF. Finally (estimation step), on one hand each patch of F is marked with white and is returned to the initial image, and, on the other hand a binary matrix of patches (BMP) is created with logical 1, if the correspondent patch belongs to F, and 0, in rest. By counting the "1"s from BMP, taking into account the total number of patches, the relative flooded area is evaluated.
The image characteristics may change as a function of distance from the ground and camera inclination with respect to the vertical axis. To avoid such issues the UAV has to respect a few additional constraints: (a) the altitude remains constant (even through some ground areas, may have different heights, we take as reference the water level, which remains constant). Floods are approximately at the same distance from the UAV, hence, if the flight plan is accurately followed, the resolution remains approximately the same for a given reference altitude; (b) the payload camera has to be oriented such that the lens are perpendicular to the surface of the Earth.
For each UAV there is a channel in GCS for image acquisition and, at the end of the mission, the images from all the UAVs are stored and processed in IPU. The methodology for flood evaluation based on patch analysis consists in the following steps:

1.
In IPU, ortho-rectified images are created and then they are combined into a single image without overlapping and without gaps (orthophotoplan).

2.
From the orthophotoplan, adjacent cropped images of dimension 6000 × 4000 pixels are investigated for flood evaluation. 3.
non-overlapping box decomposition of the tested image is made. So, a grid of boxes is created and its dimension will represent the resolution of flood segmentation. Thus, if the image dimension is R × M = 2 r × 2 m and the box dimension is 2 u × 2 v , then the resolution of segmentation The flood segmentation is made by patch classification in two regions of interest (flood-F and non flood-NF) taking into account the patch signatures and class representatives, which contain information about color and texture. As we mentioned earlier, the process has two phases: the learning phase (for feature selection and parameter adjustment) and the mission phase (for flood detection, segmentation and evaluation). Flood evaluation is made for each cropped image and, finally, the sum of partial results is calculated.

Learning Phase
Generally, the aerial images taken from UAVs have textural aspects. Moreover, the remote images for water (and particularly for flood) are characterized by high contrast in texture behavior between the flooded zones and the remaining soil. Therefore, the texture information and, in particular, texture features can be used for flood detection. The selection of effective features must group the patches with flood and differentiating them from the non-flood ones (it must increase the between-class separability and decrease the within-class variance). To this end, a set of significant texture features were analyzed in the learning phase, in order to select the most efficient ones for the classification process. The tested features are of different types: mean intensity (Im), contrast (Con), energy (En), entropy (Ent), homogeneity (Hom), correlation (Cor), variance (Var), mass fractal dimension (Dm), lacunarity (L) and histogram of Local Binary Pattern (LBP). They take into account the chromatic information as well (on R, G, B, H, S and V color components). The general formulas for the most used features in texture classification are given in Table 2, where: R is the number of rows of the image representation (matrix I), M is the columns of I and K represents the levels on color channels. C d is the normalized co-occurrence matrix [38] calculated as an average of the co-occurrence matrices  [39]. Dm (15) is calculated, based on Differential Box-Counting (DBC), for monochrome images, in [40]. A grid of boxes is created with the image divided in boxes with the factor r. For a box in position (u, v), the difference n r (u, v) between the maximum value p(u, v) and minimum value q(u, v) of the intensity are considered. Then, the sum of all the differences (17) is used to evaluate Dm. Similarly, the lacunarity L(r) is calculated as in [38]. Energy To evaluate the characteristics derived from co-occurrence matrix, besides the classical gray level co-occurrence matrix [33], applied on each color channel, we used the mean color co-occurrence matrix (CCM), between pairs of two spectral components of an input image [41]. So, in H, S, V decomposition, the image I is seen as a three-dimensional array with R rows, M columns and 3 layers (spectral bands). Each array element can take L positive integer (discrete values representing the color component's intensity of each pixel). The image I can be mathematically defined as: I ∈ N R×M×3 ).
Let H and S two components of a color space H, S, V. So, the mean CCM is considered as a square matrix, having L × L elements in N. It has two parameters: the distance d (the co-occurrence is the same as in GLCM case), and the component-pair (H, S) between which it is calculated. Each element of the mean color co-occurrence matrix CMM HS d (i, j) represents how many times a pixel of the H component, having an intensity level of i, is located near a pixel with intensity j in the spectral component S, at a d distance. Then, the elements of the mean CCM are [37]: Obviously, the next symmetry can be easily demonstrated: A simple example of calculating the mean CCM is given in Figure 5, where we consider two image components H, S ∈ N 3×4 , having 4 levels of pixel intensity, and the mean CCM computed between these two components, along a distance d = 1: A simple example of calculating the mean CCM is given in Figure 5, where we consider two image components 4 3 ,   N S H , having 4 levels of pixel intensity, and the mean CCM computed between these two components, along a distance d = 1: The algorithm for calculating CMM is presented [41] and the pseudocode in Appendix A. In order to establish the features to be selected, a cluster of 20 patches containing only flood (PF) are considered to form the representatives of the class "flood" (F) and 20 patches containing non flood elements (PNF), e.g., buildings and vegetation, are considered for the class "non flood" (NF). Each candidate feature Ti to flood signature is investigated according to the following algorithm: i.
Ti is calculated for all the learning patches (PF) and the confidence interval where, η(A) is the cardinal number of the set A. iv. The features Ti with greatest i CI are selected in decreasing order, until the fixed number of features imposed for flood signature is reached. For example, in Section 3, a signature T, with 6 elements is considered (12): v. As a consequence of the signature T, a set  of confidence intervals is created (13).  will be the representative of the class F: where: vi. For each selected feature Ti a weight wi is calculated as follows. Another set of 100 patches (50flood and 50-non flood) is considered and the confusion matrix for the feature Ti is calculated based on a preliminary classification criterion: The weight wi is established as in Equation (15): The algorithm for calculating CMM is presented [41] and the pseudocode in Appendix A. In order to establish the features to be selected, a cluster of 20 patches containing only flood (PF) are considered to form the representatives of the class "flood" (F) and 20 patches containing non flood elements (PNF), e.g., buildings and vegetation, are considered for the class "non flood" (NF). Each candidate feature T i to flood signature is investigated according to the following algorithm: i.
T i is calculated for all the learning patches (PF) and the confidence interval where m i and σ i represents, respectively, the mean and the standard deviation of T i . ii.
Similarly, T i is calculated for all the learning patches from PNF and the resulting set of values is noted as NF i . iii.
A confidence indicator for feature T i , CI i is created: where, η(A) is the cardinal number of the set A. iv.
The features T i with greatest CI i are selected in decreasing order, until the fixed number of features imposed for flood signature is reached. For example, in Section 3, a signature T, with 6 elements is considered (12): v.
As a consequence of the signature T, a set of confidence intervals is created (13). will be the representative of the class F: where: vi. For each selected feature T i a weight w i is calculated as follows. Another set of 100 patches (50-flood and 50-non flood) is considered and the confusion matrix for the feature T i is calculated based on a preliminary classification criterion: the patch B ∈ F if T i ∈ i .
The weight w i is established as in Equation (15): where F,F represents the number of patches manually selected as belonging to class F and classified to class F after feature T i . Similarly, F,NF represents the number of patches manually selected as belonging to class lF and classified to class NF after feature T i . Observations: a.
Obviously, CI i = 1 represents an ideal situation and are not encountered. b.
If λ i = λ j , then T i and T j are redundant and one can be eliminated.

Mission Phase
In the mission phase, the images from orthophotoplan are decomposed in patches with dimension of 50 × 50 pixels. Each patch (box) is indicated by a pair (row number, column number) in the squared grid of the image with an ID number. The mission phase has three steps: patch classification, image segmentation and flood estimation.
For classification of a box (B) of as flooded, a weighted vote D is considered (16), where D(B) is the sum of partial weighted vote for each selected feature (17): where: The patch B is considered as flood (18) if the weighted vote is greater than 0.8 from the sum of all weights (the maximum of D): where 0.8 is an experimentally chosen threshold. Inside of the analyzed image, a segmentation process is done with the aid of the detected flood patches. For visualization purposes, the flood boxes are marked with white. With the patches from an image, an associate matrix BMP is obtained. Each patch corresponds to an element in BMP; so, for an image dimension of 4000 × 6000 pixels and a patch of 50 × 50 pixels, then the BMP matrix dimension is dim BMP = 80 × 120. If the number of marked boxes is n, then the percentage of flood zone in the analyzed image is PF (22):

Algorithm for Flood Detection
The proposed algorithm has two phases: the Learning Phase-Algorithm 1 and the Mission Phase-Algorithm 2.

Inputs:
Learning patches (40 patches for feature selection-set 1 and 100 patches for weight establishing-set 2), set of feature to be investigated; Outputs: Selected features T i , the weights for selected features w i , and the intervals i , i = 1, . . . , 6.
For each patch of the first set: 1. Image decomposition on color channels (R, G, B, H, S, V) of patches; 2. Reject noise with median local filter; 3. Calculate the features: Im, Con, En, Hom, Ent, Var, Dm and L on color channels; 4. Until end of set 1; 5. Calculate the Confidence Indicator CI i for each feature based on Equation (11); 6. Feature selection: T i , i = 1, . . . , 6; 7. Determine the intervals for flood class representative i , Equations (13) and (14) For each Ti: 8. Calculate the confusion matrices CM i from the set 2; 9. Calculate the weights w i , i = 1, . . . , 6; Equations (15) and (17) 10. Return {T i , w i }.

Inputs:
Images to be analyzed, Selected features T i , the weights for selected features w i , and the intervals i , i = 1, . . . , 6;

Outputs: Segmented images and percent of flooded areas
For each image I: 1 Image decomposition in small non-overlapping patches (50 × 50 pixels); For each patch B 2. Calculate the selected features Im R , Con HH , En HS , Hom HH , Dm G and L R ; 3. Calculate Algorithm 1 is executed only once, at the beginning of the mission, while Algorithm 2 runs continuously throughout the mission. Both are implemented in deep neural networks (DNN). The DNN for Algorithm 2 is presented in Figure 6 and contains, besides the input and output layers, other three layers.
Layer 1 is dedicated to simultaneously calculate the features of patches and create the corresponding binary matrices of patches. Layer 2 is dedicated to local filtering of matrices from Layer 1, in order to eliminate the noise from BMP. Layer 3 creates the final BMP by voting scheme. Finally, the Output layer provides the segmented image and the relative flood size.

Experimental Results
For experimental results we used a UAV, designed, as coordinator, by University POLITEHNICA of Bucharest, MUROS project [35]. The main characteristics and technical specifications of UAV MUROS, as mobile node for image acquisition, are presented in Table 3. To evaluate the algorithms presented in Section 2, an image dataset of a flooded area was gathered with MUROS. The photographs have been captured along a path generated as in Section 2, with distances between lines d = 75 m and height of flight de = 100 m (wind strength was considered to be negligible). The portion from the orthophotoplan of an application near Bucharest, during a flood, is presented in Figure 7. The images analyzed with the algorithm described in Section 2 are marked with the specified IDs.
In the learning phase, for patch signature determination T, the first set of 40 patches of dimension 50 × 50 pixels (20 patches for flood and 20 for non-flood), manually selected, was used ( Figure 8). From this set, a cluster of 20 patches containing only flood (PF) are considered to form the prototypes

Experimental Results
For experimental results we used a UAV, designed, as coordinator, by University POLITEHNICA of Bucharest, MUROS project [35]. The main characteristics and technical specifications of UAV MUROS, as mobile node for image acquisition, are presented in Table 3. To evaluate the algorithms presented in Section 2, an image dataset of a flooded area was gathered with MUROS. The photographs have been captured along a path generated as in Section 2, with distances between lines d = 75 m and height of flight de = 100 m (wind strength was considered to be negligible). The portion from the orthophotoplan of an application near Bucharest, during a flood, is presented in Figure 7. The images analyzed with the algorithm described in Section 2 are marked with the specified IDs.
In the learning phase, for patch signature determination T, the first set of 40 patches of dimension 50 × 50 pixels (20 patches for flood and 20 for non-flood), manually selected, was used ( Figure 8). From this set, a cluster of 20 patches containing only flood (PF) are considered to form the prototypes for the class "flood" (F) and 20 patches, containing non flood elements (PNF), e.g., buildings and vegetation, are considered for the class "non-flood" (NF).
The results obtained in the learning phase (Table 4) show that the selected features (with CI criterion) are: ImR, ConHH, HomHH, EnHS, DmG and LR, where R, G, H, and S are the components of the color spaces. Thus, features on different types (first order statistics, second order statistics and fractal), on different channel color are selected. If CI falls below 0.80, then the accuracy can also decrease. It must be mentioned that the list of selected features can be changed in the learning phase, upon the requirements of the application. The fractal dimension was calculated by means of FracLac [42] plug-in of ImageJ and the features extracted from co-occurrence matrix were computed using MATLAB software. In Table 4, the values marked with * are those that are not within the corresponding confidence intervals.
Next step is the calculation of the confusion matrices for the selected features (Table 5). To this end, we used the second set (100 patches) for the learning phase, which contains 50 patches marked as flood (actually) and 50 patches marked as non-flood. From the confusion matrices we calculate the weights w i which will be used further for patch classification. The results obtained in the learning phase (Table 4) show that the selected features (with CI criterion) are: ImR, ConHH, HomHH, EnHS, DmG and LR, where R, G, H, and S are the components of the color spaces. Thus, features on different types (first order statistics, second order statistics and fractal), on different channel color are selected. If CI falls below 0.80, then the accuracy can also decrease. It must be mentioned that the list of selected features can be changed in the learning phase, upon the requirements of the application. The fractal dimension was calculated by means of FracLac [42] plug-in of ImageJ and the features extracted from co-occurrence matrix were computed using MATLAB software. In Table 4, the values marked with * are those that are not within the corresponding confidence intervals.
Next step is the calculation of the confusion matrices for the selected features (Table 5). To this end, we used the second set (100 patches) for the learning phase, which contains 50 patches marked as flood (actually) and 50 patches marked as non-flood. From the confusion matrices we calculate the weights wi which will be used further for patch classification.  In order to analyze the performances of the algorithm for flood detection, a set of 50 images with flood was investigated (see orthophotoplan from Figure 7). Random patches of flood and non flood types (Figure 9) are classified based on the voting scheme and the results are presented in Table 6.

T T T T T T R ConHH HomHH EnHS DmG LR
Here, D(B) is calculated as in (16)   So, the signature of the patch is: In order to analyze the performances of the algorithm for flood detection, a set of 50 images with flood was investigated (see orthophotoplan from Figure 7). Random patches of flood and non flood types (Figure 9) are classified based on the voting scheme and the results are presented in Table 6.
Here, D(B) is calculated as in (16) and compared with maximum value of 0.8 · ( s ∑ i=1 w i ) = 5.59 as in (17). For example, patches B6_F and B10_F with flood are wrongly classified as non flood. For the mission phase, an example of 6 images used for the algorithm test is presented in Figure 10 and the result of the segmentation, in Figure 11. Figure 12 overlaps the RGB images with masks generated by segmented images.
The random errors of the classification process are characterized by sensitivity, specificity, and accuracy [10,43] which are calculated in Table 7, where: TP is the number of true positive cases, TN is the number of true negative cases, FP is the number of false positive cases, and FN is the number of false negative cases. In [12] an accuracy of 87% is obtained using RGB information and six texture features (fixed) extracted from gray level co-occurrence matrix. Our method uses selected features (selected by a performance criterion at the beginning of the segmentation operation) on color channels (chromatic co-occurrence matrix and fractal type) and the accuracy was of 98.1%. the segmentation, in Figure 11. Figure 12 overlaps the RGB images with masks generated by segmented images. The random errors of the classification process are characterized by sensitivity, specificity, and accuracy [10,43] which are calculated in Table 7, where: TP is the number of true positive cases, TN is the number of true negative cases, FP is the number of false positive cases, and FN is the number of false negative cases. In [12] an accuracy of 87% is obtained using RGB information and six texture features (fixed) extracted from gray level co-occurrence matrix. Our method uses selected features (selected by a performance criterion at the beginning of the segmentation operation) on color channels (chromatic co-occurrence matrix and fractal type) and the accuracy was of 98.1%.    the segmentation, in Figure 11. Figure 12 overlaps the RGB images with masks generated by segmented images. The random errors of the classification process are characterized by sensitivity, specificity, and accuracy [10,43] which are calculated in Table 7, where: TP is the number of true positive cases, TN is the number of true negative cases, FP is the number of false positive cases, and FN is the number of false negative cases. In [12] an accuracy of 87% is obtained using RGB information and six texture features (fixed) extracted from gray level co-occurrence matrix. Our method uses selected features (selected by a performance criterion at the beginning of the segmentation operation) on color channels (chromatic co-occurrence matrix and fractal type) and the accuracy was of 98.1%. Figure 8. Patches for establish the flood signature (Pi_F as patch with flood and Pj_NF as non flood patch).                     Table 6. Some experimental results concerning the patch classification based on voting scheme. Gray rows mean wrong classification.

Discussion
Because we considered only complete flooded boxes, the approximation will be underestimated. Similarly, if mixed boxes are considered, an over approximation will be obtained. Table 8 presents the number of patches considered as F and the corresponding percent of flooded area for each images. The main cause was the patches from the contour of flooded which appear as mixed ones. Further studies will consider the decomposition of these patches in boxes increasingly small. On the other hand, by properly choosing textural features on color channels and patch dimension, the proposed algorithm can be extended to more classes like: road, vegetation, buildings, etc. Combining thermal camera with video, the system is able to detect possible persons in difficulty and to monitor the rescue operation. In this case, a flexible and dynamic strategy for trajectory design is necessary. Also the collaboration between the mobile nodes (UAVs) will improve the mission efficiency. The algorithms used for trajectory design minimize total path length while in the same time passing through (or within predefined distance) of a priori given way-points. In further work we plan to: (i) reconfigure trajectories on the fly such that the flooded areas are covered efficiently; and (ii) partition the workload of the UAVs such that total time/effort is minimized (for now we simply divide the area of observation into disjoint regions, one per each UAV).

Conclusions
The paper presented a comprehensive system and methodology for the detection and segmentation of flooded areas in a pre-determined zone. The contributions are focused on two important objectives: the planning of an optimal trajectory to cover the area under investigation and the image processing required to detect and to evaluate the flood spread. For the first the novelty lies in the analysis and computation of an optimal path covering the area of interest and for the second, the novelty lies in combining the information for different color channels with information about spatial pixel distribution obtained from chromatic co-occurrence matrix and mass fractal dimension. First, the paper studied a typical photogrammetry problem through the prism of control and optimization theory. That is, for a given polyhedral region which has to be covered by parallel lines (along which photographs are taken) we have given both an estimation of the required number of photographs and provided a minimum-length path covering the area. For the latter case we formulated a constrained optimization problem where various constraints and parameters were considered in order to obtain a minimum-length path. We took into account the maximum distance between consecutive lines and turn conditions (such that the UAV is guaranteed to follow the interior lines). We have also discussed the path generation problem in the presence of wind and for regions with non-convex shapes. Second, a methodology for the detection, segmentation and evaluation of flooded areas from the acquired images was presented. A color co-occurrence matrix was introduced and some efficient features. Furthermore, we illustrated that fractal type features on color component improve the local classification process on flooded and non-flooded boxes. The algorithm was tested on a large number of sub-images and the results showed good performances. We conclude that, by including the color information to texture analysis, by selection of feature based on maximum criterion and by using the fractal techniques, the accuracy of the detection of flooded boxes was increased up to 99.12%.

Conflicts of Interest:
The authors declare no conflict of interest. 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, and in the decision to publish the results. for i = 0+offset [1] to A.length+offset [1] for j = 0+offset [2]