CC BY
9
УДК 692.7.052
Coordinate and Time-frequency Support of a Spacecraft Flight by Means of Autonomic Navigation Using Sigma-point Kalman Filter Algorithm
Vladimir A. Filimonov* Vyacheslav V. Shavrin^ Vladimir I. Tislenko* Aleksey P. Kravets Vitaliy Yu. Lebedev Vadim N. Shkolniy
Research Institute of Radiotechnic Systems Tomsk State University of Control Sytems and Radioelectronics
Lenina, 40, Tomsk, 634050 Russia
Received 20.09.2015, received in revised form 26.10.2015, accepted 11.11.2015 There are investigated errors of a spacecraft current coordinate and velocity estimates, time-frequency scale bias of the board receiver reference oscillator. The receiver is a component of the autonomic navigation system using signals of global navigation systems GLONASS or GPS. Estimates are formed by the navigation calculator on the bases of pseudo range and pseudo velocity measurements taken in the regime of navigation signals delay and Doppler frequencies tracking under the condition of the receiver noise. The calculator utilizes sigma-point Kalman filter. Analysis of estimations accuracy is done by the Monte-Carlo method while spacecraft moves along a high elliptic orbit.
Keywords: autonomic navigation system, spacecraft, coordinate estimation, on board time scale bias, sigma-point Kalman filter, statistic simulation. DOI: 10.17516/1997-1397-2015-8-4-385-393
1. Introduction and preliminaries
There has been constantly growing interest to the space autonomic navigation systems (ANS). These systems utilize the signals of the global navigation space systems (GNSS). Problems and principles of the ANS design are properly described in [1, 2]. Among many other problems let's discuss some, not connected with the technology of the system production. One of them is the geometric factor reaching tens and hundreds for spacecrafts (SC) put on geostationary or high elliptical orbits (GEO or HEO); next is the low level of the ratio of the carrier power to the noise power per bandwidth unit (C/N0) up to 40 dB-Hz or less; the problem is also the deficit of the navigation spacecrafts number under observation [2]. Modern ANS receivers are designed to receive not only GLONASS signals but also GPS (Navstar), Galileo, Compass, which allows to solve some problems including the number of navigation spacecrafts within the interval of observation.
The authors of the paper suggest that ANS utilizes two-step procedure while solving the navigation problem [3]. At the first stage is formed the vector of observation z(tk), that consist
* filimonov90va@gmail.com t svv281088@ gmail.com ^ wolar1491@yandex.ru © Siberian Federal University. All rights reserved
of the results of pseudo range and pseudo velocity measurements having fluctuations. The signal z(tk) is formed at the outputs of the channels of several space vehicles signals parameters tracking. Vector z(tk) is connected by the non-linear non-inertial transformation with the information process x(tk) as follows: z(tk) = h[x(tk), nz(tk)], where nz(tk) is the white Gaussian noise vector. The noise is within the tracking channel band and is due to the own noise at the receiver input.
The subject of investigation is connected with the processing second step where filtration of measurements z(tk) is done in order to get the most accurate estimation x(tk) of the information process x(t). It consist of six components describing the SC current Cartesian coordinates, components of its velocity and two components describing the on-board reference oscillator (RO) time scale dynamic. The bases of the optimal (and quasi-optimal) algorithms of the useful information filtering synthesis is the modern Markov non-linear filtering theory (Bayesian filtering) [4,5,6].
It is known [4,6,7] that the optimal current estimation X(tk) in case of the given sequence of observations Zk = z0, z1,..., zk is defined as the conditional math expected value Xk = M[xk/Zk] defined by the posteriori probability density (PPD) W[xk/Zk]. Markov property of the process x(t) and models as first-order differential equations system allows to get the recursive procedure of PPD calculation [4,5,6]. Forming the strict optimal estimation under the condition of observations and state nonlinear models suggests calculation of PPD at every time step. It is realized by the algorithms of the particle filter [4,5]. But they are not useful in case of on-board ANS, because they need powerful calculation resource to integrate using Monte-Carlo methods. Application find quasi optimal algorithms based on Gaussian approximation of PPD, which is equivalent to the procedure of the initial nonlinear problem replace by its linear analog [6,7]. Further is usually used the extended Kalman filter (EKF) algorithm [7]. In this paper is used sigma-point Kalman filter (SPKF) algorithm which allows to make nonlinear problem as linear more correct than EKF [4,13]. The peculiarities are discussed further in the part 3.
2. Filtration problem
Let's define the state vector-column x(t) information message as follows:
(t) = [x(t) y(t) z(t) Vx(t) Vy(t) Vz(t) St(t) St(t)] . (1)
x
Navigation calculator plays the role of coordinate and time-frequency support of ANS. At some sequential discrete moments t1,t2,... ,tk it calculates estimates Xk, yk, Zk of the SC current coordinates xk, yk, zk, estimates of the vector Vk velocity components Vxk, Vyk, Vzk and
estimates Stk and Stk of process St(tk) = Stk and velocity of its drift Stk , respectively. Vector Sro(t) = [<5(t) 5(t)] describes the state of the on-board reference oscillator, while St(t) is the reference oscillator relative frequency drift, which means St = [Sf (t)/f0] [3,14,15].
Let's notice that the second of two components x7(t) = S(t) and x8(t) = S(t),connected with the reference oscillator frequency is directly controlled. Within the frames of the modern Markov non-linear filtering theory the solving of the filtration information x(t) problem using the signal that could be measured z(tk) but containing disturbances suggests that the observation statistical z(tk) and the information messages models x(t) are known.
2.1. Math model of observations z(tk)
Vector z(tk) math model is described by the system of 2m equations for the observed pseudoranges and pseudo-velocities. According to (1) and (2) for V i = 1,... ,m we'll get
(i) A;
=R(i + c ■ Stk + nV = [(xk - Xf )2 + (yk - Y?>)2 + (Zk - Z?)T5 +
(i)
(i) 2 Lfc
(i) 2
Jk
+ c ■ Stk + c ■ n
Ai) _ ¿>(i)
2k
(i) . Atk.
(i)
Rk + c ■ Stk + nfk + (Vzk - Zki))2 ■ ez]a5 + c ■ Stk + c ■ nfk ;
[(Vxk - Xki))2 ■ ex + (Vyk - Yki))2 ■ ey+
(2)
where R^ is the current real distance; X^jY^Z^ and X(\
Y (i)
Y ,
Z(l) are components of
position and velocity vectors of the i-th navigation space vehicles at the moment t within inertial geocentric coordinate system. e k = [ex ey ez] = r(t)/r(t) is the unity vector describing direction of the navigation space vehicles to the ANS receiver. Measurements of each pair of 2m vector z k elements at (2) contains additive disturbance nZk = [n^tk nf ]T as random stationary discrete sequences, which have zero averages, are not correlated in time and among themselves in the same and different channels. The noise covariance matrix in the channels of observation is block-diagonal with elements:
D
(i,i) 22
M[i
(i) lz k
(nZi)k )T ] =
2
t
0
D
(i,j) 22
M[
(i) zk
(nj k )T ]=0,
(3)
where M is operator of the mathematical expectation; a&t2 and af2 are dispersions of the Doppler bias of pseudo-delay and pseudo-frequency bias discreet measurements error at the output of the tracker systems measuring these parameters. In the matrix form equations (2), when (3) is taken into consideration, looks like follows:
Zk = h(xk) + c ■ nzk,
(4)
where h(xk) is the non-linear vector-function of the equations (2); right parts nzk is a block 2m vector of Gaussian sequences of noise with covariance matrix (3); c is speed of light.
2
a
2.2. Math model of the information message x(t)
The system of equations for x(t) is formed by six non-linear differential equations of the SC orbit motion for xsc(t). They are added by the system of two equations for the random vector Sro(t) components. In equations of the SC orbit motion we'll take into consideration the Earth gravitation field irregularity and neglect the gravitation of the Sun and the Moon. Equations of the orbit motion are given in [1,3] and looks like follows:
Xsc(t)= [rsc(t) Vsc(t)]T = fSc(rSc (t), Vsc(t)) + ng (t), (5)
wherefsc(rsc(t), Vsc(t)) is the given vector function; rsc(t) = y/(x2(t) + y2(t)+ z2(t)) is the distance from the Earth centre to the SC. Disturbances by acceleration components that were not taken into consideration are described as additive non-correlated white Gaussian noise by three components of the SC acceleration:
ng (t)= [0 0 0 ng x(t) ng y (t) ng Z (t)]T . (6)
In case of discrete Euler's approximation (5) the white noise (6) transforms into Gaussian vector sequence n*(kT) which components have dispersion Dg and are not correlated.
Estimation of state is formed within the closed loop of the dynamic object control at the step of state filtering (1) on the bases of observation results (4) with the help of quasi-optimal filter realizing algorithm SPKF [9].
In the discrete form when the reference oscillator frequency control is taken into consideration residual equations for the reference oscillator state variables look like follows [14]:
xr(k) = xr(k — 1) + x8(k — 1) • T + nr(k — 1);
(7)
x8(k) = x8(k — 1) + SU (tk-l) + n8(k — 1),
where T is the interval of time digitization in the calculator during the process of state estimates xk ; SU(tk~i) is the increase of the control function within the interval T; n7(k — 1) and n8(k — 1) are discrete independent in time and non-correlated Gaussian sequences. Covariance matrix of the vector n'T'o(k) = [n7(k) n8(k)] is described by the following relation [15]:
Qn = M[n„(k) • nTo(k)]
[Sf • T + Sg • T3/3] Sg • T2/2 Sg • T2/2 Sg • T
where Sf is the spectral density of the white noise of excitation by relative frequency with dimension (s) and Sg is the spectral density of the white noise of excitation by the relative frequency variation acceleration with dimension (s_1).
In matrix form the system of nonlinear difference equations for the vector (1) looks as follows:
x(k) = f [x(k - 1); SU(tk-i)] + nx(k - 1). (8)
Random initial conditions x(0) for (8) have Gaussian distribution of probabilities with parameters M[x(0)] = m0 and diagonal dispersion matrix V0 = M{[x(0) — m0] • [x(0) — m0]T}; nx(k) = [nT (k) nUk)]T.
3. Algorithm of filtering. Sigma-point algorithm of the Kalman filter
The problem under consideration belongs to the class of non-linear problems of dynamic filtering and control [7, 8] with additive white Gaussian noise of the model having constant states and discrete observations. The object of control is the reference oscillator.
In practical applications is widely used the approach based on non-linear functions linearization in (4) and (8) near some known point [6,7]. As such a point one can use, for instance, the current estimation of state or the point belonging to some reference (nominal) track. In the end the one should deal with the linear Gaussian problem with the square quality criteria. It was shown, that in this case the theorem of division is true [7, 8]. It says that the optimum calculation algorithm (with linear approximation) suggests forming the optimal linear state estimation (when the means of control are known) followed by optimal control U[x7(k),x8(k)] calculation as a function of optimal state estimations.
In this paper control is not optimal and is formed in the controller with proportional and integrating components. This variant needs less calculation resources. Values of the control sequence in (7) are formed according to the rule
'0, for ytk e ((k — 1)To; kTo);
SU (tk ) = { t +T t +T (9)
K7 • (1/To)J2tk x7(tk) + K8 • (1/To)J2k+To xs(tk), for tk = kTo
where K7, K8 are weighting coefficients defining the estimations x7(k), x8(k) role in the control; To is the time of observation during which the estimations are being smoothed. Algorithm of observations processing is realized according to the Kalman filtering. In the general form the recursive equation of estimation looks like follows [4, 6, 7]:
x(k) = x-(k) + K(k)[z(k) - Z-(k)], (10)
where x-(k) is the condition estimation extrapolated by one step and utilizing all observations ZQ- from t0 to tk-1 ; Z-(k) is extrapolated estimation of observation; K(k) is the matrix of the Kalman filter amplifying coefficients.
Numerous modifications of the quasi-optimal recursive filtering algorithms [4,5,9,10] differ by the ways how the estimations x-(k), Z-(k) , appropriating them covariance matrixes and matrix K(k) are calculated. The methods of local approximation, one of them is EKF, utilize the vector-function expansion into Taylor series keeping linear (Gaussian filters of the first order) or quadratic (filters of the second order) members. SPKF is also one of algorithms with local approximation. It is based on the unscented transformation [11-13]. High in general case (compared with EKF) quality of estimation when UKF (SPKF) is utilized is due to the absence of traditional (based the Taylor series) nonlinearity approximation [11,12]. It really performed more accurate approximation of the probability density for finite set of sigma-points xai; i = 1,..., La, weighted with wi. In particular the a priori probability density at the step (k — 1) looks as follows:
W[xk-i/ZQ-1] Wi • S(xk-1 — x*i,k-i), (11)
=1
where S(:) is delta-function.
In this case for the extrapolated estimation X-(k) of state x(k) taking (11) into consideration we'll get
, .„ ^ „ • t(x„
'J-œ i=i
rœ r<œ
k-1
x-(k) = ... f (xk-i)W (xk-i/ZQ-1)dxk — 1 Wi • f (xai,k-i)-
Opposed to many algorithms of "particle filtering" utilizing computational integration based on Monte-Carlo method [4,5] generation of sigma-points is determinate. The full description of the SPKF algorithm is given in [11,13].
4. Results of modeling
Probabilistic characteristics of the accuracy of estimates x(k) are found by the statistic averaging of 200 independent Gaussian noise realizations in (4) and (8) under random Gaussian initial conditions x(0) at (8). Rate of the data coming to the channel of observation is 1 ms. The SC is put at the high elliptic orbit, the geometric factor (GDOP) in case of four navigation satellites is 40.7 and only slightly changes within an hour.
Parameters t and aj at (3) recalculated into pseudo-range and pseudo-velocity are (0.32 m)2 and (0.009 m/s)2. Let's notice that the mentioned values of a\t and af received at the regime of the delay and phase of navigation signal monitoring when ratio of the navigation signal power and the white noise spectral density (C/N0) = 35 dB-Hz within the band of 0.02 Hz. The noise intensity by acceleration components in (5) is the same and equal yDg = 10-5 (m2/s2). Disturbances spectral densities of the reference oscillator model corresponds rubidium frequency standard [14] and are equal Sg « 7.9 • 10-28 (s-1), Sf « 10-20 (rad?/Hz). Parameters defining the control sequence SU(k) are: Tu = 1 s; q = 103; To = 25 s; K7 = 0.01; K8 = 1. So the rate of
the control renewal is 25 s. Diagonal elements (%/V0)jj, defining the initial estimates x(0) mean square deviation by coordinates are 105 m, by velocity 103 (m/s), by the time scale bias 10-4 c and by relative frequency 10-7 . Initial estimates are: aii(0) = 35061 km; x2(0) = 28118 km; x3(0) = 9711.4 km; x4(0) = 1.3 m/s; x5(0) = 928.4 m/s; x6(0) = -1224.8 m/s; x7(0) = 10-6 s; x8(0) = 10-7. In Fig. 1 and 2 there are shown realizations of the SC position and its velocity current estimates errors (tracking errors) when during an hour there are received signals from 4 navigation satellites.
Time, set. Time, sec
Fig. 1. The tracking error of SC position
Fig. 2. The tracking error of SC velocity
At the end of observation interval the averaged root mean square errors are 0.55 m by the SC position, 0.002 m/s by velocity, 1.9 ns by the time scale bias, 3 ■ 10-13 by relative frequency. In Fig. 3 and 4 there are shown dependences of root mean square errors of the SC position and velocity estimates over the interval of observation 7 ■ 10-3 s. Calculation is done on the condition that signal of one of the SC disappears within the intervals from 2000 s till 2900 s and from 5000 s till 5900 s. One can see that the signal reappearance provide stable estimates forming.
Fig. 3. Root mean square error of the SC Fig. 4. Root mean square error of the SC position tracking velocity tracking
In Fig. 5 and 6 there shown solitary time-domain realizations of two reference oscillator state variable changing within the closed loop of control: [c ■ x7(k)] is bias of the time scale in meters and [c ■ x8(k)] is frequency relative instability in m/s.
Fig. 5. Time scale bias
Fig. 6. Reference oscillator frequency relative bias
In Fig. 7 there is shown the signal of control SU(k) realization within the closed loop of control by the dimensionless variable x8(k). The value SU(k) = 0 corresponds to the control when the reference oscillator frequency is equal to the frequency of the navigation SC on-board reference oscillator.
Fig. 7. Variations of the signal SU(k) within the closed contour of control
5. Conclusion
The carried out investigation allows to conclude as follows.
1. Utilization of the Kalman filter sigma-point algorithm in the system of a SC autonomic navigation at the high elliptic orbit allows to get root mean square error of the position and velocity estimations of 0.55 m and 0.002 m/s; time scale bias and relative on-board reference oscillator frequency 1.9 ns and 3 ■ 10-13 correspondingly. This result is got when the regime of the navigation signal delay and phase are observed in case of GDOP equal 40.7 and C/N0 = 35 dB-Hz within the band 0.02 Hz.
2. If the signal of one of 4 navigation spacecrafts is not observed during the time up to 15 min the error grows 1.5-2 times. After reappearance of the disappeared signal the filter again works normally.
Supported by Russia Ministry of Education and Sciences, agreement 14-574-21.0101 (ID RFMEFI57414X0101).
References
[1] N.V.Michaylov, Autonomous navigation of space crafts using GNSS, St. Peterburg, Po-litechnika, 2014 (in Russian).
[2] M.C.Moreau, GPS receiver architecture for autonomous navigation in high Earth orbits, PhD thesis. Univer-sity of Colorado, USA, 2001.
[3] Ed. A.Perov, V.Kharisov, GLONASS. Principles of designing and functioning, Moscow, Radiotechnika, 2010 (in Russian).
[4] M.Simandl, Lecture notes on state estimation of nonlinear non-Gaussian stochastic systems, Pilsen, Univ. of West Bohemia, 2006.
[5] M.S.Arulampalam, S.Maskell, N.Gordon, T.Clapp, A tutorial on particle filters for online nonlinear / non-Gaussian Bayesian tracking, IEEE Trans. on Signal Processing, 50(2002), no. 2, 174-178.
[6] V.I.Tikhonov, V.N.Kharisov, Statistical analysis and synthesis of radio engineering devices and systems, Moscow, Radio i svyaz, 1991 (in Russian).
[7] A.P.Sage,J.L.Melse, Estimation theory with application to communication and control, N.Y., McGraw-Hill, 1972.
[8] J.S.Meditch, Stochastic Optimal Linear Estimation and Control, N.-Y., McGraw-Hill, 1969.
[9] V.P.Loginov, Approximate algorithms of non-linear filtering, Zarubezhnaya radioelektronika, 1(1975), no. 2, 28-48 (in Russian).
10] V.P.Loginov, N.P.Ustinov, Approximate algorithms of non-linear filtering, Zarubezhnaya radioelektronika, 2(1964), no. 3, 3-28 (in Russian).
11] S.J.Julier, J.K.Uhlmann, A New Approach for Filtering Nonlinear Systems , Proc. of Aero Science. The 11th Int. Symp. on Aerospace Defence Sensing, Simulation and Controls, 1997.
12] S.J.Julier, The Scaled Unscented Transformation, Proc. of the American Control Conference, American Automatic Control Council, 1108-1114 (2002).
13] R. Van Der Merwe, E.Wan, Sigma-Point Kalman Filters for Probabilistic Inference in Dynamic State-Space Models, N.-Y., Proc. of the Workshop on Advances in Machine Learning, 2003.
14] F.Graas, S.Craig, W.Pelgrum, S.Ugazio, Laboratory and Flight Test Analysis of Rubidium Frequency Reference Performance , Navigation, 60(2013), no. 2, 123-131.
15] R.G.Brown, P.Y.Hwang, Introduction to Random Signals and Applied Kalman Filtering , Wiley, 1992.
Координатное и частотно-временное обеспечение полета космического аппарата в режиме автономной навигации с использованием сигма-точечного алгоритма фильтра Калмана
Владимир А. Филимонов Вячеслав В. Шаврин Владимир И. Тисленко Алексей П. Кравец Виталий Ю. Лебедев Вадим Н. Школьный
В работе выполнено исследование среднеквадратичной погрешности оценок текущих координат и скорости космического аппарата (КА), смещения шкалы времени и частоты опорного генератора бортового приемника КА, который оснащен системой автономной навигации,, использующей сигналы глобальных навигационных систем ГЛОНАСС или GPS. Оценки формируются в навигационном вычислителе на основе обработки наблюдений псевдодальностей и псевдоскоростей, формируемых навигационным приемником в режиме слежения за временными задержками навигационных сигналов и их частотами Доплера при наличии шума приемника. Навигационный вычислитель реализует сигма-точечный алгоритм фильтра Калмана. Анализ точности оценок выполнен численным методом Монте-Карло при движении КА на высокоэллиптической орбите.
Ключевые слова: автономная система навигации,, космический аппарат, оценка координат, смещение бортовой шкалы времени, сигма-точечный фильтр Калмана, статистическое моделирование.
