Известия РАН. Теория и системы управления, 2020, № 2, стр. 143-155

ПРИМЕНЕНИЕ ТРЕХИМПУЛЬСНОЙ СХЕМЫ МАНЕВРИРОВАНИЯ ДЛЯ ВОЗВРАТА С ОКОЛОЛУННОЙ ОРБИТЫ К ТОЧКЕ ВХОДА В АТМОСФЕРУ ЗЕМЛИ

Н. М. Гаврикова a*, Ю. Ф. Голубев b**

a МГУ им. М.В. Ломоносова
Москва, Россия

b ИПМ им. М.В. Келдыша РАН
Москва, Россия

* E-mail: ng062974@gmail.com
** E-mail: golubev@keldysh.ru

Поступила в редакцию 11.06.2019
После доработки 29.07.2019
Принята к публикации 30.09.2019

Полный текст (PDF)

Аннотация

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

Введение. Задача выбора оптимальной (с точки зрения затрат топлива) схемы маневрирования для траекторий возврата к Земле с произвольной окололунной орбиты рассматривается с 1950-х гг., однако она до сих пор окончательно не решена. В случае если имеется возможность вывести космический аппарат (КА) на целевую орбиту ожидания, допустимо применение одноимпульсной схемы маневрирования [1]. Однако при возникновении нештатной ситуации затраты характеристической скорости в случае использования одного импульсного маневра будут значительно возрастать, что при условии ограниченного запаса топлива может привести к неблагоприятному исходу. Тогда проблема может быть решена при помощи многоимпульсной схемы маневрирования.

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

1. Системы координат. Из неинерциальных систем координат (СК) рассматривается Гринвичская СК (ГСК) – геоцентрическая СК, связанная с земным эллипсоидом, нулевым меридианом которой считается меридиан, проходящий через Гринвичскую обсерваторию [11]. В качестве инерциальных СК рассматриваются геоцентрическая СК EJ2000.0 и селеноцентрическая СК MJ2000.0 (оси этих систем параллельны осям СК J2000.0), а также Гринвичская инерциальная СК на момент времени $t$ (далее – ГИСК($t$)), представляющая собой ГСК, зафиксированную в этот момент времени.

Если некоторый радиус-вектор или вектор скорости задан в СК, связанной с Луной, он будет обозначаться как $\rho $ или ${\mathbf{u}}$ соответственно. Если же некоторый радиус-вектор или вектор скорости задан в СК, связанной с Землей, для него будет использоваться обозначение ${\mathbf{r}}$ или ${\mathbf{v}}$ соответственно.

2. Постановка задачи. Пусть в некоторый момент времени ${{t}_{0}}$ в СК MJ2000.0 задан вектор состояния КА $({{\rho }_{0}},{{{\mathbf{u}}}_{0}})$. Требуется определить последовательность импульсных маневров, выполняемых внутри СД Луны и приводящих КА в некоторую точку входа в атмосферу Земли при условии, что следующие параметры считаются заданными:

1) широта ${{\varphi }_{{{\text{вх}}}}}$ и долгота ${{\lambda }_{{{\text{вх}}}}}$ точки входа в атмосферу;

2) высота точки входа в атмосферу ${{h}_{{{\text{вх}}}}}$;

3) наклонение траектории возврата в точке входа в атмосферу ${{i}_{{{\text{вх}}}}}$;

4) высота условного перигея ${{h}_{\pi }}$.

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

В приближении ограниченной задачи трех тел уравнения движения КА в СК EJ2000.0 можно записать в следующем виде:

(2.1)
${\mathbf{\ddot {r}}} = {{{\mathbf{f}}}_{E}}({\mathbf{r}},t) + {{{\mathbf{f}}}_{L}}({\mathbf{r}},t),\quad t \in \left[ {{{t}_{0}},{{t}_{{{\text{вх}}}}}} \right],$
где ${{{\mathbf{f}}}_{E}}({\mathbf{r}},\,t)$ – ускорение, возникающее из-за действия силы гравитационного притяжения Земли, ${{{\mathbf{f}}}_{L}}({\mathbf{r}},\,t)$ – ускорение, возникающее из-за действия силы гравитационного притяжения Луны, ${{t}_{{{\text{вх}}}}}$ – момент времени достижения точки входа в атмосферу.

Если задана последовательность $n$ импульсных маневров $\Delta {{{\mathbf{v}}}_{i}}({{t}_{i}})$, $i = \overline {1,n} $ , то уравнения (2.1) последовательно интегрируются на отрезках $\left[ {{{t}_{0}},{{t}_{1}}} \right], \ldots ,\left[ {{{t}_{{n - 1}}},{{t}_{n}}} \right],\left[ {{{t}_{n}},{{t}_{{{\text{вх}}}}}} \right]$ при следующих начальных условиях:

(2.2)
$\begin{gathered} {\mathbf{r}}({{t}_{0}}) = {{{\mathbf{r}}}_{0}},\quad {\mathbf{\dot {r}}}({{t}_{0}}) = {{{\mathbf{v}}}_{0}}, \\ {\mathbf{r}}({{t}_{i}}) = {{{\mathbf{r}}}_{i}},\quad {\mathbf{\dot {r}}}({{t}_{i}}) = {{{\mathbf{v}}}_{i}} + \Delta {{{\mathbf{v}}}_{i}},\quad i = \overline {1,n} , \\ \end{gathered} $
где $\left( {{{{\mathbf{r}}}_{i}},{{{\mathbf{v}}}_{i}}} \right)$ – вектор состояния непосредственно перед совершением i-го маневра.

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

Первый этап – вычисление вектора состояния $({{{\mathbf{r}}}_{\pi }},{{{\mathbf{v}}}_{\pi }})$ в условном перигее и времени перелета $\Delta t$ от момента ${{t}_{0}}$ старта с орбиты искусственного спутника Луны до достижения условного перигея. На этом этапе используется решение задачи Ламберта в постановке Бэттина, которое подробно рассмотрено в следующем пункте.

Второй этап – определение однопараметрического семейства гипербол возврата к Земле, обеспечивающих достижение краевых условий в точке входа в атмосферу Земли (параметром является время попадания на границу СД Луны). Каждая гипербола возврата задается селеноцентрическим вектором состояния на момент времени ее пересечения с СД Луны. Необходимость определения такого семейства гипербол обусловлена следующим. В предлагаемом алгоритме по фиксированному времени старта с окололунной орбиты t0 вычисляется нулевое приближение целевого вектора состояния $({{\rho }_{{{\text{СД}}\,\,{\text{0}}}}},{{{\mathbf{u}}}_{{{\text{СД}}\,\,{\text{0}}}}})$, соответствующее моменту достижения границы СД Луны ${{t}_{{{\text{СД}}\,\,{\text{0}}}}}$. Затем строится первое приближение для маневров $\left\{ {\Delta {{{\mathbf{u}}}_{i}},{{t}_{i}}} \right\}$, $i = \overline {1,n} $. После того как эти маневры будут учтены, момент времени достижения границы СД Луны ${{t}_{{{\text{СД}}\,\,{\text{1}}}}}$ и получившийся вектор состояния $({{\rho }_{{{\text{СД}}\,\,{\text{1}}}}},{{{\mathbf{u}}}_{{{\text{СД}}\,\,{\text{1}}}}})$ уже не будут совпадать со значениями, выбранными для нулевого приближения, вследствие чего схему маневрирования необходимо пересчитывать.

Третий этап – вычисление начального приближения последовательности импульсных маневров и соответствующей им целевой гиперболы возврата, определяемой своим вектором состояния $({{\rho }_{{{\text{СД}}}}},{{{\mathbf{u}}}_{{{\text{СД}}}}})$.

Последний этап – окончательное уточнение траектории возврата, т.е. вычисление последовательности импульсных маневров $\Delta {{{\mathbf{v}}}_{i}}({{t}_{i}})$, $i = \overline {1,n} $, обеспечивающих с некоторой точностью ${{\left( {{{\varepsilon }_{\varphi }},{{\varepsilon }_{i}},{{\varepsilon }_{\pi }}} \right)}^{{\text{T}}}}$ достижение параметров ${{\left( {{{\varphi }_{{{\text{вх}}}}},{{i}_{{{\text{вх}}}}},{{h}_{\pi }}} \right)}^{{\text{T}}}}$ при последовательном интегрировании уравнений движения (2.1) с начальными условиями (2.2).

4. Задача Ламберта в постановке Бэттина. Пусть для некоторой эллиптической орбиты относительно притягивающего центра известны радиус ${{r}_{\pi }}$ ее перицентра, радиус некоторой точки на орбите r0, а также время перелета $\Delta t$ от перицентра до этой точки. Требуется найти параметры орбиты и истинную аномалию ${{\vartheta }_{0}}$ точки r0 при условии, что ${{\vartheta }_{0}} \leqslant 180^\circ $ (т.е. движение происходит в одной полуплоскости относительно линии апсид эллипса).

Решение задачи находится из уравнения (4.1), например, при помощи метода Ньютона:

(4.1)
$\Delta t = \sqrt {\frac{{{{a}^{3}}}}{{{{\mu }_{E}}}}} \left\{ {\left( {\varepsilon - \sin \varepsilon - \pi } \right) - \left( {\delta - \sin \delta } \right)} \right\},$
(4.2)
$\varepsilon = 2\arcsin \sqrt {\frac{{{{r}_{0}} + {{r}_{\pi }} + с}}{{4a}}} ,\quad \delta = 2\arcsin \sqrt {\frac{{{{r}_{0}} + {{r}_{\pi }} - с}}{{4a}}} ,$
(4.3)
$с = \sqrt {r_{0}^{2} + r_{\pi }^{2} - 2{{r}_{0}}{{r}_{\pi }}\cos {{\vartheta }_{0}}} ,\quad a = \frac{{{{r}_{\pi }}({{r}_{\pi }} - {{r}_{0}}\cos {{\vartheta }_{0}})}}{{2{{r}_{\pi }} - {{r}_{0}}(1 + \cos {{\vartheta }_{0}})}}.$

В (4.3) представлены выражения для вычисления большой полуоси $a$ орбиты и хорды c (рис. 1), соединяющей два положения на орбите (перицентр и точку r0). Решение, вообще говоря, существует не для любых соотношений параметров ${{r}_{0}}$, ${{r}_{\pi }}$ и $\Delta t$. Укажем необходимые условия для существования решения, а также область, в которой его следует искать.

Рис. 1.

Задача Ламберта в постановке Бэттина

Пусть $r_{0}^{'} = 2a - {{r}_{0}}$ и $r_{\pi }^{'} = 2a - {{r}_{\pi }}$ – радиусы соответствующих положений относительно второго фокуса. Из неравенства $r_{0}^{'} + r_{\pi }^{'} \geqslant c$ следует

(4.4)
$a \geqslant \frac{{c + {{r}_{0}} + {{r}_{\pi }}}}{4},$
которое служит необходимым условием для существования эллиптической орбиты, проходящей через точки с радиусами ${{r}_{\pi }}$ и r0.

Найдем области, в которых может лежать решение ${{\vartheta }_{0}}$. Учитывая, что ${{\vartheta }_{0}} \in \left[ {0,\pi } \right]$ и ${{{{r}_{\pi }}} \mathord{\left/ {\vphantom {{{{r}_{\pi }}} {{{r}_{0}} \leqslant 1}}} \right. \kern-0em} {{{r}_{0}} \leqslant 1}}$, представим неравенство a > 0 в следующем виде (квадратная скобка обозначает выполнение хотя бы одного из условий):

(4.5)
$\left[ {\begin{array}{*{20}{l}} {\arccos \left( {2{{{{r}_{\pi }}} \mathord{\left/ {\vphantom {{{{r}_{\pi }}} {{{r}_{0}}}}} \right. \kern-0em} {{{r}_{0}}}} - 1} \right) \leqslant {{\vartheta }_{0}} \leqslant \pi ,} \\ {0 \leqslant {{\vartheta }_{0}} \leqslant \arccos \left( {{{{{r}_{\pi }}} \mathord{\left/ {\vphantom {{{{r}_{\pi }}} {{{r}_{0}}}}} \right. \kern-0em} {{{r}_{0}}}}} \right).} \end{array}} \right.$

При изменении истинной аномалии в интервале $0 \leqslant {{\vartheta }_{0}} \leqslant \arccos \left( {{{{{r}_{\pi }}} \mathord{\left/ {\vphantom {{{{r}_{\pi }}} {{{r}_{0}}}}} \right. \kern-0em} {{{r}_{0}}}}} \right)$ большая полуось будет изменяться в интервале $a \in \left[ {0,{{{{r}_{\pi }}} \mathord{\left/ {\vphantom {{{{r}_{\pi }}} 2}} \right. \kern-0em} 2}} \right]$, уменьшаясь от rπ/2 до 0. Однако условие (4.4) не выполняется на интервале $0 \leqslant {{\vartheta }_{0}} \leqslant \arccos \left( {{{{{r}_{\pi }}} \mathord{\left/ {\vphantom {{{{r}_{\pi }}} {{{r}_{0}}}}} \right. \kern-0em} {{{r}_{0}}}}} \right)$, так что этот интервал следует исключить из рассмотрения. Решение задачи лежит в интервале $\arccos \left( {2{{{{r}_{\pi }}} \mathord{\left/ {\vphantom {{{{r}_{\pi }}} {{{r}_{0}}}}} \right. \kern-0em} {{{r}_{0}}}} - 1} \right) \leqslant {{\vartheta }_{0}} \leqslant \pi $, в нем большая полуось уменьшается от +∞ до ${{\left( {{{r}_{\pi }} + {{r}_{0}}} \right)} \mathord{\left/ {\vphantom {{\left( {{{r}_{\pi }} + {{r}_{0}}} \right)} 2}} \right. \kern-0em} 2}$, а необходимое условие (4.4) существования орбиты верно для всех значений истинной аномалии.

Наконец, исследуем ограничения на время движения $\Delta t$. Функция $\Delta t({{\vartheta }_{0}})$ монотонно возрастает на полуинтервале $\left( {\arccos \left( {2{{{{r}_{\pi }}} \mathord{\left/ {\vphantom {{{{r}_{\pi }}} {{{r}_{0}}}}} \right. \kern-0em} {{{r}_{0}}}} - 1} \right),\pi } \right]$ (в этом можно убедиться, сделав несложные аналитические преобразования). При этом в точке ${{\vartheta }_{0}} = \arccos \left( {2{{{{r}_{\pi }}} \mathord{\left/ {\vphantom {{{{r}_{\pi }}} {{{r}_{0}}}}} \right. \kern-0em} {{{r}_{0}}}} - 1} \right)$ величина большой полуоси стремится к +∞, и эту точку следует исключить из рассмотрения.

Ограничения на время перелета определяются в соответствии с формулами:

(4.6)
$\begin{gathered} \Delta {{t}_{{\min }}} \leqslant \Delta t \leqslant \Delta {{t}_{{\max }}}, \\ \Delta {{t}_{{\min }}} = \Delta t\left( {\arccos \left( {2{{{{r}_{\pi }}} \mathord{\left/ {\vphantom {{{{r}_{\pi }}} {{{r}_{0}}}}} \right. \kern-0em} {{{r}_{0}}}} - 1} \right) + \nu } \right), \\ \Delta {{t}_{{\max }}} = \Delta t(\pi ) = \pi \sqrt {\frac{{{{{\left( {{{r}_{\pi }} + {{r}_{0}}} \right)}}^{3}}}}{{8\mu }}} , \\ \end{gathered} $
где $\nu $ – некоторый малый параметр (к примеру, $\nu = 0.001^\circ $), позволяющий избавиться от неопределенности функции $\Delta t({{\vartheta }_{0}})$ в точке ${{\vartheta }_{0}} = \arccos \left( {2{{{{r}_{\pi }}} \mathord{\left/ {\vphantom {{{{r}_{\pi }}} {{{r}_{0}}}}} \right. \kern-0em} {{{r}_{0}}}} - 1} \right)$.

Таким образом, решение задачи Ламберта в постановке Бэттина существует на полуинтервале $\arccos \left( {2{{{{r}_{\pi }}} \mathord{\left/ {\vphantom {{{{r}_{\pi }}} {{{r}_{0}}}}} \right. \kern-0em} {{{r}_{0}}}} - 1} \right) < {{\vartheta }_{0}} \leqslant \pi $, при этом для существования решения должно выполняться условие (4.6), а также, очевидно, должно быть ${{r}_{0}} \geqslant {{r}_{\pi }}$.

5. Определение параметров в условном перигее. Рассмотрим задачу двух тел (Земля, КА), в которой КА движется к Земле по эллиптической орбите. Обозначим через $\Delta t$ время движения КА от начального положения в момент времени ${{t}_{0}}$ до условного перигея. Пусть известны:

1) географические широта ${{\varphi }_{{{\text{вх}}}}}$ и долгота ${{\lambda }_{{{\text{вх}}}}}$ точки входа в атмосферу в ГИСК(${{t}_{0}} + \Delta t$), соответствующей моменту времени ${{t}_{0}} + \Delta t$;

2) наклонение ${{i}_{{{\text{вх}}}}}$ орбиты возврата к точке входа в ГИСК(${{t}_{0}} + \Delta t$);

3) высота условного перигея ${{h}_{\pi }}$;

4) момент времени начала движения ${{t}_{0}}$;

5) радиус-вектор ${{\rho }_{0}}$ КА в СК MJ2000.0 в момент времени ${{t}_{0}}$.

Требуется найти соответствующий вектор состояния $({{{\mathbf{r}}}_{\pi }},\,{{{\mathbf{v}}}_{\pi }})$ в условном перигее в СК EJ2000.0, а также оценить время перелета $\Delta t$ от момента ${{t}_{0}}$ до достижения условного перигея.

По радиус-вектору ${{\rho }_{0}}$ КА в СК MJ2000.0 в момент времени ${{t}_{0}}$ определяется радиус-вектор ${{{\mathbf{r}}}_{0}}$ КА в СК EJ2000.0 на этот же момент времени. Таким образом, орбита содержит как известное положение ${{{\mathbf{r}}}_{0}}$, так и пока неизвестное положение перицентра ${{{\mathbf{r}}}_{\pi }}$. К моменту ${{t}_{0}}$ центр Луны должен находиться приблизительно в апоцентре искомой эллиптической орбиты, а начальное положение ${{{\mathbf{r}}}_{0}}$ должно быть сравнительно близким к апоцентру. Учитывая, что движение будет происходить в одной полуплоскости эллипса относительно линии апсид, а начальное положение берется близким к апоцентру искомой орбиты перелета, то (как показывают предварительные расчеты) время перелета от ${{{\mathbf{r}}}_{0}}$ до ${{{\mathbf{r}}}_{\pi }}$ окажется в интервале $\Delta t \in [4,\,5.5]$ сут. Для уточнения значения времени перелета $\Delta t$ используется метод долготной привязки концов траектории перелета [9, 10]. В этом методе сравниваются долгота восходящего узла орбиты $\Omega _{{}}^{{{\text{ГИСК}}}}$ в ГИСК(${{t}_{0}} + \Delta t$) и долгота восходящего узла орбиты $\Omega _{{}}^{{{\text{EJ2000}}}}$ в СК EJ2000.0.

Когда для эллиптической орбиты определены долгота восходящего узла и наклонение, можно найти остальные параметры. Сначала необходимо вычислить истинную аномалию ${{\vartheta }_{0}}$ точки ${{{\mathbf{r}}}_{0}}$. Для этого решается задача Ламберта в постановке Бэттина: используются заданное расстояние ${{r}_{0}}$ от центра Земли, расстояние от центра Земли до условного перигея (т.е. до перицентра этой орбиты) ${{r}_{\pi }} = {{h}_{\pi }} + {{R}_{E}}$ и полученное выше время перелета $\Delta t$. По найденной истинной аномалии ${{\vartheta }_{0}}$, согласно (4.3), вычисляется большая полуось $a$. После вычисления аргумента широты и затем аргумента периценра орбиты все параметры орбиты становятся известными, поэтому можно найти вектор состояния в перицентре $({{{\mathbf{r}}}_{\pi }},\,{{{\mathbf{v}}}_{\pi }})$.

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

1) методом долготной привязки рассчитываются время перелета и два элемента орбиты (долгота восходящего узла и наклонение);

2) решается задача Ламберта в постановке Бэттина и по времени перелета рассчитывается истинная аномалия окололунного положения, эксцентриситет и большая полуось орбиты;

3) по всем полученным параметрам определяется вектор состояния в условном перигее.

Следует отметить, что не для любого момента времени старта ${{t}_{0}}$ существует эллиптическая орбита, обеспечивающая попадание в условный перигей. Зададим некоторое опорное время t и будем изменять момент времени старта в интервале $T = \left[ {t,t + 90\;{\text{сут}}} \right]$, выполняя первые два этапа описанной выше процедуры. Ширина интервала выбрана из соображения наглядности демонстрируемых результатов. Для каждого момента времени ${{t}_{0}} \in T$ вычисляется геоцентрическое начальное положение ${{{\mathbf{r}}}_{0}}({{t}_{0}})$, после чего рассчитывается время перелета $\Delta t$.

На рис. 2 точками показаны значения времени перелета, а также границы допустимого времени перелета (нижняя граница колеблется в пределах от 1.9 до 2.3 сут, верхняя граница – в пределах от 4.5 до 5.4 сут), определенные в соответствии с формулами (4.6). В случае если время перелета превышает максимально допустимое значение, задача Ламберта в постановке Бэттина не имеет решения.

Рис. 2.

Зависимость времени перелета от момента времени старта: сплошной линией обозначено максимальное время перелета, прерывистой – минимальное время перелета

6. Целевой вектор состояния на границе СД Луны. Предположим, что известны целевой вектор состояния в условном перигее $\left( {{{{\mathbf{r}}}_{\pi }},\,{{{\mathbf{v}}}_{\pi }}} \right) = \left( {{{{\mathbf{r}}}_{\pi }}\left( {{{t}_{0}}} \right),\,{{{\mathbf{v}}}_{\pi }}\left( {{{t}_{0}}} \right)} \right)$, момент времени старта ${{t}_{0}}$ и время перелета по эллиптической орбите, описанной в разд. 5, $\Delta t = \Delta t\left( {{{t}_{0}}} \right)$. Тогда КА должен оказаться в условном перигее в момент времени ${{t}_{0}} + \Delta t$. Пользуясь методом игнорирования возмущений, представим траекторию возврата состоящей из двух частей: внутри СД Луны и вне СД Луны. Относительно Луны часть траектории, расположенная внутри СД Луны, является ветвью гиперболы с притягивающим центром в центре Луны. Часть траектории, расположенная вне СД Луны, является дугой эллипса с притягивающим центром в центре Земли. При этом для указанной эллиптической траектории вектор состояния в перицентре равен $({{{\mathbf{r}}}_{\pi }},\,{{{\mathbf{v}}}_{\pi }})$, а момент времени прохождения перицентра равен ${{t}_{0}} + \Delta t$. Требуется определить вектор состояния на границе СД Луны описанной выше гиперболической траектории (гиперболы возврата).

Проинтегрируем траекторию в обратном времени от начального момента ${{t}_{0}} + \Delta t$ и начального вектора состояния $({{{\mathbf{r}}}_{\pi }},\,{{{\mathbf{v}}}_{\pi }})$ до момента достижения границы СД Луны ${{t}_{{{\text{СД}}}}}$ (т.е. до выполнения условия ${{\rho }_{{{\text{СД}}}}}$ = 66 000 км [13]). Пусть в момент времени ${{t}_{{{\text{СД}}}}}$ КА имеет вектор состояния $({{{\mathbf{r}}}_{{{\text{СД}}}}},{{{\mathbf{v}}}_{{{\text{СД}}}}})$ в СК EJ2000.0. Вычислим соответствующий ему вектор состояния $({{\rho }_{{{\text{СД}}}}},{{{\mathbf{u}}}_{{{\text{СД}}}}})$ в СК MJ2000.0. Таким образом, для каждого момента времени ${{t}_{0}}$ будут определены вектор состояния на границе СД $\left( {{{\rho }_{{{\text{СД}}}}}\left( {{{t}_{0}}} \right),{{{\mathbf{u}}}_{{{\text{СД}}}}}\left( {{{t}_{0}}} \right)} \right)$ и момент времени попадания на границу СД ${{t}_{{{\text{СД}}}}}\left( {{{t}_{0}}} \right)$. Варьируя момент времени ${{t}_{0}}$ (при этом сохраняя селеноцентрический вектор состояния КА в начальный момент времени) в некотором интервале ${{t}_{0}} \in \left[ {t{\kern 1pt} ',t{\kern 1pt} '{\kern 1pt} '} \right]$, можно определить как зависимости $\left( {{{\rho }_{{{\text{СД}}}}}\left( {{{t}_{0}}} \right),{{{\mathbf{u}}}_{{{\text{СД}}}}}\left( {{{t}_{0}}} \right)} \right)$, ${{t}_{{{\text{СД}}}}}\left( {{{t}_{0}}} \right)$, так и зависимость $\left( {{{\rho }_{{{\text{СД}}}}}\left( {{{t}_{{{\text{СД}}}}}} \right),{{{\mathbf{u}}}_{{{\text{СД}}}}}\left( {{{t}_{{{\text{СД}}}}}} \right)} \right).$

Если КА оказался на границе СД Луны в какой-либо из моментов времени ${{t}_{{{\text{СД}}}}}\left( {{{t}_{0}}} \right)$, ${{t}_{0}} \in \left[ {t{\kern 1pt} ',t{\kern 1pt} '{\kern 1pt} '} \right]$, и в этот момент времени имеет вектор состояния $\left( {{{\rho }_{{{\text{СД}}}}}\left( {{{t}_{{{\text{СД}}}}}} \right),{{{\mathbf{u}}}_{{{\text{СД}}}}}\left( {{{t}_{{{\text{СД}}}}}} \right)} \right)$, то траектория КА удовлетворяет краевым условиям на правом конце.

7. Выбор начального приближения для импульсных маневров. Предположим, что известны зависимости ${{t}_{{{\text{СД}}}}}\left( {{{t}_{0}}} \right)$ и $\left( {{{\rho }_{{{\text{СД}}}}}\left( {{{t}_{{{\text{СД}}}}}} \right),{{{\mathbf{u}}}_{{{\text{СД}}}}}\left( {{{t}_{{{\text{СД}}}}}} \right)} \right)$ для ${{t}_{0}}$, принадлежащего некоторому отрезку $\left[ {t{\kern 1pt} ',t{\kern 1pt} '{\kern 1pt} '} \right]$. Вообще говоря, вектор состояния $\left( {{{\rho }_{{{\text{СД}}}}}\left( {{{t}_{{{\text{СД}}}}}} \right),{{{\mathbf{u}}}_{{{\text{СД}}}}}\left( {{{t}_{{{\text{СД}}}}}} \right)} \right)$ на границе СД Луны на момент времени ${{t}_{{{\text{СД}}}}}\left( {{{t}_{0}}} \right)$ может быть недостижим для одноимпульсной схемы. Действительно, если от этих вектора состояния и момента времени проинтегрировать траекторию в обратном времени, она может не пересечь орбиту ожидания. Поэтому для одноимпульсной схемы в качестве вектора целевых параметров на границе СД Луны рассматривается не вектор состояния в момент времени ${{t}_{{{\text{СД}}}}}$, а вектор ${\mathbf{y}}({{t}_{{{\text{СД}}}}}) = ({{{\mathbf{u}}}_{{{\text{СД}}}}}({{t}_{{{\text{СД}}}}}),{{{\mathbf{\hat {h}}}}_{f}}({{t}_{{{\text{СД}}}}}))$, где ${{{\mathbf{\hat {h}}}}_{f}}\left( {{{t}_{{{\text{СД}}}}}} \right)$ – направление кинетического момента гиперболы возврата, определенный вектором состояния на границе СД Луны.

Расчеты показывают, что подобный выбор вектора целевых параметров обеспечивает попадание в точку на границе СД Луны, близкую к положению ${{\rho }_{{{\text{СД}}}}}\left( {{{t}_{{{\text{СД}}}}}} \right)$. Для трехимпульсной схемы в качестве вектора целевых параметров рассматривается полный вектор состояния ${\mathbf{y}}({{t}_{{{\text{СД}}}}})$ = = $({{\rho }_{{{\text{СД}}}}}({{t}_{{{\text{СД}}}}}),{{{\mathbf{u}}}_{{{\text{СД}}}}}({{t}_{{{\text{СД}}}}}))$ на границе СД Луны в момент времени ${{t}_{{{\text{СД}}}}}$.

Зафиксируем момент времени старта ${{t}_{0}}$ и определим для него из соответствующих зависимостей разд. 6 время попадания на границу СД Луны ${{t}_{{{\text{СД 0}}}}} = {{t}_{{{\text{СД}}}}}\left( {{{t}_{0}}} \right)$ и целевой вектор ${{{\mathbf{y}}}_{0}} = {\mathbf{y}}\left( {{{t}_{{{\text{СД 0}}}}}} \right)$. Предположим также, что выбрана схема маневрирования, приведенная ниже в разд. 8, 9, т.е. по вектору ${{{\mathbf{y}}}_{0}}$ определена некоторая последовательность маневров $\left\{ {{{t}_{i}}\left( {{{{\mathbf{y}}}_{0}}} \right),\Delta {{{\mathbf{u}}}_{i}}\left( {{{{\mathbf{y}}}_{0}}} \right)} \right\}_{{i = 1}}^{n}$ в СД Луны. При интегрировании уравнений движения с учетом заданной схемы маневрирования $\left\{ {{{t}_{i}}\left( {{{{\mathbf{y}}}_{0}}} \right),\Delta {{{\mathbf{u}}}_{i}}\left( {{{{\mathbf{y}}}_{0}}} \right)} \right\}_{{i = 1}}^{n}$ до достижения границы СД Луны в нецентральном поле определяется момент времени ${{\tilde {t}}_{{{\text{СД}}}}}$ попадания КА на границу СД Луны. При этом, вообще говоря, окажется ${{\tilde {t}}_{{{\text{СД}}}}} \ne {{t}_{{{\text{СД 0}}}}}$. Следовательно, необходимо найти такой вектор ${{{\mathbf{y}}}_{*}} = {{{\mathbf{y}}}_{*}}(t_{{{\text{СД}}}}^{'})$, чтобы для определенной по нему последовательности маневров $\left\{ {{{t}_{i}}\left( {{{{\mathbf{y}}}_{*}}} \right),\Delta {{{\mathbf{u}}}_{i}}\left( {{{{\mathbf{y}}}_{*}}} \right)} \right\}_{{i = 1}}^{n}$, приводящей КА на границу СД Луны в момент времени $t_{{{\text{СД}}}}^{{''}}$, выполнялось условие $t_{{{\text{СД}}}}^{'} = t_{{{\text{СД}}}}^{{''}}$.

Процедура определения ${{{\mathbf{y}}}_{*}}$ осуществляется посредством минимизации невязки $\Delta t = \Delta t(t_{{{\text{СД}}}}^{*}) = t_{{{\text{СД}}}}^{*} - {{t}_{{{\text{СД}}}}}(\{ {{t}_{i}}({{{\mathbf{y}}}_{*}}),\Delta {{{\mathbf{u}}}_{i}}({{{\mathbf{y}}}_{*}})\} _{{i = 1}}^{n})$ (зависимость от момента времени старта ${{t}_{0}}$ и начального вектора состояния опускается, так как эти параметры зафиксированы), где ${{t}_{{{\text{СД}}}}}(\{ {{t}_{i}},\Delta {{{\mathbf{u}}}_{i}}\} _{{i = 1}}^{n})$ – момент времени, в который КА оказывается на границе СД Луны с использованием схемы $\left\{ {{{t}_{i}},\Delta {{{\mathbf{u}}}_{i}}} \right\}_{{i = 1}}^{n}$, а ${{{\mathbf{y}}}_{*}} = {\mathbf{y}}(t_{{{\text{СД}}}}^{*})$. Задачу минимизации подобной невязки можно представить как задачу нелинейного программирования (подробно описана в разд. 10), в которой нелинейной функцией, подлежащей минимизации, будет $F(t_{{{\text{СД}}}}^{*}) = {\text{|}}\Delta t(t_{{{\text{СД}}}}^{*}){\text{|}}$, $t_{{{\text{СД}}}}^{*} \in \mathbb{R}$, при этом должно выполняться ограничение $t_{{{\text{СД}}}}^{*} - {{t}_{0}} > 0$.

Таким образом, определяется целевой момент времени пересечения границы СД ${{t}_{{{\text{СД *}}}}}$ и вектор целевых параметров ${{{\mathbf{y}}}_{*}} = {\mathbf{y}}({{t}_{{{\text{СД }}*}}})$, для которого выбранная схема маневрирования (начальное приближение для последующих вычислений) в нецентральном поле Луны приводит КА на границу СД Луны в момент времени ${{t}_{{{\text{СД }}*}}}$ (с заданной точностью). Кроме того, будет найден целевой вектор состояния на границе СД $({{\rho }_{{{\text{СД }}*}}},{{{\mathbf{u}}}_{{{\text{СД }}*}}}) = ({{\rho }_{{{\text{СД}}}}}({{t}_{{{\text{СД }}*}}}),{{{\mathbf{u}}}_{{{\text{СД}}}}}({{t}_{{{\text{СД }}*}}}))$.

8. Одноимпульсная схема маневрирования. Пусть для выведения КА с орбиты ожидания на гиперболу возврата с заданным вектором целевых параметров ${\mathbf{y}} = ({{{\mathbf{u}}}_{{{\text{СД}}}}},{{{\mathbf{\hat {h}}}}_{f}})$ используется один импульсный маневр.

Известный вектор состояния $({{\rho }_{0}},{{{\mathbf{u}}}_{0}})$ полностью определяет все параметры орбиты ожидания. Приближенно вектор асимптотической скорости ${\mathbf{u}}_{\infty }^{ + }$ гиперболы возврата [13] можно найти по формуле

(8.1)
${\mathbf{u}}_{\infty }^{ + } \approx {\text{u}}_{\infty }^{ + }\frac{{{{{\mathbf{u}}}_{{{\text{СД}}}}}}}{{\left\| {{{{\mathbf{u}}}_{{{\text{СД}}}}}} \right\|}},\quad {\text{u}}_{\infty }^{ + } = \sqrt {{\text{u}}_{{{\text{СД}}}}^{2} - 2\frac{{{{\mu }_{M}}}}{{{{\rho }_{{{\text{СД}}}}}}}} ,$
где ${{\mu }_{M}}$ – гравитационный параметр Луны. Большая полуось гиперболы возврата выражается равенством ${{a}_{f}} = {{{{\mu }_{M}}} \mathord{\left/ {\vphantom {{{{\mu }_{M}}} {{{{({\text{u}}_{\infty }^{ + })}}^{2}}}}} \right. \kern-0em} {{{{({\text{u}}_{\infty }^{ + })}}^{2}}}}$.

Определим истинную аномалию места совершения маневра на орбите ожидания. Вообще говоря, плоскость гиперболы возврата может пересекать орбиту ожидания в двух точках. Единичный вектор направления радиус-вектора положения для совершения маневра вычисляется по формуле ${{\hat {\rho }}_{1}} = \pm {{{{{{\mathbf{\hat {h}}}}}_{0}} \times {{{{\mathbf{\hat {h}}}}}_{f}}} \mathord{\left/ {\vphantom {{{{{{\mathbf{\hat {h}}}}}_{0}} \times {{{{\mathbf{\hat {h}}}}}_{f}}} {\left\| {{{{{\mathbf{\hat {h}}}}}_{0}} \times {{{{\mathbf{\hat {h}}}}}_{f}}} \right\|}}} \right. \kern-0em} {\left\| {{{{{\mathbf{\hat {h}}}}}_{0}} \times {{{{\mathbf{\hat {h}}}}}_{f}}} \right\|}}$, где ${{{\mathbf{\hat {h}}}}_{0}} = {{{{\rho }_{0}} \times {{{\mathbf{u}}}_{0}}} \mathord{\left/ {\vphantom {{{{\rho }_{0}} \times {{{\mathbf{u}}}_{0}}} {\left\| {{{\rho }_{0}} \times {{{\mathbf{u}}}_{0}}} \right\|}}} \right. \kern-0em} {\left\| {{{\rho }_{0}} \times {{{\mathbf{u}}}_{0}}} \right\|}}$ – направление кинетического момента орбиты ожидания. Одно положение соответствует ситуации, когда после маневра КА будет совершать движение по гиперболе возврата, минуя ее перицентр, а второе – ситуации, когда КА не пройдет через перицентр. Дальнейшие вычисления следует производить для обоих вариантов. В случае если решение существует для обоих знаков, следует выбрать тот, для которого затраты характеристической скорости меньше.

К началу маневра КА находится на орбите ожидания, поэтому истинная аномалия определяется как ${{\vartheta }_{{1,0}}} = \operatorname{sgn} ({{{\mathbf{\hat {h}}}}_{0}} \cdot \left[ {{{{\hat {\rho }}}_{{p,0}}} \times {{{\hat {\rho }}}_{1}}} \right])\arccos \left( {{{{\hat {\rho }}}_{{p,0}}} \cdot {{{\hat {\rho }}}_{1}}} \right)$, где ${{\hat {\rho }}_{{p,0}}}$ – направление на перицентр орбиты ожидания. Радиус-вектор ${{\rho }_{1}}$ места совершения маневра и вектор скорости ${{{\mathbf{u}}}_{{1,0}}}$ рассчитываются по известным элементам орбиты ожидания и служат уточнением для радиус-вектора ${{\rho }_{0}}$.

Осталось найти истинную аномалию места совершения маневра ${{\vartheta }_{{1f}}}$ на гиперболе возврата. Вычислим угол между вектором ${{\rho }_{1}}$ и вектором ${\mathbf{\hat {u}}}_{\infty }^{ + }$:

(8.2)
$\Delta \vartheta = \operatorname{sgn} ({{{\mathbf{\hat {h}}}}_{f}} \cdot [{{\hat {\rho }}_{1}} \times {\mathbf{\hat {u}}}_{\infty }^{ + }])\arccos ({{\hat {\rho }}_{1}} \cdot {\mathbf{\hat {u}}}_{\infty }^{ + }).$

Из уравнения орбиты

(8.3)
${{\rho }_{1}} = \frac{{{{a}_{f}}(1 - e_{f}^{2})}}{{1 + {{e}_{f}}\cos {{\vartheta }_{{1f}}}}}$
получается формула для расчета эксцентриситета гиперболы возврата:
(8.4)
${{e}_{f}} = \frac{{\sqrt {2{{k}^{2}} + {{{\sin }}^{2}}\left| {\Delta \vartheta } \right| + 2k\left( {1 - \cos \left| {\Delta \vartheta } \right|} \right) \pm \sin \left| {\Delta \vartheta } \right|\sqrt {{{{\sin }}^{2}}\left| {\Delta \vartheta } \right| + 4k\left( {1 - \cos \left| {\Delta \vartheta } \right|} \right)} } }}{{\sqrt 2 k}},$
где $k = {{a}_{f}}{\text{/}}{{\rho }_{1}}$. Истинная аномалия места совершения маневра на гиперболе возврата в случае, если $\operatorname{sgn} \left( {\Delta \vartheta } \right) = 1$, вычисляется как ${{\vartheta }_{{1f}}} = {{\vartheta }_{\infty }} - \Delta \vartheta $, а в случае, если $\operatorname{sgn} \left( {\Delta \vartheta } \right) = - 1$, – как ${{\vartheta }_{{1f}}} = {{\vartheta }_{\infty }} + \left| {\Delta \vartheta } \right|$ – 2π (здесь ${{\vartheta }_{\infty }} = \arccos \left( { - {1 \mathord{\left/ {\vphantom {1 {{{e}_{f}}}}} \right. \kern-0em} {{{e}_{f}}}}} \right)$ – предельный угол гиперболы возврата [13]).

После определения угла раствора гиперболы возврата $\delta = 2\arcsin \left( {{1 \mathord{\left/ {\vphantom {1 {{{e}_{f}}}}} \right. \kern-0em} {{{e}_{f}}}}} \right)$ можно найти направление вектора асимптотической скорости, соответствующего прилетной ветви гиперболы:

(8.5)
${\mathbf{\hat {u}}}_{\infty }^{ - } = \,{\mathbf{\hat {u}}}_{\infty }^{ + }\,\cos \delta + \,\,[{\mathbf{\hat {u}}}_{\infty }^{ + } \times {{{\mathbf{\hat {h}}}}_{f}}]\sin \delta ,$
а затем – направление на перицентр гиперболы возврата и направление скорости в перицентре:

(8.6)
${{\hat {\rho }}_{{pf}}} = \frac{{{\mathbf{\hat {u}}}_{\infty }^{ - } - {\mathbf{\hat {u}}}_{\infty }^{ + }}}{{\left\| {{\mathbf{\hat {u}}}_{\infty }^{ - } - {\mathbf{\hat {u}}}_{\infty }^{ + }} \right\|}},\quad {{{\mathbf{\hat {v}}}}_{{pf}}} = \frac{{{\mathbf{\hat {u}}}_{\infty }^{ - } + {\mathbf{\hat {u}}}_{\infty }^{ + }}}{{\left\| {{\mathbf{\hat {u}}}_{\infty }^{ - } + {\mathbf{\hat {u}}}_{\infty }^{ + }} \right\|}}.$

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

(8.7)
${\mathbf{u}}_{1}^{ + } = \sqrt {\frac{{{{\mu }_{M}}}}{{{{a}_{f}}(e_{f}^{2} - 1)}}} \left[ { - \sin {{\vartheta }_{{1,f}}}{{{\hat {\rho }}}_{{p,f}}} + \left( {{{e}_{f}} + \cos {{\vartheta }_{{1,f}}}} \right){{{{\mathbf{\hat {u}}}}}_{{p,f}}}} \right],$
а вектор скорости непосредственно перед совершением маневра есть ${\mathbf{u}}_{1}^{ - } = {{{\mathbf{u}}}_{{1,0}}}$, приращение скорости во время маневра определяется как $\Delta {{{\mathbf{u}}}_{1}} = {\mathbf{u}}_{1}^{ + } - {\mathbf{u}}_{1}^{ - }$. Так как истинная аномалия места совершения маневра известна, то путем решения уравнения Кеплера для орбиты ожидания определяется момент времени ${{t}_{1}}$ совершения маневра, уточняющий параметр ${{t}_{0}}$.

9. Трехимпульсная схема маневрирования. Рассмотрим теперь ситуацию, когда для выведения КА с орбиты ожидания на гиперболу возврата с заданным вектором целевых параметров ${\mathbf{y}} = \left( {{{\rho }_{{{\text{СД}}}}},{{{\mathbf{u}}}_{{{\text{СД}}}}}} \right)$ используется трехимпульсная схема маневрирования.

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

Первый маневр совершается в перицентре орбиты 0 и преобразует ее в эллиптическую орбиту, переводя КА в перицентр орбиты 1. Это делается для того, чтобы маневр изменял только величину скорости, не изменяя при этом ее направления. В случае, если орбита 0 близка к круговой, истинная аномалия ${{\vartheta }_{{1,0}}}$ места схода с орбиты ожидания варьируется, при этом после маневра КА по-прежнему оказывается в перицентре орбиты 1. Второй маневр осуществляется на линии пересечения плоскостей орбиты ожидания и гиперболы возврата и изменяет плоскость орбиты. Третий маневр совершается в общем перицентре орбит 2 и 3, т.е. радиус-векторы перицентров этих орбит совпадают: ${{\rho }_{{p,2}}} = {{\rho }_{{p,3}}}$.

Так как вектор целевых параметров ${\mathbf{y}} = \left( {{{\rho }_{{{\text{СД}}}}},{{{\mathbf{u}}}_{{{\text{СД}}}}}} \right)$ полностью определяет гиперболу возврата, а перицентр у орбит 2 и 3 общий, то остается определить радиус апоцентра ${{\rho }_{{a,1}}}$ орбиты 1 (геометрическое ограничение). Значение этого параметра ограничивается так, чтобы период орбиты 1 находился в диапазоне 1–2 сут. Так как второй маневр, изменяющий плоскость движения, совершается на орбите 1, то за счет увеличения расстояния до апоцентра орбиты 1 уменьшается модуль скорости на ней в месте совершения маневра, что приведет к уменьшению величины импульса при маневре. Изменяя радиус апоцентра орбиты 1 и определяя последовательность маневров по алгоритму, представленному ниже, вычисляется параметр ${{\rho }_{{a,1}}}$, обеспечивающий минимальные затраты характеристической скорости для рассматриваемого алгоритма построения схемы и для заданных начальных условий в центральном поле притяжения Луны. Если орбита ожидания круговая, то помимо параметра ${{\rho }_{{a,1}}}$ варьируется также истинная аномалия места схода с орбиты ожидания ${{\vartheta }_{{1,0}}}$.

Алгоритм построения трехимпульсной схемы выполняется следующим образом.

1. Определяется направление ${{{\mathbf{\hat {h}}}}_{f}}$ кинетического момента гиперболы возврата по вектору состояния ${\mathbf{y}}$.

2. По вектору состояния ${\mathbf{y}}$ рассчитываются элементы гиперболы возврата (например, по формулам из [13]) и вектор состояния $\left( {{{\rho }_{{p,3}}},{{{\mathbf{u}}}_{{p,3}}}} \right)$ в ее перицентре (например, по формулам из разд. 7).

3. Вычисляется первый импульсный маневр. Так как орбита ожидания полностью задана, а также определено место совершения первого маневра (либо перицентр орбиты ожидания, либо, если орбита ожидания круговая, некоторое положение ${{\vartheta }_{{1,0}}}$), то несложно вычислить направление ${{\hat {\rho }}_{1}} = {{\hat {\rho }}_{{p,0}}}$ радиус-вектора места совершения первого маневра, а также скорость на орбите ожидания непосредственно перед маневром ${\mathbf{u}}_{1}^{ - }$. Осталось вычислить скорость на орбите 1 непосредственно после маневра. Расстояние до перицентра орбиты 1 равно расстоянию до перицентра орбиты ожидания, т.е. ${{\rho }_{{p,1}}} = {{\rho }_{{p,0}}}$, расстояние ${{\rho }_{{a,1}}}$ до апоцентра орбиты 1 задано (геометрическое ограничение), большая полуось орбиты 1 определяется по формуле ${{a}_{1}} = {{\left( {{{\rho }_{{a,1}}} + {{\rho }_{{p,1}}}} \right)} \mathord{\left/ {\vphantom {{\left( {{{\rho }_{{a,1}}} + {{\rho }_{{p,1}}}} \right)} 2}} \right. \kern-0em} 2}$. Учитывая, что первый маневр переводит КА из перицентра орбиты ожидания в перицентр орбиты 1, найдем скорость после совершения первого маневра по формуле

(9.1)
${\mathbf{u}}_{1}^{ + } = \sqrt {{{\mu }_{M}}\left( {\frac{2}{{{{\rho }_{{p,1}}}}} - \frac{1}{{{{a}_{1}}}}} \right)} [{{{\mathbf{\hat {h}}}}_{0}} \times {{\hat {\rho }}_{{p,1}}}].$

Таким образом, первый импульсный маневр равен $\Delta {{{\mathbf{u}}}_{1}} = {\mathbf{u}}_{1}^{ + } - {\mathbf{u}}_{1}^{ - }$.

4. Второй импульсный маневр изменяет плоскость движения. Следовательно, он осуществляется на пересечении плоскостей гиперболы возврата и орбиты ожидания. Направление радиус-вектора места его совершения выражается формулой

(9.2)
${{\hat {\rho }}_{2}} = \pm \frac{{{{{{\mathbf{\hat {h}}}}}_{0}} \times {{{{\mathbf{\hat {h}}}}}_{f}}}}{{\left\| {{{{{\mathbf{\hat {h}}}}}_{0}} \times {{{{\mathbf{\hat {h}}}}}_{f}}} \right\|}},$
где ${{{\mathbf{\hat {h}}}}_{0}}$ – направление кинетического момента орбиты ожидания. Дальнейшие вычисления следует производить для обоих вариантов знака в формуле (9.2). В случае, если решение существует для обоих знаков, выбирается та схема маневрирования, для которой затраты характеристической скорости меньше. Истинная аномалия места совершения второго маневра на орбите 1 имеет вид
(9.3)
${{\vartheta }_{{2,1}}} = \operatorname{sgn} ({{{\mathbf{\hat {h}}}}_{0}} \cdot \left[ {{{{\hat {\rho }}}_{{p,1}}} \times {{{\hat {\rho }}}_{2}}} \right])\arccos \left( {{{{\hat {\rho }}}_{{p,1}}} \cdot {{{\hat {\rho }}}_{2}}} \right),$
эксцентриситет орбиты 1 задается формулой ${{e}_{1}} = {{\left( {{{\rho }_{{a,1}}} - {{\rho }_{{p,1}}}} \right)} \mathord{\left/ {\vphantom {{\left( {{{\rho }_{{a,1}}} - {{\rho }_{{p,1}}}} \right)} {\left( {{{\rho }_{{a,1}}} + {{\rho }_{{p,1}}}} \right)}}} \right. \kern-0em} {\left( {{{\rho }_{{a,1}}} + {{\rho }_{{p,1}}}} \right)}}$, после чего из уравнения орбиты вычисляется расстояние до места совершения второго маневра:
(9.4)
${{\rho }_{2}} = \frac{{{{a}_{1}}(1 - e_{1}^{2})}}{{1 + {{e}_{1}}\cos {{\vartheta }_{{2,1}}}}}$
и его радиус-вектор ${{\rho }_{2}} = {{\rho }_{2}}{{\hat {\rho }}_{2}}$.

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

(9.5)
${\mathbf{u}}_{2}^{ - } = \sqrt {{{\mu }_{M}}\left( {\frac{2}{{{{\rho }_{2}}}} - \frac{1}{{{{a}_{1}}}}} \right)} [{{{\mathbf{\hat {h}}}}_{0}} \times {{\hat {\rho }}_{2}}].$

Чтобы найти скорость непосредственно после совершения второго маневра, необходимо знать истинную аномалию ${{\vartheta }_{{2,2}}}$ места совершения второго маневра на орбите 2:

(9.6)
${{\vartheta }_{{2,2}}} = \operatorname{sgn} ({{{\mathbf{\hat {h}}}}_{f}} \cdot \left[ {{{{\hat {\rho }}}_{{p,2}}} \times {{{\hat {\rho }}}_{2}}} \right])\arccos \left( {{{{\hat {\rho }}}_{{p,2}}} \cdot {{{\hat {\rho }}}_{2}}} \right),$
где ${{\hat {\rho }}_{{p,2}}} = {{\hat {\rho }}_{{p,3}}}$ – направление на перицентр орбиты 2 (совпадает с направлением на перицентр гиперболы возврата). Эксцентриситет орбиты 2 определяется по формуле

(9.7)
${{e}_{2}} = \frac{{{{\rho }_{2}} - {{\rho }_{{p,2}}}}}{{{{\rho }_{{p,2}}} - {{\rho }_{2}}\cos {{\vartheta }_{{2,2}}}}}$.

Если ${{e}_{2}} > 1$ и ${{\vartheta }_{{2,2}}} > 0$, то решения для рассматриваемой схемы не существует, так как в этом случае КА совершает движение по гиперболической орбите 2, удаляясь от ее перицентра, а третий маневр должен производиться именно в ее перицентре. В этом случае следует выбрать другой знак в формуле (9.2) и заново произвести все вычисления.

Так как орбита 2 может быть как эллиптической, так и гиперболической, то определение дальнейших параметров зависит от величины эксцентриситета.

Если ${{e}_{2}} < 1$, то после нахождения большой полуоси орбиты ${{a}_{2}} = {{{{\rho }_{{p,2}}}} \mathord{\left/ {\vphantom {{{{\rho }_{{p,2}}}} {\left( {1 - {{e}_{2}}} \right)}}} \right. \kern-0em} {\left( {1 - {{e}_{2}}} \right)}}$ скорость непосредственно после совершения второго маневра выражается формулой

(9.8)
${\mathbf{u}}_{2}^{ + } = \sqrt {\frac{{{{\mu }_{M}}}}{{{{a}_{2}}(1 - e_{2}^{2})}}} \left( { - \sin {{\vartheta }_{{2,2}}}{{{\hat {\rho }}}_{{p,2}}} + \left( {{{e}_{2}} + \cos {{\vartheta }_{{2,2}}}} \right){{{{\mathbf{\hat {u}}}}}_{{p,2}}}} \right),$
где ${{{\mathbf{\hat {u}}}}_{{p,2}}} = {{{\mathbf{\hat {h}}}}_{f}} \times {{\hat {\rho }}_{{p,2}}}$ – направление вектора скорости в перицентре орбиты 2. Если ${{e}_{2}} > 1$, то после нахождения большой полуоси орбиты ${{a}_{2}} = {{{{\rho }_{{p,2}}}} \mathord{\left/ {\vphantom {{{{\rho }_{{p,2}}}} {\left( {{{e}_{2}} - 1} \right)}}} \right. \kern-0em} {\left( {{{e}_{2}} - 1} \right)}}$ скорость непосредственно после совершения второго маневра имеет вид

(9.9)
${\mathbf{u}}_{2}^{ + } = \sqrt {\frac{{{{\mu }_{M}}}}{{{{a}_{2}}(e_{2}^{2} - 1)}}} \left( { - \sin {{\vartheta }_{{2,2}}}{{{\hat {\rho }}}_{{p,2}}} + \left( {{{e}_{2}} + \cos {{\vartheta }_{{2,2}}}} \right){{{{\mathbf{\hat {u}}}}}_{{p,2}}}} \right).$

Таким образом, второй импульс равен $\Delta {{{\mathbf{u}}}_{2}} = {\mathbf{u}}_{2}^{ + } - {\mathbf{u}}_{2}^{ - }$.

5. Третий маневр осуществляется в перицентре орбиты 2 и переводит КА в перицентр орбиты 3. Поэтому скорость непосредственно перед импульсным маневром вычисляется по формуле

(9.10)
${\mathbf{u}}_{3}^{ - } = \sqrt {{{\mu }_{M}}\left( {\frac{2}{{{{\rho }_{{p,2}}}}} - \frac{1}{{{{a}_{2}}}}} \right)} {{{\mathbf{\hat {u}}}}_{{p,3}}},$
а скорость после совершения маневра определяется как ${\mathbf{u}}_{3}^{ + } = {{{\mathbf{u}}}_{{p,3}}}$. Таким образом, третий импульс равен $\Delta {{{\mathbf{u}}}_{3}} = {\mathbf{u}}_{3}^{ + } - {\mathbf{u}}_{3}^{ - }$.

Так как истинные аномалии мест совершения маневров известны, а все орбиты полностью определены, то путем решения соответствующих уравнений Кеплера вычисляются моменты времени ${{t}_{1}}$, ${{t}_{2}}$, ${{t}_{3}}$. При этом момент времени ${{t}_{1}}$ уточняет параметр ${{t}_{0}},$ а радиус-вектор ρ1 – параметр ρ0.

10. Решение краевой задачи. Рассмотрим краевую задачу, для которой заданы следующие краевые условия: вектор состояния КА $({{\rho }_{0}},{{{\mathbf{u}}}_{0}})$, заданный на начальный момент времени ${{t}_{0}}$, и условия на правом конце: достижение параметров ${{\left( {{{\varphi }_{{{\text{вх}}}}}\left( {{{h}_{{{\text{вх}}}}}} \right),\,{{i}_{{{\text{вх}}}}}\left( {{{h}_{{{\text{вх}}}}}} \right),{{h}_{\pi }}} \right)}^{{\text{T}}}}$. Требуется найти последовательность маневров $\left\{ {{{t}_{i}},\Delta {{{\mathbf{u}}}_{i}}} \right\}_{{i = 1}}^{n}$, которая обеспечивает достижение условий на правом конце при последовательном интегрировании уравнений (2.1) с начальными условиями (2.2).

В качестве начального приближения будем использовать решение $\{ t_{i}^{0},\Delta {\mathbf{u}}_{i}^{0}\} _{{i = 1}}^{n}$, полученное в разд. 7–9. Использование данного начального приближения с учетом нецентральности гравитационного поля Луны и нецентральности гравитационного поля Земли приводит КА на границу СД Луны с ошибкой ${{{\mathbf{c}}}_{{{\text{СД}}}}}$ по вектору ${\mathbf{y}}$ целевых параметров (в случае одноимпульсной схемы маневрирования ${{{\mathbf{c}}}_{{{\text{СД}}}}} = {{\left( {\left| {\Delta {{t}_{{{\text{СД}}}}}} \right|,\left\| {\Delta {{{\mathbf{u}}}_{{{\text{СД}}}}}} \right\|,\left\| {\Delta {{{{\mathbf{\hat {h}}}}}_{f}}} \right\|} \right)}^{{\text{T}}}}$, а в случае трехимпульсной схемы маневрирования cСД = ${{\left( {\left| {\Delta {{t}_{{{\text{СД}}}}}} \right|,\left\| {\Delta {{\rho }_{{{\text{СД}}}}}} \right\|,\,\left\| {\Delta {{{\mathbf{u}}}_{{{\text{СД}}}}}} \right\|} \right)}^{{\text{T}}}}$). Для достижения целевых параметров в точке входа в атмосферу необходимо устранить сначала указанные выше ошибки посредством коррекции имеющегося начального приближения последовательности маневров.

После минимизации невязок ${{{\mathbf{c}}}_{{{\text{СД}}}}}$ (процедура минимизации невязок описана ниже) и получения некоторой схемы $\{ t_{i}^{{{\text{СД}}}},\Delta {\mathbf{u}}_{i}^{{{\text{СД}}}}\} _{{i = 1}}^{n}$ гипербола возврата будет близка к целевой, однако краевые условия на правом конце все же будут выполняться с некоторой ошибкой, т.е. необходимо также минимизировать невязку ${{{\mathbf{c}}}_{{{\text{вх}}}}} = {{\left( {\left| {\Delta {{\varphi }_{{{\text{вх}}}}}\left( {{{h}_{{{\text{вх}}}}}} \right)} \right|,\left| {\Delta {{i}_{{{\text{вх}}}}}\left( {{{h}_{{{\text{вх}}}}}} \right)} \right|,\left| {\Delta {{h}_{\pi }}} \right|} \right)}^{{\text{T}}}}$. Для минимизации этой невязки в качестве начального приближения принимается схема $\{ t_{i}^{{{\text{СД}}}},\Delta {\mathbf{u}}_{i}^{{{\text{СД}}}}\} _{{i = 1}}^{n}$.

В случае одноимпульсной схемы параметрами для минимизации невязок ${{{\mathbf{c}}}_{{{\text{СД}}}}}$ и ${{{\mathbf{c}}}_{{{\text{вх}}}}}$ служат параметры схемы маневрирования $\left\{ {{{t}_{1}},\Delta {{{\mathbf{u}}}_{1}}} \right\}$ (т.е. четыре параметра), в случае же трехимпульсной схемы параметрами для минимизации будут $\left\{ {{{t}_{i}},\Delta {{{\mathbf{u}}}_{i}}} \right\}_{{i = 1}}^{3}$ (т.е. 12 параметров).

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

(10.1)
${\mathbf{x}} = {{({{t}_{1}},\Delta {\mathbf{u}}_{1}^{{\text{T}}})}^{{\text{T}}}}\quad {\text{или}}\quad {\mathbf{x}} = {{({{t}_{1}},\Delta {\mathbf{u}}_{1}^{{\text{T}}},{{t}_{2}},\Delta {\mathbf{u}}_{2}^{{\text{T}}},{{t}_{3}},\Delta {\mathbf{u}}_{3}^{{\text{T}}})}^{{\text{T}}}}.$

Вектор невязок зависит от схемы маневрирования: ${\mathbf{c}} = {\mathbf{c}}\left( {\mathbf{x}} \right)$. Выберем некоторый вектор положительных коэффициентов k той же размерности, что и вектор c, определим с его помощью нелинейную функцию $F\left( {\mathbf{x}} \right) = {{{\mathbf{k}}}^{{\text{T}}}}{\mathbf{c}}\left( {\mathbf{x}} \right)$.

Тогда задача минимизации вектора невязок c сводится к поиску минимума:

(10.2)
$\mathop {\min }\limits_{\mathbf{x}} F\left( {\mathbf{x}} \right),\quad {\mathbf{x}} \in {{\mathbb{R}}^{n}},\quad n = 4\quad {\text{или}}\quad n = 12,$
при этом должны выполняться s ограничений на компоненты вектора x:

(10.3)
${{g}_{j}}\left( {\mathbf{x}} \right) \geqslant 0,\quad {\mathbf{x}} \in {{\mathbb{R}}^{n}},\quad j = \overline {1,s} .$

Для одноимпульсной (или трехимпульсной) схемы маневрирования в терминах $\left\{ {{{t}_{1}},\Delta {{{\mathbf{u}}}_{1}}} \right\}$ (или $\left\{ {{{t}_{i}},\Delta {{{\mathbf{u}}}_{i}}} \right\}_{{i = 1}}^{3}$) целесообразно задать, например, следующие три (s = 3) или девять (s = 9) ограничений:

(10.4)
$\left\{ \begin{gathered} {{t}_{1}} - {{t}_{0}} \geqslant 0, \hfill \\ {{t}_{{\max }}} - {{t}_{1}} \geqslant 0, \hfill \\ {{u}_{{\max }}} - \left\| {\Delta {{{\mathbf{u}}}_{1}}} \right\| \geqslant 0 \hfill \\ \end{gathered} \right.\quad {\text{или}}\quad \left\{ \begin{gathered} {{t}_{i}} - {{t}_{{i - 1}}} \geqslant 0, \hfill \\ {{t}_{{\max ,i}}} - {{t}_{{i - 1}}} \geqslant 0, \hfill \\ {{u}_{{\max ,i}}} - \left\| {\Delta {{{\mathbf{u}}}_{i}}} \right\| \geqslant 0, \hfill \\ \end{gathered} \right.\quad i = 1,2,3.$

Параметры ${{t}_{{\max }}}$, ${{t}_{{\max ,i}}}$, ${{u}_{{\max }}}$, ${{u}_{{\max ,i}}}$ задаются эмпирически (например, ${{t}_{{\max }}} = 6$ сут, ${{t}_{{\max ,i}}} = 0.25$ сут, ${{u}_{{\max }}} = 3$ км/с, ${{u}_{{\max ,i}}} = 2$ км/с).

Для решения задачи нелинейного программирования (10.2) с ограничениями (10.4) был применен метод последовательного квадратичного программирования (SLSQP), который является одним из наиболее распространенных и эффективных методов решения такого класса задач. Этот метод подробно описан в статье [14].

11. Результаты моделирования. Рассмотрим в СК MJ2000.0 околокруговую орбиту ожидания высоты 100 км с долготой восходящего узла 20°, аргументом перицентра 0° и эксцентриситетом, равным 0.001 [8]. Наклонение орбиты ожидания ${{i}_{0}}$ в СК MJ2000.0 будет варьироваться. Пусть широта точки входа равна ${{\varphi }_{{{\text{вх}}}}} = - 7.5^\circ $, высота условного перигея ${{h}_{\pi }} = 51.7$ км и наклонение траектории возврата в точке входа в атмосферу ${{i}_{{{\text{вх}}}}} = 54.14^\circ $.

Важным параметром, оказывающим влияние на затраты характеристической скорости, является угол между плоскостями орбиты ожидания и гиперболы возврата $\alpha $. Результаты работы алгоритма с использованием одноимпульсной и трехимпульсной схем маневрирования для склонения Луны к экваториальной плоскости Земли в начальный момент времени 0.56° приведены в табл. 1. В таблице цветом выделены строки, соответствующие ситуациям, когда использование одноимпульсной схемы лучше с точки зрения затрат характеристической скорости. Этот эффект означает, что либо из-за введения ограничений не существует трехимпульсной схемы, имеющей меньшие затраты, чем одноимпульсная схема, либо рассматриваемая трехимпульсная схема не является оптимальной. Подобный результат является ожидаемым, так как в этих ситуациях величина угла $\alpha $ для одноимпульсной схемы очень мала, т.е. изменение плоскости движения практически не требуется.

Таблица 1
${{i}_{0}}$, град Одноимпульсная схема Трехимпульсная схема
$\Delta V$, км/с $\alpha $, град $\Delta t$, сут $\Delta V$, км/с $\alpha $, град $\Delta t$, сут
0 1.17 22.95 4.18 1.25 18.59 5.27
  (1.30) (0.01) (1.22) (0.50)
30 1.43 23.13 4.28 1.30 26.93 5.16
  (1.48) (0.03) (1.31) (0.39)
60 2.22 47.85 4.70 1.68 60.13 7.82
  (2.25) (0.04) (1.61) (0.60)
90 2.43 78.21 3.16 1.72 82.20 4.95
  (2.45) (0.08) (1.74) (0.34)
120 2.52 104.46 6.04 1.83 111.16 6.61
  (2.53) (0.04) (1.88) (0.33)
150 2.80 132.61 4.73 2.04 138.67 5.38
  (2.88) (0.03) (1.99) (0.72)
180 3.05 159.32 4.39 2.14 161.56 5.39
  (3.09) (0.07) (2.16) (0.53)

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

Таблица 2
${{i}_{0}}$, град ${{t}_{1}} - {{t}_{0}}$, сут $\Delta {{V}_{1}}$, км/с ${{t}_{2}} - {{t}_{0}}$, сут $\Delta {{V}_{2}}$, км/с ${{t}_{3}} - {{t}_{0}}$, сут $\Delta {{V}_{3}}$, км/с
0 0.013 0.410 0.154 0.376 0.495 0.437
0.012 0.396 0.157 0.408 0.499 0.455
30 0.031 0.413 0.172 0.429 0.392 0.472
0.030 0.416 0.170 0.414 0.392 0.474
60 0.080 0.308 0.143 0.818 0.600 0.486
0.078 0.299 0.145 0.852 0.601 0.532
90 0.035 0.413 0.175 0.863 0.343 0.463
0.035 0.415 0.176 0.834 0.342 0.476
120 0.035 0.412 0.174 0.999 0.336 0.470
0.037 0.418 0.178 0.902 0.327 0.512
150 0.075 0.415 0.200 1.089 0.720 0.493
0.072 0.427 0.192 1.068 0.721 0.547
180 0.065 0.415 0.194 1.207 0.530 0.546
0.064 0.410 0.198 1.175 0.534 0.558

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

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

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

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

Список литературы

  1. Боровин Г.К., Голубев Ю.Ф., Грушевский А.В., Заславский Г.С., Захваткин М.В., Корянов В.В., Лавренов С.М., Морской И.М., Симонов А.В., Степаньяц В.А., Тучин А.Г., Тучин Д.А., Ярошевский В.С. Баллистико-навигационное обеспечение полетов автоматических космических аппаратов к телам Солнечной системы. Химки: АО “НПО Лавочкина”, 2018.

  2. Bean W.C. Minimum ∆V, Three-impulse Transfer onto a Trans-mars Asymptotic Velocity Vector // NASA Technical Note TN D-5757, 1970. P. 1–26.

  3. Doll J.R., Gobetz F.W. A Survey of Impulsive Trajectories // AIAA Journal. 1969. V. 7. № 5. P. 801–834.

  4. Edelbaum T. How Many Impulses? // 3rd and 4th Aerospace Sciences Meeting. New York, NY. U.S.A.: American Institute of Aeronautics and Astronautics, 1966.

  5. Gunther P.  Asymptotically Optimum Two-impulse Transfer from Lunar Orbit // AIAA Journal. 1966. V. 4. № 2. P. 346–352.

  6. Jones D., Ocampo C. Optimal Impulsive Escape Trajectories from a Circular Orbit to a Hyperbolic Excess Velocity Vector // AIAA/AAS Astrodynamics Specialist Conf. Toronto, Ontario, Canada, 2010. P. 1–25.

  7. Jones D.R., Ocampo C. Optimization of Impulsive Trajectories from a Circular Orbit to an Excess Velocity Vector // J. Guidance, Control, and Dynamics. 2012. V. 35. № 1. P. 234–244.

  8. Ocampo C., Saudemont R.R. Initial Trajectory Model for a Multi-Maneuver Moon-to-Earth Abort Sequence // J. Guidance, Control, and Dynamics. 2010. V. 33. № 4. P. 1184–1194.

  9. Егоров В.А., Гусев Л.И. Динамика перелетов между Землей и Луной. М.: Наука. Физматлит, 1980.

  10. Гаврикова Н.М., Голубев Ю.Ф. Построение траектории возврата с орбиты ИСЛ к точке входа в атмосферу Земли // Препринты ИПМ им. М.В. Келдыша. 2019. № 53. С. 1–39.

  11. Жаров В.Е. Сферическая астрономия. Фрязино: Век-2, 2006.

  12. Голубев Ю.Ф., Грушевский А.В., Корянов В.В., Тучин А.Г., Тучин Д.А. Основное свойство интеграла Якоби для гравитационных маневров в Солнечной системе // Препринты ИПМ им. М.В. Келдыша. 2019. № 34. С. 1–24.

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

  14. Kraft D. A Software Package for Sequential Quadratic Programming // Technical Report DFVLR-FB 88-28. Oberpfaffenhofen: Institut fuer Dynamik der Flugsysteme. 1988.

Дополнительные материалы отсутствуют.