Understanding the Spatial-Temporal Characteristics of Land Subsidence in Shenzhen under Rapid Urbanization Based on MT-InSAR

Land subsidence is a common geological hazard in urban areas. It can be caused by human activities such as land reclamation and subway construction. However, rapid urbanization can accelerate the progress of these activities. Therefore, spatial-temporal analysis of land subsidence can provide a guarantee for rapid urbanization progress. In this study, Shenzhen, one of the fastest urbanizing cities in China, is selected for continued land subsidence monitoring. We collected 153 Sentinel-1A images from December 2015 to July 2021 and utilized an interferometric point target analysis method to generate the land subsidence in Shenzhen. Additionally, we analyzed the relationship between land subsidence and typical urbanization progress: land reclamation, transportation network, and urban construction. The results show that land subsidence mostly occurs in the western reclamation area, with a maximum rate of −74.9 mm/a in Qianhai Bay. Areas with larger quaternary marine sediment thickness in the range of 10–20 m are more likely to experience subsidence. Obvious subsidence (<−20 mm/a) along the subway mainly occurred on Lines 11 and Line 5, which is inversely proportional to the distance from the metro. The correlation coefficient is −0.197. Areas with low road densities have a high probability of experiencing great subsidence. A relationship between land subsidence and building construction in the Guangming district is found. Obvious subsidence (<−20 mm/a) mainly occurred in the building areas, with a percentage of 43.48%.

affected by land subsidence, China is experiencing geological hazards in 17 provinces, affecting an area of 70 000 km 2 [1], [2]. However, the current situation of land subsidence in China deteriorates with urbanization progress, especially in cities with rapid urbanization. Additionally, in addition to land subsidence, other secondary disasters, such as sinkholes [3], urban flooding [4], infrastructure damage [5], and seawater intrusion in coastal areas [6], [7] can threaten people's lives and the normal operation of cities. Specific to the city of Shenzhen, it has undergone tremendous urbanization with the rapid development of the Guangdong-Hong Kong-Macao Greater Bay Area (GBA) in the past few decades [8]. Therefore, long-term monitoring of land subsidence and analysis of the relationship between land subsidence and urbanization are of great significance to the healthy development of cities.
Multitemporal InSAR (MT-InSAR) is a commonly used method for deformation monitoring that was developed from D-InSAR technology after overcoming the limitations caused by temporal-spatial decorrelation and atmospheric effects. It has been widely used in the fields of deformation monitoring for infrastructure [9], [10], [11], [12], mining [13], [14], landslides [15], [16], urban areas [17], [18], [19], earthquakes [20], [21], and volcanoes [22], [23], [24]. Among these applications, land subsidence under urbanization has been studied extensively. These studies mainly involved in spatial and temporal characteristics of ground deformation, cause analysis of land subsidence, and relationship exploration between land subsidence and building properties in the progress of urbanization. For example, the ground deformation monitoring under the project of "Mountain Excavation and City Construction" (MECC) in the Yan'an New District (YND) [25], [26], [27], [28]. Investigate the temporal and spatial patterns of land subsidence in Lanzhou New District and Ho Chi Minh City as a result of land reclamation and urban expansion [29], [30], [31], [32]. To generate the connection between the ground load and land subsidence, Chen et al. [33] and Zhou et al. [34] used the building IBI index and machine learning techniques to analyze the relationship between land subsidence and causes, such as static and dynamic load, groundwater level, and sediment thickness. In addition, other factors, including groundwater extraction, building construction, geological strata difference, and seasonal precipitation can also contribute to urban land subsidence [35], [36], [37], [38]. In Pisa, Italy, Solari et al. studied the connection between land subsidence and the age of buildings [39]. Yang et al. showed a weak correlation between ground subsidence and building volume in the Chaoyang and Tongzhou districts of Beijing [40].
In the case of Shenzhen, the hot spots of research are land reclamation and underground space excavation. For example, Xu et al. confirmed a strong connection between land reclamation and subsidence in the coastal area of Shenzhen [41], [42]. Additionally, subway construction is a factor that can contribute to land subsidence, according to Hu et al. [43] and Liu et al. [44]. He et al. employed InSAR and LiDAR data to study the connection between ground subsidence and buildings in Shenzhen [45]. The findings revealed that subsidence was primarily brought on by new buildings. However, few studies have linked surface subsidence to urbanization characteristics in Shenzhen thus far. Therefore, three typical urbanization processes, land reclamation, transportation network including subway network and ground road network, and urban construction were adopted to analyze the relationship between urbanization and land subsidence. The findings of this work can provide valuable information for urban planning and key points of potential geologic hazard monitoring of similar cities.
In this study, we calculated the time series deformation of Shenzhen in the last six years based on the IPTA-SBAS method. Then, we analyzed the characteristics of subsidence with external datasets, such as optical images, geological data, and road data. We generated the relationship between land subsidence and urbanization with GIS spatial analysis and a correlation index.

A. Regional Overview
Shenzhen, one of the largest cities in China, is located in the south of Guangdong Province (22°26 -22°51 N, 113°45 -114°37 E). As shown in Fig. 1, the entire city has a high northeast and low southwest topography. In the south and west of the city, it is surrounded by the ocean with an elevation between 1-10 m, but the average elevation in the northeast reaches 80 m because it is covered by mountains there. The geological structure of Shenzhen is complex, and the fracture structure is the main component of the city. There are three primary faults crossing the city: the Henggang fault, the Wentang-Guanlan fault, and the Wuhua-Shenzhen fault [46].
In terms of urban development, Shenzhen is one of the Special Economic Zones (SEZ) in China after the reform and opening-up policy. It has experienced a series of urbanization processes since 1980. We selected three typical processes, land reclamation, development of transportation networks, and urban construction, to analyze the evolution of urbanization in Shenzhen. Land reclamation is one of the important methods used to expand the area of the city and accelerate urbanization. According to the statistics, Shenzhen's reclaimed land area exceeds 74.77 km 2 , and it will continue to increase in accordance with the Five-Year Action Plan for Urban Infrastructure Construction in Shenzhen (2021-2025). Another aspect of urbanization is the construction of the subway network. Currently, 11 subway lines have been completed and put into use. By 2025, another five subway lines will be planned, and the accumulated subway miles will reach 647 kilometers. As for urban construction, a mass of mountains has been excavated to build infrastructure in Guangming and Longgang districts in recent years, especially when the policy "Several Opinions on Supporting Guangming Science City to Build a World-class Science City" was proposed in 2019. Besides, to get a clear understanding of the urbanization process in Shenzhen, we also collected the landmark events of each typical urbanization process, as shown in Table S1.

B. Datasets
We collected 153 Sentinel-1A images from December 2015 to July 2021 to obtain the time series deformation of Shenzhen. These SAR images are retrieved from ASF DAAC, and processed by ESA. The orbit number is 40 858 and the number of two adjacent frames is 65 and 71. The coverage of these images is shown in Fig. 1. The 1-arc Shuttle Radar Topography Mission (SRTM) is used as an external DEM. We also collected datasets from five GPS stations located in Hong Kong to evaluate the accuracy of the InSAR results. Additionally, we also collected some geological datasets [47], historical optical images (from Google Earth), and road network data of Shenzhen (https://www.ritu.cn/) to analyze the relationship between land subsidence and urbanization.

C. IPTA-SBAS Method
Interferometric point target analysis (IPTA) is a commonly used MT-InSAR method to generate time series deformation. Different from the traditional MT-InSAR methods, such as SBAS-InSAR [48] and PS-InSAR [49], IPTA-SBAS adopts an iterative method to refine the coherence points and estimate the parameters simultaneously. First, the IPTA-SBAS method generates interferometric pairs with temporal and perpendicular baseline thresholds. The differential interferograms are obtained with standard D-InSAR procedures [50]. Second, an iterative regression analysis is utilized to compute the residual topography. For each iteration, the residual phase is unwrapped with an MCF method, and the coherent candidates are refined with a given threshold for the residual phase. Then, linear deformation is calculated based on the unwrapped phase after removing the residual topographic phase. Finally, the nonlinear deformation is separated from the atmospheric signal with a spatial-temporal filtering procedure. The final time series deformation is generated by adding the nonlinear deformation and the linear deformation. Fig. 2 illustrates the flowchart of IPTA-SBAS.
Specifically, IPTA-SBAS is executed in the following three steps.
Step 1: Differential interferometric phase generation. In this study, 2017 119 was chosen as the master image, and slave images were coregistered to the master image. A total of 446 interferometric pairs were selected with temporal and perpendicular baseline thresholds of 100 days and 150 m, respectively. The topographical phase was removed with the 1-arc SRTM DEM. A multilook operation, 8 and 2 in the range and azimuth directions, respectively, was adopted to reduce phase noise.
Step 2: Iterative regression analysis. The residual topographical error phase was estimated by three separated parts. The first part, denoted as ΔZ 1 , was calculated by 1-D regression analysis. The unknown parameter is the residual topographic phase, and the observation is the original wrapped differential phase. A threshold of residual phase, 1.3 rad in this study, was adopted to improve the accuracy of the parameter and refine the coherence candidates [51]. Additionally, the wrapped residual phase was unwrapped with the MCF method, and part of the atmospheric phase was estimated by spatial filtering. The window size of the filtering is set as 1600 m empirically. The second part, denoted as ΔZ 2 , was estimated with the abovementioned iterative 1-D regression analysis. However, the observation is the remaining phase after subtracting the partial atmospheric phase and residual topographic phase from the unwrapped phase. The third part, denoted as ΔZ 3 , was calculated by the 2-D regression analysis method. The unknown parameters are the third part of the residual topographical phase and the linear deformation rate.
Step 3: Separation between nonlinear deformation and atmospheric signals. The atmospheric signal was estimated by high-pass temporal filtering following low-pass spatial filtering. The windows for temporal and spatial filtering are 180 days and 800 m [52], respectively. Nonlinear deformation was generated after removing the atmospheric phase from the residual unwrapped phase. The final time series deformation was the sum of the linear deformation and nonlinear deformation.

III. RESULTS AND ANALYSIS
The average deformation rate of Shenzhen was calculated by the IPTA-SBAS, as shown in Fig. 3. Considering the weak horizontal deformation in Shenzhen [41], we used the deformation in the line-of-sight (LOS) direction to analyze the spatialtemporal characteristics in the following sections. Positive value represents movement toward the satellite (uplift), while negative value represents movement away from the satellite (subsidence). The deformation magnitude of the inland is gentle, ranging from −15 mm/a to 5 mm/a. The subsidence is mainly distributed along the coastal areas with a maximum rate of 74.9 mm/a. Compared with the eastern area of the city, obvious subsidence is mainly distributed on the western coastline area.

A. Characteristics of Subsidence in the Land Reclamation Area
Land reclamation is a significant urbanization symbol in Shenzhen. As shown in Fig. 3(a), the changes in the coastline between 1979 and 2021 represent the land reclamation area, where land subsidence is mainly concentrated. To analyze the characteristics of subsidence, five areas with subsidence rates over 15 mm/a were selected, as shown in Fig. 3 Area S1 is located in the northwestern part of Shenzhen. The subsidence in this area is mainly distributed in the land reclamation area (zone 1) and on both sides of the Dongbao River (D.R) (zone 2), as shown in Fig. 3(b). The largest settlement occurred in the Shenzhen Convention and Exhibition Center (zone 3), with a maximum subsidence rate of 42.0 mm/a. Area S2 is located around Baoan International Airport (B.I.A), with three obvious subsidence zones [zones 1-3 in Fig. 3(c)]. Zone 1 and Zone 2 are within Baoan International Airport, whose maximum subsidence rates reach 39.8 mm/a and 23.3 mm/a, respectively. Zone 3 is located outside the airport with a maximum rate of 66.2 mm/a. In Nanshan District, subsidence is mainly located in Qianhai Bay (Q.B) and Houhai (H.H), marked as area S3 and area S4 in Fig. 3(a). Subsidence in Qianhai Bay is concentrated in three zones, as shown in Fig. 3(d). Zone 1 is located at Dachanwan wharf (D.W), and Zone 2 and Zone 3 are located between the Liantong canal (L.C) and Shuangjie River (S.R). The maximum subsidence rates are 38.2 mm/a, 29.8 mm/a, and 35.2 mm/a, respectively. Unlike the case in Qianhai Bay, subsidence in Houhai Bay is mainly located along the roads, as shown in Fig. 3(e). The maximum subsidence rates of Zone 1 and Zone 2 are 44.3 mm/a and 28.7 mm/a, respectively.
Area S5 is situated in Yantian Port (Y.P), one of the major ports in Shenzhen. Two obvious subsidence zones are identified, as shown in Fig. 3(f). Zone 1 is near the Yantianlu subway station (Y.S), with a maximum settlement rate of 22.1 mm/a. Zone 2 is within the land reclamation area, and the maximum subsidence rate is approximately 17.8 mm/a. Compared with the other four areas, the subsidence in Area S5 is the smallest.

B. Accuracy Verification
We collected five GPS stations (HKTK, HKSL, HKLT, HKKT, and HKFN) derived from the Nevada Geodetic Laboratory (NGL) to verify the results of IPTA-SBAS [53]. GNSS datasets derived from the NGL are usually used as the true value of surface deformation [54], [55]. The GPS results have been translated to the line of sight (LOS) direction of SAR images according to the equations described in [6]. Then, the relative values of InSAR were calibrated to the absolute values according to the reference HKTK station. For each GPS station, InSAR measurement points located within a circle with a radius of 100 m are selected to calculate the mean deformation value, which is a proxy for the comparison with two measurements, as shown in Fig. 4 [56]. The root mean square error (RMSE) between the two results for the four stations is 6.9 mm, 8.0 mm, 6.5 mm, and 6.8 mm, which evaluates the accuracy of the IPTA-SBAS results.

C. Characteristics of Subsidence Along the Subway Line
Underground exploitation, especially the construction and operation of subways, is a common method for urban traffic expansion and a symbol of the urbanization process. However, overexploitation in urban areas may lead to geological hazards such as land subsidence and sinkholes [57], [58]. In this study, 11 subway lines with a 500 m buffer zone of Shenzhen are analyzed, and the deformation results are shown in Fig. 5. We assumed that subsidence with a rate over 10 mm/a was mainly caused by subway tunneling. Four representative areas with subsidence rates exceeding 20 mm/a were identified; see the enlarged views in Fig. 5. The maximum subsidence rates of area S1 to area S4 are 34 mm/a, 66.2 mm/a, 35.5 mm/a. and  20.8 mm/a, respectively. Three of the four areas (S2, S3, and S4) are distributed along the coastal area, which is consistent with the overall regional subsidence characteristics. In addition, the area with the largest subsidence is located at Line 11, with a maximum subsidence rate of 71.5 mm/a.
To further analyze the subsidence characteristics of the subway lines, we counted the number of coherent points, the maximum subsidence rate, the percentage of coherent points with a subsidence rate greater than 10 mm/a, the location, the completion time, and the operation time of each line. We ranked the 11 subway lines in descending order of the coherent point percentage, as shown in Table I. Four subway lines with maximum subsidence rates exceeding 40 mm/a, that is, Lines 2, Line 5, Line 10, and Line 11. The corresponding maximum rates of the four lines are 47.3 mm/a, 53.5 mm/a, 40.5 mm/a, and 71.5 mm/a, respectively. Additionally, the top two subway lines in the coherent point percentage are Line 5 and Line 11, with values of 2.6% and 3.2%, respectively. Line 11 is located along the coastal area, while Line 5 is located inland. Considering the subsidence and location of the subway lines, we selected Line 5 and Line 11 for a more comprehensive analysis.
The subsidence along Line 5 is concentrated in two zones: one is located in Qianhai Bay, while the other is located between Shenzhen North station and Shangshuijing station, as shown in the enlarged view of Fig. 6(a). In Zone 1, we detected four areas with obvious subsidence funnels to analyze the subsidence features along the subway [see Fig. 6(b)-(e)]. In Fig. 6(b), (d), (e)), the subsidence funnels are close to the subway station, and the maximum subsidence rate of each bowl was quite similar, approximating 16 mm/a, 15.2 mm/a, and 17.5 mm/a, respectively. In Fig. 6(c), the subsidence was located in a parking lot near the Pearl River Delta Ring Expressway with a maximum subsidence rate of 29.4 mm/a. The average deformation rate of this subsidence bowl is over 20 mm/a, which is the most serious among the four funnels. It is noted that the subsidence in Qianhai Bay (Zone 2) will be analyzed in the following section due to the overlap between Line 5 and Line 11.
Line 11 was constructed from April 2012 and was put into operation in June 2016. Different from the centralized subsidence of Line 5, the subsidence distribution of Line 11 is dispersed. Five areas with subsidence rates greater than 20 mm/a were observed near Songgang station (area S1), Houting station (area S2), Jinwan Avenue (area S3), Qianhaiwan station (area S4), Hongshuwannan station, and Shenwan station (area S5), as shown in Fig. 7(a). In areas S1, S2, S4, and S5, subsidence was mainly gathered around subway stations, and the maximum subsidence rates of the four areas are 33.2 mm/a, 22.9 mm/a, 34.0 mm/a, and 20 mm/a, respectively. In addition to the subsidence near the subway station, another subsidence funnel within area S3 is distributed along Line 11. The subsidence funnel is located in an industrial park, with mean and maximum subsidence rates of 29.2 mm/a and 71.5 mm/a, respectively. Among these five areas, the subsidence in area S3 is the most severe, and the largest total subsidence area occurs in area S4. Moreover, fieldwork was carried out in area S3, which is located at the south of Baoan International Airport [point TP1 in Fig. 7(a)], as shown in Fig. 8. Obvious subsidence (approximate 8 cm) of the building was observed, which can also evaluate the reliability of InSAR results.
To analyze the temporal feature of subsidence along Line 11, we plotted the time series deformation of six typical points (P1-P6) in area S3 and area S4 by considering the magnitude  and area of subsidence. The deformation trend of P1-P3 in area S3 was nearly linear, as shown in Fig. 7(b), (d), (f). The distance from P1, P2, and P3 to Line 11 are 153 m, 248 m, and 349 m, respectively. The cumulative settlements of P1 and P2 are 16 cm and 8 cm, respectively, while the cumulative settlement of P3 was almost 0 cm, which is inversely proportional to the distance along the subway. Contrary to the characteristics of P1-P3 in area S3, the time series deformation of P4-P6 in area S4 shows an obvious periodic trend [see Fig. 7(c), (e), (g)]. The cumulative settlements of P4, P5, and P6 are 14.27 cm, 16.1 cm, and 3.6 cm, respectively. The reason for these two different deformation trends is probably caused by the intervention of human activities in the coastal land reclamation area; that is, the linear deformation in area S3 may be in the initial stage of soft soil compaction; however, the periodic deformation in area S4 probably related to the planned construction work around the commercial center of Qianhai Bay.
In addition, we calculated the correlation between the coherent point percentage and the operation time, as shown in Fig. S1. The result is −0.08, indicating that there is no obvious relationship between operation time and subsidence in the case of Shenzhen subway.

D. Characteristics of Subsidence in Guangming District
Guangming district was approved in 2018, and it is one of the newest districts of Shenzhen. As an important node of the  Guangzhou-Shenzhen Science and Technology Corridor and an important area of the Guangdong-Hong Kong-Macao Bay Area (GBA), Guangming District has carried out much work to promote urbanization, such as infrastructure construction and land excavation, in the past few years. Therefore, we choose the Guangming district as a representative for the urbanization of Shenzhen. The deformation result of the Guangming district is shown in Fig. 9(a). The subsidence is mainly distributed in the south and northeast of the Guangming district. Three areas with obvious subsidence are selected to analyze the causes of the deformation, marked as S1, S2, and S3 in Fig. 9(a). In area S1, the subsidence is mainly located in the development zone that is derived from vegetation and agricultural land, as shown in Fig. 9(b). The subsidence is probably caused by the engineering construction. Unlike the situation in area S1, the subsidence in area S3 is mainly distributed in the industrial and residential areas, as shown in Fig. 9(d). The discrete subsidence is most likely caused by groundwater exploitation for industrial and domestic water. In area S2, the subsidence is mainly concentrated on the Gongming Reservoir, as shown in Fig. 9(c). To better understand the subsidence process of the reservoir, we divided the dam into three parts according to the chronological order of construction. As we expected, the subsidence is inversely proportional to the construction time; that is, Part 3 experiences the largest settlement, Part 2 follows, and Part 1 is the last. The cumulative settlement in Part 2 and Part 3 reached 8.2 cm and 23.6 cm, respectively.

A. Causes of Land Subsidence in Qianhai Bay
Qianhai Bay is a typical reclamation area with obvious subsidence, as shown in Fig. 10(a), which is the future financial center of Shenzhen. The reclaimed land area exceeds 7.7 km 2 . To understand the causes of subsidence, we collected sediment thickness data of Qianhai Bay, as shown in Fig. 10(b). The sedimentary layer has the characteristics of large porosity, high water content, and high compressibility. Ground with a thick sedimentary layer is more prone to deformation [59]. We found that subsidence is mainly situated in the zone, where the sediment thickness is between 10 and 20 m. Some discrete settlements were also found in areas with sediment thickness between 3 and 10 m. It is noted that areas with sediment thickness below 3 m are located in mountainous areas, Nanshan Park of Shenzhen, resulting in low point density of InSAR technique. Subsidence in this area is not obvious. The abovementioned features indicate that sediment thickness may be an important cause of subsidence in Qianhai Bay. In addition to the geological conditions, we also analyzed the infrastructure construction of Qianhai Bay due to its important economic position. We collected three optical images corresponding to the acquisition time of SAR images and generated the built-up areas of Qianhai Bay with visual interpretation. The results are shown in Fig. 11. The range of the changing built-up area is consistent with the location of subsidence, which indicates that urban construction may be one of the factors affecting land subsidence in this area [45].

B. Relationship Between Land Subsidence and Subway
As described in Section III-C, obvious subsidence is observed along the subway. However, the relationship between land subsidence and subway distance is still unclear. To generate the relationship between these two factors, we adopted a proximity analysis to extract the minimum distance of each point from the subway. First, we set a distance threshold (500 m in this case) to select valid coherent points. Next, the distance from the subway was set as the independent variable, and the subsidence rate was set as the dependent variable. To better analyze the relationship, we divided the subsidence rate into four classes according to the natural breakpoint method. Then, the correlation between subsidence and subway distance was calculated based on Spearman's rank correlation coefficient method [33]. The results are shown in Table II. A negative correlation between the subsidence rate and distance from the subway is found when the subsidence rate is between 20.5 and 49.5 mm/a. The spearman's rank correlation coefficient is −0.1967. This indicates that regions closer to the subway may be more prone to land subsidence. But this correlation is quite weak [60], which means that the impact of distance from the subway on land subsidence is relatively small at present.

C. Relationship Between Land Subsidence and Road Network Density
In addition to the underground subways, Shenzhen also has abundant ground transportation. The mileage of the road network reached 7500 km in 2020, an increase of approximately 1000 km compared with 2015. To investigate the relationship between subsidence and the road network, we calculated the road network density of Shenzhen with a grid whose size was 1000 m × 1000 m. The road network density is defined as the road mileage of each grid divided by the grid area. It is noted that the road network density in this section represents the progress of urban exploitation. For example, higher road network density means less engineering construction. By contrast, lower road density means areas were under construction, which can also lead to land subsidence. Then, we divided the road network density grids into five levels according to the natural breakpoint method. Areas with higher road network density are mainly concentrated in southern and northeastern Shenzhen, as shown in Fig. 12. Similarly, we classified the subsidence rate into five classes and counted the corresponding coherent point percentage of each road network density level, as shown in Fig. 13(a). We noticed that the coherent point percentage has a negative  correlation with the subsidence rate when the road network density is higher (L4, L5). In contrast, a positive correlation between the coherent point percentage and subsidence rate is observed when the road network density is low (L1, L2, L3). Additionally, we also plot the scatter map between the road network density and the coherence point percentage for five subsidence rate classes, as shown in Fig. 13(b)-(f). We note that there is a negative correlation between the road network density level and coherent point percentage when the subsidence rate is over 5mm/a, as shown in Fig. 13(c)-(f). The above study indicates that areas with low road network density may be more likely to occur obvious subsidence, which can be caused by abundant engineering construction.

D. Evolution Feature of Land Subsidence in the Reclamation Area
Land reclamation is a typical symbol in Shenzhen's urbanization process, which can lead to land subsidence in different  degrees. To analyze the evolution feature of land subsidence in reclamation areas, we obtained the changes in reclamation areas from 1979 to 2015 with the assistance of optical images derived from Google Earth. We first divided the reclamation time into six phases, as shown in Fig. 14(a). It is noted that the coastline is stable and there is no significant new reclamation area since 2015. Then, we calculated the percentage of coherent points with subsidence rates over 5, 10, and 20 mm/a in the partitioned reclamation areas. The results are shown in Fig. 14(b). Two different subsidence stages were found for the three subsidence rate intervals. The coherent point percentage of each subsidence rate interval became decreased with the accumulation of reclamation time from 1979 to 2001, which is consistent with the settlement law of landfill or reclamation areas: the subsidence rate slows down with time [61]. However, subsidence became active with the reclamation time since 2001, indicating that the subsidence in the reclamation area is probably in an increasing trend. It will take some time for the subsidence to stabilize. Moreover, frequent human activities such as engineering construction can also accelerate land subsidence in the reclamation areas.
According to the previous analysis, land subsidence in the reclamation area of Shenzhen is mainly caused by the compaction of landfill material, underground exploitation, and ground infrastructure. Hence, three corresponding factors, distance from the coastline, distance from the subway, and road density, were selected to analyze the dominant factor that leads to land subsidence in reclamation areas. We used the Pearson correlation coefficient to calculate the correlation coefficient with a level of significance smaller than 0.1, as shown in Table Ⅲ. The correlation coefficients are −0.101, −0.084, and −0.122 for the distance from the coastline, the distance from the subway, and road density, respectively. Among these three factors, the road density and the distance from the coastline showed a modest correlation with land subsidence.

E. Relationship Between Subsidence and Urban Construction in Guangming District
Guangming District is one of the most noticeable areas with urban expansion and urban construction in Shenzhen over the past few years. To analyze the relationship between land subsidence and urban construction, we adopted building change as a proxy for urban construction events in the Guangming district. We generated the newly added buildings between 2015 and 2021 with visual interpretation based on historical optical images, as shown in Fig. 15. The changing buildings are mainly located around the edge of the Guangming district, which is consistent with the trend of urban expansion. Additionally, land subsidence areas with obvious rates are mainly located within the changing  areas. S1-S6 in Fig. 15 show the land subsidence within six typical regions of building change in the Guangming district. Fig. 16 shows a comparison of historical images of these regions during the six years. This indicates that urban construction may affect the occurrence of land subsidence.
In addition to the qualitative analysis, we also calculated the percentage of changing building areas. Specifically, we divided the subsidence rate into five grades and calculated the percentage of the changing building pixels of each grade. The histogram between these two factors is shown in Fig. 17. It is clear that the larger the subsidence rate grade is, the higher the percentage. The percentage for each grade is 5.57%, 16.4%, 21.88%, 33.59%, and 43.48%. Both qualitative and quantitative analyses prove that land subsidence in Guangming District has a high correlation with urban construction during the urbanization process.

F. Comparison With Other Studies
Many scholars have generated the deformation of Shenzhen in different periods based on SAR data [41], [43], [44], [45], as shown in Table IV. In these studies, subsidence areas were  10, respectively. The causes of land subsidence mainly focused on land reclamation, underground activities, and building construction. Besides the spatial and temporal characteristics of subsidence, Xu et al. predicted that the subsidence along the coastal areas will continue after 2010 [41]. He et al. found that the safety of a building is affected by the construction phase difference and underground activities [45]. Liu et al. found that the seasonality of time series deformation in Shenzhen is mainly related to precipitation [44].
In this study, we generated the deformation of the entire Shenzhen except for a small area in the eastern. Areas with obvious subsidence are not only distributed on the coastal reclamation areas and along the subway, but also in the Shajing community and Yantian port with a maximum subsidence rate of 42 mm/a and 22.1 mm/a, respectively. The maximum subsidence rate of Baoan International airport, Qianhai bay, and Houhai reach 66.2 mm/a, 38.2 mm/a, and 44.3 mm/a respectively, which is larger than the value (25 mm/a) before 2010 [41]. This shows that the deformation of the three coastal areas is increasing from 2010 to 2021. In addition, we analyzed the relationship between land subsidence and typical urbanization process, such as land reclamation, urban construction, underground exploitation (subway), and road network density. The results show that land subsidence has a certain relationship with urban construction and road network density. Therefore, besides the reclamation areas, more attention should be paid to areas with intensive urban construction, low road density, and subway across in the future urbanization process, aiming to provide early warning for possible geological disasters.

V. CONCLUSION
In this study, we collected 153 Sentinel-1A images from 2015 to 2021 and adopted the IPTA-SBAS method to generate the time series deformation of Shenzhen. We qualitatively and quantitatively analyzed the spatial and temporal characteristics of land subsidence in Shenzhen. Land subsidence is mainly distributed along the coastal areas, especially in the Shajing community, Baoan International Airport, Qianhai Bay, Houhai, and Yantian Port. The maximum subsidence rate is approximately 74.9 mm/a. In addition to the coastal areas, obvious subsidence has also occurred in inland districts with abundant construction, such as the Guangming district. Land subsidence is mainly concentrated in the urban construction area, industrial area, farmed area, and residential area. Additionally, we also analyzed the subsidence along the 11 subways of Shenzhen. Among these metro lines, obvious subsidence has occurred along Line 11 and Line 5 with maximum subsidence rates of 71.5 mm/a and 53.4 mm/a, respectively.
Moreover, the relationship between land subsidence and typical urbanization progress was also discussed in the fields of land reclamation, transportation network including subway construction and road network density, and urban construction: 1) Land subsidence in the reclamation area is mainly caused by soft soil geological conditions and urban construction; 2) The correlation between the subsidence rate and the distance from the subway is weak, with a correlation coefficient of −0.197 when the subsidence rate was over 20 mm/a; 3) Land subsidence is more likely to occur in areas with low road network density; 4) Taking the Guangming district as an example, a positive correlation between land subsidence and urban construction is found, and obvious subsidence mainly occurs in areas with urban construction.
Based on the abovementioned conclusions, in addition to coastal areas of Shenzhen, more attention should be paid to inland areas with abundant urban constructions, low road density, and subway across for geological disaster monitoring in the future. The findings of this study can help guide urban planning and provide suggestions for government decision-making of similar cities.