Разработка методов выявления опасных сближений космического аппарата с наблюдаемыми объектами и способов уклонения от столкновений на фазирующей орбите тема диссертации и автореферата по ВАК РФ 05.07.09, кандидат наук Каратунов Максим Олегович

• Каратунов Максим Олегович
• кандидат наук
• 2020

Оглавление диссертации кандидат наук Каратунов Максим Олегович

ВВЕДЕНИЕ

ГЛАВА 1. ВЫЯВЛЕНИЕ СБЛИЖЕНИЙ КОСМИЧЕСКИХ ОБЪЕКТОВ В ОКОЛОЗЕМНОМ КОСМИЧЕСКОМ ПРОСТРАНСТВЕ

1.1 Постановка задачи и общий подход к решению

1.2 Анализ существующих подходов к предварительной фильтрации

1.3 Предварительная фильтрация объектов, заведомо не представляющих угрозу для защищаемого космического аппарата

1.3.1 Отсев по высоте полёта

1.3.2 Отсев по компланарным элементам орбит

1.3.3 Отсев по некомпланарным элементам орбит

1.3.4 Отсев по фазам движения

1.4 Положение точки максимального сближения

1.5 Анализ работы алгоритма фильтрации

1.6 Выбор размера защищаемой области

1.7 Определение параметров сближения

1.8 Численное моделирование орбитального движения

1.8.1 Гравитационное влияние внешних небесных тел

1.8.2 Атмосферное торможение

1.8.3 Давление солнечного радиационного излучения на поверхность КА

1.8.4 Гравитационное влияние внешних небесных тел

ГЛАВА 2. ОЦЕНКА СТЕПЕНИ ОПАСНОСТИ СБЛИЖЕНИЯ КОСМИЧЕСКОГО ОБЪЕКТА С ЗАЩИЩАЕМЫМ КОСМИЧЕСКИМ АППАРАТОМ

2.1.Анализ существующих методов оценки степени опасности сближения космического объекта с защищаемым космических аппаратом

2.2.Постановка задачи и общий подход к решению

2.3.Генерация массива псевдослучайных векторов

2.4.Выявление факта пересечения трехмерных моделей КО

2.5.Оценка точности, условие автоматической остановки расчета

2.6.Примеры расчета вероятности столкновения

ГЛАВА 3. УЧЕТ МАНЕВРОВ АКТИВНЫХ КА

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

3.2. Оценка параметров импульсных маневров

3.2.1. Оценка параметров одноимпульсного маневра

3.2.2. Оценка одноимпульсного маневра с учетом ошибок определения параметров орбит

3.2.3. Оценка двухимпульсного маневра

3.3.Примеры оценки импульсных маневров

3.3.1. Численное моделирование

3.3.2. Оценка одноимпульсных маневров КА Электро-Л2 и КА Луч-5В

3.3.3. Оценка маневров космического объекта БСК

3.4. Оценка параметров продолжительных маневров

3.4.1. Оценка одиночного продолжительного маневра

3.4.2. Оценка продолжительного двухимпульсного маневра

3.4.3. Оценка продолжительного маневра с учетом ошибок определения параметров орбит

3.5. Примеры оценки продолжительных маневров

ГЛАВА 4. МАНЕВРЫ УКЛОНЕНИЯ

4.1.Расчет одноимпульсных маневров уклонения

4.2. Метод обеспечения безопасности четырехимпульсной динамической операции встречи на околокруговых орбитах

ЗАКЛЮЧЕНИЕ

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

ПРИЛОЖЕНИЕ 1. Акт о внедрении результатов диссертационной работы

152

Введение диссертации (часть автореферата) на тему «Разработка методов выявления опасных сближений космического аппарата с наблюдаемыми объектами и способов уклонения от столкновений на фазирующей орбите»

ВВЕДЕНИЕ

Актуальность темы исследования. Все космические аппараты (КА) находятся под угрозой разрушения из-за высокой вероятности столкновения с элементами космического мусора (КМ), который заполняет околоземное космическое пространство (ОКП). По последним данным, количество наблюдаемых космических объектов (КО) в ОКП составляет от 21 тысячи [1] до 25 тысяч [2], при этом доля действующих КА составляет 4.8-6.2%. Ежегодно центры управления космическими полётами регистрируют до 70 потенциально опасных сближений с отдельными активными КА [3]. При этом относительная скорость объектов достигает значений, при которых любое столкновение гарантированно приводит к разрушению. В то же время согласно отчёту европейского космического агентства [2] ежегодно регистрируется около 8 разрушений КО по различным причинам.

На основании прогнозов численности популяции КО [4, 5, 6] и роста числа запусков ракетоносителей можно сделать однозначный вывод: число объектов техногенного происхождения в ближайшие годы будет стремительно расти, а вместе с ним будет расти и вероятность взаимного столкновения. Это может привести к лавинообразному росту числа КО (так называемому синдрому Кесслера [7]) и полной непригодности ОКП к полетам. Существуют несколько методов, позволяющих снизить темпы роста количества КО:

- уклонение от столкновений;

- увод с орбиты КА по истечении срока активного существования;

- предотвращение разрушений в результате взрывов топливных баков и аккумуляторных батарей;

- очищение ОКП специализированными КА.

Регулярное уклонение от сближений и увод с орбиты хотя бы четырёх крупногабаритных фрагментов в год позволят остановить рост популяции КО [8].

Степень разработанности темы исследования. Наиболее реализуемым методом обеспечения безопасности космических полетов в ОКП на данный момент является оперативное выявление опасных сближений защищаемого КА с каталогизированными КО и последующее уклонение от столкновения. Данному вопросу посвящено более 1 50 научных работ. Исследования в данной области затрагивают различные аспекты проблемы, в числе основных можно выделить следующие группы и подгруппы:

1. Предварительная фильтрация.

1.1. Исключение из рассмотрения объектов, гарантировано не представляющих угрозу защищаемым КО (ЗКО) на всем интервале анализа.

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

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

1.4. Сравнение и тестирование методов фильтрации.

2. Поиск точки максимального сближения.

2.1.Вычисление минимального расстояния между двумя эллипсами с общим фокусом.

2.2.Нахождение локального минимума функции взаимного расстояния между КО.

3. Определение границ временного интервала сближения.

4. Учёт влияния возмущений и ошибок исходных данных.

5. Расчёт вероятности столкновения.

5.1. Оценка риска столкновения КО, сближающихся на короткий промежуток времени.

5.2. Вероятность столкновения КО при продолжительном сближении.

5.3. Применение методов Монте-Карло для расчёта вероятности столкновения.

5.4. Учёт формы и ориентации сближающихся объектов.

6. Применение алгоритмов параллельных вычислений.

В рамках обозначенной проблемной области существует ряд обзорных работ, среди которых стоит отметить [9], [10], [11]. Схожими задачами, возникшими задолго до проблемы космического мусора, являются задачи анализа столкновений естественных космических объектов в Солнечной системе [12], [13], [14].

На текущий момент разработаны, по меньшей мере, 9 систем и программных комплексов: COMBO (CspOC, USA), CRASS (ESA), STK Advanced CAT (AGI Inc. USA), ShadowCAT (AGI Inc. USA), SOCRATES (CSSI, USA), CSive (Aerospace Corp., USA), CAESAR (CNES, French Republic), CAOS-D (AIAA, USA), АСПОС ОКП (Роскосмос, РФ).

В России исследования в области мониторинга ОКП и выявления опасных сближений КО проводятся организациями: АО «ЦНИИМаш», ОАО МАК «ВЫМПЕЛ», ИПМ им. М.В.Келдыша РАН, ИКИ РАН, АО «НПК «СПП», а также АО «АНЦ».

Подход к решению задачи прогнозирования опасных сближений КО с защищаемым КА можно обобщить в виде следующих последовательных этапов: предварительная фильтрация объектов, заведомо не представляющих угрозу для защищаемого КА, вычисление координатно-временных характеристик сближения, расчёт вероятности столкновения. В тоже время, сравнивая источники можно заметить, что постановки задачи поиска сближений каталогизированных КО различаются в двух аспектах. Первое отличие заключается в составе исходных данных. В одном случае исходными данными являются векторы состояния, заданных на начальный момент времени, в другом - массив изохронных эфемерид (координат и скоростей) всех КО, рассчитанных с некоторым шагом на всём интервале анализа. Второе отличие заключается в количестве защищаемых объектов, если множество ЗКО равно множеству каталогизированных объектов, то говорят, что решается задача «все на все», в другом случае, когда указываются конкретные ЗКО, говорят о задаче «список на каталог». Далее представлен

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

В рамках постановки, при которой исходными данными являются векторы состояния КО, заданные на начальный момент интервала анализа, классической считается работа [15], где предлагается следующий набор фильтров. Апогейно-перигейный, учитывает минимальную и максимальную высоту полёта. Геометрический фильтр, учитывает минимальное расстояния между траекториями, представленными в виде эллиптических кривых. Временной фильтр, учитывает моменты времени прохождения опасных участков траектории. Идеи, предложенные в исследовании [15] получили развитие во многих последующих работах [16], [17], [18], [19].

На базе ФГУП «ЦНИИМаш» был разработан метод «обеспечения безопасности полетов околоземных космических аппаратов при условии наличия космического мусора» [20], основанный на подходе [15] на пошаговом отсеве объектов каталога, заведомо не представляющих угрозу для защищаемого КА. Данный алгоритм обладает высокими показателями быстродействия и используется при регулярном прогнозе опасных сближений элементов космического мусора с МКС [21]. Вместе с тем, в методике присутствует существенная неточность - за точку минимального расстояния принимается одна из двух точек траектории защищаемого КА, лежащих на линии пересечения плоскостей орбит КА и КО [22]. Нетрудно доказать, что точка сближения может не принадлежать линии пересечения плоскостей орбит. Данный эффект усиливается с ростом эксцентриситета рассматриваемых орбит и с уменьшением угла между плоскостями. Кроме того, третий этап отсева [23] опирается на отмеченную выше ошибочную гипотезу. В качестве критерия отсева на данном этапе принимается разность фаз в момент прохождения КА линии пересечения плоскостей орбит. Данный момент, как было отмечено, может не являться критичным с точки зрения сближения.

К методам предварительной фильтрации также можно отнести работы Лабуткиной Т.В. и др. [24], [25], [26], [27], в которых предлагается концепция

быстрого прогноза конфликтности системы орбитальных тел, основанная на поэтапном выявлении опасных участков траектории (узлов конфликтов). Скорость вычисления на порядки превосходит прямой метод за счёт отсутствия необходимости численного прогноза движения, что является безусловным преимуществом. В то же время описанные алгоритмы обладают рядом недостатков. В работе [26] предложен дополнительный этап анализа степени опасности сближения, основанный на следующем допущении: опасные участки траекторий сближающихся объектов представляют собой дуги окружности одного радиуса.

При этом, аппроксимация траекторий дугами окружности может привести к существенной погрешности в случае, когда точка максимального сближения не лежит на линии пересечения орбит. Данная ситуация характерна для объектов на геостационарных орбитах (ГСО) и траекторий с большим эксцентриситетом. Также в работе [24] присутствует математическая неточность, допущенная при рассмотрении случая совпадающих плоскостей орбит: расстояние между двумя эллипсами, ограничивающими область возможного местоположения КО, является постоянным только в случае совпадения центров эллипсов, в то время, как, общим является фокус эллипсов. В результате, в окрестности перицентра эллипсы ограничивают зону, размеры которой меньше критической сферы, что в свою очередь, может привести к пропуску потенциально опасных ситуаций. Помимо упомянутых статей Лабуткиной Т.В., различные аспекты проблемы загрязненности ОКП автором были затронуты в работах [28], [29], [30].

Этап расчёта координатно-временных характеристик сближения широко освещён в литературе. В работе Хутаровского З.Н. [31] приведён алгоритм, основанный на трёхэтапной процедуре прогноза движения КО. В статье Холшевникова К.В. и Балуева Р.В. [32] получено аналитическое решение задачи поиска экстремума функции взаимного расстояния между двумя кеплеровыми орбитами в строгой постановке. В работе [33] предложен альтернативный подход к поиску точки максимально сближения, основанный на разложении функции относительного положения в ряд Тейлора высоких порядков, с последующим

использованием стандартизированного программно-аналитического модуля поиска глобального экстремума COSY-GO [34].

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

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

Помимо расчёта координатно-временных характеристик сближения в работах [31] [35] [36] предлагается метод расчёта вероятности столкновения, который получил широкое распространение. Авторам удалось получить аналитическую формулу при введении следующих допущений:

1) КО имеют сферическую форму;

2) в пределах интервала сближения движение прямолинейно;

3) составляющие скоростей КО определены без ошибок;

4) размер эллипсоида рассеивания сближающегося КО значительно превосходит габариты защищаемого КО;

5) ошибки определения параметров движения распределены по нормальному закону.

Данный подход к расчёту вероятности хорошо работает для объектов, находящихся на низких орбитах, так как измерение траектории происходит при помощи радиолокационных средств, и ошибка определения скорости обычно не превосходит 1 м/с. Параметры движения объектов на более высоких орбитах определяются на основе оптических измерений, что обуславливает более высокую погрешность. В этом случае допущение 3 может привести к серьёзным ошибкам в расчёте вероятности столкновения. Допущение 4 не позволяет применять данный метод к крупногабаритным объектам, таким как, МКС.

В работе [37] авторы сводят задачу расчёта вероятности столкновения к нахождению вероятности попадания вектора относительного положения в критической точке опасного сближения в сферу заданного радиуса. В данной постановке вероятность определяется объёмным интегралом при помощи аппроксимации кубами малого размера, либо с использованием асимптотических рядов. К недостаткам данного подхода можно отнести следующие:

- расчёт вероятности происходит только в точке опасного сближения, окрестность данной точки не рассматривается;

- в отличие от метода Хуторовского З.Н. [31] применяются численные методы вычисления интеграла;

- форма и ориентация сближающихся объектов не учитывается;

Описание методов определения вероятности столкновения КО также встречается в работах [33], [38], [39], [40].

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

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

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

- необходимость проведения регулярных расчетов для большого числа КО, что приводит к проблеме быстродействия алгоритмов;

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

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

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

всегда имеется возможность провести дополнительную коррекцию с целью уклонения от столкновения.

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

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

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

Достижение поставленной цели потребовало решения следующих научно-технических задач:

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

2) Разработка метода расчёта вероятности столкновения с учётом формы и ориентации сближающихся объектов.

3) Разработка метода учета маневров активных КА при поиске опасных сближений.

4) Разработка метода расчета параметров безопасной четырёхимпульсной динамической операции встречи на околокруговых орбитах.

Объектом исследования является совокупность каталогизированных КО в

ОКП.

Предметом исследования являются методы выявления и оценки степени опасности сближений космического аппарата с пассивными и маневрирующими космическими объектами, а также методы расчета маневров уклонения.

Методология и методы исследования. Основным методом диссертационного исследования является математическое моделирование, основанное на теории космического полёта. При разработке аспектов маневрирования активных КО были использованы методы теории маневрирования космических аппаратов в окрестности круговой орбиты [44]. Разработка критерия оценки степени опасности сближения производилась в рамках теории вероятности с применением статистического моделирования. Программная реализация разработанных методов производилась в рамках методологии объектно-ориентированного программирования.

Положения, выносимые на защиту:

1) Методы выявления опасных сближений и расчета вероятности столкновения защищаемого КА с наблюдаемыми космическими объектами.

2) Метод учета маневров КО при выявлении факта сближения с защищаемым КА.

3) Метод расчета параметров маневра встречи, обеспечивающих уклонение от столкновений на фазирующей орбите.

Научная новизна диссертационной работы состоит в следующем:

- в рамках метода предварительной фильтрации предложен новый критерий отсева для КО на околокруговых орбитах, основанный на применении линеаризованных уравнений движения;

- разработан новый метод расчёта вероятности столкновения КО, учитывающий форму и ориентацию объектов на всем интервале сближения;

- предложен новый способ учёта маневров активных КА при выявлении опасных сближений КО, основанный на апостериорной оценке параметров динамических операций;

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

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

- более чем на порядок сократить время выявления фактов возможного сближения КО по сравнению с методом прямого расчета, не прибегая к допущению о положении точки максимального сближений на линии пересечения плоскостей орбит;

- избежать допущений классических методов расчета вероятности столкновения об отсутствии ошибок определения составляющих скоростей и о прямолинейности движения;

- определить вероятность столкновения с учетом формы и ориентации КО на всем интервале сближения;

- задать требуемую точность и достоверность оценки вероятности столкновения двух КО;

- существенно повысить точность поиска опасных сближений с маневрирующими КО в условиях отсутствия априорной информации о параметрах маневра;

- без дополнительных специальных маневров избежать столкновения на орбите фазирования при проведении четырехимпульсной динамической операции встречи.

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

Внедрение результатов работы. Алгоритмы, разработанные на основе предложенных методов выявления и оценки степени опасности сближений космического аппарата с пассивными и маневрирующими космическими объектами, используются в программно-аналитическом комплексе обработки траекторной информации в АО «Астрономический научный центр», что подтверждается соответствующими актами о внедрении (Приложение 1).

Работа выполнена при финансовой поддержке Министерства образования и науки Российской Федерации в рамках мероприятия 1.2 федеральной целевой программы «Исследования и разработки по приоритетным направлениям развития научно-технологического комплекса России на 2014-2020 годы» (Соглашение от 26 сентября 2017 года № 14.574.21.0146, уникальный идентификатор работ RFMEFI57417X0146).

Апробация работы. Основные положения и результаты диссертации докладывались и обсуждались на российских и международных конференциях:

- Международная молодёжная научная конференция «Гагаринские чтения», Москва, апрель 2012;

- XXXVIII академические чтения по космонавтике, Москва, февраль 2014;

- International workshop on "Key Topics in Orbit Propagation Applied to Space Situational Awareness" (KEPASSA). 18-30 Oct. 2015. Toulouse (France);

- XL академические чтения по космонавтике, Москва, январь 2016;

- XLI академические чтения по космонавтике, Москва, январь 2017;

- 3rd IAA Conference on Dynamics and Control of Space Systems (DYCOSS). 30 May - 1 June. 2017. Moscow (Russia);

- Шестые Репинские научные чтения, Москва, ноябрь 2018;

- Всероссийская конференция с международным участием «Космический мусор: фундаментальные и практические аспекты угрозы», Москва, 2019.

Публикации. Результаты работы изложены в 6 статьях [45], [46], [47], [48], [49], [50], [51], [52], опубликованных в рецензируемых научных изданиях, рекомендованных ВАК РФ, в том числе 3 статьи [50], [51], [52] опубликованы в журналах, входящих в базы данных Web of Science и/или Scopus.

Структура и объём диссертации. Диссертация состоит из введения, четырех глав, заключения и выводов, списка литературы и приложения. Объём диссертации - 152 страницы. Работа включает в себя 43 иллюстраций и 14 таблиц. Список литературы содержит 102 наименования.

ГЛАВА 1. ВЫЯВЛЕНИЕ СБЛИЖЕНИЙ КОСМИЧЕСКИХ ОБЪЕКТОВ В ОКОЛОЗЕМНОМ КОСМИЧЕСКОМ ПРОСТРАНСТВЕ

1.1 Постановка задачи и общий подход к решению

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

Исходными данными являются:

- границы интервала поиска сближений [¿нач, ¿кон];

- вектор состояния защищаемого КА на момент времени ¿нач;

- векторы состояния наблюдаемых КО на момент времени ¿нач;

- ковариационные матрицы векторов состояния КА и наблюдаемых КО;

- параметры математической модели, использованные при решении краевой задачи для определения каждого вектора состояния.

Выходными данными является массив, элементы которого содержат информацию по каждому отдельному сближению, а именно:

- минимальное расстояние между КА и КО;

- момент времени максимального сближения ¿сб;

- векторы состояния КА и КО на момент времени ¿сб;

- временной интервал сближения [¿1,?2], на котором взаимное расстояние между КА и КО не превышает критического значения.

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

каталоге может быть отсеяно при помощи упрощенных зависимостей на предварительном этапе. После этапа предварительного отсева производился поиск минимального расстояния между КА и КО при помощи численного интегрирования уравнений движения с учётом возмущений. На этом этапе также находится время максимального сближения tсб. Далее с учетом возможных отклонений от номинальной траектории определяется интервал сближения [^2].

1.2 Анализ существующих подходов к предварительной фильтрации

Рассмотрим особенности существующих подходов к предварительной фильтрации. В рамках постановки, при которой исходными данными являются векторы состояния КО, заданные на начальный момент интервала анализа, классической считается работа [15], где предлагается следующий набор фильтров.

Апогейно-перигейный фильтр, учитывающий требование пересечения интервалов высот ЗКО и сближающегося КО (СКО). Условие фильтрации выражается следующим отношением:

, [15] (1.1)

где q - больший из двух радиусов перигеев, Q - меньший из двух апогеев, В -размер защищаемой области.

Геометрический фильтр основан на поиске минимального расстояния между эллиптическими траекториями без учёта положения КО на данных траекториях. В рамках данного фильтра автор рассматривает только случай некомпланарных орбит. Определяется положение линии пересечения плоскостей с последующим поиском локального минимума в данной окрестности, представленной на рисунке 1.1. Взаимное расстояние между объектами ге/ выражается в виде функции углов и^, , отсчитываемых от линии пересечения

плоскостей до текущего положения ЗКО и СКО соответственно. Функция имеет вид:

ГЫ = Г1 + Г - 2грГ с°8 У, [15] (1.2)

где г, Г - модули радиус-векторов ЗКО и СКО соответственно, у - угол между радиус-векторами ЗКО и СКО.

Список литературы диссертационного исследования кандидат наук Каратунов Максим Олегович, 2020 год

/ / / /

г/

0 100

20

15 &

3

Ю °

200 300 400 500 600 700 Защищаемая область [км]

800

900

0

1000

Рисунок 1.10 - График зависимости количества объектов, представляющих угрозу от величины

защищаемой области для объекта №2

250

200

са

О &

Щ ев

ю о о и

п о

150

100

50

1 1 у —'

— -О -о .............В зъекты по □асные об эемя рабоп слс фильт ьекты гы 1ЯЦИИ

у у у * -

У ✓ ✓ ** 1* -

У * у -

*

50

100 150 200 250 300 350 Защищаемая область [км]

400

450

3 -

о ю га

3 а

1

500

Рисунок 1.11 - График зависимости количества объектов, представляющих угрозу от величины

защищаемой области для объекта №3

1.6 Выбор размера защищаемой области

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

Произведём оценку радиальной (Лг), трансверсальной (А/) и бинормальной (Лг) составляющей отклонения, вызванного ошибками определения орбиты. Прибегнув к линеаризации уравнений движения в цилиндрической системе координат (СК) [64], получим следующие выражения для отклонений:

Лг ЪУ ЪУ

— = Ъг0 +—- бШ ф + 2—L (1 - СОБ ф), (1.45)

г0 У0 У0

— = Ъ/0 - 2 (1 - соб ф) - ЪУ (3ф - 4бш ф), (1.46)

г0 У0 У0

Лг Ъ У

— = Ъг0 + —г бШ ф. (147)

г0 У0

Здесь г0 - радиус опорной круговой орбиты; У0 - модуль вектора скорости движения по опорной орбите; 8г0, 8/0, 8г0 - ошибки определения положения КО

в начальный момент времени; 8УГ, 8У(, дУ2 - совокупное отклонение составляющих вектора скорости КО от У0 , обусловленное наличием ненулевого

эксцентриситета и ошибками определения параметров орбиты в начальный момент времени; ф - угол, отсчитываемый в плоскости орбиты от радиус вектора КО в начальный момент времени в направлении его движения.

Исследуя функции (1.45), (1.46), (1.47) на экстремумы, нетрудно найти максимальное значение отклонений |Дг| , |Д/| , |Д^| на интервале

I 1тах ' I 1тах' I 1тах г

Ф <= [0, ф^ ]. Угол фу соответствует положению КО в конечный момент времени. Стоит отметить, что отсчёт угла ф у производится с учётом всех витков, пройденных КО на интервале анализа (в общем случае фу > 2 л).

Оценку составляющих Дгр, Д/р, Лгр отклонения, вызванного

возмущающими факторами, можно осуществить на основе упрощённых зависимостей [63].

Отклонения, вызванные сжатием Земли:

|Дв/\0 * — ла 2 - Збш2 I + 3вт21 соб2фу

г0 V

(1.48)

|Д

а„

в" \о

2 г

ла

|б1П2/| .

(149)

1ДпЪ = КЪ * - а

Г0

(1.50)

|Д

а„

п 2\о

2 г

а,

(1.51)

где Д 1\ , |Лв- вековые отклонения в трансверсальном и боковом

направлениях; |ДИ 1\ , |Д- периодические отклонения в трансверсальном и

боковом направлениях; ае - большая полуось Земного эллипсоида; а -коэффициент сжатия Земли.

Отклонения, вызванные аэродинамическим торможением:

К lLo * 6*2 W

|А

r *—4ncprr

в \aero г 0

М

aero

Ia

п \aero

12тс2срг02 14.8

12тс2срг02 59

(1.52)

(1.53)

(1.54)

(1.55)

где |Д l\ , Ia r\ - вековые отклонения в трансверсальном и радиальном

I в ¡aero I в ¡aero

направлениях; |ДИl| , \\Ааего - периодические отклонения в трансверсальном и

радиальном направлениях; с - баллистический коэффициент; р - плотность воздуха.

Для определения р в контексте данной задачи целесообразно использовать простейшую (изотермическую) модель атмосферы. Расчёт влияния аэродинамического сопротивления проводится для орбит r0 < 800 км. Отклонения, вызванные притяжением Солнца и Луны:

Mms * 3-25™0

/11

max

|a

вZMS * °-75™0

/11

Яо

(1.56)

(1.57)

где |Дв/| , |Дв- вековые отклонения в трансверсальном и боковом

/1

направлениях;

- отношение максимального возмущающего ускорения к

тах

ускорению земного притяжения на экваторе (определяется по известным таблицам).

Найдём суммарные вековые отклонения:

V = |Д4 + |Д„/|_ + , (1.58)

НМ^ (1.59)

Д ^=К 4 +Д Аш • (160)

Отклонения, накопленные на конечный момент времени с учётом максимальных периодических отклонений:

Д'р=Д-' Ъ+|Д"/|°+|Д-/|-, (1.61)

ДГр = ДВ в + 1-ДП, (1.62)

Д-р = Д,В + |Д„.'|е О.Ю)

2л

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

Существует несколько работ [19], [65], [66] посвящённых решению данной проблемы выбора размера защищаемой области. В исследовании [66] предлагается выбирать пороговое значение расстояния основываясь на предельно допустимой вероятности столкновения. Авторам удалось вывести аналитическую зависимость связывающую вероятность столкновения с размером сближающихся

КО, отношением сторон обобщённого эллипса рассеивания и взаимным расстоянием в точке максимального сближения.

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

Р =

1

2-я-ахау }г^

р г г\1г2 —.

х — х„

\2 г

К а х У

Л2

У — Ут

к ау У

с!ус!х,

[67] (1.64)

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

характеризующие оси суммарного эллипса рассеивания.

Стоит упомянуть, что уравнение (1.64) даёт истинное значение вероятности столкновения в рамках ряда допущений, подробное рассмотрение которых представлено в соответствующем разделе 2.1.

Выражение (1.64) не имеет аналитического решения и крайне неудобно для расчета порогового значения допустимого расстояния. В этой связи автор [66] использует подход предложенный [68], при котором двумерная задача сводится к одномерной путём приведения эллипса к окружности (сжатие по оси х). В данном случае выражение для определения вероятности имеет следующий вид:

Р = - (1 — в~и!2) + Т, [66] (1.65)

где и = г У (а х а у), [66] (1.66)

* = ( х2т/а2х ) + (у2т/а2у ), [66] (1.67)

Т < 1 - и2 - 3 - г*!2 - 4, [66] (1.68)

здесь Т - остаточный член разложения.

Далее используется необходимое условие максимума вероятности [69], согласно которому ут = 0, и вводится параметр соотношения сторон эллипса рассеивания:

ЛЯ = а ^ а, [66] (1.69)

При этом уравнения (1.66), (1.67) принимают вид:

и = г у (ЛЯ • а) [66] (1.70)

3 = й2/(ЛЯ2 •а2). [66] (1.71)

Затем, фиксируя значения ЛЯ, г и d, определяется частная производная вероятности по а . Приравнивая нулю производную, находится выражение для

максимума вероятности:

=

С \ а ( 1 ]

1(! + а)) 1(! + а))

а = (ЛЯ • гVй2).

[66] (1.72) [66] (1.73)

Для получения зависимости расстояния й от максимума вероятности Рах автор использует разложение в степенной ряд. Заключительное аналитическое соотношение имеет следующий вид:

а = р-

40 4

--ъ —,

49 •р 7

[66] (1.74)

Р =

12•е• Рт 7

144 • (е • Ртах)2 6528 • е • р^ , 2816 272

49

2401

2401 343

[66] (1.75)

В качестве примера применения метода приводится рисунок 1.12, на котором изображены линии соответствующие различным пороговым значениям

максимальной вероятности, цифры на линиях обозначают порядок величины (10й ). По оси абсцисс отложен суммарный размер объектов г, по оси ординат расстояние между объектами d. Соотношение осей ЛЯ совмещённого эллипса рассеивания для данного примера равен 3.

10° 101

Combined Object Radius (meters)

Рисунок 1.12 - Линии равной максимальной вероятности столкновения КО

ЬВюС^ах), АЯ = 3 [66]

Таким образом, задавая допустимый уровень вероятности и располагая информацией о размерах КО, можно определить минимально допустимое

пороговое значение взаимного расстояния. Например, для Рах = 10~5 и г = 30 м, на рисунке 1.12 область обозначенная окружностью, пороговое значение расстояния d = 10 км. Сближения с минимальным расстоянием, превышающим найденное значение, будет безопасным в рамках заданного уровня вероятности столкновения.

Методическая погрешность предложенного подхода обусловлена остаточными членами разложений (1.65) и (1.74). В рамках исследования [66] был проведен анализ величины данной погрешности на основе 4.24 млрд. симуляций.

Исходные данные для каждой симуляции выбирались случайным образом в диапазонах: г. е(1...100) м, е(1...100)км , АЯ е(1...100). За истинное значение вероятности принимался (1.64), затем по формулам (1.74) и (1.73) вычислялась оценка взаимного расстояния . Расхождения ¿/г - в зависимости от коэффициента аг изображены рисунке 1.13.

Рисунок 1.13 - Анализ методической погрешности [66]

Точки, формирующие огибающую линию, соответствуют максимальным отклонениям, полученным в результате симуляций. Абсолютный максимум погрешности равен 1880 м, величина минимума составила всего -2.5 м. В данном

анализе учитывался знак отклонения - ^, положительная величина означает, что пороговое значение взаимного расстояния, полученное на основе предложенного метода, будет иметь некоторый запас. Другими словами, положительная погрешность не приведёт к пропуску опасных ситуаций. Чтобы избежать пропусков из-за отрицательной погрешности достаточно ввести небольшой запас 2.5 м.

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

В рамках апогейно-перигейного фильтра авторами работы [19] предлагается анализировать значения радиусов апоцентра и перицентра на концах интервала поиска сближений с целью учёта вековых смещений. Суммарная прибавка к защищаемой области складывается из вековых уходов и максимальных короткопериодических флуктуаций. Оценка величины данной прибавки, полученной на основе анализа движения более 11 тысяч объектов каталога, составляет 30 км. При этом авторы отмечают, что оценка является приблизительной и не гарантирует отсутствие пропусков.

В результате тестирования тороидального геометрического фильтра [53] автор приходит к выводу о необходимости дополнительной оценки величин, обеспечивающих отсутствие ошибок второго рода. В рамках рассмотрения случая пропуска опасного сближения КО «SL-3 Rocket Body (NORAD ID 09904)» и КО «Fengyun 1С debris (NORADID 31921)» была выявлена существенная неустойчивость величины минимального расстояния между орбитами при изменении исходных параметров орбит. При этом различие в элементах орбит, в силу возмущающих факторов, может быть обусловлено моментом времени, на который данные элементы получены. График данной зависимости изображён на рисунке 1.14, где по оси абсцисс отложены моменты времени получения исходных данных, по оси ординат минимальное расстояния между орбитами (не объектами). Можно заметить, что минимальное расстояния варьируется от 0 до 155 км, таким образом, чтобы избежать пропуска сближения размер защищаемой области должен быть больше 155 км. Такой размер области приведёт к увеличению количества ошибок первого рода и снижению эффективности алгоритма фильтрации.

Рисунок 1.14 - Зависимость минимального расстояния между орбитами от момента времени

получения исходных данных [53]

Уменьшить необходимый размер защищаемой области можно прибегнув к идее [66] разбиения интервала анализа на более короткие временные периоды и многократного применения фильтра для каждого периода в отдельности. На рисунке 1.15 представлен график зависимости диапазона изменчивости минимально расстояния между траекториями от продолжительности периода анализа. Так, если изначальный интервал анализа разбить на полусуточные периоды, то для исключения ошибок второго рода достаточно V = 89 км, Оь = 70 км.

<D

<3

160 140 120 100 80 60 40 20

+ + + * + + + Ф + + + > О О О О О - + + + + + > О О О о о - + + +■ + + >00000

+ + + . + + + * о о

+ о* * О

+ о о

+ о D D о о 1 □

Ф о о 1 □ □ □

□ □ +++ Total 000 In Plane 1 Out of Plane

0.5

1.5

2.5

3.5

Time Sampling Interval (days)

Рисунок 1.15 - Зависимость величины диапазона изменения минимального расстояния между орбитами от продолжительности периода разбиения. [53]

В работе [70] предложен обобщённый подход, позволяющий оценить взаимосвязь между конкурирующими характеристиками, а также определить минимальную величину запаса при заданном уровне точности. Авторы применяют методы теории обнаружения сигналов и концепцию рабочей характеристики приёмника ^ОС-кривой), которая заключается в построении графика зависимости между долей истинно-положительных сигналов и долей ложно-положительных сигналов. В терминах задачи опасных сближений КО истинно-положительным сигналом является пара объектов, прошедшая отсев и сближающаяся на опасное расстояние, ложно-положительным сигналом является пара объектов, также прошедшая отсев, но при этом не представляющая опасность.

Для построения ROC-кривой осуществляется поиск опасных сближений «все на все» на заданном срезе каталога с включенной и выключенной фильтрацией. Изменяя пороговое значение и величину запаса можно построить графики количества ошибок первого и второго рода, представленные на рисунках 1.16 и 1.17 [70] соответственно. Если задать допустимый уровень ошибок второго рода, то на основе графика 1.17, можно определить минимально допустимую величину запаса. Характеристическая кривая, изображённая на рисунке 1.18, позволяет наглядно представить взаимосвязь критериев быстродействия (ось абсцисс) и точности (ось ординат). Для оператора интерес может представлять точка баланса, где производная ROC-кривой равна единице, или окрестность стопроцентной горизонтальной линии отсутствия пропусков.

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

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

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

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

Рисунок 1.16 - Доля ошибок первого рода апогейно-перигейного фильтра [70]

Рисунок 1.17 - Доля ошибок второго рода апогейно-перигейного фильтра [70]

percent false positive

Рисунок 1.18 - ROC-кривая апогейно-перигейного фильтра [70]

1.7 Определение параметров сближения

После применения описанных выше алгоритмов фильтрации для оставшихся объектов производится поиск минимального расстояния до «защищаемого» КА с учётом возмущений. Для сокращения времени расчёта целесообразно использовать двухэтапный подход.

На перовом этапе применяется численно-аналитическая теория прогнозирования движения [71], основанная на применении специальных возмущающих функций А.Р. Голикова. Численно-аналитический прогноз позволяет быстро оценить величину минимального расстояния с известной точностью. В случаях, когда найденное расстояние оказывается меньше защищаемой области, умноженной на коэффициент запаса, запускается второй этап поиска. На втором этапе осуществляется высокоточное численное интегрирование уравнений движения с учётом возмущений. При этом с заданным шагом определяется знак проекции относительной скорости на ось, соединяющую КА и КО. Интервал, на котором проекция меняет знак с положительного на отрицательный, содержит локальный минимум расстояния. Далее, методом золотого сечения внутри данного интервала, определяется минимальное расстояние и время сближения. Если найденное расстояние меньше заданного, то производится расчёт всех остальных параметров, таких как интервал опасного сближения, вектор относительной скорости, вектор относительного положения на момент сближения и др. За интервал опасного сближения принимается интервал времени, на котором эллипсоиды рассеивания двух КО имеют общие точки. Задача поиска момента начала и конца опасного интервала требует проверки факта пересечения эллипсоидов с некоторым шагом. Проверка базируется на аналитическом алгоритме, предложенном в работе [72].

Как было отмечено, разработанные алгоритмы внедрены в технологический цикл поддержания каталога КО, наблюдаемых при помощи оптических средств, и осуществляют регулярный поиск опасных сближений для ряда геостационарных спутников. Список защищаемых КА представлен в таблице 1.2. В таблице 1.3

приведены наиболее опасные сближения защищаемых КА, обнаруженные в период с февраля по сентябрь 2016 года.

Таблица 1.2 - Список защищаемых КА

№ Наименование Международное Номер

обозначение NORAD

1 EXPRESS MD1 2009-007B 33596

2 EXPRESS AM-33 2008-003A 32478

3 RADUGA 1M-1 2007-058A 32373

4 EXPRESS AM-3 2005-023A 28707

5 EXPRESS AM-2 2005-010A 28629

6 EXPRESS AM-1 2004-043A 28463

7 ELEKTRO-L 2011-001A 37344

8 COSMOS 2473 2011-048A 37806

9 AMOS 5 2011-074A 37950

10 YAMAL 300K 2012-061B 38978

11 YAMAL 402 2012-070A 39022

12 RADUGA 1M-2 2010-002A 36358

Таблица 1.3 - Выявленные сближения защищаемых КА с наблюдаемыми КО

№ Минимальное расстояние [км] Дата и время сближения UTC Номер КА NORAD Номер Ш NORAD Доля отсеянных КО Время вычислений [сек]

1 0.719 02.02.16 22:52:02.100 36358 20473 96.9 % 77

2 0.804 05.02.16 21:49:17.251 37806 19344 95.4 % 93

3 1.964 08.02.16 09:19:43.511 28629 19344 96.7 % 93

4 2.288 06.02.16 09:47:25.594 37806 19344 95.4 % 93

5 2.354 07.02.16 08:27:50.294 32478 18575 96.1 % 93

6 2.466 09.02.16 21:13:45.210 28629 19344 96.7 % 93

7 3.687 03.04.16 23:34:11.739 38978 33749 96.0 % 64

1.8 Численное моделирование орбитального движения

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

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

- численные методы интегрирования дифференциальных уравнений, учитывающие специфику функций, описывающих орбитальное движение КА;

- программная реализация численных методов;

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

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

Ошибки численной модели, определяющие её точность, могут быть разделены на следующие классы:

- параметрические, связанные с неточностью физических констант;

- ошибки, вызванные неадекватностью математической модели физическому процессу;

- ошибки округления в компьютерной арифметике;

- ошибки аппроксимации в численном методе.

В рамках данной диссертационной работы в качестве метода численного интегрирования обыкновенных дифференциальных уравнений был выбран многошаговый метод Адамса-Мультона-Коуэлла с переменным шагом. Реализация метода написана на языке программирования высокого уровня C++. Математические модели возмущающих ускорений представлены в данном разделе. Уравнения орбитального движения в системе координат, связанной с неподвижным массивным центральным телом, представимы в виде [73]:

<2 х ц

2 3-х + Р(7,х,х) = Р + Р, (1.76)

¿А |х|

где х - вектор положения исследуемого тела, t - время, / - гравитационный параметр, Г - силы центрального гравитационного поля, а Р - возмущающие силы, причем |Р| |г|. В процессе прогнозирования движения КА учитывались

следующие возмущающие факторы:

- нецентральность поля тяготения Земли. Стандарты гравитационного геопотенциала: ПЗ-90.2, EGM96, JGM3, WGS72, WGS844;

- атмосферное торможение. Стандарты моделей плотности атмосферы Земли: ГОСТ Р 25645.166-2004, DMA2004, DMA2004m, MSIS90, аКЛ90;

гравитационное влияние Солнца и Луны; световое давление с учетом теневого эффекта.

1.8.1 Гравитационное влияние внешних небесных тел

Несферичность поля тяготения Земли описывается с помощью его силовой функции в гринвичской вращающейся системе координат, представляемой в виде разложения в ряд по сферическим функциям:

II f N f R X

KNa = Nl ^ p] Pn (sin ф) + r n=2 V r )

(1.77)

M n fn \n

+ \Cnm C0S mX + dnm Sin Pnm (Sin ф)

n=2 m=1 V Г )

где r - расстояние от центра планеты, ф - широта в гринвичской системе координат, X - долгота в гринвичской системе координат, | - гравитационная постоянная Земли, R - экваториальный радиус Земли, cnm, dnm - коэффициенты

гравитационного поля, Pnm (sinф) - присоединенные функции Лежандра,

Pn0 (sinф) = Pn (sinф) - многочлены Лежандра.

Модели гравитационного поля Земли (N*M) задаются набором коэффициентов, где N - максимальное значение индекса n (число зональных гармоник) и M - максимальное значение индекса m (для тессеральных и секторальных гармоник). Здесь стоит упомянуть результаты исследования А.И.Назаренко "Погрешности прогнозирования движения спутников" [74], где показано, что для современных моделей среднеквадратические отклонения (СКО) погрешностей коэффициентов разложения гравитационного потенциала Земли

имеют порядок 10 9. Поэтому нет необходимости учитывать столь громадные модели поля Земли.

Возмущающие ускорения rNC = {Зсдг; yNC; zNC} = grad У?дг, вызываемые

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

Г N J М п

NC 2 I / ' иО п+1,1 + 1 -SS[( CnmAn+1,m+1 + dnmBn+1,m+1 ) — ^ ^

Г 1и=2 ^ п=2 т=\ (1.78)

- (n - m + 2) • (n - m + 1) • (Cnm An+1,m—1 + dnm Bn+1,m—1 )]}

|И 1 п

Г 1и=2 ^ и=2 т=1 (1.79)

+ (п - т + 2) • (п - т +1) •(с В^ ,- < А ,)]},

V /V / ^ пт п+1,т-1 пт п+1,т-1 / I '

ц ( N М п Л

= - 4 ■■ +• с*0Л+1,о + -т +• + ) к (1 - 80)

п=2 т=1

где А =

^ п,т

^ т? V-1

р» (вт фу соб тк, вии1 =

пт \ т / ? п,т

V г У

^ Т? V-1

Л®

V г У

'Рпт (^ ф)* тк

При введенных обозначениях возмущающая силовая функция записывается как

Ц [ ^ М п I

^С =_ •] ^ сп 0 Ап,0 СптАп,т + <птВп,т ) к

(181)

п=2

п=2 т=1

Для вычисления значений Аит, Вит удобно использовать рекуррентные соотношения:

Ат+1,т+1 = (2т + 1) * {хАт,т - УВт,т }* Л /Г' Вт+1 ,т+1

= ( 2 т +

1 )•{ У Ат, т + хВт,т }•Л®/ г 2

для индекса т = п и

Ап,т = [(2п - 1) •• * Ап-1т - (п + т - 1) • Л® • Ап-2,т ] •

В =[(2п -1)* г В . -(п + т -1)* Л • В , ]

п,т / п-1,т \ / ш п-2,т ^

для индекса т < п .

п - т

п - т

(1.82)

(1.83)

(1.84)

(1.85)

Они легко выводятся из рекуррентных соотношений для присоединенных функций Лежандра:

ртт(х) = (2т -1 - х2 • Рт-1,т-1(х) (для индекса т = п), Роо(х) = 1,

рт+1т(х)= (2т +1) • хРтт (х) (для индекса т = п -1),

Рпт (х) = [(2п - 1) • хРп-т (х) - (п + т -1) • РРп-2,т(х)]/(п - т) (для индекса т < п -1).

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

М п

этом меняется порядок суммирования: вместо ^^ внешнее суммирование

п=2 т=1

М М

проводится по индексу т: ^^ .

т=1 п=1

1.8.2 Атмосферное торможение

Выражение для силы атмосферного торможения записывается в виде

Р^=~рСвА.Уа.Уа = -тр8Уа.Уа, (1.86)

где р - плотность воздуха, Уа - относительная скорость набегающего потока воздуха, А - площадь миделева сечения, Сд - безразмерный аэродинамический

С А

коэффициент лобового сопротивления, £ = —---баллистический коэффициент,

2т

т - масса космического аппарата.

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

основе траекторных измерений. Относительная скорость Уа спутника зависит от

сложной атмосферной динамики, но приемлемым приближением считается

допущение о том, что атмосфера вращается вместе с Землей почти с той же угловой скоростью юф. Следовательно, вектор относительной скорости

Га=Г-[щхг], (1.87)

где в инерциальной геоцентрической системе координат: V = {х;у;г} - скорость, г = [х;у;г} - радиус-вектор спутника, юф = {0;0;ооф } - вектор угловой скорости вращения Земли.

Максимально наблюдаемым отклонением в таком допущении - порядка 40% [75], что ведет к неопределенности в силе торможения не более, чем 5%. Компонентами от подъемной и боковой силы можно пренебречь.

Зависимость плотности верхних слоев атмосферы Земли задаётся формулой:

р = р0 к -к )/я°, (1.88) где р0, Щ - плотность и шкала высот на референтной высоте И0. Высота И над общеземным эллипсоидом зависит от радиус-вектора г и широты ф. С точностью О (а^), высота определяется по формуле:

к « г - - а® (г/г)2], (1.89)

где ^ - геометрическое сжатие Земли.

Величина Щ , в свою очередь, является функцией от высоты к . Поэтому в моделях для вычисления плотности ночной атмосферы ря применяется аппроксимация. В моделях ГОСТ 22721-77 и ГОСТ 25645.115-84 [76] приводится следующая зависимость:

(1.90)

где а0, а, а, а - коэффициенты модели.

В моделях ГОСТ Р 25645.000-2001 [77] и ГОСТ Р 25645.166-2004 [78]:

Современные российские модели (ГОСТ) плотности атмосферы Земли построены с использованием базы данных обработки измерений движения спутников за несколько десятков лет, что дает высокую точность и возможность использования в баллистических центрах.

Альтернативный подход применен в американских моделях серии (М818-77, М818-83, М818-86, М818Б-90) [79], а также во французских эмпирических термосферических моделях плотности ЭТМ [80], где расчет распределения плотности и температуры осуществляется по каждому газу в термосфере. Неудобство таких моделей при интегрировании заключается в сложности и громоздкости формул, что в свою очередь, ведет к увеличению времени вычислений. Помимо основной составляющей ря, определяющей характер распределения плотности, производится учет дополнительных физических факторов: суточный эффект (К1), полугодовой эффект (К2), долгопериодические изменения индекса солнечной активности ^ с циклом 81 суток (К0), краткосрочные изменения индекса солнечной активности ^0,7 внутри цикла 81 суток (К3), зависимость плотности атмосферы от геомагнитной возмущенности кр (К4). В скобках здесь указаны обозначения соответствующих коэффициентов в ГОСТ Р 25645.166-2004. С учетом вышеперечисленных факторов плотность атмосферы определяется следующей формулой:

р = рн • К0 •(! + К + К + Къ + КА),

(1.92)

Анализ современных динамических моделей плотности атмосферы Земли, базирующихся на эффективных таблицах и формулах (до высот 1500 км), показал, что отклонения от реальных значений плотности не превосходят 5^15%.

1.8.3 Давление солнечного радиационного излучения на поверхность КА

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

E = Es-^/Д)2, (1.93)

где £5 «1373 Вт/м2 - мощность солнечной радиации на поверхности Земли (солнечная постоянная), ^ - среднее расстояние Земли от Солнца.

Давление солнечной радиации прямо пропорционально мощности Е светового потока и обратно пропорционально скорости света и равно Г8 «4,56 -10-6 Н/м2 на расстоянии 1 А.Е. В зависимости от отражательных

свойств поверхности существуют полное поглощение, зеркальное и их комбинации. Как правило, чисто зеркальное отражение или чисто полное поглощение не существуют, на практике обычно применяют коэффициент е отражения, величина которого зависит от материала поверхности. Для типичных материалов, используемых при конструкции спутников, величина е находится между 0,2 и 0,9. В этом случае возмущающее ускорение от светового давления представимо в виде:

4р •Ч'-(Асо8е)'[(1_8)'% +2всо80-й]/т, (1.94)

где (A cos 9) - площадь поперечного (миделева) сечения светового пучка, освещающего рассматриваемый фрагмент поверхности, Э - угол между нормалью к поверхности Я и вектором направления на Солнце, cos 9 = (Я, es), т -масса КА, у - теневая функция.

В рамках данной диссертационной работы использовалась цилиндрическая модель тени Земли, при этом теневая функция рассчитывалась при помощи следующих формул:

0, если i!^l<-Jl-(Rjrf; г • г

v= ,__(1-95)

1,еспи^>->/1 -{Rjr)1-

r • rs

где Яе - радиус Земли, г - радиус-вектор спутника, вектор г5 - радиус-вектор Солнца в геоцентрической системе координат.

1.8.4 Гравитационное влияние внешних небесных тел

Согласно ньютоновскому закону притяжения, ускорение спутника в инерциальной системе координат описывается уравнением:

г — г г — г

InerUonal = GMa ' 77 Ze = ^ ' Г^ , (196)

\r _ r\ \r _ r\

G G

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

= (1-97)

г Л К

О

получаем возмущающее ускорение КА в геоцентрической системе координат:

Г = и

г - Г г