# Support Vector Regression-Based Adaptive Divided Difference Filter for

Authors: Hongjian Wang; Jinlong Xu; Aihua Zhang; Cun Li; Hongfei Yao
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

We present a support vector regression-based adaptive divided difference filter (SVRADDF) algorithm for improving the low state estimation accuracy of nonlinear systems, which are typically affected by large initial estimation errors and imprecise prior knowledge of process and measurement noises. The derivative-free SVRADDF algorithm is significantly simpler to compute than other methods and is implemented using only functional evaluations. The SVRADDF algorithm involves the use of the theoretical and actual covariance of the innovation sequence. Support vector regression (SVR) is employed to generate the adaptive factor to tune the noise covariance at each sampling instant when the measurement update step executes, which improves the algorithm’s robustness. The performance of the proposed algorithm is evaluated by estimating states for (i) an underwater nonmaneuvering target bearing-only tracking system and (ii) maneuvering target bearing-only tracking in an air-traffic control system. The simulation results show that the proposed SVRADDF algorithm exhibits better performance when compared with a traditional DDF algorithm.


---

## Body

## 1. Introduction

The problem of state estimation for nonlinear systems has been a subject of considerable research interest in recent years, but there is still no single solution that outperforms all other approaches. Most proposed estimators are nonlinear extensions of the dominated Kalman filter (see [1]), and each approach provides a suboptimal trade-off between properties such as numerical robustness, computational burden, and estimation accuracy. The extended Kalman filter (EKF), which linearizes both nonlinear terms of a current estimated state trajectory, is based on a first-order Taylor series and displays poor performance if the system is highly nonlinear. The limitations of EKFs are enumerated in [2]. Another improved algorithm is the iterated extended Kalman filter (IEKF), which linearizes the nonlinear model around an updated state rather than the predicted state (see [3]). Although IEKFs have been proven to perform better than EKFs in addition to globally guaranteeing convergence, the algorithm still requires a Jacobian matrix just like EKFs. However, no solution exists for the Jacobian matrix in nonlinear systems for some situations, which limits the potential application of both EKFs and IEKFs.In recent years, a new class of filter known as sigma-point Kalman filter (SPKF) has attracted a great deal of attention. In SPKFs, the algorithm propagates a cluster of points centered on the current state instead of linearizing the system dynamics to improve the approximations of the conditional mean and covariance. Unscented Kalman filters (UKF) and divided difference filters (DDF) are two kinds of SPKFs.UKFs use a deterministic sampling technique to pick a minimal set of sample points around the mean to catch the higher order statistics of the system so as to better estimation accuracy and convergence characteristics (see [4]). In [5, 6], a UKF for a class of nonlinear discrete-time systems with correlated noises was designed to deal with the problem of nonlinear filtering failure found in conventional UKFs when system noise is correlated with measurement noise. The proposed UKF breaks the limitation of conventional UKFs that requires system noise and measurement noise to be uncorrelated Gauss white noises, thus extending the potential application of conventional UKFs. In [7], a UKF filtering algorithm with colored measurement noise was proposed. The algorithm was first derived on the basis of augmented measurement information and minimum mean square error estimation, and a filtering recursive formula of UKF with colored noise then added by applying an unscented transformation to calculate the posterior mean and covariance of the nonlinear state within the optimal framework. The proposed UKF effectively dealt with the fact that traditional UKFs fail when measurement noise is colored. In [8], a UKF was applied to multiple target tracking, with the proposed UKF shown to have improved performance versus previous EKF approaches. In [9], the paper discussed an adaptive multiuser receiver for CDMA systems in which the scaled unscented filter (SUF) and the square root unscented filter (SURF) were used for joint estimation and tracking of the code delays and multipath coefficients of the received CDMA signals. The proposed channel estimators were more near-far resistant than in conventional EKFs and presented lower complexity than conventional particle filter- (PF-) based methods. Computer simulation results demonstrated the superior performance of the proposed channel estimators, and the proposed estimators were shown to exhibit lower complexity relative to the PF-based method. Although UKFs have undergone a significant amount of meaningful theory innovation and are now used in many fields, [10–12] show that UKF accuracy is lower than that of DDF, while also having a higher computational cost.The divided difference filter (DDF) first proposed by Nøgaard et al. (see [13]) linearizes the nonlinear terms based on Stirling’s interpolation polynomial approximations formula rather than Taylor’s approximation of nonlinear terms in an EKF. Conceptually, the implementation principle resembles that of an EKF; however, the DDF is significantly simpler as it does not need to calculate the Jacobian matrix and no derivatives are required. The DDF that Nøgaard et al. developed works on general discrete-time nonlinear models in which the noises are not assumed to be additive. In [14], the paper further formulated a DDF in terms of the innovation vector approach, the additive process, and the measurement noise sources. In [15], the paper proposed a new filter named the maximum likelihood-based iterated divided difference filter (MLIDDF), which improved the low state estimation accuracy of nonlinear state estimation that results from large initial estimation errors and the nonlinearity of the measurement equations. Simulation results showed that the MLIDDF algorithm possessed better state estimation accuracy and a faster convergence rate. In [16], the authors proposed a novel adaptive version of the DDF that was applicable to nonlinear systems with a linear output equation. In order to make the filter robust to modeling errors, upper bounds on the state covariance matrix were derived. The parameters of the upper bound were then estimated using a combination of offline tuning and online optimization with a linear matrix inequality constraint, which ensured that the predicted output error covariance was larger than the observed output error covariance. Simulation results demonstrated the superior performance of the proposed filter as compared to the standard DDF. Reference [17] presented an ensemble-based approach that handled nonlinearity based on a simplified divided difference approximation through Stirling’s interpolation formula. The algorithm used Stirling’s interpolation formula to evaluate the statistics of the background ensemble during the prediction step, employing the formula in an ensemble square root filter (EnSRF) at the filtering step to update the background for analysis. In this sense, the algorithm is a hybrid of Stirling’s interpolation formula and the EnSRF method, while the computational cost of the algorithm is less than that of EnSRF.Different studies have focused on the application of DDFs to nonlinear state estimation problems. In [18], time delay and channel gain estimation for multipath fading code division multiple access (CDMA) signals using a DDF were investigated, and the simulation results showed that the DDF was simpler to implement and more resilient to near-far interference in CDMA networks compared with an EKF. In [19, 20], the relative kinematic states of a reentry vehicle obtained from noisy seeker measurements using a DDF were examined. The results were compared to those obtained using an EKF and a UKF and showed that the DDF was more accurate than estimators based on a Taylor approximation like the EKF. Reference [21] investigated the possibility of using a DDF for estimating the internal variables of a synchronous generator, such as the rotor angle where the data acquired is from a phasor measurement. The effectiveness of the method was tested on a single machine infinite bus system, a nine-bus system, and a 68-bus New England-New York interconnected system. In [22], a DDF using orientation estimation was considered. The fourth element of the quaternion error vector was removed from the system states to alleviate estimated error covariance matrix divergence. The measurement system was a MARG sensor, which consisted of a triaxial rate gyro, a triaxial accelerometer, and a triaxial magnetometer. The nonlinear measurement model was obtained based on the principals of operation of the magnetometer and accelerometer and the properties of the quaternion vector space. The performance of three filters, DDF, EKF, and UKF, was compared with different sampling frequencies. The work showed that the tested DDF and the UKF were more robust than the EKF under the same initial angle-error conditions. The DDF also performed better than the UKF, although the computational load for the UKF was less. In [23], a DDF-based data fusion algorithm was presented, which utilized the complementary noise profile of rate gyros and gravimetric inclinometers to extend their limits and achieve more accurate attitude estimates. In [24], a DDF-based ballistic target tracking system for the reentry phase was proposed. The paper compared DDF, EKF, and UKF algorithms using a Monte Carlo simulation approach, with the simulation results showing that the DDF outperformed both the EKF and UKF in terms of estimation accuracy and filtering credibility. In [25], a DDF with quaternion-based dynamic process modeling was applied to global positioning system (GPS) navigation to increase navigation estimation accuracy at high-dynamic regions while preserving precision at low-dynamic regions. Some properties and performance metrics were assessed and compared to those using EKF and UKF approaches.Despite their recent popularity, DDF algorithms require that both the system model and the stochastic information must be accurate. However, this condition cannot be satisfied in many practical situations, which forces the filter to adapt itself to changing conditions. One of the problems with this requirement is that any change in the process introduces measurement noise covariance. In this work, we make use of the theoretical and actual covariance of the innovation sequence, employing SVR to generate the scale factor to tune the noise covariance at each sampling instant when the measurement update step is executed to adapt the filtering algorithm. This paper is organized as follows. Section2 briefly introduces DDF theory and the proposed SVR-based adaptive strategy. Passive target tracking is then carried out to evaluate the performance of DDF and SVRADDF algorithms using a Monte Carlo simulation in Section 3. Finally, conclusions are provided in Section 4.


## 2. Development of the Support Vector Regression-Based Adaptive Divided Difference Filter

### 2.1. Divided Difference Filter

Consider the nonlinear function(1)y=f(x),
where x∈Rnx and y∈Rny. If the function is analytic, then the first-order Taylor series expanded about some point x=x¯ becomes
(2)y=((Δx1∂∂x1+Δx2∂∂x2+⋯+Δxn∂∂xn)f(x)f(x¯+Δx)=f(x¯)+DΔxf=f(x¯)iiiii+(Δx1∂∂x1+Δx2∂∂x2+⋯+Δxn∂∂xn)f(x))|x=x¯
with (2) truncated after the first-order term. Note that (2) can achieve a better local approximation if more terms are included. However, such an expanded Taylor series requires derivatives and cannot be fulfilled in some situations. Stirling’s interpolation formula is based on a finite number of evaluations of the function and does not require derivatives, with the first-order approximation yielding
(3)y=f(x¯+Δx)=f(x¯)+D~Δxf(x¯)ii=f(x¯)+1h(∑i=1nΔxiμiδi)f(x¯),
where h denotes a selected interval length and δ and μ are determined by
(4)δif(x¯)=f(x¯+h2ei)-f(x¯-h2ei),μif(x¯)=12(f(x¯+h2ei)+f(x¯-h2ei)),
with ei being the ith unit vector.Figure1 compares the results found by using (2) and (3). The function example is f(x)=ex, where h=0.56. From the figure, we can see that Stirling’s interpolation provides better accuracy than the Taylor series under the same order approximations.Comparison of first-order Stirling series with first-order Taylor series results.
(a)
(b)We now assume that the variablex has a Gaussian density with mean x¯ and covariance Px. We can introduce a transformation matrix Sx which we select as a square Cholesky factor of Px, such that Px=SxSxT. To illustrate how others can be derived, we introduce the linear transformation of x:
(5)z=Sx-1x.This linear transformation results in a stochastic decoupling ofx as the elements of z become mutually uncorrelated (see [13]). This changes (3) to
(6)y=f(Sxz)=f~(z)=f~(z¯)+D~Δzf~,
where D~Δzf~ is determined by
(7)D~Δzf~=1h(∑i=1n‍Δziμiδi)f(z¯).The meany¯, covariance Pyy, and cross covariance Pxy of y are obtained from
(8)y¯=E[y]=E[f~(z¯)+D~Δzf~]=E[f~(z¯)]=f(x),Pyy=E[(y-y¯)(y-y¯)T]iiiii=E[(D~Δzf~)(D~Δzf~)T]iiiii=14h2∑i=1n‍(f(z¯+hei)-f(z¯-hei))iiiiiiiiii×(f(z¯+hei)-f(z¯-hei))Tiiiii=14h2∑i=1n‍(f(x¯+hsx,i)-f((x¯-hsx,i)))iiiiiiiiii×(f(x¯+hsx,i)-f((x¯-hsx,i)))T,Pxy=E[(x-x¯)(y-y¯)T]iiiii=E[SxΔz(D~Δzf~)T]iiiii=12h∑i=1n‍sx,i(f(x¯+hsx,i)-f(x¯-hsx,i))T,
where sx,i is the ith column of the matrix Sx.Consider the following nonlinear dynamic system with states to be estimated:(9)xk+1=f(xk)+ωkiiiyk=h(xk)+νk,
where ωk and νk are assumed to be independent and identically distributed and independent of current and past states, such that ωk~N(o,Qk) and νk~N(0,Rk).The DDF takes the same predictor-corrector structures in the EKF and can be described as follows.Step 1 (initialization).
Suppose the state distribution atk instant is xk~N(x^k,Pk), where x^k and Pk are obtained by
(10)x^k=E[xk],Pk=E[(x-xk)(x-xk)T].Step 2 (square Cholesky factorizations).
Consider the following:(11)Pk=SxSxTSxx^=12h{fi(x^k+hsx,j)-fi(x^k-hsx,j)},
where sx,i is the jth column of Sx.Step 3 (state and covariance propagation).
One has(12)x^k+1-=f(x^k)(13)Pk+1-=Sxx^(Sxx^)T+Qkiiiiiii=Sx-(Sx-)T(14)Syx^=12h{hi(x^k+1-+sx,j-)-hi(x^k+1--sx,j-)},
where x^k+1- is the predicted state and Pk+1- is the predicted covariance matrix.Step 4 (observation and innovation covariance propagation).
Consider(15)y^k+1-=h(x^k+1-)(16)Pk+1νν=(Syx^)(Syx^)T+Rk(17)Pk+1xy=Sx-(Syx^)T,
where y^k+1- is the predicted observation vector, Pk+1νν is the innovation covariance matrix, and Pk+1xy is the cross correlation matrix.Step 5 (update).
Consider the following:(18)κk+1=Pk+1xy(Pk+1νν)-1Pk+1+=Pk+1--κk+1Pk+1ννκk+1Tυk+1=yk+1-y^k+1-x^k+1+=x^k+1-+κk+1υk+1,
where κk+1 is the gain, Pk+1+ is the updated covariance matrix, υk+1 is the innovation vector, and x^k+1+ is the updated estimated state.


### 2.2.  Adaptive DDF Algorithm

As stated previously, the DDF algorithm assumes a complete prior knowledge of the process and the measurement noise statisticsQk and Rk. However, Qk and Rk are unknown in most applications, and incorrect priori noise statistics can lead to performance degradation or even divergence for the solution. One of the effective ways to overcome this weakness is to use an algorithm to adapt the noise statistics. In this paper, we propose using a support vector regression adaptive scheme of the DDF to adjust Qk and Rk, respectively.


#### 2.2.1. Support Vector Regression

The principles underlying support vector regression (SVR) developed by Vapnik and are presented in several works (see [26, 27]).Given the train setTS={(x1,y1),…,(xm,ym)}∈(Rn×R)m, where xi∈Rn, yi∈R, i=1,2,…,m, the SVR problem can be defined as solving for the nonlinear function g(x) about x∈Rn to construct a relationship between the output and an arbitrary input x:
(19)y(x,ω)=∑i=1m‍ωig(xi)+b=(ω·g(x))+b.In [26], the author shows that changing the regression estimate minimizes the risk functional by using the following form:
(20)y(x,α)=∑i=1m‍(αi*-αi)(g(xi)·g(x))+biiiiiiiiiii=∑i=1m‍(αi*-αi)K(xi,x)+b,
where αi* and αi are Lagrange multipliers that satisfy the condition αi,αi*≥0, αiαi*=0 and K(xi,x) is a kernel function that satisfies Mercer’s condition. In this paper, we use a translation invariant Gaussian kernel, that is, =exp⁡(-∥xi-x∥2/2σ2).


#### 2.2.2.  Adaptive Scheme Based on SVR (Qk Is Fixed)

The covariance matrixRk represents the accuracy of the measurement instrument. If we assume that the noise covariance Qk is completely known, then we can derive the SVR algorithm to estimate the measurement noise covariance Rk by defining an adaptive factor Δrk to get the form:
(21)Rk=ΔrkR,
where R is the constant noise covariance matrix.This work uses such an SVR algorithm to derive the adaptive factor at time instantk, so as to estimate the value of Rk during the algorithm’s execution.The innovation sequenceυk+1 has a theoretical covariance Pk+1νν, as (16) shows. The actual residual covariance Pk+1νν can be approximated using its sample covariance by averaging inside a moving window of size N, where
(22)Ck=1N∑i=k-N+1k‍υk+1υk+1T.IfCk differs from Pk+1νν, then the diagonal elements of Rk can be adjusted to minimize the difference as much as possible. The size of this difference is given by
(23)DOMk=diag⁡(Pk+1νν-Ck)∈Rn,
where the function diag⁡ denotes the diagonal elements of the matrix. If the elements of DOMk are rounded with zero, the covariance matrix Rk is likely accurate; otherwise there may be a large deviation between values.To adjustRk, we define the SVR train set TSΔr as TSΔr={(DOMk-N+1,Δrk-N+1),…,(DOMk,Δrk)}. With this form, we can then create a corresponding function about the adaptive factor Δrk and DOMk based on (19), such that when DOMk+1 comes, then Δrk+1 can be resolved by the function.To avoid having any element within the train set be too large to affect the accuracy ofΔrk, we use a normalized function of the train set TS; that is, TSΔr′={(DOMk-N+1′,Δrk-N+1′)…,(DOMk′,Δrk′)}, where DOMk′ and rk′ can be expressed as follows:
(24)DOMk′=1N{∑i=1N‍DOMk-i+1,j}j=1,2,…n,Δrk′=1N∑i=1N‍Δrk-i+1.


#### 2.2.3.  Adaptive Scheme Based on SVR (Rk Is Fixed)

Assuming that the noise covarianceRk is completely known, we can derive an SVR algorithm to estimate the measurement noise covariance Qk. From (13), (14), and (16), we can deduce that a change in Qk will affect the covariance matrix Pk+1νν; if we increase Qk, then Pk+1νν also increases. We can adjust Qk in the SVR by deliberately mismatching Pk+1νν and Ck+1.We first define an adaptive factorΔrk, where Rk has the following form:
(25)Qk=ΔqkQ,
where Q is the constant noise covariance matrix.We can then define the SVR train setTSΔq as TSΔq={(DOMk-N+1,Δrk-N+1),…,(DOMk,Δrk)}. At this point, the solution process is the same as for solving Δrk.


## 2.1. Divided Difference Filter

Consider the nonlinear function(1)y=f(x),
where x∈Rnx and y∈Rny. If the function is analytic, then the first-order Taylor series expanded about some point x=x¯ becomes
(2)y=((Δx1∂∂x1+Δx2∂∂x2+⋯+Δxn∂∂xn)f(x)f(x¯+Δx)=f(x¯)+DΔxf=f(x¯)iiiii+(Δx1∂∂x1+Δx2∂∂x2+⋯+Δxn∂∂xn)f(x))|x=x¯
with (2) truncated after the first-order term. Note that (2) can achieve a better local approximation if more terms are included. However, such an expanded Taylor series requires derivatives and cannot be fulfilled in some situations. Stirling’s interpolation formula is based on a finite number of evaluations of the function and does not require derivatives, with the first-order approximation yielding
(3)y=f(x¯+Δx)=f(x¯)+D~Δxf(x¯)ii=f(x¯)+1h(∑i=1nΔxiμiδi)f(x¯),
where h denotes a selected interval length and δ and μ are determined by
(4)δif(x¯)=f(x¯+h2ei)-f(x¯-h2ei),μif(x¯)=12(f(x¯+h2ei)+f(x¯-h2ei)),
with ei being the ith unit vector.Figure1 compares the results found by using (2) and (3). The function example is f(x)=ex, where h=0.56. From the figure, we can see that Stirling’s interpolation provides better accuracy than the Taylor series under the same order approximations.Comparison of first-order Stirling series with first-order Taylor series results.
(a)
(b)We now assume that the variablex has a Gaussian density with mean x¯ and covariance Px. We can introduce a transformation matrix Sx which we select as a square Cholesky factor of Px, such that Px=SxSxT. To illustrate how others can be derived, we introduce the linear transformation of x:
(5)z=Sx-1x.This linear transformation results in a stochastic decoupling ofx as the elements of z become mutually uncorrelated (see [13]). This changes (3) to
(6)y=f(Sxz)=f~(z)=f~(z¯)+D~Δzf~,
where D~Δzf~ is determined by
(7)D~Δzf~=1h(∑i=1n‍Δziμiδi)f(z¯).The meany¯, covariance Pyy, and cross covariance Pxy of y are obtained from
(8)y¯=E[y]=E[f~(z¯)+D~Δzf~]=E[f~(z¯)]=f(x),Pyy=E[(y-y¯)(y-y¯)T]iiiii=E[(D~Δzf~)(D~Δzf~)T]iiiii=14h2∑i=1n‍(f(z¯+hei)-f(z¯-hei))iiiiiiiiii×(f(z¯+hei)-f(z¯-hei))Tiiiii=14h2∑i=1n‍(f(x¯+hsx,i)-f((x¯-hsx,i)))iiiiiiiiii×(f(x¯+hsx,i)-f((x¯-hsx,i)))T,Pxy=E[(x-x¯)(y-y¯)T]iiiii=E[SxΔz(D~Δzf~)T]iiiii=12h∑i=1n‍sx,i(f(x¯+hsx,i)-f(x¯-hsx,i))T,
where sx,i is the ith column of the matrix Sx.Consider the following nonlinear dynamic system with states to be estimated:(9)xk+1=f(xk)+ωkiiiyk=h(xk)+νk,
where ωk and νk are assumed to be independent and identically distributed and independent of current and past states, such that ωk~N(o,Qk) and νk~N(0,Rk).The DDF takes the same predictor-corrector structures in the EKF and can be described as follows.Step 1 (initialization).
Suppose the state distribution atk instant is xk~N(x^k,Pk), where x^k and Pk are obtained by
(10)x^k=E[xk],Pk=E[(x-xk)(x-xk)T].Step 2 (square Cholesky factorizations).
Consider the following:(11)Pk=SxSxTSxx^=12h{fi(x^k+hsx,j)-fi(x^k-hsx,j)},
where sx,i is the jth column of Sx.Step 3 (state and covariance propagation).
One has(12)x^k+1-=f(x^k)(13)Pk+1-=Sxx^(Sxx^)T+Qkiiiiiii=Sx-(Sx-)T(14)Syx^=12h{hi(x^k+1-+sx,j-)-hi(x^k+1--sx,j-)},
where x^k+1- is the predicted state and Pk+1- is the predicted covariance matrix.Step 4 (observation and innovation covariance propagation).
Consider(15)y^k+1-=h(x^k+1-)(16)Pk+1νν=(Syx^)(Syx^)T+Rk(17)Pk+1xy=Sx-(Syx^)T,
where y^k+1- is the predicted observation vector, Pk+1νν is the innovation covariance matrix, and Pk+1xy is the cross correlation matrix.Step 5 (update).
Consider the following:(18)κk+1=Pk+1xy(Pk+1νν)-1Pk+1+=Pk+1--κk+1Pk+1ννκk+1Tυk+1=yk+1-y^k+1-x^k+1+=x^k+1-+κk+1υk+1,
where κk+1 is the gain, Pk+1+ is the updated covariance matrix, υk+1 is the innovation vector, and x^k+1+ is the updated estimated state.


## 2.2.  Adaptive DDF Algorithm

As stated previously, the DDF algorithm assumes a complete prior knowledge of the process and the measurement noise statisticsQk and Rk. However, Qk and Rk are unknown in most applications, and incorrect priori noise statistics can lead to performance degradation or even divergence for the solution. One of the effective ways to overcome this weakness is to use an algorithm to adapt the noise statistics. In this paper, we propose using a support vector regression adaptive scheme of the DDF to adjust Qk and Rk, respectively.


### 2.2.1. Support Vector Regression

The principles underlying support vector regression (SVR) developed by Vapnik and are presented in several works (see [26, 27]).Given the train setTS={(x1,y1),…,(xm,ym)}∈(Rn×R)m, where xi∈Rn, yi∈R, i=1,2,…,m, the SVR problem can be defined as solving for the nonlinear function g(x) about x∈Rn to construct a relationship between the output and an arbitrary input x:
(19)y(x,ω)=∑i=1m‍ωig(xi)+b=(ω·g(x))+b.In [26], the author shows that changing the regression estimate minimizes the risk functional by using the following form:
(20)y(x,α)=∑i=1m‍(αi*-αi)(g(xi)·g(x))+biiiiiiiiiii=∑i=1m‍(αi*-αi)K(xi,x)+b,
where αi* and αi are Lagrange multipliers that satisfy the condition αi,αi*≥0, αiαi*=0 and K(xi,x) is a kernel function that satisfies Mercer’s condition. In this paper, we use a translation invariant Gaussian kernel, that is, =exp⁡(-∥xi-x∥2/2σ2).


### 2.2.2.  Adaptive Scheme Based on SVR (Qk Is Fixed)

The covariance matrixRk represents the accuracy of the measurement instrument. If we assume that the noise covariance Qk is completely known, then we can derive the SVR algorithm to estimate the measurement noise covariance Rk by defining an adaptive factor Δrk to get the form:
(21)Rk=ΔrkR,
where R is the constant noise covariance matrix.This work uses such an SVR algorithm to derive the adaptive factor at time instantk, so as to estimate the value of Rk during the algorithm’s execution.The innovation sequenceυk+1 has a theoretical covariance Pk+1νν, as (16) shows. The actual residual covariance Pk+1νν can be approximated using its sample covariance by averaging inside a moving window of size N, where
(22)Ck=1N∑i=k-N+1k‍υk+1υk+1T.IfCk differs from Pk+1νν, then the diagonal elements of Rk can be adjusted to minimize the difference as much as possible. The size of this difference is given by
(23)DOMk=diag⁡(Pk+1νν-Ck)∈Rn,
where the function diag⁡ denotes the diagonal elements of the matrix. If the elements of DOMk are rounded with zero, the covariance matrix Rk is likely accurate; otherwise there may be a large deviation between values.To adjustRk, we define the SVR train set TSΔr as TSΔr={(DOMk-N+1,Δrk-N+1),…,(DOMk,Δrk)}. With this form, we can then create a corresponding function about the adaptive factor Δrk and DOMk based on (19), such that when DOMk+1 comes, then Δrk+1 can be resolved by the function.To avoid having any element within the train set be too large to affect the accuracy ofΔrk, we use a normalized function of the train set TS; that is, TSΔr′={(DOMk-N+1′,Δrk-N+1′)…,(DOMk′,Δrk′)}, where DOMk′ and rk′ can be expressed as follows:
(24)DOMk′=1N{∑i=1N‍DOMk-i+1,j}j=1,2,…n,Δrk′=1N∑i=1N‍Δrk-i+1.


### 2.2.3.  Adaptive Scheme Based on SVR (Rk Is Fixed)

Assuming that the noise covarianceRk is completely known, we can derive an SVR algorithm to estimate the measurement noise covariance Qk. From (13), (14), and (16), we can deduce that a change in Qk will affect the covariance matrix Pk+1νν; if we increase Qk, then Pk+1νν also increases. We can adjust Qk in the SVR by deliberately mismatching Pk+1νν and Ck+1.We first define an adaptive factorΔrk, where Rk has the following form:
(25)Qk=ΔqkQ,
where Q is the constant noise covariance matrix.We can then define the SVR train setTSΔq as TSΔq={(DOMk-N+1,Δrk-N+1),…,(DOMk,Δrk)}. At this point, the solution process is the same as for solving Δrk.


## 2.2.1. Support Vector Regression

The principles underlying support vector regression (SVR) developed by Vapnik and are presented in several works (see [26, 27]).Given the train setTS={(x1,y1),…,(xm,ym)}∈(Rn×R)m, where xi∈Rn, yi∈R, i=1,2,…,m, the SVR problem can be defined as solving for the nonlinear function g(x) about x∈Rn to construct a relationship between the output and an arbitrary input x:
(19)y(x,ω)=∑i=1m‍ωig(xi)+b=(ω·g(x))+b.In [26], the author shows that changing the regression estimate minimizes the risk functional by using the following form:
(20)y(x,α)=∑i=1m‍(αi*-αi)(g(xi)·g(x))+biiiiiiiiiii=∑i=1m‍(αi*-αi)K(xi,x)+b,
where αi* and αi are Lagrange multipliers that satisfy the condition αi,αi*≥0, αiαi*=0 and K(xi,x) is a kernel function that satisfies Mercer’s condition. In this paper, we use a translation invariant Gaussian kernel, that is, =exp⁡(-∥xi-x∥2/2σ2).


## 2.2.2.  Adaptive Scheme Based on SVR (Qk Is Fixed)

The covariance matrixRk represents the accuracy of the measurement instrument. If we assume that the noise covariance Qk is completely known, then we can derive the SVR algorithm to estimate the measurement noise covariance Rk by defining an adaptive factor Δrk to get the form:
(21)Rk=ΔrkR,
where R is the constant noise covariance matrix.This work uses such an SVR algorithm to derive the adaptive factor at time instantk, so as to estimate the value of Rk during the algorithm’s execution.The innovation sequenceυk+1 has a theoretical covariance Pk+1νν, as (16) shows. The actual residual covariance Pk+1νν can be approximated using its sample covariance by averaging inside a moving window of size N, where
(22)Ck=1N∑i=k-N+1k‍υk+1υk+1T.IfCk differs from Pk+1νν, then the diagonal elements of Rk can be adjusted to minimize the difference as much as possible. The size of this difference is given by
(23)DOMk=diag⁡(Pk+1νν-Ck)∈Rn,
where the function diag⁡ denotes the diagonal elements of the matrix. If the elements of DOMk are rounded with zero, the covariance matrix Rk is likely accurate; otherwise there may be a large deviation between values.To adjustRk, we define the SVR train set TSΔr as TSΔr={(DOMk-N+1,Δrk-N+1),…,(DOMk,Δrk)}. With this form, we can then create a corresponding function about the adaptive factor Δrk and DOMk based on (19), such that when DOMk+1 comes, then Δrk+1 can be resolved by the function.To avoid having any element within the train set be too large to affect the accuracy ofΔrk, we use a normalized function of the train set TS; that is, TSΔr′={(DOMk-N+1′,Δrk-N+1′)…,(DOMk′,Δrk′)}, where DOMk′ and rk′ can be expressed as follows:
(24)DOMk′=1N{∑i=1N‍DOMk-i+1,j}j=1,2,…n,Δrk′=1N∑i=1N‍Δrk-i+1.


## 2.2.3.  Adaptive Scheme Based on SVR (Rk Is Fixed)

Assuming that the noise covarianceRk is completely known, we can derive an SVR algorithm to estimate the measurement noise covariance Qk. From (13), (14), and (16), we can deduce that a change in Qk will affect the covariance matrix Pk+1νν; if we increase Qk, then Pk+1νν also increases. We can adjust Qk in the SVR by deliberately mismatching Pk+1νν and Ck+1.We first define an adaptive factorΔrk, where Rk has the following form:
(25)Qk=ΔqkQ,
where Q is the constant noise covariance matrix.We can then define the SVR train setTSΔq as TSΔq={(DOMk-N+1,Δrk-N+1),…,(DOMk,Δrk)}. At this point, the solution process is the same as for solving Δrk.


## 3. Monte Carlo Simulation Results and Discussion

In this section, we report the experimental results obtained by applying SVRADDF to the nonlinear state estimation of a nonmaneuvering target in an underwater tracking control scenario and a maneuvering target in an air-traffic control scenario. To demonstrate the performance of the SVRADDF algorithm, we compare its performance against a DDF algorithm.


### 3.1. Underwater Nonmaneuvering Target Bearing-Only Tracking Control Scenario

We consider a bearing-only tracking control scenario, where an underwater target executes a uniform motion in a horizontal plane but unknown velocity, while a passive sonar platform performs a uniform circular motion in a horizontal plane. Figure2 shows a representative trajectory of the target and the passive sonar platform. The kinematics of the relative motion between the target and the platform can be modeled using the following linear process equation:
(26)xk+1=[10T0010T00100001]xk+wk.Figure 2
Target trajectory (I-initial position and F-final position) and sonar trajectory.Here, the state of the equation isx=[xyx˙y˙], where x and y denote position and x˙ and y˙ denote velocity in the x and y directions, respectively, and T is the time interval between two consecutive measurements and the process noise wk~N(0,Q) with a nonsingular covariance where
(27)Q=q1×[T22T22TT].In (27), the parameter q1 is related to process noise intensities. The measurement equation is written as follows:
(28)yk=θk=tan-1(ykxk)+υk,
where the measurement noise υk~N(0,R) with a nonsingular covariance.Given the following initial conditions:(29)T=1sq1=0.0001m2s-3iR=0.02mrad,the true initial state is(30)x0=[0m1500m0ms-10ms-1]T
and the associated covariance is
(31)P0=diag⁡[100m22000m21m2s-21m2s-2].The initial estimate statex^0 was chosen randomly from N(x0,P0) in each run, and the total number of scans per run was 1000.To provide a fair comparison, we performed 50 independent Monte Carlo runs. To track the underwater target, we used both the SVRADDF and the DDF algorithms and compared their performance. The adaptive factor was set toΔqk=Δrk=0.1. Both of the filters were initialized with the same initial conditions for each run.Performance metrics: to compare the nonlinear performance of the filters, we used the root mean square error (RMSE) of the target position and velocity. The RMSE yields a combined measure of the bias and variance of a filter estimate. The RMSE of the position at timek was found using
(32)RSMEpos(k)=1N∑i=1N‍((xki-x^ki)2+(yki-y^ki)2),
where (xki,yki) and (x^ki,y^ki) are the true and estimated positions, respectively, in the ith Monte Carlo run. The form for the RMSE of the velocity is similar.Figures3 and 4 show the estimated RMSE in target position and velocity. The SVRADDF uses SVR to adjust the adaptive factor during algorithm execution, which leads to a marginally better performance compared to the DDF, as seen in the figures.Figure 3
RMSE in position for DDF and SVRADDF.Figure 4
RMSE in velocity for DDF and SVRADDF.


### 3.2. Maneuvering Target Tracking in the Air-Traffic Control Scenario

A typical air-traffic control scenario was considered next, where an aircraft executes a maneuvering turn in a horizontal plane at a constant and known turn rateΩ. Figure 5 shows a representative trajectory of the aircraft.Figure 5
Aircraft trajectory and radar location.The kinematics of the turning motion can be modeled using the following nonlinear process equation:(33)xk=[1sinΩTΩ0-1-cos⁡ΩTΩ0cos⁡ΩT0-sinΩT01-cos⁡ΩTΩ1sinΩTΩ0sinΩT0cos⁡ΩT]iiiiiii×xk-1+wk-1.The state of the aircraft is given byx=[xx˙yy˙]T, where x and y denote position and x˙ and y˙ denote velocity in the x and y directions, respectively, and T is the time interval between two consecutive measurements and the process noise wk~N(0,Q) with a nonsingular covariance, where
(34)Q=q1×[T33T2200T22T0000T33T2200T22T].The parameterq1 related to process noise intensities. A passive radar is fixed at the origin and equipped to measure the bearing θ. The measurement equation is written as follows:
(35)yk=θk=tan-1(ykxk)+υk,
where the measurement noise υk~N(0,R).The parameters used in this simulation wereT=1s, Ω=-3°s-1, q1=0.1m2s-3, and R=10mard. The true initial state of the aircraft was x0=[1000m300ms-11000m0ms-1]T and the associated covariance matrix was P0=diag⁡[100m210m2s-2100m210m2s-2]. The initial state x^0 for the filters was chosen randomly from N(x0,P0) in each Monte Carlo run, and the simulation time per run was 1000.For a fair comparison, we performed 100 independent Monte Carlo runs for each filter. To track the maneuvering aircraft, we used both the SVRADDF and the DDF algorithms and compared their performance. The adaptive factor was set toΔqk=0.3 and Δrk=0.5. Both filters were initialized with the same initial conditions for each run.Figures6 and 7 show the estimate RMSE in target position and velocity for the SVRADDF and DDF filters. Not surprisingly, both filters exhibit divergence due to a mismatch between the initial filter design assumption and the Gaussian noise nature of the problem. The SVRADDF filter exhibits marginally better performance compared to the DDF since it was able to adjust the statistical properties of the noise during execution.Figure 6
RMSE in position for DDF and SVRADDF.Figure 7
RMSE in velocity for DDF and SVRADDF.


## 3.1. Underwater Nonmaneuvering Target Bearing-Only Tracking Control Scenario

We consider a bearing-only tracking control scenario, where an underwater target executes a uniform motion in a horizontal plane but unknown velocity, while a passive sonar platform performs a uniform circular motion in a horizontal plane. Figure2 shows a representative trajectory of the target and the passive sonar platform. The kinematics of the relative motion between the target and the platform can be modeled using the following linear process equation:
(26)xk+1=[10T0010T00100001]xk+wk.Figure 2
Target trajectory (I-initial position and F-final position) and sonar trajectory.Here, the state of the equation isx=[xyx˙y˙], where x and y denote position and x˙ and y˙ denote velocity in the x and y directions, respectively, and T is the time interval between two consecutive measurements and the process noise wk~N(0,Q) with a nonsingular covariance where
(27)Q=q1×[T22T22TT].In (27), the parameter q1 is related to process noise intensities. The measurement equation is written as follows:
(28)yk=θk=tan-1(ykxk)+υk,
where the measurement noise υk~N(0,R) with a nonsingular covariance.Given the following initial conditions:(29)T=1sq1=0.0001m2s-3iR=0.02mrad,the true initial state is(30)x0=[0m1500m0ms-10ms-1]T
and the associated covariance is
(31)P0=diag⁡[100m22000m21m2s-21m2s-2].The initial estimate statex^0 was chosen randomly from N(x0,P0) in each run, and the total number of scans per run was 1000.To provide a fair comparison, we performed 50 independent Monte Carlo runs. To track the underwater target, we used both the SVRADDF and the DDF algorithms and compared their performance. The adaptive factor was set toΔqk=Δrk=0.1. Both of the filters were initialized with the same initial conditions for each run.Performance metrics: to compare the nonlinear performance of the filters, we used the root mean square error (RMSE) of the target position and velocity. The RMSE yields a combined measure of the bias and variance of a filter estimate. The RMSE of the position at timek was found using
(32)RSMEpos(k)=1N∑i=1N‍((xki-x^ki)2+(yki-y^ki)2),
where (xki,yki) and (x^ki,y^ki) are the true and estimated positions, respectively, in the ith Monte Carlo run. The form for the RMSE of the velocity is similar.Figures3 and 4 show the estimated RMSE in target position and velocity. The SVRADDF uses SVR to adjust the adaptive factor during algorithm execution, which leads to a marginally better performance compared to the DDF, as seen in the figures.Figure 3
RMSE in position for DDF and SVRADDF.Figure 4
RMSE in velocity for DDF and SVRADDF.


## 3.2. Maneuvering Target Tracking in the Air-Traffic Control Scenario

A typical air-traffic control scenario was considered next, where an aircraft executes a maneuvering turn in a horizontal plane at a constant and known turn rateΩ. Figure 5 shows a representative trajectory of the aircraft.Figure 5
Aircraft trajectory and radar location.The kinematics of the turning motion can be modeled using the following nonlinear process equation:(33)xk=[1sinΩTΩ0-1-cos⁡ΩTΩ0cos⁡ΩT0-sinΩT01-cos⁡ΩTΩ1sinΩTΩ0sinΩT0cos⁡ΩT]iiiiiii×xk-1+wk-1.The state of the aircraft is given byx=[xx˙yy˙]T, where x and y denote position and x˙ and y˙ denote velocity in the x and y directions, respectively, and T is the time interval between two consecutive measurements and the process noise wk~N(0,Q) with a nonsingular covariance, where
(34)Q=q1×[T33T2200T22T0000T33T2200T22T].The parameterq1 related to process noise intensities. A passive radar is fixed at the origin and equipped to measure the bearing θ. The measurement equation is written as follows:
(35)yk=θk=tan-1(ykxk)+υk,
where the measurement noise υk~N(0,R).The parameters used in this simulation wereT=1s, Ω=-3°s-1, q1=0.1m2s-3, and R=10mard. The true initial state of the aircraft was x0=[1000m300ms-11000m0ms-1]T and the associated covariance matrix was P0=diag⁡[100m210m2s-2100m210m2s-2]. The initial state x^0 for the filters was chosen randomly from N(x0,P0) in each Monte Carlo run, and the simulation time per run was 1000.For a fair comparison, we performed 100 independent Monte Carlo runs for each filter. To track the maneuvering aircraft, we used both the SVRADDF and the DDF algorithms and compared their performance. The adaptive factor was set toΔqk=0.3 and Δrk=0.5. Both filters were initialized with the same initial conditions for each run.Figures6 and 7 show the estimate RMSE in target position and velocity for the SVRADDF and DDF filters. Not surprisingly, both filters exhibit divergence due to a mismatch between the initial filter design assumption and the Gaussian noise nature of the problem. The SVRADDF filter exhibits marginally better performance compared to the DDF since it was able to adjust the statistical properties of the noise during execution.Figure 6
RMSE in position for DDF and SVRADDF.Figure 7
RMSE in velocity for DDF and SVRADDF.


## 4. Conclusions

This work has proposed and developed an innovation-based SVRADDF algorithm. The algorithm introduces an adaptive factor estimated using SVR, which allows for estimation of the noise statistical characteristics of nonlinear stochastic systems. The SVRADDF algorithm avoids instability and divergence in the solution which is caused by incorrect statistical characteristics of the noise. Monte Carlo simulation results of an underwater nonmaneuvering bearing-only target tracking system and a maneuvering target bearing-only tracking system in an air-traffic control setting showed that the SVRADDF algorithm provides better state estimation accuracy than a traditional DDF algorithm.Although the SVRADDF algorithm has showed better performance under Monte Carlo simulation, there are still several challenging issues to be considered for future study. To improve its feasibility and effectiveness in the complex environment, more comprehensive and detailed studies are still needed to solve under the nonlinear and non-Gaussian noise conditions. Since the algorithm has only been tested under Monte Carlo simulation, the following work might be testing the algorithm under trial data to approve its better performance online.



---

*Source: 139503-2014-05-25.xml*