Определение состояния спускаемого аппарата по данным автономных измерений с использованием квазилинейного фильтра Калмана

Рассмотрим задачу навигации спускаемого аппарата (СА) с малым и средним аэродинамическим качеством по данным автономных измерений. Пусть информация, необходимая для ре­шения задачи навигации, поступает от датчиков перегрузок (аксе­лерометров), а для обработки информации используется бортовая цифровая вычислительная машина (БЦВМ). В качестве случайных возмущений, оказывающих влияние на процесс оценивания состоя­ния СА, рассматриваются начальные ошибки «входа» в атмосферу, разбросы аэродинамических характеристик, геометрических и кон­структивных параметров СА, порывы ветра, вариации плотности атмосферы, ошибки измерения перегрузок СА.

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

Особую сложность представляет проблема построения форми­рующих фильтров для атмосферных возмущений — ветра и вариа­ций плотности,— поскольку в настоящее время эти возмущения при­нято рассматривать как нестационарные по высоте полета случай­ные процессы.

Остановимся на построении формирующих фильтров для атмо­сферных возмущений подробнее. Ветер принято задавать в виде трех статистически независимых составляющих, одна из которых горизонтальна и лежит в плоскости Z = 0, где Z — боковое от­клонение положения центра масс СА в скоростной системе коорди­нат, другая вертикальна и также лежит в плоскости Z = 0, а третья перпендикулярна плоскости Z = 0. Каждая из состав­ляющих скорости ветра W рассматривается как нестационарный гауссовский процесс, заданный в виде отрезка канонического разложений , где h — высота полета; , i—1, 2, ...,

n - координатные функции; — случайные величины, такие, что , где — символ Кронекера. Приведенный отрезок ка­нонического разложения однозначно определяет корреляционную функцию составляющей скорости ветра:

Аналогично вариации плотности атмосферы рассматривают­ся как нестационарный по высоте случайный процесс и задаются в виде отрезка канонического разложения

— координатные функции. Это разложение определяет корреляционную функцию вариаций плотности:

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

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

где — аргументы корреляционной функции, а последовательности образуют систему линейно независимых функций. При этом порядок уравнения формирующего фильтра равен числу членов канонического разложения я, которое может быть достаточ­но велико (n= 12...14).

С целью сокращения порядка формирующего фильтра может быть поставлена задача построения формирующего фильтра поряд­ка не выше заданного. Математически задача формируется следую­щим образом. Для нестационарного гауссовского процесса y(t) с корреляционной функцией K( ) требуется подобрать линейный дифференциальный оператор заранее заданного порядка, такой, чтобы корреляционная функция K*( ) случайного процесса x(t), полученного воздействием белого шума на фор­мирующий фильтр, определяемый L(p, t), как можно меньше отли­чалась бы от K( ) в смысле некоторого критерия.

Считаем, что корреляционная функция K( ) «гладко стекает» с диагонали квадрата — отрезок не­зависимой переменной, на котором задан процесс y(t). Иными сло­вами, предполагается, что скачки производных корреляционной функции K(t, x) по ее аргументам при t = отсутствуют. Это озна­чает, что в правой части дифференциального уравнения формирую­щего фильтра будут отсутствовать производные от белого шума [34] и задача сведется к определению коэффициентов следующего дифференциального уравнения формирую­щего фильтра:

Этому уравнению соответствует система из л-уравнений 1-го поряд­ка для n-мерного вектора х, включающего переменную и ее (п—1) производную:

Дифференциальное уравнение для корреляционной матрицы век­тора х (2.70) здесь имеет вид

где — интенсивность шума , или, в скалярной форме

В (3.101) K_ij — элементы корреляционной матрицы К_х. Посколь­ку матрица К_х и производные ее элементов считаются известны­ми на основании исходных статистических данных, систему уравнений (3.101) можно рассматривать как алгебраическую относительно не­известных . Нетрудно видеть, что из п² уравнений сис­темы (3.101) только п уравнений содержат коэффициенты — всего n+1 коэффициент. Таким образом, систе­ма уравнений (3.101) является неопределенной относительно и .

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

В качестве критерия можно принять квадратичный: , причем сечения получить путем интегриро­вания уравнений

при начальных условиях

где —элементы корреляционной матрицы К_х в соответст­вующие моменты времени.

Заметим, что и, следовательно, критерий J не зависят b(t).

При заданном порядке формирующего фильтра количество на­чальных условий должно равняться порядку фильтра. Так, напри­мер, если ограничимся первым порядком, то обеспечим совпадение дисперсии случайного процесса на выходе формирующего фильтра с дисперсией заданного случайного процесса. При втором порядке формирующего фильтра, кроме совпадения дисперсий, будет совпа­дать и первая производная на диагонали квадрата заданной корреляционной функции с производной корре­ляционной функции на выходе формирующего фильтра. Если поря­док уравнения будет n, то будут совпадать дисперсия и (п—1) про­изводных сечений корреляционных функций и . Таким образом, определяем коэффициенты в результате ми­нимизации функционала . Минимизация осуществля­ется, например, методом наискорейшего спуска [28]. Коэффициент b(t) определяем из последнего уравнения системы (3.101).

Анализ корреляционных функций , по­лученных по экспериментальным данным, проведенный в [20], по­казывает, что разумным порядком формирующего фильтра для ветра и вариаций плотности атмосферы является второй, так что соответствующее возмущение характеризуется ее величиной и пер­вой производной по высоте полета: . Производные более высокого порядка вводить нецелесообразно, так как они не обладают достаточной достоверностью.

Таким образом, формирующий фильтр ветра определяется урав­нениями

Здесь — обратная высота; h₀ — начальная высота спус­ка; h — текущая высота полета.

Это означает, что уравнения формирующего фильтра (3.102) записаны по обратной высоте, поскольку движение СА происходит при убывании высоты, а поэтому корреляционная функция процесса на входе формирующего фильтра должна_ховпадать с исходной именно при убывании высоты h, т. е. при возрастании обратной вы­соты h_обр.

Аналогично записывается и определяется дифференциальное уравнение формирующего фильтра для Δр:

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

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

Поэтому возможно использование математических моделей ат­мосферных возмущений, дающих хорошее приближение лишь на малых отрезках . Такие модели, называемые в дальнейшем локальными, вообще говоря, различны для различных отрезков . Вид и способ задания локальных моделей определяется в каждом конкретном случае отдельно: это может быть отрезок канонического разложения малой размерности, формирующий фильтр первого порядка. Можно использовать на всем отрезке [0, Т] модель одного вида, периодически изменяя («настраивая») ее па­раметры. Для этого при переходе от отрезка к следующему в качестве априорных оценок параметров модели на отрезке необходимо принять апостериорные оценки этих парамет­ров, полученные на предыдущем интервале.

В качестве примера рассмотрим локальную модель вариаций плотности атмосферы. Для описания вариаций плотности в малом слое

будем основываться на простой экспоненциальной модели зависи­мости плотности от высоты:

где — значение плотности на некоторой фиксированной вы­соте —константа. Фактическая плотность p(h) будет отли­чаться от номинальной на :

где — случайная величина; — случайное отклонение показателя экспоненты от номинального значения . Таким обра­зом, значение определяется двумя случайными величинами: .

Введем величину :

Разлагая в ряд Тейлора по степеням с точностью до чле­нов первого порядка малости, получаем

Выражение (3.104) представляет собой локальную параметри­ческую модель вариаций плотности воздуха. Если для величин записать формальные уравнения формирующих фильтров, то получим векторную форму локальной модели:

Кроме атмосферных возмущений, математические модели кото­рых приведены выше, необходимо учитывать разбросы аэродинами­ческих и конструктивных параметров СА: координат центра масс, площади миделя, аэродинамических коэффициентов. В каждой реа­лизации процесса спуска значения этих параметров не известны заранее и задаются своими априорными статистическими характе­ристиками так, что в целом для совокупности реализаций спуска они — величины случайные.

Будем полагать, что эти случайные величины статистически не­зависимы и подчиняются гауссовским законам распределения. Ана­логичный характер имеют и отклонения начальных условий «вхо­да».

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

с начальным условием q(0) (случайный вектор с характеристика­ми: — диагональная матрица).

Перейдем к построению модели оцениваемого процесса в целом. С целью упрощения рассмотрим движение СА в вертикальной плос­кости. При этом с целью дальнейшего упрощения уравнений дви­жения будем считать, что влияние ветра на движение СА сводится к изменению аэродинамических сил X и У и момента M_z только за счет изменения величины скорости и возникновения дополнитель­ного угла атаки согласно соотношениям

где Vw — воздушная скорость СА; У — скорость СА; — угол на­клона вектора скорости к местному горизонту; — приращение угла атаки; W— горизонтальная скорость ветра.

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

В этих уравнениях введены следующие обозначения:

L — дальность, измеренная по дуге большого круга от точки входа в атмосферу до СА; - воздушный скоростной напор;

m — масса СА; R — радиус Земли, , — коэффициенты нормаль­ной и тангенциальной аэродинамических сил; s_M — площадь миделя; — угол крена, J_х — момент инерции СА относительно оси Ох; М_упр — управляющий момент по крену.

К системе (3.107) следует добавить уравнения формирующих фильтров для ветра и плотности воздуха (3.102) и (3.103) или, если используется локальная модель вариаций плотности, (3.105) вместо (3.103).

Однако в уравнениях формирующих фильтров (3.102), (3.103) и (3.105) независимой переменной является не время, а высота, вследствие чего эти уравнения не могут быть непосредственно ис­пользованы совместно с системой уравнений (3.107). Учитывая диф­ференциальную связь между временем и высотой dh/dt—Vsin , уравнения формирующих фильтров можно привести к виду

Здесь W — вектор 2×1, компонентами которого являются скорость ветра и ее производная по времени; — вектор 2×1, компонентами которого являются вариации плотности воздуха и ее производная по времени; A_w, , B_w, — матрицы 2×2 и 2×1 соответственно вида

— белые шумы.

Кроме уравнений (3.108), к системе (3.107) следует добавить дифференциальные уравнения (3.106) для возмущений, рассматри­ваемых как случайные величины. Вектор q для этих возмущений включает компоненты , где — раз­бросы координат центра масс СА; — разброс площади миделя; — разбросы коэффициентов нормальной и тангенциальной аэродинамических сил. Сформируем расширенный вектор состояния х (15x1) СА с компонентами

где W, p, q — векторы соответствующих размерностей, определяе­мых формирующими фильтрами для этих возмущений.

Вектору х соответствует стохастическое дифференциальное уравнение

где f — вектор-функция, образованная правыми частями диффе­ренциальных уравнений (3.106), (3.107), (3.108); — вектор белых шумов с компонентами ; В — матрица 15×2, устроенная следующим образом

причем символом 0_ij обозначены нулевые блоки соответствующих размерностей; и — управляющее воздействие;

В качестве модели измерений рассмотрим соотношения

т. е. будем считать, что измерению доступны перегрузки СА по свя­занным осям. Кроме того, будем полагать, что в нашем распоряжении имеется дополнительная информация относительно углов тан­гажа и крена, полученная с помощью соответствующих гироскопов, установленных на борту СА. С целью упрощения модели измерений, а также вследствие относительно малого времени спуска считаем, что систематические ошибки измерений отсутствуют и измерения, осуществляемые дискретно, сопровождаются ошибками, представ­ляющими собой дискретную последовательность статистически не­зависимых гауссовских случайных величин , таких, что , — диагональная матрица l×l (l — размерность век­тора измерений ). Таким образом, модель измерений имеет вид

причем формируется на основе оценок вектора . Для обра­ботки информации целесообразно использовать квазилинейный фильтр Калмана, описываемый соотношениями (3.95), (3.96), (3.98), (3.99).

Эффективность алгоритма, описываемого этими соотношениями, может быть проверена имитационным моделированием процесса навигации гипотетического СА, для которого разброс геометриче­ских и конструктивных параметров характеризуется относительным средним квадратичным отклонением в 2%. Ошибки «входа» можно считать соответствующими следующим среднеквадратичным от­клонениям:

Ошибки измерения ускорений примем соответствующими среднему квадратичному отклонению 2-10^-3 м/с²; ошибки измерения угловых величин — средним квадратичным отклонениям 0,01°.

При использовании модифицированного квазилинейного фильт­ра Калмана в данной задаче необходимо учитывать следующее. Формальное применение такого алгоритма может привести к рез­кому увеличению объема вычислений по сравнению с линейным фильтром Калмана и, как следствие, к повышенным требованиям к объему памяти, длине разрядной сетки и быстродействию БЦВМ. Поэтому можно предложить следующую последовательность дейст­вий при построении оценок компонент вектора состояния СА. Счи­тая все компоненты вектора состояния СА, кроме компонент и характеризующих угловое положение СА относительно цент­ра масс, известными, оцениваем последние, используя информацию от гироскопических датчиков. Получив эти оценки с некоторой точ­ностью, переходим к оцениванию компонент, характеризующих по­ложение и скорость центра масс, а также возмущающие воздейст­вия, полагая известными компоненты, характеризующие угловое положение СА. При этом интервал дискретности при оценке ком­понент углового положения СА относительно центра масс должен быть существенно (примерно на порядок) меньше, чем при оценке остальных компонент вектора х.

Такой подход вполне возможен, так как СА оснащен, как пра­вило, системой стабилизации углового движения, и переходные процессы по компонентам вектора состояния, характеризующим угловое положение СА, протека­ют быстрее, чем по компонентам, характеризующим положение и ско­рость центра масс СА. Следует указать также, что точность оценки уг­ла траектории сущест­венно зависит от точности оценки возмущений: ветра, вариаций плотности и разброса аэродинамических коэффициентов.

Поэтому целесообразно вначале проводить оценку составляю­щих возмущений при априорных значениях угла наклона траекто­рии, и только после уточнения возмущений переходить к совместной оценке компонент вектора х. Характер переходных процессов, воз­никающих при таком способе оценивания в квазилинейном фильтре Калмана, показан на рис. 3.9. Апостериорные дисперсии компонент, характеризующих положение и скорость центра масс, а также воз­мущения приведены на рис. 3.10.

Рис. 3.10. Апостериорные дисперсии компонент вектора состояния ЛА

Фильтр второго порядка

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

где — вектор состояния и его компонен­ты; —прогнозированная оценка компоненты , т. е. оцен­ка этой компоненты, полученная по измерениям до момента включительно; —постоянные для данного момента коэффициенты.

Опишем способ определения этих коэффициентов. Обозначим

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

где — матрицы 4×4, диагональные элементы j-й и r-й строк которых единицы, а остальные — нули.

Из условия совпадения функций в точках , найдем коэффициенты полинома:

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

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

С целью получения соотношений, удобных для реализации на ЭВМ, введем в момент t_i,(i — 0, 1, 2, ..., N) расширенный вектор фа­зовых координат в виде

Иначе говоря, будем вводить компоненты вектор-функции поочередно в пределах одного шага измерений t_i в число компонент расширенного вектора состояния оцениваемой динамической систе­мы. В результате окажется, что одна из компонент расширенного вектора состояния, а именно , измеряется непосредственно. При этом компоненты вектора измерений y_i, относящиеся к моменту t_i, будут обрабатываться поочередно; апостериорные характеристики компонент вектора состояния, полученные в результате обработки -й компоненты вектора , будут априорными при обработке (k+1)-й компоненты вектора y_i и т. д.

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

Для того чтобы найти математическое ожидание вектора , достаточно определить априорное математическое ожидание его компоненты . В соответствии с выражением (3.113) имеем

Априорная корреляционная матрица вектора определяется в виде

где — априорная дисперсия компоненты ; — матрица взаимных корреляционных моментов компонент вектора состояния ; и компоненты . В соответствии с выражением (3.113), полагая третьи моменты равными нулю, оп­ределяем

где — центральные моменты четвертого порядка, выражаемые через центральные моменты второго порядка следующим об­разом [34]:

Взаимные корреляционные моменты компонент и , j— 1, 2, ..., п, определяющие элементы матриц всоответствии с (3.113), могут быть записаны при условии, что третьи моменты равны нулю, в виде

где — компоненты вектора .

Таким образом, с помощью приведенных выше соотношений (3.115), (3.117), (3.118), (3.119) статистические характеристики вектора определяющие его априорную к моменту обработки k-й составляющей вектора у_i плотность вероятностей как гауссовскую, определены.

Для определения апостериорных характеристик вектора можно формально применить соотношения (3.42), (3.43), опреде­ляющие линейный дискретный фильтр Калмана, с учетом того, что при принятой схеме обработки матрица H, связывающая k-ю (k = 1, 2, ..., l) компоненту вектора у_i с компонентами вектора представляет собой матрицу-строку 1× (n×1), первые n-столбцов которой нули, а последний — единица:

где — дисперсия ошибки измерения k-й компоненты вектора у_i, y_ik — значение k-й компоненты вектора у_i, — апостериорная ковариационная матрица вектора , т. е. матрица, получающаяся после обработки k-й компоненты вектора у_i, — байесовская оценка вектора , полученная в результате обработки k-й компо­ненты вектора у_i.

Так как обрабатываемые измерения являются скалярными, то соотношения фильтра второго порядка удобно записать в скалярной форме:

В этих соотношениях - соответственно элементы априорной и апостериорной ковариационных матриц, полученные до и после обработки компоненты вектора — соответ­ственно априорная и апостериорная оценки компонент вектора при тех же условиях.

«Начальные условия» для рекуррентных соотношений, описы­вающих фильтр второ-го порядка, те же, что и для описанных выше байесовских алгоритмов. Переход от момента к моменту в фильтре второго порядка происходит также, как и в квазилиней­ном фильтре Калмана, т. е. с помощью соотношений (3.95), (3.96).

3.2.7. Применение фильтра второго порядка для определения орбиты стационарного ИСЗ по данным автономных навигационных средств

Автономная система навигации стационарного искусст­венного спутника Земли является весьма перспективной, поскольку она позволяет разгрузить наземный командно-измерительный комп­лекс, и ее функционирование не зависит от исправности наземных навигационных средств. Однако при реализации автономной нави­гационной системы СИСЗ возникает ряд трудностей, связанных с обеспечением требуемой точности определения орбиты. Дело в том, что автономные навигационные средства обладают сущест­венно более низкой точностью и информативностью, чем наземные.. Низкая информативность автономных (бортовых) навигационных средств выражается прежде всего в малой чувствительности изме­ряемых параметров к изменению положения СИСЗ на орбите. Кроме того, соотношения, связывающие измеряемые (навигацион­ные) параметры с компонентами вектора состояния CPIC3, сущест­венно нелинейны. Наличие больших ошибок линеаризации в соче­тании со значительными систематическими и случайными ошибка­ми измерений делает невозможным применение на борту линейного и квазилинейного фильтров Калмана. В частности, если использо­вать для навигации такие параметры, как углы Солнце—СИСЗ — Земля и Полярная Звезда — СИСЗ — Земля (рис. 3.11), то в про­цессе имитационного моделирования обработки информации с по­мощью квазилинейного фильтра Калмана можно наблюдать эффекг «расходимости» процесса оценивания в том смысле, что отклонения: оценок от истинных значений оцениваемых величин возрастают^ В дальнейшем, как это принято в настоящее время при рассмотре­нии задач автономной на­вигации, измеряемые углы: будем называть СОЗ» (Солнце — объект — Земля) и 303 (Звезда — объект—Земля). Расходи­мость процесса оценива­ния иллюстрируется рис. 3.12 и 3.13. На этих рисун­ках приведены графики-изменения апостериорных: дисперсий координат по­ложения СИСЗ в эквато­риальной декартовой си­стеме координат, а также-графики изменения оши­бок оценок, соответствую-

Рис. 3.11. Схема измерений при автономной навигации СИСЗ

Рис. 3.12. Графики измерения апостериорных дисперсий при автономной нави­гации СИСЗ

щие этим дисперсиям. Исходные данные при моделировании были приняты следующими. Априорные среднеквадратичные отклонения компонент вектора состояния в декартовой экваториальной системе координат =10 км, = 20 км,