Long-Term Effects of Fire Severity and Climatic Factors on Post-Forest-Fire Vegetation Recovery

Abstract

As a major disturbance to forest ecosystems, wildfires pose a serious threat to the ecological environment. Monitoring post-fire vegetation recovery is critical to quantifying the effects of wildfire on ecosystems and conducting forest resource management. Most previous studies have analyzed short-term (less than five years) post-fire recovery and limited the driving factors to temperature and precipitation. The lack of long-term and multi-faceted observational analyses has limited our understanding of the long-term effects of fire on vegetation recovery. This study utilized multi-source remote sensing data for a long time series analysis of post-fire vegetation recovery in China based on Google Earth Engine (GEE) cloud computing platform. Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EVI), Normalized Burn Ratio (NBR), and Normalized Difference Moisture Index (NDMI) were selected to quantify the low, moderate, and high severity of burned areas. Ridge Regression Model (RRM) was used to analyze the relationship between 15 driving factors and the vegetation regeneration process. The results show that it took at least 7–10 years for the vegetation index to recover to the pre-fire level after a forest fire. The recovery rate of high severity combustion areas was the fastest within the first two years. From the results of Ridge Regression, it came out that the overall fitting degree of the model with NDVI as the dependent variable was superior than that with EVI. The four variables of temperature, precipitation, soil temperature, and soil moisture were able to explain the change in more detail in vegetation indices. Our study enriches the research cases of global forest fires and vegetation recovery, provides a scientific basis for the sustainable development of forest ecosystems in China, and provides insight into environmental issues and resource management.

1. Introduction

Fire is a natural disturbance that occurs every few years in many forest ecosystems [1]. The distributions and characteristics of forest fuels are changing due to global warming, which impacts the occurrence of forest fires [2,3,4,5]. Fire affects atmospheric chemistry, forest ecosystem succession, carbon cycling, forests degradation, and land-use change [6,7,8,9,10]. Forest fires have occurred frequently in recent years, such as the California forest fires in the United States in October 2003 [11,12], the Australian forest fires in February 2001 and March 2002 [13,14,15], and the forest fires that occurred in the northeastern forest areas of China in 1987, 1992, and 1996 [16,17,18,19], etc., which have seriously affected the social, economic, and ecological benefits of forest areas [20,21,22,23]. On the other hand, fire also increases the diversity of ecosystem species and promotes natural regeneration [19]. Forests are renewable resources with a long growth cycle [24]. After severe fires, forests are often replaced by sparse woodland or scrub. If repeatedly burned, forests may become barren grassland or even bare ground [25,26].

The theory of vegetation recovery trends after forest fires is one of the hot spots in studies of forest fires. As early as the 1930s, Russian scholars began to study the effects of fire on the ecological environment [27]. By the 1950s, the United States and Canada began to pay attention to the impact of fire on various landscape types [21,28]. Currently, satellite remote sensing data are widely used in the study of forest fires [29]. Compared with traditional in situ monitoring methods, satellite remote sensing has the advantages of a larger observation scale and more continuity in time and spatial dimensions [30,31]. Therefore, remote sensing is becoming a primary method for studying forest fires. It helps to monitor the burned area and the severity of damage in real-time [32,33,34]. For example, forest fires can be detected by the brightness temperature of infrared (IR) radiation. In particular, the mid-infrared (MIR) and thermal infrared (TIR) bands are effective in detecting forest fires [35]. Since the spectral characteristics possessed by green vegetation are different from those of soil, water, and other features [36,37,38], hyperspectral remote sensing images can provide a good record of the spectral characteristics based on their chemistry and morphology [27], as well as provide good monitoring of the development, health, and growth conditions of the vegetation [39,40,41].

With the development of satellite observation technology, high-resolution optical image data have been used to study fire severity [42,43] and monitor the dynamic changes of vegetation after fire [12,44,45,46,47]. The study by Caccamo, G., R. Bradstock, L. Collins, T. Penman, and P. Watson of the vegetation recovery after eucalyptus fires in Australia indicated that spectral recovery was remarkably rapid in the first year after the fire, and the influence of severity was limited to the first two years after the fire [14]. Analysis of short-term (<5 years) post-fire vegetation recovery in Larix gmelinii (Larix gmelinii [Rupr.] Kuzen.) forests in northeastern China showed that climatic conditions and fire severity in the first growing season after fire play important roles in vegetation recovery [18]. The recovery rate of Normalized Burn Ratio (NBR) in the short term after a fire varies with time, severity, and forest type. Recovery rate usually reaches the highest shortly after the fire [32].

Previous studies have focused on short-term (<5 years) changes in vegetation greenness after forest fires [48,49]. The effects of soil structure and meteorological factors (temperature and precipitation) on vegetation recovery were also analyzed using ordinary least squares [12,50,51]. However, long-term vegetation recovery has seldom been studied. Although the distribution of fires at different temporal and spatial scales was summarized [52], the analysis of the number and area of fires were not comprehensive enough. The least-squares method is a widely used classical statistical algorithm. However, the negative impact of covariance on least squares estimation in a regression setting is inevitable. The presence of multiple covariance problems may lead to invalid least squares estimates [19,53,54].

Ridge Regression provides a way to solve the covariance problem without removing variables from the original set of independent variables [55,56]. Specifically, Ridge Regression makes the regression coefficients estimable by introducing a k-unit array. The regression results are more stable than simple linear regression [57]. The regression coefficients are more realistic by abandoning the unbiased nature of least squares at the expense of losing some information [58]. It is more tolerant to pathological data than the least-squares method.

In this study, we analyzed the vegetation recovery under different fire severities and vegetation types. Revealing the effects of fire on forest succession through long-term post-fire monitoring. Ridge Regression, which is more tolerant of pathological data and has more accurate regression models, was used to assess the impact of multiple driving factors on post-fire vegetation recovery. The results of this study can provide a scientific basis for forest fire management and forest resource protection departments.

The main research objectives of this paper are:

• To explore the pattern of forest recovery under different fire severities and vegetation types conditions through long-term observations;
• To explore the effect of each factor on post-fire recovery using Ridge Regression.

2. Materials and Methods

2.1. Study Areas

We used the global fire dataset (CMS_Global_Fire_Atlas_1642) to identify fire points (https://daac.ornl.gov/CMS/guides/CMS_Global_Fire_Atlas.html accessed on 27 July 2021). As shown in Figure 1a, the forest fires in China mainly distribute in the northern and southwestern regions. According to the “Regulations on Forest Fire Prevention” and related literature, we chose the study areas based on:

• Areas with a forest burned area greater than $10 \mathrm{k}{\mathrm{m}}^2$ to ensure that the recovery time would not be too fast;
• Fires that occurred between 2003 and 2010, to ensure a recovery study period of at least ten years;
• Areas in which no other fires occurred within the next five years, in the recovery process after this fire.

The western regions were not included in the analysis because of the small burned area and number of fires.

We first filtered the Chinese region by the fire dataset. Next, we filtered out forest fires. Finally, the final study areas were selected based on burned area and year. We selected four forest fires in three regions, i.e., northeast China (Zone 1, Zone 2), central China (Zone 3), and southwest China (Zone 4), as the study areas. In this study, we analyzed the recovery process of vegetation in burned areas from two aspects, different burn severities and different vegetation types. Detailed information on the study area is shown in

Table 1. Location, time, and burning area of fire points in China from 2003 to 2016.

Study Area	Location	Coordinates	Time	Burning Area $\left(\mathbf{k}{\mathbf{m}}^2\right)$	Main Vegetation Types
Zone 1	Hulun Buir	51.1312 N–121.839 E	2–20 May 2003	981.567	Deciduous Needleleaf Forests
					Deciduous Broadleaf Forests
					Mixed Forests
					Shrublands
					Grasslands
Zone 2	Da Hinggan Ling Prefecture	52.0021 N–126.125 E	30 April–16 June 2003	4596.004	Deciduous Needleleaf Forests
					Deciduous Broadleaf Forests
					Mixed Forests
					Grasslands
Zone 3	Shanxi	38.0021 N–113.329 E	20–25 April 2006	81.565	Shrublands
					Grasslands
Zone 4	Yunnan	25.681 N–100.376 E	10–21 March 2010	40.576	Mixed Forests
					Grasslands

.

Zone 1 and Zone 2 (Figure 1b–e) are the northernmost and highest latitude border regions in China. The area has a cold-temperate continental monsoon climate with rainy summers and winters and significant temperature differences [18]. The terrain consists of mesas, low mountains, hills, and intermountain basins. The average elevation is 573 m; the highest peak is 1528 m. The average annual temperatures in the south and north are 1 °C and −6 °C, and the precipitation is 442 mm and 240 mm, respectively. The natural vegetation type of the Da Hinggan Ling Prefecture is the cold-temperate coniferous and deciduous forest area of the Eurasian forest plant sub-region of the Pan-Arctic flora [16,59]. The forest cover rate is over 60%, with coniferous forests dominated by Larix gmelinii (Larix gmelinii [Rupr.] Kuzen.). The main species are Camphor pine (Pinus sylvestris var. mongholica Litv.), Red-barked spruce (Picea koraiensis Nakai), Birch (Betula platyphylla Suk.), Mongolian oak (Quercus mongolica Fisch. ex Ledeb), Populus (Populus davidiana Dode), etc. The Da Hinggan Ling Prefecture is one of the state-owned forest areas and the main area of natural forests in China. In the mountains above 1000 m of sea level, there are a large number of low-temperature plants, such as Pinus pumila (Pinus pumila [Pall.] Regel), Yuehua birch (Betula ermanii Cham.), Artemisia lagocephala (Artemisia lagocephala), Linnaea borealis (Linnaea borealis Linn.), etc. The slopes below 1000 m are widely covered with Cuckoo forests with a well-developed understory of shrubs. The Deciduous Broadleaf Forests are mainly distributed below 400 m above sea level. The forest composition is dominated by Mongolian oak (Quercus mongolica Fisch. ex Ledeb). Meadow plants are mainly found on gentle slopes of forest margins, interstices of forests, and broad valley edges. Zone 3, which corresponds to the 2006 forest fire (Figure 1d), occurs within Shanxi Province, with the burned area located at the border of Jinzhong City and Yangquan City. This area belongs to the warm-temperate semi-humid continental monsoon climate zone. The topography is high in the west and low in the east, with the highest point at 1804 m above sea level and the lowest end at only 350 m above sea level. The average annual temperature is generally 8–12 °C, and the yearly precipitation is 450–550 mm. Influenced by monsoons and complex topography, the four seasons are distinct, and the regional climate in the territory has noticeable vertical variations. The vegetation type in the study area is mainly shrubs and grasses due to topography, geomorphology, and climate. Most forests are artificially planted with Red pine (Pinus koraiensis Sieb. et Zucc.), Cypress (Cupressus funebris Endl.), Poplar (Populus tomentosa Carr), etc. The forest fire in Yunnan Province in 2010 (Zone 4) (Figure 1e) started at the junction of Dali City and Binchuan County. The topography of the study area is high in the northwest and low in the southeast. It is a northern subtropical plateau monsoon climate type with an average annual temperature of 14.9 °C. Due to the influence of topography and climate, the rainfall in the northwest is more significant than that in the southeast. The main tree species include Yunnan pine (Pinus yunnanensis Franch.), Huashan pine (Pinus armandii Franch.), Hemlock (Tsuga chinensis [Franch.] Pritz.), Fir (Abies fabri [Mast.] Craib), Horsetail fir (Phlegmariurus phlegmaria [L.] Holub), Simao pine (Pinus kesiya var. langbianensis), Cypress (Cupressus funebris Endl.), Camphor (Cinnamomum camphora [L.] presl), Tsubaki (Ailanthus altissima [Mill.] Swingle), Oak (Quercus L.) and so on. The study area is dominated by Mixed Forests and Grasslands. The annual temperature difference is slight, and the seasons are not very pronounced. The temperature at this site decreases with increasing altitude, and rainfall increases with increasing space altitude [60].

2.2. Data

Moderate-resolution imaging spectroradiometer (MODIS) products were selected to extract fire trails and were used as the primary data source for monitoring the dynamic changes of vegetation in fire trails. Three meteorological data products, the fifth generation ECMWF (European Centre for Medium-Range Weather Forecasts) Atmospheric Reanalysis of the global climate (ERA5) data, Terra Climate data, and Famine Early Warning Systems Network (FEWS NET) Land Data Assimilation System (FLDAS) data, were selected as the data sources for extracting the long time series of driving factors, and the specific information is shown in

Table 2. Dataset used in this study.

Product/Bands		Temporal and Spatial Resolution	Explain
MODIS	MCD64A1	500 m	Identify fire tracks in the study area and draw vector maps
	MOD14A1	1 km-8 d	The fire confidence is specified in pixels
	MOD13A1	500 m-16 d	Extraction and calculation of vegetation index (NDVI, EVI, NDMI)
	M0D09A1	500 m-8 d	Calculation the Normalized Burn Ratio (NBR) and the differenced Normalized Burn Ratio (dNBR)
	MCD13Q1	1 km	Land cover type classification
ERA5	temperature	2 m-1 monthly	Average air temperature at 2 m (daily mean)
	precipitation	m	Total precipitation (daily total)
	Soil temperature	K	The temperature of the soil in layer 1 (0–7 cm) of the ECMWF Integrated Forecasting System.
Terra Climate	Soil	mm	Soil moisture derived using a one-dimensional soil water balance model
	Aet	mm	Actual evapotranspiration derived using a one-dimensional soil water balance model
	Def	mm	Climate water deficit derived using a one-dimensional soil water balance model
	Pr	mm	Precipitation accumulation
	Tmmn	°C	Minimum temperature
	Tmmx	°C	Maximum temperature
FLDAS	Qg_tavg	$\raisebox{1ex}{$\mathrm{W}$}\!\left/ \!\raisebox{-1ex}{${\mathrm{m}}^2$}\right.$	Soil heat flux
	Qair_f_tavg	$\raisebox{1ex}{$\mathrm{k}\mathrm{g}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{k}\mathrm{g}$}\right.$	Specific humidity
	SoilMoi00_10cm	$\raisebox{1ex}{${\mathrm{m}}^3$}\!\left/ \!\raisebox{-1ex}{${\mathrm{m}}^3$}\right.$	Soil moisture (00–10 cm underground)
	SoilMoi10_40cm		Soil moisture (10–40 cm underground)
	SoilMoi40_100cm		Soil moisture (40–100 cm underground)
	SoilMoi100_200cm		Soil moisture (100–200 cm underground)

. We downloaded data through the Google Earth Engine (GEE) platform [61,62,63,64], and the projection coordinate system was set to “EPSG:32649”. The pre-processing of the source data was performed as follows. First, the images were cropped by the software ENVI (version 5.3. Colorado, USA: Exelis Visual Information Solutions). Next, the data were corrected using the “Band Math Batch” tool in ENVI to remove outliers and multiply by a scale factor of 0.0001 to make the reflectance values within (0,1). Finally, MATLAB (version R2021a. Massachusetts, USA: MathWorks) was used to calculate each vegetation index and synthesize the data with different temporal resolutions into annual data and export them.

2.3. Methods

Based on multi-source remote sensing data, this paper quantifies the long-term recovery process of vegetation on fire sites by monitoring post-fire vegetation recovery from forest fires for 10–20 years. The differences in vegetation recovery under different fire severities and vegetation types are also explored. Ridge Regression Method (RRM) was used to evaluate the influence of multiple driving factors on the recovery of fire sites and identify the key factors. The specific process is shown in Figure 2.

2.3.1. Vegetation Index

We used 8-day synthetic surface reflectance MODIS data at 500 m resolution to calculate indices in

Table 3. Selected remote sensing vegetation indices to detect fire severity in China.

Index	Expression	Description
Normalized Difference Vegetation Index $\left(\mathrm{NDVI}\right)$	$\mathrm{NDVI}=\frac{{\mathsf{\rho }}_{\mathrm{NIR}}-{\mathsf{\rho }}_{\mathrm{RED}}}{{\mathsf{\rho }}_{\mathrm{NIR}}+{\mathsf{\rho }}_{\mathrm{RED}}}$	NDVI and EVI are used to monitor changes in vegetation greenness (The values of NDVI and EVI are between −1 and 1, and the normal values of green vegetation)
Enhanced Vegetation Index $\left(\mathrm{EVI}\right)$	$\mathrm{EVI}=2.5\times \frac{({\mathsf{\rho }}_{\mathrm{NIR}}-{\mathsf{\rho }}_{\mathrm{RED}})}{{\mathsf{\rho }}_{\mathrm{NIR}}+6.0{\mathsf{\rho }}_{\mathrm{RED}}-7.5{\mathsf{\rho }}_{\mathrm{BLUE}}+1}$
Normalized Burn Ratio $\left(\mathrm{NBR}\right)$	$\mathrm{NBR}=\frac{{\mathsf{\rho }}_{\mathrm{NIR}}-{\mathsf{\rho }}_{\mathrm{SWIR}}}{{\mathsf{\rho }}_{\mathrm{NIR}}+{\mathsf{\rho }}_{\mathrm{SWIR}}}$	Monitor whether the fire occurs and identify the fire trace
Normalized Difference Moisture Index $\left(\mathrm{NDMI}\right)$	$\mathrm{NDMI}=\frac{{\mathsf{\rho }}_{\mathrm{NIR}}-{\mathsf{\rho }}_{\mathrm{MIR}}}{{\mathsf{\rho }}_{\mathrm{NIR}}+{\mathsf{\rho }}_{\mathrm{MIR}}}$	Extract the water content of the vegetation canopy
the differenced Normalized Burn Ratio $\left(\mathrm{dNBR}\right)$	$\mathrm{dNBR}={\mathrm{NBR}}_{\mathrm{pre}}-{\mathrm{NBR}}_{\mathrm{post}}$	Estimate the burn severity of the fire

. NDVI is the most widely used vegetation index indicating changes in vegetation cover, biomass, and ecosystem [65,66,67]. However, it can be oversaturated in areas with high above-ground biomass, such as tropical rainforests [68,69,70]. It is susceptible to soil background at low vegetation cover, which can cause changes in the NDVI time series patterns and thus result in some uncertainties in the parameter extraction [71,72,73]. Therefore the MODIS vegetation index product team developed a new enhanced vegetation index EVI. It mitigates the effects of atmospheric aerosol scattering and soil background radiation by introducing blue-band reflectance and mathematically transforming the calculation formula [28,74,75].

In this study, the differenced Normalized Burn Ratio (dNBR) was obtained by calculating the difference between the pre-fire and post-fire NBRs. We performed the k-means clustering analysis to obtain the thresholds and distributions of the three fire intensities, as detailed in

Table 4. Classification of different fire severities.

Study Area	Low Severity	Moderate Severity	High Severity
Zone 1	0 < dNBR < 0.2	0.2 ≤ dNBR < 0.4	dNBR ≥ 0.4
Zone 2	0.1 < dNBR < 0.16	0.16 ≤ dNBR < 0.3	dNBR ≥ 0.3
Zone 3	0.055 < dNBR < 0.145	0.145 ≤ dNBR < 0.22	dNBR ≥ 0.22
Zone 4	0.04 < dNBR < 0.12	0.12 ≤ dNBR < 0.24	dNBR ≥ 0.24

.

We monitor the changes in NBR and NDMI by calculating the year-over-year growth rate of the indicators and thus determine the occurrence of fires. The rate is the year-on-year growth rate, ${\mathrm{value}}_{\mathrm{i}}$ is the indicator value for the current period, and ${\mathrm{value}}_{\mathrm{i}-1}$ is the indicator value for the previous period.

$$\mathrm{rate}=\frac{{\mathrm{value}}_{\mathrm{i}}-{\mathrm{value}}_{\mathrm{i}-1}}{{\mathrm{value}}_{\mathrm{i}-1}}\times 100\%$$

(1)

The average EVI value of the three years before the fire disturbance was used as the average level of the vegetation index in this area. The recovery level of EVI each year after fire disturbance is calculated by the following formula:

$$\mathrm{recovery} \mathrm{level}=\frac{{\mathrm{EVI}}_{\mathrm{i}}}{{\mathrm{EVI}}_{\mathrm{mean}}}\times 100\%$$

(2)

where ${\mathrm{EVI}}_{\mathrm{i}}$ represents the EVI values for different fire severities in that year and ${\mathrm{EVI}}_{\mathrm{mean}}$ is the average EVI value for the three years before the fire.

2.3.2. Ridge Regression

Ridge Regression is one of the most frequently used regularization methods for regression analysis of ill-posed problems. It is commonly used to solve multicollinearity problems in multiple linear regression [54]. The least-squares method is widely used in regression analysis, which is an unbiased estimation. For a well-posed problem, $\mathrm{X}$ is usually column full rank i.e., $\mathrm{X}\mathsf{\theta }=\mathrm{y}$. Using the least-squares method, the loss function is defined as the square of the residuals, minimizing the loss function as ${‖\mathrm{X}\mathsf{\theta }-\mathrm{y}‖}^2$. The above optimization problem can be solved by the gradient descent method or directly by using the following Equation (3).

$$\mathsf{\theta }={({\mathrm{X}}^{\mathrm{T}}\mathrm{X})}^{-1}{\mathrm{X}}^{\mathrm{T}}\mathrm{y}$$

(3)

When $\mathrm{X}$ is not a column of full rank, or when the linear correlation between columns is relatively significant, the determinant of ${\mathrm{X}}^{\mathrm{T}}\mathrm{X}$ approaches 0, i.e., ${\mathrm{X}}^{\mathrm{T}}\mathrm{X}$ approaches the singular matrix, and the problem becomes an ill-posed problem. At this point, the error in calculating ${({\mathrm{X}}^{\mathrm{T}}\mathrm{X})}^{-1}$ will be large, and the traditional least-squares method lacks stability and reliability.

Adding the regularization term to the loss function transforms it into a fitness problem,

$${‖\mathrm{X}\mathsf{\theta }-\mathrm{y}‖}^2+{‖\mathrm{\Gamma }\mathsf{\theta }‖}^2$$

(4)

let $\mathrm{\Gamma }=\mathrm{kI}$, then

$$\mathsf{\theta }(\mathrm{k})={({\mathrm{X}}^{\mathrm{T}}\mathrm{X}+\mathrm{kI})}^{-1}{\mathrm{X}}^{\mathrm{T}}\mathrm{y}$$

(5)

where $\mathrm{I}$ is the unit matrix. As $\mathrm{k}$ increases, $\left|\mathsf{\theta }{(\mathrm{k})}_{\mathrm{i}}\right|$ tends to become smaller, and their deviation from the correct value ${\mathsf{\theta }}_{\mathrm{i}}$ becomes larger. $\mathsf{\theta }(\mathrm{k})$ tends to 0 as $\mathrm{k}$ tends to infinity. The trajectory of $\mathsf{\theta }(\mathrm{k})$ varying with $\mathrm{k}$ is called the ridge trace. Ridge Regression complements least squares regression, which loses unbiasedness in exchange for high numerical stability, resulting in higher computational accuracy [76]. Although Ridge Regression yields a more enormous sum of squared residuals than least squares regression, in this way, it is far more tolerant of pathological data than least squares [77].

Assuming that the data matrix $\mathrm{X}=\left\{{\mathrm{x}}_{\mathrm{ij}}\right\}$ of the independent variables is $\mathrm{n}\times \mathrm{p}$, the least-squares estimation seeks the coefficient $\mathsf{\beta }$ that minimizes the sum of squared residuals, i.e.,

$$({\widehat{\mathsf{\alpha }}}^{(\mathrm{ols})},{\widehat{\mathsf{\beta }}}^{(\mathrm{ols})})=\underset{(\mathsf{\alpha },\mathsf{\beta })}{\arg \min }{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}}{({\mathrm{y}}_{\mathrm{i}}-\mathsf{\alpha }-{\displaystyle \sum _{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{X}}_{\mathrm{ij}}{\mathsf{\beta }}_{\mathrm{j}}})}^2}$$

(6)

and the Ridge Regression requires an additional term $\mathsf{\lambda }{\displaystyle \sum _{\mathrm{j}=1}^{\mathrm{n}}{\mathsf{\beta }}_{\mathrm{j}}^2}$ to constrain the magnitude of the coefficients, such that the residual sum of squares is small and the coefficients are not inflated, i.e.,

$$({\widehat{\mathsf{\alpha }}}^{(\mathrm{ridge})},{\widehat{\mathsf{\beta }}}^{(\mathrm{ridge})})=\underset{(\mathsf{\alpha },\mathsf{\beta })}{\arg \min }{\displaystyle \sum _{\mathrm{i}=\mathrm{1}}^{\mathrm{n}}\left\{{({\mathrm{y}}_{\mathrm{i}}-\mathsf{\alpha }-{\displaystyle \sum _{\mathrm{j}=\mathrm{1}}^{\mathrm{n}}{\mathrm{X}}_{\mathrm{ij}}{\mathsf{\beta }}_{\mathrm{j}}})}^2+\mathsf{\lambda }{\displaystyle \sum _{\mathrm{j}=\mathrm{1}}^{\mathrm{n}}{\mathsf{\beta }}_{\mathrm{j}}^2}\right\}}$$

(7)

where $\widehat{\mathsf{\beta }}(\mathrm{k})={({\mathrm{X}}^{\mathrm{T}}\mathrm{X}+\mathrm{k}\mathrm{I})}^{-\mathrm{1}}{\mathrm{X}}^{\mathrm{T}}\mathrm{Y}$ is the Ridge Regression parameter estimate, $\mathrm{k}$ is the ridge parameter. The effect of multicollinearity can be reduced by the choice of the $\mathrm{k}$ value. Since the ridge parameter $\mathrm{k}$ is not uniquely chosen, the Ridge Regression parameter estimate $\widehat{\mathsf{\beta }}(\mathrm{k})$ is an estimation class of the regression parameter $\mathsf{\beta }$ with different values of $\mathrm{k}$ and different values of $\widehat{\mathsf{\beta }}(\mathrm{k})$. Ridge trace plots are used to determine $\mathrm{k}$ values. The component ${\widehat{\mathsf{\beta }}}_{\mathrm{i}}(\mathrm{k})$ of the ridge is estimated as a function of $\mathrm{k}$. When $\mathrm{k}$ varies in $\left[0,+\infty \right)$, the graph depicted by the plane rectangular coordinate system is the ridge trace plot. The principle of determining the value of $\mathrm{k}$ is to choose the smallest value of ${\mathrm{k}}_{\min }$ which makes the standardized regression coefficients of each variable stable. When $\mathrm{k}=0$ is equivalent to ordinary linear regression (OLS). The regression model is obtained by substituting the value of ${\mathrm{k}}_{\min }$ for the modeling.

In the regression analysis of various models, autocorrelation between independent variables is inevitable. If the multicollinearity between the independent variables can be reduced or eliminated, the accuracy and reliability of the index prediction can be improved. Multicollinearity refers to the existence of an approximately linear relationship between the independent variables of a multiple regression model. In this paper, the Variance Inflation Factor (VIF) is chosen to measure the covariance between the explanatory variables of the model [56,78]. The variance inflation factor is the reciprocal of tolerance $(\mathrm{1}-{\mathrm{R}}_{{}^{\mathrm{i}}}^2)$, where ${\mathrm{R}}_{\mathrm{i}}$ is the negative correlation coefficient of the independent variable ${\mathrm{x}}_{\mathrm{i}}$ for regression analysis of other independent variables.

$$\mathrm{VIF}=\frac{1}{1-{\mathrm{R}}_{\mathrm{i}}^2}$$

(8)

Usually, a VIF value greater than 10 indicates a severe covariance problem. The existence of the multicollinearity problem makes the accuracy of the estimation significantly reduced; the stability of the estimated value becomes poor. Even though the overall significance of the regression equation was high, some regression coefficients could not pass the significance test. The positive and negative signs of the regression coefficients are inverted, which cannot reasonably explain the regression equation, thus reducing the application value of the regression equation [55,79].

Common methods for overcoming multicollinearity among explanatory variables include increasing the sample size, stepwise regression, and Ridge Regression. However, stepwise regression methods solve the covariance problem by directly removing the independent variables from covariance. This may remove key variables resulting in a meaningless regression equation. Ridge Regression provides a method to address the problem without eliminating variables from the original set of independent variables [54]. Therefore, this paper will solve this problem by modifying the model using the Ridge Regression analysis method, which is more general and has higher accuracy in the regression equation [80].

The Ridge Regression is divided into two steps. First, from the ridge trace plot, the optimal $\mathrm{k}$ value is determined. Second, the $\mathrm{k}$ values are brought to the regression model. In this study, ridge trajectory plots were obtained using IBM SPSS Statistics software (version 24.0. Armonk, NY, USA: IBM Corp). The horizontal axis is the ridge parameter and the vertical axis is the standardized regression coefficient. Then, we selected 15 independent variables as input variables, including the average air temperature at 2 m (air_temperature_2m), total daily precipitation (precipitation), soil temperature at 0–7 cm (soil_temperature), monthly climatic precipitation (pr), climatic moisture gain/loss using a one-dimensional soil moisture balance model (def), actual evapotranspiration (aet), soil_moisture, minimum temperature (tmmn), maximum temperature (tmmx), soil moisture at four depths (SoilMoi100_200cm, SoilMoi40_100cm, SoilMoi10_40cm, SoilMoi00_10cm), soil heat flux (Qg), and specific humidity (Qair). Regression analysis was performed with NDVI and EVI as dependent variables, respectively. The output models were used to study the effects of 15 climatic factors, including monthly mean temperature and precipitation, on changes in vegetation indices.

3. Results

3.1. Forest Fire Vegetation Restoration

3.1.1. Forest Fire Identification

The MCD64A1 product contains per-pixel burned area and mass information to determine the burning date for each 500 m grid cell within each MODIS tile. The MOD14A1 product includes fire occurrence time (day/night), fire location, detection confidence, fire radiated power, and other layers describing fire pixel properties. We distinguish between fire, no fire, and no observation by Fire-Mask layers. Fire points are identified, and vector boundaries are mapped by combining the two products.

From Figure 3, it can be observed that the moderate and high burning areas in Zone 1 account for 90% of the fire area. The high burning areas in Zone 2 are mainly distributed in the eastern part of the whole study area, and the fire intensity gradually decreases to the southwest. Zone 3 and Zone 4 are relatively small in terms of the fire area and severity and do not vary significantly in the images.

For NBR, as shown in Figure 4a, the year-on-year growth rate fluctuated within the range of (−1~30) % in the two years before the fire. After the disturbance of the fire occurred (shaded part), the year-on-year growth rate of NBR became negative and became more prominent with the increase of the severity of the forest fire. The negative growth rate of NBR reached 90.22%, 73.64%, 50.27%, and 55.15% in the four high burning zones, respectively. The year-on-year growth rate of the moderate burning zone fire year in Zone 4 was 5.02%, and the negative growth rate of the moderate burning zone one year after the fire was 36.76%, indicating that this forest fire had a lagging effect on the moderate burning zone. The year-on-year growth rate in the high severity burn zone one year after the fire was much greater than that in the low and moderate severity burn zone, with year-on-year growth rates of 8.17 and 3.6 times (Zone 1); 2.96, and 2.04 times (Zone 2); 1.86 and 1.31 times (Zone 3); and 4.02 and 1.15 times (Zone 4) for the low and moderate severity burn zone, respectively. Then, in the absence of fire, the year-on-year growth rate gradually tends to 0 and remains stable with no abnormality.

NDMI is more effective than NDVI in extracting the moisture content of vegetation canopy; as shown in Figure 4b, the year-on-year growth rate fluctuated within (−5~50) % in the two years before the fire. Fires consume a large amount of water in the vegetation, after which the indicator decreases and becomes negative and becomes more prominent with increasing severity of forest fires. The NDMI values gradually recovered and leveled off two years after the fire. Then, the growth rate values fluctuated above and below the value of 0 without abnormalities in the absence of fire. Combining the two indices can determine whether a disturbance has occurred and the duration of the observed disturbance condition.

3.1.2. Effects of Different Fire Severities on Post-Fire Recovery

In the burned areas, vegetation index values and species richness increased substantially recovered over time. Figure 5 shows the time-series variation of NDVI and EVI for different fire intensities. Since the fire years are different, the post-fire monitoring time ranges are also different.

Under the three fire intensities, NDVI and EVI index values have similar trends. Still, NDVI values are 0.2 to 0.4 units higher than EVI values, and NDVI recovery time is longer than EVI overall. It can be seen from Figure 5 that both NDVI and EVI values drop below the normal fluctuation range after the fire, and the decline varies with the severity. The index values increased rapidly within two years after the fire (one year after the Zone 4 fire). Then the rate of increase slowed down, and gradually stabilized.

From the trend in Figure 5, the recovery rate of the high severity burning zone is the fastest, followed by the moderate severity burning area within one year after fire disturbance. The EVI increments in Zones 1–4 were 2.34, 1.99, 1.54, and 2.24 times higher than those in the low-severity burning zone and 1.23, 1.42, 1.16, and 1.27 times higher than those in the moderate-severity burning zone, respectively. The similar recovery trends in the two study areas of Zone 1 and Zone 2 are related because they are both located in the northeast and have identical vegetation type compositions. Combined with Figure 6, it can be seen that low-severity areas were always the first to recover to pre-fire levels, requiring only one to four years. Four of the study areas reaching normal EVI levels four years (107% of normal for Zone 1), two years (101% of normal for Zone 2), and one year (100% for both Zone 3 and Zone 4) after the fire, respectively.

For Zone 1, the change in NDVI growth (reduction) rate after the fire was in a smaller band than EVI (Figure 6). The negative growth rates of the low, moderate, and high fire zones as indicators of NDVI after the 2003 fire were 9.67%, 20.37%, and 32.95%, respectively. The growth rate of the moderate fire zone was 0.16%, almost no growth, and the growth rate of the high fire zone was only 4.61%, indicating a slow recovery rate; the growth rate only increased significantly in the fourth year after the fire (2007). In contrast, the negative growth rates of the low, moderate, and high burned areas were 20.87%, 40.49%, and 57.03%, respectively. The year-on-year growth rates were 9.54%, 23.24%, and 38.67% in the year after the fire, and the comparison with the growth rate of EVI shows that NDVI is a weaker indicator than EVI for monitoring the vegetation condition.

The recovery of EVI from Figure 6 shows that the recovery rate varies for different severity levels. The recovery of EVI for low to moderate burning zones reaches 100% in the second (Zone 3), sixth (Zone 4), seventh (Zone 2), and eighth (Zone 1) years after the fire, respectively. For heavily burned areas, recovery took longer, requiring eight to eleven years to reach pre-fire levels. Fire has a significant impact on the successional recovery of vegetation.

In the overall EVI recovery process, Zone 1 took eight years to reach 102.1% of the pre-fire three-year average and remained above 100% after that. Zone 2 took seven years to get 101.7% of the pre-fire level, followed by three years (2011–2013) of minor fluctuations but was still above 99.6% of the ordinary level overall. Zone 3 recovered to the normal level one year after the fire, and the subsequent EVI values have remained above 1.25 times the normal level, which positively influenced the vegetation recovery, indicating that the fire promoted vegetation regeneration. Zone 4 never reached the pre-fire level, and the best recovery was 97.37% of the pre-fire normal level, also indicating that this fire disturbance inhibited the vegetation regeneration, resulting in the decrease of EVI.

3.1.3. Effects of Different Vegetation Types on Post-Fire Recovery

Since land cover types in Zone 3 and Zone 4 are relatively simple, Zone 1 and Zone 2 were selected to analyze the restoration of different vegetation types. Five land cover types were extracted in the Northeast study area, namely Deciduous Needleleaf Forests, Deciduous Broadleaf Forests, Mixed Forests, Shrublands, and Grasslands.

Combined with Figure 7a, through 19 years, it can be seen that in Zone 1, Deciduous Needleleaf Forests and Deciduous Broadleaf Forests had the highest year-on-year growth rates in the first year after the fire, 32.23% and 32.02%, respectively. This was followed by grassland with 27.26%, while Mixed Forests and Shrublands had the slowest growth rates, 26.39%, and 24.21%, respectively. However, from the third year onwards, the growth rate slowed down and remained around 1% on average. Shrublands showed a negative growth rate, indicating a slowing down of the vegetation recovery process. In the 11th year after the fire, all five vegetation types recovered to their pre-fire mean levels. In Zone 2, the recovery was less volatile than in Zone 1. The fastest growth rate was in Mixed Forests (33.71%), followed by Deciduous Broadleaf Forests (27.31%) and Grasslands (24.61%) in the first year after the fire. The slowest growth rate was in Deciduous Needleleaf Forests (22.88%) and the EVI values of Deciduous Needleleaf Forests were kept at the lowest level, seven years after the fire (2010). The EVI values of all four land cover species returned to pre-fire levels seven years after the fire (2010) and continued to increase slowly during the following decade.

To further understand the influence of fire on the vegetation succession process, this paper also extracted the information of phenological parameters to monitor the change of vegetation type characteristics during the vegetation restoration process more intuitively. The essential level value, maximum value, and seasonal amplitude were selected as phenological change parameters. The parameter characteristics of typical land cover in the two study areas of Zone 1 and Zone 2 were obtained by screening the location information of specific land cover (

Table 5. Typical land cover phenological parameters.

Type	Baseval.	Max	Mean	Amplitude
Deciduous Needleleaf Forests	0.2643	0.4477	0.3501	0.1835
Deciduous Broadleaf Forests	0.3029	0.5137	0.3970	0.2111
Mixed Forests	0.2648	0.4449	0.3716	0.1801
Shrublands	0.2481	0.4100	0.3153	0.1763
Grasslands	0.2365	0.4248	0.3524	0.1885

).

Figure 8 shows the monthly (April to September) B-spline curves of EVI, with the red dotted line showing the fire disturbance year. Figure 8 shows that in 2003, the EVI values of Zone 1 and 2 both began to decline when the EVI value did not reach the maximum value, indicating that disturbance occurred.

As shown in

Table 6. Typical land cover phenological parameters of Zone 1.

Year	Baseval.	Max	Ampl.	Mean	Estimation of Main Vegetation Types
2002	0.1935	0.4667	0.2732	0.3577	Forest
2003	0.1739	0.2765	0.1026	0.2026	Grasslands
2004	0.1886	0.3188	0.1303	0.2536	Grasslands
2005	0.2117	0.3482	0.1364	0.2904	Grasslands
2006	0.2666	0.3550	0.0883	0.3080	Grasslands
2007	0.2760	0.3757	0.0997	0.3261	Grasslands
2008	0.2550	0.3617	0.1067	0.3109	Grasslands
2009	0.2695	0.4151	0.1456	0.3410	Grasslands
2010	0.2892	0.4117	0.1225	0.3484	Grasslands
2011	0.2695	0.4419	0.1723	0.3606	Grasslands
2012	0.2611	0.4462	0.1851	0.3509	Grasslands
2013	0.2646	0.4953	0.2307	0.3809	Forest
2014	0.2592	0.4817	0.2224	0.3697	Forest
2015	0.2749	0.4539	0.1790	0.3606	Forest
2016	0.2790	0.4613	0.1823	0.3708	Forest
2017	0.2788	0.4890	0.2101	0.3756	Forest
2018	0.2733	0.4807	0.2074	0.3691	Forest
2019	0.2624	0.4458	0.1834	0.3579	Forest
2020	0.2653	0.4872	0.2218	0.3723	Forest

and Figure 8, the maximum and the average values of EVI for Zone 1 in the first five years after the fire (2004–2008) were small. The fire was more intense and severely impacted the main vegetation types. From 2009 to 2012, the EVI maximum and mean values increased to some extent (the maximum value increased to more than 0.40, and the mean value increased to about 0.35.) After 2013, the mean value of EVI in the fire sites increased to more than 0.38, and the main vegetation types started to show forest characteristics. A relatively severe forest fire occurred in May 2003 at the Zone 1 site. After the fire, the site showed a nearly grassland-like ground cover and began to recover gradually. During the 19 years of vegetation recovery, the site went through a transitional stage from grassland to forest. Grassland dominated for about 10 years, and the process of vegetation recovery was slow.

From

Table 7. Typical land cover phenological parameters of Zone 2.

Year	Baseval.	Max	Ampl.	Mean	Estimation of Main Vegetation Types
2002	0.2449	0.5275	0.2827	0.3975	Forest
2003	0.2276	0.3487	0.1211	0.2745	Grasslands
2004	0.2447	0.4619	0.2172	0.3449	Grasslands
2005	0.2541	0.4251	0.1709	0.3346	Grasslands
2006	0.2837	0.4438	0.1600	0.3577	Grasslands
2007	0.2926	0.5170	0.2245	0.3968	Grasslands
2008	0.2712	0.4365	0.1653	0.3529	Grasslands
2009	0.2665	0.4630	0.1965	0.3564	Grasslands
2010	0.2750	0.5377	0.2627	0.3923	Forest
2011	0.2819	0.5063	0.2243	0.3912	Forest
2012	0.2850	0.5158	0.2307	0.3849	Forest
2013	0.2725	0.5012	0.2287	0.3733	Forest
2014	0.2715	0.5460	0.2746	0.4077	Forest
2015	0.2946	0.5190	0.2244	0.3950	Forest
2016	0.3081	0.5683	0.2602	0.4264	Forest
2017	0.3053	0.5578	0.2524	0.4162	Forest
2018	0.3016	0.5428	0.2413	0.4012	Forest
2019	0.2929	0.5591	0.2663	0.4196	Forest
2020	0.2905	0.5484	0.2578	0.4095	Forest

and Figure 8, it can be seen that the mean value of EVI in Zone 2 dropped to 0.27, and the maximum value was only 0.35 after the fire. From 2004, the vegetation recovery process became more apparent, and the maximum and mean values of EVI started to increase gradually. The mean and maximum values suddenly increased in 2007 to 0.5 and 0.39, respectively, then there was a slight decrease in the following two years. After 2010, the EVI maximum value reached above 0.5, and the mean value was more significant than 0.40 and remained stable, and began to show the characteristics of woodland. A relatively severe forest fire occurred in May 2003 at the Zone 2 site. After the fire, the site showed a nearly grassland-like ground cover and began to recover gradually. During the 19 years of vegetation recovery, the site went through a gradual transition from grassland to forest. The grassland stage lasted only six years, and the vegetation recovery process was rapid.

3.2. Analysis of Driving Factors of Forest Fire

Simple linear regression of the model was done using IBM SPSS Statistics software (version 24.0. Armonk, NY, USA: IBM Corp), and the covariance VIF test is shown in

Table 8. Multicollinearity test.

Independent Variable	Zone 1	Zone 2	Zone 3	Zone 4
air_temperature_2m	280.9687	441.9791	198.6482	48.5518
precipitation	5.5803	5.2305	6.7810	8.9000
soil_temperature	161.7993	217.5032	101.5689	30.9440
pr	11.8155	8.2197	18.0869	8.0774
def	4.6275	7.6822	36.3580	58.7483
aet	5.9457	5.4454	24.4846	25.9693
soil_moisture	3.8887	2.8730	3.4613	5.8624
tmmn	116.1371	156.3384	111.4930	22.2615
tmmx	103.9297	127.7386	68.0482	14.1496
SoilMoi100_200cm	1.9654	1.2187	1.5865	11.1703
SoilMoi40_100cm	2.8919	2.9092	4.0294	68.9484
SoilMoi10_40cm	3.0035	4.7844	4.1863	72.8044
SoilMoi00_10cm	4.7541	3.4573	5.6312	37.7164
Qg	5.0672	8.8173	7.2950	8.1072
Qair	33.2108	39.3646	32.3019	34.8219

. Many variables showed VIF values greater than 100, which means there is a severe covariance problem among the variables.

To reduce the effect of multicollinearity, the study was conducted using a Ridge Regression Model with better performance. Although the parameter estimates obtained by using Ridge Regression analysis are biased, the significance and stability of the parameter estimates are significantly higher than those of ordinary regression, and their likelihood of being close to the truth is greater.

From the ridge trace plots (Figure 9), it can be seen that the standardized regression coefficients of the independent variables tend to be stable after the optimal K value is selected. The regression coefficients and results for the four study areas are shown in

Table 9. Ridge Regression correlation coefficient summary, where the marked stars are the variables that passed the significance test (* p < 0.05 ** p < 0.01).

Dependent Variable	Zone 1		Zone 2		Zone 3		Zone 4
	NDVI	EVI	NDVI	EVI	NDVI	EVI	NDVI	EVI
constant	−3.803409 **	−2.088772 **	−2.882776 **	−2.068646 **	−2.136167	−0.419936	2.546419	−0.634158
air_temperature_2m	0.007401 **	0.003959 **	0.006699 **	0.004553 **	0.008268 **	0.003899 *	0.007125	0.005414
precipitation	0.047458	0.060939	−0.086581	−0.087635	−0.147393	−0.140648	−0.320704	−0.024361
soil_temperature	0.007086 **	0.003763 **	0.003380 *	0.002259 *	0.013212 **	0.009271 **	−0.016231	−0.003412
pr	0.000054	0.000230 **	0.000072	0.000257	0.000361	0.000393 *	−0.000696 *	−0.000243
def	−0.000034	0.000005	−0.000019	0.000003	0.000053	0.000055 **	−0.000027	−0.000003
aet	0.000028	0.000010	0.000031	0.000034	0.000017	0.000038	0.000055	0.000036
soil_ moisture	0.000002	−0.000004	−0.000031	−0.000039	−0.000007	0.000004	0.000022	−0.000002
tmmn	0.000542 **	0.000327 **	0.001066 **	0.000606 **	0.000749 *	0.000579 **	0.000232	0.000241
tmmx	0.000680 **	0.000284 **	0.000555 **	0.000302 **	0.000008	0.000113	0.001693	0.000290
SoilMoi100_200cm	0.491963 **	0.255826 **	0.602195 **	0.474445 *	−12.1318 **	−10.5677 **	0.744633 **	0.405073 **
SoilMoi40_100cm	0.052904	0.154634	0.015860	0.114421	0.359546	0.213440	0.062536	0.240833
SoilMoi10_40cm	−0.071292	−0.174804	−0.045349	−0.242848	0.042609	−0.022720	0.124141	0.086111
SoilMoi00_10cm	−0.231088	−0.081440	0.040574	0.051104	0.606331	0.300097	0.153566	0.111667
Qg	0.004082 **	0.004766 **	0.007151 **	0.007434 **	−0.00771 **	−0.000791	−0.02041 **	−0.005873
Qair	5.053745	7.735988 **	11.940298 **	12.668206 **	−0.281842	0.283235	−3.006526	0.951760
R-square	0.911133	0.852603	0.933542	0.88843	0.876618	0.820434	0.557328	0.675701
Adjusted R-square	0.898675	0.83194	0.924226	0.872789	0.859321	0.795261	0.495271	0.630238
F value, F (15,107)	73.14, p = 0.000	41.26, p = 0.000	100.2, p = 0.000	56.80, p = 0.000	50.68, p = 0.000	32.59, p = 0.000	8.98, p = 0.000	14.86, p = 0.000

. The marked stars are the variables that passed the significance test, $* p<0.05, ** p < 0.01$; the corresponding t-values are in parentheses. The variables that passed the significance test were retained, and the final Ridge Regression Model was obtained as shown in

Table 10. Ridge Regression results.

Study Area	Dependent Variable	Adjusted R-Square	K	Ridge Regression Equation
Zone 1	NDVI	0.898675	0.08	$\begin{array}{ll}\mathrm{NDVI}=& 0.007{\mathrm{X}}_1+0.007{\mathrm{X}}_2+0.001{\mathrm{X}}_3+0.001{\mathrm{X}}_4+0.492{\mathrm{X}}_5+\\ & 0.004{\mathrm{X}}_6-3.803\end{array}$
	EVI	0.83194	0.22	$\begin{array}{ll}\mathrm{EVI}=& 0.004{\mathrm{X}}_1+0.004{\mathrm{X}}_2+0.000327{\mathrm{X}}_3+0.000284{\mathrm{X}}_4+\\ & 0.256{\mathrm{X}}_5+0.005{\mathrm{X}}_6+7.736{\mathrm{X}}_7+0.0002{\mathrm{X}}_8-2.089\end{array}$
Zone 2	NDVI	0.924226	0.03	$\begin{array}{ll}\mathrm{NDVI}=& 0.007{\mathrm{X}}_1+0.003{\mathrm{X}}_2+0.001{\mathrm{X}}_3+0.001{\mathrm{X}}_4+0.602{\mathrm{X}}_5+\\ & 0.007{\mathrm{X}}_6+11.940{\mathrm{X}}_7-2.883\end{array}$
	EVI	0.872789	0.05	$\begin{array}{ll}\mathrm{EVI}=& 0.005{\mathrm{X}}_1+0.002{\mathrm{X}}_2+0.001{\mathrm{X}}_3+0.0003{\mathrm{X}}_4+0.474{\mathrm{X}}_5+\\ & 0.007{\mathrm{X}}_6+12.668{\mathrm{X}}_7-2.069\end{array}$
Zone 3	NDVI	0.859321	0.02	$\begin{array}{ll}\mathrm{NDVI}=& 0.008{\mathrm{X}}_1+0.013{\mathrm{X}}_2+0.001{\mathrm{X}}_3-12.132{\mathrm{X}}_5-\\ & 0.008{\mathrm{X}}_6-2.136\end{array}$
	EVI	0.795261	0.04	$\begin{array}{ll}\mathrm{EVI}=& 0.004{\mathrm{X}}_1+0.009{\mathrm{X}}_2+0.001{\mathrm{X}}_3-10.568{\mathrm{X}}_5+0.0004{\mathrm{X}}_8+\\ & 0.0001{\mathrm{X}}_9-0.420\end{array}$
Zone 4	NDVI	0.495271	0.03	$\mathrm{NDVI}=0.745{\mathrm{X}}_5-0.020{\mathrm{X}}_6-0.0007{\mathrm{X}}_8+2.546$
	EVI	0.630238	0.04	$\mathrm{EVI}=0.405{\mathrm{X}}_5-0.634$

, where, ${\mathrm{X}}_1$: temperature, ${\mathrm{X}}_2$: soil_temp, ${\mathrm{X}}_3$: tmmn, ${\mathrm{X}}_4$: tmmx, ${\mathrm{X}}_5$: SoilMoi100_200cm, ${\mathrm{X}}_6$: Qg, ${\mathrm{X}}_7$: Qair, ${\mathrm{X}}_8$: pr, ${\mathrm{X}}_9$: def.

In Zone 1, the K value for Ridge Regression analysis with NDVI as the dependent variable was taken as 0.08 and the model ${\mathrm{R}}^2=0.911$, implying that 15 independent variables such as air_temperature_2m can explain 91.11% of the variation in NDVI. When EVI is the dependent variable, the model ${\mathrm{R}}^2=0.853$, implies that air_temperature_2m and other 15 variables can explain 85.26% of the variation in EVI. Differences in the equations of different dependent variables in the same study area: the EVI equation in Zone 1 increased the variables pr and Qair. In Zone 2 study area, NDVI, EVI remained the same; Zone 3𠆒s EVI equation increased the variables pr, def, Qg; Zone 4’s EVI equation decreased the variable temperature, pr, soilmoi40_100cm, and Qg, and the variables SoilMoi100_200cm were added to the Zone 4’s EVI equation.

In the regression model, the goodness of fit ${\mathrm{R}}^2$ of the equation with NDVI as the dependent variable was better than the ${\mathrm{R}}^2$ obtained with EVI as the dependent variable, except for Zone 4. The six variables ${\mathrm{X}}_1$–${\mathrm{X}}_6$ in the equation had a greater influence on the vegetation index, indicating that temperature, precipitation, soil temperature, soil water content, and soil heat flux explained more of the variation in the vegetation index. In general, regression coefficients are mainly used for prediction, while standard regression coefficients indicate the relative importance of each independent variable under the same conditions. The degree of influence of each factor on vegetation indices (NDVI, EVI) is shown in Figure 10.

It is clear that among the nine variables that passed the significance test after Ridge Regression, mean monthly temperature (${\mathrm{X}}_1$), soil surface temperature (${\mathrm{X}}_2$), and maximum temperature (${\mathrm{X}}_4$) have a higher ability to influence changes in NDVI in the Zone 1, all reaching 20% of overall importance. In contrast, soil heat flux (${\mathrm{X}}_6$) has the weakest ability to control, with only 9% importance. For Zone 1, the role of EVI variation is relatively average, and the most significant influence is also temperature (${\mathrm{X}}_1$, ${\mathrm{X}}_2$) accounting for 15%. For the study area of Zone 2, the minimum temperature (${\mathrm{X}}_3$) can explain 27% of the NDVI variation with the greatest influence; the minimum temperature (${\mathrm{X}}_3$), soil heat flux (${\mathrm{X}}_6$), and specific humidity (${\mathrm{X}}_7$) can explain 20% of the EVI variation, respectively. Zone 3 had relatively few influences, with 27% positive influence of soil temperature (${\mathrm{X}}_2$) on NDVI and 17.8% negative influence on soil heat flux (${\mathrm{X}}_6$). The regression equation for Zone 4 is not shown in the bar chart because of the small number of variables included in the regression equation.

During post-fire recovery, managers can target human intervention to create a natural environment suitable for vegetation recovery and promote forest recovery according to the magnitude of the impact of different climatic variables on the vegetation index.

4. Discussion

Every year there are more than 200,000 forest fires worldwide, large numbers of hectares of forest are burned in China. Forest fires have a severe impact on the ecological environment [77]. In this study, we analyzed the effect of fire severity on the recovery process of vegetation by comparing the recovery rates of vegetation indices under different fire intensities.

The experimental results showed that the areas with the highest fire severity had the fastest rate of vegetation recovery within the beginning two years after the fire, which is consistent with previous studies [15,81]. The burned areas could recover to pre-fire levels over seven years. Still, due to the lack of field data, it was hard to determine whether the disturbed vegetation had fully recovered. The results of this paper also demonstrate that fire promotes the growth and development of vegetation on fire sites, keeping the vegetation index values above pre-fire levels [25,35,82]. This happens due to the fact that fire burned many tall trees and let in more sunlight in low-lying areas, the seeds of some low-lying herbaceous plants were cracked through the fire, promoting plant growth [4,30]. In addition, previous studies have shown that full recovery of vegetation after severe forest fires takes more than 20 years [38]. To gain a deeper understanding of the links between ecosystem processes and fire impacts, more areas need to be studied. In the future, more efforts are required to explore long-term competitive relationships and recovery processes among different species in burned areas using higher-resolution remote sensing data, or through stratified field sampling.

Post-fire vegetation regeneration is controlled by meteorological conditions, and recovery is more rapid in spring and summer [3,22]. The adaptation of different vegetation to fire varies is also related to temperature and soil conditions [48]. The resolution of the climatic data used in this study was relatively coarse compared to the Landsat and topographic data. This may lead to a decrease in the explanatory power of the driving factors of vegetation recovery. Due to multicollinearity among different climatic data, the analysis could not be performed by the OLS regression method [26,42,54]. Therefore, we tried to build a model using the Ridge Regression Method (RRM), which is more tolerant of pathological data. The accuracy of the regression model was improved without reducing the variables [54]. In addition, we obtained four critical driving factors, namely temperature, precipitation, soil temperature, and soil moisture. After the fire, according to the actual condition of the fire site, we can create favorable conditions for the recovery of vegetation by increasing soil moisture and artificial rainfall.

Our results contribute to a better understanding of the interactions between climate, disturbance, and vegetation dynamics. We believe that the application of forest restoration using different indicators will allow us to evaluate the progress of forest restoration from different perspectives. It will potentially help forest management authorities to make scientific and rational arrangement. This study also has the following limitations. First, due to the lack of field measurement data, we could not determine the main reason for the recovery of NDVI and EVI to pre-fire levels. It may be caused by the recovery of fire-disturbed vegetation or the role of developing vegetation. Second, the low resolution of meteorological data selected for this study may affect the explanatory power of meteorological factors for small study area. Finally, the number of study samples in this study is small, and the conclusions obtained for post-fire recovery do not apply to the whole ecosystem. Still, the research method can be used for studies in other areas. Although this paper has made preliminary achievements in the analysis of factors affecting vegetation recovery, there are still many tasks that need to be done further. The key ones are briefly discussed as follows: First, how to find a feasible method for proper planning of post-fire vegetation recovery from the macro-level based on the results of Ridge Regression analysis. Second, there is a need to study on a larger spatial scale to summarize more comprehensive statistical patterns of the post-fire vegetation restoration process. Finally, a more comprehensive prediction can be obtained by using a combination of different statistical models.

5. Conclusions

The study of post-fire recovery is a crucial determinant of building re-establishment of ecosystem services. We studied the long-term recovery processes of four forest fires in northeastern, central, and southwestern China. Considering the importance of post-fire climatic conditions in explaining the recovery process, we used the highest possible spatial and temporal resolution data for monitoring. Post-fire methods and rates were assessed by long time series analysis of multivariate remote sensing data. In combination with Ridge Regression, the covariance between variables was resolved, and attribution analysis was performed to obtain the main driving factors. This study showed that different severity levels had a significant effect on the recovery rate of vegetation, with higher severity recovering faster in the short term, and the recovery rate gradually leveling off in the middle and late stages of fire. The four main factors affecting vegetation recovery obtained from the Ridge Regression analysis were temperature, precipitation, soil temperature, and soil moisture. These results improve our understanding of vegetation dynamics after forest fires. Combining these methods makes it possible to assess the different stages of forest succession after the fire. However, additional studies are needed to include other potentially driving factors (soil properties, etc.) to obtain a more comprehensive investigation. Our results contribute to the understanding of the interactions between climate, disturbance, and vegetation dynamics. It can provide a scientific basis for forest management, post-fire ecosystem restoration, and risk assessment. Future research directions can monitor forest fire recovery on a larger scale, combine more climate products, and explore new methods to obtain more extensive results.