Characterization of Background Temperature Dynamics of a Multitemporal Satellite Scene through Data Assimilation for Wildﬁre Detection

: Detection of an active ﬁre in an image scene relies on an accurate estimation of the background temperature of the scene, which must be compared to the observed temperature, to decide on the presence of ﬁre. The expected background temperature of a pixel is commonly derived based on spatial-contextual information that can overestimate the background temperature of a ﬁre pixel and therefore results in the omission of a ﬁre event. This paper proposes a method that assimilates brightness temperatures acquired from the Geostationary Earth Orbit (GEO) sensor MSG-SEVIRI into a Diurnal Temperature Cycle (DTC) model. The expected brightness temperatures are observational forecasts derived using the ensemble forecasting approach. The threshold on the difference between the observed and expected temperatures is derived under a Constant False Alarm Rate (CFAR) framework. The detection results are assessed against a reference dataset comprised of MODIS MOD14/MYD14 and EUMETSAT FIR products, and the performance is presented in terms of user’s and producer’s accuracies, and Precision-Recall and Receiver Operating Characteristic (ROC) graphs. The method has a high detection rate when the data assimilation is implemented with an Ensemble Kalman Filter (EnKF) and a Sampling Importance Resampling (SIR) particle ﬁlter, while the weak-constraint Four-Dimensional Variational Assimilation (4D-Var) has comparatively lower detection and false alarm rates according to the reference dataset. Consideration of the diurnal variation in the background temperature enables the proposed method to detect even low-power ﬁres.


Introduction
Data acquired from satellite sensors are widely used in the automatic detection of active wildfires. Fire detection techniques supported by such data provide the most feasible and practical means for the detection and continuous monitoring of regional and global fire activities [1]. Though data used in the design of satellite-based fire detection have been acquired from sensors onboard sun-synchronous Low Earth Orbit (LEO) systems (mainly polar satellites) and onboard Geostationary Earth Orbit (GEO) satellites, thermal imagery data obtained from the global system of GEO sensors are the most appropriate satellite-based data for early fire detection in most regions of the earth susceptible to wildfires. The devastating effects of wildfires on land, atmosphere, flora, and fauna can be mitigated when wildfire outbreaks are timely detected to enable early interventions.
In thermal-based fire detection, satellites have relied on images received from the fire-sensitive Mid-Infrared (MIR) spectral window (i.e., 3-5 µm spectral range) to distinguish fire and non-fire pixels. This spectral region is highly sensitive to temperatures of vegetation fires-covering even a subpixel scale-and is less affected by solar reflections than the Shortwave Infrared (SWIR) spectral window. A pixel in a satellite image is directly or implicitly flagged as a fire pixel when the pixel's observed temperature is significantly higher than its expected background temperature-the latter defined as the expected temperature of a pixel under non-fire condition. The first satellite-based remote sensing techniques for hotspot detection were developed using data from the LEO systems. Given the LEO system's high spatial resolution advantage, the expected background temperature was extracted from spatial-contextual information (for example, see [2][3][4][5][6][7]). A pixel's expected background temperature can also be estimated using a multitemporal model for the evolution of a pixel's background temperature. It can also be defined using multispectral information either by assuming that the temperatures derived from the fire-sensitive MIR spectral window channel and temperatures derived from the earth surface-sensitive Long-Wave IR (LWIR) spectral window channel during non-fire conditions are correlated or by performing spectral unmixing [8,9]. Since GEO fire-sensitive MIR data were available, most GEO-based fire detection techniques were conceived by adapting the pre-existing LEO methods on GEO data. These techniques used the same spatial-contextual information to define the expected pixel's background temperature during fire detection [10][11][12][13][14] or fire confirmation [15][16][17]. A spatial-contextual mechanism has its drawbacks, such as Point Spread Function (PSF) effect, spatial heterogeneity, and undetected cloud-contaminated pixels that contribute to both omission and commission errors [1,18,19]. The performance of such techniques is also restricted by the low spatial resolution of GEO data. GEO sensors have a high temporal resolution, and exploiting the high temporal resolution enhances the chance of detecting a fire at ignition. Given that fires are generally small when started, the high-temporal data can then allow the detection of small and low-power fire events, typically only possible with the LEO system [20][21][22]. In this study, GEO satellite temporal MIR data are assimilated into the model of Diurnal Temperature Cycle (DTC) parameters of a pixel to describe non-fire background temperature dynamics to aid the early detection of wildfires. Temporal information incorporates land surface characteristics, but to also consider the weather and atmosphere, a combination of temporal, spatial, and spectral information must be assimilated into the background temperature model of a pixel.
Multitemporal fire detection algorithms define the expected pixel's background brightness temperature in MIR in different ways. The detection statistic termed intensity index [23,24] defines the pixel's expected background temperature as the mean of observed temperatures at the same time of the day and month of the year. Likewise, the temporal test for fire confirmation proposed by [17] defines the pixel's expected background temperature as the mean at the same time of the day in a series of thirty days, and the temporal test in the fire detection method proposed by [25,26] that exploits the polar-orbiting NOAA-AVHRR MIR spectral channel, is determined as the mean of the brightness temperatures over a number of days prior to the inspected day. The change detection algorithm proposed by [27,28] defines the expected difference between a pixel's background temperature and an auxiliary temperature at the current time as the difference of observed brightness temperature and auxiliary temperature at an anterior non-fire time. The time to an anterior time can be fixed to the minimum time between successive images up to and not exceeding twenty-four hours, and the auxiliary temperatures are modeled as temperatures that represent the natural evolution of temperatures and can be, among others, non-fire atmospherically corrected temperatures of the current day or previously observed DTC with no cloud or fire-contaminated samples. In computation of a detection statistic named 'time differential index', a fire detection method proposed by [24] predicts the expected difference between a pixel location's background temperature at the current time to its background temperature at a previous time-the previous time and the current time are at the same time of the day and month of the year. Instances of cloud cover or missing samples in the available data up to the inspection time hinder the derivation of the expected background temperature in these multitemporal fire detection methods. The problem can be solved by interpolating the invalid values due to cloud cover or fire events. The interpolation is achievable by fitting available valid brightness temperatures to a physics-based or data-driven DTC model to produce the expected DTC of the current day of a pixel location that can be used to define the expected background temperature at an inspection time. Given the presence in observed temperatures at the Top Of the Atmosphere (TOA) of outliers, mainly cloud or fire-contaminated samples, the fitting can be undertaken by first detecting and removing outliers. However, instances of missed detection of clouds or fire events negatively affect the results of the estimation of the expected background temperatures. Another option is to use methods from robust statistics to avoid the need for clearing outliers from the data. A robust fit of a DTC model on past observed DTCs has been used to represent a pixel's background temperature [29][30][31]. Some fire detection techniques apply, in one of their multiple tests, a dynamic threshold on observed MIR brightness temperature, and the threshold varies with the solar zenith angle [15,28,32]. This dynamic threshold that is also used by the current version (i.e., the third version) of the EUMETSAT FIR product [13,33] implicitly represents the pixel's background temperature dynamics. Similar to the robust fit on previous DTCs, the dynamic threshold avoids the influence of both cloud and missed fire events on the result of the fire detection, but does not consider the particularity of a day: the dynamic threshold does not change with the day, and the robust fit on previous DTCs assumes that the previous days define the whole of the current day.
To include all past temporal information on a pixel when deriving its expected brightness temperature, methods based on recursive Bayesian estimation adjust the expected brightness temperature of a pixel each time new information becomes available. The background temperature is derived by a Kalman filter-based prediction using one of the physics-based or data-driven DTC models as the forecast model that determines the time evolution of the background temperature [29,[34][35][36]. The new valid observed brightness temperature, when available, is assimilated into the background temperature model. The expected background temperature can also incorporate both multitemporal and spatial information. For example, the spatial differential index proposed by [24], regressing the inspected image in the current scene on a number of past images of the scene [37,38], expected background temperature of a pixel as a convex combination of a spatial-contextual estimate and a temporal estimate [39], and a fire detection method proposed by [40] determines the expected rise of a pixel's background radiance with respect to its neighborhood's background temperature from a recent image, acquired at the same time of the day, of a non-fire condition on the pixel location and its neighbors. The background temperature can also be derived from information extracted from both multitemporal and multispectral information. For example, the method proposed by [41] predicts MIR background temperature from the earth sensitive LWIR temperature using a linear relation derived from the principal components of a temporal sequence of non-fire MIR and LWIR brightness temperatures. The current multitemporal fire detection methods based on recursive Bayesian estimation, first derive the expected DTC of a day (i.e., the model of a day) online before the start of the day and data are assimilated at observation time into the expected DTC. The derivation of the expected DTC of a day allows only for the inclusion of recent observed DTCs with a short period of invalid data (cloud or fire-contaminated data or missing data) into the training set but observed DTCs with a long period of invalid data are not included. The model of a day is also expected to change with land cover/use change, season, and even weather. Given that the model of a day is derived from previously observed DTCs, the required number of recent valid DTCs to include in the training set might not be available for a given period of no change in a pixel location. This approach, therefore, does not adapt rapidly to non-fire change on land or in the atmosphere, for example, in the aftermath of a fire on a pixel location. This paper, with the aim of rapid adaptation to normal change, proposes a method that derives the background temperature by assimilating data into the model of DTC's parameters of a pixel location each time valid observations are available to estimate the vector of the best DTC model parameters, which is the initial value of the forecast. The ensemble forecasting approach is applied from the data assimilation time for a skillful forecast that is used to determine the expected background temperatures of a pixel location onward.
The method proposed in this paper implements an early active fire detection by subtracting the observed temperature from the expected background temperature to obtain the prediction residual at an inspection time in a given pixel location. The residual value is compared to an adaptive threshold computed using a Constant False Alarm Rate (CFAR) framework to decide whether the pixel location is burning. The dynamics of the background temperature are described using a state-space model. The variables of the forecast model are a subset of the parameters of the DTC model proposed by [42] (this model is referred to as the BER06 model in this paper), and the forecast model variable vector is modeled as a Gaussian random walk process. Two types of data are assimilated into the model of the DTC parameters at each data assimilation cycle, namely, the observed TOA brightness temperatures derived from the MIR window data acquired from the METEOSAT Second Generation (MSG)-Spinning Enhanced Visible and Infrared Imager (SEVIRI) sensor and the observed offsets of the thermal sunrise (i.e., the start time of a DTC) of a day to the sunrise of the same day. The latter is assimilated to transition from one DTC to the subsequent DTC. To produce a smooth continuous data assimilation result in the temporal dimension, three data assimilation methods-each offering different assimilation advantages-were separately used. A variational data assimilation method, which emulates the weak-constraint Four-Dimensional Variational Assimilation (4D-Var) and is implemented only in the time dimension (one-dimensional analysis), performs a point estimation of the posterior. Two sequential data assimilation methods, namely the classical Ensemble Kalman Filter (EnKF) and Sampling Importance Resampling (SIR) particle filter, estimate the posterior probability density function (pdf). The three methods are compared in terms of detection performance that depends on their ability to estimate the dynamics of the background brightness temperatures. The expected background brightness temperature is the observational forecast determined through model integration from the initial state and defined at a given False Alarm Rate (FAR). Observational forecasts are assumed to be uncorrelated across spatial and temporal dimensions, and the same for real observations. The observational forecasts and real observations are also assumed uncorrelated, and as result, each pixel is assessed independently from its neighbors and only the temporal dimension used to compute the expected background temperature. The proposed method has less computational complexity in terms of time than methods that rely on fitting past observations to a DTC model to determine the expected background temperature [29,31,[34][35][36]. The computation time is reduced because the model fits are only required when the proposed detection system is launched to determine the initial values at launch. The proposed method is expected to achieve higher accuracy than that of other temporal fire detection methods because the estimation of background temperature is based even on the most recent non-fire observation without relying on the expected DTC of a day, and the detection uses a CFAR threshold to adapt to any non-fire condition in a single pixel location and to adapt to any region. Since the MIR spectral window is affected by solar reflections, cloud cover increases both omission and commission errors. Cloud detection and the masking of cloud-contaminated samples are usually undertaken as a pre-processing step, and fire detection results depend then on the accuracy and detection threshold of the cloud detection mechanism of choice. In this study, no step is taken to clear cloud samples from data. The potential improvement to the proposed method includes the assimilation of data in more dimensions than just the temporal dimension and the assimilation of more types of data to better characterize the background temperature dynamics under cloud cover. The identification of a water body is implemented as a post-processing step. This is done to set water body pixels, which might have been falsely classified as fire or undecided, to a non-fire status. The detection is assessed in terms of user's accuracy (i.e., the performance in terms of commission error), producer's accuracy (i.e., the performance in terms of omission error), and Precision-Recall and Receiver Operating Characteristic (ROC) graphs using data on a region of Southern Africa during the dry winter season of 2007. The detection performance of the proposed method is assessed and discussed using EUMETSAT FIR and MODIS MOD14/MYD14 active fire products.

Fire Detection Algorithm
Observed data are assimilated into the background DTC model. In this study, the data assimilation is implemented at each time a new observation is flagged as a non-fire observation (i.e., fire-contaminated samples are not assimilated). The data assimilation aims to derive better initial conditions of forecasts. The data assimilation and the forecast model of the background temperatures are presented in Section S1 of the Supplementary Materials to this article. The expected background temperature of a pixel is derived from the forecasts initiated at the time of the assimilation of data. The forecast is achieved using the ensemble forecasting scheme, described in Section S1 of the Supplementary Materials of this article. The shortest-range ensemble forecasts of brightness temperature verifying at the inspection time k in a pixel location represent the distribution of the background (non-fire ) temperature of the pixel location at time step k. The mean value of this ensemble is the expected value of the background (non-fire ) temperature of the pixel location at a time step k. The initial condition of the forecast is derived using the weak constraint 4D-Var with a data assimilation window of length N w = 6, the EnKF with number of ensemble N e = 51, and the SIR particle filter with number of particles N p = 51.
The input to the fire detector comprises the observations (brightness temperatures in IR 3.9 MSG spectral channel) at the inspection time, not yet assimilated. It also comprises the expected distribution of the background temperature (non-fire temperature) of each pixel location at the inspection time. A pixel is classified as belonging to either the fire class or the non-fire class using an anomaly detection (one-class classification) mechanism. At each inspection time, only statistics on the non-fire class are available-i.e., this study models only the pixel's non-fire condition and not the fire condition. The brightness temperatures in the IR 3.9 MSG spectral channel y 1 (t k ) ≡ y 1,k are mapped into the feature space that easily enables identification of the two classes with fewer fire omissions and false alarms. The aim is to detect the start and end of a fire with minimum delay and to monitor fire in real time. The extracted feature (or derived sufficient statistic for the detection) at time step k in a given pixel location is the difference between the pixel's observed temperature and the pixel's expected background temperature. The difference is the prediction residual expressed as where x f k is the forecast state, y 1,k is the observed brightness temperature, and h 1,k (·) is the observation operator relating the state vector to the observed brightness temperature. The predictable component of the brightness temperature (the expected background temperature) is filtered out from the observation to get the prediction residual. The (conditional) distribution of the prediction residual under non-fire class (i.e., p(y 1,k |x k , condition = non-fire )) is estimated. The classifier is applied to a prediction residual and outputs either non-fire pixel (do not reject H 0 ) or fire pixel (reject H 0 ) based on a decision boundary (threshold T k ) in the feature space provided by a CFAR detector whose design is detailed in Section 2.2. The fire detection algorithm is summarized in Algorithm 1. When the analysis is implemented using the weak-constraint 4D-Var, the feature extraction (residual generation) and the fire reporting start at k = 6 (end of the first analysis cycle), and five first samples are assumed to be under non-fire condition. When the analysis is achieved by the EnKF and the SIR particle filter, the fire reporting starts at k = 1. Hereafter, we refer to the fire detection method according to the method used during the data assimilation stage. The weak constraint 4D-Var with a data assimilation window of length N w = 6, the EnKF with the number of ensemble N e = 51, and the SIR particle filter with the number of particles N p = 51 belong, respectively, to 4D-Var-CD (N w =6) method, EnKF-CD (N e =51) method, and the SIR-CD (N p =51) method.

Design of a CFAR Based Fire Detection
The optimal threshold on detection statistics (i.e., the prediction residual r k ) can only be derived when both the distribution of the fire (target) temperature and the distribution of the background temperature (non-fire temperature) are provided or estimated. Given that the background temperatures of a pixel location are modeled in this study (state-space model given by Equations (S2), (S5) and (S6)), and the distribution of background temperature can then be estimated while the distribution of fire temperature is unknown (i.e., no information on parameters of fire conditions is provided), a convenient design criterion that can be used to design fire detection system in online detection format consists of minimizing the delay for detection (that results in few omissions) for a fixed probability of false alarm (i.e., CFAR). The CFAR detector sets an adaptive threshold on the feature data r k based on the allowed level of false alarm probability (i.e., FAR). We assume that the non-fire prediction residual r k in one-pixel location over a lapse of one day long is independent and identically distributed (i.i.d) (a range of one day is taken to satisfy a compromise between the time before the land surface change and sufficient samples to characterize a time step k). Then, the CFAR window, R * ,k , under H 0 , contains 96 consecutive fire-free samples (in a temporal CFAR mechanism considering every single pixel independently) generated up to time step k − 1. Since in a given moving window of 96 samples, some samples can be fire contaminated, samples are selected from time step k − 1 down to the launch time of the detector, skipping fire-contaminated samples. Thus, the need for a censoring procedure is then lessened.
A simulation test was conducted to select the most suitable parametric model of the probability distribution of non-fire prediction residuals in a reference set R * ,k (the pdf of R * ,k under H 0 is denoted by f R * ,k (r k |H 0 )) and to model the relationship between the threshold coefficient γ and the false alarm rate P FA , given that where r 0k is the standardized variate derived from the original variate r k , and R * ,k is a random variable whose possible value is r k . The simulation was conducted by randomly drawing 2500 land pixel locations among land pixels not affected by fire (non-fire pixel) according to MODIS MOD14/MYD14 and EUMETSAT FIR products over a period of seven days, from 30 July 2007 at 12:00 a.m. UTC, and per each pixel location, randomly drawing 20-time instants from 31 July 2007. A uniform draw by a simple random sampling (with replacement) was performed using the Mersenne Twister pseudorandom number generator, which was initialized once at the start of the simulation (of each analysis method); and the simulation counts then a total of 50,000 realizations. Each pixel is described by a sequence of immediate previous 96 samples of r k of its location, which forms a set of generated prediction residuals under H 0 (i.e., one realization of the simulation). The suitable parametric model f R * ,k (r k |H 0 ) is asserted, taking into consideration that the medium was not cleared of outliers (e.g., clouds, fire-contaminated samples not reported by MODIS MOD14/MYD14 and EUMETSAT FIR products) and occurrences of possible abrupt changes in the DTC models from one day to the next (e.g., an unexpectedly hot day). After inspection of histograms of R * ,k , the distribution of the test statistic, r k , was assumed to be unimodal, and of Location-Scale (L-S) type for the detection to conform to the CFAR requirements [43]. Based on the high proportion of the sets whose distributions have an excess kurtosis (G 2 ) strictly less than 3 (i.e., a moderate kurtosis), and also based on the values of the first quartile, third quartile and the minimum of G 2 of the sets (presented in Table 1), Gaussian distribution and other L-S types, leptokurtic family of distributions with a moderate kurtosis can model the distribution of the prediction residual. The proportion of the sets with a skewness (G 1 ) between −0.5 and 0.5, and the first and third quartiles of G 1 of the sets recorded in Table 1, show that the distribution of prediction residual exhibits a considerable negative skewness. However, Anderson-Darling (A-D) goodness-of-fit test (using Best Linear Unbiased Estimate (BLUE) parameter estimates of distributions), implemented per analysis method, did not reject the null hypothesis of Gaussian distribution, at 1% nominal significance level, of a high number of the sets of prediction residuals compared to the number of the sets that were not rejected under the null hypothesis of the Gumbel distribution for minima (reverse Gumbel) distribution at the same significance level (see observed rejection rate of H 0 presented in Table 1). Student's t location-scale distribution is also generally applicable to the region of interest for non-fire pixels. However, a global constant shape parameter (i.e., the number of degrees of freedom ν) must be determined per analysis method. The prediction residuals under H 0 when the analysis is implemented using 4D-Var-CD (N w =6) can also better be modeled by reverse Gumbel distribution. In this study, Gaussian distribution was selected as a single parametric model of the detection statistic r k under H 0 when any of the three different methods of analysis is used to initialize the forecast. It is defined by its BLUE of the location parameter φ L(R * ,k ) and BLUE of the scale parameter φ S(R * ,k ) extracted from the reference set R * ,k at time step k. After the launch of a detection system, before 96 non-fire prediction residuals in the window R * ,k are available, the location parameter is assumed to be zero at time step k and the scale parameter is the standard deviation of the prediction residuals estimated beforehand via the simulation test. For each of the randomly drawn 2500-pixel locations, and each analysis method (each pixel had 569 samples starting from 30 July 2007 at 5:57 a.m. UTC (this is the Level 1.5 image time)), the variance of the predictions was found first for each pixel location. Then, the median of these variances was retained. The estimated global variances are listed in Table 1. Table 1. Characteristics of the distribution of the detection statistic r k (one MC run). Given equivariant estimators, φ L(R * ,k ) and φ S(R * ,k ) , for the location and the scale, respectively, at time step k estimated from the reference set R * ,k , the detection is declared if the prediction residual at time step k exceeds the threshold T k given by By keeping the threshold coefficient γ constant, the false alarm rate P FA is expected to be conserved in time and space. The threshold coefficient for r k (or the (1 − P FA )-quantile of the empirical distribution of r k ) under H 0 , used for this temporal detection, is modeled as where P FA is the admissible false alarm rate. The number of samples in R * ,k must also be an explanatory variable to this model, but, given that this number is constant, it is ignored. Parameters a, b, and c of the threshold coefficient model are obtained by MC simulations (50,000 runs from 2500 random pixel locations and 20 timestamps per pixel location, explained above in this section). For each MC run, 96 prediction residuals are standardized with location and scale parameters obtained by BLUE. The empirical Cumulative Distribution Function (CDF) is found for each random instant time of a random pixel location, and (1 − P FA )-quantile function is determined considering γ > 0. A Least Squares (LS) model fit on the threshold coefficient function given by Equation (4) is implemented using a Nelder-Mead simplex method with the initial parameters being the parameters obtained by fitting (1 − P FA )-quantile function of the standard Gaussian distribution on the function model for P FA ranging from 10 −6 to 0.5 with a step of 10 −6 . The coefficients found (with their 95% sample confidence interval) are depicted in Table 2. Table 2. Fitting parameters (with their 95% (sample) confidence interval) for the threshold coefficient.

Experimental Data and Detection Algorithm Setup
For our experiments, we used the MSG-SEVIRI IR 3.9 images covering a region of Southern Africa spanning from (20 S, 23 E) to (33 S, 38 E) (precisely, from pixel (778, 822) to pixel (1153, 1094); i.e., a 376 × 273 MSG Level 1.5 image data) recorded late July 2007 to early August 2007 (i.e., during the winter fire episode of 2007). Each image of the region counts 22, 934 water pixels that fall in water bodies, identified through MODerate Resolution Imaging Spectroradiometer (MODIS) land/sea mask product (MOD03) and manual visual inspection. Estimation of background temperatures and detection algorithms run on all pixels of the study area, and water-body pixels are masked out after detection to report fires exclusively on land pixels. A fire that is detected in a single pixel at a given inspection time k is counted as one fire event, then instants of fire events are scrutinized. Fire detection results are reported per FAR and per analysis method. Moreover, no cloud screening is implemented (but noting that a very cloudy DTC or very fired DTC or water pixel DTC is more likely removed from the training set based on conditions set on DTC parameters highlighted in Table S2 in the Supplementary Materials), and no sunglint-contaminated sample screening is implemented. MATLAB is used to implement and test the proposed fire detection method (The source codes are publicly available online at https://github.com/satfiretech/fire-dataAssimilation).

Performance Evaluation Metrics
Results of fire detection are assessed and compared using two performance accuracies, user's accuracy (a measure of the commission error) and producer's accuracy (a measure of the omission error), given by user's accuracy(U) = #correctly reported fire pixels #reported fire pixels (5) producer's accuracy(P) = #correctly reported fire pixels #fire pixels (6) where the #S stands for the cardinality of a set S, '#fire pixels' is the number of fire pixels reported by the reference data product, namely MODIS MOD14/MYD14 and EUMETSAT FIR , 'reported fire pixels' is the number of pixels reported by the designed (proposed) fire detection algorithm (per analysis method), and 'correctly reported fire pixels' is the number of fire pixels according to the reference data products correctly reported (classified) by the designed algorithms. The performance of the fire detection method is also assessed through two performance curves, a graph of user's accuracy against producer's accuracy (also known as Precision-Recall curve) and ROC curve. The Precision-Recall curve is sensitive to both class skews and varying unequal/asymmetric misclassification (error) costs (i.e., when the cost of misclassifying a sample from one class is much higher than the cost of misclassifying a sample from the other class) while the ROC curve is insensitive to both [44], and the class skews are evident in a fire detection application where the true instances of non-fire class occur much more often than the true instances of fire class. However, due to a high number of non-fire pixels compared to fire pixels, the Precision-Recall curves visually expose clearer the difference between the detection methods, in a given FAR range of interest, that is not noticeable in the ROC space [45]. The pixel reported as fire pixel by the validation fire products but is reported by the proposed fire detection method as a non-fire pixel or or has NaN (expected) background brightness temperature or has NaN observation (i.e., instant of completely missing data [missed pixels]), is counted as a missed detection.
If an MSG-SEVIRI pixel containing more than one MODIS fire pixel is not reported as a fire pixel by a proposed fire detection algorithm, then all MODIS MOD14/MYD14 fire pixels belonging to the MSG pixel are missed (i.e., counted only under 'fire pixels'), while if reported as a fire pixel, all MODIS fire pixels belonging to the MSG pixel are detected (i.e., counted under 'correctly reported fire pixels' and 'fire pixels'). The process of deriving the best admissible CFAR at which a fire detection system can operate can be conducted during a phase that can be referred to as the 'detection calibration phase' using an available validation set [reference data] (in this study, the validation set is made solely of MODIS MOD14/MYD14 and EUMETSAT FIR reports). The operational CFAR (or the best CFAR for a detection method according to the reference data) can be determined based on different criteria. This study determines the operational CFAR, as a FAR, that minimizes an error function J, given by where p(fire ) and p(non-fire ) are the prior probability that a fire condition and a non-fire condition, respectively, will occur -they represent the proportion of correct fire or non-fire pixels to the sum of all pixels, and can be referred to as fire or non-fire true class distribution, respectively. p(fire ) and p(non-fire ) are representations of class skews. c(non-fire |fire ) and c(fire |non-fire ) are the misclassification cost of omission and cost of commission, respectively. FAR m is the measured FAR while the variable FAR is the design FAR (known also as the nominal or the estimated FAR). Noting that the detection cost function J(·) given by Equation (7) assumes that the cost of correct classifications is zero, and J(·) accounts of class skew [46].

Fire Detection Performance
User's accuracy (U) and producer's accuracy (P) computed per estimated FAR are shown in Figure 1a and Figure  performs better than other detection methods in terms of U (gives less commission errors) except at the low estimated FARs. However, the low estimated FAR region was not considered during the design of the CFAR detector (the estimation of FAR starts from 10 −6 to 0.5 as indicated in Section 2.2). For the estimated FARs in the low FAR region (i.e., high threshold values), 4D-Var-CD (N w =6) can produce a lot of false alarms given that the assimilation window of 4D-Var is more likely to contain more missed fire samples. SIR-CD (N p =51) performs worse than both 4D-Var-CD (N w =6) and EnKF-CD (N e =51) in terms of U. At high estimated FAR, all methods give similar results in terms of U. Overall, SIR-CD (N p =51) performs better than other methods in terms of P (low omission error or high probability of detection), and at some estimated FARs, it performs more like EnKF-CD (N e =51) .    The performance of the detection methods when any of the three analysis is implemented is fairly modest (no method with a Precision-Recall above the shown diagonal/random classifier) and is not good due to a high number of false alarms according to the reference data products. EnKF-CD (N e =51) shows a better performance compared to other detection methods, as EnKF-CD (N e =51) dominates in Precision-Recall space, it does also dominate in ROC space as shown in ROC curves in Figure 2b.
A single operational FAR, when the analysis is implemented with any of the three methods considered in this study, is derived using the error (cost) function J(·) given by Equation (7). This error function considers both commission and omission errors and is a function of the estimated FAR. It is assumed in this study that c(non-fire |fire )p(fire ) = c(fire |non-fire )p(non-fire ) = 0.5 and the error function simplifies then to J = 1 2 [(1 − P) + FAR m ]. Figure 2c depicts a resulting J(·) for a single test set made of data spanning the period 31 July 2007-4 August 2007 with a reference dataset consisting of a combination of MODIS MOD14/MYD14 and EUMETSAT FIR products, and MODIS detection time is rounded to the nearest SEVIRI detection time. The optimal design FAR point is found among discrete points ranging from 10 −8 to 10 −1 (with a base-10 logarithmic step of 0.5). EnKF-CD (N e =51) gives the lowest error compared with other detection methods, and its optimal point 10 −3.5 FAR gives the operational FAR. The 10 −3.5 FAR point is closer to the optima of 4D-Var-CD (N w =6) and SIR-CD (N p =51) , and the error difference (∆J) from 10 −3.5 point to the optima of 4D-Var-CD (N w =6) and SIR-CD (N p =51) is 0.0172 and 0.0151, respectively.   When each detection not confirmed by EUMETSAT FIR product is confirmed by any MODIS MOD14/MYD14 report in 20 min before or after the detection time, the number of reported false alarms according to the validation reference data (the combination of MODIS MOD14/MYD14 and EUMETSAT FIR ) is reduced for all detection methods as shown in Figure 3, and the reported overall performance is then increased. In this case under consideration, fire events confirmed by the same MODIS MOD14/MYD14 report are counted as one 'correctly reported fire' during computation of both P and U, and the number of 'fire pixels' is the sum of the number of EUMETSAT FIR reports and the number found by subtracting the number of MODIS MOD14/MYD14 fire reports to the number of times one or more EUMETSAT FIR fire reports are found in the interval of 20 min around a MODIS MOD14/MYD14 report (i.e., 1 to 3 fire events). The results depict an instance the validation products are used to validate a whole fire course on a pixel location and not just a fire event at a given inspection time k on a given pixel location.
Results of a comparison between 4D-Var-CD (N w =6) , EnKF-CD (N e =51) and SIR-CD (N p =51) detection methods at 10 −3.5 FAR when only MODIS MOD14/MYD14 fire product is taken in the validation reference data (and MODIS detection time rounded to the nearest SEVIRI detection time) are depicted in Table 3. The number of true positives, according to a MODIS MOD14/MYD14 product, is presented (the number of true positives in Table 3, is presented in terms of the number of MODIS fire pixels). From the results presented in Table 3, we can conclude that SIR-CD (N p =51) performs better than other detection methods (throughout the 5 days) in terms of detection, according to MODIS MOD14/MYD14 product. EnKF-CD (N e =51) and SIR-CD (N p =51) at 10 −3.5 FAR achieve a true positive rate of 78.64% and 79.34%, respectively, of fire events reported in MODIS MOD14/MYD14 product. A moderate portion of MODIS MOD14/MYD14 reports is detected by EnKF-CD (N e =51) and SIR-CD (N p =51) methods, and the accuracy is restricted given that small and cool fires might not be detected due to the low spatial resolution of the SEVIRI compared to the MODIS resolution.   Table 3 for MODIS MOD14/MYD14 and true positive results presented in Table 4, the proposed fire detection algorithms detect MODIS MOD14/MYD14 fire reports at a lesser extent than they detect EUMETSAT FIR fire reports. Because, as mentioned above, of the small and cool fires that might not be detected due to the low spatial resolution of SEVIRI compared to MODIS. The degree of omission of EUMETSAT FIR fire reports by the proposed algorithms mainly depends on the level of the design FAR, and in this section, the design FAR is taken to be 10 −3.5 . However, the omission of EUMETSAT FIR reports at a given design FAR highly materializes while reporting fire events during times of extreme fires, and omission error arises in case the fire temperature dynamics at certain times are not far from the expected background temperature dynamics. Results of Table 4 show that EnKF-CD (N e =51) and SIR-CD (N p =51) perform relatively similar and better than 4D-Var-CD (N w =6) detection method. EnKF-CD (N e =51) detects EUMETSAT FIR reports at 94.10%, and detects the combination of EUMETSAT FIR fire events and MODIS MOD14/MYD14 fire events not reported by EUMETSAT FIR at 86.09%. SIR-CD (N p =51) detects EUMETSAT FIR reports at 96.48%, and detects the combination of EUMETSAT FIR fire events and MODIS MOD14/MYD14 fire events not reported by EUMETSAT FIR at 87.67%. All proposed detection methods report a significant number of false alarms according to MODIS MOD14/MYD14 and EUMETSAT FIR , and detect EUMETSAT FIR reports relatively well. Their detection performance is reduced from around t m (i.e., the time of the maximum thermal radiation) to around t s (onset of the nighttime temperature decay) with respect to EUMETSAT FIR and MODIS MOD14/MYD14 results, and such performance can be viewed to be more associated with a high number of fire events reported by both MODIS MOD14/MYD14 and EUMETSAT FIR products in the time range. During time of a high number of fire events, a proposed algorithm can miss a fire event reported by the fire products but can report a fire before or after the fire product reports. Furthermore, the observed temperature in MSG IR 3.9 spectral channel of an MSG pixel identified by MODIS to be affected by fire, might not deviate significantly from its expected background temperature. 4D-Var-CD (N w =6) has an overall low detection, especially during the nighttime-the time of high fluctuations in brightness temperatures due to different coupled land-atmosphere phenomena, such as orographic forcing effect. In the afternoon, its low detection is more explained by the DTC linearization from around t m to around t s (or likewise, based on the approximate observation operator given by Equation (S8)). The detection rate of MODIS MOD14/MYD14 report is low compared to the detection rate of EUMETSAT FIR reports, and SIR-CD (N p =51) and EnKF-CD (N e =51) perform better in detecting MODIS MOD14/MYD14 reports compared to 4D-Var-CD (N w =6) at the 4 MODIS overpass times (EnKF-CD (N e =51) is better than all detection methods during the night). Based on the MODIS MOD14/MYD14 report, all detection methods certainly detect both daytime and nighttime fires, with 4D-Var-CD (N w =6) to a lower extent. Noting that MODIS detection time is rounded to the nearest SEVIRI detection time, and the number of 'fire pixels' is counted in terms of MODIS pixels to determine the true positive rate (P), according to MODIS MOD14/MYD14 product.    The empirical distribution of MSG IR 3.9 brightness temperatures for MODIS MOD14/MYD14 reported fire pixels and EUMETSAT FIR reported fire pixels are shown in Figure 5a,b, respectively. MODIS detection time is rounded off to the nearest MSG-SEVIRI detection time. For MODIS MOD14/MYD14 reported fire pixels, the temperature value of an MSG pixel containing the centers of a number of MODIS fire pixels at a given time step k is counted as many times as the number of the centers of the MODIS pixels. For fire detection results of 4D-Var-CD (N w =6) method, producer's accuracy for each MSG brightness temperature value is shown in Figure 5c,d when the validation set is made of MODIS MOD14/MYD14 product and EUMETSAT FIR product, respectively. 4D-Var-CD (N w =6) is able to detect cooler fire compared to EUMETSAT FIR , down to approximately a brightness temperature value of 285 K (with P = 1), as shown in Figure 5c when only MODIS MOD14/MYD14 product is accommodated in the validation set. Some very hot fires are missed by 4D-Var-CD (N w =6) , especially at the SEVIRI's saturation point (which can be attributed to the linear approximation of the observation function and the data assimilation window that does not contain only clean samples). For fire detection results of the SIR-CD (N p =51) method, producer's accuracy for each brightness temperature value is shown in Figure 5g,h when the validation set is made of MODIS MOD14/MYD14 product and EUMETSAT FIR product, respectively. SIR-CD (N p =51) reports cold fires more frequently than 4D-Var-CD (N w =6) , and can detect fires with brightness temperatures below 285 K (with P = 1) as shown in Figure 5g. For fire detection results of EnKF-CD (N e =51) method, producer's accuracy for each brightness temperature value is shown in Figure 5e,f when the validation set is made of MODIS MOD14/MYD14 product and EUMETSAT FIR product, respectively. Similar to 4D-Var-CD (N w =6) and SIR-CD (N p =51) methods, the EnKF-CD (N e =51) method is able to detect cooler fire compared to EUMETSAT FIR , even for fires with a brightness temperature value below 280 K (with P = 1) as shown in Figure 5e when only the MODIS MOD14/MYD14 product makes the validation set.  An example of MSG IR 3.9 images at three MSG timestamps on 3 August 2007 is illustrated in Table 5 to portray detection reports of 4D-Var-CD (N w =6) , EnKF-CD (N e =51) , and SIR-CD (N p =51) against both MODIS MOD14/MYD14 (MODIS detection time is rounded off to the nearest MSG-SEVIRI detection time) and EUMETSAT FIR reports. In addition, the detection reports of MODIS MOD14/MYD14 and EUMETSAT FIR at the three timestamps are provided to present and describe MODIS MOD14/MYD14 and EUMETSAT FIR products. The table provides Level 1.5 observation time (expressed in Coordinated Universal Time (UTC)) in 'Time' column, the pseudocolor versions of the MSG IR 3.9 images in 'Pseudocolor image' column to better visualize the brightness temperature distribution in the images, and MODIS MOD14/MYD14 and EUMETSAT FIR fire reports overlaid on the MSG IR 3.9 images in the 'MOD14/MYD14 & FIR' column. A great number of MODIS MOD14/MYD14 reports are not reported by EUMETSAT FIR . Essentially, since in this study, we assume that a MODIS pixel belongs to an MSG pixel into which its center falls, if a burning part of the MODIS pixel location is in a different MSG pixel, some reports would be missed by a fire detection method that uses MSG-SEVIRI data. The table also presents detection reports of 4D-Var-CD (N w =6) , EnKF-CD (N e =51) and SIR-CD (N p =51) methods on the three images, overlaid on the corresponding MSG IR 3.9 images and placed under '4D-Var', 'EnKF' and 'SIR' columns, respectively. EnKF-CD (N e =51) and SIR-CD (N p =51) detect all EUMETSAT FIR reports while 4D-Var-CD (N w =6) missed, at 08:12 a.m.

Discussion
4D-Var-CD (N w =6) has the least level of false alarms while SIR-CD (N p =51) has the highest (see Figure 1a), but SIR-CD (N p =51) has the least number of missed detections (see Figure 1b). The Precision-Recall graphs shown in Figure 2a (when MODIS detection time is rounded off to the nearest SEVIRI detection time) and Figure 3 (when any MODIS fire event in 20 min around a reported fire detection is used to confirm detection), and the ROC graph shown in Figure 2b prove that both EnKF-CD (N e =51) and SIR-CD (N p =51) perform better than 4D-Var-CD (N w =6) . Both EnKF-CD (N e =51) and SIR-CD (N p =51) perform better in terms of detection but exhibit a high level of false alarms when the validation set is made of both MODIS MOD14/MYD14 and EUMETSAT FIR products. The ROC curves in Figure 2b also reveal a disparity between the measured FAR and the estimated FAR (the estimated FAR is varied from 10 −8 to 10 −1 with a base-10 logarithmic step of 0.5). A higher number of false alarms than expected can be attributed to, among others, true fire events missing from the validation fire products and the PSF effect on the MSG images with respect to MODIS MOD14/MYD14 report during the MODIS overpass. Though the proposed method registers false alarms in case of abrupt changes in temperature of a day due to cloud or smoke, they, however, minimize false alarms due to change in the model of temperatures from one day to the next-a problem inherent to previous temporal fire detection methods. False alarms in the results of the proposed method can be reduced in various ways, including the implementation of censoring mechanisms in the designed CFAR detector, lessening the PSF effect on MSG imagery, and incorporating spatial-neighborhood information and atmospheric effects into the model of the background temperature.
At an operational, best FAR found to be 10 −3.5 (see Figure 2c), a high detection rate of fire events reported by the EUMETSAT FIR product is observed, and MODIS MOD14/MYD14 reported fires are also detected but with limitations caused by the low spatial resolution of MSG-SEVIRI. EnKF-CD (N e =51) and SIR-CD (N p =51) have the lowest fire omission.
EnKF-CD (N e =51) and SIR-CD (N p =51) detect 94.10% and 96.48% of fires reported by EUMETSAT FIR , respectively, and both detect approximately 79% of fires reported by MODIS. Given that the MODIS sensor has a high spatial resolution (1 km nominal GSD between pixel centers at nadir) compared to the SEVIRI sensor (3.4100 km to 4.6800 km nominal GSD between pixel centers in the study area), the MODIS sensor is expected to detect fires covering a small area or low-temperature fires. EnKF-CD (N e =51) and SIR-CD (N p =51) detect, respectively, 86.09% and 87.67% of MODIS reports that are not reported by EUMETSAT FIR product. The EUMETSAT FIR product used in this study is produced from SEVIRI results of the EUMETSAT FIR baseline algorithm, the same as the MODIS algorithm of MODIS MOD14/MYD14 product, it considers spatial neighborhood information but not multitemporal information in deciding whether a pixel is a fire pixel. Even though the SEVIRI sensor has a low spatial resolution, the early detection of fires enhances the possibility of detecting a fire on small areas or/and with low temperatures at any time. For example, in the DTCs given in Figure 6 (and results reported in Table 6), on 2 August 2007, MODIS MOD14/MYD14 reports a brief fire during Aqua morning overpass, but did not report any fire at an earlier time in the same day (i.e., during Terra morning overpass). However, EnKF-CD (N e =51) and 4D-Var-CD (N w =6) , which reports the fire at the time of Aqua morning overpass, report the fire at an earlier time even during the time of Terra morning overpass (SIR-CD (N p =51) reports fires only during early morning and not during Aqua overpass). By visual inspection of the DTCs in Figure 6 Table 6, and in a 3 × 3 vicinity of the examined pixel, it reports fire events in two pixels, and in each  Given the diurnal variation of background temperature of a pixel location, the use of an adaptive threshold in the detection of fire enhances the accuracy of detection (see Figure 5c,e,g). The proposed algorithms are able to detect fires at any time of a day, and EnKF-CD (N e =51) tends to outperform all other detection methods in detecting night fires (see Figure 4c). An example given in Figure 7 (and results reported in Table 7) depicts the expected background temperature (i.e., predicted brightness temperature) and the state of a pixel (either fire or non-fire ) at 10 −3.5 FAR at each time step for two DTCs starting from 31 July 2007 and 1 August 2007 of an MSG pixel location (983, 901). No fire was reported in a 5 × 5 pixel array around the examined pixel by EUMETSAT FIR in the two DTCs, and MODIS MOD14/MYD14 reports a brief nighttime fire on 1 August 2007. Even though the MODIS MOD14/MYD14 report is a nominal-confidence fire (i.e., 32%) and is reported during the time of low diurnal temperature, EnKF-CD (N e =51) and SIR-CD (N p =51) detection approaches detected the fire report. EnKF-CD (N e =51) detected all MODIS MOD14/MYD14 reports and has a low level of false alarms, according to MODIS MOD14/MYD14 and EUMETSAT FIR reports. 4D-Var-CD (N w =6) missed the night MODIS MOD14/MYD14 fire report (due to its degraded detection performance during a time of high fluctuations in brightness temperature, mainly when the data assimilation window accumulates missed fire-contaminated samples during the first stage of a fire).  The proposed methods detect either brief or persistent fires. The persistent fires, throughout their duration, are well monitored at least from the time they are first detected (i.e., when the ignition is detected) to the last time the fire is detectable. The analysis methods proposed in this study adapt rapidly to the change in non-fire conditions (quick update of the DTC parameters). However, the detection methods that use these analysis methods might be slow in detecting the extinction of a real fire since the adaptation speed depends on the rate of change of the proportional between the forecast error variance and the observation error variance. Figure 8 (and results reported in Table 8) shows an example that helps analyze the performance of 4D-Var-CD (N w =6) , EnKF-CD (N e =51) , and SIR-CD (N p =51) in detecting a fire of great magnitude. Specifically, in this example, a fire that causes MSG-SEVIRI IR 3.9 saturation (i.e., a fire that evidently has a high intensity) and has a prolonged heating duration. It should be mentioned that the heating duration, together with the fire intensity, contribute to fire severity, and are among factors that are taken into consideration when quantifying the fire severity. In this example, fire extends over two DTCs and raises brightness temperatures in the MSG-SEVIRI IR 3.9 spectral channel to the saturation point of the IR 3.  For near real-time fire detection, an assessment is implemented per fire event (per pixel). However, for a general assessment, a fire can be considered in its entirety across space and time (i.e., assess per object per image or likewise per object per time step), and this can help, among others, overcome the problem of overlapping of adjacent MSG pixels in the assessment of the detection. Fire detection can also be assessed per time interval (or per a fire incident) on a single-pixel location, and this can help, among others, overcome problems induced by smoke and cloud (e.g., obscuring a burning pixel location at some time instants during a fire) during performance assessment.

Conclusions
This study defined the expected background temperature of a pixel at a given constant False Alarm Rate (FAR) as a forecast from the latest analysis. The brightness temperatures in the MSG-SEVIRI IR 3.9 band and the offset of the sunrise to the thermal sunrise time of a non-fire condition are assimilated into the forecast model of the parameters of the DTC. The paper presented three analysis methods, namely, weak constraint 4D-Var, EnKF, and SIR particle filter, and a forecast is implemented using the ensemble forecasting scheme. The expected background and the observation temperatures of a pixel are used to identify a fire pixel in image acquired at each time step. A temporal CFAR detector that provides an adaptive threshold on the detection statistics is a component of the proposed fire detection method.
The paper proposed a multitemporal fire detection method that considers the diurnal variation in the background temperatures of a pixel location in assessing whether an observed brightness temperature at a given time step is extracted from a fire or non-fire pixel. The method is then able to detect even small and cool fires at any time of a day. Each time a non-fire observation is acquired, the DTC model is updated to help track rapidly variations in background temperatures, which leads to high performance of the fire detection method. The experimental results show a high detection rate of both EnKF-CD (N e =51) and SIR-CD (N p =51) but with a high level of unconfirmed detections. 4D-Var-CD (N w =6) recorded relatively a high omission (though exhibits a low level of false alarms) because the 4D-Var data assimilation window might include cloudy or smoky contaminated samples, or missed fire samples.
The use of assimilation strategies that will improve the efficiency in terms of both quality of the analysis estimate and computational cost in high dimensional space will result in better estimation of background temperature and accuracy of the fire detection. The improvement in modeling the variations of DTC parameters over time will be beneficial in fast adaptation to non-fire background change and early detection of the start and end of fires. False alarms can be reduced by extending the forecast model and observation model to include other parameters, such as parameters related to the atmosphere, land surfaces, or other MSG spectral bands apart from IR 3.9 . The spatial dimension and pressure levels will also be considered to formulate a global model (or regional model), taking into consideration the statistical correlation between neighbor pixels. The design of a method to deal with the Point Spread Function (PSF) effect will also be a task of future work. Different ensemble forecasting methods can be exploited for a better forecast skill. Different detection statistics can be derived from ensemble forecasts to detect a fire as soon as possible and be able to track it as long as it is still active. The content of the validation reference dataset will be improved.

Acknowledgments:
The authors would like to thank EUMETSAT for making accessible MSG archived data which were used in this work-All MSG Data ©2007 and 2020 EUMETSAT.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript: