Методы моделирования гравитационного поля малых небесных тел сложной формы на примере кометы 67Р/Чурюмова — Герасименко Текст научной статьи по специальности «Физика»

УДК 550.831.015

БОТ: 10.18698/2308-6033-2023-11-2319

Методы моделирования гравитационного поля малых небесных тел сложной формы на примере кометы 67Р/Чурюмова — Герасименко

© Чжоцзинь Ли, А.Н. Клишин, В.В. Корянов, Д.С. Колесникова МГТУ им. Н. Э. Баумана, Москва, 105005, Россия

Определение параметров движения космического аппарата вблизи поверхности небесных тел несферической формы является непростой задачей для исследования ввиду трудности описания гравитационного поля таких тел. Проанализированы различные подходы к моделированию гравитационного поля тел сложной формы: метод сферических функций, метод минимальной объемлющей сферы, метод минимального объемлющего эллипсоида и метод многогранников. Проведен сравнительный анализ используемых моделей и предложен подход, основанный на применении нейронной сети для оптимизации расчета путем аппроксимации данных, полученных методом многогранников. Изложены результаты моделирования для описания гравитационного поля короткопериодической кометы 67Р/Чурюмова — Герасименко и проиллюстрированы рисунками и таблицами. Предложены рекомендации по использованию того или иного метода в процессе моделирования классической задачи двух тел.

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

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

Цель исследования — разработка методов моделирования гравитационного поля кометы 67Р/Чурюмова — Герасименко для решения задачи посадки на поверхность кометы.

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

1) получены физические параметры и модель формы кометы 67Р;

2) исследованы методы построения минимальной охватывающей сферы и минимального охватывающего эллипсоида для моделирования гравитационного поля кометы 67Р;

3) исследован метод многогранников для моделирования гравитационного поля кометы 67Р;

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

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

Метод сферических гармонических функций [1-7] широко используется при моделировании гравитационного поля небесных тел. Однако из-за ограничений сходимости ряда сферических гармоник его применяют только для описания области вне исследуемого небесного тела. Этот недостаток очень важен при решении задач моделирования посадки на поверхность кометы, когда требуется высокая точность определения гравитационного потенциала вблизи поверхности.

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

Тем не менее, оба подхода не позволяют описывать тела сложной невыпуклой формы. Более того, используя их, невозможно определить, где находится исследуемая точка — внутри или вне поверхности небесного тела, что неприемлемо при решении классической задачи двух тел [3, 4].

Метод групп частиц заключается в том, что заполняется пространство, охватываемое поверхностью малого небесного тела, N малыми сферами «масконами», тем самым заменяя общее распределение массы небесного тела совокупностью отдельных частей. Данный метод отличается сравнительной простотой и отсутствием необходимости дополнительного анализа областей, для которых замена поверхностью вращения неприемлема.

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

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

Теоретическое описание используемых методов. Поскольку исследуемое гравитационное поле кометы 67P предназначено для моделирования орбитальной динамики и посадки на небесное тело, традиционные методы сферических гармонических функций и эллипсоидальных гармонических функций неприменимы. Анализу подвергаются метод многогранников и его адаптация к вычислениям с помощью нейронных сетей [10, 11]. Для верификации метода многогранников используется метод масштабирования к минимальному охватывающему эллипсоиду по критерию равенства массы или площади поверхности.

Используемая в данной работе модель кометы 67P содержит 25 987 вершин. Набор точек выборки получается с использованием метода интервальных значений, где поверхность исследуемого тела заменяется сеткой с шагом, соответствующим принятому интервалу. Далее составляется стандартное матричное уравнение эллипсоида f (x), охватывающего набор точек выборки:

f(x) = (x - c )т A (x - c) = 1, (1)

где c — вектор центра эллипсоида; A — матрица признаков эллипсоида 3^3.

В ходе решения уравнения (1) по матрице признаков возможно определить параметры эллипсоида и далее использовать их для аппроксимации напряженности гравитационного поля [2]:

GMm

ep =--

r

( 2 2 2 0/22,7.22,22^ ^ a2 + b2 + c2 3 (a xo + b Уо +c zo )

1---1—---

10r2 10r4

У

(2)

где Ер — напряженность гравитационного поля в точке Р; G — гравитационная постоянная; М — масса центрального тела; т — масса / 2 2 2

КА; г = у1 х0 + Уо + г0 — расстояние между точкой Р и притягиваю-

Чжоцзинь Ли, А.Н. Клишин, В.В. Корянов, Д.С. Колесникова щим центром; x0, y0, z0 — координаты точки P в декартовой системе

координат (СК); a, b, c — полуоси эллипсоида.

Если принять, что масса распределена равномерно внутри эллипсоида, получим следующее выражение [2]:

Ep = -bcGmpG |l, -, , (3)

^ a a)

где p — усредненная плотность эллипсоида.

Тогда для получения гравитационного потенциала возьмем объемный интеграл следующего вида в сферической системе координат (СК):

G (1, P, q) =

= ГЦ1 i Г 2sin Ф drd ф 9,

у](r sin ф cos 9 - x)2 + (pr sin фsin 9 - y)2 + (qr cos ф - z)2

(4)

- c 9

где p = —, q =--введенные параметры; 9, ф — параметры, опре-

a a

деляющие зенитный и азимутальный углы соответственно. Приближенное решение интеграла (4):

4^a 2raa3 (p2 + q2 +1) 2na5 (p2y2 + q2z2 + x2 )

G1(1, p, q) ~ —---—3-+--i-. (5)

3r 15r3 5r5

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

Введем правостороннюю декартову СК с базисными векторами 1, ^ к и определим в ней вектор г = Ах1 + Ау$ + Агк, который задает положение элементарной массы астероида | относительно

точечной массы 0(х,у,г) (рис. 1).

В этом случае Ах = £ - х; Ау = | - у; Лг = £ - г, а единичный век-

$ Ах $ Ау $ Аг $

тор г = — 1 н--] н--к задает направление векторного поля. При-

г г г

чем нетрудно заметить, что

r

^vx +y +z )j

div(r) = divl - \ = ^---J-i+ ^---J-j +

^ r) dx dy

д

V

X + у + Г )

д2 г

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

и (гр)=с ти ^=чи-1 =2 чя, ,

где а— средняя плотность исследуемого тела.

(6)

0(х,у,г)

Рис. 1. Обозначение элементов в заданном базисе

Далее, согласно теореме Остроградского — Гаусса [2], тройной интеграл от дивергенции векторного поля г по области V, ограниченной поверхностью £, равен потоку векторного поля г через замкнутую кусочно-гладкую поверхность £ в направлении внешней нормали. Тогда гравитационный потенциал можно записать следующим образом:

1Оа\\\у ЖуГёК = 1 П • Гё£. (7)

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

I Л Л Л I ЛЛ Л Л

и = -п • гё£ = - Оа V (£6 п/ • гё£ :

2 2 ^ л /

2 2 / еФ

=2Оа! (,п/•{г 1 ,

(8)

Г еФ

где / — произвольная грань; ф — множество всех граней многогранника; п / — внешняя нормаль к грани /.

Поскольку в полученном выражении (8) двойной интеграл имеет смысл потока векторного поля через определенную грань, то выражение п/ • г = Аг является постоянной величиной, что позволяет вынести его за знак интеграла:

1 ^ 2 ( п/ (г)^ = — ^ 2 п/ • г/ (— ^, (9)

2 /е Ф Vг У 2 /еф ■> г

где г/ — вектор, определяющий положение грани /.

Интегральный член ( — в формуле (9) представляет собой

л / г

гравитационный потенциал двумерной плоскости. Преобразуем этот поверхностный интеграл следующим образом:

( =( (—^-( ^^ (—0)

Л/ \ г г-3 Л / „3 4 '

Чт л/^ Г Г j

2 2 2 2 л Поскольку r = Ax + Ay + Az , а Az = n f • rf постоянно, то:

c ^=c fr^+r!-^]d^-Azis Ald^=

cfr cf I r 3 r 3 Wfr 3

, v л J (11)

= С f-^x+ d^-AzС ^d^.

•r3f ^dAx r dAy r J •r3 f г3

Преобразуем первое слагаемое с использованием формулы Грина к криволинейному интегралу по замкнутому контуру С, ограничивающему плоскость /, а во втором слагаемом проведем замену:

( = ш/. Тогда

+ d^ -Az С ^

•r3f ^dAx r dAy r J nf r

ф -^(AxdAy -AydAx) - n f • r f •

J C т

(12)

ш f.

В такой постановке выражение для ш/ представляет собой телесный угол, образованный проекцией грани на шар единичного радиуса (рис. 2).

Согласно соотношению между соседними гранями многогранника, показанному на рис. 3, получаем

ш ^ = 2аге1аи

Г •(г / х гк)

1 ( ^ Гк ) + Г (Гк • Г ) + Гк ( • Г )'

(13)

Гг/к + Г

где г, Гj, Гк — векторы, проведенные из точки О в каждую вершину грани многогранника.

сЬ> О

Единичный

шар ¿Ь< О

Рис. 2. Телесный угол, образованный гранью

0(х,у,г)

Рис. 3. К пояснению формулы (13)

Линейный интегральный член в формуле (12) можно выразить путем упрощения внешних нормалей ребер многогранника:

-(АхёАу - АуёАх) = П£ • г£ ]£ - ёБ,

(14)

где rf = Axc i + Ajc j + Azc k — вектор от точки O до точки, принадлежащей грани f (в дальнейшем индекс f опустим).

Г 1

Интегральный член I — в формуле (14) представляет собой

•С r

потенциал бесконечной нити. Согласно [2], этот потенциал может быть выражен через отрезки a и b от точки O до концов ребра e и длину ребра e:

f ^dS = ln

JC r

a + b + e

= L„.

(15)

JC r a + b - e

Подставим результаты преобразований (13), (15) в (12). Получим

(16)

U = 2 Ге " Ee ' Ге ' Le - 2 G° S Vf ' Ff ' Гf f ^

2 eeE 2 fe®

Ee = nA ^П12 J + nв ^n21J (A, B e Ф ;12,21e E );

т

F/ = n f n f,

где E — множество всех ребер многогранника; Ee — квадратичный

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

^ ^A ^ ^B E12 = n A n12 + nB n 21.

Рис. 4. Схематичное изображение определения вектора E

В соответствии с приведенной формулой (16) расчета потенциала гравитационного поля многогранника оценим распределение потенциальной гравитационной энергии астероида и ускорение свободного падения, которое может быть выражено как частная производная функции гравитационного потенциала по радиусу-вектору:

^ (X, у, г) = 1^4 (17)

5г( х, у, г)

Тогда ускорение свободного падения в пространственной точке О (х, у, г) определяется как

Р(г) = Уи(г) = Са2 ге • Ее • ге • Ье -Са 2 г/ • Р/ • г/ •ш/;

ееЕ /еФ

V(VU(г)) = Са2 Ее • ье - Са 2 Р/ •ш/ V2U(г) = -Са 2 ш/, (18)

ее Е / е Ф / е Ф

где У(Уи (г)) — матрица Гессе; V2U (г)— оператор Лапласа.

Положение пространственной точки О (х, у, г) относительно при-

2

тягивающего центра определяется в зависимости от значения V и (г):

V 2и (г) _

—С^= 2 ш/ =

Са ГеФ

4 л - точка находится с внутренней стороны многогранника; (19)

0 - точка находится с внешней стороны многогранника.

Характеристики кометы 67Р. В качестве объекта исследования рассматривается комета 67Р, наиболее известная как комета Чурю-мова — Герасименко. Основные характеристики кометы сведены в табл. 1 [12, 13].

Таблица 1

Характеристики кометы 67Р

Характеристика Значение Характеристика Значение

Размеры, км: большая часть 4,1x3,2x1,3 Средняя плотность, г/см3 0,47

меньшая часть 2,5x2,5x2,0 Альбедо 0,06

Масса, кг 8,8022 х1012 Площадь поверхности, м2 4,7782 х107

Модель фигуры кометы в виде многогранника, используемая в исследовании, содержит 25 987 вершин, которые образуют 51 970 треугольных плоскостей. Моделирование и триангуляция фигуры кометы осуществлялась в программном пакете МЛТЬЛВ. Результат изображен на рис. 5.

г, м

1500

1000 500 0

-500 -1000 -1500

-2000 -1000 0 1000 2000 мм

Рис. 5. Модель кометы 67Р

Моделирование гравитационного поля кометы 67Р. На первом этапе получим наиболее грубое приближение гравитационного поля исследуемой кометы — модель минимальной охватывающей сферы. Результирующая сфера имеет радиус 27 003 м (рис. 6).

Рис. 6. Минимальная охватывающая сфера

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

Методы моделирования гравитационного поля малых небесных тел., z, м

1500 1000 500 0

-500 -1000 -1500

W

La

2000 1000 0 -1000 -2000 х,м Рис. 7. Минимальный охватывающий эллипсоид

Матрица признаков и вектор центра полученного эллипсоида в процессе решения матричного уравнения (1) имеют следующий вид:

1,371 х 10"7 1,883 х10"9 1,0038 х10"8 "99,59"

A = 1,883 х10"9 2,531 х10"7 1,2269 х10"8 , c = 5,91

1,003 х10"8 1,226 х10"8 2,9682 х!0"7 39,83

Параметры полученного эллипсоида: a = 2,7066-10 м, b = = 1,8235-103 м, с = 1,8480-103 м.

На следующем этапе в программном комплексе MATLAB на основе математической модели метода многогранников, представленной выше, получены параметры гравитационного поля. Сначала для каждой вершины рассчитано значение напряженности гравитационного поля [14]. Отклонение полученного вектора от нормали к поверхности в каждой точке исследуемой выборки 0г- представлено на рис. 8.

ег, град 180

100 90

Рис. 8. Зависимость отклонения вектора напряженности гравитационного поля от нормали к поверхности в каждой точке

Диапазон значений отклонений составляет [90°...180°], где 90° означает, что направление ускорения гравитационного поля параллельно поверхности, а 180° означает, что направление ускорения гравитационного поля перпендикулярно поверхности и противоположно направлению нормали в данной точке.

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

и, м2/с2 -0,12

-0,14

-0,16

-0,18

-0,20

-0,22

-0,24

-0,26

Рис. 9. Распределение потенциальной энергии гравитационного поля, полученные методом многогранников

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

Сравнение методов моделирования гравитационного поля. Для того чтобы данные гравитационного поля, полученные разными методами, были сопоставимы, необходимо масштабировать модели сферы и эллипсоида из условий равенства массы или площади поверхности [14, 15].

В предположении однородности модели эллипсоида его массу 4 13

можно оценить как т = 3 шЬерср = 1,795 -10 кг. Следовательно, для масштабирования модели эллипсоида по критерию равенства массы

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

0 „ АарЬр + Ьрср + арср с г 1п7 2 5 = 4кр-= 5,6-107 м2,р = 1,6075.

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

Для сравнения изоповерхности (поверхности равного потенциала), полученные методом эллипсоида (при условии равенства массы) и методом многогранников в пределах сферы влияния кометы, представлены на рис. 10.

и, м2/с2

х10'

х10'

Рис. 10. Изоповерхности гравитационного потенциала, полученные методом многогранников и методом охватывающего эллипсоида

Значения абсолютной погрешности определения потенциала гравитационного поля методом охватывающего эллипсоида и методом многогранников на различном расстоянии от поверхности кометы в пределах сферы влияния кометы представлены в табл. 2.

Согласно данным, представленным на рис. 10 и в табл. 2, можно заметить, что на расстоянии более 10 000 м значения гравитационного потенциала, установленные разными методами, отличаются не более чем на 1 %. Это определяет границу использования метода эллипсоида при моделировании движения КА в пределах сферы влияния кометы. В дальнейшем исключим из рассмотрения область за пределами сферы данного радиуса.

Таблица 2

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

г, км и, м2/с2 Ди2 , м2/с2 Дих, м2/с2 ДиУ , м2/с2

Эллипсоид Многогранник

10 -0,0586 -0,0584 -0,0003 0,0009 -0,0002

20 -0,0293 -0,0293 -0,00006 0,0002 -0,00005

30 -0,0196 -0,0195 -0,00002 0,00008 -0,00002

40 -0,0147 -0,0146 -0,00001 0,00005 -0,00001

50 -0,0117 -0,0117 -0,000007 0,00003 -0,000008

60 -0,0098 -0,00978 -0,000004 0,00002 -0,000005

2/ 2 М/С

х10

-0,25

х10

Рис. 11. Изоповерхности гравитационного потенциала, полученные методом многогранников (а, в) и методом охватывающего эллипсоида (б, г) в изометрии (а-б) и в проекциях на плоскость ОХУ (в-г) на расстояниях до 10 км

Во внутренней области гравитационного воздействия (на расстоянии менее 10000 м) наблюдаются существенные различия в форме изоповерхностей (рис. 11). В табл. 3 для сравнения представлены

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

Таблица 3

Погрешность модели эллипсоида относительно модели многогранников в различных точках вблизи ее поверхности

г, км и, м2/с2 Ди2 , м2/с2 дих, м2/с2 диг , м2/с2

Эллипсоид Многогранник

2 -0,28258 -0,25569 -0,02689 0,06149 0,01150

4 -0,14608 -0,14204 -0,00403 0,01013 -0,00050

6 -0,09785 -0,09655 -0,00130 0,00331 -0,00047

8 -0,07348 -0,07288 -0,00060 0,00161 -0,00031

10 -0,05880 -0,05847 -0,00033 0,00094 -0,00021

Определение параметров нейронной сети. Поскольку моделирование гравитационного поля вблизи поверхности кометы методом многогранников характеризуется большим объемом и низкой скоростью вычислений, данный метод неэффективно использовать напрямую при моделировании траектории приземления на поверхность кометы. Поэтому предлагается использовать иные подходы к моделированию гравитационного поля вблизи кометы [10, 11, 16-22], обладающие высокими точностью и скоростью расчета.

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

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

Такая задача относится к классу регрессионных [11], для ее решения используется полносвязная нейронная сеть с методом обратного распространения ошибки. На вход нейронной сети подается позиционный вектор г = [х, у, г], а на выходе имеем вектор напряженности гравитационного поля а = [ах, а2, аз ] при соответствующем позиционном векторе. В исследовании выбрана следующая архитектура: шесть скрытых слоев, каждый из которых содержит по 256 нейронов.

В качестве нелинейной функции активации в данном исследовании используется функция КеЬи(х):

(х < 0);

ЯеШ( х) = тах( х,0) =

0

х (х > 0).

(20)

В исследовании вся выборка данных представляет собой набор из 5 000 000 точек в диапазоне сферы влияния кометы. Размеры такой гравитационной сферы влияния Я можно оценить как

Я0 = Яп

т

ц. тела

тСолнца

(21)

Для кометы 67Р по формуле (21) это значение составит Я0 = 68493,528 м.

Далее необходимо рассчитать значение напряженности гравитационного поля методом многогранников по формуле (18) для всей выборки и использовать для последующего расчета лишь те точки, которые находятся вне поверхности исследуемой кометы. Отфильтровать эти значения можно согласно условию (19). В итоге набор данных составляет 4 661 276 случайных точек, среди которых 4 600 000 точек выбираются в качестве обучающей выборки, а остальные представляют собой тестовую выборку (рис. 12).

4661276 случайных точек

у, м -5 х, м

Рис. 12. Обучающая выборка

2

Поскольку гравитационное ускорение вблизи кометы 67Р невелико (порядок величин составляет 10), при отсутствии нормировки влияние ошибок вычислений будет относительно существенным, что

повлияет на окончательный результат обучения сети. Поэтому перед обучением набор данных нормируется к интервалу [0, 1] минимаксным методом, согласно которому:

*

X = •

X - X„

X — X

Лшах Лшт

(22)

где xmin, хшах — минимальное и максимальное значения среди

исходных данных.

В качестве оптимизируемого функционала рассматривается ошибка определения величины гравитационного ускорения в заданной точке. Размер батча при обучении установлен равным 128, количество эпох обучения — 200. Для того чтобы обеспечить эффективность обучения, скорость обучения п устанавливается равной 0,001 первых 10 раундах, далее — равной 0,0001. Используемая метрика — среднеквадратичная ошибка. В качестве оптимизации выбрана модель стохастического градиентного спуска (SGD) с моментом (SGD+Momentum), которая отличается от классического метода стохастического градиентного спуска возможностью преодолевать «овраги» функции ошибки [11].

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

0,0007 -

0,0001

0 50 100 150 200 Количество эпох обучения

0 50 100 150 200 Количество эпох обучения

б

Рис. 13. Результаты обучения нейронной сети: изменение ошибки (а) и изменение точности (б) в зависимости от эпохи обучения

В конце обучения ошибка составила 0,000135 м2/с4, а точность теста достигла 98,51 % Эти показатели приведены до процедуры постобработки, которая позволила бы сделать качественный вывод о результате обучения. В ходе постобработки данные денормализуются и

средняя ошибка по определению величины гравитационного ускорения составляет 4,914^ 10-6 м/c2, относительная погрешность о = 9,06 %.

Обучение нейронной сети проводилось с использованием ресурсов графической карты NVIDIA GeForce GTX 1650 (объем памяти

4 Гб, версия драйвера 531.14, версия CUDA 12.1). Обучение состояло из 200 эпох, общее время обучения — 10 ч 40 мин 12 с.

Для анализа производительности модели нейронной сети и модели многогранников была проведена серия из 20 тестов, оценивающих время получения вектора гравитационного поля в случайных 60000 точках в пределах сферы влияния кометы. Для обеспечения сопоставимости результатов расчеты проводили с использованием ресурсов центрального процессора. При выполнении тестов методом многогранников результат составил 256,2 с, а с использованием нейронной сети — 0,88 с, что в 291 раз быстрее. Следует также отметить, что метод с использованием нейронной сети легко адаптируется к высокопроизводительному языку программирования CUDA, который в дальнейшем позволит помимо обучения выполнять реализацию обученной модели с помощью ресурсов графического процессора, что также улучшит производительность модели.

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

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

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

ЛИТЕРАТУРА

[1] Werner R., Scheeres D. Exterior gravitation of a polyhedron derived and compared with harmonic and mascon gravitation representations of Asteroid 4769 Castalia. Celestial Mechanics and Dynamical Astronomy, 1996, vol. 65, pp. 313-344. https://doi.org/10.1007/BF00053511

[2] Werner R. The gravitational potential of a homogeneous polyhedron or don't cut corners. Celestial Mechanics and Dynamical Astronomy, 1994, vol. 59, pp. 253-278. https://doi.org/10.1007/BF00692875

[3] Субботин М.Ф. Введение в теоретическую астрономию. Москва, Наука, 1968, 800 с.

[4] Дубошин Г.Н. Справочное руководство по небесной механике и астродинамике. Москва, Наука,1976, 234 с.

[5] Эльясберг П.Е. Введение в теорию полета искусственных спутников Земли. Москва, Едиториал УРСС, 2020, 542 с.

[6] Чеботарев Г.А. Аналитические и численные методы небесной механики. Москва, Наука, 1965, 367с.

[7] Емельянов Н.В. Практическая небесная механика. Москва, Физический факультет МГУ, 2018, 270 с.

[8] Rossi A., Marzari F., Farinella P. Orbital evolution around irregular bodies. Earth, Planets and Space, 1999, vol. 51, pp. 1173-1180. https://doi.org/10.1186/BF03351592

[9] Geissler P., Petit J., Durda D., Greenberg R., Bottke W., Nolan M., Moore J. Erosion and ejecta reaccretion on 243 ida and its moon. Icarus, 1996, vol. 120, no. 0042, pp. 140-157.

[10] Zhao Y., Yang H., Li S. On-board modeling of gravity fields of elongated asteroids using Hopfield neural networks. Astrodynamics, 2023, vol. 7, no. 1, pp. 101-114. https://doi.org/10.1007/s42064-022-0151-3

[11] Zhang A., Lipton Z., Mu Li, Smola J. Dive into Deep Learning. Amazon Science, 2023, 1181 p.

[12] Chao B.F., Rubincam D.P. The gravitational field of Phobos. Geophysical Research Letters, 1989, vol. 16, no. 8, pp. 859-862. https://doi.org/10.1029/GL016i008p00859

[13] Scheeres D.J. Dynamics about uniformly rotating triaxial ellipsoids: applications to Asteroids. Icarus, 1994, vol. 110, no. 2, pp. 225-238. https://doi.org/10.1006/icar.1994.1118

[14] Shang H., Wei B., Lu J. Recent advances in modeling gravity field of small bodies. Journal of Deep Space Exploration, 2022, vol. 9, no. 4, pp. 359-372. https://doi.org/10.15982/jissn.2096-9287.2022.20220074

[15] Rausenberger O. Lehrbuch der Analytischen Mechanik. Leipzig, B.G. Teubner, 1888.

[16] Romain G., Jean-Pierre B. Ellipsoidal Harmonic expansions of the gravitational potential: Theory and application. Celestial Mechanics and Dynamical Astronomy, 2001, vol. 79, no. 4, pp. 235-275. https://doi.org/10.1023/A:1017555515763

[17] Shi J., Ma Y., Liang H. Research of activity of Main Belt Comets 176P/LINEAR, 238P/Read and 288P/(300163) 2006 VW139. Sci Rep, 2019, vol. 9, no. 5492. https://doi.org/10.1038/s41598-019-41880-0

[18] Turgut B. Application of back propagation artificial neural networks for gravity field modelling. ActaMontanistica Slovaca, 2016, vol. 21, no. 3, pp. 200-207.

[19] Yang H., Bai X., Baoyin H. Rapid generation of time-optimal trajectories for asteroid landing via convex optimization. Journal of Guidance, Control, and Dynamics, 2017, vol. 40, pp. 628-641. https://doi.org/10.2514/LG002170

[20] Hobson E.W. The Theory of Spherical and Ellipsoidal Harmonics. Chelsea Publishing Company, 1965, 514 p.

[21] Cui H., Zhang Z., Yu M. Computing and analysis of gravity field of Eros433 using polyhedron model. Journal of Harbin Institute of Technology, 2012, vol. 44, pp. 17-22.

[22] Jiang Yu, Baoyin H. Orbital Mechanics near a Rotating Asteroid. Journal of Astrophysics and Astronomy, 2014, vol. 35, pp. 17-38. https://doi.org/10.1007/s12036-014-9259-z

Статья поступила в редакцию 25.09.2023

Ссылку на эту статью просим оформлять следующим образом: Чжоцзинь Ли, Клишин А.Н., Корянов В.В., Колесникова Д.С. Методы моделирования гравитационного поля малых небесных тел сложной формы на примере кометы 67Р/Чурюмова — Герасименко. Инженерный журнал: наука и инновации, 2023, вып. 11. http://dx.doi.org/10.18698/2308-6033-2023-11-2319

Чжоцзинь Ли — студент кафедры «Динамика и управление полетом ракет и космических аппаратов» МГТУ им. Н.Э. Баумана. e-mail: lzj114065264@outlook.com

Клишин Алексей Николаевич — канд. техн. наук, доцент кафедры «Динамика и управление полетом ракет и космических аппаратов» МГТУ им. Н.Э. Баумана; автор более 150 научных и научно-популярных работ в области баллистики, динамики полета, управления движением летательных аппаратов и истории оружия. e-mail: alkl@mail.ru

Корянов Всеволод Владимирович — канд. техн. наук, доцент, первый заместитель заведующего кафедрой «Динамика и управление полетом ракет и космических аппаратов» МГТУ им. Н.Э. Баумана; автор более 200 научных работ в области моделирования баллистики и динамики движения космических и спускаемых аппаратов. e-mail: vkoryanov@bmstu.ru

Колесникова Дарья Сергеевна — аспирант 1-го года обучения кафедры «Динамика и управление полетом ракет и космических аппаратов» МГТУ им. Н.Э. Баумана. e-mail: dskolesnikova@gmail.ru

Methods in simulating gravitational field of the small celestial bodies of complex shapes on the 67P/Churyumov—Gerasimenko comet example

© Zhojin Li, A.N. Klishin, V.V. Koryanov, D.S. Kolesnikova Bauman Moscow State Technical University, Moscow, 105005, Russia

Determining the spacecraft motion parameters in vicinity to surface of a celestial body of the non-spherical shape is being actively studied, and it is a complex task due to a problem in describing the gravitational field of such a body. The paper analyzes various approaches to simulating gravitational field of bodies of the complex shape including the methods of spherical functions, minimum ambient sphere, minimum ambient ellipsoid and the polyhedral method. The models used were comparatively analyzed, and an approach was proposed based on introducing a neural network to optimize calculation by approximating the data obtained by the polyhedron method. Simulation results are presented to describe the gravitational field of the 67P/Churyumov — Gerasimenko short-period comet, they are illustrated with figures and tables. The paper proposes recommendations to using one or another method in simulating the classical two-body problem.

Keywords: comet gravitational field, spherical functions method, polyhedral method, mascon method, neural networks

REFERENCES

[1] Werner R., Scheeres D. Exterior gravitation of a polyhedron derived and compared with harmonic and mascon gravitation representations of Asteroid 4769 Castalia. Celestial Mechanics and Dynamical Astronomy, 1996, vol. 65, pp. 313-344. https://doi.org/10.1007/BF00053511

[2] Werner R. The gravitational potential of a homogeneous polyhedron or don't cut corners. Celestial Mechanics and Dynamical Astronomy, 1994, vol. 59, pp. 253-278. https://doi.org/10.1007/BF00692875

[3] Subbotin M.F. Vvedenie v teoreticheskuyu astronomiyu [Introduction to theoretical astronomy]. Moscow, Nauka Publ., 1968, 800 p.

[4] Duboshin G.N. Spravochnoe rukovodstvo po nebesnoy mekhanike i astrodinamike [Reference guide in celestial mechanics and astrodynamics]. Moscow, Nauka Publ., 1976, 234 p.

[5] Elyasberg P. E. Vvedenie v teoriyu poleta iskusstvennykh sputnikov Zemli [Introduction to the theory of flight of the artificial Earth satellites]. Moscow, Editorial URSS Publ., 2020, 542 p.

[6] Chebotarev G. A. Analiticheskie i chislennye metody nebesnoy mekhaniki [Analytical and numerical methods of celestial mechanics]. Moscow, Nauka Publ., 1965, 367 p.

[7] Emelyanov N.V. Prakticheskaya nebesnaya mekhanika [Practical celestial mechanics]. Moscow, Fizicheskiy Fakultet MGU Publ., 2018, 270 p.

[8] Rossi A., Marzari F., Farinella, P. Orbital evolution around irregular bodies. Earth, Planets and Space, 1999, vol. 51, pp. 1173-1180. https://doi.org/10.1186/BF03351592

[9] Geissler P., Petit J., Durda D., Greenberg R., Bottke W., Nolan M., Moore J. Erosion and ejecta recreation on 243 ida and its moon. Icarus, 1996, vol. 120, no. 0042, pp. 140-157.

[10] Zhao Y., Yang H., Li S. On-board modeling of gravity fields of elongated asteroids using Hopfield neural networks. Astrodynamics, 2023, vol. 7, no. 1, pp. 101-114. https://doi.org/10.1007/s42064-022-0151-3

[11] Zhang A., Lipton Z., Mu Li, Smola J. Dive into Deep Learning. Amazon Science, 2023, 1181 p.

[12] Chao B.F., Rubincam D.P. The gravitational field of Phobos. Geophysical Research Letters, 1989, vol. 16, no. 8, pp. 859-862. https://doi.org/10.1029/GL016i008p00859

[13] Scheeres D.J. Dynamics about uniformly rotating triaxial ellipsoids: Applications to Asteroids. Icarus, 1994, vol. 110, no. 2, pp. 225-238. https://doi.org/10.1006/icar.1994.1118

[14] Shang H., Wei B., Lu J. Recent advances in modeling gravity field of small bodies. Journal of Deep Space Exploration, 2022, vol. 9, no. 4, pp. 359-372. https://doi.org/10.15982/jissn.2096-9287.2022.20220074

[15] Rausenberge, O. Lehrbuch der Analytischen Mechanik. Leipzig, B.G. Teubner, 1888.

[16] Romain G., Jean-Pierre B. Ellipsoidal harmonic expansions of the gravitational potential: Theory and application. Celestial Mechanics and Dynamical Astronomy, 2001, vol. 79, no. 4, pp. 235-275.

https://doi.org/10.1023M.: 1017555515763

[17] Shi J., Ma Y., Liang H. Research of activity of main belt comets 176P/LINEAR, 238P/Read and 288P/(300163) 2006 VW139. Sci Rep, 2019, vol. 9, no. 5492. https://doi.org/10.1038/s41598-019-41880-0

[18] Turgut B. Application of back propagation artificial neural networks for gravity field modelling. ActaMontanistica Slovaca, 2016, vol. 21, no. 3, pp. 200-207.

[19] Yang H., Bai X., Baoyin H. Rapid generation of time-optimal trajectories for asteroid landing via convex optimization. Journal of Guidance, Control, and Dynamics,, 2017, vol. 40, pp. 628-641. https://doi.org/10.2514/LG002170

[20] Hobson E.W. The Theory of Spherical and Ellipsoidal Harmonics. Chelsea Publishing Company, 1965, 514 p.

[21] Cui H., Zhang Z., Yu M. Computing and analysis of gravity field of Eros433 using polyhedron model. Journal of Harbin Institute of Technology, 2012, vol. 44, pp. 17-22.

[22] Jiang Yu, Baoyin H. Orbital Mechanics near a Rotating Asteroid. Journal of Astrophysics and Astronomy, 2014, vol. 35, pp. 17-38. https://doi.org/10.1007/s12036-014-9259-z

Zhojin Li, Student, Department of Dynamics and Control of the Rocket and Spacecraft Flight, Bauman Moscow State Technical University. e-mail: lzj114065264@outlook.com

Klishin A.N., Cand. Sc. (Eng.), Associate Professor, Department of Dynamics and Control of the Rocket and Spacecraft Flight, Bauman Moscow State Technical University; author of more than 150 scientific and popular works in ballistics, flight dynamics, aircraft flight control and weapons history. e-mail: alkl@mail.ru

Koryanov V.V., Cand. Sc. (Eng.), Associate Professor, First Deputy Head of the Department of Dynamics and Flight Control of Rockets and Spacecraft, Bauman Moscow State Technical University; author of more than 200 scientific papers in ballistics simulation and motion dynamics of the spacecraft and descent vehicles. e-mail: vkoryanov@bmstu.ru

Kolesnikova D.S., 1st year Postgraduate, Department of Dynamics and Flight Control of Rockets and Spacecraft, Bauman Moscow State Technical University. e-mail: dskolesnikova@gmail.com