A Millimeter-Wave 3D Imaging Algorithm for MIMO Synthetic Aperture Radar

Multiple-input-multiple-output synthetic aperture radar (MIMO-SAR) is being studied and applied in more and more scenarios. However, there is still a certain distance away from real-time imaging using advanced algorithms. The traditional backpropagation algorithm (BPA) multi-accumulation integration is unsuitable for dealing with large-size scanning data, and the wavenumber domain algorithm requires the array to satisfy Nyquist sampling law in azimuth to avoid aliasing in imaging reconstruction. Based on these issues, a novel 3D imaging method is proposed for MIMO-SAR. An appropriate transformation and inverse Fourier transform (FT) is carried out for the frequency domain; thus, accumulation in the wavenumber domain is not required, which is easy to implement. The computational complexity of the algorithm is much lower than BPA and has better generalizability than the wavenumber domain algorithm. Coherence factor (CF) is also introduced to achieve sidelobe suppression. Proof-of-principle experiments were also carried out in the 92.5 GHz band based on the MIMO-SAR prototype system. Both simulation and experimental results of different distributed targets show good performance of imaging and do not lose the quality of image reconstruction.

Millimeter-wave systems were typically designed based on a wideband transceiver combined with a 2D spatial scanning scheme in the early stage [6][7][8][9][10]. This results in show data collection and are not suitable for practical applications. For real-time security scenarios, lower system costs and higher data collection efficiency are required. To meet the requirements, a multiple-input-multiple-output (MIMO) array with a synthetic aperture in alone-track is an effective solution to achieve millimeter wave 3D imaging, which can acquire high-resolution images of concealed objects [11][12][13][14][15]. Compared with conventional synthetic aperture radar, MIMO-SAR has various advantages, such as a wide field of view, smaller system volume, and lower cost hardware system. Due to these advantages, MIMO-SAR is being studied and applied in more and more scenarios.
It is easy to understand that the imaging algorithm is definitely one of the key techniques of the MIMO-SAR application. In array imaging applications, the most commonly used focusing algorithm is the backpropagation algorithm (BPA) [16]. Although BPA is plane z = 0 m to illuminate a target in front of the aperture. The coordinates of the transmit and receive antennas are defined as (x t , y l , 0) and (x r , y l , 0) while the location of the scattering point is represented by (x, y, z). To obtain a certain resolution in the z direction for 3D imaging, the transmitted wave should have a certain bandwidth. Hence, we choose a step-frequency continuous wave (SFCW) signal as the transmitted waveform. Under this MIMO-SAR imaging scheme, the corresponding received signal [20] can be denoted as: where o(x, y, z) denotes the target reflectivity function, k i = 2π f i /c(i = 1, 2, . . . , N k ) is the wavenumber, N f is the number of wavenumber domain, respectively, and where R T and R R are the distance between the target and the transmit and receive elements, and R M is the distance from the round trip.
that is parallel to the X-axis scans along the height direction Y, which is located a 0 z = m to illuminate a target in front of the aperture. The coordinates of the and receive antennas are defined as ( ) , , x y z . To obtain a certain resolution in the z d for 3D imaging, the transmitted wave should have a certain bandwidth. Hence, w a step-frequency continuous wave (SFCW) signal as the transmitted waveform. Un MIMO-SAR imaging scheme, the corresponding received signal [20] can be deno , o x y z denotes the target reflectivity function, where T R and R R are the distance between the target and the transmit and rec ments, and M R is the distance from the round trip. Assume that i k is sampled uniformly, such as k is the starting wavenumber, and Δk is the sample interval.
Then, the traditional BP imaging algorithm for MIMO-SAR systems can be e pressed as Assume that k i is sampled uniformly, such as where k 0 is the starting wavenumber, and ∆k is the sample interval. Then, the traditional BP imaging algorithm for MIMO-SAR systems can be easily expressed as For each spatial grid point in the region of interest (ROI), a complex quadruple integration is necessarily performed to obtain the reflectance of the point. Thus it is necessary to optimize the traditional BP algorithm.

The Proposed Algorithm
Through the analysis of the traditional BP algorithm in the previous section, it can be found that the time-consuming operation is a quadruple integration, so the integration step is considered reduced. Meanwhile, to ensure the applicability of the algorithm to arbitrary arrays, a certain transformation is taken to eliminate this integration layer for the wavenumber k i in this section.
Substituting (5) to (1), we have Introduce a Fourier transform pair with respect to k i : (4) can be written as: E(x t , x r , y l , where E(x t , x r , y l , k i ) in Equation (10) is the phase error obtained from R n approximating R M , small ∆R can reduce the phase error. Then, ss(x t , x r , y l , n k ) can be obtained simply by carrying a 1D inverse Fourier transform (FT) for sS(x t , x r , y l , k i ): where int(·) is the round-off operation, and the function δ(n k ) is defined as: In Equation (9), we have the approximation R n ≈ R M . Minimize the sampling interval ∆R to reduce the phase error; that is, the number of inverse FT points should be enough. Assume MN k is the number of inverse FT output, Equation (11) can be rewritten as: ss(x t , x r , y l , n k ) = IFFT k i sS(x t , x r , y l , The solution of n k is the same as that of (12) and ∆R is Sensors 2023, 23, 5979

of 15
Generally, it is believed that the E(x t , x r , y l , k i ) is negligible by satisfying |R m − R n | λ min , that is ∆R λ min , where λ min is the smallest wavelength within the frequency band. Therefore, the number of inverse FT sampling points should satisfy the following: As for a linearized inverse scattering problem, the target reflectivity function can be solved as follows: To improve the precision of the reconstructed image, the phase error term E(x t , x r , y l , k i ) can be compensated as follows: The final proposed algorithm after compensation can be expressed as follows:

The CF-Proposed Algorithm
CF is defined as the ratio of the coherent power P coh to the incoherent power P inc of the target reflectivity function [18], as follows: CF(x, y, z) = P coh (x, y, z) P inc (x, y, z) Combined with Formula (19), P coh and P inc in (20) can be defined as: P inc (x, y, z) = dx t dx r dy l · ss(x t , x r , y l , n k )| n k =int( R M ∆R ) Simply calculate the incoherence factor during image reconstruction to suppress the sidelobe. Then, the final reconstructed image can be expressed as: The proposed algorithm presents a 3D reconstruction scheme by performing IFFT on the wave number to speed up the traditional BP algorithm, then can be easily removed the coherent integrals of the k domain. The error compensation term is also analyzed to reduce the error. Based on the above formulation, the 3D reconstruction scheme can be separated into the following eight steps, as shown in Table 1: Table 1. Detailed procedures of the CF-proposed algorithm for MIMO-SAR imaging.
Input: Recorded the Original Echo for Each Channel.
Step 1. Obtain 4D echo data sS(x t , x r , y l , k i ) based on the MIMO-SAR system.
Step 2. Calculate the condition that the inverse FT sampling points MN f meet according to Formula (16).
Step 3. Perform an inverse MN f FT on the wavenumber dimension of the echo signal in Formula (1).
Step 5. Calculate the incoherent power P inc .
Step 6. Perform coherent accumulation in three dimensions x t , x r , y l .

Computation Complexity
The computational complexity of the proposed algorithm will be analyzed in this section. We use floating-point arithmetic to calculate the computational load required by the proposed algorithm. According to the above imaging formulas, the computation load of the proposed algorithm and traditional BPA can be obtained as follows: where N xt , N xr denote the numbers of transmitters and receivers, N yl represents the number of scanning points in the y-direction. Let N x , N y , and N z denote the voxel in three dimensions, and N x is the number of the voxels in the x-orientation of the 3D focused image.
To further simplify the comparison, we assume that N xt = N xr = √ N and N yl = N x = N y = N z = N k = N, then the computational complexity of (20) and (21) can be expressed more intuitive: As shown in (22) and (23), the computational complexity of the proposed algorithm is much lower than that of the conventional BPA. According to (29), it can be seen that the CF-proposed algorithm has the same quantum of computational complexity.

Error Analysis
The main error is caused by neglecting the error term E(x t , x r , y l , k i ) in Formula (9), and by the compensation operation in Formula (18), the residual maximum can be expressed as follows: Such a phase error can be controlled at a very low level. We can easily control the error to about 0.001π. To ensure imaging quality for each area, the maximum phase error should be lower than 0.25π. Thus the imaging error generated by the proposed algorithm is negligible. Additionally, considering the case of removing the compensation term. The maximum error between the proposed algorithm and the traditional BPA is about 2PHE, this will also not have a significant impact on the reconstruction results. Therefore, in practical applications, we can consider neglecting the effect brought by the error term.

Numerical Simulation and Experimental Verification
In this section, simulation experiments of point targets and letter A scattering models are carried out, and we also design a schematic prototype for experiments to verify the effectiveness of the algorithm. The advantages of the proposed algorithm over the conventional BP algorithm are demonstrated by the comparison of the computational speed, and the superiority of the proposed algorithm over the wave number domain algorithm is verified by setting up sparse array experiments. The specific parameter settings for simulation and experiments are shown in Table 2.

Simulation Results of Point Targets
The 1D MIMO array with 39 transmitters and 6 receivers applied in this simulation is shown in Figure 2. Seven reference target points are set with the center point 0.7 m away from the array in the simulation scenario. All the reconstructed images are divided into 101 × 101 × 100 voxels, and the voxel size is 1 mm in the azimuth dimension. The echo signal of the scatterer is calculated by MATLAB, which uses echo signal model simulation.

Simulation Results of Point Targets
The 1D MIMO array with 39 transmitters and 6 receivers applied in this simula is shown in Figure 2. Seven reference target points are set with the center point 0.7 m aw from the array in the simulation scenario. All the reconstructed images are divided 101 × 101 × 100 voxels, and the voxel size is 1 mm in the azimuth dimension. The e signal of the scatterer is calculated by MATLAB, which uses echo signal model simulat  Figure 3 gives the projection of BPA and the proposed algorithm in XY and XZ pla respectively, and the 3D reconstruction results with dynamic range limited to −10 dB. 1D profile along x-axis and y-axis of the center point is given in Figure 4. The simula results verify that both the proposed algorithm and the BP algorithm perform excel reconstruction of the point target.  Figure 3 gives the projection of BPA and the proposed algorithm in XY and XZ planes, respectively, and the 3D reconstruction results with dynamic range limited to −10 dB. The 1D profile along x-axis and y-axis of the center point is given in Figure 4. The simulation results verify that both the proposed algorithm and the BP algorithm perform excellent reconstruction of the point target. Figure 3 gives the projection of BPA and the proposed algorithm in XY and XZ planes, respectively, and the 3D reconstruction results with dynamic range limited to −10 dB. The 1D profile along x-axis and y-axis of the center point is given in Figure 4. The simulation results verify that both the proposed algorithm and the BP algorithm perform excellent reconstruction of the point target.   (e) (f)

Simulation Results of Letter A Scattering Models
To further evaluate the performance of the proposed algorithm, we set an ideal letter A scattering model for simulation. The echo signal of the MIMO-SAR is calculated using the following formula: where o x i , y j , z m denotes the reflectance function of the target cell x i , y j , z m , which we set to 1. The model shown in Figure 5 is placed at a distance of 0.7 m from the array. The 3D and 2D simulation results are shown in Figure 6, where the 3D dynamic range is limited to −12 dB. Both algorithms achieve a well-reconstructed ideal letter A scattering model. For some artifacts on the periphery, it is a grating-lobes problem due to the sparsity in the y-scan direction, which is also caused by the close proximity of the target. Of course, these grating lobes can be removed by certain methods, the CF method in the proposed algorithm also has some suppression effect on the grating lobes, but this is not considered in this paper.

Simulation Results of Letter A Scattering Models
To further evaluate the performance of the proposed algorithm, we set an ideal letter A scattering model for simulation. The echo signal of the MIMO-SAR is calculated using the following formula: x y z , which we set to 1. The model shown in Figure 5 is placed at a distance of 0.7 m from the array The 3D and 2D simulation results are shown in Figure 6, where the 3D dynamic range is limited to −12 dB. Both algorithms achieve a well-reconstructed ideal letter A scattering model. For some artifacts on the periphery, it is a grating-lobes problem due to the spar sity in the y-scan direction, which is also caused by the close proximity of the target. O course, these grating lobes can be removed by certain methods, the CF method in the pro posed algorithm also has some suppression effect on the grating lobes, but this is not con sidered in this paper.

Lab Experiment 1
The proposed algorithm is further verified experimentally by self-developed millimeter-wave MIMO scanning radar based on a microwave Vector Network Analyzer (VNA), which is shown in Figure 7, while the main parameters of the experiment are listed in Experiment 1 of Table 2. Due to the width limitation of the antenna front-end, the interval between the adjacent transmitter and receiver in the experiment is 2 cm. The total length of the array is 0.335 m. The fourth motor is used to control the horizontal guide rail to scan in the vertical direction to achieve a synthetic aperture. After completing a round of horizontal transmission and reception, the horizontal guide railway moves 5mm along the vertical y-direc-

Lab Experiment 1
The proposed algorithm is further verified experimentally by self-developed millimeterwave MIMO scanning radar based on a microwave Vector Network Analyzer (VNA), which is shown in Figure 7, while the main parameters of the experiment are listed in Experiment 1 of Table 2.

Lab Experiment 1
The proposed algorithm is further verified experimentally by self-developed millimeter-wave MIMO scanning radar based on a microwave Vector Network Analyzer (VNA), which is shown in Figure 7, while the main parameters of the experiment are listed in Experiment 1 of Table 2. Due to the width limitation of the antenna front-end, the interval between the adjacent transmitter and receiver in the experiment is 2 cm. The total length of the array is 0.335 m. The fourth motor is used to control the horizontal guide rail to scan in the vertical direction to achieve a synthetic aperture. After completing a round of horizontal transmission and reception, the horizontal guide railway moves 5mm along the vertical y-direc- Due to the width limitation of the antenna front-end, the interval between the adjacent transmitter and receiver in the experiment is 2 cm. The total length of the array is 0.335 m. The fourth motor is used to control the horizontal guide rail to scan in the vertical direction to achieve a synthetic aperture. After completing a round of horizontal transmission and reception, the horizontal guide railway moves 5 mm along the vertical y-direction, and the next round of data acquisition starts.
The reconstructions for the three models in Figure 8 are compared in Figure 9a-d. The reconstruction results of all the targets are shown in the dynamic range of −20 dB, which is mainly determined by the number of array transceiver units, location, and signal-to-noise ratio of the measurement.

Lab Sparse Experiment 2
To verify the advantages of the proposed algorithm over traditional wavenumber domain algorithms, comparative tests are conducted with the three targets in Figure 8.

Lab Sparse Experiment 2
To verify the advantages of the proposed algorithm over traditional wavenumber domain algorithms, comparative tests are conducted with the three targets in Figure 8.

Lab Sparse Experiment 2
To verify the advantages of the proposed algorithm over traditional wavenumber domain algorithms, comparative tests are conducted with the three targets in Figure 8. PSM [17], as a classical wavenumber domain algorithm, is used in this paper to compare with the proposed algorithm. The number of transmitting array elements is reduced from 39 to 20 by equally spaced sparsity, while the number of receiving array elements remains unchanged, as shown in Figure 10. The parameters of the experiment are consistent with Experiment 2 in Table 2.
Comparative test results are shown in Figure 11. For uniformly sparse arrays, the proposed algorithm still exhibits good image reconstruction performance, as shown in Figure 11c; compared to Figure 9c, the side lobe of the reconstructed target image has significant side lobes, which can be easily removed by certain measures. For the PSM algorithm, the presparing reconstructed target image is shown in Figure 9b, while the reconstructed target image after sparse is shown in Figure 11b. The details of the target have been severely distorted after sparsity, leading to a serious loss of detail in the reconstructed image. This is also the biggest limitation of wavenumber domain algorithms. Comparative tests indicate that the proposed algorithm has better adaptability to array element sparsity than the wavenumber domain algorithm. This can reduce the hardware cost of the system to a certain extent. Figure 11d shows the 3D reconstruction results of the CF-proposed. Compared with Figure 11c, the sidelobe suppression effect is obvious. PSM [17], as a classical wavenumber domain algorithm, is used in this paper to compare with the proposed algorithm. The number of transmitting array elements is reduced from 39 to 20 by equally spaced sparsity, while the number of receiving array elements remains unchanged, as shown in Figure 10. The parameters of the experiment are consistent with Experiment 2 in Table 2. Comparative test results are shown in Figure 11. For uniformly sparse arrays, the proposed algorithm still exhibits good image reconstruction performance, as shown in Figure 11c; compared to Figure 9c, the side lobe of the reconstructed target image has significant side lobes, which can be easily removed by certain measures. For the PSM algorithm, the pre-sparing reconstructed target image is shown in Figure 9b, while the reconstructed target image after sparse is shown in Figure 11b. The details of the target have been severely distorted after sparsity, leading to a serious loss of detail in the reconstructed image. This is also the biggest limitation of wavenumber domain algorithms. Comparative tests indicate that the proposed algorithm has better adaptability to array element sparsity than the wavenumber domain algorithm. This can reduce the hardware cost of the system to a certain extent. Figure 11d shows the 3D reconstruction results of the CF-proposed. Compared with Figure 11c, the sidelobe suppression effect is obvious PSM [17], as a classical wavenumber domain algorithm, is used in this paper to compare with the proposed algorithm. The number of transmitting array elements is reduced from 39 to 20 by equally spaced sparsity, while the number of receiving array elements remains unchanged, as shown in Figure 10. The parameters of the experiment are consistent with Experiment 2 in Table 2. Comparative test results are shown in Figure 11. For uniformly sparse arrays, the proposed algorithm still exhibits good image reconstruction performance, as shown in Figure 11c; compared to Figure 9c, the side lobe of the reconstructed target image has significant side lobes, which can be easily removed by certain measures. For the PSM algorithm, the pre-sparing reconstructed target image is shown in Figure 9b, while the reconstructed target image after sparse is shown in Figure 11b. The details of the target have been severely distorted after sparsity, leading to a serious loss of detail in the reconstructed image. This is also the biggest limitation of wavenumber domain algorithms. Comparative tests indicate that the proposed algorithm has better adaptability to array element sparsity than the wavenumber domain algorithm. This can reduce the hardware cost of the system to a certain extent. Figure 11d shows the 3D reconstruction results of the CF-proposed. Compared with Figure 11c, the sidelobe suppression effect is obvious

Discussion
To compare the performance of the algorithms more clearly, the structural similarity (SSIM) evaluation index is introduced to evaluate the focusing performance of 2D images [20][21][22]. The calculation of SSIM can be expressed by the following formula:

Discussion
To compare the performance of the algorithms more clearly, the structural similarity (SSIM) evaluation index is introduced to evaluate the focusing performance of 2D images [20][21][22]. The calculation of SSIM can be expressed by the following formula: where I 1 represents the reference image and I 2 represents the image to be evaluated. µ 1 , µ 2 , σ 1 , σ 2 and σ 12 are the average, variance, and cross-covariance of the two images, respectively. C 1 and C 2 are constants that maintain stability in order to keep the denominator from being 0, usually take C 1 = (K 1 L) 2 , SSIM is calculated for values from 0 to 1, and the larger its value, the smaller the gap between the image to be evaluated and the reference image.
Using the BPA as a reference, the difference between the proposed algorithm, PSM, and BPA are compared by calculating SSIM. Then, we evaluate the effectiveness of the algorithm using computational time and SSIM for the simulations and experiments in Sections 3 and 4; as shown in Table 3, it can be seen that the proposed algorithm has obvious efficiency advantages. Table 3 also shows the computation time of the CF-proposed algorithm in two sets of experiments to highlight the impact of CF on the computational efficiency of the proposed algorithm.

Conclusions
In this paper, a modified BPA is proposed for 3D reconstruction based on MIMO-SAR. To address the problem of the high redundancy of the traditional BPA and its iteration operations, the wavenumber dimension data are improved, and the inverse FT method is adopted to accelerate the algorithm without damaging the image quality. The Experiment 1 results in Table 3 show that, compared with BPA, the computational efficiency of the proposed algorithm is significantly improved, with an SSIM value of 0.9960 and SSIM value of 0.8570 for PSM, indicating that the proposed algorithm performs 3D reconstruction without compromising the imaging quality. The results of Experiment 2 show that the effect of array sparsity likewise does not degrade the quality of the reconstructed images by the proposed algorithm, while the quality of the images reconstructed by PSM is severely degraded. Experimentally, it is proved that although the wave number domain algorithm has obvious advantages in reconstruction speed, it loses imaging accuracy to a certain extent. In an overly sparse array, the wavenumber domain algorithm reconstruction effect is severely degraded, and the original information of the target cannot be recovered; the proposed algorithm does not have this problem with BPA.
In addition, the CF-proposed algorithm is used to reconstruct the images under each of the two sets of experiments to achieve the suppression of the sidelobe energy without losing the image quality. The main contribution of this paper is as follows: The first is to propose an algorithm that has the same reconstructed image quality as the traditional BPA but is much more efficient than the BPA. The second is to combine CF with the proposed algorithm to suppress sidelobe energy.
Author Contributions: B.L. performed the theoretical study, conducted the experiments, processed the data, and wrote the manuscript. C.L. designed the imaging system and revised the manuscript together with Y.J. helped in performing the experiments. X.L. and G.F. provided the experiment equipment and funds for the research. All authors have read and agreed to the published version of the manuscript.