Орбитальное движение искусственного спутника Земли

Модели сил гравитационного притяжения
Модели сил негравитационной природы
Программное обеспечение расчета орбитального движения спутника
Список литературы

$\newcommand{\pprime}{{\prime\prime}}$ $\newcommand{\mat}[1]{\bar{#1}}$ $\newcommand{\vecg}[1]{\vec{#1}}$ $\renewcommand{\epsilon}{\varepsilon}$ $\renewcommand{\phi}{\varphi}$ $\renewcommand{\leq}{\leqslant}$ $\renewcommand{\geq}{\geqslant}$

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

Предполагается, что векторы $\vec{r}$ и $\vec{v}$, определяющие положение и скорость спутника, задаются в инерциальной системе координат. Уравнения движения спутника имеют вид

\begin{equation} \frac{d^2 \vec{r}}{d t^2} = \vec{f}, \label{basic} \end{equation}

где $\vec{f}$ — ускорение спутника, равное

$$ \vec{f} = \vec{f}_E + \vec{f}_S + \vec{f}_M + \vec{f}_a + \vec{f}_p , $$

$\vec{f}_E$ — ускорение, обусловленное действием поля тяготения Земли;

$\vec{f}_S$ — ускорение, вызванное притяжением Солнца;

$\vec{f}_M$ — ускорение, вызванное притяжением Луны;

$\vec{f}_a$ — ускорение, обусловленное сопротивлением движению в атмосфере;

$\vec{f}_p$ — ускорение, вызванное давлением солнечного излучения.

Численное интегрирование уравнений (\ref{basic}) является основой расчета орбитального движения спутника. Для этого необходимо знать выражения сил, действующих на спутник в околоземном пространстве.

Модели сил гравитационного притяжения

Поле тяготения Земли

Потенциал поля тяготения Земли удобно определять в геоцентрической экваториальной системе координат, вращающейся вместе с Землей. Связь между декартовыми $(x, y, z)$ и полярными координатами $(r, \varphi, \lambda)$ задается формулами (рис.)

$$ \begin{aligned} x &= r \cos\varphi \cos\lambda, \\ y &= r \cos\varphi \sin\lambda, \\ z &= r \sin\varphi, \end{aligned} $$

где $\varphi$ — широта;

$\lambda$ — долгота;

$r$ — расстояние от начала координат до точки наблюдения.

$Связь между декартовыми $(x, y, z)$ и полярными координатами $(r, \varphi, \lambda)$ в геоцентрической экваториальной системе координат$

Связь между декартовыми $(x, y, z)$ и полярными координатами $(r, \varphi, \lambda)$ в геоцентрической экваториальной системе координат

Согласно [1], потенциал поля тяготения Земли записывается в виде разложения в ряд по сферическим функциям

\begin{equation} U (r, \varphi, \lambda) = \frac{\mu}{r} \sum_{n=0}^{\infty} \sum_{m=0}^n \left(\frac{R}{r}\right)^n P_{nm}(\sin\varphi) \left( C_{nm}\cos m\lambda + S_{nm} \sin m\lambda \right) , \label{U} \end{equation}

где $\mu = 398600,4415$ км$^3$/c$^2$ — гравитационная постоянная Земли;

$R = 6378,137$ км — радиус Земли;

$P_{nm}$ — присоединенная функция Лежандра степени $n$ и порядка $m$;

$C_{nm}$ и $S_{nm}$ — коэффициенты, зависящие от распределения масс внутри Земли и задаваемые таблично.

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

$$ \begin{aligned} V_{nm} &= \left( \frac{R}{r} \right)^{n+1} P_{nm} (\sin\varphi) \cos m\lambda, \\ W_{nm} &= \left( \frac{R}{r} \right)^{n+1} P_{nm} (\sin\varphi) \sin m\lambda. \end{aligned} $$

Тогда выражение для потенциала поля тяготения Земли (\ref{U}) приобретет вид

$$ U = \frac{\mu}{R} \sum_{n=0}^{\infty} \sum_{m=0}^n \left( C_{nm} V_{nm} + S_{nm} W_{nm} \right) . $$

Для вычисления присоединенных функций Лежандра степени $n$ и порядка $m$ используются рекуррентные соотношения. Так, все функции вида $P_{mm}$, начиная с $P_{00}=1$, вычисляются при помощи формулы

\begin{equation} P_{mm}(u) = (2m-1)(1-u^2)^{1/2}P_{m-1,m-1}(u), \label{Pmm} \end{equation}

где $u = \sin\varphi$, $(1-u^2)^{1/2} = \cos\varphi$.

Оставшиеся функции вычисляются из соотношений

$$ P_{m+1,m}(u) = (2m+1)uP_{mm}(u), $$

\begin{equation} P_{nm}(u) = \frac{1}{n-m} \left( (2m-1)uP_{n-1,m}(u) - (n+m-1)P_{n-2,m}(u) \right), \label{Pnm} \end{equation}

где $n > m+1$.

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

$$ \begin{aligned} \cos (m+1)\lambda &= \cos m\lambda \cos\lambda - \sin m\lambda \sin\lambda, \\ \sin (m+1)\lambda &= \sin m\lambda \cos\lambda + \cos m\lambda \sin\lambda, \end{aligned} $$

и соотношениями (\ref{Pmm})—(\ref{Pnm}), запишем рекуррентные формулы для вычисления $V_{nm}$ и $W_{nm}$

\begin{equation} \begin{aligned} V_{mm} &= (2m-1) \left[ \frac{xR}{r^2} V_{m-1,m-1} - \frac{yR}{r^2} W_{m-1,m-1}\right], \\ W_{mm} &= (2m-1) \left[ \frac{xR}{r^2} W_{m-1,m-1} + \frac{yR}{r^2} V_{m-1,m-1}\right] \end{aligned} \label{VWmm} \end{equation}

\begin{equation} \begin{aligned} V_{nm} &= \left(\frac{2n-1}{n-m}\right) \cdot \frac{zR}{r^2} V_{n-1,m} - \left(\frac{n+m-1}{n-m}\right) \cdot \frac{R^2}{r^2} V_{n-2,m}, \\ W_{nm} &= \left(\frac{2n-1}{n-m}\right) \cdot \frac{zR}{r^2} W_{n-1,m} - \left(\frac{n+m-1}{n-m}\right) \cdot \frac{R^2}{r^2} W_{n-2,m}. \end{aligned} \label{VWnm} \end{equation}

Отметим, что для вычисления $V_{m+1,m}$ и $W_{m+1,m}$ удобно использовать частный случай формулы (\ref{VWnm})

\begin{equation} \begin{aligned} V_{m+1,m} &= (2m+1) \cdot \frac{zR}{r^2} V_{mm}, \\ W_{m+1,m} &= (2m+1) \cdot \frac{zR}{r^2} W_{mm} . \end{aligned} \label{VWm_plus_1m} \end{equation}

Вычисления $V_{nm}$ и $W_{nm}$ производятся в следующем порядке.

Значения $V_{00}$ и $W_{00}$ известны и равны

$$ V_{00} = \frac{R}{r}, \quad W_{00} = 0. $$

Пользуясь ими, вычислим:

$V_{10}$, $W_{10}$ по формулам (\ref{VWm_plus_1m}) при $m=0$;
$V_{n0}$ для $n > 1$ по формулам (\ref{VWnm}) (при этом все $W_{n0}=0$);
$V_{11}$, $W_{11}$ по формулам (\ref{VWmm}).

Затем цикл вычислений повторяется:

вычислим $V_{21}$, $W_{21}$ по формулам (\ref{VWm_plus_1m});
вычислим $V_{n1}$, $W_{n1}$ для $n > 2$ по формулам (\ref{VWnm});
вычислим $V_{22}$, $W_{22}$ по формулам (\ref{VWmm})

и т. д.

Таким образом, вычисления (рис.) выполняются слева направо и сверху вниз. Стрелка $\downarrow$ означает расчет по формулам (\ref{VWnm}) (или (\ref{VWm_plus_1m}), как частному случаю (\ref{VWnm})), а $\searrow$ — по формулам (\ref{VWmm}).

$Схема вычисления коэффициентов $V_{nm}$, $W_{nm}$$

Схема вычисления коэффициентов $V_{nm}$, $W_{nm}$. Стрелка $\downarrow$ означает расчет по формулам (\ref{VWnm}) (или по их частному случаю (\ref{VWm_plus_1m})), а $\searrow$ — расчет по формулам (\ref{VWmm})

Ускорение $\ddot{\vec{r}}$, равное градиенту потенциала $U$, в декартовых координатах запишется как

$$ \ddot x = \sum_{n,m} \ddot x_{nm}, \quad \ddot y = \sum_{n,m} \ddot y_{nm}, \quad \ddot z = \sum_{n,m} \ddot z_{nm}, $$

где $\ddot x_{nm}$, $\ddot y_{nm}$ и $\ddot z_{nm}$ определяются по формулам

\begin{align} \ddot x_{nm} & \stackrel{(m=0)}{=} \frac{\mu}{R^2} \cdot \left[ -C_{n0}V_{n+1,1} \right] \notag\\ & \stackrel{(m>0)}{=} \frac{\mu}{R^2} \cdot \frac{1}{2} \cdot \left[ (-C_{nm}V_{n+1,m+1} - S_{nm}W_{n+1,m+1}) +\right. \notag\\ &\phantom{{}\stackrel{(m=0)}{=}} \left. + \frac{(n-m+2)!}{(n-m)!} (C_{nm}V_{n+1,m-1} + S_{nm}W_{n+1,m-1}) \right], \notag\\ \ddot y_{nm} & \stackrel{(m=0)}{=} \frac{\mu}{R^2} \cdot \left[ -C_{n0}W_{n+1,1} \right] \\ & \stackrel{(m>0)}{=} \frac{\mu}{R^2} \cdot \frac{1}{2} \cdot \left[ (-C_{nm}W_{n+1,m+1} + S_{nm}V_{n+1,m+1}) +\right. \notag\\ &\phantom{{}\stackrel{(m=0)}{=}} \left. + \frac{(n-m+2)!}{(n-m)!} (-C_{nm}W_{n+1,m-1} + S_{nm}V_{n+1,m-1}) \right], \notag\\ \ddot z_{nm} & = \frac{\mu}{R^2} \cdot \left[ (n-m+1) \cdot (-C_{nm}V_{n+1,m} - S_{nm}W_{n+1,m}) \right]. \notag \end{align}

Приливные деформации Земли

Согласно [2], возмущения от приливных деформаций центрального тела вносят заметный вклад в возмущения искусственных спутников Земли.

Гравитационное воздействие сторонних тел, таких как Луна и Солнце, приводит к деформации Земли, вследствие чего изменяется её гравитационное поля. Простейшей моделью представления потенциала действующих на спутник сил, обусловленных приливными деформациями Земли, является модель Лява (Love) или модель твердых приливов [2].

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

Более точные модели учитывают влияние приливных деформаций, возникающих в океане и в атмосфере Земли. В этом случае влияние приливов может быть выражено в виде поправок к коэффициентам геопотенциала $C_{nm}$, $S_{nm}$, учитывающих неупругие деформации и периодические вариации чисел Лява (модель Вара (Wahr) [2]).

Притяжение Солнца и Луны

В инерциальной системе координат ускорение спутника, вызываемое притяжением точечной массы $M$ (материальной точки), равно

\begin{equation} \ddot{\vec{r}} = GM\frac{\vec{s}-\vec{r}}{|\vec{s}-\vec{r}|^3}, \label{pointmass1} \end{equation}

где $\vec{r}$ и $\vec{s}$ — геоцентрические координаты спутника и $M$ соответственно.

Следует учесть, что под действием притяжения $M$ центр масс Земли, в свою очередь, приобретает ускорение

\begin{equation} \ddot{\vec{r}} = GM\frac{\vec{s}}{|\vec{s}|^3}. \label{pointmass2} \end{equation}

В результате, ускорение спутника в относительной системе координат связанной с центром масс Земли равно разности ускорений (\ref{pointmass1}) и (\ref{pointmass2})

\begin{equation} \ddot{\vec{r}} = GM \left( \frac{\vec{s}-\vec{r}}{|\vec{s}-\vec{r}|^3} - \frac{\vec{s}}{|\vec{s}|^3} \right). \label{pert:pointmass} \end{equation}

Возмущающие влияния Солнца и Луны принято считать независимыми друг от друга, поэтому формулой (\ref{pert:pointmass}) можно пользоваться в обоих случаях.

Для вычисления силы притяжения по формуле (\ref{pert:pointmass}) необходимо знать координаты возмущающего тела. Во многих случаях для этого достаточно использовать простые уравнения, которые следуют из аналитических теорий движения Солнца (Ньюкома) и Луны (Брауна) и обеспечивают отклонение, не превышающее $1%$ от наиболее совершенных современных теорий [1].

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

$$ \begin{aligned} a &= 149600000 \text{км},\\ e &= 0,016709,\\ i &= 0^\circ,0000,\\ \Omega+\omega &= 292^\circ,9400,\\ M &= 357^\circ,5256 + 35999^\circ,049 \cdot T, \end{aligned} $$

где $T = (JD - 2451545,0)/36525,0$ — число юлианских столетий, прошедших с эпохи J2000.

В астрономии эпоха — это момент времени, указанный в какой-либо шкале времени. Эпоха J2000 — полдень ($12^h$) 1 января 2000 г. по земному времени TT (Terrestrial Time) — стандартная эпоха современных астрономических расчетов.

С учетом малости эксцентриситета и наклонения, координаты Солнца представляются в форме рядов [1]

\begin{equation} \begin{aligned} \lambda_\odot &= \Omega + \omega + M + 6892^\pprime \sin M + 72^\pprime \sin 2M, \\ r_\odot &= (149,619 - 2,499 \cos M - 0,021 \cos 2M) \cdot 10^6 \text{км}, \end{aligned} \label{sun_coords} \end{equation}

где $\lambda_\odot$ и $r_\odot$ — геоцентрические эклиптические долгота и расстояние до Солнца соответственно. Отметим, что геоцентрическая эклиптическая широта Солнца $\beta_\odot$ равна нулю с точностью не хуже $1^\prime$ [1].

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

\begin{equation} \mat{r}_\odot = R_x(-\epsilon)\left( \begin{array}{c} r_\odot \cos\lambda_\odot \cos\beta_\odot \\ r_\odot \sin\lambda_\odot \cos\beta_\odot \\ r_\odot \sin\beta_\odot \end{array} \right), \label{cartesian_sun} \end{equation}

где $\epsilon = 23^\circ,43929111$ — наклон эклиптики относительно земного экватора;

$$ R_x(\epsilon) = \left( \begin{array}{ccc} 1 & 0 & 0 \\ 0 & +\cos\epsilon & +\sin\epsilon \\ 0 & -\sin\epsilon & +\cos\epsilon \\ \end{array} \right). $$

С учетом равенства $\beta_\odot = 0$, получим из (\ref{cartesian_sun})

\begin{equation} \mat{r}_\odot = \left( \begin{array}{c} r_\odot \cos\lambda_\odot \\ r_\odot \sin\lambda_\odot \cos\epsilon \\ r_\odot \sin\lambda_\odot \sin\epsilon \end{array} \right). \label{cartesian_sun2} \end{equation}

Координаты $\lambda_\odot$, $\beta_\odot$ и $r_\odot$ в (\ref{cartesian_sun}) относятся к среднему равноденствию и эклиптике в эпоху J2000.

Прецессия, являющаяся результатом действия возмущающих сил Солнца, Луны и планет, приводит к медленному движению как эклиптики так и равноденствия. Плоскость эклиптики меняет свою ориентацию сравнительно медленно — менее чем на $1^\prime$ в столетие. Движение равноденствия более выраженно и составляет $5030^\pprime$ в столетие. Поэтому, для того, чтобы найти координаты равноденствия в некоторую эпоху $T_\text{eqx}$ (измеряется в столетиях с начала эпохи J2000), нужно добавить к $\lambda_\odot$ следующую поправку

\begin{equation} 1^\circ,3972 \cdot T_\text{eqx} . \label{equinox_correction} \end{equation}

Эклиптическая широта $\beta_\odot$ в поправках не нуждается, так как изменяется не более чем на $1^\prime$ за столетие.

Таким образом, формулы (\ref{sun_coords}) и (\ref{cartesian_sun2}) с учетом (\ref{equinox_correction}) позволяют вычислить прямоугольные координаты Солнца в геоцентрической экваториальной системе координат.

Разложения в ряды, аналогичные (\ref{sun_coords}), существуют и для координат Луны. Однако, в связи с тем, что Луна испытывает сильные возмущения со стороны Солнца и Земли, для описания ее движения необходимо использовать большее число членов. Рассматриваемые далее ряды позволяют вычислить широту и долготу Луны с точностью до нескольких угловых минут, а расстояние до Луны — с точностью около 500 км. Расчет возмущений основывается на пяти фундаментальных аргументах: средней долготе Луны $L_0$, средней аномалии Луны $l$, средней аномалии Солнца $l^\prime$, среднем угловом расстоянии Луны от восходящего узла $F$ и разности $D$ между средними долготами Солнца и Луны. Долгота восходящего узла $\Omega$ в явном виде не используется, ее получают из разности $\Omega = L_0 - F$.

Фундаментальные аргументы вычисляются по формулам

\begin{equation} \begin{aligned} L_0 &= 218^\circ,31617 + 481267^\circ,88088 \cdot T - 1^\circ,3972 \cdot T, \\ l &= 134^\circ,96292 + 477198^\circ,86753 \cdot T, \\ l^\prime &= 357^\circ,52543 + 35999^\circ,04944 \cdot T, \\ F &= 93^\circ,27283 + 483202^\circ,01873 \cdot T, \\ D &= 297^\circ,85027 + 445267^\circ,11135 \cdot T. \end{aligned} \label{moon_fund_args} \end{equation}

Координаты Луны $\lambda_M$, $\beta_M$ и $r_M$, относящиеся к среднему равноденствию и эклиптике в эпоху J2000, записываются следующим образом [1]

\begin{align} \lambda_M &= L_0 + 22640^\pprime \sin l + 769^\pprime \sin 2l -\notag\\ &\phantom{{}=L_0}-4586^\pprime \sin(l-2D) + 2370^\pprime \sin 2D -\notag\\ &\phantom{{}=L_0}-668^\pprime \sin l^\prime - 412^\pprime \sin 2F -\notag\\ &\phantom{{}=L_0}-212^\pprime \sin(2l-2D) - 206^\pprime \sin(l+l^\prime-2D) +\notag\\ &\phantom{{}=L_0}+192^\pprime \sin(l+2D) - 165^\pprime \sin(l^\prime-2D) +\notag\\ &\phantom{{}=L_0}+148^\pprime \sin(l-l^\prime) - 125^\pprime \sin D -\notag\\ &\phantom{{}=L_0}-110^\pprime \sin(l+l^\prime) - 55^\pprime \sin(2F-2D) ,\notag\\ \beta_M &= 18520^\pprime \sin(F+\lambda_M-L_0+412^\pprime \sin 2F + 541^\pprime \sin l^\prime) - \label{eq:moon_coords}\\ &\phantom{{}=L_0}-526^\pprime \sin(F-2D) + 44^\pprime \sin(l+F-2D) -\notag\\ &\phantom{{}=L_0}-31^\pprime \sin(-l+F-2D) - 25^\pprime \sin(-2l+F) -\notag\\ &\phantom{{}=L_0}-23^\pprime \sin(l^\prime+F-2D) + 21^\pprime \sin(-l+F) +\notag\\ &\phantom{{}=L_0}+11^\pprime \sin(-l^\prime+F-2D) ,\notag\\ r_M &=(385000 - 20905\cos l - 3699 \cos(2D-l) -\notag\\ &\phantom{{}=L_0}-2956 \cos 2D - 570 \cos 2l + 246 \cos(2l-2D) -\notag\\ &\phantom{{}=L_0}-205 \cos(l^\prime-2D) - 171 \cos(l+2D) .\notag \end{align}

Сферические эклиптические координаты (\ref{eq:moon_coords}) преобразуются в прямоугольные экваториальные координаты по формулам

\begin{equation} \mat{r}_M = R_x(-\epsilon)\left( \begin{array}{c} r_M \cos\lambda_M \cos\beta_M \\ r_M \sin\lambda_M \cos\beta_M \\ r_M \sin\beta_M \end{array} \right) . \label{cartesian_moon} \end{equation}

Преобразование координат при переходе от инерциальной системы эпохи J2000 к экватору и равноденствию текущей эпохи $T_\text{eqx}$ выполняется так же, как и для солнечных координат.

Итак, прямоугольные координаты Луны в геоцентрической экваториальной системе координат заданной эпохи вычисляются по формулам (\ref{moon_fund_args}), (\ref{eq:moon_coords}) и (\ref{cartesian_moon}) с учетом поправки (\ref{equinox_correction}).

Более точное определение положений Солнца и Луны

Для определения координат Солнца и Луны широко используются разложения в ряды, основанные на аналитических теориях Ньюкома и Брауна. Точность этих теорий, как правило, вполне приемлема при описании влияния гравитационного притяжения Солнца и Луны на движение низкоорбитальных спутников Земли. Однако при прогнозировании движения спутников на высоких околоземных орбитах, использование теорий Ньюкома и Брауна вносит значительные ошибки в вычисленные положения спутника [2]. Эти ошибки можно существенно уменьшить, если использовать для вычисления координат Луны и Солнца высокоточные эфемериды — таблицы вычисленных положений небесных объектов.

В настоящее время широко употребляются высокоточные эфемериды больших планет, Солнца и Луны DE430 и DE431, построенные в Лаборатории реактивного движения (JPL) NASA. Эти эфемериды получены численными методами и дают координаты в виде полиномов Чебышева в прямоугольной барицентрической системе координат с экватором Земли и равноденствием, относящимися к эпохе J2000. В настоящее время они являются основой «Astronomical Almanac». Эфемериды DE430 охватывают период от 1550 г. до 2650 г., DE431 — от 13201 г. до н. э. до 17191 г. н. э.

Эфемериды, построенные в JPL, можно найти на ftp-сервере лаборатории. Кроме того, на веб-сайте JPL доступна программа HORIZONS для онлайн-генерации эфемерид объектов Солнечной системы (в том числе: астероидов, комет, спутников планет и некоторых космических аппаратов).

Европейский стандарт моделей описания космической среды (ECSS-E-ST-10-04C, от 15.11.2008) предполагает использование более ранней версии эфемерид JPL: DE-405/LE-405 (LE означает: эфемериды Луны — "Lunar Ephemerides").

Эффекты общей теории относительности

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

Модели гравитационного поля Земли

В настоящее время Международный центр глобальных моделей гравитационного поля (International Center for Global Gravity Field Models) насчитывает 157 моделей гравитационного поля Земли. Наиболее современные из них можно получить на веб-сайте этой организации.

Одной из наиболее распространенных моделей является Earth Gravity Model 2008 (EGM2008), содержащая сферические гармоники до 2159 степени и порядка включительно (существует расширенная версия — вплоть до гармоник 2190 степени и 2159 порядка). Она разработана Национальным управлением геопространственной разведки США и используется Международной службой вращения Земли (IERS Conventions, 2010).

Европейский стандарт моделей космической среды ECSS-E-ST-10-04C предполагает использование модели EIGEN-GL04C, разработчиками которой являются исследовательские центры из Потсдама (GFZ Potsdam) и Тулузы (GRGS Toulouse).

Модели сил негравитационной природы

Сопротивление атмосферы

Среди сил негравитационной природы наибольшее воздействие на движение низкоорбитальных спутников (т. е. спутников, движущихся на высотах от 150 до 1500 км) оказывают аэродинамические силы, вызванные влиянием атмосферы Земли. Действие этих сил преимущественно выражается в сопротивлении движению спутника, направленном противоположно его скорости $\vec{v}_r$ относительно атмосферы.

Рассмотрим малый элемент атмосферы массой $\Delta m$, сталкивающийся со спутником с площадью поперечного сечения $A$ за промежуток времени $\Delta t$

$$ \Delta m = \rho A v_r \Delta t , $$

где $\rho$ — плотность атмосферы в месте расположения спутника;

$v_r = \|\vec{v}_r\|$.

Импульс $\Delta \vec{p}$, испытываемый спутником, при этом равен

$$ \Delta \vec{p} = \Delta m \vec{v}_r = \rho A v^2_r \Delta t \cdot \frac{\vec{v}_r}{v_r}, $$

а соответствующая сила $\vec{F}$ определяется из соотношения $\vec{F} = \Delta \vec{p} / \Delta t$. В результате, ускорение спутника, обусловленное сопротивлением атмосферы равно

\begin{equation} \ddot{\vec{r}} = -\frac{1}{2} C_D \frac{A}{m} \rho v^2_r \vec{e}_v, \label{acc_atm_drag} \end{equation}

где $C_D$ — коэффициент сопротивления;

$m$ — масса спутника;

$\vec{e}_v = \vec{v}_r / v_r$.

Коэффициент сопротивления $C_D$ представляет собой безразмерную величину, которая описывает взаимодействие атмосферы с материалом поверхности спутника. Обычно, значения $C_D$ лежат в диапазоне от 1,5 до 3. Направление ускорения противоположно скорости движения спутника относительно атмосферы $\vec{e}_v$. Последняя зависит от сложной динамики атмосферы. Тем не менее, во многих практических случаях разумным приближением будет считать атмосферу вращающейся со скоростью вращения Земли. Следовательно, для $\vec{v}_r$ можно записать

\begin{equation} \vec{v}_r = \vec{v} - \vecg{\omega}_\oplus \times \vec{r}, \label{v_rel} \end{equation}

где $\vec{r}$ и $\vec{v}$ — положение и скорость спутника в инерциальной СК;

$\vecg{\omega}_\oplus$ — вектор угловой скорости вращения Земли ($\omega_\oplus = 0,7292 \cdot 10^{-4}$ рад/c).

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

Экспоненциальная атмосфера

Эта простейшая модель предполагает, что плотность атмосферы экспоненциально убывает с увеличением высоты (см., в частности, ГОСТ 4401-81)

$$ \rho = \rho_0 e^{-h/H_0} , $$

где $\rho_0$ — плотность атмосферы на поверхности земного референц-эллипсоида, то есть приближения формы поверхности Земли эллипсоидом вращения ($\rho_0 = 1,225$ кг/м$^3$ при температуре $15^{\,\circ}$C);

$h$ — высота над поверхностью референц-эллипсоида (уровнем моря);

$H_0$ — масштаб высоты по давлению, определяемый формулой

$$ H_0 = \frac{RT}{Mg} \approx 8,42\ \text{м}, $$

где $R$ — универсальная газовая постоянная;

$M$ — молярная масса сухого воздуха;

$T$ — температура воздуха на уровне моря;

$g$ — ускорение свободного падения.

Модель Харриса-Прейстера (Harris-Priester)

Эта модель, разработанная в начале 1960-х гг. [1], расширяет экспоненциальную модель, позволяя учесть суточные изменения плотности. Нагрев атмосферы солнечным излучением приводит к постепенному увеличению ее плотности — с освещенной стороны Земли возникает выпуклость, вершина которой (апекс) отстает приблизительно на 2 часа (или $30^\circ$ к востоку) от подсолнечной точки. Плотность атмосферы на заданной высоте в апексе $\rho_m(h)$ и противоположной ему относительно Земли точке (антапексе) $\rho_M(h)$ вычисляется с помощью экспоненциальной интерполяции между табулированными минимальным $\rho_m(h_i)$ и максимальным $\rho_M(h_i)$ значениями плотности на высотах $h_i$

$$ \begin{aligned} \rho_m(h) &= \rho_m(h_i)\exp\left(\frac{h_i-h}{H_m}\right) \quad (h_i \leq h \leq h_{i+1}), \\ \rho_M(h) &= \rho_M(h_i)\exp\left(\frac{h_i-h}{H_M}\right), \end{aligned} $$

где $h$ — высота над поверхностью земного референц-эллипсоида.

Соответствующие масштабы высоты по давлению равны

$$ \begin{aligned} H_m(h) &= \frac{h_i - h_{i+1}}{\ln(\rho_m(h_{i+1})/\rho_m(h_i))}, \\ H_M(h) &= \frac{h_i - h_{i+1}}{\ln(\rho_M(h_{i+1})/\rho_M(h_i))}. \end{aligned} $$

Плотность атмосферы в заданной точке $\rho(h)$ выражается через $\rho_m(h)$ и $\rho_M(h)$ следующим образом

$$ \rho(h) = \rho_m(h) + \left( \rho_M(h) - \rho_m(h) \right) \cdot \cos^n \left(\frac{\Psi}{2}\right), $$

где $\Psi$ — угол между радиусом-вектором спутника $\vec{r}$ и вершиной атмосферной выпуклости (апексом).

Широтные вариации атмосферы учитываются приблизительно, выбором показателя степени $n$, который для орбит с низким наклонением равен 2, а для приполярных орбит — 6.

$$ \cos^n \left(\frac{\Psi}{2}\right) = \left( \frac{1+\cos\Psi}{2} \right)^\frac{n}{2} = \left( \frac{1}{2}+\frac{\vec{e}_r \cdot \vec{e}_a}{2} \right)^\frac{n}{2}, $$

где $\vec{e}_r = \vec{r}/r$.

Единичный вектор апекса $\vec{e}_a$ задается как

$$ \vec{e}_a = \left( \begin{array}{l} \cos\delta_\odot \cos(\alpha_\odot + \lambda_l) \\ \cos\delta_\odot \sin(\alpha_\odot + \lambda_l) \\ \sin\delta_\odot \end{array} \right), $$

где $\alpha_\odot$ и $\delta_\odot$ — прямое восхождение и склонение Солнца соответственно;

$\lambda_l \approx 30^\circ$ — угол запаздывания апекса относительно подсолнечной точки по долготе.

Значения $\rho_m(h_i)$ и $\rho_M(h_i)$ задаются таблично (см. [1]).

Современные модели атмосферы

Современными моделями плотности верхней атмосферы является модель ГОСТ Р 25645.166-2004 и модель NASA NRLMSISE-00. Обе они построены по результатам наблюдений со спутников и геофизических ракет и рассчитаны на применение в диапазоне высот от 120 до 1500 км над поверхностью Земли. Сравнительный анализ этих моделей (см. [2]) показал, что по своим возможностям они весьма близки.

Европейский стандарт моделей описания космической среды ECSS-E-ST-10-04C предполагает использование для плотности верхней атмосферы моделей NASA NRLMSISE-00 и Jacchia-Bowman 2006 (JB-2006). Последняя подробно описана в [3].

Солнечное давление

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

Солнечное давление определяется величиной падающего на спутник потока солнечного излучения

$$ \Phi = \frac{\Delta E}{A \Delta t}, $$

то есть энергией $\Delta E$, проходящей через поверхность площадью $A$ за промежуток времени $\Delta t$.

Единичный фотон с энергией $E_\nu$ переносит импульс

$$ p_\nu = \frac{E_\nu}{c}, $$

где $c$ — скорость света.

Следовательно, полный импульс, переданный освещенному телу за время $\Delta t$, равен

$$ \Delta p = \frac{\Delta E}{c} = \frac{\Phi}{c} A \Delta t. $$

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

$$ F = \frac{\Delta p}{\Delta t} = \frac{\Phi}{c} A. $$

Соответственно, давление равно

\begin{equation} P = \frac{\Phi}{c}. \label{sun_pressure} \end{equation}

Поток солнечной энергии на расстоянии 1 а. е. от Солнца (астрономическая единица приблизительно равна среднему расстоянию от Земли до Солнца: 1 а. е. = 149597870,691 км) равен

$$ \Phi \approx 1367\ \text{Вт}\cdot\text{м}^{-2}. $$

Следовательно, согласно (\ref{sun_pressure}), давление солнечного излучения на орбите Земли равно

\begin{equation} P_\odot \approx 4,56 \cdot 10^{-6} \ \text{Н}\cdot\text{м}^{-2}. \label{P_sun} \end{equation}

При этом предполагается, что спутник расположен перпендикулярно потоку излучения и поглощает все падающие на него лучи.

На практике, часть излучения отражается от поверхности спутника. Доля энергии отраженного излучения характеризуется коэффициентом отражения $\epsilon$ ($\epsilon=1$ в случае полного отражения и $\epsilon=0$ — для полного поглощения).

Взаимодействие потока излучения с произвольно ориентированной плоской поверхностью для случаев полного поглощения и отражения показано на рис., где $\vec{n}$ — вектор нормали, определяющий ориентацию поверхности площадью $A$, $\vec{e}_\odot$ — единичный вектор, указывающий направление на Солнце, $\theta$ — угол между $\vec{n}$ и $\vec{e}_\odot$.

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

Сила, действующая на поверхность, полностью поглощающую падающее излучение, равна

\begin{equation} \vec{F}_{a} = -P_\odot \cos\theta A \vec{e}_\odot, \label{Fabs} \end{equation}

где $\cos\theta A$ — площадь поверхности, перпендикулярная потоку излучения.

В свою очередь сила, действующая на полностью отражающую поверхность, равна

\begin{equation} \vec{F}_{r} = -2P_\odot \cos\theta A \cos\theta \vec{n}. \label{Frefl} \end{equation}

Как видно из рис., импульс $\Delta \vec{p}_n$, передаваемый поверхности в направлении нормали, в случае полного отражения в два раза больше, чем в случае полного поглощения.

$Изменение импульса в направлении нормали к поверхности $\Delta \vec{p}_n$: a) полное поглощение, б) полное отражение$

Изменение импульса в направлении нормали к поверхности $\Delta \vec{p}_n$: a) полное поглощение, б) полное отражение

Для тела, отражающего часть энергии $\epsilon \Delta E$ и поглощающего оставшуюся часть $(1-\epsilon) \Delta E$, выражение для силы солнечного давления представляет собой комбинацию (\ref{Fabs}) и (\ref{Frefl})

$$ \vec{F} = -P_\odot \cos\theta A \left[ (1-\epsilon)\vec{e}_\odot + 2\epsilon\cos\theta \vec{n} \right] . $$

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

\begin{equation} \ddot{\vec{r}} = -\nu P_\odot \frac{1\,\text{а.е.}^2}{r_\odot^2} \frac{A}{m} \cos\theta \left[ (1-\epsilon)\vec{e}_\odot + 2\epsilon\cos\theta \vec{n} \right] , \label{sun_accel} \end{equation}

где $m$ — масса спутника;

$\cos\theta = \vec{n} \cdot \vec{e}_\odot$;

$\nu$ — функция тени, которую мы определим ниже.

Эллиптичность орбиты Земли приводит к тому, что расстояние между спутником, находящимся на околоземной орбите, и Солнцем изменяется в течение года в диапазоне от $147\cdot 10^6$ км до $152\cdot 10^6$ км. Это приводит к годичным вариациям солнечного давления в пределах $\pm3,3$%.

Во многих случаях (например, для спутников с большими панелями солнечных батарей) для оценки сил солнечного давления достаточно предположить что нормаль к поверхности $\vec{n}$ совпадает с направлением на Солнце. Тогда выражение (\ref{sun_accel}) упрощается до

$$ \ddot{\vec{r}} = -\nu P_\odot C_R \frac{A}{m} \frac{\vec{r}_\odot}{r_\odot^3} \text{а.е.}^2 , $$

где $C_R$ — коэффициент солнечного давления.

Функция тени $\nu$ характеризует условия освещенности спутника Солнцем и принимает значения в диапазоне от 0 до 1 ($\nu=1$ соответствует освещенному спутнику; $\nu=0$ — спутнику, находящемуся в тени; значения $0<\nu<1$ соответствуют нахождению спутника в области полутени).

Простейшую функцию тени получим, предположив, что спутник может быть либо полностью освещен Солнцем, либо скрыт от него в тени планеты. Обозначим через $s$ проекцию радиуса-вектора спутника $\vec{r}$ на направление $\vec{e}_\odot$. Тогда спутник освещен, если он находится между Солнцем и планетой или «позади» планеты, но не закрыт ею от Солнца (рис.). Этот результат можно записать в виде следующей функции

\begin{equation} \nu = \left\lbrace \begin{array}{ll} 1,& s > 0,\\ 1,& (s < 0)\ \text{и}\ |\vec{r} + s\cdot\vec{e}_\odot| > R_E,\\ 0,& \text{в остальных случаях}, \end{array} \right. \label{eq:simple_shadow} \end{equation}

где $s = (\vec{r}\cdot\vec{e}_\odot)$;

$R_E$ — радиус-Земли.

$К выводу выражения (\ref{eq:simple_shadow}) для функции тени, без учета полутени$

К выводу выражения (\ref{eq:simple_shadow}) для функции тени, без учета полутени

Для более точных расчетов необходимо учесть попадание спутника в область полутени. Значения $\nu$ в этом случае зависят от степени перекрытия Солнца затеняющим телом (Землей или Луной). Рассмотрим модель в которой Солнце и затеняющее тело представлены круговыми дисками. Схема затенения Солнца сферическим телом показана на рис..

Схема затенения Солнца сферическим телом

Согласно ей, видимый радиус Солнца $a$ равен

$$ a = \arcsin \frac{R_\odot}{|\vec{r}_\odot - \vec{r}|}, $$

где $R_\odot$ — радиус Солнца (см. пояснения на рис.).

Видимый радиус затеняющего тела $b$ равен

$$ b = \arcsin \frac{R_B}{s}, $$

где $R_B$ — радиус затеняющего тела;

$\vec{s} = \vec{r}-\vec{r}_B$, $\vec{r}_B$ — геоцентрический радиус-вектор центра затеняющего тела, $s = |\vec{s}|$ — расстояние между спутником и затеняющим телом.

$$ c = \arccos \frac{-\vec{s}\cdot (\vec{r}_\odot - \vec{r})}{s|\vec{r}_\odot - \vec{r}|} $$

— видимое расстояние между центрами Солнца и затеняющего тела.

К определению видимого радиуса Солнца

Тогда площадь перекрытого участка солнечного диска равна

$$ A = A_{CFC'} + A_{CDC'}. $$

При этом

\begin{equation} |a-b| < c < a+b. \label{occ_condition} \end{equation}

Область перекрытия можно выразить как

\begin{equation} A = 2(A_{BCF}-A_{BCE}) + 2(A_{ACD}-A_{ACF}). \label{occ_area} \end{equation}

Вводя обозначения $\overline{AE}=x$, $\overline{CE}=y$, $\angle CAE = \alpha$, найдем

$$ \renewcommand{\arraystretch}{1.3} \begin{array}{l} A_{ACD} = \frac{1}{2}\alpha a^2, \\ A_{ACE} = \frac{1}{2}xy. \end{array} $$

Вычислив аналогично остальные площади из (\ref{occ_area}), получим в результате

$$ A = a^2 \arccos(x/a) + b^2 \arccos((c-x)/b) - cy, $$

где

$$ x = \frac{c^2+a^2-b^2}{2c}, \qquad y = \sqrt{a^2-x^2}. $$

Доля незатененной части Солнца равна

$$ \nu = 1 - \frac{A}{\pi a^2}. $$

Если условие (\ref{occ_condition}) не выполняется, то реализуется одна из следующих возможностей:

затенения не происходит (при $a+b \leq c$);
затенение будет полным (когда $c < b-a$ при $a<b$);
затенение будет частичным, но максимально возможным (когда $c < a-b$ при $a>b$).

Давление излучения Земли

Помимо давления солнечного излучения, на спутник действует также давление от излучения, испускаемого Землей. В этом излучении выделяют две составляющие: коротковолновое оптическое излучение и длинноволновое инфракрасное излучение. Возмущающее ускорение от излучения Земли постепенно уменьшается с увеличением высоты полета. С одной стороны, поток излучения подчиняется закону обратных квадратов, то есть убывает с высотой. С другой стороны, уменьшение потока частично компенсируется тем, что с увеличением высоты растет площадь поверхности Земли, освещающей спутник. Амплитуда типичного альбедо-ускорения для спутников на низких околоземных орбитах составляет от 10% до 35% от ускорения, вызываемого давлением непосредственно солнечного излучения [1].

Программное обеспечение расчета орбитального движения спутника

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

Рассмотрим программное обеспечение расчета орбитального движения космических аппаратов, которое распространяется свободно и с открытым исходным кодом.

GMAT (General Mission Analysis Tool) разрабатывается коллективом сотрудников Центра им. Годдарда NASA и группой частных компаний. Он предназначен для моделирования, анализа и оптимизации траекторий космических аппаратов в различных режимах полета, начиная с низких околоземных орбит и движения вокруг Луны и заканчивая межпланетными траекториями и другими дальними космическими полетами.

Анализ предстоящего полета в GMAT начинается с создания ресурсов (resources), таких как: космический корабль, пропагатор (propagator, то есть метод расчета орбиты), оптимизатор (optimizer, способ оптимизации орбиты) и т. п. Параметры ресурсов задаются в зависимости от задач полета. Так, можно указать типы и количество реактивных двигателей, объем баков с горючим, вид горючего и мн. др.

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

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

Orekit (ORbits Extrapolation KIT) — низкоуровневая библиотека для решения задач динамики космического полета.

Orekit предназначен для создания точных и эффективных низкоуровневых компонентов, используемых в программных приложениях из области динамики космического полета. Разработчики стремились сделать из него набор «строительных блоков», который можно было бы легко использовать в самых разных ситуациях — от простых оценочных расчетов, до анализа межпланетных перелетов. Библиотека используется в Национальном центре космических исследований (CNES, Франция).

Код Orekit написан на Java, и зависит только от наличия Java Standard Edition и математической библиотеки Hipparchus.

Scilab — это пакет прикладных программ для научных и инженерных расчетов, представляющий собой наиболее полную общедоступную альтернативу проприетарному пакету MATLAB. В состав пакета входит Xcos — инструмент для создания блочных моделей, аналог Simulink из пакета MATLAB. Scilab разработан в INRIA (Франция) и в настоящее время развивается консорциумом Scilab Consortium.

Scilab поддерживает несколько десятков пакетов-расширений (toolboxes), которые можно загрузить из репозитория ATOMS (AuTomatic mOdules Management for Scilab). Среди них, для расчета орбитального движения космического аппарата предназначены CelestLab (разработки CNES) и Aerospace Blockset — набор блоков для Xcos, основанных на CelestLab. Общим недостатком этих расширений является то, что в них используются лишь простейшие модели плотности атмосферы. Одним из способов решения этой проблемы является использование моделей, имеющихся в библиотеке Orekit (использовать Java-объекты в Scilab можно, установив расширение JIMS).

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

Montenbruck O., Gill E. Satellite Orbits. Models, Methods and Applications. Heidelberg : Springer-Verlag, 2005.
Бордовицына Т. В., Авдюшев В. А. Теория движения искусственных спутников Земли. Аналитические и численные методы. Томск : Изд-во Том. ун- та, 2007.
Markley F. L., Crassidis J. L. Fundamentals of Spacecraft Attitude Determination and Control. New York : Springer, 2014.