A joint thermal and electromagnetic diagnostics approach for the inspection of thick walls

In this study, we present a numerical inversion approach to detect and localize inclusions in thick walls under natural solicitations. The approach is based on a preliminary analysis of the surface temperature field evolution with time (for instance acquired by infrared thermography); after, this analysis is improved by taking advantage of a priori information provided by ground penetrating radar reconstructions of the structure under investigation. In this way, it is possible to improve the accuracy of the images achievable with the stand-alone thermal reconstruction method in the case of quasi-periodic natural excitation. 5


Introduction
The integration of geophysical and non-invasive diagnostic methods is a topic of timely interest, and several fields can benefit of this strategy, such as monitoring of the environment, infrastructure protection and cultural heritage management.
In this scientific frame, several efforts have been made and interesting results and strategies are drawn in Eppelbaum (2014) and Alperovich et al. (2013).
Here, we propose to combine two complementary approaches by coupling thermal and electromagnetic reconstructions both based on non-invasive sensing techniques. The non-destructive diagnostic approach proposed here aims at detecting and localizing inclusions in thick walls by exploiting surface thermal and ground-penetrating radar (GPR) surveys. Here, we consider inclusions in the structure, which may represent for example delaminations or inner cavi-ties, whose prompt identification can be crucial in structural health monitoring applications.
The identification of thermal sources (heat flux) or thermal properties (thermal conductivity and/or capacity) of a material has numerous applications. Indeed, a variation of the thermal parameters may be the signature of an inclusion, and the identification process is typically performed by solving an inverse thermal problem. Several reconstruction methods have been proposed to solve the inverse thermal problem, either for source reconstruction (Beck and Blackwell, 1985) or for parameter reconstruction (Ozisik, 2000). In both cases, it is necessary to solve an optimization problem, consisting of the search for the global minimum of a cost function representing the difference between the collected measurements and the model data. Several optimization approaches based on the conjugate gradient or Levenberg-Marquardt algorithm and using the adjoint-state method have been developed in one-dimensional (Ozisik, 2000) or in multidimensional scenarios (Jarny et al., 1991).
In this frame, a method based on the adjoint-state and finite-element method has been developed in recent years at IFSTTAR, with the aim to reconstruct the thermal field over an investigated domain (Nassiopoulos, 2008;Crinière et al., 2014). Afterward, this method has been extended to the reconstruction of the thermal properties of thick walls in two-dimensional (Brouns, 2014;Brouns et al., 2014) or three-dimensional (Nassiopoulos and Bourquin, 2010) geometries. Nevertheless, the results obtained with this method suffer from large, undesired variations of the retrieved parameters close the instrument boundaries due to side effects .
In order to mitigate these spurious effects, in this work we propose a strategy where the reconstruction performance of the thermal method is improved by resorting to the information provided by the electromagnetic approach based on ground-penetrating radar. To this end, the data are collected by long-term thermal monitoring with an adapted infrared thermography system in combination with a GPR.
As is well known, GPR imaging has proved to be useful in several application fields, e.g., the detection of buried pipes, land mines, defects in structures like bridges or roads and the investigation of archaeological sites (Persico, 2014;Persico and Soldovieri, 2006;Matera et al., 2016;Rodríguez-Abad et al., 2016).
The working principle of the GPR is the same as the traditional radar save for the fact that wave propagation occurs in a lossy dielectric medium instead of free space. A wideband electromagnetic pulse is radiated in the investigated medium; owing to the presence of subsurface anomalies, a part of this wave is scattered/reflected, and this backscattered signal is collected by a receiving antenna (Daniels, 2004). Therefore, as for the thermal imaging, GPR imaging relies upon the fact that the electromagnetic signal is capable of penetrating inside opaque materials. Starting from the knowledge of the scattered field, a reconstruction of the electromagnetic properties of the investigated scene can be achieved by solving an electromagnetic inverse scattering problem, which is is nonlinear (Isernia et al., 1996(Isernia et al., , 1999 and ill-posed (Soldovieri et al., 2009).
The diagnostic approach proposed in this work couples both thermal and electromagnetic reconstruction methods and consists of three stages. First a preliminary thermal inversion is carried out to identify possible anomalies in the structure. Then, a GPR reconstruction is performed to refine the preliminary and rough localization of potential inclusions achieved by the thermal method. Finally, the information about the scenario provided by the GPR reconstructed image is exploited in the thermal inverse modeling to improve the accuracy of this method.
The paper is organized as follows. The thermal inverse model is detailed in Sect. 2. Then, Sect. 3 deals with the electromagnetic inverse scattering approach. Section 4 describes the way in which the "geometrical" information provided by GPR is exploited for the thermal data inversion. Numerical examples assess the effectiveness of the proposed inversion approach and are reported in Sect. 5. Conclusions and perspectives are addressed in Sect. 6.

Thermal modeling
The 2-D geometry relevant to the thermal reconstruction problem is sketched in Fig. 1. It is assumed that temperature measurements are collected at the upper surface of the investigated structure, which is a rectangular wall containing several anomalies. The acquisition of the surface tempera- ture time evolution can be carried out by infrared thermography coupled with environmental sensors as presented in Dumoulin et al. (2013) and Dumoulin and Boucher (2014). This outdoor experimental setup permits surveying and recording the pseudo-periodic evolution of the temperature for several days under natural solicitations. In such a measurement field configuration, also called "passive infrared thermography", no external heat source (like lamps or alternative heat excitation sources) are used to tune the system in "active infrared thermography" mode (not easy to set up in outdoor conditions); anyway monitoring of environmental conditions is mandatory for the identification procedure.
A direct thermal model is first established to get numerical temperature data at the surface of the investigated wall. Thereafter, a thermal inversion is performed with the adjointstate method in order to get the reconstruction of the shallower part of the investigated wall. The direct and inverse thermal models are described below.

Thermal state computation
As shown in Fig. 1, the measurements are performed along the upper wall surface m . A heat flux, convective exchanges with the ambient air and radiative exchanges with the sky arise at this surface. On the other wall surfaces 0 , an adiabatic constraint has been imposed. Such an approach is low in time, because heat diffusion propagation and interaction with the defective area inside the inspected thick wall are naturally governed by its thermal properties. On the other hand, an emerging smart solution for autonomous long-term thermal monitoring systems as presented in Dumoulin et al. (2013) and Dumoulin and Boucher (2014) can be set up at a real site to provide measurement data over a long time period, more adapted to depth investigation from surface measurements.
The heat modeling for the investigated wall D is given by Geosci. Instrum. Method. Data Syst., 6, 81-92, 2017 www.geosci-instrum-method-data-syst.net/6/81/2017/ where T (x, t) is the temperature at the generic location x for the time t, ρ is the density, C is the thermal capacity, k is the thermal conductivity, φ is a thermal flux, h is the convective exchange coefficient and h rad is the equivalent radiative exchange coefficient h rad = εσ (T surf + T s )(T 2 surf + T 2 s ) (temperatures here in K). This coefficient is obtained with the linearization of Stefan's law. Moreover, ε is the emissivity of the surface; σ is Stefan's constant; T surf is the temperature at the surface; T s is the sky temperature; T a is the air temperature n is the outward-pointing normal; m is the boundary located at the surface of the domain, where measurements are performed and 0 is the whole of the other boundaries.
We denote with L 2 ( ) the space of square integrable functions on and H 1 ( ) the Sobolev space: (2) The variational and Dirichlet formulations (Gartling and Reddy, 2010) allows one to write, on an orthogonal basis For a numerical solution of the thermal problem, we consider a mesh of elements { e } and nodes {N e }. The approximation T ≈ T is made (Allaire, 2005), with T ∈ V q the discrete space of piecewise polynomials of degree ≤ q in each mesh element e .
with P q being the space of piecewise polynomials of degree ≤ q. A basis of V q is the piecewise polynomials (φ i ) of degree ≤ q such that where δ ij is the Kronecker symbol. Eq.
(3) is still valid for T written on the basis of (φ i ). Therefore, to simplify the notation, T will refer to the approximated solution in the space V q .
From Eq.
(3), we get a differential equation on time for the temperature at the nodes {T }: This equation is temporally solved with an implicit Euler algorithm, so the temperature is obtained at each node of the mesh and for each time step. Finally, in order to provide a more realistic modeling, white Gaussian noise is added to the computed temperature.

Thermal reconstruction
The aim of the inverse thermal problem is to achieve a reconstruction of the thermal capacity ρC and conductivity k in the investigated domain. To achieve this goal, we introduce a functional J to be minimized (Brouns, 2014;Nassiopoulos and Bourquin, 2010). The minimization of this functional is performed with the Levenberg-Marquardt algorithm. This algorithm needs the computation, at each iteration, of the gradient of the locally linearized functional, which is carried out by means of the adjoint-state method. Note that, by using the adjoint-state method, it is possible to get the adjoint equations which have a similar structure to thermal direct Eq. (1). The functional J is expressed as with u ∈ U being the pair of thermal parameters u = {k, ρC} and U the space of unknowns: In the above formula, T (u) is the thermal state, here the temperature field in the space M, associated with u; ε is the Tikhonov parameter; u 0 is an a priori estimation of the thermal parameters u; T m are the measured values obtained with the direct model at the location m for timescale [0, t a ]. Moreover, M is the space of measures M = L 2 ( m , L 2 ([0, t a ])). The norm in M is defined as (Brouns, 2014) where R = L 2 ( ) 2 is the regularization space.
The inverse problem consists in finding u ∈ U as To solve this non-quadratic optimization problem, we apply the Levenberg-Marquardt algorithm, which is frequently used to solve nonlinear inverse problems, as thermal reconstruction (Jarny et al., 1991) or the reconstruction of electric conductivity (Bal et al., 2012). At each iteration, the Levenberg-Marquardt algorithm linearizes locally the functional in order to find a minimum point, which becomes the point around which the linearization is carried out at the next iteration (Levenberg, 1944).
At each iteration, the functional J is linearized in J u next to the value of the reconstructed thermal parameters u at this iteration. For a small variation δu of the thermal parameters, it follows that where ν is the damping parameter, which allows defining a confidence interval around u (Hanke, 2010). The minimum of J u is reached in δu as where J u is defined as To minimize J u , we use the conjugate gradient method, which needs the computation of the differential J u (δu).
The adjoint-state method initially developed in the control theory (Lions, 1971) is also applied. Applying the adjoint method allows keeping the same equation structure and also using the same resolution method for the direct and the inverse problem (Jarny et al., 1991;Brouns, 2014). This method consists in introducing and computing the adjoint operator of T (u). First, we introduce δθ, i.e., the solution of the tangent linear model (Eq. 14): Then, we introduce δθ * , i.e., the adjoint operator of δθ. This operator is the solution of the adjoint problem (Jarny et al., 1991;Brouns, 2014) The adjoint model (Eq. 15) has the same structure as the tangent linear model (Eq. 14) and the direct model (Eq. 1). Therefore, the adjoint and the tangent linear models can also be solved on the same mesh of the direct model with the same finite-element formulation.
The differential of J u can be computed with the optimal control theory by estimating (Brouns, 2014) The conjugate gradient algorithm is applied at each iteration of the Levenberg-Marquardt algorithm, with δ k and δ ρC denoting the infinitesimal variation of these parameters at each step. With J u (δu)δ u, we can get a reconstruction of the thermal parameters k and ρ C in the investigated domain.

Electromagnetic inverse modeling
The GPR imaging of an unknown wall requires recording electric field measurements at the air-wall surface. The GPR is composed of a transmitting and a receiving antenna having a spatial common offset, which are moved simultaneously to gather data along a scanning line. The transmitting antenna radiates an electromagnetic signal in the wall, and a part of this signal is reflected by buried anomalies and detected by the receiving GPR antenna.

Electric field computation
For GPR signal modeling, we consider a 2-D geometry consisting in a three-layered medium where the upper medium is free space, the central medium represents the lossy dielectric wall and the lower medium is free space (see Fig. 1). It is supposed that the upper air-wall interface is located at z = 0, having assumed an x − z coordinate system, with the x axis directed parallel to the air-wall surface and the z axis normal to it. The free space is characterized by dielectric permittivity and magnetic permeability equal to ε 0 and µ 0 , respectively. The wall has a relative dielectric permittivity equal to ε b and an electrical conductivity equal to σ b . All media are supposed to be non-magnetic (Soldovieri et al., 2007a).
In order to simulate the electric field at the measurement points along the air-wall surface, the forward problem has to be solved based on the known geometrical and electromagnetic properties of the scenario. This problem is herein solved numerically by using the popular gprMax tool developed by Giannopoulos (2005). The code discretizes Maxwell equations in space and time by means of the finite-difference time-domain (FDTD) method based on Yee's algorithm. The antennas are modeled as electric line sources assumed invariant along the y axis (Soldovieri et al., 2007a). The perfectly matched layers (PMLs) are used to terminate the computational domain in order to prevent spurious reflections of the Geosci. Instrum. Method. Data Syst., 6, 81-92, 2017 www.geosci-instrum-method-data-syst.net/6/81/2017/ electromagnetic waves from the outer boundaries. At each time step, the electric field measurements are recorded at the receiving antenna location. Moreover, a simulation is carried out for each position of the transmitting and receiving antennas in order to compute GPR data in the multi-bistatic configuration. Finally, in order to account for the effects of measurement noise, white Gaussian noise is added to the computed electric field.

Electromagnetic inverse problem
The goal of GPR imaging is to reconstruct the dielectric properties of the investigated domain D plane starting from the knowledge of the scattered field measurements taken at the surface of the structure under test. The unknown of the inverse scattering problem is described in terms of the electric contrast function χ (Persico et al., 2005;Persico, 2014;Soldovieri et al., 2009Soldovieri et al., , 2007b: where and are the equivalent complex dielectric permittivities of the target and the wall, respectively, and σ b and ε b are the conductivity and the relative dielectric permittivity of the wall. The scattered electric field is related to the unknown contrast function χ by the integral equation (Persico, 2014) where E s is the scattered field datum probed at x s , k s is the wave number in the wall, G e is the external Green's function and E is the total field (summation of the incident and scattered field) in D.
The electromagnetic inverse scattering problem is nonlinear and ill-posed. This fact implies the necessity to apply a regularization scheme in order to obtain a stable solution. Furthermore, the nonlinear nature of the problem brings additional difficulties related to the presence of false solutions (local minima).
In order to avoid the nonlinearity problem and reduce the computation complexity of the related data processing algorithms, the Born approximation is applied to linearize the problem in Eq. (20) and overcome the local minima problem. According to the Born approximation, the total field E in D is approximated as the incident field, i.e., E ≈ E inc . Accordingly, the scattering phenomenon is governed by the linear integral Eq. (20) (Leone and Soldovieri, 2003;Persico et al., 2005;Soldovieri et al., 2007a).
In order to retrieve the unknown contrast function χ, the singular value decomposition (SVD) of the linear operator defined by Eq. (21) is computed. Specifically, the singular system {σ n , u n , v n } is introduced, with {σ n } being the sequence of singular values sorted in decreasing order, {u n } being the basis function in the unknown space and {v n } being the basis functions in the data space. The regularization of the solution is achieved by means of the truncated SVD (TSVD) algorithm (Leone and Soldovieri, 2003;Persico et al., 2005;Soldovieri et al., 2007a).
where , denotes the scalar product in the data space and N is a truncation index fixed usually fixed as a good compromise between accuracy and stability of the solution.

Coupling of thermal and electromagnetic methods
The GPR imaging method presented in Sect. 3 has several drawbacks. First of all, the retrieved contrast function depends on the electrical properties of the background (Soldovieri et al., 2011). Moreover, due to the Born approximation, the reconstructions are only qualitative; i.e., they provide an indication about the position and approximate shape of the targets (Leone and Soldovieri, 2003;Persico et al., 2005). It must be further stressed that some distortions unavoidably arise in GPR images due to multipath clutter (Gennarelli and Soldovieri, 2015). The thermal reconstruction method also exhibits several drawbacks. The effusivity ratio has an important effect on the quality of the reconstructions . The defect area can be overestimated for some kinds of inclusions. The reconstruction is also sensitive to the side effects  and to the noise (Nassiopoulos and Bourquin, 2013).
In this work, we propose to reduce the area of the thermal investigation domain in order to mitigate the side effects and noise sensitivity if there is no inclusion close to the measurement points. This refinement of the investigation domain extent is performed when the GPR image indicates that there is no inclusion near the measurement points.
Accordingly, we first perform a reconstruction of the wall with the GPR. The GPR reconstruction allows defining a smaller investigation region where an inclusion is likely to be located. Thereafter, a thermal reconstruction is performed over this smaller domain.
To get a first location of the inclusions with the GPR, the spatial map defined by the normalized amplitude of the contrast functionχ is considered (Soldovieri et al., 2009).
The normalized contrast function is null where the dielectric permittivity is equal to that of the background and different from 0 over target regions. As introduced and discussed in Soldovieri et al. (2009), a threshold T m is also introduced to get the characteristic function U r . In a first approach, we chose the same value (0.5) for the threshold.
which defines a binary image of the probed scene on the basis of the amplitude of the retrieved contrast function. Based on the information retrieved from GPR images, we define three useful subdomains for the thermal reconstruction (see Fig. 2).
The first subdomain D 1 contains the elements of the thermal mesh located close to the points having a characteristic function equal to 1. This subdomain is most likely to contain the discontinuities or inclusions detected with the GPR.
The second subdomain D 2 contains the elements of the thermal mesh where the characteristic function is equal to 0; inclusions are supposed not to be located in this subdomain.
To reduce the discontinuity between the two below subdomains, we consider a third subdomain D 3 located between D 1 and D 2 . This subdomain aims at being an area in which an inclusion might be retrieved. However, D 3 is different from D 1 insofar as the thermal reconstruction performed in this area takes into account the fact that the probability of an inclusion being in this area is lower than in D 1 .  Let us define A as the set of points (x, z) such that U r (x, z) = 1 and B as the set of mesh element centroids used for the thermal reconstruction.
The subdomain D 1 is made of the elements of centroid (x, z) ∈ B such that The subdomain D 2 is made of the elements of centroid (x, z) ∈ B such that where l 1x , l 2x , l 1z+ , l 1z− , l 2z+ and l 2z− are length such that The lower part of the inclusion might be badly located (Persico et al., 2005). To take into account this effect, we add the following condition: Geosci. Instrum. Method. Data Syst., 6, 81-92, 2017 www.geosci-instrum-method-data-syst.net/6/81/2017/ The subdomain D 3 is composed of the elements of the thermal mesh which are not in D 1 or in D 2 .
The location of these subdomains is shown in Fig. 2. In the case D 3 = ∅, we define in D 3 the norm · D 3 relative to the hermitian scalar product ·, · D 3 defined as We propose to take into account the a priori information coming from the electromagnetic characteristic function U r in the following way (see Fig. 3).
The subdomain D 2 is far from the a priori location of inclusions or discontinuities. The assumption is made that no defect is located in D 2 . Consequently, we suppose that there is no inclusion in this subdomain and that the thermal properties u in this subdomain are those of the background u 0 .
To take into account the fact that the subdomain D 3 is located between D 1 , where the inclusions have been reconstructed with the GPR reconstruction, and D 2 , where we have made the assumption that there is no inclusion, we suppose that for subdomain D 3 the thermal properties are close to those of the background. Anyway, we leave open the option that thermal parameters could not be equal to the background parameters in order to be able to reveal a potential inclusion which would have been located wrong or not identified with the GPR. To do that, we add to the functional J a constraint term: where µ is a scalar allowing the intensity of the constraint to be adjusted. We define U as the subspace of U such that ∀(x, z) ∈ D 2 , u(x, z) = u 0 (x, z). The optimization problem becomes the following: find u ∈ U such that where To minimize J , the differential J u (δu)δ u, introduced in Sect. 2.2, needs to be evaluated. This differential has been evaluated, without a constraint term, in Eq. (16) and becomes This new formulation of the functional J has to be compared to the first one (Eq. 16) expressed in Sect. 2.2, as it integrates now a new constraint term, which allows us to exploit a priori information retrieved from GPR data inversion.

Numerical tests
The numerical examples reported in this study are concerned with a limestone wall containing a wood inclusion (see Fig. 4). The electric and thermal properties of these materials can be found in the pertinent literature (Dumoulin et al., 2010;Giannopoulos, 2005) and are listed in Table 1. The wall has a length of 2 m and a height of 1.5 m, and the wood inclusion is 1 m long and 10 cm thick. The top of this inclusion is located at a depth of 40 cm.
Synthetic GPR data are computed using the gprMax code, by considering a Ricker pulse radiated at the central frequency of 1300 MHz. The data are gathered along the observation line ranging from x = 0.1 m to x = 1.9 m with a step of 3 cm.
The investigation domain for GPR reconstruction has the same size as the wall, and it is discretized with a step of 1.4 cm in both directions. For the GPR reconstruction, we ex- Table 3. Thermal and electric properties of the materials.  ploit a frequency band [100, 2600] MHz, which is discretized with a frequency step equal to 30 MHz. The TSVD algorithm is applied with a threshold on the singular values equal to −30 dB; i.e., the singular values larger than 0.0316 times the maximum singular value are retained.
Then, to get the three investigation subdomains, we use the parameters summarized in Table 2.
With regard to the thermal direct model, we define the following synthetic periodic thermal loading at the monitored surface m : where T a and T s are expressed in units of degrees Celcius ( • C), and rad in watts per square meter (W m −2 ). The measurements are collected over a time interval of 5 days. The   Figure 10. Reconstruction of inclusions with the thermal method with a priori information.
length of this chosen time period is enough in terms of heat diffusion for the inclusion to have an influence on to the surface temperature evolution due to the geometry and materials considered. As the thick wall acts as a low-pass filter in the thermal domain, measurements are not influenced by high-frequency phenomena and in practice different periods can be found for which the weather conditions match quasi-periodic behavior at a natural site. The convective exchange coefficient is equal to h = 10 W m −2 K −1 . The equivalent radiative exchange coefficient is evaluated as h rad = σ ε (T surf + T s ) (T 2 surf + T 2 s ), where T surf is equal to mean air temperature and T s to the mean of sky temperature. White Gaussian noise with an amplitude equal to 0.5 • C is added to the data calculated in accordance with the direct model.
The reconstructed fields for thermal parameters without a priori information (not geometry information inferred by GPR inversion) are shown in Fig. 5a and b.
In order to analyze these images in a simple but quantitative way, we implement the following approach to identify the contour shape of the reconstructed anomaly. To this end, we analyze both the values of the reconstructed thermal parameters and their gradients. An element of the mesh is part of a contour if both reconstructed values, k and ρC, are distant by at least a % of the background parameters. This criterion aims at excluding from the inclusion all the elements having reconstructed parameters close to the background parameters. The threshold a = 10 % is chosen for this study.
We suppose that the thermal parameters have strong variability along the contour of the inclusion and assume that, at the boundaries of the inclusion, the norm of the gradient of both reconstructed parameters is higher than elsewhere in the investigated domain. On the basis of this consideration, we assume that a point (x, z) belongs to the contour of an inclusion only if After processing, closed contours that verified the above condition Eq. (36) are kept, and their contexts constitute a set of subdomains that belongs to the defective area. At that stage, no filtering on the size of the defective area is applied, to preserve multiple detection possibilities.  In order to have a first estimation of the contours, we start with values of the threshold b to 0 and select the points of the investigated domain satisfying both criteria. Then, we increase the value of b until significant variations in the size or the number of the reconstructed inclusions occur. Beyond this critical value b 0 , we observe that at least one main contour which exists with b < b 0 is not closed yet. Here, we have empirically chosen a threshold b = 0.5.
With both criteria, in the case of no a priori information, we get the reconstruction of the inclusion shown in Fig. 6.
Let us turn to consider the case of the thermal reconstruction when the information about the geometry is gained by GPR imaging. The reconstructed (normalized) contrast function is shown in Fig. 7a, and the corresponding geometry of the inclusion is highlighted by the characteristic function in Fig. 7b.
The inferred thermal subdomains are shown in Fig. 8. The thermal reconstruction applied on the subdomains D 1 , D 2 and D 3 lead to the thermal fields presented in Fig. 9a for the conductivity and in Fig. 9b for the capacity.
The reconstruction of the inclusion is shown in Fig. 10; in this case we have empirically chosen a gradient threshold b = 2.2.
To complete our study, we focus on the case in which the inclusion is comprised of air or steel instead of wood. The air is characterized by a relative dielectric permittivity equal to 1, whereas the steel can be assumed as a perfectly electricalconducting object. The thermal properties of these materials are described in Table 3 (Dumoulin et al., 2010;Giannopoulos, 2005).
As can be observed in Figs. 11 to 12, the location and the shape of the inclusion are better identified with a priori information. The area of the reconstructed inclusion is lower with a priori information and closer than the area of the initial inclusion, regardless the nature of the inclusion.
In practice, we focus the minimization of the functional J on two subdomains when we have a priori information. These two subdomains do not include here the part of the initial investigated domain located close to the surface, so the side effects are highly reduced with the GPR information. In particular, the variations of the thermal parameters close to the surface, visible in Fig. 5a and b, do not interfere if we use a priori information. In fact, with the GPR reconstruction, we get the information that there is no inclusion close to the surface, so there is no research of inclusion next to the surface. The noise added to thermal measurements has essentially no Geosci. Instrum. Method. Data Syst., 6, 81-92, 2017 www.geosci-instrum-method-data-syst.net/6/81/2017/ effect on the reconstruction because the high frequency of this noise cannot propagate far from the surface effect (the thermal diffusion from surface acts as a low-pass filter). The coupling of these two methods allows one to greatly reduce the side effects and to get an improvement of the shape of the reconstructed inclusion. However, we can notice that the reconstructed thermal properties' (thermal conductivity and heat capacity) values encompass not only the location of the inclusion but also connected areas in the non-affected part of the material matrix. The relative differences between the reconstructions and the targets are reported in Table 4. The relative errors for the heat capacity are significantly improved by using the joint approach excepted for inclusion filled with air. The relative errors for the thermal conductivity remain important, but care must be taken while making that analysis as the inclusion nature has low thermal conductivity, except for steel. For instance, air thermal conductivity is very low, and small numerical error will significantly affect the relative error. Nonetheless, the reconstructed parameters are coherent: the reconstructed values are higher than those of the background for steel, and lower for wood and air, as was expected. Finally, the inclusion shape reconstruction is enhanced, though this reconstruction phase could take advantage of more advanced post-processing approaches. In particular, we want to point out results obtained for the steel inclusion for which it can be observed that the joint thermal and electromagnetic approach allows us now to better reconstruct the inclusion and not only detect it by a signal processing approach.

Conclusions
In this numerical study, we have proposed a diagnostic approach, for thick wall structures, combining thermal and GPR imaging approaches in order to improve the fidelity of the reconstruction results obtained with only one method.
The thermal method is based on the minimization of a functional, which is achieved with the Levenberg-Marquardt algorithm by computing a gradient with the adjoint-state method.
The GPR imaging problem has been formulated by resorting to a linearized inverse scattering model based on the Born approximation. The resulting linear problem has been inverted via the truncated singular value decomposition scheme.
The combination of both reconstruction methods has then been presented. Based on GPR reconstruction, three subdomains have been identified. The first one is the most likely to contain inclusions; the second one is far from potential inclusions revealed by the GPR. The third one is located among the previous ones. A constraint term has been added to the functional to minimize it in order to take into account the a priori information provided by GPR images.
Numerical examples have been reported and have confirmed improved performance in terms of location of the inclusion, shape of the inclusion, area and thermal parameters. However, the thermal conductivity is still far from the true inclusion conductivity.
Future research work will focus on the validation of the proposed approach against experimental data at a dedicated test site in outdoor conditions.