Определение состояния спускаемого аппарата по данным автономных измерений с использованием квазилинейного фильтра Калмана
Рассмотрим задачу навигации спускаемого аппарата (СА) с малым и средним аэродинамическим качеством по данным автономных измерений. Пусть информация, необходимая для решения задачи навигации, поступает от датчиков перегрузок (акселерометров), а для обработки информации используется бортовая цифровая вычислительная машина (БЦВМ). В качестве случайных возмущений, оказывающих влияние на процесс оценивания состояния СА, рассматриваются начальные ошибки «входа» в атмосферу, разбросы аэродинамических характеристик, геометрических и конструктивных параметров СА, порывы ветра, вариации плотности атмосферы, ошибки измерения перегрузок СА.
В соответствии с постановкой задачи байесовского оценивания необходимо записать стохастическое дифференциальное уравнение для расширенного вектора состояния, подлежащего оцениванию. В компоненты этого расширенного состояния, помимо переменных, характеризующих положение и скорость СА, включаются также указанные выше возмущения. Следовательно, для возмущений .должны быть записаны стохастические дифференциальные уравнения или уравнения формирующих фильтров.
Особую сложность представляет проблема построения формирующих фильтров для атмосферных возмущений — ветра и вариаций плотности,— поскольку в настоящее время эти возмущения принято рассматривать как нестационарные по высоте полета случайные процессы.
Остановимся на построении формирующих фильтров для атмосферных возмущений подробнее. Ветер принято задавать в виде трех статистически независимых составляющих, одна из которых горизонтальна и лежит в плоскости 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) Kij — элементы корреляционной матрицы Кх. Поскольку матрица Кх и производные ее элементов считаются известными на основании исходных статистических данных, систему уравнений (3.101) можно рассматривать как алгебраическую относительно неизвестных
. Нетрудно видеть, что из п2 уравнений системы (3.101) только п уравнений содержат коэффициенты
— всего n+1 коэффициент. Таким образом, система уравнений (3.101) является неопределенной относительно
и
.
Выбор коэффициентов дифференциального уравнения формирующего фильтра будем осуществлять в процессе минимизации критерия J, характеризующего близость сечений поверхности заданной корреляционной функции и корреляционной функции
случайного процесса на выходе формирующего фильтра.
В качестве критерия можно принять квадратичный: , причем сечения
получить путем интегрирования уравнений
при начальных условиях
где —элементы корреляционной матрицы Кх в соответствующие моменты времени.
Заметим, что и, следовательно, критерий J не зависят b(t).
При заданном порядке формирующего фильтра количество начальных условий должно равняться порядку фильтра. Так, например, если ограничимся первым порядком, то обеспечим совпадение дисперсии случайного процесса на выходе формирующего фильтра с дисперсией заданного случайного процесса. При втором порядке формирующего фильтра, кроме совпадения дисперсий, будет совпадать и первая производная на диагонали квадрата
заданной корреляционной функции с производной корреляционной функции на выходе формирующего фильтра. Если порядок уравнения будет n, то будут совпадать дисперсия и (п—1) производных сечений корреляционных функций
и
. Таким образом, определяем коэффициенты
в результате минимизации функционала
. Минимизация осуществляется, например, методом наискорейшего спуска [28]. Коэффициент b(t) определяем из последнего уравнения системы (3.101).
Анализ корреляционных функций , полученных по экспериментальным данным, проведенный в [20], показывает, что разумным порядком формирующего фильтра для ветра и вариаций плотности атмосферы является второй, так что соответствующее возмущение характеризуется ее величиной и первой производной по высоте полета:
. Производные более высокого порядка вводить нецелесообразно, так как они не обладают достаточной достоверностью.
Таким образом, формирующий фильтр ветра определяется уравнениями
Здесь — обратная высота; h0 — начальная высота спуска; h — текущая высота полета.
Это означает, что уравнения формирующего фильтра (3.102) записаны по обратной высоте, поскольку движение СА происходит при убывании высоты, а поэтому корреляционная функция процесса на входе формирующего фильтра должна_ховпадать с исходной именно при убывании высоты h, т. е. при возрастании обратной высоты hобр.
Аналогично записывается и определяется дифференциальное уравнение формирующего фильтра для Δр:
Возможен и иной подход к построению формирующих фильтров для атмосферных возмущений, основанный на использовании так называемых локальных моделей. Сущность этого подхода сводится к следующему.
При вычислении оценок текущего состояния динамической системы точность статистического описания возмущений важна на том отрезке времени, где производится очередное измерение.
Поэтому возможно использование математических моделей атмосферных возмущений, дающих хорошее приближение лишь на малых отрезках . Такие модели, называемые в дальнейшем локальными, вообще говоря, различны для различных отрезков
. Вид и способ задания локальных моделей определяется в каждом конкретном случае отдельно: это может быть отрезок канонического разложения малой размерности, формирующий фильтр первого порядка. Можно использовать на всем отрезке [0, Т] модель одного вида, периодически изменяя («настраивая») ее параметры. Для этого при переходе от отрезка
к следующему
в качестве априорных оценок параметров модели на отрезке
необходимо принять апостериорные оценки этих параметров, полученные на предыдущем интервале.
В качестве примера рассмотрим локальную модель вариаций плотности атмосферы. Для описания вариаций плотности в малом слое
будем основываться на простой экспоненциальной модели зависимости плотности от высоты:
где — значение плотности на некоторой фиксированной высоте
—константа. Фактическая плотность p(h) будет отличаться от номинальной на
:
где — случайная величина;
— случайное отклонение показателя экспоненты от номинального значения
. Таким образом, значение
определяется двумя случайными величинами:
.
Введем величину :
Разлагая в ряд Тейлора по степеням
с точностью до членов первого порядка малости, получаем
Выражение (3.104) представляет собой локальную параметрическую модель вариаций плотности воздуха. Если для величин записать формальные уравнения формирующих фильтров, то получим векторную форму локальной модели:
Кроме атмосферных возмущений, математические модели которых приведены выше, необходимо учитывать разбросы аэродинамических и конструктивных параметров СА: координат центра масс, площади миделя, аэродинамических коэффициентов. В каждой реализации процесса спуска значения этих параметров не известны заранее и задаются своими априорными статистическими характеристиками так, что в целом для совокупности реализаций спуска они — величины случайные.
Будем полагать, что эти случайные величины статистически независимы и подчиняются гауссовским законам распределения. Аналогичный характер имеют и отклонения начальных условий «входа».
Возмущения типа случайных параметров образуют вектор q, для которого можно записать формальное дифференциальное уравнение
с начальным условием q(0) (случайный вектор с характеристиками: — диагональная матрица).
Перейдем к построению модели оцениваемого процесса в целом. С целью упрощения рассмотрим движение СА в вертикальной плоскости. При этом с целью дальнейшего упрощения уравнений движения будем считать, что влияние ветра на движение СА сводится к изменению аэродинамических сил X и У и момента Mz только за счет изменения величины скорости и возникновения дополнительного угла атаки согласно соотношениям
где Vw — воздушная скорость СА; У — скорость СА; — угол наклона вектора скорости к местному горизонту;
— приращение угла атаки; W— горизонтальная скорость ветра.
Вертикальная составляющая принимается равной нулю. Предполагается, что управление СА в вертикальной плоскости происходит путем разворотов его по крену вокруг вектора скорости. При этом СА балансируется на некотором номинальном угле атаки. При изменении угла крена изменяется проекция нормальной силы на вертикальную плоскость, что позволяет управлять дальностью полета СА. Наличие углов крена разных знаков используется для парирования боковых отклонений СА. Управление углом крена осуществляется с помощью газореактивной системы. С учетом этих допущений уравнения движения СА в вертикальной плоскости в скоростной системе координат имеют вид [19]
В этих уравнениях введены следующие обозначения:
L — дальность, измеренная по дуге большого круга от точки входа в атмосферу до СА; - воздушный скоростной напор;
m — масса СА; R — радиус Земли, ,
— коэффициенты нормальной и тангенциальной аэродинамических сил; sM — площадь миделя;
— угол крена, 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, компонентами которого являются вариации плотности воздуха и ее производная по времени; Aw,
, Bw,
— матрицы 2×2 и 2×1 соответственно вида
— белые шумы.
Кроме уравнений (3.108), к системе (3.107) следует добавить дифференциальные уравнения (3.106) для возмущений, рассматриваемых как случайные величины. Вектор q для этих возмущений включает компоненты , где
— разбросы координат центра масс СА;
— разброс площади миделя;
— разбросы коэффициентов нормальной и тангенциальной аэродинамических сил. Сформируем расширенный вектор состояния х (15x1) СА с компонентами
где W, p, q — векторы соответствующих размерностей, определяемых формирующими фильтрами для этих возмущений.
Вектору х соответствует стохастическое дифференциальное уравнение
где f — вектор-функция, образованная правыми частями дифференциальных уравнений (3.106), (3.107), (3.108); — вектор белых шумов с компонентами
; В — матрица 15×2, устроенная следующим образом
причем символом 0ij обозначены нулевые блоки соответствующих размерностей; и — управляющее воздействие;
В качестве модели измерений рассмотрим соотношения
т. е. будем считать, что измерению доступны перегрузки СА по связанным осям. Кроме того, будем полагать, что в нашем распоряжении имеется дополнительная информация относительно углов тангажа и крена, полученная с помощью соответствующих гироскопов, установленных на борту СА. С целью упрощения модели измерений, а также вследствие относительно малого времени спуска считаем, что систематические ошибки измерений отсутствуют и измерения, осуществляемые дискретно, сопровождаются ошибками, представляющими собой дискретную последовательность статистически независимых гауссовских случайных величин , таких, что
,
— диагональная матрица l×l (l — размерность вектора измерений
). Таким образом, модель измерений имеет вид
причем формируется на основе оценок вектора
. Для обработки информации целесообразно использовать квазилинейный фильтр Калмана, описываемый соотношениями (3.95), (3.96), (3.98), (3.99).
Эффективность алгоритма, описываемого этими соотношениями, может быть проверена имитационным моделированием процесса навигации гипотетического СА, для которого разброс геометрических и конструктивных параметров характеризуется относительным средним квадратичным отклонением в 2%. Ошибки «входа» можно считать соответствующими следующим среднеквадратичным отклонениям:
Ошибки измерения ускорений примем соответствующими среднему квадратичному отклонению 2-10-3 м/с2; ошибки измерения угловых величин — средним квадратичным отклонениям 0,01°.
При использовании модифицированного квазилинейного фильтра Калмана в данной задаче необходимо учитывать следующее. Формальное применение такого алгоритма может привести к резкому увеличению объема вычислений по сравнению с линейным фильтром Калмана и, как следствие, к повышенным требованиям к объему памяти, длине разрядной сетки и быстродействию БЦВМ. Поэтому можно предложить следующую последовательность действий при построении оценок компонент вектора состояния СА. Считая все компоненты вектора состояния СА, кроме компонент
и
характеризующих угловое положение СА относительно центра масс, известными, оцениваем последние, используя информацию от гироскопических датчиков. Получив эти оценки с некоторой точностью, переходим к оцениванию компонент, характеризующих положение и скорость центра масс, а также возмущающие воздействия, полагая известными компоненты, характеризующие угловое положение СА. При этом интервал дискретности при оценке компонент углового положения СА относительно центра масс должен быть существенно (примерно на порядок) меньше, чем при оценке остальных компонент вектора х.
Такой подход вполне возможен, так как СА оснащен, как правило, системой стабилизации углового движения, и переходные процессы по компонентам вектора состояния, характеризующим угловое положение СА, протекают быстрее, чем по компонентам, характеризующим положение и скорость центра масс СА. Следует указать также, что точность оценки угла траектории существенно зависит от точности оценки возмущений: ветра, вариаций плотности и разброса аэродинамических коэффициентов.
Поэтому целесообразно вначале проводить оценку составляющих возмущений при априорных значениях угла наклона траектории, и только после уточнения возмущений переходить к совместной оценке компонент вектора х. Характер переходных процессов, возникающих при таком способе оценивания в квазилинейном фильтре Калмана, показан на рис. 3.9. Апостериорные дисперсии компонент, характеризующих положение и скорость центра масс, а также возмущения приведены на рис. 3.10.
Рис. 3.10. Апостериорные дисперсии компонент вектора состояния ЛА
Фильтр второго порядка
Естественным развитием квазилинейного фильтра Калма-на, основанного на линеаризации уравнений движения и соотношений для измерений в окрестности оптимальной оценки, полученной на предыдущем шаге измерений, является алгоритм, обеспечивающий дальнейший более полный учет нелинейных свойств соотношений для измерений. Такой учет может быть осуществлен аппроксимацией зависимостей, связывающих компоненты вектора состояния и вектора измерений, полиномами наилучшего в среднеквадратичном смысле приближения или полиномами с узлами, определяемыми оценкой, полученной на предыдущем шаге оценивания. Рассмотрим подробнее последний способ аппроксимации. Пусть модели изменения состояния и измерений совпадают с (3.90) и (3.91) соответственно. Представим компоненты вектор-функции gi, входящей в выражение (3.91), в виде
где — вектор состояния и его компоненты;
—прогнозированная оценка компоненты
, т. е. оценка этой компоненты, полученная по измерениям до момента
включительно;
—постоянные для данного момента коэффициенты.
Опишем способ определения этих коэффициентов. Обозначим
где — диагональные элементы апостериорной корреляционной матрицы
, характеризующие точность прогнозированных оценок
. Потребуем совпадения полинома и исходной нелинейной функции
в момент
в четырех точках:
где — матрицы 4×4, диагональные элементы j-й и r-й строк которых единицы, а остальные — нули.
Из условия совпадения функций в точках
,
найдем коэффициенты полинома:
Как и при использовании квазилинейного фильтра Калмана, по мере уточнения оценок компонент вектора состояния размеры области аппроксимации функции полиномом
уменьшаются. Заметим также, что в общем случае аппроксимация функции
может быть осуществлена полиномом, порядок которого выше второго. Однако при этом объем вычислений на каждом шаге измерений существенно увеличивается.
В отдельных случаях точность оценивания может быть повышена путем использования для аппроксимации полинома наилучшего в среднеквадратичном смысле приближения.
С целью получения соотношений, удобных для реализации на ЭВМ, введем в момент ti,(i — 0, 1, 2, ..., N) расширенный вектор фазовых координат в виде
Иначе говоря, будем вводить компоненты вектор-функции поочередно в пределах одного шага измерений ti в число компонент расширенного вектора состояния оцениваемой динамической системы. В результате окажется, что одна из компонент расширенного вектора состояния, а именно
, измеряется непосредственно. При этом компоненты вектора измерений yi, относящиеся к моменту ti, будут обрабатываться поочередно; апостериорные характеристики компонент вектора состояния, полученные в результате обработки
-й компоненты вектора
, будут априорными при обработке (k+1)-й компоненты вектора yi и т. д.
С введением вектора возникает необходимость в определении статистических характеристик этого вектора, априорных по отношению к обрабатываемому измерению. Здесь ограничимся определением априорных математического ожидания и корреляционной матрицы вектора
т. е. характеристик не выше второго порядка.
Для того чтобы найти математическое ожидание вектора , достаточно определить априорное математическое ожидание его компоненты
. В соответствии с выражением (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, yik — значение 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 км,