Efficient Bi-Iterative Method for Source Position and Propagation Speed Estimation Using TDOA Measurements

Abstract

This paper develops an efficient bi-iterative source location and propagation speed estimation method utilizing time difference of arrival (TDOA) measurements. The source location and propagation speed estimation is a nonlinear problem due to the nonlinearity in the TDOA measurement equations. The newly developed bi-iterative method computes the source location and propagation speed alternately. The asymptotic convergence of the new bi-iterative method is theoretically analyzed. First-order perturbation analysis is applied to the newly developed solution to derive its bias and variance. The first-order analytical results show that the proposed method provides approximately unbiased source position and propagation speed estimates for low noise levels and the accuracy of these estimates approaches the Cramer-Rao lower bound (CRLB). The extension of the new bi-iterative method to the more general situation where the sensor locations are subject to random errors is also presented. Simulation studies are given to show the good performance of the proposed method.

1. Introduction

Passive source localization using time difference of arrival (TDOA) measurements has received considerable attention and has been widely applied in target tracking [1, 2], navigation [3], sensor networks [4, 5], and wireless communications [6, 7]. In the past decades, a number of efficient algorithms such as those in [5–9] were presented for TDOA-based source localization. Nevertheless, all these works assume that the signal propagation speed is known a priori so that the obtained TDOA measurements can be converted into range differences for source positioning. In practical applications such as seismic exploration [10], tangible interface for human-computer interaction [11], and underwater acoustic [12], the propagation speed is unknown and depends strongly on the propagation medium. In this case, the unknown propagation speed needs to be estimated jointly with the source location. For this problem, Mahajan and Walworth [13] proposed an unconstrained least-squares (LS) method, in which the nonlinear TDOA measurement equations are converted into pseudolinear ones by introducing two auxiliary variables. The source location and two auxiliary variables are then jointly estimated in LS sense. Reed et al. [14] selected the LS solution as the starting value and developed a four-step method which alternately estimates the source location and propagation speed. Recently, Zheng et al. [15] proposed a three-stage approach to simultaneously compute the source location and propagation speed. Very interestingly, the accuracy of the source location and propagation speed estimates approximates the CRLB for sufficiently small noise conditions. A disadvantage of this three-stage method is that it produces two results, where only one is the true solution. Furthermore, it may generate complex solutions when finding the square root in its last stage. Annibale and Rabenstein [16] investigated the influence of a wrongly presumed propagation speed because of temperature variations on the positioning accuracy of two well-known location methods. Oyzerman and Amar [17] extended the spherical intersection (SX) technique in [5] to the joint source position and propagation speed estimation. Annibale and Rabenstein [18] derived two closed-form methods for estimation of the propagation speed using time of arrival (TOA) and TDOA measurements.

However, these works [13–18] mentioned above are based on the assumption that the sensor locations are accurately known. This assumption may be sometimes invalid in modern practical scenarios; for example, in wireless sensor networks (WSNs), the sensors are frequently localized level by level [19]. Due to this level-by-level procedure, the sensor locations generally include random errors in which case the sensor location errors need to be taken into consideration [20]. In this paper, we propose an efficient bi-iterative source position and propagation speed estimation method using TDOA measurements and compare its performance with the LS solution [13], four-step method [14], three-step weighted least-squares (WLS) method [15], and the CRLB. Our contributions include the following.

• (1)

  We develop an efficient bi-iterative method for the source position and propagation speed estimation.

• (2)

  The asymptotic convergence of the new bi-iterative estimation method is analyzed.

• (3)

  The proposed solution is shown to approach the CRLB for low Gaussian TDOA measurement noise levels.

• (4)

  We study the CRLB of the source position and propagation speed with sensor location errors and extend the bi-iterative method to the more general case where the sensor locations are subject to random errors.

The remainder of the paper is structured as follows. Section 2 describes the problem of computing source position and propagation speed under a Gaussian TDOA noise model. Development of the new bi-iterative estimation method is presented in Section 3. Besides, the asymptotic convergence and approximate efficiency of the bi-iterative method are analyzed under low noise level conditions. Section 4 extends the bi-iterative estimation method to the general situation where the sensor locations are subject to random errors. Simulation studies are given in Section 5 to validate the theoretical analysis and to test our bi-iterative estimation method's performance. Section 6 concludes the paper.

Throughout the paper, we utilize boldface lowercase letters to stand for vectors and boldface capitals to represent matrices. ${\mathbf{A}}_{ij}$ and ${\mathbf{A}}_k$ denote the ith to jth columns and the kth column of matrix A, respectively. $\mathbf{A}⪰\mathbf{B}$ means that $\mathbf{A}-\mathbf{B}$ is positive semidefinite. ${\mathbf{0}}_{k\times l}$ stands for the $k\times l$ zero matrix and ${\mathbf{1}}_k$ denotes the $k\times \mathrm{1}$ vector of ones. Moreover, $\mathrm{diag}\left(*\right)$ denotes a diagonal matrix, ${‖*‖}_{\mathrm{2}}$ represents Euclidean norm, and $E\left[*\right]$ stands for the expectation operator.

2. Problem Formulation

Consider a two-dimensional plane where a sensor array is composed of M sensors at exactly known locations ${\mathbf{s}}_j={\left[x_j,y_j\right]}^T$ , $j=\mathrm{1,2},\dots ,M$ . It is straightforward to generalize the methodology to three-dimensional case. We postulate that all these sensors are not lying on a straight line, which guarantees that the matrix A in the proposed method has full column rank. Let $\mathbf{u}={\left[x,y\right]}^T$ be the position of the source which emits signal in free-field conditions to be determined.

Denote by $r_j$ the true distance from the source to the jth sensor:

$$\begin{array}{cc}& r_j={‖\mathbf{u}-{\mathbf{s}}_j‖}_{\mathrm{2}},\hspace{1em}j=\mathrm{1,2},\dots ,M.\end{array}$$

(1)

Without losing generality, we set the first sensor to be the reference sensor. The TDOA measurement obtained from sensor pair ${\mathbf{s}}_i$ and ${\mathbf{s}}_{\mathrm{1}}$ would be [20], after taking into account the measurement noise,

$$\begin{array}{cc}& {\tilde{t}}_{i\mathrm{1}}=t_{i\mathrm{1}}+n_{i\mathrm{1}}=\frac{r_{i\mathrm{1}}}{c}+n_{i\mathrm{1}},\hspace{1em}i=\mathrm{2,3},\dots ,M,\end{array}$$

(2)

where $t_{i\mathrm{1}}=r_{i\mathrm{1}}/c$ , $r_{i\mathrm{1}}=r_i-r_{\mathrm{1}}$ , and c is the unknown propagation speed. $n_{i\mathrm{1}}$ is the additive measurement noise. For simplicity, we presume that all the TDOA measurements used for source localization come from line-of-sight propagation.

Equation (2) can be arranged into matrix form

$$\begin{array}{cc}& \tilde{\mathbf{t}}={\left[{\tilde{t}}_{\mathrm{21}},{\tilde{t}}_{\mathrm{31}},\dots ,{\tilde{t}}_{M\mathrm{1}}\right]}^T=\mathbf{t}+\mathbf{n},\end{array}$$

(3)

where $\mathbf{t}={\left[t_{\mathrm{21}},t_{\mathrm{31}},\dots ,t_{M\mathrm{1}}\right]}^T$ and $\mathbf{n}={\left[n_{\mathrm{21}},n_{\mathrm{31}},\dots ,n_{M\mathrm{1}}\right]}^T$ . n is assumed to be Gaussian distributed with zero-mean and known covariance matrix ${\mathbf{Q}}_t$ . This assumption has been frequently adopted in literatures such as [14, 15, 18] for source location and propagation speed estimation.

Notice that TDOA measurement equations (2) are not linearly related to the source position u and propagation speed c. The problem of interest is to estimate the source position u and propagation speed c as accurately as possible via exploiting measurement equations (3).

3. Bi-Iterative Estimation Method

3.1. Bi-Iterative Source Position and Propagation Speed Estimation Method

Achieving high source position and propagation speed estimation accuracy is far from being straightforward due to the nonlinearity implied in the TDOA localization equations. According to (2), the TDOA measurement equation can be rearranged as $c{\tilde{t}}_{i\mathrm{1}}+r_{\mathrm{1}}=r_i+c{n}_{i\mathrm{1}}$ . Squaring both sides of $c{\tilde{t}}_{i\mathrm{1}}+r_{\mathrm{1}}=r_i+c{n}_{i\mathrm{1}}$ and considering (1), we obtain

$$\begin{array}{l}{\left({\mathbf{s}}_i-{\mathbf{s}}_{\mathrm{1}}\right)}^T\left(\mathbf{u}-{\mathbf{s}}_{\mathrm{1}}\right)+\mathrm{0.5}{\tilde{t}}_{i\mathrm{1}}^{\mathrm{2}}{c}^{\mathrm{2}}+{\tilde{t}}_{i\mathrm{1}}c{r}_{\mathrm{1}}-\mathrm{0.5}{‖{\mathbf{s}}_i-{\mathbf{s}}_{\mathrm{1}}‖}_{\mathrm{2}}^{\mathrm{2}}\\ \hspace{1em}=c{r}_i{n}_{i\mathrm{1}}+c^{\mathrm{2}}{n}_{i\mathrm{1}}^{\mathrm{2}}.\end{array}$$

(4)

By stacking (4) for $i=\mathrm{2,3},\dots ,M$ together, we have

$$\begin{array}{cc}& \mathbf{G}\left(\mathbf{u}-{\mathbf{s}}_{\mathrm{1}}\right)+\mathrm{0.5}\tilde{\mathbf{t}}\odot \tilde{\mathbf{t}}{c}^{\mathrm{2}}+\tilde{\mathbf{t}}c{r}_{\mathrm{1}}-\mathbf{h}=\mathbf{\eta },\end{array}$$

(5)

where $\mathbf{G}={\left[{\mathbf{s}}_{\mathrm{2}}-{\mathbf{s}}_{\mathrm{1}},{\mathbf{s}}_{\mathrm{3}}-{\mathbf{s}}_{\mathrm{1}},\dots ,{\mathbf{s}}_M-{\mathbf{s}}_{\mathrm{1}}\right]}^T$ , $\mathbf{h}=\mathrm{0.5}{\left[{‖{\mathbf{s}}_{\mathrm{2}}-{\mathbf{s}}_{\mathrm{1}}‖}_{\mathrm{2}}^{\mathrm{2}},{‖{\mathbf{s}}_{\mathrm{3}}-{\mathbf{s}}_{\mathrm{1}}‖}_{\mathrm{2}}^{\mathrm{2}},\dots ,{‖{\mathbf{s}}_M-{\mathbf{s}}_{\mathrm{1}}‖}_{\mathrm{2}}^{\mathrm{2}}\right]}^T$ , and the symbol $\odot$ represents the Schur product operator (element-by-element multiplication). On the right-hand side of (5), $\mathbf{\eta }=\mathbf{B}\mathbf{n}+c^{\mathrm{2}}\mathbf{n}\odot \mathbf{n}$ is the equation error vector and B is defined as $\mathbf{B}=\mathrm{diag}\left(c\left[r_{\mathrm{2}},r_{\mathrm{3}},\dots ,r_M\right]\right)$ .

The source position u and propagation speed c can be achieved by minimizing the following WLS problem:

$$\begin{array}{cc}& J\left(\mathbf{u},c\right)={\mathbf{\eta }}^T\mathbf{W}\mathbf{\eta }={\left(\mathbf{A}\mathbf{\phi }-\mathbf{h}\right)}^T\mathbf{W}\left(\mathbf{A}\mathbf{\phi }-\mathbf{h}\right),\end{array}$$

(6)

where W is the weighting matrix chosen as [21] $\mathbf{W}=E{\left[\mathbf{\eta }{\mathbf{\eta }}^T\right]}^{-\mathrm{1}}$ , $\mathbf{A}=\left[\mathbf{G},\mathrm{0.5}\tilde{\mathbf{t}}\odot \tilde{\mathbf{t}},\tilde{\mathbf{t}}\right]$ , and $\mathbf{\phi }={\left[{\left(\mathbf{u}-{\mathbf{s}}_{\mathrm{1}}\right)}^T,c^{\mathrm{2}},c{r}_{\mathrm{1}}\right]}^T$ . When the TDOA measurement noises $n_{i\mathrm{1}}$ , $i=\mathrm{2,3},\dots ,M$ , are relatively small, the second-order noise terms can be ignored and the weighting matrix can be approximated by $\mathbf{W}\approx {\left(\mathbf{B}{\mathbf{Q}}_t{\mathbf{B}}^T\right)}^{-\mathrm{1}}$ . It should be noted that the source position u and propagation speed c are coupled together by $c{r}_{\mathrm{1}}$ . In order to eliminate this coupling effect, we herein follow the bi-iterative scheme utilized in [22] to alternately estimate the source position and propagation speed.

Firstly, let us consider the situation where the source position u is fixed. From (5) and (6), we have

$$\begin{array}{cc}& \underset{c}{\min }\hspace{0.17em}J\left(\mathbf{u},c\right)=\underset{c}{\min }{\left({\mathbf{a}}_{\mathrm{2}}{c}^{\mathrm{2}}+{\mathbf{a}}_{\mathrm{1}}c+{\mathbf{a}}_{\mathrm{0}}\right)}^T\mathbf{W}\left({\mathbf{a}}_{\mathrm{2}}{c}^{\mathrm{2}}+{\mathbf{a}}_{\mathrm{1}}c+{\mathbf{a}}_{\mathrm{0}}\right),\end{array}$$

(7)

where ${\mathbf{a}}_{\mathrm{0}}=\mathbf{G}\left(\mathbf{u}-{\mathbf{s}}_{\mathrm{1}}\right)-\mathbf{h}$ , ${\mathbf{a}}_{\mathrm{1}}=r_{\mathrm{1}}\tilde{\mathbf{t}}$ , and ${\mathbf{a}}_{\mathrm{2}}=\mathrm{0.5}\tilde{\mathbf{t}}\odot \tilde{\mathbf{t}}$ . Thus, the cost function $J\left(\mathbf{u},c\right)$ is a quartic function with respect to the propagation speed c. By taking the derivative of $J\left(\mathbf{u},c\right)$ with respect to c and equating the result to zero, we obtain

$$\begin{array}{cc}& b_3{c}^3+b_2{c}^2+b_{1}c+b_0=\mathrm{0},\end{array}$$

(8)

where the coefficients are $b_0={\mathbf{a}}_{\mathrm{1}}^T\mathbf{W}{\mathbf{a}}_{\mathrm{0}}$ , $b_{\mathrm{1}}=\mathrm{2}{\mathbf{a}}_{\mathrm{2}}^T\mathbf{W}{\mathbf{a}}_{\mathrm{0}}+{\mathbf{a}}_{\mathrm{1}}^T\mathbf{W}{\mathbf{a}}_{\mathrm{1}}$ , $b_2=\mathrm{3}{\mathbf{a}}_{\mathrm{2}}^T\mathbf{W}{\mathbf{a}}_{\mathrm{1}}$ , and $b_{\mathrm{3}}=\mathrm{2}{\mathbf{a}}_{\mathrm{2}}^T\mathbf{W}{\mathbf{a}}_{\mathrm{2}}$ . The solution of the cubic polynomial equation in (8) can be obtained efficiently utilizing Cardano's method [23].

Secondly, consider another situation where the propagation speed c is fixed; we arrive at

$$\begin{array}{cc}& \underset{\mathbf{u}}{\min }\hspace{0.17em}J\left(\mathbf{u},c\right)=\underset{\mathbf{u}}{\min }{\left(\mathbf{H}\mathbf{y}-{\mathbf{h}}_u\right)}^T\mathbf{W}\left(\mathbf{H}\mathbf{y}-{\mathbf{h}}_u\right),\end{array}$$

(9)

where $\mathbf{H}=\left[\mathbf{G},c\tilde{\mathbf{t}}\right]$ , $\mathbf{y}={\left[{\left(\mathbf{u}-{\mathbf{s}}_{\mathrm{1}}\right)}^T,r_1\right]}^T$ , and ${\mathbf{h}}_u=\mathbf{h}-\mathrm{0.5}{c}^{\mathrm{2}}\tilde{\mathbf{t}}\odot \tilde{\mathbf{t}}$ . It is noticed that (9) is nonlinear and nonconvex with respect to u since $\left(\mathbf{u}-{\mathbf{s}}_{\mathrm{1}}\right)$ is related to $r_{\mathrm{1}}$ through (1). Here we apply the semidefinite relaxation (SDR) technique [24] to convert (9) to an approximate but convex problem. The minimization problem can be reexpressed as

$$\begin{array}{l}\underset{\mathbf{y}}{\min }\phantom{\text{iiiiiiiiiii}}{\left(\mathbf{H}\mathbf{y}-{\mathbf{h}}_u\right)}^T\mathbf{W}\left(\mathbf{H}\mathbf{y}-{\mathbf{h}}_u\right)\\ \text{subject}\mathrm{\hspace{0.17em}\hspace{0.17em}}\text{to}\hspace{1em}\mathbf{y}\left(\mathrm{3}\right)={‖\mathbf{y}\left(\mathrm{1}:\mathrm{2}\right)‖}_{\mathrm{2}},\end{array}$$

(10)

where $\mathbf{y}\left(\mathrm{1}:\mathrm{2}\right)$ and $\mathbf{y}\left(\mathrm{3}\right)$ denote the first to second elements and the third element of vector y, respectively. Expressing $\mathbf{H}\mathbf{y}-{\mathbf{h}}_u$ as $\mathbf{H}\mathbf{y}-{\mathbf{h}}_u=\left[\mathbf{H},-{\mathbf{h}}_u\right]{\left[{\mathbf{y}}^T,\mathrm{1}\right]}^T$ and writing the cost function in (10) into

$$\begin{array}{l}{\left(\mathbf{H}\mathbf{y}-{\mathbf{h}}_u\right)}^T\mathbf{W}\left(\mathbf{H}\mathbf{y}-{\mathbf{h}}_u\right)\\ \hspace{1em}=\left[{\mathbf{y}}^T,\mathrm{1}\right]\mathbf{F}\left[\begin{array}{l}\mathbf{y}\\ \mathrm{1}\end{array}\right]=\text{trace}\left\{\left[\begin{array}{ll}\mathbf{Y}& \mathbf{y}\\ {\mathbf{y}}^T& \mathrm{1}\end{array}\right]\mathbf{F}\right\},\end{array}$$

(11)

where $\text{trace}\left(*\right)$ denotes the trace operator, $\mathbf{F}=\left[\begin{array}{cc}{\mathbf{H}}^T\mathbf{W}\mathbf{H}& -{\mathbf{H}}^T\mathbf{W}{\mathbf{h}}_u\\ -{\mathbf{h}}_u^T\mathbf{W}\mathbf{H}& {\mathbf{h}}_u^T\mathbf{W}{\mathbf{h}}_u\end{array}\right]$ , and $\mathbf{Y}=\mathbf{y}{\mathbf{y}}^T$ , accordingly, the constraint in (10) can be written into $\mathbf{Y}\left(\mathrm{3,3}\right)=\text{trace}\left\{\mathbf{Y}\left(\mathrm{1}:\mathrm{2,1}:\mathrm{2}\right)\right\}$ , where $\mathbf{Y}\left(\mathrm{1}:\mathrm{2,1}:\mathrm{2}\right)$ represents the upper-left two-by-two submatrix of Y and $\mathbf{Y}\left(\mathrm{3,3}\right)$ is the $\left(\mathrm{3,3}\right)$ th entry of Y. Hence, the optimization problem (10) can be equivalently written in the following form:

$$\begin{array}{l}\underset{\mathbf{Y},\mathbf{y}}{\min }\phantom{\text{iiiiiiiiiiiii}}\text{trace}\left\{\left[\begin{array}{ll}\mathbf{Y}& \mathbf{y}\\ {\mathbf{y}}^T& \mathrm{1}\end{array}\right]\mathbf{F}\right\}\\ \text{subject}\mathrm{\hspace{0.17em}\hspace{0.17em}}\text{to}\hspace{1em}\mathbf{Y}\left(\mathrm{3,3}\right)=\text{trace}\left\{\mathbf{Y}\left(\mathrm{1}:\mathrm{2,1}:\mathrm{2}\right)\right\}\\ \phantom{\text{iiiiiiiiiiiiiiiiiii}}\mathbf{Y}=\mathbf{y}{\mathbf{y}}^T.\end{array}$$

(12)

The minimization problem (12) is just as difficult to solve as problem (10) because of the nonconvexity in the second equality constraint. Thus, replacing the equality constraint $\mathbf{Y}=\mathbf{y}{\mathbf{y}}^T$ by an inequality $\mathbf{Y}⪰\mathbf{y}{\mathbf{y}}^T$ to meet the convex specification and establishing the equivalence between $\mathbf{Y}⪰\mathbf{y}{\mathbf{y}}^T$ and $\left[\begin{array}{cc}\mathbf{Y}& \mathbf{y}\\ {\mathbf{y}}^T& \mathrm{1}\end{array}\right]⪰{\mathbf{0}}_{\mathrm{4}\times \mathrm{4}}$ utilizing the Schur complement [24], we have the following SDP:

$$\begin{array}{l}\underset{\mathbf{Y},\mathbf{y}}{\min }\phantom{\text{iiiiiiiiiiii}}\text{trace}\left\{\left[\begin{array}{ll}\mathbf{Y}& \mathbf{y}\\ {\mathbf{y}}^T& \mathrm{1}\end{array}\right]\mathbf{F}\right\}\\ \text{subject}\mathrm{\hspace{0.17em}\hspace{0.17em}}\text{to}\hspace{1em}\mathbf{Y}\left(\mathrm{3,3}\right)=\text{trace}\left\{\mathbf{Y}\left(\mathrm{1}:\mathrm{2,1}:\mathrm{2}\right)\right\}\\ \phantom{\text{iiiiiiiiiiiiiiiiiii}}\left[\begin{array}{ll}\mathbf{Y}& \mathbf{y}\\ {\mathbf{y}}^T& \mathrm{1}\end{array}\right]⪰{\mathbf{0}}_{\mathrm{4}\times \mathrm{4}}.\end{array}$$

(13)

That is to say, the nonconvex optimization problem in (9) is transformed into a convex one. The optimal value of (13) is less than or equal to the optimal value of (12), since we minimize the same cost function over a larger set. The global optimum can be obtained by applying existing numerical methods, for instance, the interior-point methods.

Based upon the above analysis, the procedures of the proposed bi-iterative method for source position and propagation speed estimation are given as follows.

• (1)

  Obtain the initial source position and propagation speed estimates ${\mathbf{u}}_{\mathrm{0}}$ and $c_0$ using the solution in [13].

•  

  For $k=\mathrm{1,2},\dots$ , repeat the following two steps until ${‖{\mathbf{u}}_{k+\mathrm{1}}-{\mathbf{u}}_k‖}_{\mathrm{2}}<\delta$ ( $\delta =\mathrm{1}{\mathrm{0}}^{-\mathrm{3}}$ ).

• (2)

  Put $c_{k-1}$ into (9) and obtain the source position estimate ${\mathbf{u}}_k$ by solving the convex optimization problem (13) utilizing the convex optimization toolbox CVX [24].

• (3)

  Substituting ${\mathbf{u}}_k$ into (8) produces a cubic equation in c, and the roots of this cubic equation can be obtained efficiently utilizing Cardano's method. We select the one that gives the smallest $J\left({\mathbf{u}}_k,c\right)$ in (7).

3.2. Asymptotic Convergence

Here, we show that the bi-iterative estimation method has asymptotic convergence. It is apparent that the cost function $J\left(\mathbf{u},c\right)$ given in (6) is continuous because it is differentiable. From the definition of $J\left(\mathbf{u},c\right)$ , we have

$$\begin{array}{cc}& J\left(\mathbf{u},c\right)={\mathit{\phi }}^T{\mathbf{A}}^T\mathbf{W}\mathbf{A}\mathbf{\phi }-2{\mathit{\phi }}^T{\mathbf{A}}^T\mathbf{W}\mathbf{h}+{\mathbf{h}}^T\mathbf{W}\mathbf{h}.\end{array}$$

(14)

Using the fact that A has full column rank and applying the Cauchy-Schwartz inequality [25] to the second term in the right-hand side of (14), we obtain

$$\begin{array}{cc}& J\left(\mathbf{u},c\right)\ge {\lambda }_{\min }{‖\mathit{\phi }‖}_2^2-2{‖\mathit{\phi }‖}_2{‖{\mathbf{A}}^T\mathbf{W}\mathbf{h}‖}_2+{\mathbf{h}}^T\mathbf{W}\mathbf{h},\end{array}$$

(15)

where ${\lambda }_{\min }$ is the smallest eigenvalue of ${\mathbf{A}}^T\mathbf{W}\mathbf{A}$ . From (15), it can be concluded that the set $\left\{\mathbf{u},c\mid {\lambda }_{\min }{‖\mathbf{\phi }‖}_{\mathrm{2}}^{\mathrm{2}}-\mathrm{2}{‖\mathbf{\phi }‖}_{\mathrm{2}}{‖{\mathbf{A}}^T\mathbf{W}\mathbf{h}‖}_{\mathrm{2}}+{\mathbf{h}}^T\mathbf{W}\mathbf{h}\le J\left(\mathbf{u},c\right)<d\right\}$ is bounded for any finite positive constant d. In addition, steps $(2)$ and $(3)$ lead to

$$\begin{array}{l}J\left({\mathbf{u}}_{k-1},c_{k-1}\right)\ge J\left({\mathbf{u}}_k,c_{k-1}\right)=\underset{\mathbf{u}}{\min }\hspace{0.17em}J\left(\mathbf{u},c_{k-1}\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\ge J\left({\mathbf{u}}_k,c_k\right)=\underset{c}{\min }\hspace{0.17em}J\left({\mathbf{u}}_k,c\right).\end{array}$$

(16)

Thus, for $k=\mathrm{1,2},\dots$ , the cost function reduces or stays the same. As a result, the cost function $J\left(\mathbf{u},c\right)$ is the so-called Lyapunov function [26], which indicates that the proposed bi-iterative method has asymptotic convergence.

3.3. Performance Analysis

We evaluate the bias and covariance matrix of the source position and propagation speed estimates $\widehat{\mathbf{u}}$ and $\widehat{c}$ through using the first-order perturbation method [27]. From (16), when the iteration procedure starting from appropriate initial source position and propagation speed estimates is terminated, the gradient of $J\left(\mathbf{u},c\right)$ at $\widehat{\mathbf{p}}={\left[{\widehat{\mathbf{u}}}^T,\widehat{c}\right]}^T$ satisfies

$$\begin{array}{cc}& {\left.\frac{\partial J\left(\mathbf{u},c\right)}{\partial \left(\mathbf{u},c\right)}\right|}_{\begin{array}{c}\mathbf{u}=\widehat{\mathbf{u}}\\ c=\widehat{c}\end{array}}={\left.{\mathbf{J}}^T\mathbf{W}\left(\mathbf{A}\mathbf{\phi }-\mathbf{h}\right)\right|}_{\begin{array}{c}\mathbf{u}=\widehat{\mathbf{u}}\\ c=\widehat{c}\end{array}}\approx {\mathbf{0}}_{\mathrm{3}\times \mathrm{1}},\end{array}$$

(17)

where $\mathbf{J}=\left[\mathbf{G}+c\tilde{\mathbf{t}}{\mathit{\rho }}_{\mathbf{u}\mathbf{,}{\mathbf{s}}_{\mathrm{1}}}^T,c\tilde{\mathbf{t}}\odot \tilde{\mathbf{t}}+r_{\mathrm{1}}\tilde{\mathbf{t}}\right]$ is the $\left(M-\mathrm{1}\right)\times \mathrm{3}$ Jacobian matrix and ${\mathit{\rho }}_{\mathbf{a},\mathbf{b}}={\left(\mathbf{a}-\mathbf{b}\right)/‖\left(\mathbf{a}-\mathbf{b}\right)‖}_{\mathrm{2}}$ denotes a unit vector. A, $\widehat{\mathbf{u}}$ , and $\widehat{c}$ can be expressed as $\mathbf{A}={\mathbf{A}}^o+\Delta \mathbf{A}$ , $\widehat{\mathbf{u}}=\mathbf{u}+\Delta \mathbf{u}\mathbf{,}$ and $\widehat{c}=c+\Delta c$ , respectively, where $\Delta \mathbf{u}$ and $\Delta c$ are the estimation bias in $\widehat{\mathbf{u}}$ and $\widehat{c}$ . ${\mathbf{A}}^o$ is obtained by replacing $\tilde{\mathbf{t}}$ with t in A, and $\Delta \mathbf{A}$ is the deviation of A from ${\mathbf{A}}^o$ .

Ignoring the perturbations that are of higher order than 1, we have

$$\begin{array}{l}\mathbf{A}\mathbf{\phi }-\mathbf{h}={\mathbf{A}}_{\mathrm{12}}^o\left(\mathbf{u}+\Delta \mathbf{u}-{\mathbf{s}}_{\mathrm{1}}\right)+\left({\mathbf{A}}_{\mathrm{3}}^o+\Delta {\mathbf{A}}_{\mathrm{3}}\right){\left(c+\Delta c\right)}^{\mathrm{2}}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\left({\mathbf{A}}_{\mathrm{4}}^o+\Delta {\mathbf{A}}_{\mathrm{4}}\right)\left(c+\Delta c\right){‖\mathbf{u}+\Delta \mathbf{u}-{\mathbf{s}}_{\mathrm{1}}‖}_{\mathrm{2}}-\mathbf{h}\\ \hspace{1em}\hspace{1em}\hspace{1em}\approx \left({\mathbf{A}}_{\mathrm{12}}^o+c{\mathbf{A}}_{\mathrm{4}}^o{\mathit{\rho }}_{\mathbf{u},{\mathbf{s}}_{\mathrm{1}}}^T\right)\Delta \mathbf{u}+\left(\mathrm{2}{\mathbf{A}}_{\mathrm{3}}^{o}c+{\mathbf{A}}_{\mathrm{4}}^o{r}_{\mathrm{1}}\right)\Delta c+\mathbf{B}\mathbf{n},\end{array}$$

(18)

where the fact that $\Delta {\mathbf{A}}_{\mathrm{12}}={\mathbf{0}}_{\left(M-\mathrm{1}\right)\times \mathrm{2}}$ is utilized. Assigning $\Delta \mathbf{p}={\left[\Delta {\mathbf{u}}^T,\Delta c\right]}^T$ and preserving only the linear perturbation terms, we obtain from (17)

$$\begin{array}{cc}& {\left.\frac{\partial J\left(\mathbf{u},c\right)}{\partial \left(\mathbf{u},c\right)}\right|}_{\begin{array}{c}\mathbf{u}=\widehat{\mathbf{u}}\\ c=\widehat{c}\end{array}}\approx {\mathbf{J}}^{oT}\mathbf{W}\left({\mathbf{J}}^o\Delta \mathbf{p}+\mathbf{B}\mathbf{n}\right)\approx {\mathbf{0}}_{\mathrm{3}\times \mathrm{1}},\end{array}$$

(19)

where ${\mathbf{J}}^o=\left[{\mathbf{A}}_{\mathrm{12}}^o+c{\mathbf{A}}_{\mathrm{4}}^o{\mathit{\rho }}_{\mathbf{u},{\mathbf{s}}_{\mathrm{1}}}^T,\mathrm{2}c{\mathbf{A}}_{\mathrm{3}}^o+r_1{\mathbf{A}}_{\mathrm{4}}^o\right]$ . Then the estimation bias $\Delta \mathbf{p}$ can be approximated as

$$\begin{array}{cc}& \Delta \mathbf{p}\approx -{\left({\mathbf{J}}^{oT}\mathbf{W}{\mathbf{J}}^o\right)}^{-\mathrm{1}}{\mathbf{J}}^{oT}\mathbf{W}\mathbf{B}\mathbf{n}.\end{array}$$

(20)

We have $E\left[\Delta \mathbf{p}\right]\approx {\mathbf{0}}_{\mathrm{3}\times \mathrm{1}}$ by using the assumption that n is a zero-mean random vector, which indicates that the source position and propagation speed estimates are approximately unbiased for low noise levels. Recalling the definition of W in (6) we have

$$\begin{array}{cc}& \mathrm{cov}\left(\widehat{\mathbf{p}}\right)=E\left[\Delta \mathbf{p}\Delta {\mathbf{p}}^T\right]\approx {\left({\mathbf{J}}^{oT}{\mathbf{B}}^{-T}{\mathbf{Q}}_t^{-\mathrm{1}}{\mathbf{B}}^{-\mathrm{1}}{\mathbf{J}}^o\right)}^{-\mathrm{1}}.\end{array}$$

(21)

Note that ${\left({\mathbf{J}}^{oT}{\mathbf{B}}^{-T}{\mathbf{Q}}_t^{-\mathrm{1}}{\mathbf{B}}^{-\mathrm{1}}{\mathbf{J}}^o\right)}^{-\mathrm{1}}$ in (21) is identical to the CRLB derived in [14], which means that the theoretical performance of the bi-iterative method closely meets the CRLB. The estimation bias and covariance matrix derived here are valid only if the TDOA measurement noises are relatively small.

4. Algorithm Extension

In practical environments, the sensor locations need to be measured and generally include random errors. In this section, we first examine the increase in the CRLB of the source position and propagation speed due to sensor location errors and then extend the bi-iterative method to the situation where the sensor locations are subject to random errors. Let ${\stackrel{-}{\mathbf{s}}}_j={\mathbf{s}}_j+\Delta {\mathbf{s}}_j$ denote the erroneous location of the jth sensor, where ${\mathbf{s}}_j$ is the true sensor location and $\Delta {\mathbf{s}}_j$ is the random location error, $j=\mathrm{1,2},\dots ,M$ . Arrange these sensor location errors $\Delta {\mathbf{s}}_{\mathrm{1}},\Delta {\mathbf{s}}_{\mathrm{2}},\dots ,\Delta {\mathbf{s}}_M$ into matrix form as $\Delta \mathbf{s}={\left[\Delta {\mathbf{s}}_{\mathrm{1}}^T,\Delta {\mathbf{s}}_{\mathrm{2}}^T,\dots ,\Delta {\mathbf{s}}_M^T\right]}^T$ . $\Delta \mathbf{s}$ is assumed to be Gaussian distributed with zero-mean and covariance matrix ${\mathbf{Q}}_{\mathbf{s}}=E\left[\Delta \mathbf{s}\Delta {\mathbf{s}}^T\right]$ and to be mutually independent of the measurement noise vector n.

4.1. CRLB

We obtain the CRLB of $\mathbf{p}={\left[{\mathbf{u}}^T,c\right]}^T$ with sensor location errors by following the same approach as in [20]:

$$\begin{array}{l}\text{CRLB}\left(\mathbf{p}\right)={\left(\mathbf{X}-\mathbf{Y}{\mathbf{Z}}^{-\mathrm{1}}{\mathbf{Y}}^T\right)}^{-\mathrm{1}}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}={\mathbf{X}}^{-\mathrm{1}}+{\mathbf{X}}^{-\mathrm{1}}\mathbf{Y}{\left(\mathbf{Z}-{\mathbf{Y}}^T{\mathbf{X}}^{-\mathrm{1}}\mathbf{Y}\right)}^{-\mathrm{1}}{\mathbf{Y}}^T{\mathbf{X}}^{-\mathrm{1}},\end{array}$$

(22)

where

$$\begin{array}{cc}& \mathbf{X}={\left(\frac{\partial \mathbf{t}}{\partial {\mathbf{p}}^T}\right)}^T{\mathbf{Q}}_t^{-\mathrm{1}}\left(\frac{\partial \mathbf{t}}{\partial {\mathbf{p}}^T}\right),\\ \mathbf{Y}={\left(\frac{\partial \mathbf{t}}{\partial {\mathbf{p}}^T}\right)}^T{\mathbf{Q}}_t^{-\mathrm{1}}\left(\frac{\partial \mathbf{t}}{\partial {\mathbf{s}}^T}\right),\\ \mathbf{Z}={\left(\frac{\partial \mathbf{t}}{\partial {\mathbf{s}}^T}\right)}^T{\mathbf{Q}}_t^{-\mathrm{1}}\left(\frac{\partial \mathbf{t}}{\partial {\mathbf{s}}^T}\right)+{\mathbf{Q}}_s^{-\mathrm{1}},\end{array}$$

(23)

where $\left(\partial \mathbf{t}/\partial {\mathbf{p}}^T\right)$ and $\left(\partial \mathbf{t}/\partial {\mathbf{s}}^T\right)$ are

$$\begin{array}{cc}& \frac{\partial \mathbf{t}}{\partial {\mathbf{p}}^T}=\left[\begin{array}{c}\frac{{\mathit{\rho }}_{\mathbf{u},{\mathbf{s}}_{\mathrm{2}}}^T}{c}-\frac{{\mathit{\rho }}_{\mathbf{u},{\mathbf{s}}_{\mathrm{1}}}^T}{c},-\frac{\left(r_2-r_1\right)}{c^2}\\ \frac{{\mathit{\rho }}_{\mathbf{u},{\mathbf{s}}_{\mathrm{3}}}^T}{c}-\frac{{\mathit{\rho }}_{\mathbf{u},{\mathbf{s}}_{\mathrm{1}}}^T}{c},-\frac{\left(r_3-r_1\right)}{c^2}\\ ⋮\\ \frac{{\mathit{\rho }}_{\mathbf{u},{\mathbf{s}}_M}^T}{c}-\frac{{\mathit{\rho }}_{\mathbf{u},{\mathbf{s}}_{\mathrm{1}}}^T}{c},-\frac{\left(r_M-r_1\right)}{c^2}\end{array}\right],\\ \frac{\partial \mathbf{t}}{\partial {\mathbf{s}}^T}=\frac{\mathrm{1}}{c}\left[\begin{array}{ccccc}{\mathit{\rho }}_{\mathbf{u},{\mathbf{s}}_{\mathrm{1}}}^T& -{\mathit{\rho }}_{\mathbf{u},{\mathbf{s}}_{\mathrm{2}}}^T& {\mathbf{0}}_{\mathrm{2}\times \mathrm{1}}^T& \cdots & {\mathbf{0}}_{\mathrm{2}\times \mathrm{1}}^T\\ {\mathit{\rho }}_{\mathbf{u},{\mathbf{s}}_{\mathrm{1}}}^T& {\mathbf{0}}_{\mathrm{2}\times \mathrm{1}}^T& -{\mathit{\rho }}_{\mathbf{u},{\mathbf{s}}_{\mathrm{3}}}^T& \cdots & {\mathbf{0}}_{\mathrm{2}\times \mathrm{1}}^T\\ ⋮& ⋮& & \ddots & \\ {\mathit{\rho }}_{\mathbf{u},{\mathbf{s}}_{\mathrm{1}}}^T& {\mathbf{0}}_{\mathrm{2}\times \mathrm{1}}^T& {\mathbf{0}}_{\mathrm{2}\times \mathrm{1}}^T& \cdots & -{\mathit{\rho }}_{\mathbf{u},{\mathbf{s}}_M}^T\end{array}\right].\end{array}$$

(24)

The term ${\mathbf{X}}^{-\mathrm{1}}$ in (22) is the CRLB of p when there exists no sensor location error [14]. The expression ${\mathbf{X}}^{-\mathrm{1}}\mathbf{Y}{\left(\mathbf{Z}-{\mathbf{Y}}^T{\mathbf{X}}^{-\mathrm{1}}\mathbf{Y}\right)}^{-\mathrm{1}}{\mathbf{Y}}^T{\mathbf{X}}^{-\mathrm{1}}$ in (22) stands for the increase in the CRLB of source position and propagation speed because of the presence of sensor location errors.

4.2. Algorithm Extension

The true sensor locations ${\mathbf{s}}_i$ and ${\mathbf{s}}_{\mathrm{1}}$ can be expressed as ${\mathbf{s}}_i={\stackrel{-}{\mathbf{s}}}_i-\Delta {\mathbf{s}}_i$ and ${\mathbf{s}}_{\mathrm{1}}={\stackrel{-}{\mathbf{s}}}_{\mathrm{1}}-\Delta {\mathbf{s}}_{\mathrm{1}}$ , respectively. Then (4) can be rewritten as

$$\begin{array}{l}{\left({\stackrel{-}{\mathbf{s}}}_i-{\stackrel{-}{\mathbf{s}}}_{\mathrm{1}}\right)}^T\left(\mathbf{u}-{\stackrel{-}{\mathbf{s}}}_{\mathrm{1}}\right)+\mathrm{0.5}{c}^{\mathrm{2}}{\tilde{t}}_{i\mathrm{1}}^{\mathrm{2}}+c{r}_1{\tilde{t}}_{i1}-\mathrm{0.5}{‖{\stackrel{-}{\mathbf{s}}}_i-{\stackrel{-}{\mathbf{s}}}_{\mathrm{1}}‖}_{\mathrm{2}}^{\mathrm{2}}\\ \hspace{1em}=c{r}_i{n}_{i1}-{\left(\mathbf{u}-{\stackrel{-}{\mathbf{s}}}_{\mathrm{1}}\right)}^T\Delta {\mathbf{s}}_{\mathrm{1}}+{\left(\mathbf{u}-{\stackrel{-}{\mathbf{s}}}_i\right)}^T\Delta {\mathbf{s}}_i,\end{array}$$

(25)

where the second-order error terms have been neglected.

It is noteworthy that the variable $r_1$ in (25) depends on the true sensor position ${\mathbf{s}}_{\mathrm{1}}$ . Hence, expand $r_1$ around the erroneous sensor location ${\stackrel{-}{\mathbf{s}}}_{\mathrm{1}}$ and preserve only the first two terms as

$$\begin{array}{cc}& r_1={\stackrel{-}{r}}_1+{\mathit{\rho }}_{\mathbf{u},{\stackrel{-}{\mathbf{s}}}_{\mathrm{1}}}^T\Delta {\mathbf{s}}_{\mathrm{1}},\end{array}$$

(26)

where ${\stackrel{-}{r}}_{\mathrm{1}}={‖\mathbf{u}-{\stackrel{-}{\mathbf{s}}}_{\mathrm{1}}‖}_{\mathrm{2}}$ . Inserting in (25) the expression of $r_1$ yields

$$\begin{array}{l}{\left({\stackrel{-}{\mathbf{s}}}_i-{\stackrel{-}{\mathbf{s}}}_{\mathrm{1}}\right)}^T\left(\mathbf{u}-{\stackrel{-}{\mathbf{s}}}_{\mathrm{1}}\right)+\mathrm{0.5}{\tilde{t}}_{i1}^2{c}^2+{\tilde{t}}_{i1}c{\stackrel{-}{r}}_1-\mathrm{0.5}{‖{\stackrel{-}{\mathbf{s}}}_i-{\stackrel{-}{\mathbf{s}}}_{\mathrm{1}}‖}_{\mathrm{2}}^{\mathrm{2}}\\ \hspace{1em}=c{r}_i{n}_{i1}-{\left(c{\tilde{t}}_{i1}{\mathit{\rho }}_{\mathbf{u},{\stackrel{-}{\mathbf{s}}}_{\mathrm{1}}}+\mathbf{u}-{\stackrel{-}{\mathbf{s}}}_{\mathrm{1}}\right)}^T\Delta {\mathbf{s}}_{\mathrm{1}}+{\left(\mathbf{u}-{\stackrel{-}{\mathbf{s}}}_i\right)}^T\Delta {\mathbf{s}}_i.\end{array}$$

(27)

Stacking (27) for $i=\mathrm{2,3},\dots ,M$ gives

$$\begin{array}{cc}& \stackrel{-}{\mathbf{A}}\stackrel{-}{\mathbf{\phi }}-\stackrel{-}{\mathbf{h}}=\stackrel{-}{\mathbf{\eta }},\end{array}$$

(28)

where $\stackrel{-}{\mathbf{A}}=\left[\stackrel{-}{\mathbf{G}},\mathrm{0.5}\tilde{\mathbf{t}}\odot \tilde{\mathbf{t}},\tilde{\mathbf{t}}\right]$ , $\stackrel{-}{\mathbf{G}}={\left[{\stackrel{-}{\mathbf{s}}}_{\mathrm{2}}-{\stackrel{-}{\mathbf{s}}}_{\mathrm{1}},{\stackrel{-}{\mathbf{s}}}_{\mathrm{3}}-{\stackrel{-}{\mathbf{s}}}_{\mathrm{1}},\dots ,{\stackrel{-}{\mathbf{s}}}_M-{\stackrel{-}{\mathbf{s}}}_{\mathrm{1}}\right]}^T$ , $\stackrel{-}{\mathbf{\phi }}={\left[{\left(\mathbf{u}-{\stackrel{-}{\mathbf{s}}}_{\mathrm{1}}\right)}^T,c^{\mathrm{2}},c{\stackrel{-}{r}}_1\right]}^T$ , and $\stackrel{-}{\mathbf{h}}=\mathrm{0.5}{\left[{‖{\stackrel{-}{\mathbf{s}}}_{\mathrm{2}}-{\stackrel{-}{\mathbf{s}}}_{\mathrm{1}}‖}_{\mathrm{2}}^{\mathrm{2}},{‖{\stackrel{-}{\mathbf{s}}}_{\mathrm{3}}-{\stackrel{-}{\mathbf{s}}}_{\mathrm{1}}‖}_{\mathrm{2}}^{\mathrm{2}},\dots ,{‖{\stackrel{-}{\mathbf{s}}}_M-{\stackrel{-}{\mathbf{s}}}_{\mathrm{1}}‖}_{\mathrm{2}}^{\mathrm{2}}\right]}^T$ . On the right-hand side of (28), $\stackrel{-}{\mathbf{\eta }}=\mathbf{B}\mathbf{n}+\mathbf{D}\Delta \mathbf{s}$ is the equation error vector, and the $\left(i-\mathrm{1}\right)$ th row, $i=\mathrm{2,3},\dots ,M$ , of D is

$$\begin{array}{cc}& \left[-{\left(c{\tilde{t}}_{i1}{\mathit{\rho }}_{\mathbf{u},{\stackrel{-}{\mathbf{s}}}_{\mathrm{1}}}+\mathbf{u}-{\stackrel{-}{\mathbf{s}}}_{\mathrm{1}}\right)}^T,{\mathbf{0}}_{\mathrm{2}(i-\mathrm{2})\times \mathrm{1}}^T,{\left(\mathbf{u}-{\stackrel{-}{\mathbf{s}}}_i\right)}^T,{\mathbf{0}}_{\mathrm{2}(M-i)\times \mathrm{1}}^T\right].\end{array}$$

(29)

Then we have the following cost function:

$$\begin{array}{cc}& \stackrel{-}{J}\left(\mathbf{u},c\right)={\stackrel{-}{\mathbf{\eta }}}^T\stackrel{-}{\mathbf{W}}\stackrel{-}{\mathbf{\eta }}={\left(\stackrel{-}{\mathbf{A}}\stackrel{-}{\mathbf{\phi }}-\stackrel{-}{\mathbf{h}}\right)}^T\stackrel{-}{\mathbf{W}}\left(\stackrel{-}{\mathbf{A}}\stackrel{-}{\mathbf{\phi }}-\stackrel{-}{\mathbf{h}}\right),\end{array}$$

(30)

where the weighting matrix is defined as $\stackrel{-}{\mathbf{W}}={\left(\mathbf{B}{\mathbf{Q}}_t{\mathbf{B}}^T+\mathbf{D}{\mathbf{Q}}_s{\mathbf{D}}^T\right)}^{-\mathrm{1}}$ .

Once the cost function $\stackrel{-}{J}\left(\mathbf{u},c\right)$ is obtained, we can use the same bi-iterative scheme presented in Section 3.1 to compute the source position and propagation speed.

5. Simulation Study

In this section, a set of Monte Carlo simulations are carried out to assess the bi-iterative method's performance by comparing with the LS solution [13], the four-step method [14], the three-step WLS method [15], the CRLB, and the theoretical bias [28]. We simply follow [15] to consider a tangible acoustic interface application of interactive displays [29]. Five sensors are placed on a 1 m × 1 m glass pane with coordinates ${\mathbf{s}}_{\mathrm{1}}={\left[\mathrm{0.5,0.5}\right]}^T$  m, ${\mathbf{s}}_{\mathrm{2}}={\left[\mathrm{0,0}\right]}^T$  m, ${\mathbf{s}}_{\mathrm{3}}={\left[\mathrm{1,0}\right]}^T$  m, ${\mathbf{s}}_{\mathrm{4}}={\left[\mathrm{1,1}\right]}^T$  m, and ${\mathbf{s}}_{\mathrm{5}}={\left[\mathrm{0,1}\right]}^T$  m, respectively; the speed of propagation is $\mathrm{1200}$  m/s; they are the same as those in [15]. The covariance matrix of TDOA measurement noise vector n is ${\mathbf{Q}}_t=\mathrm{0.5}{\sigma }_t^2\left({\mathbf{I}}_{M-\mathrm{1}}+{\mathbf{1}}_{M-\mathrm{1}}{\mathbf{1}}_{M-\mathrm{1}}^T\right)$ [6], where ${\mathbf{I}}_{M-\mathrm{1}}$ is an identity matrix of size $M-\mathrm{1}$ . In the four-step method, we use the well-known two-stage WLS technique [6] to compute the source location. The performance is assessed by root mean square error (RMSE) and norm of estimation bias, which are defined as $\sqrt{{\sum }_{l=1}^L{‖{\widehat{\mathit{\lambda }}}^{(l)}-\mathit{\lambda }‖}_2^2/L}$ and ${‖{\sum }_{l=1}^L\left({\widehat{\mathit{\lambda }}}^{(l)}-\mathit{\lambda }\right)/L‖}_{\mathrm{2}}$ , respectively, where ${\widehat{\mathbf{\lambda }}}^{(l)}$ is the estimate of $\mathbf{\lambda }$ at the lth trial and $\mathbf{\lambda }$ is either the source position u or propagation speed c. $L=\mathrm{2000}$ is the total number of Monte Carlo trials.

5.1. Localization Performance without Sensor Location Errors

We first examine the performance of the bi-iterative estimation method when there exists no sensor location error. Figure 1 plots the estimation RMSEs versus ${\sigma }_t$ when the source is at $\mathbf{u}={\left[\mathrm{0.2,0.1}\right]}^T$  m. Remarkably, the LS solution cannot achieve the CRLB performance even when ${\sigma }_t$ is as small as $\mathrm{1}{\mathrm{0}}^{-\mathrm{6}}$ . The remaining three methods closely approach the CRLB for both source position and propagation speed estimates when ${\sigma }_t\le \mathrm{4.6}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{6}}$ , which verifies the analysis in Section 3.3 on the approximate efficiency of the bi-iterative method. As ${\sigma }_t$ further increases, the bi-iterative, three-step WLS and four-step methods deviate from the CRLB. This is the so-called threshold phenomenon, which occurs in nonlinear estimation when the noises exceed a threshold value [8]. Note that the estimation RMSEs of the bi-iterative method are below the CRLB for high noise levels. A reasonable explanation for this is that the bi-iterative method gives biased estimates for high noise levels. It is possible that when an estimate is biased, its RMSE can be smaller than the CRLB [21].

Figure 1

Comparison of the source position and propagation speed estimation RMSEs of the LS solution, the three-step WLS, and four-step and bi-iterative methods versus ${\sigma }_t$ for source at $\mathbf{u}={\left[\mathrm{0.2,0.1}\right]}^T$  m.

Figure 2 compares the estimation biases of the bi-iterative method with the remaining three methods and the theoretical bias [28] for source at $\mathbf{u}={\left[\mathrm{0.2,0.1}\right]}^T$  m. As expected, the source position and propagation speed estimation biases of the bi-iterative method are all reasonably small for low noise levels, which illustrates the unbiasedness of the bi-iterative method for small measurement noises. As ${\sigma }_t$ increases, the estimation biases of the bi-iterative method deviate from the theoretical bias apparently. The bi-iterative method is superior to the remaining three methods.

Figure 2

Comparison of the source position and propagation speed estimation biases of the LS solution, the three-step WLS, and four-step and bi-iterative methods versus ${\sigma }_t$ for source at $\mathbf{u}={\left[\mathrm{0.2,0.1}\right]}^T$  m.

Figures 3 and 4 show the estimation RMSE results for two sources at $\mathbf{u}={\left[\mathrm{0.5,0.1}\right]}^T$  m and $\mathbf{u}={\left[\mathrm{0.2,0.5}\right]}^T$  m, respectively. Each of these two sources has one coordinate being identical to that of the reference sensor. More specifically, the x-coordinate of ${\left[\mathrm{0.5,0.1}\right]}^T$ and the y-coordinate of ${\left[\mathrm{0.2,0.5}\right]}^T$ are equivalent to the corresponding coordinates of ${\mathbf{s}}_{\mathrm{1}}$ . Clearly, the estimation RMSEs of the bi-iterative and four-step methods closely meet the CRLB for ${\sigma }_t\le \mathrm{1}{\mathrm{0}}^{-\mathrm{4}}$ , while those of the LS and three-step WLS methods are significantly larger than the CRLB. The most likely reason for this is that a matrix ill-conditioned problem occurs in computing the weighting matrix, which indicates that the LS and three-step WLS methods are more sensitive to the source positions.

Figure 3

Comparison of the source position and propagation speed estimation RMSEs of the LS solution, the three-step WLS, and four-step and bi-iterative methods versus ${\sigma }_t$ for source at $\mathbf{u}={\left[\mathrm{0.5,0.1}\right]}^T$  m.

Figure 4

Comparison of the source position and propagation speed estimation RMSEs of the LS solution, the three-step WLS, and four-step and bi-iterative methods versus ${\sigma }_t$ for source at $\mathbf{u}={\left[\mathrm{0.2,0.5}\right]}^T$  m.

Figures 5 and 6 plot the estimation RMSE results of the bi-iterative and four-step methods versus the iteration number for sources at $\mathbf{u}={\left[\mathrm{0.5,0.1}\right]}^T$  m and $\mathbf{u}={\left[\mathrm{0.2,0.5}\right]}^T$  m, respectively. Initial values for the bi-iterative and four-step WLS methods are the LS solutions. It is clearly demonstrated that both the bi-iterative and four-step methods can converge very well in about four iterations.

Figure 5

Comparison of the source position and propagation speed estimation RMSEs of the bi-iterative and four-step methods versus iteration number for source at $\mathbf{u}={\left[\mathrm{0.5,0.1}\right]}^T$  m.

Figure 6

Comparison of the source position and propagation speed estimation RMSEs of the bi-iterative and four-step methods versus iteration number for source at $\mathbf{u}={\left[\mathrm{0.2,0.5}\right]}^T$  m.

In the following, we evaluate the computational complexity of the bi-iterative method. In step $(2)$ , the complexity of calculating propagation speed c is about $O\left(2Mn+6M^2\right)$ , where M is the number of sensors and n is the dimension of the source position. In step $(3)$ , the complexity induced by calculating F is roughly $O\left(2\left(M^{2}n+M{n}^2\right)\right)$ . Denote by u and v the number of equality constraints and the optimization problem size, respectively; the worst case computational complexity of solving SDP is $O\left(\left(u^3+u^2{v}^2+u{v}^3\right)v^{0.5}\right)$ [30]. In the SDP (13), $u=\mathrm{2}$ and $v=n+\mathrm{2}$ . Therefore, the total complexity is about $O\left(2\left(M^{2}n+M{n}^2\right)+\left(u^3+u^2{v}^2+u{v}^3\right)v^{0.5}\right)$ . The complexity is $O\left(2\left(M^{2}n+M{n}^2\right)+\left(6n+7\right)n^2\right)$ for the four-step method. It is clear that the bi-iterative and four-step methods have comparable computational complexity since we have $n=\mathrm{2}$ or $\mathrm{3}$ .

5.2. Localization Performance with Sensor Position Errors

Here, we investigate the CRLB of source position and propagation speed and evaluate the bi-iterative method's performance with sensor location errors. The sensor location errors at different sensors are set to be unequal to make the comparison more general, where ${\mathbf{Q}}_s={\sigma }_s^{\mathrm{2}}\mathrm{diag}\left(\left[\mathrm{4,4},\mathrm{2,2},\mathrm{8,8},\mathrm{20,20,10,10}\right]\right)$ . Fix ${\sigma }_t$ at $\mathrm{1}{\mathrm{0}}^{-\mathrm{5}}$ and Figure 7 depicts the trace of the upper-left $\mathrm{2}\times \mathrm{2}$ submatrix and the third diagonal element of $\text{CRLB}\left(\mathbf{p}\right)$ as ${\sigma }_s$ increases for source at $\mathbf{u}={\left[\mathrm{0.2,0.1}\right]}^T$  m. The CRLBs of source position and propagation speed without sensor location errors are also plotted for comparison purpose. It is clear that the deviations of the CRLB of source position and propagation speed with sensor location errors from the case without sensor location errors become larger and larger as ${\sigma }_s$ increases.

Figure 7

Comparison of the CRLB of the source position and propagation speed with and without sensor location errors versus ${\sigma }_s$ for source at $\mathbf{u}={\left[\mathrm{0.2,0.1}\right]}^T$  m.

We set ${\sigma }_s$ to be $\mathrm{1}{\mathrm{0}}^{-\mathrm{2}}$ and check the influence of TDOA measurement noise on the estimation accuracy in the presence of sensor location errors. Figure 8 plots the RMSE results of the bi-iterative and LS methods as ${\sigma }_t$ increases for source at $\mathbf{u}={\left[\mathrm{0.2,0.1}\right]}^T$  m. Again, the LS solution cannot achieve the CRLB even when ${\sigma }_t=\mathrm{1}{\mathrm{0}}^{-\mathrm{6}}$ , while the bi-iterative method approaches the CRLB until it suffers from the thresholding effect at ${\sigma }_t=\mathrm{4.3}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{5}}$ . A comparison of Figures 1 and 8 reveals that a small error in the sensor locations can significantly degrade the estimation accuracy when the TDOA measurement noise is relatively small.

Figure 8

Comparison of the source position and propagation speed estimation RMSEs of the bi-iterative method and the LS solution versus ${\sigma }_t$ for source at $\mathbf{u}={\left[\mathrm{0.2,0.1}\right]}^T$  m.

We next keep ${\sigma }_t$ at $\mathrm{1}{\mathrm{0}}^{-\mathrm{5}}$ and investigate the effect of ${\sigma }_s$ on the source position and propagation speed estimation accuracy. Figure 9 plots the RMSE results against ${\sigma }_s$ for the source at $\mathbf{u}={\left[\mathrm{0.2,0.1}\right]}^T$  m. Clearly, neither the source position nor propagation speed estimation accuracy in the LS solution can achieve the CRLB performance even when ${\sigma }_s$ is as small as $\mathrm{1}{\mathrm{0}}^{-\mathrm{3}}$ , while the bi-iterative method closely approximates the CRLB when ${\sigma }_s$ is less than $\mathrm{3.6}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{2}}$ .

Figure 9

Comparison of the source position and propagation speed estimation RMSEs of the bi-iterative method and the LS solution versus ${\sigma }_s$ for source at $\mathbf{u}={\left[\mathrm{0.2,0.1}\right]}^T$  m.

6. Conclusions

The problem of estimating the source location and the propagation speed using TDOA measurements was considered in this paper. Using the bi-iterative scheme, an efficient method that alternately computes the source position and propagation speed is developed. The asymptotical convergence of the proposed bi-iterative estimation method is analyzed. The proposed method is shown by first-order perturbation analysis to approach the CRLB for low noise levels. In contrast to other existing estimation methods, the threshold effect of the new bi-iterative method occurs later as the measurement noise increases. Simulation studies are given to illustrate the effectiveness of the proposed bi-iterative method.