Quantifying Forest Fire and Post-Fire Vegetation Recovery in the Daxin’anling Area of Northeastern China Using Landsat Time-Series Data and Machine Learning

: Many post-ﬁre on-site factors, including ﬁre severity, management strategies, topography, and local climate, are concerns for forest managers and recovery ecologists to formulate forest vegetation recovery plans in response to climate change. We used the Vegetation Change Tracker (VCT) algorithm to map forest disturbance in the Daxing’anling area, Northeastern China, from 1987 to 2016. A support vector machine (SVM) classiﬁer and historical ﬁre records were used to separate burned patches from disturbance patches obtained from VCT. Afterward, stepwise multiple linear regression (SMLR), SVM, and random forest (RF) were applied to assess the statistical relationships between vegetation recovery characteristics and various inﬂuential factors. The results indicated that the forest disturbance events obtained from VCT had high spatial accuracy, ranging from 70% to 86% for most years. The overall accuracy of the annual ﬁre patches extracted from the proposed VCT-SVM algorithm was over 92%. The modeling accuracy of post-ﬁre vegetation recovery was excellent, and the validation results conﬁrmed that the RF algorithm provided better prediction accuracy than SVM and SMLR. In conclusion, topographic variables (e.g., elevation) and meteorological variables (e.g., the post-ﬁre annual precipitation in the second year, the post-ﬁre average relative humidity in the ﬁfth year, and the post-ﬁre extreme maximum temperature in the third year) jointly affect vegetation recovery in this cold temperate continental monsoon climate region.


Introduction
Fire is a natural disturbance in ecosystems and promotes diversity and natural regeneration [1]. However, warm and dry conditions increase wildfire activity [2]. At the global scale, fire is a critical factor affecting atmospheric chemistry [3], forest ecosystems succession and health [4], carbon budgets, and land use transformation [5]. At the regional scale, fire affects the socioeconomic system, soil erosion, and the hydrological cycle [6]. The diversity of fire impacts requires a comprehensive assessment of fire vulnerability of different systems [7], including the analysis of the potential damage of fire on society and ecological systems. Additionally, knowing when and where fire occurs is of great significance for forest management, carbon cycle prediction after fire events, and local or regional scale climate effects [5]. In turn, vegetation recovery influences the climate and forest management [8]. Post-fire recovery can affect surface radiation balance, carbon budgets, and water balance by changing the surface roughness, surface albedo, soil moisture, erosion, and evapotranspiration, thus affecting the climate [9][10][11]. In addition, vegetation change also affects the re-establishment of ecological functions, species habitat recovery, and biota

Study Area
The Daxing'anling area, which covers 83,000 km 2 , is located in northeastern China at latitudes of 50 • 10 to 53 • 33 N and longitudes of 121 • 12 to 127 • 00 E (Figure 1). The area is located at an average elevation of 573 m in a cold temperate zone and has a continental monsoon climate with long, cold winters and hot, short summers. According to the meteorological administration of China, the daily mean temperature is −3 • C, the minimum temperature is −48 • C, and the maximum temperature is 36 • C. The annual rainfall is 500 mm. The minimum average monthly rainfall occurs typically in April and June, leading to a relatively dry period. The details on the Daxing'anling area were obtained from the Daxing'anling District Administration Office (http://www.dxal.gov.cn/ (accessed on 24 December 2020)).

Figure 1.
Location of the study area in Northern China. The colored tiles on the right show the Landsat footprints (p119r024, p119r025, p120r023, p120r024, p120r025, p121r023, p121r022, p121r025, p122r023, p122r024, and p123r025), and the black polygon is the exact study site.

Data and Preprocessing
Annual Landsat Thematic Mapper (TM)/Enhanced Thematic Mapper (ETM+)/ Operational Land Imager (OLI) observations during 1987-2016 were used in this study. The TM/ETM+/OLI images were geometrically and atmospherically corrected by the USGS EROS data center using the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) in batch mode [33] or the Landsat 8 Surface Reflectance Code (LaSRC) system per customer requests when placing orders. The surface reflectance images were then freely downloaded from the USGS EROS Center Science Processing Architecture (ESPA) (https://espa.cr.usgs.gov/ (accessed on 24 December 2020)). The images had to conform to the following two rules to ensure sufficient accuracy of the forest change detection: (1) the dates of the image acquisition had to fall into the growing season (mid-June to mid-September) and (2) the images had to have low or no cloud cover. For images with a relatively high cloud cover in a given year, we collected other images from the same path/row tile in the same growing season to composite a new image if the clouds in the images were not spatially coincident. Additionally, the Shuttle Radar Topography Mission 3 (SRTM3) digital elevation model (DEM) dataset was downloaded from the Global Land Cover Facility (GLCF) at the University of Maryland to support the analysis.
The fire data in point format used in this study was provided by the Shenyang Institute of Applied Ecology at the Chinese Academy of Sciences. The data consisted of all fire events recorded during 1967-2004 in the Daxing'anling region, including the center latitude and longitude coordinates of the fire points, the estimated burned area, and the burn severity. The fire point data after 2004 were obtained from the "13th Five-Year Plan" reports compiled by the State Forestry Administration of China and the China Ecological Remote Sensing Information Service Network (http://rsapp.nsmc.org.cn (accessed on 24 December 2020)). Additionally, from 12 August to 17 August 2017, our team conducted a field investigation of the vegetation recovery status and progress of a catastrophic forest fire that occurred on 6 May, 1987. We obtained the locations of fire points in 1987 and interviewed the local people to glean information on vegetation recovery methods. The fire point data after 2004 and the field investigation data included the geographical coordinates and the year of fire occurrence. Four large fires that occurred in 1987, 2000, 2006, and 2010 were selected to train the fire extraction model and vegetation recovery algorithms because of the large burned areas, high burn severity, and significant ecological effects.
The meteorological data were the annual climate data set of China, provided by the China meteorological data network (http://data.cma.cn/ (accessed on 24 December 2020)). A detailed description of the meteorological data is listed in Table 2. We only downloaded observations that included the annual average temperature, the extreme maximum temperature, the extreme minimum temperature, the average maximum temperature, the average minimum temperature, the annual precipitation, and the average relative moisture from the 33 meteorological stations in Heilongjiang province. The discrete station observations were interpolated to create raster layers using the co-kriging interpolation algorithm, which is based on the theory of variograms and structure analysis to obtain unbiased estimates of the regionalized variables in a limited area [34]. In the interpolation, elevation was used as a covariable of the meteorological factors, and GS+ software was used to conduct semi-variance analysis of the meteorological station data. The optimal semi-variogram function was an exponential function, with an R-squared (R 2 ) value of above 0.9.

VCT-Based Forest Disturbance Mapping and Validation
We used the VCT algorithm to create forest disturbance time series in the study area from 1987 to 2016. The VCT algorithm considers the spectral-temporal characteristics of land cover and forest change; individual image analysis and time-series analysis are conducted. The output of the VCT algorithm includes the years of forest disturbance, forest disturbance maps, and the disturbance magnitude. The integrated forest z-score (IFZ) was used in the VCT algorithm to track forest changes. The IFZ value is defined by integrating FZ i over the spectral bands of the multi-spectral images. The calculation formulas of the IFZ and the FZ i are as follows: where b pi represent the value of a given pixel of band i, and b i , and SD i represent the mean and standard deviation of the spectral value of the forest pixel of each band i, respectively. NB represents the number of bands, which was 3 in this study (Landsat TM/ETM+ images use bands 3, 5, 7; Landsat OLI images use bands 4, 6, 7). The IFZ is an inverse measure of the likelihood of each pixel being a forest pixel. When the IFZ value of a pixel shows sustained relative stability and is lower than 0.3 (dense forests) or 0.2 (sparse forests), the pixel is identified as persisting forest. Noise observations are likely to result in an IFZ peak (a high IFZ value preceded and followed by low IFZ values), whereas disturbance observations are likely to result in consecutive high IFZ values. Consecutive high IFZ values were used to avoid noise. Forest disturbance in a year is usually accompanied by a sharp increase in the IFZ value, whereas vegetation recovery was defined as the gradual reduction in IFZ values. The disturbance maps could provide information on where and when disturbances had occurred; thus, we focused on the disturbance maps in this research. The background values in the disturbance maps of the VCT were removed. Persisting forest and previously disturbed forest that was spectrally similar to the forest in this year were defined as forests. The remaining categories were defined as non-forest, including disturbed in this year, persisting non-forest, persisting water, and post-disturbance non-forest. Additionally, we adopted the accuracy validation method proposed by Li et al. [34] to assess the accuracy of the forest disturbance maps using a spatial agreement index. Specifically, we used a spatial agreement index (SAI) based on visual interpretation and Google Earth high-resolution images to verify the accuracy of the VCT forest disturbance products. We randomly selected three 3 × 3 km squares in the study area as the validation area ( Figure 2) to derive the SAI statistics using Equation (3). In this method, S v is the disturbed area detected by VCT, S g is the disturbed area determined by visual interpretations, and S cd is the overlap area of the two detected regions. We superimposed the annual maps of aggregated forest and non-forest obtained from VCT in ArcGIS to obtain a map of the annual vegetation recovery status.

Burned Area Extraction Using the VCT Products and SVM Algorithm
We combined the SVM classifier and the disturbances detected by the VCT to extract the burned areas in 1987, 2000, 2006, and 2010. We used the SVM classifier because this algorithm has the advantages of seeking the best classification solution and dealing with a small number of samples [35]. Fires that occurred between 2012 and 2016 were excluded from the analysis because in the subsequent vegetation recovery modeling, we used the NDVI values of the fire patches in the fifth year after the fire to characterize the vegetation recovery status. Thus, the large fire events extracted in the current work included but were not limited to the "May 6th" fire in 1987, the Huma County fire in 2000, the Kandu River fire in the Songling area in 2006, and the Huzhong fire in 2010.
where ρ N IR , ρ Red ,ρ Blue , ρ SW IR1 , and ρ SW IR2 are the surface reflectance in the Landsat TM/ETM+ bands 4, 3, 1, 5, and 7 and the Landsat OLI+ bands 5, 4, 2, 6, and 7. We used the fire point data in 1987, 2000, 2006, and 2010 and the coordinates of the fire point data and grew burned pixels from the centers of 10-15 fire patches. The resulting number of burned pixels per year was around 500 (depending on the burned area in each year) ( Table 1); these pixels acted as the training data for the burned area. Next, the non-burned pixels were identified from the false-color composites of the Landsat bands using visual interpretation, including forest, grassland, cropland, logging sites, built-up, and water. The number of non-burned pixels was nearly equal to that of the burned pixels. The identified pixels were then randomly divided into two groups; one group (70% of the pixels) was used for SVM classifier training, and the other group (30% of the pixels) was used for classification accuracy validation. In the SVM classification, the Gaussian radial basis function was regarded as the kernel function, the cost parameter C was 10 (representing a penalty for miss-classification errors), and γ was 50. Finally, a pixel was considered a fire disturbance pixel if it was labeled as a disturbance pixel by the VCT algorithm and was also classified as a burned pixel by the SVM. Thus, we were able to use fire points in each year of the Landsat time-series as training sets, and the SVM algorithm was used to extract the burned patches from the VCT-detected disturbances. As a result, the final dataset of forest fire patches in the four years was created by combining the VCT products and SVM algorithm. We used the remaining 30% of the pixels to create a confusion matrix to determine the overall accuracy, user's accuracy, producer's accuracy, and Kappa coefficient to quantify the reliability of the classifications.

Post-Fire Vegetation Recovery Modeling
In this study, we focused on fire disturbances that occurred in the Daxing'anling area. The NDVI is an efficient vegetation index in characterizing post-fire vegetation recovery [42]. Previous studies indicated that the NBR could be used to characterize recovery both in the short-term and long-term [43,44]; however, NDVI has the characteristics of rapid saturation and has been primarily used for shorter-term recovery assessments [45]. The increase in the post-disturbance NDVI value may be related to an increase in the cover of trees, shrubs, and herbaceous vegetation [30]. We focused on the early stages of post-fire vegetation recovery to reduce potential saturation effects. The NDVI value of the burned area in the first year was used as an independent variable to reflect the initial vegetation condition for possible vegetation recovery, whereas the NDVI values of the fifth year were used as the dependent variable to reflect the vegetative condition after recovery. Other independent variables included many topographic features (elevation, slope) [30], climatic factors (derived from meteorological station interpolation data) and fire severity (dNBR) ( Table 2). Meteorological factors included the annual average temperature, precipitation, and moisture [46]. Three methods were selected to model post-fire vegetation recovery, including stepwise multiple linear regression (SMLR), SVM, and RF. Unlike the SMLR model, RF and SVM are non-parametric machine-learning algorithms that can deal with complex nonlinear relationships. RF has the advantages of determining variable importance, assessing the robustness of data reduction, high accuracy, and low sensitivity to parameter adjustments [47]. SVM is particularly well suited for a small sample size.
The SMLR was implemented by adding an important variable to or eliminating an unimportant variable from a collection of explanatory variables based on predetermined criteria. Typically, this step is performed using a series of F-or t-tests to evaluate the influence of a particular variable on the dependent variable [48]. This process was repeated twice until no new variables could be introduced into the model.
The RF algorithm has been widely used in data mining, ecology, and other disciplines [49]. An RF typically consists of a large number of regression trees, each of which has a set of constraints that are hierarchically organized and applied from the roots of the tree to the leaves. An RF starts with many randomly selected bootstrap samples that are iteratively replaced by those from the original training data set. Multiple regression trees are then generated to form an RF, and the final prediction result is obtained by averaging the predictions from all regression trees [50]. In the current analysis, we first defined two variable importance measures based on different response types before creating the RF.
(1) After permutating each predictor variable randomly, the out-of-bag (OOB) error was used to compute the reduction in predictive performance. (2) The reduction in the node impurity can be quantified by the Gini index. Thus, the contribution of each predictor variable to the modeling performance was assessed and ranked. The percent increase in the mean square error (MSE) (PercentIncMSE) and the increase in the NodePurity (IncN-odePurity) depend on the importance of the predictor variable. After multiple tries, we determined the following parameter values: The number of random regression trees (ntree) was 500, the number of segmentation variables (mtry) was 3, and the minimum sample size of the terminal node (node size) was 5.
The objective of the SVM is to minimize the prediction error and maximize the flatness of the function [51]. This trade-off is controlled by setting the error penalty factor C. The parameter ε is a zone defined around the regression function, and the error in the range of ±ε is ignored. If all the training tuples are in a zone with a width of 2ε, a function is output at the center of the flat channel containing all the training tuples, in which case the total error is 0. Therefore, ε controls the fitting degree between the function and the training tuples [52]. We used the same suite of independent variables derived from the stepwise regression analysis mentioned above in the SVM modeling and RF modeling to allow for the comparison of the models. A Gaussian kernel function was used; the error penalty factor C was 10, and ε-svr was the regression algorithm in the SVM regression.
We randomly selected around 100 pixels (the specific number depended on the recovery area) from the burned areas mapped by the VCT-SVM algorithm in each year and extracted the NDVI value and several impact factors for the 5 years (when the image was missing, we will use the image of next year) to fit the models. It is worth noting that collinearity usually exists among the meteorological variables. Thus, the selection of predictors used in the RF, SVM and SMLR models followed two steps: (1) we selected the predictors using the importance function in the RF model; (2) we eliminated the meteorological factors causing collinearity using SMLR [53]. An accuracy assessment was performed to assess the model's prediction performance. In this study, we used 70% of the selected pixels for modeling and the remaining 30% for validation. We used the validation dataset to derive the mean absolute error (MAE), the coefficient of determination (R 2 ), and the root mean square error (RMSE) to evaluate the prediction performance of the models.
where y i is the observed NDVI value,ŷ i is the value predicted by the model, y is the average value of the observed samples, and N is the number of samples. Table 3 shows the validation statistics derived from the three validation quadrats. It was observed that forest disturbances frequently occurred during the study period, and no disturbance was observed in seven, six, and four years in the study period in the eastern, middle, and western quadrat, respectively. Most of the spatial agreement ranged from 70% to 86%; the highest spatial agreement (92.35%) was observed in the eastern quadrat in 1987, and the lowest (60.03%) was observed in the western quadrat in 2010. Additionally, a very large disturbance patch was detected in 1987 ( Figure 2); it was the site of "the May 6th" catastrophic fire in 1987 that led to 193 deaths, one million hectares of burned land (both inside and outside the study area), and 50,000 homeless people.  Figure 3 shows the dynamics of the disturbance and recovery events (disturbed forest area, disturbance rate, vegetation recovery area, and recovery rate) during 1987-2016 in the Daxing'anling region. The largest forest disturbance area was observed in 1987, with a total forest disturbance area of 362,288.8 ha and a disturbance rate of 4.33%; this result was primarily attributed to the 1987 "May 6th" catastrophic fire. The second-and third-largest disturbance areas occurred in 2001 and 2002, with similar areas of about 143,873 ha and similar disturbance rates of about 1.72%. The smallest disturbance area was observed in 1988, with an area of 31,708.98 ha and a disturbance rate of 0.38%. The largest recovery area was 293,112.5 ha in 1994, with a recovery rate of 3.5%, followed by 1997, with a recovery area of 182,858.5 ha and a recovery rate of 2.19%. The smallest recovery area was observed in 2004, with a recovery rate of 0.46%. Overall, except for these cases, the majority of the forest disturbance rates and recovery rates were less than 2%, and they fluctuated over time.  Figure 4 shows the spatial distribution of the burned areas extracted from the VCT-SVM algorithms, and Table 4 shows the error matrix derived from the independent samples. The classification accuracy of the fire patches was high, with an overall accuracy of 94% or higher for the four years. The user's and producer's accuracies were also above 90%, except for the 2016 producer's accuracy (88.9%) of the non-burned class. The results show that Figure 4 is accurate and reliable, and the map can be used for the subsequent vegetation recovery analysis. Only two classes (burned area and non-burned area) were considered in the SVM classification. The burned areas were flagged as disturbance events by VCT, and small burnt areas (<1 ha) were excluded from the analysis. These strategies guaranteed the high extraction accuracy of the fire patches.

Model Fitting and Validation
We first used the variable importance analysis in the RF package to identify the 25 most important potential environmental variables ( Figure 5). Next, we implemented a stepwise regression analysis of the 25 potential variables to address the multicollinearity between the 25 variables to obtain the final set of 11 variables (Table 5).

Figure 5.
Variable importance ranking of all predictor variables. The acronyms of the predictive variable consist of two parts: abbreviation and number. Number means from the year after fire (the year of fire is 0) in the burned area. Abbreviation description: mean_tem (average annual temperature), ex_high (Extreme maximum temperature), ex_low (extreme minimum temperature), mean_high (average maximum temperature), mean_low (average lowest temperature), moisture (average relative humidity), preci (annual precipitation). A detailed description can be found in Table 2.). The modeling results indicated that climatic factors played an important role in postdisturbance vegetation recovery of the first few years. The relative importance of the predictor variables ( Figure 5) showed that the meteorological factors exhibited a good performance for vegetation recovery prediction. The precipitation factors in our study had high variable importance, indicating that precipitation was one of the limiting factors for vegetation recovery in this area. Although there was a good correlation between NDVI_5 and moisture_0, other moisture factors exhibited a poor performance for predicting the vegetation recovery status in this region. Extreme high temperatures, especially ex_high_3, performed well for vegetation recovery prediction. Figure 6 shows the model fitting and validation results. RF outperformed SVM and SMLR regarding modeling performance. The R 2 of the RF model fitting was 0.926, and the MAE and RMSE were lower than those of the SVM. The R 2 , MAE, and RMSE of the SVM fitting were 0.8288, 0.025, and 0.044, respectively. The results further showed that the validation R 2 of the RF was higher than that of the SVM model (0.7113), with the highest validation R 2 values of 0.7727. The model evaluation results showed that the modeling accuracy and validation accuracy were higher for the RF than the SVM and SMLR ( Figure 6). Overestimation of low values and underestimation of high values occurred in all three models, but these effects were least pronounced in the RF model.

Fire Disturbance and Recovery
Since fire is a crucial factor resulting in forest decline in the Daxing'anling area [54], we assessed fire disturbance and recovery in this region using time-series analysis of remotely sensed spectral indices. The results were consistent with Zhang's research describing a known burned forest area in the Daxing'anling area [55]. Our study identified a large disturbance event in 1987, which was consistent with the "May 6th" catastrophic fire in 1987. This was a catastrophic fire, leading to extensive loss of life and property and one million hectares of burned forest area (both inside and outside the study area) [56]. Postfire recovery is affected by fire severity, biogeographic conditions, climate factors, forest management strategy, and other factors [30]. Thus, there are various recovery methods for different fires, including no land-use change, tree planting, or natural restoration; these practices have different effects on vegetation restoration [56]. Fire and post-fire recovery are crucial processes in sustainable forest management that shape the forest age, forest structure and forest type in northeastern China [30,57].

Post-Fire Recovery Response to Climate
The climatic factors had a substantial influence on post-fire vegetation recovery in this region, especially temperature factors play an important role in vegetation recovery in this region. Minore et al. pointed out that temperature was a critical factor for forest growth when sufficient moisture is available [58]. Extremely high temperatures (e.g., ex_high_3) and extremely low temperatures (e.g., ex_low_2) may limit post-fire vegetation recovery. Extreme temperatures might be a proxy for drought stress or indicate high solar radiation and temperature [30,59]. Precipitation factors were also important in this research, especially in the early stage (e.g., preci_1 and preci_2), which was in line with the findings by Meng et al. [30]. Liu also pointed out that higher growing-season precipitation one year after a fire had a positive influence on post-fire recovery [60]. Additionally, we found that topographic variables (especially elevation) in the modeling was important, which was supported by Bright [61]. However, this study selected the 5-year post-fire NDVI value as the post-fire state of recovery vegetation [22], and there may be a relationship between the elevation and the biomass amount [50], elevation was possibly less influential than our modeling results suggested. Lacking consideration of pre-disturbance state [22] can also explain the low importance of the dNBR, which was opposite to previous research [30]. Since the dNBR represents the difference in the NBR before and after the fire, the correlation between the dNBR and NDVI_5 (the amount of greenness at that time) was not high. However, fire creates favorable conditions for fire-adapted forest types to germinate and regenerate; therefore, some areas burned at higher severities and recovered more quickly than the areas at lower severities [31].

Vegetation Recovery Modeling
Although the overall performances of the three models are acceptable, the validation performances are inferior to the fitting performances, and overestimation for low values and underestimation for high values occurred in all three models ( Figure 5). The reason may be that the fitting of the three models requires a large number of samples, and there were far more samples with intermediate values than with low values and high values, leading to inadequate or biased training of the model. In addition, the machine learning models, such as RF and SVM, performed better than SMLR. The likely reason is that most variables were non-normally distributed, causing visibly nonrandom trends in the residuals of the initial linear models [61].

Limitations and Future Improvements
Although important results were obtained in this study, the following aspects require further research: (1) Due to the limitations of VCT monitoring, ground-truth data, and the fire point data, some low-severity or small fires were not detected in the time-series data. Therefore, more sensitive spectral indices and more fixed-point field surveys are required in future research. (2) Vegetation indices are typically used for burned area extraction. However, some vegetation indices are highly correlated (such as NDVI and EVI), which may adversely affect the extraction result. (3) As mentioned above, we need to consider the pre-fire conditions to describe vegetation recovery. João used three indicators (compared with pre-fire condition) to measure different facets of the post-fire recovery process, which might be a better approach to capture recovery than NDVI_5 in our study [22]. (4) We used a small number of samples (3 quadrats) to verify the VCT disturbance results, which may have limited the ability to capture the heterogeneity and diversity of spatial, environmental, and spectral conditions in the area. Stratified random sampling will be considered to improve this problem [25]. (5) We did not consider the differences in forest type or tree species and only evaluated post-fire vegetation recovery of all forested areas. Different tree species may exhibit different recovery patterns after a fire and different responses to climate [30]. In future studies, different tree species and forest types will be considered to compare the post-fire recovery mechanism and provide comprehensive forest management suggestions after a fire.

Conclusions
Landsat time-series analysis can describe post-fire vegetation recovery at the landscape scale. In this study, a VCT-SVM method was used to extract fire disturbance areas and investigate vegetation recovery in burned areas. The RF, SVM, and SMLR models explained the relationship between post-fire vegetation recovery and different predictor variables, and the RF algorithm provided the best performance for predicting vegetation recovery after fire events. A good coupling relationship was observed between the climatic factors and post-fire vegetation recovery. The results of this study are applicable to other fire areas, where landscape-scale information on post-disturbance vegetation recovery is desired. The findings can inform management decisions and are suitable for carbon forecasting in forest ecosystems.

Data Availability Statement:
Publicly available datasets were analyzed in this study. This data can be found here: [https://espa.cr.usgs.gov/, http://data.cma.cn/ (accessed on 24 December 2020)]. Fire point data was obtained from the Shenyang Institute of Applied Ecology at the Chinese Academy of Sciences and are available with the permission of the Shenyang Institute of Applied Ecology at the Chinese Academy of Sciences.

Acknowledgments:
The authors would like to acknowledge the United States Geological Survey (USGS) and National Aeronautics and Space Administration (NASA) for providing the Landsat image data. Special thanks to the Shenyang Institute of Applied Ecology at the Chinese Academy of Sciences and the China Meteorological Data Network for sharing their fire data and climate data, respectively.

Conflicts of Interest:
The authors declare no conflict of interest.