1. Introduction

There are several ways of approaching THz imaging including raster scanning and tomographic reconstruction [1]. However, such techniques are limited by data acquisition speed. A more promising imaging technique is time-reversal imaging which has successfully demonstrated fast reconstruction of one-, two-, and even three-dimensional objects [2,3]. Time-reversal imaging is a coherent technique in which the time reversal symmetry of Maxwell’s wave equations is exploited to reconstruct an object from the scattered fields. One measures the scattered fields from an object at multiple detector positions and the object is reconstructed by numerically back-propagating the scattered fields. More specifically, an image point is reconstructed by first computing the time delays for a THz pulse to propagate from that point to every detector position. The received signal samples are then summed together in order to produce the amplitude of the image point. Each image point requires a different set of time delays and hence a different summation of the scattered signals. This method of time-reversing the scattered data is based on the “delay-sum” algorithm, commonly used in ultrasound imaging, and is nearly equivalent to back propagating the scattered fields using the time-reversed Huygens-Fresnel diffraction integral but requires less computation [3].

As is well-known, real-time imaging in ultrasound using delay-sum or similar algorithms is enabled through the use of multiple parallel emitters and detectors. The main objective of this work is to improve on the THz time-reversal technique in two ways. We first adapt a waveguide approach previously pioneered in ultrasound to increase the effective numerical aperture of the system without decreasing the data acquisition speed of the THz system [4]. Secondly, we implemented a model-based reconstruction technique that uses the actual impulse response of the THz system and is therefore better suited for reconstructing the object and eliminating spurious signals than the a priori time-reversal algorithm [5]. We demonstrate these methods for improved time-reversal imaging using a single emitter and a scanning single detector (i.e. in a form of synthetic-aperture imaging); however, the extension of these methods to multiple emitter and detector implementations is expected to be straightforward, and indicates a promising route to real-time THz imaging.

2. Experimental System

The experimental setup, as given in Fig. 1, is a typical electro-optic THz sampling system [6] with an additional stage in the pump arm to compensate for the horizontal translation of one of the imaging parabolas. A femtosecond laser pulse is split into a pump and probe pulse by a beam splitter. The delayed pump pulse illuminates a large-area photoconductive emitter (TeraSED, GigaOptics GmbH) to generate a nearly single cycle THz pulse. The THz beam is then collimated by a polyethylene lens with a focal length of 7.6 cm. The collimated THz beam is used to illuminate an object, which in this experiment consists of two metal slits with dimensions of 1 mm×8 mm and a spacing of 2.0 mm. The slits are bounded by two 12 inch flat mirrors which act as a planar waveguide redirecting the THz scattered at large angles and therefore delivering higher spatial frequencies into the detection region of the imaging parabolas. The focal plane of parabola A is imaged onto the electro-optic (EO) crystal by scanning parabola A across the exit face of the waveguide [7]. A pellicle reflects the probe pulse to propagate collinearly with the THz pulse in the EO crystal at the focus of parabola B. The EO crystal is a (1 1 0) cut ZnTe crystal which velocity matches the THz pulse and near IR probe pulse to enable coherent detection of the THz pulse. The THz pulse induces a birefringence in the crystal through the linear electro-optic effect (Pockels effect), which is probed by the linearly polarized sampling pulse. The induced phase modulation of the probe pulse is converted into an intensity modulation and detected by a differential photodiode [6].

 

Fig. 1. Experimental setup

Download Full Size | PDF

3. Waveguide-enhanced time-reversal imaging

The resolution of an imaging system is limited by its numerical aperture (NA). A higher NA can be obtained by collecting data at larger angles which requires scanning the detector over more spatial positions and hence leads to a longer acquisition time. However the waveguide technique can effectively increase the numerical aperture of the imaging setup without increasing the number of spatial scan positions. The waveguide technique described here was adapted from an ultrasound experiment conducted in a water channel bounded by two plane interfaces [4]. The ultrasound experiment showed an improvement in both the spatial and temporal compression of the time-reversed fields by simply time reversing the direct path signal and the set of multipath signals corresponding to multiple reflections of the incident wave on the interfaces. Similarly, we can intentionally introduce multipath into our THz imaging setup by bounding our object with two planar mirrors and time reversing both the direct signal and the multipath signals. We can then invoke the principle of mirror images to explain how bounding our object with the planar mirrors can effectively increase the numerical aperture of our imaging system. Each reflected pulse that is detected corresponds to a virtual detector position [4]. Hence, we can effectively double our numerical aperture with virtual detector positions by simply capturing the first set of reflections off the mirrors and accounting for their proper time delays in the time-reversal algorithm [2,4]. Furthermore, the reflected pulses which diffracted at larger angles than the direct path signals have a higher spatial frequency content and thus by using a waveguide to redirect them we can improve the resolution of our system. Thus, sampling more spatial points translates to simply scanning longer in time.

 

Fig. 2. Measured THz wavefield plots without a waveguide and with a symmetric and asymmetric placed waveguide and the corresponding reconstructed images based on time-reversal.

Download Full Size | PDF

In our experimental demonstration, we illuminated a double slit and measured the scattered THz radiation in the far field. The THz wavefield plots shown in Fig. 2(a) & 2(c) were obtained by scanning parabola A horizontally in increments of one millimeter over a range of 52 mm and a range of 59 mm respectively. At each position, the time domain THz waveform was measured by scanning the delay between pump and probe pulses over a 40 s acquisition time. Hence, the y-axis represents time delay and the x-axis represents effective detector position at the exit of the waveguide. We then carried out the same experiment with the waveguide symmetrically and asymmetrically placed about the object as shown in Fig. 2(e) and Fig. 2(g). The THz wavefield plots in the waveguide cases show, in addition to the direct path signals, pulses arriving at a later time corresponding to a single reflection from the waveguide mirrors. By the principle of mirror images, each reflected pulse that is detected corresponds to a virtual detector position. Hence, the first set of reflected signals should effectively double the numerical aperture of the THz system. In practice, as the angle of an image point increases corresponding to the arrival of reflected pulses later in time, the signal strength decreases due to the limited acceptance angle of the imaging parabola. Thus only the first set of reflections is captured in our system and it is not possible to achieve a full doubling of the numerical aperture of the system. In ultrasound, the acoustic detectors are isotropic, enabling more reflections to be captured; the use of a shorter focal length parabola and a waveguide with a smaller aspect ratio should yield more reflections and thus enable a larger numerical aperture for the THz imaging system.

The reconstructed images in Fig. 2(b), 2(d), 2(f) and 2(h) were obtained by numerically back-propagating the corresponding experimental wavefield plots using the delay-sum algorithm. By accounting for the proper time delays of the reflected signals, we were able to achieve better spatio-temporal compression of the time-reversed fields with the waveguide in place regardless of its symmetry about the object. We can quantify the temporal and spatial improvement as a result of using a waveguide by taking a horizontal slice through the reconstructed images. Figure 3(a) & 3(c) show an intensity enhancement for the symmetric and asymmetric waveguides of 2.6 and 1.9 respectively. These values were computed by taking the ratio of the peaks of the waveguide curve (red) to the non-waveguide curve (blue). The increase in intensity for the waveguide cases is attributed to an enhancement of the fields. This enhancement is the result of the coherent addition upon back-propagation of the pulses reflected from the waveguide walls to the direct pulses, resulting in a larger constructive interference of their maxima as well as a larger destructive interference of side lobes [4]. This is in contrast to the case without the waveguide, in which only the direct pulses are available for back-propagation and coherent addition. Figure 3(b) & 3(d) also shows that the waveguide in both the symmetric and asymmetric cases has led to better spatial focusing of the time-reversed fields [4]. From the resolution plot for the symmetric waveguide case, Fig. 3(b), we computed a full width half maximum (FWHM) value of 1.08 mm for the first peak of the red curve and a FWHM of 1.64 mm for the respective peak of the blue curve. The slits that we imaged, as mentioned earlier, have widths of 1 mm.

 

Fig. 3. Intensity plots for symmetric (a) and asymmetric (c) waveguide configurations. Resolution plots for symmetric (b) and asymmetric (d) waveguide configurations.

Download Full Size | PDF

Likewise, for the asymmetric waveguide case, Fig. 3(d), we computed a FWHM of 1.16 mm for the first peak of the red curve and a FWHM of 1.36 mm for the respective blue curve. In both waveguide cases, the time-reversal reconstruction of the two slits yielded a reconstructed object with dimensions much closer to the true dimensions of the two slits. However, in the asymmetric case the blue curve of Fig. 3(d) has a smaller FWHM than the blue curve in the symmetric case of Fig. 3(b) and this is consistent with the fact that in the asymmetric case we scanned more positions and hence we should have achieved better reconstruction. However, the FWHM of the red curve in the asymmetric case, Fig. 3(d), did not match the FWHM of the red curve of the symmetric case, Fig. 3(b), even though we sampled more detector positions. The reason for this discrepancy is evident in the wave field plot for the asymmetric case, Fig. 2(g). One side of the waveguide was much closer to the object than the other side. Hence we have more reflected pulses coming from the farther side than the closer side which will significantly contribute to a sharper rise on that side of the reconstructed object. Furthermore, the few reflected pulses from the closer side that are present in the wave field plot contribute very little to the reconstruction of the object because their amplitudes are small because they correspond to an image point with a large angle. Hence the right part of the red curve in the asymmetric case is overlapped with the corresponding part of the blue curve. As mentioned previously, achieving a larger angular acceptance in the THz detection will enable a higher effective NA and would alleviate this problem. In our experimental demonstration, we were still far from the ultimate diffraction limit of our system which has a λ_peak=429 µm and a λ_mean=119 µm [2].

4. Model-based image reconstruction

Ideally the time-reversal algorithm enables the realization of an optimal spatio-temporal filter as a result of the reciprocity theorem, which states that the position of a source and receiver can be interchanged without altering the resulting field [4]. We have shown that the introduction of a waveguide has effectively increased the NA of our setup and thus we have achieved a better spatio-temporal compression of our time-reversed fields than without the waveguide. However, the performance of the time-reversal algorithm is nonetheless degraded by the presence of temporal ringing on the THz pulse. Ringing arises in the system due to reflections in the THz emitter and detector, atmospheric absorption lines, and the non-ideal response of the ZnTe electro-optic crystal. Although our THz pulse is far from being a clean single cycle pulse, we can mitigate the effect of the ringing in our reconstruction algorithm by taking into account the measured impulse response of our THz imaging system. That is, we can approach image reconstruction from scattered fields as a model-based inverse problem in which we try to recover some underlying function that describes the object from the collected data in a “best fit” manner without overly fitting the noise [5]. Hence we have investigated replacing the time-reversal reconstruction algorithm with a more general statistical algorithm that estimates what the object is from the data collected [5].

At every detector position, the received signal is just a superposition of THz pulses from every point in the object plane with an appropriate delay parameter. That is the observed signal at the mth detector position can be expressed as

where θ_n denotes the unknown value of the object’s transmissivity at the nth sample position in the object plane, and h(t) is a THz pulse that is delayed by a known parameter τnm. By concatenating our observed signals into one vector, we can recast the above equation as:

or more succinctly as:

where Y is a vector consisting of observed signals, A is a known system matrix and θ is a vector of unknown parameters. We could find an estimate for θ from Y by minimizing the following least-squares criterion:

However, since our goal is not only to obtain an estimate of θ but also to reduce the presence of artifacts in our reconstructed images, we minimize instead the following regularized least-squares (RLS) cost function:

The additional term is a regularizing penalty term and its effect is to discourage disparities in neighboring pixel values while the effect of the first term is to encourage a best fit of the measured data [5]. Since these two effects are conflicting the adjustable parameter β controls the tradeoff between the two and controls the balance between spatial resolution and noise in the final estimate [5].

We implemented a one-dimensional reconstruction algorithm based on the RLS criterion given in Eq. (5) and compared its performance to the performance of the time-reversal algorithm in the case without a waveguide and in the case with an asymmetrically placed waveguide. The algorithm for both experiments took on average 5 iterations to converge. In the case without the waveguide, the RLS algorithm achieves a better reconstruction of the object than the time-reversal algorithm as evident by the 4.5 x improvement in intensity as shown in Fig. 4(a) and we have also calculated a peak SNR ratio improvement of 2.2. The RLS algorithm has also improved the resolution of the system. The FWHM for the red curve is 1.03 mm and the FWHM for blue curve is 1.36 mm. Furthermore, the accuracy of our system model can be determined by how well the simulated wavefield data, which can be obtained by multiplying the reconstructed object from the RLS algorithm by the system matrix A, matches the measured wavefield data. The measured wavefield data of Fig. 4(c) and the simulated wavefield data in Fig. 4(d) are well matched, and the presence of spurious signals due to imperfections in the imaging optics have been removed from the measured data. Furthermore, we can better show the accuracy of our system model by taking any arbitrary vertical time slice from wavefield plots of Fig. 4(c) & 4(d) and determining how well the measured THz and the simulated THz pulses match at a particular detector position. Figure 4(e) shows a comparison between the measured and the forward-propagated RLS-reconstructed THz pulse at detector position 10. The simulated pulse shows strong agreement with the measured direct path THz pulse at -25 ps, while showing strong suppression of spurious signals, most notably the one present at approximately -15 ps.

 

Fig. 4. Absolute (a) and normalized intensity (b) obtained by the RLS algorithm (red curve) versus the time-reversal algorithm (blue curve). Figure 4(c) & (d) show the similarity of the measured and simulated fields. Figure 4(e) shows the measured and forward-propagated RLS-reconstructed THz waveform.

Download Full Size | PDF

The model-based reconstruction algorithm works well in reconstructing the scattering object from the measured data. However in the waveguide case, the algorithm did not yield a substantial improvement over the time-reversal algorithm as seen in the intensity and the resolution plots of Fig. 5(a) & 5(b) respectively. We computed a peak SNR improvement of 1.47 for model-based algorithm over the time-reversal algorithm for the case with the waveguide. The lack of a substantial improvement in the waveguide case can be attributed to a combination of two factors. The first factor is the sensitivity of the RLS approach to modeling errors in the system matrix A. If there is a phase shift that we are not accounting for in the waveguide case, the impulse function used to create the A matrix will not be able to completely model the waveguide system, particularly the reflections off the mirrors. We can see in Fig. 5(c) & 5(d) that there is a discrepancy between the measured and simulated wavefield data indicating the presence of an unknown phase shift. Figure 5(e) further shows that although we modeled the direct part of the signal accurately as evident by the strong agreement between the direct parts of the measured and simulated signals, there is a mismatch between the reflected parts. The other factor preventing a substantial improvement in reconstruction for the waveguide case is the diffraction limit of the system. The presence of the waveguide has improved the NA of the system beyond the acceptance angle of the imaging parabolas and hence the RLS algorithm has very little to improve upon as it asymptotically approaches the diffraction limit of the system. Nonetheless, the model-based algorithm in conjugation with the waveguide performs better than the time-reversal algorithm with the waveguide and we anticipate better performance provided we can generate the correct system matrix.

 

Fig. 5. Absolute (a) and normalized intensity (b) obtained by the RLS algorithm (red curve) versus time-reversal with waveguide (blue curve). Both plots show the peak intensity relative to the time-reversal algorithm without the waveguide. Figure 5(c) & (d) show the slight discrepancy between the measured and simulated wavefield data. Figure 5(e) shows the measured and simulated THz signals at detector position 10.

Download Full Size | PDF

5. Conclusion

In conclusion we have presented two methods for improving the time-reversal imaging technique. We first used a waveguide to increase the effective numerical aperture of the system. The waveguide technique not only yields an improvement in the numerical aperture of the system, but more generally illustrates how techniques used in ultrasound may be fruitfully adapted to THz imaging technology. Secondly, we implemented a model-based reconstruction technique that uses the actual impulse response of the experimental THz system and is therefore better suited for reconstructing the object and eliminating spurious signals than the simple time-reversal algorithm. We have demonstrated the model-based algorithm for a THz system operating in transmission mode; however, we can easily extend this algorithm for systems operating in reflection mode provided that the impulse response for the system is known in advance in order to construct the system matrix.