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

( The influence of taking into account atmosphere density variations, caused by solar activity, on accuracy of definition of ISS motion
Preprint, Inst. Appl. Math., the Russian Academy of Science)

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

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

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

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

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

          Указанные трудности не возникают, если ограничиться анализом движения станции на тех отрезках времени в 1999 г., когда ее полет был неуправляем и проходил в режиме гравитационной ориентации вращающегося спутника [1]. В этом режиме станция вращалась с угловой скоростью 0.2 – 0.3 град./с вокруг главной центральной оси минимального момента инерции, совершавшей малые колебания относительно местной вертикали. Длительность отрезков неуправляемого полета доходила до 20 сут. Баллистический коэффициент на этих отрезках оставался практически неизменным.

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

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

В системе (III) каждая из функций ,  и  аппроксимировалась суммой линейной функции и отрезка ряда Фурье. Например, при определении движения станции на отрезке  функция  аппроксимировалась выражением

,

где  и  – удовлетворяющие неравенствам ,  ближайшие к точкам  и  средние точки интервалов постоянства функции  в системе (II),  – число таких интервалов постоянства, имеющих общие точки с отрезком . Коэффициенты выражения  находились в два этапа. Сначала они определялись из условия равенства этого выражения значениям функции  в средних точках интервалов постоянства, пересекающихся с отрезком . Таких точек – , и задача первого этапа решается единственным образом. Затем коэффициенты , , …,  умножались на сглаживающие -множители Ланцоша [2]. Аналогичным образом строились аппроксимирующие выражения  и  для функций  и .

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

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

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

,       ,

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

,        .

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

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

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

,    ,    .

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

3. Результаты построения аппроксимации. На отрезке времени с 26.III по 14.IV.1999 г. полет станции проходил в режиме гравитационной ориентации вращающегося спутника. По информации КБ «Салют» на этом отрезке никакие управляющие силы и моменты на станцию не действовали. Графики изменения индексов , ,  и  на указанном отрезке времени приведены на рис. 1. Графики задают кусочно-постоянные функции. Интервалы постоянства для индексов  и  имеют длину 1 сут, для индексов  и  – 3ч. Для удобства восприятия на графиках указаны вертикальные отрезки, соединяющие предельные значения функций в точках разрыва. Графики соответствующих аппроксимирующих выражений , и  приведены на рис. 2. На этом рисунке нижние кривые – графики аппроксимирующих выражений, верхние кривые – графики тех же выражений и ломаные с вершинами в точках, по которым эти выражения были построены. Абсциссы точек – середины отрезков постоянства функций ,  и  на рис. 1, ординаты – соответствующие значения этих функций.

При определении движения станции на интервалах времени в период с 26.III по 14.IV.1999 г. в системах (II), (III) использовались зависимости индексов ,  и  от времени, представленные на рис. 1, 2. Точнее, построение аппроксимирующих выражений , и  для системы (III) при сглаживании с ее помощью данных измерений в точках  (1, 2,…, ; ) выполнялось по описанному выше алгоритму с , . Однако отличие получаемых при этом аппроксимирующих выражений от выражений на рис. 2 незначительно. В системе (I) для указанного промежутка времени принималось ,  (). Это – данные, принятые в ЦУП. Во всех трех системах принималось  м/кг. Это значение баллистического коэффициента также использовалось ЦУП.

По двухстрочным элементам, относящимся к промежутку времени 26.III – 14.IV.1999 г., была рассчитана последовательность значений фазового вектора станции из 67 членов. Эта последовательность была разбита на 5 непересекающихся отрезков, охватывающих от 2.3 до 3.7 сут. Полученные отрезки данных были объединены в более продолжительные отрезки. На некоторых отрезках с помощью решений систем (I), (II) и (III) выполнялось сглаживание значений фазового вектора. Характеристики использованных для сглаживания отрезков данных приведены в табл. 1. Здесь  – число измерений фазового вектора на данном отрезке;  и  – начальная и конечная точки отрезка по шкале времени UTC. Интервалы 1 – 5 таблицы отвечают упоминавшимся выше 5 исходным непересекающимся отрезкам данных. Для удобства во всех случаях начальные условия сглаживающих решений задавались в момент .

Таблица 1. Интервалы аппроксимации движения МКС.

Опишем сначала результаты сглаживания, полученные с использованием системы (I). Эти результаты приведены в табл. 2 – 4 и на рис. 3, 4. Табл. 2 содержит значения параметра , стандартные отклонения ошибок задания исходных данных и стандартные отклонения определяемых параметров сглаживающих решений – начальных условий и . В табл. 3, 4 приведены характеристики рядов остатков  для . На рис. 3, 4 изображены диаграммы этих рядов, полученных при аппроксимации интервалов 2 и 1Ú2Ú3 из табл. 1. Диаграммы представляют собой ломаные с вершинами в точках . Для большей наглядности вершины указаны маркерами.

Сопоставление табл. 2 с графиком функции  на рис. 1 показывает, что изменение параметра  отслеживает изменение индекса . Этот факт легко объясним. Поскольку в системе (I) индексы  и  приняты постоянными, фактическое увеличение индекса  приводит к увеличению . Однако изменением  не удается полностью скомпенсировать резкое изменение  внутри интервала сглаживания. Например, интервалы 2Ú3 и 2Ú3Ú4, содержащие примерно посередине участок резкого возрастания индекса  (см. рис. 1), сглажены заметно менее точно, чем имеющие примерно одинаковую с ними длину соседние интервалы табл. 1. Заметим, что соседние интервалы 3Ú4, 1Ú2Ú3, 3Ú4Ú5 содержат тот же самый участок, но с краю.

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

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

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

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

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

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

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

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

Результаты сглаживания, полученные с использованием системы (II), приведены в табл. 5 – 7 и на рис. 5, 6. Эти таблицы и рисунки аналогичны описанным выше табл. 2 – 4 и рис. 3, 4. Попарное сравнение табл. 2 и 5, 3 и 6, 4 и 7 показывает, что при переходе от системы (I) к системе (II) точность аппроксимации измерений на сравнительно коротких интервалах времени 1 – 5 не улучшилась. В большинстве случаев (3 из 5) она даже незначительно ухудшилась. Точность аппроксимации на более длинных интервалах повысилась. Особенно значительно это повышение в случае интервалов 2Ú3 и 2Ú3Ú4. Такое повышение обусловлено, прежде всего, учетом изменения параметров  и  в системе (II). Именно учет изменения этих параметров позволил устранить влияние резкого скачка функции  внутри интервалов 2Ú3 и 2Ú3Ú4. Сравнение рис. 5, 6 с рис. 3 и 4 и аналогичных рисунков, не вошедших в препринт, показывает уменьшение систематической составляющей в остатках  и  в случае аппроксимации решениями системы (II) данных измерений на продолжительных интервалах времени.

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

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

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

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

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

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

Результаты сглаживания, полученные с использованием системы (III), приведены в табл. 8 – 10. Эти таблицы аналогичны таблицам, рассматривавшимся выше. Анализ таблиц показывает, что переход от системы (I) к системе (III) также заметно повышает точность аппроксимации данных измерений на продолжительных интервалах, но в целом точность системы (III) ниже точности системы (II).

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

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

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

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

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

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

Проведенное исследование показало, что на интервалах времени более 5 сут учет изменения индекса  позволяет заметно увеличить точность аппроксимации данных измерений. Специфическое поведение индекса солнечной активности в период с 26.III по 14.IV.1999 г. дает возможность проиллюстрировать это и ряд других сделанных утверждений следующим более ярким примером, чем приведенные в табл. 1 – 10. Возьмем часть рассматривавшегося выше исходного отрезка данных измерений, у которой ,  01.IV,  10.IV (использована шкала времени UTC), . Новый отрезок данных назовем интервалом 6. Характеристики точности аппроксимации этого интервала решениями систем (I), (II) и (III) указаны в табл. 11 – 13. Как видно из таблиц, учет изменения индекса  позволил повысить точность аппроксимации более чем в три раза.

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

измерений и параметров сглаживающих решений для интервала 6.

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

сглаживающих решений для интервала 6.

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

решений для интервала 6.

На рис. 7, 8 изображены диаграммы рядов остатков  для решений систем (I) и (II). Сравнение рисунков показывает значительное уменьшение систематической составляющей в остатках  и  в случае аппроксимации данных измерений решением системы (II).

4. Заключение. Проведенное исследование показало, что учет изменения параметра  позволяет в несколько раз увеличить точность аппроксимации данных измерений на интервалах времени более 5 сут. При аппроксимации данных измерений на интервалах времени менее 3 сут выигрыш в точности практически отсутствует.

Литература

1.      Ветлов В.И., Новичкова С.М., Сазонов В.В., Матвеев Н.В., Бабкин Е.В. Режим гравитационной ориентации Международной космической станции. Космические исследования, 2001, т. 39, № 4, с. 436-448.

2.      Ланцош К. Практические методы прикладного анализа. М., Физматгиз, 1961.

Приложение