A Total Lp-Norm Optimization for Bearing-Only Source Localization in Impulsive Noise with SαS Distribution

Abstract

This paper considers the problem of robust bearing-only source localization in impulsive noise with symmetric $\alpha$-stable distribution based on the Lp-norm minimization criterion. The existing Iteratively Reweighted Pseudolinear Least-Squares (IRPLS) method can be used to solve the least LP-norm optimization problem. However, the IRPLS algorithm cannot reduce the bias attributed to the correlation between system matrices and noise vectors. To reduce this kind of bias, a Total Lp-norm Optimization (TLPO) method is proposed by minimizing the errors in all elements of system matrix and data vector based on the minimum dispersion criterion. Subsequently, an equivalent form of TLPO is obtained, and two algorithms are developed to solve the TLPO problem by using Iterative Generalized Eigenvalue Decomposition (IGED) and Generalized Lagrange Multiplier (GLM), respectively. Numerical examples demonstrate the performance advantage of the IGED and GLM algorithms over the IRPLS algorithm.

1. Introduction

Bearing-Only Source Localization (BOSL) using spatially distributed sensors can be widely applied in network localization [1], vehicle localization [2], gunshot localization [3], animal behavior monitoring [4] and rigid body localization [5], to name but a few. The problem of BOSL is to estimate the source location from a set of noise-corrupted bearing measurements where its main challenge originates in the highly nonlinear nature of the angle observations with the true source position. Under the assumption of Gaussian measurement noise, many estimation approaches have been proposed to handle the nonlinearity, including the grid search method, the pseudolinear estimators [6,7], the iterative maximum likelihood methods [8] and the subspace approaches [9]. However, in the complex field environment, the sensor is vulnerable to external interference, enemy attack or node failure. The bearing measurements may suffer impulse noise [10,11,12,13] and those outlier data can degrade the localization performance dramatically. Therefore, it is necessary to develop new estimators that are robust to impulsive noise.

In fact, non-Gaussian models corresponding to impulsive noise have been extensively studied in the literature [13,14,15]. These studies have demonstrated that the Symmetric $\alpha$-Stable ($S\alpha S$) distribution is more suitable to model impulsive noise than the Gaussian distribution. The family of stable distribution is a generalization of the Gaussian distribution under the stable law, including a large range of distributions with mutable values of impulsiveness ($\alpha$), skewness ($\beta$) and dispersion ($\gamma$). In particular, two special cases of $\alpha$-Stable distribution can be obtained by letting impulsiveness parameter $\alpha$ take values of 1 and 2. One is Cauchy distribution ($\alpha =1$), and the other is Gaussian distribution ($\alpha =2$).

In the present study, we focus on the problem of robust BOSL with impulsive noise modeled as $S\alpha S$ distribution. In this situation, the methods derived using L2-norm optimization exhibit unreliable estimates because they are not robust to outliers. Therefore, robust statistics need to be used to improve the localization performance. Under the assumption of $S\alpha S$ measurement noise, a Maximum Likelihood Estimator (MLE) was proposed in [16] to produce optimal estimates. Unfortunately, no closed-form expressions exist for the general likelihood function in the cases of $1<\alpha <2$ and therefore no explicit solution for MLE is possible. Although the $S\alpha S$ distribution does not have finite variance, it has finite Fractional Lower Order Moments (FLOM) [13], which can be calculated by its dispersion $\gamma$ and its characteristic exponent $\alpha$. Moreover, minimizing the FLOM of estimation errors is equivalent to minimizing the dispersion. It is well known as the minimum dispersion criterion, which minimizes the Lp-norm of the estimation residuals. Unlike L2-norm minimization, the least Lp-norm estimator ($1<p<2$) does not have a closed-form solution and consequently needs to be solved in an iterative manner [17].

The least Lp-norm estimator belongs to the M-estimator. The main idea of M-estimate is to replace the sum of squares of least squares residuals with cost functions more robust to impulsive noise, so as to reduce the sensitivity of estimators with respect to model errors. These functions include Huber [18], Bi-square [19], the negative log-likelihood of the Cauchy distribution [20], Wilcoxon [21], L1-norm [22,23], Lp-norm [24] and L∞-norm [25], etc. Ref. [19] presented a distributed robust localization algorithm based on energy information for sensor networks. The algorithm uses Bi-square function as the cost function of M-estimate. A distributed incremental least mean square algorithm based on Wilcoxon norm was proposed in [21] for parameter estimation of sensor networks. Ref. [26] proposed a robust structure total least-squares algorithm for BOSL by using the improved Danish weight function to suppress the impact of outlier data on the localization performance. In addition to the M-estimator, there are other algorithms that can handle outlier data, such as outlier detection [27], clustering [28] and game theoretic techniques [29]. The outlier detection method [27] is to detect suspected outlier data first, and separate it from the original data set. The clustering based techniques [28] can be used to classify normal and abnormal sensors. In game theory [29], the defense strategies could be employed to detect the outlier data and adaptive threshold selection can be also introduced.

Recently, an Iteratively Reweighted Pseudolinear Least-Squares (IRPLS) method was proposed in [14] to reduce biases attributed to the impulsive noise. However, IRPLS still suffers from a major bias problem caused by the correlation between system matrix and the noise vector. This bias can be reduced by exploiting the use of an Instrumental-Variable (IV) matrix [14] that is approximately uncorrelated with the noise vector. Inconsistent with the IV method, we present a robust total least-squares method using Lp-norm minimization that can reduce biases by minimizing the errors in the system matrix and the data vector simultaneously. We first formulate the problem of BOSL subjected to impulsive noise modeled as $S\alpha S$ distribution and review the pseudo-linear measurement model for BOSL. Next, we present the Total LP-norm Optimization (TLPO) method for BOSL to minimize the the errors in the system matrix and the data vector under the minimum dispersion criterion. Two algorithms, named Iterative Generalized Eigenvalue Decomposition (IGED) and Generalized Lagrange Multiplier (GLM), are designed to solve the TLPO problem. The performance advantage of the proposed algorithms is demonstrated by numerical simulations in terms of both bias and Root-Mean-Square-Error (RMSE) performance. The main contributions of the proposed method can be summarized as follows:

• Development of a new bias reduced estimator based on TLPO for BOSL when the measurement noise is modeled as $S\alpha S$ distribution;
• Development of two algorithms for TLPO optimization using the IGED approach and the GLM method, respectively.

The rest of this paper is organized below. Section 2 briefly summarizes the $S\alpha S$ distribution, presents the measurement model for BOSL and discusses the nonlinear least Lp-norm for BOSL. In Section 3, the pseudolinear estimator and the iteratively reweighted pseudolinear least-squares algorithm are reviewed. In Section 4, a new TLPO method is presented, and two iterative algorithms are developed to solve the TLPO problem. Section 5 presents the Cramér-Rao Lower Bound and Section 6 illustrates the bias and RMSE performance of PLE, TLS, LAR and TLAR using various numerical examples. Lastly, conclusions are drawn in Section 7.

2. Lp-Norm Optimization for Robust BOSL

2.1. Symmetric Alpha-Stable Distribution

The impulsive noise is more likely to exhibit outliers than normal noise. Studies [13] have shown that $S\alpha S$ distribution can model the impulsive noise better than Gaussian distribution due to the reason that the $S\alpha S$ densities have heavier tails than the Gaussian density. The characteristic function of $S\alpha S$ distribution is described as:

$$\chi \left(\tau \right)=exp\{i\delta \tau -\gamma |\tau {|}^{\alpha }\}$$

(1)

where $\alpha$ ($0<\alpha \le 2$) denotes the characteristic exponent, $\gamma$ indicates the dispersion parameter and $\delta$ stands for the location parameter. The $S\alpha S$ distribution is completely determined by these three parameters. The value of $\alpha$ indicates the heaviness of the tails of the density function. A small positive value of $\alpha$ implies high impulsiveness, while a value of $\alpha$ close to 2 shows a type of Gaussian-like shape. The dispersion $\gamma$ performs like the variance and $\delta$ controls the location of the density function.

A $S\alpha S$ distribution is called standard if $\delta =0$, $\gamma =1$. Let x be a random variable that follows the standard $S\alpha S$ distribution. By taking the inverse Fourier transform of its characteristic function, the density function of x is of the form:

$$f_{\alpha }\left(x\right)=\frac{1}{2\pi }{\int }_{-\infty }^{+\infty }exp\left(-ix\tau -{\gamma \left|\tau \right|}^{\alpha }\right)d\tau$$

(2)

It is known that a $S\alpha S$ distribution with characteristic exponent $\alpha$ only has finite moments for orders less than $\alpha$, which are called the Fractional Lower Order Moments (FLOM). In particular, the FLOM of a standard $S\alpha S$ random variable x is given by

$${\mathbb{E}\left\{\right|x|}^p\}=C(p,\alpha ){\gamma }^{p/\alpha },0<p<\alpha$$

(3)

where $\mathbb{E}\{·\}$ denotes expectation operator,

$$C(p,\alpha )=\frac{2^{p+1}\mathsf{\Gamma }\left(\frac{p+1}{2}\right)\mathsf{\Gamma }\left(-\frac{p}{\alpha }\right)}{\alpha \sqrt{\mathsf{\pi }}\mathsf{\Gamma }\left(-\frac{p}{2}\right)}$$

(4)

is a constant that depends only on q and $\alpha$ and $\mathsf{\Gamma }\left(a\right)={\int }_0^{+\infty }{y}^{(a-1)}exp(-y)dy$ is the Gamma function. It is worth mentioning that the linear space of a $S\alpha S$ process is a Banach space for $\alpha \in [1,2)$, and it is only a metric space for $\alpha \in (0,1)$ [30]. Therefore, the tools of Hilbert space are not applicable when one solves a linear estimation problem with the $S\alpha S$ process. In the present study, we only focus on the case of $\alpha \in (1,2)$.

2.2. Measurement Model

The problem of robust BOSL is depicted in Figure 1, where $\mathit{t}={[t_x,t_y]}^T$ denotes the unknown target location vector, ${\mathit{r}}_{\mathit{m}}={[r_{x,m},r_{y,m}]}^T$ represents the mth sensor location vector and ${\theta }_m$ is the true bearing at sensor m. The nonlinear relationship between ${\theta }_m$, $\mathit{t}$ and ${\mathit{r}}_{\mathit{m}}$ is given by [31]

$${\theta }_m={tan}^{-1}(t_y-r_{y,m},t_x-r_{x,m}),{\theta }_m\in [0,2\pi )$$

(5)

where ${tan}^{-1}$ denotes the two-argument inverse tangent function and $m=1,2,\dots ,M$.

The bearing measurement taken at sensor m is given by

$${\tilde{\theta }}_m={\theta }_m+n_m$$

(6)

where $n_m$ follows independent zero mean $S\alpha S$ distribution with ${\gamma }_m$ representing the dispersion parameter and $\alpha$ accounting for the characteristic exponent. The dispersion ${\gamma }_m$ can vary with m. If the noise dispersion ${\gamma }_m$ is known a priori, then the noise term $n_m$ can be normalized. Let $e_m={\gamma }_m^{-1/\alpha }{n}_m$ denote the normalized noise that has unit dispersion. The normalized measurement equation can be written as

$${\tilde{\varphi }}_m={\varphi }_m+e_m$$

(7)

where ${\tilde{\varphi }}_m={\gamma }_m^{-1/\alpha }{\tilde{\theta }}_m$ is the normalized bearing measurement and ${\varphi }_m={\gamma }_m^{-1/\alpha }{\theta }_m$ is the normalized true bearing. Stacking (7) in the vector form yields

$$\tilde{\mathit{\varphi }}=\mathit{\varphi }+\mathit{e}$$

(8)

where $\tilde{\mathit{\varphi }}={[{\tilde{\varphi }}_1,{\tilde{\varphi }}_2,\dots ,{\tilde{\varphi }}_M]}^T$ denotes the normalized measurement vector, $\mathit{\varphi }={[{\varphi }_1,{\varphi }_2,\dots ,{\varphi }_M]}^T$ represents the true bearing vector and $\mathit{e}={[e_1,e_2,\dots ,e_M]}^T$ indicates the measurement noise vector.

2.3. Minimum Dispersion Criterion

The main difficulty for parameter estimation with the $\alpha$-stable process is that all non-Gaussian $\alpha$-stable distributions have infinite variance. As such, the traditional nonlinear least-squares techniques [32], which rely on the second order moments, are not suitable for solving BOSL problems with impulsive noise. Fortunately, we can use the minimum dispersion criterion instead of minimizing the variance. To do this, we first define the norm of the $S\alpha S$ random variable x as

$${\parallel x\parallel }_{\alpha }={\gamma }^{\frac{1}{\alpha }},\alpha \in [1,2)$$

(9)

Thus, a suitable measure of the dispersion could be $\gamma ={\parallel x\parallel }_{\alpha }^{\alpha }$. Combining (3) and (9), one can easily obtain

$${\left(\frac{\mathbb{E}\left\{\right|x{|}^p\}}{C(p,\alpha )}\right)}^{\frac{1}{p}}\equiv {\parallel x\parallel }_{\alpha },p\in (0,\alpha ),\alpha \in [1,2)$$

(10)

Therefore, we may estimate the source location by solving the following nonlinear Lp-norm optimization problem:

$$\widehat{\mathit{t}}=arg\underset{\mathit{t}}{min}J\left(\mathit{t}\right),J\left(\mathit{t}\right)=\parallel \tilde{\mathit{\varphi }}{-\mathit{\varphi }\parallel }_p^p=\sum _{m=1}^M{|{\tilde{\varphi }}_m-{\varphi }_m\left(\mathit{t}\right)|}^p$$

(11)

where $p\in (1,\alpha )$. The definition of the Lp-norm of the vector $\mathit{\zeta }$ is

$${\parallel \mathit{\zeta }\parallel }_p={\left(\sum _k{\left|{\zeta }_k\right|}^p\right)}^{1/p}$$

(12)

In the present study, we do not consider the case of $0<p\le 1$, since the corresponding Lp-norm cost function is not differentiable. For non-normalized bearing measurements, the Lp-norm objective function in (11) becomes

$$J\left(\mathit{t}\right)=\sum _{m=1}^M{\left|{\gamma }_m^{-1/\alpha }({\tilde{\theta }}_m-{\theta }_m\left(\mathit{t}\right))\right|}^p$$

(13)

The optimization problem listed in (11) can be solved numerically for a given 2D space of interest. To be more specific, we can perform a grid search over that 2D space. The global solution of maximum likelihood (ML) estimator is guaranteed for the given set of data as long as the spacing between grids is small enough. However, if the range of the parameter of interest is not confined to a relatively small interval or the dimension of unknown parameter vector is high, the grid search approach is computationally infeasible. Instead, one may resort to iterative optimization methods, such as gradient decent, and Gauss–Newton [33], etc. In particular, the Gauss–Newton algorithm is

$${\widehat{\mathit{t}}}^{(i+1)}={\widehat{\mathit{t}}}^{\left(i\right)}-{\left(\nabla J{\left(\mathit{t}\right)}^T\nabla J\left(\mathit{t}\right)\right)}^{-1}\nabla J{\left(\mathit{t}\right)}^T\left(\tilde{\mathit{\varphi }}-\mathit{\varphi }\left({\widehat{\mathit{t}}}^{\left(i\right)}\right)\right)$$

(14)

where $\nabla J\left(\mathit{t}\right)$ denotes the Jacobian matrix,

$$\nabla J\left(\mathit{t}\right)=p\sum _{m=1}^M\frac{\partial {\varphi }_m\left(\mathit{t}\right)}{\partial \mathit{t}}{\left|e_m\right|}^{p-1}\mathrm{sign}\left(e_m\right)$$

(15)

If $e_m>0$, $\mathrm{sign}\left(e_m\right)=1$ and $-1$ otherwise. The statistical properties of the nonlinear Lp-norm minimizer have been studied in the literature [24]. The theoretical covariance $\mathit{B}$ is related to the value of p, and it is given by [24]:

$$\mathit{B}\approx \frac{C(2p-2,\alpha )}{{(p-1)}^2{C}^2(p-2,\alpha )}{(\nabla {\mathit{\varphi }}^T\nabla \mathit{\varphi })}^{-1}$$

(16)

where $\nabla \mathit{\varphi }$ denotes the Jacobian matrix of $\mathit{\varphi }$ with respect to $\mathit{t}$. By minimizing the scalar term of $\mathit{B}$, we obtain the optimal choice of p

$$p^o=arg\underset{p}{min}\frac{C(2p-2,\alpha )}{{(p-1)}^2{C}^2(p-2,\alpha )}$$

(17)

where $p^o$ denotes the optimal value of p.

3. Pseudolinear Lp-Norm Minimization

3.1. Pseudolinear Estimator

The robust BOSL problem is nontrivial because the measurement Equation in (7) is nonlinearly related to the unknown source location. An attractive solution is to set up a pseudolinear equation by lumping the nonlinearities into the noise term. As illustrated in Figure 1, an orthogonal vector sum between the measured angle vector and the true angle vector can be geometrically described by

$$\begin{array}{c}\hfill {\mathit{u}}_m={\tilde{\mathit{u}}}_m+{\epsilon }_m=\mathit{t}-{\mathit{r}}_m\end{array}$$

(18)

where ${\mathit{u}}_m$ is the true angle vector between ${\mathit{r}}_m$ and $\mathit{t}$, ${\tilde{\mathit{u}}}_m$ is the measured angle vector starting from ${\mathit{r}}_m$ and produces the noisy bearing ${\tilde{\theta }}_m$ with respect to the horizontal direction and ${\epsilon }_m$ is the error vector. Let ${\mu }_m={[cos{\tilde{\theta }}_m,sin{\tilde{\theta }}_m]}^T$ and ${\nu }_m={[sin{\tilde{\theta }}_m,-cos{\tilde{\theta }}_m]}^T$ denote two orthogonal unit trigonometric vectors. Then, ${\tilde{\mathit{u}}}_m$ and ${\epsilon }_m$ are given by

$$\begin{array}{cc}\hfill {\tilde{\mathit{u}}}_m& =\parallel {\mathit{u}}_m\parallel cos{n}_m·{\mathit{\mu }}_m\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill {\mathit{\epsilon }}_m& =\parallel {\mathit{u}}_m\parallel sin{n}_m·{\mathit{\nu }}_m\hfill \end{array}$$

(20)

where $\parallel ·\parallel$ denotes Euclidean norm. Note from (19) that ${\tilde{\mathit{u}}}_m^T{\mathit{\nu }}_m=0$. Substituting (19) and (20) into (18) and multiplying (18) with ${\nu }_m^T$ yields

$$\begin{array}{c}\hfill {\xi }_m={\nu }_m^T\mathit{t}-{\nu }_m^T{\mathit{r}}_m\end{array}$$

(21)

where ${\xi }_m=\parallel {\mathit{u}}_m\parallel sin{n}_m$ is a nonlinear transformed measurement error. Collecting the pseudolinear equation errors as a vector $\mathit{\xi }={[{\xi }_1,{\xi }_2,\dots ,{\xi }_K]}^T$, we obtain

$$\begin{array}{c}\hfill \mathit{\xi }=\mathit{A}\mathit{t}-\mathit{b}\end{array}$$

(22)

where $\mathit{A}={[{\mathit{\nu }}_1,{\mathit{\nu }}_2,\dots ,{\mathit{\nu }}_M]}^T$, $\mathit{b}={[{\mathit{\nu }}_1^T{\mathit{r}}_1,\dots ,{\mathit{\nu }}_M^T{\mathit{r}}_M]}^T$ are the measurement matrix and vector, respectively.

The PLE requires that $\mathit{t}$ be estimated by minimizing ${\parallel \mathit{A}\mathit{t}-\mathit{b}\parallel }_2^2$ with respect to $\mathit{t}$ in an L2-norm optimization sense, and the solution is

$$\begin{array}{c}\hfill {\widehat{\mathit{t}}}_{\mathrm{PLE}}={\left({\mathit{A}}^T\mathit{A}\right)}^{-1}{\mathit{A}}^T\mathit{b}\end{array}$$

(23)

The above PLE solution has a large bias because the abnormal nodes can produce significantly large bearing measurement errors. As a typical robust estimation algorithm, Lp-norm minimization has been widely used for parameter estimation using measurements corrupted by impulse noise. In this paper, we aim to develop robust estimators using Lp-norm optimization to reduce the bias as much as possible.

3.2. Iteratively Reweighted Pseudolinear Least-Squares Algorithm

This subsection reviews an Iteratively Reweighted Pseudolinear Least-Squares (IRPLS) algorithm, which has been proposed in [14]. The IRPLS algorithm is derived from the following Lp-norm optimization problem:

$$\underset{\mathit{t}}{min}{\parallel \mathit{G}(\mathit{A}\mathit{t}-\mathit{b})\parallel }_p^p=\underset{\mathit{t}}{min}\sum _{m=1}^M{\left|{\mathit{g}}_m({\left(\mathit{A}\mathit{t}\right)}_m-b_m)\right|}^p$$

(24)

where $\mathit{G}=\mathrm{diag}(g_1,g_2,\dots ,g_M)$, $g_m=({\gamma }_m^{1/\alpha }\parallel \mathit{t}-{\mathit{r}}_m{\parallel )}^{-1}$, ${\left(\mathit{A}\mathit{t}\right)}_m$ and $b_m$ are the m-th element of $\mathit{A}\mathit{t}$ and $\mathit{b}$, respectively. When the measurement noise is small, we have $sin{n}_m\approx n_m$. Then, $|{\mathit{g}}_m({\left(\mathit{A}\mathit{t}\right)}_m-b_m){|}^p$ can be expressed approximately as

$$|{\mathit{g}}_m({\left(\mathit{A}\mathit{t}\right)}_m-b_m){|}^p\approx {\gamma }_m^{-p/\alpha }\parallel \mathit{t}-{\mathit{r}}_m{\parallel }^{-2}|{\tilde{\theta }}_m-{\theta }_m{|}^{(p-2)}{|{\left(\mathit{A}\mathit{p}\right)}_m-b_m|}^2$$

(25)

The above interpretation suggests an iterative reweighted pseudolinear least-squares with the weight matrix for the i-th iteration given by

$${\mathit{W}}^{\left(i\right)}=\mathrm{diag}\left({\gamma }_1^{-\frac{p}{2\alpha }}\parallel {\widehat{\mathit{t}}}^{(i-1)}-{\mathit{r}}_1{\parallel }^{-1}{\left|{\tilde{\theta }}_1-{\widehat{\theta }}_1^{(i-1)}\right|}^{\frac{p-2}{2}},\dots ,{\gamma }_M^{-\frac{p}{2\alpha }}{\parallel {\widehat{\mathit{t}}}^{(i-1)}-{\mathit{r}}_M\parallel }^{-1}{\left|{\tilde{\theta }}_M-{\widehat{\theta }}_M^{(i-1)}\right|}^{\frac{p-2}{2}}\right)$$

(26)

where ${\widehat{\mathit{t}}}^{(i-1)}$ and ${\widehat{\theta }}_m^{(i-1)}$ denote the estimated source location and the mth bearing obtained from the previous iteration $i-1$. Using the weight matrix, the weight error for the k-th iteration is given by

$${\left({\mathit{W}}^{\left(i\right)}{\mathit{\xi }}^{\left(i\right)}\right)}^T{\mathit{W}}^{\left(i\right)}{\mathit{\xi }}^{\left(i\right)}=\sum _{m=1}^M{\left|{\mathit{g}}_m({\left(\mathit{A}\mathit{t}\right)}_m-b_m)\right|}^p$$

(27)

The homotopy method can be applied in the iterations by starting with a value for p equals to 2 (the weight matrix is unit matrix) and decreasing it each iteration until it reaches the designated value. Thus,

$$p^{\left(i\right)}=max(p,\kappa p^{(i-1)})$$

(28)

where $\kappa <1$ is the homotopy parameter. The source location can be estimated by performing least-squares estimation:

$${\widehat{\mathit{t}}}^{\left(i\right)}={\left({\left({\mathit{W}}^{\left(i\right)}\mathit{A}\right)}^T{\mathit{W}}^{\left(i\right)}\mathit{A}\right)}^{-1}{\left({\mathit{W}}^{\left(i\right)}\mathit{A}\right)}^T{\mathit{W}}^{\left(i\right)}\mathit{b}$$

(29)

The IRPLS algorithm is stopped when $\parallel {\widehat{\mathit{t}}}^{(i-1)}-{\widehat{\mathit{t}}}^{\left(i\right)}\parallel \le ϵ$, where $ϵ$ is a tolerance parameter.

Finally, the whole process of the IRPLS algorithm is summarized in Algorithm 1.

Algorithm 1 The IRPLS algorithm.

1: Initialization: Set ${\widehat{\mathit{t}}}^{\left(0\right)}={\widehat{\mathit{t}}}_{\mathrm{PLE}}$. 2: for$i=0$; $i++$do 3:  Compute the weight matrix using (31). 4:  Determine the value of p for each iteration using (28). 5:  Perform weighted least squares estimation using (29). 6:  if $\parallel {\widehat{\mathit{t}}}^{(i-1)}-{\widehat{\mathit{t}}}^{\left(i\right)}\parallel \le ϵ$, we obtain the final solution, and stop the iteration. Otherwise $i=i+1$, go to step 3. 7: end for

The Lp-norm criterion ensures that the cost function in (24) gives less weight to large deviations, and therefore reduces the bias formed by pseudolinear errors with large residuals. However, the IRPLS algorithm implicitly assumes that only $\mathit{b}$ is subjected to errors. This is not the case, since the system matrix $\mathit{A}$ is corrupted with measurement noises as well. The correlation between $\mathit{A}$ and $\mathit{b}$ causes the IRPLS estimator to be inconsistent. In this subsection, we analyze such bias attributed to the correlation between $\mathit{A}$ and $\mathit{b}$.

Let ${\widehat{\mathit{t}}}^{*}$ denote the final solution of the IRPLS algorithm. This solution satisfies

$${\widehat{\mathit{t}}}^{*}={\left({\mathit{A}}^T\mathcal{W}\left({\widehat{\mathit{t}}}^{*}\right)\mathit{A}\right)}^{-1}{\mathit{A}}^T\mathcal{W}\left({\widehat{\mathit{t}}}^{*}\right)\mathit{b}$$

(30)

where

$$\mathcal{W}\left({\widehat{\mathit{t}}}^{*}\right)=\mathrm{diag}\left({\gamma }_1^{-\frac{p}{\alpha }}\parallel {\widehat{\mathit{t}}}^{*}-{\mathit{r}}_1{\parallel }^{-2}{\left|{\tilde{\theta }}_1-{\widehat{\theta }}_1^{*}\right|}^{(p-2)},\dots ,{\gamma }_M^{-\frac{p}{\alpha }}{\parallel {\widehat{\mathit{t}}}^{*}-{\mathit{r}}_M\parallel }^{-2}{\left|{\tilde{\theta }}_M-{\widehat{\theta }}_M^{*}\right|}^{(p-2)}\right)$$

(31)

and ${\widehat{\theta }}_m^{*}$ denotes the mth bearing obtained from (5) by using ${\widehat{\mathit{t}}}^{*}$.

According to (22) and (30), we can get

$$\Delta \mathit{t}={\widehat{\mathit{t}}}^{*}-\mathit{t}=-{\left({\mathit{A}}^T\mathcal{W}\left({\widehat{\mathit{t}}}^{*}\right)\mathit{A}\right)}^{-1}{\mathit{A}}^T\mathcal{W}\left({\widehat{\mathit{t}}}^{*}\right)\mathit{\xi }$$

(32)

The analytical bias of IRPLS is given by [14]

$$\mathbb{E}\{\Delta \mathit{t}\}=-\frac{1}{p-1}\mathbb{E}{\left\{\frac{{\mathit{A}}^T\mathcal{W}\left(\mathit{t}\right)\mathit{A}}{M}\right\}}^{-1}\mathbb{E}\left\{\frac{{\mathit{A}}^T\mathcal{W}\left(\mathit{t}\right)\mathit{\xi }}{M}\right\}$$

(33)

where the weighting matrix $\mathcal{W}\left(\mathit{t}\right)$ is computed from (31) using true source location parameter $\mathit{t}$.

4. Total Lp-Norm Optimization

In this section, we present a method of TLPO for BOSL. The IRPLS algorithm only considers the deviation of the system vector $\mathit{b}$. But in fact, the system matrix $\mathit{A}$ also has a measurement residual. The IRPLS algorithm will inevitably cause a large bias, because the system matrix $\mathit{A}$ and the system vector $\mathit{b}$ are statistically correlated.

4.1. Method Description

In order to improve the performance of the IRPLS algorithm and reduce the bias caused by the correlation between $\mathit{A}$ and $\mathit{b}$, we develop the TLPO method in this subsection. Unlike the IRPLS method, the TLPO method uses the correction matrix $\Delta \mathit{A}$ and the correction vector $\Delta \mathit{b}$ to compensate the system matrix $\mathit{A}$ and system vector $\mathit{b}$, respectively. The normalized equation can be written as

$$\mathit{G}(\mathit{A}+\Delta \mathit{A})\mathit{t}=\mathit{G}(\mathit{b}+\Delta \mathit{b})$$

(34)

where $\Delta \mathit{A}$ is the perturbation matrix of $\mathit{A}$ and $\Delta \mathit{b}$ is the perturbation vector of $\mathit{b}$. Let $\mathit{K}=\mathit{G}[\mathit{A},\mathit{b}]$ and $\mathit{v}={[\mathit{t},-1]}^T$ be the augmented matrix and vector, and Equation (34) can be rewritten as

$$(\mathit{K}+\mathbf{\Omega })\mathit{v}=0$$

(35)

where $\mathbf{\Omega }=\mathit{G}[\Delta \mathit{A},\Delta \mathit{b}]$ is the error matrix.

In order to reduce the bias of IRPLS, we use the TLPO approach to minimize the perturbation matrix $\Delta \mathit{A}$ and perturbation vector $\Delta \mathit{b}$ simultaneously, and the TLPO problem for BOTL can be formulated as

$$\underset{\mathit{E},\mathit{t}}{min}{\parallel \mathbf{\Omega }\parallel }_p,\mathrm{s}.\mathrm{t}.(\mathit{K}+\mathbf{\Omega })\mathit{v}=0$$

(36)

for $1<p<2$. Note that if $p=1$, (36) becomes the total least absolute residual method, and if $p=2$, it is the well-known TLS method. To solve the TLPO problem conveniently, it is necessary to develop an equivalent form of (36). The following proposition holds.

Proposition 1.

Given the TLPO problem defined in (36), the estimation of $\mathit{t}$ from the minimization of ${\parallel \mathbf{\Omega }\parallel }_p$ subject to $(\mathit{K}+\mathbf{\Omega })\mathit{v}=0$ is equivalent to

$$\underset{\mathit{v}}{min}{\parallel \mathit{K}\mathit{v}\parallel }_p,\mathrm{s}.\mathrm{t}.{\parallel \mathit{v}\parallel }_q-1=0$$

(37)

where $\mathit{v}={[\mathit{t},-1]}^T$, $\frac{1}{p}+\frac{1}{q}=1$ and $1<p<\alpha$.

Proof. 

From the constraint $(\mathit{K}+\mathbf{\Omega })\mathit{v}=0$, we have ${\parallel \mathit{K}\mathit{v}\parallel }_p={\parallel \mathbf{\Omega }\mathit{v}\parallel }_p$. Let ${\omega }_m^T$ denote the mth row of $\mathbf{\Omega }$. ${\parallel \mathbf{\Omega }\mathit{v}\parallel }_p^p$ becomes

$$\begin{array}{cc}\hfill {\parallel \mathbf{\Omega }\mathit{v}\parallel }_p^p& =\sum _{m=1}^M{\left|{\mathit{\omega }}_m^T\mathit{v}\right|}^p\hfill \\ \hfill & \le \sum _{m=1}^M{∥{\mathit{\omega }}_m∥}_p^p{∥\mathit{v}∥}_q^p\hfill \end{array}$$

(38)

for p and q satisfy the equation $\frac{1}{p}+\frac{1}{q}=1$ and $1<p<\alpha$. To derive (38), we used the properties of hölder’s inequality [34]. If the Lq-norm of $\mathit{v}$ satisfies ${\parallel \mathit{v}\parallel }_q=1$, then the inequality (38) becomes ${\parallel \mathbf{\Omega }\mathit{v}\parallel }_p\le {\parallel \mathbf{\Omega }\parallel }_p$. Therefore, we can conclude that the minimization of ${\parallel \mathbf{\Omega }\parallel }_p$ subject to $(\mathit{K}+\mathbf{\Omega })\mathit{v}=0$ is equivalent to (37). □

Proposition 1 provides a facilitated way to solve the TLPO problem. In the following two subsections, we concentrate on deriving the solution on this particular problem and developing the IGED and GLM algorithms.

4.2. The IGED Algorithm

The optimization problem in Proposition 1 is equivalent to

$$\underset{\mathit{v}}{min}{\parallel \mathit{K}\mathit{v}\parallel }_p^p,\mathrm{s}.\mathrm{t}.{\parallel \mathit{v}\parallel }_q^q-1=0$$

(39)

To solve (39), we first transform it into an unconstrained optimization problem by using the Lagrange multiplier method. The appropriate Lagrangian function of (39) can be written as

$$\begin{array}{cc}\hfill \mathcal{L}(\mathit{v},\mathit{\lambda })& ={\parallel \mathit{K}\mathit{v}\parallel }_p^p{+\lambda (1-\parallel \mathit{v}\parallel }_q^q)\hfill \\ \hfill & ={\mathit{v}}^T{\mathit{K}}^T\mathit{D}\mathit{K}\mathit{v}+\lambda (1-{\mathit{v}}^T\mathit{C}\mathit{v})\hfill \end{array}$$

(40)

where $\mathit{D}=\mathrm{diag}\left(\right|{\left(\mathit{K}\mathit{v}\right)}_m{|}^{p-2})$, ${\left(\mathit{K}\mathit{v}\right)}_m$ represents the mth element of vector $\mathit{K}\mathit{v}$ for $m=1,2,\dots ,M$, $\mathit{C}=\mathrm{diag}\left(\right|v_1^{q-2}|,|v_2^{q-2}|,|v_3^{q-2}\left|\right)$ and $\lambda$ is the Lagrangian multiplier. Taking the partial derivative of $\mathcal{L}(\mathit{v},\lambda )$ with respect to $\mathit{v}$ and setting it to zero yields

$${\mathit{K}}^T\mathit{D}\mathit{K}\mathit{v}=\lambda \mathit{C}\mathit{v}$$

(41)

After multiplying both sides of (41) with ${\mathit{v}}^T$, it can be found that only $\lambda ={\mathit{v}}^T{\mathit{K}}^T\mathit{D}\mathit{K}\mathit{v}$ is the cost to be minimized. So, the optimal solution ${\mathit{v}}^{*}$ is the eigenvector of the smallest generalized eigenvalue of $({\mathit{K}}^T\mathit{D}\mathit{K},\mathit{C})$.

Figure 2 gives a flowchart to depict the whole process of the IGED algorithm. The algorithm starts from the initialization procedure, where ${\mathit{v}}_0$ is set as ${[{\widehat{\mathit{t}}}_{\mathrm{PLE}}^T,-1]}^T$ and ${\widehat{\mathit{t}}}_{\mathrm{PLE}}$ is obtained from (23). Given p and q, matrices $\mathit{C}$ and $\mathit{D}$ can be computed, respectively. After performing generalized eigenvalue decomposition from the pair $({\mathit{K}}^T\mathit{D}\mathit{K},\mathit{C})$, the solution of $\mathit{v}$ in (41) is the generalized eigenvector of $({\mathit{K}}^T\mathit{D}\mathit{K},\mathit{C})$ that gives the minimum generalized eigenvalue. Since IGED is an iterative algorithm, it requires the convergence check. Let ${\lambda }_i$ and ${\lambda }_{i+1}$ be the minimum generalized eigenvalues at the ith and $i+1$th iterations, respectively. If there exists $\epsilon$ such that

$$\left|\frac{{\lambda }_{i+1}-{\lambda }_i}{{\lambda }_{i+1}}\right|\le \epsilon$$

(42)

then we would get the final solution ${\mathit{v}}^{*}$ and the estimate for source location is

$$\widehat{\mathit{t}}=\frac{{\mathit{v}}^{*}(1:2)}{{\mathit{v}}^{*}\left(3\right)}.$$

(43)

In summary, the IGED Algorithm 2 can be implemented by the following procedures.

Algorithm 2 The IGED algorithm.

1: Initialization:${\mathit{v}}_0={[{\widehat{\mathit{t}}}_{\mathrm{PLE}}^T,-1]}^T$. 2: for$i=0$; $i++$do 3:  Compute $\mathit{C}$ and $\mathit{D}$. 4:  Perform generalized eigenvalue decomposition of the pair $({\mathit{K}}^T\mathit{D}\mathit{K},\mathit{C})$. 5:  Select the smallest generalized eigenvalue and its corresponding generalized eigenvector. 6:  Set ${\lambda }_i$ and ${\lambda }_{i+1}$ as the minimum generalized eigenvalues for iterations i and $i+1$. If $|({\lambda }_{i+1}-{\lambda }_i)/{\lambda }_{i+1}|>\epsilon$, go to step 3. Otherwise, obtain the final solution ${\mathit{v}}^{*}$ and stop. 7:  Get the source location estimate $\widehat{\mathit{t}}$ using (43). 8: end for

4.3. The GLM Algorithm

The GLM method [35] is widely used for constrained optimization. Its basic principle is to add a penalty term to the Lagrangian function to form an augmented Lagrangian function, which can impose a larger penalty on infeasible points. Thus, the constrained optimization problem (37) would be transformed into a new unconstrained optimization problem that can be solved efficiently.

To solve (37) with the GLM method, we formulate an augmented Lagrangian function as follows:

$$\mathcal{L}(\mathit{v},\lambda ,s)=g\left(\mathit{v}\right)+\frac{s}{2}{\parallel h\left(\mathit{v}\right)\parallel }_2^2-\lambda h\left(\mathit{v}\right)$$

(44)

where $g\left(\mathit{v}\right)={\parallel \mathit{K}\mathit{v}\parallel }_p$ and $h\left(\mathit{v}\right)=1-{\parallel \mathit{v}\parallel }_q$. Compared with the Lagrangian function, the GLM cost function has an additional term $s/{2\parallel h\left(\mathit{v}\right)\parallel }_2^2$. This item is a punishment for violation of constraint $h\left(\mathit{v}\right)=0$. The punishment parameter s determines the degree of punishment, and it is generally sufficiently large. $\lambda$ is a positive scalar, and the term $\lambda h\left(\mathit{v}\right)$ is to ensure that the optimal solution is a strict local minimum point of $F\left(\mathit{v}\right)=g\left(\mathit{v}\right)+\frac{s}{2}{\parallel h\left(\mathit{v}\right)\parallel }_2^2$ under the condition of obeying the constraint $h\left(\mathit{v}\right)=0$.

Figure 3 shows the flow diagram of the GLM algorithm. For the initialization step, we set ${\mathit{v}}_0={[{\widehat{\mathit{t}}}_{\mathrm{PLE}}^T,-1]}^T$. Given initial multiplier $\lambda$, penalty factor s, admissible error $\epsilon$, scaling parameter $a>1$ and parameter $\rho \in (0,1)$, the problem (44) can be solved by using the existing unconstrained optimization methods, such as Newton, Interior point method, simplex method, etc. We use fminsearch function in MATLAB to estimate $\mathit{v}$ from (44). Then, we check whether the termination criteria are satisfied. If $\parallel h\left({\mathit{v}}_i\right)\parallel \le \delta$, the iteration terminates and ${\mathit{v}}_i$ is the near optimal solution of (37). Otherwise, go to the next step which is to control the convergence speed. If there exists $\rho \in (0,1)$ such that

$$\frac{\parallel h\left({\mathit{v}}_i\right)\parallel }{\parallel h\left({\mathit{v}}_{i-1}\right)\parallel }\ge \rho$$

(45)

we can increase the penalty factor s by setting $s=a·s$ ($a>1$). If not, then we leave the value of s unchanged. The penalty factor is used for multiplier update $\lambda =\lambda +s\parallel h\left({\mathit{v}}_i\right)\parallel$. The updated $\lambda$ and s are used to solve the unconstrained optimization problem (44) again.

Finally, the implementation process of the GLM algorithm is summarized in Algorithm 3.

Algorithm 3 The GLM algorithm.

1: Initialization:${\mathit{v}}_0={[{\widehat{\mathit{t}}}_{\mathrm{PLE}}^T,-1]}^T$, $\delta$, $\lambda$, s, $a>1$, $\rho \in (0,1)$. 2: for$i=0$; $i++$do 3:  Using $v_i$ as the initial point, find the minimum value of $\mathcal{L}(\mathit{v},\lambda ,s)$. 4:  If $\parallel h\left({\mathit{v}}_i\right)\parallel \le \delta$ holds, stop the iterations and get the final solution ${\mathit{v}}^{*}$, and go to step 7. Otherwise, go to step 5. 5:  If $\parallel h\left({\mathit{v}}_i\right)\parallel /\parallel h\left({\mathit{v}}_{i-1}\right)\parallel \ge \rho$, update the penalty coefficient $s=a·s$, and go to step 6. Otherwise, go to step 6 directly. 6:  Update $\lambda =\lambda -\rho \parallel h\left({\mathit{v}}_i\right)\parallel$ and go to step 3. 7:  Find the target position $\widehat{\mathit{t}}={\mathit{v}}^{*}(1:2)/{\mathit{v}}^{*}\left(3\right)$. 8: end for

4.4. Computational Complexity Analysis

In this subsection, we analyze the computational complexities of the proposed algorithms and compare with those of PLE and IRPLS. We only consider the asymptotic computational cost. For each algorithm, we separate the complexity calculation to small steps, e.g., multiplying two matrices of size $n\times m$ and $m\times p$ costs $\mathcal{O}\left(nmp\right)$. Let $L_1$, $L_2$ and $L_3$ denote the number of iterations with regarding to IRPLS, IGED and GLM, respectively. The iterative implementation of IRPLS, IGED and GLM, however, could be very sensitive to initialization. In this paper, we use the solution of PLE for initialization.

The computational costs are listed in

Table 1. Computational complexity.

Method	Operation	Cost
PLE	${\mathit{A}}^T\mathit{A}$, ${\mathit{A}}^T\mathit{b}$	$\mathcal{O}\left(4M\right)+\mathcal{O}\left(2M\right)$
	${\left({\mathit{A}}^T\mathit{A}\right)}^{-1}$	$\mathcal{O}\left(2^3\right)$
	${\left({\mathit{A}}^T\mathit{A}\right)}^{-1}{\mathit{A}}^T\mathit{b}$	$\mathcal{O}\left(4M\right)$
IRPLS	Compute ${\mathit{W}}^{\left(i\right)}\mathit{A}$ $L_1$ times	$\mathcal{O}\left(2L_{1}M\right)$
	Compute ${\mathit{W}}^{\left(i\right)}\mathit{b}$ $L_1$ times	$\mathcal{O}\left(L_{1}M\right)$
	Compute (29) $L_1$ times	$\mathcal{O}\left(L_1(10M+8)\right)$
IGED	Compute $\mathit{D}\mathit{K}$ $L_2$ times	$\mathcal{O}\left(3L_{2}M\right)$
	Perform GED $L_2$ times	$\mathcal{O}\left(3^3{L}_2\right)$
GLM	$L_3$ evaluations of BFGS	$\mathcal{O}\left({\sum }_{i=1}^{L_3}{N}_i{O}_i\left(M\right)\right)$

. Table 1 summarizes the computational complexities of the PLE and IRPLS, IGED and GLM algorithms. PLE requires the least amounts of computation. IRPLS and IGED exhibit the similar computational costs. For the GLM algorithm, we employ the Broyden–Fletcher–Goldfarb–Shanno (BFGS) [33] algorithm to solve the unconstrained problem presented in (44). BFGS is an iterative method for solving unconstrained nonlinear optimization problems by providing an approximation to the Hessian matrix. We assume that the number of iterations for BFGS is $N_i$ for the ith iteration of GLM. For each iteration, we assume the order of computational complexity of BFGS is $O_i\left(M\right)$. The GLM algorithm needs to estimate $\mathit{v}$ which is computationally demanding. The estimation procedure needs to repeat $L_3$ times since it is iterative.

5. Performance Bound

In this section, we present the Cramér-Rao Lower Bound (CRLB) for robust BOSL when the measurement noise is modeled as $S\alpha S$ distribution. We then discuss how the CRLB is calculated. The CRLB for any unbiased estimator of $\mathit{t}$ is derived by the inverse of the Fisher Information Matrix (FIM) [36], which is computed as

$$\mathrm{CRLB}\left(\mathit{t}\right)={\mathit{J}}^{-1}$$

(46)

where $\mathit{J}$ is the FIM given by the $2\times 2$ matrix

$$\mathit{J}=\mathbb{E}\left\{\left[\begin{array}{cc}{\left(\frac{\partial f\left(\tilde{\mathit{\varphi }}\right|\mathit{t})}{\partial t_x}\right)}^2& \frac{\partial f\left(\tilde{\mathit{\varphi }}\right|\mathit{t})}{\partial t_x}\frac{\partial f\left(\tilde{\mathit{\varphi }}\right|\mathit{t})}{\partial t_y}\\ \frac{\partial f\left(\tilde{\mathit{\varphi }}\right|\mathit{t})}{\partial t_x}\frac{\partial f\left(\tilde{\mathit{\varphi }}\right|\mathit{t})}{\partial t_y}& {\left(\frac{\partial f\left(\tilde{\mathit{\varphi }}\right|\mathit{t})}{\partial t_y}\right)}^2\end{array}\right]\right\}$$

(47)

where $f\left(\tilde{\mathit{\varphi }}\right|\mathit{t})$ denotes the Probability Density Function (PDF) of $\tilde{\mathit{\varphi }}$. Since $e_1,e_2,\dots ,e_M$ have the same standard dispersion, the PDF of ${\tilde{\varphi }}_m$ is

$$f_{\alpha }\left(e\right)=f_{\alpha }({\tilde{\varphi }}_m-{\varphi }_m)=\frac{1}{2\pi }{\int }_{-\infty }^{+\infty }exp\left(-i({\tilde{\varphi }}_m-{\varphi }_m)\tau -\gamma {\left|\tau \right|}^{\alpha }\right)d\tau$$

(48)

According to [37], the FIM $\mathit{J}$ for the parameters $t_x$ and $t_y$ is given by

$$\mathit{J}={\kappa }_{\alpha }\sum _{m=1}^M\left[\begin{array}{cc}{\left(\frac{\partial {\varphi }_m}{\partial t_x}\right)}^2& \frac{\partial {\varphi }_m}{\partial t_x}\frac{\partial {\varphi }_m}{\partial t_y}\\ \frac{\partial {\varphi }_m}{\partial t_x}\frac{\partial {\varphi }_m}{\partial t_y}& {\left(\frac{\partial {\varphi }_m}{\partial t_y}\right)}^2\end{array}\right]$$

(49)

where ${\kappa }_{\alpha }$ is the Fisher information for the location of $f_{\alpha }\left(e\right)$,

$${\kappa }_{\alpha }={\int }_{-\infty }^{+\infty }\frac{{\left(f_{\alpha }^{\prime }\left(e\right)\right)}^2}{f_{\alpha }\left(e\right)}de$$

(50)

and the expressions of $\frac{\partial {\varphi }_m}{\partial t_x}$ and $\frac{\partial {\varphi }_m}{\partial t_x}$ are given by

$$\frac{\partial {\varphi }_m}{\partial t_x}=-\frac{t_y-r_{y,m}}{\parallel \mathit{t}-{\mathit{r}}_m\parallel },\frac{\partial {\varphi }_m}{\partial t_y}=\frac{t_x-r_{x,m}}{\parallel \mathit{t}-{\mathit{r}}_m\parallel }$$

(51)

Note that the value of ${\kappa }_{\alpha }$ depends on the parameter $\alpha$. If $\alpha =1$, $f_{\alpha }\left(e\right)$ follows a Cauchy PDF and ${\kappa }_{\alpha }$ equals to $\frac{3}{5}$ for $\gamma =1$. In particular, when $\alpha =2$, $f_{\alpha }\left(e\right)$ is Gaussian with variance $2\gamma$ and ${\kappa }_{\alpha }=\gamma$. Therefore, the CRLB of $\mathit{t}$ in the case of $\alpha =2$ is consistent with the CRLB result given in [8].

However, when $1<\alpha <2$, $f_{\alpha }\left(e\right)$ has no closed-form expressions. Hence, computing $f_{\alpha }\left(e\right)$ involves numerically evaluating the integrals in (48). We use the MATLAB function STBLPDF to compute the pdf of the $S\alpha S$ distribution. The derivative of $f_{\alpha }\left(e\right)$ is also calculated by using numerical differentiation methods.

6. Simulations

In this section, we present three simulation examples to illustrate the performance of the proposed IGED and GLM algorithms and compare them with the PLE and IRPLS algorithm, as well as the CRLB. The localization performance is characterized with the bias norm and the root-mean-square-error (RMSE). The RMSE for source localization is defined as

$$\mathrm{RMSE}={\left(\frac{1}{N}\sum _{n=1}^N{\parallel {\widehat{\mathit{t}}}_n-\mathit{t}\parallel }^2\right)}^{1/2}$$

(52)

where N represents the total number of Monte Carlo runs, and ${\widehat{\mathit{t}}}_n$ denotes the source location estimate at the nth Monte Carlo run. The expression of the bias norm is given by

$$\mathrm{BIAS}=∥\frac{1}{N}\sum _{n=1}^N({\widehat{\mathit{t}}}_n-\mathit{t})∥$$

(53)

The RMSE given in (52) is bounded by the square root of the trace of the Cramér-Rao Lower Bound matrix

$$\mathrm{RMSE}\ge \sqrt{tr\left(\mathrm{CRLB}\right(\mathit{t}\left)\right)}$$

(54)

where $tr(·)$ denotes the trace of a matrix.

In the simulations, we examine the localization accuracy versus three kinds of parameters: (1) noise dispersion; (2) number of sensors; and (3) noise impulsiveness. All the sensors are uniformly placed in a $100\times 100$ ${\mathrm{m}}^2$ plane centered at $(50,50)$ m. The source is located at $(100,100)$ m. The IRPLS algorithm and the proposed IGED and GLM algorithms are iterative and they are all initialized by the PLE derived from (23) for ensuring a fair comparison. The termination parameters $ϵ$, $\epsilon$ and $\delta$ are set as $ϵ={10}^{-5}$, $\epsilon ={10}^{-10}$ and $\delta ={10}^{-5}$, respectively. For IRPLS, IGED and GLM, the maximum iteration is fixed at 200.

6.1. Various Levels of Noise Dispersion

In this example, we illustrate the RMSE and bias performance of PLE, IRPLS, IGED and GLM versus noise dispersion. For convenience, we use the generalized signal-to-noise ratio (GSNR) to characterize different noise levels. The GSNR is inversely proportional to the noise dispersion $\gamma$ and is calculated as $1/\gamma$ after normalized. The range of ${\gamma }^{1/\alpha }$ is from $2\pi /180$ radian to $6\pi /180$ radian. The characteristic exponent $\alpha$ is set to $1.5$ and the corresponding optimum value of p is $1.225$.

Figure 4 plots the RMSE and bias curves of the PLE, IRPLS, IGED and GLM algorithms together with the CRLB. For RMSE curves shown in Figure 4a, the green line with Asterisk is the RMSE value for the PLE algorithm, the blue line with Point is for IRPLS, the black line with Cross is for IGED, the red line with Square is for GLM and the blue dash line is the CRLB. The PLE method can not give accurate target location estimates in the presence of impulsive noise. The large bias norm formed by the outlier data in turn leads to a poor performance for the PLE. The least-norm estimators, including IRPLS, IGED and GLM, yield much better localization performance. However, the IRPLS method still has a large estimation bias because it cannot overcome the bias attributed to the correlation between the system matrix $\mathit{A}$ and the data vector $\mathit{b}$. On the other hand, the IGED and GLM algorithms utilizes the total Lp-norm optimization technique that can minimize the errors in $\mathit{A}$ and $\mathit{b}$ simultaneously. Therefore, the IGED and GLM algorithms are capable of outperforming the IRPLS in much lower RMSEs and biases.

Beyond the RMSE and bias performance, we need to explore the number of iterations and the computation time of IRPLS, IGED and GLM. The results are depicted in Figure 5. The algorithms run on a laptop with CPU i5-7200U @ 2.5 GHz and RAM 8 GB. The version of software is MATLAB 2017a. Collectively, the average number of iterations decreases when the GSNR increases. The GLM appears to have the least number of iterations. Its average number of iterations reduces from $1.05$ to $1.01$ when the GSNR ranges from $29.51$ to $153.3$. However, its computation time for single iteration is high (i.e., 0.1343 s), because the GLM algorithm has the steps of unconstrained optimization. On the contrary, the IGED algorithm has the least computing time since its computational complexity mainly lies in generalized eigenvalue decomposition. The IRPLS algorithm has high computing time due to the large number of iterations (over 30 iterations).

6.2. Different Number of Sensors

In this subsection, we compare the RMSE and bias performance of the PLE, IRPLS, IGED and GLM algorithms over the number of sensors M ranging from 10 to 30, five at a time, when ${\gamma }^{1/\alpha }$ is kept at $4\pi /180$ radian. The other parameters remain the same. The simulation results are shown in Figure 6 and Figure 7.

From Figure 6, we observe that the GLM algorithm has better RMSE and bias performance than that of PLE, IRPLS and IGED. When $M\ge 20$, the RMSE of GLM is close to the CRLB. The IGED algorithm slightly deviates from the CRLB, because the algorithm has not fully converged. Moreover, we also observe from Figure 7 that the average iterations is kept at $2.5$ and is almost the same as that of the GLM algorithm. Meanwhile, the bias of IGED is higher than that of GLM. Setting more iterations would help the IGED algorithm achieve better performance. As expected, the PLE exhibits unreliable results since it is not robust to the impulsive noise. The bias of IRPLE does not vanish as the number of sensors increases. This phenomenon is also validated in (33). In terms of computational complexity, the total computing time for GLM is about 0.3378 s when the number of sensors is fixed at 20. The IGED has the least amounts of computation and it only requires $0.01032$ s at $N=20$.

6.3. Various Values of Noise Impulsiveness

To further verify the effectiveness of the proposed algorithms, we examine the performance of the algorithms for different levels of noise impulsiveness. In the following figures drawing the simulation results, the value of the noise impulsiveness level $\alpha$ deviates from $1.9$ to $1.1$ and the corresponding optimum values of p are set as $1.546$, $1.430$, $1.348$, $1.282$, $1.225$, $1.174$, $1.127$, $1.083$, $1.041$ from (17). In this example, ${\gamma }^{1/\alpha }$ is fixed at $4\pi /180$ radian. The other parameters remain unchanged.

Figure 8 depicts the RMSE and bias performance with respect to various $\alpha$. As shown in Figure 8, the RMSE and the bias of PLE is much higher than that of IRPLS, IGED and GLM, which is caused by the impulsive noise. As $\alpha$ decreases, the influence of the noise impulsiveness becomes more significant. The level of noise impulsiveness also affects the bias and RMSE performance of IRPLS. The GLM algorithm has relatively small bias and its RMSE is very close to the CRLB when $\alpha \ge 1.4$. The IGED has the comparable RMSE performance with the GLM. However, the bias of IGED slightly deviates from that of GLM as $\alpha$ decreases. From Figure 9, we observe that the number of iterations of IRPLS, GLM and IGED decrease as $\alpha$ increases, i.e., the number of average iterations reduces from $38.8$ to $29.5$ for IRPLS, from $4.4$ to 1 for GLM and from $2.9$ to $2.2$ for IGED. The computation time also decreases when $\alpha$ increases. Again, the computation time of GLM is much higher than that of IGED.

6.4. Different Number of Iterations

In order to further verify the influence of the number of iterations on the proposed algorithm, we tested the performance of the algorithm for the increasing number of iterations. The simulation results are shown in Figure 10, where the number of iterations is set 50 for IRPLS and IGED, and 6 for GLM. The RMSE and bias are recorded every 10 iterations for IRPLS and IGED. The RMSE and bias performance of GLM is plotted for each iteration. In this example, ${\gamma }^{1/\alpha }$ is fixed at $4\pi /180$ radian. The number of sensors is kept at 20. The characteristic exponent $\alpha$ is set to $1.5$ and the corresponding optimum value of p is $1.225$.

Figure 10 describes the RMSE and bias performance of the IRPLS, IGED and GLM algorithms versus different iterations. As shown in Figure 10a, the RMSE result of IRPLS remains at $2.1$ m after 30 iterations. We observe that the bias of the IRPLS algorithm decreases as the number of iterations increasing, as shown in Figure 10b. However, the bias does not vanish. It still has $0.41$ m after 50 iterations. The RMSE value of IGED keeps at $1.8$ m after 30 iterations and the bias of IGED reduces to $0.12$ m using about 40 iterations. Differing from the IRPLS and IGED algorithms, the RMSE result of GLM reaches at $1.8$ m using only four iterations. After 5 iteratios, the bias of GLM reduces to 0.12 m. From Figure 10, we can observe that the RMSE and bias performance of IGED is almost the same as that of GLM if IGED runs a sufficient number of iterations.

6.5. Scalability Evaluation

We analyze the scalability of the IRPLS, IGED and GLM algorithms in this subsection for increasing number of sensors, from 20 to 100, and decreasing levels of noise impulsiveness, from $\alpha =1.1$ to $\alpha =1.9$. To examine the scalability of the iterative algorithms, the IRPLS, IGED and GLM algorithms share the same stopping condition. Let ${\widehat{\mathit{t}}}_i$ denote estimation result for the ith iteration. The criterion for stopping these iterative algorithms is given by $\parallel {\widehat{\mathit{t}}}_i-{\widehat{\mathit{t}}}_{i-1}\parallel <{10}^{-5}$.

Figure 11 compares the number of iterations, computing time for each iteration, RMSE and bias performance of the IRPLS, IGED and GLM algorithms versus different numbers of sensors. The other parameters remain the same as Section 6.4. It can be observed from Figure 11 that the number of iterations of these algorithms does not reduce as the number of sensors increases. Among them, the number of iterations for IRPLS and IGED maintains at 35, while that of the GLM algorithm keeps at 10. The computing time for each iteration, however, slightly increases as the number of sensors increasing. The RMSE and bias performance decreases as the number of sensors increases due to the fact that more normal sensors can be utilized. This result is also demonstrated in Figure 6. We also observe that the RMSE and bias curves of IGED and GLM are consistently lower than those of IRPLS, which shows the performance advantage of the proposed algorithms.

Figure 12 shows the number of iterations, computing time for each iteration, RMSE and bias performance of the IRPLS, IGED and GLM algorithms versus different levels of noise impulsiveness. The number of sensors is set at 20. Other parameters remain unchanged. From Figure 12a, we observe that the number of iterations for IRPLS scales from 88 to 11, 43 to 7 for IGED. However, the iteration value does not change for GLM and the number keeps at 7. That is because GLM performs optimization at each iteration to reduce the impact of impulsive noise. Moreover, the RMSE and bias performance of these algorithms decreases as increasing the value of $\alpha$. As expected, IGED and GLM exhibit almost the same performance. A small value of $\alpha$ results in a poorer IRPLE performance, as illustrated in Figure 12b. When $\alpha$ is greater than 1.6, the RMSE and bias curves of IRPLS reach those of IGED and GLM due to the fact that the level of noise impulsiveness has less impact on the bias performance.

7. Conclusions

In this paper, a total Lp-norm optimization method is presented to solve the BOSL problem and two algorithms, named IGED and GLM, are proposed to fulfill the optimization of the total least Lp-norm. By minimizing the errors in the system matrix and the data vector simultaneously for the least Lp-norm optimization, the proposed algorithms can overcome the bias arising from the correlation between the system matrix and the pseudolinear noise vector.

Simulation results show that the proposed IGED and GLM algorithms have much better RMSE and bias performance than the IRPLS algorithm. The number of iterations for GLM remains at a low level. With only a few iterations, the RMSEs of GLM can approach the CRLB. However, GLM expends more computation time than IGED under the same number of iterations, since GLM requires unconstrained optimization in each iteration.