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

( The influence of atmosphere density variations, caused by geomagnetic activity, on accuracy of the forecast of ISS motion
Preprint, Inst. Appl. Math., the Russian Academy of Science)

          1. Постановка задачи. Оценивается точность прогноза неуправляемого движения МКС на интервалах времени длиной до 15 ч. Прогноз считается приемлемым, если модуль ошибки предсказания радиуса-вектора станции не превышает 400 м. Из всех возможных источников ошибок прогноза учитывается только неточность знания плотности атмосферы. Оценка находится методом статистического моделирования с использованием реальных данных о состоянии атмосферы и модели плотности верхней атмосферы, определенной ГОСТ 25645.115-84 (редакция 1990 г.).

Эта модель содержит три параметра: среднесуточное значение  индекса солнечной активности , среднее взвешенное значение  индекса  за предшествующие 81 сутки и среднесуточный индекс геомагнитной активности . Значения параметров с некоторым запаздыванием относятся к суткам, на которые вычисляется искомая плотность. Характерное время изменения индекса  – несколько суток, и это изменение достаточно точно прогнозируется. При статистическом моделировании значение  считалось постоянным. Индекс  является огрубленной характеристикой. Детальное изменение геомагнитной активности во времени описывается трехчасовым индексом , значения которого служат исходной информацией при вычислении . В расчетах повышенной точности следует использовать именно их. Статистическое моделирование использует значения , относящиеся к 1998 – 2002 гг.

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

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

.                               (1)

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

          Решения системы (1) рассматривались на фиксированном отрезке времени ,  ч, при одинаковых начальных условиях и параметрах , . Решения отличались только использованными для их построения функциями . Эти функции выбирались кусочно-постоянными с отрезками постоянства , где  ч, ,  – минимальное натуральное число, удовлетворяющее условию . При  принималось . Исследование точности прогноза движения станции состояло в статистическом исследовании разброса значений  рассматриваемых решений на реалистическом множестве изменения параметров .

Исследование разброса величин  сводилось к исследованию разброса разностей . Здесь  – опорное решение системы (1), построенное при постоянном значении  и том же начальном условии, что и решения рассматриваемого семейства. Как оказалось, в представляющей интерес области изменения величин  и  с относительной погрешностью около 0.1% выполнено равенство , где  – решение линейной начальной задачи

,    .     (2)

Функция  является кусочно-постоянной: при  имеем . Введем функции , положив  при  и  в остальных случаях . Поскольку

,

решение начальной задачи (2) и вектор  можно представить в виде

,      ,                             (3)

где  – решения начальных задач

,     .

Величины  вычислялись однократно перед выполнением статистического моделирования.

В качестве значений  при моделировании использовались данные измерений трехчасового индекса  за 5 лет – с 1 января 1998 г. по 31 декабря 2002 г. Примеры использованных данных приведены на рис. 1, 2. Здесь изображены графики изменения во времени индекса  и отвечающего ему также трехчасового индекса . Последний имеет ясный физический смысл и более наглядно характеризует геомагнитную активность. Графики задают кусочно-постоянные функции с интервалами постоянства 3 ч. Для удобства восприятия на рисунках указаны вертикальные отрезки, соединяющие предельные значения функций в точках разрыва. Графики на рис. 1 показывают изменение индексов ,  в течение первых 200 суток рассматриваемого интервала. Левые графики на рис. 2 дают пример максимально возможного возмущения с , . Правые графики на этом рисунке иллюстрирует продолжительный отрезок времени слабой геомагнитной возмущенности.

Используемые данные измерений  были упорядочены во времени. Из этих данных выбиралась последовательность идущих подряд значений

  ,                                  (4)

где  – натуральное число. Первые  членов этой последовательности пересчитывались в значения индекса  и усреднялись. Результат усреднения пересчитывался назад в  и обозначался . Вычисления выполнялись по схеме:

,    ,    .

Затем находилась разность , в которой  и  определены второй формулой (3) соответственно при значениях   из последовательности (4) и при одинаковых . Фактически и  не находились, а сразу по указанной формуле (3) при  вычислялась их разность. Введение величин и  удобно для ее интерпретации.

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

Векторы  вычислялись для всех возможных последовательностей (4). На первом шаге в качестве такой последовательности использовались первые  элементов имеющегося массива значений индекса . На втором шаге использовались идущие подряд  элементов этого массива, начиная со второго, на третьем –  элементов начиная с третьего и т. д. При  таких последовательностей из имеющего массива данных можно выбрать более 14500.

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

,    ,    .

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

Составим вектор . Среднее значение этого вектора и его ковариационную матрицу [1], вычисленные по всем реализациям последовательностей (4), обозначим соответственно  и . Квадратные корни из диагональных элементов матрицы  являются стандартными отклонениями соответствующих компонент вектора . Обозначим эти стандартные отклонения . Анализ результатов моделирования состоял в анализе величин ,  и т. п. в функции параметров  и , а также в изучении гистограмм распределений компонент вектора .

3. Результаты статистического моделирования. Первым опорным решением служило фактическое движение станции 09.03.2003 г. Его начальные условия были определены по данным измерений GPS, выполненных на Американском сегменте МКС. Опорное решение рассчитывалось при  , и 12:28:38 ДМВ. Параметры атмосферы и баллистический коэффициент (0.003365 м/кг) были взяты из данных ЦУП. Результаты моделирования для некоторых значений  и  представлены в табл. 1, 2 и на рис. 3, 4. В заголовках таблиц и подписях к рисункам рассматриваемый прогноз обозначен как прогноз 1. В таблицах здесь и ниже  – число испытаний в данном варианте,  – доля испытаний, в которых было выполнено неравенство  м.

Таблица 1. Прогноз 1 при  и  ч.

Таблица 2. Прогноз 1 при  и .

Из табл. 1 видно, что точность прогноза довольно слабо зависит от числа  измерений индекса , по которым вычисляется используемое в прогнозе значение . Хотя, как и следовало ожидать, с уменьшением  точность прогноза несколько повышается. Согласно табл. 2 ошибка прогноза быстро уменьшается с уменьшением . Согласно обеим таблицам основная ошибка в прогнозе положения станции совершается вдоль орта  – по направлению вектора абсолютной скорости в точке . Ошибка в направлении орта  (по радиусу-вектору) примерно в 50 раз меньше, ошибка в направлении орта  (по нормали к плоскости орбиты) пренебрежимо мала. Ошибки прогноза абсолютной геоцентрической скорости станции наиболее велики в направлении орта  и пренебрежимо малы в направлении . Эти факты находятся в полном согласии с теорией возмущенного движения ИСЗ в атмосфере [2, 3].

Распределение ошибок в пространстве  можно пояснить, рассмотрев собственные числа и собственные векторы ковариационной матрицы . Эти величины слабо зависят от , поэтому ограничимся анализом случая . Собственные числа матрицы  при  составляют (при расчете  единицей измерения длины служили метры, скорости – мм/с)

,   ,   ,   ,   ,   .

Нормированные собственные векторы, отвечающие максимальному и второму по величине собственным числам, соответственно имеют вид

,

.

Модули компонент вектора

,

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

Из приведенного анализа собственных чисел и векторов матрицы  следует, что ошибки прогноза в пространстве  распределены в основном вдоль направления . Модули коэффициентов корреляции лишь в случае пар ,  составляют 0.88. Модули коэффициентов корреляции остальных пар переменных превосходят 0.96. Проекции вектора  на подпространства  и  характеризуют распределение ошибок прогноза положения и скорости станции в базисе .

          Детальное представление о полученном при статистическом моделировании распределении ошибок прогноза по величине дают гистограммы [1] переменных , , , ,  и . Гистограммы строились единообразно, поэтому ограничимся описанием построения гистограммы координаты . Пусть  – натуральное число. Отрезок  разбивался на  отрезков равной длины, и для каждого отрезка разбиения подсчитывалось число оказавшихся в нем полученных при моделировании значений . Пусть в отрезок с номером  попало  точек. Гистограммой называется кусочно-постоянная функция , которая на отрезке с номером  принимает значение . Напомним,  – общее число испытаний в данном варианте моделирования (при данном ). Графики гистограмм , , , ,  и , построенные для варианта с  при , изображены на рис. 3. Эти графики выглядят как ступенчатые линии, поскольку для наглядности на рисунках указаны вертикальные отрезки, соединяющие предельные значения гистограмм в точках разрыва.

          Полученные гистограммы достаточно точно аппроксимируются плотностью вероятности двухстороннего показательного распределения

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

          Для сравнения на рис. 4 приведены результаты аппроксимации полученных гистограмм плотностью вероятности нормального распределения

.

В этом случае при аппроксимации  принималось ,  и т. д. Как показывает рисунок, такая аппроксимация является существенно менее точной.

          Рис. 3, 4 построены для варианта  в случае . Визуально примерно такое же качество аппроксимации достигается во всех вариантах из табл. 1 при . На рис. 5 приведены гистограммы и их аппроксимации плотностью вероятности  для варианта  в случае .

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

Таблица 3. Прогноз 1 при  и  ч.

Таблица 4. Прогноз при  и .

Результаты расчетов в случае  приведены в табл. 5, 6 и на рис. 7, 8. Таблицы аналогичны предыдущим. Они показывают, что погрешность прогноза заметно возросла по сравнению со случаем . На рис. 7, 8 изображены гистограммы , …, , построенные при  для варианта , и аппроксимации этих гистограмм плотностями вероятностей  (рис. 7) и  (рис. 8). Параметры ,  и  аппроксимирующих функций определялись выборочными средними значениями и дисперсиями найденных ошибок прогноза. По сравнению со случаем  аппроксимация гистограмм плотностью  стала точнее. Точность аппроксимации плотностью  уменьшилась, но в отличие от случая  это уменьшение точности вызвано увеличением разброса переменных , …, . Распределение ошибок прогноза приблизилось к нормальному.

Таблица 5. Прогноз 1 при  и  ч.

Таблица 6. Прогноз при  и .

Поведение матрицы , ее собственных чисел и векторов в случаях  и  в точности такое же, как в случае .

4. Прогноз при неизменном . Для более полной оценки полученных результатов сравним прогноз движения станции, выполняемый при выборе  по схеме раздела 3, с более простым прогнозом, который всякий раз рассчитывается при одном и том же значении этого индекса. Статистическое моделирование такого прогноза проводилось следующим образом. Все варианты рассчитывались при  ч . Задавалось значение . На первом шаге из имеющегося упорядоченного массива данных измерений трехчасового индекса  (см. раздел 2) выбирались первые пять элементов. Эта упорядоченная пятерка принимались в качестве параметров . Для таких значений  и  способом, описанным выше, рассчитывался вектор . На втором шаге в качестве параметров  использовались взятые подряд 5 элементов этого массива начиная со второго. Для прежнего значения  и новых  снова рассчитывался вектор . На третьем шаге из массива значений  выбирались 5 элементов начиная с третьего и т. д. Всего таких шагов выполнялось 14605. Использовались прежние опорные решения. Обработка результатов вычислений выполнялась по описанной выше схеме. В качестве значений  для определенности использовались значения 1,1+ (1.3333), (1.6667) и т.д. Результаты приведены лишь для некоторых значений  и обязательно содержат вариант с минимальным значением . В заголовках таблиц и подписях к рисункам новый прогноз обозначен как прогноз 2.

Результаты, полученные для опорного решения в случае , представлены в табл. 7 и на рис. 9. Таблица аналогична табл. 1, 3 и 5, только в ней вместо значений  указаны значения , и отсутствует столбец со значениями  – во всех вариантах моделирования теперь . Из таблицы видно, что изменение  приводит только к смещению среднего значения прогноза, а его стандартные отклонения остаются неизменными. Указанное смещение приводит к значительному увеличению . При  смещение весьма мало, и описанный простейший способ прогноза практически не уступает по точности, способу рассмотренному выше.

Матрица  во втором способе прогноза практически не зависит от . При  ее собственные числа составляют

,   ,   ,   ,   ,   ;

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

Модули компонент вектора

,

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

Таблица 7. Прогноз 2 при .

На рис. 9 приведены гистограммы значений переменных , , …, , полученных в варианте статистического моделирования с  (ступенчатые линии). Способ построения этих гистограмм описан в предыдущем разделе. В данном случае . Эти гистограммы достаточно точно аппроксимируются плотностью вероятности гамма распределения с показателем формы 2:

.

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

Результаты, полученные при , представлены в табл. 8 и на рис. 10. Эта таблица в точности аналогична табл. 7. Из нее видно, что погрешность прогноза заметно уменьшилась. На рис. 10 изображены гистограммы , …,  (ступенчатые линии), построенные при  для варианта , и аппроксимации этих гистограмм плотностью вероятности  (плавные линии) при описанном выше выборе параметров  и . Из рисунка видно, что в данном случае аппроксимация не так точна, как в случае . Распределение ошибок прогноза стало похоже на одностороннее показательное распределение. Используя термины гамма распределения можно сказать, что коэффициент формы заметно уменьшился.

Таблица 8. Прогноз 2 при .

Результаты расчетов в случае  приведены в табл. 9 и на рис. 11. Таблица аналогична табл. 4, 5. Из нее видно, что погрешность прогноза заметно возросла по сравнению со случаем .

На рис. 11 изображены гистограммы , …, , построенные при  для варианта , и аппроксимации этих гистограмм плотностью вероятности . Из рисунка видно, что и в данном случае аппроксимация не так точна, как в случае . Распределение ошибок прогноза стало похоже на гамма-распределение с коэффициентом формы, заметно превышающем 2.

Для второго и третьего опорных решений также имеет место независимость стандартных отклонений прогноза от . При высоком уровне солнечной активности прогноз 2 заметно уступает по точности прогнозу 1. При низком и среднем уровне солнечной активности оба прогноза обеспечивают в среднем примерно одинаковую точность.

Таблица 9. Прогноз 2 при .

5. Заключение. Статистическое моделирование показало, что влияние геомагнитной активности на точность предсказания положения станции существенно зависит от уровня солнечной активности. В случае длинных интервалов эта зависимость может оказаться критической. Например, в конце 15-часового интервала прогнозирования при малой солнечной активности вероятность неудачного прогноза составляет несколько процентов (4% при . В случае средней солнечной активности (при . эта вероятность возрастает до 14%. При высокой солнечной активности (при  вероятность неудачного прогноза превосходит 35%. С уменьшением длины интервала вероятность неудачного прогноза снижается. На интервалах с длиной менее 7.5 ч эта вероятность не превосходит 1%. На интервалах менее 4.5 ч неудачных прогнозов не оказалось.

Литература

1.     Ивченко Г.И., Медведев Ю.И. Математическая статистика. М., Высшая шклола, 1984.

2.     Охоцимский Д.Е., Сихарулидзе Ю.Г. Основы механики космического полета. М., Наука, 1990.

3.     Эльясберг П.Е. Введение в теорию полета искусственных спутников Земли. М., Наука, 1965.

Приложение