Phase Shift Migration With SIMO Superposition for MIMO-Sidelooking Imaging at Terahertz Band

In this paper, the enhanced phase shift migration (PSM) algorithm for multi-input-multi-output (MIMO) side-looking (SL) scheme was proposed for the 2-dimension (2D) fast image reconstruction in the terahertz (THz) band. By decomposing the MIMO imaging problem into a superposition of several single-input-multiple-output (SIMO) sub-arrays, the proposed algorithm can be applied to systems with randomly distributed transmitting array. By deriving and applying the modified phase shift factor and necessary phase compensation, zero-padding in the transmitting array and data reorganization in the spatial-frequency domain of the conventional wavenumber domain MIMO algorithm can be totally avoided, which will hopefully further reduce the amount of calculation. A comprehensive MIMO-SL-EPSM algorithm was established in the wavenumber domain. Proof-of-principle experiments in the 0.1-THz band were performed based on a MIMO-SL prototype system. The reconstructed images for different targets with high efficiency and quality demonstrate the theoretical results and the effectiveness of the proposed algorithm.


I. INTRODUCTION
Generally, THz waves are referred to the spectrum from 0.1 to 10 THz, which lie in the gap between the microwave and infrared. In recent years, many THz imaging systems have been designed and developed for various applications, particularly for nondestructive testing (NDT) and security inspection [1]- [15].
The imaging frame rate and system cost are generally the key factors in the practical application which should be considered in the development of terahertz imaging technology. Early terahertz antennas were difficult and expensive to manufacture, therefore, imaging systems were often developed based on a single transceiver combined with a mechanical scanning scheme, such as a quasi-optical element The associate editor coordinating the review of this manuscript and approving it for publication was Wei Wang . combined with a rotatable mirror for beam focusing and spatial scanning [6], [12]- [15]. In addition, the method combining the concept of synthetic aperture with quasi-optical scanning was proposed to integrate the fan-beam imaging system in the THz band [2].With the reduction of the cost of terahertz devices and the development of MIMO technology [16], [17], the active THz-MIMO technology has caught more and more attention [18]- [21].
Compared with conventional synthetic aperture radar (SAR), spatially distributed multiple transceivers applied in the MIMO radar system can not only significantly reduce the used transceivers and the data-acquisition time but also capitalize on the diversity of target illumination and scattering by viewing the target from multiple aspects [22]. In the MIMO scheme, due to hardware limitations, the maximum number of transmitting and receiving channels is usually limited. In order to further improve the azimuth resolution through a given number of hardware channels, the bistatic sparse array topology is usually employed [22]. Due to the complexity of the MIMO system with sparse array geometry and the coupling of the transceiver array in the wavenumber domain, the traditional FFT based SAR mono-static imaging algorithms, such as the fast algorithm employed in [14], cannot be directly applied to the MIMO system with the bi-static scheme. Therefore, the imaging result of the MIMO systems is conventionally reconstructed based on the back projection (BP) algorithm [7], [24], which is based on time-domain signal phase compensation and multi-channel coherent summation. However, the BP algorithm is always time-consuming and not applicable to real-time imaging applications.
In recent years, researchers are devoted to the study of fast BP algorithms [25]- [27], however, the approximation is usually introduced, which will increase the limitation of BP algorithm and reduces the imaging quality. The range migration algorithm (RMA) applied for MIMO systems [28] based on the Fourier Transform (FT) fully utilizes the efficiency of the FFT algorithm and had been widely used. However, accurate interpolation and data reorganization are necessary for practical applications; moreover, the quality of the interpolation directly determines the final image quality. As we all know, the interpolation process usually takes up most of the computing resources and affects the imaging speed. And the truncated interpolation kernels may cause distortion, especially for edge targets. In addition, RMA imposes strict requirements on the geometry of array, which limits its range of applications. Therefore, it is still a challenging task to develop a fast imaging algorithm suitable for MIMO systems.
Considering the limitations of the RMA, in our previous work [3], [29], the phase shift migration algorithm originated from reflection seismology was firstly extended for the MIMO imaging scheme. MIMO-PSM algorithm based on the explosion reflection model was proposed for 2-D imaging and achieved imaging quality comparable to BP algorithm with imaging efficiency far faster than BP algorithm. However, in order to prepare the data for interpolation where spatial frequencies over transmitting and receiving arrays must be fused into the same grid of data, the sampling steps of spatial frequencies for both arrays shall be the same. That means, the transmitting and receiving arrays must have the same spatial length and uniform spatial steps. Therefore, the zeropadding operation of both transmitting and receiving arrays is always required prior to the FFT operation [28].
In this paper, we further research and optimize the MIMO-PSM imaging technology in the THz band. By decomposing the MIMO imaging problem into a superposition of several SIMO sub-arrays, the SIMO compensation operator is proposed to avoid data reorganization in the wavenumber domain. In addition, combined with the PSM algorithm, frequency-wavenumber domain interpolation is avoided, further eliminating the influence of the interpolation kernel on the imaging results. The proposed methods and the algorithms were verified by the proof-of-principle experiments in THz band.  This paper is organized as follows. Details about the imaging scheme are presented in Section II. In Section III, the principle and formula of the Enhanced PSM algorithm are introduced and gives the implementation steps of the algorithm in practical application in detail. In Section IV, some simulation results in the 0.1-THz band are given to demonstrate the performance of the proposed imager. And proof-of-concept experimental results are given in Section V. Finally, a conclusion is drawn in Section VI.

II. DESCRIPTION OF THE MIMO-SL IMAGING SCHEME
The SL imaging scheme is shown in Fig.1, and the geometry of the typical MIMO 1-D linear array with both sides transmitting and middle receiving is shown in Fig.2. The targets in the scene are located in the x-y plane, and the MIMO array with height h above the ground lies in the azimuth direction (x-direction). The beams illuminate the imaging scene with the angle θ look of SL mode, which represents the angle between the center direction of the beam and the negative z-direction. The swath width represents the range of the imaging scene along the y-axis, which depends on the beam width and the SL Angle. The beamwidth in the azimuth direction determines the range of the imaging scene along the x-axis.

A. EXPLODING REFLECTOR MODEL FOR IMAGE RECONSTRUCTION
In this paper, the phase shift migration originated from reflection seismology is the basis of the imaging algorithm. VOLUME 8, 2020 The fundamental theory is based on the exploding reflector model [30], [31] in geophysics, in which all the targets are assumed to 'explode' at the time t = 0 and the 'exploded fields' propagate with velocity c 2 toward the receiver. The literature [29] pointed out that the explosion field of the 2-D imaging problem can be expressed as (1), (2), (3): where U (x, y, ω) is the exploded field at (x, y) with the angular frequency ω andŨ (k x , y, ω) is the spatial wavenumber domain function of U (x, y, ω) . Through the explosion reflection model we can find that the explosion field at t = 0 has a similar form as the objective function to be recovered [11], and the 2-D objective function can be obtained by the integration of angular frequency ω , specifically expressed as: (4) The traditional PSM algorithm records the exploded field in a monostatic array, which cannot be directly applied to the bistatic MIMO radar. In our previous work [29], based on the bistatic echo data, a reorganization method in the wavenumber domain was developed to derive the explode fields in the ''equivalent frequencies'' at time t = 0. Finally, the PSM algorithm was successfully applied to imaging reconstruction in the bistatic MIMO array.

B. DERIVATION OF THE ENHANCED PSM ALGORITHM
In this Section, the relationship between the ''explosive field'' and echo data of the SIMO sub-array was firstly derived. By decomposing the MIMO imaging problem into a superposition of several SIMO sub-arrays, it is feasible to make the receiving array satisfy the Nyquist space sampling criterion only. Necessarily, the SIMO compensation factor must be applied to the explosion field of each sub-array. Then, the corresponding explosion fields at time t = 0 can be successfully recovered. Finally, the image reconstruction of MIMO-SL is realized. Based on the above theories and methods, a comprehensive method of Enhanced PSM for MIMO-SL imaging is finally developed, which can realize MIMO-SL imaging of arbitrary transmitting element positions. Details are as follows.
Based on a realistic THz imaging geometry given in Fig.1 the echoed signal model s(x t , x r , k) can be obtained. The side-looking scene is shown in Fig.3 and the MIMO array is located at the plane y = 0 . The arbitrary transmitters of the MIMO array are located at (x t , 0) ; the uniform receivers are located at (x r , 0) . The following derivation process takes the first sub-array as an example. Under the first-order Born approximation, the received wave field corresponding to the transmitter x t1 can be extracted where O(x , r ) represents the objective function of the target over the x-r plane. k = 2π f /c represents the frequency wavenumber that is directly related to the frequency f . And c represents the speed of the electromagnetic wave, which is generally equal the speed of light. R T 1 , R R are the distances from the target located at (x , r ) to the first transmitter and all receivers, can be expressed as (6) and (7), respectively By ignoring the amplitude decay factor which causes negligible influence on reconstructed image [11], (5) can be simplified as Fourier transform was firstly implemented on (8) over x r , to transfer the extracted echo data s 1 (x r , k) into the spatial frequency s 1 (x r , k) , as given in (9) (9) where F (k xr , k) represents the Fourier transform of the exponential exp (−jkR R ) , can be expressed as The integral in (10) can be solved by the Method of Stationary Phase (MSP), as given in (11) Then, substitute (11) into (9) From (12), the 2-D dispersion relations are defined as By comparing (12) and (2), we can interestingly find that a similar form of the 'exploding fields' can be derived by applying the above dispersion relationship for variable substitution. It can be seen from (13), the spatial frequency k x , k r are only related to the spatial frequency k xr of the receiving array, which indicated that SIMO decomposition achieves the decoupling of the spatial-frequency domain of the transmitting and receiving array.
Based on the dispersion relations in (13), the ''explosive field'' in the frequency-wavenumber domain with the bistatic SIMO sub-array data can be derived as The 'exploding fields' derived from the bistatic imaging topology have the extrapolation property in r-coordinate [29]. That means, the 'exploding fields' at different distances r i = r 0 , r 1 , . . . , r n−1 can be obtained by phase-shifting the explosion field at r = 0 . The specific operation is given as where exp(jk r r i ) is the extrapolated phase shift factor of the PSM algorithm along the propagation direction of the electromagnetic wave. This is similar to the monostatic imaging scheme [11], So it can be used for more complex MIMO-SL imaging. After the 1-D Inverse Fourier Transform over k x , the ''approximate exploding field'' in each distance can be obtained. To get an accurate exploding field, the residual phase has to be compensated. The specific operation can be expressed as The image in the x-r coordinate system can be reconstructed by the integral of (16) over frequency to recover the ''explosive field'' at time t = 0. To facilitate understanding, the wavenumber k can be substituted by the angular frequency ω u 1 (x, r, t)| t=0 = u 1 (x, r, ω) · exp(jωt)dω (17) The scale relationship between r i and y i in the Cartesian coordinate system is shown in (18). Then, the integrated formula for imaging with the proposed algorithm can be expressed as After scale-transformation, the imaging result is mapped from x-r plane to x-y plane where targets are located. Finally, the imaging result can be obtained by accumulating results of the SIMO sub-arrays: where, N xt represents the actual number of transmitters. We can find that the above formulas do not assume the positions of the transmitters. Therefore, (19) and (20) can be used for the MIMO linear array imaging system with arbitrarily distributed transmitters. According to the principle of reciprocity, the proposed algorithm can also be applied to the array with arbitrarily distributed receivers and uniformly distributed transmitters. For the convenience of discussion, we assume that the receivers are uniformly distributed in this paper.

C. SAMPLING CRITERIA & SPATIAL RESOLUTION
Theoretically, to obtain the accurate MIMO imaging results using SIMO superposition, it is required that any aliasing of the target should not arise in the imaging results of each SIMO sub-array. The analysis of the azimuth performance of the SIMO aperture is actually similar to a one-way aperture [32]. Defining that the distance between the target and the center of the receiving array is R, the length of the target along the range is L r , and the length of the target along azimuth is L a , the frequency sampling interval and the spatial sampling interval of the receiving array must satisfy the Nyquist criterion are given by where L xr is the aperture length of the receiving array, λ min is the wavelength corresponding to the highest transmitting VOLUME 8, 2020 frequency, and it can be found that there is no constraint on the aperture of the transmitting array in the above formula. The theoretical spatial resolutions are shown in (23) and (24), where B represents the bandwidth of the transmitted signal. With a fixed system frequency and detection range, the resolution of the array direction (x-direction) is related to the transmitting and receiving arrays aperture [22]. For the resolution of the range direction (y-direction), under the SL regime, the angle between each horizontal distance line and the array is inconsistent, so it depends on the target distance in the y-direction and the height h of the array.

D. IMPLEMENTATION
The detailed processing of the proposed EPSM algorithm is shown in Fig.4. By applied FFT and IFFT, the algorithm can achieve high computation efficiency. The proposed algorithm avoids the zero-padding of the transmitting array, we can easily extract the echo data corresponding to the SIMO sub-array for processing. After the FFT over x r , the wavenumber domain data U 1 (k x , r i , k) at the specific distance r i can be estimated by applying the phase shift factor exp(jk r r i ),where k x and k r are defined in (13).
Then the sub-images at each frequency point can be obtained after 1D-IFFT over k xr and introducing the SIMO compensation factor which is used to compensate the residual phase. Obviously, the 2-D focused image in x-r coordinates can be obtained by integration over the frequency domain, as shown in (16).
Finally, the reconstructed image on the x-y plane can be obtained after scale-transformation and accumulating results of all sub-arrays. And the imaging result needs to be uniform in both (x, y) dimensions, which means the imaging result in the space (x, r) is nonuniform along the r-coordinate according to (18). However, thanks to the independence of the PSM algorithm at different distances as shown in (15), we can apply the nonuniform samples in r-coordinate to get the uniform samples in y-direction perfectly and precisely without any interpolation processing.

E. COMPUTATIONAL COST
In order to quantitatively evaluate the computational complexity of the proposed reconstruction algorithm, it is assumed that the echo data is recorded at sampling points in the frequency domain, the phase shift factor can be stored in a register in advance for easy recall, and the calculation cost of the phase shift operation can be ignored in application. Through a comprehensive analysis of the imaging formulas, the calculation cost of the algorithm can be expressed as the total number of real multiplications and real additions: (25) where N f is the frequency steps number of the SFCW signal, and N r is the number of samples along r-coordinate. t mul and t add are the computation time of real multiplication and real addition for a specific computer, respectively. The detailed computation cost is shown in Table 1.
In our previous work [29], the computational cost of the MIMO-PSM algorithm had been compared with the BP algorithm. In order to show the computational efficiency of the proposed algorithm, we compare it with the MIMO-PSM algorithm. The computational cost of the MIMO-PSM algorithm can be expressed as (26). (26) where N x = N xt + N xr − 1, it should be noted that N xt and N rt in the formula no longer represent the actual number of transmitters and receivers, but the number after zeropadding. In the MIMO array that transmitters at both sides and receivers in the middle, such as the array shown in Fig.5, the numbers of transmitters and receivers after zero-padding are always much larger than the actual numbers. Therefore, the proposed EPSM algorithm can obtain much faster computational efficiency than the MIMO-SL-PSM algorithm. The great improvement of the efficiency and will be validated  further with simulations and experiments in the following parts.

IV. NUMERICAL SIMULATIONS
The numerical simulations based on the proposed algorithms are demonstrated in this section. The simulation results of BP algorithm and the proposed algorithm are compared by point spread function (PSF). After that, the computational costs of the proposed algorithm and the MIMO-PSM algorithm are quantitatively compared with the same simulation parameters. The geometry of MIMO array shown in Fig.5 is employed in this section. The linear MIMO array composes of 6 transmitters and 39 receivers. The total length of the linear array is 0.3 m with 2.5 mm transmitting array sampling interval and 7.5 mm receiving array sampling interval. The frequency of the signal ranges from 75 GHz to 110 GHz with 201 points sampled in a step of 175 MHz. In the current simulation, the θ look shown in Fig.2 is set to 60 • . Through the theoretical model, bistatic MIMO simulated echo data can be obtained. Then, the MIMO-SL-EPSM algorithm developed in Section III. can be used to efficiently reconstruct 2-D images.
The simulated target consists of an array of 17-point scatters distributed within an area of 0.2m × 0.5m in the x-y plane. The center point of the target is located at y = 3.03m. The 2-D reconstructed images are given in Fig.6. For comparison, the images reconstructed from the proposed algorithm and the conventional BP algorithm are both given, which all get fully focused. The red circles represent the real positions of the targets.
To illustrate more clearly, the 1-D imaging results of the center target along the line y = 3.03m and x = 0m are shown in Fig.7(a)-(b), and the 1-D imaging results of the edge target along the line y = 2.79m and x = − 0.09m are shown in Fig.7(c)-(d), from which a well consistency between the proposed algorithm and BP algorithm can be clearly observed. Taking the width of the main lobe of -3dB as the resolution index, the simulation results show that the azimuth resolution is 1.42cm at y = 2.79m, which is near with the calculation result of (22) of 1.44cm, indicating the correctness of the theoretical analysis.
While as shown in Table 2, Under the same simulation parameters, the computational cost of the proposed algorithm is less than that of the MIMO-SL-PSM algorithm and much less than that of BP algorithm.

V. EXPERIMENTAL RESULTS
To illustrate the performance of the proposed MIMO-SL-EPSM algorithm in practice, as shown in Fig.8, a prototype imager with a 0.1 THz band was developed. The complete system exhibits a bandwidth from 90 to 110 GHz with a frequency step of 100MHz. The system was developed based on a microwave Vector Network Analyzer (VNA) combined  with two transmitters. The transmitters generate a 0.1 THz transmit signal, and then the receiver converts the 0.1 THz signals to the VNA intermediate frequency (IF) to extract the amplitude and phase information of the echo data. The triaxial motion controller as shown in Fig.8 can control the motion of three motors independently and precisely.
In this paper, the controller is implemented as a receiver controlled by a central motor, and two transmitters controlled by the motors on both sides to realize the equivalent one-dimensional MIMO array. We controlled the central receiver to cross 39 different positions in the center of the equivalent array, and control the two transmitters to cross 3 different positions on each side, to achieve an equivalent MIMO 1-D array with 39 receivers and 6 transmitters consistent with simulation. After recording all the data, the obtained scattering parameters in the 0.1-THz band are transmitted to the computer for data processing and image processing.
As shown in Fig.9, two experimental results are given. The target used in the experiment was made of foam material and covered with copper foil on the surface. The center of the target is located at 1.2 m away from the plane of y = 0m. By placing several thin foam materials under the end of the target, the target is placed at an angle of approximately 35 • to the ground to form an equivalent SL angle of approximately 60 • . Obviously, both algorithms successfully reconstruct the targets with high quality and the results of different algorithms exhibits good agreement with each other.
As shown in Table 3, the proposed EPSM algorithm requires only about 5.074 × 10 8 real multiplications and 3.132 × 10 8 real additions while MIMO-SL-PSM requires 13.679 × 10 8 real multiplications and 14.515 × 10 8 real additions. For a Lenovo PC with 4 Intel core i5-2400 CPUs, the entire reconstruction requires less than 0.7 s which is less than the time consumed by MIMO-SL-PSM, and undoubtedly, much less than the time consumed by BP.
Limited by the preliminary experimental conditions of the laboratory, we set up a small-scene experimental scheme to verify the high efficiency and high quality of the proposed algorithm. Theoretically, by increasing the length of array and  the number of frequency steps, the proposed algorithm can be applied to the security inspection of human body. And by modifying the N r , N xr and N f in (25), the calculation cost of the proposed algorithm under the specified imaging scene can be obtained.
Here we add the computation cost analysis for the practical scene to cover an entire human body. The human body scene is set to 0.7m × 2m and the imaging parameters are set according to the above-mentioned resolution requirements. After evaluation, the following modified parameters can be obtained.
For the MIMO-SL-PSM algorithm, the following modified parameters can be obtained. N xt = 301; N xr = 101; N r = 301; N f = 301 (28) In the human body imaging scene, the computational costs computed by (25) and (26) are shown in Table 4. Obviously, the computational cost of the proposed algorithm is quite lower than that of the MIMO-SL-PSM algorithm.

VI. CONCLUSION
In this paper, the enhanced PSM algorithm for MIMO-SL imaging in THz band is studied. By applying the SIMO compensation factor to the explosion field of each transmitter, the analytical expression of the ''explosion field'' of the corresponding SIMO subarray is derived, and the corresponding explosion field at time t = 0 is successfully recovered. Finally, a comprehensive MIMO-SL-EPSM algorithm is developed. Compared with the conventional BP algorithm, the proposed algorithm takes advantage of the fast wavenumber domain transformation, reducing the time consumption greatly without deteriorating the imaging quality. Since the PSM algorithm does not require interpolation, there will be no edge loss due to truncated interpolation kernels in the imaging results compared with the RM algorithm. Compared with MIMO-PSM, the proposed algorithm avoids the time-consuming data reorganization in the spatial-frequency domain and transmitting array zeropadding, which further reduces the amount of calculation.
Furthermore, the proposed algorithm can be applied to the MIMO system with a uniform receiving array and nonuniform transmitting array, which enhances the flexibility of the system greatly. A prototype MIMO-SL imager was designed for principle verification experiments in 0.1 THz band. The imaging results of different targets are given to verify the theoretical results obtained in this paper and the effectiveness and efficiency of the proposed algorithm.