Constrained total least squares localization using angle of arrival and time difference of arrival measurements in the presence of synchronization clock bias and sensor position errors

Abstract

Based on measurements of angle of arrival and time difference of arrival, a method is proposed to improve the accuracy of localization with imperfect sensors. A derivation of the Cramér–Rao lower bound and the root mean square error is presented aimed at demonstrating the significance of taking synchronization errors into consideration. Subsequently, a set of pseudo-linear equations are constructed, based on which the constrained total least squares optimization model has been formulated for target localization and the Newton iteration is applied to obtain the source position and clock bias simultaneously. The theoretical performance of the constrained total least squares localization algorithm subject to sensor position errors and synchronization clock bias is derived, and a framework for the performance analysis is developed. In addition, the first-order error analysis illustrates that the proposed method can achieve the Cramér–Rao lower bound under moderate Gaussian noises by a mathematic derivation. Finally, simulation results are presented that verify the validity of the theoretical derivation and superiority of the new algorithm.

Introduction

Passive source localization based on the angle of arrival (AOA) and the time difference of arrival (TDOA) has been of considerable interest as a research topic in regard to sensor networks in the field of signal processing and for its broad applications in target tracking, surveillance, and navigation. ¹ This kind of hybrid localization not only yields a better accuracy than using TDOA alone ² but also overcomes the drawback that the accuracy is sensitive to the baseline between target and sensors when using only AOA localization.

A constrained weighted least squares (CWLS) positioning approach using AOA and TDOA measurements has been presented in the study by Cheung et al., ³ which achieves zero bias and is approximately within the Cramér–Rao lower bound (CRLB) under small measurement noise. A new AOA- and TDOA-based localization algorithm was formulated in the study by Bishop et al., ⁴ which exploits a geometrically derived constraint on the measurement errors in estimating the maximum likelihood values of errors given an appropriate cost function. AOA and TDOA measurements employing a single sensor have been jointly implemented thereafter in localizing radar emitters without a priori knowledge of the reflecting environment. ⁵ By exploiting the hybrid TDOA and AOA measurements, a closed-form solution has been presented in the study by Yin et al. ⁶ for the three-dimensional (3D) source localization problem using only two observation stations.

However, most existing AOA and TDOA localization algorithms assume that the sensor position is actually known, which is very unlikely to be satisfied in practice because of the possible variations in or inaccurate knowledge of sensor positions. ⁷ For instance, the sensors employed to receive target signals are fixed in planes or vessels, and their positions need to be estimated by self-positioning systems, from which true positions are compromised by either system errors or sensor position errors. The localization accuracy has been shown in the study by Ho et al. ⁸ to be sensitive to accurate knowledge of the sensor position and a slight error will create a sharp degradation in localization accuracies. Hence, the uncertainty associated with sensor positions needs to be taken into account in AOA- and TDOA-based localization algorithms to improve localization accuracy.

Another challenge for AOA and TDOA localization is the strict synchronization between sensors and the reference. To obtain high accuracy in AOA estimates, the sensors are always distributed far apart. Nevertheless, this is impractical if synchronized sampling needs to be obtained. In addition, a large number of sensor nodes and their non-ideal clocks will lead to a non-synchronization of sensors and result in time offsets. As a result, another system error, the so-called synchronization error occurring in TDOA measurements, can degrade the localization accuracy considerably. ⁹

As mentioned earlier, considering that strict synchronization of a sensor’s clock and precise positioning of sensors are not available in practice, both system errors should be taken into account together. The localization scenario considered in the study by Wang and Ho ¹⁰ consisted of a set of uncorrelated sensors with timing synchronization offsets and random position errors existing among different sensor groups, in which sensors are synchronized with each other. To some degree, this method can eliminate the synchronization errors, but fails in determining source position due to an insufficient number of sensors or a number of partially synchronized sensors. Sun and Yang ¹¹ derived the localization mean square error (MSE) for the case where the source clock bias is absent and achieved time synchronization and source localization. However, this article used only time of arrival (TOA) measurements, and a new method utilizing multiple measurements should be exploited. This article investigates an AOA and TDOA localization method that handles both sensor position errors and clock bias and is not restricted to partially synchronized sensors. Hence, it has universal applicability in a number of scenarios.

It is well-known that the total least squares (TLS) algorithm can overcome the influence of noise in the observation matrix when estimates of the source position are obtained using the linear localization model. ¹² When noise is independent and identically distributed, this method can achieve an asymptotic optimal performance. Nevertheless, because pseudo linearization in source localization is used, noise hidden in the observation matrix is in the structured components and subsequent statistical correlations. There is no doubt that the TLS method performs poorly with high noise levels and cannot approach the CRLB.

Fortunately, inspired by the CWLS algorithm proposed in the study by Markovsky and Huffel, ¹² which exploits known relationships between intermediate variables and the source location coordinates, the constrained total least squares (CTLS) algorithm extends the TLS to solve the statistical correlations in noise from the viewpoint of their numerical structural information. Yu et al. ¹³ and Wang et al. ¹⁴ proposed CTLS methods based on measurements that were only AOA and only TDOA, respectively. The corresponding performance curves attained CRLB and possessed a higher threshold compared with the closed-form method. Wang extended the CTLS method to localization scenarios based on an external illuminator. ¹⁵ Subsequently, this method was implemented to determine the source position with sensor location uncertainties using only TDOA measurements in the study by Yang et al. ¹⁶ and then joint TDOA/frequency difference of arrival (FDOA) measurements in the study by Chen et al. ¹⁷ Hence, this article develops a CTLS method based on a hybrid AOA and TDOA localization considering both sensor position errors and synchronization errors.

The study starts by deriving the CRLB for errors in sensor position and synchronization and then exhibits the root mean square error (RMSE) of a source position estimate acquired from the optimal estimator but not accounting for synchronization errors that actually exist. A comparison of the CRLB and RMSE significantly illustrates the need to find a novel method taking synchronization into account to perfectly estimate performance. The theoretical analysis forms the framework for AOA and TDOA localization in the presence of a non-ideal factor with good universality.

The study proceeds with the CTLS method based on AOA and TDOA measurements. We first transform the non-linear algebraic equation into a pseudo-linear equation depending on the geometry of the source and sensors. Following this, we construct the CLTS optimum localization model and apply the iterative Newton’s method to solve this problem and obtain the source position. Moreover, this article develops a unified theoretical framework in which the CTLS localization algorithm performance is derived and strictly proves that the theoretical performance reaches the CRLB to first order in error analysis. The theoretical framework of this analysis is not limited to a specific observation equation and can extend to more general localization scenarios. Therefore, to examine the robustness and practicality of the CTLS localization algorithm, the algorithms and theoretical analysis are developed with errors modeled and finally tested in simulations to verify positioning performance.

The rest of article is organized as follows. The pseudo-linear formulation based on AOA and TDOA measurements is described in section “Measurement model and problem formulation.” In section “CRLB and RMSE,” the CRLB with errors modeled and the RMSE assuming no synchronization errors are developed in detail. Section “CTLS localization method” provides a system of derivation of the CTLS method and its theoretical performance. Simulation results are presented in section “Simulations,” and conclusions are given in section “Conclusion.” The proofs and details of the main results are shown in Appendixes 1–3.

We use the following notations throughout the article: (1) boldface lower-case letters and boldface upper-case letters denote vectors and matrices, respectively; (2) ${\mathit{O}}_{N\times n}$ denotes the $N\times n$ matrix with all-zero entries, and similarly ${\mathit{0}}_n$ a n-dimensional zero vector. (3) ${{\mathit{1}}_M^{\left(1\right)}}^{\mathrm{T}}$ denotes the $M$ column vector of zeros except the ith element is 1, and ${\mathit{I}}_N$ denotes the N-dimensional identity matrix. (4) · represents the Hadamard product (multiplied by the corresponding elements); (5) ⊗ represents the Kronecker product; (6) ${\left[·\right]}^{†}$ represents the Moore–Penrose inverse of a matrix; (7) ${〈·〉}_n$ represents the $n$ elements of a matrix; (8) $\mathrm{diag}[·]$ represents a diagonal matrix equipped with the elements of a vector; and (9) $\mathrm{blkdiag}[·]$ represents a block diagonal matrix consisting of matrices and/or vectors.

Measurement model and problem formulation

Consider a 3D scenario composed of $M$ sensors each of which is capable of determining both AOA and TDOA of a signal source to jointly determine the unknown radiation source located at $u^0=[x^0,y^0,z^0{]}^T$ . Figure 1 shows that the localization geometry, and the unit for the positions is meter. The ith sensor whose true position denoted by ${\mathit{s}}_i^o=[x_i^o,y_i^o,z_i^o{]}^{\mathrm{T}},i=1,2,\dots ,M$ are assumed to be known and provides the AOA measurements consisting of the azimuth angle ${\theta }_i^o$ and elevation angle ${\beta }_i^o$ , whose relationship with the target and sensor i can be described as

$$\{\begin{array}{c}{\theta }_i^o=\arctan \frac{y^o-y_i^o}{x^o-x_i^o}\hfill \\ \multicolumn{1}{c}{{\beta }_i^o=\arctan \frac{z^o-z_i^o}{\sqrt{{(x^o-x_i^o)}^2+{(y^o-y_i^o)}^2}}}\end{array},i=1,2,\dots ,M$$

(1)

Figure 1.

3D localization scenario with a source and one sensor to show the geometry relationship of the AOA/TDOA measurements and source position when the sensor position is uncertain.

As sensors are not synchronized, the constant clock bias existing between the sensor pair $i$ and 1 (the reference sensor) is ${\tau }_{i1}^o$ , which is so-called synchronization error in nanoseconds. Hence, the sets of equations relating the range difference $r_{i1}^o$ with target and sensor positions are given by

$$r_{i1}^o=\Vert {\mathit{u}}^o-{\mathit{s}}_i^o\Vert -\Vert {\mathit{u}}^o-{\mathit{s}}_1^o\Vert +{\delta }_{i1}^o,i=2,\dots ,M$$

(2)

where ${\delta }_{i1}^o=c{\tau }_{i1}^o$ and $c$ is the signal propagation speed in meter per second.

We define the vectors ${\mathit{\theta }}^o=[{\theta }_1^o,{\theta }_2^o,\dots ,{\theta }_M^o{]}^{\mathrm{T}}$ , ${\mathit{\beta }}^o=[{\beta }_1^o,{\beta }_2^o,\dots ,{\beta }_M^o{]}^{\mathrm{T}}$ , and ${\mathit{r}}^{\mathit{o}}=[r_{21}^o,r_{31}^o,\dots ,r_{M1}^o{]}^{\mathrm{T}}$ are the pure angle and range difference measurements, which have not been contaminated by noises. Besides, ${\mathit{\delta }}^{\mathit{o}}=[{\delta }_{21}^o,{\delta }_{31}^o,\dots ,{\delta }_{M1}^o{]}^{\mathrm{T}}$ is a synchronization error vector.

With the aid of trigonometric identity that $\tan {\theta }_i^o=\sin {\theta }_i^o/\cos {\theta }_i^o$ in Figure 1, equation (1) can be formatted as the pseudo-linear equations

$$\{\begin{array}{c}(x^o-x_i^o)\sin \left({\theta }_i^o\right)=(y^o-y_i^o)\cos \left({\theta }_i^o\right)\hfill \\ \multicolumn{1}{c}{(z^o-z_i^o)\cos \left({\beta }_i^o\right)=((x^o-x_i)\cos \left({\theta }_i^o\right)+(y^o-y_i^o)\sin \left({\theta }_i^o\right))\sin \left({\beta }_i^o\right)}\end{array},i=1,2,\dots ,M$$

(3)

Furthermore, rewrite equation (3) into vector format

$$\{\begin{array}{c}{\mathit{\omega }}_1^{\mathrm{T}}\left({\theta }_i^o\right){\mathit{u}}^o={\mathit{\omega }}_1^{\mathrm{T}}\left({\theta }_i^o\right){\mathit{s}}_i^o\hfill \\ \multicolumn{1}{c}{{\mathit{\omega }}_2^{\mathrm{T}}({\theta }_i^o,{\beta }_i^o){\mathit{u}}^o={\mathit{\omega }}_2^{\mathrm{T}}({\theta }_i^o,{\beta }_i^o){\mathit{s}}_i^o}\end{array},i=1,2,\dots ,M$$

(4)

where

$$\{\begin{array}{c}{\mathit{\omega }}_1\left({\theta }_i^o\right)={[\sin \left({\theta }_i^o\right),-\cos \left({\theta }_i^o\right),0]}^{\mathrm{T}}\hfill \\ \multicolumn{1}{c}{{\mathit{\omega }}_2({\theta }_i^o,{\beta }_i^o)={[\cos \left({\theta }_i^o\right)\sin \left({\beta }_i^o\right),\sin \left({\theta }_i^o\right)\sin \left({\beta }_i^o\right),-\cos \left({\beta }_i^o\right)]}^{\mathrm{T}}}\end{array}$$

(5)

Different from the AOA equations, the TDOA equations cannot be transformed into the pseudo-linear equations directly. However, it can be solved validly in coordination with angle measurements. We resort to the geometry relationship of target and sensors as follows

$$\Vert {\mathit{u}}^o-{\mathit{s}}_i^o\Vert ={\mathit{\omega }}_3^{\mathrm{T}}({\theta }_i^o,{\beta }_i^o)({\mathit{u}}^o-{\mathit{s}}_i^o)$$

(6)

where

$${\mathit{\omega }}_3({\theta }_i^o,{\beta }_i^o)={[\cos \left({\theta }_i^o\right)\cos \left({\beta }_i^o\right),\sin \left({\theta }_i^o\right)\cos \left({\beta }_i^o\right),\sin \left({\beta }_i^o\right)]}^{\mathrm{T}}$$

(7)

Substitute equation (6) into equation (2) and obtain the TDOA equations

$$\begin{array}{c}{r}_{i1}^o+{\mathit{\omega }}_3^{\mathrm{T}}({\theta }_i^o,{\beta }_i^o){\mathit{s}}_i^o-{\mathit{\omega }}_3^{\mathrm{T}}({\theta }_1^o,{\beta }_1^o){\mathit{s}}_1^o\hfill \\ =({\mathit{\omega }}_3^{\mathrm{T}}({\theta }_i^o,{\beta }_i^o)-{\mathit{\omega }}_3^{\mathrm{T}}({\theta }_1^o,{\beta }_1^o)){\mathit{u}}^o+{\delta }_{i1}^o,i=2,\dots ,M\hfill \end{array}$$

(8)

We define the measurement vector ${\mathit{\xi }}^{\mathit{o}}=[{\mathit{\theta }}^{\mathit{o}\mathrm{T}},{\mathit{\beta }}^{\mathit{o}\mathrm{T}},{\mathit{r}}^{\mathit{o}\mathrm{T}}{]}^{\mathrm{T}}$ and sensor position vector ${\mathit{s}}^{\mathit{o}}=[{\mathit{s}}_1^{\mathrm{T}},{\mathit{s}}_2^{\mathrm{T}},\dots ,{\mathit{s}}_M^{\mathrm{T}}{]}^{\mathrm{T}}$ and note the true estimated vector ${\mathit{\mu }}^o=[{{\mathit{u}}^o}^{\mathrm{T}},{{\mathit{\delta }}^o}^{\mathrm{T}}{]}^{\mathrm{T}}$ . Then simultaneous equations (4) and (8) give the joint localization equation set based on AOA and TDOA measurements as follows

$$\mathit{A}\left({\mathit{\xi }}^{\mathit{o}}\right){\mathit{\mu }}^o=\mathit{b}({\mathit{\xi }}^{\mathit{o}},{\mathit{s}}^{\mathit{o}})$$

(9)

where

$$\mathit{A}\left({\mathit{\xi }}^{\mathit{o}}\right)={\left[\begin{array}{cc}{\mathit{\omega }}_1^{\mathrm{T}}\left({\theta }_1^o\right)\hfill & {\mathit{0}}_{M-1}^{\mathrm{T}}\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill \\ \multicolumn{1}{c}{{\mathit{\omega }}_1^{\mathrm{T}}\left({\theta }_M^o\right)}& {\mathit{0}}_{M-1}^{\mathrm{T}}\hfill \\ \multicolumn{1}{c}{{\mathit{\omega }}_2^{\mathrm{T}}({\theta }_1^o,{\beta }_1^o)}& {\mathit{0}}_{M-1}^{\mathrm{T}}\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill \\ \multicolumn{1}{c}{{\mathit{\omega }}_2^{\mathrm{T}}({\theta }_M^o,{\beta }_M^o)}& {\mathit{0}}_{M-1}^{\mathrm{T}}\hfill \\ \multicolumn{1}{c}{({\mathit{\omega }}_3^{\mathrm{T}}({\theta }_2^o,{\beta }_2^o)-{\mathit{\omega }}_3^{\mathrm{T}}({\theta }_1^o,{\beta }_1^o))}& {{\mathit{1}}_{M-1}^{\left(1\right)}}^{\mathrm{T}}\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill \\ \multicolumn{1}{c}{({\mathit{\omega }}_3^{\mathrm{T}}({\theta }_M^o,{\beta }_M^o)-{\mathit{\omega }}_3^{\mathrm{T}}({\theta }_1^o,{\beta }_1^o))}& {{\mathit{1}}_{M-1}^{(M-1)}}^{\mathrm{T}}\hfill \end{array}\right]}_{3M-1\times K+M-1}$$

(10)

$$\mathit{b}({\mathit{\xi }}^{\mathit{o}},{\mathit{s}}^{\mathit{o}})=\left[\begin{array}{c}{\mathit{\omega }}_1^{\mathrm{T}}\left({\theta }_1^o\right){\mathit{s}}_1^o\hfill \\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{{\mathit{\omega }}_1^{\mathrm{T}}\left({\theta }_M^o\right){\mathit{s}}_M^o}\\ \multicolumn{1}{c}{{\mathit{\omega }}_2^{\mathrm{T}}({\theta }_1^o,{\beta }_1^o){\mathit{s}}_1^o}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{{\mathit{\omega }}_2^{\mathrm{T}}({\theta }_M^o,{\beta }_M^o){\mathit{s}}_M^o}\\ \multicolumn{1}{c}{r_{21}^o+{\mathit{\omega }}_3^{\mathrm{T}}({\theta }_2^o,{\beta }_2^o){\mathit{s}}_2^o-{\mathit{\omega }}_3^{\mathrm{T}}({\theta }_1^o,{\beta }_1^o){\mathit{s}}_1^o}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{r_{M1}^o+{\mathit{\omega }}_3^{\mathrm{T}}({\theta }_M^o,{\beta }_M^o){\mathit{s}}_M^o-{\mathit{\omega }}_3^{\mathrm{T}}({\theta }_1^o,{\beta }_1^o){\mathit{s}}_1^o}\end{array}\right]$$

(11)

We would like to illustrate the physical significance of the aforementioned algebraic expressions that $\mathit{A}({\mathit{\xi }}^{\mathit{o}})$ represents the coefficient matrix, $\mathit{b}({\mathit{\xi }}^{\mathit{o}},{\mathit{s}}^{\mathit{o}})$ represents the pseudo-linear measurement vector.

In practice, ${\mathit{\xi }}^{\mathit{o}}$ derives from the coherent processing of a radio signal, which will be interfered by the noises inevitably. Hence, we assume that the AOA and TDOA measurements can be described by the additive noise model as

$$\mathit{\xi }={\mathit{\xi }}^{\mathit{o}}+{\mathit{e}}_1$$

(12)

where ${\mathit{e}}_1=[e_{11},\dots ,e_{1,3M-1}{]}^{\mathrm{T}}$ is the measurement noises, which is submitted to zero mean Gaussian distribution with covariance matrix ${\mathit{Q}}_1=\mathrm{E}[{\mathit{e}}_1{\mathit{e}}_1^{\mathrm{T}}]$ .

On the contrary, although ${\mathit{s}}^{\mathit{o}}$ is easy to access through actual measurement, it may be perturbed by noise inevitably (such as sensors are fixed at moving airborne and marine platforms, or random positions). As a result, we assume the available sensor position vector is

$$\mathit{s}={\mathit{s}}^{\mathit{o}}+{\mathit{e}}_2$$

(13)

where ${\mathit{e}}_2=[e_{21},e_{22},\dots ,e_{2,\mathrm{MK}}{]}^{\mathrm{T}}$ is sensor agitation noise, which is zero mean distributed with covariance matrix ${\mathit{Q}}_2=\mathrm{E}[{\mathit{e}}_2{\mathit{e}}_2^{\mathrm{T}}]$ .

It is an easily provable fact that measurement noises ${\mathit{e}}_1$ and sensor agitation noises ${\mathit{e}}_2$ are independent of one another for originates in different.

We can gather the two kinds of noises as the vector $\mathit{e}=[{\mathit{e}}_1^{\mathrm{T}},{\mathit{e}}_2^{\mathrm{T}}{]}^{\mathrm{T}}$ and the covariance matrix can be expressed as $\mathit{Q}=\mathbf{cov}(\mathit{e})=\mathrm{E}[\mathit{e}{\mathit{e}}^{\mathrm{T}}]=\mathrm{blkdiag}[{\mathit{Q}}_1,{\mathit{Q}}_2]$ , whose block structure rooting in the statistic independence between ${\mathit{e}}_1$ and ${\mathit{e}}_2$ .

CRLB and RMSE

CRLB is regarded as the lower bound of an arbitrary unbiased estimation. This section first exhibits the CRLB with and without synchronization errors based on AOAs and TDOAs, intending to illustrate that the localization accuracy will be declined with respect to synchronization errors. Second, we derive the localization MSE based on an optimum algorithm (in the sense of achieving CRLB) that assumes the synchronization errors are absent but exist factually. The comparison of MSE assuming no synchronization errors and CRLB with synchronization errors is designed to explain the necessity to explore a novel localization algorithm considering the synchronization errors to improve the localization accuracy.

CRLB

Given that the measurement noises vector ${\mathit{e}}_1$ and sensor position noises vector ${\mathit{e}}_2$ are Gaussian distributed and independent of each other, the logarithm of the joint probability density function (PDF) of $\mathit{\xi }$ and $\mathit{s}$ parameterized on $\mathit{\rho }=[{\mathit{\mu }}^{\mathrm{T}},{\mathit{s}}^{\mathrm{T}}{]}^{\mathrm{T}}$ is

$$\begin{array}{cc}\multicolumn{1}{c}{\ln f(\mathit{\xi },\mathit{s};\mathit{\rho })}& =k_1-\frac{1}{2}{(\mathit{\xi }-{\mathit{\xi }}^o)}^{\mathrm{T}}{{\mathit{Q}}_1}^{-1}(\mathit{\xi }-{\mathit{\xi }}^o)\hfill \\ \multicolumn{1}{c}{}& +k_2-\frac{1}{2}{(\mathit{s}-{\mathit{s}}^o)}^{\mathrm{T}}{{\mathit{Q}}_2}^{-1}(\mathit{s}-{\mathit{s}}^o)\hfill \end{array}$$

(14)

where $k_1=-1/2\ln ((2\pi {)}^{2(3M-1)}\left|{\mathit{Q}}_1\right|)$ and $k_2=-1/2\ln ((2\pi {)}^{6M}\left|{\mathit{Q}}_2\right|)$ are constants. $\mathit{\mu }$ is an estimated vector consisting of source position and clock bias.

The CRLB of $\mathit{\rho }$ is equal to

$$\begin{array}{cc}\multicolumn{1}{c}{\mathrm{CRLB}\left(\mathit{\rho }\right)}& =-E{\left[{\frac{{\partial }^2\ln f(\mathit{\xi },\mathit{s};\mathit{\rho })}{\partial \mathit{\rho }\partial {\mathit{\rho }}^{\mathrm{T}}}|}_{\genfrac{}{}{0ex}{}{\mathit{\mu }={\mathit{\mu }}^o}{\mathit{\delta }={\mathit{\delta }}^o}}\right]}^{-1}\hfill \\ \multicolumn{1}{c}{}& ={\left[\begin{array}{cc}{\mathit{X}}_{11}\hfill & {\mathit{X}}_{12}\hfill \\ \multicolumn{1}{c}{{\mathit{X}}_{12}^{\mathrm{T}}}& {\mathit{X}}_{22}\hfill \end{array}\right]}^{-1}\hfill \end{array}$$

(15)

where

$$\{\begin{array}{c}{\mathit{X}}_{11}={\mathit{\Omega }}_1^{\mathrm{T}}{\mathit{Q}}_1^{-1}{\mathit{\Omega }}_1\hfill \\ \multicolumn{1}{c}{{\mathit{X}}_{12}={\mathit{\Omega }}_1^{\mathrm{T}}{\mathit{Q}}_1^{-1}{\mathit{\Omega }}_2}\\ \multicolumn{1}{c}{{\mathit{X}}_{22}={\mathit{\Omega }}_2^{\mathrm{T}}{\mathit{Q}}_1^{-1}{\mathit{\Omega }}_2+{\mathit{Q}}_2^{-1}}\end{array}$$

(16)

$$\{\begin{array}{c}{\mathit{\Omega }}_1=\frac{\partial {\mathit{\xi }}^o}{\partial {{\mathit{\mu }}^o}^{\mathrm{T}}}=[\frac{\partial {\mathit{\xi }}^o}{\partial {{\mathit{u}}^o}^{\mathrm{T}}},\frac{\partial {\mathit{\xi }}^o}{\partial {{\mathit{\delta }}^o}^{\mathrm{T}}}]={[{\mathit{\Omega }}_{11}^{\mathrm{T}},{\mathit{O}}_{(M-1)\times 2M}^{\mathrm{T}},{\mathit{I}}_{M-1}]}^{\mathrm{T}}\hfill \\ \multicolumn{1}{c}{{\mathit{\Omega }}_2=\frac{\partial {\mathit{\xi }}^o}{\partial {{\mathit{s}}^o}^{\mathrm{T}}}={[{\left(\frac{\partial {\mathit{\theta }}^o}{\partial {{\mathit{s}}^o}^{\mathrm{T}}}\right)}^{\mathrm{T}},{\left(\frac{\partial {\mathit{\beta }}^o}{\partial {{\mathit{s}}^o}^{\mathrm{T}}}\right)}^{\mathrm{T}},{\left(\frac{\partial {\mathit{r}}^o}{\partial {{\mathit{s}}^o}^{\mathrm{T}}}\right)}^{\mathrm{T}}]}^{\mathrm{T}}}\\ \multicolumn{1}{c}{{\mathit{\Omega }}_{11}={[{\left(\frac{\partial {\mathit{\theta }}^o}{\partial {{\mathit{u}}^o}^{\mathrm{T}}}\right)}^{\mathrm{T}},{\left(\frac{\partial {\mathit{\beta }}^o}{\partial {{\mathit{u}}^o}^{\mathrm{T}}}\right)}^{\mathrm{T}},{\left(\frac{\partial {\mathit{r}}^o}{\partial {{\mathit{u}}^o}^{\mathrm{T}}}\right)}^{\mathrm{T}}]}^{\mathrm{T}}}\end{array}$$

(17)

The detailed derivatives of $\partial {\mathit{\xi }}^o/\partial {{\mathit{\mu }}^o}^{\mathrm{T}}$ and $\partial {\mathit{\xi }}^o/\partial {{\mathit{s}}^o}^{\mathrm{T}}$ are addressed in Appendix 1, and the CRLB of $\mathit{\mu }$ is given as

$$\begin{array}{cc}\multicolumn{1}{c}{\mathrm{CRLB}\left(\mathit{\mu }\right)}& ={({\mathit{X}}_{11}-{\mathit{X}}_{12}^{\mathrm{T}}{\mathit{X}}_{22}^{-1}{\mathit{X}}_{12})}^{-1}\hfill \\ \multicolumn{1}{c}{}& ={\left({\mathit{\Omega }}_1^{\mathrm{T}}{({\mathit{Q}}_1+{\mathit{\Omega }}_2{\mathit{Q}}_2{\mathit{\Omega }}_2^{\mathrm{T}})}^{-1}{\mathit{\Omega }}_1\right)}^{-1}\hfill \end{array}$$

(18)

The square root of first K diagonal elements sum in equation (18) is the minimum possible source location MSE that any linear unbiased estimator can achieve. In the same way, the minimum possible MSE of synchro errors can be acquired from the later $M-1$ diagonal elements.

RMSE of source position without accounting for synchronization errors

This subsection discusses the minimum possible RMSE of source position estimator assuming the synchronization errors are absent but existing actually. The derivation is based on Taylor series linearization expansion, which can achieve the CRLB after convergence to the global minimum when sensors are synchronized.

Define three quantities parameterized on source position $\mathit{u}$ as follows

$$\{\begin{array}{c}{\stackrel{⌣}{\theta }}_i(\mathit{u},{\mathit{s}}_i)=\arctan \frac{y-y_i}{x-x_i}\hfill \\ \multicolumn{1}{c}{{\stackrel{⌣}{\beta }}_i(\mathit{u},{\mathit{s}}_i)=\arctan \frac{z-z_i}{\sqrt{{(x-x_i)}^2+{(y-y_i)}^2}}}\end{array},i=1,2,\dots ,M$$

(19)

$${\stackrel{⌣}{r}}_{i1}(\mathit{u},{\mathit{s}}_i)=\Vert \mathit{u}-{\mathit{s}}_i\Vert -\Vert \mathit{u}-{\mathit{s}}_1\Vert ,i=2,\dots ,M$$

(20)

Define the vector

$$\{\begin{array}{c}\stackrel{⌣}{\mathit{\theta }}(\mathit{u},\mathit{s})={[{\stackrel{⌣}{\theta }}_1(\mathit{u},{\mathit{s}}_i),{\stackrel{⌣}{\theta }}_2(\mathit{u},{\mathit{s}}_i),\dots ,{\stackrel{⌣}{\theta }}_M(\mathit{u},{\mathit{s}}_i)]}^{\mathrm{T}}\hfill \\ \multicolumn{1}{c}{\stackrel{⌣}{\mathit{\beta }}(\mathit{u},\mathit{s})={[{\stackrel{⌣}{\beta }}_1^o(\mathit{u},{\mathit{s}}_i),{\stackrel{⌣}{\beta }}_2^o(\mathit{u},{\mathit{s}}_i),\dots ,{\stackrel{⌣}{\beta }}_M(\mathit{u},{\mathit{s}}_i)]}^{\mathrm{T}}}\\ \multicolumn{1}{c}{\stackrel{⌣}{\mathit{r}}(\mathit{u},\mathit{s})={[{\stackrel{⌣}{r}}_{21}(\mathit{u},{\mathit{s}}_i),{\stackrel{⌣}{r}}_{31}(\mathit{u},{\mathit{s}}_i),\dots ,{\stackrel{⌣}{r}}_{M1}(\mathit{u},{\mathit{s}}_i)]}^{\mathrm{T}}}\\ \multicolumn{1}{c}{{\mathit{\delta }}^{\mathit{o}}={[{\delta }_{21}^o,{\delta }_{31}^o,\dots ,{\delta }_{M1}^o]}^{\mathrm{T}}}\\ \multicolumn{1}{c}{\mathit{s}={[{\mathit{s}}_1^{\mathrm{T}},{\mathit{s}}_2^{\mathrm{T}},\dots ,{\mathit{s}}_M^{\mathrm{T}}]}^{\mathrm{T}}}\end{array}$$

(21)

where $\stackrel{⌣}{\mathit{r}}(\mathit{u},\mathit{s})$ is supposed to be only contaminated by the measurement noises and sensor position noises in spite of the actual occurrence of synchronization errors.

We define $\stackrel{⌣}{\mathit{\xi }}(\mathit{u},\mathit{s})=[{\stackrel{⌣}{\mathit{\theta }}}^{\mathrm{T}}(\mathit{u},\mathit{s}),{\stackrel{⌣}{\mathit{\beta }}}^{\mathrm{T}}(\mathit{u},\mathit{s}),{\stackrel{⌣}{\mathit{r}}}^{\mathrm{T}}(\mathit{u},\mathit{s}){]}^{\mathrm{T}}$ and apply Taylor series expansion up to the first-order term

$$\stackrel{⌣}{\mathit{\xi }}(\mathit{u},\mathit{s})\stackrel{·}{=}\stackrel{⌣}{\mathit{\xi }}({\mathit{u}}^o,{\mathit{s}}^o)+{\frac{\partial \stackrel{⌣}{\mathit{\xi }}(\mathit{u},\mathit{s})}{\partial {\mathit{u}}^{\mathrm{T}}}|}_{{\mathit{u}}^o,{\mathit{s}}^o}(\mathit{u}-{\mathit{u}}^o)+{\frac{\partial \stackrel{⌣}{\mathit{\xi }}(\mathit{u},\mathit{s})}{\partial {\mathit{s}}^{\mathrm{T}}}|}_{{\mathit{u}}^o,{\mathit{s}}^o}(\mathit{s}-{\mathit{s}}^o)$$

(22)

Since ${\frac{\partial \stackrel{⌣}{\mathit{\xi }}(\mathit{u},\mathit{s})}{\partial {\mathit{u}}^{\mathrm{T}}}|}_{{\mathit{u}}^o,{\mathit{s}}^o}$ and ${\frac{\partial \stackrel{⌣}{\mathit{\xi }}(\mathit{u},\mathit{s})}{\partial {\mathit{s}}^{\mathrm{T}}}|}_{{\mathit{u}}^o,{\mathit{s}}^o}$ have the same form as $\partial {\mathit{\xi }}^o/\partial {{\mathit{u}}^o}^{\mathrm{T}}$ and $\partial {\mathit{\xi }}^o/\partial {{\mathit{s}}^o}^{\mathrm{T}}$ given in Appendix 1, equation (22) can be rewritten as

$$\stackrel{⌣}{\mathit{\xi }}(\mathit{u},\mathit{s})=\stackrel{⌣}{\mathit{\xi }}({\mathit{u}}^o,{\mathit{s}}^o)+{\mathit{\Omega }}_{11}(\mathit{u}-{\mathit{u}}^o)+{\mathit{\Omega }}_2(\mathit{s}-{\mathit{s}}^o)$$

(23)

By subtracting the observation vector $\mathit{\xi }$ to $\stackrel{⌣}{\mathit{\xi }}(\mathit{u},\mathit{s})$ in equation (23) and replacing $(\mathit{s}-{\mathit{s}}^o)$ by ${\mathit{e}}_2$ , the measurement noises term is

$${\mathit{\epsilon }}_1=\mathit{\xi }-\stackrel{⌣}{\mathit{\xi }}(\mathit{u},\mathit{s})=\mathit{\xi }-\stackrel{⌣}{\mathit{\xi }}({\mathit{u}}^o,{\mathit{s}}^o)-{\mathit{\Omega }}_{11}(\mathit{u}-{\mathit{u}}^o)-{\mathit{\Omega }}_2{\mathit{e}}_2$$

(24)

We gather the measurement noises and sensor position errors as $\mathit{\epsilon }=[\begin{array}{cc}{\mathit{\epsilon }}_1^{\mathrm{T}}\hfill & {\mathit{e}}_2^{\mathrm{T}}\hfill \end{array}{]}^{\mathrm{T}}$ and equation (24) can be formulated as

$$\mathit{\Xi }\mathit{\epsilon }=\mathit{\xi }-\stackrel{⌣}{\mathit{\xi }}({\mathit{u}}^o,{\mathit{s}}^o)-{\mathit{\Omega }}_{11}(\mathit{u}-{\mathit{u}}^o)$$

(25)

where $\mathit{\Xi }=[\begin{array}{cc}{\mathit{I}}_{3M-1}\hfill & {\mathit{\Omega }}_2\hfill \end{array}]$ is full row rank.

We find $\mathit{u}$ by minimizing

$${\mathit{\epsilon }}^{\mathrm{T}}{\mathit{Q}}^{-1}\mathit{\epsilon }$$

(26)

with the constraints in equation (25), which can be reformulated as

$$\begin{array}{c}\underset{\mathit{u}}{\min }{(\mathit{\xi }-\stackrel{⌣}{\mathit{\xi }}({\mathit{u}}^o,{\mathit{s}}^o)-{\mathit{\Omega }}_{11}(\mathit{u}-{\mathit{u}}_o))}^{\mathrm{T}}{\left(\mathit{\Xi }\mathit{Q}{\mathit{\Xi }}^{\mathrm{T}}\right)}^{-1}\hfill \\ (\mathit{\xi }-\stackrel{⌣}{\mathit{\xi }}({\mathit{u}}^o,{\mathit{s}}^o)-{\mathit{\Omega }}_{11}(\mathit{u}-{\mathit{u}}_o))\hfill \end{array}$$

(27)

Taking derivation in equation (27) with respect to u and making the gradient zero yields

$$\begin{array}{c}\mathit{u}={\mathit{u}}_o+{\left({\mathit{\Omega }}_{11}^{\mathrm{T}}{\left(\mathit{\Xi }\mathit{Q}{\mathit{\Xi }}^{\mathrm{T}}\right)}^{-1}{\mathit{\Omega }}_{11}\right)}^{-1}\hfill \\ {\mathit{\Omega }}_{11}^{\mathrm{T}}{\left(\mathit{\Xi }\mathit{Q}{\mathit{\Xi }}^{\mathrm{T}}\right)}^{-1}(\mathit{\xi }-\stackrel{⌣}{\mathit{\xi }}({\mathit{u}}^o,{\mathit{s}}^o))\hfill \end{array}$$

(28)

Based on equation (24), let $\mathit{u}={\mathit{u}}^o$ and take the actual synchronization errors into consideration, we have the following expression

$$\mathit{\xi }-\stackrel{⌣}{\mathit{\xi }}({\mathit{u}}^o,{\mathit{s}}^o)={\mathit{e}}_1+{\mathit{\Omega }}_2{\mathit{e}}_2+\mathrm{\Delta }\mathit{\xi }$$

(29)

where $\mathrm{\Delta }\mathit{\xi }$ contains the synchronization errors.

Moreover, the estimate error of source position is given by

$$\begin{array}{c}\mathit{u}-{\mathit{u}}_o={\left({\mathit{\Omega }}_{11}^{\mathrm{T}}{\left(\mathit{\Xi }\mathit{Q}{\mathit{\Xi }}^{\mathrm{T}}\right)}^{-1}{\mathit{\Omega }}_{11}\right)}^{-1}\hfill \\ {\mathit{\Omega }}_{11}^{\mathrm{T}}{\left(\mathit{\Xi }\mathit{Q}{\mathit{\Xi }}^{\mathrm{T}}\right)}^{-1}({\mathit{e}}_1+{\mathit{\Omega }}_2{\mathit{e}}_2+\mathrm{\Delta }\mathit{\xi })\hfill \end{array}$$

(30)

Taking the expectation of equation (30) indicates that the estimator $\mathit{u}$ is biased for the presence of $\mathrm{\Delta }\xi$ , whereas ${\mathit{e}}_1$ and ${\mathit{e}}_2$ are zero mean. Multiplying equation (30) with its transpose and taking expectation yields the covariance matrix

$$\begin{array}{cc}\multicolumn{1}{c}{\mathbf{cov}(\mathit{u}-{\mathit{u}}^o)}& =E\left[(\mathit{u}-{\mathit{u}}^o){(\mathit{u}-{\mathit{u}}^o)}^{\mathrm{T}}\right]\hfill \\ \multicolumn{1}{c}{}& ={\left({\mathit{\Omega }}_{11}^{\mathrm{T}}{\left(\mathit{\Xi }\mathit{Q}{\mathit{\Xi }}^{\mathrm{T}}\right)}^{-1}{\mathit{\Omega }}_{11}\right)}^{-1}\hfill \\ \multicolumn{1}{c}{}& +{\left({\mathit{\Omega }}_{11}^{\mathrm{T}}{\left(\mathit{\Xi }\mathit{Q}{\mathit{\Xi }}^{\mathrm{T}}\right)}^{-1}{\mathit{\Omega }}_{11}\right)}^{-1}{\mathit{\Omega }}_{11}^{\mathrm{T}}{\left(\mathit{\Xi }\mathit{Q}{\mathit{\Xi }}^{\mathrm{T}}\right)}^{-1}\hfill \\ \multicolumn{1}{c}{}& \mathrm{\Delta }\mathit{\xi }\mathrm{\Delta }{\mathit{\xi }}^{\mathrm{T}}{\left(\mathit{\Xi }\mathit{Q}{\mathit{\Xi }}^{\mathrm{T}}\right)}^{-1}{\mathit{\Omega }}_{11}{\left({\mathit{\Omega }}_{11}^{\mathrm{T}}{\left(\mathit{\Xi }\mathit{Q}{\mathit{\Xi }}^{\mathrm{T}}\right)}^{-1}{\mathit{\Omega }}_{11}\right)}^{-1}\hfill \end{array}$$

(31)

Apparently, equation (31) is composed of two terms. The first term is the RMSE with inaccurate sensor positions but accurate clocks. The second exhibits the additional error resulting from the sensor clock non-synchronization, the trace of which is the increase in RMSE.

In view that the synchronization errors have not been taken into consideration, that is to say it is deemed to be zero in the estimator. Consequently, the RMSE of synchronization error can be obtained from

$$\mathbf{cov}(\mathit{\delta }-{\mathit{\delta }}^{\mathit{o}})={\mathit{\delta }}^{\mathit{o}}{\mathit{\delta }}^{\mathit{o}\mathrm{T}}$$

(32)

Comparison of CRLB and RMSE

In order to obtain some insight, we compare the CRLB with and without synchronization errors with the RMSE not accounting for synchronization errors. A sensor network consisting of three sensors whose locations are listed in Table 1 is employed to determine the source at ${\mathit{u}}^o=[550,500,400{]}^{\mathrm{T}}$ . The AOA and TDOA measurements are contaminated by zero-mean Gaussian noises with the covariance matrices ${\mathit{Q}}_{\mathrm{AOA}}={\sigma }_{\mathrm{AOA}}^2\mathit{I}$ and ${\mathit{Q}}_{\mathrm{TDOA}}={\sigma }_{\mathrm{TDOA}}^2\mathit{R}$ , where $\mathit{I}$ is an identity matrix and $\mathit{R}$ is equal to unity in the diagonal elements and 0.5 otherwise. Measurement noise power of AOA and TDOA are ${\sigma }_{\mathrm{AOA}}^2={10}^{-3.5}(\mathrm{ra}{\mathrm{d}}^2)$ and ${\sigma }_{\mathrm{TDOA}}^2={10}^{0.5}({\mathrm{m}}^2)$ , respectively, which are modified to represent the different noise levels. Therefore, the covariance matrix of measurement noises can be expressed as ${\mathit{Q}}_1=\mathrm{blkdiag}[{\mathit{Q}}_{\mathrm{AOA}},{\mathit{Q}}_{\mathrm{TDOA}}]$ . In addition, the covariance matrix of sensor position perturbation noises is ${\mathit{Q}}_{\mathrm{s}}={\sigma }_{\mathrm{s}}^2\mathit{R}$ , where ${\sigma }_{\mathrm{s}}^2=100({\mathrm{m}}^2)$ is the noise power. We set clock bias ${\mathit{\delta }}^{\mathit{o}}={\sigma }_{\mathrm{c}}[0.8,2.6{]}^{\mathrm{T}}(\mathrm{ns})$ and ${\sigma }_{\mathrm{c}}$ is the standard deviation of synchronization errors and varies from 10 to 30 dBns.

Table 1.

Positions of sensors (unit: m).

Sensor no.	1	2	3
$x_i^o$	550	−200	−700
$y_i^o$	−350	750	−900
$z_i^o$	−300	250	−550

Numerical experiment results about position accuracy are shown in Figure 2, where Figure 2(a) shows a 6 m difference in the CRLB with and without synchronization errors for that the decline of synchro errors acting on CRLB is a fixed coefficient resulting from operation of derivation. The RMSE curve assuming no synchro errors can achieve the CRLB without synchro when ${\sigma }_{\mathrm{c}}$ is very small and gradually exceeds the CRLB with synchro errors when ${\sigma }_{\mathrm{c}}$ increases to 21 dBns approximately. Afterward, it deviates from the two CRLB farther and farther even up to 100 m errors in position estimation. Figure 2(b) shows the difference in RMSE assuming no synchro errors and CRLB with synchro errors to prove the performance improvement when synchro errors are considered. One can see that when ${\sigma }_{\mathrm{c}}<21\mathrm{dBns}$ the RMSE is lower than CRLB, which indicates that the localization model by ignoring the small synchro errors can achieve a little higher accuracy. However, the introduction of synchro errors will bring a more and more improvement as the ${\sigma }_{\mathrm{c}}$ increases. For example, about 35 m promotion appears in position precision just when ${\sigma }_{\mathrm{c}}$ is up to 26 dBns. Hence, it is worth introducing the synchro errors in the localization model for it brings little performance reduction in small synchro errors and large promotion in big errors.

Figure 2.

Comparison of CRLBs with and without synchronization errors and the RMSE assuming no synchro errors in the view of position accuracy. The dotted line in (b) represents the negative difference of RMSE and CRLB with synchro error in (a) and solid line represents the active difference.

As the clock bias can be estimated by introducing synchro errors, we compare the RMSE assuming no synchro errors and CRLB with synchro errors to illustrate the estimate accuracy promotion of clock bias, which is exhibited in Figure 3 based on equations (18) and (32). The CRLB without synchro errors is not plotting for its inability to estimate clock bias. A similar appearance in Figure 3(a) is worth investigating that the RMSE is lower than CRLB before ${\sigma }_{\mathrm{c}}=21\mathrm{dBns}$ , which suggests that the clock bias would prefer to be ignored in small synchro errors. It can be interpreted that the estimation accuracy of clock bias is trade by the degree of freedom (DOF) of estimate parameters with limited measurements. From Figure 3(b), we can see that there are at most 100 ns accuracy decrease when ${\sigma }_{\mathrm{c}}<21\mathrm{dBns}$ , and 700 ns improvement when ${\sigma }_{\mathrm{c}}=30\mathrm{dBns}$ , even higher promotion when sensors fail to synchronize clocks in actual project. From the aspect of the clock bias estimate, synchro errors should be taken into consideration for performance improvement.

Figure 3.

Comparison of CRLB with synchro errors and the RMSE assuming no synchro errors in the view of synchro clock bias. The dotted line in (b) represents the negative difference of RMSE and CRLB with synchro error in (a) and solid line represents the active difference.

CTLS localization method

It is well known that the TLS localization method is able to effectively restrain the measurement noises with the noisy coefficient matrix. Nevertheless, the TLS method can achieve the optimization performance when all disturbance components are subject to independent and identical distribution in the coefficient matrix and the observation vector. Considering that the disturbance in $\mathit{A}(\mathit{\xi })$ and $\mathit{b}(\mathit{\xi },\mathit{s})$ must be structuralized for $\mathit{A}({\mathit{\xi }}^{\mathit{o}})$ and $\mathit{b}({\mathit{\xi }}^{\mathit{o}},s^{\mathit{o}})$ are the function matrix with respect to ${\mathit{\xi }}^{\mathit{o}}$ collectively, the CTLS method can enhance the location accuracy relative to the TLS method, which cannot attain the optimal statistical performance.

Construction of the CTLS localization model

As aforementioned, only noisy measurements and imprecise sensor positions can be utilized to determine the source position. We first expand $\mathit{A}({\mathit{\xi }}^{\mathit{o}})$ and $\mathit{b}({\mathit{\xi }}^{\mathit{o}},s^{\mathit{o}})$ with respect to $\mathit{\xi }$ and $\mathit{s}$ through Taylor series and ignore the second and higher order terms

$$\{\begin{array}{c}\mathit{A}\left({\mathit{\xi }}^{\mathit{o}}\right)\stackrel{.}{=}\mathit{A}\left(\mathit{\xi }\right)+\sum _{k=1}^{3M-1}{\stackrel{·}{\mathit{A}}}_{1k}\left(\mathit{\xi }\right)·{〈{\mathit{\xi }}^{\mathit{o}}-\mathit{\xi }〉}_k\hfill \\ \multicolumn{1}{c}{=\mathit{A}\left(\mathit{\xi }\right)-\sum _{k=1}^{3M-1}{\stackrel{·}{\mathit{A}}}_{1k}\left(\mathit{\xi }\right)·{〈{\mathit{e}}_1〉}_k}\\ \multicolumn{1}{c}{\mathit{b}({\mathit{\xi }}^{\mathit{o}},{\mathit{s}}^{\mathit{o}})\stackrel{.}{=}\mathit{b}(\mathit{\xi },\mathit{s})+{\mathit{B}}_1(\mathit{\xi },\mathit{s})({\mathit{\xi }}^{\mathit{o}}-\mathit{\xi })+{\mathit{B}}_2(\mathit{\xi },\mathit{s})({\mathit{s}}^{\mathit{o}}-\mathit{s})}\\ \multicolumn{1}{c}{=\mathit{b}(\mathit{\xi },\mathit{s})-{\mathit{B}}_1(\mathit{\xi },\mathit{s}){\mathit{e}}_1-{\mathit{B}}_2\left(\mathit{\xi }\right){\mathit{e}}_2}\end{array}$$

(33)

where

$$\begin{array}{c}{\stackrel{·}{\mathit{A}}}_{1k}\left(\mathit{\xi }\right)={\frac{\partial \mathit{A}\left({\mathit{\xi }}^{\mathit{o}}\right)}{\partial {〈{\mathit{\xi }}^{\mathit{o}}〉}_k}|}_{{\mathit{\xi }}^{\mathit{o}}=\mathit{\xi }},\hfill \\ {\mathit{B}}_1(\mathit{\xi },\mathit{s})={\frac{\partial \mathit{b}({\mathit{\xi }}^{\mathit{o}},{\mathit{s}}^{\mathit{o}})}{\partial {\mathit{\xi }}^{\mathit{o}\mathrm{T}}}|}_{{\mathit{\xi }}^{\mathit{o}}=\mathit{\xi },{\mathit{s}}^{\mathit{o}}=\mathit{s}},{\mathit{B}}_2\left(\mathit{\xi }\right)={\frac{\partial \mathit{b}({\mathit{\xi }}^{\mathit{o}},{\mathit{s}}^{\mathit{o}})}{\partial {\mathit{s}}^{\mathit{o}\mathrm{T}}}|}_{{\mathit{\xi }}^{\mathit{o}}=\mathit{\xi }}\hfill \end{array}$$

(34)

Noting that estimated vector $\mathit{\mu }=[{\mathit{u}}^{\mathrm{T}},{\mathit{\delta }}^{\mathrm{T}}{]}^{\mathrm{T}}$ , based on equations (9) and (33), the CLTS method can be modeled as

$$\{\begin{array}{c}min{\mathit{e}}^{\mathrm{T}}{\mathit{Q}}^{-1}\mathit{e}\hfill \\ \multicolumn{1}{c}{\mathrm{s}.\mathrm{t}.(\mathit{A}\left(\mathit{\xi }\right)-\sum _{k=1}^{3M-1}{\stackrel{·}{\mathit{A}}}_{1k}\left(\mathit{\xi }\right){〈{\mathit{e}}_1〉}_k)\mathit{\mu }=\mathit{b}(\mathit{\xi },\mathit{s})-{\mathit{B}}_1(\mathit{\xi },\mathit{s}){\mathit{e}}_1-{\mathit{B}}_2\left(\mathit{\xi }\right){\mathit{e}}_2}\end{array}$$

(35)

Furthermore, equation (35) can be simplified as

$$\{\begin{array}{c}\min {\mathit{e}}^{\mathrm{T}}{\mathit{Q}}^{-1}\mathit{e}\hfill \\ \multicolumn{1}{c}{\mathrm{s}.\mathrm{t}.\mathit{b}(\mathit{\xi },\mathit{s})-\mathit{A}\left(\mathit{\xi }\right)\mathit{\mu }=({\mathit{B}}_1(\mathit{\xi },\mathit{s})-{\mathit{E}}_1(\mathit{\xi },\mathit{\mu })){\mathit{e}}_1+{\mathit{B}}_2\left(\mathit{\xi }\right){\mathit{e}}_2=\mathit{\Phi }(\mathit{\xi },\mathit{s},\mathit{\mu })\mathit{e}}\end{array}$$

(36)

where

$$\{\begin{array}{c}{\mathit{E}}_1(\mathit{\xi },\mathit{\mu })=[{\stackrel{·}{\mathit{A}}}_{11}\left(\mathit{\xi }\right)\mathit{\mu },{\stackrel{·}{\mathit{A}}}_{12}\left(\mathit{\xi }\right)\mathit{\mu },\dots ,{\stackrel{·}{\mathit{A}}}_{1,3M-1}\left(\mathit{\xi }\right)\mathit{\mu }]\hfill \\ \multicolumn{1}{c}{\mathit{\Phi }(\mathit{\xi },\mathit{s},\mathit{\mu })=[{\mathit{\Phi }}_1(\mathit{\xi },\mathit{s},\mathit{\mu }),{\mathit{\Phi }}_2\left(\mathit{\xi }\right)]=[{\mathit{B}}_1(\mathit{\xi },\mathit{s})-{\mathit{E}}_1(\mathit{\xi },\mathit{\mu }),{\mathit{B}}_2\left(\mathit{\xi }\right)]}\end{array}$$

(37)

The details of algebraic translation in equation (37) are cited in Appendix 2.

We dispose equation (36) as an optimization problem, which is described in Appendix 3 for the aim to achieve global optimal solutions just for $\mathit{\mu }$ . The objective function is given by

$$\begin{array}{c}\underset{{\mathit{u}}^o,{\mathit{\delta }}^o}{\min }{(\mathit{A}\left(\mathit{\xi }\right)\mathit{\mu }-\mathit{b}(\mathit{\xi },\mathit{s}))}^{\mathrm{T}}\hfill \\ {\left(\mathit{\Phi }(\mathit{\xi },\mathit{s},\mathit{\mu })\mathit{Q}{\mathit{\Phi }}^{\mathrm{T}}(\mathit{\xi },\mathit{s},\mathit{\mu })\right)}^{-1}(\mathit{A}\left(\mathit{\xi }\right)\mathit{\mu }-\mathit{b}(\mathit{\xi },\mathit{s}))\hfill \end{array}$$

(38)

Newton iteration solution to CTLS method

As the closed solution of equation (38) is hardly derived to solve the complicated problem in equation (38), we introduce Newton iteration solution equipped with repaid convergence rate in this subsection, which has been regarded as a classic tool to solve the CTLS problem.

Note the cost function in equation (38) as

$$J\left(\mathit{\mu }\right)={\mathit{\phi }}^{\mathrm{T}}(\mathit{\xi },\mathit{s},\mathit{\mu })\mathit{\Psi }(\mathit{\xi },\mathit{s},\mathit{\mu })\mathit{\phi }(\mathit{\xi },\mathit{s},\mathit{\mu })$$

(39)

where

$$\{\begin{array}{c}\mathit{\phi }(\mathit{\xi },\mathit{s},\mathit{\mu })=\mathit{A}\left(\mathit{\xi }\right)\mathit{\mu }-\mathit{b}(\mathit{\xi },\mathit{s})\hfill \\ \multicolumn{1}{c}{\mathit{\Psi }(\mathit{\xi },\mathit{s},\mathit{\mu })={\left(\mathit{\Phi }(\mathit{\xi },\mathit{s},\mathit{\mu })\mathit{Q}{\mathit{\Phi }}^{\mathrm{T}}(\mathit{\xi },\mathit{s},\mathit{\mu })\right)}^{-1}}\end{array}$$

(40)

The derivation of $J(\mathit{\mu })$ with respect to $\mathit{\mu }$ yields the gradient vector

$$\begin{array}{c}\mathit{h}\left(\mathit{\mu }\right)=2{\mathit{A}}^{\mathrm{T}}\left(\mathit{\xi }\right)\mathit{\Psi }(\mathit{\xi },\mathit{s},\mathit{\mu })\mathit{\phi }(\mathit{\xi },\mathit{s},\mathit{\mu })\hfill \\ +{\left(\frac{\partial \mathrm{vec}\left(\mathit{\Psi }(\mathit{\xi },\mathit{s},\mathit{\mu })\right)}{\partial {\mathit{\mu }}^{\mathrm{T}}}\right)}^{\mathrm{T}}(\mathit{\phi }(\mathit{\xi },\mathit{s},\mathit{\mu })\otimes \mathit{\phi }(\mathit{\xi },\mathit{s},\mathit{\mu }))\hfill \end{array}$$

(41)

where

$$\begin{array}{c}\frac{\partial \mathrm{vec}\left(\mathit{\Psi }(\mathit{\xi },\mathit{s},\mathit{\mu })\right)}{\partial {\mathit{\mu }}^{\mathrm{T}}}\hfill \\ =(\left(\mathit{\Psi }(\mathit{\xi },\mathit{s},\mathit{\mu })\mathit{\Phi }(\mathit{\xi },\mathit{s},\mathit{\mu })\mathit{Q}\right)\otimes \mathit{\Psi }(\mathit{\xi },\mathit{s},\mathit{\mu }))\stackrel{·}{\overline{\mathit{A}}}\left(\mathit{\xi }\right)\hfill \\ (\mathit{\Psi }(\mathit{\xi },\mathit{s},\mathit{\mu })\otimes \left(\mathit{\Psi }(\mathit{\xi },\mathit{s},\mathit{\mu })\mathit{\Phi }(\mathit{\xi },\mathit{s},\mathit{\mu })\mathit{Q}\right))\mathit{\Pi }\stackrel{·}{\overline{\mathit{A}}}\left(\mathit{\xi }\right)\hfill \end{array}$$

(42)

Matrix $\mathit{\Pi }$ satisfies $\partial \mathrm{vec}({\mathit{\Psi }}^{\mathrm{T}}(\mathit{\xi },\mathit{s},\mathit{\mu }))/\partial {\mathit{\mu }}^{\mathrm{T}}=\mathit{\Pi }(\partial \mathrm{vec}(\mathit{\Psi }(\mathit{\xi },\mathit{s},\mathit{\mu }))/\partial {\mathit{\mu }}^{\mathrm{T}})$ and $\stackrel{·}{\overline{\mathit{A}}}(\mathit{\xi })$ can be expressed as

$$\begin{array}{c}\stackrel{·}{\overline{\mathit{A}}}\left(\mathit{\xi }\right)=\hfill \\ {[{\stackrel{·}{\mathit{A}}}_{11}^{\mathrm{T}}\left(\mathit{\xi }\right),{\stackrel{·}{\mathit{A}}}_{12}^{\mathrm{T}}\left(\mathit{\xi }\right),\dots ,{\stackrel{·}{\mathit{A}}}_{1,3M-1}^{\mathrm{T}}\left(\mathit{\xi }\right),{\mathit{O}}_{(K+M-1)\times \left(\left(3M\right)(3M-1)\right)}]}^{\mathrm{T}}\hfill \end{array}$$

(43)

Beyond that, the Hessian matrix of $J(\mathit{\mu })$ is given by

$$\begin{array}{cc}\multicolumn{1}{c}{\mathit{H}\left(\mathit{\mu }\right)}& =\frac{{\partial }^{2}J\left(\mathit{\mu }\right)}{\partial \mathit{\mu }\partial {\mathit{\mu }}^{\mathrm{T}}}=\frac{\partial \mathit{h}\left(\mathit{\mu }\right)}{\partial {\mathit{\mu }}^{\mathrm{T}}}\hfill \\ \multicolumn{1}{c}{}& =2{\mathit{A}}^{\mathrm{T}}\left(\mathit{\xi }\right)\mathit{\Psi }(\mathit{\xi },\mathit{s},\mathit{\mu })\mathit{A}\left(\mathit{\xi }\right)\hfill \\ \multicolumn{1}{c}{}& +2{(\mathit{\phi }(\mathit{\xi },\mathit{s},\mathit{\mu })\otimes \mathit{A}\left(\mathit{\xi }\right))}^{\mathrm{T}}\left(\frac{\partial \mathrm{vec}\left(\mathit{\Psi }(\mathit{\xi },\mathit{s},\mathit{\mu })\right)}{\partial {\mathit{\mu }}^{\mathrm{T}}}\right)\hfill \\ \multicolumn{1}{c}{}& +{\left(\frac{\partial \mathrm{vec}\left(\mathit{\Psi }(\mathit{\xi },\mathit{s},\mathit{\mu })\right)}{\partial {\mathit{\mu }}^{\mathrm{T}}}\right)}^{\mathrm{T}}(\mathit{\phi }(\mathit{\xi },\mathit{s},\mathit{\mu })\otimes \mathit{A}\left(\mathit{\xi }\right))\hfill \\ \multicolumn{1}{c}{}& +{\left(\frac{\partial \mathrm{vec}\left(\mathit{\Psi }(\mathit{\xi },\mathit{s},\mathit{\mu })\right)}{\partial \mathit{\mu }}\right)}^{\mathrm{T}}({\mathit{I}}_{3M-1}\otimes \mathit{\phi }(\mathit{\xi },\mathit{s},\mathit{\mu })\mathit{A}\left(\mathit{\xi }\right))\hfill \end{array}$$

(44)

Based on the gradient vector and Hessian matrix in equations (41) and (44), the specific realization of the Newtown iteration method is listed as

• Step 1: compute a rough estimate ${\mathit{\mu }}^{(0)}$ as the initial value, which is needed for iteration, then note $k:=0$ and set the threshold $\epsilon ={10}^{-8}$ .

• Step 2: make a judgment if $\Vert \mathit{h}({\mathit{\mu }}^{(k)}){\Vert }_2\le \epsilon$ and pause the iteration with YES, otherwise turn to step 3.

• Step 3: update the estimate by ${\mathit{\mu }}^{(k+1)}={\mathit{\mu }}^{(k)}-{\lambda }^k{\mathit{H}}^{-1}({\mathit{\mu }}^{(k)})\mathit{h}({\mathit{\mu }}^{(k)})$ , then let $k:=k+1$ and transform to step 2.

In order to ensure the convergence of Newton iteration solution, we select the successive step-size reduction Armijo rule. The main idea is that choose an initial step size ${\lambda }_0$ and if the corresponding vector ${\mathit{\mu }}^{(k+1)}$ does not yield an improved value of $J$ , the step size is reduced, perhaps repeatedly, by a certain factor, until the value of $J$ is improved. Then, we present this rule in detail. Here, fixed scalars $\beta$ and $\rho$ with $0<\beta <1$ and $0<\rho <1$ . Then, we set ${\lambda }^k={\beta }^{m_k}{\lambda }_0$ , where $m_k$ is the first non-negative integer $m$ for which

$$\begin{array}{c}J\left({\mathit{\mu }}^{\left(k\right)}\right)-J\left({\mathit{\mu }}^{(k+1)}\right)\hfill \\ >-\rho {\beta }^{m_k}{\lambda }_0\mathit{h}\left({\mathit{\mu }}^{\left(k\right)}\right){\mathit{H}}^{-1}\left({\mathit{\mu }}^{\left(k\right)}\right)\mathit{h}\left({\mathit{\mu }}^{\left(k\right)}\right)\hfill \end{array}$$

(45)

Usually, $\rho$ is chosen close to zero, for example, $\rho \in [{10}^{-5},{10}^{-1}]$ . The reduction factor $\beta$ is chosen from 0.5 to 0.1 depending on the confidence we have on the quality of the initial step size ${\lambda }_0$ . We always take ${\lambda }_0=1$ and multiply the direction ${\mathit{H}}^{-1}({\mathit{\mu }}^{(k)})\mathit{h}({\mathit{\mu }}^{(k)})$ by a scaling factor. In this article, we set $\rho ={10}^{-3}$ and $\beta =0.3$ .

It is generally known that the convergence of Newton iteration algorithm is related to the initial estimate. We can easily obtain the initial estimate by the TLS method or the close-form weighted least squares (WLS) method resorting to the pseudo-linear equation (equation (9)). Hence, the convergence can be ensured.

In the following, the computational complexity of the proposed method is analyzed. Assuming that the initial estimate has been known, we only need to discuss the calculation in Newton iteration solution. The calculation complexity will be evaluated by the times of multiplication. Before analysis, we have to restate that M denotes the sensor amount, and the calculation of all unknowns is presented in Table 2.

Table 2.

Calculation of all unknowns.

Unknowns	Calculation
$\mathit{\phi }(\mathit{\xi },\mathit{s},\mathit{\mu })$	$3M^2+5M-2$
$\mathit{\Psi }(\mathit{\xi },\mathit{s},\mathit{\mu })$	$216M^3-144M^2+45M-2$
$\partial \mathrm{vec}(\mathit{\Psi }(\mathit{\xi },\mathit{s},\mathit{\mu }))/\partial {\mathit{\mu }}^{\mathrm{T}}$	$162M^5+135M^4-105M^3+30M^2-2M$
$\mathit{h}(\mathit{\mu })$	$18M^3+44M^2-34M+6$
$\mathit{H}(\mathit{\mu })$	$36M^4+138M^3+76M^2-102M+20$
${\mathit{H}}^{-1}(\mathit{\mu })$	$O((M+2{)}^3)$
${\mathit{H}}^{-1}(\mathit{\mu })\mathit{h}(\mathit{\mu })$	$(M+2{)}^2$

Therefore, the calculation of one iteration is sum of all unknowns in the list, that is $162M^5+171M^4+267M^3+10M^2-84M$ $+26+O((M+2{)}^3)$ . If the iteration times $N_{\mathrm{iter}}$ are defined, the total calculation is $(162M^5+171M^4+267M^3+10M^2-84M+26$ $+O((M+2{)}^3))N_{\mathrm{iter}}$ , which is much larger than the close-form method like WLS and TLS methods and slightly higher than the Taylor series method. However, the proposed method derived from a pseudo-linear equation can achieve high position accuracy and equips with a higher. The value of $N_{\mathrm{iter}}$ and the convergence of Newton iteration are examined in section “Convergence speed of the proposed method.”

Theoretical analysis of CTLS method

In this subsection, we would demonstrate that the optimal estimate performance of the CTLS method with the aid of small error analysis. The estimator is supposed to yield a solution ${\mathit{\mu }}_{\mathrm{CTLS}}$ after convergence to the global minimum, that is ${\mathit{\mu }}_{\mathrm{CTLS}}=\underset{k\to +\infty }{lim}{\mathit{\mu }}^{(k)}$ . Substituting ${\mathit{\mu }}_{\mathrm{CTLS}}$ into gradient vector in equation (41) results in

$$\underset{k\to +\infty }{lim}\mathit{h}\left({\mathit{\mu }}^{\left(k\right)}\right)=\mathit{h}\left({\mathit{\mu }}_{\mathrm{CTLS}}\right)={\frac{\partial J\left(\mathit{\mu }\right)}{\partial \mathit{\mu }}|}_{\mathit{\mu }={\mathit{\mu }}_{\mathrm{CTLS}}}={\mathit{0}}_K$$

(46)

Using equation (41), the gradient vector is equal to

$$\begin{array}{c}2{\mathit{A}}^{\mathrm{T}}\left(\mathit{\xi }\right)\mathit{\Psi }(\mathit{\xi },\mathit{s},{\mathit{\mu }}_{\mathrm{CTLS}})\mathit{\phi }(\mathit{\xi },\mathit{s},{\mathit{\mu }}_{\mathrm{CTLS}})\hfill \\ +{\left({\frac{\partial \mathrm{vec}\left(\mathit{\Psi }(\mathit{\xi },\mathit{s},\mathit{u})\right)}{\partial {\mathit{\mu }}^{\mathrm{T}}}|}_{\mathit{\mu }={\mathit{\mu }}_{\mathrm{CTLS}}}\right)}^{\mathrm{T}}\hfill \\ (\mathit{\phi }(\mathit{\xi },\mathit{s},{\mathit{\mu }}_{\mathrm{CTLS}})\otimes \mathit{\phi }(\mathit{\xi },\mathit{s},{\mathit{\mu }}_{\mathrm{CTLS}}))={\mathit{0}}_K\hfill \end{array}$$

(47)

The estimate can be expressed as the true value plus error

$$\mathrm{\Delta }{\mathit{\mu }}_{\mathrm{CTLS}}={\mathit{\mu }}_{\mathrm{CTLS}}-{\mathit{\mu }}^o$$

(48)

We substitute equation (48) to $\mathit{\phi }(\mathit{\xi },\mathit{s},\mathit{\mu })$ in equation (40) and ignoring the second-order error term gives

$$\begin{array}{c}\mathit{\phi }(\mathit{\xi },\mathit{s},{\mathit{\mu }}_{\mathrm{CTLS}})\stackrel{·}{=}\mathit{A}\left({\mathit{\xi }}^o\right){\mathit{\mu }}^o-\mathit{b}({\mathit{\xi }}^o,{\mathit{s}}^o)\hfill \\ +\mathit{A}\left({\mathit{\xi }}^o\right)\mathrm{\Delta }{\mathit{\mu }}_{\mathrm{CTLS}}+({\mathit{B}}_1({\mathit{\xi }}^o,{\mathit{s}}^o)-{\mathit{E}}_1({\mathit{\xi }}^o,{\mathit{\mu }}^o)){\mathit{e}}_1\hfill \\ +{\mathit{B}}_2\left({\mathit{\xi }}^o\right){\mathit{e}}_2\hfill \\ =\mathit{A}\left({\mathit{\xi }}^o\right)\mathrm{\Delta }{\mathit{\mu }}_{\mathrm{CTLS}}-\mathit{\Phi }({\mathit{\xi }}^o,{\mathit{s}}^o,{\mathit{\mu }}^o)\mathit{e}\hfill \end{array}$$

(49)

Applying equation (49) and omitting the second term in equation (47) for the second-order error terms are small to be neglected, we have the approximation

$${\mathit{0}}_K\approx {\mathit{A}}^{\mathrm{T}}\left({\mathit{\xi }}^o\right)\mathit{\Psi }({\mathit{\xi }}^o,{\mathit{s}}^o,{\mathit{\mu }}^o)(\mathit{A}\left({\mathit{\xi }}^o\right)\mathrm{\Delta }{\mathit{\mu }}_{\mathrm{CTLS}}-\mathit{\Phi }({\mathit{\xi }}^o,{\mathit{s}}^o,{\mathit{\mu }}^o)\mathit{e})$$

(50)

Then, the estimation error can be decomposed as

$$\begin{array}{cc}\multicolumn{1}{c}{\mathrm{\Delta }{\mathit{\mu }}_{\mathrm{CTLS}}}& \approx {\left({\mathit{A}}^{\mathrm{T}}\left({\mathit{\xi }}^o\right)\mathit{\Psi }({\mathit{\xi }}^o,{\mathit{s}}^o,{\mathit{\mu }}^o)\mathit{A}\left({\mathit{\xi }}^o\right)\right)}^{-1}\hfill \\ \multicolumn{1}{c}{}& {\mathit{A}}^{\mathrm{T}}\left({\mathit{\xi }}^o\right)\mathit{\Psi }({\mathit{\xi }}^o,{\mathit{s}}^o,{\mathit{\mu }}^o)\mathit{\Phi }({\mathit{\xi }}^o,{\mathit{s}}^o,{\mathit{\mu }}^o)\mathit{e}\hfill \end{array}$$

(51)

Multiplying equation (51) with its transpose and taking expectation yields the covariance matrix as

$$\begin{array}{c}\mathbf{cov}\left(\mathrm{\Delta }{\mathit{\mu }}_{\mathrm{CTLS}}\right)=\mathrm{E}\left[\mathrm{\Delta }{\mathit{\mu }}_{\mathrm{CTLS}}\mathrm{\Delta }{\mathit{\mu }}_{\mathrm{CTLS}}^{\mathrm{T}}\right]\hfill \\ ={\left({\mathit{A}}^{\mathrm{T}}\left({\mathit{\xi }}^o\right)\mathit{\Psi }({\mathit{\xi }}^o,{\mathit{s}}^o,{\mathit{\mu }}^o)\mathit{A}\left({\mathit{\xi }}^o\right)\right)}^{-1}\hfill \\ ={\left({\mathit{A}}^{\mathrm{T}}\left({\mathit{\xi }}^o\right){\left(\mathit{\Phi }({\mathit{\xi }}^o,{\mathit{s}}^o,{\mathit{\mu }}^o)\mathit{Q}{\mathit{\Phi }}^{\mathrm{T}}({\mathit{\xi }}^o,{\mathit{s}}^o,{\mathit{\mu }}^o)\right)}^{-1}\mathit{A}\left({\mathit{\xi }}^o\right)\right)}^{-1}\hfill \end{array}$$

(52)

In the following pages, we would observe the asymptotic performance of the CTLS method. Applying ${\mathit{\xi }}^o={\mathit{\xi }}^o({\mathit{s}}^o,{\mathit{\mu }}^o)$ to equation (9) gives

$$\mathit{A}\left({\mathit{\xi }}^o({\mathit{s}}^o,{\mathit{\mu }}^o)\right){\mathit{\mu }}^o=\mathit{b}({\mathit{\xi }}^o({\mathit{s}}^o,{\mathit{\mu }}^o),{\mathit{s}}^o)$$

(53)

Computing the Jacobi matrix of both sides with respect to ${\mathit{\mu }}^o$ , equation (53) becomes

$$\mathit{A}\left({\mathit{\xi }}^o\right)+{\mathit{E}}_1(\mathit{\xi },{\mathit{\mu }}^o){\mathit{\Omega }}_1={\mathit{B}}_1({\mathit{\xi }}^o,{\mathit{s}}^o){\mathit{\Omega }}_1$$

(54)

Based on equations (37) and (54), $\mathit{A}({\mathit{\xi }}^o)$ is given by

$$\mathit{A}\left({\mathit{\xi }}^o\right)=({\mathit{B}}_1({\mathit{\xi }}^o,{\mathit{s}}^o)-{\mathit{E}}_1({\mathit{\xi }}^o,{\mathit{\mu }}^o)){\mathit{\Omega }}_1={\mathit{\Phi }}_1({\mathit{\xi }}^o,{\mathit{s}}^o,{\mathit{\mu }}^o){\mathit{\Omega }}_1$$

(55)

Again calculating the Jacobi matrix with respect to ${\mathit{s}}^o$ , equation (53) is rewritten as

$${\mathit{E}}_1(\mathit{\xi },{\mathit{\mu }}^o){\mathit{\Omega }}_2={\mathit{B}}_1({\mathit{\xi }}^o,{\mathit{s}}^o){\mathit{\Omega }}_2+{\mathit{B}}_2\left({\mathit{\xi }}^o\right)$$

(56)

With some algebraic manipulations, we have the equality relationship

$$-{\mathit{\Phi }}_1^{-1}({\mathit{\xi }}^o,{\mathit{s}}^o,{\mathit{\mu }}^o){\mathit{\Phi }}_2\left({\mathit{\xi }}^o\right)={\mathit{\Omega }}_2$$

(57)

Multiplying equation (51) with its transpose, taking expectation and using equations (55) and (57) to simplify yields

$$\begin{array}{c}\mathbf{cov}\left(\mathrm{\Delta }{\mathit{\mu }}_{\mathrm{CTLS}}\right)={\left({\mathit{A}}^{\mathrm{T}}\left({\mathit{\xi }}^o\right){({\mathit{\Phi }}_1({\mathit{\xi }}^o,{\mathit{s}}^o,{\mathit{\mu }}^o){\mathit{Q}}_1{\mathit{\Phi }}_1^{\mathrm{T}}({\mathit{\xi }}^o,{\mathit{s}}^o,{\mathit{\mu }}^o)+{\mathit{\Phi }}_2\left({\mathit{\xi }}^o\right){\mathit{Q}}_2{\mathit{\Phi }}_2^{\mathrm{T}}\left({\mathit{\xi }}^o\right))}^{-1}\mathit{A}\left({\mathit{\xi }}^o\right)\right)}^{-1}\hfill \\ ={\left({\mathit{\Omega }}_1^{\mathrm{T}}{\mathit{\Phi }}_1^{\mathrm{T}}({\mathit{\xi }}^o,{\mathit{s}}^o,{\mathit{\mu }}^o){({\mathit{\Phi }}_1({\mathit{\xi }}^o,{\mathit{s}}^o,{\mathit{\mu }}^o){\mathit{Q}}_1{\mathit{\Phi }}_1^{\mathrm{T}}({\mathit{\xi }}^o,{\mathit{s}}^o,{\mathit{\mu }}^o)+{\mathit{\Phi }}_2\left({\mathit{\xi }}^o\right){\mathit{Q}}_2{\mathit{\Phi }}_2^{\mathrm{T}}\left({\mathit{\xi }}^o\right))}^{-1}{\mathit{\Phi }}_1({\mathit{\xi }}^o,{\mathit{s}}^o,{\mathit{\mu }}^o){\mathit{\Omega }}_1\right)}^{-1}\hfill \\ ={\left({\mathit{\Omega }}_1^{\mathrm{T}}{({\mathit{Q}}_1+{\mathit{\Phi }}_1^{-1}({\mathit{\xi }}^o,{\mathit{s}}^o,{\mathit{\mu }}^o){\mathit{\Phi }}_2\left({\mathit{\xi }}^o\right){\mathit{Q}}_2{\mathit{\Phi }}_2^{\mathrm{T}}\left({\mathit{\xi }}^o\right){\mathit{\Phi }}_1^{-\mathrm{T}}({\mathit{\xi }}^o,{\mathit{s}}^o,{\mathit{\mu }}^o))}^{-1}{\mathit{\Omega }}_1\right)}^{-1}\hfill \\ ={\left({\mathit{\Omega }}_1^{\mathrm{T}}{({\mathit{Q}}_1+{\mathit{\Omega }}_2{\mathit{Q}}_2{\mathit{\Omega }}_2^{\mathrm{T}})}^{-1}{\mathit{\Omega }}_1\right)}^{-1}\hfill \\ =\mathrm{CRLB}\left(\mathit{\mu }\right)\hfill \end{array}$$

(58)

As a result, the proposed method is equipped with the asymptotic optimal performance and achieves the CRLB over the small error region.

Simulations

In the following simulations, the proposed method (labeled as “CTLS”) is compared with the total least square method (labeled as “TLS”) and the CRLB with synchronization errors is also plotted in picture. It is needed to be illustrated that both the measurement noises and sensor position errors are Gaussian distributed. In addition, synchronization errors rooting in sensor clock biases are constant.

Estimate performance of position and clock biases

In this subsection, we examine the estimate accuracy of source position in 3D scenario which consists of three sensors, and their locations are listed in Table 1. There exist two radiation sources located at ${\mathit{u}}^o$ . In fact, sensors 2 and 3 have clock biases relative to sensor 1 (the reference sensor), which can be collected as ${\mathit{\delta }}^{\mathit{o}}$ . We first fix the noise power of sensor position as ${\sigma }_{\mathrm{s}}^2$ , and claim ${\sigma }_{\mathrm{AOA}}^2$ and ${\sigma }_{\mathrm{TDOA}}^2$ , where ${\sigma }_{\mathrm{m}}^2$ represents the measurement noises. Conduct $N$ Monte Carlo experiments and estimation accuracy is evaluated in terms of $\mathrm{RMSE}(\mathit{\mu })=(\sum _{n=1}^N{\Vert \mathit{\mu }-{\mathit{\mu }}^0\Vert }^2/N{)}^{1/2}$ , which is shown in log scale as the noise power increases. (The $10·\log 10(·)$ operator is utilized for explicit when both position and velocity errors are curved in same plot.) Detailed parameter settings are shown in Table 3.

Table 3.

Simulation parameter settings when measurement errors change.

Parameter	Simulation value
${\mathit{u}}^o$	$[550,500,400{]}^{\mathrm{T}}(\mathrm{m})$ (near-field)
	$[1200,2000,500{]}^{\mathrm{T}}(\mathrm{m})$ (far-field)
${\mathit{\delta }}^{\mathit{o}}$	$[80,260{]}^{\mathrm{T}}(\mathrm{ns})$
${\sigma }_{\mathrm{s}}^2$	$100({\mathrm{m}}^2)$
${\sigma }_{\mathrm{AOA}}^2$	${10}^{-4}{\sigma }_{\mathrm{m}}^2(\mathrm{ra}{\mathrm{d}}^2)$
${\sigma }_{\mathrm{TDOA}}^2$	${\sigma }_{\mathrm{m}}^2({\mathrm{m}}^2)$
${\sigma }_{\mathrm{m}}^2$	−20 to 20 dB
$N$	500

The RMSEs of position and synchronization errors versus measurement noises are shown in Figures 4 and 5, respectively, from which we can see that both TLS and CTLS methods attain the CRLB at the low noise level. As the measurement of noises increases, the TLS method performs much worse and is the first to departure the CRLB, which can be illuminated as the threshold effect generally appearing early in the closed method. It is evident that the CTLS method shows superior performance than the TLS method especially for a near-field source, which achieves CRLB in the presence of synchronization errors and its theoretical RMSE coincides with CRLB perfectly before the thresholding effect occurs, which verifies the former analysis.

Figure 4.

Comparisons of source position RMSEs of the proposed method and the TLS method with the CRLB as measurement noises vary: (a) near-field localization scenario and (b) far-field localization scenario.

Figure 5.

Comparisons of synchronization error RMSEs of the proposed method and the TLS method with the CRLB as measurement noises vary: (a) near-field localization scenario and (b) far-field localization scenario.

In the following, we examine the estimate performance of source position and sensor clock bias by comparing Figure 4 with Figure 5, whose ordinates suggest that the position accuracy exceeds the estimate accuracy of sensor clock bias. One can derive that the CTLS method focuses more efforts on localization rather than estimating the synchronization errors to correct the sensor clock for the synchronization errors only hide in TDOA measurement and its real value is dissolved by one-order derivation.

Next, simulation parameters settings are shown in Table 4, and results are depicted in Figures 6 and 7, which screen that there are big gaps between the RMSE of TLS estimator and CRLB in low sensor position errors no matter how far the source occurs. This difference elucidates that the estimation performance is sensitive to measurement noises or synchro errors under small sensor position errors. As ${\sigma }_{\mathrm{s}}^2$ increases, the RMSE of TLS estimator is close to CRLB for sensor position errors play more influential role in the estimator. As for the CTLS method, in the near-field scenario, the simulation results agree with the theoretical RMSE curve, which coincides with CRLB very well. In far-field scenario, the proposed method cannot achieve CRLB when sensor position errors are small but performs better than the TLS method. It can be interpreted that clock bias is small relative to TDOA measurements and cannot be accurately identified in the localization model. With the increase in ${\sigma }_{\mathrm{s}}^2$ , the CTLS reaches CRLB for the sensor position errors dominant the position accuracy.

Table 4.

Simulation parameter settings when sensor position errors change.

Parameter	Simulation value
${\mathit{u}}^o$	$[550,500,400{]}^{\mathrm{T}}(\mathrm{m})$ (near-field)
	$[1200,2000,500{]}^{\mathrm{T}}(\mathrm{m})$ (far-field)
${\mathit{\delta }}^{\mathit{o}}$	$[80,260{]}^{\mathrm{T}}(\mathrm{ns})$
${\sigma }_{\mathrm{s}}^2$	−10 to 30 dB
${\sigma }_{\mathrm{AOA}}^2$	${10}^{-3.5}{\sigma }_{\mathrm{m}}^2(\mathrm{ra}{\mathrm{d}}^2)$
${\sigma }_{\mathrm{TDOA}}^2$	${10}^{0.5}({\mathrm{m}}^2)$
${\sigma }_{\mathrm{m}}^2$	5 dB
$N$	500

Figure 6.

Comparisons of source position RMSEs of the proposed method and the TLS method with the CRLB as sensor position noises vary: (a) near-field localization scenario and (b) far-field localization scenario.

Figure 7.

Comparisons of synchronization error RMSEs of the proposed method and the TLS method with the CRLB as sensor position noises vary: (a) near-field localization scenario and (b) far-field localization scenario.

Sensitivity to the initial estimate

In this subsection, the sensitivity to the initial estimate of the proposed method is demonstrated. We consider a 3D scenario with the same simulation conditions in part A. We assume that the exact position of radiation source is unknown but distributed randomly around the true value. Thus, we choose the initial estimate as ${\mathit{\mu }}^{(0)}={\mathit{\mu }}^o+\mathit{\Delta }$ , where $\mathit{\Delta }$ represents the initial estimate error and is submitted to zero-mean Gaussian distribution with covariance ${\sigma }_0^2$ . Parameter settings are shown in Table 5. Figures 8 and 9 are the results of RMSEs versus ${\sigma }_0^2$ of the proposed method. It is evident that the proposed method coincides with CRLB in theory when the noise power is moderate and its simulation results fluctuate at the curve of CRLB before the threshold emerging at ${\sigma }_0^2=20\mathrm{dB}$ . This phenomenon reveals that the proposed method is robust to the initial estimate, and we can make an interpretation that this robustness is derived from the joint of AOA and TDOA-based localization. More analytically, the non-linearization in the wake of TDOA is commendably transforming to linearization problem resort to AOA. Therefore, the global optimum can be attained by iterative Newton’s method and the initial estimate error is alleviated.

Table 5.

Simulation parameter settings when synchro errors change.

Parameter	Simulation value
${\mathit{u}}^{\mathit{o}}$	$[550,500,400{]}^{\mathrm{T}}(\mathrm{m})$
${\mathit{\delta }}^{\mathit{o}}$	${[80,260]}^{\mathrm{T}}(\mathrm{ns})$
${\sigma }_{\mathrm{s}}^2$	$100({\mathrm{m}}^2)$
${\sigma }_{\mathrm{AOA}}^2$	${10}^{-4}{\sigma }_{\mathrm{m}}^2(\mathrm{ra}{\mathrm{d}}^2)$
${\sigma }_{\mathrm{TDOA}}^2$	${\sigma }_{\mathrm{m}}^2({\mathrm{m}}^2)$
${\sigma }_{\mathrm{m}}^2$	5 dB
$N$	300
${\sigma }_0^2$	−20 to 25 dB

Figure 8.

Sensitivity of simulation and theoretical position accuracy of the CTLS method compared with CRLB as initial estimate errors vary.

Figure 9.

Sensitivity of simulation and theoretical synchronization error accuracy of the CTLS method compared with CRLB as initial estimate errors vary.

Convergence speed of the proposed method

We concentrate on the effect of the convergence speed of the proposed method in this subsection. The localization scenario and simulation conditions are same as section “Sensitivity to the initial estimate.” We let ${\sigma }_0^2$ vary from 50 to 200 with 50 intervals and each of them corresponds to an initial estimate. We execute the proposed method proceeding from each initial estimate with 300 runs and compute the mean of the cost function in iterations. The convergence performance is evaluated in terms of this mean, and simulation results are exhibited in Figure 10, where the number of iteration times is varying from 1 to 10. It is obvious that the proposed method generally converges to the global minimum with three iterations. It has been proved that iterative Newton’s method is a powerful tool to solve the CTLS method for the rapid convergence rate.

Figure 10.

Convergence rate of iterative Newton’s method employed to solve the CTLS method.

Performance improvement when introducing synchro errors

In this subsection, we compare the position performance of the proposed method considering the synchro errors with the WLS method in the study by Chen and Ho ⁷ not considering the errors in the non-synchronization scenario to demonstrate the improvement. CRLBs with and without synchro errors are also plotted for comparison except the RMSE derived in section “CRLB and RMSE,” which make no sense as a contrast. It can be explained that the RMSE just assumes no synchro errors in the localization model but the WLS method does ignore not only synchro errors but also sensor position errors. In order to meet the condition of the WLS method, we consider a 3D localization scenario consisting of the first two sensors in Table 1, which are non-synchronization with a clock bias ${\delta }_0$ . Other parameter settings are shown in Table 3, and Figure 11 shows the results when ${\delta }_0=260\mathrm{ns}$ . In the near-field localization scenario, the theoretical RMSE curve is consistent with CRLB. The simulation RMSE slightly fluctuates around CRLB and departures up to ${\sigma }_{\mathrm{m}}^2=20\mathrm{dB}$ . The TLS method achieves the CRLB in low noise levels and begins to separate when ${\sigma }_{\mathrm{m}}^2=5\mathrm{dB}$ . The WLS method assuming no synchro errors performs weakly in low measurement noises, which can be explained that the clock bias is much larger than the measurement noises and this method suffers mismatch for the dominate synchro errors. As the ${\sigma }_{\mathrm{m}}^2$ increases, the synchro errors may offset be the large measurement noises and the WLS method attains higher accuracy than the CTLS and TLS methods introducing synchro errors into the localization model. In the far-field localization scenario, the CTLS and TLS methods have similar performances and both of them achieve CRLB except a little divergence in large noises. The WLS method performs better than the other two methods in small measurement noises and is even lower than CRLB. It can be interpreted that synchro error among the two sensors is very small compared with the TDOA measurements and can be ignored in a far-field localization scenario. Hence, the WLS method with smaller DOF of estimate parameters (only position no clock bias) can obtain higher accuracy. With the increase in ${\sigma }_{\mathrm{m}}^2$ , the measurement noises are dominant in localization accuracy and the methods perform similarly with and without considering synchro errors.

Figure 11.

Comparisons of source position RMSEs of the method considering and not considering the synchro errors with the CRLB as measurement noises vary with a 260 ns clock bias: (a) near-field localization scenario and (b) far-field localization scenario.

Figure 12 shows the results when ${\delta }_0=1600\mathrm{ns}$ , where the WLS method fails to determine the near-field source position. However, the TLS and CTLS methods are still practicable in such a large clock bias. In far-field scenario, this clock bias cannot be ignored and played a major role in localization. The simulation results and the reasons are same as Figure 11(a) for the role of ${\delta }_0=1600\mathrm{ns}$ in far-field scenario is similar to ${\delta }_0=260\mathrm{ns}$ in near-field. As a result, only when clock bias is large enough compared with TDOA measurements, the model introducing synchro errors can bring evident performance improvement and the proposed method can alleviate the performance loss caused by synchro errors.

Figure 12.

Comparisons of source position RMSEs of the method considering and not considering the synchro errors with the CRLB as measurement noises vary with a 1600 ns clock bias: (a) near-field localization scenario and (b) far-field localization scenario.

Performance of random distributed source

In the actual scenario, the source position is uncertain. Hence, we conduct this simulation to prove that the proposed method is valid wherever the source locates. We arrange the near-field source and far field distributed at a limited area with a radius of 600 and 6000 m, respectively. Other parameters are same as Table 1, and simulation results are shown in Figures 13 and 14. From the results, we can see that the proposed method still performs better than the TLS method in both near- and far-field localization scenarios.

Figure 13.

Comparisons of source position RMSEs of the proposed method and the TLS method with the CRLB as measurement noises vary: (a) near-field localization scenario and (b) far-field localization scenario.

Figure 14.

Comparisons of synchronization error RMSEs of the proposed method and the TLS method with the CRLB as measurement noises vary: (a) near-field localization scenario and (b) far-field localization scenario.

Conclusion

A system of CTLS localization method based on AOA and TDOA measurements in the presence of sensor position errors and synchronization errors has been exhibited in this article. The study begins with the derivation of CRLB and MSE aiming to demonstrate the significance of taking synchronization errors into consideration. Subsequently, we construct a set of pseudo-linear equations, based on which a CTLS localization model has been formulated and Newton iteration is resorted to solve it. Simulation results illustrate that the proposed method can achieve the CRLB under the moderate noises, which verify the validity of theoretical derivation and the superiority of the new algorithm.