Remote Sensing Detection of Vegetation and Landform Damages by Coal Mining on the Tibetan Plateau

Abstract

In order to satisfy the needs of constant economic growth, the pressure to exploit natural resources has been increasing rapidly in China. Particularly with the implementation of the National Western Development Strategies since 1999, more and more mining activities and related infrastructure constructions have been conducted on the Tibetan Plateau (TP). Mining activities are known to have substantial impacts on plant dynamics and hence the water and energy cycles. Identifying mining activities and quantifying their effects on vegetation cover are critical to the monitoring and protection of the pristine TP environment. Thus, this study aims to develop an automated approach that detects the timing of initial mining development and assess the spatial distribution of mining-ruined vegetation. The Breaks for Additive Seasonal and Trend (BFAST) algorithm was used to decompose the signal in the normalized difference vegetation index (NDVI) time series derived from high-frequency MODIS images, and to detect abrupt changes of surface vegetation. Results show that the BFAST algorithm is able to effectively identify abrupt changes in vegetation cover as a result of open-mining development on the studied alpine grassland. The testing study in Muli Town of Qinghai Province shows that the mining development began in 2003 and massive destructions of vegetation cover followed between 2008 and 2012. The integrated use of Landsat imagery and multi-temporal DEMs further reveals detailed areal and volumetric changes in the mining site. This study demonstrates the potential of applying multi-mission satellite datasets to assess large-scale environmental influences from mining development, and will be beneficial to environmental conservation and sustainable use of natural resources in remote regions.

1. Introduction

The Tibetan Plateau (TP, hereafter) is well known as the Earth’s “Third Pole”, with an average altitude exceeding 4000 m [1,2]. Due to its broad geographic extent and high altitude, the plateau plays an important role in regulating climate, ecosystem and water budgets both in East Asia and globally [3,4]. The vegetation cover is a fundamental tie that links the carbon uptake of ecosystem processes, energy exchange, and water cycles on the TP [5,6]. Alpine grassland is the primary vegetation cover on the TP, accounting for ~61% of the plateau area [7,8,9]. Evidence from field observations, remote sensing measurements, and ecosystem model simulations consistently suggest that the alpine vegetation growth is particularly sensitive and vulnerable to the rapidly changing climate [10,11,12,13]. In recent decades, substantial climate changes have been observed on the TP [4,14,15,16], which acted on a diverse range of land surface types and have resulted in complex responses from the vegetation, including enhanced greenness, productivity [8,13,17], and respiration [18,19], loss of species richness [20], and advanced spring phenology [9,21,22].

Beyond the overall vegetation recovery, there were also serious grassland degradations at local scales during recent years [13,23,24,25,26]. These could be largely attributed to a combination of various factors, including climate variability, overgrazing, unscientific livestock management, herbal medicine collection, tourisms, and bioturbation from small mammals [5,23,24,25,26,27]. Most existing studies emphasized large-scale climate variability and overgrazing activities associated with the vegetation growth. The influences of mining development on the TP grasslands and related environmental impacts remain poorly addressed. The mining activities not only destroy the original vegetation cover, but also disturb the freezing/thawing cycle of the permafrost’s active layers, leading to possible deduction in rangeland productivity. Thus, it is necessary to develop an automated approach that can detect and assess mining-induced vegetation destructions on the remote TP.

Remote sensing-based normalized difference vegetation index (NDVI) has been widely used to monitor phenological dynamics and land cover conversions [28]. However, as many disturbance sources (e.g., climate, grazing, volcano events) often lead to similar surface spectral responses [29,30], identifying the impacts of mining activities on the vegetation cover remains a challenging task. An effective approach requires frequent observations with long-term temporal coverage in order to differentiate natural changes from those associated with human activities. Changes within long-term NDVI time series can be partitioned into three major signals: seasonal, gradual, and abrupt changes [31]. While seasonality depicts the natural variability within the vegetation’s phenological cycle [32], the latter two reveal possible changes in trend due to external forces. A gradual change is usually caused by slowly acting environmental processes, such as climate change or land degradation. Over time, the gradual change may stall or reverse [33,34], manifested as an abrupt change in magnitude and/or direction which is defined as a trend break [31]. In this study, we expect that the initiation of a mining activity could cause such an abrupt change in the NDVI time series. The Breaks for Additive Seasonal and Trend (BFAST) tool has been developed to distinguish among different signals in the NDVI time series [35]. Therefore, we here test the BFAST performance on decomposing NDVI time series derived from high-temporal-resolution Moderate Resolution Imaging Spectroradiometer (MODIS) imagery, and aim to develop an automated approach of detecting the timing of mining activities and the spatial distributions of mining-ruined vegetation covers on the TP. In addition, the medium-resolution Landsat images and SRTM/ALOS digital elevation models (DEMs) were used to map the fine-detailed mining extents and to quantify the excavated soil and mineral volumes.

2. Study Area

The study area is located in Muli Town in the northwestern part of the Haixi Mogolian and Tibetan Autonomous Prefecture, Qinghai Province, China (Figure 1). Geographically, this region is south of the Qilian Mountains in the northern TP and northwest of Lake Qinghai, the largest lake in China. The study area covers approximately 400 km² with altitudes ranging from 3950 m to 4150 m. The mineral reserve in Muli Town is estimated over 3.5 billion tons and may account for more than 67% of the entire mineral reserve in Qinghai Province [36,37]. Climate of this region belongs to the plateau continental regime that is cold and dry throughout the year, with low annual precipitation from 280 mm to 770 mm, very low air temperature, large diurnal temperature, long sunshine time and strong solar radiation [38]. The vegetation cover is mainly composed of alpine meadows and meadow steppe communities, which are dominated by Gramineae and Asteraceae, and accompanied by Ophiocordyceps sinensis, Armillaria luteo-virens, Nostoc commune and other probiotic species [39,40]. This region is also the headwater area for Qinghai Lake and Datong River. Therefore, it is considered as one of the most critical regions for wetland and ecological conservations in China, especially Qinghai Province.

3. Data and Methodology

3.1. Data and Processing

The MODIS instrument onboard the Terra and Aqua satellites provide globally-covered daily multispectral images that are ideal for monitoring large-scale phenological dynamics [41,42]. In this study, spectral reflectance data from MODIS Terra were used since they cover a longer time period than that from MODIS Aqua. Specifically, we used the 8-day composite, 500-m resolution Surface Reflectance Product MOD09A1 [43], downloaded from the MODIS/VIRS Subset website (https://modis.ornl.gov/cgi-bin/MODIS/global/subset.pl). This product contains reflectance observations at both red and near-infrared (NIR) channels which are required for calculating normalized differential vegetation index (NDVI), the index we employed for detecting vegetation changes. This level-3 data were atmospherically corrected using lower-level data, and composited from consecutive daily observations within each 8-day interval. This temporal composition reduced cloud contaminations, and thus are desirable for continuous time series analysis. We confined our study period from 27 July 2000 to 14 March 2018, in order to start from the vegetation’s peak growing season while ensuring the longest possible monitoring duration. This results in a total of 812 MODIS scenes. Each scene fully contains the study area with a spatial range of 30 km × 30 km and 9760 pixels (160 columns × 61 lines).

To delineate the extent of the mining-ruined area at higher spatial resolutions, Landsat 5 Thematic Mapper (TM), Landsat 7 Enhanced Thematic Mapper Plus (ETM+), and Landsat 8 Operational Land Imager (OLI) Level-1 terrain corrected data at 30 m resolution have been used. Landsat images acquired from 2003 to 2016 were downloaded from the USGS Global Visualization Viewer (GloVis) website (https://glovis.usgs.gov). In addition, digital elevation models (DEMs) were used to reveal the topographic changes due to mining activities. The Shuttle Radar Topography Mission DEM at 1″ (~30-m) resolution (SRTM-1), which collected the land elevations using space-borne radar measurements in February 2000, was utilized for representing the topographic condition before mining activities. To understand the topography condition after mining, we used the ALOS World 3D–30 m resolution (AW3D30) DEM [44] which was derived from stereo images acquired during January 2006 to April 2011 by the Panchromatic Remote-sensing Instrument for Stereo Mapping (PRISM).

3.2. Methods

Figure 2 illustrates the schematic framework of automated identification and extraction of mining-ruined vegetation areas. The procedure includes four modules: (1) preparing pixel-by-pixel NDVI time series for the BFAST input, including calculating NDVI for each pixel, and temporally filtering and filling NDVI values, (2) decomposing NDVI time series by BFAST and identifying breakpoint pixels, (3) mapping mining-ruined areas from Landsat imagery in the identified mining site, and (4) calculating the excavated mass in the mining area.

3.2.1. Preprocessing of NDVI Time Series

The NDVI index is widely used to indicate vegetation health conditions [28,34,41]. It is here computed by the reflectance of red band (${\rho }_{RED}$, band 1) and the near-infrared band (NIR) (${\rho }_{NIR}$, band 2) from the MOD09 product:

$$\mathrm{NDVI}=\frac{{\rho }_{NIR}-{\rho }_{RED}}{{\rho }_{NIR}+{\rho }_{RED}}$$

(1)

Due to instrument malfunction and extreme weather conditions (e.g., cloud cover or snowfall), the NDVI time series may be incomplete and can contain extreme values beyond the normal range of [−1,1]. We defined extreme values larger than 0.95 and smaller than −0.05 (mostly during winter seasons) as noises and replaced them with the value 0. The reason for selecting value 0 for this filling is to avoid introducing high or low values that may be mistaken as “abrupt changes” in the time series analysis. After filling the abnormal values, the complete NDVI time series at each pixel was output to text format that was later used as the input of the BFAST module. The calculation, filtering and output of NDVI were done with scripts in ENVI/IDL and Python/GDAL environments.

All the valid pixels in a single scene contributing to time-series were available for BFAST-based time series analysis.

3.2.2. Identification of Mining-Ruined Areas Based on BFAST-Detected NDVI Time Series

BFAST decomposes the NDVI time series into three components including trend, seasonal, and remainder with methods of detecting time-series changes [35]. It can be used to detect and characterize abrupt changes in trend and seasonal components to analyze many kinds of indexes or parameters, such as those from hydrology, climatology, and econometrics. In this study, BFAST package in Version 1.5.7 (in R tool) [31,35] was used to perform time series decomposition of the NDVI time series. Our purpose is to detect the abrupt changes that are associated with coal mining activities. The result of such activities is usually a land cover change from vegetation to barren land.

To effectively detect such changes, we focused on determining three key parameters in the BFAST application. The first one is the setting of ‘break’, which defines the maximum number of breakpoints (abrupt changes) that can be found in the trend component. We set this parameter to 1, as we intended to detect the most prominent change in each NDVI time series, which is most likely associated with the initiation of coal mining. The second parameter is the setting or ‘h’ value, which defines the minimal segment size that will be separated by the breakpoint. This parameter was set to 0.1, meaning that the minimal segment size is about two years. The third parameter is to select the seasonal model, which the options of dummy, harmonic and none. The harmonic seasonal model was used here, as it is less sensitive to inherent noise like clouds and atmospheric scatter than the dummy seasonal model [34]. The default setting was used for other parameters, such as max.iter which means the maximum iterative number (default value 2) in the trend analysis, level—the threshold value (default value 0.05) for the generalized fluctuation tests, hpc—the high performance computing support, and a type argument to empirical fluctuation processes.

After detecting the breakpoint in each time series, the next step was to determine which kind of breakpoint probably represents the coal mining activities. By analyzing different patterns of NDVI time series associated with different land surface status or processes, we used three conditions to retrieve the ‘true’ breakpoint:

(1)

There is one breakpoint detected in the time series. Assuming that one pixel is located in the new mining area, the NDVI time series of the pixel will show structural decreases after the mining activity. Based on the BFAST algorithm, one breakpoint should be detected if the mining-ruined areas have not yet recovered ecologically.

(2)

There is a decrease in NDVI magnitude in the trend components after the breakpoint. As the mining destruction of the land surface is an irreversible process in a short term and the value of NDVI after mining is generally lower than that in the beginning.

$$A_{starting}-A_{ending}>0$$

(2)

where $A_{starting}$ represents the last value in the trend component of the first segment, and $A_{ending}$ denotes the first value in the trend component of the second segment.

(3)

The magnitude of the post-breakpoint NDVI decrease is larger than the threshold 0.1. We used the difference between the mean NDVI value in the first segment $A_{starting-mean}$ and the mean NDVI value in the last segment $A_{ending-mean}$ to measure the level of NDVI decrease before and after the breakpoints. With a large NDVI decrease, we assume the change is probably associated with coal mining. The optimum threshold value for the change magnitude was set to 0.1 based on our tests on the threshold sensitivity, which will be discussed in the next section.

$$A_{starting-mean}-A_{ending-mean}>0.1$$

(3)

Based on the above three conditions, pixels that represent a change from vegetation cover to mining area during the study period were identified. Such examples are shown at the 4880th pixel and the 4083th pixel in Figure 3a,b. There are obvious trend breakpoints in both trend components. Meanwhile, the absolute value of the difference between two average values ($A_{starting-mean}$ and $A_{ending-mean}$) was greater than or equal to 0.1. The two examples have been validated in Google Earth history images.

By contrast, the 3743th pixel and 3733th pixel do not show such cover changes. Their time series and BFAST analysis results are illustrated in Figure 3c,d. For Pixel 3743, there is no breakpoint detected. For Pixel 3733, there is a positive trend in the trend component, which means the NDVI value has an abrupt increase probably because of some natural or artificial activity. The two pixels are located in areas with absence of mining activities, as confirmed from Google Earth images.

3.2.3. Extraction of Mining Areas from Landsat Imagery

In order to obtain a more accurate boundary of the mining area, Landsat 8 OLI imagery were used to map different land cover types using several unsupervised classification methods. The classified land cover types include mining area, water body, vegetation, bare land, and snowfield.

3.2.4. Estimation of Mining Landmass Volume

The ~30-m-resolution SRTM-1 DEM and AW3D were used to detect the topographic changes due to the detected mining activities. Before the estimation of mining landmass volume, the vertical accuracy was investigated for the two DEMs by referring to high-accuracy altimetry elevations measured by the ICESat satellite. Both SRTM-1 and AW3D achieved fairly high accuracy with mean absolute errors of 1.83 m and 2.64 m, respectively. Hence, the influence of different DEM sources on the volume change estimation is considered to be small. The fine-detailed mining extent mapped from the Landsat image was used as a spatial mask to calculate the landmass volume. This was conducted by computing the differences between the fills and cuts volume in ArcGIS software.

4. Results and Analyses

4.1. Abrupt Changes of NDVI Time Series Due to Mining Activities

According to the method in Section 3.2.1, there are a total of 152 pixels detected as mining-area breakpoints from all examined pixels. As shown in Figure 4, all the detected breakpoints are featured by centralized distribution which are consistent with the real situation of mining area as displayed in the high-resolution Google Earth imagery. Furthermore, based on the BFAST method, the years of dramatic change occurrence for each pixels are also identified. It seems that the mining area appear as a gradual expansion from 2003 to 2013 (Figure 4). The earliest mining operation was conducted in the center of our detected breakpoint cluster in 2003. Massive destructions of vegetation cover were mainly observed between 2008 and 2012.

In order to analyze the temporal changes in NDVI due to mining activities, the spatially averaged NDVI time series over BFAST-detected breakpoints from 2000 to 2018 is illustrated in Figure 5. There is a slight upward shift from 2000 to 2003, but an evident decline occurred after 2006 and the NDVI value became relatively stable with small seasonal fluctuations after 2011.

The applied MODIS imagery identified the initial time and general extent of the mining area. From these results, we next applied Landsat imagery to reveal how the boundary of this mining area evolved in finer spatial detail. We here compared three commonly-used unsupervised classification methods including maximum likelihood, minimum distance and Mahalanobis distance. As shown in Figure 6, the classification results based on the maximum likelihood method achieved the best performance, given our manual inspection against high-resolution Google Earth imagery. The main problem of the other two methods is the commission error surrounding the mining area, especially in the mountain regions. In this paper, the classification results based on the maximum likelihood method were adopted and converted into the final extents of the mining area. The post-processing is was conducted, including removal of small sliver polygons.

Although the NDVI value has proven to be sensitive to mining activities, the overall inter-annual variability within the mining area still needs to be investigated. The averaged NDVI values derived from Landsat of the mining area on 13 dates between 2003 and 2016 are illustrated in Figure 7. It should be noted that the NDVI value of 2013 is missing due to the lack of the cloud-free image in this year. In general, the NDVI values exhibit a decreasing trend which is particularly apparent during 2010 and 2015. However, many prior studies have indicated vegetation greening and advanced phenological green-up dates in the context of warming and wetting climate changes during the recent two decades [22,45,46]. The downward trend, especially with an abrupt drop around the year 2011, implies that mining activities have a significant negative impact on the local surface vegetation. However, NDVI values on June 23 2006 and August 31 2008 were slightly anomalous. This phenomenon was possibly explained by the image acquisition time when the vegetation was not so luxuriant in early summer and late summer.

4.2. Mining-Induced Landmass Volume Changes

In order to reveal the influences of open-mining activities on geomorphologic landform in the study area, we compared the topographic features of the study area between pre-mining and post-mining periods based on two different DEMs. As shown in Figure 8, the SRTM-1 DEM acquired in February 2000 represents the surface morphology before mining activities. By contrast, the terrain information provided by AW3D30 DEM around 2010 indicates that the mining activities had substantially changed the local topography. The post-mining elevation values in the central location are much lower than those of pre-miming values, implying a large amount of earth excavation induced by mineral/coal exploitation. As a result, the excavated sand/soil was deposited in the surrounding area, which led to evident elevation increase (Figure 8b) and might have suppressed the natural vegetation growth.

Figure 8c reveals the spatially-explicit changes in surface elevation after the detected mining activities. The elevation change for each pixel has been calculated by computing the difference between SRTM-1 DEM and AW3D30 DEM (Figure 8c). The deepest elevation drop (down to 175 m) is located in the central mining pit, and the major elevation aggradation (up to 128 m) is distributed right surrounding the mining pit.

Eventually, we applied the Cut Fill tool in ArcGIS to estimate the surface volume change due to the mining activities. The result shows that the excavated volume is 0.3231 km³ and the earth deposition volume is 0.3211 km³. These two values are highly comparable to each other, which further testifies the accuracy of our mining extent delineation and volume estimations. According to the statistical record (http://www.qhtjj.gov.cn/), the total coal reserve within and surrounding the identified mining area is up to ~3.54 billion tons, which is the largest open-pit coal mining area in Qinghai province.

5. Discussions

5.1. Sensitivity of BFAST Detection of Abrupt Changes

5.1.1. Sensitivity Test of the h Parameter in BFAST

In order to improve the robustness of detecting abrupt changes, we further test the sensitivity of key parameters for BFAST. In this section, the h parameter, which has a significant impact on the accuracy of BFAST detection, is investigated with different values. Pixel 5368 and pixel 5860 located in the mining area are chosen to detect significantly different break dates according to different h values, so they were taken as samples to acquire time series images processed by the BFAST algorithm from different h values of 0.05, 0.08, 0.10, 0.12, and 0.15.

When the h value is set as 0.05, that is, the minimal segment size from detected breaks in the trend component is about one year, there is no abrupt changes to be detected as shown in the first picture for Sample 5368 (Figure 9). That is because h value is too low to find a breakpoint within one-year cycle. When the h value is set as 0.08, there is one breakpoint and an abrupt change and the decomposed seasonality item also has an abrupt change. The trend components when h value is 0.10, 0.12 or 0.15 are similar and work well in detecting the abrupt change. The dating information of breakpoints (

Table 1. The dates that trend breakpoints and seasonal breakpoints occurred when h value is 0.05, 0.08, 0.1, 0.12, and 0.15, respectively. “None” means there is no breakpoints detected. 334 means the 333th eight-day of MODIS image slices which is based on the time resolution of MOD09A1, and 2007(12) means the 11th eight-day in 2007.

h	Trend Breakpoints			Seasonal Breakpoints
	Corresponding to Breakdates:			Corresponding to Breakdates:
	2.5%	Breakpoints	97.5%	2.5%	Breakpoints	97.5%
0.05		None			None
0.08	29 Marth 2008	24 November 2008	1 January 2009	16 November 2008	26 December 2008	02 February 2009
0.10	30 Marth 2007	26 December 2008	06 Marth 2009		None
0.12	30 Marth 2007	26 December 2008	06 Marth 2009		None
0.15	30 Marth 2007	26 December 2008	06 Marth 2009		None

) acquired from the result of BFAST further confirmed that the time when the abrupt change happened is consistent. Therefore, h value of 0.10, 0.12 or 0.15 is alternative. In general, assuming that the restoration measurements are not taken, the mining operation-induced vegetation destruction will last in much longer period than the periods of vegetation variations caused by climate variability. Thus the h value can be recommended for the segmented time window larger than 2–3 years or longer term.

As for Sample 5680, there is a breakpoint and an abrupt change in both trend and seasonal components when h value is equal to 0.05, 0.10 or 0.12. (Figure 10) By contrast, there is only one breakpoint in the trend component when h value is set to 0.08 or 0.15. From the breakpoints information in

Table 2. The date that trend breakpoints and seasonal breakpoints occurred when h value is 0.05, 0.08, 0.10, 0.12 and 0.15.

h	Trend Breakpoints			Seasonal Breakpoints
	Corresponding to Breakdates:			Corresponding to Breakdates:
	2.5%	Breakpoints	97.5%	2.5%	Breakpoints	97.5%
0.05	3 December 2006	27 December 2006	6 Marth 2007	2 February 2008	5 March 2008	6 April 2008
0.08	11 December 2006	17 January 2007	2 June 2007	None
0.10	3 December 2006	27 December 2006	06 Marth 2007	2 February 2008	5 March 2008	6 April 2008
0.12	3 December 2006	27 December 2006	06 Marth 2007	2 February 2008	5 March 2008	6 April 2008
0.15	11 December 2006	17 January 2007	2 June 2007	None

, the time of the abrupt changes of seasonal and trend component is the same when h value is 0.05, 0.10 and 0.12, and different from h value of 0.08 and 0.15. Thus, the h value of 0.05, 0.10 and 0.12 is alternative in this case. Consequently, h value of 0.10 can be an appropriate choice.

5.1.2. Sensitivity Test of the NDVI Different Threshold

It is necessary to optimize the parameters involved in the BFAST tool so that only the breakpoints due to mining activities can be detected. As shown in Figure 11, the detected breakpoints with different threshold values are overlaying with Landsat images. The mapping accuracies for different threshold setting were also validated with Landsat-based result (in

Table 3. The overall accuracy and kappa index of the confusion matrix based on Landsat image and each superimposed mask of breakpoints by different thresholds.

Threshold	0.05	0.08	0.10	0.12	0.15
Overall accuracy	0.9670	0.9736	0.9760	0.9783	0.9779
Kappa index	0.5320	0.5560	0.5552	0.5029	0.2741

).

It is obvious that the detected mining areas exaggerate the extent when the threshold value is set to 0.05 or 0.08. When the threshold is small, the second judging condition (Equation (3)) is easy to satisfy so that more pixels around the mining areas or in non-mining area with other disturbances may be included. By contrast, the extracted area is smaller than the reference extent when the threshold value is set to 0.12 or 0.15. By comparing with the result derived from Landsat mapping, the threshold value setting as 0.10 shows relatively more matched delineation. In the process of threshold analysis, it is easy to understand that although the different thresholds will affect the number of the NDVI time-series pixels identified with abrupt changes, the mining locations can be detected despite that the extents may be underestimated or overestimated. It indicates that the BFAST method can effectively detect the existence of mining-induced vegetation destruction based on time-series MODIS images. However, it is still necessary to combine the Landsat or other higher spatial-resolution imagery to assist the accurate extent extraction of mining areas over the BFAST-identified locations.

5.2. Influences of Mining Activities on Environmental Damages on the TP

The TP, also regarded as “Asian Water Tower”, is the source of the seven greatest rivers in Asia, which feed more than one-third of the world’s population [47]. The plateau is traditionally considered to have a pristine environment, and most of Tibetan rivers are considered to be free from contamination [48]. Recent studies suggest that pollutants originating from outside the region are transported through the atmosphere and deposited on the plateau [49,50]. Apart from this natural processes, several mineral deposit belts have been discovered or being explored, thus the mining developments are likely to become important sources of water contamination in the high-altitude areas. Within the study area, the mine’s water inlets and outlets both lie on the adjacent rivers that are mainly supplied with rainfalls and snow meltwater in the Qilian Mountains and that originally provided drinking water for local herders. The large volume of sediment deposited in rivers can cause aggradation in channels, and uncontrolled mine water discharges become a significant source of surface water pollution.

There are also serious concerns over the destruction of permafrost by open-cast mining over the TP where the permafrost is rather important for maintenance of regional water sources. The identified coalfield is located in the critical zone of permafrost regions in the TP [51,52]. It belongs to the alpine permafrost belt of the Altun–Qilian Mountains in the TP [51]. The permafrost thickness is 50–90 m, and the annual average temperature is −1.35 °C to −0.24 °C [52]. The mining activities will consequently change the heat-exchange conditions and water-heat transport process between the ground and the atmosphere [11,53]. Once the overlying vegetation is excavated, the underground ice will begin to melt. This will lead to a series of frozen-soil environment problems, including hot-melt collapse, slope mudflow, vegetation and soil degradation, and mining-area desertification. What is worse, the destruction of permafrost caused by open pit mining is often irreversible. According to the field survey carried out by China National Administration of Coal Geology and China Geological Survey, there is a great potential reserve of minerals, including coal and gas hydrate, located in broad permafrost regions of the TP, such as Wuli- Kaixinling region, Tuotuo River and Qiangtang Basin. In the future, a larger-scale investigation needs to be conducted based on multi-source remote sensing observations and more attentions should be paid on the mining development impacts in the environmental issues of the TP and other parts of northern China.

6. Conclusions

Identifying mining activities and assessing their effects on surface vegetation and landforms are of critical importance to local sustainable development and natural resources management. This is particularly true in vulnerable and pristine environments such as the Tibetan Plateau. In this study, we developed a fully automated approach that uses multiple remote sensing datasets to identify the timing and scale of mining activities. This approach was successfully demonstrated for the recent mining development in Muli Town situated in northeastern TP. Our results show that the BFAST algorithm is able to detect abrupt changes in vegetation cover from multi-temporal images, as a result of open-mining development for the studied alpine-grassland region. The design of a combined parameter setting in the BFAST algorithm can effectively differentiate mining-induced vegetation destruction from other types of vegetation changes caused by climate or grazing. A deeper investigation of the algorithm parameters (e.g., h and threshold of NDVI difference) suggests that the parameter setting may cause slight biases in the timing and number of detected pixels of abrupt changes, but does not impair the overall effectiveness of the proposed method in identifying mining locations. The testing study in Muli (northeastern TP) reveals that a mining development began in 2003 and massive destructions of vegetation cover followed between 2008 and 2012. The detection results from Landsat images and multi-temporal DEMs explicitly indicate that the mining activities have substantial destructive influences on local grassland and landform. The assessment of eco-environmental and geomorphologic damages caused by open mining development is expected to attract attentions from local environmental agencies and policy-makers in the TP. The demonstrated method may also be applicable to detecting and evaluating mining activities in other regions of the world.