## 1.

## Introduction

Flooding events are the most common natural disaster worldwide, and their frequency may increase in the future due to global climate change;^{1} thus, flood monitoring is a national policy issue of increasing importance, requiring rapid access to accurate information that identifies changes induced by floods. Multitemporal remote sensing imagery has proven particularly useful in addressing computer-assisted change-detection applications related to flood monitoring.^{2}3.^{–}^{4} To precisely extract the flood extent information from multitemporal satellite imagery, it is necessary to carry out radiometric correction, which minimizes the unfavorable impact of radiometric differences on change detection caused by variations in imaging conditions. Two types of radiometric corrections, absolute and relative, are commonly employed to normalize remote-sensing images for the comparison of multitemporal satellite images.^{5} The absolute radiometric correction extracts the absolute reflectance of scene targets at the time of data acquisition. Most methods for absolute radiometric correction require the input of simultaneous atmospheric conditions and sensor calibration parameters, which are difficult to acquire in many cases, especially when historical data are used for change-detection analysis.^{6} Relative radiometric normalization is preferred because it does not require *in situ* atmospheric data at the time of satellite overpass.^{7}^{,}^{8} This method involves normalizing the intensities of multitemporal images, band-by-band, to a reference image selected by the analyst. Well-normalized images would appear as if they were acquired with the same sensor under similar atmospheric and illumination conditions to those of the reference image.^{9} In performing the relative radiometric normalization, it is assumed that the relationship between the radiance obtained by the sensors at two different times can be approximated with linear functions.^{10} In this method, the critical issue is determining suitable time-invariant features that can be used as the basis for normalization.^{11}^{,}^{12} A variety of relative radiometric normalization methods, such as pseudoinvariant features (PIFs), dark and bright set, simple regression, no-change determined from scattergrams, histogram matching, and multivariate alteration detection (MAD), have been investigated extensively from a theoretical and practical aspect over the last few decades.^{13}14.^{–}^{15} A new method based on PIFs, which includes automatic selection and optimization of PIFs for the radiometric normalization of multisensor images^{16}, has been introduced. A hierarchical regression method based on spectral difference has been proposed recently to reduce the radiation difference for multitemporal images, which extracts PIFs and optimizes the normalization parameters.^{17} To generate a mosaic image using multitemporal airborne thermal infrared images, a polynomial regression method was recently applied to improve the radiometric agreement between adjacent stitched images.^{18} A recent study used a parallelization method and iterative MAD to meet the demands of rapid radiometric normalization of remote sensing images for mosaicking.^{19}

One of the widely used methods is the MAD transformation because it is invariant to linear transformations of the original image intensities, which indicates that it is insensitive to differences in atmospheric conditions or sensor calibrations. For this reason, it is considered a more robust method over traditional methods for change detection.^{20} The iteratively reweighted (IR)-MAD was proposed by Nielsen to improve the robustness of the MAD transformation with the iterative updating of weights.^{21} The conceptual basis for IR-MAD is simply an iterative scheme to assign high weights on pixels that exhibit little change over time. The chi-square distribution is used to model the probability of no-change at every pixel. These probabilities are then used as weights for each iteration. The IR-MAD procedure is superior to the ordinary MAD transformation in isolating the no-change pixels suitable for use in relative radiometric normalization, particularly for data sets that exhibit large seasonal changes. However, the efficiency and robustness of IR-MAD has not yet been verified in data sets that exhibit substantial changes in land cover due to floods. The objective of this study is to extract reliable PIFs from the images acquired before and after flooding and provide an accurate radiometric normalization of a series of bitemporal very high-resolution (VHR) satellite images for flood change detection. To accomplish this, an MAD-based method that can consider the influence of open water at the study sites in the computation of the covariance matrices of the MAD transform is presented. The MAD method aims at finding a linear relationship between two sets of variables. The target image is therefore corrected by the linear relationship on the PIFs. If image pixels affected by the flood are extracted as PIFs, linear correlation information on the PIFs becomes unreliable. To accurately extract the flood-induced change information, it is necessary to carry out radiometric corrections that minimize the unfavorable impact of radiometric differences caused by image pixels affected by the flood. To address this issue, the current study uses the normalized difference water index (NDWI) to estimate open water features in the study site. In addition, the current study develops a weighting function to adaptively assign weights to the pixels based on the NDWI difference value for the radiometric normalization of multitemporal images acquired before and after flooding.

## 2.

## Image Preparation

In this study, two bitemporal images acquired by the KOMPSAT-2 satellite over the city of N’djamena, Chad, and the Atbara River, Sudan, were used to evaluate the performance and feasibility of our methodology.^{22} The topography of the two regions is relatively flat and intersected by rivers. In Chad, flooding is a frequent consequence of heavy rainfall caused by tropical cyclones. Despite serious water shortages, flash floods caused by torrential rainfall and run-off are common in Sudan. Major floods in Chad and Sudan occurred due to heavy rainfall episodes on October 12, 2012, and August 2, 2013, respectively. The technical specifications of the KOMPSAT-2 datasets are described in Table 1. In general, a radiometric difference exists between the bitemporal images acquired under different view directions at different solar hours. Such a radiometric difference is clearly observed in the vegetated area on the left bottom side of the bitemporal images of Chad in Fig. 1. There is also a large time difference between the bitemporal images for each site, as reported in Table 1.

## Table 1

KOMPSAT-2 satellite data characteristics.

Chad (site 1) | Sudan (site 2) | |||
---|---|---|---|---|

Before the flood event | After the flood event | Before the flood event | After the flood event | |

Acquisition date | June 22, 2010 | October 14, 2012 | January 9, 2012 | August 9, 2013 |

Spatial resolution | PAN: 1 m | PAN: 1 m | PAN: 1 m | PAN: 1 m |

MS: 4 m | MS: 4 m | MS: 4 m | MS: 4 m | |

Radiometric resolution | 10 bit | 10 bit | 10 bit | 10 bit |

Number of band (spectral resolution) | PAN: 1 | PAN: 1 | PAN: 1 | PAN: 1 |

MS: 4 | MS: 4 | MS: 4 | MS: 4 | |

Wavelength information | PAN: 500 to 900 nm | |||

MS 1 (blue): 450 to 520 nm | ||||

MS 2 (green): 520 to 600 nm | ||||

MS 3 (red): 630 to 690 nm | ||||

MS 4 (NIR): 760 to 900 nm | ||||

Off-nadir angle (deg) | 2 | 24 | −2 | −19 |

Processing level | Level 1R | Level 1R | Level 1R | Level 1R |

Bitemporal images just before and after the flood event are more appropriate for this application. Unfortunately, bitemporal images with short time differences could not be acquired in this experiment because these images are generally difficult to obtain due to cloud cover, long revisiting cycles of high-resolution data gathering satellites, and scarcity of data of appropriate quality. The images for each site exhibit a high proportion of changes due to significant flooding, as shown in Fig. 1. As indicted in Table 1, VHR satellite images, such as KOMPSAT-2, provide a high spatial resolution panchromatic (PAN) image and a set of multispectral (MS) images with lower spatial resolution but higher spectral resolution. To take advantage of both high spatial and spectral resolution in the change-detection process, a pan-sharpening process that merges the PAN and MS images to create a set of MS images with both a high spectral resolution and enhanced spatial resolution is required. In this study, the Gram–Schmidt adaptive (GSA) pan-sharpening method was used to generate single high-resolution MS images for both images of each study site. The GSA method provided by ENVI software comprises two steps: high-frequency spatial information is extracted from the PAN image and then injected into the resized MS images.^{23} The images were taken with different off-nadir look angles, as shown in Table 1. Therefore, it was necessary to georeference the datasets to a common coordinate system using an image registration technique. The bitemporal pan-sharpened images of each site were coregistered using the manual image-to-image registration module provided in the ENVI image processing software; the accuracy of the coregistration evaluated using 10 checkpoints gave a positional accuracy within a root mean square error of 0.5 pixels for each study site.

## 3.

## Methodology

The schematic diagram in Fig. 2 illustrates the concept and procedure of the proposed method. Our procedure for relative radiometric normalization based on the MAD transformation is conducted using the following six steps. (1) Generate the NDWI difference between the two images taken before and after the flood event. (2) Compute the weight value using the proposed weighting function defined in Eq. (6), which used the NDWI difference value and near-infrared (NIR) reflectance. (3) Apply the weight value to the bitemporal image to determine its mean and covariance matrices defined in Eq. (8). (4) Perform canonical correlation analysis (CCA) to construct the MAD variates. (5) Select pixels with no-change probability [Eq. (10)] exceeding a predefined threshold as PIFs. (6) Perform the orthogonal linear regression based on the selected PIFs to normalize the image taken after the flood band-by-band to the one taken before the flood. In this study, we compared our method with the IR-MAD method to evaluate its performance and feasibility. The quality of the normalized images generated from each method was evaluated in terms of the paired $t$-test and $F$-tests. The accuracy of flood change detection is also used as another way to examine the performance of the proposed method. Under the same condition, change vector analysis (CVA)- and MAD-based change detections are applied to the normalized images obtained by each method.

## 3.1.

### MAD and IR-MAD Transformation

The MAD transformation is an orthogonal transformation based on the CCA between two groups of variables; the transformation identifies the linear combinations that provide a set of mutually orthogonal difference images (MAD components) of decreasing variance.^{20} Let us consider two $k$-band MS images, $\mathbf{F}$ and $\mathbf{G}$, acquired in the same geographical area at two different times. We can represent the observations in different bands of the multispectral images as random vectors $\mathbf{F}=[{F}_{1},{F}_{2},\cdots ,{F}_{k}]$ and $\mathbf{G}=[{G}_{1},{G}_{2},\cdots ,{G}_{k}]$. The MAD transformation can be formulated as follows:

## (1)

$$\mathbf{U}={\mathbf{a}}^{\mathrm{T}}\mathbf{F},\phantom{\rule[-0.0ex]{2.0em}{0.0ex}}\mathbf{V}={\mathbf{b}}^{\mathrm{T}}\mathbf{G},\phantom{\rule{0ex}{0ex}}\mathbf{M}=\mathbf{U}-\mathbf{V}={\mathbf{a}}^{\mathrm{T}}\mathbf{F}-{\mathbf{b}}^{\mathrm{T}}\mathbf{G},$$## (2)

$$\sum _{fg}\sum _{gg}^{-1}\sum _{gf}\mathbf{a}={\rho}^{2}\sum _{ff}\mathbf{a}\text{\hspace{0.17em}\hspace{0.17em}},\phantom{\rule{0ex}{0ex}}\sum _{gf}\sum _{ff}^{-1}\sum _{fg}\mathbf{b}={\rho}^{2}\sum _{gg}\mathbf{b}\text{\hspace{0.17em}\hspace{0.17em}},$$^{21}

## 3.2.

### Determining Weights for the Weighted Covariance Matrices

From the above section, it can be seen that the MAD and IR-MAD transformations intrinsically project data containing the total difference between two images into uncorrelated $k$ components ${M}_{i}(i=\mathrm{1,2},\cdots ,k)$ to detect the difference between them, subject to the constraint of maintaining the total difference information as much as possible. However, these methods lead to an incorrect projection of the MAD variates in data sets with a large number of pixels in the scene that change over time, such as flood disaster data sets. The problem arises because the covariance matrix is greatly influenced by the image pixels affected by the flood; the PIF extraction result is unreliable to a considerable degree.^{24} To circumvent this problem and achieve reliable extractions of PIFs in such data sets, we devised a weighting function for the robust estimation of the covariance matrices. The basic concept for calculating weights is to assign a small weight value to pixels over open water features in calculating the covariance matrices in Eq. (2) because these pixels have a high probability of belonging to the area affected by the flood. In this study, the NDWI is employed to generate weights according to the presence of an open water feature in the image pixel. The NDWI has been developed to delineate open water features in remotely sensed digital imagery.^{25}^{,}^{26} There are two popular versions of NDWI, one using NIR and short-wave infrared (SWIR) bands proposed by Gao^{25} and the other using green and NIR bands proposed by McFeeters.^{26} In this study, we used the NDWI version proposed by McFeeters because the KOMPSAT-2 satellite does not provide SWIR band data. The NDWI uses reflectance information from the green and NIR spectral bands. The open water condition influences the interaction between these two spectral regions. The NDWI is expressed as follows:

## (3)

$$\mathrm{NDWI}=\frac{{B}_{\text{Green}}-{B}_{\mathrm{NIR}}}{{B}_{\text{Green}}+{B}_{\mathrm{NIR}}},$$Water usually has higher green reflectance than NIR reflectance. As a result, the NDWI value generally takes positive values in regions with open water within the range $[-1;1]$. In this study, we propose a weighting function that allows different weights to be assigned according to the NDWI difference value and NIR reflectance. The proposed weighting function consists of the product of two simple probability functions. In the first function, which is designed using Gaussian function, each pixel is weighted according to the difference in NDWI value. The function is defined as follows:

where $d$ is the difference in NDWI values before and after the flood at a pixel and $\sigma $ is the standard deviation that determines the height and width of the bell-shaped curve. A large standard deviation creates a bell that is short and wide while a small standard deviation creates a tall and narrow curve. This function assigns higher weights when the differences in the NDWI values are small. In the second function, which is designed using a logistic function, each pixel is weighted according to the NIR reflectance value of the image taken after the flood. The function is defined as follows: where $r$ is the NIR reflectance value of the image taken after the flood event, $k$ controls the steepness of the sigmoid curve, and ${r}_{0}$ is the midpoint of the sigmoid curve at which the curvature changes from concave to convex. In this study, this midpoint was set to the third quartile of the NIR reflectance values in the NIR image, where the third quartile is the central value between the median and the highest value of the data set. This function assigns higher weights when the NIR reflectance is high. The final weighting function made with these two functions is defined as follows:## (6)

$$\omega (d,r)={\omega}_{1}(d)\times {\omega}_{2}(r)=\mathrm{exp}(-\frac{{d}^{2}}{2\sigma})\times \frac{1}{1+\mathrm{exp}[-k(r-{r}_{0})]}.$$We performed the experiment using various values of standard deviation $\sigma $ and steepness $k$ to find the most acceptable value of these parameters. The best results were acquired when we set the standard deviation to 0.0001 and the steepness to 3 in the Chad dataset (site 1). The same parameters were applied to the Sudan dataset (site 2) to check whether these parameters are suitable for the extraction of PIFs.

These weights enter the calculation of the mean and covariance matrix as feature vectors during the MAD method procedure. For example, the covariance matrix ${\sum}_{ff}$ in Eq. (2) is simply recalculated considering the weights as follows:

## (7)

$${\overline{F}}_{i}=\frac{{\sum}_{j=1}^{N}{\omega}_{j}{F}_{ji}}{{\sum}_{j=1}^{N}{\omega}_{j}},$$## (8)

$${S}_{ik}=\frac{{\sum}_{j=1}^{N}{\omega}_{j}}{{({\sum}_{j=1}^{N}{\omega}_{j})}^{2}-{\sum}_{j=1}^{N}{\omega}_{j}^{2}}\sum _{j=1}^{N}{\omega}_{j}({F}_{ji}-{\overline{F}}_{i})({F}_{jk}-{\overline{F}}_{k}).$$If all the elements of vector $\omega $ are set to 1, the proposed method is identical to the original MAD transformation.

## 3.3.

### Relative Radiometric Normalization Using MAD Combined with NDWI

The most important step is the determination of suitable PIFs in the relative radiometric normalization because the normalization performance can vary depending on the quality and quantity of PIFs selected from the image. As described in Sec. 3.1, the MAD variates are uncorrelated (orthogonal) and invariant under affine transformations of the bitemporal images. This invariance can be exploited to determine PIFs suitable for relative radiometric normalization. To choose the PIFs using the MAD transformation, the random variable $Z$ represents the sum of the squares of the standardized MAD variates

## (9)

$$Z=\sum _{i=1}^{k}{\left(\frac{{M}_{i}}{{\sigma}_{{M}_{i}}}\right)}^{2},\phantom{\rule{0ex}{0ex}}{\sigma}_{{M}_{i}}^{2}=\mathrm{var}({M}_{i})=\mathrm{var}({U}_{i}-{V}_{i})=2(1-{\rho}_{k-i+1}),$$^{21}The calibration parameters for radiometric normalization were determined using the orthogonal linear regression and the selected PIFs.

^{24}

## 4.

## Experimental Results and Discussion

## 4.1.

### Results of Relative Radiometric Normalization

To test the performance of the proposed method, a comparison with the IR-MAD method was conducted. The PIFs obtained from each method were used to perform an orthogonal regression for relative radiometric normalization. To compare these two methods, we visually inspected the exact geometrical position of the PIFs selected in each method. The PIF extraction results obtained using these two methods are shown in Fig. 3. The pixel positions of the PIFs are shown on the flooded images with green arrows for visual inspection in Fig. 3. As shown in Fig. 3, the number of the extracted PIFs using the IR-MAD method is smaller than that of the proposed method, and the majority of extracted PIFs using the IR-MAD method are located in the area affected by flooding, which are unsuitable for radiometric normalization. The proposed method produces relatively more PIFs in the nonflooded area, with more homogeneous surface characteristics compared with the IR-MAD methods. This improvement is due to using the proposed weighting method in calculating the covariance matrices for the MAD transformation.

The selected PIFs were subsequently used to normalize the target image to the reference image using orthogonal linear regression; we designated the image taken before the flood event as the reference image and the image taken after the flood event as the target image to be normalized. Figures 4 and 5 show the results of orthogonal regression analysis using the PIFs obtained by each method for each site. As shown in Figs. 4 and 5, it is clear that the proposed method produced a better result than the IR-MAD method in terms of the linear relationship for both sites.

Figure 6 shows the radiometric normalized images generated from the orthogonal linear regression using the PIFs selected from each method. In the Chad scenes, the normalized image using the IR-MAD method provides a visually bad result, as shown in Fig. 6(a); there are clear radiometric distortions throughout all regions when compared with the proposed method. This difference is due to the quality of the extracted PIFs; most of the extracted 150 PIFs in the IR-MAD method are located in the area affected by flooding, as shown in Fig. 3(a), which are unreliable PIFs for the radiometric normalization.

## 4.2.

### Discussion Using Statistical Evaluation Method

To quantify the comparison between the proposed and IR-MAD methods, the quality of the normalized images generated from each method was evaluated in terms of the paired $t$-test and $F$-tests for equal means and variance, respectively. In general, paired $t$-test values close to zero and $F$-test values close to one indicate good matches. In both statistical tests, for a significance level of 0.05, $P$-values close to one are desirable. The null hypothesis of equal mean and variance is accepted when the $P$-values are greater than a predefined significance level, which is traditionally set to 0.05.^{24} The statistical comparisons of hold-out test pixels for the normalized images obtained from the two methods for each site are listed in Tables 2–5. The hold-out test pixels are those that are not used in the estimation of orthogonal regression parameters and are only used for the assessment of accuracy. Figure 7 shows hold-out test pixels to estimate the accuracy of each method for each site.

## Table 2

Comparison of means and variance for 44 hold-out test pixels, with paired t-tests and F-tests for equal means and variances, and normalization using the IR-MAD method for the Chad scenes.

Band 1 | Band 2 | Band 3 | Band 4 | |
---|---|---|---|---|

Reference mean | 694.522 | 706.545 | 510.523 | 564.886 |

Normalized mean | 719.765 | 756.451 | 571.468 | 574.684 |

Difference | 25.243 | 49.905 | 60.946 | 9.798 |

$t$-stat | 5.685 | 7.096 | 7.729 | 2.432 |

$P$-value | 0.000 | 0.000 | 0.000 | 0.019 |

Reference var. | 596.674 | 1822.858 | 3608.023 | 973.312 |

Normalized var. | 38.317 | 19.926 | 68.403 | 25.953 |

F-stat. | 0.064 | 0.011 | 0.019 | 0.0267 |

P-value | 0.000 | 0.000 | 0.000 | 0.000 |

## Table 3

Comparison of means and variance for 44 hold-out test pixels, with paired t-tests and F-tests for equal means and variances, and normalization using the proposed method for the Chad scenes.

Band 1 | Band 2 | Band 3 | Band 4 | |
---|---|---|---|---|

Reference mean | 694.522 | 706.545 | 510.523 | 564.886 |

Normalized mean | 697.139 | 708.777 | 513.030 | 567.389 |

Difference | 2.616 | 2.231 | 2.507 | 2.503 |

$t$-stat | 1.103 | 0.712 | 0.863 | 1.199 |

$P$-value | 0.276 | 0.480 | 0.393 | 0.237 |

Reference var. | 596.674 | 1822.858 | 3608.023 | 973.312 |

Normalized var. | 422.146 | 1476.537 | 3352.496 | 846.695 |

F-stat. | 0.708 | 0.810 | 0.929 | 0.869 |

P-value | 0.261 | 0.493 | 0.811 | 0.649 |

## Table 4

Comparison of means and variance for 326 hold-out test pixels, with paired t-tests and F-tests for equal means and variances, and normalization using the IR-MAD method for the Sudan scenes.

Band 1 | Band 2 | Band 3 | Band 4 | |
---|---|---|---|---|

Reference mean | 562.349 | 578.798 | 428.497 | 442.193 |

Normalized mean | 576.394 | 601.476 | 448.912 | 447.687 |

Difference | 14.045 | 22.678 | 20.415 | 5.494 |

t-stat | 13.785 | 15.198 | 16.58 | 4.519 |

P-value | 0.000 | 0.000 | 0.000 | 0.000 |

Reference var. | 563.656 | 1593.424 | 1815.285 | 2008.267 |

Normalized var. | 35.013 | 194.875 | 542.063 | 708.274 |

F-stat. | 0.062 | 0.122 | 0.299 | 0.353 |

P-value | 0.000 | 0.000 | 0.000 | 0.000 |

## Table 5

Comparison of means and variance for 326 hold-out test pixels, with paired t-tests and F-tests for equal means and variances, and normalization using the proposed method for the Sudan scenes.

Band 1 | Band 2 | Band 3 | Band 4 | |
---|---|---|---|---|

Reference mean | 562.349 | 578.798 | 428.497 | 442.193 |

Normalized mean | 561.959 | 577.956 | 427.051 | 442.193 |

Difference | −0.391 | −0.841 | −1.446 | −0.128 |

t-stat | −0.732 | −1.187 | −1.631 | −0.131 |

P-value | 0.465 | 0.236 | 0.104 | 0.896 |

Reference var. | 563.656 | 1593.424 | 1815.285 | 2008.267 |

Normalized var. | 657.418 | 1793.870 | 2166.284 | 2518.705 |

F-stat. | 1.1663 | 1.126 | 1.193 | 1.254 |

P-value | 0.1659 | 0.286 | 0.112 | 0.042 |

Comparing data for the Chad scenes (Tables 2 and 3), it is clear that the proposed method produced a better result than the IR-MAD method in both statistical tests; none of the $P$-values for bandwise tests for equal means and variance are acceptable in the IR-MAD method results. In the Sudan scene, as shown in Tables 4 and 5, the proposed method also generated a better result than the IR-MAD method in both statistical tests; the difference between the reference and normalized mean is much smaller than that of the IR-method. From these results, the weighting scheme in the proposed method allowed a more precise identification of PIFs prior to the relative radiometric normalization for bitemporal images exhibiting a significant amount of change due to flooding.

## 4.3.

### Results of Change Detection

Because the significance of the statistical evaluation using the paired $t$-test and $F$-tests crucially depends on the test pixels, it is difficult to conclude, using these tests alone, that the proposed method is better than the IR-MAD method. Another approach to assessing the performance of the proposed method is to compare the accuracy of flood change detection using the proposed method with that using the IR-MAD method. Under the same conditions, changes are detected in the radiometric normalized images produced from each method. Several change-detection methods have been proposed for remote sensing change detection.^{27}28.^{–}^{29} Among them, we used the CVA- and MAD-based change detection-methods for flood change detection.

The MAD approach, originally designed by Nielsen et al.,^{20} can be used directly for change detection. The thresholds for deciding between change and no-change can be set in terms of the standard deviation about the mean for each MAD component. All pixels in an MAD component ${M}_{i}$ whose intensities are within $\pm 2{\sigma}_{{M}_{i}}$, where ${\sigma}_{{M}_{i}}$ is the standard deviation, are no-change pixels. The final change pixels are obtained by applying a union operator on the change-detection results of each MAD component.

The CVA is one of the simplest and most widely used change-detection methods in the literature.^{30} The CVA is applied to the radiometric normalized image generated from each method. The basis for CVA is that a particular pixel with different values over time resides at substantially different locations in the feature space. The spectral change vector (SCV) is calculated from the vector difference of spectral feature vectors associated with pairs of corresponding pixels in two images acquired at two different times.^{31} The magnitude of SCV is used to establish a simple criterion for identifying the changed area.^{32} Due to the properties of the magnitude operator, it is possible to assert that pixels showing a magnitude higher than a given threshold value are changed, while pixels showing a magnitude lower than the threshold value are unchanged.^{33} This method performed best using Landsat TM data in a comparative evaluation of change-detection techniques for detecting areas associated with flood events.^{34}

A threshold, indicating the changed area, needs to be determined on the SCV magnitude image. Selecting an appropriate threshold value to identify change is difficult.^{35} Too low a threshold will exclude areas of change, and too high will include too many areas of change. In the literature, several threshold-selection methods have been proposed for identifying the threshold value, which separates changed areas in bitemporal satellite images.^{36} Among them, we applied the expectation maximization (EM)-based thresholding method to the SCV magnitude image to assign each image pixel to one of two opposing classes: changed and unchanged areas. The EM-based threshold algorithm requires estimates of the statistical parameters of classes, i.e., the class prior probabilities and class-conditional probabilities. The estimated class-statistical parameters are then used with the Bayes decision rule for minimum-error in the automatic determination of an optimal decision threshold.^{37} Figure 8 shows the change-detection results obtained by the MAD-based change-detection method for each site; whereas, Fig. 9 shows the change-detection results using the CVA-based change-detection method for each site. In general, the results from both change-detection methods suggest that change has been overdetected, as shown in Figs. 5 and 6.

## 4.4.

### Discussion Using Change-Detection Accuracy

To evaluate and compare the proposed algorithm, reference images for each site were produced from the original image by manually digitizing the flooded areas,^{38} as shown in Fig. 7. In the construction of the reference image, we only considered the visually salient flooded area along the river; it is challenging to visually identify all changes in urban residential districts, and we were focusing on changes due to floods. By comparing this reference image with the results from each change-detection method, we obtained a measure of detection accuracy. Comparisons between the reference and the results from each change-detection method were quantified using constructed error matrices; commission error (CE), omission error (OE), and overall accuracy (OA) were calculated to assess the whole accuracies of each result image. The OA is the sum of the correctly classified pixels divided by the total number of reference pixels.^{39}

Table 6 shows detailed quantitative results using the MAD-based change-detection method for both sites. When visually compared with the reference image from each site, there are many omissions in the results for both sites, but the OEs of the proposed method are slightly lower than those of the IR-MAD method, as shown in Table 6. From Table 6, it is also observed that the proposed method yields better accuracy in measuring OA for both sites. Table 7 represents detailed quantitative results using the CVA-based change-detection method for both sites. The CVA-based method provided a better visual and quantitative result than the MAD-based method for both sites.

## Table 6

Accuracy assessment results of the MAD-based change-detection method for each site: (OA) overall accuracy, (CE) commission error, (OE) omission error, (F) flood (in pixels), and (NF) no flood (in pixels).

Classified change | Reference change | |||||
---|---|---|---|---|---|---|

F | NF | Sum | CE (%) | |||

Chad (site 1) | IR-MAD | F | 1332 | 20,138 | 21,470 | 93.8 |

NF | 102,326 | 126,204 | 228,530 | 44.8 | ||

Sum | 103,658 | 146,342 | 250,000 | |||

OE (%) | 98.7 | 13.7 | ||||

OA = 51.01% | ||||||

Proposed method | F | 3590 | 22,313 | 25,903 | 86.14 | |

NF | 100,068 | 124,029 | 224,097 | 44.6 | ||

Sum | 103,658 | 146,342 | 250,000 | |||

OE (%) | 96.5 | 15.2 | ||||

OA = 51.04% | ||||||

Sudan (site 2) | IR-MAD | F | 44,820 | 103,149 | 147,969 | 69.7 |

NF | 189,287 | 662,744 | 852,031 | 22.2 | ||

Sum | 234,107 | 765,893 | 1,000,000 | |||

OE (%) | 80.8 | 13.4 | ||||

OA = 70.75% | ||||||

Proposed method | F | 89,421 | 93,269 | 182,690 | 51.1 | |

NF | 144,686 | 672,624 | 817,310 | 17.7 | ||

Sum | 234,107 | 765,893 | 1,000,000 | |||

OE (%) | 61.8 | 12.1 | ||||

OA = 76.20% |

## Table 7

Accuracy assessment results of the CVA-based change-detection method for each site: (OA) overall accuracy, (CE) commission error, (OE) omission error, (F) flood (in pixels), and (NF) no flood (in pixels).

Classified change | Reference change | |||||
---|---|---|---|---|---|---|

F | NF | Sum | CE (%) | |||

Chad (site 1) | IR-MAD | F | 93,079 | 34,933 | 128,012 | 27.2 |

NF | 10,579 | 111,409 | 121,988 | 8.7 | ||

Sum | 103,658 | 146,342 | 250,000 | |||

OE (%) | 10.2 | 23.8 | ||||

OA = 81.79% | ||||||

Proposed method | F | 89,839 | 27,116 | 116,955 | 23.18 | |

NF | 13,819 | 119,226 | 133,045 | 10.39 | ||

Sum | 103,658 | 146,342 | 250,000 | |||

OE (%) | 13.3 | 18.5 | ||||

OA = 83.62% | ||||||

Sudan (site 2) | IR-MAD | F | 220,160 | 353,501 | 573,661 | 61.62 |

NF | 13,947 | 412,392 | 426,339 | 3.27 | ||

Sum | 234,107 | 765,893 | 1,000,000 | |||

OE (%) | 5.95 | 46.16 | ||||

OA = 63.25% | ||||||

Proposed method | F | 195,556 | 202,251 | 397,807 | 50.84 | |

NF | 38,551 | 563,642 | 602,193 | 6.40 | ||

Sum | 234,107 | 765,893 | 1,000,000 | |||

OE (%) | 16.47 | 26.41 | ||||

OA = 75.91% |

Results obtained for both sites with both the IR-MAD and proposed methods were consistent with actual flood changes. Upon close inspection of the change-detection results using the reference images (Fig. 10), the flood extent extracted using the IR-MAD method overestimated changes in comparison to the results from the proposed method. In the results from the IR-MAD method, more pixels were identified as flooded areas than were identified by visual inspection; there was also a considerable CE, particularly in the forest and permanent water body region, relative to the proposed method. As shown in Table 7, the proposed method produced a better result than the IR-MAD method with an OA of 83.62% and 75.91%, respectively, for each site. The OA of the proposed method is higher by 1.8% and 12.6% for the two sites than the OA of the IR-MAD method. These results demonstrate the feasibility and effectiveness of employing the NDWI difference and NIR reflectance for the relative radiometric correction of remote sensing imagery to extract flooded areas. The proposed method has the advantages of not only being able to extract more precise PIFs but also performing well in differentiating the flooded area in comparison to the IR-MAD method. Our method is computationally efficient since it does not require an iterative procedure. However, there are disadvantages in that it focuses only on flood-related issues and is sensitive to the NIR band because it is designed based on the NIR band of the post-flood image.

In future work, to increase the robustness of the proposed method, we will explore new strategies to fuse flood-related information of VHR synthetic aperture radar (SAR) images when SAR images of the same area acquired during and after the flood are available.

## 5.

## Conclusions

In this study, we presented a new approach combining MAD transformation and open water features for the relative radiometric normalization of bitemporal VHR satellite imagery. The proposed method constructs a weighting function based on differences in NDWI and NIR reflectance; this function is used to estimate the covariance matrix for the MAD transformation and reliably extract PIFs from images acquired at different times. The effectiveness of this approach was verified with the experimental results, showing the extraction of PIFs from two KOMPSAT-2 VHR satellite images acquired before and after flooding. To test the performance of the proposed approach, the results were compared with those of IR-MAD-based radiometric normalization techniques. Both statistical tests and actual performances of the flood change detection were evaluated, and results demonstrated that the proposed method can extract PIFs suitable for use in relative radiometric normalization for flood change detection. Both the statistical paired $t$-test and $F$-test on hold-out test pixels also convincingly indicate that, for bitemporal scenes exhibiting a large amount of change due to flooding, the proposed method produces a better result than the IR-MAD method. Based on OA, in actual experiments of flood change detection using the MAD- and CVA-based methods, the proposed method also produced better results than the IR-MAD method. The CVA-based method produces better change-detection results than the MAD-based method for both sites. Using the CVA-based change-detection method, the OA achieved with the proposed method was 1.8% and 12.6% better than those obtained with the IR-MAD method for the two study sites. To improve the accuracy and effectiveness of the proposed method, our future research will focus on developing a strategy to utilize VHR SAR images.

## Acknowledgments

This research was supported by a grant (18RDRP-B076564-05) from the Regional Development Research Program funded by the Ministry of Land, Infrastructure, and Transport of Korean Government.

## References

## Biography

**Younggi Byun** received his MS degree and PhD in civil and environmental engineering from Seoul National University, Seoul, South Korea, in 2004 and 2011, respectively. His major research interests include UAV image processing, change-detection, and multisensor image matching and fusion.

**Dongyeob Han** received his MS degree in civil and environmental engineering and his PhD in civil, urban, and geosystem engineering from Seoul National University, Seoul, South Korea, in 1998 and 2007, respectively. His research interests include UAV image processing, laser scanning data processing, and multisensor image registration.