domain: math sde

text:
Igor Yu Gejadze

Victor Shutyaev
email: shutyaev@inm.ras.ru

Francois-Xavier Le Dimet

Hessian-based covariance approximations in variational data assimilation

Keywords: variational data assimilation, optimality system, Hessian, estimation error covariance MSC: 65K10, 86A22, 93B40

The problem of variational data assimilation (estimation) for a nonlinear model is considered in general operator formulation. Hessian-based methods are presented to compute the estimation error covariances. The importance of dynamic formulation and the role of the Hessian and its inverse are discussed.

Introduction

Data assimilation (state or parameter estimation) for high-dimensional distributed parameter dynamical systems has become a powerful analysis tool in recent decades. It is a fusion process of incomplete and, possibly, indirect observations of state variables with a mathematical model governed by partial differential equations, complemented by a priori information. The applications include model initialization in meteorology and oceanography, air and water quality monitoring, 'calibration' of groundwater and reservoir models, discharge estimation and forecasting in river hydraulics and hydrology, flow estimation and control in aerospace engineering, process control in chemical and nuclear engineering, etc. In different applications these estimation problems are also referred as 'inverse problems', 'data assimilation'(DA), and 'calibration'. Variational methods were introduced in meteorology by Sasaki [START_REF] Sasaki | An objective analysis based on the variational method[END_REF]. These methods consider the equations governing the flow as constraints and the problem is closed by using a variational principle, e.g. the minimization of the discrepancy between the model prediction and the observations. Optimal Control Approach (Lions [19], Marchuk [START_REF] Marchuk | Perturbation theory and the statement of inverse problems[END_REF]) for data assimilation was proposed by Marchuk and Penenko [START_REF] Marchuk | Application of optimization methods to the problem of mathematical simulation of atmospheric processes and environment // Modelling and Optimization of Complex Systems[END_REF], Penenko and Obraztsov [START_REF] Penenko | A variational initialization method for the fields of the meteorological elements[END_REF], Le Dimet [START_REF] Dimet | A general formalism of variational analysis[END_REF], Le Dimet and Talagrand [START_REF] Dimet | Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects[END_REF], Talagrand and Courtier [START_REF] Talagrand | Variational assimilation of meteorological observations with the adjoint vorticity equation[END_REF]. Therefore, variational estimation/DA is a method based on the optimal control theory, which can also be understood as a special case of the maximum a-posteriory probability (MAP) estimator [START_REF] Dashti | Map estimators and their consistency in bayesian nonparametric inverse problems[END_REF]. This method is preferred for weather and ocean forecasting in major operational centers around the globe, particularly in the form of the incremental 4D-Var [START_REF] Courtier | A strategy for operational implementation of 4d-var, using an incremental approach[END_REF], and in the form of the ensemble 4D-Var [START_REF] Clayton | Operational implementation of a hybrid ensemble/4d-var global data assimilation system at the met office[END_REF]. Variational estimation is widely used in other scientific and engineering applications, such as aerospace engineering [START_REF] Alifanov | Extreme Methods for Solving Ill-Posed Problems with Applications to Inverse Heat Transfer Problems[END_REF] and astrophysics [START_REF] Bardsley | Total variation-penalized poisson likelihood estimation for ill-posed problems[END_REF], to mention a few.

Uncertainty quantification is an important topic closely associated with estimation and data assimilation. An overview of the original work by the authors on the Hessian-based advanced methodology for computing the estimation error covariances in the framework of variational estimation is presented in the coming sections. The main results are given in a general operator formulation. The accent is made on feasibility of the suggested methods for the uncertainty quantification in high-dimensions, where the statistical methods (Monte Carlo involving associated tricks, e.g. localization, importance sampling, etc.) may not produce a sensible outcome due to a very small sample being available. Some of the approaches could equally be useful for improving (in terms of acceleration, memory savings and robustness) the estimation/DA technique itself.

Variational data assimilation: basic concept

Let us consider a mathematical (or numerical) model which describes behavior of a natural system in terms of its state variables X ∈ X . Let U ∈ U be the set of the model inputs (controls), then the model can be considered as a 'control-to-state' mapping M : U → X , such that

X = M(U), (2.1)
where U and X are the input and state spaces, correspondingly. For modeling the system behavior the true input vector Ū must be specified. Under the 'perfect model' assumption the following can be postulated: X = M( Ū). In reality, some components of Ū contain uncertainties ε ∈ U. Thus, instead of Ū we use its best available approximation (background/prior)

U * = Ū + ε, (2.2)
where ε is also called the background error. Because of the presence of ε, the predicted state X|U * = M(U * ), that is, X evaluated (or conditioned) on U * , also contains an error δX = M(U * ) -M( Ū).

The state observing tools are represented by an observation operator

C : X → Y in the form Y = C(X) = C(M(U)) := G(U), (2.3)
where G : U → Y is a generalized input-to-observations mapping and Y is the 'observation' space. The true observations would be Ȳ = G( Ū), however the actual observations usually contain noise ξ (observation uncertainty), i.e.

Y * = Ȳ + ξ. (2.4)
The aim of data assimilation is to obtain Û = U|Y * , i.e. an estimate of U conditioned on observations Y * , which should be better than the prior U * in the sense Û -Ū < U * -Ū .

In the Bayesian framework the posterior probability density of U conditioned on observations Y * is given by the Bayes formula [START_REF] Stuart | Inverse problems: a Bayesian perspective[END_REF] p

(U|Y * ) = p(Y * |U)p(U) p(Y * ) . (2.5)
Looking for the mode of the posterior density p(U|Y * ), i.e. maximizing p(U|Y * ), is the essence of variational data assimilation (estimation). Under the Gaussian assumption on the prior and observation uncertainties, i.e. ε ∼ N(0, B), ξ ∼ N(0, R), where B is the background error covariance and R -the observation error covariance, maximizing p(U|Y * ) is equivalent to minimizing the cost-function [START_REF] Tarantola | Inverse Problem Theory and Methods for Model Parameter Estimation[END_REF] J

(U) = 1 2 R -1/2 (G(U) -Y * ) 2 Y + 1 2 B -1/2 (U -U * ) 2 U . (2.6)
Thus, the estimate Û is obtained from the necessary optimality condition

J ′ U ( Û) = 0. (2.7)
For the operator G(U) we define the tangent linear operator G ′ (U) (Gateaux derivative) and its adjoint (G ′ (U)) * [START_REF] Marchuk | Adjoint Equations and Perturbation Algorithms in Nonlinear Problems[END_REF] as follows:

G ′ U (U)w = lim t→0 G(U + tw) -G(U) t , w ∈ U, (2.8)
(w, (G ′ U (U)) * w * ) U = (G ′ U (U)w, w * ) Y , w * ∈ Y.
(2.9)

Given the above operator definitions, the full gradient of J(u) in (2.7) can be expressed in the form:

J ′ U (U) = (G ′ U (U)) * R -1 (G(U) -Y * ) + B -1 (U -U * ). (2.10)
Thus, the estimate Û is the solution to the operator equation

(G ′ U ( Û)) * R -1 (G( Û) -Y * ) + B -1 ( Û -U * ) = 0. (2.11)
The above presented general operator formulation of DA/estimation problems is common in nonlinear regression and in the inverse problems theory.

In the theory of dynamical systems it is usual to formulate DA/estimation problems in terms of non-stationary partial differential equations (PDE) [START_REF] Lions | Contrôle optimal des systèmes gouvernés par des équations aux dérivées partielles[END_REF]. Let X be the state space such that X = L 2 (0, T ; Ω), with a norm

• X = (•, •) 1/2
X , where Ω is a bounded domain of the natural space R d (d = 1, 2 or 3). For the initial state control problem considered in [START_REF] Gejadze | On analysis error covariances in variational data assimilation[END_REF] one can write

   ∂ϕ ∂t = F (ϕ) + f, t ∈ (0, T ) ϕ| t=0 = u, (2.12)
where ϕ = ϕ(x, t) (analog of X) is the unknown function belonging for any time instant t to L 2 (Ω), x ∈ R d is the spatial variable, u(x) ∈ L 2 (Ω), f (x) ∈ X are the initial state and the source term, respectively, and F (x) is a nonlinear operator mapping L 2 (Ω) into L 2 (Ω). Suppose that for a given pair (u, f ) there exists a unique solution ϕ ∈ X to (2.12). Next, we introduce the functional

J(u) = 1 2 R -1/2 (C(ϕ) -Y * ) 2 Y + 1 2 B -1/2 (u -u * ) 2 L 2 (Ω) , (2.13)
where u * is a prior (background) of u. Thus, operator G is defined by the formula G(u) = C(ϕ), where ϕ is the solution to equation (2.12). Now, consider the following data assimilation problem with the aim to identify the initial condition: find u ∈ L 2 (Ω) and ϕ ∈ X such that they satisfy (2.12), and on the set of solutions to (2.12), the functional J(u) takes the minimum value, i.e.

         ∂ϕ ∂t = F (ϕ) + f, t ∈ (0, T ) ϕ| t=0 = u, J(u) = inf v∈L 2 (Ω) J(v).
(2.14)

The necessary optimality condition reduces the problem (2.14) to the following optimality system [START_REF] Lions | Contrôle optimal des systèmes gouvernés par des équations aux dérivées partielles[END_REF], [START_REF] Agoshkov | On solvability and numerical solution of data assimilation problems[END_REF], [START_REF] Marchuk | Approaches to the solution of data assimilation problems[END_REF]:

   ∂ϕ ∂t = F (ϕ) + f, t ∈ (0, T ) ϕ| t=0 = û, (2.15)    -∂ϕ * ∂t -(F ′ (ϕ)) * ϕ * = -(C ′ (ϕ)) * R -1 (C(ϕ) -Y * ), t ∈ (0, T ) ϕ * | t=T = 0, (2.16
)

B -1 (û -u * ) -ϕ * | t=0 = 0 (2.17)
with the unknowns ϕ, ϕ * , û, where (F ′ (ϕ)) * is the adjoint to the Frechet derivative of F , and (C ′ (ϕ)) * is the adjoint to the Frechet derivative of C defined by

(C ′ (ϕ)η, ψ) Y = (η, (C ′ (ϕ)) * ψ) X , η ∈ X , ψ ∈ Y.
We assume that the system (2.15)-(2.17) has a unique solution. It is easy to see that this system is a detalized representation of the estimator equation (2.11) in the case the model M in (2.1) is a dynamical model and the input U = {u} is the initial state ϕ(x, 0). Obviously, the general results obtained for (2.11) are valid for the system (2.15)-(2.17). However, the latter is used for better understanding of the structure of adjoint operators, which is sometimes vital for practical implementation of the variational DA method.

3 Estimation error covariance and the inverse Hessian

Clarification of the existing theory

Here, we use the main results from [START_REF] Gejadze | On analysis error covariances in variational data assimilation[END_REF] and [START_REF] Gejadze | On optimal solution error covariances in variational data assimilation problems[END_REF] using the general operator formulation, whereas in papers the dynamic formulation is used. The latter has its own value (see Section 3.3).

Let us consider an estimation error δU = Û -Ū and equation (2.11). We notice that

G( Û ) -Y * = G( Û ) -(G( Ū) + ξ) = G ′ U ( Ũ)δU -ξ,
where Ũ = Ū + τ δU, τ ∈ [0, 1], and

Û -U * = ( Û -Ū ) -(U * -Ū ) = δU -ε.
Then, equation (2.11) yields the error equation

(G ′ U ( Û )) * R -1 (G ′ U ( Ũ )δU -ξ) + B -1 (δU -ε) = 0, (3.18)
from where the estimation error can be expressed:

δU = (G ′ U ( Û )) * R -1 G ′ U ( Ũ ) + B -1 -1 (G ′ U ( Û )) * R -1 ξ + B -1 ε . (3.19)
Using the first order approximations for operators

G ′ U ( Û ) ≈ G ′ U ( Ũ ) ≈ G ′ U ( Ū ) we express δU as follows: δU ≃ H -1 ( Ū ) (G ′ U ( Ū)) * R -1 ξ + B -1 ε , (3.20)
where

H(•) = (G ′ U (•)) * R -1 G ′ U (•) + B -1 . (3.21)
Here, (•) denotes a placeholder for the argument of operators G and H, which shall be called the 'origin'.

It has been shown that H(•) represents the first-order term in the Hessian of the cost-function J(U) in (2.6) or, otherwise, is the Hessian of the auxiliary cost-function as follows:

J (•, δU) = 1 2 R -1/2 (G ′ U (•)δU -f 1 ) 2 Y + 1 2 B -1/2 (δU -f 2 ) 2 U , (3.22)
where f 1 ∈ Y and f 2 ∈ U are trial functions, including f 1 = 0 and f 2 = 0. We assume that H is positive definite and, hence, invertible. If the errors ξ and ε indeed satisfy the conditions ε ∼ N(0, B), ξ ∼ N(0, R), then the estimation error covariance is

P = E[δUδU T ] ≃ H -1 ( Ū ). (3.23)
The above relationship is exact for linear G; for nonlinear G it is valid if the tangent linear hypothesis is valid. Since in reality the true value Ū may not be known, one uses the particular estimate (event) Û instead, i.e. P ≃ H -1 ( Û ).

(

3.24)

There are several papers presenting essentially the same result for dynamical models [START_REF] Rabier | Four-dimensional assimilation in the presence of baroclinic instability[END_REF][START_REF] Thacker | The role of the hessian matrix in fitting models to measurements[END_REF]. However, there has also been some confusion. For example, H(•) is often perceived as a Hessian of the cost-function J(U) in (2.6), whereas it is the Hessian of the auxiliary cost-function J (U) in (3.22). Another important point is that H -1 ( Ū) is the best possible approximation to P using H, whereas H( Û) is only an approximation to H( Ū ). The latter fact has never been clearly underlined before.

Computing the inverse Hessian using the LBFGS

First, we define the preconditioned (projected) Hessian in the form

H = (B 1/2 ) * HB 1/2 = (B 1/2 ) * (G ′ U (•)) * R -1 G ′ U (•)B 1/2 + I,
which is usually much better conditioned that H. Having computed H-1 in the limited-memory form H -1 is recovered as follows:

H -1 = B 1/2 H-1 (B 1/2 ) * .
It is easy to see that H is the Hessian of the following auxiliary cost-function

J (•, δU) = 1 2 R -1/2 (G ′ U (•)B 1/2 δU -f 1 ) 2 Y + 1 2 δU -f 2 2 U . (3.25)
We use the LBFGS method [START_REF] Liu | On the limited memory BFGS method for large scale minimization[END_REF] for generating the limited-memory representation of H-1 as the by-product of minimization of J (•, δU). A set of secant pairs of a given (prescribed) size is accumulated during minimization, then the product H-1 v is recovered by the BFGS recursion.

The key improvement of the quality of approximation (for a given number of secant pairs) has been achieved by using the exact step search, as described in Section 5.1 in [START_REF] Gejadze | On analysis error covariances in variational data assimilation[END_REF]. The optimized (in terms of computational expenses) version of this algorithm is presented in our later paper [START_REF] Gejadze | Computation of the analysis error covariance in variational data assimilation problems with nonlinear dynamics[END_REF], §6.2.

The Lanczos method is routinely used for computing H-1 , whereas the BFGS or LBFGS are usually seen as inferior to the Lanczos. Indeed, for the same accuracy of approximation achieved, representation using the BFGS update formula is less compact than the eigenvalue representation. The eigenvalue analysis, however, could be noticeably more expensive to compute (in terms of function calls). Besides, control of computational costs in the LBFGS algorithm is easier than in the Lanczos method.

We can cite only [START_REF] Veerse | Variable-storage quasi-newton operators as inverse forecast/analysis error covariance matrices in variational data assimilation[END_REF] (preceding to [START_REF] Gejadze | On analysis error covariances in variational data assimilation[END_REF][START_REF] Gejadze | On optimal solution error covariances in variational data assimilation problems[END_REF] in time) where the same technique is considered. The main focus in this paper is, however, on methods for computing the inverse of the analysis error covariance P , given the BFGS representation of H -1 obtained in the inner loop of the incremental approach. At the same time, matrix H -1 itself is not of interest and its accuracy is not assessed.

On the importance of dynamic formulation

Let us note that in papers [START_REF] Gejadze | On analysis error covariances in variational data assimilation[END_REF] and [START_REF] Gejadze | On optimal solution error covariances in variational data assimilation problems[END_REF] the dynamic formulation is used. For example in the initial state control problem, the Hessian-vector product is defined by the successive solution of the tangent linear (TL) and adjoint models as follows

   ∂ψ ∂t -F ′ (ϕ)ψ = 0, t ∈ (0, T ), ψ| t=0 = v, (3.26)
   -∂ψ * ∂t -(F ′ (ϕ)) * ψ * = -(C ′ (ϕ)) * R -1 C ′ (ϕ)ψ, t ∈ (0, T ) ψ * | t=T = 0, (3.27)
H(u)v = B -1 v -ψ * | t=0 , (3.28)
where ϕ(u) is the solution of equation (2.12). These equations correspond to equation (3.21) with U = {u}. The dynamic formulation is important because it shows different implementation options. For example, if operator G is considered as a black box, then operators G ′ U (•) and (G ′ U (•)) * can be obtained by means of the Automatic Differentiation (AD) following the 'discretize-then-optimize' approach. In one hand, the forward, tangent linear and adjoint models obtained this way are mutually 'consistent'. This is a vital property required for minimization of realistic high-dimensional models. On the other hand, the models generated by means of the AD are often computationally cumbersome and should be optimized manually.

For computing H -1 one must compute the product Hv several times at the same origin point (model trajectory). Thus, the trajectory may be computed only once. In practice, computing the solution and the perturbation in the AD produced TL model is done 'simultaneously', i.e. one line of the TL solver follows the corresponding line of the direct solver. Moreover, in case of nonlinear iterations (to converge on ϕ), the TL solver follows step by step the iterative procedure, computing the perturbation variable by increments. This is highly inefficient. For computing the Hessian, however, it might be sufficient to have full consistency between the TL and adjoint models only. The forward and TL models can be partially consistent, i.e. consistent in terms of the discretized operator F (•), however the time integration method could be different. In this case the model's trajectory once computed and stored in memory can be supplied into the TL solver. Since the TL model is linear, the nonlinear treatment involved with the forward solver is not involved with the TLM. Thus, one can obtain a computationally efficient TL model. By applying the AD to the auxiliary cost-function (3.22), one can get the code which computes its gradient. The latter defines the product Hv according to (3.25). The code involves the adjoint model consistent with the modified TLM. Let us note that this adjoint may not be suitable for solving the DA problem itself.

For the parameter estimation (control) problem considered in [START_REF] Gejadze | On optimal solution error covariances in variational data assimilation problems[END_REF], the mathematical model is described by the evolution problem

   ∂ϕ ∂t = F (ϕ, λ) + f, t ∈ (0, T ) ϕ| t=0 = u, (3.29)
where u is the known initial condition, λ ∈ X p is an unknown model parameter and X p is the parameter space.

Let us introduce the functional

J(λ) = 1 2 R -1/2 (C(ϕ) -Y * ) 2 Y + 1 2 B -1/2 (λ -λ * ) 2 Xp , (3.30)
where λ * ∈ X p is a prior (background) of λ, and B is the covariance of δλ = λ * -λ, and consider the following DA problem with the aim to estimate the parameter λ: for given (u, f ), find λ ∈ X p and ϕ ∈ X such that they satisfy (3.29), and on the set of solutions to (3.29), the functional J(λ) takes the minimum value, i.e.

         ∂ϕ ∂t = F (ϕ, λ) + f, t ∈ (0, T ) ϕ| t=0 = u, J(λ) = inf v∈Zp J(v).
(3.31)

The Hessian-vector product is defined by the successive solution of the tangent linear and adjoint models as follows: That is, there exists a basic TL equation (system of PDE-based equations) which describes the evolution of perturbations through the system. However, the perturbations are applied at different entries: at the initial condition (case of the initial value control), and as a source term (case of parameter control). The adjoint equations coincide, the expression for the gradient expressed via the adjoint variable ψ * is different. This structure suggests an easy way of modifying the system (3.26)-(3.28) into (3.32)- (3.34). It is important to keep consistency between operator F ′ λ (•, •) and its adjoint. Since F ′ λ (•, •) is not acting on ψ, it may not necessarily be consistent with F ′ ϕ (•, •). All these details are important to design a computationally efficient TL and adjoint models for computing the Hessian-vector product (possibly not suitable for computing the gradient of J(U)).

   ∂ψ ∂t -F ′ ϕ ( φ, λ)ψ = F ′ λ ( φ, λ)v, t ∈ (0, T ), ψ| t=0 = 0, (3.32)    -∂ψ * ∂t -(F ′ ϕ ( φ, λ)) * ψ * = -(C ′ ( φ)) * R -1 C ′ ( φ)ψ, t ∈ (0, T ) ψ * | t=T = 0, (3.33)
H( λ)v = B -1 v -(F ′ λ ( φ, λ)) * ψ * . ( 3

On the role of the Hessian and its inverse

The importance of the Hessian matrix and its inverse in variational DA for geophysical applications is underlined in [START_REF] Thacker | The role of the hessian matrix in fitting models to measurements[END_REF], although this has been a well-known fact in statistics for decades (see, for example, [START_REF] Amemiya | Handbook of Econometrics[END_REF]). The previous Section illustrates that for linear and moderately non-linear DA/estimation problems H -1 (•) can serve as an approximation of the estimation (analysis) error covariance matrix. In particular, confidence intervals for the components of the estimate Û can be defined by the corresponding diagonal elements (variance) of H -1 ( Û ). A column c i of H -1 (•) which includes the i th diagonal element can be obtained by solving the equation H(•)c i = e i (where e i is a Euclidean unit vector). If the number of requested diagonal elements is significant, it would be much less expensive to evaluate H -1 (•) once and keep it in some limited-memory form, than to retrieve necessary diagonal elements using the Hessian-vector product rule.

In addition, H -1 (•) is involved in several other aspects of statistical pre-and post-processing and design of DA systems. Firstly, as the Hessian H(•) is equivalent to the Fisher information matrix (up to a constant multiplier), the diagonal elements of the inverse Hessian can also be used in the context of optimal experimental design involving such optimality criteria as loptimality [START_REF] Ucinski | Optimal Measurement Methods for Distributed-Parameter System Identification[END_REF], for example. Secondly, the analysis probability density function (pdf) is defined by the estimate Û and the estimation error covariance. Random functions from the Gaussian distribution N ( Û, H -1 ( Û )) can therefore be used as 'particles' of the ensemble of initial states, which may be useful for ensemble forecasting [START_REF] Toth | Ensemble forecasting at NMC: The generation of perturbations[END_REF]. These functions can be generated using U = U * + H -1/2 ( Û )ξ, where ξ ∼ N (0, I), or using the eigenvalues of H( Û) [START_REF] Ehrendorfer | Optimal prediction of forecast error covariances through singular vectors[END_REF]. However, in highly non-linear cases the 'particles' generated using this approach are unlikely to belong to the true posterior distribution, thus one must solve perturbed DA problems. This approach is referred to as the fully non-linear ensemble method [START_REF] Gejadze | On analysis error covariances in variational data assimilation[END_REF], or randomised maximum likelihood method [START_REF] Chen | Ensemble randomized maximum likelihood method as an iterative ensemble smoother[END_REF]. In these cases, an approximation of H -1/2 ( Û) can be used for preconditioning the non-linear minimization process to accelerate convergence, often with impressive results. Thirdly, the analysis error δU and the data errors ε and ξ are related via the approximate error equation (3.20), which can be used as a meta-model for investigating the effects of non-Gaussian data errors on the estimation error pdf. Lastly, if the model depends on parameters θ ∈ Θ, where Θ is the parameter space, an important problem is to evaluate the sensitivity of the estimation error to uncertainty ξ θ in these parameters. This can be done using the relationship

H(•)δU = (G ′ U (•)) * R -1 G ′ θ (•)ξ θ .
Once again, H must be inverted to obtain δU.

Another application is related to the use of H as a coefficient matrix (and preconditioner) in incremental 4D-Var [START_REF] Courtier | A strategy for operational implementation of 4d-var, using an incremental approach[END_REF]. Here each step of an outer iterative Gauss-Newton process is of the form U j+1 = U j + α j δU j , with a discrete approximation U j of the unknown initial state at iteration j, descent step α j and update (descent direction) δU j . As the update satisfies

H(U j )δU j = -G(U j ), (3.35)
where G(U j ) is the gradient of the cost function, a system of linear equations involving H has to be solved at each step. Given a Hessian-vector product evaluation routine, the systems in (3.35) are usually solved iteratively using, for example, the conjugate gradient (CG) algorithm.

An approximation of H -1 or H -1/2 , if available at a reasonable cost, can therefore be used to precondition equation (3.35) to accelerate convergence of this inner iteration.

Effective inverse Hessian method 4.1 Theory of the method

Here, we use the main results from [START_REF] Gejadze | Computation of the analysis error covariance in variational data assimilation problems with nonlinear dynamics[END_REF] and [START_REF] Shutyaev | Optimal solution error covariance in highly nonlinear problems of variational data assimilation[END_REF]. As mentioned above, the DCA error contains two components: linearization error component and the origin error component. First, we will see how the linearization component can be reduced.

Let us first consider a related abstract problem. Let g be a random vector of dimension n defined as follows

g = A(u, η)ξ,
where ξ is a Gaussian random vector of dimension n with zero mean and the covariance matrix B ξ , η(ξ) is a random vector with elements which are linearly or nonlinearly dependent on ξ, u is a given non-random vector and A is a n × n matrix with elements a i,j nonlinearly dependent on u and η. Let us assume that the matrix A can be presented in the form

A(u, η) = A(u, 0) + δA(u, η),
where A(u, 0) is the non-random part of A and δA(u, η) is the random part. The covariance matrix of the vector g is

P g = E (A(u, η)ξ -E[g])(A(u, η)ξ -E[g]) T . (4.36)
One can see that if η does not depend on ξ, i.e. vectors η and ξ are not correlated, then E[g] = 0, and

E A(u, η)ξξ T A T (u, η) = E A(u, η)E[ξξ T ]A T (u, η) .
For correlated vectors η and ξ (since we assume that η = η(ξ)), the following proposition can be put forward. Matrix

E A(u, η)E[ξξ T ]A T (u, η) = E A(u, η)B ξ A T (u, η
) is a better approximation to P g in (4.36) than A(u, 0)B ξ A T (u, 0), both in terms of the Frobenius and Riemann distances. The approximation accuracy increases when the correlation level between η and ξ decreases. The justification of this proposition is provided in [START_REF] Gejadze | Computation of the analysis error covariance in variational data assimilation problems with nonlinear dynamics[END_REF], and its validity has been analyzed numerically.

Let us consider equation (3.19) in the form

δU = H( Ū , δU, τ ) -1 (G ′ U ( Ū + δU)) * R -1 ξ + B -1 ε , (4.37)
where

H( Ū , δU, τ ) = (G ′ U ( Ū + δU)) * R -1 G ′ U ( Ū + τ δU) + B -1 . (4.38)
Since the nonlinear least-squares estimator is asymptotically unbiased, we assume that E[δU] is small and, therefore, for the covariance P we have an expression as follows:

P := E[δUδU T ] = E H -1 (G ′ U ) * R -1 ξξ T R -1 G ′ U (H -1 ) * + E H -1 B -1 εε T B -1 (H -1 ) * + E H -1 B -1 εξ T R -1 G ′ U (H -1 ) * + E H -1 (G ′ U ) * R -1 ξε T B -1 (H -1 ) * . (4.39)
As discussed in the beginning of this Section, we approximate the products ξξ T , εε T , ξε T and εξ T in (4.39) by E[ξξ T ] = R, E[εε T ] = B, E[ξε T ] = 0 and E[εξ T ] = 0, respectively. Thus, we write an approximation of P as follows:

P = E H -1 ( Ū, δU, τ ) (G ′ U ( Ū + δU)) * R -1 G ′ U ( Ū + δU) + B -1 (H -1 ( Ū , δU, τ )) * , ( 4
.40) where the expression in the round brackets is the Hessian (3.21) at the point Ū +δU. Therefore, the equation (4.40) can be written in the form

P = E H -1 ( Ū, δU, τ ) H( Ū + δU) (H -1 ( Ū , δU, τ )) * .
(4.41)

The value τ = 1/2 is the optimal one to achieve the best approximation of a difference

E[G( Ū + δU) -G( Ū)] by E[G ′ U ( Ū + τ δU)δU],
for δU ∼ N(0, P ). In this case, however, one must deal with operator H( Ū, δU, τ = 1/2), which is neither symmetric, nor positive definite, which may seriously complicate its eigenvalue analysis and the subsequent limited-memory representation. Besides, the double product formula is sensitive to the errors which result from the already accepted approximations. By assuming τ = 1 in H( Ū, δU, τ ) and keeping in mind that H( Ū , δU, τ = 1) := H( Ū + δU) this formula can be further simplified by the following expression:

P = E H -1 ( Ū + δU) . (4.42)
The right-hand side of (4.42) may be called the 'effective inverse Hessian', hence the name of the suggested method. In order to compute P directly using this equation, the expectation must be substituted by the ensemble mean:

P = 1 L L l=1 H -1 ( Ûl ), (4.43)
where Ûl are elements from the ensemble of estimates { Ûl }, l = 1, . . . , L, being obtained using perturbed data. Obviously, having such ensemble evaluated, the covariance matrix can be computed according to its definition as follows:

P = 1 L s Ls l=1 ( Ûl -Ū )( Ûl -Ū) T . (4.44)
The advantage of (4.43) is, however, that for L = L s it gives much better approximation of P than (4.44), particularly for small L s .

Key implementation details 4.2.1 Preconditioning

Implementation of formula (4.43) implies that L inverse Hessians have to be computed. This looks like an extremely laborious task, however, computational costs can be drastically reduced by preconditioning. We notice that the inverse Hessians are evaluated at the origin points perturbed around the true input value Ū , which means that the difference between H -1 ( Ū +δU l ) and H -1 ( Ū) can be described by just a few low-rank updates. Thus, we consider a projected Hessian in the form

H(•) = (B 1/2 ) * H(•)B 1/2 = (B 1/2 ) * (G ′ U (•)) * R -1 G ′ U (•)B 1/2 + I, (4.45)
and use H-1/2 ( Ū ) for the second-level preconditioning of H( Ū + δU l ) as follows:

H( Ū + δU l ) = H-1/2 ( Ū) H( Ū + δU l ) H-1/2 ( Ū ). (4.46)
One can see that H(•) is the Hessian of an auxiliary cost-function

J (•, δU) = 1 2 R -1/2 (G ′ U (•) H-1/2 ( Ū)B 1/2 δU -f 1 ) 2 Y + 1 2 δU -f 2 2 U . (4.47)
The limited memory approximation of H(•) is evaluated as a by-product of minimization of this cost-function by the LBFGS method, then, H -1 ( Ū + δU l ) is recovered by the formula

H -1 ( Ū + δU l ) = B 1/2 H-1/2 ( Ū ) H-1 ( Ū + δU l ) H-1/2 ( Ū)(B 1/2 ) * . (4.48)
Taking into account (4.43) and (4.48), the expression for the optimal solution error covariance P reads

P = B 1/2 H-1/2 ( Ū ) 1 L L l=1 H-1 ( Ū + δU l ) H-1/2 ( Ū)(B 1/2 ) * . (4.49)
The first step is, therefore, to find a limited-memory approximation of H-1 ( Ū ). This is done by minimizing the auxiliary cost-function (3.25) using the LBFGS algorithm. Next, having the product H-1 ( Ū)•v defined, one can evaluate the leading eigenpairs ({λ

(0) k , W (0)
k }, k = 1, . . . , K 0 ) of H-1 ( Ū ) using the Lanczos algorithm.

Given the leading (maximum magnitude) eigenpairs ({λ k , W k }, k = 1, . . . , K) of any symmetric operator A, the limited-memory representation of A in power β can be constructed as follows:

A β • v = I • v + K k=1 (λ β k -1)W k (W k ) * • v, (4.50)
In particular, for the square-root H-1/2 ( Ū ) we get:

H-1/2 ( Ū ) • v = I • v + K 0 k=1 ((λ (0) k ) 1/2 -1)W (0) k (W (0) k ) * • v, (4.51)
Let us note that at this factorization stage we do not run PDE models, so the Lanczos method is affordable. The product H-1/2 ( Ū) • v is involved in (4.47).

Quasi-random approach

The implementation of formulas (4.43) or (4.49) requires a set of estimates Ûl = Ū + δU l , l = 1, . . . , L to be computed. Each estimate is a solution to the original data assimilation problem for the cost-function (2.6) with perturbed data U * l = Ū + ξ l and Y * l = Ȳ + ε l , where ξ l and ε l are random events from ξ ∼ N(0, R) and ε ∼ N(0, B), respectively. Evaluating such a set could be fairly expensive.

Here, an alternative approach is suggested. If we denote by f δU the multivariate probability density of the estimation error δU, then the equation (4.42) can be re-written in the form

P = +∞ -∞ H -1 ( Ū + v)f δU (v)dv.
(4.52)

Since we assume that the estimation error is approximately Gaussian with zero mean and covariance P we obtain

P = 1 (2π) M/2 |P | 1/2 +∞ -∞ H -1 ( Ū + v) exp - 1 2 v T P -1 v dv. (4.53)
In contrast to formula (4.43), the above expression is a nonlinear matrix integral equation (deterministic) with respect to P , while v is a dummy variable. This equation can be solved by the fixed point iterative process

P k+1 = 1 (2π) M/2 |P k | 1/2 +∞ -∞ H -1 ( Ū + v) exp - 1 2 v T (P k ) -1 v dv, (4.54)
for k = 0, 1, ..., starting with P 0 = H -1 (ū).

The evaluation of the multi-dimensional integral in (4.54) using the quasi-Monte Carlo method means returning to formula (4.43). Taking into account (4.54), the iterative process takes the form

       P k+1 = 1 L L l=1 H -1 ( Ū + δU k l ), P 0 = H -1 (ū), k = 0, 1, ... (4.55)
where δU k l ∼ N(0, P k ). One can see that for each k the last formula looks similar to (4.43) with one key difference: δU k l in (4.55) is not an estimation error itself, but a vector having the statistical properties of this error. It is generated as follows:

δU k l = (P k ) 1/2 η l , η l ∼ N(0, I).
where η l is a quasi-random sequence from N (0, I). This approach is similar to the inflation approach used in statistical estimation methods (ensemble Kalman and particle filtering). If the preconditioning (4.45)-(4.46) is involved, the process (4.55) reads as follows:

       P k+1 = 1 L L l=1 H-1 ( Ū + δU k l ), P 0 = I, k = 0, 1, ..., (4.56)
where δU k l ∼ N(0, P k ), and the perturbation is defined as follows:

δU k l = ( P k ) 1/2 H-1/2 ( Ū )(B 1/2 ) * η l , η l ∼ N(0, I).
Finally, the covariance P in the original space is recovered by the formula

P = B 1/2 H-1/2 ( Ū ) P H-1/2 ( Ū )(B 1/2 ) * . (4.57)
5 Estimation error covariance versus Bayesian posterior covariance

General theory

Here, we use the main results from [START_REF] Gejadze | Analysis error covariance versus posterior covariance in variational data assimilation[END_REF], which point out the difference between the classical estimation (analysis) error covariance and the Bayesian posterior covariances, both in theoretical and computational aspects.

The estimation error covariance is associated with trying to find an approximation around the truth Ū , whereas the data is also assumed to come from the truth: Y * = G( Ū) + ξ, U * = Ū + ε, where ξ ∼ N (0, R) and ε ∼ N(0, B) are the observation and background error, respectively. The estimation error is defined as δU = U -Ū and its covariance is given by

P = E a [(U -Ū)(U -Ū) T ] = E a [δUδU T ],
where E a denotes averaging (expectation) with respect to the probability density function (pdf) ρ a which, taking into account the definitions of the data Y * and U * , can be defined as follows:

ρ a (U) = const • exp - 1 2 R -1/2 (G(U) -G( Ū)) 2 Y - 1 2 B -1/2 (U -Ū ) 2 U .
(5.58)

The Bayesian posterior covariance is given by

P = E p [(U -E p [U])(U -E p [U]) T ], (5.59)
where E p denotes averaging (expectation) with respect to the pdf

ρ p (U) = const • exp - 1 2 R -1/2 (G(U) -Y * 0 ) 2 Y - 1 2 B -1/2 (U -U * 0 ) 2 U .
(5.60)

In the above formula Y * 0 stands for the actual observations made by an instrument. It could be, in principle, seen in the form Y * 0 = G( Ū ) + ξ 0 , where ξ 0 is no longer a random variable, but a particular event (number) from N (0, R). However, this representation is not used in the Bayesian approach. Similarly, U * 0 stands for the actual background. Then, the estimate Û0 satisfies the equation

(G ′ U ( Û0 )) * R -1 (G( Û0 ) -Y * 0 ) + B -1 ( Û0 -U * 0 ) = 0. (5.61)
The posterior covariance is often approximated by trying to find the second moment of ρ p centered around Û instead of E p [U], which is natural because Û is the output of 4D-Var. In the has to be small as compared to K l because the difference between H( Ûl ) and H( Ûl ) is only due to the presence of the second-order term. Thus, the expenses of computing the eigenpairs {λ

(l) k , W (l) k } and { λ(l) k , W (l)
k } for given l can make only a fraction of those associated with computing {λ

(0) k , W (0) k }.
Given all eigenpairs are computed, the posterior covariance in (5.73) is evaluated as follows:

P = B 1/2 H-1/2 ( Û0 ) A( Ûl ) H-1/2 ( Û0 )(B 1/2 ) * , (5.76)
where

A( Ûl ) = 1 L L l=1 H-1/2 ( Ûl ) H-α ( Ûl ) H-1/2 ( Ûl ).
(5.77)

In (5.77), cases α = 2, α = 1 and α = 0 correspond to approximations (5.70), (5.71) and (5.72), respectively. The quasi-random implementation described in Sect.4.2.2 is applied to substitute Ûl by Û0 + δU l .

In the dynamic formulation the Hessian-vector product H(•)v is defined by the successive solutions of the following problems [START_REF] Gejadze | On analysis error covariances in variational data assimilation[END_REF]:

   ∂ψ ∂t -F ′ (ϕ)ψ = 0, t ∈ (0, T ), ψ| t=0 = v, (5.78)
   - ∂ψ * ∂t -(F ′ (ϕ)) * ψ * = (F ′′ (ϕ)ψ) * ϕ * -C * V 2 Cψ, t ∈ (0, T ) ψ * | t=T = 0, (5.79)
H(u)v = V 1 v -ψ * | t=0 .
(5.80)

Here ϕ and ϕ * are involved, being taken from (2.15)-(2.16). The problem (5.79) is the socalled second-order adjoint problem [START_REF] Dimet | Second-order information in data assimilation[END_REF]. One can see that it is a standard adjoint problem with a specially defined source term (F ′′ (ϕ)ψ) * ϕ * . For certain F , this term is relatively expensive to compute as compared to (F ′ (ϕ)) * ψ * , which could be a difficulty if the explicit time integration scheme is used. This is not an issue, however, when the implicit or semi-implicit schemes are used.

Conclusions

The main issue considered in this paper is the relationship between the Hessian of the costfunction and the estimation error covariance matrix in variational DA, whereas the discussion is focused on the essentially nonlinear case. In the classical (frequentist) statistical approach the estimation error δU is considered as a difference between the estimate Û and the 'truth' Ū . The corresponding covariance matrix can be approximated by the inverse of the Hessian H of an auxiliary cost-function computed at the 'truth' point. If the 'truth' is known (the identical twin experiment framework only), this is the first-order approximation with respect to δU. In practice, the Hessian is computed at the estimate Û0 , which is obtained after assimilating the real data. The difference between Û0 and Ū yields an additional error in the covariance (the origin error), which can not be removed in principle. Evaluating the inverse Hessian using the LBFGS algorithm with exact step search has been a novel computational approach.

In the Bayesian approach the estimation error δU is considered as a difference between the estimates Û and Û0 , conditioned on a perturbed and unperturbed data, respectively. Due to different centering of data, the first-order (with respect to δU) approximations of the estimation (analysis) error covariance and the Bayesian posterior covariance are different and should not be confused. The latter is computed via the Hessian-product formula (5.70), which involves both the Hessian of the original cost-function H, and the Hessian of the auxiliary cost-function H (often called the Jacobian). In the dynamic formulation the Hessian-vector product H(•)v is defined by the second-order adjoint problem. The formulas of the type (5.70) occasionally appear in the literature on statistics and nonlinear regression in the general context of nonlinear least-squares (see e.g. [START_REF] Donaldson | Computational experience with confidence regions and confidence intervals for nonlinear least squares[END_REF]).

Introducing the 'effective' Hessian-based covariance approximations (4.43) and (5.73), as well as the relatively feasible methods for their evaluation, is our main contribution in terms of novelty. The concept of such approximations stems from the exact error equations (3.18) and (5.65), the idea of which, in turn, has been first presented in [START_REF] Dimet | On deterministic error analysis in variational data assimilation[END_REF]. In practice, evaluating (4.43) or (5.73) requires an ensemble of the inverse Hessians H -1 ( Ûl ), or the products H -1 ( Ûl )H( Ûl )H -1 ( Ûl ) to be computed, however the size of such ensemble can be small enough (L = 25 -50). The suggested method can be regarded as an alternative to the standard sampling-based covariance evaluation method (involving localization and other possible tricks). The latter requires much larger ensemble size to achieve a similar approximation quality. The ways of improving the accuracy of the 'effective' Hessian-based covariance approximations are the matter of future research.

. 34 )

 34 These equations correspond to equation (3.21) with U = {λ}. From (3.26)-(3.28) and (3.32)-(3.34) one can understand the difference between G ′ u (•) and G ′ λ (•).

Author-produced version of the article published in Russian Journal of Numerical Analysis and Mathematical Modelling, 2018, N°33(1), p.25-39. The original publication is available at https://www.degruyter.com Doi: 10.1515/rnam-2018-0003

Acknowledgements

This work was carried out within the Russian Science Foundation project 17-77-30001(studies in Section 3), and the project 15-01-01583 of the Russian Foundation for the Basic Research.

Author-produced version of the article published in Russian Journal of Numerical Analysis and Mathematical Modelling, 2018, N°33(1), p.25-39. The original publication is available at https://www.degruyter.com Doi: 10.1515/rnam-2018-0003 linear Gaussian setup E p [U] and Û coincide. This is not true in general, but can be expected to be a good approximation if the volume of data is large and /or noise is small (due to the asymptotic properties of the nonlinear least-square estimator, see [START_REF] Amemiya | Handbook of Econometrics[END_REF], Vol. 1, Ch.6). Most importantly, due to different centering of Gaussian data, the Bayesian posterior covariance and the estimation error covariance are different objects and should not be confused. However, they are equal in the linear case.

In order to simulate the distribution (5.60) we define perturbed data as follows: Y * = Y * 0 +ξ and U * = U * 0 + ε. Then, the perturbed estimate Û satisfies equation

(5.62)

One way to derive the error equation is to subtract (5.61) from (5.62). For the difference in the residual terms we obtain

Using the Taylor-Lagrange formula [START_REF] Marchuk | Adjoint Equations and Perturbation Algorithms in Nonlinear Problems[END_REF] the expressions in the brackets take the form

For the difference in the penalty terms we obtain

(5.64) Combining (5.63) and (5.64) we obtain the error equation as follows:

In the operator form this equation can be rewritten as

where

Let us note that, generally, H is neither symmetric, nor positive definite. The estimation error is expressed from (5.66) as follows: Since operators H( Û, Ũ1 , Ũ2 ) and G ′ U ( Û) nonlinearly depend on δU, it is not possible to consider them as constant multipliers when applying the expectation operator. One way is to fix Û , Ũ1 , Ũ2 at the value Û0 . In this case H(•) becomes the Hessian of the cost-function (2.6) at the origin point Û0 . Next, by applying (5.69) we obtain

(5.70) (compare this equation to (3.24)).

Since the condition number of the double product H -1 HH -1 is usually larger than the condition number of its components H and H, the above formula could be quite sensitive to the approximation errors. One way is to simplify this formula assuming H and H are not very different. These simplification is

and, subsequently,

Another approach is to use the 'effective inverse Hessian' method presented in Sect. 4. In this case equation (5.70) becomes

H -1 ( Ûl )H( Ûl )H -1 ( Ûl ).

(5.73)

Implementation

As before, the computationally efficient and, therefore, feasible implementation is achieved using preconditioning. The first-level preconditioning yields

and the second-level preconditioning yields

(5.75)

At the first step we compute the leading eigenpairs ({λ

k }, k = 1, . . . , K 0 ) of H( Û0 ). This allows H-1/2 ( Û0 ) to be defined in the limited-memory form. Then, inside the ensemble loop (index l) we compute: a) the leading eigenpairs ({λ

k }, k = 1, . . . , Kl ) of H( Ûl ). For a reasonably small δU l = Ûl -Û0 , the number K l of eigenpairs describing the Riemann distance between H( Ûl ) and H( Û0 ) has to be small as compared to K 0 . Similarly, Kl

Author-produced version of the article published in Russian Journal of Numerical Analysis and Mathematical Modelling, 2018, N°33(1), p.25-39. The original publication is available at https://www.degruyter.com Doi: 10.1515/rnam-2018-0003


Thierry Sentenac
email: sentenac@mines-albi.fr

Rémi Gilblas

Noise effect on the interpolation equation for near infrared thermography

This paper investigates the performance of interpolation equations for a near infrared thermal imager operating over wavelengths from 0.9 µm to 1.7 µm with various filter bandwidths and a broad temperature range from 300 • C to 1000 • C. The equations are based on a general formulation of the effective wavelength as a function of the temperature. The quality of the interpolation is assessed in relation to the order of the effective wavelength. However, the noise induced by the imperfections of the thermal imager significantly disturbs the signal, and this phenomenon is enhanced as the bandwidth of the filter increases (i.e. for low-temperature applications). The main purpose of this paper is to establish the right choice of the filter bandwidth and the expression and order of the interpolation equation in relation to the noise level on the thermal imager and the desired accuracy. This paper first outlines the background on interpolation equations and then tests them on synthetic data from signals delivered first by an ideal thermal imager (i.e. free from noise) and then from noisy signals. This simulation study provides a framework for users to select an interpolation equation with an adequate order for near infrared thermal imagers. The performances of the selected interpolation equations are finally demonstrated on real images performed by a near infrared thermal imager.

Introduction

An accurate on-line temperature measurement is the key parameter for the optimization and quality assurance of many industrial processes. Non-contact infrared systems, such as pyrometers and thermal imagers, are suitable in many industrial applications or research projects.

For temperatures ranging from 300 • C to 1000 • C, near infrared thermal imagers, operating at short wavelengths (from 0.9 µm to 1.7 µm), provide high performances in terms of temperature sensitivity, and also have the advantage of being less sensitive to emissivity. Furthermore, they deliver images with a high spatial resolution which is suitable for studies on materials with non-uniform surfaces. These thermal imagers are then widely used for thermal imaging applications [START_REF] Mayrhofer | Optical characterization of thin layers grown on metal components[END_REF], for radiance temperature measurements [START_REF] Rotrou | Thermographie courtes longueurs d'onde avec des caméras silicium: contribution à la modélisation radiométrique PhD Thesis Ecole nationale supérieure de l'aéronautique et de l'espace[END_REF], for multi-spectral measurements [START_REF] Mansouri | Optical calibration of a multispectral imaging system based on interference filters[END_REF][START_REF] Meriaudeau | Real time multispectral high temperature measurement: application to control in the industry Image[END_REF] and for true temperature field measurements [START_REF] Sentenac | Bi-color near infrared thermoreflectometry: a method for true temperature field measurement[END_REF].

For these applications, it is essential that a radiometric calibration of the thermal imager is carried out with low temperature error [START_REF] Ferrero | Low-uncertainty absolute radiometric calibration of a CCD[END_REF]. The signal of the thermal imager is the integral of its spectral responsivity multiplied by the spectral radiance of the black body. An ideal radiometric calibration requires the measurement of the spectral responsivity of the thermal imager and the integration of Planck's law over its spectral bandwidth. Unfortunately, the spectral responsivity provided by the manufacturer is often given as a relative value.

When the spectral responsivity of the thermal imager is unknown, an algebraic interpolation equation approximates the relationship between the output signal of the thermal imager and the temperature of the black body. The calibration parameters of the interpolation equation are estimated by a common multi-point calibration method based on leastsquares fitting of the measured output signals versus the different temperatures of the black body. Then, the difference between the estimated temperature obtained by inverting the interpolation equation and the true black-body temperature is known as the interpolation error. It should ideally be smaller than about 1/3 to 1/5 of the quadrature sum of all the other uncertainty components of the radiometric calibration.

Hence, it should ideally be comparable to the accuracy of the black-body temperature (around 0.1 • C). Therefore, the interpolation equation should not introduce additional errors and the interpolation error should not impact the total uncertainty. In this secondary radiometric calibration approach to thermal imagers, performed with a multi-point method using interpolation equations, choosing a proper interpolation equation is a fundamental step to minimize the interpolation error.

The easiest interpolation equation is the monochromatic approximation of the product of black-body radiance in Planck or Wien formulation by the thermal imager responsivity around a mean operating wavelength. The suitability of this equation depends on the temperature range, the spectral bandwidth of the thermal imager and the required accuracy of the application. As shown in [START_REF] Kostkowski | Theory and Methods of Optical Pyrometry[END_REF], the monochromatic approximation introduces large interpolation errors in the near infrared spectral band, which is not acceptable for many applications. This error can be reduced by adding correction factors used with a reference wavelength [START_REF] Bezemer | Spectral sensitivity corrections for optical standard pyrometers[END_REF] and at high temperatures (above 1500 • C) [START_REF] Hahn | Interpolation equation for the calibration of infrared pyrometers[END_REF]. However, the estimated temperature calculated from this interpolation equation cannot be derived analytically. In contrast, the Sakuma-Hattori [START_REF] Sakuma | Establishing a practical temperature standard by using a narrow-band radiation thermometer with a silicon detector Temperature: its Measurement and Control in Science and Industry[END_REF] interpolation equation, based on a temperature-dependent effective wavelength, which is a polynomial function at first order of the temperature, is analytically invertible. However, this equation is fitted using non-linear least-squares techniques. An interpolation equation with a general formulation of the inverse effective wavelength as a polynomial function at various orders of the inverse temperature, suggested in [START_REF] Saunders | General interpolation equations for the calibration of radiation thermometers[END_REF], can easily be fitted using linear least-squares techniques when this interpolation equation is expressed with the Wien approximation. Moreover, a physical interpretation of the calibration parameters is provided.

This paper addresses the issue of the choice of an interpolation equation with this final formulation of a temperature-dependent effective wavelength for a near infrared thermal imager equipped for various filters representing different applications of near infrared thermography and operating over a broad temperature range (300 • C-1000 • C). With respect to [START_REF] Saunders | Interpolation errors for radiation thermometry[END_REF], the paper assesses the interpolation errors for an ideal thermal imager (i.e. free from noise) to study the error due to the different orders of the effective wavelength and the influence of the filter characteristics. In addition, the interpolation error is investigated from the interpolation equation fitting to noisy data. The interpolation error is then the addition of an error due to the interpolation equation (evaluated from noiseless data) and an error due to the noise in the thermal imager signal. The main objective of this paper is to evaluate the performance of the interpolation error in relation to the level of noise.

This paper is organized as follows: in section 2 the main interpolation equations required for the study are recalled. In section 3, interpolation errors are calculated for simulated data from an ideal near infrared thermal imager equipped with Gaussian filters. The interpolation errors are also given in relation to the noise of the thermal imager signal. Finally, section 4 is devoted to the assessment of the performances of interpolation equations with real data from a near infrared thermal imager with a different bandwidth filter.

Background to interpolation equations

This section recalls the fundamental radiometric equation which links the black-body temperature to the thermal imager's signal through the integration of Planck's law and the spectral responsivity of the thermal imager. When this spectral responsivity is unknown, interpolation equations provide a useful alternative by approximating the integral using an analytic function with fitting parameters.

This section discusses various definitions of interpolation equations.

Fundamental radiometric equation

Let us consider a thermal imager, in a dark room, viewing a black-body source at temperature T , with an atmosphere of unitary transmission and without any other heat sources. Then, the output signal I D , provided by the thermal imager with a linear response, is the integral of the spectral responsivity (W (λ)) multiplied by the spectral radiance of the black body (defined on a finite bandwidth λ):

I D = k λ L 0 (λ, T ) W (λ) dλ, ( 1
)
where k is a constant dependent on geometrical, optical, electrical and digital properties of the thermal imager which is expressed in DL m 2 sr W -1 (DL: digital level). This constant can be computed up to a scale factor and for this study it will be set at one. The integration bandwidth λ depends on the thermal imager. The relative spectral responsivity W (λ) includes the detector sensitivity and the transmittances of the filter and lens. The spectral radiance of the black body, L 0 , is given by Planck's law at wavelength λ and temperature T . According to equation ( 1), an ideal radiometric calibration requires the measurement of the spectral responsivity W of the thermal imager and the integration of Planck's law over this spectral responsivity [START_REF] Campos | Spectral responsivity uncertainty of silicon photodiodes due to calibration spectral bandwidth[END_REF]. However, its measurement remains difficult to obtain.

Principle of the interpolation equation

As a first approach the interpolation equation can be expressed as a monochromatic approximation, based on the mean operating wavelength of the thermal imager (which is the mean wavelength value of the product of filter transmittance and imager responsivity), λ 0 , as follows:

k λ L 0 (λ, T ) W (λ) dλ ≈ ÎD = k λ L 0 (λ 0 , T ) W (λ 0 ). (2)
The next step can be the choice between Planck's form and Wien's approximation for the expression of the black-body spectral radiance L 0 (λ 0 , T ). Wien's approximation is accurate to within 1% if the condition λ T < 3000 µm K is observed. This choice depends on the temperature range and the spectral band of the thermal imager. For the temperature range 300 • C-1000 • C, considered in our application, the wavelength should be shorter than 2.4 µm. Wien's approximation is then available for our wavelengths (0.9 µm to 1.7 µm) and the simplest interpolation equation can be expressed as follows:

ÎD = A exp -C 2 λ 0 T , (3)
where A is a fitting parameter related to the spectral responsivity given by

A = k C 1 λ W (λ)λ -5 dλ ≈ k λ C 1 λ -5 0 W (λ 0 ), (4)
C 1 and C 2 are the first and the second radiation constants (C 1 = 1.191 042 × 10 8 W µm 4 m -2 and C 2 = 14 388 K µm).

Parameter A is assumed to be a property of the thermal imager that depends only on its spectral characteristics and responsivity. The third step is the calculation of the interpolation error, which is the difference between the estimated temperature provided by the interpolation equation and the black-body temperature. This error is compared with the required accuracy of the application. For example, in the near infrared spectral band, with Wien's formulation, this monochromatic interpolation equation provides an interpolation error greater than 1 • C [START_REF] Rotrou | Near infrared thermography with silicon FPA-comparison to MWIR and LWIR thermography Quant[END_REF] in the temperature range 300 • C-1000 • C. A correction factor is then usually applied to decrease the interpolation error. A first method is the reference wavelength method [START_REF] Coates | Wavelength specification in optical and photoelectric pyrometry[END_REF][START_REF] Hahn | Interpolation equation for the calibration of infrared pyrometers[END_REF] which uses a temperature-dependent factor A(T ) and a fixed reference wavelength, λ r , in the exponential term (λ r can exhibit a value close to λ 0 ). However, the estimated temperature calculated from this interpolation equation cannot be derived analytically but only by iteration. In contrast, the effective wavelength method [START_REF] Sakuma | Establishing a practical temperature standard by using a narrow-band radiation thermometer with a silicon detector Temperature: its Measurement and Control in Science and Industry[END_REF] is based on a fixed parameter A and a temperature-dependent wavelength, λ x (T ). Some forms of the interpolation equation provided by this method are analytically invertible to calculate the estimated temperature. The effective wavelength method is thus widely used and is detailed in the next paragraph.

Effective wavelength method

The effective wavelength method uses a varying wavelength, called the extended effective wavelength, and is denoted by λ x , which is calculated as a function of temperature. The interpolation equation based on Wien's law is given by the following equation:

ÎD = A exp -C 2 λ x (T ) × T , (5)
where

1 λ x = T C 2 ln λ C 1 W (λ)λ -5 dλ λ W (λ)L 0 (λ, T ) dλ .
Without the knowledge of W , an analytical extended effective wavelength versus the temperature should be suggested. Reference [START_REF] Sakuma | Establishing a practical temperature standard by using a narrow-band radiation thermometer with a silicon detector Temperature: its Measurement and Control in Science and Industry[END_REF] suggests an extended effective wavelength which is expressed as a first order polynomial function of the inverse temperature. The Sakuma-Hattori interpolation equation is then given by

ÎD = A exp -C 2 b 0 × T + b 1 , (6)
where b 0 and b 1 are fitting parameters.

Reference [START_REF] Saunders | General interpolation equations for the calibration of radiation thermometers[END_REF] generalizes this formulation and the inverse of the extended effective wavelength is a polynomial function of the inverse temperature with different orders, as follows:

1 λx = N i=0 a i T i (7)
where a i are fitting parameters and N is the order of the polynomial function, which depends on the operating wavelength of the thermal imager, the spectral bandwidth of the filter, the temperature range and the required accuracy of the application. The interpolation equation formed using equations ( 5) and ( 7) is difficult to invert for high orders.

However, when used with Wien's approximation, it can be fitted, to all orders, using linear least-squares techniques, unlike the Sakuma-Hattori equation which requires non-linear least-squares techniques. As shown in [START_REF] Saunders | Physical basis of interpolation equations for radiation thermometry[END_REF], for relatively narrow bandwidth thermal imagers, it is sufficient to use a polynomial function at first order (N = 1 and a 2 = 0) and the parameters a 0 and a 1 have a physical meaning and can be well approximated by equations ( 8) and ( 9) as follows:

a 0 = 1 λ 0 1 + 6 σ f λ 0 2 ( 8
)
a 1 = - C 2 2 λ 2 0 σ f λ 0 2 , (9)
where λ 0 and σ f are the mean wavelength and standard deviation of the spectral responsivity of the thermal imager.

If the ratio r = σ f λ 0 ≪ 1, the parameter a 0 will be close to the inverse of the mean wavelength of the spectral responsivity. The term a 1 describes the temperature dependence of the extended effective wavelength. The parameters a 0 and a 1 depend only on the mean wavelength and standard deviation of the spectral responsivity and they are independent of its shape (Gaussian, uniform, etc). Reference [START_REF] Saunders | Interpolation errors for radiation thermometry[END_REF] provides a comprehensive study of the influence of r on the interpolation error. For wide bandwidth thermal imagers, the order N of equation ( 7) can be higher and at minimum the term a 2 is required.

Simulation of the interpolation error

A simulation study of the interpolation error is carried out with an interpolation equation ( 5) based on an extended effective wavelength equation [START_REF] Kostkowski | Theory and Methods of Optical Pyrometry[END_REF] expressed as first and second order. The interpolation error is calculated, in the case of thermal imagers operating in NIR, for narrow and large bandwidths and for a wide temperature range from 300 • C to 1000 • C. The first aim is to highlight the variation of interpolation error in relation to the extended effective wavelength's order and the wide variation of the relevant parameter r = σ f λ 0 , as in [START_REF] Saunders | Interpolation errors for radiation thermometry[END_REF]. The second goal is to understand the influence of the noise due to experimental errors (non-linearity, non-uniformity, temporal noise, etc) on the interpolation error. This simulation analysis is carried out with noiseless and noisy data.

Simulation procedure

The simulated output signal I D of the thermal imager is obtained according to equation ( 1) and therefore by integrating the product of the spectral responsivity of our thermal imager (W ) with Planck's law. The parameter W is equal to the product of the spectral responsivity of the detector (set between 0 and 1), provided by the manufacturer, with the spectral distribution of the filters. The filters are Gaussian with a mean wavelength of λ 0 = 1310 nm and different full-width at half-maximum (FWHM), from 10 nm to 600 nm (which corresponds to a standard deviation σ f = FWHM √ 8 ln 2 varying from 4.24 nm to 254.8 nm). The relevant selected parameter, the ratio r = σ f λ 0 , varies from 3.24 × 10 -3 to 0.194. Finally, equation ( 1) is calculated over the temperature range from 300 • C to 1000 • C (respectively for a bandwidth λ from 900 nm to 1700 nm).

For the second simulation, noise is optionally added to the output signal by defining ǏD j = I D j + ν j , where ν is the j th trial value of M trials (j = 1 . . . M) for a random variable following a Gaussian distribution N (0, σ 2 I D ) of mean zero and standard deviation σ I D .

The parameters of interpolation equation ( 5) are fitted to the previous simulated output signal, I D or ǏD , and to the reference temperature values T of Planck's law. The fit method is based on weighted least-squares fit where the weights are inversely proportional to the output signal values.

Finally, the estimated temperature T is calculated by the inverse of equation ( 5) and then by solving the polynomial equation ( 10) with Cardan's formula at third order.

a 2 T 3 + a 1 T 2 + a 0 T - 1 C 2 (ln( ÎD ) -ln(A)) = 0. ( 10
)
The absolute interpolation error E T = |T -T | and the maximum absolute interpolation error over all temperatures, E T peak , are then calculated. For noisy simulations, the mean value ( ĒT ) and the standard deviation (σ E T ) of interpolation errors are calculated over M trials. The mean value over all temperatures is denoted ¯Ē T .

Influence of extended effective wavelength order

This first series of simulations is intended to study the influence of the extended effective wavelength's order in relation to the relevant parameter r = σ f λ 0 . The simulation is carried out without noise on the output signal I D . This approach clearly defines the dependence of the interpolation error on the order of the extended effective wavelength.

Figure 1 shows that the maximum absolute interpolation error, E T peak , decreases when the order of extended effective wavelength increases and when the ratio r decreases. At order zero, an interpolation error lower than 0.1 • C is achieved with a value of r lower than 9 × 10 -3 , which corresponds to a FWHM of about 35 nm. At first order, the same value of the interpolation error is achieved for a value of r lower than 7 × 10 -2 , or a FWHM value around 220 nm. Finally, at second order, the interpolation error is lower than 0.1 • C for all the spectral bandwidths of the thermal imager. As shown in table 1, and according to equation ( 4), the fitting parameter A increases with the bandwidth of the spectral responsivity of the camera. The increase, in absolute value, of the parameters a 0 and a 1 is also consistent with equations ( 8) and [START_REF] Hahn | Interpolation equation for the calibration of infrared pyrometers[END_REF]. The difference between the fitting values and the theoretical values remains low (a few per cent) as long as the value of r is lower than 0.05. The value of parameter a 2 is small with a narrow bandwidth filter and increases quickly with a wide bandwidth filter.

As a first conclusion and with noiseless data, a high order always decreases the interpolation error, and this decrease is even more significant for large values of r. For wide bandwidth filters, an interpolation equation with an extended effective wavelength at second order is the best choice to achieve an interpolation error of less than 0.1 • C.

Influence of noise

The second series of simulations are performed to study the influence of noise on the choice of the extended effective wavelength order.

The noisy output signal, ǏD , is generated with a variation of relative standard deviation σ I D from 0% to 2% (for an output signal digitalized on 12 bits, the maximum relative standard deviation represents 80 DL). The number of trials is set to M = 100. This simulation is only performed with a wide bandwidth filter and with a value of r of 0.194 (for a mean wavelength of λ 0 = 1310 nm and a Gaussian filter, it corresponds to a FWHM value around 600 nm). This case induces a maximum difference of the interpolation error between an extended effective wavelength at first and second order.

For each M value of ǏD , the temperatures are estimated with an interpolation equation with an extended effective wavelength at first and second order. Figure 2 shows the mean value of the interpolation error over the trials and temperatures, called ¯Ē T , for both orders. ¯Ē T is very different when the relative standard deviation σ I D is lower than 1.5%. After this value, the quantity ¯Ē T becomes similar for both cases.

E T peak / E T peak / A/ a 0 / a 1 / a 2 / E T peak / r • C • C (DL s -1 ) µm -1 (K µm -1 ) (K 2 µm -1
A second-order extended effective wavelength significantly reduces the interpolation error only if the noise standard deviation of the output signal is lower than 1.5%. With a very noisy output signal, when the relative standard deviation σ I D is above 1.5%, an increase in the order of the extended effective wavelength moderately improves the interpolation error. As a second conclusion, for low-noise data, the largest part of the interpolation error is due to the interpolation equation. An extended effective wavelength at second order significantly reduces the interpolation error. In contrast, for high-noise data, a second order does not seem to provide improvements. The next paragraph studies the contribution of the second-order extended effective wavelength for high-noise signals.

Coupled influence of order and noise

The third series of simulations extends the previous results to the entire range of variation of parameter r (from 3.24 × 10 -3 to 0.194) and with a variation of relative standard deviation σ I D of the output signal from 0% to 1%. The number of trials is set at M = 100.

Figures 3 and4 confirm that the quantity ¯Ē T increases with the ratio r and with the noise of the output signal. The mean interpolation error obtained with an extended effective wavelength at first or second order is always greater than 0.1 • C for a relative standard deviation higher than 0.15%, whatever the value of r. Above this value, the interpolation error is mainly due to the noise in the output signal.

However, the quantity ¯Ē T increases faster with the ratio r for an extended effective wavelength at first order than it does at second order. It reaches 0.79 • C at order 1 for a relative standard deviation σ I D = 1% and for a ratio r = 0.194. For the same configuration, the quantity ¯Ē T for an extended effective wavelength at second order rises to 0.59 • C. Figure 5 demonstrates that the difference between quantity ¯Ē T for an Figure 5. Difference of interpolation errors averaged over the trials and temperatures for an extended effective wavelength at first and second order versus the ratio r of a Gaussian filter (mean wavelength of λ 0 = 1310 nm and FWHM from 10 nm to 600 nm) and for six levels (from 0% to 1%) of the relative standard deviation of the noise of the output signal.

interpolation equation with an extended effective wavelength at first and second order, denoted ¯Ē T = ¯Ē T (order 1) -¯Ē T (order 2), is always positive.

As a third conclusion, an interpolation equation with an extended effective wavelength at second order is always more efficient than it is with one at first order. As this positive difference increases with the ratio r, an extended effective wavelength at second order is always advisable with wide bandwidth filters. However, this benefit decreases when the noise standard deviation on the output signal increases. For example, for the highest values of parameter r = 0.194, the second order brings an additional accuracy of about 0.43 • C compared with the first order. When the noise reaches 1%, this performance decreases to only 0.2 • C. For very noisy signal, and for applications involving high time resolution, a secondorder parametrization of the extended effective wavelength may not be the best compromise.

Experimental validation of the interpolation error

To compare simulation with experimental results, the interpolation equation with an extended effective wavelength at first or second order is fitted from measured data provided by a near infrared thermal imager.

The calibration methodology, involving a black body and the thermal imager, is described in the first paragraph. Compared with the simulation step, the thermal imager can be seen as a matrix of detectors and each detector as a thermometer. Then, the interpolation equation is fitted from the output signal of each detector. However, this signal is not ideal and some departures from ideal behaviour are characterized in the second paragraph. Finally, the interpolation error of each detector is determined in two configurations of calibration: the thermal imager equipped with a narrow Gaussian filter and without it (i.e. on the entire spectral band of the thermal imager). The mean value of this interpolation error field is compared with simulation results.

Experimental set-up

The calibration set-up is composed of a thermal imager aligned with a black body on a translation stage.

The easiest way to calibrate the thermal imager is to illuminate it uniformly so that each pixel receives the same flux. A flat-plate calibrator or a system based on an integrating sphere would offer such an illumination compatible with the field of view of the imager. Unfortunately, the flatplate calibrator cannot provide a reference temperature up to 500 • C, which is necessary to calibrate the near infrared thermal imager, and an integrating sphere would complicate the experimental set-up. The calibration of the thermal imager is then carried out with a cavity black body with relatively small aperture size in relation to the field of view of the thermal imager. It provides reference temperatures over the range 50 • C to 1200 • C with a cavity diameter of 50 mm and an emissivity of 0.995. The NIR lens used has a 50 mm focal length and the distance between the thermal imager and the black body is around 1 m. The advantage of this calibration is the reliable, high and constant emissivity of the cavity black body. The major drawback is that the aperture only fills a relatively small part of the image (25% or 20 000 pixels). In order to assess the entire image it is necessary to carry out multiple measurements with the source positioned at different points within the image. In this paper, only one position is shown. The temperature of the black body is measured by a reference S-type thermocouple. The experimental set-up is located in a black room to minimize surrounding radiation.

The thermal imager is a 320×256 InGaAs array operating in the near infrared spectral band from 0.9 µm up to 1.7 µm and with a Peltier cooler. The set-up includes a filter wheel with three Gaussian filters, a black screen and an empty hole. The wheel is motorized to rotate and is piloted by the software. For this experiment, the mean wavelength is λ 0 = 1310 nm with a FWHM of 50 nm for the narrow filter considered and a FWHM of 600 nm without filters. The thermal imager is equipped with a lens of focal length 50 mm and the F-number is set to 4. The exposure time range available is between 1 µs and 400 s. However, the minimum time is limited to 100 µs to minimize the non-linearity of the signal with the exposure time, and the maximum time is limited to 1 s to minimize the duration of the calibration. Consequently, the calibration temperature range is limited between 500 • C and 850 • C, which is lower than the range used for the simulations. As shown in reference [START_REF] Saunders | Interpolation errors for radiation thermometry[END_REF], for the Sakuma-Hattori equation and when the temperatures are equally spaced, the peak interpolation error increases as the cube of the temperature range. The values of the simulated and measured interpolation errors will not be the same. However, trends are comparable, which validates the noise impact on the interpolation equation.

Evaluation of thermal imager signal

The thermal imager is an array device composed of many detectors. Its behaviour is a composition of the individual behaviour of each detector, and of their juxtaposition.

As a single detector, the output signal of each detector is given by the sum of three terms: signal due to the radiation source, noise signal from the dark current and offsets and the signal produced by temporal noises. The ambient radiation entering the black body can be neglected above 500 • C.

The dark current and the offsets are part of the signal acquired in darkness. This signal, denoted I dark D , can be subtracted from the signal obtained during the illumination with the radiation source, denoted I D .

Moreover, the exposure time ti is managed to avoid the thermal imager signal becoming lower than the noise level and to avoid overexposure. When a detector is illuminated with sufficient radiance, the exposure time is then calculated to obtain a signal level around 75% of the maximum value. Around this value and for the limited exposure time chosen, the signal of each detector is linearly dependent on the exposure time.

Under these conditions, the corrected signal for each detector of coordinates u and v, denoted I c D (u, v), is calculated from the difference between the signal under illumination, I D , and the dark signal, I dark D , and is normalized by the exposure time, as given in equation [START_REF] Saunders | General interpolation equations for the calibration of radiation thermometers[END_REF].

I c D (u, v) = I D (u, v) -I dark D (u, v) ti ( 11
)
The effects of temporal noise are minimized by averaging the corrected signal over several frames (in practice 50 frames). The histogram plotted in figure 6 shows the dispersion of the corrected signal for a central pixel. The temporal mean Ī c D and its associated standard deviation, denoted σ I c D , are calculated. The values are recorded for a black-body temperature of T = 500 • C and without filter. Considering the noise as independent of the temperature, this measurement gives the minimum signal-to-noise ratio and represents the worst situation. The standard deviation of the corrected signal is low (around σ I c D = 0.5%), and therefore the mean value represents an accurate representation of each detector's value. It should be noted that this standard deviation, with a λ 0 = 1.3 µm monochromatic thermal imager, corresponds to a temperature uncertainty of T = T 2 λ 0 C 2 σ I c D = 0.27 • C at a temperature of T = 500 • C. As our thermal imager is not monochromatic, the measured error would be higher than this value.

As a device array, additional factors need to be investigated, such as image non-uniformity and size-of-source effect (SSE). Firstly, the calculation of the corrected signal compensates the non-uniformity due to dark current and offsets. Secondly, the residual non-uniformity is due to the difference in responsivity of each detector, denoted W (u, v) in equation ( 1). This can be taken into account by a calibration of each detector across the calibration parameter A. To do so, the responsivity of each detector has to be assumed to be independent and therefore the SSE and the cross-talk problem must be negligible. It should be noted that this calibration is performed under the same conditions as our application [START_REF] Sentenac | Bi-color near infrared thermoreflectometry: a method for true temperature field measurement[END_REF], for which the temperature measurement is carried out on a target with the same size as the black-body aperture and at the same distances. These two phenomena are therefore assessed over a range of distances in the depth of field of the thermal imager and with the image well focused. The central area is then well overfilled by the aperture of the black body at a fixed temperature of 700 • C. Under these conditions of a small change in angular black-body aperture, the signal indicated by the detectors within the central area of the black body slightly increases as the black-body distance decreases. The variation in the level of mean value of the corrected value is less than 0.5%. As a consequence, correction of SSE is not applied, and the imperfections of optical components and the crosstalk problem can be neglected. The response of each detector is then assumed to be independent.

Table 2. Mean spatial and standard deviations values of fitting parameters with an extended effective wavelength at first and second order for a narrow bandwidth filter with a FWHM of 50 nm and a mean wavelength of λ 0 = 1310 nm. As a conclusion, the radiometric calibration of the thermal imager is equivalent to the calibration of each detector. Each detector produces a corrected signal that can be calibrated independently of the others, and all can be calibrated together at the same time. The calibration signal of each detector is the temporally averaged signal, denoted Ī c D :

1 λx Ã/(DL s -1 ) 100 × σ A ã0 /µm -1 100 × σ a 0 ã1 /(K µm -1 ) 100 × σ a 1 ã2 /(K 2 µm -1 )
Ī c D (u, v) = A exp C 2 λ x (T ,u,v)×T (12) with 1 λ x (T , u, v) = a 0 (u, v) + a 1 (u, v) T + a 2 (u, v) T 2 . ( 13
)
The calibration depends on three or four images of fitting parameters, A(u, v), a 0 (u, v), a 1 (u, v) and a 2 (u, v), that are determined during the radiometric calibration process.

Interpolation error of the thermal imager with a Gaussian filter with a FWHM of 50 nm

For a thermal imager equipped with a Gaussian filter with a FWHM of 50 nm and a mean wavelength of λ 0 = 1310 nm, the first calibration procedure is carried out measuring the corrected signal from equation [START_REF] Saunders | General interpolation equations for the calibration of radiation thermometers[END_REF] and knowing the black body temperatures ranging from 500 • C to 850 • C with an increment of 25 • C. Hence, three or four images of the calibration parameters (A(u, v), a 0 (u, v), a 1 (u, v) and a 2 (u, v)) of interpolation equation ( 12) with an extended effective wavelength at first and second order are fitted. The parameter values are then computed for each detector of coordinates (u, v) in the area illuminated by the black body. It should be noted that the detectors outside the area of the black body are not calibrated. The mean spatial value, denoted Ã, ã0 , ã1 and ã2 and the standard deviations, denoted σ A , σ a 0 , σ a 1 and σ a 2 of calibration parameters are tabulated in table 2.

The relative standard deviation of the parameter A is around 2%. This value shows the low non-uniformity of the detector's responsivity (see equation ( 4)). For the calibration parameters ã0 and ã1 , the experimental values are close to the theoretical values (see table 1). According to equation [START_REF] Bezemer | Spectral sensitivity corrections for optical standard pyrometers[END_REF], the inverse of the parameter ã0 is homogeneous to the mean wavelength of the filter, λ 0 (when σ λ 0 ≪ 1), which is equal to 1307 nm. This value can be compared with the value of the mean wavelength of the filter, which is equal to 1310 nm. Even in calibration experiments with a noisy signal of the thermal imager, this comparison demonstrates the physical consistency of the behaviour of the interpolation equation with two orders of the extended effective wavelength. The low values of standard deviations, 0.38% for σ a 0 and 5.41% for σ a 1 , again show the small spatial dispersion of the detector's sensitivity.

Next, for each detector, the field of estimated temperatures T (u, v) is calculated by the resolution of the polynomial equation [START_REF] Sakuma | Establishing a practical temperature standard by using a narrow-band radiation thermometer with a silicon detector Temperature: its Measurement and Control in Science and Industry[END_REF]. The field of interpolation errors E T (u, v) = T -T (u, v) is calculated. Its mean spatial value, denoted ẼT , and its standard deviation, denoted σ E T , across the calibrated detectors, are calculated. Figure 7 shows the mean spatial value versus the black-body temperature T . Error bars represent the standard deviation of interpolation error with an extended effective wavelength at first and second order.

The mean spatial value of interpolation errors, ẼT , calculated with an interpolation equation with an extended effective wavelength at first order is similar to the one calculated for an extended effective wavelength at second order. For an interpolation equation with an extended effective wavelength at first order (respectively at second order), the standard deviation σ E T is between 2.7% and 12% (respectively between 2.3% and 10%).

Finally, the absolute mean value over all temperatures, denoted ĒT , of the spatial mean value ẼT is slightly lower for an extended effective wavelength at second order ( ĒT = 0.33 • C) than for an extended effective wavelength at first order ( ĒT = 0.407 • C). These values are consistent with those calculated in the previous section with simulated noisy output signal, Ǐ D , with a standard deviation of noise of σ I D = 0.5%. This value is also consistent with the temporal standard deviation of the corrected signal (see figure 6).

For a thermal imager with a standard deviation of noise around of σ I D = 0.5% equipped with a narrow bandwidth Gaussian filter, and even for a large temperature range, an interpolation equation with an extended effective wavelength at first order can be a good approximation to achieve a low interpolation error. A higher order of the extended effective wavelength moderately improves the interpolation error. These experimental results confirm the simulated results of figure 5 in which, for r ≪ 1, the difference in interpolation errors with an extended effective wavelength at first and second order is very low.

Interpolation error of the thermal imager over its entire bandwidth

The second calibration was carried out for a thermal imager over its entire bandwidth (0.9 µm to 1.7 µm), without any filter, and with the same calibration procedure as previously used.

The spatial mean and standard deviation values of the image's calibration parameters (A(u, v), a 0 (u, v), a 1 (u, v) and a 2 (u, v)) of the interpolation equation with an extended effective wavelength at first and second order are tabulated in table 3.

Compared with the case of the narrow Gaussian filter, the value of à is high. According to equation ( 4), the value of the integral over the spectral bandwidth increases with the bandwidth, as it does for ã0 and ã1 , due to the fact that the standard deviation of the spectral responsivity of the thermal imager is wider. The thermal imager has a rectangular spectral responsivity from 0.9 µm to 1.7 µm. In this case of a wide bandwidth, the value of the calibration parameter ã2 is also very high. This value justifies the choice of an interpolation equation at a higher order. It should be noted that for this calibration, the standard deviation values are also very low.

The mean spatial values, denoted ẼT , of the interpolation error field are shown versus the black body temperature T in figure 8 for an interpolation equation with an extended effective wavelength at first and second orders. The error bars represent the standard deviation, denoted σ E T .

An interpolation equation with an extended effective wavelength parametrized at first order involves a higher interpolation error than for the first order. It reaches an interpolation error up to one degree for a few temperatures. In contrast, an interpolation equation with an extended effective wavelength at second order limits the interpolation error around a mean value over all temperatures of ĒT = 0.36 • C. The second order shows a significant improvement in the value of the temperature error ( ĒT = 0.73 • C for first order versus ĒT = 0.36 • C for second order). In this case, the interpolation error depends on the order of the extended effective wavelength. This result has already been noted in figure 5 when the value of r increases.

Comparison of the simulated and measured interpolation errors

The two previous experiments confirm the trends observed in figure 5, and in this paragraph, our approach is to compare the values of simulated and measured interpolation errors. The simulation is performed with simulated noisy output signal, ǏD , which is generated with the measured relative standard deviation σ I D on the signal from the thermal imager which is around 0.5% (see figure 6). The standard deviation of the imager's spectral responsivity σ f varies and the variation of the ratio r = σ f λ 0 is from 3.24 × 10 -3 to 0.194. This simulation is intended to represent the interpolation error in a configuration similar to that of the thermal imager. Figure 9 shows the quantity ¯Ē T for simulated interpolation errors with an interpolation equation with an extended effective wavelength at first and second order. In this figure, the measured interpolation errors, ĒT , are also reported.

For an extended effective wavelength at second order, the measured interpolation errors are consistent with the simulated ones for a standard deviation of noise in the output around σ I D = 0.5%. For an extended effective wavelength at first order, the measured interpolation errors are higher than the simulated ones with the same standard deviation of noise in the output signal. The measured values should correspond to a higher standard deviation of noise in the output around σ I D = 0.6%. This result demonstrates that for middle-noise on the output signal, an extended effective wavelength at second order provides lower residual fits on the fitting parameters A, a 0 and a 1 .

As a conclusion, for our thermal imager with a middlenoise on the output signal, for a ratio r < 0.07, an extended effective wavelength at second order slightly improves the interpolation error (around 0.07 • C) with measured data, unlike the conclusions reached with simulated data in the previous section. For a ratio up to r > 0.07, as demonstrated with simulated data, only an interpolation equation with an extended effective wavelength at second order involves a significant improvement in the value of the interpolation error.

Conclusion

This paper has described the application of an interpolation equation with a temperature-dependent effective wavelength for the calibration of near infrared thermal imagers (see equation [START_REF] Sentenac | Bi-color near infrared thermoreflectometry: a method for true temperature field measurement[END_REF]). This dependence is defined as follows: the inverse extended effective wavelength is a polynomial function of the inverse temperature (see equation [START_REF] Kostkowski | Theory and Methods of Optical Pyrometry[END_REF]). The order of the function mainly depends on the mean wavelength and the standard deviation of the imager's bandwidth, the temperature range and the accuracy required.

For temperatures ranging from 300 • C to 1000 • C, a simulation performed with a free noise signal of the thermal imager showed that the interpolation error was lower than 0.1 • C if the filter's FWHM was lower than 220 nm (for a temperature-dependent effective wavelength at first order and a mean wavelength of 1310 nm). Using an effective wavelength at second order, the FWHM could be extended to 600 nm.

A second simulation with a noisy signal of the thermal imager gave the interpolation error in relation to the percentage of noise and the characteristics of the imager's bandwidth. This simulation showed that an interpolation equation with an effective wavelength at second order decreases the interpolation error with noise. However, the interpolation error was not negligible and was lower than 0.1 • C only for a filter bandwidth lower than 30 nm.

The simulations were compared with experimental results obtained by a thermal imager with a noise standard deviation around 0.5% and equipped with a narrow bandwidth filter (FWHM of 50 nm). An interpolation equation with an effective wavelength, at first or second order, provided a comparable interpolation error which was around 0.3 • C. However, for wider spectral bandwidth filters (FWHM higher than 50 nm), an effective wavelength at second order was required to keep a low value of interpolation error.

The interpolation equation with an extended effective wavelength, at first and second orders, was implemented in an original method [START_REF] Sentenac | Bi-color near infrared thermoreflectometry: a method for true temperature field measurement[END_REF] for measuring the true temperature field by fusing the radiance temperature and the reflectivity measurement fields.

Figure 1 .

 1 Figure 1. Logarithmic representation of the maximum absolute interpolation error, E T peak , versus the ratio r = σ f λ 0 for Gaussian filters with a mean wavelength of λ 0 = 1310 nm and standard deviation values from σ f = 4.2 × 10 -2 to 0.25.

Figure 2 .

 2 Figure 2. Logarithmic representation of ¯Ē T versus the relative standard deviation of the noise of the output signal σ I D for a Gaussian filter with a ratio r = 0.194 (mean wavelength of λ 0 = 1310 nm and FWHM of 600 nm).

Figure 3 .

 3 Figure 3. Mean (over the trials and temperatures) of interpolation errors for an extended effective wavelength at first order versus the ratio r of a Gaussian filter (mean wavelength of λ 0 = 1310 nm and FWHM from 10 nm to 600 nm) and for seven levels (from 0% to 1%) of the relative standard deviation of the noise of the output signal.

Figure 4 .

 4 Figure 4. As for figure 3 but with an interpolation equation with an extended effective wavelength at second order.

Figure 6 .

 6 Figure 6. Histogram over fifty samples of corrected signal of equation (11) for a black-body temperature at 500 • C.

Figure 7 .

 7 Figure 7. Values of mean ( ẼT ) and standard deviation (σ E T ) of the interpolation error field versus the temperature for a Gaussian filter with a FWHM of 50 nm and a mean wavelength of λ 0 = 1310 nm.

Figure 8 .

 8 Figure 8. As for figure 7 but for a thermal imager without filter.

Figure 9 .

 9 Figure 9. Quantities ¯Ē T or ¯Ẽ T versus the ratio r for the noise's relative standard deviation of the output signal σ I D = 0.5% and σ I D = 0.6%.

Table 1 .

 1 Maximum temperature interpolation errors, E T peak , for each extended effective wavelength's order; fitting calibration parameter values for an extended effective wavelength at second order.

	1/ λ x	0th order 1st order	2nd order

  100 × σ a 2

	1st order 2nd order 1.84 × 10 11 2.264 × 10 11 1.91 2.14	0.765 0.761	0.38 1.85	-4.48 -34.42	5.41 3.2	-1.95 × 10 4	-3.2

Table 3 .

 3 As for table 2 but for a thermal imager without filter. DL s -1 ) 100 × σ A ã0 /µm -1 100 × σ a 0 ã1 /(K µm -1 ) 100 × σ a 1 ã2 /(K 2 µm -1 ) 100 × σ a 2

	1 λx Ã/(1st order 3.78 × 10 13 0.85 2nd order 8.84 × 10 13 5.26	1.16 1.39	0.091 0.47	-217.39 -446.62	0.207 0.81	-8.211 × 10 4	-0.97

Acknowledgments

This work was carried out in the context of the R3T (Real Time True Temperature) project, supported by the French National Research Agency.