Flood mapping through principal component analysis of multitemporal satellite imagery considering the alteration of water spectral properties due to turbidity conditions

ABSTRACT We propose a methodology for flood mapping by remote sensing considering the alteration of the spectral response of water due to turbidity conditions using principal component analysis (PCA) of multitemporal satellite imagery limited to the visible and the near IR region of the electromagnetic spectrum. Crosta technique is applied to the resulting load matrix to select the PCs of interest, and to interpret the nature of the changes. Finally, it is proposed the use of linear unmixing, using the selected PCs as input, to classify the categories required for flood mapping. The proposed methodology is applied to mapping the flood areas after the tropical storm Manuel during September 2013 in Acapulco, Mexico. The proposed procedure is a more robust alternative to quotients or index-orientated approaches based on fixed spectral response of the water, e.g. the Normalized Difference Water Index. This approach can be useful to authorities, civil protection and other organizations dedicated to risk management during natural contingences to assess quickly the dimension of affected areas, without the expensive and complicated mobilization of recourses to the site, and to give a more efficient response.


Introduction
Floods are the natural phenomenon that cause more losses around the world (Amarnath 2013). As a result of quickly, unequal and disorganized urban sprawl, there are highly populated areas exposed to this phenomenon, increasing the potential damage to assets and population. Flood monitoring should be a priority topic in national and international agendas, and in many places the actions of mitigation, prevention and response have been strengthened under the paradigm of climate change. Nevertheless, these goals are far from an integral and effective flood risk management. It is very common, in any part of the world that authorities cannot prevent and respond quickly enough to assist the population affected by floods or to manage the contingency in large areas. This situation is even more dramatic and has more relevance in development countries.
Earth observation satellites (EOS) and remote sensing (RS) are ideal instruments to follow natural events like floods in large areas, due to its high spatial and temporal resolution and its capacity for acquiring quickly multispectral or radar images of the affected areas, allowing an early diagnostic of critic situation without the mobilization of recourses to the site. The use of radar satellites in recent years has gained advantages for flood monitoring due to its active sensor capabilities CONTACT Marco A. Torres mtorresp@iingen.unam.mx; Eduardo Reinoso ereinosoa@iingen.unam.mx (Pulvirenti et al. 2011). However, a quickly and precise analysis is still the principal dare for change detection through RS (Adams et al. 1995;Womble et al. 2006). Images from EOS have great advantages in the study of natural hazards, such as continental or coastal floods. Even though multispectral satellite images are limited by their wavelength and adverse meteorological effects that affect the quality and continuity of the records. There are available references in the literature on the use of RS for flood mapping (Jain et al. 2005;Xiao et al. 2005;Sakamoto et al. 2007), flood monitoring (Jain et al. 2006), flood damage analysis (Dewan et al. 2007), flood forecasting (Van Der Knijff et al. 2008;Patro et al. 2009), validation of mathematic models and runoff analysis (Machado & Ahmad 2007).
In multitemporal studies, it is very important that satellites provide images with periodic coverture and similar conditions, such as altitude and angle of incidence. It is recommended that images of different dates (scenes) were obtained from the same satellite and sensor with the same spectral range (Chuvieco 2002). An increase in spatial resolution of the main satellite constellations has motivated the development, improvement and use of more accurate methods of change detection through RS.
Early methods were generally pixel based. As the resolution in the image increased, more sophisticated approaches emerged (such as feature-based methods) for change detection (Ilsever & € Unsalan 2012). Some authors have proposed quotients or index-orientated approaches to delimitate the temporal or permanent waterbodies using the ranges of the electromagnetic spectrum where water has its particular spectral characteristics, the best known and most cited is the NDWI (Normalized Difference Water Index; Gao 1996), which usually uses the range of blue, green and near IR. The first studies of flood areas with Landsat imagery showed that the near IR was the more suitable band for the identification of flood areas (Frazier & Page 2000). However, as it will be explained, the analysis of waterbodies throughout RS is very sensitive in their spectral characteristics owing to changes in the composition and suspended sediments. So, these variables must be considered in the multitemporal studies and change detection analyses through RS.
The multitemporal changes in satellite imagery are mainly related to transformations in their spectral or spatial features (Robin 1998). Changes in spectral features are related to the alterations in the reflexed or emitted signal of the element of interest, and are the most common in RS. Furthermore, spatial changes imply transformations of shape and size elements, even though the spectral features are the same (Chuvieco 2002).
In the dynamic of imagery due to natural events, it is necessary to distinguish two aspects related to the changes: the velocity in which the event produces the changes in the scene and the persistence or duration of those changes. Both aspects are not necessarily correlated: events with short duration can have long-term effects, and vice versa.
When visible light passes through a water column, its intensity decreases exponentially while the depth increases; this process is known as attenuation. The intensity of attenuation differs in accordance with the range of the electromagnetic spectrum. In the visible region, the red is attenuated faster, followed by the blue. When water has not suspended sediments present a peak in the green that decreases while the water depth increases. In contrast, near IR presents a second peak of reflectance; as depth increases, the spectral response declines. Although the turbidity of the water can change the spectral signature, the exponential decay of light intensity with increasing depth results from two processes, absorption and scattering. The absorption involves the conversion of electromagnetic energy into other forms such as heat or chemical energy (e.g. photosynthesis). The principal elements that involve absorption are: algae (phytoplankton), inorganic and organic particles in suspension (excluding algae), dissolved organic compounds (yellow substances), water itself, which strongly absorbs red light and has a smaller effect on shorter wavelength blue light (Mumby & Edwards 2000).
In multispectral satellite images, red band is located in the strong chlorophyll absorption region, while the near IR band is located in the high reflectance plateau of vegetation canopies. But, the water spectral signature also has a higher response in the near IR. In the visible region, the water may have different responses due to different conditions such as: the clarity of water column, depth, turbidity, dissolved organic compounds and eutrophication. Originally, the spectral response of water presents a peak in the green band, but this characteristic can change among different waterbodies. Wavelength of electromagnetic radiation may interact with suspended particles in the water column and change direction (see Figure 1). This process of scattering is caused by inorganic and organic particles and increases with the suspended sediment load (turbidity) of the water. These changes in the spectral response of water must be taken into account for a correct multitemporal analysis in flood mapping or the monitoring of biophysical conditions.
The main aim of this study is to propose a methodology for flood mapping applying principal component analysis (PCA) to a multitemporal dataset (the scenes before and after the flood event). PCA is useful for change interpretation when spectral signatures are affected severely by the own event, in this case by the water's turbidity in the second time. This particular form of PCA generates new bands uncorrelated among them with a new data arrangement. Then, Crosta technique is applied to the resulted load matrix and is used to select the PCs that contain the changes of interest, and to interpret the nature of those changes. Finally, it is proposed the use of linear unmixing, using the selected PCs as input, to classify the categories required for flood mapping. The proposed methodology is applied to map the flood areas in Acapulco's outskirts produced after the pass of tropical cyclone Manuel in 2013.

Methodology
Adjacent bands in multispectral imagery often contain similar information which produces redundancies in a multispectral dataset. The statistical technique known as PCA is used to reduce the data Source: ASTER spectral library (Baldridge et al. 2009). redundancy to perform more efficient data analysis and easier interpretation (Chuvieco 2002). This process typically results in fewer principal components uncorrelated among them, which are ordered in terms of their explanatory capacity (higher order components contain more information). The PCA has been used frequently as a spectral technique to enhance visual interpretation or as a feature extraction technique for subsequent use in digital classification procedures (Weng 2012).
Change detection analysis requires a multitemporal scene as input (i.e. a scene before the change and another scene after). The aim of PCA in change detection analysis is not to retain common information between two or more scenes, but rather to highlight the information that changed. In this sense, the resulting primary components are not very interesting because they contain the common information in both scenes, privileging the albedo response in the first PC. On the other hand, the secondary components offer the lesser common information, in this context this is precisely the object of interest, the information that changed. However, this transformation generally is complex and difficult to interpret directly. So, it is proposed to use the Crosta technique (Crosta et al. 1979) to select the PCs where the changes in the waterbodies are contained. Finally, a linear classification, as unmixing, is used to mapping the flood areas. The structure of the proposed methodology is shown in Figure 2, and it is discussed in what follows.

Image pre-processing
The image pre-processing may be different according to the quality level of each image. In this study the following processes were performed to both images: geometric adjustment, conversion of the digital numbers (DN) to reflectance values, and the atmospheric correction.
Superposition errors in satellite images are more representative during multitemporal interpretation; therefore, different scenes should adjust with higher level of precision. Authors suggest a precision higher than 0.2 pixels to guarantee errors below 10% in multitemporal comparison  (Townshend et al. 1992;Dai & Khorram 1998). It is necessary that images present a good correlation among the analysed scenes. If the geometrical adjustment of the images were not corrected or improved, there can be omission or commission errors; this is, changes should be assumed where do not exist, and vice versa.
When waterbodies are studied throughout RS techniques, it is necessary to convert the DN into a geophysical variable. When reflectance values are utilized, a homogeneous parameters are adjusted between the images, and the multitemporal analysis will guarantee that the elements of interest be in a natural and standardized scale. The absolute method DOS (Dark Object Subtraction), also known as Chavez method, was applied, obtaining a new histogram to a pristine form, diminishing the minimum value to each wavelength to obtain images independent of any kind of atmospheric effect (Chavez 1998).

Principal component analysis (PCA)
To perform the PCA is necessary to count with multitemporal and multiband images. The input images must share the same dimensions (rows and columns), pre-processing level, number of bands for each scene, geographical extension, format and preferably the same angle of incidence.
The purpose of using PCA is to reduce the dimensionality of the data to maximize the amount of information from the original bands into the least number of PC, in this case, the number of original bands. A set of correlated variables (original bands) is transformed in other uncorrelated variables (PC) which contain the maximum original information with a physical meaning that needs to be explored.
Assuming a multitemporal image with eight input bands (four for each time scene), it can be expressed in matrix format in the following way: where n represents the number of pixels and b the number of bands (eight in a SPOT multitemporal composition). Considering each band as a vector, the above matrix can be simplified as follows: . .
where k is the number of bands. To reduce the dimensionality of the original bands the eigenvalues of the covariance matrix must be calculated. This matrix can be calculated as follows: where s i;j is the covariance of each pair of different bands.
where DN p;i is the digital number of the pixel p in the band i; DN p;i is the digital number of the pixel p in the band j; m i and m j are the averages of the DN for the bands i and j, respectively. From the variance-covariance matrix, the eigenvalues (λ) are calculated as the roots of the characteristic equation: where C is the covariance matrix of the bands and I is the diagonal identity matrix. The eigenvalues indicate the original information that is retained. From these values the percentage of original variance explained by each PC can be obtained calculating the ratio of each eigenvalue in relation to the sum of all them. The PC can be expressed in matrix from: . .
where Y is the vector of the principal components, w the transformation matrix, and x the vector of the original data. The coefficients of the transformation matrix w are the eigenvectors that diagonalized the covariance matrix of the original bands. These values provide information on the relationship of these bands with each PC (Figure 3). From these values it is possible to link a main PC with a real variable. The eigenvectors can be calculated from the vector-matrix equation for each eigenvalue l k .
where C is the covariance matrix, l k is the k eigenvalues (eight in a SPOT multitemporal composition), I is the diagonal identity matrix, and W k is the k eigenvectors (Estornell et al. 2013).

Crosta technique
An inherent problem of PCA is the difficulty of interpreting a priori the content of the PCs, because the transformation is purely statistical and hence is high-dependent on the numerical characteristics of the image. Some authors have proposed that PC1 indicates the general bright and the PC2 the greenery (Ingebritsen & Lyon 1985). However, these assumptions imply that the image contains a sufficient quantity of vegetal cover (Vald es 2015). Therefore, general rules cannot be used for PC interpretation. The Crosta technique (Crosta et al. 1979) is also known as a feature oriented to principal component selection. The analysis of the eigenvector values allows the identification of the PCs that contain spectral information of specific targets, as well as the contribution of the original bands to the components of the spectral response of interest. This technique indicates the materials that are represented as bright or dark pixels in the PCs according to the magnitude and sign of the eigenvector loadings (Ranjbar et al. 2004).
The load matrix analysis, as a result of the application of multispectral transformation technique, is very useful to spatially locate and to understand the spectral nature of the most relevant changes of the input. When PCA is used in studies of change detection, the input raster is a multiband and multitemporal file. In this way, depending on the weight and direction of the eigenvectors, the principal components will be generated, concentrating most of the variance in the PC 1 and decreasing the variance in the following principal components. Each PC has a preponderant weight in some of the bands which are used like input in the process (Fraser 1991;Loughlin 1991).

Linear unmixing classification
Linear unmixing determines the relative abundance of materials that are depicted in multispectral imagery based on their spectral characteristics. The reflectance at each pixel is assumed to be a linear combination of each material. Linear Unmixing classification assumes that each surface component within a pixel is sufficiently large so that no multiple scattering exists between the components (Singer & McCord 1979). The use of linear mixture modelling relies on four assumptions (Settle & Drake 1993): (1) There is no significant occurrence of multiple scattering between the different surface components.
(2) Each surface component within the image has sufficient spectral contrast to allow their separation.
(3) For each pixel the total land cover is the unity.
(4) Each surface component is known.
The linear mixture model is expressed by the next equation: where R i is the reflectance for the ith pixel, r j is the spectral reflectance of the jth surface component, f ij is the fraction of the jth surface component in the ith pixel.

Application
The aim of this study is mapping flood areas from satellite multispectral images with high spatial resolution and low spectral resolution (only four bands of the SPOT-6 images). The proposed methodology is applied to mapping flood after the pass of the tropical storm Manuel during September 2013 in Acapulco, Mexico. During this month, the Mexican territory was affected by the simultaneous pass of two tropical cyclones, one in the Pacific Ocean (Manuel) and another one in the Gulf of Mexico (Ingrid). The damages produced by both hurricanes were the largest since the earthquake of 1985 in Mexico City (CENAPRED 2014). Due to the dimensions of the contingency, the capacity of the authorities was overwhelmed. Although Manuel does not reach the coast of Guerrero, the rainfall over the mountains triggered thousands of landslides (Rodr ıguez 2014) and caused considerable floodplains in the principal basins and coastal areas (Figure 4).

Area of study
Acapulco is a coastal city located %300 km to the south of Mexico City. Its urban sprawl during the last decades has overstepped the bay where it was originally founded, up to the lagoons (Coyuca de Benitez to the east and Tres Palos to the west). The lack of planning in the outskirts of the city transformed the river and coastal ecosystems of the region, increasing the vulnerability of the population to hydrometeorological and geological events which are very common in the Mexican Pacific coast. During the last 25 years, Acapulco has been largely affected by hurricane Pauline in 1997 (Garc ıa 2006) and tropical storm Manuel in 2013 ( Figure 5).

Data set
Two time scenes of SPOT-6 data were acquired, before and after the occurrence of tropical storm Manuel (February and October of 2013) throughout ERMEX NG (Spanish acronym for Reception Station Mexico New Generation). The images corresponding to different seasons of the year (spring and summer). The obtained SPOT-6 data were ortho processing level, with four bands: blue (0.450-0.520 mm), green (0.530-0.590 mm), red (0.625-0.695 mm) and near IR (0.760-0.890 mm). The characteristics from the available satellite images are shown in Table 1. The correlation of bands can be appreciated in Table 2, where it is shown that the ranges corresponding to the visible spectrum are highly correlated. On the other hand, the near IR ranges are less correlated to others ranges of the electromagnetic spectrum.
Some authors have treated the temporality in change detection studies using RS (e.g. Robin 1998;Chuvieco & Huete 2002). Most of them agree that the time scenes employed must correspond to the same annual season due to the homogeneity of vegetation conditions among time scenes. However, this is very difficult to satisfy in change detection analysis due to adverse weather during hydrometeorological events.
In multispectral imagery, fluvial floods are often accompanied by sediments in the second time scene, complicating the characterization of waterbodies in images with only three channels in the visible range and one in the near IR. The spectral response of the water in second time scene changes significatively respect the first time scene, increasing the reflectance in the red band of the second  time scene due to the high concentration of clays in the water column. The characteristics of reflectance and absorbance of the elements in the image are strongly linked among them (Table 3). The identification of water pixels at the water/soil interface may be complicated because it depends on several factors. Deep waterbodies have quite distinct representation in the imagery than shallow waterbodies. However, very shallow water/turbid water can be mistaken with soils while saturated soils can be mistaken with water pixels. It is also possible that a pixel presents mixed conditions (some part of water and other part of soils). Therefore, very shallow water will have the same colour and brightness as very turbid water due to the high concentration of suspended sediments ( Figure 6). Besides, there is an important change in the greenery of vegetal cover among the two scenes.

Results
The interpretation procedure for PC needs more elaboration and classification, especially in terms of effects of different combination of negative/positive eigenvectors with high/low input image DN values (Gupta et al. 2013).
In the load matrix the magnitude of the new components is measured by the variance or the percentage of eigenvalue stored in each PC. As it can be seen in Figure 7 and in Table 4, the first five PCs collect 98.21% of total variance, as follows manner: PC 1 has 63.84%, PC 2 has 16.73%, PC 3 has 8.81%, PC 4 has 7.37% and PC 5 has 1.44%. Finally, PC 6, PC 7 and PC 8 have percentages under 1%. So, PC 1 collects common information of both scenes, whereas PC 2-PC 4 have the less common information (denoting principal changes) and PC 6-PC 8 have noise and residual information.
In this study, the changes of interest are related to the increment of the area of water at the second time scene respect to the first time scene. In addition, the sediments in the second time scene have a great repercussion in the response of water from a peak shift from the green to red and blue bands. As the bands are in the multitemporal composition input, the results obtained through PCA allow to integrate all possible changes in the scene: an increase or decrease of waterbodies, the sediments accumulated or missed and the geomorphological changes in the flood plain. The ability to detect all changes and to distribute them to all different PCs is the principal advantages respect to other methods or techniques of change detection. The approaches based on the comparison of thematic classifications only focus on the spatial detection of the changes but do not analyse the spectral changes among time scenes. The techniques based on indexes or multitemporal quotients use specific bands (e.g. green or near IR in the case of NDWI) which implies a higher probability of obtaining false positives. The PCA integrates a major quantity of data, using a higher number of bands, which increase the possibility to detect more changes and to isolate them into one or more of the generated components. This technique requires only a basic understanding of the spectral properties of the presented elements (Loughlin 1991) as the key for the correlation matrix. If the second time scene corresponds to the maximum flood and it does not have an obvious change in the water colour, the increment of the depth of water would be only detected by the use of green and near IR bands. However, the differences between waterbody in the first time scene and the second time scene imply a complex behaviour between the spectral responses of the same target. Water concentrates its spectral response in the near IR region due to its reflectance characteristics which is independent of external factors, instead of the visible region where the response depends highly on other variables from water composition. On the other hand, with the presence of sediments and organic material, it is valid to assume that the response can be identified in some of the PCs, and this can be done through the load matrix. For this reason, it is considered for flood Figure 6. Change in the spectral response of water bodies between the two times of the scene (before and after the tropical cyclone Manuel) in five different locations compared with the spectral response of distilled water. Note the increase in the reflectance value in the red region due to the increment of the sediments after the tropical cyclone Manuel. mapping those PCs with near IR, red and blue inputs, which correspond to PC 2, PC 3 and PC 5, respectively (Table 4).
PC 1 contains the effect of topography and albedo; as it can be expected that it is dominated by a positive contribution of all the input bands with a higher response of the green of the first time scene (green 1st). PC 2 presents a positive contribution of the near IR of the first and second time scene (near IR 1st and near IR 2nd), which reflect the effect of vegetation and water bodies. The negative values in the load matrix for the second time scene made possible a more precise delimitation of waterbodies in this time, in other words the flood areas. PC 3 is dominated by a positive contribution of near IR in the first time scene (Near IR 1st) and contrasted by the negative loads of the visible region of the first time; this allows the delimitation of the waterbodies in the first time scene, i.e. perennial waterbodies. PC 4 has a major response of the near IR 2nd and of the visible region in the same time scene which is related to the structure of the vegetation; the negative values of the visible region of the first time scene allow the identification of banks of sediments along the principal rivers. However, the utility of PC 4 was rejected because the characteristics of waterbodies in the second  time scene are best reflected in PC 2. The response of PC 5 is higher in the red 1st and red 2nd, and the lower response corresponds to the blue 1st which can be interpreted as the result of the change in the response of the water to the red band due to the increase of sediments in the second time scene (Figure 8). The rest of PCs (PC 6-8) contain information that is not useful for flood mapping.
New datasets were generated from these PCs, corresponding to waterbodies in the first time scene and waterbodies in the second time scene. Using these new data sets as input, four categories were classified by means of unmixing: (1) waterbodies of first time scene, (2) waterbodies of second time scene, (3) vegetation and (4) urban areas. Each category was represented by 50 samples. The results are shown in Figures 9 and 10. From these results, it is concluded that analysis of change detection through PCA is useful for the identification of changes due to floods when the correct variables are integrated and is performed a correct interpretation of the load matrix knowing the response of the target in the spectral channels used as input.

Discussion
Change detection using RS often implicates a dilemma between errors by omission or by commission; in other words, it implies to accept or reject changes that exist or detect changes which have not occurred. Commission and omission errors in change detection through PR are caused by temporality and spatial resolution of the images. It is very difficult to count on of useful images of the exact moment or just after the flood. This problem compromises the identification of areas that could be or could not be affected by the flood. Nowadays, there is a major number of multispectral satellites in orbit compared to the number of radars, which allows to count on of images from different constellations for the same event. Therefore, it is necessary to continue developing methodologies which use multispectral images from satellites with remote passive sensors. These images can be easily acquired by any person or institution. Multispectral images, in accordance with its spatial    resolution level and its origin, can be free and have historical data which can be useful for many types of studies. The analysis of the spectral response of waterbodies reveals a greater response in the near IR band rather than the blue or the green bands due to the increase of turbidity of the water. The use of multispectral sensors with more bands in the IR or in the visible range could improve the detail of the classification of waterbodies. Besides, the meteorological conditions can reduce the capacity of the sensors in the satellites (especially passive sensors) to detect the changes among the scenes.
At present, there are various multispectral sensors with high spectral resolutions but limited by their spatial or temporal resolutions (Landsat-8, Sentinel-2A, CBERS-4 and MODIS) which can be used with relative success depending on the characteristics of the event. More bands/channels could be useful to a more detailed delimitation of waterbodies.
As it was explained, the increase of sediments or other suspended agents after the detonating event represents an important factor in flood mapping due to the change of the typical spectral response of water. However, most of the existent approaches of flood mapping do not include this variation. The proposed approach based on PCA analysis takes into account this effect for a more complete and robust analysis.
The utility of the PCA resides in the mathematic transformation, which allows to obtain noncorrelated bands from spectral information highly correlated. However, the interpretation of the load matrix and the eigenvectors generally is confuse and produce wrong conclusions; additionally, the PCA could present dependency among scenes, and even depend on the software used. The changes in the scene are interpreted through the contribution of the original bands in the load matrix by means of Crosta technique. The PCs taken as indicators of flood areas are those with dominant contribution in the near IR and red bands, which are PC 2, PC 3 and PC 5 ( Table 4).
The results of change detection analysis are conditioned by the temporality of the images and by the persistence or duration of the effects of the event; it is much better that the images were taken in times closer between them, and closer to the moment when the changes are produced. For example, it is very complicated to obtain the floodplain of flash floods, so typical in arid zones, which are characterized by its short duration.
The PCA is focused on the identification of specific elements or characteristics which are difficult to recognize or to delimit by its spectral attributes. This technique was proposed by different authors for mineral prospection. Although these studies used only one image corresponding to a single time scene, the analysis of load matrix is useful in studies of change detection. Besides, these studies limited the application of Crosta technique to specific geographic conditions, e.g. desert conditions where the response of vegetation is very rare or null. The proposed approach can be used in different geographic conditions, even if the spectral response of first time scene differs considerably of the second time scene since the changes in waterbodies will be reflected in any of the PCs of the load matrix.

Conclusions
Conventionally, the application of PCA to change detection analysis has been focused mainly in vegetation cover and land use (Maldonado & dos Santos 2005;Ponce 2008;Rojas 2011) because of the difficulty to interpreting the load matrix and the contribution of principal components of the original bands (Lira 2010). The proposed methodology, based on PCA and unmixing classification, uses a multitemporal and multiband combination as the input that highlight the change of the spectral response of water, as a result of the presence of suspended sediments. This approach represents a more robust alternative to quotients or index-orientated approaches based on fixed spectral response of the water, such as the NDWI.
This kind of methodologies can be applied by authorities, civil protection and other organisms dedicated to risk management during natural contingences to assess quickly the dimension of affected areas and to carry out a more efficient response. RS allows to get in almost real-time the effects of a natural hazard, and to help in a rapid decision-making.