A time series decomposition approach to detect coal fires in parts of the Gondwana coalfields of India from VIIRS data

ABSTRACT This paper proposes an algorithm to detect coal fires using time-series Land Surface Temperature (LST) maps generated from data acquired by a wide-swath Visible Infrared Imaging Radiometer Suite (VIIRS) satellite sensor. VIIRS provides remotely sensed thermal images daily during daytime and nighttime. The proposed algorithm primarily uses the trend component extracted at pixel level from time-series LST maps by using the Seasonal Trend decomposition based on Loess (STL) model. The proposed algorithm is suitable for detecting coal fires at the regional scale. The obtained results were validated with coal fire observations collected during a field survey.


Introduction
India has third ranking in global coal production and second ranking in global coal consumption (BP 2021).Commercial coal mining in India was reportedly started in 1774 in Raniganj Coalfield, situated in West Bengal state (Singh and Yadav 1995).Due to unscientific mining techniques, Indian coal mining companies faced the problem of coal fires.The first coal fire incident occurred in India in 1916 in the Jharia coalfield situated in Jharkhand state (Michalski et al. 1997), and coal fires are still occurring.Most coal fires occur due to the spontaneous combustion of coal (Misra and Singh 1994).When coal gets in contact with oxygen, the heat is liberated.If there is no passage of heat in the atmosphere, it accumulates and raises the temperature at the surface of coal.The process of increase in temperature continues and the coal gets ignited.Coal fires are notoriously recognised for destroying prime coking coal, causing land subsidence, which results in lethal accidents to the people residing in the coalfield area (Chatterjee et al. 2015), making agricultural land infertile, releasing harmful gases in the atmosphere, which subsequently become root cause of health problems to the local residents (Finkelman 2004, Kolker et al. 2009), and creating obstacles for new coal production.The knowledge of underground coal fires enables the government and agencies involved in mining operations to take action on securing the unburnt coal resources and resettlement of the affected population in coal fire affected areas.It has been reported that extinguished coal fires may start burning again within a few weeks (Kuenzer et al. 2013).Therefore, their monitoring at appropriate intervals can help the concerned departments and agencies make timely action plans to tackle potentially dangerous coal fires.Most previous studies have detected coal fires with an interval of more than 1 year, e.g.Martha et al. (2010) detected coal fires in 2003 and 2006 at 3-year interval, Pandey et al. (2017) detected coal fires during 1988-2013 with consecutive five-year intervals, and Singh et al. (2021) detected coal fires for 2006, 2009, and 2015.The present study monitors coal fires at consecutive 1year time interval from wide-swath data.However, with the availability of wide-swath data on a daily basis, the proposed algorithm can be used to monitor coal fires at shorter intervals.
Surface and sub-surface coal fires raise the temperature at the place of occurrence as compared to their surroundings by liberating heat.These locations are termed coal fireinduced thermal anomalies.Spaceborne remote sensing technology has been widely used to detect and monitor coal fires (Syed et al. 2018).It was used by Cracknell and Mansor (1992) for the first time to detect sub-surface coal fires in the Jharia coalfield using Landsat Thematic Mapper (TM) data.Surface coal fire can be detected from short-wave infrared and mid-wave infrared data, whereas subsurface coal fire can be detected from long-wave TIR data (Hecker et al. 2007).Spaceborne thermal sensors ensure data continuity and cover a larger spatial extent required for detecting and monitoring coal fires over a large geographical region.Most thermal anomaly detection methods use manual thresholding, also called density slicing (Prakash et al. 1995, Guha et al. 2008).This approach separates coal fires and background pixels using a manually defined threshold.This technique requires field knowledge of some coal fire locations.A threshold is arbitrarily defined such that the known coal fire pixels fall in one part of the threshold and the background pixels fall in the other part.The threshold is adjusted until the desired results are achieved.The threshold can also be derived by using statistical techniques like the one developed by Kuenzer et al. (2007).In this technique, a thermal image is subjected to a moving window filter.Suppose that the moving window size is M1× N1 and the step size to move the window over the image is X × Y.Then, each pixel in the window is sampled (M1× N1)/(X × Y) times.The window size usually varies from 11 × 11 to 35 × 35, and the step size is kept as 1 × 1 to maximise the number of times a pixel is sampled.The histogram of each subset around the central pixel is analysed.The first local minimum after the global maximum is set as the threshold to separate thermally anomalous and background pixels in the subset.This operation is done more than 1000 times, and the pixels labelled as thermally anomalous in more than 70% of cases are considered coal fire pixels.Chatterjee (2006) determined the threshold temperature for surface coal fire pixels by field-based modelling of the pixelintegrated temperature for surface coal fire mixed pixels in thermal infrared Landsat TM data.Roy et al. (2015) proposed a new statistical approach in which they picked up some random blocks from the study area.A scatter plot between mean and maximum temperature was plotted.A coal fire and a non-coal fire cluster were formed in the scatter plot.The maximum temperature in the non-coal fire cluster was taken as the threshold.
Most of the studies related to coal fire detection have used medium-resolution TIR data from sensors like Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER), Landsat Thematic Mapper (TM), Enhanced Thematic Mapper Plus (ETM+), Thermal Infrared Sensor (TIRS), etc. Medium-resolution TIR data has its limitations in spatial coverage and revisit time.The swath of ASTER and Landsat is 60 km and 185 km, respectively, which is inadequate for regional mapping of coal fires.The revisit time of ASTER and Landsat is 16 days, which is very coarse and results in the availability of a smaller number of observations in the same area over the study period.Since the winter season is very good for detecting coal fires (Roy et al. 2015), the availability of goodquality observation further reduces due to thin cloud cover and bad weather conditions (Mujawdiya et al. 2019).The Visible Infrared Imaging Radiometer Suite (VIIRS) sensor onboard the Suomi National Polar-orbiting Partnership (SNPP) satellite platform overcomes both limitations with a wide swath of 3060 km and full daily coverage of the earth's surface in the daytime and nighttime.Therefore, the availability of good-quality observations substantially increases.
Furthermore, the daytime and nighttime thermal images from VIIRS are available on the same day.This helps minimise false alarms produced by various elements, e.g.water, riverbed sand, built-up land, and anthropogenic activities.VIIRS Land Surface Temperature (LST) data is available through VNP21 products.Coal fire detection using LST maps has an advantage over brightness temperature maps.LST increases the subtle difference between the coal fire and background pixels (Chatterjee et al. 2017).
In this study, we employ Seasonal Trend decomposition based on Loess (STL) model to wide-swath VIIRS satellite images available at high temporal resolution to detect coal fire induced thermal anomalies over the extensive region of Gondwana coal fields of India.The STL model was proposed by Cleveland et al. (1990).STL is one of the statistical approaches aiming to decompose a time series of atmospheric and environmental variables into different components of land surface dynamics.These components include seasonal variations, gradual trends, and abrupt changes.The seasonal component possesses regular and periodic fluctuations in the LST on an annual scale.The trend component reflects gradual changes in the land surface dynamics, and the residual component shows abrupt changes occurring on the land surface.Breaks For Additive Seasonal and Trend (BFAST) (Verbesselt et al. 2010) and Detecting Breakpoints and Estimating Segments in Trend (DBEST) (Jamali et al. 2015) are other models used to decompose remotely sensed time series of atmospheric and environmental variables.These models are specifically designed to detect breakpoints in the remotely sensed time-series datasets.However, we prefer STL over other models since this study does not aim at identifying breakpoints, and in addition, STL is computationally efficient with large datasets (Ben Abbes et al. 2018).STL has been widely used for the decomposition of time series of environmental, geological, and demographic parameters, i.e. normalised difference vegetation index (Lu et al. 2003), groundwater level and weekly discharge (Duy et al. 2021), daily crime data (Yang et al. 2021), and LST (Mujawdiya et al. 2022).
This research aims to detect the spatial distribution of coal fires in parts of the Gondwana coalfields of India over a multi-year period from 2013 to 2019.Detection of coal fires in many of the Gondwana coalfields has rarely been attempted except for the Jharia and Raniganj coalfields.Therefore, this research will be helpful for future researchers who will work on these unexplored coalfields to detect coal fires.

Study area
Gondwana coalfields are situated in Jharkhand, West Bengal, Madhya Pradesh, Chhattisgarh, Odisha, Maharashtra and Telangana states of India.The study area is bounded by geographic coordinates 77.74 N to 87.6 N and 16.68 E to 25.13 E. Talcher, Mand-Raigarh, Raniganj, Ib-River, Godavari, Jharia, North Karanpura, Rajmahal, Singrauli, Korba, East Bokaro, Sohagpur, South Karanpura, and Wardha Valley are some of the major Gondwana coalfields in India in terms of available coal resources and coal production.These coalfields are the source of various types of coal (based on quality), such as prime coking, medium coking, semi coking, superior non-coking, high sulphur, inferior noncoking, and ungraded non-coking.Jharia coalfield is the only source of prime coking coal, which is 0.54% of overall coal reserves, while 97.25% coal of all coal reserves is non-coking type (GSI 2018).Prime coking coal is highly prone to the coal fire phenomenon.Therefore, most coal fire incidents have been reported in the Jharia coalfield.However, incidents of coal fire have also been reported in other coalfields (Sahu and Pal 1998).The study area is shown in Figure 1.Field photographs taken during field surveys in Jharia and Raniganj coalfields are demonstrated in the supplementary file, Figure S1.

Data used
VIIRS onboard the SNPP satellite scans the earth's surface in 22 spectral bands ranging from 412 nm to 12 µm.VIIRS observes the earth's surface in seven bands in the thermal infrared region of the electromagnetic spectrum.Of the seven bands, five have 750 m spatial resolution, and two have 375 m spatial resolution.The main data source of this study is the VIIRS VNP21A2 8-day LST product.The VNP21A2 LST product is produced at 1 km spatial resolution using the Temperature Emissivity Separation (TES) method, which uses three moderate resolution thermal infrared bands at 750 m spatial resolution, namely M14, M15, and M16 (Islam et al. 2017).The wavelength ranges of these thermal bands are given in Table 1.The data used in this research were collected from a web user interface of level-1 and atmosphere archive and distribution system (LAADS) distributed active archive centre for the period January 2013 to December 2019.

Methodology
The trend component of coal fire pixels deviates significantly from the general background trend of surrounding pixels (Mujawdiya et al. 2022).Therefore, this deviation can be quantified and used to segregate coal fire pixels.The LST Pixel Time Series (LPTS) vectors were generated for each pixel in the study area.An LPTS vector represents the time-series LST observations at a specific pixel.Separate LPTS vectors were generated for daytime and nighttime data.The trend components were extracted from each LPTS vector.The working of the STL model used for the decomposition of LPTS vectors is detailed in section 4.2.Yearly deviation maps were generated for both daytime and nighttime by the procedure described in section 4.3.The deviation maps contain information on the difference between the LST values from their surroundings.The deviation maps were generated at different neighbourhood sizes (3 × 3, 5 × 5, 7 × 7, 9 × 9, and 11 × 11) around the pixels.Yearly product deviation maps were generated by taking the product of corresponding yearly daytime and nighttime deviation maps.The presence of water, built-up land, bare soil, river bed sand, industrial heat sources, and so on, may produce false alarms.The product of yearly daytime and nighttime deviation maps significantly reduces false alarms since the radiations from false alarms vary during their diurnal temperature cycles.The procedure employed for detecting coal fire pixels in the yearly product deviation maps is described in section 4.4.The flowchart of the proposed algorithm is shown in Figure 2.

Pre-processing of satellite data
The VNP21A2 data files are provided in Hierarchical Data Format version 5 (HDF5).A single file contains 11 bands.These bands include emissivity maps in wavelength ranges 8.4-8.7 µm, 10.26-11.26µm, and 11.54-12.49µm, daytime and nighttime LST maps, as well as quality control, view angle, and view time bands for LST maps.The required bands were extracted from HDF files and converted into Tag Image File Format (TIFF).The LST maps have sinusoidal projection; therefore, their projection was changed to the WGS84 UTM zone 45 projected coordinate system.The data is

STL model
STL is an additive time-series decomposition model that decomposes a time series into three components: seasonal, trend, and residue.These three components are added to the original time series (Equation 1).
Where, Y t is a non-linear and non-stationary discrete time series defined on time t from 1 to N. N is the total number of observations in the time series.S t , T t ; and R t are seasonal, trend, and residue components, respectively.
STL applies LOcally wEighted regreSsion Smoother (LOESS) to the input time-series data.Suppose we have a dependent variable as y and an independent variable as x.A LOESS smoothing curve g x ð Þ is fitted on y.For calculating g x ð Þ at every value x i (i ¼ 1toN), q observations of y are taken around the x i such that q � N. The weight of each observation in q is calculated by using a tricube weight function.More distant observations in q will have low weights as compared to weights assigned to nearer observations.The weights for observations outside q are zero.STL has an inner and outer loop.The inner loop is nested inside the outer loop.Trend and seasonal components are updated with each pass of the inner loop, while robustness weights are updated with each pass of the outer loop.All robustness weights are set to 1 for first pass of the outer loop.This model needs six parameters to be defined, which are: The values of parameters were selected based on the characteristics of the time series (Cleveland et al. 1990).Defining the value of n p ð Þ is straightforward as it is equal to the number of observations per year in the LPTS.n s ð Þ is recommended to be equal to the n p ð Þ to make seasonal length equal to periodic length (Jiang et al. 2010).The value of The LST time series were decomposed by using 'stl' function available in the R programming language (R Core Team 2019).An example of decomposition of a VIIRS nighttime LST pixel time-series vector at location (23.69 N, 86.33 E) is shown in Figure S2 (see supplementary file).

Generation of deviation maps
Yearly deviation maps were generated using a neighbourhood-based concept.An N*N neighbourhood was created around the focal pixel.The background trend for the focal pixel is obtained by averaging the trend components of all pixels in the neighbourhood except the focal pixel.The background trend is subtracted element-wise from the trend component of the focal pixel.All elements of the resultant vector are annually aggregated by average.A positive deviation shows that the focal pixel has a higher value than its neighbourhood.In contrast, a negative deviation shows that the focal pixel has a lower value than its neighbourhood.This process is repeated for all the pixels in the study area by considering each pixel as a focal pixel.Yearly deviation maps for the years from 2013 to 2019 were generated with different neighbourhood sizes of 3 × 3, 5 × 5, 7 × 7, 9 × 9, and 11 × 11 for daytime and nighttime data separately.Therefore, a total of 70 deviation maps were generated.The product of the corresponding daytime and nighttime deviation maps is produced, leading to 35 product deviation maps.This process was implemented in the R programming language (R Core Team 2019).Daytime and nighttime deviation maps were produced to take advantage of the diurnal temperature variation phenomenon.For most features, thermal radiation emission is time-dependent, i.e. they emit higher thermal radiation during some parts of the diurnal cycle and lower thermal radiation in some other parts.Features like built-up land emit high thermal radiation during nighttime (Meng et al. 2018) and pose challenges in the segregation of coal fire pixels since the solar radiation is trapped in the built-up during daytime and this trapped radiation is emitted at nighttime, making total emitted radiation during nighttime equivalent to that of coal fires.However, the thermal emission is very low in the daytime from the built-up land.Such false alarms can be eliminated by combining their diurnal temperature characteristics.Therefore, nighttime and daytime deviation maps were produced, and their product substantially reduced the false alarms.The process for generating the deviation maps is shown in Figure 3.

Segregation of coal fire pixels
A threshold is determined for each product deviation map to separate coal fire pixels from non-coal fire pixels.A statistical method based on mean and standard deviation is employed to determine the threshold.Pixels with values greater than the limit (mean + 2*standard deviation) are flagged as anomalous.Pixels in product deviation maps are eliminated if their corresponding pixels in the daytime and nighttime deviation maps have values less than one.Such deviation values cannot belong to a potential coal fire pixel.This criterion was formed by analysing the daytime and nighttime deviation values of known coal fire pixels in the Jharia coalfield.Thermal anomaly maps are generated from all product deviation maps by employing the above-mentioned procedure.Therefore, a total of 35 anomaly maps are generated for the years 2013 to 2019, with five maps for each year belonging to different neighbourhood sizes of 3 × 3, 5 × 5, 7 × 7, 9 × 9, and 11 × 11.A pixel appearing anomalous for a particular year in at least four out of five anomaly maps was considered a coal fire pixel.

Results and discussion
Coal fire maps obtained for the years 2013 and 2019 are shown in Figures 4 and 5, respectively.A total of 123 and 127 coal fire locations were found in 2013 and 2019, respectively.Coal fire maps for the years from 2014 to 2018 are also provided in the supplementary file (Figures S3-S7).Coalfield-wise, the numbers of coal fire locations obtained for 2013 and 2019 are provided for major Gondwana coalfields (see the supplementary file, Figure S8).The individual coal fire maps for some coalfields are illustrated in Figure 6.Moreover, a map has been prepared to show new, dormant, and continued coal fire pixels (Figure 7).In this map, the coal fire pixels were categorised using coal fire maps of 2013 and 2019.Coal fire pixels that appeared in the 2013 map but disappeared in the 2019 map were termed as dormant coal fire pixels; coal fire pixels that did not exist in the 2013 map but appeared in the 2019 map were termed as new coal fire pixels.Coal fire pixels that existed in both the 2013 and 2019 maps were termed as continued coal fire pixels.Field surveys were carried out, and coal fire locations were visited in various collieries of Jharia and Raniganj coalfields to validate the results of the present study.However, for other coalfields, due to the unavailability of ground-based survey data, validation was done using high spatial resolution images available on Google Earth.Spatial association of detected coal fires and collieries was observed.Most of the detected coal fire pixels were observed over the opencast mining projects, which assured the reliability of the obtained results.
The product of corresponding daytime and nighttime deviation maps substantially weakened the false alarms like water, sand, built-up land, bare soil, etc., due to their varied diurnal temperature.Daytime, nighttime, and product deviation values for some of the false alarms are shown in Table 2.It is observed that water has a positive deviation during nighttime which is attributed to the high temperature of the water as compared to the background.In comparison, the same water pixel has a negative deviation during daytime due to the low temperature of the water as compared to the background.For pixels with major land cover as built-up, this technique is more effective as built-up land acts as an urban heat island during nighttime and temperature may increase as high as coal fire pixels, which may cause a false interpretation of built-up land as coal fire.However, low thermal radiations from built-up land during the daytime keep their temperature low; subsequently, the product of nighttime and daytime deviations becomes low.River bed sand is also a major false alarm; however, due to its negative deviation during nighttime and positive deviation during the daytime, the resulting product deviation value becomes negative.Such pixels no longer remain anomalous.Therefore, it becomes easy to identify and eliminate such pixels in the product deviation map.Threshold values for all yearly product deviation maps with different neighbourhood sizes are provided in Table 3.It was    interpreted that the threshold value increased with the neighbourhood size.This increase in threshold can be attributed to an increase in the deviation of the focal pixels from the background since increasing the neighbourhood size leads to a weak background trend.

Validation
The details of sites visited during the field survey in 2018-19 are given in Table 4.A total of 17 sites were visited in various collieries of Jharia (13 sites) and Raniganj (4 sites) coalfields.Since one pixel of the VIIRS VNP21A2 product covers a 1 km 2 area on the ground, the obtained results were validated within a 1 km radius around the visited coal fire location.
In other words, the obtained results were considered valid if the obtained coal fire locations fall within a 1 km radius of the visited coal fire locations.The proposed algorithm achieved an accuracy of 88%.Some detected coal fires were not temporally consistent; therefore, they are not found in all the coal fire maps generated from 2013 to 2019.Such behaviour can be attributed to the decreasing intensity of coal fires or diminishing coal fires during the observation period (Mujawdiya et al. 2022).

Conclusions
A satellite-sensor independent contextual coal fire detection algorithm was proposed to detect coal fires in parts of the Gondwana coalfields of India annually from 2013 to 2019 using time series daytime and nighttime VIIRS LST maps.The proposed algorithm derives a threshold inherently from the data and avoids the problems faced by fixed-thresholding approaches.Automatic thresholding makes this algorithm the best choice when quick detection of coal fires at the regional scale is needed.Furthermore, most of the previous research has been limited to Jharia and Raniganj coalfields; therefore, the present study attempts to detect coal fires in all major Gondwana coalfields, which will be helpful for future researchers to study the unexplored Gondwana coalfields.The proposed algorithm is independent of the location and size of the study area; therefore, this algorithm can be used globally.This algorithm uses time series of images in contrast to limited images used in most of the previous coal fire detection algorithms available in published literature.It minimises the chances of false identification of coal fires due to poor data quality.LST time-series data provide the long-term LST trend at coal fire pixels, which are robust to external disturbances, such as cloud cover, sensor artefacts, bad weather conditions, anthropogenic activities, etc., compared to a single observation of coal fire pixels.LPTS vectors of all pixels in the study area were decomposed by employing the STL model to extract the trend components.Yearly daytime and nighttime deviation maps were generated by using the extracted trend components.The product of daytime and nighttime maps significantly weakened false alarms.Moreover, the neighbourhood-based technique used for generating yearly deviation maps makes the proposed algorithm successful in detecting coal fire-induced thermal anomalies at the regional scale.Using multiple neighbourhood sizes for generating yearly deviation maps makes the proposed algorithm more robust against false alarms.The accuracy of the proposed algorithm was calculated using ground-truth data collected during field surveys.The proposed algorithm achieved good accuracy of 88%.However, the present algorithm has some limitations in removing false alarms completely.In the Godavari Valley coalfield, false alarms were discovered on the riverbed of the Godavari River.Such false alarms can be removed by employing the present algorithm to high spatial resolution time-series LST maps.An inventory of industrial plants can be prepared to mask false alarm pixels produced by them.

Figure 1 .
Figure 1.Map of Gondwana coalfields located in India.Background Image is a true colour high resolution satellite image used as a base map.
available in rescaled digital numbers that must be multiplied by 0.02 to obtain actual LST values.The LAADS web interface provides VIIRS images in tiles.Due to the large size of the study area, a total of three tile locations were covered by the study area.Each tile location had 322 daytime and 322 nighttime LST maps.Therefore, the corresponding tiles were mosaicked to form full images covering the entire study area.LST maps may have poor quality LST pixels due to cloud cover, haze, aerosols, and so on.Such LST pixels in LST maps were filtered based on quality bands accompanying the LST maps.Pixel values in the quality band range from 0 to 255.The quality of the retrieved LST can be measured based on the pixel value in the corresponding quality band(Hulley et al. 2017).The LST is not generated for some values (e.g. 2, 3, 6, 10, and 14) of the quality band.Poor or marginal quality LST is generated for quality bands with pixel values like 13, 29, 77, and 93; such LST values are identified and removed from LST maps.

Figure 2 .
Figure 2. Flow-chart of the proposed algorithm.

Figure 3 .
Figure 3. Algorithm for generating daytime and nighttime deviation maps.

Figure 4 .
Figure 4. Coal fire map obtained using proposed algorithm for the year 2013.The base map is a VIIRS NIR band of October 2013.

Figure 5 .
Figure 5. Coal fire map obtained using proposed algorithm for the year 2019.The base map is a VIIRS NIR band of October 2019.

Figure 7 .
Figure 7. Location of new, dormant, and continued coal fire pixels.Coal fire pixels were categorised into new, dormant, and continued categories using coal fire maps of 2013 and 2019.

Table 1 .
Thermal bands used in the TES algorithm for retrieval of LST.

Table 3 .
Threshold values for yearly product deviation maps with different neighbourhood sizes.

Table 2 .
Deviation values for the year 2016 using 7 × 7 neighbourhood size.