Estimating the area burned by agricultural fires from Landsat 8 Data using the Vegetation Difference Index and Burn Scar Index

Additional keywords: remote sensing of environment.

Farmers worldwide often remove excess crop residue from fields by burning farmland. The economic benefits of crop residue utilisation are not high, and burning is a quick and cheaper management method. Fire can also remove weeds, pests and diseases, enabling farmers to prepare the next crop (McCarty et al. 2009). In the early ripening season (mostly early June) in winter wheat planting areas such as the Haihe River Basin in northern China, it is easy for farmers to lose control of fires owing to dry weather, wind and other natural factors. This can result in the loss of unharvested crops. Frequent farmland fires have serious consequences, such as fire-induced trace gas emissions (Andreae and Crutzen 1997; Andreae and Merlet 2001; Yang et al. 2008; Hao and Larkin 2014) and grain loss. Farmland fires that spread to neighbouring grasslands and forests are also one of the main causes of wildland fires (Houghton et al. 2000). Some of these fires are caused by stubble and straw burning, whereas others involve the burning of mature crops and are made worse by dry conditions (Maingi and Henry 2007; Vadrevu and Lasko 2015). Because of the randomness and rapid spread of fires on farmland such as the winter wheat planting area in the Haihe River Basin, fire managers find it difficult to determine the location of fires and measure the extent of burned area. Retrieval of accurate burned area data is essential for determining the source of wildland fire, and modelling air pollution and grain loss (França et al. 2014).

Most previous studies estimating burned area through remote sensing have classified images using principal component analysis and vegetation indices (Carlson and Ripley 1997; Chuvieco et al. 2002; Domenikiotis et al. 2002; Hudak and Brockett 2004; Mitri and Gitas 2004; Kučera et al. 2005; Loboda et al. 2007; Maingi and Henry 2007; Smith et al. 2007; Chuvieco et al. 2008; Palandjian et al. 2009; Stroppiana et al. 2009; Boschetti et al. 2010; Bastarrika et al. 2011; Parker et al. 2015). Many studies have relied on the Normalized Difference Vegetation Index (NDVI), including modified versions to reduce the sensitivity of the index to different atmospheric and soil conditions (Chuvieco et al. 2002; Domenikiotis et al. 2002; Kučera et al. 2005; Stroppiana et al. 2009; Veraverbeke et al. 2011a). Some of these modified NDVI indices, such as the Enhanced Vegetation Index (EVI), Soil Adjusted Vegetation Index (SAVI) and Modified Soil-Adjusted Vegetation Index (MSAVI), have been used effectively for burned area estimation in a range of habitats such as forests and grasslands (Huete 1988; Pinty and Verstraete 1992; Qi et al. 1994; Huete et al. 2002). Some researchers have also proposed other spectral indices, such as the Burned Area Index (BAI) and Normalized Burned Ratio (NBR) (Chuvieco et al. 2002; French et al. 2008; Veraverbeke et al. 2010b; Araújo and Ferreira 2015), which are particularly sensitive to the spectral features of a burned area. However, the existence of different land types can easily lead to spectral confusion for burned areas and spectrally similar areas, such as water bodies, roads and buildings, which results in uncertainties in burn scar models, largely due to the difficulty in determining threshold values for these models (Lasaponara 2006; Stroppiana et al. 2009; Boschetti et al. 2010; Veraverbeke et al. 2011a, 2011b; Boschetti et al. 2015). To deal with this type of spectral complexity (spectral mixing), various models and methods have been used to highlight fire-induced changes in pre- and post-fire imagery; however, this can potentially create additional constraints in relation to image-to-image normalisation (Conghe Song and Woodcock 2003; Verbyla et al. 2008; Veraverbeke et al. 2010a, 2010b, 2011a, 2011b).

A few remote sensing methods have been used to extract burn scar information, with a focus mainly on forests, grasslands and other non-farmland land types. Only a few models and methods have addressed the issue of identifying burned areas within farmland. Cultivated land is often distributed unevenly and characterised by a mixture of different land-use types, such as rural settlements, mines, water conservation and irrigation facilities, rivers, lakes and asphalt roads. As a result, different land types often share the same spectra, which creates interference when burn scar information is extracted from remote sensing data. In addition, because of differences in crop types and planting times, the distribution of crop stubble, as well as the presence of mature and non-mature crops, is likely to cause interference, leading to uncertainty in estimations of burned crop area.

The main purpose of the present study was to develop a method for effectively extracting information on the burned crop area caused by agricultural fires from remotely sensed data. Considerable spectral confusion exists in relation to distinguishing burned areas from areas with similar surface features in monotemporal imagery (Lasaponara 2006; Stroppiana et al. 2009; Veraverbeke et al. 2011b). The following issues are of considerable importance regarding the extraction of precise information on burned areas: (1) the method should allow rapid identification of the spatial distribution of mature crops in such a way that researchers can determine the possible range of a burned area, and can reduce the complexity of background information; (2) furthermore, the method can clearly highlight burned areas and integrate different types of background information, for example on mature crops and stubble.

To test the accuracy of the proposed method in determining the locations of burned areas, we selected the winter wheat planting area in the Haihe River Basin. The growth curve features extrapolated from a time series of Moderate Resolution Imaging Spectroradiometer (MODIS) data; based on this, two NDVI images extracted from Landsat 8 data were selected for this area. Burned areas with different spectral features, based on the band features of Landsat 8 data, were analysed and a burn scar index (BSI) was then constructed for extracting information on the burned area. Model application and validation were conducted for the winter wheat planting area of the Haihe River Basin. The method greatly improved the accuracy of satellite image-based burn scar estimates of burnt areas resulting from farmland fires. It can also be used to determine the local and regional contributions of particulate and trace gas emissions, which affect both air quality and public health (McCarty et al. 2009; Li et al. 2014; Chen et al. 2017).

The Haihe River Basin, China’s major grain production base, is located in central China, and covers a total area of 318 200 km². The study area is located west of the Haihe River Basin and the main grain crops include winter wheat and maize, as shown in Fig. 1. The mean annual temperature in the study area ranges from 12.7 to 13.7°C, whereas the mean annual precipitation is about 600 mm, based on data from the Anyang Meteorological Bureau. Only sparse precipitation occurs during spring and the air is generally dry; therefore, winter wheat normally requires irrigation. After the winter wheat reaches maturity or is harvested, prescribed fires or wildfires in farmland areas may cause serious air pollution and a loss of grain harvesting opportunities. On 9 June 2015, an agricultural wildfire disaster resulted in the loss of many human lives and large amounts of grain in the study area.

Fig. 1.  Location of the Haihe River Basin, China. The study area is located to the south of the Haihe River Basin and is shown by the red box.

Soil, winter wheat and burned area spectral data were collected through field experiments during winter wheat growth and maturity periods. Winter wheat spectral data were collected from 4 April to 3 June in both 2001 and 2002, and stubble and burn spectral data were collected from 9 to 14 June 2015.

An ASD FieldSpec3 spectrometer (Analytical Spectral Devices, Boulder, CO, USA) and a PSR spectrometer (Spectral Evolution Co., Lawrence, MA, USA) were used for field data collection. These devices have a bandwidth between 350 and 2500 nm, a viewing angle of 25° and a height of 20 cm to the measured sample. Measurements were conducted on sunny days with favourable visibility according to the criterion of 3 days in a row without precipitation prior to spectral measurement. The areas surrounding the measurement points were broad, with no large obstructions, and the measurements were conducted at the local time of 1000–1400 hours. Each sample was measured 10 times, and the average reflectance of each sample was then calculated. The total collected sample number was 212. Portions of the spectrum near 1900 nm were removed because of noise.

A time series of MODIS 09A data products was selected for crop growth curve analysis. Landsat 8 reflectance data were used for the vegetation difference index (VDI) and BSI models, and Gaofen-1 satellite (GF-1) (Jia et al. 2016) data were used for model validation (Table 1). A topographic map with a scale of 1 : 100 000, and some latitude and longitude coordinates, were derived from a Google Earth map of the study area, enabling geometric rectification to be conducted.

Table 1.  Sources for the remote sensing data employed in the present study
Abbreviations: MODIS, Moderate Resolution Imaging Spectroradiometer; GF-1, Gaofen-1 satellite; NIR, near-infrared; SWIR, short-wave infrared

ENVI software (ver. 4.8; ITT Visual Solutions, Boulder, CO, USA) was used for spatial geometric precision correction. The correction was accurate to within half a pixel, which enabled transformation of the data into a universal transverse Mercator (UTM) projection.

All analyses of the burned area estimates for the study region were conducted using ENVI software (ver. 4.8). Fig. 2 is a flowchart of the steps implemented to satisfy the study objectives. The methodology can be divided into four parts: (i) deriving a crop phenology curve from a time series of MODIS NDVI data, which is used to select the date of the TM image. The reference TM image was obtained during the vigorous growth period of the winter wheat; the monitored TM image was obtained just after the fire; (ii) development of a burned area extraction method based on VDI and BSI; (iii) cross-comparison with the results from a previous model (Table 2); and (iv) burned area simulation using VDI and BSI.

Fig. 2.  Overview of the methodology implemented in this study. Abbreviations: BSI, Burn Scar Index; VDI, Vegetation Difference Index; MODIS, Moderate Resolution Imaging Spectroradiometer; NDVI_ref, Normalized Difference Vegetation Index of the reference remote sensing image; NDVI_moni, Normalized Difference Vegetation Index of the monitored remote sensing image.

Table 2.  Spectral indices used in this study
Abbreviations: R, red; NIR, near-infrared; SWIR, short-wave infrared; LSWIR, longer short-wave infrared; L, adjusted parameter

Ground spectral measurements of the typical underlying surfaces in the study area were conducted using a portable ground object spectrometer (Model PSR-3500m Spectral Evolution Co.), with the spectrum ranging from 400 to 2500 nm (Figs 3 and 4). Moisture absorption bands at 1360–1420 and 1780–1986 nm were removed from the curves.

Fig. 3.  Photographs of mature winter wheat (a), stubble (b) and burned areas (c) within the study area.

Fig. 4.  Spectrum curves of the typical underlying surfaces in the study area overlapped with Thematic Mapper (TM) bands. The continuous spectrum curves were measured using a spectrometer, and the short horizontal lines are the second to seventh TM bands of the corresponding underlying surfaces.

As shown in Fig. 4, the winter wheat reflectance spectrum for the filling stage peaked at 530 and 1610 nm, with absorption troughs at 680 and 1420 nm and a high reflection region between 760 and 1300 nm. The trend in the variation of the spectrum for mature winter wheat was similar to that of green winter wheat, although the reflection peaks and absorption valleys were more obvious in the former. The reflection peaks and absorption troughs of mature winter wheat in the infrared region (760–1100 nm) were lower than those at the filling stage, whereas in the other bands, they were higher.

The trend in the variation of the burned area spectrum was similar to that of the stubble and soil spectra throughout the entire spectral range. The reflectance of the burned area was generally weak in the 400–1800-nm band range, while it was higher than that of mature winter wheat and the filling stage at 1800–2500 nm. The curve in the 1800–2500-nm region was similar to those of soil and stubble. This result was similar to that reported in previous studies (Lasaponara 2006; Smith et al. 2007; Chuvieco et al. 2008). The purpose of the ground spectral measurements was to obtain the features of different underlying surfaces and then use them to select bands from Thematic Mapper (TM) images. For certain features, the second to seventh TM bands overlapped the spectrum curves measured by the spectrometer, as shown by the short horizontal lines in Fig. 4.

For the VDI parameters and , in the subtrahend represents the NDVI value during the growth period, and represents the NDVI value during the mature period. An analysis of the NDVI curves for different growth stages of the typical underlying surfaces was conducted in the Haihe River Basin (Fig. 5). The NDVI curve from the time series of the MODIS data was used to analyse trends in the NDVI with respect to differences between crop area and other typical underlying surfaces, such as forests and grasslands, and cities and towns. The time series of NDVI values for winter wheat began to increase from Julian day 40 and peaked between Julian days 90 and 120, when the curve began to decline, reaching its lowest value at approximately Julian day 170 (Fig. 5).

Fig. 5.  Changes in wheat phenology over 1 year obtained from MODIS time series data.

In comparison with the NDVI curves of forests and grasslands, and of cities and towns, the NDVI curves of the winter wheat were higher before Julian day 140, but declined before Julian day 170. From the above analysis, it is clear that the NDVI curve of winter wheat changed over time, and was also different from that of forests and grasslands as well as cities and towns.

Based on the above time series of NDVI curve analysis, a VDI was developed, and specific models and parameters were configured as below:

where the NDVI is given by . R_Nir and R_Red are the fifth (near-infrared (NIR)) and fourth (red) bands of the Landsat 8 sensor respectively. in the equation represents the NDVI value during the growth period and represents the NDVI value during the mature period.

Previous observations and analyses (Veraverbeke et al. 2011a) demonstrated a low discriminatory power for the visible spectral region, and indicated that the highest sensitivity is in the short-wave infrared (SWIR) spectral region.

The specific parameters of the BSI are described below. Background features such as soil, stubble and mature crop, can be defined as:

where , R_Red, R_Nir, and are Landsat 8 bands 7 (SWIR 2), 4 (red), 5 (NIR) and 3 (green) respectively. The equation is proposed based on spectral analysis of the relationship between the burned area and the background information. The value of a burned area is similar to those of soil and stubble, and is higher than that of a mature crop. However, the value of a burned area is lower than those of soil, stubble and mature winter wheat. In the denominator, m is an adjustment factor. An appropriate value of m for can highlight burned areas, reduce the influence of background information, and reduce differences between the different types of background information (thereby reducing uncertainty in the estimation of burned areas). A value of 4 for m was considered appropriate for the present study according to the simulation analysis (Tables 3 and 4).

Table 3.  Comparison of Burn Scar Index (BSI) values of typical underlying surfaces and C() values for different m values(the bold number is the biggest C() value of each m value)
See text for definition of terms

Table 4.  Comparison of the Burn Scar Index (BSI) values of typical underlying features and ND() (normalized C()) values for different values of m
See text for definition of terms

C denotes comparison between various types of background information and the information valuable for the BSI.

where is the mean BSI of the burned stubble area retrieved from a Landsat 8 image, is the mean BSI of the burned mature crop area in the same image, is the mean BSI of the unburned stubble area, and is the mean BSI of the unburned mature crop area. C() is the difference between the two compared parameters; the larger the value of C(), the greater the difference between the two compared parameters. The proportions of these different land-cover types within burned and unburned areas were validated by field sampling and high-resolution satellite images (i.e. GF-1 data). Based on the classification results for the Landsat 8 image, the locations of the samples were obtained and the mean value of each cover type was calculated.

To facilitate the analysis, the range of C() was normalized to (0,1) (Wang et al. 2015). is the maximum value of , , and , calculated as follows:

where , , , and are the normalized values of , , , and respectively.

By retrieving the values of , , and for the study area from a Landsat 8 image, a set of m values from 0 to 10 was obtained in increments of 0.5. The m values of 0, 2, 4 and 6 were selected and used in example calculations. Then, the BSI and C() values were calculated, as shown in Table 3. The maximum values of , , , and were calculated for different values of m. , , , and were calculated as shown in Table 4. Fig. 6 shows the four ND() curves with m values ranging from 0 to 10, in increments of 0.5. An ND() value of 1 indicates that the two parameters being compared are totally distinct, whereas the closer the value is to 0, the less distinct are the parameters. In Fig. 6, in a comparison of the different m values for the BSI model, each of the four ND() curves reaches a saturation value when m is close to 4. A BSI model in which m = 4 has a high compression capability and can provide detailed background information.

Fig. 6.  Comparison of the four types of ND() value for the Burn Scar Index (BSI) model according to different values of m.

The VDI model was prepared using two NDVI images extracted from Landsat 8 data (Fig. 7) and 950 groups of BSI values were obtained from field samples and high-spatial-resolution GF-1 satellite data. The 2-m-resolution GF-1 remote sensing fusion data were based on the multispectral and panchromatic bands images by using the Nearest Neighbor Decision (NND) fusion method. The overall precision of the model exceeded 95% with respect to estimation of mature winter wheat and stubble areas.

Fig. 7.  Normalized Difference Vegetation Index (NDVI) values on Julian days 145 (a), and 161 (b) using Landsat 8 reflectance data.

Mature winter wheat data could be extracted with the VDI model by setting a threshold, with the results indicating that an obvious difference between and values allows more effective extraction. Therefore, a large NDVI value as and small NDVI value as should be selected to present a very different trend in the change in forest and grassland and cities and towns, from that in winter wheat cover between Julian days 140 and 170 (Fig. 8).

Fig. 8.  Estimation of the area of mature winter wheat according to the VDI model. (a), (b), (d), and (e) were derived from Landsat 8 bands 6 (short-wave infrared (SWIR) 1), 5 (near-infrared (NIR)) and 4 (red) on Julian days 145 and 161. In Fig. 8, (a) is the burned area, (b) and (e) are stubble, and (c) is mature winter wheat. In Fig. 8b, e, the letter ‘a’ stands for the underlying surface of the city and town, ‘b’ and ‘e’ for stubble, ‘d’ for water body, ‘g’ for water body and ‘h’ for cloud; the blue part in (c) and (f) indicates mature winter wheat or stubble areas. The red square is generated by the image processing software, no meaning here.

The accuracy of the BSI model was verified through an intercomparison of 286 groups of Landsat-based BSI values with objects observed in high-resolution satellite imagery and field surveys. The accuracy assessment was performed within a short period of time after Landsat data acquisition. The validated burned area, roads and crop residue were obtained by field surveys and high-resolution remote sensing images (GF-1). The specific features corresponding to each class were digitised in the high-resolution images. The digitised features representing burned area, crop residue and roads were then used to produce a random stratified point sample for each class. Finally, this random sample of points was used as the ground truth for accuracy assessment of the Landsat-based classification. The false-positive, false-negative, missing alarm and false alarm (Table 5) rates were found to be acceptable according to evaluation of the target detection accuracy (Burke et al. 1988; Lienhart and Maydt 2002; Dumitrescu et al. 2003; Tyre et al. 2003).

Table 5.  Precision analysis of burned areas based on the Burn Scar Index (BSI) model
Abbreviations: SAVI, Soil Adjusted Vegetation Index; CSI, Char Soil Index; BAI, Burned Area Index; NBR, Normalized Burned Ratio; VDI, Vegetation Difference Index; BSI, Burn Scar Index

Compared with other models, such as SAVI, the Char Soil Index (CSI) and BAI, the new method presented herein is clearly more effective for estimating burn scar areas.

Based on the BSI data in Fig. 9, a number of conclusions can be drawn. (1) The BSI model can estimate burned areas, such as in Fig. 9c, but shadows from clouds may result in misidentification (as seen in Fig. 9f) due to spectral differences between shallow cloud, smoke and burned areas. (2) Burned areas differed markedly in shape and size (Fig. 10). Large areas represent severe wildfires in mature winter wheat or stubble areas; within these areas, some of the fires were natural wildfires, whereas others were caused by straw burning. (3) A large number of burned areas were apparent in the images, indicating that agricultural wildfires generally occur during the winter wheat harvesting season, and these fires may cause air pollution and grain loss.

Fig. 9.  Burned area estimates for different underlying surfaces (the letter ‘a’ and the other large red patches in (c) indicate areas burned by wildfires; the letter ‘b’ in (f) denotes cloud shadows, whereas ‘c’ in (i) indicates smoke from straw burning).

Fig. 10.  Information extracted from the study area using the BSI model.

To summarise, the following conclusions can be drawn. (1) Accurate identification of a crop planting area requires highly precise information on burned areas. During the harvest period, wildfires are often caused by prescribed fires that were set to burn crop straw or stubble. These wildfires are a direct cause of air pollution and grain loss, especially in north and north-east China (Shi et al. 2014; Long et al. 2016). The occurrence of wildfires is related to local climatic conditions and land-management practices. Therefore, local fire management personnel need to monitor both prescribed fires and areas affected by wildfires. (2) The VDI model may effectively simplify background information during monitoring of burned areas. Using to differences between the growth curves of crop areas and sources of background information, the VDI model was developed, which can be used to effectively estimate mature crop areas in regions where wildfires may occur. Through the use of the VDI model, certain types of background information, such as the locations of water bodies, grasslands and roads can be filtered out, with the result that the BSI model may be more effective for estimating burned areas. (3) The simulations and experiments conducted during this study showed that the BSI model was moderately reliable in determining crop areas. The parameters used in the BSI model were applied to the winter wheat planting area in the Haihe River Basin. The results showed that the BSI model may be sufficiently precise for monitoring burned areas; the VDI model was also shown to be effective.

However, when a comprehensive analysis was performed using a combination of on-site investigation, simulation analysis and remote sensing mapping, it was evident that the VDI model was more suitable for monitoring early harvest crops. The early harvest crop areas showed a different trend in cover changes to those of forests and grasslands, and cities and towns during the harvest period. When the growth curve trend changes, the VDI model should be adjusted according to crop type. With the BSI model, it was also difficult to differentiate burned areas of mature crops from burned areas of stubble because of the complexity of the soil spectrum in of burned stubble areas.

The authors declare they have no conflicts of interest.