18 June 2013 Invariant embedding technique for medium permittivity profile reconstruction using terahertz time-domain spectroscopy
Author Affiliations +
Optical Engineering, 52(6), 068203 (2013). doi:10.1117/1.OE.52.6.068203
Abstract
Reconstructing the dielectric permittivity profile (depth dependence of sample dielectric permittivity) is an important inverse problem. We present a new method for permittivity profile reconstruction based on terahertz time-domain spectroscopy signal processing. Reconstruction is accomplished in two steps. First, the sample pulse function is reconstructed using sample time-domain reflection data. Low- and high-frequency noise filtering and the interpolation of the pulse function at low frequencies are then applied. Second, an invariant embedding technique is used to calculate the dielectric permittivity profile based on the sample pulse function. Samples with known permittivity profiles have been studied experimentally using this procedure in order to verify this algorithm. This algorithm is stable to additive Gaussian white noise as shown using mathematical modeling based on the finite-difference time-domain technique. Possible applications of this permittivity profile reconstruction technique are discussed.
Zaytsev, Karasik, Fokina, and Alekhnovich: Invariant embedding technique for medium permittivity profile reconstruction using terahertz time-domain spectroscopy

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.1 For instance, THz imaging and spectroscopy systems are used in the detection of concealed weapons, drugs, and explosives,2,3 medical diagnostics,45.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.9 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,11 study of skin structure,12 study of tooth demineralization,13 and cornea tissue14) and methods of art inspection.15

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 ε(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.18 The THz beam is focused on the sample. Moreover, the aperture ratio of lens is rather small: D/f=12.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.

OE_52_6_068203_f001.png

Two signals need to be obtained for sample permittivity profile reconstruction: Es(t) is the signal reflected from the sample of interest, and Er(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 Es(t) and Er(t) and the reconstruction of the sample permittivity profile ε(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 R˜(νt). Pulse function reconstruction is a signal deconvolution process.

The first problem in signal deconvolution is that simply dividing the sample signal Fourier spectrum E˜s(νt) by the base signal Fourier spectrum E˜r(νt) to reconstruct the sample complex amplitude reflectivity R˜(νt)

(1)

R˜(νt)=E˜s(νt)E˜r(ν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

(2)

R˜(νt)=E˜s(νt)E˜r(νt)F˜(νt).

A technique utilizing a Wiener filter F˜(νt) (see Ref. 19) was applied in previous research20,21 to process TDS system data. The formulation of a derived Wiener filter F˜(νt) with several modifications can be described as follows:

(3)

F˜(νt)=[E˜r(νt)maxνt{E˜r(νt)}]2[E˜r(νt)maxνt{E˜r(νt)}]2+N(νt)S(νt),
where N(νt) is a noise power spectrum model and S(νt) is a signal power spectrum model. Equation (3) is negligible in a frequency range where the noise power spectrum model N(νt) exceeds the signal power spectrum model S(νt). Conversely, if S(νt) exceeds N(νt), the function F˜(ν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(νt):

(4)

N(νt)=K˜,
where K˜ is a constant such that 0.0<K˜<1.0. The signal power spectrum S(νt) is based on a Gaussian monopulse function fGMp(t):

(5)

fGMp(t)=2·e12·(π·νtC·t)·exp[2(π·νtC·t)2],
where νtC is the most powerful harmonic in the monopulse spectrum f˜GMp(νt). The value ν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:

(6)

S(νt)=[Ϝt{fGMp(t)}]2maxνt{[Ϝt{fGMp(t)}]2},
where Ϝt is the direct Fourier transform operator.

Choosing a central monopulse frequency νtC and noise spectral power K˜ in Eqs. (3)–(6), it is possible to provide a qualitative reconstruction of the sample reflectivity for different conditions of waveform Er(t) and Es(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 F˜(νt) generated for a typical reference signal Er(t) with filter parameters: νtC=0.8THz, K˜=0.002. Luckily for our TDS system, the center frequency ν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 E˜r(νt) and Wiener filter F˜(νt)

OE_52_6_068203_f002.png

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 R˜(νt) based on known frequency domain information within the ranges [νt_max,νt_min] and [νt_min,νt_max], and to obtain low-frequency data from the interpolation results. Since the pulse function spectrum R˜(ν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 |R˜(νt)| and phase φ{R˜(ν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 |R˜(νt)|. Since the sample pulse function R(t) is a real function, the series has the form:

(7)

|R˜(νt)|interp=a0|R˜|+2n=1Nan|R˜|.cos(2πnΔνtνt),
where Δνt=νt_maxνt_min, the size of the determined sample reflectivity spectrum, with boundaries νt_min=0.1THz and νt_max=2.0THz. a0|R˜| and an|R˜| are the interpolation coefficients:

(8)

a0|R˜|=1Δνtνt_minνt_max|R˜(νt)|dνt,an|R˜|=1Δνtνt_minνt_max|R˜(νt)|·cos(2πnΔνtνt)dνt.

A trigonometric series is also used to interpolate the phase of the complex amplitude reflectivity φ{R˜(νt)}. φ{R˜(νt)} is an odd function, and hence, the interpolating series has the following form:

(9)

φ{R˜(νt)}interp=a0φ{R˜}νt+2n=1Nanφ{R˜}·sin(2πnΔνtνt).
The interpolation coefficients a0φ{R˜} and anφ{R˜} can be written as

(10)

a0φ{R˜}=1Δνtνt_minνt_maxνtφ{R˜(νt)}dνt,anφ{R˜}=1Δνtνt_minνt_maxφ{R˜(νt)}·sin(2πnΔνtνt)dνt.
The first order moment a0φ{R˜} should be used in Eq. (9) instead of the zero order moment, because the sample reflectivity phase φ{R˜(νt)} always has a linearly decreasing character.

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

(11)

R˜(νt)interp={|R˜(νt)|interpexp(jφ{R˜(νt)}interp),νt[νtmin,νtmin]R˜(νt),νt[νtmin,νtmin],R(t)=Fνt1{R˜(νt)interp}.
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, ε(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 ε(z) can be described by the following integral equation:16

(12)

Er(z0,t)=0tR+(z0,tt)·Ei(z0,t)dt,
where Ei(z0,t) is the incident plane wave amplitude registered at the point z=z0, Er(z0,t) is the reflected plane wave amplitude registered at the point z=z0, R+(z0,t) is the scattering kernel of the integral transformation for signals registered at the point z=z0, 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+(z0,t) describes the region of the investigated medium in the depth range [z0,L], and depends only on the dielectric permittivity ε(z) and effective conductivity σ(z) of the sample.

Fig. 3

Dielectric permittivity profile function.

OE_52_6_068203_f003.png

The first surface of the object is assumed to be the origin of the spatial coordinate OZ. Signals Ei(z0,t), Er(z0,t), and the integral kernel R+(z0,t) can be determined if z0<0, because the field amplitude can be detected only outside the object of interest. The dielectric permittivity is considered to remain constant ε1 before the medium z<0 and take a different constant value ε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)]:

(13)

R+(z0=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:16

(14)

l=0Lε(z)ε0μ0dz,x=x(z)=0Zdzlc(z),s=tl,R(x,s)=lR+(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:16

(15)

R(x,s)x=2R(x,s)sB(x)R(x,s)12[A(x)+B(x)]0sR(x,s)R(x,ss)ds,s>0,R(1,s)=0,s>0R(x,0)=14[A(x)B(x)],
where A(x) and B(x) are coefficients that depend on medium properties ε(z) and σ(z):

(16)

A(x)=ddx[ln(1ε(z(x))ε0μ0)],B(x)=lσ(z(x))ε(z(x))ε0.
Equations (15) and (16) follow from dissipative wave equation.

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

(17)

R(x,s)x2R(x,s)s=12A(x)0sR(x,s)R(x,ss)ds,s>0,R(1,s)=0,s>0,R(x,0)=14A(x).
To determine sample dielectric permittivity ε(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):

(18)

R(x,s2x)x=12A(x)0s2xR(x,s)R(x,s2xs)ds,s>0,R(1,s)=0,s>0,R(x,0)=14A(x).

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

(19)

z(x)=c(0)l0xexp{0x'A(x)dx}dx,0<x<1,ε[z(x)]=ε1exp{20xA(x)dx},0<x<1,
where c(0) is the speed of light at the point z=0. The values c(0) and ε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)=0xA(x)dx becomes equal to a constant at xL

(20)

Akorr(x)=A(x)+C,C=[ddxf(x)]|xL.
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.

OE_52_6_068203_f004.png

The algorithm can be modified to take media absorbance into account. At first sight, it is impossible to reconstruct both an arbitrary dielectric permittivity ε(z) and arbitrary conductivity σ(z) [or absorption α(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 σ(z) or some approximation of the conductivity distribution during reconstruction. In addition, it is possible to connect the sample conductivity σ(z) with sample permittivity ε(z) utilizing a mathematical model of the complex dielectric permittivity ε˜. In this case, reconstruction of some ε˜ parameter profile will take place instead of a direct permittivity ε(z) and conductivity σ(z) reconstruction. Note that the function A(x) in Eq. (16) depends on ε(z) and the function B(x) depends on both ε(z) and σ(z). We also need to include a complex dielectric permittivity model ε˜ in the description of these functions to deal with absorbing media. After parameter profile reconstruction, one can obtain both profiles of interest from ε˜(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 R˜(ν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.

HRFZ-SIQuartzHDPETPX
Permittivity (arb. units)11.679±0.0104.536±0.0502.372±0.0202.135±0.050
Absorption (cm−1)0.025±0.0051.000±0.5000.400±02000.300±0.100
Average penetration depth (mm)400.00010.00025.00033.333

The steps for algorithm implementation for the quartz flat in air are as follows. An example of sample Es(t) and reference Er(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 Es(t). There are no satellite pulses, which could exist due to etalon effects that arise from multiple reflections within the sample, in signal Es(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)22 of reference time-domain data were 2000 and 25, respectively.

Fig. 5

Reference Er(t) and sample (1.0-mm thick quartz window) Es(t) waveforms.

OE_52_6_068203_f005.png

Figure 6(a) illustrates the absolute value of the sample reflectivity |R˜(νt)|, obtained after signal deconvolution and interpolation. Function |R˜(νt)| has a strongly modulated character due to the presence of two waves in the sample waveform Es(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 |R˜(νt)| (a) and pulse function R(t) (b) of 1.0-mm thick quartz window.

OE_52_6_068203_f006.png

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 Es(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.

OE_52_6_068203_f007.png

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 ε(z) model.

Fig. 8

Flow chart for the study of algorithm stability against noise.

OE_52_6_068203_f008.png

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).23 The assumptions made during numerical modeling were:

  • All media of interest were completely nondispersive and nonabsorptive.

  • No nonlinearities of media optical properties were present.

  • The structure and properties of the media did not change in the lateral directions OX and OY, and vary only with depth direction OZ. Therefore, the light interaction with ε(z) (incident and reflected waves) can be described with plane waves running along the OZ axis.

  • The time dependence of the electric field of the incident wave Einc(t) was described with a Gaussian monopulse function [Eq. (5)].

  • There are no satellite pulses caused by multiple reflections in the layers of the object present in the reflected signals.

After FDTD modeling, the incident Einc(t) and reflected Erefl(t) waveforms were obtained for each studied profile ε(z). White Gaussian additive noise with a standard deviation ση was added to both incident Einc(t) and reflected Erefl(t) waveforms. Moreover, different noise standard deviations were implemented for each type of permittivity profile model ε(z). A permittivity profile reconstruction procedure was applied to the waveforms with noise added for estimating the permittivity profile ε(z,ση). Note, that there was no noise filtering procedure applied before permittivity reconstruction. The estimating permittivity profile ε(z,ση) was compared with the initial function ε(z) using the standard deviation of the profiles σε:

(21)

σε(ση)=1Ni=1N[ε(zi,ση)ε(zi)]2,
where N is the number of permittivity profile grid points.

Figure 9 contains examples of (a) the reference Einc(t), and (b) the sample Erefl(t) waveforms with different additive noise, as well as (c) results of permittivity profile reconstruction for various waveform noise ση. Figure 9(c) contains both reconstruction results ε(z,ση) for different noise levels and the initial permittivity profile function ε(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 ση is less than 1.9%. Dependence of the standard deviation σε for the profiles on waveform noise standard deviation ση for the same function ε(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.

OE_52_6_068203_f009.png

Fig. 10

Profile standard deviation σε versus waveform noise standard deviation ση.

OE_52_6_068203_f010.png

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 ε(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 ση (Fig. 10). Reconstruction results for all examined functions ε(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 ση=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.

Acknowledgments

We thank Stanislav O. Yurchenko, Alexei M. Khorokhorov, Vladimir M. Orlov (Bauman Moscow State Technical University, Moscow, Russia) for valuable discussions. This work was partially supported by Grants from the Ministry of Education and Science of Russian Federation (Grant Agreement Nos 14.В37.21.0898 and 14.В37.21.1282) and from the Russian Foundation for Basic Research (Grant Agreement Nos RFBR 12-08-31104 and RFBR 12-08-33112).

References

1. 

Y.-S. Lee, Principles of Terahertz Science and Technology, Springer, New York, NY (2009).Google Scholar

2. 

J. F. Federiciet al., “THz imaging and sensing for security applications: explosives, weapons and drugs,” Semicond. Sci. Technol. 20(7), 266–280 (2005).SSTEET0268-1242http://dx.doi.org/10.1088/0268-1242/20/7/018Google Scholar

3. 

S. R. Murrillet al., “Enhanced terahertz imaging system performance analysis and design tool for concealed weapon identification,” Proc. SPIE 8188, 81880J (2011).PSISDG0277-786Xhttp://dx.doi.org/10.1117/12.898371Google Scholar

4. 

V. P. Wallaceet al., “Terahertz pulsed imaging of cancers,” Proc. SPIE 4949, 352–359 (2003).PSISDG0277-786Xhttp://dx.doi.org/10.1117/12.500121Google Scholar

5. 

P. C. Ashworthet al., “Terahertz pulsed spectroscopy of freshly excised human breast cancer,” Opt. Express 17(15), 12444–12454 (2009).OPEXFF1094-4087http://dx.doi.org/10.1364/OE.17.012444Google Scholar

6. 

C. Reid, “Spectroscopic methods for medical diagnosis at terahertz wavelength,” Ph.D. Thesis, University College London (2009).Google Scholar

7. 

C. D. Stoik, “Nondestructive evaluation of aircraft composites using terahertz time domain spectroscopy,” Ph.D. Thesis, Air Force Institute of Technology (2008).Google Scholar

8. 

N. Karpowiczet al., “Fire damage on carbon fiber materials characterization by THz waves,” Proc. SPIE 6212, 62120G (2006).PSISDG0277-786Xhttp://dx.doi.org/10.1117/12.665852Google Scholar

9. 

D. M. Mittlemanet al., “T-ray tomography,” Opt. Lett. 22(12), 904–906 (1997).OPLEDP0146-9592http://dx.doi.org/10.1364/OL.22.000904Google Scholar

10. 

W. L. ChanJ. DeibelD. M. Mittleman, “Imaging with terahertz radiation,” Rep. Prog. Phys. 70(8), 1325–1379 (2007).RPPHAG0034-4885http://dx.doi.org/10.1088/0034-4885/70/8/R02Google Scholar

11. 

M. H. Arbabet al., “Terahertz reflectometry of burn wounds in a rat model,” Biomed. Opt. Express 2(8), 2339–2347 (2011).BOEICL2156-7085http://dx.doi.org/10.1364/BOE.2.002339Google Scholar

12. 

E. PickwellV. P. Wallace, “Biomedical applications of terahertz technology,” J. Phys. D: Appl. Phys. 39(17), R301–R310 (2006).JPAPBE0022-3727http://dx.doi.org/10.1088/0022-3727/39/17/R01Google Scholar

13. 

E. Pickwellet al., “A comparison of terahertz pulsed imaging with transmission microradiography for depth measurement of enamel demineralisation in vitro,” Caries Res. 41(1), 49–55 (2007).CAREBK0008-6568http://dx.doi.org/10.1159/000096105Google Scholar

14. 

D. B. Bennettet al., “Terahertz sensing in corneal tissues,” J. Biomed. Opt. 16(5), 057003 (2011).JBOPFO1083-3668http://dx.doi.org/10.1117/1.3575168Google Scholar

15. 

J.-P. Caumeset al., “Terahertz tomographic imaging of XVIIIth Dynasty Egyptian sealed pottery,” Appl. Opt. 50(20), 3604–3608 (2011).APOPAI0003-6935http://dx.doi.org/10.1364/AO.50.003604Google Scholar

16. 

G. KristenssonR. J. Krueger, “Direct and inverse scattering in the time domain for a dissipative wave equation: I. Scattering operators,” J. Math. Phys. 27(6), 1667–1682 (1986).JMAPAQ0022-2488http://dx.doi.org/10.1063/1.527083Google Scholar

17. 

G. KristenssonR. J. Krueger, “Direct and inverse scattering in the time domain for a dissipative wave equation: II. Simultaneous reconstruction of dissipation and phase velocity profiles,” J. Math. Phys. 27(6), 1683–1693 (1986).JMAPAQ0022-2488http://dx.doi.org/10.1063/1.527084Google Scholar

18. 

B. Schulkin, “Miniature terahertz time-domain spectrometry,” Ph.D. Thesis, Rensselaer Polytechnic Institute (2008).Google Scholar

19. 

R. C. GonzalezR. E. Woods, Digital Image Processing, Prentice Hall, Upper Saddle River, NJ (2007).Google Scholar

20. 

D. M. Mittlemanet al., “Gas sensing using terahertz time-domain spectroscopy,” Appl. Phys. B 67(3), 379–390 (1998).APBOEM0946-2171http://dx.doi.org/10.1007/s003400050520Google Scholar

21. 

Y. ChenS. HuangE. Pickwell-MacPherson, “Frequency-wavelet domain deconvolution for terahertz reflection imaging and spectroscopy,” Opt. Express 18(2), 1177–1190 (2010).OPEXFF1094-4087http://dx.doi.org/10.1364/OE.18.001177Google Scholar

22. 

M. NaftalyR. Dudley, “Methodology for determining the dynamic ranges and signal-to-noise ratios of terahertz time-domain spectrometers,” Opt. Lett. 34(8), 1213–1215 (2009).OPLEDP0146-9592http://dx.doi.org/10.1364/OL.34.001213Google Scholar

23. 

A. TafloveS. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, Arteck House, Norwood, MA (2000).Google Scholar

Biography

Kirill I. Zaytsev is currently a PhD student at Bauman Moscow State Technical University (BMSTU), Moscow, Russia. His research interests include THz technology, THz spectroscopy, inverse problems in optics, digital signal processing, and computational electrodynamics. Since 2011 he has worked at the Research and Educational Center “Photonics and Infrared Technology” as an associate researcher. He is a member of SPIE.

Valeriy E. Karasik completed his PhD at Bauman Moscow State Technical University (BMSTU), Moscow, USSR. He received DrSc degree in 2002 (at BMSTU). Since 2002 Valeriy Karasik became a professor of BMSTU. He is the author of over 100 publications and conference presentations. He is a fellow of SPIE and OSA, as well as head of Research and Educational Center “Photonics and Infrared Technology.”

Irina N. Fokina is currently a PhD student at Bauman Moscow State Technical University (BMSTU), Moscow, Russia. Her research interests include THz technology, THz imaging systems, and problems of electromagnetic waves scattering. Since 2011 she has worked at the Research and Educational Center “Photonics and Infrared Technology” as an associate researcher. She is a member of SPIE and OSA.

Valentin I. Alekhnovich received his PhD degree at Bauman Moscow State Technical University (BMSTU), Moscow, Russia. He is currently an assistant professor of applied mathematics chair (Fundamental Sciences Department) and an assistant professor of laser and optoelectronic systems chair (Department of Radio Electronics and Laser Technology) at BMSTU. He is the author of over 100 publications and conference presentations. His research interests include laser spectroscopy, inverse and ill-posed problems in optics, optical, and laser devices.

Kirill I. Zaytsev, Valery E. Karasik, Irina N. Fokina, Valentin I. Alekhnovich, "Invariant embedding technique for medium permittivity profile reconstruction using terahertz time-domain spectroscopy," Optical Engineering 52(6), 068203 (18 June 2013). http://dx.doi.org/10.1117/1.OE.52.6.068203
Submission: Received ; Accepted
JOURNAL ARTICLE
11 PAGES


SHARE
KEYWORDS
Terahertz radiation

Reconstruction algorithms

Dielectrics

Reflectivity

Terahertz spectroscopy

Quartz

Spectroscopy

Back to Top