Influence of Open-Pit Coal Mining on Ground Surface Deformation of Permafrost in the Muli Region in the Qinghai-Tibet Plateau, China

The Qinghai-Tibet Plateau (QTP) is the largest mid-to low latitude and high-altitude permafrost. Open-pit coal mining and other activities have caused serious damage to the alpine ecological environment and have accelerated the degradation of permafrost on the QTP. In this study, the influence of open-pit coal mining on the time series ground surface deformation of the permafrost in the Muli region of the QTP was analyzed from 19 January 2018 to 22 December 2020 based on Landsat, Gaofen, and Sentinel remote sensing data. The primary methods include human-computer interactive visual interpretation and the small baseline subsets interferometric synthetic aperture radar (SBAS-InSAR) method. The results showed that the spatial distribution of displacement velocity exhibits a considerably different pattern in the Muli region. Alpine meadow is the main land use/land cover (LULC) in the Muli region, and the surface displacement was mainly subsidence. The surface subsidence trend in alpine marsh meadows was obvious, with a subsidence displacement velocity of 10–30 mm/a. Under the influence of changes in temperature, the permafrost surface displacement was characteristics of regular thaw subsidence and freeze uplift. Surface deformation of the mining area is relatively severe, with maximum uplift displacement velocity of 74.31 mm/a and maximum subsidence displacement velocity of 167.51 mm/a. Open-pit coal mining had resulted in the destruction of 48.73 km2 of natural landscape in the Muli region. Mining development in the Muli region had increased the soil moisture of the alpine marsh meadow around the mining area, resulting in considerable cumulative displacement near the mining area and the acceleration of permafrost degradation.


Introduction
As an important and large part of the cryosphere, permafrost is significantly affected by global warming [1,2]. The Qinghai-Tibet Plateau (QTP) is the largest area of permafrost in the mid-to low latitude and with the highest elevation across the world [3]. With increasing global warming, the permafrost in the QTP has been degrading at a large scale [4][5][6][7][8][9][10]. The landscape structure and hydrothermal regimes of underlying permafrost have changed due to extensive human activities such as mining development, engineering construction, and urban expansion, which affect the soil water content and heat exchange

Study Region
The Muli region is located in the upper reaches of the Datong River, a tributary of the Yellow River, in the Qilian Mountains of the QTP. The geographical extent of the study region is 98 • 56 42"-99 • 42 27"E and 37 • 49 29"-38 • 13 26"N, and its altitude ranges from 3683 to 4844 m. Muli is home to four drainage basins from east to west, namely, Jiangcang, Azigou, Duosuo, and Cuokamori ( Figure 1). The study region has an alpine continental climate [43], with the mean annual air temperature of −4.3 to 3.6 • C and annual precipitation of 466.9-656.2 mm from 2017 to 2019. All of the river systems originate from the Datong River Basin, the major tributaries are Duosuo and Jingcang Rivers, and there are many small lakes in the area. The Muli region is an important water source of Datong River Basin and one of the most important wetland and ecological environment protection areas in China [42]. The coal mines in the Muli region are the largest open-pit coal mines on the QTP, and this region is rich in coal resources. At present, the two main mines are Juhugeng (JHG) mining area (mining began in 2003) and Jiangcang (JC) mining area (mining began in 2004) [44]. JHG mining area is located in an intermontane basin of Duosuo River Basin with an altitude of 4000-4200 m. The terrain is relatively transitional, and the ground surface is covered by dense alpine meadow. JC mining area is located in the intermontane basin of Jiangcang River Basin with an altitude of 2800-3900 m, and it is dominated by alpine marsh meadow [43].
Remote Sens. 2022, 14,2352 3 of 15 monitoring using SBAS-InSAR, the influence of mining development on the surrounding permafrost was evaluated. This study reported the degradation of permafrost under the background of global warming. Moreover, it may provide scientific basis for the protection, restoration, and sustainable development of alpine ecology in the Qilian Mountains.

Study Region
The Muli region is located in the upper reaches of the Datong River, a tributary of the Yellow River, in the Qilian Mountains of the QTP. The geographical extent of the study region is 98°56′42″-99°42′27″E and 37°49′29″-38°13′26″N, and its altitude ranges from 3683 to 4844 m. Muli is home to four drainage basins from east to west, namely, Jiangcang, Azigou, Duosuo, and Cuokamori ( Figure 1). The study region has an alpine continental climate [43], with the mean annual air temperature of −4.3 to 3.6 °C and annual precipitation of 466.9-656.2 mm from 2017 to 2019. All of the river systems originate from the Datong River Basin, the major tributaries are Duosuo and Jingcang Rivers, and there are many small lakes in the area. The Muli region is an important water source of Datong River Basin and one of the most important wetland and ecological environment protection areas in China [42]. The coal mines in the Muli region are the largest open-pit coal mines on the QTP, and this region is rich in coal resources. At present, the two main mines are Juhugeng (JHG) mining area (mining began in 2003) and Jiangcang (JC) mining area (mining began in 2004) [44]. JHG mining area is located in an intermontane basin of Duosuo River Basin with an altitude of 4000-4200 m. The terrain is relatively transitional, and the ground surface is covered by dense alpine meadow. JC mining area is located in the intermontane basin of Jiangcang River Basin with an altitude of 2800-3900 m, and it is dominated by alpine marsh meadow [43].  We collected 75 C-band (wavelength = 5.63 cm) single look complex (SLC) images from Sentinel-1A, which were acquired in the interferometric wide (IW) swath mode with VV polarization from 19 January 2018 to 22 December 2020. All scenes were acquired in an ascending orbit, and the incidence angles were approximately 36.83 • . The time interval passed between capturing the two images of the same area was 108 days from 26 April 2020 to 12 August 2020 and 60 days from 19 June 2019 to 18 August 2019. The frequency of revisiting other images was 12 or 24 days. To improve the registration accuracy of Sentinel-1A images and remove systematic errors caused by orbit errors, precise orbit ephemerides (POD) data corresponding to the time of Sentinel-1A image data were used. The positioning accuracy of the POD data was >5 cm. Sentinel-1A data were obtained from the Alaska Satellite Facility (ASF) data search Vertex (https://search.asf.alaska.edu (accessed on 20 October 2021)).

Landsat and Gaofen Data
To analyze the historical development process of the mining area, Landsat thematic mapper (TM) and operational land imager (OLI) remote sensing images with a spatial resolution of 30 m were used to extract the boundary of the mining area from 2003 to 2020. The images of Landsat in the Muli region were obtained from the United States Geological Survey (USGS) (https://earthexplorer.usgs.gov (accessed on 1 September 2021)). LULC in the Muli region in 2020 was extracted from Gaofen-1 remote sensing images with a spatial resolution of 2 m. The images of Gaofen were obtained from the China Centre for Resources Satellite Data and Application (http://www.cresda.com (accessed on 1 September 2021)).

Other Data
Data for the digital elevation model (DEM) were obtained from the Shuttle Radar Topography Mission (SRTM) (spatial resolution of 30 m) and the USGS (https://earthexplorer. usgs.gov (accessed on 20 October 2021)). DEM data were used to simulate the topographic phase and geocode displacement results.
The meteorological information consisting of daily temperature and precipitation obtained from the meteorological station in Muli Town. It was used to analyze the relationship between daily mean temperature and daily cumulative precipitation and time series displacement.

Human-Computer Interactive Visual Interpretation
Based on Landsat and Gaofen remote sensing images, in combination with the shape, spectrum, texture, and other features of land class, the mining development process from 2003 to 2020 and the LULC of the Muli region in 2020 were extracted based on humancomputer interactive visual interpretation. LULC was divided into natural and human activities. Natural land was divided into grassland (alpine meadow), water bodies (river, lake, and pit pond), and unused land (bottomland, bare land, and bare rock). Human activity-related land was divided into mining areas (mine building, mine, finished coal, slag, mine water bodies, and grassland restoration), roads, and settlements.

SBAS-InSAR
SBAS-InSAR is based on interferometric pairs of multiple master images. The deformation information of time series in the study area was recovered based on high coherence points [31]. SBAS-InSAR considers the temporal and perpendicular baseline distribution of images in the selection of interferometric pairs. It can overcome the aforementioned strong temporal decorrelation over permafrost terrain and calculate high-precision time series surface deformation [32].
SBAS-InSAR used m Sentinel-1A images from the same orbit in the same region at different times. According to the set temporal and perpendicular thresholds, n interferograms were generated. A relatively stable target point was selected as the reference pixel point to phase unwrapping of n multi-view differential interferograms. A pixel point (x, r) was indicated in the i differential interferometry phase diagram, and the unwrapped differential interferometry phase can be expressed using Equation (1): where t B and t A are the acquisition time of Sentinel-1A master and slave images, respectively; X and r are the pixel azimuth and range direction coordinates, respectively. We assumed that the surface displacement of radar line-of sight (LOS) in this area at reference time t 0 was equal to 0. d(t B ,x,r) and d(t A ,x,r) are the LOS cumulative displacement of the master and slave images relative to t 0 , respectively; λ is the central wavelength of Sentinel-1A; ∆ϕ i,arb (x,r) is the orbit error phase; ∆ϕ i,res (x,r) is the residual phase, including atmospheric delay, phase noise, high-pass deformation component, and so on; and ∆ϕ i,topo (x,r) is the terrain residual phase, and it is expressed using Equation (2): where ∆h is the elevation error; R i is the slope distance of Sentinel-1A master image; θ i is the incident angle of Sentinel-1A master image; and B ⊥ is the vertical baseline distance of the interferometric pairs. SBAS-InSAR was utilized to calculate the mean displacement velocity and time series displacement in the Muli region. The image from 19 May 2018 was selected as the master image. After image clipping, sub-pixel precision level registration was completed based on the master image. To avoid low robustness and accuracy of the final singular value decomposition (SVD) solution owing to the isolation of small baseline sets, the temporal baseline threshold was set to 120 days, and the perpendicular baseline threshold was set to 2% of the critical threshold. Interferometric pairs with low coherence or poor unwrapping were deleted and 249 groups were retained. To suppress the decoherence noise in the interferogram and improve the coherence of the interferogram, 1 × 4 multi-view time and Goldstein adaptive filtering were used. Delaunay minimum cost flow (Delaunay MCF) method was used to complete phase unwrapping. A pixel with coherence coefficient > 0.2 was selected to estimate the displacement velocity and residual topography. After estimating and removing atmospheric delay, the filter window in the perpendicular and temporal domains was set for 1200 and 365 days, respectively.

LULC Structure and Mining Development Process in the Muli Region
The LULC of the Muli region in 2020 interpreted by Gaofen-1 image is shown in Figure 2a; the alpine meadow area is the largest (904.55 km 2 ) in the natural class, accounting for 75.69% of the total area in the study region. The bare rock area is 170.80 km 2 , accounting for 14.29% of the total area, and it is mainly distributed in mountainous areas with high elevations. The Bare land, bottomland, and water bodies comprise small areas of 40.03 km 2 , 13.53 km 2 , and 12.07 km 2 , respectively ( Figure 2b). Mining area is the largest area affected by human activities. The area of natural surface damage due to coal mining is 48.73 km 2 , which comprises 14.22 km 2 mine area and 11.01 km 2 slag area. Mine building, finished coal, and mine water bodies have relatively small areas. In recent years, the ecological environment of the mining area has been recovered in the Muli region, with 18.84 km 2 mining area recovered to grassland. Areas occupied by roads and settlements are less, with an area of 4.71 km 2 and 0.64 km 2 , respectively ( Figure 2c).
Among the mining development processes from 2003 to 2020 (Figure 2d), the development of JHG and JC mining areas was mainly concentrated from 2003 to 2015, with development rates of 2.72 km 2 /a and 1.05 km 2 /a, respectively. JHG mining caused more damage to the natural surface. In 2015, the area of JHG mining was 23.61 km 2 more than that of JC mining. Under the grass planting ecological restoration and management of the ecological environment of the mining area after 2015, some mines have been backfilled. This is also attributed to the reduced rate of mining.

Surface Displacement Velocity and Time Series Displacement in the Muli Region
Deformation area was 707.93 km 2 , accounting for 59.28% of the total area of the Muli region. Considerable incoherence was observed in the alpine marsh meadow in the Among the mining development processes from 2003 to 2020 (Figure 2d), the development of JHG and JC mining areas was mainly concentrated from 2003 to 2015, with development rates of 2.72 km 2 /a and 1.05 km 2 /a, respectively. JHG mining caused more damage to the natural surface. In 2015, the area of JHG mining was 23.61 km 2 more than that of JC mining. Under the grass planting ecological restoration and management of the ecological environment of the mining area after 2015, some mines have been backfilled. This is also attributed to the reduced rate of mining.

Surface Displacement Velocity and Time Series Displacement in the Muli Region
Deformation area was 707.93 km 2 , accounting for 59.28% of the total area of the Muli region. Considerable incoherence was observed in the alpine marsh meadow in the intermontane basin and in some mining areas. The surface displacement in the Muli region was mainly attributed to subsidence, accounting for 76.36% of the total deformation area. There were obvious differences in surface deformation between natural and human activities ( Figure 3). intermontane basin and in some mining areas. The surface displacement in the Muli region was mainly attributed to subsidence, accounting for 76.36% of the total deformation area. There were obvious differences in surface deformation between natural and human activities ( Figure 3).

Surface Displacement Velocity and Time Series Displacement of Natural Type
The displacement velocity of natural type was between −30 and 10 mm/a. The main surface displacement was subsidence, accounting for 75.73% of the natural deformation area. The subsidence of alpine marsh meadow in the intermountain basin was the most significant, and the mean surface displacement velocity was between −30 and −10 mm/a.

Surface Displacement Velocity and Time Series Displacement of Natural Type
The displacement velocity of natural type was between −30 and 10 mm/a. The main surface displacement was subsidence, accounting for 75.73% of the natural deformation area. The subsidence of alpine marsh meadow in the intermountain basin was the most significant, and the mean surface displacement velocity was between −30 and −10 mm/a. The surface of bare rock in the middle and high mountains in the southern area of Muli was relatively stable, and the mean surface displacement velocity was between −10 and 10 mm/a. In Jiangcang River Basin, the alpine meadow distributed along the river showed a tendency of uplifting, and it was surrounded by the alpine marsh meadow. According to the characteristic point time series displacement of three typical regions, namely, the stable region of mountain bare rock (A, mean displacement velocity −0.87 mm/a), the subsidence region of alpine marsh meadow (B, mean displacement velocity −12.40 mm/a), and the uplift region of alpine meadow (C, mean displacement velocity 2.85 mm/a) (Figure 4), the surface displacement in bare rock region fluctuated around zero, and there was no obvious periodic deformation trend. The time series displacement of alpine marsh meadow presented a trend of periodic and linear slow subsidence. The time series displacement of alpine meadow showed a trend of periodic and linear slow uplift. Both alpine marsh meadow and alpine meadow are characteristics of obvious freeze-thaw cycles on permafrost surface. Combined with the temperature and precipitation data of Muli Town meteorological station, the annual periodic thaw process in the permafrost region in 2018-2020 started from the end of April or early May to the end of September or early October, and the surface subsidence reached the maximum value. The freeze process started in early October and continued until early May of the following year. The change of permafrost displacement was observed to lag behind climate change. The maximum temperature and precipitation appeared at the end of July or early August, and the maximum time lag of surface thaw subsidence was approximately 60 days. The minimum temperature and precipitation appeared in the end of January or early February, and the maximum time lag of surface freeze uplift was approximately 90 days.
The surface of bare rock in the middle and high mountains in the southern area of Muli was relatively stable, and the mean surface displacement velocity was between −10 and 10 mm/a. In Jiangcang River Basin, the alpine meadow distributed along the river showed a tendency of uplifting, and it was surrounded by the alpine marsh meadow.
According to the characteristic point time series displacement of three typical regions, namely, the stable region of mountain bare rock (A, mean displacement velocity −0.87 mm/a), the subsidence region of alpine marsh meadow (B, mean displacement velocity −12.40 mm/a), and the uplift region of alpine meadow (C, mean displacement velocity 2.85 mm/a) (Figure 4), the surface displacement in bare rock region fluctuated around zero, and there was no obvious periodic deformation trend. The time series displacement of alpine marsh meadow presented a trend of periodic and linear slow subsidence. The time series displacement of alpine meadow showed a trend of periodic and linear slow uplift. Both alpine marsh meadow and alpine meadow are characteristics of obvious freeze-thaw cycles on permafrost surface. Combined with the temperature and precipitation data of Muli Town meteorological station, the annual periodic thaw process in the permafrost region in 2018-2020 started from the end of April or early May to the end of September or early October, and the surface subsidence reached the maximum value. The freeze process started in early October and continued until early May of the following year. The change of permafrost displacement was observed to lag behind climate change. The maximum temperature and precipitation appeared at the end of July or early August, and the maximum time lag of surface thaw subsidence was approximately 60 days. The minimum temperature and precipitation appeared in the end of January or early February, and the maximum time lag of surface freeze uplift was approximately 90 days.

Surface Displacement Velocity and Time Series Displacement of Mining Area
Due to frequent human activities, most areas in JHG mining area were incoherent. The displacement in the center of the mine was mainly uplifted (Figure 3). The point with the maximum uplifted displacement was located at the center of the mine in coal mine No. 3, with a mean displacement velocity of 74.31 mm/a. The surface uplift displacement velocity decreased from the center of the mine to the surrounding area. There was an obvious subsidence displacement at the mine edges of coal mines No. 4, No. 5, and No. 6. The slag of all mine showed obvious subsidence displacement. The maximum subsidence point was located in coal mine No. 5 with a mean displacement velocity of −167.51 mm/a. The displacement of mine building and finished coal was mainly subsidence, but its displacement was smaller than that of slag.
The incoherence area of JC mining area was mainly located in the slag of coal mines No. 1 and No. 3, and an obvious ecological restoration effect was observed in these mines. There was no obvious surface uplift displacement in all mine areas, and a small area at the edge of the mine was slightly uplifted. The maximum uplift displacement was

Surface Displacement Velocity and Time Series Displacement of Mining Area
Due to frequent human activities, most areas in JHG mining area were incoherent. The displacement in the center of the mine was mainly uplifted (Figure 3). The point with the maximum uplifted displacement was located at the center of the mine in coal mine No. 3, with a mean displacement velocity of 74.31 mm/a. The surface uplift displacement velocity decreased from the center of the mine to the surrounding area. There was an obvious subsidence displacement at the mine edges of coal mines No. 4, No. 5, and No. 6. The slag of all mine showed obvious subsidence displacement. The maximum subsidence point was located in coal mine No. 5 with a mean displacement velocity of −167.51 mm/a. The displacement of mine building and finished coal was mainly subsidence, but its displacement was smaller than that of slag.
The incoherence area of JC mining area was mainly located in the slag of coal mines No. 1 and No. 3, and an obvious ecological restoration effect was observed in these mines. There was no obvious surface uplift displacement in all mine areas, and a small area at the edge of the mine was slightly uplifted. The maximum uplift displacement was located in the slag area in the north of coal mine No. 1, with a mean displacement velocity of 56.70 mm/a. The other slags were dominated by subsidence displacement. The maximum subsidence displacement velocity was −118.83 mm/a in the southern slag area of coal mine No. 1 (Figure 3).
The typical points of mine, slag, mine building, and finished coal area are shown in Figure 5. The times series displacement at the typical points in the mining area exhibited no obvious periodic change rules of freeze uplift and thaw subsidence. However, in winter, when the temperature was the lowest, the typical points of the BC1-BC4 showed that the subsidence displacement reached the maximum value and then started decreasing slowly. All of the displacements are in the state of subsidence. An increase in fluctuation was observed, and the maximum displace was −80.4 mm. The displacement of the typical points of Mi1-Mi4 showed uplift and linearly increased, and the maximum displacement of uplift was 23.9-106.5 mm. The displacement of the typical points of Sl1-Sl4 showed subsidence and linearly increase. The displacement was higher than the displacement of other typical points, and the maximum displacement was −157.0 to −250.3 mm. located in the slag area in the north of coal mine No. 1, with a mean displacement velocity of 56.70 mm/a. The other slags were dominated by subsidence displacement. The maximum subsidence displacement velocity was −118.83 mm/a in the southern slag area of coal mine No. 1 (Figure 3).
The typical points of mine, slag, mine building, and finished coal area are shown in Figure 5. The times series displacement at the typical points in the mining area exhibited no obvious periodic change rules of freeze uplift and thaw subsidence. However, in winter, when the temperature was the lowest, the typical points of the BC1-BC4 showed that the subsidence displacement reached the maximum value and then started decreasing slowly. All of the displacements are in the state of subsidence. An increase in fluctuation was observed, and the maximum displace was −80.4 mm. The displacement of the typical points of Mi1-Mi4 showed uplift and linearly increased, and the maximum displacement of uplift was 23.9-106.5 mm. The displacement of the typical points of Sl1-Sl4 showed subsidence and linearly increase. The displacement was higher than the displacement of other typical points, and the maximum displacement was −157.0 to −250.3 mm.

Mean Cumulative Displacement in Permafrost Region Outside the Mining Area
Further, the influence of open-pit coal mining on permafrost surface deformation in the Muli region was analyzed and evaluated. Based on the boundaries of JHG and JC mining areas, a total of 50 buffer gradient regions (within 5 km outside the mining area) were set at 100 m intervals. The mean value of cumulative displacement was calculated in the buffer gradient region from January 2018 to December 2020. The mean value of cumulative displacement varying with buffer distance was analyzed.
The mean cumulative displacement fluctuated between −23.6 mm (0-100 m) and −19.3 mm (300-400 m), with an average value of −21.4 mm in the 0-700 m buffer zones of JHG mining area. When the buffer distance exceeded 700 m, the mean cumulative displacement decreased with increasing elevation from −21.4 mm (600-700 m) to −8.1 mm (2100-2200 m), and the changing rate was 0.8 mm per 100 m. When the buffer distance exceeded 2200 m, the mean cumulative displacement changed slightly and fluctuated between −14.0 mm and −7.5 mm (Figure 6a). The mean cumulative displacement in the 0-2200 m buffer zones of JC mining area fluctuated between −20.6 mm (1600-1700 m) and −7.5 mm (500-600 m), with an average value of −15.0 mm. When the buffer exceeded 2200 m, the mean cumulative displacement decreased with increasing altitude from −17.3 mm (2100-2200 m) to 0.1 mm (3000-3100 m), with an increase rate of 1.8 mm per 100 m. When the buffer exceeded 3100 m, the mean cumulative displacement varied only slightly and fluctuated from −5.9 mm to 0.1 mm (Figure 6b). Field investigations revealed that the JHG and JC mining areas are located in the valley basin region, and the surrounding area is mainly alpine marsh meadow. The development of open-pit coal mines transformed and changed the runoff path of surface water. The water discharged from the slag mountain now flows directly to the alpine marsh meadow nearby. Furthermore, waterretaining walls were set up around the mine, which led to accumulation of water around the mine into lakes and meadows, increasing soil moisture around the mine (Figure 7). This phenomenon led to a relatively high mean cumulative displacement around the mining area and accelerated the degradation rate of permafrost.
Remote Sens. 2022, 14, 2352 10 of 15 Further, the influence of open-pit coal mining on permafrost surface deformation in the Muli region was analyzed and evaluated. Based on the boundaries of JHG and JC mining areas, a total of 50 buffer gradient regions (within 5 km outside the mining area) were set at 100 m intervals. The mean value of cumulative displacement was calculated in the buffer gradient region from January 2018 to December 2020. The mean value of cumulative displacement varying with buffer distance was analyzed.
The mean cumulative displacement fluctuated between −23.6 mm (0-100 m) and −19.3 mm (300-400 m), with an average value of −21.4 mm in the 0-700 m buffer zones of JHG mining area. When the buffer distance exceeded 700 m, the mean cumulative displacement decreased with increasing elevation from −21.4 mm (600-700 m) to −8.1 mm (2100-2200 m), and the changing rate was 0.8 mm per 100 m. When the buffer distance exceeded 2200 m, the mean cumulative displacement changed slightly and fluctuated between −14.0 mm and −7.5 mm (Figure 6a). The mean cumulative displacement in the 0-2200 m buffer zones of JC mining area fluctuated between −20.6 mm (1600-1700 m) and −7.5 mm (500-600 m), with an average value of −15.0 mm. When the buffer exceeded 2200 m, the mean cumulative displacement decreased with increasing altitude from −17.3 mm (2100-2200 m) to 0.1 mm (3000-3100 m), with an increase rate of 1.8 mm per 100 m. When the buffer exceeded 3100 m, the mean cumulative displacement varied only slightly and fluctuated from −5.9 mm to 0.1 mm (Figure 6b). Field investigations revealed that the JHG and JC mining areas are located in the valley basin region, and the surrounding area is mainly alpine marsh meadow. The development of open-pit coal mines transformed and changed the runoff path of surface water. The water discharged from the slag mountain now flows directly to the alpine marsh meadow nearby. Furthermore, water-retaining walls were set up around the mine, which led to accumulation of water around the mine into lakes and meadows, increasing soil moisture around the mine (Figure 7). This phenomenon led to a relatively high mean cumulative displacement around the mining area and accelerated the degradation rate of permafrost.

Precision Analysis of Surface Deformation Monitoring
The surface deformation obtained through SBAS-InSAR technology was improved by increasing the temporal and perpendicular baseline thresholds. This approach can effectively reduce the incoherence and high phase error of images during calculation, and thus, the monitoring accuracy of surface deformation can reach the millimetre scale [31,46]. This method can be effectively applied to the analysis of mining areas and

Precision Analysis of Surface Deformation Monitoring
The surface deformation obtained through SBAS-InSAR technology was improved by increasing the temporal and perpendicular baseline thresholds. This approach can effectively reduce the incoherence and high phase error of images during calculation, and thus, the monitoring accuracy of surface deformation can reach the millimetre scale [31,46]. This method can be effectively applied to the analysis of mining areas and permafrost surface deformation. From January 2018 to December 2020, the surface displacement velocity in the Muli region in the QTP ranged from −30 mm/a to 10 mm/a, and its displacement velocity was similar to that of Yeniugou (−30 mm/a to 10 mm/a, October 2014 to October 2016) [38], of Beiluhe (−20 mm/a to 20 mm/a, June 2007 to December, 2010) [47], and of Qinghai-Tibet Railway (−30 mm/a to 15 mm/a, March 2017 to March 2020) [48] in the permafrost regions of the QTP.
Surface deformation in permafrost region is complexly affected by vegetation, topography, soil texture, ground ice content, unfrozen water content, and human activities [40,49,50]. The displacement velocity and time series of displacement of permafrost varied with different LULC [40]. In alpine meadow area with dense vegetation cover, permafrost is stable and surface displacement is small [51]. In plain areas where many rivers and lakes dominate in alpine marsh meadow, the soil moisture content is relatively high, and the corresponding surface freeze uplift and thaw subsidence are relatively large [49]. In this study, the mean displacement velocity of the alpine marsh meadow in the Muli region was more than −10 mm/a, and the maximum value was −30 mm/a. Bare rock is mostly located in mountainous areas, whose complex topography and large slope are not conducive to water retention. The unfrozen water content of soil is low, and the surface deformation is small [50]. The mean displacement velocity of bare rock in the Muli region was less than −1 mm/a, and the surface displacement fluctuated within a small range during the year.

Characteristics of Surface Deformation in Open-Pit Mining Area and Its Influence on Surrounding Permafrost
Human activities exacerbate the surface deformation of permafrost; in particular, open-pit coal mining has more direct and severe influences than linear engineering such as highway and railway [18,43,52,53]. Open-pit coal mining in the Muli region has caused considerable surface deformation as well as the formation of many unstable surfaces. Surface subsidence after mining also disturbs the hydrothermal balance of the surface, which causes the dry-out of surrounding alpine meadow. Thus, the surface thermal radiation changes, and the degradation of permafrost around the mining area is aggravated [54]. In the Muli region, surface deformation in each coal mine was strong, and it was mainly distributed in mine and slag areas. In JC mining area, the slag mountains were mainly characterized by subsidence displacement, and the displacement velocity reached −118.83 mm/a of coal mine No. 1. In JHG mining area, the subsidence displacements were the most significant in slag mountain area, and the maximum displacement velocity was located at the top of the slag mountain, which were −164.87 mm/a and −167.51 mm/a and appeared in coal mines No. 4 and No. 5, respectively. Combined with the changes in the surface landscape in the slag mountain area, backfilling of mines and grass restoration occurred in the slag mountain area. The mining and backfilling of the mine changed the hydrothermal condition of the original natural stratum and the physical properties of the soil [55]. Under the combined influence of the compaction of the original active layer, the compaction of the slag mountain itself, and the insufficient compaction of the surface after backfilling, along with the creep deformation of the permafrost, there was uneven subsidence of the slag mountain and the mine [56,57]. Therefore, mine backfilling, slag hill greening, and natural slope deformation are important factors that lead to large surface deformation of slag hills. The uplift of the bottom of the mine is mainly attributed to the elastic rebound of the overlying soil after excavation. This phenomenon is slow, and it commonly occurs when there are considerable changes in the surface [58].
The development of mines in the Muli region has directly changed the natural geomorphic landscape, the structure, and the freeze-thaw balance of the original permafrost [56]. Mining causes the ground ice to melt, which leads to significant degradation or even disappearance of permafrost in the mining area. The permafrost table continues to deepen under the influence of global warming, and the degradation of permafrost near the table leads to continuous surface subsidence displacement [4,59]. In addition, the destruction of alpine meadow, the change of river system, slope collapse phenomenon, and other related factors, exacerbate the degradation of permafrost around the mining area [60].

Conclusions
The Muli region is located in the permafrost region of the QTP. Influenced by global climate change and extensive human activities, the spatial distribution of surface deformation varies significantly. Based on Landsat, Gaofen, and Sentinel satellite remote sensing data, this study analyzed the impact of open-pit coal mine development on surface deformation in the permafrost. Open-pit coal mining in the Muli region has resulted in significant destruction of 48.73 km 2 of natural land surface. Based on the analysis of deformation results calculated by SBAS-InSAR, it is determined that the surface deformation of alpine marsh meadow in intermontane basin area is the largest, with the largest contribution of subsidence displacement and the mean displacement velocity between −30 and −10 mm/a. Due to climate change, there are obvious periodic characteristics of thaw settlement and frost heave. The degree of surface deformation is the least in mountain bare rock area. Compared with natural slow deformation, the surface deformation of the mining area is severe, and it is mainly distributed in slag and mining areas. The maximum displacements of typical points in the mining area and slags were 23.9 to 106.5 mm and −157.0 to −250.3 mm, respectively. Open-pit coal mining has significantly exacerbated the surface deformation and accelerated the degradation of permafrost in some surrounding areas. C-band SAR has significant incoherence in areas with high vegetation cover and high soil moisture content, and the deformation in areas with severe outburst in the mining area cannot be monitored. Therefore, there is absence of deformation points in 40.72% of the entire study area. This increased the error and uncertainty of regional deformation results. The deformation calculation methods and data sources will be improved in subsequent studies to obtain deformation results with higher accuracy.
Author Contributions: H.W., Y.Q. and J.Z (Juan Zhang). designed the experiments and wrote the paper; J.Z. (Jinlong Zhang), R.Y., J.G. and S.Z. processed the satellite images and performed the experiments; and J.W. and D.L. reviewed the paper and gave constructive advice. All authors contributed to the redaction of the paper. All authors have read and agreed to the published version of the manuscript.