sustainability

Article
The Monitoring and Analysis of Land Subsidence in Kunming
(China) Supported by Time Series InSAR
Bo Xiao 1,2 , Junsan Zhao 1 , Dongsheng Li 3, *, Zhenfeng Zhao 1,2 , Wenfei Xi 4 and Dingyi Zhou 1

1 Faculty of Land Resources Engineering, Kunming University of Science and Technology,

Kunming 650093, China
2 Faculty of Road and Construction Engineering, Yunnan Communications Vocational and Technical College,
Kunming 650500, China
3 International Cooperation Department, Kunming Metallurgy College, Kunming 650033, China
4 Faculty of Geography, Yunnan Normal University, Kunming 650500, China
* Correspondence: 311705000620@home.hpu.edu.cn; Tel.: +86-183-8719-0341

Abstract: As urban construction has been leaping forward recently, large-scale land subsidence has
been caused in Kunming due to the special hydrogeological conditions of the city; the subsidence
scope has stretched out, and the subsidence rate has been rising year by year. As a consequence,
Kunming’s sustainable development has seriously hindered. The PS-InSAR (Persistent Scatterer In-
terferometric Synthetic Aperture Radar) and the SBAS-InSAR (Small Baseline Subsets Interferometric
Synthetic Aperture Radar) technologies were adopted to process the descending Sentinel-1A data
stacks from July 2018 to November 2020 to monitor the land subsidence of Kunming, so as to ensure
the sustainable development of the city. Moreover, the causes were analyzed. As revealed by the
results, (1) the overall subsidence trend of Kunming was large in the south (Dian lakeside), whereas
it was relatively small in the north. The significant subsidence areas showed major distributions
in Xishan, Guandu and Jining district. The maximal average subsidence rates of PS-InSAR and
Citation: Xiao, B.; Zhao, J.; Li, D.;
SBAS-InSAR were −78 mm/a and −88 mm/a, respectively. (2) The ground Subsidence field of
Zhao, Z.; Xi, W.; Zhou, D. The
Monitoring and Analysis of Land
Kunming was analyzed, and the correlation coefficient R2 of the two methods was reported as 0.997.
Subsidence in Kunming (China) In comparison with the leveling data of the identical period, the root mean square error (RMSE) is
Supported by Time Series InSAR. 6.5 mm/a and 8.5 mm/a, respectively. (3) Based on the urban subway construction data, geological
Sustainability 2022, 14, 12387. structure, groundwater extraction data and precipitation, the causes of subsidence were examined. As
https://doi.org/10.3390/ revealed by the results, under considerable urban subways construction, special geological structures
su141912387 and excessive groundwater extraction, the consolidation and compression of the ground surface
Academic Editor: George D.
could cause the regional large-area subsidence. Accordingly, the monthly average precipitation in
Bathrellos Kunming in the identical period was collected for time series analysis, thereby indicating that the land
subsidence showed obvious seasonal variations with the precipitation. The results of this study can
Received: 27 July 2022
provide data support and facilitate the decision-making for land subsidence assessment, forecasting
Accepted: 26 September 2022
and construction planning in Kunming.
Published: 29 September 2022

Publisher’s Note: MDPI stays neutral Keywords: Sentinel-1A; time series InSAR; land subsidence; attribution analysis
with regard to jurisdictional claims in
published maps and institutional affil-
iations.

1. Introduction
Land subsidence is a local subsidence movement due to consolidation and compres-
Copyright: © 2022 by the authors.
sion of underground loose strata on the basis of natural factors and human factors, thus
Licensee MDPI, Basel, Switzerland. resulting in the reduced ground elevation. It is a loss of environment and resources that
This article is an open access article is difficult to recover [1–3]. As urbanization has rapidly expanded, land subsidence has
distributed under the terms and turned out to be a geological disaster hindering urban development and agglomerations [4].
conditions of the Creative Commons Since the introduction of the 13th Five-Year Plan, Yunnan Province has expedited the con-
Attribution (CC BY) license (https:// struction journey of building a regional international center city based in southwest China
creativecommons.org/licenses/by/ that will face the whole nation and impact South and Southeast Asia as part of the active
4.0/). rollout of the national “Belt and Road” development strategy. Land subsidence turned out

Sustainability 2022, 14, 12387. https://doi.org/10.3390/su141912387 https://www.mdpi.com/journal/sustainability

Sustainability 2022, 14, 12387 2 of 21

to be an urgent issue in urban construction. In particular, the main Kunming urban area,
i.e., the core area of construction, is considered the most serious area of land subsidence.
Accordingly, the effective monitoring of land subsidence in the Kunming city main urban
area and its surrounding areas, as well as the spatial-temporal variation causes of land
subsidence, can be referenced for the subsequent urban construction and scientific land
subsidence prevention and control.
Early studies initiated a few years ago using leveling, total station triangulation, and
Global Navigation Satellite System (GNSS) monitoring as well as discrete point-based
geodesy methods [5–7]. The mentioned conventional geodetic surveys can be used for high-
precision subsidence monitoring. However, these surveys have several limitations (e.g.,
high cost, small monitoring scope, and rigorous requirements for monitoring environment),
so large-area monitoring results are difficult to obtain efficiently. In contrast with conven-
tional deformation monitoring methods, Differential Interferometric Synthetic Aperture
Radar (DInSAR) [8–10] developed rapidly in the past 30 years. It is widely used in urban
surface subsidence monitoring [11–13], landslide disaster identification [14,15], earthquake
deformation monitoring [16–18] and other geological disaster fields [19] because of its large
space coverage, short revisit period, high spatial resolution, high deformation monitoring
accuracy, all-day, all-weather, no personnel to reach the disaster or subsidence area for
direct measurement and so on [20,21]. This technology has gradually become one of the
indispensable technical means of ground subsidence monitoring. Over the past few years,
with several radar satellites launching and data sharing, the monitoring accuracy of InSAR
technology has continuously increased. Moreover, data processing methods have been
optimized, and coherent target-based temporal InSAR technologies have been proposed
successively, PS-InSAR technology [22,23], SBAS-InSAR technology [24], etc. Temporal
InSAR technology, by exploiting SAR data stacks for temporal analysis and processing, can
effectively overcome the conventional D-InSAR in space-time incoherence errors, atmo-
spheric delay, terrain effect and other aspects, and realize inverse surface deformation in
high precision.
Land subsidence has been detected and monitored over the past few years using
time-series InSAR [25–29]. Orhan [30] investigated the land subsidence in the Konya city
center and its surroundings based on SAR and global navigation satellite system (GNSS)
data. In addition, the correlation between the land subsidence observed in the region was
examined with land use/land cover changes and groundwater data. Solano-Rojas [31]
presented a novel method for extracting signals relating to geotechnical and infrastructure
monitoring at individual-building scale in fast-subsiding metropolises through band-pass
filtering. Moreover, the usefulness of their method was verified by applying the method
to Mexico City. Orhan et al. [32] expanded the InSAR investigation of land subsidence
in Karapınar region using a variety of data; GNSS observation data, monthly groundwa-
ter level data and InSAR data were adopted to gain more insights into the correlations
between groundwater extraction, land subsidence and sinkhole formation in the Kara-
pinar basin. The above research obtained a surface velocity field based on time-series
InSAR. The representative studies on land subsidence monitoring in Kunming’s main
urban area by adopting time-series InSAR technology are presented below. Zhu et al. [33]
utilized a multi-temporal InSAR technique to obtain a decade of times series deformation
of Kunming. Shao et al. [34] employed PS-InSAR technique to monitor the land subsidence
in Kunming. As revealed by the results, the land subsidence in Kunming was mainly
concentrated along the Lake Dian and the southern urban area, and the trend in land
subsidence was slowing down. Xiong et al. [35] used SBAS to process 28 sentinel-1A dual
polarization images explore the land subsidence in Kunming. As revealed by the results,
there was a relatively significant subsidence phenomenon in the Kunming main urban
area, and significant subsidence funnel was seriously formed in some areas. The maximal
annual average sedimentation rate is −54 mm/a (VV polarization) and −54 mm/a (VH
polarization). Jiang et al. [36] employed PS-InSAR and SBAS-InSAR technique to study
the surface subsidence of the Kunming main urban area respectively. The results showed
Sustainability 2022, 14, 12387 3 of 21

that the correlation between the two methods is high, and the maximal annual subsidence
rate in 2018–2019 reached −39 mm/a. Time series InSAR technique can conduct high
precision land subsidence monitoring in Kunming, whereas the study also showed that the
researchers only use temporal InSAR technique to detect subsidence funnel in Kunming
urban area. The causes of the spatial-temporal variation in land subsidence in Kunming
were comprehensively explored without the inversion of regional subsidence rate, urban
subway construction, geological structure, rainfall data and field survey data.
Due to current situation of large-scale urban construction in Kunming city, this study
used PS-InSAR and SBAS-InSAR technologies through 69 C-band Sentinel-1A data stacks
from July 2018 to November 2020 to study the characteristics of large-scale land surface
subsidence in Kunming main urban area, which obtained subsidence scope and accumu-
lated subsidence amount of Kunming main urban areas. It revealed spatial distribution
characteristics of ground deformation in Kunming. Then, the accuracy of InSAR results is
verified from three aspects of internal precision, external precision and field verification.
Lastly, the causes of surface subsidence in Kunming main urban area were explored com-
prehensively from the perspectives of urban subway construction, fault tectonic activity,
groundwater level change and rainfall. A subsidence monitoring and analysis system was
proposed. On the one hand, this study overcame the conventional surface subsidence
monitoring technology disadvantages and improved the spatial and temporal resolution,
which reduced the economic cost significantly. It could realize small deformation accurate
monitoring in a wide range of urban areas and provide strong technical support for moni-
toring surface subsidence in urban areas. On the other hand, this study analyzed the causes
of Kunming city ground subsidence from different aspects. A subsidence monitoring and
analysis system was established in this study to solve the problems of subsidence reason is
not specific and the subsidence trend is not clear in urban areas, which is beneficial to the
inversion of the subsidence trend in urban areas.
This study is organized as follows: basic information in the study area and the datasets
are given in Section 2. The methodologies are elaborated in Section 3. Subsequently, the re-
sult and analysis are elucidated in Section 4, followed by the discussions in Sections 5 and 6
in the end.

2. Study Area and Datasets

2.1. Study Area
Kunming (102◦ 330 1300 E~102◦ 590 1600 E; 24◦ 440 1200 N~25◦ 110 2500 N) is situated in the
central Yunnan-Guizhou Plateau, bordering Lake Dian in the south, as well as surrounded
by mountains on the other three sides. Kunming acts as the frontier and gateway of China
to Southeast Asia, South Asia and even the Middle East, Southern Europe and Africa,
exhibiting unique geographical advantages. It connects the coastal areas of Guizhou and
Guangxi in the east, facing Sichuan and Chongqing in the north, Thailand and Cambodia
in the south, Myanmar, India and Pakistan in the west. The altitude of Kunming is approxi-
mately 1891 m. The overall terrain is high in the north and low in the south, progressively
descending form north to south. Besides, the region has a subtropical-plateau monsoonal
climate. As impacted by the warm and humid airflow in the southwest of the Indian
Ocean, the sunshine is long, the frost period is short, and the annual average temperature
is 15 ◦ C. The annual precipitation is 1450 mm, and the precipitation is concentrated in May
to October, taking up nearly 85% of the annual precipitation. From November to the next
April is the dry season, taking up only about 15% of the annual precipitation.
The Kunming main urban area is located in what was a fault basin of the Late Ceno-
zoic, where loose sedimentary layers are broadly developed. On the whole, it comprises
lacustrine silt and clay with multi-layer silt, peat and lignite layers [37]. Given the geologi-
cal boreholes, the quaternary sedimentary phase transition is complicated vertically and
horizontally. Delta deposits are dominant in the northeast and lacustrine deposits in the
southwest. The soil structure stratum comprised alternating soft and hard layers, interbed-
ded aquifer and aquifer, interlaced high compressibility soil and low compressibility soil.
 geological boreholes, the quaternary sedimentary phase transition is complicated verti-
cally and horizontally. Delta deposits are dominant in the northeast and lacustrine depos-
Sustainability 2022, 14, 12387
its in the southwest. The soil structure stratum comprised alternating soft and hard layers,
4 of 21
interbedded aquifer and aquifer, interlaced high compressibility soil and low compressi-
bility soil.
As revealed from the hydrogeological conditions of Kunming main urban area, the
KunmingAs revealed from the
soil stratum hydrogeological
is characterized conditions
by soft of Kunming
structure, main urban
low consolidation area,and
degree the
Kunming soil stratum is characterized by soft structure, low consolidation degree and
high compressibility. Considerable urban construction in this area can easily cause foun- high
compressibility. Considerable urban construction in this area can easily cause foundation
dation vertical displacement variations, thereby posing the risk of large-scale engineering
vertical displacement variations, thereby posing the risk of large-scale engineering con-
construction collapse. Thus, this study selected Chenggong, Guandu, Xishan, Panlong,
struction collapse. Thus, this study selected Chenggong, Guandu, Xishan, Panlong, and
and Wuhua Districts and the northern part of Jinning District as the study area (Figure 1).
Wuhua Districts and the northern part of Jinning District as the study area (Figure 1).

Figure 1.
Figure Schematic diagram
1. Schematic diagram of
of the
the study
study area.
area.

2.2. Datasets
2.2. Datasets
In this study, 69 C-band Sentinel-1A images downloaded by ESA acted as research
In this study, 69 C-band Sentinel-1A images downloaded by ESA acted as research
data. The imaging lasted from July 2018 to November 2020, under the polarization mode
data. The imaging lasted from July 2018 to November 2020, under the polarization mode
of VV polarization, the wavelength of 5.63 cm, as well as the revisit period of 12 days.
of VV polarization,
Furthermore, DEM the
datawavelength
with 30 mof 5.63 cm, as
resolution well as by
provided theJapan
revisitSpace
period of 12 days.
Agency were
Furthermore, DEM data with 30 m resolution provided by Japan Space Agency were se-
selected for eliminating the effect of topographic relief. Table 1 lists the basic parameters of
lected fordata.
satellite eliminating the effect of topographic relief. Table 1 lists the basic parameters of
satellite data.
Table 1. Basic parameters of satellite data.
Table 1. Basic parameters of satellite data.
Dataset Orbit Mode Incidence Wavelength Polarization Mode Number Time Span
Dataset Orbit Mode Incidence
◦ Wavelength Polarization Mode Number Time Span
Sentinel-1A Descending 39.5 C VV 69 20 July 2018~18 November 2020
20 July 2018~18 November
Sentinel-1A Descending 39.5° C VV 69
2020
3. Methodology
3.1. PS-InSAR Method
3. Methodology
The basic principle of PS-InSAR technology [23] is to adopt SAR datasets (generally
3.1. PS-InSAR Method
over 20) covering the identical area to detect stable point targets not affected by temporal
The basic
and spatial principle
baseline of on
based PS-InSAR technology
a statistical analysis [23] is amplitude
of the to adopt SARanddatasets (generally
phase information
over 20) covering
stability the identical
in time series. area to detect
The mentioned stablestable
pointpoint targets
targets not affected
are likely by temporal
to include artificial
and
buildings, exposed rocks, artificial angle reflectors and others, which are termedinformation
spatial baseline based on a statistical analysis of the amplitude and phase PS points as
stability
almost notin affected
time series. The mentioned
by speckle stable point
noise in time-series SARtargets
imagesareand
likely
stilltomaintaining
include artificial
stable
buildings, exposed rocks, artificial angle reflectors and others, which are
scattering after a long time interval. The geometric size of PS points is significantlytermed PSsmaller
points
than the SAR datasets’ spatial resolution unit, whereas their backscattering coefficients
dominate the echo signal in the whole resolution unit, which are independent of the spatial
baseline and show good coherence [38]. According to time-series datasets, even if no
interference fringes can be seen on a single interferogram, the whole interference phase can
be obtained from PS points phase information by interpolation or fitting.
Sustainability 2022, 14, 12387 5 of 21

The PS-InSAR technology primarily focuses on the spatially dispersed PS points, and
the image in the time-series dataset with N single-look complex SAR images is expressed
as Si (i = 1, 2, . . . , N ). The interferogram comprises the i and k SAR images, expressed
I i,k = si s∗k . * represents complex conjugate. The target p0 is assumed as the reference
point, so the interference phase ∆∅i,k i,k
p,p0 = ∠ I p,p0 will be determined by noise differences
between the target p and the reference point p0 , atmospheric disturbance, deformation, and
elevation. Phase difference attributed to elevation and linear deformation difference of the
two targets is denoted in the following:
4π 1
∆∅i,k
H,p,p0 = ∆h p,p0 Bni,k (1)
λ Rsinθ
and
4π
∆∅i,k
D,p,p0 = ∆v p,p0 Bti,k (2)
λ
In Equations (1) and (2), ∆h p,p0 and v p,p0 are the difference between target p and p0 in
elevation and deformation rate, respectively. Bni,k and Bti,k represent the space baseline and
time baseline between SAR image i and k, respectively, perpendicular to the radar sight
line. θ represents the SAR incident angle; R denotes the oblique distance between the target
and the satellite; λ represses the radar signal wavelength.
In the classical PS-InSAR technique, the elevation and subsidence rate of the target are
estimated by maximizing the temporal coherence of ξ p [39]:

(∆ĥ p , ∆v̂ p ) = arg max( ξ p )


(3)

1 j(∆∅i,k i,k i,k

p − ∆∅ H,p − ∆∅ D,p )
ξp =
M ∑e (4)
i,k

In Equations (3) and (4), M is the number of interferogram; ∆∅i,k

p is the obtained differ-
ential interference phase (the flat part and rough terrain phase have been removed); ∆∅i,k
H,p
is the elevation-related part phase of Equation (1); ∆∅i,kD,p is the topographic deformation
phase related part of Equation (2).
Surface deformation ∆∅i,kD,p consists of linear and nonlinear deformation phases:

∆∅i,k i,k i,k

D,p = ∆∅linear + ∆∅non−linear
i,k
(5)
λ ( ∆v p Bt
4π
= + De f onon−linear )

In Equation (5), ∅i,k non−linear is a nonlinear form variable. Equation (4), the maximum
factor is ξ̂ p = max ( ξ p ) known as the temporal coherence, the time coherence of target can
−δ2 /2
approximate expression for phase deviation factor of residual function ξ̂ p = e φ . Lastly,
the variance of elevation and deformation rate estimated with the use of PS-InSAR, i.e., the
accuracy index δ∆h2 and δ2 with phase deviation δ2 , the difference between interferogram
∆v φ
geometry and time baseline (δφ2 and δB2 t ), as well as the number of interferometers M satisfy
the following [40]:
2
2 λRsinθ 2 δφ
δ∆h ≈( ) (6)
4π MδB2 n
and
2
2 λ 2 δφ
δ∆v ≈( ) (7)
4π MδB2
t

As expressed in Equations (6) and (7), the accuracy of the elevation and deforma-
tion rates obtained from PS-InSAR is considered a function of temporal coherence and
spatio-temporal baseline distribution. Such a model is not correlated with interference
combination, so its expression applies to all time-series InSAR image analysis techniques.
Sustainability 2022, 14, 12387 6 of 21

3.2. SBAS-InSAR Method

SBAS-InSAR refers to time-series InSAR analysis developed by Berardino, Lanari et al.
in accordance with a variety of strategies of PS-InSAR [24,41]. First, the time-series radar in-
terference datasets were integrated to generate several short-space baseline interferograms
that are capable of overcoming spatial de-correlation [24]. Subsequently, the chosen GCP
points are applied in terms of denoising and orbit refining. Lastly, the deformation rate was
obtained with the use of the singular value decomposition (SVD) method, and the SAR
datasets separated using large spatial baselines were connected to up-regulate the observed
data’s time sampling rate. SAR data processing is elucidated below:
(1) It obtained N + 1 SAR images in the identical area arranged in chronological order
t0 , . . . , tn , selecting one image as the master image, and registered other SAR images
on the master image. N + 1 SAR image generates M differential interferogram, and
N ( N +1)
M satisfies N2+1 ≤ M ≤ 2 .
(2) In terms of the j differential interferogram obtained according to SAR image obtained
from slave image t A and mater image t B (t B > t A ), the pixel interference phase with
azimuth direction coordinate x and range direction coordinate r is written as follows:

4π j j j
δφj ( x, r ) = φB ( x, r ) − φ A ( x, r ) ≈ [d(t B , x, r ) − d(t A , x, r )] + ∆φtopo ( x, r ) + ∆φ APS (t B , t A , x, r ) + ∆φnoise ( x, r ) (8)
λ
In Equation (8), j ∈ (1, . . . M); λ represents the signal central wavelength; d(t B , x, r )
and d(t A , x, r ) denote the cumulative shape variables of the phase at t B and t A for the radar
j
sight direction line d(t0 , x, r ) = 0. ∆φtopo ( x, r ) denotes the participating phase within the
differential interferogram. The DEM added to differential interference is assumed to have
high accuracy and can largely remove the topographic phase, so the differential interfero-
j
gram has small residual topographic information, which is negligible. ∆φ APS (t B , t A , x, r )
j
denotes the atmospheric delay phase; ∆φnoise ( x, r ) represents coherent noise. Without
considering noise phase, residual terrain phase, or atmospheric delay phase, Equation (8)
is simply expressed in the following:

δφj ( x, r ) = φB ( x, r ) − φ A ( x, r )
(9)
≈ 4πλ [ d ( t B , x, r ) − d ( t A , x, r )]

(3) For generating the subsidence sequence that achieves physical significance, the prod-
uct of the average phase velocity and time of two acquisition times expresses the
phase in Equation (9):
φ j − φ j −1
vj = (10)
t j − t j −1
The j interferogram’s phase value is:
t B,j

∑ (tk − tk−1 )vk = δφj (11)

k =t A,j+1

The velocity integral in the respective period on the time interval of the main and
auxiliary images is written as the matrix form:

Bv = δφ (12)

Equation (12) represents an M × M matrix. Matrix B may have rank deficiency since
the small baseline sets employ multi-master image strategy. The generalized inverse matrix
of matrix B is established with SVD method, as an attempt to yield the minimal velocity
vector norm solution. Lastly, the shape variable of the respective period is acquired by
velocity integration in the respective period.
Sustainability 2022, 14, 12387 7 of 21

3.3. Data Processing

This study was primarily based on SARscape software and 69 Sentinel-1A images.
The specific experimental processing included SBAS interference processing and PS inter-
ference processing.

3.3.1. SBAS Interference Processing

(1) Data preprocessing: reading 69 original images and clipping the study area using
vector files.
(2) Connection diagram and connection pair image generation: the super master image
was the image dated 5 November 2018. The max normal baseline was set to 2%, and
the max temporal baseline threshold was set to 60 days. Lastly, 316 image pairs were
generated, as depicted in Figure 2a.
(3) Interferometric workflow: 316 image pairs were conducted interferometric proces-
sion. Coregistration, interferogram generation, interferogram flattening, filtering and
phase unwrapping were carried out. Coregistration all image pairs was conducted
with the super master image to prepare for orbit refinement, re-flattening and SBAS
inversion. The unwrapping method selected in this study was minimum cost flow.
The advantage of this method was that better results could be obtained when the
unwrapping was difficult due to large areas of low coherence or other factors limiting
growth. Goldstein was selected as the filtering method.
(4) Orbit refinement and re-flattening: This was adopted for estimating and removing
the residual phase. Ground control points (GCPs) were employed for correcting SAR
data because of inaccurate DEM geographic positioning or inaccurate satellite orbit.
Criteria for GCP selection are elucidated below: (1) Selecting the area without residual
terrain fringe; (2) Selecting a region without deformation, unless the deformation rate
of that point is known; (3) Selecting the point without phase jump; (4) The number of
selected points should be more than 20. A total of 50 GCPs were selected in this study.
(5) First inversion of SBAS: this step serves as the core of inversion, estimating deforma-
tion rate and residual terrain phase, and optimizing the input interferograms through
secondary unwrapping.
(6) The second inversion of SBAS: The residual phase of the SBAS second unwrapping
interferograms was separated into low-pass phase and high-pass phase. The low-pass
filter and high-pass filter are employed to filter the residual phase to obtain nonlinear
Sustainability 2022, 14, x FOR PEER REVIEW 8 of 22
deformation. The linear deformation phase and non- linear deformation phase are
added together to obtain accurate surface deformation information.
(7) Geocoding: the results of surface deformation were projected onto the geographic co-
A total ofsystem.
ordinate 2,535,786 stable
In this points,
study, the namely SDFPprojected
results were points, were
onto obtained in direction.
the vertical this experi-A
ment.
total of 2,535,786 stable points, namely SDFP points, were obtained in this experiment.

(a) (b)
Figure 2.2. Spatiotemporal
Figure spatiotemporal baseline
baselinedistribution.
distribution.(a,b)
(a) and
are (b) are SBAS-InSAR
SBAS-InSAR and PS-InSAR
and PS-InSAR spatio-
spatiotemporal
temporal baseline distribution, respectively.
baseline distribution, respectively.

3.3.2. PS Interference Processing

(1) Data preprocessing: Data preprocessing was the same as SBAS processing.
(2) Connection diagram generation: Likewise, the super master image was the image on
5 November 2018, and the connection diagram was depicted in Figure 2b. The max
Sustainability 2022, 14, 12387 8 of 21

3.3.2. PS Interference Processing

(1) Data preprocessing: Data preprocessing was the same as SBAS processing.
(2) Connection diagram generation: Likewise, the super master image was the image on
5 November 2018, and the connection diagram was depicted in Figure 2b. The max
temporal baseline was 60 days and the min temporal baseline was 12 days.
(3) Interferometric workflow: 69 SAR images were processing by coregistration, flattening
and differential interferometric. 1,295,860 effective stable PS points were seleteded out
by amplitude deviation index. Most of these PS points were distributed in buildings,
Bridges, etc.; Lake and green areas were almost no PS points.
(4) PS first inversion: the residual topography and the displacement rate were obtained to
flatten the synthetic interferogram. The main method was based on the identification
of a certain number of Persistent scatterers, which should be stable (change less than
1 mm) and can be detected by SAR antenna.
(5) PS second inversion: in this step, the atmospheric phase was removed for obtaining
the final rate of deformation.
(6) Geocoding.

4. Results and Analyses

4.1. Time-Series Land Subsidence
Since the deformation rate in the Line of Sight (LOS) direction is the surface deforma-
tion rate projection in the radar sight line, any surface deformation can be expressed by the
deformation in the North-South (SN), East-West (EW) and vertical (U) components [42]. As
revealed by the geometric relationship between radar side-view imaging, LOS deformation
and surface deformation observed by InSAR, the vertical deformation contribution rate to
LOS direction exceeds 90% regardless of orbit ascending or descending [43]. Since InSAR
observation results are most sensitive to vertical surface deformation, the LOS deformation
rate field determined in this study was employed as the analysis data of Kunming deforma-
tion trend. As decipted in Figure 3, the Sentinel-1A radar data covering the study area were
processed using SBAS-InSAR and PS-InSAR. Subsequently, the results obtained by the
two InSAR technologies are underwent Kriging interpolation by ArcGIS. The continuous
spatial deformation rate diagram of LOS to the same point in Kunming main urban area
from July 2018 to November 2020 was obtained. In Figure 3, the deformation rate was
positive, thereby demonstrating that the surface displacement moved in the direction of the
satellite, and the surface was uplifted. Under the displacement, the ground moved away
from the direction of the satellite and showed subsidence.
As depicted in Figure 3, the land subsidence covered a large area with uneven spatial
distribution in the study area. In addition, the land subsidence centers in the study
area were mainly concentrated in the northeast, east, southeast of Lake Dian and central
Kunming. The areas with relatively serious subsidence refer to A, B, C, D and E. Area A
was located near Daguan Park in Xishan District, including Zongshuying Street and Ma
Street. The maximal annual subsidence rate monitored by PS and SBAS was −78 mm/a
and −88 mm/a, respectively. Area B was located in Guandu District along the eastern
coast of Lake Dian, which has developed into a continuous subsidence area stretched out
to the urban center, with Liujia Street, Guandu Street and Yiliu Street involved. In the
above area, the maximal annual subsidence rate was −59 mm/a and −67 mm/a. Area C
was located in Xiaobanqiao Street near Guangwei Overpass in Guandu District, and the
maximal annual subsidence rate was −31 mm/a and −27 mm/a, respectively. Area D
was located in Luojia Village, Tujia Village and Zhonghe Village, Dayu Street, Chenggong
District, and the maximal annual subsidence rate was −34 mm/a and −51 mm/a. The
difference between the maximal annual subsidence rate obtained by the two monitoring
methods in this area was −18 mm/a, probably because the subsidence area was largely
distributed in farmland. The PS points were relatively sparse. Area E was located in Xinjie
Street near No.8 Middle School of Jinning County, with the maximal annual subsidence
rate of −54 mm/a and −56 mm/a, respectively. In the study area, only Panlong District
 Sustainability 2022, 14, 12387 9 of 21

Sustainability 2022, 14, x FOR PEER REVIEW 9 of 22

was relatively stable with relatively small subsidence, and the average annual Subsidence
rate was lower than −10 mm/a.

Figure 3. 3.
Figure Deformation
Deformation rate ofthe
rate of thestudy
study area.
area. (a)areand
(a,b) (b) of
results are results and
PS-InSAR of PS-InSAR
SBAS-InSAand SBAS-InSA
processing
processing respectively.
respectively. ①~④ fault,
1 ~ 4 are Puduhe are Puduhe fault, Chenggong-Fumin
Chenggong-Fumin fault, Heilongtan-Guandu,
fault, Heilongtan-Guandu, Baiyi-Hechong
Baiyi-Hechong fault, respectively.
fault, respectively.

4.2.
AsInSAR Accuracy
depicted Assessment
in Figure 3, the land subsidence covered a large area with uneven spatial
This study
distribution in theverified
study the reliability
area. of InSAR
In addition, thesubsidence monitoring
land subsidence results
centers in from three area
the study
aspects.
were mainly First, internal coincidence
concentrated accuracyeast,
in the northeast, [44] southeast
was used to of verify the consistency
Lake Dian and central of Kun-
experimental PS and SBAS results. Second, the reliability of InSAR was verified and
ming. The areas with relatively serious subsidence refer to A, B, C, D and E. Area A was
compared in accordance with the leveling results. Third, field sampling survey verified the
located nearofDaguan
accuracy monitoringPark in Xishan District, including Zongshuying Street and Ma Street.
results.
The maximal annual subsidence rate monitored by PS and SBAS was −78 mm/a and −88
(1) Internal coincidence accuracy verification analysis
mm/a, respectively. Area B was located in Guandu District along the eastern coast of Lake
In Figure
Dian, which has 3, the maximal
developed intoannual subsidence
a continuous rate mapsarea
subsidence obtained by PS
stretched and
out to SBAS
the urban
were, respectively, − 78 mm/a and − 88 mm/a, with a difference of 10
center, with Liujia Street, Guandu Street and Yiliu Street involved. In the above area, themm/a, thereby
demonstrating that serious surface subsidence exists in the study area. The location and
maximal annual subsidence rate was −59 mm/a and −67 mm/a. Area C was located in
range of the subsidence area identified by the two methods are consistent. To further
Xiaobanqiao Street near Guangwei Overpass in Guandu District, and the maximal annual
verify the reliability of the two methods based on the experimentally achieved results, the
subsidence rate wasmethod
cross-comparison −31 mm/a was and −27 to
adopted mm/a,
verifyrespectively.
the results. InArea D was located
that process, in Luojia
532 identical
Village,
pointsTujia
wereVillage
randomly and Zhonghe
extracted Village,
from Dayu Street,
the subsidence rate Chenggong
maps generated District, and the max-
by employing
imal annual subsidence rate was −34 mm/a and −51 mm/a. The
the PS-InSAR and the SBAS-InSAR respectively. The correlation coefficient diagram difference between
was the
maximal
drawn annual subsidence
by applying the averagerateannual
obtained by the rate
subsidence twoofmonitoring methods
PS as the horizontal in and
axis thiswith
area was
the average annual subsidence rate of SBAS as the vertical axis (Figure
−18 mm/a, probably because the subsidence area was largely distributed in farmland. The4). The correlation
R was 2 was 0.997, thereby demonstrating that the experimentally
0.998, Rsparse.
PS coefficient
points were relatively Area E was located in Xinjie Street near No.8 Middle
achieved results in this study have high correlation and consistency.
School of Jinning County, with the maximal annual subsidence rate of −54 mm/a and −56
mm/a, respectively. In the study area, only Panlong District was relatively stable with rel-
atively small subsidence, and the average annual Subsidence rate was lower than −10
mm/a.

4.2. InSAR Accuracy Assessment

This study verified the reliability of InSAR subsidence monitoring results from three
aspects. First, internal coincidence accuracy [44] was used to verify the consistency of ex-
perimental PS and SBAS results. Second, the reliability of InSAR was verified and com-
 were randomly extracted from the subsidence rate maps generated by employing the PS-
InSAR and the SBAS-InSAR respectively. The correlation coefficient diagram was drawn
by applying the average annual subsidence rate of PS as the horizontal axis and with the
average annual subsidence rate of SBAS as the vertical axis (Figure 4). The correlation
Sustainability 2022, 14, 12387 coefficient R was 0.998, R2 was 0.997, thereby demonstrating that the experimentally
10 of 21
achieved results in this study have high correlation and consistency.

Correlationcoefficient
Figure4.4.Correlation
Figure coefficientdiagram
diagramofofannual
annualaverage
averageSubsidence
Subsidencerate
rateofofPS
PSand
andSBAS.
SBAS.

Notably, the
Notably, the number
number ofof PS
PSpoints
pointsininthe
theresearch area
research obtained
area by by
obtained PS-InSAR technique
PS-InSAR tech-
nique were 1,295,860, and the number of SDFP points obtained by SBAS-InSAR wereIn
were 1,295,860, and the number of SDFP points obtained by SBAS-InSAR were 2,535,786.
this study, 532 corresponding points were selected randomly, which had the characteristics
2,535,786. In this study, 532 corresponding points were selected randomly, which had the
of small number and centralized distribution. Therefore, the calculated correlation coeffi-
characteristics of small number and centralized distribution. Therefore, the calculated cor-
cient was close to 1. It was not representative, which only indicates that the results obtained
relation coefficient was close to 1. It was not representative, which only indicates that the
by PS-InSAR and SBAS-InSAR were relatively consistent. In the experiment, it was found
results obtained by PS-InSAR and SBAS-InSAR were relatively consistent. In the experi-
that the correlation coefficient decreases nonlinearly with the increase of the corresponding
ment, it was found that the correlation coefficient decreases nonlinearly with the increase
points number involved in the calculation. We will be discussed in the future experiments.
of the corresponding points number involved in the calculation. We will be discussed in
(2) future
the Comparison and verification with leveling results
experiments.
(2) Comparison and verification
InSAR measurements withofleveling
show line results
sight deformation, whereas conventional leveling
measurements show vertical deformation. In order
InSAR measurements show line of sight deformation, whereas to verify theconventional
monitoring results
levelingof
measurements show vertical deformation. In order to verify the monitoring results ofwere
PS-InSAR and SBAS-InSAR, the monitoring data of 12 level points at the same time PS-
obtained,
InSAR andas depicted in Figure
SBAS-InSAR, 1. With the
the monitoring dataleveling pointpoints
of 12 level as theatcenter of thetime
the same circle, the
were
InSAR subsidence point in the 30 m radius buffer was extracted. After the
obtained, as depicted in Figure 1. With the leveling point as the center of the circle, the extraction
subsidence
InSAR was divided
subsidence point inbytheCOSθ
30 m (θ as the
radius incident
buffer was angle), theAfter
extracted. average
the subsidence
extraction sub-was
determined
sidence and compared
was divided by 𝐶𝑂𝑆𝜃 with(𝜃the
asleveling measurement
the incident result,
angle), the and subsidence
average the root mean square
was de-
error (RMSE) was introduced to assess the verification accuracy. The RMSE
termined and compared with the leveling measurement result, and the root mean square is determined
by Equation (13) [20]:
error (RMSE) was introduced to assess thesverification accuracy. The 𝑅𝑀𝑆𝐸 is deter-
n
mined by Equation (13) [20]: RMSE = ∑ ∆v2i /n (13)
i =1

where, ∆vi denotes the deformation rate difference between the deformation rate at the i-th
leveling point and InSAR monitoring results; n expresses the number of test sites.
Table 2 lists the comparison results of the leveling survey and InSAR average deforma-
tion rate. According to Table 2, the deformation rate root mean square error of the leveling
monitoring point and the PS and SBAS subsidence point are 8.5 mm/a and 6.7 mm/a
respectively, and the maximal deviation is less than 15 mm/a. The deviation ratio in
(10, 15) mm/a was nearly 10%. The above results suggest that PS and SBAS monitoring
results exhibit high reliability.

Table 2. The average subsidence rate comparison between leveling data and InSAR.

Number of Deviation in Deviation in Deviation in

RMSE Maximum Minimum
Stations (0 5) (%) (5 10) (%) (10 15) (%)
PS 8.5 13.0 1.4 14 66.3 23.5 10.2
SBAS 6.7 11.2 0.9 14 71.7 19.2 9.1
 Table 2. The average subsidence rate comparison between leveling data and InSAR.

Number of Deviation in Deviation in Deviation in

𝑹𝑴𝑺𝑬 Maximum Minimum
Stations (0 5) (%) (5 10) (%) (10 15) (%)
PS 8.5 13.0 1.4 14 66.3 23.5 10.2
SBAS 2022,
Sustainability 6.714, 12387 11.2 0.9 14 71.7 19.2 9.1 11 of 21

(3) Field investigation

(3) To more
Field specifically verify the location and scope of the identified subsidence area,
investigation
subsidence area B, close to the downtown of Kunming, was selected to conduct field in-
To more specifically verify the location and scope of the identified subsidence area,
vestigation for verifying the accuracy of the experimentally achieved results. Figure 5 is a
subsidence area B, close to the downtown of Kunming, was selected to conduct field
photo taken on
investigation for the spot. In
verifying thethe area connected
accuracy by subsidence
of the experimentally centersresults.
achieved in Figure 3, house
Figure 5 is a
wall cracking (Figure 5a), house tilt (Figure 5b) and road foundation deformation
photo taken on the spot. In the area connected by subsidence centers in Figure 3, house cracking
(Figure
wall 5c) were
cracking caused
(Figure in the tilt
5a), house field, and 5b)
(Figure the and
trend of foundation
road deterioration continued,cracking
deformation thereby
demonstrating that InSAR results exhibit high monitoring accuracy.
(Figure 5c) were caused in the field, and the trend of deterioration continued, thereby
demonstrating that InSAR results exhibit high monitoring accuracy.

Figure 5. Field photos. (a–c) A house wall crack, house tilt and road foundation deformation cracking
Figure by
caused 5. Field photos. (a–c)
the subsidence, A house wall crack, house tilt and road foundation deformation crack-
respectively.
ing caused by the subsidence, respectively.
5. Discussion
5. Discussion
5.1. Impact Analysis of Ground Subsidence and Urban Subway Construction
5.1. Impact
Subway Analysis
is knownof Ground Subsidence
as another landmark and of
Urban Subway
urban Construction
modernization. Kunming metro was
built Subway
in 2008. Table 3 lists
is known asthe construction
another landmark andofoperation time of each subway
urban modernization. Kunming line. In 2017,
metro was
all six lines running through the Kunming main urban area and gradually expanded
built in 2008. Table 3 lists the construction and operation time of each subway line. In 2017, to
Anning, Songming, Yiliang, Jining, Yangzong and other counties and urban
all six lines running through the Kunming main urban area and gradually expanded to areas. The rapid
subway
Anning,construction and the line
Songming, Yiliang, network
Jining, coverage
Yangzong andmay aggravate
other countiesthe
andground
urbansubsidence
areas. The
according
rapid subwayto subway line. and
construction Basedtheon SBAS-InSAR
line experimentally
network coverage achieved
may aggravate results,sub-
the ground the
monitoring results to
sidence according within
subway500line.
m along
Based the
onsubway line were
SBAS-InSAR extracted (Figure
experimentally 6), and
achieved the
results,
evolution characteristics of land subsidence along the subway network under construction
in Kunming were analyzed. As depicted in Figure 6, there were different degrees subsidence
along Line 1 and branch line. The maximal annual subsidence rate was −14 mm/a, which
was located at Xiaodong Village Station. The first phase of Metro Line 2 was opened
in 2014 (from the Northern coach Station to South Ring Road Station), with a maximal
annual subsidence rate of −14 mm/a, which was located between Sijiaying station and
the northern coach station. The second phase of Line 2 (South Ring Road station to
Baofeng Village Station) was fully constructed in 2015 and is expected to open in 2022. The
phase II project of Line 2 under construction was 12.7 km long and built along Lake Dian,
passing through continuous subsidence areas, among which the most serious subsidence
is the section from Baofeng Village station to Guangfu Road, with the maximal annual
subsidence rate of −73 mm/a. Line 3 represents a completed running through Kunming
city’s east-west direction, with 23.4 km length in total and relatively stable overall. There is
a slight subsidence between Liangjiahe Station and municipal gymnasium Station, with
the maximal annual subsidence of −12 mm/a. The maximal subsidence area of Line 3 is
located near the eastern coach station, and the maximal annual subsidence is −16 mm/a.
Metro Line 4 connects northwest and southeast of Kunming. There are four significant
subsidence areas along the line, from northwest to southeast, located at the starting stations
Jinchuan Road Station, Guangwei Station, Gucheng Station and Kunming South Railway
Station. The corresponding maximal annual subsidence rate is −24 mm/a, −26 mm/a,
−11 mm/a and −14 mm/a. Line 5 is under construction. As impacted by construction
and train cyclic load operation, uneven ground surface subsidence occurs as well. The
Sustainability 2022, 14, 12387 12 of 21

subsidence area of Line 5 is mainly located between Dianchi International Convention and
Exhibition Center and Fuhai in Liujia Street, with the maximal annual subsidence rate of
−67 mm/a. Metro Line 6 in the study area has been stable since it was put into operation.
Only one point of subsidence was detected in the eastern coach station, and the maximal
annual settlement was −16 mm/a.

Table 3. Timetable for construction and operation of Kunming subway lines.

Subway Line Time (20-) Stations Length/km

08 09 10 11 12 13 14 15 16 17 18 19 20 21 22
Line 1 Phase I 18 41.4
Line 1 Branch 5 5.3
Northwest Extension of Line 1 8 7.6
Line 2 Phase I 14 22.8
Line 2 Phase II 10 12.7
Line 3 20 23.4
Line 4 29 43.4
Line 5 22 25.9
Line 62022,
Sustainability Phase14,
I x FOR PEER REVIEW 4 26.6
13 of 22
Line 6 Phase II 5 7.3
Note: Yellow indicates the subway construction stage; green indicates the subway operation stage.

Figure 6. Deformation
Figure 6. Deformation rate
rate within
within 500
500 m
m buffer
buffer of
of Kunming
Kunmingmetro
metroline.
line.

5.2.
5.2. Analysis
Analysis of
of the
the Land
Land Subsidence
Subsidence Impact
Impact and
and Fault
Fault Tectonic
TectonicActivity
Activity
Kunming was located in late Cenozoic downfaulted basin with a complicated structure.
Kunming was located in late Cenozoic downfaulted basin with a complicated struc-
It covers four major fault zone surrounding distribution, i.e., the Kunming Baiyi-Hengchong
ture. It covers four major fault zone surrounding distribution, i.e., the Kunming Baiyi-
faults in eastern and central Heilongtan-Guandu faults, as well as west Pudu river and
Hengchong faults in eastern and central Heilongtan-Guandu faults, as well as west Pudu
north west of Chenggong-Fuming fault (Figure 3). The fault is the surrounding potential
river and north west of Chenggong-Fuming fault (Figure 3). The fault is the surrounding
earthquake source area in Kunming, as well as one of the factors of land subsidence
potential earthquake source area in Kunming, as well as one of the factors of land subsid-
that cannot be ignored. According to existing research, compared with linear urban
ence that cannot be ignored. According to existing research, compared with linear urban
land subsidence, subway and expressway, the crustal deformation in fault zone has its
land subsidence, subway and expressway, the crustal deformation in fault zone has its
particularity, thereby suggesting that the small deformation magnitude and the remarkably
particularity, thereby suggesting that the small deformation magnitude and the remarka-
slow change. The active intensity of most faults is several million meter per year, and the
bly slow change. The active intensity of most faults is several million meter per year, and
maximal exceeds a dozen million meter [43]. The Chenggong-Fumin fault with serious land
the maximal exceeds a dozen million meter [43]. The Chenggong-Fumin fault with serious
subsidence was selected as the research object, and the cross-sectional lines P1 P10 , P2 P20 , P3 P30
land subsidence
perpendicular was
to the selected as the fault
Chenggong-Fumin research
zone object, and (Figure
were made the cross-sectional lines
3b), and the cross-
𝑃 𝑃 , 𝑃 𝑃 , 𝑃 𝑃 perpendicular to the Chenggong-Fumin fault zone were made (Figure
sectional map perpendicular to the fault zone was created (Figure 7). The intersection of the
3b), and the cross-sectional map perpendicular to the fault zone was created (Figure 7).
The intersection of the blue vertical line and the section line in the figure represent the
position of the fault line and the subsidence rate on the cross-sectional map. As depicted
in Figure 7, the current vertical crustal deformation rate was relatively high in the near
Chenggong-Fumin fault area and the maximal annual subsidence rate was −70 mm/a,
Sustainability 2022, 14, 12387 13 of 21

blue vertical line and the section line in the figure represent the position of the fault line and
the subsidence rate on the cross-sectional map. As depicted in Figure 7, the current vertical
crustal deformation rate was relatively high in the near Chenggong-Fumin fault area and
the maximal annual subsidence rate was −70 mm/a, thereby demonstrating significant
tectonic movement in the study area is, especially controlled by the Chenggong-Fumin
fault. According to existing studies on land subsidence, seismicity and tectonic movement
in Kunming city [26,45], Since the late Cenozoic, the stress field in the Kunming urban area
was dominated by the northward compressive stress, and the tensile stress was generated
in the basin area bounded by the Puduhe fault. On that basis, the five blocks of Xishan,
Puji, Sheshan, Longtan and Baiyi, are controlled and pulled northward successively by the
Heilongtan–Guandu fault and Baiyi Hengchong fault. As impacted by the different degree
in the consolidation fault zone, the distance of fault block pulling out was not identical.
The longest block was Heilongtan, and the nearest block was Sheshan, thereby forming a
tectonic type area of uplift and subsidence in the basin. The fault subsidence extent was
inconsistent, and the extent of fault subsidence was large in the south and relatively small
in the north. From the distribution of surface subsidence, the subsidence areas were all
located in the southern urban area accordance with the Lake Dian, i.e., areas A, B, C, D and
E marked in the blue rectangle in Figure 3.

5.3. Impact Analysis of Ground Subsidence and Groundwater Extraction

Given the survey, the groundwater extracted in the area mainly contained pore water,
karst water and fissure water. It was primarily employed for industrial, construction,
domestic, etc. As April 2021, Kunming has been using overall 1711 underground wells (i.e.,
1498 cold water wells and 213 hot water wells). Cold water wells were largely distributed
in five districts (i.e., Chenggong, Guandu, Xishan, Panlong and Wuhua districts), and hot
water wells had the major distributions in Guandu District and Xishan District. Since
2008, Kunming city has issued a series of regulations and documents to strengthen the
protection, utilization and management of groundwater. On the whole, 441 cold water
wells have been closed (e.g., 369 in the main urban area of Dianchi Basin). Although the
groundwater utilization has slowed down to a certain extent, groundwater is still over-
utilized. The statistical table of groundwater level changes in the study area (Table 4) was
obtained the collation of historical information from the environmental monitoring stations
in Yunnan Province, and the corresponding groundwater exploitation sites in the blocks
are shown in Figure 8a. In the early 1980s, the burial depth of the hot water level was
4~6 m, and the single well gushing water volume was 1000~1700 m3 /d. However, the
current water level is generally lower than the water level when the well was completed,
and the average burial depth of the water level has increased to 11~25 m; the maximum
decline value is more than 35 m, the average decline rate is 1.80 m/a, and the maximum is
4.59 m/a [46]. The groundwater exploitation and water level landing funnels in Kunming
over the years confirmed that there was a clear correlation between ground subsidence
and the continuous decline of groundwater levels [37]. Since the groundwater exploitation
data could not be obtained, we used the groundwater data released by NASA to obtain
groundwater data for the same period, and analyzed four regional subsidence feature
points. The point distribution of the four feature points ( J1 ∼ J5 ) is illustrated in Figure 8a.
The results of the analysis are shown in Figure 8b, and the cumulative subsidence values of
the four points have obvious correlation with the seasonal changes of groundwater level.
In addition, Kunming was located in the crustal active zone, covering the intermountain
basin, valley basin pores, fractures, karst water area with a wide range of carbonate rock
and karst caves distributed and with funnels scattered all over the city. The underground
confined water in the bedrock exhibits a certain supporting function. Under the excessive
exploitation of groundwater and the soil weight pressure, the bedrock is easy to lose its
support, thereby causing several geological disasters (e.g., ground collapse, foundation
subsidence and earthquake).
Sustainability 2022, 14, x FOR PEER REVIEW 14 of 22

Sustainability 2022, 14, 12387

accordance with the Lake Dian, i.e., areas A, B, C, D and E marked in the blue14rectangle
of 21
in Figure 3.

Figure
Figure 7. 7. Velocityand
Velocity and cross-section
cross-section diagram
diagramof fault in the
of fault instudy area. (a–c)
the study area.are plots of
(a)–(c) aredisplacements
plots of displace-
0 , P P0 and P P0 , respectively.
along profiles P
ments along profiles𝑃
1 P1 𝑃 2 , 𝑃 𝑃 and
2 3 3 𝑃 𝑃 , respectively.

5.3.Table 4. The
Impact statisticsoffor
Analysis the exploitation
Ground Subsidenceof groundwater for different
and Groundwater blocks.
Extraction
NumberGiven the survey,
of Mining theSpacing
Minimum groundwater
between extracted inDrawdown
Cumulative the areaofmainlyAverage
contained
Drop inpore wa-
Block
Wells
ter, (Eyes)
karst Production
water and fissure WellS (m)
water. It was Wateremployed
primarily Level (m) Water Levelconstruction,
for industrial, (m/a)
Urban area 28
domestic, etc. As April 2021, 150 Kunming has been20.22~25.10 2.26
using overall 1711 underground wells
Guanshang-Jinmasi 34 100 9.18~24.88 1.32
Yangfangao-Paomashan (i.e., 1498
21 cold water wells and
200 213 hot water wells). Cold
5.70~9.38 water wells were
1.18 largely dis-
Haigeng Sanatorium 46 150 11.50~21.30
tributed in five districts (i.e., Chenggong, Guandu, Xishan, Panlong and Wuhua districts), 1.62

andNote:
hotThe data cited
water wellsfromhad
the literature
the major[37,46].
distributions in Guandu District and Xishan District.
Since 2008, Kunming city has issued a series of regulations and documents to strengthen
the protection, utilization and management of groundwater. On the whole, 441 cold water
wells have been closed (e.g., 369 in the main urban area of Dianchi Basin). Although the
groundwater utilization has slowed down to a certain extent, groundwater is still over-
utilized. The statistical table of groundwater level changes in the study area (Table 4) was
 Sustainability
Sustainability 2022,
2022, 14,14,12387
x FOR PEER REVIEW 1615ofof2221

Figure8.8.Analysis
Figure Analysis diagram
diagram of
of ground
ground subsidence
subsidence and
and groundwater extraction. (a)
groundwater extraction. (a)isisthe
thelocation
location
map of Groundwater exploitation site and four feature points ( J ∼
map of Groundwater exploitation site and four feature points (𝐽 1~𝐽 ). 5(b) is the analysis results.
J ). (b) is the analysis results.
Sustainability 2022, 14, x FOR PEER REVIEW 17 of 22
Sustainability 2022, 14, 12387 16 of 21

5.4. Impact Analysis of Land Subsidence and Rainfall

5.4. Impact Analysis of Land Subsidence and Rainfall
Kunming
Kunming is is in
inaalow-latitude
low-latitudeplateau
plateauwith
witha acomplex
complexand anddiverse
diversetopography
topography andanda
large terrain
a large height
terrain heightdifference. Obvious
difference. Obviousvertical and and
vertical horizontal differences
horizontal in climate
differences have
in climate
been identified. The annual precipitation can be divided into dry and wet
have been identified. The annual precipitation can be divided into dry and wet seasons in seasons in terms
of timeofdistribution.
terms The rainy
time distribution. The season is from
rainy season May to
is from May October, withwith
to October, about 85% 85%
about of the
of an-
the
nual
annual precipitation. The dry season from November to April accounts for only 15% of the
precipitation. The dry season from November to April accounts for only 15% of the
annual
annual precipitation.
precipitation. ToTo explore
explore the
the correlation
correlation between
between land land subsidence
subsidence rule
rule and
and rainfall
rainfall
in
in Kunming,
Kunming, this
this study
study selected
selected five
five feature
feature points
points distributed
distributed in in vital
vital subsidence
subsidence areasareas in
in
the subsidence rate map obtained by SBAS technology for analysis and
the subsidence rate map obtained by SBAS technology for analysis and statistics. The point statistics. The
point distribution
distribution of theof thefeature
five five feature ( S1 ∼(𝑆S5~𝑆
pointspoints ) is )illustrated
is illustrated in Figure
in Figure 9. 9.

Figure
Figure 9.
9. Distribution
Distribution map of 𝑆S1~𝑆∼ feature
map of points.
S5 feature points.

The precipitation
The precipitationdata datafrom fromJuly
July2018
2018toto November
November 2020
2020 in the
in the studystudy
areaarea
werewere
ac-
acquired
quired from
from China
China Meteorological
Meteorological DataData Network
Network (http://data.cma.cn/,
(http://data.cma.cn/, accessed accessed
on 21 Au- on
21 August
gust 2021).2021). The time-series
The time-series for 𝑆 for
valuesvalues ~𝑆 S1feature
∼ S5 feature
points points and monthly
and monthly precip-
precipitation
itation
data datastudy
in the in the studyare
period period
listedare listed 5.
in Table inThe
Table 5. The time-series
time-series subsidencesubsidence
values of thevalues
five
of the five feature points in the study period were compared with
feature points in the study period were compared with the monthly average precipitation the monthly average
precipitation
(Figure (Figure 10).
10). According to theAccording to the
figure, the figure, the
subsidence subsidence
time-series time-series
of the of the
five feature five
points
feature points showed a nonlinear decline with seasonal variations,
showed a nonlinear decline with seasonal variations, and there was continuous ground and there was contin-
uous ground
subsidence at subsidence
𝑆 point during at S1 thepoint during
whole theperiod.
study whole study period. 2018,
In November In November 2018,
the precipita-
the precipitation in Kunming city was 0, and the subsidence value
tion in Kunming city was 0, and the subsidence value of that month achieved an extreme of that month achieved
an extreme
value of −34.4 value −34.4time,
mm.ofOver mm.the Over time, thesubsidence
cumulative cumulativerate subsidence
reached ratethe reached
maximalthe of
maximal of − 160.7 mm in August 2020, and then the subsidence
−160.7 mm in August 2020, and then the subsidence slowed down. From July 2018 slowed down. From July
to Jan-
2018 to
uary January
2019, 𝑆 , 𝑆 2019, S2 , S𝑆3 , Srose
, 𝑆 and S5 rose
4 andfirst andfirst
thenand then subsided,
subsided, duringduring
whichwhich the overall
the overall sub-
siding amount was relatively small. Subsequently, they continued to subside. In the In
subsiding amount was relatively small. Subsequently, they continued to subside. the
rainy
rainy season of 2019, the subsiding rate slowed down and then continued
season of 2019, the subsiding rate slowed down and then continued to subside. From June to subside. From
June
to to November
November 2020,2020, the subsiding
the subsiding tended
tended to be
to be gentle.
gentle. AsAs depictedFigure
depicted Figure10,10,the
theperiod
period
when the selected subsidence feature points decreased in the five vital subsidence areas
when the selected subsidence feature points decreased in the five vital subsidence areas
corresponded to the rainy season with the maximal precipitation, thereby demonstrating
corresponded to the rainy season with the maximal precipitation, thereby demonstrating
that precipitation is inversely correlated with land subsidence. This is because rainfall is
that precipitation is inversely correlated with land subsidence. This is because rainfall is
capable of penetrating into the ground through the soil and replenishing the groundwater
capable of penetrating into the ground through the soil and replenishing the groundwater
level. With the rise of the groundwater level and the increase in the soil water content, the
level. With the rise of the groundwater level and the increase in the soil water content, the
support of bedrock was improved, thereby effectively down-regulating the subsidence rate
support of bedrock was improved, thereby effectively down-regulating the subsidence
of the ground. In spring and winter, however, drought and little rain, soil water content
Sustainability 2022, 14, 12387 17 of 21

decreased, groundwater extraction increased, groundwater level declined, bedrock lost

support again, thereby intensifying land subsidence.

Table 5. Subsidence times series of S1 ~S5 feature points and monthly average rainfall data.

Displacement (mm) Precipitation

Date
S1 S2 S3 S4 S5 (mm)

20 July 2018 0 0 0 0 0 172.8

13 August 2018 −2 8 4 5 6 289.6
18 September 2018 −8 6 1 −1 11 82.8
12 October 2018 −13 4 2 −4 8 53.3
17 November 2018 −34 −3 −3 −11 8 0.0
11 December2018 −23 −8 −6 −17 4 39.3
16 January 2019 −32 −13 −4 −22 3 31.7
21 February 2019 −40 −20 −8 −29 −12 5.7
29 March 2019 −53 −28 −15 −38 −28 4.5
22 April 2019 −60 −32 −18 −45 −39 7.1
28 May 2019 −75 −43 −26 −61 −58 20.1
21 June 2019 −86 −48 −28 −68 −63 94.1
15 July 2019 −90 −53 −26 −70 −62 279.5
20 August 2019 −99 −59 −26 −75 −71 120.7
25 September 2019 −105 −64 −29 −79 −69 186.1
19 October 2019 −113 −72 −31 −84 −77 74.2
24 November 2019 −116 −79 −34 −86 −84 7.5
30 December 2019 −124 −86 −37 −91 −91 9.0
23 January 2020 −127 −88 −38 −95 −94 61.5
28 February 2020 −134 −96 −42 −99 −99 24.2
23 March 2020 −140 −102 −44 −103 −101 5.3
28 April 2020 −141 −109 −48 −109 −106 45.3
22 May 2020 −148 −113 −53 −113 −112 *
27 June 2020 −156 −120 −54 −119 −118 122.0
21 July 2020 −156 −122 −54 −121 −120 379.8
26 August 2020 −161 −126 −54 −125 −126 258.0
19 September 2020 −158 −125 −55 −126 −128 96.2
25 October 2020 −155
Sustainability 2022, 14, x FOR PEER REVIEW −133 −57 −129 −130 26.9 19 of 22
18 November 2020 −155 −134 −57 −128 −128 6.6
* Indicates that the current month data is unmeasured.

Figure 10.
Figure Correlationbetween
10.Correlation betweenthe
thetime
time sequence
sequence ofof Subsidence
Subsidence feature
feature points
points andand monthly
monthly aver-
average
rainfall.
age rainfall.

6. Conclusions
With the development of urbanization, land subsidence has become one of the most
prominent geological disasters in Kunming city. Land subsidence has caused significant
damage to houses, roads, pipelines and other types of infrastructure. In this study, two
Sustainability 2022, 14, 12387 18 of 21

6. Conclusions
With the development of urbanization, land subsidence has become one of the most
prominent geological disasters in Kunming city. Land subsidence has caused significant
damage to houses, roads, pipelines and other types of infrastructure. In this study, two time-
series InSAR techniques were used to invert the time-series deformation of the Kunming
main urban area from July 2018 to November 2020. Then, the accuracy of the time-series
InSAR technology was verified using three methods: cross-comparison, comparison with
the level data of the same period and field investigation. Finally, combined with time-
series deformation, groundwater use survey data and rainfall data, this study analyzed
the causes of land subsidence in the Kunming main urban area from four aspects of urban
subway construction, geological structure, groundwater use and rainfall comprehensively.
It established the monitoring and analysis system of urban land subsidence. Based on this
study, the following conclusions are drawn:
(1) The deformation results achieved by applying PS-InSAR and SBAS-InSAR techniques
were consistent. There were numerous significant subsidence areas, mainly dis-
tributed in Xishan District, Guandu District and Jinning County. In the study period,
the maximal annual subsidence rates identified by the two methods were −78 mm/a
and −88 mm/a, respectively. The overall subsidence trend was large in the south
(along the Lake Dian) and relatively small in the north. Moreover, cross-comparison
verification was adopted to verify the InSAR monitoring results reliability, and the
correlation coefficient R2 of the two methods was 0.997. Compared with the identical
period leveling data, root mean square error (RMSE) reached 8.5 mm/a and 6.7 mm/a,
respectively. Through the field sampling investigation, the two InSAR techniques
were verified with high correlation and reliability.
(2) The ground subsidence of Kunming city is correlated with subway construction,
geological structure and groundwater extraction. Kunming tectonic movement is
significant, especially under the control by Chenggong-Fuming fault in North West.
Under the special geological structure condition, considerable urban construction
altered the original foundation stress state. With the increase in the load of the
foundation and the excessive exploitation underground confined water in bedrock,
the underground water level decreased, which lost their support for the foundation.
Compression consolidation of soil under the double action of self-weight and external
load induced regional ground subsidence.
(3) The surface subsidence of Kunming city showed seasonal variations with rainfall.
The flood season lasts from May to October. Moreover, the rainfall could replenish
the soil water content and groundwater level effectively, and the land subsidence rate
decreased. From November to following year April, the drought and lack of rain
increased water consumption and increased land subsidence rate.
In this study, the land subsidence mechanism in Kunming city during the study period
was proved qualitatively and quantitatively. The application value can be summarized
according to two aspects. First, we overcame traditional surface subsidence monitoring
technology disadvantages and improved the spatial and temporal resolution, which signifi-
cantly reduced the economic costs. It can realize small-deformation accurate monitoring in
various urban areas and provide strong technical support for monitoring surface subsidence
in urban areas. Second, the causes of Kunming city ground subsidence were analyzed
in this study from different aspects. A subsidence monitoring and analysis system was
established in this study to solve the problems of subsidence reason is not specific and
the subsidence trend is not clear in urban areas, which is beneficial to the inversion of
the subsidence trend in urban areas. However, two problems in this study that should be
discussed in subsequent studies. First, when the PS-InSAR and SBAS-InSAR technologies
were adopted to verify the internal coincidence accuracy, too few corresponding points
were selected. Distribution of corresponding points were concentrated relatively and not
distributed evenly in the whole research area. The calculated correlation coefficient was
Sustainability 2022, 14, 12387 19 of 21

close to 1, which should be corrected in subsequent studies. Second, in the process of

data processing, atmospheric correction is not carried out. The influence of atmosphere on
experimental results is ignored, which should be taken into account in future studies.

Author Contributions: Investigation, B.X., D.Z. and Z.Z.; formal analysis, B.X.; validation, B.X.,
J.Z. and Z.Z.; software, B.X. and D.Z.; methodology, B.X. and W.X.; Conceptualization, B.X. and
J.Z.; resources, J.Z.; project administration, D.L.; supervision, J.Z. and D.L.; visualization, B.X. and
D.Z.; writing—review and editing, B.X. and D.L.; writing—original draft preparation, B.X.; data
curation, J.Z.; funding acquisition, J.Z. All authors have read and agreed to the published version of
the manuscript.
Funding: This study was funded by Natural Science Foundation of China, grant number 41761081,
42161067.
Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.
Data Availability Statement: Not applicable.
Acknowledgments: Three arc-seconds SRTM DEM were freely downloaded from the website http:
//rmw.recordist.com/index.html, accessed on 10 May 2021. Thanks to ESA for its Sentinel-1A data
from July 2018 to November 2020 and China Meteorological Data Network (http://data.cma.cn/,
accessed on 1 April 2020 ) for its monthly precipitation data from July 2018 to November 2020.
Conflicts of Interest: The authors declare no conflict of interest.

References
1. Strozzi, T.; Wegmuller, U.; Tosi, L.; Bitelli, G.; Spreckels, V. Land subsidence monitoring with differential SAR interferometry.
Photogramm. Eng. Remote Sens. 2001, 67, 1261–1270.
2. Xu, J.Q.; Ma, T.; Lu, Y.K.; Bai, W.M.; Zhao, S.H. Land subsidence monitoring in north Henan plain based on SBAS-InSAR
technology. J. Jilin Uni. (Earth Sci.) 2019, 49, 1182–1191. [CrossRef]
3. Dong, S.C.; Zhong, Y.H.; Hu, H.; Huang, L.L. Ground subsidence monitoring during 2015–2018 in Changzhou based on times
series InSAR method. J. NanJing Unvi. (Natrual Sci.) 2019, 55, 370–380. [CrossRef]
4. Xu, Q.; Pu, C.H.; Zhao, K.Y.; He, P.; Zhang, H.Y.; Liu, J.L. Time series InSAR monitoring and analysis of spatiotemporal evolution
characteristics of land subsidence in Yan’an new district. Geom Inf Sci Wuhan Univ. 2021, 46, 957–969. [CrossRef]
5. Gong, L.X.; Zhang, J.F.; Guo, Q.S. Measure groundwater pumping induced subsidence with D-InSAR. IEEE Int Symp Geosci
Remote Sens. 2005, 7, 25–29. [CrossRef]
6. Luo, X.G.; Wang, J.; Xu, Z.; Zhu, S.; Meng, L.S.; Liu, J.K.; Cui, Y. Dynamic analysis of urban ground subsidence in Beijing based on
the permanent scattering InSAR technology. J. Appl. Remote Sens. 2018, 12, 026001. [CrossRef]
7. Wu, Q.; Jia, C.T.; Chen, S.B.; Li, H.Q. SBAS-InSAR Based Deformation Detection of Urban Land, Created from Mega-Scale
Mountain Excavating and Valley Filling in the Loess Plateau: The Case Study of Yan’an City. Remote Sens. 2019, 14, 1673.
[CrossRef]
8. Zebker, H.A.; Goldstein, R.M. Topographic mapping from interferometric synthetic aperture radar observations. J. Geophys. Res.
1986, 91, 4993–4999. [CrossRef]
9. Perna, S.; Wimmer, C.; Moreira, J.; Fornaro, G. Airborne DInSAR time series at X-band. IEEE Int. Symp. Geosci. Remote Sens. 2010,
2868–2871. [CrossRef]
10. Buonanno, S.; Zeni, G.; Fusco, A.; Manunta, M.; Marsella, M.; Carrara, P.; Lanari, R. A GeoNode-Based platform for an effective
exploitation of advanced DInSAR measurements. Remote Sens. 2019, 11, 2133. [CrossRef]
11. Carnec, C.; Fabriol, H. Monitoring and modeling land subsidence at the Cerro Prieto geothermal Field, Baja California, Mexico,
using SAR interferometry. Geophys. Res. Lett. 1999, 26, 1211–1214. [CrossRef]
12. Xia, Y.; Kaufmann, H.; Guo, X.F. Landslide monitoring in the Three Gorges area using D-InSAR and corner reflectors. Photogramm.
Eng. Remote Sens. 2004, 70, 1167–1172. [CrossRef]
13. Lubitz, C.; Motagh, M.; Wetzel, H.U.; Kaufmann, H. Remarkable urban uplift in Staufen im Breisgau, Germany: Observations
from TerraSAR-X InSAR and leveling from 2008 to 2011. Remote Sens. 2013, 5, 3082. [CrossRef]
14. Dong, J.; Zhang, L.; Tang, M.G.; Liao, M.S.; Xu, Q.; Gong, J.Y. Mapping landslide surface displacements with time series SAR
interferometry by combining persistent and distributed scatterers: A case study of Jiaju landslide in Danba, China. Remote Sens.
Environ. 2018, 205, 180–198. [CrossRef]
15. Li, M.H.; Zhang, L.; Shi, X.G.; Liao, M.S.; Yang, M.S. Monitoring active motion of the Guobu landslide near the Laxiwa
hydropower station in China by time-series point-like targets offset tracking. Remote Sens Environ. 2019, 221, 80–93. [CrossRef]
Sustainability 2022, 14, 12387 20 of 21

16. Bayik, C. Deformation analysis of 2020 mw 5.7 Karliova, Turkey, earthquake using DInSAR method with different incidence
angle SAR data. Arab J. Geosci. 2021, 14, 273. [CrossRef]
17. Tang, X.W.; Guo, R.M.; Xu, J.Q.; Sun, H.P.; Chen, X.D.; Zhou, J.C. Probing the Fault Complexity of the 2017 Ms7.0 Jiuzhaigou
earthquake based on the InSAR Data. Remote Sens. 2021, 13, 1573. [CrossRef]
18. Li, Y.J.; Li, Y.S.; Hu, X.P.; Liu, H.Q. Fault geometry and mechanism of the Mw 5.7 Nakchu earthquake in Tibet inferred from
InSAR observations and stress measurements. Remote Sens. 2021, 13, 5142. [CrossRef]
19. Yang, C.S.; Zhang, Q.; Zhao, C.Y.; Ji, L.Y. Small baseline bubset InSAR technology used in Datong ground subsidence, fissure and
fault zone monitoring. Geomat. Inf. Sci Wuhan Univ. 2014, 39, 945–950. [CrossRef]
20. He, X.F.; Gao, Z.; Xiao, R.Y.; Luo, H.B.; Feng, C. Monitoring and analysis of subsidence along Lian-Yan railway using multi-
temporal Sentinel-1A InSAR. Acta Geodaetica et Cartographica Sinica 2021, 50, 600–611. [CrossRef]
21. Zhu, B.Y.; Yao, F.Y.; Sun, J.W.; Wang, X. Attribution analysis on land subsidence feature in Hexi area of Nanjing by InSAR and
geological data. Geomat. Inf. Sci. Wuhan Univ. 2020, 45, 442–450. [CrossRef]
22. Ferretti, A.; Prati, C.; Rocca, F. Nonlinear subsidence rate estimation using permanent scatterers in differential SAR interferometry.
IEEE Trans Geosci Remote Sens. 2000, 5, 2202–2212. [CrossRef]
23. Ferretti, A.; Prati, C.; Rocca, F. Permanent scatterers in SAR interferometry. IEEE Trans Geosci Remote Sens. 2001, 39, 8–20.
[CrossRef]
24. Berardino, P.; Fornaro, G.; Lanari, R.; Sansosti, E. A new algorithm for surface deformation monitoring based on small baseline
differential SAR interferograms. IEEE Trans. Geosci. Remote Sens. 2002, 40, 2375–2383. [CrossRef]
25. Zhou, C.F.; Gong, H.L.; Chen, B.B.; Li, J.W.; Gao, M.L.; Zhu, F.; Chen, W.F.; Liang, Y. InSAR time-series analysis of land subsidence
under different land use types in the eastern Beijing plain, China. Remote Sens. 2014, 9, 380. [CrossRef]
26. Hu, B.; Chen, J.Y.; Zhang, X.F. Monitoring the land subsidence area in a coastal urban area with InSAR and GNSS. Sensors 2019,
19, 3181. [CrossRef]
27. Yalvac, S. Validating InSAR-SBAS results by means of different GNSS analysis techniques in medium- and high-grade deformation
areas. Environ Monit Assess. 2020, 192, 120. [CrossRef]
28. Hakim, W.L.; Achmad, A.R.; Lee, C.W. Land subsidence susceptibility mapping in Jakarta using functional and meta-ensemble
machine learning algorithm based on time-series InSAR data. Remote Sens. 2020, 12, 3627. [CrossRef]
29. Hussain, M.A.; Chen, Z.L.; Shoaib, M.; Shah, S.U.; Khan, J.K.; Ying, Z. Sentinel-1A for monitoring land subsidence of coastal city
of Pakistan using Persistent Scatterers InSAR technique. Sci. Rep. 2022, 12, 5294. [CrossRef]
30. Orhan, O. Monitoring of land subsidence due to excessive groundwater extraction using small baseline subset technique in
Konya, Turkey. Environ Monit Assess. 2021, 193, 174. [CrossRef]
31. Solano-Rojas, D.; Wdowinski, S.; Cabral-Cano, E.; Osmanoğlu, B. Detecting differential ground displacements of civil structures
in fast-subsiding metropolises with interferometric SAR and band-pass filtering. Sci. Rep. 2020, 10, 15460. [CrossRef]
32. Orhan, O.; Oliver-Cabrera, T.; Wdowinski, S.; Yalvac, S.; Yakar, M. Land subsidence and its relations with sinkhole activity in
Karapınar region, Turkey: A multi-sensor InSAR time series study. Sensors 2021, 21, 774. [CrossRef]
33. Zhu, W.; Li, W.L.; Zhang, Q.; Yang, Y.; Zhang, Y.; Qu, W.; Wang, C.S. A decade of ground deformation in Kunming (China)
revealed by multi-temporal synthetic aperture radar interferometry (InSAR) technique. Sensors 2019, 19, 4425. [CrossRef]
34. Shao, J.M.; Li, J.P.; Lu, H.J. Application of PS-InSAR technology in ground subsidence monitoring of Kunming. J. Henan Univ.
Urb. Constr. 2017, 26, 45–49. [CrossRef]
35. Xiong, P.; Zuo, X.Q.; Li, Y.F.; You, H. Application of dual-polarized Sentinel-1 data to subsidence monitoring in Kunming. Prog.
Geophys. 2020, 4, 1317–1322. [CrossRef]
36. Jiang, N.Q.; Zuo, X.Q.; Wang, Z.H.; You, H.; Zhao, Y.H. Monitoring land settlement in the main urban area of Kunming based on
PS-InSAR and SBAS-InSAR. J. Guizhou Univ. (Natural Sci). 2020, 37, 72–78. [CrossRef]
37. Xue, C.D.; Liu, X.; Li, B.Z.; Meng, G.T.; Li, F. Mechanism analysis of land subsidence in Kunming city area. Chin. J. Geol. Hazard
Control 2004, 15, 47–53. [CrossRef]
38. Perissin, D.; Prati, C.; Engdahl, M.E.; Desnos, Y.L. Validating the SAR wavenumber shift principle with the ERS Envisat PS
coherent combination. IEEE Trans. Geosci. Remote Sens. 2006, 44, 2343–2351. [CrossRef]
39. Foroughnia, F.; Nemati, S.; Maghsoudi, Y.; Perissin, D. An iterative PS-InSAR method for the analysis of large spatio-temporal
baseline data stacks for land subsidence estimation. Int. J. Appl. Obs. Geoinformation. 2019, 74, 248–258. [CrossRef]
40. Colesanti, C.; Ferretti, A.; Novali, F.; Prati, C.; Rocca, F. SAR monitoring of progressive and seasonal ground deformation using
the permanent scatterers technique. IEEE Trans. Geosci. Remote Sens. 2003, 41, 1685–1701. [CrossRef]
41. Lanari, R.; Mora, O.; Manunta, M.; Mallorquí, J.; Berardino, P.; Sansosti, E. A small-baseline approach for investigating
deformations on full-resolution differential SAR interferograms. IEEE Trans. Geosci. Remote Sens. 2004, 42, 1377–1386. [CrossRef]
42. Yu, S.Y.; Yang, Y.Y.; Zhang, P.F.; Sun, J.; Luo, J.J. Monitoring Land Subsidence and Fault Activity in Hefei City Based on MT-InSAR.
J. Geod. Geodyn. 2020, 41, 398–402. [CrossRef]
43. Qu, C.Y.; Dan, X.J.; Zhang, G.H.; Xu, X.B.; Song, X.G.; Zhang, G.F.; Liu, Y.H. The research progress in measurement of fault
activity by time series InSAR and discussion of related issues. Seismol. Geol. 2014, 36, 731–748.
44. Liu, Q.; Yue, G.S.; Ding, X.B.; Yang, K.; Feng, G.C.; Xiong, Z.Q. Temporal and Spatial Characteristics Analysis of Deformation
Along Foshan Subway Using Time Series InSAR. Geomat. Inf. Sci. Wuhan Univ. 2019, 7, 1099–1106.
Sustainability 2022, 14, 12387 21 of 21

45. Ding, X.F.; Zhang, G.D. An analysis of the geologic condctions of Haigeng geothermal well in Kunming, Yunnan. Yunnan Geol.
2019, 38, 238–243.
46. Zhou, D.Y.; Zuo, X.Q.; Zhao, Z.F. Constructing a large-scale urban land subsidence prediction method based on neural network
algorithm from the perspective of multiple factors. Remote Sens. 2022, 14, 1803. [CrossRef]