Идентификация маневров, выполняемых двигателями малой тяги космического аппарата Текст научной статьи по специальности «Математика»

УДК 629.7

Г. К. Б о р о в и н, М. В. З а х в а т к и н, В. А. С т е п а н ь я н ц, А. Г. Т у ч и н, Д. А. Т у ч и н, В. С. Я р о ш е в с к и й

ИДЕНТИФИКАЦИЯ МАНЕВРОВ, ВЫПОЛНЯЕМЫХ ДВИГАТЕЛЯМИ МАЛОЙ ТЯГИ КОСМИЧЕСКОГО АППАРАТА

Рассмотрена задача идентификации параметров маневра космического аппарата на основе обработки данных оптических измерений. Сложность этой задачи определяется большой длительностью работы двигателей и создаваемыми малыми ускорениями. Решена задача разделения интервала наблюдений на участки «пассивного» движения и определения параметров орбиты на них; оценены интервалы включения и выключения двигателя малой тяги; уточнены создаваемые двигателем ускорения. Предложенный алгоритм решения задачи используется при обработке результатов измерений Научной сети оптических инструментов для астрометрических и фотометрических наблюдений.

E-mail: bc@kiam1.rssi.ru

Ключевые слова: двигатель малой тяги, оптические измерения, орбитальные параметры, идентификация и определение маневров.

Идентификация параметров маневров космического аппарата (КА), выполняемых двигателями малой тяги, на основе обработки результатов оптических измерений является сложной задачей ввиду большой длительности работы двигателей КА и малых ускорений, создаваемых этими двигателями. Например, КА, выполняющий на геостационарной орбите (ГСО) маневр для изменения точки стояния, должен получить характеристическую скорость 2...3 м/с. Такая скорость достигается в результате работы двигателя малой тяги в течение одних или нескольких суток. Для удержания КА на ГСО в окрестности заданной точки стояния проводятся коррекции с затратами характеристической скорости 0,1.0,2 м/с. Двигатель малой тяги для выполнения таких коррекций включается на один или несколько часов.

При идентификации параметров маневров необходимо выполнить:

• разделение интервала наблюдений на участки «пассивного» движения и определение параметров орбиты на этих участках;

• оценку интервалов включения и выключения двигателей малой тяги;

• уточнение ускорений, создаваемых двигательной установкой.

Для определения параметров движения КА используют:

1) метод наименьших квадратов;

2) расширенный фильтр Калмана.

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

Расширенный фильтр Калмана — эмпирическое распространение фильтра Калмана для линейных систем на нелинейный случай. Метод предполагает наличие неизвестного шума, воздействующего на систему. Характеристики шума задаются его ковариационной матрицей. Особенностью этого метода является то, что текущие значения оцениваемых параметров должны находиться в достаточно близкой окрестности относительно их истинных значений. Кроме того, требуется, чтобы не было длительных интервалов времени, в которых нет измерений. Указанные обстоятельства не позволяют использовать расширенный фильтр Калмана для решения задачи разделения интервала наблюдений на участки «пассивного» движения и определения параметров орбиты на этих участках.

Для разделения интервала наблюдений на участки «пассивного» движения и определения параметров орбиты на этих участках предлагается применить метод, изложенный в работах [1—3]. Суть метода состоит в том, что оценку параметров орбиты КА получают из условия минимизации функционала, зависящего как от невязок расчетных и измеренных значений, так и от определяемых возмущений. Минимизация функционала выполняется итерационно. На каждом шаге итерационного процесса определяется поправка к искомым параметрам исходя из условия минимума функционала для линейной системы. Поиск минимума функционала для линейной системы приводит к двум последовательностям рекуррентных формул. Одна из них продвигается от первого измерения к последнему и позволяет определить поправку на момент последнего измерения; она эквивалентна рекуррентным формулам фильтра Калмана для линейной системы. Другая последовательность, которая идет от последнего значения к первому и позволяет восстановить возмущения, называется сглаживанием [4].

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

Основная идея предлагаемого алгоритма состоит в использовании результатов сглаживания для определения возможных моментов возмущений.

Пусть движение КА описывается дифференциальным уравнением

Л

^ = F х) + В (t), (1)

где х — вектор состояния; F (7, х) — вектор-функция; В (7) — матрица, описывающая воздействие шума на систему; \ (7) — белый шум с нулевым математическим ожиданием и заданной матрицей интенсивности Р (7) порядка 3 х 3.

Начальные условия для системы (1) задаются априорным вектором х0 и его ковариационной матрицей Р0.

В моменты времени t2,..., проводятся измерения функций Т1, Т2,..., Тм, значение которых обозначим через ¥гнаб . Для каждого момента времени ti справедливо

*/наб = Ъ , X (■)) + %, (2)

где % — случайный вектор, имеющий нулевое математическое ожидание и ковариационную матрицу .

В линейном случае (2) приобретает вид

гг = Нг )х) + ц,, (3)

где г, — вектор параметров размерности г, измеряемых в момент времени Н, {ti) — матрица размера г, х п.

Запись в качестве параметра х функции означает, что функция зависит не от мгновенного значения вектора состояния, а от функции х (V), которая является решением уравнения (1).

Решение уравнения (1) на интервале tN] будем аппроксимировать функциями вида

хЛ О) = хБ 0) + хР ^X

при этом

dx

D

dt

F (t, x D);

•n

xp (t) = Ф(t, tN)xP (tN)-J Ф(t,r)B(r)dT.

t

Матричная функция Ф (t, т) удовлетворяет уравнению dФ (t, tN )_dF

dt

dx

=x D (t)

Ф(t, ^N )

при начальных условиях Ф , tN ) = Е; хр ) — вектор, параметризующий семейство функций хр (t).

Покажем, что хА ) аппроксимирует решения системы (1). Дей-

ствительно,

dx. w ч öF

~t=F (t,x D )+äF

x=x D (t)

Ф (t, tN )x p (tN ) +

+B(t W) -

d¥

dx

x=xD (t) t

N

J Ф (t, r)B(r)£,(r)dz.

Сгруппировав члены, содержащие

отношение (4), получим dx

sf

dx

и использовав со-

X=XD (t)

w , dF dt F (t, x"

x=xd (t)

x p (t)+B(t) ^(t).

ÖF

Поскольку F (t, xA ) = F (t, xD + xp ) « F (t, xD ) н--

Sx

x=x d (t)

то ха (t) аппроксимирует решение уравнения (1).

Отметим, что хр (7) удовлетворяет линейному стохастическому дифференциальному уравнению

dx p _ dF dt c)x

x=x d (t)

xp (t) + B(t)£(t).

Зависимость между векторами состояния хА (ч,) в дискретные моменты времени t0, Ч1,..., tN можно выразить разностным уравнением

XА (Ч+1) = XВ (Ч,+1) + Ф (Ч+1, Ч ) (хА (Ч ) - Xд (Ч )) +

Чг+1

| Ф(^.+1,г)Б(г)§(г) йт.

+

г+1

Пусть у, = | Ф(чг+1, г)Б(г)^(г) йг. Случайный вектор у, имеет

Ч

нулевое математическое ожидание и ковариационную матрицу

О, = | Ф(^.+1,г)Б(Г)0(Г)Бт(г)Фт (чг+1,г)йг.

ч

Положим Xр (tN ) = 0. Тогда вектор начальных условий хВ (tN ) = хА (tN ) = х (tN ) однозначно определяет значения вектор-функции хА [ч) в дискретных точках Ч0, ..., tN. Однако при вычислении значений функции (7г-, хА требуется знать зависимость ха (ч) в окрестности каждого момента времени Ч,:

X A (0 =

хв ) = XN, если Ч = tN;

хв (Ч) + Ф(Ч Ч)[ха (Ч)-хв (Ч)] + У-1, если Ч <Ч < Ч^ хв (Ч0 ) + Ф (Ч Ч0 )[хА (Ч0 )- хв (Ч0 )], еСЛИ Ч0 < Ч < Ч1.

Таким образом, построена параметрическая зависимость хА (ч, q), где q — вектор уточняемых параметров, состоящий из компонент векторов х ^ ), у 0,..., у N-V

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

J =

2

(XA (^0, q) - X0 Г P0_1 (XA (?0, q) - X0 )

1 jN / \ т / \

+2£(*»наб - , XA (i, q))) R^ (TiHBß - (ti, XA (f, q))) +

1 N-1

+21 ■

2 i=0

Если минимум функционала (5) искать методом Ньютона, то поправки каждого шага итерации минимизируют квадратичную форму, полученную из (5) заменой нелинейных зависимостей линейными членами ряда Тейлора. Это означает, что на шаге итерации 5 решается задача оптимальной оценки состояния линейной системы (1),

(3). Матрицы Л^5) и И^^), I = 1,..., N, этой системы вычисляют по следующим формулам:

v J 8x

_1))!

H

(s

dx

a (t ,q(s-1))'

Значения переменных, для которых достигается минимум функционала

J (X0,N, X1,N, ■■■, XN,N, V0,N, ■■■, VN-1,N) = T(X0,N _ X0 ) P0 (X0,N _ X0 ) +

N-1

i=0 z L

H(Jw j Ry.,.1 (z, i+1 Hi+1Xi+1,N ) + Qi v

i+1 "i+1A i+1, N

i, N

при ограничениях

0(s)X,N + v i n , i = 0,1,..., N -1,

N+1,N ~ wz Ai,N 1 t/,N

удовлетворяет системе уравнении

X0,N _ Р0Ф0^0 + X0;

4+1,N - Ф/Xi,N + Qi h;

v i,N = Qi ь i;

^i = Ф L^iX, + H

¿+1 г+i

z i H i 3х

i+1Ai+1,N 5

i = 0,..., N -1;

Эта система уравнений распадается на две системы рекуррентных соотношений. Одна из них связывает априорную оценку с конечной поправкой на момент последнего измерения. Происходит продвижение по измерениям вперед. Эта рекуррентная система совпадает с алгоритмом фильтра Калмана:

хN,N = ФN-1хN-1,N-1 + р,NНN[2N ~ НNФN XN-1,N-1 ].

Другая рекуррентная система продвигается от последнего измерения к первому и осуществляет сглаживание. Она позволяет определить оценки векторов состояния (поправок) на момент времени Ч, по информации на интервале от Ч1 до N а также восстановить значение векторов шума:

Рассмотрим работу алгоритма на примере обработки результатов оптических измерений, получаемых с помощью КА 2000-031A (Express 3A).

Научная сеть оптических инструментов для астрометрических и фотометрических наблюдений [5] провела наблюдения трассы орбиты указанного КА в интервале от 18 июня до 29 сентября 2009 г.

На интервале измерений КА выполнил ряд маневров, в результате которых сначала был совершен перелет из точки стояния с долготой 349° в точку стояния с долготой 38°, а затем КА был переведен на орбиту с высотой более 400 км относительно высоты ГСО.

Трасса орбиты КА показана на рис. 1. Зависимость долготы подспутниковой точки от времени приведена на рис. 2 (по оси абсцисс — время в сутках от 0 ч 18 июня 2009 г, по оси ординат — долгота подспутниковой точки в градусах).

X i = Ф

v i,N = Qi i;

Рис. 1. Трасса орбиты КА Express 3A

О 10 20 30 40 50 60 70 80 90 100

Рис. 2. Зависимость долготы подспутниковой точки КА от времени

На первом этапе работы алгоритма были обнаружены моменты выполнения маневров по шкале иТС, указанные в табл. 1.

Таблица 1

Моменты включения двигателей малой тяги

Дата UTC Время UTC Дата UTC Время UTC

2009/06/19 05:08:07.000 2009/08/22 06:41:51.490

2009/08/03 08:20:28.030 2009/08/24 02:17:05.470

2009/08/08 20:44:06.530 2009/08/25 03:57:02.040

2009/08/15 19:49:34.990 2009/08/27 03:59:12.710

2009/08/17 19:54:41.340

Рассмотрим отношение невязки между расчетным и измеренным значениями к априорному СКО ошибки измерения. Это отношение

является безразмерной величиной и при ошибках измерений, соответствующих априорным значениям СКО ошибок, и правильных моделях движения не должно превышать 3. Будем называть это отношение относительной невязкой. Зависимости относительных невязок прямого восхождения и склонения от времени показаны на рис. 3 (по оси абсцисс — время в сутках, по оси ординат — относительная невязка).

Рис. 3. Зависимость относительных невязок по прямому восхождению и склонению от времени

Далее моменты выполнения импульсов и значения импульсов уточнялись по соответствующим интервалам измерений. При этом одновременно определялись вектор состояния на начало интервала измерений, время маневра и значение импульса. Результаты представлены в табл. 2: wMR.SE — это СКО относительных невязок; шБУ — уточненное значение модуля импульса; vR и уК — проекции вектора уточненного импульса на радиальное и трансверсальное направления соответственно; уВ — проекции вектора уточненного импульса на направление, ортогональное к плоскости орбиты.

Параметры, приведенные в табл. 2, позволяют оценить достоверность работы алгоритма. Оценку импульса на 2009/08/04 05:01:38 иТС нельзя назвать удовлетворительной, так как СКО относительных невязок превышает значение 3 и найденный импульс имеет большую радиальную составляющую. Поэтому об этом маневре можно говорить только как об обнаруженном маневре. Оценка параметров этого маневра недостоверна. Маневр 2009/08/18 0:19:18 имеет низкую характеристическую скорость. По-видимому, этого маневра не было. Маневр 2009/08/27 3:32:36 имеет высокую радиальную составляющую, что указывает на неправильность определения времени маневра.

Таблица 2

Результаты уточнения параметров маневров

Дата UTC Время wMRSE mDV, м/с vR, м/с vN, м/с vB, м/с

2009/06/19 09:27:14 1,80 1,280 -0,647 -1,069 -0,279

2009/08/04 05:01:38 5,26 10,915 10,373 3,317 -0,719

2009/08/15 20:17:17 0,25 0,154 0,154 -0,001 0,005

2009/08/18 00:19:18 0,29 0,019 0,019 0 0

2009/08/20 08:43:07 0,30 1,812 -0,358 1,715 -0,462

2009/08/22 08:48:06 0,26 1,579 0,487 1,484 -0,231

2009/08/23 08:44:16 0,22 1,419 -0,051 1,399 -0,230

2009/08/24 12:04:46 0,22 1,711 -0,845 1,447 0,347

2009/08/27 03:32:36 0,35 3,810 -3,807 0,059 -0,124

2009/09/02 10:09:57 0,39 71,085 -5,313 -1,163 70,877

Таким образом, в существующем состоянии из 10 найденных фактов маневров достоверными оказались девять. Из девяти оценок маневров достоверны семь.

Вывод. Пример обработки данных КА Express 3A показал эффективность предложенной схемы и разработанных методов обработки оптических измерений для определения параметров маневров и орбиты движения КА.

СПИСОК ЛИТЕРАТУРЫ

1. Kamensky S., Tuchin A., Stepanyants V., Alfriend K. T. Algorithm of automatic detection and analysis of non-evolutionary changes in orbital motion of geocentric objects // AAS/AIAA Astrodinamics Specialist Conf. Paper AAS 09-103.

2. Тучин А. Г. Определение параметров движения КА по результатам измерений при наличии шума в динамической системе. Препр. / ИПМ им. М.В. Келдыша РАН. - 2004. - № 2.

3. Аким Э. Л., Астахов А. П., Бакитько Р. В. и др. Автономная навигационная система околоземного космического аппарата // Изв. РАН. ТиСУ. -2009. - № 2.

4. Haupt G., Kasdin N., Keiser G., Parkinson B. An optimal recursive iterative algorithm for discrete nonlinear least-square estimation // AIAA-95-3218. AIAA Guidance, navigation, and control conf. Aug. 7-10, 1995. Baltimore.

5. Молотов И. E., Агапов В. М., Куприянов В. В. и др. Научная сеть оптических инструментов для астрометрических и фотометрических наблюдений // http://lfvn.astronomer.ru/report/0000042/index.htm.

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