Космические исследования, 2020, T. 58, № 2, стр. 138-148

ВВЕДЕНИЕ

Исследование малых тел Солнечной системы (астероидов, комет) с помощью автоматических межпланетных станций, которые позволяют изучать эти тела на близком расстоянии, а также контактным методом, стало одним из главных направлений исследования дальнего космоса. Уже осуществлен ряд миссий для изучения малых тел Солнечной системы, например, “Вега”, “Near”, “Deep Space-1”,“Stardust”,“Rosetta”,“Deep Impact”, “Hayabusa” и др. К основным целям этих исследований можно отнести выявление механизма происхождения и эволюции Солнечной системы, поиск ценных ресурсов, изучение проблемы астероидно-кометной опасности, а также демонстрацию мер по предотвращению угрозы Земле.

Сейчас возрастает роль экспедиций к малым небесным телам с возвращением космического аппарата от этих тел к Земле. К данному моменту реально разработаны 4 космической миссии к малым небесным телам с возвращением к Земле. Из них две миссии (“Stardust”, NASA, и “Hayabusa”, Япония) уже совершили возврат спускаемого аппарата (СА) на Землю [1, 2]. И еще 2 миссии (“Hayabusa-2”, Япония, и “OSIRIS-REx” [3], NASA) осуществляются сейчас.

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

Астероид Апофис представляет потенциально серьезную угрозу для безопасности Земли [4, 5]. В России НПО Лавочкина разрабатывает проект полета к этому “опасному” астероиду Апофис [6]. Настоящая работа посвящена разработке алгоритма построения оптимальных орбит КА для экспедиции Земля–астероид–Земля, главными целями которой являются исследование астероида, взятие и доставка образцов его грунта обратно на Землю. В работе сначала рассматривается задача оптимизации траекторий полета к произвольному астероиду с возвратом к Земле, затем результаты анализа будут применены для экспедиции к астероиду Апофис. Данная работа является развитием [7, 8].

Рассмотрена следующая схема экспедиции Земля–астероид–Земля. Ракета-носитель (РН) выводит КА с разгонным блоком (РБ) большой тяги на опорную орбиту искусственного спутника Земли. После пассивного движения по ней в некоторый оптимальный момент t₀ двигательная установка РБ со скоростью истечения c₁ сообщает КА импульс скорости ∆V₁, производится разгон КА, и КА переводится на орбиту полета к астероиду. Затем РБ отделяется от КА, и дальнейшие маневры осуществляются с помощью второй двигательной установки большой тяги, ДУ2, со скоростью истечения газов c₂. В момент t_1сд аппарат выходит из сферы действия Земли. Далее, в момент t₂ КА подлетает к астероиду. С помощью ДУ2 сообщается импульс скорости ∆V₂, осуществляется торможение КА, и КА переходит на круговую орбиту искусственного спутника астероида (ИСА) радиусом r₂. В окрестности астероида КА пребывает некоторое время ∆t₂₃, это – “время ожидания” [9]. В течение этого времени возможны посадка на поверхность астероида, взятие образцов его грунта и другие исследования. Затем, в момент t₃ основному КА сообщается импульс скорости ∆V₃, КА разгоняется и переходит на траекторию возвращения к Земле. В момент t_4сд КА подлетает к сфере действия Земли. При подлете к Земле от КА отделяется спускаемый аппарат, и в момент t_f происходит его гиперболический вход в атмосферу Земли, затем – торможение, посадка. Для уменьшения энергетических затрат в качестве основной принята схема полета, когда тормозной импульс скорости двигателем не сообщается при подлете к Земле, как предлагал К.Э. Циолковский [10] и как сделано в проектах “Stardust”, “Hayabusa”. В этом случае энергетические затраты на экспедицию в номинале определяются тремя величинами импульсов скорости ΔV₁, ΔV₂, ΔV₃.

1. ПОСТРОЕНИЕ ОПТИМАЛЬНЫХ ПО ПОЛЕЗНОЙ МАССЕ ТРАЕКТОРИЙ КА ДЛЯ ЭКСПЕДИЦИИ ЗЕМЛЯ–АСТЕРОИД–ЗЕМЛЯ

Разработана и использовалась двухэтапная методика определения оптимальных межпланетных траекторий полета КА от Земли к астероиду и возврата от астероида к Земле в классе трехимпульсных перелетов. На первом этапе оптимальные гелиоцентрические траектории перелета КА Земля–астероид и астероид–Земля определяются в модели точечных сфер действия Земли и астероида, в импульсной постановке. В основном варианте оптимизации максимизируем полезную массу экспедиции m_p, определяемую как конечная масса КА при возвращении к Земле m_f минус масса ДУ2 вместе с ее топливными баками, масса которых зависит от массы топлива, т.е. импульсов ΔV₁, ΔV₂, ΔV₃. На втором этапе анализа уточняются характеристики полученных оптимальных межпланетных траекторий в более точной модели с учетом параметров траекторий на геоцентрических и астероидоцентрических участках, возмущений, эфемерид небесных тел и гравитационных потерь на участках маневров.

1.1. ПЕРВЫЙ ЭТАП – ПОСТРОЕНИЕ ОПТИМАЛЬНЫХ ТРАЕКТОРИЙ КА В ПРИБЛИЖЕННОЙ МОДЕЛИ

На первом этапе гелиоцентрические траектории перелета КА Земля–астероид и астероид–Земля определяются в модели точечных сфер действия Земли и астероида, поэтому орбиты этих перелетов строятся в центральном ньютоновском поле притяжения Солнца. Схема решения задачи будет следующей. При задании граничных времен экспедиции t₁ (отлет с орбиты Земли), t₂ (подлет к орбите астероида), t₃ (отлет с орбиты астероида), t₄ (подлет к орбите Земли) гелиоцентрические орбиты перелета между небесными телами определяются путем двукратного решения задачи Эйлера–Ламберта (с учетом возможности совершения одного пассивного витка хотя бы по одной орбите). Это позволяет найти скорости “на бесконечности” V_∞1, V_∞2, V_∞3, V_∞4 в граничные времена t_i и требуемые импульсы скорости для перелета ΔV₁, ΔV₂, ΔV₃. По этим скоростям можно определить конечную массу m_f и полезную m_p массу КА [11], с учетом отделяемых масс РБ и ДУ2:

где m₀ – начальная масса КА на опорной орбите ИСЗ, Δm₁ – отделяемая после разгона у Земли масса РБ, m₂₀ – постоянная часть массы ДУ2, a_T2 – коэффициент пропорциональности массы топливных баков ДУ2 массе топлива.

При этом скорости истечения газов с₁, с₂ из двигательных установок РБ и ДУ2, вообще говоря, различны.

Обычно при оптимизации траекторий рассматривают минимизацию характеристической скорости V_хар, равной сумме величин импульсов скорости (см., например, [12, 13]):

или максимизацию конечной массы:

Нами при построении оптимальных межпланетных траекторий перелета в основном варианте анализа максимизируется полезная масса экспедиции m_p:

Этот функционал, вообще говоря, лучше отражает требование энергетической эффективности траектории, чем характеристическая скорость V_хар и конечная масса m_f.

Тогда задача заключается в выборе времен t₁, t₂, t_3,t₄ (при заданных областях для этих времен) для нахождения оптимальных траекторий с максимальной полезной массой.

Классическим способом поиска окна запуска для полета КА к другим планетам является построение графика изолиний функционала с помощью метода прямого перебора (см., например, [14]). При этом время вычисления быстро увеличивается с ростом размерности пространства поиска. Поэтому для уменьшения времени расчета и применения в более сложных задачах разработан комплексный метод, сочетающий несколько методов: метод И.М. Соболя [15–17], генетический алгоритм (ГА) [18] и квазиньютоновский BFGS (Broyden-Fletcher-Goldfarb-Shanno) метод [19].

Сначала для глобальной оптимизации применяется метод И.М. Соболя. В рамках этого метода используется детерминированный подход для построения точек ЛП_τ – последовательностей, которые очень равномерно расположены в области поиска. Используя сетку, построенную точками ЛП_τ – последовательностей, можно получить много различных решений и одновременно оценить максимумы и (или) минимумы функционала. В итоге, точки ЛП_τ – последовательностей позволяют дать хорошее представление о поведении функционала со сравнительно меньшим количеством вычислений, чем в случае прямого перебора. Детальное описание метода построения и свойств точек ЛП_τ – последовательностей дано в работах [15–17]. В частности, метод Соболя позволяет быстро найти области, где расположены локальные оптимумы и глобальный оптимум.

Последующий запуск генетического алгоритма (ГА) в этих областях позволяет определить глобальный оптимум с точностью до суток по граничным временам экспедиции. ГА использует принципы и терминологию, заимствованные у биологической науки – генетики. В ГА каждая особь представляет потенциальное решение некоторой задачи. Каждой особи предоставляется определенная пригодность путем вычисления значения фитнесс-функции. Множество этих особей называется популяцией. В процессе эволюции популяции с помощью генетических операторов – репродукции, кроссинговера и мутации, генерируют новые популяции. При этом особи с лучшими значениями целевой функции (пригодности) c большей вероятностью переходят в следующее поколение, т.е. в следующую итерацию. Таким образом, после нескольких поколений (итераций) решения сходятся к оптимальному. Детальное описание метода ГА дано, например, в работе [18]. Для реализации простого ГА использованы следующие основные операторы: репродукция с помощью рулеточного отбора для выбора “родителей”, одноточечный кроссинговер с вероятностью P_c = 0.8, одноточечная мутация с вероятностью P_m = 0.01 и элитарный отбор “особей” в следующее поколение. Метод Соболя и ГА удобны также при решении задач оптимизации с ограничениями, в частности, при ограничениях на граничные времена экспедиции, на скорость входа в атмосферу Земли при возвращении.

В случае внутреннего, по ограничениям, оптимума при необходимости более точного (в пределах суток) его определения, используется квазиньютоновский BFGS метод [19]. В этом методе, как и в других квазиньютоновских методах, матрица Гессе G_k целевой функции f заменяется на ее аппроксимацию А_k с учетом информации о градиенте функции f(x): g_k = grad f(x_k). При этом производные $\frac{{\partial f}}{{\partial x}}$ определяются численно. Применение BFGS метода позволяет быстро (за 2–8 итераций) уточнить оптимум. Применение этого метода также позволяет снизить требования к точности оптимизации в ГА.

Использование данного комплексного метода позволило быстро и точно проводить оптимизацию траекторий для экспедиции Земля–астероид–Земля. При этом рассмотрено 5 задач:

1) Основная задача оптимизации: при заданной общей продолжительности экспедиции ∆t_∑ = = t₄ – t₁ и заданном времени пребывания КА (времени ожидания) у астероида ∆t₂₃ = t₃ – t₂ оптимизируются время старта t₁ и время перелета от Земли до астероида ∆t₁₂ = t₂ – t₁, чтобы выполнялось (5). Это близко к постановке, данной в работе [13]. В этой задаче, таким образом, есть два параметра оптимизации.

2) При заданном времени ожидания ∆t₂₃ и при ограничении на общую продолжительность экспедиции ∆t_∑ (например, ∆t_∑ ≤ 2 года), оптимизируются времена t₁, ∆t₁₂ и ∆t_∑. Здесь три параметра оптимизации.

3) При заданном суммарном времени экспедиции ∆t_∑ оптимизируются времена t₁, ∆t₁₂, и ∆t₂₃. Здесь также три параметра оптимизации.

4) Оптимизируются все времена ∆t_∑, t₁, ∆t₁₂ и ∆t₂₃. Здесь четыре параметра оптимизации.

5) Полная четырехпараметрическая оптимизация времен ∆t_∑, t₁, ∆t₁₂ и ∆t₂₃ с учетом ограничения на скорость входа КА в атмосферу Земли при возврате от астероида V_вх: V_вх ≤ V_{вх max}.

1.2. ПРОВЕРКА ОПТИМАЛЬНОСТИ ПОЛУЧЕННЫХ ТРАЕКТОРИЙ КА В КЛАССЕ МНОГОИМПУЛЬСНЫХ ПЕРЕЛЕТОВ

После построения на первом этапе оптимальных гелиоцентрических траекторий перелета КА эти траектории проверяем на выполнение необходимых условий оптимальности в классе многоимпульсных перелетов с помощью сопряженных функций. Получены выражения для базис-вектора Лоудена p (вектора λ_V, сопряженного к гелиоцентрической скорости КА V) в граничные времена t₁, t₂, t₃, t₄ для всех трех указанных выше (3–5) функционалов – как для минимума скорости V_хар, так и для максимума конечной массы m_f, и для максимума полезной массы КА m_p [20]. При этом было учтено наличие двух двигательных установок с разными скоростями истечения c₁, c₂ продуктов сгорания, а также наличие отделения массы Δm₁ между сообщением импульсов скорости ΔV₁, ΔV₂. При минимизации характеристической скорости V_хар (3) граничные значения базис-вектора ${{{\mathbf{p}}}_{i}} = {\mathbf{p}}({{t}_{i}})$ определяются выражениями:

Здесь V_p1 – с корость в перигее орбиты отлета от Земли; V_p2, V_p3 – скорости в перицентрах астероидоцентрических орбит подлета к астероиду и отлета от астероида. Радиусы r₂, r₃ соответствуют круговым орбитам спутника астероида после подлета к астероиду и при отлете от него. Отметим, что притяжение астероида часто очень мало. Так для астероида Апофис μ_А ~ 2 м³/с². В то же время скорости на бесконечность V_∞2, V_∞3 велики (километры в секунду – сотни метров в секунду). Поэтому вторые слагаемые под знаком радикала для V_p2, V_p3 (4а) малы по сравнению с первыми слагаемыми. Если мы этими вторыми слагаемыми пренебрегаем, т.е. притяжение астероида не учитывается, то в (6) V_p2 ≈ V_∞2, V_p3 ≈ V_∞3. Для случая максимизации конечной массы (4):

При максимизации полезной массы КА m_p (5):

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

где r, V – гелиоцентрические радиус-вектор и вектор скорости КА, индексы “0” и “f” соответствуют началу и концу каждой дуги перелета,

Эти матрицы задаются, например, формулами для изохронных производных [12, 13, 21–23 и др.]. Поскольку переходная матрица для (p, p') идентична матрице для вариаций (δr, δV), производная базис-вектора p' = dp/dt в начальный момент:

Тогда можно определить базис-вектор на всей траектории:

при этом ${{{\Phi }}_{1}},$ ${{\Phi }_{2}}$ пересчитываются по (9а) в каждый момент t. Можно также определить базис-вектор p(t) интегрированием системы уравнений [13]:

где μ_S – гравитационный параметр Солнца. Для оптимальности траекторий в классе многоимпульсных перелетов необходимо выполнение условия [20]:

Если это условие, играющее здесь роль “принципа максимума”, нарушается на некотором участке, то траекторию можно улучшить введением дополнительных импульсов или вариацией граничных времен [20, 24, 25 ].

Также получены производные от функционалов (3–5) по граничным временам траектории [20], они могут быть использованы для проверки выполнения условий трансверсальности и для улучшения траектории по функционалу, например, градиентным или квазиньютоновским методом. Для функционала F₃ = –m_p (5) имеем [20]:

Если времена t_i (i = 1, 2, 3, 4) не зависят друг от друга и лежат внутри допустимых областей, то для оптимальной траектории должно быть ∂F₃/∂t_i = 0. Если между временами t_i есть связи, то переходим к независимым переменным x_j.

1.3. ВТОРОЙ ЭТАП – УТОЧНЕНИЕ ХАРАКТЕРИСТИК ПОЛУЧЕННЫХ ТРАЕКТОРИЙ КА

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

где r(x, y, z)– радиус-вектор КА в прямоугольной невращающейся системе координат; r_i – радиус-вектор i-го небесного тела (из эфемериды DE421); μ_⊙, μ_Gi – гравитационные параметры центрального и i-го небесных тел соответственно; Δ₁, Δ₂ – ускорения за счет сжатия Земли (J₂) и давления солнечного света. Данная система (14) рассматривается на двух отрезках t ∈ [t₀, t₂], t ∈ [t₃, t_f]. Здесь t₀ – момент отлета с опорной орбиты спутника Земли, t₂ – момент перехода КА на орбиту спутника астероида, t₃ – момент отлета КА с орбиты спутника астероида для полета к Земле, t_f – момент вход в атмосферу Земли при возврате c астероида.

Траектория перелета КА уточняется решением двух задач. Внутренняя задача – решение краевых задач. Для перелета от Земли до астероида варьируем три параметра в начальной точке траектории отлета от Земли, например Ω, ω и |V_∞1|, с которыми одновременно варьируются и t₀, r₀, V₀. Первоначальные приближения t₀, r₀, V₀ определяются с помощью t_1СДЗ (= t₁ из первого этапа) и V_∞1. Добиваемся выполнения условия, что в момент t₂ прицельная дальность в картинной плоскости равняется 1 км. Для перелета от астероида до Земли задача решается в два этапа. Сначала варьируется вектор скорости V_∞3, чтобы обеспечить в момент t_4СДЗ (в качестве начального приближения t_4СДЗ, взято оптимальное время t₄ из первого этапа) вход КА в сферу действия Земли с вектором прицельной дальности, соответствующим входу в атмосферу Земли с необходимой высотой условного перигея. Через эту высоту условного перигея орбиты подлета к Земле также обеспечим угол входа в атмосферу Земли. Для численного анализа полагаем, что эта высота равна 50 км. Затем учитываются возмущения до входа в атмосферу, и получается момент t_f . Первые решения внутренней задачи называем “квазиоптимальными” траекториями. Внешняя задача – оптимизация возмущенных траекторий по граничным временам t₀, t₂, t₃, t_f . В связи с переходом от упрошенной модели к полной модели, оптимальные решения будут слегка отличаться от решений, полученных на первом этапе. Поэтому на втором этапе также выполняется оптимизация – методом покоординатного спуска [30] при малых вариациях. Таким образом, на каждой итерации решаются краевые задачи. После этого уже получаются оптимальные траектории в полной модели.

При коррекции массо-энергетических характеристик экспедиции, прежде всего, был учтен неимпульсный характер разгона КА возле Земли. В этом случае при расчете характеристической скорости отлета КА от Земли $\Delta V_{{\text{1}}}^{*}$ принимаются во внимание “гравитационные потери” δV_1gr. Они оценивались численно и аналитически, на основе формулы [26–29]:

где t_e – время работы двигателя для создания импульса скорости ∆V₁ за счет конечной тяги на расстоянии r от центра с гравитационным параметром μ, k ≈ 0.019 [28]. Эта константа немного уточнялась на основе численных расчетов. Для уменьшения этих потерь рассматривается режим разгона с двумя-тремя включениями двигателя. Для коррекции массо-энергетических характеристик также предусмотрены дополнительные импульсы скорости на коррекцию траекторий и т.д.

2. ИССЛЕДОВАНИЕ ОПТИМАЛЬНЫХ ПО ПОЛЕЗНОЙ МАССЕ ТРАЕКТОРИЙ КА ДЛЯ ЭКСПЕДИЦИИ ЗЕМЛЯ–АПОФИС–ЗЕМЛЯ

С помощью разработанного метода, для экспедиции “Земля–Апофис–Земля” построены и исследованы энергетически оптимальные, по максимуму полезной массы КА, межпланетные траектории при полете в течение 2019–2022 гг. Исследованы варианты с использованием РН “Союз-2.1a”, “Союз-2.1б”, “Зенит” и РБ “Фрегат”. Для двигательной установки ДУ2 удельная тяга 304 с; в (1а) постоянная составляющая массы m₂₀ = 100 кг; коэффициент массы топливных баков а_T2 = 0.15. Начальная опорная орбита пребывания основного КА у Апофиса взята круговой с радиусом 500 м.

При коррекции массовых характеристик гравитационные потери импульсов скорости ΔV₂, ΔV₃ не учитываются, т.к. гравитация Апофиса мала (μ_А ~ 2–3 м³/с²). Для надежности анализа при уточнении траектории предусмотрены дополнительные импульсы скорости на коррекцию: до подлета к астероиду – 50 м/с, возле астероида – 10 м/с и после отлета КА от астероида – 25 м/с. По рекомендации НПО им. С.А. Лавочкина были скорректированы некоторые параметры РБ “Фрегат” и КА, в частности уточнены массовые характеристики отделяемых частей блока; предложены гарантийные запасы топлива; введены отделяемые массы – для небольшого спутника астероида и посадочного устройства – 10 и 20 кг.

Сначала приведем данные в случае использования РН “Союз-2.1a”. Для задачи оптимизации 1) фиксировано время ожидания ∆t₂₃ = t₃ – t₂ = 7 сут. При этом суммарное время ∆t_∑ выбиралось из множества T_S = [390; 420; 450; 510; 540; 570; 600; 630; 660; 690; 730] сут. Решение этой задачи на первом этапе анализа дало оптимальные траектории для разных времен ∆t_∑. Среди них максимальная полезная масса КА m_p = 272 кг получена для траектории № 1, для которой: ∆t_∑ = 690 сут, t₁ = 24.V.2019, Δt₁₂ = t₂ – t₁ = 335 сут, t₄ = 13.IV.2021, V_хар = 6.618 км/с, m_f = 527 кг, см. рис. 1, где представлены зависимость полезной массы от времен t₁ и ∆t₁₂ для варианта Δt_∑ = 690 сут, а также изолинии полезной массы в окрестности глобального оптимума. На рис. 2 приведены картины перелета по траектории № 1: здесь рис. 2а – для полета от Земли (P₁) до Апофиса (P₂); рис. 2б – для полета от Апофиса (P₃) до Земли (P₄). На второй дуге полета делается пассивный виток, прилет к Земле – у восходящего узла орбиты Апофиса.

Решена задача 2) трехмерной оптимизации при заданном времени ожидания Δt₂₃ ∈ [7; 30; 60; 90; 120; 130] сут, при ∆t_∑ ≤ 2 года. Получено, что оптимальное время ожидания КА у Апофиса ∆t_23 opt достигается на интервале 90–120 сут, см. рис. 3. Так, для траектории с Δt₂₃ = 120 сут: t₁ = 6.V.2020, Δt₁₂ = 297 сут, Δt_Σ = 716 сут, V_хар = 6.35 км/с, m_p = = 328 кг.

В задаче (3) трехмерной оптимизации при задании ∆t_∑ = 690 сут, для оптимального решения (траектория № 3): t₁ = 23.V.2019, Δt₁₂ = 336 сут, Δt_23opt = 93 сут, V_хар = 6.519 км/с, m_f = 544 кг, m_p = = 293 кг. Картинки траектории № 3 в плоскости эклиптики показаны на рис. 4.

Решена также задача 4) полной четырехмерной оптимизации при условиях ∆t_∑ ≤ 2 года; Δt₂₃ ≥ 7 сут. Для нее получена траектория № 4 с характеристиками: ∆t_∑ = 716 сут, t₁ = 05.V.2020, Δt₁₂ = 300 сут, Δt_23opt = 112 сут, V_хар = 6.343 км/с, m_f = 545 кг, m_p = = 329 кг. Картинки траектории перелетов № 4 в плоскости эклиптики показаны на рис. 5. Возврат к Земле происходит также вблизи восходящего узла орбиты Апофиса, как и для траекторий № 1 и 3.

Для полученных траекторий построены сопряженные функции, в частности, модуль базис-вектора. На рис. 6 дано изменение модуля базис-вектора p(t) для основной траектории № 1. Здесь p(t) > 1 на некотором начальном отрезке второй дуги, после приложения импульса ΔV₃, условие (12) не выполняется, поэтому траекторию можно улучшить. Поскольку $p_{3}^{'}$ = dp₃/dt > 0, из выражений (13) для частных производных от функционала видно, что если dt₃ > 0, то уменьшим функционал, dF₃ < 0. Это подтверждает улучшение характеристик траектории № 3, у которой большее время t₃. И время ожидания у нее, Δt₂₃ = 93 сут, как раз соответствует интервалу оптимального времени ожидания 90–120 сут. На рис. 7 представлено изменение p(t) для этой траектории. Здесь p ≤ 1, условие оптимальности выполняется. Для траектории № 4 условие p(t) ≤ 1 также выполняется.

Далее, на втором этапе анализа, выполнено уточнение полученных на первом этапе характеристик оптимальных траекторий. Проводится интегрирование уравнений движения КА на участках перелета от Земли к Апофису и от Апофиса к Земле, и решаются соответствующие краевые задачи. Система уравнений движения КА (14) интегрируется методом Рунге-Кутта-Фельберга 7-го порядка с контролем ошибки 8 порядка. Затем в окрестности граничных времен полученных траекторий выполнена оптимизация задачи на множестве уточненных траекторий – методом покоординатного спуска. Также сделана коррекция массово-энергетических характеристик экспедиции. Полученные характеристики оптимальных траекторий № 1а, 3а, 4а приведены в табл. 1, где Δt_s = t_f – t₀. Уточнение на втором этапе привело к некоторому уменьшению полезной массы. При численном анализе получено также, что если о-птимальное время t₀ варьировать в окрестности [t₀ – 3 сут, t₀ + 3 сут], то полезная масса экспедиции меняется мало, в пределах 1 кг. Это, с одной стороны, обеспечивает большое значение полезной массы КА, и с другой стороны, позволяет при необходимости (плохой погоде, отказе в работе бортовой системы и т.д.) переносить дату запуска в небольшом диапазоне.

Кроме РН “Союз-2.1а”, рассмотрено использование РН “Союз-2.1б” и “Зенит” с РБ “Фрегат”. При этом для РН “Союз” при разгоне у Земли РБ делается два включения, для РН “Зенит” – три включения РБ, с целью уменьшения гравитационных потерь. В табл. 2 приведены значения масс m_f и m_p для траекторий № 1a, 3а, 4а при использовании этих РН.

В основном варианте анализа полагали, что при подлете к Земле тормозной импульс скорости двигателем не сообщается, аналогично [1, 2]. КА в момент t_4СД подлетит к сфере действия Земли, и дальше входит в сферу действия Земли со скоростью V_СДЗ по гиперболической орбите с некоторым расстоянием условного перигея r_p4. Затем КА входит в атмосферу Земли. Для оценки скорости входа V_вх в качестве границы атмосферы Земли берем h_a = 125 км [1, 2]. Определена скорость V_вх для полученных оптимальных траекторий. Получается, что V_вх = 12.74, 12.32, 13.26 км/c, для траекторий № 1a, 3а, 4а, соответственно. Для оценки их технологической реализуемости сравниваются их характеристики с теми для миссий “Stardust”, “Hayabusa” и “OSIRIS-REx”, табл. 3. Опыт этих миссий показывает, что современные технологии позволяют успешно приземляться КА, который входит в атмосферу Земли с гиперболической скоростью до 12.9 км/с [1, 2]. Но самая лучшая траектория, полученная нами, № 4а не удовлетворяет этому условию. Поэтому рассмотрена еще одна задача – полная четырехпараметрическая оптимизация с ограничением на скорость входа в атмосферу Земли, V_вх ≤ V_{вх max}. В табл. 4 приведены результаты анализа зависимости полезной массы от V_вх (12.8–13.26 км/с). Отметим, что при этом при использовании РН “Союз-2.1б” и РН “Зенит” полезная масса остается большей, чем сухая масса КА “Stardust” (300 кг).

Таким образом, выполненный анализ показывает, что существует принципиальная возможность осуществить экспедицию к астероиду Апофис в течение 2019–2022 гг. на основе существующих ракет типа “Союз”, “Зенит”.

ЗАКЛЮЧЕНИЕ

Разработана методика построения энергетически оптимальных, с максимальной полезной массой КА, траекторий для экспедиции Земля–астероид–Земля при использовании обычных двигательных установок большой тяги. Разработанная методика применена к построению оптимальных траекторий космической экспедиции Земля–Апофис–Земля, предназначенной для изучения этого опасного астероида. Получены оптимальные траектории данной экспедиции при запуске КА в течение 2019–2022 гг. с общей продолжительностью полета до двух лет. Для анализа оптимальности полученных межпланетных траекторий в классе многоимпульсных перелетов разработан алгоритм построения базис-вектора, т.е. вектора, сопряженного к вектору скорости. Построен базис-вектор для ряда решений. Выявлены качественные характеристики оптимальных траекторий: (1) возврат к Земле происходит вблизи восходящего узла орбиты Апофиса; (2) на номинальной траектории нужно не более трех импульсов скорости; (3) получен оптимальный диапазон времени ожидания у Апофиса – в пределах 3–4 мес. Результаты исследования показывают, что существует принципиальная возможность реализации космической экспедиции Земля–Апофис–Земля с использованием обычных двигателей большой тяги и РН типа “Союз”, “Зенит”.

В заключение авторы выражают искреннюю признательность сотрудникам “НПО им. С.А. Лавочкина” к. т. н. В.Г. Полю и к. т. н. А.В. Симонову, а также к.т.н. И.В. Крылову за поддержку и полезные обсуждения работы, а также участникам семинара им. В.А. Егорова в МГУ им. М.В. Ломоносова за очень полезные обсуждения проблемы и настоящей работы. Авторы признательны аспиранту Гуо Пэн за помощь в совершенствовании вычислительной программы и проведении расчетов. В.В. Ивашкин признателен проф. J. Martinez-Garcia, д-ру M. Bello-Mora, д-ру E. Revilla-Pedreira, проф. P. Sanz-Arangues и проф. T. Elices (GMV, Madrid Polytechnic University) за их поддержку исследований проблемы полета к малым небесным телам.