CC BY
9
УДК 621.396.969.11
A-priory Error of the Initial Conditions While Solving the Problem of the Space Vehicle Navigation Using Pulsar Signals
Aleksey S. Konakov Vladimir I. Tislenko* Vyacheslav V. Shavrin
Research Institute of Radiotechnic Systems Tomsk State University of Control Sytems and Radioelectronics
Lenina, 40, Tomsk, 634050
Russia
Received 21.01.2016, received in revised form 27.03.2016, accepted 10.05.2016 The problem of a space vehicles (SV) navigation while processing the moments of signals from several space sources of Roentgen illumination incoming taking into consideration random error of measurements is discussed in the paper. The device of processing realizes Kalman filter quasi-optimal sigma-point algorithm and forms at the output estimates of SV current coordinates and velocity. It is done estimates error analysis when a-priory input and measured date about the SV position and velocity differ. There is also analyzed intensity of measurement of incoming moments errors due to the observation properties compared with the initial values. Method of statistical tests allowed to estimate the influence of a-priory inaccuracy on the root-mean-square error of the SV coordinate and velocity estimation.
Keywords: autonomic navigation in space, pulsar, non-linear filtration, sigma-point Kalman filter, error of estimation, a-priory data, Monte-Carlo method. DOI: 10.17516/1997-1397-2016-9-3-310-319.
Introduction
Nowadays in order to solve the problem of space vehicles (SV) navigation in near-the-Earth space there are used either ground stations or global navigation satellite systems (GNSS) [1,2]. Some advantages have alternative navigation variants connected with application of inertial navigation systems (INS) or astronomic navigation systems (ANS) [3-5], using the maps of the star heaven. But sometimes they are also not attractive. INS accumulate errors, ANS have problems of optics durability under the condition of the open space. In order to get the acceptable navigation errors the initial signals of outside beacons or onboard sensors (in case of INS) should have stable parameters [1]. It looks as if the Nature "has foreseen" the problems the Man would have while traveling in space. There are in space the needed beacons, this astronomic sources (AS) of electromagnetic illumination with sufficient stability level. They are far away star groups or galaxies called quasars and also some star systems or separate stars-pulsars [6-10]. One of the first publications [9] where errors of SV position estimation using pulsars was discussed belong to the 80-th of the last century. During the following years the interest to the problem grew [10].
* wolar1491@yandex.ru © Siberian Federal University. All rights reserved
In [6,7] there is described the history and evolution of the space navigation using pulsars, there are given results of illumination parameter stability investigation. Signals are received in the Roentgen band of electromagnetic illumination where there is maximal spectra density of the power flow [8]. The most stable are the signals of AS with periods of illumination from parts to tens of milliseconds. Stability of these signals is compared to the stability of nuclear frequency standard [6,10]. The problem of a SV autonomic navigation system (ANS) has many aspects. The authors activity is a continuation [11] and is connected with investigation of navigation processor algorithm forming the navigation vector (NV) current estimations and including SV position vector r(t) = [x(t) y(t) z(t)]T and vector of its velocity v(t) = [vx(t) vy(t) vz(t) ]T relative the Earth within heliocentric coordinate system.
The theoretical basis of processing algorithms synthesis in modern navigation systems is the Markov theory of non-linear filtration (MTNF) [12,13]. Application of MTNF suggests application of two kinds models. The first of them describes dynamic properties of NV condition x(t) = [r(t) v(t)]T with the noise nx(t) that appears due to the necessity to take into consideration random field perturbation near the Earth. This model is described by six stochastic differential equations of the first order. The second model made by algebraic equations shows the relation of the observation vector z(t) = h[x(t), nz(t)], NV elements x(t) and errors nz(t) always existing in real signal at the processor input and appear due to the receiver noise.
It is known that MTNF is based on the Bayes theory of errors [12] where the initial conditions x(0) for the SV dynamics are considered random with given math expectation mx(0) and covariance matrix Vx(0). So exist ensemble of SV real tracks. Observations provide information about real tracks of SV movement. The quality of SV x(t) estimates x(t) depends on on the covariance matrix of filtration errors Vx(t) = M{[x — x(t)] • [x — x]T}, where M is the operator of math expectation. Diagonal elements [Vx(t)]j¿ define root-mean-square error (RMSE) of SV estimates when initial conditions x(0) and many perturbations nx(t) and nz(t) are taken into consideration.
Let's mention that realization MTNF processor suggests introduction into the filtration algorithm a-priory information about statistical vector x(0) properties and random perturbations
nx(t) and nz(t) intensity [12]. In the filtration algorithm there are introduced x(f )(0) and its (f)
covariance matrix Vx^) as well as the data about the perturbations nx (t) and nz (t) intensity. For the optimal in sense RMS Bayesian estimate its initial values should be x(f)(0) = mx(0) and Vxf0) = Vx(0). It is obvious, that a-priory data introduced into processing algorithm differ
from the real ones having place due to the nature. This define the properties of the condition
(f)
and observation en-semble. Elements of the matrix Vx^) depends on the way and technique the a-priory data about the suggested orbit are received on board SV. In particular, the navigation of SV to orbit altitude to (1 ^ 3) • 105 km based on GNSS solutions provided radiochannel energy challenges, provides a position error standard deviation is not more units km. Star astronaviga-tion systems have error also comparable with this value.
Misalignment of the initial condition is especially important in the transient filter regime. De-
(f)
crease of the matrix Vx^) elements compared to their real values in the matrix Vx(0) describing tracks dispersion at the initial moment drag out the transient process of NV estimates RMSE. The filter "optimistically" accepts the information about initial estimates of condition and not effectively takes into consideration the real observations z(t) providing information about real conditions x(0). Decrease of the noise introduce in the filter compared to its real level make the filter too "trustful" concerning the current observation (in reality they are nor very reliable). At the opposite condition the filter is rather slow to accept the current observations. As the result
in both cases errors of estimates in the transient regime grow and becomes possible the unstable regime of processing.
In the paper the authors describe how the difference between the a-priory data introduced for filtration and they real values influence RMSE estimates of SV position and velocity. Sigma-point Kalman filter algorithm is used in the process of analysis.
1. Mathematical model of SV movement in the near the Earth space
Here there given the equations of SV movement in the Solar system near the Earth space. It is considered [3] that gravitation of the Sun, the Moon and the Earth is stronger than of other space objects. On order to make the navigation solution with this model more accurate non-sphericity of the Earth gravitation field is taken into consideration. Tidal forces, atmosphere resistance the Sun illumination pressure considered small compared to the effects mentioned above. Their influence is taken into consideration by introduction into the equations right parts random disturbances nK(t) as the vector of the white Gaussian non-correlated noise [3]. Equations of condition for the vector x(t) describing SV dynamics within heliocentric inertial coordinate system looks like follows [3]:
X(t)
m
v(t)
f(r(t), v(t)) + nx(t),
(1)
where components fi(r(t), v(t)) of the vector functions f(r(t), v(t)) in the (1) right part answer the relations:
fi(-) = vx(t); f2(-)= vy (t); f3(-)= vz (t); h(-) = ^ fo[r(t)]+ Vs
h(-) = fe(-) =
Ve y
Ve Z
fo[r(t)] + Vs
fo [r(t)]+ Vs
\ rs - r\ 3 \rs\3
ys - y ys
\ rs - r\ 3 \Zs\3
Zs - z Zs
K - r\3 |rs|3
+ Vm
+ Vm
+ Vm
xm -x xm
\ rm - r\ 3 \rm \3
ym - y ym
\ rm - r\ 3 \rm\3
zm -z zm
(2)
\rm - r\3 \rm\3
where r = \Jx2 + y2 + z2 is the distance from the Earth centre to the SV; Ve, Vs, Vm are gravitation constants of the Earth, the Sun and the Moon; rs = [xs ys ZsY, rm = [xm ym zm]T are vectors of the Earth, the Sun and the Moon position. General in (2) function
*<•>=' + '2 )) ((-5 ( Z )2) + J3 ( ^ )1 (3 - 7 ( Z )2) Z-■
- (^)4 8 (3 - 42 (Z)2 + 63 (Z)4) .
where Re is the Earth radius; J2, J3, J4 are coefficients of the zonal harmonics of the Earth field of gravitation. The random initial conditions x(0) for the system (1) have Gaussian probability distribution with parameters M [x(0)] = m0 and diagonal covariance matrix
Vo = M {[x(0) - mo] ■ [x(0) - mo]T } ■
3
r
3
r
2. Mathematical model of the observed signals
Observations onboard SV are formed as described in [14]. Equation (1) that defines dynamic of the SV movement allows to calculate the suggested position of the object in space and time. For some particular moment there is calculated the moment tbc of the AS signal coming into the Solar system barycenter (Fig. 1). The same moment is used to find the time of the signal from AS tsv coming. It is considered that within the Solar system the unit vector nas describing the direction to AS is the same at any point of the near space. So in heliocentric coordinate system the difference At = (tbc — tsv) is defined by the projection of vector rsv on vector nas. The first approximation of the difference (relative correction trel taken into consideration) is [14]
1
1
1
At tbc tsv (nas , rsv ) + trel (nas , rsv ) + ^ j [(nas , rsv )
c c 2cL
2^s
|rsv I2 — 2(rbc, rsv ) + 2(nas, rbc)(nas, rsv )] + ' ln
(nas, rsv ) + |rs
(nas , rbc) + |rbc
+ 1
(3)
where (a, b) is the vectors scalar product; L is the distance between the Sun and the AS; c = 2.99792458 -108 (m/s) is the light speed; rbc is the vector defined by the barycenter position. The second component in the right part of (3) keeps Remer's delay which is the geometry delay taking into consideration Doppler effect and the early parallax. The last component keeps the temporal change of the electro-magnetic illumination at the presence of the Sun gravitation field (Shapiro effect). It was shown in [10] that the level of the relative correction of the delay time (equal to 482.13600481109), when SV receive ARGOS pulsar of Crab Nebula signals, is less than 10-7
Fig. 1. Propagation of the AS illumination in heliocentric coordinate system
The math model of the observed i-th AS signal is described by (3) if At(i) is written with rsv = r(t) + r_E(t) seen in Fig. 1. Dimension of the vector z(t) corresponds the number m of pulsars whose signals are used for calculation SV coordinates and velocity. The vector components should be calculated in the distance units. Supposing the random errors additive (which looks correct from the physics point of view) the vector of the discreet in time observations would have the general view
z(tn ) = c - [At(1) At(2) ... At(i) ... At(m)]T + c - nz (tn ) = h(x(t„),t„) + c - nz (tn), (4)
the additive error nz(tn) is the m-order Gaussian vector of discreet statistically independent and in time errors of measurement. They have zero average and the given value of root-mean-
square deviation az. The signal illuminated by pulsars is periodical, it means that the calculated difference between delays of the foretold time of impulse incoming into barycenter and the time registered on the SV board has ambiguity. The foretold SV position can help to escape this ambiguity. As the period of the main part of AS applicable for SV navigation is about several milliseconds the error of SV, used to escape the ambiguity, position should not be more than 300 km, which is possible. So in order to find NV estimate it is necessary at every step to predict the position according to (1). Further the suggested coordinates are made more accurate according to the measured time of delay.
3. Algorithm of SV navigation parameters filtration
Filtration problem of NV x(t) with dynamics according to (1) and observations (3) belongs to the class of non-linear problems MTNF because of equations (1), (3) non-linearity. Accumulation of signals within interval of processing AT up to several hundreds of seconds is used to find for sure the signal at the output of the photons detector and provide the needed accuracy of the time of incoming [6,7]. The concrete interval depends on AS illumination spectrum intensity, intensity of the outer phone illumination, square of the photons detector, loses during the process of accumulation and intensity of the detector own noise [7,15]. In particular, RMSE of the pulsar signal temporary position fixation recalculated into the distance, when the time of processing is about 100 sec and the photons detector square 1 m2, according to the data [7, Fig. 8] is about 1 km. So the level of noise in the channels of observation is high and the influence of non-linearity, as it is known, grows [12]. In order to solve the given non-linear problem quasi-optimal extended Kalman filter (EKF) algorithm [3,4,6] and sigma-point Kalman filter algorithm [5,11,14] are widely used. Both algorithms allow to realize in navigation processor recursive regime of current estimates NV X(t) forming and are oriented on Gaussian approximation of probability density W [x(tn)/Z(tn)], where Z(tn) = [z(t1 ),■■■, z(tn)] is the vector of the whole complex of observations [4,16,17].
SPKF algorithm [5,16,17], differs from EKF algorithm by the fact that in order to calculate the needed estimates of the condition and observation extrapolation is realized without application expansion of f(x(t)) and h(x(t)) into the Taylor series and the following linearization. Algorithm SPKF utilize probability density approximation. In order to do that at any step there is formed of the sigma-points set [xa\tn); i = 1, ■ ■ ■, N}, making simplex with the center at the point of the current estimate x(t). Sigma-points dissipation depends of the current calculated covariation of errors. To do so there is used the procedure of covariation estimates matrix square root (Cholesky decomposition). As the result, extrapolative condition estimates, observation prognosis and covariation of the prognosis errors are calculated as weighted mean value on the sample set Xa which dimension is NX = 2nx +1, where nx = 6 is the dimension of NV. It allows to make linear approximation more correct and provide high accuracy at strong non-linearity and high level of disturbances. So it becomes not necessary to calculate Jakobi matrix at any step as in case of EKF.
In the paper in navigation processor there is used algorithm SPKF. The full description of algorithm is done in [5,11,14,17]. In [11] there is described stability of filtration at the same condition. The structure circuit of the system, at which output NV of SV is formed, is shown in Fig. 2.
Fig. 2. The structure circuit of the autonomous system measuring SV coordinates and velocity using pulsar illumination
4. Conditions of modeling and the results
From the physic's point of view it is obvious that signals from several AS (minimum three) are needed to measure SV coordinates and velocity. Calculations were performed for pulsars whose parameters according to [9] are given in the Tab. 1.
In the Tab. 1 can be seen that the main part of AS are in one plane orthogonal the Earth equator. It makes geometric factor (GDOP) [15] more important. So in the paper five signals from AS, mentioned in the table, were processed together while modeling the algorithm. The SV orbit is geostationary, its parameters are given in Tab. 2.
Calculation of real current NV is done by the system (1) integration in MATLAB using the Runge-Kutta method of the 5-th order with automatic choice of integration step and relative error 0.001. Elements of the movement model are found by the method of ephemerid from the theory DE-405. Ephemerides are available in binary format at the NASA sight (ftp://ssd.jpl.nasa.gov/pub/eph/planets/bsp/). Initialization of the real NV current values suggests introduction in the filter the initial estimates X(f)(0) and covariance matrix of its error (f)
x(o). Gaussian sequence nx(tj) corresponding the white noise in (1) is formed by the standard sensor with RMSE equal 10 m by coordinates and 0.01 m/s by velocity.
Analysis of NV probable estimates is done by the Monte-Carlo method averaging by N=75 runs. Statistical averaging is done by Gaussian values x(t0) and random disturbances nx(tj). Simulation is done within the time interval 104 s. Analysis of time difference (4) RMSE az influence on the NV estimates convergence character is done for az = 20 ¡s; 50 ¡s; 100 ¡s.
The quality of Estimates of position and velocity SV are used statistically average realization as root-mean-square error (RMSE)
RM S EPoS (tn)
1
N
—J2Mtn) -r(tn)] ; RMSEvel(tn)
N- 1
" i=1
1
N
ZTîYl [Vi (tn) - vi(tn)] '
N1
" i=1
Inaccuracy of a-priory data introduced in the filtration algorithm is described be three parameters: a = aXf)(0)/ax(0); 3 = oVf)(0)/av(0); 7 = a{J](0)/az(0). They define the mismatch of values introduced into the filter compared to the real ones. Values a, 7 change when the numerator valued are constant.
Table 1. Parameters of pulsars
Parameter PSR XTE PSR B 1617-155 PSR J0437-
B1937+21 J1751-305 J0218+4232 4715
Wavelength, m 467676 689523 695519 96833 1723810
Period, s 0.00156 0.00230 0.00232 0.00323 0.00575
Latitude, deg 57.51 7.80 184.56 Unknown -47.25
Longitude, deg -0.29 5,58 5.58 Unknown 69.31
Spectral inten- 4.99 • 10~5 6.65 • 10~5 6.65 • 10~5 1.62 • 10~5 0.93 • 10~5
sity, ph^sm2 •Hz
Signal shape pa- 7.41 1.18 22.55 76.2 50.8
rameter rs, s1/2
Signal shape pa- 7.41 1.18 22.55 76.1 31.4
rameter r&, s1/2
Table 2. Parameters of modeled orbit
x, m y, m z, m Vx, m Vy, m Vz, m
26214220335.009 33024715796.160 0 2408.201 -1911.571 0
Fig. 3 shows dependence of the error within the interval 104 s for position and velocity (one realization) and RMSE for position and velocity when RMS of the measurement error is
[ (f)]
az = 20 ¡s. Elements of the diagonal matrices VXf0)) = [V^o)]«. It means that the filter
L ()Jii
introduced data coinside coincide with the real values describing real SV tracks. Dispersions of the initial estimates are considered (9 km)2 by three coordinates and (100 m/s)2 by velocity vector components. Results of RMS calculation for the moment t = 5 • 103 s are shown in Tab. 3.
Table 3. Variation of RMSE of the SV position and velocity estimates depending on az when t = 5 • 103 s
RMSE of the difference of the incoming moments into barycenter and on SV board, az ¡s 20 50 100
RMSE of the SV velocity estimate, m/s 1.70 3.12 4.95
RMSE of the SV position estimate, m 2615.70 4297.71 9902.6
Fig. 4 shows the results of filtration algorithm susceptibility to the a-priory inaccuracy due to the matrices V(f)x(o) and Vx(o) mismatch. Calculation is done for af) = az = 20 ¡s, results are given for t = 5 • 103 s.
Tab. 4 provide information about RMS inaccuracy of the measured moments of incoming influence on the SV position and velocity RMSE. A-priory RMS inaccuracy is described by 7 = aZf Vz.
Calculation is done for a = /3 = 1 and, for the values increased with comparison with the previous values [vX^o)]ii = [Vx(0)](25km)2 by coordinates, (1 km/s)2 by velocity.
Fig. 3. Error of the SV position and velocity estimates depending on az = 20 ¡s xlOJ
11.5 10.5
9.5 8.5 1 7.5 6.5 5.5 4.5
----a=l
...... a=0J3
•t .............
J-.....
4 6
13
(a)
10
(b)
Fig. 4. (a) Variation of SV position RMSE at the moment t = 5 • 103 s from a-priory inaccuracy of initial errors estimates covariation; (b) Variation of SV velocity RMSE at the moment t = 5 • 103 s from a-priory inaccuracy of initial errors estimates covariation
Table 4. RMSE of SV position and velocity at the moment az when t = 5 • 103 s from a-priory inaccuracy of RMSE of time of incoming measurements, when a = 3 =1
7 = ]/°z 0.01 0.1 7 = 1; a{J) = 50 lis 10 100
RMSE of the SV velocity estimate, m/s 1113.7 907.3 860.4 973.7 1042.5
RMSE of the SV position estimate, m 981.4 697.4 482.9 791.6 1093.7
Conclusion
In the paper is performed analysis of the Kalman filter sigma-point algorithm sensitivity to the a-priory data about RMSE of initial estimates about SV position and velocity incorrectness. There are also investigated RMSE in the channel of observations which are introduced in filter while calculations of the SV navigation vector current estimates are performed. It is shown high sensitivity of the NV current estimates RMSE to the initial conditions RMSE that define ensemble of SV track along the given orbit. For instance, in case of parameters, that define SV NV filtration, deviation of a-priory RMS of the initial position and velocity estimates (parameter a) from the optimal value (a = 1) increases RMSE of SV position up to 30 percents and of SV velocity up to 20 percents.
The work is performed within the framework of a project of the Federal Special Purpose Program "Research and development in accordance with the priority development fields of science and industrial complex of Russia for 2014-2020", agreement no. 14-574-21.0101 (unique identifier RFMEFI57414X0101).
References
[1] N.V.Mikhaylov, 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, University of Colorado at Boulder, 2001.
[3] Z.Xiong, GEO satellite autonomous navigation using X-ray pulsar navigation and GNSS measurements, International Journal of Innov. Comp. Inform. and Control, 8(2012), no. 5, 2965-2977.
[4] X.Kai, W.Chunling, L.Liangdong, The use of X-ray pulsars for aiding navigation of satellites in constellations, Acta Astronautica, 64(2009) no. 4, 427-436.
[5] J.Ali, F.Jiancheng, SINS/ANS integration for augmented performance navigation solution using unscented Kalman filtering, Aerospace Science and Technology, 10(2006), 233-238.
[6] S.I.Sheikh, D.J.Pines, Recursive Estimation of Spacecraft Position and Velocity Using X-ray Pulsar Time of Arrival Measurements, Navigation, 53(2006), no. 3, 149-166.
[7] S.I.Sheikh, Spacecraft navigation using X-ray pulsars, Journal of Guidance, Control and Dynamics, 29(2006), no. 1, 49-63.
[8] J.Hanson, Noise analysis for X-ray navigation systems, Location and Navigation Symposium (IEEE/ION), 2008, 704-713.
[9] T.J.Chester, Navigation Using X-ray Pulsars, NASA Techn. Rep. N81-27129, 1981, 22-25.
[10] P.J.Buist, Overview of Pulsar Navigation: Past, Present and Future Trends, Navigation, 58(2011), no. 2, 153-164.
[11] A.S.Konakov, V.V Shavrin , V.I.Tislenko, Non-linear filtration algorithm synthesis in the problem of space devices navigation using pulsar and quasar signals, Trudy XII Vseross. konf o problemah otvetstvennosti t VSPU-2014, Moscow, 2014, 3646-3656 (in Russian).
[12] A.P.Sage, J.L.Melse, Estimation theory with application to communication and control, New York, McGraw-Hill, 1972.
[13] M.S.Iarlykov, Statistical theory of radio navigation, Radio i Sviaz, 1985 (in Russian).
[14] J.Liu, P.White, J.Ma, Doppler/XNAV-integrated navigation system using small-area X-ray sensor, IET Radar Sonar Navigation, 2011, 1010-1017.
[15] J.Wu, Z.Yang, N.Yang, The accuracy analysis of the spacecraft autonomous navigation system based on X-ray pulsars, Journal Systems and Control Engineering, 227(2010), 121-128.
[16] J.Julier, J.K.Uhlmann, A new approach for filtering nonlinear systems, Proceedings of the American Control Conference, 1995, 1628-1632.
[17] R. van der Merwe, E.Wan, Sigma point Kalman filters for integrated navigation, Proceedings of the 60th Annual Meeting of the Institute of Navigation, 2004.
Априорная погрешность начальных условий в задаче навигации космического аппарата по сигналам пУльсаров
Алексей С. Конаков Владимир И. Тисленко Вячеслав В. Шаврин
Радиотехнический факультет Томский государственный университет систем управления и радиоэлектроники
Ленина, 40, Томск, 634050 Russia
В статье рассмотрена задача навигации космического аппарата (КА) на основе обработки наблюдений моментов времени прихода сигналов от нескольких космических источников рентгеновского излучения с учетом наличия случайной ошибки измерений. Устройство обработки реализует квазиоптимальный сигма-точечный алгоритм фильтра Калмана и формирует на выходе оценки текущих координат и скорости КА. Выполнен анализ погрешности оценок при наличии расхождения вводимых в алгоритм априорных данных о начальных оценках положения и скорости КА, а также об интенсивности ошибок измерений моментов времени с их фактическими значениями, обусловленными свойствами наблюдений. Методом статистических испытаний определено влияние априорной неточности на величину среднеквадратических погрешностей оценок координат и скорости КА.
Ключевые слова: автономная навигация в космосе, пульсары, нелинейная фильтрация, сигма-точечный фильтр Калмана, погрешность оценок, априорные данные, метод Монте-Карло.
