Thermographic detection of internal defects using 2D photothermal super resolution reconstruction with sequential laser heating

One of the main factors that govern the resolution limit of thermographic nondestructive testing methods is the diffusive nature of heat conduction in solids. This has an especially severe impact on resolving defect/inhomogeneity structures that lie deep below the surface of the object under test (OuT). As a general rule, using well-established conventional thermographic testing methods, defects can only be fully resolved as long as their lateral extension is as large as they lie deep within the OuT. With the application of photothermal super resolution (SR) reconstruction techniques, it has already been proven that this limit can be overcome.

Super resolution imaging techniques have been successfully applied in several fields where there are resolution barriers to overcome. While geometrical super resolution techniques enhance the spatial resolution of modern detectors to a sub-pixel accuracy,^1,2 optical super resolution techniques allow us to overcome the classical Abbe diffraction limit of optical imaging systems.^3,4 However, currently established photothermal super resolution reconstruction methods within the field of active thermographic materials testing are restricted to the reconstruction of one-dimensional defect structures^5–7 or only approximate two-dimensional super resolution by combining several one-dimensional illumination patterns.⁸ 

Furthermore, in this context, it should be noted that in contrast to the aforementioned methods relying on the spatially structured heating of the OuT, there also exist reconstruction methods working with temporally structured illumination, which allow for a three-dimensional reconstruction of internal structures while also trying to eliminate thermal diffusion effects.^9,10

Within this work, we expand the experimental approach behind laser-based photothermal SR reconstruction toward the detection of two-dimensional defects. This experimental approach is characterized by taking multiple independent measurements for a set of equidistantly distributed positions where the OuT is heated using a single round laser spot with a high power laser. For each of those illuminations, an independent measurement of the temperature response of the OuT is recorded with an infrared camera. In order to then reconstruct the internal defect/inhomogeneity pattern from the acquired set of independent measurements and gaining true two-dimensional information about the defect shape stripped from the blurring-effects of thermal diffusion, two different numerical methods for inverting the underlying severely ill-posed inverse problem are proposed. Subsequently, both methods are experimentally validated on a purpose-made additively manufactured sample with internal defect structures suitable to benchmark the resolution capabilities of each method.

The temporal progression of the front surface (⁠$z=0$⁠) temperature $Tmeas$ of an OuT exposed to a localized external heat flux $Q$ can be modeled as the sum of the initial temperature distribution on the OuT front surface at $t=0s$ and the two-dimensional spatial convolution (denoted as $∗x,y$⁠) of the thermal point spread function (PSF) $ΦPSF$ and a heat source distribution $a$⁠,

The thermal PSF resembles the kernel function in Green’s-function like representation of the thermal diffusion differential equation. The PSF can be calculated analytically for the special case of a defect-free thin plate with thickness $L$ as follows:^6,11

where

• $cp$

  is the specific heat capacity,

• $ρ$

  is the bulk density,

• $α$

  is the thermal diffusivity,

• $ndim$ is the dimensionality of the heat flow:

• point-wise excitation: $ndim:=3$⁠,

• line-wise excitation: $ndim:=2$⁠,

• full-surface excitation: $ndim:=1$⁠,

• $(x¯,y¯)$ is the coordinate centroid of the excitation area,

• $R$

  is the thermal wave reflection coefficient (for typical metals $R≈1$⁠),

• $∗t$

  denotes convolution in time,

• $Q^$

  is the amplitude of the external heat flux $Q$⁠, and

• $It$

  is the temporal structure of $Q$⁠.

For the derivation of Eq. (2), it is assumed that the external heat flux $Q$ can be separated into its temporal structure given by $It$ and its spatial structure $Ix,y∗x,yaext$ according to Eq. (3),

Typically, $It$ with $It(t)∈[0,1]$ is a boxcar function with an active interval from $t∈[0,tpulse]$⁠, while $Ix,y(x,y)∈[0,1]$ encodes the intensity distribution of the external photothermal heating and $aext$ constitutes a distribution of Dirac delta impulses $δi$ with a unit impulse at each excitation position,

While the temporal structure of the external heat flux is already considered within the PSF, the spatial structure is incorporated into the heat source distribution $a$⁠, which is defined as follows:

Here, the additional term $aint$ represents the spatial distribution of the internal “apparent” heat sources. In typical OuTs, even though there are no active internal heat sources present, the fact that regions with internal defects/inhomogeneities show up as higher temperature regions in $Tmeas$ lets them appear in the data as if they were independent heat sources (see Fig. 1). This is due to the fact that they impede the local heat flow at their respective positions. It should be noted that there exist also scenarios in which defective regions show a lower temperature for which $aint(x,y)$ is negative (apparent heat sinks). This is, for example, the case if the OuT includes inhomogeneities with a higher thermal effusivity than its bulk material. Both variants are equally covered by the presented modeling approach.

Therefore, if the heat source distribution $a$ can be reconstructed for the OuT, then, in turn, also the defect/inhomogeneity structure of the OuT is known, which, in general, makes the reconstruction of $a$ from Eq. (1) the overall goal of photothermal SR reconstruction. $aint$ can be described as a Dirac delta pulse distribution as follows:

where $ζi∈[0,1[$ serve as numerical weights in order to incorporate different defect signal strengths. These weights originate from the model assumption that all defects/inhomogeneities viewed as apparent heat sources feature the same PSF as the external heating but depending on their cross section and depth attenuated by the factor $ζi$⁠. Since $ζi$ cannot be easily separated from $aint$ within this modeling approach, by reconstructing the heat source distribution $a$ and, therefore, $aint$⁠, no information about the defect cross section in the depth plane can be gained. $aint$⁠, therefore, only contains the in-plane information on the cross sections of the contained defects/inhomogeneities.

In order to solve the ill-posed problem as stated in Eq. (1) for the heat source distribution $a$⁠, $nm$ measurements with $m∈{1,…,nm}$ are performed. For each measurement $m$⁠, the photothermal heating pattern $aext$ is varied, leading to a system of $nm$ equations to be solved simultaneously,

with $Tdiffm(x,y,t)=Tmeasm(x,y,t)−T0m(x,y)$⁠.

Even though the external component $aext$ from Eq. (5) is in principle known as prior information, it was shown to be beneficial to also simultaneously reconstruct the external component $aext$ besides $aint$⁠, leading to a blind reconstruction. If conversely a non-blind reconstruction approach is chosen, it is critical that the position and intensity profile of the external heating for each measurement is precisely known in order to not introduce reconstruction artifacts. In contrast, for a blind reconstruction, this requirement does not apply, reducing the overall experimental complexity. However, in order to make it possible to separate $aext$ and $aint$ without artifacts in a blind reconstruction context, it must hold true that the sum of the external excitation over all measurements is constant as demanded by Eq. (8),

This condition is based on the fact that it can otherwise not be differentiated if a deviation from the mean is caused by an interference of a defect/inhomogeneity at that location or if it is caused by the non-uniformity of the sum of the external excitation.

Furthermore, to decrease the computational complexity, it is advantageous to eliminate the time dependency of Eq. (7) by only considering a single time step $teval$⁠. When choosing a suitable time step, several factors have to be considered. In order to minimize measurement artifacts from the photothermal heating, a time step during the cooling phase $teval>tpulse$ should be chosen. Due to the exponential nature of the cooling process, an early $teval$ leads to higher SNR. Additionally, the thermal diffusion length of the thermal wave $Ldiff∝t1/2$ increases with time with the thermal diffusivity $α$ of the bulk material acting as the proportionality constant,

Since photothermal SR is only able to gain in-plane information of the contained defects/inhomogeneities, only the planar projection of the defect/inhomogeneity distribution within the testing volume spanned by the region of interest (ROI) and $Ldiff$ can be reconstructed. Therefore, choosing $teval$ is a balance act between achieving sufficient SNR and maximum detection depth for internal defects as is emphasized by Fig. 2.

Currently, well-established laser-based photothermal SR reconstruction methodology is restricted to reconstruction of defects along a single spatial dimension.^6,12 Here, the OuT is heated using line shape laser spots or patterns consisting of several laser lines either by scanning continuously over the ROI or by stepwise scanning. Even though in this way only defects that vary only in a single dimension can be detected, the method has proven to highly increase the resolution limit of thermographic testing. One factor of its success is based on the fact that one of the spatial dimensions can be collapsed by summation leading to a dramatic increase in SNR since all measured data in this collapsed dimension are redundant, which greatly reduces the measurement noise, which is not possible for typical two-dimensional defects.

Other approaches on achieving two-dimensional super resolution have been based on using movable slit-masks in conjunction with photothermal heating with the use of flashlamps. Here, the slit mask needs to not only be moved but also rotated to gain super resolution in the second dimension and the information gain about the second spatial dimension is discretized to the different rotation angles of the slit mask.⁸ In order to expand the method to a two-dimensional reconstruction approach and to be sensitive in both spatial dimensions, a different heating strategy needs to be chosen.

If instead of scanning the OuT with laser lines a sequential point-wise scanning approach is used, where the sample is iteratively scanned on an equally-spaced grid by a single laser spot, 2D-super resolution can be achieved. In this scenario, $Ix,y$ can be modeled as a top-hat round laser spot with a diameter of $dspot$⁠. An equal spacing between all neighboring measurement positions can be achieved by arranging all measurements on a grid consisting of equilateral triangles of side length $rd$ (cf. Fig. 3). A square ROI with area $AROI,2D$ and side length $LROI,1D$ can, therefore, be covered by approximately $nm,2D≈3/2×AROI,2D/rd2$ independent measurements:

In order to achieve sufficient uniformity according to Eq. (8), $rd$ has to be chosen sufficiently small compared to the spatial full width at half maximum (FWHM) extension $dFWHM$ of the PSF, which has a direct effect on the amounts of independent measurements needed to cover the whole ROI. If for a reconstruction along a single dimension $nm,1D=200$ measurements are to be taken, for a 2D-reconstruction, now $nm,2D≈34800$ measurements are necessary to achieve a similar resolution for a fixed spacing of measurements $rd$⁠.

The FWHM spot diameter $dFWHM(t)$ for a PSF with $ndim=2$ is given by the following term considering the standard deviation $σ(t)=2α×t$ of the PSF,

Here, the time dependency can be exploited such that choosing a late-enough $teval$ will lead to a larger $dFWHM$ and, therefore, a larger necessary $rd$⁠, decreasing the amount of measurements to be performed to cover the ROI and increase the viability of the method. However, this comes at the cost of lower SNR and a reduced resolution limit.

The severely ill-posed fundamental equation of photothermal SR reconstruction as stated in Eq. (7) is formulated as a set of convolution problems $Φ∗x,yam=Tdiffm$⁠. Consequently, a suitable deconvolution algorithm needs to be applied to obtain a reconstruction of the internal heat source distribution $arec$ and, therefore, of the inner defect structure of the OuT. For a similar problem in matrix vector product form $Φ⋅am=Tdiffm$⁠, this can be achieved by applying the iterative convex optimization alternating direction method of multipliers (ADMM) algorithm¹³ in conjunction with $ℓ2,1$-norm and $ℓ2$-norm regularizing terms to constrain the solution space,^14,15

where $‖arec‖2,1$ is the $ℓ2,1$-norm defined as

$λ2$ and $λ2,1$ are numerical weights controlling the strength of the regularizing terms in comparison to the least squares minimization term $1/2‖Φ⋅arecm−Tdiffm‖22$ and have to be empirically determined for optimal reconstruction results. Currently, it is part of an ongoing research, how to efficiently select the regularizer weights in the context of photothermal SR reconstruction using machine learning techniques.^16,17

While the $ℓ2$-norm constrains the amplitude of the reconstruction result, the $ℓ2,1$-norm incentivizes the reconstruction of (joint-)sparse defect/inhomogeneity structures, which effect multiple measurements. Since the defect/inhomogeneity structure of the OuT does not change for all measurements and defects/inhomogeneities are typically scarcely distributed within the ROI, this regularization approach greatly increases the reconstruction quality.

Additionally, the ADMM algorithm introduces a penalty parameter $ρ$ (see detailed ADMM-implementations in Algorithms 1 and 2), which balances the effect of the least squares minimization term against the regularizing terms. This parameter can be determined programmatically using the L-curve method.¹⁸ 

Finally, the internal heat source distribution $aint$ widened by $Ix,y$ can be extracted by summing over all measurements where the external component $aext$ is eliminated due to the homogeneity constraint as stated in Eq. (8),

Furthermore, the constant term in Eq. (15) originating from Eq. (8) vanishes due the application of the $ℓ2$-norm within the ADMM algorithm as it will converge to zero over the performed $niter$ iterations.

Within Secs. IV A and IV B, two methods are proposed for transforming the spatial convolution problem into a multiplicative form, such that the aforementioned inversion via the ADMM algorithm can be performed.^19,20

In order to transform the fundamental equation from a convolution problem into a multiplicative form to be able to use the ADMM algorithm for inversion, first the spatial dimensions have to be combined into one single dimension $r∈{1,…,nx⋅ny}$ by vectorizing,

After that, a convolution matrix operator $h(⋅)$ can be introduced. This operator constructs for an input matrix of shape $Φr$ a convolution matrix $h(Φr)∈R2nx⋅ny−1×nx⋅ny$ such that

Here, $Tr0m∈R2nx⋅ny−1$ comprises the measured surface temperature data $Trm$ symmetrically padded with $nx⋅ny/2$ zeros. The convolution matrix $h(Φr)$ constitutes a lower triangular matrix with Toeplitz-structure and is, therefore, quite sparsely populated. Since all $nm$ measurements need to be considered simultaneously for the reconstruction, the $nm$ equations can be stacked on top of each other as follows:

with $h=h(Φr)$⁠, $H∈R(2nx⋅ny−1)⋅nm×nx⋅ny⋅nm$⁠, $A∈Rnx⋅ny⋅nm$⁠, and $TR0∈R(2nx⋅ny−1)⋅nm$⁠. Even though the dimensional size of the matrix $H$ is orders of magnitudes larger than the original data size of $nx⋅ny⋅nm$ data points, it is as a diagonal matrix of sparse component matrices very sparse itself (only $≈1/2nm$ of the entries in $H$ contain a non-zero value), making it possible to be still efficiently stored and handled on modern computer hardware.

In summary, the following minimization problem emerged:

for which the algorithm presented as Algorithm 1 can be applied to solve for $A$ and ultimately reconstruct the internal defect/inhomogeneity map as encoded in $aint$⁠.

The computational complexity of Algorithm 1 is mainly dominated by repeated matrix vector multiplications. Naïvely, this will run in at most $O((2nx⋅ny−1)2⋅nx⋅ny⋅nm3)$ time but due to the high degree of sparsity involved, this establishes only a very conservative upper bound.²¹ 

Another possible way to transform a spatial convolution problem into a multiplicative form is to switch into the spatial frequency domain. This can be achieved by means of a two-dimensional fast Fourier transform (FFT) as follows:

where $fft$ indicates the two-dimensional FFT in $x$ and $y$⁠, $ifft$ its inverse, $ifftshift$ a function that swaps the first quadrant with third and the second with the fourth quadrant of the input matrix, overlines designate taking the complex conjugate, and $⊙$ indicates element-wise multiplication (Hadamard product). In order to get rid of the element-wise multiplication and transform the problem to a matrix vector multiplication problem, the following substitution can be applied:

Since the constraints governing the selection of the regularization terms for the problem inversion, like the sparse distribution of defects/inhomogeneities, are only applicable in the spatial domain, the problem has to be transformed back from the spatial frequency domain in order to apply their effect. This leads to one invocation of the $fft$ and one of the $ifft$ operation necessary per iteration.

The computational complexity of Algorithm 2 is mainly dominated by several invocations of the FFT per iteration. It will run in $O(nx⋅ny⋅nm⋅log⁡(nx⋅ny⋅nm))$ time, which is significantly less than the sparse matrix stacking minimization approach presented before for most circumstances.

In order to test the experimental approach in conjunction with both proposed inversion algorithms as an extension to laser-based photothermal SR reconstruction, a purpose-made sample additively manufactured from 316L stainless steel has been examined. Within the ROI as shown in Fig. 4 lie multiple internal defects in the form of cubical voids with side length 2 mm, filled with unfused metal powder from the manufacturing process. The square OuT features a side length of 58.5 mm and a thickness of $L=4.5mm$⁠, while its fused bulk material comes with a thermal conductivity of $k=15Wm−1K$⁠, a density of $ρ=7950kgm−3$⁠, a specific heat capacity of $cp=502Jkg−1K−1$⁠, and, therefore, an inferred thermal diffusivity of $α=3.76×10−6m2s−1$⁠.^22,23

The ROI is chosen such that it covers several pairs of defects with a spacing of 0.5, 1, 2, and 4 mm between them. This allows not only for testing if all defects are detected by the algorithm but also if all defects can be detected as individual defects resulting in a measure of the resolution capabilities of the method.

As a photothermal heating source, a $Q^max=500W$ fiber-couple diode laser with a wavelength of $λ=940nm$ is used. With this, the OuT front surface can be illuminated via reflecting off a dichroitic mirror (see Fig. 5), which is highly reflective for the laser wavelength but transparent for midwave infrared radiation. The surface temperature is subsequently measured with a midwave infrared camera with a detector size of $npix=1280×1024pixel$ at a spatial resolution of $Δx,Δy=52μm$ per pixel with a frame rate of $fcam=100Hz$⁠. Two linear translational stages are utilized to move the sample with respect to the laser and allow for scanning of the whole surface of the OuT.

For a regular measurement grid of equilateral triangles of side length $rd=0.743mm$ across the ROI with $nm=403$ vertices/independent measurement positions arranged in $7$ rows with $54$ or $53$ positions each and subsequent reconstruction using the sparse matrix stacking approach, the resulting reconstruction of $arec$ is shown in Fig. 7. At each grid point, a single measurement has been recorded after an illumination with a laser spot with spot size $dspot=0.6mm$ at $Q^=15W$⁠. To eliminate the time dimension, an evaluation time of $teval=500ms$ has been taken into account. This corresponds to a thermal diffusion length of $Ldiff≈1.5mm$⁠.

The reconstruction result for applying the sparse matrix stacking algorithm is shown in Fig. 7. For this reconstruction, the ADMM-parameters $λ2,1=1570$⁠, $λ2=100$⁠, $ρ=16$⁠, and $niter=400$ have been used. As can be seen, all defects have been resolved as independent defects and the measurement noise has been successfully suppressed. The reconstruction of the exact defect geometry is still to be improved. For reference, in Fig. 6, the measured data for illuminating the whole OuT front surface homogeneously using the same setup as for the point-wise illumination with $Q^=500W$ and a pulse length of $tpulse=500ms$ are presented.

The reconstruction result for reconstructing $arec$ within the spatial frequency domain is shown in Fig. 8. Here, the ADMM-parameters $λ2,1=27$⁠, $λ2=500$⁠, $ρ=16$⁠, and $niter=400$ have been used to achieve the shown reconstruction results. With reconstruction within the spatial frequency domain, a speedup of approximately 50 times has been found for the reconstruction within the given ROI, while this method is also much more sensitive to perturbation, making it harder to find a suitable set of parameters for the reconstruction.

In contrast to the sparse matrix stacking reconstruction result, this method also reconstructs the edges of the illuminated area as defects as can be seen near the left and right borders of the ROI. This can be explained by the local violation of the necessary condition for super resolution reconstruction given in Eq. (8). The sparse matrix stacking inversion method has been experienced to show much less dependence on the smoothness of the external excitation pattern. Even though both numerical inversion methods mathematically approximate the true solution to the same reconstruction problem, both methods show considerably different results. This can not only be explained due to performing the inversion either in the flattened spatial domain (sparse matrix stacking) or the spatial frequency domain (FFT-based method) but also by the high non-linearity and filtering properties of the used solver in form of the ADMM algorithm and the applied regularizing terms. How to exploit those differences to improve the reconstruction quality is still part of an ongoing research.

In order to better assess the resolution gain of the proposed photothermal SR reconstruction methods compared to well-established conventional methods, a qualitative comparison of different methods can be found in Fig. 9. Here, the reconstruction results of both proposed photothermal SR reconstruction methods is shown next to raw data acquired after homogeneous illumination (cf. Fig. 6) as-is and after different post-processing steps. The homogeneous illumination data are recorded within the same setup as shown in Fig. 5. The additional post-processing steps shown contain a difference thermogram where from the measured raw data $Tdiff$⁠, the temperature evolution of a region of the same OuT without defects is subtracted. As a second post-processing method, the amplitude and phase data for a pulsed phase thermography (PPT) of $Tdiff$ for an evaluation frequency of $fPPT=0.1Hz$ (found to feature the highest contrast for this dataset) are shown.²⁴ 

Since quantitatively comparing the increase in resolution capability comparing the sparse reconstruction results gained from photothermal SR reconstruction and dense measurement data acquired from conventional methods is not trivially possible, we limit ourselves to a qualitative comparison only where our main criterion comprises the separability of closely spaced defects by the method. As can be seen in Fig. 9, in the difference thermogram, as well as a PPT analysis of $Tdiff$⁠, all defects are detected, but the smallest gap between defects (⁠$dgap=0.5mm$⁠) cannot be clearly identified with the available SNR of the measurement. In this direct comparison, the advantage of the sparse nature of the reconstruction of the photothermal SR reconstruction results can be clearly identified. Here, all defects can be clearly separated.

Within this work, an experimental approach for expanding the capabilities of laser-based photothermal SR reconstruction to the detection of two-dimensional defects has been successfully introduced. Two numerical inversion methods for the reconstruction of the internal defect/inhomogeneity structure of an OuT have been proposed and their reconstruction performance has been tested experimentally. Both proposed numerical methods have proven to make it possible to detect even closely positioned defects individually and both methods have provided an approximation of the defect shape. While the proposed sparse matrix stacking inversion method led to qualitatively better shape reconstruction, it took almost 50 times more time to compute on modern computer hardware, making it less suited for the reconstruction of large ROIs. The reconstruction in the frequency domain has proven to be more sensitive to the choice of regularization parameters and to imperfections in the measurements but takes less time to compute and is more memory efficient. Both methods are still equally affected by the lack of a suitable method for finding the best suited set weights for the regularizing terms programmatically in a feasible amount of time. Solving this issue and the improvement of the reconstruction quality of the defect shape is still part of an ongoing research.¹⁶ 

As the presented experimental approach can be recreated using pre-existing standard equipment, a simple adaptation for existing industrial test processes is possible. Its capability of achieving super resolution further allows for the detection of defects beyond the spatial resolution limit of the conventional methods presented. In addition, the sparse nature of the reconstruction results highly contrasts defective regions and allows for easy automation of the test result evaluation process. However, the experimental prerequisite to perform a large amount of sequential independent measurements with single laser-spot excitation to cover larger ROIs is still to be improved for this method. While for testing scenarios demanding high-defect detection resolution at larger depths, the increased testing effort can be justified, e.g., for high-performance parts for aerospace applications or for medical diagnostics, for high-volume/low-cost testing scenarios, the experimental complexity still needs to be improved. Here, the projection of two-dimensional pixel patterns using the latest laser-projector technology shows promising results to significantly reduce the necessary number of measurements.²⁵ 

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.