# A New Sparse Gauss-Hermite Cubature Rule Based on Relative-Weight-Ratios for Bearing-Ranging Target Tracking

Authors: Lijun Peng; Xiaojun Duan; Jubo Zhu
Journal: Modelling and Simulation in Engineering
 (2017)

---

## Abstract

A new sparse Gauss-Hermite cubature rule is designed to avoid dimension explosion caused by the traditional full tensor-product based Gauss-Hermite cubature rule. Although Smolyak’s quadrature rule can successfully generate sparse cubature points for high dimensional integral, it has a potential drawback that some cubature points generated by Smolyak’s rule have negative weights, which may result in instability for the computation. A relative-weight-ratio criterion based sparse Gauss-Hermite rule is presented in this paper, in which cubature points are kept symmetric in the input space and corresponding weights are guaranteed to be positive. The generation of the new sparse cubature points set is simple and meaningful for practice. The difference between our new sparse Gauss-Hermite cubature rule and other cubature rules is analysed. Simulation results show that, compared with Kalman filter with those types of full tensor-product based Gauss-Hermite rules, our new sparse Gauss-Hermite cubature rule based Kalman filter can lead to a substantially reduced number of cubature points, more stable computation capability, and maintaining the accuracy of integration at the same time.


---

## Body

## 1. Introduction

Bayesian recursive estimation is commonly used in target tracking, positioning, and signal processing [1]. Generally, a Bayesian recursive estimation algorithm requires a state model and a measurement model. The a posteriori density function can describe the behaviour of the estimated state [2]. Since a closed form solution to the Bayesian recursive estimation is available only for a few special cases [3], such as the linear Gaussian system (which leads to the classical standard Kalman filter), a suboptimal solution is a preferable choice in the general case [4, 5]. Ito and Xiong [4] suggest that a local Gaussian filter can be used to approximate the general Bayesian recursive estimation suboptimally. The core problem of local Gaussian filters is in fact a high dimensional Gaussian weighted integration, which has been studied in numerical analysis; see [6–9] and the references therein. Lots of researches concentrate on the numerical approximation methods to solve the Gaussian weighted integral problem, and different approaches result in different local Gaussian filters, such as the cubature Kalman filter [10], the quadrature Kalman filter [11], and related variants [12–14]. Gauss-Hermite filter, introduced by [2, 4], makes use of the Gauss-Hermite quadrature rule and has the highest accuracy among all the above filters. However, it suffers from the curse of dimensionality since the number of cubature points increases exponentially with the dimension of the state.To avoid the curse of dimensionality, several sparse rules are developed. Smolyak rule [15] is one of the useful tools to generate a small number of cubature points for high dimensional integral. The computational cost does not increase exponentially by using this sparse grid method. Reference [16] compares the approximation accuracies of various sparse grid methods, including trapezoidal rule, Crenshaw Curtis rule, and Gauss Patterson rule. Reference [17] combines Smolyak’s rule and Gaussian Quadrature rule for high dimensional likelihood function in economic models. Jia et al. [18, 19] make a tracking comparison among Smolyak’s rule based Gauss-Hermite filter of different levels and the traditional cubature filters in the context of determination of the spacecraft attitude and the lower-earth orbit satellite orbit. Simulations prove the effectiveness of Smolyak’s rule based Gauss-Hermite filter. Reference [20] proposes the multiple sparse grid Gauss-Hermite filter based on sparse grid Gauss-Hermite filter and state-space partitioning. It is claimed that the computational burden is further reduced with respect to the Gauss-Hermite filter and the sparse grid Gauss-Hermite filter.However, as Heiss and Winschel [17] mentioned, there exists a potential drawback that some of the cubature points have negative weights, which may result in instability for computation. So we come up with the idea to avoid this potential drawback. In this paper, a classification of full tensor-product Gauss-Hermite cubature points is obtained by using equivalence classification by position permutation and signum function. There exist plenty of cubature points with low weights, which can be deleted from the cubature scheme directly considering a relative-weight-ratio threshold. Comparing with the Smolyak-based rule, our construction of the new Gauss-Hermite rule is simple and practically meaningful. The corresponding weights of the cubature points are all positive; meanwhile, the full symmetry property still remains. There also exist some interesting relationships among the Gauss-Hermite filter, the 3rd embedding cubature filter, the sparse Gauss-Hermite filter, and the Unscented Kalman filter.The remainder of this paper is organized as follows: In Section2, a brief review of the nonlinear system and its Bayesian recursive estimation framework is presented. In Section 3, a cycle of general local Gaussian filter is presented, which offers six kinds of Gaussian weighted integrals. This leads to various specific local Gaussian filters by various choice of the cubature points set and the corresponding weights. In Section 4, a full tensor-product Gauss-Hermite integral cubature rule is presented. By introducing a simple and elegant operator, the full tensor-product based cubature points can be divided into different categories. A small set of positive weighted cubature points is generated by a threshold of relative-weight-ratio. Relationships between the new sparse Gauss-Hermite filter and the traditional cubature filters are analysed as well. In Section 5, a typical Bearing-Ranging tracking problem with a general 2D manoeuvring target motion is demonstrated to test the performance of our new sparse Gauss-Hermite filter. Some conclusions are given in Section 6.


## 2. Bayesian Recursive Estimation Algorithm

Consider the following discrete nonlinear system:(1)xk=fxk-1+wk-1,(2)zk=hxk+vk,where (1) and (2) are the motion model and the measurement model of the system, respectively; xk∈Rn is the state of the system and zk∈Rm is the measurement; wk∈N(0,Qk) is the process noise and vk∈N(0,Rk) is the measurement noise. xk is independent with wk and vk. By the discrete time Chapman-Kolmogoroff equation,(3)pxk∣Z1:k-1=∫Rnpxk∣xk-1pxk-1∣Z1:k-1dxk-1,where Z1:k-1={z1,…,zk-1}. By the Bayesian formula,(4)pxk∣Z1:k=pzk∣xkpxk∣Z1:k-1∫Rnpzk∣xkpxk∣Z1:k-1dxk,where Z1:k={z1,…,zk}. Equations (3) and (4) are the time update formula and the measurement update formula. Figure 1 shows the recursive Bayesian estimation algorithm. When pxk∣Z1:k is obtained, we can estimate the state x^k and its covariance matrix Pk by minimizing the mean square error. The results are presented as follows:(5)x^k=Exk∣Z1:k=∫Rnxkpxk∣Z1:kdxk,Pk=Exk-x^kxk-x^kT∣Z1:k=∫Rnxk-x^kxk-x^kTpxk∣Z1:kdxk.However, when pxk∣Z1:k becomes an arbitrary probability density distribution, it is difficult to calculate the high dimensional integrations in (5) directly. A suboptimal way is to place normal distributions on pxk-1∣Z1:k-1, pxk∣Z1:k-1, and Pxk∣Z1:k; thus the first-order moment and second-order moment can be used to describe the posterior probability. In other words, at time epoch k-1, the posterior density is(6)pxk-1∣Z1:k-1=Nxk:x^k-1,Pk-1.At time epoch k, the predicted probability density of the state is(7)pxk∣Z1:k-1=Nxk:x^k∣k-1,Pk∣k-1.Correspondingly, the posterior density of the state is(8)pxk∣Z1:k=Nxk:x^k,Pk.Figure 1
One cycle of Bayesian recursive estimation algorithm [18].


## 3. Local Gaussian Filter

Section2 points out that, under the Gaussian assumption, the first-order moment and second-order moment can be used to describe the property of the state. Hence, a suboptimal local Gaussian filter algorithm can replace the Bayesian recursive estimation (3) and (4). Details of general Gaussian filter can be found in [4].(1)
Time update.
(i)
Prediction of statex^k∣k-1:(9)x^k∣k-1=Efxk-1+wk∣Z1:k-1=∫Rnfxk-1Nxk-1:x^k-1,Pk-1dxk-1.(ii)
Prediction of covariance matrixPk∣k-1 given state x^k∣k-1:(10)Pk∣k-1=Exk-x^k∣k-1·xk-x^k∣k-1T∣Z1:k-1=Efxk-1+wk-1-x^k∣k-1·fxk-1+wk-1-x^k∣k-1T∣Z1:k-1=Qk+∫Rnfxk-1-x^k∣k-1·fxk-1-x^k∣k-1TNxk-1:x^k-1,Pk-1dxk-1.(2)
Measurement update.
(i)
Prediction of measurementz^k∣k-1:(11)z^k∣k-1=Ehxk+vk∣Z1:k-1=∫RmhxkNxk:x^k∣k-1,Pk∣k-1dxk.(ii)
Covariance matrixPk∣k-1zz given z^k∣k-1: (12)Pk∣k-1zz=Ezk-z^k∣k-1zk-z^k∣k-1T∣Z1:k-1=Evk+hxk-z^k∣k-1vk+hxk-z^k∣k-1T∣Z1:k-1=Rk+∫Rnhxk-z^k∣k-1hxk-z^k∣k-1T·Nxk:x^k∣k-1,Pk∣k-1dxk.(iii)
Covariance matrixPk∣k-1xz between state x^k∣k-1 and measurement z^k∣k-1:(13)Pk∣k-1xz=Exk-x^k∣k-1zk-z^k∣k-1T∣Z1:k-1=Exk-x^k∣k-1hxk+vk-z^k∣k-1T∣Z1:k-1=∫Rnxk-x^k∣k-1hxk-z^k∣k-1TNxk:x^k∣k-1,Pk∣k-1dxk.(iv)
Gain of the filter:(14)Kk=Pk∣k-1xzPk∣k-1zz-1.(v)
State updatex^k and its covariance matrix Pk:(15)x^k=x^k∣k-1+Kkzk-z^k∣k-1,Pk=Pk∣k-1-KkPk∣k-1zzKkT.We can see that (9), (10), (14), and (15) are similar to those corresponding formulae of standard Kalman filter. The only difference between local Gaussian filter and standard Kalman filter is the calculation of x^k∣k-1, z^k∣k-1, Pk∣k-1, Pk∣k-1xz, and Pk∣k-1zz. For standard Kalman filter, there exist specific formulae, while for the local Gaussian filter, numerical approximation methods need to calculate the higher dimensional integrals. A general high dimensional integral with Gaussian type weighting function can be expressed as follows:(16)I=∫RnFt12πndet⁡Σ1/2×exp12t-x¯T×Σ-1t-x¯dt.Suppose that Σ=STS, where S is a lower triangular matrix. By a linear transform t=STS+x¯, we can get that(17)I=∫RnFSTS+x¯12πn/2e-1/2s2ds≈∑iNwiFSTξi+x¯.Various methods can approximate (17), such as the spherical radial integral rule [7, 21], Gauss-Hermite rule [2], central differential rule [4], Stochastic Integral Rule [22], Gaussian localized cumulative distribution rule [23], and Smolyak rule based sparse grid rule [18]. All of them use specific ways to generate {ξi}i=1N and their corresponding weights {wi}i=1N, where N is the number of cubature points set.Remark 1.
Generally, the choice of{ξi,wi}i=1N has no relationship with the estimation of state x^k and its covariance matrix Pk, but uniquely determined by the specific cubature integral rule. In fact, almost all cubature rules based {ξi,wi}i=1N can be computed and stored offline.


## 4. Sparse Gauss-Hermite Integral Rules

### 4.1. Symmetry Set Based on Full Position Permutation

Definition 2.
In order to utilize the symmetry property of the cubature points,[•]n is introduced to achieve a more concise denotation of the cubature points set. Denote • as a nonnegative seed; we have its corresponding seed vector v(•) as a n-dimensional vector (•,0)T, where n is the state dimension. Then we have [•]n={FP(v(•)),-FP(v(•))}, where FP(v(•)) means that the set of elements consists of the full permutation vectors of the seed vector and -FP(v(•)) is generated by randomly setting one or several nonzero values to their opposites in FP(v(•)).Let the state dimension ben and take [1,1]n as an example, where {1,1} is the seed. Since two values are given, then the n-dimensional seed vector is (1,1,0,…,0)T. It has total n2 full permutation, so [1,1]n has n2+21n2+22n2 elements. For simplicity, suppose n=3; then [1]3 has 31+31 elements, which is presented as the following set:(18)13=100,010,001-100,0-10,00-1.Remark 3.
With the help of full position permutation, it becomes much easier to describe the cubature points set{ξi,i=1,2,…,N}. Furthermore, it brings some interesting properties for the cubature points.


### 4.2. Full Tensor-Product Based Gauss-Hermte Integral Rule

The classical one-dimension Gauss-Hermite rule can be easily obtained according to [24]. Table 1 presents us with cubature points set {ξ} and corresponding weights {w} from orders 1 to 5. Here the order m also represents the number of cubature points in each dimension. We can see that the points and weights are both symmetric in R. For details about the construction of Gauss-Hermite quadrature points, one can refer to [25]. Furthermore, [26] provides a better and more stable way to construct cubature points and weights than moment matching method.Table 1
Cubature points and corresponding weights for one dimension(m=1,2,3,4,5).
m
{
ξ
}
{
w
}
1
{
0
}
{
1
}
2
{
-
1,1
}
1
2
,
1
2
3
{
-
3
,
0
,
+
3
}
1
6
,
2
3
,
1
6
4
{
-
2.334
,
-
0.742,0.742,2.334
}
{
0.046,0.454,0.454,0.046
}
5
{
-
2.857
,
-
1.356,0
,
1.356,2.857
}
{
0.0113,0.2221,0.5333,0.2221,0.0113
}Forn dimension integral problem, a full tensor-product based Gauss-Hermite rule can be obtained directly by the following formula [2]:(19)∫RnfxNx,0,Indx≈∑i=1mwin⋯∑i1=1mwi1fξi1,…,ξin=∑i=1mnwifξi.As mentioned above, n is the dimension of the state, m is the number of the cubature points for each dimension, and 1≤i≤mn as presented in (19). Each point of the cubature set is a n-dimensional vector; that is, ξi=[ξi1,…,ξin]T. The weight of point i can be calculated as wi=∏p=1n‍wip. For example, let m=3 and n=2; we can easily get the cubature points and the corresponding weights as follows:(20)ξ=0,0,0,-3,0,3,3,0,3,3,3,-3,-3,0,-3,3,-3,-3,w=49,19,19,19,136,136,19,136,136.Recall Definition 2; the cubature points in (20) can be categorized as three sets {ξ0,ξ1,ξ2}; that is,(21)ξ0=02,ξ1=32,ξ2=3,32,and the weights can be summarized as three types as well: (22)w0=49,w1=19,w2=136.The number of each category can be easily computed as(23)n0=1,n1=21+21=4,n2=22+2122+2222=4.The total number of cubature points N=n0+n1+n2=9=32. Equations (21), (22), and (23) show that such a classification is an equivalent notation of the full tensor-product rule. Now suppose that m=3 remains unchanged; we would like to consider the general case of n dimension. According to Definition 2, there exist n+1 categories. The corresponding {ξi,i=0,1,2,…,n}, {wi,i=0,1,2,…,n}, and the number of each category {ni,i=0,1,2,…,n} are shown as follows:(24)ξ0=0nξ1=3nξ2=3,3n⋮ξn=3,…,3nw0=23nw1=23n-116w2=23n-2162⋮wn=16nn0=1n1=n1+n1=2nn2=n220+21+22=2nn-1n3=n330+31+32+33=4nn-1n-23⋮nn=nnn0+n1+⋯+nn=2n.The total number of cubature pointsN=n0+n1+⋯+nn=3n. It also shows that such an expression is another form of the full tensor-product based rule. When the dimension of state n increases, the number of cubature points N increases exponentially with n. The curse of dimensionality limits the usage of Gauss-Hermite integral rule.Now, we are going to give the cubature points set and their corresponding weights withm=1,2,4,5. It can be easily conducted that when m=1, there is only one cubature point; that is, ξ=(0,0,…,0), and its weight is 1. When m=2, we can refer to Table 1 to find that there is only one nonnegative value (i.e., 1), so there exists only one category as well: (25)ξ0=1,1,…,1n;w0=12n,n0=2n.When m=2, it is interesting to find that the GHKF is similar to the 3rd embedding cubature Kalman filter [27]. We will elaborate their connections in the following Remark.Remark 4.
The main difference between GHKF and the 3rd embedding cubature Kalman filter is the central point(0,0,…,0)T and its weight w0 utilized. The classical standard 3rd cubature rule based on spherical radial rule [10, 28], whose cubature points set is [n], may exceed the range of state. Both the GHKF and the 3rd embedding cubature rule are stable for high dimensional integral rule. However, since the number of cubature points in 3rd embedding cubature Kalman filter is 2n+1, we still face the threat of dimension explosion when dealing with high dimensional system.Whenm=4, there still exists n+1 type of category which is analogous to the case m=3. Here we can simply regard p4,1 as 0 of case m=3 only to remember that p4,1 is positive. We have p4,1=0.7442 and  p4,2=1.3556 according to Table 1. The total number of cubature points becomes 4n. More specifically, the cubature points set {ξ}, the corresponding weights {w}, and the number {ni} for each category are shown as follows:(26)ξ0=p4,1,p4,1,…,p4,1nξ1=p4,1,…,p4,1,p4,2n⋮ξn=p4,2,p4,2,…,p4,2nw0=w4,1nw1=w4,1n-1w4,2w2=w4,1n-2w4,22⋮wn=w4,2nn0=2nn1=2nn1+n1=2n·2n2=2n·2nn-1⋮n4=2n,where w4,1=0.4541,w4,2=0.0459. It can be seen that the number of each cubature category has a factor 2n, which will result in the dimension explosion as well.WhenN=5, the total number of cubature points is N=5n. It can be classified as follows:(27)ξ0=0nξ1=p5,1nξ2=p5,2nξ3=p5,1,p5,1nξ4=p5,1,p5,2nξ5=p5,2,p5,2ξ6=p5,1,p5,1,p5,1nξ7=p5,1,p5,1,p5,2nξ8=p5,1,p5,2,p5,2nξ9=p5,2,p5,2,p5,2n⋮.The corresponding weights are presented as(28)w0=w5,0nw1=w5,0n-1w5,1w2=w5,0n-1w5,2w3=w5,0n-2w5,12w4=w5,0n-2w5,1w5,2w5=w5,0n-2w5,22w6=w5,0n-3w5,13w7=w5,0n-3w5,12w5,2w8=w5,0n-3w5,1w5,22w9=w5,0n-3w5,23⋮,where p5,1=1.3556,p5,2=2.857. w0=0.5333, w4,1=0.2221,w4,2=0.1123. We can calculate the number of each category to obtain the following equations:(29)n0=1n1=2nn2=2nn3=2nn-1n4=2nn-1n5=2nn-1n6=4nn-1n-23n7=4nn-1n-23n8=4nn-1n-23n9=4nn-1n-23⋮.Here we can see that, for m=2,4, the number of each category has a factor 2n, so those cubature points sets are not a good choice for the high dimensional system. However, the number of each category is a polynomial of the dimension of state when m=3,5. In the face of the dimension explosion, it is rather tempting that we can just choose several categories to approximate the high dimensional integration. We have introduced the relative-weight-ratio as a criterion to achieve this goal. Detailed analysis is presented in the following subsection.


### 4.3. Sparse Gauss-Hermite Rule Based on Relative-Weight-Ratio

Takem=3 as an example and recall (21) and (22) which present the cubature points set and corresponding weights set, respectively. The relative-weight-ratio between category ξ0=[0]n and category ξk=[3,…,3]n can be easily obtained:(30)Ratiok=w0wk=4kk≥0.Equation (30) shows that as the subindex k of the cubature points ξk set increases, the Ratiok decreases exponentially. It indicates that when considering a certain degree of accuracy, there is no need to use all Gauss-Hermite cubature points. When a threshold value K is settled, only the categories satisfying Ratiok≤K will be sufficient for the high dimensional integral. For example, set K=100; then we have(31)w0w3=64<100,w0w4=128>100.According to the criterion, the categories including [0]n, [3]n, [3,3]n, and [3,3,3]n are enough to approximate the integral. The total number of the cubature points is(32)N=1+2n+2nn-1+4nn-1n-13.The number of cubature points has a polynomial relationship with the dimension of state. If we set threshold value as K=50, the remaining categories will become [0]n, [3]n, and [3,3]n. The corresponding number of the cubature points is(33)N=1+2n+2nn-1.Remark 5.
By setting the ratio thresholdK, most low weighted cubature points will be eliminated. Such a strategy can efficiently eliminate the problem of dimension explosion when dealing with high dimensional system. Bear in mind that since lots of low weighted cubature points are eliminated, the sum of remaining points’ weights does not equal to 1. We can deal with this problem by renormalizing the remaining weights and meanwhile maintaining the relative-weight-ratio. Unlike the Unscented Transform rules [29] or the higher order spherical radial rule cubature rule [28], our transformation keeps all weights positive.Remark 6.
Whenm=3, considering [0]n and [3]n as selected cubature points, the total number is 2n+1 which equals the number of the UKF with parameter κ=3-n [29]. The difference lies in the determination of weights of the cubature points. In the context of sparse GHKF, the weight of [0]n is w0,GHKF=4/2n+4. However under the background of UKF, the weight becomes w0,UKF=1-n/3. Similarly, we have weights of [3]n as w1,GHKF=1/2n+4 and w1,UKF=1/6 for sparse GHKF and UKF, respectively. When dimension n>3, the weight of [0]n becomes negative in the UKF context, which is a potential unstable factor for integral approximation. The sparse Gauss-Hermite cubature points for m=5 can also be obtained by repeating the above discussed process with a preset relative-weight-ratio. Comparing [p5,1,p5,1]n with [p5,2]n, the weight of the former category is larger than the latter one.


## 4.1. Symmetry Set Based on Full Position Permutation

Definition 2.
In order to utilize the symmetry property of the cubature points,[•]n is introduced to achieve a more concise denotation of the cubature points set. Denote • as a nonnegative seed; we have its corresponding seed vector v(•) as a n-dimensional vector (•,0)T, where n is the state dimension. Then we have [•]n={FP(v(•)),-FP(v(•))}, where FP(v(•)) means that the set of elements consists of the full permutation vectors of the seed vector and -FP(v(•)) is generated by randomly setting one or several nonzero values to their opposites in FP(v(•)).Let the state dimension ben and take [1,1]n as an example, where {1,1} is the seed. Since two values are given, then the n-dimensional seed vector is (1,1,0,…,0)T. It has total n2 full permutation, so [1,1]n has n2+21n2+22n2 elements. For simplicity, suppose n=3; then [1]3 has 31+31 elements, which is presented as the following set:(18)13=100,010,001-100,0-10,00-1.Remark 3.
With the help of full position permutation, it becomes much easier to describe the cubature points set{ξi,i=1,2,…,N}. Furthermore, it brings some interesting properties for the cubature points.


## 4.2. Full Tensor-Product Based Gauss-Hermte Integral Rule

The classical one-dimension Gauss-Hermite rule can be easily obtained according to [24]. Table 1 presents us with cubature points set {ξ} and corresponding weights {w} from orders 1 to 5. Here the order m also represents the number of cubature points in each dimension. We can see that the points and weights are both symmetric in R. For details about the construction of Gauss-Hermite quadrature points, one can refer to [25]. Furthermore, [26] provides a better and more stable way to construct cubature points and weights than moment matching method.Table 1
Cubature points and corresponding weights for one dimension(m=1,2,3,4,5).
m
{
ξ
}
{
w
}
1
{
0
}
{
1
}
2
{
-
1,1
}
1
2
,
1
2
3
{
-
3
,
0
,
+
3
}
1
6
,
2
3
,
1
6
4
{
-
2.334
,
-
0.742,0.742,2.334
}
{
0.046,0.454,0.454,0.046
}
5
{
-
2.857
,
-
1.356,0
,
1.356,2.857
}
{
0.0113,0.2221,0.5333,0.2221,0.0113
}Forn dimension integral problem, a full tensor-product based Gauss-Hermite rule can be obtained directly by the following formula [2]:(19)∫RnfxNx,0,Indx≈∑i=1mwin⋯∑i1=1mwi1fξi1,…,ξin=∑i=1mnwifξi.As mentioned above, n is the dimension of the state, m is the number of the cubature points for each dimension, and 1≤i≤mn as presented in (19). Each point of the cubature set is a n-dimensional vector; that is, ξi=[ξi1,…,ξin]T. The weight of point i can be calculated as wi=∏p=1n‍wip. For example, let m=3 and n=2; we can easily get the cubature points and the corresponding weights as follows:(20)ξ=0,0,0,-3,0,3,3,0,3,3,3,-3,-3,0,-3,3,-3,-3,w=49,19,19,19,136,136,19,136,136.Recall Definition 2; the cubature points in (20) can be categorized as three sets {ξ0,ξ1,ξ2}; that is,(21)ξ0=02,ξ1=32,ξ2=3,32,and the weights can be summarized as three types as well: (22)w0=49,w1=19,w2=136.The number of each category can be easily computed as(23)n0=1,n1=21+21=4,n2=22+2122+2222=4.The total number of cubature points N=n0+n1+n2=9=32. Equations (21), (22), and (23) show that such a classification is an equivalent notation of the full tensor-product rule. Now suppose that m=3 remains unchanged; we would like to consider the general case of n dimension. According to Definition 2, there exist n+1 categories. The corresponding {ξi,i=0,1,2,…,n}, {wi,i=0,1,2,…,n}, and the number of each category {ni,i=0,1,2,…,n} are shown as follows:(24)ξ0=0nξ1=3nξ2=3,3n⋮ξn=3,…,3nw0=23nw1=23n-116w2=23n-2162⋮wn=16nn0=1n1=n1+n1=2nn2=n220+21+22=2nn-1n3=n330+31+32+33=4nn-1n-23⋮nn=nnn0+n1+⋯+nn=2n.The total number of cubature pointsN=n0+n1+⋯+nn=3n. It also shows that such an expression is another form of the full tensor-product based rule. When the dimension of state n increases, the number of cubature points N increases exponentially with n. The curse of dimensionality limits the usage of Gauss-Hermite integral rule.Now, we are going to give the cubature points set and their corresponding weights withm=1,2,4,5. It can be easily conducted that when m=1, there is only one cubature point; that is, ξ=(0,0,…,0), and its weight is 1. When m=2, we can refer to Table 1 to find that there is only one nonnegative value (i.e., 1), so there exists only one category as well: (25)ξ0=1,1,…,1n;w0=12n,n0=2n.When m=2, it is interesting to find that the GHKF is similar to the 3rd embedding cubature Kalman filter [27]. We will elaborate their connections in the following Remark.Remark 4.
The main difference between GHKF and the 3rd embedding cubature Kalman filter is the central point(0,0,…,0)T and its weight w0 utilized. The classical standard 3rd cubature rule based on spherical radial rule [10, 28], whose cubature points set is [n], may exceed the range of state. Both the GHKF and the 3rd embedding cubature rule are stable for high dimensional integral rule. However, since the number of cubature points in 3rd embedding cubature Kalman filter is 2n+1, we still face the threat of dimension explosion when dealing with high dimensional system.Whenm=4, there still exists n+1 type of category which is analogous to the case m=3. Here we can simply regard p4,1 as 0 of case m=3 only to remember that p4,1 is positive. We have p4,1=0.7442 and  p4,2=1.3556 according to Table 1. The total number of cubature points becomes 4n. More specifically, the cubature points set {ξ}, the corresponding weights {w}, and the number {ni} for each category are shown as follows:(26)ξ0=p4,1,p4,1,…,p4,1nξ1=p4,1,…,p4,1,p4,2n⋮ξn=p4,2,p4,2,…,p4,2nw0=w4,1nw1=w4,1n-1w4,2w2=w4,1n-2w4,22⋮wn=w4,2nn0=2nn1=2nn1+n1=2n·2n2=2n·2nn-1⋮n4=2n,where w4,1=0.4541,w4,2=0.0459. It can be seen that the number of each cubature category has a factor 2n, which will result in the dimension explosion as well.WhenN=5, the total number of cubature points is N=5n. It can be classified as follows:(27)ξ0=0nξ1=p5,1nξ2=p5,2nξ3=p5,1,p5,1nξ4=p5,1,p5,2nξ5=p5,2,p5,2ξ6=p5,1,p5,1,p5,1nξ7=p5,1,p5,1,p5,2nξ8=p5,1,p5,2,p5,2nξ9=p5,2,p5,2,p5,2n⋮.The corresponding weights are presented as(28)w0=w5,0nw1=w5,0n-1w5,1w2=w5,0n-1w5,2w3=w5,0n-2w5,12w4=w5,0n-2w5,1w5,2w5=w5,0n-2w5,22w6=w5,0n-3w5,13w7=w5,0n-3w5,12w5,2w8=w5,0n-3w5,1w5,22w9=w5,0n-3w5,23⋮,where p5,1=1.3556,p5,2=2.857. w0=0.5333, w4,1=0.2221,w4,2=0.1123. We can calculate the number of each category to obtain the following equations:(29)n0=1n1=2nn2=2nn3=2nn-1n4=2nn-1n5=2nn-1n6=4nn-1n-23n7=4nn-1n-23n8=4nn-1n-23n9=4nn-1n-23⋮.Here we can see that, for m=2,4, the number of each category has a factor 2n, so those cubature points sets are not a good choice for the high dimensional system. However, the number of each category is a polynomial of the dimension of state when m=3,5. In the face of the dimension explosion, it is rather tempting that we can just choose several categories to approximate the high dimensional integration. We have introduced the relative-weight-ratio as a criterion to achieve this goal. Detailed analysis is presented in the following subsection.


## 4.3. Sparse Gauss-Hermite Rule Based on Relative-Weight-Ratio

Takem=3 as an example and recall (21) and (22) which present the cubature points set and corresponding weights set, respectively. The relative-weight-ratio between category ξ0=[0]n and category ξk=[3,…,3]n can be easily obtained:(30)Ratiok=w0wk=4kk≥0.Equation (30) shows that as the subindex k of the cubature points ξk set increases, the Ratiok decreases exponentially. It indicates that when considering a certain degree of accuracy, there is no need to use all Gauss-Hermite cubature points. When a threshold value K is settled, only the categories satisfying Ratiok≤K will be sufficient for the high dimensional integral. For example, set K=100; then we have(31)w0w3=64<100,w0w4=128>100.According to the criterion, the categories including [0]n, [3]n, [3,3]n, and [3,3,3]n are enough to approximate the integral. The total number of the cubature points is(32)N=1+2n+2nn-1+4nn-1n-13.The number of cubature points has a polynomial relationship with the dimension of state. If we set threshold value as K=50, the remaining categories will become [0]n, [3]n, and [3,3]n. The corresponding number of the cubature points is(33)N=1+2n+2nn-1.Remark 5.
By setting the ratio thresholdK, most low weighted cubature points will be eliminated. Such a strategy can efficiently eliminate the problem of dimension explosion when dealing with high dimensional system. Bear in mind that since lots of low weighted cubature points are eliminated, the sum of remaining points’ weights does not equal to 1. We can deal with this problem by renormalizing the remaining weights and meanwhile maintaining the relative-weight-ratio. Unlike the Unscented Transform rules [29] or the higher order spherical radial rule cubature rule [28], our transformation keeps all weights positive.Remark 6.
Whenm=3, considering [0]n and [3]n as selected cubature points, the total number is 2n+1 which equals the number of the UKF with parameter κ=3-n [29]. The difference lies in the determination of weights of the cubature points. In the context of sparse GHKF, the weight of [0]n is w0,GHKF=4/2n+4. However under the background of UKF, the weight becomes w0,UKF=1-n/3. Similarly, we have weights of [3]n as w1,GHKF=1/2n+4 and w1,UKF=1/6 for sparse GHKF and UKF, respectively. When dimension n>3, the weight of [0]n becomes negative in the UKF context, which is a potential unstable factor for integral approximation. The sparse Gauss-Hermite cubature points for m=5 can also be obtained by repeating the above discussed process with a preset relative-weight-ratio. Comparing [p5,1,p5,1]n with [p5,2]n, the weight of the former category is larger than the latter one.


## 5. Numerical Investigation

In this section, a Nearly Constant Turning (NTC) manoeuvring target model is settled to verify the tracking performance of our new sparse Gauss-Hermite Kalman filter. NCT model is commonly used in aircraft flying control. We consider the 2-dimensional movement. Suppose the unknown turning rate isω. Since ω is unknown, the NCT model is a time-varying nonlinear model [10]. The discrete state model can be settled as follows:(34)xk=1sinωk-1Δtωk-10cosωk-1Δt-1ωk-100cosωk-1Δtωk-10-sinωk-1Δt001-cosωk-1Δtωk-11sinωk-1Δtωk-100sinωk-1Δt0cosωk-1Δt000001xk-1+wk,where xk=xkx˙kyky˙kωkT. xk and x˙k represent the position and velocity of axis x, respectively. yk and y˙k represent the position and velocity of axis y, respectively. ΔT is the sampling time interval. ωk is the turning rate at time epoch k. The corresponding process noise wk≈N(0,Q), where Q can be rewritten as follows:(35)Q=qxΔT33qxΔT22000qxΔT22qxΔT00000qyΔT33qyΔT22000qyΔT22qyΔT00000qwΔT.q
x and  qy represent the Power Spectrum Density for positions x and y. qw is the Power Spectrum Density of the turning rate. Suppose that the sensor platform is settled at the origin points. Bearing-Ranging measurement is used for the target tracking. The measurement equation can be expressed as follows:(36)zk=x1,k2+x3,k2tan-1x1,kx3,k+vk,where vk is the measurement noise, which satisfies the Gaussian distribution with zero mean. The corresponding covariance matrix R can be expressed as follows:(37)R=diagσr2σθ2.The corresponding parameters are settled by Table 2 [10].Table 2
Parameter setting for the experiment.
ΔT
1s
w0
-
3
d
e
g
/
s
q
x, qy
0.1
m
2
/
s
3
qw
1.75
×
1
0
-
4
s
-
3
σr
10
m
σθ
10mradWe set the initial state of the targetx0 as(38)x0=1000m300ms-11000m0ms-1-3degs-1. Then the corresponding covariance matrix P0 of the initial state is(39)P0=diag100m210m2s-2100m210m2s-2100mrad2s-2.The estimate of initial state of the filter x^0 is randomly generated by the distribution N(x:x0,P0). Our new sparse Gauss-Hermite Kalman filter (m=3) with the cubature sets {[0],[3],[3,3]} is used in the experiment. In order to make a reliable comparison between different filters, we implement our experiment with the (5th) CKF [28], full tensor-product based Gauss-Hermite Kalman filter (m=3) [2], and Smolyak rule based sparse Gauss-Hermite Kalman filter [18] (m=3).The above experiment is repeated 100 times independently. Root mean square errorsRMSEpos, RMSEvel, and RMSEω are used to evaluate the tracking performance of the position, velocity, and the turning rate. The detailed formulations are shown as follows: (40)RMSEposk=1Nmc∑i=1Nmcx^1,ki-x1,ki2+x^3,ki-x3,ki2,RMSEvelk=1Nmc∑i=1Nmcx^2,ki-x2,k2+x^4,ki-x4,ki2,RMSEωk=1Nmc∑i=1Nmcx^5,ki-x5,ki2.Similarly, CRMSE are used to evaluate the average tracking performance:(41)CRMSEposk=1Nmc∑i=1Nmc1Nsim∑j=1Nsimx^1,ji-x1,ji2+x^3,ji-x3,ji2,CRMSEvelk=1Nmc∑i=1Nmc1Nsim∑j=1Nsimx^2,ji-x2,ji2+x^4,ji-x4,ji2,CRMSEωk=1Nmc∑i=1Nmc1Nsim∑j=1Nsimx^5,ji-x5,ji2,wherex^k,ji represents the kth estimate element of state at time epoch j for the ith experiment. xk,ji represents the kth theoretical element of state at time epoch j for the ith experiment. Nmc represents the Monte Carlo experiment times. Nsim represents the state number of each experiment. We set Nmc=100. Nsim=100 in this experiment. For a single simulation experiment, the tracking performance is shown in Figures 2, 3, and 4. For Nmc repeated experiments, the tracking performance is shown in Figures 5, 6, and 7. The CRMSEs of different filters are shown in Table 3.Table 3
CRMSE of different filters.
Filter
CKF (5th)
GHKF (full)
SGHKF (Smolyak)
SGHKF (novel)
2
n
2 points
3
n points
4
n
2
+
1 points
2
n
2
+
1 points
CRMSEpos(m)
11.4217
11.4215
15.8146
14.2034
CRMSEvel(m/s)
6.6320
6.6320
14.2974
10.4270
CRMSEω(rad)
0.0088
0.0088
0.0207
0.0145Figure 2
Estimation error of position (one simulation).Figure 3
Estimation error of velocity (one simulation).Figure 4
Estimation error of turn rate (one simulation).Figure 5
RMSEs of position of different filters.Figure 6
RMSEs of velocity of different filters.Figure 7
RMSEs of turn rate of different filters.From Table3, it can be seen that the full tensor-product based Gauss-Hermite filter and the (5th) CKF have almost exactly the same tracking performance, while the Smolyak-based sparse Gauss-Hermite filter has relatively lower tracking accuracy. The new sparse Gauss-Hermite filter has better tracking accuracy comparing with sparse Gauss-Hermite filter and has less tracking accuracy comparing with the full tensor-product based Gauss-Hermite filter and the (5th) CKF. Our new sparse Gauss-Hermite filter is a good choice to track this manoeuvring target problem with a relative sparse cubature points set where all corresponding weights are kept positive. Furthermore, our method burdens low computational costs and maintains stable numerical accuracy.


## 6. Conclusion and Discussion

Gauss-Hermite cubature rule is an effective way to approximate nonlinear Gaussian weighted type integral. However, the full tensor-product based Gauss-Hermite cubature rule may cause a dimension explosion and is rarely helpful for high dimensional case. Smolyak rule based sparse Gauss-Hermite cubature rule is a way to generate small cubature set by specific generation design, which has a potential drawback that some cubature points may have negative weights. This novel Gauss-Hermite rule is easily generated by using our notation in Definition2 to classify the full tensor-product based Gauss-Hermite cubature points. Based on a proper relative-ratio-weight criterion, a sparse cubature set is regenerated. Simulations show that our sparse Gauss-Hermite rules will have comparable accuracy with the full tensor-product based Gauss-Hermite filter and higher precision comparing with Smolyak’s rule based Gauss-Hermite filter. Furthermore, our method can lead to a substantially reduced number of cubature points and a more stable high dimensional integration. Our method can maintain accuracy of integration at the same time and could be a good choice for practice problems in engineering.



---

*Source: 2783781-2017-09-11.xml*