Enhanced Moving Source Localization with Time and Frequency Difference of Arrival: Motion-Assisted Method for Sub-Dimensional Sensor Networks

Abstract

Localizing a moving source by Time Difference of Arrival (TDOA) and Frequency Difference of Arrival (FDOA) commonly requires at least $N+1$ sensors in N-dimensional space to obtain more than N pairs of TDOAs and FDOAs, thereby establishing more than $2N$ equations to solve for $2N$ unknowns. However, if there are insufficient sensors, the localization problem will become underdetermined, leading to non-unique solutions or inaccuracies in the minimum norm solution. This paper proposes a localization method using TDOAs and FDOAs while incorporating the motion model. The motion between the source and sensors increases the equivalent length of the baseline, thereby improving observability even when using the minimum number of sensors. The problem is formulated as a Maximum Likelihood Estimation (MLE) and solved through Gauss–Newton (GN) iteration. Since GN requires an initialization close to the true value, the MLE is transformed into a semidefinite programming problem using Semidefinite Relaxation (SDR) technology, while SDR results in a suboptimal estimate, it is sufficient as an initialization to guarantee the convergence of GN iteration. The proposed method is analytically shown to reach the Cramér–Rao Lower Bound (CRLB) accuracy under mild noise conditions. Simulation results confirm that it achieves CRLB-level performance when the number of sensors is lower than $N+1$, thereby corroborating the theoretical analysis.

1. Introduction

Source localization has been fundamental across various applications, including radar, sonar, and navigation [1,2,3], and has further expanded its utility into diverse commercial domains such as indoor sensing, robotics, and sensor networks [4,5]. This trend underscores the growing significance of precise positioning in modern society. Localization traditionally involves two stages. The first stage acquires parameters through signals transmitted from a source to several spatially distributed sensors. Subsequently, the second stage formulates a nonlinear optimization problem to determine the source position. Typical parameters for stationary source localization include Time Difference of Arrival (TDOA) [6,7,8,9], Time of Arrival (TOA) [10,11,12], Angle of Arrival (AOA) [13,14], and Received Signal Strength (RSS) [15,16]. Among them, TDOA stands out as one of the most popular choices. Distinguished from TOA-based technologies, TDOA localization avoids the requirement of synchronization between the source and sensors, and its performance commonly exceeds that of AOA and RSS-based localization methods.

When relative motion exists between the source and sensors, the Frequency Difference of Arrival (FDOA) can additionally be obtained. Combining FDOA with TDOA enhances the estimation accuracy of both the position and velocity of the source. However, moving source localization using TDOA and FDOA measurements presents challenges due to the highly nonlinear and nonconvex nature of the Maximum Likelihood Estimation (MLE) problem. The conventional approach to handling the nonlinear MLE problem typically involves linearizing it through Taylor-series expansion [17]. This results in an iterative algorithm that can reach the optimal estimate. However, this method requires an accurate initial guess to ensure convergence; otherwise, the iteration may be trapped by local minima. This limitation prompted researchers to seek better solutions, one of which is the closed-form solution (CFS). Two-Step Weighted Least Squares (2SWLS) [18] is perhaps the most popular CFS for solving TDOA-FDOA localization [19,20]. This approach first linearizes the measurement equations by introducing two redundant parameters after squaring, forming a pseudo-linear Weighted Least Squares (WLS) problem that can straightforwardly find the solution. The second stage leverages constraints between the redundant parameters and the source position and velocity to enhance the location accuracy. Studies such as [21,22,23] took account of the sensor uncertainties of sensor position and velocity errors, developing new CFSs to reduce estimation errors. Since sensor uncertainties are a type of common mode error that is identical for different sources, an asymptotically efficient algebraic solution for multiple disjoint sources was proposed to better calibrate the error [24]. Ref. [25] formulated the problem as a quadratic optimization with two quadratic equality constraints and solved it using iterative Constrained Weighted Least Squares (CWLS), significantly improving upon previous methods. Although the aforementioned CFSs achieve the Cramér–Rao Lower Bound (CRLB) of Mean-Square Error (MSE) and avoid possible divergence problems, they are not noise resilient if the noise is relatively large. Therefore, Multidimensional Scaling (MDS) analysis was incorporated to improve robustness against noise [26]. Another approach to enhancing noise robustness is resorting to direct position determination [27], but this comes with demanding computational complexity.

Compromising between performance and complexity, convex optimization [28,29,30] is a competitive technology superior to the aforementioned approaches. It typically formulates TDOA-FDOA localization as a semidefinite programming (SDP) problem by utilizing Semidefinite Relaxation (SDR) to relax the nonconvex problem to a nearly equivalent convex one. The SDP solution can mitigate the divergence issue that may occur in MLE implemented by iteration. While SDR sometimes relaxes the problem with insufficient tightness, resulting in a suboptimal estimate, it is sufficiently close to the true value and can serve as an appropriate initialization to guarantee the convergence of iterative MLE [31,32]. The other strategy to improve the accuracy of SDP is error correction [33], which is capable of achieving the CRB level of performance. The aforementioned studies focus on developing feasible and optimal solutions for estimating the source position, but they did not take into account the influence of sensor placement [34]. The optimal localization performance is primarily determined by the Geometric Dilution of Precision (GDOP). Optimizing the sensor placement can improve the GDOP so a better performance is achievable.

In general, the minimum number of sensors required for unambiguously localizing a moving source in an N-dimensional space using TDOA and FDOA is more than $N+1$, as the TDOA and FDOA equations are nonlinear with respect to the source position and velocity. Some studies [19,26] introduce additional instrumental parameters for pseudo-linearization, thus necessitating more sensors to ensure the proposed closed-form algorithms are solvable. Other researchers have explored localization with the minimum number of sensors. For instance, $N+1$ sensors are required for TDOA localization [35], while elliptic localization in Multiple-Input Multiple-Output (MIMO) systems requires at least one transmitter and three receivers or two transmitters and two receivers [36,37]. These methods leverage measurements from sensors or receivers extracted from one frame of snapshots to establish the equations; thus, they are limited by the actual number of sensors or receivers. Typically, to guarantee that the optimization problem is overdetermined and the solution for source location is unique, it is necessary to incorporate more than $N+1$ sensors into the localization system in most scenarios [38,39].

When the number of sensors is limited, one possible method to localize a moving source is known as Target Motion Analysis (TMA). TMA is a technique used in various fields, notably in military and maritime operations, to determine the movement characteristics of a target based on observations from one or more sensors. The primary goal of TMA is to estimate the position, velocity, and possibly other relevant parameters of a target, such as its course and acceleration. TMA incorporates a kinetic model, allowing it to handle moving source localization with fewer sensors. The early development of TMA evolved from bearing-only (BO) measurements [40]. The team led by K. Doğançay made significant contributions to BO-TMA closed-form solutions, addressing issues such as two-dimensional single-station TMA [41], bias reduction in target motion parameter estimation [42], and the effect of receiver station location errors [43]. These research findings were further extended to three-dimensional targets [44]. TMA is not only designed for bearing-only measurements but is also applicable to TDOA [45] and FDOA [46]. In a scenario where the number of sensors is smaller than the dimension of space, leveraging TMA can help handle constant velocity moving source localization using TDOA and FDOA measurements.

This paper focuses on the TDOA-FDOA localization problem when the number of sensors is less than $N+1$ for an N-dimensional scenario. By incorporating the kinetic model of the source, which is an effective method involving TMA, it becomes feasible to solve for $2N$ unknowns, such as position and velocity. The sensors observe multiple times while the source is moving, obtaining multiple sets of TDOAs and FDOAs. This paper first formulates the moving source localization as a MLE problem, which is nonlinear and nonconvex. The MLE is then implemented by Gauss–Newton (GN) iteration. To find an initialization to start the GN iteration, the original MLE is then relaxed to a convex optimization problem through SDR technology, resulting in an SDP problem with penalties. The solution can be found straightforwardly using mature tools such as the CVX Toolbox [47]. Although the SDP solution is suboptimal, it is sufficient to guarantee that the GN iteration converges to an optimal point if the noise is relatively mild. Mean-square analysis shows that the performance of the proposed estimator achieves the Cramér–Rao Lower Bound (CRLB), and the bias is close to zero in the small noise region. Simulations validate the analytical results and additionally demonstrate the performance behavior while the observation time, source range, and source velocity vary.

The paper begins with the measurement model of TDOA and FDOA in Section 2. Section 3 introduces the formulation of the MLE problem for TDOA-FDOA localization and proposes a solution by GN iteration, where the initial solution for both velocity and position is subsequently found by SDP. Section 4 provides the CRLB and analyzes the performance of the proposed MLE. Section 5 presents the simulation results, and Section 6 concludes the paper.

Notations: The following notations are used throughout the paper. Bold uppercase and bold lowercase letters denote matrices and vectors, respectively. ${\mathbf{I}}_m$ is an $m\times m$ identity matrix, and ${\mathbf{1}}_m$ is an $m\times 1$ vector whose elements are all ones. ${\mathbf{O}}_{m\times n}$ is an $m\times n$ zero matrix, and $\mathbf{0}$ is an all-zero vector with appropriate size. ${\ast }^T$ is the transposition of a matrix or vector. $E(·)$ denotes expectation, and $\parallel ·\parallel$ is the $l_2$-norm. $\mathbf{A}(:,i)$ denotes the i-th column of $\mathbf{A}$, and $\mathbf{A}(i,j)$ denotes the $(i,j)$-th element of $\mathbf{A}$. $\mathbf{a}\left(i\right)$ denotes the i-th element of $\mathbf{a}$. $\mathbf{A}⪰\mathbf{B}$ means that $\mathbf{A}-\mathbf{B}$ is positive semidefinite. $tr\left\{\mathbf{A}\right\}$ represents the trace of $\mathbf{A}$. ⊗ is the Kronecker product between two matrices. The symbols with a superscript ${ }^o$ represent the true value corresponding to the noisy one. The vectors in this paper are column vectors.

The main acronyms used in this paper are summarized in

Table 1. Acronym list.

Acronym	Description
TDOA	Time Difference of Arrival
TOA	Time of Arrival
AOA	Angle of Arrival
RSS	Received Signal Strength
FDOA	Frequency Difference of Arrival
MLE	Maximum Likelihood Estimation
CFS	Closed-form solution
2SWLS	Two-Step Weighted Least Squares
WLS	Weighted Least Squares
CWLS	Constrained Weighted Least Squares
CRLB	Cramér–Rao Lower Bound
MSE	Mean-Square Error
MDS	Multidimensional Scaling
SDP	Semidefinite programming
SDR	Semidefinite Relaxation
MIMO	Multiple-Input Multiple-Output
TMA	Target Motion Analysis
BO	Bearing-only
GN	Gauss–Newton
QN	Quasi-Newton
LM	Levenberg–Marquardt
LMI	Linear Matrix Inequality
SOC	Second-order cone
FIM	Fisher Information Matrix
RMSE	Root mean-squared error
UAV	Unmanned aerial vehicle
GDOP	Geometric Dilution of Precision

.

Table 2. Categories of the related works.

Category	Method	Paper
TDOA and FDOA	CFS	[18,19,20,21,22,23,24,25,26]
TDOA and FDOA	Iterative solution	[27]
TDOA and FDOA	SDP solution	[28,29,30]
BO-TMA	Particle filter	[40]
BO-TMA	CFS	[41,42,43,44]
TMA for TDOA only	Iterative solution	[45]
TMA for FDOA only	CFS and SDP	[46]
TMA for TDOA and FDOA	SDP and GN	Proposed

summarizes the related works and categorizes them.

2. Localization Scenario

Consider a scenario as shown in Figure 1 where M moving sensors are distributed in an N-dimensional space ($N=3$ or $N=2$). The known initial positions and velocities of these sensors are denoted as ${\mathbf{s}}_{i,1}\in {\mathbb{R}}^N$ and ${\dot{\mathbf{s}}}_i\in {\mathbb{R}}^N$, respectively. These sensors are deployed to determine the position and velocity of a moving source. Assuming synchronization in time and frequency among the sensors, they can produce TDOA and FDOA measurements between receiving sensors and a reference sensor. Without loss of generality, the reference sensor is chosen as ${\mathbf{s}}_1$. The unknown initial position and velocity of the source are denoted as ${\mathbf{u}}^o\in {\mathbb{R}}^N$ and ${\dot{\mathbf{u}}}^o\in {\mathbb{R}}^N$, respectively, which will be estimated using TDOA and FDOA measurements. The velocities of the source and sensors are constant. During the motion of the sensors, the source is observed a total of K times, with a time interval of $\tau$ between the adjacent observation frame. Assuming the observation time is short, the sensors and source are considered as quasi-static during one observation frame. At the k-th observation frame, the positions of the sensors are denoted as

$${\mathbf{s}}_{i,k}={\mathbf{s}}_{i,1}+{\dot{\mathbf{s}}}_i(k-1)\tau ,k=1,2,\cdots ,K,$$

(1)

and the position of the source is

$${\mathbf{u}}_k^o={\mathbf{u}}^o+{\dot{\mathbf{u}}}^o(k-1)\tau .$$

(2)

Collecting the initial position and velocity of the source together, the unknown vector is ${\mathit{\beta }}^o={[{\mathbf{u}}^{oT},{\dot{\mathbf{u}}}^{oT}]}^T$.

The TDOAs between the i-th and 1st sensors are interchangeable with range differences by multiplying with the signal propagation speed,

$$r_{i,k}^o=d_{i,k}^o-d_{1,k}^o,i=2,\cdots ,M.$$

(3)

where

$$d_{i,k}^o=\parallel {\mathbf{u}}_k^o-{\mathbf{s}}_{i,k}\parallel ,i=1,\cdots ,M,$$

(4)

is the true distance from the source to the i-th sensor during the k-th observation frame. The rate of the time derivative of the true range difference in (3) is

$${\dot{r}}_{i,k}^o={\dot{d}}_{i,k}^o-{\dot{d}}_{1,k}^o,i=2,\cdots ,M,$$

(5)

where ${\dot{d}}_{i,k}^o$ is the rate of the range between the source and the i-th sensor in the k-th observation frame

$${\dot{d}}_{i,k}^o=\frac{{(\dot{\mathbf{u}}-{\dot{\mathbf{s}}}_i)}^T({\mathbf{u}}_k^o-{\mathbf{s}}_{i,k})}{\parallel {\mathbf{u}}_k^o-{\mathbf{s}}_{i,k}\parallel },i=1,\cdots ,M,$$

(6)

The FDOA between the i-th sensor and the reference is the range difference rate divided by the signal propagation speed c,

$$f_{i,k}^o={\dot{r}}_{i,k}^o/c.$$

(7)

When the noise exists, the measured TDOA and FDOA are

$$\begin{array}{cc}\hfill r_{i,k}& =r_{i,k}^o+n_{i,k},\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill {\dot{r}}_{i,k}& ={\dot{r}}_{i,k}^o+e_{i,k},\hfill \end{array}$$

(9)

where $n_{i,k}$ is the measurement noise of range difference and $e_k$ is the measurement noise of range difference rate. The collections of TDOA and FDOA measurements form

$$\begin{array}{cc}\hfill \mathbf{r}& ={\mathbf{r}}^o+\mathbf{n}={[r_{2,1}^o,r_{3,1}^o,\cdots ,r_{M,1}^o,\cdots ,r_{M,K}^o]}^T+\mathbf{n},\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill \dot{\mathbf{r}}& ={\dot{\mathbf{r}}}^o+\mathbf{e}={[{\dot{r}}_{2,1}^o,{\dot{r}}_{3,1}^o,\cdots ,{\dot{r}}_{M,1}^o,,\cdots ,{\dot{r}}_{M,K}^o]}^T+\mathbf{e},\hfill \end{array}$$

(11)

where

$$\begin{array}{cc}\hfill \mathbf{n}& ={[n_{2,1},n_{3,1},\cdots ,n_{M,1},\cdots ,n_{M,K}]}^T,\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill \mathbf{e}& ={[e_{2,1},e_{3,1},\cdots ,e_{M,1},\cdots ,e_{M,K}]}^T,\hfill \end{array}$$

(13)

are the noise vectors obeying a zero-mean Gaussian distribution. Putting the two sets of measurements together yields the total measurements vector $\mathit{\alpha }={[{\mathbf{r}}^T,{\dot{\mathbf{r}}}^T]}^T$. The corresponding measurement error vector $\Delta \mathit{\alpha }={[{\mathbf{n}}^T,{\mathbf{e}}^T]}^T$ is also zero-mean and has the covariance matrix

$$E[\Delta \mathit{\alpha }\Delta {\mathit{\alpha }}^T]={\mathbf{Q}}_{\mathit{\alpha }}.$$

(14)

3. Localization Method

MLE is widely utilized for parameter estimation in statistical modeling due to its desirable properties like consistency, efficiency, and asymptotic normality. It offers accurate, optimal, and unbiased estimation for source localization, effectively leveraging observation data while retaining flexibility to different types of source localization problems. Supported by theoretical foundations and mathematical derivations, MLE plays a crucial role in the development and application of source localization technologies. Maximizing the likelihood function of ${\mathbf{u}}^o$ and ${\dot{\mathbf{u}}}^o$ is equivalent to minimizing the objective function under Gaussian noise conditions [48],

$$\begin{array}{cc}\hfill J\left(\mathit{\beta }\right)& ={(\mathit{\alpha }-\overline{\mathit{\alpha }}\left(\mathit{\beta }\right))}^T{\mathbf{Q}}_{\mathit{\alpha }}^{-1}(\mathit{\alpha }-\overline{\mathit{\alpha }}\left(\mathit{\beta }\right))\hfill \\ \hfill & =\sum _{k_1,k_2=1}^K\sum _{i_1,i_2=2}^M{\omega }_{a(i_1,k_1),b(i_2,k_2)}\left(r_{i_1,k_1}-{\overline{r}}_{i_1,k_1}\right)\left(r_{i_2,k_2}-{\overline{r}}_{i_2,k_2}\right)\hfill \\ \hfill & +\sum _{k_1,k_2=1}^K\sum _{i_1,i_2=2}^M{\omega }_{c(i_1,k_1),d(i_2,k_2)}\left({\dot{r}}_{i_1,k_1}-{\overline{\dot{r}}}_{i_1,k_1}\right)\left({\dot{r}}_{i_2,k_2}-{\overline{\dot{r}}}_{i_2,k_2}\right),\hfill \end{array}$$

(15)

where $a(i_1,k_1)=(k_1-1)(M-1)+i_1-1$, $b(i_2,k_2)=(k_2-1)(M-1)+i_2-1$, $c(i_1,k_1)=K(M-1)+a(i_1,k_1)$, $d(i_2,k_2)=K(M-1)+b(i_2,k_2)$ and ω_ı,ȷ is the (ı,ȷ)-th element of ${\mathbf{Q}}_{\mathit{\alpha }}^{-1}$. $\overline{\mathit{\alpha }}\left(\mathit{\beta }\right)$, as well as ${\overline{r}}_{i_1,k_1}$ and ${\overline{\dot{r}}}_{i_1,k_1}$, are the reconstructed measurements from $\mathit{\beta }$ according to (3) and (5).

The MLE problem above is nonconvex and nonlinear, which means it is nontrivial to find the optimal solution. Besides the exhaustive method, an effective strategy to handle the MLE is iterative solution. It should be noted that iteration methods, such as GN, Newton–Raphson, Quasi-Newton (QN) and Levenberg–Marquardt (LM), involve asymptotic MLE. Their solutions can approach the CRLB when the noise is mild but diverge or are trapped by local minima if there is a relatively large amount of noise. The GN method offers advantages over others as it balances the convergence and computational costs by computing Jacobian matrices instead of Hessians, suitable for large-scale and high-dimensional parameter spaces, making it particularly effective for highly nonlinear problems while being less computationally demanding and more scalable, although it may be slightly sensitive to initial conditions.

Therefore, the proposed MLE consists of two parts. One part is an iterative equation, which was developed for minimizing (15) by means of the GN implementation [48]. Considering the initialization sensitivity of GN iteration, the other part provides an suboptimal solution for ${\mathit{\beta }}^o$ by SDR and SDP. This solution utilized to initialize the iterative process is sufficient for guaranteeing the convergence.

3.1. Gauss–Newton Iteration

Taking the Taylor-series expansion of $\mathit{\alpha }\left(\mathit{\beta }\right)$ at ${\mathit{\beta }}^{\left(0\right)}$ and retaining the first-order terms give a quadratic approximation of $J\left(\mathit{\beta }\right)$

$$J\left(\mathit{\beta }\right)\approx {\left[\mathit{\alpha }-\overline{\mathit{\alpha }}\left({\mathit{\beta }}^{\left(0\right)}\right)-{\left.\frac{\partial \mathit{\alpha }}{\partial \mathit{\beta }}\right|}_{\mathit{\beta }={\mathit{\beta }}^{\left(0\right)}}(\mathit{\beta }-{\mathit{\beta }}^{\left(0\right)})\right]}^T{\mathbf{Q}}_{\mathit{\alpha }}^{-1}\left[\mathit{\alpha }-\overline{\mathit{\alpha }}\left({\mathit{\beta }}^{\left(0\right)}\right)-{\left.\frac{\partial \mathit{\alpha }}{\partial \mathit{\beta }}\right|}_{\mathit{\beta }={\mathit{\beta }}^{\left(0\right)}}(\mathit{\beta }-{\mathit{\beta }}^{\left(0\right)})\right].$$

(16)

Differentiating it with respect to $\mathit{\beta }$ and setting the result to zero yield the GN iterative equation for asymptotically minimizing $J\left(\mathit{\beta }\right)$, after replacing the superscript,

$${\mathit{\beta }}^{(g+1)}={\mathit{\beta }}^{\left(g\right)}+{\left({\mathbf{G}}^{\left(g\right)T}\mathbf{W}{\mathbf{G}}^{\left(g\right)}\right)}^{-1}{\mathbf{G}}^{\left(g\right)T}\mathbf{W}\left(\mathit{\alpha }-{\mathit{\alpha }}^{\left(g\right)}\right),g=0,1,\cdots ,L.$$

(17)

where g is the iteration count and L is the number of iterations. ${\mathit{\beta }}^{\left(0\right)}$ is the initial solution guess whose choice will be elaborated later, which should not be confused with ${\mathit{\beta }}^o$ that represents the true value. In (17), $\mathbf{W}={\mathbf{Q}}_{\mathit{\alpha }}^{-1}$ and ${\mathit{\alpha }}^{\left(g\right)}$ are the re-generated measurement vectors with ${\mathit{\beta }}^{\left(g\right)}$ replacing ${\mathit{\beta }}^o$ in ${\mathit{\alpha }}^o$. ${\mathbf{G}}^{\left(g\right)}$ is the Jacobian matrix given by

$${\mathbf{G}}^{\left(g\right)}={\left.\frac{\partial \mathit{\alpha }}{\partial {\mathit{\beta }}^T}\right|}_{\mathit{\beta }={\mathit{\beta }}^{\left(g\right)}}=\left[\begin{array}{cc}{\displaystyle \frac{\partial \mathbf{r}}{\partial {\mathbf{u}}^{\left(g\right)T}}}& {\displaystyle \frac{\partial \mathbf{r}}{\partial {\dot{\mathbf{u}}}^{\left(g\right)T}}}\\ {\displaystyle \frac{\partial \dot{\mathbf{r}}}{\partial {\mathbf{u}}^{\left(g\right)T}}}& {\displaystyle \frac{\partial \dot{\mathbf{r}}}{\partial {\dot{\mathbf{u}}}^{\left(g\right)T}}}\end{array}\right]=\left[\begin{array}{cc}{\mathbf{R}}^{\left(g\right)}& {\mathbf{B}}^{\left(g\right)}\\ {\dot{\mathbf{R}}}^{\left(g\right)}& {\dot{\mathbf{B}}}^{\left(g\right)}\end{array}\right].$$

(18)

The $\left(\right(k-1\left)\right(M-1)+(i-1\left)\right)$-th row of ${\mathbf{R}}^{\left(g\right)}$ is

$$\begin{array}{cc}\hfill & {\mathbf{R}}^{\left(g\right)}\left((k-1)(M-1)+(i-1),:\right)=\frac{{\left({\mathbf{u}}_k^{\left(g\right)}-{\mathbf{s}}_{i,k}\right)}^T}{d_{i,k}^{\left(g\right)}}-\frac{{\left({\mathbf{u}}_k^{\left(g\right)}-{\mathbf{s}}_{1,k}\right)}^T}{d_{1,k}^{\left(g\right)}},\hfill \\ \hfill & k=1,\cdots ,K,i=2,\cdots ,M,\hfill \end{array}$$

(19)

where $d_{i,k}^{\left(g\right)}$ and $d_{1,k}^{\left(g\right)}$ are from (4) with the true value ${\mathbf{u}}^o$ replaced by ${\mathbf{u}}^{\left(g\right)}$, and similarly the $\left(\right(k-1\left)\right(M-1)+(i-1\left)\right)$-th row of ${\dot{\mathbf{R}}}^{\left(g\right)}$ is

$$\begin{array}{cc}\hfill & {\dot{\mathbf{R}}}^{\left(g\right)}((k-1)(M-1)+(i-1),:)\hfill \\ \hfill & =\frac{{\dot{\mathbf{u}}}^{\left(g\right)T}-{\dot{\mathbf{s}}}_i^T}{d_{i,k}^{\left(g\right)}}-\frac{{\dot{d}}_{i,k}^{\left(g\right)}}{d_{i,k}^{\left(g\right)2}}{\left({\mathbf{u}}_k^{\left(g\right)}-{\mathbf{s}}_{i,k}\right)}^T-\frac{{\dot{\mathbf{u}}}^{\left(g\right)T}-{\dot{\mathbf{s}}}_1^T}{d_{1,k}^{\left(g\right)}}+\frac{{\dot{d}}_{1,k}^{\left(g\right)}}{d_{1,k}^{\left(g\right)2}}{\left({\mathbf{u}}_k^{\left(g\right)}-{\mathbf{s}}_{1,k}\right)}^T,\hfill \end{array}$$

(20)

The $\left(\right(k-1\left)\right(M-1)+(i-1\left)\right)$-th row of ${\mathbf{B}}^{\left(g\right)}$ is

$$\begin{array}{cc}\hfill {\mathbf{B}}^{\left(g\right)}((k-1)(M-1)+(i-1),:)& =\tau (k-1){\mathbf{R}}^{\left(g\right)}((k-1)(M-1)+(i-1),:)\hfill \\ \hfill & =\tau (k-1)\left(\frac{{\left({\mathbf{u}}_k^{\left(g\right)}-{\mathbf{s}}_{i,k}\right)}^T}{d_{i,k}^{\left(g\right)}}-\frac{{\left({\mathbf{u}}_k^{\left(g\right)}-{\mathbf{s}}_{1,k}\right)}^T}{d_{1,k}^{\left(g\right)}}\right),\hfill \end{array}$$

(21)

and the $\left(\right(k-1\left)\right(M-1)+(i-1\left)\right)$-th row of ${\dot{\mathbf{B}}}^{\left(g\right)}$ is

$$\begin{array}{cc}\hfill & {\dot{\mathbf{B}}}^{\left(g\right)}((k-1)(M-1)+(i-1),:)\hfill \\ \hfill & =\frac{{(\tau (k-1)({\dot{\mathbf{u}}}^{\left(g\right)}-{\dot{\mathbf{s}}}_i)+({\mathbf{u}}_k^{\left(g\right)}-{\mathbf{s}}_{i,k}))}^T}{d_{i,k}^{\left(g\right)}}-\frac{\tau (k-1){({\mathbf{u}}_k^{\left(g\right)}-{\mathbf{s}}_{i,k})}^T{\dot{r}}_{i,k}^{\left(g\right)}}{d_{i,k}^{\left(g\right)2}}\hfill \\ \hfill & -\frac{{(\tau (k-1)({\dot{\mathbf{u}}}^{\left(g\right)}-{\dot{\mathbf{s}}}_1)+({\mathbf{u}}_k^{\left(g\right)}-{\mathbf{s}}_{1,k}))}^T}{d_{1,k}^{\left(g\right)}}-\frac{\tau (k-1){({\mathbf{u}}_k^{\left(g\right)}-{\mathbf{s}}_{1,k})}^T{\dot{r}}_{1,k}^{\left(g\right)}}{d_{1,k}^{\left(g\right)2}}.\hfill \end{array}$$

(22)

Ensuring that the initial solution ${\mathit{\beta }}^{\left(0\right)}$ is sufficiently close to the true value is crucial for initializing the iteration process. Failure to achieve this closeness may result in either divergence or convergence to a local minimum for (17). Our attention is now directed towards obtaining this essential initial solution.

3.2. Initial Solution

To guarantee convergence of the GN iteration, this paper introduces an SDR method that transforms the nonconvex ML problem into a convex one, thereby forming an SDP problem. Let us start by considering the ML cost function (15) and reformulating it as a constrained optimization problem,

$$\begin{array}{cc}\hfill & \underset{\mathbf{u},\dot{\mathbf{u}},\mathbf{d},\dot{\mathbf{d}}}{min}{(\mathbf{r}-\mathbf{A}\mathbf{d})}^T{\mathbf{Q}}^{-1}(\mathbf{r}-\mathbf{A}\mathbf{d})+{(\dot{\mathbf{r}}-\mathbf{A}\dot{\mathbf{d}})}^T{\dot{\mathbf{Q}}}^{-1}(\dot{\mathbf{r}}-\mathbf{A}\dot{\mathbf{d}})\hfill \\ \hfill & \mathrm{s}.\mathrm{t}.d_{i,k}=∥{\mathbf{u}}_k-{\mathbf{s}}_{i,k}∥,\hfill \\ \hfill & {\dot{d}}_{i,k}=\frac{{(\dot{\mathbf{u}}-{\dot{\mathbf{s}}}_i)}^T({\mathbf{u}}_k-{\mathbf{s}}_{i,k})}{∥{\mathbf{u}}_k-{\mathbf{s}}_{i,k}∥},i=1,\cdots ,M,k=1,\cdots ,K,\hfill \end{array}$$

(23)

where

$$\mathbf{d}={[d_{1,1},d_{2,1},\cdots ,d_{M,1},\cdots ,d_{M,K}]}^T,$$

(24)

$$\dot{\mathbf{d}}={[{\dot{d}}_{1,1},{\dot{d}}_{2,1},\cdots ,{\dot{d}}_{M,1},\cdots ,{\dot{d}}_{M,K}^o]}^T,$$

(25)

and $\mathbf{A}={\mathbf{I}}_K⨂\mathbf{B}$, $\mathbf{B}=[-{\mathbf{1}}_{M-1},{\mathbf{I}}_{M-1}]$, $\mathbf{Q}=E\left[\mathbf{n}{\mathbf{n}}^T\right]$, $\dot{\mathbf{Q}}=E\left[\mathbf{e}{\mathbf{e}}^T\right]$. Denote $\mathbf{h}={[{\mathbf{d}}^T,{\dot{\mathbf{d}}}^T]}^T,{\mathbf{A}}_1=\mathbf{A}[{\mathbf{I}}_{MK},{\mathbf{O}}_{MK\times MK}]$ and ${\mathbf{A}}_2=\mathbf{A}[{\mathbf{O}}_{MK\times MK},{\mathbf{I}}_{MK}]$. The problem given in (23) can be rewritten as

$$\begin{array}{ccc}\hfill & \underset{\begin{array}{c}\mathbf{u},\dot{\mathbf{u}},\mathbf{h}\end{array}}{min}{(\mathbf{r}-{\mathbf{A}}_1\mathbf{h})}^T{\mathbf{Q}}^{-1}(\mathbf{r}-{\mathbf{A}}_1\mathbf{h})+{(\dot{\mathbf{r}}-{\mathbf{A}}_2\mathbf{H})}^T{\dot{\mathbf{Q}}}^{-1}(\dot{\mathbf{r}}-{\mathbf{A}}_2\mathbf{h})\hfill & \hfill \left(26\mathrm{a}\right)\\ \hfill & \mathrm{s}.\mathrm{t}.\mathbf{h}((k-1)M+i)=∥{\mathbf{u}}_k-{\mathbf{s}}_{i,k}∥,\hfill & \hfill \left(26\mathrm{b}\right)\\ \hfill & \mathbf{h}((k-1+K)M+i)=\frac{{(\dot{\mathbf{u}}-{\dot{\mathbf{s}}}_i)}^T({\mathbf{u}}_k-{\mathbf{s}}_{i,k})}{∥{\mathbf{u}}_k-{\mathbf{s}}_{i,k}∥},\hfill & \hfill \left(26\mathrm{c}\right)\\ & i=1,\cdots ,M,k=1,\cdots ,K.\hfill & \hfill \end{array}$$

The objective function (26a) for minimization is nonconvex. The convex equivalent form is

$$tr\left[({\mathbf{A}}_1^T{\mathbf{Q}}^{-1}{\mathbf{A}}_1+{\mathbf{A}}_2^T{\dot{\mathbf{Q}}}^{-1}{\mathbf{A}}_2)\mathbf{H}\right]-2{\mathbf{H}}^T({\mathbf{A}}_1^T{\mathbf{Q}}^{-1}\mathbf{r}+{\mathbf{A}}_2^T{\dot{\mathbf{Q}}}^{-1}\dot{\mathbf{r}})$$

(27)

where $\mathbf{H}=\mathbf{h}{\mathbf{h}}^T$ and $tr(·)$ represents the trace of a matrix. Two constant terms (26b) and (26c) are discarded without affecting the results in (27).

Acknowledging that the constraints (26b) and (26c) are nonconvex, they must be relaxed into convex constraints while maintaining a tight connection with the original ones. Let $\mathbf{X}=\left[\mathbf{u},\tau \dot{\mathbf{u}}\right]$ and $\mathbf{Y}={\mathbf{X}}^T\mathbf{X}$. The constraint (26b) takes on an equivalent form as follows:

$$\mathbf{H}((k-1)M+i,(k-1)M+i)=\mathbf{h}{((k-1)M+i)}^2={\mathbf{q}}^T\mathbf{Y}\mathbf{q}-2{\mathbf{q}}^T{\mathbf{X}}^T{\mathbf{s}}_{i,k}+{\mathbf{s}}_{i,k}^T{\mathbf{s}}_{i,k},$$

(28)

for $i=1,\cdots ,M$, $k=1,\cdots ,K$, where $\mathbf{q}={\left[1,k-1\right]}^T$ thus ${\mathbf{u}}_k=\mathbf{X}\mathbf{q}$. Applying the Cauchy–Schwartz inequality to (26b) obtains

$$\mathbf{H}(l_1,l_2)\ge |{\mathbf{q}}_1^T\mathbf{Y}{\mathbf{q}}_2-{\mathbf{q}}_1^T{\mathbf{X}}^T{\mathbf{s}}_{i_2,k_2}-{\mathbf{q}}_2^T{\mathbf{X}}^T{\mathbf{s}}_{i_1,k_1}+{\mathbf{s}}_{i_1,k_1}^T{\mathbf{s}}_{i_2,k_2}|,$$

(29)

where ${\mathbf{q}}_1={\left[1,k_1-1\right]}^T$, ${\mathbf{q}}_2={\left[1,k_2-1\right]}^T$. The variables $i_1$, $i_2$, $k_1$, and $k_2$ satisfy $1\le l_1<l_2\le MK$ where $l_1=(k_1-1)M+i_1$, $l_2=(k_2-1)M+i_2$.

Considering the difficulty in formulating constraint (26c) in a convex form due to its fraction, multiplying it by (26b) yields

$$\mathbf{h}((k-1)M+i)\mathbf{h}((k-1+K)M+i)={\dot{\mathbf{u}}}^T\mathbf{u}-{\dot{\mathbf{u}}}^T{\mathbf{s}}_i-{\mathbf{s}}_i^T\mathbf{u}+{\dot{\mathbf{s}}}_i^T{\mathbf{s}}_i.$$

(30)

Thus, the nonconvex constraint above can be expressed as

$$\begin{array}{cc}\hfill & \mathbf{H}\left(\right(k-1)M+i,(k-1+K)M+i)\hfill \\ \hfill & =\mathbf{Y}(2,1)/\tau +(k-1)\mathbf{Y}(2,2)/\tau -\mathbf{X}{(:,2)}^T{\mathbf{s}}_{i,k}/\tau -{\mathbf{q}}^T\mathbf{X}{(:,1)}^T{\dot{\mathbf{s}}}_i+{\mathbf{s}}_{i,k}^T{\dot{\mathbf{s}}}_i.\hfill \end{array}$$

(31)

Additionally, applying the Cauchy–Schwartz inequality to (26c) obtains an additional constraint,

$$\left|\mathbf{h}((k-1+K)M+i)\right|\le \parallel \dot{\mathbf{u}}-{\dot{\mathbf{s}}}_i\parallel .$$

(32)

Squaring both sides of (32) results in

$$\mathbf{H}((k-1+K)M+i,(k-1+K)M+i)\le \mathbf{Y}(2,2)/{\tau }^2-2\mathbf{X}{(:,2)}^T{\dot{\mathbf{s}}}_i/\tau +{\dot{\mathbf{s}}}_i^T{\dot{\mathbf{s}}}_i.$$

(33)

At this stage, two nonconvex constraints remain: $\mathbf{H}=\mathbf{h}{\mathbf{h}}^T$ and $\mathbf{Y}={\mathbf{X}}^T\mathbf{X}$. The SDR method can be employed [49] to relax these constraints into convex inequalities $\mathbf{H}⪰\mathbf{h}{\mathbf{h}}^T$ and $\mathbf{Y}⪰{\mathbf{X}}^T\mathbf{X}$, which can be expressed as Linear Matrix Inequalities (LMIs):

$$\begin{array}{cc}\hfill \left[\begin{array}{cc}1& {\mathbf{h}}^T\\ \mathbf{h}& \mathbf{H}\end{array}\right]& ⪰\mathbf{0},\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill \left[\begin{array}{cc}{\mathbf{I}}_N& \mathbf{X}\\ {\mathbf{X}}^T& \mathbf{Y}\end{array}\right]& ⪰\mathbf{0},\hfill \end{array}$$

(35)

where the rank-1 constraints are ignored.

Since the rank of $({\mathbf{A}}_1^T{\mathbf{Q}}^{-1}{\mathbf{A}}_1+{\mathbf{A}}_2^T{\dot{\mathbf{Q}}}^{-1}{\mathbf{A}}_2)$ is equal to $2K(M-1)$, the convex optimization formulation (27)–(29), (31), and (32) may exhibit ambiguities. To address this, two additional penalties ${\eta }_{1}tr\left(\mathbf{H}(1:MK,1:MK)\right)$ and ${\eta }_{2}tr\left(\mathbf{H}(MK+1:2MK,MK+1:2MK)\right)$ need to be introduced to tighten the constraints, where ${\eta }_1$ and ${\eta }_2$ are penalty factors. Additionally, second-order cone (SOC) constraints [50]

$$∥\mathbf{X}(:,1)\mathbf{q}-{\mathbf{s}}_{i,k}∥\le \mathbf{h}((k-1)M+i),i=1,\cdots ,M,$$

(36)

are incorporated to enhance the location accuracy. The discussions above lead to the proposed SDP problem as follows,

$$\begin{array}{cc}\hfill & \underset{\mathbf{h},\mathbf{H},\mathbf{X},\mathbf{Y}}{min}tr\left[({\mathbf{A}}_1^T{\mathbf{Q}}^{-1}{\mathbf{A}}_1+{\mathbf{A}}_2^T{\dot{\mathbf{Q}}}^{-1}{\mathbf{A}}_2)\mathbf{H}\right]-2{\mathbf{H}}^T({\mathbf{A}}_1^T{\mathbf{Q}}^{-1}\mathbf{r}+{\mathbf{A}}_2^T{\dot{\mathbf{Q}}}^{-1}\dot{\mathbf{r}})\hfill \\ \hfill & +{\eta }_{1}tr\left(\mathbf{H}(1:MK,1:MK)\right)+{\eta }_{2}tr\left(\mathbf{H}(MK+1:2MK,MK+1:2MK)\right),\hfill \\ \hfill & \mathrm{s}.\mathrm{t}.\mathbf{H}((k-1)M+i,(k-1)M+i)=\mathbf{h}{((k-1)M+i)}^2={\mathbf{q}}^T\mathbf{Y}\mathbf{q}-2{\mathbf{q}}^T{\mathbf{X}}^T{\mathbf{s}}_{i,k}+{\mathbf{s}}_{i,k}^T{\mathbf{s}}_{i,k},\hfill \\ \hfill & \mathbf{H}(l_1,l_2)\ge |{\mathbf{q}}_1^T\mathbf{Y}{\mathbf{q}}_2-{\mathbf{q}}_1^T{\mathbf{X}}^T{\mathbf{s}}_{i_2,k_2}-{\mathbf{q}}_2^T{\mathbf{X}}^T{\mathbf{s}}_{i_1,k_1}+{\mathbf{s}}_{i_1,k_1}^T{\mathbf{s}}_{i_2,k_2}|,\hfill \\ \hfill & \mathbf{H}\left(\right(k-1)M+i,(k-1+K)M+i)\hfill \\ \hfill & =\mathbf{Y}(2,1)/\tau +(k-1)\mathbf{Y}(2,2)/\tau -\mathbf{X}{(:,2)}^T{\mathbf{s}}_{i,k}/\tau -{\mathbf{q}}^T\mathbf{X}{(:,1)}^T{\dot{\mathbf{s}}}_i+{\mathbf{s}}_{i,k}^T{\dot{\mathbf{s}}}_i.\hfill \\ \hfill & \mathbf{H}((k-1+K)M+i,(k-1+K)M+i)\le \mathbf{Y}(2,2)/{\tau }^2-2\mathbf{X}{(:,2)}^T{\dot{\mathbf{s}}}_i/\tau +{\dot{\mathbf{s}}}_i^T{\dot{\mathbf{s}}}_i.\hfill \\ \hfill & \left[\begin{array}{cc}1& {\mathbf{h}}^T\\ \mathbf{h}& \mathbf{H}\end{array}\right]⪰\mathbf{0},\hfill \\ \hfill & \left[\begin{array}{cc}{\mathbf{I}}_N& \mathbf{X}\\ {\mathbf{X}}^T& \mathbf{Y}\end{array}\right]⪰\mathbf{0},\hfill \\ \hfill & ∥\mathbf{X}(:,1)\mathbf{q}-{\mathbf{s}}_{i,k}∥\le \mathbf{h}((k-1)M+i),i=1,\cdots ,M.\hfill \end{array}$$

(37)

This problem can be solved using the CVX Toolbox [47]. Therefore, the initial solution for the GN iteration described above is ${\mathit{\beta }}^{\left(0\right)}={[\mathbf{X}{(:,1)}^T,\mathbf{X}{(:,2)}^T/\tau ]}^T$.

Suitable values of ${\eta }_1$ and ${\eta }_2$ are crucial for achieving good estimation performance. However, obtaining the optimal values analytically is challenging, as they depend on various factors such as the distance, relative velocity between the sensor and the source nodes, and the noise level [51]. In each scenario, the first set of measurements is utilized to select suitable parameters. Specifically, $C_1$ values of ${\eta }_1$ are denoted as ${\eta }_{1,c_1}$, $c_1=1,2,\dots ,C_1$, and $C_2$ values of ${\eta }_2$ are denoted as ${\eta }_{2,c_2}$, $c_2=1,2,\dots ,C_2$. For each pair of penalty factors, the (37) is computed in parallel. The $({\eta }_{1,c_1},{\eta }_{2,c_2})$ giving the estimation ${\mathbf{X}}_{c_1,c_2}$ that minimizes $J\left(\mathit{\beta }\right)$ in (15) were selected for following simulations.

4. Analysis

In this section, the CRLB for position and velocity estimation will be provided, and the performance of the proposed MLE will be investigated. The covariance analysis will be limited to the first order, disregarding second- and higher-order noise terms in the estimate, which is valid only within a small error region [8,18]. The bias evaluation will consider up to second-order noise terms [6,52].

4.1. CRLB

The CRLB is given by [48]

$$error\left({\mathit{\beta }}^o\right)={\mathbf{FIM}}^{-1}\left({\mathit{\beta }}^o\right),$$

(38)

where ${\mathbf{FIM}}^{-1}\left({\mathit{\beta }}^o\right)$ is the Fisher Information Matrix (FIM),

$${\mathbf{FIM}}^{-1}\left({\mathit{\beta }}^o\right)=\frac{\partial {\mathit{\alpha }}^{oT}}{\partial {\mathit{\beta }}^o}{\mathbf{Q}}_{\alpha }^{-1}\frac{\partial {\mathit{\alpha }}^o}{\partial {\mathit{\beta }}^{oT}},$$

(39)

and

$$\frac{\partial {\mathit{\alpha }}^o}{\partial {\mathit{\beta }}^{oT}}=\left[\begin{array}{cc}{\displaystyle \frac{\partial \mathbf{r}}{\partial {{\mathbf{u}}^o}^T}}& {\displaystyle \frac{\partial \mathbf{r}}{\partial {\dot{\mathbf{u}}}^{oT}}}\\ {\displaystyle \frac{\partial \dot{\mathbf{r}}}{\partial {{\mathbf{u}}^o}^T}}& {\displaystyle \frac{\partial \dot{\mathbf{r}}}{\partial {\dot{\mathbf{u}}}^{oT}}}\end{array}\right],$$

(40)

where the derivatives are similar to (18)–(22) by replacing ${\mathit{\beta }}^{\left(g\right)}$ with ${\mathit{\beta }}^o$.

4.2. Bias and Covariance

The performance of the proposed solution is ultimately determined by the GN iteration. Therefore, attention will be directed towards the solution provided in Section 3.1. The bias and covariance of a column vector estimation derived from minimizing a cost function are summarized in [53]. This approach will be followed to evaluate the performance of the proposed MLE.

The estimate of $\mathit{\beta }$ implies

$$\mathbf{0}={\left.\frac{\partial J\left(\mathit{\beta }\right)}{\partial \mathit{\beta }}\right|}_{\mathit{\beta }=\widehat{\mathit{\beta }}},$$

(41)

where $\widehat{\mathit{\beta }}$ is the estimated value. Taking the Taylor-series expansion of the right side around the true value ${\mathit{\beta }}^o$ yields

$$\mathbf{0}\approx {\left.\frac{\partial J\left(\mathit{\beta }\right)}{\partial \mathit{\beta }}\right|}_{\mathit{\beta }={\mathit{\beta }}^o}+{\left.\frac{{\partial }^{2}J\left(\mathit{\beta }\right)}{\partial \mathit{\beta }\partial {\mathit{\beta }}^T}\right|}_{\mathit{\beta }={\mathit{\beta }}^o}(\widehat{\mathit{\beta }}-{\mathit{\beta }}^o).$$

(42)

The gradient vector and Hessian matrix evaluated at the true value are

$$\begin{array}{cc}\hfill ▽\mathcal{J}\left(\mathit{\beta }\right)=& {\left.\frac{\partial J\left(\mathit{\beta }\right)}{\partial \mathit{\beta }}\right|}_{\mathit{\beta }={\mathit{\beta }}^o}=2{\left(\frac{\partial {\mathit{\alpha }}^o}{\partial {\mathit{\beta }}^{oT}}\right)}^T{\mathbf{Q}}_{\alpha }^{-1}(\mathit{\alpha }-{\mathit{\alpha }}^o),\hfill \end{array}$$

(43)

$$\begin{array}{cc}\hfill \mathcal{H}\left(\mathit{\beta }\right)=& {\left.\frac{{\partial }^{2}J\left(\mathit{\beta }\right)}{\partial \mathit{\beta }\partial {\mathit{\beta }}^T}\right|}_{\mathit{\beta }={\mathit{\beta }}^o}=2{\left(\frac{{\partial }^2{\mathit{\alpha }}^o}{\partial {\mathit{\beta }}^o\partial {\mathit{\beta }}^{oT}}\right)}^T{\mathbf{Q}}_{\alpha }^{-1}(\mathit{\alpha }-{\mathit{\alpha }}^o)-2{\left(\frac{\partial {\mathit{\alpha }}^o}{\partial {\mathit{\beta }}^{oT}}\right)}^T{\mathbf{Q}}_{\alpha }^{-1}\frac{\partial {\mathit{\alpha }}^o}{\partial {\mathit{\beta }}^{oT}}.\hfill \end{array}$$

(44)

Using (43) and (44) to (42) yields

$$\widehat{\mathit{\beta }}-{\mathit{\beta }}^o\approx -▽\mathcal{J}\left(\mathit{\beta }\right)\mathcal{H}{\left(\mathit{\beta }\right)}^{-1}$$

(45)

Taking the expectation on both sides, we arrive at the bias

$$E[\widehat{\mathit{\beta }}-{\mathit{\beta }}^o]\approx -E[▽\mathcal{J}\left(\mathit{\beta }\right)]{\{E\left[\mathcal{H}\left(\mathit{\beta }\right)\right]\}}^{-1}=\mathbf{0}.$$

(46)

And the covariance matrix is

$$\begin{array}{cc}\hfill & E\left[(\widehat{\mathit{\beta }}-{\mathit{\beta }}^o){(\widehat{\mathit{\beta }}-{\mathit{\beta }}^o)}^T\right]\hfill \\ \hfill & \approx {\{E\left[\mathcal{H}\left(\mathit{\beta }\right)\right]\}}^{-1}E[▽\mathcal{J}\left(\mathit{\beta }\right)▽\mathcal{J}{\left(\mathit{\beta }\right)}^T]{\{E\left[\mathcal{H}\left(\mathit{\beta }\right)\right]\}}^{-1}.\hfill \end{array}$$

(47)

Putting (43) and (44) to the right side of (47) verifies that

$$E\left[(\widehat{\mathit{\beta }}-{\mathit{\beta }}^o){(\widehat{\mathit{\beta }}-{\mathit{\beta }}^o)}^T\right]\approx \mathbf{CRLB}\left({\mathit{\beta }}^o\right)$$

(48)

holds when the small noise conditions

$$\begin{array}{cc}\hfill & n_{i,k}/r_{i,k}\approx 0,e_{i,k}/r_{i,k}\approx 0,\hfill \end{array}$$

(49)

$$\begin{array}{cc}\hfill & \mathrm{or}{n}_{i,k}\approx 0,e_{i,k}\approx 0,\hfill \end{array}$$

(50)

are satisfied.

4.3. Complexity

This section analyzes the complexity of the proposed methods. The worst-case complexity of solving an SDP, as mentioned in [54], can be expressed as:

$$\sqrt{\mu }\left(m\sum _{i=1}^{N_{soc}}{\left(n_i^{soc}\right)}^2+m^2\sum _{i=1}^{N_{sd}}{\left(n_i^{sd}\right)}^2+m\sum _{i=1}^{N_{sd}}{\left(n_i^{sd}\right)}^3+m^3\right)\times ln(1/ϵ)$$

(51)

where m is the number of variables, $N_{soc}$(resp. $N_{sd}$) is the number of second-order cone (resp. semidefinite cone) constraints, and $n_i^{soc}$ (resp. $n_i^{sd}$) is the dimension of the i th second-order cone (resp. semidefinite cone),

$$\mu =\sum _{i=1}^{N_{sd}}{n}_i^{sd}+2N_{soc}$$

(52)

is the so-called barrier parameter that measures the geometric complexities of the cones involved, and $ϵ>0$ is the solution precision. The proposed SDP algorithm (37) has ${\left(2MK\right)}^2/2+\left(MK\right)+2N+3$ variables, $4.5KM+K^2{M}^2/2$ semidefinite cone constraints of size 1 ($2MK$ linear equality constraints in (28) and (31); $M^2{K}^2/2+KM/2$ linear inequality constraints in (29) and (33)); and one semidefinite cone constraint of size $2MK+1$ in (34), one semidefinite cone constraint of size $N+2$ in (35), and $MK$ second-order cone constraints of size N in (36). Since $MK$ is significantly larger than N in practice, the worst-case complexity is in the order of $O\left(\left(\sqrt{4.5K^2{M}^2+N^2}\right)\left(M^6{K}^6\right)(ln(1/ϵ))\right)$.

Regarding the computational complexity of the GN iteration, updating ${\mathit{\alpha }}^{\left(g\right)}$ in (17) requires $K(2M+4N+5)$ additions, $K(2N+9)$ multiplications, and K square root operations. Computing ${\mathbf{G}}^{\left(g\right)}$ needs $K(M-1)(12N+5)$ additions and $K(M-1)(15N+13)$ multiplications, putting them into the (17) iteration and obtaining a new ${\mathit{\beta }}^{(k+1)}$ requires $16N{K}^2{M}^2-32N{K}^{2}M+16N{K}^2-8NKM+8NK+16N^{2}KM-8K{N}^2-4N^2-8N^{2}K+2KM-2K$ additions and $16N{K}^2{M}^2-32N{K}^{2}M+16N^{2}KM-16N^{2}K+4NKM-4NK$ multiplications. The computational load of inverting ${\mathbf{G}}^{\left(g\right)T}\mathbf{W}{\mathbf{G}}^{\left(g\right)}$ is $O\left(N^3\right)$. As a result, the amount of computation in each GN iteration is $O\left(N{K}^2{M}^2\right)$.

Following the similar analysis above, the computational loads of the QN and LM methods are also $O\left(N{K}^2{M}^2\right)$ and are dominated by updating the Jacobian and the correction amount. The QN needs to approximate Hessian by an additional iterative optimization, e.g., BFGS [55], while the LM method introduces a term with a damping factor that will be updated in each iteration. They both require more computational load. The processing times tabulated in Table 5 validate the above analysis.

5. Simulation

Simulations were utilized to examine the performance of the proposed solution. The proposed SDP initial solution and the optimal GN solution are denoted by “SDP” and “GN” in the following figures. The TDOA and FDOA measurement noises are assumed to be independent Gaussian random variables, so their covariance matrices are set to $\mathbf{Q}={\sigma }^2\left({\mathbf{I}}_{M-1}+{\mathbf{1}}_{M-1}{\mathbf{1}}_{M-1}^T\right)/2$ and $\dot{\mathbf{Q}}=0.1\mathbf{Q}$, where ${\sigma }^2$ is the measurement noise level in ${\mathrm{m}}^2$. The observation interval $\tau$ is fixed at 1 s. The number of ensemble runs is $T=200$. The proposed SDP algorithm is implemented using the CVX toolbox with the solver SeDuMi, where the precision is set to the best. The performance is evaluated by the root mean-squared error (RMSE), defined as

$$\left\{\begin{array}{cc}RMSE\left(\widehat{\mathbf{u}}\right)\hfill & =\sqrt{\frac{1}{T}{\sum }_{t=1}^T{\parallel {\mathbf{u}}^o-\widehat{\mathbf{u}}\parallel }^2},\hfill \\ RMSE\left(\widehat{\dot{\mathbf{u}}}\right)\hfill & =\sqrt{\frac{1}{T}{\sum }_{t=1}^T{\parallel {\dot{\mathbf{u}}}^o-\widehat{\dot{\mathbf{u}}}\parallel }^2},\hfill \end{array}\right.$$

(53)

where $\widehat{\mathbf{u}}$ and $\widehat{\dot{\mathbf{u}}}$ are the estimates at the t-th test. The values of ${\eta }_1$ and ${\eta }_2$ are selected from $\{{10}^{-2},{10}^{-3},\cdots ,{10}^{-10}\}$ and $\{1,{10}^{-1},\dots ,{10}^{-8}\}$ through the strategy elaborated at the end of Section 3.2.

The proposed method is applicable for both 2D and 3D cases. If the sensors and source are unmanned ground vehicles (UGVs) that move on a plane, the localization scenario is 2D, while for unmanned aerial vehicles (UAVs) or underwater moving sensors and source deployed with different depths, the localization scenario is 3D. Thus, separate simulations were conducted to evaluate the algorithm’s performance in both 2D and 3D scenarios. Of the large number of papers investigated, only a few studied TMA using TDOA or FDOA [45,46]. However, they are neither applicable to hybrid measurements nor the same scenario or problem. Thus, the proposed SDP and GN methods are compared with the CRB, the optimal theoretical lower bound for evaluating the performance, to validate the effectiveness and feasibility.

5.1. 3D Scenario

Let us begin by testing the 3D case. As shown in Figure 2, the scenario comprises 3 sensors utilized to locate a moving source, where the noise varies from ${10}^{-1}$${\mathrm{m}}^2$ to ${10}^3$${\mathrm{m}}^2$. The initial positions and velocities of the sensors are listed in

Table 3. Positions and velocities of the sensors in the 3D scenario.

Sensor No.	$\mathbf{x} \left(\mathbf{m}\right)$	$\mathbf{y} \left(\mathbf{m}\right)$	$\mathbf{z} \left(\mathbf{m}\right)$	$\dot{\mathbf{x}} (\mathbf{m}/\mathbf{s})$	$\dot{\mathbf{y}} (\mathbf{m}/\mathbf{s})$	$\dot{\mathbf{z}} (\mathbf{m}/\mathbf{s})$
1	300	100	150	30	−20	20
2	400	150	100	−30	10	20
3	300	500	200	10	−20	10

. The source is initially observed from $\mathbf{u}={[285,325,275]}^T$ m with a velocity of $\dot{\mathbf{u}}=[$ ${-20,15,40]}^T$ m/s. The number of observing times is $K=16$. The RMSEs of the SDP and GN methods and the corresponding CRLB are shown in Figure 3. The GN method achieves a CRLB accuracy if the noise is not larger than ${10}^2$${\mathrm{m}}^2$. However, when the noise is at ${10}^3$${\mathrm{m}}^2$, the GN method encounters divergence problems since the estimation results from SDP are too far from the true values. The RMSE of these estimations using the SDP method fails to attain the CRLB since the effectiveness of the penalty terms is limited when only three sensors exist. This limitation hinders the transformed problem from accurately approximating the original problem and renders the SDP method suboptimal.

Considering a practical unmanned aerial vehicle (UAV) object, the aim is to localize the UAV by intercepting its communication and control signal at 2.4 GHz with 10 MHz bandwidth. The observation time is ${10}^{-4}$ s, and the SNR is 0 dB. Following the CRLB of TDOA estimation given by Equation (57) in [56], the theoretical noise power in TDOA is approximately 0.1 ${\mathrm{m}}^2$. Thus, the noise is reasonably fixed at ${\sigma }^2={10}^{-1}$${\mathrm{m}}^2$ in following simulations.

In the second numerical experiment, the RMSE performance is assessed as the number of observing times K increases over the range $[2,16]$. The noise is fixed at ${\sigma }^2={10}^{-1}$${\mathrm{m}}^2$, and other parameters remain consistent with those in Figure 3. As depicted in Figure 4, the performance of the proposed method improves with larger K. This enhancement stems from the relative motion between the source and sensors, which increases the effective baseline length and accumulates more information about the source’s position and velocity. Even with only two observations of the source, the GN algorithm remains feasible, achieving an RMSE at the CRLB level. This feasibility arises from the establishment of 8 equations to solve for 6 unknowns.

Let us further examine the RMSE of the proposed SDP and GN methods as the source range increases. The initial position of the source is oriented with azimuth and elevation $\theta =48.75$ deg and $\varphi =32.46$ deg relative to the coordinate origin. The source’s initial position is calculated as

$$\mathbf{u}=d·{[cos\theta cos\varphi ,sin\theta cos\varphi ,sin\varphi ]}^T,$$

(54)

where d represents the distance between the source’s initial position and the origin. Other settings are the same as above. The results are shown in Figure 5. The performance trends of the proposed algorithms align with the variations in the CRLB. As expected, the GN method can closely match the CRLB regardless of whether the source is near or far. Although SDP fails to reach the CRLB, it closely approximates the GN method in estimating the source’s position. Moreover, it does not encounter the divergence problem, even when the source range increases to 7000 m. While the SDP yields an RMSE approximately 5 dB higher than the CRLB, it still ensures that the GN method converges to the optimal estimate.

Figure 6 illustrates the RMSE results as the source velocity increases. The noise power is set to ${\sigma }^2={10}^{-1}$ ${\mathrm{m}}^2$, with other parameters remaining consistent with those in Figure 3. For clarity in displaying the performance as the velocity magnitude changes, the source velocity is defined as

$$\dot{\mathbf{u}}=vs.·{[cos\tilde{\theta }cos\tilde{\varphi },sin\tilde{\theta }cos\tilde{\varphi },sin\tilde{\varphi }]}^T,$$

(55)

where $\tilde{\theta }=143$ deg and $\tilde{\varphi }=58$ deg represent the direction of the source velocity. The SDP method consistently provides a reliable initial estimate as the source velocity increases, thereby ensuring the optimal performance of the GN algorithm in approaching the CRLB.

5.2. 2D Scenario

The proposed SDP and GN methods are not limited to 3D scenarios; they are also applicable to 2D cases. In a 2D scenario, the SDP and GN methods require at least 2 sensors, with the initial positions and velocities listed in

Table 4. Positions and velocities of the sensors in the 2D scenario.

Sensor No.	$\mathbf{x} \left(\mathbf{m}\right)$	$\mathbf{y} \left(\mathbf{m}\right)$	$\dot{\mathbf{x}} (\mathbf{m}/\mathbf{s})$	$\dot{\mathbf{y}} (\mathbf{m}/\mathbf{s})$
1	400	150	−50	−30
2	150	−100	−20	20

. The geometric scenario is depicted in Figure 7.

Similar to the simulations in the 3D scenario, this section evaluates the performance for the 2D case as the noise power increases. The initial position and velocity of the source are set to $\mathbf{u}={[300,200]}^T$ m and $\dot{\mathbf{u}}={[20,-15]}^T$ m/s, respectively. As illustrated in Figure 8, the GN algorithm can achieve CRLB performance when the noise power does not exceed 10 ${\mathrm{m}}^2$, while it encounters divergence issues when the noise power reaches ${10}^3$ ${\mathrm{m}}^2$. SDP provides 3 dB worse RMSE compared to the GN method, which is sufficient as an initialization.

The simulation also examines the RMSE of both the SDP and GN methods for the 2D scenario with respect to the number of observation times. As shown in Figure 9, the GN algorithm achieves CRLB performance as the number of observation times increases. The noise power is set to ${\sigma }^2={10}^{-1}$ ${\mathrm{m}}^2$. It requires at least three observations to establish 6 equations, which is overdetermined for solving 4 unknowns.

Figure 10 illustrates the RMSEs of the proposed solutions as the source range increases in 2D. Here, the arrival angle with respect to the origin is set to $\theta =80$ deg. Thus, the source’s initial position is denoted by

$$\mathbf{u}=d·{[cos\theta ,sin\theta ]}^T,$$

(56)

where d is the distance between the source’s initial position and the origin. The noise power is set to ${\sigma }^2={10}^{-3}$ ${\mathrm{m}}^2$, and $\tau$ is set to 4 s. SDP is not as effective as in the 3D scenario. However, fortunately, it provides an initialization for the GN algorithm, which can also achieve the CRLB when the source is within 2100 m. It slightly deviates from the CRLB when the source’s initial position is farther away.

Figure 11 depicts the results as the source velocity increases. The noise power is specified as ${\sigma }^2={10}^{-1}$ ${\mathrm{m}}^2$. The source velocity is expressed as

$$\dot{\mathbf{u}}=d·[cos\tilde{\theta },sin\tilde{\theta }],$$

(57)

where $\tilde{\theta }=$$-36.87$ deg denotes the direction of the source’s motion. The GN algorithm consistently achieves CRLB performance, while the SDP method provides a suboptimal estimate. However, it is sufficiently close to the true value such that the GN method can converge.

5.3. Geometric Dilution of Precision

The GDOP contour for the 2D scenario is presented to evaluate the theoretical performance of the proposed method. The sensor configuration and the velocity of the source are similar to those depicted in Figure 7. As illustrated in Figure 12, in the 2D scenario, the localization accuracy remains high when the source is nearby but deteriorates significantly as the source moves farther away.

5.4. Convergence and Complexity Comparisons

Section 1 has clarified that the GN method advances the QN and LM methods since it balances the performance and efficiency. Considering that the Newton–Raphson method is designed for a single unknown estimation problem only, the proposed GN method will be compared with the QN and LM methods that are implemented by the MATLAB built-in functions fminunc and lsqnonlin. The convergence curves and processing times are shown in Figure 13 and

Table 5. Processing times.

Method	GN	LM	QN
Time (ms)	0.615	3.3	8.926

.

It is demonstrated in Figure 13 that the GN and LM methods converge to the CRB after only two iterations on average. Additionally, these two methods exhibit similar performance since their convergence curves coincide. However, the QN method requires 26 iterations to converge. Table 5 demonstrates that the processing times of the proposed GN method are significantly lower than those of the LM and QN methods.

6. Conclusions

The paper investigates moving source localization utilizing TDOA and FDOA measurements, particularly addressing scenarios where the number of sensors is limited to the dimensionality of the problem. It begins by formulating the MLE problem and proposes a two-stage approach for optimal solution. The first stage involves an SDP formulation via SDR, yielding a suboptimal solution utilized as an initialization for the GN iteration in the second stage. Simulations demonstrate that the SDP provides an initial estimate close to the true value, enabling the GN method to achieve CRLB performance. The proposed MLE is shown to be effective with a minimum of 3 sensors in 3D and 2 sensors in 2D, requiring at least 2 and 3 observation times, respectively. The method opens possibilities for simultaneous estimation of position and velocity with few sensors, potentially benefiting applications such as UAV surveillance, vehicle localization, and person tracking. The proposed method of this paper includes an SDP stage that demands high complexity. One future research subject is improving efficiency and developing closed-form solutions. The other limitation is the study assumes the sensor positions and velocities are accurately known, and the measurements do not contain outliers. How to handle sensor uncertainties and address outliers in measurements is perhaps a valuable research field. Moreover, the deployment and movement of sensors have a significant effect on source localization performance, making it an interesting subject for future research.