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

This paper considers the problem of robust bearing-only source localization in impulsive noise with symmetric α-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.


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 α-Stable (Sα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 (α), skewness (β) and dispersion (γ). In particular, two special cases of α-Stable distribution can be obtained by letting impulsiveness parameter α take values of 1 and 2. One is Cauchy distribution (α = 1), and the other is Gaussian distribution (α = 2).
In the present study, we focus on the problem of robust BOSL with impulsive noise modeled as SαS distribution. In this situation, the methods derived using L2-norm opti-• 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α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.

Symmetric Alpha-Stable Distribution
The impulsive noise is more likely to exhibit outliers than normal noise. Studies [13] have shown that SαS distribution can model the impulsive noise better than Gaussian distribution due to the reason that the SαS densities have heavier tails than the Gaussian density. The characteristic function of SαS distribution is described as: where α (0 < α ≤ 2) denotes the characteristic exponent, γ indicates the dispersion parameter and δ stands for the location parameter. The SαS distribution is completely determined by these three parameters. The value of α indicates the heaviness of the tails of the density function. A small positive value of α implies high impulsiveness, while a value of α close to 2 shows a type of Gaussian-like shape. The dispersion γ performs like the variance and δ controls the location of the density function. A SαS distribution is called standard if δ = 0, γ = 1. Let x be a random variable that follows the standard SαS distribution. By taking the inverse Fourier transform of its characteristic function, the density function of x is of the form: It is known that a SαS distribution with characteristic exponent α only has finite moments for orders less than α, which are called the Fractional Lower Order Moments (FLOM). In particular, the FLOM of a standard SαS random variable x is given by where E{·} denotes expectation operator, is a constant that depends only on q and α and Γ(a) = +∞ 0 y (a−1) exp(−y)dy is the Gamma function. It is worth mentioning that the linear space of a SαS process is a Banach space for α ∈ [1, 2), and it is only a metric space for α ∈ (0, 1) [30]. Therefore, the tools of Hilbert space are not applicable when one solves a linear estimation problem with the SαS process. In the present study, we only focus on the case of α ∈ (1, 2).

Measurement Model
The problem of robust BOSL is depicted in Figure 1, where t = [t x , t y ] T denotes the unknown target location vector, r m = [r x,m , r y,m ] T represents the mth sensor location vector and θ m is the true bearing at sensor m. The nonlinear relationship between θ m , t and r m is given by [31] θ where tan −1 denotes the two-argument inverse tangent function and m = 1, 2, . . . , M. The bearing measurement taken at sensor m is given bỹ where n m follows independent zero mean SαS distribution with γ m representing the dispersion parameter and α accounting for the characteristic exponent. The dispersion γ m can vary with m. If the noise dispersion γ m is known a priori, then the noise term n m can be normalized. Let e m = γ −1/α m n m denote the normalized noise that has unit dispersion. The normalized measurement equation can be written as whereφ m = γ −1/α mθm is the normalized bearing measurement and φ m = γ −1/α m θ m is the normalized true bearing. Stacking (7) in the vector form yields

Minimum Dispersion Criterion
The main difficulty for parameter estimation with the α-stable process is that all non-Gaussian α-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αS random variable x as Thus, a suitable measure of the dispersion could be γ = x α α . Combining (3) and (9), one can easily obtain Therefore, we may estimate the source location by solving the following nonlinear Lp-norm optimization problem: where p ∈ (1, α). The definition of the Lp-norm of the vector ζ is In the present study, we do not consider the case of 0 < p ≤ 1, since the corresponding Lp-norm cost function is not differentiable. For non-normalized bearing measurements, the Lp-norm objective function in (11) becomes 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 iŝ where ∇J(t) denotes the Jacobian matrix, If e m > 0, sign(e m ) = 1 and −1 otherwise. The statistical properties of the nonlinear Lp-norm minimizer have been studied in the literature [24]. The theoretical covariance B is related to the value of p, and it is given by [24]: where ∇φ denotes the Jacobian matrix of φ with respect to t. By minimizing the scalar term of B, we obtain the optimal choice of p where p o denotes the optimal value of p.

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 where u m is the true angle vector between r m and t,ũ m is the measured angle vector starting from r m and produces the noisy bearingθ m with respect to the horizontal direction and ε m is the error vector. Let µ m = [cosθ m , sinθ m ] T and ν m = [sinθ m , − cosθ m ] T denote two orthogonal unit trigonometric vectors. Then,ũ m and ε m are given bỹ where · denotes Euclidean norm. Note from (19) thatũ T m ν m = 0. Substituting (19) and (20) into (18) and multiplying (18) where ξ m = u m sin n m is a nonlinear transformed measurement error. Collecting the pseudolinear equation errors as a vector ξ = [ξ 1 , ξ 2 , . . . , ξ K ] T , we obtain where A = [ν 1 , ν 2 , . . . , ν M ] T , b = [ν T 1 r 1 , . . . , ν T M r M ] T are the measurement matrix and vector, respectively.
The PLE requires that t be estimated by minimizing At − b 2 2 with respect to t in an L2-norm optimization sense, and the solution iŝ 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.

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: where G = diag(g 1 , g 2 , . . . , The above interpretation suggests an iterative reweighted pseudolinear least-squares with the weight matrix for the i-th iteration given by 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 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, where κ < 1 is the homotopy parameter. The source location can be estimated by performing least-squares estimation: The IRPLS algorithm is stopped when Finally, the whole process of the IRPLS algorithm is summarized in Algorithm 1.

4:
Determine the value of p for each iteration using (28).

5:
Perform weighted least squares estimation using (29). 6: if t (i−1) −t (i) ≤ , 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 b is subjected to errors. This is not the case, since the system matrix A is corrupted with measurement noises as well. The correlation between A and b causes the IRPLS estimator to be inconsistent.
In this subsection, we analyze such bias attributed to the correlation between A and b.
Lett * denote the final solution of the IRPLS algorithm. This solution satisfieŝ andθ * m denotes the mth bearing obtained from (5) by usingt * .
According to (22) and (30), we can get The analytical bias of IRPLS is given by [14] E{∆t} where the weighting matrix W (t) is computed from (31) using true source location parameter t.

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 b. But in fact, the system matrix A also has a measurement residual. The IRPLS algorithm will inevitably cause a large bias, because the system matrix A and the system vector b are statistically correlated.

Method Description
In order to improve the performance of the IRPLS algorithm and reduce the bias caused by the correlation between A and b, we develop the TLPO method in this subsection. Unlike the IRPLS method, the TLPO method uses the correction matrix ∆A and the correction vector ∆b to compensate the system matrix A and system vector b, respectively. The normalized equation can be written as where ∆A is the perturbation matrix of A and ∆b is the perturbation vector of b. Let and v = [t, −1] T be the augmented matrix and vector, and Equation (34) can be rewritten as where Ω = G[∆A, ∆b] is the error matrix. In order to reduce the bias of IRPLS, we use the TLPO approach to minimize the perturbation matrix ∆A and perturbation vector ∆b simultaneously, and the TLPO problem for BOTL can be formulated as 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 t from the minimization of Ω p subject to (K + Ω)v = 0 is equivalent to Proof. From the constraint (K + Ω)v = 0, we have Kv p = Ωv p . Let ω T m denote the mth row of Ω. Ωv p p becomes for p and q satisfy the equation 1 p + 1 q = 1 and 1 < p < α. To derive (38), we used the properties of hölder's inequality [34]. If the Lq-norm of v satisfies v q = 1, then the inequality (38) becomes Ωv p ≤ Ω p . Therefore, we can conclude that the minimization of Ω p subject to (K + Ω)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.

The IGED Algorithm
The optimization problem in Proposition 1 is equivalent to 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 3 |) and λ is the Lagrangian multiplier. Taking the partial derivative of L(v, λ) with respect to v and setting it to zero yields After multiplying both sides of (41) with v T , it can be found that only λ = v T K T DKv is the cost to be minimized. So, the optimal solution v * is the eigenvector of the smallest generalized eigenvalue of (K T DK, C). Figure 2 gives a flowchart to depict the whole process of the IGED algorithm. The algorithm starts from the initialization procedure, where v 0 is set as [t T PLE , −1] T andt PLE is obtained from (23). Given p and q, matrices C and D can be computed, respectively. After performing generalized eigenvalue decomposition from the pair (K T DK, C), the solution of v in (41) is the generalized eigenvector of (K T DK, C) that gives the minimum generalized eigenvalue. Since IGED is an iterative algorithm, it requires the convergence check. Let λ i and λ i+1 be the minimum generalized eigenvalues at the ith and i + 1th iterations, respectively. If there exists ε such that then we would get the final solution v * and the estimate for source location iŝ In summary, the IGED Algorithm 2 can be implemented by the following procedures. Compute C and D.

4:
Perform generalized eigenvalue decomposition of the pair (K T DK, C).

5:
Select the smallest generalized eigenvalue and its corresponding generalized eigenvector. 6: Set λ i and λ i+1 as the minimum generalized eigenvalues for iterations i and i + 1. Get the source location estimatet using (43). 8: end for

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: where g(v) = Kv p and h(v) = 1 − v q . Compared with the Lagrangian function, the GLM cost function has an additional term s/2 h(v) 2 2 . This item is a punishment for violation of constraint h(v) = 0. The punishment parameter s determines the degree of punishment, and it is generally sufficiently large. λ is a positive scalar, and the term λh(v) is to ensure that the optimal solution is a strict local minimum point of 2 2 under the condition of obeying the constraint h(v) = 0. Figure 3 shows the flow diagram of the GLM algorithm. For the initialization step, we set v 0 = [t T PLE , −1] T . Given initial multiplier λ, penalty factor s, admissible error ε, scaling parameter a > 1 and parameter ρ ∈ (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 v from (44). Then, we check whether the termination criteria are satisfied. If h(v i ) ≤ δ, the iteration terminates and 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 ρ ∈ (0, 1) such that 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 λ = λ + s h(v i  Finally, the implementation process of the GLM algorithm is summarized in Algorithm 3.

Algorithm 3
The GLM algorithm.
Using v i as the initial point, find the minimum value of L(v, λ, s).

4:
If h(v i ) ≤ δ holds, stop the iterations and get the final solution v * , and go to step 7. Otherwise, go to step 5.

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 × m and m × p costs O(nmp). 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. 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 (M). The GLM algorithm needs to estimate v which is computationally demanding. The estimation procedure needs to repeat L 3 times since it is iterative.

Method
Operation Cost

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αS distribution. We then discuss how the CRLB is calculated. The CRLB for any unbiased estimator of t is derived by the inverse of the Fisher Information Matrix (FIM) [36], which is computed as where J is the FIM given by the 2 × 2 matrix where f (φ|t) denotes the Probability Density Function (PDF) ofφ. Since e 1 , e 2 , . . . , e M have the same standard dispersion, the PDF ofφ m is According to [37], the FIM J for the parameters t x and t y is given by where κ α is the Fisher information for the location of f α (e), Note that the value of κ α depends on the parameter α. If α = 1, f α (e) follows a Cauchy PDF and κ α equals to 3 5 for γ = 1. In particular, when α = 2, f α (e) is Gaussian with variance 2γ and κ α = γ. Therefore, the CRLB of t in the case of α = 2 is consistent with the CRLB result given in [8].
However, when 1 < α < 2, f α (e) has no closed-form expressions. Hence, computing f α (e) involves numerically evaluating the integrals in (48). We use the MATLAB function STBLPDF to compute the pdf of the SαS distribution. The derivative of f α (e) is also calculated by using numerical differentiation methods.

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 where N represents the total number of Monte Carlo runs, andt n denotes the source location estimate at the nth Monte Carlo run. The expression of the bias norm is given by The RMSE given in (52) is bounded by the square root of the trace of the Cramér-Rao Lower Bound matrix RMSE ≥ tr(CRLB(t)) (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 × 100 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 , ε and δ are set as = 10 −5 , ε = 10 −10 and δ = 10 −5 , respectively. For IRPLS, IGED and GLM, the maximum iteration is fixed at 200.

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 γ and is calculated as 1/γ after normalized. The range of γ 1/α is from 2π/180 radian to 6π/180 radian. The characteristic exponent α 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 A and the data vector b. On the other hand, the IGED and GLM algorithms utilizes the total Lp-norm optimization technique that can minimize the errors in A and 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).

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 γ 1/α is kept at 4π/180 radian. The other parameters remain the same. The simulation results are shown in Figures 6 and 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 ≥ 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.

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 α deviates from  (17). In this example, γ 1/α is fixed at 4π/180 radian. The other parameters remain unchanged. Figure 8 depicts the RMSE and bias performance with respect to various α. 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 α 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 α ≥ 1.4. The IGED has the comparable RMSE performance with the GLM. However, the bias of IGED slightly deviates from that of GLM as α decreases. From Figure 9, we observe that the number of iterations of IRPLS, GLM and IGED decrease as α 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 α increases. Again, the computation time of GLM is much higher than that of IGED.

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, γ 1/α is fixed at 4π/180 radian. The number of sensors is kept at 20. The characteristic exponent α 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.

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 α = 1.1 to α = 1.9. To examine the scalability of the iterative algorithms, the IRPLS, IGED and GLM algorithms share the same stopping condition. Lett i denote estimation result for the ith iteration. The criterion for stopping these iterative algorithms is given by t i −t i−1 < 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 α. As expected, IGED and GLM exhibit almost the same performance. A small value of α results in a poorer IRPLE performance, as illustrated in Figure 12b. When α 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.

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.