# Terahertz nondestructive stratigraphic reconstruction of paper stacks based on adaptive sparse deconvolution

Min Zhai, Alexandre Locquet, D.S. Citrin

► To cite this version:

Min Zhai, Alexandre Locquet, D.S. Citrin. Terahertz nondestructive stratigraphic reconstruction of paper stacks based on adaptive sparse deconvolution. Journal of the European Optical Society : Rapid publications, 2024, 20 (1), pp.4. 10.1051/jeos/2024001. hal-04651528

HAL Id: hal-04651528

https://hal.science/hal-04651528v1

Submitted on 17 Jul 2024

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 ARTICLE

OPEN ACCESS

# Terahertz nondestructive stratigraphic reconstruction of paper stacks based on adaptive sparse deconvolution

M. Zhai<sup>1,2</sup>, A. Locquet<sup>1,2</sup>, and D.S. Citrin<sup>1,2,*</sup>

Received 14 September 2023 / Accepted 18 January 2024

**Abstract.** Characterizing the number of sheets in a stack of paper typically involves mechanical separation of the individual sheets. Here, we explore a nondestructive method that can be applied to the intact paper stack. Namely, terahertz time-of-flight tomography, together with post signal-processing technique sparse deconvolution based on a two-step iterative shrinkage-thresholding algorithm (SD/TWIST), is employed to reconstruct the stratigraphy of stacks of sheets of paper with multilayered structure in a nondestructive and noncontact manner. The double-Gaussian mixture model (DGMM) is also incorporated to suppress dispersion in the reflected THz echoes. The effectiveness and accuracy of the proposed adaptive sparse-deconvolution method are verified experimentally and numerically. Compared with the commonly used frequency wavelet-domain deconvolution (FWDD) method and previous implementations of sparse deconvolution based on an iterative-shrinkage and thresholding algorithm (SD/IST), the proposed sparse-deconvolution approach can provide a clearer and rapid stratigraphic reconstruction of the paper stacks studied, while ensuring accurate thickness information for each paper sheet in the presence of noise, revealing the potential usage of real-time THz tomographic-image processing.

**Keywords:** Terahertz time-domain spectroscopy, Paper handling, Super-resolution characterization, Dispersion, Sparse deconvolution.

## 1 Introduction

Nondestructive evaluation (NDE) approaches are needed to avoid wastage, for quality control, and to verify various production processes prior to the next step. An area that could benefit from one-step NDE approaches is paper handling. For example, to count the number of sheets in a paper stack, particularly when more than one paper type is present, the individual sheets must in some fashion be mechanically separated for counting. An approach that avoids such laborious mechanical separation would be welcome. Our work presented here applies terahertz (THz) NDE to achieve this aim.

The THz band of the electromagnetic spectrum, lying at the boundary between photonics and electronics, is of interest in part because THz waves can penetrate numerous electrically insulating materials that may be opaque in the visible and infrared bands or provide little contrast in the X-rays enabling the characterization of structure interior to various objects. Moreover, in contrast to X-ray imaging, THz waves present no known health risks to human tissue

and in contrast to ultrasonic inspection, do not require a coupling medium (e.g., water) [1]. Therefore, THz waves have been considered as a candidate NDE technique for a range of commercial applications, such as steel production [2–4], cultural heritage conservation [5–7], plastics [8–11], medical diagnosis [12–14], and security [15–17].

A method explored in numerous contexts to obtain structural information as a function of depth is THz time-of-flight tomography (TOFT) [18]. In a stratified dielectric medium, the physical thickness of a layer is determined by the optical delay between two echoes reflected from successive interfaces multiplied by the speed of light in that layer, i.e., the in-vacuum speed of light *c* divided by the refractive index *n* in question. However, owing to the temporal overlap of THz echoes reflected from optically thin layers, pulse-spreading associated with frequency-dependent attenuation and dispersion during the propagation of THz waves in materials, as well as low signal-to-noise ratio (SNR) at region of interests because of the low Fresnel coefficients at those interfaces, it may be difficult to distinguish the reflected THz echoes. Our aim in this paper is to demonstrate an efficient and accurate signal-processing approach to nondestructively characterize layer thicknesses in

<footnote n="1"><sup>1</sup>Georgia Tech-CNRS IRL2958, Georgia Tech-Europe, 2 Rue Marconi, 57070 Metz, France</footnote>

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

<footnote>* Corresponding author: david.citrin@ece.gatech.edu</footnote>
complex multilayer structures in the presence of material dispersion. The test case on which we focus are stacks of paper.

To date, several advanced signal-processing methods have been proposed and employed to analyze the complex reflected THz echoes, such as frequency wavelet-domain deconvolution (FWDD) [19] and autoregressive (AR) spectral extrapolation [20, 21]. Based on the assumption that the reflected THz echoes are the time-shifted amplitude-scaled replicas of the THz incident signals, both FWDD and AR methods are mainly limited to deconvolve reflected signals with time-invariant pulses. In practice, however, pulse spreading associated with dispersion (frequency dependence of the refractive index) occurs, and will lower the performance of these deconvolution methods dramatically.

Sparse deconvolution (SD) has attracted attention to retrieve the impulse response function $h(t)$ owing to its robustness, relative insensitivity to noise, and resolvability for overlapped echoes. As a pure time-domain method, SD does not introduce erroneous features into the low- and high- frequency region of the transfer function thus in effect further reducing the axial resolution. The conventional SD approach employs iterative shrinkage and a thresholding algorithm (IST) [22], where the $l_1$-norm regularization problems were solved by adjusting the components of the sparse vector with a fixed step size. Despite the excellent denoising performance of SD/IST, in excess of several hundred iterations is typically required to reach the stopping criterion, which is extremely time-consuming rendering the approach impractical, especially when results are needed in real time [23, 24]. In addition, the waveform changes caused by material dispersion are not dealt with in a straightforward manner by the algorithm. Thus, an approach that both accounts for dispersion in thick samples as well as can identify thin layers is needed.

In this paper, leveraging off the simplicity and the capabilities of SD/IST, an adaptive SD method based on a two-step iterative shrinkage-thresholding algorithm (SD/TWIST) is proposed and employed to quantitatively reconstruct the stratigraphy of a multi-layer sample. As a test case, we focus on a stack of ordinary copy paper. An appropriate dispersion strategy is incorporated into the SD/TWIST to compensate for the dispersion. The criterion of selecting the optimal regularization parameters to analyze time-varying THz signals is also discussed. A series of numerical simulations and experiments are performed to validate the robustness, reliability, and efficiency of the proposed SD/TWIST approach on quantitative thickness measurements compared to conventional THz signal-processing methods, such as FWDD. A clearer stratigraphic analysis with limited erroneous spike-like features in $h(t)$ as well as lower computational cost is demonstrated.

The sequel of this paper is organized as follows. Section 2 contains the description of the experimental setup for THz reflection measurements as well as the theory of the adaptive SD method. Section 3 presents the application of the proposed SD method to address the simulated and experimental dispersive reflected THz signals, respectively. The conclusions are summarized in Section 4.

# 2 Theory and equipment

## 2.1 Experimental setup

A commercial THz system (TPS Spectra 3000, TeraView Ltd., Cambridge, UK) is employed for the THz TOFT measurements. Figure 1 shows the schematic diagram of the THz system in reflection mode. The generation of pulsed THz radiation is based on a photoconductive switch. Quasi-single-cycle THz pulses with a bandwidth ranging from 60 GHz to 3 THz are generated in a biased gallium arsenide (GaAs) antenna after excitations by an Er-doped fiber which emits 780 nm pulses with sub-100 fs pulse duration at a repetition rate of 100 MHz and has an averaged output power >65 mW. By gating the photoconductive gap with an ultrafast optical pulse synchronized to the THz emission, a current proportional to the THz electric field is measured. A delay line is also incorporated into the probe beam to change the difference in the optical delay between the incoming THz pulse and the probe laser pulse at the receiver. A bias is also applied across the emitter to generate a time-gated output signal. The corresponding spectrum is obtained after Fourier transform.

In this study, THz measurements were carried out in an air-conditioned laboratory at 22 °C with humidity <48% at near normal incidence (~3°). The reference THz pulse produced by the apparatus, was recorded by setting a bare metal plate (i.e., an excellent THz reflector) at the sample position. According to Fresnel's law, the reflections happened at the conductive metal plate give phase shifts of $\pi$ respectively. Therefore, the incident THz pulse has the same phase shifts as the reference pulse. The ringing after the main THz pulses and the sharp absorption lines in the measured spectra correspond to the excitation of ro-vibrational modes of atmospheric H$_2$O molecules during the propagation of THz pulses [25]. Measurements can also be made in a dry atmosphere of N$_2$; however, by collecting data under ambient conditions, our work is closer to what might be expected in a production environment. The recorded reflected temporal signal contains 4096 data points (for a sampling period $T_s = 0.0116$ ps) and is averaged over 10 shots per pixel to reduce the abnormal fluctuations associated with system noises.

## 2.2 Model estimation strategy

In the time domain, the measured reflected THz signal $r(t)$ is assumed to be the convolution of the incident signal $i(t)$ with the impulse response function $h(t)$, which represents the structure and the dielectric properties at a given transverse position,

$$r(t) = i(t) \otimes h(t) = \int_{-\infty}^{\infty} i(t - \tau)h(\tau)d\tau, \quad (1)$$

where the symbol $\otimes$ denotes the convolution operation and $\tau$ is the time delayed operation satisfying $t \ge \tau$. Accounting for the noise $e(t)$ originating from the measurement [26], the discrete format of equation (1) can be written as a matrix multiplication,
<figure><img src="image_2.png" /><figcaption>Fig. 1. Schematic diagram of the TPS Spectra 3000 apparatus in reflection mode from TeraView Ltd., Cambridge, UK.</figcaption></figure>

$$r = Ih + e \quad (2)$$

where $I$ is the convolution dictionary matrix whose rows are delayed versions of the discrete form of $i(t)$. $r$ and $e$ are the discretized matrix form of THz reflected signal $r(t)$ and noise $e(t)$. $h$ accounts for the matrix form of the impulse response function $h(t)$ that contain important structural information of sample studied.

Due to the aforementioned sparsity of $h$, the sparse deconvolution aims at approximating the received THz signal $r$ with $Ih$, equation (2) can be solved by $l_1$-norm regularization optimization, which is defined as

$$\min_{h} \frac{1}{2} \|Ih - r\|_{2}^{2} + \lambda \|h\|_{1}, \quad (3)$$

where $\|I'h - r\|_{2}^{2}$ denotes $\sum_{i}^{N}(Ih_{i} - r_{i})^{2}$, $\|h\|_{1}$ denotes $\sum_{i}|h_{i}|$, $N$ is the length of the vector, and $\lambda$ is the non-negative regularization parameter that controls the trade-off between the sparsity of $h$ and the residue norm.

The iterative-shrinkage thresholding (IST) algorithm is an effective method to solve the $l_1$-norm regularization problems. Despite the simplicity and excellent deconvolution performance, the extremely slow convergence rate resulted from the matrix multiplication regarding $I$ and $I^T$ during each iteration, hinders the practical application of the SD/IST method seriously. In this study, a two-step iterative shrinkage-thresholding algorithm (TWIST) is utilized to retrieve $h$, which can be expressed as

$$h_1 = S_T(h_0 + I^T(r - I h_0)), \quad (4)$$

$$h_{k+1} = (1 - \alpha)h_{k-1} + (\alpha - \beta)h_k + \beta S_T(h_k + I^T(r - I h_k)) \quad (5)$$

where $h_0$ is the initial vector, $\alpha$ and $\beta$ are two convergence parameters of TWIST, which are associated with the eigenvalues of matrix multiplication $I^T I$. In addition,

$$S_T(h) = \begin{cases} \operatorname{sgn}(h_i)(|h_i| - \lambda), & |h_i| > \lambda \\ 0, & \text{otherwise}, \end{cases} \quad (6)$$

where $S_T(h)$ is a soft-thresholding or shrinkage operator with respect to each component of the current vector $h_i$. Superior to SD/IST mentioned above, because each iteration of SD/TWIST is nonlinear, it can produce a more robust solution. More details on SD/TWIST method can be found in Refs. [27, 28].

# 3 Results and discussion

## 3.1 Numerical simulation

In this section, before considering experimental data, numerical simulations based on synthetic data are performed to validate the effectiveness and accuracy of the proposed SD/TWIST for the analysis of the reflected THz signal $r(t)$ from a complex layered structure. Theoretically, the reflected signal $r(t)$ in practice is the convolution of the reference THz pulse $i(t)$ and the assumed impulse response function $h_0(t)$, together with the dispersion factors. In order to construct the synthetic reflected THz signal $r(t)$, a double-Gaussian mixture model (DGMM) [29], which is expressed as

$$\beta_1 \exp\left(-\frac{(t - t_1)^2}{a_1^2}\right) + \beta_2 \exp\left(-\frac{(t - t_2)^2}{a_2^2}\right) \quad (7)$$

is implemented to estimate the original THz reference signal $i(t)$ first. Here, $\beta_{1,2}$ is the amplitude, $t_{1,2}$ is the initial location ($t_1 \neq t_2$), and $a_{1,2}$ is the bandwidth factor. The Levenberg-Marquardt algorithm is used to address the nonlinear least square procedure to obtain the optimal coefficients [30]. In order to qualify the evaluating performance of DGMM utilized, root-mean-square error (RMSE)

$$\text{RMSE} = \sqrt{\frac{1}{N} \sum_{n=1}^{N} [\hat{i}(n) - i(n)]^2} \quad (8)$$

is utilized. Here $N$ is the length of the THz reference signal $i(t)$, $\hat{i}$ and $i$ are the compensated THz signal and raw THz
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** (a) The comparison between the fitted result by DGMM and the measured THz reference signal $i(t)$. The power spectrum in (b) corresponding to the time-domain reference signals in (a) after Fourier transform.</figcaption></figure>

<table><caption>**Table 1.** Parametric vectors for the synthetic reflected THz echo components with dispersion.</caption><thead><tr><th>Parameters</th><th>Incident</th><th>atom1</th><th>atom2</th><th>atom3</th></tr></thead><tbody><tr><td>β<sub>1</sub></td><td>1.4296</td><td>0.4289</td><td>0.7148</td><td>0.1430</td></tr><tr><td>T</td><td>0</td><td>0</td><td>0.35</td><td>13.92</td></tr><tr><td>a<sub>1</sub></td><td>0.3203</td><td>0.3203</td><td>0.3203</td><td>1.6015</td></tr><tr><td>β<sub>2</sub></td><td>-0.4959</td><td>-0.1488</td><td>-0.2480</td><td>-0.0496</td></tr><tr><td>a<sub>2</sub></td><td>0.9147</td><td>0.9147</td><td>0.9147</td><td>4.5735</td></tr></tbody></table>

signal, respectively, and $n$ represents the index of the data points. It can be clearly observed in Figure 2 that DGMM fit the THz reference signal $i(t)$ well in both time and frequency domains, and the calculated RMSE is 0.0203. Minor mismatch at low ($<$0.1 THz) and high frequencies ($>$2.5 THz) might originate from the existence of noises near the upper and lower limits of the experimental system bandwidth.

The reflected THz signal $r(t)$ can be simulated using the fitted parameter vector, and the corresponding parametric vectors of each echo components are summarized in Table 1. Figure 3 shows the schematic diagram of the synthetic dispersive reflected THz signal. The time interval (optical delays) between the atom1 and atom2, and atom2 and atom3, which is associated with the thickness of layer I and II, are 0.35 ps and 13.57 ps, corresponding to an air gap with a thickness of 52.5 µm and 2.035 mm, respectively. Temporal pulse broadening associated with dispersion is accounted for atom3 at $t = 26.95$ ps. In order to simulate the reflected signal $r(t)$ obtained in the actual noisy environment, additive Gaussian white noise $e(t)$ with SNR = 5 dB, 10 dB, and 20 dB is included in the synthetic reflected signal $r(t)$, as shown in Figures 4a–4c.

Our first attempt to obtain the thicknesses of layers I and II employs FWDD. In this work, accounting for the noise present in the synthetic reflected signal $r(t)$, a Wiener

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Schematic diagram of sparse representation for synthetic dispersive reflected THz signal $r(t)$. The SNR of reflected signal $r(t)$ is set to 5 dB, 10 dB, and 20 dB, respectively.</figcaption></figure>

filter is applied with noising desensitizing factor 0.2 max ($|I(v)|^2$), where $I(v)$ is the Fourier transform of the reference signal $i(t)$ after compensation, and the sym4 wavelet function with decomposition level 5 is set for wavelet decomposition. Figures 4d–4f show the FWDD deconvolved result for the synthetic reflected signal $r(t)$. Because of the effective suppression of low- and high-frequency noise after band-passed frequency-domain filtering and stationary wavelet shrinkage process, the SNR of the deconvolved signal is significantly enhanced compared to the raw synthetic signal $r(t)$ (black). Owing to the narrowing of the bandwidth after frequency-domain filtering, however, the width of THz pulses in impulse response function $h(t)$ broadens after FWDD, leading to the failure of separating echoes reflected from layer I in the deconvolved result. In addition,
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Comparison between the raw synthetic reflected signal under various noise levels and the corresponding deconvolved results. (a)–(c) the synthetized reflected signals with 5 dB, 10 dB, and 20 dB, respectively. (d)–(f) the deconvolved results $h(t)$ by FWDD. (g)–(i) the deconvolved results $h(t)$ by adaptive SD/TWIST. Inset in (g)–(i) is the Zoom-in version of the deconvolved result between 12 ps and 16 ps.</figcaption></figure>

the position of resolved feature at ~26.94 ps have somewhat shifted, as shown in Figures 4d–4f, which might originate from the pulse spreading of atom3, demonstrating that the resolution in depth of FWDD is not high enough to retrieve the impulse response function $h(t)$ from dispersive THz signals reflected from complex structure. Note that FWDD is simply a linear filter.

To improve the axial resolution and successfully reconstruct the impulse response function $h(t)$, we now apply the adaptive SD method proposed. Superior to FWDD, adaptive SD is a pure time-domain technique that represents the impulse response function $h(t)$ by a few nonzero spikes. Because the selection of an appropriate regularization parameter is a crucial step to avoid the overfitting and underfitting, the fidelity item $||I'h - r||_2$ versus the preliminary regularization parameter $\lambda$ is employed to estimate the optimal $\lambda$ that guarantee the accuracy and efficiency of the reconstructed solution of the regularization problem. Figures 4g–4i present $h(t)$ reconstructed by adaptive SD method. Three peaks (Two positive echoes and one negative echo) are found, and the corresponding optical delay between the 1st and 2nd echo, and the 2nd and 3rd echo are 0.39 ps and 13.8 ps, confirming the significant enhancement in axial resolution.

Table 2 shows the CPU times of SD/TWIST, SD/IST, and SD/MM (one common SD method based on majorization–minimization optimization [31]). As expected, owing to the implementation of a nonlinear two-step iterative version of IST, the convergence rate of TWIST is clearly faster than SD/IST and SD/MM with the same parameters and stop criterion, $||h^{(k+1)} - h^k|| < \text{tolerance assumed to be } 10^{-6}$. For a raster scan over 10 cm by 10 cm area with a 1 mm step size, the corresponding number of pixels is

<table><caption>**Table 2.** Computational cost in seconds of three SD algorithms. @Intel(R) Core (TM) i7-9750H CPU @ 2.60 GHz.</caption><thead><tr><th>Algorithm</th><th>TWIST</th><th>MM</th><th>IST</th></tr></thead><tbody><tr><td>CPU time (s)</td><td>2.56</td><td>304</td><td>189.9</td></tr></tbody></table>

10,000, which might only take TWIST 35 min to analyze through parallelizing the computation into the 12 threads of a 6-core Intel i7 processor. A significant improvement of the computational speed is achieved to meet the requirements by practical applications.

## 3.2 Experimental verification

In the above numerical simulations, the proposed SD/TWIST algorithm exhibits better performance in characterizing specimens with complex structure than FWDD, and less computational cost than classic SD/IST. In this section, a paper stack comprised of standard A4 copy paper (Paperbox Brand multifunction white laser-printer and photocopy paper, manufactured by LECTA, Milton Keynes, UK [32]) is employed in the experiment. Note that air gaps of varying thickness also occur between the paper sheets.

In order to obtain the stratigraphic reconstruction of the paper stack studied, THz TOFT experiments were carried out, and a typical reflected signal $r(t)$ shown in Figure 5a. The largest pulse corresponds to the echo reflected from the surface of paper stack on the side from which the THz pulse is incident, and the subsidiary peaks at optical
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Comparison among the measured reflected THz signal $r(t)$ and the deconvolved signal $h(t)$ by FWDD and adaptive SD/TWIST. (a) Received temporal THz signal $r(t)$ from the paper stack; (b) Deconvolved signal $h(t)$ with FWDD; (c) Deconvolved signal $h(t)$ with adaptive SD/TWIST; (d) Schematic diagram of the seventeen-layered structure, which is composed of nine copy papers with eight air gaps between them.</figcaption></figure>

delay >10 ps corresponding to echoes from various subsequent interfaces. The relatively low local SNR makes it difficult to discriminate the effective interfacial features based on the raw reflected signal.

Adaptive SD/TWIST and FWDD were also implemented to extract the interface positions, and the corresponding reconstructed results are shown in Figures 5b and 5c. Even though a significant improvement in SNR is observed based on FWDD deconvolved results, several slight fluctuations associated with residual noises can still be found. Moreover, due to the presence of dispersion-induced pulse spreading in time-domain reflected THz signal, unexpected errors might be introduced when locating the echo position based on deconvolved signal by FWDD. Instead, the proposed adaptive SD method provides a more accurate structural characterization of the paper stack studied after incorporating DGMM dispersion model, and the broadened reflected echoes are not treated as a sum of several narrow pulses. The sign peak in Figure 5c is attributed to the refractive indices of successive adjacent layers. The thickness information of all individual copy paper in the paper stack can be expressed as [33]

$$d_{\text{Layer}} = \frac{\Delta t}{2} \frac{c}{n \cos \alpha}, \qquad (9)$$

where c is the speed of the light in vacuum, Δt corresponds to the optical delay of echoes reflected from copy layer, and the factor of one half arises since in reflection, the THz pulse passes through any given layer twice, α corresponds to the emitter angle, which is ~3°, n corresponds to the index of refraction of copy paper across the THz band. Note that multiple reflections occur during the propagation of THz waves in the paper stack are not

<table><caption>**Table 3.** The estimated physical thickness of copy layer in the paper stack from different methods.</caption><thead><tr><th>Sample</th><th>Nominal thickness (μm)</th><th>FWDD method (μm)</th><th>SD method (μm)</th></tr></thead><tbody><tr><td>Copy paper</td><td>102 ± 4</td><td>81 ± 7</td><td>98 ± 7</td></tr></tbody></table>

accounted for. The dielectric properties of the paper stack have been characterized using THz time-domain spectroscopy in ref. [32], and the index of refraction of copy paper varies little in the frequency range from 0.2 THz to 2.2 THz, and ~1.55 at 1 THz. Table 3 shows the thickness information of copy paper reconstructed from FWDD and adaptive SD/TWIST. The nominal physical thickness of copy paper is measured using caliper as a reference. Compared to the standard FWDD approach, adaptive SD/TWIST is more adept at locating dispersive THz pulses, demonstrating the effective and applicability of the proposed SD method in addressing reflected THz echoes from complex stratified structures.

## 4 Conclusions

In this study, an adaptive SD method has been employed to accurately reconstruct the stratigraphy of a paper stack. Balancing the computational cost against robustness, a two-step iterative shrinkage-thresholding (TWIST) algorithm is proposed to solve the residual ill-posed inverse problem. A double-Gaussian mixture model (DGMM) is also incorporated into the SD approach to compensate for the
broadening of THz pulses during the propagation of THz
waves in materials, which is associated with frequency-
dependent attenuation and dispersion. Both numerical sim-
ulations and experiments are implemented to validate the
efficiency and accuracy of SD/TWIST in characterizing
samples that are optically thick, but contain optically thin
layers. Compared to FWDD, the SD/TWIST provides a
more accurate and clearer representation of the impulse
response function $h(t)$ of the paper stack, with calculated
thicknesses of copy paper consistent with actual values.
Moreover, SD/TWIST presents a significant decrease in
the computation time compared with SD/IST. To con-
clude, SD/TWIST is shown to be a rapid and effective tool
in the structural characterization and layer thickness mea-
surement of multilayered structures nondestructively for
industrial applications.

Funding

The work has been supported financially Conseil Régional Grand Est.

Conflicts of Interest

The authors declare no conflict of interest.

Data availability statement

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

Author contribution statement

Conceptualization, A.L. and D.S.C.; methodology, A.L.; software, M.Z.; validation, A.L.; formal analysis, M.Z.; investigation, M.Z.; resources, M. Z. and D.S.C.; data curation, M.Z.; writing – original draft preparation, M.Z. and D.S.C.; writing – review and editing, A.L. and D.S.C.; visualization, M.Z.; supervision, D.S.C.; project administration, A.L. and D.S.C.; funding acquisition, D.S.C. All authors have read and agreed to the published version of the manuscript.

References

1 Zhong S., Nsengiyumva W. (2022) Ultrasonic testing techniques for nondestructive evaluation of fiber-reinforced composite structures, *Nondestructive Testing and Evaluation of Fiber-Reinforced Composite Structures*, Springer, Singapore, pp. 133–195.

2 Zhai M., Locquet A., Roquelet C., Alexandre P., Daheron L., Citrin D.S. (2020) Nondestructive measurement of mill-scale thickness on steel by terahertz time-of-flight tomography, *Surf. Coat. Technol.* **393**, 125765.

3 Xu Y., Jiang X. (2022) Nondestructive testing and imaging of corrosion thickness of steel plates using THz-TDS, *Infrared Phys. Technol.* **127**, 104467.

4 Zhai M., Locquet A., Roquelet C., Borean J.L., Meilland P., Citrin D. S. (2023) *Nondestructive tertiary mill-scale thickness measurement on commercial hot-rolled steel strip: terahertz time-of-flight tomography*, Steel Research International.

5 Krugener K., Ornik J., Schneider L.M., Jäckel A., Koch-Dandolo C. L., Castro-Camus E., Riedl-Siedow N., Koch M., Viöl W. (2020) Terahertz inspection of buildings and architectural art, *Appl. Sci.* **10**, 15, 5166.

6 Fukunaga K. (2023) Nondestructive evaluation of lined paintings by THz pulsed time-domain imaging, *Heritage* **6**, 4, 3448–3460.

7 Tornari V., Andrianakis M., Duchêne S., Nowik W., Brissaud D., Giovannacci D., Küppers M., Rehorn C., Blümich B., Ricci G., Artioli G. (2023) A combined ND diagnostic investigation by DHSPI, SIRT, NMR, THZ, on Giotto fresco, *J. Cult. Herit.* **63**, 206–216.

8 Zhai M., Locquet A., Citrin D.S. (2020) Pulsed THz imaging for thickness characterization of plastic sheets, *NDT&E Int.* **116**, 102338.

9 Xu Y., Hao H., Citrin D.S., Wang X., Zhang L., Chen X. (2021) Three-dimensional nondestructive characterization of delamination in GFRP by terahertz time-of-flight tomography with sparse Bayesian learning-based spectrum-graph integration strategy, *Compos. Part B Eng.* **225**, 109285.

10 Lu X., Shen Y., Xu T., Sun H., Zhu L., Zhang J., Chang T., Cui H.L. (2022) Accurate detectyion of porosity in glass fiber reinforced polymers by terahertz spectroscopy, *Compos. Part B Eng.* **242**, 110058.

11 Calvo-de la Rosa J., Pomarède P., Antonik P., Mertaghni F., Citrin D.S., Rontani D., Locquet A. (2023) Determination of the process-induced microstructure of woven glass fabric reinforced polyamide 6.6/6 composite using terahertz pulsed imaging, *NDT&E Int.* **136**, 102799.

12 P.H. Siegel, Therahertz technology in biology and medicine, IEEE Trans. Microw. Theory Tech. 52(10), 2438-2447.

13 Brun M.A., Formanek F., Yasuda A., Sekine M., Ando N., Eishii Y. (2010) Terahertz imaging applied to cancer diagnosis, *Phys. Med. Biol.* **55**, 16, 4615.

14 Zaytsev K.I., Dolganova I.N., Chernomyrdin N.V., Katyba G.M., Gavdush A.A., Cherkasova O.P., Komandin G.A., Shchedrina M.A., Khodan A.N., Ponomarev D.S., Reshetov I.V., Karasik V.E., Skorobogatiy M., Kurlov V.N., Tuchin V.V. (2019) The progress and perspectives of terahertz technology for diagnosis of neoplasms: a review, *J. Optics* **22**, 1, 013001.

15 Davies A.G., Burnett A.D., Fan W., Linfield E.H., Cunningham J.E. (2008) Terahertz spectroscopy of explosives and drugs, *Mater. Today* **11**, 3, 18–26.

16 Yan Z., Shi W. (2022) Detection of aging in the common explosive RDX using terahertz time-domain spectroscopy, *J. Opt. Soc. Am. B* **39**, 3, A9–A12.

17 Sharma M., Sharma B., Gupta A.K., Pandey D. (2023) Recent developments of image processing to improve explosive detection methodologies and spectroscopic imaging techniques for explosive and drug detection, *Multimed. Tools. Appl.* **82**, 6849–6865.

18 Mittleman D.M., Hunsche S., Boivin L., Nuss M.C. (1997) T-ray tomography, *Opt. Lett.* **23**, 12, 904–906.

19 Chen Y., Huang S., Pickwell-MacPherson E. (2010) Frequency-wavelet domain deconvolution for terahertz reflection imaging and spectroscopy, *Opt. Exp.* **18**, 2, 1177–1190.

20 Dong J., Locquet A., Citrin D.S. (2017) Depth resolution enhance-
ment of terahertz deconvolution by autoregressive spectral extrap-
olation, *Opt. Lett.* **429**, 1828–1831.

21 Zhai M., Locquet A., Roquelet C., Citrin D.S. (2020) Terahertz Time-of-flight tomography beyond the axial resolution limit: Autoregressive spectral estimation based on the modified covariance method, *J. Infrared Millim. Terahertz Waves* **41**, 926–939.

22 Dong J., Wu X., Locquet A., Citrin D.S. (2017) Terahertz superres-
olution stratigraphic characterization of multilayered structures
using sparse deconvolution, *IEEE Trans. Terahertz Sci. Technol.*
**7**, 3, 260–267.

23 Chang Y., Zi Y., Zhao J., Yang Z., He W., Sun H. (2017) An adaptive sparse deconvolution method for distinguishing the overlapping echoes of ultrasonic guided waves for pipeline crack inspection, *Meas. Sci. Technol.* **28**, 035002.

24 Zhai M., Citrin D.S., Locquet A. (2021) Terahertz nondestructive stratigraphic analysis of complex layered structures: reconstruction techniques, *J. Infrared Millim. Terahertz Waves* **42**, 929–946.

25 Huang Y., Sun P., Zhang Z., Jin C. (2017) Numerical method based on transfer function for eliminating water vapor noise from terahertz spectra, *Appl. Opt.* **56**, 20, 5698–5704.

26 Neu J., Schmuttenmaer C.A. (2018) Tutorial: An introduction to terahertz time domain spectroscopy, *J. Appl. Phys.* **124**, 23, 231101.

27 Bioucas-Dias J.M., Figueiredo M.A.T. (2007) A new TwIST: Two-step iterative shrinkage/thresholding algorithms for image restoration, *IEEE Trans. Image Process.* **16**, 12, 2992–3004.
28 Huang C., Ji H., Qiu J., Wang L., Wang X. (2022) TWIST sparse regularization method using cubic B-spline dual scaling functions for impact force identification, *Mech. Syst Signal Process.* **167**, 108451.

29 Xu Y., Fang X., Fan S., Zhang L., Yan R., Chen X. (2022) Double Gaussian mixture model-based terahertz wave dispersion compensation method using convex optimization technique, *Mech. Syst. Signal Process.* **164**, 108223.

30 Moré J.J. (1977) The Levenberg-Marquardt algorithm: implementation and theory, in *Numerical analysis: Proceedings of the Biennial Conference Held at Dundee, June 28-July 1, 1977*, Springer, Berlin, Heidelberg, pp. 105–116.

31 Abdessalem B., Farid C. (2020) Resolution improvement of ultrasonic signals using sparse deconvolution and variational mode decomposition algorithms, *Russ. J. Nondestruct. Test.* **56**, 6, 479–489.

32 Zhai M., Locquet A., Citrin D.S. (2021) Terahertz imaging for paper handling of legacy documents, *Sensors* **21**20, 6756.

33 Vandrevala F., Einarsson E. (2018) Decoupling substrate thickness and refractive index measurement in THz time-domain spectroscopy, *Opt. Exp.* **26**, 2, 1697–1702.