Mapping the 2021 October Flood Event in the Subsiding Taiyuan Basin by Multitemporal SAR Data

A flood event induced by heavy rainfall hit the Taiyuan basin in north China in early October of 2021. In this study, we map the flood event process using the multitemporal synthetic aperture radar (SAR) images acquired by Sentinel-1. First, we develop a spatiotemporal filter based on low-rank tensor approximation (STF-LRTA) for removing the speckle noise in SAR images. Next, we employ the classic log-ratio change indicator and the minimum error threshold algorithm to characterize the flood using the filtered images. Finally, we relate the flood inundation to the land subsidence in the Taiyuan basin by jointly analyzing the multitemporal SAR change detection results and interferometric SAR (InSAR) time-series measurements (preflood). The validation experiments compare the proposed filter with the Refined-Lee filter, Gamma filter, and a statistically homogeneous pixels selection based multitemporal SAR filter. The results demonstrate the effectiveness and advantage of the proposed STF-LRTA method in SAR despeckling and detail preservation, and the applicability to change scenes. The joint analyses reveal that land subsidence might be an important contributor to the flood event, and the flood recession process linearly correlates with time and subsidence magnitude.

we analyze multitemporal spaceborne synthetic aperture radar (SAR) images to restore this flood event.
In the past decades, the spatial flood extent was usually modeled as spatial feature extraction or dual-temporal change detection by one or a pair of SAR images [1], [2], [3], which, however, is not reasonable, as flood events are long-term dynamic processes. Disaster assessment should consider the spatiotemporal change of floods. The temporal evolution of the flood could be inferred by a Bayesian network using multitemporal SAR intensity, coherence, and other ground information [4]. Another Bayesian approach was applied to generate time-series probabilistic flood maps rather than binary flood information to constrain the correlation between the distribution of flooded areas and the rainfall [5]. The water persistence estimated by the Otsu method was used to evaluate the area of the flooded agricultural land [6]. These studies assess the long-term impacts of floods considering the temporal and spatial features of flood events using more than two SAR images, to achieve disaster monitoring and mitigation.
The studies above focus on advanced algorithms for change extraction but ignore the SAR speckle noise reduction. Speckle noise affects the change detection and prevents accurate extraction of flood changes [7]. Therefore, as the first step of change detection, SAR despeckling is essential to ensure the accuracy of change detection. Classical filters for single SAR images, such as Refined-Lee and Gamma filters with boxcar window [8], [9], have been used for data preprocessing, but they may blur image details, resulting in misclassification of the boundary of change patches and missing the areas with small changes. Meanwhile, single-image filters ignore the temporal features. The multitemporal SAR filters with statistically homogeneous pixels selection (SHPS) take the temporal features of SAR data stack into account and can maintain spatial details [10], [11], but they require a lot of images with unchanged ground objects when identifying homogeneous pixels, so they cannot work well in the regions with changes. A robust multitemporal SAR filter applicable for an area with changes is needed to improve the change detection accuracy. Tensor decomposition is often used in high-dimensional feature extraction. The feature of polarimetric SAR local tensor is applied to improve the classification accuracy [12]. The redundancy of spatial and polarimetric information reduces the impact of speckle noise on classification. Similarly, multitemporal SAR images with low-rank property in spatial and temporal dimensions contribute to despeckling and are not affected by the changes of time series, which, however, is rarely studied. The Taiyuan basin, hit worst by this flood event, has serious surface subsidence due to long-term excessive groundwater exploitation, according to interferometric SAR (InSAR) timeseries analysis [13], [14], [15]. The distribution of subsidence funnels and belts suggests the impacts of human activities on this region. Similarly, the subsidence of New Orleans in the U.S. is also related to groundwater withdrawal [16], [17], and it was stricken by a flood event brought by Hurricane Katrina in August 2005. The InSAR measurements in this area revealed that rapid subsidence may have led to levee failures during the storm surge [18]. However, the evolution process of the flood was rarely mentioned. Although time series InSAR analysis and multitemporal change detection are two different methods, their results can be jointly applied to disaster analysis. The correlation between land subsidence and flood event should be studied by jointly analyzing the InSAR displacement time series and the spatiotemporal change detected using both SAR amplitude and phase data.
In this study, we first propose a spatiotemporal filter based on low-rank tensor approximation (STF-LRTA) and apply it to the SAR dataset acquired by Sentinel-1. Next, we use the log-ratio change indicator and minimum error threshold algorithms, to extract flood inundation from the filtered images. A reservoir and cultivated land are selected to validate the method for simple and complicated scenes, respectively. We compare the proposed STF-LRTA with the Refined-Lee filter, Gamma filter, and an SHPS-based multitemporal SAR filter to assess the filtered images and change detection results. Finally, we investigate the evolution process of the flood event and measure the annual surface deformation velocity in the whole Taiyuan basin before the flood event, and discuss their correlation.

II. STUDY AREA AND DATASETS
The Taiyuan basin is located between the mountain ranges of Taihang and Lvliang, with an average altitude of 800-900 m (see Fig. 1). The Fen River and some small rivers run through the central part of the basin, and the secondary terraces are widely developed along the bank. Except the Fen River, the rivers in  the Taiyuan basin are relatively small and have high sediment content. The riverbed is wide and shallow. So this region is vulnerable to floods. Due to long-term excessive groundwater exploitation, many regions in this basin subside, though the velocity has reduced by up to ∼70% from 2017 to 2020 to the period of 2007 and 2010 [19]. From October 2 to 6, 2021, this area was stricken by continuous rainstorms, resulting in the water level rise of tributaries of Fen River, which are the Wuma River and Ciyao River, and flood outburst. A large area of agricultural land and bare land was flooded and become open water as the water could not be discharged timely. Since both agricultural land and bare land are permeable surface, the process of flood recession is affected by the evaporation and infiltration.
The data used in this study are listed in Table I. The ground range detected (GRD) and single look complex (SLC) data of VV polarization over the study area acquired by Sentinel-1 were collected for extracting flood change and time-series deformation, respectively. The accuracy of flood detection in VV polarization is higher than that in VH polarization [20], [21], so we selected the SAR images in VV polarization. The L2A products obtained by Sentinel-2 were processed to generate the modified normalized difference water index (MNDWI) [22] for validation. Fig. 2 shows the two test areas, the Wenyuhe Reservoir (test area A) and cultivated land (test area B), which are located in or near the Taiyuan basin. The manually labeled ground truths are from MNDWI and SAR images.

III. METHODOLOGY
The data processing consists of three steps (see Fig. 3).
Step one maps the flood changes, taking the GRD datasets spanning the flood event as input. The GRD data are geometrically coregistered and radiometrically normalized.
Step two maps the preflood land subsidence through InSAR time-series analysis of the SLC data stack. The above two steps are conducted with different datasets.
Step three jointly analyzes the flood change map and the land subsidence map to explore their relationship.

A. Multitemporal SAR Images Despeckling by STF-LRTA 1) Principle of Low-Rank Tensor Approximation (LRTA):
Tensor A ࢠ R L1×L2 × …×Ln is the high-order equivalent of vector (first-order tensor) and matrix (second-order tensor). Each order is the tensor's nth mode. And a stack of SAR amplitude images, which is a three-dimensional data cube, can be expressed as a third-order tensor X ࢠ R I×J ×K , in which I and J are the spatial dimensions, and K is the number of temporal acquisitions. It follows the basic operation rule of tensors and the solution of tensor decomposition [23], [24], [25].
Generally, a third-order tensor can be decomposed into a linear combination of R rank-1 tensors, as follows: where a r , b r , and c r are the vectors on three modes, λ r is the weight of each rank-1 tensor, and ࢪ denotes the outer product operation. In this equation, R is called the rank of tensor X. The value of R cannot be solved accurately, and it needs to be determined in advance to solve (1). Assuming that a given rank of X is equal to m and m is less than R, the low-rank decomposition of X can be solved approximately by the iterative method until the convergence condition is satisfied, and the low-rank tensor X should be are the factor matrices formed by combining the vectors of the same mode with m rank-1 tensors. The process of solvingX is called LRTA, andX contains the main information of X.
For a tensor constructed by SAR image stacks, speckle noises are considered to be independent of the signal. On the premise that rank m is limited, redundant information is compressed in a low-rank tensor, and the noise can be removed directly by LRTA, which can be written as From formulas (2) and (3), the LRTA can preserve the main features in different modes of a third-order tensor. Therefore, it can reduce the spatial noise, and also detect the features in the images with temporal changes. Thus, we propose the spatiotemporal filter, STF-LRTA, which is described below.
2) Solution for STF-LRTA: The size of the original tensor X is determined by the size and number of SAR images. Therefore, when solving the low-rank tensorX, the optimal rank m should be changed with tensor X, which affects the efficiency and robustness of the method. STF-LRTA improves the solution strategy by denoising SAR stack in small patches without determining the rank value. First, a SAR image is divided into super-pixel patches. Each patch is composed of adjacent pixels with similar brightness and textures [26], [27]. Proper segmentation considers the edge features of ground objects and rejects complicated scenes and changes. The super-pixel patches are processed separately, and the average equivalent filtering window size for the whole image is where I × J is the image size and N is the number of super-pixel patches. W can be empirically set according to image resolution and ground complexity, so N is easy to calculate. Afterward, the minimum third-order tensor is constructed for each super-pixel patch to perform the LRTA procedure. Because the tensor is small, uncomplicated, and even homogeneous, we set very low ranks, m = 1, 2, or 3, to solve the approximate tensor without considering the size of the original tensor. The change of rank at such a low level does not affect the results, so STF-LRTA does not need to find the optimal rank. For a given rank m, the alternating least squares (ALS) method is used to solve formula (2), which fixes two of factor matrices A, B, and C and solves the rest one. This operation is repeated until the convergence conditions are met. For example, if B and C are where X (1) is the mode-1 matrix of tensor X and diag(λ) represents the diagonal matrix formed by [λ 1 , λ 2 , …, λ m ]. and T is the Khatri-Rao product and the transpose operators, respectively. In such a case, formula (5) is a linear least-squares problem, which can be solved as follows: whereÂ should be normalized to obtain the new mode-1 factor matrix. Table II presents the solution process of the low-rank tensor for each patch.
Finally, after all, super-pixel patches are processed, the filtered SAR images can be obtained by stitching according to the spatial position.

B. Flood Inundation Mapping
After the SAR dataset filtering, the flood inundation is mapped by extracting change pixels. Change extraction includes difference map generation and threshold segmentation. We choose the log-ratio change indicator and the minimum error threshold algorithm for these two substeps.
The log-ratio [28] is employed to generate the time-series difference maps by comparing the reference image and each postflood image. The change indicator, (8), is suitable for the SAR amplitude images that are subject to a gamma distribution Another classical algorithm, the minimum error threshold [29] evolved from the Bayes minimum error probability criterion, is used for difference map binarization. It minimizes the where h(g) represents the probability density histogram and ε(g, T) denotes the cost function calculated by the posterior probability in the changed or unchanged situation with a given pixel value g and threshold T. It has good performance in detecting large-scale changes, which present a two-peaks distribution in the difference map histogram. We choose these two algorithms for two reasons. First, these two methods apply to different change scenes and support rapid mapping. Second, to emphasize the change detection improvement resulting from filtering, we use these traditional algorithms that have few effects on the results.

C. Time-Series InSAR Processing
In the past two decades, many time-series InSAR methods have been widely used to obtain surface subsidence. Among them, persistent scatterer (PS) and small baseline subset (SBAS) InSAR techniques estimate surface deformation based on point targets and distributed targets [30], [31], respectively. Given the rural environment of the Taiyuan basin, we select SBAS to estimate the large-scale land subsidence to ensure high spatial density of coherent points. Using a digital elevation model data, we first perform the batch processing of differential interferometry with the SLC dataset (see Fig. 4). The interferometric pairs are selected according to the given temporal and spatial baselines. Then the interferogram stack is used as the input to execute the StaMPS-SBAS method [32]. After coherent points selection, phase unwrapping, terrain residual phase, and atmospheric delay phase estimation, the time-series deformation will be obtained finally.
The result of time-series InSAR processing will be presented as the annual average deformation velocity of point targets. We will interpolate the values of discrete points into a raster image with the same resolution as the change detection map. The raster map presents the long-term ground deformation distribution before the flood, and it will be divided into several levels as the base map for the joint analysis with the dynamic flood change.

A. Evaluation of the Despeckling and Change Detection Results
The experiments were carried out according to the processing flow in Fig. 3. First, the GRD data with a size of 5000×7500 pixels were empirically divided into 240 000 super-pixel patches for LRTA with rank m = 2, so the average equivalent filtering tensor is about 13×13×5. Then, we extracted the changes between each postflood image and the reference image (preflood) to extract the flood changes, which can be expressed as a flood duration map.
To assess the performance of the proposed filter, we compared it with two classical filters, Gamma filter and Refined-Lee filter, and an advanced multitemporal SHPS-based filter [11]. The window size was set as 5×5 for the Gamma filter and the Refined-Lee filter to ensure the details can be maintained while reducing noises. The search window of 13×13 pixels was set for the SHPS-based method to identify the homogeneous pixels in five time-series GRD images. The SAR images of two test areas acquired on September 6 and November 17, close to the acquisition time of clear optical images, were selected for evaluation. As the filtered images in Fig. 5 show, all methods reduced the noises, but the Refined-Lee filter and Gamma filter got excessively smoothed results, and the SHPS-based and STF-LRTA methods can preserve more linear and point features (e.g., narrow rivers and small buildings). Two indicators were employed for the comparison. The speckle strength index (SSI) measures the noise level, and the edge sharpness index (ESI) estimates the edge preservation in the vertical and horizontal directions [33]. Their expressions are as follows: where var(), E[], and mean() are variance, expectation, and mean operators, respectively.x is the pixel value in the region of interest. x 1f and x 2f are values on either side of the edge. The smaller the SSI, the lower the noise level, and the larger the ESI, the better the edge detail enhancement. We calculated the SSI of the flood inundation region and the ESI of the linear features in the whole patch in postevent images. As listed in Table III, in both the two test areas, the Refined-Lee filter performed best in speckle noises reduction, but worst in edge information preservation as the images were excessively smoothed. The Gamma filter also smoothed most areas of the images, but it overemphasized the highlight features, so its SSI value is larger than that of the Refined-Lee filter. Both SSI and ESI suggest that the SHPS-based method has poor performance, because of the identification error of homogeneous pixels caused by the small number of available images and images with significant changes. The STF-LRTA method preserved the most edge information and ranked second in noise reduction. After denoising, we mapped the flood inundation by the methods described in Section III. Fig. 6 and Table IV present the binary change detection results and evaluation indicators, respectively. Precision and recall were calculated to indicate the false alarm and missed detection, whose expressions are where true positive (TP) is the number of changed pixels correctly classified, false positive (FP) denotes the number of unchanged pixels misclassified as changed ones, and false negative (FN) is the number of changed pixels misclassified into unchanged ones. All four methods detected the main change patterns, but have different levels of false alarm and missed detection. As shown in Table IV, the proposed method has the highest change detection precisions, 0.735 in test area A and 0.902 in test area B, which means lower noise level and fewer FP detections than the other methods. All methods have similar recall rates in test  area A, the simple test scene, because the change areas were concentrated and easy to be detected. In the complicated test scene, area B, small change patterns were dispersed, which increased the difficulty of change extraction. But our method had stable performance and extracted flood inundation areas with clear boundaries. The results have the least missed detection, confirming the obvious advantage of the proposed method.
The above experimental results indicate that the proposed STF-LRTA retains more image details while achieving good results in despeckling. And the images filtered by STF-LRTA are suitable for change detection.

B. InSAR Time-Series Measurements and Flood Dynamic Change
We employed 128 SAR SLC images spanning from February 2016 to December 2020 over the Taiyuan basin for InSAR time-series processing. The coregistration and differential interferometry were processed with the GAMMA software. And the thresholds of temporal intervals and spatial baselines were set as 180 days and 300 m, respectively, to ensure a good coherence of interferometric pairs. Then we applied 617 interferograms to generate the average deformation velocity in the line of sight direction using the StaMPS software as shown in Fig. 7.
As shown in Fig. 7, the ground deformation in the Taiyuan basin is relatively concentrated, showing the feature of largescale regional subsidence. In theory, PS method can also reconstruct the surface deformation [34]. However, the Taiyuan basin is mainly covered by distributed targets, such as bare land or farmland, which is unsuitable for the PS method. Therefore, SBAS is chosen to estimate the surface deformation of the experimental area.
Based on the results of test areas A and B, we mapped the timeseries changes of flood inundation in the Taiyuan basin using four postevent images and one pre-event image. Since there was little rainfall from October 7 to November 17 (a total of 42 days), this result can be regarded as the flood duration. We applied four grey levels to denote the duration of the flood, and superposed them on the deformation velocity map (see Fig. 8). Most regions  of the Taiyuan basin had subsidence before the flood event, with the highest subsidence velocity of more than 30 mm/year. Two concentrated subsiding regions are distributed in the northeast and southwest, where the positions of dike-break were located. In addition, the inundated areas are irregularly distributed on the whole map, relatively dense near the positions of dike-break, but have no obvious correlation with the land subsidence, which is confirmed by the statistical histograms of flood distribution on October 12 (6 days after the flood event) and November 17 (42 days after the flood event) (see Fig. 9). Then, we calculated the ratio of the flood inundation area on November 17 to that on October 12 in the subsidence regions with the deformation levels of (−Ý, −30], (−30, 20], (−20, −10], and (−10, +Ý). We found that in the area with a subsidence velocity greater than 30 mm/year, the flood persistence rate was the smallest, only 11.1%, and it increased almost linearly with the decrease of deformation velocity. This means that the water recession may be related to surface deformation.

A. Spatiotemporal Features of Images Based on Different Filters
To understand the change detection results, we analyzed the texture features of the flood inundation area in multitemporal  Fig. 10(c)-(g). In the unchanged part, all four methods reduce the speckle noise and process images smoothly. However, on the boundary of the flood inundation area outlined by dotted lines in Fig. 10(d) and (e), the Refined-Lee filter and Gamma filter smoothed the sudden change to be a gradual change, which may hinder the change extraction. Therefore, although the SAR images filtered by Refined-Lee had the optimal SSI, the fuzzy details may cause more noisy false alarms. In addition, due to the selection error of homogeneous pixels caused by flood changes, the SHPS method cannot remove much noises in the flooded area as shown in Fig. 10(f).
We also studied the temporal features of the unchanged ground object in multi-temporal SAR images, as shown in Fig. 11. We selected all five GRD SAR images over test area B [see Fig. 11(a)-(e)], and analyzed the statistical box plots of a column of unchanged pixels next to the flooded area. The results are shown in Fig. 11(f). The statistical medians of the box plots in the original SAR images indicate that the pixel values have no temporal change trend. However, the box plots of the images filtered by the Refined-Lee and Gamma methods show obvious changes from the second image. As these two methods can smooth details, ambiguous spatial features will lead to errors in temporal features extraction. And such errors will be further introduced into the process of change detection, resulting in false alarm or missed detection. The SHPS method extracts correct temporal features from the five images, and the box plots are relatively concentrated, indicating that the filtered images are smooth. But there are many outliers caused by homogeneous  pixels selection errors due to the small number of images. The box plots of STF-LRTA also show correct temporal features. The value range of the five filtered images is similar, which helps the change detection.
The above-mentioned analysis proved that the proposed STF-LRTA method has advantages in despeckling and feature extraction, and it is more suitable for multitemporal SAR change detection than the other three methods.

B. Joint Analysis of Flood Dynamic Change and Land Subsidence
We selected the dike break on the Ciyao River and Wuma River (see Fig. 12), to analyze the effect of land subsidence on dike-break. The dike-break position on the Ciyao River is located in the center of a subsidence funnel [see Fig. 12(a)]. As the red curve in Fig. 12(c) shows, the total subsidence of this region has reached 127 mm from 2016 to 2020, which may cause unstable segments of the river bank. Fig. 12(b) presents the dike break on the Wuma River. The surface subsidence magnitudes on the north bank and south bank of the river are different. The broken dike is just located near the junction of different subsidence levels. Such inconsistent deformation may damage the levee. Besides, the two dike breaks were found near the confluences of rivers. After the heavy rainfall, greater water flow struck the unstable bank and finally broke the dike.
To further explore the correlation between flood changes and subsidence, we calculated the flood recession rates of four deformation levels (see Fig. 13). The flood recession rate at temporal O is where FI is the number of flood inundation pixels. In Fig. 13, three of the four deformation levels, i.e., (−30, 20], (−20, −10], and (−10, +Ý), show that the flood recession rates decreased linearly with time. The water recession rates are linearly correlated with the subsidence velocity during the first 30 days of the flood event and tend to be similar 42 days later. This is consistent with the result in Fig. 9, confirming that the flood recession process is related to surface deformation. The possible reason is that the subsidence of the Taiyuan basin was caused by groundwater exploitation, and the flood may penetrate the underground cavity. The more serious the subsidence, the faster the penetration. As the groundwater level rose, the flood recession rate reduced and became stable. But the recession rate of the largest subsidence level (−Ý, −30] is different. It is the highest on the 18th day and 42nd day, but the lowest on the 30th day among the four defamation levels. In the regions with the subsidence level of (−Ý, −30], the water penetrates the underground cavity quickly after the rainfall, so the recession rate is as high as 0.76 on the 18th day. The left water has lower pressure and the cavity is gradually filled, so the recession rate is reduced. The low value on the 30th day may further affect the recovery of the flood recession rate on the 42nd day. To get the exact reasons we need detailed groundwater and geological data, and SAR data with a higher temporal resolution, which is unavailable in this study. But this does not affect the overall trend of flood recession and its relationship with surface subsidence.

VI. CONCLUSION
In this study, we monitored the 2021 October flood event in the Taiyuan basin, Shanxi Province, China using Sentinel-1 SAR images. First, we proposed a multitemporal SAR despeckling method, STF-LRTA, to obtain accurate spatiotemporal flood changes. Then, we employed the classical log-ratio change indicator and the minimum error threshold algorithm to extract the flood inundation features from the filtered GRD SAR images, and obtained the surface deformation velocity by the StaMPS-SBAS method from a 5-year preflood SLC dataset. Finally, we analyzed the interaction between the flood event and land subsidence by superposing the change detection results of flood inundation and time-series deformation.
Comparative experiments prove that the proposed STF-LRTA method has a good performance in reducing noises and maintaining spatiotemporal features. The images filtered by this method are quite suitable for multitemporal change detection, as the precision and the recall rate of change detection results are both improved, either in simple or complicated scenes. The joint analyses indicate that the land subsidence in the Taiyuan basin might be an important contributing factor to this flood event. The underground cavities in the subsidence area might be filled with the flood, resulting in the flood recession rate slowing down almost linearly and becoming stable within 42 days after the flood breakout. There is also a correlation between the flood recession process and surface subsidence level. The more serious the subsidence, the greater the flood recession rate, and vice versa. These findings, combined with geological data, may contribute to the emergency response to flood disasters.
Further work will focus on characterizing flood changes using multiple features of SAR images, e.g., amplitude, coherence, and polarization. And the filtering method applies to other applications.