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 is reconstructed based on time-domain signals. Second, the sample permittivity profile is reconstructed utilizing the pulse function . 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.
Materials and Methods
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: . 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.
Two signals need to be obtained for sample permittivity profile reconstruction: is the signal reflected from the sample of interest, and 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.
Dielectric Permittivity Profile Reconstruction Algorithm
As previously noted, sample dielectric permittivity profile reconstruction is accomplished in two steps: reconstruction of the sample pulse function based on TDS system signals and and the reconstruction of the sample permittivity profile based on a sample pulse function .
Sample pulse function reconstruction
The sample pulse function 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 . Pulse function reconstruction is a signal deconvolution process.
The first problem in signal deconvolution is that simply dividing the sample signal Fourier spectrum by the base signal Fourier spectrum to reconstruct the sample complex amplitude reflectivity1) 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
A technique utilizing a Wiener filter (see Ref. 19) was applied in previous research20,21 to process TDS system data. The formulation of a derived Wiener filter with several modifications can be described as follows:3) is negligible in a frequency range where the noise power spectrum model exceeds the signal power spectrum model . Conversely, if exceeds , the function is close to unity. Smooth filter edges are also provided. For simplicity, we assume a Gaussian white noise model for the noise power spectrum :
Choosing a central monopulse frequency and noise spectral power in Eqs. (3)–(6), it is possible to provide a qualitative reconstruction of the sample reflectivity for different conditions of waveform and 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 generated for a typical reference signal with filter parameters: , . Luckily for our TDS system, the center frequency 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.
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 based on known frequency domain information within the ranges and , and to obtain low-frequency data from the interpolation results. Since the pulse function spectrum 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 and phase of the sample complex reflectivity can be interpolated separately. It is convenient to use a trigonometric interpolation series to reconstruct the amplitude reflectivity modulus . Since the sample pulse function is a real function, the series has the form:
A trigonometric series is also used to interpolate the phase of the complex amplitude reflectivity . is an odd function, and hence, the interpolating series has the following form:9) instead of the zero order moment, because the sample reflectivity phase always has a linearly decreasing character.
After reflectivity interpolation, the sample pulse function can be obtained using the following expressions:
Sample permittivity profile reconstruction
The first surface of the object is assumed to be the origin of the spatial coordinate OZ. Signals , , and the integral kernel can be determined if , because the field amplitude can be detected only outside the object of interest. The dielectric permittivity is considered to remain constant before the medium and take a different constant value from a certain medium depth (Fig. 3). The following kernel of the integral transformation corresponds to the sample pulse function , obtained at the signal deconvolution step [Eq. (2)]:16
All normalized medium kernels should satisfy the following nonlinear integro-differential equation:1615) and (16) follow from dissipative wave equation.
As mentioned above, absorption of the studied sample is negligible; that is, , and therefore, Eq. (15) takes the form: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 function. For numerical computation, it is convenient to use the following notation of Eq. (17):
As was determined, the dielectric permittivity profile can be found using an inverse transformation of Eq. (16):
Prior to permittivity profile calculation, the function 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 . To fulfill this condition, a constant should be added to so that the integral becomes equal to a constant atFig. 4.
The algorithm can be modified to take media absorbance into account. At first sight, it is impossible to reconstruct both an arbitrary dielectric permittivity and arbitrary conductivity [or absorption ] profile simultaneously, since there is only one algorithm input function . However, one can use some a priori information about the sample conductivity profile or some approximation of the conductivity distribution during reconstruction. In addition, it is possible to connect the sample conductivity with sample permittivity 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 and conductivity reconstruction. Note that the function in Eq. (16) depends on and the function depends on both and . 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 .
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 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.
Optical properties of THz materials.
|Permittivity (arb. units)||11.679±0.010||4.536±0.050||2.372±0.020||2.135±0.050|
|Average penetration depth (mm)||400.000||10.000||25.000||33.333|
The steps for algorithm implementation for the quartz flat in air are as follows. An example of sample and reference 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 . There are no satellite pulses, which could exist due to etalon effects that arise from multiple reflections within the sample, in signal . 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.
Figure 6(a) illustrates the absolute value of the sample reflectivity , obtained after signal deconvolution and interpolation. Function has a strongly modulated character due to the presence of two waves in the sample waveform . 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.
The quartz window sample pulse function is presented in Fig. 6(b). The pulse function contains information about the sample interfaces as also for the sample waveform . Obviously, the pulse function contains two types of powerful distortions. First, there is the 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.
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.
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 model.
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 (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 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 and reflected waveforms were obtained for each studied profile . White Gaussian additive noise with a standard deviation was added to both incident and reflected waveforms. Moreover, different noise standard deviations were implemented for each type of permittivity profile model . A permittivity profile reconstruction procedure was applied to the waveforms with noise added for estimating the permittivity profile . Note, that there was no noise filtering procedure applied before permittivity reconstruction. The estimating permittivity profile was compared with the initial function using the standard deviation of the profiles :
Figure 9 contains examples of (a) the reference , and (b) the sample 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 for different noise levels and the initial permittivity profile function . 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 is shown in Fig. 10. This figure also contains the border between stable (1) and unstable (2) regions.
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 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 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.
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 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 . 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.
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.
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).
Y.-S. Lee, Principles of Terahertz Science and Technology, Springer, New York, NY (2009).Google Scholar
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
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
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
C. Reid, “Spectroscopic methods for medical diagnosis at terahertz wavelength,” Ph.D. Thesis, University College London (2009).Google Scholar
C. D. Stoik, “Nondestructive evaluation of aircraft composites using terahertz time domain spectroscopy,” Ph.D. Thesis, Air Force Institute of Technology (2008).Google Scholar
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
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
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
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
B. Schulkin, “Miniature terahertz time-domain spectrometry,” Ph.D. Thesis, Rensselaer Polytechnic Institute (2008).Google Scholar
R. C. GonzalezR. E. Woods, Digital Image Processing, Prentice Hall, Upper Saddle River, NJ (2007).Google Scholar
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
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
A. TafloveS. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, Arteck House, Norwood, MA (2000).Google Scholar
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.