The Trend of Coal Mining-Disturbed CDR AVHRR NDVI (1982–2022) in a Plain Agricultural Region—A Case Study on Yongcheng Coal Mine and Its Buffers in China

Abstract

The destruction of arable land caused by coal mining in coal grain compound areas is a major bottleneck restricting grain production increase. The spatiotemporal correlation between the decline in cultivated land quality and crop growth deterioration due to mining subsidence still needs to be clarified. This study employed the CDR AVHRR NDVI dataset and applied correlation and trend analysis methods to extract vegetation cover information from 1982 to 2022. It also explored the relationships between vegetation cover and temperature and precipitation. The study found the following: (1) Over the past 41 years, the NDVI in the study area showed a significant upward trend. Specifically, the average annual NDVI growth rate in the mining area was 51.85%, while the corresponding growth rates for the 10 km buffer area, 20 km buffer area, and check area (CK) were 65.91%, 65.86%, and 68.09%, respectively. The start of the growing season (SOS) for winter wheat in the mining area and control area advanced by 49 ± 1.5 days and 65 ± 1.5 days, respectively, while the length of the growing season (LOS) extended by 59 ± 1.5 days and 72 ± 1.5 days, respectively. For summer maize, the SOS advanced by 11 ± 1.5 days and 15 ± 1.5 days, respectively, and the LOS extended by 17 ± 1.5 days and 19 ± 1.5 days, respectively. The study area exhibited a significant positive correlation between the NDVI and temperature. Specifically, the correlation coefficient for the mining area was 0.6865 (p < 0.01); for the 10 km buffer zone, it was 0.5937 (p < 0.01), for the 20 km buffer zone, it was 0.6775 (p < 0.01), and for the control check area (CK), it was 0.6591 (p < 0.01). The results of this study can provide data support for the collaborative rehabilitation of and source reduction in coal grain compound areas, as well as for the restoration of damaged farmland.

1. Introduction

Surface subsidence caused by natural or anthropogenic factors has become a significant global development issue [1]. Coal, as a key fundamental energy source, occupies an important position in the global energy structure [2]. However, with the increase in coal mining activities, problems such as surface subsidence [3], land destruction [4,5], and decreased arable land productivity [6] have become increasingly severe. These issues not only have negative impacts on the environment but also directly threaten food security and sustainable agricultural development, becoming one of the critical bottlenecks limiting the growth of food production.

Against this backdrop, satellite remote sensing technology has been widely adopted due to its capability to efficiently monitor changes in mining areas and their surrounding environments [7,8]. Studies have shown that the continuous increase in coal production has caused irreversible damage to China’s ecological environment [9,10]. The primary manifestation of ecological destruction in mining areas is its impact on vegetation health. A time series of the normalized difference vegetation index (NDVI), one of the most widely used vegetation indices in remote sensing studies, can reflect the trend of vegetation in the ecological environment in coal-mining areas [11].

Over the past three decades, global warming has led to increased NDVI values, an advanced start of the growing season (SOS), a postponed end of the growing season (EOS), and an extended length of the growing season (LOS), which have been recorded in many parts of the world [12]. These changes have also been confirmed by related studies at domestic and international scales: in the context of global warming, the SOS in the Nordic Fennoscandia region became 11.8 ± 2.0 days earlier (p < 0.01) from 1982 to 2011 [13]. The LOS in the northern hemisphere was extended due to climate change (1982–2008), and the vegetation index increased significantly [14]. Approximately 56.30% of surface vegetation showed a significant seasonal increase from 1982 to 2011 [15].

Vegetation indices have been proved to be an effective method of capturing vegetation growth [16,17]. The NDVI is one of the most widely used vegetation evaluation indices in remote sensing images. The construction of a NDVI time series can quantify the ecological dynamics of different vegetation [18,19]. China’s national average annual NDVI increased by 7.4%, and the annual NDVI of Henan Province increased by 14.5%. These values were affected by driving factors such as air temperature and increases in summer precipitation [20,21]. The maximum NDVI of the vegetation cover in north China has increased significantly by 0.011/10a over the past 25 years [22]. Mining areas have been affected by climate-driven factors and human mining activities, and the findings of related studies of vegetation-based ecological change have been as follows: Wu et al. [23] used the SPOT 4/5 NDVI (1999–2008) to study the vegetation cover index of the Shendong mining area. The results showed that most of the vegetation cover in the mining area improved significantly (by 82.13%). Hu et al. [24] used Landsat TM images to describe the Shenfu mining area from 1986 to 2006. The vegetation cover index of the mining area improved overall. Quan et al. [25] analyzed the relationship between land use and ecological vulnerability by developing an ecological vulnerability model. The results showed that the exploitation of four open pit mines in the Shengli Coalfield resulted in a significant increase in the vulnerability of the mining area, and a 300–2000 m region on the periphery of the mine had evolved into a highly ecologically fragile area. Coal resource exploitation was the main factor leading to the increase in regional vulnerability. Zhang et al. [26] analyzed the spatial heterogeneity and driving mechanism of vegetation cover by multitemporal Landsat TM images. The results showed that the vegetation cover in the open area, the dump site, and the surrounding areas within 1.5 km had significantly decreased, and these decreases tended to be consistent. The vegetation restoration exhibited a “high-low-high” cycle over time. Ma et al. [27], via the GIMMS AVHRR NDVIg (1981–2006) and MODIS NDVI (2000–2010) datasets, found that the increase in NDVI values in the Shendong mining area was lower than that in the check area (CK). With the increase in mining intensity, the NDVI growth rate slowly decreased, and it was lower than that in the CK. However, the temporal scale of inference was short or intermittent, owing to a lack of multiscale spatial comparison and climate-related analysis. Using the GIMMS AVHRR NDVI3g (1982–2013) dataset, Ma et al. [28] found that in the Lu’an mining area, the annual NDVI growth rate during the study period ranked among different areas as follows: mining area (1.09%) < 10 km buffer area (2.16%) < 20 km buffer area (8.86%) < CK (9.87%). Thus, mining activity indisputably had a significant impact on regional ecological processes. Karan, S.K. et al. [29] assessed reclamation in the Jharia coal field using Landsat imagery from 2000 to 2015. Vegetation cover on reclaimed sites grew by 21%, from 213.88 to 258.9 hectares, showing effective reclamation. Blachowski, J. et al. [30] conducted an integrated spatiotemporal analysis of vegetation in the rehabilitated areas of the Przyjaźń Narodów-Szyb Babina brown coal mine, western Poland, using Sentinel-2 data. The study found significant vegetation changes related to pre-mining activities even 50 years post-mining. Buczyńska, A. et al. [31] analyzed vegetation development and post-coal mining activities in the Babina mine area, western Poland, from 1989 to 2019 using Landsat TM data. The results indicated overall vegetation improvement, particularly near coal waste dumps and former open-pit mines. Juanda, E.T. et al. [32] analyzed vegetation changes in a coal mining restoration area in South Kalimantan, Indonesia, using Sentinel-2 imagery and UAV photogrammetry. The NDVI significantly increased to 0.369 and 0.417, with healthy vegetation percentages of 68.13% and 81.39%, respectively, in the third and fourth years post-restoration.

In this paper, we compare the effects of mining activities on vegetation health among different climatic types from the perspectives of temporal, spatial, and climate change. The Yongcheng mining area, a farming plain with a healthy ecological environment and vegetation with a high capacity for self-repair, was considered a study case. The intra-annual and inter-annual variability, vegetation cover successional patterns, and LOS in the mining area, indirectly affected areas, and an unexploited area (i.e., the CK) will be discussed. Moreover, NDVI changes, which may be attributed to climate change or to human activities, will be identified. The results of this study will provide experimental support and a theoretical basis for the comprehensive management of coal mine environments.

2. Materials and Methods

2.1. Study Area

(1)

Mining area

The Yongcheng mining area is located in Yongcheng City, which is in the eastern part of Henan Province, and adjacent to western Anhui Province. The geographical position of the mining area is 33°42′–34°18′ E, 115°58′–116°40′ N, with a total area of approximately 85,000 hm² [33]. The area consists of five coal mines, namely, the Chensilou mine, Chengjiao mine, Cheji mine, Shunhe mine, and Xinqiao mine, as shown in Figure 1.

The Yongcheng mining area is located in the central part of the Huang-Huai Plain in China. The geomorphic landscape of the area is characterized by flatlands, with a semi-humid and semi-arid continental monsoon climate. The annual average temperature is 14.3 °C, and the annual average precipitation is 874 mm [34]. The region has a simple soil type and comprises mainly agricultural land. The region is an important grain production base for China, producing mainly agricultural plants such as wheat and maize. Because of the good ecological environment and high degree of cultivation in this region, its environmental adaptability and resistance are strong.

(2)

Buffers and the CK

To further reveal the response to coal-mining activities in the mining and surrounding areas, this study drew on previous studies of grassland ecology [35]. Based on analysis of the mining area, two buffers (10 km and 20 km), and CK (Figure 1), the vegetation ecology of mining-disturbed and non-mining-disturbed areas was studied. Buffer areas are types of geographical spatial targets that represent areas of influence or service ranges around point, line, or polygon entities. Specifically, they refer to the automatic establishment of areas with certain widths around these entities. Using the “Proximity” analysis tool in GIS, a “buffer” function was employed to create buffers of 10 km and 20 km around the mining area boundary line. To reflect the process of vegetation growth in natural ecological environments, the CK was selected based on the principle of pseudo-invariant feature areas (PIFs) (Jensen’s Book) [36]:

• The CK should be a vegetation-covered area.
• The CK and the experimental area should be at the same latitude. The coordinates of the central point of the mining area are (115°5′48″ E, 34°11′8″ N), and the coordinates of the central point of the verification area are (116°23′19″ E, 33°56′53″ N). Here, their longitudinal difference is 27′31″ and latitudinal difference is 14′15″.
• The geometric center of the CK is 50 km away from that of the mining area, which is less affected by mining activities.
• The slope and elevation in the CK are similar to those in the mining area. The average slope of the mining area is 1.71°, while the average slope of the verification area is 2.07°.

2.2. Data Use

The long-term CDR AVHRR NDVI dataset (1982–2022) was acquired from NOAA’s National Centers for Environmental Information (http://www.NCEI.noaa.gov, accessed on 5 March 2023). The temporal resolution is 1 d, with a spatial resolution of 0.05 degrees per pixel. Bimonthly data were calculated by the maximum value composite (MVC) procedure. The dataset covered 984 half-months of synthetic NDVI images.

The PKU GIMMS NDVI long-term dataset (1982–2020) is from the zenodo website (https://zenodo.org, accessed on 10 April 2023). The dataset is available in geographic latitude/longitude projection and TIFF formats. This dataset has a temporal resolution of half a month and a spatial resolution of 1/12. The dataset contains the global vegetation area from 1982 to 2020 [37].

The meteorological data consisted of a 1 km monthly mean temperature dataset for China (1901–2022) and 1 km monthly precipitation dataset for China (1901–2022), which came from the National Tibetan Plateau Data Center (http://data.tpdc.ac.cn, accessed on 15 May 2023). In order to obtain the annual mean meteorological data within the study area, annual mean synthesis was performed on the data within the study area.

2.3. Data Processing Method

(1)

Average and cumulative NDVI

NDVI is the best indicator of vegetation growth states [38,39]. The annual average NDVI can be used to analyze changes in the vegetation index [40,41]. The annual accumulated NDVI can reveal the trend of regional vegetation biomass during the study period and provide a basis for vegetation management. In this study, half-monthly composite values were calculated into annual means and cumulative values using an algorithm based on the IDL programming language, and the NDVI trend curve was then constructed in Origin to analyze monthly and annual variations in the vegetation index in the area. The calculation equations were as follows:

$$a={\displaystyle \sum _{i=1}^m{b}_i}\left\{\begin{array}{c}Sum={\displaystyle \sum _{j=1}^n{a}_j}\\ Average=Sum/m\cdot n\end{array}\right.$$

(1)

Equation (1) expresses a progressive calculation process. Here, $b_i$ is a single-pixel NDVI value composited bimonthly; m is the cumulative number of composites and is equal to 2 (when implementing a monthly MVC), 6 (when implementing a quarterly MVC), or 24 (when implementing an annual MVC); Sum is a cumulative value; and n is the number of pixels over the area.

(2)

Gaussian fitting

The Gaussian function is a probability density function that represents a continuous random variable. It is often used to express the periodic evolution of phenology in nature [42]. In this study, based on the characteristics of the dual growing seasons of winter wheat and summer maize in the Yongcheng mining area, a Gaussian fitting method (double-peak type) was adopted to smooth the NDVI time series. This approach allowed for the better identification of trend changes within the NDVI time series. The one-dimensional expression of the function (bimodal type) is as follows:

$$y=y_0+{\displaystyle \frac{A_1}{w_1\cdot \sqrt{\pi /2}}}{e}^{-{\displaystyle \frac{2{(x-x_{{\mathrm{c}}_1})}^2}{w_1^2}}}+{\displaystyle \frac{A_2}{w_2\cdot \sqrt{\pi /2}}}{e}^{-{\displaystyle \frac{2{(x-x_{{\mathrm{c}}_2})}^2}{w_2^2}}}$$

(2)

In Equation (2), $y_0$ is an offset of the baseline. Ai is the integral area on the underside of the bell curve (here, A1 and A2 represent the NDVI-equivalent biomass of winter wheat and summer maize, respectively). $x_{\mathrm{ci}}$ is the height in the central peak (equivalent mathematical expected value). $w_{\mathrm{i}}$ = 2σ is approximately equal to 0.849 times the width at the peak half-height (the width of the bell curve when the peak height is 1/2). The value implies the length of growth (LOS) during prosperous vegetation.

According to the characteristics of the Gauss function, the coordinates of the four inflection points of the fitted curve can be obtained:

$$(x_{c1}\pm {\displaystyle \frac{w_1}{2}},y_0+{\displaystyle \frac{A_1}{w_1\cdot \sqrt{\pi e/2}}}+{\displaystyle \frac{A_2}{w_2\cdot \sqrt{\pi e/2}}})(x_{c2}\pm {\displaystyle \frac{w_2}{2}},y_0+{\displaystyle \frac{A_1}{w_1\cdot \sqrt{\pi e/2}}}+{\displaystyle \frac{A_2}{w_2\cdot \sqrt{\pi e/2}}})$$

(3)

That is, the time when the start of the season (SOS) and the end of the season (EOS) occur (the day of the year, DOY). The length of the growing season (days) (LOS) can then be calculated.

(3)

Linear regression

The greenness rate of change (GRC) for vegetation was quantified through linear regression analysis, represented by the slope of the normalized difference vegetation index (NDVI) regression line over a specified period [43]. This study aimed to elucidate the temporal trends in the NDVI across different regions during this interval. Specifically, the inter-annual NDVI time series for the mining areas, buffer areas, and CK were fitted using linear regression models, with their respective slopes being indicative of ecological quality dynamics within these areas. Furthermore, the study conducted linear regression fitting of the inter-annual variations in the SOS, EOS, and LOS to reflect phenological changes in the vegetation across different regions. The equation for the slope is presented as follows:

$${\Theta }_{slope}={\displaystyle \frac{{\displaystyle \sum _{i=1}^n(x_i-\overline{x})(y_i-\overline{y})}}{{\displaystyle \sum _{i=1}^n{(xi-\overline{x})}^2}}}$$

(4)

Here, $\overline{x}={\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^n{x}_i}$ and $\overline{y}={\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^n{y}_i}$. $x_i$ is the independent variable, $y_i$ is the dependent variable, and ${\Theta }_{slope}$ represents the slope of the linear regression equation.

(4)

NDVI change rate

To elucidate the inter-annual NDVI differences between the mining-disturbed and non-mining-disturbed areas, the annual NDVI change rate was unified and calculated using the following equation:

$$r_{NDVI}(\%)={\displaystyle \frac{{\Theta }_{slope}}{Average}}\times k\times 100\%$$

(5)

Here, ${\Theta }_{slope}$ is obtained by the regression line between the annual mean NDVI over 41 years (1982~2022) (determined by Equation (4)). k is the statistical period; Average represents the arithmetic mean of the cumulative value of NDVI over k years (determined by Equation (1)). The parameter r is used to compare the trend of vegetation activity in the mining-disturbed zone with that in the non-mining-disturbed zone.

3. Results

3.1. The Intra-Annual Variation in the NDVI

The mean (from January to December) values of the half-monthly NDVI for the whole 41-year period (1982~2022) were used to measure the intra-annual variation characteristics. These values were charted and implemented via an intra-annual timing feature analysis (Figure 2a–d, light gray polyline). The variation in the NDVI in the four areas during the 41 years was characterized by a double-peaked Gauss curve, which confirmed the phenological characteristics of the winter wheat and summer maize in northern China. The NDVI was lowest at the end of October, and the first peak was approximately 0.61 (the winter wheat growing season) at the end of April. The first valley was approximately 0.32 at the end of June (the NDVI value decreased before and after the summer harvest). The second peak of approximately 0.65 appeared at the end of August (the summer maize growing season). In late February and June, the NDVI value increased significantly. There were two rapid declines at the end of May and September, and the second valley was approximately 0.28 at the end of October. In the four areas, as shown in Figure 2a–d, the NDVI changed slightly in January and February. These slight changes indicated that the termination of vegetation growth was affected by low temperatures. In late February, the NDVI increased steeply, and the natural vegetation in the experimental area entered the first SOS. The change in the NDVI became steep in mid-June, and the natural vegetation in the experimental area began to enter the first EOS. At the end of June, the NDVI increased sharply, and the natural vegetation in the experimental area entered a second SOS. In mid-October, the NDVI decreased rapidly, and the natural vegetation in the experimental area entered the second EOS. These phenological variabilities were consistent with the previously reported trend of vegetation cover in Henan Province [44]. Based on the spatial analysis, there was no obvious spatial difference in the trend of the NDVI among the four areas, which indicated that the vegetation types, ecological environment, and farming modes were consistent in the four areas. This similarity indicated that the sites of the calibration had a high consistency. Bimodal models of the four areas were obtained by multi-peak Gaussian fitting (Equation (2)), and the goodness of fit (R²) was 0.9231 (the mining area), 0.9216 (the 10 km buffer), 0.9360 (the 20 km buffer), and 0.9275 (the CK).

The Gaussian fitting functions are expressed as follows:

In the mining area,

$$y_a=0.3311+{\displaystyle \frac{1.3081}{3.4903\cdot \sqrt{\pi /2}}}{e}^{-{\displaystyle \frac{2{(x-7.8650)}^2}{{3.4903}^2}}}+{\displaystyle \frac{1.3734}{3.0622\cdot \sqrt{\pi /2}}}{e}^{-{\displaystyle \frac{2{(x-15.6663)}^2}{{3.0622}^2}}},R^2=0.9232$$

(6)

In the 10 km buffer,

$$y_b=0.3202+{\displaystyle \frac{1.4921}{3.6234\cdot \sqrt{\pi /2}}}{e}^{-{\displaystyle \frac{2{(x-7.8261)}^2}{{3.6234}^2}}}+{\displaystyle \frac{1.5127}{3.0608\cdot \sqrt{\pi /2}}}{e}^{-{\displaystyle \frac{2{(x-15.6526)}^2}{{3.0608}^2}}},R^2=0.9216$$

(7)

In the 20 km buffer,

$$y_c=0.3124+{\displaystyle \frac{1.4206}{3.5901\cdot \sqrt{\pi /2}}}{e}^{-{\displaystyle \frac{2{(x-7.8616)}^2}{{3.5901}^2}}}+{\displaystyle \frac{1.5397}{3.1729\cdot \sqrt{\pi /2}}}{e}^{-{\displaystyle \frac{2{(x-15.6589)}^2}{{3.1729}^2}}},R^2=0.9361$$

(8)

In the CK,

$$y_d=0.3265+{\displaystyle \frac{1.3875}{3.5810\cdot \sqrt{\pi /2}}}{e}^{-{\displaystyle \frac{2{(x-7.9974)}^2}{{3.5810}^2}}}+{\displaystyle \frac{1.4864}{3.3106\cdot \sqrt{\pi /2}}}{e}^{-{\displaystyle \frac{2{(x-15.6436)}^2}{{3.3106}^2}}},R^2=0.9275$$

(9)

3.2. Analysis of Vegetation Growth Period

Phenological changes are the result of complex interactions between temperature, precipitation, solar radiation, and many other factors. An effective method for phenology identification or modeling is to fit long time series information and then analyze its periodic variation [45]. In the agricultural area, to analyze the Gauss multi-peak fitting curve of the average NDVI of the winter wheat and summer maize can well reveal the phenology of the crops (Figure 3a–d). Here, A1 and A2 represent the NDVI-equivalent biomass of the winter wheat and summer maize, respectively. Width1 and Width2 are approximately equal to 0.845 times the width at the peak half-height (the width of the bell curve when the peak is at half-height). The value implies the length of growth during prosperous vegetation (unit: month). Heigth1 and Heigth2 are the heights in the central peak and represent the maximum NDVI (NDVImax) of the winter wheat and summer maize, respectively. The x_c1 and x_c2 are the time of occurrence of the NDVImax of the winter wheat and summer maize (unit: half-month).

There are many ways to monitor the season of vegetation growth by means of the NDVI, such as the NDVI threshold method [46,47,48], the moving average method [49,50], and the NDVI midpoint method [51]. In the study, the vegetation growing season was monitored and analyzed by the NDVI threshold method.

The half-monthly NDVI values in the experimental areas were interpolated in MATLAB (Figure 4). The vegetation in Yongcheng is mainly composed of winter wheat and summer maize. In order to better distinguish the growth status of the winter wheat and summer maize, a threshold value was selected according to the change in the NDVI and divided into two LOSs. In the winter wheat LOS, the vegetation started to change rapidly when the NDVI was greater than 0.56, and when the NDVI was less than 0.44, the vegetation was out of the rapid change state. In the summer maize LOS, when the NDVI was above 0.55, the vegetation began to rapidly change, and when the NDVI was less than 0.40, the vegetation left the state of rapid change. According to the NDVI threshold method, in the winter wheat LOS, when the NDVI increased (>0.56) or decreased (<0.44), it corresponded to the beginning and end of vegetation photosynthesis. In the second LOS, when the NDVI increased (>0.55) or decreased (<0.40), it marked the beginning and end of vegetation photosynthesis. This means that in the winter wheat LOS, the vegetation began the SOS when the NDVI value was greater than 0.56 and entered the EOS when it was less than 0.44. In the summer maize LOS, the vegetation began the SOS when the NDVI value was greater than 0.55 and entered the EOS when it was less than 0.40.

The trend of the SOS and EOS: With a linear regression analysis of the SOS and EOS (1982–2022), we determined the trend of inter-annual phonological change (Figure 5). A negative “regression scissors” pattern was observed between the SOS and EOS from 1982 to 2022. Long-term changes in the beginning and end of growing seasons have been observed based on trends observed in phenological indicators, resulting in an increase in the length of growing seasons [52]. Therefore, the LOS of the vegetation showed an increasing trend over the 41 years. The average LOS of the winter wheat was 79.16 d (the mining area), 86.27 d (the 10 km buffer), 81.33 d (the 20 km buffer), and 80.96 d (the CK). The average LOS of the summer maize was 67.88 d (the mining area), 66.90 d (the 10 km buffer), 81.33 d (the 20 km buffer), and 69.33 d (the CK).

Historically, the types of local farming have varied dramatically; therefore, it was more reasonable to select the 1998–2022 data to study the vegetation phenology.

During the period from 1998 to 2022, the SOS of the winter wheat was advanced by 49 ± 1.5 (mining area), 53 ± 1.5 (10 km buffer area), 62 ± 1.5 (20 km buffer area), and 65 ± 1.5 (CK) days, while the EOS was delayed by 11 ± 1.5 (mining area), 10 ± 1.5 (10 km buffer area), 12 ± 1.5 (20 km buffer area), and 15 ± 1.5 (CK) days, respectively.

For the summer maize, the SOS was advanced by 11 ± 1.5 days (mining area), 3 ± 1.5 days (10 km buffer), 9 ± 1.5 days (20 km buffer), and 15 ± 1.5 days (CK), while the EOS was postponed by 6 ± 1.5 days (mining area), 7 ± 1.5 days (10 km buffer), 8 ± 1.5 days (20 km buffer), and 5 ± 1.5 days (CK), respectively.

The phenological variability differed among the four regions. For the winter wheat, the LOS of the mining area experienced a 59 ± 1.5 d extension, and the LOSs of the buffers (10 km and 20 km) were extended by 62 ± 1.5 d and 72 ± 1.5 d. In addition, the LOS of the CK underwent a 71 ± 1.5 d expansion. The linear growth rate in the mining area was 2.35 d/a, and the growth rate in the CK was 2.85 d/a.

For the summer maize, the LOS of the mining area was extended by 17 ± 1.5 d, while the buffers (10 km and 20 km) experienced extensions of 10 ± 1.5 d and 17 ± 1.5 d. Furthermore, the CK was extended by 19 ± 1.5 d. The linear growth rate in the mining area was 0.67 d/a, and the growth rate in the CK was 0.75 d/a. Using the dot matrix interpolation method, the time resolution was improved from the original 15 days to 1.5 days.

3.3. The Inter-Annual Variability in the NDVI

The influence of mining disturbance on the NDVI of the mining area was analyzed by comparing the NDVI between the mining-disturbed and non-mining-disturbed areas. Figure 6 shows that the four experimental areas were highly similar. The Pearson correlation coefficients between the mining area and the other three areas were 0.9906 (10 km buffer), 0.9905 (20 km buffer), and 0.9838 (CK). The average annual NDVI value of the CK was lower than that of the mining area in the 1980s. However, in recent years, the NDVI in the CK has begun to overtake that in the mining area. Under the same ecological environment, the vegetation activity in the CK was more active than that in the mining area.

The θ_slope in the mining area obtained by the regression equation (Equation (3)) was 0.0056, and the regression slopes in the other three areas were 0.0072 (10 km buffer area), 0.0070 (20 km buffer area), 0.0074 (check area). The mean NDVI of the four areas over the past 41 years was 0.4428 (the mining area), 0.4454 (the 10 km buffer), 0.4339 (the 20 km buffer), and 0.4462 (the CK). According to Equation (4), the annual mean (1982–2022) NDVI growth rates r_NDVI were 51.85% (in the mining area), 65.91% (in the 10 km buffer), 65.86% (in the 20 km buffer), and 68.09% (in the CK). Spatially, the differences between the validation area and the mining area are significant. In the mining area, 38.09% of the region has a regression slope below 0.0055, while in the validation area, it consistently remains above 0.0055. The average regression slopes across the four major regions are 0.0054 (the mining area); 0.0070 (the 10 km buffer); 0.0067 (the 10 km buffer); 0.0070 (the CK) (Figure 7a–d). The growth rate of the mining area was the lowest of the four experimental areas, which may have been related to the long-term impact of coal mining. Assuming that the vegetation growth rate in the CK is natural, the deficits in the growth rate of the other regions can be partly attributed to the impact of mining activities. The contribution rates of mining disturbance to the regional NDVI changes were calculated to be −16.24% in the mining area, −2.18% in the 10 km buffer, and −2.22% in the 20 km buffer. Thus, mining disturbance had a certain impact on the vegetation growth in the mining area. The Yongcheng mining area has superior climatic conditions, superior soil that is suitable for fine-scale farming, and vegetation with a strong self-repairing ability. Therefore, vegetation growth in the area and its surrounding areas is less affected by coal mining than that in other areas [24,25,26].

3.4. Climate Correlation Analysis

Characteristics of temperature changes in the mining area: The lowest average temperature in the experimental area was 14.43 °C (in 1984), the highest average temperature was 16.31 °C (in 2007), and the average temperature over the 41 years was 15.47 °C. Over the past 41 years, the change in the annual average temperature in the Yongcheng mining area was consistent with the trend of the change in global temperature. The overall trend exhibited fluctuations (Figure 8a,b). The rate of temperature increase was 0.27 °C/10a, which was lower than the national rate of increase of 0.29 °C/10a over the past 50 years (1951–2004) [53] and even higher than the rate of change in the annual average temperature of 0.17 °C/10a (1961–2013) in Henan Province [54]. Moreover, the rate was far higher than the rate of 0.12 °C/10a of global temperature change for the same period [55,56].

Characteristics of precipitation changes in the mining area: As shown in Figure 8c,d, contrary to the trend of precipitation change across the country, the precipitation in the Yongcheng mining area over the 41 years exhibited a slight increasing trend, but the annual precipitation in the mining area has decreased significantly in recent years [57]. The average annual precipitation in the experimental area in 1999 was the lowest, at 534.56 mm. In 2003, the precipitation increased steeply to 1113.66 mm, which was the highest precipitation during the study period. The average precipitation in the mining area was 759.06 mm during the 41 years, and the rate of increase in the precipitation was low, i.e., only 1.78 mm/a.

Correlation between the NDVI and climatic factors: Temperature, precipitation, moisture, light, and other climatic factors are important factors affecting the phenology of vegetation [58,59]. Based on the correlation analysis between the average NDVI and the temperature or precipitation (Figure 9a–d), the NDVI of both the mining area and the buffers was positively correlated with the precipitation: 0.20043 (the mining area, p > 0.05), 0.18085 (the 10 km buffer, p > 0.05), 0.23015 (the 20 km buffer, p > 0.05), and 0.17704 (the CK, p > 0.05). In spatial terms, the correlation between the NDVI and precipitation in both the mining and comparison areas consistently remained below 0.25, with average correlations of 0.22924 (the mining area, p > 0.05) and 0.18530 (the CK, p > 0.05) (Figure 9c,d). However, these correlations were not statistically significant. There was a significant positive correlation between the NDVI and the temperature-related driving factors. The correlation coefficients were 0.68652 (the mining area, p < 0.01), 0.59372 (the 10 km buffer, p < 0.01), 0.67759 (the 20 km buffer, p < 0.01), and 0.65912 (the CK, p < 0.01). Spatially, the NDVI in both the mining area and CK showed a significant positive correlation with the temperature, ranging from a minimum of 0.53281 to a maximum of 0.71994 (Figure 9a,b). Therefore, the sensitivity of the vegetation to temperature was much higher than that to precipitation; this result is consistent with the findings of previous studies [60,61,62,63,64]. Against the background of global warming, vegetation activities increased. The mining areas were the most sensitive to temperature, which was also a response to global climate change.

Through correlation analysis, this study determined that the vegetation index was significantly correlated with the temperature, and the vegetation index of the mining areas responded strongly to global warming. This finding was in contrast to the finding that the process of desertification in the mining area has been more affected by climate change [65]. In addition to precipitation and temperature, other climatic factors, such as humidity, wind, potential evaporation, sunshine days, and other factors, affect vegetation. The comprehensive analysis of vegetation cover variations in different ecoregions of a typical system demonstrates that vegetation cover variations are affected by the combined or synergetic action of hydrothermal conditions such as temperature, precipitation, and humidity. Therefore, future studies will provide a more comprehensive discussion about the impact of the above-mentioned climatic factors on vegetation cover in mining areas.

4. Discussion

The mining area ecosystem, which is impacted by many factors, including artificial mining activities and climate change, is a fragile mixed ecosystem. Mining is a strong disturbance that affects the structure and functioning of ecosystems. Therefore, the strength of a mine’s ability to recover from a disturbance is related to the natural resilience of the mine’s ecosystems [66]. Long-term and high-intensity mining are bound to have an impact on vegetation activities. Many scholars believe that the vegetation damage caused by exploitation under certain conditions can be improved by artificial vegetation reclamation [67,68]; however, an increase in the vegetation index after artificial restoration within a short period of time does not necessarily imply the normal evolution of the original vegetation in the mining area. Based on the addition of the buffers and the CK, the influence of coal mining on vegetation was analyzed, and the results were determined. However, discovering how to quantify the differences between the effects of artificial disturbance and climatic factors on the evolution of the vegetation in the mining area represents the next step that must be overcome.

(1)

Comparison of Vegetation Affected by Mining Disturbance with Naturally Grown Vegetation

Coal mining activities can lead to large-scale, high-intensity ecological disturbances that can destroy large amounts of land resources [69]. In addition, mining activities have a significant impact on vegetation during the growing season [70]. Many studies have shown that mining has a greater impact on vegetation in the mining area than on that outside the area [24,25,26]. However, the true impact of coal mining on the vegetation in the mining area can be determined only if the natural growth rate is separated from the total growth rate [27,28,71]. This study assumed that the vegetation growth rate of the CK was the natural growth rate. Although the impact of coal mining on the vegetation ecology may be very subtle (mainly because the Yongcheng mining area has fertile soil, suitable climatic conditions, and good field management, such that its vegetation (mainly crops) has a strong self-repairing ability), long-term analysis still detected some changes. We furthered the contribution rate of human activities and natural changes from the NDVI trend by cartographic analysis. Under the same climatic conditions, we believe that the differences (approximately 8.0–9.3% of the total NDVI) between the three affected areas and the CK are mainly caused by influencing mining factors (Figure 10a–c).

Because of the multiple noise components in the initially generated NDVI time series, smoothing was required to obtain a smoother and more stable time series [72]. In Figure 10a–c, the yellow line is the original NDVI. The blue line was generated by smoothing to remove irregularities and noise from the NDVI dataset using a Svaitzky–Golay filter. The red or blue bars are the residual NDVI between the three affected areas and the CK. The smoothed and residual NDVI were generated by linear fitting, respectively. The slopes of the NDVI and residual NDVI were obtained (

Table 1. Slope and precision of linear regression in time series NDVI. Note: The upward arrow indicates an increase, while the downward arrow indicates a decrease.

Area	NDVI Type	Slope	Standard Error	Fitting Adj. R-Square
Mining area	Original	0.0002325	0.0000164	0.16953
	Smoothed	0.0002326	0.0000156 ↑	0.18458 ↑
	Residual	−0.0000754 ↓	0.0000058	0.14666
10 km buffer	Original	0.0002971	0.0000179	0.21853
	Smoothed	0.0002972	0.0000171 ↑	0.23583 ↑
	Residual	−0.0000108 ↓	0.0000029	0.00348
20 km buffer	Original	0.0002895	0.0000174	0.22028
	Smoothed	0.0002896	0.0000166 ↑	0.23713 ↑
	Residual	−0.0000184 ↓	0.0000055	0.01125

).

As shown in Table 1, the linear regression slope of the NDVI, standard error of the slope, and R-square were calculated. Curve smoothing processing was performed on the original NDVI data before the linear regression analysis. The method of the curve smoothing was a Svaitzky–Golay filter, and the step size in this method was 3. It can be seen that the accuracy (standard error) and the goodness of fit (R-square) for the NDVI slope calculation were improved after the smoothing process. Most importantly, the variation in the residual NDVI (the NDVI difference between the affected area and the CK) was obtained. That is, there was a negative impact of human activities on regional vegetation growth. The slope of the residual NDVI of three research areas (mining area, 10 km buffer, and 20 km buffer) are −0.0000754, −0.0000108, and −0.0000184, respectively. The sums of the residual NDVI values are −3.3316 (the mining area), −0.7825 (the 10 km buffer), and −10.3479 (the 20 km buffer), which shows that there is no obvious spatial pattern in the sums of the residual NDVI. The sum of absolute values of residual NDVI is 40.5880 (the mining area), 36.2022 (the 10 km buffer), and 34.3726 (the 20 km buffer). It is obvious that the rank of the sum of the absolute values of the residual NDVI values is the mining area, the 10 km buffer, and then the 20 km buffer. The farther away they are from the mining area, the smaller the fluctuation of the residual values are.

(2)

Comparison of Different NDVI Datasets

The NDVI is an important tool for monitoring changes in surface vegetation, and NDVI datasets from different sources and at different scales are used differently and can be cross-checked with each other [73]. The PKU GIMMS NDVI dataset was selected for comparison with the CDR AVHRR NDVI dataset, and both of them used half-monthly maximum synthetic data to obtain the annual mean values. A comparison of the trends over the study period showed that both sets of data showed a significant increase in the four study areas: the mining area, 10 km buffer area, 20 km buffer area, and check area (Figure 11 and Figure 12).

In terms of the correlation between the two datasets, in the four regions, the correlation between the annual averages of the PKU GIMMS NDVI and CDR AVHRR NDVI is as follows: 0.66912 in the mining area (p < 0.01), 0.841 in the 10 km buffer (p < 0.01), 0.8346 in the 20 km buffer (p < 0.01), and 0.76864 in the check area (p < 0.01), with the highest correlation in the 10 km buffer area and the lowest correlation in the mining area. Figure 12a–d present scatter plots generated from the monthly PKU GIMMS NDVI and CDR AVHRR NDVI data. Linear fitting of the scatter plots reveals a significant correlation between the two datasets, with correlation coefficients of 0.82 (mining area, p < 0.01), 0.84 (10 km buffer area, p < 0.01), 0.85 (20 km buffer area, p < 0.01), and 0.78 (check area, p < 0.01). Analysis of the scatter distribution indicates that the NDVI values from the two datasets fluctuate over time, with the PKU GIMMS NDVI generally showing higher values than the CDR AVHRR NDVI. The slope of the fitted line is 1.38, reinforcing this observation. The slopes of the fitted lines for annual and bi-monthly scales are provided in

Table 2. PKU GIMMS NDVI and CDR AVHRR NDVI fitted straight line slopes.

Scale	Dataset	Mining Area	10 km Buffer	20 km Buffer	Check Area
Yearly	CDR AVHRR NDVI	0.00649	0.00742	0.007	0.00737
	PKU GIMMS NDVI	0.00155	0.0025	0.00253	0.00219
Monthly	CDR AVHRR NDVI	0.81	0.84	0.90	0.87
	PKU GIMMS NDVI

.

(3)

Comparison of Different NDVI Datasets

By comparing the four mining areas, it can be found that the ecological environment and vegetation types are different, and their vegetation index and growth period are also different. Figure 13 provides a brief comparative introduction of the geographical location and basic information of the four mining areas.

There are differences in the ecological characteristics of the vegetation. The annual vegetation growth pattern of the Shendong mining area is obviously “unimodal”, and the fitting curve is similar to the “Voigt” curve with wider bottom. Because it is located in an arid and semi-arid area, the environment of the mine area is bad and the vegetation growth is relatively depressed, so the vegetation biomass is the lowest. The annual vegetation growth pattern of the Binchang mining area is like a “trapezoid”; about 30% of the land in the mining area is used for agriculture, which causes the vegetation index to drop slightly in June after the peak in May, but it returns to a peak in August. The fitted curve presents a “Lorentz” line with a wider base and no peak. The annual vegetation growth in the Lu’an mining area shows a “single peak type with small shoulders”, and the variation in the “main and secondary double peaks” appears in the area. There are two peaks in May (secondary peak) and August (main peak), while the vegetation index begins to decline in late May (early spring crop harvest) and September (corn harvest). This is in line with the growth cycle of small-scale winter wheat and mainly summer maize and also confirms the local farming situation of two crops per year. The fitted curve presents a narrow “Doppler” type of curve at the bottom. The vegetation change in the Yongcheng mining area shows the characteristic of “double peaks”. The mining area is mainly composed of agricultural land and is an important grain production base in China, mainly growing crops such as wheat and maize. The fitting curve of the vegetation index shows a “Gaussian bimodal” line with a wider bottom.

There are differences in the phenological characteristics of the vegetation. During the period of 1981–2006 in the Shendong mining area, due to the influence of human activities, the SOS was earlier than in the natural ecoregion, the EOS lagged behind that of the natural ecoregion, and the LOS was extended again [27]. The vegetation LOS in the Binchang mining area showed a trend of extension during 32 years from 1982 to 2013, and the LOS of the mining area was extended by 11 ± 3 d, with a linear growth rate of 0.34 d/a [71]. From 1982 to 2013, the calculated vegetation SOS in the Lu’an mining area showed a lagging trend, the EOS showed an advancing trend, and the LOS showed a shortening trend. The SOS lagged by 3 d, the EOS advanced by 30 d, and the LOS was shortened by 33 d [28]. The Yongcheng mining area was divided into two growth periods within a year based on the growth conditions of local winter wheat and summer maize. During the period from 1998 to 2022, the SOS for the winter wheat advanced by 49 ± 1.5 d, the EOS was delayed by 11 ± 1.5 d, and the growth period was extended by 59 ± 1.5 d. The linear growth rate was 2.35 d/a. The summer maize SOS was advanced by 11 ± 1.5 d, the EOS was postponed by 6 ± 1.5 d, the LOS was extended by 17 ± 1.5 d, and the linear growth rate was 0.67 d/a. The LOSs of the Shendong, Binchang, and Yongcheng mining areas were extended during the study period. However, there were obvious differences between the Lu‘an mining area and the other three mining areas, with its SOS lagging, its EOS advancing, and its LOS shortening.

5. Conclusions

(1)

In the analysis of the monthly variation in the NDVI in the Yongcheng mining area during 1982–2022, the change in the NDVI in the experimental areas presented a significant bimodal shape, which reflected the phenology associated with the local farming practice of two-crop cultivation in northern China. The peaks in the second half of April and August corresponded to the growing season of the local winter wheat and summer maize, and the lowest value in June was related to the summer harvest of the winter wheat. The second decrease occurred at the end of September, and the lowest value was reached in the February of the next year. This decrease was closely related to the autumn harvest of the summer maize.

(2)

Research on the vegetation-growing season in the experimental areas during 1982–2022 showed that the four areas entered the SOS of the winter wheat from the end of February to the beginning of March and entered the EOS in the middle of June. The SOS of the summer maize was entered at the end of June and the EOS was entered in mid-October. The rate of advancement in the SOS was significantly higher than that in the EOS. Ultimately, both LOSs were extended. The extension rate of the LOS in the winter wheat and summer maize in the mining area was lower than that in the CK. Apparently, the changes in the phenology of the vegetation in the mining-disturbed area differed from those in the natural growing area.

(3)

The inter-annual variation analysis of the NDVI revealed that, regardless of whether they have been in the mining area, buffer zone, or control area, there has been an increase in the NDVI values. Notably, the vegetation activity under natural growth conditions (CK) was higher than that in the coal mining area (mining area), further confirming that coal mining activities indeed have a negative impact on vegetation growth.

(4)

The analysis of relationships between climatic parameters revealed that the NDVI in the Yongcheng mining area, buffers, and CK was significantly correlated with temperature and weakly correlated with precipitation. The vegetation in the mining area was more sensitive to temperature. In the context of global warming, the vegetation activity increased, but the growth rate of the NDVI in the mining areas declined. This result fully demonstrates that coal mining has had an impact on the growth of vegetation in the mining area.