Je souhaite remercier avant tout mes encadrants qui m'ont accompagné durant ces trois années. En premier lieu, merci à Jean Dumoulin pour l'ensemble des remarques scientifiques qui ont concrétisé cette thèse. Au-delà des échanges professionnels, merci pour ta compréhension, ta patience ainsi que tes conseils avisés.

Je remercie également Laurent Mevel pour nos discussions très pertinentes qui ont beaucoup apporté dans la thèse et ce jusqu'aux derniers instants. Mes différentes visites à Rennes ont à chaque fois permis une grande progression dans mes travaux. Ainsi, l'implication des différents membres d'I4S fût très bénéfique. En particulier, je tiens à remercier Qinghua Zhang pour les échanges que nous avons eu qui ont fait mûrir ma thèse. Même si je n'ai pas toujours compris l'ensemble des remarques immédiatement, celles-ci se sont avérées pertinentes et m'ont fait de prendre du recul sur mes travaux. De la même manière je souhaite remercier Frédéric Gillot pour le temps passé à réfléchir à mes équations. Merci notamment pour ta pédagogie et pour m'avoir présenté le Krigeage qui fût finalement un point-clé de ces travaux. Je remercie l'IFSTTAR, l'Inria ainsi que la région Bretagne pour avoir financé cette thèse. Je remercie Guillaume Moreau pour avoir bien voulu présider mon jury. Des remerciements aussi aux rapporteurs Xavier Maldague et Christophe Pradère pour leur examen minutieux du manuscrit ainsi que leurs remarques qui ont permis de finaliser et d'éclaircir certains points. J'adresse également des remerciements à Alexia Gorecki qui a accepté de faire partie du jury et Laurent Ibos pour avoir assisté à la soutenance et m'avoir suivi lors du comité de suivi de thèse. Je souhaite aussi remercier Sébastien Bourguignon de l'École Centrale Nantes pour m'avoir offert l'opportunité de donner des cours mais aussi de m'avoir guidé dans ce processus. D'une manière moins académique, ces trois années de thèses furent une expérience personnelle épanouissante, notamment grâce à l'ambiance qui règne au sein du laboratoire SII et de l'équipe I4S. Les différents collègues (permanents, stagiaires, doctorants et post-docs) rencontrés ont été une source d'enrichissement. Merci ainsi à Louis-Marie Cottineau à qui je souhaite une bonne retraite et à Vincent Le Cam (aérobic) les directeurs successifs du laboratoire. Des remerciements aussi à l'ensemble des chercheurs avec qui j'ai pu intéragir, Vincent Baltazar, David Bétaille, Xavier Chapeleau et Alex Coiret. Je me dois de saluer aussi l'ensemble des collègues qui contribuent à cette bonne entente : Yveline pour son travail de l'ombre, Ivan pour l'animation des débats autour du café, Quentin pour ses goûts musicaux, Jean-Luc pour ses connaissances en champignons et ses connaissances tout court, Jean-Marc pour nos échanges sur le traitement d'images et Jean-Philippe pour ses conseils de bricoleur. Je pense aussi à ceux qui ont été de passage et qui l'ont marqué, Ludo pour ne citer personne. Une pensée à tous les autres doctorants qui vont bientôt soutenir. Dernière ligne droite Shreedhar ! Courage Antoine, tu es le prochain sur la liste, bonne chance pour tes futures fonctions. Salut à toi Arthur à qui je n'aurais jamais dû dire que je voulais faire un triathlon, quelle idée ! Merci au futur docteur David pour tous les échanges que nous avons pu avoir et que nous aurons sur Kalman et l'informatique en général (enfin quelqu'un qui comprend mes délires informatiques !). Courage pour la suite ! Merci à tous mes collègues de bureau, notamment Guillaume qui a eu le courage de se reconvertir en architecte Java. Une grande pensée pour Nico qui m'a emmené dans des drôles d'aventures (encore désolé pour la boule noire) et qui a répondu à mes nombreuses questions mathématiques. Profite bien de la campagne, je publierais une version du jeu python un jour. Enfin, merci à Martin maître-pneu investisseur, le petit dernier du bureau qui regorge d'idées innovantes. Accroche-toi, la route est longue mais le jeu en vaut la chandelle ! Merci à Emacs pour me simplifier la vie. Mes remerciements les plus chaleureux s'adressent à toute ma famille et à ma belle-famille qui n'ont cessé de me soutenir et de m'aider durant ces années, notamment pour que la maison avance.

Je termine ces pensées avec toi Mélina. Tu as été mon soutien le plus fort grâce à ta présence à mes côtés au jour le jour et tous les sacrifices que cela a représenté. Il y a une part de toi dans l'aboutissement de ces travaux. En cette fin de thèse, tu m'as offert la plus belle rencontre qui soit et il est ainsi temps de démarrer notre nouvelle vie à trois.

Conclusions et perspectives

Les travaux effectués dans cette thèse mettent en perspective l'utilisation de la thermographie infrarouge pour la surveillance de structures comme outil quantitatif. Une étude bibliographique a permis de mettre en avant trois moyens d'améliorer la précision de la mesure de la température.

Premièrement, l'exploitation de données issues de multiples capteurs permet de corriger les effets environnementaux et d'améliorer ainsi l'estimation de la température. Lorsque les données locales ne sont pas disponibles, des solutions en ligne et libres d'accès sont proposées. Ensuite, le développement d'un simulateur d'échanges radiatifs diffus sur des scènes 3D complexes permet d'optimiser de nouvelles solutions d'instrumentation et ce, en amont de toute exploitation. Ce simulateur peut être amélioré pour intégrer d'autres phénomènes physiques tels que la diffusion ou la convection, ainsi que des environnements participatifs. De plus amples modèles de caméras et détecteurs peuvent aussi être ajoutés. Enfin, un ensemble de méthodes utilisant l'inférence statistique a été proposé. Chacune de ces méthodes s'appuie sur un ensemble d'hypothèses et des méthodes de résolutions différentes. Si les méthodes MCMC et CMA-ES offrent de bons résultats, leur temps d'exécution reste long et un suivi dynamique des paramètres n'est pour l'instant pas possible. A contrario, la méthode KIKF offre une résolution en un temps relativement faible. Des travaux complémentaires, notamment sur l'optimisation des paramètres intrinsèques au filtre pourraient être conduits. De plus, d'autres approches d'optimisation pourraient être étudiée en vue d'améliorer la convergence du filtre.

Dans tous les cas, les méthodes proposées offrent un résultat prometteur. En effet, une quantification de la mesure de la température par thermographie infrarouge est réalisée grâce à une estimation de l'incertitude, inhérente à ces méthodes statistiques.

Résumé étendu en Français

Introduction L e travail réalisé dans cette thèse s'inscrit dans un contexte de développement croissant de l'utilisation de nouvelles technologies pour la surveillance de structures. Le contrôle de santé des infrastructures est primordial afin d'assurer la pérennité des réseaux de transports et des constructions. Au regard des récents évènements en Italie, la prévention et la réparation des infrastructures est un enjeu à la fois social et économique. Pour pallier la vétusté des structures, des moyens de mesures afin de quantifier leur état sont nécessaires. En particulier, l'utilisation de caméras infrarouges bas-coûts pour la surveillance long-terme d'infrastructures est prometteuse grâce aux récentes avancées technologiques du domaine. En effet, les derniers développements technologiques ainsi que la miniaturisation des détecteurs infrarouges non-refroidis permettent d'avoir des caméras bas-coûts et faciles à mettre en oeuvre. La thermographie infrarouge vise alors à fournir une mesure sans-contact et plein champs de la température. Cependant, une mesure précise de la température des surfaces observées in-situ se heurte au manque de connaissance des propriétés radiatives de la scène. En effet, la thermographie infrarouge permet de mesurer un flux radiatif arrivant à la caméra. Cependant, ce flux est la résultante de différentes contributions radiatives : l'émission propre de l'objet, permettant de remonter -en théorie -à sa température, mais aussi l'atmosphère, le soleil, les autres objets environnants de la scène, l'optique de la caméra, etc.

Une étude bibliographique de la mesure par thermographie infrarouge ainsi que l'écriture des équations du transfert radiatif met en exergue ce phénomène complexe multi-varié. De cette étude se dégage trois voies d'amélioration de la précision de la mesure par thermographie infrarouge.

Surveillance long-terme d'infrastructures par thermographie infrarouge

Dans un premier temps, il est montré que la correction des paramètres environnementaux inhérents à la mesure ainsi que l'exploitation de données issues d'une instrumentation multi-capteurs permettent d'accroître la précision de la mesure. Dans l'optique d'une mesure précise, les processus de calibration spatiale et thermique sont tout d'abord présentés. Afin d'exploiter l'ensemble des données provenant d'une instrumentation multi-capteurs, un logiciel de visualisation et de traitement est présenté. En particulier, un module YAML permet un traitement d'un large nombre de données par lots ; facilitant ainsi l'ensemble du processus. L'interface graphique du logiciel permet une vue d'ensemble des mesures et assiste les utilisateurs dans l'exploitation des données. Le modèle de conversion des im-Figure 1: Principe de la mesure par thermographie infrarouge: le flux résultant dépend de l'environnement et des différentes réflections de la scène.

ages infrarouges en niveaux numériques vers la température est parallélisé par une approche GPGPU (General-purpose processing on graphics processing units), afin d'obtenir un logiciel de traitement efficace.

Ensuite, l'exploitation des données multi-capteurs montre qu'il est possible d'améliorer l'estimation de la température. De plus, des solutions faisant usage de bases de données climatiques en ligne et libres d'accès sont proposées pour pallier tout manque d'instrument de mesure.

Enfin, une étude de sensibilité des différents paramètres est menée. Il apparaît alors que l'émissivité, propriété radiative propre à tout matériaux, est un facteur prépondérant dans l'estimation de la température ; nonobstant l'usage de données multi-capteurs. Par conséquent, des méthodes d'estimation conjointe d'émissivité et de température sont proposées dans la suite.

Simulateur multi-spectral de scènes 3D

Dans un second temps et dans l'optique d'évaluer les méthodes d'estimation conjointe d'émissivité et de température, un simulateur de scènes 3D complexes dans l'infrarouge a été étudié et développé. Suite à une étude bibliographique, la méthode des radiosités progressives est sélectionnée pour son approche par éléments finis. Cet algorithme a été implémenté en utilisant l'accélération matérielle des machines ainsi que les dernières bibliothèques graphiques. Une approche GPGPU permet une parallélisation des calculs pour un rendu accéléré des radiosités. Afin de pallier ces précédentes limitations, une dernière méthode reposant sur un filtre de Kalman en interaction krigé (KIKF) est proposé. Cette méthode combine la précédente méthode du filtre de Kalman en intéraction avec une méthode de Krigeage dans le but d'inférer un modèle de covariance spatiale à l'émissivité. L'hypothèse de température homogène locale est aussi considérée. L'ensemble des solutions possibles est alors réduit permettant de résoudre le problème en mono-spectral. L'avantage de cette méthode est sa capacité à suivre l'évolution des paramètres et de la température en même temps. Ainsi, lorsque le krigeage dégénère, il est possible de remarquer une modification des propriétés spatiales. De plus, par comparaison aux précédentes métaheuristiques, l'exécution est réalisée en un temps raisonnable qui est potentiellement parallélisable (moins d'une minute d'exécution entre deux mesures sans parallélisation). Aerosol related

• atm Atmospheric contribution

• env Environmental contribution

• g Mixed gas related (mainly O 2 and CO 2 )

• λ Spectral quantity (wavelength dependent) µm

• ∆λ i Defined on the median of the spectral band

• ∆λ i Defined on a given spectral band (

∫ λ i +∆λ/2 λ i -∆λ/2
• (i,j)
Defined for a given pixel (i, j)

• ↑ Up-welling quantity

I Radiant intensity W • sr -1 Ĩ Directional emissive power W • sr -1 J Radiosity W • m -2 L Radiance W • m -2 • sr -1 l Longitude M Radiant exitance W • m -2

Introduction

I

nfrared thermography is a non-contact and full-field measurement technique used in numerous infrastructure diagnosis applications thanks to its non-invasive nature and large scale implementation potential. In the context of energy savings, the thermal characterization of buildings envelopes plays an essential role. Determining the thermal properties of a structure helps at refining and integrating various phenomena in diagnostics tools in order to measure its thermal efficiency but also its aging. Obtaining refined data on large infrastructures makes it possible to optimize repair costs, and thus save both energy and economic expenses. The latest improvements on uncooled infrared (IR) cameras have brought new opportunities for inexpensive thermal diagnostics in the Civil Engineering's field. Although the technological maturity is ready, the methods to retrieve an accurate temperature measurement or at least with a known uncertainty are missing. In fact, the radiative flux that arrives at the camera's sensors is not directly related to the observed target temperature. Instead, this flux depends on a combination of many parameters that varies with the 3D geometry of the scene, the meteorological conditions, the self radiative properties of the objects, etc. In such context, the estimation of the temperature from those measurements is not trivial. If some parameters can be estimated from other local or online measurements, the most influential radiative quantity that will affect the temperature estimation is locally difficult to evaluate. This quantity is called the emissivity and depends on the object's own characteristics (material, roughness), but also the angle of view and the temperature.

As a consequence, in-situ infrared thermography for the thermal inspection of infrastructures is mainly used as a qualitative tool.

The work done in this thesis aims at providing methods to estimate simultaneously the emissivity and the temperature during long-term and in-situ multi-spectral infrared thermography measurements. One of the advantage of in-situ measurements is the presence of the sun as a natural heating source to excite the observed surfaces. Those various conditions coupled with long-term monitoring and multi-spectral measurements provide a set of data from which the emissivity and temperature can be extracted under some assumptions and statistical priors. Even more, a confidence interval is deducted to make the infrared thermography usable as a quantitative tool. In order to develop such method, the equations that govern the radiative heat transfer are introduced and presented in Chapter 2. The radiative phenomenon is complex and the equations are written at different observations scales to understand fully the assumptions used in the final model equation. In particular, it is shown that even if multi-spectral data brings more information to the system it is still insufficient to estimate 1 Chapter 1 -Introduction simultaneously the temperature and the emissivity.

A bibliographical review on the IR measurement process and radiative balance, the numerical simulation of 3D thermal infrared scenes and the simultaneous estimation of emissivity and temperature through multi-spectral infrared thermal measurements has been conducted. This bibliographical study, presented in Chapter 3, highlights three ways to improve the temperature estimation accuracy.

The first way is to correct the environmental effects and exploit in-situ complementary environmental measurements. In Chapter 4, a description of how coupled environmental data measurements can improve the estimation result is presented. Moreover, when local measurements from multi-sensor acquisitions are not available, it is shown that open-data available online can also be exploited.

Then, in order to optimize new instrumentation solutions for in-situ studies and test developed emissivity and temperature estimation methods, a 3D simulation software tool has been developed

and studied. This software, presented in Chapter 5 is implemented on graphic hardware to reduce the computational time. It enables the numerical simulation of the radiative diffuses exchanges of 3D complex scenes. Thanks to an atmospheric spectral solar irradiance model and sun modelization, insitu dynamic scenes can be simulated. Furthermore, the acquisition chain of the irradiance received at the camera sensor is also modeled and complete the measurement process simulation.

Developed mathematical models for the simultaneous estimation of the temperature and emissivity are presented in Chapter 6. Two methods from the literature have been implemented and three new methods are then proposed. The first one is an interacting Kalman filter which is the combination of a particle filter on a Kalman filter. While the Kalman filter tracks the temperature evolution through time, the particle filter estimates the model's parameters. However, satisfying results are obtained only for small spectral bandwidth, which is suitable for hyper-spectral applications only. The second one is based on the covariance matrix adaptation evolution strategy (CMA-ES) meta-heuristic optimization algorithm. By assuming that on a given pixel's neighborhood the temperature is homogeneous, the equation system can be solved on a single-band image, at every moment. However, this robust method induces important computational costs. To overcome those previous issues, a combination of an interacting Kalman filter with a spatial approximation is proposed. By using a priori on the spatial distribution of the emissivity through Kriging, a single-band or multi-band emissivity and temperature estimation method that tracks the parameters through time is proposed. Those methods are applied to a study case, obtained by the previous introduced simulation tool and compared.

Finally, conclusions and perspectives on in-situ infrared thermography for the monitoring of infrastructures as well as obtained results are proposed and discussed.

Chapter 2

Context and problem positioning

P

roposed methods in this thesis aim at making the temperature measurement of infrared thermal structure monitoring more accurate, or at least with a quantifiable error. This can only be achieved by identifying the surface properties of the observed object and in particular its emissivity, as exposed in this chapter. After introducing the concepts of the radiative transfer theory, the main equations to describe radiative exchanges in this thesis will be derived. Then, the main scientific challenges for in-situ infrared thermal monitoring are presented as well as the technological context on which this thesis is part of.

Structural Health Monitoring (SHM)

Structural health monitoring aims at detecting and characterizing damages on engineering structures. It has become a major public concern due to aging infrastructures and their intensive use [START_REF]L'entretien du réseau routier en France[END_REF].

A damage can be any change in the geometry or the material's properties of the infrastructure system that may affect its performance. The system is observed over time thanks to various sensors. This gathered information is then exploited to prevent catastrophic events, save maintenance costs but also for prognostic. Therefore, SHM relies on multiple pillars (see Fig. 2.1). The first one concerns the diagnosis of the infrastructure through the observed physical phenomenon, which can be seen as an extension of Non Destructive Evaluation (NDE) [START_REF] Balageas | Introduction to Structural Health Monitoring[END_REF]. The second one is the monitoring system in itself and all the parameters that defines it: sensors technology, acquisition architecture, data processing, communication, etc. Therefore, SHM needs to consider the monitoring as a whole sub-system inside the structure: embedded and connected monitoring system integration, communication layers, power and even realtime data processing ability. All those elements need to reconsider the design and management of the structure to make it "smarter". Such complete systems have been studied for the monitoring of road transport infrastructures for example in [START_REF]Integrated System of Structural Health Monitoring and Intelligent Management for a Cable-Stayed Bridge[END_REF], [START_REF] Tennyson | Structural Health Monitoring of Innovative Bridges in Canada with Fiber Optic Sensors[END_REF] or [START_REF] Crinière | Contribution au développement d'outils d'analyse de séquences d'images infrarouges: Application au contrôle non destructif de structures de Génie Civil[END_REF]. In those SHM applications, the monitoring is long-term, which needs to have a fully autonomous acquisition systems. Those systems are therefore critical for the monitoring and need to be scalable to any specific application. This implies to handle the various sensors [START_REF] Crinière | Management of Local Multi-Sensors Applied to SHM and Long-Term Infrared Monitoring: Cloud2IR Implementation[END_REF], the amount of data to store, send [START_REF] Le | Applications Des Réseaux de Capteurs Intelligents et de La Communication sans Fil à l'instrumentation Des Structures de Génie Civil[END_REF] and process [START_REF] Dumoulin | The Detection and Thermal Characterization of the Inner Structure of the 'Musmeci' Bridge Deck by Infrared Thermography Monitoring[END_REF][START_REF] Dumoulin | Infrared Thermography System for Transport Infrastructures Survey with Inline Local Atmospheric Parameter Measurements and Offline Model for Radiation Attenuation Evaluations[END_REF]. One concern about most civil engineering structures is their subject to environmental parameters due to their outdoor characteristic.

The thermal monitoring of infrastructure is an interesting component of SHM. In fact, outdoor ral heating excitation source for the system [START_REF] Lascoup | On the Feasibility of Defect Detection in Composite Material Based on Thermal Periodic Excitation[END_REF]. Thermo-physical properties of the structure can therefore been followed through time. In this thesis, emphasis is given to the monitoring of infrastructures through infrared cameras. In fact, infrared thermography cameras have experienced significant technological development over the last decade. Uncooled cameras with good resolution, relatively inexpensive and even multi-spectral are now possible [START_REF] Yon | 320 x 240 Uncooled IRFPA with Pixel Wise Thin Film Vacuum Packaging[END_REF]. Infrared thermography offers an interesting solution to monitor structure's thermal evolution through time. First, infrared cameras provide a non-contact and non-destructive technique which simplifies the setting up on field. Then, infrared thermography offers a multi-point measurement with a large field of view, suitable for large scale structures [START_REF] Clark | Application of Infrared Thermography to the Non-Destructive Testing of Concrete and Masonry Bridges[END_REF]. However, in-situ infrared measurements are often used as qualitative data only [START_REF] Grinzato | Quantitative Infrared Thermography in Buildings[END_REF].

In fact, the measured flux that arrives at the sensor depends on the surrounding and meteorological environment (see Section 2.4.1 and Chapter 4). As a consequence, complementary data measurements are needed to estimate the temperature of the observed object. The objective of this thesis is to use this multiple sensors configuration and propose a method to estimate the temperature with a known uncertainty. The sensitivity of the measurement model to its different parameters will be studied on actual measurements.

Some elements of the thermal radiative transfers theory are given in the following section. This theoretical background will be used to derive the equations that govern the main phenomena that occur during infrared thermography.

Thermal radiative transfers

Radiation transfer is one of the three means for thermal energy to be transferred with conduction and convection. The diffusion represents the transfer of heat via molecular interactions whereas radiation represents the transfer of heat via photons/electromagnetic waves. The radiative transfer theory (RT) aims at studying the absorption, emission and scattering of electromagnetic radiation as it passes through a medium. Such interactions can be described mathematically by the equation of radiative transfer. RT is involved in a wide variety of fields, including remote sensing, atmospheric science, astrophysics, optics etc. The purpose of this section is to provide the main principles of the RT theory and to introduce its concepts through the equations. Those definitions and a more detailed description of RT theory can be found in standard textbooks [START_REF] Howell | Thermal Radiation Heat Transfer. 5th[END_REF][START_REF] Chandrasekhar | Radiative Transfer[END_REF][START_REF] Vriendt | La transmission de la chaleur : introduction au rayonnement thermique[END_REF]

Radiative transfer theory

For the purposes of this study, radiation is viewed as the transport of energy in electromagnetic waves. All substances continuously emit electromagnetic radiation due to the different interactions and motions of the molecules and atoms they are made. The term thermal radiation is reserved for the visible and infrared portion of the spectrum since the related radiation should be detected by either heat or light [START_REF] Buglia | Introduction to the Theory of Atmsopheric Radiative Transfer[END_REF]. The infrared part of the electromagnetic spectrum is given in Fig. 2.2. In addition to emitting their own radiation, the radiation can be scattered, absorbed, reflected and transmitted (also, those effects can affect the frequency of the radiation beam). All those different interactions make the thermal radiation a complex phenomenon. Absorption and scattering effects can be studied on their own: the size of the particle relative to the wavelength of the radiation, the index of refraction (complex) of the material and also the shape of the particle will have an impact on the scattered proportion of the incident radiation. All of this leads to the study of different parameters such as the particle's phase function, the scatter cross section and the absorption cross section. The purpose of this study is not to go too deeply in the description of all the physical effects that occur since it will make the approach too complicated. However, it is still important to know that other physical phenomena exist at lower scales that cannot be taken into account in the presented models. For more readings on this subject, see [START_REF] Van De Hulst | Light Scattering by Small Particles[END_REF] and [START_REF] Bohren | Absorption and Scattering of Light by Small Particles[END_REF]. Two approaches can be used for describing the radiation theory: the classical electromagnetic theory and quantum mechanics. Whereas quantum mechanics will take into account the microscopic interactions of the radiation with matter, the classical theory such as Maxwell's equations will only derive the macroscopic behavior. As explained before, we cannot go too deeply onto lower scales since the mathematical complexity of the approach will increase and become unusable for practical study cases. Therefore, a combination of the two approaches is generally used in the litterature. For example, the spectral intensity is derived from the photon model whereas the absorption, transmission and emissivities are considered from a more classical approach as constant of proportionnality in the equations. This is why most of the time those parameters are derived from experimental studies even though they could be developed by the quantum mechanics through Einstein's transition probability coefficients. In such context, the main physical quantities involved into the RT theory and mainly described in [START_REF] Howell | Thermal Radiation Heat Transfer. 5th[END_REF] will be presented. Then, based on those definitions, the different radiative transfer equations will be derived.

Basic radiation properties

Let first start by giving some definitions on temperature. One can define temperature differently [START_REF] Becker | Surface Temperature and Emissivity at Various Scales: Definition, Measurement and Related Problems[END_REF]:

• Thermodynamic temperature: T is defined according to the second principle of thermodynamics, for a medium in thermal equilibrium, which can be measured directly by a thermometer. For nonisothermal bodies, the temperature T (x, y, z) is defined as the temperature of an elementary isothermal volume at the location (x, y, z). The surface temperature is then defined by the limit of T (x, y, 0) with a width approaching zero. • Radiometric temperature: This time, the temperature is defined from the radiance L λ emitted by a surface and defined later on in this chapter. The radiance L λ measured by a radiometer may be written with appropriate approximations and the temperature retrieved from those equations.

It is important to note that the surface temperature is well-defined for homogeneous materials at thermal equilibrium. However, it becomes much more complicated when considering non isothermal and heterogeneous materials. Most of the time in this thesis, the materials will be considered homogeneous and that thermal equilibrium is reached.

Radiation can be viewed as the transport of energy by discrete photons for which the relationship between its energy and frequency ν (or wavelength λ = c/ν) is given by Planck's relation:

ϵ = hν = hc λ (2.1)
where h is Planck's constant and c the speed of light in vacuum. Eq. (2.1) shows an interesting fact:

the energy of the radiation is inversely proportional to the radiation wavelength. This property can be easily observed by looking at the change of a burning flame as it goes hotter. Hotter temperatures mean higher energy which leads to a shift in frequency and therefore, in its color. In the XIX th century, the Stefan-Boltzmann law was established, deriving the maximum emission value for a given surface at a given temperature.

Definition 2.2.1 (Stefan-Boltzmann Law). We can describe the power radiated from the black-body in terms of its temperature:

M • (T ) = σT 4 (2.2)
This relation describes the power radiated from a black body surface (denoted with the superscript • ) which depends only on the temperature of the surface and is expressed in W/m 2 . All real surfaces will emit radiation at a rate smaller than M • which represents an ideal surface called black body (more details on black body and real surfaces radiation will be given in the next section). Then, at the beginning of the XX th century, Planck showed that the radiative energy emitted by a black body could be expressed over the spectrum a relation now called the Planck's law. Definition 2.2.2 (Planck's law). Planck's law describes the radiative energy emitted by a black body to its temperature (T ) over the spectrum through the Planck's function denoted L • λ (T ). An increase of temperature will lead to an increase of the emitted frequency.

L • λ (T ) = 2hc 2 λ 5 1 e hc λk B T -1 (2.3)
Two constants are often used to simplify the writing of Planck's law

C 1 = 2hc 2 and C 2 = hc k B : L • λ (T ) = C 1 e C 2 λT -1 (2.4)
By integrating Planck's Law over the spectrum, the Stefan-Boltzmann law is found: At the same time, two other relations linking the radiative energy emitted by a black body were found the Rayleigh-Jeans formula and the Wien's approximation law. Both of those formulae can be retrieved by Planck's law and may be used for particular applications as an approximation of Planck's law.

∫ ∞ 0 L • λ (T )dλ = σT 4 (2.
Definition 2.2.3 (Rayleigh-Jeans formula). By expanding the exponential of the denominator of Eq. (2.21) in a series:

e hc λkT = ∞ ∑ n=0 ( hc λkT ) n n! , ∀ hc λkT > 0 (always the case) (2.6)
For λT >> hc k , we can approximate the denominator by using the development until n = 1 which gives:

L • λ (T ) = 2ckT λ 4 (2.7)
Definition 2.2.4 (Wien's approximation law). This was the model originaly proposed by Wien to describe the complete spectrum of thermal radiation until Planck's one.

Wien's approximation can now be derived from Planck's law when hν ≫ kT :

L • λ (T ) = 2hc 2 λ 5 1 e hc λk B T -1 ≈ 2hc 2 λ 5 e -hc λk B T (2.8)
The temperature of a black body can therefore be deducted from Planck's law: In the next section, the main quantities of the radiative transfer theory will be defined in order to derive -at the end -the radiative transfer equations and get a better understanding of this physical phenomenon.

T = C 2 ln ( (C 1 +λ 5 M • λ ) λ (λ 5 M • λ ) λ ) (2.9)
Wien's approximation is often used to approximate Planck's law and particularly to linearize the equation system as shown on the emissivity -temperature section. This approximation can be justified in common civil engineering applications: since hc k = 14388µm • K, in the LWIR band (λ ∈ [8µm; 14µm]), we must have 959.2 >> T .

Note

Radiative transfer quantities definitions

Radiative transfer represents exchanges between surfaces and inside medium. Due to the complexity of the phenomenon it is important to keep in mind that the radiative transfer deals with waves that propagate among surfaces and through a medium. Therefore, the surfaces particularities and the geometry between the surfaces will have an effect on this heat transfer. To describe those exchanges, it is necessary to remind the solid angle definition. Definition 2.2.5 (Solid angle). The solid angle is an equivalent of the planar angle but in three-dimensional space. It represents a part of 3D space delimited by a cone (the vertex of the cone is the vertex of the solid angle). It is a measure of the field of view from a particle point that an object covers and therefore used in radiometry to determine the exchanges between two bodies. The solid angle is defined as the ratio of its base area to the square of chord length. Its dimension is the steradian sr:

dΩ = A r 2 (2.10)
The solid angle can also be defined on an unit hemisphere, which will be used later on (see Fig. 2.5):

dΩ = rdθr sin(θ)dΦ r 2 = sin θdθdΦ (2.11)
In radiometry, the quantities are defined for a given solid angle, to represent the beam direction.

However, for approximation and simplicity needs or sometimes due to uniformity, the quantities may be referred to as hemispherical. In that case, it means that the quantity is integrated on the unit hemisphere over the surface. Hemispherical properties get rid of the angle component and propose a more global representation of the quantity, as an average over the surface.

Definition 2.2.6 (Radiance). The analysis of radiation field deals with the analysis of the amount of energy dQ λ (x, Ω, t) that is transported across an element of area dA [START_REF] Chandrasekhar | Radiative Transfer[END_REF]. The radiation will then be highly dependent on the spatial geometry and the location of this area element. Let D ⊂ R d (d ∈ {2, 3}) and A be the unit sphere of R d . This energy is defined by unit of projected surface, dA cos(θ)

dQ λ (x, Ω, t) = L λ (x, Ω, t) cos θdλdAdΩdt (2.12)
where θ is the angle which the direction considered makes the outward normal to d⃗ n. The quantity

L λ (x, Ω, t) is called spectral radiance (W • sr -1 • m -2 • µm -1 )
We have to note the fact that L λ (x, Ω, t) depends on seven variables in total:

• 3 variables giving the position in space x,

• 2 variables giving the direction of the radiation Ω, in which we use the inclination angle θ and the azimuth angle φ, • The time variable t,

• The wavelength λ or sometimes the energy hν or the frequency ν = c λ .

This dimensionality makes the solving of the general radiative transfer equation difficult both analytically and numerically. Moreover, most of the time, the radiation has an impact on the crossed material, which also has an impact on the incomming radiation. Therefore, the resolution of the problem must be coupled, iteratively.

For simplification in the writting, and regarding the amount of dependent variables that L λ (x, Ω, t)

depends on, we may abusively write 1 : where Q is the radiant energy emitted, reflected, transmitted or received (J). From Eq. (2.1) we can define the spectral flux in frequency ϕ ν = ∂ϕ ∂ν and the spectral flux in wavelength ϕ λ = ∂ϕ ∂λ 1 .

L λ ≡ L λ (x,
Finally, one can combine Eq. (2.4) and Eq. (2.1) to get the spectral radiance in terms of flux which represents how much of the power is emitted, reflected, transmitted or received by a surface when an optical system is looking at the surface from a specified angle of view as unit.

L λ = ∂ 2 ϕ λ cos(θ)∂Ω∂A (2.15)
where:

• ϕ λ is the radiant spectral flux (emitted, reflected, transmitted or received),

• Ω is the solid angle,

• ∂A cos(θ) is the projected surface. Therefore:

J = ∂ϕ ∂A = M + J reflected + J transmitted (2.18)
When the surface is opaque, J transmitted = 0 =⇒ J = M + J reflected . 1 The same notations are applied to the other quantities.

Definition 2.2.11 (Radiant intensity). The radiant intensity is the ratio of the flux ϕ e emitted per unit solid angle dΩ:

I = ∂ϕ e ∂Ω (2.19)
As explained in [START_REF] Howell | Thermal Radiation Heat Transfer. 5th[END_REF], the intensity of radiation from a black-body is defined on the basis of the normal area. It means that I is independent of the direction of emission, by definition which enables the description of a black-body intensity without defining the normal of the surface or the angles (θ, Φ).

Definition 2.2.12 (Directional emissive power). The radiant intensity is defined by its normal which means that the maximum energy emitted by a surface is reached if the surface is normal to the receiver.

When there is an angle between the two surfaces, then the energy reaching the receiver is not the radiant intensity anymore. By assuming that the emission is uniform in all directions, the directional spectral emissive power for a black surface can be derived:

Ĩ(θ, Φ) = I cos(θ) = Ĩ(θ) (2.20)
where Ĩ stands for the directional emissive power. This equations is valid for Lambertian surfaces which represent a perfect diffusion: the brightness of the surface appears the same no matter the observer's angle of view (Fig. Based on this law, one can integrate the emission of a black-body through a unit hemisphere over the surface. The spectral emission from an infinitesimal area dA per unit of time and unit surface area passing through the element on the hemisphere is given by: Ĩλ cos(θ)dΩdλ = Ĩλ cos θ sin θdθdΦdλ (2.21)

By integrating over the hemisphere and using Eq. (2.20)

Ĩλ dλ = I λ dλ ∫ 2π Φ=0 ∫ π 2 θ=0 cos θ sin θdθdΦ = πI λ dλ (2.22)
The black-body hemispherical emissive power is therefore π times the black-body intensity.

We also have the equivalent in luminance. For a flux leaving the surface in an isotropic way:

M λ = πL λ (2.23)
Definition 2.2.13 (Bouguer's law). Finally, Bouguer's law defines the relation between the irradiance E of a receiving surface, due to a source S and the intensity of that source in the direction of the receiver lying at a distance d. In the case of a nonparticipating media.

E = I cos(θ) d 2 (2.24)
The Lambert's cosine law can be recognized in this expression as well as the inverse square law.

Now that the main quantities that describe the radiation theory have been introduced, the next section aims at developing the principles behind the actual phenomena that occur at a given surface.

Surface properties

Until now, the black body has been said to be a body that emits the maximum of energy for a given temperature. To be more accurate, a black body is an idealized physical body that absorbs all incident electromagnetic radiation, regardless of the frequency or the angle of incidence [START_REF] Planck | The Theory of Heat Radiation[END_REF]. At the opposite of a black body, a white body will reflect all incident rays uniformly and entirely in all directions.

Whereas black body represents an idealized object, real materials are considered to emit energy at a fraction of a black body energy levels. The ratio that links the two materials is called emissivity (ϵ, more rigorously defined later on). Finally, surfaces that have their properties independent of wavelength are called grey surfaces. When dealing with real radiating materials, multiple interactions at interfaces may occur and are schematized in Fig. 2.7. The emission has been defined previously for a black body, as well as the incident radiation. However, the reflectivity, absorption and transmission have not been mentioned yet and are addressed in this section.

To understand well the complexity of those phenomena, it is necessary to view those quantities as local properties of the material. As for the previous quantities, they will depend on the solid angle, the wavelength and the temperature but also on the intrinsic characteristics of the surface's material.

Wood and metal will not absorb or reflect the same way; as well as polished iron will not have the same properties as unpolished one. From a practical point of view, giving values to those quantities is difficult. In fact, when the reflection or the absorption of a given material is given in the literature (usually the hemispherical one), it cannot take into account the local roughness of the material or its detailed atomic composition. However, some analytic results exist for specific materials, in particular for dielectric ones. Determining the material's properties represents a field on its own. for particular materials such as dielectric ones, derived from the electromagnetic theory, analytical formulations are most of the time unknowns for other materials. The reflectivity can be defined as a bidirectional function and denoted by ρ λ (θ 1 , ϕ 1 , θ 2 , ϕ 2 ). In that case, it means that the reflected beam depends on the incoming beam solid angle.

The directional hemispherical spectral reflectivity can be derived from this definition:

ρ λ (θ, ϕ) = d 2 Q λ , r(θ, ϕ) d 2 Q λ , i(θ, ϕ) = ∫ ρ λ (θ r , ϕ r , θ, ϕ) cos(θ r )dΩ (2.25)
Similarly, the transmissivity can be defined as a bidirectional spectra quantity τ λ (θ 1 , ϕ 1 , θ 2 , ϕ 2 ).

For an opaque material, τ λ (θ 1 , ϕ 1 , θ 2 , ϕ 2 ) = 0.

Definition 2.2.15 (Emissivity). The emissivity represents the effectiveness of a surface of a material in emitting energy as thermal radiation. Based on the Planck's law and as explained previously, the directional spectral emissivity ϵ λ (θ, ϕ, T ) can be defined mathematically as the ratio of the radiance L λ (T ) emitted by a body at temperature T to the radiance L • λ (T ) emitted by a black body at the same temperature T :

ϵ λ (θ, ϕ, T ) = L λ (T ) L • λ (T ) (2.26)
Definition 2.2.16 (Absorptivity). The absorptivity represents the part of the incoming energy that is absorbed by the material. It is denoted as α λ (θ, ϕ, T ) and is strongly related to absorbing surface's temperature property. If the incident surrounding of the surface is a black body at temperature T • and the total absorbed energy is Q λ,a (θ, ϕ, T ), the absorptivity can be written as:

α λ (θ, ϕ, T ) = d 2 Q λ,a (θ, ϕ, T )dλ I • λ (T • )dA cos(θ)dΩdλ (2.27)
Definition 2.2.17 (Kirchhoff's law). Kirchhoff's law links the absorptivity and the emissivity capacities of a given body. It actually depends on which of the spectral, directional or hemispherical quantities are considered. A general statement is that at thermal equilibrium, for a given body in a radiative enclosure, the body needs to return the same amount of energy its has absorbed to the enclosure in order to keep the thermal equilibrium.

By using the definition of the emissivity, the energy emitted per unit time by an element dA in a wavelength interval dλ and solid angle dΩ is

d 2 Q λ,e (θ, ϕ, T )dλ = ϵ λ (θ, ϕ, T )I • λ (T )dA cos(θ)dΩdλ (2.28)
To ensure thermal equilibrium and isotropy within the black enclosure:

d 2 Q λ,e (θ, ϕ, T )dλ = d 2 Q λ,a (θ, ϕ, T )dλ (2.29)
which leads by using the definition of emissivity and absorptivity to:

α λ (θ, ϕ, T ) = ϵ λ (θ, ϕ, T ) (2.30)
Please note that when considering the polarization of the radiation, this equality holds only for the components that have the same polarization.

For an opaque material, we can write the directional spectral emissivity and

reflectivitiy ϵ λ (θ, ϕ, T ) = 1 -ρ λ (θ, ϕ, T ).

Note

The emissivity is one of the most important parameters in the physical process of radiation heat transfer. The emissivity depends on the view angle, wavelength, temperature and intrinsic characteristics of the material. ). The alumina is from previous studies [START_REF] Monchau | Mesure d'émissivité pour la thermographie infrarouge appliquée au diagnostic quantitatif des structures[END_REF].

Most of the time, the emissivity of a material is given on a particular spectral band and integrated over the directions (called hemispherical emissivity). The material is then considered optically thick.

As a consequence, emissivity values found in the literature must be taken with care. In fact, the variation of the emissivity of a surface with the observation's angle means that surfaces that are not flat will have a locally varying apparent emissivity, even for a given material. The roughness of a material will therefore impact the emissivity value, leading to important uncertainties (see Fig. 2.9).

Radiative Transfer Equation (RTE)

Differential form

The RTE represents what affects a beam radiation as it travels. Three main effects can occur:

• Absorption ⇒ Energy is lost • Emission ⇒ Energy is gained • Scattering ⇒ Energy is redistributed
We will see in that section how to consider those different effects and how it affects the final equation.

We can construct the radiative transfer equation by summing the different effects that will affect This physical phenomenon can be formulated mathematically through the differential form of the radiative transfer equation:

1 c ∂L ν ∂t + Ω.∇L ν = j ν -(κ ν,s + κ ν,a )L ν + 1 4π ∫ Ω κ ν,s L ν dΩ (2.32)
where:

• c is the speed of light,

• L ν is the spectral radiance,

• Ω is the propagation direction,

• k ν,a is the absorption coefficient,

• k ν,s is the scattering coefficient,

• j ν is the emission coefficient,

• κν 4π ∫
Ω L ν dΩ is the radiation scattered from other directions onto the surface. Most of the time the term 1 c ∂Lν ∂t is neglected since it would involve fast responsive systems compared to the involved distances. 2 More details on this equation can be found in [START_REF] Howell | Thermal Radiation Heat Transfer. 5th[END_REF]. Analytical solutions of the RTE exist only under some assumptions and numerical ones based on finite elements 2 Fast here means systems for which the variation of Lν in time is comparable to the light of speed. For most of the systems, we have 1 c ≪ ∂Lν ∂t which justifies the use of this approximation. However, when the distances become bigger like in astronomy for example, this approximation should be re-considered. models are more often used. In the thermography field, this equation is replaced by a simplified radiometric equation, often partially integrated in space that represents the radiative exchanges at a more macroscopic level, presented in section Section 2.4.1.

Radiative transfer in nonparticipating medium

The previous equation is suited for participating mediums. However, obtaining a solution for such equation is non trivial especially in 3-dimensions. This equation is often simplified and used for solving complex systems through finite elements methods by using physical approximations. Eq. (2.32) describes a change in radiative energy during the propagation of a beam of radiative energy in a participating medium. However, for a numerous number of systems, the scattering from a participating medium can be neglected without introducing any significant error to calculations. When it is the case and that there is no scattering particles involved, the RTE is reduced to a partial differential equation.

Exchange between two surfaces

Let now consider a simpler case of nonparticipating (i.e. no scattering, no absorption and no emission) medium. Instead of looking at a single beam propagation phenomenon, let make the balance on a given surface since the medium is nonparticipating. Let consider two infinitesimal surfaces of area dA 1 and dA 2 and separated by a distance r. In this simplified case, the infinitesimal surfaces emit uniformly in all directions and reflect the energy received from other surfaces. Properties of emission and reflections will be more described later on. The energy dQ λ,d1 leaving dA 1 in direction x 1 is given by: dQ λ,d1 = Ĩλ,d1 cos θ 1 dλdA 1 dΩ 1

(2.33)

The solid angle can be written by:

dΩ 1 = cos(θ 2 ) r 2 dA 2 (2.34)
where dA 2 cos(θ 2 ) is the projection of dA 2 on the direction x 1 . When no absorption is considered in the medium, the radiative energy reaching dA 2 is therefore given by:

dQ λ,d1→d2 = Ĩλ,d1 ( cos(θ 1 ) cos(θ 2 ) r 2 ) dλdA 1 dA 2 (2.35)
When Eq. (2.32) was integrated on the direction, Eq. (2.34) is now integrated on the area only. For finite element surface, the integration yields to:

Q λ,1→2 = ∫ ∫ A 1 ,A 2 dQ λ,d1→d2 = ∫ ∫ A 1 ,A 2 Ĩλ,d1 cos(θ 1 ) cos(θ 2 ) r 2 dλdA 1 dA 2 (2.36)
By considering the previous stated assumptions that the radiation leaving A 1 is directionally uniform, Ĩλ,d1 = Ĩλ,1 = I λ,1 π and can be taken outside the integral:

Q λ,1→2 = Ĩλ,1 ∫ ∫ A 1 ,A 2 cos(θ 1 ) cos(θ 2 ) r 2 dλdA 1 dA 2 = I λ,1 (∫ ∫ A 1 ,A 2 cos(θ 1 ) cos(θ 2 ) πr 2 dA 1 dA 2 ) dλ (2.37)
The term inside the parenthesis is called a form factor (sometimes configuration or view factor) and represents the geometrical exchanges interaction coefficient between the two surfaces A 1 and A 2 . It depends only on the geometry of the scene and is denoted F 1→2 :

Q λ,1→2 = I λ,1 A 1 F 1→2 dλ (2.38)
where

A 1 F 1→2 = ∫ ∫ A 1 ,A 2 cos(θ 1 ) cos(θ 2 ) πr 2 dA 1 dA 2 (2.39)

Radiosity equation

Let now consider an enclosure S ∈ R 3 with a nonparticipating medium inside and opaque surface.

Let x be a parametric point of such enclosure and an infinitesimal surface dA (see Fig. 2.10). Based on its definition, the radiance at a given point x of the enclosure and given time can be derived as:

L λ,leaving (x, Ω) = L λ,emitted (x, Ω) + L λ,reflected (x, Ω) (2.40)
The reflected term L λ,reflected (x) is the radiation part that is received at x (by all the other radiated parts) and reflected by x. Based on the definition of the bidirectional reflectivity and Eq. (2.37), this term can be written as a sum of the energy received from all other directions.

L λ,reflected (x, Ω) = ∫ (x ′ ,Ω ′ )∈S,cos(θ)>0 ρ λ (x, Ω, Ω ′ )V (x, x ′ )L λ,leaving (x ′ , Ω ′ ) cos(θ)dΩ ′ (2.41)
where cos(θ) = n x .(x-x ′ )/||x-x ′ || 2 (Lambert's cosine law) and V (x, x ′ ) represents a visibility factor which is equal to 1 if x sees x ′ and 0 otherwise, as illustrated in Fig. 2.10. By using the expression of dΩ ′ as previously, the integral can be integrated with the areas and expressed through the form factors: Finally, the radiance at x and solid angle Ω is:

L λ,reflected (x, Ω) = ∫ (x ′ ,Ω ′ )∈S,cos(θ)>0 ρ λ (x, Ω, Ω ′ )V (x, x ′ )L λ,leaving (x ′ , Ω ′ ) cos(θ) cos(θ ′ ) ||x -x ′ || 2
L λ,leaving (x, Ω) = L λ,emitted (x, Ω) + ∫ (x ′ ,Ω ′ )∈S,cos(θ)>0 ρ λ (x, Ω, Ω ′ )F (x, x ′ )L λ,leaving (x ′ , Ω ′ )dA ′ (2.43)
where

F (x, x ′ ) = V (x, x ′ ) cos(θ) cos(θ ′ ) ||x -x ′ || 2 2 (2.44)
Let suppose now that the surfaces are lambertians, and therefore does not depend on the direction.

The Eq. (2.24) can be used for diffuse surfaces for which the reflectivity is independent of (θ r , ϕ r ):

ρ λ,h = ρ λ (θ, ϕ) ∫ cos(θ r )dΩ = πρ λ (θ, ϕ) (2.45)
which is analogous to the Eq. (2.22). As pointed by [START_REF] Howell | Thermal Radiation Heat Transfer. 5th[END_REF], since ρ λ,h = 1α λ , the reflectivity cannot depend on the angle of incidence:

ρ λ,h = πρ λ (2.46)
Therefore, for diffuse surfaces:

L λ,leaving (x) = L λ,emitted (x) + ρ λ (x) π ∫ S F (x, x ′ )L λ,leaving (x ′ )dA ′ (2.47)
Instead of dealing with punctual quantities, it is possible to integrate those quantities over the hemisphere:

∫ ∩ L λ,leaving (x) cos(θ)dΩ J λ (x) = ∫ ∩ L λ,emitted (x) cos(θ)dΩ M λ (x) + ρ λ (x) π ∫ S F (x, x ′ ) ∫ ′ ∩ L λ,leaving (x ′ ) cos ( θ ′ ) dΩ ′ ∫ S F (x,x ′ )J λ (x ′ )dA ′ dA ′
(2.48)

J λ (x) = M λ (x) + ρ λ (x) π ∫ S F (x, x ′ )J λ (x ′ )dA ′ (2.

Infrared thermography

Simplified radiometric equation

The main quantities involved in the radiation theory have been presented and two radiative transfer equations have been derived. In this part, a simplified model that governs the radiative exchanges for in-situ infrared thermography measurements will be exposed. This equation will rely on looking at the exchanges as fluxes on the scene instead of beams and focuses on the irradiance received at the camera's detectors. Thermography aims at determining the spatial distribution of heat in the observed scene as well as its time dependence. However, the actual radiation received at the sensors will not only depend on the observed object thermal properties but it will also depend on environmental conditions (e.g.

weather), the thermo-optics and the geometrical properties of the measurement scene [START_REF] Chang | An Investigation Into The Factors Affecting Infrared Temperature Measurements For Building Applications[END_REF][START_REF] Chrzanowski | Influence of Measurement Conditions and System Parameters on Accuracy of Remote Temperature Measurement with Dualspectral IR Systems[END_REF][START_REF] Van | The Influence of Environmental Parameters on the Thermographic Analysis of the Building Envelope[END_REF].

Therefore, the irradiance received at one sensor is the result of many radiometric interactions, different from scene to scene. Since most of the quantities derived in the previous sections are unknown when using an infrared camera, a simpler model is needed to be able to describe the quantity received at the sensor. The main interactions involved in such radiometric balance are represented in Fig. 2.11. To take into account the flux arriving at an infrared detector, a mathematical model is necessary and can be derived from a flux balance that reaches the detector [START_REF] Minkina | Simulation Analysis of Uncertainty of Infrared Camera Measurement and Processing Path[END_REF]. Therefore, we can derive the simplified radiometric equation at pixel (i, j) for a given time [START_REF] Gaussorgues | Infrared Thermography[END_REF]. In this equation, multiple approximations are implicitly done. First, the surfaces are considered as greybody in the considered band. Usually, surfaces are assumed to be diffuse (Lambertian hypothesis). The atmospheric, environment and sun contributions are assumed to be black-body equivalent radiative surfaces. The environmental contribution is considered to be an equivalent contribution of the sum of all the radiation that come from the surroundings. The sun contribution represents the radiance that arrives at the sensor, from the sun emittance that crosses the atmosphere both direct and indirect. The optical, atmospheric, sun and environmental contributions are supposed to be constant over the image as a first approximation. The contribution from the scene can be expressed as:

L (i,j) ∆λ,total (T ) =τ opt [ τ atm ( ϵ (i,j) ∆λ,obj L •(i,j) ∆λ,obj (T (i,j) obj ) + (1 -ϵ (i,j) ∆λ,obj )(L • ∆λ,env + L • ∆λ,sun ) +(1 -τ atm )L • ∆λ,atm ] + (1 -τ opt )L • ∆λ,opt (2.51)
Where τ atm and τ opt are respectively the atmospheric and optic transmissions.

From Eq. (2.50) we can find the radiative contribution of the observed object:

L •(i,j)
∆λ,obj (T

(i,j) obj ) = L (i,j) ∆λ,total (T ) -(1 -ϵ (i,j) ∆λ,obj )τ atm τ opt (L • ∆λ,env + L • ∆λ,sun ) ϵ (i,j) ∆λ,obj τ atm τ opt -(1 -τ atm )τ opt L • ∆λ,atm -(1 -τ opt )L • ∆λ,opt ϵ (i,j) ∆λ,obj τ atm τ opt (2.52)
By using the contribution of the observed object and Planck's law, it is possible to find back the object's temperature. However, this can be achieved if and only if the other parameters of the model are well known. This is the topic of the following section.

Needs for emissivity-temperature simultaneous estimation methods

Once the contribution of our object has been derived, we still need to convert this quantity to temperature. We usually use a conversion model based on a thermal camera calibration to achieve this. Such step will be detailed in Chapter 4. Nevertheless, this equation shows how difficult it is to estimate the temperature from IR measurements without prior knowledge of the observed scene. In fact, even if the involved quantities can be determined through local weather reports and phenomenological laws with some errors (see Chapter 4), the system still needs to get the knowledge of the observed emissivity to retrieve the temperature from the measurements. Assuming that all the parameters are known except for the temperature and emissivity, the Fig. 2.12 shows the solution space for a given measurement value. Needs for methods that would estimate simultaneously the emissivity and the temperature can be emphasized from this figure. Without any knowledge about the emissivity, the temperature cannot be derived. Even more, an error made on the emissivity will lead to an error made on the temperature estimation.

Analytic values of spectral emissivities are often unknowns in practical applications. For an opaque material dielectric or metallic, typical directional and spectral emissivity can be derived from the Hagen-Rubens relation and based on the electro-magnetic theory. Another possibility is to obtain spectral emissivity values from databases. However, such approach is generally not enough for in-situ thermal monitoring of infrastructures to get an accurate estimate because of the lack of information about the observed material, the unknown roughness of the observed surface or even the angle of view with To overcome those issues, temperature-emissivity separation methods appeared in order to get rid of the emissivity in the thermal measurement procedure. Finally, with the development of Bayesian inference, new separation methods appeared in order to estimate conjointly the emissivity and the temperature from statistical a priori, having the potential to give an estimate with a known standard deviation error. In this thesis, different methods that separates emissivity and temperature that are used in the pyrometry or remote sensing fields are presented and compared to new methods involving spatial and temporal information usage.

The evaluation of emissivity and temperature can be performed through active methods, on which a surface is heated by a laser for example in the case of noncontact measurements. Those methods are often used for materials characterization, defects detection and are well-developed in those fields.

However, active measurement for long-term infrared monitoring of infrastructure may not be possible for a couple of reasons. First of all, when dealing with large infrastructures, it may not be possible to perform a characterization on the entire structure. Such investigation will have an impact on the cost of the operation and may not be feasible for security purposes again. Furthermore, retrieving local information may not be useful if the structure is observed in its entirety through a particular point of view. In this thesis, the focus is done on noncontact and passive temperature / emissivity retrieval methods for in-situ infrared thermal monitoring. This topic will have a dedicated bibliographical study in Chapter 3 and proposed method will be developed in Chapter 6

Measurements bias due to environmental and spatial conditions

The wide range of variables involved in the balance have to be identified and considered in any in-situ temperature estimation through IR thermography. We will review here the different impacts of the environment, the geometry of the scene and the digitization on our temperature estimate.

Spatial resolution and angle of view

The measurement set up will influence the obtained image. As we have seen previously, the emissivity depends on the material and its roughness but it also depends on the angle of view. As a consequence, the positioning and the angle of view of the camera on the scene will have an impact on the measurement (see Fig. 2.13b). Moreover, the perspective projection induced by the camera field of view will result on a non-constant spatial discretization, that will have to be considered in order to estimate the final temperature estimation error.

Perspective projection

The image acquired by the IR camera will result in a projection of the 3D world coordinates to the 2D frame of the image. Therefore, the resulting image and projection will depend on the position of the camera onto the scene, the focal length, the aperture, the field of view etc. As said previously, the radiation measurement depends on the distance and angle of the sensor to the object. Most of the time, for in-situ civil engineering applications, large field of view are used to get a full view of the structure and the complex geometry of those 3D structures leads to a wide variety of angles between the object and the camera view (see Fig. 2.13a).

Induced discretization

The position of the camera and its optical properties has to be considered in order to achieve a proper estimation of the temperature of the observed object. However, the angle of view is not the only geometrical quantity that will impact the measurement. The discretization and in particular the non-constant spatial sampling size will give a maximum accuracy reachable that will variate within the image, as represented in Fig. 2.14. In Chapter 4, a spatial calibration procedure will be presented and the interpolation due to such resectioning will be discussed.

Mixed pixel effect

Last but not least, the discretization will lead to the so called "mixed pixel effect". In fact, one sensor may see multiple materials for which the measurement will lead to a mixture of those materials. Fig. 2.15 illustrates this effect. In a mixed pixel, the material property of the observed area can be considered as a combination of the properties of the mixed materials whereas in a pure pixel, we can consider only the material properties of the observed area. Therefore, most of the time, it is difficult to identify the emissivity of a particular pixel / material (especially when dealing with angles).

As a consequence, the radiative properties are often estimated in terms of pixels i.e. an equivalent temperature and emissivity of the real observed surface.

Infrared cameras for thermal radiative measurements

As seen in the previous sections, the estimation of temperature from infrared measurements is a challenging task in which many factors are involved, from theoretical approximation in the radiative model to the actual measurements system and configuration. In order to understand the concept behind infrared thermography, IR cameras technologies principles are briefly introduced. Then a simple use-case measurement equation in single-band and multi-band thermography are derived.

History and basic principles

Since the discovery of infrared radiation by Herschel first experiment in 1800 with a prism and measuring the temperature of each color, needs for quantitative detection of radiation for many different application field, and in particular for Defense purposes appeared [START_REF] Rogalski | History of Infrared Detectors[END_REF]. Over the past 80 years, many developments have been done in the infrared detectors which are used now in many fields, ranging from satellite defense applications to building insulation inspection. In particular, the latest improvements have lead to the development of small infrared uncooled cameras with the size of a smartphone camera: inexpensive and easy to deploy on site.

Because of the two bands of good transmission in the infrared range of the atmosphere, most detectors are divided into short-wave (SW -1.4-3µm) and long-wave (LW -8-15µm) camera devices.

However, we can still find detectors that will work in the near IR (NIR -0.75-1.4µm) or even the far IR (FIR -15-1000µm)

Cameras can then be classified into two different groups: cooled or uncooled systems. Cooled systems have to maintain a low temperature at ranges that may vary between 4K to 203.15K depending on type and performance levels. At the opposite, uncooled systems, appeared more recently, operate at ambient temperature. Such systems are cheaper but measurements may not provide the same quality.

We aimed at using our developed model on such devices, for their simplicity of exploitation and price.

Uncooled cameras Uncooled cameras come with thermal detectors (microbolometer) that when heated through IR radiation, measure a temperature change. When an infrared radiation with wavelength between 8-15µm strikes the IR absorbing metal, it heats it, thus changing its electrical resistance. The change of the resistance is measured and processed through a CMOS (Complementary metaloxidesemiconductor) ROIC (readout integrated circuit) [START_REF] Wang | Fabricating Microbolometer Array on Unplanar Readout Integrated Circuit[END_REF]. Additionally, a reflector can be put beneath the absorbing metal to send back to the absorption layer the energy that it did not absorbed previously, thus enhancing the maximum energy absorption. Since the system is uncooled, the absorbing material must be thermally isolated from the bottom. This is achieved by performing a sacrificial layer.

This micro-bridge structure makes an optical resonant cavity possible at λ/4 (i.e. 2 µm for the 8 µm to 14 µm spectral band) between the probe and the substrate. An example of a microbolometer structure is shown on Fig. 2.17. The drift of such system may be important and the calibration is an essential part of the trust you can accord to the temperature measurement [START_REF] Machin | High-Quality Blackbody Sources for Infrared Thermometry and Thermography between -40 and 1000řC[END_REF].

Main advantages

• Low power consumption

• Lower cost • Moderate sensitivity • No mechanical cooling

Main drawbacks

• Fixed and long-time constants

• Relatively slow frame rate

• No control over timing or triggering

• Sensitivity to ambient temperature =⇒ needs for thermal calibration

Cooled cameras

In order to count the number of photons, such cameras use quantum detectors that need to be cooled for cutting thermal noise. The materials are mainly semiconductors with narrow band gaps. Incoming photons will bring energy into those semiconductors causing excitation that will be measured. The necessity of cooling those detectors makes the manufacturing more complicated but the sensitivity of those cameras is much better than the uncooled systems.

Main advantages

• Very short integration times • Fast frame rates

• Superior sensitivity

• Ability to see and measure small targets

Main drawbacks

• High cost

• Power consumption

• Requires cryogenic cooling

Thermograms creation

In order to create an image from an array of detectors, two main technologies exist. The first one is to use a single detector or a ruler of detectors and perform a scan of the observed area. In a single detector camera, the image is built point by point, in a ruler of detectors, it is constructed line by line (see Fig. 2.18). On a single point detector, one can note that all points of a thermogram have identical parameters because the temperature measured at every point is measured by the same detector. Such property can be useful in detecting the temperature difference at two points of a homogeneous object. Historically, the point scanner cameras was the first system used. Then, rulers of detectors and finally arrays of detectors were created. In the last systems, there is no mechanical scanning parts: the array observes the scene through the camera optics. Such system enables the creation of cameras that can record fast thermal processes such as 100 thermograms per second. Each pixel is read 25 times per second (PAL system -Europe) or 30 times (NTSC system -USA) [START_REF] Minkina | Infrared Thermography: Errors and Uncertainties[END_REF].

Single-band pyrometry

In order to emphasis again the needs for emissivity and temperature separation, let consider a simple case in which the irradiance at a sensor located at a distance r of an object is derived. Such case of study needs some assumptions:

Assumptions

• Flat surfaces are small enough to be considered punctuals.

• Only the diffuse interaction between the sensor and the object is considered.

• The scene geometry is known Under those assumptions, the irradiance received at the sensor can be derived from Bouguer's law:

E sensor = I object cos(θ sensor ) r 2 (2.53)
Where E sensor is the irradiance at the sensor, I object , the radiant intensity leaving the object, r and θ sensor are defined in Fig. 2.19. With the infinitesimal surfaces hypothesis, one can get the luminance from the intensity:

dI object = L object cos ( θ object ) dS object (2.54)
By using the fact that the surfaces are punctuals, the irradiance received at the sensor is finally:

E sensor = cos ( θ object ) cos(θ sensor )S object r 2 L object (2.55)
By using Planck's law and defining g :

(θ sensor , θ object , r) → cos(θ object) cos(θsensor)S object r 2
, assumed to be known in this case:

E sensor = g(θ sensor , θ object , r)ϵ object C 1 λ 5 (e C 2
λT object -1)

(2.56) with C 1 = 2hc 2 λ 5 and C 2 = hc k .
By using Wien's approximation and taking the logarithm of the previous expresion, the emissivity and temperature can be separated to derive the following expression:

ln(E sensor ) + ln ( λ 5 C 1 g(θ sensor , θ object , r) ) = ln(ϵ object ) - C 2 λT object (2.57)
The first term in Eq. (2.56) is the physical quantity that will be measured at the sensor. The second term depends on the geometry of the scene and is assumed to be known here. Finally, the last two terms represent the unknowns of the system: the temperature and the emissivity of the object. Therefore, one can see that a passive infrared thermography measurement system is facing an undetermined problem with one equation and two unknowns.

Multi-band thermography

Multiband sensors provide a measurement set of spectrally integrated radiance in various finitebandwidth spectral channels. An idea that could come into mind is to perform a measurement at another wavelength in order to bring new information to the system. However, new unknowns are also added when working in multi-band thermography, the spectral emissivity at those supplementary wavelengths [START_REF] Barducci | Temperature and Emissivity Retrieval from Remotely Sensed Images Using the" Grey Body Emissivity" Method[END_REF]. Let N be the number of spectral bands, i ∈ N, i ∈ [1, N ] and ∆λ i be the median value of the i th spectral band:

                     ln ( E ∆λ 1 ,sensor ) + ln ( ∆λ 5 1 C 1 g(θsensor,θ object ,r) ) = ln ( ϵ ∆λ 1 ,object ) - C 2 ∆λ 1 T object ln ( E ∆λ 2 ,sensor ) + ln ( ∆λ 5 2 C 1 g(θsensor,θ object ,r) ) = ln ( ϵ ∆λ 2 ,object ) - C 2 ∆λ 2 T object ... ln ( E ∆λ N ,sensor ) + ln ( ∆λ 5 N C 1 g(θsensor,θ object ,r) ) = ln ( ϵ ∆λ N ,object ) - C 2 ∆λ N T object (2.58)
As a consequence the system is still undetermined and is composed of N equations for N + 1 unknowns. Even if the equation system is unknown, multispectral measurements did brought information. Indeed, the infrared energy radiated by an object is proporional to its absolute temperature T object and emissivity ϵ object . However, sources of stray radiation (background, sky, sun and so on) having a temperature T env and emissivity ϵ env can add, after reflection a disturbing radiative contribution to the previous intrinsic emission and reflectance of the object. Of course, an infrared imager will be sensitive to the sum of all those energies and will become a function of all the parameters that will interfere. In the case of the grey-body assumption, a measurement on a single spectral band will give one relation between T obj and ϵ obj . Multispectral measurements can therefore provide more information and particularly when ϵ obj is a function of wavelength. For example, if the same radiance is produced in one band (∆λ 2 in Fig. 2.20) because of a lower emissivity compensated by a higher temperature; it will probably not be the case for another band (∆λ 1 ):

Synthesis

Infrared thermography for structural health monitoring has been introduced in this chapter. Equations that govern the physical phenomena of radiative transfers have been derived and studied for the particular case of IR thermography measurements. It has been shown that the IR thermography measurement process involves many variables that need to be taken into account in the hope of getting accurate temperature measurements. In particular, the emissivity plays an important role in the equations. Without emissivity, the temperature estimation uncertainty may be high and therefore, emissivity-temperature simultaneous estimation methods are needed. The objective of this thesis is to develop such method, applied to in-situ long-term thermal monitoring of infrastructures. To achieve this, the understanding of the full measurement process is needed. In particular and as shown in the former sections, the knowledge of the environmental parameters helps at improving the accuracy of the temperature estimation. Moreover, the infrared camera needs to be calibrated both spatially and thermally. These topics will be discussed in Chapter 4. Furthermore, it has been shown that the set-up of the instrumentation will have an impact on the obtained measurement and therefore the potential temperature estimate. In order to optimize such instrumentation but also to test the emissivity and In the next chapter, the methods developed in the literature for separating the emissivity and the temperature will be presented. Then, a bibliographical study about solving the radiosity equation will also be done in view of implementing a custom simulation software.

Chapter 3

Bibliographical study I n this chapter we will review the problem of extracting the temperature and the emissivity information from thermal infrared multispectral radiance data. We will present the different methods that exist in the litterature and discuss them. Being able to identify the emissivity and the temperature at the same time while knowing the approximation made on the measurement is a challenge that could enhance infrared data exploitation and infrared cameras applications possibilities. We will focus here on the simultaneous evaluation of emissivity and surface temperature of objects observed by infrared thermography, for in-situ applications.

Moreover, in order to test the efficiency, and robustness of those methods but also to optimize the in-situ instrumentation, a simulation software tool using the forward radiosity equation solution has been developed. Thanks to literature's models, infrared images can be simulated from 3D models. A brief state of the art of illumination algorithms to render such images is therefore presented in this chapter. As explained previously, the temperature cannot be deduced from IR measurements without a good emissivity estimate.

Emissivity and temperature separation methods

As seen before, the lack of knowledge on radiation properties of the real scene makes the thermal measurement through infrared thermography complicated. In first approximation, the irradiance received by a sensor for a given wavelength depends on the objects' temperature and its emissivity at this wavelength, as explained in the previous chapter. However, the radiative system of such configuration leads to an undetermined system. Therefore, retrieving the temperature from such a measurement needs the knowledge of the object, or by default, the use of inferences and approximations on the observed system. The methods are classified into two categories. A first category regroup the optimization based methods that are derived from a mathematical approximation at some point to reduce the number of variables or to increase the redundancy of data to get more equations. The second one gathers the Bayesian methods for which a prior statistical information is added to the system in order to be able to solve it. The literature has different models for estimating the emissivity from remotely sensed images and most of those models are used for computing satellite data. The main difficulty here is to separate the two unknowns of the equation: the emissivity and the temperature. Even if the equations behind airborne applications and in-situ applications are close, the conditions and assump-tions to be applied may differ, in particular when considering the atmospheric contribution. However, the methods used in the remote field sensing can be inspiring for in-situ long-term thermal monitoring. [START_REF] Krapez | Radiative Measurements of Temperature[END_REF] is a good introduction for having a global point of view on general noncontact temperature measurement by radiative methods.

In the field of Remote Sensing

Land surface temperature (LST) is an important parameter for various applications. Such indicator provides information on the temporal and spatial variations of the surface equilibrium state and is used for vegetation monitoring, hydrological cycle, urban and environmental studies or climate change for example. Many methods have been developed to retrieve the LST from satellite infrared data. The radiative transfer equation for measurements performed at the top of the atmosphere is different than the one used for in-situ measurements and presented in the previous chapter. Those measurements have to take into account the atmospheric path that the beam has taken until it reaches the satellite sensor. Even though the model is different, it is still interesting to take a look at the developed emissivity / temperature separations methods and see if it can be applied to in-situ long-term IRT on ground.

The measured radiance at sensor level can be written [START_REF] Li | Satellite-Derived Land Surface Temperature: Current Status and Perspectives[END_REF]:

L λ (θ, ϕ, T ) = τ λ (θ, ϕ) ( ϵ λ (θ, ϕ)L • λ,land (T ) + (1 -ϵ λ (θ, ϕ))(L ↓ λ,a + L ↓ λ,sun ) + Lλ,sun ) +L ↑ λ,atm + L↑ λ,atm (3.1)
where (see Fig. 3.1):

• τ λ (θ, ϕ) is the atmospheric transmission,

• ϵ λ (θ, ϕ) is the land surface emissivity,

• T is the land surface temperature,

• L • λ,land (T ) is the land surface emission, • L ↓
λ,atm is the down-welling atmospheric radiance reflected by the surface, • L ↓ λ,sun is the down-welling diffuse solar radiance reflected by the surface, -L λ,sun is the direct solar radiance,

• Lλ,sun is the direct solar radiance reflected by the surface,

• L ↑ λ,atm is the atmospheric emission, • L↑ λ,atm is the atmospheric scattering.
and Lλ,sun is defined by:

Lλ,sun = ρ λ (θ, ϕ, θ s , ϕ s )L λ,sun cos(θ s )τ λ (θ s , ϕ s ) (3.2)
For simplification in the writting, the formalism of the Remote Sensing field will be used in this section. In fact, due to the narrow bands involved, λ notation is used in place of ∆λ. Moreover, when dealing with satellite data, the atmospheric effects along the path from the ground to the sensor are assumed to be properly estimated. This case makes the emissivity retrieval close to the problem that occurs on land measurements (N measurements for N + 1 unknowns). In order to get rid of those parameters, atmospheric radiative transfer models such as SMARTS2 ( [START_REF] Gueymard | Simple Model for the Atmospheric Radiative Transfer of Sunshine (SMARTS2) Algorithms and Performance Assessment[END_REF]), MODTRAN ( [START_REF] Berk | MODTRAN: A Moderate Resolution Model for LOWTRAN[END_REF]) or MATISSE ( [START_REF] Fauqueux | MATISSE-v2.0 Infrared Sea Images: Sensitivity Analysis and Experimental Validation Using MIRAMER Campaign Measurements in Solar Glint Configurations[END_REF]) can be used. In this field, the unknowns are referred to LSE and LST for land surface emissivity and land surface temperature reciprocally. Many of the temperature and emissivity separation methods presented here are introduced and discussed in [START_REF] Li | Satellite-Derived Land Surface Temperature: Current Status and Perspectives[END_REF] state of the art.

Note

Known emissivity

The following methods are temperature -emissivity separation methods but those methods do not actually estimate conjointly the temperature and emissivity. They need prior knowledge on LSE [START_REF] Li | Satellite-Derived Land Surface Temperature: Current Status and Perspectives[END_REF].

Those methods are still presented here for their data information extraction from IR measurements.

Thermal log residuals By using Wien's approximation of Planck's law Eq. (2.8) the log can be expressed as:

ln   ϵ × C1 πλ 5 exp ( C 2 λT )   = ln(ϵ) + ln(C 1 ) -5 ln(λ) -ln(π) - C 2 λT (3.3)
By multiplying this equation by λ, we can dissociate T from λ. Let X pc be the resulting quantity for the pixel p and the channel c of a multispectral infrared sensor.

X pc = λ c ln(ϵ pc ) + λ c ln(C 1 ) -5λ c ln(λ) -λ c ln(π) - C 2 T p (3.4)
Then, the "residuals" (Y pc ) are defined as the remaining of this quantity after having substracted the average value over the image ( Xc ), the average value over the channels ( Xp ) and the average value over both images and channels ( X):

Y pc = X pc -Xp -Xc + X (3.5)
With N the number of spectral channels and M the number of pixels in one image:

Xp = 1 N N ∑ j=1 X pj Xc = 1 M M ∑ i=1 X ic X = 1 N M M ∑ i=1 N ∑ j=1 X ij (3.6)
From equations Eq. (3.4) and Eq. (3.6):

Y pc = λ c ln(ϵ pc ) - 1 N N ∑ j=1 λ j ln(ϵ pj ) - λ c M M ∑ i=1 ln(ϵ ic ) + 1 N M N ∑ j=1 M ∑ i=1 λ j ln(ϵ ij ) (3.7)
Finally, let l c = λc λ and ϵ pc = e -lc pc where λ is the mean of the central wavelengths of all channels, then:

Y pc λ = ln(e pc ) - 1 N N ∑ j=1 ln(e pj ) - 1 M M ∑ i=1 ln(e ic ) + 1 N M N ∑ j=1 M ∑ i=1 ln(e ij ) = ln(e pc ) - 1 N ln   N ∏ j=1 e pj   - 1 M ln ( M ∏ i=1 e ic ) + 1 N M ln   N ∏ j=1 M ∏ i=1 e ij   (3.8)
Therefore, the thermal log residual enables to separate the variation of temperature and emissivity to a variation of combined emissivity.

Alpha residuals

Similarly to thermal log residuals, the alpha log residuals performs a mean substraction but only with the wavelength-weighted log pixel this time [START_REF] Simon | A Comparison of Techniques for Extracting Emissivity Information from Thermal Infrared Data for Geologic Studies[END_REF].

α pc = X pc -Xp (3.9)
Therefore, let d :

x c → x c -1 N ∑ N c=1
x c , the alpha log residuals equals to:

α pc = d(λ c ln(ϵ pc )) + (ln(C 1 ) -ln(π))d(λ c ) -5d(λ c ln(λ c )) Kc (3.10)
Where K c is a constant known for a given channel and depends only on the wavelength.

Let r be one reference channel. Therefore, it is possible to write:

ϵ pc ϵ nc pr = e 1 λc [(αpc+Kc)-(αpr+Kr)] , n c = λ c λ r (3.11)
The same authors showed that

α pc = d(λ c ln(B pc )) + K c (3.12)
Meaning that the alpha residuals can be directly computed from the measured radiance.

Those two methods (alpha log and thermal log residuals) have the advantages of being simple and applicable to one channel sensors only. However, some priors on the emissivity are required and the important source of errors in those methods will be the accuracy of the emissivity pixel's estimate.

Moreover, the uncertainty on atmospheric estimated parameters will have a strong effect on LST retrieval [START_REF] Li | Satellite-Derived Land Surface Temperature: Current Status and Perspectives[END_REF].

Other methods with known emissivity

Other methods with known emissivity rather than the log residuals have been developed. Those methods use different assumptions to improve the accuracy of the estimated LST or, more generally to reduce the strong knowledge requirement of known emissivity at each pixel. Multi-channel algorithms for example use different atmospheric absorption in adjacent channels and reduce the need for accurate atmospheric profiles [START_REF] Becker | Towards a Local Split Window Method over Land Surfaces[END_REF][START_REF]A Generalized Split-Window Algorithm for Retrieving Land-Surface Temperature from Space[END_REF][START_REF] Sobrino | Improvements in the Split-Window Technique for Land Surface Temperature Determination[END_REF][START_REF] Larry | Estimation of Sea Surface Temperatures from Two Infrared Window Measurements with Different Absorption[END_REF].

Similarly, another way is to use the angle dependence of the atmospheric parameters. A change of the angle will change the atmospheric absorption due to the different path-lengths from different views. Again, such method reduces the need for accurate atmospheric profiles. However, the angular emissivity variation with angle needs to be known and also geometric registration becomes fundamental [START_REF] Sòria | ENVISAT/AATSR Derived Land Surface Temperature over a Heterogeneous Region[END_REF][START_REF] Sobrino | Multi-Channel and Multi-Angle Algorithms for Estimating Sea and Land Surface Temperature with ATSR Data[END_REF][START_REF]Land Surface Temperatures Derived from the Advanced Very High Resolution Radiometer and the Along-Track Scanning Radiometer 1. Theory[END_REF].

Unknown emissivity -two steps-methods

The previous presented methods were assuming that the LSE was known at some point. However, in practice, knowing LSE is challenging due to the heterogeneity of the surfaces and its angular and spectral variations. Moreover, an uncertainty on the emissivity leads to an important uncertainty on the temperature estimation (see next chapter for a sensitivity study regarding in-situ long-term monitoring). Therefore, methods for estimating both LSE and LST have been developed in the literature. In [START_REF] Li | Satellite-Derived Land Surface Temperature: Current Status and Perspectives[END_REF] those methods are classified into two groups: the two-steps methods and the simultaneous estimation of LSE and LST. The first group methods estimate first the LSE and then use one of the former method to estimate the LST. At the opposite, the simultaneous estimation methods aims at providing a-priori on the system to estimate at the same time LSE and LST.

Temperature Independent Spectral Indices (TISI) T ISI coefficients have been proposed for spectral analysis in the TIR (thermal infrared, centered at 8.9µm) band. Those coefficients aim at getting rid of surface temperature in their expression to only express a variation of emissivity through a reference channel. Those coefficients are called Temperature-Independent Spectral Indices (T ISI). Please note that T ISI theory and algebra have been developed in the original article and different application examples are provided in [START_REF] Becker | Temperature-Independent Spectral Indices in Thermal Infrared Bands[END_REF]. Based on those coefficients, a new two-steps method has been proposed by assuming that T ISI coefficients on particular channels will not vary substantially between day and night [START_REF] Li | Feasibility of Land Surface Temperature and Emissivity Determination from AVHRR Data[END_REF][START_REF] Li | A Physically Based Algorithm for Land Surface Emissivity Retrieval from Combined Mid-Infrared and Thermal Infrared Data[END_REF]. The T ISI ij are referred by two channels i and j, usually i is the MIR (mid infrared, centered at 3.8µm) channel and j the TIR (thermal infrared, centered at 8.9µm) reference channel.

T ISI ij are defined by using the power law approximation of Planck's function:

L • i (T ) = a i T n i (3.13)
where a i and n i are channel-specific constants. T ISI ij indices are therefore defined by:

T ISI ij = L • i (T gi ) -L ↓ i,a L • i (T gj ) -L ↓ i,a (3.14)
which they shown mathematically that [START_REF] Becker | Temperature-Independent Spectral Indices in Thermal Infrared Bands[END_REF]:

T ISI ij ≈ ϵ i ϵ n i /n j j (3.15)
Those spectral indices are therefore independent of ground surface temperature and are an indicator of the LSE evolution in the two channels. By assuming that the T ISI variation can be neglected between daytime without the contribution of solar illumination and nighttime then the bidirectional reflectivity in MIR channel is found back by removing the emitted radiance during the day and comparing it with night T ISI ij coefficients. LSEs in the TIR channels are then deducted based on the T ISI and emissivity in the MIR channel. Former methods are then used to estimate the LST.

T ISI indices have been used in different methods and have shown satisfying overall results [START_REF] Li | Evaluation of Six Methods for Extracting Relative Emissivity Spectra from Thermal Infrared Images[END_REF].

However, due to the constant emissivity ratio during day and night assumption, a strong requirement is to have the same observation view during night and day with accurate geometry registration. Furthermore, this method works only for one MIR channel and one TIR channel at least.

Unknown emissivity -simultaneous estimation of LSE and LST

The following methods aim at estimating conjointly LSE and LST. This time, the methods provide directly both estimates, without estimating first the LSE and then the LST as in two-steps methods.

Normalized-Emissivity separation method (NES)

Let assume that the maximum emissivity ϵ max is known, without knowing especially its associated wavelength or band. Firstly, a relatively high temperature T is chosen. The following radiometric model is then solved for each spectral band at each pixel to get a first estimate of the emissivity :

ελ = L λ,land (T ) -L ↓ λ L • λ,land ( T ) -L ↓ λ (3.16)
Where L ↓ λ = L ↓ λ,a + L ↓ λ,s and L λ,land is the ground-leaving radiance. The temperature T is then decreased respectively increased until ελ = ϵ max if L λ,land (T ) > L ↓ λ , respectively L λ,land (T ) < L ↓ λ . Another way to achieve such separation is to make the gray body assumption where ϵ λ = ϵ max . In this case, the temperature can be found when L λ,land (T ) > L ↓ λ :

T = max λ [ L •-1 ( L λ,land (T ) -L ↓ λ ϵ max + L ↓ λ )]
(3.17)

When L λ,land (T ) < L ↓ λ , a minimum must be found instead. Once T has been estimated, the emissivity spectrum is deduced from Eq. (3.16).

According to [START_REF] Krapez | Radiative Measurements of Temperature[END_REF], in the case of bare soils with rocks, this method is quite precise (0.02 emissivity uncertainty, ϵ max = 0.96). However, when the surface is unknown, the underlying hypothesis is not satisfied and the inacurracy increase and becomes finally high (±3K).

Temperature-Emissivity separation method (TES)

In order to overcome the high inaccuracies induced by the NES method, the TES method has been developed. Proposed by [START_REF]A Temperature and Emissivity Separation Algorithm for Advanced Spaceborne Thermal Emission and Reflection Radiometer[END_REF], the TES method is based on the fact that authors have observed that the relative emissivity of β λ = ελ / ε is rather insensitive to the error made on the temperature estimation [START_REF] Krapez | Radiative Measurements of Temperature[END_REF]. ελ is derived from Eq. (3.16) and ε its average on wavelength. In order to obtain the absolute spectrum from this relative one, [START_REF]A Temperature and Emissivity Separation Algorithm for Advanced Spaceborne Thermal Emission and Reflection Radiometer[END_REF] proposed a correlation:

ϵ λ = β λ ϵ min β min , ϵ min ≈ A -B(β max -β min ) C (3.18)
The coefficients A, B and C are obtained through a regression depending on a material database. In their article, A = 0.994, B = 0.687 and C = 0.737 but other values can be found in the literature.

Finally, the temperature is derived from:

T = L •-1   L λ M (T ) -L ↓ λ M ϵ max + L ↓ λ M   , λ M = arg max λ (ϵ λ ) (3.19)
TES method is used to calclulate surface temperature for ASTER satelite. The TES has been validated on different sites, showing that TES performs well for ASTER data (+1.5K and +0.015 of correct values) [START_REF] Sabol | Field Validation of the ASTER TemperatureEmissivity Separation Algorithm[END_REF]. However, the correlation is highly dependent on the chosen database materials and thus the generalization to any site and any scale of this method seems tricky [START_REF] Gillespie | Residual Errors in ASTER Temperature and Emissivity Standard Products AST08 and AST05[END_REF].

Classification based methods Classification methods are widespread for identifying parameters in images. In order to derive LSE from the measurements, a classification algorithm is performed to classify a given pixel to a given class. By using the right database these methods offer good accuracy results (about 0.01) [START_REF] Snyder | Classification-Based Emissivity for Land Surface Temperature Measurement from Space[END_REF]. Once the LSE has been estimated, LST is derived from previous methods. However, the classification algorithm needs a strong knowledge of the emissivity of each class as well as robust and reliable classification map which becomes complicated when considering soil moisture, seasoning, snow etc. Moreover, classification is also more difficult to perform for coarse pixels, due to the mixed materials observation and then strong correlation within classes [START_REF] Gillespie | Temperature/Emissivity Separation Algorithm Theoretical Basis Document[END_REF].

Other methods Please note that this method list is not comprehensive, and other methods could be quoted. For example, the spectra smoothness method which states that if LST is not well estimated, then saw-teeth will appear on the estimated LSE signal due to atmospheric absorption lines [START_REF] Li | Satellite-Derived Land Surface Temperature: Current Status and Perspectives[END_REF].

When the estimated LST comes closer to the real temperature, the estimated LSE becomes progressively smoother. The estimate is then found when LSE is the smoothest. This method offers interesting results for satelite and airborne applications. However, the underlying hypothesis is to use the atmospherical spectrum properties to have a better estimate of the temperature and emissivity. Furthermore, this method requires hyperspectral instruments and a good high spectral resolution in order to be able to distinguish the atmospheric spectral features from the emissivity features. Such method is therefore out of scope for in-situ multispectral measurements.

General separation methods

Whereas the previous section introduced Remote Sensing field specific methods, this section presents a more general overview of the temperature and emissivity separation problem, from a mathe-matical basis. Whatever the field of application, once the environmental or atmospheric parameters have been identified, it is possible to bring the model back to the case:

γ ∆λ i = f (∆λ i , ϵ ∆λ i , T )
where f is known without accounting for specific channels or atmospheric properties to estimate the temperature and emissivity. Therefore, some presented methods such as the grey body and multi-temperature methods have been actually developed in the Remote Sensing field. In the following, the conjoint estimation of emissivity and temperature will be studied, based on the sensor -object Eq. (2.56).

Linearization of the problem

Let γ ∆λ i be g -1 (θ sensor , θ object , r)C -1 1 ∆λ i 5 E ∆λ i ,sensor . This quantity is assumed to be known in the following, representing the measurements for which we have prior knowledge about the geometry.

As a reminder, we want to solve the system:

ln(γ ∆λ i ) = ln(ϵ ∆λ i ) - C 2 ∆λ i T , ∀i ∈ [1; N ] (3.20)

Monomial model

In order to solve this underdetermined problem, a first approach is to estimate ln(ϵ ∆λ i ) by a polynom. Let (α 0 , α 1 , ..., α m ) ∈ R m+1 where m ≤ N -2 such as:

ln(ϵ ∆λ i ) = m ∑ k=0 α k ∆λ i k , ∀i ∈ [|1; N |] (3.21)
Such approximation leads to the following equation system:

Γ = M X =        1 ∆λ 1 ∆λ 1 2 . . . ∆λ 1 m -C 2 ∆λ 1 1 ∆λ 2 ∆λ 2 2 . . . ∆λ 2 m -C 2 ∆λ 2 . . . . . . . . . . . . . . . . . . 1 ∆λ N ∆λ N 2 . . . ∆λ N m -C 2 ∆λ N                     α 0 α 1 α 2 . . . α m 1 T ′              (3.22)
Where Γ = (ln(γ ∆λ i )) 1≤i≤N and M is a [N × (m + 2)] matrix and T ′ is the estimated temperature obtained from the polynomial approximation and measurement. In the following, we will remove the ∆ notation inferring that the sensor integrates on a small band. Then, the λ i notation will refer to a band and not a particular wavelength.

Monomial model (m = N -2)

Let study the case where we have exactly N equations for N unknowns. This is possible by using a N -2 order polynomial. By rewriting Eq. (3.22), we get a square matrix on which we recognize the shape of a Vandermonde matrix and we can compute the determinant:

det(M ) = C 2 (-1) N +1 ∏ 1≤i<j≤N (λ i -λ j ) ∏ N i=1 λ i ≤ C 2 (-1) N +1 ∆λ N max ∏ N i=1 λ i (3.23)
where ∆λ max is the maximum difference of two bands. We can see that the absolute value of this determinant decreases exponentially as N increases, meaning that M is approaching singularity when the number of bands increases. The literature shows that the polynomial approximation leads to unrealistics temperatures for N > 3 [START_REF] Coates | Multi-Wavelength Pyrometry[END_REF]. Such behavior can be explained by the fact that even if the Taylors serie development of the ln function enables us to estimate the logarithm in Eq. (3.21) because of the boundaries of ϵ, we need to get a N large enough to get a good approximation as shown on Fig. 3.2. Therefore, for a N too small, the variables will be considered sufficiently independents but the approximation is rough. At the opposite, when N increases, the approximation of the logarithm becomes satisfying but we induce correlation between the different variables which makes the system difficult to inverse.

λ i ln(ϵ i ) = N -2 ∑ j=0 a j-1 λ j i (3.24)
where a -1 is the constant term evaluated on λ i = 0, giving an idea of the error made on the temperature, from Eq. (3.20):

a -1 = C 2 ( 1 T - 1 T ′ ) (3.25)
Since the polynom will pass exactly through the N values of λ i ln(ϵ i ), we can use Lagrange polynomials to find the values

p(λ) = N -2 ∑ j=0 a j-1 λ j = N ∑ i=1   λ i ln(ϵ i ) N ∏ j=1,j̸ =i λ j -λ λ j -λ i   (3.26)
With Eq. (3.18), we get the value of a -1 [START_REF] Coates | Multi-Wavelength Pyrometry[END_REF]:

a -1 = C 2 ( 1 T - 1 T ′ ) = N ∏ j=1 λ j N ∑ i=1 [ ln(ϵ i ) ∏ N j=1,j̸ =i (λ j -λ i ) ] (3.27)
Therefore, only a polynom that will exactly match the N values of ln(ϵ i ) will have a null extrapolation error and be the real emissivity; which is unlikely to happen.

Based on those equations, a particular study of the different pyrometry techniques (two, three and four colors) can be done ( [START_REF] Mansoor | Noncontact Temperature Measurement. I. Interpolation Based Techniques[END_REF][START_REF] Coates | Multi-Wavelength Pyrometry[END_REF][START_REF] Helcio | Thermal Measurements and Inverse Techniques, Heat Transfer[END_REF]). We can observe the conditions that are necessaries to get a null error (id est solving α -1 = 0).

Let λ i ̸ = λ j , ∀i, j ∈ [|1, N |], i ̸ = j N = 2 C 2 ( 1
T - 1 T ′ ) = λ 1 λ 2 λ 1 -λ 2 ln( ϵ 2 ϵ 1 ) = 0 ⇐⇒ ϵ 1 = ϵ 2 N = 3 C 2 ( 1 T - 1 T ′ ) = 0 ⇐⇒ ϵ 1 ϵ 3 = ϵ 2 2 (equidistant bands) (3.28)
Those results show that strong hypothesis are needed on the system we use that may lead to uncertain measurements without prior knowledge on emissivities. Finally, it is well known that overfitting of the orignal data on interpolation based methods leads to errors (Runge's phenomenon). Such methods can be retained only for simple pyrometers (two or three wavelengths at most) [START_REF] Coates | Multi-Wavelength Pyrometry[END_REF].

To overcome the previous overfitting issue, one may try to reduce the order of the polynom used to fit the data. In that case, the estimated temperature value is obtained by a least-square method, to solve the overdetermined problem.

Monomial model (m < N -2)

To overcome some numerical issues, it is necessary to normalize the system matrix in order to reduce its condition number. This can be achieved by replacing the wavelength band in the previous polynomial expression by its reduced centered value such as:

λ * i ∈ [-1; 1]
[135]:

λ * i = 2 λ i -λ min λ max -λ min (3.29)
For the same reasons, we normalize the temperature with T ref such as C 2 λ i T ref is close to one. Therefore, we can write the new normalized system:

Γ = M X =         1 λ * 1 λ * 1 2 . . . λ * 1 m -C 2 λ * 1 T ref 1 λ * 2 λ * 2 2 . . . λ * 2 m -C 2 λ * 2 T ref . . . . . . . . . . . . . . . . . . 1 λ * N λ * N 2 . . . λ * N m -C 2 λ * N T ref                      α * 0 α * 1 α * 2 . . . α * m T ref T ′              (3.30)
In that case, the least-square method resolution leads to solve the normal equations, to minimize the sum of the square differences .

(

M T M )X = M T Γ (3.31)
Many methods exist in the literature to solve such equation system such as the Singular Value Decomposition (SVD), the QR decomposition, etc. [START_REF] Nath | Numerical Methods for Linear Control Systems: Design and Analysis[END_REF] Orthogonal basis (m < N -2) The approximation of taking the logarithm as a polynomial could be improved by choosing the right polynomial basis instead of the previous presented basis function {P 0 , P 1 , ..., P N }. In fact, it is clear that the approximation of a function to a polynomial is like a projection on the space spanned by the basis function. Using orthogonal polynomials helps to have uncorrelated explanatory variables. Let {Φ 0 , Φ 1 , ..., Φ m } an orthogonal set of functions of R + * and scalars (α 0 , α 1 ..., α m ) such as:

Γ = M X =         Φ 0 (∆λ * 1 ) Φ 1 (∆λ * 1 ) . . . Φ m (∆λ * 1 ) -C 2 λ 1 T ref Φ 0 (∆λ * 2 ) Φ 1 (∆λ * 2 ) . . . Φ m (∆λ * 2 ) -C 2 λ 2 T ref . . . . . . . . . . . . . . . Φ 0 (∆λ * N ) Φ 1 (∆λ * N ) . . . Φ m (∆λ * N ) -C 2 λ N T ref                   α 0 α 1 . . . α m T ref T ′           (3.32)
By combining equations Eq. (3.31) and Eq. (3.32), one can see that the orthogonal basis will increase diagonal coefficients and reduce other coefficients :

M T M =         ∑ N i=1 Φ 0 2 (∆λ * i ) ∑ N i=1 Φ 0 (∆λ * i )Φ 1 (∆λ * i ) ... ∑ N i=1 Φ 0 (∆λ * i )Φm(∆λ * i ) - ∑ N i=1 Φ 0 (∆λ * i ) C 2 λ i T ref ∑ N i=1 Φ 0 (∆λ * i )Φ 1 (∆λ * i ) ∑ N i=1 Φ 1 2 (∆λ * i ) ... ∑ N i=1 Φ 1 (∆λ * i )Φm(∆λ * i ) - ∑ N i=1 Φ 1 (∆λ * i ) C 2 λ i T ref . . . . . . . . . . . . . . . - ∑ N i=1 Φ 0 (∆λ * i ) C 2 λ i T ref - ∑ N i=1 Φ 1 (∆λ * i ) C 2 λ i T ref ... - ∑ N i=1 Φm(∆λ * i ) C 2 λ i T ref ∑ N i=1 ( C 2 λ i T ref ) 2         (3.33)
In order to illustrate the effect of the orthogonalization onto the final estimation, a simple usecase has been experimented. The Eq. (3.20) is computed for ∆λ i = {2.5µm, 3µm, 4µm, 5µm, 8µm, 10µm, 12µm}, a given temperature profile (represented in Fig. 3.4), random emissivity values, θ sensor = θ object = 0 and r = 1. Once the data have been simulated, the emissivities and temperatures are found back with the polynomial model presented in Eq. (3.30) and the Tchebychev polynomials of first kind applied to the system Eq. (3.32). The systems are solved with the constraint : ln(ϵ ∆λ i ) ≤ 0, ∀i ∈ [1, N ]. Moreover, the emissivity profile is generated randomly. The purpose of this is to test the robustness of the fitting on a random profile. The inversion becomes more difficult when the standard deviation of the emissivity profile is larger, due to the important discontinuity that appears and the small amount of band measurements. The Fig. 3.3 shows the simulated irradiance in comparison to the irradiance computed with the estimated parameters. One can see that this reprojection error can be important since the algorithm stops when a local minimum has been found. However, the Fig. 3.4 shows that the estimation is slightly different between the Tchebychev basis and the monomial basis. The orthogonal basis tends to be more robust than the monomial's one, as expected. However, the results are not quite satysfying in the simple test case, even when using an orthogonal basis. In fact, the polynomial approximation is also depending on the temperature estimation sensitivity [START_REF] Krapez | Radiative Measurements of Temperature[END_REF].

Splines

In the previous approximation procedures, the emissivity was evaluated to a single identical polynomial. It is possible instead to split the domain into different region and define a polynomial's function piecewise. Such function is called a spline and is often prefered to the previous methods for interpolation problems. Let now S j (λ), j ∈ [1, N -1] be the polynoms that interpolate the ln(ϵ ∆λ i ), i ∈ [1, N ] knots values and m the degree of the polynoms.

S j = m ∑ k=0 α j,k Φ k (λ * ), ∀j ∈ [1, N -1]
(3.34) To satisfy the knot's equality :

(C 0 ) =                                S 1 (∆λ * 1 ) = ln ( ϵ ∆λ * 1 ) S 1 (∆λ * 2 ) = ln ( ϵ ∆λ * 2 ) S 2 (∆λ * 2 ) = ln ( ϵ ∆λ * 2 ) S 2 (∆λ * 3 ) = ln ( ϵ ∆λ * 3 ) . . . S (N -1) (∆λ * N ) = ln ( ϵ ∆λ * N ) (3.35)
The Eq. (3.35) implies 2(N -1) equations, to satisfy the knots equalities. In order to be able to solve the system and reduce the number of unknowns, one needs to impose conditions to the system. In the case of cubic interpolation (m = 3), the number of unknowns is : 4(N -1) + 1 (4 coefficients for N -1 polynoms and the temperature). The constraint of equality of the first derivatives at knot's points can be added, to smooth the curves and avoid discontinuities. Therefore, we have the second equation system:

(C 1 ) =                  S ′ 1 (∆λ * 2 ) = S ′ 2 (∆λ * 2 ) S ′ 2 (∆λ * 3 ) = S ′ 3 (∆λ * 3 ) . . . S ′ N -2 (∆λ * N -1 ) = S ′ N -1 (∆λ * N -1 ) (3.36)
Moreover, one can also add the same constraint for the second derivatives:

(C 2 ) =                  S ′′ 1 (∆λ * 2 ) = S ′′ 2 (∆λ * 2 ) S ′′ 2 (∆λ * 3 ) = S ′′ 3 (∆λ * 3 ) . . . S ′′ N -2 (∆λ * N -1 ) = S ′′ N -1 (∆λ * N -1 ) (3.37)
Finally, (C 0 ) + (C 1 ) + (C 2 ) constraints lead to 4N -6 equations, for 4(N -1) + 1 unknowns, an undetermined system. To overcome this, boundaries constraints can be added. For example, a quadratic boundary constraint and a second derivative equal to zero for the first and last polynomials will result in a overdetermined system. Let modify the first and last polynomial to be quadratic. Eq. (3.34) becomes :

S j (λ * ) =    ∑ 2 k=0 α j,k Φ k (λ * ) if j ∈ {1, N -1} ∑ 3 k=0 α j,k Φ k (λ * ) if j ∈ [2, N -2] (3.38)
We also add their first and second derivatives to be equal at the boundaries (called "Periodic Spline"):

   S ′ 1 (∆λ * 1 ) = S ′ N -1 (∆λ * N ) S ′′ 1 (∆λ * 1 ) = S ′′ N -1 (∆λ * N ) =⇒    α 1,1 = α N -1,1 α 1,2 = α N -1,2 (3.39)
In that case, the number of unknowns becomes 4N -5 for 4(N -1) equations. The system Eq. (3.32) to be solved can be updated with the constraints :

                                                                                                     Γ * =         ln ( γ ∆λ 1 ) ln ( γ ∆λ 2 ) ln ( γ ∆λ 2 ) . . . ln ( γ ∆λ N )         =          Φ 0 (∆λ * 1 ) Φ 1 (∆λ * 1 ) Φ 2 (∆λ * 1 ) 0 0 . . . 0 -C 2 ∆λ 1 T ref Φ 0 (∆λ * 2 ) Φ 1 (∆λ * 2 ) Φ 2 (∆λ * 2 ) 0 0 . . . 0 -C 2 ∆λ 2 T ref 0 0 0 Φ 0 (∆λ * 2 ) Φ 1 (∆λ * 2 ) . . . 0 -C 2 ∆λ 2 T ref . . . . . . . . . . . . . . . . . . . . . . . . 0 0 0 0 0 . . . Φ 2 (∆λ * N ) -C 2 ∆λ N T ref          M 0         α 1,0 α 1,1 . . . α N -1,2 T ref T ′         0 =       Φ ′ 0 (∆λ * 2 ) Φ ′ 1 (∆λ * 2 ) Φ ′ 2 (∆λ * 2 ) -Φ ′ 0 (∆λ * 2 ) . . . 0 0 0 0 0 Φ ′ 0 (∆λ * 3 ) . . . 0 0 . . . . . . . . . . . . . . . 0 . . . 0 0 0 0 . . . -Φ ′ 3 (∆λ * N -1 ) 0       M 1         α 1,0 α 1,1 . . . α N -1,2 T ref T ′         0 =       Φ ′′ 0 (∆λ * 2 ) Φ ′′ 1 (∆λ * 2 ) Φ ′′ 2 (∆λ * 2 ) -Φ ′′ 0 (∆λ * 2 ) -Φ ′′ 1 (∆λ * 2 ) . . . 0 0 0 0 0 Φ ′′ 0 (∆λ * 3 ) Φ ′′ 1 (∆λ * 3 ) . . . 0 0 . . . . . . . . . . . . . . . . . . 0 . . . 0 0 0 0 0 . . . -Φ ′′ 2 (∆λ * N -1 ) 0       M 2         α 1,0 α 1,1 . . . α N -1,2 T ref T ′         0 = ( 0 1 0 . . . 0 -1 0 0 0 0 1 . . . 0 0 -1 0 ) M 3         α 1,0 α 1,1 . . . α N -1,2 T ref T ′         (3.40)
The matrices M 0 , M 1 and M 2 and M 3 can be merged together to form the following linear system:

    Γ * 0 0     =        M 0 M 1 M 2 M 3                  α 1,0 α 1,1 . . . α N -1,2 T ref T ′           (3.41)
Using splines interpolation has the advantage of adding continuity constraints to the interpolating function. Due to the fact that the emissivity is a continuous function in wavelength, splines interpolation should be preferred to the simple polynomial fitting. approach (i.e. space-time) but applies only a temporal constraint in their article. The Kalman filter is used conjointly with the σ -SEVIRI model based on look-up tables. Such method is specific to the SEVIRI measurements and cannot be applied directly to the in-situ IRT measurements case. However, the Kalman Filter appears to be an interesting approach to perform parameters estimation on the emissivity/temperature retrieval problem. In particular, the temporal approach for long-term monitoring through a Kalman filter is used in Chapter 6 on which more details on the Kalman Filter are given. Likewise, particle filters which are also Bayesian filters have been used in the Remote sensing field. However, particle filter has not been particularly used for LST / LSE simultaneous retrieval to our knowledge. This can be particularly explained by the fact that satellite measurements represent a huge amount of data for which particle filtering may not be suitable due to its complexity. Yet, particle filtering and ensemble kalman filters have been used for data assimilation [START_REF] Li | Estimating Net Radiation Flux in the Tibetan Plateau by Assimilating MODIS LST Products with an Ensemble Kalman Filter and Particle Filter[END_REF] or data downscaling purposes [START_REF] Mechri | Genetic Particle Filter Application to Land Surface Temperature Downscaling[END_REF]. The particle filter is also more described in Chapter 6.

Monte-Carlo Markov Chains (MCMC)

The Monte-Carlo Markov Chains (MCMC) methods are meta-heuristic methods that sample probability distributions. The basic idea is to be able to browse a Markov Chain that will have the targeted distribution as equilibrium distribution. Such method can be employed for the temperature and emissivity retrieval [START_REF] Ash | Temperature-Emissivity Separation for LWIR Sensing Using MCMC[END_REF]. This method has been implemented for comparison purposes in this thesis and is therefore more detailed here. In this case, knowing the measurement vector γ, the objective is to determine the probability distributions of p(T |γ) and p(ϵ|γ). For this, let suppose that the emissivity ϵ = (ϵ λ i ) 1≤i≤N and temperature T are independent, and follow a given prior:

p(ϵ) = N (µ ϵ , Σ ϵ ) p(T ) = U(T min , T max ) (3.42)
The radiative equation used is derived from Eq. (2.50) where the optical transmission is assumed equal to 1 [START_REF] Berrett | A Bayesian Nonparametric Model for Temperature-Emissivity Separation of Long-Wave Hyperspectral Images[END_REF] and the sun contribution is included in the environmental contribution. A difference is made by considering the variables as stochastic and therefore adding a measurement noise η:

L (i,j) λ, ⃗ d,tot (T ) = τ atm ( ϵ λ L • λ (T (i,j) obj ) + (1 -ϵ λ )(L λ,env ) ) + (1 -τ atm )L • λ,atm + η (3.43)
This equation can be written:

γ = D T ϵ + c + η (3.44)
where

D T = diag(τ atm ( L • λ (T (i,j) obj ) -L λ,env ) and c = τ atm L λ,env + (1 -τ atm )L • λ,atm .
Finally, it is possible to derive the conditional distributions with this equation. Since γ is a linear function of a Gaussian variable ϵ plus noise η, then the distribution of ϵ knowing γ has also to be Gaussian [START_REF] Ash | Temperature-Emissivity Separation for LWIR Sensing Using MCMC[END_REF]. By using Gaussian conditioning:

p(ϵ|T, γ) = N (μ ϵ , Σϵ ) μϵ = µ ϵ + Σ ϵ D T (D T Σ ϵ D T + Σ η ) -1 (γ -D T µ ϵ -c) Σϵ = Σ ϵ -Σ ϵ D T (D T Σ ϵ D T + Σ η ) -1 D T Σ ϵ (3.45)
The temperature distribution is determined up-to proportionality:

p(T |ϵ, γ) ∝ p(γ|T, ϵ)p(T |ϵ) = p(γ|T, ϵ)p(T ) (3.46)
Since the joint density a posteriori of ϵ and T are difficult to derive and interpret, a Gibbs sampler will allow to sample the probability of p(ϵ, T |γ) based on the conditional probabilities p(ϵ|T, γ) and p(T |ϵ, γ).

Let π a distribution over a universe Ω. Drawing randomly an element of Ω following π can be done through a Gibbs sampler [START_REF]Explaining the Gibbs Sampler[END_REF]. This type of sampling is called a slice-within-Gibbs sampling. Even though this method provides quite accurate results and its convergence is ensured in this case, the required computation is heavy and processing time important. Such method can therefore be used for post-process or handling a small amount of data but the complexity of it makes it difficult to use on full resolution IR images and real-time processing. A comparison of MCMC method with other ones is done in Chapter 6.

Artificial Neural Network (ANN)

Artificial neural network can robustly perform non-linear and parallel computations and have therefore been largely used for remote sensing [START_REF] Aires | A New Neural Network Approach Including First Guess for Retrieval of Atmospheric Water Vapor, Cloud Liquid Water Path, Surface Temperature, and Emissivities over Land from Satellite Microwave Observations[END_REF][START_REF] Li | Satellite-Derived Land Surface Temperature: Current Status and Perspectives[END_REF][START_REF] Aires | Remote Sensing from the Infrared Atmospheric Sounding Interferometer Instrument 2. Simultaneous Retrieval of Temperature, Water Vapor, and Ozone Atmospheric Profiles[END_REF][START_REF]Retrieval of Atmospheric and Land Surface Parameters from Satellite-Based Thermal Infrared Hyperspectral Data Using a Neural Network Technique[END_REF]. A neuron represents a particular transfer function that will transform the inputs to outputs according to a particular rule (which can be summation, thresholding etc.). Neurons are then connected between them with different possible patterns.

Through statistical (usually Bayesian) learning methods based on a particular data-set, interconnected neurons store their acquired knowledge by updating their weight value (this behavior replicates the synaptic plasticity of biological networks). However, such methods are facing multiple difficulties such as overfitting which occurs when a given model is too dependent on training data and may not correspond to the actual targeted model. As a consequence, the selection of training data and their characteristics are essential for the method to give accurate results. Furthermore, the architecture used for the network is also one of the key-point [START_REF]The Application of Artificial Neural Networks to the Analysis of Remotely Sensed Data[END_REF]. Unfortunately, there is no standard way to help in the choice of those parameters. Most of the time, a subjective trial and error approach is performed to set the network. In addition, the retrieval LST/LSE process cannot be well controlled and improving the output is complex due to the difficulty to interpret the network weights.

Synthesis

A review of the different emissivity / temperature separation methods have been made. Generally speaking, the algorithms can be distinguished into three main families:

• Purely mathematics algorithms where a mathematical model aims at simplifying the equation • Derived from physical characteristics of the measurements or the observed object • Derived from a statistical prior on the measurement or the observed object It has been shown that adding more wavelengths to the system does not help to solve the equations' system. The sole introduction of an emissivity polynomial model introduces high correlations and therefore leads to poor results [START_REF] Krapez | Radiative Measurements of Temperature[END_REF]. However, using some prior knowledge about the observed materials is facing the mixed pixel issue: the sensor does not see a unique material but a combination of materials with their own properties. The resulting measurement is a new equivalent material with its own equivalent thermal properties. It is important to note also that some methods that works well in the Remote Sensing field cannot be applied for in-situ infrared thermography measurements. As a conclusion, the Bayesian approaches seem to offer interesting prospects, particularly due to the low priors that can be put on the system. One difficulty of those methods is to adjust well the priors. One drawback also can be the computation time of such methods, particularly when meta-heuristics are involved.

Other methods than the presented one exist. For example, active methods but out-of-scope in this application (see previous chapter). It is also possible to put a reference target with a known emissivity on the field. However, for longterm practical applications of infrastructures, the reference emissivity may drift through time leading to inaccuracies. It may also been difficult to place the target so that the reference measurement is not too much influenced by its surroundings or modifies the dynamic of the measure.

Note

Solving the radiative transfer equation

The equations presented in Chapter 2 can be used as a forward model to simulate the radiative exchanges that occur in a 3D scene. Such simulation tool can be useful for sizing infrared thermography on site before any exploitation and also generate data to test or train algorithms. As a matter of fact, recent buildings constructions come with a building information modeling (BIM) that contains a model of the constructed building, usable for simulation purpose. The presented RTE (Eq. (2.32)) can be solved numerically as in [START_REF] Le | Conception et étude d'infrastructures de transports à énergie positive: de la modélisation thermomécanique à l'optimisation de tels systèmes énergétiques[END_REF] where diffusion, convection and radiative transfers are coupled on a discretized domain in space and angles. However, even if such methodology provides accurate results relatively to its discrete domain, it is also computationally intensive and unfortunately not suitable for fast computation on large refined 3D models. Instead of using the full equation governing the radiation beam in a participating medium Eq. (2.32), the radiosity equation reminded in Eq. (2.48) for the particular case of nonparticipating mediums and diffuse surfaces is used.

J λ (x) = M λ (x) + ρ λ (x) π ∫ S F (x, x ′ )J λ (x ′ )dS(x ′ ) (3.47)
Eq. (3.47) is called a Fredholm integral equation of the second kind.

Note

Finite element approach

Analytical solutions to Eq. (2.48) exist in particular cases only. Instead, a numerical solution can be computed by decomposing the radiosity function on a given basis. Such approximation is called a Galerkine method, on which the domain of study is meshed and the estimated function is interpolated on this mesh. The Galerkine methods are a set of methods that transform a continuous problem to a discrete one and is often used in finite elements methods. Finite element method needs to make the study domain discrete. In the radiosity equation (nonparticipating medium), it leads to cut the domain into a finite number of elementary surfaces (mesh). Let S ⊂ R 3 be the study domain and L 2 (S) the space of square-integrable functions (Lebesgue space):

L 2 (S) = { f ∈ S ∫ S |f (x)| 2 dS(x) < ∞ } (3.48)
That is to say, the function is (almost) finite everywhere in the domain. Let define the inner product

on L 2 (S) ⟨•, •⟩: ∀u, v ∈ L 2 (S), ⟨u, v⟩ = ∫ S uvdS (3.49)
Let now V = H m (S) be an Hilbert space of functions of L 2 (S) with the inner product Eq. (3.49). Let

also ∥•∥ V = ∥•∥ H m .
A weak form can be derived from Eq. (3.48). Let v ∈ V :

∫ S J λ (x)v(x)dS(x) = ∫ S M λ (x)v(x)dS(x) + ∫ S ρ λ (x) π ∫ S F (x, x ′ )J λ (x ′ )v(x)dS(x ′ )dS(x) (3.50)
A bilinear form a(•, •) and a linear form L(•) defined respectively over L 2 (S) × L 2 (S) and L 2 (S) such as:

a(J λ , v) = ⟨J λ , v⟩ V - ∫ S ρ λ (x) π ∫ S F (x, x ′ )J λ (x ′ )v(x)dS(x ′ )dS(x) L(v) = ⟨M λ , v⟩ V (3.51)
Such writing gives the following variationnal form: finding v ∈ L 2 (S) such as:

a(J λ , v) = L(v), ∀v ∈ L 2 (S) (3.52)
In order to use the Lax-Milgram theorem, let show that a(•, •) is continuous on V × V and V -elliptical and that L(•) is a linear continuous form on V. F (x, x ′ ) is square integrable. By using Cauchy-Schwartz and triangular inequality, it follows that:

∀u, v ∈ V, ∥a(u, v)∥ ≤ 2∥u∥ V ∥v∥ V ∀u ∈ V, ∥L(u)∥ ≤ ∥M λ ∥ V ∥u∥ V (3.53)
and therefore a(•, •) is continuous on V ×V and L(•) is continuous on V . From a physical point of view, since the emissivity ranges in the interval ]0; 1[, thus there exists a positive number ϵ 0 such as [START_REF] Mbiock | Radiation in Enclosures: Elliptic Boundary Value Problem, en, Scientific Computation[END_REF]:

∀x ∈ S, ϵ(x) ≥ ϵ 0 > 0 (3.54) Which shows that a(•, •) is V -elliptical: ∀v ∈ V, a(v, v) ≥ ϵ 0 ∥v∥ 2 V (3.55)
As a conclusion, Lax-Milgram theorem ensures that a solution v ∈ L 2 (S) of Eq. (3.52) exists and is unique.

Based on this weak formulation, let V h ⊂ L 2 (S) a finite-dimensional subspace. And let suppose (ϕ i ) 1≤i≤N a basis function of V h . The study domain is now a subset of H 1 (S) made of polynomial functions defined on every element of the mesh. Based on the Galerkine formulation, it is possible to approximate the continuous function J λ to a discrete one Jλ ∈ V h expressed as:

a( Jλ , v) = b(v), ∀v ∈ V h (3.56)
Since Jλ ∈ V h , Jλ can be expressed with the basis functions of V h :

J λ ≈ Jλ = n ∑ j=1 Jλ,j a(ϕ j , ϕ i ) = b(ϕ i ) (3.57)
With Eq. (3.51) and Eq. (3.57):

n ∑ j=1 Jλ,j ( ⟨ϕ j , ϕ i ⟩ V h - ∫ S ρ λ (x) π ∫ S F (x, x ′ )ϕ j (x ′ )ϕ i (x)dS(x ′ )dS(x) ) = ∫ S ϕ i (x)M λ (x)dS(x) (3.58)
For orthonormal basis functions:

⟨ϕ j , ϕ i ⟩ V h = δ ij ∫ S ρ λ (x) π ∫ S F (x, x ′ )ϕ j (x ′ )ϕ i (x)dS(x ′ )dS(x) = ρ i F i→j ∫ S ϕ i (x)M λ (x)dS(x) = M λ,i (3.59)
where F i→j are the form factor defined in the previous chapter, ρ i is the average reflectance of element i and M λ,i the average emission of element i. For an opaque material it follows that:

J λ,i = M λ,i + (1 -ϵ λ,i ) n ∑ j=1 F i→j J λ,j (3.60)
Eq. (3.60) is the discrete version of the radiosity equation. In the following, methods for solving this equation will be reviewed. In particular, even if the origin of this radiosity equation is derived from the thermal transfer theory, it has been widely used in a very different field: graphics rendering [START_REF]Modeling the Interaction of Light Between Diffuse Surfaces[END_REF].

In fact, if radiative heat transfer is interesting in the propagation of thermal radiation through sur-faces, the equations are analogous for light propagation to render diffuses surfaces. Therefore, the following methods are presented from this former field, even though those methods could be applied from a more finite element point of view. In this thesis, the mathematical quantities refer to radiometry and not photometry. The photometry aims at measuring the light in terms of human's perception whereas radiometry measure radiant energy (including light) in terms of absolute power. Therefore, terms radiance, radiant flux and radiant intensity are used in place of luminance, luminous flux and luminous intensity for consistency. Please note that the two fields aim at different goals. When thermal specialists look for an accurate thermal solution of radiative exchanges equations, computer graphics specialists look for a fast solution as real as possible (through the following presented methods). In this thesis, the aim is to develop a simulation tool that gathers those two approaches.

3D image synthesis

Image synthesis history starts during the end of the 1950s for automotive and aeronautic researches purposes. First electronic pens appeared in the first 1960s into an IBM computer so drafting the first mechanical parts through a computer for General Motors was possible [START_REF]IASLonline NetArt: History of Computer Art IV.2 Computer Animation[END_REF]. but with the apparition of home computers, the 1980s are an upsurge time for computer graphics. But it is only in the 1990s that image synthesis and 3D become trendy and now used in many application thanks to the hardware development of graphical processing units (GPU): scientific simulations, CAD, art but also video games. Nowadays, the many algorithms and various applications exist for image synthesis and it is a constant evolving field mixing video games, scientific simulation and animation movies. In this thesis, the objective is to develop a scientific simulation software to render infrared images from in-situ 3D scenes. The light transport equations used in the image synthesis field are similar to the radiation transfer equation. However, one of the main differences between the radiation theory and the 3D image synthesis is the initialization condition. In fact, in most image synthesis algorithms, the emissivity of the surfaces are initialized as follows:

ϵ =    0, if it is not a light object 1, if it is a light object (3.61)
Then, the reflectivity of the surface is used to reflect light. Even if emissive objects are sometimes considered, it represents only a -usually -small part of the scene. At the opposite, in the radiative transfer field all surfaces are supposed to emit which is a non negligible fact. A bibliographic overview of image synthesis is therefore presented in this section. Most recent developed algorithms are now called PBR for physically based rendering (the 2004 reference book has now its 3 rd free online edition [START_REF]Physically Based Rendering: From Theory to Implementation[END_REF])

and tend to respect as most as possible physics' laws to make the rendering as real as possible.

Note

GI algorithms are used in 3D image synthesis and take into account not only the direct illumination of light sources but also the scattered and reflected light form other patches present in the scene.

Therefore, each elementary surface's illumination cannot be computed individually and the system that modelize all the radiative exchanges can only be rendered globally. Many applications are involved such as animation movies (e.g. Disney, Pixar), special effects or video games for texture rendering. Two main types of algorithms exists [START_REF] François | Radiosity and Global Illumination[END_REF]:

• Ray tracing based methods (ray tracing, photon mapping, metropolis light transport, etc.)

• Direct computation of the radiosity's coefficients equation (radiosity method, finite elements, optimization, etc.)

A review of those methods are presented in this section.

Moreover, only a geometrical description of the scene is required, and may not need to be discretized. Due to the independence of the rays, the algorithm can be parallelized. However the main drawbacks are that the rendering equation is not used in basic ray tracing methods and is therefore not considered as photo-realistic since the physic may not be respected. Also, the ray tracing algorithm produces noise due to the random sampling which reduces with the number of particles thrown per pixel (see Fig. 3.7). Finally, the rendered image is view-dependent for standard ray tracing which means that for any movement of the camera in static scenes, the algorithm needs to be computed again. This behavior enables the methods to compute the solution quickly. In fact, the rays are thrown from the pixels of the camera only. View independent ray tracing methods does exist, especially by rendering the scene from multiple views. However, such fix induces a high increase of computation time [START_REF] Dutré | Advanced Global Illumination[END_REF].

Monte-Carlo ray tracing

Monte-Carlo methods are a set of methods that compute a given numerical value based on stochastic processes. The numerical value is approximated by running multiple simulation.

A simple example of this method is to approximate π value. Let C be a circle of R 2 of radius r. Let now S be the square in which C is inscribed, x and y drawn from a uniform distribution in the square. The probability that the point p(x, y) lies into C can be given by the ratio between the area of the square to the area of the circle: The Monte-Carlo ray tracing method is the basic ray tracing that was used to throw a lot of random particles in random directions over the 3D scene and observe their evolution. This method has the advantage of being light and simple to develop and to consider all radiative transfer effects without approximation. However, as shown with the simple π approximation example, the convergence may take a long time and an important particles number to be thrown [START_REF] Shirley | Time Complexity of Monte Carlo Radiosity[END_REF]. As a consequence, this method may lead to an important computational cost.

P (p ∈ C) = A S A C = 4r 2 πr 2 (3.
Path tracing Instead of throwing randomly the rays, other methods were developed to optimize the path taken by the rays. In path tracing, the thrown rays may not have randomness in its process. This is the most common behavior used in most ray tracing algorithms. A ray can be drawn from a deterministic way even if some randomness can be used, particularly for taking into account advanced physical phenomena [START_REF] Dutré | Advanced Global Illumination[END_REF].

Radiosity methods

In order to get a more realistic global illumination algorithm, many developments have been made

to solve the radiosity equation. This equation faces multiple challenges: computing the visibility between the surfaces, computing their form factor and finally solving the equation system. Even if the discrete version of the equation is a linear system, this system can become very large for refined meshes and therefore computationally intensive to solve. On the contrary of ray tracing algorithm, once the radiosity have been computed, the solution is global and does not depend on camera's view.

In real-time 3D computer graphics, radiosity is used to compute the so-called lightmaps. Lightmaps are precomputed textures that aim at caching light data of a 3D scene. Once the radiosity of the scene have been computed, the brightness of the surfaces is stored on such texture. To render the scene, the texture is then mapped against the surfaces.

Note

Progressive radiosity To overcome this issue, Cohen et al. [START_REF] Michael | A Progressive Refinement Approach to Fast Radiosity Image Generation[END_REF] proposed an iterative approach to the solution system. The advantage of this method is its ability to progressively solve the equation system and display the result at each iteration. The original algorithm takes into account diffuses surfaces on nonparticipating media with occlusion. Please note that improvements have been made on this algorithm over time to take into account transmission, absorption, specular reflection, participating media or bidirectional reflection [START_REF] François | A General Two-Pass Method Integrating Specular and Diffuse Reflection[END_REF][START_REF] François | A Global Illumination Solution for General Reflectance Distributions[END_REF][START_REF] Cerezo | A Survey on Participating Media Rendering Techniques[END_REF][START_REF] Michael | Radiosity and Realistic Image Synthesis[END_REF]. Thanks to the technological development of computer's graphical processing units and related software, this algorithm can be optimized to render more efficiently the radiative transfer exchanges inside 3D scenes. The parallelization task is however challeng-ing due to the high dependence that exists between the mesh elements [START_REF] Coombe | Radiosity on Graphics Hardware[END_REF][START_REF] Wallner | An Extended GPU Radiosity Solver: Including Diffuse and Specular Reflectance and Transmission[END_REF][START_REF]AN EMPIRICAL EVALUATION OF A GPU RADIOSITY SOLVER[END_REF]. This algorithm has been retained in this thesis for its physical and view-independent approach. Therefore, more details on this algorithm and implementation will be given on the dedicated Chapter 5.

Improvements on the radiosity method

As stated previously, improvements for computing the radiosity have been done in the literature, in particular to increase the speed of the rendering or to inject more phenomena in the radiosity equation. Some possible improvements are described here.

As for example, general reflectance functions can be approximated by spherical harmonics [START_REF] François | A Global Illumination Solution for General Reflectance Distributions[END_REF].

This method consists in decomposing the radiance into a sum of spherical harmonics, linked to Legendre polynomials, and has been developed as well in the radiative transfer field [START_REF] Modest | Radiative Heat Transfer[END_REF].

Meshing is an important part of the radiosity process which will have an impact on both accuracy and speed computation. In fact, shadows can fall between mesh vertices and can be entirely missed for example. A common extension to the progressive radiosity method is to refine the mesh as the solution is computed [START_REF] Heckbert | Adaptive Radiosity Textures for Bidirectional Ray Tracing[END_REF][START_REF] Vedel | A Testbed for Adaptive Subdivision in Progressive Radiosity[END_REF][START_REF] Paulin | Adaptive Mesh Generation for Progressive Radiosity: A Ray-Tracing Based Algorithm[END_REF] or in post process, which is common in finite element methods [START_REF] Shephard | Approaches to the Automatic Generation and Control of Finite Element Meshes[END_REF].

Hierarchical methods are a meshing strategy used to reduce the computational cost of the radiosity solution. A hierarchical subdivision method is based on the previous adaptive subdivision but this time, the radiosity solution is approximated at different scales [START_REF] Cohen | An Efficient Radiosity Approach for Realistic Image Synthesis[END_REF][START_REF] Michael | Radiosity and Realistic Image Synthesis[END_REF]. The same number of elements is kept but if a local patch needs to be refined due to an important gradient for example, a subdivision can be made on this part only of the solution and a hierarchy is performed which implies to keep only a hierarchical interaction (see Fig. 3.9). This method provides important improvements if the original Figure 3.9: Hierarchical subdivision mesh can be refined without any need of new tesselation (mainly for planes and square shapes). When dealing with concave and convex geometries, the improvements falls drastically: a curved surface already needs to be represented as small polygons and thus, the hierarchical radiosity will not bring a lot of improvements. Nevertheless, some techniques can be added to the radiosity method in the particular case of curved shapes [START_REF]Hierarchical Radiosity On Curved Surfaces[END_REF][START_REF] Cerezo | A Survey on Participating Media Rendering Techniques[END_REF]. This algorithm has been used recently for the simulation of thermal radiative transfer inside houses and implemented on a central processing unit (CPU) [START_REF] Le | Contribution Du Rayonnement Au Confort Thermique et Aux Économies d'énergie Dans l'habitat[END_REF].

When parallelized, the adaptive subdivision becomes more complex; fortunately this can be overcome on current graphic hardwares with the use of quadtree textures [START_REF] Coombe | Chapter 39. Global Illumination Using Progressive Refinement Radiosity[END_REF].

Note

For participating media, the classical radiosity method can be extended through voxels (i.e. pixels in 3D) to consider the radiative exchanges in 3D [START_REF]The Zonal Method for Calculating Light Intensities in the Presence of a Participating Medium[END_REF][START_REF] Arquès | Proximity Radiosity : Exploiting Coherence to Accelerate Form Factor Computations[END_REF][START_REF] François | A Unified Hierarchical Algorithm for Global Illumination with Scattering Volumes and Object Clusters[END_REF].

To complete the previous note on quadtree textures to represent subdivided patches, the voxel representation in recent hardware is made possible by the use of octree [START_REF] Lefebvre | Chapter 37. Octree Textures on the GPU[END_REF]

Note

Finally, please note that other methods also exists. For example meshless methods have been developed [START_REF] Lehtinen | Meshless Finite Elements for Hierarchical Global Illumination[END_REF][START_REF] Lehtinen | A Meshless Hierarchical Representation for Light Transport[END_REF] on which only a mathematical description of the domain is needed, without any need for meshes. In this case, a basis function is used to approximate the solution on the 3D objects at some points. Those points are then interpolated to get the final estimate.

Hybrid methods

To overcome the radiosity difficulty to take into account specular reflection, participating media and also dynamic scenes, hybrid methods that combine both ray tracing and radiosity have also been developed [START_REF] Shirley | Hybrid Radiosity/Monte Carlo Methods[END_REF]. Two common methods called photon mapping and instant radiosity are presented first. Then more recent and advanced ones will be briefly presented.

Photon mapping

The photon mapping approach introduced by Jensen and Chris-tensen [START_REF]Efficient Simulation of Light Transport in Scenes with Participating Media Using Photon Maps[END_REF] is a twopass algorithm. First, the photons are traced from the light sources into the scene. Each time a photon intersects a surface, the intersection point and incoming direction are stored on a cache called photon map. Then, based on this photon map, the rendering can be performed using the radiosity equation. A ray tracing method is used where the ray is divided into categories: direct illumination, specular reflection, transmission and indirect illumination. The direct illumination, specular reflection and transmission are ray-traced but the indirect illumination is computed using the previously cached photon map.

For a given point, N nearest photons are gathered. The photon density is then computed by adding the flux of those photons and divided by the projected area of the sphere that contains those N photons [START_REF] Dutré | Advanced Global Illumination[END_REF]. Photon mapping is a simple and efficient method able to deal with complex geometries and multiple light phenomena such as caustics. However, as stated in [START_REF] Cerezo | A Survey on Participating Media Rendering Techniques[END_REF], it suffers from various artifacts such as blurred shadows and requires a high amount of memory when a lot of light are present on the scene. In fact the photon map is a shortcut that allows the interpolation of photon density. However, if all surfaces are considered to emit such as in the radiation field, then the photon map may not be efficient.

Instant radiosity Instant radiosity is also a two-pass hybrid approach. The aim of the algorithm is to provide a real-time radiosity solution [START_REF] Keller | Instant Radiosity[END_REF]. In this method, the indirect illumination is replaced by the direct illumination from new light source points (called virtual point light (VPL)). As for the photon mapping, a first set of photons are thrown on the scene and new point sources are placed where a trajectory intersects a surface. Then the direct illumination due to these point light sources is accumulated, resulting in an approximation of both direct and indirect illumination. This method provides actually fast results. However, the resulting image may have bias due to the VPL introduced into the scene. First, a too much intensity at the location of the VPL can occur. Then, if the VPL number is too small or not well sampled regarding the geometry, some interactions may be missed. A solution to this issue is given in [START_REF] Segovia | Metropolis Instant Radiosity[END_REF].

Recent developments Based on the VPL, more recent developments have been made, in particular to improve the computation speed of the process and adapt it to modern hardware [START_REF] Walter | Lightcuts: A Scalable Approach to Illumination[END_REF][START_REF] Per | Point-Based Approximate Color Bleeding[END_REF][START_REF] Laurent | Forward Light Cuts: A Scalable Approach to Real-Time Global Illumination[END_REF].

Also, efficient and nowadays often used methods for real-time global illumination algorithms are the light propagation volume method which is an extent of VPL where a hierarchical storing data is performed: every VPL stores which points in the world they light up on a 3D grid. Then the light is spread by using spherical harmonics for the directional reflectance function approximation. This method was developed by Crytek for their rendering engine CryEngine.

Finally, another method is the voxel cone tracing. Again the scene is represented by volume elements (voxels) and stored in an octree texture. Then, instead of tracing a ray represented by a line, the ray is represented by a cone. This cone can then propagates its energy through the octree across the different refinements scales [START_REF] Crassin | Interactive Indirect Illumination Using Voxel Cone Tracing[END_REF].

Those advanced algorithms render a highly looking realistic scene in real-time but still resides in a view-dependent solution and approximation on the radiosity equation. The number of bounces or VPL in those algorithms will refine the solution but increase the computational time. Moreover, those algorithms are also based on the fact that not all the surfaces on the scene are a source of light.

Softwares for radiative simulations

Softwares have been developed in the literature for computing radiative exchanges of 3D scenes.

SOLENE notably, initially developed for architectural purposes, takes into account the radiative exchanges in urban scenes [START_REF] Groleau | Simulation of the Radiative Behavior of an Urban Quarter of Marseille with the SOLENE Model[END_REF][START_REF] Hénon | High Resolution Thermo-Radiative Modeling of an Urban Fragment in Marseilles City Center during the UBL-ESCOMPTE Campaign[END_REF]. It proposes thermal balance, building surface temperature, effect of vegetation on streets etc. This software has multiple modules to make different thermal interactions based on finite element models. The method used for solving the radiative exchanges is the progressive radiosity refinement method, for its rigorous physical approach [START_REF] Lalanne | Development and Validation of a Numerical Tool for Simulating the Surface Temperature Field and the Infrared Radiance Rendering in an Urban Scene[END_REF]. Although the results are interesting and the models used are complete to simulate a micro-climate urban scene, this software is not available freely and suffers from a lack of online documentation.

Synthesis

The radiosity continuous equation have been discretized through a finite element method, inspired from the radiative transfer field. This equation has then been taken to the image synthesis field for rendering realistic 3D scenes. Then, a short overview of different techniques used in image synthesis have been shown. If the ray tracing based algorithms offer a way to render complex phenomena, those methods suffer from their view-point dependence, and the need for denoising the rendered image in post-process. Moreover, some of those methods may not strictly respect the physical equations. At the opposite, methods that aim at solving the radiosity equation are view-independent and respect the physical equations (relative to the spatial discretization). However, radiosity based methods are computationally intensive and are not able to take into account all physical phenomena easily. To conclude, many methods have been developed to simulate the light transport. Some algorithms will therefore perform better in specific situations: from a generic perspective, radiosity based algorithms will behave better in pure diffuse environment and ray tracing in highly specular scenes. In order to simulate the radiative exchange between surfaces for in-situ infrared thermography measurements, the progressive radiosity method has been chosen. It has been implemented on recent graphic hardware and slightly adapted with a GPGPU approach inspired from the image synthesis field and presented in Chapter 5.

Multiples solutions exist to improve the temperature estimation on in-situ infrared thermography.

As shown in this chapter, the first way is to develop emissivity and temperature separation methods to get a more accurate estimate of the temperature, or at least a known uncertainty on this estimation.

This can be achieved through the use of Bayesian methods. Secondly, the optimization of in-situ setup due to the angle and view dependency of the measurement will also have an important impact on the final temperature estimation. As a consequence, a 3D radiative exchanges simulation tool can be used to optimize and visualize the measurement, before any setting up. This tool, can also be used

Chapter 4

In-situ long-term thermal monitoring of structures: environmental measurements bias compensation I n the previous chapters, the different equations that govern the thermal radiative transfer have been exposed. It has been shown that the simplified radiometric equation, used to derive the temperature from infrared measurements depends on many factors. To overcome the lack of knowledge about most of those parameters, one solution is to use other in-situ measurements, or at least find a way to estimate those secondaries variables. Practical instrumentation considerations for in-situ infrared thermal monitoring will be exposed in this chapter. In particular, the need for coupling multiple sensors with an infrared camera will be presented. Handling multi-sensors data over long-term thermal monitoring can lead to an important amount of information to store. Those data need to be well registered. In this case, this is done by using standardized file formats for the records that allow the writing of metadata (a description of the recorded data). Moreover, a specific tool has been developed during this thesis to visualize and process such amount of data.

Like any measuring tool, an infrared camera needs to be calibrated. In order to retrieve the temperature from the digital levels returned by the camera, the process of thermal calibration is shown.

Also, due to the geometrical dependency between the camera and the observed surfaces, a spatial calibration is proposed. Such spatial calibration allows the retrieval of some geometrical parameters such as the angle and distance between the object and the camera.

Moreover, a study of the effects of the different parameters that influence the measurement are also exposed. Finally, it will be shown that it is still possible to estimate some parameters of the model when no coupling has been done. This is achieved by using open data available online.

Infrared camera calibration

This section introduces the thermal and spatial calibration processes, needed to make as trustful measurements as possible. The thermal calibration process is performed prior to the measurements in laboratory conditions. The spatial calibration can be done in different ways, exposed below.

Thermal calibration

Thermal calibration functions

The equations behind the measurements have been presented previously. Once the object radiation L

•(i,j) ∆λ (T (i,j) obj ) has been estimated, it is still necessary to convert this radiation to temperature. The most popular calibration function used for radiation thermometers is the Sakuma-Hattori equation [START_REF] Minkina | Infrared Thermography: Errors and Uncertainties[END_REF][START_REF] Lane | Uncertainty of Temperature Measurements by Infrared Thermography for Metal Cutting Applications[END_REF][START_REF]The Ultimate Infrared Handbook for R&D Professionals[END_REF]. Introduced in 1982 [START_REF] Sakuma | Establishing a practical temperature standard by using a narrow-band radiation thermometer with a silicon detector[END_REF], the Sakuma-Hattori equation provide the electromagnetic signal of thermal radiation from an object's temperature. This equation is then used during the thermal calibration process to relate the blackbody source temperature to the signal received at sensors. Multiple equations have been derived then from this main equation to fit different bandwidth and needs [START_REF] Saunders | General Interpolation Equations for the Calibration of Radiation Thermometers[END_REF]. The Planckian form of the Sakuma-Hattori function is given in the following, where A, B and C are constant values and C 2 is the second radiation constant.

S(T ) = C exp( C 2 AT +B ) -1 (4.1)
where S(T ) is the measured quantity, in digital levels. An equivalent function is often found in the IR thermography field called the RBF-function.

S(T ) = R exp( B T ) + F (4.2)
This equations are mainly used due to their close form to the integrated over spectral responsivity's Planck's law and the fact that it can be inverted.

Thermal calibration process

In order to calibrate thermally the camera, a black body calibration source is used. The infrared camera is placed in a climatic chamber where the temperature is controlled. Different measurements at different source temperature are made. Finally, the RBF function is fitted to the measurements through an optimization algorithm, deriving the R, B and F constants. An example of such calibration process is presented in the next part.

Thermal calibration example

Spatial calibration

The necessity of knowing the geometrical information has been presented earlier. Some methodologies and algorithms to perform camera resectioning (i.e. calibrate spatially the camera) are proposed. Many algorithms and methods have been developed in the computer vision field for many different applications such as machine vision, mobile robot navigation, biomedical, etc. Some differences from the traditional algorithms that we can find in the literature can be noted. In the case of infrared structure monitoring, the input images comes from other spectral ranges than the visible spectrum which leads to a different pixel range dynamics. Moreover, in some cases, it can be difficult to put a geometric reference on a test site for safety purposes. However, we usually use a 3D model of the structure we are monitoring that is assumed to be available.

Camera resectioning

Introduction The objective of camera resectioning is to find the new coordinates of the projected points into the image from the real scene (see Fig. A camera resectioning model can be defined through two types of parameters, represented in Fig. 4.4:

• Extrinsic parameters, that transform the 3D coordinates of the scene into the 3D coordinates of the camera. These parameters mainly depend on the scene and the camera disposition relative to it. • Intrinsic parameters, that transform the 3D coordinates of the camera into pixels coordinates.

The intrisic parameters depend on the camera properties due to the lens used for example.

The simplest model of a perspective matrix is given below. We will review in this section how to extend this model and how to estimate the different parameters.

    f 0 0 0 0 f 0 0 0 0 1 0            x y z 1        =     f x f y z     (4.

Extrinsic parameters

The first step for converting the real world coordinates into pixel coordinates is to translate the information of the camera position.

The camera can be translated and/or rotated from the real world. Therefore, a combination of rotation and translation is enough to define the camera position from the real world:

-→ p c = R w→c -→ p w + -→ T w→c (4.4)
where -→ p c , -→ p w respectively define the point from the camera frame and a point from the scene frame. R w→c is the rotation matrix from the world to the camera and T w→c is the translation one.

This transformation creates 5 degrees of freedom (DOF), due to the need of setting 3 points in R 3 .

• First point p 1 : 3 degrees of freedom (x 1 , y 1 , z 1 )

• Second point p 2 : 1 degree of freedom d = ∥p 1 , p 2 ∥

• Third point p 3 : 1 degree of freedom θ, angle around the line created by p 1 , p 2 .

In order to get from this new camera 3D frame to the pixel frame, we now need to get the intrinsic parameters.

Intrinsic parameters With Eq. (4.3), the simple perspective transformation can be written as:

u = f x z + u 0 v = f y z + v 0 (4.5)
where u 0 and v 0 are the camera origin coordinates in the pixel frame. However, by consediring the fact that pixels are in some arbitrary spatial units, and not always square:

u = α x z + u 0 v = β y z + v 0 (4.6)
Moreover, we can also consider a coefficient of non-orthogonality between x-axis and y-axis of our frame, called "skew" and represented by Fig. 4.5:

u = α x z -α cot(θ) y z + u 0 v = β sin(θ) y z + v 0 (4.7) u v v ′ u ′ θ v ′ sin(θ) = v u ′ = u -cos(θ)v ′ = u -cot(θ)v
    zu zv z     =     α -α cot(θ) u 0 0 0 β sin(θ) v 0 0 0 0 1 0            x y z 1        -→ p ′ =K -→ p
In order to have a more readable and understable matrix, let s be the skew, (c x , c y ) the camera center and a the aspect ratio:

K =     f s c x 0 af c y 0 0 1     (4.8)
Parameters combination By combining Equation 4.4 and Equation 4.8:

-→ p ′ = K(R w→c -→ T w→c ) -→ p = M T -→ p (4.9)
With I 3 the identity matrix, M T is:

M T =     f s c x 0 af c y 0 0 1         1 0 0 0 0 1 0 0 0 0 1 0     ( R 3x3 0 3x1 0 1x3 1
) ( I 3 T 3x1 0 1x3 1 ) (4.

Camera resectioning algorithms

The methods can be divided into two main categories, traditional methods and self-calibration methods [START_REF] Qi | Review on Camera Calibration[END_REF]. We use in the next part the non-linear based method that can take into account the distorsion of the camera. We assume that some reference are known in the image. At least, the Computeraided design (CAD) is known, which is the case for most recent applications of civil engineering.

Classic camera calibration

In traditional methods, a reference on the scene or a set of correspondances points between the scene and the image are known. Then, an optimization algorithm is used to find the camera model parameters. Different methods exists such as direct linear transformation, non-linear methods, two-steps method (Tsai's algorithm), Zhang's method, etc..

Camera auto-calibration

When there is no reference on the scene or that nothing on the scene is known, self-calibration techniques determines the internal camera parameters directly. Such calibration techniques can use multiple views or when the camera is moving perform an active camera selfcalibration.

Direct Linear Transformation

The direct linear transformation (DLT) was introduced in 1974 by I.E Sutherland [START_REF] Sutherland | Three-Dimensional Data Input by Tablet[END_REF] for digitizing 2D drawings to 3D model with a pen. This algorithm solves a set of variables links to a set of similarity equations. Based on this algorithm, [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF] have provided a "Gold standard" algorithm for calibration We have seen in the previous section that at least 6 points were needed to solve the equation system. If we have at least such points, we can create a homogeneous set of equations that defines a least squares problem when over constrained.

With M T = (m ij ) 1≤i≤2 1≤j≤3 and by writting the system for n points, we have:

A 2nx12             m 00 m 01 m 02 ... m 22 m 23             = 0 2n (4.

11)

Where

A =          X 1 Y 1 Z 1 1 0 0 0 0 -u 1 X 1 -u 1 Y 1 -u 1 Z 1 -u 1 0 0 0 0 X 1 Y 1 Z 1 1 -v 1 X 1 -v 1 Y 1 -v 1 Z 1 -v 1 ... .
X n Y n Z n 1 0 0 0 0 -u 1 X n -u 1 Y n -u 1 Z n -u 1 0 0 0 0 X n Y n Z n 1 -v 1 X n -v 1 Y n -v 1 Z n -v 1          (4.12)
Therefore, we can minimize ∥Am∥ to get a least squares problem when over constrained (6 or more points). Since the scale is arbitrary, we can also add the constraint ∥m∥ = 1. This can finally be formulated as a homogeneous least square problem:

min ∥x∥=1 ∥Ax∥ (4.13)
(Get the eigenvector of A T A with the smallest associated eigenvalue).

Singular Value Decomposition

Theorem 1. Let M be a m × n matrix for which coefficients belongs to the field K where K = R or K = C. M can be decomposed such as:

M = U ΣV * (4.14)
where U is a m × m unitary matrix on K, Σ is a m × n matrix for which the diagonal coefficients are reals positives or zero and all others are zeros. V * is the conjugate transpose to V , a n × n unitary matrix on K. This decomposition is called the singular value decomposition of M .

The first way to get the parameters is to solve 4.13 by eigenvalue decomposition. This can be achieved by a singular value decomposition of A. Let y = V T x, therefore we want to minimize ∥Dy∥ subject to ∥x∥ = ∥V T x∥ = ∥y∥ = 1. However, since D is diagonal with decreasing values, the minimum of ∥Dy∥ is found when y = (0, 0, 0, ..., 0, 1) T (4.16)

Therefore, x = V y is the last column in V .

Moreover,

A T A = (U DV T ) T (U DV T ) = V D T U T U DV T = V D T DV T (4.17)
Since D T D is diagonal, this means that the singular value decomposition corresponds to the eigenvalue decomposition. D T D contains the square eigenvalues of A and V contains the eigenvectors. Nonlinear methods Those methods are preferred to the previous one. The idea is to define an error function between the two sets of points (projected 3D points and image positions). Then minimize E with nonlinear optimization methods such as Levenberg-Marquardt. This concept has been introduced by [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF] with the direct linear transformation (DLT) algorithm. The objective is to estimate the transformation matrix M by minimizing an error function f :

min M T f (x i , M T X i ) (4.18)
Prior to this minimization, a normalization is recommended since it provides an invariant behavior regarding the scale and coordinate origins of the input points. Moreover, it also improves the accuracy of the results (see section 4.4.4 of [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]). Different normalization are possible. We will present here the isotropic scaling that we apply to the 3D and 2D dataset of points:

• The points are translated so that their centroid will be at the origin.

• The points are scaled so that in average, they will be around the unit ball of there space in euclidian norm.

Moreover, multiple error functions can be chosen even if they are obviously based on a geometric error function.

Error in one image

The simplest way is to minimize the geometric distance of the projeted points with the objectives points:

∑ i d(x i , M T X i ) 2 (4.19)
2. Error in both images -symmetric transfer error However, in a more realistic case where image measurement errors occur in both the images, we can minimize the error for the forward (M T )

and backward (M -1 T ) transformations:

∑ i d(x i , M T X i ) 2 + d(X i , M -1 T x i ) 2 (4.20)

Error in both images -reprojection error

In this case, we want to get a perfect match of the forward and backward projected points. To do so, we determine new correspondences of our input points x i ↔ xi , X i ↔ Xi so that M will perfectly match the two points and minimize the distance between the correspondances.

∑ i d(x i , xi ) 2 + d(X i , Xi ) 2 , with xi = M T Xi (4.21)
The illustration in Fig. 4.6 inspired from [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF] shows the differences between the symmetric transfer error and the reprojection error.

Levenberg-Marquadt

By using the linear estimate as a starting-point, we minimize the geometric error:

∑ 1≤i≤N d( xi , MT Xi ) Denormalization
The matrix is finally obtained by: M T = V -1 MT U 6. Improving the estimation It has been shown that at least 6 correspondences were needed to solve the system, meaning that the only source of error is in the measurement's point position.

Most of the time this error is said to follow a Gaussian distribution. However, in practical, this assumption may not be verified since mismatch points will be outliers for the Gaussian distribu- Therefore, the profile extraction will induce errors from which an upper bound can be derived from Taylor's series in the case of polynomial interpolation. Of course, those errors are only from a mathematical point of view. Physically, it is even more complicated. In fact, considering the fact that one pixel in the image represents a varying surface in reality make the process of sampling more complex. A square pixel will potentially represent a non-square surface and all the pixels within the image may not represent the same surface in reality. Furthermore it is likely to have the mixed pixel effect, on which the pixel value does not represent a single surface material properties but multiple ones [START_REF] Jones | Scaling of Thermal Images at Different Spatial Resolution: The Mixed Pixel Problem[END_REF] (see also Fig. 2.15). This is actually a challenge often addressed in Remote Sensing [START_REF] Chen | The Mixed Pixel Effect in Land Surface Phenology: A Simulation Study[END_REF].

Recent infrared thermography applications make use of unmanned aerial vehicles (UAV) to inspect large-scale structures [START_REF] Omar | Remote Sensing of Concrete Bridge Decks Using Unmanned Aerial Vehicle Infrared Thermography[END_REF]. The use of UAVs in the field of civil engineering has been investigated years an advanced inspection of the structure and can therefore overcome the former sampling issues. However, those applications are facing other technical difficulties such as motion blur, image stabilization, georeferencing, etc. Furthermore, UAVs may not be suitable for the live monitoring of running systems for security purposes such as bridges with operating cars and most of the employed methods are actives. Nevertheless, IR thermography combined with UAVs is an emerging promising method for many different applications [START_REF] Deane | Application of NDT Thermographic Imaging of Aerospace Structures[END_REF].

Multi-sensor data exploitation

In order to perform a complete thermal monitoring of an infrastructure, one needs to set up a complex instrumentation based on various sensors. The aim of such instrumentation is to overcome the lack of knowledge about environmental and meteorological parameters. A reminder of the simplified radiometric equation to solve is given in Eq. (4.22). Thus, a multi-sensor instrumentation helps for estimating the different involved parameters.

L (i,j) ∆λ,total (T ) =τ opt [ τ atm ( ϵ (i,j) ∆λ,obj L •(i,j) ∆λ,obj (T (i,j) obj ) + (1 -ϵ (i,j) ∆λ,obj )(L • ∆λ,env + L • ∆λ,sun ) +(1 -τ atm )L • ∆λ,atm ] + (1 -τ opt )L • ∆λ,opt (4.22)
Therefore, a need for connected sensors and a flexible architecture with different types of measurements (punctual, in-situ full field thermal images, distributed measurements, etc.) appeared. Even if some tools have already been studied and developed for gathering and recording all the data from such instrumentation [START_REF] Crinière | Cloud2IR: Infrared Thermography and Environmental Sensors Integrated in an Autonomoussystem for Long Term Monitoring of Structures[END_REF][START_REF] Crinière | Management of Local Multi-Sensors Applied to SHM and Long-Term Infrared Monitoring: Cloud2IR Implementation[END_REF], a software for exploiting and helping to process all this amount of data was lacking. In such context, a thermal infrared images post-processing tool has been developed under MATLAB ® in order to visualize infrared images and process data from various sensors. The developed software had to be user-friendly to help for importing data, to navigate between different time stamps of the experiment and also export data to other formats.

Input data and standard formats

The need for interoperability between previously developed tools from the community is primary.

In this sense, the developed software makes use of OGC standards. The Open Geospatial Consortium (OGC) is an international not for profit organization committed to developing open standards. The advantage of using standards is a certain guarantee of the interoperability with other developed tools. Therefore, one of the main feature of the developed software is its ability to read and convert infrared images from To develop such CUDA kernel on MATLAB ® , a C program linked with the CUDA kernel is developed and interfaced with MATLAB ® through a MEX compilation [START_REF]GPU CUDA and MEX Programming -MATLAB & Simulink -MathWorks France[END_REF]. In order to be able to compare the three implementations, the processing time is calculated from the beginning of the conversion call process

Graphical User's Interface (GUI)

To be able to visualize the data efficiently, an interface has been developed to facilitate the processing and be able to detect easily errors or aberration. Fig. 4.12 shows the user's interface.

Batch processing

Even though a GUI is user friendly, it is not efficient to process a large amount of data at once.

Instead, in order to apply some specific processes on a large amount of data, a batch process module have been created. This module relies on YAML Ain't Markup Language (YAML) language [START_REF]The Official YAML Web Site[END_REF]. It is a markup language that has the advantage of being readable for humans, representing data through lists. Furthermore, its non-hierarchical structure but relational based one avoid redundancy: identical data can be referenced multiple times in the same file without the need to repeat the data. The process works as follow, as shown in the schematic in Fig. 4.13:

1. The user describe the input data,

• Folders / files paths,

• Begin date and end date,

• Optionally how to extract the date from the file.

2. The output that will be stored into a global structure with its data type.

The actions to be executed,

• Name of the action.

Definition of the action,

• Name of the action (as in 3.),

• When the action is executed (Begin / For every file or given frame rate / End),

• Function to call for executing the action.

An example of such YAML configuration file is given in Appendix A.

Tests sites and use cases

In order to illustrate the bias through practical use-cases, long-term in-situ measurements made on two mock-up are used: an instrumented wood house and an instrumented road section that are presented in the following. to our multi-sensor data acquisition system Cloud2IR. Due to the particular angle of view, distance and planar properties of the scene, those data are interesting to process to observe spatial sampling issues, mixed pixels and transmissions effects. Moreover, the presence of thermocouples on the road enables us some comparison with infrared measurements.

Instrumented road section

Instrumented wood house

Parameters sensitivity 4.4.1 Sensitivity to emissivity

Different analysis of emissivity sensitivity have been previously done in the literature [START_REF] Chrzanowski | Comparison of Shortwave and Longwave Measuring Thermal-Imaging Systems[END_REF][START_REF] Chrzanowski | Influence of Measurement Conditions and System Parameters on Accuracy of Remote Temperature Measurement with Dualspectral IR Systems[END_REF][START_REF] Pfänder | Pyrometric Temperature Measurements on Solar Thermal High Temperature Receivers[END_REF][START_REF] Krapez | Radiative Measurements of Temperature[END_REF][START_REF] Gaussorgues | Infrared Thermography[END_REF] and [START_REF] Minkina | Infrared Thermography: Errors and Uncertainties[END_REF] especially. Those analyses are extended to the case of the in-situ model presented previously and focus on the emissivity, environment temperature and atmospheric transmission effects. In particular, the use of multi-sensors data for practical use cases to improve the temperature estimation is discussed. By considering that the object's temperature and emissivity are independent:

               ∂L (i,j)
∆λ,total (θ,ϕ) ∂ϵ

(i,j) ∆λ,obj (θ,ϕ) = τ atm τ opt ( L (i,j) ∆λ,obj (θ, ϕ, T (i,j) obj ) -L ∆λ,env -L ∆λ,sun ) ∂L (i,j) ∆λ,total (θ,ϕ) ∂Tenv = ( 1 -ϵ (i,j) ∆λ,obj (θ, ϕ) ) τ atm τ opt ∂L ∆λ,env Tenv = ( 1 -ϵ (i,j) ∆λ,obj (θ, ϕ) ) τ atm τ opt RBe B/Tenv T 2 env (e B/Tenv +C) 2 ∂L (i,j) ∆λ,total (θ,ϕ) ∂τatm = ϵ (i,j) ∆λ,obj (θ, ϕ)τ opt L (i,j) ∆λ,obj (θ, ϕ, T (i,j) obj ) + (1 -ϵ (i,j) ∆λ,obj (θ, ϕ))τ opt ( L ∆λ,env + L ∆λ,sun ) -τ opt L • ∆λ,atm (4.23)
Those equations will be used in the following in order to observe the impact of the different model's parameters to the temperature retrieval.

Emissivity measurements have been made on the wood house (test site #2) in order to get an emissivity map, represented in Fig. 4.16. Figure 4.17: Effect of adjusting the emissivity for different material compared to a constant emissivity However, one can see that adjusting the emissivity will make the temperature estimation more accurate. In fact, in this example, we have chosen the best candidate for a constant emissivity value over all the image because the average emissivity value is 0.93 with a standard deviation of 0.025 (derived from Table 1). As a consequence, we can see that even when the best constant emissivity value is chosen, this is still leading to, at most, 0.5K difference compared to the use of an emissivity map.

Sensitivity to sky temperature and atmospheric transmission

The sky temperature (T sky ) will also have an impact on the final temperature estimation ( Tobj ), as shown in Fig. To overcome this, one can use the measurements made in-situ to get the sky temperature. In Such correlations can be found in the literature [START_REF] Adelard | Sky Temperature Modelisation and Applications in Building Simulation[END_REF]. Fig. 4.20 shows the effect of the Swinbanks correlation for two different periods: in december and august. Moreover, one can see that the sky temperature derived from the correlation will have an impact of the same magnitude's order as if using in-situ measurements.

If local measurements data are not available, the Swinbanks' formula can be used from meteorological air temperature data. The following figure shows the sky temperature derived from air temperature from three different data sources:

• Local weather station

• Copernicus' data store [START_REF]Copernicus Data Store[END_REF] • METAR data Even if differences exists in term of sky temperature (see Fig. 4.21), it represents a small quantity of thermal radiation, regarding the involved temperatures. Therefore, it is always possible to enrich the Moreover, the atmospheric transmission is also influenced by the involved distances. In particular, when dealing with infrastructures monitoring applications, great distances differences may occur due to large field of views (e.g. during the monitoring of a bridge). In such case, a spatial calibration is necessary to compute a distance map and take into account those distances. As for example, the Fig. 4.23

shows the effect of the distance to the atmospheric transmission in the 7.5µm -13µm band and the need for considering the distances due to the camera view projection within the image. As shown in this figure, considering the same distance between the observed surface and the camera may lead to inaccuracies.

Sensitivity to sun contribution

The sun's irradiance will also have an impact on the final temperature estimation. In fact, depending on the emissivity of the observed object, the irradiance at camera's sensor will contain the sun's reflected contribution, both specular and diffuse. An estimation of the solar irradiance can be obtained thanks to a pyrgeometer for example. In the case of the LWIR (7.5µm -13µm), even though the sun contribution may be small compared to the NIR, SWIR and MWIR bands, it stills represents a source of error in the conversion process. In particular for low object's temperatures as illustrated in Fig. 4.24, we can see that the direct sun irradiance contribution can be responsible for up-to 3% of the object's emittance, leading to incorrect temperature estimation if not considered. In the case of SWIR-MWIR band the sun contribution should not be neglected, in particular when observing surfaces with low temperatures (Fig. 4.25).

Other parameters

Finally, rain and fog conditions will have an important impact on the flux received at cameras sensors. Fig. 4.27 illustrates the weather variation during a measurement campaign. Such variation will produce measurements that are difficult to exploit for accurate temperature estimation (see Fig. 4.28).

This may not lead to accurate temperature measurements. However, such effect could be compensated by using a black-body source and perform automatic thermal calibration in case of rain or using adequate radiative transmission model (for instance in case of fog [START_REF] Dumoulin | Infrared Thermography System for Transport Infrastructures Survey with Inline Local Atmospheric Parameter Measurements and Offline Model for Radiation Attenuation Evaluations[END_REF]). 4.24). In this particular case, the final temperature estimation difference with the multi-sensors conversion and without it can be up-to 2K and could be even more in a NIR band due to the sun influence. In this example, the error made on the sky contribution and the emissivity

Synthesis

The long-term infrared thermography measurement process has been reviewed in this chapter:

from the calibration of the sensors to the actual data processing tools. Moreover, the sensitivity of the simplified radiometric equation to its parameters has been studied. It was shown that estimating the parameters through multi-sensors data is efficient. When local sensors data are not available, it is still possible to use online data from various sources (METAR, Copernicus, etc.). After such compensation, one of the most important source of temperature retrieval error is the emissivity. This is why temperature / emissivity retrieval methods are fundamental in the infrared thermal monitoring process. The objective of the last chapter is to propose simultaneous temperature / emissivity separation methods for long-term infrared thermal monitoring. However, finding back the presented in-situ data is a challenging task, particularly because there is no ground truth to trust. As a consequence, an infrared multi-spectral images simulation tool has been developed. This simulation tool is a way to test the developed methods on simple cases. Then, noise can be added by the simulator by considering various environmental parameters. The next chapter will present this simulator based on the progressive radiosity algorithm. Then, developed methods will be presented and finally tested via simulated images.

Chapter 5

Study and development of an infrared multispectral images simulator

A forward solution simulation software has been implemented in order to generate data based on 3D object models. The objectives of this software are to simulate radiative diffuses exchanges between surfaces on a 3D scene as well as integrating atmospheric and infrared camera models for in-situ infrared thermography. The spectral bands used in the software can be adjusted and the computation time has been optimized by using graphical hardware acceleration. Different applications can be done with such simulation tool: simulation's generated data are an efficient and quick way to test temperature and emissivity separation methods robustness or infrared thermal monitoring systems sizing for example. The proposed software gives a forward solution to the diffuse radiative transfer equation through the radiosity method, more detailed in this chapter. The radiosity equations and the form factors computation will be firstly presented. The progressive radiosity method used in the developed software and its implementation on graphic hardware is introduced and discussed.

Then, models to simulate solar spectral irradiance at ground and camera's detectors are conferred to finish with some practical applications examples.

Radiosity method

The radiosity method has been chosen for its physical equations solving approach. The idea is to solve a simplified linear system representing all the contributions of the elements present in the scene.

Such algorithm can be combined with Monte-Carlo methods [START_REF] Shirley | Hybrid Radiosity/Monte Carlo Methods[END_REF] or the equations to be solved can be extended to consider multiple physical phenomena [START_REF] François | A General Two-Pass Method Integrating Specular and Diffuse Reflection[END_REF][START_REF] Wallace | A Two-Pass Solution to the Rendering Equation: A Synthesis of Ray Tracing and Radiosity Methods[END_REF] for considering advanced reflections, participating medium, Fresnel effect, etc. Contrary to Monte-Carlo algorithm, such approach renders the 3D scene globally, id est the generated image does not depend on the viewpoint. Once the scene has been rendered, the camera can be moved inside to inspect each surface. Moreover, Monte-Carlo algorithms suffer from sampling issues and noise. The advantage of Monte-Carlo methods are their ability to consider participating medium, specular reflections and directional reflectance functions easily.

However, in order to get accurate results, multiple bounces are needed which means casting a lot of rays which make the complexity and the time processing more important.

Radiosity equations

The different radiosity methods are all based on the radiosity equation. The aim of the different methods is to solve this linear equation as fast as possible either by providing improved algorithms or by giving an approximation of the solution. The assumptions made here are nonparticipating medium, with orthotropy and lambertian diffuses materials. The form factor F 1→2 which characterizes the geometric interaction between those two diffuse elements is defined in Eq. ( 5.

1) [93]

F 1→2 = ∫ A 1 ∫ A 2 cos(θ 1 ) cos(θ 2 ) πr 2 dA 1 dA 2 (5.1)
For a given spectral band ∆λ i , the radiosity B k,∆λ i of a patch k is given by its self emission E k,∆λ i added to the contribution of all others patches j present on the scene. Therefore, after discretization (see Eq. (3.60)):

J k,∆λ i = M k,∆λ i + (1 -ϵ k,∆λ i ) j=N elements ∑ j=1,j̸ =k V kj F k→j J j,∆λ i (5.2)
where V kj represents an occulting factor

V kj =    1 if k sees j 0 otherwise (5.3)
One can see that three main steps appear from equation 5.3. First, the visibility between the different patches must be solved. Then, the form factors must be computed and finally the system has to be solved. The next part will expose the different methods used for solving those three main steps.

Form factors computation

The direct calculus of the integral of the form factors is computationally intensive, in particular when an important amount of patches is used. Therefore, different methods have been developed

for improving the computation speed of the form factors. Many tables of pre-computed form factors can be found in the literature and used as-is in lookup tables [START_REF] Howell | Thermal Radiation Heat Transfer. 5th[END_REF]. However, such technique cannot be generalized and can only be used on simple cases. In the following, some of the main methods to compute or approximate the form factors are presented.

Contour integration

Instead of being evaluated through the areas, the integral can be evaluated through a contour integration, by applying Stokes' theorem. As a consequence, the multiple integration over a surface area is reduced to a single integration around the boundary of the area [START_REF] Hiram | The Scientific Basis of Illuminating Engineering[END_REF][START_REF] Sparrow | Radiation Heat Transfer[END_REF]. In such context, the form factor from an infinitesimal element to a finite area F d1→2 can be expressed as a sum of single integrations. Such method improves both accuracy and computing time as shown in [START_REF] Shapiro | Computer Implementation, Accuracy, and Timing of Radiation View Factor Algorithms[END_REF].

Nusselt analog

The form factors can also be computed through the unit sphere method introduced by Nusselt [START_REF] Nusselt | Graphische bestimmung des winkelverhaltnisses bei der wärmestrahlung[END_REF].

The idea is to build a hemisphere around the differential element dA 1 and project onto the hemisphere plane the area to be integrated. In such case, the form factor from dA 1 to A 2 becomes:

F d1→2 = 1 π ∫ A 2 cos(θ 1 ) cos(θ 2 )dA 2 r 2 = 1 π ∫ A 2 cos(θ 1 )dΩ (5.4) (dΩ = dA p = cos(θ 2 )dA 2 r 2
). Moreover, cos(θ 1 )dA p is the projection of dA p onto the hemisphere plan (see Fig. 5.2). By integrating the previous quantities, we can derive the form factor such as: The relation Eq. (5.5) can be used to improve the computation speed of the form factors by using the hardware acceleration [START_REF] Coombe | Chapter 39. Global Illumination Using Progressive Refinement Radiosity[END_REF]. In practical, the Z-Buffer hardware algorithm can be exploited to compute the form factors and patches visibility at the same time. However, one drawback of this method is that although polygon edges project to curves on a hemisphere, rasterization produces only straight edges.

F d1→2 = 1 π ∫ A 2 cos(θ 1 )dΩ = A n π (5.
As a consequence, tessellation should be important in the image. More details for such projection are given in the next section.

Hemicube method

One of the most common method used for computing the form factors is the hemicube approach.

Rather than using a unit-sphere, the projection is done on a cube centered over the receiving element [START_REF] Michael | The Hemi-Cube: A Radiosity Solution for Complex Environments[END_REF]. Such method has been developed for its simplicity and efficiency to be implemented on hardware.

However, this method needs to render each face of the hemicube (5 rendering) and depends on the resolution of the chosen hemicube. Small elements that fall into a too large hemicube's resolution will have an important form factor approximation. The Z-Buffer algorithm can also be used for determining the visibility and form factors at the same time. Moreover, in recent hardwares, a special hemicube buffers are available to render hemicubes efficiently. One drawback of the hemicube method is its linear projection that induces issues at edges of the hemicube, particularly for coarse meshing [START_REF] Michael | Radiosity and Realistic Image Synthesis[END_REF].

Single plane method

Different methods have been developed through years for increasing the computation speed such as the planar projection proposed by F. Sillion [START_REF] François | A General Two-Pass Method Integrating Specular and Diffuse Reflection[END_REF]. In this method, a projection plane just above the surface of interest on which other elements are projected is used. This plane is then subdivided adaptively depending on the size of the projected elements to improve sampling. The advantage of this method is to reduce the number of rays to be thrown for the computation of form factors. In particular, in this study, F. Sillion introduces extended form factors for computing specular component.

Monte-Carlo method

Similarly to the ray tracing algorithms for rendering presented in Chapter 3, a ray casting approach is possible to compute the form factors. One of the main advantages of this technique is its ability to handle many surface types, including curved ones. Moreover, it is also possible to get a priori upper bound as described in [START_REF] Pellegrini | Monte Carlo Approximation of Form Factors with Error Bounded a Priori[END_REF]. However, one drawback is that the expense can be more important depending on the distribution to approximate [START_REF] Michael | Radiosity and Realistic Image Synthesis[END_REF]. In Monte-Carlo methods, the full area-to-area form factor integral can be approximated or the Nusselt's analog can be derived [START_REF] Michael | Radiosity and Realistic Image Synthesis[END_REF].

Approximate solution to the linear system

Some iterative methods that could be used for solving the linear radiosity system of equations are presented in this section and translated to the radiosity problem [START_REF] Michael | Radiosity and Realistic Image Synthesis[END_REF]. Most of these methods are available in standard numerical texts books. In this part, the idea is to find the solution to the linear system Ax = b, with A = (a ij ) 1≤i,j≤n , x = (x i ) 1≤i≤n . Then at the end of the proposed method, the algorithm is developed for the radiosity problem.

Relaxation methods

The basic idea behind relaxation methods is to compute iteratively the result of the system Ax = b for large systems.

Fixed point methods

In order to introduce the Jacobi and Gauss-Seidel methods, let decompose A such as A = G -H with G invertible:

Ax = b ⇔ x = G -1 Hx + G -1 b (5.8)
The previous linear equation system is therefore transformed to a problem fixed point search1 :

   x 0 ∈ R n x (k+1) = G -1 Hx (k) + G -1 b
(5.9)

If this sequence converges, then its limit is the solution to the previous linear system Ax = b.

Jacobi method

The Jacobi method has been introduced by C.G.J. Jacobi in the XIX th century. Based on the previous equation, A can be decomposed into two components, a diagonal matrix G (which is invertible provided that its coefficients are also invertibles, assumed here) and a remainder H constituted of the others elements of A: The matrix G -1 H is then called the Jacobi matrix and the system can be solved iteratively:

A = G -H with G =       
x (k+1) i = b i - ∑ n j=1,j̸ =i a ij x (k) j a ii ( 5.11)
A sufficient condition for this method to converge for all x (0) is to have a diagonally dominant matrix. Let now r be the residual vector: r = Axb. In that case, we have:

r (k) i = b i -a ii x (k) i - n ∑ j=1,j̸ =i a ij x (k) j ⇒ x (k+1) i = x (k) i + r (k) i a ii
(5.12)

Where the residuals are updated after each iteration.

Algorithm 6: Pseudo-code for Jacobi radiosity solver

Initialization x (0) i = emission i , r = b -Ax while not converged for i ∈ Faces x (k+1) i = x (k) i - r (k) i a ii r (k+1) = b -Ax (k+1) return x
Gauss-Seidel The Gauss-Seidel method is a slight variation of the Jacobi one. Instead of taking the diagonal matrix as G, the Gauss-Seidel method uses the lower triangular part of the matrix:

A = G -H with G =        a 11 0 . . . 0 a 21 a 22 . . . 0 . . . . . . . . . . . . a n1 0 . . . a nn        and H = -        0 a 12 . . . a 1n 0 0 . . . a 2n . . . . . . . . . . . . 0 0 . . . 0        (5.13)
By using forward substitution, a sequential formula is obtained:

x (k+1) i = b i - ∑ i-1 j=1 a ij x (k+1) j - ∑ n j=i+1 a ij x (k) j
a ii (5.14) The differences between the Gauss-Seidel and the Jacobi methods reside in the amount of elements that are needed to be stored. For a very large problem, Gauss-Seidel is preferred due to the fact that only one storage vector is required due to the fact that elements can be overwritten as they are computed.

Similarly to the Jacobi method, one can express the iteration based on the residuals:

x (k+1) i = x (k) i + r (k) i a ii (5.15)
Where this time the residuals are updated after each step.

Algorithm 7: Pseudo-code for Gauss-Seidel radiosity solver

Initialization x (0) i = emission i while not converged for i ∈ Faces x (k+1) i = x (k) i + r (k) i a ii r (k) i = b i - ∑ n j=1 A ij x (k) j return x
By looking more closely to the Gauss-Seidel algorithm and its physical interpretation in terms of radiosities, one can see that at each step a surface radiosity value i is updated by summing all the contribution of the others surfaces on the scene. Physically speaking, this step is equivalent to take a single patch and consider the radiation from all the surrounding patches which is called a gathering step in the global illumination field, as illustrated in Fig. 5.3. Southwell method A variant of the Gauss-Seidel method is the Southwell algorithm. Developed during the 40', the Southwell algorithm was made for solving linear systems for humans. In this method, instead of relaxing in the i th order, the relaxation is performed on the element with the greatest residual.

Such ordering is called a Southwell iteration [START_REF] Gortler | Radiosity and Relaxation Methods[END_REF][START_REF] Gastinel | Linear Numerical Analysis[END_REF]. However, in order to know which element has the greatest residual at iteration k it seems that all the r i should be computed. To overcome this, a trick is to make use of the changes ∆x (k) between one approximation to the other:

x (k+1) = x (k) + ∆x (k) (5.16)
The residuals are given by:

r (k+1) = b -A ( x (k) + ∆x (k) ) = r (k) -A∆x (k) (5.17)
Since for a given iteration, only ∆x arg max r (k) i ,∀i∈ [1,n] is not zero, we have:

∀j ∈ [1, n], r (k+1) j = r (k) j - a ji a ii r (k) i (5.18)
In such case, the point fixed algorithm is applied to the residuals. The initialization is straightforwardly adapted:

r (0) = b x (0) = 0 (5.19)
The pseudo-code for the Southwell relaxation is given in Algorithm 8

Algorithm 8: Pseudo-code for Southwell radiosity solver

Initialization r i = emission i x i = 0 while not converged i = arg max r j j∈[1,N ] x i = x i + r i a ii temp = r i for j ∈ Faces r j = r j - a ji a ii return x
If each step of the Gauss-Seidel algorithm is O(n), it is more complicated to characterize the complexity of the Southwell algorithm since one given element may be visited multiple times before visiting some other element. The physical comparison with the Gauss-Seidel method can be made. In the Southwell solver, the greatest residual is taken and then the residual of the other surfaces are updated, at each step. On the contrary of Gauss-Seidel, a given element updates the others: it is called a shooting step.

be distributed. At each iteration, the current approximation can be displayed.

Algorithm 9: Pseudo-code for progressive radiosity

Initialization x i = emission i ∆x i = x i while not converged i = arg max (x j Area j ) j∈[1,N ]
//Maximum of energy for j ∈ Faces ∆u i = ∆x i a ji //Unshot energy to be distributed to j ∆x j = ∆x j + ∆u i //j will have to distribute this part of energy x j = x j + ∆u i //Update the approximation of j ∆x i = 0 //i has shot its unshot energy dispImage(x) //Display current approximation return x By looking at the algorithm, one can see that by changing x in the progressive radiosity algorithm by x+r in the Southwell one gives the same algorithm. One difference though is that in the progressive radiosity, the largest energy to shoot is taken instead of the largest residual.

Implementation on accelerated hardware

The progressive refinement approach described previously has been implemented in C++. Parallel radiosity is difficult to have in practice since the computation of the visibility and the form factors in the scene implies many links between the different patches of the scene and therefore to be able to know at a given time the value of other patches. The objective of the software is to provide a fast way to render the radiosities of a 3D scene thanks to recent hardware and software technologies as long as returning physical quantities that can be exploited for infrared thermography. The OpenGL ® 's API [START_REF]OpenGL -The Industry Standard for High Performance Graphics[END_REF] has been used for the graphics rendering, Assimp's library [START_REF]ASSIMP -The Open-Asset-Importer-Lib[END_REF] for loading 3D objects and Nuklear's library for rendering the GUI (Graphical User Interface) elements [START_REF] Mettke | Nuklear: A Single-Header ANSI C Gui Library[END_REF].

Hardware acceleration and OpenGL ®

Figure 5.6: OpenGL ® pipeline. Blue boxes are programmable steps and dashed lines correspond to optional stages [START_REF]Rendering Pipeline Overview -OpenGL Wiki[END_REF].

OpenGL ® (Open Graphics Library) is an API for rendering 2D and 3D vector graphics. This library defines a set of abstract functions that are normalized among the Khronos Group consortium. Graphics processing unit vendors provides an implementation of OpenGL ® (and additional functionalities specifics to each vendors/cards) with their graphic cards so that softwares can interact with the GPU thanks to OpenGL ® . OpenGL ® is used in many types of applications from video games (see Fig. 5.7) to CAD.

Initially, the rendering of any object in 3D was normalized through a rendering pipeline (see Fig. 5.6) that includes automatically and efficiently the tesselation, the clipping, the rasterization, the memory management, the buffers management, etc. Such architecture was designed to be handy, so that the user only had to create 2 type of scripts : one for handling the geometry and manipulating vertices (called vertex shader) and one for handling the rendered image (called the fragment shader). However, with the development of new technologies and new needs for imaging and rendering, the pipeline has been extended so that any part of the pipeline could be customized through scripts. More recently, with the apparition of GPGPU needs for scientist and computing resources, compute shaders appeared with OpenGL ® version 4.3. This time, the compute shaders replace entirely the pipeline and users can perform data manipulation on the GPU with OpenGL ® . This feature has been used in our implementation in order to optimize the computation time of the algorithm.

Progressive radiosity GPGPU implementation

The compute shader enable a GPGPU approach within the OpenGL ® API which appears to be an interesting architecture for the radiosities. The computation of the radiosities is done in a texture through a compute shader. This texture is read for the rendering through a classical OpenGL ® pipeline by using a vertex shader and a fragment shader. This architecture minimizes CPU/GPU transfers. The size of the texture depends on the number of bands and number of patches. The actual texture stores radiosities values as-is and should be improved in the future by using mipmap textures or even hierarchical storage. A flowchart of the developed architecture is presented in Fig. 5.8.

The main steps of the developed software are summarized below:

1. The 3D geometry of the scene is loaded thanks to a 3D file model description 2. Textures are initialized (radiosities values, residuals, areas, etc.) and stored on the GPU memory through 32bits 2D textures 3. A first projection/rendering is done on a hemisphere onto a buffer 4. The highest patch is selected on the GPU 5. The progressive radiosity algorithm is executed with all the previous computed data on a compute shader

The algorithms presented previously are derived from realistic images synthesis that mainly work in the visible spectrum (generally 3 bands). In this particular case, sources (that emits light) must be distinguished from other elements. As a consequence, in the classical approach the radiosities are initialized to a floating-point (usually between 0 and 1), representing the light intensity if it is a source and to 0 otherwise. The values are then kept on [0; 1] interval since only the rendering view aspect really matters. However, in the infrared spectrum any patch will emit radiation and in order to retrieve

M k,∆λ i = ϵ k,∆λ i ∫ ∆λ i L • (λ, T k )dλ (5.21)
Those computed spectral emittances are the initialization of the radiosity and unshot radiosity textures, which corresponds to the first picture of Fig. 5.5.

The initialization step is performed at the beginning, once the 3D model has been loaded. If the meshing is finely refined, the process can be quite long.

The arrays are filled on the CPU and then transfered to the GPU textures. To overcome this, one possible way (not yet implemented) is to parallelize the texture filling. In any case, the textures are saved onto OpenEXR textures files to avoid computing those values again. OpenEXR is a widespread format developed by George Lucas' ILM company for saving 32-bits values textures up-to 64 different channels.

Note

Once the textures have been initialized, some parameters such as atmospherics, solar, or camera parameters and described in the next sections can be tuned. Then, the parallelized progressive radiosity algorithm can be started.

Form factors computation

Instead of using the classical hemicube algorithm, the Nusselt's Analog hemisphere has been used

to reduce the number of renderings and avoid cubemaps usage. In fact, hemicube rendering needs to render the scene from the five sides of the hemicube. Instead, a non-linear projection is performed and used with the Nusselt equations to derive an approximation of the form factors [START_REF] Coombe | Chapter 39. Global Illumination Using Progressive Refinement Radiosity[END_REF]. A square texture is used in offline rendering to render the view from the current's element of interest. Surrounding elements are then projected onto a unit hemisphere by using a non-linear transformation, analogous to a fisheye lens as shown in Fig. 5.9. Even though a non-linear projection is more costly than a linear one, this extra-cost is manageable over 5 renderings.

A comparative study of form factors computation with stereographic projections is done in [START_REF]AN EMPIRICAL EVALUATION OF A GPU RADIOSITY SOLVER[END_REF].

This article points out the fact that using a stereographic projection in the vertex shader accelerates the creation of visibility textures but has many drawbacks. Among those drawbacks, the information of overlapped surfaces is distorted and leads to incorrect shadow boundaries. In fact, since only vertices are affected by the projection, straight edges in the scene are mapped onto straight edges in the projected surface instead of curved segments. This issue becomes critical when the meshes are too coarse.

Extreme deformations may lead to incorrect visibility classification and therefore incorrect distribution of energy. This last drawback leads to inconsistent values and increases the algorithm computation time [START_REF]AN EMPIRICAL EVALUATION OF A GPU RADIOSITY SOLVER[END_REF] due to the fact that the algorithm may not converge correctly if the visibility is not well computed. Instead, a real non-linear projection is shown to overcome most of those issues and to decrease the rendering time considerably. This algorithm has been developed by [START_REF] Gascuel | Fast Non-Linear Projections Using Graphics Hardware[END_REF] and propose actual nonlinear projection for coarse meshes. However, this method has not been implemented in the software.

First, most of the 3D scenes used can be refined thanks to the important amount of available memory on recent hardware, reducing the approximation of elliptic edges to straight ones. The drawback is that the meshing has to be refined everywhere and adaptive subdivision is therefore not possible otherwise distortion will lead to wrong refinement. It seems that adaptive refinement would be possible by using a real non-linear projection (including edges) as described in [START_REF] Gascuel | Fast Non-Linear Projections Using Graphics Hardware[END_REF].

When loading the 3D mesh, each surface is given a color-id. This unique id is then used to retrieved back corresponding elements on the projection texture.

In the actual implementation a limit of 256 3 -2 faces can be used due to the color-id limit (2 colors (black and white) are kept special to check the background for example). This number could be extended by using a hierarchical method or a specific texture containing the ids. In such case, no limit on the number of faces would exist.

Note

The Eq. (5.5) is then used by counting the number of pixels for each element seen on the projection texture. By dividing by the total number of pixels in the hemispherical base, an approximation of the form factors can be obtained. The algorithm is implemented onto the GPU. The OpenGL ® rendering pipeline is used for projecting the scene at every iteration from the point of view of the element that needs to distribute its energy. To keep the texture onto the GPU and avoid back and forths exchanges between the CPU and GPU, the counting of the pixels is performed by using compute shaders, memory barriers and atomic counters [START_REF] Shreiner | OpenGL Programming Guide: The Official Guide to Learning OpenGL[END_REF]. Even though memory barriers break the parallelism essence of the method, it enables the process to be kept onto the GPU without any additional memory transfer.

By taking the equations of the form factors and using the Nusselt analog:

F 1→2 = 1 A 1 ∫ A 1 ∫ A 2 dF d1→d2 dA 1 = 1 A 1 ∫ A 1 F d1→2 dA 1 = 1 A 1 ∫ A 1 A n π dA 1
By assuming that the meshing is refined enough and that the surfaces dk are small enough regarding the involved scene distances; the quantity F d1→2 can be considered constant over A 1 :

F 1→2 ≈ A n π 1 A 1 ∫ A 1 dA 1 = A n π

Note

The radiosity equation is therefore adapted to the implemented Nusselt hemisphere form factor algorithm. Let N b be the number of pixels in the hemisphere basis of the projection texture and el(i) the function that returns the element index based on the pixel color id. Since

A n = N n /N b J k,∆λ i = M k,∆λ i + (1 -ϵ k,∆λ i ) 1 πN b N b ∑ p=1 J el(p),∆λ i (5.

Progressive radiosity implementation

Now that the form factors can be computed, the full algorithm can be implemented and its more detailed pseudo-code is given in Algorithm 10. The distribution of energy is made through a GPGPU approach by using compute shaders. The rendered projection texture is the input of the compute shader which processes each pixel separately on a single grid element. The purpose of this method is to parallelize the For each loop in Algorithm 10 and is applicable thanks to the coherent memory access of OpenGL ® 's Image Load Store functions property.

Next element function

As in the standard progressive radiosity algorithm, the next shooting element is the one with the greatest energy. From a practical point of view, it means that the largest unshot radiosity times area element textures must be found. If finding the maximum value in an array through the CPU is a standard thing it is more difficult on a GPU. In fact, looking into a GPU array the maximum value can be done in parallel but may fail due to concurrent writings and readings.

Concurrency: occurs when two asynchronous processes try to write / read to the same location at the "same" time which can lead to memory corruptions and wrong data written or read. Concurrency

Post-process

Finally, once the radiosity have been computed, the rasterized triangles cannot be directly used for texturing so missing values are linearly interpolated.

Camera model

Camera parameters such as focal length, position, resolution or even acquisition frequency for example are inputs to the software that are used to visualize the 3D scene through this specific camera from one side and also to compute the infrared image from the other part. For instance, the focal length, position and camera resolution are used for the projection of the scene through the pinhole model for the visualization. In the current version, there is no optical transfer function to simulate advanced effects such as diffraction, opening shape and various aberrations but such feature could be added in the future. However, the camera is also defined by a matrix of detectors that can be customized to take into Figure 5.12: Signal received at the camera (left). Generated image by the camera (impulse response), the digitization can be observed as well as the noises impact on the image.

Solar spectral irradiance

In order to consider the in-situ characteristic of the infrared thermography structure monitoring, environmental models have been implemented to simulate the sun contribution, in particular in the NIR band. Such simulation needs the knowledge of the location at a given time to get the solar position in the sky and an atmospheric model to take into account the different absorption, reflection and diffusion effects of the extraterrestrial radiation through the atmosphere. Many radiative transfer codes have been developed through years (a Wikipedia page with a non exhaustive list even exists [START_REF]Atmospheric Radiative Transfer Codes[END_REF]) for various applications such as climate models or solar photo-voltaic systems efficiency optimization for example. In this thesis, the concern was about finding a solar spectral model that could give the solar spectral irradiance at ground and fully documented to be integrated into the software. Those requirements were met by the SMARTS2 (Simple Model for the Atmospheric Radiative Transfer of Sunshine) algorithm described in [START_REF] Gueymard | Simple Model for the Atmospheric Radiative Transfer of Sunshine (SMARTS2) Algorithms and Performance Assessment[END_REF]. This model is derived from SPCTRAL2 and SMARTS1, improved by the comparison and fitting of MODTRAN model for adjusting look-up tables coefficient values [START_REF] Utrillas | A Comparative Study of SPCTRAL2 and SMARTS2 Parameterised Models Based on Spectral Irradiance Measurements at Valencia, Spain[END_REF]. By using look-up tables, the SMARTS2 algorithm provides a fast way to simulate the solar spectral irradiance at ground by considering various atmospheric factors.

Sun position

The sun position can be derived thanks to the knowledge of the longitude (l), latitude (φ), date and time considered [START_REF] Le | Conception et étude d'infrastructures de transports à énergie positive: de la modélisation thermomécanique à l'optimisation de tels systèmes énergétiques[END_REF]. The declination angle δ that represents the angular position of the sun at solar midday from the equatorial plan can be written, from [START_REF] Cooper | The Absorption of Radiation in Solar Stills[END_REF]:

δ = 23.45 sin ( 284 + n 365 × 360 ) (5.23)
Where n is the n th day of the year. Solar hour can be expressed as a function of the legal hour that depends on the time zone UTC+X by:

h sol = 4 60 × (15 × X -l) + E 60 + h leg (5.24)
where the minute correction E that comes from the variation of the Earth self-rotating speed is given by [START_REF] Spencer | Fourier Series Representation of the Position of the Sun[END_REF] (see footnote in [START_REF] Le | Conception et étude d'infrastructures de transports à énergie positive: de la modélisation thermomécanique à l'optimisation de tels systèmes énergétiques[END_REF]).

Therefore, the zenith angle is obtained from [START_REF] Duffie | Solar Engineering of Thermal Processes[END_REF]:

cos(θ z ) = sin(δ) sin(φ) + cos(δ) cos(φ) cos(ω) (5.25)
Finally, the azimuth angle, from the same reference: based on altitude variation. Also, 4 aerosol models (urban, rural, maritime and tropospheric) can be chosen. Fig. 5.14 shows the configuration window in the software with the possibility to plot the resulting transmittance and solar spectrum.

Extraterrestrial solar spectrum

SMARTS2 uses a reference extraterrestrial solar spectrum with 1nm resolution in [280nm; 1700nm] range and 5nm in [1705nm; 4000nm] range. This spectrum was derived from previous studies and adjusted to satellite data. Since this spectrum is interpolated in the current implementation to fit the user's requested resolution, this input spectrum could be changed or adapted to more recent data; if needed.

Direct beam radiation

The direct beam radiation by a surface normal to the sun is obtained by using the following formula:

E bn,λ = E on,λ τ R,λ τ o,λ τ n,λ τ g,λ τ w,λ τ a,λ (5.27)
with

• E on,λ the extraterrestrial irradiance reaching the Earth,

• τ R,λ transmittance due to rayleigh scattering,

• τ o,λ transmittance due to ozone absorption,

• τ n,λ transmittance due to nitrogen dioxide absorption,

• τ g,λ transmittance due to mixed gases absorption (mainly O 2 and CO 2 ),

• τ w,λ transmittance due to water vapor absorption,

• τ a,λ transmittance due to aerosol extinction.

Some examples are provided in the following figures.

Diffuse beam radiation

The diffuse component of the sun spectral radiation is approximated by using the same transmittance functions used in the direct component in order to simplify the computation. A physical justification can be that the photons that are not directly transmitted are scattered in all directions for which a roughly predictable fraction will be directed downwards [START_REF] Gueymard | Simple Model for the Atmospheric Radiative Transfer of Sunshine (SMARTS2) Algorithms and Performance Assessment[END_REF]. The diffuse irradiance is then given by:

E d,λ = E d,R,λ + E a,R,λ + E b,R,λ (5.28)
with • E d,R,λ due to Rayleigh scattering,

• E a,R,λ due to aerosol scattering,

• E b,R,λ due to backscattering (interreflections between the ground and the atmosphere).

Implementation and units tests

The SMARTS2 algorithm and its corresponding look-up tables have been implemented in C++ providing a model to simulate the radiation from the sun. The different parameters can be customized in a specific window as shown in Fig. 5.14. The sun and sky are considered as a hemisphere over the scene and contribute to the radiation through the standard radiosity algorithm for which the patches only emit. The diffuse component is spread over the entire sphere whereas the direct one affects the patches based on the zenith angle. Currently, the circumsolar radiation part of the SMARTS2 algorithm is not implemented yet.

Again, units tests have been implemented to check the reliability between the expected article results and the implemented ones. The results in the article have been reproduced (see Fig. 5.17

Results and examples

The different algorithms and models used to simulate the radiative exchanges on a 3D scene have been exposed and their implementation commented in the previous section. In this section, results are provided and some examples of usage are shown to illustrate the user-friendly part of the developed software and its possible applications.

Graphical User Interface (GUI)

Once the radiosity algorithm has been computed -or during the computation -, the user can switch between multiple rendering shaders for displaying different available results. The radiosites values on every surface on the scene can therefore be observed and more information can be obtained by clicking on a specific surface. Another rendering shader enables the visualization of the energy received at camera's sensor that simulates the irradiance at the camera. This is achieved through Bouguer's law Eq. (2.24) and a model of the detectors of the camera as explained in the next section.

I k,∆λ i = J k,∆λ i cos(θ) d 2 (5. 29
)
Where I is the irradiance received at the camera, d is the distance and θ the angle between the pixel of the camera and the observed surface of radiosity B.

YAML files configuration

The simulation settings can be customized by a YAML (YAML Ain't Markup Language) configuration file. The YAML syntax [START_REF]The Official YAML Web Site[END_REF] allows readable and easy to write syntax. A user can therefore customize the simulation settings and such configuration file is a good way to save and restore different simulation parameters. Variables can be defined locally on the YAML file, according to the YAML 1.2 specification.

Moreover, a specific feature has been added to the YAML parser of the software so that it understands formulas through the ExprTk library [START_REF]C++ Mathematical Expression Library (ExprTk) -By Arash Partow[END_REF] and the tag !!formula. An example and more details about those configuration files are given in Appendix B

Python interpreter and module

Even if the graphical user interface enables interactions with the software, it may be complicated to batch process multiple simulations. Moreover, extending the software or writing C++ code in it to fits the expectations may be difficult for most users. That is why a Python interpreter (CPython) has been included and a specific Python module has been developed so that the user can interface with the software from Python. Python is a widespread high-level programming language among scientist with a strong and dynamic typing. CPython offers a C implementation of Python. Interacting with the software is made easier thanks to the Python interpreter. The user has access to the 3D scene and the radiosity renderer programmatically, which can be useful to perform multiple simulations with different parameters or even animate the 3D scene (see Fig.

Static results

Dynamic results

The simulation enables to have a first sense of the infrared measurement on site, which can be viewed as a decision making tool for infrared measurements systems sizing. Simulating data can help to choose the right measurement system, prior to the exploitation. In the following use-case, this tool has been used in this philosophy to study the feasibility of measuring train's hot boxes [START_REF] Toullier | Comparative Study of Moving Train Hot Boxes Pre-Detection and Axles Counting by in-Situ Implementation of Two Infrared Cameras[END_REF] assuming that those bands correspond approximately to blue, green and red, respectively. Therefore,

Synthesis

An infrared imaging software that simulates radiative exchanges of 3D scenes statics and dynamics have been produced. The implementation of the progressive radiosity method on GPU by using OpenGL ® enables a parallelized rendering and visualization at the same time. Thanks to detectors and solar irradiance spectral data models from the literature, advanced in-situ multi-band infrared thermography simulations are possible. Settings different parameters onto the models through multiple simulations helps in setting infrared measurements systems before the exploitation. Moreover, thanks to the Python module, batch processing and advanced dynamic use-cases can be automatically rendered. However, improvements can be made in the simulation tool. In fact, the parallelized radiosity method could be improved by using an adaptive meshing that makes some refinement based on the energy gradient [START_REF] Coombe | Chapter 39. Global Illumination Using Progressive Refinement Radiosity[END_REF] and therefore modify the projection algorithm [START_REF]AN EMPIRICAL EVALUATION OF A GPU RADIOSITY SOLVER[END_REF][START_REF] Gascuel | Fast Non-Linear Projections Using Graphics Hardware[END_REF]. Moreover, an overshooting method is not performed but could be incorporated to speed-up the convergence [START_REF] Gortler | Radiosity and Relaxation Methods[END_REF]. Also, to overcome the orthotropy, a BRDF (Bidirectional reflectance distribution function) could be used [START_REF] François | Radiosity and Global Illumination[END_REF] and specular reflection as well as transmission could be incorporated as in [START_REF] Wallner | An Extended GPU Radiosity Solver: Including Diffuse and Specular Reflectance and Transmission[END_REF].

The developed software provides a way to study different parameters for in-situ infrared thermography measurements. Its parallel implementation on graphics hardware uses advances from the image synthesis and computer science fields to enables fast interaction and visualization of in-situ 3D spectral scene which has not been done previously according to our knowledge. Moreover, the use of solar and camera's detectors models from the literature and parallelized through the graphics hardware offer a tool to simulate IR in-situăacquisition from the radiative exchanges to the generated image in digitial levels.

A flowchart that summarizes the developed software is given in Fig. 5.25. The software will be used in the following to simulate data for the evaluation of different emissivity and temperature retrieval methods. Proposed and studied methods for the simultaneous estimation of temperature and emissivity

I

n the previous chapter, different emissivity and temperature estimation methods were presented.

Those methods are either restricted to a range of applications and materials or computationally intensive. In particular, no specific method for in-situ IRT structure monitoring was introduced. In this chapter, a presentation of the developed methods for the simultaneous estimation of temperature and emissivity through IRT is made.

Firstly, the state-space representation of the problem is derived from the equations presented in Chapter 3. In order to track the dynamic component of the measured signal (i.e. the temperature under our assumptions), a Kalman filter has been implemented. However, the lack of knowledge about the emissivity and the non-linear characteristics of the model lead us to extend this Kalman filter to a particle filter. Such approach exploits the temporal and multi-bands aspects of the signal but does not take into account the spatial one. Therefore, a Kriging step is finally added to the particle filter which serves two purposes. The first purpose is to integrate into the model a spatial description for the measured signal. The second purpose is to use the Kriging as an interpolator and thus avoid having the particle filter to run on all the pixels of the image but only on a selection of the pixels. This last step reduces the required computation needs.

Results of the different implementations are presented and discussed around a theoretical test case simulation scenario.

Linear time varying (LTV) systems in state-space representation

The methods presented in the following make use of a state-space representation in their formu- where:

• x ∈ R n is the state vector,

• y ∈ R o is the output vector,

• u ∈ R i is the control vector,

• A ∈ M n (R) is the state matrix,

• B ∈ M n,i (R) is the control matrix,

• C ∈ M o,n (R) is the observation matrix,

• D ∈ M o,i (R) is the feedthrough matrix.
with M n (R) the set of matrices with coefficients in R. The first equation is usually called the dynamic equation and the second one the measurement equation.

A discrete version of the previous state-space model can be derived from Eq. (6.1):

x k+1 = A k x k + B k u k y k = C k x k + D k u k (6.2)
Based on the linearization of the temperature-emissivity problem for IRT (Eq. (3.20)), a state-space representation of the emissivity and temperature estimation problem for IRT can be derived. By assuming that the emissivity is constant through time and using the notation defined in Chapter 3, the measurement equation can be defined as:

y k =        1 ∆λ 1 ∆λ 1 2 . . . ∆λ 1 m -C 2 ∆λ 1 1 ∆λ 2 ∆λ 2 2 . . . ∆λ 2 m -C 2
1 ∆λ N ∆λ N 2 . . . ∆λ N m -C 2 ∆λ N                     α 0 α 1 α 2 . . . α m 1 T ′ k              (6.3)
for a given discrete time k. The associated state-space representation can be defined as: The state matrix A k that represents the dynamic evolution of the state vector is undefined, yet. In fact, the evolution of the temperature remains difficult to predict. In the following, different dynamic equations will be considered for characterizing the evolution of the temperature. Even if some models can be found in the literature, a perfect dynamic clear sky ground's temperature evolution may not be suitable for having a good A k estimate. On the other hand, choosing a complex model will lead to add more information onto the system that may not be always available. For the Kalman filter, the dynamic model will be taken as A k = 1, leading to a random walk model. In the particle filtering based methods, the A k will be considered as a parameter to be estimated as well as the emissivity. However, even if it has not been implemented in this work, an auto-regressive based dynamic model is possible. In fact, by using daily temperature measurements, an auto-regressive model can be fitted and then applied to the dynamic equation. Such model is proposed in Appendix D. Finally, such notations lead to this simplified state-space representation:

x k = 1 T k y k = (ln (γ ∆λ i )) 1≤i≤N u k = u =           α 0 α 1 α 2 . . . α m           B k = 0 C k = C =        -C 2 ∆λ 1 -C 2 ∆λ 2 . . . -C 2 ∆λ N        D k = D =        1 ∆λ
x k+1 = A k x k y k = Cx k + Du (6.5)
The C matrix can therefore be defined as in Eq. (3.32) for polynomial approximation or in Eq. (3.41) for spline approximation.

This state-space representation hides many assumptions. The first one is that only the temperature and its evolution are supposed to vary through time (or the emissivity is assumed to vary slowly regarding the variation of temperature). The emissivity and the temperature are assumed to be independent. Those assumptions will be discussed later on. Please note that the presented model will evolve through this chapter to fit the developed methods. In the following, Bayesian filters that will be used to develop the estimation methods are presented.

Optimal Bayes filter

Prediction

The idea is to estimate the next state of the system given past measurements (next and past are related to current time). Let note Y t the measurement up to time t, including the past t -1 measurements. Consider the state-space system (x t , y t ). The objective is to find the prediction p(x t |Y t-1 )

given the belief p(x t-1 |Y t-1 ). By using the law of total probability, we have:

p(x t |Y t-1 ) = ∫ p(x t , x t-1 |Y t-1 )dx t-1 = ∫ p(x t |x t-1 , Y t-1 )p(x t-1 |Y t-1 )dx t-1 = ∫ p(x t |x t-1 )p(x t-1 |Y t-1 )dx t-1 (6.6)
The last line of the equation is obtained by assuming that the phenomena are independent: state x t depends only on state x t-1 . Such assumption means that only the immediate past matters: p(x t |x t-1 , x t-2 , . . . , x 0 ) = p(x t |x t-1 ).

Update

The which appears to be a normalization factor in Eq. (6.7).

Synthesis

Assuming that p(x 0 ) is known, the recursive computation of p(x t , Y t ) is based on the optimal filter Eq. (6.6) and Eq. (6.7): To conclude, the Bayesian filter is based on two main steps:

       p(x t |Y t-1 ) = ∫ p(
• Make prediction based upon previous belief

• Update belief based upon Measurement Three algorithms (Kalman filter, particle filter and interacting Kalman filter) based on those equations are now introduced and then written for the particular case of the conjoint estimation of temperature and emissivity.

Kalman Filter

The Kalman Filter is an infinite impulse response filter that estimates the current state of a linear system based on the previous one and the current measurements. The filter can therefore track linear dynamical models in Gaussian noise. This filter involves two main steps, the prediction and the correction based on the observations.

Equations

White noises are now added to the dynamic and measurement equations and the state-space model is considered as a stochastic process with stochastic variables. The purpose of the Kalman filter is to reconstruct the state variable through time by taking into account this state-space model when the system to observe is assumed to be linear. Initial conditions are assumed to be relatively well-known in order to allow a reliable linearization.

The estimate of the state at time k is denoted by xk|k and xk|k-1 represents the estimate of the system's state at time step k before k -th measurement. Let v k and w k be random independent Gaussian noises white in time (i.e. v i and v j are independent one another for all i ̸ = j, same for w k ). The covariance matrices of xk|k , v k and w k are respectively denoted by Σ k|k , Σ v and Σ w .

Usually normal distributions are taken as a priori. Apart from the fact that this distribution is widely observed in nature, such assumptions comes mainly from its mathematical properties:

• Stable with finite variance

N (x, Σ x ) + N (y, Σ y ) ∼ N (x + y, Σ x + Σ y )
• Convenient to generate isotropic search points

• Maximum entropy relative to all probability distributions covering R but having finite mean and variance [START_REF] Cover | Elements of Information Theory[END_REF].

Note

Prediction

Given the measurements (y i ) 0≤i<k , the expectation xk|k-1 and covariance matrix Σk|k-1 can be written:

xk|k-1 = A k xk-1|k-1 + B k u k Σ k|k-1 = A k Σ k-1|k-1 A T k + Σ v (6.10)

Correction

The Kalman filter alternates two stages: the prediction and the correction of the prediction based on input measurements at step k: y k . The expectation xk|k and covariance matrix Σ k|k are given by using the equations of the best unbiased linear estimator.

ỹk = y k -C k xk|k-1 -D k u k Innovation S k = C k Σ k|k-1 C T k + Σ w Covariance of the innovation K k = Σ k|k-1 C T k S -1 k Kalman gain xk|k = xk|k-1 + K k ỹk Corrected estimation Σ k|k = Σ k|k-1 -K k C k Σ k|k-1 Corrected covariance (6.11)
The Kalman filter aims at correcting the vector's state based on measurements. By taking Eq. (6.9):

• The prediction is obtained by using the previous corrected estimate p(x t-1 |Y t-1 ) and the dynamic model linking the previous state to the next one in time: p(x t |x t-1 ). • The correction at step t is obtained by using the prediction and an observation model p(y t |x t ).

As a consequence, on a Kalman filter, the Gaussian is propagated through time. The mean and standard deviation associated to the dynamic noisy model can be predicted and corrected. The Kalman filter is therefore a feedback system on which the state vector expectation and covariance matrices are updated as new measurements are made. The Kalman filter gives an estimation of x k but also an appreciation of the result through the estimation of the covariance matrix Σ k .

Test case

In order to test this method, emissivity profiles have been taken from the database provided in [START_REF]Determination of Toxicity of Spoil Substrates After Brown Coal Mining Using a Laboratory Reproduction Test with Enchytraeus Crypticus (Oligochaeta) | SpringerLink[END_REF]. The emissivity profile is decimated to get a smaller number of samples. A standard temperature 1 st -order time evolution has been considered (see Fig. Firstly, the model defined in Eq. (3.32), by using a 6-order Tchebychev polynomials basis has been taken. The results are presented in Fig. 6.2. One can see that the polynomial does not fit the emissivity profile. This is mainly due to a local minimum that is found in the process at the beginning and that will not change through the Kalman evolution. However, the temperature evolution does follow the original time profile.

To overcome the issue of the emissivity profile, the Kalman filter is combined with a gradient descent. As a consequence, the emissivity profile is updated at each step, following the temperature evolution as in a Levenberg-Marquardt algorithm. The prediction step in the Kalman is replaced by a gradient descent (with a dumping factor like in the Levenberg-Marquardt algorithm). Let J be the Jacobian matrix of the vectorial function defined by Eq. (3.20) and β the dumping factor. The prediction

J T J + βdiag(J T J ) ) -1 J T (y k -C xk-1|k-1 ) ) + Bu Σ k|k-1 = A k Σ k-1|k-1 A T k + Σ v (6.12)
Let µ ∈ R, the dumping factor is then updated once the correction has been made:

β = β µ , if y k -C xk|k-1 2 < y k-1 -C xk-1|k-1 2 (6.13)
Fig. 6.3 shows the effect of the gradient descent. The first approximations are rough but once more data arrive into the system, the polynom fits more and more to the original profile. Of course, the same input parameters have been used for a fair comparison. As a consequence, the difference between the temperature estimate and the original one is reduced. Even if those results are encouraging, two main issues still remain. The first one is the fact that the emissivity profile is not bounded inside the [0; 1]

interval. The second one is that the 6-order polynomial seems insufficient for fitting the emissivity profile. The first issue can been tackled by using a modified parameterization to fit [START_REF] Masiello | Kalman Filter Physical Retrieval of Surface Emissivity and Temperature from Geostationary Infrared Radiances[END_REF]:

e ∆λ i = ln ( ϵ ∆λ i 1 -ϵ ∆λ i ) (6.14)
The advantage is that e ∆λ i is defined on R while ϵ ∆λ i is kept into [0; 1]. The second part could be tackled with the previous presented splines system. Even though the convergence occurs to a better estimate, the algorithm is still leading to a local solution which is not satisfying even with the splines equation.

Synthesis

One of the main advantages of the Kalman filter is that updates are computationally efficient. However, the Kalman filter is restricted to unimodal distributions with single hypothesis and it can be sensitive to process noise. In the presented Kalman filter, the emissivity, temperature are estimated as state variables of our system which is not accurate regarding Eq. (6.5). In fact, the Kalman applied on the previous defined system leads to a local solution of the actual system. The Kalman filter is not well suited since the previous model is not linear (an extended Kalman filter would also in this particular case converge to a local solution due to a lack of information on the system parameters). Instead, the A k scalar and emissivity should be considered as model parameters and therefore be estimated as is. In the next part, a particle filter in addition to the Kalman filter is presented in order to estimate conjointly the model's parameters and the temperature. Therefore, a particle filter is running in supplement to the Kalman one.

Particle filter

As for the Kalman filter, the particle filter is a Bayes filter that uses a Monte-Carlo process to estimate a probability density function. The density is represented by a set of particles (one possible state) and their associated weights. In the Kalman filter, the prediction is made through the dynamics and uncertainty is added based upon the process noise (dynamic noise). In the particle filter, a prediction distribution is generated by sampling particles based upon their probability, which enables the usage of the Bayes filter in the case that the state-space system is nonlinear. In fact, in that particular case the optimal filter equations cannot be solved explicitly. Then those particles states are propagated through dynamics and perturbed with noise. By tacking back the state system previously presented and formulating it as a non-linear state-system:

x k+1 ∼ f k (x k+1 |x k ) + v k y k ∼ g k (y k |x k ) + w k (6.15) Let B(R n ) be the Borel σ-algebra on R n . The objective is to get p(x t ∈ A) ≈ p(x t |Y t ), where A ∈ B(R n ). If N particles is the number of particles: p(x t ∈ A) = N particles →∞ p(x t |Y t ).
To make those particles evolve through time, a genetic mutation-selection algorithm is performed.

In the particle filter, only the posterior density p(x t |Y t ) is followed by the estimates. In contrast, the estimates obtained with the MCMC method presented in Chapter 3 follows the full posterior density p(x t , x t-1 , . . . , x 0 |Y t ).

Note

Let (ζ i k ) 1≤i≤N particles a set of independent random samples particles with their associated weight w i k , following a probability density p(x 0 ). Two stages are then applied at each iteration: the mutation and the selection. Those two steps are similar to the prediction and update steps performed in the Bayes filter. More details of those equations and convergence results can be found in many textbooks such as [START_REF] Del | Mean Field Simulation for Monte Carlo Integration[END_REF].

Mutation

A random move is locally applied to every particle at the beginning of the iteration according to p(x k |x k-1 ). Such move allows the set of particles to browse the field of solutions. The term

∑N particles i=1 w i k 1 ζ i k
is therefore an approximation of the conditional probability density p(y k |x k ). Once the random move has been done, an approximation of p(y k |x k ) is given by

∑N particles i=1 w i k 1 ζi k .

Selection

The weights of every particle ζi k are evaluated using the Bayes filter equations :

w i k = w i k-1 p(y k |x k = ζi k ) ∑N particles j=1 w j k-1 p(y k |x k = ζj k ) (6.16)

Resampling

As the weights are updated, some particles may be irrelevant due to their low weight. A resampling is therefore done to replace this particle by a new one that will bring more significant information.

Different sampling strategies can be used such as sequential importance resampling (SIR) [START_REF] Gordon | Novel Approach to Nonlinear/Non-Gaussian Bayesian State Estimation[END_REF] or sequential importance sampling (SIS). In order to estimate the moment to perform resampling, the standard effective particles number N eff is used. The resampling is performed when N eff < 0.9N particles , according to [START_REF] Liu | Sequential Monte Carlo Methods for Dynamic Systems[END_REF] N eff = 1

∑N particles i=1 ( w i k ) 2 (6.17)
Resampling consists in sampling new particles based on their weights and removing the ones whose weight is less than 1/N particles . Then all weights are sets to 1/N particles .

Synthesis

The particle filter can be used in the temperature -emissivity estimation by modifying the Eq. (6.5)

and including the emissivity in the state vector [START_REF] Ching | Bayesian State and Parameter Estimation of Uncertain Dynamical Systems[END_REF]. Let xk be this new state vector at time k:

xk = ( x k u T k ) T ũk = 0 Ck = ( D k C k ) Dk = 0 Ãk = ( A k 0 0 I ) (6.18)
In such method, the particles have components that correspond to both parameters and state approximations. Thus, the computed likelihood will take into account the overall approximation, which may not give a good indication about the approximation of the state or the parameters. This is also link with the expression of Ãk for which the parameters are only following a random walk. Instead, another possibility is to use the Kalman filter for the state estimation and combine it with a particle filter for parameters estimation, which is called the interacting Kalman filter (IKF).

Interacting Kalman filter

In this section, a Bayesian parameter recursive estimation based on an interacting Kalman filter is presented. The parameters to be estimated are the emissivity and the scalar A k . Interacting Kalman filter is a combination of two widely used Bayesian estimation methods: the particle filter and the Kalman filter [START_REF] Zghal | Modal parameter estimation using interacting Kalman filter[END_REF]. Please note that in this case, the particle filter is used to estimate the parameters of the model and not the state vector, which is estimated by a Kalman filter [START_REF] Kantas | On Particle Methods for Parameter Estimation in State-Space Models[END_REF].

Test case

A first approximation of f is to take the previous linear function

f k : x k → A k x k with A k a scalar.
Such function could be extended to a more complicated one. As for example, [START_REF] Kawashima | Relations between Surface Temperature and Air Temperature on a Local Scale during Winter Nights[END_REF] shows a regression model of surface temperature depending on wind speed and air temperature, [START_REF]Development of Statistical Temperature Prediction Models for a Test Road in Edmonton[END_REF] presents other models to predict road's surface and depth temperature based on the air temperature. In Appendix D, a seasonal autoregressive integrated moving average model (SARIMA) is presented in order to model the f k function.

Depending on the different sensors available on the test site, the g function could also be extended.

In a theoretical case, the g function could be taken as the material's emittance (used in the following test case). However, in a more realistic study case and as done later on in Section 6.8, it should be integrated over the sensor's measurement band and derived from the simplified radiative transfer model, presented in Eq. (2.51).

A given particle is defined through its parameter vector ζ k :

ζ k = ( a ζ k α ζ 0,k α ζ 1,k . . . α ζ m,k ) T = ( a ζ k α ζ k T ) T (6.19)
Where a ζ k ∈ R is a scalar that describes the temperature evolution through time and

α ζ k = (α ζ i,k = ln(ϵ ζ k (∆λ i ))) T 1≤i≤N (6.20)
Please note that the linearized emissivity model could also be used for which emissivity would be given by ϵ

ζ k (∆λ) = ∑ m i=0 α ζ i,k Φ i (∆λ) for a given (Φ i ) 1≤i≤m function basis (see Eq. (3. 32
)). In this case, the linearized state-space system could still be used. The state-space model for a given particle is then defined as:

   x ζ k+1 = A k (ζ k )x ζ k + v ζ k y ζ k = g(x k , ζ k ) + w ζ k (6.21)
Based on the previous measurement equation, and Eq. (6.20), g can be defined as:

g : (x k , ζ k ) → α ζ k - C 2 ∆λ x k = α ζ k + Cx k (6.22)
The likelihood criterion p(y k |ζ k|k-1 ) is therefore given by [START_REF] Zghal | Modal parameter estimation using interacting Kalman filter[END_REF]:

p(y k |ζ k|k-1 ) = N ( α ζ k + C xζ k|k-1 , CΣ ζ k|k-1 + Σ v ) (6.23)
One can see that the weights are adjusted based on a combination of the Kalman filter and the particle one. In order to improve the result, an additional criterion is added to the particles weights computation, based on the first derivative of the measurement equation in wavelength. By using the linear continuous equation measurement of the system:

∂y(t, λ) ∂λ = 1 ϵ(λ) dϵ(λ)
dλ (6.24)

Synthesis

Fig. 6.5 shows the result on the theoretical test case. One can see that the emissivity profile shape is restored leading to a satisfying temperature estimation. However, it is important to note that the current algorithm implementation is computationally intensive. Although the algorithm can be parallelized with a GPGPU approach, the number of particles needed for a high image resolution would make the algorithm complexity explode. Furthermore, due to the additional wavelength first derivative criterion, the algorithm is suitable for hyper-spectral measurement data only, with short distances between the spectral bands.

y(s) = µ(s) + Z(s) (6.29)
Let s * ∈ S * and S = (s i ) 1≤i≤n , the Simple Kriging equations are [START_REF] Roustant | Re)Dice Packages for Computer Experiments[END_REF]:

µ(s * |S) = µ(s * ) + C(s * , S) T C -1 (y -µ(S)) σ 2 (s * |S) = C(s * , s * ) -C(s * , S) T C -1 C(s * , S) (6.30)
where C = C(s i , s j ) 1≤i,j≤n is the covariance matrix of y(S)

Let consider now that Z is a Gaussian random vector N (0, C) where the covariance function C is known:

[ Y (S) = µ(S) + Z(S) Y (s * ) = µ(s * ) + Z(s * ) ] ∼ N ([ µ(S) µ(s * ) ] , [ C C(s * , s) T C(s * , s) C(s * , s * ) ]) (6.31)
The conditioning on a Gaussian vector leads to

p(Y (s * )|Y (S) = y) ∼ N (µ(s * ) + C(s * , S) T C -1 (y -µ(S)) , C(s * , s * ) -C(s * , S) T C -1 C(s * , S) (6.32)
In the previous equation, we recognize the SK equations. In fact, in the case that Z is assumed Gaussian, the SK equations coincide with the classical conditionning equations. From the SK equations, one can note the following properties of the Kriging: the SK interpolates exactly the experimental data.

The SK variance is independant of the measurements y (homoscedasticity), whatever the chosen kernel C. For a spatially stationary field, the kriging equations are a simple linear minimum mean square error (LMMSE) interpolator.

Universal Kriging

In SK, the trend function µ is assumed to be known. However, in most practical applications, this is not the case. Universal Kriging's (UK) objective is to estimate on the fly the best linear predictions of Y on S and the trend µ. In this case, the trend is of the form µ(s) = ∑ q i=1 β i f i (s) where the vector β = (β i ) 1≤i≤q is unknown and f i : R 2 → R are given basis functions. Therefore, we have:

Y (s) = q ∑ i=1 β i f i (s) + Z(s) (6.33)
Let note F (s * ) = (f i (s * )) 1≤i≤q and abusively, F (S) = (f i (s j )) 1≤i≤q,1≤j≤n . In the case that Z ∼ N (0, C) and β ∼ N (b, B):

[ Y (S) Y (s * ) ] = [ F (S) T β + Z(S) F (s * ) T β + Z(s * ) ] ∼ N ([ F (S) T b F (s * ) T b ] , [ C + F (S)BF (S) T C(S, s * ) + F (S)BF (s * ) T C(S, s * ) T + F (s * )BF (S) T C(s * , s * ) + F (s * )BF (s * ) T ]) (6.34)
By using the conditionning on Eq. (6.34), we get the kriging estimator's equations: In order to simplify the expression of the previous equation, the two following matrix inversion lemma are used [START_REF] Henderson | On Deriving the Inverse of a Sum of Matrices[END_REF]:

µ(s * |S) =F (s * ) T b + ( C(S, s * ) T + F (s * )BF (S) T ) ( C + F (S)BF (S) T ) -1 ( y -F (S) T b ) σ 2 (s * |S) =C(s * , s * ) + F (s * )BF (s * ) T - ( C(S, s * ) T + F (s * )BF (S) T ) ( C + F (S)BF (S) T ) -1 ( C(S, s * ) T + F (s * )BF (S) T ) T (6.
Lemme 6.6.1 ((A + U BV ) -1 = A -1 -A -1 U ( I + BV A -1 U ) -1 BV A -1 ). Lemme 6.6.2 ((A -1 +B -1 ) -1 = A-A(A+B) -1 A).
Let K be a commutative field. Let A and B be two inversible matrices of K N and suppose that their sum and the sum of their inverse are also inversible.

Let x ∈ K N , I be the identity matrix of K N and look for a solution such as

(A -1 + B -1 ) -1 = A + x. Therefore, (A + X)(A -1 + B -1 ) = I ⇐⇒ AA -1 + AB -1 + X(A -1 + B -1 ) = I ⇐⇒ X = -AB -1 (A -1 + B -1 ) -1 ⇐⇒ X = -AB -1 [ (A -1 B + I)B -1 ] -1 ⇐⇒ X = -AB -1 B(A -1 B + I) -1 ⇐⇒ X = -A(A -1 B + I) -1 ⇐⇒ X = -A(A + B) -1 A (6.37)
By using successively the first and second lemma, we can derive the expression of (C+F (s)BF (s) T ) -1 :

( C + F (S)BF (S) T ) -1 = C -1 -C -1 F (S) ( B -1 + F (S) T C -1 F (S) ) -1 F (S) T C -1 = C -1 -C -1 F (S) [ ( F (S) T C -1 F (S) ) -1 - ( F (S) T C -1 F (S) ) -1 ( ( F (S) T C -1 F (S) ) -1 + B ) -1 ( F (S) T C -1 F (S) ) -1 ] F (S) T C -1 (6.
C -1 F (S)) -1 + B) -1 (F (S) T C -1 F (S)) -1 F (S) T C -1 y = F (s * )((F (S) T C -1 F (S)) -1 B -1 + I) -1 (F (S) T C -1 F (S)) -1 F (S) T C -1 y → B→∞ F (s * )(F (S) T C -1 F (S)) -1 F (S) T C -1 y (6.39) Let β = (F (S) T C -1 F (S)) -1 F (S) T C -1 y: F (s * )BF (S) T (C + F (S)BF (S) T ) -1 y → B→∞ F (s * ) β (6.40)
This time, by multiplying the Eq. (6.38) on the left by C (S, s * ) T and on the right by y:

C (S, s * ) T (C + F (S)BF (S) T ) -1 y = C (S, s * ) T [ C -1 -C -1 F (S) ( B -1 + F (S) T C -1 F (S) ) -1 F (S) T C -1 ] y → B→∞ C (S, s * ) T C -1 [ y -F (S) ( F (S) T C -1 F (S) ) -1 F (S) T C -1 y ] = C (S, s * ) T C -1 [ y -F (S) β] (6.41)
Eq. (6.40) and Eq. (6.41) are the two terms in Eq. (6.36):

μ(s * ) = F (s * ) β + C (S, s * ) T C -1 [ y -F (S) β] (6.42) Similarly, σ2 (s * ) = C(s * , s * ) -C(S, s * ) T C -1C(S, s * ) + U T ( F (S) T C -1 F (S) ) -1 U (6.43)
where U = F (S) T C -1 C(S, s * ) -F (s * ) T .

Kriging examples

In order to get an overview of Kriging, here are simple examples made with the DiceKriging [START_REF] Roustant | Two Packages for the Analysis of Computer Experiments by Kriging-Based Metamodeling and Optimization[END_REF] library in R. Those examples are inspired from the DiceKriging help package.

1D interpolation case: covariance kernel type influence

Let 5 measurement points taken every 0.5s starting at 1s. We want to make the interpolation on t ∈ [0, 4] with a second order polynomial trend. It is clear on this simple case that the interpolation will perform differently depending on the chosen covariance. The choice of the covariance kernel will mainly depend on the shape and nature of the input data. As described earlier, the Matérn 5/2 is often used by default, and it is the standard option in the DiceKriging package. Moreover, the confidence interval is larger at the boundaries, due to the lack of information. At the opposite, the measurements points are exactly interpolated.

1D interpolation case: trend influence

Similarly to the choice of the covariance kernel, the choice of the trend function will also have an impact on the prediction. By using the same example but this time with a Matérn 5/2 covariance kernel and different trend models, Fig. 6.8 shows the influence of the trend model, for three different cases: constant, effing and sinusoidal.

2D interpolation case: experimental variogram

In order to illustrate the importance of choosing the right model for the variogram fitting, a simple example of 2D data have been taken from the Meuse R package [START_REF] Pebesma | The Meuse Data Set: A Brief Tutorial for the Gstat R Package[END_REF]. Zinc concentration over a region has to be estimated based on some measurements. Fig. 6.9 shows the differences between a Spherical and a Matérn model for this specific case. The final concentration estimation will therefore be different depending on the chosen covariance model.

Synthesis

The equations of the Kriging have been derived in this section. The simple examples at the end show that the variogram model choice is a fundamental parameter in Kriging estimation. In particular, the fit of the experimental variogram has to be done cautiously in order to be consistent results. In the following, a Kriging step will be combined with the previously presented IKF method. In order to exploit most information from the recorded images, a combination of both Kriging and Kalman methods has been developed. The idea behind this approach is to use the Kriging to conserve the spatial correlation between the objects in the scene and the Kalman filter to use the long-term characteristic of the measurements. Kriged-Kalman Filters (KKF) have been particularly used in the climatic domain, to establish temporal and spatial evolving models with few a-priori [START_REF] Noel | Statistics for Spatio-Temporal Data[END_REF][START_REF] Pierluissi | New Lowtran Band Model for Water Vapor[END_REF]. Such filters can also be used for optimal sensors placement as in [START_REF]Spatio-Temporal Field Estimation Using Kriged Kalman Filter (KKF) with Sparsity-Enforcing Sensor Placement[END_REF] for example.

Main problem statement

Let S ⊂ R 2 be the image space domain and T ⊂ R the time domain. Let suppose that an object has a slow varying emissivity value through time. After removing the geometric coefficients, the measurement equation can be written as:

y(s, t, λ) = ϵ(s, λ)M (s, t, λ), s ∈ S, t ∈ T, λ ∈ R (6.44)
with y the measured value at the sensor and wavelength λ after considering the geometry of the scene, ϵ the emissivity at location s and M the emittance at location s and time t.

A separation of the dynamic and stationary parts of the measure is performed:

log(y(s, t, λ)) = log(ϵ(s, λ)) stationnary + log(M (s, t, λ)) dynamic (6.45)
As done previously by using Wien's approximation and putting a linearized model on ϵ(s, λ), let v s and v d be the stationary component and the dynamic one respectively.

ỹ(s, t, λ) = log(ϵ(s, λ)) - C 2 λT (s, t) = C t (λ) (v s (s) + v d (s, t)) (6.46)
In the following, one important assumption will be that the quantities follow a Gaussian law. This assumption can be discussed regarding the involved log functions in the previous equations. However, since the sampling of the IRT is generally more important than the involved thermal time constants of the observed surfaces, this oversampling is used so that the measurements are averaged over N time values. We finally get the spatio-temporal evolution with the following state model, with a process noise vector ν s (t) ∼ N (0, σ s (t)):

ỹ(s, t, λ) = C t (λ) (v s (s) + v d (s, t)) + ν s (t) (6.47)

Tracking of the dynamic component -Kalman

Let the spatially stationary field v s (s) follows a random process with µ ϵ its spatial expectation and Σ ϵ its spatial covariance matrix:

v s (s) ∼ N (µ ϵ , Σ ϵ ) (6.48)
v s (s) is firstly interpreted as a noise term to get the measurement equation that will be used in the Kalman filter [START_REF] Wikle | A Dimension-Reduction Approach to Space-Time Kalman Filtering[END_REF]:

ỹ(s, t, λ) -C t (λ)µ ϵ = C t (λ)v d (s, t) + C t (λ) (v s (s) -µ ϵ ) + ν s (t) (6.49)
The measurement model equation is therefore derived from Eq. (6.49):

ŷ(s, t, λ) = C t (λ)v d (s, t) + e s (t), with            ŷ(s, t, λ) = ỹ(s, t, λ) -C t (λ)µ ϵ e s (t) = C t (λ) (v s (s) -µ ϵ ) + ν s (t) e s (t) ∼ N (0, C t (λ)Σ ϵ C t (λ) T + σ 2 ν ) (6.50)
With the dynamic equation defined in the previous section and process noise η s (t) ∼ N (0, σ 2 η ), the resulting state-space model is:

   ŷ(s, t, λ) = C t (λ)v d (s, t) + e s (t) v d (s, t) = A t v d (s, t) + η s (t) (6.51)
The Kriging needs therefore an experimental variogram to fit a model to the emissivity component.

In order to get an initial estimate of the emissivity, it is necessary to give a first estimate of the temperature at locations s. This is achieved by choosing a measurement for which the temperature is uniform or at least some knowledge about the temperature gradient is known for the considered image.

Particles update

The weights of the particles are then updated by assuming that locally (on a given neighborhood) the temperature is the same. The likelihood is then derived by reprojecting the estimated emissivity and the given local temperature to a measurement value and computing the mean square error. Contrary to the particle filter presented in section Section 6.4, the particles are thrown for only a given number of pixel over the entire image and not for a single pixel, which reduces the algorithm complexity. The assumption that the emissivity and temperature are constant over this window is made:

∀k ∈ w,    ϵ k = ϵ p T k = T p (6.56)
Therefore, we can formulate the following optimization problem, where Y is the image data and g : (x) → R the observation model :

∀p ∈ (N -2) × (M -2), x p = (ϵ p , T p ) = arg min ⃗ x   ∑ k∈wp ∥Y k -g(x)∥ 2 2   (6.57)
Such optimization problem can be solved by using the standard algorithms such as lsqnonlin function from MATLAB. However, the lsqnonlin does not guarantee a global minimum, but a local one. As a consequence, we compare here the lsqnonlin function with the Trust-Region-Reflective Least Squares algorithm to the CMA-ES metaheuristic.

Trust-Region-Reflective Least Squares algorithm

In order to browse the solution space, the idea is to update the current point to x + h if f (x + h) < f (x). If f (x+h) ≥ f (x) then the region is shrunk and a new h is computed. This algorithm is based on an approximation of the function f to be minimized to a quadratic approximation. Such approximation can be defined with the first two terms of the Taylor approximation of f at x. Then, the problem is formulated as follow:

min ( 1
2 h T Hh + h T g such that∥Dh∥ 2 < ∆ ) (6.58)
where H is the Hessian matrix, g the gradient at x, D is a diagonal scaling matrix and ∆ a positive scalar.

Multi-level approach

In order to solve this system through an optimization problem, a multi-level approach has been used. The algorithm is first applied to a high level (low-resolution) image. Then, the result is used as input to the algorithm on the next lowest level. Therefore, this multi-resolution approach makes the window size defined previously varying automatically. Also, it provides an estimate to the next level.

Of course, the result obtained at low resolution looks blurry, but this result is refined at each iteration, as shown in Fig. 6.11.

Results

Whereas the lsqnonlin can give a local optimum, the CMA-ES has the advantage of providing a global optimum, making it more robust. One other advantage of the CMA-ES is its quasi parameterfree aspect: only the initial solution point, the initial step-size and population size have to be provided.

However such algorithm has a cost and the computation of the result is quite long and expensive, see In this algorithm, the emissivity map makes more sense in terms of physics. In fact, the emissivity has slight changes in this image with the angle of view.

In particular, if we focus on the cuboid at the middle-right part of the scene, we can see that the CMA-ES and lsqnonlin emissivity take into account the angle of view, leading to a more homogeneous surface's temperature, as shown in the Fig. 6.13. An homogeneous temperature surface in this case makes more sense. However, no ground truth information on those data are available, so the final interpretation of this result is difficult. One can still note that the CMA-ES emissivity map is quite blurry, this detail could be improved in a post-processing step.

In order to perform a comparison as fair as possible, it is necessary to compare not only the result but also the processing time. Since the CMA-ES is a stochastic global optimization algorithm, the number of iteration is quite important compared to the trust reflective algorithm.

Synthesis

An overview of different methods for the simultaneous emissivity and temperature estimation have been proposed. In the following section, those methods will be compared to a test case obtained by numerical simulation.

Results analysis

In order to test the developed methods, a study case presented in the next section has been set up.

A comparison of the different introduced methods is performed: MCMC, CMA-ES and KIKF. Due to the fact that the IKF needs small spectral bandwidth, this method has not been tested but will be part of the final results summary proposed at the end of this section.

Study case

The comparison of the methods is done through the creation of a material target (see Fig. The material properties are derived from the ECOSTRESS spectral library [START_REF] Baldridge | The ASTER Spectral Library Version 2.0[END_REF][START_REF] Meerdink | The ECOSTRESS Spectral Library Version 1.0[END_REF] except from the alumina spectrum, derived from previous studies [START_REF] Monchau | Mesure d'émissivité pour la thermographie infrarouge appliquée au diagnostic quantitatif des structures[END_REF]. The four materials have been chosen for their particular emissivity spectral profile as shown in Fig. 6.16. The gold offers a reflective material, while black paint from which the frame is made of, an emissive one. The alumina is particularly interesting for its particular emissive window in the thermal IR band. Wood and concrete materials are commonly used in the civil engineering field. Five different spectral bands are used corresponding to possible infrared detectors spectral bands:

• Band 1: 1.5µm -3µm

• Band 2: 3µm -5µm

• Band 3: 8µm -10µm

• Band 4: 10µm -12µm

• Band 5: 8µm -12µm

In the proposed scenario, the frame temperature remains constant at 293.15K. The temperature of the different materials is assumed to be controlled by a Peltier device and varies according to real temperature measurements at ground taken from the instrumented road section, presented in Chapter 4 and shown in Fig. 6.17. The temperature profile has been decimated to have a sampling rate every 16min.

Monte-Carlo Markov Chain (MCMC)

The Monte-Carlo Markov Chain algorithm, as presented in Chapter 3 has been used on the simulated data. The prior used in this test for the emissivity is weak: N (0.5, 0.5). For the temperature, the initial prior has been chosen with N (273.15, 5) which is reasonable for practical applications where a first guess could be taken by using the air temperature for example. estimation is satisfying on this test case. In fact, one can see that the emissivity spectral shape is retrieved and that the temperature mean tends to the actual temperature. This meta-heuristic provides interesting results on single-time spectral data, with weak priors. However, one drawback of the MCMC method is that it is difficult to know when the convergence has been reached. Moreover, the method is computationally intensive and may need many iterations to have a good estimate. The result has been computed in approximately 4 × 5 = 20minute on a i7-6820HQ CPU @ 2.70GHz on MATLAB ® ), for a single pixel at a single timestamp. As a consequence, only the temperature and emissivity of one pixel at one time for each of the material have been estimated. The implementation could be parallelized through a GPGPU approach to reduce the computation time. Nevertheless, even with a parallel implementation the time needed for a full resolution image over multiple days would take a long time.

Instead, the method could be improved. For instance, a change of the proposal distribution and / or going onto an Ensemble MCMC method or other more advanced methods could be investigated to get better results and faster convergence [START_REF] Drovandi | Ensemble MCMC: Accelerating Pseudo-Marginal MCMC for State Space Models Using the Ensemble Kalman Filter[END_REF]. In this study, a SSIM SSIM m metrics is used instead of the SSIM index, for its mathematical properties [START_REF] Brunet | On the Mathematical Properties of the Structural Similarity Index[END_REF]: A small error on materials with low emissivity leads to a more important error in temperature than materials with higher emissivity. The material average RMSE is up to 4K against 0.05 for the emissivity, which is still lower in average than for a constant emissivity case. Band: 10µm -12µm.

SSIM = √ 1 -SSIM index ( 6
As stated in Chapter 4, a small error on a small emissivity value will lead to an important one on the temperature. Such behaviour is shown with the gold material (top right material). Even if the RMSE and SSIM have low values, the final RMSE and SSIM on the temperature are more important. As a comparison to the MCMC method, the CMA-ES provides interesting results due to its constant local temperature assumption only. No other priors are given. In the current implementation (MATLAB ® ), the full image is processed in approximately 10h. As a consequence, the CMA-ES seems to be a good way to get an emissivity map of the scene without any knowledge. As shown in Chapter 4, an emissivity map can reduce the temperature estimation error drastically. Another possibility for this algorithm is the possibility to identify the different materials within the scene, for classification purposes for example.

Moreover, the algorithm could be improved by using the multi-band measurements, which has not been done yet. Such improvement could lead to disambiguate the solutions and therefore improve the emissivity estimation. In any case, the estimated temperature error is still lower in average than in a constant emissivity case as in Fig. 6.20.

KIKF

Finally, the developed KIKF method has been tested on the study case. However, some adjustments need to be performed in order to have a fair comparison with the other algorithms. First, in the simulation, the materials are assumed to be heated by Peltier modules whereas the frame is at constant temperature of 20 • C. One of the assumption of the KIKF is to suppose that the scene temperature is heated by the sun which does not apply here. As a consequence, the references locations (denoted previously by s) are taken within the materials areas, outside the target's frame. The initial temperature is taken at T 0 = 270K. Then, from this first estimate, the KIKF is run as new data are incoming. During the morning (corresponding to A k < 1), the KIKF filter is applied while a simple Kalman filter is applied during the rest of the day, with the last emissivity approximation. The final estimates are shown in Fig. 6.26. For one given time, the estimation of the emissivity can be more or less accurate. Such behavior could be improved by modifying the particle filter to mutate accordingly to other particles solutions, similarly to particle swarm optimization algorithms. Moreover, it appears that the estimation is more difficult for reflective materials which is an expected result in reality, particularly when dealing with direct solar illumination (not considered here). Fig. 6.27 is an example of one pixel's time evolution. Again, depending on the selected pixel, the estimation will be more or less accurate. In any case, the filter provides information about how the estimates can be trusted: if the estimated temperature is not exact, the error made is at least quantifiable. A small error on materials with low emissivity leads to a more important error in temperature than materials with higher emissivity. The material average RMSE is less than 3K against 0.04 for the emissivity, which is still lower in average than for a constant emissivity case. Band: 10µm -12µm.

This method provides an interesting option to the previous meta-heuristics. In fact, the model is able to track in time the evolution of the system's parameters and could, at term, indicates of a change of properties in the observed structure. Such change can be noticed by the Kriging prediction that will be compromised. A reset factor could be introduced to catch-up any change. Thanks to the Kriging, only part of the pixels are tracked through time, reducing the data storage and computation time. The current implementation in R has estimated the results in 10min for two days of single-band data over the entire image. Moreover, the method is designed for a single-band camera. Introducing the other bands could, as for the CMA-ES disambiguate the solutions. Furthermore, the Kriging model has been put on the pixel spatial covariance. In Chapter 4, the spatial calibration process has been introduced.

This calibration process could be used to get a spatial covariance model based on the world scale and not the image scale anymore which would lead to a more physical meaning. The KIKF method could be extended to take into account the previous remarks but more particularly to improve its robustness.

In fact, the current method needs to set-up the parameters of the Kalman filter, the particle filter and the Kriging which is not always trivial.

Synthesis

In this thesis, methods have been studied and developed to estimate conjointly the temperature and emissivity for in-situ infrared multi-spectral measurements. Bayesian methods offer a manor of including prior information that could be used when additional knowledge is available such as emissivity material specific prior [START_REF] Ash | Temperature-Emissivity Separation for LWIR Sensing Using MCMC[END_REF]. Furthermore, Bayesian methods offer a way to measure the uncertainty made on the estimation which is valuable for practical applications. Table 6.1 shows a summary of the physical assumptions and the application domain of the algorithms. Full image KIKF T 0 at thermal equilibrium Single-band T homogeneous on a given neighborhood Needs temporal data ϵ time evolution slow compared to T Table 6.1: Summary of the compared simultaneous emissivity and temperature methods: physical assumptions If the MCMC and CMA-ES methods provide interesting results, the computation energy and time needed are an obstacle to a practical long-term monitoring solution. At the opposite, the KIKF method proposes a time tracking solution suitable for long-term thermal monitoring of infrastructures. However, this method needs to be deepened to be more robust and easier to setup. Table 6.2 summarizes the method's pros and cons.

Method Physical hypothesis

Chapter 7

Conclusion and future work

The work done in this thesis allowed to review the infrared thermography for in-situ infrastructures long-term thermal monitoring. The necessity and advantages of infrared thermography on a structural health monitoring context has been presented. In order to understand fully the measurement process of infrared thermography, the different equations and quantities used in the radiative transfer theory have been introduced as well. Then, needs for emissivity and temperature separation methods have been emphasized based on those equations.

A bibliographical study has been driven. The temperature and emissivity separation methods from the literature and particularly in the Remote Sensing field have been listed and analysed first. Secondly, with an aim of developing a 3D in-situ infrared scene simulator, the main algorithms of the image synthesis field have been described. Based on the study of literature's algorithms, it has been shown that the progressive radiosity method allows an iterative resolution of the radiative equation system that respect the physical equations.

As a consequence, an infrared simulation tool that uses the progressive radiosity method has been developed. This tool uses the latest computer graphics technologies and has a graphical interface to visualize the 3D scene and the rendering. Based on image synthesis and derived from the thermal field, the progressive radiosity method solves the radiosity equation in a computationally efficient way. Due to the fact that all the surfaces emit in the infrared spectrum, the algorithm has been parallelized and more particularly the shooting step through a GPGPU approach. The software makes use of the graphical libraries and hardware algorithms to render the scene. In addition, a simple model of atmospheric transfer has been transcribed in the software to simulate the spectral solar irradiance at ground. Furthermore, a camera and detector models are used. Please note that the camera and detector implementations are parallelized and rendered on the GPU in real-time by benefiting from hardware acceleration.

However, the detector model is based on quantum detectors and should be modified to other type of detector so that it would actually fit the infrared thermography measurement process. Moreover, the camera model optic does not consider the optical aberrations yet. A transfer function of the camera could therefore be added.

A sensitivity analysis of the parameters that influence the infrared thermography measurement has also been conducted. In particular, the use of coupled measurements refines the radiative model and helps at reducing the uncertainty on the temperature estimate. If for some reason sensors and local data are not available, it is possible to retrieve those informations online thanks to recent open databases. Especially, the use of the European Copernicus climate data store and their API coupled with empirical models from the literature provides a way to estimate most of the parameters of the equation. Moreover, a tool has been developed to visualize and process the large amount of data gathered by long-term measurements. This tool and its graphical interface helps users to get an overview of the spatial, spectral and time data. To improve the display of the images, the radiative conversion process from digital levels to temperatures have been implemented onto the graphical processing unit (GPU). Some extensions could be added to this software. At the present time, the data from online sources are firstly saved onto files and then read and parsed via the software. However, it would be interesting to incorporate the gathering of online data from different sources directly onto the software's interface. Furthermore, a batch processing tool that allows the process of many data at once in a simplified configuration file was developed. This tool is helpful for repeating tasks on different data and could be generalized to any other application.

Finally, new methods for estimating simultaneously the temperature and emissivity have been presented. The Monte-Carlo Markov chains method from the literature has been implemented. By using statistical priors on the temperature and emissivity distributions, the simultaneous estimation of temperature and emissivity is performed, pixel by pixel, at every moment. If this method provides satisfying results, its current implementation suffers from a long processing time. However, improvements could be made by using Ensemble MCMC for example and some parallelization process. A second method based on an interacting Kalman filter has been proposed. It enables the tracking of the temperature through time and the estimation of the emissivity. However, this estimation relies on a hypothesis of small bandwidth which is suitable for hyper-spectral applications only. Then, the CMA-ES optimization method has been used, relying on a local homogeneous hypothesis. As for the MCMC method, if CMA-ES provides interesting results, its computation is intensive. To overcome those issues, a Kriged Interacting Kalman filter (KIKF) has been developed. By using a prior on the spatial emissivity distribution, an interacting kalman filter tracks the temperature and estimates the system parameters. Preliminary results are interesting: on simulation data, the temperature profile is retrieved and the emissivity map values restoration satisfying even if differences with actual values exist. The algorithm needs to be improved to be more reliable, particularly in setting the different filters parameters.

Moreover, tests on laboratory with real measurements data and then on-field measurements need to be performed to fully validate the method.

YAML batch processing

This appendix shows an example of a YAML batch processing file used in the multi-sensors data visualization and processing tool.

Title: Simultaneous characterization of objects temperature and radiative properties through multispectral infrared thermography

Keywords : infrared thermography, thermal monitoring, thermo-optical properties, Bayesian methods, Kalman filters, progressive radiosities Abstract : The latest technological improvements in low-cost infrared cameras have brought new opportunities for long-term infrastructures monitoring. The accurate measurement of surfaces' temperatures is facing the lack of knowledge of radiatives properties of the scene. By using multi-sensors instrumentation, the measurement model can be refined to get a better estimate of the temperature. To overcome a lack of sensors instrumentation, it is shown that online and free available climatic data can be used. Then, Bayesian methods to estimate simultaneously the emissivity and tem-perature have been developed and compared to literature's methods. A radiative exchange simulator of 3D scenes have been developed to compare those different methods on numerical data. This software uses the hardware acceleration as well as a GPGPU approach to reduce the computation time. As a consequence, obtained numerical results emphasized an advanced use of multi-spectral infrared thermography for the monitoring of structures. This simultaneous estimation enables to have an estimate of the temperature by infrared thermography with a known uncertainty.

L

  'importation de scènes 3D issues d'outils de CAO (Conception assistée par ordinateurs) ainsi que la visualisation et le rendu sont faits en temps réel via une interface graphique. L'intégration d'un interpréteur Python au sein du logiciel le rend programmable par l'utilisateur. Ainsi, des scripts Python peuvent être utilisés pour interagir avec la scène 3D et le moteur de rendu. Des scènes statiques ou dynamiques peuvent alors être générées. De plus, l'intégration d'un modèle de rayonnement solaire spectral et de transmission atmosphérique permet une simulation accrue de scènes in-situ. De même, l'introduction d'un premier modèle de caméra infrarouge équipée d'un détecteur quantique dont les performances métrologiques sont paramétrables permet d'optimiser de nouvelles solutions d'instrumentations pour des essais in-situ.

Figure 2 :

 2 Figure 2: Exemple de rendu dynamique avec des boogies de trains, sur deux bandes spectrales différentes. Les données simulées sont comparées à des données terrain.

  (a) Méthode MCMC: estimation de l'émissivité, ici avec un profile correspondant à l'or. (b) Méthode MCMC: histogramme permettant d'obtenir l'estimé en température.

Figure 3 :

 3 Figure 3: Exemple de résultat pour la méthode MCMC pour un pixel à un temps donné. (a) Estimation de l'émissivité. La valeur attendue (continue) est comparée à la valeur estimée (par bande). (b) Estimation de la température. La valeur attendue est représentée en pointillés.

Figure 4 :

 4 Figure 4: Exemple de résultat pour la méthode CMA-ES. Une erreur relativement faible sur une valeur faible d'émissivité peut induire une erreur importante dans l'estimation de la température

Figure 5 :

 5 Figure 5: Estimé final après traitement de l'ensemble des mesures sur une bande spectrale avec la méthode KIKF, pour deux jours de mesures simulées.

Figure 6 :D

 6 Figure 6: Température estimée sur un pixel avec la méthode KIKF: evolution de la température au cours du temps. Comparaison avec la valeur attendue.

  Form factor (or view factor) from element i to element j viiList of Figures

4. 4 75 4. 5 81 4. 7 84 4. 8 131 6. 1 141 6. 4 142 6. 5 146 6. 6 153 xi 6 . 8

 4755817848131114141425146615368 Extrinsics and intrinsics parameters over the camera resectioning. . . . . . . . . . . . . Frame skew . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 4.6 A comparison between the symmetric and reprojection error, as proposed in [85] . . . . Result of interpolation after calibration and ROI recovering . . . . . . . . . . . . . . . . Spatial re-sampling influence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 4.9 HDF file recording format: different data type from various sensors can be registered and then read back fastly due to the HDF internal architecture. . . . . . . . . . . . . . . . 86 4.10 Any text file format can be imported and formated easily . . . . . . . . . . . . . . . . . . 87 4.11 Parallelization of the temperature conversion process onto the GPU. Comparison of the different implementations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 4.12 GUI of the tool . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 4.13 YAML batch processing flowchart . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 4.14 Road section instrumentation: IR image (a) and test site overview (b) . . . . . . . . . . . 90 4.15 Wood house instrumentation: IR image (a) and test site overview (b) . . . . . . . . . . . 91 4.16 Emissivity map for the wood house, colors represents one particular emissivity measurement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 4.17 Effect of adjusting the emissivity for different material compared to a constant emissivity 92 4.18 Difference between true object temperature and estimated one depending on the error made on the environment temperature and emissivity. In this case, the true object is at 293.15K and its emissivity is 0.93. How to read:on the left the true object's temperature and emissivity. A point on the image represents the difference between the true object temperature and the estimated one with Tenv (x-axis) and εobj (y-axis) values. The images represent different environment temperatures values. . . . . . . . . . . . . . . . 93 ix 4.19 Comparison of ground truth temperature values to temperature estimation with and without sky correction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.20 Effect of using a correlation for the sky temperature based on the air temperature . . . . 4.21 Sky temperature derived from different sources of air temperature . . . . . . . . . . . . 4.22 Difference between true object temperature and estimated one depending on the error made on the atmospheric transmission. Here is the difference map for an object temperature of 293.15K and emissivity of 0.93. . . . . . . . . . . . . . . . . . . . . . . . . . 4.23 Different distance values for the atmospheric transmission, through time. . . . . . . . . 4.24 Sun irradiance reflection on 7.5µm-13µm band for a normal object and different emissivity values . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.25 Sun irradiance reflection on 1.5µm -5.4µm band for a normal object and different emissivity values . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.26 Comparison of sun's irradiance measured and simulated on may 5th in the 0.3µm -2.8µm band . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.27 Rain and wind in-situ measurements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.28 IR image view: in dry environment (a), rain with optic issues (b) . . . . . . . . . . . . . . 4.29 A set of environmental data for long-term thermal monitoring . . . . . . . . . . . . . . . 4.30 Comparison between the temperature estimation by considering M tot = ϵ obj M obj or using the presented model and multi-sensor data . . . . . . . . . . . . . . . . . . . . . . 5.1 Form factors for two infinitesimal elements . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Nusselt's analog . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Gathering step as in the Gauss-Seidel algorithm, figure inspired from [40] . . . . . . . . . 5.4 Shooting step as in the Southwell algorithm, figure inspired from [40] . . . . . . . . . . . 5.5 Progressive radiosity examples with a Cornell's Box where all surfaces emit in [1.5µm -5.4µm] band. (a) Initial state. The ceiling source is at 353.15K and the objects at 293.15K with different emissivities. (b) After 10 iterations (≈ 0.15s). (c) After 100 iterations (≈ 0.3s). (d) After 1000 iterations with interpolation to smooth the result (≈ 4.4s). Almost all surfaces have shot their unshot radiosity (red texture). . . . . . . . . . . . . . 5.6 OpenGL ® pipeline. Blue boxes are programmable steps and dashed lines correspond to optional stages [191]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.7 Mirror's Edge Catalyst driven by Unreal Engine 3, which can interface OpenGL ® -Credits to deadendthrills.com . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.8 Principle of the progressive radiosity GPGPU implementation . . . . . . . . . . . . . . . 5.9 Projection image of one of the plane from an element of the other plane . . . . . . . . . 5.10 Concurrency example. At first, the maximum value is 10 at index 0. The concurrency makes the orders in the queues not the same. The maximum value is always guarantee but not its corresponding index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x 5.11 Energy shot at a given time of the algorithm as a percentage of the first shot . . . . . . . 120 5.12 Signal received at the camera (left). Generated image by the camera (impulse response), the digitization can be observed as well as the noises impact on the image. . . . . . . . . 120 5.13 Sun position over the sky vault for a 3D wood house scene . . . . . . . . . . . . . . . . . 122 5.14 Atmospheric model parameters in the software . . . . . . . . . . . . . . . . . . . . . . . 122 5.15 Spectrum of the solar irradiance received by a surface normal to the sun for different visibility and altitude values . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 5.16 Transmittance values for different zenith angles . . . . . . . . . . . . . . . . . . . . . . . 124 5.17 Effect of the pathlengths of ozone and N O 2 to their transmittance and comparison with the original figure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 5.18 Animation of an object in the 3D scene through Python: an object on the scene is rotated of 4 degrees and the scene is displayed on screen at each iteration. . . . . . . . . . . . . 126 5.19 Radiosity map of a building for a given day, at different time (from left to right: 6AM, 12PM and 4PM) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 5.20 3D model of the boogie in its scene environment . . . . . . . . . . . . . . . . . . . . . . . 127 5.21 Boogie and environment in the simulation software. The preview of the rendered infrared image is done in real-time at the bottom left corner. . . . . . . . . . . . . . . . . . 128 5.22 Comparison of simulation data with experimental ones. The motion blur is well represented in the simulation. The temperature levels however are not well represented due to a lack of information about the initial temperatures values of the boogie. . . . . . . . 128 5.23 Original and approximated reflectivity spectrum used for the validation . . . . . . . . . 129 5.24 Cornell's box result for two different lights. Such result is computed on a professional laptop -M1000M Nvidia graphic card, Intel i7-6820HQ CPU @ 2.70GHz. The computation time is approximately 60 seconds with a progressive rendering and aliasing on screen every 50 iterations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 5.25 Flowchart overview of the developed software showing the different steps of the numerical simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Temperature and emissivity profiles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 6.2 6-order Tchebychev approximation through Kalman filtering . . . . . . . . . . . . . . . 140 6.3 6-order Tchebychev polynomial approximation through combined Kalman filter and gradient descent. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Kalman filter flowchart overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Temperature and emissivity estimation through the interacting Kalman filter (8000 particles) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Interacting Kalman filter flowchart overview . . . . . . . . . . . . . . . . . . . . . . . . 147 6.7 Influence of the covariance kernel on the prediction (in black). The 95% confidence interval in light gray and the input measurements points in red. . . . . . . . . . . . . . . Influence of the trend on the prediction (in black). The 95% confidence interval in light gray and the input measurements points in red. . . . . . . . . . . . . . . . . . . . . . . . 6.9 Difference between a spherical and Matérn variogram models for a 2D problem . . . . . 6.10 Kriged Interacting Kalman filter flowchart overview . . . . . . . . . . . . . . . . . . . . 6.11 Multi-resolution approach for solving the optimization problem. From left to right: coarse level to highest final resolution details. . . . . . . . . . . . . . . . . . . . . . . . . 6.12 Comparison of the trust-region-reflective algorithm and the CMA-ES . . . . . . . . . . . 6.13 Cuboid's result comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.14 Comparison of processing time, for one given scale with a given number of pixels . . . . 6.15 Target used for the study case. (a) Rendering of the scene in the visible spectrum. The camera is represented, the target and the hemisphere for the environment. (b) Result example for a frame at 20 • C and other materials at 5.2 • C . . . . . . . . . . . . . . . . . 6.16 Emissivity spectrum of the materials of the target. . . . . . . . . . . . . . . . . . . . . . 6.17 Target material's temperature evolution through time. The measurements correpond to the temperature at ground during two days of january 2017. The decimated temperature profile is shifted by 1K for drawing purposes. . . . . . . . . . . . . . . . . . . . . . . . . 6.18 Focus on the emissivity retrieval on one pixel of each of the four materials for the 8µm -12µm spectral band. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.19 Results of the temperature / emissivity retrieval on one pixel of each of the four materials. (a) (c) (e) (g): Emissivity retrieval. The ground truth profile (continous) is compared to the estimated one (per spectral band).(b) (d) (f) (h): Temperature estimation. The ground truth temperature is represented by the dashed line. . . . . . . . . . . . . . . . . 6.20 Initial prior for the CMA-ES optimization algorithm. A constant emissivity is chosen leading to important differences with the ground truth values. Band: 10µm -12µm. . .

Figure 2 . 1 :

 21 Figure 2.1: Main steps of structural health monitoring process

Figure 2 . 2 :

 22 Figure 2.2: Focus on the IR part of the electromagnetic spectrum. Common spectral bands of IR are shown: Near-infrared (NIR: 0.75 -1.4µm); Short-wavelength infrared (SWIR 1.4 -3µm); Mid-wavelength infrared (MWIR 3 -8µm); Long-wavelength infrared (LWIR 8 -15µm); Far infrared (FIR 15 -1000µm). The atmospheric transmission is also represented.

5 )Fig. 2 .

 52 Fig. 2.3 shows the Planck's law as a function of wavelength, for different temperatures values. The maximum of this distribution is shifted as the temperature changes. This is known as the displacement

Figure 2 . 3 :

 23 Figure 2.3: Planck's law as a function of wavelength at different temperature values

Figure 2 . 4 :

 24 Figure 2.4: Temperature as a function of exitance at different wavelengths.

Fig. 2 .

 2 Fig. 2.4 shows the temperature value as a function of radiant exitance, for different wavelength values. This figure shows two main trends regarding the temperature and exitance dependence. In the MWIR-LWIR bands, highest energy levels correspond to a more important increase of the temperature than in the NIR-SWIR bands.

Figure 2 . 5 :

 25 Figure 2.5: Solid angle on hemisphere. Note that the solid angle is not necessary a circular cone

  2.6). The emission is therefore isotropic and the intensity follows Eq. (2.20) called Lambert's cosine law.

Figure 2 . 6 :

 26 Figure 2.6: Emission rate in a normal direction and off direction with the Lambert's cosine law on the right.

Figure 2 . 7 :

 27 Figure 2.7: Interactions that occur at the interface of an object

Fig. 2 .

 2 [START_REF] Appel | Some Techniques for Shading Machine Renderings of Solids[END_REF] shows some examples of different materials' emissivity profiles.

Figure 2 . 8 :

 28 Figure 2.8: Spectral emissivity of a series of various materials. (Data from reflectance spectra in ASTER spectral library http://speclib.jpl.nasa.gov, Copyright: Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA; see also Balridge, A.M. et al., Remote Sens. Environ., 113, 711, 2009.). The alumina is from previous studies[START_REF] Monchau | Mesure d'émissivité pour la thermographie infrarouge appliquée au diagnostic quantitatif des structures[END_REF].

Figure 2 . 9 :

 29 Figure 2.9: Effect of the state of the surface on emissivity (from [71])

Figure 2 . 10 :

 210 Figure 2.10: Enclosure for which the radiosity is computed at the surface. Point x 3 and x 2 cannot see them each other.

Fig. 2 .Figure 2 . 11 :

 2211 Fig. 2.11 shows the different flux (ϕ emitted, Φ received) that are considered: • ϕ obj , the flux emitted by the object • ϕ env , the flux emitted by the ambient environment and reflected by the object

Figure 2 . 12 :

 212 Figure 2.12: Possible solutions for a given measurement when everything but temperature and emissivity are known

Figure 2 . 13 :

 213 Figure 2.13: The geometry of the scene leads to different angle of view between the object and the camera in the image

Figure 2 . 14 :

 214 Figure 2.14: Different views of the same scene, the spatial sampling is not constant within the images

Figure 2 . 15 :

 215 Figure 2.15: Mixed pixel effect illustration with three different materials.

Figure 2 . 16 :Figure 2 . 17 :

 216217 Figure 2.16: Small uncooled IR camera with 3D-printed case. A coin of 1e is placed at the top as a scale factor.

Figure 2 . 18 :

 218 Figure 2.18: Different types of detectors reading

Figure 2 . 19 :

 219 Figure 2.19: Sensor / Object configuration

Figure 2 .

 2 Figure 2.20: Multi-band measurements bring more information

Figure 3 . 1 :

 31 Figure 3.1: Radiative transfer equation for satellite measurements, done at the top of the atmosphere.

Figure 3 . 2 :

 32 Figure 3.2: ln approximation with its Taylor serie

Figure 3 . 3 :

 33 Figure 3.3: Original irradiance at sensor compared to the irradiance with the estimated emissivity and temperature

Figure 3 . 4 :

 34 Figure 3.4: Found temperature and emissivity profiles

Figure 3 . 5 :

 35 Figure 3.5: Cubic spline interpolation with quadratic boundaries conditions.

Algorithm 1 :Algorithm 2 :

 12 Pseudo code for a Gibbs sampler Function Gibbs(p(T ), p(ϵ|T i-1 , y), p(T |ϵ i , γ), N mcmc , N burn ) ϵ(0..N mcmc ) ← 0 ; T (0..N mcmc ) ← 0 ; Draw T (0) ∼ p(T ); For i ← 1..N mcmc Draw ϵ(i) ∼ p(ϵ|T (i -1), γ); Draw T (i) ∼ p(T |ϵ(i), γ); Remove N burn samples from T and ϵ; However, since the distribution p(T |ϵ, γ) is difficult to compute, a slice-sampling algorithm [72] is used to sample from it. Pseudo code for a slice sampling Function Slice-sampling(f ) Choose x such as f (x) > 0; Repeat once Draw y ∼ U (0, f (x)); Plot an horizontal line at y; Sample the point (x, y) with the edge created at the intersection of f (x) and the horizontal line;

  Further development have been performed during the 1960s in order to develop computer aided design (CAD) programs for 2D and 3D production drawings. The Computer Graphics term appeared in 1961 and is attributed to Boeing art's CAD director William Fetter. In the late 1960s, Utah university specialized in the 3D modeling which aimed at forming a computer model of an object's shape. In beginning of 1970s, Martin Newell at Utah university proposed a 3D model of what will become the famous Utah teapot, which was an object challenging to model but still possible with curves and reflective surfaces. In the late 1970s, 3D computer graphics software appeared for home computers such as 3D Art Graphics in the Apple II. The development of 3D computer graphics were slowdown during the 1970s due to the high hardware cost

  3D image synthesis algorithms can be classified into two categories: local illumination and global illumination. Local illumination methods as their name stipulate consider only the local parameters of the scene (usually the probe lights sources, meshes material properties and normal) to compute the apparent luminosity of the meshes. Those algorithms when executed on material hardware are rendered in real-time really efficiently. However those algorithms (Lamberts, Gouraud, Phong, Cook-Torrance, Kajiya, etc.) are only based on the aesthetic of the rendering and computed empirically. Those methods take into account direct lightning only and therefore exchanges between the elements other than the light → object interactions are not considered. This can even lead to inconsistencies with the first thermodynamic law. At the opposite, global illumination (GI) based algorithms compute both direct and indirect lightning that take into account the incoming radiation (or light) and also the different recursives reflections with the surrounding environment. Those algorithms are mostly based on the radiosity equation presented below.

  (a) Ray-traced Sponza sample scene with Cycles renderer (zoom). (b) Ray-traced Sponza sample scene with Cycles renderer after denoising (zoom).

Figure 3 . 7 :

 37 Figure 3.7: Example of ray-traced rendering with the Sponza scene and Cycles open source renderer.

62 )NoteFigure 3 . 8 :

 6238 Figure 3.8: Monte-Carlo π approximation after 500 iterations and 10000 iterations.

A

  black body calibration source Mikron M310 (ϵ ef f = 1, ∀λ ∈ [8µm, 14µm]) has been used for thermally calibrating the infrared camera Fig. 4.1. Experimental calibration at 10K sampling has been made for T ∈ [293.15K, 343.15K]. The quantity ( T object -T calibration ) 2 is minimized through the Levenberg-Marquardt algorithm. The Fig. 4.2 shows the optimized parameters improvement on the conversion process from the measurements. The distance between the objective temperature and the measure one is reduced and (a) Mikron M310 (b) Infrared view of the source

Figure 4 . 1 :

 41 Figure 4.1: Black-body calibration source

Figure 4 . 2 :

 42 Figure 4.2: Result of thermal calibration

Figure 4 . 3 :

 43 Figure 4.3: Principle of camera resectioning.

Figure 4 . 4 :

 44 Figure 4.4: Extrinsics and intrinsics parameters over the camera resectioning.

Figure 4 . 5 :

 45 Figure 4.5: Frame skew

10 )Extrinsics••

 10 The equation 4.10 is called the pin-hole model. Therefore, a camera matrix is described by several parameters: Translation T from the origin of world coordinates • Rotation R from the world coordinates Intrinsics Focal length f • Principle points (c x , c y ) • Aspect ratio a • Skew s Finally, we have 11 = 5(Projection) + 3(Rotation) + 3(Translation) degrees of freedom.Therefore, we can estimate the transformation matrix if we have a prior knowledge of our observed scene. Which means having at least 6 points for which we know the position in the scene frame and in the image.

Proof.

  Let A = U DV T be the singular decomposition of A (with real coefficients) such as D is sorted by decreasing values. ∥Ax∥ = ∥U DV T x∥ = ∥DV T x∥ (4.15)

Algorithm 3 :∥Ax∥Initialization

 3 Pseudo-code of the SVD computation Objective: Given n ≥ 6 image point correspondences x i ↔ X i , compute the singular value decomposition of A that will minimize the homogeneous least square problem: min ∥x∥=1 Generate matrix A with x i and X i inputsCompute SVD Perform eigenvalues decompositionA T A = V DV TFind the minimum value of D diagonal corresponding to the smallest eigenvalue return the corresponding eigenvector from V Output: Estimate of the calibration matrix M T

Figure 4 . 6 : 4 . 4 Algorithm 4 :

 46444 Figure 4.6: A comparison between the symmetric and reprojection error, as proposed in[START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]

  to be able to recover the region of interest (ROI). However, such operation induces a new sampling obtained by an interpolation and illustrated in Fig.4.7. The reference image in this figure shows the non-square and non-constant spatial sampling in the image.

Figure 4 . 7 :

 47 Figure 4.7: Result of interpolation after calibration and ROI recovering

  (a) The re-sampled profile is dilated and smoothed due to the interpolation (b) The dilation coefficient depends on pixels position

Figure 4 . 8 :

 48 Figure 4.8: Spatial re-sampling influence

  various files formats; from proprietary ones to OGC ones such as HDF. Hierarchical Data Format (HDF) is a digital container format used to structure and handle large amount of data. It is supported by numerous programming languages and offers possibilities to compress data. The data is structured around groups that handles different type of data and datasets that represents a single type of data. Such type of format is necessary for in-situ large scale monitoring. In fact, it allows recording and gathering multiple type of data in one single file (see Fig. 4.9).

  (a) Description of one HDF5 file record (b) View of the recorded file through HDFView software

Figure 4 . 9 :

 49 Figure 4.9: HDF file recording format: different data type from various sensors can be registered and then read back fastly due to the HDF internal architecture.

  to the final converted image (it means that it includes the transfer communication CPU ↔ GPU). Moreover, this amount has been averaged over 60 different images conversions. The standard gpuArray represents an improvement of a factor 10 compared to the CPU conversion whereas the custom CUDA kernel represents an improvement of a factor 40. The difference between the gpuArray and the CUDA kernel can be explained due to the fact that gpuArray offers an easy way for developers to parallelize their algorithm. As a result, it also need to perform multiple checking to be able to handle any user code. Please note that with the built-in functions, developers only need to declare their variables as gpuArray, the code does not change. At the opposite, the CUDA kernel development needs some knowledge on C, CUDA and compilation. Fig.4.11 shows the comparison of those implementations. (a) CUDA kernels in MATLAB ® (b) Time performances for the different methods

Figure 4 . 11 :

 411 Figure 4.11: Parallelization of the temperature conversion process onto the GPU. Comparison of the different implementations

Figure 4 .

 4 Figure 4.12: GUI of the tool

Figure 4 .

 4 Figure 4.13: YAML batch processing flowchart

( a )

 a Image of the road in the infrared spectrum (b) Image of the road in the visible spectrum

Figure 4 . 14 :

 414 Figure 4.14: Road section instrumentation: IR image (a) and test site overview (b)

AFigure 4 . 15 :

 415 Figure 4.15: Wood house instrumentation: IR image (a) and test site overview (b)

Figure 4 . 16 :

 416 Figure 4.16: Emissivity map for the wood house, colors represents one particular emissivity measurement

  4.18 where the difference between the estimated sky temperature ( Tsky ) is compared to the difference made on the estimated object temperature, for different values of T obj and T sky . The R, B and F values used in the model to compute this figure are derived from a thermal calibration, as explained previously[START_REF] Le Touz | Infrared Thermography Applied to the Study of Heated and Solar Pavement: From Numerical Modeling to Small Scale Laboratory Experiments[END_REF]. One can note that those coefficients are only valid for a limited range of object temperature (283.15K -343.15K in this case) and could expect some more important error out of this range. The following figures compare the differences between the actual object temperature and the estimated one by using estimates on T env ( Tenv ) and ϵ obj (ε obj ).

Figure 4 . 18 :

 418 Figure 4.18: Difference between true object temperature and estimated one depending on the error made on the environment temperature and emissivity. In this case, the true object is at 293.15K and its emissivity is 0.93. How to read: on the left the true object's temperature and emissivity. A point on the image represents the difference between the true object temperature and the estimated one with Tenv (x-axis) and εobj (y-axis) values. The images represent different environment temperatures values.

Fig. 4 .

 4 Fig.4.19, the ground truth temperature obtained from thermocouples data is compared to the temperatures estimated with and without sky correction. In the case of the sky correction, the sky temperature is derived from a pyrgeometer data whereas a constant sky temperature is taken in the other case. The Fig.4.[START_REF] Berk | MODTRAN: A Moderate Resolution Model for LOWTRAN[END_REF] shows that taking into account the sky temperature from local measurements gives a better estimate of the object temperature. In this particular example, we compare the estimation of the temperature with and without the sky temperature correction to thermocouples measurements at the surface of a concrete pavement structure. The sky correction reduces in average the difference between the thermocouples measurements and the estimated temperature from infrared data.

Figure 4 . 19 :

 419 Figure 4.19: Comparison of ground truth temperature values to temperature estimation with and without sky correction

Figure 4 . 20 :

 420 Figure 4.20: Effect of using a correlation for the sky temperature based on the air temperature

Figure 4 . 21 :

 421 Figure 4.21: Sky temperature derived from different sources of air temperature

Figure 4 . 22 :

 422 Figure 4.22: Difference between true object temperature and estimated one depending on the error made on the atmospheric transmission. Here is the difference map for an object temperature of 293.15K and emissivity of 0.93.

Figure 4 . 23 :

 423 Figure 4.23: Different distance values for the atmospheric transmission, through time.

Figure 4 . 24 :Figure 4 . 25 :

 424425 Figure 4.24: Sun irradiance reflection on 7.5µm -13µm band for a normal object and different emissivity values

Figure 4 . 26 :Figure 4 . 27 :

 426427 Figure 4.26: Comparison of sun's irradiance measured and simulated on may 5th in the 0.3µm -2.8µm band

Figure 4 .

 4 Figure 4.28: IR image view: in dry environment (a), rain with optic issues (b)

Figure 4 .Fig. 4 .

 44 Figure 4.29: A set of environmental data for long-term thermal monitoring

Figure 4 . 30 :

 430 Figure 4.30: Comparison between the temperature estimation by considering M tot = ϵ obj M obj or using the presented model and multi-sensor data

Figure 5 . 1 :

 51 Figure 5.1: Form factors for two infinitesimal elements

5 )Figure 5 . 2 :

 552 Figure 5.2: Nusselt's analog

Figure 5 . 3 :

 53 Figure 5.3: Gathering step as in the Gauss-Seidel algorithm, figure inspired from[START_REF] Michael | Radiosity and Realistic Image Synthesis[END_REF]

Fig. 5 .

 5 [START_REF] Aires | Remote Sensing from the Infrared Atmospheric Sounding Interferometer Instrument 2. Simultaneous Retrieval of Temperature, Water Vapor, and Ozone Atmospheric Profiles[END_REF] shows an example of the progressive radiosity algorithm, for which the implementation is detailed in the next section.

Figure 5 . 5 :

 55 Figure 5.5: Progressive radiosity examples with a Cornell's Box where all surfaces emit in [1.5µm -5.4µm] band. (a) Initial state. The ceiling source is at 353.15K and the objects at 293.15K with different emissivities. (b) After 10 iterations (≈ 0.15s). (c) After 100 iterations (≈ 0.3s). (d) After 1000 iterations with interpolation to smooth the result (≈ 4.4s). Almost all surfaces have shot their unshot radiosity (red texture).

Figure 5 . 7 :

 57 Figure 5.7: Mirror's Edge Catalyst driven by Unreal Engine 3, which can interface OpenGL ® -Credits to deadendthrills.com

Figure 5 . 8 :

 58 Figure 5.8: Principle of the progressive radiosity GPGPU implementation

22 )Figure 5 . 9 :

 2259 Figure 5.9: Projection image of one of the plane from an element of the other plane

Figure 5 . 10 :

 510 Figure 5.10: Concurrency example. At first, the maximum value is 10 at index 0. The concurrency makes the orders in the queues not the same. The maximum value is always guarantee but not its corresponding index

Figure 5 . 11 :

 511 Figure 5.11: Energy shot at a given time of the algorithm as a percentage of the first shot

γ)( 5 . 26 )

 526 s = sign(ω) arccos ( cos(θ z ) sin(φ)sin(δ) sin(θ z ) cos(φ) An example of the change of the sun position over the sky vault in the software is shown in Fig. 5.13.

Figure 5 . 13 :Figure 5 . 14 :

 513514 Figure 5.13: Sun position over the sky vault for a 3D wood house scene

Fig. 5 .

 5 15a shows the impact of the zenith angle to the received direct irradiance component whereas Fig.5.15b shows the impact of the visibility.

Figure 5 . 15 :

 515 Figure 5.15: Spectrum of the solar irradiance received by a surface normal to the sun for different visibility and altitude values

Figure 5 . 16 :

 516 Figure 5.16: Transmittance values for different zenith angles

  ) and units tests for the different look-up tables interpolations algorithms performed.

Figure 5 . 17 :

 517 Figure 5.17: Effect of the pathlengths of ozone and N O 2 to their transmittance and comparison with the original figure.

  5.18).

Figure 5 . 18 :

 518 Figure 5.18: Animation of an object in the 3D scene through Python: an object on the scene is rotated of 4 degrees and the scene is displayed on screen at each iteration.

Fig. 5 .Figure 5 . 19 :

 5519 Fig. 5.[START_REF] Berk | MODTRAN: A Moderate Resolution Model for LOWTRAN[END_REF] shows an example of radiosity map that can be obtained in a static configuration with a building (near infrared band to view the effect of the sun). The sun contribution can be observed among with the created shadow due to this illumination. The meshing impact on the final resolution and the interpolation quality is also visible, and very documented in the literature[START_REF] Michael | Radiosity and Realistic Image Synthesis[END_REF]. The Python script used to generate the simulation is given in Appendix C

20 Figure 5 . 20 :

 20520 Figure 5.20: 3D model of the boogie in its scene environment

Figure 5 . 21 :

 521 Figure 5.21: Boogie and environment in the simulation software. The preview of the rendered infrared image is done in real-time at the bottom left corner.

Figure 5 . 22 : 5 . 6 . 6

 522566 Figure 5.22: Comparison of simulation data with experimental ones. The motion blur is well represented in the simulation. The temperature levels however are not well represented due to a lack of information about the initial temperatures values of the boogie.

Figure 5 . 25 :

 525 Figure 5.25: Flowchart overview of the developed software showing the different steps of the numerical simulation

  6.1).

Figure 6 . 1 :

 61 Figure 6.1: Temperature and emissivity profiles

Figure 6 .

 6 Figure 6.2: 6-order Tchebychev approximation through Kalman filtering

Figure 6 .

 6 Figure 6.3: 6-order Tchebychev polynomial approximation through combined Kalman filter and gradient descent.

Figure 6 . 4 :

 64 Figure 6.4: Kalman filter flowchart overview

Figure 6 . 5 :Figure 6 . 6 :

 6566 Figure 6.5: Temperature and emissivity estimation through the interacting Kalman filter (8000 particles)

35 )

 35 In order to get the UK formula, we take the limits of µ(s * |S) and σ 2 (s * |S) when the prior on the trend weights β ∼ N (b, B = λ B), λ ∈ R becomes non informative. That is to say when λ → ∞ and b = 0. Therefore, we define the estimates as: μ(s * ) = lim λ→∞ b=0 µ(s * |S) σ2 (s * ) = lim λ→∞ b=0 σ 2 (s * |S) (6.36)

38 )

 38 By multiplying the previous equation on the left by F (s * )BF (S) T and on the right by y: F (s * )BF (S) T (C + F (S)BF (S) T ) -1 y = F (s * )B((F (S) T

Figure 6 . 7 :

 67 Figure 6.7: Influence of the covariance kernel on the prediction (in black). The 95% confidence interval in light gray and the input measurements points in red.

Figure 6 . 8 :

 68 Figure 6.8: Influence of the trend on the prediction (in black). The 95% confidence interval in light gray and the input measurements points in red.

Figure 6 . 9 : 6 . 7 . 1 A

 69671 Figure 6.9: Difference between a spherical and Matérn variogram models for a 2D problem

Figure 6 . 10 :

 610 Figure 6.10: Kriged Interacting Kalman filter flowchart overview

Figure 6 . 11 :

 611 Figure 6.11: Multi-resolution approach for solving the optimization problem. From left to right: coarse level to highest final resolution details.

Fig. 6 .

 6 Fig. 6.14.

Fig. 6 .

 6 Fig.6.12 shows the comparison of the two algorithms. The first row shows the result for a constant emissivity over the entire image. The second one shows the result of the trust-region-reflective algorithm and the last one the result obtained with the CMA-ES algorithm. One can see that the result of the CMA-ES differs from the two others. Whereas in the lsqnonlin, a local solution is computed where the emissivity is almost constant, the CMA-ES tends to a more homogeneous temperature field.

Fig. 6 .

 6 [START_REF] Baldridge | The ASTER Spectral Library Version 2.0[END_REF] shows the processing time for each algorithm, depending on the number of pixels in the image, for a 3 × 3 window size. One can see that depending on image's size, the processing time is quite large for the CMA-ES algorithm. on a 41 × 36 image, the CMA-ES takes over 2 minutes and 40 seconds (just for this image, not for the entire multi-scale pyramid's process). For the full image (480 × 640), it represents more than 10 hours of process, compared to 1 hour and a half for the lsqnonlin function.

Figure 6 . 12 :

 612 Figure 6.12: Comparison of the trust-region-reflective algorithm and the CMA-ES

Figure 6 . 14 :

 614 Figure 6.14: Comparison of processing time, for one given scale with a given number of pixels

Figure 6 . 15 :

 615 Figure 6.15: Target used for the study case. (a) Rendering of the scene in the visible spectrum. The camera is represented, the target and the hemisphere for the environment. (b) Result example for a frame at 20 • C and other materials at 5.2 • C

Figure 6 . 16 :

 616 Figure 6.16: Emissivity spectrum of the materials of the target.

Figure 6 . 17 :

 617 Figure 6.17: Target material's temperature evolution through time. The measurements correpond to the temperature at ground during two days of january 2017. The decimated temperature profile is shifted by 1K for drawing purposes.

Fig. 6 .

 6 18 and Fig. 6.19 show the result obtained with the algorithm after 10k iterations, for a single pixel value of each material for each spectral band. The histograms Fig. 6.19b, Fig. 6.19d, Fig. 6.19f and Fig. 6.19h show a shift of the mean in the posterior distribution of the temperature. The temperature and emissivity simultaneous

Figure 6 . 18 :

 618 Figure 6.18: Focus on the emissivity retrieval on one pixel of each of the four materials for the 8µm -12µm spectral band.

Figure 6 . 19 :

 619 Figure 6.19: Results of the temperature / emissivity retrieval on one pixel of each of the four materials. (a) (c) (e) (g): Emissivity retrieval. The ground truth profile (continous) is compared to the estimated one (per spectral band).(b) (d) (f) (h): Temperature estimation. The ground truth temperature is represented by the dashed line.

Figure 6 . 21 :

 621 Figure 6.21: Coarse level preliminary CMAES result. The image is blurry due to the level used. Band: 10µm -12µm.

Fig. 6 .

 6 Fig.6.[START_REF] Brunet | On the Mathematical Properties of the Structural Similarity Index[END_REF] shows the final result at fine resolution. Due to the window size of 5 × 5 pixels, the result looks noisy. This noise could be removed either by increment the window size to an extent that would not compromise the results or, as done here, by a post-process filtering (median filter).

Figure 6 . 22 :

 622 Figure 6.22: Full result of the CMAES algorithm at a given time and given band. At the opposite of Fig. 6.21, the image is noisy due to the small neighborhood involved. Band: 10µm -12µm.

. 60 )

 60 Fig.6.23, the different emissivity materials levels are found. However, some differences exist between the expected value and the estimated one. By looking at the RMSE of the emissivity, an important gap exists between the bottom right material (alumina) emissivity expected and estimated values. However, the SSIM shows that this estimated solution is found over the entire material and that the spatial structure of the emissivity map corresponds to the expected values.

Figure 6 .

 6 Figure 6.24: CMAES: SSIM and RMSE local values between estimated emissivity and temperature and ground truth values.A small error on materials with low emissivity leads to a more important error in temperature than materials with higher emissivity. The material average RMSE is up to 4K against 0.05 for the emissivity, which is still lower in average than for a constant emissivity case. Band: 10µm -12µm.

Fig. 6 .

 6 [START_REF] Cerezo | A Survey on Participating Media Rendering Techniques[END_REF] shows the result of the first kriging step, with this first T 0 . The first image in this figure shows the estimated emissivity at T 0 and the second and third one the kriging interpolation from this emissivity map and the references points.

Figure 6 . 25 :

 625 Figure 6.25: Initial estimate by using a constant temperature value T = T 0 = 270K. The emissivity at T 0 is shown in the first image. The Kriged model result from this emissivity map is shown in the two other images with the mean estimate and standard deviation. Band: 10µm -12µm.

Figure 6 . 26 :

 626 Figure 6.26: Final estimate over the entire period. The ground truth values are compared to the estimates. As for the CMA-ES the materials emissivities are distincts. Band: 10µm -12µm.

Figure 6 . 27 :

 627 Figure 6.27: Estimated temperature time evolution at one pixel

Figure 6 .

 6 Figure 6.28: KIKF: SSIM and RMSE local values between estimated emissivity and temperature and ground truth values.A small error on materials with low emissivity leads to a more important error in temperature than materials with higher emissivity. The material average RMSE is less than 3K against 0.04 for the emissivity, which is still lower in average than for a constant emissivity case. Band: 10µm -12µm.

  Application ϵ and T independent Per pixel MCMC Single-band Single sample ϵ and T independent Per pixel CMA-ES T homogeneous on a given neighborhood Single-band Single sample ϵ and T independent Per pixel IKF Small bandwidth Multi-band only Needs temporal data ϵ and T independent

  2.1 Main steps of structural health monitoring process . . . . . . . . . . . . . . . . . . . . . 2.2 Focus on the IR part of the electromagnetic spectrum. Common spectral bands of IR are shown: Near-infrared (NIR: 0.75 -1.4µm); Short-wavelength infrared (SWIR 1.4 -3µm); .13 The geometry of the scene leads to different angle of view between the object and the camera in the image . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.14 Different views of the same scene, the spatial sampling is not constant within the images 2.15 Mixed pixel effect illustration with three different materials. . . . . . . . . . . . . . . . . 2.16 Small uncooled IR camera with 3D-printed case. A coin of 1e is placed at the top as a scale factor. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.17 Cross-sectional view of a bolometer pixel structure principle . . . . . . . . . . . . . . . . 2.18 Different types of detectors reading . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Radiative transfer equation for satellite measurements, done at the top of the atmosphere. 37 3.2 ln approximation with its Taylor serie . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 3.3 Original irradiance at sensor compared to the irradiance with the estimated emissivity and temperature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.4 Found temperature and emissivity profiles . . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.5 Cubic spline interpolation with quadratic boundaries conditions. . . . . . . . . . . . . . 49 3.6 Ray tracing images can be realistics, example with the POV-Ray software rendering . . . 61 3.7 Example of ray-traced rendering with the Sponza scene and Cycles open source renderer. 62 3.8 Monte-Carlo π approximation after 500 iterations and 10000 iterations. . . . . . . . . . 63 3.9 Hierarchical subdivision . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

	Mid

-wavelength infrared (MWIR 3 -8µm); Long-wavelength infrared (LWIR 8 -15µm); Far infrared (FIR 15 -1000µm). The atmospheric transmission is also represented. . . . . . . 2.3 Planck's law as a function of wavelength at different temperature values . . . . . . . . . 2.4 Temperature as a function of exitance at different wavelengths. . . . . . . . . . . . . . . 2.5 Solid angle on hemisphere. Note that the solid angle is not necessary a circular cone . . 2.6 Emission rate in a normal direction and off direction with the Lambert's cosine law on the right. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7 Interactions that occur at the interface of an object . . . . . . . . . . . . . . . . . . . . . 2.8 Spectral emissivity of a series of various materials. (Data from reflectance spectra in ASTER spectral library http://speclib.jpl.nasa.gov, Copyright: Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA; see also Balridge, A.M. et al., Remote Sens. Environ., 113, 711, 2009.). The alumina is from previous studies [130]. . . . 2.9 Effect of the state of the surface on emissivity (from [71]) . . . . . . . . . . . . . . . . . . 2.10 Enclosure for which the radiosity is computed at the surface. Point x 3 and x 2 cannot see them each other. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.11 Main radiatives contributions received by the infrared camera . . . . . . . . . . . . . . . 2.12 Possible solutions for a given measurement when everything but temperature and emissivity are known . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

22.19 Sensor / Object configuration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.20 Multi-band measurements bring more information . . . . . . . . . . . . . . . . . . . . . viii 3.1 4.1 Black-body calibration source . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.2 Result of thermal calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.3 Principle of camera resectioning. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

  The emissivity at T 0 is shown in the first image. The Kriged model result from this emissivity map is shown in the two other images with the mean estimate and standard deviation. Band: 10µm -12µm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173 6.26 Final estimate over the entire period. The ground truth values are compared to the estimates. As for the CMA-ES the materials emissivities are distincts. Band: 10µm -12µm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 174 6.27 Estimated temperature time evolution at one pixel . . . . . . . . . . . . . . . . . . . . . 174 6.28 KIKF: SSIM and RMSE local values between estimated emissivity and temperature and ground truth values. A small error on materials with low emissivity leads to a more

6.21 Coarse level preliminary CMAES result. The image is blurry due to the level used. Band: 10µm -12µm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.22 Full result of the CMAES algorithm at a given time and given band. At the opposite of Fig. 6.21, the image is noisy due to the small neighborhood involved. Band: 10µm -12µm.170 6.23 Filtered CMAES result after a 2D median filter. The materials have distinct emissivities values on the emissivity map. As stated in Chapter 4 and observed on the temperature absolute difference image, a small error on low emissivity values can lead to important error in the temperature estimate. Band: 10µm -12µm. . . . . . . . . . . . . . . . . . . 6.24 CMAES: SSIM and RMSE local values between estimated emissivity and temperature and ground truth values. A small error on materials with low emissivity leads to a more important error in temperature than materials with higher emissivity. The material average RMSE is up to 4K against 0.05 for the emissivity, which is still lower in average than for a constant emissivity case. Band: 10µm -12µm. . . . . . . . . . . . . . . . . . xii 6.25 Initial estimate by using a constant temperature value T = T 0 = 270K. important error in temperature than materials with higher emissivity. The material average RMSE is less than 3K against 0.04 for the emissivity, which is still lower in average than for a constant emissivity case. Band: 10µm -12µm. . . . . . . . . . . . . . . . . . 175 D.1 A day forecasting with the SARIMA model used (2, 0, 0) × (0, 1, 1) 1 and its comparison to actual measurements with the 95% confidence interval . . . . . . . . . . . . . . . . . 203 D.2 A day forecasting with the SARIMA model used (2, 0, 0) × (0, 1, 1) 1 and its comparison to actual measurements with the 95% confidence interval . . . . . . . . . . . . . . . . . 203 xiii Chapter 1

Definition 2.2.9 (Exitance

  ). The radiant exitance denoted M is defined as the local value of the ratio of the flux ϕ e emitted by an object and the area of this object. It is the power emitted per unit area

	Definition 2.2.8 (Irradiance). The irradiance denoted E is defined as the local value of the ratio of
	the flux ϕ r received by an object and the area of this object. It is the power received per unit area
	(W • m -2 ):
	E =	∂ϕ r ∂A	(2.16)
	(W • m -2 ):
	M =	∂ϕ e ∂A	(2.17)
	Definition 2.2.10 (Radiosity). The radiosity of a surface denoted J is the radiant flux leaving (emitted,
	reflected and transmitted by) a surface per unit area. It is the power left per unit area (W • m -2 ).

  tion. Therefore, if we have n > 6 correspondences, a more robust estimation can be made by performing a Random sample consensus (RANSAC)[START_REF] Fischler | Random Sample Consensus: A Paradigm for Model Fitting with Applications to Image Analysis and Automated Cartography[END_REF] algorithm. A pseudo-code of a generalized RANSAC algorithm for this purpose is given on Algorithm 5.

	Algorithm 5: Pseudo-code for RANSAC algorithm where the error function and the model are
	defined appart.
	Function RANSAC(data, fcnModel, fcnError, nData, nIteration, thres, nMinData) data ObservationsB set of observations ( ObservationsA )
	fcnModel	model to be adjusted to the data, returns the parameters p
		of the model (p = fcnModel(inputP oints))
	fcnError	estimate error between a point and the parameters: e = fcnError(point, p)
	nData	number of data used to adjust the model
	nIteration	number of maximum iterations
	thres	a threshold value to adjust the data to the model
	nMinData	number of minimal data necessary to consider this model
	Initialization iIteration = 0, finalError = ∞, other allocations here
	while iIteration < nIteration
	rPoints = select nData random points from data;
	modelParameters = fcnModel(rPoints) ;	/* Get those parameters */
	setOfPoints = rPoints ;	/* We initialize the set of points */
	foreach iPoint ∈ data and iPoint / ∈ rPoints /* For all other points	*/
	if fcnError(iPoint, modelParameters) < thres
	add iPoint to setOfPoints;
	/* If the number of elements is satisfying	*/
	if numberOfElements(setOfPoints) > nMinData
	e = fcnError(setOfPoints, modelParameters) ; /* Compute the error with all
	the points */
	if e < finalError /* Best model yet	*/
	/* Save the data	*/
	finalModelParameters = modelParameters;
	finalSetOfPoints = setOfPoints;
	finalError = e;
	iIteration++;
	Result: finalModelParameters, finalSetOfPoints, finalError

Table 4 .

 4 1 shows the associated emissivity measurements.

	Material	Color	Emissivity
	Wood wall (west) Blue	0.95
	Wood wall (east) Yellow	0.90
	Grass	Green	0.95
	Sidewalk	Red	0.89
	Tar	Grey	0.95
	Roof	Purple	0.92

Table 4 .

 4

1: Emissivity measurements values in the 8µm -14µm band.

  lation of the problem. Let define the state-space model of a linear system with i ∈ N inputs, o ∈ N outputs and n ∈ N so called state variables. A general state-space representation is given in the fol-

	lowing equation:
	ẋ(t) = A(t)x(t) + B(t)u(t)
	(6.1)
	y(t) = C(t)x(t) + D(t)u(t)

  update step involves computing p(x t |Y t ) given a prediction p(x t |Y t-1 ) and a measurement y t . By applying Bayes rule: Again, the last line is obtained through independence assumptions. That is to say, the observations only depends on the current state: p(y t |x t , x t-1 , . . . , x 0 , Y t-1 ) = p(y t |x t ). By applying Chapman-

	p(x t |Y t ) = =	p(x t |Y t-1 )p(y t |x t , Y t-1 ) p(y t |Y t-1 ) p(y t |Y t-1 ) p(x t |Y t-1 )p(y t |x t )	(6.7)
	Kolmogorov equations, one can get
		∫
	p(y t |Y t-1 ) =	p(y t |x t )p(x t |Y t-1 )dx t	(6.8)

  x t |x t-1 )p(x t-1 |Y t-1 )dx t-1 prediction p(x t |Y t ) = p(x t |Y t-1 )p(y t |x t )

			(6.9)
	∫	p(y t |x t )p(x t |Y t-1 )dx t	update

Caractérisation conjointe de la température et des propriétés ra- diatives des objets par thermographie infrarouge multispectrale Mot clés : thermographie

  infrarouge, monitoring thermique, propriétés thermo-optiques, approche Bayésienne, filtres de Kalman, radiosités progressives Résumé : L'utilisation de caméras infrarouges bas coûts pour la surveillance long-terme d'infrastructures est prometteuse grâce aux dernières avancées technologiques du domaine. Une mesure précise de la température des surfaces observées in-situ se heurte au manque de connaissance des propriétés radiatives de la scène. L'utilisation d'une instrumentation multicapteurs permet d'affiner le modèle de mesure afin d'obtenir une estimation plus précise de la température. A contrario, il est montré qu'il est toujours possible d'exploiter des données climatiques en ligne pour pallier un manque de capteur. Des méthodes Bayesiennes d'estimation conjointe d'émissivité et de température sont ensuite développées et comparées aux méthodes de la littérature. Un simulateur d'échanges radiatifs diffus de scènes 3D a été implémenté afin de tester ces différentes méthodes. Ce logiciel utilise l'accélération matérielle de la machine pour réduire les temps de calcul. Les résultats numériques obtenus mettent en perspective une utilisation avancée de la thermographie infrarouge multi-spectrale pour la surveillance de structures. Cette estimation conjointe permet alors d'obtenir un estimé de la température par thermographie infrarouge avec une incertitude connue.

	Titre:
	#=================================================================
	# This is a sample YAML file for running batch processes
	# [o] stands for optionnal
	properties:	# [o] -> properties of the simulation
	recursive:	#	-> Either to look onto subfolder or not.
	showprogress:	#	-> show progress of the process
	logfile:	# [o] -> save the output to a log file
	input:	#	-> tag for defining the IR data
	folders:	#	-> can be a path to a file or a folder
	date:	# [o] -> defines the range of time to consider
	beginDate:	#	->
	endDate:	#	->
	fileFormat:	#	-> the format of date in the filename
	fileDateExtract: #	-> the way to extract the date in the filename (regexp)
	output:	#	-> Initialize the fields of outputData
	actionList:	#	-> which actions should be executed in which order
	myaction:	#	-> the action name
	...
	myaction:	#	-> The actual definition of the action "myaction"
	func:	#	-> The function that will be executed
	type:	#	-> batch: will be executed for every single file
		#	-> before/after: will be executed in the main program
	frameRate:	# [o] -> FrameRate if batch
	beginDate:	# [o] -> Begin date if batch
	endDate:	#	-> End date if batch
	args:	# [o] -> arguments that will be given to the function
	extra:	# [o] -> extra arguments. For clarity, you can pass
	name1: value1	#	arguments with the tag 'extra'. Its child tags
	name2: value2	#	will be parsed as 'name', 'value' arguments.
	#=================================================================

dA ′(2.42)

The superscript(k) denotes the k th complete iteration. The superscript (0) corresponds to a prior given as input.

Remerciements

Multi-temperature optimization

One way to overcome the undetermined system is to perform multiple measurements, at multiple temperatures, assuming that the emissivity will not change. In principle, such assumption makes the equations' system solvable by removing one temperature variable. This multi-temporal method seems attractive in theory. However, this method is highly sensitive to measurement noise due to the constant emissivity hypothesis between the two measurements [START_REF] Li | Satellite-Derived Land Surface Temperature: Current Status and Perspectives[END_REF][START_REF] Krapez | Radiative Measurements of Temperature[END_REF].

Similarly, multi-temperature has been used in the Remote Sensing field by assuming that the emissivity is invariant at two times. However, this method requires accurate geometry registration to get observation under similar angles (which favors geostrationary satellite data). Moreover, it requires multitemporal data and important temperature gradient between the measurements.

Note

Grey body hypothesis

The grey body emissivity method tackles also the equation system unknown issue by assuming that on two given area of the spectrum, the emissivity is constant. Such flat region in the emissivity profile is actually a strong requirement and is not applicable to most of the observed materials.

Similarly, grey body methods have been used in the Remote Sensing field [START_REF] Barducci | Temperature and Emissivity Retrieval from Remotely Sensed Images Using the" Grey Body Emissivity" Method[END_REF] but for the same reason, the assumption is not satisfyed in most cases.

Note

Bayesian methods

In this section are presented the Bayesian inference based methods. Those statistical methods are based on the Bayes' theorem, used to update the probability for a hypothesis as more information becomes available. With the increase of computation means, those methods have become common. More details on Bayesian inference are given on Chapter 6.

Dynamic filtering with Kalman and particle filters

The Kalman Filter applied to the emissivity / temperature retrieval has been used in the literature, in [START_REF] Masiello | Kalman Filter Physical Retrieval of Surface Emissivity and Temperature from Geostationary Infrared Radiances[END_REF][START_REF] Venafra | RETRIEVAL OF SURFACE TEMPERATURE AND EMISSIVITY FROM SEVIRI DATA: IMPLEMENTATION OF KALMAN FILTER METHODOLOGY TO FULL DISK AND VALIDATION WITH IN-SITU OBSERVATIONS[END_REF] with SEVIRI data for example. In particular [START_REF] Masiello | Kalman Filter Physical Retrieval of Surface Emissivity and Temperature from Geostationary Infrared Radiances[END_REF] formulates a four-dimensional

Ray tracing

Ray tracing methods principle is to browse the light path backward from reality: a ray is thrown from the camera toward the 3D scene and bounces until it reaches the light. This behavior is the opposite than reality where the ray goes from the light to the camera. Each time a ray encounter an obstacle, its energy can be distributed according to the different physical phenomena that occur on the object (absorption, reflection, transmission, diffraction, etc.). The first algorithm was developed in 1968 by Arthur Appel from IBM [START_REF] Appel | Some Techniques for Shading Machine Renderings of Solids[END_REF]. The idea was to throw a ray from the eye, one per pixel and stop once the closest object blocked the ray. This led to determine the shading of an object. Then, in 1979, the algorithm was improved by performing multiple bounces and casting new rays from an incoming ray that represents the reflection, refraction and shadow [START_REF] Whitted | An Improved Illumination Model for Shaded Display[END_REF]. Then, improvements on the method and the consideration of many other phenomena were made on the literature [START_REF] Dutré | Advanced Global Illumination[END_REF]. A few methods for ray tracing based rendering are presented. Basic ray tracing algorithm are straightforward to implement nowadays (for simple considerations of physical phenomena) and can provide realistic results (see Fig. 3.6 for an advanced result). Figure 3.6: Ray tracing images can be realistics, example with the POV-Ray software rendering to test the previous temperature and emissivity separation methods. Finally, the estimation of the parameters other than emissivity and temperature can also help to improve the estimation accuracy.

The next chapter aims at introducing the concepts of practical infrared thermography measurements and the impact of the environmental parameters on the final temperature estimation. tioning. Some articles can be found in the literature for matching textured 3D models and 2D images ( [START_REF] Wu | 3D Model Search and Pose Estimation from Single Images Using VIP Features[END_REF], [START_REF] Wu | 3D Model Matching with Viewpoint-Invariant Patches (VIP)[END_REF]). [START_REF] Dubreuil | Calibration Sur Modèle CAO Pour La Mesure Par Vision In-Situ[END_REF] have proposed an algorithm for performing a calibration based on the CAD of the observed scene. This approach seems interesting but unfortunately has not been tested in this thesis.

Summary of main resectioning algorithms

Note

Spatial sampling

Another topic that is also challenging is the extraction of profiles from infrared images. In fact, it can be interesting to extract thermal profiles from infrared images to compare them to theoretical models for example. In particular, the error made on the profile compared to the theoretical thermal one can be complicated to compute.

Monitoring a civil engineering structure requires a large angle of view which forces users to use a mast to place the camera. This situation results in images with varying spatial resolution and important angle between the object and the camera (e.g. Fig. 2.14). The underlying consequences are non-square pixels in the real-world and varying angles values along the same plane which has an impact on the emissivity values to be considered. One way to overcome this is to perform spatial calibration on images

The software can therefore read and import different data through multiple file formats (.hdf; proprietary IR: .seq, .fff; binaries IR: .irl; common ones: .csv, .txt etc.) and make the conversion between any of those formats. The importation and formatting is facilitated by a custom importation module (see Fig.

Data processing

Conversion process

The tool supports different data processing algorithm. The conversion model from digital levels to temperature can be customized and the parameters derived from the thermal calibration can be sets.

The defaults values are the one stored into the metadata of the input file. Even though the conversion from digital levels to temperature is not a difficult computation to perform for MATLAB ® , it may be too long for fast visualization purposes. In fact, when displaying the images and browsing from them, any latency can make the user's experience complicated. As a consequence, the computation speed has been improved by the parallelization of the conversion model. The digital levels to temperature conversion is then quickly done and can be customized depending on user's hardware. Three choices are available, from the slowest to the quickest: CPU only, GPU through gpuArray's built-in MATLAB ® functions and CUDA using custom GPU kernel. Compute Unified Device Architecture (CUDA) is a General-Purpose Computing on Graphics Processing Units (GPGPU) technology. Due to the constant improvements of graphic hardware for 3D graphics compute-intensive real-time applications, the GPUs became really efficient for performing highly parallel tasks on large blocks of data. Therefore, this design led to increase the number of applications and therefore derive the original purpose of 3D graphics to scientific general computations. CUDA was therefore developed by NVidia for their hardware to achieve such goals and interface directly with the GPU.

Numerical solution to the linear radiosity system

Once the form factors have been computed, the radiosity system can be solved. The Eq. (5.2) shows that the radiosity system is a set of linear equations which can be summarized by:

where A is the interactions matrix, x a vector of radiosity values and b a vector of emission terms. By combining Eq. (5.2) and Eq. (5.6), we can get the shape of A:

Most of the time, the diagonal coefficient of A are equals to 1. In fact, one assumption is to consider that one element cannot "see" itself (in practical, the meshing is supposed to be refined enough and the involved patches are often planar surfaces). Therefore, in the following for all i ∈ N, F i→i = 0.

Note

Many numerical methods for solving linear equation's systems can be found in the literature. It should be noted here that the system involved in the radiosity is a large matrix most of the time. Furthermore, when dynamic scenes are involved, the main difficulty is that the matrix will be changed and that it should be inversed again. As a consequence fast algorithms are needed. Solving such system is a large active field on its own, so only a review of the main algorithms used in the literature and applied to the radiosity's equation is proposed here. The next section focuses on methods for solving the linear system Eq. (5.6) of unknown x.

Solving the linear system

Once the system has been solved, the radiosity values are found and the rendering is completed. As a consequence, finding the exact solution in a fast way seems to be the best option. However, finding the inverse of A is non trivial and simple inversions methods (e.g. Gauss-Jordan elimination or LUdecomposition) will not work to solve the system in a reasonable amount of time. Those algorithms are usually in order of O(n 3 ) which is not suitable for the large matrix that is involved in the radiosity equation system. Instead, iterative methods that provide an approximation of the solution will be preferred. The Southwell method is actually a successive over-relaxation (SOR) method.

A SOR method is a variant of the Gauss-Seidel one on which a relaxation factor is incorporated and the A matrix is defined by A = D + L + U where D is diagonal, L lower and U upper triangulars matrices.

For a positive definite matrix, the algorithm converges for ω ∈ [0, 2]. For ω = 1, the equation is equivalent to Gauss-Seidel.

Note

Progressive Radiosity

Based on the Southwell algorithm, [START_REF] Michael | A Progressive Refinement Approach to Fast Radiosity Image Generation[END_REF] proposed a similar approach to render the radiosities which enables both computation efficiency and rendering called progressive radiosity or sometimes progressive refinement. As for the Southwell method, the progressive radiosity method processes the patches in a sorted order according to their energy contribution to the environment. This approach avoid computing all the form factors at the beginning of the algorithm. Thus, each iteration needs the computation of the form factor from the current patch to all other patches which can, however, lead to redundant computations. The pseudo-code for the progressive refinement method is presented in Algorithm 9.

The Southwell method is slightly changed so that the result can be shown progressively. To achieve this, an additional variable is introduced: unshot radiosities which can be stored during the process.

Taking into account the unshot radiosities is a way to take into account what part of energy remains to

Algorithm 10: Proposed code for progressive radiosities algorithm on hardware Function RADIOSITY()

While stopCriteria unsatisfied hemiProjection = nonLinearProjection(i);

For each pixels j ∈ hemiProjection ; //performed on the compute shader

A j ; radiosity j = radiosity j + ∆B; residuals j = residuals j + ∆B;

can be avoided by different methods such as message-passing where the requested are queued.

Many possibilities can be investigated to overcome the concurrency on the GPU. The most efficient one would probably be to interface with CUDA or OpenCL libraries and a specific kernel. One drawback of this method is that the code becomes more complicated and an overload of includes libraries just to find a maximum value. Another possibility is to get the texture back to the CPU and find the maximum on the CPU but the efficiency of this method may vary between hardwares and breaks the GPU-only process. Instead, the chosen solution is to use again atomic counters thanks to the GL_ARB_shader_atomic_counter_ops extension which incorporates an atomic max operation.

Atomic counter: a specific openGL variable type which can be seen as a limited buffer image variable on which atomic memory operations can be performed on them. Therefore, no memory barriers are needed to synchronize their access: each read/modify or write operation is completed in its entirety before other operations are permitted.

The current implementation uses two atomic counters, one for keeping the maximum value and one for keeping the index of the maximum. Such design leads to a lack of synchronization in-between the counters:

1. atomicCounterMax function replaces the value of the max atomic counter if the new value is greater than the oldest one.

atomicCounterExchange exchanges the old value in the index atomic counter by the new index

Therefore, if the atomicCounterMax function always guarantees to store the maximum, the atomicCounterExchange does not guarantees that the index actually corresponds to the max found.

In fact, even if the queues of the counters are synchronized by themselves, the maximum and exchange queues are not synchronized together, as shown in Fig. 5.10.

The algorithm is different from the progressive refinement on which the maximum of energy is not always chosen but one of the maximum patch energy is instead. Concurrency can occur when two another one that goes onto the red. In Fig. 5.24, the first top left zoom shows the diffusion that occurs between the red wall to the tall box, as expected. However, the classical issues of the radiosity method can also be found and are well addressed in the literature [START_REF] François | Radiosity and Global Illumination[END_REF]. An example is the wrong shadows, on the second zoomed box at the lower right corner of Fig. 5.24. Such issue is due to the too coarse mesh in this region and the aliasing algorithm. It can be prevented by modifying the aliasing method and/or modify the mesh so that it is more refined where the gradient is important [START_REF] Heckbert | Adaptive Radiosity Textures for Bidirectional Ray Tracing[END_REF][START_REF] Vedel | A Testbed for Adaptive Subdivision in Progressive Radiosity[END_REF][START_REF] Paulin | Adaptive Mesh Generation for Progressive Radiosity: A Ray-Tracing Based Algorithm[END_REF].

Figure 5.24: Cornell's box result for two different lights. Such result is computed on a professional laptop -M1000M Nvidia graphic card, Intel i7-6820HQ CPU @ 2.70GHz. The computation time is approximately 60 seconds with a progressive rendering and aliasing on screen every 50 iterations.

Kriging

In the previous section, the temporal and multi-band parts of the data have been used to get an estimate of the temperature and emissivity from infrared image. However, one interesting part in those data has not been taken into acount yet: the spatial information of the image. Spatial issues have been presented in Chapter 4, in terms of discretization, angle of view and calibration. This section focuses on extracting spatial information on the images. Moreover, even if the interactive Kalman filter has shown interesting results on a simple sensor data, the required computation time to process an entire image may not be convenient for in-situ practical applications. This section focuses on exploiting the spatial information to serve two purposes:

• Bringing spatial information to the temperature/emissivity retrieval,

• Get rid of the short spectral band distance requirement,

• Reducing the required computation time.

Such goals have been reached by using a Kriging model on top of the interactive Kalman filter. The kriging method was developed in the 60's by Georges Matheron, from Danie G. Krige works in the geostatistical field. Kriging is a spatial interpolation method that provides estimates at unknown locations based on sampled ones. Under some assumptions, kriging gives the best linear unbiased prediction for the intermediate values. Based on the modeling of the experimental variogram, the spatial interpolation is made by computing the expected value of a random variable. We will introduce first the Kriging equations and in the next section present the final temperature / emissivity retrieval method KIKF (Kriged Interactive Kalman Filter).

Let S ⊂ R 2 be the studied spatial domain, S ⊂ S be the measurement domain and S * ⊂ S be the spatial domain where lies the points on which the estimation has to be made such as S ∪ S * = S .

Variogram

Definition

Let s ∈ S a spatial point with its associated measurement y(s). Consider now the measure-

The variogram is defined as the expectation of the random variable (y(s)y(s + h)) 2 [START_REF] Howarth | Mining Geostatistics[END_REF]. Since in practice only one realization of the couple (y(s), y(s+h)) is known, the hypothesis behind the modelization of the variogram is to assume that the variogram only depends on the modulus and direction h and not its location s. Therefore, the variogram is defined as: 2 (6.25)

where N (h) represents the number of points separated by a vector h in S. The variogram can be viewed as the estimation of the variance of the resulting estimation prediction error when the mea-148 6.6. Kriging surement in s + h is estimated from the measurement in s. Modeling the variogram during kriging is one of the most important part since the final estimation will depend on this spatial model. In the particular case of a spatially stationnary field, this variogram can be expressed as a spatial covariance [START_REF] Noel | Statistics for Spatio-Temporal Data[END_REF].

Variogram models

The variogram can be deduced from the experimental data based on Eq. (6.

2 ) (6.27)

Matérn

The Matérn covariance function is commonly used for its mathematicals properties, in particular due to its finite differentiability (for finite ν) [START_REF] Michael | Interpolation of Spatial Data: Some Theory for Kriging[END_REF] and the fact that it includes the modified Bessel function of second kind makes it respect physical equations, depending on the choice of ν.

where Γ is the gamma function, and K ν the modified Bessel function of second kind, ρ ∈ R + * , ν ∈ R + *

Simple Kriging

Depending on the knowledge that one can get on its data, or the knowledge that one wants to infere into its model, the kriging method and formula has to be adapted. In Simple Kriging (SK), the measurements y are assumed to be the sum of a known trend function µ : s ∈ S → µ(s) ∈ R and a

This state-space model is updated each time new incoming data are gathered through a Kalman Filter, as explained in Section 6.3.

One can see from Eq. (6.50) that the presented method needs the priors µ ϵ and Σ ϵ . The prior on Σ ϵ can be roughly estimated by the input data assuming that the space covariance emissivity matrix will remain the same through time (as long as the scene is static). However, the prior on µ ϵ is a strong one that will have an impact on the final temperature and emissivity estimation. If reference emissivity points are available and spread in the image, then those references should be taken as locations for which the Kalman filter must be executed due to a fair estimation of µ ϵ .

More generally, to overcome the lack of knowledge on µ ϵ , a possibility is to replace the Kalman filter by an interacting Kalman filter. In that case, the weights of the particles are updated by assuming that for each s's neighbors measurements locations, the temperature is the same as in s.

Furthermore, the A t dynamic evolution of the temperature is also unknown. The A t value is therefore added to the unknowns of the model and added to the components of the particles.

A t represents the temperature's evolution from one time to the next one. As a consequence, A t can be greater than 1. This is actually an issue for the Kalman recursive filter. If A t is greater than 1, the initial temperature value may not be forgotten as the algorithm is computed. To overcome this issue, and for outdoor data process, the algorithm is started during the morning to guarantee that A t < 1 most of the time (the dynamic evolution is in 1 T ).

Note

Spatial data information -Kriging

The Kriging step is used in order to exploit the spatial information by using the assumption that the spatial covariance matrix model is constant through time. The Kriging interpolation is used to retrieve the emissivity values on the entire image. As a consequence, only a few points on the image are tracked through time. This methodology helps to reduce the computation time of the particle filter algorithm on the entire image. By writing the measurement equation with the stationary component Eq. ( 6.52), the emissivity at any new location s ∈ S within the image can be derived from the Kriging Eq. (6.53) (see Section 6.6):

CMA-ES applied to the temperature / emissivity estimation

CMA-ES (covariance matrix adaptation evolution strategy) belongs to stochastic and derivative-free methods, based on an evolution strategy [START_REF] Hansen | The CMA Evolution Strategy: A Comparing Review[END_REF]. It has been used to address the temperature and emissivity separation problem and compare it to the KIKF method. The idea behind an evolution strategy is to choose λ children from µ parents population. In order to produce one child, ρ parents are combined.

Then children are mutated, usually based on a multivariate normal distribution. Finally, the µ best children based on the evaluation of a minimization criteria are selected and become the new parents population.

The CMA-ES algorithm relies on this principle. However, instead of mutating the children with a multivariate normal distribution a priori, the covariance matrix of this distribution is updated.

Theory

As for the particle filter presented previously, the CMA-ES algorithm will have a sampling step, a mutation and then a selection. However, instead of updating only the mean of the search points, the covariance matrix (i.e. the shape of the distribution ellipsoid) is also updated. Let λ be the population size and i ∈ [1, λ]. A sample x i is therefore normally distributed with its mean m and covariance matrix C

Where σ ∈ R + is a step-size that controls the particles evolution.

From Eq. (6.54) different evolution strategies (ES) can be applied to update µ i , Σ i and σ. The (µ/µ, λ)-ES chooses the best µ points from new solutions which is a non-elitistic selection. Then, recombination is applied. Let f be the objective function to be minimized. The solutions points are ranked according to their objective value f (x 1:λ ) ≤ . . . ≤ f (x λ:λ ). The mean is then updated according to their weights:

where N eff has been defined in the particle filter section.

In order to control the algorithm evolution, a path length control quantity p is introduced to update both the covariance matrix and the step-size factor. This is achieved by using a cumulative step-size adaptation (CSA) as described in [START_REF] Alexandre | Cumulative Step-Size Adaptation on Linear Functions: Technical Report[END_REF]:

Spatial optimization problem

As for the KIKF method, the approximation that the temperature is the same on a given pixel's window has been used.

CMA-ES

Then, the CMA-ES described in this chapter has been tested on the simulated data. As explained previously, the temperature is assumed to be constant over a given pixel's neighborhood. Such assumption can be discussed. In fact, at thermal equilibrium, for a short distance and small angle of view this assumption remains valid. However, when important distances are involved, a given pixel's represents a larger area. In that case, two pixels may have very different temperature. The method has been applied over the entire image at one timestamp and the results are presented for the band (10µm-12µm).

Boundaries are applied to the optimization problem: 0 < ϵ < 1 and 260.15K < T < 300.15K.

An initial emissivity constant value is given as a starting point of the optimization algorithm. This first guess shown in Fig. 6.20 leads to important difference between initial temperature and ground truth temperature. The result at level 3 is shown in Fig. 6.21. The emissivity map starts to appear in this coarse preliminary result. However, the temperature estimation suffers from the coarse approximation, due to the assumption that it is constant over a given neighborhood.

The post process of Fig. 6.22 is shown in Fig. 6.23. Two metrics have been used to compare the result given by the CMAES method with the ground truth images. The first one is the root-mean-square error (RMSE) that evaluates the distance between the estimated values and the expected ones. The RMSE is performed locally on a given window of 5 × 5 pixels. The second one is called the structural similarity (SSIM) which aims at measuring the similarity between two images, developed for image reconstruction and compression. The SSIM index is defined between two window u and v as:

where

• σ u is the variance of u,

• cov uv is the covariance of u and v,

• c 1 , c 2 and c 3 are constants that aim to characterize the saturation effects and ensure numerical stability for denominators close to zero. Finally, in order to improve the estimation process, one idea is to use hybrid methods by combining for example the CMA-ES and KIKF. The CMA-ES could provide a first estimate to the KIKF method that would track dynamically the evolution of the system's parameters. When the KIKF tracking is made difficult, the CMA-ES could be called again.

Method

Appendices YAML example configuration file

An example of a YAML configuration file used in the infrared multi-spectral images simulator is given. Appendix C

Simulator Python script example

Temperature auto-regressive model

In this appendix, the dynamic model that will describe the temperature's evolution through time is discussed. Since the temperature evolution through time has to be predicted in the Kalman filter, a model for this prediction is needed. The evolution of daily temperature follows some patterns that can be exploited.

D.1 Seasonal ARIMA model

In • The variance γ is the same at all time points

• The covariance between any two time points t, tk depends only on k and not on the value of t

The definition of a weakly stationary process is important for the ACF to make sense and therefore for an ARIMA model. An ARIMA model is defined by the combination of an autoregressive model (AR), a moving average one (MA) and a differencing part (I).

An autoregressive model of order p, denoted AR(p) is given by: Definition D.1.2 (AR(p)).

with (Φ i ) 1≤i≤p the model parameters, c a constant and w t ∼ N (0, σ 2 w )

A moving average model of order q, denoted M A(q) is a weighted past error, given by: Definition D.1.3 (MA(q)).

X t = µ + w t + θ 1 w t-1 + θ 2 w t-2 + ... + +θ q w t-q (D.2)

201

Let defines the lag operator, L which will be used in the ARIMA model definition.

Definition D.1.4 (Lag operator). Given the time serie X t , t ∈ Z, the lag operator L produces the previous element of a time serie:

An ARM A(p, q) model can be defined as:

By using the lag operator L:

In the case that the polynomial

has a unit root of multiplicity d, this defines an ARIMA(p', d, q) model:

Please note that an ARIMA(0,1,0) leads to a random walk. Moreover, seasonality can be added to the ARIMA model when high values or low values tend always occur in particular pattern. In that case, a seasonal ARIMA model (SARIMA) can be added for which an ARIMA model is multiplied by the seasonality lag:

Definition D.1.5 (Seasonal ARIMA model). In a SARIMA model, seasonal and non-seasonal components are mixed up. Let S be the seasonality (S = 12 for monthly seasonality for example). We define the SARIMA model such as: (p, d, q) × (P, D, Q) S where (p, d, q) and (P, D, Q) are the non-seasonal and seasonal components respectively.

By removing the seasonal component of the signal, the SARIMA model enables the use of lower polynomial orders than an ARIMA equivalent would have used in order to improve the computation performances.

D.2 SARIMA model for daily surface temperature

The SARIMA model have been fitted to actual measurements by using the astsa (Applied Statistical Time Series Analysis) R package [START_REF] Shumway | Time Series: A Data Analysis Approach Using R, English, 1 edition[END_REF]. A daily seasonal model has been chosen to take into account the