Multi-Target Localization of MIMO Radar with Widely Separated Antennas on Moving Platforms Based on Expectation Maximization Algorithm

: This paper focuses on multi-target parameter estimation of multiple-input multiple-output (MIMO) radar with widely separated antennas on moving platforms. Aiming at the superimposed signals caused by multi-targets, the well-known expectation maximization (EM) is used in this paper. Target’s radar cross-section (RCS) spatial variations, different path losses and spatially-non-white noise appear because of the widely separated antennas. These variables are collectively referred to as signal-to-noise ratio (SNR) ﬂuctuations. To estimate the echo delay/Doppler shift and SNR, the Q function of EM algorithm is extended. In addition, to reduce the computational complexity of EM algorithm, the gradient descent is used in M-step of EM algorithm. The modiﬁed EM algorithm is called generalized adaptive EM (GAEM) algorithm. Then, a weighted iterative least squares (WILS) algorithm is used to jointly estimate the target positions and velocities based on the results of GAEM algorithm. This paper also derives the Cramér-Rao bound (CRB) in such a non-ideal environment. Finally, extensive numerical simulations are carried out to validate the effectiveness of the proposed algorithm.


Introduction
Recent years have witnessed a rapid development in signal processing technologies for array or multi-node radar systems, such as space-time adaptive processing radar [1,2], distributed coherent radar [3], multiple-input multiple-output (MIMO) radar with colocated antennas [4] and MIMO radar with widely separated antennas [5][6][7][8][9]. Among them, MIMO radar with widely separated antennas simultaneously observes the target from multiple different angles, thus reducing the probability that all antennas measure small radar cross-section (RCS) and low Doppler shift [10]. Furthermore, putting these antennas on moving platforms can ensure the radar system gets superior sensing opportunities. For example, in an urban area where targets may be obscured by buildings, the moving radars can arrive in more suitable positions to obtain the target information [11]. Hence, this paper concerns multi-target localization of MIMO radar with widely separated antennas on moving platforms.
In the field of target localization based on MIMO radar with widely separated antennas, Ref. [7] analyzes the possibility of high-precision target localization through the derivation of Cramér-Rao bound (CRB) and ambiguity function. Ref. [10] studies the optimized system/configuration design based on CRB. Ref. [12] uses the best linear unbiased estimator to analyze the geometric dilution of precision (GDOP) of target location for MIMO radar non-ideal environment is derived. Finally, numerical simulations are conducted to verify the proposed algorithm.
The rest of the paper is organized as follows: Section 2 firstly establishes the signal model which exists the superimposed signal and SNR fluctuations. Then, the GAEM algorithm is proposed, and CRB is derived. Section 3 verifies the feasibility of the algorithm through numerical simulations. Sections 4 and 5 draw the discussion and conclusion, respectively. Figure 1 shows an application scenario of MIMO radar with widely separated antennas on moving platforms. The radar system contains multiple radar nodes on moving platforms, and the distance between radar nodes is much larger than half wavelength. As shown in Figure 1, multiple radar nodes are mounted on the unmanned aerial vehicle (UAV) and shuttle between buildings to approach the target, thereby completing the estimation of target parameters. However, due to the large distance between radar nodes, the propagation paths of the echoes received by different radar nodes are different. When the building blocks the echo signal of a certain propagation path, the echo SNR of this propagation path will be much lower than the echo SNR of other propagation path. As a result, the reliability of the target information carried by the echo in this propagation path decreases. At the same time, the targets (such as cars) in the urban scene are densely distributed. The generated superimposed signals have an impact on the target positioning. Topology diagram of MIMO radar and targets is shown in Figure 2, which shows two practical problems concerned in this paper: (1) target echo overlap, (2) different SNR caused by different radar receiver noises and different propagation paths. Let K, L, and Q 0 denote the total number of transmitting radar nodes, receiving radar nodes, and targets, respectively. Further, τ q lk and f q lk represent the delay and Doppler shift transmitted by the kth radar node, reflected by the qth target and received by the lth radar node, respectively. Suppose that the coordinates of the kth transmitter and lth receiver are p t k = x t k , y t k T and p r l = x r l , y r l T , respectively, and the velocities of the kth transmitter and lth receiver are

Signal Model
, v y r l ] T , respectively. we denote superscript T as the transpose.
denote the location and velocity of the qth target, respectively. Then, the echo delay is related only to the distance between the radar node and the target, i.e., where c denotes the speed of light in free space; · 2 represents the l 2 norm. The echo Doppler shift is related only to the projection of the velocity (radar velocity and target velocity) in the radar-target direction, i.e., where λ denotes the wavelength of the carrier frequency. Let x r l ,ỹ r l T , ∆x r l , ∆y r l T , [ṽ x r l ,ṽ y r l ] T and [∆v x r l , ∆v y r l ] T be the nominal position, position deviation, nominal velocity, and velocity deviation of the lth receiving radar, respectively. Then, we define x r l =x r l + ∆x r l , y r l =ỹ r l + ∆y r l , v x r l =ṽ x r l + ∆v x r l , v y r l =ṽ y r l + ∆v y r l . (3) The system deviation definition of the kth transmitting radar node is similar to that in (3). Let the number of slow times be M. Let f s denote the sampling frequency of radar system. Then, the received waveform at the lth receiver, time n f s and mth pulse is where E is the total energy transmitted by MIMO radar; ξ q lk represents the target reflection coefficient multiplied by path loss; s k [n] represents the waveform emitted by the kth transmitter; PRI is the pulse repetition interval; φ q lk is the phase deviation; w lm [n] is the complex white Gaussian noise with zero mean and variance of (σ wl ) 2 . Suppose s k [n] is a unit energy waveform that satisfies ∑ N−1 n=0 |s k [n]| 2 f s = 1, where N = PRI · f s . In order to simplify the solution process of the gradient descent algorithm, the echo signal can be transformed to the frequency domain. This means that the envelope movement among the radar nodes is converted to a phase shift in the frequency domain. The Fourier transform of (4) is Based on the definition of w lm [n], ς lm [u] is a Gaussian noise with a mean of 0 and a variance of (σ l ) 2 . Decomposing Y lm [u] into its signal components: where ς q lm [u] is the Gaussian noise with zero mean and variance of σ For convenience, an alternative parameter vector is defined as where The subvectors of γ in (9) are Every radar node is assumed to be a transceiver, which means L equals K. Collect the parameter of interest in where β q is a subvector composed of the qth target position and velocity; other subvectors in (11) represent radar position deviations and radar velocity deviations.
In this study, echo delays, Doppler shifts, η, and σ 2 (represented by vector ψ) are estimated first. Then based on these estimates, the target locations, velocities and radar system deviations (represented by vector α) are estimated.

Stage 1: Delay-Doppler-SNR Estimation
In this section, a fast EM algorithm is proposed to estimate ψ. The proposed GAEM algorithm is illustrated in Figure 3.

Q function and Derivation of Complete Data
The , which denotes the complete data selected by the EM algorithm in this study, and we obtain According to the definition of the Q function in the EM algorithm, the Q function of the (i + 1)th iteration is where E X [·] represents the expectation of the formula with respect to the random vector X. (15) denotes the estimated result of the ith iteration, which is known in the (i + 1)th iteration. The parameters to be estimated in the Q function here are echo delay/Doppler shift γ, echo amplitude η and noise variance σ 2 . Substituting (14) into (15), we can get where (·) * in (16) and (18) (18) and (19) satisfies σ (16), the Q function can be reformulated as Unlike the use of the EM algorithm to estimate the parameters of the superimposed targets in [15], this study extends the parameters to be estimated of the Q function in the EM algorithm, where the echo amplitude η and noise variance σ 2 are estimated in each iteration of EM algorithm. The estimations of echo amplitude η and noise variance σ 2 conform to the actual scenario of SNR fluctuations in different propagation paths of the MIMO radar, which ensures the accuracy of multi-target parameter estimation results.

Estimation of Parameters in SNR
Since the signals transmitted from different radar nodes are assumed to be orthogonal, we obtain Based on (20) and (21), the maximum likelihood (ML) estimator of η q lk is given aŝ By taking the partial differentiation of the formula in arg max{·} with respect to η q lk and setting the derivative equal to zero, the explicit expression of the ML estimateη q lk can be written asη Equation (20) can be simplified based on (21). Next, by substituting (23) into simplified (20), it can be seen from (8) that we can change ψ in (20) into γ; σ 2 , thus obtaining a new Q function, which is named as Q 1 in this paper, as shown below: Since the derivation of (24) involves the estimation of η q lk , and the parameter estimation accuracy is affected by the SNR, (24) is an approximation of (20). By taking the partial derivative of (24) with respect to σ q2 l and setting the derivative equal to zero, the explicit expression for the ML estimateσ q2 l is obtained aŝ Now, the explicit expression for the ML estimateη q lk andσ q2 l are derived.

Echo Delay and Doppler Shift Estimation
Substituting (25) into (24), we simplify γ; σ 2 in (24) to γ, and Q 1 is updated and named as Q 2 : Similar to the derivation of Q 1 , compared with (24), the derivation process of (26) includes the estimation of σ q2 l . The parameter estimation accuracy will be affected by the SNR, so (26) is an approximation of (24). Since σ (26) are composed of the estimation results of the ith iteration or echo data, the echo delay and Doppler parameter vector γ can be estimated by the following optimization criterion: For one kqlth propagation path, the dimension of γ can be reduced to 2. Then, we can draw an objective function of (27) as below: From Figure 4, we can find the objective function is convex when the distance varies from 3970 m to 4030 m and the speed varies from 0 m/s to 10 m/s. Hence, different from the traditional grid point search method, this study employs the gradient descent algorithm to solve the echo delay and Doppler shift. That is, instead of computing the maximizer of Q 2 ·, ψ (i) , we can find a vector γ (i+1) that satisfies Q 2 γ (i+1) , ψ (i) ≥ Q 2 γ (i) , ψ (i) (i.e., an improvement). This is conducive for reducing the computational complexity of each M-step in EM algorithm. To achieve Q 2 γ (i+1) , ψ (i) ≥ Q 2 γ (i) , ψ (i) , one iterative gradient algorithm can be formulated as where µ is the step size that depends on ∇ γ Q 2 γ, ψ (i) γ=γ (i) . So far, the delay, Doppler and parameters related to SNR are estimated. In addition, we provide a computational complexity analysis in terms of complex multiplications here. In the traditional EM algorithm, the computational time of the algorithm is mainly spent in solving (27). For the grid point search method, it requires approximately O(LKQ 0 MN M 1 N 1 ) operations, where M 1 and N 1 denote the search points of the Doppler and echo delay, respectively. Moreover, the estimation precision of the delay and Doppler is restricted by the grid interval. That is, for the grid point search method, the computational complexity must be increased to improve the estimation accuracy. For the GAEM algorithm in this study, since the gradient descent algorithm is used to solve (27), each iteration requires approximately O(LKQ 0 MN) operations. In addition, the convergence of GAEM algorithm can be achieved in a finite number of iterations. Therefore, compared with the traditional EM algorithm, the GAEM algorithm proposed in this study is faster.

Stage 2: Target Parameters and System Deviations Estimation
The estimated echo delay and Doppler in Section 2.2 is denoted by the vectorγ. Since the GAEM algorithm proposed in this paper is essentially an ML algorithm, the ML estimateγ is asymptotically distributed according to N γ, J −1 γ . J γ represents the Fisher information matrix (FIM) ofγ. Section 2.4 derives the expression of J γ and the SNR required to calculate J γ are estimated in Section 2.2.2.
Let γ = g(α). Due to the asymptotic statistical property ofγ, the estimator of target parameters and system deviations is arg min α (γ − g(α)) T J γ (γ − g(α)), which can be reformulated as arg min Equation (29) can be solved as follows. Let J γ γ i = J γ g(α i ) and J γγ = J γ g(α i+1 ), where α i is the result of the ith iteration. Then, a Taylor series expansion of J γ g(α i+1 ) around α i is given as Set ∂g(α) ∂α α=α i = P i . Then, we have Based on (29), a rough estimation result of the target position/velocity can be calculated when the radar system deviation is set to 0. The rough estimation result of the target position/velocity and the zero system deviation are integrated into the initial value α 0 of (31).
At this point, a robust and fast self-calibration algorithm is realized. Robustness implies that the algorithm can be applied to non-ideal scenarios such as superimposed signals and SNR fluctuations. Quickness implies that we reduce not only the computational complexity within each iteration of the EM algorithm, but also the overall number of iterations of the EM algorithm.

Cramér-Rao Bound in the Non-Ideal Environment
From the definition of CRB [31], we can get Using the chain rule, From the derivation of [32], The elements of J(ψ) are calculated in Appendix A. (σ wl ) 2 in (36) is the variance of Gaussian noise w lm [n]. Based on the definition of Section 2.1, ς lm [u] is the Fourier transform of w lm [n], and (σ l ) 2 is the variance of Gaussian noise ς lm [u]. From (7), the variances mentioned above have the following functional relationship: For clarity, J(ψ) is expressed as From the definition of Schur complement [33], J γ proposed in Section 2.3 can be expressed as

Results
In this section, the effectiveness of the proposed GAEM and WILS algorithm is verified by numerical simulations. In order to show the performance of the algorithm briefly, we define four radar nodes, and the distance between two radar nodes is 500 m. The angles of the four targets are −45 • , 60 • , 30 • and 35 • , respectively, and the distances between the four targets and first radar node is 2000 m, 2010 m, 2020 m and 2030 m, respectively. Both radar nodes and targets are moving.

Estimation of Time Delay and Doppler
In this subsection, the SNR is 20 dB. The bandwidth of radar signal is 10 MHz, so the 3 dB width of the sinc function in the time domain is 15 m. However, the closest distance between targets is 10 m, which generates the superimposed signal. The simulation of the superimposed signal is shown in Figure 5. Figure 5 shows the signal envelope in time domain after pulse compression. The dotted curve is the echo signal reflected from the first target. The dashed curve is the echo signal reflected from the second target. The solid curve is the superimposed signal of two targets. As can be seen from Figure 5a, when the distance between two targets is 60 m, the target distance corresponding to the superimposed signal peak is shifted from the real target distance. According to Figure 5b, when the distance between two targets is 30 m, the number of superimposed signal peaks becomes 1. The distance corresponding to the superimposed signal peak is still offset compared to the real target distance. Therefore, if the traditional ML algorithm is used to estimate the echo delay, the superimposed signals will affect the target parameter estimation.
The proposed GAEM is used to solve this problem. We compare GAEM algorithm with a gradient descent method based on traditional ML [34], and SQUAREM algorithm [24] is added to make GAEM converge faster. The root mean square error (RMSE) is used to evaluate the algorithm performance, as shown below where MC represents the number of Monte Carlo simulations; p (m) esti represents the estimated value; p real represents the real value. Due to page limitations, only the delay and Doppler estimation results of the first target are shown. Figure 6 shows the RMSE of the estimation results against the number of iterations.  The line with square markers represents the echo delay estimate only using gradient descent algorithm. The line with diamond markers represents the echo delay estimate using GAEM algorithm proposed in this paper. It can be seen from the Figure 6a that GAEM algorithm is more robust than traditional "hill climbing" algorithm when the target echoes are superimposed. This is because compared with gradient descent, GAEM algorithm has E-step, that is, the adaptive decomposition of multi-target echoes based on the estimation results of the previous M-step. The line with square markers represents the estimation results after combining the proposed GAEM algorithm with SQUAREM algorithm. The simulation results indicate that the SQUAREM algorithm can further reduce the number of iterations of the GAEM algorithm. Figure 6b shows the estimation results of Doppler shift. The curve definition and conclusions in Figure 6b are consistent with those in Figure 6a.
The computational complexity curves versus the grid resolutions are depicted in Figure 7. According to Figure 7, the GAEM algorithm is much more efficient than the existing EM method.

Target Parameters and System Deviations Estimation
This subsection simulates the WILS algorithm proposed in Section 2.3. Since MIMO radar with widely separated antennas uses different radar nodes as receivers and has a large baseline, the noise variance and the echo coefficient in different propagation paths are considered to be different. In other words, the SNRs of signal echo in different propagation paths are different. In the simulation of this section, the SNRs of the three transmitting and three receiving radar nodes with the same target are set as SNR 0 [0.01, 0.6, 1, 0.8, 0.02, 1, 0.4, 0.06, 1]. The SNR mentioned in the simulation figures refers to the largest item (i.e., SNR 0 ) in the above vector. When the SNR is 20 dB, the estimation results of target parameters and radar system deviations by GAEM and WILS algorithm are shown in Figure 8. Figure 8 is also the deployment of radar nodes and targets.
We compare the proposed WILS algorithm with the traditional LS algorithm [25], and add CRB simulation to show the performance of the proposed WILS algorithm. The simulation results of target parameter estimation are shown in Figure 9. The RMSE curves of WLIS algorithm are close to CRB curves. In contrast, the RMSE of the ILS algorithm lies above the CRB and the RMSE of the WILS algorithm. This is because the WILS algorithm adds the use of the SNR in different propagation paths, which ensures WILS has a more robust performance than ILS for non-ideal environment such as different noise variances of radar nodes and anisotropic targets. The noise variances and echo amplitudes are estimated in GAEM algorithm, which ensures the use of SNR in WILS algorithm.
The simulation results of radar system deviation estimation are shown in Figure 10. The conclusions of Figure 10 are identical to those in Figure 9a.
In summary, the proposed technique enables efficient and robust multi-target parameter estimation for MIMO radar with widely separated antennas on moving platforms.

Discussion
This paper deals with the problems of dense multi-target, SNR fluctuations of different propagation paths and high computational complexity of MIMO radar with widely separated antennas on moving platforms. In order to reduce the computational complexity, the gradient descent method is used to optimize the M-step in the EM algorithm. The simulation results indicate that the combination of the gradient descent and EM algorithm is feasible. Further, the echo amplitudes and noise variances are estimated and the WILS algorithm is used to ensure the robustness of the target parameter estimation algorithm. Moreover, the asymptotic property of the ML estimator ensures the realization of the notion that an echo with a high SNR should be assigned a large weight.
This study focuses on solving the practical problems encountered in the parameter estimation of MIMO radar with widely separated antennas on moving platforms. However, the work in this paper still has potential for further development. Firstly, inspired by Ref. [9], co-located MIMO antennas can be deployed in each platform. For this novel radar system, the advantages of both MIMO radar with widely separated antennas and co-located MIMO radar can be exploited, and waveform optimization can be used to further improve the performance of target parameter estimation. Secondly, for MIMO radar with widely separated antennas, target localization accuracy is affected by the topology of MIMO radar and targets [12]. MIMO radar with widely separated antennas on moving platforms naturally has the advantage of changing the radar topology in real time. How to optimize the radar topology configuration based on the estimated multi-target positions in this paper to further improve the target positioning accuracy is the area that needs to be studied in the future. Finally, it can be seen from [35] that when the number of parameters to be estimated is determined, adding more information will improve the estimation accuracy. Ref. [26] combines the target echo delay and direct wave delay to complete the estimation of radar system deviations. In the future, we can model the fast-time and slow-time of target echo and direct wave, and further improve the parameter estimation accuracy based on the increased direct wave information.

Conclusions
In this paper, a robust parameter estimation algorithm is proposed to achieve target localization in scenarios such as superimposed signals and SNR fluctuations of different propagation paths. Firstly, in order to estimate the echo delay/Doppler and SNR, this paper derives the Q function of EM algorithm. Secondly, in order to reduce the computational complexity of M-step in EM algorithm, the iterative gradient algorithm is used to realize the idea of the generalized EM algorithm. Furthermore, based on the estimation results of the previous method, the WILS algorithm is used to estimate the target positions and velocities. In addition, this paper deduces the CRB in the non-ideal environment. Finally, the simulation results verify the effectiveness of the algorithm proposed in this paper. Possible future works include the combination of MIMO radar with widely separated antennas and co-located MIMO radar, optimization of radar topology, and use of direct wave signals.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author, F.L., upon reasonable request.

Acknowledgments:
The authors would like to thank the editor and anonymous reviewers for their helpful comments and suggestions.

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

Appendix A. Derivation of the Fisher Information Matrix
Based on (35)-(37), the elements of J(ψ) can be calculated as follows. respectively. β k in (A12) is the effective bandwidth, which satisfies where B is the signal bandwidth; S k ( f ) denotes the Fourier transform of s k (t). In this paper, the SNR relative to the kqlth propagation path is defined as where N 0l is the noise power spectral density, which equals σ 2 wl f s . In the scenario of MIMO radar with widely separated antennas, the target reflection coefficients and path losses of different propagation paths are different. At the same time, the noise power spectral densities of different radar receivers are different. Therefore, it can be seen from (A20) that the SNRs of different propagation paths are different.