Cooperative Target Localization and Tracking with Incomplete Measurements

Abstract

This study investigates a problem on target localization and tracking for two cases where either the slant range information of dual stations is lost or the slant range information of one station and the pitch angle information of the other one are missing. The models of cooperative localization with incomplete measurements are presented and the Kalman filtering algorithm is applied for target tracking. For improving tracking precision, a strategy of observers path planning based on the gradient of circular error probability (CEP) is integrated into the Kalman filtering algorithm. Several numerical examples are used to illustrate the tracking performance of the proposed algorithm with the corresponding root mean square error (RMSE) and Cramer-Rao lower bound (CRLB). The Monte Carlo simulation results validate the effectiveness of the presented algorithm.

1. Introduction

In practical target detection systems, measurement signals are usually subject to missing information due to a variety of reasons, such as the high maneuverability of the tracked target, a fault in the measurement, intermittent sensor failures, network congestion, accidental loss of some collected data, or the fact that some of the data may be jammed or coming from a very noisy environment. This missing measurement phenomenon is usually called incomplete target information or incomplete measurement [1–3]. Such an imperfect target detection attracted considerable attention during the past few years and various approaches have been presented in the literature to model this incomplete measurement phenomenon and many filtering algorithms and strategies are proposed to handle this problem; see [4–7] and the references cited therein. However, the above mentioned target localization and tracking methods for incomplete measurements have been investigated only for the case of single observation station while the case of dual or multiple stations with incomplete measurements is rarely considered. In fact, such methods for target localization and tracking using single observation station with incomplete measurements may have the problem that the system is unobservable. This can be illustrated as the following special situation: if the slant range measurement is missing and the single observation station is still, the localization and tracking system is unobservable. In such situation, only when the single observation station is moving along an appropriate path that the system can be observable. Therefore, this paper is dedicated to handling the problem of imperfect target detection with dual observation stations by cooperative localization or cross localization methods.

The cooperative localization or cross localization [8–10] can integrate incomplete measurements information from different observation stations and estimate the state parameters of targets. In particular, dual observation stations are adopted for cooperative localization and tracking. The fact is that dual observation stations can offer enough measurements for localization and tracking. In this paper, the target localization and tracking system is completed by dual observation stations with incomplete measurements, and the problem of target localization is done by approach of geometric analysis and the system is observable in this situation. In [11–13], the problem of trajectory parameters estimation is solved with dual or multiple observation stations. All of the methods mentioned in such research are for passive tracking systems with bearings-only measurements. In [14], the authors propose a case where not only the bearings-only measurements but also some other measurements including Doppler shift and distance between the target and each of the sensors are available. But, as described in [14], at least 4 stationary sensors with Doppler-shift measurements can achieve a finite number of solutions for location of target in two-dimensional plane. The number of stationary sensors is usually more than 2 even if there are other types of measurements in addition to the Doppler-shift measurements. What is more, [15] points out that the measurement equation is nonlinear when measurements are Doppler shift, time of arrival, and angle of arrival. In [12], the interactive multiple Kalman filters based on 3D target tracking for two airborne platforms with bearings-only measurements are given. Lots of researches also are concerned with the key point of how to improve localization and tracking accuracy when doing cooperative localization and tracking with dual observation stations. For example, in [16, 17], the interacting multiple models approaches are provided to achieve a superior accuracy in state estimation for manoeuvring target tracking. The analysis of error model and positioning accuracy which helps the layout or placement of observation stations are proposed in [18, 19]. And [15] gives optimality analysis of sensor-target localization geometries. The main idea is that such geometric configurations are likely to result in accurate localization (at least for efficient estimation algorithms). These researches provide different methods to improve localization and tracking accuracy, but all of these researches are based on the situation in which the observation stations are still.

In this paper, the observation stations are assumed to be photoelectric sensors and can provide three types of measurements: azimuth, pitch, and slant range. And the cooperative localization is fulfilled by using the approach of geometric analysis, the azimuth, pitch, and slant range measurements convert to x, y, and z in Cartesian coordinate system and the redundant measurements detected by dual observation stations can be used to eliminate false location point. At last, the measurement equation becomes linear when taking x, y, and z as “new measurements.” What is more, a strategy of cooperative localization and tracking is proposed to handle the problem of imperfect target detection with dual mobile observation stations. For improving tracking precision, a path planning strategy is integrated into the Kalman filtering algorithm. This strategy uses the gradient-descent method which is based on the CEP function of target localization accuracy. Furthermore, the constraints such as communication distance and topography in practical applications are also considered and the corresponding path planning strategy of dual mobile stations is presented.

The main contribution of this paper can be described as follows.

• (1)

  The cooperative localization model of dual stations with incomplete target measurements is established using geometric constraints and the corresponding target localization accuracy is derived by error propagation.

• (2)

  The path planning strategy is integrated into Kalman filtering algorithm to improve target tracking accuracy. The path planning strategy of dual mobile stations is fulfilled in terms of the gradient-descent method which is based on the circular error probability (CEP) function of target localization accuracy and the obstacle avoidance of dual mobile stations is achieved during the target tracking process.

The work considered in this paper can be applied to target localization and tracking in a distributed sensor network with little modification. Distributed sensor networks (DSNs) can use lots of physical sensors (e.g., acoustic, seismic, and visual) to achieve high-quality measurement information. All of the sensors can communicate via the ad hoc wireless network. The localization and tracking system is a part of distributed sensor networks. The data detected by the sensor nodes are processed by processing unit and then a path planning strategy can also be designed for intelligent sensors. What is more, the localization and tracking system has some merits of distributed sensor networks. For example, it has the diversity and redundancy. In [20], diversity means that different sensors can provide complementary perspectives while redundancy denotes the reliability and high-resolution of measurements. Furthermore, the problem of asynchronous time between dual sensors is also solved in this paper and this problem is existing in a distributed sensor network. So the work is applicable in a distributed sensor network.

This paper is organized as follows. Section 2 presents two cooperative localization models as well as the corresponding localization accuracy analysis with incomplete measurements. The target tracking algorithms in several typical different situations with incomplete target information and the path planning strategies of dual mobile stations are described in Section 3, respectively. Section 4 proposes several numerical simulation results while the conclusions of this study are drawn in Section 5.

2. Cooperative Localization Models of Dual Observation Stations

In the process of localization, some measurement information (slant range, azimuth, and pitch angle) of targets may be lost or inaccurate. The cooperative localization of dual or multiple stations can attenuate or even avoid the influence of this imperfect detection. In this paper, the target location is the intersection of localization surfaces and/or lines which is determined by spatial analytic geometry principle. It is often only needed to measure three parameters’ information to obtain the target location in practical applications. For dual or multiple observation stations, the redundant measurement information can be used to improve the localization accuracy. In this section, two different cooperative localization models are proposed and the corresponding localization accuracy parameters are also derived based on the theory of error propagation.

2.1. Case with Missing Slant Range Information

In this case, the slant range information of dual observation stations is missing, and we only use the azimuth and pitch angle information measured by the dual stations in different observation locations to determine the location of target in Cartesian coordinates. As shown in Figure 1, two observation stations are located in A $(x_a,y_a,z_a)$ and C $(x_c,y_c,z_c)$ , respectively, to observe the target T $(x_t,y_t,z_t)$ . Two groups of measurements, the azimuth and pitch angle information of $({\beta }_{\mathrm{1}},{\epsilon }_{\mathrm{1}};{\beta }_{\mathrm{2}},{\epsilon }_{\mathrm{2}})$ can be obtained, which are used to compute the location of target by cooperative localization methods. The measurements $({\beta }_{\mathrm{1}},{\epsilon }_{\mathrm{1}};{\beta }_{\mathrm{2}},{\epsilon }_{\mathrm{2}})$ are detected by the observation stations A and C. These measurements use coordinate systems which take stations as the origin and x, y, and z are parallel to that of the reference coordinate as shown in Figure 1.

Figure 1

Locations of observation stations and target.

For the uniqueness of the solution of target location with the measurements of $({\beta }_{\mathrm{1}},{\epsilon }_{\mathrm{1}};{\beta }_{\mathrm{2}},{\epsilon }_{\mathrm{2}})$ , we have the following assumption.

Assumption 1.

Observation stations A, C and target T are not collinear, and A, C do not overlap.

Two localization surfaces as shown in Figure 2 can be achieved by the azimuth and pitch angle measurements obtained by the station A, and then the corresponding localization line of $AT$ is generated by the above two localization surfaces. In a similar way, the localization line of $CT$ can be determined by the measurements of the station C. The equations of $AT$ and $CT$ are described as (1) and (2), respectively. Then, the intersection of lines of $AT$ and $CT$ , that is, the target T location coordinate $(x_t,y_t,z_t)$ , can be calculated by spatial analytic geometry approach with the form of (3). In this paper, target T and stations A and C are in the same reference Cartesian coordinates system

$$\begin{array}{cc}& \frac{x_t-x_a}{m_{\mathrm{1}}}=\frac{y_t-y_a}{n_{\mathrm{1}}}=\frac{z_t-z_a}{l_{\mathrm{1}}}\end{array}$$

(1)

$$\begin{array}{cc}& \frac{x_t-x_c}{m_{\mathrm{2}}}=\frac{y_t-y_c}{n_{\mathrm{2}}}=\frac{z_t-z_c}{l_{\mathrm{2}}}\end{array}$$

(2)

$$\begin{array}{cc}& x_t=\frac{m_{\mathrm{1}}(n_{\mathrm{2}}(x_c-x_a)-m_{\mathrm{2}}(y_c-y_a))}{m_{\mathrm{1}}{n}_{\mathrm{2}}-m_{\mathrm{2}}{n}_{\mathrm{1}}}+x_a\\ y_t=\frac{n_{\mathrm{1}}(n_{\mathrm{2}}(x_c-x_a)-m_{\mathrm{2}}(y_c-y_a))}{m_{\mathrm{1}}{n}_{\mathrm{2}}-m_{\mathrm{2}}{n}_{\mathrm{1}}}+y_a\\ z_t=\frac{l_{\mathrm{1}}(n_{\mathrm{2}}(x_c-x_a)-m_{\mathrm{2}}(y_c-y_a))}{m_{\mathrm{1}}{n}_{\mathrm{2}}-m_{\mathrm{2}}{n}_{\mathrm{1}}}+z_a,\end{array}$$

(3)

where $m_{\mathrm{1}}=\cos {\epsilon }_{\mathrm{1}}\cos {\beta }_{\mathrm{1}}$ , $n_{\mathrm{1}}=\cos {\epsilon }_{\mathrm{1}}\sin {\beta }_{\mathrm{1}}$ , $l_{\mathrm{1}}=\sin {\epsilon }_{\mathrm{1}}$ , $m_{\mathrm{2}}=\cos {\epsilon }_{\mathrm{2}}\cos {\beta }_{\mathrm{2}}$ , $n_{\mathrm{2}}=\cos {\epsilon }_{\mathrm{2}}\sin {\beta }_{\mathrm{2}}$ , and $l_{\mathrm{2}}=\sin {\epsilon }_{\mathrm{2}}$ .

Figure 2

Localization surfaces of azimuth and pitch angle.

Remark 2.

In terms of Assumption 1, there exists a unique solution of target T location of $(x_t,y_t,z_t)$ in (3). In fact, the direction cosines of the localization lines of $AT$ and $CT$ are $(m_{\mathrm{1}},n_{\mathrm{1}},l_{\mathrm{1}})$ and $(m_{\mathrm{2}},n_{\mathrm{2}},l_{\mathrm{2}})$ , respectively. Therefore, the following inequality in (4) is satisfied

$$\begin{array}{cc}& \frac{m_{\mathrm{1}}}{m_{\mathrm{2}}}\ne \frac{n_{\mathrm{1}}}{n_{\mathrm{2}}}\ne \frac{l_{\mathrm{1}}}{l_{\mathrm{2}}}.\end{array}$$

(4)

The above inequality can be easily obtained due to the noncollinearity of A, C, and T and the nonoverlapping condition of A, C according to Assumption 1 and the principle of spatial analytic geometry. Therefore, $m_{\mathrm{1}}{n}_{\mathrm{2}}-m_{\mathrm{2}}{n}_{\mathrm{1}}\ne \mathrm{0}$ is satisfied.

The corresponding standard deviations of $x_t$ , $y_t$ , and $z_t$ can be derived from (3) using the theory of error propagation, as shown in (5), (6), and (7), respectively

$$\begin{array}{cc}& {\sigma }_{x_t}=\sqrt{\left(a_{\mathrm{1}}^{\mathrm{2}}+a_{\mathrm{3}}^{\mathrm{2}}\right){\sigma }_{\beta }^{\mathrm{2}}+{\sigma }_{x_o}^{\mathrm{2}}},\end{array}$$

(5)

$$\begin{array}{cc}& {\sigma }_{y_t}=\sqrt{\left(b_{\mathrm{1}}^{\mathrm{2}}+b_{\mathrm{3}}^{\mathrm{2}}\right){\sigma }_{\beta }^{\mathrm{2}}+{\sigma }_{y_o}^{\mathrm{2}}},\end{array}$$

(6)

$$\begin{array}{cc}& {\sigma }_{z_t}=\sqrt{\left(c_{\mathrm{1}}^{\mathrm{2}}+c_{\mathrm{3}}^{\mathrm{2}}\right){\sigma }_{\beta }^{\mathrm{2}}+(c_{\mathrm{2}}^{\mathrm{2}}+c_{\mathrm{4}}^{\mathrm{2}}){\sigma }_{\epsilon }^{\mathrm{2}}+{\sigma }_{z_o}^{\mathrm{2}}},\end{array}$$

(7)

where ${\sigma }_{\beta }$ and ${\sigma }_{\epsilon }$ are the standard deviations of azimuth and pitch angle. The standard deviations of ${\sigma }_{x_o}$ , ${\sigma }_{y_o}$ , and ${\sigma }_{z_o}$ are the measurement accuracy of the location of observation stations. The coefficients in (5), (6), and (7) are listed as follows:

$$\begin{array}{l}{a}_{\mathrm{1}}=\frac{m_{\mathrm{2}}[n_{\mathrm{2}}(x_c-x_a)-m_{\mathrm{2}}(y_c-y_a)]{\cos }^{\mathrm{2}}{\epsilon }_{\mathrm{1}}}{{\left(m_{\mathrm{1}}{n}_{\mathrm{2}}-m_{\mathrm{2}}{n}_{\mathrm{1}}\right)}^{\mathrm{2}}},\\ a_{\mathrm{3}}=\frac{m_{\mathrm{1}}[m_{\mathrm{1}}(y_c-y_a)-n_{\mathrm{1}}(x_c-x_a)]{\cos }^{\mathrm{2}}{\epsilon }_{\mathrm{2}}}{{\left(m_{\mathrm{1}}{n}_{\mathrm{2}}-m_{\mathrm{2}}{n}_{\mathrm{1}}\right)}^{\mathrm{2}}},\\ b_{\mathrm{1}}=\frac{n_{\mathrm{2}}[n_{\mathrm{2}}(x_c-x_a)-m_{\mathrm{2}}(y_c-y_a)]{\cos }^{\mathrm{2}}{\epsilon }_{\mathrm{1}}}{{\left(m_{\mathrm{1}}{n}_{\mathrm{2}}-m_{\mathrm{2}}{n}_{\mathrm{1}}\right)}^{\mathrm{2}}},\\ b_{\mathrm{3}}=\frac{n_{\mathrm{1}}[m_{\mathrm{1}}(y_c-y_a)-n_{\mathrm{1}}(x_c-x_a)]{\cos }^{\mathrm{2}}{\epsilon }_{\mathrm{2}}}{{\left(m_{\mathrm{1}}{n}_{\mathrm{2}}-m_{\mathrm{2}}{n}_{\mathrm{1}}\right)}^{\mathrm{2}}},\\ c_{\mathrm{1}}=\frac{l_{\mathrm{1}}[n_{\mathrm{2}}(x_c-x_a)-m_{\mathrm{2}}(y_c-y_a)](n_{\mathrm{1}}{n}_{\mathrm{2}}+m_{\mathrm{1}}{m}_{\mathrm{2}})}{{\left(m_{\mathrm{1}}{n}_{\mathrm{2}}-m_{\mathrm{2}}{n}_{\mathrm{1}}\right)}^{\mathrm{2}}},\\ c_{\mathrm{2}}=\left([\cos {\epsilon }_{\mathrm{1}}(m_{\mathrm{1}}{n}_{\mathrm{2}}-m_{\mathrm{2}}{n}_{\mathrm{1}})+\sin {\epsilon }_{\mathrm{1}}(w_{\mathrm{1}}{n}_{\mathrm{2}}-w_{\mathrm{2}}{m}_{\mathrm{2}})]\right.\\ \hspace{1em}\phantom{\text{iiid}}\times \left.[n_{\mathrm{2}}(x_c-x_a)-m_{\mathrm{2}}(y_c-y_a)]\right)\\ \hspace{1em}\hspace{1em}\times {\left({\left(m_{\mathrm{1}}{n}_{\mathrm{2}}-m_{\mathrm{2}}{n}_{\mathrm{1}}\right)}^{\mathrm{2}}\right)}^{-1},\\ \hspace{1em}{c}_{\mathrm{3}}=\frac{l_{\mathrm{1}}[m_{\mathrm{1}}(y_c-y_a)-n_{\mathrm{1}}(x_c-x_a)](n_{\mathrm{2}}^{\mathrm{2}}+m_{\mathrm{2}}^{\mathrm{2}})}{{\left(m_{\mathrm{1}}{n}_{\mathrm{2}}-m_{\mathrm{2}}{n}_{\mathrm{1}}\right)}^{\mathrm{2}}},\\ \hspace{1em}{c}_{\mathrm{4}}=\frac{l_{\mathrm{1}}\left[n_{\mathrm{1}}(x_c-x_a)m_{\mathrm{1}}(y_c-y_a)\right](m_{\mathrm{2}}{v}_{\mathrm{2}}-n_{\mathrm{2}}{v}_{\mathrm{1}})}{{\left(m_{\mathrm{1}}{n}_{\mathrm{2}}-m_{\mathrm{2}}{n}_{\mathrm{1}}\right)}^{\mathrm{2}}},\end{array}$$

(8)

where $w_{\mathrm{1}}=\sin {\epsilon }_{\mathrm{1}}\cos {\beta }_{\mathrm{1}}$ , $w_{\mathrm{2}}=\sin {\epsilon }_{\mathrm{1}}\sin {\beta }_{\mathrm{1}}$ , $v_{\mathrm{1}}=\sin {\epsilon }_{\mathrm{2}}\cos {\beta }_{\mathrm{2}}$ , and $v_{\mathrm{2}}=\sin {\epsilon }_{\mathrm{2}}\sin {\beta }_{\mathrm{2}}$ .

2.2. Case with Missing Slant Range and Pitch Angle Information

In this case, the cooperative localization information is the azimuth and pitch angle measured by the station A and the slant range and azimuth measured by the station C. As shown in Figure 3, azimuth ${\beta }_{\mathrm{1}}$ and pitch angle ${\epsilon }_{\mathrm{1}}$ of the target T are obtained by the station A, and slant range $d_{\mathrm{2}}$ and azimuth ${\beta }_{\mathrm{2}}$ of the target T are measured by the station C. The measurements of $({\beta }_{\mathrm{1}},{\epsilon }_{\mathrm{1}};d_{\mathrm{2}},{\beta }_{\mathrm{2}})$ can be used to determine the location of the target T by cooperative localization methods. The measurements such as $({\beta }_{\mathrm{1}},{\epsilon }_{\mathrm{1}};d_{\mathrm{2}},{\beta }_{\mathrm{2}})$ are detected by observation stations with corresponding sensors. And these measurements also use coordinate systems which take stations as the origin and x, y, and z are parallel to that of the reference coordinate, respectively, as shown in Figure 1.

Figure 3

Locations of observation stations and target.

Similarly, for the uniqueness of the solution of target location with the measurements of $({\beta }_{\mathrm{1}},{\epsilon }_{\mathrm{1}};d_{\mathrm{2}},{\beta }_{\mathrm{2}})$ , the following assumption is proposed.

Assumption 3.

The projection points $A^{\prime }$ , $C^{\prime }$ , and $T^{\prime }$ of observation stations A, C and target T in the plane of $xoy$ are not collinear, and $A^{\prime }$ , $C^{\prime }$ do not overlap.

As shown in Figure 4, the localization line of $T{T}^{\prime }$ can be obtained by the two localization surfaces of azimuths of ${\beta }_{\mathrm{1}}$ , ${\beta }_{\mathrm{2}}$ . First, we calculate the line of $T{T}^{\prime }$ ; then we determine the intersection of the line of $T{T}^{\prime }$ and the localization spherical surface of $d_{\mathrm{2}}$ . The target location coordinates of $(x_t,y_t,z_t)$ can be calculated and given by

$$\begin{array}{l}{x}_t\\ =\frac{x_c\sin {\beta }_{\mathrm{2}}\cos {\beta }_{\mathrm{1}}-x_a\sin {\beta }_{\mathrm{1}}\cos {\beta }_{\mathrm{2}}+\cos {\beta }_{\mathrm{2}}\cos {\beta }_{\mathrm{1}}(y_a-y_c)}{\sin ({\beta }_{\mathrm{2}}-{\beta }_{\mathrm{1}})}\\ y_t\\ =\frac{y_a\sin {\beta }_{\mathrm{2}}\cos {\beta }_{\mathrm{1}}-y_c\sin {\beta }_{\mathrm{1}}\cos {\beta }_{\mathrm{2}}+\sin {\beta }_{\mathrm{2}}\sin {\beta }_{\mathrm{1}}(x_c-x_a)}{\sin ({\beta }_{\mathrm{2}}-{\beta }_{\mathrm{1}})}\\ \hspace{1em}\hspace{1em}\hspace{1em}{z}_t=\pm \sqrt{d_{\mathrm{2}}^{\mathrm{2}}-{\left(x_t-x_c\right)}^{\mathrm{2}}-{\left(y_t-y_c\right)}^{\mathrm{2}}}+z_c.\end{array}$$

(9)

Figure 4

Localization surfaces of azimuths and slant range.

Remark 4.

In terms of Assumption 3, there exists a unique solution of target T location of $(x_t,y_t,z_t)$ in this case. The intersection of localization line of $T{T}^{\prime }$ and the localization spherical surface of slant range of $d_{\mathrm{2}}$ measured by the station C can be determined using spatial analytic geometry method. If the intersection is two points, we can use the pitch angle information of ${\epsilon }_{\mathrm{1}}$ measured by station A to exclude the false point. If the intersection is just only one point, then the point is the location of target T. Moreover, according to Assumption 3, the measurements of ${\beta }_{\mathrm{1}}$ , ${\beta }_{\mathrm{2}}$ satisfy the following relationship:

$$\begin{array}{cc}& {\beta }_{\mathrm{1}}\ne {\beta }_{\mathrm{2}}\pm k\pi ,\hspace{1em}k=\mathrm{0,1},\mathrm{2},\dots \end{array}$$

(10)

which yields $\sin ({\beta }_{\mathrm{2}}-{\beta }_{\mathrm{1}})\ne \mathrm{0}$ .

Using the theory of error propagation, we can determine the corresponding localization standard deviations of $x_t$ , $y_t$ , and $z_t$ from (9), which are described by (11), (12), and (13), respectively. Consider

$$\begin{array}{cc}& {\sigma }_{x_t}=\sqrt{(a_{\mathrm{2}}^{\mathrm{2}}+a_{\mathrm{4}}^{\mathrm{2}}){\sigma }_{\beta }^{\mathrm{2}}+{\sigma }_{x_o}^{\mathrm{2}}},\end{array}$$

(11)

$$\begin{array}{cc}& {\sigma }_{y_t}=\sqrt{(b_{\mathrm{2}}^{\mathrm{2}}+b_{\mathrm{4}}^{\mathrm{2}}){\sigma }_{\beta }^{\mathrm{2}}+{\sigma }_{y_o}^{\mathrm{2}}},\end{array}$$

(12)

$$\begin{array}{cc}& {\sigma }_{z_t}=\sqrt{c_{\mathrm{6}}^{\mathrm{2}}{\sigma }_d^{\mathrm{2}}+(c_{\mathrm{5}}^{\mathrm{2}}+c_{\mathrm{7}}^{\mathrm{2}}){\sigma }_{\beta }^{\mathrm{2}}+{\sigma }_{z_o}^{\mathrm{2}}}.\end{array}$$

(13)

The coefficients in (11), (12), and (13) are listed as follows:

$$\begin{array}{cc}& a_{\mathrm{2}}=\frac{(x_c-x_a)\sin {\beta }_{\mathrm{2}}\cos {\beta }_{\mathrm{2}}+{\cos }^{\mathrm{2}}{\beta }_{\mathrm{2}}(y_a-y_c)}{\text{si}{\text{n}}^{\mathrm{2}}({\beta }_{\mathrm{2}}-{\beta }_{\mathrm{1}})},\\ a_{\mathrm{4}}=\frac{(x_a-x_c)\sin {\beta }_{\mathrm{1}}\cos {\beta }_{\mathrm{1}}+{\cos }^{\mathrm{2}}{\beta }_{\mathrm{1}}(y_c-y_a)}{\text{si}{\text{n}}^{\mathrm{2}}({\beta }_{\mathrm{2}}-{\beta }_{\mathrm{1}})},\\ b_{\mathrm{2}}=\frac{(x_c-x_a)\text{si}{\text{n}}^{\mathrm{2}}{\beta }_{\mathrm{2}}+\sin {\beta }_{\mathrm{2}}\cos {\beta }_{\mathrm{2}}(y_a-y_c)}{\text{si}{\text{n}}^{\mathrm{2}}({\beta }_{\mathrm{2}}-{\beta }_{\mathrm{1}})},\\ b_{\mathrm{4}}=\frac{(x_a-x_c)\text{si}{\text{n}}^{\mathrm{2}}{\beta }_{\mathrm{1}}+\sin {\beta }_{\mathrm{1}}\cos {\beta }_{\mathrm{1}}(y_c-y_a)}{\text{si}{\text{n}}^{\mathrm{2}}({\beta }_{\mathrm{2}}-{\beta }_{\mathrm{1}})},\\ c_{\mathrm{5}}=\frac{(x_c-x_t+y_c-y_t)b_{\mathrm{2}}}{\sqrt{d_{\mathrm{2}}^{\mathrm{2}}-{\left(x_t-x_c\right)}^{\mathrm{2}}-{\left(y_t-y_c\right)}^{\mathrm{2}}}},\\ c_{\mathrm{6}}=\frac{d_{\mathrm{2}}}{\sqrt{d_{\mathrm{2}}^{\mathrm{2}}-{\left(x_t-x_c\right)}^{\mathrm{2}}-{\left(y_t-y_c\right)}^{\mathrm{2}}}},\\ c_{\mathrm{7}}=\frac{(x_c-x_t+y_c-y_t)b_{\mathrm{4}}}{\sqrt{d_{\mathrm{2}}^{\mathrm{2}}-{\left(x_t-x_c\right)}^{\mathrm{2}}-{\left(y_t-y_c\right)}^{\mathrm{2}}}}.\end{array}$$

(14)

3. Target Tracking Algorithm of Dual Observation Stations

The location information of targets in Cartesian coordinates can be obtained by the cooperative localization models of dual stations with incomplete measurements in Section 2. In this section the state model and the measurement model of targets are presented and the covariance matrixes of measurement noise for each case in Section 2 are derived. The Kalman filtering algorithm is used to track targets.

The incomplete measurement information obtained by stations A and C is collected to the data processing center during the tracking process. Different sampling frequencies of A and C may cause asynchronous time. Therefore, the time alignment of asynchronous measured data should be handled first. The virtual fusion method is used for time synchronization. We assume that the ratio of different sampling periods of stations A, C is integer; then the virtual fusion method is applied to complete the time alignment [21]. In this tracking system, the interpolation method is used for data synchronization. Then, time alignment for sampling time synchronization is done with virtual fusion method. All of these measures can guarantee that the parameters obtained from the two stations are not of phase. The ratio is supposed to be $n=\Delta \tau /\Delta T$ , where $\Delta T$ and $\Delta \tau$ are the sampling periods of stations A and C, respectively. The next updating time of station A then is

$$\begin{array}{cc}& k\Delta \tau =\left[\left(k-\mathrm{1}\right)\Delta \tau +\Delta \tau \right]=\left[\left(k-\mathrm{1}\right)\Delta \tau +n\Delta T\right],\end{array}$$

(15)

where $(k-\mathrm{1})\Delta \tau$ is the latest updating time. Then, the n measurement values of station A are fused into a virtual measurement value as the measurement value of the k time by the least squares method.

After the time alignment for measurement data, the Kalman filtering algorithm is used for target tracking. First, we consider the problem of target tracking with dual still observation stations, and then the corresponding target tracking algorithm is proposed. Then, for improving the tracking accuracy, this paper studies the strategies of target tracking with mobile observation stations. The strategies use the gradient-descent method based on the CEP function of target localization. Here, CEP is an assessment index which is adopted to assess target localization performance and is defined as

$$\begin{array}{cc}& \text{CEP}=\mathrm{0.5887}\left({\sigma }_{x_t}+{\sigma }_{y_t}\right),\end{array}$$

(16)

where ${\sigma }_{x_t}$ and ${\sigma }_{y_t}$ are the standard deviations of target location of $x_t$ and $y_t$ coordinates in the plane of $xoy$ , respectively.

3.1. Still Case

The still case is defined as such situation where both observation stations do not move. The discrete linear state equation and the measurement equation of constant velocity (CV) dynamic system are considered in this paper which can be described as

$$\begin{array}{cc}& {\mathbf{x}}_k=F{\mathbf{x}}_{k-\mathrm{1}}+{\mathbf{v}}_{k-\mathrm{1}},\\ {\mathbf{z}}_k=H{\mathbf{x}}_k+{\mathbf{w}}_k,\end{array}$$

(17)

where F is a deterministic transition matrix, H is a known deterministic matrix, F and H are all nonsingular matrices, and $\mathbf{v}(k)~\mathcal{N}(\mathrm{0},Q)$ is white Gaussian process noise. $\mathbf{w}(k)~\mathcal{N}(\mathrm{0},R_k)$ is a zero-mean independent Gaussian noise. We adopt

$$\begin{array}{cc}& F=\left[\begin{array}{cccccc}\mathrm{1}& \Delta T& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \Delta T& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \Delta T\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right],\hspace{1em}\hspace{1em}{\mathbf{x}}_k=\left[\begin{array}{c}x\left(k\right)\\ \dot{x}\left(k\right)\\ y\left(k\right)\\ \dot{y}\left(k\right)\\ z\left(k\right)\\ \dot{z}\left(k\right)\end{array}\right],\\ H=\left[\begin{array}{cccccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\end{array}\right],\hspace{1em}\hspace{1em}{\mathbf{z}}_k=\left[\begin{array}{c}{x}_t\left(k\right)\\ y_t\left(k\right)\\ z_t\left(k\right)\end{array}\right],\\ E\left({\mathbf{w}}_k{\mathbf{w}}_l^T\right)={\delta }_{kl}Q,\hspace{1em}\hspace{1em}E\left({\mathbf{v}}_k{\mathbf{v}}_l^T\right)={\delta }_{kl}{R}_k,\end{array}$$

(18)

where $\Delta T$ is the sampling period and the Kronecker delta of ${\delta }_{kl}$ is only nonzero for $k=l$ and is used to illustrate that process and measurement noises are white sequences. $R_k$ is furthermore assumed to be positive definite, which implies that all components of the measurement ${\mathbf{z}}_k$ are uncertain.

The measurement noise covariance matrix $R_k$ in the case where the slant range information of both stations A and C is missing can be calculated by

$$\begin{array}{cc}& R_k=\left[\begin{array}{ccc}{R}_{\mathrm{1}}\left(\mathrm{1,1}\right)& R_{\mathrm{1}}\left(\mathrm{1,2}\right)& R_{\mathrm{1}}\left(\mathrm{1,3}\right)\\ R_{\mathrm{1}}\left(\mathrm{2,1}\right)& R_{\mathrm{1}}\left(\mathrm{2,2}\right)& R_{\mathrm{1}}\left(\mathrm{2,3}\right)\\ R_{\mathrm{1}}\left(\mathrm{3,1}\right)& R_{\mathrm{1}}\left(\mathrm{3,2}\right)& R_{\mathrm{1}}\left(\mathrm{3,3}\right)\end{array}\right],\end{array}$$

(19)

where

$$\begin{array}{cc}& R_{\mathrm{1}}\left(\mathrm{1,1}\right)=\left[a_{\mathrm{1}}^{\mathrm{2}}\left(k\right)+a_{\mathrm{3}}^{\mathrm{2}}\left(k\right)\right]{\sigma }_{\beta }^{\mathrm{2}},\\ R_{\mathrm{1}}\left(\mathrm{1,2}\right)=R_{\mathrm{1}}\left(\mathrm{2,1}\right)=\left[a_{\mathrm{1}}\left(k\right)b_{\mathrm{1}}\left(k\right)+a_{\mathrm{3}}\left(k\right)b_{\mathrm{3}}\left(k\right)\right]{\sigma }_{\beta }^{\mathrm{2}},\\ R_{\mathrm{1}}\left(\mathrm{1,3}\right)=R_{\mathrm{1}}\left(\mathrm{3,1}\right)=\left[a_{\mathrm{1}}\left(k\right)c_{\mathrm{1}}\left(k\right)+a_{\mathrm{3}}\left(k\right)c_{\mathrm{3}}\left(k\right)\right]{\sigma }_{\beta }^{\mathrm{2}},\\ R_{\mathrm{1}}\left(\mathrm{2,2}\right)=\left[b_{\mathrm{1}}^{\mathrm{2}}\left(k\right)+b_{\mathrm{3}}^{\mathrm{2}}\left(k\right)\right]{\sigma }_{\beta }^{\mathrm{2}},\\ R_{\mathrm{1}}\left(\mathrm{2,3}\right)=R_{\mathrm{1}}\left(\mathrm{3,2}\right)=\left[b_{\mathrm{1}}\left(k\right)c_{\mathrm{1}}\left(k\right)+b_{\mathrm{3}}\left(k\right)c_{\mathrm{3}}\left(k\right)\right]{\sigma }_{\beta }^{\mathrm{2}},\\ R_{\mathrm{1}}\left(\mathrm{3,3}\right)=\left[c_{\mathrm{1}}^{\mathrm{2}}\left(k\right)+c_{\mathrm{3}}^{\mathrm{2}}\left(k\right)\right]{\sigma }_{\beta }^{\mathrm{2}}+\left[c_{\mathrm{2}}^{\mathrm{2}}\left(k\right)+c_{\mathrm{4}}^{\mathrm{2}}\left(k\right)\right]{\sigma }_{\epsilon }^{\mathrm{2}}.\end{array}$$

(20)

The measurement noise covariance matrix $R_k$ in the case where the slant range information of the station A and the pitch angle information of the station C is missing can be computed by

$$\begin{array}{cc}& R_k=\left[\begin{array}{ccc}{R}_{\mathrm{2}}\left(\mathrm{1,1}\right)& R_{\mathrm{2}}\left(\mathrm{1,2}\right)& R_{\mathrm{2}}\left(\mathrm{1,3}\right)\\ R_{\mathrm{2}}\left(\mathrm{2,1}\right)& R_{\mathrm{2}}\left(\mathrm{2,2}\right)& R_{\mathrm{2}}\left(\mathrm{2,3}\right)\\ R_{\mathrm{2}}\left(\mathrm{3,1}\right)& R_{\mathrm{2}}\left(\mathrm{3,2}\right)& R_{\mathrm{2}}\left(\mathrm{3,3}\right)\end{array}\right],\end{array}$$

(21)

where

$$\begin{array}{cc}& R_{\mathrm{2}}\left(\mathrm{1,1}\right)=\left[a_{\mathrm{2}}^{\mathrm{2}}\left(k\right)+a_{\mathrm{4}}^{\mathrm{2}}\left(k\right)\right]{\sigma }_{\beta }^{\mathrm{2}},\\ R_{\mathrm{2}}\left(\mathrm{1,2}\right)=R_{\mathrm{2}}\left(\mathrm{2,1}\right)=\left[a_{\mathrm{2}}\left(k\right)b_{\mathrm{2}}\left(k\right)+a_{\mathrm{4}}\left(k\right)b_{\mathrm{4}}\left(k\right)\right]{\sigma }_{\beta }^{\mathrm{2}},\\ R_{\mathrm{2}}\left(\mathrm{1,3}\right)=R_{\mathrm{2}}\left(\mathrm{3,1}\right)=\left[a_{\mathrm{2}}\left(k\right)c_{\mathrm{5}}\left(k\right)+a_{\mathrm{4}}\left(k\right)c_{\mathrm{7}}\left(k\right)\right]{\sigma }_{\beta }^{\mathrm{2}},\\ R_{\mathrm{2}}\left(\mathrm{2,2}\right)=\left[b_{\mathrm{2}}^{\mathrm{2}}\left(k\right)+b_{\mathrm{4}}^{\mathrm{2}}\left(k\right)\right]{\sigma }_{\beta }^{\mathrm{2}},\\ R_{\mathrm{2}}\left(\mathrm{2,3}\right)=R_{\mathrm{2}}\left(\mathrm{3,2}\right)=\left[b_{\mathrm{2}}\left(k\right)c_{\mathrm{5}}\left(k\right)+b_{\mathrm{4}}\left(k\right)c_{\mathrm{7}}\left(k\right)\right]{\sigma }_{\beta }^{\mathrm{2}},\\ R_{\mathrm{2}}\left(\mathrm{3,3}\right)=\left[c_{\mathrm{5}}^{\mathrm{2}}\left(k\right)+c_{\mathrm{6}}^{\mathrm{2}}\left(k\right)\right]{\sigma }_{\beta }^{\mathrm{2}}+c_{\mathrm{7}}^{\mathrm{2}}\left(k\right){\sigma }_{\epsilon }^{\mathrm{2}}.\end{array}$$

(22)

By the cooperative localization models proposed in Section 2, the measurement ${\mathbf{z}}_k$ in Cartesian coordinate system can be obtained by (3) and (9), respectively. Then, the Kalman filtering algorithm is used for target tracking.

The detailed steps of the tracking algorithm for this still case are as follows.

Algorithm 5.

Target tracking algorithm for still dual stations is as follows.

• ①

Filter initialization: initialize ${\widehat{\mathbf{x}}}_{\mathrm{0}\mid \mathrm{0}}$ and $P_{\mathrm{0}\mid \mathrm{0}}^{\mathbf{x}\mathbf{x}}$ .

②

State vector prediction:

$$\begin{array}{cc}& {\widehat{\mathbf{x}}}_{k|k-\mathrm{1}}=F{\widehat{\mathbf{x}}}_{k-\mathrm{1}\mid k-\mathrm{1}},\\ P_{k|k-\mathrm{1}}^{\mathbf{x}\mathbf{x}}=F{P}_{k-\mathrm{1}\mid k-\mathrm{1}}^{\mathbf{x}\mathbf{x}}{F}^{\text{T}}+Q.\end{array}$$

(23)

③

Observation-related prediction:

$$\begin{array}{cc}& {\widehat{\mathbf{z}}}_{k|k-\mathrm{1}}=H{\widehat{\mathbf{x}}}_{k\mid k-\mathrm{1}},\\ P_{k|k-\mathrm{1}}^{\mathbf{z}\mathbf{z}}=H{P}_{k\mid k-\mathrm{1}}^{\mathbf{x}\mathbf{x}}{H}^{\text{T}}+R_k,\\ P_{k|k-\mathrm{1}}^{\mathbf{x}\mathbf{z}}=P_{k\mid k-\mathrm{1}}^{\mathbf{z}\mathbf{z}}{H}^{\text{T}}.\end{array}$$

(24)

④

Kalman filter update:

$$\begin{array}{cc}& K_k\triangleq P_{k|k-\mathrm{1}}^{\mathbf{x}\mathbf{z}}{\left(P_{k\mid k-\mathrm{1}}^{\mathbf{z}\mathbf{z}}\right)}^{-\mathrm{1}},\\ {\widehat{\mathbf{x}}}_{k\mid k}={\widehat{\mathbf{x}}}_{k\mid k-\mathrm{1}}+K_k\left({\mathbf{z}}_k-{\widehat{\mathbf{z}}}_{k\mid k-\mathrm{1}}\right),\\ P_{k\mid k}^{\mathbf{x}\mathbf{x}}=P_{k\mid k-\mathrm{1}}^{\mathbf{x}\mathbf{x}}-K_k{P}_{k\mid k-\mathrm{1}}^{\mathbf{z}\mathbf{z}}{K}_k^{\text{T}}.\end{array}$$

(25)

⑤

Store results: store ${\widehat{\mathbf{x}}}_{k\mid k}$ and $P_{k\mid k}^{\mathbf{x}\mathbf{x}}$ to a Track File.

⑥

Return to ②.

3.2. Moving Case

The moving case is defined as such situation where dual observation stations can move autonomously. Changing the placement of observation stations can improve the localization and tracking accuracy when the observation stations are still during the localization and tracking process. But this method cannot respond to the changes of target route in time. Here, we propose a path planning strategy of the observation stations for target tracking on the basis of the gradient of CEP. This approach can significantly improve the tracking accuracy.

In this case, we consider two typical scenarios of moving dual stations. One is the case in which dual stations can move in the plane of $xoy$ without any additional constraints to track targets, and the other is the case in which dual stations’ moving paths are subjected to the constraints of topography (mountain, hollow), obstacle, and communication distance. For each case, the corresponding target tracking algorithms are presented in subsequence sections.

3.2.1. Without Constraints

In the process of target tracking, the predicted value at the current time can be achieved with the last filtering value and the predicted CEP value at the current time can also be obtained. The analytical expression of the predicted CEP about $(x_a,y_a)$ , $(x_c,y_c)$ , and $(x_t,y_t,z_t)$ is derived from the analysis of the localization accuracy and it is described as

$$\begin{array}{cc}& \text{CEP}=\mathrm{0.5887}\left({\sigma }_{x_t}+{\sigma }_{y_t}\right)=f\left(x_a,y_a,x_c,y_c,x_t,y_t,z_t\right).\end{array}$$

(26)

Then, use the predicted value of the target location to calculate the gradient of CEP by

$$\begin{array}{cc}& \text{CEP}{A}_{\text{grad}}=\left[\frac{\partial f}{\partial x_a}\frac{\partial f}{\partial y_a}\right],\\ \text{CEP}{C}_{\text{grad}}=\left[\frac{\partial f}{\partial x_c}\frac{\partial f}{\partial y_c}\right].\end{array}$$

(27)

The detailed formulas of (27) can be found in the Appendix. Then, the direction cosines of the above gradients can be determined by

$$\begin{array}{cc}& \text{CEP}{A}_{dc}=\frac{\text{CEP}{A}_{\text{grad}}}{\left|\text{CEP}{A}_{\text{grad}}\right|},\\ \text{CEP}{C}_{dc}=\frac{\text{CEP}{C}_{\text{grad}}}{\left|\text{CEP}{C}_{\text{grad}}\right|}.\end{array}$$

(28)

Station A and station C can move along the opposite direction of CEP gradient from $(x_a(k-\mathrm{1}),y_a(k-\mathrm{1}))$ or $(x_c(k-\mathrm{1}),y_c(k-\mathrm{1}))$ to $(x_a(k),y_a(k))$ or $(x_c(k),y_c(k))$ , respectively, according to (29) to improve tracking accuracy. Consider

$$\begin{array}{cc}& x_a\left(k\right)=x_a\left(k-\mathrm{1}\right)-\text{CEP}{A}_{dc}\left(\mathrm{1}\right)\times v_a\times \Delta T,\\ y_a\left(k\right)=y_a\left(k-\mathrm{1}\right)-\text{CEP}{A}_{dc}\left(\mathrm{2}\right)\times v_a\times \Delta T,\\ x_c\left(k\right)=x_c\left(k-\mathrm{1}\right)-\text{CEP}{C}_{dc}\left(\mathrm{1}\right)\times v_c\times \Delta T,\\ y_c\left(k\right)=y_c\left(k-\mathrm{1}\right)-\text{CEP}{C}_{dc}\left(\mathrm{2}\right)\times v_c\times \Delta T,\end{array}$$

(29)

where $v_a$ and $v_c$ are the velocities of stations A and C, respectively. Here, we present an algorithm of target tracking by integrating path planning strategy of observation stations into Kalman filtering algorithm. The detailed steps of the proposed target tracking algorithm are listed as follows.

Algorithm 6.

Target tracking algorithm for moving dual stations without constraints is as follows.

• ①

Filter initialization: initialize ${\widehat{\mathbf{x}}}_{\mathrm{0}\mid \mathrm{0}}$ and $P_{\mathrm{0}\mid \mathrm{0}}^{\mathbf{x}\mathbf{x}}$ .

②

State vector prediction:

$$\begin{array}{cc}& {\widehat{\mathbf{x}}}_{k\mid k-\mathrm{1}}=F{\widehat{\mathbf{x}}}_{k-\mathrm{1}\mid k-\mathrm{1}},\\ P_{k\mid k-\mathrm{1}}^{\mathbf{x}\mathbf{x}}=F{P}_{k-\mathrm{1}\mid k-\mathrm{1}}^{\mathbf{x}\mathbf{x}}{F}^{\text{T}}+Q.\end{array}$$

(30)

③

Observation-related prediction:

$$\begin{array}{cc}& {\widehat{\mathbf{z}}}_{k\mid k-\mathrm{1}}=H{\widehat{\mathbf{x}}}_{k\mid k-\mathrm{1}},\\ P_{k\mid k-\mathrm{1}}^{\mathbf{z}\mathbf{z}}=H{P}_{k\mid k-\mathrm{1}}^{\mathbf{x}\mathbf{x}}{H}^{\text{T}}+R_k,\\ P_{k\mid k-\mathrm{1}}^{\mathbf{x}\mathbf{z}}=P_{k\mid k-\mathrm{1}}^{\mathbf{z}\mathbf{z}}{H}^{\text{T}}.\end{array}$$

(31)

④

Path planning strategies are as follows.

• (1)

Substitute ${\widehat{\mathbf{x}}}_{k\mid k-\mathrm{1}}(\mathrm{1})$ , ${\widehat{\mathbf{x}}}_{k\mid k-\mathrm{1}}(\mathrm{3})$ , ${\widehat{\mathbf{x}}}_{k\mid k-\mathrm{1}}(\mathrm{5})$ and $(x_a(k-\mathrm{1}),y_a(k-\mathrm{1}))$ , $(x_c(k-\mathrm{1}),y_c(k-\mathrm{1}))$ into (27) to calculate $\text{CEP}{C}_{\text{grad}}$ .

(2)

Station C moves to $(x_c(k),y_c(k))$ by (29).

(3)

Substitute the ${\widehat{\mathbf{x}}}_{k\mid k-\mathrm{1}}(\mathrm{1})$ , ${\widehat{\mathbf{x}}}_{k\mid k-\mathrm{1}}(\mathrm{3})$ , ${\widehat{\mathbf{x}}}_{k\mid k-\mathrm{1}}(\mathrm{5})$ and $(x_a(k-\mathrm{1}),y_a(k-\mathrm{1}))$ , $(x_c(k),y_c(k))$ into (27) to calculate $\text{CEP}{A}_{\text{grad}}$ .

(4)

Station A moves to $(x_a(k),y_a(k))$ by (29).

⑤

Kalman filter update:

$$\begin{array}{cc}& K_k\triangleq P_{k\mid k-\mathrm{1}}^{\mathbf{x}\mathbf{z}}{\left(P_{k\mid k-\mathrm{1}}^{\mathbf{z}\mathbf{z}}\right)}^{-\mathrm{1}},\\ {\widehat{\mathbf{x}}}_{k\mid k}={\widehat{\mathbf{x}}}_{k\mid k-\mathrm{1}}+K_k\left({\mathbf{z}}_k-{\widehat{\mathbf{z}}}_{k\mid k-\mathrm{1}}\right),\\ P_{k\mid k}^{\mathbf{x}\mathbf{x}}=P_{k\mid k-\mathrm{1}}^{\mathbf{x}\mathbf{x}}-K_k{P}_{k\mid k-\mathrm{1}}^{\mathbf{z}\mathbf{z}}{K}_k^{\text{T}}.\end{array}$$

(32)

⑥

Store results: store ${\widehat{\mathbf{x}}}_{k|k}$ and $P_{k|k}^{\mathbf{x}\mathbf{x}}$ to a Track File.

⑦

Return to ②.

3.2.2. With Constraints

In practical engineering applications, the constraints such as topography (mountain, hollow), obstacle, and communication distance during observation stations moving process need to be taken into account. The distance between stations A and C can be calculated by

$$\begin{array}{cc}& l=\sqrt{{\left(x_a-x_c\right)}^{\mathrm{2}}+{\left(y_a-y_c\right)}^{\mathrm{2}}+{\left(z_a-z_c\right)}^{\mathrm{2}}}\end{array}$$

(33)

and then the maximum admissible communication distance is defined as $l_c$ .

In the process in which the stations move along the optimal paths determined by the gradient of CEP, the data processing center needs to judge whether the stations will encounter some obstacles or the distances among the stations are out of the maximum admissible communication distance. When the station C meets obstacles between $(x_c(k-\mathrm{1}),y_c(k-\mathrm{1}))$ and $(x_c(k),y_c(k))$ or the station A meets obstacles between $(x_a(k-\mathrm{1}),y_a(k-\mathrm{1}))$ and $(x_a(k),y_a(k))$ , the desired direction angles of ${\varphi }_d^A$ and ${\varphi }_d^C$ [22] need to be found by (34). In Figure 5, $d_r$ and $d_l$ are tangents. The angles of ${\beta }_l$ and ${\beta }_r$ can be obtained by direction detection sensors. Then, such direction cosines are given in (35). Consider

$$\begin{array}{cc}& {\varphi }_d^A={\tan }^{-\mathrm{1}}\left[\frac{{\text{CEP}A}_{dc}\left(\mathrm{2}\right)}{{\text{CEP}A}_{dc}\left(\mathrm{1}\right)}\right]+{\beta }_l^A,\hspace{1em}{\beta }_l^A\le {\beta }_r^A,\\ {\varphi }_d^A={\tan }^{-\mathrm{1}}\left[\frac{{\text{CEP}A}_{dc}\left(\mathrm{2}\right)}{{\text{CEP}A}_{dc}\left(\mathrm{1}\right)}\right]-{\beta }_r^A,\hspace{1em}{\beta }_l^A>{\beta }_r^A,\\ {\varphi }_d^C={\tan }^{-\mathrm{1}}\left[\frac{{\text{CEP}C}_{dc}\left(\mathrm{2}\right)}{{\text{CEP}C}_{dc}\left(\mathrm{1}\right)}\right]+{\beta }_l^C,\hspace{1em}{\beta }_l^C\le {\beta }_r^C,\\ {\varphi }_d^A={\tan }^{-\mathrm{1}}\left[\frac{{\text{CEP}C}_{dc}\left(\mathrm{2}\right)}{{\text{CEP}C}_{dc}\left(\mathrm{1}\right)}\right]-{\beta }_r^C,\hspace{1em}{\beta }_l^C>{\beta }_r^C,\end{array}$$

(34)

$$\begin{array}{cc}& C_{dc}=\left[\begin{array}{cc}\cos {\varphi }_d^C& \sin {\varphi }_d^C\end{array}\right],\\ A_{dc}=\left[\begin{array}{cc}\cos {\varphi }_d^A& \sin {\varphi }_d^A\end{array}\right].\end{array}$$

(35)

Then, station A and station C can move from $(x_a(k-\mathrm{1}),y_a(k-\mathrm{1}))$ or $(x_c(k-\mathrm{1}),y_c(k-\mathrm{1}))$ to $(x_a(k),y_a(k))$ or $(x_c(k),y_c(k))$ by (36) to avoid the obstacle. Consider

$$\begin{array}{cc}& x_a\left(k\right)=x_a\left(k-\mathrm{1}\right)+A_{dc}\left(\mathrm{1}\right)\times v_a\times \Delta T,\\ y_a\left(k\right)=y_a\left(k-\mathrm{1}\right)+A_{dc}\left(\mathrm{2}\right)\times v_a\times \Delta T,\\ x_c\left(k\right)=x_c\left(k-\mathrm{1}\right)+C_{dc}\left(\mathrm{1}\right)\times v_c\times \Delta T,\\ y_c\left(k\right)=y_c\left(k-\mathrm{1}\right)+C_{dc}\left(\mathrm{2}\right)\times v_c\times \Delta T.\end{array}$$

(36)

Here, an algorithm of target tracking combining the path planning strategies of moving observation stations A and C with constraints is proposed. The detailed steps of such tracking algorithm are as follows.

Figure 5

Desired direction.

Algorithm 7.

Target tracking algorithm for moving dual stations with constraints is as follows.

• ①

Filter initialization: initialize ${\widehat{\mathbf{x}}}_{\mathrm{0}\mid \mathrm{0}}$ and $P_{\mathrm{0}\mid \mathrm{0}}^{\mathbf{x}\mathbf{x}}$ .

②

State vector prediction:

$$\begin{array}{cc}& {\widehat{\mathbf{x}}}_{k\mid k-\mathrm{1}}=F{\widehat{\mathbf{x}}}_{k-\mathrm{1}\mid k-\mathrm{1}},\\ P_{k\mid k-\mathrm{1}}^{\mathbf{x}\mathbf{x}}=F{P}_{k-\mathrm{1}\mid k-\mathrm{1}}^{\mathbf{x}\mathbf{x}}{F}^{\text{T}}+Q.\end{array}$$

(37)

③

Observation-related prediction:

$$\begin{array}{cc}& {\widehat{\mathbf{z}}}_{k\mid k-\mathrm{1}}=H{\widehat{\mathbf{x}}}_{k\mid k-\mathrm{1}},\\ P_{k\mid k-\mathrm{1}}^{\mathbf{z}\mathbf{z}}=H{P}_{k\mid k-\mathrm{1}}^{\mathbf{x}\mathbf{x}}{H}^{\text{T}}+R_k,\\ P_{k\mid k-\mathrm{1}}^{\mathbf{x}\mathbf{z}}=P_{k\mid k-\mathrm{1}}^{\mathbf{z}\mathbf{z}}{H}^{\text{T}}.\end{array}$$

(38)

④

Path planning strategies are as follows.

• (1)

Substitute ${\widehat{\mathbf{x}}}_{k\mid k-\mathrm{1}}(\mathrm{1})$ , ${\widehat{\mathbf{x}}}_{k\mid k-\mathrm{1}}(\mathrm{3})$ , ${\widehat{\mathbf{x}}}_{k\mid k-\mathrm{1}}(\mathrm{5})$ and $(x_a(k-\mathrm{1}),y_a(k-\mathrm{1}))$ , $(x_c(k-\mathrm{1}),y_c(k-\mathrm{1}))$ into (27) to calculate $\text{CEP}{C}_{\text{grad}}$ .

(2)

Calculate $(x_c(k),y_c(k))$ by (29).

(3)

Substitute $(x_a(k-\mathrm{1}),y_a(k-\mathrm{1}))$ , $(x_c(k),y_c(k))$ into (33) to calculate l.

(4)

Station C stops moving if $l>l_c$ and go to step $(7)$ ; otherwise go to step $(5)$ .

(5)

Station C moves to $(x_c(k),y_c(k))$ by (29) if there is no obstacle and go to step $(7)$ ; otherwise go to step $(6)$ .

(6)

Station C moves to $(x_c(k),y_c(k))$ by (35).

(7)

Substitute the ${\widehat{\mathbf{x}}}_{k\mid k-\mathrm{1}}(\mathrm{1})$ , ${\widehat{\mathbf{x}}}_{k\mid k-\mathrm{1}}(\mathrm{3})$ , ${\widehat{\mathbf{x}}}_{k\mid k-\mathrm{1}}(\mathrm{5})$ and $(x_a(k-\mathrm{1}),y_a(k-\mathrm{1}))$ , $(x_c(k),y_c(k))$ into (27) to calculate $\text{CEP}{A}_{\text{grad}}$ .

(8)

Calculate $(x_a(k),y_a(k))$ by (29).

(9)

Substitute $(x_a(k),y_a(k))$ , $(x_c(k),y_c(k))$ into (33) to calculate l.

(10)

Station A stops moving if $l>l_c$ and go to step ⑤; otherwise go to step $(11)$ .

(11)

Station A moves to $(x_a(k),y_a(k))$ by (29) if there is no obstacle and go to step ⑤; otherwise go to step $(12)$ .

(12)

Station A moves to $(x_a(k),y_a(k))$ by (35).

⑤

Kalman filter update:

$$\begin{array}{cc}& K_k\triangleq P_{k\mid k-\mathrm{1}}^{\mathbf{x}\mathbf{z}}{\left(P_{k\mid k-\mathrm{1}}^{\mathbf{z}\mathbf{z}}\right)}^{-\mathrm{1}},\\ {\widehat{\mathbf{x}}}_{k\mid k}={\widehat{\mathbf{x}}}_{k\mid k-\mathrm{1}}+K_k\left({\mathbf{z}}_k-{\widehat{\mathbf{z}}}_{k\mid k-\mathrm{1}}\right),\\ P_{k\mid k}^{\mathbf{x}\mathbf{x}}=P_{k\mid k-\mathrm{1}}^{\mathbf{x}\mathbf{x}}-K_k{P}_{k\mid k-\mathrm{1}}^{\mathbf{z}\mathbf{z}}{K}_k^{\text{T}}.\end{array}$$

(39)

⑥

Store results: store ${\widehat{\mathbf{x}}}_{k\mid k}$ and $P_{k\mid k}^{\mathbf{x}\mathbf{x}}$ to a Track File.

⑦

Return to ②.

4. Simulation and Analysis

The Monte Carlo simulations have been done for the tracking strategies and the RMSE curves of x, y, and z are given, respectively. The RMSE is defined as

$$\begin{array}{cc}& \text{RMSE}\left(k\right)=\sqrt{\frac{\mathrm{1}}{N}{\displaystyle \sum }_{i=\mathrm{1}}^N‍{\left[{\widehat{\mathbf{x}}}_i\left(k\right)-\mathbf{x}\left(k\right)\right]}^{\mathrm{2}}},\end{array}$$

(40)

where N is the number of Monte Carlo simulation times and $\mathbf{x}(k)$ is the true value of target state at the kth sampling time. ${\widehat{\mathbf{x}}}_i(k)$ is the corresponding estimated value of target state at time k in the ith simulation.

For comparison, the corresponding Cramer-Rao lower bound (CRLB) is also given, which can be calculated in a recursive way by [23, 24]

$$\begin{array}{cc}& J_{k+\mathrm{1}}={\left({\Gamma }^T{Q}_k\Gamma +\Phi J_k^{-\mathrm{1}}{\Phi }^T\right)}^{-\mathrm{1}}+H^T{R}_{k+\mathrm{1}}^{-\mathrm{1}}H,\\ C_k\triangleq J_k^{-\mathrm{1}},\end{array}$$

(41)

where $J_{k+\mathrm{1}}$ is the Fisher information matrix and $C_k$ is the corresponding error covariance matrix.

What is more, in this paper, simulation results are also compared with the results that use single observation station with complete measurements: azimuth, pitch angle, and slant range.

4.1. Simulation without Constraints

In the Cartesian coordinate system, the initial coordinates of stations A and C are located at $(\mathrm{0,0},\mathrm{0})$ and $(\mathrm{1000,1000,0})$ , respectively. The target is assumed to fly along a horizontal beeline without any posture change at a height of $z=\mathrm{500}$  m and the velocity is set as $v=\mathrm{80}$  m/s. The sample period of station A is $\Delta T=\mathrm{20}$  ms and the sampling period of station C is $\Delta \tau =\mathrm{40}$  ms. The total sampling time is denoted as t. The standard deviations of slant range, azimuth, and pitch angle of both station A and station C are chosen as ${\sigma }_d=\mathrm{5}$  m, ${\sigma }_{\beta }=\mathrm{2}$  mil, and ${\sigma }_{\epsilon }=\mathrm{2}$  mil, respectively.

The localization accuracy of stations A and C is defined as ${\sigma }_{x_o}={\sigma }_{y_o}={\text{CEP}}_{\text{GPS}}/\sqrt{\mathrm{2}\ln \mathrm{2}}$ , ${\sigma }_{z_o}=\mathrm{20}$  m, and CEP_GPS is the localization accuracy of observation stations using GPS. Here, we choose $\text{CE}{\text{P}}_{\text{GPS}}$ as 15 meters. The velocities of observation stations A and C are both set as $v_a=v_c=\mathrm{10}$  m/s. Finally, we define $T_{\mathrm{0}}$ as the initial location of target and θ as the course angle of target flight. The simulation results are compared not only with Cramer-Rao lower bound (CRLB) but also with results of single station, and the corresponding simulation parameters are set to be the same as that in the situation of dual stations. The single observation station $A^{\prime }$ is also located at $(\mathrm{0,0},\mathrm{0})$ :

• (A)

  measurement information of station A: ${\beta }_{\mathrm{1}}$ , ${\epsilon }_{\mathrm{1}}$ ; measurement information of station C: ${\beta }_{\mathrm{2}}$ , ${\epsilon }_{\mathrm{2}}$ ; measurement information of single station $A^{\prime }$ : ${\beta }^{\prime }$ , ${\epsilon }^{\prime }$ , $d^{\prime }$ ; $\theta =-\mathrm{4}{\mathrm{0}}^{\circ }$ , $T_{\mathrm{0}}=(-\mathrm{3500,1500,500})$ , $t=\mathrm{40}$  s (see Figure 6),

• (B)

  measurement information of station A: ${\beta }_{\mathrm{1}}$ , ${\epsilon }_{\mathrm{1}}$ ; measurement information of station C: $d_{\mathrm{2}}$ , ${\beta }_{\mathrm{2}}$ ; measurement information of single station $A^{\prime }$ : ${\beta }^{\prime }$ , ${\epsilon }^{\prime }$ , $d^{\prime }$ ; $\theta =-\mathrm{4}{\mathrm{0}}^{\circ }$ , $T_{\mathrm{0}}=(-\mathrm{2000,3000,500})$ , $t=\mathrm{40}$  s (see Figure 7).

Figure 6

The simulation results of (A).

Figure 7

The simulation results of (B).

The still case is defined as such situation where both observation stations are still while the moving case describes the scene where dual stations move along the planned paths; the single observation station legend means that the tracking system uses only one observation station. The results of still case are achieved with Algorithm 5 while the results of moving case are achieved with Algorithm 6, and the corresponding simulation results are shown in Figures 6 and 7, respectively. The trajectories of target T and stations A, C are shown in Figures 6(a) and 7(a), respectively. The RMSE curves of x, y, and z, respectively, are given in Figures 6(b), 6(c), and 6(d) and Figures 7(b), 7(c), and 7(d).

The RMSE curves of x, y, and z in the cases of (A) and (B) show that Algorithm 6 using path planning strategy possesses smaller tracking errors than Algorithm 5 where stations are still. Furthermore, the RMSE curves are closer to the ideal CRLB curves, and it is shown that the tracking accuracy has a significant improvement with the proposed path planning strategy. The simulation results are found to be the same as what is expected and generally speaking both the proposed algorithms have better performance than the existing technique using single observation station. The reason is that dual stations have redundant measurements for cooperative localization and tracking. In Algorithm 6, the state vector prediction is used to calculate optimal path planning and a higher localization accuracy will be achieved if stations move along the planned path. Such planning process is conducted repeatedly during the whole tracking process. So tracking accuracy will become higher than that of Algorithm 5 obviously. Figures 6 and 7 provide intuitive verification.

4.2. Simulation with Constraints

The simulation environment is set to be the same as that in Section 4.1 and we use the rectangular and circular to simulate different obstacles (mountain or pothole) that the observation stations may encounter (see Figures 8(a) and 9(a)). The maximum admissible communication distance is set as $l_c=\mathrm{2500}$  m. The observation stations stop moving if the distance between stations A and C is out of the maximum admissible communication distance:

• (A)

  measurement information of station A: ${\beta }_{\mathrm{1}}$ , ${\epsilon }_{\mathrm{1}}$ ; measurement information of station C: ${\beta }_{\mathrm{2}}$ , ${\epsilon }_{\mathrm{2}}$ ; $\theta =-\mathrm{4}{\mathrm{0}}^{\circ }$ , $T_{\mathrm{0}}=(-\mathrm{3500,1500,500})$ , $t=\mathrm{40}$  s (see Figure 8),

• (B)

  measurement information of station A: ${\beta }_{\mathrm{1}}$ , ${\epsilon }_{\mathrm{1}}$ ; measurement information of station C: $d_{\mathrm{2}}$ , ${\beta }_{\mathrm{2}}$ ; $\theta =-\mathrm{4}{\mathrm{0}}^{\circ }$ , $T_{\mathrm{0}}=(-\mathrm{2000,3000,500})$ , $t=\mathrm{40}$  s (see Figure 9).

Figure 8

The simulation results of (A).

Figure 9

The simulation results of (B).

The simulation results in cases of (A) and (B) are shown in Figures 8 and 9. The trajectories of target T and stations A, C are given in Figures 8(a) and 9(a). The RMSE curves of x, y, and z are shown in Figures 8(b), 8(c), and 8(d) and Figures 9(b), 9(c), and 9(d). The results of situation without constraints are achieved with Algorithm 6 while the results of situation with constraints are achieved with Algorithm 7.

Figures 8(a) and 9(a) show that the observation stations can effectively avoid obstacles in the process of movement and the observation stations will stop if the distance between them is out of the maximum admissible communication distance. The RMSE curves in Figures 8(b), 8(c), and 8(d) and Figures 9(b), 9(c), and 9(d) show that the tracking accuracy has some deterioration in the situation with constraints. In fact, Algorithm 6 is the optimal method to achieve high tracking accuracy and there is no doubt that the RMSE curves of Algorithm 6 are closer to the ideal CRLB curves. Algorithm 7 is concerned with actual situations and can offer rational strategy to avoid constraints. The negative impact of constraints is that the tracking accuracy will decline more or less. Generally speaking, Algorithm 7 is more comprehensive and can be used for target localization and tracking in actual distributed sensor networks.

5. Conclusions

This paper studies the cooperative localization of dual observation stations with incomplete measurements and gives path planning strategies of the observation stations for target tracking. The cooperative localization method can locate the target in the situation where the slant range measurements of both observation stations are missing or the situation where the slant range measurement of station A and the pitch angle measurement of station C are missing. The optimal path planning strategies of dual mobile stations for tracking have been proposed. Furthermore, this paper studies the strategies for tracking with the constraints such as communication distance and topography. At last several simulation experiments have been presented and the simulation results show that the target tracking accuracy has a significant improvement when the optimal path planning strategies are adopted. The mobile observation stations can move normally under constraints, but the constraints may lead to a decline in tracking accuracy. The proposed strategy of target tracking can provide technical support for high-precision cooperative target tracking with multiple mobile observation stations under modern electronic countermeasure environment.