Issue 
J. Eur. Opt. SocietyRapid Publ.
Volume 20, Number 1, 2024
THz imaging



Article Number  4  
Number of page(s)  8  
DOI  https://doi.org/10.1051/jeos/2024001  
Published online  08 March 2024 
Research Article
Terahertz nondestructive stratigraphic reconstruction of paper stacks based on adaptive sparse deconvolution
^{1}
Georgia TechCNRS IRL2958, Georgia TechEurope, 2 Rue Marconi, 57070 Metz, France
^{2}
School of Electrical and Computer Engineering, Georgia Institute of Technology, Atlanta, GA 303320250, USA
^{*} Corresponding author: david.citrin@ece.gatech.edu
Received:
14
September
2023
Accepted:
18
January
2024
Characterizing the number of sheets in a stack of paper typically involves mechanical separation of the individual sheets. Here, we explore an nondestructive method that can be applied to the intact paper stack. Namely, terahertz timeofflight tomography, together with post signalprocessing technique sparse deconvolution based on a twostep iterative shrinkagethresholding 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 doubleGaussian mixture model (DGMM) is also incorporated to suppress dispersion in the reflected THz echoes. The effectiveness and accuracy of the proposed adaptive sparsedeconvolution method are verified experimentally and numerically. Compared with the commonly used frequency waveletdomain deconvolution (FWDD) method and previous implementations of sparse deconvolution based on an iterativeshrinkage and thresholding algorithm (SD/IST), the proposed sparsedeconvolution 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 realtime THz tomographicimage processing.
Key words: Terahertz timedomain spectroscopy / Paper handling / Superresolution characterization / Dispersion / Sparse deconvolution
© The Author(s), published by EDP Sciences, 2024
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
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 onestep 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 Xrays enabling the characterization of structure interior to various objects. Moreover, in contrast to Xray 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 timeofflight 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 invacuum 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, pulsespreading associated with frequencydependent attenuation and dispersion during the propagation of THz waves in materials, as well as low signaltonoise 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 signalprocessing approach to nondestructively characterize layer thicknesses in complex multilayer structures in the presence of material dispersion. The test case on which we focus are stacks of paper.
To date, several advanced signalprocessing methods have been proposed and employed to analyze the complex reflected THz echoes, such as frequency waveletdomain deconvolution (FWDD) [19] and autoregressive (AR) spectral extrapolation [20, 21]. Based on the assumption that the reflected THz echoes are the timeshifted amplitudescaled replicas of the THz incident signals, both FWDD and AR methods are mainly limited to deconvolve reflected signals with timeinvariant 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 timedomain 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 timeconsuming 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 twostep iterative shrinkagethresholding algorithm (SD/TWIST) is proposed and employed to quantitatively reconstruct the stratigraphy of a multilayer 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 timevarying 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 signalprocessing methods, such as FWDD. A clearer stratigraphic analysis with limited erroneous spikelike 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. Quasisinglecycle 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 Erdoped fiber which emits 780 nm pulses with sub100 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 timegated output signal. The corresponding spectrum is obtained after Fourier transform.
Fig. 1 Schematic diagram of the TPS Spectra 3000 apparatus in reflection mode from TeraView Ltd., Cambridge, UK. 
In this study, THz measurements were carried out in an airconditioned 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 π 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 rovibrational 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,$$\begin{array}{c}r\left(t\right)=i\left(t\right)\u2a02h\left(t\right)=\underset{\infty}{\overset{\infty}{\int}}i\left(t\tau \right)h\left(\tau \right)\mathrm{d}\tau ,\end{array}$$(1)where the symbol ⨂ denotes the convolution operation and τ is the time delayed operation satisfying t ≥ τ. Accounting for the noise e(t) originating from the measurement [26], the discrete format of equation (1) can be written as a matrix multiplication,$$\begin{array}{c}\mathit{r}=\mathit{Ih}+\mathit{e},\end{array}$$(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$$\begin{array}{c}\underset{h}{\mathrm{min}}\frac{1}{2}{\left\right\mathit{Ih}\mathit{r}\left\right}_{2}^{2}+\lambda {\left\right\mathit{h}\left\right}_{1},\end{array}$$(3)where ${\left\right\mathit{I\text{'}h}\mathit{r}\left\right}_{2}^{2}$ denotes ${\sum}_{i}^{N}{(\mathit{I}{\mathit{h}}_{\mathit{i}}{\mathit{r}}_{\mathit{i}})}^{2}$, h_{1} denotes ${\sum}_{i}^{N}\left{\mathit{h}}_{\mathit{i}}\right$, N is the length of the vector, and λ is the nonnegative regularization parameter that controls the tradeoff between the sparsity of h and the residue norm.
The iterativeshrinkage 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 twostep iterative shrinkagethresholding algorithm (TWIST) is utilized to retrieve h, which can be expressed as$$\begin{array}{c}{\mathit{h}}_{1}={S}_{T}\left({\mathit{h}}_{0}+{\mathit{I}}^{T}\left(\mathit{r}\mathit{I}{\mathit{h}}_{0}\right)\right),\end{array}$$(4) $$\begin{array}{c}{\mathit{h}}_{k+1}={\left(1\alpha \right){\mathit{h}}_{k1}+\left(\alpha \beta \right){\mathit{h}}_{k}+\beta S}_{T}({\mathit{h}}_{k}+{\mathit{I}\mathbf{\text{'}}}^{T}(\mathit{r}\mathit{I}{\mathit{h}}_{k}\left)\right)\end{array}$$(5)where h_{0} is the initial vector, α and β are two convergence parameters of TWIST, which are associated with the eigenvalues of matrix multiplication I^{T} I. In addition,$${S}_{T}\left(\mathit{h}\right)=\{\begin{array}{cc}\mathrm{sgn}\left({h}_{i}\right)\left(\left{h}_{i}\right\lambda \right),& \left{h}_{i}\right\lambda \\ 0,& \mathrm{otherwise}\end{array},$$(6)where S_{T}(h) is a softthresholding 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 doubleGaussian mixture model (DGMM) [29], which is expressed as$${\beta}_{1}\mathrm{exp}\left(\frac{{\left(t{t}_{1}\right)}^{2}}{{a}_{1}^{2}}\right)+{\beta}_{2}\mathrm{exp}\left(\frac{{\left(t{t}_{2}\right)}^{2}}{{a}_{2}^{2}}\right)$$(7)is implemented to estimate the original THz reference signal i(t) first. Here, β_{{1,2}} is the amplitude, t_{{1,2}} is the initial location (t_{1} ≠ 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, rootmeansquare error (RMSE)$$\begin{array}{c}\mathrm{RMSE}=\sqrt{\frac{1}{N}\sum _{n=1}^{N}{[\widehat{i}\left(n\right)i\left(n\right)]}^{2}}\end{array}$$(8)is utilized. Here N is the length of the THz reference signal i(t), $\widehat{i}$ and $i$ are the compensated THz signal and raw THz 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.
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 timedomain reference signals in (a) after Fourier transform. 
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.
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. 
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 Zoomin version of the deconvolved result between 12 ps and 16 ps. 
Parametric vectors for the synthetic reflected THz echo components with dispersion.
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 filter is applied with noising desensitizing factor 0.2 max(I(ν)^{2}), where I(ν) 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 highfrequency noise after bandpassed frequencydomain 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 frequencydomain 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, 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 timedomain 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 λ is employed to estimate the optimal λ 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 twostep 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 } < 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 10,000, which might only take TWIST 35 min to analyze through parallelizing the computation into the 12 threads of a 6core Intel i7 processor. A significant improvement of the computational speed is achieved to meet the requirements by practical applications.
Computational cost in seconds of three SD algorithms. @Intel(R) Core (TM) i79750H CPU @ 2.60 GHz.
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 laserprinter 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 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.
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 seventeenlayered structure, which is composed of nine copy papers with eight air gaps between them. 
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 dispersioninduced pulse spreading in timedomain 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]$$\begin{array}{c}{d}_{\mathrm{Layer}}=\frac{\u2206t}{2}\frac{c}{n\mathrm{cos}\alpha},\end{array}$$(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 accounted for. The dielectric properties of the paper stack have been characterized using THz timedomain 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.
The estimated physical thickness of copy layer in the paper stack from different methods.
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 twostep iterative shrinkagethresholding (TWIST) algorithm is proposed to solve the residual illposed inverse problem. A doubleGaussian 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 frequencydependent attenuation and dispersion. Both numerical simulations 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 conclude, SD/TWIST is shown to be a rapid and effective tool in the structural characterization and layer thickness measurement 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 corresponding 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
 Zhong S., Nsengiyumva W. (2022) Ultrasonic testing techniques for nondestructive evaluation of fiberreinforced composite structures, Nondestructive Testing and Evaluation of FiberReinforced Composite Structures, Springer, Singapore, pp. 133–195. [CrossRef] [Google Scholar]
 Zhai M., Locquet A., Roquelet C., Alexandre P., Daheron L., Citrin D.S. (2020) Nondestructive measurement of millscale thickness on steel by terahertz timeofflight tomography, Surf. Coat. Technol. 393, 125765. [CrossRef] [Google Scholar]
 Xu Y., Jiang X. (2022) Nondestructive testing and imaging of corrosion thickness of steel plates using THzTDS, Infrared Phys. Technol. 127, 104467. [NASA ADS] [CrossRef] [Google Scholar]
 Zhai M., Locquet A., Roquelet C., Borean J.L., Meilland P., Citrin D.S. (2023) Nondestructive tertiary millscale thickness measurement on commercial hotrolled steel strip: terahertz timeofflight tomography, Steel Research International. [Google Scholar]
 Krugener K., Ornik J., Schneider L.M., Jäckel A., KochDandolo C.L., CastroCamus E., RiedlSiedow N., Koch M., Viöl W. (2020) Terahertz inspection of buildings and architectural art, Appl. Sci. 10, 15, 5166. [CrossRef] [Google Scholar]
 Fukunaga K. (2023) Nondestructive evaluation of lined paintings by THz pulsed timedomain imaging, Heritage 6, 4, 3448–3460. [CrossRef] [Google Scholar]
 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. [CrossRef] [Google Scholar]
 Zhai M., Locquet A., Citrin D.S. (2020) Pulsed THz imaging for thickness characterization of plastic sheets, NDT&E Int. 116, 102338. [Google Scholar]
 Xu Y., Hao H., Citrin D.S., Wang X., Zhang L., Chen X. (2021) Threedimensional nondestructive characterization of delamination in GFRP by terahertz timeofflight tomography with sparse Bayesian learningbased spectrumgraph integration strategy, Compos. Part B Eng. 225, 109285. [CrossRef] [Google Scholar]
 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. [CrossRef] [Google Scholar]
 Calvode la Rosa J., Pomarède P., Antonik P., Mertaghni F., Citrin D.S., Rontani D., Locquet A. (2023) Determination of the processinduced microstructure of woven glass fabric reinforced polyamide 6.6/6 composite using terahertz pulsed imaging, NDT&E Int. 136, 102799. [Google Scholar]
 P.H. Siegel, Therahertz technology in biology and medicine, IEEE Trans. Microw. Theory Tech. 52(10), 2438–2447. [Google Scholar]
 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. [NASA ADS] [CrossRef] [Google Scholar]
 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. [Google Scholar]
 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. [CrossRef] [Google Scholar]
 Yan Z., Shi W. (2022) Detection of aging in the common explosive RDX using terahertz timedomain spectroscopy, J. Opt. Soc. Am. B 39, 3, A9–A12. [CrossRef] [Google Scholar]
 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. [CrossRef] [Google Scholar]
 Mittleman D.M., Hunsche S., Boivin L., Nuss M.C. (1997) Tray tomography, Opt. Lett. 23, 12, 904–906. [NASA ADS] [CrossRef] [Google Scholar]
 Chen Y., Huang S., PickwellMacPherson E. (2010) Frequencywavelet domain deconvolution for terahertz reflection imaging and spectroscopy, Opt. Exp. 18, 2, 1177–1190. [NASA ADS] [CrossRef] [Google Scholar]
 Dong J., Locquet A., Citrin D.S. (2017) Depth resolution enhancement of terahertz deconvolution by autoregressive spectral extrapolation, Opt. Lett. 42 9, 1828–1831. [NASA ADS] [CrossRef] [Google Scholar]
 Zhai M., Locquet A., Roquelet C., Citrin D.S. (2020) Terahertz Timeofflight tomography beyond the axial resolution limit: Autoregressive spectral estimation based on the modified covariance method, J. Infrared Millim. Terahertz Waves 41, 926–939. [NASA ADS] [CrossRef] [Google Scholar]
 Dong J., Wu X., Locquet A., Citrin D.S. (2017) Terahertz superresolution stratigraphic characterization of multilayered structures using sparse deconvolution, IEEE Trans. Terahertz Sci. Technol. 7, 3, 260–267. [CrossRef] [Google Scholar]
 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. [NASA ADS] [CrossRef] [Google Scholar]
 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. [NASA ADS] [CrossRef] [Google Scholar]
 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. [NASA ADS] [CrossRef] [Google Scholar]
 Neu J., Schmuttenmaer C.A. (2018) Tutorial: An introduction to terahertz time domain spectroscopy, J. Appl. Phys. 124, 23, 231101. [NASA ADS] [CrossRef] [Google Scholar]
 BioucasDias J.M., Figueiredo M.A.T. (2007) A new TwIST: Twostep iterative shrinkage/thresholding algorithms for image restoration, IEEE Trans. Image Process. 16, 12, 2992–3004. [CrossRef] [Google Scholar]
 Huang C., Ji H., Qiu J., Wang L., Wang X. (2022) TWIST sparse regularization method using cubic Bspline dual scaling functions for impact force identification, Mech. Syst Signal Process. 167, 108451. [CrossRef] [Google Scholar]
 Xu Y., Fang X., Fan S., Zhang L., Yan R., Chen X. (2022) Double Gaussian mixture modelbased terahertz wave dispersion compensation method using convex optimization technique, Mech. Syst. Signal Process. 164, 108223. [NASA ADS] [CrossRef] [Google Scholar]
 Moré J.J. (1977) The LevenbergMarquardt 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. [Google Scholar]
 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. [CrossRef] [Google Scholar]
 Zhai M., Locquet A., Citrin D.S. (2021) Terahertz imaging for paper handling of legacy documents, Sensors 21 20, 6756. [NASA ADS] [CrossRef] [Google Scholar]
 Vandrevala F., Einarsson E. (2018) Decoupling substrate thickness and refractive index measurement in THz timedomain spectroscopy, Opt. Exp. 26, 2, 1697–1702. [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Parametric vectors for the synthetic reflected THz echo components with dispersion.
Computational cost in seconds of three SD algorithms. @Intel(R) Core (TM) i79750H CPU @ 2.60 GHz.
The estimated physical thickness of copy layer in the paper stack from different methods.
All Figures
Fig. 1 Schematic diagram of the TPS Spectra 3000 apparatus in reflection mode from TeraView Ltd., Cambridge, UK. 

In the text 
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 timedomain reference signals in (a) after Fourier transform. 

In the text 
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. 

In the text 
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 Zoomin version of the deconvolved result between 12 ps and 16 ps. 

In the text 
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 seventeenlayered structure, which is composed of nine copy papers with eight air gaps between them. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.