Tracking an Underwater Object with Unknown Sensor Noise Covariance Using Orthogonal Polynomial Filters

In this manuscript, an underwater target tracking problem with passive sensors is considered. The measurements used to track the target trajectories are (i) only bearing angles, and (ii) Doppler-shifted frequencies and bearing angles. Measurement noise is assumed to follow a zero mean Gaussian probability density function with unknown noise covariance. A method is developed which can estimate the position and velocity of the target along with the unknown measurement noise covariance at each time step. The proposed estimator linearises the nonlinear measurement using an orthogonal polynomial of first order, and the coefficients of the polynomial are evaluated using numerical integration. The unknown sensor noise covariance is estimated online from residual measurements. Compared to available adaptive sigma point filters, it is free from the Cholesky decomposition error. The developed method is applied to two underwater tracking scenarios which consider a nearly constant velocity target. The filter’s efficacy is evaluated using (i) root mean square error (RMSE), (ii) percentage of track loss, (iii) normalised (state) estimation error squared (NEES), (iv) bias norm, and (v) floating point operations (flops) count. From the simulation results, it is observed that the proposed method tracks the target in both scenarios, even for the unknown and time-varying measurement noise covariance case. Furthermore, the tracking accuracy increases with the incorporation of Doppler frequency measurements. The performance of the proposed method is comparable to the adaptive deterministic support point filters, with the advantage of a considerably reduced flops requirement.


Introduction
Target tracking of an underwater object in a passive mode is a challenging problem owing to the nonlinear and low information content of the passive measurements, coupled with the complexities of the ocean environment [1][2][3][4][5][6]. A typical passive measurement used in such problems is bearings, and the resulting tracking process is termed bearingsonly tracking (BOT), where the objective is to estimate the position and velocity of a moving target using bearings-only measurements. Another type of underwater tracking utilises both bearings and Doppler measurements received from the target and is known as Doppler-bearing tracking (DBT). Due to their importance in underwater target tracking, both BOT and DBT have received considerable attention in the literature [7][8][9][10][11]. Since the main aim in the solution to this problem is to estimate the position and velocity of a moving target, this problem is often known as target motion analysis (TMA) [2,7].
In this paper, we assume that the object is approaching towards the ownship following a nearly constant velocity path. The sonar mounted on the ownship measures only While the percentage of track loss of the developed method is comparable to adaptive deterministic support point filters, the flops requirement of our approach is considerably less. The robustness of the proposed method is tested by assuming a target maneuver with a small and constant turn rate, while the estimators wrongly modelled it as constant velocity. To compensate for this, the estimators use a higher process noise covariance, and it has been shown that the proposed filters work with acceptable accuracy, particularly for DBT.
The remaining part of the paper is organised as follows. An underwater tracking problem is formulated in Section 2. In Section 3, we derive the proposed filter. The simulation results are discussed in Section 4, followed by a brief discussion and conclusions.

System Model
The state vector of a target at time index k is given by X t k = X t 1,k X t 2,kẊ t 1,kẊ t 2,k , where (X t 1,k , X t 2,k ) are the true positions and (Ẋ t 1,k ,Ẋ t 2,k ) are the true velocities of the target along the x and y axis, respectively, and represents the transpose of a matrix. Similarly, the state vector of the observer (which is ownship) can be expressed as where F is the state transition matrix. The process noise, η k−1 , follows a Gaussian distribution with zero mean and covariance, whereq is the process noise intensity, and T is the sampling time. We define the relative state vector as Using Equation (1), we can write the state equation as where U k−1,k = X 0 k − FX 0 k−1 is known as the vector of deterministic inputs [7]. Here, we consider a single target which follows (i) a nearly constant velocity path described by a constant velocity (CV) model, or (ii) a maneuvering trajectory described by a constant turn (CT) model with constant and known turn rate. It is interesting to note that, for both cases, the process equation is linear.

Constant Velocity Model
The state transition matrix for the CV model is given by and the vector of deterministic inputs [2,7] becomes

Constant Turn Model
In this case, it is assumed that the target maneuvers with a constant and known turn rate, ω. The state transition matrix of this model is expressed as [23,29,30] and the expression of U k−1,k becomes It is to be noted that in the limiting condition ω → 0, Equations (6) and (7) reduce to Equations (4) and (5), respectively.

Measurement Model
Measurement sensors are mounted on the ownship and all the sensors are passive so that some specific tactical advantages can be achieved. It is assumed that the sonar mounted on the ownship is only capable of measuring the bearing angle. We assume the following two cases.

Case I-Bearings-Only Measurement
The sensors measure the direction of the target location with respect to true north. Such measurements can be expressed as where γ(X k ) = tan −1 (X 1,k /X 2,k ), and ν θ k is noise in bearing measurement, which is assumed to be zero mean, white Gaussian with standard deviation, σ θ , i.e., ν θ k ∼ N (0, σ 2 θ ). For such measurements, the system is unobservable [2,7] and the estimation of the target state is only possible when the ownship starts maneuvering. It is to be noted that even if, sometimes, the sensor measures the bearing angle with respect to its own local coordinates, it is easy to transform the measurement with respect to true north.

Case II-Doppler-Bearing Measurement
As the target is moving, a Doppler shift occurs on the tonal frequency emitted from the target, which can also be measured. In such cases, both the noise-corrupted bearing angles and Doppler-shifted frequencies are available. This is popularly known as Doppler-bearing (DB) measurement and the overall measurement function γ(X k ) becomes where f T is the constant tonal frequency emitted from the target, r 2 k = X 2 1,k + X 2 2,k , c is the speed of the sound in water, and the measurement noise vector in this case is Note that ν f k is the sensor noise in the Doppler measurement and assumed to be zero mean, white, Gaussian with variance σ 2 f , i.e., ν f k ∼ N (0, σ 2 f ). Further, we assume that the noises in bearing and Doppler measurements is uncorrelated, with overall measurement noise covariance R k = diag(σ 2 θ , σ 2 f ). It is to be noted that for Doppler-bearing (DB) measurements, the system is observable and maneuvering of the ownship is not required for estimators to converge [9]. Further, it is expected that if noise variance for Doppler measurement is reasonably low, the Doppler-bearing tracking (DBT) should provide more accurate results than bearings-only tracking (BOT).

Tracking Methodology
In this section, we formulate a tracking filter for underwater TMA when sensor noise covariance is unknown. The proposed filter linearises the nonlinear measurement function using a first-order orthogonal polynomial approximation with a Gaussian measure as a weight function and estimates the measurement noise covariance R k at each step.

Orthogonal Polynomial Filter
The proposed orthogonal polynomial filter is a kind of linearised Kalman filter, where, unlike the extended Kalman filter, the linearisation is not done using the Taylor series; rather, it is done with a first-order Hermite polynomial. The method was developed for a general nonlinear state estimation problem in our earlier publication [28]. Here, we provide a more straightforward derivation of the filter for an underwater tracking problem.
The problem considered here consists of a linear process model, so the time update step can be done by the Kalman filter (KF). The available measurements are assumed to be (i) bearings, (ii) bearings with Doppler-shifted frequencies. Thus, our measurement Y k ∈ R n y , where n y is 1 or 2 for bearings and Doppler-bearing measurements, respectively. As the measurement equation is nonlinear, it is linearised with the first-order Hermite polynomial, and subsequently the KF scheme is applied for the measurement update.
Let us assume an arbitrary function of Gaussian random variable (r.v.) X , f (X ), which maps X to a real space of arbitrary dimension n p , i.e., f (X ) : R n x → R n p is approximated with a first-order Hermite polynomial [28,[31][32][33]; so, where the coefficients A 0 and A 1 are matrices with dimension n p × 1 and n p × n x , respectively. H 0 (X ) and H 1 (X ) are the zero and first-order Hermite polynomial, respectively. Please note that in Equation (10), we neglect higher-order terms so that the function is approximated by a linear equation. All the terms with higher-order Hermite polynomials that are neglected in Equation (10) combined together form the residual error. We choose the Hermite polynomial since its weighting function expansion is the same as the Gaussian [31,32]. As the Hermite polynomial is orthogonal [33], where <, > represents the dot product, X ∼ N (X , P), and I n x is the identity matrix of dimension n x . The coefficient A i can be calculated as It is easy to verify that H 0 (X ) = 1 and H 1 (X ) = S −1 (X −X ), where P = SS . Substituting the values of H 0 (X ) and H 1 (X ) in Equation (10), the approximation of Following Equation (13), we can linearise the measurement function γ(X k ) as where the state X k ∈ R 4 follows a Gaussian distribution with meanX k|k−1 and covariance P k|k−1 , and the coefficients A k and B k are n y × 1 and n y × 4 matrices, respectively. S k|k−1 is the square root of the covariance matrix P k|k−1 (i.e., P k|k−1 = S k|k−1 S k|k−1 ), which can be calculated by the Cholesky decomposition. Using Equation (12), the coefficient matrices A k and B k can be expressed as The above integrals cannot be solved analytically for any arbitrary nonlinear function [34]. Various methods such as the Gauss-Hermite (GH) quadrature rule [21], unscented transform [19,20], cubature quadrature (CQ) rule [23,24], and high-degree cubature quadrature (HDCQ) rule [35,36] are available to solve such integrals. Here, in brief, we discuss the GH quadrature rule.
Gauss-Hermite quadrature rule: A single-dimensional Gauss-Hermite (GH) quadrature rule [21,34] is given by where ξ i are the quadrature points with associated weights ω i for i = 1, . . . , m. The above integral can be extended for a multidimensional variable using the product rule [21]. To calculate these quadrature points, we take a symmetric tridiagonal matrix J, with zero diagonal elements and J i,i+1 = √ i/2, 1 ≤ i ≤ (m − 1) [21,37]. The quadrature points are located at 1 is the first element of the i-th normalised eigenvector of J. Thus, a multidimensional integral (18) can be approximately evaluated as [38] The total number of required points in the above method for any n x dimensional system is m n x . Here, we see that the number of points in the GH quadrature rule increases exponentially with the system's dimension. For more about the point and weight generation methods of the UKF, CQKF, HDCQKF, readers are referred to [19,21,24,36].
Calculation of A k and B k : To calculate A k and B k (Equations (15) and (16)) using the above-discussed GH quadrature rule, we need to transform Equations (15) and (16) into standard Gaussian integrals. Transforming the Gaussian r.v. X k into a standard Gaussian r.v. x k following X k =X k|k−1 + S k|k−1 x k , the integrals (Equations (15) and (16)) become Using the points (ξ i ) and weights (ω i ), the coefficients A k (Equation (20)) can be expressed as where For bearings-only measurements given in Equation (8), where χ 1,i represents the first element of the i-th transformed quadrature point. For the DB measurements given in Equation (9), the expressions of A k and B k become

Remark 1.
Even with the assumption of Gaussian pdf of states, the approximation arises due to (i) neglecting the terms containing second-and higher-order Hermite polynomials; (ii) approximate evaluation of the coefficients A 0 and A 1 by the summation method. However, if the nonlinear function is in polynomial form, the exact evaluation of coefficients can be done with the Gaussian integral [17].

Remark 2.
On neglecting higher-order Hermite polynomial terms in the approximation (Equation (10)), an accumulative error may occur. It can be compensated by the following approaches [12] (pp. 395-402): (i) increasing the process noise artificially; (ii) multiplying the error covariance by a factor slightly greater than unity. In both cases, the Kalman gain increases and more weight is given to the present innovation, thus compensating for the error which occurs during linearisation.

Remark 3.
From the above discussion, we have seen that A k and B k can be evaluated using the weighted sum of deterministic sample points. If we use the GH quadrature rule [21,34] to obtain the points, the estimation method is called Gauss-Hermite orthogonal polynomial-based EKF (GHO-EKF). The extension EKF is used because the proposed method linearises the nonlinear function (but not with Taylor series expansion). Similarly, when the unscented transform [19] and CQ rule [23,24] are used to generate the points, the filters are called unscented orthogonal polynomial-based EKF (UO-EKF) and cubature orthogonal polynomial-based EKF (CO-EKF), respectively.

Measurement Update
As we have linearised the measurement equation, the measurement update step can be performed with the Kalman filter (KF) algorithm. The estimated value of the measurement and its covariance can be calculated aŝ and The cross-covariance between the state and measurement can be calculated as Remark 4. Although, throughout our work, we consider the measurement noises to be Gaussian, it may not always be so [30,39,40]. In the case of non-Gaussian noise, it can be approximated by a Gaussian distribution and we can apply the proposed method. However, in this case, since there is a mismatch of noise statistics with the assumed one, robust filters might be helpful. Further, the noise can also be represented as a sum of Gaussians, and Gaussian sum filters [41] with the proposed method may help to achieve the desired results.

Estimation of Measurement Noise Covariance
The measurement noise considered here is time-varying, zero mean white Gaussian with unknown covariance. In the literature, filtering for such problems is known as adaptive filtering [42][43][44][45][46][47][48]. Several adaptive filters have been proposed based on the KF [42,44], EKF [48], UKF [46,49], CKF [50], and GHF [46]. A few works on the TMA using BOT for unknown time-varying measurement noise are found in [51,52]. In this paper, we have developed a noise adaptation technique using the orthogonal polynomial filter for both BOT and DBT. Here, we derive an expression for online R k estimation compatible with the orthogonal polynomial filter and presented in Lemma 1.

Lemma 1.
The covariance of measurement noise, R k can be expressed as where P Y Y k|k is the covariance of residual measurement, P k|k is the posterior error covariance, and b k = B k S −1 k|k−1 Proof. The measurement equation using Equation (14) can be written as The expression for the posterior estimate of measurement iŝ Subtracting Equation (29) from Equation (28), we obtain the residual of measurement As X k ,X k|k−1 , and k|k−1 are uncorrelated with measurement noise ν k , Using these relations, we can write the following where the Kalman gain Now, Substituting Equations (31)- (33) in Equation (30), we have Equation (27).
From the measurement data sequence, P Y Y k|k is calculated using the following equation: where i 0 is the user's choice to fix a window length of (k − i 0 + 1). Substituting Equation (34) into Equation (27), The detailed algorithm to implement the proposed adaptive orthogonal polynomialbased filter is presented in Algorithm 1.

Algorithm 1 : Pseudo code for adaptive polynomial filter for target tracking
Step 1: Initialisation 1 Initialise the filter withX 0|0 , P 0|0 , and R 0 . 2 Compute the sample points ξ i , and their respective weights ω i .
Step 2: Time update 3 Compute the predicted mean and covariance aŝ Step 3: Measurement update 4 ComputeŶ k|k−1 , P Y Y k|k−1 , and P X Y k|k−1 using Equations (24)-(26). 5 Calculate the Kalman gain, 7 Compute the posterior error covariance, Step 4: Estimate unknown R k 8 Calculate the posterior estimate of the measurement,Ŷ k|k using Equation (29). 9 Calculate the residual of measurement as Y k −Ŷ k|k . 10 Calculate the covariance of measurement residual (P Y Y k|k ) using Equation (34). 11 Estimate R k using Equation (35).

Alternative to the Square Root Filters
In sigma point filters, to calculate the sigma point at each time step, the Cholesky factorisation of the covariance matrix (P), i.e., P = SS , is required. The accumulated round-off error associated with the limited word length arithmetics of the processing software [25][26][27] sometimes leads to an asymmetric, non-positive definite covariance matrix, which forces the filter to stop. To circumvent such a problem, the square root sigma point filters such as square root UKF (SRUKF), square root quadrature Kalman filter (SRQKF), square root CKF (SRCKF), etc., were developed [23,[25][26][27]. In square root filtering, at each time step, the square root of P, i.e., S, is propagated directly using QR decomposition. The QR decomposition guarantees positive definite covariance matrices at the expense of increased computational burden.
The orthogonal polynomial filters are free from such problems and less sensitive to the round-off error and preserve the symmetry and positive definiteness property [28]. Furthermore, the computational cost of this filter is less than the square root filters, which we show using the flops count in the next subsection.

Flops Count
We analyze the computational complexity of the filters by counting the number of floating point operations (flops) [26,28,53]. A flop is defined as one basic algebraic operation, such as addition, subtraction, multiplication, and division, between any two floating point numbers [26,53]. The filter algorithm consists of the following matrix operations: (i) addition, (ii) multiplication, (iii) inverse, and (iv) Cholesky decomposition. The addition of two matrices with dimension m × n requires mn flops. The multiplication of an m × n matrix with n × p matrix requires mp(2n − 1) number of flops [53]. The inverse and Cholesky operation of any matrix A ∈ R n×n require n 3 and n 3 /3 + 2n 2 flops, respectively [53]. The QR decomposition of any matrix with dimension m × n using the Householder method requires 2mn 2 − 2/3n 3 flops.
We calculate the flops count for all the filters (both ordinary and adaptive) for BOT and DBT measurements, and these are summarised in Table 1. From the table, we see that the filters with DBT measurements require more flops than those with the BOT. The adaptation of the unknown measurement noise covariance brings an extra computational burden on the filtering algorithm. Among all filters, the EKF has the lowest flops count. The orthogonal polynomial filters such as CO-EKF, UO-EKF, and GHO-EKF have lower flops counts than their square root counterparts (SRCKF, SRUKF, and square root GHF (SRGHF)).

Tracking Scenarios
In this paper, we consider two tracking scenarios, as shown in Figure 1. From the figure, we see that in both scenarios, the target follows a nearly constant velocity path. In the first engagement scenario (Figure 1a), the ownship maneuvers smoothly during the period 61-420 s with a turn rate of 0.4 • /s. In the second scenario (Figure 1b), the observer sharply maneuvers at t = 181 s and changes its course to 320 • from 180 • . In the second scenario, the rate of change for the bearing angle during the maneuver is high, which means that it is more difficult for any suboptimal filter to track. It is well known that the maneuvering of ownship makes the system observable. We took two scenarios from the earlier literature [2,7] as benchmark problems so that we could demonstrate the efficacy of our proposed algorithm against the background of existing available results. For both scenarios, the estimators use (i) bearings-only and (ii) bearings and Dopplershifted frequency measurements, and tracking is performed for 15 min. The speed of sound in water (c) is taken to be 1.5 km/s, and the tonal frequency ( f T ) of the target is 385 Hz. Although, in an underwater scenario, the tonal frequency may not be available continuously and it may change after some time interval, here, we assume that the tonal frequency is fixed and available throughout the simulation. The parameters used in the simulation are provided in Table 2.

Performance Metrics
We evaluate the performance of the tracking filters using the following performance metrics.

Root Mean Square Error
The position root mean square error (RMSE) at any time index k can be expressed as where M is the number of ensembles, and (X i,t 1,k , X i,t 2,k ) is the position at the k-th time step of i-th ensemble with corresponding estimated value (X i,t 1,k ,X i,t 2,k ).

Normalised (State) Estimation Error Squared
The normalised (state) estimation error squared (NEES) at the k-th time index can be defined as [12] (p. 234) If we assume that the estimators are consistent, and the system is linear with Gaussian noise, then the NEES follows a chi-square distribution with n x degree of freedom, and E[NEES k ] = n x , the dimension of the state. The average value of the NEES (ANEES) for M Monte Carlo (MC) runs can be calculated as [12] (p. 234) [54] The estimators are called consistent if the ANEES k ∈ [l b , u b ], where l b and u b are the lower and upper bounds of the acceptance interval, respectively. Furthermore, the estimator is considered optimistic [12] (p. 245) if ANEES k > u b , because, in this case, the value of P k|k is too small, whereas the estimator is pessimistic [12] (p. 245) when the value of ANEES k is lower than l b . The calculation of l b and u b can be found in [28,55] and references therein.

Bias Norm
The position bias norm at the k-th time step for M number of MC runs is expressed as [56] Bias where · 2 denotes the vector norm, X p,k = (X 1,k , X 2,k ) is the true position, and X p,k = (X 1,k ,X 2,k ) is its posterior estimate.

Track Loss
The terminal error of estimation is defined as where X t 1,kmax andX t 1,kmax are the true and estimated x-position at the final time step, with similar definitions for these variables with subscript 2, which correspond to the y-position. A track is deemed lost if the terminal error (e b ) is above a threshold, e T .

Performance Comparison
All the filters, namely the orthogonal polynomial filters and deterministic sample point filters, are implemented with BOT and DBT measurements for Scenario I and Scenario II in the following cases: (i) known and time-invariant measurement noise (we call it R known ), (ii) known and time-varying measurement noise (we call it R k,known ), (iii) unknown and time-invariant measurement noise (we denote it R unknown ), (iv) unknown and time-varying measurement noise (we denote it R k,unknown ). It is needless to say that we replace the filters with their adaptive counterparts whenever we encounter unknown measurement noise covariance. For the R k,unknown case, we vary R k from (1.5 • ) 2 to (4 • ) 2 linearly based on the range of target from ownship. Results are summarised in figures and tables showing the RMSE, ANEES, bias norm, track loss, and computational time for various cases. Throughout the simulation, for the SRUKF, we set κ = 1 as the choice κ = 3 − n x was found to be problematic due to the negative weights. However, such problems did not appear in UO-EKF and so, for this filter, we set κ = −1. For SRGHF and GHO-EKF, we used three univariate points (i.e., m = 3), which resulted in a total of 3 4 = 81 support points.  Figure 3a with 95% probability regions, for which the lower bound and upper bound become l b = 0.939 and u b = 1.0629, respectively [28,55]. From this figure, we see that all the filters' ANEES values fall inside the concentration region only by the end of the simulation time. Position bias norms obtained from 500 MC runs are plotted in Figure 3b. From this figure, we see that the bias norms of all filters are sufficiently low, and among them, the SRGHF and GHO-EKF attain the lowest bias norm.  To show the efficacy of the proposed measurement noise adaptation technique, in Figure 4, we plot R unknown and R k,unknown and their estimates obtained from the adaptive GHO-EKF (AD-GHO-EKF). From this figure, we see that the estimated noise covariances converge to their true values. Very similar estimates are obtained from other adaptive polynomial filters and are not plotted in the figure. To observe the effect of unknown R in state estimation (both BOT and DBT cases), position and velocity RMSEs (excluding track loss cases, with error threshold e T = 500 m) obtained from a single representative filter, say AD-GHO-EKF, are plotted in Figure 5a,b, respectively. Here, we also provide the RMSE results for the DBT case in order to show the improvement in RMSE that we could achieve in comparison to BOT if Doppler shifts are incorporated as measurements. The findings of the result for the DBT case are explicitly discussed in Scenario I: Case II. The results are compared with GHO-EKF for both the time-varying and time-invariant unknown measurement noise covariance. From these figures, we see that the GHO-EKF provides lower RMSEs compared to AD-GHO-EKF in this situation. Furthermore, the RMSE for the R unknown case is slightly less than R k,unknown , as expected.  ANEES metrics are calculated for various adaptive filters for e T ∈ {200 m, 500 m}. For e T = 500 m, none of the filters' ANEES is inside the concentration region. In Figure 6a, the ANEES of various adaptive filters are plotted from 500 MC runs e T = 200 m, and for this case, we notice that the ANEES of the GHO-EKF (DBT, R known ) is consistently inside the probability region, whereas the same obtained from other filters fall inside the region only after 650 s. The poor ANEES results are mainly due to low observability and unknown timevarying measurement noise covariance. Furthermore, we provide the position bias norm of GHO-EKF and AD-GHO-EKF for unknown measurement noise covariance obtained from 500 MC runs in Figure 6b. From this figure, we observe that GHO-EKF (DBT) has the lowest bias norm, whereas the AD-GHO-EKF for R k,unknown attains the highest.
Next, we study the effect of sampling time on the accuracy of estimation. In Figure 7, we plot the position and velocity RMSE (excluding the track divergence with e T = 500 m) of AD-GHO-EKF obtained from 500 MC runs for different sampling times, T, ranging from 1 s to 60 s. From these figures, we observe that for lower sampling time T, the AD-GHO-EKF exhibits better performance, as expected. However, it is to be noted that the sampling time of the sensor cannot be very small as certain time is required to accumulate the energy received in order to measure the bearing angle.  (Figures 8 and 9), we can conclude that the filters with DBT measurements track more effectively than those with BOT measurements.    The efficacy of the proposed adaptive estimator is also checked for R unknown using DBT measurements. We provide the position and velocity RMSEs of GHO-EKF and AD-GHO-EKF (R unknown ) obtained from 500 MC runs (excluding the track divergence cases, e T = 500 m)) in Figure 5. From this figure, we see that GHO-EKF attains lower RMSEs than AD-GHO-EKF. We also see that the incorporation of the Doppler-shifted frequency improves the filters' estimation accuracy. Now, we compare the performance of filters (EKF, orthogonal polynomial filters, deterministic square root sample point filters, and their adaptive counterparts) for BOT and DBT in terms of percentage track loss obtained from 10,000 MC runs, presented in Table 3. The track divergence bound is set to e T = 500 m. From this table, we see that the filters with DBT measurements have lower track loss than those with BOT. As expected, filters with known R show better performance than those with unknown R, i.e., adaptive filters. The orthogonal polynomial filters provide almost similar track loss as their square root counterparts. We also observe that the filters with time-varying R k (known and unknown) show higher track loss than time-invariant R, and these results are not provided. Furthermore, we observe that the track loss count of all the filters increases with the increase in sampling time. This happens due to the fact that with the increase in sampling time, the observability of the system decreases, which further deteriorates the filters' performance. For BOT, the position and velocity RMSE of different filters for R known from 500 MC runs are plotted in Figure 10a,b, respectively. From these figures, we see that except for the EKF, all other filters attain similar RMSEs. Figure 11a shows that except EKF, the ANEES of all filters (with e T = 500 m) fall within the provided probability region. The position bias norm plot (Figure 11b) shows that all the filters attain a sufficiently low bias norm.   Figure 11. (a) ANEES; (b) position bias norm of the BOT for Scenario II.

Scenario II : Case II-Doppler-Bearing Tracking
For DBT, the position and velocity RMSE of different filters for R known obtained from 500 MC runs are plotted in Figure 12a,b, respectively. From this figure, we see that the RMSEs of all filters are comparable and lower than the RMSEs obtained from BOT.
The estimation accuracy of the proposed AD-GHO-EKF for R unknown for BOT and DBT is examined in Figure 13 in terms RMSEs (obtained from 500 MC runs). For comparison, we also include the RMSEs of GHO-EKF for R known . From this figure, it can be seen that the GHO-EKF provides better estimation accuracy than the AD-GHO-EKF for both cases of measurements. We also see that the RMSE with DBT measurements is lower than the RMSE obtained from the BOT.  Now, we calculate the percentage track loss of the EKF, square root filters, orthogonal polynomial filters, and their adaptive counterparts out of 10,000 MC runs for BOT and DBT, the results of which are presented in Table 3. The track loss bound is set to be e T = 500 m. From this table, it can be seen that the filters with DBT measurements perform better than those with BOT. We also see that for known measurement noise covariance, the filters have lower track loss. It can be also observed that the orthogonal polynomial filters have almost similar track loss to the square root filters.
Finally, we compare the filtering performance in terms of their execution time. The relative execution times of all filters with respect to the EKF for both BOT and DBT are reported in Table 4. From this table, we see that the execution times of the GHO-EKF and SRGHF are almost three to four times those of the EKF. The CO-EKF and UO-EKF take nearly 50% more execution time than the EKF, whereas the computation time of the SRCKF and SRUKF is nearly twice that of the EKF. We also observe that for both BOT and DBT measurements, the orthogonal polynomial filters have lower computation times than their respective square root deterministic sample filters. As expected, the adaptation of measurement noise covariance further increases the execution time of the filters. To assess how our proposed algorithms work in the presence of a target maneuver, we consider a target that is maneuvering with the constant turn rate (CT) model, as described in Section 2.1.2, with known ω = −0.015 • /s, which is shown in Figure 1a. Initially, we assume that the estimators know the target maneuvering model and check the performance of all the filters. We see that all the filters perform with similar accuracy, as described in Figures 2-13. Further, to study the robustness of these filters, we assume that the target is following a maneuver as described above, but the estimators are unaware of this dynamics and they use the constant velocity model to estimate the position and velocity of the target. To cope with such model mismatch, the process noise intensity of the filtering algorithm is increased to 100q. The target and all the estimators are initialised as described in Section 4.3. Initially, we experiment with a 500 m track loss threshold, i.e., e T = 500 m, and it has been observed that the final track error in all the filters in this model mismatched case was above this threshold in most of the runs. Further, we calculate the percentage of track loss for all the filters considering e T = 1500 m and list them in Table 5. From this table, we note that DBT shows much less track loss than BOT and the adaptive filters in DBT show less than 1% track loss, while their non-adaptive counterparts exhibit zero track loss. However, for BOT, the percentage of track loss obtained from adaptive filters is almost double that of their non-adaptive counterparts. In Figure 14, we plot the RMSE of position and velocity obtained from the adaptive and ordinary GHO-EKF for both the bearings-only and Doppler-bearing tracking. Comparing Figure 5 with Figure 14, it has been observed that the position RMSE for the model mismatch case is converging to a slightly higher value than that of the the constant velocity scenario. Further, the RMSE obtained from DBT is lower compared to BOT, and the same obtained with GHO-EKF is lower than its adaptive counterpart.

Discussion and Conclusions
In this paper, we have proposed an adaptive orthogonal polynomial-based filter that can estimate the target trajectories along with the measurement noise error covariance. It is based on the linearisation of the measurement function using the first-order orthogonal polynomial approximation, and the coefficients of these polynomials are evaluated using the numerical integration technique. The proposed method is used to estimate a constant velocity target trajectory for two different tracking scenarios, each with bearings-only and Doppler-bearing measurements. The estimation accuracy of the filters is studied in terms of RMSE. We also use the percentage track loss, computation time, average NEES, and bias norm to compare the performance of the filters. The robustness of the proposed algorithm is studied when there is a mismatch between the target and the filter-assumed dynamics. From the simulation results, it can be deduced that the proposed filters consistently outperform the Taylor series-based EKF and provide comparable performance with their square root counterparts of the deterministic sample point filters with a reduced computational cost. Moreover, the proposed filters with CV model target dynamics can estimate the track of a lightly maneuvering target with acceptable accuracy. Target motion analysis with unknown sensor misalignment remains under the scope of future works.