# Detecting, Monitoring, and Analyzing the Surface Subsidence in the Yellow River Delta (China) Combined with CenterNet Network and SBAS-InSAR

Authors: Zhenjin Li; Zhiyong Wang; Wei Liu; Xing Li; Maotong Zhou; Baojing Zhang
Journal: Journal of Spectroscopy
 (2022)

---

## Abstract

Long-term industrial activities tend to cause surface subsidence and damage to ground facilities and local ecological environment. Monitoring and analyzing surface subsidence is of great significance to prevent potential disasters. The surface type of the Yellow River Delta in China is complex and there are many industrial activities, so it is necessary to monitor the surface subsidence in this area. Small Baseline Subset InSAR (SBAS-InSAR) can monitor the surface subsidence with millimeter-level accuracy, but it takes a long time to process wide images (Sentinel-1) and is seriously affected by atmospheric errors. To avoid these limitations, we constructed a method combining the CenterNet network and SBAS-InSAR (CNSBAS-InSAR). Firstly, the CenterNet network is used to automatically detect the subsidence areas from the wide differential interferogram formed by two SAR satellite images and determine the location of the subsidence area. Then, the SBAS-InSAR monitoring is performed on the detected multiple subsidence areas. Finally, the small-scale subsidence results are obtained. In this study, based on 24 Sentinel-1A satellite images acquired from 10 January 2018 to 24 December 2018, nine subsidence areas in Yellow River Delta were detected. Three of them had long-term surface subsidence. They were located in Zhanhua District, Xianhe Town, and Hongguang Village, respectively. This paper focuses on analyzing these three areas. The maximum subsidence rate of Zhanhua District, Xianhe Town, and Hongguang Village were −135.21 mm/a, −330.91 mm/a, and −209.68 mm/a, respectively. In addition, the analysis showed that precipitation in the Zhanhua District could effectively slow down the subsidence rate of the area. The subsidence of Xianhe Town threatened the safety of the Shugang Expressway. The subsidence of Hongguang Village caused the safety risks of buildings. The results of this study prove that CNSBAS-InSAR method is reliable for monitoring subsidence areas and it can provide a reference for local construction and protection of Yellow River Delta.


---

## Body

## 1. Introduction

Surface subsidence is a complex geological disaster affected by natural and artificial factors [1]. The main causes of surface subsidence include mining, urban groundwater extraction, and industrial activities [2–4]. Large-scale subsidence can easily damage ground infrastructure, causing foundation collapse and road cracks. In addition, rapid subsidence in a short time may lead to surface collapse, ecological damage, landslides, and debris flow [5–7]. With the change in the ecological environment and the continuous exploitation of resources, it is becoming more important to monitor the surface subsidence. It can give early warning of danger and disaster in time to ensure human safety and reduce economic loss.In recent years, the interferometric synthetic aperture radar (InSAR) technique has been widely applied to monitor ground subsidence in various situations, such as airports and oil fields [8, 9]. Compared with the traditional global navigation satellite system (GNSS) and precise leveling [10, 11], the monitoring accuracy of InSAR can reach centimeter-level, and it has the characteristics of saving time and labour, low cost, and wide monitoring range [12]. However, in the long-term monitoring, InSAR is limited by signal incoherence, which is caused by atmospheric phase error and noise [13], resulting in reduced monitoring accuracy and a loss of ability to monitor tiny deformations. In view of this problem, several scholars improved the phase filtering [14] and phase unwrapping [15] methods in InSAR or established a refined model for a single subsidence basin [16] to improve the monitoring accuracy. In 2002, Small Baseline Subset InSAR (SBAS-InSAR) was proposed by Berardino et al. [17]. This technique combined many differential interferograms produced by data pairs to limit the spatial decorrelation phenomena and enabled the monitoring accuracy to reach a millimeter level.The Yellow River Delta has developed industry, agriculture, aquaculture, etc. Shengli Oilfield, the second largest oil field in China, is located in the area. A large number of mining work, special geographical locations, and complex industrial structures have caused serious damage to the surface in many areas of the Yellow River Delta [18]. Therefore, it is necessary to monitor the surface subsidence in this area. At present, the InSAR technique has been widely used in surface subsidence monitoring in this area. Liu et al. [19] combined the SBAS-InSAR with a geotechnical model and geological, geomorphological, hydrogeological, and geotechnical data to assess the long-term spatiotemporal subsidence in the Yellow River Delta. Zhang et al. [20] developed a multitemporal InSAR method to map ground subsidence over the Yellow River Delta area by exploiting both persistent scatterers (PS) and distributed scatterers (DS) and obtained surface subsidence from 2007 to 2010 and 2015 to 2018. The above methods improve the accuracy of monitoring results by introducing parameters and combining different algorithms, but they require the operation of a wide image, including many areas without subsidence.Due to the large amount of data required by SBAS-InSAR, it takes a long time to process the original wide image (Sentinel-1), and it is seriously affected by atmospheric error during large-scale monitoring [21], resulting in lower accuracy of monitoring results. At present, several scholars have proposed new methods to reduce the errors in the SBAS-InSAR monitoring process. For example, Duan et al. [13] took incoherent noise and atmospheric turbulence as random variables and proposed a new weighting method to reduce the impact on the monitoring results. Havazli and Wdowinski [22] proposed a method to estimate the detection threshold of SBAS-InSAR to reduce tropospheric delay to improve monitoring accuracy. Most of these methods are improvements to the algorithm, less for the image. In most cases, there are few subsidence areas in a wide image like Sentinel-1, and we pay more attention to the subsidence area, while monitoring the area without subsidence has not much significance and will consume a lot of time. Therefore, these problems can be effectively solved by finding suitable SBAS-InSAR monitoring locations. Wang et al. [23] used the histogram of oriented gradients and a support vector machine model in machine learning to extract the subsidence basins of the mining area from the wide differential interferogram, and the accuracy rate reached 85%. Bata et al. [24] proposed the automatic subsidence area detection in SAR interferograms by the method of circlet transform and tested it on the Upper Silesian Coal Basin located in Southern Poland, and the detection efficiency of this method was improved by 20% compared with the Hough transform. However, traditional mathematical and machine learning methods are easily limited by the amount of data and need to add features manually [25]. At present, most of the methods can only perform manual clipping of the original wide image after determining the geographic position of the subsidence area through field investigation [26], which takes a long time and is prone to omission.Object detection technology [27] under the framework of deep learning has made practical progress in many fields, such as detection of coastline garbage [28], animal species [29], and vehicles [30], which greatly facilitates people’s work. In the application of SAR remote sensing, object detection technology is mainly used to identify ships [31] and offshore oil spills [32]. By combining the object detection technology, the subsidence areas could be automatically identified in the wide differential interferogram, and the small-scale subsidence areas could be located. Then the SBAS-InSAR monitoring is performed for the small subsidence areas, which saves a lot of monitoring time and reduces the influence of atmospheric phase error on monitoring accuracy. The CenterNet network is an object detection technology proposed by Zhou et al. [33] in 2019. It had no anchor boxes and modeled the object as a single point detection mode so that the detection results were not limited by the object size, which was very suitable for detecting subsidence areas. In this study, we constructed a method combining CenterNet and SBAS-InSAR technology (CNSBAS-InSAR) for monitoring surface subsidence. Applying this method to the Yellow River Delta, the location of subsidence was determined and the millimeter-scale subsidence results were obtained successfully.The sections of this paper are arranged as follows: In Section2, detailed information about the study area and data is introduced; Section 3 is the method of this paper. It mainly introduces the method of detecting subsidence area automatically by CenterNet network and obtaining subsidence rate and time series deformation by SBAS-InSAR technology; Section 4 is experimental results and validation, including subsidence areas detection results, subsidence rate results, and validation by different SAR tracks; Section 5 is about the discussion and analysis of subsidence area in the experimental area. In Section 6, some important conclusions are presented.


## 2. Materials and Data

### 2.1. Study Area

The Yellow River Delta is located in Shandong Province, China, between 36°55′–38°16′N, 117°31′–119°18′E, bordering Laizhou Bay in the east and Bohai Bay in the north and occupying an area of over 10,000 km2. The terrain is flat and the altitude is between 0 m and 15 m. It is located in the warm temperate zone and has a temperate continental monsoon climate. The annual average precipitation is about 551.6 mm. The Yellow River Delta is rich in biodiversity and natural resources and has been extensively studied by scholars [34, 35].The study area in this paper is the coverage of the Sentinel-1A image, as shown in the red box in Figure1(a), and the coverage of Sentinel-1B images used for validation is shown in the blue box. The background is DEM, and the altitude is denoted by the color scale from white to jade. Figure 1(b) shows the location of the Yellow River Delta in Shandong Province, China, Asia. Figure 1(c) shows the Landsat 8 image of the part of the study area for analysis in this paper, and the imaging time is 25 March 2018.Figure 1
Location of the study area. (a) SRTM DEM 90 m image of the study area, the red box denotes the Sentinel-1A range, and the blue box denotes the Sentinel-1B range; (b) the study area in Shandong Province, China, Asia; (c) the Landsat 8 image of part of the study area.
(a)(b)(c)


### 2.2. Data

To monitor the surface subsidence in the Yellow River Delta in 2018, 24 ascending Sentinel-1A images from the ESA were acquired from 10 January 2018 to 24 December 2018. The 24 descending Sentinel-1B images were used to verify the accuracy of subsidence results. The parameters of the Sentinel-1 data are shown (see Table1).Table 1
Parameters of Sentinel-1.
ParametersSentinel-1ASentinel-1BBandCCTrack6976Beam modeIW∗IWIncidence angle (°)38.939.2PolarizationVVVVOribit directionAscendingDescendingResolution (m)2.3 × 13.92.3 × 13.9∗IW: Interferometric wide swath.In addition, precise orbit determination (POD) data released by ESA was used for the orbital refinement and phase reflattening. Topographic phases were removed using shuttle radar topography mission (SRTM) DEM data with 90 m resolution provided by the National Aeronautics and Space Administration (NASA).


## 2.1. Study Area

The Yellow River Delta is located in Shandong Province, China, between 36°55′–38°16′N, 117°31′–119°18′E, bordering Laizhou Bay in the east and Bohai Bay in the north and occupying an area of over 10,000 km2. The terrain is flat and the altitude is between 0 m and 15 m. It is located in the warm temperate zone and has a temperate continental monsoon climate. The annual average precipitation is about 551.6 mm. The Yellow River Delta is rich in biodiversity and natural resources and has been extensively studied by scholars [34, 35].The study area in this paper is the coverage of the Sentinel-1A image, as shown in the red box in Figure1(a), and the coverage of Sentinel-1B images used for validation is shown in the blue box. The background is DEM, and the altitude is denoted by the color scale from white to jade. Figure 1(b) shows the location of the Yellow River Delta in Shandong Province, China, Asia. Figure 1(c) shows the Landsat 8 image of the part of the study area for analysis in this paper, and the imaging time is 25 March 2018.Figure 1
Location of the study area. (a) SRTM DEM 90 m image of the study area, the red box denotes the Sentinel-1A range, and the blue box denotes the Sentinel-1B range; (b) the study area in Shandong Province, China, Asia; (c) the Landsat 8 image of part of the study area.
(a)(b)(c)


## 2.2. Data

To monitor the surface subsidence in the Yellow River Delta in 2018, 24 ascending Sentinel-1A images from the ESA were acquired from 10 January 2018 to 24 December 2018. The 24 descending Sentinel-1B images were used to verify the accuracy of subsidence results. The parameters of the Sentinel-1 data are shown (see Table1).Table 1
Parameters of Sentinel-1.
ParametersSentinel-1ASentinel-1BBandCCTrack6976Beam modeIW∗IWIncidence angle (°)38.939.2PolarizationVVVVOribit directionAscendingDescendingResolution (m)2.3 × 13.92.3 × 13.9∗IW: Interferometric wide swath.In addition, precise orbit determination (POD) data released by ESA was used for the orbital refinement and phase reflattening. Topographic phases were removed using shuttle radar topography mission (SRTM) DEM data with 90 m resolution provided by the National Aeronautics and Space Administration (NASA).


## 3. Methods

To solve the problem of rapid and accurate monitoring of surface subsidence in the Yellow River Delta, we constructed a novel method combining the CenterNet network and SBAS-InSAR technology (called CNSBAS-InSAR) for monitoring surface subsidence. Firstly, the sample datasets of subsidence areas were made based on the characteristics of the subsidence areas. Then, the CenterNet network was used to detect the subsidence areas from the wide interferogram. Finally, SBAS-InSAR processing was conducted on the multiple selected subsidence areas of the Yellow River Delta, and the small-scale accurate subsidence results were obtained. The flow chart of CNSBAS-InSAR is shown in Figure2, which includes the generation of a wide interferogram, detecting the subsidence areas with the CenterNet network, locating the subsidence areas, cropping the image, and obtaining the small-scale subsidence results.Figure 2
The flow chart of CNSBAS-InSAR.


### 3.1. CenterNet Network

#### 3.1.1. Technical Principle

The CenterNet network is one of the object detection methods of deep learning, which is widely used in industry [36], transportation [37], and other fields. The principle [33] is modeling the object to be detected as a single point, namely the center point of the bounding box, and then determining the center point through the keypoint heatmap. Other features of the object can be regressed according to the image features of the center point, such as object size, 3D extent, orientation, and pose. These features are input into the network to directly predict the height and weight of the object box. The CenterNet supports Hourglass [38], DLA [39], and ResNet [40] as the backbone. Due to a large number of Hourglass and DLA network parameters, it is not convenient for practical application. Therefore, ResNet50 was adopted as the backbone in this paper.The Yellow River Delta is a complex industrial area, and the size of the subsidence areas caused by different industrial activities varies greatly. CenterNet network is not limited by the size of the object to be detected because it has no anchor box and adopts the method of heatmap regression prediction. Therefore, it is very suitable for detecting the subsidence areas in wide differential interferograms. In this paper, CenterNet was used to automatically detect the subsidence areas from wide interferograms before SBAS-InSAR processing.


#### 3.1.2. Datasets Making

The sample datasets for training the CenterNet network were obtained from differential interferograms of real subsidence areas. The two-pass approach [41] was used to carry out DInSAR processing on some areas with many subsidence areas. The DInSAR process only needs to obtain the differential interferogram after phase noise filtering, so as to accurately obtain the characteristics of the subsidence area on the interferogram. Then, the interferograms were cropped into some subimages with a size of 416 × 416 and the image data annotation software “LabelImg” was used to draw the outer rectangular boxes of subsidence areas on each sample image. The subsidence areas which lost features such as interferometric fringes due to atmospheric phase errors and decorrelation noise were removed. Data augmentation was used in this paper with rotation, translation, and flipping to expand SAR data samples [42–44]. In this study, by operating on 12 differential interferograms in the Huaibei mining area, a total of 611 subsidence area sample datasets were produced. We used 80% of the samples for training and 20% for testing. The parameters of interferograms are listed in Table 2. The sample datasets of some subsidence areas are shown in Figure 3.Table 2
Parameters of interferograms for making sample datasets.
Master imageSlave imagePathSpatial baseline (m)Temporal baseline (d)22/12/2017∗27/01/2018142−80.563615/01/201808/02/2018142−101.252427/01/201804/03/2018142−32.343608/02/201804/03/201814228.722411/11/201817/12/201814238.833617/12/201822/01/2019142−50.793617/12/201827/02/2019142−84.143603/02/201927/02/2020142−105.332406/11/201930/11/201914239.522430/11/201924/12/201914272.782405/01/202024/12/2019142−41.112429/01/202022/02/202014242.4824∗Note. Date/month/year in Table 1.Figure 3
Partial sample datasets for training the CenterNet network.


#### 3.1.3. Training and Evaluation of Detection Capability

The model training was carried out using the prepared sample datasets. The platform was Windows 10. The processor was an Intel Core I5-8400H (8 GB memory). The deep learning framework was TensorFlow 2.2.0, with CPU for training and CUDA for acceleration.The loss function of CenterNet is divided into three parts: heatmap loss, offset loss, and size loss. In CenterNet, each center point corresponds to an object location, and overlapping judgment is not required. In practical training, the loss of other points around the center point is the loss after attenuation, while the length and width of the object are obtained by regression and correspond to the current center point. We have recorded the loss curve during the training of CenterNet, as shown in Figure4. The loss curve decreased significantly during the training process and then tended to be stable after epoch = 30 when the loss value remained at about 1.8. The results showed that the CenterNet network could be a good fit to detect subsidence areas.Figure 4
Loss curve trained by CenterNet.For the detection results,T and F, respectively, indicated that the sample was correctly classified and incorrectly classified. P and N indicated that the sample was detected to be a positive sample and a negative sample. There were four types of detection results of the model: true positive (TP) was the positive sample that was correctly detected; false positive (FP) was the positive sample that was incorrectly detected; true negative (TN) was the negative sample that was correctly detected; false negative (FN) was the negative sample that was incorrectly detected. Precision (P) indicated the proportion of all detected samples that were correct. Recall (R) indicated the proportion of the objects recognized by the network among all the objects that are required to be recognized. The equations of precision (P) and recall (R) are as follows:(1)P=TPTP+FP,(2)R=TPTP+FN.TheP-R curve took the precision P as the ordinate and the recall as the abscissa. The equation for AP is shown in the following formula:(3)AP=∑PM​N,where PM is the detection accuracy of subsidence areas in each image and N is the number of images.In this paper, AP measured the detection performance of subsidence areas; that is, the precision values on theP-R curve were averaged. The calculated AP value was 93.69%, indicating that CenterNet can better detect the subsidence areas. The P-R curve is shown in Figure 5.Figure 5
P-R curve of the CenterNet network. Note that the AP value is marked at the top of the figure.


### 3.2. CNSBAS-InSAR

Before SBAS-InSAR processing, the two Sentinel-1A images were registered first, and then the DInSAR was processed [41] to obtain the wide differential interferogram. Goldstein filtering [45] was used to reduce phase noise on the image to make the interferometric fringe clear. In the process, the interferogram was cut into 36 subimages according to 6 × 6 clipping. Grey scales were added to the blank part of the image to make the pixel values of width and height multiples of 32, which met the detection requirements of CenterNet. Then, the 36 subimages were detected by the CenterNet network in turn to obtain the detection results of the subsidence areas, and several subsidence subregions were located according to the detection results.According to the pixel position of the subsidence subregions in the original image, theN+1 Sentinel-1A images requiring SBAS-InSAR processing were cropped and constructed into a time series [3, 46]. The acquisition time of SAR images covering the subregions with N+1 scenes was T1,T2,…,TN and assuming each scene image formed an interferogram with another scene image. Then, M interferograms were generated based on spatial and temporal baseline threshold value. M satisfied the following equation:(4)N+12≤M≤NN+12.Assuming interferogramj was generated by combining SAR acquisitions at times tA and tBtB>tA, and then the differential interferometric phase of j in pixel of azimuth and range coordinates x,r was mainly composed of five parts, which is expressed in the following equation:(5)φjx,r=φjtB,x,r−φjtA,x,r≈φdisp+φtopo+φorb+φatm+φnoise,where x and r are azimuth and range coordinates, respectively. φdisp is the phase caused by changes in the line of sight (LOS) distance between the object and the radar; φtopo is the residual phase due to inaccuracies in reference DEM; φorb is the satellite orbit phase error; φatm is the atmospheric phase error; and φnoise is the phase caused by other noise.In order to obtain accurate deformation results,φtopo, φorb, φatm, and φnoise should be removed, then a system of M equations in N unknowns is obtained from equation (5). The matrix form w expressed in the following equation:(6)Aφ=δφ,where A is the M×N coefficient matrix; M and N denotes the number of interferograms and SAR acquisitions, respectively; φ is N unknown phase values related to high-coherence pixels; and δ is the vector of unwrapped phase values associated with differential interferograms.In order to obtain the deformation rate of highly coherent pixels, the following equation was formed:(7)Bv=δφ,where B is the M×N coefficient matrix and vT can be expressed as equation:(8)v=v1=φ1t1−t0,…,vN=φN−φN−1tN−tN−1.Finally, the deformation rate was obtained by singular value decomposition (SVD) in equation (7) and the time series deformation was obtained according to the time span between SAR image acquisitions.


## 3.1. CenterNet Network

### 3.1.1. Technical Principle

The CenterNet network is one of the object detection methods of deep learning, which is widely used in industry [36], transportation [37], and other fields. The principle [33] is modeling the object to be detected as a single point, namely the center point of the bounding box, and then determining the center point through the keypoint heatmap. Other features of the object can be regressed according to the image features of the center point, such as object size, 3D extent, orientation, and pose. These features are input into the network to directly predict the height and weight of the object box. The CenterNet supports Hourglass [38], DLA [39], and ResNet [40] as the backbone. Due to a large number of Hourglass and DLA network parameters, it is not convenient for practical application. Therefore, ResNet50 was adopted as the backbone in this paper.The Yellow River Delta is a complex industrial area, and the size of the subsidence areas caused by different industrial activities varies greatly. CenterNet network is not limited by the size of the object to be detected because it has no anchor box and adopts the method of heatmap regression prediction. Therefore, it is very suitable for detecting the subsidence areas in wide differential interferograms. In this paper, CenterNet was used to automatically detect the subsidence areas from wide interferograms before SBAS-InSAR processing.


### 3.1.2. Datasets Making

The sample datasets for training the CenterNet network were obtained from differential interferograms of real subsidence areas. The two-pass approach [41] was used to carry out DInSAR processing on some areas with many subsidence areas. The DInSAR process only needs to obtain the differential interferogram after phase noise filtering, so as to accurately obtain the characteristics of the subsidence area on the interferogram. Then, the interferograms were cropped into some subimages with a size of 416 × 416 and the image data annotation software “LabelImg” was used to draw the outer rectangular boxes of subsidence areas on each sample image. The subsidence areas which lost features such as interferometric fringes due to atmospheric phase errors and decorrelation noise were removed. Data augmentation was used in this paper with rotation, translation, and flipping to expand SAR data samples [42–44]. In this study, by operating on 12 differential interferograms in the Huaibei mining area, a total of 611 subsidence area sample datasets were produced. We used 80% of the samples for training and 20% for testing. The parameters of interferograms are listed in Table 2. The sample datasets of some subsidence areas are shown in Figure 3.Table 2
Parameters of interferograms for making sample datasets.
Master imageSlave imagePathSpatial baseline (m)Temporal baseline (d)22/12/2017∗27/01/2018142−80.563615/01/201808/02/2018142−101.252427/01/201804/03/2018142−32.343608/02/201804/03/201814228.722411/11/201817/12/201814238.833617/12/201822/01/2019142−50.793617/12/201827/02/2019142−84.143603/02/201927/02/2020142−105.332406/11/201930/11/201914239.522430/11/201924/12/201914272.782405/01/202024/12/2019142−41.112429/01/202022/02/202014242.4824∗Note. Date/month/year in Table 1.Figure 3
Partial sample datasets for training the CenterNet network.


### 3.1.3. Training and Evaluation of Detection Capability

The model training was carried out using the prepared sample datasets. The platform was Windows 10. The processor was an Intel Core I5-8400H (8 GB memory). The deep learning framework was TensorFlow 2.2.0, with CPU for training and CUDA for acceleration.The loss function of CenterNet is divided into three parts: heatmap loss, offset loss, and size loss. In CenterNet, each center point corresponds to an object location, and overlapping judgment is not required. In practical training, the loss of other points around the center point is the loss after attenuation, while the length and width of the object are obtained by regression and correspond to the current center point. We have recorded the loss curve during the training of CenterNet, as shown in Figure4. The loss curve decreased significantly during the training process and then tended to be stable after epoch = 30 when the loss value remained at about 1.8. The results showed that the CenterNet network could be a good fit to detect subsidence areas.Figure 4
Loss curve trained by CenterNet.For the detection results,T and F, respectively, indicated that the sample was correctly classified and incorrectly classified. P and N indicated that the sample was detected to be a positive sample and a negative sample. There were four types of detection results of the model: true positive (TP) was the positive sample that was correctly detected; false positive (FP) was the positive sample that was incorrectly detected; true negative (TN) was the negative sample that was correctly detected; false negative (FN) was the negative sample that was incorrectly detected. Precision (P) indicated the proportion of all detected samples that were correct. Recall (R) indicated the proportion of the objects recognized by the network among all the objects that are required to be recognized. The equations of precision (P) and recall (R) are as follows:(1)P=TPTP+FP,(2)R=TPTP+FN.TheP-R curve took the precision P as the ordinate and the recall as the abscissa. The equation for AP is shown in the following formula:(3)AP=∑PM​N,where PM is the detection accuracy of subsidence areas in each image and N is the number of images.In this paper, AP measured the detection performance of subsidence areas; that is, the precision values on theP-R curve were averaged. The calculated AP value was 93.69%, indicating that CenterNet can better detect the subsidence areas. The P-R curve is shown in Figure 5.Figure 5
P-R curve of the CenterNet network. Note that the AP value is marked at the top of the figure.


## 3.1.1. Technical Principle

The CenterNet network is one of the object detection methods of deep learning, which is widely used in industry [36], transportation [37], and other fields. The principle [33] is modeling the object to be detected as a single point, namely the center point of the bounding box, and then determining the center point through the keypoint heatmap. Other features of the object can be regressed according to the image features of the center point, such as object size, 3D extent, orientation, and pose. These features are input into the network to directly predict the height and weight of the object box. The CenterNet supports Hourglass [38], DLA [39], and ResNet [40] as the backbone. Due to a large number of Hourglass and DLA network parameters, it is not convenient for practical application. Therefore, ResNet50 was adopted as the backbone in this paper.The Yellow River Delta is a complex industrial area, and the size of the subsidence areas caused by different industrial activities varies greatly. CenterNet network is not limited by the size of the object to be detected because it has no anchor box and adopts the method of heatmap regression prediction. Therefore, it is very suitable for detecting the subsidence areas in wide differential interferograms. In this paper, CenterNet was used to automatically detect the subsidence areas from wide interferograms before SBAS-InSAR processing.


## 3.1.2. Datasets Making

The sample datasets for training the CenterNet network were obtained from differential interferograms of real subsidence areas. The two-pass approach [41] was used to carry out DInSAR processing on some areas with many subsidence areas. The DInSAR process only needs to obtain the differential interferogram after phase noise filtering, so as to accurately obtain the characteristics of the subsidence area on the interferogram. Then, the interferograms were cropped into some subimages with a size of 416 × 416 and the image data annotation software “LabelImg” was used to draw the outer rectangular boxes of subsidence areas on each sample image. The subsidence areas which lost features such as interferometric fringes due to atmospheric phase errors and decorrelation noise were removed. Data augmentation was used in this paper with rotation, translation, and flipping to expand SAR data samples [42–44]. In this study, by operating on 12 differential interferograms in the Huaibei mining area, a total of 611 subsidence area sample datasets were produced. We used 80% of the samples for training and 20% for testing. The parameters of interferograms are listed in Table 2. The sample datasets of some subsidence areas are shown in Figure 3.Table 2
Parameters of interferograms for making sample datasets.
Master imageSlave imagePathSpatial baseline (m)Temporal baseline (d)22/12/2017∗27/01/2018142−80.563615/01/201808/02/2018142−101.252427/01/201804/03/2018142−32.343608/02/201804/03/201814228.722411/11/201817/12/201814238.833617/12/201822/01/2019142−50.793617/12/201827/02/2019142−84.143603/02/201927/02/2020142−105.332406/11/201930/11/201914239.522430/11/201924/12/201914272.782405/01/202024/12/2019142−41.112429/01/202022/02/202014242.4824∗Note. Date/month/year in Table 1.Figure 3
Partial sample datasets for training the CenterNet network.


## 3.1.3. Training and Evaluation of Detection Capability

The model training was carried out using the prepared sample datasets. The platform was Windows 10. The processor was an Intel Core I5-8400H (8 GB memory). The deep learning framework was TensorFlow 2.2.0, with CPU for training and CUDA for acceleration.The loss function of CenterNet is divided into three parts: heatmap loss, offset loss, and size loss. In CenterNet, each center point corresponds to an object location, and overlapping judgment is not required. In practical training, the loss of other points around the center point is the loss after attenuation, while the length and width of the object are obtained by regression and correspond to the current center point. We have recorded the loss curve during the training of CenterNet, as shown in Figure4. The loss curve decreased significantly during the training process and then tended to be stable after epoch = 30 when the loss value remained at about 1.8. The results showed that the CenterNet network could be a good fit to detect subsidence areas.Figure 4
Loss curve trained by CenterNet.For the detection results,T and F, respectively, indicated that the sample was correctly classified and incorrectly classified. P and N indicated that the sample was detected to be a positive sample and a negative sample. There were four types of detection results of the model: true positive (TP) was the positive sample that was correctly detected; false positive (FP) was the positive sample that was incorrectly detected; true negative (TN) was the negative sample that was correctly detected; false negative (FN) was the negative sample that was incorrectly detected. Precision (P) indicated the proportion of all detected samples that were correct. Recall (R) indicated the proportion of the objects recognized by the network among all the objects that are required to be recognized. The equations of precision (P) and recall (R) are as follows:(1)P=TPTP+FP,(2)R=TPTP+FN.TheP-R curve took the precision P as the ordinate and the recall as the abscissa. The equation for AP is shown in the following formula:(3)AP=∑PM​N,where PM is the detection accuracy of subsidence areas in each image and N is the number of images.In this paper, AP measured the detection performance of subsidence areas; that is, the precision values on theP-R curve were averaged. The calculated AP value was 93.69%, indicating that CenterNet can better detect the subsidence areas. The P-R curve is shown in Figure 5.Figure 5
P-R curve of the CenterNet network. Note that the AP value is marked at the top of the figure.


## 3.2. CNSBAS-InSAR

Before SBAS-InSAR processing, the two Sentinel-1A images were registered first, and then the DInSAR was processed [41] to obtain the wide differential interferogram. Goldstein filtering [45] was used to reduce phase noise on the image to make the interferometric fringe clear. In the process, the interferogram was cut into 36 subimages according to 6 × 6 clipping. Grey scales were added to the blank part of the image to make the pixel values of width and height multiples of 32, which met the detection requirements of CenterNet. Then, the 36 subimages were detected by the CenterNet network in turn to obtain the detection results of the subsidence areas, and several subsidence subregions were located according to the detection results.According to the pixel position of the subsidence subregions in the original image, theN+1 Sentinel-1A images requiring SBAS-InSAR processing were cropped and constructed into a time series [3, 46]. The acquisition time of SAR images covering the subregions with N+1 scenes was T1,T2,…,TN and assuming each scene image formed an interferogram with another scene image. Then, M interferograms were generated based on spatial and temporal baseline threshold value. M satisfied the following equation:(4)N+12≤M≤NN+12.Assuming interferogramj was generated by combining SAR acquisitions at times tA and tBtB>tA, and then the differential interferometric phase of j in pixel of azimuth and range coordinates x,r was mainly composed of five parts, which is expressed in the following equation:(5)φjx,r=φjtB,x,r−φjtA,x,r≈φdisp+φtopo+φorb+φatm+φnoise,where x and r are azimuth and range coordinates, respectively. φdisp is the phase caused by changes in the line of sight (LOS) distance between the object and the radar; φtopo is the residual phase due to inaccuracies in reference DEM; φorb is the satellite orbit phase error; φatm is the atmospheric phase error; and φnoise is the phase caused by other noise.In order to obtain accurate deformation results,φtopo, φorb, φatm, and φnoise should be removed, then a system of M equations in N unknowns is obtained from equation (5). The matrix form w expressed in the following equation:(6)Aφ=δφ,where A is the M×N coefficient matrix; M and N denotes the number of interferograms and SAR acquisitions, respectively; φ is N unknown phase values related to high-coherence pixels; and δ is the vector of unwrapped phase values associated with differential interferograms.In order to obtain the deformation rate of highly coherent pixels, the following equation was formed:(7)Bv=δφ,where B is the M×N coefficient matrix and vT can be expressed as equation:(8)v=v1=φ1t1−t0,…,vN=φN−φN−1tN−tN−1.Finally, the deformation rate was obtained by singular value decomposition (SVD) in equation (7) and the time series deformation was obtained according to the time span between SAR image acquisitions.


## 4. Results

### 4.1. Detection Results

In this paper, two interferometric pairs of the Yellow River Delta from 10 January 2018 to 23 March 2018 and from 06 November 2018 to 24 December 2018 were selected. Due to the serious decorrelation noise generated by vegetation at other times, the subsidence areas lost their characteristics, so the detection effect was not good. Therefore, we no longer showed and listed their results. The parameters of interferograms are shown in Table3.Table 3
Parameters of interferograms for detecting the subsidence areas.
Master imageSlave imagePathSpatial baseline (m)Temporal baseline (d)22/12/2017∗27/01/2018142−80.563615/01/201808/02/2018142−101.2524Subsidence areas were detected in 9 subareas. The interferogram from 10 January 2018 to 23 March 2018 detected 7 subsidence areas as shown in Figure6(a). They were located in Shuinan Village, Zhangqiu District, Jinan ①, Maqiao Town, Huantai County, Zibo ②, Hanting District, Weifang ③, Yangkou Town, Shouguang City, Weifang ④, Hongguang Village, Kenli District, Dongying ⑤, Zhanhua District, Binzhou ⑥, and Xianhe Town, Hekou District, Dongying ⑦. The interferogram from 06 November 2018 to 24 December 2018 detected 6 subsidence areas as shown in Figure 6(b). They were located in Yangkou Town, Shouguang City, Weifang ④), Hongguang Village, Kenli District, Dongying, ⑤, Zhanhua District, Binzhou ⑥), Xianhe Town, Hekou District, Dongying ⑦, south of the Yellow River estuary in Kenli District, Dongying ⑧, and Mazishan Town, Wudi County, Binzhou ⑨. Among them, Hongguang Village ⑤, Zhanhua District ⑥, and Xianhe Town ⑦ had the same location of the subsidence areas in the interferogram of both time periods, indicated the surface subsidence lasted a long time. The detection results of each subregions are shown in Figures 6(a) and 6(b). The lines in Figures 6(a) and 6(b) denote image cropping mode, and the subsidence areas are marked by red boxes. Hongguang Village, Zhanhua District, and Xianhe Town are shown in the enlarged Figures 6(c), 6(d), and 6(e).Figure 6
Detection results of subsidence areas. (a) Detection results from 10/01/2018 to 23/03/2018 interferometric pair; (b) detection results from 13/10/2018 to 24/12/2018 interferometric pair. The results of (c) Hongguang village ⑤, (d) Zhanhua district ⑥, and (e) Xianhe town ⑦ are enlarged, respectively.
(a)(b)(c)(d)(e)


### 4.2. SBAS-InSAR Results

#### 4.2.1. Subsidence Rate

The detected subsidence areas were cropped from Sentinel-1A images from January 2018 to December 2018 to form a time series. A maximum temporal baseline was less than 180 days and spatial baseline was shorter than 200 m generated. External registration was carried out with SRTM DEM 90M provided by NASA. Goldstein filter [45] was used to eliminate the noise phase. The minimum cost flow method (MCF) [47] was used for phase unwrapping. The Yellow River Delta was covered with vegetation in many places, so there was a wide range of low coherence areas after SAR interferometry. The unwrapping decomposition level was set to 2 to improve the quality and efficiency of unwrapping and reduce unwrapping errors. The interferometric pairs with low coherence and poor unwrapping were removed. A smooth ground control point (GCP) without an obvious phase transition was used to remove the residual terrain phase. The influence of atmospheric phase error was reduced by the spatiotemporal filter. The subsidence rate and time series subsidence results were obtained. Finally, the results were converted to the World Geodetic System 1984 (WGS84) coordinate system by geocoding. The spatial and temporal baselines are shown in Figure 7.Figure 7
Spatial and temporal baselines in this study. The yellow diamond denotes the super master image, the green diamonds denote the slave image, and each blue line denotes a connection from the master image to the slave image.In Hongguang Village ⑤, Zhanhua District ⑥, and Xianhe Town (⑦), the same location of the subsidence areas appeared in two time periods, indicating that the surface subsidence was caused by long-term damage, so this paper focused on these three areas. The subsidence rate results of these three areas are shown in Figure8. We superimposed results on Landsat 8 images with 50% transparency to better display subsidence and surface information. The value is denoted by the color scale from blue to red.Figure 8
Subsidence rate results. (a) Subsidence rate of Zhanhua district ⑥ in Figure5; (b) subsidence rate of Xianhe town ⑦ in Figure 5; (c) subsidence rate of Hongguang village ⑤ in Figure 5. Negative values indicate an increase in the distance along the LOS (subsidence) and positive values present a decrease in the distance along the LOS (uplift).
(a)(b)(c)


#### 4.2.2. Cross Track Validation

In order to evaluate the accuracy of SBAS-InSAR monitoring results, cross orbital validation was used in this study. The subsidence value was obtained by using the 24-scene Sentinel-1B data when monitoring in the same way without changing the parameters. Ten validation points evenly distributed in large subsidence areas were randomly selected in Xianhe Town, as shown in Figure9(a), to check whether the monitoring results of SBAS-InSAR at different orbits were consistent, and root mean square error (RMSE) was adopted as the evaluation index. The subsidence curves of 10 points in different orbits are shown in Figure 9(b). Due to the 6-day time interval for Sentinel-1 image acquisition and the different look direction of ascending orbit and descending orbit, a slight difference was acceptable. After calculation, the RMSE was 14.66 mm, and the experimental results are considered reliable. The equation of RMSE is as follows:(9)RMSE=1n∑i=1ndas,i−ddes,i,where n denotes the number of verification points; i denotes the point of the serial number; das denotes the subsidence rate at the verification point of the ascending image; and ddes denotes the subsidence rate at the point of the descending image.Figure 9
Validation results. (a) Location distribution of 10 validation points. Red circle indicates the validation point, and the number of the point is marked on the right side of circle. (b) Subsidence rate at each point. Black square denotes ascending Sentinel-1 datasets, and red circle denotes descending Sentinel-1 datasets.
(a)(b)


## 4.1. Detection Results

In this paper, two interferometric pairs of the Yellow River Delta from 10 January 2018 to 23 March 2018 and from 06 November 2018 to 24 December 2018 were selected. Due to the serious decorrelation noise generated by vegetation at other times, the subsidence areas lost their characteristics, so the detection effect was not good. Therefore, we no longer showed and listed their results. The parameters of interferograms are shown in Table3.Table 3
Parameters of interferograms for detecting the subsidence areas.
Master imageSlave imagePathSpatial baseline (m)Temporal baseline (d)22/12/2017∗27/01/2018142−80.563615/01/201808/02/2018142−101.2524Subsidence areas were detected in 9 subareas. The interferogram from 10 January 2018 to 23 March 2018 detected 7 subsidence areas as shown in Figure6(a). They were located in Shuinan Village, Zhangqiu District, Jinan ①, Maqiao Town, Huantai County, Zibo ②, Hanting District, Weifang ③, Yangkou Town, Shouguang City, Weifang ④, Hongguang Village, Kenli District, Dongying ⑤, Zhanhua District, Binzhou ⑥, and Xianhe Town, Hekou District, Dongying ⑦. The interferogram from 06 November 2018 to 24 December 2018 detected 6 subsidence areas as shown in Figure 6(b). They were located in Yangkou Town, Shouguang City, Weifang ④), Hongguang Village, Kenli District, Dongying, ⑤, Zhanhua District, Binzhou ⑥), Xianhe Town, Hekou District, Dongying ⑦, south of the Yellow River estuary in Kenli District, Dongying ⑧, and Mazishan Town, Wudi County, Binzhou ⑨. Among them, Hongguang Village ⑤, Zhanhua District ⑥, and Xianhe Town ⑦ had the same location of the subsidence areas in the interferogram of both time periods, indicated the surface subsidence lasted a long time. The detection results of each subregions are shown in Figures 6(a) and 6(b). The lines in Figures 6(a) and 6(b) denote image cropping mode, and the subsidence areas are marked by red boxes. Hongguang Village, Zhanhua District, and Xianhe Town are shown in the enlarged Figures 6(c), 6(d), and 6(e).Figure 6
Detection results of subsidence areas. (a) Detection results from 10/01/2018 to 23/03/2018 interferometric pair; (b) detection results from 13/10/2018 to 24/12/2018 interferometric pair. The results of (c) Hongguang village ⑤, (d) Zhanhua district ⑥, and (e) Xianhe town ⑦ are enlarged, respectively.
(a)(b)(c)(d)(e)


## 4.2. SBAS-InSAR Results

### 4.2.1. Subsidence Rate

The detected subsidence areas were cropped from Sentinel-1A images from January 2018 to December 2018 to form a time series. A maximum temporal baseline was less than 180 days and spatial baseline was shorter than 200 m generated. External registration was carried out with SRTM DEM 90M provided by NASA. Goldstein filter [45] was used to eliminate the noise phase. The minimum cost flow method (MCF) [47] was used for phase unwrapping. The Yellow River Delta was covered with vegetation in many places, so there was a wide range of low coherence areas after SAR interferometry. The unwrapping decomposition level was set to 2 to improve the quality and efficiency of unwrapping and reduce unwrapping errors. The interferometric pairs with low coherence and poor unwrapping were removed. A smooth ground control point (GCP) without an obvious phase transition was used to remove the residual terrain phase. The influence of atmospheric phase error was reduced by the spatiotemporal filter. The subsidence rate and time series subsidence results were obtained. Finally, the results were converted to the World Geodetic System 1984 (WGS84) coordinate system by geocoding. The spatial and temporal baselines are shown in Figure 7.Figure 7
Spatial and temporal baselines in this study. The yellow diamond denotes the super master image, the green diamonds denote the slave image, and each blue line denotes a connection from the master image to the slave image.In Hongguang Village ⑤, Zhanhua District ⑥, and Xianhe Town (⑦), the same location of the subsidence areas appeared in two time periods, indicating that the surface subsidence was caused by long-term damage, so this paper focused on these three areas. The subsidence rate results of these three areas are shown in Figure8. We superimposed results on Landsat 8 images with 50% transparency to better display subsidence and surface information. The value is denoted by the color scale from blue to red.Figure 8
Subsidence rate results. (a) Subsidence rate of Zhanhua district ⑥ in Figure5; (b) subsidence rate of Xianhe town ⑦ in Figure 5; (c) subsidence rate of Hongguang village ⑤ in Figure 5. Negative values indicate an increase in the distance along the LOS (subsidence) and positive values present a decrease in the distance along the LOS (uplift).
(a)(b)(c)


### 4.2.2. Cross Track Validation

In order to evaluate the accuracy of SBAS-InSAR monitoring results, cross orbital validation was used in this study. The subsidence value was obtained by using the 24-scene Sentinel-1B data when monitoring in the same way without changing the parameters. Ten validation points evenly distributed in large subsidence areas were randomly selected in Xianhe Town, as shown in Figure9(a), to check whether the monitoring results of SBAS-InSAR at different orbits were consistent, and root mean square error (RMSE) was adopted as the evaluation index. The subsidence curves of 10 points in different orbits are shown in Figure 9(b). Due to the 6-day time interval for Sentinel-1 image acquisition and the different look direction of ascending orbit and descending orbit, a slight difference was acceptable. After calculation, the RMSE was 14.66 mm, and the experimental results are considered reliable. The equation of RMSE is as follows:(9)RMSE=1n∑i=1ndas,i−ddes,i,where n denotes the number of verification points; i denotes the point of the serial number; das denotes the subsidence rate at the verification point of the ascending image; and ddes denotes the subsidence rate at the point of the descending image.Figure 9
Validation results. (a) Location distribution of 10 validation points. Red circle indicates the validation point, and the number of the point is marked on the right side of circle. (b) Subsidence rate at each point. Black square denotes ascending Sentinel-1 datasets, and red circle denotes descending Sentinel-1 datasets.
(a)(b)


## 4.2.1. Subsidence Rate

The detected subsidence areas were cropped from Sentinel-1A images from January 2018 to December 2018 to form a time series. A maximum temporal baseline was less than 180 days and spatial baseline was shorter than 200 m generated. External registration was carried out with SRTM DEM 90M provided by NASA. Goldstein filter [45] was used to eliminate the noise phase. The minimum cost flow method (MCF) [47] was used for phase unwrapping. The Yellow River Delta was covered with vegetation in many places, so there was a wide range of low coherence areas after SAR interferometry. The unwrapping decomposition level was set to 2 to improve the quality and efficiency of unwrapping and reduce unwrapping errors. The interferometric pairs with low coherence and poor unwrapping were removed. A smooth ground control point (GCP) without an obvious phase transition was used to remove the residual terrain phase. The influence of atmospheric phase error was reduced by the spatiotemporal filter. The subsidence rate and time series subsidence results were obtained. Finally, the results were converted to the World Geodetic System 1984 (WGS84) coordinate system by geocoding. The spatial and temporal baselines are shown in Figure 7.Figure 7
Spatial and temporal baselines in this study. The yellow diamond denotes the super master image, the green diamonds denote the slave image, and each blue line denotes a connection from the master image to the slave image.In Hongguang Village ⑤, Zhanhua District ⑥, and Xianhe Town (⑦), the same location of the subsidence areas appeared in two time periods, indicating that the surface subsidence was caused by long-term damage, so this paper focused on these three areas. The subsidence rate results of these three areas are shown in Figure8. We superimposed results on Landsat 8 images with 50% transparency to better display subsidence and surface information. The value is denoted by the color scale from blue to red.Figure 8
Subsidence rate results. (a) Subsidence rate of Zhanhua district ⑥ in Figure5; (b) subsidence rate of Xianhe town ⑦ in Figure 5; (c) subsidence rate of Hongguang village ⑤ in Figure 5. Negative values indicate an increase in the distance along the LOS (subsidence) and positive values present a decrease in the distance along the LOS (uplift).
(a)(b)(c)


## 4.2.2. Cross Track Validation

In order to evaluate the accuracy of SBAS-InSAR monitoring results, cross orbital validation was used in this study. The subsidence value was obtained by using the 24-scene Sentinel-1B data when monitoring in the same way without changing the parameters. Ten validation points evenly distributed in large subsidence areas were randomly selected in Xianhe Town, as shown in Figure9(a), to check whether the monitoring results of SBAS-InSAR at different orbits were consistent, and root mean square error (RMSE) was adopted as the evaluation index. The subsidence curves of 10 points in different orbits are shown in Figure 9(b). Due to the 6-day time interval for Sentinel-1 image acquisition and the different look direction of ascending orbit and descending orbit, a slight difference was acceptable. After calculation, the RMSE was 14.66 mm, and the experimental results are considered reliable. The equation of RMSE is as follows:(9)RMSE=1n∑i=1ndas,i−ddes,i,where n denotes the number of verification points; i denotes the point of the serial number; das denotes the subsidence rate at the verification point of the ascending image; and ddes denotes the subsidence rate at the point of the descending image.Figure 9
Validation results. (a) Location distribution of 10 validation points. Red circle indicates the validation point, and the number of the point is marked on the right side of circle. (b) Subsidence rate at each point. Black square denotes ascending Sentinel-1 datasets, and red circle denotes descending Sentinel-1 datasets.
(a)(b)


## 5. Discussion

### 5.1. Analysis of Zhanhua District

Figure10 shows the time series of surface deformation in Zhanhua District from 10 January 2018 to 24 December 2018, which reflects the evolution process of this area. 10 January 2018 is considered as the reference time of the time series deformation map, and a total of 23 deformation maps are obtained. Here, 11 of them are selected for display. There was one image every month (except January) to show the deformation characteristics in different periods. The date is in the upper left corner of each map. The deformation value is denoted by the color scale from blue to red.Figure 10
Cumulative deformation map of time series in 2018 in Zhanhua district (with reference to 10 January 2018). The date of obtaining the deformation map is marked in the upper left corner of each map, and the number of the subsidence area is marked in the ellipse. The triangle up in the last image denotes the point where the subsidence curve is plotted.Five major subsidence zones were observed in Zhanhua District, which were defined as ZoneA, Zone B, Zone C, Zone D, and Zone E, respectively. Five zones were marked by red dashed ellipses in Figure 10. By calculating the amount of subsidence in each area, the maximum subsidence values of Zone A, B, C, D, and E were −125.8 mm, −56.3 mm, −76.4 mm, −95.0 mm, and −60.7 mm, respectively. According to the color in the time series cumulative deformation map, except for the 5 subsidence areas, most of the other areas tended to be stable, and the subsidence range was between −10 mm and 20 mm. With the passage of time, the amount of subsidence value of the five subsidence zones increased, and the subsidence range expanded continuously. A point in each of the five subsidence zones was selected to draw the subsidence variation curve as shown in Figure 11. The location of the subsidence point was triangled up in the deformation map of 24 December 2018 in Figure 10.Figure 11
Time series of subsidence at selected points.According to Figure11, we can obtain the cumulative subsidence and the subsidence rate according to the value and the slope change of the line. From January 10 to March 1, the surface subsidence of each zone always kept pace, and the subsidence in Zone A grew the fastest, reaching about 40 mm. From March 11 to March 23, Zone A and C were briefly uplifted, and the subsidence rate slowed down in Zone B, D, and E. From March 23 to August 2, the subsidence rate of Zone A remained fast. Zone B was uplifted from June 3 to June 27 and subsided slowly at other times. The subsidence rate of Zone C was always stable. In Zone D, the subsidence rate accelerated from April 16 to April 28, and then the surface uplift appeared from April 28 to June 4, and the subsidence rate was stable in other times. Zone E kept a stable subsidence rate. From August 2 to September 7, Zone A, B, and E showed obvious surface uplift, while the subsidence rate of Zones C and D slowed down. After September 7, subsidence continued in all areas, with only one uplift occurring in Zone B from October 25 to November 6. It is worth noting that the subsidence points in Zones A, B, and E had obvious surface uplift from August 2 to September 7.From July to August, surface uplift and subsidence rate slowed down in many zones. The climate of the Yellow River Delta determined that the precipitation from July to August was generally greater than other months, and the precipitation could recharge groundwater to slow subsidence. To verify this possibility in this area, the precipitation data of 2018 were acquired from the National Meteorological Science Data Center (https://data.cma.cn/) for analysis, as shown in Figure 12(a). The average precipitation in August was more than 300 mm in 2018, much higher than in other months.Figure 12
Analysis of subsidence in Zhanhua district combined with precipitation. (a) Monthly average precipitation curve in Zhanhua district, red square is precipitation in August. The deformation of 353 points in (b) January, (c) May, (d) August, and (e) November. The red line is the average deformation of 353 points in each month.
(a)(b)(c)(d)(e)When the amount of subsidence is large, it can better reflect the characteristics of subsidence. Therefore, 353 points with a subsidence rate greater than −50 mm/a were selected in the five zones, and the deformation of these points in different months was acquired to draw a scatter diagram. Figures12(b)–12(e), respectively, show the deformation of these points in January, May, August, and November. The red line denotes the average deformation of these points. It was −24.31 mm in January, −9.22 mm in May, −3.57 mm in August, and −7.34 mm in November. In August, the average subsidence value was the smallest, and only in August, there were many points with deformation greater than 0 mm, which was significantly different from other months, confirming that heavy rainfall in August slowed down the subsidence.


### 5.2. Analysis of Xianhe Town

The maximum subsidence rate of Xianhe Town was −330.9 mm/a. In order to further explore the characteristics of subsidence, two profile lines AA′ and BB′ were established for the subsidence area. The cumulative surface subsidence of 30 uniformly distributed points in each profile line was extracted to draw deformation curves, as shown in Figure13(c). It can be seen that the subsidence of the two profile lines is not uniform. The subsidence on the side close to Dongying Shugang Expressway was more serious than that on the side far from the expressway.Figure 13
Subsidence analysis map of Xianhe town. (a) Enlarged plot of settlement rate. AA′ and BB′ are profile lines, the black and yellow lines are Shugang expressway, and the areas with subsidence rates greater than −300 mm/a, −250 mm/a, and −150 mm/a are denoted by red, green, and blue vectors, respectively. (b) The vector denotes the range of cumulative subsidence greater than 20 mm (with reference to 10 January 2018). The vectors of brown, dark green, purple, and light green are the cumulative subsidence in February 15, April 10, August 14, and December 24. (c) Profile maps based on AA′ and BB′. (d) The fitted line of the distance between the subsidence boundary and the Shugang expressway and month.
(a)(b)(c)(d)The Shugang Expressway was an important hub between Hekou District and the port. The traffic flow was large, so it was very important to ensure the safety of the Shugang Expressway. In Figure13(a), we set three vector boundaries of subsidence rate that the value is greater than −300 mm/a, −250 mm/a, and −150 mm/a, which are represented in wine red, green, and blue. It can be seen that some traffic positions of Shugang Expressway have crossed the blue boundary, denoting that the subsidence rate of this part has been greater than −150 mm/a.According to the standard for dangerous building appraisal in China, when the surface subsidence of infrastructure was more than 20 mm, the ground facilities would suffer safety risks. So we used the time series deformation to draw the boundary of 20 mm subsidence at different times and measured the vertical distance between the boundary and the expressway by taking the point at the outer end of the boundary once a month. The linear fitting function between distance and month was established. Distance measurement is shown in Figure13(b). Four of these measurements are shown here. Brown, dark green, purple, and light green vector, respectively, denote the boundary of 20 mm subsidence on February 15, April 10, August 14, and December 24, respectively.The black line denotes distance, and values are marked in the figure. The linear fitting function is shown in Figure13(d), and the R2 (Relevance) is 0.9739. The distance and time almost meet the linear relationship. This proved that the expansion rate of the subsidence range was stable and there was no obvious slowing trend. Subsequently, it was necessary to strengthen governance to slow down the continuous occurrence of surface subsidence so as to reduce the expressway’s potential risk.


### 5.3. Analysis of Hongguang Village

There were three subsidence areas in Hongguang Village, and the maximum subsidence rate was -209.68 mm/a. Combined with the Landsat 8 image, two of them were located in the north, about 5 km and 1.5 km away from the Hongguang Village building, respectively, which were defined as ZoneA and Zone B. The other was located in the south, about 3 km away from the building, which was defined as Zone C. We marked the geographical locations of buildings with red boxes in the enlarged Landsat 8 image, as shown in Figure 14(a) and overlaid the red boxes on the enlarged subsidence rate map, as shown in Figure 14(b), to reflect the location relationship between subsidence and Hongguang village buildings.Figure 14
Subsidence analysis map of Hongguang village. (a) Enlarged landsat 8 image of Hongguang village, the buildings are marked by red line; (b) the enlarged cumulative map of subsidence. The three subsidence areas are marked withA, B, and C. The profile line is DD′. (c) Profile map. The location of the points in the building area is boxed on the map.
(a)(b)(c)The profile line DD′ through Hongguang Village was established as shown in Figure14(b), and the cumulative surface subsidence of 30 uniformly distributed points on the profile line was extracted to draw a profile map. The profile map is shown in Figure 14(c).We marked the points in the location of Hongguang Village buildings with a red box in Figure14(c). The surface subsidence value of these buildings was small. However, because it was located between subsidence Zone B and Zone C, when the subsidence range on both sides expanded continuously, the surface around Hongguang village buildings would sink. Before reaching a stable state, the subsidence would gradually extend to Hongguang Village, causing safety risks to the building.


## 5.1. Analysis of Zhanhua District

Figure10 shows the time series of surface deformation in Zhanhua District from 10 January 2018 to 24 December 2018, which reflects the evolution process of this area. 10 January 2018 is considered as the reference time of the time series deformation map, and a total of 23 deformation maps are obtained. Here, 11 of them are selected for display. There was one image every month (except January) to show the deformation characteristics in different periods. The date is in the upper left corner of each map. The deformation value is denoted by the color scale from blue to red.Figure 10
Cumulative deformation map of time series in 2018 in Zhanhua district (with reference to 10 January 2018). The date of obtaining the deformation map is marked in the upper left corner of each map, and the number of the subsidence area is marked in the ellipse. The triangle up in the last image denotes the point where the subsidence curve is plotted.Five major subsidence zones were observed in Zhanhua District, which were defined as ZoneA, Zone B, Zone C, Zone D, and Zone E, respectively. Five zones were marked by red dashed ellipses in Figure 10. By calculating the amount of subsidence in each area, the maximum subsidence values of Zone A, B, C, D, and E were −125.8 mm, −56.3 mm, −76.4 mm, −95.0 mm, and −60.7 mm, respectively. According to the color in the time series cumulative deformation map, except for the 5 subsidence areas, most of the other areas tended to be stable, and the subsidence range was between −10 mm and 20 mm. With the passage of time, the amount of subsidence value of the five subsidence zones increased, and the subsidence range expanded continuously. A point in each of the five subsidence zones was selected to draw the subsidence variation curve as shown in Figure 11. The location of the subsidence point was triangled up in the deformation map of 24 December 2018 in Figure 10.Figure 11
Time series of subsidence at selected points.According to Figure11, we can obtain the cumulative subsidence and the subsidence rate according to the value and the slope change of the line. From January 10 to March 1, the surface subsidence of each zone always kept pace, and the subsidence in Zone A grew the fastest, reaching about 40 mm. From March 11 to March 23, Zone A and C were briefly uplifted, and the subsidence rate slowed down in Zone B, D, and E. From March 23 to August 2, the subsidence rate of Zone A remained fast. Zone B was uplifted from June 3 to June 27 and subsided slowly at other times. The subsidence rate of Zone C was always stable. In Zone D, the subsidence rate accelerated from April 16 to April 28, and then the surface uplift appeared from April 28 to June 4, and the subsidence rate was stable in other times. Zone E kept a stable subsidence rate. From August 2 to September 7, Zone A, B, and E showed obvious surface uplift, while the subsidence rate of Zones C and D slowed down. After September 7, subsidence continued in all areas, with only one uplift occurring in Zone B from October 25 to November 6. It is worth noting that the subsidence points in Zones A, B, and E had obvious surface uplift from August 2 to September 7.From July to August, surface uplift and subsidence rate slowed down in many zones. The climate of the Yellow River Delta determined that the precipitation from July to August was generally greater than other months, and the precipitation could recharge groundwater to slow subsidence. To verify this possibility in this area, the precipitation data of 2018 were acquired from the National Meteorological Science Data Center (https://data.cma.cn/) for analysis, as shown in Figure 12(a). The average precipitation in August was more than 300 mm in 2018, much higher than in other months.Figure 12
Analysis of subsidence in Zhanhua district combined with precipitation. (a) Monthly average precipitation curve in Zhanhua district, red square is precipitation in August. The deformation of 353 points in (b) January, (c) May, (d) August, and (e) November. The red line is the average deformation of 353 points in each month.
(a)(b)(c)(d)(e)When the amount of subsidence is large, it can better reflect the characteristics of subsidence. Therefore, 353 points with a subsidence rate greater than −50 mm/a were selected in the five zones, and the deformation of these points in different months was acquired to draw a scatter diagram. Figures12(b)–12(e), respectively, show the deformation of these points in January, May, August, and November. The red line denotes the average deformation of these points. It was −24.31 mm in January, −9.22 mm in May, −3.57 mm in August, and −7.34 mm in November. In August, the average subsidence value was the smallest, and only in August, there were many points with deformation greater than 0 mm, which was significantly different from other months, confirming that heavy rainfall in August slowed down the subsidence.


## 5.2. Analysis of Xianhe Town

The maximum subsidence rate of Xianhe Town was −330.9 mm/a. In order to further explore the characteristics of subsidence, two profile lines AA′ and BB′ were established for the subsidence area. The cumulative surface subsidence of 30 uniformly distributed points in each profile line was extracted to draw deformation curves, as shown in Figure13(c). It can be seen that the subsidence of the two profile lines is not uniform. The subsidence on the side close to Dongying Shugang Expressway was more serious than that on the side far from the expressway.Figure 13
Subsidence analysis map of Xianhe town. (a) Enlarged plot of settlement rate. AA′ and BB′ are profile lines, the black and yellow lines are Shugang expressway, and the areas with subsidence rates greater than −300 mm/a, −250 mm/a, and −150 mm/a are denoted by red, green, and blue vectors, respectively. (b) The vector denotes the range of cumulative subsidence greater than 20 mm (with reference to 10 January 2018). The vectors of brown, dark green, purple, and light green are the cumulative subsidence in February 15, April 10, August 14, and December 24. (c) Profile maps based on AA′ and BB′. (d) The fitted line of the distance between the subsidence boundary and the Shugang expressway and month.
(a)(b)(c)(d)The Shugang Expressway was an important hub between Hekou District and the port. The traffic flow was large, so it was very important to ensure the safety of the Shugang Expressway. In Figure13(a), we set three vector boundaries of subsidence rate that the value is greater than −300 mm/a, −250 mm/a, and −150 mm/a, which are represented in wine red, green, and blue. It can be seen that some traffic positions of Shugang Expressway have crossed the blue boundary, denoting that the subsidence rate of this part has been greater than −150 mm/a.According to the standard for dangerous building appraisal in China, when the surface subsidence of infrastructure was more than 20 mm, the ground facilities would suffer safety risks. So we used the time series deformation to draw the boundary of 20 mm subsidence at different times and measured the vertical distance between the boundary and the expressway by taking the point at the outer end of the boundary once a month. The linear fitting function between distance and month was established. Distance measurement is shown in Figure13(b). Four of these measurements are shown here. Brown, dark green, purple, and light green vector, respectively, denote the boundary of 20 mm subsidence on February 15, April 10, August 14, and December 24, respectively.The black line denotes distance, and values are marked in the figure. The linear fitting function is shown in Figure13(d), and the R2 (Relevance) is 0.9739. The distance and time almost meet the linear relationship. This proved that the expansion rate of the subsidence range was stable and there was no obvious slowing trend. Subsequently, it was necessary to strengthen governance to slow down the continuous occurrence of surface subsidence so as to reduce the expressway’s potential risk.


## 5.3. Analysis of Hongguang Village

There were three subsidence areas in Hongguang Village, and the maximum subsidence rate was -209.68 mm/a. Combined with the Landsat 8 image, two of them were located in the north, about 5 km and 1.5 km away from the Hongguang Village building, respectively, which were defined as ZoneA and Zone B. The other was located in the south, about 3 km away from the building, which was defined as Zone C. We marked the geographical locations of buildings with red boxes in the enlarged Landsat 8 image, as shown in Figure 14(a) and overlaid the red boxes on the enlarged subsidence rate map, as shown in Figure 14(b), to reflect the location relationship between subsidence and Hongguang village buildings.Figure 14
Subsidence analysis map of Hongguang village. (a) Enlarged landsat 8 image of Hongguang village, the buildings are marked by red line; (b) the enlarged cumulative map of subsidence. The three subsidence areas are marked withA, B, and C. The profile line is DD′. (c) Profile map. The location of the points in the building area is boxed on the map.
(a)(b)(c)The profile line DD′ through Hongguang Village was established as shown in Figure14(b), and the cumulative surface subsidence of 30 uniformly distributed points on the profile line was extracted to draw a profile map. The profile map is shown in Figure 14(c).We marked the points in the location of Hongguang Village buildings with a red box in Figure14(c). The surface subsidence value of these buildings was small. However, because it was located between subsidence Zone B and Zone C, when the subsidence range on both sides expanded continuously, the surface around Hongguang village buildings would sink. Before reaching a stable state, the subsidence would gradually extend to Hongguang Village, causing safety risks to the building.


## 6. Conclusion

In this paper, we constructed the CNSBAS-InSAR method to monitor and analyze the surface subsidence in the Yellow River Delta. This method first uses the CenterNet network to automatically detect the subsidence areas from the wide differential interferogram and determine the subsidence area. The main purpose of this step was to determine the location to crop the image, so as to save time and reduce the atmospheric phase error. Then we conducted SBAS-InSAR processing on the small-scale subsidence area to obtain the subsidence results. Based on 24 scene Sentinel-1A images from 10 January 2018 to 24 December 2018, the subsidence location and small-scale subsidence rate of the Yellow River Delta were obtained. Finally, based on precipitation data, Landsat 8 image, highway vector, and building, we analyzed the causes and potential risks. Through the study in this paper, some valuable conclusions can be drawn as follows:(1)
9 subsidence areas were detected in the Yellow River Delta. Three of them had long-term surface subsidence. They were located in Zhanhua District, Xianhe Town, and Hongguang Village.(2)
There were 5 obvious subsidence areas in Zhanhua District, and the maximum subsidence rate was −135.21 mm/a. The surface subsidence in this area was closely related to precipitation. Heavy rainfall can reduce the amount of surface subsidence.(3)
The maximum subsidence rate was −330.91 mm/a in Xianhe Town. It was worth noting that the subsidence rate of part of the Shugang Expressway exceeded −150 mm/a. The distance between the subsidence area boundary and the expressway and time almost met the linear function (R2 = 0.9739), which indicated that the subsidence area was still expanding without a slowing downtrend.(4)
There were three subsidence zones in Hongguang Village, and the maximum subsidence rate was −209.68 mm/a. In 2018, the subsidence value of Hongguang village buildings was small. However, through the analysis of the profile line of subsidence, the subsidence on both sides of the building posed a great threat to its safety.Overall, CNSBAS-InSAR accurately detected the subsidence location of the Yellow River Delta, and then the small-scale SBAS-InSAR processing was conducted to obtain a millimeter-level surface subsidence value of the subsidence area. The method in this paper plays a positive role in short time and efficient SBAS-InSAR monitoring. And the analysis results can also provide a reference for local construction and protection in Yellow River Delta, but this study also has some limitations. Some error detection and leak detection occurred when using CenterNet to detect subsidence areas. This requires continuous improvement of the network and data set. In addition, the detection capability of CenterNet at different noise levels is also worth exploring to meet the requirements of different spatial and temporal baselines for SBAS-InSAR processing. We will solve these above problems in the future work.



---

*Source: 2672876-2022-05-24.xml*