Evaluating the Best Spectral Indices for the Detection of Burn Scars at Several Post-Fire Dates in a Mountainous Region of Northwest Yunnan, China

Abstract

Remote mountainous regions are among the Earth’s last remaining wild spots, hosting rare ecosystems and rich biodiversity. Because of access difficulties and low population density, baseline information about natural and human-induced disturbances in these regions is often limited or nonexistent. Landsat time series offer invaluable opportunities to reconstruct past land cover changes. However, the applicability of this approach strongly depends on the availability of good quality, cloud-free images, acquired at a regular time interval, which in mountainous regions are often difficult to find. The present study analyzed burn scar detection capabilities of 11 widely used spectral indices (SI) at 1 to 5 years after fire events in four dominant vegetation groups in a mountainous region of northwest Yunnan, China. To evaluate their performances, we used M-statistic as a burned-unburned class separability index, and we adapted an existing metric to quantify the SI residual burn signal at post-fire dates compared to the maximum severity recorded soon after the fire. Our results show that Normalized Burn Ratio (NBR) and Normalized Difference Moisture Index (NDMI) are always among the three best performers for the detection of burn scars starting 1 year after fire but not for the immediate post-fire assessment, where the Mid Infrared Burn Index, Burn Area Index, and Tasseled Cap Greenness were superior. Brightness and Wetness peculiar patterns revealed long-term effects of fire in vegetated land, suggesting their potential integration to assist other SI in burned area detection several years after the fire event. However, in general, class separability of most of the SI was poor after one growing season, due to the seasonal rains and the relatively fast regrowth rate of shrubs and grasses, confirming the difficulty of assessment in mountainous ecosystems. Our findings are meaningful for the selection of a suitable SI to integrate in burned area detection workflows, according to vegetation type and time lag between image acquisitions.

1. Introduction

With the opening of the Landsat archive in 2008, nearly four decades of moderate resolution images of the Earth have been made available to the global public. Scientists from a vast range of fields are using these data for research purposes and to develop all sort of applications at local to global scales. In recent years, thanks to increasing storage capabilities, local computing power, and cloud computing technology, the analysis of dense time series has become accessible to ordinary users at reduced costs, catalyzing the development of algorithms to systematically process the data and extract targeted information [1,2]. In order to monitor specific land surface dynamics such as natural or human-induced disturbances, land use/land cover change, or urban expansion, a common approach involves the processing and transformation of the raw spectral information contained in remotely sensed images, resulting in the enhancement of the target phenomena. For this purpose, spectral indices (SI) based on the most sensitive spectral bands while minimizing noise have been designed and integrated in change detection workflows to highlight particular land change events and evaluate their temporal trends [3,4,5,6]. In the case of fire disturbance, several SI initially intended to monitor vegetation productivity, such as the Normalized Difference Vegetation Index or the Enhanced Vegetation Index, as well as specially tailored SI such as Normalized Burn Ratio and Burn Area Index, have been employed to detect burned areas and assess the degree of impact of the disturbance using different satellite sensors (some examples: [7,8,9,10]). More recently, these SI have been incorporated in time series analysis in order to systematically detect burned areas and monitor long-term vegetation recovery [11,12,13,14,15,16,17].

To help researchers in the selection of the most suitable SI for a given situation, several authors performed comparative studies in different ecosystems and using different satellite sensors. In Mediterranean-like ecosystems characterized by chaparral, shrubs, and pine vegetation, several SI were tested on multitemporal sets of Landsat Thematic Mapper (TM), Enhanced Thematic Mapper (ETM+), NOAA Advanced Very High Resolution Radiometer (AVHRR), Moderate Resolution Imaging Spectroradiometer/Advanced Spaceborne Thermal Emission and Reflection Radiometer (MODIS/ASTER) Airborne Simulator, and the Chinese Huanjing (HJ) [7,18,19,20,21,22], while in boreal forests, grasslands, and tundra, comparative assessments have been done mainly using Landsat [23,24,25] and Sentinel-2 [26]. The latter study included other vegetation types, such as tropical and subtropical woody savannas, which have also been assessed using TM [27]. A particular study from Schepers et al. [28] tested SI capabilities in the heathland shrubs of Belgium using Airborne Imaging Spectroscopy. All the referenced studies evaluated SI capabilities for the discrimination of burned areas and for mapping burn severity short time after the fires occurred. With a focus on post-fire vegetation recovery patterns, other authors compared different SI responses at multiple post-fire dates using time-series approaches [29,30,31,32]. Summarizing the results of this comprehensive volume of literature, performance divergences were found in different ecosystems and in the time lag between date of the burn and the post-fire image used for the assessment. However, certain patterns have been observed and interpretations suggested. For example, the SI that integrate the shortwave infrared domain of the electromagnetic spectrum, corresponding roughly to wavelengths between 1100 and 3000 nanometers, tended to be more suitable for long-term vegetation recovery than for immediate burn-scar mapping, thanks to their sensitivity to forest structure and moisture content. On the contrary, SI relying on visible and near infrared light could accurately map burned patches if employed immediately after the fire event, when the magnitude of change in vegetation’s chlorophyll content was the highest [29,31,32]. In general, even if the famous Normalized Burn Ratio [33] was often considered the best SI, our literature review founds that there’s no single SI that constantly excels and over performs the others in every ecosystem, scale and time lag conditions, both for burned area mapping and severity assessments, as well as for vegetation recovery characterization.

Of the tested vegetation types, Mediterranean pine/shrubs, semi-dry savannas, and boreal forests are among the most fire prone ecosystems on Earth, explaining the higher amount of published research targeting these regions. Much less can be found about impacts of fire in mountain environments, although they are extremely important for their ecological value while being highly vulnerable to climate change and natural disasters (http://www.fao.org/mountain-partnership, last accessed on 9 March 2018). Remoteness, contrasting topography and patchy landscapes that characterize alpine regions dictate harsh living conditions and access difficulties, limiting field research opportunities and presenting particular challenges for the monitoring of disturbances using remote sensing approaches [34]. Medium to high resolution imagery is required to assess changes in such complex ecosystems and frequent cloud cover fairly affects the quality of the images. Topographic effects are known to introduce distortions and artifacts that need to be corrected during the preprocessing stage of image analysis, together with routine radiometric and atmospheric corrections. Moreover, rugged terrain and sun angle create shades that obscure some portions of the image, highly affecting reflectance values. Several topographic correction methods have been developed and tested, and it is in general suggested to include them in preprocessing operations [35,36], as well as using spectral band ratios which are less affected by topographic effects [34]. However, in very rugged terrain and low sun elevation conditions, the performance of these correction methods decreases drastically [37,38].

In this study, we focus on the mountainous region of northwest Yunnan, China, where wildland fires occur frequently during the dry and windy season, between December and May [39,40]. According to official statistics, more than 99% of forest fires in southwest China including Yunnan are caused by anthropogenic ignition, often due to the widespread use of fire in agriculture [41]. The main peculiarity of forest fires in this region is their relatively small size, in general 100 to 300 hectares and rarely above 1000 hectares, which heavily impact the performance of the existing burned area products commonly used in regional and global assessments of fire impacts, such as MODIS MCD64A1 and ESA Fire_CCI [42]. Therefore, efficient burned area extraction algorithms require image data of a finer resolution, such as 30 m Landsat imagery, and the use of suitable SI. However, cloud-free images are not always available in this region to a point that dense time series are impossible to construct. Time lag between acceptable scene acquisitions can reach several months or even years, further complicating burn scar detection. In this context, our study aims to address the common difficulties encountered by analysts working in mountainous areas who need to deal with low image frequencies and several vegetation types coexisting in the same region. We selected eleven among the most popular SI widely used for burned area and burn severity assessments and analyzed the spectral trajectories of burned and unburned pixels over time in northwest Yunnan’s four major vegetation types. The permanence of the spectral signal that characterize burned land was assessed following a 1-year temporal sampling design and the SI that better highlight the signal were identified at each post-fire year. Our approach can be associated with vegetation recovery evaluations because the detectability of the residual signal of burned areas is strictly related to vegetation regrowth. The next section explains in detail how we adapted existing conceptual frameworks and tailored the metrics to better serve our research question.

2. Data and Methods

2.1. Description of the Study Region

Northwest Yunnan (NWY) is situated in the Yunnan province of China, between approximately latitude 24.5°N and 29.5°N, and longitude 98°E and 101.5°E (see Figure 1). Administratively, it is divided into four prefectures: Dali, Lijiang, Nujiang and Diqing, covering a total area of more than 86,000 km². NWY lies in a peculiar geological region at the transition between the Qinghai-Tibet Plateau and the Yunnan-Guizhou Plateau. It is a mountainous region characterized by a very rugged terrain, forming high mountain peaks that reach more than 6000 m, and deep valleys carved by four major Asian rivers; the Dulongjiang (a tributary of the Irrawaddy), the Nujiang (upper Salween), the Lancangjiang (upper Mekong) and the Jinshajiang (upper Yangtze). These exceptional bio-physical features interact with the seasonal East Asian monsoon, conferring the region a wide variety of climates, landscapes and natural habitats. Major vegetation types include alpine needleleaved forests dominated by species of Pinus, Picea and Abies; subtropical evergreen broadleaved forests composed of Lithocarpus, Castanopsis and Quercus species; varied shrublands with succulent Euphorbia royleana and Opuntia monacantha as indicator species but also hosting young Pinus yunnanensis, Rhododendron and Magnolia species. Moreover, the region hosts alpine meadows, grasslands and arid savannas in dry valleys [43,44].

2.2. Selection of Burn and Reference Plots and Sampling Methodology

Official and accurate data about forest fires such as location and date are unfortunately not available for our study region. Consequently, the burned plots were manually drawn within natural or anthropogenic burned areas (no experimental burns) by the means of visual interpretation of Landsat TM scenes and MODIS fire products (MCD45A1 and MCD64A1) at various dates (http://modis-fire.umd.edu/, last accessed on 9 March 2018). Burned plots followed the burned area perimeter using particular care in avoiding peripheral pixels to keep the core of the burned area. Consequently, shape and size were different among plots. In a first stage, a total of 50 fires were identified from Landsat images (path/row 131/042 and 132/041), 25 of which were confirmed by the MODIS products. A second round of selection was performed in order to obtain 12 burn plots equally representing 4 dominant land cover classes (3 plots per class). The land cover classes were extracted and aggregated from the European Space Agency Climate Change Initiative (ESA CCI) Landcover product (https://www.esa-landcover-cci.org/, last accessed on 9 March 2018), and were: needleleaved forests, broadleaved forests, shrublands, and grasslands. The basic requirements for the 12 final burn plots consisted in having a relatively narrow time lag between the pre-fire and the immediate post-fire image dates, and having good quality, 1-year temporal interval Landsat scenes in all of the five years following the burn. Images affected by clouds were also considered as long as they didn’t cover the selected burned plots. The final dataset was composed of a 7-images time series for each burn plot. The frequent presence of clouds in the study region was the main problem in both selection phases and reduced drastically the number of acceptable burn plots. All Landsat scenes were retrieved from the United States Geological Survey EarthExplorer repository (https://earthexplorer.usgs.gov/, last accessed on 9 March 2018). The scenes were already preprocessed by the provider according to Level-2 surface reflectance standards using the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) radiometric calibration and atmospheric correction algorithms [45]. Level-2 images account for most of the artifacts introduced by rugged terrain but topographic shades persist in the images. Our tests of existing topographic correction methods gave poor results due to the excessive steepness of most of the slopes in the study region. We therefore chose to rely on the slight normalizing effect of band ratios, which is used by most of the SI, and on the selection of yearly Landsat scenes at same periods of the year, where sun elevation, sun azimuth and vegetation phenology are similar. Average fire severity for each burn plot was estimated using Delta NBR between the pre-fire and post-fire images. According to the burn severity assessment protocol used by the USDA Forest Service (FireMon documentation: https://www.fs.fed.us/rm/pubs/rmrs_gtr164/rmrs_gtr164_13_land_assess.pdf, last accessed on 9 March 2018), five plots were high severity burns, three plots were moderate-high severity burns, two plots were moderate-low severity burns, and two plots were low severity burns. Details of the 12 final burn plots are described in

Table 1. Landsat TM time series for the 12 selected burn plots. Locations of the burn plots are indicated with the geographic coordinates of their center, WGS84, decimal degrees. ESA LC: land cover class extracted from the European Space Agency Landcover product; DOY: day of year (Julian day); LS_pre, LS_post, LS_1, …: Landsat scene date in yyyymmdd format at pre-burn, post-burn, 1,2, … years after burn, respectively.

Plot ID	Lat	Long	ESA LC	Size Sq. km	MODIS Year, DOY	LS_pre	LS_post	LS_1	LS_2	LS_3	LS_4	LS_5
Plot_01	27.74	98.69	broad	2.15	2006 007	20051028	20060201	20070220	20080310	20090209	20100212	20110303
Plot_02	27.54	99.31	broad	1.26		19990214	19990302	20000217	20010118	20020105	20030108	20040111
Plot_03	27.41	99.02	shrub	0.7	2006 005	20051028	20060201	20070103	20080310	20090209	20100212	20110303
Plot_04	27.46	99.85	needle	0.74	2005 334	20051113	20060201	20070103	20080310	20090209	20100212	20110303
Plot_05	26.69	99.12	shrub	7.25	2006 011	20051028	20060201	20070103	20080310	20090209	20100212	20110303
Plot_06	26.55	100.35	broad	1.1		19980204	19980409	19990223	20000414	20010212	20020303	20030306
Plot_07	25.87	100.08	shrub	2.08	2003 065	20030218	20030306	20040221	20050207	20060226	20070301	20080216
Plot_08	25.87	100.77	needle	0.33		20030306	20030509	20040425	20050428	20060517	20070301	20080404
Plot_09	26.42	100.84	needle	3.43		19950212	19950316	19960302	19970201	19980204	19990223	20000210
Plot_10	26.11	100.37	grass	2.28		20021029	20030101	20040104	20050106	20060125	20061227	20080216
Plot_11	26.28	100.44	grass	1.67		20030306	20030322	20040221	20050207	20060226	20070301	20080216
Plot_12	26.28	100.58	grass	1.44		20030101	20030218	20040104	20050106	20060125	20061227	20080216

.

Thirty random points (samples) were created inside each burn plot and thirty control points of similar land cover classes and similar slope/aspect were created in the near proximity of each burn plot, to be used as unburned reference samples. Slope/aspect similarity was considered in order to further minimize the bias introduced by topography. Control points were within 4 km except for burn plot number 5 for whom, due to the presence of snow and clouds, suitable sample points could be found only 35 km away. In summary, every plot was characterized by a dominant land cover class, was composed of 30 burn samples and 30 external unburned samples, and each sample was assigned its own land cover class. In total, 720 sample points were used in our analysis.

2.3. Description of Spectral Indices and Processing

Among the wide range of existing spectral indices and image transformation techniques described in the literature, we selected 11 different approaches developed for the Landsat sensors and commonly used for the assessment of burn severity and the extraction of burned areas. SI equations were applied to each Landsat scene in the time series (see Table 1). Finally, SI values from the newly created layers were extracted at the location of each of the 720 sample points. All operations were performed using the ENVI software from Harris Geospatial (https://www.harris.com/solution/envi, last accessed on 9 March 2018) and QGIS (https://www.qgis.org/en/site, last accessed on 9 March 2018).

SI equations, literature references and all abbreviations used in this section are listed in

Table 2. List of spectral indices and image transformation techniques. Blue = visible blue Landsat Thematic Mapper (TM) band 1; Green = visible green Landsat TM band 2; Red = visible red Landsat TM band 3; NIR = near infrared Landsat TM band 4; sSWIR = shorter shortwave infrared Landsat TM band 5; lSWIR = longer shortwave infrared Landsat TM band 7.

Index Full Name	Abbreviation	Equation	Reference
Normalized Difference Vegetation Index	NDVI	$\frac{Red-NIR}{Red+NIR}$	[46]
Global Environmental Monitoring Index	GEMI	$\gamma \left(1-0.25\times \gamma \right)-\frac{Red-0.125}{1-Red}$ where $\gamma =\frac{2\left(NI{R}^2-Re{d}^2\right)+1.5\times NIR+0.5\times Red}{NIR+Red+0.5}$	[47]
Burn Area Index	BAI	$\frac{1}{{\left(0.1-Red\right)}^2+{\left(0.06-NIR\right)}^2}$	[9]
Normalized Burn Ratio	NBR	$\frac{NIR-lSWIR}{NIR+lSWIR}$	[33,48,49]
Normalized Difference Moisture Index	NDMI	$\frac{NIR-sSWIR}{NIR+sSWIR}$	[50]
Burned Area Index Modified–LSWIR	BAIML	$\frac{1}{{\left(NIR-0.05\times NIR\right)}^2+{\left(lSWIR-0.2\times lSWIR\right)}^2}$	[51]
Burned Area Index Modified-sSWIR	BAIMs	$\frac{1}{{\left(NIR-0.05\times NIR\right)}^2+{\left(sSWIR-0.2\times sSWIR\right)}^2}$
Mid Infrared Burn Index	MIRBI	$10\times lSWIR-9.8\times sSWIR+2$	[52]
TassCap Brightness	BRI	$0.2043\times Blue+0.4158\times Green+0.5524\times Red$ $+0.5741\times NIR+0.3124\times sSWIR+0.2303\times lSWIR$	[53,54]
TassCap Greenness	GRE	$-0.1603\times Blue-0.2819\times Green-0.4934\times Red$ $+0.794\times NIR-0.0002\times sSWIR-0.1446\times lSWIR$
TassCap Wetness	WET	$0.0315\times Blue+0.2021\times Green+0.3102\times Red$ $+0.1594\times NIR-0.6806\times sSWIR-0.6109\times lSWIR$

and in its legend. The following sections describe the selected indices, divided in four groups according to the spectral bands used for their design.

2.3.1. SI Using Visible and Near Infrared Domains of the Electromagnetic Spectrum

One of the most known and widely used spectral indices, NDVI takes advantage of green vegetation’s strong absorption of visible red light, in contrast with its high reflection of near-infrared light. The damage to vegetation caused by fire results in a significant decrease of the NDVI. Because of its applicability to a wide range of sensors, including those lacking SWIR bands, its simplicity and relatively good performance, NDVI has been used extensively to map burned areas, assess burn severity and monitor vegetation recovery [7,14,55,56,57]. Similar to NDVI, GEMI uses the visible red and near infrared domains but was designed following a non-linear approach in an attempt to reduce undesired atmospheric effects. Although its application to burn scar detection was mainly done with MODIS and SPOT sensors [58,59,60], we decided to include it in our Landsat TM study and evaluate its potential. BAI was initially designed for NOAA-AVHRR images in Mediterranean environments but has also been tested with Landsat TM [18]. It focuses on the specific spectral signal of charcoal in areas affected by fire and is computed from the spectral distance from each pixel to a reference spectral point, where recently burned areas tend to converge [18]. Reference reflectance values of the convergence point, set to 0.1 for red and 0.06 for NIR, were defined based on literature and several sets of satellite sensor images analyzed by the authors [9].

2.3.2. SI Using Near Infrared and Shortwave Infrared Domains of the Electromagnetic Spectrum

Initially designed for burned area extraction, NBR is the most popular spectral index used for burn severity assessments with different sensors in several ecosystems around the world. In numerous comparative analyses, NBR proved to be one of the most efficient SI [14,19,21,23,28]. In time series analyses, NBR showed good correlations with field-based composite burn index scores several years after a fire, thus representing an efficient tool for vegetation recovery monitoring [29,30]. Unlike the three previously described SI, NBR uses the band pair SWIR and NIR instead of the visible red, which has shown to be sensitive to reflectance changes caused by fire. NDMI is very similar to NBR, but uses the shorter SWIR band instead of the longer one, which has revealed to be equally sensitive to wet/dry content of soils and vegetation. However, it has been rarely tested for burned area discrimination [11,19,29]. BAIML and BAIMs are a modified version of the BAI tailored for the MODIS sensors which offer SWIR channels, but they are also used with Landsat images [11,27,61]. Following the same logic as the one defined for BAI, convergence point values for NIR and SWIR bands were chosen by their creators based on the analysis of MODIS time series, and were set to 0.05 for NIR and 0.2 for SWIR. These values were kept unchanged for the long and short versions of the SWIR used by the Landsat TM instrument.

2.3.3. SI Using Different Shortwave Infrared Domains of the Electromagnetic Spectrum

Designed for arid savannah and shrub ecosystems, MIRBI combines the two SWIR bands because of their better separability in burned areas. Several authors working in these vegetation types have reported its good performance, often superior to the widely used NBR and other SI [25,28].

2.3.4. Image Transformation Techniques Using Multiple Bands

The Tasseled Cap Transformation (TassCap) is a conversion of the original image data to a new coordinate system with a new set of orthogonal axes by applying different weights to the bands. Studies have shown how using TassCap in replacement or in addition to other indices can improve the accuracy of the results when detecting fire scars and fire severity [14,29,62,63]. The added value of integrating TassCap derived indices Brightness, Greenness, and Wetness lies in their integration of a larger range of spectral information instead of limited band ratios. The band coefficients for surface reflectance of Landsat TM are given by Crist [54].

2.4. Analysis and Metrics

Figure 2 illustrates a conceptual model of the different temporal trajectories of burned and unburned samples. The proposed model has been adapted from previous models developed initially by Key [33] and then modified by Chompuchan & Lin [64]. We assume that the ecological condition of vegetation, which could be measured by a given spectral index, follows a regular pattern through time. In addition to natural seasonal fluctuations, the trajectory is also subject to a certain degree of interannual variability due to multiple factors such as different interannual climate/weather conditions, differences in phenological timing, or residual radiometric inaccuracies. Our analysis follows a temporal resolution scheme of 1 year but cloud-free Landsat scenes are rarely available at the same date for several consecutive years. For this reason, instead of analyzing the single trajectory of burned samples, we decided to compare it with the trajectory of reference unburned pixels of same vegetation type. Considering the spectral fluctuations of vegetation not affected by fire gives a more realistic picture on the ecological condition of the disturbed vegetation and presents some interesting advantages. For example, a simple comparison of solely burned samples at two different dates, such as between pre-fire and post-fire, doesn’t account for potential variations of vegetation phenology and unusual weather events which may affect the overall vegetation. Consequently, an undefined bias in the quantification of change is introduced. According to specific purposes, this bias may need to be quantified, especially when working with long-term time series [17,65]. Another advantage of this approach lies on the possibility to analyze differences in the pre-fire condition of burned and unburned vegetation, illustrated in Figure 2 as the distance $t_{burn}^0$ − $t_{ref}^0$. However, this is not the focus of the present study.

According to the conceptual model (Figure 2), the physical effect of burning on vegetation will result in an abrupt deviation from the regular trajectory, reaching the maximum distance short time after the fire or, if a delayed vegetation mortality effect exists, several months or years after the fire [64,65]. This distance corresponds to the magnitude of maximum severity and is quantified as $t_{burn}^m$ − $t_{ref}^m$. At any time after the point of maximum severity, residual severity can be calculated using the distance $t_{burn}^a$ − $t_{ref}^a$.

In this paper, the capabilities of the previously described SI to detect burned areas at several dates after the fire events were evaluated and compared using the following three approaches:

• Visual analysis of temporal plots of the SI

Using Figure 2 as a model, temporal trajectories of each SI were plotted to visually analyze the spectral response of vegetation after a fire event, compared to unburned reference samples. This kind of graphic representation allows for a quick preliminary assessment of SI behavior and fire disturbance patterns. Specifically, the trajectories of burned and unburned samples show the effect of fire and the natural temporal fluctuation of unburned vegetation. Maximum severity and the existence of a delayed mortality effect can be identified. Moreover, recovery and residual severity dynamics can be followed along the timeline.

• Analysis of burned-unburned separability capabilities of SI, the M-statistic

The ability of each SI to separate burned from unburned pixels can be compared using statistical metrics that quantify the spectral distance between two distributions. The M-statistic is a measure of class separability defined by Kaufman and Remer [66] as the difference between the means (μ) of two classes normalized by the sum of their standard deviations (σ). It can be interpreted as a signal-to noise-ratio, where the distance of the means represents the signal and the sum of the standard deviations represents the noise. This metric has been previously used in several studies related to burned area discrimination [27,60,67,68]. The M-statistic is expressed as:

$$M=\frac{\left({\mu }_1-{\mu }_2\right)}{\left({\sigma }_1+{\sigma }_2\right)}.$$

A value of M > 1.0 indicates good separation between the burned and the unburned classes, while M < 1 indicates poorer separation, given by closer means and/or larger standard deviations. Figure 3 shows a graphic representation of a burned area and the M value a pre-fire and several post-fire dates.

• SI enhancement of residual severity

Recovery and residual severity can be quantified by dividing their magnitudes at a given post-fire date with the magnitude of the disturbance at the time of maximum severity. This concept has been used by Lin et al. [69] to assess vegetation recovery after a landslide. The authors defined the Vegetation Recovery Rate which uses NDVI as input. The same structure was adopted by Chompuchan & Lin [64] to assess recovery after fire disturbances, but NDVI was replaced with NBR. We adapted this methodology to compare the capabilities of each selected SI to enhance residual severity at several post-fire dates. Our adjustments consist in the quantification of the residual severity instead of the recovery (see Figure 2), and the use of SI values of unburned plots as a reference for healthy vegetation instead of pre-fire values of burned samples. We name this metric Residual Severity Ratio (RSR), and using the terminology presented in Figure 2, we define it as the ratio between Residual Severity and Maximum Severity. RSR is expressed with the following equation:

$$RS{R}_{SI}=\frac{SI\left(t_{burn}^a\right)-SI\left(t_{ref}^a\right)}{SI\left(t_{burn}^m\right)-SI\left(t_{ref}^m\right)}.$$

RSR values close to 1, indicate low recovery from the stage of maximum severity and may indicate that an SI shows good capabilities to enhance residual severity. However, when interpreting this metric, particular attention should be given to the pattern of the SI trajectories of the burned and the unburned classes displayed on the graphic plots and to the separation capabilities of the SI at the time of maximum severity.

Following the three approaches mentioned above, the analysis was initially performed on the overall dataset, using the mean of all samples of all fires. To investigate possible differences in SI responses between dominant ecosystems in northwest Yunnan, sample points were grouped according to their vegetation class, and the same analysis was performed.

3. Results and Discussion

3.1. SI General Response to Fire in NWY

A comprehensive representation of the temporal trajectories of the 11 SI is illustrated in Figure 4. Most of the SI present similar response patterns following fire, with a pronounced deviation right after the fire (point of maximum severity), a relatively fast recovery during the first year after the fire, and a slow recovery during the following years. Visually, all SI except the TassCap Wetness show a well-marked response at the immediate post-fire assessment which confirms their wide employment for burn severity evaluation. However, a less sharp recovery slope following the maximum severity point suggests that NDMI and NBR could better enhance residual severity than the other SI several years after the fire, while MIRBI appears to be the least performant. The observation of the graphic plots also reveals singular behaviors of Brightness and Wetness. The former has a similar response than the other indices right after fire: the effect of burning on the biomass and the deposit of ash and charcoal significantly decrease the values of Brightness. Though, this effect seems to last very shortly. In fact, one year after the fire, the index’s values increase drastically, becoming higher than those of unburned vegetation (reference and burned lines crossed). This is probably due to the effect of the rainy season washing away dark deposits of ash and charcoal and exposing bare soil, while new green matter (leaves, grass, etc.) hasn’t had the time to recover properly. The Wetness index follows a different pattern. The immediate effect of burn seems to have only a slight decreasing effect on Wetness values, but the magnitude of this effect grows in the following years, showing a clear difference from the unburned vegetation. This suggests that the long-term effects of fire on vegetation include a decrease of moisture content in the affected site in favor of more arid conditions, which, together with post-fire climate conditions, may affect vegetation regrowth dynamics several years after the fire [70].

The results of separability and residual severity assessments are shown in

Table 3. Separability index (M-statistic) and Residual Severity Ratio (RSR) for the 11 SI. Highlighted cells indicate the three best performers at the given post-fire time in the following rank/color order: 1st, 2nd, 3rd. Pre-fire M-statistic values are shown as reference. MS = maximum severity from which RSR is calculated. Negative values of RSR indicate that the trajectory of the burned class crossed the trajectory of the unburned control class.

	M-Statistic-All Fire Plots (Mean)							RSR-All Fire Plots (Mean)
Spectral Index	Pre	Post	1y	2y	3y	4y	5y	Post	1y	2y	3y	4y	5y
BAI	0.08	1.12	0.28	0.28	0.16	0.17	0.15	MS	0.21	0.17	0.10	0.09	0.08
GEMI	0.12	1.08	0.35	0.30	0.20	0.18	0.16	MS	0.36	0.31	0.21	0.19	0.16
NBR	0.07	0.88	0.60	0.55	0.53	0.48	0.39	MS	0.59	0.52	0.49	0.48	0.38
NDVI	0.05	0.82	0.49	0.43	0.37	0.34	0.28	MS	0.60	0.53	0.46	0.46	0.36
BRI	0.03	0.39	0.13	0.14	0.16	0.16	0.14	MS	−0.39	−0.37	−0.50	−0.49	−0.39
GRE	0.12	1.24	0.36	0.34	0.23	0.23	0.20	MS	0.38	0.33	0.24	0.22	0.18
WET	0.05	0.11	0.31	0.36	0.37	0.36	0.29	0.32	0.93	0.95	MS	0.98	0.78
BAIML	0.13	0.74	0.50	0.47	0.40	0.43	0.32	MS	0.31	0.20	0.15	0.16	0.11
BAIMs	0.09	0.75	0.41	0.34	0.27	0.31	0.23	MS	0.36	0.23	0.17	0.18	0.14
MIRBI	0.12	1.29	0.20	0.07	0.06	0.07	0.03	MS	0.20	0.06	−0.06	−0.06	−0.03
NDMI	0.04	0.64	0.55	0.54	0.54	0.51	0.38	MS	0.77	0.73	0.72	0.71	0.55

. According to the ideal value of M > 1, the separability values are in general not very optimal. At the first post-fire assessment, only MIRBI, Greenness, BAI, and GEMI obtained M values higher than 1. However, if we compare post-fire M values with the pre-fire ones, NBR and NDVI still show a significant difference (pre-fire M: 0.07 and 0.05, post-fire M: 0.88 and 0.82, respectively). Brightness and Wetness appear to be the less suitable for burned area mapping shortly after fire. These results don’t necessarily mean that the indices cannot detect burned areas, but that relying on a static threshold to separate burned from unburned patches may be subject to a higher rate of error. It must be noted that these results are based on the average values of all burn plots which, taken individually, have different time lags between pre-fire and immediate post-fire images and include different vegetation types (see Table 1). The influence of time lag and vegetation type on the quality of post-fire assessments by the SI has been pointed out in numerous research results [71,72,73]. Although several comparative studies mentioned in the introduction recognize NBR as the best SI for immediate post-fire assessment, our results found that MIRBI has the greatest separability capabilities, as confirmed by Lu et al. [25]. However, starting from 1 year after the fire, separability radically decreases for the four best performers mentioned before, while NBR, NDMI, NDVI, BAIML, and BAIMs proved to be more resistant. Among them, NBR, NDMI, and BAIML are the three best performers in all five post-fire years. In terms of RSR and as shown in the graphic plots, NDMI is the most suitable SI to enhance the residual severity several years after fire, even if its separability potential decreases. NDVI shows good RSR at 1 and 2 years after fire, while Brightness takes over at 3, 4, and 5 years after fire. Nevertheless, the separability capability of Brightness at these time frames remains poor. The case of Wetness needs further clarification. Maximum severity and best separability of Wetness, albeit rather poor, is reached at the 3rd post-fire year. The trajectory of the burned samples, as showed in the graphic plot, indicates that the biggest shift from unburned samples appears at the 1 year post-fire assessment. Afterwards, RSR remains relatively stable, which explains the resulting high values. The meaning of this pattern is that the perturbations in Wetness caused by fire persists for a few years before the recovery phase begins. The graphic plots and RSR also indicate that, except for the Wetness index whose trajectory reveals persisting moisture stress, no delayed mortality effect seem to exist regarding fire disturbance in northwest Yunnan. Maximum severity appears always at the immediate post-fire assessment (MS value in Table 3).

3.2. SI Response by Vegetation Type

The four main vegetation types of northwest Yunnan react differently to fire and show dissimilar recovery patterns. These patterns can be observed in the SI trajectories, and can be quantified by comparing the values of M-statistic and RSR between vegetation types. The temporal trajectories of burned and unburned samples for the four major vegetation classes in NWY are represented in Figure 5. Among them, we can notice that needleleaved and broadleaved forests tend to recover slower than shrublands and grasslands. This pattern is found in all the SI. Grasslands samples are characterized by burned and unburned trajectories that quickly converge after 1 year. Recovery of shrublands seems to be a little slower than grasslands, with some residual severity visible at 1 year post-fire in NDMI, NBR, NDVI, Brightness and Wetness. A particular behavior is observed in the shrublands graphic plots. At the 1 year post-fire assessment, the trajectory of unburned pixels shows an important shift similar but not as severe as the burned samples (GEMI, NBR, NDVI, Greenness, Wetness, NDMI). A valid explanation lies in the dates of burn of the shrub burn plots. As shown in Table 1, two of the three plots were burned during the dry season between 2005 and 2006, which is described as one of the most severe droughts in Yunnan and other regions of China in the last 50 years [74]. This extreme climate conditions had a significant impact on vegetation and soils, such as pronounced decreases in moisture content and primary productivity, which are visible in the SI. The trajectories show that the unburned vegetation remains under stress for several years after the drought event, while the burned trajectory joins the unburned one after a few years. This suggests that the shrubs have recovered quickly from the fire but they are still under stress from the general drought condition, showing similar absolute values of SI. Another isolated behavior is found in BAIMs which is the only SI to detect a delayed mortality effect at the 1 year post-fire assessment in needleleaved vegetation. The peculiarity of BAIMs is its combined use of the near infrared and the shorter shortwave infrared bands of the Landsat TM sensor, an approach also adopted by the NDMI. The latter, although not showing a delayed mortality effect in needle leaved forest, tends to remain stable between the post-fire and the 1 year assessment compared to the other SI. Such a pattern may reveal higher sensitivity of BAIMs and NDMI for the detection of delayed mortality. Finally, the Wetness index draws an unusual trajectory in grassland vegetation. Instead of decreasing constantly after the fire like in other vegetation classes, it increases slightly at the immediate post-fire assessment and rapidly rejoins the same conditions as the unburned grass. A possible explanation for this distinct behavior could be that fire boosts recovery of grassland by eliminating dead plant material, competing weeds and small shrubs. Furthermore, with adequate soil moisture conditions, the darkening effect of fire increases soil surface temperature, stimulating earlier growth [75]. In support to this hypothesis, the values of Brightness in grasslands fall drastically in this time frame compared to other vegetation types.

More precise information is given by the M-statistic and RSR results shown in

Table 4. Separability index (M-statistic) and Residual Severity Ratio (RSR) for the 11 SI in the four major vegetation types of northwest Yunnan. Highlighted cells indicate the three best performers at the given post-fire time in the following rank/color order: 1st, 2nd, 3rd. Pre-fire M-statistic values are shown as reference. MS = maximum severity from which RSR is calculated. Negative values of RSR indicate that the trajectory of the burned class crossed the trajectory of the unburned control class.

	M-Statistic						Residual Burn Ratio
	Coniferous Needleleaved Forest
Spectral Index	Post	1y	2y	3y	4y	5y	Post	1y	2y	3y	4y	5y
BAI	1.13	0.70	0.55	0.38	0.43	0.33	MS	0.65	0.41	0.32	0.29	0.24
GEMI	1.11	0.56	0.42	0.30	0.31	0.24	MS	0.57	0.41	0.33	0.31	0.24
NBR	1.32	1.12	1.13	1.05	0.97	0.81	MS	0.79	0.67	0.63	0.67	0.53
NDVI	1.36	0.82	0.78	0.88	0.87	0.61	MS	0.72	0.64	0.58	0.63	0.43
BRI	0.26	0.16	0.29	0.30	0.36	0.25	−0.69	0.44	0.81	0.91	MS	0.70
GRE	1.27	0.63	0.54	0.36	0.41	0.34	MS	0.63	0.50	0.39	0.39	0.30
WET	0.48	0.71	0.70	0.74	0.75	0.70	0.58	0.94	0.92	0.95	MS	0.82
BAIML	0.75	0.85	0.98	0.77	0.85	0.70	MS	0.78	0.47	0.43	0.46	0.35
BAIMs	0.75	0.82	0.87	0.66	0.74	0.60	0.87	MS	0.63	0.59	0.59	0.51
MIRBI	1.31	0.39	0.26	0.06	0.09	0.04	MS	0.35	0.22	0.06	0.08	0.03
NDMI	1.09	1.10	0.98	0.99	0.92	0.78	MS	0.96	0.81	0.80	0.84	0.68
	Broadleaved Forest
BAI	1.09	0.46	0.53	0.29	0.27	0.26	MS	0.32	0.31	0.17	0.15	0.15
GEMI	1.13	0.59	0.60	0.39	0.34	0.35	MS	0.52	0.63	0.45	0.39	0.37
NBR	0.81	0.81	0.76	0.88	0.69	0.60	MS	0.81	0.81	0.81	0.74	0.59
NDVI	1.13	0.77	0.81	0.78	0.74	0.71	MS	0.75	0.76	0.64	0.61	0.57
BRI	0.30	0.21	0.09	0.17	0.20	0.17	MS	−0.69	−0.32	−0.72	−0.79	−0.63
GRE	1.27	0.68	0.71	0.47	0.43	0.43	MS	0.58	0.66	0.49	0.43	0.41
WET	0.21	0.39	0.44	0.50	0.45	0.38	0.41	0.84	0.84	MS	0.94	0.74
BAIML	0.75	0.73	0.72	0.77	0.73	0.59	MS	0.51	0.35	0.25	0.23	0.15
BAIMs	0.77	0.64	0.56	0.50	0.52	0.34	MS	0.54	0.40	0.25	0.23	0.15
MIRBI	1.07	0.41	0.34	0.09	0.03	0.06	MS	0.37	0.30	0.10	0.03	0.05
NDMI	0.64	0.73	0.77	0.81	0.69	0.55	0.94	0.93	0.99	MS	0.94	0.74
	Shrubland
BAI	1.02	0.00	0.00	0.09	0.11	0.11	MS	0.00	0.00	−0.04	−0.05	−0.05
GEMI	1.00	0.23	0.11	0.01	0.05	0.07	MS	0.23	0.10	0.01	−0.05	−0.06
NBR	1.02	0.49	0.38	0.29	0.22	0.15	MS	0.48	0.34	0.26	0.20	0.13
NDVI	0.71	0.32	0.22	0.14	0.05	0.03	MS	0.45	0.32	0.20	0.07	0.05
BRI	0.44	0.20	0.20	0.20	0.17	0.17	MS	−0.47	−0.41	−0.46	−0.37	−0.34
GRE	1.15	0.25	0.12	0.02	0.05	0.05	MS	0.25	0.11	0.02	−0.04	−0.05
WET	0.08	0.34	0.42	0.32	0.30	0.23	0.24	MS	0.98	0.84	0.82	0.57
BAIML	0.77	0.32	0.24	0.10	0.09	0.01	MS	0.14	0.08	0.02	0.02	0.00
BAIMs	0.78	0.17	0.03	0.05	0.06	0.09	MS	0.10	0.01	−0.02	−0.02	−0.04
MIRBI	1.42	0.21	0.11	0.24	0.30	0.18	MS	0.16	−0.06	−0.15	−0.18	−0.09
NDMI	0.69	0.43	0.42	0.34	0.29	0.17	MS	0.65	0.58	0.46	0.42	0.24
	Grassland
BAI	1.90	0.05	0.11	0.08	0.15	0.17	MS	0.02	0.05	0.03	0.05	0.06
GEMI	1.54	0.13	0.07	0.04	0.04	0.07	MS	−0.10	−0.06	−0.03	0.03	0.04
NBR	1.12	0.37	0.32	0.21	0.07	0.04	MS	−0.28	−0.23	−0.14	−0.07	−0.03
NDVI	0.72	0.12	0.13	0.06	0.00	0.05	MS	−0.17	−0.19	−0.08	0.00	0.05
BRI	1.18	0.03	0.07	0.03	0.09	0.05	MS	0.04	0.08	0.04	0.09	0.04
GRE	1.69	0.22	0.15	0.09	0.05	0.01	MS	−0.21	−0.12	−0.07	0.04	0.01
WET	0.47	0.16	0.10	0.04	0.06	0.04	MS	0.43	0.24	0.10	0.14	0.09
BAIML	0.87	0.17	0.12	0.09	0.05	0.01	MS	−0.05	−0.04	−0.02	0.01	0.00
BAIMs	0.91	0.11	0.03	0.05	0.06	0.07	MS	−0.05	−0.01	−0.02	0.02	0.03
MIRBI	1.94	0.02	0.15	0.20	0.11	0.01	MS	−0.01	−0.11	−0.14	−0.07	−0.01
NDMI	0.57	0.40	0.19	0.04	0.01	0.01	MS	−0.63	−0.29	−0.05	0.02	−0.01

. At the immediate post-fire assessment, MIRBI, GEMI, Greenness and BAI have very good separability (M > 1) in all vegetation types. NBR follows with good separability in three vegetation types (M lower than 1 in broadleaved forest only) and NDVI in two vegetation types (M lower than 1 in shrubland and grassland). Brightness proves good separability in grassland (M = 1.18) but very poor separability in other vegetation types. Overall, the best separability is found in needleleaved forest, with all SI having M values above 1 at the immediate post-fire assessment, except for BAIML, BAIMs, Brightness, and Wetness. The last one tends to respond to fire later than other indices, as a result of decreased biomass and gradual loss of moisture content in soils some years after a fire event [62]. The colored highlights in Table 4 clearly tell us that, in all vegetation types and in most of the post-fire dates, NBR and NDMI score among the three best performers in both separability and enhancement of residual severity capabilities. Slightly lower scores are observed in grassland and shrubs ecosystems, where at some dates other SI such as MIRBI and Wetness performed similarly or better. Moreover, NDVI, which is normally used as an indicator of general vegetation condition, obtained separability and RSR scores comparable to the best index in broadleaved and needleleaved forests. RSR values are in general higher in needleleaved and broadleaved forests than it is in shrublands and grasslands, indicating slower recovery and longer persistence of fire scars. Except for NDMI that quantifies RSR at 63% (0.63) with the highest separability (M = 0.40) 1 year after the burn, all other indices perform significantly worse in grass-dominated landscapes, and from 2 years after the burn, scars are almost undistinguishable by any SI. A similar tendency though less pronounced is observable in shrubland ecosystems.

The interpretation of the performance of Brightness, especially regarding RSR values which are often higher than other SI, is prone to confusion. First, it is important to note that its separability abilities are quite low in all land cover types and at any post-fire date, except for a noticeable M = 1.10 in grassland at the immediate post-fire assessment. Second, abrupt shifts of the burn trajectory versus the unburned reference trajectory as indicated by RSR negative values call for particular caution and require further clarification. At the 1 year post-fire assessment, the trajectory of burned samples shifts to values higher than the unburned trajectory in most of the vegetation classes. In the following years, the two trajectories progress almost in parallel with small variations of both reference and burned samples. These variations coupled with the smaller spectral distance at the immediate post-burn assessment (initial maximum severity point) explain the inflated RSR values. In other words, the RSR of Brightness could highlight a persistent disturbed status of the vegetation, however, the M-statistic values being relatively low indicate high overlap between burned and unburned pixels. These results advise against the employment of Brightness for burned area mapping without other supporting tools.

3.3. Summary of Results, Which SI Should We Choose?

Our results in a mountainous region characterized by frequent small fires and a fast recovering vegetation due to abundant summer rains shows that all SI are valid tools to separate burned from unburned vegetation at the immediate post-fire assessment, except for Wetness and Brightness that follow different patterns. Among them, MIRBI, BAI, and Greenness are the best when averaging all vegetation types, while NDVI and NBR get better scores in isolated vegetation types. Other studies obtained similar results, e.g., better discrimination capabilities of BAI [18,22,28], Greenness [24,62], or MIRBI in shrub-savannah ecosystems [25,27,52]. Although NBR is often considered the best SI for burned area mapping short time after fire and therefore widely used for burn severity assessments [7,19,21,23,27,28,29], our results demonstrated that in certain vegetation types other SI would offer a better option. The limitations of NBR for immediate post-fire assessment were already pointed out by Roy et al. [76].

From one year after the fire events, burned areas in shrublands and grasslands are very difficult to discriminate. Despite low discrimination capabilities, NDMI and NBR are the best choice, while Wetness could offer a valuable support in shrublands starting from two years after the fires. In more vegetated areas such as broadleaved and needleleaved forests, fire scars are more persistent and NBR has often shown to perform better than in sparsely vegetated land [20,77]. From one to five years after the fire, the best SI are NBR and NDMI, followed by NDVI in broadleaved vegetation and BAIML in needle leaved vegetation. Brightness and Wetness have lower separability capabilities but can be useful to analyze residual severity and recovery. The other SI analyzed in this study performed consistently worse and proved to be less suitable for long term monitoring of fire disturbance. The spectral response of BAI strongly depends on the persistence of charcoal deposits [18]. Both of its modified versions, which were initially designed for the MODIS sensor, couldn’t reach the same performance as the best SI, as also found by [27], and therefore represent a second choice for long term monitoring of post-fire vegetation recovery. The spectral behavior of GEMI is similar to Greenness but results of metrics are almost always worse. Finally, MIRBI resulted very unsuitable for long term burned area detection.

Our results are in line with recent studies on SI capabilities for mapping vegetation recovery, which found that the SI that integrated the shortwave infrared domain of the electromagnetic spectrum such as NBR, NDMI, and Wetness are more suitable for long term vegetation recovery, while indices relying on visible and near-infrared domains were more suitable for immediate post-fire assessment [29,30,31,32]. Our graphic plots show that this is particularly true in less vegetated ecosystems (shrubs, grass). However, we also found that NDVI is as robust as NBR and NDMI in broadleaved forest several years after fire.

Finally, it is important to note that M-statistic and RSR only quantify separability between burned and unburned samples and residual severity magnitude. These metrics do not quantify the quality of burned area mapping. Further analysis of the mapping capabilities of SI should include change detection and classification methodologies with detailed accuracy assessments quantifying omission and commission errors. Moreover, to improve our knowledge of fire in mountainous regions, additional research on burned area detectability by SI in different slope/aspect and hill shading conditions as well as their capabilities to discriminate fire disturbances from other types of change is highly suggested.

4. Conclusions

The wide-ranging amount of literature on Landsat spectral indices performance for the detection of burned areas suggests that there’s no unique model that over performs the others in any spatio-temporal circumstances. Instead, different vegetation types and different time lags between image acquisition dates greatly influence the spectral indices behavior. In this study, the spectral trajectories of eleven among the most common spectral indices widely employed in burned area extraction have been analyzed in four major vegetation types in a mountainous area of northwest Yunnan, China. As previous literature suggested, two main performance variability patterns can be highlighted. First, a temporal pattern: significant performance differences were observed between immediate post-fire assessment and from one to several years after fire. Second, a vegetation type pattern: a clear difference between more vegetated (needleleaved and broadleaved forest) versus less vegetated landscapes (shrublands and grasslands) exists. Therefore, in order to optimize accuracy of burned area detection algorithms, the best spectral indices should be selected accordingly. In general, when time lag is short (<1 year), MIRBI and the Tasseled Cap Greenness are the most suitable indices, followed by the BAI. When time lag is longer (>1 year), NBR and NDMI are the most suitable to detect residual severity and more realistically track vegetation recovery. However, fire scars in grassland ecosystems are almost undetectable by any spectral index after one growing season. The Tasseled Cap Wetness revealed to be very useful for monitoring long-term residual severity starting from 2 years after a fire.

Delicate mountain ecosystems are particularly sensitive to climate variations which, together with patchy land cover, frequent clouds, and relatively fast vegetation regeneration, make it challenging to systematically extract burned areas and assess impacts and recovery. To better isolate fire disturbance from the overall variation of the vegetation’s ecological condition, we designed a conceptual model that relativizes the spectral response of burned pixels by integrating the spectral response of unburned pixels. The metrics used to quantify class separability (M-statistic) and magnitude of residual severity (RSR), were adapted to this model. The high inter-annual variability of unburned vegetation observed in the resulting graphic plots confirms the validity of this approach. Therefore, relying solely on the pre-fire condition of the burned vegetation to assess recovery, may introduce bias due to other disturbances. Instead, keeping track of general vegetation conditions over time allows for the identification of other environmental factors influencing recovery, such as general drought. Thus, this approach is particularly suitable to assess vegetation vulnerability to fire under changing climate and weather conditions.