An Efficient Estimator for Source Localization Using TD and AOA Measurements in MIMO Radar Systems

3-D target localization based on multi-input and multioutput radar systems is studied in this letter by utilizing time delay and angle of arrival measurements, which can be obtained by some local sensor signal processing methods. For such a positioning problem, the typical closed-form solution first establishes pseudolinear equation and then, applies the weighted least squares (WLS) technique to determine the final position. However, the WLS solution often suffers from performance degradation when the noise is large. To alleviate this problem, the cost function of the WLS is represented, and a quadratic constraint is imposed such that the expectation of the cost function can attain the minimum value at the true position. Moreover, simulation experiments show that the proposed method performs better than the state-of-art algorithms.


I. INTRODUCTION
Target positioning based on multiple-input and multiple-output (MIMO) radar systems with distributed antennas [1] has attracted considerable attention in recent years. Localization methods [2]- [7] in MIMO systems can be attributed to use different measurements, including time delay (TD), bistatic range (BR), angle of arrival (AOA), and their combinations. The Doppler shift (DS) and BR rate (BRR) can be also utilized when the source and/or receivers are moving [8]- [11].
Similar to the hyperbolic location [12], the maximum likelihood (ML) estimator of target positioning in MIMO radar systems leads to some nonlinear and nonconvex optimization problems, whose global solutions are difficult to obtain. Pioneers have researched many methods for solving these problems. By introducing nuisance parameters, [2] and [3] give final estimate solutions based on multistage WLS technique. In [4], the localization problem is formulated as a constrained weighted least squares (CWLS) problem, which is a quadratically constrained quadratic programming (QCQP) problem, and the QCQP problem is then converted into the linear constrained quadratic programming (LCQP) problem, which can provide the closed-form solution. Recently, some efforts toward joint transmitter and target localization, joint synchronization and localization, and localization under uncertainty in distributed MIMO radars for both stationary and moving target are also made [8]- [11]. These closed-form methods can reach the Cramer-Rao lower bound (CRLB) [13] performance at low noise levels.
Compared to the methods using only TD measurements, hybrid TD-AOA methods perform better. In [5], a unit vector is introduced to eliminate the extra variables and then the problem is solved by WLS. In [6], first, the pseudolinear equation is established with nuisance parameters. Then, the second step refines the final solution by estimating the error of the first step to improve the estimated accuracy. Recently, the localization problem is recast into a convex optimization problem and solved by the semiclosed form solution in [7]. All these methods can achieve CRLB at low noise levels. However, performance degradation occurs to method [5] at a high noise level due to the presence of measurement noises in both the regressor and the regressand. Similarly, the squaring operations of the second stage in [6] and the rooting operations of method [7] may also lead to increased bias as well as poor positioning accuracy when the noise is large. Thus, it is necessary to solve these problems.
In this letter, we propose an explicit solution based on TD and AOA measurements for target localization in MIMO radar systems. To the best of our knowledge, the expectation of the cost function is not considered in the traditional WLS solution [5], [6]. To better handle the randomness of the regressor and regressand introduced by noise in the WLS solution, the proposed method aims to minimize the cost function by imposing a quadratic constraint such that the expectation of the cost function will attain its minimum value in the case of the true solution. Moreover, the proposed method can provide a closedform solution. Numerous simulations confirm the effectiveness of the proposed method.
We shall follow the convention that bold upper and lower case letters denote matrix and column vector. a(i, 1) is the ith element of a, and A(i, :) is the ith row of A. diag(a) denotes a diagonal matrix formed by elements of the vector a, and blkdiag(A, . . . , B) represents the block-diagonal matrix formed by A, …, B. In addition, ||a|| is the Euclidean norm of a and E ( * ) is the expectation of * . O M×N denotes the M × N matrix of zero. a o is the true value of a noisy quantity a , and the symbol ⊗ stands for Kronecker product. and receivers are precise and available. Besides, the signals obtained directly and indirectly by each receiver can produce one AOA pair and M TDs. Let the true AOA pair at the nth receiver be denoted by (θ o n , φ o n ) , where θ o n represents the azimuth angle and φ o n denotes the elevation angle. These true angle values are related to the coordinates of the source and receiver n through the following equations:

II. LOCALIZATION MODEL
The measured versions of AOA pairs at the nth receiver are modeled as θ n = θ o n + θ n and φ n = φ o n + φ n , where θ n and φ n are the error terms. By combing all the AOA measurements introduced above, the AOA measurement vectors can be constructed by where θ and φ are assumed to be zero-mean Gaussian noise vectors with covariance matrices After some simple operations, each true TD can be converted into a BR r o m,n , and it is expressed by where R o sn = ||u − s n || and R o tm = ||u − t m ||. Collecting all the true BRs into a matrix form, we have

A. Traditional Closed-Form Solution
The closed-form solution [5] introduces a unit vector of the actual target position with respect to the nth receiver for eliminating the extra variables. This unit vector is denoted as In addition, it also gives two vectors orthogonal to d o n through the following equations: Thus, a set of linear equations are established below by utilizing (5) and (6) and the details of the derivation can be found in [5] where h = [h T r , h T aoa ] T and G = [G T r , G T aoa ] T . G aoa and h aoa are expressed as The kth entry of h r and the kth row of G r are given by where k = (m − 1)N+n, n = 1, . . . , N, and m = 1, . . . , M. Substituting TD and AOA measurements into (9) yieldŝ where the regressandĥ and regressorĜ are the same as the previous defined in (7) except that the true values are replaced by their measurement versions. η = [ r T , θ T , φ T ] T represents the total noise vector with the covariance matrix Q η = E ( η T η), and B is given by Notably, calculating the weight matrix W requires the true target position. Method [5] first approximates W = Q −1 η and uses (11) to obtain an initial position. This initial value is then used to update W so that a better estimated position can be found.

B. Proposed Method
The closed-form algorithm summarized above can achieve CRLB under a low noise level. However, this solution suffers from performance degradation at high noise levels due to the presence of noise disturbances in bothĥ andĜ. Therefore, inspired by the BiasRed introduced in [14], a new approach for improving the estimation accuracy is proposed in this subsection.
The proposed method starts from introducing the augmented matrix A = [−Ĝ,ĥ] and vector v = [u T , 1] T . To the end, the cost function J can be rewritten as The matrix A can be decomposed as is the true augmented matrix and A is the noise matrix. By subtracting the true value A o from A and considering only the linear noise terms, A can be reformulated as A = TF, where T = blkdiag(T 1 , T 2 ). The matrices T 1 and T 2 can be obtained by The matrix F can be decomposed into three parts as follows: with The ith row of f 2 is f 2 (i, : . . . , N. Besides, the ith row of f 1 and f 3 , respectively, are computed by Obviously, substituting A = A o + A into (12) yields the following expression: The second term on the right side of (19) is a mean of zero. Therefore, by taking the expectation of J, we have Therefore, E (J ) will take its minimum value at the true value v o when the last term on the right side of (20) is constrained to a constant κ [14]. Thus, the solution for v can be expressed as where the constant κ can be any positive value and its value causes only a different scaling of v [14] and = E ( A T W A). Substituting A = TF into yields = F T PF, where P = E (T T WT). After that, by dividing the weighting matrix W into (M + 2) 2 N × N blocks with the (i, j) block denoted as W i j , i, j = 1, 2, . . . , M + 2, the matrix P can be expressed by P = blkdiag(P 1 , P 2 ) ( 2 2 ) where Similarly, P 2 is partitioned into (M + 2) block matrices denoted as P 2 = [P T 21 , . . . , P T 2(M+1) , P T 2(M+2) ] T , and the ith block matrix is expressed as where i = 1, . . . M, M + 2 and P 2(M+1) is also given by A feasible approach for solving (21) is to exploit Lagrange multiplier λ. After taking its derivative with respect to v and setting it to zero, we have the following equation: Multiplying v T on both sides of (26) and using the constrained Consequently, the smallest value of the cost function J is the minimum generalized eigenvalue λ min of the pair (A T WA, ) when κ = 1, and the related parameter vector v is the generalized eigenvector corresponding to λ min . Thus, the source location is given by Another point to note is that the true values are needed to construct W. We first approximate the true values in F by using their noisy versions and set the weighting matrix W = Q −1 to obtain an initial estimate. Then, this initial value can be employed to get a better estimated position of u o . Simulations show that the performance loss due to this operation is negligible for such problems.

IV. SIMULATION
In this section, we exhibit the performance of the proposed method and compare it with methods [5]- [7]. The simulation results are averaged by L = 10000 Monte Carlo experiments. Two indicators are used to evaluate the estimated performance. One is the root-mean-square error (RMSE), defined as RMSE(u) = L l=1 ||u l − u o || 2 /L, and the other is estimated bias, defined by bias(u) = || L l=1 u l /L − u o || 2 , where u l is the estimated value of u o at the lth run. Meanwhile, the CRLB (the details of the derivation can be found in [5]) is also considered as a benchmark in the RMSE comparison. The positions of transmitters and receivers are given in Table 1. One hundred targets are randomly generated on a circle centered at the origin of the coordinate system with a radius of 1000 m to obtain the average RMSE and bias. Consistent with [5] and [6], the covariance matrices, respectively, are set to be Q r = σ 2 r I MN , Q θ = σ 2 θ I N , and Q φ = σ 2 φ I N , where σ 2 r = 1600σ 2 and σ 2 θ = σ 2 φ = 0.1σ 2 . σ 2 is the noise-related multiplier that varies from −60 to 10 dB. For simplicity, the constant κ is set to 1. Fig. 1 compares the RMSE of different algorithms, from which we find that all the methods perform well at low noise levels. However, the method [5] suffers performance loss when σ 2 > −30 dB, which is mainly because the measurement noises dominate the performance. The method [6] and method [7] also start to deviate from the CRLB when σ 2 > −10 dB. In contrast, our proposed method achieves near-CRLB performance over small to relatively high noise levels, and the RMSE performance of the proposed method improves, averaging 5.7 dB compared to that of the others methods.
Next, we focus on the bias performance in Fig. 2. It is clear that   the bias of the method [5] is significantly higher than that of the other methods when σ 2 > −40 dB. The squaring operations of the second stage in [6] and the rooting operations of the method [7] lead to increased estimation bias when σ 2 > −10 dB. Furthermore, the proposed method always has the smallest bias over a wide range of noise levels. Finally, the complexity analysis of different methods is also roughly evaluated via average CPU run times in Table 2. The proposed method does not significantly increase the computational complexity while maintaining good performance.

V. CONCLUSION
Aiming at the target localization in MIMO radar systems, we present an efficient estimator using TD and AOA measurements in this letter. The cost function is reformulated and subjected to a quadratic constraint so that the expectation of the cost function can reach its ideal state at the true position. Moreover, the proposed method is a closed-form solution. Simulation results demonstrate that the proposed method outperforms the existing methods under the Gaussian noise model.