1.

Introduction

Terahertz (THz) radiation is located between the infrared and microwave regions of the electromagnetic spectrum, from about 0.1 to 10.0 THz. THz waves have some specific properties, such as high penetration depth into dielectric materials and the ability to interact with both vibrational and rotational molecular levels. THz radiation is nonionizing in nature. Because of these properties, there are many potential applications of THz technology.¹ For instance, THz imaging and spectroscopy systems are used in the detection of concealed weapons, drugs, and explosives,^2,3 medical diagnostics,^4–6 and nondestructive evaluation of construction materials.^7,8

THz time-domain spectroscopy (TDS) utilizes short pulses of THz radiation with a wideband spectrum (from 0.1 to 3.0 THz) to measure the THz optical properties of various media. The electric field of a THz pulse is detected with a very high time resolution (about 50.0 fs) after the transmission of the pulse through the sample or the reflection of the pulse from the surface. By applying a fast Fourier transformation to the detected time-domain signals, it is possible to analyze the complex transmission or reflection coefficients of the sample in a wide frequency range. Time-domain data can also be used to reconstruct the permittivity profile in a sample (i.e., the dependence of the sample dielectric permittivity on depth). This reconstruction method aids in the study of the internal structure of a sample, known as THz tomography or T-ray tomography.⁹ THz tomography is expected to be useful in medical diagnosis, nondestructive evaluation of construction materials, inspection of art objects, and other applications. This paper describes a new way to solve this inverse problem.

Many ways have been presented to solve the inverse scattering problem.^9,10 However, several aspects of the inverse problem have not previously been considered, such as sample pulse function low-frequency interpolation, and sample pulse function filtering during the deconvolution process. The reconstruction results significantly depend on the particular method of solution for these problems. In addition, the stability of different algorithms with respect to noise has not been previously studied. The present work does not critically analyze previous solutions to the inverse problem. Instead, in this paper, an alternative way to solve the inverse scattering problem is presented. The results of implementing this new algorithm and its stability with respect to noise are shown in detail.

Several assumptions concerning media properties are made in the present work. The medium of interest is assumed to be nondispersive and nonabsorbing. Sample optical properties are assumed to be constant in the lateral direction and vary only with depth. The media of interest must possess no optical nonlinearity. When the THz penetration depth is much larger than the full optical width, a medium is considered to be low-absorbing. The dispersion of optical properties can be neglected when the variation of the medium refractive index has small fluctuations with frequency. The medium is assumed to be thick such that multiple reflections in sample layers do not impact the reflected signal data (i.e., the number of pulses in the sample response). All media considered in this paper are assumed to satisfy these assumptions.

Restricting our analysis to nondispersive, nonabsorbing media, we apply the present algorithm for nondestructive evaluation of polymer construction materials. Future modifications of this method would help generalize this technique to absorbing media and media with dispersion. For example, dispersion can be described with a spectral dielectric permittivity model, such as the Debye, Lawrence, or Drude models. This improvement would aid in the development of new methods of medical diagnosis with THz TDS systems (diagnosis of skin burns,¹¹ study of skin structure,¹² study of tooth demineralization,¹³ and cornea tissue¹⁴) and methods of art inspection.¹⁵

Most inverse scattering algorithms used in different electromagnetic spectral ranges, such as in the radio frequency range, consist of two steps for sample permittivity profile reconstruction. First, a sample pulse function $R(t)$ is reconstructed based on time-domain signals. Second, the sample permittivity profile $\epsilon (z)$ is reconstructed utilizing the pulse function $R(t)$. In this paper, an invariant embedding technique is suggested as the basis for the permittivity profile reconstruction procedure. This technique was previously used outside of the THz range.^16,17

Several problems should be solved in order to generalize this technique to THz TDS. These problems include the low- and high-frequency noise filtering of the sample pulse function, sample pulse function interpolation at low frequencies, and the correction of the dielectric permittivity profile using a priori information. Permittivity profile reconstruction based on an invariant embedding technique was implemented and tested experimentally and using mathematical modeling. We consider all stages of permittivity profile reconstruction using this algorithm. We discuss all methods for solving the problems of invariant embedding technique generalization to THz spectroscopy.

2.

Materials and Methods

2.1.

THz Time-Domain Spectrometer

First, a THz TDS system used to collect experimental data is described. Figure 1 presents a THz time-domain spectrometer working in reflection mode. Generation and detection of THz pulses is made using optical femtosecond pulses of a fiber laser with full width at half maximum of about 80.0 fs. The generation of THz pulses is produced in a photoconductive (PC) antenna. The laser repetition rate is 50.0 MHz and the THz pulses are modulated at 102.4 kHz using PC voltage modulation.¹⁸ The THz beam is focused on the sample. Moreover, the aperture ratio of lens is rather small: $D/f^{\prime }=1:2.5$. Rays in THz beam transmitted through the lens have different incident angles. The maximum incident angle is equal to the aperture angle of the lens. Since the aperture angle is rather small, then the incident angle is also small. This condition is an important one for correct implementation of the algorithm. THz radiation reflects from the surface of the sample during sample signal measurements or from the reference surface during reference signal registration and is incident on the electrooptical THz field detector based on ZnTe crystal. Using a probe femtosecond optical pulse transmitted through the delay stage, the THz electrical field magnitude is detected with very high time resolution. A lock-in amplifier is synchronized with the PC antenna modulating voltage to provide a high signal-to-noise ratio.

Fig. 1

Schematic representation of the TDS system.

Two signals need to be obtained for sample permittivity profile reconstruction: $E_s(t)$ is the signal reflected from the sample of interest, and $E_r(t)$ is the reference signal reflected from the reference surface with a very high and homogeneous reflectivity in a wide spectral range. A planar gold mirror was used to obtain the reference signal. As the method of signal detection is known, all stages of dielectric permittivity profile reconstruction can be considered.

2.2.

Dielectric Permittivity Profile Reconstruction Algorithm

As previously noted, sample dielectric permittivity profile reconstruction is accomplished in two steps: reconstruction of the sample pulse function $R(t)$ based on TDS system signals $E_s(t)$ and $E_r(t)$ and the reconstruction of the sample permittivity profile $\epsilon (z)$ based on a sample pulse function $R(t)$.

2.2.1.

Sample pulse function reconstruction

The sample pulse function $R(t)$ is the medium response to the influence of an incident electromagnetic wave with a delta function amplitude time dependence. The sample pulse function is also the result of an inverse Fourier transformation applied to the sample amplitude reflectivity $\tilde{R}({\nu }_t)$. Pulse function reconstruction is a signal deconvolution process.

The first problem in signal deconvolution is that simply dividing the sample signal Fourier spectrum ${\tilde{E}}_s({\nu }_t)$ by the base signal Fourier spectrum ${\tilde{E}}_r({\nu }_t)$ to reconstruct the sample complex amplitude reflectivity $\tilde{R}({\nu }_t)$

Eq. (1)

$$\tilde{R}({\nu }_t)=\frac{{\tilde{E}}_s({\nu }_t)}{{\tilde{E}}_r({\nu }_t)}$$

does not work. Since both the sample and reference signals contain no useful information at low frequencies ($<0.05$–0.1 THz) and high frequencies ($>2.0$–3.0 THz), Eq. (1) only results in noise in these spectral regions. Any sample reflectivity reconstruction procedure must contain some low- and high-frequency noise-filtering operation. The filtering procedure should allow us to obtain as much data on the sample reflectivity as possible, and should contain smoothed filter edges to prohibit the emergence of Gibbs noise in the sample pulse function. With a filter, Eq. (1) is modified as

Eq. (2)

$$\tilde{R}({\nu }_t)=\frac{{\tilde{E}}_s({\nu }_t)}{{\tilde{E}}_r({\nu }_t)}\tilde{F}({\nu }_t).$$

A technique utilizing a Wiener filter $\tilde{F}({\nu }_t)$ (see Ref. 19) was applied in previous research^20,21 to process TDS system data. The formulation of a derived Wiener filter $\tilde{F}({\nu }_t)$ with several modifications can be described as follows:

Eq. (3)

$$\tilde{F}({\nu }_t)=\frac{{\left[\frac{{\tilde{E}}_r({\nu }_t)}{\underset{{\nu }_t}{\max }\{{\tilde{E}}_r({\nu }_t)\}}\right]}^2}{{\left[\frac{{\tilde{E}}_r({\nu }_t)}{\underset{{\nu }_t}{\max }\{{\tilde{E}}_r({\nu }_t)\}}\right]}^2+\frac{N({\nu }_t)}{S({\nu }_t)}},$$

where $N({\nu }_t)$ is a noise power spectrum model and $S({\nu }_t)$ is a signal power spectrum model. Equation (3) is negligible in a frequency range where the noise power spectrum model $N({\nu }_t)$ exceeds the signal power spectrum model $S({\nu }_t)$. Conversely, if $S({\nu }_t)$ exceeds $N({\nu }_t)$, the function $\tilde{F}({\nu }_t)$ is close to unity. Smooth filter edges are also provided. For simplicity, we assume a Gaussian white noise model for the noise power spectrum $N({\nu }_t)$:

Eq. (4)

$$N({\nu }_t)=\tilde{K},$$

where $\tilde{K}$ is a constant such that $0.0<\tilde{K}<1.0$. The signal power spectrum $S({\nu }_t)$ is based on a Gaussian monopulse function $f_{\text{GMp}}(t)$:

Eq. (5)

$$f_{\text{GMp}}(t)=-2·e^{-\frac{1}{2}}·(\pi ·{\nu }_{tC}·t)·\exp [-2{(\pi ·{\nu }_{tC}·t)}^2],$$

where ${\nu }_{tC}$ is the most powerful harmonic in the monopulse spectrum ${\tilde{f}}_{\mathrm{GMp}}({\nu }_t)$. The value ${\nu }_{tC}$ depends on the methods that are used for generation and detection of THz pulses. The signal power spectrum model can be obtained using the following equation:

Eq. (6)

$$S({\nu }_t)=\frac{{[{Ϝ}_t\{f_{\text{GMp}}(t)\}]}^2}{\underset{{\nu }_t}{\max }\{{[{Ϝ}_t\{f_{\text{GMp}}(t)\}]}^2\}},$$

where ${Ϝ}_t$ is the direct Fourier transform operator.

Choosing a central monopulse frequency ${\nu }_{tC}$ and noise spectral power $\tilde{K}$ in Eqs. (3)–(6), it is possible to provide a qualitative reconstruction of the sample reflectivity for different conditions of waveform $E_r(t)$ and $E_s(t)$ registration, such as the type of THz emitter and detector, the number of signal averages, and other factors. Figure 2 demonstrates an example of a Wiener filter function $\tilde{F}({\nu }_t)$ generated for a typical reference signal $E_r(t)$ with filter parameters: ${\nu }_{tC}=0.8\mathrm{THz}$, $\tilde{K}=0.002$. Luckily for our TDS system, the center frequency ${\nu }_{tC}$ matches a window in the atmospheric water vapor absorption spectrum. Outside this window the filtering procedure obviously leads to noise suppression in frequency ranges corresponding to high atmospheric absorption.

Fig. 2

An example of reference spectrum ${\tilde{E}}_r({\nu }_t)$ and Wiener filter $\tilde{F}({\nu }_t)$

The second problem in the deconvolution procedure is pulse function interpolation at low frequencies. To reconstruct a sample permittivity profile with the invariant embedding technique, it is essential to know the sample pulse function spectral components at low frequencies, excluding the zero point of the discrete Fourier-domain. It is possible to interpolate the sample reflectivity $\tilde{R}({\nu }_t)$ based on known frequency domain information within the ranges $[-{\nu }_{t\_\max },-{\nu }_{t\_\min }]$ and $[{\nu }_{t\_\min },{\nu }_{t\_\max }]$, and to obtain low-frequency data from the interpolation results. Since the pulse function spectrum $\tilde{R}({\nu }_t)$ character is determined both by the internal structure of the sample and by the dispersion of media optical properties, interpolation can be easily accomplished if the investigated medium has a negligible dispersion of optical characteristics.

The modulus $|\tilde{R}({\nu }_t)|$ and phase $\phi \{\tilde{R}({\nu }_t)\}$ of the sample complex reflectivity can be interpolated separately. It is convenient to use a trigonometric interpolation series to reconstruct the amplitude reflectivity modulus $|\tilde{R}({\nu }_t)|$. Since the sample pulse function $R(t)$ is a real function, the series has the form:

Eq. (7)

$$|\tilde{R}({\nu }_t){|}_{\mathrm{interp}}=a_0^{|\tilde{R}|}+2\sum _{n=1}^N{a}_n^{|\tilde{R}|}.\cos (2\pi \frac{n}{\mathrm{\Delta }{\nu }_t}{\nu }_t),$$

where $\mathrm{\Delta }{\nu }_t={\nu }_{t\_\max }-{\nu }_{t\_\min }$, the size of the determined sample reflectivity spectrum, with boundaries ${\nu }_{t\_\min }=0.1\mathrm{THz}$ and ${\nu }_{t\_\max }=2.0\mathrm{THz}$. $a_0^{|\tilde{R}|}$ and $a_n^{|\tilde{R}|}$ are the interpolation coefficients:

Eq. (8)

$$a_0^{|\tilde{R}|}=\frac{1}{\mathrm{\Delta }{\nu }_t}{\int }_{{\nu }_{t\_\min }}^{{\nu }_{t\_\max }}|\tilde{R}({\nu }_t)|\mathrm{d}{\nu }_t,a_n^{|\tilde{R}|}=\frac{1}{\mathrm{\Delta }{\nu }_t}{\int }_{{\nu }_{t\_\min }}^{{\nu }_{t\_\max }}|\tilde{R}({\nu }_t)|·\cos \left(2\pi \frac{n}{\mathrm{\Delta }{\nu }_t}{\nu }_t\right)\mathrm{d}{\nu }_t.$$

A trigonometric series is also used to interpolate the phase of the complex amplitude reflectivity $\phi \{\tilde{R}({\nu }_t)\}$. $\phi \{\tilde{R}({\nu }_t)\}$ is an odd function, and hence, the interpolating series has the following form:

Eq. (9)

$$\phi {\{\tilde{R}({\nu }_t)\}}_{\mathrm{interp}}=a_0^{\phi \{\tilde{R}\}}{\nu }_t+2\sum _{n=1}^N{a}_n^{\phi \{\tilde{R}\}}·\sin \left(2\pi \frac{n}{\mathrm{\Delta }{\nu }_t}{\nu }_t\right).$$

The interpolation coefficients $a_0^{\phi \{\tilde{R}\}}$ and $a_n^{\phi \{\tilde{R}\}}$ can be written as

Eq. (10)

$$a_0^{\phi \{\tilde{R}\}}=\frac{1}{\mathrm{\Delta }{\nu }_t}{\int }_{{\nu }_{t\_\min }}^{{\nu }_{t\_\max }}{\nu }_t\phi \{\tilde{R}({\nu }_t)\}\mathrm{d}{\nu }_t,a_n^{\phi \{\tilde{R}\}}=\frac{1}{\mathrm{\Delta }{\nu }_t}{\int }_{{\nu }_{t\_\min }}^{{\nu }_{t\_\max }}\phi \{\tilde{R}({\nu }_t)\}·\sin \left(2\pi \frac{n}{\mathrm{\Delta }{\nu }_t}{\nu }_t\right)\mathrm{d}{\nu }_t.$$

The first order moment $a_0^{\phi \{\tilde{R}\}}$ should be used in Eq. (9) instead of the zero order moment, because the sample reflectivity phase $\phi \{\tilde{R}({\nu }_t)\}$ always has a linearly decreasing character.

After reflectivity interpolation, the sample pulse function $R(t)$ can be obtained using the following expressions:

Eq. (11)

$$\begin{gathered} \tilde{R}{({\nu }_t)}_{\mathrm{interp}}=\{\begin{array}{ll}{\begin{array}{c}|\tilde{R}({\nu }_t)|\end{array}}_{\mathrm{interp}}\exp \left({j\phi \{\tilde{R}({\nu }_t)\}}_{\mathrm{interp}}\right),& {\nu }_t\in [-{\nu }_{t\min },{\nu }_{t\min }]\\ \tilde{R}({\nu }_t),& {\nu }_t\notin [-{\nu }_{t\min },{\nu }_{t\min }]\end{array}, \\ R(t)={\mathrm{F}}_{{\nu }_t}^{-1}\left\{\tilde{R}{({\nu }_t)}_{\mathrm{interp}}\right\}. \end{gathered}$$

The interpolation type applied in this work provides good results if one deals with sharp sample permittivity profile variations or with a layered structure. Since the pulse function is known, $R(t)$, sample permittivity profile, $\epsilon (z)$, reconstruction is possible.

2.2.2.

Sample permittivity profile reconstruction

Now consider a plane wave incident normally on the sample surface (Fig. 3). A wave reflected from $\epsilon (z)$ can be described by the following integral equation:¹⁶

Eq. (12)

$$E^r(z_0,t)=\underset{0}{\overset{t}{\int }}{R}^{+}(z_0,t-t^{\prime })·E^i(z_0,t^{\prime })\mathrm{d}{t}^{\prime },$$

where $E^i(z_0,t^{\prime })$ is the incident plane wave amplitude registered at the point $z=z_0$, $E^r(z_0,t)$ is the reflected plane wave amplitude registered at the point $z=z_0$, $R^{+}(z_0,t)$ is the scattering kernel of the integral transformation for signals registered at the point $z=z_0$, and $z$ is the depth coordinate. This linear representation of the reflection process is valid when the reflected signal has no components, arising from the multiple reflections in the object layers. The kernel of the integral transformation $R^{+}(z_0,t)$ describes the region of the investigated medium in the depth range $[z_0,L]$, and depends only on the dielectric permittivity $\epsilon (z)$ and effective conductivity $\sigma (z)$ of the sample.

Fig. 3

Dielectric permittivity profile function.

The first surface of the object is assumed to be the origin of the spatial coordinate OZ. Signals $E^i(z_0,t^{\prime })$, $E^r(z_0,t)$, and the integral kernel $R^{+}(z_0,t)$ can be determined if $z_0<0$, because the field amplitude can be detected only outside the object of interest. The dielectric permittivity is considered to remain constant ${\epsilon }_1$ before the medium $z<0$ and take a different constant value ${\epsilon }_L$ from a certain medium depth $z>L$ (Fig. 3). The following kernel of the integral transformation corresponds to the sample pulse function $R(t)$, obtained at the signal deconvolution step [Eq. (2)]:

Eq. (13)

$$R^{+}(z_0=0,t)=R(t).$$

This kernel characterizes the entire sample. For convenience, it is better to use nondimensional spatial and temporal coordinates and to normalize the integral kernels:¹⁶

Eq. (14)

$$l={\int }_0^L\sqrt{\epsilon (z){\epsilon }_0{\mu }_0}\mathrm{d}z,x=x(z)={\int }_0^Z\frac{\mathrm{d}{z}^{\prime }}{lc(z^{\prime })},s=\frac{t}{l},R(x,s)=l{R}^{+}(z,t),$$

where $l$ is the time needed for the wavefront to travel through the investigated medium (from 0 to $L$), $x$ is the normalized optical depth ($0<x<1$), $s$ is the normalized temporal coordinate ($0<s<2$), and $c(z)$ is the dependence of the speed of light on depth.

All normalized medium kernels $R(x,s)$ should satisfy the following nonlinear integro-differential equation:¹⁶

Eq. (15)

$$\begin{gathered} \frac{\partial R(x,s)}{\partial x}=2\frac{\partial R(x,s)}{\partial s}-B(x)R(x,s)-\frac{1}{2}[A(x)+B(x)]{\int }_0^{s}R(x,s^{\prime })R(x,s-s^{\prime })\mathrm{d}{s}^{\prime },s>0,R(1,s)=0,s>0 \\ R(x,0)=-\frac{1}{4}[A(x)-B(x)], \end{gathered}$$

where $A(x)$ and $B(x)$ are coefficients that depend on medium properties $\epsilon (z)$ and $\sigma (z)$:

Eq. (16)

$$A(x)=-\frac{\mathrm{d}}{\mathrm{d}x}\left[\ln \left(\frac{1}{\sqrt{\epsilon (z(x)){\epsilon }_0{\mu }_0}}\right)\right],B(x)=-\frac{l\sigma (z(x))}{\epsilon (z(x)){\epsilon }_0}.$$

Equations (15) and (16) follow from dissipative wave equation.

As mentioned above, absorption of the studied sample is negligible; that is, $\sigma (z)=0$ $[B(x)=0]$, and therefore, Eq. (15) takes the form:

Eq. (17)

$$\frac{\partial R(x,s)}{\partial x}-2\frac{\partial R(x,s)}{\partial s}=-\frac{1}{2}A(x){\int }_0^{s}R(x,s^{\prime })R(x,s-s^{\prime })\mathrm{d}{s}^{\prime },s>0,R(1,s)=0,s>0,R(x,0)=-\frac{1}{4}A(x)\mathrm{.}$$

To determine sample dielectric permittivity $\epsilon (z)$, the specified initial boundary problem in Eqs. (14), (16), and (17) needs to be solved. The integro-differential Eq. (17) is first solved numerically with the initial condition in Eq. (13) to find the $A(x)$ function. For numerical computation, it is convenient to use the following notation of Eq. (17):

Eq. (18)

$$\frac{\partial R(x,s-2x)}{\partial x}=-\frac{1}{2}A(x){\int }_0^{s-2x}R(x,s^{\prime })R(x,s-2x-s^{\prime })\mathrm{d}{s}^{\prime },s>0,R(1,s)=0,s>0,R(x,0)=-\frac{1}{4}A(x).$$

As $A(x)$ was determined, the dielectric permittivity profile can be found using an inverse transformation of Eq. (16):

Eq. (19)

$$z(x)=c(0)l{\int }_0^x\exp \{-\underset{0}{\overset{x\text{'}}{\int }}A(x^{\prime \prime })\mathrm{d}{x}^{\prime \prime }\}\mathrm{d}{x}^{\prime },0<x<1,\epsilon [z(x)]={\epsilon }_1\exp \{2\underset{0}{\overset{x}{\int }}A(x^{\prime })d{x}^{\prime }\},0<x<1,$$

where $c(0)$ is the speed of light at the point $z=0$. The values $c(0)$ and ${\epsilon }_1$ should be known a priori and are often simply the values for air.

Prior to permittivity profile calculation, the function $A(x)$ can be corrected using some a priori information regarding the sample permittivity. For example, the sample dielectric permittivity has to be constant from a certain depth $z=L$. To fulfill this condition, a constant should be added to $A(x)$ so that the integral $f(x)={\int }_0^{x}A(x^{\prime })\mathrm{d}{x}^{\prime }$ becomes equal to a constant at $x\ge L$

Eq. (20)

$$\begin{gathered} A_{\mathrm{korr}}(x)=A(x)+C, \\ C=[\frac{d}{dx}f(x)]{|}_{x\ge L}. \end{gathered}$$

This correction helps to significantly reduce any error in permittivity profile reconstruction. The permittivity profile reconstruction algorithm is resumed in Fig. 4.

Fig. 4

Scheme of the algorithm.

The algorithm can be modified to take media absorbance into account. At first sight, it is impossible to reconstruct both an arbitrary dielectric permittivity $\epsilon (z)$ and arbitrary conductivity $\sigma (z)$ [or absorption $\alpha (z)$] profile simultaneously, since there is only one algorithm input function $R(t)$. However, one can use some a priori information about the sample conductivity profile $\sigma (z)$ or some approximation of the conductivity distribution during reconstruction. In addition, it is possible to connect the sample conductivity $\sigma (z)$ with sample permittivity $\epsilon (z)$ utilizing a mathematical model of the complex dielectric permittivity $\tilde{\epsilon }$. In this case, reconstruction of some $\tilde{\epsilon }$ parameter profile will take place instead of a direct permittivity $\epsilon (z)$ and conductivity $\sigma (z)$ reconstruction. Note that the function $A(x)$ in Eq. (16) depends on $\epsilon (z)$ and the function $B(x)$ depends on both $\epsilon (z)$ and $\sigma (z)$. We also need to include a complex dielectric permittivity model $\tilde{\epsilon }$ in the description of these functions to deal with absorbing media. After parameter profile reconstruction, one can obtain both profiles of interest from $\tilde{\epsilon }(z)$.

3.

Experimental Results

In order to validate this algorithm, several media with known permittivity profiles were studied. These test media were built from a set of flats with known thickness and THz optical properties. A list of THz materials used in this work and their optical properties are given in Table 1. The optical thicknesses of the test samples were about 5 mm. The THz wave penetration depth in each of the materials listed in Table 1 is much higher than the sample thickness. As a consequence, material absorption can be neglected. The materials selected for the study have small dispersion for the most part of spectral range. Some of them have dispersion at frequencies above 2.0 THz. To neglect it we should reconstruct sample reflectivity $\tilde{R}({\nu }_t)$ only in frequency range under 2.0 THz. Notice, while studying materials with rather high dispersion of optical properties we are reconstructing sample permittivity profile, which corresponds to the average of the sample permittivity in the entire bandwidth of the THz pulse spectrum.

Table 1

Optical properties of THz materials.

The steps for algorithm implementation for the quartz flat in air are as follows. An example of sample $E_s(t)$ and reference $E_r(t)$ waveforms is given in Fig. 5. The reference signal was reflected by a planar gold mirror and the sample amplitude was reflected by a 1.0-mm thick crystalline quartz window. Therefore, two pulses reflected from air/quartz and quartz/air interfaces are present in $E_s(t)$. There are no satellite pulses, which could exist due to etalon effects that arise from multiple reflections within the sample, in signal $E_s(t)$. Appearance of such pulses would lead to distortion of reconstruction results. The noise standard deviation was less than 0.5% relative to the reference peak pulse amplitude. The dynamic range and signal-to-noise ratio (SNR)²² of reference time-domain data were 2000 and 25, respectively.

Fig. 5

Reference $E_r(t)$ and sample (1.0-mm thick quartz window) $E_s(t)$ waveforms.

Figure 6(a) illustrates the absolute value of the sample reflectivity $|\tilde{R}({\nu }_t)|$, obtained after signal deconvolution and interpolation. Function $|\tilde{R}({\nu }_t)|$ has a strongly modulated character due to the presence of two waves in the sample waveform $E_s(t)$. The complexity of reflectivity modulation increases with an increasing number of sample pulses. Such an oscillating function can be represented as a trigonometric series during the low-frequency interpolation procedure in Eqs. (7)–(11). The sample reflectivity is distorted with noise as caused by both detector noise and changes in air composition along the THz beam path during the measurement process. Particularly strong noise lines are due to the fluctuation of water vapor content along the THz beam path.

Fig. 6

Absolute value of reflectivity $|\tilde{R}({\nu }_t)|$ (a) and pulse function $R(t)$ (b) of 1.0-mm thick quartz window.

The quartz window sample pulse function $R(t)$ is presented in Fig. 6(b). The pulse function contains information about the sample interfaces as also for the sample waveform $E_s(t)$. Obviously, the pulse function $R(t)$ contains two types of powerful distortions. First, there is the $R(t)$ distortion caused by a sample reflectivity low-frequency interpolation error. Second, there is a Gibbs noise, which is caused by the rather sharp edges of the Wiener filter. With increasing smoothness of the filter, Gibbs noise decreases but the depth resolution of the permittivity profile reconstruction also correspondingly decreases.

The quartz plate permittivity profile reconstruction is shown in Fig. 7(a). Figure 7 also contains the results for the other test samples: Fig. 7(b) corresponds to the 1.0-mm thick quartz plate and a thick, high-resistivity floating zone silicon plate; Fig. 7(c) corresponds to the sample made of a 1.0-mm thick layer of quartz and a high-density polyethylene layer; and Fig. 7(d) corresponds to the sample made of the same quartz window and a thick polymethylpentene plate. In addition to the algorithm implementation results, any a priori information concerning sample permittivity (ideal profiles) is presented for all four figures.

Fig. 7

Permittivity profiles of test samples: (a) corresponds to 1.0-mm thick crystalline quartz window located in air; (b) corresponds to 1.0-mm thick crystalline quartz window and thick HRFZ-SI plate; (c) corresponds to 1.0-mm thick crystalline quartz window and a thick HDPE plate; and (d) corresponds to 1.0-mm thick crystalline quartz window and a thick TPX plate.

The results for algorithm implementation [Fig. 7(a)–7(d)] and the a priori data for sample permittivity allow us to conclude that profile reconstruction was produced accurately enough for all samples of interest. All reconstructed profiles contain both low- and high-frequency noise, but low-frequency noise caused by interpolation errors leads to a larger distortion of the permittivity profile. Any small Gibbs noise present in the reconstruction results in a much smaller distortion.

4.

Algorithm Stability in the Presence of Noise

The stability of the permittivity profile reconstruction algorithm was considered by mathematical modeling of algorithm implementation in the presence of additive white Gaussian noise in TDS system signals (Fig. 8). Models of layer medium permittivity profiles as well as models of permittivity profiles with slightly smoothed edges were considered. The interaction of THz light with matter was examined and reflected waveforms were obtained. White Gaussian noise with various standard deviations was added to the waveforms and the permittivity reconstruction procedure was applied. Results of algorithm implementation were compared with the initial $\epsilon (z)$ model.

Fig. 8

Flow chart for the study of algorithm stability against noise.

First, a model for the sample permittivity profile was generated and the interaction of THz pulses with the model was considered using a finite-difference time-domain technique for solving Maxwell’s equations (FDTD).²³ The assumptions made during numerical modeling were:

After FDTD modeling, the incident $E^{\text{inc}}(t)$ and reflected $E^{\text{refl}}(t)$ waveforms were obtained for each studied profile $\epsilon (z)$. White Gaussian additive noise with a standard deviation ${\sigma }_{\eta }$ was added to both incident $E^{\text{inc}}(t)$ and reflected $E^{\text{refl}}(t)$ waveforms. Moreover, different noise standard deviations were implemented for each type of permittivity profile model $\epsilon (z)$. A permittivity profile reconstruction procedure was applied to the waveforms with noise added for estimating the permittivity profile ${\epsilon }^{\prime }(z,{\sigma }_{\eta })$. Note, that there was no noise filtering procedure applied before permittivity reconstruction. The estimating permittivity profile ${\epsilon }^{\prime }(z,{\sigma }_{\eta })$ was compared with the initial function $\epsilon (z)$ using the standard deviation of the profiles ${\sigma }_{\epsilon }$:

Figure 9 contains examples of (a) the reference $E^{\text{inc}}(t)$, and (b) the sample $E^{\text{refl}}(t)$ waveforms with different additive noise, as well as (c) results of permittivity profile reconstruction for various waveform noise ${\sigma }_{\eta }$. Figure 9(c) contains both reconstruction results ${\epsilon }^{\prime }(z,{\sigma }_{\eta })$ for different noise levels and the initial permittivity profile function $\epsilon (z)$. Profiles obtained after implementation of the inverse scattering algorithm to the numerical modeling results are close to the initial profile when the additive noise standard deviation ${\sigma }_{\eta }$ is less than 1.9%. Dependence of the standard deviation ${\sigma }_{\epsilon }$ for the profiles on waveform noise standard deviation ${\sigma }_{\eta }$ for the same function $\epsilon (z)$ is shown in Fig. 10. This figure also contains the border between stable (1) and unstable (2) regions.

Fig. 9

Reference (a) and sample (b) waveforms with additive white Gaussian noise, and results of test media permittivity profile reconstruction (c) for various waveform noise.

Fig. 10

Profile standard deviation ${\sigma }_{\epsilon }$ versus waveform noise standard deviation ${\sigma }_{\eta }$.

The dependence of the standard deviation of permittivity profile reconstruction on the noise standard deviation is less than 15.0% for rather high SNR of time-domain data. The quality of $\epsilon (z)$ reconstruction decreases with increasing depth. Therefore, the thinner the investigated medium is, the higher the average accuracy of permittivity profile reconstruction. Some noise-filtering procedure (such as Fourier-domain or wavelet-domain noise filtering) could be implemented to the waveforms to increase the stability border to a higher standard deviation ${\sigma }_{\eta }$ (Fig. 10). Reconstruction results for all examined functions $\epsilon (z)$ were stable at the same waveform noise power for noise standard deviation less than 1.9%. Certainly, the stability and uniqueness of the algorithm should be proven in analytical form, but the results shown in this chapter let us conclude that the algorithm is stable to additive Gaussian white noise.

5.

Discussion

Experimental algorithm implementation and mathematical modeling shows high reconstruction accuracy when layered samples and samples with slightly smoothed borders between layers were studied. The average reconstruction error was smaller than 15.0% for all tested media if the waveform noise standard deviation was in the region of stability (Fig. 10). The border of stability is equal to the additive noise standard deviation ${\sigma }_{\eta }=1.9\%$ and can be raised by applying a noise-filtering procedure. The theoretical limit of depth resolution is considered to be half of the TDS signal minimum wavelength, which exhibits an impulse response $R(t)$. Resolution in the range of 40 to 50 μm can be achieved for the TDS system used in this work.

The accuracy of reconstruction also depends on other factors: interpolation technique, waveform noise filtering method, media dispersion, and media absorption. The accuracy might be increased by implementing modifications to this algorithm, which may include increasing the quality of complex amplitude reflectivity interpolation, developing new methods of reconstructed permittivity profile correction, or applying a noise-filtering procedure. As mentioned above, future algorithm modifications would help to study absorbing samples and samples with dispersion of optical properties described by some complex dielectric permittivity profile model. One of the possible ways to take media absorption into account was discussed above.

6.

Conclusions

In the present work, the ability of THz TDS applied to media dielectric permittivity profile reconstruction was shown. An algorithm for solving the inverse scattering problem based on an invariant embedding technique was developed and implemented. Problems appearing during sample pulse function reconstruction were solved. A pulse function low- and high-frequency filtering procedure based on a Wiener filter was used and trigonometric interpolation of the sample pulse function was applied. The present algorithm can be used in various fields of THz technology, where media with negligible absorption are of interest, such as for nondestructive evaluation of polymer construction materials. The results of algorithm implementation show high reconstruction accuracy for all test samples. Result of algorithm stability modeling aids in the determination of the region of TDS system noise in which reconstruction can be applied successfully. Different methods for modifying the algorithm were discussed. Future work will involve providing the ability to study absorbing media with a dispersion in optical properties.