Tsunami-induced Change Detection Using SAR Intensity and Texture Information Based on the Generalized Gaussian Mixture Model

Abstract

The remote sensing technique using SAR data have many advantages when applied to the disaster site due to its wide coverage and all-weather acquisition availability. Although a single-pol (polarimetric) SAR image cannot represent the land surface better than a quad-pol SAR image can, single-pol SAR data are worth using for disaster-induced change detection. In this paper, an automatic change detection method based on a mixture of GGDs (generalized Gaussian distribution) is proposed, and usability of the textural features and intensity is evaluated by using the proposed method. Three ALOS/PALSAR images were used in the experiments, and the study site was Norita City, which was affected by the 2011 Tohoku earthquake. The experiment results showed that the proposed automatic change detection method is practical for disaster sites where the large areas change. The intensity information is useful for detecting disaster-induced changes with a 68.3% g-mean, but the texture information is not. The autocorrelation and correlation show the interesting implication that they tend not to extract agricultural areas in the change detection map. Therefore, the final tsunami-induced change map is produced by the combination of three maps: one is derived from the intensity information and used as an initial map, and the others are derived from the textural information and used as auxiliary data.

1. Introduction

After the occurrence of a severe disaster, rapid relief activities are extremely important for reducing the casualties. For the urgent relief activities, analysis of the disaster site is needed, but it is usually very difficult to immediately obtain information regarding the site through a field survey. Therefore, remote sensing data play an important role in the event of a disaster because remote sensors can observe the disaster site widely and rapidly. Among the various types of remote sensing data, SAR (Synthetic Aperture Radar) data are very efficient because of their all-time and all-weather acquisition ability. Because the SAR sensor is an active sensor, which means that it generate its own signal unlike optic sensor, it can observe the surface during day or night. In addition, as the signal of SAR sensor can penetrate clouds, SAR images can be acquired under all weather conditions (Zyl and Kim, 2011).

The previous research on detecting disaster-induced damage used the change detection method (Matsuoka and Yamazaki, 2004; Bovolo and Bruzzone, 2007; Gamba et al., 2007; Park et al., 2013). It compared the pre- and post-event images and extracted the changed areas as the damage areas. Matsuoka and Yamazaki (2004) calculated the backscattering coefficient and the intensity correlation between two ERS (European Remote Sensing) satellite images and extracted areas with lower values as the damage areas. Gamba et al. (2007) detected them using the intensity and phase information of multitemporal ASAR (Advanced SAR) images with GIS data. Park et al. (2013) extracted tsunami damage areas with the various polarimetric parameters of ALOS (Advanced Land Observation Satellite)/PALSAR (Phased Array type L-band SAR) images, and revealed that quad-pol (polarimetric) SAR data can detect damage areas more accurately than single-pol SAR data can. The four types of scattering information (HH, HV, VH, and VV) in quad-pol SAR data are decomposed into various polarimetric parameters, and such parameters reflect the various characteristics of the land surface. Therefore, quadpol data can well describe the features of the land surface.

The single-pol imaging mode, however, covers a greater area because of its wider swath width compared to the quadpol mode (Ainsworth et al., 2009). Furthermore, much more pre-event images are available for comparison with the postevent images than quad-pol SAR data when a disaster occurs because a good number of single-pol SAR images have been obtained over the past decades. From this point of view, the change detection method using single-pol data is appropriate for disaster-induced change detection. Many disaster-induced damage detection cases using single-pol SAR data compared simple factors such as intensity or correlation between pre- and post-event images. To compensate for the insufficient information on such factors, texture information can be applied. Texture information has been used in many previous studies for classification, feature extraction, and change detection (Soh and Tsatsoulis, 1999; Zhu et al., 2012; Kang et al., 2015). In particular, Kang et al. (2015) showed the potential use of texture information for urban change detection from VHR (very high resolution) SAR imagery.

Change detection for disaster-caused damage detection from SAR data should be automatic and unsupervised because prior knowledge such as ground truth information hardly exists in the disaster site. Unlike the use of optical data, automatic change detection using SAR data has been less exploited due to the SAR data’s speckle noise and geometry distortions (Bolvoro et al., 2013). Furthermore, because the statistical distribution is affected by the sensor type, the land cover type, etc., SAR data do not always follow the Gaussian distribution, which is traditionally envisaged in automatic thresholding algorithms like the K&I (Kittler and Illingworth, 1986) and Otsu (Otsu, 1975) methods. SAR intensity images generally follow various gamma distributions such as Rayleigh, Nakagami, Weibull, etc. (Li et al., 2007). The major impediment, however, to using a generalized gamma distribution is its complexity and difficulty in parameter estimation (Gomes et al., 2008). Therefore, automatic change detection based on a GGD (Generalized Gaussian distribution) can be an alternative due to its high flexibility. The change detection method based on the GGD showed the proper results when it applied to SAR images (Bazi et al., 2009).

In this paper, an automatic change detection method based on the GGD is suggested for detecting disaster-induced damage areas. By applying the proposed change detection to intensity and texture images, a change map is produced. The only two single-pol SAR images are needed in the proposed method as it is considered that various data such as DEM (Digital Elevation Model) and optic images are not always available after the disaster. Texture images are first generated by GLCM (Gray Level Co-occurrence Matrix). Statistical analysis based on the GGD is performed to DI (Difference Image) of each factor (intensity and textural features) through the EM (Expectation Maximization) algorithm. The final change detection maps are derived from the statistical analysis results. In the following section, the methods of generating texture information and of automatically detecting change based on the GGD through the EM algorithm are explained in detail. The experiment results are reported and compared with one another in section 3, using the ALOS/PALSAR data of the 2011 Tohoku earthquake. Finally, conclusions are drawn in section 4.

 

2. Methodology

The proposed method consists of three steps: preprocessing, extracting the textural features, and automatic change detection. The preprocessing step includes geocoding, coregistration, and reducing speckle noise, in that sequence. A transformed slant-range image, which is the natural result of radar-range measurement systems like the SAR system, is geocoded to a ground-range image, which is projected to a specific coordinate system. Co-registration and speckle filtering are also performed for better change detection performance. In this paper, geo-coding and co-registration were performed by the SNAP (Sentinel Application Platform) which was downloaded through STEP (Science Toolbox Exploitation Platform) of ESA (European Space Agency) (http://step.esa.int/main/download/). The SNAP is the common architecture for all Sentinel Toolboxes which are the processing tools to support the diverse data including Sentinel-1, ERS-1 &2, ALOS/PALSAR, etc. The enhanced Lee filter with a 5×5-sized window was used for reducing the speckle noise. The other steps are explained below.

2.1 Extracting the textural features

Texture information is about the spatial distribution of the pixel values in an image. GLCM is the most efficient method of representing the texture of an image (Lee et al., 2005). Three fundamental parameters must be defined when calculating the GLCM: the quantization levels, displacement (d), and orientation (α). GLCM is created by calculating the number of pixels with the gray value i, which are in the specific relationship to the pixels with the gray value j. The specific relationship is defined by distance (d) and orientation (α). Fig. 1 shows the creation of GLCM when the distance is set to 1 pixel and the orientation is set to 0°(vertical). The size of the resultant GLCM is affected by the size of the gray level of an image. In Fig. 1, the gray value varies from 0 to 3, and the size of the resultant GLCM is 4×4. If an image has n gray values from 0 to n-1, the size of the resultant GLCM is n×n. Therefore, the size of the gray level must be quantized for efficient calculation and for obtaining a desirable resultant texture information.

Fig. 1.Calculation of GLCM (distance=1; orientation=0°): (a) objective image; (b) GLCM

Haralick et al. (1973) suggested a set of 28 textural features using GLCM. In this paper, seven popular textural features were tested: autocorrelation, contrast, correlation, dissimilarity, energy (angular second moment), entropy, and homogeneity. The equations below define these features.

where p(i, j) is the value of the (i, j)th cell in the normalized GLCM.

where μi, μj, σi and σj are the means and standard deviations of row (i) and column (j).

If some relative frequencies p(i, j) are zero, log (0) is not defined. Therefore, the set of (i, j) where the value of p(i, j) is 0 is excluded during the summation.

2.2 Automatic change detection

A change detection method traditionally consists of two parts: DI generation and thresholding. As it was assumed, however, that the DI of SAR data follows the GGD, statistical modeling of the DI based on a mixture of GGDs using the EM algorithm was additionally performed before thresholding.

2.2.1 DI generation

The DI of intensity was directly calculated using the logratio method (Eq. (8)).

where Ipre, Ipost are the intensity values of the pre- and post-event images, respectively. The ratio method (r=I1/ I1) is desirable for SAR intensity images because unlike the difference method (d=I1−I0), it does not depend on the intensity level of the pixels and is very robust against calibration errors (Rignot et al., 1993). For the DI of textural images, log calculation was first applied to the pre-processed images.

The textural features that were used in this paper were calculated by the GLCMs. The difference method is desirable for texture feature images.

2.2.2 Statistical modeling of DI

(1) A mixture of GGDs

The Gaussian distribution, which is known as a normal distribution, is the most popular distribution in statistics (Nadarajah, 2005). When the Gaussian distribution is naturally generalized, it becomes a GGD. The pdf (probability density function) of GGD is given by (Elguebaly and Bouguila, 2011)

where x (x ∈ ℝ) is the random variable which is a gray value in this paper, and θ is the set of statistical parameters, including μ, α and β which denote the mean, scale, and shape parameters, respectively. The Γ(∙) is the Gamma function. The shape parameter, β, determines the flatness of the pdf (Fig.2). By changing the value of β, the GGD is made to include a variety of statistical distributions, such as the uniform, Gaussian, Laplacian, and impulsive distributions. For example, when β = 2, the pdf is the Gaussian distribution, and when β = 1, the pdf is the Laplacian distribution. The flexibility of the GGD is one of the reasons that many previous works adopted it.

Fig. 2.GGDs with different shape parameter β (Elguebaly and Bouguila, 2011)

It was assumed that the DI follows a mixture of two GGDs: one for a class of changed pixels and the other for a class of unchanged pixels. If the random variable, x, follows a mixture of k GGDs, then pdf is represented as

where Pi is the prior probabilities and θi is the set of statistical parameters associated with the ith distribution. It was noticed that an ith GGD has four unknown parameters: Pi, μi, αi and βi. Therefore, modeling an image with a mixture of k GGDs means that 4k parameters should be estimated. This problem can be solved by the maximizing the log-likelihood using the EM algorithm.

(2) EM algorithm

The EM algorithm (Dempster et al., 1977) is well known for the estimation of unknown parameters for a mixture model and used in several studies (Bazi et al., 2007; Yang et al., 2012; Nguyen et al., 2014). The EM algorithm for the parameters of the GGD was well explained in Bazi et al. (2007). If an image, X, is a mixture of k GGDs and is composed of N pixels, the likelihood of X is

where Θ refers to the entire set of parameters to be estimated. The EM algorithm considers the data incomplete and having a missing part. Therefore, each pixel, xj, is associated with zj for applying the EM algorithm to the image, X, which means a hidden variable that indicates which distribution includes xj. For example, if xj belongs to the ith distribution, zji is equal to 1. Otherwise, zji is equal to 0. The likelihood of a complete X is given by

The log-likelihood of Eq. (12) is expressed as

where Z is the set of zj. Using Eq. (9), the log-likelihood of X is given by

Each zij can be replaced by its conditional expectation associated with xj and the parameters of the ith distribution (Eq(15)).

The EM algorithm has an iterative construction consisting of two steps: the expectation step and the maximization step. In the expectation step, zij is computed using the current statistical parameter Θ In the maximization step, Θ is updated using Eqs. (16) - (19). These formulations are computed based on the partial derivatives of the log-likelihood of X (Eq. (14)) with respect to each unknown parameter.

where t denotes the number of iteration steps. is obtained using Eq. (15) using and the parameter at the tth step. Ψ(∙) is the digamma function and ξ(∙) is expressed as Eq. (20).

As Eq. (18) and (20) are nonlinear equations, the Newton-Raphson method was used to solve them. The EM algorithm stops when the number of iterations reaches the maximum or when the difference between Θt and Θt+1 is less than the threshold. The maximum iteration number was set to 100, and the threshold was set to 1.0×10−5 for the EM algorithm.

The limitation of the EM algorithm is that its result is not always the global maximum and is sometimes a local maximum. In this paper, the EM algorithm was initialized with automatic thresholding algorithms such as the K&I and Otsu methods to avoid the local maximum problem.

(3) Thresholding

It was assumed that the distribution of DI was a mixture of two GGDs. The two GGDs indicated the changed and unchanged classes. As a result of the EM algorithm, the unknown parameters (P, β, α and μ) were determined, and the probability of x associated with each class could be calculated. T in Eq. (21) is the Bayes minimum error threshold, where the image should be divided (Kittler and Illingworth, 1986).

 

3. Results

3.1 Study site and data

The March 11, 2011 Tohoku earthquake, also referred to as the Great East Japan earthquake, was the fourth most powerful earthquake in the world to date (magnitude 9.0 Mw). The earthquake triggered tsunami waves and severely destroyed the eastern part of Japan. The study site is located in Norita in Miyagi prefecture, whose eastern area was seriously affected by the disaster. The study site has diverse land cover types, such as agricultural, mountain, and built-up areas. Fig. 3 shows the pre-processed ALOS/PALSAR images that were projected to UTM (Universal Transverse Mercator) Zone 54 North projection with the WGS-84 (World Geodetic System 1984). The spatial resolution of all the images was 15 m. The data were acquired on May 20, 2010, November 20, 2010, and April 7, 2011. The datasets are suitable for detecting the changes created by a tsunami as their acquisition parameters are very similar (Table 1).

Fig. 3.ALOS/PALSAR intensity images: (a) pre-event (2010.05.20); (b) pre-event (2010.11.20); (c) post-event (2011.04.07)

Table 1.Data specifications

A tsunami inundation map is the only available reference data for the study site (Fig. 4). The map was downloaded from the World Map (http://worldmap.havard.edu). Japan Society of GeoInformatics produced this map through a field survey. Although tsunami waves do not devastate all the ground features of an inundated area, most of the scattering characteristics will be changed due to the moisture, salinity, and debris of the destroyed features.

Fig. 4.Tsunami inundation map of Norita (white)

3.2 Experimental results and discussion

There are some parameters to be selected in the proposed method. The optimal values for the parameters were derived through trial and error. As for the quantization level, it is known that 8 or 16 normally does not cause the loss of texture information in SAR data (Kang et al., 2015), both values were used for GLCM, and 16 was selected as the parameter. The values with more accurate results were chosen in this paper. The window size of GLCM was 5×5, and the displacement (d) was 1 pixel. For the orientation (α), four kinds of values (0, 45, 90, 135°) were applied, and the mean of their results was extracted as the textural value.

3.2.1 Modeling statistical distributions of DI

The proposed statistical modeling method was separately applied to the reduced and enhanced DIs. To evaluate the usefulness of the proposed change detection method, the results of the proposed method were compared with the results based on the Gaussian distribution. In this paper, a g-mean (geometric mean) was used for accuracy evaluation, as in Park et al. (2013). The g-mean is the rate of accuracy of both the changed and unchanged areas (Eq. (22); Kubat et al., 1988).

Table 2 shows the estimated shape parameters of changed and unchanged distributions for the enhanced and the reduced areas in each DI. The proposed method failed to estimate the parameters of the enhanced area in the contrast DI because there were few pixels in the area. Apart from this area, the estimated shape parameters of DIs varied of which the values were from 0.753 to 3.808. It presented that the proposed change detection method included various distributions, and the DIs in the study site followed the diverse distributions.

Table 2.Estimated shape parameters

For comparison, two comparable detection maps were created using the Otsu method and an EM algorithm based on the Gaussian distribution. The EM algorithm based on the Gaussian distribution is different from the EM algorithm based on the GGD in terms of the value of shape parameter β. The shape parameter for the former was fixed to 2. The g-means of the proposed method, Otsu method, and the EM algorithm based on the Gaussian distribution were 0.6833, 0.037, and 0.6573, respectively, in terms of detecting the tsunami inundation area in the study site. The comparison with the Otsu method showed that it failed to describe the distribution of the DI. Through the Otsu method, a few pixels were extracted as the changed class. This was because two distributions of the changed and unchanged classes were highly overlapping as the tsunami changed a large area. The g-mean with the EM algorithm based on the Gaussian distribution was lower than that with the proposed method. The K-S (Kolmogorov-Smirnov) test was also performed for the statistical models from the two kinds of EM algorithm and the proposed method better described the DI than the method based on the Gaussian distribution. The K-S statistics of the models from the proposed method and the method based on the Gaussian distribution were 0.0370 and 0.0503, respectively, when modeling the reduced DI of the study site. The flexibility of the GGD contributed to the resolution of the thresholding problem when the two classes intricately overlapped with each other, when changes happened in a large area, as in the case of a tsunami. Therefore, the statistical modeling method based on GGD is useful for tsunami-induced change detection when the changed and unchanged classes were within the statistical distribution.

3.2.2 Tsunami-induced change detection

Fig. 5 represents the change detection results using the intensity and textural features. The sea surface was masked out in Fig. 5 to highlight the changes in the land surface. Although speckle noise filtering was applied to the SAR data, an effect of speckle noise still appeared in all the results. It was expected that using a bigger window size for filtering could reduce the effect, but it could also cause the loss of the detailed information. The results of the use of two intensity DIs were acceptable considering that most of the changed areas were close to the sea. In particular, reduced change was detected in the entire seashore that was directly struck by the tsunami waves. Excessive changes were detected in the agricultural areas even though some paddy fields were not affected by the tsunami waves (Fig. 5(a) and (b) dottedred- lined areas). It was also noticeable that the agricultural areas were determined to be enhanced-change areas based on the 2010.05.20 and 2011.04.07 images; on the contrary, they were determined to be reduced-change areas based on the 2010.11.20 and 2011.04.07 images. The rapid and dynamic seasonal difference in crops caused these results.

Fig. 5.Tsunami-induced change detection map: (a) intensity (2010.05.20-2011.04.07); (b) intensity (2010.11.20-2011.04.07); (c) autocorrelation; (d) contrast;(e) correlation; (f) dissimilarity; (g) energy; (h) entropy; (i) homogeneity (black: reduced change; white: enhanced change;gray: no change;solid-red-lined area: agricultural area without inundation; dotted-red-lined area: agricultural area with inundation)

Table 3 shows the g-means of tsunami-induced change detection using intensity and texture DIs. The g-means derived from the combination of the images obtained on May 20, 2010 and April 7, 2011, and the combination of the images obtained on November 20, 2010 and April 7, 2011 were 68.3% and 64.1%, respectively. The difference in the g-means was caused from the seasonal differences between the images. The former combination had less seasonal difference than the latter. The major reason for the low g-means was the detected change pixels in the agricultural area that were not reached by the tsunami waves. The g-mean value of 68.3% implied, however, that the intensity of single-pol SAR data is worth using for disaster-induced change detection considering that the highest g-mean in Park et al. (2013) was 63.4% when the reference data was a tsunami inundation area.

Table 3.g-means of tsunami-induced change mapping using intensity

The textural features were generated only from the combination of the images taken on May 20, 2010 and April 7, 2011 because its g-mean using intensity images was higher than the value obtained using the combination of intensity images taken on November 20, 2010 and April 7, 2011 (Table 3). As the highest g-mean derived from textural features was just 51.0%, the use of textural information brought about no meaningful result in terms of tsunami-induced change detection. However, textural features showed the possibility of detecting an agricultural area. The results that were obtained from autocorrelation and correlation tend not to extract changed pixels in agricultural areas (Fig. 5 (c) and (e) red-lined areas). Texture information of SAR data, in general, provides significant information about land cover and has been used for land cover mapping (Dekker, 2003; Zhu et al., 2012; Kumar et al., 2015). Zhu et al. (2012) showed that texture information improved the classification result. In this respect, it was expected that change detection results using texture information also gave the information about the land cover types.

To confirm the possibility of using the two textural features to extract agricultural areas, the results of derived from autocorrelation and correlation were combined, and a mask for agricultural areas was produced. First, the large unchanged areas were extracted from each result. The common areas between the extracted areas became the initial agricultural areas. After applying a morphological filter for smoothing the mask image, the final mask for the agricultural area was obtained. Fig. 6(a) represents the mask that was obtained from change maps using autocorrelation and correlation. The mask includes some areas with tsunami inundation, but a large percentage of agricultural areas without inundation. This implied that change detection results derived from the textural features such as autocorrelation and correlation could give the information about the land cover type, particularly about the agricultural areas. The final tsunami-induced change detection map is shown in Fig. 6(b). It was noticed that the changed areas extracted from the agricultural areas without inundation were reduced.

Fig. 6.Change detection with textural information and intensity: (a) agricultural area (white); (b) final tsunami-induced change detection map (2010.11.20-2011.04.07) (solid-lined area: agricultural area without inundation; dotted-lined area; agricultural area with inundation)

 

4. Conclusion

Although quad-pol SAR data are currently receiving much attention, single-pol SAR data are still worth using for disaster-induced change detection due to its wide imaging mode and the amount of data that have been acquired over the past decades. In this paper, single-pol SAR data were used to detect tsunami-induced changes and to evaluate their usability. For the change detection method, an automatic change detection method based on a mixture of GGDs was proposed. The usability of eight kinds of SAR factors: intensity and seven textural features derived from GLCM, was evaluated to detect tsunami-induced change using the proposed algorithm.

The experiments were performed using three ALOS/PALSAR intensity images near Natori, which was seriously affected by the tsunami waves that were generated by the 2011 Tohoku earthquake. The results showed that the proposed change detection method is more practical than the other methods based on the Gaussian distribution for detecting a tsunami-induced area. Therefore, the automatic change detection method proposed in this paper can be applied to other SAR data that have complicated statistical distributions, and will provide fine change detection results. Applying the proposed method, tsunami-induced changes were detected using intensity information, and a desirable result with a 68.3% g-mean was obtained, which was comparable with that obtained in the previous study. On the contrary, the textural features produced meaningless results in terms of tsunami-induced change detection. In terms of the land cover type, however, autocorrelation and correlation were found usable for detecting tsunami-induced change in agricultural areas. Consequently, the final change detection map was derived from the detection maps of intensity, autocorrelation, and correlation. The final change map included less pixels in agricultural areas without inundation than the initial map. Based on the experiment results, it can be concluded that single-pol SAR data have a potential for disaster-induced change detection, and that texture information can improve the change detection result when used as auxiliary data.

Eight kinds of factors were evaluated for their usability in this paper, but there are much more developed textural information and other factors such as phase information. Therefore, evaluation of the usability of other factors will be conducted for the future work. Additionally, to refine the effects of the speckle noise that appeared in the change detection map, MRF (Markov random field) will be used.