SBAS-InSAR-Based Analysis of Surface Deformation in the Eastern Tianshan Mountains, China

Due to the unique geographical characteristics of cold alpine and high-altitude regions, glaciers, permafrost, ground ice, rock glaciers, and other periglacial geomorphology have developed with fragile habitats, and these areas are often the birthplaces of many river basins and natural hazards. With global warming and the extensive cryogenesis and physical weathering, the thermal state of permafrost and the mass balance of glaciers have been changed, and thus it can be deduced that the surface deformation is of great concern. To obtain ground subsidence or uplift over a large area to understand local surface changes, the small baseline subset interferometric synthetic aperture radar (SBAS-InSAR) technique was applied to process 89-scene of Sentinel-1A images ranging from December 25, 2017 to January 2, 2021 to obtain surface deformation for these 3 years for the eastern Tianshan Mountains, China. The surface deformation characteristics of the area were analyzed to provide a basic dataset for environmental protection policies and mitigation or reduction of natural hazards in this region, and to verify the applicability of SBAS-InSAR technology in alpine and high-altitude areas. The results show that the SBAS-InSAR technique processing with sentinel-1A dataset cannot be effectively used to acquire ground deformation in areas covered by trees, scrub/shrub, glaciers, snow, and ground ice, where the decohered phenomenon is serious. In other regions, SBAS-InSAR can effectively measure surface subsidence or uplift. Surface deformation is significant throughout the study area, with rates ranging from −70.7 to 50.8 mm/a and with an average rate of 1.1 mm/a. There are obvious regions of uplift in the northwest, northeast, and central sections of the study area, with uplift greater than 155.73 mm in 3 years, and three obvious regions of subsidence in the northeast and west sections of the study area, with subsidence of at least −125.20 mm in 3 years. The remaining areas of deformation are scattered, with smaller amounts of settlement and uplift and with an isolated and sporadic distribution. Areas with elevations of 3,150 to 4,275 m.a.s.l., slopes of 15°–50°, and southwest, west, and northwest aspects are geologic disaster-prone regions and should receive more attention and more field monitoring. The results of this study have important implications for local environmental protection and hazard prevention.


INTRODUCTION
Cold alpine and high-altitude regions are characterized by fragile ecosystem habitats where periglacial geomorphology such as glaciers, permafrost, rock glaciers, and subsurface ice is richly developed due to the perennial low-temperature environment . In recent years, with global warming (Hansson et al., 2021), the periglacial geomorphology of high alpine-altitude regions has been put at great risk. This warming will accelerate the degradation of glaciers, permafrost, rock glaciers, snow cover, and subsurface ice (Zhu et al., 1996;Shan et al., 2015;Cheng et al., 2019;Zhao et al., 2019), exposing their subsurface land and touching off continuous uneven subsidence (Qin et al., 2018). This may lead to geologic hazards such as landslides, collapses, and debris flows (Metternicht et al., 2005;Shan et al., 2014), which may entail huge economic losses.
The eastern Tianshan Mountains in China are located in a typical alpine and high-altitude region with a developed periglacial geomorphology and an abundance of minerals (Du et al., 2020a;Du et al., 2021c). To ensure safe conduct of mining and to understand the subsidence of the regional surface, 89 acquisitions of Sentinel-1A images covering the study area from December 25, 2017 to January 2, 2021 were selected and processed by SBAS-InSAR to obtain the ground deformation of the study area during this period. The results can provide a theoretical basis for the formulation of regional development policies, as well as safety recommendations for mining activities within the region.

OUTLINE OF STUDY AREA
The study area is situated in the eastern section of the Tianshan Mountains in China ( Figure 1) and belongs to the Xinjiang Uyghur Autonomous Region, with a geographical location of 43.25-43. 62°N, 84.87-85.49°E, an area of 1,590.78 km 2 , and an altitude of 2,549-5,250 m above sea level (m.a.s.l.), with an average altitude of 3,839 m.a.s.l. (Figure 1B), making this a typical cold-alpine, high-altitude region (Du et al., 2020a;Du et al., 2021a). According to a land-cover dataset, Esri 2020 Global Land Cover (Karra et al., 2021) (Figure 1D), and site inspection, the ground surface cover of the study area is mainly bare land, followed by snow/ice and glaciers (Figures 1C,D) and grass ( Figure 1D). In addition, there are also some other land cover types including water, scrub/shrub, trees, crops, and built areas. The terrain is generally high in the north and low in the south, with snow cover and ice developing in the alpine-mountain areas in the north (3,971-5,250 m.a.s.l.) and grassland widely distributed in the lower, flatter areas in the south (2,549-3,628 m.a.s.l.) and mostly accompanied by rivers, wetlands (scrub/shrub) and trees in the southeast section.
The amount of thawing of snow cover and ice affects the succession and growth of the local grassland ecosystem (Du et al., 2020b), which has a significant influence on the standard of living of pastoralists living in this area. Snow cover and ice are also a source of water resource recharge for many rivers in the adjacent low-altitude areas, and their phase behavior change is a guarantee of water supply for people living and working downstream along the rivers. As global warming accelerates the thawing of snow and ice (Qin et al., 2018), the exposed subsurface land will accelerate subsidence of the regional ground surface. Excessive melting in till sedimentation (moraine) can alter local topography, thereby affecting regional runoff and the water cycle. The amount of surface subsidence or uplift is also related to the frequency of natural disasters such as landslides, collapses, and debris flows. Du et al. (2021b) obtained deformation information from the open-pit mining zone in the lower right corner of the study area (43.28-43.35°N, 84.95-85.12°E) using the SBAS-InSAR technique and compared the results with field GNSS monitoring data. This effort demonstrated the outstanding measurement capability of SBAS-InSAR in this region and obtained credible results. Based on this work, the parameters involved in the SBAS-InSAR processing for this paper were consistent with their settings in the literature (Mora et al., 2002;Hu et al., 2021a;Du et al., 2021b;Kursah et al., 2021), ensuring the reliability of the results, and therefore it was unnecessary to verify the accuracy of the measurement results in this study. This paper focuses on obtaining surface settlement information through SBAS-InSAR processing and analyzing settlement or uplift characteristics to understand the deformation features of the study area over a 3-year period, which will provide basic data for regional development policy formulation. The obtained surface deformation information can be used as an input risk assessment dataset for natural disasters and other studies, thus offering great significance and usefulness.

DATASET AND METHODOLOGY Dataset
The 89-scene of Sentinel-1A ascending (satellite travelling south to north) acquisitions from December 25, 2017 to January 2, 2021 were downloaded from the ASF website (https://vertex.daac. asfalaska.edu/), along with the precise orbital data corresponding to each acquisition as obtained from the ESA website (https://qc.sentinel1.eo.esa.int/). In general, longwavelength (low-frequency) bands have good long-range performance and easy access to high-power transmitters and huge antennas (Ferretti et al., 2007). Short-wavelength (highfrequency) bands generally give precise distances and positions, but have a short range of action . Sentinel-1A has been acquiring data since October 2014 and Sentinel-1B since September 2016 with c-band (4-8 GHz) (Zhu et al., 2019), which has some penetration while maintaining good monitoring capability (high spatial accuracy). Each sentinel-1 constellation is in a near-polar, sun-synchronous orbit, with a 12-day repeat cycle and the two-satellite constellation offers a 6-day repeat cycle, which ensures excellent monitoring capability (high temporal resolution). They are very commonly used as SAR image data and are free and open-source, with multiple modes Frontiers in Earth Science | www.frontiersin.org November 2021 | Volume 9 | Article 729454 (Stripmap, SM; Interferometric Wide swath, IW; Extra Wide swath, EW; Wave, WV) and multi-polarization (HH, VV, HH + HV, VV + VH) (Yang et al., 2015). The data used in this paper were the single-look complex (SLC) Level 1 product in IW mode, which acquires 250 km of surface data at a spatial resolution of 5 × 20 m (single scene) (Yang et al., 2015) and contains both phase and amplitude information. Phase is a function of time, and distance measurements can be obtained based on phase and velocity information . Based on this, SAR images can be used for distance measurement and deformation observation (Ye, 2010;Hu et al., 2021b). To reduce the time spent on data processing, only the extent of images containing the study area was selected for processing in this paper. DEM data were provided by the Shuttle Radar Topography Mission (SRTM) 1 Arc-Second Global DEM from the USGS Earth Explorer website (https://earthexplorer.usgs.gov/) with a spatial resolution of 30 m, which can be used to produce aspect and slope dataset products ( Figure 2) and extract contour lines. Glacier distribution data were derived from the Second Glacier Inventory of China, which can be obtained from the National Tibetan Plateau/Third Pole Environment Data Center (TPDC) (https://data.tpdc.ac.cn/) in vector format.
A ground-surface land coverage dataset ( Figure 1D), Esri 2020 Global Land Cover (Karra et al., 2021) is based on the dataset produced for the Dynamic World Project by National Geographic Society in partnership with Google and the World Resources Institute. The map is derived from ESA Sentinel-2 imagery at 10 m resolution. It is a composite of LULC predictions for 10 classes throughout the year in order to generate a representative snapshot of 2020 and was produced by a deep learning model trained using over five billion hand-labeled Sentinel-2 pixels, sampled from over 20,000 sites distributed across all major biomes of the world. This dataset is now officially available for global public access and can be downloaded freely through the Land Cover Downloader Map website page (https://www.arcgis.com/apps/instant/media/index. html?appid fc92d38533d440078f17678ebc20e8e2 or https:// ai4edataeuwest.blob.core.windows.net/io-lulc/io-lulc-model-001-v01-composite-v03-supercell-v02-clip-v01.zip). Compared to FROM-GLC10 dataset (Gong et al., 2019), Esri 2020 Global Land Cover is closer to the date of acquisition of the SAR images (December 25, 2017 to January 2, 2021) of the study area and has a higher accuracy of ground features classification.

SBAS-InSAR
Since its introduction in 2002 (Mora et al., 2002), SBAS-InSAR technology has been used in many surface deformation monitoring studies. The principles, processing procedures, and other features of SBAS-InSAR have been extensively described in the literature (A. J. Hooper, 2008;Liu et al., 2019;Hu et al., 2021b) and will not be repeated in this paper. In this study, the SBAS-InSAR process was implemented on the SARscape 5.2.1 software platform, and its parameter settings and process were consistent with those in the literature (Chen et al., 2013;Du and Zhao, 2020;Reinosch et al., 2020;Du et al., 2021a;Du et al., 2021b;Du et al., 2021d). Table 1 gives the parameter settings of each step. The annual rate (mm/a) for the study area during the 3-year period from December 25, 2017 to January 2, 2021 and the surface cumulative deformation information (mm) will be available once the SBAS-InSAR processing has been completed. Deformation includes subsidence and uplift and represents cumulative deformation information describing ground surface changes for each relevant period starting on December 25, 2017 and ending with the date of image acquisition. Many important intermediate files used during processing will also be available, such as the coherence coefficient graphs, deformation maps of the first estimate, deformation maps after the second estimate, and other important reference files. In the InSAR process, the coherence coefficient is an important parameter ranging from 0 to 1 that determines the accuracy of the measurement results. The magnitude of the coherence coefficient indicates the extent of decohered, and the coherence coefficient can also be used as a threshold to remove low-coherence areas to improve the accuracy of the monitoring results.
Generally, if the coherence is less than 0.2 Hu et al., 2021b), there may be few decohered regions, but the accuracy of the measurement results will be reduced. The coherence threshold will mask the areas where the coherence is less than the threshold value. Masked regions will be shown as blank in the result maps. In this study, the threshold was set to 0.35 to ensure the precision of the results.

Deformation Characterization
Deformation velocity maps and cumulative deformation maps were analyzed using the ArcGIS 10.6 software platform to identify the spatial and temporal deformation variation characteristics of the study region, focusing on the consistency and heterogeneity of the study area at overall and local scales. Regions with intense subsidence or uplift changes within the study area were analyzed and zoned in detail, and a set of observation points (OPs) was selected to extract the cumulative deformation for every period. Based on the Land Cover dataset, we have statistically analyzed the surface deformation characteristics of each type of ground feature, including the velocity range, decohered area, and decohered area percentage. This involved reclassifying the deformation velocity and transforming the raster results to vector format to extract the areas of each class, followed by reanalysis and counting of the land-cover area in each class.
The Raster to Point tool was used to transfer the acquired deformation velocity result to vector format points which were used as a baseline to extract other corresponding values including slope, aspect, elevation, and surface feature using the Extract Values to Points tool. The velocity map can also be transferred to vector polygon for area calculation using the Raster to Polygon tool. However, as the Raster to Polygon tool only supports data in integer, it is necessary to convert the results of the data type float to integer before raster to polygon processing. The extraction and processing of the relevant data corresponding to each ground cover data is similar.
The overlay method was used to analyze the relationship between deformation and land cover, as well as to collect deformation information for each land-cover type. More attention was paid to areas with intense variation and steep slopes to explore the intrinsic relevance between them. Slope data were obtained by calculating the DEM.
During the analysis, all the processes just described were not separate, but rather intertwined with each other. These analyses provided an understanding of the ground surface change characteristics of the study area and a preliminary analysis of deformation factors. The results provide a theoretical basis for formulating regional development policies as well as basic data for physical geological hazard prevention.

Deformation Velocity and Cumulative Deformation
A velocity map ( Figure 3) and cumulative deformation maps (Figures 5-7) for the 3 years were created once SBAS-InSAR processing had been completed. Warm colors indicate surface deformation away from the sensor line-of-sight (LOS) direction, whereas cool colors indicate surface deformation toward the sensor LOS direction. Blank regions represent missing data due to decohered. All the deformation descriptions information including the uplift and subsidence about the velocity and cumulative deformation in this paper are based on the sensor (LOS) direction. Figure 3 shows regions of significant subsidence in the northeastern, southwestern, and southeastern sections of the study area, with rates of subsidence of −30.2 mm/a or more. The northern, central, and southern parts of the study area contain regions of significant uplift, with uplift rates greater than 15.6 mm/a. Areas with obvious deformation are isolated and sporadic. In the southwestern section, there is an open-pit mine, where Du et al. (2021b) have analyzed the associated ground deformation and obtained results in accord with the velocity results comparing with the situ global navigation satellite system (GNSS) measurements. The processing and parameter settings in this paper are identical to those of the literature (Du et al., 2021b), which ensures the reliability and scientific validity of the SBAS-InSAR measurements results in this paper. The maximum settlement rate in the study area was −71.7 mm/a, and the maximum uplift rate was 50.8 mm/a. Settlement rates in other valuable areas mostly ranged from −15.4 to 15.6 mm/a (Figure 4).
The areas of cumulative deformation change in 2018 are concentrated at P11, P10, P9, and P8 (marked with small rectangles, pentagrams, ellipses, and large rectangles Frontiers in Earth Science | www.frontiersin.org November 2021 | Volume 9 | Article 729454 6 respectively in Figure 5, respectively). All four sites are in the subsidence zone and the changes are mainly reflected in an increase in the area of subsidence and an acceleration of subsidence rate. However, there are also some seasonal fluctuations, with an increase in subsidence from January to June and a slowdown in subsidence from June to December. The reason for this phenomenon may be due to the lifting of the ground surface caused by freeze-thaw action.
The change from 2018 to 2019 is most evident at P14, P12, P9, and P13, (marked with a triangle, large ellipses, small ellipses, and large circle on Figure 6, respectively). By January 1, 2019, subsidence at P14 continued to increase until December 27, 2019 and had expanded in area. Other regional changes were characterized by similar cyclical fluctuations as in 2018. However, between June 30, 2019 and September 22, 2019, small landslides occurred in P13, P9, and the south-central section through the  investigation of historical imageries and information, resulting in greater subsidence during this period, shown as a number of small dark blue dots in Figure 6. By the end of 2019, deformation in these three regions had stabilized.
Changes in 2019 are mainly in P14 (marked by a red triangle on Figure 7), with continued subsidence and an increase in the area of subsidence. Other areas show little change, mainly due to cyclical variations caused by freeze-thaw action. However, due to landslides during 2018 in P13, P9 and the south-central region, freeze-thaw uplift caused by loose accumulation is greater than that in the past and is shown on the map as a more concentrated darker area (marked with blue ovals and circles on Figure 7, respectively).
From Figures 5-7, we find that the cumulative deformation for each period was based on the starting date of December 25, 2017 and ended with the SAR image acquisition date. It is obvious that as time goes by, the cumulative deformation of the area near P14 (marked with a red triangle in Figures 6, 7)  According to velocity maps ( Figure 2) and cumulative deformation graphs (Figures 5-7), the regions with violent settlement or uplift are relatively scattered. Settlement regions are distributed in the northeast, southwest, and southeast sections of the study area, whereas uplift is obvious in the north and central areas. All these zones are isolated and sporadic.

Area Statistics and Quantitative Analysis
In order to understand the deformation characteristics corresponding to each land cover type, we took the land cover raster dataset as a baseline and converted it to vector data using the Raster to Polygon tool. The vector data were then used to calculate the area of each ground land cover and the area of decohered regions (blank area in Figures 3, 5-7). Statistical indicators were calculated for each land cover, including the proportion of decohered area, velocity range, mean, media number, and standard deviation (SD). The results are shown in Table 2.
After the calculation,the entire area of the entire research region is 1,590.78 km 2 . Bare ground occupied the largest area with 772.63 km 2 or 45.45%. This is followed by snow/ice and grass with 35.18% (559.47 km 2 ) and 18.82% (299.21 km 2 ), respectively. The other 0.55% were scrub/shrub (6.23 km 2 ), Frontiers in Earth Science | www.frontiersin.org November 2021 | Volume 9 | Article 729454 8 water (1.83 km 2 ), built-up areas (0.35 km 2 ), trees (0.25 km 2 ), and crops (0.06 km 2 ) respectively. From Table 2 and statistical calculations, the results indicate that the decohered phenomenon is serious within the entire study area, with the area of decorrelation being 984.11 km 2 and accounting for 61.87% of total area. The effective measurement area is only 38.13%, with an area of 606.67 km 2 . And bare ground, grass, and snow/ice are 398.38, 158.44, and 48.45 km 2 , respectively. The areas of other land cover types are all less than 0.6 km 2 .
Of the eight land cover types from Table 2, bare ground has the greatest range of velocity variation at −70.3-50.8 mm/a and crops have the smallest range of velocity variation at −9.0-6.9 mm/a. This result is due to the large area of bare ground and landslides occurring within. The area of crops is small and does not vary much. The settlement rate of grassland is significantly greater than the uplift rate, which is due to the fact that part of the area where the landslide is close to grassland, resulting in an increase in the settlement rate of grassland. Snow/ ice deformation rates ranged from −27.4-27.2 mm/a relatively symmetrically, with sedimentation rates slightly greater than uplift rates due to seasonal freeze-thaw. But the mean and median rates are positive, indicating the presence of localised widespread subsidence at snow/ice. Scrub/shrub rates varied mainly due to environmental succession caused by seasonal  freeze-thaw and snow/ice thawing, ranging from −16.1-24.9 mm/ a with an average rate of 1.7 mm/a. The built area is dominated by a subsidence trend with an average rate of −1.0 mm/a, mainly due to the thawing of the ground ice and permafrost, and more care should be taken to prevent the risk of collapse in the future. Water changes are somewhat fortuitous and the results are not very meaningful. The region of trees is completely decohered and cannot be effectively monitored for surface deformation.
The statistical analysis revealed that the study area was severely decohered. The severity of the decorrelation was found to be in trees, scrub/shrub, snow/ice, water, grass, crops, bare ground, and built area with the area decohered proportion of 100, 91.75, 91.34, 54.14, 47.05, 45.35, 44.87, and 38.7%, respectively.

Time-Series Variation Characteristics of Cumulative Deformation
To determine the variation characteristics of cumulative deformation, 14 observation points (OPs) (Table 3, Figures 3,   6A) were selected in areas with dramatic variation, and the cumulative deformation was extracted for each point. Figure 8 shows the results.
Due to the posture limitations of the SAR satellite when observing, the deformation information obtained through InSAR processing is usually not the real deformation of surface ground, instead the LOS direction. The effectiveness of the deformation results obtained is closely related to the terrain (slope and aspect) (see Figure 8), especially in mountainous areas. Specifically, with a few exceptions, the direction of surface movements fits expectations: the parallel movement, caused by widespread slope processes (e.g., general creep), was dominant in the middle section of slopes; for other portions, the rotational motion was prevalently caused by the alluvial accumulation or the combination dynamics of permafrost and the overlaid active layer (Chen et al., 2013). In addition, the occurrence of landslide will form a serious deformation area in a short time. Figure 8 shows a high agreement between the real surface deformation and the motion in the LOS direction, both at the top of the mountain and at the foot of the mountain, fore and back the slope (a, aa, b, bb, c 1 , cc, d, dd, e, and ee in Figure 8). Ground heave shows as uplift in LOS and surface settlement as subsidence in LOS. It is only when the slope faces satellite and slope gradient is greater than satellite incidence angle that the inconsistency is shown. In Figure 8, c2 marks a slope process (blue line), and an uplift signal is manifested in the LOS direction (red line). The incidence angle of Sentinel-1A images used in this study is 39.08°. For ascending orbit, the west aspect is the satellite-facing slope and the east aspect is the satellite-away slope (Schaufler et al., 2018). Combining the slope and aspect dataset (Figure 2), the deformation velocity map (Figure 3), and the slope-velocity and aspect-velocity statistics ( Figures 11A,B), we recognize that most uplift signals are located on satellite-facing slopes (west) and extensive subsidence signals are located on the east aspect. The signals (uplift and subsidence) are more likely caused by slope processes/land sliding.
To further understand the relationship between deformation rate and slope, slope orientation and vegetation, we counted the percentage of each element in the area with deformation rate greater than 15.0 mm/a. The results showed that, except for trees, the deformation rates of bare ground, grassland, snow, and wetland were all greater than 15.0 mm/a among the remaining seven types of ground features, with the percentages of 75.99, 17.91, 5.98, and 0.12%, respectively, and the deformation rates were concentrated between 16.0 and 21.5 mm/a. More statistical information is shown in Table 4.
It can be found from Table 4 that the areas with deformation velocity greater than 15.0 mm/a in the study area are mostly concentrated in the west, southwest, and northwest sides of the slope direction, accounting for 37.71, 29.54, and 18.83%, respectively. The west aspect is the satellite-facing slope and ascending orbits show a positive shift in backscatter for slopes facing west (Schaufler et al., 2018), so that the landslide accumulation phenomenon can be well observed in this area. Moreover, the slope of the hills sloping to the west is also relatively large, concentrated between 17.1 and 35.3°, and the average slope is all greater than 24.3°. In addition, the area is mostly covered with bare land, followed by grassland, with little snow/ice and scrub/shrub coverage, which is more likely to cause landslides and should be paid special attention.
The cumulative deformation of the OPs shows that points P1-P7 are situated in a settlement area, whereas P9-P14 are situated in an uplift area. Over time, the cumulative deformation becomes greater, involving both settlement and uplift. The slope of the cumulative deformation curves shows the velocity of deformation, and P3 and P8 have the maximum values of settlement and uplift, respectively, which agrees with the velocities in Table 3. Compared with settlement curves, uplift curves show more volatility, especially those for P13 and P14. However, some anomalies were observed in settlement curves, such as P2 on August 22, 2018 and September 16, 2019 and P7 on April 19, 2019. The sudden jumps that appeared in Figure 9 were very likely due to unwrapping errors or other noises.
By investigating historical imagery data, mainly via comparing Esri satellite maps and Google searching historical images, we found that the location where P7 is located originally had significant amounts of snow cover and ice, which has now largely thawed with global warming (Figure 10). The seasonal freezing and thawing of snow cover and ice, as well as the effects of solid winter precipitation, have induced small seasonal fluctuations changes in the cumulative deformation chart at point P7, but the overall trend is settlement.
According to close-up shots of the 14 OPs based on Esri satellite image maps (Figure 9), the settlement areas are mostly at the foot of mountains, which have an abundance of broken rocks or areas covered with snow and ice, especially on steep ridges where the landslides are prone to occur. Uplift regions occur mostly in the basins between ridges with relatively flat terrain, near the runoff points, so that changes in river discharge may modify the cumulative deformation process.

Relationship Between Deformation and Topographical Elements
Raster to Point tool was used to create the monitoring points (MPs) consistent with the deformation velocity raster grids within the area and extract the velocity, DEM, aspect, and slope of each point. Aspect and slope dataset for the study area ( Figure 2) were produced from DEM data. Decohered areas without a deformation rate were rejected, and a total of 242,778 points were generated. The statistical analysis addressed the distribution of effective MPs and the relationships between velocity and aspect, slope, and elevation. The results are shown in Figure 11. Most of the MPs are distributed in the areas toward the southeast, south, and southwest (112.5-247.5), with slopes ranging from 15°to 50°and altitudes ranging from 3,150 m.a.s.l. to 4,275 m.a.s.l. MPs are less prevalent in areas with slopes greater than 60°and altitudes higher than 4,275 m.  Frontiers in Earth Science | www.frontiersin.org November 2021 | Volume 9 | Article 729454 12 Therefore, most of the deformation occurs in areas with slopes of 15-50°and elevations of 3,150-4,275 m.a.s.l., especially in regions with southwest, west, and northwest aspects. Ground areas with severe deformation may induce some secondary geological disasters, such as landslides, collapses, and mudrock flows.

DISCUSSION
Due to the inherent limitations of InSAR measurements, certain areas lack deformation data due to decorrelation. In this study, significant incoherent regions covered with glaciers, snow cover, and ice were observed within the study area. In addition to the shortcomings of InSAR itself, this phenomenon is also related to the fact that degeneration of glaciers and snow cover not only occurs abruptly, but also in large amounts. Effective solutions for areas with missing values include field surveys, which can then be supplemented by spatial interpolation (Zhang, 2010). Alternatively, InSAR measurements can be processed from multiple SAR image sources and the results combined (Wang, 2019;Zhu et al., 2019). Sources of SAR images include products taken by different sensors, at different orientations, or in different time periods. In addition, reducing the coherence threshold is an effective method.
The deformation data obtained are of great significance, despite the large number of decoherent zones in the study area. The measurement results can provide basic advice and information for later continuous observations, and use of a sequential adjustment method  ensures that the long time series of SBAS-InSAR measurement data can be supplemented at a later stage. At the same time, current deformation values can provide a theoretical basis for environmental development by the local government. More importantly, the deformation data can provide a general picture over large areas and long time series of the distribution of mines and the degeneration of snow cover and ice in the area.
There are some differences regarding land cover between the close-up shots of OPs P5-P7 in Figure 9 and the land-cover data in Table 3; these may have been caused by the different dates of the images obtained. Land-cover data is derived from ESA Sentinel-2 imagery at 10 m resolution. And it is a composite of LULC predictions for 10 classes throughout the year in order to Frontiers in Earth Science | www.frontiersin.org November 2021 | Volume 9 | Article 729454 generate a representative snapshot of 2020. Whereas the close-up shots were based on the 2020 Esri image maps, which were the integrated product of various satellite images. The question of whether the deformation in the LOS direction reflects the true displacement of the ground surface has been addressed in many papers (Chen et al., 2013;Liu et al., 2019;Hu et al., 2021b;Ge et al., 2021). As the monitoring results in this paper are obtained based on a single data source and a single InSAR processing method, it is difficult to obtain the real 3-D ground surface deformation. However, combining the available relevant theories (Chen et al., 2013;Reinosch et al., 2020) and the geometric relationship between SAR satellites and monitoring sites (mainly concerning slope and aspect), we found that there are almost no anomalies (the deformation in the LOS direction does not coincide with the true surface displacement) in this study area and the deformation values of the monitoring points are highly consistent with the real surface displacement. The reliability of results in this paper is also supported by the high agreement between the InSAR deformation monitoring values and the actual monitoring values for the mine sites located within the study area (Du et al., 2021b). The correlation between slope and deformation velocity needs to be investigated in depth in the follow-up to propose slope threshold of slope when the anomalies occurred, which will allow us to better analyze the real surface displacements.

CONCLUSION
1) The maximum subsidence velocity in the study area was -70.7 mm/a in the western, northeastern, and southeastern regions, and the maximum uplift rate was 50.8 mm/a in the northern and central regions. The areas with intense subsidence and uplift were isolated and sporadic. 2) The maximum cumulative subsidence in the study area for the 3-year period from December 25, 2017 to January 2, 2021 was −211.58 mm, and the maximum cumulative uplift was 149.03 mm, which was consistent with the deformation velocity.
3) The study area was severely decohered, with a decorrelation area of 984.11 km 2 , accounting for 61.87% of the total area. The severity of the decorrelation is trees, scrub/shrub, snow/ ice, water, grass, crops, bare ground and built area,with decohered percentage of own area 100, 91.75, 91.34, 54.14, 47.05, 45.35, 44.87, and 38.7%, respectively. 4) There were correlations between deformation and slope, slope direction, and elevation. Specifically, deformation was obvious in areas with elevations between 3,150 and 4,275 m.a.s.l., slopes of 15°-50°, and southwest, west, and northwest aspects, which are disaster-prone regions.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.