РАКЕТНО-КОСМИЧЕСКОЕ ПРИБОРОСТРОЕНИЕ И ИНФОРМАЦИОННЫЕ СИСТЕМЫ
2017, том 4, выпуск 1, c. 3–7
КОСМИЧЕСКИЕ НАВИГАЦИОННЫЕ СИСТЕМЫ И ПРИБОРЫ.
РАДИОЛОКАЦИЯ И РАДИОНАВИГАЦИЯ
УДК 692.7.052
Сигма-точечныйалгоритм фильтра Калмана
в задаче автономнойнавигации космического аппарата
В. А.Филимонов1, В.И.Тисленко2, В.Ю.Лебедев3,А.П.Кравец4
1,4аспирант, 2д. т. н., проф., 3к. т. н.
Томскийуниверситет систем управления ирадиоэлектроники
e-mail: Svv281088@gmail.com, wolar1491@yandex.ru
Аннотация. В работе выполнен анализ среднеквадратичной погрешности (СКП) оценок вектора состояния, который содержит
текущие координаты, скорость космического аппарата (КА), а также смещения шкалы времени и частоты бортового
опорного генератора (ОГ). Оценки формируются в навигационном вычислителе, на вход которого поступают наблюдения в виде
псевдодальностей и псевдоскоростей по нескольким навигационным аппаратам ГНСС. Математическая модель состояния
учитывает наличие случайных возмущений по компонентам ускорения КА, обусловленных наличием вариаций орбиты КА,
каналы наблюдений содержат случайные ошибки в виде белого гауссовского шума. Модель вариаций частоты ОГ соответствует
свойствам рубидиевого стандарта. Показано, что в условиях движения КА на высокоэллиптической орбите при больших
величинах геометрического фактора GDOP и слабых сигналах алгоритм сигма-точечного фильтра Калмана обеспечивает СКП
местоположения менее 1 м, модуля скорости менее 0,01 м/с, смещения шкалы времени 1,9 нс и относительной нестабильности
частоты 3 · 10−13. Алгоритм сохраняет работоспособность при нарушении полноты созвездия НКА до нескольких минут.
Ключевые слова: автономная система навигации, космический аппарат, оценка координат, смещение бортовой шкалы времени,
нестабильность частоты, сигма-точечный фильтр Калмана
Sigma Point Algorithm of the Kalman Filter
in Spacecraft Autonomous Navigation
V.A.Filimonov1,V. I.Tislenko2,V.Yu.Lebedev3,A.P.Kravets4
1,4post-graduate student, 2doctor of engineering science, professor, 3candidate of engineering science
Tomsk State University of Control Systems and Radioelectronics
e-mail: Svv281088@gmail.com, wolar1491@yandex.ru
Abstract. This paper describes the analysis of a RMS error of the state vector estimations, which contain the current coordinates
and velocity of a spacecraft (SC).Moreover, it includes shift of the time scale and frequency of the onboard reference generator (RG).
Theestimates aretransferred to theoutputof the navigation computer,where they areformed. Theprocessing algorithmuses
observation data in the form of pseudoranges and pseudovelocities according to several GNSS SC. The mathematical state model
takes into account the presence of random fluctuations based on the SC acceleration components caused by SC orbit variations.
Observation channels contain random errors in the form of white Gaussian noise. The variation model of RG frequency corresponds
to the model of the rubidium frequency standard. The article shows that the sigma point Kalman filter provides a positioning RMS
error of less than 1 m, a velocity model error of less than 0.01 m/s, a time scale shift of less than 1.9 ns, and a frequency instability
of less than 3 · 10−13, in conditions when SC moves in a highly elliptical orbit (HEO) at high values of GDOP and weak signals.
The algorithm does not lose its efficiency even at incomplete constellations of GNSS SC up to several minutes.
Keywords: autonomous navigation system, spacecraft, coordinates estimation, onboard time scale bias, frequency instability, sigma
point Kalman filter
Стр.1
4
В.А.ФИЛИМОНОВ, В.И.ТИСЛЕНКО, В.Ю.ЛЕБЕДЕВ, А.П.КРАВЕЦ
Введение
Проблемы разработки и научно-технические
принципы проектирования систем автономной навигации
(САН) систематически изложены в [1–3].
В работе предполагается, что структура САН соответствует
схеме с двухэтапной процедурой решения
навигационной задачи [3]. На первом этапе
формируется вектор наблюдений z(tk), состоящий
из измерений псевдодальностей и псевдоскоростей,
содержащих возмущения. Вектор z(tk) связан
нелинейным безынерционным преобразованием
с информативным процессом x(t) соотношением
вида: z(tk)= h[x(tk),nz(tk)],где nz(tk) —вектор
гауссовских шумов, которые попадают в полосу каналов
слежения и обусловлены наличием собственного
шума на входе приемника.
На втором этапе обработки выполняется фильиз
шести компонент, определяющих текущие декартовы
координаты x(t), y(t), z(t) КА, составляющие
его скорости Vx(t), Vy(t), Vz(t), и дополнительно
двух компонент δT
трация наблюдений z(tk) с целью получения оценки
ог(t)=[δ(t) ˙δ(t)],определяющих
динамику вариаций шкалы времени δ(t)
бортового опорного генератора (ОГ).
Известно [4–6], что оптимальную (по квадx(tk)
при заданной последовательности наблюдений
Zk
ратичному критерию качества) текущую оценку
= Mxk/Zk
процесса x(t) и задание моделей в виде системы
дифференциальных уравнений (СДУ) первого порядка
позволяют получить рекурсивную процедуру
вычисления АПВ [4–6]. Численный алгоритм
фильтра частиц (particle filter) реализует вычисление
оптимальных текущих оценок состояния [4].
Практическое применение находят квазиоптиятностей
(АПВ) W xk/Zk
0. Свойство марковости
мальные алгоритмы, основанные на гауссовской аппроксимации
АПВ, что эквивалентно замене исходной
нелинейной задачи ее линейным аналогом
[4–6]. Далее, как правило, применяют известный
алгоритм расширенного фильтра Калмана
(extended Kalman filter — EKF) [4,6].
В данной работе используется сигма-точечный
алгоритм фильтра Калмана (sigma-point Kalman
тор условного математического ожидания
0 по апостериорной плотности веро0
= z0, z1, ... , zk определяет операxk
=
x(t) информативного процесса x(t), состоящего
filter — SPKF), позволяющий более корректно,
по сравнению с EKF, выполнить сведение нелинейной
задачи к линейному варианту [4, 8–10]. Применение
алгоритма SPKF для решения задачи навигации
показано в [12].
Постановка задачи и алгоритм
фильтрации
Вектор-столбец состояния x(t) определим
ввиде
xT (t)=
=[x(t) y(t) z(t) Vx(t) Vy(t) Vz(t) δt(t) ˙δt(t)] =
= xка(t) δог(t). (1)
Навигационный вычислитель выполняет функ
δ(t) и x8(t) ˙δ(t), последняя, связанная с частотой
ОГ, является непосредственно управляемой.
цию координатного и частотно-временного обеспечения
САН.
Отметим, что в (1) из двух компонент x7(t)
Математическая модель вектора z(tk) определена
заданием системы 2m (i = 1, . . . ,m) уравнений для
псевдодальностей и псевдоскоростей, которые имеют
вид
z1(i)
k = R(i)
k +n(i)
= (xk −X(i)
z2(i)
k = ˙R(i)
∆tk =
k )2 +(yk −Y (i)
k +n(i)
+c · δtk +c · n(i)
fk =
=(Vxk − ˙X(i)
где R(i)
Z(i)
+(Vzk − ˙Z(i)
векторов положения и скорости i-го НКА на моk
и ˙X(i)
k — текущая истинная дальность; X(i)
k , ˙Y (i)
k , ˙Z(i)
k — соответственно компоненты
мент времени tk в инерциальной геоцентрической
системе координат; ek =[ex ey ez]= r(t)/r(t) —
единичный вектор, определяющий направление визирования
по линии НКА–КА. Аддитивные возмущения
n(i)
zk = n(i)
∆tk n(i) T
fk
в виде случайных стационарных
гауссовских дискретных последовательРАКЕТНО-КОСМИЧЕСКОЕ
ПРИБОРОСТРОЕНИЕ И ИНФОРМАЦИОННЫЕ СИСТЕМЫ т. 4 вып. 1 2017
k )2 +(zk −Z(i)
∆tk;
k ) · ex +(Vyk − ˙Y (i)
k )20,5
+
k ) · ey+
k ) · ez +c · ˙δtk +c · n(i)
fk,(2)
k , Y (i)
k ,
Стр.2
СИГМА-ТОЧЕЧНЫЙ АЛГОРИТМ ФИЛЬТРА КАЛМАНА В ЗАДАЧЕ АВТОНОМНОЙ НАВИГАЦИИ 5
ностей имеют нулевые средние значения, некоррелированы
во времени и между собой в одном
ивразных каналах.
Система уравнений для x(t) определена заданием
шести нелинейных дифференциальных уравнений
орбитального движения КА для xка(t),которые
дополняются СДУ для компонент случайного
вектора δог(t). В уравнениях движения КА
учитывалась только неравномерность гравитационного
потенциала Земли. Уравнения орбитального
движения КА соответствуют приведенным в [1, 3]
и отличаются от них введением аддитивных некоррелированных
между собой БГШ по всем компонентам
ускорения, что позволяет учесть влияние
возмущающих орбиту КА факторов.
Оценка
туре управления ОГ на шаге фильтрации состояния
(1), на основе обработки наблюдений (2) с помощью
квазиоптимального фильтра, реализующего
алгоритм SPKF.
В дискретной форме с учетом управления ОГ
x(tk) формируется в замкнутом конпо
частоте разностные уравнения для переменных
состояния ОГ принимают следующий вид [11]:
x7(k)= x7(k −1)+ x8(k −1) · T +n7(k −1),
x8(k)= x8(k −1)+ δU(tk−1)+ n8(k −1),
(3)
где T — интервал временной дискретизации в вычислителе
при формировании оценок состояния
независимые во времени и коррелированные между
собой гауссовские последовательности. Ковариационная
матрица вектора nT
δU(tk−1) — приращение функции управления на
интервале T; n7(k − 1) и n8(k − 1) —дискретные
xk;
кретных наблюдений. Объектом управления является
ОГ.
Известно [7], что оптимальный алгоритм обнок
состояния.
Вданнойработеуправлениенеявляется опx7(k),
x8(k) как функции оптимальных оцеработки
(в линейном приближении) предполагает
формирование оптимальной линейной оценки
состояния (при заданном управлении) с последующим
вычислением оптимального управления
Uопт
тимальным и формируется в ПИ-регуляторе, который
содержит пропорциональную и интегральную
компоненты. Этот вариант требует меньших вычислительных
ресурсов. Значения δU(k) управляющей
последовательности в (3) формировались по следующему
правилу:
δU(tk)=
0для ∀ tk ∈ ((k −1)T0;kT0),
K7 · (1/T0)
tk+T0
-
tk
Ч
tk+T0
-
tk
x7(tk)+K8 · (1/T0)×
x8(tk) для tk = kT0,
(4)
время усреднения, в течение которого выполнялось
сглаживание оценок.
где K7, K8 — весовые коэффициенты, определяющие
вклад оценок
x7(k),
x8(k) вуправление; T0 —
Результаты моделирования
Вероятностные характеристики точности оце×10−28
(с−1) Sf ≈ 10−20 (рад2/Гц). Случайные начальные
условия x(0) имеют гауссовское распреответствует
виду, приведенному в [13]. При этом
величины спектральных плотностей возмущений
в модели ОГ соответствуют значениям для рубидиевого
стандарта частоты [11] и равны: Sg ≈ 7,9 ·×
деление вероятностей с параметрами M[x(0)] =
= m0 и диагональной дисперсионной матрицей
V0 =M[x(0)−m0] · [x(0)−m0]T.
Поставленная задача относится к классу нелинейных
задач динамической фильтрации и управления
[6, 7] с аддитивными белыми гауссовскими
шумами вмоделинепрерывных состоянийидисог(k)=
n7(k) n8(k) сонения
по ансамблю из 200 независимых реализаций
гауссовских шумов в модели состояний, наблюдений
и случайным гауссовским начальным
условиям x(0). Временной темп поступления данных
в канале наблюдений 1 мс. КА располагался на
высокоэллиптической орбите, геометрический фактор
(GDOP) при приеме сигналов 4 НКА равен 40,7
и в течение часа изменялся незначительно.
Интенсивности дискретных шумов в наблюнок
дениях σ2
∆t и σ2
∆t и σ2
и псевдоскорость равны соответственно 0,32 м2
f в пересчете на псевдодальность
и 0,009 м/с2. Отметим, что указанные значения
σ2
f достигаются в режиме слежения
РАКЕТНО-КОСМИЧЕСКОЕ ПРИБОРОСТРОЕНИЕ И ИНФОРМАЦИОННЫЕ СИСТЕМЫ т. 4 вып. 1 2017
x(k) определены путем статистического усред
Стр.3
6
В.А.ФИЛИМОНОВ, В.И.ТИСЛЕНКО, В.Ю.ЛЕБЕДЕВ, А.П.КРАВЕЦ
за задержкой и фазой навигационного сигнала
при отношении мощности навигационного сигнала
к спектральной плотности белого шума (C/N0)=
= 35 дБ-Гц в полосе 0,02 Гц. Интенсивности шумов
в модели состояния по компонентам ускорения
одинаковы и равны Dg = 10−5 м/с2.
Параметры, определяющие приращение δU(k)
управляющей последовательности: Tu = 1с; q =
= 103; T0 = 25 с; K7 = 0,01; K8 = 1. Таким образом,
временной темп обновления управления составлял
25 с.
Диагональные элементы V0
динатам 105 м, по компонентам скорости 103 м/с,
по смещению шкалы времени 10−4 сипо относительной
частоте 10−7. Начальные оценки:
щие СКО начальных оценок
x1(0)=35 061 км,
= 9711,4 км,
x6(0)= −1224,8 м/с,
x4(0)= 1,3 м/с,
x2(0)= 28 118 км,
На рис.1и2наинтервале наблюдения 1чпоx7(0)=
10−6 с,
x8(0)= 10−7.
казаны несколько реализаций модуля ошибок текущих
оценок положения КА и его скорости при
приеме сигналов 4 НКА.
Рис. 3. Зависимость СКП модуля положения КА
x5(0)= 928,4 м/с,
x3(0)=
x(0),равны:по коорii
,
определяюРис.
2. Реализация модуля ошибок текущих оценок скорости
КА
Рис. 1. Реализация модуля ошибок текущих оценок положения
КА
При этом в конце интервала наблюдения
среднеквадратические погрешности (СКП) по ансамблю
реализаций равны: СКП положения КА —
0,55 м, СКП скорости — 0,002 м/с, СКП оценки
смещения шкалы времени — 1,9 нс, СКП оценки
смещения относительной частоты — 3 · 10−13.
На рис.3и4приведены зависимостиСКП
модуля соответственно положения и скорости КА
на интервале наблюдения 7 · 103 с. Расчет выполнен
для случая пропадания сигналов одного
Рис. 4. Зависимость СКП модуля скорости КА
НКА на интервале времени от 2000 с до 2900 с
и от 5000 с до 5900 с. На рис. 5 и 6 показаны одиночные
реализации оценок [c ·
в замкнутом контуре управления.
x7(k)] и c ·
x8(k)
Заключение
Полученные в работе результаты позволяют
сделать следующие выводы.
1. Применение сигма-точечного алгоритма
фильтра Калмана в задаче координатного и часРАКЕТНО-КОСМИЧЕСКОЕ
ПРИБОРОСТРОЕНИЕ И ИНФОРМАЦИОННЫЕ СИСТЕМЫ т. 4 вып. 1 2017
Стр.4
СИГМА-ТОЧЕЧНЫЙ АЛГОРИТМ ФИЛЬТРА КАЛМАНА В ЗАДАЧЕ АВТОНОМНОЙ НАВИГАЦИИ 7
Список литературы
1. Михайлов Н.В. Автономная навигация космических
аппаратов при помощи спутниковых радионавигационных
систем. СПб.: Политехника, 2014.
362 с.
2. Moreau M.C. GPS receiver architecture for autonomous
navigation in high Earth orbits. 2001.
P. 207.
Рис. 5. Оценка [c ·
x7(k)] в замкнутом контуре управления
3.
ГЛОНАСС. Принципы построения и функционирования
/ Под ред. А. И.Перова, В.Н.Харисова.
4-е изд., перераб. и доп. М.: Радиотехника, 2010.
800 с.
4. Simandl M. Lecture notes on state estimation of nonlinear
non-Gaussian stochastic systems // Department
of Cybernetics, Faculty of Applied Sciences, University
of West Bohemia, Pilsen. 2006. 154 p.
5. Тихонов В.И., Харисов В.Н.Статистическийанализ
и синтез радиотехнических устройств и систем.
М.: Радио и связь, 1991. 608 с.
6. Сейдж Э., Мелс Дж. Теория оценивания и ее применение
в связи и управлении / Пер. с англ. под ред.
Б.Р.Левина. М.: Связь, 1976. 496 с.
Рис. 6. Оценка c ·
x8(k) в замкнутом контуре управления
тотно-временного
обеспечения САН при движении
КА на ВЭО с величиной GDOP 40,7 и (C/N0)=
= 35 дБ-Гц в полосе 0,02 Гц приводит к величинам
СКП оценок, равным: 0,55 м и 0,002 м/с, соответственно,
по положению и скорости КА; 1,9 нс
и3 · 10−13, соответственно, по смещению шкалы
времени и относительной частоте бортового ОГ.
2. Изменение количества НКА с 4 до 3 на интервале
до 15 мин в условиях рассмотренной задачи
приводит к расходимости оценок положения с увеличением
СКП положения КА до 1,5–2 раз с последующим
восстановлением устойчивой работы алгоритма
при появлении дополнительного НКА.
Работа выполнена в рамках проекта ФЦП
«Исследования и разработки по приоритетным
направлениям развития научно-технологического
комплекса России на 2014–2020 годы», соглашение
№14.574.21.0101 (уникальный идентификатор
RFMEFI57414X0101).
7. Медич Дж. Статистически оптимальные линейные
оценки и управление / Пер. с англ. под ред.
А. С.Шаталова. М.: Энергия, 1973. 440 с.
8. Julier S. J., Uhlmann J.K.ANewApproach for Filtering
Nonlinear Systems // In Proc. of AeroSence:
The 11th Intern. Symp. on Aerospace Defence Sensing,
Simulation and Controls, Orlando FL, USA 1997.
9. Julier S. J. The Scaled Unscented Transformation //
Proceedings of the American Control Conference,
American Automatic Control Council, Evanston, IL.
2002. May. Р. 1108–1114.
10. Van Der Merwe R., Wan E.Sigma-Point Kalman
Filters for Probabilistic Inference in Dynamic StateSpace
Models // In Proceedings of the Workshop
on Advances in Machine Learning. 2003.
11. Graas F., Craig S., Pelgrum W., Ugazio S.Laboratory
and Flight Test Analysis of Rubidium Frequency
Reference Performance // Navigation, 2013. V. 60,
№2. P. 123–131.
12. Filimonov V., Shavrin V., Tislenko V., Kravets A.,
Lebedev V., Shkolniy V. Coordinate and Time-Frequency
Support of a Spacecraft Flight by Means
of Autonomic Navigation Using Sigma-Point Kalman
Filter Algorithm // Журнал Сибирского федерального
университета. Серия Математика и Физика,
2015, т. 8, №4. С. 385–393.
13. Brown R.G. and Hwang P.Y.C.IntroductiontoRandom
Signals and Applied Kalman Filtering, 2rd Edition,
Wiley, 1992. P. 502.
РАКЕТНО-КОСМИЧЕСКОЕ ПРИБОРОСТРОЕНИЕ И ИНФОРМАЦИОННЫЕ СИСТЕМЫ т. 4 вып. 1 2017
Стр.5