Robust Localization in Distributed MIMO Radar Using Delay and Angle Measurements with Impulsive Noise Robust TD/AOA Localization in Impulsive Noise

Abstract

This paper considers target localization using time delay (TD) and angle of arrival (AOA) measurements in distributed multiple-input multiple-output (MIMO) radar. Aiming at the problem that the localization performance of existing algorithms degrades sharply in the presence of impulsive noise, we propose a novel localization algorithm based on ${\ell }_p$-norm minimization and iteratively reweighted least squares (IRLS). Firstly, the TD and AOA measurement equations are established in the presence of zero-mean symmetric α-stable noise; then, the localization problem is transformed to a ${\ell }_p$-norm minimization problem by linearizing the measurement equations; and finally, the ${\ell }_p$-norm minimization problem is solved using IRLS by which the target position estimate is obtained, and the optimal choice of norm order $p$ is deduced. Moreover, the Cramér–Rao bound (CRB) for target position estimation in impulsive noise is also derived, generalizing the Gaussian CRB. Simulation results demonstrate that the proposed algorithm outperforms existing algorithms in terms of localization accuracy and robustness in impulsive noise.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Multiple-input multiple-output (MIMO) radar is a new kind of sensing system, which sends mutually orthogonal waveforms from multiple transmit antennas and extracts these waveforms from each of the receive antennas by a set of matched filters. This kind of radar system has an enlarged virtual receive aperture and a finer spatial resolution compared with the conventional radar systems [1]. According to the transmit/receive antenna configuration, MIMO radars can be divided into two categories: colocated MIMO radar with closely spaced antennas [2–4] and distributed MIMO radar with widely separated antennas [5–7]. Both architectures have their respective advantages, and this paper considers the target localization in the latter case.

Time delay (TD) and angle of arrival (AOA) are commonly used types of measurements for target localization in distributed MIMO radar. The TD measurement traces out an ellipsoidal surface for the possible target positions with foci located at the transmit and receive antennas. The AOA measurement induces a line from the receive antenna to the target. Theoretically, the target position can be estimated as the intersection of the lines and ellipsoids corresponding to the TD and AOA measurements. However, it is far from straightforward to estimate the target position from the TD and AOA measurements since both TD and AOA are highly nonlinear with respect to the target position.

In recent years, some effort has been devoted to this challenging problem. Borrowing the well-known two-stage weighted least squares (2WLS) idea by Chan and Ho [8], A. Noroozi et al. presented an algebraic algorithm for 3D hybrid TD/AOA localization in [9], where a pseudolinear set of equations is established by introducing nuisance parameters and the dependencies of nuisance parameters on target position are employed to yield a final estimate. Unlike Noroozi’s algorithm in [9], which is multistage estimators, R. Amiri et al. developed a different algebraic solution in [10] which uses the AOA measurements to linearize the TD measurement equations and identifies the target position in only one WLS stage. Both Noroozi’s algorithm and Amiri’s algorithm are shown analytically and confirmed by numerical simulations to attain the Cramér–Rao bound (CRB) under small measurement noise conditions. However, they perform unsatisfactorily at large noise levels and suffer from “threshold effect.” More recently, in order to improve the target localization performance at large noise levels, S. A. R. Kazemi et al. proposed an efficient convex solution for target position estimation in [11], where the associated localization problem is formulated as a nonconvex constrained quadratic problem and then recast as a convex problem, from which the target position estimate is determined by using polynomial root-finding. Kazemi’s algorithm is shown numerically to outperform previous algorithms and reach the CRB up to relatively large measurement noise levels. Nevertheless, the abovementioned algorithms are designed based on such assumption that the TD/AOA measurement noises are subject to Gaussian distribution, implicitly or explicitly. Gaussian distribution is a noise model widely accepted and employed in radar and communication fields. Its second-order statistics bring significant convenience to least squares (LS)-based algorithm development. In reality, the measurement noise is not always Gaussian distributed but often presented in a more impulsive nature [12–15]. Under such impulsive measurement noise, the abovementioned algorithms based on Gaussian noise assumption and LS approach will suffer performance degradation and even invalidation because the LS approach based on ${\ell }_2$-norm minimization of errors is very sensitive to outliers. Hence, there is yet a need for developing a robust algorithm for target localization in distributed MIMO radar using TD and AOA measurements with impulsive noise.

The α-stable distribution based on the generalized central limit theorem is a generalization of the Gaussian distribution and provides a better model for the measurement noise with a more impulsive nature [16, 17]. Therefore, it has been widely used to describe impulsive noises such as burst noise in indoor localization applications, clutter returns in radar applications, and man-made noise in acoustic localization [18–21]. On the other hand, when the TD/AOA measurement noise is α-stable distributed, existing algorithms based on the Gaussian noise assumption and LS approach will produce unreliable estimate since impulsive noise has no second-order moments. The ${\ell }_p$-norm minimization with $p<2$ is an attractive alternative to the LS approach and can be used to achieve robust estimation in impulsive noise because it is less sensitive to outliers. However, it is unrealistic to obtain an analytical solution for the ${\ell }_p$-norm minimization problem since the ${\ell }_p$-norm minimization is based on the $p$-order statistics. Iteratively reweighted least squares (IRLS) which converts the ${\ell }_p$-norm minimization problem into ${\ell }_2$-norm minimization problem for solving [22] is an attractive solver to deal with the ${\ell }_p$-norm minimization and has been successfully applied to parameter estimation in impulsive noise [23] in recent years. For the target localization in distributed MIMO radar using TD and AOA measurements, however, IRLS approach cannot be applied directly because the target position is nonlinearly related to TD and AOA measurements.

Motivated by the above facts, we investigate in this paper the problem of target localization in distributed MIMO radar using TD and AOA measurements with impulsive noise. Based on the ${\ell }_p$-norm minimization and IRLS approach, a robust algebraic solution for target position is proposed. Firstly, the TD and AOA measurement equation is established in the presence of zero-mean symmetric α-stable noise; then, the target localization problem is converted to a ${\ell }_p$-norm minimization problem by linearizing the TD and AOA measurement equations; finally, by using IRLS approach, the ${\ell }_p$-norm minimization problem is solved and the target position estimate is obtained. The best choice of $p$ and the CRB in impulsive noise will also be derived. Numerical simulations are also performed to verify the superiority of the proposed algorithm over existing algorithms [9–11] in impulsive measurement noise.

The paper is organized as follows. Section 2 presents the localization scenario and introduces the symbols involved. Section 3 presents a robust algebraic solution for the localization problem, and Section 4 derives the CRB. Section 5 contains the simulation results to evaluate the localization performance of the proposed algorithm, and Section 6 is the conclusion.

2. Problem Formulation

We consider finding the position of a target in 3D space with a distributed MIMO radar. The basic idea is the transmit antennas produce electromagnetic waves which are reflected or scattered by a single target; the receive antennas are utilized to collect the target echo, from which the TD and AOA measurements can be extracted to estimate the target position.

2.1. Measurement Equation

As illustrated in Figure 1, assume the distributed MIMO radar system is equipped with $M$ geographically separated receive antennas at positions ${\mathbf{s}}_{t,m}={x_{t,m},y_{t,m},z_{t,m}}^{\text{T}}$, $m=\mathrm{1,2},\dots ,M$, and $N$ widely spaced receive antennas at positions ${\mathbf{s}}_{r,n}={x_{r,n},y_{r,n},z_{r,n}}^{\text{T}}$, $n=\mathrm{1,2},\dots ,N$. The transmit antennas send out a set of waveforms, which are then reflected by the target at unknown position $\mathbf{u}={x,y,z}^{\text{T}}$.

[figure omitted; refer to PDF]

Based on the above geometry, the true AOA pair for receive antenna $n$, i.e., the elevation angle denoted by ${\theta }_n$ and the azimuth angle denoted by ${\varphi }_n$, is given by$$\begin{array}{c}\text{(1)}& {\theta }_n=\arctan \frac{y-y_{r,n}}{x-x_{r,n}},{\varphi }_n=\arctan \frac{z-z_{r,n}}{\sqrt{{x-x_{r,n}}^2+{y-y_{r,n}}^2}}.\end{array}$$

The range between transmit antenna $m$ and the target is $R_{t,m}=\mathbf{u}-{\mathbf{s}}_{t,m}$, the range between receive antenna $n$ and the target is $R_{r,n}=\mathbf{u}-{\mathbf{s}}_{r,n}$, and the baseline distance between transmit antenna $m$ and the receive antenna $n$ is $R_{t,m,r,n}={\mathbf{s}}_{t,m}-{\mathbf{s}}_{r,n}$. Thus, the TD between the direct path signal from transmit antenna $m$ and the corresponding reflected signal arriving at receive antenna $n$ can be expressed as$$\begin{array}{}\text{(2)}& {\tau }_{m,n}=\frac{1}{c}{R}_{t,m}+R_{r,n}-R_{t,m,r,n}.\end{array}$$

Considering the unavoidable measurement noises in reality, we have the erroneous TD and AOA measurements as$$\begin{array}{}\text{(3)}& {\tilde{\theta }}_n={\theta }_n+e_{\theta ,n},\text{(4)}& {\tilde{\varphi }}_n={\varphi }_n+e_{\varphi ,n},\text{(5)}& {\tilde{\mathbf{\tau }}}_{m,n}={\tau }_{m,n}+e_{\tau ,m,n},\end{array}$$where ${\tilde{\theta }}_n$, ${\tilde{\varphi }}_n$, and ${\mathbf{\phi }}_n$ represent the noisy TD and AOA measurements and $e_{\theta ,n}$, $e_{\varphi ,n}$, and $e_{\tau ,m,n}$ represent the corresponding measurement noises.

There are $N$ elevation angle measurements, $N$ azimuth angle measurements, and $MN$ TD measurements altogether. For easier manipulation, collect these TD and AOA measurements in column vector form as follows:$$\begin{array}{}\text{(6)}& \tilde{\mathbf{\alpha }}=\mathbf{\alpha }+\mathbf{e},\end{array}$$where $\tilde{\mathbf{\alpha }}={{\tilde{\mathbf{\theta }}}^{\text{T}},{\tilde{\mathbf{\varphi }}}^{\text{T}},{\tilde{\mathbf{\tau }}}^{\text{T}}}^{\text{T}}$, $\mathbf{\alpha }={{\mathbf{\theta }}^{\text{T}},{\mathbf{\varphi }}^{\text{T}},{\mathbf{\tau }}^{\text{T}}}^{\text{T}}$, and $\mathbf{e}={{\mathbf{e}}_{\theta }^{\text{T}},{\mathbf{e}}_{\varphi }^{\text{T}},{\mathbf{e}}_{\tau }^{\text{T}}}^{\text{T}}$, with$$\begin{array}{c}\text{(7)}& \tilde{\mathbf{\theta }}={{\tilde{\theta }}_1,{\tilde{\theta }}_2,\dots ,{\tilde{\theta }}_N}^{\text{T}},\mathbf{\theta }={{\theta }_1,{\theta }_2,\dots ,{\theta }_N}^{\text{T}},\\ {\mathbf{e}}_{\theta }={e_{\theta ,1},e_{\theta ,2},\dots ,e_{\theta ,N}}^{\text{T}},\\ \tilde{\mathbf{\varphi }}={{\tilde{\varphi }}_1,{\tilde{\varphi }}_2,\dots ,{\tilde{\varphi }}_N}^{\text{T}},\\ \mathbf{\varphi }={{\varphi }_1,{\varphi }_2,\dots ,{\varphi }_N}^{\text{T}},\\ {\mathbf{e}}_{\varphi }={e_{\varphi ,1},e_{\varphi ,2},\dots ,e_{\varphi ,N}}^{\text{T}},\\ \tilde{\mathbf{\tau }}={{\tilde{\mathbf{\tau }}}_1^{\text{T}},{\tilde{\mathbf{\tau }}}_2^{\text{T}},\dots ,{\tilde{\mathbf{\tau }}}_M^{\text{T}}}^{\text{T}},\\ {\tilde{\mathbf{\tau }}}_m={{\tilde{\tau }}_{m,1},{\tilde{\tau }}_{m,2},\dots ,{\tilde{\tau }}_{m,N}}^{\text{T}},\\ \mathbf{\tau }={{\mathbf{\tau }}_1^{\text{T}},{\mathbf{\tau }}_2^{\text{T}},\dots ,{\mathbf{\tau }}_M^{\text{T}}}^{\text{T}},\\ {\mathbf{\tau }}_m={{\tau }_{m,1},{\tau }_{m,2},\dots ,{\tau }_{m,N}}^{\text{T}},\\ {\mathbf{e}}_{\tau }={{\mathbf{e}}_{\tau ,1}^{\text{T}},{\mathbf{e}}_{\tau ,2}^{\text{T}},\dots ,{\mathbf{e}}_{\tau ,M}^{\text{T}}}^{\text{T}},\\ {\mathbf{e}}_{\tau ,m}^{\text{T}}={e_{\tau ,m,1},e_{\tau ,m,2},\dots ,e_{\tau ,m,N}}^{\text{T}}.\end{array}$$

2.2. Measurement Noise

It is assumed in previous studies [9–11] that the TD and AOA measurement noises are subject to zero-mean Gaussian distribution. However, in practical applications, the TD and AOA measurement noises in distributed MIMO radar are often presented in a more impulsive nature, which cannot be simply modelled as Gaussian distribution. As a generalization of the Gaussian distribution, α-stable distribution provides a more suitable model for the measurement noise with a more impulsive nature. Thus, we assume, in this paper, the TD and AOA measurement noises are subject to zero-mean symmetric α-stable (SαS) distribution, whose probability density function (PDF) cannot be written analytically but the general characteristic function (CF), that is, the Fourier transform of its PDF, can be expressed explicitly as follows [21]:$$\begin{array}{}\text{(8)}& \varsigma t=\exp -\gamma t^{\alpha },\end{array}$$where $\alpha \in \mathrm{0,2}$ known as the characteristic exponent determines the impulsiveness of the noises and the tail of the distribution. A smaller value of $\alpha$ will result in a heavier tail of PDF, a more impulsive nature, and more samples deviating from the mean. As $\alpha$ increases, the impulsiveness of the noises and the tail heaviness of the PDF reduces, with the upper bound $\alpha =2$ corresponding to the Gaussian distribution. Out of consideration for mathematical and practical concerns, it is usually assumed in SαS distribution processing that $\alpha \in \mathrm{1,2}$. We also follow the same assumption in this work. $\gamma >0$ known as the dispersion parameter determines the spread of the distribution around the mean. When $\alpha =2$, $\gamma$ is similar to the variance of the Gaussian distribution.

Figure 2 compares the noise samples generated from Gaussian distribution and those from SαS distribution under the same dispersion parameter. It can be seen that, unlike Gaussian distribution, SαS distribution has outliers far from the mean. These outliers cause the LS-based estimators to generate unreliable parameter estimate because the performance of the ${\ell }_2$-norm minimizer is very sensitive to outliers.

[figure omitted; refer to PDF]

For SαS distribution, there exists no finite second-order moment but only finite moments for orders less than $\alpha$ which are known as fractional lower-order moment (FLOM). Supposing $\vartheta$ is a random variable subject to zero-mean SαS distribution, the FLOM of $\vartheta$ is given as$$\begin{array}{}\text{(9)}& E{\vartheta }^q=C_{\alpha }q,\alpha {\gamma }^{q/\alpha },\hspace{1em}q\in -\mathrm{1,1}\cup 1,\alpha ,\end{array}$$where $\cup$ represents the union operator and$$\begin{array}{}\text{(10)}& C_{\alpha }q,\alpha =\frac{\Gamma q+1/2\Gamma -q/\alpha }{\alpha \sqrt{\pi }\Gamma -q/2}{2}^{q+1}\end{array}$$with$$\begin{array}{}\text{(11)}& \Gamma x={\displaystyle {\int }_0^{\infty }\exp -t{t}^{x-1}\text{d}t},\end{array}$$being the gamma function. Furthermore, it can be derived that$$\begin{array}{}\text{(12)}& E{\vartheta }^{q-1}\text{sign}\vartheta =0.\end{array}$$

In this work, we are interested in robustly identifying the unknown target position $\mathbf{u}$ from impulsive noise-contaminated TD and AOA measurements. Nevertheless, this is a potentially challenging task due to the nonlinearity of the TD/AOA measurements and the impulsive nature of the measurement noise.

3. Proposed Algorithm

3.1. Localization Objective Function

The derivation begins with converting the localization to a ${\ell }_p$-norm minimization problem. To achieve this, we first normalize the TD and AOA measurements in (3)∼(5) as$$\begin{array}{}\text{(13)}& {\tilde{\psi }}_{\theta ,n}={\psi }_{\theta ,n}+{\epsilon }_{\theta ,n},\text{(14)}& {\tilde{\psi }}_{\varphi ,n}={\psi }_{\varphi ,n}+{\epsilon }_{\varphi ,n},\text{(15)}& {\tilde{\psi }}_{\tau ,m,n}={\psi }_{\tau ,m,n}+{\epsilon }_{\tau ,m,n},\end{array}$$where$$\begin{array}{c}\text{(16)}& {\tilde{\psi }}_{\theta ,n}={\gamma }_{\theta ,n}^{-1/{\alpha }_{\theta ,n}}{\tilde{\theta }}_n,{\psi }_{\theta ,n}={\gamma }_{\theta ,n}^{-1/{\alpha }_{\theta ,n}}{\theta }_n,\\ {\epsilon }_{\theta ,n}={\gamma }_{\theta ,n}^{-1/{\alpha }_{\theta ,n}}{e}_{\theta ,n},\\ {\tilde{\psi }}_{\varphi ,n}={\gamma }_{\varphi ,n}^{-1/{\alpha }_{\varphi ,n}}{\tilde{\varphi }}_n,\\ {\psi }_{\varphi ,n}={\gamma }_{\varphi ,n}^{-1/{\alpha }_{\varphi ,n}}{\varphi }_n,\\ {\epsilon }_{\varphi ,n}={\gamma }_{\varphi ,n}^{-1/{\alpha }_{\varphi ,n}}{e}_{\varphi ,n},\\ {\tilde{\psi }}_{\tau ,m,n}={\gamma }_{\tau ,m,n}^{-1/{\alpha }_{\tau ,m,n}}{\tilde{\tau }}_{m,n},\\ {\psi }_{\tau ,m,n}={\gamma }_{\tau ,m,n}^{-1/{\alpha }_{\tau ,m,n}}{\tau }_{m,n},\\ {\epsilon }_{\tau ,m,n}={\gamma }_{\tau ,m,n}^{-1/{\alpha }_{\tau ,m,n}}{e}_{\tau ,m,n}.\end{array}$$

After the above processing, the measurement noises $e_{\theta ,n}$, $e_{\varphi ,n}$, and $e_{\tau ,m,n}$ are converted to normalized noises ${\epsilon }_{\theta ,n}$, ${\epsilon }_{\varphi ,n}$, and ${\epsilon }_{\tau ,m,n}$, which are independent identically distributed with unit noise dispersion.

By defining the following vectors,$$\begin{array}{c}\text{(17)}& \tilde{\mathbf{\psi }}={{\tilde{\mathbf{\psi }}}_{\theta }^{\text{T}},{\tilde{\mathbf{\psi }}}_{\varphi }^{\text{T}},{\tilde{\mathbf{\psi }}}_{\tau }^{\text{T}}}^{\text{T}},\mathbf{\psi }={{\mathbf{\psi }}_{\theta }^{\text{T}},{\mathbf{\psi }}_{\varphi }^{\text{T}},{\mathbf{\psi }}_{\tau }^{\text{T}}}^{\text{T}},\\ \mathbf{\epsilon }={{\mathbf{\epsilon }}_{\theta }^{\text{T}},{\mathbf{\epsilon }}_{\varphi }^{\text{T}},{\mathbf{\epsilon }}_{\tau }^{\text{T}}}^{\text{T}},\end{array}$$where$$\begin{array}{c}\text{(18)}& {\tilde{\mathbf{\psi }}}_{\theta }={{\tilde{\psi }}_{\theta ,1},{\tilde{\psi }}_{\theta ,2},\dots ,{\tilde{\psi }}_{\theta ,N}}^{\text{T}},{\mathbf{\psi }}_{\theta }={{\psi }_{\theta ,1},{\psi }_{\theta ,2},\dots ,{\psi }_{\theta ,N}}^{\text{T}},\\ {\mathbf{\epsilon }}_{\theta }={{\epsilon }_{\theta ,1},{\epsilon }_{\theta ,2},\dots ,{\epsilon }_{\theta ,N}}^{\text{T}},\\ {\tilde{\mathbf{\psi }}}_{\varphi }={{\tilde{\psi }}_{\varphi ,1},{\tilde{\psi }}_{\varphi ,2},\dots ,{\tilde{\psi }}_{\varphi ,N}}^{\text{T}},\\ {\mathbf{\psi }}_{\varphi }={{\psi }_{\varphi ,1},{\psi }_{\varphi ,2},\dots ,{\psi }_{\varphi ,N}}^{\text{T}},\\ {\mathbf{\epsilon }}_{\varphi }={{\epsilon }_{\varphi ,1},{\epsilon }_{\varphi ,2},\dots ,{\epsilon }_{\varphi ,N}}^{\text{T}},\\ {\tilde{\mathbf{\psi }}}_{\tau }={{\tilde{\mathbf{\psi }}}_{\tau ,1}^{\text{T}},{\tilde{\mathbf{\psi }}}_{\tau ,2}^{\text{T}},\dots ,{\tilde{\mathbf{\psi }}}_{\tau ,M}^{\text{T}}}^{\text{T}},\\ {\tilde{\mathbf{\psi }}}_{\tau ,m}={{\tilde{\psi }}_{\tau ,m,1},\dots ,{\tilde{\psi }}_{\tau ,m,N}}^{\text{T}},\\ {\mathbf{\psi }}_{\tau }={{\mathbf{\psi }}_{\tau ,1}^{\text{T}},\dots ,{\mathbf{\psi }}_{\tau ,M}^{\text{T}}}^{\text{T}},\\ {\mathbf{\psi }}_{\tau ,m}={{\psi }_{\tau ,m,1},{\psi }_{\tau ,m,2},\dots ,{\psi }_{\tau ,m,N}}^{\text{T}},\\ {\mathbf{\epsilon }}_{\tau }={{\mathbf{\epsilon }}_{\tau ,1}^{\text{T}},\dots ,{\mathbf{\epsilon }}_{\tau ,M}^{\text{T}}}^{\text{T}},\\ {\mathbf{\epsilon }}_{\tau ,m}^{\text{T}}={{\epsilon }_{\tau ,m,1},\dots ,{\epsilon }_{\tau ,m,N}}^{\text{T}},\end{array}$$and we can reorganize (13)∼(15) in vector form as$$\begin{array}{}\text{(19)}& \tilde{\mathbf{\psi }}=\mathbf{\psi }+\mathbf{\epsilon }.\end{array}$$

By now, the TD/AOA-based localization problem can be expressed as solving the following ${\ell }_p$-norm minimization:$$\begin{array}{}\text{(20)}& J\mathbf{u}={\tilde{\mathbf{\psi }}-\mathbf{\psi }}_p^p,\end{array}$$where $1<p<\alpha$ and ${\ast }_p$ represents the ${\ell }_p$-norm. Note that the ${\ell }_p$-norm objective function $J\mathbf{u}$ is essentially weighted by the normalization in (3)∼(5) according to the characteristic exponent parameter and the dispersion parameter of the TD/AOA measurement noises. To be more specific, the ${\ell }_p$-norm objective function $J\mathbf{u}$ can be further rewritten as$$\begin{array}{}\text{(21)}& J\mathbf{u}={\mathbf{\Lambda }\tilde{\mathbf{\alpha }}-\mathbf{\alpha }}_p^p,\end{array}$$where$$\begin{array}{c}\text{(22)}& \mathbf{\Lambda }=\mathrm{diag}{\mathbf{\Lambda }}_{\theta },{\mathbf{\Lambda }}_{\varphi },{\mathbf{\Lambda }}_{\tau },{\mathbf{\Lambda }}_{\theta }=\mathrm{diag}{\gamma }_{\theta ,1}^{-1/{\alpha }_{\theta ,1}},{\gamma }_{\theta ,2}^{-1/{\alpha }_{\theta ,2}},\dots ,{\gamma }_{\theta ,N}^{-1/{\alpha }_{\theta ,N}},\\ {\mathbf{\Lambda }}_{\varphi }=\mathrm{diag}{\gamma }_{\varphi ,1}^{-1/{\alpha }_{\varphi ,1}},{\gamma }_{\varphi ,2}^{-1/{\alpha }_{\varphi ,2}},\dots ,{\gamma }_{\varphi ,N}^{-1/{\alpha }_{\varphi ,N}},\\ {\mathbf{\Lambda }}_{\tau }=\mathrm{diag}{\mathbf{\Lambda }}_{\tau ,1},{\mathbf{\Lambda }}_{\tau ,2},\dots ,{\mathbf{\Lambda }}_{\tau ,M},\\ {\mathbf{\Lambda }}_{\tau ,m}=\mathrm{diag}{\gamma }_{\tau ,m,1}^{-1/{\alpha }_{\tau ,m,1}},{\gamma }_{\tau ,m,2}^{-1/{\alpha }_{\tau ,m,2}},\dots ,{\gamma }_{\tau ,m,N}^{-1/{\alpha }_{\tau ,m,N}}.\end{array}$$

By using the weighting matrix $\mathrm{\Lambda }$, the measurements with higher accuracy are given more weight. Thus, the optimization problem in (22) can be seen as a weighted least ${\ell }_p$-norm estimation problem, analogous to the weighted least squares estimation in the Gaussian noise.

3.2. Iteratively Reweighted Least Squares Estimator

In this subsection, we focus on solving the ${\ell }_p$-norm minimization presented in (22) by using IRLS approach. As mentioned above, part of the problem with applying IRLS approach to solve (22) is that the TD and AOA measurement equations are nonlinear with respect to target position. To circumvent this, we shall transform the TD and AOA measurement equations to linear equations as follows:$$\begin{array}{c}\text{(23)}& x\sin {\tilde{\theta }}_n-y\cos {\tilde{\theta }}_n=x_{r,n}\sin {\tilde{\theta }}_n-y_{r,n}\cos {\tilde{\theta }}_n+{\mathbf{b}}_{\theta ,\theta ,n}^{\text{T}}\mathbf{u}-{\mathbf{s}}_{r,n}{e}_{\theta ,n},\text{(24)}& x\cos {\tilde{\theta }}_n\sin {\tilde{\phi }}_n+y\sin {\tilde{\theta }}_n\sin {\tilde{\phi }}_n-z\cos {\tilde{\phi }}_n=x_{r,n}\cos {\tilde{\theta }}_n\sin {\tilde{\phi }}_n+y_{r,n}\sin {\tilde{\theta }}_n\sin {\tilde{\phi }}_n-z_{r,n}\cos {\tilde{\phi }}_n+{\mathbf{b}}_{\phi ,\theta ,n}^{\text{T}}\mathbf{u}-{\mathbf{s}}_{r,n}{e}_{\theta ,n}+{\mathbf{b}}_{\phi ,\phi ,n}^{\text{T}}\mathbf{u}-{\mathbf{s}}_{r,n}{e}_{\phi ,n},\\ \text{(25)}& 2{{\mathbf{s}}_{r,n}-{\mathbf{s}}_{t,m}+c{\tilde{\tau }}_{m,n}{\mathbf{\rho }}_n+R_{t,m,r,n}{\mathbf{\rho }}_n}^T\mathbf{u}=2{\mathbf{s}}_{r,n}^T{\mathbf{s}}_{r,n}+2c{\tilde{\tau }}_{m,n}{R}_{t,m,r,n}+c^2{\tilde{\tau }}_{m,n}^2+2c{\tilde{\tau }}_{m,n}+R_{t,m,r,n}{\mathbf{\rho }}_n^T{\mathbf{s}}_{r,n}+2c{\tilde{\tau }}_{m,n}+R_{t,m,r,n}{\mathbf{b}}_{\tau \text{,}\theta \text{,}n}^T\mathbf{u}-{\mathbf{s}}_{r,n}{e}_{\theta ,n}\\ +2c{\tilde{\tau }}_{m,n}+R_{t,m,r,n}{\mathbf{b}}_{\tau \text{,}\phi \text{,}n}^T\mathbf{u}-{\mathbf{s}}_{r,n}{e}_{\phi ,n}-2c{R}_{t,m}^o{e}_{\tau ,m,n},\end{array}$$where$$\begin{array}{c}\text{(26)}& {\mathbf{b}}_{\theta \text{,}\theta \text{,}n}={\cos {\tilde{\theta }}_n,\sin {\tilde{\theta }}_n,0}^{\text{T}}{\mathbf{b}}_{\phi \text{,}\theta \text{,}n}={-\sin {\tilde{\phi }}_n\sin {\tilde{\theta }}_n,\sin {\tilde{\phi }}_n\cos {\tilde{\theta }}_n,0}^{\text{T}},\\ {\mathbf{b}}_{\phi \text{,}\phi \text{,}n}={\cos {\tilde{\phi }}_n\cos {\tilde{\theta }}_n,\cos {\tilde{\phi }}_n\sin {\tilde{\theta }}_n,\sin {\tilde{\phi }}_n}^{\text{T}}\\ {\mathbf{b}}_{\tau \text{,}\theta \text{,}n}={-\cos {\tilde{\phi }}_n\sin {\tilde{\theta }}_n,\cos {\tilde{\phi }}_n\cos {\tilde{\theta }}_n,0}^T,\\ {\mathbf{b}}_{\tau \text{,}\phi \text{,}n}={-\sin {\tilde{\phi }}_n\cos {\tilde{\theta }}_n,-\sin {\tilde{\phi }}_n\sin {\tilde{\theta }}_n,\cos {\tilde{\phi }}_n}^T,\\ {\mathbf{\rho }}_n={\cos {\tilde{\phi }}_n\cos {\tilde{\theta }}_n,\cos {\tilde{\phi }}_n\sin {\tilde{\theta }}_n,\sin {\tilde{\phi }}_n}^T,\end{array}$$for $m=\mathrm{1,2},\dots ,M$ and $n=\mathrm{1,2},\dots ,N$.

Stack (23)∼(25) for $m=\mathrm{1,2},\dots ,M$ and $n=\mathrm{1,2},\dots ,N$ in vector form as$$\begin{array}{}\text{(27)}& \mathbf{B}\mathbf{e}=\mathbf{G}\mathbf{u}-\mathbf{h},\end{array}$$where$$\begin{array}{c}\text{(28)}& \mathbf{G}=\begin{array}{c}{\mathbf{G}}_{\theta }\\ {\mathbf{G}}_{\phi }\\ {\mathbf{G}}_{\tau }\end{array},\mathbf{h}=\begin{array}{c}{\mathbf{h}}_{\theta }\\ {\mathbf{h}}_{\phi }\\ {\mathbf{h}}_{\tau }\end{array},\\ \mathbf{B}=\begin{array}{ccc}{\mathbf{B}}_{\theta ,\theta }& {\mathbf{O}}_{N\times N}& {\mathbf{O}}_{N\times MN}\\ {\mathbf{B}}_{\phi ,\theta }& {\mathbf{B}}_{\phi ,\phi }& {\mathbf{O}}_{N\times MN}\\ {\mathbf{B}}_{\tau ,\theta }& {\mathbf{B}}_{\tau ,\phi }& {\mathbf{B}}_{\tau ,\tau }\end{array},\end{array}$$with the inner elements of the submatrix in (28) given by$$\begin{array}{c}\text{(29)}& {{\mathbf{G}}_{\theta }}_{n,1:3}=\sin {\tilde{\theta }}_n,-\cos {\tilde{\theta }}_n,0,{{\mathbf{G}}_{\phi }}_{n,1:3}=\cos {\tilde{\theta }}_n\sin {\tilde{\phi }}_n,\sin {\tilde{\theta }}_n\sin {\tilde{\phi }}_n,-\cos {\tilde{\phi }}_n,\\ {{\mathbf{G}}_{\tau }}_{m-1\times N+n,1:3}=2{{\mathbf{s}}_{r,n}-{\mathbf{s}}_{t,m}+c{\tilde{\tau }}_{m,n}{\mathbf{\rho }}_n+R_{t,m,r,n}{\mathbf{\rho }}_n}^T,\\ {{\mathbf{h}}_{\theta }}_{n,1}=x_{r,n}\sin {\tilde{\theta }}_n-y_{r,n}\cos {\tilde{\theta }}_n,\\ {{\mathbf{h}}_{\phi }}_{n,1}=x_{r,n}\cos {\tilde{\theta }}_n\sin {\tilde{\phi }}_n+y_{r,n}\sin {\tilde{\theta }}_n\sin {\tilde{\phi }}_n-z_{r,n}\cos {\tilde{\phi }}_n,\\ {{\mathbf{h}}_{\tau }}_{m-1\times N+n,1}=2{\mathbf{s}}_{r,n}^T{\mathbf{s}}_{r,n}+2c{\tilde{\tau }}_{m,n}{R}_{t,m,r,n}+c^2{\tilde{\tau }}_{m,n}^2+2c{\tilde{\tau }}_{m,n}+R_{t,m,r,n}{\mathbf{\rho }}_n^T{\mathbf{s}}_{r,n},\\ {{\mathbf{B}}_{\theta ,\theta }}_{n,n}={\mathbf{b}}_{\theta \text{,}\theta \text{,}n}^T\mathbf{u}-{\mathbf{s}}_{r,n},\\ {{\mathbf{B}}_{\phi ,\theta }}_{n,n}={\mathbf{b}}_{\phi \text{,}\theta \text{,}n}^T\mathbf{u}-{\mathbf{s}}_{r,n},\\ {{\mathbf{B}}_{\phi ,\phi }}_{n,n}={\mathbf{b}}_{\phi \text{,}\phi \text{,}n}^T\mathbf{u}-{\mathbf{s}}_{r,n},\\ {{\mathbf{B}}_{\tau ,\theta }}_{m-1\times N+n,n}=2c{\tilde{\tau }}_{m,n}+R_{t,m,r,n}{\mathbf{b}}_{\tau \text{,}\theta \text{,}n}^T\mathbf{u}-{\mathbf{s}}_{r,n},\\ {{\mathbf{B}}_{\tau ,\phi }}_{m-1\times N+n,n}=2c{\tilde{\tau }}_{m,n}+R_{t,m,r,n}{\mathbf{b}}_{\tau \text{,}\phi \text{,}n}^T\mathbf{u}-{\mathbf{s}}_{r,n},\\ {{\mathbf{B}}_{\tau ,\tau }}_{m-1\times N+n,m-1\times N+n}=-2c{R}_{t,m}^o,\end{array}$$for $m=\mathrm{1,2},\dots ,M$ and $n=\mathrm{1,2},\dots ,N$.

It can be deduced from (27) that $\mathbf{e}={\mathbf{B}}^{-1}\mathbf{G}\mathbf{u}-\mathbf{h}=\tilde{\alpha }-\alpha$. Substituting this fact into (21) leads to$$\begin{array}{}\text{(30)}& J\mathbf{u}={{\mathbf{G}}_1\mathbf{u}-{\mathbf{h}}_1}_p^p,\end{array}$$where ${\mathbf{G}}_1=\mathbf{\Lambda }{\mathbf{B}}^{-1}\mathbf{G}$ and ${\mathbf{h}}_1=\mathbf{\Lambda }{\mathbf{B}}^{-1}\mathbf{h}$. For the ${\ell }_p$-norm minimization problem presented in (30), we can transform it to ${\ell }_2$-norm minimization problem for solving. To accomplish this, we rewrite (30) in ${\ell }_2$-norm form as$$\begin{array}{}\text{(31)}& J\mathbf{u}={\mathbf{D}{\mathbf{G}}_1\mathbf{u}-{\mathbf{h}}_1}^2,\end{array}$$where$$\begin{array}{}\text{(32)}& \mathbf{D}={{\mathbf{G}}_1\mathbf{u}-{\mathbf{h}}_1}^{p-2/2}.\end{array}$$

Clearly, the weighted least squares solution for the ${\ell }_2$-norm minimization problem presented in (31) can be obtained as$$\begin{array}{}\text{(33)}& \mathbf{u}={{\mathbf{G}}_1^T\mathbf{W}{\mathbf{G}}_1}^{-1}{\mathbf{G}}_1^T\mathbf{W}{\mathbf{h}}_1,\end{array}$$where $\mathbf{W}$ is the weighting matrix given by$$\begin{array}{}\text{(34)}& \mathbf{W}={\mathbf{D}}^T\mathbf{D}.\end{array}$$

It is noted that $\mathbf{W}$ is related to the unknown parameter $\mathbf{u}$ and thus cannot be evaluated properly. To circumvent this problem, we can estimate the weighted matrix $\mathbf{W}$ by iteration and use iteratively the refined weighting matrix $\mathbf{W}$ to compute the target position. The process can be reformulated with more detail as follows:

(i) Start with an initial weighting matrix $\mathbf{W}=\mathbf{I}$, and then we have $\widehat{\mathbf{u}}={{\mathbf{G}}_1^T{\mathbf{G}}_1}^{-1}{\mathbf{G}}_1^T{\mathbf{h}}_1$

(ii) Repeat the following:

  Substitute $\widehat{\mathbf{u}}$ into (34) to compute a more accurate weighting matrix $\mathbf{W}$

  Substitute $\mathbf{W}$ into (33) to produce a more accurate estimate $\widehat{\mathbf{u}}$

  Stop the above iteration when the difference of target position estimate between two iterations is less than a specified threshold or the number of iterations reaches the specified value.

3.3. Choice of p

The norm order $p$ is an important parameter with great influence on the performance of the least ${\ell }_p$-norm estimator. In this subsection, we shall deduce the optimal choice of $p$.

Using (9) and (12), we have$$\begin{array}{}\text{(35)}& E\nabla J\mathbf{u}{\nabla }^{T}J\mathbf{u}=p^2{C}_{\alpha }2p-2,\alpha {\gamma }^{2p-2/\alpha }{\mathbf{G}}_1^T{\mathbf{G}}_1,\text{(36)}& E\mathbf{H}J\mathbf{u}=pp-1C_{\alpha }p-2,\alpha {\gamma }^{p-2/\alpha }{\mathbf{G}}_1^T{\mathbf{G}}_1,\end{array}$$where $\mathbf{H}J\mathbf{u}$ is the Hessian matrix.

By using (35) and (36), we can compute the covariance matrix of $\mathbf{u}$ as$$\begin{array}{}\text{(37)}& {\mathbf{C}}_{\alpha }\mathbf{u}=\frac{C_{\alpha }2p-2,\alpha }{{p-1}^2{C}_{\alpha }^{2}p-2,\alpha }{\gamma }^{2/\alpha }{{\mathbf{G}}_1^T{\mathbf{G}}_1}^{-1}.\end{array}$$

The scalar term of parameter variance of $\mathbf{u}$, characterized by $p$, is then obtained from (37) as$$\begin{array}{}\text{(38)}& {\Phi }_{\alpha }p=\frac{C_{\alpha }2p-2,\alpha }{{p-1}^2{C}_{\alpha }^{2}p-2,\alpha }.\end{array}$$

To determine the value $p$ that minimizes ${\mathrm{\Phi }}_{\alpha }p$, we take the derivative of (34) with respect to $p$ and set it to zero, which yields$$\begin{array}{}\text{(39)}& {\Phi }_{\alpha }^{\prime }p={\Phi }_{\alpha }p{\lambda }_{\alpha }p=0,\end{array}$$where$$\begin{array}{c}\text{(40)}& {\lambda }_{\alpha }p=\Psi 2-p-\Psi 2-\frac{p}{2}+\Psi p-\frac{1}{2}-\Psi \frac{p+1}{2}+\frac{2}{\alpha }\Psi 1-\frac{p-2}{\alpha }-\Psi 1-\frac{2p-2}{\alpha },\end{array}$$with $\mathrm{\Psi }x$ being the digamma function given by the following [24]:$$\begin{array}{}\text{(41)}& \Psi x=\frac{{\Gamma }^{\prime }x}{\Gamma x}=-\kappa +{\displaystyle \sum _{n=1}^{+\infty }\frac{1}{n}-\frac{1}{x+n-1}},\end{array}$$where $\kappa$ being the Euler–Mascheroni constant.

It can further be deduced from (39) that ${\lambda }_{\alpha }p=0$ because ${\mathrm{\Phi }}_{\alpha }p=0$. Applying Taylor series expansion on (40) and exploiting (41), we can rewrite (40) as$$\begin{array}{}\text{(42)}& {\lambda }_{\alpha }p={\displaystyle \sum _{k=0}^{+\infty }{a}_k{p}^k+b_k{p}^{-k-1}}=0,\end{array}$$where$$\begin{array}{c}\text{(43)}& a_k=\frac{{-1}^k}{k!}\begin{array}{c}\frac{2-2^{k+1}}{a^{k+1}}{\Psi }^k\frac{2+\alpha }{\alpha }+\frac{2^k-1}{2^k}{\Psi }^{k}2\\ -{-\frac{1}{2}}^k-{-1}^k{\Psi }^k\frac{3}{2}\end{array},b_k=2{-1}^k-\frac{{-1}^k+1}{2^k},\end{array}$$with ${\mathrm{\Psi }}^k\cdot$ being the polygamma function given by the following [24]:$$\begin{array}{}\text{(44)}& {\Psi }^{k}x={\displaystyle \sum _{n=1}^{+\infty }\frac{{-1}^{k+1}k!}{{n+x-1}^{k+1}}},\end{array}$$$k!$ representing the factorial of $k$.

In (42), $k$ tends to infinity, which is not implementable for numerical calculation. To this end, we can approximate (42) as$$\begin{array}{}\text{(45)}& {\lambda }_{\alpha }p={\displaystyle \sum _{k=0}^K{a}_k{p}^k+b_k{p}^{-k-1}}=0,\end{array}$$where $K$ is sufficiently large. By this, we establish a polynomial whose roots can be determined by using for the numerical calculation. As reported from the previous study [25], $p=\alpha +1/2$ is close to the optimal value and thus can be selected as the initial guess for the numerical calculation to guarantee convergence.

4. Cramér–Rao Bound for Impulsive Noise

In this section, the CRB on the accuracy of estimating the target position is derived for impulsive measurement noise. It should be pointed out that the well-known CRB derived by [9, 10] for Gaussian noise is inapplicable to the target localization problem in impulsive noise. Hence, we focus on deriving a general expression of the CRB for target localization in impulsive noise. It is generally known that the CRB for any unbiased estimator of target position $\mathbf{u}$ is given as follows [26]:$$\begin{array}{}\text{(46)}& \text{CRB}\mathbf{u}={\text{FIM}}^{-1}\mathbf{u},\end{array}$$where $\text{FIM}\mathbf{u}$ is the Fisher information matrix (FIM) given by$$\begin{array}{}\text{(47)}& \text{FIM}\mathbf{u}=E\frac{\partial \ln f\tilde{\mathbf{\psi }}|\mathbf{u}}{\partial \mathbf{u}}{\frac{\partial \ln f\tilde{\mathbf{\psi }}|\mathbf{u}}{\partial \mathbf{u}}}^T\end{array}$$where $f\tilde{\mathbf{\psi }}|\mathbf{u}$ is the PDF of $\tilde{\psi }$ given $\mathbf{u}$. By using the fact that the PDF of $\tilde{\psi }$ satisfies $f_{\tilde{\psi }}\tilde{\mathbf{\psi }}=f_{\epsilon }\tilde{\mathbf{\psi }}-\mathbf{\psi }=f\tilde{\mathbf{\psi }}|\mathbf{x}$, we can rewrite the $\text{FIM}\mathbf{u}$ as$$\begin{array}{}\text{(48)}& \text{FIM}\mathbf{u}=I_c\alpha ,\gamma {\frac{\partial \mathbf{\psi }}{\partial \mathbf{u}}}^{\text{T}}\frac{\partial \mathbf{\psi }}{\partial \mathbf{u}}=I_c\alpha ,\gamma {\frac{\partial \mathbf{\alpha }}{\partial \mathbf{u}}}^{\text{T}}{\mathbf{\Lambda }}^{\text{T}}\mathbf{\Lambda }\frac{\partial \mathbf{\alpha }}{\partial \mathbf{u}},\end{array}$$where $\partial \mathbf{\alpha }/\partial \mathbf{u}$ denotes the partial derivative of $\alpha$ with respect to $\mathbf{u}$. The inner elements of $\partial \mathbf{\alpha }/\partial \mathbf{u}$ have been derived in [9, 10], and interested reader can refer to [9, 10] for details.$$\begin{array}{}\text{(49)}& I_c\alpha ,\gamma ={\displaystyle {\int }_{-\infty }^{+\infty }\frac{{f_{\mathbf{\epsilon }}^{\prime }\epsilon }^2}{f_{\mathbf{\epsilon }}\epsilon }\text{d}\epsilon },\end{array}$$with $f_{\epsilon }\epsilon$ being the PDF of normalized noise and $f_{\epsilon }^{\prime }\epsilon$ the derivative of $f_{\epsilon }\epsilon$. As can be seen from (48) and (49), the noise PDF $f_{\epsilon }\epsilon$ affects the CRB through the scalar multiplier $I_c\alpha ,\gamma$. For $f_{\epsilon }\epsilon$ Gaussian with zero-mean and unit variance, we have $I_c\alpha ,\gamma =1$ and can recover the Gaussian FIM and hence the Gaussian CRB presented in [9, 10]. However, for impulsive noise subject to SαS distribution, as mentioned above, there is no general closed-form PDF. Fortunately, the numerical method in [27] can be adopted to compute the PDF of the α-stable distribution. Hence, the CRB for impulsive noise can be numerically computed using (46), (48), and (49) [28].

5. Simulation Results

This section contains some Monte Carlo simulations to evaluate the performance of the proposed algorithm. The localization scenario is set as follows: a distributed MIMO radar system with $N=6$ geographically separated receive antennas and $M=4$ transmit antennas is deployed to locate a target at position $\mathbf{u}={\mathrm{1000,2000,1000}}^T$m. The positions of the transmit/receive antennas are enumerated in Table 1. $MN=24$ TDs and $N=6$ AOA pairs are extracted to determine the target position. In order to simulate a practical localization scenario, zero-mean SαS distributed noises with specified characteristic exponent parameter $\alpha$ and dispersion parameter $\gamma$ are added to actual true TDs/AOAs. By default, the unit of AOA measurement noise is deg and TD measurement noise is microsecond (μs). The localization accuracy is assessed using root mean squares error (RMSE) and bias defined as follows:$$\begin{array}{c}\text{(50)}& \mathrm{RMSE}\mathbf{u}=\sqrt{\frac{1}{5000}{\displaystyle \sum _{l=1}^{5000}{{\widehat{u}}^l-\mathbf{u}}_2^2}},\text{b}ias\mathbf{u}={\frac{1}{5000}{\displaystyle {\sum }_{l=1}^{5000}{\widehat{u}}^l}-\mathbf{u}}_2,\end{array}$$where ${\widehat{\mathbf{u}}}^l$ is the estimation of $\mathbf{u}$ at the $l\text{th}$ trial, and the number of trials for each case is set as $L=5000$.

Table 1

Positions of the transmit/receive antennas.

Transmitter	TX1	TX2	TX3	TX4	Receiver	RX1	RX2	RX3	RX4	RX5	RX6
xt,m(m)	2000	–2000	2000	–2000	xr,m(m)	4500	–4500	0	6000	–6000	0
yt,m(m)	3000	3000	–3000	–3000	yr,n(m)	4500	–4500	6000	0	0	–6000
zt,m(m)	2000	1000	800	1200	zr,n(m)	2000	1000	2000	1000	1500	1000

The localization RMSEs and biases of the proposed algorithm are evaluated via comparison with existing algorithms including Noroozi’s algorithm (two-stage estimator) in [9], Amiri’s algorithm (one-stage estimator) in [10], Kazemi’s algorithm (Convex estimator) in [11], and the root CRB, under different noise conditions. In order to achieve a more comprehensive insight on the performance of the proposed algorithm, factors including the noise dispersion, noise impulsiveness, target distance, and computation complexity are considered.

5.1. Performance versus Noise Impulsiveness

In the first simulation, the influence of noise characteristic exponent $\alpha$ on localization performance is evaluated. Figure 3 shows the RMSE of the algorithms versus $\alpha \in \mathrm{1,2}$ with noise dispersion parameter $\gamma =1$. The value of $\alpha$ controls the noise impulsiveness. The smaller the value of $\alpha$ is, the more impulsive the noise is. As can be seen from Figure 3(a), under different noise impulsiveness conditions, the RMSE of Amiri’s algorithm is lower than that of Noroozi’s algorithm. This is because Amiri’s algorithm requires only one-stage WLS minimization, and fewer second- and higher-order error terms are discarded. The sensitivity of Kazemi’s algorithm to measurement noise is further reduced. Compared with Noroozi’s algorithm and Amiri’s algorithm, Kazemi’s algorithm achieves smaller RMSE. However, because the above three algorithms are based on the Gaussian noise assumption, the localization RMSEs of the three existing algorithms rise sharply and far larger than the CRB when $\alpha$ is small. This demonstrates that the three existing algorithms are not robust in the presence of impulsive noise and can be considered basically invalid under strong impulsive noise. The proposed algorithm, by contrast, achieves localization RMSEs about 2 to 3 orders of magnitude lower. Under different $\alpha$, the RMSE of the proposed algorithm is basically within 103 m and slightly above the CRLB, which verifies that it is robust to noise impulsiveness. It is worth noting that with the increase in $\alpha$, the RMSE of the proposed algorithm tends to be close to that of the existing algorithms. This is because the noise distribution reduces to Gaussian distribution when $\alpha$ approaches 2.0. Figure 3(b) plots the localization biases of the algorithms. When $\alpha$ approaches 2.0, the biases of the algorithms are comparable and close to zero, which is consistent with the research conclusions in [9–11]. However, when $\alpha$ approaches 1.0, the localization RMSEs of Noroozi’s algorithm, Amiri’s algorithm, and Kazemi’s algorithm are extremely large. The proposed algorithm, by contrast, can still provide approximate unbiased estimate, which reflects the robustness of the proposed algorithm in terms of localization bias.

[figures omitted; refer to PDF]

5.2. Performance versus Noise Dispersion

In the second simulation, we analyze the influence of noise dispersion parameter on the localization performance of the algorithms. Figure 4 presents the RMSE and bias of the target position estimate versus noise dispersion parameter $\gamma$ at $\alpha =1.5$. As mentioned earlier, the dispersion parameter is similar to the standard deviation of Gaussian distribution. It can be seen from Figure 4(a) that with the increase in noise dispersion parameter, the RMSEs of the algorithms increase. As expected, the proposed algorithm has the best performance under different noise dispersion conditions: its RMSE is at least two orders of magnitude lower than that of the existing algorithms and reaches the CRB under different noise dispersion parameters. It can be seen from Figure 4(b) that under different dispersion parameters, there are significant biases in the estimation results of Amiri’s algorithm, Noroozi’s algorithm, and Kazemi’s algorithm. Moreover, due to the nonlinearity of the localization problem, the biases of the three existing algorithms rise sharply with the increase in the noise dispersion parameter. In comparison, the bias of the proposed algorithm is close to zero under different dispersion parameters, which proves the advantage of the proposed algorithm in terms of bias performance.

[figures omitted; refer to PDF]

5.3. Performance versus Target Distance

In the third simulation, we assess the localization performance over the target distance, using the simulation setup as follows: the noise characteristic exponent $\alpha =1.5$, the dispersion parameter $\gamma =1$, the positions of the transmit/receive antennas are listed in Table 1, and the target position is given by $\mathbf{x}=R{\cos {30}^{\mathrm{°}}\cos {45}^{\mathrm{°}},\cos {30}^{\mathrm{°}}\sin {45}^{\mathrm{°}},\sin {30}^{\mathrm{°}}}^{\text{T}}m$ with $R$ varied from 1 km to 100 km. Figure 5 presents the RMSE and bias performance of the algorithms under different target distance $R$. Unsurprisingly, we observe that the proposed algorithm again provides a good localization performance for different target distances, with RMSE close to the CRB and bias about zero. This once again verifies the performance superiority of the proposed algorithm over existing algorithms. Moreover, it can be observed that the localization RMSE and bias performance of the algorithms degrades as the target distance increases. This is consistent with the previous studies [9–11] for the dependence of the localization performance on the geometry of the target and the distributed MIMO radar.

[figures omitted; refer to PDF]

5.4. Performance versus Norm Order p

Next, we evaluate how the variation on norm order $p$ affects the localization performance, using the simulation setup as follows: the noise characteristic exponent $\alpha =1.8$, the dispersion parameter $\gamma =1.8$, the positions of the transmit/receive antennas are listed in Table 1, and the target is located at position $\mathbf{u}={\mathrm{1000,2000,1000}}^T$m. According to the derivation of norm order above, the initial guess $p=1.4$ and the global optimum $p=1.43$. Figure 5 presents the RMSE and bias performance of the algorithms for norm order $p\in \mathrm{1,2}$. As it can be seen from Figure 6, the proposed algorithm achieves the best localization performance at the theoretical global optimum value of $p$, with its RMSE attaining the CRLB. The initial guess of $p$ is very close to the global optimum, which justifies the rationality of choosing $p=\alpha +1/2$ as the initial guess for the numerical calculation of $p$ to guarantee convergence. If the norm order $p$ is too large or too small, the localization performance of the proposed algorithm will degrade. However, by contrast, the performance decline caused by a too small $p$ is relatively not significant. Furthermore, a useful conclusion that can be deduced is that for the impulsive noise with unknown distribution parameters, a small $p$ is recommended to achieve relatively decent localization performance.

[figures omitted; refer to PDF]

5.5. Computation Complexity Comparison

Finally, to evaluate the proposed algorithm in terms of computational complexity, we count the average running time of the algorithms from 5000 independent Monte Carlo runs. The main configuration of the computer is shown as follows: Intel(R) Core(TM) CPU [email protected] GHz; 8.00 G RAM; Windows 10 64 bit Operating System; Matlab 2019a Software. The comparison results are given in Table 2.

Table 2

Time cost of the algorithms.

Algorithms	Average run time (ms)
Noroozi1’s algorithm	1.76
Amiri’s algorithm	0.97
Kazemi’s algorithm	2.36
Proposed algorithm	1.24

As presented in Table 2, the time cost of Noroozi1’s algorithm is almost twice higher than that of Amiri’s algorithm. This is because Noroozi1’s algorithm requires two WLS stages while Amiri’s algorithm determines the target position in only one WLS stage. Kazemi’s algorithm incurs the highest computation complexity among the algorithms. By contrast, the proposed algorithm has the time cost comparable with Amiri’s algorithm. Combining with the performance comparison in Sections 5.1∼5.3, we can conclude that the proposed algorithm significantly improves the target localization performance in impulsive noise, without apparent increase in computation complexity.

6. Conclusions

We have proposed a novel algebraic solution for TD/AOA-based target position estimation distributed MIMO radar. A significant distinction of our study is that the presence of impulsive measurement noise is considered. The proposed algorithm replaces the ${\ell }_2$-norm minimization by the ${\ell }_p$-norm minimization of the error terms and solves the ${\ell }_p$-norm minimization by linearizing the TD and AOA measurement equations and adopting the IRLS approach. Some theoretical analyses, including the optimal choice of the norm order $p$ and the CRB for target position estimation in impulsive noise, are also performed. Simulations results demonstrate that, in impulsive noise, the proposed algorithm outperforms the existing algorithms in terms of localization accuracy and robustness.

Acknowledgments

This study was supported by the Foundation for University Key Teacher by Ministry of Education of Henan Province (No. 2018GGJS234) and the National Natural Science Foundation of China (No. 62003313).

References

[1] S. Sun, A. P. Petropulu, H. V. Poor, "MIMO radar for advanced driver-assistance systems and autonomous driving: advantages and challenges," IEEE Signal Processing Magazine, vol. 37 no. 4, pp. 98-117, DOI: 10.1109/msp.2020.2978507, 2020.

[2] J. Li, P. Stoica, "MIMO radar with colocated antennas," IEEE Signal Processing Magazine, vol. 24 no. 5, pp. 106-114, DOI: 10.1109/msp.2007.904812, 2007.

[3] E. Tohidi, M. Coutino, S. P. Chepuri, H. Behroozi, M. M. Nayebi, G. Leus, "Sparse antenna and pulse placement for colocated MIMO radar," IEEE Transactions on Signal Processing, vol. 67 no. 3, pp. 579-593, DOI: 10.1109/tsp.2018.2881656, 2019.

[4] S. Imani, M. M. Nayebi, S. A. Ghorashi, "Colocated MIMO radar SINR maximization under ISL and PSL constraints," IEEE Signal Processing Letters, vol. 25 no. 3, pp. 422-426, DOI: 10.1109/lsp.2018.2796603, 2018.

[5] A. M. Haimovich, R. S. Blum, L. J. Cimini, "MIMO radar with widely separated antennas," IEEE Signal Processing Magazine, vol. 25 no. 1, pp. 116-129, 2007.

[6] H. Shin, W. Chung, "Distributed MIMO radar target altitude estimation for ground‐based systems," IET Radar, Sonar & Navigation, vol. 13 no. 4, pp. 627-637, DOI: 10.1049/iet-rsn.2018.5159, 2019.

[7] J. Wang, J. Ye, G. Hua, "The optimal detector of distributed MIMO radar under Swerling-Chi scattering models," IEEE Geoscience and Remote Sensing Letters, vol. 17 no. 7, pp. 1129-1133, DOI: 10.1109/lgrs.2019.2943185, 2020.

[8] Y. T. Chan, K. C. Ho, "A simple and efficient estimator for hyperbolic location," IEEE Transactions on Signal Processing, vol. 42 no. 8, pp. 1905-1915, DOI: 10.1109/78.301830, 1994.

[9] A. Noroozi, M. A. Sebt, "Algebraic solution for three‐dimensional TDOA/AOA localisation in multiple‐input-multiple‐output passive radar," IET Radar, Sonar & Navigation, vol. 12 no. 1, pp. 21-29, DOI: 10.1049/iet-rsn.2017.0117, 2018.

[10] R. Amiri, F. Behnia, H. Zamani, "Efficient 3-D positioning using time-delay and AOA measurements in MIMO radar systems," IEEE Communications Letters, vol. 21 no. 12, pp. 2614-2617, DOI: 10.1109/lcomm.2017.2742945, 2017.

[11] S. A. R. Kazemi, R. Amiri, F. Behnia, "Efficient convex solution for 3-D localization in MIMO radars using delay and angle measurements," IEEE Communications Letters, vol. 23 no. 12, pp. 2219-2223, DOI: 10.1109/lcomm.2019.2948175, 2019.

[12] L. Li, T.-S. Qiu, D.-R. Song, "Parameter estimation based on fractional power spectrum under alpha-stable distribution noise environment in wideband bistatic MIMO radar system," AEU-International Journal of Electronics and Communications, vol. 67 no. 11, pp. 947-954, DOI: 10.1016/j.aeue.2013.05.006, 2013.

[13] H. Sadreazami, M. O. Ahmad, M. N. S. Swamy, "A study of multiplicative watermark detection in the contour let domain using alpha-stable distributions," IEEE Transactions on Image Processing, vol. 23 no. 10, pp. 4348-4360, DOI: 10.1109/tip.2014.2339633, 2014.

[14] X. Shen, H. Zhang, Y. Xu, S. Meng, "Observation of alpha-stable noise in the laser gyroscope data," IEEE Sensors Journal, vol. 16 no. 7, pp. 1998-2003, DOI: 10.1109/jsen.2015.2506120, 2016.

[15] G. Zhang, J. Wang, G. Yang, Q. Shao, S. Li, "Nonlinear processing for correlation detection in symmetric alpha-stable noise," IEEE Signal Processing Letters, vol. 25 no. 1, pp. 120-124, DOI: 10.1109/lsp.2017.2776317, 2018.

[16] Y. Chen, J. Chen, "Novel S α S PDF approximations and their applications in wireless signal detection," IEEE Transactions on Wireless Communications, vol. 14 no. 2, pp. 1080-1091, DOI: 10.1109/twc.2014.2364181, 2015.

[17] G. Laguna-Sanchez, M. Lopez-Guerrero, "On the use of alpha-stable distributions in noise modeling for PLC," IEEE Transactions on Power Delivery, vol. 30 no. 4, pp. 1863-1870, DOI: 10.1109/tpwrd.2015.2390134, 2015.

[18] E. E. Kuruoglu, "Nonlinear least ℓp-norm filters for nonlinear autoregressive α -stable processes," Digital Signal Processing, vol. 12 no. 1, pp. 119-142, 2002.

[19] X. Jiang, W.-J. Zeng, H. C. So, S. Rajan, T. Kirubarajan, "Robust matched filtering in ℓ p -Space," IEEE Transactions on Signal Processing, vol. 63 no. 23, pp. 6184-6199, DOI: 10.1109/tsp.2015.2464179, 2015.

[20] P. G. Georgiou, C. Kyriakakis, "Robust maximum likelihood source localization: the case for sub-Gaussian versus Gaussian," IEEE Transactions on Audio Speech and Language Processing, vol. 14 no. 4, pp. 1470-1480, DOI: 10.1109/tsa.2005.860846, 2006.

[21] W. Wen-Jun Zeng, H. C. So, L. Lei Huang, "ℓ p -MUSIC: robust direction-of-arrival estimator for impulsive noise environments," IEEE Transactions on Signal Processing, vol. 61 no. 17, pp. 4296-4308, DOI: 10.1109/tsp.2013.2263502, 2013.

[22] C. S. Burrus, "Iterative re-weighted least squares," Communications on Pure and Applied Mathematics, vol. 44 no. 6, 2009.

[23] Q. Liu, Y. Gu, H. C. So, "DOA estimation in impulsive noise via low-rank matrix approximation and weakly convex optimization," IEEE Transactions on Aerospace and Electronic Systems, vol. 55 no. 6, pp. 3603-3616, DOI: 10.1109/taes.2019.2909728, 2019.

[24] M. Abramowitz, I. A. Stegun, Hand Book of Mathematical Functions, 1964.

[25] Y. Chen, H. C. So, E. E. Kuruoglu, "Variance analysis of unbiased least ℓ p -norm estimator in non-Gaussian noise," Signal Processing, vol. 122, pp. 190-203, DOI: 10.1016/j.sigpro.2015.12.003, 2016.

[26] B. M. Sadler, R. J. Kozick, T. Moore, "Performance analysis for direction finding in non-Gaussian noise," Proceedings of the IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), vol. 5, pp. 2857-2860, DOI: 10.1109/icassp.1999.761358, .

[27] J. P. Nolan, "Numerical calculation of stable densities and distribution functions," Communications in Statistics-Stochastic Models, vol. 13 no. 4, pp. 759-774, DOI: 10.1080/15326349708807450, 1997.

[28] R. J. Kozick, B. M. Sadler, "Maximum-likelihood array processing in non-Gaussian noise with Gaussian mixtures," IEEE Transactions on Signal Processing, vol. 48 no. 12, pp. 3520-3535, 2000.