Signal Source Localization with Long-Term Observations in Distributed Angle-Only Sensors

Abstract

Angle-only sensors cannot provide range information of targets and in order to determine accurate position of a signal source, one can connect distributed passive sensors with communication links and implement a fusion algorithm to estimate target position. To measure moving targets with sensors on moving platforms, most of existing algorithms resort to the filtering method. In this paper, we present two fusion algorithms to estimate both the position and velocity of moving target with distributed angle-only sensors in motion. The first algorithm is termed as the gross least square (LS) algorithm, which takes all observations from distributed sensors together to form an estimate of the position and velocity and thus needs a huge communication cost and a huge computation cost. The second algorithm is termed as the linear LS algorithm, which approximates locations of sensors, locations of targets, and angle-only measures for each sensor by linear models and thus does not need each local sensors to transmit raw data of angle-only observations, resulting in a lower communication cost between sensors and then a lower computation cost at the fusion center. Based on the second algorithm, a truncated LS algorithm, which estimates the target velocity through an average operation, is also presented. Numerical results indicate that the gross LS algorithm, without linear approximation operation, often benefits from more observations, whereas the linear LS algorithm and the truncated LS algorithm, both bear lower communication and computation costs, may endure performance loss if the observations are collected in a long period such that the linear approximation model becomes mismatch.

1. Introduction

For some passive sensors, such as infrared sensors, photoelectric sensors and cameras, they can detect targets by receiving electromagnetic signals. As they do not emit signals, they can probe targets in a stealth manner [1]. However, such sensors can measure only angles of signals and thus are termed as angle-only sensors subsequently. The signal position information, which is of great concern in many situations, cannot be obtained with a sensor. To determine the position of signal sources, one can connect distributed sensors with communication links and then estimate the position through a fusion algorithm. This is a hot topic in recent years and gains wide attentions of scholars in different fields [2,3,4,5].

In the 3-dimensional (3D) scenario, the angle information measured by each passive sensor includes the azimuth and elevation of the signal. From a mathematical perspective, each angle observation can be represented by a straight line passing the sensor and a target in space. If no error occurs in this process, all the lines will intersect in a point in space, which is the location of a signal source. In practice, both sensor location measures and angle measures are inevitably contaminated by measurement noises and then the lines may not intersect a point in space. However, as if the signals are from the same target, the lines will intersect in a small volume, whose center can be deemed as the location of a target. Following this concept, an angle-only positioning algorithm is presented in [6] and a closed-form solution is derived.

In real applications, if the targets of interest are static, or if the sampling frequency to the signals is too high in contrast to the velocities of possible targets, the algorithms can be developed under an assumption that the velocity of the target is static. The least squares (LS) algorithm is applied in the target position estimation based on angle-only measurements by linearizing the angle observation equations [7,8,9,10]. The intersection localization algorithm is obtained by considering that the straight lines formed by the angle observations will intersect in a small volume in space [11,12,13]. Real sensor often makes observations in an asynchronous manner, namely the observations are not obtained at the same instants. The stationary target assumption will also make the sensor synchronization problem easier, because we can totally drop the timing information of the observations. If the target is stationary, even if the sensors are moving, the straight lines formed by the angle measurements of multiple sensors at different times will converge to a small area near the target location. Therefore, in this scenario, one just needs to solve a target positioning algorithm in an asynchronous manner.

Once the target motion should be considered at different observation instants, target location estimation will face greater biases and then one has to take the target motion issue into consideration. Meanwhile, for moving targets, the observation instants should be taken into account and then as the distributed fusion algorithm should take instant information into account, the fusion algorithm becomes more complicated. There are mainly two strategies available so far. The first strategy is to use the filtering algorithms, such as the Kalman filter that can estimate target velocity through observations from different instants. In the target tracking theoretical framework, the angle-only observations can be described by a measurement equation, although it is heavily nonlinear. Therefore, a nonlinear filtering algorithm should be used [14,15]. For instance, the extended Kalman filtering (EKF) algorithm linearizes the angle measurement function through the first-order Taylor approximation, and then uses the standard Kalman filtering algorithm for the angle-only target tracking problem [16]. The cubature Kalman filter (CKF) [17], the unscented Kalman filter (UKF) [18], the pseudo linear Kalman filter (PLKF) [19,20], the particle filter [21] and a series of sigma-point based algorithms can also be used in the target tracking problem with distributed angle-only sensors.

Although the tracking algorithms have been widely used to estimate moving target positions, it requires the noise distribution parameters known a priori. It also faces the convergence problem if the initial state is set improperly [22,23]. In distributed sensor networks, if each observation undergoes a tracking process, the computation cost will also be high since the data amount of observations are often intensive in practice. Therefore, a good positioning algorithm should be implemented before filtering. For instance, in [15], short-term angle-only observations are fused by a distributed positioning algorithm, whose outputs are then processed by a tracking algorithm.

In the other strategy, the target position and velocity can be estimated together and then the result is valid in a longer period. In this case, the tracking operation can be performed in a longer period, so that the computation cost can be further reduced. However, if the velocity is estimated, more optimization variables are involved and then the optimization problem is more complicated. Meanwhile, in a distributed sensor configuration, the communication cost between the sensors may be high if all the observations are transmitted to a fusion center. In this paper, we study the distributed positioning of moving targets with distributed asynchronous angle-only sensors. We consider the scenario where multiple asynchronous passive sensors are linked with the fusion center through communication links. First, we formulate an algorithm, termed as gross LS algorithm, that takes all angle observations of multiple sensors together with their positions in certain period to estimate the position and velocity of the target. Different observations contribute different lines and with many lines available, both the target position and its velocity can be estimated. The classical LS algorithm is formulated such that the computation cost is reduced a lot.

Due to the huge amount of data, this algorithm still has high computational complexity and high communication cost. In order to reduce the communication and computation cost, we further present a distributed positioning algorithm, termed as linear LS algorithm, that can implement the fusion algorithm in a parallel computation manner. In detail, both positions and angle observations of local sensors are processed by LS operations, whose outputs are zero and first orders of the Taylor series of corresponding parameters. The outputs are then transmitted to a fusion center for which we derive a fusion algorithm to efficiently combine position and velocity estimates for a higher parameter estimation accuracy. The later algorithm can greatly compress the data rate from local sensors to the fusion center, such that the communication cost is greatly reduced. Meanwhile, local observations are represented by a few parameters and thus the fusion algorithm also needs a lower computation cost. The sensor location can be recorded asynchronously with the angle observations and thus can make the algorithm easier in applications. Meanwhile, a truncated LS algorithm, which replaces the velocity estimation of the linear LS algorithm by a simply average operation, is also presented.

Numerical results are obtained with distributed asynchronous angle-only sensors measuring a moving target with certain velocity. The convergence performance of both the algorithms are presented first, in order to examine the impact of the number of observations on the positioning performance. Then the impact of the linear approximation of position and angle measures on the estimation accuracy is analyzed. It will be found that the gross LS algorithm often benefits from more observations. However, although the linear LS algorithm and the truncated LS algorithm will perform good if the number of observations is small, as the number of observations increase, their performances will degrade, as a resulting of the linear model mismatching. The truncated LS algorithm will perform better in a short period than the linear LS algorithm but worse in a longer time. To an extreme, the estimation performance of the linear LS algorithm may deteriorate with more observations if the model mismatching is severe. We also verify that the linear LS algorithm has a lower communication cost in most situations and examine the performance loss due to inaccurate platform velocity estimates. Numerical angle distortion errors under the linear approximations are also analyzed.

2. Localization with Angle-Only Passive Sensors

2.1. Signal Model of Passive Observations

Consider a passive sensor network with N widely separated sensors and M targets in the surveillance volume. All the N passive sensors can measure only direction of arrival (DOA) of signals, based on which real position of a signal emitter can be estimated. Assume that all the sensors operate in the same coordinate system through some inherent position and attitude measurement devices, such as the Global Positioning System (GPS) and inertial sensors. A typical coordinate system is the earth-centered earth-fixed (ECEF) of the World Geodetic System 84 (WGS84). Both the targets and the sensors are in motion by assumption. The real position of the nth sensor at instant t is denoted by ${\mathbf{p}}_n^{\mathrm{o}}\left(t\right)={[x_{n,s}^{\mathrm{o}}\left(t\right),y_{n,s}^{\mathrm{o}}\left(t\right),z_{n,s}^{\mathrm{o}}\left(t\right)]}^{\mathrm{T}}$, $n=1,2,\dots ,N$, where ${(·)}^{\mathrm{T}}$ denotes the transpose operation, and $x_{n,s}^{\mathrm{o}}\left(t\right),y_{n,s}^{\mathrm{o}}\left(t\right),z_{n,s}^{\mathrm{o}}\left(t\right)$ denote the $x,y,z$ coordinates of the nth sensor in the common coordinate system at instant t, respectively. The real position of the mth target at instant t is denoted by ${\mathbf{g}}_m^{\mathrm{o}}\left(t\right)={[x_{m,g}^{\mathrm{o}}\left(t\right),y_{m,g}^{\mathrm{o}}\left(t\right),z_{m,g}^{\mathrm{o}}\left(t\right)]}^{\mathrm{T}},m=1,\cdots ,M$, where $x_{m,g}^{\mathrm{o}}\left(t\right),y_{m,g}^{\mathrm{o}}\left(t\right),z_{m,g}^{\mathrm{o}}\left(t\right)$ denote the $x,y,z$ coordinates of the mth target at instant t, respectively. The topology of the passive sensors and targets are shown in Figure 1.

For the nth sensor, signals are detected and their DOAs are measured at instants denoted by $t_{k,n},k=1,\cdots ,N_n$, where $N_n$ denotes the number of observations of the nth sensor. At the instant $t_{k,n}$, assume that the position of the nth sensors is measured as

$${\mathbf{p}}_{k,n}\left(t_{k,n}\right)={\mathbf{p}}_n^{\mathrm{o}}\left(t_{k,n}\right)+\mathsf{\Delta }{\mathbf{p}}_n\left(t_{k,n}\right)={[x_{n,s}\left(t_{k,n}\right),y_{n,s}\left(t_{k,n}\right),z_{n,s}\left(t_{k,n}\right)]}^{\mathrm{T}},k=1,\cdots ,N_n$$

(1)

where $\mathsf{\Delta }{\mathbf{p}}_n\left(t_{k,n}\right)$ denotes the sensor self-positioning error. For simplicity, we assume that the sensor self-positioning error follows zero mean Gaussian distributions with covariance matrices ${\mathbf{C}}_{\mathrm{s}}(k,n)=\mathbb{E}\left(\mathsf{\Delta }{\mathbf{p}}_n\left(t_{k,n}\right)\mathsf{\Delta }{\mathbf{p}}_n^{\mathrm{T}}\left(t_{k,n}\right)\right)$, where $\mathbb{E}$ denotes the expectation operation.

At instant $t_{k,n}$, the real position of the mth signal source is denoted by

$${\mathbf{g}}_m^{\mathrm{o}}\left(t_{k,n}\right)={[x_{m,g}^{\mathrm{o}}\left(t_{k,n}\right),y_{m,g}^{\mathrm{o}}\left(t_{k,n}\right),z_{m,g}^{\mathrm{o}}\left(t_{k,n}\right)]}^{\mathrm{T}},m=1,\cdots ,M.$$

(2)

Assume that all the observations regarding the same target are obtained in a short period $\mathcal{T}=[T_1,T_2]$. In this period, assume that the location of the mth target can be expressed by

$${\mathbf{g}}_m^{\mathrm{o}}\left(t\right)\approx {\mathbf{g}}_m^{\mathrm{o}}\left(t_0\right)+{\mathbf{v}}_m^{\mathrm{o}}(t-t_0)$$

(3)

where ${\mathbf{g}}_m^{\mathrm{o}}\left(t_0\right)$ denotes the location of the target at the reference instant $t_0$, and ${\mathbf{v}}_m^{\mathrm{o}}$ denotes the velocity over $\mathcal{T}$. The signal model in use depends on the velocity of the target and the period of observations. If all the observations are collected in a short period and the velocity is small, then one can simply assume ${\mathbf{g}}_m^{\mathrm{o}}\left(t\right)\approx {\mathbf{g}}_m^{\mathrm{o}}\left(t_0\right)$ as [6]. Under the signal model (3), more observations can be used to make an estimation of the target space locations. If the observations are obtained in a long period and the velocity is huge, then this model may also mismatch and higher order approximations may be used.

For the nth angle-only passive sensors, the lth observation at $t_{k,n}$ is indexed by a triple $(l,k,n),l=1,\cdots ,L_{k,n}$. For simplicity, we also encode all the triples available, corresponding to all the observations available, with a one-to-one function $\Omega (l,k,n)\to i$. Then we define a set ${\mathcal{L}}_{k,n}$ by

$${\mathcal{L}}_{k,n}=\left\{i\right|i=(l,k,n),l=1,\cdots ,L_{k,n}\},k=1,\cdots ,N_n,n=1,\cdots ,N$$

(4)

denotes a set of signal indices detected at the instant $t_{k,n}$ by the nth sensor. Therefore, $L_{k,n}=\left|{\mathcal{L}}_{k,n}\right|$, where $|·|$ over a set denotes the cardinality of the set. As the possibility of miss detection, false alarms and overlapping of signal sources, $L_{k,n}$ may not be equal to M. Denote

$${\mathcal{L}}_n={\cup }_{k=1}^{N_n}{\mathcal{L}}_{k,n},\mathcal{L}={\cup }_{n=1}^N{\mathcal{L}}_n$$

(5)

where ∪ denotes the union operation. The total number of observations by N sensors is denoted by

$$N_{\mathrm{s}}=\left|\mathcal{L}\right|=\sum _{n=1}^N\sum _{k=1}^{N_n}{L}_{k,n}.$$

(6)

Each observation is associated with one of M targets or the false alarm indexed by 0, represented by a set $\mathcal{M}=\{0,1,\cdots ,M\}$. It can be considered as a mapping $\psi :\mathcal{L}\to \mathcal{M}$, which is a correct mapping and is thus typically unknown in practice. According to our setting, the index set $\mathcal{L}$ can be partitioned into $M+1$ disjoint sets ${\mathcal{A}}_0,{\mathcal{A}}_1,\cdots ,{\mathcal{A}}_M$, and ${\mathcal{A}}_m$ is defined by

$${\mathcal{A}}_m=\left\{i\right|\psi \left(i\right)=m,i\in \mathcal{L}\},$$

(7)

where ${\mathcal{A}}_0$ denotes the index of observations corresponding to false alarms, and ${\mathcal{A}}_m$ denotes the index set of observations from the mth signal source. As a partition of $\mathcal{A}$, we have ${\mathcal{A}}_i\cap {\mathcal{A}}_j=\varnothing$, $i,j\in \mathcal{M}$, $i\ne j$, and $\mathcal{A}={\cup }_{i=0}^M{\mathcal{A}}_i$, where ∩ denotes the intersection operation of sets. Assume that $|{\mathcal{A}}_m|=L_{\mathrm{m}}$ and there are totally $L={\sum }_{m=0}^M{L}_{\mathrm{m}}$ observations available.

The signal indices in ${\mathcal{A}}_m$ are composed of signal indices from all the sensors and the sub set for the nth sensor is denoted by

$${\mathcal{A}}_{n,m}={\mathcal{L}}_n\cap {\mathcal{A}}_m$$

(8)

which indicates the observations from the nth sensor probing the mth target. Denote $K_{n,m}=\left|{\mathcal{A}}_{n,m}\right|$ and then we have $L_{\mathrm{m}}={\sum }_{n=1}^N{K}_{n,m}$ and ${\mathcal{A}}_m=\cup {\mathcal{A}}_{n,m}$.

For simplicity, we first assume that the mapping $\psi$ is exactly known and then observations associated with ${\mathcal{A}}_m$ is exactly known. For observation $i\in {\mathcal{A}}_m$, real azimuth angle and elevation angle, regarding the nth sensor at $t_{k,n}$, can be expressed by

$$\begin{array}{cc}\hfill {\theta }_i^{\mathrm{o}}& ={tan}^{-1}\left(y_{m,g}^{\mathrm{o}}\left(t_{k,n}\right)-y_{n,s}^{\mathrm{o}}\left(t_{k,n}\right),x_{m,g}^{\mathrm{o}}\left(t_{k,n}\right)-x_{n,s}^{\mathrm{o}}\left(t_{k,n}\right)\right)\hfill \\ \hfill {\phi }_i^{\mathrm{o}}& =arctan\frac{z_{m,g}^{\mathrm{o}}\left(t_{k,n}\right)-z_{n,s}^{\mathrm{o}}\left(t_{k,n}\right)}{\sqrt{{(x_{m,g}^{\mathrm{o}}\left(t_{k,n}\right)-x_{n,s}^{\mathrm{o}}\left(t_{k,n}\right))}^2+{(y_{m,g}^{\mathrm{o}}\left(t_{k,n}\right)-y_{n,s}^{\mathrm{o}}\left(t_{k,n}\right))}^2}}\hfill \end{array}$$

(9)

respectively, where ${\theta }_i^{\mathrm{o}}\in (-\pi ,\pi )$, ${\phi }_i^{\mathrm{o}}\in (-\frac{\pi }{2},\frac{\pi }{2})$, ${tan}^{-1}(\ast )$ is called the two-argument inverse tangent function [24,25] and $arctan(\ast )$ is the inverse tangent function. Denote ${\mathit{\eta }}_i^{\mathrm{o}}={[{\theta }_i^{\mathrm{o}},{\phi }_i^{\mathrm{o}}]}^{\mathrm{T}}$. The azimuth angle and elevation angle measures can be written as

$$\begin{array}{cc}\hfill {\mathit{\eta }}_i& ={[{\theta }_i,{\phi }_i]}^{\mathrm{T}}={\mathit{\eta }}_i^{\mathrm{o}}+\mathsf{\Delta }{\mathit{\eta }}_i\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill {\theta }_i& ={\theta }_i^{\mathrm{o}}+\mathsf{\Delta }{\theta }_i\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill {\phi }_i& ={\phi }_i^{\mathrm{o}}+\mathsf{\Delta }{\phi }_i\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill \mathsf{\Delta }{\mathit{\eta }}_i& ={[\mathsf{\Delta }{\theta }_i,\mathsf{\Delta }{\phi }_i]}^{\mathrm{T}}\hfill \end{array}$$

(13)

where $\mathsf{\Delta }{\theta }_i$ and $\mathsf{\Delta }{\phi }_i$ represent the measurement noise of the azimuth angle and elevation angle, respectively.

For simplicity, we assume that observation noises $\mathsf{\Delta }{\theta }_i$ and $\mathsf{\Delta }{\phi }_i$ are statistically independent and follow zero-mean Gaussian distribution. The covariance matrices of $\mathsf{\Delta }\mathit{\eta }$ are denoted by ${\mathbf{C}}_{\eta }\left(i\right)=\mathbb{E}\left(\mathsf{\Delta }\mathit{\eta }\mathsf{\Delta }{\mathit{\eta }}^{\mathrm{T}}\right)\in {\mathbb{C}}^{2\times 2}$, namely $\mathsf{\Delta }\mathit{\eta }\sim \mathcal{N}(\mathbf{0},{\mathbf{C}}_{\mathit{\eta }}\left(i\right))$, which is typically affected by the SNR of the signal, where $\mathcal{N}(\mathbf{0},{\mathbf{C}}_{\eta }\left(i\right))$ denotes the zero-mean Gaussian distribution with mean $\mathbf{0}$ and covariance matrix ${\mathbf{C}}_{\eta }\left(i\right)$.

2.2. Estimation of Target Track

Each angle-only observation contributes a line in 3D space and without measurement error, a target will be present at the line. With many angle-only observations, real position of the target can be determined. The line associated with the ith observation can be expressed by

$$\begin{array}{c}\hfill L_i:\mathbf{x}\left(t_i\right)={\mathbf{p}}_i+{\alpha }_i{\mathbf{e}}_i,{\alpha }_i\in \mathbb{R}\end{array}$$

(14)

where ${\mathbf{p}}_i$ denotes the sensor location regarding the ith observation, ${\alpha }_i$ is a parameter indicating the distance to the origin ${\mathbf{p}}_i$, ${\mathbf{e}}_i=\mathbf{e}\left({\mathit{\eta }}_i\right)={[e_{i,x},e_{i,y},e_{i,z}]}^{\mathrm{T}}\in {\mathbb{R}}^{3\times 1}$ is the normalized direction vector associated with the angle observation ${\mathit{\eta }}_i$, namely $\parallel {\mathbf{e}}_i\parallel =1$, $\parallel ·\parallel$ over a vector denotes the ${\ell }_2$-norm, and

$$\begin{array}{cc}\hfill e_{i,x}& =cos\left({\theta }_i\right)cos\left({\phi }_i\right)\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill e_{i,y}& =sin\left({\theta }_i\right)cos\left({\phi }_i\right)\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill e_{i,y}& =sin\left({\phi }_i\right).\hfill \end{array}$$

(17)

In what follows, we consider the observations in ${\mathcal{A}}_m,m\ne 0$. From (3), we can rewrite

$$\mathbf{x}\left(t_i\right)={\mathbf{p}}_i+{\alpha }_i{\mathbf{e}}_i={\mathbf{g}}_{0,m}+{\mathbf{v}}_{\mathrm{g},m}(t_i-t_0)+{\mathit{ϵ}}_{\mathrm{e}}\left(i\right)={\mathbf{A}}_{\mathrm{e}}\left(i\right){\mathbf{q}}_m+{\mathit{ϵ}}_{\mathrm{e}}\left(i\right)$$

(18)

where ${\mathbf{q}}_m={[{\mathbf{g}}_{0,m}^{\mathrm{T}},{\mathbf{v}}_{\mathrm{g},m}^{\mathrm{T}}]}^{\mathrm{T}},{\mathbf{g}}_{0,m}={\mathbf{g}}_m^{\mathrm{o}}\left(t_0\right),{\mathbf{v}}_{\mathrm{g},m}={\mathbf{v}}_{\mathrm{m}}^{\mathrm{o}}\left(t_0\right)$, ${\mathit{ϵ}}_{\mathrm{e}}\left(i\right)$ denotes the bias term,

$${\mathbf{A}}_{\mathrm{e}}\left(i\right)=[\mathbf{I},(t_i-t_0)\mathbf{I}]\in {\mathbb{R}}^{3\times 6}$$

(19)

and $\mathbf{I}$ denotes the identity matrix. In (18), there are totally 7 unknown parameters and an observation can provide 3 equations. In addition to an observation, one can obtain another 3 equations and the number of unknown parameters will increase by 1. Unless specified, we always refer to the mth target and drop the subscripts m in situations without ambiguity subsequently, e.g., denote ${\mathbf{g}}_{0,m}\to {\mathbf{g}}_0,{\mathbf{v}}_{\mathrm{g},m}\to {\mathbf{v}}_{\mathrm{g}},{\mathbf{q}}_m\to \mathbf{q}$.

In order to determine the location and velocity of the target, the optimization problem can be formulated as

$$\underset{\alpha ,\mathbf{q}}{min}\parallel \mathbf{P}+\mathbf{E}diag\left(\alpha \right)-{\mathbf{A}}_{\mathrm{h}}({\mathbf{1}}_6\otimes \mathbf{q})\parallel$$

(20)

where $\parallel ·\parallel$ refers to the ${\ell }_2—$norm subsequently unless explicitly specified, $\mathit{\alpha }={[{\alpha }_1,\cdots ,{\alpha }_{L_{\mathrm{m}}}]}^{\mathrm{T}}$, $\mathbf{P}$ is a matrix whose columns are ${\mathbf{p}}_i,i\in {\mathcal{A}}_m$, $\mathbf{P}=[{\mathbf{p}}_{,}\cdots ,{\mathbf{p}}_{L_{\mathrm{m}}}]$, ${\mathbf{1}}_6$ denotes a $6\times 1$ all-one vector of length $L_{\mathrm{m}}$, $\mathbf{E}$ is a matrix whose columns are ${\mathbf{e}}_i,i\in {\mathcal{A}}_m$, ⊗ denotes the Kronecker product operation,

$${\mathbf{A}}_{\mathrm{h}}=[{\mathbf{A}}_{\mathrm{e}}\left(i\right),\cdots ]\in {\mathbb{R}}^{3\times 6L_{\mathrm{m}}},i\in {\mathcal{A}}_m,$$

(21)

and $diag(·)$ with a vector entry denotes a diagonal matrix with the vector as diagonal elements.

2.3. The Gross LS Algorithm

In practice, the observations are generally contaminated by measurement noise and then the lines often do not intersect into a point in space. With $L_{\mathrm{m}}$ observations, there are totally $L_{\mathrm{m}}+6$ unknown parameters and $3L_{\mathrm{m}}$ equations. Therefore, if we have at least 3 observations from angle-only sensors, we can find 9 unknown parameters together. For that purpose, let $\overline{\mathbf{q}}={[{\mathbf{q}}^{\mathrm{T}},{\mathit{\alpha }}^{\mathrm{T}}]}^{\mathrm{T}}$ and then we can reformulate (18). In order to minimize the mismatch, the optimization problem can be rewritten as

$$\underset{\mathit{\alpha },\mathbf{q}}{min}\parallel \mathbf{p}+\mathbf{D}\mathit{\alpha }-{\mathbf{A}}_{\mathrm{d}}\mathbf{q}\parallel .$$

(22)

The combination of equations for all observation indices in ${\mathcal{A}}_m$ can be formulated as

$${\mathbf{A}}_{\mathrm{a}}\overline{\mathbf{q}}=\mathbf{p}+{\mathit{ϵ}}_{\mathrm{e}}$$

(23)

where $\mathbf{p}=\mathrm{vec}\left(\mathbf{P}\right)$, $\mathrm{vec}(·)$ denotes the vectorization operation.

$$\begin{array}{cc}\hfill {\mathit{ϵ}}_{\mathrm{e}}& ={[{\mathit{ϵ}}^{\mathrm{T}}\left(1\right),\cdots ,{\mathit{ϵ}}^{\mathrm{T}}\left(L_{\mathrm{m}}\right)]}^{\mathrm{T}}\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill {\mathbf{A}}_{\mathrm{a}}& =[{\mathbf{A}}_{\mathrm{d}},-\mathbf{D}]\in {\mathbb{R}}^{3L_{\mathrm{m}}\times (6+L_{\mathrm{m}})}\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill {\mathbf{A}}_{\mathrm{d}}& ={[{\mathbf{A}}_{\mathrm{e}}^{\mathrm{T}}\left(1\right),\cdots ,{\mathbf{A}}_{\mathrm{e}}^{\mathrm{T}}\left(L_{\mathrm{m}}\right)]}^{\mathrm{T}}\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill & =[\mathbf{1},\mathbf{t}]\otimes \mathbf{I}\in {\mathbb{R}}^{3L_{\mathrm{m}}\times 3}\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill \mathbf{D}& =diag({\mathbf{e}}_1,\cdots ,{\mathbf{e}}_{L_{\mathrm{m}}})\in {\mathbb{R}}^{3L_{\mathrm{m}}\times L_{\mathrm{m}}}\hfill \end{array}$$

(28)

where $\mathbf{t}={[t_1-t_0,t_2-t_0,\dots ,t_{L_{\mathrm{m}}}-t_0]}^{\mathrm{T}}$, and $diag(·)$ with some matrices inputs denotes a block diagonal matrix with the input matrices as block diagonal elements.

For this optimization problem, we can find the classical LS solution as

$${\widehat{\mathbf{q}}}_{\mathrm{a}}={\left({\mathbf{A}}_{\mathrm{a}}^{\mathrm{T}}{\mathbf{A}}_{\mathrm{a}}\right)}^{-1}{\mathbf{A}}_{\mathrm{a}}^{\mathrm{T}}\mathbf{p}.$$

(29)

It can be proved that for symmetric matrices $\mathbf{A}$ and $\mathbf{B}$, and $\mathbf{C}$, all of appropriate sizes,

$${\left[\begin{array}{cc}\mathbf{A}& {\mathbf{C}}^{\mathrm{T}}\\ \mathbf{C}& \mathbf{B}\end{array}\right]}^{-1}=\left[\begin{array}{cc}{(\mathbf{A}-{\mathbf{C}}^{\mathrm{T}}{\mathbf{B}}^{-1}\mathbf{C})}^{-1}& -{(\mathbf{A}-{\mathbf{C}}^{\mathrm{T}}{\mathbf{B}}^{-1}\mathbf{C})}^{-1}{\mathbf{C}}^{\mathrm{T}}\mathbf{B}\\ -{\mathbf{B}}^{-1}\mathbf{C}{(\mathbf{A}-{\mathbf{C}}^{\mathrm{T}}{\mathbf{B}}^{-1}\mathbf{C})}^{-1}& {(\mathbf{B}-\mathbf{C}{\mathbf{A}}^{-1}{\mathbf{C}}^{\mathrm{T}})}^{-1}\end{array}\right]$$

(30)

Consequently, with a fact that ${\mathbf{D}}^{\mathrm{T}}\mathbf{D}=\mathbf{I}$, we have

$${\left({\mathbf{A}}_{\mathrm{a}}^{\mathrm{T}}{\mathbf{A}}_{\mathrm{a}}\right)}^{-1}={\left[\begin{array}{cc}{\mathbf{A}}_{\mathrm{d}}^{\mathrm{T}}{\mathbf{A}}_{\mathrm{d}}& -{\mathbf{A}}_{\mathrm{d}}^{\mathrm{T}}\mathbf{D}\\ -{\mathbf{D}}^{\mathrm{T}}{\mathbf{A}}_{\mathrm{d}}& {\mathbf{D}}^{\mathrm{T}}\mathbf{D}\end{array}\right]}^{-1}=\left[\begin{array}{cc}\mathbf{X}& {\mathbf{Z}}^{\mathrm{T}}\\ \mathbf{Z}& \mathbf{Y}\end{array}\right]$$

(31)

where

$$\begin{array}{cc}\hfill \mathbf{X}& ={({\mathbf{A}}_{\mathrm{d}}^{\mathrm{T}}{\mathbf{A}}_{\mathrm{d}}-{\mathbf{A}}_{\mathrm{d}}^{\mathrm{T}}\mathbf{D}{\left({\mathbf{D}}^{\mathrm{T}}\mathbf{D}\right)}^{-1}{\mathbf{D}}^{\mathrm{T}}{\mathbf{A}}_{\mathrm{d}})}^{-1}\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill & ={({\mathbf{A}}_{\mathrm{d}}^{\mathrm{T}}{\mathbf{A}}_{\mathrm{d}}-{\mathbf{A}}_{\mathrm{d}}^{\mathrm{T}}\mathbf{D}{\mathbf{D}}^{\mathrm{T}}{\mathbf{A}}_{\mathrm{d}})}^{-1}\in {\mathbb{R}}^{6\times 6}\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill \mathbf{Y}& ={({\mathbf{D}}^{\mathrm{T}}\mathbf{D}-{\mathbf{D}}^{\mathrm{T}}{\mathbf{A}}_{\mathrm{d}}{\left({\mathbf{A}}_{\mathrm{d}}^{\mathrm{T}}{\mathbf{A}}_{\mathrm{d}}\right)}^{-1}{\mathbf{A}}_{\mathrm{d}}^{\mathrm{T}}\mathbf{D})}^{-1}\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill & ={(\mathbf{I}-{\mathbf{D}}^{\mathrm{T}}{\mathbf{A}}_{\mathrm{d}}{\left({\mathbf{A}}_{\mathrm{d}}^{\mathrm{T}}{\mathbf{A}}_{\mathrm{d}}\right)}^{-1}{\mathbf{A}}_{\mathrm{d}}^{\mathrm{T}}\mathbf{D})}^{-1}\in {\mathbb{R}}^{L_{\mathrm{m}}\times L_{\mathrm{m}}}\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill \mathbf{Z}& ={\left({\mathbf{D}}^{\mathrm{T}}\mathbf{D}\right)}^{-1}{\mathbf{D}}^{\mathrm{T}}{\mathbf{A}}_{\mathrm{d}}\mathbf{X}\hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill & ={\mathbf{D}}^{\mathrm{T}}{\mathbf{A}}_{\mathrm{d}}{({\mathbf{A}}_{\mathrm{d}}^{\mathrm{T}}{\mathbf{A}}_{\mathrm{d}}-{\mathbf{A}}_{\mathrm{d}}^{\mathrm{T}}\mathbf{D}{\mathbf{D}}^{\mathrm{T}}{\mathbf{A}}_{\mathrm{d}})}^{-1}\in {\mathbb{R}}^{L_{\mathrm{m}}\times 3}\hfill \end{array}$$

(37)

and the following equation is used in above formulations.

It can also be proved that

$$\begin{array}{cc}\hfill {\mathbf{A}}_{\mathrm{d}}^{\mathrm{T}}{\mathbf{A}}_{\mathrm{d}}& =\left[\begin{array}{cc}L\mathbf{I}& T_{\mathrm{d}}\mathbf{I}\\ T_{\mathrm{d}}\mathbf{I}& T_{\mathrm{D}}\mathbf{I}\end{array}\right]=\mathbf{M}\otimes \mathbf{I}\hfill \end{array}$$

(38)

$$\begin{array}{cc}\hfill {\mathbf{A}}_{\mathrm{d}}^{\mathrm{T}}\mathbf{D}& =\left[\begin{array}{c}\mathbf{E}\\ {\mathbf{E}}_{\mathrm{t}}\end{array}\right],\hfill \end{array}$$

(39)

where $T_{\mathrm{d}}={\sum }_{i=1}^L(t_i-t_0),T_{\mathrm{D}}={\sum }_{i=1}^L{(t_i-t_0)}^2$, ${\mathbf{E}}_{\mathrm{t}}=\mathbf{E}diag\left(\mathbf{t}\right)$

$$\mathbf{M}=\left[\begin{array}{cc}L& T_{\mathrm{d}}\\ T_{\mathrm{d}}& T_{\mathrm{D}}\end{array}\right]\iff {\mathbf{M}}^{-1}=\frac{1}{L{T}_{\mathrm{D}}-T_{\mathrm{d}}^2}\left[\begin{array}{cc}{T}_{\mathrm{D}}& -T_{\mathrm{d}}\\ -T_{\mathrm{d}}& L\end{array}\right]$$

(40)

and ⊗ denotes the Kronecker product operation.

Still evoke (30) and we have

$${\left({\mathbf{A}}_{\mathrm{d}}^{\mathrm{T}}{\mathbf{A}}_{\mathrm{d}}\right)}^{-1}={\mathbf{M}}^{-1}\otimes \mathbf{I}$$

(41)

and thus

$$\begin{array}{cc}\hfill {\mathbf{Y}}^{-1}& =\mathbf{I}-\frac{1}{L{T}_{\mathrm{D}}-T_{\mathrm{d}}^2}[{\mathbf{E}}^{\mathrm{T}},{\mathbf{E}}_{\mathrm{t}}^{\mathrm{T}}]\left[\begin{array}{c}{T}_{\mathrm{D}}\mathbf{E}-T_{\mathrm{d}}{\mathbf{E}}_{\mathrm{t}}\\ -T_{\mathrm{d}}\mathbf{E}+L{\mathbf{E}}_{\mathrm{t}}\end{array}\right]\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill & =\mathbf{I}-\frac{1}{L{T}_{\mathrm{D}}-T_{\mathrm{d}}^2}(T_{\mathrm{D}}{\mathbf{E}}^{\mathrm{T}}\mathbf{E}-T_{\mathrm{d}}{\mathbf{E}}_{\mathrm{t}}^{\mathrm{T}}\mathbf{E}-T_{\mathrm{d}}{\mathbf{E}}_{\mathrm{t}}^{\mathrm{T}}\mathbf{E}+L{\mathbf{E}}_{\mathrm{t}}^{\mathrm{T}}{\mathbf{E}}_{\mathrm{t}}).\hfill \end{array}$$

(43)

Meanwhile,

$${\mathbf{A}}_{\mathrm{a}}^{\mathrm{T}}\mathbf{p}=\left[\begin{array}{c}{\mathbf{A}}_{\mathrm{d}}^{\mathrm{T}}\mathbf{p}\\ -{\mathbf{D}}^{\mathrm{T}}\mathbf{p}\end{array}\right]=\left[\begin{array}{c}{\mathbf{q}}_{\mathrm{t}}\\ -{\mathbf{p}}_{\mathrm{e}}\end{array}\right]$$

(44)

where

$$\begin{array}{cc}\hfill {\mathbf{q}}_{\mathrm{t}}& =\mathbf{P}{[\mathbf{1},{\mathbf{t}}^{\mathrm{T}}]}^{\mathrm{T}}\hfill \end{array}$$

(45)

$$\begin{array}{cc}\hfill {\mathbf{p}}_{\mathrm{e}}& =diag\left({\mathbf{E}}^{\mathrm{T}}\mathbf{P}\right).\hfill \end{array}$$

(46)

Therefore, the estimates for ${\mathbf{q}}_m$ and $\mathit{\alpha }$ can be written separately as

$$\begin{array}{cc}\hfill {\widehat{\mathbf{q}}}_m& =\left[\begin{array}{c}{\widehat{\mathbf{g}}}_0\\ {\widehat{\mathbf{v}}}_{\mathrm{g}}\end{array}\right]=\mathbf{X}{\mathbf{q}}_{\mathrm{t}}-\mathbf{X}\left[\begin{array}{c}\mathbf{E}\\ {\mathbf{E}}_{\mathrm{t}}\end{array}\right]{\mathbf{p}}_{\mathrm{e}}\hfill \end{array}$$

(47)

$$\begin{array}{cc}\hfill {\widehat{\mathit{\alpha }}}_m& ={\mathbf{A}}_{\mathrm{d}}^{\mathrm{T}}\mathbf{D}\mathbf{X}{\mathbf{q}}_{\mathrm{t}}-\mathbf{Y}{\mathbf{p}}_{\mathrm{e}}\hfill \end{array}$$

(48)

where $\mathbf{X}$ can be written in a concise form as

$${\mathbf{X}}^{-1}=\left[\begin{array}{cc}L\mathbf{I}-\mathbf{E}{\mathbf{E}}^{\mathrm{T}}& T_{\mathrm{d}}\mathbf{I}-\mathbf{E}{\mathbf{E}}_{\mathrm{t}}^{\mathrm{T}}\\ T_{\mathrm{d}}\mathbf{I}-{\mathbf{E}}_{\mathrm{t}}{\mathbf{E}}^{\mathrm{T}}& T_{\mathrm{D}}\mathbf{I}-{\mathbf{E}}_{\mathrm{t}}{\mathbf{E}}_{\mathrm{t}}^{\mathrm{T}}\end{array}\right].$$

(49)

Under the assumption that all the sample under consideration is from the same target, the track, parameterized by $\mathbf{q}$, is identical for all observations hereafter.

2.4. The Linear LS Algorithm

In a long period, one can obtain a sequence of $\mathit{\eta }=(\theta ,\varphi )$ observations. If we take all the observations into consideration for optimization, a huge computation cost may be required. In some cases, it is also unnecessary at all. We can extract information from local observations and then transmit estimated parameters to the fusion center for target location estimation. The target location has been approximated by a linear model. Next, we express the DOA and sensor location by a linear model as well.

For DOA measures from a sensor, we can approximate a series of angle measures by

$$\mathit{\eta }=\left\{\begin{array}{c}\theta \approx {\theta }_0+{\dot{\theta }}_0(t-t_0)\hfill \\ \phi \approx {\phi }_0+{\dot{\phi }}_0(t-t_0).\hfill \end{array}\right.$$

(50)

Now consider the nth sensor and let

$$\left\{\begin{array}{c}{\theta }_0\to {\theta }_{0,n}\hfill \\ {\dot{\theta }}_0\to {\dot{\theta }}_{0,n}\hfill \end{array}\right.\Rightarrow {\overline{\mathit{\theta }}}_n=\left[\begin{array}{c}{\theta }_{0,n}\\ {\dot{\theta }}_{0,n}\end{array}\right],\left\{\begin{array}{c}{\phi }_0\to {\phi }_{0,n}\hfill \\ {\dot{\phi }}_0\to {\dot{\phi }}_{0,n}\hfill \end{array}\right.\Rightarrow {\overline{\mathit{\phi }}}_n=\left[\begin{array}{c}{\phi }_{0,n}\\ {\dot{\phi }}_{0,n}\end{array}\right].$$

(51)

For all observations in ${\mathcal{A}}_{n,m}$, we can write a polynomial regression problem as

$$\left\{\begin{array}{c}{\mathit{\theta }}_n={\mathbf{A}}_n{\overline{\mathit{\theta }}}_n+{\mathit{ϵ}}_{\theta }\hfill \\ {\mathit{\phi }}_n={\mathbf{A}}_n{\overline{\mathit{\phi }}}_n+{\mathit{ϵ}}_{\phi }\hfill \end{array}\right.$$

(52)

where ${\mathit{\theta }}_n={\left[{\theta }_i\right]}_{i\in {\mathcal{A}}_{n,m}},{\mathit{\phi }}_n={\left[{\phi }_i\right]}_{i\in {\mathcal{A}}_{n,m}}$, ${\mathit{ϵ}}_{\theta }$ denotes bias for ${\mathit{\theta }}_n$, ${\mathit{ϵ}}_{\varphi }$ denotes bias for ${\mathit{\phi }}_n$, and

$${\mathbf{A}}_{\mathrm{a}}\left(n\right)=\left[\begin{array}{c}1,t_i-t_0\\ ⋮\end{array}\right]\in {\mathbb{R}}^{K_{n,m}\times 2},i\in {\mathcal{A}}_{n,m}.$$

(53)

It can be proved that the LS estimate of the directions for the nth sensor can be directly written as

$$\begin{array}{cc}\hfill {\widehat{\mathit{\theta }}}_n& ={\left({\mathbf{A}}_{\mathrm{a}}{\left(n\right)}^{\mathrm{T}}{\mathbf{A}}_{\mathrm{a}}\left(n\right)\right)}^{-1}{\mathbf{A}}_{\mathrm{a}}{\left(n\right)}^{\mathrm{T}}{\mathit{\theta }}_n=\frac{1}{L{T}_{\mathrm{D}}-T_{\mathrm{d}}^2}\left[\begin{array}{c}{T}_{\mathrm{D}}{\mathit{\theta }}_n^{\mathrm{T}}\mathbf{1}-T_{\mathrm{d}}{\mathit{\theta }}_n^{\mathrm{T}}\mathbf{t}\\ -T_{\mathrm{d}}{\mathit{\theta }}_n^{\mathrm{T}}\mathbf{1}+L{\mathit{\theta }}_n^{\mathrm{T}}\mathbf{t}\end{array}\right]\hfill \end{array}$$

(54)

$$\begin{array}{cc}\hfill {\widehat{\mathit{\phi }}}_n& ={\left({\mathbf{A}}_{\mathrm{a}}{\left(n\right)}^{\mathrm{T}}{\mathbf{A}}_{\mathrm{a}}\left(n\right)\right)}^{-1}{\mathbf{A}}_{\mathrm{a}}{\left(n\right)}^{\mathrm{T}}{\mathit{\phi }}_n=\frac{1}{L{T}_{\mathrm{D}}-T_{\mathrm{d}}^2}\left[\begin{array}{c}{T}_{\mathrm{D}}{\mathit{\phi }}_n^{\mathrm{T}}\mathbf{1}-T_{\mathrm{d}}{\mathit{\phi }}_n^{\mathrm{T}}\mathbf{t}\\ -T_{\mathrm{d}}{\mathit{\phi }}_n^{\mathrm{T}}\mathbf{1}+L{\mathit{\phi }}_n^{\mathrm{T}}\mathbf{t}\end{array}\right].\hfill \end{array}$$

(55)

In this case, ${\widehat{\mathit{\theta }}}_n$ and ${\widehat{\mathit{\phi }}}_n$, instead of ${\mathit{\theta }}_n$ and ${\mathit{\phi }}_n$, will be transmitted to the fusion center, such that the communication cost will be greatly reduced.

In (18), ${\widehat{\mathit{\theta }}}_n$ and ${\widehat{\mathit{\phi }}}_n$ affects the positioning accuracy through the normalized direction vector ${\mathbf{e}}_i,i\in {\mathcal{A}}_{n,m}$. With a linear approximation model, $\mathbf{e}$ can be rewritten as

$$\mathbf{e}\left(t\right)\approx {\mathbf{e}}_0+{\dot{\mathbf{E}}}_0^{\mathrm{T}}{\dot{\mathit{\eta }}}_0(t-t_0)$$

(56)

where ${\mathbf{e}}_0=\mathbf{e}\left(t_0\right)$,

$$\begin{array}{cc}\hfill \dot{\theta }& =\frac{\partial \theta }{\partial t},\dot{\phi }=\frac{\partial \phi }{\partial t}\hfill \end{array}$$

(57)

$$\begin{array}{cc}\hfill {\mathit{\eta }}_0& ={[{\theta }_0,{\phi }_0]}^{\mathrm{T}},\mathit{\eta }={[\theta ,\phi ]}^{\mathrm{T}}\hfill \end{array}$$

(58)

$$\begin{array}{cc}\hfill {\dot{\mathit{\eta }}}_0& =\dot{\mathit{\eta }}{|}_{t=t_0}={[{\dot{\theta }}_0,{\dot{\phi }}_0]}^{\mathrm{T}},\dot{\mathit{\eta }}={[\dot{\theta },\dot{\phi }]}^{\mathrm{T}}\hfill \end{array}$$

(59)

$$\begin{array}{cc}\hfill {\dot{\mathbf{E}}}_0& =\dot{\mathbf{E}}{\left(\mathit{\eta }\right)|}_{t=t_0}=\dot{\mathbf{E}}\left({\mathit{\eta }}_0\right)\hfill \end{array}$$

(60)

and $\dot{\mathbf{E}}\left(\mathit{\eta }\right)$ denotes the Jacobi matrix defined by

$${\dot{\mathbf{E}}}^{\mathrm{T}}\left(\mathit{\eta }\right)=\frac{\partial \mathbf{e}}{\partial {\mathit{\eta }}^{\mathrm{T}}}=\left[\begin{array}{cc}-sin\left(\theta \right)cos\left(\phi \right)& -cos\left(\theta \right)sin\left(\phi \right)\\ cos\left(\theta \right)cos\left(\phi \right)& -sin\left(\theta \right)sin\left(\phi \right)\\ 0& cos\left(\phi \right).\end{array}\right]$$

(61)

For the nth sensor, denote ${\mathbf{e}}_0\to {\mathbf{e}}_{0,n},{\dot{\mathbf{E}}}_0\to {\dot{\mathbf{E}}}_{0,n}$, ${\dot{\mathit{\eta }}}_0\to {\dot{\mathit{\eta }}}_{0,n}$, ${\mathit{\eta }}_0\to {\mathit{\eta }}_{0,n}$ and so on.

In practice, the position of the platform is also measured by a device, such as an inertial system or a positioning system. In either case, if the nth sensor is moving with speed ${\mathbf{v}}_{\mathrm{s}}$ at $t=t_0$, the position can be expressed by

$${\mathbf{p}}_n\left(t\right)={\mathbf{p}}_{0,n}+{\mathbf{v}}_{\mathrm{s},n}(t-t_0)+{\mathit{ϵ}}_{\mathrm{s}}\left(n\right)$$

(62)

where ${\mathbf{p}}_{0,n}$ denotes the position of the nth sensor, ${\mathbf{v}}_{\mathrm{s},n}$ denotes the velocity, both at $t=t_0$, and ${\mathit{ϵ}}_{\mathrm{s}}$ denotes the bias of position estimation error.

For simplicity, we assume that the sensor location is measured at $t={\overline{t}}_{k,n},k=1,\cdots ,{\overline{N}}_n$. In practice, in this configuration, the sensor location can be measured at instants other than $t_{k,n},k=1,\cdots ,N_n$. Now we can construct equations as

$${\mathbf{p}}_{k,n}={\mathbf{p}}_{0,n}+{\mathbf{v}}_{\mathrm{s},n}({\overline{t}}_{k,n}-t_0)+{\mathit{ϵ}}_{\mathrm{s}}(k,n),k=1,\cdots ,{\overline{N}}_n,n=1,\cdots ,N$$

(63)

or in another form as

$${\mathbf{p}}_{\mathrm{s},n}={\mathbf{A}}_{\mathrm{s}}{\overline{\mathbf{p}}}_n+{\mathit{ϵ}}_{\mathrm{s}}\left(n\right)$$

(64)

where ${\overline{\mathbf{p}}}_n={[{\mathbf{p}}_{0,n}^{\mathrm{T}},{\mathbf{v}}_{\mathrm{s},n}^{\mathrm{T}}]}^{\mathrm{T}}$, ${\mathbf{p}}_{\mathrm{s},n}={[{\mathbf{p}}_{1,n}^{\mathrm{T}},\cdots ,{\mathbf{p}}_{{\overline{N}}_n,n}^{\mathrm{T}}]}^{\mathrm{T}}$, and

$${\mathbf{A}}_{\mathrm{s}}\left(n\right)={\left[\begin{array}{c}\mathbf{I},({\overline{t}}_{k,n}-t_0)\mathbf{I}\\ ⋮\end{array}\right]}_{k=1,\cdots ,{\overline{N}}_n}=[\mathbf{1},{\mathbf{t}}_{\mathrm{s}}]\otimes \mathbf{I}$$

(65)

for which

$${\left({\mathbf{A}}_{\mathrm{s}}^{\mathrm{T}}\left(n\right){\mathbf{A}}_{\mathrm{s}}\left(n\right)\right)}^{-1}={\mathbf{M}}_{\mathrm{s}}^{-1}\otimes \mathbf{I}$$

(66)

where ${\mathbf{M}}_{\mathrm{s}}^{-1}$ is similar to $\mathbf{M}$ in (40).

The LS estimate of $\overline{\mathbf{p}}$ is

$$\begin{array}{cc}\hfill {\widehat{\mathbf{p}}}_n& =\left[\begin{array}{c}{\widehat{\mathbf{p}}}_{0,n}^{\mathrm{T}}\\ {\widehat{\mathbf{v}}}_{\mathrm{s},n}\end{array}\right]={\left({\mathbf{A}}_{\mathrm{s}}^{\mathrm{T}}\left(n\right){\mathbf{A}}_{\mathrm{s}}\left(n\right)\right)}^{-1}{\mathbf{A}}_{\mathrm{s}}^{\mathrm{T}}\left(n\right){\mathbf{p}}_{\mathrm{s},n}\hfill \end{array}$$

(67)

$$\begin{array}{cc}\hfill & =({\mathbf{M}}_{\mathrm{s}}^{-1}\otimes \mathbf{I})\left[\begin{array}{c}\mathbf{P}\mathbf{1}\\ {\mathbf{P}}_{\mathrm{t}}\mathbf{t}\end{array}\right]\hfill \end{array}$$

(68)

$$\begin{array}{cc}\hfill {\widehat{\mathbf{p}}}_{0,n}& =\frac{1}{L{T}_{\mathrm{D}}-T_{\mathrm{d}}^2}(T_{\mathrm{D}}\mathbf{P}\mathbf{1}-T_{\mathrm{d}}{\mathbf{P}}_{\mathrm{t}}\mathbf{t})\hfill \end{array}$$

(69)

and

$$\begin{array}{cc}\hfill {\widehat{\mathbf{v}}}_{\mathrm{s},n}& =\frac{1}{L{T}_{\mathrm{D}}-T_{\mathrm{d}}^2}(-T_{\mathrm{d}}\mathbf{P}\mathbf{1}+L{\mathbf{P}}_{\mathrm{t}}\mathbf{t}).\hfill \end{array}$$

(70)

With above operations, we can obtain a linear sensor location parameter ${\widehat{\mathbf{p}}}_n$, linear DOA parameters ${\widehat{\mathit{\theta }}}_n,{\widehat{\mathit{\varphi }}}_n$, and linear target location parameters ${\mathbf{q}}_m$. The accuracy depends on the interval of observations, the speed of the target, and the speed of the sensors. A series of observations can now be approximated by two parameters and it is now unnecessary to transmit all the local observation to a fusion center anymore.

With only one observation available, a fusion center can now estimate the target position and velocity with the following equation

$$\begin{array}{cc}\hfill \mathbf{x}\left(t\right)& ={\mathbf{g}}_0+{\mathbf{v}}_{\mathrm{g}}(t-t_0)+{\mathit{ϵ}}_{\mathrm{a}}\left(n\right)={\mathbf{p}}_{0,n}+{\mathbf{v}}_{\mathrm{s},n}(t-t_0)+{\alpha }_n\left(t\right){\mathbf{e}}_n\left(t\right)\hfill \\ \hfill & \approx {\mathbf{p}}_{0,n}+{\mathbf{v}}_{\mathrm{s},n}(t-t_0)+({\alpha }_{0,n}+{\dot{\alpha }}_{0,n}(t-t_0))({\mathbf{e}}_{0,n}+{\dot{\mathbf{E}}}_{0,n}^{\mathrm{T}}{\dot{\mathit{\eta }}}_{0,n}(t-t_0))\hfill \end{array}$$

(71)

where ${\dot{\alpha }}_{0,n}$ represents the change rate of ${\alpha }_{0,n}$. By expanding Equation (71) and ignoring the second-order term, ${\mathit{ϵ}}_{\mathrm{a}}\left(n\right)$ can be reformulated as

$$\begin{array}{cc}\hfill {\mathit{ϵ}}_{\mathrm{a}}\left(n\right)& ={\mathbf{p}}_{0,n}+{\alpha }_{0,n}{\mathbf{e}}_{0,n}-{\mathbf{g}}_0+({\alpha }_{0,n}{\dot{\mathbf{E}}}_{0,n}^{\mathrm{T}}{\dot{\mathit{\eta }}}_{0,n}+{\mathbf{v}}_{\mathbf{s},n}+{\dot{\alpha }}_{0,n}{\mathbf{e}}_{0,n}-{\mathbf{v}}_g)(t-t_0)\hfill \end{array}$$

(72)

where $n=1,2\dots ,N$ and ${\mathit{ϵ}}_{\mathrm{a}}\left(n\right)$ denotes the bias term in approximation target position and direction of arrival by the linear models.

The following equation can be obtained using (72) for N sensors,

$$\begin{array}{c}\hfill {\mathit{ϵ}}_{\mathrm{a}}={\mathbf{p}}_0+{\mathbf{D}}_0{\mathit{\alpha }}_0-\overline{\mathbf{I}}{\mathbf{g}}_0+({\mathbf{X}}_{\mathrm{d}}{\mathit{\alpha }}_0+{\mathbf{v}}_s+{\mathbf{D}}_0{\dot{\mathit{\alpha }}}_0-\overline{\mathbf{I}}{\mathbf{v}}_g)(t-t_0)\end{array}$$

(73)

where

$$\begin{array}{cc}\hfill {\mathit{\alpha }}_0& ={[{\alpha }_{0,1},{\alpha }_{0,2},\cdots ,{\alpha }_{0,N}]}^{\mathrm{T}}\hfill \end{array}$$

(74)

$$\begin{array}{cc}\hfill {\dot{\mathit{\alpha }}}_0& ={[{\dot{\alpha }}_{0,1},{\dot{\alpha }}_{0,2},\cdots ,{\dot{\alpha }}_{0,N}]}^{\mathrm{T}}\hfill \end{array}$$

(75)

$$\begin{array}{cc}\hfill {\mathbf{p}}_0& ={[{\mathbf{p}}_{0,1}^{\mathrm{T}},{\mathbf{p}}_{0,2}^{\mathrm{T}},\dots ,{\mathbf{p}}_{0,N}^{\mathrm{T}}]}^{\mathrm{T}}\hfill \end{array}$$

(76)

$$\begin{array}{cc}\hfill {\mathbf{D}}_0& =\mathrm{diag}({\mathbf{e}}_{0,1},{\mathbf{e}}_{0,2},\dots ,{\mathbf{e}}_{0,N})\hfill \end{array}$$

(77)

$$\begin{array}{cc}\hfill {\mathbf{X}}_{\mathrm{d}}& =\mathrm{diag}({\dot{\mathbf{E}}}_{0,1}^{\mathrm{T}}{\dot{\mathit{\eta }}}_{0,1},{\dot{\mathbf{E}}}_{0,2}^{\mathrm{T}}{\dot{\mathit{\eta }}}_{0,2},\dots ,{\dot{\mathbf{E}}}_{0,N}^{\mathrm{T}}{\dot{\mathit{\eta }}}_{0,N})\hfill \end{array}$$

(78)

$$\begin{array}{cc}\hfill {\mathbf{v}}_{\mathrm{s}}& ={[{\mathbf{v}}_{\mathrm{s},1}^{\mathrm{T}},{\mathbf{v}}_{\mathrm{s},2}^{\mathrm{T}},\dots ,{\mathbf{v}}_{\mathrm{s},N}^{\mathrm{T}}]}^{\mathrm{T}}\hfill \end{array}$$

(79)

$$\begin{array}{cc}\hfill \overline{\mathbf{I}}& ={\mathbf{1}}_N\otimes {\mathbf{I}}_3\in {\mathbb{R}}^{3N\times 3}\hfill \end{array}$$

(80)

and

$$\begin{array}{cc}\hfill {\mathit{ϵ}}_{\mathrm{a}}& ={[{\mathit{ϵ}}_{\mathrm{a}}^{\mathrm{T}}\left(1\right),{\mathit{ϵ}}_{\mathrm{a}}^{\mathrm{T}}\left(2\right),\cdots ,{\mathit{ϵ}}_{\mathrm{a}}^{\mathrm{T}}\left(N\right)]}^{\mathrm{T}}.\hfill \end{array}$$

(81)

In order to minimize the total bias ${\mathit{ϵ}}_a$, we can minimize

$$\begin{array}{c}\hfill \underset{{\mathbf{q}}_{\mathrm{c}},{\mathbf{q}}_{\mathrm{v}}}{min}\parallel {\mathbf{p}}_0+{\mathbf{D}}_0{\mathit{\alpha }}_0-\overline{\mathbf{I}}{\mathbf{g}}_0+({\mathbf{X}}_{\mathrm{d}}{\mathit{\alpha }}_0+{\mathbf{v}}_s+{\mathbf{D}}_0{\dot{\mathit{\alpha }}}_0-\overline{\mathbf{I}}{\mathbf{v}}_g)(t-t_0)\parallel ,t\in \mathcal{T}.\end{array}$$

(82)

where ${\mathbf{q}}_{\mathrm{c}}={[{\mathbf{g}}_0^{\mathrm{T}},{\mathit{\alpha }}_0^{\mathrm{T}}]}^{\mathrm{T}}$ and ${\mathbf{q}}_{\mathrm{v}}={[{\mathbf{v}}_{\mathrm{g}}^{\mathrm{T}},{\dot{\mathit{\alpha }}}_0^{\mathrm{T}}]}^{\mathrm{T}}$.

To ensure the bias is minimized for $t\in \mathcal{T}$, both the initial position bias, ${\mathbf{p}}_0+{\mathbf{D}}_0{\mathit{\alpha }}_0-\overline{\mathbf{I}}{\mathbf{g}}_0$, and the speed bias, ${\mathbf{X}}_{\mathrm{d}}{\mathit{\alpha }}_0+{\mathbf{v}}_s+{\mathbf{D}}_0{\dot{\mathit{\alpha }}}_0-\overline{\mathbf{I}}{\mathbf{v}}_g$ should be minimized. Therefore, we can solve the optimization problem through solving the following two optimization problems of smaller scale,

$$\begin{array}{cc}\hfill & \underset{{\mathbf{q}}_{\mathrm{c}}}{min}\parallel {\mathbf{p}}_0+{\mathbf{D}}_0{\mathit{\alpha }}_0-\overline{\mathbf{I}}{\mathbf{g}}_0\parallel \hfill \end{array}$$

(83)

$$\begin{array}{cc}\hfill & \underset{{\mathbf{q}}_{\mathrm{v}}}{min}\parallel {\mathbf{X}}_{\mathrm{d}}{\mathit{\alpha }}_0+{\mathbf{v}}_{\mathrm{s}}+{\mathbf{D}}_0{\dot{\mathit{\alpha }}}_0-\overline{\mathbf{I}}{\mathbf{v}}_{\mathrm{g}}\parallel .\hfill \end{array}$$

(84)

The solutions to the problems can be found directly through the LS algorithm as

$$\begin{array}{cc}\hfill {\widehat{\mathbf{q}}}_{\mathrm{c}}& ={[{\widehat{\mathbf{g}}}_0^{\mathrm{T}},{\widehat{\mathit{\alpha }}}_0^{\mathrm{T}}]}^{\mathrm{T}}={\left({\mathbf{A}}_{\mathrm{c}}^{\mathrm{T}}{\mathbf{A}}_{\mathrm{c}}\right)}^{-1}{\mathbf{A}}_{\mathrm{c}}^{\mathrm{T}}{\mathbf{p}}_0\hfill \end{array}$$

(85)

$$\begin{array}{cc}\hfill {\widehat{\mathbf{q}}}_{\mathrm{v}}& ={[{\widehat{\mathbf{v}}}_{\mathrm{g}}^{\mathrm{T}},{\widehat{\dot{\alpha }}}_0^{\mathrm{T}}]}^{\mathrm{T}}={\left({\mathbf{A}}_{\mathrm{c}}^{\mathrm{T}}{\mathbf{A}}_{\mathrm{c}}\right)}^{-1}{\mathbf{A}}_{\mathrm{c}}^{\mathrm{T}}({\mathbf{X}}_{\mathrm{d}}{\mathit{\alpha }}_0+{\mathbf{v}}_{\mathrm{s}})\hfill \end{array}$$

(86)

where

$$\begin{array}{cc}\hfill {\mathbf{A}}_{\mathrm{c}}& =[{\mathbf{1}}_N\otimes {\mathbf{I}}_3,-{\mathbf{D}}_0].\hfill \end{array}$$

(87)

It can be proved that

$$\begin{array}{c}\hfill {\mathbf{A}}_{\mathrm{c}}^{\mathrm{T}}{\mathbf{A}}_{\mathrm{c}}=\left[\begin{array}{cc}{\overline{\mathbf{I}}}^{\mathrm{T}}\overline{\mathbf{I}}& -{\overline{\mathbf{I}}}^{\mathbf{T}}{\mathbf{D}}_0\\ -{\mathbf{D}}_0^{\mathrm{T}}\overline{\mathbf{I}}& {\mathbf{D}}_0^{\mathrm{T}}{\mathbf{D}}_0\end{array}\right]=\left[\begin{array}{cc}N{\mathbf{I}}_3& -{\mathbf{E}}_0\\ -{\mathbf{E}}_0^{\mathrm{T}}& \mathbf{I}\end{array}\right]\end{array}$$

(88)

and then

$${\left({\mathbf{A}}_{\mathrm{c}}^{\mathrm{T}}{\mathbf{A}}_{\mathrm{c}}\right)}^{-1}=\left[\begin{array}{cc}{(N\mathbf{I}-{\mathbf{E}}_0{\mathbf{E}}_0^{\mathrm{T}})}^{-1}& {(N\mathbf{I}-{\mathbf{E}}_0{\mathbf{E}}_0^{\mathrm{T}})}^{-1}{\mathbf{E}}_0\\ {\mathbf{E}}_0^{\mathrm{T}}{(N\mathbf{I}-{\mathbf{E}}_0{\mathbf{E}}_0^{\mathrm{T}})}^{-1}& \mathbf{I}+{\mathbf{E}}_0^{\mathrm{T}}{(N\mathbf{I}-{\mathbf{E}}_0{\mathbf{E}}_0^{\mathrm{T}})}^{-1}{\mathbf{E}}_0\end{array}\right].$$

(89)

where ${\mathbf{E}}_0=[{\mathbf{e}}_{0,1},{\mathbf{e}}_{0,2},\dots ,{\mathbf{e}}_{0,N}]$ and ${\mathbf{P}}_0=[{\mathbf{p}}_{0,1},{\mathbf{p}}_{0,2},\dots ,{\mathbf{p}}_{0,N}]$.

Meanwhile,

$$\begin{array}{c}\hfill {\mathbf{A}}_{\mathrm{c}}^{\mathrm{T}}{\mathbf{p}}_0=\left[\begin{array}{c}{\mathbf{P}}_0\mathbf{1}\\ -{\mathbf{D}}_0^{\mathrm{T}}{\mathbf{p}}_0\end{array}\right]\end{array}$$

(90)

and thus,

$$\begin{array}{cc}\hfill {\widehat{\mathbf{q}}}_0& ={[{\widehat{\mathbf{g}}}_0^{\mathrm{T}},{\widehat{\mathit{\alpha }}}_0^{\mathrm{T}}]}^{\mathrm{T}}={\left({\mathbf{A}}_{\mathrm{c}}^{\mathrm{T}}{\mathbf{A}}_{\mathrm{c}}\right)}^{-1}{\mathbf{A}}_{\mathrm{c}}^{\mathrm{T}}{\mathbf{p}}_0\hfill \\ \hfill {\widehat{\mathbf{g}}}_0& ={(N\mathbf{I}-{\mathbf{E}}_0{\mathbf{E}}_0^{\mathrm{T}})}^{-1}{\mathbf{P}}_0\mathbf{1}-{(N\mathbf{I}-{\mathbf{E}}_0{\mathbf{E}}_0^{\mathrm{T}})}^{-1}{\mathbf{E}}_0{\mathbf{D}}_0^{\mathrm{T}}{\mathbf{p}}_0\hfill \end{array}$$

(91)

$$\begin{array}{cc}\hfill & ={(N\mathbf{I}-{\mathbf{E}}_0{\mathbf{E}}_0^{\mathrm{T}})}^{-1}({\mathbf{P}}_0\mathbf{1}-{\mathbf{E}}_0{\mathbf{D}}_0^{\mathrm{T}}{\mathbf{p}}_0)\hfill \\ \hfill {\widehat{\mathit{\alpha }}}_0& ={\mathbf{D}}_0^{\mathrm{T}}\overline{\mathbf{I}}{\widehat{\mathbf{g}}}_0-{\mathbf{D}}_0^{\mathrm{T}}{\mathbf{p}}_0\hfill \end{array}$$

(92)

$$\begin{array}{cc}\hfill & ={\mathbf{D}}_0^{\mathrm{T}}\overline{\mathbf{I}}{(N\mathbf{I}-{\mathbf{E}}_0{\mathbf{E}}_0^{\mathrm{T}})}^{-1}({\mathbf{P}}_0\mathbf{1}-{\mathbf{E}}_0{\mathbf{D}}_0^{\mathrm{T}}{\mathbf{p}}_0)-{\mathbf{D}}_0^{\mathrm{T}}{\mathbf{p}}_0\hfill \end{array}$$

(93)

which is identical to the solution of (82). A minor difference is that the matrix inverse operation is over an $N\times N$ matrix $N\mathbf{I}-{\mathbf{E}}_0{\mathbf{E}}_0^{\mathrm{T}}$.

The solution to ${\mathbf{q}}_{\mathrm{v}}$ can be expressed by

$$\begin{array}{cc}\hfill {\mathbf{A}}_{\mathrm{c}}^{\mathrm{T}}({\mathbf{X}}_{\mathrm{d}}{\mathit{\alpha }}_0+{\mathbf{v}}_{\mathrm{s}})& =\left[\begin{array}{c}{\mathbf{1}}_N^{\mathrm{T}}\otimes {\mathbf{I}}_3\\ -{\mathbf{D}}_0^{\mathrm{T}}\end{array}\right](\mathrm{vec}(\mathbf{X}diag\left({\mathit{\alpha }}_0\right))+\mathrm{vec}\left({\mathbf{V}}_{\mathrm{s}}\right))\hfill \\ \hfill & =\left[\begin{array}{c}\mathbf{X}{\mathit{\alpha }}_0+{\mathbf{V}}_{\mathrm{s}}\mathbf{1}\\ -{\mathbf{v}}_{\mathrm{a}}-{\mathbf{v}}_{\mathrm{e}}\end{array}\right]\hfill \end{array}$$

(94)

where

$$\begin{array}{cc}\hfill \mathbf{X}& =[{\dot{\mathbf{E}}}_{0,1}^{\mathrm{T}}{\dot{\mathit{\eta }}}_{0,1}{\dot{\mathbf{E}}}_{0,2}^{\mathrm{T}}{\dot{\mathit{\eta }}}_{0,2},\dots ,{\dot{\mathbf{E}}}_{0,N}^{\mathrm{T}}{\dot{\mathit{\eta }}}_{0,N}]\hfill \end{array}$$

(95)

$$\begin{array}{cc}\hfill {\mathbf{V}}_{\mathrm{s}}& =[{\mathbf{v}}_{\mathrm{s},1},{\mathbf{v}}_{\mathrm{s},2},\dots ,{\mathbf{v}}_{\mathrm{s},N}]\hfill \\ \hfill {\mathbf{v}}_{\mathrm{a}}& =diag\left({\mathit{\alpha }}_0\right){\mathbf{E}}_0^{\mathrm{T}}\mathbf{X}\hfill \end{array}$$

(96)

$$\begin{array}{cc}\hfill & ={[{\alpha }_{0,1}{\mathbf{e}}_{0,1}^{\mathrm{T}}{\dot{\mathbf{E}}}_{0,1}^{\mathrm{T}}{\dot{\mathit{\eta }}}_{0,1},\cdots ,{\alpha }_{0,N}{\mathbf{e}}_{0,N}^{\mathrm{T}}{\dot{\mathbf{E}}}_{0,N}^{\mathrm{T}}{\dot{\mathit{\eta }}}_{0,N}]}^{\mathrm{T}}\hfill \end{array}$$

(97)

and

$$\begin{array}{cc}\hfill {\mathbf{v}}_{\mathrm{e}}& ={[{\mathbf{e}}_{0,1}^{\mathrm{T}}{\mathbf{v}}_{\mathrm{s},1},\cdots ,{\mathbf{e}}_{0,N}^{\mathrm{T}}{\mathbf{v}}_{\mathrm{s},N}]}^{\mathrm{T}}.\hfill \end{array}$$

(98)

Consequently, we can obtain

$$\begin{array}{cc}\hfill {\widehat{\mathbf{v}}}_{\mathrm{g}}& ={(N\mathbf{I}-{\mathbf{E}}_0{\mathbf{E}}_0^{\mathrm{T}})}^{-1}(\mathbf{X}{\mathit{\alpha }}_0+{\mathbf{V}}_{\mathrm{s}}\mathbf{1}-{\mathbf{E}}_0{\mathbf{v}}_{\mathrm{a}}-{\mathbf{E}}_0{\mathbf{v}}_{\mathrm{e}})\hfill \end{array}$$

(99)

$$\begin{array}{cc}\hfill {\widehat{\dot{\mathit{\alpha }}}}_0& ={\mathbf{E}}_0^{\mathrm{T}}{(N\mathbf{I}-{\mathbf{E}}_0{\mathbf{E}}_0^{\mathrm{T}})}^{-1}(\mathbf{X}{\mathit{\alpha }}_0+{\mathbf{V}}_{\mathrm{s}}\mathbf{1})-{\mathbf{v}}_{\mathrm{a}}-{\mathbf{v}}_{\mathrm{e}}-{\mathbf{E}}_0^{\mathrm{T}}{(N\mathbf{I}-{\mathbf{E}}_0{\mathbf{E}}_0^{\mathrm{T}})}^{-1}{\mathbf{E}}_0({\mathbf{v}}_{\mathrm{a}}+{\mathbf{v}}_{\mathrm{e}})\hfill \end{array}$$

(100)

$$\begin{array}{cc}\hfill & ={\mathbf{E}}_0^{\mathrm{T}}{(N\mathbf{I}-{\mathbf{E}}_0{\mathbf{E}}_0^{\mathrm{T}})}^{-1}(\mathbf{X}{\mathit{\alpha }}_0+{\mathbf{V}}_{\mathrm{s}}\mathbf{1}-{\mathbf{E}}_0{\mathbf{v}}_{\mathrm{a}}-{\mathbf{E}}_0{\mathbf{v}}_{\mathrm{e}})-{\mathbf{v}}_{\mathrm{a}}-{\mathbf{v}}_{\mathrm{e}}\hfill \end{array}$$

(101)

$$\begin{array}{cc}\hfill & ={\mathbf{E}}_0^{\mathrm{T}}{\widehat{\mathbf{v}}}_{\mathrm{g}}-{\mathbf{v}}_{\mathrm{a}}-{\mathbf{v}}_{\mathrm{e}}.\hfill \end{array}$$

(102)

One can also derive in another way. According to (93), the change rate ${\dot{\mathit{\alpha }}}_0$ of ${\mathit{\alpha }}_0$ can be expressed as

$$\begin{array}{c}\hfill {\dot{\mathit{\alpha }}}_0={\mathbf{X}}_{\mathrm{d}}^{\mathrm{T}}\overline{\mathbf{I}}{\widehat{\mathbf{g}}}_0+{\mathbf{D}}_0^{\mathrm{T}}\overline{\mathbf{I}}{\mathbf{v}}_{\mathrm{g}}-{\mathbf{X}}_{\mathrm{d}}^{\mathrm{T}}{\mathbf{p}}_0-{\mathbf{D}}_0^{\mathrm{T}}{\mathbf{v}}_{\mathrm{s}}.\end{array}$$

(103)

Take (103) into (84), which can be rewritten as

$$\begin{array}{c}\hfill \underset{{\mathbf{v}}_{\mathrm{g}}}{min}\parallel {\mathbf{X}}_{\mathrm{d}}{\widehat{\mathit{\alpha }}}_0+{\mathbf{v}}_{\mathrm{s}}+{\mathbf{D}}_0{\mathbf{X}}_{\mathrm{d}}^{\mathrm{T}}\overline{\mathbf{I}}{\widehat{\mathbf{g}}}_0-{\mathbf{D}}_0{\mathbf{D}}_0^{\mathrm{T}}{\mathbf{v}}_s-{\mathbf{D}}_0{\mathbf{X}}_{\mathrm{d}}^{\mathrm{T}}{\mathbf{p}}_0+({\mathbf{D}}_0{\mathbf{D}}_0^{\mathrm{T}}-\mathbf{I})\overline{\mathbf{I}}{\mathbf{v}}_{\mathrm{g}}\parallel .\end{array}$$

(104)

The solution of (104) can also be found through the LS algorithm, which can be expressed as

$$\begin{array}{c}\hfill {\widehat{\mathbf{v}}}_{\mathrm{g}}={\left({\mathbf{A}}_{\mathrm{v}}^{\mathrm{T}}{\mathbf{A}}_{\mathrm{v}}\right)}^{-1}{\mathbf{A}}_{\mathrm{v}}^{\mathrm{T}}{\mathbf{b}}_{\mathrm{v}}\end{array}$$

(105)

where

$$\begin{array}{cc}\hfill {\mathbf{A}}_{\mathrm{v}}& =(\mathbf{I}-{\mathbf{D}}_0{\mathbf{D}}_0^{\mathrm{T}})\overline{\mathbf{I}}\hfill \end{array}$$

(106)

$$\begin{array}{cc}\hfill {\mathbf{b}}_{\mathrm{v}}& ={\mathbf{X}}_{\mathrm{d}}{\widehat{\mathit{\alpha }}}_0+{\mathbf{v}}_{\mathrm{s}}+{\mathbf{D}}_0{\mathbf{X}}_{\mathrm{d}}^{\mathrm{T}}\overline{\mathbf{I}}{\widehat{\mathbf{g}}}_0-{\mathbf{D}}_0{\mathbf{D}}_0^{\mathrm{T}}{\mathbf{v}}_s-{\mathbf{D}}_0{\mathbf{X}}_{\mathrm{d}}^{\mathrm{T}}{\mathbf{p}}_0.\hfill \end{array}$$

(107)

2.5. The Truncated LS Algorithm

For a better performance, it is necessary to estimate the change rate $\dot{\mathit{\alpha }}$ and if we ignore this term, the optimization problem becomes

$$\underset{{\mathbf{v}}_{\mathrm{g}}}{min}\parallel {\mathbf{X}}_{\mathrm{d}}{\mathit{\alpha }}_0+{\mathbf{v}}_{\mathrm{s}}-\overline{\mathbf{I}}{\mathbf{v}}_{\mathrm{g}}\parallel$$

(108)

whose solution, termed as truncated LS algorithm subsequently, is

$${\widehat{\mathbf{v}}}_{\mathrm{g}}=\frac{1}{N}({\mathbf{V}}_{\mathrm{s}}\mathbf{1}+\mathbf{X}{\widehat{\mathit{\alpha }}}_0)$$

(109)

which is an average operation. Note that the truncated LS algorithm shares the same position estimate with the linear LS algorithm.

The estimate of the target location at $t\in \mathcal{T}$ can be written as

$${\widehat{\mathbf{g}}}_m\left(t\right)={\widehat{\mathbf{g}}}_0+{\widehat{\mathbf{v}}}_{\mathrm{g}}(t-t_0),t\in \mathcal{T}.$$

(110)

In (99) and (109), the velocity terms ${\mathbf{v}}_{\mathrm{s},k},k=1,\cdots ,N$ are unknown and should be replaced by their estimates, typically ${\widehat{\mathbf{v}}}_{\mathrm{s}}$ estimated in the LS algorithm as in (70). In practice, besides the linear regression method performed at local sensors, there may be other methods that can output more accurate velocity and angle difference information. For instance, some inertial devices can measure the velocity more accurately than the LS algorithm in use. With a more accurate velocity estimate, it is possible to obtain a better positioning performance.

Both the linear LS algorithm and the truncated LS algorithm estimate the velocity of the target and thus can make the time-consuming nonlinear filtering operation update in a longer time interval. Subsequently, the performances of these algorithms will be analyzed in numerical results.

3. Numerical Results

In order to evaluate the performance of the concerned positioning algorithms, we first consider a scenario where four angle-only sensors are estimating the position of a target with their angle-only observations. Both the sensors and the target are moving with a constant speed during the period of observation by assumption. The initial position and the constant speed of the sensors and the target are shown in

Table 1. Positions and velocities of sensors and the target.

	Position (m) at t = 0 s	Velocity (m/s)
Sensor #1	${\mathbf{p}}_1^{\mathrm{o}}\left(0\right)={[1000,1000,0]}^{\mathrm{T}}$	${[-100,0,0]}^{\mathrm{T}}$
Sensor #2	${\mathbf{p}}_2^{\mathrm{o}}\left(0\right)={[1000,2000,0]}^{\mathrm{T}}$	${[-100,-80,0]}^{\mathrm{T}}$
Sensor #3	${\mathbf{p}}_3^{\mathrm{o}}\left(0\right)={[2000,1000,0]}^{\mathrm{T}}$	${[-100,-50,0]}^{\mathrm{T}}$
Sensor #4	${\mathbf{p}}_4^{\mathrm{o}}\left(0\right)={[1500,1500,0]}^{\mathrm{T}}$	${[-100,-60,0]}^{\mathrm{T}}$
Target #1	${\mathbf{g}}_1^{\mathrm{o}}\left(0\right)={[0,100,1000]}^{\mathrm{T}}$	${[200,100,0]}^{\mathrm{T}}$

. The scenario is illustrated in Figure 2.

All the sensors output observations at a frequency of 50 Hz, i.e., with a period of 20 ms. But they operate on an asynchronous manner, namely the sensors record the observations at independent instants. The differences of the sampling instants are randomly generated within 20ms. This assumption is important in real situations because it allows distributed sensors to operate asynchronously. We also assume that there is no error in recording the instants of the observations and for all the sensors, no signal is missed in detection during the observation period.

Assume that the self-positioning error is distributed with zero-mean normal distribution, whose variance is 1m for all the sensors, namely

$${\mathbf{C}}_{\mathrm{s}}(k,n)=\mathbb{E}\left(\mathsf{\Delta }{\mathbf{p}}_n\left(t_{k,n}\right)\mathsf{\Delta }{\mathbf{p}}_n^{\mathrm{T}}\left(t_{k,n}\right)\right)=\mathbf{I},n=1,\cdots ,N.$$

(111)

The angle measurement error also follows zero-mean normal distribution with variance of $0.5$ degree for all the observations, namely

$${\mathbf{C}}_{\mathit{\eta }}\left(i\right)=\mathbb{E}\left(\mathsf{\Delta }{\mathit{\eta }}_i\mathsf{\Delta }{\mathit{\eta }}_i^{\mathrm{T}}\right)=0.5\mathbf{I},i\in {\mathcal{A}}_m.$$

(112)

At the current stage, we do not consider the measurement errors from the gyroscopes installed on the platforms along with the sensors. Therefore, the angle measurement error is caused by the sensors only.

In order to evaluate the performance of the algorithms, we run $N_{\mathrm{e}}=20$ random experiments and take the root mean square error (RMSE) as the resulting performance metric. The RMSE of position, RMSE of velocity and the gross RMSE at instant t in scale and in dB are defined by

$$\begin{array}{ccc}\hfill {\mathrm{RMSE}}_{\mathrm{p}}\left(t\right)& =\sqrt{\frac{1}{N_{\mathrm{e}}}\sum _{k=1}^{N_{\mathrm{e}}}{|{\widehat{\mathbf{g}}}_0(t;k)-{\mathbf{g}}_1^{\mathrm{o}}\left(0\right)|}^2},\hfill & \hfill {\mathrm{RMSE}}_{\mathrm{p}}\left(t\right):\mathrm{dB}=20{log}_{10}{\mathrm{RMSE}}_{\mathrm{p}}\left(t\right)\end{array}$$

(113)

$$\begin{array}{ccc}\hfill {\mathrm{RMSE}}_{\mathrm{v}}\left(t\right)& =\sqrt{\frac{1}{N_{\mathrm{e}}}\sum _{k=1}^{N_{\mathrm{e}}}{|{\widehat{\mathbf{v}}}_{\mathrm{g}}(t;k)-{\mathbf{v}}_{\mathrm{g},1}^{\mathrm{o}}\left(0\right)|}^2},\hfill & \hfill {\mathrm{RMSE}}_{\mathrm{v}}\left(t\right):\mathrm{dB}=20{log}_{10}{\mathrm{RMSE}}_{\mathrm{v}}\left(t\right)\end{array}$$

(114)

$$\begin{array}{ccc}\hfill {\mathrm{RMSE}}_{\mathrm{g}}\left(t\right)& =\sqrt{{\mathrm{RMSE}}_{\mathrm{p}}^2\left(t\right)+{\mathrm{RMSE}}_{\mathrm{v}}^2\left(t\right)},\hfill & \hfill {\mathrm{RMSE}}_{\mathrm{g}}\left(t\right):\mathrm{dB}=20{log}_{10}{\mathrm{RMSE}}_{\mathrm{g}}\left(t\right)\end{array}$$

(115)

respectively, where ${\widehat{\mathbf{g}}}_0(t;k)$ denotes the initial position of the target at the kth experiment, ${\mathbf{g}}_1^{\mathrm{o}}\left(0\right)$ and ${\mathbf{v}}_1^{\mathrm{o}}\left(0\right)$ are constants during experiments, and ${\widehat{\mathbf{v}}}_{\mathrm{g}}(t;k)$ denotes the estimate of the target speed at the kth experiment. At each random experiment, the position error and the angle measurement error are generated randomly.

3.1. The Convergence Curves

As the number of observations increase, the localization performance will improve. Figure 3 shows the mean RMSE of position and velocity and the gross RMSE of the gross LS algorithm, the linear LS algorithm and the truncated LS algorithm.

From Figure 3a,b, it can be seen that with observations in a short while, roughly in about $0.8$ s corresponding to 40 observations, all the algorithms have close gross RMSE curves. However, as more observations are available, the linear LS algorithm and the truncated LS algorithm will perform worse and the position RMSE even increase with sample number. It is a predictable result, because the linear approximations of the target and platform motion will be inaccurate gradually, resulting a deteriorated positioning performance. The gross LS algorithm will always benefit from the increase of the observations, because it does not rely on the linear approximation, and more observations will contribute more information of the target position.

From Figure 3c,d, with some initial observations, the truncated LS algorithm performs the best and the linear LS algorithm performs the worst and close to the gross LS algorithm. As more observations are involved, the truncated LS algorithm converges to a level much higher than that can be achieved by the gross LS algorithm and the linear LS algorithm. Therefore, ignoring the term ${\dot{\mathit{\alpha }}}_0$ will cause performance loss for long term observations. The gross LS algorithm is still benefitting from the increase of observations and it is slightly better than the linear LS algorithm for short term observations. The linear LS algorithm can reach a lower RMSE level but will still suffer performance degradation due to the linear model mismatch. With about 40 snapshots of observations, corresponding to 160 observations and $0.8$ s period, the velocity estimation performances of two algorithms will depart.

The gross RMSE of all the algorithms are shown in Figure 3e,f, which have very close appearances to Figure 3c,d. That is because the velocity estimation errors are much greater than the positioning errors. Therefore, although the algorithm can estimate the velocity of targets, the accuracy is low due to a short observation period. In order to estimate the velocity in a higher accuracy, one needs to use observations from a longer period, which can be achieved through a filtering operation.

In order to show the way in which the algorithms converge to the real value, Figure 4a,b are presented to show estimated target positions and velocities at the first 20 snapshots and the latest 20 snapshots, respectively. It can be seen that with a few observations, the linear LS algorithm will converge to the real position of the target to a high accuracy. However, as more observations are available, the gross LS algorithm is closer to the real target position and the linear LS algorithm converges to other locations. Therefore, the gross LS algorithm is more robust in real applications.

3.2. Computation Cost

The advantage of the linear LS algorithm and the truncated LS algorithm lies in its computation cost and communication cost. In applications of the LS algorithms, the locations of the sensor will be approximated by a linear model, which is described by an initial position and a velocity term, with totally 6 parameters. Therefore, it is unnecessary to transmit all observations to the fusion center anymore and thus the communication cost will be reduced. If 100 position estimates are described by 6 parameters, the data to transmit will be reduced to $2\%$. Of course, there is a limit to which the data can be reduced and the limit depends on the platform speed of the sensors, the positions of sensors, the position of the target, and the periods of the observations.

The computation cost reduction, for both the linear LS algorithm and the truncated LS algorithm, stems from reduced number of multiplication and summation operations at the fusion center. The linear LS algorithm and the truncated LS algorithm can be implemented in a structure like parallel computation, namely, the linear regression of the platform position and the local DOA measures are performed at local sensors, and the fusion center just operates on the results of local sensors. To illustrate this fact, we record the computation times of the 20 random experiments for both the algorithms and show the computation times in Figure 5a,b, in scale mode and dB mode respectively. It can be seen that as the number of observations increases, the gross LS algorithm requires a longer computation time, but the linear LS algorithm and the truncated LS algorithm have much plain slopes. Meanwhile, the linear LS algorithm needs more computation cost, as a result of estimating ${\mathbf{v}}_{\mathrm{a}}$ and ${\mathbf{v}}_{\mathrm{e}}$. In fact, the computation cost of the linear LS algorithm does not vary with the sample number too much because it always computes with the same number of parameters, namely the number of sensors N. The computation cost increase due to more observations is imposed over local sensors now.

3.3. The Impact of Velocity Estimation Error

In theoretical derivations, we assumed that the velocities of sensors are estimated by position measures from a device on a platform. In practice, the platform may provide other means to measure the velocity in a higher accuracy. Meanwhile, in order to check whether the performance degradation of the linear LS algorithm in a long period is a result of inaccurate estimation of the platform velocity, we perform a simulation in a way that the estimated velocity ${\widehat{\mathbf{v}}}_{\mathrm{g}}$ is replaced by its real value ${\mathbf{v}}_{\mathrm{g}}$. In this case, there is no velocity error and the only measurement bias is from position measurement. Note subsequently that the linear LS algorithm shares the same position estimate with the truncated LS algorithm.

The RMSE of position estimation is shown in Figure 6. In Figure 6a,b, the sensor location uncertainty is zero-mean normal distributed with ${\mathbf{C}}_{\mathrm{s}}(k,n)=\mathbf{I}$. It can be clearly seen that for the linear LS algorithm, the RMSE curve with real sensor velocity is very close to the RMSE curve using estimated sensor velocity. In order to examine whether a higher position estimation error will make a difference, we make another simulation with ${\mathbf{C}}_{\mathrm{s}}(k,n)=10\mathbf{I}$ and the results are shown in Figure 6c,d. Two RMSE curves are still very close. After some experiments with other position measurement errors, we find that the platform location and velocity regression algorithms can reach a high accuracy and thus will not cause too much performance degradation. This conclusion depends heavily on a fact that the numbers of observations under consideration is often huge according to our configurations.

In fact, from (92), the position estimate of the target does not depend on the velocity estimation of the sensor platform too much. However, there still an insignificant impact, because in our simulation configuration, the sensor location at $t=0$ is obtained by an interpolation operation and if the platform velocity is exactly known a priori, the position estimation will be more accurate.

From (105), the target velocity estimation depends on the sensor velocity more. In order the examine the impact of the sensor velocity estimation on the target velocity estimation performance, we run a simulation with ${\mathbf{C}}_{\mathrm{s}}(k,n)=10\mathbf{I}$ and the results are shown in Figure 7a,b, in scale and dB respectively. It can be seen that it makes a little difference to use real sensor velocity instead of estimated velocity, especially in few earlier observations. As more observations are taken into account, it makes a minor different to replace by real platform velocities. That is because more observations make the velocity estimation more accurate. However, accurate sensor velocity information does not make the target velocity estimation better necessarily, and sometimes, its impact is a bit negative. As the target also moves in a constant velocity, it is reasonable to infer that the linear approximation of the signal DOA has a great impact on the target position and velocity estimation accuracy, as will be analyzed in the subsequent results.

3.4. Nonlinearity of the DOA Approximation

In order to examine the impact of DOA nonlinearity on the final performance, Figure 7c,d show the azimuth angles and elevation angles of the target in the four sensors. The azimuth angle and elevation angle change by about ${15}^{\mathrm{o}}$ at most during 100 snapshots. Over 100 observations, corresponding to 2 s, the nonlinearity of both DOA angles becomes obvious. One should refer to explicit numerical quantities to evaluate the nonlinearity acceptable.

4. Conclusions

This paper studies the target position and velocity estimation problem with distributed passive sensors. The problem is formulated with distributed asynchronous sensors connected to a fusion center with communication links. We first present a gross LS algorithm that takes all angle observations from distributed sensors into account to make a LS estimation. The algorithm is simplified after some matrix manipulations, but as it needs local sensors transmitting all local observations to a fusion center, the communication cost is high. Meanwhile, the computation cost at the fusion center is also high. The communication cost is mainly a result of high-dimensional received data. In order to reduce the communication cost and computation cost, we present a linear LS algorithm that approximates local sensor locations and angle observations with linear models and then estimate target position and velocity with the parameters of the linear models. In order to simplify the velocity estimation, we also present a truncated LS algorithm that just take an average operation to estimate target velocity. In this manner, both the communication cost and the computation cost at the fusion center are reduced significantly. However, the linear LS algorithm and the truncated LS algorithm faces the model mismatching problem, namely, if the linear approximation is not accurate anymore, the performance may degrade greatly. That is a difference from the gross LS algorithm, which always benefits from more observations, as if the linear target position model holds.

The performance of the concerned algorithms is verified with numerical results. It is found that with less observations, the truncated LS algorithm performs the best. As the number of observations increase, the linear LS algorithm and the gross LS algorithm perform better. With more observations available, the linear model mismatch and then the gross LS algorithm perform the best. The gross LS algorithm always benefits from more local observations, which is a difference from the other two algorithms. The cost is a higher communication cost and a higher computation cost at fusion center. We also examined the angle distortion problem that is the only nonlinear term in the simulation configurations. Our matrix operations often make the estimation need less computation costs.

Compared to localization and tracking framework, the algorithms with velocity estimation needs a much lower rate of tracking operations, whose matrix inverse operation often need huge computation cost. Meanwhile, it can provide more accurate measures of target states and the tracking algorithm will also benefit from that. In our simulations, the sensor location error is not taken into account. In practice, this is inevitable. If the self-positioning error is non zero-mean Gaussian distributed, one may incorporate this goal in a distributed angle-only based positioning algorithm, which will be considered in our future works.