Астрономический вестник, 2022, T. 56, № 5, стр. 344-355

ВВЕДЕНИЕ

По данным службы контроля околоземного космического пространства (ОКП) на околоземных орбитах находится в настоящее время более 500 тыс. объектов искусственного происхождения, причем 20 тысяч из них имеют размеры, большие 10 см, и доступны для позиционных наблюдений (Оголев, Морозов, 2019). Вся эта масса объектов создает угрозу действующим аппаратам (Klinkrad, 2006), что приводит к необходимости постоянного отслеживания большого количества объектов. А это, в свою очередь, требует постоянного совершенствования программно-математического обеспечения, предназначенного для обработки результатов оптических наблюдений, определения и улучшения орбит, а также высокоточного численного моделирования движения (Андрианов и др., 2019).

Коллектив авторов настоящей работы состоит из двух групп, одна из которых имеет большой опыт в разработке методов высокоточного численного прогнозирования движения ИСЗ (Александрова и др., 2017; Авдюшев, 2020), а вторая группа (Левкина и др., 2013) ведет позиционные наблюдения на уникальной научной установке Цейсс–2000 в ЦКП “Терскольская обсерватория” Института астрономии РАН.

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

Разработанный ранее программный комплекс “Численная модель движения систем ИСЗ” (Александрова и др., 2017) был усовершенствован для определения параметров движения и сил из обработки результатов высокоточных позиционных наблюдений. Были уточнены модели сил в соответствии с (IERS Conventions 2010), использован более эффективный интегратор (Авдюшев, 2020), введено точное вычисление изохронных производных путем численного интегрирования соответствующих уравнений. В новом программном комплексе совместно с параметрами орбиты определяется коэффициент парусности, представляющий собой отношение площади миделевого сечения объекта к его массе. Этот параметр определяет величину влияния светового давления на динамику объекта.

Начальное определение параметров орбит в рамках данного программного комплекса находится нетрадиционным способом, который в применении к астероидно-кометной задаче был описан в (Батурин, Чувашов, 2006). Как отмечается в работе (Андрианов и др., 2019), в связи с тем, что классические способы первоначального определения орбит мало приемлемы при обработке наблюдений объектов космического мусора, интерес к этой задаче в последнее время возрос (2009; Milani и др., 2004; Tommei и др., 2007; Maruskin и др., Kolessa и др., 2014; 2019).

Для апробации были выбраны объекты с номерами 90 008, 90 031, 90 214 и 97 149. Номера объектов даны в соответствии с нумерацией в динамической базе данных космических объектов Института прикладной математики им. М.В. Келдыша РАН (Молотов и др., 2020). Полученные параметры использованы для исследования долговременной орбитальной эволюции объектов.

ОПИСАНИЕ ПРОГРАММНОГО КОМПЛЕКСА “ЧИСЛЕННАЯ МОДЕЛЬ ДВИЖЕНИЯ ИСЗ”

Численную модель орбитального движения ИСЗ можно представить в виде

где ${{{\mathbf{p}}}^{{\text{C}}}}$–$2N$-мерный вектор модельных представлений наблюдений: $N$ – число моментов наблюдений; $\alpha _{i}^{{\text{C}}}$ и $\delta _{i}^{{\text{C}}}$ $(i = 1, \ldots ,N)$ – прямое восхождение и склонение спутника: топоцентрические угловые координаты относительно экватора на эпоху J2000; ${\mathbf{q}} = {{({{{\mathbf{x}}}_{0}},{{{\mathbf{\dot {x}}}}_{0}},\gamma )}^{T}}$ – 7-мерный вектор модельных параметров: ${{{\mathbf{x}}}_{0}}$ и ${{{\mathbf{\dot {x}}}}_{0}}$ – векторы положения и скорости соответственно на начальную эпоху ${{t}_{0}}$; $\gamma = {A \mathord{\left/ {\vphantom {A m}} \right. \kern-0em} m}$ – коэффициент парусности, представляющий собой отношение площади миделевого сечения $A$ спутника к его массе m.

Угловые координаты на каждый момент наблюдения вычисляются из топоцентрического вектора положения спутника ${{{\mathbf{x}}}_{{\text{T}}}} = {{({{x}_{{{\text{T}}1}}},{{x}_{{{\text{T}}2}}},{{x}_{{{\text{T}}3}}})}^{T}}$ как

В свою очередь, топоцентрический вектор определяется как разность геоцентрических векторов положения спутника ${\mathbf{x}}$ и наблюдателя ${{{\mathbf{x}}}_{{\text{O}}}}$: ${{{\mathbf{x}}}_{{\text{T}}}} = {\mathbf{x}} - {{{\mathbf{x}}}_{{\text{O}}}}$. Положение наблюдателя корректируется за подвижку земной коры вследствие тектонических смещений и приливных деформаций, вызванных, главным образом, притяжением Луны и Солнца. В представлении наблюдений также учитывается эффект запаздывания света.

Движение спутника относительно геоцентра формализуется дифференциальными уравнениями 2-го порядка:

которые интегрируются численно коллокационным методом 16-го порядка (Авдюшев, 2020) в компьютерной арифметике с двойной точностью. Здесь ${\mathbf{P}}$ – равнодействующая сил, которая включает в себя (согласно IERS Conventions 2010): притяжение Земли как протяженного тела (с учетом приливных деформаций) до гармоник 360-го порядка; притяжение Луны и Солнца, рассматриваемых как материальные точки; влияние светового давления в рамках сферической модели спутника.

Для приведения времени наблюдений UTC к эфемеридному времени ТТ ($t$), относительно которого формализуется движения ИСЗ (4), используется поправка

Поправки ${\text{UT1}} - {\text{UTC}}$ и ${\text{UT1}} - {\text{TAI}}$ ежедневно определяются в “Международной службе вращения Земли” (https://hpiers.obspm.fr/eop-pc/index.php), а 32.184 с – это разность между двумя эфемеридными временами ${\text{TAI}} - {\text{TT}}$.

Параметры модели ${\mathbf{q}}$ определяются из наблюдений спутника ${{{\mathbf{p}}}^{{\text{O}}}}$ в рамках задачи наименьших квадратов:

которая решается итерационно методом Гаусса–Ньютона (7): где ${\mathbf{p}}_{{\mathbf{q}}}^{'} = {{\partial {{{\mathbf{p}}}^{{\text{C}}}}} \mathord{\left/ {\vphantom {{\partial {{{\mathbf{p}}}^{{\text{C}}}}} {\partial {\mathbf{q}}}}} \right. \kern-0em} {\partial {\mathbf{q}}}}$ – матрица частных производных от модельных представлений наблюдений по параметрам для приближения ${{{\mathbf{q}}}_{k}}$; ${{{\mathbf{p}}}^{{\text{O}}}}$ – вектор наблюдений спутника: Для каждого i-го момента наблюдения частные производные вычисляются как

Частные производные ${{\partial \alpha _{{}}^{{\text{C}}}} \mathord{\left/ {\vphantom {{\partial \alpha _{{}}^{{\text{C}}}} {\partial {{{\mathbf{x}}}_{T}}}}} \right. \kern-0em} {\partial {{{\mathbf{x}}}_{T}}}}$ и ${{\partial \delta _{{}}^{{\text{C}}}} \mathord{\left/ {\vphantom {{\partial \delta _{{}}^{{\text{C}}}} {\partial {{{\mathbf{x}}}_{T}}}}} \right. \kern-0em} {\partial {{{\mathbf{x}}}_{T}}}}$ получаются путем прямого дифференцирования формул (4). Частные производные ${{\partial {\mathbf{x}}} \mathord{\left/ {\vphantom {{\partial {\mathbf{x}}} {\partial {\mathbf{q}}}}} \right. \kern-0em} {\partial {\mathbf{q}}}}$ определяются из дифференциальных уравнений

которые интегрируются численно совместно с уравнениями движения (4).

В результате обработки наблюдений формируется ковариационная матрица параметров ${\mathbf{q}},$ которая характеризует их неопределенность вследствие ошибок наблюдений:

где ${{\sigma }^{2}} = {{S({\mathbf{\hat {q}}})} \mathord{\left/ {\vphantom {{S({\mathbf{\hat {q}}})} {2N - 6}}} \right. \kern-0em} {2N - 6}}$ – несмещенная оценка дисперсии ошибок ($\sigma $ – среднеквадратическая ошибка).

Поскольку определяемые параметры ${\mathbf{q}}$ имеют разную размерность, целесообразно прибегать к их масштабированию (нормализации), что может улучшить обусловленность задачи и тем самым уменьшить вычислительные ошибки в поправках к параметрам в итерационной схеме Гаусса–Ньютона. Масштабирование предполагает умножение частных производных ${\mathbf{p}}_{{\mathbf{q}}}^{'}$ и поправок в итерационной схеме на диагональную матрицу $\operatorname{diag} ({{s}_{1}}, \ldots ,{{s}_{7}}).$ Оптимальный вариант масштабирования – ${{s}_{i}} = {{q}_{i}}$ $(i = 1, \ldots ,7)$. Обращение матрицы осуществляется методом SVD, что позволяет следить за обусловленностью процесса улучшения орбит.

Начальное приближение для итерационной схемы (7) определяется нетрадиционным способом (Батурин, Чувашов, 2006). Из двух наблюдений на соседние близкие моменты времени ${{t}_{i}}$ и ${{t}_{{i{\kern 1pt} + {\kern 1pt} 1}}}$ вычисляются топоцентрические векторы положения ${{({{{\mathbf{x}}}_{{\text{T}}}})}_{i}}$ и ${{({{{\mathbf{x}}}_{{\text{T}}}})}_{{i{\kern 1pt} + {\kern 1pt} 1}}}$ по формулам

где топоцентрическое расстояние ИСЗ $\left| {{{{\mathbf{x}}}_{{\text{T}}}}} \right|$ задается приближенно из сторонних предположений (см., например, Tommei и др., 2007). Тогда топоцентрический вектор скорости ${{{\mathbf{\dot {x}}}}_{{\text{T}}}}$ оценивается как и в качестве начального приближения динамического состояния в эпоху ${{t}_{0}} = {{t}_{i}}$ выбираем

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

Новая версия программного комплекса “Численная модель движения ИСЗ” написана на процедурном высокоуровневом языке программирования Fortran. Интерфейсом программного комплекса служит входной текстовый файл. В этом файле задаются следующие данные: положение наблюдателя; начальная эпоха; параметры движения объекта (при наличии); коэффициенты масштабирования определяемых параметров; выборка наблюдений (время наблюдения, прямое восхождение, склонение, блеск); параметры интегратора (шаг, порядок, доверительная точность); переключатели учета возмущающих факторов (порядок и степень гармоник геопотенциала и селенопотенциала, Луна, Солнце, планеты, световое давление, релятивистские эффекты, приливы, атмосфера); параметры спутника (масса, площадь миделевого сечения, коэффициенты лобового сопротивления и отражения).

Таким образом, новый программный комплекс существенно отличается от своего предшественника, он позволяет не только прогнозировать движение, но и определять параметры движения и сил из обработки высокоточных позиционных наблюдений. Кроме того, в нем уточнены модели сил в соответствии с (IERS Conventions 2010), использован новый более эффективный интегратор (Авдюшев, 2020), который позволяет работать с системой уравнений, состоящей из уравнений первого и второго порядка, что удобно при совместном интегрировании уравнений движения, уравнений в вариациях и уравнений параметров MEGNO, введено точное вычисление изохронных производных путем численного интегрирования соответствующих уравнений. В новом программном комплексе совместно с параметрами орбиты определяется коэффициент парусности объекта.

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

РЕЗУЛЬТАТЫ ОПРЕДЕЛЕНИЯ ПАРАМЕТРОВ ДВИЖЕНИЯ ПО ДАННЫМ ИЗМЕРЕНИЙ

В качестве объектов наблюдения была выбрана группа, состоящая из четырех неуправляемых геосинхронных спутников с номерами 90 008, 90 031, 90 214 и 97 149. Форма представления наблюдений показана в табл. 1.

С помощью описанной выше численной модели движения систем ИСЗ, содержащей блоки вычисления изохронных производных и определения параметров движения по данным измерений, было проведено определение и уточнение орбит перечисленных выше объектов. В процессе численного моделирования движения учитывались возмущения от гармоник геопотенциала до 20-го порядка и степени, возмущения от Луны и Солнца, рассматриваемых как материальные точки, а также световое давление. В табл. 2 приведены данные о процессе улучшения орбит по наблюдениям, полученным в сентябре 2020 г.

Далее в табл. 3–6 приведены координаты и скорости объектов, элементы их орбит на выбранные эпохи, а также средние квадратические ошибки процесса улучшения орбит и параметры парусности объектов, полученные по данным наблюдений. На рис. 1 и 2 показано изменение блеска объектов в процессе наблюдения и распределения невязок в сферических координатах.

Как показывают приведенные результаты, орбиты объектов 90 008, 90 031 и 97 149 определяются с высокой точностью, а точность определения орбиты объекта 90 214 недостаточна. В то же время методическая точность прогнозирования движения всех объектов высокая. Данные приведены в табл. 7.

Здесь в качестве оценок точности $\Delta {\mathbf{x}}$ даны максимальные расхождения численных решений на интервале времени 1 год, полученные при двух соседних порядках задаваемой точности интегрирования ${{e}_{{{\text{tol}}}}} = 10$ и ${{e}_{{{\text{tol}}}}} = 11$. В процессе интегрирования учитывался полный набор возмущающих сил.

Наблюдения двух объектов 90 008 и 90 031 были продолжены в сентябре 2021 г. Данные о результатах улучшения орбит по этим измерениям приведены в табл. 7.

Приведенные в табл. 8 данные говорят о том, что наблюдения в сентябре 2021 г. не дают той точности, которая была получена по наблюдениям в 2020 г. Тем не менее была предпринята попытка объединить наблюдения объектов, выполненные в 2020 и 2021 годах, которая для объекта 90 031 увенчалась успехом.

В табл. 10 введены следующие обозначения: N – количество наблюдений, ON – отбраковано наблюдений, σ" – в конце процесса улучшения, контр. σ", во всех строках, кроме первой, – контрольное значение, полученное по 11 наблюдениям в октябре 2021 г., не участвовавших в улучшении орбиты, Т – число обусловленности Тодда.

Результаты, полученные по новому ПО, сравнивались с результатами, полученными по разработанному ранее программному комплексу “Численная модель движения систем ИСЗ”. Модели показали хорошее согласие (до третьей значащей цифры) в невязках и среднеквадратических ошибках. Использование точных значений изохронных производных в новом ПО позволяет повысить точность сходимости итерационного процесса улучшения орбит на два порядка.

ИССЛЕДОВАНИЕ ОРБИТАЛЬНОЙ ЭВОЛЮЦИИ ОБЪЕКТОВ

Перейдем к рассмотрению долговременной орбитальной эволюции выбранных объектов. Для исследования орбитальной эволюции будем использовать разработанный пакет программ, дополненный интегрированием уравнений для вычисления текущего и усредненного параметров MEGNO (Cincotta, Simo, 2000; Cincotta и др., 2003). Эволюция во времени усредненного параметра MEGNO дает представление о характере движения объекта: для квазипериодических (регулярных) орбит параметр MEGNO осциллирует около 2, для таких орбит усредненное значение MEGNO всегда стремится к 2, а для устойчивых орбит типа гармонического осциллятора усредненное значение MEGNO равно нулю. При значении усредненного параметра MEGNO больше 2 имеет место хаотизация движения, что не позволяет точно прогнозировать эволюцию элементов орбиты. А для исследования резонансных характеристик объектов будем использовать методики и программы, разработанные ранее (Томилова и др., 2019; Александрова и др., 2020).

Представление о точности интегрирования на 10-летнем интервале времени дают графики, приведенные на рис. 3.

Оценка точности интегратора была проведена методом сравнения с эталонной орбитой. На графике представлены результаты для двух рассматриваемых в данной работе объектов. В программном комплексе для выбора шага интегрирования задается параметр LL, который связан с задаваемой точностью следующим образом: ${{\left\| {\mathbf{e}} \right\|}_{{{\text{tol}}}}} = {{10}^{{ - LL}}}$. В качестве эталонной орбиты была выбрана орбита, полученная интегрированием уравнений движения методом 16-го порядка и параметром LL, равным 10, на 128-битной разрядной сетке. Остальные расчеты были выполнены с помощью арифметики двойной точности. Прогнозирование движения осуществлялось на 10-летнем интервале времени. Как показывают оценки, субсантиметровую точность прогнозирования можно получить при значении параметра LL = 5. Дальнейшее увеличение параметра LL нецелесообразно на 64-битной разрядной сетке, поскольку приводит к ненужному дроблению шага и накоплению ошибки округления. Орбитальные характеристики объектов отличаются величинами больших полуосей. Как показали, приведенные ниже исследования, орбитальная эволюция объектов также очень различна: объект 90 008 захвачен в резонанс, и его движение обладает высокой устойчивостью, а объект 97 149 движется вне зоны устойчивого резонанса, имеет регулярное, но менее устойчивое движение.

Движение каждого из объектов было промоделировано численно на интервале времени 10 лет. В процессе моделирования учитывались возмущения от геопотенциала до V_15,15, притяжение Луны, Солнца, влияние светового давления и приливных деформаций в теле Земли. Кроме того, был проведен анализ резонансных возмущений по методике Р. Алана (Allan, 1967a; 1967b), уточненной в работах (Кузнецов и др., 2012; Томилова и др., 2019).

Полученные из наблюдений значения больших полуосей говорят, что объекты находятся в окрестности резонанса 1 : 1 со скоростью вращения Земли. Указанная выше методика в этом случае дает следующие формулы для резонансных соотношений и резонансных (критических) аргументов:

причем первые три компоненты совпадают.

Результаты исследования приведены на графиках ниже (рис. 4–7). Здесь в блоках (а) на единой временной шкале показаны изменения наклонения орбиты $i$, эксцентриситета $e$, большой полуоси $a$ и усредненного параметра MEGNO. В блоках (б) показана эволюция орбит объектов в плоскости $\left\{ {a,\,\,\lambda } \right\}$. В блоках (в) дана эволюция во времени долготы подспутниковой точки. На графиках, обозначенных (г), (д) и (е), показана эволюция во времени резонансных соотношений и критических аргументов компонент ${{\Phi }_{1}},$ ${{\Phi }_{4}},$ ${{\Phi }_{5}}$ мультиплета орбитального резонанса со скоростью вращения Земли. Влияние светового давления выделено красным цветом.

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

Общий анализ результатов, приведенных на рис. 4–7, позволяет сделать ряд выводов о динамике исследуемых объектов:

− все изучаемые объекты принадлежат области геосинхронного движения; полученные результаты хорошо согласуются с результатами (Breiter и др., 2005);

− основными возмущающими факторами в движении объектов являются: влияние гармоник геопотенциала, гравитационное влияние Луны и Солнца, а также влияние светового давления;

− движение объектов 90 008, 90 031 и 90 214 является исключительно устойчивым, близко к гармоническому осциллятору, а движение объекта 97 149 является регулярным;

− объекты 90 008, 90 031 и 90 214 захвачены в резонанс 1 : 1 со скоростью вращения Земли (резонансные соотношения проходят через нулевые значения, а критические аргументы устойчиво либрируют);

− объект 97 149 не является резонансным;

− фазовые портреты в плоскости $\left\{ {a,\,\,\lambda } \right\}$ показывают, что во вращающейся системе координат объекты 90 008, 90 031, 90 214 либрируют около устойчивой точки либрации $\lambda $ = 75°;

− в плоскости $\left\{ {a,\,\,\lambda } \right\}$ во вращающейся системе координат объект 97149 находится в области регулярного движения;

− вековые резонансы низких порядков в движении объектов отсутствуют.

ЗАКЛЮЧЕНИЕ

Таким образом, в настоящей работе представлены результаты численного моделирования движения группы геосинхронных спутников, параметры которых получены по результатам позиционных наблюдений, выполненных на уникальной научной установке Цейсс–2000 в ЦКП “Терскольская обсерватория” Института астрономии РАН.

Моделирование выполнялось с помощью модифицированного программного комплекса “Численная модель движения ИСЗ”. Представлено описание программно-математического комплекса, который позволяет решать следующие задачи:

− численно с высокой точностью моделировать движение околоземных объектов;

− выполнять определение и улучшение орбит околоземных объектов по данным позиционных наблюдений;

− производить отбраковку наблюдений с контролем обусловленности задачи;

− определять из наблюдений параметр парусности объекта;

− проводить полное исследование долговременной орбитальной эволюции объекта с выявлением резонансных характеристик его динамики;

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

Использование данного программно-математического комплекса при обработке наблюдений группы геосинхронных спутников, полученных на телескопе Цейсс–2000, полностью подтвердило работоспособность и высокие характеристики точности программного комплекса.

Анализ орбитальной эволюции исследованных объектов позволил разделить объекты на захваченные в резонанс и нерезонансные.

Объекты 90 008, 90 031 и 90 214 захвачены в резонанс 1 : 1 со скоростью вращения Земли, их резонансные соотношения проходят через нулевые значения, а критические аргументы устойчиво либрируют. Фазовые портреты в плоскости $\left\{ {a,\,\,\lambda } \right\}$ показывают, что во вращающейся системе координат объекты 90 008, 90 031, 90 214 либрируют около устойчивой точки либрации $\lambda $ = 75°.

Объект 97 149 не является резонансным, в плоскости $\left\{ {a,\,\,\lambda } \right\}$ во вращающейся системе координат объект находится в области регулярного движения.

Вековые резонансы низких порядков в движении объектов отсутствуют.

Работа выполнена в рамках государственного задания Министерства науки и высшего образования Российской Федерации (тема № FSWM-2020-0049).