Research on multi-moving target location algorithm based on improved TDOA/FDOA

Abstract

Wireless positioning and tracking technology are key technologies in applying wireless sensor networks and have vital research significance and application value. In this paper, in terms of positioning multiple non-coincident dynamic targets, aiming at the problem of low precision and large computational complexity when locating multiple dynamic targets by the two-step least squares algorithm and the constrained total least squares algorithm, an improved constrained total least squares algorithm is proposed. This algorithm fully considers the constraints, introduces the Lagrangian multiplier technology and the quasi-Newton BFGS iterative formula, avoids the calculation of the Hessian matrix, reduces the amount of calculation, and improves the positioning accuracy. Simulation experiments show that when the measurement error and sensor position error are moderate, the ICTLS positioning algorithm has a smaller RMSE than the TSWLS and CTLS positioning algorithms, showing higher positioning accuracy and stronger robustness. Secondly, aiming at the problem that the target is close to the reference node or any coordinate axis causes the positioning error of the traditional positioning algorithm to increase sharply, an optimized two-step least squares algorithm is proposed. This algorithm corrects the defects of the TSWLS algorithm by selecting the reference station and rotating the coordinate system again, and improves the positioning performance of the algorithm while reducing the amount of calculation.

Introduction

In real life, Wireless Sensor Network (WSN) is widely used and it plays an increasingly important role, and location-based devices are the key core part of all applications, while the positioning target generally occurs in a noisy environment. The measurement error is large, and the error is also large when a single TDOA (Time Difference Of Arrival) ranging technology is used to achieve target positioning. ¹ To improve the algorithm’s positioning accuracy, hybrid positioning technology is usually used to estimate the position information of the target. ² When carrying out position information estimation for multiple moving targets simultaneously, it is usually necessary to use a positioning model combined with TDOA and FDOA (Frequency Difference Of Arrival) to estimate the target speed and position, which effectively improves the positioning accuracy of moving targets. ³ However, when positioning a dynamic target, the position error of the sensor and the measurement error will bring a large positioning deviation. ⁴ Sun et al. used the TSWLS (Two Step Weighted Least Squares) algorithm to locate multiple moving targets based on the TDOA/FDOA measurement technology combined with the sensor’s position error. The experimental results show that the RMSE (Root Mean Squares Error) of the position estimation can reach CRLB (Cramer-RaoLower Bound) ⁵ only under the condition of the high signal-to-noise ratio. Since the TSWLS algorithm has poor adaptability to noise errors, the positioning deviation is large, and the positioning accuracy is limited. Chen Shaochang et al. proposed in the literature ⁶ that the constrained total least squares algorithm considers all coefficient matrices’ noise in the positioning equation and obtains the target’s position coordinates through the Newton iteration method. The experimental results show that this algorithm has better performance than the TSWLS algorithm. However, this algorithm estimates the location information of fixed targets based on TDOA ranging technology and is not suitable for locating moving targets. By analyzing the factors affecting the positioning performance of the TSWLS algorithm, Liu Yang et al. proposed a moving target positioning algorithm based on TDOA/FDOA to correct the positioning error. Without increasing the computational complexity, the positioning error was reduced and the ability to adapt to measurement noise is improved. Yanbin Zou et al. proposed a semi-definite programming (SDP) method, which transforms the common MLE (maximum likelihood estimator) problem into a convex optimization problem, which iterates with the position and velocity obtained by the SDP method as the initial values and updates the velocity using weighted least squares to update the position using SDP, the main advantage of this scheme is that the localization performance is more obvious under medium to high noise conditions, ⁷ but the position error of the sensor and the sharp decrease in localization performance when the target is close to any axis of the set reference sensor are not considered. In contrast, this paper is mainly based on TSWLS and CTLS, taking into account the sensor position error and measurement error, as well as the localization problem in special scenarios.

Based on the above analysis, this paper mainly studies the positioning of multiple moving targets from two aspects. The first aspect of the research is: Aiming at the positioning of multiple non-coincident dynamic targets, ICTLS (Improved CTLS, ICTLS) positioning is proposed on the basis of the TSWLS and CTLS positioning algorithms. The ICTLS positioning algorithm corrects two defects in the TSWLS algorithm: ① First, the deviation of WLS estimation results increases with the increase of noise error; ② Second, the nonlinear operation introduced in WLS causes large calculation errors, thereby reducing the accuracy of the positioning algorithm. The ICTLS algorithm proposed in the article corrects the connection between the additional variables introduced in the CTLS algorithm and the target position coordinates, and establishes global constraints based on the connection between the target and the additional variables; At the same time, the Lagrangian multiplier technology is introduced to solve the positioning equation. The BFGS iterative formula of the quasi-Newton method avoids the calculation of the Hessian matrix, reduces the amount of calculation and speeds up the convergence speed. Simulation experiments show that when the measurement error and sensor position error are moderate, the ICTLS positioning algorithm has a smaller RMSE than the TSWLS and CTLS positioning algorithms, showing higher positioning accuracy and stronger robustness.

The second aspect of the research is: in view of the fatal flaw in a variety of existing positioning algorithms, that is, when the target approaches a reference sensor or a certain coordinate axis, it will cause a sudden increase in the estimation error of the target position based on the TSWLS algorithm. For this problem, a revised TSWLS positioning algorithm was proposed. The correction algorithm first uses the first step of the TSWLS algorithm to initially estimate the target position, and then re-selects the reference sensor (in principle, the new reference sensor farthest from the target) and the rotation axis according to the estimated target position. Even if the target is approaching the reference sensor or a certain coordinate axis, the correction algorithm also shows superior positioning performance, which makes up for the loopholes in the classic TSWLS series of algorithms.

Algorithm model and positioning principle

Positioning scene

K non-coincident moving targets are randomly distributed in the three-dimensional space, and the true position and speed of the ith target are recorded as ${\mathit{u}}_{\mathit{i}}$ and ${\stackrel{·}{\mathit{u}}}_{\mathit{i}},i=1,2,\cdots ,K$ , respectively. M sensors are deployed in the same three-dimensional space, ${\mathit{s}}^{\mathit{0}}={[{\mathit{s}}_{\mathit{1}}^{0\mathrm{T}}{,s}_{\mathit{2}}^{0\mathrm{T}},\dots {,s}_M^{0\mathrm{T}}]}^{\mathrm{T}}$ is the real position of the sensor, ${\stackrel{·}{\mathit{s}}}^0={[{\stackrel{·}{\mathit{s}}}_1^{0\mathrm{T}}{,\stackrel{·}{\mathit{s}}}_{\mathit{2}}^{0\mathrm{T}},\dots {,\stackrel{·}{\mathit{s}}}_M^{0\mathrm{T}}]}^{\mathrm{T}}$ is the real speed of the sensor, where ${\mathit{s}}_j^0$ and ${\stackrel{·}{\mathit{s}}}_j^0$ are the real position and speed of the jth sensor, $[·{]}^0$ is the real value of $[·]$ , and $[·{]}^{\mathrm{T}}$ is the real value of $[·]$ transpose operation. In the article, it is assumed that the true position information of the sensor is unknown, but the measured values ${\mathit{s}}_j$ and ${\stackrel{·}{\mathit{s}}}_j$ with position error information are known, and $\mathit{s}={[{\mathit{s}}_1^{\mathrm{T}},{\mathit{s}}_2^{\mathrm{T}},\cdots ,{\mathit{s}}_M^{\mathrm{T}}]}^{\mathrm{T}}={\mathit{s}}^0+\mathrm{\Delta }\mathit{s}$ and $\stackrel{·}{\mathit{s}}={[{\stackrel{·}{\mathit{s}}}_1^{\mathrm{T}},{\stackrel{·}{\mathit{s}}}_2^{\mathrm{T}},\cdots ,{\stackrel{·}{\mathit{s}}}_M^{\mathrm{T}}]}^{\mathrm{T}}={\stackrel{·}{\mathit{s}}}^0+\mathrm{\Delta }\stackrel{·}{\mathit{s}}$ are the known sensor position and velocity vector, and $\mathrm{\Delta }\mathit{s}$ and $\mathrm{\Delta }\stackrel{·}{\mathit{s}}$ are the sensor’s position and velocity, respectively, the error vector. The vector of real information about sensor position and velocity is ${\mathit{\beta }}^0={[{\mathit{s}}^{0\mathrm{T}},{\stackrel{·}{\mathit{s}}}^{0\mathrm{T}}]}^{\mathrm{T}}$ . Similarly, the measurement vector containing position error is $\mathit{\beta }={[{\mathit{s}}^{\mathrm{T}},{\stackrel{·}{\mathit{s}}}^{\mathrm{T}}]}^{\mathrm{T}}={\mathit{\beta }}^0+\mathbf{\Delta }\mathit{\beta }$ , and $\mathrm{\Delta }\mathit{\beta }={[{\mathit{s}}^{\mathrm{T}},{\stackrel{·}{\mathit{s}}}^{\mathrm{T}}]}^{\mathrm{T}}$ is the position error vector. Also assume that the position error vector is a zero-mean Gaussian random process and the covariance matrix is $Q_{\beta }$ . For the convenience of calculation, the sensor with j = 1 is generally regarded as the reference sensor. Therefore, the true distance between target i and sensor j is:

$$r_{j,i}^0=\Vert {\mathit{u}}_i-{\mathit{s}}_j^0{\Vert }_2=\sqrt{{({\mathit{u}}_i-{\mathit{s}}_j^0)}^{\mathrm{T}}({\mathit{u}}_i-{\mathit{s}}_j^0)}$$

(1)

In the formula, $\Vert •{\Vert }_2$ represents the two norms. The real and measured values of the distance difference between target i and sensors 1 and j are

$$r_{j1,i}^0=r_{j,i}^0-r_{1,i}^0,r_{j1,i}=r_{j1,i}^0+\mathrm{\Delta }{r}_{j1,i}$$

(2)

Where $\mathrm{\Delta }{r}_{j1,i}$ is the measurement error. Therefore, the measurement vector in the TDOA form of the position information of the target i is recorded as:

$${\mathit{r}}_i={[r_{21,i},r_{31,i},\cdots ,r_{M1,i}]}^{\mathrm{T}}={\mathit{r}}_i^0+\mathrm{\Delta }{\mathit{r}}_i$$

(3)

Where $\mathrm{\Delta }{\mathit{r}}_i={[\mathrm{\Delta }{r}_{21,i},\mathrm{\Delta }{r}_{31,i},\cdots ,\mathrm{\Delta }{r}_{M1,i}]}^{\mathrm{T}}$ is the measurement error vector.

When the target and the sensor move relative to each other, the FDOA measurement value of the target speed information can be obtained from equation (2):

$${\stackrel{·}{r}}_{j1,i}={\stackrel{·}{r}}_{j1,i}^0+\mathrm{\Delta }{\stackrel{·}{r}}_{j1,i}={\stackrel{·}{r}}_{j,i}^0-{\stackrel{·}{r}}_{1,i}^0+\mathrm{\Delta }{\stackrel{·}{r}}_{j1,i}$$

(4)

Among them, ${\stackrel{·}{r}}_{j1,i}^0$ is the true distance change rate between target i and sensor j, and ${\stackrel{·}{r}}_{j,i}^0={\rho }_{u_i,s_j^0}^{\mathrm{T}}({\stackrel{.}{\mathit{u}}}_i-{\stackrel{·}{\mathit{s}}}_j^0)({\rho }_{\mathit{a},\mathit{b}}=(\mathit{a}-\mathit{b})/\parallel \mathit{a}-\mathit{b}\parallel$ is the unit vector pointing from b to a ). According to formula (3), In the same way, the FDOA measurement vector about the target speed information is obtained: ${\stackrel{•}{\mathit{r}}}_i={[{r^{•}}_{21,i},{r^{•}}_{31,i},\cdots ,{r^{•}}_{M1,i}]}^{\mathrm{T}}=\stackrel{•}{{\mathit{r}}_i^0}+\mathrm{\Delta }\stackrel{•}{{\mathit{r}}_i}$ .

For the convenience of the following description, the TDOA and FDOA measurement vectors of the target information are often expressed as a vector equation.^8,9 Since the measurement vector of target i can be denoted as ${\mathbf{\alpha }}_{\mathbf{i}}={[{\mathit{r}}_i^{\mathrm{T}},{\stackrel{·}{\mathit{r}}}_i^{\mathrm{T}}]}^{\mathrm{T}}={\mathbf{\alpha }}_{\mathbf{i}}^0+\mathbf{\Delta }{\mathbf{\alpha }}_{\mathbf{i}}$ and $\mathrm{\Delta }{\mathbf{\alpha }}_{\mathbf{i}}={[\mathrm{\Delta }{\mathit{r}}_i^{\mathrm{T}},\mathrm{\Delta }{\stackrel{·}{\mathit{r}}}_i^{\mathrm{T}}]}^{\mathrm{T}}$ . Therefore, the measurement equation is denoted as $\mathbf{\alpha }={\left[\begin{array}{cc}{\mathbf{\alpha }}_1,{\mathbf{\alpha }}_2,\hfill & \cdots ,{\mathbf{\alpha }}_{\mathbf{K}}\hfill \end{array}\right]}^{\mathrm{T}}={\mathbf{\alpha }}^0+\mathbf{\Delta }\mathbf{\alpha }$ and the dimension is $2K(M-1)\times 1$ . Also assume that $\mathrm{\Delta }\mathbf{\alpha }$ is a Gaussian random process with zero mean and the covariance matrix is the measurement error vector of $Q_{\alpha }$ . The position error vector $\mathrm{\Delta }\mathit{\beta }$ and the measurement error vector $\mathrm{\Delta }\mathbf{\alpha }$ are independent of each other.

Positioning model

In order to estimate the position and velocity of the target, it is necessary to construct a TDOA/FDOA positioning equation set about the position information of the target and the sensor. To this end, an additional variable $r_{1,i}^0=\Vert {\mathit{u}}_i-{\mathit{s}}_1^0{\Vert }_2$ is introduced, $r_{1,i}^0$ and $r_{j,i}^0$ are substituted into equation (2) and both sides are squared, ¹⁰ and the TDOA equations are obtained after sorting:

$$\begin{array}{c}{\left(r_{j1,i}^0\right)}^2+2r_{j1,i}^0{r}_{1,i}^0={({\mathit{s}}_j^0-{\mathit{s}}_1^0)}^{\mathrm{T}}({\mathit{s}}_j^0-{\mathit{s}}_1^0)\hfill \\ -2{({\mathit{s}}_j^0-{\mathit{s}}_1^0)}^{\mathrm{T}}({\mathit{u}}_i-{\mathit{s}}_1^0)\hfill \end{array}$$

(5)

Since the introduced additional variable $r_{1,i}^0=\Vert {\mathit{u}}_i-{\mathit{s}}_1^0{\Vert }_2$ contains the real position coordinate information of sensor 1, we expand $r_{1,i}^0$ at the sensor noise position ${\mathit{s}}_1$ as a first-order Taylor series, and ignore the higher-order terms above the second order, and get:

$$r_{1,i}^0\approx \Vert {\mathit{u}}_i-{\mathit{s}}_1{\Vert }_2+{\mathit{\rho }}_{u_i,s_1}^{\mathrm{T}}\mathbf{\Delta }{\mathit{s}}_1={\hat{r}}_{1,i}+{\mathit{\rho }}_{{\mathit{u}}_i,s_1}^{\mathrm{T}}\mathbf{\Delta }{\mathit{s}}_1$$

(6)

Where ${\hat{r}}_{1,i}=\Vert {\mathit{u}}_i-{\mathit{s}}_1{\Vert }_2$ is the estimated value of $r_{1,i}^0$ and assume ${\mathit{d}}_{1,i}={\mathit{\rho }}_{{\mathit{u}}_{\mathbf{i}},{\mathit{s}}_1}^{\mathrm{T}}$ . Substituting formula (6) and ${\mathit{s}}_1={\mathit{s}}_1^0+\mathrm{\Delta }{\mathit{s}}_1$ into formula (5), we get:

$$\begin{array}{c}2{({\mathbf{s}}_j-{\mathbf{s}}_1)}^{\mathrm{T}}({\mathit{u}}_i-{\mathbf{s}}_1)+2r_{j1,i}{\hat{r}}_{1,i}-2{(\mathrm{\Delta }{\mathbf{s}}_j-\mathrm{\Delta }{\mathbf{s}}_1)}^{\mathrm{T}}({\mathbf{u}}_i-{\mathbf{s}}_1)\hfill \\ -2\mathrm{\Delta }{r}_{j1,i}{\hat{r}}_{1,i}\hfill \\ ={({\mathbf{s}}_j-{\mathbf{s}}_1)}^T({\mathbf{s}}_j-{\mathbf{s}}_1)+2r_{j1,i}\mathrm{\Delta }{r}_{j1,i}-2{({\mathbf{s}}_j-{\mathbf{s}}_1)}^{\mathrm{T}}\mathrm{\Delta }{\mathbf{s}}_j\hfill \\ -2r_{j1,i}{\mathbf{d}}_{1,i}\mathrm{\Delta }{\mathbf{s}}_1-{\left(r_{j1,i}\right)}^2\hfill \end{array}$$

(7)

Differentiate equation (7) with time to obtain the FDOA equations:

$$\begin{array}{c}({\stackrel{·}{\mathit{s}}}_j-{\stackrel{·}{\mathit{s}}}_1{)}^{\mathrm{T}}({\mathit{u}}_i-{\mathit{s}}_1)+{\stackrel{·}{\mathit{r}}}_{j1,i}{\hat{r}}_{1,i}+{({\mathit{s}}_j-{\mathit{s}}_1)}^{\mathrm{T}}({\stackrel{·}{\mathit{u}}}_i-{\stackrel{·}{\mathit{s}}}_1)+r_{j1,i}{\hat{\overrightarrow{r}}}_{1,i}\hfill \\ -{(\mathrm{\Delta }{\stackrel{·}{\mathrm{s}}}_j-\mathrm{\Delta }{\stackrel{·}{\mathrm{s}}}_1)}^{\mathrm{T}}({\mathit{u}}_i-{\mathit{s}}_1)-\mathrm{\Delta }{\stackrel{·}{\mathit{r}}}_{j1,i}{\hat{r}}_{1,i}-{(\mathrm{\Delta }{\mathit{s}}_j-\mathrm{\Delta }{\mathit{s}}_1)}^{\mathrm{T}}({\stackrel{·}{\mathit{u}}}_i-{\stackrel{·}{\mathit{s}}}_1)\hfill \\ -\mathrm{\Delta }{r}_{j1,i}{\hat{\overrightarrow{r}}}_{1,i}={({\stackrel{·}{\mathrm{s}}}_j-{\stackrel{·}{\mathrm{s}}}_1)}^{\mathrm{T}}({\mathit{s}}_j-{\mathit{s}}_1)-r_{j1,i}{\stackrel{·}{\mathit{r}}}_{j1,i}+{\stackrel{·}{\mathit{r}}}_{j1,i}\mathrm{\Delta }{r}_{j1,i}\hfill \\ +r_{j1,i}\mathrm{\Delta }{\stackrel{·}{\mathit{r}}}_{j1,i}-{({\stackrel{·}{\mathrm{s}}}_j-{\stackrel{·}{\mathrm{s}}}_1)}^{\mathrm{T}}\mathrm{\Delta }{\mathit{s}}_j-{({\mathit{s}}_j-{\mathit{s}}_1)}^{\mathrm{T}}\mathrm{\Delta }{\stackrel{·}{\mathrm{s}}}_j-{\mathit{d}}_{2,i}^{\mathrm{T}}\mathrm{\Delta }{\mathit{s}}_1-{\mathit{d}}_{1,i}^{\mathrm{T}}\mathrm{\Delta }{\stackrel{·}{\mathrm{s}}}_1\hfill \end{array}$$

(8)

Where, ${\mathit{d}}_{2,i}=\frac{({\stackrel{·}{\mathit{u}}}_i-{\stackrel{·}{\mathit{s}}}_1)}{{\hat{r}}_{1,i}}-\frac{({\stackrel{·}{\mathit{u}}}_i-{\stackrel{·}{\mathit{s}}}_1)}{{\hat{r}}_{1,i}}{\rho }_{{\mathit{u}}_i,s_1}{\rho }_{{\mathit{u}}_{i,}}^{\mathrm{T}}$ . Synthesize (7) and (8) to obtain the positioning equations in matrix form:

$$(\mathit{A}+\mathrm{\Delta }\mathit{A})\theta =\mathit{b}+\mathrm{\Delta }\mathit{b}$$

(9)

In the formula, $\theta ={[{\mathit{w}}^{\mathrm{T}},{\hat{r}}_{1,i},{\stackrel{·}{\mathit{w}}}^{\mathrm{T}},{\widehat{\stackrel{·}{\mathit{r}}}}_{1,i}]}^{\mathrm{T}},$ $\mathit{w}={\mathit{u}}_i-{\mathit{s}}_1,\mathit{w}={\stackrel{·}{\mathit{u}}}_i-{\stackrel{·}{\mathit{s}}}_1$ , and the expression of $\mathit{A}$ , $\mathrm{\Delta }\mathit{A}$ , $\mathit{b}$ , $\mathrm{\Delta }\mathit{b}$ is consistent with the corresponding matrix vector expression in the literature. ¹¹

Optimization

For the above method of solving the target positioning equation, the TSWLS positioning algorithm only considers the error of part of the coefficient matrix in the equation, ¹² which causes the target estimated position error to increase; secondly, the TSWLS algorithm considers that the introduced variables ${\hat{r}}_{1,i}$ and ${\widehat{\stackrel{·}{\mathit{r}}}}_{1,i}$ are mutually related to $\mathit{w}$ and $\stackrel{·}{\mathit{w}}$ . Independence, without considering the constraint relationship between them, causing distortion of the positioning model, resulting in increased position estimation deviation. For this reason, the error of all coefficient matrices in the equation is considered in this paper, combined with the constraints between additional variables ${\hat{r}}_{1,i}$ and ${\widehat{\stackrel{·}{\mathit{r}}}}_{1,i}$ and $\mathit{w}$ and $\stackrel{·}{\mathit{w}}$ , ¹³ the CTLS positioning algorithm is used to estimate the position information of the target. Therefore, the value function of the vector $\theta$ to be sought can be constructed as:

$$J(\theta )=(\mathit{A}\theta -b{)}^{\mathrm{T}}{\mathit{W}}_{\theta }(\mathit{A}\theta -b)$$

(10)

In the formula, ${\mathit{W}}_{\theta }$ is the weighting matrix. Therefore, the least squares solution of vector $\theta$ is:

$$\hat{\theta }=\arg \underset{\tilde{\theta }}{min}(\mathit{A}\tilde{\theta }-b{)}^{\mathrm{T}}{\mathit{W}}_{\tilde{\theta }}(\mathit{A}\tilde{\theta }-b)$$

(11)

In the formula, $\hat{\theta }$ is the estimated value of $\theta$ , and $\tilde{\theta }$ is the optimized value of $\theta$ . Since the introduced variables ${\hat{r}}_{1,i}$ and ${\widehat{\stackrel{·}{\mathit{r}}}}_{1,i}$ have a certain constraint relationship with $\mathit{w}$ and $\stackrel{·}{\mathit{w}}$ , it is necessary to solve equation (11) in combination with this constraint. The constraint conditions of the comprehensive formulas (1) and (4) can be expressed as:

$${\hat{r}}_{1,i}=\sqrt{{\mathit{w}}^{\mathrm{T}}\mathit{w}},{\widehat{\stackrel{·}{\mathit{r}}}}_{j,i}={\stackrel{·}{\mathit{w}}}^{\mathrm{T}}\mathit{w}/{\hat{r}}_{1,i}$$

(12)

Substituting the deformed formula (12) into the formula (10), the constraint conditions in matrix form are obtained after finishing:

$${\theta }^{\mathrm{T}}{\Sigma }_1\theta =0,{\theta }^{\mathrm{T}}{\Sigma }_2\theta =0$$

(13)

In the formula, ${\Sigma }_1=\left[\begin{array}{cc}{\Sigma }_{11}\hfill & 0\hfill \\ \multicolumn{1}{c}{\mathit{O}}& \mathit{O}\hfill \end{array}\right],{\Sigma }_2=\left[\begin{array}{cc}\mathit{O}\hfill & {\Sigma }_{12}\hfill \\ \multicolumn{1}{c}{{\Sigma }_{21}}& \mathit{O}\hfill \end{array}\right]$ , $\mathit{O}$ is a 4 × 4 all-zero matrix, ${\mathit{\Sigma }}_{11}={\mathit{\Sigma }}_{12}={\mathit{\Sigma }}_{22}=\mathrm{diag}\{1,1,1,-1\}$ .

Therefore, the total least squares solution of vector $\theta$ becomes the CTLS solution:

$$\begin{array}{c}\hat{\theta }=\arg \underset{\tilde{\theta }}{min}(\mathit{A}\tilde{\theta }-b{)}^{\mathrm{T}}{\mathit{W}}_{\tilde{\theta }}(\mathit{A}\tilde{\theta }-b)\hfill \\ \mathrm{s}.t.{\theta }^{\mathrm{T}}{\Sigma }_1\theta =0,{\theta }^{\mathrm{T}}{\Sigma }_2\theta =0\hfill \end{array}$$

(14)

For solving the problem of constrained optimization minimum of formula (14), the Lagrangian multiplier technique is introduced, and the value function of A (10) becomes:

$$\begin{array}{c}L(\theta ,{\lambda }_1,{\lambda }_2)=(\mathit{A}\theta -\mathit{b}{)}^{\mathrm{T}}{\mathit{W}}_{\theta }(\mathit{A}\theta -\mathit{b})+{\lambda }_1{\theta }^{\mathrm{T}}{\Sigma }_1\theta \hfill \\ +{\lambda }_2{\theta }^{\mathrm{T}}{\Sigma }_2\theta \hfill \end{array}$$

(15)

Where ${\lambda }_1$ and ${\lambda }_2$ are the introduced Lagrange multipliers. Therefore, equation (15) is a quadratic constrained non-convex optimization problem. To solve this problem, first obtain a partial differential of equation (15) with respect to $\theta$ ,^14,15 and then set $\frac{\partial L[\theta ,{\lambda }_1,{\lambda }_2]}{\partial \theta }$ to obtain:

$$\begin{array}{c}2[{(\mathit{A}-\mathit{D})}^{\mathrm{T}}{\mathit{W}}_{\theta }\mathit{A}+{\lambda }_1{\mathit{\Sigma }}_1+{\lambda }_2{\Sigma }_2]\theta -2(\mathit{A}-\mathit{D}{)}^{\mathrm{T}}\hfill \\ {\mathit{W}}_{\theta }\mathit{b}=0\hfill \end{array}$$

(16)

In the formula, $\mathit{D}=[{\mathit{EQF}}_1^{\mathrm{T}}{\mathit{W}}_{\theta }(\mathit{A}\theta -\mathit{b}),\cdots ,{\mathit{EQF}}_{L+1}^{\mathrm{T}}{\mathit{W}}_{\theta }(\mathit{A}\theta -\mathit{b})]$ , L  = 8 is the dimension of vector $\theta$ , $\mathit{E}=\sum _{l=1}^{L+1}{\theta }_l{\mathit{F}}_l-{\mathit{F}}_{L+1},\mathit{Q}=E\left(\mathit{n}{\mathit{n}}^{\mathrm{T}}\right)$ is the initial weighting matrix and $\mathit{n}={[\mathrm{\Delta }{\mathit{\beta }}^{\mathrm{T}},\mathrm{\Delta }{\mathbf{\alpha }}^{\mathrm{T}}]}^{\mathrm{T}}$ is random noise, and ${\mathit{F}}_1,\dots ,{\mathit{F}}_{\mathit{L}}$ is the coefficient matrix of the disturbance matrix $\mathit{\Delta }\mathit{A},\mathit{\Delta }\mathit{b}$ . Therefore, the final estimate of vector $\theta$ is:

$$\hat{\theta }={[{(\mathit{A}-\mathit{D})}^{\mathrm{T}}{\mathit{W}}_{\theta }\mathit{A}+{\lambda }_1{\Sigma }_1+{\lambda }_2{\Sigma }_2]}^{-1}(\mathit{A}-\mathit{D}{)}^{\mathrm{T}}{\mathit{W}}_{\theta }\mathit{b}$$

(17)

Substituting formula (17) into the constraint formula (13) again, we get:

$${\hat{\theta }}^{\mathrm{T}}{\Sigma }_1\hat{\theta }=0,{\hat{\theta }}^{\mathrm{T}}{\Sigma }_2\hat{\theta }=0$$

(18)

Equation (18) composes a high-order polynomial equation set about ${\lambda }_1$ and ${\lambda }_2$ . ${\lambda }_1$ and ${\lambda }_2$ are obtained by solving the equation set and substituted into equation (17) to obtain the final estimated position information of the target. The system of high-order polynomial equations is solved by the quasi-Newton iteration method, ¹⁶ and the formula (18) can be changed to

$$f_k({\lambda }_1,{\lambda }_2)={\mathit{p}}^{\mathrm{T}}{\mathit{G}}_f^{-\mathrm{T}}{\mathit{\Sigma }}_k{\mathit{G}}_f^{-1}\mathit{p}=0,k=1,2$$

(19)

Where, ${\mathit{G}}_f=(\mathit{A}-\mathit{D}{)}^{\mathrm{T}}{\mathit{W}}_{\theta }\mathit{A}+{\lambda }_{\mathrm{I}}{\mathit{\Sigma }}_1+{\lambda }_2{\mathit{\Sigma }}_2$ , $\mathit{p}=(\mathit{A}-\mathit{D}{)}^{\mathrm{T}}{\mathit{W}}_{\theta }\mathit{b}$ , and define: $\mathit{\lambda }={\left[\begin{array}{cc}{\lambda }_1\hfill & {\lambda }_2\hfill \end{array}\right]}^{\mathrm{T}},\mathit{f}={\left[\begin{array}{cc}{f}_1\hfill & f_2\hfill \end{array}\right]}^{\mathrm{T}}$ . Solve $\lambda$ using the quasi-Newton iterative formula:

$${\lambda }^{(k+1)}={\lambda }^{(k)}-{\mathbf{\alpha }}^{(\mathbf{k})}{\mathit{H}}^{(k)}{\mathit{g}}^{(k)}$$

(20)

In the formula, ${\lambda }^{(k)}$ is the kth iteration of $\lambda$ , ${\mathbf{\alpha }}^{(\mathbf{k})}$ is the step size, ${\mathit{g}}^{(k)}=\nabla f\left({\mathit{\lambda }}^{(k)}\right)$ . ${\mathit{H}}^{(k)}$ is the symmetric matrix of the corresponding dimension, updated by the BFGS formula to:

$$\begin{array}{c}{\mathit{H}}^{(k+1)}={\mathit{H}}^{(k)}-\frac{{\mathit{H}}^{(k)}{\mathit{y}}^{(k)}{\mathit{y}}^{(k)\mathrm{T}}{\mathit{H}}^{(k)}}{{\mathit{y}}^{(k)\mathrm{T}}{\mathit{H}}^{(k)}{\mathit{y}}^{(k)}}+\frac{{\mathit{s}}^{(k)}{\mathit{s}}^{(k)\mathrm{T}}}{{\mathit{y}}^{(k)\mathrm{T}}{\mathit{s}}^{(k)}}\hfill \\ +{\mathit{w}}^{(k)}{\mathit{w}}^{(k)\mathrm{T}}\hfill \end{array}$$

(21)

The specific steps of using the quasi-Newton BFGS formula to iterate are as follows^17–20:

• Set parameter $\delta \in (0,1),\sigma \in (0,0.5),{\mathit{W}}_{\theta }={\mathit{Q}}^{-1}$ . Set the precision $\epsilon =10^{-3}$ , take the initial point ${\lambda }^{(1)}=[0,0{]}^{\mathrm{T}}$ , the initial matrix ${\mathit{H}}^{(1)}$ is the unit matrix I, k = 0;

• Calculate ${\mathbf{g}}^{(k)}=\nabla f\left({\lambda }^k\right)$ , if $\Vert {\mathbf{g}}^{(k)}\Vert \le \epsilon$ , stop the iteration and output ${\lambda }^k$ as the approximate solution;

• Calculate the search direction ${\mathit{p}}^{(k)}=-{\mathit{H}}^{(k)}{\mathit{g}}^{(k)}$ ;

• Let $m_k$ be the smallest non-negative integer that satisfies $f\left({\mathit{\lambda }}^{(k)}\right)+{\delta }^m{\mathit{p}}^{(k)}\le f\left({\mathit{\lambda }}^{(k)}\right)+\sigma {\delta }^m{\mathit{g}}^{(k)T}{\mathit{p}}^{(k)}$ ;

• Use the BFGS formula of the quasi-Newton method to update ${\mathit{H}}^{(k+1)}$ , set k = k + 1, and go to step (b).

Finally, the obtained optimal value of ${\lambda }^k$ is substituted into equation (18) to obtain the optimal estimated value of the vector $\theta$ to be sought, that is, the position information of the target.

Experimental verification and result analysis

In order to verify the feasibility of the proposed ICTLS algorithm to locate multiple dynamic targets, this section designs three simulation experiments for target location scenarios. In order to compare and analyze the positioning performance of the ICTLS algorithm in detail, the simulation experiment in this article compares the simulation results of the three positioning algorithms TSWLS, CTLS, and ICTLS. The experiment is based on the simulation scheme in Ho et al. ⁴ It is assumed that there are six observation sensors in the three-dimensional space. The specific positions and speeds of the sensors are shown in Table 1. In this paper, two targets are taken as examples to carry out the simulation experiment of multi-target positioning, of course, it can also be extended to more targets. The real positions of the two targets to be positioned are ${\mathit{u}}_1^0=[600,650,550{]}^T$ and ${\mathit{u}}_2^0=[310,480,245{]}^{\mathrm{T}}$ respectively, and the speeds of the two targets are ${\stackrel{·}{\mathit{u}}}_1^0=[-20,15,40{]}^{\mathrm{T}}$ and ${\stackrel{·}{\mathit{u}}}_2^0=[40,15,-20{]}^{\mathrm{T}}$ respectively. Suppose the position error of the sensor is a Gaussian random process with zero mean, and the covariance matrix is ${\mathit{Q}}_{\beta }=\mathrm{diag}[{\sigma }_s^2\mathit{R}\left({\mathrm{m}}^2\right),{\stackrel{·}{\sigma }}_s^2\mathit{R}{(\mathrm{m}/\mathrm{s})}^2]$ , ${\sigma }_s^2$ and ${\stackrel{·}{\sigma }}_s^2$ are the position and velocity error variances of the sensor, and ${\stackrel{·}{\sigma }}_s^2=0.1{\sigma }_s^2$ , matrix $\mathit{R}=\mathrm{diag}[1,1,1,2,2,2,10,10,10,40,40,40,20,20,20,3,3,3]$ . The measurement errors based on TDOA and FDOA are all Gaussian random vectors with zero mean, and the covariance matrix is ${\mathit{Q}}_{\alpha }=\mathrm{diag}[{\sigma }_t^2\mathit{J}\left({\mathrm{m}}^2\right),{\sigma }_f^2\mathit{J}{(\mathrm{m}/\mathrm{s})}^2]$ , ${\sigma }_t^2$ and ${\sigma }_f^2$ are the measurement time difference and frequency difference variance respectively, set ${\sigma }_t^2={10}^{-4}(m^2)$ , and J is a (M − 1) × (M + 1) square matrix, whose main diagonal element is 1 and other elements are 0.5. Run 9999 Monte Carlo experiments on a PC to compare the RMSE estimated by the ICTLS algorithm and the TSWLS and CTLS methods for target position and velocity. The RMSE of the estimated target position and velocity is defined as follows:

$$\{\begin{array}{c}\mathrm{P}-\mathrm{RMSE}\left({\mathit{u}}_i\right)=\sqrt{\frac{1}{10000}\sum _{l=1}^{10000}\Vert {\tilde{\mathit{u}}}_i(l)-{\mathit{u}}_i^0{\Vert }_2^2}\hfill \\ \multicolumn{1}{c}{\mathrm{V}-\mathrm{RMSE}\left({\stackrel{·}{\mathit{u}}}_i\right)=\sqrt{\frac{1}{10000}\sum _{l=1}^{10000}\Vert {\widetilde{\stackrel{·}{\mathit{u}}}}_i(l)-{\stackrel{·}{\mathit{u}}}_i^0{\Vert }_2^2}}\end{array}$$

(22)

Table 1.

Sensor position and velocity coordinates.

Sensor number i	$x_i$ (m)	$y_i$ (m)	$z_i$ (m)	${\stackrel{·}{x}}_i$ (m/s)	${\stackrel{·}{y}}_i$ (m/s)	${\stackrel{·}{z}}_i$ (m/s)
1	299	99	151	30	−20	20
2	399	151	101	−30	10	20
3	299	501	199	10	−20	10
4	348	201	101	10	20	30
5	−101	−99	−99	−20	10	10
6	199	−299	−199	20	−10	10

Among them: P-RMSE and V-RMSE respectively represent the RMSE of the estimated position and velocity of the target, and ${\tilde{\mathit{u}}}_i(l)$ is the estimated result of the first Monte Carlo experiment. Three positioning scenarios are designed in the article: ① To locate a single dynamic target, the real position and speed of the target are ${\mathit{u}}_1^0=[600,650,550{]}^{\mathrm{T}}$ and ${\stackrel{·}{\mathit{u}}}_1^0=[-20,15,40{]}^{\mathrm{T}}$ respectively. ② Locate two different dynamic targets. In addition to target ${\mathit{u}}_1^0$ , the real position and velocity of the other target are ${\mathit{u}}_2^0=[310,480,245{]}^{\mathrm{T}}$ and ${\stackrel{·}{\mathit{u}}}_2^0=[40,15,-20{]}^{\mathrm{T}}$ respectively. ③ Locate two different static targets, assume that the speed of the two targets ${\mathit{u}}_1^0$ and ${\mathit{u}}_2^0$ are zero, and select the sensor ${\mathit{s}}_3^0$ close to ${\mathit{u}}_2^0$ as the reference sensor. The location distribution of targets and sensors in space is shown in Figure 1.

Figure 1.

Spatial location distribution of targets and sensors.

Figure 2 shows the positioning performance of the three positioning algorithms TSWLS, CTLS, and ICTLS when locating a single dynamic target, as well as the comparison results of RMSE and CRLB of the target position velocity estimation deviation. It can be seen from the curve in the figure: For target position estimation, when $10\lg \left({\sigma }_s^2\right)\le 5\mathrm{dB}{\mathrm{m}}^2$ , the RMSE of the position estimation deviation of the three positioning algorithms can all reach CRLB. However, the RMSE of the position estimation error of the three positioning algorithms of TSWLS, CTLS, and ICTLS at $10\lg \left({\sigma }_s^2\right)\ge 5{\mathrm{dBm}}^2,10{\mathrm{dBm}}^2,12.5{\mathrm{dBm}}^2$ are gradually higher than CRLB. The reason for the different thresholds of the three positioning algorithms of TSWLS, CTLS, and ICTLS is that the deviation of the WLS results in the first step of the TSWLS positioning algorithm increases with the increase of the noise error, resulting in poor positioning accuracy and weak robustness; Although the CTLS positioning algorithm considers the measurement error and the position error of the sensor as a whole, it does not consider the relationship between the introduced additional variable ${\mathit{r}}_{1,i}^0$ and the target position coordinate ${\mathit{u}}_i^0$ , and does not consider the correlation between the measurement noise and the sensor position error. Therefore, the CTLS positioning algorithm cannot reach the optimal solution; the ICTLS positioning algorithm in the article is based on the CTLS algorithm, not only considering the correlation between the measurement noise and the sensor position error, but also approximating processing of the introduced additional variables reduces the estimation error in the first step, and also consider the correlation between the introduced variables and the target position into the positioning algorithm, and the optimal solution can be reached under constraints.

Figure 2.

Comparison of RMSE and CRLB of the three positioning algorithms under scenario ①.

Therefore, this paper’s ICTLS positioning algorithm has the best performance in positioning dynamic targets among the three positioning algorithms.

In terms of target speed estimation, when $10\lg \left({\sigma }_s^2\right)\le 2.5\mathrm{dB}{\mathrm{m}}^2$ , the RMSE of the speed estimation of the three positioning algorithms of TSWLS, CTLS, and ICTLS can all reach CRLB, and the three positioning algorithms’ thresholds are, respectively, 2.5, 5, and 7.5 dBm². It can be seen that the three positioning algorithms correspond to The position estimation thresholds are respectively greater than their speed estimation thresholds. Therefore, when the three positioning algorithms locate dynamic targets, the ability of speed estimation to adapt to noise is worse than that of position estimation, and the randomness of estimation errors is also increased.

Figure 3 depicts the situation where the RMSE and CRLB of the three positioning algorithms of TSWLS, CTLS, and ICTLS are consistent with the target position and velocity estimation errors when multiple dynamic targets are located. Figure 3 shows the RMSE curve of position and velocity estimation deviation. The upper curve is the RMSE curve of positioning target ${\mathit{u}}_1^0$ , and the lower curve is the RMSE curve of positioning target ${\mathit{u}}_2^0$ . As can be seen from the curve in the figure, when positioning two dynamic targets ${\mathit{u}}_1^0$ and ${\mathit{u}}_2^0$ , the ICTLS positioning algorithm always has higher accuracy and stronger robustness than the other two algorithms in terms of position and velocity estimation, even if the measurement error and sensor position error CRLB can be achieved even larger. In addition, when the TSWLS and CTLS positioning algorithms are obtained from the curve in the figure, when locating a short-range target ${\mathit{u}}_2^0$ , the thresholds in the estimation of target position and speed are smaller than those of the long-range target ${\mathit{u}}_1^0$ . This shows that TSWLS and CTLS The poor performance of positioning algorithms in locating close-range targets is mainly due to the large influence of noise between close-range targets without considering the correlation of noise errors.

Figure 3.

Comparison of RMSE and CRLB of three positioning algorithms under scenario ②.

Figure 4 depicts the RMSE and CRLB comparison results of the position estimation errors of the three positioning algorithms TSWLS, CTLS and ICTLS when positioning two static targets. For long-distance target ${\mathit{u}}_1^0$ , when $10\lg \left({\sigma }_s^2\right)⩽5\mathrm{dB}{\mathrm{m}}^2$ , the RMSE of the three positioning algorithms can reach CRLB. Even in $10\lg \left({\sigma }_s^2\right)>5\mathrm{dB}{\mathrm{m}}^2$ , the RMSE of the three positioning algorithms will increase with the increase of the error noise power, but at the same noise Under the condition of power increase, the increase in RMSE of static target position estimation error is less than the increase in RMSE of dynamic target position estimation error. The reason is the lack of some parameters such as speed and frequency when positioning the static target, which weakens the interaction between the participating variables. Influence, improve the calculation efficiency and positioning accuracy. However, when locating the target ${\mathit{u}}_2^0$ close to the reference sensor, the TSWLS and CTLS positioning algorithms’ threshold is obviously smaller than that of the target ${\mathit{u}}_1^0$ , and the positioning algorithm has poor positioning performance when locating the target close to the reference sensor.

Figure 4.

Comparison of RMSE and CRLB of the three positioning algorithms under scenario ③.

Simulation experiments show that when the measurement error and sensor position error are moderate, the ICTLS positioning algorithm has a smaller RMSE than the TSWLS and CTLS positioning algorithms, showing higher positioning accuracy and stronger robustness. However, when the measurement error and sensor position error are large, the performance of the algorithm in this paper is not as good as expected, which is the limitation of this algorithm.

Math positioning algorithm in specific scenarios

In the previous chapter, based on the estimation of multiple targets’ position by the TSWLS algorithm and the CTLS algorithm, the ICTLS positioning algorithm was proposed. This algorithm fully considered the noise of all coefficient matrices in the positioning equation, introduced the Lagrangian multiplier technology, and used the quasi-Newton The BFGS formula of the method effectively avoids the calculation of the Hessian matrix, reduces the amount of calculation, and comprehensively improves the performance of the positioning algorithm. However, whether it is the classic TSWLS, CTLS positioning algorithm, or the ICTLS positioning algorithm proposed in the article, there is a fatal flaw in locating the target, which leads to a large positioning deviation and a sharp decline in positioning performance. When the target is located at or close to any coordinate axis of the set reference sensor, the above positioning algorithm will have a large positioning deviation, and the positioning performance will drop sharply. In order to avoid the algorithm defects caused by special positioning scenarios, we carefully analyzed the principle of the above algorithm to solve the positioning equation, combined with the measurement error and the position error of the sensor, and introduced a simple and effective strategy based on the TSWLS algorithm, namely, reselecting the reference sensor and the rotating coordinate system avoids the fatal shortcomings of the above positioning algorithm and reduces the target position estimation error in a specific scene. The principle and flow of the improved algorithm to locate the target will be introduced in detail below.

Positioning scene and model establishment

To be better suitable for practical applications, the article considers positioning the target in three-dimensional space. Combining the targets and sensors set in the previous chapter, this section also assumes that M randomly distributed sensors are used to locate a static target ${\mathit{u}}^0={[x^0,y^0,z^0]}^{\mathrm{T}}$ , the real position coordinate of the ith sensor is ${\mathit{s}}_i^0={[x_i^0,y_i^0,z_i^0]}^{\mathrm{T}}$ and i = 1, 2,…, M. The measured position vector of the sensor is $\mathit{s}={[{\mathit{s}}_1^{\mathrm{T}},{\mathit{s}}_2^{\mathrm{T}},\cdots ,{\mathit{s}}_M^{\mathrm{T}}]}^{\mathrm{T}}={\mathit{s}}^0+\mathrm{\Delta }\mathit{s}$ , and the real position vector is ${\mathit{s}}^0={[{\mathit{s}}_1^{0\mathrm{T}},{\mathit{s}}_2^{0\mathrm{T}},\cdots ,{\mathit{s}}_M^{0\mathrm{T}}]}^{\mathrm{T}}$ , and $\mathrm{\Delta }\mathit{s}={[\mathrm{\Delta }{\mathit{s}}_1^{\mathrm{T}},\mathrm{\Delta }{\mathit{s}}_2^{\mathrm{T}},\cdots ,\mathrm{\Delta }{\mathit{s}}_M^{\mathrm{T}}]}^{\mathrm{T}}$ is the sensor position error vector. Suppose $\mathrm{\Delta }\mathbf{s}$ is a Gaussian random process with zero mean and the covariance matrix is ${\mathit{Q}}_{\mathrm{s}}$ .

To be more representative, let be the reference sensor and be the origin of the coordinate system. Refer to the positioning scenario and model establishment in Chapter 2, the measurement equation based on TDOA is:

$$\mathit{r}={\mathit{r}}^0+\mathrm{\Delta }\mathit{r}$$

(23)

Where $\mathit{r}={[r_{2,k},r_{3,k},\cdots ,r_{M,k}]}^{\mathrm{T}}$ , $\mathrm{\Delta }\mathit{r}={[\mathrm{\Delta }{r}_{2,k},\mathrm{\Delta }{r}_{3,k},\cdots ,\mathrm{\Delta }{r}_{M,k}]}^{\mathrm{T}}$ is the Gaussian random measurement error vector with zero mean and the covariance matrix is ${\mathit{Q}}_{\mathrm{t}}$ . The additional variable $r_k^0$ is also introduced, and $r_k^0$ is determined by the target position coordinate ${\mathit{u}}^0$ and the reference sensor ${\mathit{s}}_k^0$ . Therefore, the additional variable $r_k^0$ is expanded by a first-order Taylor series at ${\mathit{s}}_k^0$ , and only the linear term of $\mathrm{\Delta }{\mathit{s}}_k$ is retained. After sorting, we get:

$$r_k^0\approx \Vert {\mathit{u}}^0-{\mathit{s}}_k\Vert +{\mathit{\rho }}_{{\mathit{u}}^0,{\mathit{s}}_k}^{\mathrm{T}}\mathrm{\Delta }{\mathit{s}}_k={\hat{r}}_k^0+{\mathit{\rho }}_{{\mathit{u}}^0,{\mathit{s}}_k}^{\mathrm{T}}\mathrm{\Delta }{\mathit{s}}_k$$

(24)

In the formula, ${\hat{r}}_k^0=\Vert {\mathit{u}}^0-{\mathit{s}}_k\Vert$ is the estimated value of and $r_k^0$ is the unit vector (from ${\mathit{s}}_k$ to ${\mathit{u}}^0$ ). Substituting equation (23) into equation (24) and finishing the pseudo-linear positioning equation for TDOA:

$$\begin{array}{c}{r}_i^0{n}_{i,k}+{\mathit{d}}_i^{\mathrm{T}}\mathrm{\Delta }{\mathit{s}}_k+{({\mathit{u}}^0-{\mathit{s}}_i)}^T\mathrm{\Delta }{\mathit{s}}_i\hfill \\ =0.5(r_{i,k}^2+{\mathit{s}}_k^{\mathrm{T}}{\mathit{s}}_k-{\mathit{s}}_i^{\mathrm{T}}{\mathit{s}}_i)-{({\mathit{s}}_k-{\mathit{s}}_i)}^{\mathrm{T}}{\mathit{u}}^0+r_{i,k}{\hat{r}}_k^0\hfill \end{array}$$

(25)

Algorithm optimization

Let ${\mathit{\phi }}_1^0={[{\mathit{u}}^{0\mathrm{T}},{\hat{r}}_k^0]}^{\mathrm{T}}$ , synthesize all the pseudo-linear equations of formula (25) to obtain the positioning equation in matrix form of the vector to be sought:

$${\mathit{\epsilon }}_1={\mathit{h}}_1-{\mathit{G}}_1{\phi }_1^0={\mathit{B}}_1\mathit{n}+{\mathit{D}}_1\mathrm{\Delta }\mathit{s}$$

(26)

Where, $\mathit{h}=\frac{1}{2}\left[\begin{array}{c}{r}_{2,k}^2+{\mathit{s}}_k^{\mathrm{T}}{\mathit{s}}_k-{\mathit{s}}_2^{\mathrm{T}}{\mathit{s}}_2\hfill \\ \multicolumn{1}{c}{r_{3,k}^2+{\mathit{s}}_k^{\mathrm{T}}{\mathit{s}}_k-{\mathit{s}}_3^{\mathrm{T}}{\mathit{s}}_3}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{r_{M,k}^2+{\mathit{s}}_k^{\mathrm{T}}{\mathit{s}}_k-{\mathit{s}}_M^{\mathrm{T}}{\mathit{s}}_M}\end{array}\right],\mathit{G}=\left[\begin{array}{cc}{({\mathit{s}}_k-{\mathit{s}}_2)}^{\mathrm{T}}\hfill & -r_{2,k}\hfill \\ \multicolumn{1}{c}{{({\mathit{s}}_k-{\mathit{s}}_3)}^{\mathrm{T}}}& -r_{3,k}\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill \\ \multicolumn{1}{c}{{({\mathit{s}}_k-{\mathit{s}}_M)}^{\mathrm{T}}}& -r_{M,k}\hfill \end{array}\right],$ ${\mathit{B}}_1=\mathrm{diag}\left\{{[r_2^0,r_3^0,\cdots ,r_M^0]}^{\mathrm{T}}\right\}$ , ${\mathit{D}}_1=\left[\begin{array}{cccc}{\mathit{d}}_2^{\mathrm{T}}\hfill & {\mathit{0}}_{3(2-1)\times 1}^{\mathrm{T}}\hfill & {({\mathit{u}}^0-{\mathit{s}}_2)}^{\mathrm{T}}\hfill & {\mathit{0}}_{3(M-2)\times 1}^{\mathrm{T}}\hfill \\ \multicolumn{1}{c}{{\mathit{d}}_3^{\mathrm{T}}}& {\mathit{0}}_{3(3-1)\times 1}^{\mathrm{T}}\hfill & {({\mathit{u}}^0-{\mathit{s}}_3)}^{\mathrm{T}}\hfill & {\mathit{0}}_{3(M-3)\times 1}^{\mathrm{T}}\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & ⋮\hfill & ⋮\hfill \\ \multicolumn{1}{c}{{\mathit{d}}_M^{\mathrm{T}}}& {\theta }_{3(M-1)\times 1}^{\mathrm{T}}\hfill & {({\mathit{u}}^0-{\mathit{s}}_M)}^{\mathrm{T}}\hfill & {\mathit{0}}_{3(M-M+1)\times 1}^{\mathrm{T}}\hfill \end{array}\right]$ .

Therefore, equation (26) becomes a linear equation about ${\phi }_1^0$ , and the WLS estimation result of ${\phi }_1^0$ is:

$${\phi }_1={[{\hat{\mathit{u}}}^{\mathrm{T}},{\hat{r}}_k]}^{\mathrm{T}}={\left({\mathit{G}}_1^{\mathrm{T}}{\mathit{W}}_1{\mathit{G}}_1\right)}^{-1}{\mathit{G}}_1^{\mathrm{T}}{\mathit{W}}_1{\mathit{h}}_1$$

(27)

Where $\hat{u}$ and ${\hat{r}}_k$ are the estimated values of $\mathit{u}$ and $r_k$ respectively. The weighting matrix ${\mathit{W}}_1$ is

$${\mathit{W}}_1=\mathrm{E}{\left({\mathit{\epsilon }}_1{\mathit{\epsilon }}_1^{\mathrm{T}}\right)}^{-1}={({\mathit{B}}_1{\mathit{Q}}_t{\mathit{B}}_1^{\mathrm{T}}+{\mathit{D}}_1{\mathit{Q}}_s{\mathit{D}}_1^{\mathrm{T}})}^{-1}$$

(28)

Define $\mathrm{\Delta }{\phi }_1={\phi }_1-{\phi }_1^0$ as the estimation error of the first WLS, then the estimation error $\mathrm{\Delta }{\phi }_1$ can be expressed as

$$\mathrm{\Delta }{\phi }_1={[\mathrm{\Delta }{\mathit{u}}^{\mathrm{T}},\mathrm{\Delta }{\hat{r}}_k]}^{\mathrm{T}}={\left({\mathit{G}}_1^{\mathrm{T}}{\mathit{W}}_1{\mathit{G}}_1\right)}^{-1}{\mathit{G}}_1^{\mathrm{T}}{\mathit{W}}_1{\mathit{\epsilon }}_1$$

(29)

So the covariance matrix of the estimated vector ${\phi }_1$ is

$$\mathrm{cov}\left({\phi }_1\right)={\left({\mathit{G}}_1^{\mathrm{T}}{\mathit{W}}_1{\mathit{G}}_1\right)}^{-1}$$

(30)

When the estimated deviation of vector ${\phi }_1$ is relatively small, ${\phi }_1$ can be regarded as a random vector fluctuating around ${\phi }_1^0$ . Therefore, the estimated vector ${\phi }_1$ can be expressed as

$$\{\begin{array}{cc}{\phi }_1(1)=x^0+e_1\hfill & {\phi }_1(2)=y^0+e_2\hfill \\ \multicolumn{1}{c}{{\phi }_1(3)=z^0+e_3}& {\phi }_1(4)=r_k^0+e_4\hfill \end{array}$$

(31)

In the formula, ${\phi }_1(i),i=1,2,3,4$ is the ith element of the estimated vector ${\phi }_1$ , and $e_i$ is the estimated deviation of the ith element in ${\phi }_1$ . The first three elements in ${\phi }_1$ are respectively subtracted from the position coordinate ( $x_k$ , $y_k$ , and $z_k$ ) of the reference sensor, and the two sides are squared to obtain another set of linear positioning equations:

$${\mathit{\epsilon }}_2={\mathit{h}}_2-{\mathit{G}}_2{\phi }_2$$

(32)

Where, ${\mathit{h}}_2=\left[\begin{array}{c}{({\phi }_1(1)-x_k)}^2\hfill \\ \multicolumn{1}{c}{{({\phi }_1(2)-y_k)}^2}\\ \multicolumn{1}{c}{{({\phi }_1(3)-z_k)}^2}\\ \multicolumn{1}{c}{{({\phi }_1(4))}^2}\end{array}\right],{\mathit{\phi }}_2=\left[\begin{array}{c}{(x^0-x_k)}^2\hfill \\ \multicolumn{1}{c}{{(y^0-y_k)}^2}\\ \multicolumn{1}{c}{{(z^0-z_k)}^2}\end{array}\right],$ ${\mathit{G}}_2=\left[\begin{array}{c}\mathrm{diag}(1,1,1)\hfill \\ \multicolumn{1}{c}{{\mathit{1}}_{1\times 3}}\end{array}\right]$ , ${\epsilon }_2$ is the error vector of the first WLS estimation B. Substituting formula (31) into formula (32), after sorting, we get:

$$\{\begin{array}{c}{\epsilon }_2(1)=2(x^0-x_k)e_1+e_1^2\approx 2(x^0-x_k)e_1\hfill \\ \multicolumn{1}{c}{{\epsilon }_2(2)=2(y^0-y_k)e_2+e_2^2\approx 2(y^0-y_k)e}\\ \multicolumn{1}{c}{{\epsilon }_2(3)=2(z^0-z_k)e_3+e_3^2\approx 2(z^0-z_k)e_3}\\ \multicolumn{1}{c}{{\epsilon }_2(4)=2r_k^0{e}_4+e_4^2\approx 2r_k^0{e}_4}\end{array}$$

(33)

Where ${\epsilon }_2(i),i=1,2,3,4$ is the ith element of ${\mathit{\epsilon }}_2$ . When the target ${\mathit{u}}^0$ is located at or close to the reference sensor ${\mathit{s}}_k$ , $(x^0-x_k)\to 0,(y^0-y_k)\to 0\mathrm{or}(z^0-z_k)\to 0$ appears, then the approximate error equation (33) will become meaningless. Both the classic TSWLS and CTLS positioning algorithms, or the ICTLS mentioned in the article, have this hidden danger. Therefore, in order to avoid the impact of this malpractice on target positioning, we need to re-select the reference sensor, usually based on the first estimation result ${\phi }_1$ , select the sensor ${\mathit{s}}_{k^{\prime }}$ farthest from the target as the new reference sensor.^21,22

In order to further analyze the reason why the re-selectors refer to the sensor, define ${\phi }_{{\mathit{u}}^0}=\left[\begin{array}{c}2(x^0-x_k){\mathit{e}}_1+{\mathit{e}}_1^2\hfill \\ \multicolumn{1}{c}{2(y^0-y_k){\mathit{e}}_2+{\mathit{e}}_2^2}\\ \multicolumn{1}{c}{2(z^0-z_k){\mathit{e}}_3+{\mathit{e}}_3^2}\\ \multicolumn{1}{c}{2r_k^0{\mathit{e}}_4+{\mathit{e}}_4^2}\end{array}\right]$ and the approximate error of equation (33) can be rewritten as:

$$\begin{array}{c}\sum _{i=1}^4\left|\frac{{\mathit{\phi }}_u(i)-{\mathit{\epsilon }}_2(i)}{{\mathit{\epsilon }}_2(i)}\right|\hfill \\ =\frac{1}{2}(\left|\frac{{\mathit{e}}_1}{x^0-x_k}\right|+\left|\frac{{\mathit{e}}_2}{y^0-y_k}\right|+\left|\frac{{\mathit{e}}_3}{z^0-z_k}\right|+\left|\frac{{\mathit{e}}_4}{r_k^0}\right|)\hfill \end{array}$$

(34)

In order to facilitate the calculation of the coordinates, the position distribution of the target and the reference sensor is shown in Figure 5. From the geometric relationship in the figure, $x^0-x_k=r_k\sin \theta \cos \varphi ,y^0-y_k=r_k\sin \theta \sin \varphi$ , and $z=r_k\cos \theta$ can be obtained. Substituting the above relationship into equation (34), we get:

$$\begin{array}{c}\sum _{i=1}^4\left|\frac{{\mathit{\phi }}_u(i)-{\mathit{\epsilon }}_2(i)}{{\mathit{\epsilon }}_2(i)}\right|\hfill \\ =\frac{1}{2r_k}[\left|\frac{1}{\sin \theta }\right|(\left|\frac{{\mathit{e}}_1}{\cos \varphi }\right|+\left|\frac{{\mathit{e}}_2}{\sin \varphi }\right|)+\left|\frac{{\mathit{e}}_3}{\cos \theta }\right|+\left|{\mathit{e}}_4\right|]\hfill \end{array}$$

(35)

Figure 5.

Target and reference sensor distribution.

The mathematical expectation of the approximate error formula (35) is:

$$E\left(\sum _{i=1}^4\left|\frac{{\mathit{\phi }}_u(i)-{\mathit{\epsilon }}_2(i)}{{\mathit{\epsilon }}_2(i)}\right|\right)=\frac{P}{2r_k}$$

(36)

Where $P=\left|\frac{1}{\sin \theta }\right|(\frac{E\left(\left|e_1\right|\right)}{|\cos \varphi |}+\frac{E\left(\left|e_2\right|\right)}{|\sin \varphi |})+\frac{E\left(\left|e_3\right|\right)}{|\cos \theta |}+E\left|e_4\right|$ , $r_k$ represents the measured distance between the kth sensor and the target.

In order to minimize the mathematical expectation formula (36) of the estimated deviation, two steps are required: ① let $r_k$ take the maximum value, ② let $P$ take the minimum value. ²³ In order to get the maximum value of $r_k$ , a new sensor ${\mathit{s}}_{k^{\prime }}$ needs to be selected as the reference sensor, and ${\mathit{s}}_{k^{\prime }}$ is the sensor farthest from the target among all sensors, so as to ensure that ${\mathit{s}}_{k^{\prime }}$ gets the maximum value. In order to get the minimum value of P, we carefully analyzed the expression of P. In view of the relatively small estimation error, $E\left(\left|{\mathit{e}}_1\right|\right)=E\left(\left|{\mathit{e}}_2\right|\right)=E\left(\left|{\mathit{e}}_3\right|\right)$ can be obtained, and P can get the minimum value when $\theta \approx {55}^{°},\varphi ={45}^{°}$ . So we can make $\theta \approx {45}^{°},\varphi ={55}^{°}$ by rotating the coordinate axis twice. ²⁴

Rotate the coordinate axis for the first time so that $\varphi ={45}^{°}$ , take the z axis as the rotation axis, rotate the coordinate axis clockwise, and the rotation angle is $\mathbf{\alpha }=\mathbf{\pi }/4$ . Therefore, after rotating the coordinate axis for the first time, the ith sensor in the new coordinate system can be expressed as:

$${\mathit{s}}_i^{\prime }=\left[\begin{array}{ccc}\cos \alpha \hfill & \sin \alpha \hfill & 0\hfill \\ \multicolumn{1}{c}{-\sin \alpha }& \cos \alpha \hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill \end{array}\right]{\mathit{s}}_i$$

(37)

Rotate the coordinate axis for the second time so that $\theta \approx {55}^{°}$ , take the line passing through the origin of the coordinate and perpendicular to the plane $\overline{s_k\mathrm{ug}}$ as the rotation axis, rotate the coordinate axis clockwise, the rotation angle is $\mathit{\beta }=55\mathbf{\pi }/180$ , and the cosine direction of the rotation axis is $n_x=n_y=\sqrt{2}/2,n_z=0$ . ²⁵ Therefore, after the second rotation of the coordinate axis, the coordinate system of the second rotation of the ith sensor can be expressed as

$${\mathit{s}}_i^{\prime \prime }=\left[\begin{array}{ccc}\cos \beta \hfill & \frac{1-\cos \beta }{2}\hfill & \frac{\sqrt{2}\sin \beta }{2}\hfill \\ \multicolumn{1}{c}{\frac{1-\cos \beta }{2}}& \frac{1-\cos \beta }{2}+\cos \beta \hfill & \frac{\sqrt{2}\sin \beta }{2}\hfill \\ \multicolumn{1}{c}{\frac{\sqrt{2}\sin \beta }{2}}& \frac{\sqrt{2}\sin \beta }{2}\hfill & \cos \beta \hfill \end{array}\right]{\mathit{s}}_i^{\prime }$$

(38)

In the same way, the target ${\mathit{u}}^0$ becomes $({\mathit{u}}^0{)}^{\prime \prime }$ after two rotations of the coordinate axis.

The WLS estimation result of the program (32) under the new coordinate system is

$${\mathit{\phi }}_2^{\prime \prime }={\left({\mathit{G}}_2^{\mathrm{T}}{\mathit{W}}_2{\mathit{G}}_2\right)}^{-1}{\mathit{G}}_2^{\mathrm{T}}{\mathit{W}}_2{\mathit{h}}_2^{\prime \prime }$$

(39)

Where ${\mathit{W}}_2$ is the weighting matrix:

$$\{\begin{array}{c}{\mathit{W}}_2=\mathrm{E}{\left[{\mathit{\epsilon }}_{\mathit{2}}{\mathit{\epsilon }}_2^{\mathrm{T}}\right]}^{-1}={\left[{\mathit{B}}_2\mathrm{cov}\left({\mathit{\phi }}_1^{\prime \prime }\right){\mathit{B}}_2^{\mathrm{T}}\right]}^{-1}\hfill \\ \multicolumn{1}{c}{{\mathit{B}}_2=\mathrm{diag}\{{\left(x^0\right)}^{\prime \prime }-x_{k^{\prime }}^{\prime \prime },{\left(y^0\right)}^{\prime \prime }-y_{k^{\prime }}^{\prime \prime }{\left(z^0\right)}^{\prime \prime }-z_{k^{\prime }}^{\prime \prime },{\left(r_{k^{\prime }}^0\right)}^{\prime \prime }\}}\end{array}$$

(40)

The covariance matrix of the estimated vector ${\phi }^{\prime \prime }$ is approximate:

$$\mathrm{cov}\left({\mathit{\phi }}_2^{\prime \prime }\right)={\left({\mathit{G}}_2^{\mathrm{T}}{\mathit{W}}_2{\mathit{G}}_2\right)}^{-1}$$

(41)

Therefore, the estimated position of the target under the new coordinate system is:

$${\hat{\mathit{u}}}^{\prime \prime }=\mathbf{\Pi }\sqrt{{\mathit{\epsilon }}_2}+{\left[\begin{array}{ccc}{\mathbf{x}}_{{\mathbf{k}}^{\prime }}^{\prime \prime }\hfill & {\mathbf{y}}_{{\mathbf{k}}^{\prime }}^{\prime \prime }\hfill & {\mathbf{z}}_{{\mathbf{k}}^{\prime }}^{\prime \prime }\hfill \end{array}\right]}^{\mathrm{T}}$$

(42)

In the formula, $\mathbf{\Pi }=\mathrm{diag}\left\{\mathrm{sgn}({\mathit{u}}^{\prime \prime }(1:3)-{\mathit{s}}_{{\mathbf{k}}^{\prime }}^{\prime \prime })\right\}$ , and $\mathrm{sgn}(*)$ is a symbolic function. The position estimate $\hat{u}$ of the original coordinate system’s target can be obtained by formulas (37) and (38). The improved TSWLS localization algorithm is shown in Algorithm 1.

Algorithm 1: Improved TSWLS localization algorithm

Input: Number of sensors M, sensor measurement vector $\mathit{s}$ , error vector $\mathrm{\Delta }\mathit{s}$ . Output: Estimated value of location $\hat{\mathit{u}}$ . 1: Introduce additional variable $r_k^0$ and expand it at ${\mathit{s}}_k$ according to (24); 2. According to (23) and (25), a pseudo-linear positioning equation based on TDOA is established; 3. Set ${\phi }_1^0={[{\mathit{u}}^{0T},{\hat{r}}_k^0]}^{\mathrm{T}}$ , through (26)–(29), use WLS to make a preliminary estimate of vector ${\phi }_1$ ; 4. Rotate the coordinate axis twice to make $\varphi ={45}^{°}$ , $\theta \approx {55}^{°}$ , get a new coordinate system, and calculate the coordinates ${\mathit{s}}_i^{\prime }$ and ${\mathit{s}}_i^{\prime \prime }$ of the ith sensor in the new coordinate system; 5. According to (39)–(40), the error vector ${\phi }_2$ is estimated by WLS in the new coordinate system; 6. Calculate the estimated target position ${\hat{\mathit{u}}}^{\prime \prime }$ under the new coordinate system according to (42); 7. The estimated position value $\hat{\mathit{u}}$ of the target in the original coordinate system can be obtained by (37) and (38).

Algorithm performance analysis

In order to facilitate the calculation, we analyze the performance of the optimized algorithm in the paper under the original coordinate system. The following will evaluate the positioning performance of the optimized algorithm by comparing $\mathrm{cov}\left({\mathit{\phi }}_2^{\prime \prime }\right)$ in equation (41) with the CRLB of the position estimation. The detailed derivation process of CRLB ( ${\mathit{u}}^0$ ) for target position estimation can be found in literature, ⁷ that is, $\mathrm{CRLB}{\left({\mathit{u}}^0\right)}^{-1}={\mathit{X}}^{-1}-\mathit{Y}{\mathit{Z}}^{-1}{\mathit{Y}}^{\mathrm{T}}$ , and it satisfies:

$$\{\begin{array}{c}X={\left(\frac{\partial {\mathit{r}}^0}{\partial {\mathit{u}}^0}\right)}^{\mathrm{T}}{\mathit{Q}}_t^{-1}\left(\frac{\partial {\mathit{r}}^0}{\partial {\mathit{u}}^0}\right)\hfill \\ \multicolumn{1}{c}{Y={\left(\frac{\partial {\mathit{r}}^0}{\partial {\mathit{u}}^0}\right)}^{\mathrm{T}}{\mathit{Q}}_t^{-1}\left(\frac{\partial {\mathit{r}}^{\prime }}{\partial {\mathit{s}}^0}\right)}\\ \multicolumn{1}{c}{Z={\left(\frac{\partial {\mathit{r}}^0}{\partial {\mathit{s}}^0}\right)}^{\mathrm{T}}{\mathit{Q}}_t^{-1}\left(\frac{\partial {\mathit{r}}^0}{\partial {\mathit{s}}^0}\right)+{\mathit{Q}}_s^{-1}}\end{array}$$

(43)

Substituting ${\mathit{W}}_1$ in equation (28), $\mathrm{cov}\left({\phi }_1\right)$ in equation (30), and ${\mathit{W}}_2$ in equation (40) into equation (41), we get:

$$\begin{array}{c}\mathrm{cov}{\left({\phi }_2\right)}^{-1}=\hfill \\ G_3^{\mathrm{T}}{Q}_t^{-1}{G}_3-G_3^{\mathrm{T}}{Q}_t^{-1}{G}_4{(Q_s^{-1}+G_4^{\mathrm{T}}{Q}_t^{-1}{G}_4)}^{-1}{G}_4^{\mathrm{T}}{Q}_t^{-1}{G}_3\hfill \end{array}$$

(44)

In the formula, ${\mathit{G}}_3={\mathit{B}}_1^{-1}{\mathit{G}}_1{\mathit{B}}_2^{-1}{\mathit{G}}_2,{\mathit{G}}_4={\mathit{B}}_1^{-1}{\mathit{D}}_1$ . Comparing $\mathrm{cov}{\left({\phi }_2\right)}^{-1}$ and CRLB ( ${\mathit{u}}^0$ )⁻¹, it is not difficult to find that the two have the same expression form. When the position error of the sensor is small, then ${\mathit{G}}_4\approx (\partial {\mathit{r}}^0/\partial {\mathit{s}}^0)$ is established. It is only necessary to prove that ${\mathit{G}}_3\approx (\partial {\mathit{r}}^0/\partial {\mathit{u}}^0)$ is established,then $\mathrm{cov}(\mathrm{\Delta }\mathit{u}{)}^{-1}\approx \mathrm{CRLB}{\left({\mathit{u}}^0\right)}^{-1}$ is established. The certification process is as follows:

Let the ith row element of $(\partial {\mathit{r}}^0/\partial {\mathit{u}}^0)$ be:

$$\frac{\partial r_{i,k}^0}{{(\partial {\mathit{u}}^0)}^{\mathrm{T}}}={\mathit{\rho }}_{u^0,s_i^0}^{\mathrm{T}}-{\mathit{\rho }}_{u^0,s_k^0}^{\mathrm{T}}$$

(45)

In the same way, if $B_1,G_1,B_2\mathrm{and}{G}_2$ are brought into $G_3$ , then the ith row element of $G_3$ is

$$G_3(i:)=\frac{{({\mathit{s}}_i-{\mathit{s}}_k)}^{\mathrm{T}}}{r_i^0}+\frac{r_{i,k}}{r_i^0}{\mathit{\rho }}_{u,s_k}^{\mathrm{T}}$$

(46)

When the measurement error and the sensor position error are small, $\mathrm{\Delta }{s}_i/r_i^0\approx 0$ , $\mathrm{\Delta }{s}_k/r_i^0\approx 0$ , $n_{i,k}/r_i^0\approx 0$ , and ${\phi }_2/r_1^0\approx 0$ is established. Therefore, the following approximate formula is obtained:

$$\frac{r_{i,k}}{r_i^0}=\frac{r_i^0-r_k^0}{r_i^0}+\frac{n_{i,k}}{r_i^0}\approx \frac{r_i^0-r_k^0}{r_i^0}$$

(47)

$$\begin{array}{c}\multicolumn{1}{c}{{\mathit{\rho }}_{u,s_k}^{\mathrm{T}}=\frac{{({\mathit{u}}^0+{\phi }_2-{\mathit{s}}_k)}^{\mathrm{T}}}{\Vert {\mathit{u}}^0+{\phi }_2-{\mathit{s}}_k\Vert }\approx \frac{{({\mathit{u}}^0-{\mathit{s}}_k)}^{\mathrm{T}}}{\Vert {\mathit{u}}^0-{\mathit{s}}_k\Vert }}\\ +\frac{{\phi }_2^{\mathrm{T}}}{\Vert {\mathit{u}}^0-{\mathit{s}}_k\Vert }\approx {\mathit{\rho }}_{u^0,s_k^0}^{\mathrm{T}}\hfill \end{array}$$

(48)

Substituting formula (47) and (48) into formula (46), we get:

$$\begin{array}{c}{G}_3(i:)\approx \frac{{({\mathit{u}}^0-{\mathit{s}}_k)}^{\mathrm{T}}-{({\mathit{u}}^0-{\mathit{s}}_i)}^{\mathrm{T}}}{r_i^0}+(1+\frac{r_k^0}{r_i^0}){\mathit{\rho }}_{u^0,s_k^0}^{\mathrm{T}}\hfill \\ \approx \frac{{({\mathit{u}}^0-{\mathit{s}}_k^0)}^{\mathrm{T}}}{r_k^0}-\frac{{({\mathit{u}}^0-{\mathit{s}}_i^0)}^{\mathrm{T}}}{r_i^0}\approx -\frac{\partial r_{i,k}^0}{\partial {\left({\mathit{u}}^0\right)}^{\mathrm{T}}}\hfill \end{array}$$

(49)

Therefore, $\mathrm{cov}(\mathrm{\Delta }\mathit{u}{)}^{-1}\approx \mathrm{CRLB}{\left({\mathit{u}}^0\right)}^{-1}$ is established, indicating that the position estimation accuracy of the proposed optimization algorithm can reach CRLB.

In order to comprehensively evaluate the performance of the optimized algorithm in the paper, we will compare the computational complexity of the optimized algorithm in the paper with the TSWLS algorithm in Sun and Ho. ⁵ Generally, RM(X) is used to represent the number of multiplications required to calculate the variable X. Usually, N³RMs are required to inverse a N × N matrix, and M × N× PRMs is required to multiply the matrix of M × N and the matrix of N×P. Taking locating a short-range target as an example, the computational complexity of the TSWLS algorithm in Sun and Ho ⁵ and the optimization algorithm in the paper are shown in Tables 2 and 3, respectively.

Table 2.

The computational complexity of the TSWLS algorithm in Sun and Ho. ⁵

Step	Equations in the literature	RMs
First step
Initialization ${\mathit{W}}_1={{\mathit{Q}}_{\alpha }}^{-1}$	Formula (37)	0
Calculate ${\phi }_1$	Formula (17)	4M2 + 28M2 + 32
Calculate ${\mathit{W}}_1$	Formula (18)	2 × (5M3 − 15M2 + 15M − 5)
Calculate ${\phi }_1$	Formula (17)	2 × (4M2 + 28M2 + 32)
Second step
Calculate $\mathrm{cov}\left({\phi }_1\right)$	Formula (19)	0
Form ${\mathit{W}}_2$	Formula (32)	192
Calculate ${\phi }_2$	Formula (31)	320
Calculate $\hat{\mathit{u}}$	Formula (34)	3
Calculate B 2	Formula (28)	0
Calculate ${\mathit{W}}_2$	Formula (32)	2 × 192
Calculate ${\phi }_2$	Formula (31)	2 × 320
Calculate $\hat{\mathit{u}}$	Formula (34)	2 × 3
Total number of multiplications: $10M^3-18M^2+114M+1631\mathrm{RMs}$

Table 3.

The computational complexity of the optimized algorithm in this paper.

Step	The formula in this paper	RMs
First step		$4M^2+28M+32$
Repeat twice	Formulas (7) and (8)	$10M^3-22M^2+86M+54$
Second step
Calculate ${\mathit{s}}_i^{\prime }$ and ${\mathit{s}}_i^{\prime \prime }$	Formulas (11) and (12)	9 + 9
Use ${\phi }_1$ to calculate ${\phi }_2$		18
Third step		415
Repeat twice		2 × 415
Use ${\hat{\mathit{u}}}^{\prime \prime }$ to calculate $\hat{\mathit{u}}$		18
Total number of multiplications: $10M^3-18M^2+114M+1331\mathrm{RMs}$

From Tables 2 and 3, it can be seen that the optimization algorithm in the paper requires less total computation than the TSWLS algorithm in literature, ⁵ and the optimization algorithm in this paper can show better positioning performance, especially when the target is at or close to the reference sensor.

Experimental verification and result analysis

This section compares the proposed optimization algorithm’s performance and the TSWLS algorithm for positioning targets through simulation experiments in two positioning scenarios. The two positioning scenes are ① special positioning scene, that is, the target is located at or close to the reference sensor; ② random positioning scene, that is, the target and sensors are randomly distributed in three-dimensional space. The experiments in the two positioning scenarios were run 5000 Monte Carlo experiments on a PC, and the positioning performance of the algorithm was analyzed by comparing the MSE and deviation of the position estimation of the two algorithms. The MSE and bias of the position estimate are defined as ²⁶ :

$$\{\begin{array}{c}{R}_{\mathrm{MSE}}(\hat{\mathit{u}})=\sum _{l=1}^L\frac{{\Vert {\hat{\mathit{u}}}_l-{\mathit{u}}^0\Vert }^2}{L}\hfill \\ \multicolumn{1}{c}{R_{\mathrm{BIAS}}(\hat{\mathit{u}})={\Vert \sum _{l=1}^L\frac{\mid {\hat{\mathit{u}}}_l-{\mathit{u}}^0\parallel }{L}\Vert }^2}\end{array}$$

(50)

Where ${\hat{\mathit{u}}}_l$ is the estimated value of ${\mathit{u}}_l$ ’s lth Monte Carlo experiment and L = 5000. The position error variance of the sensor is ${\sigma }_s^2$ , and the covariance matrix is ${\mathit{Q}}_s={\sigma }_s^2·\mathrm{diag}[10,10,10,2,2,2,10,10,10,40,40,40,20,20,20,3,3,3]$ ; the variance of the average measurement error is ${\sigma }_t^2={10}^{-4}{m}^2$ , its covariance matrix is ${\mathit{Q}}_t={\sigma }_t^2\mathit{R}$ , and $\mathit{R}$ is a main diagonal element of 1 and other elements of 0.5 (M − 1) × (M − 1) matrix; assume that the measurement error and the position error of the sensor are independent of each other. ²¹

Target positioning in a specific scene

There are two targets to be located in the monitoring area, one is the short-range target ${\mathit{u}}_1^0=[600,80,300]\mathrm{m}$ , and the other is the long-range target ${\mathit{u}}_2^0=[2000,3000,1500]\mathrm{m}$ . The two targets to be located are used to estimate the target position information from the data measured by the six sensors in the area. The real position coordinates of the sensors are shown in Table 4.

Table 4.

Real position coordinates of the sensor.

Sensor number i	xi/m	yi/m	zi/m
1	299	99	151
2	399	151	101
3	299	501	199
4	348	201	101
5	−101	−99	−99
6	199	−299	−199

The distribution of targets and sensors in space is shown in Figure 6, assuming that A is the reference sensor. We compare the position estimation deviation and MSE to evaluate the positioning performance of the TSWLS algorithm in Sun and Ho ⁵ and the optimized algorithm in the paper.

Figure 6.

Target and sensor position distribution in a specific scenario.

Figure 7 depicts the relationship between the MSE and the deviation of the short-distance target position estimated by the two algorithms and the sensor position error variance. As can be seen from the figure, $10\lg \left({\sigma }_s^2\right)\le -40\mathrm{dB}{\mathrm{m}}^2$ ; and the threshold of the optimized algorithm in the text is −10 dBm², which is 30 dBm² higher than the threshold of the TSWLS algorithm, which has strong robustness. In addition, the optimization algorithm in the article at $10\lg \left({\sigma }_s^2\right)=-20\mathrm{dB}{\mathrm{m}}^2$ time reduces the estimation deviation by about 8 dB compared with the TSWLS algorithm. Therefore, the optimization algorithm shows high positioning accuracy and strong robustness when locating short-range targets.

Figure 7.

Performance comparison of two algorithms when locating close-range targets.

Figure 8 shows the relationship between the MSE and the deviation of the long-distance target position estimated by the two algorithms and the sensor position error variance. It can be seen from the curve in the figure that when $10\lg \left({\sigma }_s^2\right)\le -10\mathrm{dB}{\mathrm{m}}^2$ , the MSE estimated by the two positioning algorithms can reach CRLB. The noise thresholds of the TSWLS algorithm and the optimized algorithm are −10 and −5 dBm², respectively. Obviously, the optimized algorithm is more robust than the TSWLS algorithm. In addition, the optimization algorithm also has a smaller position estimation deviation than the TSWLS algorithm. In short, when locating long-distance targets, the optimization algorithm also showed better positioning performance.

Figure 8.

Performance comparison of two algorithms when locating long-distance targets.

Figure axis labels target positioning in random acenes

Hundred target sources are randomly distributed in the setting space, and the target sources are randomly distributed in a cube with a side length of 600 m. Similarly, 100 sensors are randomly distributed in a cube with a side length of 200 m, and the centers of the two cubes coincide. Randomly select six sensors and one target as a positioning scene, and run 5000 Monte Carlo experiments on the PC to locate each scene’s target. The settings of other parameters are the same as those of target positioning in a specific scene. The spatial distribution of targets and sensors is shown in Figure 9.

Figure 9.

Target and sensor location distribution in random scenarios.

Figure 10 depicts the relationship between the MSE and deviation of the target position estimated by the two algorithms and the variance of the sensor position error when the target and sensors are randomly distributed. It can be seen from the graph curve that when $10\lg \left({\sigma }_s^2\right)=-30\mathrm{dB}{\mathrm{m}}^2$ or $10\lg \left({\sigma }_s^2\right)=0\mathrm{dB}{\mathrm{m}}^2$ , the MSE estimated by the TSWLS algorithm for the target position suddenly deviates from the CRLB and the deviation increases sharply, this shows that the target to be located at this time may be just at or close to the reference sensor, resulting in an increase in estimation error. On the contrary, the optimization algorithm in the article can reach CRLB at the error threshold. Compared with the TSWLS algorithm, the optimization algorithm reduces the position estimation MSE by 4 and 7 dBm, respectively under the above two sensor position error variance points, and the position estimation deviation is reduced by 7 and 6 dBm, respectively. This is mainly because the optimization algorithm corrects the defects in the TSWLS algorithm by reselecting the reference sensor and rotating the coordinate system. However, in the same positioning scenario, the optimization algorithm takes a slightly longer time to locate the target than the TSWLS algorithm. This is because the optimization algorithm introduces the reference sensor’s reselection and the rotating coordinate system, which causes the amount of calculation to increase. Considering comprehensively, when randomly distributed sensors are used to locate randomly distributed targets, the optimization algorithm is better than the TSWLS algorithm in the positioning targets’ overall performance.

Figure 10.

Performance comparison of the two algorithms when targets and sensors are randomly distributed.

Figure 11 shows the relationship between the average MSE and the deviation of the target position estimated by the two positioning algorithms with the number of sensors. In the experiment, the number of sensors is gradually increased from 3 to 8. Simultaneously, set the average position error variance of the sensor to $10\lg \left({\sigma }_s^2\right)=-30\mathrm{dB}{\mathrm{m}}^2$ , and the settings of other parameters are the same as the above experiment. The experimental results show that the MSE and deviation of the two algorithms’ target position estimation decrease with the number of sensors. Especially when the number of sensors increases from 3 to 5, the MSE and deviation of the position estimation decrease sharply, but the number of sensors decreases from 6. This change will become slower when it is increased to eight, which is why six sensors are set to locate the target in the above experiment.

Figure 11.

The relationship between position estimation MSE and deviation and the number of sensors.

In the presence of both sensor position error and measurement error, this section carefully analyzes the solution process of the TSWLS algorithm and corrects the fatal defects in TSWLS by re-selecting the reference sensor and rotating the coordinate system, thus optimizing the TSWLS algorithm and improving the positioning performance of the TSWLS algorithm in general. However, the introduction of new operation steps in the optimization algorithm leads to a slight increase in the amount of operations and a slight increase in the time used to locate the target, which is a reflection of the limitations of the algorithm in this paper, and it is evident that any high-precision localization algorithm comes at the expense of the amount of operations.

Conclusion

Aiming at the positioning problem of multiple non-coincident targets, the ICTLS positioning algorithm is proposed. By analyzing the shortcomings of the TSWLS algorithm when positioning multiple targets, the proposed ICTLS algorithm establishes global constraints based on the relationship between the target and additional variables, and introduces Lagrangian multiplier technology to solve the positioning equation using the quasi-Newtonian BFGS iterative formula. Avoiding the operation of the Hessian matrix and reducing the calculation amount of the algorithm. The proposed ICTLS algorithm shows better positioning performance and strong robustness.

The optimized positioning algorithm of TSWLS solves the hidden dangers in many positioning algorithms; when the target is located at or close to the reference sensor, the estimation error increases suddenly. Based on retaining the closed-form solution, the optimization, the algorithm effectively reduces the position estimation error by reselecting the reference sensor and the rotating coordinate system, corrects the TSWLS algorithm’s defects, and comprehensively improves the positioning performance of the TSWLS algorithm.