© The Author(s), published by EDP Sciences, 2024

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 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₁-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.

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 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 π 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₂O molecules during the propagation of THz pulses [25]. Measurements can also be made in a dry atmosphere of N₂; 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)⨂h\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₁-norm regularization optimization, which is defined as$$\begin{array}{c}\underset{h}{\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||₁ denotes ${\sum }_i^N\left|{\mathit{h}}_{\mathit{i}}\right|$, N is the length of the vector, and λ is the non-negative regularization parameter that controls the tradeoff between the sparsity of h and the residue norm.

The iterative-shrinkage thresholding (IST) algorithm is an effective method to solve the l₁-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$$\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}}_{k-1}+\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₀ 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 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₀(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{{\left(t-t_1\right)}^2}{a_1^2}\right)+{\beta }_2\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₁ ≠ t₂), 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)$$\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 time-domain 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 Zoom-in version of the deconvolved result between 12 ps and 16 ps.

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(ν)|²), 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 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, 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||₂ 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 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 || < tolerance assumed to be 10⁻⁶. 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 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 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 seventeen-layered 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 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]$$\begin{array}{c}{d}_{\mathrm{Layer}}=\frac{∆t}{2}\frac{c}{n \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 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 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

All Tables

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 time-domain 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 Zoom-in 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 seventeen-layered structure, which is composed of nine copy papers with eight air gaps between them.

In the text