Signal Source Positioning Based on Angle-Only Measurements in Passive Sensor Networks

Abstract

Some passive sensors can measure only directions of arrival of signals, but the real positions of signal sources are often desirable, which can be estimated by combining distributed passive sensors as a network. However, passive observations should be correctly associated first. This paper studies the multi-target data association and signal localization problem in distributed passive sensor networks. With angle-only measurements from distributed passive sensors, multiple lines in a 3-dimensional (3D) scenario can be built and then those that will intersect in a small volume in 3D are classified into the same source. The center of the small volume is taken as an estimate of the signal source position, whose statistical distributions are formulated. If the minimum distance is less than an association threshold, then two lines are considered to be from the same signal source. In numerical results, the impacts of angle measurement accuracy and platform self-positioning accuracy are analyzed, indicating that this method can achieve a prescribed data association rate and a high positioning performance with a low computation cost.

1. Introduction

Unlike active sensors such as radars, passive sensors do not transmit signals and thus have no anti-jamming problem [1,2]. However, some passive sensors, such as infrared sensors, photoelectric sensors, electronic counter measurement (ECM) and cameras, can estimate only angles of signal sources. Therefore, their signal source positioning performances are typically poor since they have no accurate range information of signal sources. In order to estimate the positions of signal sources, passive sensors with angle-only observations can be connected with communication links into a network to measure signals sources from different spatial locations. In this case, an algorithm to combine the angle-only observations is needed [3,4,5,6,7]. Compared with the time of arrival (TOA) [8] and the time difference of arrival (TDOA) localization, angle-only localization does not require accurate time synchronization between distributed passive sensors for signal sources with low speeds [6].

In a passive sensor network, there may be multiple signal sources, and before accurate localization, one should first correctly associate observations regarding the same sources [9]. The multi-dimension assignment model is a classical method for data association in passive sensor networks [10,11,12], but it needs the locations of signal sources, which is unavailable before correct data association. A geometry-based localization algorithm for a distributed sensor network is presented in [13], which constructs a test statistic based on the minimal distance between the lines of sight for data association. The measurement errors are considered, but the platform self-positioning errors are not considered.

In this paper, how to perform data association and signal source localization in a 3-dimensional (3D) scenario for distributed passive sensors with only angle measurements is studied. We improve the intersection localization algorithm for passive sensor networks in a multi-target scenario. We first consider the data association problem and then the signal source position estimation problem. The basic concept is to construct a set of lines in a 3D scenario according to angle measurements of signal sources. In data association, measurement lines that will intersect within a small space volume are categorized into the same group. The statistical distribution of the minimal distances of the lines are formulated and the minimum distance between any two observation lines is a random variable, proved to follow the Chi-square distribution. The threshold for correct association is formulated by the misassociation probability; namely, two observations are from the same signal source but are classified into two groups. In the test statistics, not only the measurement errors but also the platform positioning errors are considered, which makes the association performance robust when the platform positioning errors exist.

After data association, observations regarding the same signal sources are grouped, based on which the location of signal sources can be estimated. Three positioning algorithms are considered. It is known that angle measurements are nonlinear functions of coordinates of signal sources. With the Taylor expansion, the least square (LS) algorithm linearizes the nonlinear angle measurements about the target position and then uses the LS method to obtain the target position estimate [6,14,15,16,17,18,19]. In real applications, different passive sensors may obtain observations of different signal-to-noise ratios (SNRs) and then the weighted least squares algorithm (WLS) [6,14] and total least square (TLS) algorithm [20] can be used to obtain a better estimate. Another source location method is the intersection localization algorithm [1,4,5,21,22,23]. The basic concept is that if multiple passive sensors simultaneously measure the signal sources without measurement error, these measurement lines of sight will intersect to the target position. The geometric method and algebraic solution method use this property to estimate the positions of signal sources.

The data association process and target location process of this method are closely combined, which ensures a lower algorithm complexity and a better positioning performance. In numerical results, the improvement of data association and signal-source positioning are analyzed. The impact of the target-sensor geometry on the localization accuracy is also studied, indicating that the localization performance will be better if the lines associated with different observations are perpendicular to each other.

We follow the convention that bold lower and upper case letters denote column vectors and matrices, respectively. A symbol with an upper script o denotes the true value. For instance, ${\mathbf{a}}^o$ denotes the true value of $\mathbf{a}$. $diag(·)$ with a vector entry denotes a diagonal matrix with the entry vector as diagonal elements. The notation diag(${\mathbf{A}}_1,{\mathbf{A}}_2,\dots ,{\mathbf{A}}_N$) stands for the block-diagonal matrix formed by the matrices ${\mathbf{A}}_1,{\mathbf{A}}_2,\dots ,{\mathbf{A}}_N$.

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. Assume that a coordinate system is available for all the sensors, such as earth-centered earth-fixed (ECEF) of the World Geodetic System 84 (WGS84). The real position of the nth sensor at instant t is denoted as ${\mathbf{s}}_n^o\left(t\right)={[x_{n,s}^o\left(t\right),y_{n,s}^o\left(t\right),z_{n,s}^o\left(t\right)]}^{\mathrm{T}}$, $n=1,2,\dots ,N$, where ${(·)}^{\mathrm{T}}$ denotes the transpose operation, t denotes time, and $x_{n,s}^o\left(t\right),y_{n,s}^o\left(t\right),z_{n,s}^o\left(t\right)$ denote the $x,y,z$ coordinates of the nth sensor at instant t, respectively. The real position of the mth target at instant t is denoted by ${\mathbf{g}}_m^o\left(t\right)={[x_{m,g}^o\left(t\right),y_{m,g}^o\left(t\right),z_{m,g}^o\left(t\right)]}^{\mathrm{T}},m=1,\cdots ,M$, where $x_{m,g}^o\left(t\right),y_{m,g}^o\left(t\right),z_{m,g}^o\left(t\right)$ denote the $x,y,z$ coordinates of the mth target at instant t, respectively. Assume that there is a self-positioning device to measure the position of each sensor. There are different kinds of instruments that can measure the position of a platform, such as the Global Positioning System (GPS) and inertial sensors. The topology of the passive sensors and targets are shown in Figure 1.

For the nth sensor, signals are detected 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 instant $t_{k,n}$, assume that the position of the nth sensors is measured as

$${\mathbf{s}}_{k,n}={\mathbf{s}}_n^o\left(t_{k,n}\right)+\Delta {\mathbf{s}}_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 $\Delta {\mathbf{s}}_n\left(t_{k,n}\right)$ denotes the sensor self-positioning error. For simplicity, we assume that the positioning errors follow zero mean Gaussian distributions with covariance matrices ${\mathbf{C}}_{k,n}=\mathbb{E}(\Delta {\mathbf{s}}_n\left(t_{k,n}\right)\Delta {\mathbf{s}}_n^{\mathrm{T}}\left(t_{k,n}\right))$, where $\mathbb{E}$ denotes the expectation operation. In practice, the self-positioning error of the sensor is approximately subject to zero-mean Gaussian distribution, so this assumption is reasonable and widely used.

For the angle-only sensors, the observations are directions of the signals and the lth signal at instant $t_{k,n}$ is denoted by ${\theta }_{l,k,n}={[{\theta }_{l,k,n},{\phi }_{l,k,n}]}^{\mathrm{T}}$, where $(l,k,n)\in {\mathcal{M}}_{k,n}$, $k=1,\cdots ,N_n$ and $n=1,\cdots ,N$, and ${\mathcal{M}}_{k,n}$ denotes a set of triples of signals detected at the kth measurement by the nth sensor. Assume that $|{\mathcal{M}}_{k,n}|=M_{k,n}$, where $|·|$ over a set denotes the cardinality of the set. As the existence of miss detection, false alarms and overlapping of signal sources, $M_{k,n}$ may not be equal to M. Denote ${\mathcal{M}}_n={\cup }_{k=1}^{N_n}{\mathcal{M}}_{k,n}$ and $\mathcal{A}={\cup }_{n=1}^N{\mathcal{M}}_n$, where ∪ denotes the union operation. The total number of observations by N sensors is denoted by

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

(2)

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

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

(3)

For the mth signal source, the real azimuth angle and elevation angle regarding the nth sensor can be expressed by

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

(4)

respectively, where ${\theta }_{m,k,n}^o\in (-\pi ,\pi )$, ${\phi }_{m,k,n}^o\in (-\frac{\pi }{2},\frac{\pi }{2})$, ${tan}^{-1}(*)$ is called the two-argument inverse tangent function [24,25] and $arctan(*)$ is the inverse tangent function. Denote ${\theta }_{m,k,n}^o={[{\theta }_{m,k,n}^o,{\phi }_{m,k,n}^o]}^{\mathrm{T}}$.

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{A}\to \mathcal{M}$. According to our setting, the index set $\mathcal{A}$ can be partitioned into $M+1$ disjoint sets ${\mathcal{A}}_0,{\mathcal{A}}_1,\cdots ,{\mathcal{A}}_M$ defined by

$${\mathcal{A}}_i=\left\{x\right|\psi \left(x\right)=i,x\in \mathcal{A}\},$$

(5)

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.

If the mth signal source is detected and indexed as the $(l,k,n)$ observation, then the azimuth angle and elevation angle measurements can be written as

$$\begin{array}{cc}\hfill {\theta }_{l,k,n}& ={[{\theta }_{l,k,n},{\phi }_{l,k,n}]}^{\mathrm{T}}\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill & ={\theta }_{m,k,n}^o+\Delta {\theta }_{l,k,n}\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill {\theta }_{l,k,n}& ={\theta }_{m,k,n}^o+\Delta {\theta }_{l,k,n}\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill {\phi }_{l,k,n}& ={\phi }_{m,k,n}^o+\Delta {\phi }_{l,k,n}\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill \Delta {\theta }_{l,k,n}& ={[\Delta {\theta }_{l,k,n},\Delta {\theta }_{l,k,n}]}^{\mathrm{T}}\hfill \end{array}$$

(10)

where $\Delta {\theta }_{l,k,n}$ and $\Delta {\phi }_{l,k,n}$ represent the measurement noise of the azimuth angle and elevation angle, respectively. For simplicity, we assume that observation noises $\Delta {\theta }_{l,k,n}$ and $\Delta {\phi }_{l,k,n},l=1,\cdots ,M_{k,n},k=1,\cdots ,N_n,n=1,\cdots ,N$ are statistically independent and follow zero-mean Gaussian distribution by assumption. The covariance matrices of $\Delta \theta$ are denoted by ${\mathbf{C}}_{l,k,n}=\mathbb{E}(\Delta \theta \Delta {\theta }^{\mathrm{T}})\in {\mathbb{C}}^{2\times 2}$, namely $\Delta \theta \sim \mathcal{N}(\mathbf{0},{\mathbf{C}}_{l,k,n})$, which is typically affected by the SNR of the signal, where $\mathcal{N}(\mathbf{0},{\mathbf{C}}_{l,k,n})$ denotes the zero-mean Gaussian distribution with mean $\mathbf{0}$ and covariance matrix ${\mathbf{C}}_{l,k,n}$. It should be noted that modeling the angle measurement noise as a zero-mean Gaussian distribution is a commonplace assumption.

2.2. The Distance of Observation Lines

The angle-only observations provide information on the directions of signal sources, and in theory, the directions can be expressed by 3D lines. Therefore, it is important to study the properties of the lines. Consider a simple scenario where the target can be deemed to be static when the observations are recorded and the location of the mth signal is simplify denoted by ${\mathbf{g}}_m^o$. For passive observations, each observation, say $(l,k,n)$, will contribute a line in space and the line can be written as $L_{l,k,n}$

$$\begin{array}{c}\hfill L_{l,k,n}:\mathbf{x}={\mathbf{s}}_{k,n}+{\alpha }_{l,k,n}{\mathbf{e}}_{l,k,n},{\alpha }_{l,k,n}\in \mathbb{R}\end{array}$$

(11)

where ${\alpha }_{l,k,n}$, a parameter indicating the distance to the origin ${\mathbf{s}}_{k,n}$, ${\mathbf{e}}_{n,m}={[e_{l,k,n,x},e_{l,k,n,y},e_{l,k,n,z}]}^{\mathrm{T}}\in {\mathbb{R}}^{3\times 1}$, is the normalized direction vector associated with the angle observation ${\theta }_{l,k,n}$, and

$$\begin{array}{cc}\hfill e_{l,k,n,x}& =cos\left({\theta }_{l,k,n}\right)cos\left({\phi }_{l,k,n}\right)\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill e_{l,k,n,y}& =sin\left({\theta }_{l,k,n}\right)cos\left({\phi }_{l,k,n}\right)\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill e_{l,k,n,y}& =sin\left({\phi }_{l,k,n}\right).\hfill \end{array}$$

(14)

It can be seen that the subscripts of the denotations are complicated. Therefore, for two observations $i,j\in {\mathcal{A}}_m$, we simplify the expression of lines by

$$\begin{array}{cc}\hfill L_{i,m}:\mathbf{x}& ={\mathbf{s}}_i+{\alpha }_{i,m}{\mathbf{e}}_{i,m}\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill L_{j,m}:{\mathbf{x}}^{\prime }& ={\mathbf{s}}_j+{\alpha }_{j,m}{\mathbf{e}}_{j,m}\hfill \end{array}$$

(16)

where the locations of two sensors regarding the two observations are denoted by ${\mathbf{s}}_i,{\mathbf{s}}_j$, ${\alpha }_{i,m},{\alpha }_{j,m}$ are two scalars indicating the distances to two origins, and ${\mathbf{e}}_{i,m}$ and ${\mathbf{e}}_{j,m}$ are two normalized vectors associated with two observations.

The difference between two points over the two lines are

$$\begin{array}{cc}\hfill \Rightarrow {\mathbf{x}}^{\prime }-\mathbf{x}& ={\mathbf{s}}_j-{\mathbf{s}}_i+{\alpha }_{j,m}{\mathbf{e}}_{j,m}-{\alpha }_{i,m}{\mathbf{e}}_{i,m}\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill & ={\mathbf{s}}_j-{\mathbf{s}}_i+(-{\mathbf{e}}_{i,m},{\mathbf{e}}_{j,m})\alpha ,\hfill \end{array}$$

(18)

where

$$\alpha =\left[\begin{array}{c}{\alpha }_{i,m}\\ {\alpha }_{j,m}\end{array}\right].$$

(19)

Without measurement error, then there will be two ${\alpha }_{i,m}$ and ${\alpha }_{j,m}$ such that ${\mathbf{x}}^{\prime }-\mathbf{x}=0$, i.e.,

$${\mathbf{g}}_m={\mathbf{s}}_i+{\alpha }_{i,m}{\mathbf{e}}_{i,m}={\mathbf{s}}_j+{\alpha }_{j,m}{\mathbf{e}}_{j,m}.$$

(20)

In general, due to inevitable measurement errors, the lines even regarding the same signal source may not coincide to each other. Therefore, we calculate minimal distance between those lines. The distance between two points over two lines can be expressed by

$$\begin{array}{cc}\hfill d=& \parallel {\mathbf{x}}^{\prime }{-\mathbf{x}\parallel }^2\hfill \\ \hfill =& \parallel {\mathbf{s}}_j-{\mathbf{s}}_i+(-{\mathbf{e}}_{i,m},{\mathbf{e}}_{j,m}){\alpha \parallel }^2\hfill \\ \hfill =& {\alpha }^{\mathrm{T}}{\mathbf{E}}_{i,j}^{\mathrm{T}}{\mathbf{E}}_{i,j}\alpha +2{\alpha }^{\mathrm{T}}{\mathbf{E}}_{i,j}^{\mathrm{T}}({\mathbf{s}}_j-{\mathbf{s}}_i)+{({\mathbf{s}}_j-{\mathbf{s}}_i)}^{\mathrm{T}}({\mathbf{s}}_j-{\mathbf{s}}_i)\hfill \end{array}$$

(21)

where $\parallel ·\parallel$ denotes the ${\ell }_2$-norm, $\alpha ={[{\alpha }_{i,m},{\alpha }_{j,m}]}^{\mathrm{T}}$, and ${\mathbf{E}}_{i,j}=[-{\mathbf{e}}_{i,m},{\mathbf{e}}_{j,m}]$. For simplicity, the distance d is actually the squared distance, instead of the distance.

In particular, if ${\mathbf{e}}_{i,m}={\mathbf{e}}_{j,m}$, namely two lines are parallel to each other, then

$$d=\parallel {\mathbf{s}}_j-{\mathbf{s}}_i+({\alpha }_{j,m}-{\alpha }_{i,m}){\mathbf{e}}_{i,m}{\parallel }^2,$$

(22)

and the minimal distance between points in two lines will be achieved if

$${\alpha }_{j,m}-{\alpha }_{i,m}={\mathbf{e}}_{i,m}^{\mathrm{T}}({\mathbf{s}}_i-{\mathbf{s}}_j).$$

(23)

It can be proved that the minimal distance is

$$d=\parallel ({\mathbf{e}}_{i,m}{\mathbf{e}}_{i,m}^{\mathrm{T}}-\mathbf{I})({\mathbf{s}}_i-{\mathbf{s}}_j){\parallel }^2.$$

(24)

If ${\mathbf{e}}_{i,m}\ne {\mathbf{e}}_{j,m}$ or $|{\mathbf{E}}_{i,j}|\ne 0$, the distance is a second order function of ${\alpha }_{i,m}$ and ${\alpha }_{j,m}$, thus the minimal value of d is unique, where $|·|$ over a matrix denotes the determinant of the input matrix. To obtain the minimal value, we take a derivative of d with respect to $\alpha$,

$$\begin{array}{cc}\hfill \frac{dd}{d\alpha }& ={\left[\frac{dd}{d{\alpha }_{i,m}},\frac{dd}{d{\alpha }_{j,m}}\right]}^{\mathrm{T}}\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill & =2{\mathbf{E}}_{i,j}^{\mathrm{T}}{\mathbf{E}}_{i,j}\alpha +2{\mathbf{E}}_{i,j}^{\mathrm{T}}({\mathbf{s}}_j-{\mathbf{s}}_i)\hfill \end{array}$$

(26)

where

$$\begin{array}{cc}\hfill \frac{dd}{d{\alpha }_{i,m}}& =-2{\mathbf{e}}_{i,m}^{\mathrm{T}}({\mathbf{s}}_j-{\mathbf{s}}_i+\mathbf{E}\alpha )\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill \frac{dd}{d{\alpha }_{j,m}}& =2{\mathbf{e}}_{j,m}^{\mathrm{T}}({\mathbf{s}}_j-{\mathbf{s}}_i+\mathbf{E}\alpha ).\hfill \end{array}$$

(28)

Let the derivative be zero and then we obtain a solution

$$\frac{dd}{d\alpha }=0\Rightarrow {\mathbf{E}}_{i,j}^{\mathrm{T}}{\mathbf{E}}_{i,j}\alpha =-{\mathbf{E}}_{i,j}^{\mathrm{T}}({\mathbf{s}}_j-{\mathbf{s}}_i).$$

(29)

Under the assumption that $|{\mathbf{E}}_{i,j}|\ne 0$, ${\mathbf{R}}_{i,j}={\mathbf{E}}_{i,j}^{\mathrm{T}}{\mathbf{E}}_{i,j}$ has a reverse matrix and then the solution can be immediately obtained by

$${\alpha }_{\mathrm{opt}}=-{\left({\mathbf{E}}_{i,j}^{\mathrm{T}}{\mathbf{E}}_{i,j}\right)}^{-1}{\mathbf{E}}_{i,j}^{\mathrm{T}}({\mathbf{s}}_j-{\mathbf{s}}_i).$$

(30)

The minimal distance can be expressed by

$$\begin{array}{cc}\hfill d_{min}& ={({\mathbf{s}}_j-{\mathbf{s}}_i)}^{\mathrm{T}}(\mathbf{I}-{\mathbf{E}}_{i,j}{\mathbf{R}}_{i,j}^{-1}{\mathbf{E}}_{i,j}^{\mathrm{T}})({\mathbf{s}}_j-{\mathbf{s}}_i)\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill & =\parallel (\mathbf{I}-{\mathbf{E}}_{i,j}{\mathbf{R}}_{i,j}^{-1}{\mathbf{E}}_{i,j}^{\mathrm{T}})({\mathbf{s}}_j-{\mathbf{s}}_i){\parallel }^2.\hfill \end{array}$$

(32)

In particular, if

$${\mathbf{E}}_{i,j}{\mathbf{R}}_{i,j}^{-1}{\mathbf{E}}_{i,j}^{\mathrm{T}}({\mathbf{s}}_j-{\mathbf{s}}_i)={\mathbf{s}}_j-{\mathbf{s}}_i,$$

(33)

namely, ${\mathbf{s}}_j-{\mathbf{s}}_i$ is an eigenvector of ${\mathbf{E}}_{i,j}{\mathbf{R}}_{i,j}^{-1}{\mathbf{E}}_{i,j}^{\mathrm{T}}$ and the eigenvalue is 1, then

$$d_{min}=0.$$

(34)

It means that two lines will intersect at a point. Without measurement errors, observations regarding the same target will form lines intersecting at the target position.

In practice, the real mapping $\psi$ should be estimated through an association algorithm. Due to measurement errors, the estimated mapping may not be correct, and then, the positioning error may raise. Therefore, an accurate data association method is important, which will be studied subsequently.

3. Data Association Based on Minimal Distance

3.1. Data Association Model

In order to perform data association, we first need to build the statistical model of the minimal distance of the observation lines. Because of both the platform positioning error $\Delta {\mathbf{s}}_i,\Delta {\mathbf{s}}_j$ and the angle measurement error $\Delta {\theta }_{i,m}$ and $\Delta {\theta }_{j,m},i,j\in {\mathcal{A}}_m$, the minimal distance $d_{min}$ is not zero and follows a certain distribution depending on the measurement errors, where $\Delta {\mathbf{s}}_i$, $\Delta {\mathbf{s}}_j$ denote the sensor self-positioning errors of the ith and the jth observations, respectively, and $\Delta {\theta }_{i,m},\Delta {\theta }_{j,m}$ denote the angle measurement error of the ith and the jth observations on the mth target.

In practice, the ith and the jth observations may belong to the same target or not, and the data association problem is treated as a test to make a decision here. For $i,j\in \mathcal{A}$, the data association problem for a target can be formulated as the following hypothesis problem

$$\left\{\begin{array}{cc}{H}_0:d_{min}=0\hfill & \psi \left(i\right)=\psi \left(j\right)\hfill \\ H_1:d_{min}>0\hfill & \psi \left(i\right)\ne \psi \left(j\right).\hfill \end{array}\right.$$

(35)

To determine the statistical distribution of $d_{min}$ under the $H_0$ hypothesis, we first define

$${\mathbf{K}}_{i,j}=\mathbf{I}-{\mathbf{E}}_{i,j}{\mathbf{R}}_{i,j}^{-1}{\mathbf{E}}_{i,j}^{\mathrm{T}}$$

(36)

and then

$$d_{min}={\parallel {\mathbf{K}}_{i,j}({\mathbf{s}}_j-{\mathbf{s}}_i)\parallel }^2.$$

(37)

Next, we discuss the eigenvalues of ${\mathbf{K}}_{i,j}$. It can be proved that ${\mathbf{K}}_{i,j}^{\mathrm{T}}{\mathbf{K}}_{i,j}={\mathbf{K}}_{i,j}$ and ${\mathbf{K}}_{i,j}^{\mathrm{T}}={\mathbf{K}}_{i,j}$. Therefore, ${\mathbf{K}}_{i,j}$ is a positive semi-definite matrix and its possible eigenvalues are either 0 or 1. As ${\mathbf{K}}_{i,j}$ is not an all-zero matrix, then there is at least an eigenvalue of 1.

The trace of the matrix $\mathbf{K}\left(i,j\right)$ satisfies

$$\begin{array}{cc}\hfill \mathrm{tr}\left({\mathbf{K}}_{i,j}\right)& =\mathrm{tr}(\mathbf{I}-{\mathbf{E}}_{i,j}{\mathbf{R}}_{i,j}^{-1}{\mathbf{E}}_{i,j}^{\mathrm{T}})\hfill \end{array}$$

(38)

$$\begin{array}{cc}\hfill & =\mathrm{tr}\left(\mathbf{I}\right)-\mathrm{tr}\left({\mathbf{E}}_{i,j}{\mathbf{R}}_{i,j}^{-1}{\mathbf{E}}_{i,j}^{\mathrm{T}}\right)\hfill \end{array}$$

(39)

$$\begin{array}{cc}\hfill & =3-\mathrm{tr}\left({\mathbf{E}}_{i,j}^{\mathrm{T}}{\mathbf{E}}_{i,j}{\mathbf{R}}_{i,j}^{-1}\right)\hfill \end{array}$$

(40)

$$\begin{array}{cc}\hfill & =3-\mathrm{tr}\left({\mathbf{R}}_{i,j}{\mathbf{R}}_{i,j}^{-1}\right)\hfill \end{array}$$

(41)

$$\begin{array}{cc}\hfill & =3-\mathrm{tr}({\mathbf{I}}_2\in {\mathbb{R}}^{2\times 2})=1\hfill \end{array}$$

(42)

where $\mathrm{tr}(·)$ with a matrix input denotes the trace of the matrix.

Therefore, it can be inferred that three eigvenvalues of ${\mathbf{K}}_{i,j}$ are $1,0,0$ and the rank of ${\mathbf{K}}_{i,j}$ is 1. In other words, ${\mathbf{K}}_{i,j}$ can be written as

$${\mathbf{K}}_{i,j}={\mathbf{e}}_{\mathrm{s}}{\mathbf{e}}_{\mathrm{s}}^{\mathrm{T}}$$

(43)

where ${\mathbf{e}}_{\mathrm{s}}$ is a unity direction vector perpendicular to both ${\mathbf{e}}_i$ and ${\mathbf{e}}_j$, defined by

$${\mathbf{e}}_{\mathrm{s}}={\mathbf{e}}_i\times {\mathbf{e}}_j/\sqrt{1-{\left({\mathbf{e}}_i^{\mathrm{T}}{\mathbf{e}}_j\right)}^2}$$

(44)

and ${\mathbf{e}}_i$ and ${\mathbf{e}}_j$ denote the unity direction vectors associated with the ith and the jth observations, respectively.

With this fact, the minimal distance can be rewritten as

$$d_{min}={\left|{\mathbf{e}}_{\mathrm{s}}^{\mathrm{T}}({\mathbf{s}}_j-{\mathbf{s}}_i)\right|}^2.$$

(45)

In practice, both ${\mathbf{s}}_j-{\mathbf{s}}_i$ and ${\mathbf{e}}_{\mathrm{s}}$ may be inaccurate. A statistical distribution is necessary to determine the impact of the measurement errors. In order to formulate the statistical distribution of $d_{min}$ under the $H_0$ hypothesis, we define

$$r={\mathbf{e}}_{\mathrm{s}}^{\mathrm{T}}({\mathbf{s}}_j-{\mathbf{s}}_i)$$

(46)

and its relationship to $d_{min}$ is $d_{min}={\left|r\right|}^2$.

The minimal distance varies with a total of 10 parameters; namely ${\theta }_{i,m},{\theta }_{j,m},{\mathbf{s}}_{i,m}$ and ${\mathbf{s}}_{j,m}$, and the corresponding partial derivatives can be written as

$$\begin{array}{cc}\hfill \frac{\partial r}{\partial {\theta }_{i,m}}& ={\left[\frac{\partial r}{\partial {\theta }_{i,m}},\frac{\partial r}{\partial {\phi }_{i,m}}\right]}^{\mathrm{T}}\hfill \end{array}$$

(47)

$$\begin{array}{cc}\hfill \frac{\partial r}{\partial {\theta }_{j,m}}& ={\left[\frac{\partial r}{\partial {\theta }_{j,m}},\frac{\partial r}{\partial {\phi }_{j,m}}\right]}^{\mathrm{T}}\hfill \end{array}$$

(48)

$$\begin{array}{cc}\hfill \frac{\partial r}{\partial {\mathbf{s}}_i}& =-{\mathbf{e}}_{\mathrm{s}}\hfill \end{array}$$

(49)

$$\begin{array}{cc}\hfill \frac{\partial r}{\partial {\mathbf{s}}_j}& ={\mathbf{e}}_{\mathrm{s}}.\hfill \end{array}$$

(50)

Consequently, under the $H_0$ hypothesis, with denotation $\mathbf{v}={[{\theta }_{i,m}^{\mathrm{T}},{\theta }_{j,m}^{\mathrm{T}},{\mathbf{s}}_{i,m}^{\mathrm{T}},{\mathbf{s}}_{j,m}^{\mathrm{T}}]}^{\mathrm{T}}$, we have an approximation

$$\begin{array}{cc}\hfill r& \approx \frac{\partial r}{\partial {\theta }_{i,m}^{\mathrm{T}}}\Delta {\theta }_{i,m}+\frac{\partial r}{\partial {\theta }_{j,m}^{\mathrm{T}}}\Delta {\theta }_{j,m}+\frac{\partial r}{\partial {\mathbf{s}}_i^{\mathrm{T}}}\Delta {\mathbf{s}}_{i,m}+\frac{\partial r}{\partial {\mathbf{s}}_j^{\mathrm{T}}}\Delta {\mathbf{s}}_{j,m}\hfill \end{array}$$

(51)

$$\begin{array}{cc}\hfill & =\frac{\partial r}{\partial {\mathbf{v}}^{\mathrm{T}}}\Delta \mathbf{v}\hfill \end{array}$$

(52)

where $\Delta \mathbf{v}={[\Delta {\theta }_{i,m}^{\mathrm{T}},\Delta {\theta }_{j,m}^{\mathrm{T}},\Delta {\mathbf{s}}_{i,m}^{\mathrm{T}},\Delta {\mathbf{s}}_{j,m}^{\mathrm{T}}]}^{\mathrm{T}}$, and

$$\frac{\partial r}{\partial {\mathbf{v}}^{\mathrm{T}}}=\left[\frac{\partial r}{\partial {\theta }_{i,m}^{\mathrm{T}}},\frac{\partial r}{\partial {\theta }_{j,m}^{\mathrm{T}}},\frac{\partial r}{\partial {\mathbf{s}}_i^{\mathrm{T}}},\frac{\partial r}{\partial {\mathbf{s}}_j^{\mathrm{T}}}\right].$$

(53)

Under the assumption that the measurement errors $\Delta {\theta }_{i,m},\Delta {\theta }_{j,m},\Delta {\mathbf{s}}_{i,m},\Delta {\mathbf{s}}_{j,m}$ are statistically independent of each other and follow a zero-mean Gaussian distribution, r follows the Gaussian distribution and then $d_{min}$ follows the central weighted Chi-square distribution with 1 degree of freedom, namely $d_{min}\sim {\chi }_1^2\left(\lambda \right)$, where $\lambda$ is the variance

$$\lambda =\mathbb{E}(\frac{\partial r}{\partial {\mathbf{v}}^{\mathrm{T}}}\Delta \mathbf{v}\Delta {\mathbf{v}}^{\mathrm{T}}\frac{\partial r}{\partial \mathbf{v}})=\frac{\partial r}{\partial {\mathbf{v}}^{\mathrm{T}}}{\mathbf{C}}_{\mathbf{v}}\frac{\partial r}{\partial \mathbf{v}}$$

(54)

and

$${\mathbf{C}}_{\mathbf{v}}=\mathbb{E}(\Delta \mathbf{v}\Delta {\mathbf{v}}^{\mathrm{T}}).$$

(55)

The probability density function (PDF) and cumulative distribution function (CDF) can be written as

$$\begin{array}{cc}\hfill p_{d_{min}}\left(d\right)& =\frac{1}{2\lambda \sqrt{\pi }}{\left(\frac{d}{\lambda }\right)}^{-\frac{1}{2}}exp(-\frac{d}{2\lambda }),d\ge 0\hfill \end{array}$$

(56)

$$\begin{array}{cc}\hfill F_{d_{min}}\left(d\right)& =\frac{1}{\sqrt{\pi }}\gamma (\frac{d}{2\lambda },\frac{1}{2}),d\ge 0\hfill \end{array}$$

(57)

respectively, where $\Gamma (·)$ denotes the Gamma function and $\gamma (·,·)$ denotes the incomplete Gamma function. If the decision rule is to keep the misassociation probability $P\left(H_1\right|H_0)$ as a constant, say $p_{\mathrm{f}}$, then the decision threshold can be obtained as

$$\begin{array}{c}\hfill \rho =2\lambda {\gamma }^{-1}(\sqrt{\pi }(1-p_{\mathrm{f}}),\frac{1}{2})\end{array}$$

(58)

where ${\gamma }^{-1}(·,\frac{1}{2})$ denotes the inverse incomplete Gamma function with $\frac{1}{2}$ a degree of freedom.

Therefore, the key is to derive the variance $\lambda$. In practice, as the measurements are carried out by different sensors, it is reasonable to assume that the measurement errors are statistically independent of each other. In this case, ${\mathbf{C}}_{\mathbf{v}}$ is a block diagonal matrix and the variance can be formulated conveniently.

3.2. Measurement Errors and Association Threshold

For simplicity, we first consider the self-positioning error $\Delta {\mathbf{s}}_i,\Delta {\mathbf{s}}_j$, whose covariance matrices are assumed to be

$$\begin{array}{cc}\hfill {\mathbf{C}}_{i,s}& =\mathbb{E}(\Delta {\mathbf{s}}_i\Delta {\mathbf{s}}_i^{\mathrm{T}})\hfill \end{array}$$

(59)

$$\begin{array}{cc}\hfill {\mathbf{C}}_{j,s}& =\mathbb{E}(\Delta {\mathbf{s}}_j\Delta {\mathbf{s}}_j^{\mathrm{T}}).\hfill \end{array}$$

(60)

Consequently, from (49) and (50), the variances due to the two terms are

$$\begin{array}{cc}\hfill {\lambda }_{i,s}& =\mathbb{E}(\frac{\partial r}{\partial {\mathbf{s}}_i}\Delta {\mathbf{s}}_{i,m}\Delta {\mathbf{s}}_{i,m}^{\mathrm{T}}{\left(\frac{\partial r}{\partial {\mathbf{s}}_i}\right)}^{\mathrm{T}})={\mathbf{e}}_{\mathrm{s}}^{\mathrm{T}}{\mathbf{C}}_{i,s}{\mathbf{e}}_{\mathrm{s}}\hfill \end{array}$$

(61)

$$\begin{array}{cc}\hfill {\lambda }_{j,s}& =\mathbb{E}(\frac{\partial r}{\partial {\mathbf{s}}_j}\Delta {\mathbf{s}}_{j,m}\Delta {\mathbf{s}}_{j,m}^{\mathrm{T}}{\left(\frac{\partial r}{\partial {\mathbf{s}}_j}\right)}^{\mathrm{T}})={\mathbf{e}}_{\mathrm{s}}^{\mathrm{T}}{\mathbf{C}}_{j,s}{\mathbf{e}}_{\mathrm{s}}.\hfill \end{array}$$

(62)

Next, we consider the variance caused by the angle measurement error $\Delta {\theta }_{i,m}$ and $\Delta {\theta }_{j,m}$. Assume that the covariance matrices of $\Delta {\theta }_{i,m}$ and $\Delta {\theta }_{j,m}$ are

$$\begin{array}{cc}\hfill {\mathbf{C}}_{i,\theta }& =\mathbb{E}(\Delta {\theta }_{i,m}\Delta {\theta }_{i,m}^{\mathrm{T}})\hfill \end{array}$$

(63)

$$\begin{array}{cc}\hfill {\mathbf{C}}_{j,\theta }& =\mathbb{E}(\Delta {\theta }_{j,m}\Delta {\theta }_{j,m}^{\mathrm{T}})\hfill \end{array}$$

(64)

respectively. The variances caused by the two terms are

$$\begin{array}{cc}\hfill {\lambda }_{i,\theta }=& \mathbb{E}(\frac{\partial r}{\partial {\theta }_{i,m}^{\mathrm{T}}}\Delta {\theta }_{i,m}\Delta {\theta }_{i,m}^{\mathrm{T}}\frac{\partial r}{\partial {\theta }_{i,m}})=\frac{\partial r}{\partial {\theta }_{i,m}^{\mathrm{T}}}\mathbb{E}(\Delta {\theta }_{i,m}\Delta {\theta }_{i,m}^{\mathrm{T}})\frac{\partial r}{\partial {\theta }_{i,m}}=\frac{\partial r}{\partial {\theta }_{i,m}^{\mathrm{T}}}{\mathbf{C}}_{i,\theta }\frac{\partial r}{\partial {\theta }_{i,m}}\hfill \end{array}$$

(65)

$$\begin{array}{cc}\hfill {\lambda }_{j,\theta }=& \mathbb{E}(\frac{\partial r}{\partial {\theta }_{j,m}^{\mathrm{T}}}\Delta {\theta }_{j,m}\Delta {\theta }_{j,m}^{\mathrm{T}}\frac{\partial r}{\partial {\theta }_{j,m}})=\frac{\partial r}{\partial {\theta }_{j,m}^{\mathrm{T}}}\mathbb{E}(\Delta {\theta }_{j,m}\Delta {\theta }_{j,m}^{\mathrm{T}})\frac{\partial r}{\partial {\theta }_{j,m}}=\frac{\partial r}{\partial {\theta }_{j,m}^{\mathrm{T}}}{\mathbf{C}}_{j,\theta }\frac{\partial r}{\partial {\theta }_{j,m}}\hfill \end{array}$$

(66)

It can be seen that the distance is not a linear function of $\theta$. As a cross product of ${\mathbf{e}}_{i,m}$ and ${\mathbf{e}}_{j,m}$, ${\mathbf{e}}_{\mathrm{s}}$ is a complicated function. From (47) and (48), $\frac{\partial r}{\partial {\theta }_{i,m}}$ and $\frac{\partial r}{\partial {\theta }_{j,m}}$ are proved in Appendix A.

Under these assumptions, ${\mathbf{C}}_{\mathrm{v}}=diag({\mathbf{C}}_{i,\theta },{\mathbf{C}}_{j,\theta },{\mathbf{C}}_{i,s},{\mathbf{C}}_{i,s})$. Consequently, the variance $\lambda$ can be expressed by

$$\lambda ={\lambda }_{i,s}+{\lambda }_{j,s}+{\lambda }_{i,\theta }+{\lambda }_{j,\theta }.$$

(67)

It should be noted that under the $H_1$ hypothesis, the distance may be arbitrary, and for simplicity, we assume that the distance is uniformly distributed over the surveillance volume. In this case, we can determine whether two observations are from the same target by the following decision rule

$$\begin{array}{c}\hfill d_{min}\underset{H_0}{\overset{H_1}{\gtrless }}\rho .\end{array}$$

(68)

The association algorithm is a mapping ${\psi }^{\prime }:\mathcal{A}\to \mathcal{M}$, which partitions $\mathcal{A}$ into $M+1$ groups and this mapping may disagree with real $\psi$ due to inevitable errors. In practice, some knowledge about the targets of interest may be available and thus can be used for better performance. For instance, if only targets on the ground are of interest, then we may use this information and discard observations in which the cross points are obviously apart from the ground, which will be studied in the future.

3.3. Localization Algorithms

With mapping ${\psi }^{\prime }$, we obtain another partition $\mathcal{A}={\cup }_{m=0}^M{\mathcal{A}}_m^{\prime }$, where ${\mathcal{A}}_m^{\prime }=\left\{x\right|{\psi }^{\prime }\left(x\right)=m,x\in \mathcal{A}\}$. For observations in ${\mathcal{A}}_m^{\prime }$, we can conduct signal source positioning, and three target positioning estimation methods will be introduced subsequently.

We first consider the intersection method. In presence of measurement errors, we often have $d_{min}\ne 0$ even if two observations are from the same target. In this case, over two lines, the two points with minimal distance are

$$\begin{array}{cc}\hfill {\mathbf{x}}_i& ={\mathbf{s}}_i+{\alpha }_{\mathrm{opt}}\left(1\right){\mathbf{e}}_{i,m}\hfill \end{array}$$

(69)

$$\begin{array}{cc}\hfill {\mathbf{x}}_j& ={\mathbf{s}}_j+{\alpha }_{\mathrm{opt}}\left(2\right){\mathbf{e}}_{j,m}.\hfill \end{array}$$

(70)

In this case, it is reasonable to take the middle of two points as the estimate of the signal position, namely

$$\begin{array}{cc}\hfill {\widehat{\mathbf{g}}}_{i,j}& =\frac{1}{2}({\mathbf{x}}_i+{\mathbf{x}}_j)\hfill \\ \hfill & =\frac{1}{2}({\mathbf{s}}_i+{\mathbf{s}}_j)-\frac{1}{2}({\mathbf{e}}_{i,m},{\mathbf{e}}_{j,m}){\mathbf{R}}_{i,j}^{-1}{\mathbf{E}}_{i,j}^{\mathrm{T}}({\mathbf{s}}_j-{\mathbf{s}}_i).\hfill \end{array}$$

(71)

With more observations available, there will be an estimate of the target location for each observation pair, and a simple estimate of the target location can be expressed by their average, namely

$${\widehat{\mathbf{g}}}_m=\frac{1}{N_m(N_m-1)}\sum _{i,j\in {\mathcal{A}}_m}{\widehat{\mathbf{g}}}_{i,j}$$

(72)

where $N_m$ denotes the number of observations associated with the mth signal source.

In practice, another widely used localization algorithm is the LS algorithm. It stems from the delta method concerned in [26]. From (4), the ith angle observation denoted by ${\theta }_{i,m}$ actually contributes a geometric relationship, which can be expressed by

$${\mathbf{G}}_i{\mathbf{g}}_m={\mathbf{G}}_i{\mathbf{s}}_i,i\in {\mathcal{A}}_m$$

(73)

where

$$\begin{array}{cc}\hfill {\mathbf{G}}_i& =\left[\begin{array}{ccc}sin{\theta }_{i,m}& -cos{\theta }_{i,m}& 0\\ cos{\theta }_{i,m}sin{\phi }_{i,m}& sin{\theta }_{i,m}sin{\phi }_{i,m}& -cos{\phi }_{i,m}\end{array}\right].\hfill \end{array}$$

(74)

In (74), we used the equality

$$\begin{array}{cc}\hfill & (x_{m,g}^o-x_{i,s}^o)cos{\theta }_{i,m}^o+(y_{m,g}^o-y_{i,s}^o)sin{\theta }_{i,m}^o\hfill \end{array}$$

(75)

$$\begin{array}{cc}\hfill =& \sqrt{{(x_{m,g}^o-x_{i,s}^o)}^2+{(y_{m,g}^o-y_{i,s}^o)}^2}\hfill \\ \hfill =& d_{i,m}^{o}cos{\phi }_{i,m}^o\hfill \end{array}$$

(76)

where the value $d_{i,m}^o$ is the true range of ${\mathbf{g}}_m^o$ to ${\mathbf{s}}_i^o$ and can be expressed by $d_{i,m}^o=\parallel {\mathbf{g}}_m^o-{\mathbf{s}}_i^o\parallel$.

Next, we can combine the observations into an equation

$$\mathbf{G}{\mathbf{g}}_m=\mathbf{y}$$

(77)

where

$$\begin{array}{cc}\hfill \mathbf{G}=\left[\begin{array}{c}{\mathbf{G}}_1\\ ⋮\\ {\mathbf{G}}_{M_m}\end{array}\right],\mathbf{y}& =\left[\begin{array}{c}{\mathbf{G}}_1{\mathbf{s}}_1\\ ⋮\\ {\mathbf{G}}_{M_m}{\mathbf{s}}_{M_m}\end{array}\right].\hfill \end{array}$$

(78)

It should be noted that after the data association operation, the number of samples associated with a target may not agree with the real number. For simplicity, we still use $M_m=\left|{\mathcal{A}}_m^{\prime }\right|$ to denote the number of observations associated with the mth target.

With $M_m$ observations, there are in total $2M_m$ equations and three unknown parameters in ${\mathbf{g}}_m$. Therefore, as if $M_m\ge 2$, we can use the LS method to obtain an optimal solution in the sense of the mean square error, as

$$\begin{array}{c}\hfill {\widehat{\mathbf{g}}}_m={\left({\mathbf{G}}^{\mathrm{T}}\mathbf{G}\right)}^{-1}\mathbf{G}\mathbf{y}.\end{array}$$

(79)

In practice, different observations may have different SNRs, and then, the WLS algorithm can also be used in this framework. Then, the angle measurement error and sensor positioning error are equally weighted in the process of the LS algorithm. Assume that the distribution of the angle measurement error and the positioning error of sensors are known a priori. With this information, we can impose different weights over the observations, which is the WLS algorithm.

With the following approximations

$$\begin{array}{cc}\hfill sin{\theta }_{i,m}& \approx sin\left({\theta }_{i,m}^o\right)+\Delta {\theta }_{i,m}^{\mathrm{T}}cos{\theta }_{i,m}^o\hfill \end{array}$$

(80)

$$\begin{array}{cc}\hfill cos{\theta }_{i,m}& \approx cos\left({\theta }_{i,m}^o\right)-\Delta {\theta }_{i,m}^{\mathrm{T}}sin{\theta }_{i,m}^o\hfill \end{array}$$

(81)

we can put (80) and (81) into (73) and then write (73) as

$$\begin{array}{cc}\hfill {ϵ}_i& ={\mathbf{G}}_i({\mathbf{g}}_m^o-{\mathbf{s}}_i)\hfill \end{array}$$

(82)

where

$$\begin{array}{cc}\hfill {ϵ}_i& ={[{ϵ}_{i,1},{ϵ}_{i,2}]}^{\mathrm{T}}\hfill \end{array}$$

(83)

$$\begin{array}{cc}\hfill {ϵ}_{i,1}& =\Delta {\theta }_{i,m}{d}_{i,m}^{o}cos{\phi }_{i,m}^o+{\mathbf{a}}_{i,m}^{o\mathrm{T}}\Delta {\mathbf{s}}_i\hfill \end{array}$$

(84)

$$\begin{array}{cc}\hfill {ϵ}_{i,2}& =\Delta {\phi }_{i,m}{d}_{i,m}^o+{\mathbf{b}}_{i,m}^{o\mathrm{T}}\Delta {\mathbf{s}}_i\hfill \end{array}$$

(85)

$$\begin{array}{cc}\hfill {\mathbf{a}}_{i,m}^o& ={[-sin{\theta }_{i,m}^o,cos{\theta }_{i,m}^o,0]}^{\mathrm{T}}\hfill \end{array}$$

(86)

$$\begin{array}{cc}\hfill {\mathbf{b}}_{i,m}^o& ={[-cos{\theta }_{i,m}^{o}sin{\phi }_{i,m}^o,-sin{\theta }_{i,m}^{o}sin{\phi }_{i,m}^o,cos{\phi }_{i,m}^o]}^{\mathrm{T}}.\hfill \end{array}$$

(87)

In (84), we have used the equality (76). In (85), we have used the equality

$$\begin{array}{cc}\hfill (x_{m,g}^o-x_{i,s}^o)sin{\theta }_{i,m}^o-(y_{m,g}^o-y_{i,s}^o)cos{\theta }_{i,m}^o& =0\hfill \end{array}$$

(88)

$$\begin{array}{cc}\hfill {({\mathbf{g}}_m^o-{\mathbf{s}}_i^o)}^{\mathrm{T}}{[cos{\theta }_{i,m}^{o}cos{\phi }_{i,m}^o,sin{\theta }_{i,m}^{o}cos{\phi }_{i,m}^o,sin{\phi }_{i,m}^o]}^{\mathrm{T}}& =d_{i,m}^o\hfill \end{array}$$

(89)

which can be easily verified in the angle measurement equations.

The vector ${ϵ}_i$ can be rewritten as

$${ϵ}_i={\mathbf{B}}_i^o{\delta }_i$$

(90)

where

$$\begin{array}{cc}\hfill {\mathbf{B}}_i^o& =\left[\begin{array}{ccc}{d}_{i,m}^{o}cos{\phi }_{i,m}^o& 0& {\mathbf{a}}_{i,m}^{o,\mathrm{T}}\\ 0& d_{i,m}^o& {\mathbf{b}}_{i,m}^{o,\mathrm{T}}\end{array}\right]\hfill \end{array}$$

(91)

$$\begin{array}{cc}\hfill {\delta }_i& ={[\Delta {\theta }_{i,m},\Delta {\phi }_{i,m},\Delta {\mathbf{s}}_i^{\mathrm{T}}]}^{\mathrm{T}}.\hfill \end{array}$$

(92)

Under the assumption that the self-positioning error and the angle measurement error are decorrelated, it can be proven that the covariance matrix of ${ϵ}_i$ can be written as

$${\mathbf{Q}}_n=\mathbb{E}\left(ϵ{ϵ}^{\mathrm{T}}\right)=diag({\mathbf{C}}_{n,\theta },{\mathbf{C}}_{n,s}).$$

(93)

Putting (90) into (82) and combining the observations into an equation yield

$$\begin{array}{c}\hfill {\mathbf{B}}^o\delta =\mathbf{G}{\mathbf{g}}_m^o-\mathbf{y}.\end{array}$$

(94)

where

$$\begin{array}{cc}\hfill {\mathbf{B}}^o=\left[\begin{array}{c}{\mathbf{B}}_1^o\\ ⋮\\ {\mathbf{B}}_{M_m}^o\end{array}\right],\delta & =\left[\begin{array}{c}{\delta }_1\\ ⋮\\ {\delta }_{M_m}\end{array}\right].\hfill \end{array}$$

(95)

In the WLS algorithm, the goal is to minimize the objective function $J\left({\mathbf{g}}_m\right)$ as

$$\begin{array}{c}\hfill J\left({\mathbf{g}}_m\right)={(\mathbf{G}{\mathbf{g}}_m-\mathbf{y})}^{\mathrm{T}}\mathbf{W}(\mathbf{G}{\mathbf{g}}_m-\mathbf{y})\end{array}$$

(96)

where $\mathbf{W}$ is the weighting matrix with the expression

$$\begin{array}{c}\hfill \mathbf{W}=\mathbb{E}{\left[{\mathbf{B}}^o\delta {\delta }^{\mathrm{T}}{\mathbf{B}}^{o\mathrm{T}}\right]}^{-1}={\left({\mathbf{B}}^o\mathbf{Q}{\mathbf{B}}^{o\mathrm{T}}\right)}^{-1}\end{array}$$

(97)

and $\mathbf{Q}$ is the error covariance matrix with an expression $\mathbf{Q}=diag({\mathbf{Q}}_1,{\mathbf{Q}}_2,\dots ,{\mathbf{Q}}_N)$.

The variable ${\mathbf{g}}_m$ to minimize the objective function $J\left({\mathbf{g}}_m\right)$ can be calculated by the least square method and the estimate of the target positioning can be expressed by

$$\begin{array}{c}\hfill {\widehat{\mathbf{g}}}_m={\left({\mathbf{G}}^{\mathrm{T}}\mathbf{W}\mathbf{G}\right)}^{-1}{\mathbf{G}}^{\mathrm{T}}\mathbf{W}\mathbf{y}\end{array}$$

(98)

which is the WLS algorithm.

In WLS, the covariance matrices of different observations should be known a priori. If the observations can be deemed to have close SNRs, then one can simply use the LS algorithm. If the covariance matrices are not known with certain accuracy, some performance loss may occur, which should be analyzed with numerical analysis.

4. Numerical Results

We first study the localization performance in the presence of a platform of self-positioning errors and angle measurement errors. Then, the data association performance will be analyzed, followed by the analysis of the impact of sensor-target geometry.

Consider a scenario where two sensors are installed on two aircraft and three targets of interest are in the scope. Assume that two aircraft move at the same speed. The parameter settings are shown in

Table 1. Positions and velocities of sensors and targets.

	Position (m) at $\mathit{t}=0$ s	Velocity (m/s)	Position (m) at $\mathit{t}=10$ s
Sensor #1	${\mathbf{s}}_1^{\mathrm{o}}\left(0\right)={[0,0,0]}^{\mathrm{T}}$	${[50,100,0]}^{\mathrm{T}}$	${\mathbf{s}}_1^{\mathrm{o}}\left(10\right)={[500,1000,0]}^{\mathrm{T}}$
Sensor #2	${\mathbf{s}}_2^{\mathrm{o}}\left(0\right)=$ [12,000, 10,000, −800${]}^{\mathrm{T}}$	${[50,100,0]}^{\mathrm{T}}$	${\mathbf{s}}_2^{\mathrm{o}}\left(10\right)=$ [12,500, 11,000, −800${]}^{\mathrm{T}}$
Target #1	${\mathbf{g}}_1^{\mathrm{o}}\left(0\right)=$ [18,000, 12,000, 8000]T	${[20,30,0]}^{\mathrm{T}}$	${\mathbf{g}}_1^{\mathrm{o}}\left(10\right)=$ [18,200, 12,300, 8000${]}^{\mathrm{T}}$
Target #2	${\mathbf{g}}_2^{\mathrm{o}}\left(0\right)=$ [15,000, 13,000, 7000${]}^{\mathrm{T}}$	${[20,30,0]}^{\mathrm{T}}$	${\mathbf{g}}_2^{\mathrm{o}}\left(10\right)=$ [15,200, 13,300, 7000${]}^{\mathrm{T}}$
Target #3	${\mathbf{g}}_3^{\mathrm{o}}\left(0\right)=$ [13,000, 12,000, 5000${]}^{\mathrm{T}}$	${[20,30,0]}^{\mathrm{T}}$	${\mathbf{g}}_3^{\mathrm{o}}\left(10\right)=$ [13,200, 12,300, 5000${]}^{\mathrm{T}}$

. For simplicity, we assume that two passive sensors collect their observations on the same instants, namely, $t_{k,i}=t_{k,j}$ for all k and $i,j\in \{1,\cdots ,N\}$. Meanwhile, the sampling interval is $0.1$ s and the simulation runs for 10 s. Assume that the covariance matrices of angle measurement error for all sensors and all targets are the same as ${\mathbf{C}}_{k,\theta }={\sigma }_{\theta }^2{\mathbf{I}}_2,\forall k$. The self-positioning error covariance is ${\mathbf{C}}_{n,s}={\sigma }_{\mathrm{s}}^2{\mathbf{I}}_3,n=1,\cdots ,N$.

4.1. Impact of Self-Positioning Error and Angle Measurement Error

We first study the impacts of sensor positioning error and angle measurement error in the proposed intersection method. In order to make a comparison between the intersection method, the LS algorithm [19] and the WLS algorithm [6], we first consider the positioning accuracy of target #1 with the two sensors at $t=0$ s.

In order to analyze the positioning performance of the concerned algorithms, we measure the performance with the average root mean square error (RMSE), defined for a target, e.g., target #1, by

$$\mathrm{RMSE}=\sqrt{\frac{1}{N_{\mathrm{s}}}\sum _{k=1}^{N_{\mathrm{s}}}{\parallel {\widehat{\mathbf{g}}}_1\left(k\right)-{\mathbf{g}}_1^{\mathrm{o}}\parallel }^2}$$

(99)

where ${\widehat{\mathbf{g}}}_1\left(k\right)$ denotes the estimate of the target position at the kth simulation run for the target #1, and $N_{\mathrm{s}}$ denotes the number of simulation runs.

With $N_{\mathrm{s}}=$ 10,000 Monte Carlo runs, the RMSE of different self-positioning errors are shown in Figure 2, under the assumption that the covariance matrix is ${\mathbf{C}}_{1,\mathrm{s}}={\sigma }_{\mathrm{s}}^2{\mathbf{I}}_3$, where we always set ${\sigma }_{\theta }=0.{03}^{\circ }$ in this simulation. It can be seen that as ${\sigma }_{\mathrm{s}}$ rises from 0 to 20 m, the three algorithms have close performances and the WLS performs better, under the assumption that the covariance matrices of ${\mathbf{C}}_{k,\theta },k=1,2,3$ and ${\mathbf{C}}_{n,\mathrm{s}},n=1,2$ are known exactly.

The impacts of the angle measurement error are shown in Figure 3, where ${\sigma }_{\theta }$ raises from $0^{\circ }$ to $1.5^{\circ }$, where we set ${\sigma }_{\mathrm{s}}=5$ m in this simulation. It shows that the three algorithms have close performances and the intersection method and the LS method perform a little worse than the WLS method. The rise of both the self-positioning error and the angle measurement error will cause the rise of the target positioning error. However, with a better weighting, the WLS algorithm often performs better than the intersection method and the LS method.

4.2. Data Association Performance

Consider the parameters of the two sensors and the targets, as shown in Table 1. As the targets are moving in this scenario, we set the sampling interval to $0.1$ s. At each sampling instant, each sensor has three measurements corresponding to the three targets. Therefore, the two sensors will totally have nine measurement combinations at each sampling instant, of which three combinations are correct. In data association, we set $p_{\mathrm{f}}=1\%$ and then the probability of correct association $p_{\mathrm{c}}=99\%$.

With $N_s=2000$ Monte Carlo runs, the averaged correct association probability is shown in Figure 4, where the self-positioning error is ${\sigma }_{\mathrm{s}}=5$ m, and the angle measurement error is ${\sigma }_{\theta }=0.{03}^{\circ }$. It can be seen that during the whole sampling period, the correct association probability of the three targets is close to the present value $99\%$, namely 20 wrong combinations, on average, for each instant.

Next, we consider the location performance with observations after the data association operation. The RMSEs of the algorithms over the three targets based on the proposed method are shown in Figure 5. In Figure 5, the RMSEs of the three targets are $23.43$, $30.05$ and $14.45$ m. Note that the data association error may affect the localization performance in this case. From Figure 5, the relative position of the sensor and the target will affect the positioning performance, so the impact of the geometry of sensors and targets on the positioning performance will be studied subsequently.

4.3. Impact of Target-Sensor Geometry

The geometry will play an important role in the positioning performance. Consider the scenario shown in Figure 6, where a target is probed with two passive sensors and the intersection angle of two azimuth lines is denoted by $\varphi$. In the 3D scenario, we also define a new elevation angle $\eta$ for convenience and the distance of the target is the same for both the sensors, i.e., 8000 m.

It is assumed that the positions of the two sensors are accurately measured and target #1 can be detected by both two sensors. The two sensors are symmetrically distributed on both sides of target #1. We explore the impact of geometry on the positioning accuracy by changing $\varphi$ and $\eta$. In order to illustrate that the errors caused by geometry on the three coordinate axes are different, we define $\Delta x$, $\Delta y$ and $\Delta z$ to represent the RMSE on the three coordinate axes, respectively, which can be expressed by

$$\begin{array}{c}\hfill \Delta x=\sqrt{\frac{1}{N_{\mathrm{s}}}\sum _{k=1}^{N_{\mathrm{s}}}{\parallel {\widehat{x}}_{1,g}\left(k\right)-x_{1,g}^{\mathrm{o}}\parallel }^2}\end{array}$$

(100)

$$\begin{array}{c}\hfill \Delta y=\sqrt{\frac{1}{N_{\mathrm{s}}}\sum _{k=1}^{N_{\mathrm{s}}}{\parallel {\widehat{y}}_{1,g}\left(k\right)-y_{1,g}^{\mathrm{o}}\parallel }^2}\end{array}$$

(101)

$$\begin{array}{c}\hfill \Delta z=\sqrt{\frac{1}{N_{\mathrm{s}}}\sum _{k=1}^{N_{\mathrm{s}}}{\parallel {\widehat{z}}_{1,g}\left(k\right)-z_{1,g}^{\mathrm{o}}\parallel }^2}.\end{array}$$

(102)

With $N_s$ = 10,000 Monte Carlo runs, the $\Delta x$, $\Delta y$, $\Delta z$ and the spatial error $sum$ denoted by $sum=\sqrt{\Delta x^2+\Delta y^2+\Delta z^2}$ are shown in Figure 7, under the assumption that $\eta ={15}^{\circ }$, and the angle measurement error ${\sigma }_{\theta }=0.{03}^{\circ }$. It can be seen that when the intersection angle $\varphi$ is changed, the positioning error $\Delta x$, $\Delta y$, $\Delta z$ are different. When the intersection angle $\varphi$ is equal to $82.{58}^{\circ }$, the spatial error $sum$ is the best, about $5.58$ m, and $\Delta x$, $\Delta y$ and $\Delta z$ are approximately equal to each other. In fact, the angle around $\varphi ={90}^{\mathrm{o}}$ will all lead to high accuracy.

The impacts of intersection $\varphi$ and $\eta$ on spatial error $sum$ are shown in Figure 8. For $\eta$ in the range of $0^{\circ }$ to ${45}^{\circ }$ and $\varphi$ in the range of ${25}^{\circ }$ to ${175}^{\circ }$, the spatial error sum is less than 20 m. Therefore, in practice, one can look for geometry with sensors and targets nearly perpendicular to each other to improve the positioning performance.

5. Conclusions

This paper studies the data association and signal source localization problems with distributed passive sensors with angle-only observations. A geometry-based data association method is considered, and the concept is that real targets will contribute observations with a small minimal distance. The statistical distribution of the minimal distance of two lines associated with the same target is formulated, based on which a data association method based on hypothesis testing is also developed. The decision threshold is formulated. Meanwhile, for observations that are classified into the same class, three positioning algorithms are studied, namely the intersection method, the LS method and the WLS method. Two kinds of measurements errors are considered, namely sensor self-positioning error and angle measurement error.

In numerical results, we analyze the data association performance of the concerned positioning algorithms and the signal source positioning performance in different scenarios, indicating that the data association algorithm works well and the positioning performances of the algorithms are very close to each other. Meanwhile, if the observation lines are approximately perpendicular to each other, then the localization performance is more accurate.

This algorithm can be used in laser, infrared, and other passive sensors with angle-only measurements. In practice, other information, such as range, ground surface and sea surface, may be available, which can be incorporated into the positioning and association algorithms to improve performance.