Spatiotemporal Assessment of Forest Fire Vulnerability in China Using Automated Machine Learning

Abstract

Frequent forest fires cause air pollution, threaten biodiversity and spoil forest ecosystems. Forest fire vulnerability assessment is a potential way to improve the ability of forests to resist climate disasters and help formulate appropriate forest management countermeasures. Here, we developed an automated hybrid machine learning algorithm by selecting the optimal model from 24 models to map potential forest fire vulnerability over China during the period 2001–2020. The results showed forest aboveground biomass (AGB) had a vulnerability of 26%, indicating that approximately 2.32 Gt C/year of forest AGB could be affected by fire disturbances. The spatiotemporal patterns of forest fire vulnerability were dominated by both forest characteristics and climate conditions. Hotspot regions for vulnerability were mainly located in arid areas in western China, mountainous areas in southwestern China, and edges of vegetation zones. The overall forest fire vulnerability across China was insignificant. The forest fire vulnerability of boreal and temperate coniferous forests and mixed forests showed obviously decreasing trends, and cultivated forests showed an increasing trend. The results of this study are expected to provide important support for the forest ecosystem management in China.

1. Introduction

Forests are key components in terrestrial systems that provide many ecosystem services contributing to climate mitigation [1,2] and human well-being [3,4]. Forest fire remains one of the main natural and widespread disturbance factors. Although forests can resist some surrounding environmental changes, fires damage forest ecosystems and increase their vulnerability to environmental changes [5,6,7]. As a result, quantifying forest vulnerability to fire disturbances and understanding the underlying mechanisms are critical for projecting climate impacts and developing effective adaptation strategies [8].

Forest ecosystem vulnerability results from both natural and anthropogenic activities [9,10], and high vulnerability generally results in reduced forest area, dramatic mass reduction in biomass, gradual loss of biodiversity, decline of forest ecological functions and deterioration of system stability [11,12]. However, a widely accepted and operational framework for evaluating forest ecosystem vulnerability to disturbances, including exotic pests, habitat destruction, and fire [13,14,15], remains to be outlined [16]. Forest fire has been viewed as the most common and serious hazard in forest ecosystems [17] because of the warmer and drier climate conditions [18,19,20]. Forest fire vulnerability is often defined as the potential loss from fire [21]. Several scholars have investigated forest fire vulnerability at regional scale. Amalina et al. [22] developed a fire vulnerability map in Way Kambas National Park based on biophysical and human disturbance factors and classified the vulnerability into high, medium, and low levels. Dupire et al. (2019) explored alpine forest fire vulnerability by simulating fire behavior and ultimately post-fire tree mortality using surface fuel measurements, forest inventory and fire weather data. They found that summer fires occurring in extremely dry conditions promote high mortality across all forest communities with different climatic conditions. Buotte et al. [23] calculated and predicted recent forest vulnerability to drought and fire in the western United States based on simulated future changes in burned area relative to historical burned area. They indicated that the proposed drought vulnerability metrics can be incorporated into Earth system models as probabilistic mortality rates. Forzieri et al. [8] defined forest vulnerability as the relative biomass loss following the occurrence of a given disturbance and quantified the vulnerability of European forests to fires, windthrows and insect outbreaks during the period 1979–2018 by the random forest regression method. They showed that about 33.4 billion tonnes of forest biomass could be seriously affected by these disturbances, with higher relative losses when exposed to windthrows (40%) and fires (34%) compared with insect outbreaks (26%). Overall, the current methods for vulnerability estimation can be summarized into two categories: knowledge-based and data-based methods. The knowledge-based method is determined by expert systems and a priori knowledge to capture and assess vulnerability, for example, the AHP-fuzzy integrated evaluation method [24], landscape ecology method [25], gray cluster analysis method [26], and principal component analysis method [27]. Such methods typically adopt a priori knowledge to identify the functional relationships and the possible amplification or dampening effects that may emerge at local scale from interactions among multiple factors [8]. The data-based method applies in vulnerability assessment with more complete data, mainly including machine learning or deep learning methods [28]. With the advancement of vulnerability concepts, the enrichment of indicator systems, and the expansion of data availability in recent years, popular methods have shifted from mathematical knowledge-based methods to spatiotemporal analysis in view of machine learning models. Machine learning can effectively identify patterns and trends of complex problems without human intervention. The newly developed automated machine learning (AutoML) framework is more efficient and flexible than a single machine learning algorithm. It has the advantages of accurately selecting specific algorithms for hyperparameters, using cross validation to adjust the model to avoid over fitting, and training stacked models [29,30]. As a comprehensive issue, forest fire vulnerability is involved in forest vegetation characteristics, climate and landscape characteristics. Therefore, it is necessary to assess forest fire vulnerability using AutoML at large scales.

For a long time, forest fire prevention and control have been very important factors included in China’s ecological planning. However, due to natural and historical factors, forest fires have occurred frequently in China. According to the China Forestry Statistical Yearbook, a total of 58,545 forest fire accidents occurred in China from 2008 to 2018 [31]. Although they displayed an overall downward trend, the frequency of forest fires is still high [32]. Currently, forest fire vulnerability in China is primarily estimated using data from provincial field surveys with a coarse resolution (around 300,000 km²) [33,34,35,36], and detailed forest fire vulnerability under different vegetation types and climatic zones is still lacking. The increasing high-spatial-resolution satellite datasets could support large-scale comprehensive efforts for regional forest parameters and disturbances. Compared with traditional forest surveys, remote sensing has the advantages of large detection range, fast data acquisition, fewer ground restrictions, and a large amount of information acquisition [37,38]. Technologically, the combination of machine learning and massive computation assists large-scale ecosystem vulnerability assessment effectively [39,40]. Thus, more systematic and in-depth research on forest vulnerability assessment in China is necessary [41,42].

In this paper, forest fire vulnerability was investigated using AutoML over China from 2001 to 2020. The specific objectives of this study were to: (1) develop a novel framework for vulnerability assessment using multisource datasets and AutoML; (2) explore the spatial and temporal patterns of vulnerability and detect possible forest hotspots that are highly sensitive to fire disturbance; and (3) discover the driving factors affecting forest vulnerability to fire. This work will shed fresh light on the quantification of the spatial and temporal vulnerability of natural disturbances on ecosystems.

2. Study Area

Forests cover approximately 21.63% (2009–2013) of the land surface in China and have a carbon stock of approximately 84.27 billion tons [43]. China’s forests range from boreal forests in the northeast to tropical rainforests in south and southwest, which include almost all major forest biomes in Northern Hemisphere [44]. They provide a set of ecosystem services that contribute to human well-being and climate regulation.

Forest distribution and types (Figure 1) in China were obtained by integrating the forest cover data from the Terra MODIS Vegetation Continuous Fields (VCF) product [45] with the 1:1,000,000 vegetation type map from The Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (RESDC, http://www.resdc.cn, accessed on 6 November 2021). We extracted the areas with tree cover greater than 10% and patches area greater than 0.5 hectares as forests from annual VCF product from 2001 to 2020 [46]. Then, the forests were further classified into alpine vegetation, broad-leaved forest, mixed forest, coniferous forest, shrub forest and artificial forest based on the 1:1,000,000 vegetation type map. To investigate the spatial differences of fire vulnerability, we used the Vegetation Regionalization Map of China from RESDC. From east to west, vegetation distribution was determined by moisture conditions, ranging from humid to dry, the area was divided into forest zones (I, II, III, IV and V), temperate grasslands zone (VI) and temperate desert zone (VII) [47]. From north to south, vegetation distribution was determined by the temperature, forest zones were further categorized into cold-temperate needleleaf forest zone (I), temperate coniferous and broadleaf mixed forests zone (II), warm-temperate deciduous broadleaf forests zone (III), subtropical evergreen broadleaf forests zone (IV) and tropical monsoon rainforest and rainforest zone (V) [48]. Considering the influence of vertical zonality or other non-zonal factors (such as geomorphic structure), the Qinghai-Tibetan Plateau alpine vegetation zone was divided.

3. Data and Methodology

Forest fire vulnerability was expressed as the relative loss of aboveground biomass (AGB) following the occurrence of a given fire disturbance. The workflow of the annual potential forest fire vulnerability map in China is presented in Figure 2. In this section, we provide a list of data collection and preprocessing (Section 3.1), dimensionality reduction in explanatory variables (Section 3.2.1), vulnerability models evaluation and selection (Section 3.2.2 and Section 3.2.3), and spatiotemporal patterns and key drivers of forest fire vulnerability (Section 3.2.4).

3.1. Data Collection and Preprocessing

Fire data were retrieved from the Global Wildfire Information System (GWIS, https://gwis.jrc.ec.europa.eu/ (accessed on 9 October 2022)) using the MODIS dataset MCD64A1, and this product was computed by an algorithm relying on encoding in a graph structure a space–time relationship among patches of burned areas [49]. Each fire record was represented by a vector feature describing the spatial range of the damaged forest patch. The results were verified by visual aerial photointerpretation, satellite imagery, and field surveys [49]. We obtained a fire record database containing 107,168 forest fire records in China from 2001 to 2019 by overlying the observed fire area with the annual forest area.

Time series annual AGB are required to assess the expected AGB loss for a given fire disturbance event. To reconstruct annual AGB changes over 2001 and 2020, we integrated 100 m aboveground AGB map (AGB₂₀₁₀) in 2010 as the reference data from the European Space Agency-Data User Element Global Biomass Project (https://globbiomass.org/products/globalmapping/ (accessed on 15 September 2021)) [50] and MODIS VCF product [45]. Time series AGB maps were reproduced under the assumption that AGB changes are entirely dependent on the change in forest cover, similar to the approach used by Forzieri et al. [8]. For each fire-disturbed forest area at year t in the fire record database, the relative AGB loss (AGBL_rel) was quantified based on the difference of the AGB before and after disturbance. We also computed one-year-lagged AGB (AGB_lag) estimates to retrieve information on carbon stocks of the recent past. The equations are as follows:

$$AGB{L}_{rel}=\frac{\mathit{max}\left(AG{B}_{t-1},AG{B}_t\right)-min\left(AG{B}_t,AG{B}_{t+1}\right)}{\mathit{max}\left(AG{B}_{t-1},AG{B}_t\right)}$$

(1)

where $AGB{L}_{rel}$ is the relative AGB loss caused by fire disturbance. $AG{B}_t$ is the AGB map of fire disturbances in year t.$AG{B}_{t-1}$ is the AGB map before fire disturbances for year t − 1. $AG{B}_{t+1}$ is the AGB map after fire disturbance for year t + 1. To focus on effective fire disturbance events, records with $AGB{L}_{rel}$ > 5% were selected, and yielding 90,208 events in total.

We collected 26 variables, including forest features, climate and landscape factors, which were used as potential predictors in the vulnerability models.

Table 1. Explanatory variables used in vulnerability models.

Category	Explanatory Variables	Abbreviation	Temporal Resolution
Forest features	Aboveground biomass	AGB	Yearly
	one-year-lagged aboveground biomass	AGBlag	Yearly
	Leaf area index	LAI	Yearly
	Fraction of absorbed photosynthetically active radiation	FAPAR	Yearly
	Specific leaf area	SLA	Static
	Leaf dry matter content	LDMC	Static
	Leaf nitrogen content per leaf dry mass	LNC	Static
	Leaf phosphorus content per leaf dry mass	LPC	Static
Climate features	Annual cumulated precipitation	Pcum	Yearly
	Annual average temperature	Tavg	Yearly
	Annual maximum temperature	Tmax	Yearly
	Annual values of downward surface Shortwave radiation	SRad	Yearly
	Palmer drought severity index	PDSI	Yearly
	Soil moisture	SM	Yearly
	Long-term mean precipitation	Pcum_long	Static
	Short-term average anomaly in cumulated precipitation	avg aPcum	Yearly
	Long-term mean temperature	Tavg_long	Static
	Short-term average anomaly in average temperature	avg aTavg	Yearly
	Moisture index	MI	Yearly
Landscape features	Coefficient of variation	CV	Static
	Evenness index	EI	Static
	Shannon Index	SI	Static
	Homogeneity index	HI	Static
	Elevation	Ele	Static
	Population density	PD	Yearly

shows the explanatory variables used in the vulnerability models. These variables were collected from multiple sources, including satellite images and reanalysis products. Data collection and processing were completed on the Google Earth Engine cloud platform. Please see Appendix A for specific information on explanatory variable extraction and processing.

We extracted the explanatory variables for each of the fire records to train the vulnerability models. Specifically, the explanatory variables were spatially averaged according to each fire-burned forest area. Considering the loss of forest features before and after fire, we deleted the fire records if the variables had null values or the change in LAI was less than 5%. Finally, we built a forest fire record database with 41,188 fire records. Each fire record consisted of multiple attributes, including AGBLrel and multiple explanatory variables. Spatial distribution maps of all explanatory variables were obtained and converted to raster with 1000 m resolution for the fire vulnerability prediction in China.

3.2. Methodology

3.2.1. Vulnerability Indicator Construction

The variance inflation factor (VIF) was widely used to quantify the collinearity among the explanatory variables [51,52,53]. VIF is calculated from the correlation determination coefficient in the regression analysis:

$$VIF=\frac{1}{1-R_i^2}$$

(2)

where $R_i^2$ is the square of the correlation determination coefficient of independent variable xi for the regression analysis of other independent variables. The VIF values were determined using the “car” package in R version 4.1.1. The results showed that the VIF values of FPAR, T_avg and T_{avg_long} were greater than 10 (Appendix A: Figure A1). Thus, we removed these 3 variables and finally obtained 23 explanatory variables for vulnerability modeling.

3.2.2. Vulnerability Modeling Based on the Regression Method of Auto ML

Vulnerability expresses the degree to which a forest ecosystem is affected when exposed to a given fire disturbance. In this study, we quantified the vulnerability in terms of potential relative AGB losses (PBL_rel) following the occurrence of a specific fire (0 means a forest is not vulnerable to the given fire events, and 1 means a forest is completely damaged when exposed to the given fire events).

We developed and tested an optimal AutoML model to predict the observed AGBL_rel (response variable) based on 23 explanatory variables (predictors/features). We randomly selected 75% (30,888 fire records) of the forest fire record database for training the vulnerability model and 25% (10,300 fire records) for verifying the results.

AutoML makes the machine learning available for non-machine learning experts, improves efficiency of machine learning and accelerates research on machine learning [54,55]. H₂O AutoML interface was used to generate a regression model for fire disturbance vulnerability estimation. H₂O AutoML is an machine learning framework that enables the automated training of a large selection of candidate models on a dataset [39]. The framework included many popular machine learning algorithms, such as the simple/deep neural network (NN) [56,57], generalized linear model (GLM) [58,59], gradient boosting machine (GBM) [60,61], distributed random forest (DRF) [62,63], and stacked ensemble learning of these models [54,64]. The H₂O AutoML interface performs random grid searches to select machine learning models and automatically completes hyperparameter optimization for each model. Finally, 24 models were developed based on the training set.

3.2.3. Performance Metrics

To evaluate the performances of the AutoML models, the estimated PBL_rel was evaluated against the AGB_rel by the root mean square error (RMSE), mean square error (MSE), mean absolute error (MAE) and root mean squared logarithmic error (RMSLE). The coefficient of determination (R²) between the simulation results of the optimal model and the AGB_rel in the test set was also calculated. These statistical indexes are calculated using the following equations:

$$RMSE=\sqrt{\frac{1}{N}{{\displaystyle \sum }}_{i=1}^N{\left({\theta }_i^{Est.}-{\theta }_i^{Obs.}\right)}^2}$$

(3)

$$MSE=\frac{1}{N}{{\displaystyle \sum }}_{i=1}^N{\left({\theta }_i^{Est.}-{\theta }_i^{Obs.}\right)}^2$$

(4)

$$MAE=\frac{1}{N}{{\displaystyle \sum }}_{i=1}^N|{\theta }_i^{Est.}-{\theta }_i^{Obs.}|$$

(5)

$$RMSLE=\sqrt{\frac{1}{N}{{\displaystyle \sum }}_{i=1}^N{\left(\mathit{log}\left({\theta }_i^{Obs.}+1\right)-\mathit{log}\left({\theta }_i^{Est.}+1\right)\right)}^2}$$

(6)

$$R^2=1-\frac{{\sum }_{i=1}^n{({\theta }_i^{Est.}-{\theta }_i^{Obs.})}^2}{{\sum }_{i=1}^n{({\theta }_i^{Est.}-\overline{{\theta }^{Obs.}})}^2}$$

(7)

where ${\theta }_i^{Est.}$ is the i-th estimated PBL_rel and ${\theta }_i^{Obs.}$ is the i-th sampled AGB_rel in the test set. $\overline{{\theta }^{Obs.}}$ is the average value of all sampled AGBrel in the test set, and n is the number of sampled AGB_rel in the test set.

3.2.4. Spatiotemporal Patterns and Key Drivers of Forest Fire Vulnerability

The best-performing model in AutoML was used to evaluate the vulnerability of forests annually between 2001 and 2020 for each grid cell (1 × 1 km²) in China. The potential absolute AGB loss (PBL_obs) was constructed by multiplying the estimated PBL_rel by the available biomass. Therefore, the vulnerability value in each grid cell reflects the biomass (relative or absolute) that will be affected when exposed to fire disturbance based on its specific local and temporal environmental conditions. The temporal and spatial vulnerability patterns and driving mechanisms were analyzed in R.

To better understand fire disturbance vulnerability in various forest biomes, we compared the vulnerability in different vegetation zones by regional statistical methods using the “raster” and “rgdal” packages in R.

Long-term trends in vulnerability were quantified for 2001–2020 by forest types and grid level, and their trend and significance were measured by 2-sided Mann–Kendall tests using the “Kendall” package in R. Kendall’s rank correlation coefficient (τ) indicated the strength of the vulnerability trend, and positive or negative values indicated an upward or downward trend, respectively. Kendall’s p represents the significance level of the trend statistics, with the significance set to p < 0.1.

The variable importance tool describes the extent to which covariables contribute to the results of the prediction model [65]. We quantified and ranked how individual explanatory variables affected forest fire vulnerability through “caret” package and “iml” package in R. For each of the explanatory variables, we also performed a correlation analysis comparing vulnerability against each of the datasets at the pixel level using Spearman’s correlation coefficient by gridcorts in R [66].

4. Results

4.1. Vulnerability Model Evaluation

Table 2. The statistical error results of the vulnerability models of PBLrel and AGBrel in the training set (root mean squared error (RMSE), mean square error (MSE), mean absolute error (MAE) and root mean squared logarithmic error (RMSLE)).

Number	Model ID	RMSE	MSE	MAE	RMSLE
1	StackedEnsemble_AllModels	0.095	0.009	0.067	0.074
2	StackedEnsemble_BestOfFamily	0.097	0.009	0.069	0.075
3	GBM_4_AutoML	0.098	0.010	0.069	0.076
4	GBM_3_AutoML	0.098	0.010	0.070	0.076
5	GBM_grid__1_AutoML_model_1	0.099	0.010	0.070	0.077
6	GBM_2_AutoML	0.099	0.010	0.070	0.077
7	GBM_1_AutoML	0.100	0.010	0.071	0.077
8	GBM_grid__1_AutoML_model_4	0.100	0.010	0.071	0.077
9	GBM_5	0.101	0.010	0.072	0.078
10	GBM_grid__1_AutoML_model_8	0.102	0.010	0.075	0.079
11	GBM_grid__1_AutoML_model_10	0.104	0.011	0.076	0.080
12	GBM_grid__1_AutoML_model_7	0.104	0.011	0.076	0.080
13	GBM_grid__1_AutoML_model_3	0.104	0.011	0.076	0.080
14	GBM_grid__1_AutoML_model_6	0.104	0.011	0.076	0.081
15	GBM_grid__1_AutoML_model_11	0.105	0.011	0.078	0.081
16	DRF_1	0.106	0.011	0.080	0.082
17	XRT_1	0.106	0.011	0.080	0.082
18	GBM_grid__1_AutoML_model_9	0.108	0.012	0.082	0.084
19	GBM_grid__1_AutoML_model_5	0.108	0.012	0.082	0.084
20	GLM_1	0.139	0.019	0.110	0.107
21	DeepLearning_1_AutoML	0.149	0.022	0.108	0.116
22	DeepLearning_grid__1	0.153	0.024	0.118	0.119
23	DeepLearning_grid__2	0.175	0.031	0.128	0.133
24	DeepLearning_grid__3	0.201	0.040	0.165	0.154

shows the statistical error results of the vulnerability models of PBL_rel and AGB_rel in the training set. The accuracy of different types of models ranked as follows: stacked ensemble > GBM > DRT > XRT > GLM > deep learning, according to the model accuracy from high to low. Through comparative analysis, we found that the accuracy evaluation results of the above machine learning models are similar. The optimal model explained 0.095 of the RMSE, with values of 0.009, 0.009, 0.067 and 0.074 for the mean residual deviance, MSE, MAE and RMSLE, respectively. The RMSE, MSE, MAE and RMSLE of the last ranked model were 0.201, 0.040, 0.165 and 0.154, respectively. Figure 3 shows R² between the optimal model simulation PBLrel and AGBrel in the test set. The optimal model explained 63% of the correlation between the estimated relative biomass loss and the samples. This model was used to further generate the vulnerability distribution over China.

4.2. Spatiotemporal Patterns of Vulnerability

4.2.1. Spatial Distribution of Vulnerability

Figure 4 shows the average potential vulnerability of Chinese forests to fire disturbances from 2001 to 2020. The vulnerability to fires exhibited a potential relative AGB loss (PBL_rel) of 26% (25.9–26%, 95% confidence interval). Considerable spatial variations in vulnerability emerged across China. Vulnerability to fire outbreaks generally occurred in dry and sparsely forested regions, with the PBL_rel increasing from east to west and reaching as much as 80% (Figure 4a). Vulnerability to fires appeared to be high in the temperate grassland zone (VI), tropical monsoon rainforest and rainforest (V), temperate desert zone (VII) and subtropical evergreen broadleaf forest zone (IV), especially in the mountainous regions (such as the Kunlun and Tianshan Mountains, Qilian and Helan Mountains, Himalayas and Hengduan Mountains), with the median value of PBLrel exceeding 25% and reaching 75% locally. However, forests in the boreal and temperate coniferous forest zone (I) and temperate coniferous and broadleaf mixed forest zone (II) with cold and wet climates generally had lower values (Figure 4b).

Figure 5 shows PBL_rel, PBL_abs and the regional average AGB in the different vegetation zones. The lowest forest vulnerability to fire occurred in zones I and II, with values of 18.66% and 19.11%, respectively, and the mean AGB of these two regions was relatively high. However, they may not have been sufficiently disaster-prone environments (i.e., the annual rainfall in these areas was usually more than 600 mm, and the accumulated temperature was also low). The PBL_abs values in these regions were not the highest, with values of approximately 0.14 and 0.21 Gt C/year. The vulnerability of VII and VIII was the highest, with values of 33.36% and 33.01%, respectively. Although the AGB in these places was low, because of the harsh climate conditions, their vulnerability was the highest, and the PBL_abs was the lowest. The vulnerability of IV was also relatively high, which may be because the region provided sufficient fuel (the regional average AGB was the highest, at 39.79 mg C/hm²) and suitable climatic conditions for fire. In summary, climate and combustible conditions may be the main factors affecting fire vulnerability. As climate gradients such as precipitation and temperature become harsher, vulnerability to fire disturbance will increase.

4.2.2. Spatiotemporal Trend of Vulnerability

Figure 6 indicated the spatial distribution of the forest vulnerability trend from 2001 to 2020. In all, 89% of vegetated grids in China exhibited a nonsignificant trend. A significant decrease in forest vulnerability occurred in zone I and zone II (overall 7% of the forest vegetated grids), and a significant increase (4% of the forest vegetated grids) was mainly distributed in vegetation zone IV and zone VII.

The temporal trends of forest vulnerability (PBL_rel) to fire with different forest types from 2001 to 2020 are shown in Figure 7. Generally, alpine vegetation showed the highest mean vulnerability, followed by cultivated forests and shrub forests; in contrast, natural forests, especially mixed forests, had the lowest mean vulnerability. No significant temporal trend of forest fire vulnerability was found in terms of the mean vulnerability across China (τ = 0.126, p value = 0.656). In detail, only the cultivated forests increased significantly, with τ = 0.316 and p value = 0.056; mixed forests (τ = −0.379, p value = 0.021) and broad-leaved forests (τ = −0.284, p value = 0.086) showed obvious reductions in vulnerability to fire.

4.3. Key Drivers of Forest Fire Vulnerability

We quantified and ranked how individual explanatory variables affected forest fire vulnerability (Figure 8). The results showed that the forest features (e.g., AGB, AGB_lag) influenced forest vulnerability the most, followed by the climate features and landscape features, with mean importance scores of 60.06, 4.83 and 3.20, respectively. Plant water stress indicated by low precipitation (P_cum) and the PDSI, and temperature stress indicated by the SRad, avg T_avg, and T_max associated with very large fuel loads acted to further increase vulnerability. The HI, SI, EI and CV had the lowest importance scores, which was possibly due to the lower resistance and susceptibility to fire spread.

Correlations between several important features and forest vulnerability are shown in Figure 9. Generally, the AGB and LAI were negatively correlated with vulnerability, with mean Spearman’s coefficients of −0.24 and −0.07, respectively. In particular, a negative correlation occurred in regions with a high coverage of coniferous forest (e.g., Greater Khingan, Changbai Mountains). The PDSI (Figure 9h), P_cum (Figure 9c) and avg aP_cum (Figure 9f) were negatively correlated with forest vulnerability, with mean correlations of −0.04, −0.03 and −0.01, respectively. The temperature factors T_avg (Figure 9d), T_max (Figure 9e) avg aT_avg (Figure 9g), and SRad (Figure 9i) with mean correlations of 0.02, 0.04, 0.05 and 0.05, respectively, were positively correlated with forest vulnerability.

5. Discussion

In this study, we investigated forest vulnerability to fire in China during 2001–2020, and multisource satellite products and explanatory variables were integrated into AutoML to acquire the most optimized model. Forest vulnerability to fire varied across China, and hotspot regions were concentrated in arid areas in western China, mountainous areas in southwestern China, and edges of vegetation zones. These results were consistent with the high fire frequency reported in previous studies in Tibet [67], northeast China’s Greater Khingan Range and the subtropical dry forest in southern Yunnan [68]. Concerning the forest vulnerability trend, few regions expressed significant reductions, especially in zones I and II of China, where the Three-North Shelterbelt [69,70] has been implemented to protect natural temperate forests. Additionally, decreased vulnerability indicated that forest susceptibility to fire improved due to strict forest protection policies, such as the complete prohibition of commercial logging of natural forests [71,72]. The forest area of zones I and II accounted for 6% of China’s forest area, but the regions with significantly reduced vulnerability accounted for 43% of China’s region.

Integrating all of the information regarding forest vulnerability, as we did in this study, has rarely been previously proposed at the national scale. The solution was a model with global representativeness, resistance to forest variability and environmental heterogeneity, and strengthened applicability and reliability. The bottom-up evaluation method is broadly used to estimate fire emissions [73,74], which are caused by burned area, fuel load, combustion and emission factors [67,75]. In this study, the potential affected AGB (PBL_abs) was calculated by multiplying biomass by the carbon loss ratio (AGBL_rel), while the fuel load, combustion and emissions were not determined. Therefore, the overestimation of AGB losses existed, which could lead to uncertainties in forest vulnerability evaluation. Quantifying the uncertainties from input data, model configuration, and other sources to optimize the model is necessary in future research.

Ecological projects, such as afforestation, returning farmland to forests and natural forest conservation, have provided substantial and long-lasting benefits in boosting the vegetation carbon sink of China’s terrestrial ecosystem [76,77]. Forest areas in China increased by 1.61 × 104 km² from 2000 to 2020 [78]. The forest vulnerability results among various forest types demonstrated that the forest vulnerability of plantations had a significant increasing trend, while the fire vulnerability of mixed forest and broad-leaved forest had a significant decreasing trend (Figure 6). Due to the differences in forest structure (e.g., forest age, forest size distribution, and canopy cover), plantations are more vulnerable to disasters [79,80]. Overall, our strategy revealed that forest vulnerability varied with forest structural characteristics and environmental conditions. In view of the expected persistence of future climate warming and the possible intensification of natural disturbances [14,81], forest ecosystem services may be challenged by these disturbances in the future. Our strategy can provide a potential way to understand the current vulnerability of forest carbon stocks and improve the ability of forests to resist disasters.

6. Conclusions

In this study, we proposed a new framework using AutoML to estimate forest vulnerability in China from 2001 to 2020. Multiple data were used to train the most optimized model and evaluate the estimations. The results showed that forest AGB had a vulnerability of 26%, corresponding to approximately 2.32 Gt C/year of forest AGB that could be affected by fire disturbances. The overall accuracy of the best vulnerability AutoML model expressed in terms of the R² was approximately 63% of the correlation across the considered disturbances. The spatial pattern of forest vulnerability was strongly controlled by the interaction between forest features and climate. Hotspot regions for forest vulnerability were mainly located in the arid areas in western China, the mountainous areas in southwestern China, and at the edges of vegetation zones. No clear trend in overall forest vulnerability was found, except in boreal and temperate coniferous forest zones at high latitudes. Furthermore, the abnormal increase in temperature in the short term and water shortage significantly increased forest vulnerability to fire. Our study sheds light on the vulnerability of China’s forests to fire disturbances and the ongoing trends in response to changing climate conditions. The proposed methodology can provide new tools and insights for forest management.