Three-Dimensional Imaging of Terahertz Circular SAR with Sparse Linear Array.

Due to the non-contact detection ability of radar and the harmlessness of terahertz waves to the human body, three-dimensional (3D) imaging using terahertz synthetic aperture radar (SAR) is an efficient method of security detection in public areas. To achieve high-resolution and all aspect imaging, circular trajectory movement of radar and linear sensor array along the height direction were used in this study. However, the short wavelength of terahertz waves makes it practically impossible for the hardware to satisfy the half-wavelength spacing condition to avoid grating lobes. To solve this problem, a sparse linear array model based on the equivalent phase center principle was established. With the designed imaging geometry and corresponding echo signal model, a 3D imaging algorithm was derived. Firstly, the phase-preserving algorithm was adopted to obtain the 2D image of the ground plane for each sensor. Secondly, the sparse recovery method was applied to accomplish the scattering coefficient reconstruction along the height direction. After reconstruction of all the range-azimuth cells was accomplished, the final 3D image was obtained. Numerical simulations and experiments using terahertz radar were performed. The imaging results verify the effectiveness of the 3D imaging algorithm for the proposed model and validate the feasibility of terahertz radar applied in security detection.

1. Introduction

Terrorist attacks have increased substantially in recent years on a global scale, and they pose a great threat to public safety. Security detection in public areas is an important approach to guard against attacks in advance [1,2]. Compared with detection methods like metal detectors and manual checking, three-dimensional (3D) imaging techniques can provide more precise information, rendering automatic detection and identification possible. Recently, 3D synthetic aperture radar (SAR) imaging has attracted increasing research interest because of its ability to obtain the real spatial position and scattering properties of the targets  of interest [3,4,5,6,7]. On the other hand, terahertz (THz) waves have proved superior in their ability to penetrate clothes with inappreciable harm to the human body, and it is thereby suitable for concealed weapon detection. Besides, terahertz waves have a short wavelength and are able to generate wide bandwidths, which can achieve high-resolution imaging [8,9,10]. Therefore, 3D imaging using terahertz radar is a promising method in security detection.

Among various kinds of 3D imaging patterns, circular SAR (CSAR) is popular for its 3D imaging ability and full observation with a single pass [11,12,13]. CSAR 3D imaging needs a pitch angle to realize resolution in the height direction, which may limit its application in some cases. Besides, because of the pitch angle, resolutions between range and height directions are coupled, which makes it difficult to improve both resolutions simultaneously. To address this issue, much research has been conducted by scholars, including studies on multipass SAR, linear array SAR, etc. [14,15,16,17,18]. In [16], holographic SAR (HoloSAR) was applied to form 3D images, where it combines the height reconstruction with the 2D imaging of multipass CSAR. Compared with single-pass CSAR 3D imaging, HoloSAR can reduce sidelobes and provide more height details. However, HoloSAR may face the nonuniform height-sampling problem, which reduces the image quality. The sensor array is another choice for 3D SAR imaging, where the array distributes along the needed direction to provide a sampling aperture. However, to satisfy the half-wavelength spacing condition to avoid grating lobes, a large number of sensors are needed, which leads to a complicated and costly system, especially for terahertz radar. Therefore, we need to reduce the number of sensors to make the radar system realizable.

A sparse array was considered for 3D SAR imaging, however, the traditional Fourier-based imaging method may be invalid in this condition [19,20]. Fortunately, the targets of interest in security detection behave with spatial sparsity generally, thus, the sparse recovery method can be employed to form the image. Among the multiple sparse recovery algorithms, sparse Bayesian learning (SBL) has superiority to deterministic algorithms due to its characteristic of incorporating prior knowledge [21,22]. Without knowing the signal sparsity, SBL estimates the distribution parameters and the posterior with iteration, which can obtain more stable solutions. In Ref. [23], SBL was applied to through-wall imaging with missing data measurements, and in [24], a local structural SBL was adopted to form ISAR images of maneuvering targets with compressively sampled data. The desired results were obtained in both of the studies. Therefore, SBL was used for height reconstruction in this study due to its favorable performance.

In this work, the imaging geometry of CSAR with sparse linear sensor array was first established, and the echo signal model in the terahertz band was developed. Then, a 2D ground plane imaging procedure using phase-preserving algorithm was deduced for each sensor. After that, the phase features of the 2D images were analyzed, and reconstruction along the height direction was transformed into a sparse recovery problem. After completing the reconstructions of all the range-azimuth cells, the 3D image was obtained. The rest of this paper is organized as follows. In Section 2, the 3D imaging geometry with sparse linear sensor array is provided and the imaging procedure is described. In Section 3, simulations and experiments are presented and analyzed to demonstrate the effectiveness of the 3D imaging method for the proposed imaging model. Conclusions are given in the last section.

2. 3D Imaging Methodology

2.1. Imaging Geometry

There are several different transmitting-receiving models for a sensor array, such as single-transmitter multiple-receiver (STMR) and multiple-transmitter multiple-receiver (MTMR). Based on the equivalent phase center principle, it can be considered equivalent to the model that treats each sensor transmits signal and receives echo independently. The specific transmitting-receiving pattern design is not discussed in this paper, and the equivalent distribution of the sensor array is presented directly. The imaging geometry of CSAR with a sparse linear sensor array is shown in Figure 1, where the radar system consists of a vertical sparse linear array operating at a terahertz band. The sensor array spreads parallel to the Z-axis and moves along a circular trajectory with radius . Assuming that there are a total of M equivalent sensors whose Z-axis coordinates are , respectively, when the radar system moves to the azimuth angle , the instantaneous coordinate of the equivalent sensor m is . It is assumed that the target scattering is isotropic and the sensors can illuminate the region of interest during the entire sampling time.

Figure 1

Imaging geometry of circular synthetic aperture radar (CSAR) with sparse linear sensor array.

In general, a target can be approximated as a sum of several independent and non-directional scattering centers in SAR imaging [8]. Therefore, the echo of the targets can be expressed as the superposition of echoes of all scattering centers. Considering a general scattering center P with coordinate , which is plotted in Figure 1, the instantaneous distance from the equivalent sensor m to the scattering center is given as

With the assumption that the target is in the far-field, the fourth and higher-order terms can be neglected for the Taylor series expansion of the cosine term. The distance expression in Equation (1) can be approximated as where is the distance from the equivalent sensor m to the origin .

2.2. 2D Phase-Preserving Imaging Method

For terahertz radar, broadband can be realized, and the pulsed linear frequency modulated (LFM) signal using the dechirp-on-receive technique is appropriate for high-resolution imaging [6,8]. The transmitted radar signal is presented as where is the fast time, T is the pulse duration, is the center frequency, and is the modulated rate. The rectangular window is defined as for and otherwise.

For scattering center P, the echo received by sensor m is expressed as where is the scattering coefficient of P, and is time delay, where c is light speed.

For broadband signal, the matched filtering (MF) method will need a high sampling rate to satisfy the Nyquist sampling criterion, which could increase the data size remarkably. Nevertheless, the dechirp technique is able to overcome this difficulty and is used extensively in terahertz radar signal processing. The dechirped signal of Equation (4) is

In Equation (5), the second exponential term is the residual video phase (RVP), and it can be removed before the next steps [25].

Next, the signal will be transformed into wavenumber domain. Let , , and represent the signal wavenumber, center wavenumber, and bandwidth wavenumber, respectively. Substituting (2) into (5), we get where represents the pitch angle between the equivalent sensor and the plane and satisfies , .

In Equation (6), for a determined equivalent sensor m, the first phase term can be compensated with compensating factor . The wavenumber of K can be decomposed as follows

By substituting (7) into (6), the signal can be represented as

As mentioned in Section 2.1, the echo of the entire region of interest is the superposition of all the scattering centers. Thus, the radar echo received by the mth sensor is expressed as where is the scattering coefficient of a scattering center located at .

The triple integral in Equation (9) can be decomposed into two parts, one is a double integral about , the other is single integral about . Since the window function in (9) has no direct impact on the phase analysis, it can be neglected for convenience in the following. Then, the signal can be represented as

Note that the double integral in Equation (10) is the 2D Fourier transform, thus, the 2D inverse Fourier transform can be performed to both sides of (10). After the 2D inverse Fourier transform, the scattering distribution for each cell is obtained:

Equation (11) indicates that for a determined cell , the 2D imaging result is the correlation stack of the scattering coefficients of all the scattering centers with different heights z. Therefore, to obtain the 3D image, we need to reconstruct the height location of the scattering centers.

2.3. Scattering Coefficient Reconstruction along the Height Direction

To avoid the grating lobes, the dense linear sensor array needs to satisfy the half-wavelength spacing condition. However, the wavelength of a terahertz wave is as short as less than 1 mm generally, thus, integration in practical implementation is almost impossible. Therefore, the sparse sensor array along the Z-axis is implemented. However, this makes the Fourier-based signal processing invalid. Considering that there are limited dominated scattering centers in the height direction, the sparse reconstruction method can provide a solution to recover the scattering coefficient distribution. The data corresponding to the same pixel cell of the multiple 2D complex images are extracted and rearranged in a vector as the measurement data. This manipulation is shown in Figure 2.

Figure 2

Data corresponding to the same pixel are extracted and rearranged in a vector.

In discrete signal processing, Equation (11) can be rewritten as where Q is the number of grids divided along the height direction, and the values of the grids are , where denotes the transpose. Furthermore, Equation (12) can be presented as the vector form , where , , and represent the scattering coefficients of the targets at grids to be recovered. Therefore, for all , we have , where , . Since we will reconstruct the height information pixel by pixel, the subscript will be omitted for brevity in the following. With additive white Gaussian noise whose mean is zero and variance is , we obtain

Equation (13) can be exhibited in Figure 3.

Figure 3

Graphical representation of Equation (13). The blank squares indicate zeros, while the colored ones indicate non-zero values.

Since , the solution of Equation (13) is underdetermined. Generally, the targets of interest in security detection behave with spatial sparsity, which means that only a few values of shown in Figure 3 are meaningful. Thus, the solving of Equation (13) can be transferred to an optimization problem as where is a proper bound.

To solve the problem in Equation (14), sparse recovery algorithms can be adopted. In SAR imaging, the number of dominated scattering centers are generally unknown, thus, the algorithms based on the sparsity of the signal are impracticable. As analyzed in [22,23,24], sparse Bayesian learning (SBL) is an effective approach for sparse recovery. Compared with other -norm minimization algorithms, the SBL-based algorithm adopts iteration during the procedure to avoid convergence to the local minimum and produces a full posterior distribution as the solution. Based on the Gaussian scale mixture (GSM), the prior probability of the distribution of scattering coefficient is represented as [23]: where is the hyperparameter that governs .

Once the scattering coefficients are given, the probability of the measurement vector is determined by the noise, which is given by

With the prior probability of and the conditional probability of , the marginal probability of is yielded by where , is the identity matrix, .

Based on the Bayesian Theorem, the posterior probability of is where

From (18), we can see that the posterior probability is a Gaussian distribution with mean and variance , which are represented in (19). For Gaussian distribution, the mean value can be treated as the estimated scattering coefficient, i.e., the estimate . It can be seen that the estimate of the scattering coefficient is relative to the noise variance and hyperparameter . For traditional SBL, and can be estimated with expectation–maximization (EM) steps [22]. However, the noise variance is a nuisance parameter and the estimate may be highly inaccurate [22]. According to [23], can be integrated out by introducing a prior distribution, which is a gamma distribution given by where a and b are deterministic small value, and is a gamma function defined as . With the prior probability of , the posterior of is yielded as

The posterior distribution given by (21) is a multivariate complex Student’s t-distribution, whose mean and covariance are given by

Like Gaussian distribution, the mean value of Student’s t-distribution can also be treated as the estimated scattering coefficient, i.e., the estimate . We can see that the estimate is merely a function of the hyperparameter rather than the noise variance . Therefore, in order to obtain the estimate of the scattering coefficient, we have to estimate the hyperparameter . Using the EM approach in SBL, can be obtained by maximizing the marginal probability of , which is . The marginal probability is yielded by where .

Taking the logarithm of and neglecting the terms that are independent of hyperparameter , the cost function is represented as

Differentiating with respect to , and setting the results to zero, we can obtain where is the qth component of and is the qth diagonal component of .

From Equation (25), we can see that the hyperparameter is a function of mean and covariance . Meanwhile, Equation (22) indicates that and are functions of . Therefore, once the initial values of , a, and b are given, an iterative procedure can be implemented between (22) and (25) to obtain the estimates of , , and . Furthermore, the scattering coefficient estimate is obtained by . The processing procedure of the SBL mentioned above is shown in Algorithm 1.

Finally, combining the 2D imaging procedure with the reconstruction of the scattering coefficient along the height direction, the overall flowchart for 3D imaging is shown in Figure 4.

Figure 4

3D imaging flowchart.

3. Simulations and Experiments

3.1. Simulation Results and Analysis

In this section, we present numerical simulations used to verify the effectiveness of the proposed imaging model and sparse recovery method. In the simulations, the radar center frequency is set as 340 GHz based on the atmospheric transmission window theory. A benefit pf the high center frequency is that a large bandwidth can be achieved for terahertz radar, which is set as 28 GHz here. The other main simulation system parameters are also listed in Table 1.

Table 1

Main system parameters

Parameters	Values
Carrier frequency	340 GHz
Bandwidth	28 GHz
Pulse duration time	1 ms
Radar radius	5 m
Equivalent receiving antenna number	9
Aperture length in height direction	0.2 m

The simulated imaging scene is shown in Figure 5a, where there contains 10 point-scattering centers in three different planes. In Figure 5a, the points with the same color are located in the same plane. The 2D imaging method is first applied to each equivalent sensor to obtain the 2D image. After that, the sparse reconstruction algorithm is adopted to recover the scattering coefficients along the height direction for each range-azimuth cell. The final 3D imaging result using the sparse recovery method is presented in Figure 5b. Compared with Figure 5a, it can be seen that the point scattering centers are well reconstructed in the 3D space. The 3D imaging result using the MF method with sparse linear sensor array is also presented, which is shown in Figure 5c. It can be seen that the points with the same range-azimuth cell cannot be distinguished along the height direction. Moreover, the imaging result using the MF method with dense sensor array that satisfies the half-wavelength spacing condition is given in Figure 5d. It shows that the points are recovered with acceptable sidelobes.

Figure 5

Imaging scene and 3D imaging results. Dotted lines link the points in the same plane. (a) Original imaging scene; (b) imaging result using sparse recovery with sparse sensor array; (c) imaging result using MF with sparse sensor array; (d) imaging result using MF with dense sensor array.

To further exhibit the imaging results, the height profiles at the range-azimuth cell (0,0) are provided in Figure 6, where the red line, black line, and blue line correspond to the imaging results of spare recovery with sparse sensor array, MF method with sparse sensor array, and MF method with dense sensor array, respectively. In the original simulated scene, there are two points in the cell (0,0) with different heights. Figure 6 shows that the blue line has two main lobes at m and m, while the red line contains two evident peaks that conform to the main lobes. It should be declared that the sparse recovery is a super-resolution imaging method, thus the height profile has no sidelobes. On the other hand, the imaging result using MF with sparse sensor array is messy and has no prominent peaks. Therefore, the simulations demonstrate that the sparse recovery can reconstruct the height location of scattering centers with sparse linear sensor array.

Figure 6

Height profile of the range-azimuth cell (0,0).

3.2. Experiment Results and Analysis

With 0.34 THz radar, the equivalent experiment was conducted to validate the proposed imaging model and 3D imaging method. The system parameters are the same with those of the simulation shown in Table 1. The experimental scenario is presented in Figure 7. In Figure 7a, the radar is fixed while the targets are placed on a rotating platform, which is equal to the radar moving in a circular trajectory while the targets stand still. Since our radar system has no sensor array, we raise the platform step-by-step to form the relative height changing between the sensor and targets. The platform is precisely controlled to make sure that it rotates from the same original angle of each height step. Figure 7b is the optical photo of the terahertz radar with a single pair of transmitting-receiving antennas. More specified descriptions for this radar system, such as the schematic diagram and front-end setup, are presented in [8,9]. Figure 7c shows the optical photo of the targets, which contain five metallic balls in three different planes.

Figure 7

Experimental scenario. (a) Diagram of imaging scene; (b) photo of the terahertz radar system; (c) photo of targets.

The experimental results are shown in Figure 8, which contains three 2D slices corresponding to the planes, whose heights are m, 0 m, and m, respectively. The yellow circles in each figure mark the locations of the metallic balls. Comparing Figure 7c with Figure 8a, we can see that only the two balls in Figure 7c, whose relative heights are in accord with the height m, appear remarkably in Figure 8a. Similarly, this phenomenon is also presented in Figure 8b,c. It should be noted that the regions in red squares are noise generated by the crabsticks. Therefore, the experimental results demonstrate the 3D imaging ability of the proposed imaging model and sparse recovery algorithm.

Figure 8

2D imaging slices of heights at (a) m; (b) m; (c) m. The yellow circles mark the locations of the balls for each slice, while the regions in red squares are noise generated by the crabsticks.

4. Conclusions

For the application of security detection of the human body in public areas, a 3D imaging model that combines CSAR with sparse linear sensor array was established in this study. Based on the CSAR model, aspect 2D high-resolution imaging was implemented using a phase-preserving algorithm. Due to the short wavelength of terahertz waves, the sparse linear sensor array was used to form the resolution in height direction. With the obtained complex 2D images, the sparse recovery method, which is an SBL-based algorithm, was applied to obtain scattering coefficient reconstruction along the height direction. After accomplishing the reconstructions of all the range-azimuth cells, the final 3D image was obtained. Numerical simulation was performed to verify the high-resolution imaging ability of the proposed imaging model and algorithm. Furthermore, we used 0.34 THz radar to carry out the equivalent experiments, and the results demonstrate that the terahertz radar can be applied in high-resolution imaging, which indicates that terahertz radar imaging is suitable for security detection of humans in public areas. However, the imaging process described in this paper used the isotropic assumption of the targets, which may not always be true in practice, and the imaging method for aspect-dependent targets needs to be further researched.