Adjacent-Track InSAR Processing for Large-Scale Land Subsidence Monitoring in the Hebei Plain

: Large-scale land subsidence has threatened the safety of the Hebei Plain in China. For tens of thousands of square kilometers of the Hebei Plain, large-scale subsidence monitoring is still one of the most difﬁcult problems to be solved. In this paper, we employed the small baseline subset (SBAS) and NSBAS technique to monitor the land subsidence in the Hebei Plain (45,000 km 2 ). The 166 Sentinel-1A data of adjacent-track 40 and 142 collected from May 2017 to May 2019 were used to generate the average deformation velocity and deformation time-series. A novel data fusion ﬂow for the generation of land subsidence velocity of adjacent-track is presented and tested, named as the fusion of time-series interferometric synthetic aperture radar (TS-InSAR) results of adjacent-track using synthetic aperture radar amplitude images (FTASA). A cross-comparison analysis between the two tracks results and two TS-InSAR results was carried out. In addition, the deformation results were validated by leveling measurements and benchmarks on bedrock results, reaching a precision 9 mm/year. Twenty-six typical subsidence bowls were identiﬁed in Handan, Xingtai, Shijiazhuang, Hengshui, Cangzhou, and Baoding. An average annual subsidence velocity over − 79 mm/year was observed in Gaoyang County of Baoding City. Through the cause analysis of the typical subsidence bowls, the results showed that the shallow and deep groundwater funnels, three different land use types over the building construction, industrial area, and dense residential area, and faults had high spatial correlation related to land subsidence bowls.


Introduction
Groundwater depletion and related land subsidence is a severe problem in many countries around the world. Several major cities worldwide suffer from aquifer depletion and land subsidence [1] due to overexploitation of groundwater resources [2]. Land subsidence may be one of the most obvious environmental impacts of groundwater extraction [3]. Globally, the demand for fresh water is the main cause of this phenomenon because of the rapid population growth [4]. Industrial fresh water is usually supplied directly by pumping from aquifer systems, which leads to compaction of the compressible aquifer [5]. It is easy to cause damage to the surface and underground infrastructure and increase the risk of flood in coastal cities [5,6].
Hebei Plain has experienced serious land subsidence (see Figure 1) and is an area of serious groundwater depletion in China [7][8][9][10][11]. The amount of water resources in the Hebei Plain is less than 1% of China, while the cultivated area accounts for 5.28% [12]. Hebei Plain covers more than 80% of towns and more than 70% of the population of Hebei Province, and groundwater extraction accounts for 80%~85% of the total water Plain covers more than 80% of towns and more than 70% of the population of Hebei Province, and groundwater extraction accounts for 80%~85% of the total water resources [13]. Since 1970, the groundwater level has been falling under the influence of natural and manmade factors in the Hebei Plain [14]. Due to the overexploitation of groundwater and the high population density, it has caused large-scale land subsidence and other geological hazards [15]. In addition to groundwater, the mining of liquid mineral resources and increased surface loads have also aggravated land subsidence [11,16,17]. Interferometric synthetic aperture radar (InSAR), which appeared in the late 1980s [18], has been widely used to monitor ground deformation information [19,20], but it is extremely vulnerable to temporal and spatial decorrelation and atmospheric phase delay [21]. The TS-InSAR technology, represented by SBAS [22], permanent scatters InSAR (PS-InSAR) [23,24], can monitor the surface centimeter and even millimeter-level slow deformation for a long time. TS-InSAR has high accuracy, a wide monitoring range, and low cost compared with earlier geological disaster monitoring methods, which have their own advantages. Multi-track data cannot be used directly to produce interferograms because of different datum and geometry. The research of 3-D and 2-D deformation monitoring using multi-source and multi-track synthetic aperture radar (SAR) data has become a hot spot [25][26][27][28], while there are relatively few studies the fusion of one-dimensional multi-track SAR data. In 2007, for the first time, Ketelaar et al. proposed a mathematical framework for the datum connection of multiple independent overlapping tracks, and discussed the SAR data fusion method of an adjacent-track under the radar coordinate system [29]. Since then, research studying the fusion of one-dimensional multi-track results has gradually begun [30][31][32][33][34][35][36][37][38]. Some methods aim to increase the spatiotemporal sampling rate or to connect time-series results having time-gapped through multi-track monitoring [32][33][34][35]. Others have only considered the offset between two track deformation rates to merge multitrack TS-InSAR results [36,37]. The time-series analysis of multi-track interferograms may cause the loss of the transient component of the deformation time-series [38].
Some studies have used TS-InSAR to monitor land subsidence in different cities in the Hebei Plain [39][40][41], but the first full coverage monitoring was completed by Zhang Yonghong et al. in 2016 [42]. This research addressed several challenges in full coverage Multi-track data cannot be used directly to produce interferograms because of different datum and geometry. The research of 3-D and 2-D deformation monitoring using multisource and multi-track synthetic aperture radar (SAR) data has become a hot spot [25][26][27][28], while there are relatively few studies the fusion of one-dimensional multi-track SAR data. In 2007, for the first time, Ketelaar et al. proposed a mathematical framework for the datum connection of multiple independent overlapping tracks, and discussed the SAR data fusion method of an adjacent-track under the radar coordinate system [29]. Since then, research studying the fusion of one-dimensional multi-track results has gradually begun [30][31][32][33][34][35][36][37][38]. Some methods aim to increase the spatiotemporal sampling rate or to connect time-series results having time-gapped through multi-track monitoring [32][33][34][35]. Others have only considered the offset between two track deformation rates to merge multi-track TS-InSAR results [36,37]. The time-series analysis of multi-track interferograms may cause the loss of the transient component of the deformation time-series [38].
Some studies have used TS-InSAR to monitor land subsidence in different cities in the Hebei Plain [39][40][41], but the first full coverage monitoring was completed by Zhang Yonghong et al. in 2016 [42]. This research addressed several challenges in full coverage and long-term (from 1992 to 2014) land subsidence monitoring, but only the offset between the two tracks is considered to merge the deformation rate maps. Recently, some studies have monitored the subsidence of the North China Plain [16,43,44]. However, this article focuses on the study of the fusion of adjacent-track TS-InSAR results and qualitatively analyzes the causes of subsidence from five aspects in the Hebei Plain. Therefore, based on the method proposed by Sergey Samsonov [45], this paper proposes a novel data fusion flow for the generation of land subsidence velocity of adjacent-track to merge the TS-InSAR results of two adjacent-tracks with overlapping regions. In addition, this method does not need to interpolate and correct the persistent scatterer (PS) point, which further reduces the error and improves the accuracy.
The main objectives of the paper are as follows: • Monitor land subsidence in the Hebei Plain (45,000 km 2 ) by using two TS-InSAR techniques with 166 Sentinel-1A images of two adjacent-track from May 2017 to May 2019; • Propose a novel data fusion flow for the generation of land subsidence vertical velocity of adjacent-track; • Analyze the spatial and temporal evolution characteristics of subsidence features in the Hebei Plain and qualitatively analyze the causes of typical land subsidence bowls' response to groundwater funnels, three land use types, and faults.
The rest of the paper is organized as follows: In Section 2, we introduce the study area and SAR data. In Section 3, TS-InSAR technologies are briefly introduced, and we propose a novel data fusion flow for the generation of the land subsidence vertical velocity of the adjacent-track. Then, data processing is briefly explained. In Section 4, several results and validations are presented. Finally, the discussion and conclusions of this study are summarized in Sections 5 and 6, respectively.

Study Area
The Hebei Plain is located in the Hebei Province of China and forms the northern part of the North China Plain. The study area is outlined by a red border in Figure 2, covering most of the Hebei Plain and involving seven cities, namely Handan, Xingtai, Shijiazhuang, Hengshui, Cangzhou, Baoding, and Langfang. The study area is approximately 45,000 km 2 (114.280-117.610E, 36.492-38.911N). The Hebei Plain has a flat topography, and the terrain is high in the northwest and low in the southeast. The ground elevation changes from 100 m in front of the mountain to 2-3 m along the Bohai Sea. There are two river basins flowing through the Hebei Plain, namely the Haihe River and the Luanhe River. The Hebei Plain belongs to a continental monsoon climate. The annual average temperature is 10-25 • C. The average annual precipitation from 2018 to 2019 is about 475 mm in Hebei Province [46,47]. The annual precipitation is mainly concentrated in summer, and there is less in spring than in winter. The total amount of water resources were 4065 and 9119 million cubic meters in 2018 and 2019, respectively, which are less than the multi-year average water resources in Hebei Province [46,47]. The shallow groundwater level in the Hebei Plain dropped by an average of 0.13 and 0.58 m between 2018 and 2019, respectively [46,47]. The detailed information of groundwater level changes in 2018 and 2019 are shown in Figure 3 and Table 1. During 2018 and 2019, 55% of the water consumption in Heibei Province came from groundwater extraction [46,47].

Data
In total, 166 scenes of ascending Sentinel-1A [48] TOPS SAR images (C-band) along two adjacent-track four frames from May 2017 to May 2019 were used to monitor subsidence. As shown in Table 2 and Figure 2, tracks 40 and 142 each had two frames. Track 40 (T40) included 86 scene images, and track 142 (T142) included 80 scene images. The dataset details are shown in Tables 2 and 3, and Figure 4. As can be seen from Figure 4, there were time-lags between the black lines, and the time was between 06/2017 and 09/2017 as well as 06/2018 and 09/2018. This is mainly due to the poor coherence of interferograms in summer. The one arc-second SRTM DEM (30*30 m) provided by the National Aeronautics and Space Administration (NASA) was used to remove topographic phases. Precise orbit determination (POD) [49] data were used for the orbital refinement and phase re-flattening.

Data
In total, 166 scenes of ascending Sentinel-1A [48] TOPS SAR images (C-band) along two adjacent-track four frames from May 2017 to May 2019 were used to monitor subsidence. As shown in Table 2 and Figure 2, tracks 40 and 142 each had two frames. Track 40 (T40) included 86 scene images, and track 142 (T142) included 80 scene images. The dataset details are shown in Tables 2 and 3, and Figure 4. As can be seen from Figure 4, there were time-lags between the black lines, and the time was between 06/2017 and 09/2017 as well as 06/2018 and 09/2018. This is mainly due to the poor coherence of interferograms in summer. The one arc-second SRTM DEM (30*30 m) provided by the National Aeronautics and Space Administration (NASA) was used to remove topographic phases. Precise orbit determination (POD) [49] data were used for the orbital refinement and phase re-flattening.

Methodology
Section 3.1 provides a brief introduction to the basic theories of two TS-InSAR technologies, the SBAS and NSBAS [50][51][52] technologies. Section 3.2 introduces a novel data fusion flow for the generation of adjacent-track land subsidence velocity. Section 3.3 briefly introduces the data processing.

TS-InSAR Technologies
Here, we briefly describe the theories of the TS-InSAR techniques. A detailed description of the SBAS algorithm is provided by Berardino et al. [22]. First, the appropriate combination of differential interferograms are generated based on small temporal and spatial baseline threshold values. The least square method and singular value decomposition

Methodology
Section 3.1 provides a brief introduction to the basic theories of two TS-InSAR technologies, the SBAS and NSBAS [50][51][52] technologies. Section 3.2 introduces a novel data fusion flow for the generation of adjacent-track land subsidence velocity. Section 3.3 briefly introduces the data processing.

TS-InSAR Technologies
Here, we briefly describe the theories of the TS-InSAR techniques. A detailed description of the SBAS algorithm is provided by Berardino et al. [22]. First, the appropriate combination of differential interferograms are generated based on small temporal and spatial baseline threshold values. The least square method and singular value decomposition (SVD) method are used to calculate the line-of-sight (LOS) direction average deformation velocity and deformation time-series. For a given pixel l, there are N l interferograms and M l SAR images. There is a linear equation [22,[50][51][52]: where d l is the phase vector of N l interferograms for pixel l; m l is a vector include M l − 1 phase delay increments for pixel l; G l is a N l × M l − 1 matrix composed of zeros and ones, constructed based on the phase of an interferogram; and φ l ij is the sum of the successive phase delays between SAR images i and j: φ l ij = j−1 ∑ k=i m l k [50]. For a relatively large number of inversion pixels, G T l G l is singular. Use the SVD method instead of least square method to solve d l = G l m l .
We also followed NSBAS as described by Lopez-Quiroz et al. [50] to build the deformation time-series and deformation velocity maps. As time went by, the SBAS showed a clearly linear deformation, which caused a deviation in the time and thus a deviation in the deformation rate. In order to solve this problem, constraints were added to the TS-InSAR inversion instead of using the SVD method. The cumulative phase delay φ l k for pixel l at time t k is defined by the following formula [50][51][52]: where m l k is the phase delay increment of the two successive SAR images at pixel l, and the following constraints are added in the inversion [50][51][52]: where e l B k ⊥ represents the DEM error phase related to the vertical baseline B k ⊥ of each acquisition. We then solve, by least square inversion, the system d c = G c m c [50]. The weight, γ, scales the additional matrix. If the SBAS network is complete, this method does not affect the inversion and is suitable for SAR data. If the SBAS network is incomplete and disconnected subsets appear, this method links these subsets [50].

Fusion of TS-InSAR Results of Adjacent-Track Using SAR Amplitude Images (FTASA)
The single-track SAR data do not cover the Hebei Plain completely, which limits the large-scale subsidence monitoring, therefore, multi-track SAR data processing has come into being. However, there are at least three problems that need to be solved in the fusion of the adjacent-track TS-InSAR results. First, due to the different incidence angles of the adjacent-track, the imaging geometries are different. Second, the unwrapping reference points selected by TS-InSAR results from adjacent-track are inconsistent. Third, there is no connection between adjacent-track when we try to merge the adjacent-track results.
This paper proposes a fusion processing flow of adjacent-track TS-InSAR results based on SAR amplitude images ( Figure 5). First, vertical deformation of every single track is retrieved. Second, adjacent-track deformation results are registered based on the SAR amplitude image. Third, the datum is unified by calculating the offset of the adjacent-track. Finally, adjacent-track deformation time-series fusion is solved by the SVD method. The detailed description of fusion processing flow is as follows:  Step 1: Vertical deformation calculation InSAR is a side-view observation technique, and its direct observation result is LOS deformation. The LOS deformation observed by InSAR can be expressed as Equation (6): sin cos sin sin cos where e d , n d , u d represent the east, north, and vertical deformations, respectively. ϑ is the incident angle, and ϕ is the angle between the satellite heading and north in a clockwise direction. Different tracks have different incidence angles. To eliminate this difference, we converted merely to vertical deformation by assuming the east and north land subsidence of Hebei is negligible [42][43][44]. Vertical deformation can be generated with Equation (7). / cos Step 2: Registration of deformation results of adjacent-track based on SAR amplitude image However, due to the different incident angles and other factors of the adjacent-track, the encoded deformation results of the adjacent-track are spatially deviated. Generally, remote sensing images are registered using the texture information of ground objects. However, the deformation results had no texture information, which makes it difficult to Step 1: Vertical deformation calculation InSAR is a side-view observation technique, and its direct observation result is LOS deformation. The LOS deformation observed by InSAR can be expressed as Equation (6): where d e , d n , d u represent the east, north, and vertical deformations, respectively. ϑ is the incident angle, and ϕ is the angle between the satellite heading and north in a clockwise direction. Different tracks have different incidence angles. To eliminate this difference, we converted merely to vertical deformation by assuming the east and north land subsidence of Hebei is negligible [42][43][44]. Vertical deformation can be generated with Equation (7).
Step 2: Registration of deformation results of adjacent-track based on SAR amplitude image However, due to the different incident angles and other factors of the adjacent-track, the encoded deformation results of the adjacent-track are spatially deviated. Generally, remote sensing images are registered using the texture information of ground objects. However, the deformation results had no texture information, which makes it difficult to register. Therefore, we propose a registration method of adjacent-track deformation results based on SAR amplitude images. The affine transformation method is used for registration with Equations (8) and (9). This method first registers the encoded SAR amplitude images of adjacent-tracks within one pixel, and the registration transfer function is obtained to register deformation results of the adjacent-track. Pre-registration and post-registration comparison images of tracks 40 and 142 are shown in Figure 6.
where X mas and Y mas are the coordinates of the master image; X sla and Y sla are the coordinates of the slave image; a ij and b ij are the polynomial coefficients to be solved, and N is the degree of the polynomial. Linear registration was used in this paper. register. Therefore, we propose a registration method of adjacent-track deformation results based on SAR amplitude images. The affine transformation method is used for registration with Equations (8) and (9). This method first registers the encoded SAR amplitude images of adjacent-tracks within one pixel, and the registration transfer function is obtained to register deformation results of the adjacent-track. Pre-registration and postregistration comparison images of tracks 40 and 142 are shown in Figure 6.
where mas X and mas Y are the coordinates of the master image; sla X and sla Y are the coordinates of the slave image; ij a and ij b are the polynomial coefficients to be solved, and N is the degree of the polynomial. Linear registration was used in this paper. Step 3: Datum unification of adjacent-track The reference points of the adjacent-track are usually difficult to keep consistent. In addition, the positions of the PS point of adjacent-track overlapping area are also not consistent. Here, we propose a unified reference method. First, filter the PS points using the overall coherence image for the single track. Second, extract the corresponding PS points of the adjacent-track based on geographic coordinates. Third, calculate the datum offset using Equation (10) [42] by analyzing the distribution of the homonymous point difference of the adjacent-track. The advantage of this method is that there is no need to interpolate the PS point of the adjacent-track. The offset is added to the adjacent-track to unify the datum.
Step 4: Fusion deformation time-series of adjacent-track Suppose that ( )  Step 3: Datum unification of adjacent-track The reference points of the adjacent-track are usually difficult to keep consistent. In addition, the positions of the PS point of adjacent-track overlapping area are also not consistent. Here, we propose a unified reference method. First, filter the PS points using the overall coherence image for the single track. Second, extract the corresponding PS points of the adjacent-track based on geographic coordinates. Third, calculate the datum offset using Equation (10) [42] by analyzing the distribution of the homonymous point difference of the adjacent-track. The advantage of this method is that there is no need to interpolate the PS point of the adjacent-track. The offset is added to the adjacent-track to unify the datum.
Step 4: Fusion deformation time-series of adjacent-track Suppose that d 1 are the vertical deformation time-series of an arbitrary pixel in two tracks overlapping, respectively. Use the SVD method to obtain the merged time-series results with Equation (11), and the schematic Remote Sens. 2021, 13, 795 9 of 23 diagram of the fusion is shown in Figure 7. Purple font T 1 0 , T 1 1 , · · · , T 1 N 1 are the times of track 1, and red font T 2 0 , T 2 1 , · · · , T 2 N 2 is the times of track 2, respectively. where is the matrix composed of time intervals, and D N 1 +N 2 +1 is the unknown merged vertical deformation velocity time-series, and Φ N 1 +N 2 is the adjacenttrack vertical deformation time-series combinations.
is the matrix composed of time intervals, and

Data Processing
The SBAS and NSBAS technique are applied to calculating the land deformation time-series using the four frames of tracks 40 and 142 of the 166 Sentinel-1A images by GAMMA and GIAnT [53,54]. Each SAR image is accurately registered to the master image. A total of 646 and 600 Sentinel-1A interferograms were generated for tracks 40 and 142, considering a maximum temporal and perpendicular baseline value of 360 days and 220 m, respectively. Concerning the interferogram coherence, 182 and 202 interferograms of tracks 40 and 142 were selected, respectively. Differential SAR interferograms were generated by using an SRTM DEM 1″. Furthermore, we also performed a 10 × 2 multilook operation for every S1A single interferogram with a 30 m × 30 m resolution. Each generated interferogram was first noise-filtered, and then unwrapped. The unwrapped phase signals were also calculated to the same reference pixel of each track. Subsequently, the LOS land deformation time-series were obtained from unwrapped DInSAR interferograms by using the SBAS and NSBAS technique. The atmospheric phase was filtered out from the obtained land deformation time-series. Moreover, assuming a simple linear relationship with topography, the stratified atmospheric phase contributions can be estimated from the data itself using a multi-scale approach [55]. We converted the LOS deformation time-series into the vertical time-series. We followed the process in Section 3.2 to get the merged subsidence velocity. First, the deformation time-series of tracks 40 and 142 were registered based on the SAR amplitude image. Second, the deformation timeseries of tracks 40 and 142 eliminated the offset 0.9 mm. Finally, the deformation timeseries of tracks 40 and 142 were merged using the SVD method.

Data Processing
The SBAS and NSBAS technique are applied to calculating the land deformation timeseries using the four frames of tracks 40 and 142 of the 166 Sentinel-1A images by GAMMA and GIAnT [53,54]. Each SAR image is accurately registered to the master image. A total of 646 and 600 Sentinel-1A interferograms were generated for tracks 40 and 142, considering a maximum temporal and perpendicular baseline value of 360 days and 220 m, respectively. Concerning the interferogram coherence, 182 and 202 interferograms of tracks 40 and 142 were selected, respectively. Differential SAR interferograms were generated by using an SRTM DEM 1". Furthermore, we also performed a 10 × 2 multilook operation for every S1A single interferogram with a 30 m × 30 m resolution. Each generated interferogram was first noise-filtered, and then unwrapped. The unwrapped phase signals were also calculated to the same reference pixel of each track. Subsequently, the LOS land deformation time-series were obtained from unwrapped DInSAR interferograms by using the SBAS and NSBAS technique. The atmospheric phase was filtered out from the obtained land deformation time-series. Moreover, assuming a simple linear relationship with topography, the stratified atmospheric phase contributions can be estimated from the data itself using a multi-scale approach [55]. We converted the LOS deformation time-series into the vertical time-series. We followed the process in Section 3.2 to get the merged subsidence velocity. First, the deformation time-series of tracks 40 and 142 were registered based on the SAR amplitude image. Second, the deformation time-series of tracks 40 and 142 eliminated the offset 0.9 mm. Finally, the deformation time-series of tracks 40 and 142 were merged using the SVD method.

Fusion of Adjacent-Track PS Results
The vertical deformation velocity maps of tracks 40 and 142 are shown in Figure 8. The black line in Figure 8 is the overlapping area of tracks 40 and 142. It is clear to see that the NSBAS technique can generally provide more coherent pixels than SBAS techniques. Meanwhile, the NSBAS technique can provide entirely consistent results with SBAS techniques, despite some differences in the detailed characteristics. Therefore, we used the results of NSBAS for subsequent analysis. From Figure 8c,d, it can be clearly seen that the deformation trend of the overlapping area of tracks 40 and 142 were basically consistent. The subsidence velocity map of Hebei Plain was generated as shown in Figure 9.  Figure 8 is the overlapping area of tracks 40 and 142. It is clear to se the NSBAS technique can generally provide more coherent pixels than SBAS techn Meanwhile, the NSBAS technique can provide entirely consistent results with SBAS niques, despite some differences in the detailed characteristics. Therefore, we use results of NSBAS for subsequent analysis. From Figure 8c,d, it can be clearly seen th deformation trend of the overlapping area of tracks 40 and 142 were basically cons The subsidence velocity map of Hebei Plain was generated as shown in Figure 9.   From Figure 9, we can see that most regions of Hebei Plain were relatively stable. However, there exist about 26 apparent land subsidence bowls (labeled with letters from d0 to d25), and they tend to be connected in the administrative area. The maximum subsidence area was located in region d22 in Gaoyang County, where the deformation rate was up to −79 mm/year. As shown in Figure 10, six land subsidence bowls were found in Handan. Seven land subsidence bowls were found in Xingtai. Two land subsidence bowls were found in Shijiazhuang. Six land subsidence bowls were found in Hengshui. Three land subsidence bowls were found in Cangzhou, and two land subsidence bowls were found in Baoding. The detailed locations of 26 subsidence bowls are shown in Table 4.  From Figure 9, we can see that most regions of Hebei Plain were relatively stable. However, there exist about 26 apparent land subsidence bowls (labeled with letters from d0 to d25), and they tend to be connected in the administrative area. The maximum subsidence area was located in region d22 in Gaoyang County, where the deformation rate was up to −79 mm/year. As shown in Figure 10, six land subsidence bowls were found in Handan. Seven land subsidence bowls were found in Xingtai. Two land subsidence bowls were found in Shijiazhuang. Six land subsidence bowls were found in Hengshui. Three land subsidence bowls were found in Cangzhou, and two land subsidence bowls were found in Baoding. The detailed locations of 26 subsidence bowls are shown in Table 4. From Figure 9, we can see that most regions of Hebei Plain were relatively stabl However, there exist about 26 apparent land subsidence bowls (labeled with letters fro d0 to d25), and they tend to be connected in the administrative area. The maximum su sidence area was located in region d22 in Gaoyang County, where the deformation ra was up to −79 mm/year. As shown in Figure 10, six land subsidence bowls were found Handan. Seven land subsidence bowls were found in Xingtai. Two land subsidence bow were found in Shijiazhuang. Six land subsidence bowls were found in Hengshui. Thr land subsidence bowls were found in Cangzhou, and two land subsidence bowls we found in Baoding. The detailed locations of 26 subsidence bowls are shown in Table 4.

Validation of the InSAR Results
In order to evaluate the accuracy of FTASA, we undertook a comparative analysis between the merged result acquired by only considering the offset (FO) and the merged result acquired by FTASA. Two merged results were also compared with the leveling results. The overlap area of tracks 40 and 142 had only two benchmarks. As shown in Figure 11, the average differences between the leveling results and FTASA, FO at the two benchmarks were 2 mm/year and 10 mm/year, respectively. This shows that FTASA has a higher accuracy than FO. We also plotted the difference between the two results as a histogram (see Figure 12  To quantitatively verify the accuracy of land subsidence monitoring using 1A data by NSBAS, not only was a comparative analysis of the differences bet NSBAS, leveling results, and benchmark on the bedrock results performed, but t 40 and 142 deformation results of overlapping were also comprehensively a Given the limited leveling data and benchmark on the bedrock data, the positio benchmark and the InSAR results were difficult to coincide. We extracted the n point corresponding to the benchmark. The leveling results were acquired from M to May 2019 and the two benchmarks on the bedrock results were acquired from  To quantitatively verify the accuracy of land subsidence monitoring using Sentinel-1A data by NSBAS, not only was a comparative analysis of the differences between the NSBAS, leveling results, and benchmark on the bedrock results performed, but the tracks 40 and 142 deformation results of overlapping were also comprehensively analyzed. Given the limited leveling data and benchmark on the bedrock data, the positions of the benchmark and the InSAR results were difficult to coincide. We extracted the nearest PS point corresponding to the benchmark. The leveling results were acquired from May 2017 to May 2019 and the two benchmarks on the bedrock results were acquired from December 2017 to June 2018. From Figure 13, we can see that the subsidence velocities obtained by InSAR were generally consistent with the leveling results and benchmark on the bedrock results. The overall RMSE between the leveling results and the InSAR results was 9 mm/year. The maximum and minimum errors were 17.5 mm/year and 1 mm/year, respectively. Overall, the InSAR results were highly correlated with the leveling results. To quantitatively verify the accuracy of land subsidence monitoring using Sentinel-1A data by NSBAS, not only was a comparative analysis of the differences between the NSBAS, leveling results, and benchmark on the bedrock results performed, but the tracks 40 and 142 deformation results of overlapping were also comprehensively analyzed. Given the limited leveling data and benchmark on the bedrock data, the positions of the benchmark and the InSAR results were difficult to coincide. We extracted the nearest PS point corresponding to the benchmark. The leveling results were acquired from May 2017 to May 2019 and the two benchmarks on the bedrock results were acquired from December 2017 to June 2018. From Figure 13, we can see that the subsidence velocities obtained by InSAR were generally consistent with the leveling results and benchmark on the bedrock results. The overall RMSE between the leveling results and the InSAR results was 9 mm/year. The maximum and minimum errors were 17.5 mm/year and 1 mm/year, respectively. Overall, the InSAR results were highly correlated with the leveling results. The deformation results of tracks 40 and 142 within the overlapped region were utilized. The histogram of the vertical deformation velocity difference within the overlapped region between tracks 40 and 142 is shown in Figure 14a. Statistical analysis showed that the 96% vertical deformation velocity difference of overlapping regions was within [−13, 13] mm/year, basically conforming to the normal distribution. The standard deviation was 7.37 mm/year. Figure 14b shows the linear fitting result of the vertical deformation velocity between tracks 40 and 142. The overlapping area of tracks 40 and 142 was approximately 14,000 km 2 , with a total of 3,400,000 PS points. The linear fitting formula of tracks 40 and 142 was T40 = 0.96 × T142 + 0.60. The correlation coefficient was 0.89. We also selected four points (see Figure 9) in the overlapping area of tracks 40 and 142, and their deformation time-series are shown in Figure 15. It is clear that the time-series of tracks 40 and 142 were highly consistent. All these prove that the vertical deformation velocity of the two tracks had a good correlation. The deformation results of tracks 40 and 142 within the overlapped region were utilized. The histogram of the vertical deformation velocity difference within the overlapped region between tracks 40 and 142 is shown in Figure 14a. Statistical analysis showed that the 96% vertical deformation velocity difference of overlapping regions was within [− 13,13] mm/year, basically conforming to the normal distribution. The standard deviation was 7.37 mm/year. Figure 14b shows the linear fitting result of the vertical deformation velocity between tracks 40 and 142. The overlapping area of tracks 40 and 142 was approximately 14,000 km 2 , with a total of 3,400,000 PS points. The linear fitting formula of tracks 40 and 142 was T40 = 0.96 × T142 + 0.60. The correlation coefficient was 0.89. We also selected four points (see Figure 9) in the overlapping area of tracks 40 and 142, and their deformation time-series are shown in Figure 15. It is clear that the time-series of tracks 40 and 142 were highly consistent. All these prove that the vertical deformation velocity of the two tracks had a good correlation. The deformation results of tracks 40 and 142 within the overlapped region were utilized. The histogram of the vertical deformation velocity difference within the overlapped region between tracks 40 and 142 is shown in Figure 14a. Statistical analysis showed that the 96% vertical deformation velocity difference of overlapping regions was within [− 13,13] mm/year, basically conforming to the normal distribution. The standard deviation was 7.37 mm/year. Figure 14b shows the linear fitting result of the vertical deformation velocity between tracks 40 and 142. The overlapping area of tracks 40 and 142 was approximately 14,000 km 2 , with a total of 3,400,000 PS points. The linear fitting formula of tracks 40 and 142 was T40 = 0.96 × T142 + 0.60. The correlation coefficient was 0.89. We also selected four points (see Figure 9) in the overlapping area of tracks 40 and 142, and their deformation time-series are shown in Figure 15. It is clear that the time-series of tracks 40 and 142 were highly consistent. All these prove that the vertical deformation velocity of the two tracks had a good correlation.

Groundwater Funnels and Land Subsidence
Shallow groundwater funnels like the Ningbolong funnel and Gaoliqing-Suning funnel, and deep groundwater funnels such as the Jizaoheng funnel and Nangong funnel are relatively large groundwater funnels in the Hebei Plain. The details of the groundwater funnels are shown in Table 5. All the groundwater data were obtained from the Department of Water Resources of Hebei Province [46,47]. Figure 16 shows that the land subsidence bowls have a strong connection with the shallow or deep groundwater funnels. The subsidence bowl d10 (see Figure 16a)

Groundwater Funnels and Land Subsidence
Shallow groundwater funnels like the Ningbolong funnel and Gaoliqing-Suning funnel, and deep groundwater funnels such as the Jizaoheng funnel and Nangong funnel are relatively large groundwater funnels in the Hebei Plain. The details of the groundwater funnels are shown in Table 5. All the groundwater data were obtained from the Department of Water Resources of Hebei Province [46,47]. Figure 16 shows that the land subsidence bowls have a strong connection with the shallow or deep groundwater funnels. The subsidence bowl d10 (see Figure 16a) Figure 16c, subsidence bowls d9 and d25 were within the Nangong funnel, and the center of the Nangong funnel was located in Jiaowang in Nangong City. In 2018, the funnel center was buried at a depth of 99.20 m, and the groundwater depth contour was 90 m surrounded the funnel area with a 664 km 2 area. The burial depth of the funnel center was 108.70 m in 2019, and the 95 m groundwater depth contour enclosing the funnel area was 947 km 2 . All in all, we found that there were six ground subsidence bowls in the shallow groundwater funnel, and two ground subsidence bowls in the deep groundwater funnel. There was a positive correlation between eight ground subsidence bowls and three groundwater funnels. 947 km 2 . All in all, we found that there were six ground subsidence bowls in the shallow groundwater funnel, and two ground subsidence bowls in the deep groundwater funnel. There was a positive correlation between eight ground subsidence bowls and three groundwater funnels.   Figure 17a shows the vertical deformation velocity map in three bowls (d0, d14, and d15) from May 2017 to May 2019. Three remarkable subsidence bowl centers (d0c, d14c, and d15c) associated with building construction were successfully detected. Two groups of Google Earth images over three subsidence bowl centers (Figure 17b) acquired in 2017 and 2019, respectively, show the detailed changes in the building construction process. Either large-scale buildings were built in these three bowls, or there were rebuilt residential areas after demolishing. Bowl d0c is located in Nanhu Wenyuan, Hanshan District of Handan City, where large new residential areas and a school were built. Subsidence bowl center d14c is located in Mingshi Mansion and Star River Garden, Gucheng County of Hengshui City. Bowl center d15c is located in Lantian Apartment and Bishuilanwan of Jingxian County. The centers of subsidence bowls d14c and d15c were rebuilt the community after demolishing low-rise houses. To study the temporal deformation evolution, three feature points over three subsidence areas were extracted and are given in Figure  17c. The maximum subsidence velocity in d0c, d14c, and d15c reached around −44 mm/year, −54 mm/year, and −57 mm/year, respectively. All these phenomena indicate that building constructions increase surface loads, aggravating land subsidence. For the first  Figure 17a shows the vertical deformation velocity map in three bowls (d0, d14, and d15) from May 2017 to May 2019. Three remarkable subsidence bowl centers (d0c, d14c, and d15c) associated with building construction were successfully detected. Two groups of Google Earth images over three subsidence bowl centers (Figure 17b) acquired in 2017 and 2019, respectively, show the detailed changes in the building construction process. Either large-scale buildings were built in these three bowls, or there were rebuilt residential areas after demolishing. Bowl d0c is located in Nanhu Wenyuan, Hanshan District of Handan City, where large new residential areas and a school were built. Subsidence bowl center d14c is located in Mingshi Mansion and Star River Garden, Gucheng County of Hengshui City. Bowl center d15c is located in Lantian Apartment and Bishuilanwan of Jingxian County. The centers of subsidence bowls d14c and d15c were rebuilt the community after demolishing low-rise houses. To study the temporal deformation evolution, three feature points over three subsidence areas were extracted and are given in Figure 17c. The maximum subsidence velocity in d0c, d14c, and d15c reached around −44 mm/year, −54 mm/year, and −57 mm/year, respectively. All these phenomena indicate that building constructions increase surface loads, aggravating land subsidence. For the first time, we tried to analyze the subsidence causes from the perspective of change detection in the Hebei Plain.

Industrial Areas and Land Subsidence
The land use may affect the groundwater balance and lead to the evolution of land subsidence [44]. As shown in Figure 18a, severe subsidence areas were mainly orange or red areas. Six remarkable subsidence bowl centers (d2c, d4c, d12c, d16c, d22c-1, and d22c-2) associated with industrial areas were found. The maximum subsidence velocity in d2c, d4c, d12c, d16c, and d22c reached around −36 mm/year, −48 mm/year, −50 mm/year, −45 mm/year, −48 mm/year, and −79 mm/year, respectively. Six feature points were extracted and their deformation time-series are given in Figure 18b. The subsidence bowl centers and the industrial areas were consistent to some extent, although the distribution range of the subsidence bowl centers were different from the industrial areas. Two groups of Google Earth images over six subsidence bowl centers (Figure 18c) acquired in 2019, respectively, show the details in the industrial area. Subsidence bowl d2c is located at the Congtai Liquor Factory, Feixiang Economic Development Zone of Handan City. It is suspected to be related to the large amount of water used for liquor making. Bowl center d4c is a gathering place of industrial areas, located in Zhaolizhuang Village of Quzhou County, Handan City. Bowl center d12c is located in Xinleitou Village of Xinji City, Shijiazhuang City, and there are 30 to 40 chemical fiber, filter paper, and food factories that require a lot of water. Figure 19 shows the accumulative time-series deformation of d12c. Bowl center d16c is located in Renhe Street and Guanzhou Road, Dongguang County of Cangzhou City. There are dozens of factory gathering areas in d16c, especially machinery factories, chemical factories, and plastic products. Bowl center d22c is located in Nansha Village and Yuejiazuo Village of Gaoyang County, Baoding City. Hundreds of companies gather in d22, especially textile companies related to oil exploration. Gaoyang County is known as the hometown of textiles, with textile enterprises accounting for 91.8%. The two horizontal and vertical deformations along the distance through the subsidence center d22c-1 are plotted in Figure 18d. In general, the larger the industrial plants, the more groundwater consumed for industrial production, leading to land subsidence. A relatively high spatial correlation also existed between locations of land subsidence bowls and industrial areas.

Industrial Areas and Land Subsidence
The land use may affect the groundwater balance and lead to the evolution of land subsidence [44]. As shown in Figure 18a, severe subsidence areas were mainly orange or red areas. Six remarkable subsidence bowl centers (d2c, d4c, d12c, d16c, d22c-1, and d22c-2) associated with industrial areas were found. The maximum subsidence velocity in d2c, d4c, d12c, d16c, and d22c reached around −36 mm/year, −48 mm/year, −50 mm/year, −45 mm/year, −48 mm/year, and −79 mm/year, respectively. Six feature points were extracted and their deformation time-series are given in Figure 18b. The subsidence bowl centers and the industrial areas were consistent to some extent, although the distribution range of the subsidence bowl centers were different from the industrial areas. Two groups of Google Earth images over six subsidence bowl centers (Figure 18c) acquired in 2019, respectively, show the details in the industrial area. Subsidence bowl d2c is located at the Congtai Liquor Factory, Feixiang Economic Development Zone of Handan City. It is suspected to be related to the large amount of water used for liquor making. Bowl center d4c is a gathering place of industrial areas, located in Zhaolizhuang Village of Quzhou County, Handan City. Bowl center d12c is located in Xinleitou Village of Xinji City, Shijiazhuang City, and there are 30 to 40 chemical fiber, filter paper, and food factories that require a lot of water. Figure 19 shows the accumulative time-series deformation of d12c. Bowl center d16c is located in Renhe Street and Guanzhou Road, Dongguang County of Cangzhou City. There are dozens of factory gathering areas in d16c, especially machinery factories, chemical factories, and plastic products. Bowl center d22c is located in Nansha Village and Yuejiazuo Village of Gaoyang County, Baoding City. Hundreds of companies gather in d22, especially textile companies related to oil exploration. Gaoyang County is known as the hometown of textiles, with textile enterprises accounting for 91.8%. The two horizontal and vertical deformations along the distance through the subsidence center d22c-1 are plotted in Figure 18d. In general, the larger the industrial plants, the more groundwater consumed for industrial production, leading to land subsidence. A relatively high spatial correlation also existed between locations of land subsidence bowls and industrial areas.

Dense Residential Areas and Land Subsidence
In addition, there a high spatial correlation exists between the locations of land subsidence bowls and dense residential areas. Figure 20a shows the deformation velocity map in the vertical direction in four subsidence bowls (d5, d7, d8, and d13) from May 2017 to May 2019. Four remarkable subsidence bowl centers (d5c, d7c, d8c, and d13c) associated with dense residential areas were successfully found. Google Earth images over four areas acquired in 2019 are shown in Figure 20b. The subsidence bowl center of d5 is located in Dawuzhuang, Xiangdu District of Xingtai City. Subsidence bowl center d7 is located in Shengfeng Times and Zhongxiadijingcheng in Wei County of Xingtai City. Subsidence bowl center d8 is located in the Shengshi Capital, Julu County of Xingtai City. Subsidence bowl center d13 is located in No. 195 Yongping Street, Hengshui City of Shenzhou City. To study the temporal deformation evolution, four feature points over four subsidence

Dense Residential Areas and Land Subsidence
In addition, there a high spatial correlation exists between the locations of land subsidence bowls and dense residential areas. Figure 20a shows the deformation velocity map in the vertical direction in four subsidence bowls (d5, d7, d8, and d13) from May 2017 to May 2019. Four remarkable subsidence bowl centers (d5c, d7c, d8c, and d13c) associated with dense residential areas were successfully found. Google Earth images over four areas acquired in 2019 are shown in Figure 20b. The subsidence bowl center of d5 is located in Dawuzhuang, Xiangdu District of Xingtai City. Subsidence bowl center d7 is located in Shengfeng Times and Zhongxiadijingcheng in Wei County of Xingtai City. Subsidence bowl center d8 is located in the Shengshi Capital, Julu County of Xingtai City. Subsidence bowl center d13 is located in No. 195 Yongping Street, Hengshui City of Shenzhou City. To study the temporal deformation evolution, four feature points over four subsidence

Dense Residential Areas and Land Subsidence
In addition, there a high spatial correlation exists between the locations of land subsidence bowls and dense residential areas. Figure 20a shows the deformation velocity map in the vertical direction in four subsidence bowls (d5, d7, d8, and d13) from May 2017 to May 2019. Four remarkable subsidence bowl centers (d5c, d7c, d8c, and d13c) associated with dense residential areas were successfully found. Google Earth images over four areas acquired in 2019 are shown in Figure 20b. The subsidence bowl center of d5 is located in Dawuzhuang, Xiangdu District of Xingtai City. Subsidence bowl center d7 is located in Shengfeng Times and Zhongxiadijingcheng in Wei County of Xingtai City. Subsidence bowl center d8 is located in the Shengshi Capital, Julu County of Xingtai City. Subsidence bowl center d13 is located in No. 195 Yongping Street, Hengshui City of Shenzhou City. To study the temporal deformation evolution, four feature points over four subsidence area were extracted and are given in Figure 20c. Figure 21 shows the accumulative time-series deformation of d12c. The maximum subsidence velocity in d5c, d7c, d8c, and d13c reached around −58 mm/year, −56 mm/year, −59 mm/year, and −51 mm/year, respectively. To sum up, land subsidence is highly consistent with dense residential areas.
area were extracted and are given in Figure 20c. Figure 21 shows the accumulative timeseries deformation of d12c. The maximum subsidence velocity in d5c, d7c, d8c, and d13c reached around −58 mm/year, −56 mm/year, −59 mm/year, and −51 mm/year, respectively. To sum up, land subsidence is highly consistent with dense residential areas.
Meanwhile, we analyzed subsidence with monthly precipitations acquired from the China Meteorological Data Service Center [56] from May 2017 to May 2019 at the p5 point. As we can see from Figure 20c, the precipitations declined noticeably from November 2017 to March 2018 and from September 2018 to March 2019 while subsidence velocity increased significantly. Precipitations were the highest from May 2017 to October 2017 and from April 2018 to August 2018 while the subsidence velocity slowed down slightly. All these phenomena suggest the following: that surface subsidence at the p5 point varies with the seasons and precipitation is one of the factors influencing land subsidence.    area were extracted and are given in Figure 20c. Figure 21 shows the accumulative timeseries deformation of d12c. The maximum subsidence velocity in d5c, d7c, d8c, and d13c reached around −58 mm/year, −56 mm/year, −59 mm/year, and −51 mm/year, respectively. To sum up, land subsidence is highly consistent with dense residential areas.

Faults and Land Subsidence
Meanwhile, we analyzed subsidence with monthly precipitations acquired from the China Meteorological Data Service Center [56] from May 2017 to May 2019 at the p5 point. As we can see from Figure 20c, the precipitations declined noticeably from November 2017 to March 2018 and from September 2018 to March 2019 while subsidence velocity increased significantly. Precipitations were the highest from May 2017 to October 2017 and from April 2018 to August 2018 while the subsidence velocity slowed down slightly. All these phenomena suggest the following: that surface subsidence at the p5 point varies with the seasons and precipitation is one of the factors influencing land subsidence.    Meanwhile, we analyzed subsidence with monthly precipitations acquired from the China Meteorological Data Service Center [56] from May 2017 to May 2019 at the p5 point. As we can see from Figure 20c, the precipitations declined noticeably from November 2017 to March 2018 and from September 2018 to March 2019 while subsidence velocity increased significantly. Precipitations were the highest from May 2017 to October 2017 and from April 2018 to August 2018 while the subsidence velocity slowed down slightly. All these phenomena suggest the following: that surface subsidence at the p5 point varies with the seasons and precipitation is one of the factors influencing land subsidence. Figure 22a demonstrates a strong correlation between the land subsidence and fault distribution in the Hebei Plain, which indicates that the main deformation direction is basically parallel to the faults. From Figure 22a, a total of 21 subsidence bowls were distributed on both sides of the major faults. Figure 22b-d shows an example of the average subsidence velocity in three fault zones along profile aa', bb', and cc'. Sharp gradients in subsidence were observed across the faults, especially profile aa' with subsidence gradients up to 40 mm/year. These high-velocity gradients correlate with faults, indicating that these faults coincide with the subsidence bowls.

Faults and Land Subsidence
Remote Sens. 2021, 13, 795 20 of 23 basically parallel to the faults. From Figure 22a, a total of 21 subsidence bowls were distributed on both sides of the major faults. Figure 22b-d shows an example of the average subsidence velocity in three fault zones along profile aa', bb', and cc'. Sharp gradients in subsidence were observed across the faults, especially profile aa' with subsidence gradients up to 40 mm/year. These high-velocity gradients correlate with faults, indicating that these faults coincide with the subsidence bowls.

Conclusions
Through using two TS-InSAR techniques, SBAS and NSBAS, we obtained a largescale deformation map (about 45,000 km 2 ) for Hebei Plain, with 166 Sentinel-1A data of tracks 40 and 142 from May 2017 to May 2019. We also identified 26 notable subsidence bowls in Handan, Xingtai, Shijiazhuang, Hengshui, Cangzhou, and Baoding, with maximum subsidence rates of −65 mm/year, −71 mm/year, −50 mm/year, −62 mm/year, −54 mm/year, and −79 mm/year, respectively. Subsidence bowls tend to be connected in the administrative area. A novel data fusion flow for the generation of land subsidence vertical velocity of adjacent-track was proposed, named FTASA. In addition, this method does

Conclusions
Through using two TS-InSAR techniques, SBAS and NSBAS, we obtained a largescale deformation map (about 45,000 km 2 ) for Hebei Plain, with 166 Sentinel-1A data of tracks 40 and 142 from May 2017 to May 2019. We also identified 26 notable subsidence bowls in Handan, Xingtai, Shijiazhuang, Hengshui, Cangzhou, and Baoding, with maximum subsidence rates of −65 mm/year, −71 mm/year, −50 mm/year, −62 mm/year, −54 mm/year, and −79 mm/year, respectively. Subsidence bowls tend to be connected in the administrative area. A novel data fusion flow for the generation of land subsidence vertical velocity of adjacent-track was proposed, named FTASA. In addition, this method does not need to interpolate and correct the PS point, which further reduces the error and improves the accuracy. FTASA was successfully applied to the deformation monitoring in the Hebei Plain. FTASA had a higher accuracy than FO. The InSAR results with leveling measurements show good agreement with a precision of 9 mm/year. The InSAR results within tracks 40 and 142 and benchmark on the bedrock results also had a good correlation. The standard deviation of the InSAR results velocity calculated between tracks 40 and 142 was 7.37 mm/year. The correlation coefficient of track 40 and 142 was 0.89.
We found that eight subsidence bowls were positively correlated with Ningbolong and Gaoliqing-Sunin shallow groundwater funnels and Nangong deep groundwater funnel. Furthermore, we found that three, five, four, and three subsidence bowls were positively correlated with building constructions, industrial areas, dense residential areas, and faults, respectively. In our future work, we will try to combine change detection by using deep learning to analyze the influence of changes in different land types on the deformation results. A more detailed summary and future work are as follows:

•
The land subsidence bowls have a positive correlation with the shallow or deep groundwater funnels. We found that there were six ground subsidence bowls in the shallow ground-water funnel, and two ground subsidence bowls in the deep groundwater funnel. In the future, the mechanisms underlying groundwater-induced land surface displacement deserves further analysis. • Building constructions increase surface loads, aggravating land subsidence. Three remarkable subsidence bowl centers associated with building construction were successfully detected. For the first time, we tried to analyze the subsidence causes from the perspective of change detection in the Hebei Plain. Whether it is the load of the building that provokes the subsidence, or the required groundwater pumping to have a dry environment to build the foundations of those buildings, the subsidence mechanism needs to be further studied.

•
Six remarkable subsidence bowl centers associated with an industrial area were found. The larger the industrial plants, the more groundwater consumed for industrial production, leading to land subsidence. A relatively high spatial correlation also existed between locations of land subsidence bowls and industrial areas. • A high spatial correlation existed between the locations of four land subsidence bowls and dense residential areas. Precipitation was one of the factors influencing land subsidence. We can analyze the causes of land subsidence from the hydrogeological properties of rocks in our future work.

•
There was a strong correlation between three land subsidence and fault distributions where the main deformation direction is basically parallel to the faults. A total of 21 subsidence bowls were distributed on both sides of major faults. These highvelocity gradients correlate with faults in three subsidence bowls, indicating that these faults coincide with the subsidence bowls.