Определение параметров движения космического аппарата с помощью упрощенного фильтра Калмана Текст научной статьи по специальности «Физика»

УЧЕНЫЕ ЗАПИСКИ Ц А Г И Том VII : 1976

№ 1

УДК 629.7.051.853

ОПРЕДЕЛЕНИЕ ПАРАМЕТРОВ ДВИЖЕНИЯ КОСМИЧЕСКОГО АППАРАТА С ПОМОЩЬЮ УПРОЩЕННОГО ФИЛЬТРА КАЛМАНА

В. М. Поединок

Рассматривается применение оптимального и упрощенного (суб-оптимального) вариантов фильтра Калмана для определения параметров движения космического летательного аппарата на заключительном участке входа в атмосферу Земли по данным инерциальных и радиолокационных (позиционных) измерений. Упрощенный фильтр позволяет значительно сократить объем вычислений, выполняемых бортовой вычислительной машиной. По результатам численного моделирования процесса фильтрации сравнивается точность оценки параметров движения с помощью оптимального и упрощенного алгоритмов.

На заключительном участке спуска космического летательного аппарата (КЛА) в атмосфере ошибки инерциальных измерений в основном из-за дрейфа гироскопов могут превышать допустимые пределы. Поэтому возникает потребность в получении дополнительной информации и совместной обработке всей имеющейся информации. Задача рационального объединения и обработки информации от различных измерителей может быть решена с помощью фильтра Калмана [1]. В качестве дополнительной может использоваться, в частности, информация о положении, поступающая от радиолокационной станции [2—4]. Применение фильтра Калмана требует обработки измерительной информации на ЭЦВМ. При этом для обеспечения большей автономности целесообразно обрабатывать информацию на бортовой вычислительной машине (БЦВМ). Поскольку как память, так и быстродействие БЦВМ ограничены, возникает необходимость в создании упрощенных алгоритмов обработки измерительной информации.

1. Уравнения фильтра Калмана в случае совместной обработки инерциальных и радиолокационных измерений. Изменение за малый промежуток времени Д^ = 4+1 — 4 параметров движения КЛА в

прямоугольной системе координат, связанной с местом посадки, можно описать уравнением

есть вектор фазовых координат КЛА, Т — знак транспонирования

Дт = £а+1 — -с; а1 — компоненты вектора ускорений (отличных от гравитационного), 1 = х, у, г. Ось х в выбранной системе координат направлена по меридиану на Север, ось у — по радиусу Земли вверх, ось г — на Восток.

Информация о негравитационных ускорениях поступает от трех акселерометров, ориентированных по трем взаимно перпендикулярным направлениям. Предполагается, что основным источником ошибок инерциальных измерений является „уход“ осей акселерометров из-за дрейфа гироскопов. При этом выходные сигналы акселерометров формируются согласно модели:

где а1 — измеренные значения ускорений, аг — истинные значения ускорений, ?г — углы „ухода“ осей акселерометров (г —лу, г), которые удовлетворяют следующим уравнениям:

где

2к+1 — ~Т ^к ^к>

^к = £ (¿к) = [хк, ук, гк, хк, ук, 2Й]

(1.1)

(1.2)

О 0 Ми о о

к

ф

0 10 О Д4 о

0 0 10 о ык 0 0 0 1 о о

0 0 0 0 1 о

0 0 0 0 0 1

(1.3)

есть матрица перехода,

(к + \

1к +1

к

к

£

*к + 1

1к+\

1к +1

*к

(1.4)

-ал, + ауТ,; ау = ау + аг - ах т2;

(1.5)

Лг О'г' Му Iл 7у>

здесь Qj, 2у — проекции угловой скорости вращения Земли на соответствующие оси, —скорости „ухода“ осей акселерометров, обусловленные дрейфом гироскопов, 1 — Х, у, Z.

Используя (1.1), можно получить следующее выражение для априорной оценки вектора состояния КЛА:

Zk+\ — -j- A k — Gk> (1-6)

где знаком „Д“ обозначены оценки соо тветствуюших величин штрихом обозначаются апостериорные оценки.

Л, Л,

Используемые в Ak значения оценок ускорений щ можно поЛ _ Л,

лучить из (1.5), полагая а^щ, ъ малыми и пренебрегая произ-

А' Л' 1 ■ • \ ведениями ц X; (l, J — X, у, 2):

А, ~ ~ л, ~ л.

ах — ах-\-аг іу- —

А, ~ ~ л, ~ А ,

ау — ау az~\х~\- о-х Тг!

л, ~ ~ л, _ л,

аг = аг + а„и~ -а* Ту.

л,

здесь "¡і — апостериорные оценки углов ухода.

Л,

Таким образом, для определения вектора Ак требуется опреЛ

делить оценки углов ухода Эти оценки получим, дополнив вектор фазовых координат КЛА углами ухода и применяя для оценки расширенного вектора фильтр Калмана. В дальнейшем для краткости указанный расширенный вектор назовем вектором состояния фильтра.

Будем считать, что радиолокационная станция смещена относительно начала координат по оси х на величину х° и измеряет одновременно наклонную дальность сі, угол азимута ф и угол возвышения X. Поскольку й, ф, X нелинейно связаны с фазовыми координатами, а в уравнениях фильтра Калмана используются измерения, связанные с оцениваемыми фазовыми координатами линейно, в качестве вектора измерений будем использовать линейную часть разложения в ряд Тейлора величин сі, ф, X в окрестности оцененных значений фазовых координат [2]. Таким образом, вектор состояния фильтра будет иметь вид:

Хт — [Д*, Ду, Дг, Ах, Ду, Дг, Ду,, Дуг],

где

л ...

Д/ = і — і, і = X, у, г, х, у, г, у, тг-

Учитывая это, представим вектор измерений в виде

¥л^НкХ„ + ^

где

У1 = \аь-йь, фА-ф„ хк-%]

есть вектор измерений

я*

Hd(tk) 0 0 0 0 0 0

H^{tk) 0 0 0 0 0 0

Hx(tk) оооооо

Ял (tk) =

НА*н) Я*(**) =

л* л

л,

хь

L ик

Л

Ук

¿к

Л

Л

гц

dk

л*

л*2

Л*2 , Л2 ■** + Zk

*2

+ •?

Л2

Л Л

-ук гк

Т-У х'!-' -2

У

Í*2 4-?2 ■*ь +гк

№

К.

^—-измеренная наклонная дальность; х\-=хк-\- х°\ ^ = [^а^к),Ц(Ь)> М^)]т ~ вектор ошибок измерений, который имеет нормальное распределение с М[£*]=0; М[%кЩ — 8*,г/?; М ского ожидания; 8*, г — символ Кронекера; /? рица ошибок радиолокационных измерений

символ математиче-корреляционная мат-

R

0 0

о

о о

где

Учитывая (1.1), для вектора состояния фильтра можно записать: **+1 = Ф, ** + ДВ*. (1.8)

ДВк = Вк-Вк, Вк = [А1 0 0 0];

Ф1==

ф ®63

Oh ^33

Ф имеет вид (1.3), Оез — нулевая матрица размерности 6X3, /33 — единичная матрица размерности 3X3.

Учитывая (1.5) и (1.7), имеем Вк — Вк— 8Вк, Вк — Вк-\-ЪВк. По-

А

лагая 7,- на промежутках времени Дtk постоянными, представим

А

ЬВь и ЬВк в виде:

А

Тогда

где

ЬВк — Ск ГА; ьВк Д Вк — Ск Д1\,

С1=[С1 Оза]; rT* = [T;(íft); b(tky, уг (tk)];

г* — [тлг(^); ту(^кУ’ тz{tk)\\ — г*

А,

г,,

Матрица Ck размерности 3X6 имеет следующие элементы:

с1\ №к) — cî2 (?к) — C3S (tк) — ck\ (¿к) ~ СЪ2 {tк) = свЗ (tk) — 0)

*к +1 _

^ls (tk) := ^21 (tk) ~ i U,Arzd'

c13 (tk) — c3i (tk) c23 itк) ~

c!2 (tk) ~ ~ c5i (tk

(k +1 _

^32 (**) - j âtA-df,

‘k

*к + 1 _ i

^*+1 _ ,

ct3 (tk)~ c6i(tk)== ~ j* o-yd.%

h

^A+l

сьз№ь) ~ съîi^»)“* J ax^vi Ai — tk^i

Учитывая, что априорная оценка ^*+1=0, а также (1.8) и полученное выше выражение для ДВА, уравнение для априорной корреляционной матрицы запишем в виде:

/>*+, = М [Хк+! А'!,! ] = Ф! р; Ф1 + Ф, Ок С1 + Ск 01Ф1 + Ск Рз С1, (1.9) где

P'k” ' P\(tk) P2(tk) ; Dk = Pktk)

_P?(t„) P'3(tk). .Рз (tk)

/>!(/*) = М[ АЛ"(/*)ДЛГТ (**)]; АХт=*[Ах, Ау, Аг, Ах, А у, Аг)\ Р'^к) = М[АХ{1,!)АГ(1к)У, РН<*) = ^[ДГ(/*)ДГт(/,)1. ..

С учетом (1.9) и принятых ранее обозначений, уравнения девятимерного фильтра Калмана для случая совместной обработки инер-циальных и радиолокационных измерений будут иметь вид:

Х'к=Рк Hl(HkPkHl + Rk)-'Yk-

Л, Л А,

Z-z = ¿* + Хк\

К = Рк - Рк Hl(Hk Рк Н1 + R„)~' Hk Рк,

Zk+\ = Zk + Bk-\ Ck Г* — Gk,

Pk+l = Ф1 P'k Ф1 + Ф1 Dk Cl + CkDl Ф1 + C* PiCl ;

(1.10)

где

zl = ¡X (tk)у (tk) z(tk)x (tn)y(tk) z (tk) ^ (tk) (tk) í2 (Í*)].

Начальные условия:

Л

Xa

= 0: Z,,

P(t0)-

2. Упрощенный (субоптимальный) фильтр. Как видно из уравнений (1.10), необходимость учета ошибок инерциальных измерений приводит к увеличению размерности фильтра и, следовательно, к увеличению объема вычислений. Использование в тех же условиях шестимерного фильтра, оценивающего только вектор фазовых координат КЛА, приводит к большим ошибкам оценки, которые (вследствие расходимости фильтрации Калмана [5]) могут даже возрастать по мере обработки новой информации. Для устранения расходимости в шестимерном фильтре возможно использование метода, учитывающего потерю информационной ценности прошлых измерений [6]. В этом случае априорную корреляционную матрицу ошибок оценки определяем по формуле

'*+i

s4>/>; фт,

(2.1)

где

Xk)(Xk-

Xk)T], XI = [Axk Дyk Д zk Axk Дyk Дzk\.

Pk = M[(Xk

Матрица Ф имеет вид (1.3), а вектор Zk — вид (1.2). Величина скаляра 5 подбирается эмпирически и зависит как от величины ошибок инерциальных измерений, так и от интервала обработки информации. При этом величина s не должна быть слишком большой, так как с увеличением s учитывается все меньшее число прошлых измерений.

Априорное значение вектора фазовых координат КЛА определяется по формуле:

2*+1=Ф Zk + Ak-Qk, (2.2)

где Ак имеет вид (1.4), но вместо истинных ускорений at используются измеренные значения ускорений а;.

Матрица Нк имеет вид:

Hd(tk) ООО нь = Щ (t k) ООО Н\ (tk) ООО

В целях дальнейшего сокращения объема вычислений в урав-

нениях фильтра используется оценки вида [7|:

корреляционная матрица ошибок

Рхх (t*) 0 0 Рхх (t„) 0 0

0 Pn^k) 0 0 Руу (¿к) 0

0 0 Pzz (tk) 0 0 P*z (tk)

Рхх (tk) 0 0 Pxx(tk) 0 0

0 Руу (¿k) 0 0 Руу (tk) 0

0 0 Ріг i.tk) 0 0 pa (tk)

(2.3)

Остальные уравнения шестимерного фильтра имеют тот же вид, что и в девятимерном фильтре.

Как видно, в матрице (2.3) учитываются лишь диагональные элементы и элементы рXX, руу, Ргг, СООТВеТСТВуЮЩИе В Среднем наибольшим по модулю коэффициентам корреляции по сравнению с другими недиагональными элементами. Точность результатов, полученных при применении матрицы типа (2.3), может быть оценена путем численного моделирования.

Уменьшение размерности фильтра, приближенный учет ошибок инерциальных измерений и применение упрощенной корреляционной матрицы (2.3) позволяют снизить количество операций сложения, выполняемых при одном шаге оценки, с 477 до 36, а операций умножения — с 570 до 51. Если предположить, что на выполнение операции умножения затрачивается в пять раз больше времени, чем на выполнение операции сложения [8J, то объем вычислений на БЦВМ при использовании упрощенного фильтра снижается более чем в 10 раз. .

3. Результаты численного моделирования процесса фильтрации. Численное моделирование по определению параметров движения KJIA с помощью оптимального девятимерного и упрощенного шестимерного фильтров производилось для участка траектории снижения с высоты примерно 30 км до высоты 6 км с постоянным углом крена 30°. Движение КЛА рассматривается в прямоугольной системе координат, начало которой расположено на поверхности Земли, ось х направлена по меридиану на Север, ось у — по радиусу Земли вверх. Началу рассматриваемого участка соответствует вектор фазовых координат КЛА:

Zl = [—38282 м, 33040 м, 97 739 м, -337 м/с, —22 м/с, 1265 м/с]; концу — вектор:

Zion = [1024 м, 6350 м, —109 м, —158 м/с, —56 м/с, 48 м/с].

Продолжительность полета на рассматриваемом участке составляет 256 с.

Численное моделирование производилось при следующих условиях.

1) Начальная корреляционная матрица ошибок оценки имеет диагональный вид, ее диагональные элементы соответствуют дисперсиям ошибок оценки начальных фазовых координат. Значения этих дисперсий равны:

<4(i0)=16 км2; Qy (t0) — 4 км2, О* (/0) = 36 км2; \

4 (*о) — °1 (¿о) — 0,324-10~3 км2/с2; (t0) = 0,9-10~5 км2/с2; (3.1;

°т-* (^о) — aiy(to) — °тг(^о) — 0»2• 10 2. I

Оценка вектора начального состояния имеет вид

л л л л Л Л Л

Zo = [х0, Уо> *0» *0. Уо, zo> 0, 0, 0],

где ненулевые компоненты вектора определяются по формуле:

А . л ...

*0 — h + 1 “ X, у, z, х, у, z

(¿0 — истинное значение компонента вектора начального состояния; Дг0 — ошибки в оценке компонента начального вектора—гауссовы случайные величины с нулевыми математическими ожиданиями

и с дисперсиями (3.1). Конкретные значения этих ошибок, как и ошибок радиолокационных измерений, моделируются с помощью стандартного датчика случайных чисел).

2) Дисперсии ошибок радиолокационных измерений [8]:

Srf = 0,196XlO~3 км2; 4 - о2х = 0,ЗХЮ-5 .

Смещение радиолокационной станции от начала координат принято равным дг° = 20 км.

3) Скорости уходов осей акселерометров из-за дрейфа гироскопов г; = 1 град/час. Измеряемые значения ускорений а, (г = х, у, г) моделируются с помощью (1.5).

4) Интервалы обработки информации Д^ = 4с, проекции угловой скорости вращения Земли Qx = 0,67ХЮ~5 l/c, = 0,95Х10~5 l/c.

5) Величина скаляра s в шестимерном упрощенном фильтре при рассматриваемых интервале обработки информации и максимальной величине скорости ухода $¡—1 град/час в результате проведенного численного моделирования выбрана равной 1,5.

Поскольку корреляционная матрица ошибок измерений R имеет диагональный вид, апостериорные характеристики оценок вычислялись последовательно для каждого компонента вектора радиолокационных измерений [2]. В этом случае апостериорные значения оценки и корреляционной матрицы, полученные после обработки одного компонента вектора измерений, используются как априорные при обработке другого компонента. Такой прием позволяет избежать операции обращения матрицы и, следовательно, сократить объем вычислений.

На фиг. 1 и 2 приводится сравнение среднеквадратических ошибок оценки положения и скорости для девятимерного (оптимального) и шестимерного (упрощенного) фильтров, полученных осреднением результатов моделирования десяти реализаций процесса фильтрации. Каждой реализации соответствует определенная совокупность случайных ошибок в определении начальных данных и случайных ошибок радиолокационных измерений.

Из сравнения ошибок видно, что они уменьшаются по мере обработки новой измерительной информации в случае девятимер-

Фиг. 1

м/с 3-мернь/й фильтр б-меаный фильтр

1

1

1 1

v / 10

\ \ \ л-число ммерении

1

А \

’ Л к

Л ■*Ч

N

V ч. N ' > Л <■'* -s

-i 1

О 10 20 30 МО SO 60 п

Фиг. 2

ного фильтра. В случае шестимерного фильтра такое уменьшение ошибок происходит лишь до некоторого момента, после которого они остаются в основном на одном уровне. Это объясняется тем, что для упрощенного учета ошибок инерциальных измерений в шестимерном фильтре, используется метод, учитывающий потерю информационной ценности прошлых измерений: при определении оценки в данный момент учитываются фактически не все предыдущие измерения, а лишь некоторое их количество.

При этом среднеквадратическая ошибка в определении положения КЛА в конце рассмотренного при численном моделировании участка траектории (на высоте 6 км) увеличивается с 20 до 30 м, а в определении скорости — от 0,2 до 2 м/с.

Автор выражает благодарность В. А. Ярошевскому за полезное обсуждение настоящей работы и ценные замечания.

ЛИТЕРАТУРА

1. Kalman R. Е. A new approach to linear filtering and prediction problems. Transaction of the ASME: Journal of basic engineering, vol. 82, I960.

2. Богуславский И. А. Методы навигации и управления по неполной статистической информации. М., „Машиностроение“, 1970.

3. Hollister W. М., Brayard М. С. Optimum mixing of inertial navigator and position fix data. Journal of Aircraft, vol. 8, N 9, 1971.

4. Carter J. P. Space shuttle terminal navigation with conventional navigation aids. AIAA paper N 72-832.

5. S chi ее F. п., St an dish C. J., Tod a N. F. Divergence in the Kalman filler. AIAA Journ., vol. 5, N 6, 1967.

6. Tarn T. J. Zaborszky J. A practical nondiverging filter.

AIAA Journ., vol. 8, N 6, 1970.

7. Merrick R. B. A simplified estimator for an aircraft landing display. Journal of Aircraft, vol. 8, N l, 1971.

_ 8. Singer R. A., Behnke K. W. Real-time tracking filter evaluation and selection for tactical applications. IEEE transactions on aerospace and electronic systems, vol. AES—7, N l, 1971.

Рукопись поступила 25jll 1974 г.