Определение параметров движения МКС по данным навигационной системы GPS
( Determination of motion parameters of ISS by use of GPS measurements
Preprint, Inst. Appl. Math., the Russian Academy of Science)

1. GPS-приемник на Российском сегменте МКС. Для выполнения высокоточного навигационного обеспечения полета МКС на ее Российском и Американском сегментах установлены несколько GPS-приемников. К настоящему времени эти приемники успешно функционируют, и полученные с их помощью GPS-измерения фазового вектора МКС в составе телеметрической информации (ТМИ) поступают на Землю.

На Российском сегменте (РС) находится аппаратура спутниковой навигации АСН-2401, опытная эксплуатация которой проводилась в 2001 – 2002 гг. и возобновилась в июне 2003 г. Аппаратура размещается внутри и вне герметического объема служебного модуля РС. Антенны установлены так, чтобы свести к минимуму их затенение солнечными батареями и другими элементами конструкции МКС при полете в орбитальной ориентации – предполагаемой основной ориентации МКС на этапе ее штатной эксплуатации. Аппаратура АСН-2401 входит в состав системы управления движением РС и управляется командами Земли через бортовой компьютер.

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

2. Данные GPS-измерений и метод их обработки. Данные определений фазового вектора МКС, получаемые с помощью аппаратуры АСН-2401, представляют собой последовательности значений радиуса-вектора и скорости станции в гринвичской системе координат (точнее, в системе WGS84): , . Здесь , , ,  – время. Компоненты этих векторов указываются в единые моменты времени с шагом 1 с в сеансах непосредственной передачи ТМИ на Землю и шагом несколько минут () в сеансах передачи ТМИ из бортовых запоминающих устройств. Ниже используются данные только второго типа. Поскольку точность определения скорости станции с помощью аппаратуры АСН невысока (см. ниже), сглаживались только измерения радиуса-вектора.

Применявшаяся математическая модель орбитального движения станции представляла собой дифференциальные уравнения движения ее центра масс, записанные в гринвичской системе координат с учетом нецентральности гравитационного поля Земли и сопротивления атмосферы. Нецентральность поля учитывалась с точностью до членов порядка (36,36) включительно в разложении гравитационного потенциала Земли в ряд по шаровым функциям. Атмосфера считалась неподвижной относительно земной поверхности, плотность атмосферы рассчитывалась согласно модели ГОСТ 25645.115-84 (редакция 1990 г.).

Эта модель содержит три параметра: среднесуточное значение  индекса солнечной активности , среднее взвешенное значение  индекса  за предшествующие 81 сутки и среднесуточный индекс геомагнитной активности . Значения этих параметров с определенным запаздыванием относятся к суткам, для которых вычисляется искомая плотность. В общем случае геомагнитная активность меняется быстрее солнечной, поэтому в расчетах повышенной точности вместо индекса  используют трехчасовой индекс .

При сглаживании данных измерений GPS использовались две системы уравнений движения станции, отличающиеся способами задания индексов ,  и  в формулах ГОСТа. Эти системы обозначим (I) и (II). В системе (I) индексы ,  и  оставались неизменными на всем интервале сглаживания. В системе (II) учитывалась зависимость этих индексов от времени, причем функции ,  и  выбирались кусочно-постоянными. Интервалы постоянства функций  и  имели длину 1 сут, интервалы постоянства функции  – 3 ч. Граничные точки интервалов постоянства определялись условиями ГОСТ, значения функций принимались равными значениям индексов, взятым в Internet. При сглаживании данных измерений на коротких интервалах времени (менее трех витков) значение баллистического коэффициента  станции было фиксировано. На длинных интервалах времени (более трех витков) этот коэффициент служил параметром согласования.

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

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

,         .

где  – вычисленная в точке  матрица системы нормальных уравнений, возникающей при минимизации  методом Гаусса-Ньютона, . Точность аппроксимации данных измерений найденным решением уравнений движения и точность определения этого решения будем характеризовать стандартным отклонением ошибок измерений , стандартными отклонениями   компонент вектора  и некоторыми статистическими характеристиками рядов остатков

,        .

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

При минимизации  удобно варьировать не баллистический коэффициент , а безразмерное отношение , где  – фиксированное номинальное значение этого коэффициента. Ниже под  всюду понимается стандартное отклонение параметра .

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

,    ,    .

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

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

На первом этапе при значении баллистического коэффициента, взятого из данных радиоконтроля орбиты, функционал  минимизировался по начальным условиям . Для найденного решения вычислялась медиана  модулей остатков . Из обработки исключались точки , для которых . При сделанных выше предположениях относительно ошибок измерений , и вероятность появления остатка с указанным значением модуля – менее 0.001 (модули остатков имеют распределение Максвелла с параметром ). На втором этапе минимизация  по  и  повторялась для модифицированного набора измерений. Полученные оценки начальных условий считались окончательными, если интервал с измерениями был коротким. В случае продолжительного интервала проводились еще два этапа. Третий этап состоял в вычислении оценки  по измерениям, допущенным ко второму этапу. Затем проводилась описанная выше отбраковка измерений, причем проверялись и ранее отвергнутые точки. Для пересмотренного набора значений  находилась новая оценка , считавшаяся окончательной.

3. Результаты сглаживания данных GPS-измерений. Ниже описываются результаты сглаживания данных измерений, полученных в основном на отрезке времени с 25.V по 30.V.2002 г. Графики изменения индексов , ,  и  на этом отрезке приведены на рис. 1. В системе (I) для указанного промежутка времени принималось ,  (). Это – данные, использовавшиеся ЦУП. В обеих системах принималось  м/кг. Такое значение баллистического коэффициента также использовалось ЦУП.

Выбранные для сглаживания данные измерений были разнесены по 10 временным интервалам, основные характеристики которых приведены в табл. 1. Первые 9 интервалов относятся к указанному выше промежутку времени, десятый интервал – более поздний. На этом интервале в системе (I) ,  (), м/кг. Все эти интервалы считались продолжительными и обрабатывались в четыре этапа. В таблице:  – число точек включенных в обработку на первом этапе;  – число точек, включенных в обработку на заключительном этапе;  и  – начальная и конечная точки исходного измерительного интервала (здесь использовано декретное московское время – ДМВ). Для удобства во всех случаях начальные условия задавались в момент . Точка  на всех этапах обработки не менялась, несмотря на случавшуюся иногда отбраковку измерения . В некоторых клетках таблицы для числа  указаны два значения, первое отвечает системе (I), второе – системе (II). В случаях совпадения этого числа для обеих систем в таблице указано одно значение. Интервалы 6, 7, 8 и 9 наполнены измерениями из предыдущих интервалов таблицы с исключением части измерений, признанных ошибочными.

Таблица 1. Обработанные интервалы

Результаты сглаживания данных измерений на интервалах табл. 1 решениями систем (I) и (II) приведены в табл. 2 – 7 и на рис. 2 – 10. Табл. 2 – 4 и рис. 2 – 8 получены с использованием системы (I), табл. 5 – 7 и рис. 9, 10 – с использованием системы (II). Табл. 2 и 5 содержит значения параметра , стандартные отклонения ошибок в исходных данных и стандартные отклонения определяемых параметров сглаживающих решений – начальных условий и . В табл. 3, 4, 6, 7 приведены указанные выше характеристики рядов остатков  для . На рис. 2 – 10 изображены диаграммы этих рядов, а также диаграммы модулей остатков  и . Диаграммы представляют собой ломаные с вершинами в точках ,  и . Для большей наглядности вершины указаны маркерами.

Полученные результаты относятся к интервалам слабо возмущенного движения станции. В настоящее время такие интервалы непродолжительны. В частности, между интервалами 2 и 3 табл. 1 лежит интервал с сильными возмущениями. Из-за малой длины интервалов слабо возмущенного движения GPS-измерения пока не удается в полной мере использовать для верификации моделей влияния сопротивления атмосферы на движение станции.

Анализ результатов обработки интервалов 1 – 9 подтверждает вывод [1] о том, что на сравнительно коротких (менее 3 сут) интервалах времени система (II) не имеет преимуществ перед системой (I). В данном случае система (I) также оказалась немного точнее. Но такой вывод справедлив «в среднем» – при отсутствии значительной геомагнитной активности. В данном случае эта активность была невелика (рис. 1). К тому же, индекс  в рассматриваемый период времени оставался стабильным. Как следствие, значения параметра  от интервала к интервалу менялись мало (см. табл. 2, 5).

Таблица 2. Параметр , стандартные отклонения ошибок данных измерений

и параметров сглаживающих решений системы (I).

Таблица 3. Средние квадратические значения рядов остатков

сглаживающих решений системы (I).

          Точность определения движения станции по данным GPS-измерений соответствует паспортной точности измерения координат аппаратурой АСН-2401 (50 м). Характеристики точности аппроксимации данных измерений координат – величины , ,  и  – в своем большинстве не превышают геометрических размеров станции. Вместе с тем, погрешность аппроксимации несколько возрастает с увеличением длины интервала сглаживания, что свидетельствует о некоторой неадекватности используемых уравнений движения.

Таблица 4. Медианы рядов модулей остатков сглаживающих

решений системы (I).

Таблица 5. Параметр , стандартные отклонения ошибок данных измерений

и параметров сглаживающих решений системы (II).

Следует отметить устойчивый вид систематических ошибок в измерениях скорости станции. В остатках  (невязки аппроксимации составляющей абсолютной скорости по радиусу-вектору) имеется постоянная составляющая около 8 м/с. В остатках  (невязки по нормали к плоскости орбиты) доминирует периодическая составляющая с орбитальным периодом и амплитудой примерно 0.8 м/с. Если эти систематические составляющие устранить, то модифицированные остатки  и , по-видимому, будут выглядеть примерно также как остатки  (невязки вдоль вектора абсолютной скорости).

Таблица 6. Средние квадратические значения рядов остатков

сглаживающих решений системы (II).

Таблица 7. Медианы рядов модулей остатков сглаживающих

решений системы (II).

          4. Прогноз движения МКС по данным GPS-измерений. Данные измерений на интервалах 6 – 8 табл. 1 позволяют оценить возможную точность такого прогноза на отрезке времени длиной около 1 сут. Каждый из указанных интервалов объединяет два сеанса измерений, разнесенных примерно на сутки. В такой ситуации интересно оценить погрешность аппроксимации данных измерений второго, более позднего по времени сеанса решением системы (I) или (II), найденным по измерениям на первом сеансе. Эту аппроксимацию назовем прогнозом 1. Его результаты для указанных интервалов представлены на рис. 11 – 16. По виду новые рисунки в точности аналогичны рис. 2 – 10, но интерпретировать их следует по-другому.

Рассмотрим, например, рис. 11. На этом рисунке изображены диаграммы остатков решения системы (I), построенного по измерениям на интервале 1 и продолженного на интервал 6, т.е. до конца интервала 2. Начальные отрезки этих диаграмм совпадают с диаграммами на рис. 2. Изменено только начало отсчета времени из-за того, что первые два измерения интервала 1 были признаны ошибочными и в интервал 6 не попали (ср. значения  для интервалов 1 и 6 в табл. 1). Заключительные части диаграмм на рис. 11 характеризуют погрешность аппроксимации данных измерений интервала 2 решением, построенным для интервала 1. Значения остатков в этой заключительной части существенно больше, чем в случае аппроксимации тех же данных решениями системы (I), построенными специально для интервала 2 (рис. 3) и для интервала 6 (рис. 5). Сравнение рис. 11 с рис. 3 и 5 дает наглядное представление о том, к какой точности прогноза в данном случае следовало бы стремиться, т. е. характеризует потерю точности при прогнозе.

          Аналогичным образом сравнение рис. 6 и 13 дает представление о погрешности прогноза на интервал 5 движения станции, построенного по измерениям на интервале 4; сравнение рис. 10 и 15 дает представление о прогнозе на интервал 4 движения станции, построенного по измерениям на интервале 3. Первый из этих прогнозов построен с использованием системы (I), второй – с использованием системы (II). Анализ рис. 11 – 16 показывает, что для целей прогноза система (II) в целом оказалась более точной, чем система (I). Отметим, что во всех вариантах за исключением варианта прогноза, представленного на рис. 12, ошибка прогноза положения станции на 15 ч =  с вперед не превышает 400 м. Для варианта на рис. 12 эта ошибка не превышает 600 м.

          Сделанные оценки точности прогноза относились к рафинированной ситуации. Сеансы измерений, по которым строился прогноз, были достаточно продолжительными, и на них удалось выполнить уточнение баллистического коэффициента станции. Рассмотрим, что может дать прогноз без уточнения этого коэффициента. Выделим на интервалах 6, 7 и 8 начальные отрезки данных, охватывающие примерно по 3 ч. Эти отрезки содержат соответственно 20, 35 и 18 точек . Построим по этим отрезкам сглаживающие решения систем (I) и (II) и оценим точность аппроксимации такими решениями данных на полных интервалах. Баллистический коэффициент при построении сглаживающего решения на коротком отрезке и последующем прогнозе будем считать неизменным. Значения этого коэффициента на интервале 6 (соответственно 7 или 8) возьмем из результатов обработки данных измерений на интервале 1 (соответственно 3 или 4). Описанную процедуру назовем прогнозом 2. Результаты такого прогноза представлены на рис. 17 – 22, которые аналогичны рис. 11 – 16. Анализ результатов, показывает, что во всех рассмотренных вариантах ошибка прогноза положения станции на 15 ч вперед не превышает 400 м. Правда, точки с абсциссой  с лежат вне областей с измерениями, но по виду диаграмм остатков (сделав их “интерполяцию”) можно судить, что указанная оценка в действительности выполняется.

Узким местом оценок, полученных в рамках прогноза 2, является выбор баллистического коэффициента. Здесь его значения выбраны наилучшим образом. С другой стороны, ясно, что совершив ошибку при выборе этого коэффициента, можно получить неприемлемо большую ошибку прогноза. По-видимому, для построения надежного точного прогноза движения станции по GPS-измерениям следует предусмотреть достаточно продолжительный (до 10 ч) интервал сбора этих измерений с целью уточнения баллистического коэффициента.

Прогноз 2 был построен и для интервала 10 (рис. 23). В этом случае баллистический коэффициент был взят из результатов обработки данных на всем интервале. Для построения прогноза использовалась система (I) и 61 точка , попавшая в трех часовой начальный отрезок. Прогноз оказался приемлемым, что не удивительно, поскольку он выполнялся по 75 % измерений.

5. Заключение. Точность определения слабо возмущенного движения станции по данным GPS-измерений ее гринвичских координат соответствует паспортным данным аппаратуры АСН-2401. Средние квадратические значения ошибок аппроксимации этих данных на интервалах времени менее 1 сут не превышают 30 м.

Данные GPS-измерений скорости станции содержат устойчивые систематические ошибки. Ошибки рассчитанных по измеренному фазовому вектору компонент абсолютной скорости в естественном трехграннике имеют следующий вид. Компонента по нормали к траектории (эта нормаль практически коллинеарна радиусу-вектору) содержит постоянную ошибку около 8 м/с. Компонента по бинормали к траектории (нормали к плоскости орбиты) содержит периодическую ошибку с орбитальным периодом и амплитудой примерно 0.8 м/с. Для выяснения природы этих ошибок необходим анализ выполняемых аппаратурой АСН-2401 измерений псевдодальности и псевдоскорости.

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

Литература

1. Беляев М.Ю., Медведев Е.С., Рулев Д.Н., Сазонов В.В. Влияние учета вариаций плотности атмосферы, вызванных солнечной активностью, на точность определения движения МКС. Препринт ИПМ им. М.В. Келдыша РАН, №    , 2004.

Приложение