MIMO Radar Accurate 3-D Imaging and Motion Parameter Estimation for Target with Complex Motions

In this paper, three-dimensional (3-D) multiple-input multiple-output (MIMO) radar accurate localization and imaging method with motion parameter estimation is proposed for targets with complex motions. To characterize the target accurately, a multi-dimensional signal model is established including the parameters on target 3-D position, translation velocity, and rotating angular velocity. For simplicity, the signal model is transformed into three-joint two-dimensional (2-D) parametric models by analyzing the motion characteristics. Then a gridless method based on atomic norm optimization is proposed to improve precision and simultaneously avoid basis mismatch in traditional compressive sensing (CS) techniques. Once the covariance matrix is obtained by solving the corresponding semi-definite program (SDP), estimating signal parameters via rotational invariance techniques (ESPRIT) can be used to estimate the positions, then motion parameters can be obtained by Least Square (LS) method, accordingly. Afterwards, pairing correction is carried out to remove registration errors by setting judgment conditions according to resolution performance analysis, to improve the accuracy. In this way, high-precision imaging can be realized without a spectral search process, and any slight changes of target posture can be detected accurately. Simulation results show that proposed method can realize accurate localization and imaging with motion parameter estimated efficiently.


Introduction
Owing to accurate localization and imaging performance, radar is widely applied in many imaging fields. Unfortunately, the localization errors will increase so the image will be distorted and worsened when target complex motions are taken into account, particularly for slow time-varying motions containing translations and rotations [1,2].
Under this circumstance, relevant studies mostly concentrate on synthetic aperture radar (SAR), inverse synthetic aperture radar (ISAR), and 3-D interference inverse synthetic aperture radar (3-D InISAR) [3,4]. Nevertheless, the imaging accuracy and real-time performance will be worse if the synthetic aperture time cannot be optimized reasonably. Following this, some relevant improvements have been applied, containing optimization of imaging accuracy [5][6][7][8][9], selection of optimal imaging time [10], improvements of imaging efficiency by CS and sparse sampling techniques [11][12][13]. However, there are still inevitable defects. On one hand, the platform is usually required to ensure baselines unchanged during synthetic aperture time, which is impractical in actual. On the other hand, the geometric models are usually one or two-dimensional which are insufficient to accurately ⊗ denotes Kronecker product. (·) T , (·) * , (·) H , (·) −1 and (·) + denote the transpose, conjugate, conjugate transpose, inverse, and pseudo-inverse operations, respectively.

Multi-Dimensional Echo Model
Appropriate array structure is necessary for MIMO radar 3-D imaging and parameter estimation [17], the layout of array in the paper is shown in Figure 1. The paper chooses a uniform linear array consists M transmitters along X axis direction, T m represents the m-th transmitter where m = 0, 1, · · · , M − 1. Signal s m (t) = p m (t) exp(j2π f c t + ϕ m ) is Hadamard orthogonal-phase encoded signal of T m , where p m (t) and f c are envelope and carrier frequency, signal orthogonality is ensured by adjusting phase ϕ m . Then a N × L uniform planar array is designed as the receiving array. Axial directions of Y and Z are considered to be the array line directions and R nl is used to represent the receiver in n-th row, l-th column where n = 0, 1, · · · , N − 1 and l = 0, 1, · · · , L − 1. In this paper, a large target with translations and rotations is considered in the model, thus all the scattering points have same motion states. Moreover, they all obey Swerling II distribution and scattering coefficients remain unchanged in one pulse period. Following the array model, the echo signal at R nl receiver can be expressed as where σ k is scattering coefficient of the k-th scattering point, τ k m,n,l = (T m k + R nl k) /c is delay time and c is propagation speed of electromagnetic signal. T m k and R nl k represent the distance from T m to the k-th scattering point and the distance from k-th scattering point to R nl , respectively. f d = 2V d /λ is Doppler frequency caused by target motions, λ is signal wavelength and V d is the synthesis of translation velocity and rotational angular velocity. After removing carrier, the receiving signal at R nl from T m can be written as Due to V d c, the model in (2) can be simplified as Taking P as the reference center of the region, then the echo d P m,n,l (t) reflected from P can be used as reference signal to compensate the target echo signal D m,n,l (t) = d m,n,l (t) · d P m,n,l (t) * Then according to the geometrical model shown in Fig.1, taking ∆R to represent the range deviation term in (4) and it can be finally written as following, the proof is shown in Appendix A.1.
where ( −→ T m P) and ( −→ R nl P) respectively represents the unit direction vector from T m and R nl to center P. Then, in order to intuitively describe the echo, we set (P X , P Y , P Z ) and (P X + x k , P Y + y k , P Z + z k ) as the coordinates of P and point k, d X is internal spacing of transmitting array, d Y and d Z respectively denote the row and column spacing inside the receiving array. Considering the coordinates of T 0 transmitter is (r X , 0, 0), the coordinates of R 00 receiver is (0, r Y , r Z ), thus the coordinates of T m and R nl are (r X + md X , 0, 0) and (0, r Y + nd Y , r Z + ld Z ), respectively. Following this, we can get where R 0 is the reference distance from P to coordinate center O. A pulse transmitted by radar is divided into Q samplings, t q = q · T p /Q(q = 0, 1, . . . , Q − 1) is the sampling time and T p is pulse width. Thus, the echo model can be written as where However, it is difficult to directly extract parameters from the model in (6) due to the large dimension. Moreover, V d is also hard to be dealt because it is the synthesis of translation velocity and 3-D rotation velocity. Therefore, we transform this model into 2-D parametric models for simplicity.

Joint 2-D Parameter Models
Due to f d t q = 2V d t q /λ = 2R q /λ, the range term produced by target motions can be expressed as following, with translation velocity V and 3-D rotating angular velocities ω x , ω y , ω z are considered where ∆R ω is caused by target rotations and its projections are ∆x, ∆y and ∆z, respectively. ∆R X , ∆R Y and ∆R Z are the projections of ∆R in three dimensions. Obviously, these deviation terms are only related to their own dimensions and do not affect each other.
It is noted that the velocity parameters are considered invariable due to the short processing interval of MIMO radar, i.e., V, ω x , ω y and ω z are all constant during Q samplings in one pulse. As for the roll, pitch, and yaw rotations of target shown in Figure 1, the following rotation matrices are supported by basic navigation theory where θ r (t), θ p (t) and θ y (t) represent the time-varying angles caused by roll, pitch, and yaw rotations. Based on this, ∆x, ∆y and ∆z in Equation (8) can be final expressed as following, the proof is presented Finally, the model in (6) can be transformed into following joint three 2-D parametric models based on Equations (8) and (10), with partial complex variables represented by new parameters Hence, we can get the parameter models of X, Y, and Z in low-dimensional space through above decoupling process of Doppler shift term. Compared with the impossibility of estimating motion parameters directly from the Doppler frequency in (6), it becomes feasible to obtain the estimation results of target location and motion parameters from the models in (11). Simultaneously, the problem of large complexity caused by large dimension can also be solved. Then according to the consistency of three models in (11), the processing in X direction will be taken as an example.

2-D Parameters Estimation without Basis Mismatch
For radar imaging, traditional sparse recovery method such as Orthogonal Matching Pursuit (OMP) will result in basis mismatch due to the construction of sparse dictionary, which depends on the discretization process of continuous variable. In view of this, it is a feasible method to avoid this phenomenon by taking SVD decomposition of echo covariance matrix and then extracting the eigenvalues from the signal subspace accordingly. However, it should be noted that it is impractical to extract all K-column eigenvectors to construct the signal subspace because the echo in (11) is one-dimensional. Therefore, to solve this problem and avoid basis mismatch, a gridless method based on atomic norm optimization is proposed in this section.
As a penalty function for convex optimization problem, atomic norm shows its convenience in solving underdetermined linear inverse problems [27][28][29][30]. For X dimension echo model in (11), we define atoms where . So, the model can be written as Then the atomic set is defined as x max ] and [α min , α max ] represent the range of x and α. Obviously, the basic components in A construct the full echo D x . In this model, A and σ are considered continuous. Therefore, the atomic norm of the echo can be expressed as where σ k = |σ k | e jφ k , σ k and φ k are amplitude and initial phase. However, it is impractical to directly construct the covariance matrix of the echo because it is not Toeplitz. Considering this, the dimension of echo is extended: , so that σ can be replaced by an equivalent diagonal matrix. Thus, we can get the rank-K matrix where A(x, α) = [a(x min , α min ), · · · , a(x max , α max )] and it is obviously a Vandermonde matrix with infinite columns. Then, with Gaussian noise considered in D x , the covariance matrix can be construct based on the optimization problem where µ is denoised truth echo and ρ is the regularization coefficient. Accordingly, this atomic norm minimization problem can be transformed into an approximate semi-definite program (SDP) as |σ k | 2 , T (s) denotes the Toeplitz covariance matrix and its first column is s. Furthermore, it has been proved that the noise can be efficiently suppressed by the optimization condition, i.e., the obtained s mainly includes the scattered sampling data even at low SNR. According to the Toeplitz covariance matrix T (s), the estimation of parameters can be realized based on ESPRIT algorithm. K large singular values can be found and relevant signal subspace U s can be obtained by taking singular value decomposition of T (s), then two subspace of U s can be obtained as where Then the parameters x and α can be directly obtained according to (13) with the K eigenvalues extracted from Ψ by eigenvalue decomposition. Following the results, scattering coefficients distributions of the points in X-α plane can be expressed as where A 0 ∈ C MQ×K is constructed by selecting corresponding columns in A(x, α), according to the parameters extraction results x 1 x 2 · · · x K and α 1 α 2 · · · α K from (20). Following this, all estimation results of three models in (11) can be obtained without any spectral search process

Target 3-D Imaging with Motion Parameters Estimated
According to (22), F X , F Y and F Z only represent the distributions of scattering points in X-α, Y-β and Z-γ planes, but the orders of the points are different in these matrices. For example, the k-th point may locates on the a-th row in F X , but the b-th row in F Y and the c-th row in F Z , where a, b, c are different because it maybe not one-to-one with the rows among F X , F Y and F Z . So the true 3-D imaging result cannot be obtained unless the coordinates in F X , F Y and F Z are paired accurately. Based on this, a set is constructed by combining all of X, Y and Z coordinates in (22), which surely contains the coordinates of K true scattering points. Γ = [g 0,0,0 , · · · , g ς,ε,ζ , · · · , g K,K,K ] g ς,ε,ζ = a ς x ⊗ a ε y ⊗ a ζ z , ς = 1, 2, · · · , K, ε = 1, 2, · · · , K, ζ = 1, 2, · · · , K} where Then a set is defined as χ = [σ 0,0,0 , · · · , σ ς,ε,ζ , · · · , σ K,K,K ]. As for the elements in σ ς,ε,ζ , σ ς 0 ,ε 0 ,ζ 0 = (σ x ς +σ y ε +σ z ζ ) 3 only if ς 0 = ς, ε 0 = ε and ζ 0 = ζ, otherwise σ ς 0 ,ε 0 ,ζ 0 = 0. Then the true positions can be obtained by following optimization problem where D = [D 1,1,1 (t 1 ), · · · , D M,N,L (t 1 )] T . All the parameters are all known in this problem, so it is actually a process of finding K minimum values from a finite set consists of K 3 elements. Therefore, 3-D imaging result is obtained T = [XŶẐσ ] (26) where T ∈ C K×4 , each row of T denotes the 3-D positions and scattering coefficients of every scattering point.
can be regarded as position selection matrices according to the results of (25). Nevertheless, the obtained results are not the most accurate because the effects of Doppler frequency shift are ignored. For this paper, accurate motion estimation plays a significant role in position compensation and target identification. Thus, efficient estimation of the motion parameters is carried out according to the relationships among the parameters in (11) and (26), we can get whereα = J 1 α,β = J 2 β,γ = J 3 γ, α, β, γ are the second column in F X , F Y , F Z , and Therefore, Least Square (LS) method can be directly used to solve the problem in (27), then motion parameter vector can be estimated

Pairing Correction
The estimation precision of the motion parameters is mostly determined by (27), so it must be guaranteed that the results of (22) and (26) are correct. However, the existence of model errors will lead to the failure of parameter estimation. After pairing, the coordinates of two points in one resolution cell may be disorder. For example, (x a , y a , z a ) and (x b , y b , z b ) are coordinates of point a and point b in same X resolution cell, with y a , z a far away from y b , z b and x a close to x b . In this way, it is hard to distinguish them in the pairing process due to the strong coherence between them. So the final pairing result maybe (x b , y a , z a ) or (x a , y b , z b ), which will break the structure of J 1 . As a result, the estimation of ω y and ω z in (29) will be inaccurate due toα = J 1 α in (27). Similarly, when y or z coordinates are difficult to distinguish, the estimation of motion parameters will also be affected. Therefore, pairing correction is carried out in this section to remove this registration error and improve the accuracy of the method, particularly the estimation accuracy of motion parameters.
First, the imaging resolution performance of the model is studied by analyzing the point spread function theoretically. Here we define the point spread function of the X model as where κ is a normalized parameter to ensure the maximum value of Ps f (k, k 0 ) is 1. Following this, we take the mathematical simplification where T p = Q × t q is pulse width. Then the limit resolution in X-α plane is that is, when the distance between two points in X direction is less than ρ x , the main lobes of them will overlap and makes their X coordinates difficult to be distinguished. As a result, the estimation of motion parameters will be inaccurate.
In a similar way, the resolutions of Y-β plane and Z-γ plane can also be obtained In the following, a rough estimation is developed by judging and removing the points in the same cell. First, we need to determine whether there are points difficult to be distinguished. Taking X dimension as example, for any point a and point b, a judging condition is set as if it does not meet the condition, there must be more than one point in one cell, then these points are taken out from the pairing result. Therefore, it is actually a process of constantly eliminating indistinguishable target points by pairing judgment. The convergence condition of this process is that all the remaining target points satisfy the judgment condition in (34). Then, the parametric coarse estimation is developed with the target points satisfying (34) where ϕ ∈ C 1×K , ω ∈ C 1×4 , Φ ∈ C 4×K , K denotes the number of the points satisfying (34). With the process of pairing judgment and correction in (34) and (35), the coordinates between these points in same cell can be distinguished according to the coarse estimation of ω . Then the position selection matrix J 1 can be corrected and Equations (26) and (27) are updated. Accordingly, Equation (29) will be carried out based on the corrected results. In this way, the registration error can be removed, and the high accuracy of localization and motion parameter estimation can be guaranteed.
As a summary, the flow chart of the whole method is shown in Figure 2. It is noted that the main computational load of the algorithm includes three parts. The first one is caused by the SDPT3 or ADMM [29] when the atomic norm optimization problem is solved. The second one is caused by singular value decomposition and eigenvalue decomposition when ESPRIT algorithm is used to calculate the target positions, which is O(2( The third one is brought by the process of pairing and motion parameter estimation, which is O(MNLK 3 ). Thus, the computational complexity of the algorithm is mainly affected by the number of antennas, samplings, and targets.
The method is applicable for many scenarios such as ships, airplanes, and accurate imaging can be realized with stable application environment. However, the performance will be deteriorated once the application environment becomes worse, such as the sea surface full of clutters and sea waves, which will be further studied in future research.

Simulation Results and Discussion
In this section, relevant simulation results are shown to verify that proposed method can realize accurate imaging and motion estimation for target with complex motions. In the simulations, a ship target with radial translation and 3-D rotations is taken into account. Scattering points are set on the ship hull shown in Figure 3 and scattering coefficients are all set to 1. Following this, imaging results and motion estimation results are shown in Figures 4-9. First, in order to verify the feasibility of proposed method, ship target are imaged at two moments t A and t B with different motion parameters, meanwhile radar parameters are set as Table 1 and motion parameters of t A and t B are set as Table 2.

Values of Motion Values of Motion Parameters at t A Parameters at t B
Radial translation velocity (m/s) 4.     6 show their 2-D projections. In these figures, accurate localization and imaging can be intuitively presented. The deviation between true scattering points and the estimation results are small and some coincidence points appear in the figures when the estimated results of these points are more close to the real situation, which proves the feasibility and high accuracy of proposed method. Moreover, Table 3 shows the motion estimation results at t A and t B . High estimation accuracy of proposed method can be verified by comparing Table 3 with Table 2, which indicates that any slight changes of target posture can be efficiently detected with the proposed imaging method.  Tables 4-6 show the time performance with different M, Q, and L which represent the number of transmitters, samplings, and targets, respectively. For simplicity, the effects caused by receivers are ignored because they are similar to the transmitters according to the computation analysis, so the number N, L of receivers are set as constant as shown in Table 1. Then, the rest of the parameters follow the settings in Table 1. As is shown in the tables, with the increase of M, Q, and L, the running time of the algorithm increases obviously. Therefore, the time efficiency decreases accordingly, which coincides with the computation analysis of the method.   Figure 7 shows the relationship between the imaging error performance and SNR. In simulation Figure 7a, we take MIMO-ISAR method [21] and modified OMP method [13] for comparison. Apparently, our method has minimum imaging errors in the figure, which verifies its high accuracy over other imaging methods. In simulation Figure 7b, the SACR-iMAP method [25] and S-TLS method [26] are simulated for references. It is obvious that these methods cannot further improve the imaging accuracy because they only focus on the improvement of 2-D resolution and the basis mismatch errors cannot be completely removed. In contrast, proposed method has higher precision. On one hand, it benefits from the denoising performance of SDP. On the other hand, target positions can be directly calculated without any spectral search process, which can avoid basis mismatch.  Figure 8 shows the comparison results of motion estimation accuracy, with MIMO-ISAR method and modified OMP method for comparison. It can be seen that proposed method has more precise motion estimation performance. In fact, the high precision is guaranteed by the gridless method and the pairing correction, where the former avoids the basis mismatch and the latter removes registration errors.
In addition, to further show the performance of proposed method under basis mismatch circumstance, Figure 9 presents the error curve of localization and motion estimation in more detail. Take 500 Monte Carlo experiments and 6 random scattering points with random motion parameters, then SNR step size is set to 1dB from −10 dB to 14 dB in the simulations. As a result, the high precision of the method can be efficiently illustrated from the simulation results Figure 9a,b.

Conclusions
In this paper, we have presented a 3-D MIMO radar imaging method with motion parameter estimation for target with complex motions. The method can reduce process difficulty by building joint 2-D parameter models. Then efficient imaging and accurate motion parameter estimation can be guaranteed by the gridless method and pairing correction process, which can eliminate basis mismatch and remove registration errors. Simulation results show that proposed method is more suitable for accurate localization and imaging of the target with complex motions.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
Appendix A.1. Proof of (5) As is shown, ∆R represents the range deviation term in (4), so it can be further processed as In the formula, the modulus values of the vectors −→ T m k, −→ R nl k, −→ T m P and −→ R nl P are used to represent the corresponding range terms, then we can get following approximations where ( −→ T m P) and ( −→ R nl P) respectively represents the unit direction vector from T m and R nl to center P, then Equation (A1) can be expressed as