Vertical Cracks Excited in Lock-in Vibrothermography Experiments: Identification of Open and Inhomogeneous Heat Fluxes

  Citation: Mendioroz, A.; Castelo, A.; Celorrio, R.; Salazar, A. Vertical Cracks Excited in Lock-in Vibrothermography Experiments: Identification of Open and Inhomogeneous Heat Fluxes. Sensors 2022,22, 2336. https://doi.org/ 10.3390/s22062336 Academic Editor: Giacomo Oliveri Received: 25 February 2022 Accepted: 16 March 2022 Published: 17 March 2022 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). sensors Article Vertical Cracks Excited in Lock-in Vibrothermography Experiments: Identification of Open and Inhomogeneous Heat Fluxes Arantza Mendioroz 1,*, Alazne Castelo 1, Ricardo Celorrio 2and Agustín Salazar 1 1Departamento de Física Aplicada, Escuela de Ingeniería de Bilbao, Universidad del País Vasco UPV/EHU, Plaza Ingeniero Torres Quevedo 1, 48013 Bilbao, Spain; [email protected] (A.C.); [email protected] (A.S.) 2 Departamento de Matemática Aplicada, EINA/IUMA, Universidad de Zaragoza, Campus Río Ebro, Edificio Torres Quevedo, 50018 Zaragoza, Spain; [email protected] *Correspondence: arantza.mendior[email protected] Abstract: Lock-in vibrothermography has proven to be very useful to characterizing kissing cracks producing ideal, homogeneous, and compact heat sources. Here, we approach real situations by addressing the characterization of non-compact (strip-shaped) heat sources produced by open cracks and inhomogeneous fluxes. We propose combining lock-in vibrothermography data at several modulation frequencies in order to gather penetration and precision data. The approach consists in inverting surface temperature amplitude and phase data by means of a least-squares minimization algorithm without previous knowledge of the geometry of the heat source, only assuming knowledge of the vertical plane where it is confined. We propose a methodology to solve this ill-posed inverse problem by including in the objective function penalty terms based on the expected properties of the solution. These terms are described in a comprehensive and intuitive manner. Inversions of synthetic data show that the geometry of non-compact heat sources is identified correctly and that the contours are rounded due to the penalization. Inhomogeneous smoothly varying fluxes are also qualitatively retrieved, but steep variations of the flux are hard to recover. These findings are confirmed by inversions of experimental data taken on calibrated samples. The proposed methodology is capable of identifying heat sources generated in lock-in vibrothermography experiments. Keywords: crack characterization; lock-in vibrothermography; ultrasound-excited thermography; sonic-infrared; inverse problems; nondestructive testing 1. Introduction Thermographic non-destructive testing (NDT) methods have demonstrated a high potential for surface and subsurface defect detection and characterization [ 1 ]. Thermographic techniques consist in generating a thermal unbalance in the material and recording the evolution of the surface temperature distribution by means of an infrared camera. The thermal perturbation can be carried out by exciting the material with light (optically excited infrared thermography (IRT)), ultrasounds (vibrothermography, thermosonics, sonic IR), or electromagnetically (inductive thermography). The most popular modality of infrared thermography uses light to heat the material surface. The presence of defects perturbs the subsequent heat diffusion, giving rise to anomalies in the surface temperature distribution with respect to a sound material. Consequently, the signature of the defect needs to be identified in a pre-existent temperature field caused by the excitation. In this regard, vibrothermography has attracted a great deal of interest in recent times due to its defect-selective nature. In vibrothermography, the material is excited with high-amplitude ultrasounds. In non-viscoelastic materials, the bulk dissipation is small and the mechanical energy is converted into heat at cracks, mainly due to friction between the crack lips. This Sensors 2022,22, 2336. https://doi.org/10.3390/s22062336 https://www.mdpi.com/journal/sensors Sensors 2022,22, 2336 2 of 20 thermal energy diffuses in the material and eventually reaches the sample surface, producing a hot region above the defect, in a cold environment. For a comprehensive description of vibrothermography, see [ 2 ]. If compared to optically excited thermography, there is a double advantage in generating an internal heat source at the defect. First, the resulting surface temperature distribution is background-free and only due to the heat generated at the defect. Second, the heat travels only one way to reach the surface, which allows sensing deeper regions in the material. These advantages apply to any NDT method generating heat at defects, for instance, the identification of metallic inclusions embedded in an electrical insulator when the parts are excited by eddy currents (inductive thermography). The heat generated at cracks in vibrothermography is generally non-uniform. Actually, in surface-breaking open cracks, the region where the crack lips are not in contact does not produce heat (unless an induced breathing mode brings the two surfaces into contact [ 3 ]) and, close to the crack border, the closure stresses might lock the crack asperities, thus preventing heat production [ 4 ]. Intermediate regions where the lips are in contact and in relative motion produce heat, generally with a non-uniform distribution. Accordingly, the geometry of this flux distribution is the information accessible from temperature data measured at the surface, rather than the crack geometry. The identification of the shape of internal heat sources from surface temperature data is a severely ill-posed inverse problem due to the diffusive nature of heat propagation. The strategies to solve this problem in a general form can be roughly categorized into leastsquares minimization methods, statistical methods, and the new “virtual wave concept” method. In least-squares minimization, the cause of the observed temperature distribution (here, the inner heat source distribution) is identified by minimizing the squared L2-norm of the difference between the data and the prediction of the model (residual). The ill-posed character of the inverse problem makes this minimization unstable, and in order to find a sensible solution, the inversion needs to be regularized. An efficient strategy to stabilize the inversion and incorporate information on the characteristics of the solution consists in adding one or several terms to the residual that provide stability to the minimization. The minimization can be carried out by either global methods (neural networks [ 5 ], genetic algorithms [ 6 ], particle swarm optimization [ 7 ]), which search for the solution over large ranges of parameter values, or local methods (Gauss-type or conjugate gradient [ 8 ]), which modify the starting parameter values in a controlled way. Global methods are aimed at finding the rough global minimum but are less precise in finding the optimum solution and entail a high computational cost, whereas local methods may find the minimum more precisely but risk getting trapped at local minima. In statistical methods [ 9 – 11 ], the solution is characterized by featuring the highest probability from a statistical point of view. Knowledge of the statistical uncertainty of the data set is required, as well as having a forward model in order to calculate the probability distribution to find the solution. Lastly, the recently developed virtual wave concept [ 12 ] is configured as a two-step problem. The first problem consists in calculating the so-called virtual wave, which can be understood as the wave equation solution counterpart of the true heat diffusion problem. Once found, in the second step, back-projection techniques allow finding the heat source distribution. The main difference between least-squares minimization and statistical methods versus the virtual wave concept is that the former need a physical model to describe the direct or forward process (calculation of the surface temperature from knowledge of the heat sources), whereas the later does not need modelization of the direct problem. So far, statistical methods and the virtual wave concept have been applied to characterize volumetric heat sources [ 12 , 13 ]. Least-squares minimization approaches have been implemented to characterize ideal, compact, and homogeneous vertical planar heat sources from lock-in vibrothermography data [ 14 , 15 ]. However, the heat generated by real cracks does not follow ideal, compact, and homogeneous distributions, unlike the sources treated in these previous works [ 14 , 15 ]. With the idea of approaching practical situations, in this work, we address the characterization of heat sources typically generated by real Sensors 2022,22, 2336 3 of 20 surface-breaking vertical cracks with half-penny shape, as well as inhomogeneous heat sources in vibrothermography experiments. We confine our study to the thermal diffusion problem, leaving aside the mechanisms that give rise to the heat generation. We focus on amplitude-modulated excitation and lock-in detection, as this modality is aimed at reducing the noise in the data, which is crucial in ill-posed inverse problems. In Section 2, we present the solution of the direct problem for the geometries addressed. In Section 3, we present a comprehensive overview of a regularized least-squares minimization approach, in order to give some insight on the meaning of regularization, and we describe the inversion algorithm. The potential and limitations of Lasso (L1) [ 16 , 17 ] and Total Variation (TV) [ 18 , 19 ] regularizations to identify open and inhomogeneous heat sources is shown in Section 4by inverting synthetic data with added noise. In Section 5, we present the experimental set-up and inversions of experimental data, discussing the results. Finally, in Section 6, we summarize and conclude. 2. Direct Problem The direct problem consists in calculating the surface temperature distribution generated by a certain distribution of modulated heat sources (at frequency f, ω = 2 π f) located in plane Π (x= 0) perpendicular to the sample surface (z= 0). We consider that the sample is semi-infinite in the zdirection and infinite in xand ydirections, with thermal conductivity Kand diffusivity D. The geometry is depicted in Figure 1a. Figure 1. ( a ) Geometry of the problem, with heat sources in red; ( b ) detail of the geometry of the heat source, representing an open half-penny crack; (c) geometry of a rectangular heat source. Neglecting heat losses by convection and radiation, the complex temperature at the surface due to the thermal waves launched at frequency ffrom Ω can be calculated by integrating the contribution of point-like modulated heat sources in plane Π (confined in area Ω) [20]: Tf(x,y, 0) = x Π Q(→ r0) 4πK e−qf|→ r−→ r0| → r−→ r0 dS0=x Ω Q(→ r0) 4πK e−qf|→ r−→ r0| → r−→ r0 dS0(1) where Q(→ r0) is the position-dependent flux amplitude (null outside Ω ) and qf=p2πi f /D is the thermal wave vector. In order to describe the heat produced by half-penny surfacebreaking cracks, we focus on heat sources featuring the shape of semi-circular bands of radii r 1 and r 2 (r 2 >r 1 ). For the sake of generality, we allow the heat source to be slightly buried, with the upper side located at a depth dwith respect to the sample surface (Figure 1b). The complex surface temperature distribution for this case is written as follows: Tf(x,y, 0)= r2 Z r1 π Z 0 Q(r0,ϕ0) 4πK e−qfqx2+(y−r0cos ϕ0)2+(d+r0sin ϕ0)2 qx2+(y−r0cos ϕ0)2+(d+r0sin ϕ0)2r0dr0dϕ0(2) This expression also includes the case of kissing half-penny cracks, by making r1= 0. Sensors 2022,22, 2336 4 of 20 For the sake of comparison, we also present inversions corresponding to other geometries. Just to give an example, we deal with rectangular heat sources of width wand height hburied at a depth dbelow the surface (Figure 1c). In this case, the expression of the surface temperature distribution is written as follows: Tf(x,y, 0)= w/2 Z −w/2 −d Z −(d+h) Q(x0,y0) 4πK e−qfqx2+(y−y0)2+z02 qx2+(y−y0)2+z02 dy0dz0(3) In Section 4, we present inversions of synthetic surface temperature data (amplitude and phase) calculated using Equations (2) and (3). For the inversion, we combine data obtained at modulation frequencies f k = 0.05, 0.1, 0.2, 0.4, 0.8, 1.6, 3.2, 6.4, and 12.8 Hz, corresponding to thermal diffusion lengths µf=pD/πf ranging from 0.3 to 5 mm: high frequencies provide sharp details, whereas low frequencies penetrate deeper in the material. 3. Inverse Problem The general solution of the inverse problem consists in finding the heat flux distribution Q(→ r0) in plane Π , responsible for the observed (noisy) surface temperature data Tδ fk , (k= 1, . . . , kmax), δbeing the noise level in the data (L2-norm of the noise). This approach entails that: 1. The heat sources are known to be confined in a plane perpendicular to the surface (prior knowledge). 2. No specific geometry of the heat source is supposed. 3. The thermal properties of the material are known. 4. The shape of the spatial distribution of heat sources is unaffected by the modulation frequency. Accordingly, even if the heat sources are known to be uniform within region Ω , the inversion is not a mere parameter estimation problem (Q,r 1 ,r 2 , and din Equation (1); Q,w, h, and din Equation (2)) but entails meshing plane Π and determining the value of Qat each mesh node. This gives generality to the solution and is of practical interest, because the shape of the heat source is not known beforehand, but increases the difficulty of solving the problem. In this context, the formulation of the inverse problem in a least-squares sense consists in finding the Qdistribution in plane Π that minimizes the L2-norm of the difference between the data and the calculated temperatures at each frequency, summed for all the modulation frequencies fk,k= 1, . . . , kmax: R2= kmax ∑ k=1  Tfk(Q)−Tδ fk   2= kmax ∑ k=1  AfkQfk−Tδ fk   2 2= kmax ∑ k=1  IfkAfkQ−Tδ fk   2 2(4) Here, Afk is the integral operator in Equation (1), and we have introduced a frequencydependent heat source distribution Qfk(r0)=IfkQ(r0) expressing Qfk(r0) as the product of two factors: a normalized heat source distribution, Q(r0) , which, according to assumption 4, is common to all modulation frequencies, and a set of intensities, Ifk , that only depend on the modulation frequency. This allows using different ultrasound amplitudes depending on the modulation frequency (typically, higher amplitude at high frequency, for which the signal is weaker). In this framework, in the inversion, the temperatures are not calculated using Equations (2) or (3) (or the corresponding expression for a particular geometry) but are obtained as the superposition of the point-like contributions of each mesh node in plane Π (Equation (1)). Accordingly, the number of unknowns in the inversion is significantly high (number of mesh nodes in plane Π ). Given the ill-posed character of the inverse problem, the minimization of R 2 is very unstable, and solving the problem requires stabilizing the Sensors 2022,22, 2336 5 of 20 inversion. A very popular method to stabilize ill-posed inverse problems is truncated singular-value decomposition (SVD). We opt for a different solution, which consists in minimizing a modified version of R2by adding stabilizing terms to the right hand side of Equation (4), because this strategy allows introducing in the inversion prior information about the solution. In the next sub-section, following [ 8 ], we present a comprehensive and progressive introduction to the penalty terms that we incorporate in our inversion, taking truncated SVD as the starting point: from the well-known zero-order Tikhonov to more sophisticated functionals such as Lasso (L1-norm) or Total Variation (TV). We start with a quite general formulation, and later on, we particularize for the problem we are addressing. We have prioritized the smoothness of an intuitive description over rigor in formalism and notation. 3.1. Regularization Functionals 3.1.1. Truncated Singular-Value Decomposition We start by writing the direct problem in an operator form: AQ = T (5) where A is a linear matrix operator that maps the discretized heat source distribution Q in plane Π into the surface temperature data T . The least-squares problem is written as follows: A∗AQ = A∗T(6) where A∗stands for a complex conjugate of A. The solution is: Q = (A∗A)−1A∗T(7) If A has full column rank, ( A*A ) −1 exists, but if it was rank-deficient, ( A*A ) −1 would not exist and Q could not be calculated using Equation (7). The SVD method allows solving Equation (6) for rank-deficient matrices. Just as a reminder, in SVD, matrix A (mby n) is factored into 3 matrices: A = USV∗(8) where U is an mby mmatrix whose columns are orthogonal vectors spanning the data space, V is an nby nmatrix whose columns are orthogonal vectors spanning the model space, and S is an mby ndiagonal matrix whose diagonal elements s i (singular values) are arranged in decreasing order. If only the first psingular values are non-zero (p<m), S can written as S=Sp0 0 0 (9) and Equation (8) can be simplified to A = UpSpVp∗ , where Up and Vp denote the matrices whose columns are the first pcolumns of U and V , respectively. The SVD can be used to compute a generalized inverse of A, the so-called Moore–Penrose pseudoinverse, A†, A†=(A∗A)−1A∗= VpS−1 pU∗ p(10) which always exists. The pseudoinverse solution is then: Q†= A†T = VpS−1 pU∗ pT(11) In an explicit form, the pseudo-inverse is written as follows: Q†= p ∑ i=1 U∗ .,iT si V.,i(12) where U.,iand V.,irepresent each of the pcolumns of Upand Vp, respectively. Sensors 2022,22, 2336 6 of 20 Equation (12) presents the solution as a linear combination of model space vectors, multiplied by factors containing the corresponding singular value s i at the denominator. The summation may include terms with very small singular values that give rise to very large coefficients for the corresponding high-frequency model space vectors V.,i , which may eventually dominate the solution, acting as noise amplifiers. A natural way to stabilize the solution consists in discarding Equation (12) model space vectors with very small associated singular values. This is so-called truncated SVD regularization. However, this stability comes at the expense of reducing the accuracy of the solution. Therefore, the criterion to discard model space vectors must be a trade-off between stability and accuracy of the solution. 3.1.2. Zero-Order Tikhonov Regularization A successful solution of an inverse problem generally involves reformulation as an approximate well-posed problem. The zero-order Tikhonov regularization [ 21 ] modifies the least-squares equation by adding a smoothing term in order to reduce the unstable effects of noise in the data. When the data are noisy, there might be many solutions that adequately fit the data, so that || AQ −T || 2 is small enough. In zero-order Tikhonov regularization, the solutions are sought among those that meet || AQ −T || 2≤δ ( δ being a specific residual misfit value), selecting the one that minimizes the L2-norm of Q: minkQk2, subject to kAQ −Tk2≤δ(13) Introducing zero-order Tikhonov regularization (for a specific regularization parameter αTK), the problem formulated in Equation (13) can be written as the minimization of: R2=kAQ −Tk2 2+αTKkQk2 2(14) Equation (14) is the so-called objective function, and the first and second terms on the right hand side are the so-called discrepancy term and regularization term, respectively. The regularization term is the product of a regularization parameter, αTK , and a regularization functional, kQk2 2 in this case. The larger the αTK , the more powerful the regularization and the larger the error in the solution. We describe our strategy to determine the optimum value of the regularization parameter in Section 3.2.1. The zero-order Tikhonov solution is equivalent to an ordinary least-squares problem augmented according to: QαTK =arg min Q∈Rn   A √αTKIQ−T 0    2 2 =arg min Q∈Rn    AaugQ−T 0    2 2 (15) The size of A remains mby n, and I is the nby nidentity matrix. As long as αTK is nonzero, the last nrows of matrix Aaug are linearly independent, so Equation (15) represents a full-rank least-squares problem that can be solved by its normal equations: A∗ augAaugQαTK =A∗ augT(16) Using the SVD of A and following the steps indicated in Section 3.1.1, the solution can be written as: QαTK = k ∑ i=1 si si2+αTK U∗ ·,iT V·,i(17) where k= min (m,n), and all non-zero singular values and vectors are included. Equation (17) can be rewritten as: QαTK = k ∑ i=1 si2 si2+αTK U∗ .,iT si V.,i= k ∑ i=1 fi U∗ .,iT si V.,i(18) Sensors 2022,22, 2336 7 of 20 where fi=s2 i/(s2 i+αTK) are the so-called filter factors, which control the contribution of the different terms to the sum, in the fashion of a low-pass filter. Comparison of Equations (12) and (18) shows that the penalization of different model space vectors depends on the relation between αTK and their associated singular values. Accordingly, the degree of regularization varies between two limiting cases: for s i >> αTK ,f i≈ 1, and the contribution of the corresponding model space vectors in Equation (18) remains the same as in Equation (12), whereas for s i << αTK ,f i≈ 0, i.e., the associated model space vectors are highly damped. For intermediate singular values, as s i decreases, f i produces a decreasing contribution of the corresponding model space vectors. The result is a filtering of model space vectors with small singular values softer than applying truncated SVD. As a consequence, zero-order Tikhonov regularization produces a smooth solution, since sharp, high-frequency model space vectors are filtered out. Finally, let us mention that it is also possible to apply penalty terms that minimize the L2-norm of the first or second derivatives of the solution, rather than the L2-norm of solution itself. These are the so-called firstand second-order Tikhonov functionals, which are mentioned in the next section. 3.1.3. Lasso and Total Variation Regularizations Focusing now on the particular inverse problem that we are addressing, we come back to Equation (4). As mentioned at the beginning of this section, our goal is to retrieve the vertical heat source distribution Qthat minimizes a regularized version of the squared L2-norm in Equation (4). In practice, this is carried out by meshing plane Πwith nnodes. If a zero-order Tikhonov penalty term is applied, the regularized version of Equation (4) is written as follows: R2= kmax ∑ k=1  IfkAfkQ−Tδ fk   2 2+αTKTK(Q)with TK(Q)=x Π|Q|2dS ≈ n ∑ i=1|Qi|2∆S(19) Zero-order Tikhonov regularization penalizes all nodes in plane Π equally, as it applies the same regularization parameter to each one, with no further information regarding possible locations of the heat sources. However, in order to optimize the degree of regularization, other non-linear regularization procedures based on local information can be implemented, aimed at performing a position-dependent penalization. Lasso (L1) [ 16 , 17 ] and total variation [ 18 , 19 ] regularization methods allow performing a position-dependent penalization by assigning a different regularization parameter to each node in plane Π , which, in turn, is made feasible by implementing iterative methods that make use of the heat source distribution retrieved in a previous iteration. This way, it is possible to have an idea of which nodes need to be penalized more in a following iteration, in order to force some of them to remain damped and keep others dominating the solution. Let us consider a penalty term based on a zero-order Tikhonov functional, as the one considered in Equation (19), but with a regularization parameter that takes into account the solution in a previous iteration: αTKi=αL1 1 QαL1 i,k−1 (20) where idenotes the node in plane Π ,kis the iteration, and we assume that QαL1 i,k−16= 0. The explicit expansion of this new discretized penalty term for all nodes is written as follows: n ∑ i=1 αTKiQαL1 i,k 2∆S=n ∑ i=1 αL1QαL1 i,k 2 QαL1 i,k−1 ∆S= αL1 1 QαL1 1,k−1QαL1 1,k 2+1 QαL1 2,k−1QαL1 2,k 2+. . . +1 QαL1 n,k−1QαL1 n,k 2!∆S. (21) Sensors 2022,22, 2336 8 of 20 As can be seen, despite αL1 being common for all terms, each term is affected by a different penalization, because the QαL1 values are divided by the local values obtained in the previous iteration. In this way, if QαL1 i,k−1 is small and thus 1 /QαL1 i,k−1 is large, then QαL1 i,k is forced to remain small. Otherwise, 1 /QαL1 i,k−1 is small and Qδ,αL1 i,k is free to increase or vary. Over iterations, eventually QαL1 i,k−1≈QαL1 i,k , and the penalty term in Equation (21) approaches: n ∑ i=1 αTKi,jQαL1 i,k 2∆S≈αL1QαL1 1,k+QαL1 2,k+. . . +QαL1 n,k∆S(22) which represents the L1-norm of QαL1 multiplied by the regularization parameter αL1 . Thus, penalizing the least-squares minimization with a penalty term based on the lasso (L1) functional: L1(Q)=x Π|Q|dS =kQk1≃lim k→∞x Π |Qk|2 qε+|Qk−1|2dS (23) can be interpreted as performing a position-dependent penalization of zero-order Tikhonov penalization. The presence of a small constant ε in the denominator of Equation (23) is aimed at avoiding computing errors when |Qk−1|≈ 0. Equations (21) and (22) describe the lagged fix-point iterations algorithm that can be used to approximate the non-quadratic L1 penalty term defined in Equation (23). Regularization with a total variation penalty term: TV(Q)=x Π|∇Q|dS =k∇Qk1(24) is based on the same principle as L1, but acting over |∇Q| instead of |Q| . It can be interpreted as the implementation of a first-order Tikhonov functional with a positiondependent regularization parameter. The lasso functional penalizes the L1 norm of the solution, and TV penalizes the L1 norm of the gradient of the solution. In practice, the main difference between L1 and TV for the solution of the inverse problem is that L1 favours sparse solutions in plane Π (compressive sensing effect), whereas TV favours solutions with areas of null derivatives, which yields blocky solutions. The combination of both is appropriate to characterize the confined heat sources representing cracks that we are seeking. Similarly to Equation (22), which approximates the L1-norm of the solution, since TV is a non-quadratic operator, it can be approximated from first-order Tikhonov penalty functional using lagged fix-point iterations: TV(Q)≃lim k→∞x Π |∇Qk|2 qε+|∇Qk−1|2dS =lim k→∞x Π (∂yQk)2+ (∂zQk)2 qε+ (∂yQk−1)2+ (∂zQk−1)2dS (25) Throughout this section, we have seen that particular regularization functionals produce specific types of solutions: zero-order Tikhonov yields smooth solutions, TV generates blocky functions, and lasso produces a compressive sensing effect. This indicates that, in ill-posed inverse problems, given some prior knowledge of the properties of the solution, the mere selection of the penalty functional is a tool to incorporate this prior information in the inversion. According to the previous results, we stabilize our inversion by penalizing the minimization with two functionals based on TV and L1, plus an auxiliary zero-order Tikhonov penalty term. The properties of TV and L1 motivate this selection, as we seek confined heat Sensors 2022,22, 2336 9 of 20 sources produced at cracks in well-defined areas. The regularized version of Equation (4) to be minimized is written as follows: R2 α=kmax ∑ k=1  Iα fkAfkQα−Tδ fk   2 2+αTKTk(Qα)+αL1L1(Qα)+αTV TV(Qα), with α=(αTK,αL1,αTV) (26) 3.2. Inversion Algorithm The regularization parameters αTK , αL1 , and αTV in Equation (26) determine how large the different regularization terms are with respect to the discrepancy term. The degree of regularization can be varied by modifying the values of the regularization parameters: large values increase the stability of the inversion process, in the sense that the solution becomes less sensitive to noise in the data, but this stability comes at the expense of introducing an error in the solution. 3.2.1. Regularization Parameters In order to find the optimum regularization parameters, our choice is to start iterations with rather high initial values, αTK0 , αL10 , and αTV0 , and reduce them in each iteration according to different decay factors: γTK = 0.3, γL1= 0.75, and γTV = 0.75, respectively. The Tikhonov regularization parameter αTK0 decays much faster than αL10 and αTV0 , so the effect of Tikhonov regularization is basically significant in the first iteration (iteration zero). Tikhonov provides smooth solutions, which is beneficial at the beginning of the inversion and guarantees that the first solution does not get dominated by noise, but sharper solutions are then sought. Moreover, L1 and total variation cannot be implemented at the beginning, because they make use of the solution in a previous iteration. Theoretical results [ 21 ] suggest that it is prudent to stop minimization iterations before achieving the noise level δ . Keeping this in mind, in this problem, we have found that stopping iterations when the minimum discrepancy term is found delivers good results. This is a heuristic stopping criterion, which probably works because we are solving a highly overdetermined problem with quite uncorrelated data noise and gives us optimum results for the retrieved normalized heat source distribution. An important aspect that is worth mentioning about the chosen stopping criterion is that there is no over-fitting of the data, i.e., fitting the noise rather than the underlying function. Regarding the optimum values of the decay factors, there is a lack of theoretical results on this subject. Small values decrease the number of iterations needed to reach the solution, but reduction factors below 0.5 may lead to steps in the discrepancy term being too large for the solution to bet retrieved accurately. The initial values of the regularization parameters as well as their decay factors are chosen by performing systematic batteries of inversions until achieving solutions in a reasonable number of iterations, about 20. Next, we describe the iterative process implemented to find the solution. 3.2.2. Iterations For the inversion procedure, we use domain decomposition iterations to retrieve the normalized heat source distribution, Qα , and the set of intensities, Iα fk , in successive iterations, known as non-linear Gauss–Seidel iterations by blocks. It is a local minimization method used in bi-linear problems such as this one. Coming back to our problem, Tδ fk≈AfkhQα fki=Iα fkAfk[Qα]for k=1, . . . , kmax (27) Sensors 2022,22, 2336 16 of 20 Figure 9. ( a ) Picture of the excitation system; ( b ) detail of the sample closed and in contact with the sonotrode for excitation. Figure 10. ( a ) Experimental natural logarithm of amplitude (left) and phase (right) obtained for a sample containing a semi-circular Cu strip of inner radius r 1 = 1.2 mm and outer radius r 2 = 2 mm buried at d= 0.32 mm, obtained at 0.2 Hz; (b) fitted thermograms. The reconstruction obtained by combining data taken in the whole frequency set (0.05 up to 12.8 Hz) is depicted in Figure 11, together with a reconstruction of the same slab buried at d= 0.71 mm and reconstructions corresponding to other open heat source geometries. Sensors 2022,22, 2336 17 of 20 Figure 11. Grey-level representation of the normalized heat source distribution in inversions from experimental data corresponding to ( a ) semicircular Cu bands of inner and outer radii r 1 = 1.2 mm and r 2 = 2 mm, respectively, buried at depths d= 0.32 and 0.71 mm; ( b ) a square Cu band of outer width 2.8 mm, outer height 1.7 mm, and thickness 0.7 mm buried at a depth d= 0.55 mm; ( c ) triangular Cu bands of outer width 3.6 mm, outer height 2 mm, and thickness 0.9 mm buried at depths d= 0.36 and 0.71 mm. Real contours depicted in red, and values of the depth of the heat sources and quality factor Fon top and under of each reconstruction, respectively. The results confirm some of the features observed in the inversion of synthetic data. On the one hand, rounded contours dominate the reconstructions, which, as explained in Section 4.1, is due to the presence of a TV term in the regularization penalty. Furthermore, the shadowing effect is visible, due to the stronger contribution of the shallowest heat sources that dominate the reconstruction. Nevertheless, in all geometries, the deeper central areas correctly show the path the bands follow, and all depths are well-recovered. The quality factors are in all cases above the cutoff value of F= 0. Next, we tried to obtain experimental data corresponding to inhomogeneous heat sources. As mentioned above, our samples are intended to produce homogeneous heat sources, so we decided to take data combining in the same experiment two strips with the shape of a quarter of a circle to form a semi-circular band with two homogeneous but different heat fluxes on its two halves. In Figure 12, we present the experimental amplitude and phase thermograms obtained by combining two quarters of circular strips made of stainless steel and W (both 25 µ m thick) with inner and outer radii r 1 = 4.2 mm and r2= 5.1 mm , respectively, buried at a depth of d= 0.16 mm below the surface, at a modulation frequency of 1.6 Hz. Unfortunately, we do not have an independent estimate of the ratios of fluxes generated by the two halves in these combinations. The reconstruction obtained by combining amplitude and phase data in the whole frequency set is depicted in Figure 13a (right), together with two more reconstructions, from data obtained using other combinations of materials: on the left, annealed and hard Cu foils, both 38 µ m thick, and at the center, 25 µ m thick hard Cu and stainless steel foils. In Figure 13b, we display the reconstructions obtained for the same material combinations but with triangular geometries. As may be noted, for either geometry, similar results are obtained regarding the heat flux generated by each material combination: the annealed and hard Cu halves (left) act as a homogeneous heat source, whereas differences in the retrieved heat source distribution are more significant for the other two material combinations: Cu–stainless steel (center) and stainless steel–W (right). These differences in the retrieved fluxes are similar for both geometries, which proves the consistency of the inversions. Although the shadowing effect makes the retrieved areas miss the contribution of the central deeper positions in the deepest cases, the overall geometry and the depths of all heat sources are well-recovered. These results prove that differences in the heat flux distributions can be qualitatively characterized with the proposed algorithm. Sensors 2022,22, 2336 18 of 20 Figure 12. ( a ) Experimental natural logarithm of amplitude (left) and phase (right) obtained at a modulation frequency of 1.6 Hz in a sample containing two quarters of circular strips made of stainless steel and W (both 25 µ m thick) with inner and outer radii r 1 = 4.2 mm and r 2 = 5.1 mm, respectively, buried at a depth of d= 0.16 mm below the surface; (b) fitted thermograms. Figure 13. Grey-level representation of the normalized heat source distribution of inversions from experimental data corresponding to ( a ) two quarters of circular bands of inner radius r 1 = 4.2 mm and outer radius r 2 = 5.1 mm buried at depths d= 0.2, 0.27, and 0.16 mm and ( b ) two halves of triangular bands of outer width 5.6 mm, outer height 2.4 mm, and thickness 0.9 mm, buried at depths 0.35 and 0.36 mm. For both geometries, the material combinations for the left and right halves of the bands are the following: 38 µ m thick annealed Cu and hard Cu foils (left), 25 µ m thick Cu and stainless steel foils (centre), and 25 µ m thick stainless steel and W foils (right). Real contours depicted in red and values of the depth of the heat sources on top of each reconstruction. Sensors 2022,22, 2336 19 of 20 6. Summary and Conclusions In this work, we have demonstrated that multi-frequency lock-in vibrothermography data in combination with a least-squares minimization algorithm regularized by TV and lasso functionals allows characterizing ”hollow” non-compact vertical heat sources typically generated by real open cracks in vibrothermography experiments. We have obtained semi-analytical expressions of the surface temperature distribution generated by vertical heat sources with the shape of semi-circular stripes, representing the behavior of open half-penny cracks excited with ultrasounds. A detailed description of the regularization strategies (starting from truncated SVD to Tikhonov, total variation, and lasso) as well as of the inversion algorithm has been presented, and we have proposed a criterion to evaluate the quality of the reconstructions. The inversions of synthetic data with added noise show that the algorithm is able to identify “hollow” uniform heat fluxes and reveal that when the heat source spans a large range of depths, the reconstructions are affected by the shadowing effect, which blurs the deepest part of the heat source, due to the stronger contribution of shallow locations. Inhomogeneities in the heat flux are qualitatively identified except in the case of radial dependence of the flux. The predictions of the reconstructions with synthetic data were confirmed by inversions of experimental data taken on calibrated samples. The results confirm that it is possible to characterize the shape of heat sources generated by open cracks is lock-in vibrothermography experiments. The lock-in processing of modulated data allows detecting signals below the NETD of the camera. The possibility of identifying the regions of the crack that produce heat and the distribution of these heat sources in lock-in vibrothermography open the way to understanding the configuration and dynamics of cracks in this kind of experiment. Author Contributions: Conceptualization, A.S. and A.M.; methodology, R.C. and A.M.; software, R.C.; validation, A.C.; formal analysis, A.M. and A.C.; investigation, A.S.; resources, A.S.; data curation, A.C.; writing—original draft preparation, A.M.; writing—review and editing, A.S.; visualization, A.C.; supervision, A.M.; project administration, A.S.; funding acquisition, A.M. and A.S. All authors have read and agreed to the published version of the manuscript. Funding: This research is part of a project with grant number PID2019-104347RB-I00 funded by MCIN/AEI/10.13039/501100011033. The research was also funded by Universidad del País Vasco, UPV/EHU, grant number GIU19/058. Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. Data Availability Statement: The data are available under reasonable request to the corresponding author. Acknowledgments: The authors are thankful for technical and human support provided by SGIker Computing Services (UPV/EHU/ ERDF, EU). Conflicts of Interest: The authors declare no conflict of interest. References 1. Maldague, X.P.V. Nondestructive Evaluation of Materials by Infrared Thermography; John Wiley & Sons: Hoboken, NJ, USA, 2001. 2. Mendioroz, A.; Celorrio, R.; Salazar, A. Ultrasound excited thermography: An efficient tool for the characterization of vertical cracks. Meas. Sci. Technol. 2017,28, 112001. [CrossRef] 3. Rothenfusser, M.; Homma, C. Acoustic thermography: Vibrational modes of cracks and the mechanism of heat generation. AIP Conf. Proc. 2005,760, 624–631. 4. Renshaw, J.; Chen, J.C.; Holland, S.D.; Thompson, R.B. The sources of heat generation in vibrothermography. NDT & E Int. 2011 , 44, 736–739. 5. Glorieux, C.; Li Voti, R.; Thoen, J.; Bertolotti, M.; Sibilia, C. Photothermal depth profiling: Analysis of reconstruction errors. Inverse Probl. 1999,15, 1149–1163. [CrossRef] 6. Li Voti, R.; Sibilia, C.; Bertolotti, M. Photothermal depth profiling by thermal wave backscattering and genetic algorithms. Int. J. Thermophys. 2005,26, 1833–1848. [CrossRef] Sensors 2022,22, 2336 20 of 20 7. Chen, Z.-J.; Zhang, S.-Y. Thermal Depth Profiling Reconstruction by Multilayer Thermal Quadrupole Modeling and Particle Swarm Optimization. Chin. Phys. Lett. 2010,27, 026502. 8. Aster, R.C.; Borchers, B.; Thurber, C.H. Parameter Estimation and Inverse Problems; Elsevier: Amsterdam, The Netherlands; Academic Press: Cambridge, MA, USA, 2013. 9. Franklin, J.N. Well-posed stochastic extensions of ill-posed linear problems. J. Math. Anal. Appl. 1970,31, 682–716. [CrossRef] 10. Orlande, H.R.; Colaço, M.J.; Dulikravich, G.S. Approximation of the likelihood function in the bayesian technique for the solution of inverse problems. Inv. Probl. Sci. Eng. 2008,16, 677–692. [CrossRef] 11. Stuart, A.M. Inverse problems: A bayesian perspective. Acta Numer. 2010,19, 451–559. [CrossRef] 12. Burgholzer, P.; Thor, M.; Gruber, J.; Mayr, G. Three-dimensional thermographic imaging using a virtual wave concept. J. Appl. Phys. 2017,121, 105102. [CrossRef] 13. Groz, M.M.; Abisset-Chavanne, E.; Meziane, A.; Sommier, A.; Pradere, C. Bayesian inference for 3D volumetric heat sources reconstruction from surfacic IR imaging. Appl. Sci. 2020,10, 1607. [CrossRef] 14. Mendioroz, A.; Castelo, A.; Celorrio, R.; Salazar, A. Characterization and spatial resolution of cracks using lock-in vibrothermography. NDT & E Int. 2014,66, 8–15. 15. Castelo, A.; Mendioroz, A.; Celorrio, R.; Salazar, A. Optimizing the inversion protocol to determine the geometry of vertical cracks from lock-in vibrothermography. J. Nondestr. Eval. 2017,36, 3. [CrossRef] 16. Daubechies, I.; Defrise, M.; De Mol, C. An iterative thresholding algorithm for linear inverse problems with a sparsity constraint. Commun. Pur. Appl. Math. 2004,11, 1413–1457. [CrossRef] 17. Tibshirani, R. Regression shrinkage and selection via the lasso. J. R. Stat. Soc. Ser. B Stat. Methodol. 1996 ,58, 267–288. [CrossRef] 18. Vogel, C.R. Computational Methods for Inverse Problems; SIAM: Philadelphia, PA, USA, 2002. 19. Brune, C.; Sawatzky, A.; Burger, M. Primal and dual Bergman methods with application to optical nanoscopy. Int. J. Comput. Vis. 2010,92, 211–229. [CrossRef] 20. Mendioroz, A.; Castelo, A.; Celorrio, R.; Salazar, A. Characterization of vertical buried defects using lock-in vibrothermography: I. Direct problem. Meas. Sci. Technol. 2013,24, 065601. [CrossRef] 21. Engl, H.W.; Hanke, M.; Neubauer, A. Regularization of Inverse Problems; Kluwer Academic: Dordrecht, The Netherlands, 2000. 22. Breitenstein, O.; Warta, W.; Langenkamp, M. Lock-In Thermography: Basics and Use for Evaluating Electronic Devices and Materials; Series in Advanced Microelectronics; Springer: Berlin/Heidelberg, Germany, 2003.