Comparison of digital beamforming algorithms for 3-D terahertz imaging with sparse multistatic line arrays

In this contribution we compare the backprojection algorithm with our recently developed modified range migration algorithm for 3-D terahertz imaging using sparse multistatic line arrays. A 2-D planar sampling scheme is generated using the array’s aperture in combination with an orthogonal synthetic aperture obtained through linear movement of the object under test. A stepped frequency continuous wave signal modulation is used for range focusing. Comparisons of the focusing quality show that results using the modified range migration algorithm reflect these of the back-projection algorithm except for some degradation along the array’s axis due to the operation in the array’s near-field. Nevertheless the highest computational efficiency is obtained from the modified range migration algorithm, which is better than the numerically optimized version of the back-projection algorithm. Measurements have been performed by using an imaging system operating in the W frequency band to verify the theoretical results.


Introduction
The use of the effective aperture of a sparse multistatic array reduces the number of required transmitters (T x ) and receivers (R x ) and preservers the array's imaging quality (Lockwood et al., 1996).Hence, this approach is widely spread in the fields of radar and ultrasonic imaging (Wiesbeck and Sit, 2014;Thomenius, 1996).The principle idea of the effective aperture concept is to illuminate the target object with an array of transmitters (T x -array) and record backscattered radiation using an array of receivers (R x -array).
The T x -array and the R x -array are designed in a way that the convolution of their apertures results in a dense effective aperture.In a subsequent step digital beam forming (DBF) algorithms are applied to the recorded data for 3-D image reconstruction of the target object.Furthermore, the use of DBF techniques does not pose any constraints on the kind of aperture sampling.So one can replace a physical array by an equivalent synthetic one, or combine both or even generate a sampling aperture that corresponds to the object shape.This approach has been adopted in different system designs.In Zhuge and Yarovoy (2011) for instance, the authors reported on the development of an imaging system at a center frequency of ca.11 GHz, which generates a 2-D sampling aperture through combining a sparse multistatic line array with a synthetic aperture generated by the movement of the array.However, in Ahmed et al. (2011) the authors used a 2-D planar array for 3-D imaging at a center frequency of 76 GHz.In the field of terahertz imaging there is also a growing interest for using sparse arrays, since imaging using quasi-optics in combination with scanning stages doesn't fulfill real-time operation requirements and focal-plane arrays (FPAs) pose a trade-off between the system's field of view and the resolution (Friederich et al., 2011).Hence, a major concern about such systems is the acceleration of data acquisition processes as well as faster image generation.
In order to achieve real-time operation, computational efficient DBF algorithms are mandatory.The range migration algorithm (RMA) is a very computational efficient algorithm since it is mainly implemented using the fast Fourier transform (FFT) algorithm.It is widely used for synthetic aperture radar (SAR) imaging and is established for monostatic and bistatic configurations (Soumekh, 1991;Lopez-Sanchez and Fortuny-Guasch, 2000).Under the assumption of a fully populated multistatic array an implementation of the RMA is presented in Zhuge and Yarovoy (2011).However, the use of sparse arrays hinders a direct implementation of the RMA because of the violation of the Nyquist sampling criterion along the sampling aperture.
In this contribution we compare a modified RMA with the back-projection (BP) algorithm and the fast-factorized backprojection (FFBP) algorithm by taking the focusing quality and the asymptotic computing complexity as criteria.
This paper is organized as follows.In Sect. 2 we formulate the imaging problem.The imaging algorithms are briefly described in Sect.3. Required sampling criteria and resulting resolution are discussed in Sect. 4. The algorithms are compared in Sect. 5.In Sect.6 terahertz image reconstructions using the mentioned algorithms from measurement data are presented.Conclusions are drawn in Sect.7.

Formulation of the imaging problem
A schematic of the exemplary imaging scenario is depicted in Fig. 1.Along the y axis a line array of N T transmitters illuminates the measurement scene.The transmitters are operated sequentially.A line array of N R receivers records backscattered radiation from the scene.The T x -array and R xarray are placed on the same distance Z o from the center of the volume to be imaged (for the sake of clarity they are illustrated separately).The coordinates of the transmitters and receivers are then given by r T = (0, u T , Z o ) and r R = (0, u R , Z o ), respectively.The elements of each array are uniformly distributed along the y axis.Both arrays are combined to generate an equivalent effective aperture.Under far-field conditions the resulting effective aperture is the convolution of the apertures of the T x -array and the R x -array (Lockwood et al., 1996).Hence the resulting number of effective sampling ele-ments is N E = N T ×N R .The main idea of this approach is to design the aperture of one array, here the T x -array, in a sparse fashion of the required effective aperture and to design the other array, here respectively the R x -array, in a denser fashion so that its aperture is used as an interpolation function.Thus, the element spacing in the effective aperture is equal to the R x -array element spacing d R .To generate a uniform effective aperture, the element spacing in the T x -array should be d T = N R × d R .The coordinates of an effective aperture element are given by, (1) The overall extent of the effective aperture is given by L E = (N E −1)×d R .The target object is translated along the x axis, thus mechanically creating a synthetic aperture orthogonally to the array's aperture.The translation along the x axis is described by the vector r v = (v, 0, 0).The signal measured at the effective aperture position u and the synthetic aperture position v is given by, Where o(r) is the object reflectivity function, V is the illuminated volume and k = ω/c is the wavenumber composed of the temporal radial frequency ω and the vacuum speed of light c.Amplitude variations due to antenna patterns and wave propagation attenuation of the volume elements (voxels) specifying the object are included in o(r).r mu is the path of the electromagnetic waves travelling from an arbitrary T x to a an arbitrary voxel located at r = (x, y, z) and back to an arbitrary R x , The aspect angle between a voxel and an effective aperture element is denoted by θ y .Through combining the effective aperture of the physical multistatic line array and the synthetic aperture generated by the movement of the object we generate a 2-D sampling aperture containing N E × N S sampling positions with N S the number of the synthetic aperture sampling positions.The required N S is estimated using the spectral support along the synthetic aperture as shown in Sect. 4. Range focusing is obtained using a stepped frequency continuous wave (SFCW) modulation with N F frequency points.The goal of an imaging algorithm is to retrieve the object function o(r).A common method to solve this problem is to use a spatial variant Matched-Filter (Pastorino, 2010), where now u comprises all N E measurement positions, v all N S synthetic aperture positions and k all N F SFCW frequency points.The reconstruction of 3-D images using a direct implementation of Eq. ( 4) has a high computational load, which makes it unsuitable for time critical imaging applications.By setting M = N E = N S = N F , the reconstruction of N 3 voxels using a direct implementation of Eq. ( 4) has a computational burden of O(N 3 M 3 ).The computational costs of the implementation of Eq. ( 4) can be reduced using the following algorithms.

The BP algorithm
Assuming that the measured signal stems from a set of point sources distributed within a volume to be inspected, this algorithm projects the measurement data back to their sources.Details on the algorithm can be found in Ulander et al. (2003).Here only the main implementation steps are mentioned.Starting from a SFCW signal the first step is to inverse Fourier transform the measurement data along the frequency axis Where ρ is the range variable and numerically performed using the FFT.In this way the data are range focused.The image is reconstructed by estimating the reflectivity of each voxel through linear interpolating the range focused data according to its relative coordinates to the sampling aperture, with a = (v, u, Z o ) at N E positions u and N S positions v.

The FFBP algorithm
This is a numerically optimized implementation of the BP algorithm.The optimization is based on reducing the computational costs of the summation in Eq. ( 6).Through iteratively dividing the total volume in sub-volumes and assigning samples from adjacent aperture positions to the centers of the sub-volumes the total number of summations in Eq. ( 6) is reduced.The performance of this algorithm is governed by a so called factorization factor n, which denotes the number of adjacent aperture positions taken to form a single new aperture position.A more detailed description can be found in Ulander et al. (2003).

The modified RMA
From Eq. ( 4) we see that o(r) is estimated through 3-D convolution in the aperture coordinates (u T , u R , v), which can be implemented in the wavenumber-domain using a complex multiplication.This is the basic idea of the RMA.However the thinning of the T x -array violates the Nyquist sampling criterion along the axis of the T x -array and leads to aliasing in the wavenumber-domain.This can be modified to be suitable for imaging with sparse arrays.In the following the main steps of the modified RMA are summarized.Since in the array's near-field the effective aperture approach is approximately valid due the spherical phase front, a phase compensation factor is obtained by calculating the difference between the measured phase and the intended far-field phase in the volume center.Then the measured signal is Fourier transformed along the sampling aperture, The next step is a regridding of S(k x , k y , k) known as the Stolt interpolation (Cafforio et al., 1991), with This can be obtained by evaluating the Fourier transform in Eq. ( 7) using the method of stationary phase, which is discussed in more detail in Appendix A. A first estimate of the reflectivity function is obtained by inverse Fourier transforming the 3-D spectrum with respect to k x and k y , The last step is to correct the range curvature in the first estimation due to near-field operation, 4 Sampling constraints and resolution For a correct reconstruction of the object sampling constraints in the space domain have to be fulfilled.These constraints and the expected spatial resolution are discussed in the following paragraphs.For an extended discussion we refer to Soumekh (1994).

Spatial sampling with a physical aperture
The spacing d R is determined with the knowledge of the spatial frequency (wavenumber) support y of the recorded signal with respect to y, Where k y,max and k y,min are the extrema of the y component k y of the wavenumber k and are given by, . (15) The expected spatial resolution along the y axis is also a function of the spatial frequency support y .The bandwidth of the spatial spectral support depends on the relative position of a voxel to the array.For the center of the object the spatial resolution is given by 4.2 Spatial sampling with a synthetic aperture Along the x axis the translation movement is used to generate a synthetic aperture in this direction.Sampling constraints are also derived from the spatial frequency x of the recorded signal with respect to x.The maximum step size is then given by Where the half-power beamwidth of a T x /R x antenna is denoted by θ 3 dB .Using a synthetic aperture we can ensure a constant spatial spectral bandwidth for all parts of the object.
The expected spatial resolution in this direction is given by

Frequency sampling and range resolution
The range resolution is determined by the signal modulation bandwidth B. For free space wave propagation it is given by Since we use a SFCW waveform, a maximum temporal frequency step size f is required in order to avoid range ambiguities.For a required range unambiguity R un , 5 Comparison of the algorithms

Focusing quality
Firstly, we compare the focusing quality of the BP algorithm, FFBP algorithm and the RMA using numerical simulations.Assuming that the required spatial resolution is 4 mm, according to Eqs. ( 12) and ( 16), an effective aperture of 50 cm is required at an imaging range of Z o = 60 cm.The spacing between the receiver elements has to be at most 2 mm in order to avoid aliasing.Hence we assume a R x -array with 10 elements spaced by 2 mm and a T x -array containing 25 T x with an element spacing of 2 cm to obtain the required effective aperture.The corresponding number of effective elements is 250 elements.We also assume a synthetic aperture of 25 cm.According to Eq. ( 17) the sampling step along the synthetic aperture has to be at most 1 mm.We use an SFCW waveform with 300 frequency steps from 75 to 110 GHz.We simulate the signal from a 3-D grid of ideal point reflectors with a uniform spacing of 5 cm using Eq. ( 2) and reconstruct it using the proposed RMA, the BP algorithm and the FFBP algorithm.
The implementation of the two latter algorithms requires the definition of a discrete 3-D rectangular volume, which covers the expected extent of the target.The distances between each two points of the volume is set to half of the expected resolution along each dimension.Two perpendicular sections from the image reconstruction obtained using the BP algorithm are depicted in Fig. 2. Figure 2a shows an xysection of the volume, which describes the estimated object reflectivity at a constant depth z = Z o .This section shows a successful focusing of the data along the x axis and the y axis with respect to the expected resolution.The section depicted in Fig. 2b is a yz-section obtained at x = 0 and demonstrates the focusing quality along the z axis.The corresponding xy-  section and yz-section, which result from the FFBP algorithm with a factorization factor of n = 2 are shown in Fig. 3a  and b, respectively.They also demonstrate the successful 3-D focusing using this algorithm.Finally the same sections, which result from the proposed RMA are provided in Fig. 4 to illustrate its 3-D focusing quality.
For a quantitative assessment of the imaging quality, we take a closer look at the focusing along the y axis and x axis of each algorithm.The upper part of Fig. 5 shows three y profiles of the estimated reflectivity function at x = 0 and z = Z o .The profiles resulting from the BP and FFBP agree very well.The 3 dB-width of the main lobes yield a width of ca. 4 mm.The side lobes level is ca.−14 dB due the rectangular shape of the aperture.The width of the main lobes of the profile resulting from the RMA fulfill the expected resolution, yet with increasing y-coordinate the RMA main lobes yield a slight shift accompanied by a broadening.An increase of the side lobes level of the RMA by ca. 3 dB is also observable.The x-profiles of the algorithms at y = 0 and z = Z o are depicted in the lower part of Fig. 5. Along with the x axis the 3 dB-width of the main lobe of the three algorithms agree very well.All algorithms achieve the desired resolution of 4 mm.Along this axis the RMA yield solely a slight shift with increasing x coordinate.The behavior of the RMA is expected, due to the approximately validness of the effective aperture approach in the near-field of the array.It is worth mentioning, that although the synthetic aperture is as large as half of the effective aperture, both achieve the same spatial resolution.This is due to the spectral support of the synthetic aperture, which is two times higher than the one of the effective aperture, compare Eqs. ( 16) and (18).

Computational complexity
Secondly, we compare the computing costs of the three discussed algorithms.Using the asymptotic computational complexity of the implementation steps of each algorithm we can estimate its computational burden.For the RMA the highest computational burden is caused by the 3-D Stolt interpolation.Since this is implemented using a complex 3-D linear interpolation, its asymptotic operation count C p is Additionally we use in this algorithm 5 complex FFT operations with a total asymptotic operation count C fft of, Hence the asymptotic computing complexity of the RMA is, The computational burden of the BP algorithm is estimated from Eqs. ( 5) and ( 6), For an aperture factorization factor n the computational load of the FFBP algorithm is given as function of the computational load of the BP algorithm (Ulander et al., 2003), For large volumes N 3 M 1 2 log 2 (M) , Comparing the FFBP with the RMA we obtain a computational saving in the order of By setting N = M we obtain the computational saving as function of amount of measurement samples, Furthermore, we have compared the performance of the mentioned algorithms on a commercial i7 CPU with 3.5 GHz and 16 GB of RAM.The scientific Python modules (van der Walt et al., 2011) have been used for the implementation of the algorithms.Using numerical simulation we investigate the dependency of the reconstruction time on the amount of voxels to be reconstructed for each algorithm.The results are depicted in Fig. 6.The x axis of Fig. 6 represents the total number of the image voxels (N 3 ).As expected the highest reconstruction duration is caused by the BP algorithm.This increases rapidly with increasing number of voxels.A significant improvement is obtained using the FFBP algorithm.
For instance, the reconstruction of 100 Kilovoxels becomes four times faster.However, the modified RMA yields the best computational performance.The 100 Kilovoxels have been reconstructed 22 times faster than the FFBP algorithm and 90 times faster than the BP algorithm.Significantly faster reconstruction times can be obtained by parallelizing the implementation of the algorithms as reported in Baccouche et al. (2017), where the mentioned algorithms are executed on a graphics processing unit.

Measurement results
We experimentally investigated the imaging performance of the algorithms using measurement data from the imaging system presented in Baccouche et al. (2015), which combines a sparse line array with a band-conveyor for 3-D tera-hertz imaging as sketched in Fig. 1.This system operates in the W -band and provides a modulation bandwidth of 35 GHz.As Fig. 7 shows, the multistatic sparse line arrays used in this system contains 12 transmitters, which are linearly distributed along the y axis.Standard-gain horn antennas are used for the transmitters as well as for the receivers.The transmitters are sequentially operated using a switching matrix.The R x -array contains also 12 elements.For the ease of realization, the elements of the R x -array are grouped in two sub-arrays of 6 elements and placed slightly shifted to each other on both sides of the T x -array with an offset of 4.5 cm.Measurement data are obtained using the full system bandwidth from 75 to 110 GHz and an effective aperture of 144 sampling positions with an effective spacing of 5 mm.This effective elements spacing, which is imposed by the mechanical dimensions of the used components, violates the Nyquist sampling criterion according to Eq. ( 15).Hence an increase of the grating lobes level is expected.For the sake of clarity, in the ideal case, the violation of the Nyquist sampling criterion within the sampling aperture is compensated by the effective aperture approach.As measurement sample we used an A-sandwich (20 cm × 10 cm × 0.5 cm) made of glass fiber reinforced plastic (GFRP) with a Rohacell core, depicted in Fig. 8a.The sample was prepared with different insertions at different locations with varying depths on the top of the Rohacell core. Figure 8b shows a schematic of the sample with the different insertions and defects, which are made of adhesive, teflon and polyethylene (PE).Furthermore, a step wedge was formed at the top of the sample through removing surface layers.The step size is less than 1 mm.The corner of the sample was marked with an aluminum quadrant.The imaging range is ca.65 cm. Figure 9 shows different 2-D layers from the 3-D image reconstruction of the sample using the BP algorithm as a reference (Fig. 9a, c and e) and the modified RMA (Fig. 9b, d and f).By comparing Fig. 9e and f, the same amount of defects at this layer can be recognized using both algorithms even if the modified RMA yields the expected degradation with respect to the array's axis (y axis).Figure 9c and d shows that the modified RMA is also sensitive to the tiny thickness variations of the sample as well as the BP algorithm.The upper and the lower layers of the sam-ple are resolved by both algorithms (Fig. 9a and b).These two figures show also the expected increase of grating lobes along the array's axis, yet confirming that the modified RMA approximates well the focusing quality of the BP algorithm.

Conclusions
We have compared the back-projection algorithm (BP) and its numerical optimized implementation, the fastfactorization back-projection (FFBP) algorithm, with a modified range migration algorithm (RMA) for 3-D terahertz imaging using a sparse multistatic line array in combination with a synthetic aperture.As criteria for comparison we took the focusing quality and the computational complexity.We have shown that the computational costs can be significantly reduced using the FFBP algorithm with competitive imaging quality to the BP algorithm.On the other side the modified RMA yields the highest computational efficiency at the costs of a minor focusing degradation along the array's axis due to the operation in the array's near-field.We used an undersampled sparse multistatic line array to inspect the inner structure of a glass fibre reinforced plastics (GFRP) sample in the frequency range 75 to 110 GHz.Using the modified RMA we can recognize the same relevant features within the sample as by using the BP algorithm with dramatically reduced computational effort.
Figure 1.Schematic of a generic imaging setup using the effective aperture concept in combination with a synthetic aperture.

Figure 5 .
Figure 5. Focusing quality comparison of the modified RMA with the BP and FFBP algorithms.

Figure 6 .
Figure 6.Reconstruction time comparison of the modified RMA with the BP and FFBP algorithms.

Figure 7 .
Figure 7. Multistatic sparse line array composed of 12 T x and 12 R x .

Figure 8 .
Figure 8. A-sandwich GFRP with a Rohacell core with a size of (20 cm × 10 cm × 0.5 cm): (a) Photograph of the sample.(b) Schematic and distribution of defects.