CC BY
27
УДК 629.78 : 681.51
ПОЛЕТНАЯ КАЛИБРОВКА И ЮСТИРОВКА АСТРОИНЕРЦИАЛЬНОЙ СИСТЕМЫ ОПРЕДЕЛЕНИЯ УГЛОВОГО ПОЛОЖЕНИЯ МАНЕВРИРУЮЩЕГО СПУТНИКА ЗЕМЛЕОБЗОРА
© 2016 Е.И. Сомов, С .А. Бутырин
Самарский научный центр Российской академии наук
Статья поступила в редакцию 06.12.2016
Представляются дискретные алгоритмы фильтрации измерений, полетной калибровки и юстировки астроинерциальной системы определения ориентации маневрирующего спутника землеобзора. Ключевые слова: спутник, определение ориентации, астроинерциальная система, полетная калибровка.
Работа поддержана РФФИ (грант 14-08-01091) и Отделением энергетики, механики, машиностроения и процессов управления РАН (программа фундаментальных исследований № 13).
ВВЕДЕНИЕ
Рассматривается астроинерциальная система определения углового положения (СОУП) маневрирующего космического аппарата (КА). Данная система содержит инерциальный измерительный модуль (ИИМ) в виде жестко закрепленных на корпусе КА гироскопических измерителей приращений угловых квазикоординат КА и астрономическую систему (АС) на основе звёздных датчиков с широким полем зрения, также закрепленных на корпусе спутника. СОУП входит в состав бесплатформенной инерциальной навигационной системы (БИНС), которая решает общую задачу навигации - определения как ориентации, так и местоположения подвижного объекта [1 - 4], в данном случае спутника. Далее считаются известными понятия и методы обработки информации в БИНС, включая модели погрешностей ИИМ и АС [3 - 6].
Предложенный в [7, 8] подход к обработке сигналов в астроинерциальной СОУП развивается для маневрирующих КА радиоэлектронного и оптико-электронного наблюдения Земли. Проблемы такой обработки связаны с интегрированием кинематических уравнений по информации только о векторе дискретных приращений угловых квазикоординат, получаемого ИИМ при наличии шумов измерений, с калибровкой (идентификацией и компенсацией вектора дрейфа ИИМ bg и погрешности m масштабного коэффициента измерения вектора угловой скорости ш ) и с юстировкой [9] - идентификацией
Сомов Евгений Иванович, кандидат технических наук, доцент, ведущий научный сотрудник отдела «Динамика и управление движением» СамНЦ РАН. E-mail: [email protected]
Бутырин Сергей Анфимович, кандидат технических наук, старший научный сотрудник отдела «Динамика и управление движением» СамНЦ РАН. E-mail: [email protected]
и компенсацией погрешности взаимной угловой установки систем координат ИИМ и АС на основе сигналов АС с периодом дискретности To, существенно большего периода дискретности Tp выходных сигналов СОУП. В процессе решения таких задач многими авторами использовались кинематические параметры в виде кватерниона Л = (Л,0, , матрицы ориентации C, вектора Эйлера ф = e 9 , вектора конечного поворота 0 = 2e tg(9/2) и др. При этом для медленных угловых движений КА с малой вариацией угла 9 и практически неизменным направлением орта Эйлера e на периоде дискретности To интегрирование кинематических соотношений для вектора Эйлера и получение значений кватерниона Л k - Л (tk ) выполнялось по схеме
= i Г = fi+1 Ш (T)dT - Int(tk To, <D (t) ,
Jtk
Фk + % = Фм ^ Ck+1 ^ Лk+1 , tk+1 = tk + To ,
где 5фk = 59kek и k e N0 - [0,1,2,...). Угловое движение маневрирующего КА землеобзора осуществляется на последовательности временных интервалов выполнения маршрутов сканирующей съемки и быстрых поворотных маневров (ПМ) с переменным направлением вектора Ш , модуль которого может достигать 5 град/сек. Для вектора Эйлера дифференциальное уравнение имеет вид
ф = Ш + 2ф хш + f (9)ф х (ф хш), где скалярная функция f (9) = (1 -9 sin 0/(2(1 - cos 9)) / 0¿ [10], при малых значениях угла 9 функция f (9) « 1/12 . На основе этого уравнения нетрудно понять физический смысл возникающего конинга (coning) - известного эффекта некоммутативного поворота, обусловленного изменением именно направления вектора Ш в процессе выполнения ПМ [10 - 12].
В статье рассматривается классическая схема СОУП в составе БИНС, где основным измерителем является ИИМ с гироскопическими датчиками угловой скорости (ДУС) при выходном векторе приращений угловых квазикоординат - интегралов от проекций измеренного вектора и на периоде дискретности Т , а сигналы АС используются для его калибровки и юстировки с периодом дискретности Т0, кратном периоду Т . Здесь применяются методы полиномиальной аппроксимации и интерполяции, выполняется численное интегрирование кинематического уравнения для вектора модифицированных параметров Родрига (МПР) С = е 1§(0/4) [13], который связан с кватернионом Л прямыми С = Х/(1 + ^0) ( Л^С ) и обратными
А, = 2а /(1 + а2); ^ = (1 -а 2)/(1 + а 2)( с^Л ) соотношениями. Кинематические прямые и обратные уравнения для вектора МПР имеют вид
С = Га (а,ш) = 1(1 -а 2)ш+1ахш
+ 1а <С,и >; ш = Г и (а, С) = 4[(1 -а 2)С - 2(ахС) + 2а < С,а >]/(1 + а2)2
при стандартном обозначении < а, Ь > скалярного произведения векторов а и Ь.
ПОСТАНОВКА ЗАДАЧИ
Аналогично [7, 8] вводятся инерциальный базис I; базис В, связанный с корпусом КА; базис 8 сенсора (телескопа, радиолокационной антенны); виртуальный базис АС А, вычисляемый на основе информации от звёздных датчиков, и виртуальный базис ИИМ О, вычисляемый на основе информации от ДУС с выходным вектором приращений угловых квазикоординат. Будем для простоты считать, что базисы В и 8 совпадают, известна информация о взаимном угловом положении базисов А и 8, полученная одним из известных методов, в частности для КА оптико-электронного землеобзора см. [14, 15]. Пусть исходная измерительная информация обрабатывается в ИИМ с частотой ~ 3 кГц и в результате от ИИМ поступают измеренные значения вектора
приращений квазикоординат 1 , ^ е N с периодом Т < Т0, а от АС - измеренные значения кватерниона Л лт к с периодом То :
1 тш=, т, ш т(о +5 п;
ш 1(1) = (1 + т) 8 Л (ш () + Ьё) ; (1)
л m * =л t °
л n.
Здесь ш m(t) представляет измеряемый
ДУС вектор угловой скорости КА в базисе G с учетом неизвестных малых и медленных вариаций вектора bg = bg (t) дрейфа ИИМ по угловой скорости; при обозначениях Са = cos а , Sa = sin а ортогональной матрице SЛ (t) погрешности взаимной угловой установки систем координат ИИМ и АС соответствует кватернион
ЛЛ = (C д eлS Д ) с неизвестными ортом
„А (~ДФ~Д Ф о-
e = {ex, ey, ez } и углом Эйлера ф ; скалярная
функция m = m (t) представляет неизвестную медленную вариацию масштабного коэффициента ИИМ, учтены центрированные дискретные
гауссовские шумы 8 П со среднеквадратичным
отклонением (СКО) аь и ЛП с СКО а3 в выходных сигналах ИИМ и АС соответственно.
Будем считать, что матрица Sд (t) близка к единичной матрице 13 и приближенно представляется в виде Sд = 13 + [Ах], где вектор «скосов» А = {Д x ,Д y ,Д z} систем координат ИИМ
и АС имеет компоненты Дi = e¿ ф , i = x, y, z. Предполагается малая вариация масштабного коэффициента ИИМ, например | m (t) |< 0.01, что позволяет считать справедливым 1 - m2 = 1.
Задача состоит в разработке алгоритмов дискретной фильтрации измерений ИИМ с периодом Tq, калибровки и юстировки СОУП по сигналам
АС с периодом To с получением оценок bg, mk и SД, а также цифровой фильтрации выходных сигналов СОУП с периодом Tp = t¡+1 - tl, кратном периоду To, в общем случае с фиксированным запаздыванием Td относительно моментов времени tk, с получением дискретных значений
кватерниона Лf и вектора угловой скорости ш^ , l е N0, которые используются для управления ориентацией маневрирующего КА землеобзора, при наличии быстрых ПМ спутника с существенно переменным направлением вектора его угловой скорости на ограниченных временных интервалах.
ПОДХОД К РЕШЕНИЮ ЗАДАЧИ
В модели (1) вектор дрейфа bg (t) представляет аддитивное возмущение, а ортогональная
матрица S Д (t) и вариация масштабного коэффициента m (t) - мультипликативные возмущения, идентификация которых является наиболее сложной проблемой. В известных методах решения таких нелинейных задач для корректируемых БИНС [9,13] используются нестационарные фильтры Калмана, которые «работают» в малой окрестности тестовых пространственных угловых движений КА, которые необходимы для обеспечения условий наблюдаемости. Применяемый
авторским подход основывается на дискретном фильтрации доступных измерений и следующих положениях [7, 8]:
• для калибровки ИИМ с периодом To в отношении дрейфа bg применяется дискретный стационарный нелинейный фильтр Луенбергера с параметрами, при которых обеспечиваются показатели качества оценки bg , близкие к аналогичным показателям соответствующего фильтра Калмана, причём такая калибровка выполняется на временных интервалах, где модуль угловой скорости КА не превышает 1 град/с и измерения АС являются вполне точными;
• при ПМ спутника с модулем угловой скорости от 1 град/с до 5 град/с на временных интервалах длительностью до 120 с коррекция СОУП временно отключается и оценка углового положения КА продолжает выполняться с использованием прогноза вариации вектора дрейфа ИИМ, который формируется на основе анализа тренда
оценки bg этой вариации в «скользящем окне» на предшествующем интервале работы СОУП с астрономической коррекцией;
• юстировка СОУП (оценка SА) и её калибровка в отношении вариации масштабного коэффициента с определением оценки m выполняются off-line на основе сопоставления значений вектора ш , которые с применением дискретной фильтрации оцениваются автономно и одновременно по сигналам ИИМ (значений вектора Шg) и АС (значений вектора шa).
Отметим, что в предлагаемом подходе полётная калибровка и юстировка СОУП выполняются на фоне решения целевых задач маневрирующего спутника, без каких-либо тестовых пространственных угловых движений КА.
ФИЛЬТРАЦИЯ ИЗМЕРЕНИИ И КОМПЕНСАЦИЯ ДРЕЙФА ИИМ
Будем предусматривать формирование оценок bg, Sk и mk, постоянных на периоде To, когда оценка bg обновляется on-line, т.е. в
каждый момент времени tk, а оценки S ^ и mk регулярно формируются off-line, т.е. на основе обработки доступной измерительной информации, накопленной на длительных временных интервалах. При дискретной фильтрации измерений вектора приращений квазикоординат с целью подавления дискретного шума ИИМ впервые применяется оригинальная двухпроходная технология фильтрации - сочетание аппроксимации измеренных данных ig ffl полиномом 3-го порядка
1
g ш
СО
к 4 в скользящем окне 9 измерений по методу наименьших квадратов (МНК) и сплайновой интерполяции по центрам двух смежных скользящих окон полиномом 5-го порядка 1 (т) для локального времени т = t - кТ0 е [0, Т0 ].
Разработанная технология иллюстрируется схемой на рис. 1, где синими «звездочками» отме-
1 5 дискретного измерения квазикоординаты ИИМ в моменты времени ts (здесь указаны значения индекса 5), зелеными пунктирными линиями выделены полиномы
51 8Ш (т) 3-го порядка, которые отражают погрешности аппроксимации в скользящем окне 9 измерений, и бордовыми линиями представлены гладко сопряженные полиномы 518И (т) 5-го порядка, которые соответствуют погрешностям процедуры интерполяции по центрам двух смежных скользящих окон. В нижней части
х10
0.1
; -0.1
* si0)0 — sTmg -3>а -
+
+ + I я + +
+. / * л .......к\..... ...А.... + \ + V У ......./. \ * t ........\........ kiV
у 1 J 1+ / t \ Г
+ + и*/ + V + г* + + +
+ • +
+ *
12
16 20 24 28
32 36 40 44 48 52 56 60
64
63
Рис. 1. Погрешности двухпроходной техники фильтрации измерений ИММ
этого рисунка приведена погрешность 5ш 8(т) оценки ш8(т) угловой скорости в базисе О , которая получается по явным аналитическим соотношениям и поэтому строго согласована с
оценкой ¡8И (т) квазикоординаты ИИМ.
При компенсации влияния вектора дрейфа
Ь8 и погрешности масштабного коэффициента т в базисе О получается оценка векторного полинома (сплайна) приращений квазикоординат
1Г (т) = (1 - т к )(~8 * (т) - Ь 8 т).
В базисе А астрономической системы оценка непрерывного вектора приращений квазикоординат 1 © (т) вычисляется по формуле
(т) = ( 8 £ )118 © (т), (2)
а непрерывная оценка © к(т) значений вектора угловой скорости ш (^) на к -ом интервале времени t е [¿к, tk+1 ] с получается аналитическим дифференцированием этого векторного сплайна, что обеспечивает полное соответствие векторных
полиномов 1 к (т) и © к(т) ■
Идентификация дрейфа Ь8 ИИМ выполняется с применением дискретного наблюдателя Луенбергера при периоде дискретности То. Для определения структуры этого наблюдателя рассмотрим сначала упрошенную линеаризованную модель ИИМ, где 8А = 13, т ^) = 0, орт е к (т) вектора МПР о к (т) = е к (т)18(6к (т)/4) неизменен в базисе I, т.е. ек (т) = ек на к -ом временном интервале с локальным временем т = t - кТо е [0, То], а вектор дрейфа Ь8(т) считается постоянным на этом временном интервале, т.е. ё Ь8 (т) / ё т = Ь8 (т) = 0 Vт е [0, То ] . Если при этом шум измерений ИИМ и АС не учитывается, то дискретное изменение значений вектора рк = 4 о к « фк с периодом То описывается соотношением
Рк+1 = Рк +1© + Ь| То. (3)
При отсутствии шума измеряемые значения
р а = рк вектора рк аналитически определяются по соотношениям А ^ о и р = 4 о на основе измеряемых дискретных сигналов АС в виде кватерниона А атк = Ак . Здесь отсутствует какая-либо
измеряемая информация о значении вектора Ь \ , поэтому решается локальная задача синтеза структуры асимптотически устойчивого линейного дискретного наблюдателя Луенбергера [16] для
получения дискретных оценок рк и Ь| векторов рк и Ь8 соответственно. При обозначениях {•} и [•] для столбца и строки представим модель (3) в стандартной форме дискретной системы с векторами состояния хк = {рк,Ь8}, входа ик = 1 © и измеряемого выхода у к = р к в виде
х
к+1
где матрицы
А ё =
= А ёх к+в ёи к; у к = х к, (4)
о
ТI I,
вё ={1з,0};
Сё = [13,0].
При символе транспонирования (•)' матрица
наблюдаемости Ц = [Сё, А Сё] этой системы
имеет гап8 (Ц ) = 6, так как ёй (Ц) = Т'3 ■ Поэтому линейная дискретная система (4) полностью наблюдаема [16] и структура линейного дискретного наблюдателя Луенбергера для получения
оценки Xк = {рк, Ь8 } вектора хк имеет вид
х к+1 = А ё х к + в ёи к + О ё (у к - С ёх к).
Здесь постоянная матрица Оё назначается из условия асимптотической устойчивости тривиального решения 8 к = х к - ххк = 0 дискретной системы
8к+1 = Аё8к - (Аё - ОёСё ) 8к, (5) что соответствует принадлежности спектра матрицы Аё внутренности единичного круга, т.е. условию 11 р |< 1 Vp = 1 ^ 6, где - корни характеристического уравнения ёе1;(II6 - А°л ) = 0 . Рационально выбрать структуру матрицы Оё в виде О л = 13, g 213}, где скалярные коэффи-
_о _о
циенты g1 и g2 определяются аналитически по каноническим представлениям линейных дискретных систем управления [17] с обеспечением одинаковой степени устойчивости системы (5) по каналам идентификации, т.е. с обеспечением одинаковых показателей сходимости оценок компонентов векторов рк и Ь8. В этом случае характеристическое уравнение каждого канала идентификации имеет вид 12 + а^ + а9) = 0,
где а', = -(2 - gl0) и а*0 = 1 - - Т^о°).
При желаемом спектре канала в виде г*,2 = ехр(УТ0 (-£, ± )
(7 = V-! ) коэффициенты характеристического уравнения вычисляются по фор-
мулам Виета
а0*ё = 4 = ехР(-2^уТ0) .
= -(1* +1*) = -2 ехр( -ОД ^(уТо^/Л^)
Это позволяет в явном виде определить коэффициенты дискретного наблюдателя Луенбергера
„о о « „ё _ ё _ „
g1 и g2 из условий а1 = а1 и а0 = а0 . Так, при выборе значений у = 3/(^ТГ) и Е, = 0.95 при отсутствии шума измерений гарантируется время Тг «регулирования» - входа оценок координат рк и компонентов Ь8к дрейфа ИИМ в соответствующие стандартные «трубки точности» размером ± 5% от начального рассогласования.
I
3
Искомая структура линейного дискретного наблюдателя Луенбергера для упрощенной линейной модели (4) представляется в виде
р к+1 = р к +1 © + 8о ар к;
Ь 8+! = Ь| + g о Ар к;
к+1
ар к+1 = р к+1
(6)
к+1
где постоянные коэффициенты g10 и g0 вычисляются по явным соотношениям.
При учете шумов измерений АС наблюдатель типа (6), где вектор р а вычисляется по измерениям АС в виде А т к = А к о Ап с кватернионом Ап дискретного центрированного гауссовского шума, некоррелированного по каналам наблюдения, может представляться также в форме фильтра Калмана. При этом только для установившегося режима в каналах идентификации аналитически определяются постоянные матрицы ковариации и коэффициентов фильтра Калмана, в переходных же режимах идентификации приходится численно решать матричное уравнение Риккати. Вычислительные ресурсы СОУП ограничены, поэтому здесь применяется наблюдатель дрейфа ИИМ в форме Луенбергера и его параметры подбираются из условий обеспечения показателей качества оценки Ь8 , близких к аналогичным показателям фильтра Калмана в установившемся режиме.
Возвратимся к исходной задаче идентификации дрейфа ИИМ Ь8 с помощью наблюдателя Луенбергера в СОУП маневрирующего КА,
когда орт Эйлера ек (т) вектора МПР ок (т) изменяет свое направление в базисе I на к -ом временном интервале с локальным временем т = t - кТо е [0, То ]. Оценка изменения углового положения КА на этом интервале времени достигается численным интегрированием векторного дифференциального уравнения
о к (т) = Р о (о к (т),ш к (т) (7)
методом ОБЕ45 [18,19] с формированием оценки вектора МПР ок (т) . Вычислительный дрейф такой процедуры интегрирования был детально исследован в [7], например при коническом движении вектора угловой скорости с модулем 0.5 град/сек и шаге интегрирования 0,015625 сек вычислительный дрейф не превосходит 0.00025 угл. сек за время 720 сек. Начальные условия для уравнения (7) формируется непосредственно по сигналам АС только в момент времени t0 (при включении СОУП), в остальных случаях начальные условия такого интегрирования на очередном к -ом временном интервале определятся по сигналам дискретного наблюдателя Луенбергера.
Не указывая индекс к временного интервала, приведем соотношения для вычисления непрерывных погрешностей при получении оценок углового положения и угловой скорости КА по предлагаемой методике. Здесь фактическое угловое движение КА задано кватернионом А(т) и вектором угловой скорости ш (т), а соответствующие оценки этих переменных определяются
кватернионом А(т) и вектором ш (т), причем выше указано, что непрерывная оценка ш (т) получается аналитическим дифференцированием
векторного сплайна 1ш (т). В результате интегрирования (7) получаются значения о(т) векторной функции МПР, на основе которых с помощью
преобразования о(т) ^ А(т) формируется
оценка А(т) изменения кватерниона. Непрерывные погрешности оценки углового положения
определяются кватернионом А = А о А = (е0, е) с вектором е = {ег.}, вектором параметров Эйлера А = {е0, е} , матрицей погрешности ориентации Се = 13 - 2[ех]Ре , где матрица
Рис. 2. Погрешности непрерывной оценки и цифровой фильтрации выходных сигналов СОУП
Q e = 13e0 + [ex], вектором модифицированных параметров Родрига (МПР) oe = e /(1 + e0) и
столбцом угловых ошибок 5ф = {5фД = 4сe. Непрерывная погрешность оценки вектора угловой скорости 5ш(х) определяется соотношением 5ш(х) = ш(т) - Сеш(т). На рис. 2 бордовыми линиями схематично представлены непрерывные погрешности при оценке углового положения 5ф и угловой скорости 5ш КА в базисе A астрономической системы.
Итак, пусть в момент времени t = tk имеются оценки матрицы SA и погрешности масштабного коэффициента mk измерений вектора угловой скорости, постоянные на периоде To (методика получения таких оценок off-line
представлена ниже), оценки кватерниона Лk
и вектора дрейфа b g ИИМ, а также измеренное АС фактическое угловое положение КА в виде
кватерниона Лm k . Пусть в момент времени tk сформированы кватернион AP, (gO, Qk ) и вектор Apk(gO,Qk) , где кватернион Qk - (?0k,qk) - (^/2>eq S,^) - Qk(el, 9k) = k ° Лk . В этот же момент времени tk с помощью преобразования Лk ^ сk определяется начальное условие 6 k (0) - 6 k численного интегрирования уравнения (7) для получения оценки вектора МПР 6k (т) на k -ом интервале. Такое интегрирование выполняется для указанной выше последовательности гладко сопряженных векторных сплайнов СО k (т) . Здесь рационально применять параллельные вычисления для аппроксимации и интерполяции при получении в базисе G
векторного сплайна ig® (т) , его представления iш (т) (2) в базисе A астрономической системы с аналитическим формированием оценки СО k(т) и собственно численного интегрирования (7).
В момент локального времени Т = TO в результате интегрирования получается значение
вектора МПР с k+1 = 6 k (To) , по которому с помощью преобразования СС k+1 ^ Rk+1 определяется
значение кватерниона R k+1. Разработанный нелинейный дискретный наблюдатель Луенбергера, где используются численное интегрирование (7) в сочетании с кватернионными и векторными соотношениями, в конечном итоге имеет вид
Л k+1 = R k ° APk (go, Q k);
b g+i = b g + Ap k (g o, Q k);
_o
APk+i( go, Qk+1) = Q k+i (e l+i, g°tyk+1); Ap k+1( g o, Q k+1) = 4 g 2°С l+1,
(8)
где значение вектора МПР С q+1 определяется по кватерниону Qk+1 явными соотношениями.
На завершающем этапе фильтрации дискретных измерений ИИМ с компенсацией его дрейфа выполняется формирование цифровых выходных
сигналов СОУП - значений кватерниона Лj, которые получаются по значениям вектора МПР С j с помощью преобразования 6 j ^ Лj, и вектора угловой скорости СО j, I е N0 с заданным периодом Tp > Tq, кратном периоду Tq. Погрешности применяемой фильтрации выходных сигналов СОУП представлены на рис. 2 черными линиями.
Несмотря на то, что описанная дискретная фильтрация измеряемых сигналов ИИМ может выполняться по схеме параллельных вычислений, здесь возникает фиксированное временное запаздывание Td = 4 Tq относительно моментов времени tk. Авторский метод аналитического учета такого запаздывания при синтезе цифровых алгоритмов управления ориентацией КА представлен в [20].
ЮСТИРОВКА И КАЛИБРОВКА МАСШТАБНОГО КОЭФФИЦИЕНТА
Как отмечено выше, юстировка СОУП (вычисление оценки SA) и определение оценки гп погрешности m масштабного коэффициента выполняются off-line на основе сопоставления значений вектора угловой скорости С спутника, которые восстанавливаются автономно по сигналам ИИМ (векторы шg) и АС (векторы ОСa ), полученных в одинаковые моменты времени.
На фоне выполнения целевых задач спутника формируется группа трёх временных «скользящих окон» с номерами j = 1,2,3, с началом в моменты времени m]T0 и длительностью UjTo, где Uj е [12... 18], причем каждое окно содержит временной интервал класса МТ а, где | ш (t) |< 1 град/сек, выполняются рекуррентная калибровка по сигналам астрономической системы и компенсация дрейфа bg ИИМ . В каждом j -ом окне длительностью UjTo:
(i) явным дифференцированием полиномов
i |ffl (т) в базисе G получаются оценки вектора
угловой скорости С jg с периодом Tp ;
(ii) по измеренным АС значениям кватерниона Лmk формируется набор значений вектора МПР С лтk и выполняется его аппроксимация по МНК полиномом 5-го порядка 6* (т j) для соответствующего локального времени т ; = t - m]To е [0, n}To ], j = 1,2,3 ;
(iii) по обратному уравнению для вектора МПР определяется аналитическая зависимость
оценки ©а(т) вектора угловой скорости внутри ] -го окна и формируются значения векторов ш га с периодом Тр.
Далее с помощью переиндексации формируются два массива значений векторов ш 8 и ш га в моменты времени t = ti с периодом Тр. Для этих векторов вычисляются значения модулей
ш 8 =|ш81, ев а=| ш аI и ортов еш. =(0 8/ш 8, еш. = ш а/со а, если ш8 > 0 и ша > 0.
Задача юстировки, которая сводится к определению взаимной ориентации ортогональных базисов A и G на основе вычисленных значений ортов ai = е ш в базисе A и ортов g. = е ш. в базисе G, формулируется так: най-
А А
ти оптимальную оценку S ортогональной матрицы SА с определителем, равным +1, которая минимизирует квадратичный критерий
Ь = 1 Еа. | gi - SАai |2, где неотрицательные числа а. являются весовыми коэффициентами с условием Еа. = 1. Как известно, решением этой
А
задачи является оптимальный кватернион А , который взаимно-однозначно связан с искомой
А
ортогональной матрице S и определяется как собственный вектор матрицы K с наибольшим собственным значением qm, т.е. кватернион Аа определяется по алгоритму
B = Е а. g .а;; д = б + Б1; z = Еа. g. х ;
K
1r B z1 z Q -13t B
; KAA= qm ЛА. (9)
ИИМ для моментов времени tk > (т3 + п3)То вплоть до начала очередного временного интервала класса МТ а, когда включается астрономическая коррекция ИИМ.
ИМИТАЦИЯ РАБОТЫ СОУП
Пусть для КА землеобзора на солнечно-синхронной орбите высотой 720 км имеется задание на съемку окрестностей российских столиц (рис. 3) в такой последовательности временных интервалов при отсчете времени от момента прохождения восходящего узла орбиты:
• ориентация КА в ОСК Vt е [0,870) с;
• ПМ1 длительностью 40 с Vt е [870,910) с;
• маршрут М1 съемки с выравниванием продольной СДИ длительностью 40 с Vt е [910,950) с при начальной точке в гринвичской системе координат с долготой Ь = 37.8 град, широтой В = 54.3 град и геодезическим азимутом А = -12 град, ширина полосы захвата 41 км, длина 271 км;
• ПМ2 длительностью 55 с Vt е [950,1005) с;
• маршрут М2 трассовой съемки длительностью 15 с Vt е [1005,1020) с начальными значениями углов крена ф1 = 8 град и тангажа ф3 = 0, ширина полосы захвата 40 км, длина 102 км.
На рис. 3 представлен закон наведения КА в инерциальном базисе для выполнения указанного задания, где цветами выделены компоненты векторов о (), ш ^) и 8 ) - синим цветом по крену, зеленым по рысканию и красным по тангажу, а модули векторов ш ($) и 8 ($) отмечены черным цветом.
Соотношения (9) представляют известный алгоритм QUEST (Quaternion Estimation) для по- ЛА
лучения оптимальной оценки кватерниона ЛА
по векторным измерениям.
Для калибровки погрешности m масштабного коэффициента вычисляется последовательность значений = 1 - (cog / со a), искомая оценка m получается их осреднением и применяется как fhk до завершения следующей калибровки. На временных интервалах, не принадлежащих классу МТ а, астрономическая коррекция выключается и оценка углового положения КА продолжает выполняться с использованием прогноза
вариации оценки lb8 вектора дрейфа bg . Такой прогноз формируется на основе анализа тренда этой вариации на предыдущем интервале работы СОУП с астрономической коррекцией. При этом выполняется аппроксимация по МНК набора оценок вектора дрейфа ИИМ, сформированных в последнем «скользящем окне» (j = 3), векторным полиномом 2-го порядка, на основе которого реализуется прогноз вариации вектора дрейфа
Рис. 3. Маршруты съемки двух столиц на карте
При компьютерной имитации работы СОУП в процессе выполнения указанных маршрутов
съемки были приняты СКО оа = 0.3 угл. сек
/л/Гц на частоте 1 Гц для шума измерения АС
и СКО оъ = 0.001 угл. сек /д/Гц на частоте 128 Гц для шума измерения вектора псевдокоординат 1, а параметры наблюдателя Луенбер-гера выбраны из условия близости его качества
Рис. 4. Векторный сплайновый закон наведения
к качеству оценки дрейфа ИИМ фильтром Калмана, когда тестовый вектор дрейфа ИИМ Ь8 = {1,- 0.8, 0.3} угл.сек / сек восстанавливается за время « 60 секунд. На рис. 5 представлены
погрешности 5Ь8 = Ь8 - Ь8 оценки вектора дрейфа ИИМ в базисе О. В моменты времени tl с периодом Тр = 0.125 с (частота 8 Гц) погрешность оценки ориентации представляется
вектором 5г = 5ф1 = 4а ^. Цифровая фильтрация
вектора 5{ с частотой 8 Гц дает вектор 5í (t), его компоненты изображены на рис. 6. Погрешности СОУП при определении вектора угловой скорости и цифровой фильтрации его значений с частотой 8 Гц приведены на рис. 7 и 8 для компонентов
вектора , см. рис. 2.
Прогноз вариации оценки Ь8 без астрономической коррекции СОУП выполнялся по набору
оценок со 8 и соа, сформированных в последнем «скользящем окне» длительностью 20 секунд.
С применением разработанных и описанных выше процедур юстировки и калибровки СОУП получены оценка компонентов вектора «скосов»
А с точностью « 3 угл. сек и оценка т погрешности масштаба измерений вектора угловой скорости с точностью ~ 0.025 % .
Рис. 5. Погрешности оценки вектора дрейфа ИИМ в базисе О
Рис. 6. Угловые погрешности СОУП при определении ориентации КА
Рис. 7. Погрешности СОУП при определении вектора угловой скорости
Рис. 8. Погрешности СОУП при определении угловых скоростей, фрагмент рис. 7. ЗАКЛЮЧЕНИЕ
Кратко представлены новый оригинальный метод дискретной двухпроходной фильтрации измерений ИИМ, методики полетной калибровки и юстировки астроинерцильной системы для определения ориентации и угловой скорости маневрирующего спутника землеобзора. Приведены результаты компьютерной имитации, демонстрирующие эффективность разработанных алгоритмов.
СПИСОК ЛИТЕРАТУРЫ
1. Ишлинский А.Ю. Ориентация, гироскопы и инер-циальная навигация. М.: Наука. 1976.
2. Теоретические основы проектирования информационно-управляющих систем космических аппаратов / В.В. Кульба, Е.А. Микрин, Б.В. Павлов, В.Н. Платонов. М.: Наука. 2006.
3. Бранец В.Н., Шлыглевский И.П. Введение в теорию бесплатформенных инерциальных навигационных систем. М.: Наука. 1992.
4. Titterton D.H., Weston J.L. Strap-down inertial navigation technology. 2nd ed. Cornwall: MPG Books
Ltd. 2004.
5. Анучин О.Н., Комарова И.Э., Порфирьев Л.Ф. Бортовые системы навигации и ориентации искусственных спутников Земли. СПб: Концерн «ЦНИИ Электроприбор». 2004.
6. Матвеев В.В., Распопов В.Я. Основы построения бесплатформенных инерциальных навигационных систем. СПб: Концерн «ЦНИИ Электроприбор». 2009.
7. Somov Ye.I. Multiple algorithms for filtration, integration and calibration of a strap-down inertial system for a spacecraft attitude determination // Proceedings of 16th Saint Petersburg International Conference on Integrated Navigation Systems. 2009. P. 110-112.
8. Somov Ye.I., Butyrin S.A. Digital signal processing, calibration and alignment of a strap-down inertial system for attitude determination of an agile spacecraft // Proceedings of 17th Saint Petersburg International Conference on Integrated Navigation Systems. 2010. P. 81-83.
9. Pittelkau M.E. Kalman filtering for spacecraft system alignment calibration // Journal of Guidance, Control and Dynamics. 2001. Vol. 24. No. 6. P. 1187-1195.
10. Bortz J.E. A new mathematical formulation for strapdown inertial navigation // IEEE Transactions on Aerospace and Electronic Systems. 1971. Vol. 7. No. 1.
P. 61-66.
11. New procedure for deriving optimized strap-down attitude algorithms / V.Z. Gusinsky, M.V., Lesyuchevsky Yu.A Litmanovich, H. Mussoff, G.T. Schmidt // Journal of Guidance, Control and Dynamics. 1997. Vol. 20. No. 4. P. 673-680.
12. Бабич О.А. Исследование некоммутативных поворотов в алгоритмах ориентации БИНС методом аксоидов // Труды МИЭА. 2013. Вып. 6. С. 18-34.
13. Yong K., Jo S., BangH. A modified Rodrigues parameter-based nonlinear observer design for spacecraft gyroscope parameters estimation // Transaction of the Japan Society for Aeronautical and Space Science. 2012. Vol. 55. No. 5. P. 313-320.
14. Somov Ye. I., Butyrin S.A., Skirmunt V.K. In-flight alignment calibration a space telescope and a star tracker cluster // Proceedings of the jubilee 15th Saint Petersburg International Conference on Integrated Navigation Systems. 2008. P. 139-143.
15. Somov Ye.I.,Butyrin S.A. In-flight alignment of a space telescope and a star tracker cluster at a scanning observation of the Earth marked objects // Proceedings
of 19th Saint Petersburg International Conference on Integrated Navigation Systems. 2012. P. 242-244.
16. Kwakernaak H., Sivan R. Linear Optimal Control Systems. Wiley-Interscience. 1972.
17. Сомов Е.И. Структурные свойства и канонические представления линейных систем управления. Изд-во Иркутского гос. ун-та. 1993.
18. Shampine L.F. Some practical Runge-Kutta formulas // Mathematics of Computation. 1986. Vol. 46. No. 173. P. 135-150.
19. Mathews J.H., Fink K.D. Numerical methods using MATLAB. 3-rd Ed. Prentice Hall. 1999.
20. Сомов Е.И. Робастная стабилизация упругих космических аппаратов при неполном дискретном измерении и запаздывании в управлении // Известия РАН. Теория и системы управления. 2001. № 2. С. 124-143.
21. In-flight calibration of the attitude determination systems for information mini-satellites / Ye. Somov, S. Butyrin, H. Siguerdidjane, C. Hajiyev, V. Fedosov // Proceedings of 19th IFAC Symposium on Automatic Control in Aerospace. 2013. P. 393-398.
IN-FLIGHT CALIBRATION AND ALIGNMENT OF ASTROINERTIAL SYSTEM FOR ATTITITUDE DETERMINATION OF A MANEUVERING LAND-SURVEY SATELLITE
© 2016 Ye.I. Somov, S.A. Butyrin
Samara Scientific Center, Russian Academy of Sciences
We have presented discrete algorithms for a filtering of measurements, in-flight calibration and alignment of astroinertial system for the attitude determination of a maneuvering land-survey satellite. Keywords: satellite, attitude determination, astroinertial system, in-flight calibration.
Yevgeny Somov, Candidate of Technics, Associate Professor, Leading Research Fellow at the Dynamics and Motion Control Department. E-mail [email protected] Sergey Butyrin, Candidate of Technics, Senior Research Fellow at the Dynamics and Motion Control Department. E-mail [email protected]
