Bayesian Inference for 3D Volumetric Heat Sources Reconstruction from Surfacic IR Imaging

: The domain of non-destructive testing (NDT) or thermal characterization is currently often done by using contactless methods based on the use of an IR camera to monitor the transient temperature response of a system or sample warmed by using any heat source. Though many techniques use optical excitation (ﬂash lamps, lasers, etc.), some techniques use volumetric sources such as acoustic or induction waves. In this paper, we propose a new inverse processing method, which allows for the estimation of 3D ﬁelds of heat sources from surface temperature measurements. This method should be associated with volumetric heat source generation. To validate the method, a volumetric source was generated by the Joule effect in a homogeneous PVC sample using an electrical thin cylindrical wire molded in the material. The inverse processing allows us to retrieve the depth of the wire and its geometrical shape and size. This tool could be a new procedure for retrieving 3D defects on NDT.


Introduction
The methods for non-destructive testing (NDT) and structural health monitoring (SHM) are numerous, exploiting physical phenomena such as optics, thermodynamics, electromagnetism, acoustics, and thermal effects. However, it is important to emphasize that despite this multidisciplinary nature, there is no optimal or unique method; each method has its advantages and disadvantages. In this work, it is the use of thermal energy as a means of NDT that will be mainly addressed.
Infrared thermography (IRT) is a technique that enables the thermal image of a sample observed in a spectral range of the infrared to be obtained. However, while tomography is now prevalent in various domains, which enables volumetric measurements (x-ray, magnetism resonance imaging, etc.), it is still not possible to measure the temperature in a 3D volume; in other words, it is not possible to realize thermal tomography. The challenge is then to reconstruct 3D heat sources in a volume from the 2D field measured at the surface by thermal imaging cameras. This inverse situation is an ill-posed problem mainly due to the diffusive nature of temperature. To circumvent this problem, classical regularization methods or statistical ones can be used. Interested readers can refer to [1,2], which provides some examples with a rather general overview of the existing inverse methods.
Among the different methods used to tackle heat source reconstruction, the work of Castelo et al. [3][4][5][6][7][8] proposes detecting defects present in a vertical plane on the surface of the material by exciting them with modulated (lock-in) or pulsed (burst) acoustic waves, creating a thermal source at the spot of the defect; the characterization of the defect takes place by the characterization of the resulting thermal source. This reconstruction uses an original regularization term. More recently, a method involving the concept of virtual waves was developed by Burgholzer et al. [9,10] to go back to the initial conditions (heat flow, temperature values, etc.) in a surface temperature field. These works show that the reconstruction methods that exist today lead to the reconstruction of volumetric heat sources (volumetric data) from temperature fields measured at the surface (surface data). However, two main limitations emerge: depth and measurement noise. Indeed, regardless of the methods used, deep sources are difficult to locate; they are blurred and spread out and their intensities are underestimated. The depth problem comes from the physics, so there is no way to overcome that. Finally, measurement noise, inevitable during experiments, is a problematic factor during inversion. Regularization methods enable approaching the solution, but the choice of the parameters of regularization is delicate because it is necessary to find a good compromise between stabilization and accuracy of the solution.
In this work, a reconstruction process based on Bayesian inference, a method largely used in other scientific domains, is proposed to improve the influence of the noise in the inverse processing. In Section 2, the method is applied to a numerical 1D problem. It is then extended to a 3D problem in Section 3. Finally, an experimental validation is performed on a hot wire molded in a plastic sample and warmed by the Joule effect.

Description of the Method through a 1D Example
Let us consider the 1D transient thermal problem depicted in Figure 1. The inverse problem can then now defined as: knowing the temperature evolution θ(t) at z = 0, the question is to determine the position of the source(s) in the thickness, i.e., the vector Ω.
In the paper, the following assumptions are made: (i) The source intensity is considered known; (ii) any geometry or shape of sources can be approximated by a sum of variable intensity source points; (iii) the material is considered orthotropic.
In [11,12], the authors tested different deterministic inversion methods, highlighting some crucial tricks: • For a regular thickness discretization, the number of reconstructable sources is directly linked to the Fourier coefficient Fo z . Indeed, an analysis of the Singular Value Decomposition (SVD) of the matrix A demonstrates that the number of significant modes is limited and linked to the Fourier coefficient. Then, one could not reconstruct more sources than this number.

•
It has also been shown that this limitation is linked to the shape of the temperature evolution at the surface depending on the depth of the source. For the deepest sources, their effect on the thermal evolution on the surface is very similar, and then, the inversion method sensitivity is not good enough to discriminate between them. The reconstruction is then trickier.

•
The reconstruction is also very sensitive to the noise that appears in the temperature measurement at the surface. A Signal to Noise Ratio (SNR) of approximately 50 is sufficient to prevent the inversion method from reconstructing even the closest source to the surface.
Then, to overcome the problem of noise sensitivity, the use of Bayesian inference to reconstruct the source(s) in the thickness is proposed.

Inverse Processing Method Based on Bayesian Inference
In the probabilistic framework, the idea is to reconstruct the repartition of the probability η of finding a source at each position i knowing the temperature evolution at the surface θ (j) , j = 1...m using Bayesian inference based on the Bayes relations [13]: where η(S|θ) is called the posterior-the probability of source repartition S knowing the surface temperature θ-P(Ω) is the prior, P(θ|Ω) is the likelihood-the probability of obtaining θ knowing S-and p(θ) is the model evidence.
In our case, the prior, which allows an a priori knowledge of η to be taken into account, is assumed to be irrelevant; the probability of finding a source is the same for all depths. This could be changed, for example, in the case of defect detection in a composite, where we could assume that the defects are more likely to be present at the layer interfaces.
Therefore, in our case, the posterior can be written as follows [14]: where θ f orward (Ω) is the surface temperature evolution depending on the source repartition Ω. This temperature can come from experimental observations made for different source repartitions, or can be computed for different Ω from a model that can be analytical or numerical. The choice of this model-called forward or direct model as it is representative of the direct problem-brings some remarks: (i) the model used does not need to mimic perfectly the reality as it is used in a probabilistic framework and only the maximum probability is considered but it has to be pertinent enough to bring information on the sources position, (ii) in order to be used efficiently, θ f orward has to be known for all the possible source repartition Ω, then, its computation has to be less time consuming as possible. Equation (2) literally writes: if the sources position is well represented by Ω then θ − θ f orward (Ω) 2 ≈ 0 and then the probability of Ω, P(Ω|θ) ≈ 1; if not, θ − θ f orward (Ω) 2 >> 1 and then P(Ω|θ) ≈ 0. Then, in order to identify the source repartition, one has to identify the Ω maximizing the likelihood in Equation (2).

Direct Problem
In this example, the computation of the forward temperature is done by using the analytical solution of a 1D thermal problem in an infinite media.
The analytical dimensionalized solution of a 1D thermal problem in an infinite media can be written as follows: is the Fourier coefficient. Then, the evolution of the temperature at the surface is given by the following: To compute the temperature evolution at the surface for any position of sources in the material thickness and then compute the likelihood of Equation (2), the problem is discretized in space and time. A source vector Ω is defined by discretizing the material thickness with n − 1 elements, and m time steps are defined between t 0 and t f to sample θ, the measured thermal field at z = 0. Then, the reconstruction problem can be written as in Equation (5): 1D Ω = θ (5) where each point of the matrix A is calculated with Equation (4) [11]. The choice of the thickness discretization requires some remarks: • The size of the discretization is directly correlated with the intensity of the source. Indeed, for a regular spatial discretization of size dz, the intensity of a source at position z is the integration of the intensity of the sources present in the space [z − dz/2; z + dz/2]. Then, the intensity of a source is associated with the chosen discretization; • In the first approximation, the discretization is chosen to be of regular thickness.

Reconstruction of an In-Depth Point Source
To illustrate the use of the Bayesian inference, let us consider the case of a unique point source whose intensity is assumed to be known located at z = z 0 , as illustrated in Figure 2. The surface temperature is built by noising the analytical solution with a white noise ε: The intensity of ε will be modulated to test the sensitivity of the method. The Bayesian inference is used to determine the vector probability P Ω = [P 1 P 2 ...P n ], where P i is the probability of finding a source at position i, knowing the measured surface temperature θ meas . At each time step t j and for each possible thickness position i, the likelihood L ij is computed by the following Equation (7): and a vector P Ω is determined as follows: As Equation (8) shows, the final probability takes into account the information in the whole time interval [0, t f ]. Figure 3 illustrates the evolution of (a) the likelihood L ij and (b) the cumulative probability P Ω (t j ) at each time step j for SNR = 100. One can then note that for early time steps, the deepest sources are still possible (likelihood= 1). This is coherent with the fact that at an early time, only the closest sources can be detected. This explains the fact that for these sources, the probability switches to 0 quickly. With the evolution of time, the deepest sources can be detected.
Then, Figure (4) enables the obtained probability to be determined for three different signal-to-noise ratios: SNR = 10 (a), SNR = 50 (b), and SNR = 100 (c). One can note that the method is not so sensitive to the noise compared to literature methods [8,12,15]: for SNR = 50 and 100, the probability at the source position is superior to 0.95, while for SNR = 10, the probability is 0.30, which can be questionable in absolute terms but not in the case where the other source position probabilities are equal to 0.

Reconstruction of Multiple Sources Using Bayesian Inference
Previously, it was assumed that a unique source was present in the thickness. Then, the likelihood was computed using the response of a unique source at different depths. Now, if n potential sources are assumed in the thickness, the number of possible temperature fields N Ω to be used in the likelihoods calculation is as follows: N Ω = 2 n .
Then, following the same idea described previously, for each likelihood, one would have to compare the measured temperature with the 2 n possible temperature fields. For example, for 25 sources, this corresponds to 2 25 ≈ 10 7 possibilities.
To overcome this problem, Markov chain Monte-Carlo methods are commonly used [14,[16][17][18]. These stochastic methods allow for the simulation of a distribution via a Markov chain: a random source repartition is initially proposed to calculate the likelihood. Among its neighboring repartitions, the most probable one is retained.
Step by step, the method converges toward the solution that corresponds to the data. However, these methods are cumbersome, as the number of tested cases remains high.
To circumvent this, using the superposition principle of the thermal response is proposed, stemming from the linear characteristic of the problem. If θ 1 (t) and θ 2 (t) are, respectively, the resulting temperature field of sources S 1 and S 2 , then the resulting temperature field of S 1 + S 2 is θ(t) = θ 1 (t) + θ 2 (t). Moreover, the sources are detected from the one closest to the surface to the deepest one, allowing for sequential detection.
The principle of the method is then to work by time steps; as soon as the position of a source is located Ω(S 1 ), its analytical contribution is subtracted from the measured temperature field θ mes using Equation (10): The source detection then goes on using the corrected surface temperature θ (t) until the time interval has elapsed and all the sources have been detected. At the end of the algorithm, if all the sources have been correctly detected, the residual temperature R = θ mes − ∑ N sources i=1 AΩ(S i ) must be on the order of magnitude of the noise. Figure 5 illustrates the reconstructed sources obtained with the proposed algorithm starting from the noisy temperature for two different 1D examples. Figure 5a,b shows the results for the first example, whereas Figure 5c,d correspond to the second example. Graphs in Figure 5a,c represent the theoretical sources as well as their reconstructions with the Bayes algorithm. Graphs in Figure 5b,d simultaneously illustrate the theoretical, the noisy, and the residual temperatures at the end of the algorithm. If the algorithm fails in reconstructing a source, as is the case for the second example (Figure 5c), one can see that the residual temperature presents a discrepancy that provides information regarding the error (Figure 5d).
By defining this residue using a floating window, as the time measurement is directly linked to the depth in the material, one can detect the position of the error and correct it in an iterative fashion until reaching a residue that is on the order of magnitude of the noise, as is the case in the graph in Figure 5b.

Generalization of the Bayesian Approach to 3D Source Reconstruction
In this part, the idea is to go from the 1D problem to the 3D one, where the objective is to reconstruct the 3D thermal source(s) from the 2D surface temperature measurement.
To do so, one can reproduce the same scheme as in Section 2.1. In the 3D orthotropic problem, each heat source diffuses in all directions, x, y, and z. The temperature evolution can then be written as follows: with Fo x , Fo y , and Fo z as the Fourier coefficients in each direction. The 3D volume can then be discretized in n x × n y × n z nodes, each one being a possible source whose density represents the density of energy in the elementary volume surrounding it. The surface temperature is also discretized in space and time. Finally, the 3D problem is given by the following: (12) In this case, the classical inversion method to reconstruct the sources is not tractable. Concerning the Bayesian inference, one could reproduce the same scheme as in Part 1 by looking for the probability repartition of sources in the 3D volume by calculating the likelihood P(η(x, y, z)|θ(x, y)). Nevertheless, the size of the problem can be prohibitive when the size of the volume and the surface become larger.
To bypass this complexity, treating the 2D temperature field as n x × n y independent 1D problems is proposed. Indeed, looking at the temperature only above the source allows us to limit the quantity of data to process without loosing much information as almost all the information is contained in it. However, to take into account the 3D component and the interaction between the sources, in this case, once a source is detected, its whole 3D analytical effect given by Equation (11) is subtracted from the measured temperature field θ mes (x, y): As an example, let us consider a problem of 3 point sources p 1 (5, 3, 27), p 2 (11,11,11) and p 3 (18,17,20) placed in a volume discretized in n x = 20, n y = 20, and n z = 30, as illustrated in Figure 6.   Figure 7a represents the 400 measured temperature as function of time, which constitutes the time evolution of the surface temperature T sur f (x, y, t). Figure 7b represents the points where the probability of finding a source is superior to 99%.
As explained earlier, the first sources that are detected are the ones closer to the surface. Therefore, the first source that will be detected for this example should be p 1 , then p 3 , and finally p 2 . In Figure 7b, the closest point to the surface is indeed placed at x = 5, y = 3, and z = 27, where the probability almost reaches 1. This point corresponds exactly to p 1 . We then subtract its effect on the surface temperature evolution, as described by Equation (14): Figure 7c represents the new temperature evolutions on each pixel of the surface. Then, the same procedure is applied until the final time is reached. Figure 7d represents the reconstruction of the second source point p 3 , the newest point that is closest to the surface.
In this model case, the three sources have been reconstructed perfectly. The residual temperature, calculated with the help of Equation (15), is null at the end of the algorithm.

Experimental Validation of a Volumetric Hot Wire Warmed by the Joule Effect
To demonstrate the efficiency of the method, the proposed 3D inverse processing is applied on an experimental case based on volumetric heat source generated by the Joule effect. For that, a chromel wire is molded in a PVC block. The thermal source is created via the Joule effect to warm up the sample; electric current travels on the chromel (nickel chrome) wire over 7.5 s. The temperature was measured by the IR camera on the black-painted rear face of the PVC block, as depicted in Figure 8. In this experiment, the diameter of the chromel wire was 200 µm, and its geometric shape was similar to an "M". It was placed on the rear face of a PVC sample 1 mm thick. To avoid thermal losses on the same side that the sample, an insulating foam was bonded to the surface. The chromel wire was on a small area of the PVC sample; therefore, it can be considered semi-infinite along the x and y directions. The chromel wire was bonded to the back of the PVC sample. To estimate the whole 3D effect, the method of the images was used; the sample was then considered to be 2 mm thick with the wire located not on its back but in the middle of the plane.
The PVC was assumed to be homogeneous and orthotropic with well-known thermal properties given as follows: a density of ρ = 1180 kg.m −3 , a heat capacity of C p = 1000 J.kg −1 .K −1 , and a diffusivity of a = 1.2 × 10 −7 m 2 .s −1 . The voltage of the generator was U = 3.5 V with an intensity of I = 1.13 A. Measurements were carried out by an IR camera with an acquisition frequency of 200 Hz, and the size of the pixel at the surface sample was 290 µm × 290 and µm. Figure 9a shows the thermogram (the mean of the temperature field on all the pixels of the surface as a function of time) as well as the measured temperature at the surface at t = 0 s Figure 9b and at t = 7.5 s Figure 9c. At t = 0 s, i.e., at the very beginning of the experiment, only noise was measured by the IR camera, as shown in Figure 9b. Indeed, the diffusion of the temperature had not yet taken place. During the experiment, the temperature rose throughout the material; at the surface, the geometrical figure formed by the wire is clearly visible. In theory, all the source points formed by the wire are of the same intensity. In addition, since the wire is placed on a plane parallel to the surface, the temperature rise at the line of the wire should be identical in all respects. It can be noted in Figure 9c that this is the case, except in two places where the temperature rise is higher. At the end of the experiment, we can note in Figure 9c that the maximum variation of temperature is approximately 20 degrees.
As explained in the previous section, all the pixels of the surface (containing a temperature as a function of time) are treated simultaneously and considered as independent 1D problems, where the goal is to reconstruct the source in depth using the temperature at the surface point. At each time step, a probability is calculated on each point of the sample using the theoretical model given by Equation (16): where L z (m) is the position at the depth where the thermal source is located, and Q (W) is the intensity (power) tested (W.m −3 ). Equation (16) corresponds to the analytical solution of the thermal response at the level measured at the surface just above a source located at depth L z for a 3D model. Hypothesis: For this study, it is assumed that the intensities of the thermal sources generated by the wire are identical and constant at all points.
Under this hypothesis, the problem can be written in such a way as to search only for the position of the thermal sources (and not their intensities). To do this, we normalize the temperature fields, the measured one (spatio-temporal) and the theoretical one given by the Equation (16) to have an evolution between 0 and 1. The time interval is between t = 0 and t f = 2 s. All temperature values are weighted by the maximum temperature value indicated at time t f . It is important to note that the temperature profiles in each pixel are conserved. The discretization of the sample along the x and y directions is given by the IR camera (pixel of 290 µm × 290 µm), which leads to n x = 126, n y = 72. The discretization along the z direction is defined as n z = 16 [11,12]. The probabilistic approach reconstruction algorithm is then applied, and the solution obtained is illustrated in Figure 10, representing the probability of the presence of the source at each point of the sample. The solution given by the program is a 3D matrix, which represents the volume of the sample, where each point of the matrix contains the probability (>98%) of the presence of a source at this point. The depths of the sources are well found with this method on average; the wire is reconstructed on the plane parallel to the surface, in the center of the plate. Regarding the reconstruction in the x and y directions, the wire is rebuilt thicker than it actually is, especially in areas where the temperature variation is greater. This error is due to the fact that the algorithm considers a constant intensity at any point of the wire. It can also be noted in Figure 10a that in these areas, the wire is reconstructed at a slightly higher position.
Finally, an important point to highlight is the execution time of the program. The calculation, run on a laptop computer (the least powerful of the computers tested) lasted 8.35 s. To compare, in [12], the proposed reconstruction method based on SVD and regularization allows estimating the sources for the same experimental problem. Even if a reduction of the study domain is applied, the execution time of the program is particularly long, more than 6 h are necessary to reconstruct the same sources. The Bayesian approach is therefore definitely faster than a deterministic approach.

Conclusions
In this paper, a methodology of the 3D reconstruction of heat sources generated by any multiphysics wave (electrical, acoustics, optics, etc.) based on a Bayesian approach has been proposed. It constitutes another way to overcome the ill-posed problem of the inversion of a heat equation rather than regularization. The main advantage of this proposed method is the low noise sensitivity and the calculation time compared to other methods. This was validated with an experimental example performed with a volumetric Joule effect heat source. It can be noted that the method can be developed to simultaneously estimate the location and the amplitude of the source and a completely non-destructive test can be performed from such methods. This constitutes the main perspective of this work. Due to the very fast time of calculation of this algorithm, the method can operate on-line and in real time in many industrial systems, such as induction inspection.