Three-Dimensional Reconstruction of Thermal Volumetric Sources from Surface Temperature Fields Measured by Infrared Thermography

Non-destructive testing (NDT) of materials and structures is a very important industrial issue in the fields of transport, aeronautics and space as well as in the medical domain. Active infrared thermography is an NDT method that consists of providing an external excitation to cause an elevation of the temperature field in the material, consequently allowing evaluation of the resulting temperature field at the surface. However, thermal exciters that are used (flash lamps, halogen, lasers) act only on the surface of the sample. On the other hand, several energy conversion systems can lead to the generation of volumetric sources; the phenomena of thermo-acoustics, thermo-induction, thermomechanics or thermochemistry can be cited. For instance, ultrasonic waves can generate volumetric heat sources if the material is viscoelastic or if there is a defect. The reconstruction of these sources is the initial process for the quantification of parameters responsible for the heating. Characterizing a heat source means reconstructing its geometry and the supplied power. Identification of volumetric heat sources from surface temperature fields is a mathematically ill-posed problem. The main cause of the issue is the diffusive nature of the temperature. In this work, 3D reconstruction of the volumetric heat sources from the resulting surface temperature field, measured by infrared thermography, is studied. An analysis of the physical problem enables specifying the limits of the reconstruction. In particular, a criterion on the achievable spatial resolution is defined, and a reconstruction limitation for in-depth sources is highlighted.


Introduction
Active infrared thermography is a widely used method in non-destructive testing (NDT). The process involves exciting the environment studied with the help of an external energy supply, and then measuring the resulting temperature on the surface. However, one of the disadvantages of thermal NDT is that the excitation and the measurement are on the surface. Indeed, exciters today used in thermal NDT are usually flash lamps, halogens or lasers, which act only on the surface of the material, making access to information in the volume difficult. Several energy conversion systems can, however, lead to the appearance of volumetric sources in a material: phenomena of thermo-acoustics, thermo-induction, thermomechanics or thermochemistry can be cited. For example, an excitation by ultrasonic waves can lead to heat sources in a material if it is viscoelastic or if it contains a defect. Thus, to characterize this defect or the parameter responsible for the heating, the first step is to characterize the volumetric heat source from the surface temperature data provided by the IR camera. This is a mathematically ill-posed problem, mainly because of the diffusive characteristic of the heat transfer problem as well as the influence of noise. This subject has first been studied as a 2D problem, where the objective was to reconstruct the size and the shape of planar subsurface defects from the blurred temperature measured at the surface [1].
Then, the study has been extended on 3D problems; several pioneering teams have proposed methods based on thermal and acoustic coupling [2][3][4][5][6] or by introducing the concept of virtual waves [7,8]. These methods of reconstruction that exist today enable reconstructing volume data (sources, initial conditions, etc.) from surface data (temperature fields). However, two main limitations emerge from this work: depth and measurement noise. Indeed, whatever the methods used, the deep sources are poorly determined; they are blurred and etched, and their intensities are under-estimated. Finally, measurement noise, unavoidable during experiments is a limiting factor and drawback for the inverse processing method. To overcome this drawback, regularization methods allow approaching the solution, but the choice of the regularization parameters is delicate. The temporal shape of excitation used (modulated, burst, pulse, etc.) does not seem to be a decisive factor for reconstruction: the methods do indeed lead to similar results. Each of these methods has its advantages and disadvantages: methods that use pulsed excitation are, among others, fast but very sensitive to noise, unlike modulated excitation methods that are less sensitive to noise but slower to implement. In this paper, it is proposed to reconstruct 3D heat sources using a method based on an impulse excitation (Dirac type). The impulse thermal response has the advantage of being completely determined analytically.The main objective is to explain, using a physical criterion, the reason for the depth limitation and to propose results on a complete 3D shape generated by the Joule effect.

Direct Problem
From the global point of view, the goal of the proposed method is to be able to reconstruct 3D volumetric sources Q(x, y, z), as illustrated Figure 1, from the measured temperature fields at the sample surface opaque to IR radiation. With this constraint, we will assume that any geometry or shape of sources can be approximated by a sum of variable intensity source points. The impulse thermal response θ Dirac at each point r = (x, y, z) at time t from a heat source located at r = (x , y , z ) delivering, at a previous time t , a pulse of heat quantity of Q(r ) (J·m −3 ) in an infinite homogeneous medium, possibly anisotropic, but of non-thermo-dependent thermophysical properties, is expressed analytically by Equation (1) [9]: where ρ (kg·m −3 ) and C p (J·K −1 ·kg −1 ) are, respectively, density and specific heat of the material, while G 3D (m −3 ) is the 3D-Green function of the material. This function can be written as: where G x , G y and G z are the 1D-Green functions along the three directions of space [9]. If the domain is infinite in all three directions, each 1D-Green function expresses itself explicitly in the same way, here, with spatio-temporal Gaussian functions, and the thermal impulse response becomes: where a x , a y and a z (m 2 ·s −1 ) are the diffusivities along each of the x, y or z directions of space, respectively. If the domain is not infinite in all three directions, but only along the directions x and y, and is only semi-infinite in the third direction (z), the 1D-Green function functions G x and G y are the same as in Equation (3), but G z is different and can be written as [9] : Now, if the domain is infinite along the directions x and y and finite along the direction z, the domain is therefore a plane sample of thickness L z . In this case, the 1D-Green functions G x and G y are the same as in Equation (3), but G z is different and can be written as [9] : For this example of a sample finite along the z direction, the quadripole method is another technique that allows the access to the temperature response [10]. It enables calculating the analytical solution G z of G z in the Laplace domain: The previous Equation (6) can be inverted numerically to find its original in temporal space. The numerical algorithms, such as Gaver-Stehfest, or a return of Laplace by Fourier [10] or using De Hoog's Invlap algorithm, can indeed be used.
For the inversion, several models, including those defined by Equations (3), (4), (5) or (6), can be used. To compare the influence of each model along the z direction, three different configurations of heat sources are tested. As the x and y direction are considered infinite, the heat sources considered on the examples are an infinite plane along these dimensions, located at a given depth in the sample, as illustrated by Figure 2a.
The first example is computed considering this heat source at z = L z /4 (near the measuring surface). On the second example, the source is located at z = L z /2 and for the last one, it is located at z = 3L z /4 (deep source). The resulting temperatures θ(z = 0, t), are calculated by three different models (finite, semi-infinite and infinite along z), measured at the surface, and given by (b), (c) and (d), respectively, in Figure 2. For the three examples, Q 0 = 0.01 J·m −3 , a = 2.5 × 10 −7 m 2 ·s −1 and ρC p = 10 6 J·K −1 ·m −3 . The study is made between t 0 = 0 and t f = L 2 z /a. Whatever the depth of the heat sources, it can be seen that for short times, the temperature responses calculated by the three different models are exactly the same. The differences occur only for the long times.
In this paper, all the works were done at short times. Therefore, the model used is the infinite model, given by Equation (3), because it is more convenient (analytically). However, the choice of the model is not decisive: the methodology that leads to the results is the same, and any models or numerical simulations can be used for the inversion.
By introducing the following dimensionless data t * = t Then for a given shape of a source the resulting temperature will be obtained by: In the case of the 3D problem, the whole material is discretized and each point of discretization is the representation of a voxel considered as a potential source, with a level of energy Q(x 0 , y 0 , z 0 ) = Q 0 Ω(x 0 , y 0 , z 0 ), which corresponds to the mean of the source energy on the defined voxel. The size of the voxels depends on the number of discretizations along each direction. The problem is then written in matrix form, given by Equation (9). The resulting surface temperature at z = 0 is a 2D field as a function of time. For the writing of the problem, each θ(i) is a vector containing the 2D temperature reshaped into a 1D transient temperature: (1) θ (2) . . .
where I is an m by n operator matrix and θ is a vector of size m, and m and n are, respectively, the number of time steps and the depth discretization of the sample. n can then be seen as the number of point heat sources present in the material. For the 1D problem, each point of the operator matrix I described above must be calculated using formula (8) since each of its values is independent of the other. However, for the 3D problem, calculation of matrix I can be optimized as it is composed of several sub-matrices that are repeated (see Figure 3). Thus, it is not necessary to calculate each value of I.x To illustrate this phenomenon, Figure 3 gives the profile of the matrix I, obtained for the case where n x = 3, n y = 4, n z = 5 and n t = 10. In the case of the 3D problem, the operator matrix can indeed be written in the form of a block symmetric Toeplitz matrix [11]: where each element A i is a matrix of size n t n y by n z n y , n z and n y corresponding to the number of discretization points in z and y, and n t to the number of time steps. All the matrices A i are themselves symmetrical in blocks, defined by: where each block A ij is of size n t by n z . The global matrix I is therefore n x n y n t by n z n y n x . Thus, instead of calculating all the points separately, just calculate the first row (by block) of each matrix A i using the formula (7) to reconstruct the entire matrix. The time savings by doing so are considerable. Indeed, for the 3D problem, the number of computations is (n x n y ) 2 n z .

Inverse Problem
The inverse problem given by Equation (9) is written as follows: In practice, I is not square. Its number of lines is much greater than its number of columns. Therefore, I is not invertible. The algorithm used for the studied problem is therefore based on the Moore-Penrose pseudo-inversion method [12,13]. The pseudo-inverse I is obtained by using singular value decomposition (SVD) method [14][15][16][17]. The principle of SVD is to decompose the matrix I of size m by n into a product of three matrices where U is an orthogonal matrix of size m by m, orthogonal V of size n by n and S a "diagonal" matrix of size m by n. The diagonal elements of S are called singular values. They are positive and ranked from the largest to the smallest.
If only the first p singular values are different from zero (p < n), the SVD can then be further reduced. Thus, the pseudo-inversion of Equation (12) gives: The reconstruction problem studied here responds well to the existence and uniqueness requirements of a well-posed inverse problem (cf. Hadamard [18]). Unfortunately, the condition of stability is not respected. One of the methods to study this stability is to observe singular values of operator matrix I: the closer they are to zero, the more the problem is ill-posed. If, in addition, their decay follows a power law such that if s i ≈ i −n , with n strictly positive, the problem is described as ill-posed. On the other hand, if their decay follows an exponential law such as s i ≈ e −i (which is the case here), then it is strongly ill-posed. Then to overcome such problems, it becomes necessary to use regularization methods.

Regularization Methods
Mathematically, in the presence of noise, Equation (12) is no longer valid: strict equality becomes an approximation and inversion can be seen as a minimization problem by the least squares. The unknown Ω is sought in such a way as to minimize the residual R by the introduction of a term of penalization: where α is a scalar called the "coefficient of regularization" and J(Ω) is an independent function of the unknown Ω. The higher the coefficient of regularization is, the more stable the inversion will be. However, in return, the term "penalization" adds a bias to the initial problem. Thus, the greater the term of regularization is, the farther the final solution is from the good one. The difficulty therefore lies in the choice of the adjustment parameters in order to have the best compromise between stability and speed of the solution. Different methods enable choosing the α parameter, for example, the UPRE (unbiased predictive risk estimator), GCV (generalized cross validation), L-Curve and discrepancy principle or Picard methods [17,[19][20][21][22][23]. Then, several types of regularization are possible depending on the choice of J. In this paper, the following types will be studied: (i), regularization by the L 2 standard (Tikhonov regularization [24]), and (ii), regularization by the L 1 standard (Lasso method [25]).
In the Tikhonov case, Equation (13) will take the following form: The comparison between the theoretical solution, given by Equation (13), and the solution of the Tikhonov regularization, given by Equation (15), shows that Tikhonov regularization has the effect of weighting each mode k of the SVD by a coefficient. In the L 1 regularization, Equation (13) will take the following form: where for each iteration i, J (i) is a diagonal matrix that contains the L 1 -norm of the solution of the previous iteration Ω (i−1) . As can be seen in the definition of this matrix, each iteration uses the results from the source Ω (i) from the previous iteration. For the first iteration, it is therefore necessary to define the matrix J (1) . The choice of this initialization does not matter, since from the first iteration, there will be correction according to the data of the inversion. However, it is unadvisable to take the null matrix since that would mean that the first iteration is not regularized. In this study, we take J (1) = Id, the identity matrix. The problem given by Equation (16) is a simple linear equation, given by: The matrix I L 1 is entirely determined since it involves only the data of SVD of the matrix of the basic problem, as well as the solution Ω (i−1) of the previous iteration (so known). Moreover, by its construction, this matrix is square and invertible. The resolution of the problem (17) can therefore be done via any inversion algorithm. Unlike Tikhonov's regularization and truncation of the spectrum of SVD, the regularization by the L 1 -norm allows penalizing each position of the source (node of discretization) with a different parameter. As it is defined, regularization by the L 1 -norm strongly penalizes the "zero" zones (no source), and conversely, it lowly weights non-zero areas (source). This approach favours a solution containing a maximum of null values, which results in strong intensity contrasts on the source reconstruction in contrast to the Tikhonov regularization, which tends to smooth the solution.

Analysis of the Criterion
To introduce the optimal criterion an example of the 1D problem shown in Figure 4, the following equation to calculate the temperature response is presented here: with θ 1D defined by: The heat source repartition Ω consists of three "source points" located at three different depths. In Figure 4, the intensity (normalized) of the nearest source of the surface is 0.5, that of the source of the medium is 1 and the last is 0.75. Here, L z = 0.01 m and Fo z = 0.1. Moreover, no noise will be taken in this part. In this part, the data are synthetic and computed as a consequence of the direct problem defined previously. Therefore, two steps are computed for the algorithm. The direct problem supposes the heat sources are known (location, intensity) and enables acquiring the resulting temperature at the surface (these data are produced by the IR-camera in a real experiment). The second step of the algorithm is the inversion process.
As presented in the previous part, inversion is based on the SVD of the operator matrix I. This matrix is defined for a discretization of every possible source excited by a Dirac excitation at time t = 0. It is actually independent of the real distribution of sources for this given discretization. In this work, it appears that the study of S (singular values matrix) is of prime importance for the performance of the inversion method.
I and S have the same size (m by n), where m is the number of time steps (between t 0 and t f ) and n the number of point of discretizations of the sample. For the computation of the algorithm, these two parameters must be defined. As this inverse problem can be seen as a linear least squares minimization, we must have m > n. Technically, the number of time steps is limited by the acquisition frequency of the IR camera (ac): it cannot exceed M = t f * ac. Therefore, m can be chosen between 1 and M. The discretization of the depth n is not submitted to the same limitations. Two cases of inversion for the problem illustrated by Figure 5 have been tested for two different discretizations of the material: the first inversion has been made with n = 20 (a), the second with n = 50 (b). The number of time steps was taken as m = n t = 300.
For the chosen parameters, if n = 20, the three heat sources are perfectly reconstructed, even the deepest source. However, if n = 50, all sources are poorly found. These results imply that there is a maximal number n of sources (which corresponds to an optimal resolution for the voxel) that we can reconstruct. To determine this limit, the inversion algorithm was carried out one hundred times for n, varying from 1 to 100. For each discretization, the error made between the reconstruction and the real source is calculated, and the obtained results are given in Figure 5d. For n < 30, the errors are negligible, meaning that the reconstructions are faithful to the real sources. However, from n = 30, the error increases rapidly with the number of disctretizations, which means that the reconstructions are erroneous. Thus, for a given problem, there is a maximum number of discretizations n φ to reconstruct the sources with accuracy. This number n φ can be determined by the study of the singular values of the operator matrix I. The SVD of the operator matrix I with a high spatial discretization n = 100 is computed and illustrated in Figure 5c. As can be seen, from a certain value n φ , the singular values are constant and very low in comparison to the first values. As each mode value provides particular information, it means that all modes after n φ do not contribute to the solution of the problem and are therefore non-significant for the reconstruction of heat sources. The comparison between Figure 5c,d enable observing that there is a direct correlation between the significant number of singular values and the optimal number of discretizations along the in-depth direction. This optimal number n φ depends on the parameters of the model, such as the temporal discretization, linked to n t and the Fourier number Fo z . A parametric study has been computed to determine n φ according to these parameters. The obtained results are illustrated in Figure 6a. It appears that the parameter n φ depends strongly on the Fourier number. The higher n φ is, the thinner the size is of the voxel in the depth. Therefore, for Fourier numbers below 0.01, as well as high Fourier numbers, it will be very difficult to reconstruct the sources with a high resolution. At most, it will be possible to reconstruct 10 sources. For a given Fourier number, increasing the number of time steps n t enables increasing the in-depth resolution (with the n φ parameter). However, this increase is made with a logarithmic rate, which means that a large number of temporal data will be necessary in order to gain a few numbers of discretization in depth. This abacus enables defining the optimal parameters (n t , n z , t f ...) in order to correctly reconstruct the heat sources. It is important to note that the criterion given in Figure 6a depends neither on the material nor on the thickness of sample. This criterion depends only on the Fourier number, which is universal, and on the number of time steps, which is numerical, for a Dirac thermal impulse response.
To see the influence of temporal excitation, the same study has been performed with a continuous temporal heat source. The resulting thermal response is simply the temporal integral of the one resulting from Dirac given by Equation (7): The methodology used for the direct and inverse problem of the algorithm is the same as that for the Dirac excitation. However, the operator matrix I is different. The results of the parametric study for Heaviside excitation are presented in Figure 6b. These results have the same tendency as those obtained with a Dirac response.

Analysis of the Influence of Each Mode
The criterion given in Figure 6 enables choosing the optimal discretization in depth (n z = n φ ). Once this parameter is defined, the operator matrix I, given by Equation (10), can be computed. There are therefore n φ modes for the inversion algorithm, where the heat sources can be reconstructed as: with 1 ≤ n r ≤ n φ . If n r = n φ , all the modes of the operator matrix are used and the solution is determined using the sum of each of its modes. The first mode gives a result to which is added the information given by the second mode and so on until the last mode n φ . One can then study the evolution of the reconstruction according to the number of modes n r used for the reconstruction. Figure 7 shows in (a) the global standard error obtained between the reconstruction obtained Ω (n r ) and the real source distribution Ω according to the number of used modes for the example of the 3 sources illustrated by Figure 4. On the same graph, the standard error at the localization of each one of the three sources is also detailed. The sources are numbered from 1 to 3, the first being the source closest to the surface and the third being the deepest one. On the same figure, the reconstructions obtained at different stages of the algorithm are illustrated, corresponding to an increasing number of used modes n r : on (b) n r = 15, (c) n r = 27 and (d) n r = 28 = n φ .
The analysis of the global error enables obtaining an overview of the relevance of the reconstruction. At the beginning of the algorithm, (i.e., with very few modes), the global error is high, reflecting a poor reconstruction of the heat source: the used modes do not provide enough information. The more the number of modes used increases, the more the error decreases, representing a convergence of the algorithm towards a single solution. In addition, the latest mode n φ appears very important because the global error decreases much faster with the addition of the latest modes than with the first modes. The errors calculated for each source, illustrated in red (source 1), yellow (source 2) and violet (source 3), enable visualizing the influence of the modes of the SVD depending on the depth. It appears that the algorithm starts to estimate the sources closest to the surface and finishes with the deepest sources. The error obtained on the source closest to the surface decreases much more rapidly than those of the deepest ones: for n r = 15, the error on source 1 is approximately 10 −2 , approximately one hundred times less than the errors on sources 2 and 3. The reconstruction obtained with n r = 15, given in (b) proves this assessment: source 1 is well reconstructed, which is not the case for the others.
The analysis of graph (a) allows estimating that the source 2, moderately deep, will be fairly well reconstructed with the first 23 modes since the obtained error decreases sharply from n r = 23. Source 3 requires all the modes of the SVD to be reconstructed. Graphs (c) and (d) illustrate the fact that the deepest source requires all modes to be perfectly estimated.

3D Generalisation
For the 3D problem, all Fourier numbers (Fo x , Fo y and Fo z ) have to be considered. Since the x and y directions are symmetrical, the analysis of the problem represented by Figure 8 is sufficient to have access to the 3D criterion. In this study, the heat sources Ω are distributed on the plane π(x = 0), perpendicular to the surface (a). This plane is meshed along the y and z directions, with n y and n z discretization points, respectively (b). As explained in Section 3.1, n z is limited and can be determined via the analysis of the SVD of the matrix operator. In (a) is illustrated the singular values obtained for Fo z = 0.1, Fo y = 1 × 10 −3 and L y = 2L z . Other parameters remain unchanged. The disposition of singular values is peculiar: they are disposed in decreasing steps. Each step contains exactly n y singular values. As in the 1D problem, the last singular values reach a plateau and are then non-significant, and the number of steps corresponds to the maximal number n z of discretizations. This distribution of singular values can be explained by the block repetitiveness in the operator matrix I. The limitation therefore does not come from n y but from n z .
Two principal observations can be made. First, the evolution of the criterion in the function of Fo z is exactly the same as in the 1D problem. Second, contrary to what we could have expected due to separability of Equation (7), the evolution of the criterion with Fo y is not the same as that with Fo z . This evolution is represented in Figure 9b, for Fo z chosen at 0.1. In Figure 9c, the same criterion is represented, but for a number of time steps set at n t = 300, with varying Fo y and Fo z . The 3D reconstruction is optimal for the smallest Fo y and a reasonable number of time steps. The more Fo y increases, the less the discretization can be thin. Physically, a high Fo y means a fast diffusion along the y direction. Therefore, the temperature at the surface due to a buried point source will be largely spread for high Fo y . It can be compared to blurred information: the sources can be estimated with a low discretization. In contrast, for a small Fo y , the surface temperature will be more focused because diffusion along the y direction is almost non-existent. In this case, heat sources can be estimated with a thin discretization.
Thus, optimal reconstruction of 3D buried heat sources is obtained for small Fourier numbers along the x and y directions and a Fourier number between 0.01 and 1 along the z direction.

Application of the Criterion on a 3D Example
To visualize the impact of the criterion, an example is illustrated by Figure 10. A representation of the π-plane vertical to the measuring surface that contains the theoretical heat sources is given in (a) and the resulting temperature at the surface in (b) at a certain time. In this example, all the sources have the same intensity, normalized to 1. The considered dimensions of the sample L x , L y and L z are all normalized in order to have x and y varying between [−0.5, 0.5] and z on the interval [0, 1]. Others parameters are defined such as n x = n y = 20, n z = 25 and n t = 300. According to previous criteria, Fourier numbers are key parameters for the estimation of heat sources. Two cases are tested for this example in order to illustrate this criterion: the first one with Fo x = Fo y = Fo z = 0.1 and the second with Fo x = Fo y = 0.01 and Fo z = 0.1.
The criterion given in Figure 9 enables comparing and verifying the results obtained for both cases. The first of the simulations is tested with Fourier numbers that, according to the criterion illustrated in Figure 9c, do not allow estimating the sources with a discretization higher than 20 in depth. The reconstruction obtained, illustrated by Figure 10d confirms this assertion: even if the sources closest to the surface are reconstructed, it is not the case for the deeper sources. Likewise, still according to this criterion, the Fourier numbers of the second simulation enable estimating the sources with thirty points of discretization in depth, which is confirmed by Figure 10c: all the sources are perfectly reconstructed.

Coupling the Fourier Criterion and the Regularization Parameter
In the previous section, data were synthetic and without noise. In the presence of noise, even with optimized parameters, regularization is necessary. According to Section 2.3, different types of regularization can be made. Figure 11 illustrates the estimation made by the inversion algorithm without regularization (a), with Tikhonov regularization (b) and with L 1 regularization for the same example as depicted in previous section, but with synthetic noised data (SNR = 100). As expected, the solution obtained without regularization (Figure 11a) is completely wrong, while those with the regularization approach the real solution. For each type of regularization, the upper part (the eyes of the smiley) is reconstructed. However, the difference is clear for the deep sources: the mouth of the smiley is completely blurred (smoothed intensities) on the solution obtained with Tikhonov's regularization, while it is much more distinct with the regularization by the L 1 -norm. Regularization by the L 1 norm thus seems to be the favoured in this example, as the heat sources are highly contrasted.

Experimental Results
As depicted in Figure 12, the experiment in this part uses the Joule effect: a Chromel wire (nickel-chrome), traversed by a current transmitted by a generator, constitutes the heat source (a). The principle of the method is to insert this Chromel thread into PVC (polyvinyl chloride) material to form a particular geometrical figure. In this experiment, the diameter of the Chromel wire was 200 µm, and its geometric shape was similar to an "M", placed on the rear face of a PVC sample of 1 mm thickness. To avoid thermal losses on the side of the sample where the wire is located, a thickness of mousse has been bonded, thereby fixing the heating wire (b). The Chromel wire lies on a small area of the PVC sample, as illustrated in (c); therefore it can be considered semi-infinite along the x and y directions. The Chromel wire is bonded to the back of the PVC sample. To estimate the whole 3D effect, the method of the images is used: the sample is then considered to be 2 mm thick with the wire located not on its backside but on the middle of the plane. Thus, at the time of reconstruction, one should find the Chromel wire centred on the middle plane of the sample along the z direction.
The PVC is supposed to be isotropic, with a density of ρ = 1180 kg·m −3 , a heat capacity of C p = 1000 J·kg −1 ·K −1 , and diffusivity a = 1.2 × 10 −7 m 2 ·s −1 . The voltage of the generator was U = 3.5 V with an intensity I = 1.13 A. The measures were carried out by an IR camera, with an acquisition frequency of 200 Hz, and the size of the pixel at the surface was 290 µm × 290 µm. In this experiment, the wire is traversed by the electric current for a duration of approximately 7 s. The area of study on which the temperature is measured as a function of time (with the IR camera) is a rectangle of dimensions 3.65 × 2 cm 2 , divided into 126 × 72 pixels. The resulting field (Digital Level) at the surface is presented in Figure 13 at t = 0 s when the power supply is turned on (a), after 2 s (b), after 7.2 s, where the maximum of the temperature is measured (c), and after 70 s (d). As explained in Section 3.1, the choice of the time interval (of which the Fourier number depends) is very important in order to be able to reconstruct the thermal sources as well as possible. In this experiment, we choose t f = 2 s, allowing to have 400 time steps and a corresponding Fourier number Fo z = 0.2. Since the excitation time of the generator is of the order of 7 s, the temperature response between t = 0 s and t = 2 s is not a Dirac response but a Heaviside response. According to the criterion depicted in Figure 9 the maximum number of discretizations along z for the selected Fourier number is approximately 20. The number of discretization points taken in z (depth direction) is here set to 16 in order to have a number of points close to this limit while remaining low enough to limit the computation time of the inverse problem.
The parameters of the model are then n t = 400, n y = 72, n x = 126 and n z = 16. The size of the operator matrix I to be inverted would therefore be n t n y n x × n z n y n x ≈ 10 7 × 10 5 . To avoid excessive computation time, the domain is divided into 18 equal parts (region of interest) of size 21 × 24 × 16, as illustrated in Figure 14. In each of the eighteen domains, the operator matrix is thus of size 201, 600 × 8064. Each of these domains is treated in a completely independent way. This division is possible because the heat sources are small and the time interval is short (2 s), thus limiting the distribution in the plane. The only problematic case in this division is when the thermal sources are at the boundary between two zones. Edge effects could possibly appear.
In each of the eighteen domains, the solution computed by the program is a cube (3D matrix) of values close to zero at points that do not contain sources and high values (which correspond to the intensity of the source) whenever there is a heat source. The juxtaposition of all the domains enables estimating the totality of the studied volume. For this experiment, the inversion with a regularization L 1 was computed on a dual processor with RAM = 128 GB: the duration of the algorithm was 5 h 45 min. Figure 15 illustrates the obtained results: view along the thickness (a) and along the x and y plane (b).
It can be seen that shape of the heat sources are well estimated: being positioned and centred in the middle of the sample, on a plane situated on the middle plane along the z direction. Then, in a predictable manner, "nodes" can be located at the intersections of each ROI domain. These nodes are artefacts due to the division of the work space and because each domain is treated independently of the others.

Conclusions
The subject of this work is to characterize 3D heat sources buried in a material from the temperature field measured at the surface. An inverse method has been presented in this paper that allows estimating the geometry and the intensity of several volumetric sources. A theoretical analysis of the reconstruction method, based on the thermal response to an impulse excitation, allowed the estimation of a criterion depending on the Fourier number. The resultant abacus of this criterion enables anticipating the quality of the reconstruction of the sources and explaining the limitation related to the depth.
Measurement noise is an important problem for the reconstruction of heat sources because the problem is mathematically ill-posed. The commonly used solution in order to bypass this problem is to add terms of regularization in the method. Two types of regularization were presented in this paper. The abacus depending on the Fourier number remains valid, but the deepest sources will remain strongly impacted by the noise.
An experimental application has been performed to validate the results. The reconstruction method with the help of the abacus and the regularization allows estimating the sources. However, as the problem is three-dimensional, the number of data to inverse is very high, which requires reduction of the study domain. Nevertheless, the execution time of the program is particularly long: several hours are necessary to rebuild sources in a volume of a few cubic centimetres with a resolution of a few hundred microns. To overcome this problem, it is possible to consider other estimation methods, such as probabilistic approach methods.