Simulating wildfire spread based on continuous time series remote sensing images and cellular automata

Study area

The research area was selected within localised burn extents of the 2020 MuLi Forest Fire, where time sequential multi-source remote sensing images were collected, as shown in Fig. 1. The geographical coordinates of MuLi County are (100.139859E, 28.942699N), encompassing a total area of 13,252 km², with an administrative district boundary measuring 134.81 km. MuLi exhibits typical mountainous, hilly and canyon-like landforms, with three primary topographical features: northwestern mountainous terrain, southeastern middle-height mountain terrain and southwestern high-mountain terrain. According to meteorological data from the China Meteorological Administration (CMA), the county experiences distinct dry and wet seasons, with an annual average temperature of 14°C. The fire started at 19:30 hours on 28 March 2020, with the ignition point located at coordinates (101.379492E, 27.900137N). The prevailing wind directions were from the northeast and northwest, while the average wind speed ranged from 2.5 to 3.0 m s⁻¹.

Fig. 1.

Location of the study area and the fire site: (a) the location of the study area in Sichuan Province, (b) the fire area in MuLi County, and (c) the fire scene in MuLi County.

Preparation of the data

The present study incorporates various spatiotemporal factors influencing fire spread, encompassing wind speed, wind direction, distribution of vegetation, thermophysical parameters of vegetation and topographic conditions. The collected data include fire ignition points (fire ignition lines), a digital elevation model (DEM), wind speed and direction, and Sentinel-2A satellite images, as well as the vegetation map of China. Table 1 presents a summary of the data sources.

We use ArcGIS 10.1 software to conduct data manipulations and analyses for pre-processing purposes. The DEM is clipped to the study area boundary and converted into the raster data representing slope and aspect. The spatial distribution information of vegetation is derived by an instance segmentation method (Zhao et al. 2022) based on Sentinel-2A remote sensing images, aiming to assess the diversity of vegetation within the study area. Further, thermo-physical parameters of various vegetation fuels are obtained through systematic sampling and measurement. Dynamic meteorological data every 3 h from 29 to 31 March 2020 are collected within the study area to consider the influence of meteorological parameters. Finally, all the data are incorporated into the geospatial information database, defined as a ‘quartel grid’. The process of data preparation is schematically illustrated in Fig. 2.

Fig. 2.

Data preparation process. Topography encompasses slope and aspect; vegetation encompasses distribution and parameters; meteorological factors encompass wind speed, wind direction, temperature and humidity; the data store comprises two pivotal aspects: multi-precision and multiple forest regions.

Vegetation inversion

For the study area, we introduce multi-source and multi-temporal remote sensing images to construct comprehensive datasets for vegetation inversion. To establish ground truth, we adopt the vegetation map of China (1:1,000,000) (Wu et al. 2022), the Global Land Cover Product with Fine Classification System at 30 m for 2020 (GLC_FCS30-2020) (Meng et al. 2022), and up-to-date data from Google Digital Image. Subsequently, we utilise Labelme (Russell et al. 2008) to construct labels. The train and test datasets are constructed using Sentinel-2A remote images, encompassing infrared spectral bands across different seasons. Specifically, the vegetation map of China (1:1,000,000) and the GLC_FCS30-2020 provide information on vegetation types within various compartments, while Google Digital Image and labels offer the boundary delineations of vegetation. Moreover, we divide the built training and test dataset into three subsets: 70% for the training dataset, 25% for the test dataset and 5% for the validation dataset.

We obtain the spatial distribution information of vegetation by combining the maximum likelihood (ML) technique (Otukei and Blaschke 2010) and the Hybrid Task Cascade (HTC) algorithm (Chen et al. 2019). The ML technique is a non-linear supervised classification method based on Bayes’ theorem, which necessitates the derivation of prior probability and conditional probability density functions from provided training samples. The HTC algorithm is a multi-task and multi-stage hybrid cascade instance segmentation approach that enhances information flow by integrating cascaded stages and multi-tasking capabilities at each stage.

The ML method integrated within ArcGIS and the HTC model are utilised for vegetation inversion, with the fusion of results achieved by incorporating ground truth data and manual priors. Moreover, non-combustible objects such as rivers and villages are identified and categorised using the abovementioned approach to serve as natural barriers for fire spread within forested regions. These natural barriers are excluded from the outcomes to enhance the reliability of the inversion.

Fire spread models

The fire spread modelling incorporates the Rothermel (1972) and the Time-WZF rate models (Sun et al. 2021), considering the topographical slope factor in MuLi County. The study area is partitioned into a grid of square-shaped cells with uniform dimensions (i.e. sides measuring 30, 50 and 100 m). The ROS R is determined using Eqn 1 when slope ≤25° as recommended in the Rothermel model (1972):

where the reaction intensity I_R (kJ m⁻² min⁻¹) represents the heat release per unit area of the fire front, ξ denotes the no-wind propagating flux ratio, Φ_s is a slope factor, Φ_w is the wind coefficient, ρ_b (kg m⁻³) is the oven-dry bulk density and Q_ig (kJ kg⁻¹) is the heat required to ignite a unit weight of fuel. The wind coefficient Φ_W is calculated using:

where

The slope factor Φ_S is calculated according to Eqn 4:

where β is the filling rate, ϕ is the slope and σ is a function coefficient.

The ROS R is calculated according to Eqn 5 when 25°< slope ≤ 60° as recommended in the Time-WZF model (Sun et al. 2021):

where R₀ is the initial ROS (m min⁻¹), K_φ is the wind coefficient, K_θ is a terrain factor, K_s is the combustible index (obtained from a look-up table (Wang 1992)); K_r is the time correction coefficient, set to 0.04 as suggested in Rui et al. (2018); a = 0.03, b = 0.05, c = 0.01, d = 0.3; T is the temperature (°C), W is the wind level, Int a set integer, RH is the air humidity (%), v is the wind speed (m s⁻¹), θ is the slope angle, and g is the direction of the hill (1 for uphill and −1 for downhill). The model takes into account the effects of combustibles, wind, temperature, humidity and slope on the ROS.

Beyond the two models above, there is scope for incorporating additional models. The selection of the most appropriate model at each simulation step depends on the precise conditions of applicability, including vegetation type, slope and wind speed. As pointed out previously, the Rothermel model is employed for a slope less than 25°, whereas the Time-WZF model is selected when the slope falls within the range of 25°–60°.

Time-adaptive cellular automata

The CA model used in each grid incorporates eight directions (east, west, south, north, southeast, southwest, northwest and northeast). The rule for fire propagation encompasses both the interior and exterior of the grid, as illustrated in Fig. 3. The state of cells either changes or remains constant according to the local evolution rules. The state of a cell at a discrete time step t_step is determined by whether the propagation displacement L_spread exceeds the cell size along the spread direction L_size. The L_size is calculated with the grid size L, L_size = L/2 in east, west, south and north directions, Lsize=L/2 in southeast, southwest, northwest and northeast directions. Two local evolution rules are taken into account as follows:

Fig. 3.

Fire propagation outside the grid and within the grid. Rule 1 corresponds to the outside grid, and Rule 2 to the grid inside.

Rule 1: the fire spreads outside the grid if the propagation displacement L_spread exceeds the cell size along the spread direction L_size. For the coordinate (x₁, y₁) of the fire point in current grid, the coordinate (x₂, y₂) in next grid during t_step can be calculated with:

where R is the ROS calculated by the model within the grid, and the coordinates (x_i, y_i) are the central coordinates of the grid considered as the coordinates of fire points. As illustrated in Fig. 3, when the propagation displacement L_spread along the northeastern direction exceeds the grid size L_size = L/2, this indicates that the fire point has extended from the centre of one grid to the centre of an adjacent grid.

Rule 2: the fire spreads within the grid if the propagation displacement L_spread is smaller than the size of the cell along the spread direction L_size. The propagation displacement L_spread is first calculated, followed by computation of the remaining time for spreading over the grid t_rest. Subsequently, the new time step is introduced into the spread calculation with Eqn 7 until the rule changes from 2 to 1:

where L_spread is the propagation displacement, L_size is the grid size along the spread direction, t_rest is the extra time required to spread out of the current grid and tstepupdate is the new time step for a new simulation within the grid in this direction. According to Rule 2, the propagation displacement L_spread is evaluated against the grid size along the spread direction (L/2) to adjust the duration within this grid, as shown in Fig. 3. This process continues until the fire point has successfully exceeded the current grid in the corresponding direction. During this phase, the coordinates of the fire point (x_i, y_i) are no longer fixed at the grid centre and change over time. Owing to the total simulating time t_total being constant, Rule 2 may lead to inconsistent output times for simulation outcomes when compared with Rule 1. To address this discrepancy, the present study employs computational programming methodologies to ensure a consistent visualisation of simulation results across different time frames.

Assuming S is the state of the grid, according to Eqns 6 and 7, the transition from one state to another is calculated by:

where Si,jt is the state of the grid (i, j) at the current moment, Si±1,j±1t+Δt is that of the grid (i ± 1, j ± 1) at the next moment, Ri±1,j±1t is the ROS for each direction, which is approximately equal to Ri,jt in the ROS model within the grid approximately, L_size is the grid size along the spread direction, Δt is the value of time step, t is the current time.

A critical condition judgment factor C_f reflects the relationship between propagation displacement and grid size within the current cell based on Eqn 8:

As time progresses discretely, the state of the neighbouring grid is determined by summing the state of the burning grid at the previous moment and the critical condition judgment factor C_f, adhering to the two local evolution rules. This process is depicted in Fig. 4. The critical condition judgment factor C_f represents the state of the current grid, defined as: C_f = 0, unburned; 0 < C_f < 1, burning, following Rule 2; and C_f ≥ 1, fully burned, following Rule 1. The variable S denotes the state of the adjacent grid at different time points, where S = 0 or S = 1 indicates that the adjacent grid state remains unburned, 1 < S < 2 indicates that the current grid has caused the adjacent grid to enter a burning state, and S > 2 indicates that the current grid has led the adjacent grid to reach a fully burned state. Wildfire propagation is often influenced by external factors such as fluctuating wind and active fire suppression efforts, resulting in slow propagation in the local direction. The proposed TCA model incorporates localised rules both within and outside the grid to achieve a controlled and deliberate fire spread in the local direction. The details of the input and output parameters of the TCA model are presented in Table 2.

Fig. 4.

Grid transition dynamics: fire propagation model, the adjacent grid states: unburned (a–b), burning (c) and fully burned (d). The centre S_{(i, j)} represents the current grid with a fire point, surrounded by adjacent grids. The light green, orange and grey grids indicate the unburned, burning and fully burned states, respectively. The black bold arrows represent the transition of the grid state, while the red and grey arrows indicate the directions of wildfire spread.

Extraction of fireline

The location of the fireline is determined by employing the Normalized Burn Ratio (NBR) to identify burned areas through remote sensing images, as Epting et al. (2005) recommended.

where ρ_nir denotes the near-infrared band, ρ_swir represents the short-wave infrared band and the NBR index has a theoretical value ranging from −1 to 1.

The initial spatial information of the fireline is corrected using time series of remote sensing images through an adaptive simulation method, which involves integrating the extracted fireline coordinates with the proposed TCA model. Specifically, the NBR method is initially utilised to extract the coordinates of the fireline contour from t₀ to t_n−1, and subsequently calculate the average ROS in each direction, denoted R₀, R₁, R₂, …, R₇ based on the extracted information. Moreover, the fireline coordinates from t₀ to t_n−1 are used as the initial spatial information of fire spread at t_n, and R₀, R₁, R₂, …, R₇ are adopted as the initial ROS instead of the results from the rate models. Finally, the proposed TCA model is integrated into a simulation system for wildfire propagation. The algorithm effectively addresses the challenges posed by non-linear factors such as dynamic wind and manual fire suppression.

Development of simulation system

We evaluated the performance of the proposed method by developing a wildfire propagation simulation system based on ArcGIS Engine by C++, C# and XAML. A comprehensive static database that stores crucial topography, meteorology and vegetation parameters is incorporated into the system. The Rothermel model and Time-WZF model are embedded to calculate the ROS within the grid. The system integrates the proposed TCA model and utilises geographic information system (GIS) technology to achieve the three-dimensional and dynamic visualisation of fire propagation. Fig. 5 illustrates the technology pathway for the system.

Fig. 5.

The technology route for the simulation of wildfire spread.

The input conditions of the simulation system encompass the time step, total time and coordinates of a single fire point and fireline. The output parameters comprise the ROS, fireline length, burned area and fireline coordinates. Further, the fireline coordinates extracted from multi-source remote sensing images could be imported into the system as new initial spatial information to achieve adaptive fire spread simulation in subsequent time steps. The system also has a fireline identification module, which provides real-time spatial coordinates of firelines to simulate wildfire spread.

Vegetation inversion results

By utilising the ML method integrated within ArcGIS and the HTC model, spatial distribution information on vegetation can be acquired with comprehensive datasets generated for training and testing purposes. The accuracy of predicting the probability of rivers reached 95%, while that for villages approached approximately 90% when compared with the GLC_FCS30-2020 and Google Earth System Digital Image. Moreover, the accuracy in predicting the probability of dominant vegetation was ~75% compared with ground truth labels. Importantly, inversion analysis provides more comprehensive information on the fire environment than previous research (Tavakol Sadrabadi and Innocente 2023; Xu et al. 2024), encompassing dominant vegetation, rivers and villages, as depicted in Fig. 6.

Fig. 6.

Vegetation inversion result: (a–b) are the inversion results with different grid sizes. The white areas indicate rivers and villages, while the other coloured areas represent vegetation.

The inversion results provide crucial spatial distribution information on key ground objects, thereby highlighting the importance of fine-grained inversion in wildfire spread simulation.

Simulation results without and with the TCA model

This section compares the simulation results without and with the proposed TCA model. The proposed model integrates two rules for higher and lower fire spread rates. The primary focus of the comparative analysis lies within the designated area highlighted in yellow shown in Fig. 7. The simulation results with the presented model effectively capture the gradual fire propagation trend in local regions. The proposed TCA model effectively addresses the limitations imposed by the time step and grid size, resulting in simulation outcomes that exhibit greater consistency with real wildfire environments. The local fire spread is specifically limited in the leeward direction of the fire front, thereby effectively preventing the fire front from surpassing its current grid within a given time step. Comparative analysis indicates that the proposed method effectively addresses the constraint imposed by grid size and time step on the spreading algorithm, thereby enabling a localised slow spreading simulation.

Fig. 7.

Simulation results without and with the proposed model, (a) 2D simulation, (b) Local details from (a), (c) 3D simulation. The first row represents the simulation result without the TCA model, whereas the second row corresponds to the simulation result with the TCA model. The red point denotes the ignition source, and the black points are generated by the spread of the wildfire. The blue line delineates the simulation result of the fire line. The yellow squares represent the local details of wildfire spread and the grey areas represent burned regions and the orange line is the wildfire boundary in 3D terrain.

Results of fireline extraction

The time series of fire point distribution is derived from the GF4 and Sentinel 2A remote sensing images in the MuLi Forest Fire, as shown in Fig. 8. The ignition occurred at 19:30 hours on 28 March, as determined through investigation (Fig. 8a), with the initial coordinate of the fire point recorded as (101.379492E, 27.900137N). By 13:15 hours on 30 March, localised areas along the west started experiencing fire spread, contradicting the overall propagation trend of the fire (Fig. 8b), which can be attributed to local dynamic winds within the canyon. At 18:18 hours, there was a significant expansion in the fire-affected area (Fig. 8c). However, at 11:47 hours on 31 March, the fire intensified along the east (Fig. 8d), aligning with prevailing wind conditions. By 13:57 hours, the smoke plume from the fire displayed an upward trajectory (Fig. 8e), indicating a resurgence in fire intensity. Subsequently, at 15:12 hours, the fire spread towards the north–northeast direction with local weather conditions (Fig. 8f).

Fig. 8.

Forest fire monitoring in MuLi County based on GF-4 satellite data: (a–f) represent the fire spots at 19:30 hours on 28 March, 13:15 hours on 30 March, 18:18 hours on 30 March, 11:47 on 31 March, 13:57 hours on 31 March, and 15:12 hours on 31 March, respectively.

According to the analysis above, uncontrollable factors such as dynamic wind patterns and human intervention in fire suppression often pose significant challenges for quantifying localised fire spread outcomes during fire propagation, which presents substantial obstacles for fire propagation simulation. Therefore, the extraction of fireline coordinate information based on continuous time series remote sensing images (Fig. 9) and research on parameter initialisation for spreading simulations demonstrate significant importance.

Fig. 9.

Multi-source satellite remote sensing image extraction of the firelines: (a–c) represent the fire spots at 11:47 hours on 31 March, 13:57 hours on 31 March, and 15:12 hours on 31 March, respectively.

The location and shape of the fireline were extracted from multi-source remote sensing data (Fig. 9). The time series remote sensing images at 11:47, 13:57 and 15:12 hours on 31 March were selected based on the extracted fireline coordinates to achieve adaptive simulation. Considering the prevailing local weather conditions, particular attention was focused on the northeast fireline along the wind direction. The coordinates of the firelines at 11:47 hours (Fig. 9a) and 13:57 hours (Fig. 9b) were acquired using ArcGIS based on WGS_1984 UTM Zone 48N, as presented in Table 3. Specifically, the time denotes temporal data regarding distinct fire perimeters, the fire points represent a sequence of fire points along the current time on the fireline, direction indicates the primary direction of fire spread, and the coordinates provide the location information of the aforementioned fire points.

Based on the results in Table 3, the ROS is calculated with the fire points coordinates at t₀, including t_0–1, t_0–2 and t_0–3, as well as the fire points at t₁, namely t_1–2, t_1–7 and t_1–13. The coordinates of fire points from t_1–1 to t_1–14 are also utilised as new fireline information for fire spread simulation at t₂.

Simulation results based on parameters initialisation and TCA

Based on the above extracted fireline results, the fire spread trend observed in the remote sensing images from t₁ to t₂, combined with directional information provided by the TCA model show that the fire spread towards the north, northeast and east during this period. Consequently, three group fire points at t₁ and t₂ were selected to calculate the ROS in these directions.

After calculations, the spread distances in each direction are 805 m to the north, 991 m to the northeast and 934 m to the east. The rates of fire spread from t₁ to t₂ are determined as 0.103, 0.127 and 0.120 m s⁻¹, respectively. To initialise the input parameters for simulation at t₂, we assign the ROS as 0.103 m s⁻¹ in the northern direction, 0.127 m s⁻¹ in the northeast direction and 0.120 m s⁻¹ in the eastern direction for all fire points ranging from t_1–1 to t_1–14. The initial fireline for simulation at t₂ is selected based on the fire points within coordinates from t_1–1 to t_1–14, as presented in Table 4. This experiment successfully initialised the coordinates of the firelines and the ROS in each direction for simulating the fire spread from t₁ to t₂, overcoming the limitations associated with calculating ROS based on geographical parameters within a grid.

The dynamic simulation of wildfire spread is achieved through the initialisation of parameters and the application of the TCA model. The simulation of MuLi Forest Fire is presented over a duration of 75 min from t₁ to t₂, with a time step interval of 25 min. The primary wind direction is set as eastward. The coordinates of fire points and their corresponding spread rates in each direction are derived from Table 4. The simulation results are shown in Fig. 10. A comparison between the dynamic simulation results and multi-source remote sensing data, as presented in Fig. 9, indicates that the dynamic simulation accurately captures the observed trends in the burning area, with consistent shape changes reflecting those of the actual burning area. This consistency serves to validate the efficacy of the proposed method for simulating real wildfire events. According to Table 4, it is evident that the rates of spread in the eastern and northern directions are comparable. Consequently, the propagation outcomes should be similar within the same time step for both directions. However, observational data suggest a greater extent of eastward propagation compared with northward propagation. This behaviour may be attributed to the amalgamation of multiple simulation results along the wind direction, resulting in an expanded dispersion area. To accurately integrate simulation results from multiple fire points, this study calculates the propagation correction factor along the wind direction at t₂ based on the firelines extracted from t₀ to t₁.

Fig. 10.

Dynamic simulation results by initialisation parameters. (a) Fire points by fireline. Simulation results (b) after 25 min, (c) after 50 min, (d) after 75 min. The simulated results coalesce on intersection.

For a qualitative assessment of the fire spread simulation outcomes, the simulation results were compared with the ground truth derived from time series remote sensing images, and the performance of the proposed method was evaluated based on the burned area and fireline length. As shown in Fig. 11, the overall spread trend, shape of the burned area and location of the fire front exhibit remarkable consistency with ground truth, thereby substantiating that time series remote sensing images demonstrated exceptional correction performance for the simulation. The effectiveness of the proposed TCA model is further confirmed by successfully simulating the gradual propagation scenario within a localised region.

Fig. 11.

Qualitative comparison of simulation results and real information. The grey area represents the simulation results, while the area bounded by the red line depicts the actual data. The period of the simulation is from t₁ to t₂.

To further demonstrate the advantages of the proposed method, we conducted a comparative analysis between the simulation results of burned area and fireline length and their corresponding ground truth values, as presented in Table 5. The findings reveal that prediction accuracy of the burned area reached an impressive 92.2%, while the accuracy for fireline length prediction achieved 87.3%. The experimental results demonstrate that the proposed methodology effectively enhances the simulation precision and successfully addresses significant prediction errors caused by objective factors such as local dynamic wind patterns and human intervention in firefighting efforts.

Subscripts