Фильтр Калмана в задаче наведения Текст научной статьи по специальности «Математика»

УДК 517.983:519.2:629.76:681.5

Фильтр Калмана в задаче наведения

© Ю.В. Журавлев МГТУ им. Н.Э. Баумана, Москва, 105005, Россия

Изучается простейшая кинематическая модель системы самонаведения для встречно-догонных режимов. Маневр цели моделируется стационарно-коррелированным случайным процессом. Для оценки координат фазового состояния кинематического звена, а именно угловой скорости линии визирования и нормального ускорения цели, применяется фильтр Калмана. Система дифференциальных уравнений Риккати для ковариаций ошибок ввиду нестационарности и нелинейности в замкнутом виде не поддается исследованию, однако в трех специальных случаях найдено замкнутое решение ее графическим или аналитическим способом. Указывается на необходимость имитационного моделирования с целью установления рамок применимости полученных в работе результатов.

Ключевые слова: система самонаведения, фильтр Калмана, уравнение Риккати.

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

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

Для практики наиболее ценны выводы и рекомендации, получаемые из математического и имитационного моделирования. Модель системы может быть представлена в виде нелинейных нестационарных стохастических систем конечно-разностных и дифференциальных уравнений, сменяющих одна другую в случайные моменты времени. Вообще говоря, принципиальных преград на этом пути иссле-

дований нет. Однако анализ качества этих сложных систем требует разработки специальных методов расчета характеристик устойчивости и точности.

Одно из направлений — анализ асимптотических свойств процессов с учетом нелинейности и нестационарности. Будем рассматривать калмановскую фильтрацию текущих измерений. Фильтр Кал-мана является байесовым оценщиком состояния. Поэтому в качестве априорной плотности вероятностного распределения фазового состояния принимаем статистические параметры маневренных характеристик цели. Можно использовать целый ряд статистических моделей маневра цели (см., например, [2]):

а) абсолютно случайностная. Нормальная перегрузка цели пц = v — это белый шум известной интенсивности . При этом

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

б) стационарно-коррелированная. Нормальной перегрузкой цели пц служит стационарный случайный процесс с заданной функцией

корреляции, например

_ы

^ц ( Т) = ^ц ^ ^

Соответствующий формирующий фильтр

• 1 2 пц =--пц + ап I-V,

ц Т ц ТЖ,

на вход которого подается белый шум V интенсивности Ж^,. Эту модель маневра используем в настоящей работе; в) предельная. При Т

Кц (Т) = а 2ц,

т. е. нормальной перегрузкой цели является гауссова случайная величина с заданным математическим ожиданием и дисперсией аПц .

Формирующий фильтр

пц = V.

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

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

противоположный через случайные промежутки времени Дt с показательной функцией плотности распределения:

f (*)

x=At

где ц — число переключений знака за единицу времени.

Кинематика наведения на цель. Используемые далее понятия и обозначения кинематической схемы наведения в земной вертикальной плоскости показаны на рис. 1 и соответствуют работе [1, гл. 8].

JIA -1—1-1--

Рис. 1. Кинематическая схема наведения ЛА на цель

Пусть L — расстояние между летательным аппаратом (ЛА) и целью (Ц); V и Уц — линейная скорость ЛА и цели соответственно; s — угол места цели, отсчитываемый от фиксированного инерциального направления Xg до линии визирования ЛА — Ц; п и пц — угол упреждения и угол убегания соответственно, отсчитываемые от направлений векторов скоростей ЛА и цели до направления линии визирования; 0, 0ц — траекторные углы ЛА и цели соответственно, отсчитываемые от фиксированного инерциального направления Xg до направлений векторов скоростей; n-n — подразумеваемое направление нормали к линии визирования; g — ускорение свободного падения.

Скорость сближения L и мгновенная угловая скорость s вращения линии визирования относительно ЛА имеют вид

L = -У cos ^ + Уц cos ; (1)

s = (Vsinц-Уц sinцц)/L. (2)

Угол места цели

8 = Л + 0 = ЛЦ +©ц. (3)

Продифференцируем (3) и результат домножения (2) на L по t: Li + Ls = V sin ц - Уц sin цц + V cos - ©) - Уц cos цц (s- ©ц),

откуда с учетом (1) получим

Ls + 2Ls = Уц©ц cos цц - V© cos ц + F (t), (4)

где

F (t ) = V sin ц - sin цц. (5)

Здесь F — проекция ускорений на линию n-n.

Запишем первое уравнение кинематического звена, исходя из (4), в виде

Ls + 2Ls= gnil cos цц - gn cos ц + F (t), (6)

Vц©ц V © ПА

где пц = ———, n =- — нормальные перегрузки цели и ЛА соот-

g g

ветственно; пц cos цц, n cos ц — проекции нормальных перегрузок на

нормаль к линии визирования.

Теперь по t продифференцируем (1) и с учетом (2), (3) получим второе кинематическое уравнение 2

L - Ls = -V© sin ц + Vn©ц sin цц - V cos ц + Vn cos цц. (7)

Уравнения (6), (7) представляют полную математическую модель кинематических связей в процессе наведения. На ранних этапах проектирования эту модель упрощают. Так, если направления векторов скоростей ЛА и цели таковы, что sin цц « sin ц « 0, cos цц « cos ц « 1, т. е.

процесс наведения развивается в окрестностях встречно-догонных ракурсов с узким сектором вариаций направления линии визирования, то F (t) « 0, и первое уравнение кинематического звена (6) принимает вид

Ls + 2Ls = £Иц -gn. (8)

Далее воспользуемся моделью (8).

Постановка задачи. Динамическая система с двумерным вектором состояния x(t) и скалярным измерением y(t) только первой координаты вектора состояния рассматривается как стохастическая модель наведения ЛА на цель:

x (t ) = A(t ) x(t ) + d (t )u (t ) + w(t ), y (t ) = hт x(t ) + ). (9)

Здесь

x(t) = (xi(t),x2(t))т, d(t) = (-L-1(t),0)т, w(t) = (0,^(t))\ h = (1,0)т,

(a(t) b(fj\ 0 -a

A(t ) =

причем

a(t) = -2L (t ), Ь(1 ) = L-1 (t), L (t) = dL|dt,

а — независимые центрированные случайные процессы

гауссова белого шума:

M ^) = M ф) = M ^ ) = 0,

M^(t)£(т) = 2аа28^ -т), а >0, а >0, Mц^)т) = г5(t -т), г > 0.

Верхним индексом «т» обозначается операция матричного транспонирования.

Угловая скорость вращения линии визирования

x1 ^) = 8,

нормальное к линии визирования ускорение цели

^) = Пц.

Для V t > 0 априори задаются следующие параметры: Mx(0) — среднее значение начального состояния; L(t) — текущая дальность от ЛА до цели; L ^) — относительная скорость сближения; а2 — интенсивность маневрирования цели; а — частота маневров цели по нормали п-п к линии ЛА — Ц; г — интенсивность шумового фона на выходе измерителя угловой скорости линии визирования; u — управляющее нормальное ускорение ЛА.

При вышеуказанных допущениях модель (9) используется на начальных этапах проектирования систем самонаведения.

Для выработки текущего управления нормальной перегрузкой ЛА необходимо оценивать вектор состояния x(t) по полной информации Y = {y(^), 0^^ t}, что реализуется фильтром Калмана — Бьюси [3, 4]:

x = A(t)x + d(t)u(t) + K(t) (y(t) - hx), < K (t) = P(t)hr-1, (10)

P(t) = A(t) P(t) + P(t) Aт (t) - P(t )hhт P(t )r 1 + 0,

где

P(t) = M (x(t) - x (t))(x(t) - x (t ))т; Q = x(0) = Mx(0), P(0) = P0.

(0 0 ^

v0 2aa y

К достоинствам такой схемы фильтрации относится возможность получения оценки Х(г) состояния х(г) в режиме текущего времени, а также то, что эта оценка является эффективной в смысле минимума следа ковариационной матрицы Р(г) и несмещенной, т. е. МХ(г) = Мх (г).

Процесс фильтрации реализуется в разомкнутом офлайн-режиме либо в замкнутом онлайн-режиме интеграцией матричного дифференциального уравнения Риккати для ковариационной матрицы Р(г). В силу симметричности матрицы Р(г) матричное дифференциальное уравнение Риккати эквивалентно системе трех обыкновенных дифференциальных уравнений с известными переменными коэффициентами а(г) и Ь(г) и постоянными параметрами а, а, г:

ри = 2а(г)рп + 2Ь(г)р12 - г-1р2, 4 р12 = (а(г)-а)Р12 + Ъ(г)Р22 -г^РпР^ (11)

.р22 = -2«Р22 -гР\2 +2«°2.

Интегрирование системы (11), далее называемой дифференциальной системой Риккати, должно вестись численно. В следующем разделе займемся не численной, а качественной стороной данной проблемы.

Исследование системы Риккати. Целью исследования будет нахождение решений системы (11) при наложении принципиально упрощающих допущений. Рамки применимости методики необходимо увязывать с особенностями конкретной ситуации и анализировать

в процессе имитационного моделирования. Имитационному моделированию будет посвящена отдельная работа. Здесь будут изучены три специальных случая.

Случай 1 (графическое решение при «замороженных» параметрах). На достаточно коротком временном промежутке [0, 7] рассмотрим систему (11) как систему с «замороженными» коэффициентами. Соответственно, коэффициенты а(^) и Ь(¿) системы (11) теперь будем считать постоянными.

Одним из актуальнейших вопросов калмановской фильтрации в стационарных динамических системах является изучение асимптотических свойств ковариационной матрицы при I ^ ю. Легко проверяется, что система (10) «объект - измеритель» при постоянных А и ё управляема по н и наблюдаема по и потому в соответствии с теорией [3] у дифференциальной системы Риккати (11) существует единственное стационарное решение (511, 512, э22), формально удовлетворяющее следующей нелинейной алгебраической системе Риккати:

0 = 2аэ11 + 2Ьэ12 - г -152, <0 = (а-а)э12 + Ьэ22 -г~ls11s12, (12)

0 = -2аs22 - г-152 + 2аа2.

Утверждение. В системе (12) э\2 Ф 0.

Действительно, в случае эй = 0 из второго уравнения системы (12) вытекало бы 522 = 0, но тогда из ее третьего уравнения вытекало бы ао = 0, что противоречит исходным условиям, и потому Ф 0.

В силу приведенного утверждения можно переписать систему (12) в эквивалентном виде:

1 2 а

'12 = 2ЬТ'П - Г11,

Ь гЬа

э11 = г(а-а) - э12 +-

2а 5

(13)

12

1 52 , _2 2га

522 = „ 512 .

Отобразим зависимости системы (13) графиками на рис. 2. В прямоугольной декартовой системе с осью абсцисс 5ц и осью ординат первому уравнению системы (13) соответствует парабола 1 на рис. 2

с вершиной А (га, га) и точками 0(0, 0) и В(2га, 0) на оси абсцисс. Второму уравнению системы (13) соответствует гипербола с ветвями 2 и 3 и асимптотами 4 и 5 (см. рис. 2); асимптота 5 совпадает с осью абсцисс; асимптота 4 наклонная с отрицательным угловым коэффициентом -2аЬ; центр гиперболы в точке С(г(а - а), 0).

^12

2 1 в ы

5

\

Тх в 6

Н 3

Рис. 2. Графическое отображение системы (13)

Парабола 1 может пересекаться с ветвями гиперболы 2 и 3 максимум в четырех точках Т1,2,3,4. Абсцисса точки Т4 отрицательна, что не отвечает решению системы (13), так как дисперсия 511 неотрицательна. При достаточно большом значении Ь асимптота 4 будет иметь большую отрицательную крутизну -2аЬ и пересечет ось абсцисс в точке С левее О, тогда точек Т1,2 вовсе не будет. Из соображений непрерывности лишь точка Т3 будет соответствовать искомому единственному стационарному решению уравнений системы (13) с постоянными а и Ь. Координаты точки Т3 дают искомые значения 511 и 512.

Искомое значение 522 можно найти с помощью третьего уравнения системы (13). На рис. 2 парабола, соответствующая этому уравнению, изображена кривой 6 с вершиной Е(а , 0) на оси абсцисс 522 и точками 0(0, аV2га ) и Н(0, - а у!2га ) пересечения с осью ординат. Точку Т3 спроецируем параллельно оси абсцисс на кривую 6 и получим точку М(522, 512) с абсциссой, равной искомому значению 522 дисперсии погрешности оценки второй координаты состояния.

Отметим, что с ростом Ь ветви 2 и 3 гиперболы все плотнее прижимаются к асимптотам 4 и 5, причем асимптота 4 перемещается влево с одновременным разворотом по часовой стрелке, стремясь занять параллельное оси ординат положение, а центр гиперболы (точка С) может перейти на отрицательную часть оси абсцисс, при этом Т3 ^ В и М ^ ¥, т. е. 5ц ^ (-4гЬ /Ь), S22 ^ а2, 512 ^ 0.

Случай 2 (медленный перехват). Изучим случай достаточной

малости коэффициентов а(^), Ь(^), когда мала величина |Ь)|/Ь(г).

Отбросив в системе (11) члены с а(() и Ь((), придем к упрощенной системе дифференциальных уравнений Риккати:

5П - Г

512 - а512 г 511512,

&22 - -2ал22 - г 1л122 + 2аа2.

(14)

Систему (14) можно решить аналитически. Действительно, находим точное решение сначала из первого уравнения системы (14):

зи(*) -

1

[ 5П(0)]-1 + г-1г

г > о,

(15)

затем из второго уравнения системы:

512(г) -

^12 (0)

-а/

лп(0) [лп(0)]-1 + г-1/'

/ > 0,

(16)

и, наконец, из третьего уравнения системы (14) с учетом выражения (16):

5И(/)-а2 -е (а2-5.2(0))-

^12 (0)

г-1/е-а/

5ц(0) [5П(0)]-1 + г-1/

г >0.

(17)

В выражениях (15)-(17) 5ц (0) = р^ (0) при /,] = 1, 2.

Пределы ковариаций при / ^ да следующие: 511(/) ^ 0, 512(/) ^ 0, 522(г) ^ О2.

Случай 3 (равномерное сближение). Рассмотрим наведение в режиме Ь(/) = Ь0 — V/ с постоянной скоростью V > 0. Дополнительно пренебрежем величиной 1/£(г), что будет соответствовать начальному участку процесса наведения. Тогда при 0 < / < будем иметь а(/) = 2V/(L0 - V/) > 0 и Ь(г) = 0. Вместо системы (11) будет исследоваться система

511 - 2а(г)511 - г~51Р

512 - (а(/) -а)512 - Г

12'

522 - - 2а522 - г 15122 + 2аа2.

(18)

В системе (18) первое уравнение является уравнением Бернулли первого порядка, общее решение которого

511 -(2 а(/} )) С1 + -1 е '2а(/} ЛЛ ].

Тогда на промежутке 0 < / < общее решение для 511 принимает следующий вид:

3гV

5ц(/) -•

(( -V/)(1 + 3гV(( -V/)

или

(19)

5„(/) -

V А0)

3гV

1 А

511(0) 3^

V А))

эЛ"

(20)

Из выражения (19) имеем

^ -

5П(0)А40 3гЩ

3 '

(21)

Если считать функцию 511(/) известной из (19), то по отношению к 512(/) второе уравнение системы (18) представляет собой уравнение первого порядка с разделяющимися переменными и общим решением вида

512(/) - С2ехР 1 |[а(/) -а

5ц(/)

ё/ к

(22)

При а(/) = 2^0 - V/)

512(/) - С2

3гV (А0 - V/) 1 + 3гV (( - V/) С1

Отсюда

^2 -

12

(0)

Фп (0)'

(23)

(24)

По отношению к 522 третье уравнение системы (18) является линейным неоднородным уравнением с общим решением

522

(/) - е-2а/ VС3 +1 |(2ага2 -4(/))е2аё/\ (25)

или, после подстановки (23)

C - 3

Ci 1 + 3rV (L0 - Vt)3 С1

C2 1

s21(t) = a2 +C3 ---—. (26)

Очевидно, что

,22(0) = a2 + C3 -CC21 'С =^2 + C3 -,ц(0)CC2 C1 1 + 3гУЬ30С1 C 3rV

= a2 +C3-512(0)C 1

C1 3rVL0 (27)

Произвольные постоянные С1, С2, С3 находятся без труда из (21), (24) и (27). Итак, уравнения (20), (23) и (26) — это аналитические выражения компонент ковариационной матрицы для случая 3.

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

ЛИТЕРАТУРА

[1] Козлов В.И. Системы автоматического управления летательными аппаратами. Москва, Машиностроение, 1979, 216 с.

[2] Федосов Е.А., ред. Динамическое проектирование систем управления автоматических маневренныхЛА. Москва, Машиностроение, 1997, 336 с.

[3] Kalman R.E., Bucy R.S. New results in linear filtering and predictions theory. J. Basic Eng., vol. 83, Ser. D, I, 1961, pp. 95-108.

[4] Bucy R.S., Joseph P.D. Filtering for stochastic processes with applications to guidance. New York, Wiley-Interscience, 1968, 217 p.

Статья поступила в редакцию 03.10.2014

Ссылку на эту статью просим оформлять следующим образом: Журавлев Ю.В. Фильтр Калмана в задаче наведения. Инженерный журнал: наука и инновации, 2014, вып. 07. URL: http://engjournal.ru/catalog/mathmodel/it/1240.html

Журавлев Юрий Васильевич родился в 1947 г., окончил МАИ в 1974 г., МГУ им. М.В. Ломоносова в 1981 г. Старший преподаватель кафедры «Вычислительная математика и математическая физика» МГТУ им. Н.Э. Баумана. Автор более 20 научных публикаций в области моделирования, идентификации, управления динамическими системами. e-mail: zhurjurwas270747@yandex.ru

Kalman filter in the problem of homing

© Yu.V. Zhuravlev Bauman Moscow State Technical University, Moscow, 105005, Russia

We studied a simple kinematic model of guidance system for counter-overtaking modes. The maneuver target is modeled by stationary correlated random process. We use Kalman filter to estimate the coordinates of the phase state of the kinematic link, namely, the angular velocity of the line of sight and the normal acceleration of the target. The system of differential Riccati equations for the covariances of the errors due to the nonstationarity and nonlinearity in closed form is not amenable to study. However, in three special cases we found its closed solution in graphical or analytical forms. The need of simulation with the purpose of determining applicability of the obtained results is shown.

Keywords: homing system, the kalman filter, the riccati equation.

REFERENCES

[1] Kozlov V.I. Sistemy avtomaticheskogo upravleniya letatelnymi apparatami [Aircraft automatic control systems]. Moscow, Mashinostroenie Publ., 1979, 216 p.

[2] Fedosov E.A., red. Dinamicheskoe proektirovanie system upravleniya avto-maticheskikh manevrennykh LA [Dynamic design of automatic control systems of maneuverable aircraft]. Moscow, Mashinostroenie Publ., 1997, 336 p.

[3] Kalman R.E., Bucy R.S. New results in linear filtering and theory predictions. J. Basic Eng., vol. 83, Ser. D, I, 1961, pp. 95-108.

[4] Bucy R.S., Joseph P.D. Filtering for stochastic processes with applications to guidance, New York, Wiley-Interscience, 1968, 217 p.

Zhuravlev Yu.V. (b. 1947) graduated from the Moscow Aviation Institute in 1974 and Lomonosov Moscow State University in 1981. Senior lecturer of the Department of Computational Mathematics and Mathematical Physics at Bauman Moscow State Technical University. Author of more than 20 scientific publications in the field of modeling, identification, control of dynamic systems. e-mail: zhurjurwas270747@yandex.ru