Study Area
The Sulabulak fire area is located in the north-central region of Xinjiang, China, approximately 86 km south of Urumqi(Fig. 1). This area lies between latitudes 87°27′ and 87°32′ and longitudes 43°20′ and 43°23′. The terrain is characterized by mid-low mountain hills, with elevations ranging from 1995 to 2211 meters, and relative height differences of 50 to 110 meters. The region has a continental climate with annual precipitation between 170.4 and 201.1 mm, and predominantly experiences southwest winds. The geological structures of the Sulabulak fire area include the Jurassic Xishanyao Formation (J2x) and Quaternary sediments. The coal seams, which are highly flammable due to their low sulfur, low phosphorus, high calorific value, and high melting temperature, consist mainly of black, blocky, or layered coal with well-developed joints and fractures. Coal fires were discovered during field surveys in 2016, primarily caused by incomplete stripping and underground mining that exposed shallow coal seams to spontaneous combustion. These fires have progressively worsened, causing significant environmental damage and economic losses.
Datasets
This study utilized 29 Landsat Collection 2 Level 2 images from 2013 to 2023. Released by the US Geological Survey (USGS), this dataset comprises enhanced Earth observation products from the Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS) on Landsat 8 and Landsat 9 satellites (Masek et al., 2020). These datasets, featuring surface reflectance and temperature products, are processed with advanced radiative transfer models and supplementary atmospheric data to enhance accuracy and usability (Claverie et al., 2018). Available in image formats of approximately 185 km by 180 km, the data utilize Universal Transverse Mercator (UTM) projection, or Polar Stereographic (PS) projection for polar regions (Dwyer et al., 2018). Stored as 16-bit unsigned integers, the datasets include metadata for converting digital numbers (DN) to radiance, reflectance, and brightness temperature. These datasets, available via the USGS Earth Explorer portal, USGS EROS Machine-to-Machine (M2M) API, and AWS Simple Storage Service (S3), offer direct access for detailed analysis (Crawford et al., 2023).
This study utilizes Sentinel-1 data from the European Space Agency (ESA). Sentinel-1, a C-band satellite with a 12-day revisit cycle (Torres et al., 2012), provided 135 images acquired between January 13, 2018, and December 13, 2023. The standard acquisition mode was Interferometric Wide (IW) with dual-polarization (VV/VH). Precise orbit corrections used data from the ESA Sentinel-1 Quality Control website. Additionally, 30-meter spatial resolution SRTM data served as external reference digital elevation models (DEMs) for InSAR processing. For SBAS-InSAR processing, 428 interferograms were generated, with temporal baselines ranging from 12 to 60 days and the longest spatial baseline being 122.1 meters. The average baseline length for PS-InSAR was − 157.2 meters. Temporal and spatial baseline graphs from SBAS-InSAR and PS-InSAR processing are shown in Fig. 2.
Vegetation cover retrieve
Accurate vegetation cover estimation is crucial for understanding the impact of coal fires on the local environment, particularly in identifying smoke fugitive channels. Vegetation cover inversion leverages the spectral absorption and reflection characteristics of vegetation across different spectral bands (Chen et al., 2024). This study employs the pixel dichotomy method to extract the Fractional Vegetation Cover (FVC) for the study area over the period from 2013 to 2023, reflecting the vegetation coverage status in the region.
The Normalized Difference Vegetation Index (NDVI) is calculated using the near-infrared and red bands of remote sensing data. Based on this, the Fractional Vegetation Cover (FVC) is calculated using the pixel dichotomy model. The formulas are as follows:
$$\:\begin{array}{c}NDVI=\frac{{R}_{NI}-R}{{R}_{NI}+R} \left(1\right)\end{array}$$
$$\:\begin{array}{c}FVC=\frac{NDVI-NDV{I}_{soil}}{NDV{I}_{veg}-NDV{I}_{soil}} \left(2\right)\end{array}$$
where NDVI is the Normalized Difference Vegetation Index, RNI is the spectral value of the near-infrared band, and R is the spectral value of the red band; NDVIsoil represents the NDVI of bare soil; and NDVIveg represents the NDVI of fully vegetated areas. The values for NDVIsoil and NDVIveg are selected from the 5%-95% confidence interval.
After obtaining the FVC values, the vegetation cover for each period is reclassified using ArcGIS software. This reclassification is based on the field vegetation conditions in the area, ensuring that the remote sensing data aligns with ground truth observations. The number of pixels and the area for each interval are then statistically analyzed to understand the temporal dynamics of vegetation cover.
This approach enables a detailed analysis of vegetation cover changes over time, providing insights into the environmental impact of coal fires and aiding in the identification of smoke fugitive channels.
Temperature retrieve
In this study, we used ENVI software to monitor land surface temperature (LST) with products from the United States Geological Survey (USGS) Landsat-8 Collection 2 Level 2. The LST data is captured by the Thermal Infrared Sensor (TIRS) on the Landsat 8 satellite and undergoes preprocessing, including atmospheric correction and surface emissivity correction (Galve et al., 2022). The LST calculation process follows these steps:
$$\:\begin{array}{c}L\left(\lambda\:\right)={M}_{L}\times\:{Q}_{cal}+{A}_{L} \left(3\right)\end{array}$$
$$\:\begin{array}{c}T=\frac{{K}_{2}}{\text{ln}\left(\frac{{K}_{1}}{L\left(\lambda\:\right)}+1\right)} \left(4\right)\end{array}$$
$$\:\begin{array}{c}{T}_{s}=\frac{T}{1+\left(\lambda\:\times\:\frac{T}{\rho\:}\right)\times\:\text{ln}\left(ϵ\right)} \left(5\right)\end{array}$$
The digital number (DN) is first converted to spectral radiance L(λ). In formula 3, ML is the radiance multiplicative factor, Qcal is the calibrated pixel value (DN), and AL is the radiance additive factor. Using the Planck equation in formula 4, we calculate the brightness temperature. T is the absolute temperature (in Kelvin), and K1 and K2 are pre-launch calibration constants. Formula 5 uses a modified Planck equation to convert the radiance back to actual temperature, where TS is the land surface temperature (in Kelvin), λ is the wavelength (in meters), and ρ is defined as h*c/σ, with h being the Planck constant, c the speed of light, and σ the Stefan-Boltzmann constant; ϵ is the surface emissivity (Meghraj et al., 2023).
InSAR deformation retrieval
SBAS-InSAR is an advanced time-series Synthetic Aperture Radar (InSAR) analysis method that constructs short-time baseline differential interferograms by selecting and combining collected data. This approach efficiently overcomes spatial decorrelation by selecting target points based on spatial coherence coefficients and building a linear model. Using Singular Value Decomposition (SVD) or the Least Squares Method (LSM), linear deformation rates and elevation error values for the target points can be obtained, resulting in a comprehensive deformation time series over the entire observation period (Selvakumaran et al., 2018). SBAS-InSAR technology excels in processing long-term observational data and monitoring large-area deformation, making it especially suitable for detailed ground deformation studies.
Persistent Scatterer Interferometry (PS-InSAR) is a technique used for monitoring ground subsidence and deformation across various regions. It identifies stable scatterers, known as Persistent Scatterers (PS), which exhibit consistent phase behavior over long periods. PS-InSAR combines high precision and efficiency, enabling continuous monitoring of ground deformation, particularly useful for detecting slow-moving deformations. This technique provides high temporal resolution deformation information, offering detailed measurements of ground subsidence and deformation (P. Zhang et al., 2023).
When used together, PS-InSAR and SBAS-InSAR allow tracking of various subsidence rates, significantly enhancing monitoring density and accuracy (Ramzan et al., 2022). The technical workflows for SBAS-InSAR and PS-InSAR are illustrated in Fig. 3. These complementary methods provide thorough and precise monitoring of ground deformation, essential for comprehending and managing the dynamics and effects of underground coal fires.
Detailed processing methods for SBAS and PS techniques can be found in (Tian et al., 2023; Z. Zhang et al., 2023). These methods provide robust technical support for studying smoke fugitive channels and surface deformation in underground coal fire areas, essential for precise monitoring and developing effective mitigation measures.
Integrated Smoke Channel Detection
Previous research shows some overlap between temperature anomalies and surface deformation (Zhou et al., 2013; L. Jiang et al., 2011). However, temperature anomalies caused by underground coal fires exhibit temporal lag, non-linear characteristics, and significant seasonal variation (Song & Kuenzer, 2014). Directly overlaying these results can lead to data omissions and errors.
To address this issue, this study integrates high-temperature threshold extraction and analysis of sparse to moderately sparse vegetation coverage, supplemented with surface deformation data, to verify and accurately locate smoke fugitive channels in coal fire areas.
Figure 4. Frequency superposition for time series images to extract coal fire pixels
The high-temperature threshold (Th) is derived from temperature inversion data. It is defined as the sum of the mean surface temperature (Tm) and twice the standard deviation (Tsd), serving as the optimal threshold to distinguish between coal fire areas and non-coal fire areas (W. Jiang et al., 2011), as shown in formula 6. To obtain more accurate temperature threshold information for coal fire areas, we incorporated thermal interference data from non-coal fire regions collected through field surveys, along with the location and thermal information of actual coal fire boundaries. We conducted a comparative filtering experiment to extract fire area ranges at different time series. As illustrated in Fig. 4, the high-temperature threshold attribute was set to 1, while other areas were set to 0, and then the data were superimposed. A pixel was identified as an abnormal temperature region when its frequency of appearance in the same location and within the set period reached the predefined high-frequency filtering threshold.
$$\:\begin{array}{c}{T}_{h}={T}_{m}+2{T}_{Sd} \left(6\right)\end{array}$$
By employing this method, this study can more accurately identify areas affected by underground coal fires, thereby providing a more reliable scientific basis for the monitoring and management of coal fire disasters.