Remote Sensing of Burn Severity Using Coupled Radiative Transfer Model: A Case Study on Chinese Qinyuan Pine Fires

: Burn severity mapping is critical to quantifying ﬁre impact on key ecological processes and post-ﬁre forest management. Satellite remote sensing has the advantages of high spatial-temporal resolution and large-scale monitoring and provides a more e ﬃ cient way to evaluate forest ﬁre burn severity than traditional ﬁeld or aerial surveys. However, the proportion of tree canopy cover (TCC) a ﬀ ects the spectral signal received by remote sensing sensors from the background charcoal and ash. Consequently, not considering this factor normally leads a spectral confusion in burn severity retrieval. In this study, the burn severity of two Qinyuan forest ﬁres was estimated using a coupled Radiative Transfer Model (RTM) and Sentinel-2A Multi-Spectral Instrument (MSI) reﬂectance data. A two-layer Canopy Reﬂectance Model (ACRM) RTM was coupled with the GeoSail RTM by replacing the spectra of the background input of GeoSail RTM to simulate the spectra of the three-layered forests for burn severity retrieval measured as the Composite Burn Index (CBI). The TCC data was then served to RTM parameterization and constrain the backward inversion procedure of the coupled RTM to alleviate spectral confusion. Finally, the inversion retrievals were evaluated using 163 ﬁeld measured CBI. The coupled RTM can simulate the radiative transfer characteristics of three-layer vegetation and has greater potential to accurately estimate burn severity worldwide. To evaluate the merit of our proposed method, the CBI was estimated through coupled RTM inversion with TCC constraint (CP_RTM + TCC), coupled RTM inversion with global optimal search (CP-RTM + GOS), Forest Reﬂectance and Transmittance (FRT) RTM inversion with TCC constraint (FRT + TCC), and random forest (RF) algorithm. The results showed that the method proposed in this study (CP_RTM + TCC) yielded the highest estimation accuracy (R 2 = 0.92, RMSE = 0.2) among the four methods used as benchmark, indicating its reasonable ability to assist forest managers to better understand post-ﬁre vegetation regeneration and forest management.


Introduction
Wildfire is a pervasive natural disturbance, which affects the global forest climate services by regulating the spatial distribution of forests, as well as the exchange of carbon, water, and energy between the land and atmosphere [1][2][3]. Burn severity has been defined as the impact of fire on cannot directly simulate the radiative transfer characteristics of three-layer vegetation structure except for the Forest Reflectance and Transmittance (FRT) model [27][28][29]. The two-layer vegetation structure generally includes tree canopy and lower grassland on the soil surface, while the three-layer vegetation structure includes upper tree canopy, middle shrub, and lower grassland. Quan et al. [30] has proved that the coupled RTM has better performance than the single RTM as the coupled RTM can better simulate the radiative transfer characteristics of multi-layer vegetation structure forest.
All of the previously presented methods that estimate burn severity by use of a single optical post-fire remotely sensed image have a defect in that there will spectral confusion between different burn severity levels with the variation of tree canopy cover (TCC). This is because the post-fire dead leaf litter and charcoal on the soil surface will have different contributions to the spectral signal received by remote sensing sensors with the variation of TCC [12]. The charcoal and dead litter on the soil surface will have a stronger contribution to the observed spectral signal for the sparse tree-covered area, while for the dense tree-covered area, the spectral signal of remote sensing sensors received from the soil surface will be shielded by the green tree canopy. The variation of TCC will lead the same burn severity will produce different spectral characteristics and the same spectrum will correspond to different burn severity values. Yin et al. [9] significantly improved the CBI estimation accuracy by using prior TCC information into the forward simulations and backward inversions of FRT model. Because there is no TCC parameter in FRT model, the stand density is used as a proxy parameter to constrain TCC. However, TCC is not only determined by stand density but also canopy diameter. Moreover, the correlation between stand density and TCC has not been validated other than the Yatir forest in Israel [31], and therefore the established empirical relationship cannot be generalized, limiting the application of FRT RTM-based inversion for forest burn severity somewhere else.
Within this context, our objectives are to (1) develop an improved RTM-based method to retrieve burn severity by coupling A two-layer Canopy Reflectance Model (ACRM) and GeoSail models at canopy level, and (2) achieve more accurately burn severity mapping. ACRM is used to model the spectra of vegetation in the middle and low stratum, while GeoSail is used to model the spectrum of the discontinuous upper tree canopy layer. The tree coverage parameter FCOV is one of the input parameters of GeoSail model, therefore it can explicitly account for variations in TCC without a proxy parameter. Our hypothesis is the coupled RTM can simulate the radiative transfer characteristics of three-layer vegetation and has greater potential to accurately estimate burn severity worldwide given it directly accounts for variation in TCC to avoid misclassify burn severity without the need of determining a site-specific stand density-TCC relationship. More accurate estimation of burn severity will contribute to better guide post-fire forest management.

Study Area
The study area is located at Qinyuan County, central China (111 • 58 30"E-112 • 32 30"E, 36 • 20 20"N-37 • 00 42"N) ( Figure 1). The study area is dominated by overstory vegetation of Pinus tabuliformis while the middle vegetation layer is dominated by low Quercus acutissima and Rosa xanthine and the grass over the soil surface are dominated by Gramineous weeds. Under the management of the local government's policy of closing hills for reforestation, the vegetation growth in the study area is very dense, and the TCC in most areas is more than 50%. The height of Pinus tabulaeformis in the study area ranged from 3 to 5 m, which has a positive correlation with fire behavior [32]. Pinus tabulaeformis will become very flammable in dry and hot weather for a long time. At the same time, pine balls will help the fire spread faster with the help of wind. The topography of the study area is rugged, and the altitudes range from 1020 to 1639 m. The average annual rainfall in this region is 656.7 mm. The average annual relative humidity is 65% and the average annual temperature is 8.7 • C. Maximum and minimum precipitations were recorded in July-September and November-March, respectively. From November to March there is little precipitation, and the temperature will continue to rise from about −10 • C to 15 • C. During this period, sparse precipitation and temperature rise will help forest ecosystems accumulate large amounts of dry and flammable fuels. Generally, the fire season begins in March and ends in May.

Data Collection
In March 2019, two fires that caught the attention nationwide broke out in south (fire 1, 14 March) and north of (fire 2, 29 March) Qinyuan County ( Figure 1). Fire 1 was small in scale (approximately burnt 48 ha), but six firefighters died during fire suppression because of sudden changes in wind direction. There were no casualties in fire 2 even though the burn area was larger (approximately 16,788 ha,). The post-fire field survey was conducted from 18 to 23 July 2019. The CBI was measured based on the field criterion proposed by Key and Benson (2006) and used as the ground truth to validate the accuracy of burn severity estimated by the methodology proposed in this study. There was a total of 20 and 143 field CBIs recorded for fire 1 and fire 2, respectively ( Figure 1). Field plots were selected within areas with homogeneous burn severity levels, as recommended by Key and Benson [7].

Sentinel-2A MSI Data
In this study, the Sentinel-2A Multi-Spectral Instrument (MSI) data were used as the satellite data source to estimate burn severity. The Sentinel-2A MSI data has 13 spectral channels in the visible (VIS), near-infrared (NIR), and shortwave infrared (SWIR) spectral ranges. To keep the spatial resolution (20 m) consistency of surface reflectance data, nine out of the 13 spectral bands (bands 2-7, 8a, 11, and 12) were used to estimate burn severity ( Table 1). Considering that NBR can provide more additional spectral information, the post-fire NBR, normalized of bands 8a and 12 were also used for burn severity estimation. To ensure that it is closest to the field survey time, the post-fire imagery was acquired on 30th June 2019. This image is the closest to the sampling date and has the best imaging quality. The Sentinel-2A MSI Level-1C data used in this study were downloaded from the Copernicus Open Access Hub (https://scihub.copernicus.eu/). The atmospheric correction of Sentinel-2A MSI data was conducted by use of the Sen2Cor Tool (version 2.3.1) [33,34]. The reflectance of the Level-1C Bottom-of-Atmosphere (BOA) was converted to Level-2A by using the default parameters. Sen2Cor consists of two main modules: the Scene Classification (SCL) module and the Atmospheric Correction (AC) module [35]. The AC is executed using a set of Look Up Tables (LUTs) generated by libRadtran [36]. AC includes three configuration parameters: aerosol type, atmosphere type, and ozone content. The aerosol type can be set as rural or maritime (the default setting is rural). There are two atmosphere types (mid-latitude summer or mid-latitude winter) that can be set, which will be automatically selected by Sen2Cor according to the geographic location and climatology of the scene [35]. The Ozone content is also automatically selected by Sen2Cor according to the geographic location and climatology.

Landsat Vegetation Continuous Fields Product
The Landsat Vegetation Continuous Fields (VCF) tree cover product contains estimates of the TCC in each 30 m pixel covered by woody vegetation greater than 5 m in height [37]. The product extracts all seven bands of Previous versions of TCC were collected using the Landsat Global Land Survey, but version 4 was created by mining and processing entire Landsat archives, providing improved and accurate representation for specific years (https://landsat.gsfc.nasa.gov/global-30m-landsat-treecanopy-version-4-released/). In this study, the TCC layer of 2015 was used to parameterize and constrain the RTM. The forests in the study area are natural forests with little human intervention, and there has been no forest fire in the study area since 2015, thus the TCC layer in 2015 can be used to represent the situation before the fire. For more information about this dataset please visit http://landcover.org/data/landsatTreecover. The TCC layer of Landsat VCF data was resampled to 20 m matching the spatial resolution of Sentinel-2A MSI data using nearest-neighbor interpolation.

Burn Severity Retrieval Using Coupled RTM
The methodology to retrieve burn severity using coupled RTM is carried out in four steps ( Figure 2): (i) remote sensing data preprocessing, including the atmospheric correction of the Sentinel 2A MSI data and resampling of the Landsat VCF product; (ii) couple RTMs of leaf and canopy level for burn severity retrieval; (iii) forward modeling, including the parameterization of the coupled RTM based on the sensitivity analysis results of previous studies and construction of the LUT [24]; and (iv) Coupled RTM backward inversion procedure by finding the simulated surface reflectance (saved in the LUT) that most similar to the extracted spectrum from satellite remotely sensed data by use of a cost function [38]. The coupled RTM can simulate the radiative transfer characteristics of three-layer vegetation and has greater potential to accurately estimate burn severity worldwide given it directly accounts for variation in TCC to avoid misclassify burn severity without the need of determining a site-specific stand density-TCC relationship.

Model Selection and Coupling
Several vegetation layers are considered for CBI evaluation. Therefore, the selected RTM should have the capacity to simulate the spectrum of multi-layered forests. In this study, the leaf-level RTM PROSPECT [39], canopy-level RTM ACRM [40,41], and GeoSail RTM [42] were coupled to simulate the spectra of the multi-layered forests for burn severity retrieval. The coupled RTM has the capacity of simulating the radiative transfer characteristics of three-layer vegetation and has greater potential to accurately estimate burn severity worldwide given it directly accounts for variation in TCC to avoid misclassifying burn severity without the need of determining a site-specific stand density-TCC relationship.
The PROSPECT model was used to simulate the spectrum of a green and brown needle leaf. Although other RTMs may be more suitable for simulating the reflectance and transmittance of needles than PROSPECT (e.g., LIBERTY model proposed by Dawson et al. [43]), they are not as operational as PROSPECT because it requires more input variables with more complex parameters. Besides, PROSPECT has been successfully applied to the simulation of needle leaves' spectra in many studies [24,44,45]. The PROSPECT model is a leaf optical model, which calculates the reflectance and transmittance of the leaf hemisphere in the range of 400 to 2500 nm in steps of 1 nm. The calculated reflectance and transmittance are expressed as a function of six input parameters: leaf structure parameter, N (unit-less); C ab (µg cm −2 ); Equivalent Water Thickness, C w (g cm −2 ); Dry Matter Content, C m (g cm −2 ); carotenoid content, C ar (µg cm −2 ); and leaf brown pigment, C s (unit-less).
ACRM assumes that the canopy consists of a homogeneous vegetation layer and a thin layer of vegetation on the soil surface. The model is an extension of the homogeneous multispectral canopy reflectance model MSRM [46] and Markov chain canopy reflectance model MCRM [47]. The model has a step size of 1 nm in the 400 to 2400 nm spectral range. This RTM is thus suitable for simulating the spectra of middle and lower canopies. Here, the ACRM was coupled with the GeoSail model by replacing the spectra of the bare soil input of GeoSail and considering the vertical stratification characteristics of vegetation in the study area. The ACRM was used to simulate the spectrum of middle vegetation layers and grass over the soil surface.
The GeoSail model was used to simulate the spectrum of the upper tree canopy layer since it is a geometric model that considers the geometry and shadowing effects of forest stands. The Geosail model described the canopy reflectance of heterogeneous and discontinuous vegetation types with relatively simple input parameters [42] and has been successfully used for burn severity retrieval [24]. The tree coverage parameter FCOV in the GeoSail model can be directly constrained by satellite-derived TCC without proxy parameters such as stand density, therefore the coupled RTM has greater potential than FRT RTM in burn severity retrieval. All the canopy vegetation parameters required by the coupled RTM are listed in Table 2. For the range and default values of model parameters, refer to the user's guide and previous studies [24,27,30]. For the coupled RTM, the diffuse radiation was set as the dominant radiation for the grassland layer by the hypothesis that most of the direct radiation intercepted by the tree canopy. Therefore, the bi-hemispherical reflectance rather than the BRDF of the understory was modeled in this case. Besides, the spectral reflectance of the background originally in the GeoSail is assumed to be the reflectance of Lambertian while the spectra modeled by the PROSAIL (PROSPECT and Scattering by Arbitrarily Inclined Leaves) is non-Lambertian.

Model Parameterization and Forward Modeling
The parameterization of the coupled RTM was based on prior information obtained from field quantitative measurements, qualitative survey, previous studies [24,30], and Sentinel 2A MSI metadata. Only the parameters most sensitive to the change of the target bands according to previous analysis from De Santis et al. [24] and Yin et al. [9] were used, i.e., the fractional cover (FCOV), leaf area index (LAI), and water content (C w ) are most sensitive to the short-wave infrared (SWIR) band. In detail, the LAI and FCOV are the most sensitive parameters for the NIR band. The C ab , LAI, and brown pigments content (C s ) are the most sensitive parameters for VIS bands. Therefore, the sensitive parameters selected for the GeoSail model were FCOV and LAI at canopy level and C w , C ab , and C s at the leaf level, as illustrated in Table 3. In the CBI field criterion, burn severity is evaluated using a set of variables for each vegetation vertical strata, and the CBI increases with the proportion of brown over green leaves. Therefore, the PROSPECT model was used at the leaf level to simulate only one representative spectrum for green and brown reference leaf types. The parameterization for PROSPECT model were N = 2.5, C ab = 70 µg cm −2 , C w = 0.048 g cm −2 , and C s = 0.2 for green leaf [48] while to N = 2.5, C ab = 20 µg cm −2 , C w = 0.008 g cm −2 , and C s = 1.5 for brown leaf [49]. The simulated green and brown leaf spectrum based on the PROSPECT model was shown in Figure 3. The ACRM model was calibrated following Yin et al. [9]. For the middle vegetation layer, the sensitive parameters are LAI of middle vegetation and a series of biochemical parameters determining leaf color (N, C ab , C w , and C s ); for the lower vegetation layer, the sensitive parameters are LAI of lower vegetation and biochemical parameters determining leaf color (N, C ab , C w , and C s ) ( Table 3).  Considering the uncertainty caused by unrealistic combinations of the input model variables, a supervised simulation scheme was executed to reduce the error because it is better than the full range variations [24,26]. The parameter FCOV (range from 0 to 0.8) in the GeoSail model was set as the only free variable to generate the LUT. When the coupling model is forward simulated, the spectral curves corresponding to different TCC conditions can be generated with the variation of FCOV. A combination example of the model input parameters used to generate LUT is listed in Table 3. The background shortly after wildfires usually consists mainly of soil and charcoal [24]. Therefore, three reference spectra were considered as the background layer for the model input: soil, dark charcoal (DCH), and light charcoal (LCH; a mixture of charcoal and ash) (Figure 4). The soil spectra were measured on sandy soil with moderate moisture, while the spectrum of dark and light charcoal was measured with GER 2600 field spectro-radiometer (Geophysical & Environmental Research Corporation, Millbrook, NY, USA) by De Santis et al. [24]. The simulated CBI value variation from 0.5-3 with steps of 0.1. The burn severity can be divided into low, moderate, and high level by use of CBI threshold values of 2 and 2.5 [9]. For low burn severity level, only the grassland and shrub on the soil surface were burned. For moderate burn severity level, the underlying grass and shrubs were completely burned and the understory of the tree was slightly burned, but the trees still retained the green canopy. For high burn severity level, the low and middle vegetation strata were completely burned and the crown was severely burned. For the part of CBI lower than 2.5, only the DCH scenario was simulated. For the part of CBI higher than 2.5, the organic matter in the soil and all vegetation strata is burned, the substratum can be charcoal or ash, therefore both DCH and LCH scenarios are simulated. Finally, there are 352 combinations of input parameters used to construct the LUT for burn severity retrieval.

Coupled RTM Inversion
In this study, the LUT was used as the inversion algorithm to retrieve burn severity using Sentinel-2A MSI satellite remotely sensed data. We derived burn severity using the information on TCC to constrain the backward inversion procedure. For each pixel in the burnt area, the Landsat VCF product was served to select the corresponding LUT. The model simulated spectral reflectance in the LUT corresponding to a specific TCC range was subsequently compared to the extracted surface reflectance of the pixel by use of a cost function. The Spectral Angle Mapper (SAM) was used as the cost function to find the model simulated spectral reflectance that was most similar to the spectra reflectance extracted from remotely sensed data. SAM is a supervised classification technology based on pixels. SAM solves spectral similarity by calculating the spectral angle (SA) between two spectral vectors [50,51]. This angle is independent of the length of vectors and is therefore insensitive to the illumination or albedo effects (Equation (1)). Consequently, SAM can eliminate most of the errors introduced by terrain relief [23].
The i represents the band number, nb represents the number of the bands, and t and r represent the simulated and observed spectrum, respectively. SA between the simulated and observed reflectance was calculated as the cost function to estimate the biophysical and biochemical variables of the vegetation that corresponding to certain CBI values. The R square (R 2 ), root mean square error (RMSE), and slope between the measured and estimated CBI were used as the accuracy evaluation index. The backward inversion process in this study was executed by use of all nine Sentinel-2A MSI bands and the coupled RTM with TCC constraint (labeled as CP-RTM+TCC). In addition to the nine spectral bands of Sentinel-2A MSI data mentioned above, the post-fire NBR was also used as an additional vector to provide more spectral information and improve the accuracy of burn severity retrieval.

Methodological Comparison
To compare and evaluate the merit of our proposed method (labeled as CP-RTM+TCC), we derived burn severity using three other existing methods. In the first method, the coupled RTM is inverted using a global optimal search and without TCC constraint (labeled as CP-RTM+GOS). The second one is the method used by Yin et al. [9] (labeled as FRT+TCC). In this method, the burn severity was estimated using the FRT model with TCC constraint. The third method uses a random forest (RF) algorithm. The RF classifier is particularly useful for burn severity mapping with promising accuracy [18,19,52]. In the RF method, the independent variables are the nine spectral bands and the post-fire NBR index of Sentinel 2A MSI data, and the dependent variable is CBI. First, the field measured plots were randomly divided into three samples (each sample comprised about 54 CBI measurements). Second, two samples (i.e., the training sample) were used for RF model calibration, and the other sample (i.e., the testing sample) was used for the validation. Third, the process was repeated three times to adequately test the performance of the RF method [53]. To minimize the statistical error caused by each run of the RT model, we set the number of trees to 1000 and run 100 times to train each RF model.

Influence of TCC on the Spectral Response of Burn Severity
The surface reflectance and post-fire NBR of all the field plots were extracted from Sentinel-2A MSI data from locations corresponding to our sampling locations and averaged by burn severity level and TCC (<20%, 20%-30%, 30%-50%, and >50%) to prove the influence of TCC on the spectral response and the confusion between different burn severity levels ( Figure 5). Under this circumstance, the burn severity was divided into low, moderate, and high level by use of CBI threshold values of 2 and 2.5 [9]. Figure 5(a1,b1), and Figure 5(c1) represent the extracted mean reflectance from all field plots without constraining TCC. The averaged reflectance of all field plots obscured the reflectance variation between different TCC ranges. By analyzing the spectrum of three different levels of burn severity under different TCC, it can be found that the high-level burn severity is not sensitive to the change of tree coverage, while the medium level burn severity shows little sensitivity to the change of tree coverage, and the low-level burn severity is most sensitive to the change of tree coverage. This was also observed from the spectral angle (SA) between the spectra of lowest and highest TCC within the same burn severity. The SA was 0.68 ( Figure 5(a2,a5)), 0.33 ( Figure 5(b2,b5)), and 0.27 ( Figure 5(c1,c5)) for low, moderate, and high burn severity, respectively. For TCC over 50%, the spectral characteristic of the high burn severity exhibited similar characteristics to the moderate and low burn severity with the lower TCC. For a detailed analysis of the influence of TCC on the inversion of burn severity, please refer to Yin et al. [9]. The results given in Figure 5 indicate that if the TCC is not considered in the inversion of burn severity, spectral confusion will be caused by the change of TCC between different burn severity levels.

Evaluation of Burn Severity Estimates
A total of 163 CBI values recorded at two sample sites were used to validate the retrieved CBI. The CBI was estimated based on the coupled RTM constrained by TCC information provided by satellite estimates (CP_RTM+TCC). The validation results at both sites showed a high correlation between the observed and simulated CBI (R 2 equal to 0.96 for study site 1 (Figure 6a), 0.91 for study site 2 (Figure 6b), and 0.92 for two study sites (Figure 6c)) and an almost 1:1 linear fit (slope:~1; constant:~0). The RMSE is also very low (0.15 for study site 1, 0.21 for study site 2, and 0.20 for two study sites). When comparing the method proposed in this study with the burn severity estimations obtained with CP_RTM+GOS, FRT+TCC, and RF methods (Table 4 and Figure 7), the method proposed in this study but without TCC constrains (CP_ RTM+GOS method) presents slightly worse performance than FRT+TCC and RF. Figure 7b also shows that the estimated results of the FRT+TCC method yield a relatively higher accuracy (R 2 equal to 0.82, RMSE equal to 0.32, and slope equal to 0.81). As can be indicated from Figure 7c that the validation accuracy (R 2 = 0.76, RMSE = o 0.35 and slop = 0.7) based on the RF method is higher than CP_RTM+GOS but lower that of the FRT+TCC method. Overall, the CP_RTM+TCC method ( Figure 6) performed best (R 2 equal to 0.92, RMSE equal to 0.2 and slop equal to 0.99) among these four methods.

Burn Severity Mapping
The focus of this study is burn severity estimation, so the dNBR as well as visual interpretation were used to extract burn area. Burn severity maps for the two sample sites of Qinyuan country were produced based on the burn area extraction and the CP_RTM+TCC method ( Figure 8). The high-level burn severity range (CBI ≥ 2.5) in site 1 is mainly distributed in the western area because the elevation of the western area of the site is relatively high and its main vegetation type is dense and dry Pinus tabulaeformis forest, which is quite flammable. The eastern part corresponds to the foot and hillside with relatively low elevation, the vegetation types are mainly yellow rose, low bush, and dense herbaceous vegetation, and the burn severity type is mainly of medium (2 ≤ CBI < 2.5) and low level (CBI < 2). For study site 2, the main vegetation type is also Pinus tabulaeformis. The middle and high-level burn severity are mainly distributed at the top of the mountain with relatively high altitude, while the low-level burn severity is mainly distributed at the foot of the mountain with relatively low and flat terrain. From the spatial distribution map of burn severity in study area 2, we can also see that there are more areas with CBI over 2.8 and below 2.0, but fewer areas with CBI between 2.0 and 2.8 than in study site 1. Figure 8. Burn severity maps for study sites 1 (left) and 2 (right) of Qinyuan fires achieved using the CP_RTM+TCC method proposed in this study. In order to show the spatial distribution characteristics of burn severity in the two study sites more finely, CBI was divided into six grades which vary from 0 to 3.
In the field measurement, in addition to the photos used to estimate the CBI value of the sample plots, we also selected some areas with a wide view to take panoramic photos for quantitative verification of burn severity inversion results. The qualitative validation was carried out by comparing the field landscape panoramic photos and the spatial distribution map of burn severity (Figure 9). A total of four validation points were selected: two (Figure 9b,c) of them point to the valley, and the other two (Figure 9a,d) point to the hillside. The qualitative accuracy validation results indicated that the burn severity map and the landscape panoramic photos from the field survey showed good spatial distribution consistency. This provided further evidence for the reliability of our burn severity estimates.

Discussion
In this study, the burn severity of two forest fires was estimated using a coupled RTM and Sentinel-2A MSI reflectance data. ACRM assumes that the canopy consists of a homogeneous vegetation layer and a thin layer of vegetation on the soil surface, which is thus suitable for simulating the spectra of middle and lower canopies. The GeoSail model is a geometric model that considers the geometry and shadowing effects of forest stands and was used to simulate the spectrum of the upper tree canopy layer. The ACRM was coupled with the GeoSail model by replacing the spectra of the bare soil input of GeoSail and considering the vertical stratification characteristics of vegetation in the study area. The coupled RTM can simulate the radiative transfer characteristics of three-layer vegetation and has greater potential to accurately estimate burn severity worldwide given it directly accounts for variation in TCC to avoid misclassify burn severity without the need of determining a site-specific stand density-TCC relationship.
Compared with the three existing methods in other studies, the proposed method (CP_ RTM+TCC) in this study performed best. The CP_ RTM+GOS method performed worst because if TCC is not constrained, there is serious spectral confusion between the different burn severity levels. This is similar to the conclusion of Yin et al. [9], who found spectral confusion between different burn severity levels when the influence of TCC was not considered in the RTM forward and backward inversion procedure. Accuracy comparison results also showed that the FRT+TCC method yields a relatively higher accuracy. This is because the influence of tree canopy cover is considered in the FRT+TCC method. However, the accuracy was still lower than that of this study. This is because the stand density was used as a proxy parameter to constrain the TCC in the FRT+TCC method, but stand density is not the only determining parameter of tree coverage. Furthermore, the statistical correlation between stand density and the lack of universality in TCC occurs because the empirical correlation has not been validated other than the Yatir forest in Israel. Therefore, the indirect constraint of tree canopy cover through stand density will produce matching errors, which will affect the estimation accuracy of CBI. In essence, the RF algorithm estimates the burn severity according to the characteristics of spectral reflectance, so it cannot get rid of the influence of TCC. Therefore, the estimated accuracy based on the RF method is still lower than the method proposed in this study.
The spatial distribution of burn severity in the two study sites shows obvious boundary characteristics. Most regions have CBI values above 2.8 or below 2.0, but few regions between 2.0 and 2.8. Combined with field sampling, we analyzed two main reasons for this phenomenon. First, the interval between the fire and the field sampling was nearly 4 months and there was rainfall during this period. The vegetation recovery ability in the low and medium burn severity area was stronger than that in the high burn severity area. Most of the trees are burned completely under the high burn severity level, which hinders quick recovery, and the low and medium burn severity is only affected by the medium degree, supporting recovery faster speeds, which is consistent with the conclusion of Meng et al. [54]. On the other hand, the relationship between vegetation recovery rate after the fire and burn severity is not linear, but is a saddle curve [54]. This curve shows that the forest has certain fire resistance. For example, some trees under fire stress have a short-term loss of canopy leaves but do not die. Therefore, their canopies will recover in a short time after the fire. However, with the increase of forest burn severity, tree survival will decrease as canopies cannot recover from the damage. The areas affected by a moderate burn severity did not cause fatal damage to trees, therefore the recovery rate was the fastest, with quick recovery to the state of lower burn severity level after the fire. In most cases, the damage to trees caused by high burn severity is irreversible, implying slow recovery. Second, field sampling saw that the high burn severity level was mainly distributed in areas with high altitude and high slope. During the period from the fire (29 March) to the field sampling (18 July), there were several precipitation events (>83 mm), and the temperature continued to rise (from 6.9 to 26.3 • C) ( Figure 10); this created good conditions for the post-fire recovery of vegetation in the study area. However, because of the fire-related decrease of water storage capacity in mountain forests, there was a slower recovery rate of trees at the top of mountains and in steep areas compared to valleys and relatively flat terrain. In the process of field survey, we also found that the spatial distribution of burn severity has a high correlation with terrain characteristics. For example, the areas with high burn severity are mostly distributed on the top of the mountain and steep hillside, which indicates that terrain factor is the important drive factor of burn severity. Future work will be focus on analyzing the relationship between post-fire burn severity and pre-fire environmental drive factors. The meteorological data before the forest fire 2 in Qinyuan county were also analyzed. Figure 10 gives the daily average temperature and precipitation data of Taiyuan meteorological station nearest to Qinyuan county (longitude: 112.35 • E, latitude: 37.37 • N) from 1 October 2018 to 31 July 2019. There was almost no precipitation in the six months before the fire and the temperature continued to rise in the three months before the fire, which provided sufficient external conditions for the fire. This is also evident in the driving effect of meteorological factors on fire. Furthermore, due to the local policy of closing the mountain for forest cultivation, many areas in sample site 2 were completely closed, resulting in the accumulation of a large number of fuel load, including live fuel, such as Pinus tabulaeformis and other shrubs; herbaceous vegetation; and dead fuel represented by litter and humus on the ground. The accumulation of fuel load is also one of the important causes of fire ignition. Meteorological variables are one of the most important factors affecting post-fire vegetation regeneration. Vegetation regeneration and growth are generally highly depending on favorable postfire meteorological conditions, such as conditions with higher moisture and proper temperature [55]. The several precipitation events and the average daily temperature continued to rise in the study area after the fire created good conditions for vegetation regeneration.

Conclusions
Burn severity is of vital importance for post-fire forest management and will assist in post-fire vegetation regeneration. In this study, the burn severity of two forest fires in Qinyuan was estimated using coupled RTM and Sentinel 2A MSI data. Traditional optical satellite remotely sensed data-based methods that do not take the variation in TCC into account will lead to spectral confusion in burn severity retrieval. Yin et al. [9] significantly improved the CBI estimation accuracy by using prior TCC information into the forward simulations and backward inversions via the stand density parameter of the FRT model. However, TCC is not only determined by stand density but also canopy diameter. Moreover, the correlation between stand density and TCC is not universal, limiting the application of FRT RTM-based inversion for forest burn severity somewhere else. In this study, ACRM was coupled with the GeoSail model, by replacing the spectra of the bare soil input of GeoSail and considering the vertical stratification characteristics of vegetation in the study area. The tree coverage parameter FCOV in the GeoSail model was directly constrained by TCC without proxy parameters, and therefore has greater potential to accurately estimate burn severity worldwide. The verification of the accuracy of results using four methods (CP_RTM+TCC, CP_RTM+GOS, FRT+TCC, and RF) showed that the method proposed in this study (CP_RTM+TCC) has the highest burn severity estimation accuracy. The more accurate estimation of the burn severity obtained from the approach presented in this study will help land managers to better understand post-fire vegetation regeneration and help in post-fire forest management. Now that a new method for accurately estimate fire severity has been found, future work will be focus on analyzing the relationship between post-fire burn severity and pre-fire environmental drive factors.