Fusion of Angle Measurements from Hull Mounted and Towed Array Sensors

Abstract

Two sensor arrays, hull-mounted array, and towed array sensors are considered for bearings-only tracking. An algorithm is designed to combine the information obtained as bearing (angle) measurements from both sensor arrays to give a better solution. Using data from two different sensor arrays reduces the problem of observability and the observer need not follow the S-maneuver to attain observability of the process. The performance of the fusion algorithm is comparable to that of theoretical Cramer–Rao lower bound and with that of the algorithm when bearing measurements from a single sensor array are considered. Different filters are used for analyzing both algorithms. Monte Carlo runs need to be done to evaluate the performance of algorithms more accurately. Also, the performance of the fusion algorithm is evaluated in terms of solution convergence time.

1. Introduction

The demand for fast and accurate tracking of the targets is increasing day by day. Under water, noise corrupted angle measurements are available from noise radiated from the target and processed from a single observer platform operating in passive mode [1]. Passive target tracking is more popular as it helps the tracker from being tracked. In the ocean environment, target motion analysis is generally performed in a two-dimensional plane using angles-only measurements. The angles generated by the sensors are well known as bearing angles. The observer with the single sensor array, say hull-mounted array (HMA), in a passive mode, determines parameters of target motion like its range, course, bearing angle, speed, etc. As the range of the target is not available and bearing being nonlinearly related to target motion parameters, nonlinearity in the process is high [2]. In addition to this, measurements from a single sensor array are only considered to be available, which makes the process unobservable and the observer must perform a proper maneuver to make the process observable [3,4,5].

The ocean environment is unstable all the time and this makes the bearing measurements from a single sensor array unreliable in obtaining all the required information about the target. Moreover, an increase in system complexity increases the need for more input from different sensor arrays to obtain complete information on the system. According to the application and need, a system with multiple sensors can choose different sensors with the required characteristics and based on the environmental conditions for obtaining the target information. Information from these sensors is then consistently combined coherently to obtain the parameters like range, course, velocity, and acceleration of the target during its movement. This is known as data fusion problem. As the bearing measurements are taken from two sensors, in this work, the maneuver of the observer to obtain the observability of the target is reduced in most of the scenarios.

Angles-only tracking, commonly known as bearings-only tracking, used classical filters like maximum likelihood estimator, pseudo linear estimator, etc. in the early 19th century [1,6] and these were batch processing filters. In real-time scenarios, these batch processing estimators need a lot of time for obtaining the estimated results. This led to the evolution of recursive estimators like least square estimator, Kalman filter (KF), etc. KF being a linear filter was not suitable for nonlinear bearings-only tracking [3]. The nonlinearity in the process leads to the evolution of extended KF (EKF), unscented KF (UKF), particle filter (PF), etc. to deal with the nonlinearities in the process and measurements [7]. EKF reduces the nonlinearities by linearizing the measurement equation using Taylor series expansion and updates the state and covariance matrices and calculates the filter gain. The filter is mostly dependent on the initializations and due to the linearization process; information is lost making the solution unreliable. In UKF, the probability density function of the state vector is propagated through a set of carefully chosen sample points, called sigma points, using unscented transform [8,9]. This helps in possession of posterior mean and covariance up to second order for any type of nonlinearity. Hence, unlike EKF, UKF can provide improved stability, accuracy, and better convergence in solution [10,11,12].

Y. Bar-Shalom, in [13], gave a deep insight regarding data association and data fusion of multiple sensors in multiple target environment using radar measurements. Several filtering algorithms were utilized including EKF, UKF, and PF. Data fusion for underwater surveillance was presented in [14,15,16], in which the technique was evaluated only using MGBEKF and UKF algorithms. In [17], the author used a data fusion technique to evaluate the performance of the fusion technique for two and more than two sensors of data fusion. It is an extensive work for the work presented by Bar-Shalom on data fusion. In [18], the authors applied a data fusion technique using a PF for estimation of robot motions. In [19], the authors analyzed the fusion technique for batch processing using the maximum likelihood estimator and analyzed the same with recursive processing with an unscented Gauss–Helmert Kalman filter. The observability of the system is further discussed in the paper. In [18], a data fusion algorithm was used along with a data association algorithm to detect multiple targets. It deals with the hardware complexity in target detection. In [19], data fusion was used to detect the target in a properly designed underwater network. In [20], authors used interacting multiple model techniques combined with EKF and UKF for tracking the single as well as multiple targets and compared two different data fusion techniques that involve fusion using nearest neighborhood calculation and fusion using probabilistic data association. The data fusion techniques in literature are extended to robotics, in air surveillance, and underwater surveillance applications. However, there is a necessity to explore the fusion techniques in detail for underwater surveillance applications.

The filters considered for evaluation of fusion algorithm are further compared with theoretical Cramer–Rao lower bound (CRLB). The mathematical modelling of CRLB for tracking algorithms using a single sensor array is in [9,21] and that for tracking algorithms using two sensor arrays is represented in Section 2.6.

In this paper, a single target following the constant velocity model is considered for simplicity. The fusion of measurement information from two sensor arrays, HMA and towed array (TA) is considered without maneuver. TA is towed behind the observer at a certain distance away from the observer. This helps in the reduction of self-noise from the observer and allows it to operate at very low frequencies. Hence, the TA will be capable of even tracking the target at great distances. The algorithm is evaluated for several scenarios in Monte Carlo simulations.

In this paper, the performance of BOT with a single sensor (HMA) and two sensors (HMA and TA) using different filters is evaluated. The mathematical modelling of all the algorithms is described in Section 2. The results obtained from the simulation of the algorithms and discussion on the results is given in Section 3. This paper is concluded in Section 4.

2. Mathematical Modeling

2.1. Tracking Using a Single Sensor Array

Tracking using a single sensor array (SST) measurement considers the bearing measurements from HMA only. The observer is assumed to be following constant speed with S—the maneuver in the course and the target follows persistence speed and course [20]. At time instant, say ‘$\tau$’, state vector of the observer [9,10] is given as Equation (1).

$$S_o\left(\tau \right)={\left[\begin{array}{cccc}{v}_{xo}\left(\tau \right)& v_{yo}\left(\tau \right)& r_{xo}\left(\tau \right)& r_{yo}\left(\tau \right)\end{array}\right]}^T.$$

(1)

Here $v_{xo}\left(\tau \right)$, $v_{yo}\left(\tau \right)$, $r_{xo}\left(\tau \right)$, $r_{yo}\left(\tau \right)$ are the respective speeds and ranges of the observer in $x$ and $y$ coordinates. With the known speed and course of the observer, change in the observer’s $x$ and $y$ position is obtained as Equations (2) and (3).

$$d{r}_{xo}\left(\tau \right)=v_{xo}\left(\tau \right)\sin \left(ocr\right)\tau ,$$

(2)

$$d{r}_{yo}\left(\tau \right)=v_{yo}\left(\tau \right)\cos \left(ocr\right)\tau .$$

(3)

Here $d{r}_{xo}\left(\tau \right)$ is the change in $x$ coordinate and $d{r}_{yo}\left(\tau \right)$ is the change in $y$ coordinates of observer. $ocr$ is the observer course with $\tau$ being the time of a second. The state vector of the target is given as Equation (4).

$$S_t\left(\tau \right)={\left[\begin{array}{cccc}{v}_{xt}\left(\tau \right)& v_{yt}\left(\tau \right)& r_{xt}\left(\tau \right)& r_{yt}\left(\tau \right)\end{array}\right]}^T.$$

(4)

Here $v_{xt}\left(\tau \right)$, $v_{yt}\left(\tau \right)$, $r_{xt}\left(\tau \right)$, $r_{yt}\left(\tau \right)$ represent target’s speed and range values in $x$ and $y$ coordinates respectively. Target position, when its course and speed are known, changes according to the following Equations (5) and (6).

$$d{r}_{xt}\left(\tau \right)=v_{xt}\left(\tau \right)\sin \left(tcr\right)\tau ,$$

(5)

$$d{r}_{yt}\left(\tau \right)=v_{yt}\left(\tau \right)\cos \left(tcr\right)\tau .$$

(6)

Here $d{r}_{xt}\left(\tau \right), d{r}_{yt}\left(\tau \right)$ are the progress in $x$ and $y$ coordinates of the target. $tcr$ is the target course and $\tau$ is the time of one second. The target’s state vector relative to the observer [9,12,13] is as in Equation (7).

$$S_s\left(\tau \right)={\left[\begin{array}{cccc}{v}_x\left(\tau \right)& v_y\left(\tau \right)& r_x\left(\tau \right)& r_y\left(\tau \right)\end{array}\right]}^T.$$

(7)

Here $v_x\left(\tau \right)$, $v_y\left(\tau \right)$, $r_x\left(\tau \right)$, $r_y\left(\tau \right)$ are relative values of respective speed and range in $x$ and $y$ coordinates. The state vector of the target relative to the observer for the proceeding time, concerning the relative target state vector at present is given as Equation (8).

$$S_s\left(\tau +1\right)=A\left(\tau \right)S_s\left(\tau \right)+b\left(\tau +1\right)+\omega C\left(\tau \right).$$

(8)

Here $A\left(\tau \right)$ is the system dynamics matrix given as Equation (9) and $b\left(\tau \right)$ is the deterministic matrix given as in Equation (10).

$$A\left(\tau \right)=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ \tau & 0& 1& 0\\ 0& \tau & 0& 1\end{array}\right],$$

(9)

$$b\left(\tau +1\right)={\left[0 0 -\left(r_{xo}\left(\tau +1\right)-r_{xo}\left(\tau \right)\right) -\left(r_{yo}\left(\tau +1\right)-r_{yo}\left(\tau \right)\right)\right]}^T.$$

(10)

Here $C\left(\tau \right)$ is the system noise and its gain matrix, ω is taken as in Equation (11).

$$\omega =\left[\begin{array}{cc}\tau & 0\\ 0& \tau \\ {\tau }^2/2& 0\\ 0& {\tau }^2/2\end{array}\right]$$

(11)

Let $Q\left(\tau \right)$ be the covariance of system noise. Then it is given as in Equations (12) and (13).

$$Q\left(\tau \right)=E\left[\left(\omega C\left(\tau \right)\right){\left(\omega C\left(\tau \right)\right)}^T\right],$$

(12)

$$Q\left(\tau \right)={\sigma }^2\left[\begin{array}{cccc}{\tau }^2& 0& {\tau }^3/2& 0\\ 0& {\tau }^2& 0& {\tau }^3/2\\ {\tau }^3/2& 0& {\tau }^4/4& 0\\ 0& {\tau }^3/2& 0& {\tau }^4/4\end{array}\right].$$

(13)

Here ${\sigma }^2$ is taken as the variance of system noise multiplied with some constant value as given in (13). The measurement equation for SST has bearing angles from single sensor array (HMA) and the bearing angle $\beta \left(\tau \right)$ is given as in Equation (14).

$$\beta \left(\tau \right)={\tan }^{-1}\left(r_x\left(\tau \right)/r_y\left(\tau \right)\right).$$

(14)

Measured bearing, ${\beta }_m\left(\tau \right)$, degraded with noise is obtained and is mathematically expressed as in Equation (15).

$${\beta }_m\left(\tau \right)=\beta \left(\tau \right)+\eta \left(\tau \right).$$

(15)

Here $\eta \left(\tau \right)$ is the noise in the measurement. Here zero-mean additive process noise and measurement noise are considered. The system measurement equation is given as in Equation (16).

$$Z\left(\tau \right)=H\left(\tau \right)S_s\left(\tau \right)+\gamma \left(\tau \right).$$

(16)

Here $H\left(\tau \right)$ matrix gives the relation between measurement and state of the system and the measurement noise matrix is represented by $\gamma \left(\tau \right)$.

2.2. Tracking Using Two Sensor Arrays

For target tracking using two sensor arrays (TST), the bearing measurements from HMA and TA are considered. Consider the true bearing angles from HMA and TA be ${\beta }_{ha}\left(1\right), {\beta }_{ha}\left(2\right), {\beta }_{ha}\left(3\right), \dots , {\beta }_{ha}\left(\tau \right)$ and ${\beta }_{ta}\left(1\right), {\beta }_{ta}\left(2\right), {\beta }_{ta}\left(3\right), \dots , {\beta }_{ta}\left(\tau \right)$ respectively where $\tau$ represents the time. Measurement vector, ${\beta }_m\left(\tau \right)$ can be written as in Equation (17).

$${\beta }_m\left(\tau \right)={\left[\begin{array}{cc}{\beta }_{ham}\left(\tau \right)& {\beta }_{tam}\left(\tau \right)\end{array}\right]}^T,$$

(17)

$\tau$ represents the time for which the measurement vector is considered. The relative target state vector, that is taken concerning observer be $S_S\left(\tau \right)$, is given as in Equation (18).

$$S_S\left(\tau \right)={\left[\begin{array}{cccccc}{\dot{v}}_x\left(\tau \right)& {\dot{v}}_y\left(\tau \right)& v_x\left(\tau \right)& v_y\left(\tau \right)& r_x\left(\tau \right)& r_y\left(\tau \right)\end{array}\right]}^T.$$

(18)

Here, ${\dot{v}}_x\left(\tau \right)$, ${\dot{v}}_y\left(\tau \right)$, $v_x\left(\tau \right)$, $v_y\left(\tau \right)$, $r_x\left(\tau \right)$, and $r_y\left(\tau \right)$ are corresponding acceleration, corresponding speed and corresponding range components of the target concerning observer. The dynamic equation of the target state is given as in Equations (19) and (20).

$$S_S\left(\tau +1\right)=A\left(\tau +1\right)S_S\left(\tau \right)+b\left(\tau +1\right)+\omega \mathrm{C}\left(\tau \right),$$

(19)

$$b\left(\tau +1\right)={\left[0 0 0 0 -\left(r_{xo}\left(\tau +1\right)-r_{xo}\left(\tau \right)\right) -\left(r_{yo}\left(\tau +1\right)-r_{yo}\left(\tau \right)\right)\right]}^T.$$

(20)

Here $A$ is the transition matrix and is given as in Equation (21) and $\omega$ is taken as in Equation (22).

$$A\left(\tau +1\right)=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ \tau & 0& 1& 0& 0& 0\\ 0& \tau & 0& 1& 0& 0\\ {\tau }^2/2& 0& \tau & 0& 1& 0\\ 0& {\tau }^2/2& 0& \tau & 0& 1\end{array}\right],$$

(21)

$$\omega ={\left[\begin{array}{cccccc}1& 0& \tau & 0& {\tau }^2/2& 0\\ 0& 1& 0& \tau & 0& {\tau }^2/2\end{array}\right]}^T.$$

(22)

The plant noise $\omega \left(\tau \right)$ is assumed to follow Gaussian distribution. It is given by Equation (23).

$$C\left(\tau \right)={\left[\begin{array}{cc}{C}_{\ddot{x}}& C_{\ddot{y}}\end{array}\right]}^T.$$

(23)

The plant noise covariance matrix is given Equations (24) and (25).

$$Q\left(\tau \right)=E\left[\omega \mathrm{C}\left(\tau \right){\mathrm{C}}^{\mathrm{T}}\left(\tau \right){\omega }^T\right],$$

(24)

$$Q\left(\tau \right)=\left[\begin{array}{cccccc}1& 0& \tau & 0& {\tau }^2/2& 0\\ 0& 1& 0& \tau & 0& {\tau }^2/2\\ \tau & 0& {\tau }^2& 0& {\tau }^3/2& 0\\ 0& \tau & 0& {\tau }^2& 0& {\tau }^3/2\\ {\tau }^2/2& 0& {\tau }^3/2& 0& {\tau }^4/4& 0\\ 0& {\tau }^2/2& 0& {\tau }^3/2& 0& {\tau }^4/4\end{array}\right]\times E\left[C\left(\tau \right)C^T\left(\tau \right)\right].$$

(25)

Initially, HMA is assumed to be at origin and TA is $L$ meters apart from HMA and towed behind the ship. Figure 1 represents the initial positions of HMA, TA, and target. Let $\left(R_{x1}, R_{y1}\right)$ and $\left(R_{x2}, R_{y2}\right)$ be $x$ and $y$ coordinates of target range corresponding to HMA and TA. The observer is assumed to be moving with a course of $\alpha$ concerning the Y-axis. $B_{ta}$ and $B_{ha}$ are the true bearing angles concerning TA and HMA. $\left(R_{x3}, R_{y3}\right)$ is the location of TA concerning HMA and is assumed to be 500m apart from HMA.

From Figure 1, $R_{x3}$ and $R_{y3}$ are given by Equations (26) and (27).

$$R_{x3}=L\sin \left(\alpha +180\right),$$

(26)

$$R_{y3}=L\cos \left(\alpha +180\right).$$

(27)

From Figure 1, the relation between range coordinates can be written as in Equations (28) and (29).

$$R_{x2}=R_{x1}+R_{x3},$$

(28)

$$R_{y2}=R_{y1}+R_{y3}.$$

(29)

The HMA bearing, ${\beta }_{ham}$ is modelled as in Equation (30).

$${\beta }_{ham}\left(\tau \right)={\tan }^{-1}\left(\frac{R_{x1}\left(\tau \right)}{R_{y1}\left(\tau \right)}\right)+{\gamma }_{ha\beta }\left(\tau \right),$$

(30)

where ${\gamma }_{ha\beta }\left(\tau \right)$ is an error in the HMA measurement and is assumed to follow zero-mean Gaussian distribution with variance ${\sigma }_{ha\beta }^2\left(\tau \right)$. The measurement and plant noises are assumed to be uncorrelated to each other. Measured bearing from TA, ${\beta }_{tam}\left(\tau \right)$ is given as in Equation (31).

$${\beta }_{tam}\left(\tau \right)={\tan }^{-1}\left(R_{x2}\left(\tau \right)/R_{y2}\left(\tau \right)\right)+{\gamma }_{ta\beta }\left(\tau \right).$$

(31)

Here, error in TA measurement is represented as ${\gamma }_{ta\beta }\left(\tau \right)$ and is assumed to follow Gaussian distribution with zero mean and having variance of ${\sigma }_{ta\beta }^2\left(\tau \right)$.

The covariance matrix of the measurements, $\varnothing \left(\tau \right)$, at time $\tau$ is taken as in Equation (32).

$$\varnothing \left(\tau \right)=\left[\begin{array}{cc}{\sigma }_{ha\beta }^2\left(\tau \right)& 0\\ 0& {\sigma }_{ta\beta }^2\left(\tau \right)\end{array}\right].$$

(32)

2.3. Algorithm of Extended Kalman Filter for TST

The main idea of EKF is to linearize the existing non-linearities in measurement as well as state equations and then perform Kalman filtering. The nonlinearity in this application is considered only in the measurement Equation (16) [20]. Therefore, a partial derivation is applied to the measurement model matrix for linearization and is linearized measurement model matrix for HMA, $H_{ha}$ is obtained as follows in Equation (33).

$$H_{ha}\left(\tau +1\right)=\left[0 0 0 0 \frac{{\widehat{R}}_{y1}\left(\tau +1,\tau \right)}{{\widehat{R}}_1^2\left(\tau +1,\tau \right)} \frac{-{\widehat{R}}_{x1}\left(\tau +1,\tau \right)}{{\widehat{R}}_1^2\left(\tau +1,\tau \right)}\right].$$

(33)

Here $R_1^2=R_{x1}^2+R_{y1}^2$

The true parameters of the target are not known in real-time scenarios, so the estimated values, ${\widehat{R}}_{x1}$, and ${\widehat{R}}_{y1}$ are used in Equation (33). Measurement matrix of TA, $H_{ta}$ is given by Equation (34).

$$H_{ta}\left(\tau +1\right)=\left[0 0 0 0 \frac{{\widehat{R}}_{y2}\left(\tau +1,\tau \right)}{{\widehat{R}}_2^2\left(\tau +1,\tau \right)} \frac{-{\widehat{R}}_{x2}\left(\tau +1,\tau \right)}{{\widehat{R}}_2^2\left(\tau +1,\tau \right)}\right].$$

(34)

Here $R_2^2=R_{x2}^2+R_{y2}^2$.

Measurement fused matrix, $H$ is given by Equation (35).

$$H\left(\tau +1\right)={\left[\begin{array}{cc}{H}_{ha}\left(\tau +1\right)& H_{ta}\left(\tau +1\right)\end{array}\right]}^T.$$

(35)

The priori state vector estimation is calculated as follows in Equation (36).

$$S_s\left(\tau +1|\tau \right)=A\left(\tau \right)S_s\left(\tau |\tau \right).$$

(36)

The estimated priori state covariance matrix estimation [2] is done as follows in Equation (37).

$$P\left(\tau +1|\tau \right)=A\left(\tau \right)P\left(\tau |\tau \right){\left(A\left(\tau \right)\right)}^T+Q\left(\tau \right).$$

(37)

The gain [2] is calculated for EKF as in Equation (38).

$$G\left(\tau \right)=P\left(\tau +1|\tau \right)H^T\left(\tau \right){\left(H\left(\tau \right)P\left(\tau +1|\tau \right)H^T\left(\tau \right)+\varnothing \left(\tau \right)\right)}^{-1}.$$

(38)

The posteriori state and its covariance matrices, based on the measurements, are updated respectively as follows in Equations (39) and (40).

$$S_s\left(\tau +1|\tau +1\right)=S_s\left(\tau +1|\tau \right)+G\left(\tau \right)Z\left(\tau \right),$$

(39)

$$P\left(\tau +1|\tau +1\right)=\left(I-G\left(\tau \right)H\left(\tau \right)\right)P\left(\tau +1|\tau \right){\left(I-G\left(\tau \right)H\left(\tau \right)\right)}^T+G\left(\tau \right)\varnothing \left(\tau \right){\left(G\left(\tau \right)\right)}^T.$$

(40)

2.4. Algorithm of Modified Gain Bearings-Only Extended Kalman Filter for TST

The modified gain function, $g_{ha}$ [12] is given by Equation (41).

$$g_{ha}=\left[0 0 0 0 \left(\frac{\cos {\beta }_{ham}}{{\widehat{R}}_{x1}\sin {\beta }_{ham}+{\widehat{R}}_{y1}\cos {\beta }_{ham}}\right) \left(\frac{-\sin {\beta }_{ham}}{{\widehat{R}}_{x1}\sin {\beta }_{ham}+{\widehat{R}}_{y1}\cos {\beta }_{ham}}\right)\right],$$

(41)

and TA’s modified gain function, $g_{ta}$ can be written as in Equation (42).

$$g_{ta}=\left[0 0 0 0 \left(\frac{\cos {\beta }_{tam}}{{\widehat{R}}_{x2}\sin {\beta }_{tam}+{\widehat{R}}_{y2}\cos {\beta }_{tam}}\right) \left(\frac{-\sin {\beta }_{tam}}{{\widehat{R}}_{x2}\sin {\beta }_{tam}+{\widehat{R}}_{y2}\cos {\beta }_{tam}}\right)\right].$$

(42)

Finally, $g$ is written as in Equation (43).

$$g={\left[\begin{array}{cc}{g}_{ha}& g_{ta}\end{array}\right]}^T.$$

(43)

The predicted covariance matrix is calculated as in Equation (44).

$$P\left(\tau +1\right)=\left(A\left(\tau +1\right)P\left(\tau \right)A^T\left(\tau +1\right)\right)+\omega C\left(\tau +1\right){\mathsf{\omega }}^{\mathrm{T}}.$$

(44)

The Kalman gain is calculated as in Equation (45).

$$G\left(\tau +1\right)=P\left(\tau +1\right)H^T\left(\tau +1\right){\left[\varnothing \left(\tau \right)+H\left(\tau +1\right)P\left(\tau +1\right)H^T\left(\tau +1\right)\right]}^{-1}.$$

(45)

The updated state matrix is calculated as in Equation (46).

$$S_s\left(\tau +1\right)=S_s\left(\tau +1\right)+G\left(\tau +1\right)\left[Z\left(\tau +1\right)-Z\left(\tau +1,S_s\left(\tau +1\right)\right)\right].$$

(46)

Here $Z\left(\tau +1, S_s\left(\tau +1\right)\right)$ is the bearing measurement obtained from predicted estimate at time index $\left(\tau +1\right)$. The updated covariance matrix is given in Equation (47).

$$P\left(\tau +1\right)=\left[I-G\left(\tau +1\right)g\left({\beta }_m\left(\tau +1\right),S_s\left(\tau +1\right)\right)\right]P\left(\tau +1\right){\left[I-G\left(\tau +1\right)g\left({\beta }_m\left(\tau +1\right),S_s\left(\tau +1\right)\right)\right]}^T+\varnothing \left(\tau \right)G\left(\tau +1\right)G^T\left(\tau +1\right).$$

(47)

2.5. Particle Filter Combined with Other Filters

PF [12] is an advanced nonlinear filter, suitable for systems that are nonlinear and applicable in the non-Gaussian environment. In PF, estimation is done by considering the probability density function (PDF) of the state vector. Initially, a set of $N$ random state vectors are generated based on the initial PDF of the state, $p\left(S_S\left(0|0\right)\right)$, which is assumed to be a known value [21]. These randomly generated state vectors are called particles and are represented as $S_S\left(n,\tau |\tau \right)\left(n=1,2,\dots ,N\right)$. These particles are propagated to the next time using (48).

$$S_S\left(n,\tau +1|\tau \right)=f\left(S_S\left(n,\tau |\tau \right),C\left(\tau +1\right)\right), n=1,2,\dots , N$$

(48)

Here, $C$ is the plant noise that is generated randomly based on the PDF of the particle. Once the measurements are obtained, the conditional relative likelihood of each particle, $p\left(Z\left(\tau \right), S_S\left(n,\tau +1|\tau \right)\right)$, is computed. For an m-dimensional measurement equation, a relative likelihood $O\left(\tau \right)$ can be computed as follows in Equation (49) [12].

$$\begin{array}{c}O\left(\tau \right)=P\left[Z\left(\tau \right)=Z,S_S\left(n,\tau \right)=S_S\left(n,\tau +1|\tau \right)\right]\\ =P\left[\gamma \left(\tau \right)=Z-H\left(S_S\left(n,\tau +1|\tau \right)\right)\right]\\ ~\left(1/{\left(2\pi \right)}^{m/2}{\left({\sigma }_B^2\right)}^{m/2}\right)\times exp\left[\frac{\left(-{\left[Z-H\left(S_S\left(n,\tau +1|\tau \right)\right)\right]}^T\left[Z-H\left(S_S\left(n,\tau +1|\tau \right)\right)\right]\right)}{2{\sigma }_B^2}\right]\end{array}$$

(49)

‘$~$’ symbol in (47) represents the proportionality relation rather than the equality relation between the probability and expression on the right side. Now the relative likelihoods obtained in (46) are normalized as follows in Equation (50).

$$O\left(\tau \right)=O\left(\tau \right)/\left({\displaystyle \sum }_{i=1}^{N}O\left(i\right)\right)$$

(50)

Now, the particles using the computed likelihoods are re-sampled. This means a new set of particles are randomly generated based on the relative likelihoods $O\left(\tau \right)$. As the number of particles increases to infinity, the accuracy in estimation increases and thereby complexity of computation. At any point in time, these sets of particles can approximate the PDF of the state accurately [21]. The advanced microprocessors, available these days can easily manage the complexity in computation regarding PF.

It is a known fact that PF has the disadvantage of high computational complexity, but this disadvantage can be reduced by the utilization of fast computing processors with high storage memory (which are available now a days). PF also has the disadvantage of sample degeneracy and sample impoverishment. This is reduced by several resampling techniques available in literature [9,22]. However, these techniques are not suitable for all types of applications. Based on application, proper sampling technique must be selected.

One approach that has been proposed for improving particle filtering is to combine it with another filter such as the EKF, UKF, or MGBEKF [12]. Here, each posterior particle will be updated using the posterior updating equations of EKF, UKF, or MGBEKF algorithms and then re-sampling (if required) is performed. This is like working with an $N$ number of Kalman filters (one for each particle), that is initialized with randomly generated state vectors and then re-sampling them (if required) after updating them with each measurement obtained. In this research work, PF combined with the EKF algorithm is explained for fused measurement data in the form of a flowchart in Figure 2. Similarly, PF is combined with MGBEKF. In this paper, the particle filter is combined with the EKF, MGBEKF and the algorithms are named as particle filter coupled with extended Kalman filter (PFEKF) and particle filter coupled with modified gain bearings-only extended Kalman filter (PFMGBEKF). $S_S\left(n,\tau +1|\tau \right)$ is updated to $S_S\left(n,\tau +1|\tau +1\right)$ according to the EKF and MGBEKF equations [12,13] in PFEKF and PFMGBEKF respectively.

2.6. Cramer–Rao Lower Bound

The Cramer–Rao lower bound (CRLB) provides a lower bound on the variance of an unbiased estimator for comparing the performance of any estimator [21]. It is defined that for a nonlinear discrete-time system, the covariance of the estimated state follows the condition in Equation (51).

$$E\left[\left({\widehat{S}}_S\left(\tau \right)-S_S\left(\tau \right)\right){\left({\widehat{S}}_S\left(\tau \right)-S_S\left(\tau \right)\right)}^T\right]\ge J_{\tau }^{-1},$$

(51)

where $J_{\tau }$ is an information matrix, which is calculated recursively [11], using Equation (52).

$$J_{\tau }=D_{\tau -1}^{22}-D_{\tau -1}^{21}{\left(J_{\tau -1}+D^{11}\right)}^{-1}{D}_{\tau -1}^{12},$$

(52)

where $D^{11}, D^{12}, D^{21}and D^{22}$, for additive Gaussian noise, are given as follows in Equations (53)–(56).

$$D^{11}=A^T{Q}^{-1}$$

(53)

$$D^{21}=-A^T{Q}^{-1}$$

(54)

$$D^{12}={\left(D^{21}\right)}^T$$

(55)

$$D_{\tau -1}^{22}=Q^{-1}+E\left[H_{\tau }^T{\left({\sigma }_B^2\right)}^{-1}{H}_{\tau }\right]$$

(56)

Here ${\mathrm{H}}_{\mathsf{\tau }}$ is given by the Jacobian of the measurement function for TST as in Equation (57).

$$H_{\tau }=\left[\begin{array}{c}0 0 0 0 \cos {\widehat{\beta }}_{ham}/R -\sin {\widehat{\beta }}_{ham}/R\\ 0 0 0 0 \cos {\widehat{\beta }}_{tam}/R -\sin {\widehat{\beta }}_{tam}/R\end{array}\right]$$

(57)

and ${\widehat{\beta }}_{ham}$ is the predicted HMA measurement and ${\widehat{\beta }}_{tam}$ is the predicted TA measurement. For SST the equations are as given in [23]. The CRLB of range error is calculated for TST as in Equation (58).

$$\mathrm{CRLB}\_{\mathrm{R}}_{\mathsf{\tau }}=\sqrt{{\mathrm{J}}_{\mathsf{\tau }}^{-1}\left(5,5\right)+{\mathrm{J}}_{\mathsf{\tau }}^{-1}\left(6,6\right)}$$

(58)

Similarly, the CRLB of speed error for TST is defined as in Equation (59).

$$CRLB\_S_{\tau }=\sqrt{J_{\tau }^{-1}\left(3,3\right)+J_{\tau }^{-1}\left(4,4\right)}$$

(59)

CRLB of course error for TST is calculated as in Equation (60).

$$CRLB\_C_{\tau }={\tan }^{-1}\left(\sqrt{J_{\tau }^{-1}\left(3,3\right)}/\sqrt{J_{\tau }^{-1}\left(4,4\right)}\right)$$

(60)

3. Implementation and Simulation

Scenarios that are considered for evaluating the algorithms are as shown in

Table 1. Scenarios for the evaluation of algorithms.

Parameter	Scenario No.
	1	2	3
Initial Range (m)	3000	3000	3000
Initial Bearing (deg)	0	0	0
Target Speed (m/s)	12	12	12
Observer Speed (m/s)	8	8	8
Target course (deg)	155	145	110
Observer course (deg)	90	90	90

. The bearing measurements from HMA and TA are assumed to be available continuously for every second for TST and SST.

In an underwater environment, a submarine in a passive mode of operation can track a target at 4 to 5 kilometers. However, this changes with the environmental conditions underwater. Scenarios with different target courses are considered to track the target moving in different directions. Scenario with an initial bearing of 20° is like the scenario with an initial bearing of 0° with coordinates frame turned by 20°. Hence the bearing value is chosen as 0° for simplicity. Target is also assumed to be a submarine, so, the speed of the target is chosen accordingly [24].

For SST, the observer is presumed to maneuver following the famous S-maneuver in its course. Therefore, for two minutes, the observer follows the current course it is i.e, 90°. Now it slowly turns 180° at the rate of 0.5° per second to attain 270°. With 0.5° per second turn rate, the observer takes 6 minutes to complete the maneuver. For TST, the observer is presumed to be traveling at a constant speed and with a constant course of 90° i.e., without any maneuver in course or speed. The noise in measurements is assumed to follow Gaussian distribution and is additive. The initializations of the target state vector and its covariance matrices are chosen as in Equations (33)–(35). The speed of the target is difficult to estimate and is assumed to be 5m/s, as the submarine travels at very low speeds. Sonar range of the day, which gives the minimum range of target detection with a 50 percent probability of detection, of 5000 m is considered in target state vector initializations.

For SST, the initial state vector is given in Equation (61) and for TST it is given in Equation (62).

$$S_s\left(0|0\right)={\left[5 5 5000\sin {\beta }_m 5000\cos {\beta }_m\right]}^T$$

(61)

$$S_S\left(0|0\right)={\left[0.1 0.1 5 5 5000\sin \left({\beta }_{ham}\left(1\right)\right) 5000\sin \left({\beta }_{ham}\left(1\right)\right)\right]}^T$$

(62)

Initial target state estimate $S(0|0)$ is assumed to be distributed uniformly and elements of the state covariance matrix are initialized accordingly as Equation (63).

$$P\left(0|0\right)=Diagonal{\left(4S_S^2\left(0|0\right)/12\right)}_i$$

(63)

Here $i=1, 2,\dots , 6$ for TST and $i=1, 2,\dots , 4$ for SST

The particles of the target state vector for PFEKF and PFMGBEKF filters are generated by adding random noise to the initial state vector as in (58) for TST. Similarly, for SST the initial state vector is added with a column matrix of random numbers with four rows as in Equation (64) [22].

$${\left(S_S\left(0|0\right)\right)}_j=S_S\left(0|0\right)+{\left[\begin{array}{cccccc}randn& randn& randn& randn& randn& randn\end{array}\right]}^T$$

(64)

Here, $j$ represents the particle number and $j=1, 2, 3, \dots , 1000$. The covariance matrix for each particle of the target state vector for PFEKF and PFMGBEKF filters are generated based on each particle generated.

It is not genuine to evaluate the performance of the algorithm on a run basis, as randomness exists in the experiment. Therefore, Monte Carlo simulation of 100 runs is performed for each filtering in MATLAB [6] for both SST and TST algorithms. The maximum (3σ) error acceptable in estimated range, course and speeds are 8% of true range, 3° and 1m/s, the root-mean-squared (RMS) (1σ) error allowed for the acceptance of the solution are 2.66% of true range, 1° and 0.33m/s respectively. The errors in the estimation of target parameters are root-mean-squared and the solution is assessed on the error acceptance levels mentioned above.

RMSE is nothing but root mean squared error and it is, in general, called standard deviation. It is implemented as follows. The difference between the true target parameter and the estimated target parameter is taken as error and it is squared. The mean of the squared error for all the 100 runs is then calculated and square rooted. This provides the RMSE of the estimated target parameters. Let $i$ be the number of Monte Carlo runs and $j$ be the total time in seconds for which the simulation is carried out. Then the error is mathematically given as Equation (65) and the RMSE as Equation (66).

$$E_{\mathit{ij}}={\left(\mathit{true} \mathit{target} \mathit{parameter}\right)}_j-{\left(\mathit{estimated} \mathit{target} \mathit{parameter}\right)}_{\mathit{ij}}$$

(65)

$${\left(\mathit{RMSE}\right)}_j=\sqrt{\mathit{mean}\left({\left(E_{\mathit{ij}}\right)}^2\right)}$$

(66)

The error in the estimated target parameters must be within the vicinity range of the weapon that has to be fired onto the target. The error in estimated target parameters is calculated based on Monte Carlo runs. Therefore, the RMSE values of the estimated target parameters are considered. With Monte Carlo simulations with 100 different noise sequences (maintaining same PDF), the outputs are root mean squared to obtain a single value. The acceptance criteria are chosen based on the weapon homing zone. The solution is said to be obtained only when the error in the estimated target parameters is less than the given values in acceptance criteria. Once the solution is obtained, the estimated target parameters can be used to release weapon onto the target.

The convergence time of the solutions for the three scenarios based on the above-mentioned acceptance criteria for 100 runs is tabulated in

Table 2. Convergence times of target parameters in seconds using the SST algorithm.

Filter	Scenario No.	Target Parameter				Simulation Time
		R	C	S	TC
EKF	1	273	301	330	330	462.82
	2	NC	NC	NC	NC	239.49
	3	339	416	339	416	408.42
MGBEKF	1	NC	NC	NC	NC	287.53
	2	NC	NC	NC	NC	371.97
	3	249	321	240	321	380.34
UKF	1	216	216	216	216	329.3
	2	244	239	329	329	345.6
	3	320	338	318	338	441.23
PFEKF	1	247	247	247	247	263.85
	2	249	246	245	249	122.39
	3	343	434	340	434	236.39
PFMGBEKF	1	273	320	269	320	398.29
	2	236	337	301	337	389.29
	3	270	457	288	457	228.24

R: Range; C: Course; S: Speed; TC: Total Convergence; NC: No Convergence.

for different filtering algorithms. In a Monte Carlo simulation, it is observed that UKF converges faster than other filters for BOT. Moreover, with single bearing measurement, the observer needs to follow S-maneuver, which can be avoided while using bearing measurements from two sensors.

The simulation is carried out for 30 minutes. Results obtained with UKF for scenario 1 in Table 1 with SST and TST, are presented in Figure 3 and Figure 4 respectively. For SST, the observer needs to perform S-maneuver to obtain the observability of the process and hence to obtain a solution. With TST, the observability and solution are obtained even without observer maneuver as it uses bearing measurements from two sensors.

Table 3. Convergence times of target parameters in seconds using the TST algorithm.

Filter	S. No	Target Parameter				Simulation Time
		R	C	S	TC
EKF	1	155	200–253; 347–442	213–257; 301–532	200–253; 347–442	301.81
	2	158–447	227–447	182–447	227–447	255.68
	3	109–148; 187–828	NC	226–676	NC	441.86
MGBEKF	1	18–676; 679–730	206–296; 377–401	159–357	206–296; 386–410	1160.08
	2	17–296	NC	124–296	NC	315.04
	3	218	333–365; 465–485; 550–611; 655–675	206–246; 260–280	NC	385.01
UKF	1	193	235	230	235	427.36
	2	175	248	221	248	430.81
	3	202	333	140	333	325.57
PFEKF	1	116	337	242	337	141.40
	2	119	437	259	437	159.59
	3	114	713	200	713	198.75
PFMGBEKF	1	116	237	150	237	436.61
	2	119	256	152	256	1038.03
	3	130	712	85	712	163.65

S. No: Scenario Number; R: Range; S: Speed; C: Course; TC: Total Convergence; NC: No Convergence.

gives the total convergence times of the solution with various filtering algorithms for TST. The solution is said to be converged if the acceptance criteria are achieved in all the three parameters at a time and the solution does not diverge for the rest of the simulation time. It can be observed from the table data that the solution convergence time is only for a short period with some of the scenarios and not for the entire simulated time for EKF. For scenario 2, the solution is converged from 227 to 447th second and then the solution diverged for the rest of the simulation time. The solution is not obtained for some of the scenarios with EKF. With MGBEKF, the same trend of solution-divergence persists. With UKF, the solution is obtained earlier than other filters for most of the scenarios. As shown in Table 3, UKF for scenario 1 has convergence time of 235 seconds whereas PFEKF and PFMGBEKF have convergence times of 337 and 237 seconds respectively. It can also be observed from Table 2 and Table 3 data that the convergence times for TST are obtained faster than SST.

As the measurements are taken from two different sensor arrays and both the measurements are nonlinearly related to the target state vector elements, the process nonlinearity will be very high. Therefore, sub-optimal nonlinear filters, EKF and MGBEKF, failed to give proper solutions. Particle filter deals with complete non-linearity and is suitable for highly nonlinear systems. However, to avoid the sample degeneracy problem it is combined with other filters and hence gives better results than EKF and MGBEKF.

RMS errors in range, course, and speed estimates of the target are represented in Figure 5, Figure 6 and Figure 7. It can be observed from the figures that the RMS errors in estimated target parameters are less in case of TST than that with SST. RMS errors in estimated parameters are obtained within the acceptance criteria for TST much earlier than that of SST. However, as the range between observer and target increases, the RMS errors in estimated target parameters increases slightly for TST than that of SST.

Figure 8 shows the acceptable error in range along with CRLB values for scenario 3, according to the acceptance criteria mentioned above, which is 2.667% of true range. RMS errors in the range of the target for different filters are also shown in Figure 8. It can be observed from the figure that lower RMS error is obtained earlier with PFEKF and PFMGBEKF filters than that of EKF, MGBEKF, and UKF filters. It can be seen from Figure 8b that the lowest RMS error is obtained with PFEKF. Though the error is below the acceptable error values, there is an effect of EKF that makes the error to vary in an unstable manner, which can be seen in the form of spikes in the plot. Similarly, with PFMGBEKF, the instability in error is reduced when compared to PFEKF.

Figure 9 shows the acceptable error in the course for the same scenario 3, according to the acceptance criteria mentioned above, which is less than 1°. RMS errors in the course of the target for different filters are also shown in Figure 9. From Figure 9a, it can be observed that the error in the course is not obtained within acceptance criteria with EKF and MGBEKF. However, the fastest convergence in the target course was observed to be obtained with UKF, which can be seen from Figure 9b. PFEKF and PFMGBEKF also have convergence in the target course at 562 and 550 s respectively.

Figure 10 shows the acceptable error in speed for the same scenario 3, according to the acceptance criteria mentioned above, which is less than 0.33 m/s. RMS errors in the speed of the target for different filters are also shown in Figure 10. It can be observed from the figure that RMS error in speed with EKF and MGBEKF are varying in an unstable manner and are out of bounds after a certain time. With UKF, the convergence of speed error within acceptable criteria takes a shorter time than other filters and has a stable pattern. With PFMGBEKF, the lowest values in error are obtained. The overall convergence times were observed to be shorter for UKF than other filters and hence UKF performs well for TST than other filters.

4. Conclusions

In an underwater scenario, the observer and target are considered as submarines. Accordingly, scenarios are considered covering different target courses. SST and TST algorithms are considered for the performance evaluation concerning the convergence of the solution using different filters. Simulation is carried out and it is confirmed that TST generates the solution faster when compared to SST and even without observer maneuver. Hence, for passive target tracking, tracking with two bearing measurements is recommended. For TST, it was observed from simulated results that UKF provides a faster solution than other filtering algorithms.