# Enhanced terahertz nondestructive evaluation and stratigraphic reconstruction with filtered deconvolution

Haolian Shi, Grace Ham, Philip Meilland, D.S. Citrin, Alexandre Locquet

► To cite this version:

Haolian Shi, Grace Ham, Philip Meilland, D.S. Citrin, Alexandre Locquet. Enhanced terahertz non-destructive evaluation and stratigraphic reconstruction with filtered deconvolution. NDT & E International, 2025, 156, pp.103440. 10.1016/j.ndteint.2025.103440. hal-05176204

HAL Id: hal-05176204

https://cnrs.hal.science/hal-05176204v1

Submitted on 22 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Research Paper

# Enhanced terahertz nondestructive evaluation and stratigraphic reconstruction with filtered deconvolution

Haolian Shi<sup>a,b,*</sup>, Grace Ham<sup>c</sup>, Philip Meilland<sup>c</sup>, D.S. Citrin<sup>a,b@</sup>, Alexandre Locquet<sup>a,b</sup>

ARTICLE INFO

Keywords:

THz-TDS

THz time-of-flight tomography

Sparse deconvolution

Orthogonal matching pursuit

Stratigraphic reconstruction

Nondestructive evaluation

ABSTRACT

Terahertz pulsed imaging (TPI) is a valuable tool for nondestructive evaluation (NDE), enabling high-resolution stratigraphic characterization and structural inspection of various materials. However, its accuracy — defined by both the precise localization of peaks and the reduction of spurious peaks — can be limited by challenges such as low signal-to-noise ratio (SNR) and temporally overlapping echoes when layers are thin. While advanced signal-processing methods such as sparse deconvolution (SD) and orthogonal matching pursuit (OMP) have shown promise in addressing these issues, their effectiveness often hinges on careful parameter selection, which typically requires significant expertise. This paper presents a novel approach that enhances the accuracy and simplifies parameter determination in THz TPI by applying SD to a filtered impulse response function in the time domain. This method also aids in OMP parameter determination, further advancing the precision and usability of THz TPI for NDE applications. We demonstrate the effectiveness of the proposed approach by applying it to experimental data obtained from Willow glass on silicon and mill scale on hot-rolled steel strip, showcasing its versatility across different material systems.

1. Introduction

Terahertz (THz) waves, falling within the frequency range 0.1–10 THz, bridge the gap between microwaves and infrared light in the electromagnetic spectrum. Radiation in the THz band can penetrate many electrical insulators, allowing nondestructive inspection of an object’s internal structure. Imaging at THz frequencies can provide superior resolution, both transverse and axial—compared with several well-established NDE techniques such as ultrasonic imaging where high frequencies may be strongly attenuated in some materials [1]. Furthermore, the lack of a couplant simplifies the inspection process. THz imaging is also noncontact. Its capacity to penetrate various electrically insulating materials, as noted above, positions THz technology as a powerful emerging tool for subsurface NDE. This technology has found diverse applications, including measuring coating thickness [2], analyzing polymer–matrix composites [3], inspecting glass-reinforced polymer composites [4–6], residual stress distribution [7] and characterizing cultural-heritage objects [8–11].

Among THz-based techniques, THz pulsed imaging (TPI), also called THz time-of-flight tomography, stands out for its ability to investigate material structures in the axial (depth) direction. By analyzing time delays between reflected pulses, or *echoes*, of an incident THz pulse

in the reflected THz signal, TPI can accurately gauge the thickness of each layer within a sample enabling stratigraphic reconstruction. The deconvolution method proposed by [12] attempts to reconstruct the impulse response function by applying an inverse Fourier transform to the ratio of the reflected and incident THz signals. In ideal cases, the impulse consists of sharp pulses. However, factors such as low signal-to-noise ratio (SNR) and the temporal overlap of echoes, especially when layers are optically thin, pose challenges to the accuracy of TPI. From a hardware perspective, advanced THz sources, especially broadband sources made from various materials, can enhance signal quality and improve performance [13,14]. Alternatively, advanced signal processing methods can provide a more direct and rapid solution for improving TPI's utility in stratigraphic reconstruction [15–17]. Although machine learning and deep learning-based techniques have gained significant attention in recent years [18,19], they have limitations, including reduced interpretability and a reliance on large amounts of high-quality training data. Consequently, classical signal processing techniques remain highly relevant. Frequency-wavelet domain deconvolution (FWDD) [20] can enhance SNR by initially applying a Wiener filter, followed by residual noise reduction using stationary wavelet decomposition. It must be noted that FWDD relies on filters and thus cannot reconstruct information that is missing from the signal.

<footnote><sup>a</sup> Georgia Tech-CNRS IRL 2958, Georgia Tech-Europe, 2 Rue Marconi, 57070 Metz, France</footnote>

<footnote><sup>b</sup> School of Electrical and Computer Engineering, Georgia Institute of Technology, Atlanta, GA, 30332-0250, USA</footnote>

<footnote><sup>c</sup> ArcelorMittal Maizière Research S.A., Voie Romaine, Maizière-lès-Metz, 57280, France</footnote>

<footnote>* Corresponding author at: School of Electrical and Computer Engineering, Georgia Institute of Technology, Atlanta, GA, 30332-0250, USA.
E-mail address: hshi94@gatech.edu (H. Shi).</footnote>
Autoregressive spectral extrapolation (AR) is another signal processing tool that has been employed, particularly for optically thin layers [21,22]. It uses tools similar to those of time-series prediction, but in the frequency domain, and aims at extending the high-SNR band of the transfer function into regions with lower SNR. AR can provide axial super-resolution capability.

Another approach leverages the sparsity of the impulse response function, which in a stratified medium consists of time-delayed sharp peaks corresponding to echoes at material interfaces, whose relative amplitudes are described by Fresnel equations. Sparse deconvolution (SD) has shown great promise in identifying echo positions, sometimes achieving super-resolution in the axial direction. SD assumes that the reflected signal contains only a few echoes, making the impulse response function sparse in the time domain.

Specifically, SD aims to minimize the difference between the reconstructed and measured signals, with regularization methods employed to enforce sparsity. Applying the L₀-norm typically leads to an NP-hard problem (see below). While orthogonal matching pursuit (OMP) [23–25] can sometimes provide an acceptable solution, especially for simpler structures [26], it lacks a guaranteed global optimum. By relaxing the regularization constraint to the L₁-norm, the problem transforms into a convex optimization, which is the standard approach in SD. Numerous successful applications of SD have been reported [27], though its effectiveness often hinges on carefully selecting parameters, a process requiring substantial expertise.

THz techniques are also employed for analyzing complex structures in transmission mode. For example, the normalized time-domain integration (NTDI) method [28] has been used to detect delamination in wood materials by enhancing the signal from multiple reflections within the delaminated region. The effectiveness of NTDI has also been demonstrated in characterizing low-velocity impact damage in composites [28].

The TPI experiment can be viewed as a linear time-invariant (LTI) system. The input is the THz electric-field pulse generated by the apparatus $x(t)$, the output is the reflected THz electric-field signal $y(t)$, and the impulse response $h(t)$ is determined by the stratigraphy of the structure at the location interrogated. Because a structure with discrete layers will lead to well-defined reflections, $h(t)$ is composed of a sequence of narrow peaks, but at most times $t$, $h(t) = 0$, i.e., $h(t)$ is sparse. Identifying the peaks in $h(t)$, including their number and temporal positions is always important. In the context of defect detection, the presence of an anomalous (unexpected) peak tends to indicate a structural defect within the sample. In the context of stratigraphic reconstruction, the delay between successive peaks determines the value of thickness of the layer. Therefore, identifying peaks in $h(t)$ is crucial for accurately determining the sample’s axial structure and detecting potential anomalies. While the LTI system assumption is widely used in THz-TDS, its applicability may be limited in scenarios involving samples with strong dispersion, significant nonlinear attenuation, or high structural complexity. Despite SD’s success in various contexts, hitherto it operates on raw, noisy THz signals without benefiting from the filtering of low-SNR spectral regions, as seen in FWDD. To address this, we propose a new method that combines the strengths of filtering with SD. In this paper, we present, for the first time, an approach that applies SD to a filtered impulse response function rather than directly to the raw time-domain signal. This not only improves accuracy but also simplifies the selection of SD parameters. This enables effective automation of SD analysis for TPI data in stratigraphic reconstruction. Additionally, we find that this method facilitates parameter determination for OMP.

## 2. Methodology

### 2.1. Problem definition

The TPI measurement is typically considered as a convolutional problem,

$$y(t) = h(t) * x(t) + e(t) \qquad (1)$$

where $y(t)$ is the measured (reflected), $x(t)$ the incident or reference signal, $h(t)$ the impulse response, and $e(t)$ additive noise. In practice, measurement of the signal involves sampling, which makes all the processing techniques based on discrete-time signals. With an assumed sampling time $t_0$, the continuous-time signal can be represented as $x[n] = x(n t_0)$. Thus the process above can be rewritten as

$$y[n] = h[n] * x[n] + e[n] = \sum_{k=-\infty}^{\infty} h[k] x[n-k] + e[n]. \qquad (2)$$

Furthermore, Eq. (2) can be written in matrix form as

$$y = Xh + e \qquad (3)$$

where $X$ is an $N \times M$ matrix each column of which represents a delayed version of the reference signal, $y$ denotes the measured signal in discrete time, and $h$ represents the system’s impulse response function. The structural information of the sample is usually encapsulated in this impulse response function.

### 2.2. Conventional approaches in the frequency and time domains

In the following, a circumflex above the impulse response function indicates that this is a reconstruction from the data.

1. Neglecting noise, we can determine an estimate of $h[n]$ as

$$\hat{h}_1[n] = F^{-1}\left(\frac{Y(\omega)}{X(\omega)}\right), \qquad (4)$$

where $Y(\omega)$ and $X(\omega)$ are the Fourier transforms of the reflected and reference signals $y[n]$ and $x[n]$ and $F^{-1}(\cdot)$ is the inverse Fourier transform. This estimated $\hat{h}_1[n]$ obtained from Eq. (4) is often poor due to the presence of noise.

2. In the time domain, reconstruction is carried out by solving

$$y[n] = \hat{h}[n] * x[n] \qquad (5)$$

and introducing a regularization term to favor sparsity, i.e., finding $\hat{h}_{1,s}[n]$ that solves

$$\min_{\hat{h}_{1,s}} \|X\hat{h}_{1,s} - y\|_2^2 + \lambda \|\hat{h}_{1,s}\|_1 \qquad (6)$$

where the first term in the minimization represents accuracy, while the second quantifies sparsity. The subscripts 2 and 1 refer to the $L_2$ and $L_1$ norm, respectively. This $\hat{h}_{1,s}[n]$ can be regarded as an estimation of the impulse response function $h[n]$ by SD. There are other options to regulate the sparsity such as the $L_0$ norm; however, in this case, the problem becomes NP-hard; a suboptimal solution can be found with OMP. In both cases, noise is neglected in the reconstruction, though it is present in the measured signals $x[n]$ and $y[n]$.

### 2.3. Filtering the transfer function

Better estimates could be obtained by filtering out spectral regions where the SNR is low.

1. In the frequency domain, filtering is applied before taking the inverse Fourier transform, viz.

$$\hat{h}_2[n] = F^{-1}\left[\frac{Y(\omega)}{X(\omega)}Fil(\omega)\right]. \qquad (7)$$

This estimate is what is obtained in the FWDD technique, before wavelet denoising is applied. $Fil(\omega)$ is the transfer function of the filter applied.

2. In the time domain by solving the convolutional equivalent of Eq. (7),

$$\hat{h}_2[n] = h_1[n] * fil[n] \qquad (8)$$

or, in matrix form where $fil[n]$ is the filter impulse response

$$\hat{h}_2 = F\hat{h}_1 \qquad (9)$$
with F a matrix where each column is a delayed version of the impulse response of the filter fil[n]. At this point, the problem exhibits the same structure as Eq. (5), making it amenable to resolution using SD. As is typical in SD, we introduce a regularization term, based on the L₁ norm, to favor sparsity, i.e., we determine $\hat{h}_{2,s}$ that solves

$$ \min_{\hat{h}_{2,s}} \left\| F\hat{h}_{2,s} - \hat{h}_{2} \right\|_{2}^{2} + \lambda \left\| \hat{h}_{2,s} \right\|_{1}. \quad (10) $$

The objective of deconvolution is to retrieve the function $\hat{h}_{2,s}$ with sparse regularization. $\hat{h}_{2,s}$ here can be considered as an estimation of $\hat{h}_1$. This filtering approach offers a new interpretation of the entire process and introduces benefits in parameter determination, which will be discussed in the following sections. The choice of the filter, including its type and cutoff frequency, can vary depending on the signal characteristics and specific applications.

## 2.4. Techniques for solving the optimization problems

The optimization problems of Eqs. (6) and (10), also known as the LASSO problem, can be solved using an iterative shrinkage algorithm [29,30]. The iterative process can be described as

$$ h^{k+1} = \text{Shrinkage} [h^k - \lambda \nabla f(h^k), \lambda \tau] \quad (11) $$

where k denotes the iteration index and τ is the step size. The gradient of the reconstruction error term is given by $\nabla f(h^k) = X^T(Xh^k - y)$. The shrinkage operator is defined element-wise as

$$ \text{Shrinkage}(f[n]) = \begin{cases} f[n] + \tau\lambda, & f[n] \le -\tau\lambda, \\ 0, & |f[n]| < \tau\lambda, \\ f[n] - \tau\lambda, & f[n] \ge \tau\lambda. \end{cases} \quad (12) $$

In practical applications, the presence of noise and dispersion can complicate the reconstruction process. For SD, iterative methods may cause peak splitting (the appearance of more than one peak in the reconstructed impulse response function for a given single real peak) due to peak broadening, especially in the presence of strong dispersion. Dispersion is typically significant when layers are thick. For OMP, determining the number of nonzero elements in the impulse response function a priori is challenging. Therefore, selecting an appropriate sparsity coefficient $\lambda$ in Eq. (6) is crucial.

Of note, the sparsity regularization term can also be expressed with the $L_0$ norm. As mentioned before, this problem is NP-hard; however; it can be partially addressed using OMP. In OMP, the algorithm calculates the correlation between the measured signal and each column of the matrix X. The column with the highest correlation is selected and the corresponding component is removed from the matrix. This process will be repeated iteratively until the residual approaches zero.

In work to date, the value of $\lambda$ in both SD and OMP is typically determined using the L-curve method [31], which involves plotting the trade-off between sparsity and reconstruction error and identifying the point where both metrics are minimized. Specifically, each $\lambda$ corresponds to a unique optimized impulse response function and thus to a distinct reconstruction error and sparsity level representing a single point on the L curve. To construct the L curve, a series of $\lambda$ values are tested with each point on the curve defined by the associated reconstruction error (horizontal axis) and sparsity (vertical axis). The optimal $\lambda$ is typically located at the corner of the L-curve, where its second derivative reaches a maximum. This point represents the inflection where the balance between sparsity and reconstruction error is optimized, minimizing both. However, this approach can be problematic due to the low SNR in real measurements, which may prevent the formation of a clear L-shaped curve. The model vacillates between reducing reconstruction error and enforcing sparsity. Consequently, finding the optimal number of iterations remains difficult.

<figure><img src="image_2.png" /><figcaption>Fig. 1. THz spectrum of the impulse response function with and without filtering. The black line represents the raw spectrum $F[\hat{h}_1]$ without filtering, while the red line corresponds to the spectrum $F[\hat{h}_2]$ after applying the filter described in Eq. (13). The vertical axis is displayed on a logarithmic scale. An abrupt transition is clearly observed. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

Enhancing the SNR would significantly help address the challenge of determining the optimal parameter mentioned above. This corresponds to the approach leading to the optimization problem expressed as Eq. (10). For the filter, we choose a Hanning window to suppress noise outside the band in which we have meaningful information,

$$ \text{Fil}(\omega) = \begin{cases} e^{i\omega t_0} \cos^2\left(\frac{\omega}{4f_c}\right) & |\omega| \le 2\pi f_c, \\ 0 & |\omega| > 2\pi f_c, \end{cases} \quad (13) $$

where $t_0$ corresponds to the optical delay for the peak in the reference signal and $f_c$ is the cutoff frequency that defines the band containing the important information. The typical value for $f_c$ for our experiments is 3 THz, and is set by a simple examination of the power spectrum of the signals where $f_c$ is the frequency beyond which the spectrum falls below the noise floor. For other systems and signals, the cutoff frequency and filter type can be adjusted accordingly. The effect of adding the filtering to the impulse response function can be observed in Fig. 1. The red line represents the spectrum of the impulse response function after filtering, while the black line shows the unfiltered spectrum. Filtering has significantly reduced the noise floor.

## 3. Experiment

To demonstrate the efficacy of the proposed method, several samples with various numbers of layers and layer thicknesses were tested. Additionally, the samples featured a diverse range of material types. All signals were measured using a TeraPulse Lx system (TeraView Ltd., Cambridge, UK) in reflection mode. The first sample is a two-layer structure with Willow alkaline earth boro-aluminosilicate glass (Corning) as the first layer and silicon as the substrate. This sample exhibits both positive and negative reflections, corresponding to interfaces between materials with relatively low-to-high and high-to-low refractive indices. The thickness of the glass layer is 125 µm and the silicon substrate has a thickness of 675 µm. The second group of samples consists of hot-rolled steel strip coupons covered by mill-scale films of thicknesses 5.1, 7.9, and 14 µm. The 5.1 µm sample was also studied in [2].

The thickness $d$ of a layer is obtained from the time delay $\Delta t$ between the two associated peaks in the reconstructed impulse response
<figure><img src="image_3.png" /><figcaption>Fig. 2. (a) A-scan for the Willow glass on silicon sample. The orange curve is the measured reference signal and the blue curve is the reflected signal. The inset schematically represents the structure of the sample under measurement; multiple reflections are weak and not considered. The thickness of glass is 125 µm and the silicon substrate is 675 µm. (b) The photograph for the Willow glass on silicon sample. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. The L-curve for Willow glass on silicon employing the unfiltered SD method. This curve illustrates the trade-off between reconstruction error (horizontal axis) and sparsity (vertical axis), calculated according to Eq. (6) (the first term for error and the second for sparsity).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. The impulse response function $\hat{h}_{1,1}$ of the Willow glass on silicon sample for the conventional SD method. Note the erroneous splitting of some of the peaks.</figcaption></figure>

function by

$$d = \frac{c \Delta t}{2n} \qquad (14)$$

where c is the speed of light in vacuum and n is the refractive index of the layer. The factor of 2 in the denominator is due to the fact that the echo makes a round trip through the layer.

<figure><img src="image_6.png" /><figcaption>Fig. 5. L-curve of the Willow glass on silicon sample with filtering. This curve illustrates the trade-off between reconstruction error (horizontal axis) and sparsity (vertical axis) in the solution, calculated according to Eq. (10) (the first term for error and the second for sparsity). The curve appears smoother compared to the conventional SD method shown in Fig. 3.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. The impulse response function $\hat{h}_{2,1}$ of the Willow glass on silicon sample using the filtering method. Peak splitting is suppressed compared to the unfiltered method shown in Fig. 4.</figcaption></figure>

### 3.1. Willow glass

Fig. 2 shows the measured signal reflected from the sample. The orange curve is the reference signal while the blue curve is the measured reflected signal. The reference signal can be seen to be a roughly single-cycle pulse. It is followed by weak ringing caused by resonances
<figure><img src="image_8.png" /><figcaption>Fig. 7. Reconstructed impulse responses (SD) along a cross-section: unfiltered (left, $\hat{h}_{1,s}$) vs. filtered (right, $\hat{h}_{2,s}$). Filtering enhances interface delineation and contrast.</figcaption></figure>

in atmospheric water vapor and background noise. The structure of the sample is shown schematically in the inset. There are three main reflections as indicated in the inset. Each reflection resembles a time-delayed version of $x[n]$ or $-x[n]$. The first echo is from the top air/glass interface, the second from glass/substrate interface, and the third from the bottom substrate/air interface as can be seen in the blue curve.

Each point in the curve is related to a $\lambda$. For each $\lambda$, we run the iteration and obtain an optimized impulse response function. This impulse response function corresponds to a reconstruction error and sparsity. First, we test SD with the raw signal; this is the standard approach hitherto employed. Fig. 3 displays the L-curve when filtering is not used. Each point on the L-curve is related to one specific $\lambda$ which balances the reconstruction error and sparsity of the reconstruction in Eq. (6). Ideally the optimal parameter $\lambda$ should be chosen when both the sparsity and reconstruction errors are low, corresponding to the elbow point of the L-curve. However, in the conventional method shown in Fig. 3, identifying such a point is challenging due to the complex, spiky shape of the L-curve. The curve significantly deviates from an ideal L-curve, complicating the task of determining the optimal $\lambda$. We decide on a value of $\lambda = 0.02$ (corresponding to the location of the abrupt change in Fig. 3) and plot the reconstructed impulse response function by SD in Fig. 4. Although we can identify the three main peaks from the reconstructed impulse response function, there are several small peaks that can be misleading in the identification of real reflections. It is clearly desired to suppress such spurious peaks.

By contrast, Fig. 5 is the L-curve for the filtered method introduced here. The L-curve is far smoother than in the unfiltered case of Fig. 3. The elbow (maximum of the second derivative) of the L-curve is also easily found. The reconstruction error for the optimal $\lambda$ ranges from 0.03 to 0.04, while the model’s sparsity lies between 0.4 and 0.6.

Fig. 6 shows the deconvolved impulse response function $\hat{h}_{2,s}$ using the filtering method. Compared to Fig. 4, the three genuine peaks of the impulse response function are clearly identified, while spurious peaks are notably absent. Fig. 7 presents a comparative visualization of the reconstructed impulse response along a cross-section, contrasting the unfiltered (left) and filtered (right) SD results. The filtered approach demonstrates improved clarity, with sharper definition of the interface responses and enhanced contrast throughout the cross-sectional image.

The filtering approach also improves parameter determination for OMP. Fig. 8 shows the relationship between the reconstruction error and number of iterations used in OMP on the same Willow glass-silicon data. The black curve represents the direct deconvolution of

the time-domain signal, while the red curve corresponds to the new filtering method. Based on the sample’s structure, the true number of peaks is three. This can be observed in both curves; before three iterations, the error drops quickly and after three iterations, the error stabilizes. However, in the case of filtering, the rate of change becomes more distinct compared to conventional SD; beyond the optimal point, the curve is straighter. This characteristic can assist in identifying the optimal number of iterations when applying the OMP method to complex samples.

## 3.2. Mill scale

Mill scale is a thin layer of oxide on the steel surface as hot-rolled steel strip cools. The three main iron-oxide phases present are wüstite, magnetite, and hematite. The majority phase wüstite is electrically insulating. The mill scale is subsequently removed on the pickling line before further processing, as its presence can hinder operations such as coating deposition. Pickling produces considerable toxic waste, and thus information on the mill-scale thickness is of great importance to optimize the mill-scale removal and to reduce waste while preventing over pickling, which might damage the steel surface.

The raw reflected signal is shown in Fig. 9. The orange curve is the reference signal and the blue curve is the reflected signal. The inset shows a schematic diagram of the sample and the measurement. Unlike the previous glass-silicon sample, the mill-scale/steel interface exhibits total reflection, resulting in the absence of a final sample/air interface reflection. Three different mill-scale thicknesses (5.1, 7.9, and 14.8 µm) were tested, with these values determined using the dissolution method.

The L-curve for conventional SD is shown in the first row of Fig. 10. Although the shape of the L-curve is preserved, an abrupt transition remains clearly visible for the thicker cases (7.9 and 14.8 µm). In contrast, with the proposed filtered SD method, the L-curve becomes significantly clearer, making it relatively easier to identify the optimal sparsity regularization parameter. More importantly, the reconstructed impulse response function $\hat{h}_{2,s}$ demonstrates a higher SNR compared to $\hat{h}_{1,s}$ obtained through conventional SD with optimal parameter choice, as shown in the last row of Fig. 10. For the 5.1 µm case, two peaks can be distinguished clearly for the filtered SD. Additionally, the reconstructed impulse response function exhibits fewer small spurious peaks, enhancing the reliability of the deconvolution process. Table 1
<figure><img src="image_9.png" /><figcaption>Fig. 8. The relationship between the reconstruction error and number of iterations for OMP. The red curve represents the filtered method, while the black curve corresponds to the unfiltered method. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9. Typical time domain signals for the mill-scale samples. The inset shows a schematic representation of the setup. The orange curve represents the measured reference signal and the blue curve corresponds to the reflected signal. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

compares the thickness estimated by the various deconvolution methods. The reference thickness is obtained by the dissolution method, as mentioned above. In all cases, the thickness reconstructed using the new method is closer to the reference thickness. The assumed refractive index is 4.7, as in the previous work [32]. It is important to note that the filtered SD approach not only excels in estimating the thickness of thin mill-scale layers but also eliminates false peaks in the impulse response function, resulting in a clearer and more accurate impulse response function compared to unfiltered SD, as shown in the last row of Fig. 10.

## 4. Conclusion

In this work, we propose a new method to enhance the accuracy of sparse deconvolution (SD) and simplify the determination of optimal parameters for both SD and orthogonal matching pursuit (OMP). Rather than directly deconvolving the reference and measured reflected signals in the time domain, the proposed approach achieves higher accuracy and easier parameter determination by deconvolving the filtered impulse response function with the filter itself in the time domain. The

<table><caption>**Table 1**
Comparison of mill-scale thickness estimation *d* by SD with and without filtering compared with the reference thickness measured by the dissolution method.</caption><thead><tr><td>Reference thickness (μm)</td><td>Conventional SD (μm)</td><td>Filtered SD (μm)</td></tr></thead><tbody><tr><td>5.1</td><td>9.9</td><td>6.6</td></tr><tr><td>7.9</td><td>10.8</td><td>7.7</td></tr><tr><td>14.8</td><td>17.8</td><td>14.0</td></tr></tbody></table>

effectiveness of this method has been validated through real measurements involving Willow glass on silicon and mill-scale layers of varying thicknesses on steel. We present representative examples demonstrating this improvement. This approach offers a practical and efficient solution for enhancing signal reconstruction quality, particularly in scenarios with high noise levels or closely spaced peaks (optically thin layers) in the impulse response function. It should be emphasized that, until now, considerable expertise has been required to manually
<figure><img src="image_11.png" /><figcaption>Fig. 10. L-curves and SD deconvolution results for mill scale of various thicknesses (L-curves are plotted with the same convention as Figs. 3 and 5). The three columns correspond to mill-scale thickness of 5.1, 7.9, and 14.8 µm. (a)–(c) The L-curve with SD of raw signal. (d)–(f) The L-curve with filtered SD. (g)–(i) The impulse response function reconstructed by the deconvolution method. The black curve is with conventional SD and the red curve is with the filtered SD. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

optimize the parameter when implementing SD. Consequently, the proposed approach is well-suited for NDE technicians, reducing reliance on highly specialized academic training. Furthermore, this method shows potential for future applications in more complex signal environments and integration with other advanced deconvolution techniques.

### CRediT authorship contribution statement

**Haolian Shi:** Writing – review & editing, Writing – original draft, Conceptualization.
**Grace Ham:** Writing – review & editing, Supervision.
**Philip Meilland:** Writing – review & editing, Supervision.
**D.S. Citrin:** Writing – review & editing, Supervision, Conceptualization.
**Alexandre Locquet:** Writing – review & editing, Supervision, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgments

We thank Shelby Nelson of Mosaic Microsystems for the Willow glass on silicon. We gratefully acknowledge funding from ArcelorMittal Mazière Research, France. This work was supported by the French PIA project "Lorraine Université d'Excellence", reference ANR-15-IDEX-04-LUE.

### Data availability

Data will be made available on request.

### References

[1] Dong J, Kim B, Locquet A, McKeon P, Declercq N, Citrin DS. Nondestructive evaluation of forced delamination in glass fiber-reinforced composites by terahertz and ultrasonic waves. Compos Part B: Eng 2015;79:667-75.

[2] Zhai M, Locquet A, Roquelet C, Alexandre P, Daheron L, Citrin DS. Nondestructive measurement of mill-scale thickness on steel by terahertz time-of-flight tomography. Surf Coat Technol 2020;393:125765.

[3] Zhu P, Zhang H, Robitaille F, Maldague X. Terahertz time-domain spectroscopy for the inspection of dry fibre preforms. NDT E Int 2024;145:103133.

[4] Kim H-S, Park D-W, Kim S-I, Oh G-H, Kim H-S. Non-destructive detection of thin micro-defects in glass reinforced polymer composites using a terahertz electromagnetic wave based on a convolution neural network. Compos Part B: Eng 2023;257:110694.

[5] Zhai M, Locquet A, Citrin DS. Terahertz nondestructive layer thickness measurement and delamination characterization of GFRP laminates. NDT E Int 2024.

[6] Pomarède P, de la Rosa JC, Antonik P, Meraghni F, Citrin DS, Rontani D, et al. THz-based monitoring of the deformation of the inner woven structure of glass-fiber reinforced polymers. IEEE Trans Terahertz Sci Technol 2023;13:396-404.

[7] Zhu P, Zhang H, Santulli C, Sfarra S, Usamentiaga R, Vavilov VP, et al. Contactless and nondestructive evaluation of residual stress distribution in modified and pure HDPE materials using a novel terahertz method and line-scan thermographic technique. Compos Part A: Appl Sci Manuf 2024.

[8] Dong J, Locquet A, Melis M, Citrin D. Global mapping of stratigraphy of an old-master painting using sparsity-based terahertz reflectometry. Sci Rep 2017;7(1).

[9] Fukunaga K. Nondestructive evaluation of lined paintings by THz pulsed time-domain imaging. Heritage 2023.
[10] Cui Y, Xu Y, Han D, Wang X, Shen Z, Hou Y, et al. Hidden-information extraction from layered structures through terahertz imaging down to ultralow SNR. Sci Adv 2023;9.

[11] Citrin DS, Locquet A, Dong J. Unearthing a prayer for the dead: Terahertz waves uncover the hidden inscription on a 16th-century lead cross. Am Sci 2023;111(5):280-6.

[12] Walker GC, Bowen JW, Labaune J, Jackson JB, Hadjiloucas S, Roberts J, et al. Terahertz deconvolution. Opt Express 2012;20 25:27230-41.

[13] E Y, Zhang L, Tcypkin A, Kozlov S, Zhang C, Zhang X-C. Broadband THz sources from gases to liquids. Ultrafast Sci 2021.

[14] Mansourzadeh S, Nielson MF, Omar A, Vogel T, Michaelis DJ, Johnson JA, et al. High-power, ultra-broadband THz generation in organic crystal MNA. In: 2023 48th international conference on infrared, millimeter, and terahertz waves (IRMMW-THz). 2023, p. 1-2.

[15] Sun F, Fan M, Cao B, Liu L. THzResNet: A physics-inspired two-stream residual network for thermal barrier coating thickness measurement. IEEE Trans Ind Inform. 2023;19:8328-39.

[16] Zhang Z, Huang Y, Zhong S, Lin T, Zhong Y, Zeng Q, et al. Time of flight improved thermally grown oxide thickness measurement with terahertz spectroscopy. Front Mech Eng 2022;17:1-11.

[17] Yasui T, Yasuda T, i. Sawanaka K, Araki T. Terahertz paintmeter for noncontact monitoring of thickness and drying progress in paint film. Appl Opt 2005;44 32:6849-56.

[18] Wang B, Qin X, Li W, ren Li Z, Zhu L-G. Learned iterative shrinkage and thresholding algorithm for terahertz sparse deconvolution. Opt Express 2022;30 11:18238-49.

[19] Xu Z, Ye D, Chen J, Zhou H. Novel terahertz nondestructive method for measuring the thickness of thin oxide scale using different hybrid machine learning models. Coatings 2020.

[20] Chen Y, Huang S, Pickwell-MacPherson E. Frequency-wavelet domain deconvolution for terahertz reflection imaging and spectroscopy. Opt Express 2010;18(2):1177-90.

[21] Dong J, Locquet A, Citrin DS. Depth resolution enhancement of terahertz deconvolution by autoregressive spectral extrapolation. Opt Lett 2017;42(9):1828-31.

[22] Zhai M, Locquet A, Roquelet C, Citrin D. Terahertz time-of-flight tomography beyond the axial resolution limit: autoregressive spectral estimation based on the modified covariance method. J Infrared Millim Terahertz Waves 2020;41:926-39.

[23] Mallat S, Zhang Z. Matching pursuits with time-frequency dictionaries. IEEE Trans Signal Process 1993;41(12):3397-415.

[24] Tropp JA. Greed is good: algorithmic results for sparse approximation. IEEE Trans Inform Theory 2004;50:2231-42.

[25] Mor E, Azoulay A, Aladjem ME. A matching pursuit method for approximating overlapping ultrasonic echoes. IEEE Trans Ultrason Ferroelectr Freq Control 2010;57.

[26] Shi H, Calvelli S, Zhai M, Ricci M, Laureti S, Singh P, et al. Terahertz nondestructive characterization of conformal coatings for microelectronics packaging. IEEE Trans Components Packag Manuf Technol 2024;14:3-9.

[27] Reyes-Reyes ES, Carriles-Jaimes R, D'Angelo E, Nazir S, Koch-Dandolo CL, Kuester F, et al. Terahertz time-domain imaging for the examination of gilded wooden artifacts. Sci Rep 2024;14.

[28] Zhu P, Zhang H, Pivarčiová E, Sfarra S, Maldague XPV. Characterization of water content and inspection of delamination in spruce, oak and meranti woods after pyrolysis processing using a new terahertz time-domain spectroscopy method. NDT E Int 2023;139:102938.

[29] Elad M, Matalon B, Shtok J, Zibulevsky M. A wide-angle view at iterated shrinkage algorithms. In: Wavelets XII, vol. 6701, SPIE; 2007, p. 15-33.

[30] Wagner CW, Semper S, Kirchhof J. Fastmat: Efficient linear transforms in python. SoftwareX 2022;18:101013.

[31] Yafei X, Wang X, Zhang L, Yan R, Yan R, Yan R, et al. Terahertz nondestructive quantitative characterization for layer thickness based on sparse representation method. NDT E Int 2021.

[32] Henning T, Begemann B, Mutschke H, Dorschner J. Optical properties of oxide dust grains. Astron Astrophys Suppl 1995;112:143.