# New joint estimation method for emissivity and temperature distribution based on a Kriged Marginalized Particle Filter: application to simulated infrared thermal image sequences

Thibaud Toullier, Jean Dumoulin, Laurent Mevel

<figure><img src="image_1.png" /></figure>

To cite this version:

Thibaud Toullier, Jean Dumoulin, Laurent Mevel. New joint estimation method for emissivity and temperature distribution based on a Kriged Marginalized Particle Filter: application to simulated infrared thermal image sequences. Science of Remote Sensing, 2025, 11, pp.100209. 10.1016/j.srs.2025.100209. hal-04966544

HAL Id: hal-04966544

https://inria.hal.science/hal-04966544v1

Submitted on 25 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
New joint estimation method for emissivity and
temperature distribution based on a Kriged
Marginalized Particle Filter : application to simulated
infrared thermal image sequences

Thibaud Toullier<sup>a,*</sup>, Jean Dumoulin<sup>a</sup>, Laurent Mevel<sup>b</sup>

<sup>a</sup>Université Gustave Eiffel, Inria, COSYS-SII, I4S Team, Allée des Ponts et Chaussées, Bouguenais CEDEX, CS 5004 - 44344, , France

<sup>b</sup>Inria, Université Gustave Eiffel, COSYS-SII, I4S Team, 263 Avenue Général Leclerc, Rennes, 35700, , France

Abstract

This paper addresses the challenge of simultaneously estimating temperature and emissivity for infrared thermography in natural environment, aiming for near real-time performance. Existing methods, mainly in satellite observation field, rely on restrictive physical assumptions unsuitable for ground-based application context (Structures and Infrastructures monitoring). Other generic methods are nonetheless computationally intensive, making them impractical for real-time use. Our objective is to provide a method with effective real-time calculation performance while still giving results comparable to those reference methods under the same hypotheses, finally achieving both good accuracy and performance. The proposed method is based on a dynamical state-space modeling for the temperature, where the state vector is assumed to be split into a dynamic component for the temperature and a stationary component representing the emissivity. Then the dynamical component is estimated by a Kalman filter approach, whereas the parameterized model and the emissivity component are estimated through a particle filtering framework resulting in a bank of Kalman filters, also called marginalized particle

<footnote>*Corresponding author
Email addresses: thibaud.toullier@univ-eiffel.fr (Thibaud Toullier),
jean.dumoulin@univ-eiffel.fr (Jean Dumoulin), laurent.mevel@inria.fr
(Laurent Mevel)</footnote>
filter. A spatial assumption of homogeneity for the temperature yields to the addition of a Kriging step to the Marginalized Particle Filter to overcome the ill-posed nature of the problem and to compute the necessary physical estimates in a reasonable amount of time while providing fair results compared to reference methods from the literature.

A comparison with two state-of-the-art methods, MCMC and CMA-ES, is presented. The results indicate that the proposed method estimates the true value within a maximum deviation of 3K, similar to CMA-ES, while MCMC achieves a more accurate estimate with a maximum deviation of 0.5K. However, the computational efficiency of the proposed method is significantly improved, reducing the processing time by seven orders of magnitude compared to MCMC and three orders of magnitude compared to CMA-ES. This remarkable efficiency highlights the method's feasibility for real-time monitoring of temperature and emissivity.

*Keywords:* infrared thermography, temperature and emissivity estimation, natural environment, kriged Kalman filtering, marginalized particle filter

# 1. Introduction

## 1.1. Infrared thermography for structural health monitoring

The latest improvements on uncooled infrared detectors for thermal cameras have brought new opportunities towards remote surface temperature measurements of observed objects. In particular, civil engineering structures can be thermally monitored through multispectral infrared thermography during a long period of time in outdoor conditions, enabling new potential applications such as structures energy consumption optimization[1], train axles bearing hot boxes detections[2], etc. The literature on the topic of absolute infrared thermography in natural environment is still sparse since it is an emerging application with many bottlenecks, shown by the complexity of the modeling and the monitoring before validation of the involved phenomenon [3, 4].

Therefore, the temperature retrieval faces the lack of knowledge of the thermo-optical properties of the observed objects and more particularly of the emissivity [5]. The emissivity denoted by $\epsilon_{obj} \in [0; 1]$ is an intrinsic parameter of a given surface that represents its capacity to emit a radiative flux. The emissivity is a difficult parameter to measure or approximate due to its dependence on material's nature, roughness, angle of observation,
wavelength and temperature. In most applications, only the wavelength dependence of the emissivity is considered. Fig. 1a shows the complexity of measuring surface temperature of an object in natural environment.

<figure><img src="image_3.png" /><figcaption>Figure 1: Schematics representing the simplified problem and the principle of infrared thermography measurement in natural environment.</figcaption></figure>

In this paper and to simplify the problem, the focus is on the sole interaction of a camera and an object with emissivity $\epsilon$ and temperature $T$. We therefore assume that all others contributions of Fig. 1a (the sky radiance $L_{sky}$, atmospheric one $L_{atm}$, etc.) and geometric parameters are known or estimated otherwise [6]. The surfaces are assumed to follow a diffuse emission law and the emissivity dependence to temperature can be neglected for the domain of application considered. The simplified system is then represented in Fig. 1b.

To derive the measurement equation from this simplified system, let $S \subset \mathbb{R}^2$ be the image space domain and $\mathcal{T} \subset \mathbb{R}^+$ the time domain. By using the previous assumptions, the irradiance $\gamma(T, \mathbf{s}, t, \lambda_i)$ received at a given pixel $\mathbf{s} \in S$, a time $t \in \mathcal{T}$ and a spectral band $\lambda_i$ can be written as:

$$ \gamma(T, \mathbf{s}, t, \lambda_i) = \epsilon(\mathbf{s}, t, \lambda_i)L^{\circ}(T, \lambda_i) \quad (1) $$

with temperature $T$ that implicitly depends on time $t$ and pixel location $\mathbf{s}$ and $L^{\circ}(T, \lambda_i)$ (W m$^{-2}$ sr$^{-1}$ µm$^{-1}$) represents the self emission governed by Planck's law that depends on temperature and wavelength, [7]. We further assume that the observed object has a slow varying emissivity value through
time, leading to a stationary emissivity state<fnref n="1" />

By taking the logarithm of (1), the measurement equation can be then split into two parts, one stationary and the other one dynamic, thanks to the previous assumption:

$$ \ln(\gamma(T, \mathbf{s}, t, \lambda_i)) = \underbrace{\ln(\epsilon(\mathbf{s}, t, \lambda_i))}_{\text{stationary}} + \underbrace{\ln(L^\circ(T, \lambda_i))}_{\text{dynamic}} \quad (2) $$

By using Wien's approximation and Planck's law:

$$ \ln\left(\frac{\gamma(T, \mathbf{s}, t, \lambda_i)\lambda^5}{C_1}\right) = \ln(\epsilon(\mathbf{s}, t, \lambda_i)) - \frac{C_2}{\lambda_i T(\mathbf{s}, t)} \quad (3) $$

Equation (3) leads to an undetermined system of one equation and two unknowns ($\epsilon$ and $T$).

When dealing with multi-spectral data, one can think that measurements in $N_{\text{bands}}$ spectral bands will help for solving this equation. However, (3) shows that $N_{\text{bands}} + 1$ unknowns lie in this system of $N_{\text{bands}}$ equations due to the dependence of the emissivity on wavelength. As a consequence, the estimation of the temperature of the observed object cannot be done directly by solving this equation system; more information and/or hypotheses are needed. Such information can include prior on the probability densities of the parameters to estimate, physical assumptions, additional measurements, etc.

This ill-posed problem is known as the temperature and emissivity separation (TES) problem [7]. Although different methods have been developed (see the reviews from [8] or [9] for the comparison of the estimation of the temperature of the land surface temperature), using Bayesian inference with Monte-Carlo Markov Chains (MCMC) [10] or Kalman filtering with [11] and [12], the Maximum Entropy estimator [13] or even neural networks [14], the TES problem is still challenging and often addressed in the field of remote sensing. In fact, the developed methods often rely on particular application hypothesis or computations based on a particular materials' database, for neural network training for instance [15]. In particular, different methods for

<footnote n="1">¹This assumption helps us to simplify the equation system, by splitting the two quantities into a dynamic one and a stationary one. Having such stationary component gives us the opportunity to represent the emissivity as a noise term in the equations (see (11))</footnote>
the land-surface temperature estimation have been developed in the satelli-
tal domain. One can cite for instance methods developed with the NASA
airborne Thermal Infrared Multispectral Scanner that assume prior knowl-
edge on the land surface emissivity such as thermal log residuals, alpha log
residuals, and many derived methods using multi-channel informations for in-
stance [16, 17, 18, 19]. Other methods, without any prior on emissivity have
also been developed, making use of specific spectral bands to allow physicial
assumptions such as Temperature-Independent Spectral Indices (TISI) [20]
and its variants, assuming that TISI coefficients on particular channels will
not vary substantially between day and night [21, 22]. However, this meth-
ods need specific mid infrared and thermal channels, observing Earth with
satellite images, channels that may not be available in ground-based appli-
cations. For instance, in our applications cases, off-the-shelf mono-spectral
or multi-spectral uncooled infrared cameras (based on microbolometers) are
used, which may not match the capabilities of satellite-embedded detectors.
Similarly, the Normalized-Emissivity Separation (NES) method and then its
extension the Temperature Emissivity Separation (TES) method have been
developed and applied to ASTER images [23]. However, those methods rely
on a correlation whose coefficients are obtained from a material database.
As a consequence, this high dependency on the chosen database materials
makes the generalization to any site and any observation scale tricky [24].
Other methods developed for the LANDSAT-8 images make use of particu-
lar spectral bands allowing specific physical assumptions. They rely on band
10 brightness temperature, band 11 brightness temperature, red band and
NIR band, which are not available in our application case (see [25] or [26]).
Finally, one generic method can be found in the litterature, that can be ap-
plied to any test site and observation scale, using MCMC [10]. This method,
however suffers from being computationally intensive.

The bibliographical review highlights several challenges in directly apply-
ing temperature and emissivity estimation methods from satellite imagery
to ground-based applications. First, off-the-shelf infrared cameras used in
ground-based settings may not match the capabilities of satellite-embedded
detectors. Then, key physical assumptions, such as atmospheric properties
considered in satellite-based methods, are not directly applicable at ground
level. Similarly, methods that rely on material databases may also be af-
fected by this scale difference, as material properties observed at satellite
altitudes may not translate accurately to close-range measurements. Finally,
targeted applications demand a low-cost, near-real-time solution, which in-
troduces both computational and time constraints that further complicate
the adaptation of satellite-based methods.

Therefore, those methods may not be appropriate for the long-term ther-
mal monitoring of civil engineering structures through infrared thermogra-
phy. However, since the measurements are made *in-situ* with the sun heat-
ing the surface temperature daily, we propose a method that exploit such
characteristic by using Bayesian methods to track in time the evolution of
the materials’ surface temperature and estimate the different parameters on
each spectral band. Particular attention is given to the necessary hypotheses
for obtaining a generic filter but also to the computational time burden to
achieve on-the-fly estimations and parameter tracking.

1.2. Paper’s objectives

In this paper, emphasis will be made on the simultaneous estimation of state and parameters with the state-space discrete-time formulation of spatio-temporal systems. The numerical application focuses on infrared thermography measurements in a long-term and outdoor context where environmental and meteorological effects are corrected. The observed surfaces are thermally excited by the daily contribution of the Sun.

The results of our method are then compared to a MCMC algorithm from
the literature and to the Covariance Matrix Adaptation Evolution Strat-
egy (CMA-ES). To our knowledge, MCMC is the most recent and only
method that attempts to address the same problem as ours. As a conse-
quence, it can be considered as a benchmark regarding the existing littera-
ture. In comparison, CMA-ES is a well-known method in global optimization,
a field on which such generic problems may be solved and can tackle MCMC
formulations in some cases. The CMA-ES has been empirically success-
ful in hundreds of applications and is considered to be useful in particular
on non-convex, non-separable, ill-conditioned, multi-modal or noisy objec-
tive functions. Moreover, in contrast to most other evolutionary algorithms,
the CMA-ES is, from the user’s perspective, quasi-parameter-free, which en-
sure that its implementation results can not be worsened due to bad tuning
choices. It is a big advantage for benchmark and comparison purposes. Please
note that we applied the CMA-ES to get a global optimization method result
in comparison with MCMC and our developed method.

The state-space formulation will be introduced first. Then, the Marginal-
ized Particle Filter and Kriging steps will be presented. A direct applica-
tion of infrared thermography along with numerical results will be proposed.
Finally, outcomes on the newly presented Kriged Marginalized Particle Filter (KMPF) method for the simultaneous estimation of emissivity and temperature are addressed.

## 2. Spatio-temporal modeling

Spatio-temporal modeling is used in a wide range of varied domains, from environmental monitoring with climate change analysis (temperature trends tracking, precipitation patterns, etc.), air quality management, wildfire predictions, to public health with disease spread modeling, epidemiological studies or even urban planning including traffic management and land use planning. One of the main objectives of such modeling is to estimate the state of a spatial system that changes over time using a sequence of noisy measurements. Unfortunately, system parameters may also be unknown and thus require to be estimated as well. Therefore, it is important to use approaches that combine parameter and state estimation on the fly.

The most well known Bayesian state-estimation is the Kalman Filter (KF) [27]. By injecting the parameters into the state vector, the Extended Kalman Filter (EKF) allows a simultaneous estimation of states and parameters of the corresponding non-linear system but is generally biased or divergent [28]. Another method, the Particle Filter (PF), is therefore well suited for both parameter and state estimation, especially when non-linearity is involved and system parameter fluctuations are not negligible [29]. However, the computation of the weight of a particle is biased due to the induced nested dependencies between the state and the parameters. Furthermore, the computational burden explodes when the dimension of the state increases for a fair approximation [30]. As a consequence, a Marginalized Particle Filter (MPF) — also called Interacting Particle Filter — that combines encapsulated KF for state estimation with a PF for parameter estimation has been developed (see the note [31] and the article [32] for an introduction, the book [33] for more details). Decoupling the estimation of states and parameters helps to keep the system dimension within limits while taking advantage of the linearity of the dynamic part. A comparison of different estimation strategies is made in [34].

For spatio-temporal systems, Kriged Kalman Filter (KKF) have previously been used in the literature to estimate temporal varying fields based on observed data (see the book [35] and the paper [36] for a robust version). It has been used in various applications, such as large-scale data in pollution
and global climate monitoring [37], interpolation of medical images [38], or placement of sensors [39]. In those methods, the Kriging is performed to spatially interpolate sparse data. Moreover, the parameters of the system are known or estimated through an autoregressive model.

The proposed filter combines the two approaches: a MPF to estimate the parameters and the state of the underlying system as well as a Kriging step that will fit a spatial covariance model to the input image. This particular step aims at reducing the computational burden and takes advantage of the spatial characteristics of the signal to refine the particle selection.

## 2.1. State-Space Formulation

Let consider the linear evolution of the following spatio-temporal state sequence {$x_k, k \in \mathbb{N}$}.

$$x_k = A_k(\theta_k)x_{k-1} \quad (4)$$

where $x_k \in \mathbb{R}^n$ is the state vector taken at $n$ locations, $A_k(\theta_k) \in \mathbb{R}^n \times \mathbb{R}^n$ is the matrix that defines the linear dynamics of the system.

Similarly, the measurement equation is written as:

$$y_k = C_k(\theta_k)x_k \quad (5)$$

Where $y_k \in \mathbb{R}^n$ are the measured quantities, $C_k(\theta_k) \in \mathbb{R}^n \times \mathbb{R}^n$ and $\theta$ is a time-varying parameter vector. Finally, (4) and (5) form the so-called discrete-time state-space model of the spatio-temporal system.

The simplified system studied in this paper can now be written as a state-space system. Looking at (3), two main components of the signal can be extracted: a time-varying variable $T(s, t)$ and an assumed stationary part $\epsilon(s, t, \lambda_i)$. The system can therefore be discretized for a given band $\lambda_i$ as:

$$y_k(s) = C_k(\theta_k)(x_{s,k} + x_{d,k}), \quad \left\{ \begin{array}{l} x_{s,k} = \frac{\ln(\epsilon_k(s, \lambda_i))}{C}, \\ x_{d,k} = \frac{T_{\text{ref}}}{T_k(s)} \end{array} \right. \quad (6)$$

where $T_{\text{ref}}$ is a normalization factor for numeric purposes [5], $C = \frac{-C_2}{\lambda_i T_{\text{ref}}}$ is known in this application, $C_k(\theta_k) = C I$ and $y_k(s) = \ln(\frac{\gamma(T, s, t, \lambda_i)\lambda_i^5}{C_1})$.

The state-space vector to be estimated is therefore given by

$$x_k = x_{s,k} + x_{d,k} \quad (7)$$
where $x_{s,k}$ represents the stationary part of the state vector and $x_{d,k}$ the dynamic one.

Finally, the derived equation system is a time varying state-space equation system, has an ill-posed measurement equation with $2n$ unknowns for $n$ measurements and has a state vector that can be divided into a dynamic and a stationary component.

### 3. Marginalized Particle Filter

Proposed in [31], the Marginalized Particle Filter (MPF) has been applied to address mechanical system monitoring issues in [40] and further improved by [41] for noise robustness. MPF relies on two Bayesian methods: the particle and the Kalman filters. Whereas the Kalman filter is a key tool for state estimation based on the a priori knowledge of the underlying dynamical model, the particle filter is well suited for non-linear parameter estimation. The main idea is to run a particle filter where each particle embeds an instance of the model with its own set of parameter values and the corresponding Kalman filter. The idea of decoupling the parameter estimation through a PF and the state estimation through a KF is to enhance the stability of the algorithm. In fact, the relatively less expensive KF is used to estimate the linear state while the costly PF is used for the possibly non-linear parameters' estimation.

#### 3.1. Envelop Particle filter

Similarly to a standard PF, a set of $N_p \in \mathbb{N}$ independent particles enable the estimation of the posterior density of the parameter variables ($A_k(\theta_k)$, $C_k(\theta_k)$). Whereas different assumptions on the system evolution can be inferred through the time evolution of the different particles, only a random walk for the parameter evolution is assumed in this paper. Let $\xi_k^p$ be the parameter vector of a given particle $p \le N_p$ at the instant $k$, the random walk can be derived as:

$$ \xi_k^p = \xi_{k-1}^p + r_k^\xi \quad (8) $$

where $r_k^\xi \sim N(\delta\xi_k; \sigma_k^\xi)$ corresponds to a Gaussian blurring of $\xi_{k-1}^p$ with a shift of $\delta\xi_k$ and a spread of $\sigma_k^\xi$.
3.2. Embedded Kalman filters

To estimate the linear state variable, each particle embeds its own Kalman
filter with its associated parameters. As a consequence, once the parameters
of the numbered particle $p$ have been derived, the prediction of the Kalman
filter state $\mathbf{x}_{k|k}^p$ based on $\mathbf{x}_{k|k-1}^p$ can be done using the Kalman equations
[40]$^2$. The predicted estimate is therefore improved by using the current
measurement $\mathbf{y}_k$. Once the Kalman update has been applied to all the par-
ticles, their weight is updated using the likelihood of the innovation and the
Bayes rule updates the prior of the particles. Fig. 2 synthesizes the basic
principle of the MPF state and parameters estimation method.

<figure><img src="image_4.png" /><figcaption>Figure 2: Representation of the 3 main steps of the MPF method. The particles and their associated weights are represented with the blue circles and the circles sizes, respectively. The posterior probability density function of the parameters is represented in red.</figcaption></figure>

<footnote n="2">²Here, $\boldsymbol{x}_{i|j}$ denotes the random variable at $i$th instant provided the measurements up-to instant $j$</footnote>
# 4. Kriged Marginalized Particle Filter

The purpose of this section is to derive the equations that govern the presented KMPF.

## 4.1. Spatial Variability Modelization

In order to make full use of the spatio-temporal data, a spatial variability model is inferred onto the system by using the Kriging method. Kriging has been introduced in the geostatistical field to perform spatial interpolation on sparse data [42]. It provides the best linear unbiased prediction at unknown locations under suitable assumptions on the priors. In this paper, a kriging step is incorporated to the previously presented MPF to achieve two purposes. The first purpose is to bring a spatial covariance model to the system, and the second one is to reduce the computation time. Therefore, the stationary part of the state vector is obtained and modeled by using the Kriging equation, and detailed in Appendix A.

Kriging has already been used and studied along with Kalman filtering, known as the KKF [35, 36]. Most of the time, kriging is performed on sparse data to estimate the state vector at unknown locations [39, 37]. However, in this paper, measured data are full: kriging is applied differently to infer a spatial model onto the system. Moreover, the stationary $x_{s,k}$ component can not be extracted directly from the data since it is an ill-posed measurement equation problem (i.e. $x_{s,k}$ is mixed with the vector $x_{d,k}$ in the measurement equation). In other terms, it will be difficult to fit the trend at each Kalman step as it is done in a usual Kriging procedure because the fitting will then be highly dependent on the current Kalman solution leading to convergence toward a local solution. Alternatively, instead of estimating the trend using the data directly, the trend will be estimated through the particle filter to browse a broader solution space, and avoid local minima.

## 4.2. Focus on the dynamic part

Let now focus on the dynamic part of the state vector. Let $S \subset \mathbb{R}^2$ be the studied spatial domain, Let $\mathcal{S} \subset S$ be a subset of $S$ on which the algorithm will estimate the parameters and state vector. The estimates on $\mathcal{S}^* = S \setminus \mathcal{S}$ are given by the kriging step<fnref n="3" />.

<footnote n="3">³To accelerate the method, only a subset of pixels is used, the estimates on the full domain are then obtained with the Kriging. As such, $\text{card}(S) = m$ and $\text{card}(\mathcal{S}) = n$</footnote>
Let $\nu_k$ be the process noise $\nu_k \sim \mathcal{N}(0, \Sigma_{\nu,k})$. With (5) and (7), the measurement equation is

$$y_k = C_k(\theta_k)(x_{d,k} + x_{s,k}) + \nu_k \quad (9)$$

One can therefore retrieve a standard LTV system with $x_{d,k}$ state vector only. As proposed in [39], the stationary component can be considered as a noise term in (9). Denote $\mu_k = (\mu_k(s_1), \dots, \mu_k(s_n))^T$ the trend function of the stationary component (see Appendix A).

$$y_k - C_k(\theta_k)\mu_k = C_k(\theta_k)x_{d,k} + C_k(\theta_k)(x_{s,k} - \mu_k) + \nu_k \quad (10)$$

The measurement model equation becomes

$$y_k = C_k(\theta_k)(\mu_k + x_{d,k}) + e_k \quad (11)$$

$$\text{with} \quad \left\{ \begin{array}{l} e_k = C_k(\theta_k)(x_{s,k} - \mu_k) + \nu_k, \\ \Sigma_e = C_k(\theta_k)\Sigma_s C_k(\theta_k)^T + \Sigma_{\nu,k}. \end{array} \right. \quad (12)$$

where the covariance matrix $\Sigma_s$ of $x_{s,k}$ can be derived in Appendix B.

Please note that depending on the Kriging model used, the random variable $e_k$ may not follow a Gaussian process. When the variogram is chosen as a Gaussian model, then the classical Kalman results apply. However, in the case that the variogram is modeled differently, then the Kalman equations could be modified to consider non-Gaussian noises. In this study, the state vector is estimated linearly and the variance of such estimation is evaluated.

With $\eta_k \sim \mathcal{N}(0, \Sigma_\eta)$ the process noise, the resulting modified state-space model for the dynamic component is

$$\left\{ \begin{array}{l} x_{d,k} = A_k(\theta_k)x_{d,k-1} + \eta_k \\ y_k = C_k(\theta_k)(\mu_k + x_{d,k}) + e_k \end{array} \right. \quad (13)$$

### 4.3. Kalman equations

The Kalman equations can now be derived based on the modified state-space model (13). The objective here is to track the evolution of the dynamic component through time $x_{d,k}$ with unknown parameters $\theta_k$. However, if the dynamic component is tracked through time, the spatial expectation $\mu_k$ is unknown. It cannot be inferred from the data due to the ill-posed characteristic of the considered measurement equation. The spatial covariance matrix
$\Sigma_s$ is assumed to be known and constant through time, obtained for example from an initial guess at instant $k = 0$.

The numbered particle $p$ can therefore be defined through its vector as $\xi_k^p = (\theta_k^p, \mu_k^p)$. As a consequence, the state-space model for a given particle, considering dynamic and observation noises is

$$ \left\{ \begin{aligned} \mathbf{x}_{d,k}^p &= \mathbf{A}_k(\theta_k^p)\mathbf{x}_{d,k-1}^p + \boldsymbol{\eta}_k \\ \mathbf{y}_k^p &= \mathbf{C}_k(\theta_k^p)(\boldsymbol{\mu}_k^p + \mathbf{x}_{d,k}^p) + \mathbf{e}_k^p \end{aligned} \right. \qquad (14) $$

The prediction and correction steps are given in Appendix C.

### 4.3.1. Particle update

From [40], the particle approximation $\hat{\xi}_k$ can be derived. The likelihood criterion $p(\mathbf{y}_k | \xi_k^p)$ is given by

$$ p(\mathbf{y}_k | \xi_k^p) = \mathcal{N}(\mathbf{C}_k(\theta_k^p)(\boldsymbol{\mu}_k^p + \mathbf{x}_{d,k}^p), \mathbf{S}_k^p) \qquad (15) $$

Similarly to a classical PF algorithm, the particle approximation is obtained with

$$ \hat{\xi}_k = \sum_{p=1}^{N_p} w_k^p \xi_k^p, \quad \text{with } w_k^p = \frac{w_{k-1}^p p(\mathbf{y}_k | \xi_k^p)}{\sum_{p=1}^{N_p} w_{k-1}^p p(\mathbf{y}_k | \xi_k^p)} \qquad (16) $$

Resampling can be done according to [41]. This algorithm is performed on each pixel $\mathbf{s}$ of $\mathcal{S}$. Kriging equation (A.10) is then used to retrieve the estimates on $\mathcal{S}^*$. To retrieve $\mathbf{x}_{d,k}^*$ on $\mathcal{S}^*$, the measurement equation is then applied with the kriged estimates $\mathbf{x}_{s,k}^*$ on $\mathcal{S}^*$ and the measurements $\mathbf{y}_k^*$.

This algorithm has the property of reducing the number of variables and using the kriging to speed up the computation. In the next section, a numerical application of the presented method for the simultaneous estimation of temperature and emissivity in infrared imaging is performed.

## 5. Numerical Application and Benchmark

In order to compare our proposed method with two other referenced methods from literature, a numerical test case has been created thanks to a dedicated radiosity simulation software tool [43].
## 5.1. *Description of the Numerical Test Case*

A target made of four distinct materials with known radiative properties has been modeled, as shown in Fig. 3. The frame's target (black coating) is kept at constant temperature whereas the four other materials (alumina, concrete, wood and gold) are assumed to have a varying controlled temperature inside the squares. The emissivity profiles are taken from a database of real materials and the temperature evolution is taken from outside ground thermocouple temperature measurements during January 2017 [6]. Four different spectral bands have been used for the simulation: [1.5; 3]µm, [3; 5]µm, [8; 10]µm, [10; 12]µm. The grey-body hypothesis is performed: the emissivity is considered constant, equals to the mean value on that band. For the alumina and gold, data from [1.5; 2]µm were missing, so the mean value of the first band was computed on [2; 3]µm instead. The simulated data represent noisy data recorded by an infrared camera.

<figure><img src="image_5.png" /><figcaption>Figure 3: Use-case illustrations. From top-left to bottom: rendering of the 3D simulation scene in the visible, result example of the target's simulated radiosity and emissivity profiles of the four different materials.</figcaption></figure>

## 5.2. *Evaluated methods*

The implementations are made in R for the KMPF method to make use of the DiceKriging library [44] and MATLAB® for the two others methods.
### 5.2.1. References methods from literature

Two methods (MCMC and CMA-ES) have been implemented and tested conjointly on the simulated data for which preliminary results were shown in [45, 46].

The MCMC metaheuristic method has already been used in the literature for temperature and emissivity estimation [47, 10]. A more detailed description of this implementation can be found in [48].

The second method is the CMA-ES. The CMA-ES [49] is a global optimization algorithm, necessary in this case where the solution space makes other optimization methods converge towards a local solution (see [48]). Even though this specific method has never been used for this purpose to our knowledge, optimization formulations for the TES problem are common [5, 9]. Furthermore, this method is parameters-free and one of the references in the global optimization field [50, 51]. The optimization problem is defined by considering that locally the emissivity and temperature are homogeneous. The details of the objective function is given in Appendix D.

Regarding the initialization of the algorithms, broad priors have been taken. For the MCMC: $N(0.5, \sigma_e^2 = 0.5)$ for the emissivity and $N(T_{\text{air}}, \sigma_T^2 = 10)$ for the temperature where $T_{\text{air}}$ is the air / environment temperature set in the simulation software. For the CMA-ES one, the initial point was set to $T_{\text{air}}$ for the temperature, giving a prior for the emissivity using the inverse model at this temperature. Upper and lower bounds are set so that the emissivity lies in $[0; 1]$ and the temperature in $[T_{\text{air}} - 2\sigma_T; T_{\text{air}} + 2\sigma_T]$.

The evaluation dataset remains consistent across methods, but MCMC results are reported for only one pixel per material at a given time due to its high computational cost and single-pixel inference approach. In contrast, CMA-ES processes the entire image for each timestamp separately.

### 5.2.2. KMPF Formulation for the application

The dynamic evolution of the temperature is assumed to be homogeneous on the image, linear and unknown: $A_k = a_k I$. Similarly, $\Sigma_{\eta}$ is a diagonal matrix and $C_k = \frac{-C_2}{\lambda_i T_{\text{ref}}} I$. Please note that this model could be extended by adding relationship between the selected pixels and therefore modifying $A_k$ and $C_k$. This state-space model is updated over time. To ensure the convergence of the Kalman filter, the algorithm is started at the beginning of the day (when $0 < a_t < 1$). The algorithm is applied for single spectral band images. The parameters to be estimated are therefore $a_k$ and $\mu$, the trend of the emissivity. KMPF processes the entire image across all timestamps.
Since it is assumed that locally, at thermal equilibrium, the temperature in a given pixel's neighborhood is the same (local homogeneity), the particle's selection criterion as defined in (16) can be refined by comparing the obtained Kriged values to the actual measurements around the pixels of interest. This criterion is suitable in this specific case due to its physical nature. Finally, an estimation of the confidence on the temperature and emissivity estimates can be deduced.

Regarding the initialization of the algorithm, the same broad priors have been taken for KMPF and MCMC while all three algorithms share the same basic assumptions for benchmark comparison.

### 5.3. Benchmark results

To ensure a fair comparison, the same timestamp was used for all methods, and results were compared per material. Results are presented in Table 1, Table 2, Table 3 which share the same true values, facilitating method comparison.

#### 5.3.1. MCMC

The simultaneous estimation of temperature and emissivity is done at one given instant, for a given pixel. The algorithm stops after ten thousand Monte Carlo realizations (which represents approximately 20 min). Fig. 4 and Fig. 5 show the results obtained for the alumina material.

The temperature of 272.16 K is estimated at 277.19 K with a RMSE of $3.42 \times 10^{-2}$ K, the best result of the four materials (see Table 1).

The MCMC method retrieves the original profiles of emissivity (with the gray body hypothesis) as well as a convergence towards the original temperature. Please note that this method has an important computational burden and that it is sensible to the number of Monte Carlo realizations, which has to be fixed in advance. A larger number of Monte Carlo realizations would eventually lead to better results.

#### 5.3.2. CMA-ES

The result is computed on an image at a given instant and spectral band. A median filter is applied to the final image to smooth the result, shown in Fig. 6.

The CMA-ES method enlightens the different materials through the emissivity retrieval. It can be seen that a relatively small difference on small emissivity values implies a relative important error on temperature, which is
MCMC: result example for alumina

<figure><img src="image_6.png" /><figcaption>Figure 4: MCMC: result example for alumina material emissivity estimation for a given pixel at a given instant. The true value (dashed) is compared to the estimation (solid).</figcaption></figure>

MCMC : result example for alumina (μ=272.19, σ=4.59)

<figure><img src="image_7.png" /><figcaption>Figure 5: MCMC: result example for alumina material temperature estimation for a given pixel at a given instant. The true value (dashed) is compared to the estimation (solid).</figcaption></figure>
<table><caption>Table 1: MCMC: results obtained for the simultaneous estimation of temperature and emissivity for the four different materials at the four different bands. ewline * Spatial average on January 18th at 7:59 PM.</caption><thead><tr><th>Parameter</th><th>Ground Truth</th><th>Estimation</th><th>RMSE</th><th>σ</th></tr></thead><tbody><tr><td>Temperature (K)</td><td>272.16</td><td>271.69</td><td>4.68e-01</td><td>3.85e+00</td></tr><tr><td>ε<sub>[1.5;3]μm</sub></td><td>0.09</td><td>0.10</td><td>1.01e-02</td><td>2.26e-01</td></tr><tr><td>ε<sub>[3;5]μm</sub></td><td>0.08</td><td>0.09</td><td>6.33e-03</td><td>2.25e-01</td></tr><tr><td>ε<sub>[8;10]μm</sub></td><td>0.12</td><td>0.12</td><td>4.80e-03</td><td>2.25e-01</td></tr><tr><td>ε<sub>[10;12]μm</sub></td><td>0.11</td><td>0.11</td><td>4.32e-03</td><td>2.25e-01</td></tr><tr><td>ε on all bands</td><td>-</td><td>-</td><td>1.35e-02</td><td>-</td></tr><tr><td>Temperature (K)</td><td>272.16</td><td>271.76</td><td>4.03e-01</td><td>5.09e+00</td></tr><tr><td>ε<sub>[1.5;3]μm</sub></td><td>0.64</td><td>0.69</td><td>5.46e-02</td><td>2.92e-01</td></tr><tr><td>ε<sub>[3;5]μm</sub></td><td>0.85</td><td>0.89</td><td>3.96e-02</td><td>2.60e-01</td></tr><tr><td>ε<sub>[8;10]μm</sub></td><td>0.96</td><td>0.98</td><td>1.24e-02</td><td>2.28e-01</td></tr><tr><td>ε<sub>[10;12]μm</sub></td><td>0.92</td><td>0.93</td><td>8.06e-03</td><td>2.26e-01</td></tr><tr><td>ε on all bands</td><td>-</td><td>-</td><td>6.91e-02</td><td>-</td></tr><tr><td>Temperature (K)</td><td>272.16</td><td>272.19</td><td>3.42e-02</td><td>4.59e+00</td></tr><tr><td>ε<sub>[1.5;3]μm</sub></td><td>0.30</td><td>0.30</td><td>1.48e-03</td><td>2.39e-01</td></tr><tr><td>ε<sub>[3;5]μm</sub></td><td>0.51</td><td>0.50</td><td>7.58e-03</td><td>2.38e-01</td></tr><tr><td>ε<sub>[8;10]μm</sub></td><td>0.95</td><td>0.94</td><td>1.20e-02</td><td>2.35e-01</td></tr><tr><td>ε<sub>[10;12]μm</sub></td><td>0.70</td><td>0.68</td><td>1.58e-02</td><td>2.29e-01</td></tr><tr><td>ε on all bands</td><td>-</td><td>-</td><td>2.13e-02</td><td>-</td></tr><tr><td>Temperature (K)</td><td>272.16</td><td>272.21</td><td>5.14e-02</td><td>5.03e+00</td></tr><tr><td>ε<sub>[1.5;3]μm</sub></td><td>0.70</td><td>0.73</td><td>2.39e-02</td><td>3.09e-01</td></tr><tr><td>ε<sub>[3;5]μm</sub></td><td>0.90</td><td>0.90</td><td>3.93e-03</td><td>2.65e-01</td></tr><tr><td>ε<sub>[8;10]μm</sub></td><td>0.87</td><td>0.85</td><td>1.52e-02</td><td>2.27e-01</td></tr><tr><td>ε<sub>[10;12]μm</sub></td><td>0.95</td><td>0.93</td><td>1.80e-02</td><td>2.27e-01</td></tr><tr><td>ε on all bands</td><td>-</td><td>-</td><td>3.38e-02</td><td>-</td></tr></tbody></table>
<figure><img src="image_8.png" /><figcaption>Figure 6: CMA-ES: simultaneous estimation of temperature and emissivity at a given instant and spectral band ([10; 12]µm). From left to right: true value, estimated value and absolute difference. From up to bottom: emissivity, temperature.</figcaption></figure>

inherent to the physical equation. Although the results are promising, the computational burden is still high, particularly for long-term and on-the-fly applications (≈ 0.08 s for one pixel at one instant). The numerical results of the estimations are presented in Table 2.

### 5.3.3. KMPF

One of the main assumptions done in this application case is that the dynamic evolution of the temperature governed by the $a_k$ parameter is the same over the entire image. However, in this simulation, the black coated frame remains at a constant temperature. To overcome this, the kriging reference points (belonging to $\mathcal{S}$) are forced to be in the 4 materials only and not in the frame. Fig. 7 and Fig. 8 show the results of the estimation for one spectral band and over the entire image and time domains. Table 3 shows the numerical results on the 4 bands at a given time.

As for the CMA-ES method, the KMPF one reveals the different materials through the emissivity retrieval (Fig. 7). By comparison with the CMA-ES, the error in temperature is more homogeneous between the materials due to the use of a spatial model. Fig. 8 shows that the standard deviation of
<table><caption>Table 2: CMA-ES: results obtained for the simultaneous estimation of temperature and emissivity for the four different materials at the four different bands.</caption><thead><tr><th>Parameter</th><th>Ground Truth</th><th>Estimation*</th><th>Min. RMSE*</th><th>Mean RMSE*</th><th>Max. RMSE*</th></tr></thead><tbody><tr><td>GOLD Temperature (K)</td><td>272.16</td><td>269.02</td><td>1.15e+00</td><td>3.15e+00</td><td>4.81e+00</td></tr><tr><td>ε[1.5;3]μm</td><td>0.09</td><td>0.13</td><td>2.61e-02</td><td>3.80e-02</td><td>1.89e-01</td></tr><tr><td>ε[3;5]μm</td><td>0.08</td><td>0.10</td><td>4.74e-03</td><td>1.47e-02</td><td>9.88e-02</td></tr><tr><td>ε[8;10]μm</td><td>0.12</td><td>0.12</td><td>2.31e-03</td><td>6.98e-03</td><td>5.10e-02</td></tr><tr><td>ε[10;12]μm</td><td>0.11</td><td>0.11</td><td>2.92e-04</td><td>4.33e-03</td><td>4.53e-02</td></tr><tr><td>ε on all bands</td><td>-</td><td>-</td><td>2.92e-04</td><td>6.99e-02</td><td>2.87e-01</td></tr><tr><td>WOOD Temperature (K)</td><td>272.16</td><td>268.98</td><td>1.04e+00</td><td>3.18e+00</td><td>3.96e+00</td></tr><tr><td>ε[1.5;3]μm</td><td>0.64</td><td>0.85</td><td>1.57e-01</td><td>2.16e-01</td><td>2.87e-01</td></tr><tr><td>ε[3;5]μm</td><td>0.85</td><td>0.94</td><td>5.12e-02</td><td>8.18e-02</td><td>1.08e-01</td></tr><tr><td>ε[8;10]μm</td><td>0.96</td><td>0.98</td><td>7.15e-03</td><td>1.76e-02</td><td>2.64e-02</td></tr><tr><td>ε[10;12]μm</td><td>0.92</td><td>0.98</td><td>4.21e-02</td><td>5.55e-02</td><td>6.81e-02</td></tr><tr><td>ε on all bands</td><td>-</td><td>-</td><td>2.92e-04</td><td>6.99e-02</td><td>2.87e-01</td></tr><tr><td>ALUMINA Temperature (K)</td><td>272.16</td><td>268.72</td><td>2.25e+00</td><td>3.44e+00</td><td>4.54e+00</td></tr><tr><td>ε[1.5;3]μm</td><td>0.30</td><td>0.46</td><td>1.30e-01</td><td>1.57e-01</td><td>2.14e-01</td></tr><tr><td>ε[3;5]μm</td><td>0.51</td><td>0.63</td><td>9.02e-02</td><td>1.20e-01</td><td>1.74e-01</td></tr><tr><td>ε[8;10]μm</td><td>0.95</td><td>0.98</td><td>1.37e-02</td><td>2.81e-02</td><td>3.81e-02</td></tr><tr><td>ε[10;12]μm</td><td>0.70</td><td>0.74</td><td>2.95e-02</td><td>4.58e-02</td><td>6.26e-02</td></tr><tr><td>ε on all bands</td><td>-</td><td>-</td><td>2.92e-04</td><td>6.99e-02</td><td>2.87e-01</td></tr><tr><td>CONCRETE Temperature (K)</td><td>272.16</td><td>270.35</td><td>3.46e-01</td><td>1.81e+00</td><td>2.35e+00</td></tr><tr><td>ε[1.5;3]μm</td><td>0.70</td><td>0.88</td><td>1.40e-01</td><td>1.78e-01</td><td>2.32e-01</td></tr><tr><td>ε[3;5]μm</td><td>0.90</td><td>0.94</td><td>1.67e-02</td><td>4.51e-02</td><td>7.15e-02</td></tr><tr><td>ε[8;10]μm</td><td>0.87</td><td>0.95</td><td>5.20e-02</td><td>7.73e-02</td><td>1.01e-01</td></tr><tr><td>ε[10;12]μm</td><td>0.95</td><td>0.98</td><td>1.59e-02</td><td>3.08e-02</td><td>3.97e-02</td></tr><tr><td>ε on all bands</td><td>-</td><td>-</td><td>2.92e-04</td><td>6.99e-02</td><td>2.87e-01</td></tr></tbody></table>

<footnote>* Spatial average on January 18th at 7:59 PM.</footnote>
<figure><img src="image_9.png" /><figcaption>Figure 7: KMPF: qualitative example of the simultaneous estimation of temperature and emissivity at a given spectral band ([10; 12]µm). From left to right: true value, estimated value and absolute difference. From up to bottom: emissivity, temperature.</figcaption></figure>
<table><caption>Table 3: KMPF: results obtained for the simultaneous estimation of temperature and emissivity for the four different materials at the four different bands.</caption><thead><tr><th>Parameter</th><th>Ground Truth</th><th>Estimation*</th><th>Min. RMSE*</th><th>Mean RMSE*</th><th>Max. RMSE*</th></tr></thead><tbody><tr><td>Temperature (K)</td><td>272.16</td><td>271.00</td><td>2.25e-03</td><td>1.04e+00</td><td>6.37e+00</td></tr><tr><td>ε[1.5;3]μm</td><td>0.09</td><td>0.10</td><td>2.32e-05</td><td>9.99e-03</td><td>6.46e-02</td></tr><tr><td>ε[3;5]μm</td><td>0.08</td><td>0.09</td><td>1.44e-06</td><td>1.37e-02</td><td>1.39e-01</td></tr><tr><td>ε[8;10]μm</td><td>0.12</td><td>0.12</td><td>4.77e-05</td><td>1.66e-02</td><td>1.84e-01</td></tr><tr><td>ε[10;12]μm</td><td>0.11</td><td>0.11</td><td>1.25e-05</td><td>9.11e-03</td><td>8.25e-02</td></tr><tr><td>ε on all bands</td><td>-</td><td>-</td><td>1.44e-06</td><td>1.23e-02</td><td>1.84e-01</td></tr><tr><td>Temperature (K)</td><td>272.16</td><td>273.69</td><td>2.53e-03</td><td>1.54e+00</td><td>5.43e+00</td></tr><tr><td>ε[1.5;3]μm</td><td>0.64</td><td>0.66</td><td>1.95e-04</td><td>2.23e-02</td><td>7.97e-02</td></tr><tr><td>ε[3;5]μm</td><td>0.85</td><td>0.80</td><td>9.99e-05</td><td>5.54e-02</td><td>1.98e-01</td></tr><tr><td>ε[8;10]μm</td><td>0.96</td><td>0.90</td><td>2.84e-04</td><td>6.13e-02</td><td>1.44e-01</td></tr><tr><td>ε[10;12]μm</td><td>0.92</td><td>0.90</td><td>1.86e-04</td><td>2.97e-02</td><td>1.32e-01</td></tr><tr><td>ε on all bands</td><td>-</td><td>-</td><td>9.99e-05</td><td>4.22e-02</td><td>1.98e-01</td></tr><tr><td>Temperature (K)</td><td>272.16</td><td>275.02</td><td>1.48e-02</td><td>3.08e+00</td><td>8.82e+00</td></tr><tr><td>ε[1.5;3]μm</td><td>0.30</td><td>0.31</td><td>1.25e-05</td><td>1.28e-02</td><td>5.55e-02</td></tr><tr><td>ε[3;5]μm</td><td>0.51</td><td>0.48</td><td>3.81e-03</td><td>2.95e-02</td><td>7.98e-02</td></tr><tr><td>ε[8;10]μm</td><td>0.95</td><td>0.89</td><td>2.81e-04</td><td>6.20e-02</td><td>1.74e-01</td></tr><tr><td>ε[10;12]μm</td><td>0.70</td><td>0.68</td><td>1.76e-04</td><td>2.30e-02</td><td>6.08e-02</td></tr><tr><td>ε on all bands</td><td>-</td><td>-</td><td>1.25e-05</td><td>3.18e-02</td><td>1.74e-01</td></tr><tr><td>Temperature (K)</td><td>272.16</td><td>273.84</td><td>1.66e-03</td><td>1.87e+00</td><td>1.41e+01</td></tr><tr><td>ε[1.5;3]μm</td><td>0.70</td><td>0.73</td><td>5.35e-05</td><td>2.85e-02</td><td>6.47e-02</td></tr><tr><td>ε[3;5]μm</td><td>0.90</td><td>0.82</td><td>5.65e-04</td><td>7.93e-02</td><td>2.38e-01</td></tr><tr><td>ε[8;10]μm</td><td>0.87</td><td>0.80</td><td>6.35e-03</td><td>6.47e-02</td><td>1.66e-01</td></tr><tr><td>ε[10;12]μm</td><td>0.95</td><td>0.92</td><td>3.44e-05</td><td>3.24e-02</td><td>2.23e-01</td></tr><tr><td>ε on all bands</td><td>-</td><td>-</td><td>3.44e-05</td><td>5.12e-02</td><td>2.38e-01</td></tr></tbody></table>

<footnote>* Spatial average on January 18th at 7:59 PM.</footnote>
<figure><img src="image_10.png" /><figcaption>Figure 8: KMPF: temperature through time for a given pixel, input data (decimated from true value) in dashed orange, true value in dot-dashed blue and estimation in solid blue. A cursor shows the time at which the other results of this paper are presented.</figcaption></figure>

the underlying Kalman Filter is increasing during the afternoon and night. This is due to the fact that this part corresponds to $A_k(\theta_k) > 1$ and thus a divergent recursive sequence.

Furthermore, Fig. 9 shows promising results for emissivity estimation, here with a focus for alumina material at a given time step. Numerical results are given in Table 3 for which the values are obtained by doing a spatial average on the area of each material. As a consequence, this average take into account the pixels that have been interpolated by the Kriging step.

Finally, thanks to the kriging that allows the tracking of a subset of pixels, the time is approximately of 0.06 ms per pixel at a given instant, so approximately 1000 times faster than the CMA-ES algorithm.

# 6. Discussion

Three different methods have been tested on this Temperature and Emissivity Separation (TES) problem for thermal measurements by infrared thermography. It is important to note that those methods do not address the same condition measurements: the two first ones range from pixel level to im-
KMPF: result example for alumine

<figure><img src="image_11.png" /><figcaption>Figure 9: KMPF: result example for alumina material emissivity estimation for a given pixel of $S$ on January 18$^{th}$ at 7:59 PM. The true value (dashed) is compared to the estimation (solid).</figcaption></figure>

age level at a given time and the last one addresses full image time evolution (dynamic scenes). All methods are studied and built to calculate emissivity and temperature in different spectral bandwidths of various lengths. The first one, MCMC, is a pixel level method for input data at a given time, whereas the second one, the CMA-ES method, is concerned by a neighborhood of pixels up to image spatial scale to get an estimate. The third and last one, the KMPF method, addresses full frame image evolution over time. By its construction KMPF will be more adapted for long term measurements on dynamic scenes, whereas MCMC and CMA-ES are more suited for static ones. Performances of KMPF rely on the temperature kinetic evolution. Moreover, performances are also affected by the kriging part of the method that interpolates the emissivity over the image and also acts as a drift corrector in interaction with the filter. As Kalman and particle filters are tracking algorithms, they need to adapt to a fast changing parameter evolution as illustrated in Fig. 8.

Fig. 10 and Fig. 11 show a comparison of the performance of the presented methods regarding the estimation bias for both the emissivity and the temperature, as well as a comparison of the required computation time for achieving those results.

In Fig. 10, for all studied methods, estimation bias is presented for the
considered spectral band. Band 1 and band 2 yield to lower performance in comparison to band 3 and band 4 located in the long wave infrared spectrum. The more emissive and regular the material is in the considered spectral band, the lower the estimation bias is. Globally all methods yield to equivalent performance for the emissivity estimation.

In Fig. 11, left graph, for a given time step we observe less dispersion for temperature estimation by MCMC than by KMPF and finally CMA-ES. Considering their performance for each material, it can be inferred that KMPF yields to similar estimates compared to the other methods.

In Fig. 11, right graph, the developed method KMPF appears to be the only one that is suited for the live estimation of temperature and emissivity through IR measurements, thanks to its relatively low computation time regarding the others methods while providing comparable results as stated above.

Finally, the quality of the method has been investigated in terms of bias for the emissivity, but spread is also another matter of evaluation that is now discussed. A comparison of the obtained emissivity for the different materials is shown in Fig. 12, where the estimation of the emissivity for the compared methods is performed. Notice that different confidence information are available for the 3 methods. For MCMC, a tolerance interval is defined, which is decreasing with the number of Monte Carlo realizations at the expense of computational burden. For others, the confidence interval is provided as a by-product of the statistical properties of the algorithm. A large spread for MCMC shows that the estimation could be eventually improved by allowing more time to the algorithm, whereas no improvement in the solution could be achieved from the implemented CMA-ES algorithm. Then, the green whiskers in Fig. 12 show the estimated standard deviation of the *a priori* emissivity distribution of the KMPF. For all KMPF estimates, the true value belongs to the corresponding statistical confidence intervals, unlike CMA-ES. This is also true for MCMC at the expense of very large tolerance regions.

It should also be noted that none of the algorithms have been optimized in terms of computing speed, even though a GPU-acceleration could have been considered. Moreover, the quality of the estimation may vary depending on the considered time instant and chosen band for the different estimations methods. However, the computational time burden gives an interesting indicator for the feasibility of an *in-situ* and live estimation of the temperature and emissivity.
<figure><img src="image_12.png" /><figcaption>Figure 10: Comparison of the emissivity normalized root mean square error for the different methods. The normalized root mean square error is defined as the ratio of the root mean square error to its true value. The values correspond to the numerical results presented in the previous tables.</figcaption></figure>
Implemented methods comparison - temperature and time burden

<figure><img src="image_13.png" /><figcaption>Figure 11: On the left, comparison of the temperature root mean square error for the different methods, corresponding to the numerical results presented in the previous tables. On the right, the computation time brought to one pixel at a given band and instant in log-scale.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 12: Comparison of the emissivity value obtained with the different methods. The blue whiskers represent the tolerance interval for the MCMC and CMA-ES methods while the green ones represent the estimated standard deviation for the KMPF.</figcaption></figure>

## 7. Conclusion

In this paper, a new method that allows the simultaneous estimation of states and parameters in spatio-temporal based systems called Kriged Marginalized Particle Filter is proposed. The general formulation of the problem is firstly introduced, then the equations are derived for the temperature and emissivity separation problem in infrared thermography. A numerical application and a comparison of the proposed method with state-of-the-art ones are finally presented. If the MCMC and CMA-ES methods offer even results, the required computation time and energy are an obstacle for on-the-fly or long-term data analysis (see Fig. 11). In contrast, the KMPF and the inference of a spatial model in the equation system, enable to reduce the computational time burden. To be applied, our method relies on the local homogeneity of the temperature, as well as having a dynamic excitation of the temperature (i.e. a set of images evolving through time). Nevertheless, the simultaneous estimation of the temperature and emissivity of the observed object on simulated data yields promising results toward absolute temperature estimation by multispectral infrared thermography. By its generic aspect, this filter could be extended to others applications, but
further work on robustness is needed for application on real data. Finally, potential future directions include enhancing the physical modeling to derive more accurate state-space equations and fine-tuning the variance evolution, such as incorporating the trend slope of parameter changes into the analysis.

Acknowledgment

The authors would like to thank the French Region Bretagne and BRIGHTER project for their support. BRIGHTER has received funding from the Chips Joint Undertaking (JU) under grant agreement No 101096985. The JU receives support from the European Union's Horizon Europe research and innovation program and France, Belgium, Portugal, Spain, Turkey.
## Appendix A. Derivation of the stationary part using Kriging equations

Let $s \in \mathcal{S}$ be a spatial point with its measurement $Z(s)$ and $\tilde{s} = s + h \in \mathcal{S}$. To model the spatial variability of the observed quantities and thus infer a spatial model into the equation system, the idea is to fit a spatial covariance function on the measurement. Such step is accomplished by modelling the variogram, a function that describes the degree of spatial dependence of the spatial random field $(Z(s) - Z(s + h))^2$ [52]. The empirical variogram is defined as:

$$2\gamma(h) = \frac{1}{N(h)} \sum_{i=1}^{N(h)} (Z(s_i) - Z(s_i + h))^2 \quad (A.1)$$

where $N(h)$ represents the number of points separated by a vector $h$. Based on (A.1), the main assumption of such modelling is to suppose that the variogram depends on the modulus and direction of the vector $h$ and not the location $s$. Once the experimental variogram has been computed from the measurements points of $\mathcal{S}$, a theoretical model is fitted to it, as shown as an illustrating example in Fig. A.13.

<figure><img src="image_15.png" /><figcaption>Figure A.13: Variogram fitting example: the dots represent the experimental variogram values and the lines represent two different variogram theoretical models.</figcaption></figure>

The fitted variogram is then used in the following to give a spatial covariance model to the stationary component of the system.
As mentioned above, the state vector in this study is assumed to be decomposable into a stationary component and a dynamic one. Let focus on the stationary one first and assume that $x_{s,k}(\mathbf{s})$ can be decomposed into a deterministic trend function $\mu_k : \mathbf{s} \in \mathcal{S} \to \mu_k(\mathbf{s}) \in \mathbb{R}$ and a square-integrable process $Z(\mathbf{s})$ with known spatial covariance model (learned from the variogram) and zero mean:

$$x_{s,k}(\mathbf{s}) = \mu_k(\mathbf{s}) + Z(\mathbf{s}) \qquad (A.2)$$

$$\mathbb{E}[x_{s,k}(\mathbf{s})] = \mathbb{E}[\mu_k(\mathbf{s})] + \mathbb{E}[Z(\mathbf{s})] = \mu_k(\mathbf{s}) \qquad (A.3)$$

$$\gamma_Z(\mathbf{h}) = \frac{1}{2} \mathbb{E} [ |Z(\mathbf{s} + \mathbf{h}) - Z(\mathbf{s})|^2 ] \qquad (A.4)$$

The objective of the universal kriging (UK) is to estimate both the best linear prediction of $x_{s,k}(\mathbf{s})$ on $\mathcal{S}$ and the trend function $\mu_k(\mathbf{s})$ with a non-constant spatial mean. In this case, the trend function is of the form $\mu_k(\mathbf{s}) = \sum_{i=1}^q \beta_i f_i(\mathbf{s})$ where $\beta = (\beta_i)_{0 \le i \le q}$ is the unknown drift coefficient vector and $f_i : \mathbb{R}^2 \to \mathbb{R}$ are given basis function:

$$x_{s,k}(\mathbf{s}) = \sum_{i=0}^{q} \beta_i f_i(\mathbf{s}) + Z(\mathbf{s}) \qquad (A.5)$$

By denoting $\mathbf{F} = (f_i(\mathbf{s}_j))_{0 \le i \le q, 1 \le j \le n}$, $\mathbf{x}_{s,k}$, the random vector of the stationary component on $\mathcal{S}$ can be written as:

$$\mathbf{x}_{s,k} = \begin{pmatrix} x_{s,k}(\mathbf{s}_1) \\ \vdots \\ x_{s,k}(\mathbf{s}_n) \end{pmatrix} \qquad (A.6)$$

$$= \begin{pmatrix} f_0(\mathbf{s}_1) & \dots & f_q(\mathbf{s}_1) \\ \vdots & \ddots & \vdots \\ f_0(\mathbf{s}_n) & \dots & f_q(\mathbf{s}_n) \end{pmatrix} \begin{pmatrix} \beta_0 \\ \vdots \\ \beta_q \end{pmatrix} + \begin{pmatrix} Z(\mathbf{s}_1) \\ \vdots \\ Z(\mathbf{s}_n) \end{pmatrix} \qquad (A.7)$$

$$= \mathbf{F}\beta + \mathbf{Z} \qquad (A.8)$$

The UK prediction $x_{s,k}^*(\mathbf{s}_0^*)$ of the value of $x_{s,k}(\mathbf{s})$ at the location $\mathbf{s}_0^* \in \mathcal{S}^*$ is given by the linear sum:

$$x_{s,k}^*(\mathbf{s}_0^*) = \sum_{i=1}^{n} \omega_i x_{s,k}(\mathbf{s}_i) = \omega^T (\mathbf{F}\beta + \mathbf{Z}) \qquad (A.9)$$
where $\omega_i$ corresponds to the kriging interpolation weights and $\omega = (\omega_1, \dots, \omega_n)^T$.

The kriging equations give the values of the weight vector $\omega$ — after satisfying unbiasedness condition and minimal observation variance of the predictions [53, 35] — by using the input measurement data $\mathbf{x}_{s,k}$, the basis function and after fitting the variogram theoretical model and trend to the data:

$$ \forall \mathbf{s}^* \in \mathcal{S}^*, \quad x_{s,k}^*(\mathbf{s}^*) = \qquad (A.10) $$

$$ [\gamma_{Z,s^*} - \mathbf{F}(\mathbf{F}^T \Gamma_Z^{-1} \mathbf{F})^{-1}(\mathbf{F}^T \Gamma_Z^{-1} \gamma_{Z,s^*} - \mathbf{f}^*)]^T \Gamma_Z^{-1} \mathbf{x}_{s,k} \qquad (A.11) $$

where

$$ \left\{ \begin{array}{l} \Gamma_Z = (\gamma_Z(\mathbf{s}_i - \mathbf{s}_j)_{1\le i,j \le n}, \\ \mathbf{f}^* = (f_i(\mathbf{s}^*))_{0\le i \le q}, \\ \gamma_{Z,s^*} = (\gamma_Z(\mathbf{s}_1 - \mathbf{s}^*), \gamma_Z(\mathbf{s}_2 - \mathbf{s}^*), \dots, \gamma_Z(\mathbf{s}_n - \mathbf{s}^*)) \end{array} \right. \qquad (A.12) $$

## Appendix B. Covariance Matrix Derivation On $\mathcal{S}$

Let compute the covariance matrix $\Sigma_s$ of the vector $\mathbf{x}_{s,k}$ denoted $\mathbf{x}_s$ for convenience. For a given instant $k$:

$$ \Sigma_s = (\operatorname{cov}(x_s(\mathbf{s}_i), x_s(\mathbf{s}_j)))_{1\le i,j \le n} \qquad (B.1) $$

To compute such covariance, one need to link the variogram of $Z$ to its covariance ($Z$ is assumed second-order stationary with zero mean):

$$ \gamma_Z(\mathbf{h}) = \frac{1}{2} \mathbb{E} \left[ (Z(\mathbf{s} + \mathbf{h}) - Z(\mathbf{s}))^2 \right] \qquad (B.2) $$

$$ = \frac{1}{2} \mathbb{E} [Z(\mathbf{s} + \mathbf{h})^2] + \frac{1}{2} \mathbb{E} [Z(\mathbf{s})^2] - \mathbb{E} [Z(\mathbf{s} + \mathbf{h})Z(\mathbf{s})] \qquad (B.3) $$

Since $Z$ the second moment $\mathbb{E}[Z(\mathbf{s})^2] = \operatorname{var}(Z(\mathbf{s}))$ is constant and therefore independent of $\mathbf{s}$ for all $\mathbf{s}$ of $\mathcal{S}$,

$$ \gamma_Z(\mathbf{h}) = \mathbb{E}[Z(\mathbf{s})^2] - \mathbb{E}[Z(\mathbf{s} + \mathbf{h})Z(\mathbf{s})] \qquad (B.4) $$

$$ = \operatorname{var}(Z(\mathbf{s})) - \operatorname{cov}(Z(\mathbf{s}), Z(\mathbf{s} + \mathbf{h})) \qquad (B.5) $$

With $\mathbf{h} = \mathbf{s}_i - \mathbf{s}_j$. From this equation, we can deduce the covariance matrix of $\mathbf{x}_s$:

$$ \Sigma_s = (\operatorname{var}(Z(\mathbf{s})) - \gamma_Z(\mathbf{s}_i - \mathbf{s}_j))_{1\le i,j \le n} \qquad (B.6) $$
# Appendix C. Kalman Derivations

## Appendix C.1. Kalman Prediction and Correction Steps

Based on (14), let note $\boldsymbol{x}_{d,k|k}^p$ the estimate of the dynamic component at instant $k$ and $\boldsymbol{x}_{d,k|k-1}^p$ the predicted estimate at instant $k$ before the correction for the considered particle. The estimate $\boldsymbol{x}_{d,k|k-1}^p$ and covariance matrix $\boldsymbol{\Sigma}_{d,k|k-1}^p$ can be written

$$ 
\left\{ 
\begin{aligned}
\boldsymbol{x}_{d,k|k-1}^p &= \mathbf{A}_k(\boldsymbol{\theta}_k^p)\boldsymbol{x}_{d,k-1|k-1}^p \\
\boldsymbol{\Sigma}_{d,k|k-1}^p &= \mathbf{A}_k(\boldsymbol{\theta}_k^p)\boldsymbol{\Sigma}_{d,k-1|k-1}^p\mathbf{A}_k(\boldsymbol{\theta}_k^p)^T + \boldsymbol{\Sigma}_\eta
\end{aligned}
\right. 
\quad (C.1) 
$$

By using the Kalman equations, the correction step can therefore be derived, using the optimal Kalman gain (see Appendix C.2)

$$ 
\left\{ 
\begin{aligned}
\dot{\boldsymbol{i}}_k^p &= \boldsymbol{y}_k - \mathbf{C}_k(\boldsymbol{\theta}_k^p)(\boldsymbol{\mu}_k^p + \boldsymbol{x}_{d,k|k-1}^p) \\
\boldsymbol{S}_k^p &= \mathbf{C}_k(\boldsymbol{\theta}_k^p)(\boldsymbol{\Sigma}_{d,k|k-1}^p + \boldsymbol{\Sigma}_s)\mathbf{C}_k(\boldsymbol{\theta}_k^p)^T + \boldsymbol{\Sigma}_{\nu,k} \\
\boldsymbol{K}_k^p &= \boldsymbol{\Sigma}_{d,k|k-1}^p\mathbf{C}_k(\boldsymbol{\theta}_k^p)^T(\boldsymbol{S}_k^p)^{-1} \\
\boldsymbol{x}_{d,k|k}^p &= \boldsymbol{x}_{d,k|k-1}^p + \boldsymbol{K}_k^p\dot{\boldsymbol{i}}_k^p \\
\boldsymbol{\Sigma}_{d,k|k}^p &= (\boldsymbol{I} - \boldsymbol{K}_k^p\mathbf{C}_k(\boldsymbol{\theta}_k^p))\boldsymbol{\Sigma}_{d,k|k-1}^p
\end{aligned}
\right. 
\quad (C.2) 
$$

## Appendix C.2. Kalman Gain Derivation

For a given particle, let compute first the *posterior* estimate covariance matrix $\boldsymbol{\Sigma}_{d,k|k}^p$.

$$ \boldsymbol{\Sigma}_{d,k|k}^p = \operatorname{cov}(\boldsymbol{x}_{d,k} - \boldsymbol{x}_{d,k|k}^p) \quad (C.3) $$

$$ = \operatorname{cov}(\boldsymbol{x}_{d,k} - (\boldsymbol{x}_{d,k|k-1}^p + \boldsymbol{K}_k^p \dot{\boldsymbol{i}}_k^p)) \quad (C.4) $$

$$ = \operatorname{cov}(\boldsymbol{x}_{d,k} - (\boldsymbol{x}_{d,k|k-1}^p + \boldsymbol{K}_k^p(\boldsymbol{y}_k - \mathbf{C}_k(\boldsymbol{\theta}_k^p)\boldsymbol{x}_{d,k|k-1}^p))) \quad (C.5) $$

Since $\mathbf{C}_k(\boldsymbol{\theta}_k^p)\boldsymbol{\mu}_k^p$ is a real-valued vector. Then,
$$ \Sigma_{d,k|k}^{p} = \operatorname{cov}(\mathbf{x}_{d,k} - \mathbf{x}_{d,k|k-1}^{p}) \qquad (C.6) $$

$$ - \mathbf{K}_{k}^{p}[\mathbf{C}_{k}(\boldsymbol{\theta}_{k}^{p})(\mathbf{x}_{d,k} + \mathbf{x}_{s,k}) + \boldsymbol{\nu}_{k}] \qquad (C.7) $$

$$ - \mathbf{C}_{k}(\boldsymbol{\theta}_{k}^{p})\mathbf{x}_{d,k|k-1}^{p}] \qquad (C.8) $$

$$ = \operatorname{cov}((\mathbf{I} - \mathbf{K}_{k}^{p}\mathbf{C}_{k}(\boldsymbol{\theta}_{k}^{p}))(\mathbf{x}_{d,k} - \mathbf{x}_{d,k|k-1}^{p}) - \mathbf{K}_{k}^{p}\boldsymbol{\nu}_{k} \qquad (C.9) $$

$$ - \mathbf{K}_{k}^{p}\mathbf{C}_{k}(\boldsymbol{\theta}_{k}^{p})\mathbf{x}_{s,k} \qquad (C.10) $$

$$ = (\mathbf{I} - \mathbf{K}_{k}^{p}\mathbf{C}_{k}(\boldsymbol{\theta}_{k}^{p}))\Sigma_{d,k|k-1}^{p}(\mathbf{I} - \mathbf{K}_{k}^{p}\mathbf{C}_{k}(\boldsymbol{\theta}_{k}^{p}))^{T} \qquad (C.11) $$

$$ + \mathbf{K}_{k}^{p}\boldsymbol{\Sigma}_{\nu}\mathbf{K}_{k}^{pT} + \mathbf{K}_{k}^{p}\mathbf{C}_{k}(\boldsymbol{\theta}_{k}^{p})\boldsymbol{\Sigma}_{s}(\mathbf{K}_{k}^{p}\mathbf{C}_{k}(\boldsymbol{\theta}_{k}^{p}))^{T} \qquad (C.12) $$

* by assuming that $\mathbf{x}_{s,k}$, $\mathbf{x}_{d,k}$ and $\boldsymbol{\nu}_{k}$ are independent and using the definition of $\Sigma_{s} = \operatorname{cov}(\mathbf{x}_{s,k})$.

The Kalman filter aims at minimizing the quantity $\mathbb{E}[\|\mathbf{x}_{d,k} - \mathbf{x}_{d,k|k}^{p}\|^{2}]$ which is equivalent to the minimization of the trace of the a posteriori estimate covariance matrix $\Sigma_{d,k|k}^{p}$.

$$ \frac{\partial \operatorname{tr}(\Sigma_{d,k|k}^{p})}{\partial \mathbf{K}_{k}^{p}} = -2(\mathbf{C}_{k}(\boldsymbol{\theta}_{k}^{p})\Sigma_{d,k|k-1}^{p})^{T} + 2\mathbf{K}_{k}^{p}\mathbf{S}_{k}^{p} = 0 \qquad (C.13) $$

The Kalman gain is therefore given by:

$$ \mathbf{K}_{k}^{p} = \Sigma_{d,k|k-1}^{p} \mathbf{C}_{k}(\boldsymbol{\theta}_{k}^{p}) \mathbf{S}_{k}^{p-1} \qquad (C.14) $$

## Appendix D. CMA-ES optimization problem

Let $\mathbf{w}_s$ be a window of $N_w$ pixels centered over a pixel $\mathbf{s} \in \mathcal{S} \setminus \partial\mathcal{S}$. The local homogeneity of temperature and emissivity can be stated as:

$$ \forall \tilde{\mathbf{s}} \in \mathbf{w}_s, \quad \begin{cases} \epsilon(\tilde{\mathbf{s}}) = \epsilon(\mathbf{s}) \\ T(\tilde{\mathbf{s}}) = T(\mathbf{s}) \end{cases} \qquad (D.1) $$

The optimization problem is therefore derived, $\forall \mathbf{s} \in \mathcal{S} \setminus \partial\mathcal{S}$:

$$ (\epsilon(\mathbf{s}), T(\mathbf{s})) = \arg\min_{(\epsilon, T)} \left( \sum_{\tilde{\mathbf{s}} \in \mathbf{w}_s} \|\mathbf{y}_k(\tilde{\mathbf{s}}) - \epsilon L^\circ(T, \lambda_i)\| \right) \qquad (D.2) $$
## References

[1] N. Le Touz, T. Toullier, J. Dumoulin, Study of an optimal heating command law for structures with non-negligible thermal inertia in varying outdoor conditions, Smart Structures and Systems 27 (2) (2021) 379–386. doi:10.12989/sss.2021.27.2.379.
URL https://inria.hal.science/hal-03145348

[2] B. Merainani, T. Toullier, S. Sriranjan, J. Dumoulin, Monitoring of moving rail-road cars through infrared thermal vision: Study of hot box detection on reconstructed scenes with and without Deep Learning approaches, in: QIRT 2024 - 17th Quantitative Infrared Thermography Conference, Zagreb, Croatia, 2024, pp. 1–2.
URL https://inria.hal.science/hal-04798560

[3] N. Lalanne, Modélisation des échanges thermiques et radiatifs en environnement urbain à très haute résolution spatiale : aide à l'interprétation des mesures par télédétection infrarouge, Ph.D. thesis, thèse de doctorat dirigée par Le Niliot, Christophe et Krapez, Jean-Claude Energétique Aix-Marseille 2015 (2015).
URL http://www.theses.fr/2015AIXM4731

[4] N. Lalanne, J. Krapez, C. Le Niliot, X. Briottet, J. Pierro, L. Labarre, Development and validation of a numerical tool for the simulation of the temperature field and infrared radiance rendering in an urban scene, in: 12th International Conference on Quantitative Infrared Thermography, BORDEAUX, France, 2014.
URL https://hal.science/hal-01087770

[5] J.-C. Krapez, Radiative Measurements of Temperature, Heat Transfer, CRC Press, Boca Raton, FL, 2011.

[6] T. Toullier, J. Dumoulin, Bias and bottlenecks study in outdoor long term thermal monitoring by infrared thermography: Leveraging opportunistic data for temperature estimation, Infrared Physics and Technology 141 (2024) 105471. doi:10.1016/j.infrared.2024.105471.
URL https://inria.hal.science/hal-04677950

[7] J. Howell, R. Siegel, M. Pinar, Thermal Radiation Heat Transfer., 5th Edition, CRC Press, 2010.
[8] S. B. A. Rolim, A. Grondona, C. L. Hackmann, C. Rocha, A Review of Temperature and Emissivity Retrieval Methods: Applications and Restrictions, American Journal of Environmental Engineering 6 (2016) 119–128.

[9] Z.-L. Li, B.-H. Tang, H. Wu, H. Ren, G. Yan, Z. Wan, I. F. Trigo, J. A. Sobrino, Satellite-derived land surface temperature: Current status and perspectives, Remote Sensing of Environment 131 (2013) 14–37. doi:10.1016/j.rse.2012.12.008.

[10] J. N. Ash, J. Meola, Temperature-emissivity separation for LWIR sensing using MCMC, in: Algorithms and Technologies for Multispectral, Hyperspectral, and Ultraspectral Imagery XXII, Vol. 9840, International Society for Optics and Photonics, 2016, p. 98401O. doi: 10.1117/12.2223263.

[11] G. Masiello, C. Serio, I. De Feis, M. Amoroso, S. Venafra, I. F. Trigo, P. Watts, Kalman filter physical retrieval of surface emissivity and temperature from geostationary infrared radiances, Atmospheric Measurement Techniques 6 (12) (2013) 3613–3634. doi:10.5194/amt-6-3613-2013.

[12] M. Pivovarník, S. J. S. Khalsa, J. C. Jiménez-Muñoz, F. Zemek, Improved Temperature and Emissivity Separation Algorithm for Multispectral and Hyperspectral Sensors, IEEE Transactions on Geoscience and Remote Sensing 55 (4) (2017) 1944–1953. doi:10.1109/TGRS.2016.2631508.

[13] A. Barducci, D. Guzzi, C. Lastri, P. Marcoionni, V. Nardino, I. Pippi, Emissivity and Temperature Assessment Using a Maximum Entropy Estimator: Structure and Performance of the MaxEnTES Algorithm, IEEE Transactions on Geoscience and Remote Sensing 53 (2) (2015) 738–751. doi:10.1109/TGRS.2014.2327218.

[14] N. Wang, Z.-L. Li, B.-H. Tang, F. Zeng, C. Li, Retrieval of atmospheric and land surface parameters from satellite-based thermal infrared hyperspectral data using a neural network technique, International Journal of Remote Sensing 34 (9-10) (2013) 3485–3502. doi: 10.1080/01431161.2012.716536.
[15] J. F. Mas, J. J. Flores, The application of artificial neural networks to the analysis of remotely sensed data, International Journal of Remote Sensing 29 (3) (2008) 617-663. doi:10.1080/01431160701352154.

[16] J. Sobrino, Zhao-Liang Li, M. Stoll, F. Becker, Improvements in the split-window technique for land surface temperature determination, IEEE Transactions on Geoscience and Remote Sensing 32 (2) (1994) 243-253. doi:10.1109/36.295038.

[17] Z. Wan, J. Dozier, A generalized split-window algorithm for retrieving land-surface temperature from space, IEEE Transactions on Geoscience and Remote Sensing 34 (4) (1996) 892-905. doi:10.1109/36.508406.

[18] F. Becker, Z.-L. Li, Towards a local split window method over land surfaces, International Journal of Remote Sensing 11 (3) (1990) 369-393. doi:10.1080/01431169008955028.

[19] L. M. McMillin, Estimation of sea surface temperatures from two infrared window measurements with different absorption, Journal of Geophysical Research (1896-1977) 80 (36) (1975) 5113-5117. doi:10.1029/JC080i036p05113.

[20] F. Becker, Z.-L. Li, Temperature-independent spectral indices in thermal infrared bands, Remote Sensing of Environment 32 (1) (1990) 17-33. doi:10.1016/0034-4257(90)90095-4.

[21] Z.-L. Li, F. Becker, Feasibility of land surface temperature and emissivity determination from AVHRR data, Remote Sensing of Environment 43 (1) (1993) 67-85. doi:10.1016/0034-4257(93)90065-6.

[22] Z. Li, F. Petitcolin, R. Zhang, A physically based algorithm for land surface emissivity retrieval from combined mid-infrared and thermal infrared data, Science in China Series E: Technological Sciences 43 (1) (2000) 23-33. doi:10.1007/BF02916575.

[23] A. Gillespie, S. Rokugawa, T. Matsunaga, J. Cothern, S. Hook, A. Kahle, A temperature and emissivity separation algorithm for advanced spaceborne thermal emission and reflection radiometer (aster) images, IEEE Transactions on Geoscience and Remote Sensing 36 (4) (1998) 1113-1126. doi:10.1109/36.700995.
[24] A. R. Gillespie, E. A. Abbott, L. Gilson, G. Hulley, J.-C. Jiménez-Muñoz, J. A. Sobrino, Residual errors in ASTER temperature and emissivity standard products AST08 and AST05, Remote Sensing of Environment 115 (12) (2011) 3681–3694. doi:10.1016/j.rse.2011.09.007.

[25] U. Avdan, G. Jovanovska, Algorithm for automated mapping of land surface temperature using landsat 8 satellite data, Journal of Sensors 2016 (1) (2016) 1480307. arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1155/2016/1480307, doi:https://doi.org/10.1155/2016/1480307. URL https://onlinelibrary.wiley.com/doi/abs/10.1155/2016/1480307

[26] J. C. Jiménez-Muñoz, J. A. Sobrino, D. Skoković, C. Mattar, J. Cristóbal, Land surface temperature retrieval methods from landsat-8 thermal infrared sensor data, IEEE Geoscience and Remote Sensing Letters 11 (10) (2014) 1840–1843. doi:10.1109/LGRS.2014.2312032.

[27] R. E. Kalman, A New Approach to Linear Filtering and Prediction Problems, Journal of Basic Engineering 82 (1) (1960) 35–45. doi:10.1115/1.3662552.

[28] L. Ljung, Asymptotic behavior of the extended Kalman filter as a parameter estimator for linear systems, IEEE Transactions on Automatic Control 24 (1) (1979) 36–50. doi:10.1109/TAC.1979.1101943.

[29] J. Ching, J. L. Beck, K. A. Porter, Bayesian state and parameter estimation of uncertain dynamical systems, Probabilistic Engineering Mechanics 21 (1) (2006) 81–96. doi:10.1016/j.probengmech.2005.08.003.

[30] C. Snyder, T. Bengtsson, P. Bickel, J. Anderson, Obstacles to High-Dimensional Particle Filtering, Monthly Weather Review 136 (12) (2008) 4629–4640. doi:10.1175/2008MWR2529.1.

[31] P. Del Moral, Nonlinear filtering: Interacting particle resolution, Comptes Rendus de l'Académie des Sciences - Series I - Mathematics 325 (6) (1997) 653–658. doi:10.1016/S0764-4442(97)84778-7.
[32] T. Schon, F. Gustafsson, P. Nordlund, Marginalized particle filters for mixed linear/nonlinear state-space models, IEEE Transactions on Signal Processing 53 (7) (2005) 2279-2289. doi:10.1109/TSP.2005.849151.

[33] P. D. Moral, Feynman-Kac Formulae: Genealogical and Interacting Particle Systems with Applications, Probability and Its Applications, Springer-Verlag, New York, 2004. doi:10.1007/978-1-4684-9393-1.

[34] Z. Chen, Bayesian Filtering: From Kalman Filters to Particle Filters, and Beyond, Statistics 182 (Jan. 2003). doi:10.1080/02331880309257.

[35] N. A. C. Cressie, C. K. Wikle, Statistics for Spatio-Temporal Data, Wiley Series in Probability and Statistics, Wiley, Hoboken, N.J, 2011.

[36] B. Baingana, E. Dall'Anese, G. Mateos, G. B. Giannakis, Robust kriged Kalman filtering, in: 2015 49th Asilomar Conference on Signals, Systems and Computers, 2015, pp. 1525-1529. doi:10.1109/ACSSC.2015.7421400.

[37] K. V. Mardia, C. Goodall, E. J. Redfern, F. Alonso, The Kriged Kalman Filter, TEST: An Official Journal of the Spanish Society of Statistics and Operations Research 7 (1998) 217-282. doi:10.1007/BF02565111.

[38] W. S. Kerwin, The Kriging Update Model and Recursive Space-Time Function Estimation, IEEE Transactions on Signal Processing 47 (11) (1999) 11.

[39] V. Roy, A. Simonetto, G. Leus, Spatio-Temporal Field Estimation Using Kriged Kalman Filter (KKF) with Sparsity-Enforcing Sensor Placement, Sensors 18 (6) (2018) 1778. doi:10.3390/s18061778.

[40] M. Zghal, L. Mevel, P. Del Moral, Modal parameter estimation using interacting Kalman filter, Mechanical Systems and Signal Processing 47 (1-2) (2014) 139-150. doi:10.1016/j.ymssp.2012.11.005.

[41] S. Sen, A. Crinière, L. Mevel, F. Cérou, J. Dumoulin, Correntropy based IPKF filter for parameter estimation in presence of non-stationary noise process, IFAC-PapersOnLine 51 (24) (2018) 420-427. doi:10.1016/j.ifacol.2018.09.611.
[42] G. Matheron, F. Blondel, Bureau de recherches géologiques et minières (France), Traité de géostatistique appliquée. Tome I Tome I, Technip, Paris, 1962.

[43] T. Toullier, J. Dumoulin, L. Mevel, Study and development of an image synthesis software for the simulation of in-situ 3D scenes monitored by infrared cameras, in: QIRT 2020 - 15th Quantitative InfraRed Thermography Conference, 2020.

[44] O. Roustant, Y. Deville, (Re)Dice Packages for Computer Experiments, in: Troisièmes rencontres R, Montpellier, France, 2014.

[45] T. Toullier, J. Dumoulin, L. Mevel, Étude comparative de quatre méthodes statistiques pour l'estimation conjointe de l'émissivité et de la température par thermographie infrarouge multispectrale, in: SFT 2020 - Congrès annuel de la Société Française de Thermique, Belfort, France, 2020, pp. 1-9. doi:10.25855/SFT2020-087. URL https://hal.inria.fr/hal-02986521

[46] T. Toullier, J. Dumoulin, L. Mevel, A kriging-based interacting particle kalman filter for the simultaneous estimation of temperature and emissivity in infra-red imaging, IFAC-PapersOnLine 53 (2) (2020) 456-461, 21th IFAC World Congress. doi:https://doi.org/10.1016/j.ifacol.2020.12.219.

[47] C. Berrett, G. P. Williams, T. Moon, J. Gunther, A Bayesian Non-parametric Model for Temperature-Emissivity Separation of Long-Wave Hyperspectral Images, Technometrics 56 (2) (2014) 200-211. doi: 10.1080/00401706.2013.869262.

[48] T. Toullier, Simultaneous characterization of objects temperature and radiative properties through multispectral infrared thermography, Ph.D. thesis, Université de Rennes 1, IFSTTAR - Nantes (2019).

[49] N. Hansen, The CMA Evolution Strategy: A Comparing Review, Studies in Fuzziness and Soft Computing, Springer, 2006. doi:10.1007/3-540-32494-1_4.

[50] S. García, D. Molina, M. Lozano, F. Herrera, A study on the use of non-parametric tests for analyzing the evolutionary algorithms' be-
haviour: A case study on the CEC'2005 Special Session on Real Pa-
rameter Optimization, Journal of Heuristics 15 (6) (2009) 617-644.
doi:10.1007/s10732-008-9080-4.

[51] N. Hansen, A. Auger, R. Ros, S. Finck, P. Pošík, Comparing results of 31 algorithms from the black-box optimization benchmarking BBOB-2009, in: Proceedings of the 12th Annual Conference Comp on Genetic and Evolutionary Computation - GECCO 10, ACM Press, Portland, Oregon, USA, 2010, p. 1689. doi:10.1145/1830761.1830790.

[52] R. J. Howarth, (A. G.) Journel and (Ch. J.) Huijbregts. Mining Geo-statistics., Mineralogical Magazine 43 (328) (1979) 563-564. doi: 10.1180/minmag.1979.043.328.34.

[53] G. Matheron, Le Krigeage Universel (Universal Kriging), Vol. 1 of Les Cahiers Du Centre de Morphologie Mathématique, École Nationale Supérieure des Mines de Paris, Fontainebleau, 1969. doi:10.1016/j.geoderma.2003.08.018.