Анализ модели движения космического аппарата на околоземной орбите

Полный текст

Введение

Спутниковая навигация (СН) стала неотъемлемой частью современной жизни. От мобильных приложений для ориентации в городе до точного позиционирования в морской, авиационной и космической отраслях СН обеспечивает надежное и точное определение местоположения, скорости и времени [1]. Эти достоинства особенно полезны для автономной навигации, автоматического пилотирования, картографии и других приложений, которые только начинают активно развиваться [2].

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

1. сигнальные помехи: спутниковые сигналы, передаваемые навигационными системами, могут подвергаться различным помехам, таким как многолучевое распространение, шумы и блокирование сигнала. Это может вызывать искажение сигналов и снижать точность определения местоположения;
2. географические ограничения: в некоторых местах на Земле, таких как глубокие долины, горные районы или густо застроенные города, сигналы от спутников могут быть ослаблены или полностью блокированы. Это приводит к проблемам с получением надежного сигнала навигации. Похожая проблема наблюдается в космическом пространстве для космических аппаратов (КА) на геостационарной орбите (ГСО), когда имеющегося числа видимых навигационных космических аппаратов (НКА) недостаточно для точного определения состояния КА;
3. зависимость от состояния НКА: СН полностью зависит от функционирования глобальных навигационных спутниковых систем (ГНСС) [3]. Любые погрешности, технические сбои, сбои в спутниковых сетях или намеренные акты вмешательства могут привести к недоступности или неадекватной работе СН.

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

Мгновенное (одномоментное) КСВ-решение по ГНСС – это метод, позволяющий определять положение объекта (вектор состояния объекта – координаты, скорость и время) в реальном времени. Этот метод удовлетворяет требованиям большинства приземных пользователей, а совместно с инерциальной системой позволяет преодолеть большинство проблем СН.

Мгновенное КСВ-решение для объектов на околоземных орбитах (ОО) аналогично «приземному», при этом особенностью орбитального движения, отличающей его от движения приземных объектов, является высокая предсказуемость – высокая точность прогноза движения, определяемого через параметры гравитации Земли, Луны, Солнца, сопротивления атмосферы, светового давления и др.

Для навигации на ОО прогноз дает возможность:

1. рассчитывать в реальном времени движение на интервалах потери слежения за сигналами НКА (подобно инерциальной навигационной системе), а также прогнозировать (предсказывать) движение в будущем;
2. рассчитывать по мгновенным координатам и скоростям (КС) среднюю на интервале орбиту (СрО), которая может быть точнее каждого из исходных мгновенных КС-решений;
3. использовать рекурсивные алгоритмы [4–7], позволяющие в режиме реального времени учитывать различные источники погрешностей, такие как ошибки измерений, шумы, межпространственные помехи, систематические ошибки в координатах НКА.

Современные исследования в области моделирования орбитальной динамики КА на околоземных орбитах сталкиваются с необходимостью балансировки между точностью прогнозирования и вычислительной сложностью, обусловленной спецификой прикладных задач. В работах [8–9] представлены упрощенные модели, ориентированные на краткосрочное прогнозирование траекторий объектов на низкой круговой (НКО) и средневысокой околоземной орбитах (СОО), что актуально для задач наведения антенн с широким углом охвата. В частности, в [8] продемонстрировано применение модели SGP (Simplified General Perturbations) для расчета азимута и угла места КА на основе данных TLE (Two-Line Elements), с оценкой погрешностей относительно более точной модели. Исследование [9] предлагает алгоритм определения положения и ориентации Земли относительно КА, что критично для симуляции процессов стыковки. В отличие от указанных работ, исследование [10] фокусируется на траекториях межпланетных миссий (на примере экспедиции к астероиду Апофис), где набор возмущающих факторов включает гравитационные влияния удаленных небесных тел, несферичность астероида и давление солнечного света. В работах [11; 12] рассматриваются методы оптимизации траекторий, включая псевдоспектральный метод Гаусса и релятивистские модели распространения, однако их вычислительная сложность ограничивает применимость в режиме реального времени для ресурсоограниченных систем.

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

Следует отметить, что большинство существующих исследований игнорируют ГСО, ограничиваясь НКО и СОО. Целью данной работы является комплексное описание модели движения КА, реализованной в навигационных модулях производства АО «КБ НАВИС», с акцентом на анализ влияния погрешностей параметров движения на точность прогнозирования для различных классов орбит. Представлены результаты численного моделирования отклонений траекторных параметров, а также данные летных испытаний, подтверждающие необходимость и возможность учета динамических возмущений в режиме реального времени.

1. Физическая модель прогноза вектора состояния КА

В качестве системы координат в целях геодезического обеспечения орбитальных полётов и решения навигационных задач были выбраны WGS 84 и ПЗ-90. Обе характеризуются тем, что начало системы расположено в центре масс Земли. Ось координат Z в соответствии с рекомендациями международной службы вращения Земли (IERS) направлена к среднему северному полюсу. Ось координат X лежит в плоскости земного экватора этой же эпохи, образуя пересечение с плоскостью начального меридиана, установленного все той же IERS, и определяет положение нуль-пункта принятой системы счета. Ось Y дополняет систему координат до правой.

Разница двух систем координат состоит в дальнейшем только в значениях набора констант. В основном используется WGS 84, поэтому приведен пример для этой системы координат:

• угловая скорость вращения Земли ${\mathrm{\omega }}_{\text{з}}$ = 7,292115×10^–5 рад/с;
• радиус Земли $R_{\text{з}}$ (большая полуось $a_{\text{e}}$) = 6378137 м;
• гравитационная постоянная µ = 398600,4418×10⁹ м³/с²;
• эксцентриситет орбиты Земли 6,69437999013×10^–3.

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

1. сила притяжения сферической и однородной модели Земли;
2. центробежная сила и сила Кориолиса, обусловленные вращением Земли и движением КА;
3. аномальная сила притяжения Земли, вызванная отличием Земли от сферической и однородной массы [14];
4. сопротивление атмосферы;
5. притяжение Солнца;
6. притяжение Луны;
7. давление солнечного света.

2. Математическая постановка задачи прогноза вектора состояния КА

2.1. Общая система уравнений

Прогнозирование параметров движения КА на ОО от момента времени t0 до момента t выполняется путем численного интегрирования следующих дифференциальных уравнений движения в прямоугольной геоцентрической системе координат (ГСК) для определения координат положения (x, y, z) и составляющих скорости (V_x, V_y, V_z):

$\frac{dx}{dt}=V_x,$ (1)

$\frac{dy}{dt}=V_y,$ (2)

$\frac{dz}{dt}=V_z,$ (3)

$\frac{d{V}_x}{dt}=-\frac{\mathrm{\mu }}{r^3}\cdot x+{\mathrm{\omega }}_3^2\cdot x+2\cdot {\mathrm{\omega }}_3\cdot V_y+{\ddot{x}}_a+{\ddot{x}}_s+{\ddot{x}}_c+{\ddot{x}}_m+{\ddot{x}}_l,$ (4)

$\frac{d{V}_y}{dt}=-\frac{\mathrm{\mu }}{r^3}\cdot y+{\mathrm{\omega }}_3^2\cdot y+2\cdot {\mathrm{\omega }}_3\cdot V_x+{\ddot{y}}_a+{\ddot{y}}_s+{\ddot{y}}_c+{\ddot{y}}_m+{\ddot{y}}_l,$ (5)

$\frac{d{V}_z}{dt}=-\frac{\mathrm{\mu }}{r^3}\cdot z+{\ddot{z}}_a+{\ddot{z}}_s+{\ddot{z}}_c+{\ddot{z}}_m+{\ddot{z}}_l,$ (6)

где $r=\sqrt{x^2+y^2+z^2}$ – радиус-вектор положения КА; ${\ddot{x}}_a,{\ddot{y}}_a,{\ddot{z}}_a$ – составляющие вектора ускорения, вызванные отличием модели Земли от сферической; ${\ddot{x}}_s,{\ddot{y}}_s,{\ddot{z}}_s$ – составляющие вектора ускорения, вызванные аэродинамическим торможением КА в атмосфере Земли; ${\ddot{x}}_c,{\ddot{y}}_c,{\ddot{z}}_c$ – составляющие вектора ускорения, вызываемые притяжением Солнца; ${\ddot{x}}_m,{\ddot{y}}_m,{\ddot{z}}_m$ – составляющие вектора ускорения, вызываемые притяжением Луны; ${\ddot{x}}_l,{\ddot{y}}_l,{\ddot{z}}_l$ – составляющие вектора ускорения, обусловленные давлением солнечного света.

Сила притяжения сферической и однородной модели Земли соответствует первому слагаемому в уравнениях (4)–(6). Центробежная сила соответствует второму, а сила Кориолиса третьему слагаемому формул (4)–(5).

2.2. Определение ускорений КА, обусловленных геопотенциалом

Составляющие вектора ускорения КА, вызванные отличием модели Земли от сферической, сначала определяются в топоцентрической системе координат (ТСК) через коэффициенты геопотенциала [14], после чего преобразуются в ГСК-составляющие ${\ddot{x}}_a,{\ddot{y}}_a,{\ddot{z}}_a$:

$\left(\begin{array}{c}{\dot{x}}_a\\ {\dot{y}}_a\\ {\dot{z}}_a\end{array}\right)=\left(\begin{array}{ccc}\frac{-x}{r}& -\frac{x\cdot z}{r\cdot r_1}& -\frac{y}{r_1}\\ -\frac{y}{r}& \frac{y\cdot z}{r\cdot r_1}& \frac{x}{r_1}\\ -\frac{z}{r}& \frac{r_1}{r_1}& 0\end{array}\right)\times \left(\begin{array}{c}{\dot{U}}_r\\ {\dot{U}}_{\varphi }\\ {\dot{U}}_{\lambda }\end{array}\right)$, (7)

где ${\dot{U}}_r,{\dot{U}}_{\varphi },{\dot{U}}_{\lambda }.$ – составляющие вектора ускорения в ТСК: по радиус-вектору r, по широте ϕ и долготе λ; $r_1=\sqrt{x^2+y^2}$.

Составляющие вектора ускорения ${\ddot{U}}_r,{\ddot{U}}_{\varphi },{\ddot{U}}_{\lambda }$ определяются как суммы тригонометрического ряда:

$\begin{array}{l}{\ddot{U}}_r=\mathrm{\mu }/r^2\sum _{n=2}^N{(a_e/r)}^n(n+1)\sum _{m=0}^n{A}_{nm}{\overline{P}}_{nm},\\ {\ddot{U}}_{\varphi }=\mathrm{\mu }/r^2\sum _{n=2}^N{(a_e/r)}^n\sum _{m=0}^n{A}_{nm}\cdot {{\overline{P}}^{\text{'}}}_{nm},\\ {\ddot{U}}_{\lambda }=\frac{\mathrm{\mu }}{r^2\cdot \cos \varphi }\sum _{n=2}^N{(a_e/r)}^n\sum _{m=0}^{n}m{B}_{nm}{\overline{P}}_{nm},\end{array}$

где

$\begin{array}{l}\sin \varphi =z/r,\\ \cos \varphi =r_1/r,\\ A_{nm}={\overline{C}}_{nm}\cdot \cos (m\cdot \lambda )+{\overline{S}}_{nm}\cdot \sin (m\cdot \lambda ),\\ B_{nm}={\overline{S}}_{nm}\cdot \cos (m\cdot \lambda )-{\overline{C}}_{nm}\cdot \sin (m\cdot \lambda ).\end{array}$

${\overline{C}}_{nm},{\overline{S}}_{nm}$ – нормированные коэффициенты разложения потенциала Земли.

Тригонометрические функции углов, кратных l, вычисляются по формулам

$\begin{array}{l}\sin (m\cdot \lambda )=\sin \left(\lambda \right)\cdot \cos \left(\right[m-1]\cdot \lambda )+\cos \left(\lambda \right)\cdot \sin \left(\right[m-1]\cdot \lambda ),\\ \cos (m\cdot \lambda )=\cos \left(\lambda \right)\cdot \cos \left(\right[m-1]\cdot \lambda )-\sin \left(\lambda \right)\cdot \sin \left(\right[m-1]\cdot \lambda ),\\ \cos \lambda =x/r_1,\sin \lambda =y/r_1.\end{array}$

Нормированные полиномы Лежандра ${\overline{P}}_{nm}$ являются функциями от $\sin \left(\mathrm{\varphi }\right)$, $\cos \left(\mathrm{\varphi }\right)$ и определяются по следующим рекуррентным зависимостям:

$\begin{array}{l}{\overline{P}}_{nm}=\sqrt{\frac{2\cdot n+1}{2\cdot n}}\cdot \cos \left(\varphi \right)\cdot {\overline{P}}_{n-\mathrm{1,}m-1},\text{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}n}=\text{m}>\text{1,}\\ {}_{\text{nm}}=\sqrt{\frac{2\cdot n+1}{n^2-m^2}}[\sqrt{2\cdot n-1}\cdot \sin (\varphi )\cdot {\overline{P}}_{n-\mathrm{1,}m}-\sqrt{\frac{(n-m-1)\cdot (n+m-1)}{2\cdot n-3}}\cdot {\overline{P}}_{n-\mathrm{2,}m}],\text{\hspace{1em}n}>\text{m,}\\ {}_{\text{nm}}=0,\text{\hspace{1em}m}>\text{n\hspace{0.33em}.}\end{array}$

Начальные значения: ${\overline{P}}_{00}=\mathrm{1,}{\overline{P}}_{10}=\sqrt{3}\mathrm{sin\varphi },P_{11}=\sqrt{3}\mathrm{cos\varphi }$.

Производные от нормированных полиномов Лежандра ${\overline{P}}_{nm}$ по $\mathrm{\varphi }$ определяются по следующим рекуррентным зависимостям:

$\begin{array}{l}{{\overline{P}}^{\text{'}}}_{nm}=\sqrt{\frac{2\cdot n+1}{2\cdot n}}\cdot \left(\cos \right(\varphi )\cdot {{\overline{P}}^{\text{'}}}_{n-\mathrm{1,}m-1}-\sin (\varphi )\cdot {\overline{P}}_{n-\mathrm{1,}m-1}),\text{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}n}=\text{m}>\text{1,}\\ {}_{\text{nm}}=\sqrt{\frac{2\cdot n+1}{n^2-m^2}}\cdot \left[\sqrt{2\cdot n-1}\right(\sin \left(\varphi \right)\cdot {{\overline{P}}^{\text{'}}}_{n-\mathrm{1,}m}+\cos \left(\varphi \right)\cdot {\overline{P}}_{n-\mathrm{1,}m})-\sqrt{\frac{(n-m-1)\cdot (n+m-1)}{2\cdot n-3}}\cdot {{\overline{P}}^{\text{'}}}_{n-\mathrm{2,}m}],\text{\hspace{1em}n}>\text{m,}\\ {}_{\text{nm}}=0,\text{\hspace{1em}m}>\text{n\hspace{0.33em}.}\end{array}$

Начальные значения: ${{\overline{P}}^{\text{'}}}_{00}=\mathrm{0,}{{\overline{P}}^{\text{'}}}_{10}=\sqrt{3}\cdot \cos \left(\mathrm{\varphi }\right),{P^{\text{'}}}_{11}=-\sqrt{3}\cdot \sin \left(\mathrm{\varphi }\right).$

2.3. Определение ускорений КА, обусловленных торможением в атмосфере

Составляющие вектора ускорения ${\ddot{x}}_s,{\ddot{y}}_s,{\ddot{z}}_s$, вызванные аэродинамическим торможением КА в атмосфере Земли [14], вычисляются по формулам

$\begin{array}{l}{}_{\text{s}}=-{\text{S}}_{\text{b}}\cdot \mathrm{\rho }\cdot \text{V}\cdot {\text{V}}_{\text{x}}\text{, }\\ {}_{\text{s}}=-{\text{S}}_{\text{b}}\cdot \mathrm{\rho }\cdot \text{V}\cdot {\text{V}}_{\text{y}}\text{,}\\ {}_{\text{s}}=-{\text{S}}_{\text{b}}\cdot \mathrm{\rho }\cdot \text{V}\cdot {\text{V}}_{\text{z}}\text{,}\end{array}$ (8)

где ρ – плотность верхней атмосферы из национальной стандартной модели (ГОСТ Р 25645.166–2004); $V=\sqrt{{V_x}^2+{V_y}^2+{V_z}^2}$ – модуль вектора скорости КА; ${\text{S}}_{\text{b}}=C_x{S}_m/2m$ – баллистический коэффициент; m – масса КА; $C_x$ – аэродинамический коэффициент лобового сопротивления; $S_m$ – площадь Миделева сечения КА.

Значения аэродинамического коэффициента $C_x$ и площади Миделева сечения $S_m$ определятся размером и формой конкретного КА.

2.4. Определение ускорений КА, обусловленных притяжением Солнца

Составляющие вектора ускорения ${\overrightarrow{\ddot{r}}}_c$ = (${\ddot{x}}_c$, ${\ddot{y}}_c$, ${\ddot{z}}_c$), вызванные возмущениями от Солнца [14], определяются по следующим формулам:

${\overrightarrow{\ddot{r}}}_c=\frac{{\mathrm{\mu }}_c\cdot ({\overrightarrow{r}}_c-\overrightarrow{r})}{{\left|{\overrightarrow{r}}_c-\overrightarrow{r}\right|}^3}-\frac{{\mathrm{\mu }}_c\cdot {\overrightarrow{r}}_c}{{\left|{\overrightarrow{r}}_c\right|}^3}$,

где ${\mathrm{\mu }}_c$ – произведение гравитационной постоянной на массу Солнца; ${\overrightarrow{r}}_c$ – вектор положения Солнца относительно Земли; $\overrightarrow{r}$ – вектор положения объекта относительно Земли.

Расчет эклиптических и геоцентрических координат Солнца [15] использует в качестве констант средние элементы орбиты Солнца, определенные на эпоху J2000. Оцениваемая погрешность расчета – не более 10 угловых мин.

Прямоугольные экваториальные координаты Солнца вычисляются при помощи радиус-вектора, долготы и наклона эклиптики к экватору. При расчете координат Солнца на текущую дату учитывается его эллиптическое движение. Радиус-вектор Солнца и долгота Солнца на текущую дату определяются через ряды по средней аномалии. Средняя аномалия на текущую дату определяется от средней аномалии на эпоху J2000 с учетом среднего движения на эпоху J2000 и времени в юлианских столетиях от эпохи J2000. Наклон эклиптики к экватору на текущую дату определяется от наклона эклиптики на J2000 и изменения наклона за время в юлианских столетиях от эпохи J2000.

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

Экваториальные координаты Солнца определяются по следующим формулам:

$\begin{array}{c}x=r\cdot \text{cos}\left(L\right),\\ y=r\cdot \text{sin}\left(L\right)\cdot \text{cos}\left(EPS\right),\\ z=r\cdot \text{sin}\left(L\right)\cdot \text{sin}\left(EPS\right),\end{array}$

где r – расстояние до Солнца:

$r=[\mathrm{149,619}-\mathrm{2,499}\cdot \text{cos}(M)-\mathrm{0,021}\cdot \text{cos}(2\cdot M\left)\right]\cdot {10}^9,$

M – средняя аномалия:

$M=357°\mathrm{,528}+35999°\mathrm{,05}\times T$,

L – долгота Солнца с учетом прецессии:

$L=282°\mathrm{,94}+M+[6892"\cdot \text{sin}(M)+72"\times \text{sin}(2\cdot M\left)\right]/3600+1°3972\cdot T$,

EPS – наклон эклиптики к экватору:

$EPS=23°\mathrm{,43929111}–46"\mathrm{,8150}\cdot T–0"\mathrm{,00059}\cdot T^2+0"\mathrm{,001813}\cdot T^3$,

T – время в юлианских столетиях от эпохи J2000.

Координаты Солнца (вектор Солнца ${\overline{R}}_{\text{S}}=\left[x_S{y}_S{z}_S\right]$) в прямоугольной геоцентрической СК определяются по формулам

$\begin{array}{c}{x}_S=x\cdot \text{cos}\left(Tgr\right)+y\cdot \text{sin}\left(Tgr\right),\\ y_S=y\cdot \text{cos}\left(Tgr\right)–x\cdot \text{sin}\left(Tgr\right),\\ z_S=z,\end{array}$

где $Tgr$ – гринвичский часовой угол от эпохи J2000:

$Tgr={280}^{\circ }\mathrm{,46061837504}+{360}^{\circ }\mathrm{,9856473662862}\cdot d$,

$i=Tgr/360$,

$RAD=180/\mathrm{\pi }$,

$Tgr=(Tgr–i\cdot 360)/RAD$, если $Tgr<0$, то $Tgr=Tgr+2\mathrm{\pi }$,

$d=\left(MJD–51544\right)+(FMJD\cdot 86400–dUTC)/86400–\mathrm{0,5};$

MJD – число дней от начала эпохи J2000,

FMJD число секунд от начала эпохи,

dUTC поправка времени к шкале UTC.

Расчет экваториальных координат Солнца проверялся по Астрономическому ежегоднику СССР на 1983 г. и контрольным точкам с сайта NASA. Примеры расчета представлены в табл. 1.

 

Таблица 1

Оценка ошибки расчета координат Солнца

Расчет экваториальных координат Солнца на 21.11.2010 12:00:00
Экваториальные координаты Солнца	X, км	Y, км	Z, км
Сайт NASA	–7,6318610 E+07	–1,1610898 E+08	–5,0335183 E+07
Расчет по формулам	–7,6038571 E+07	–1,1626342 E+08	–5,0402949 E+07
Ошибка расчета	0,0280039 E+07	0,0015444 E+08	0,0067766 E+07
Max ошибка	0,3269033 E+06 км=0,00218=7,5 угловых мин
Расчет экваториальных координат Солнца на 21.06.2011 12:00:00
Экваториальные координаты Солнца	X, км	Y, км	Z, км
Сайт NASA	9,9620932 E+05	1,3947997 E+08	6,0467823 E+07
Расчет по формулам	6,4872767 E+05	1,3948265 E+08	6,0468820 E+07
Ошибка расчета	3,4748165 E+05	0,0000268 E+08	0,0000997 E+07
Max ошибка	0,347493 E+06 км = 0,00232 = 8,0 угловых мин

 

2.5. Определение ускорений КА, обусловленных притяжением Луны

Составляющие вектора ускорения ${\overrightarrow{\ddot{r}}}_m$ = (${\ddot{x}}_m$, ${\ddot{y}}_m$, ${\ddot{z}}_m$), вызванные возмущениями от Луны [14], аналогично Солнцу, определяются по следующим формулам:

${\overrightarrow{\ddot{r}}}_m=\frac{{\mathrm{\mu }}_m\cdot ({\overrightarrow{r}}_m-\overrightarrow{r})}{{\left|{\overrightarrow{r}}_m-\overrightarrow{r}\right|}^3}-\frac{{\mathrm{\mu }}_m\cdot {\overrightarrow{r}}_m}{{\left|{\overrightarrow{r}}_m\right|}^3}$,

где ${\mathrm{\mu }}_m$ – произведение гравитационной постоянной на массу Луны; ${\overrightarrow{r}}_m$ = ($x_m$, $y_m$, $z_m$) – вектор положения Луны относительно Земли; $\overrightarrow{r}$ – вектор положения КА относительно Земли.

Для расчета координат Луны [15] надо знать текущее время в формате год, месяц, день, час, минута, секунда, поправка к ШВ UTC (US).

Порядок расчета и формулы следующие:

перевод текущей даты и времени во время в юлианских столетиях Т от эпохи J2000 аналогично Солнцу;

расчет средней долготы Луны el0 (в градусах):

$el0=\mathrm{218,31617}+\mathrm{481267,88088}\cdot T–\mathrm{1,3972}\cdot T^2$;

средняя аномалия Луны el (в градусах):

$el=\mathrm{134,96292}+477198/86753\cdot T$;

средняя аномалия Солнца elp (в градусах):

$elp=\mathrm{357,528}+\mathrm{35999,04944}\cdot T$;

среднее угловое расстояние f Луны от восходящего узла (в градусах):

$f=\mathrm{93,27283}+\mathrm{483202,01873}\cdot T;$

разность средних долгот Солнца и Луны (в градусах):

$d=\mathrm{297,85027}+\mathrm{445267,11135}\cdot T$

истинная эклиптическая долгота Луны (эпоха J2000) отличается от средней долготы периодическими членами Lon и dlon (в градусах):

$Lon=el0+dlon$;

$\begin{array}{l}dlon=[22640\cdot \text{sin}\left(el\right)\text{ + }769\cdot \text{sin}\left(2\cdot el\right)–4586\cdot \text{sin}\left(el–2\cdot d\right)+\\ +2370\cdot \text{sin}\left(2\cdot d\right)–668\cdot \text{sin}\left(elp\right)–412\cdot \text{sin}\left(2\cdot f\right)–212\cdot \text{sin}\left(2\cdot el–2\cdot d\right)\\ –206\cdot \text{sin}\left(el+elp–2\cdot d\right)+192\cdot \text{sin}\left(el+2\cdot d\right)-165\cdot \text{sin}\left(elp–2\cdot d\right)+148\cdot \text{sin}\left(el-elp\right)\\ –125\cdot \text{sin}\left(d\right)–110\cdot \text{sin}\left(el+elp\right)-55\cdot \text{sin}\left(2\cdot f–2\cdot d\right)]/3600\end{array}$

истинная эклиптическая широта Луны Lat (эпоха J2000) (в градусах):

$\begin{array}{l}Lat=[18520\cdot \text{sin}\left(f+Lon–el0+q\right)–526\cdot \text{sin}\left(f–2\cdot d\right)+44\cdot \text{sin}\left(el+f–2\cdot d\right)–\\ –31\cdot \text{sin}\left(–el+f–2\cdot d\right)–25\cdot \text{sin}\left(–2\cdot el+f\right)–23\cdot \text{sin}\left(elp+f–2\cdot d\right)\\ +21\cdot \text{sin}\left(–el+f\right) +11\cdot \text{sin}\left(–elp+f–2\cdot d\right)]/3600,\end{array}$

где $\begin{array}{l}q=412/3600\cdot \text{sin}\left(2\cdot f\right)+541/3600\cdot \text{sin}\left(elp\right)\end{array}$;

расстояние Rse от центра Земли до Луны (в метрах):

$\begin{array}{l}Rse=[38500–20905\cdot \text{cos}\left(el\right)–3699\cdot \text{cos}\left(2\cdot d–el\right)–2956\cdot \text{cos}\left(2\cdot d\right)–\\ -570\cdot \text{cos}\left(2\cdot el\right)+46\times \text{cos}\left(2el–2d\right)–20\times 5\text{cos}\left(elp–2d\right)-\\ -171\times \text{cos}\left(el+2d\right)+2–152\times \text{cos}\left(el+elp–2d\right)]\times 1000;\end{array}$

истинная эклиптическая долгота Луны с учетом прецессии (в градусах):

$Lon=Lon+\mathrm{1,3972}\times T ;$

наклон эклиптики к экватору Obe (в градусах):

$Obe=\mathrm{23,43929111} ;$

эклиптические координаты Луны xse, yse, zse:

$\begin{array}{l}xe=Rse\cdot \text{cos}\left(Lon\right)\cdot \text{cos}\left(Lat\right),\\ ye=Rse\cdot \text{sin}\left(Lon\right)\cdot \text{cos}\left(Lat\right),\\ ze=Rse\cdot \text{sin}\left(Lat\right);\end{array}$

экваториальные координаты Луны xse, yse, zse:

$\begin{array}{l}xse=xe\\ yse=ye\cdot \text{cos}\left(Obe\right)–ze\cdot \text{sin}\left(Obe\right), \\ zse=ze\cdot \text{cos }\left(Obe\right)+ye\cdot \text{sin}\left(Obe\right);\end{array}$

координаты Луны в прямоугольной геоцентрической СК:

$\begin{array}{l}{x}_m=xse\cdot \text{cos(}Tgr)+yse\cdot \text{sin}(Tgr),\\ y_m=yse\cdot \text{cos}\left(Tgr\right)–xse\cdot \text{sin}\left(Tgr\right),\\ z_m=zse,\end{array}$

где гринвичский часовой угол Tgr от эпохи J2000 рассчитывается по формулам аналогичным Солнцу.

2.6. Определение ускорений КА, обусловленных давлением солнечного света

Для вычисления ускорений КА, обусловленных действием светового давления [14], надо знать коэффициент эффективного отражения Crefl, зависящий от площади поперечного сечения, массы и других свойств поверхности. Составляющие вектора ускорения ${\overrightarrow{\ddot{r}}}_l$ = (${\ddot{x}}_l$, ${\ddot{y}}_l$, ${\ddot{z}}_l$), вызванные давлением солнечного света определяются по следующей формуле:

${\overrightarrow{\ddot{r}}}_l=Crefl\cdot {\left(\frac{1.4959787\cdot {10}^8}{\left|{\overrightarrow{r}}_c\right|}\right)}^2\cdot \frac{(\overrightarrow{r}-{\overrightarrow{r}}_c)}{\left|{\overrightarrow{r}}_c\right|}$, (9)

где $\left|{\overrightarrow{r}}_c\right|$ – расстояние от Солнца до КА.

2.7. Метод Рунге – Кутта четвертого порядка

Численное интегрирование системы обыкновенных нелинейных дифференциальных уравнений движения КА (1)–(6) вида ${y^{\text{'}}}_j=f_j(t,y_1,y_2,y_3,y_4,y_5,y_6)$ осуществляется методом Рунге – Кутта четвертого порядка с постоянным шагом [16] (j = 1,…, 6 – число уравнений, описывающих движение КА).

Для k-го шага и j-го уравнения обозначим

$y_{jk}=y_j\left(t_k\right),f_{jk}=f_j(t_k,{y_1}_k,{y_2}_k,{y_3}_k,{y_4}_k,{y_5}_k,{y_6}_k)$.

На шаге k+1 интегрирования значение искомых функций рассчитывается по формуле:

$y_{jk+1}=y_{jk}+h/6\cdot (k_{1j}+2\cdot k_{2j}+2\cdot k_{3j}+k_{4j})$,

где h – шаг интегрирования,

$\begin{array}{l}{k}_{1j}=f_j(t_k,{y_1}_k,{y_2}_k,{y_3}_k,{y_4}_k,{y_5}_k,{y_6}_k),\\ k_{2j}=f_j(t_k+\frac{h}{2},{y_1}_k+\frac{h}{2}\cdot k_{11},{y_2}_k+\frac{h}{2}\cdot k_{12},{y_3}_k+\frac{h}{2}\cdot k_{13},{y_4}_k+\frac{h}{2}\cdot k_{14},{y_5}_k+\frac{h}{2}\cdot k_{15},{y_6}_k+\frac{h}{2}\cdot k_{16}),\\ k_{3j}=f_j(t_k+\frac{h}{2},{y_1}_k+\frac{h}{2}\cdot k_{21},{y_2}_k+\frac{h}{2}\cdot k_{22},{y_3}_k+\frac{h}{2}\cdot k_{23},{y_4}_k+\frac{h}{2}\cdot k_{24},{y_5}_k+\frac{h}{2}\cdot k_{25},{y_6}_k+\frac{h}{2}{k}_{26}),\\ k_{4j}=f_j(t_k+h,{y_1}_k+h\cdot k_{31},{y_2}_k+h\cdot k_{32},{y_3}_k+h\cdot k_{33},{y_4}_k+h\cdot k_{34},{y_5}_k+h\cdot k_{35},{y_6}_k+h\cdot k_{36}).\end{array}$

При численном интегрировании на большом числе шагов происходит значительная потеря точности вследствие накопления ошибок округления. Считается, что накопление ошибок округления в координатах пропорционально количеству шагов интегрирования в степени 3/2:

${\epsilon }_n^i\approx 3n^{3/2}\cdot \epsilon \left(0\right)$,

где $\mathrm{\epsilon }\left(0\right)$ – точность вычисления на каждом шаге.

3. Погрешности прогноза движения КА по ОО при наличии ошибок

Приведенные выше уравнения прогноза были запрограммированы в среде MATLAB и выполнены расчеты:

• для НКО с периодом обращения вокруг Земли 1,5 ч и высотой орбиты над уровнем Мирового океана 500 км;
• СОО с периодом 12 ч и высотой 20000 км;
• ГСО с периодом 24 ч и высотой 35777 км.

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

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

3.1. Погрешности прогноза от ошибки начального положения

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

 

Таблица 2

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

НКО. Период орбиты: Т = 1,5 ч. Высота 500 км. Погрешность прогноза, м:
Модуль ошибки, м	Интервал прогноза, ч
	0,5	1,0	1,5			3,0		6,0
0,5	2,1	5,5	8,0			15,7		30,8
2	9,7	27,3	36,6			73,1		145,6
10	39,2	101,7	144,6			286,5		566,9
30	118,5	267,2	301,0			619,4		1256,9
СОО. Период орбиты: Т = 12 ч. Высота 20000 км. Погрешность прогноза, м:
Модуль ошибки, м	Интервал прогноза, ч
	0,5	1,0		1,5	3,0		6,0	12,0	24,0
0,5	0,5	0,5		0,7	1,8		5,6	8,7	17,7
2	2,0	2,2		2,5	6,0		18,2	38,6	76,3
10	10,4	12,1		15,7	38,3		120,9	202,6	407,0
30	32,2	39,4		51,8	122,9		399,0	730,2	1461,4
ГСО. Период орбиты: Т = 24 ч. Высота 35777 км. Погрешность прогноза, м:
Модуль ошибки, м	Интервал прогноза, ч
	0,5	1,0		1,5	3,0		6,0	12,0	24,0
0,5	0,5	0,5		0,6	0,9		2,6	10,0	18,7
2	2,0	2,0		2,2	2,9		7,3	29,4	55,9
10	10,1	10,6		11,3	16,3		43,3	166,2	310,6
30	30,6	32,4		35,5	53,9		145,1	555,2	1036,8

 

Таблица 3

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

НКО. Период орбиты: Т = 1,5 ч. Высота 500 км. Погрешность прогноза, м:
Модуль ошибки, м/с:	Интервал прогноза, ч
	0,5		1,0			1,5		3,0		6,0
0,0001	0,3		1,4			1,7		3,4		6,8
0,0005	1,9		6,9			8,6		17,0		33,6
0,0020	7,7		28,1			34,7		68,7		135,8
0,0050	19,8		67,2			86,1		170,2		336,3
СОО. Период орбиты: Т = 12 ч. Высота 20000 км. Погрешность прогноза, м:
Модуль ошибки, м/с	Интервал прогноза, ч
	0,5	1,0		1,5	3,0		6,0		12,0		24,0
0,0001	0,1	0,3		0,6	1,4		6,7		12,4		24,5
0,0005	0,9	1,9		3,1	9,3		34,5		65,7		130,1
0,0020	3,6	7,7		12,7	34,3		131,4		261,4		517,3
0,0050	9,1	19,2		30,6	94,5		352,4		621,9		1232,9
ГСО. Период орбиты: Т = 24 ч. Высота 35777 км. Погрешность прогноза, м:
Модуль ошибки, м/с	Интервал прогноза, ч
	0,5	1,0		1,5	3,0		6,0		12,0		24,0
0,0001	0,1	0,3		0,5	1,2		3,6		14,7		25,8
0,0005	0,9	1,8		2,8	6,4		19,6		74,6		129,4
0,0020	3,6	7,3		11,1	24,9		76,5		268,1		507,2
0,0050	9,0	18,3		28,3	64,9		180,2		625,7		1277,6

 

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

3.2. Ошибки прогноза от разного числа учтенных гармоник геопотенциала

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

 

Таблица 4

Ошибки прогноза от числа учтенных гармоник геопотенциала

НКО. Период орбиты: Т = 1,5 ч. Высота 500 км. Ошибки прогноза, м:
Число гармоник:	Интервал прогноза, ч
	0,5		1,0		1,5		3,0			6,0
32	1,1		1,4		1,4		4,4			14,4
16	5,3		14,8		24,4		19,0			25,2
8	20,8		65,1		66,7		140,2			197,2
2	100,2		236,6		236,6		451,9			626,1
СОО. Период орбиты: Т = 12 ч. Высота 20000 км. Ошибки прогноза, м:
Число гармоник:	Интервал прогноза, ч
	0,5	1,0		1,5		3,0		6,0	12,0		24,0
32	0,0	0,0		0,0		0,0		0,0	0,0		0,0
16	0,0	0,0		0,0		0,0		0,0	0,0		0,0
8	0,0	0,0		0,0		0,0		0,1	0,3		0,8
2	0,1	0,7		1,5		5,6		12,4	27,6		32,0
ГСО. Период орбиты: Т = 24 ч. Высота 35777 км. Ошибки прогноза, м:
Число гармоник:	Интервал прогноза, ч
	0,5	1,0		1,5		3,0		6,0	12,0		24,0
32	0,0	0,0		0,0		0,0		0,0	0,0		0,0
16	0,0	0,0		0,0		0,0		0,0	0,0		0,0
8	0,0	0,0		0,0		0,0		0,0	0,0		0,1
2	0,1	0,1		0,1		0,1		0,6	3,3		21,3

 

Для достижения минимальных ошибок прогноза для НКО нужно учитывать наибольшее число гармоник потенциала Земли. С ростом высоты орбиты необходимое число гармоник уменьшается, однако 8 гармоник является необходимым минимумом.

3.3. Влияние сопротивления атмосферы

Состояние атмосферы учитывалось по ГОСТ Р 25645.166–2004. Варьирование баллистического коэффициента Sb в (8) сильнее влияет на низкой орбите [17] в силу более плотной атмосферы. Поэтому дополнительно рассмотрены несколько НКО (табл. 5). Sb линейно связан с ускорением, а значит в (1)–(6) ошибка координат будет иметь такую же зависимость, несмотря на нелинейность орбиты, что показано для одной НКО.

 

Таблица 5

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

Sb	Интервал прогноза, ч
	0,5	1,0	1,5	3,0	6,0
НКО. Период орбиты: Т = 1,45 ч. Высота 200 км. Ошибки прогноза, м:
6,25е-05	1,3	11,8	36,6	143,5	570,6
6,25е-04	13,6	118,9	366,5	1437,7	5722,7
6,25е-03	136,2	1191,2	3683,4	>10000	>10000
6,25е-02	1364,6	>10000	>10000	>10000	>10000
НКО. Период орбиты: Т = 1,49 ч. Высота 288 км. Ошибки прогноза, м:
6,25е-03	7,8	67,8	190,0	766,8	3098,1
НКО. Период орбиты: Т = 1,50 ч. Высота 500 км. Ошибки прогноза, м:
6,25е-03	0,1	0,8	3,1	12,2	48,9
НКО. Период орбиты: Т = 1,68 ч. Высота 600 км. Ошибки прогноза, м:
6,25е-03	0,1	0,2	0,7	2,7	10,9
НКО. Период орбиты: Т = 1,76 ч. Высота 1000 км. Ошибки прогноза, м:
6,25е-03	0,0	0,0	0,1	0,1	0,4

 

В табл. 5 видно ослабевание влияния баллистического коэффициента с ростом высоты орбиты. Также заметно, что увеличение баллистического коэффициента на порядок приводит к увеличению ошибки прогноза на порядок.

В силу отсутствия атмосферы выше 3000 км отсутствует и ее влияние на ошибку прогноза для СОО и ГСО. Однако для НКО учет баллистического коэффициента необходим, а для эллиптических орбит, совмещающих свойства НКО и СОО, необходим в пределах указанной высоты.

3.4. Ошибки прогноза от притяжения Солнца

В оценке влияния притяжения Солнца на прогноз положения КА рассматривается 2 параметра (табл. 6):

• ошибка углового определения положения Солнца. Для того чтобы задать угловую ошибку в координатах Солнца (ошибку влияния притяжения), координаты смещались на равномерно распределенную случайную величину так, чтобы угол смещения был равен указанному значению, затем был выбран максимальный результат погрешности прогноза среди равных значений угла смещения;
• ошибка определения расстояния до Солнца из-за периодического (год) изменения радиуса орбиты Земли от 147098 до 152098 тыс. км.

 

Таблица 6

Ошибки прогноза от притяжения Солнца

НКО. Период орбиты: Т = 1,5 ч. Высота 500 км. Погрешности прогноза, м:
Параметр ошибки					Интервал прогноза, ч
					0,5		1,0		1,5		3,0		6,0
Угол, градусы			0,1		0,0		0,0		0,0		0,0		0,1
			0,5		0,0		0,1		0,1		0,1		0,1
			2		0,1		0,2		0,2		0,2		0,3
			10		0,5		1,0		1,0		1,4		2,5
Расстояние, тыс. км			20		0,0		0,0		0,0		0,0		0,0
			100		0,0		0,0		0,0		0,0		0,0
			500		0,0		0,0		0,0		0,1		0,1
			2500		0,0		0,1		0,1		0,2		0,4
Без учета ${\overrightarrow{\ddot{r}}}_c$					0,6		2,2		2,3		4,1		7,6
СОО. Период орбиты: Т = 12 ч. Высота 20000 км. Погрешности прогноза, м:
Параметр ошибки		Интервал прогноза, ч
		0,5		1,0		1,5		3,0		6,0		12,0		24,0
Угол, градусы	0,1	0,0		0,0		0,1		0,2		1,0		2,2		2,6
	0,5	0,1		0,2		0,5		1,1		4,6		8,5		14,9
	2	0,2		0,6		1,4		4,9		18,5		33,7		58,3
	10	1,0		3,9		8,5		25,2		74,2		173,4		222,3
Расстояние, тыс. км	20	0,0		0,0		0,0		0,0		0,1		0,2		0,3
	100	0,0		0,0		0,0		0,2		0,6		0,8		1,6
	500	0,0		0,1		0,2		0,8		3,2		3,9		7,9
	2500	0,1		0,4		0,9		4,1		16,4		19,8		40,6
Без учета ${\overrightarrow{\ddot{r}}}_c$		2,0		7,9		17,6		78,9		314,4		379,3		778,3
ГСО. Период орбиты: Т = 24 ч. Высота 35777 км. Ошибки прогноза, м:
Параметр ошибки		Интервал прогноза, ч
		0,5		1,0		1,5		3,0		6,0		12,0		24,0
Угол, градусы	0,1	0,021		0,084		0,2		0,7		2,2		9,9		43,3
	0,5	0,114		0,455		1,0		4,0		13,4		52,4		210,3
	2	0,401		1,578		3,5		13,1		39,9		189,2		755,1
	10	1,756		6,871		15,4		62,6		216,0		940,6		3986,0
Расстояние, тыс. км	20	0,0		0,0		0,0		0,1		0,2		0,9		0,976
	100	0,0		0,0		0,1		0,3		0,8		4,5		4,9
	500	0,0		0,2		0,4		1,3		4,5		22,8		24,6
	2500	0,2		0,8		1,8		6,5		21,6		109,5		118,1
Без учета ${\overrightarrow{\ddot{r}}}_c$		4,5		17,4		38,2		135,2		451,0		2284,5		2464,0

 

Также приведена ошибка прогноза в случае, когда влияние притяжения Солнца не учтено.

Из данных табл. 6 следует, что вариацию расстояния объект – Солнце из-за положения объекта на орбите практически можно не учитывать.

3.5. Ошибки прогноза от притяжения Луны

В оценке влияния притяжения Луны на прогноз положения КА рассматривается 2 параметра (табл. 7):

• ошибка углового определения Луны аналогично Солнцу;
• ошибка определения расстояния до Луны из-за периодического изменение радиуса лунной орбиты от 356410 до 406740 км.

Также приведена ошибка прогноза в случае, когда влияние притяжения Луны не учтено.

 

Таблица 7

Погрешности прогноза от притяжения Луны

НКО. Период орбиты: Т = 1,5 ч. Высота 500 км. Погрешности прогноза, м:
Параметр ошибки				Интервал прогноза,
				0,5		1,0		1,5					3,0			6,0
Угол, градусы		0,1		0,0		0,0		0,0					0,0			0,0
		0,5		0,1		0,1		0,1					0,1			0,2
		2		0,2		0,3		0,3					0,4			0,6
		10		1,2		1,6		1,9					3,6			6,5
Расстояние, км		200		0,0		0,0		0,0					0,0			0,0
		1000		0,0		0,1		0,1					0,1			0,2
		5000		0,1		0,3		0,3					0,5			0,9
		25000		0,3		1,4		1,5					2,6			4,9
Без учета ${\overrightarrow{\ddot{r}}}_m$				1,3		5,9		6,1					10,8			20,4
СОО. Период орбиты: Т = 12 ч. Высота 20000 км. Погрешности прогноза, м:
Параметр ошибки			Интервал прогноза, ч
			0,5		1,0		1,5			3,0		6,0			12,0			24,0
Угол, градусы	0,1		0,0		0,1		0,2			0,7		3,0			5,4			6.0
	0,5		0,2		0,6		1,1			3,0		11,9			21,9			31.8
	2		0,6		2,1		4,4			14,0		69,4			121,1			133,1
	10		2,6		9,7		20,3			73,0		255,4			570,3			570,3
Расстояние, км	200		0,0		0,0		0,1			0,4		1,2			1,4			3,0
	1000		0,0		0,2		0,4			1,8		6,2			6,9			15,2
	5000		0,2		0,8		1,9			8,6		30,2			34,0			74,5
	25000		1,2		5,0		11,4			51,3		178,9			200,1			442,6
Без учета ${\overrightarrow{\ddot{r}}}_m$			5,1		20,4		46,6			210,8		758,7			885,9			1855,5
ГСО. Период орбиты: Т = 24 ч. Высота 35777 км. Ошибки прогноза, м:
Параметр ошибки			Интервал прогноза, ч
			0,5		1,0		1,5		3,0		6,0			12,0			24,0
Угол, градусы	0,1		0,0		0,1		0,2		0,7		3,0			18,6			24,4
	0,5		0,1		0,3		0,8		2,7		12,6			92,4			140,7
	2		0,5		2,0		4,3		15,8		71,1			452,0			590,0
	10		2,8		11,1		24,5		89,3		413,5			2373,7			4052,2
Расстояние, км	200		0,0		0,0		0,1		0,3		1,2			5,0			15,7
	1000		0,0		0,2		0,4		1,5		6,1			25,1			78,0
	5000		0,2		0,9		1,9		7,7		30,7			127,2			394,8
	25000		0,3		1,2		2,6		10,3		41,2			170,8			530,3
Без учета ${\overrightarrow{\ddot{r}}}_m$			5,9		22,9		50,6		204,8		831,2			3355,3			>10000

 

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

3.6. Влияние давления солнечного света

Варьирование коэффициента эффективного отражения Crefl относительно нулевого значения дает следующие результаты (табл. 8).

 

Таблица 8

Ошибки прогноза от Crefl

НКО. Период орбиты: Т = 1,5 ч. Высота 500 км. Ошибки прогноза, м:
Crefl	Интервал прогноза, ч
	0,5	1,0		1,5		3,0		6,0
10–8	0,0	0,1		0,3		0,6		0,9
10–7	0,2	0,9		3,0		5,7		9,3
10–6	1,6	9,4		30,4		57,0		93,1
10–5	16,0	94,3		303,8		570,3		931,0
СОО. Период орбиты: Т = 12 часов Высота 20000 км. Ошибки прогноза, м:
Crefl	Интервал прогноза, часы
	0,5	1,0	1,5		3,0		6,0	12,0	24,0
10–8	0,0	0,1	0,1		0,55		2,6	8,4	16,9
ГСО. Период орбиты: Т = 24 часа Высота 35777 км. Ошибки прогноза, м:
10–8	0,0	0,1	0,1		0,6		2,0	13,2	61,2

 

При прогнозах движения по НКО следует не забывать исключать ускорения из-за светового давления на участках затенения Солнца Землей.

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

4. Использование модели прогноза при обработке реальных данных

В ходе проведения летных испытаний объекта на НКО с помощью спутниковой навигации был получен набор мгновенных КС-решений. В силу, прежде всего, погрешностей эфемеридно-временного обеспечения ГНСС последовательность мгновенных КС-решений должна представлять собой пилообразную ломаную вокруг реальной траектории. Поэтому в качестве эталона была использована СрО, рассчитанная по мгновенным КС-решениям методом МНК с использованием прогноза (для сведения отдельных КС-решений к общему моменту) и последующим «размножением» рассчитанной СрО на все исходные мгновенные КС-решения.

Параметры прогноза были выбраны по вышеприведенным в разделе 3 данным для интервала прогноза 4,5 ч, ожидаемой точности мгновенного положения 15 м (3σ) и допустимой погрешности прогноза не более 2 м:

• гармоники потенциала Земли до 32 порядка;
• учет ускорений от притяжения Солнца и Луны без погрешностей;
• баллистический коэффициент 0,000625;
• коэффициент эффективного отражения 10^-7.

График отклонений мгновенных КС-решений от рассчитанной СрО представлен на рисунке.

 

График отклонений реальных мгновенных КС-решений объекта на НКО от рассчитанной СрО (3 витка – 4,5 ч)

Graph of deviation of change of instantaneous KS-solutions of the object on NKO from the calculated SRO (3 seconds – 4.5 hours)

 

Из рисунка видно, что характер мгновенных КС-решений соответствует ожидаемому, отсутствует явное ухудшение точности со временем, что свидетельствует о корректном выборе параметров прогноза движения. Погрешность реальных КС-решений относительно рассчитанной СрО 15 м и 4 см/с (3 σ).

Заключение

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

Основными результатами работы являются:

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

Эта статья служит основой для дальнейшего изучения методов и алгоритмов, направленных на определение местоположения на основе актуальных навигационных данных. Следует отметить, что модель движения является базой, на которой в качестве навигационного алгоритма может эффективно использоваться фильтр Калмана [18], позволяя улучшить точность и надежность навигационных систем.

В перспективе возможно дальнейшее совершенствование модели и оптимизации вычислительных процессов, например за счет исследования и перехода к регулярным кватернионным уравнениям в других переменных [19].

Об авторах

Иван Александрович Соколов

АО “Конструкторское бюро навигационных систем”; Санкт-Петербургский государственный университет

Автор, ответственный за переписку.
Email: ivansokolof1997@mail.ru
ORCID iD: 0009-0008-2709-7080

аспирант, Санкт-Петербургский государственный университет; ведущий инженер лаборатории базового программного обеспечения навигационной аппаратуры пользователей космического назначения; АО «Конструкторского бюро навигационных систем»

Россия, 121170, г. Москва, ул. Кульнева, 3; 199034, г. Санкт-Петербург, Университетская наб., 7-9

Геннадий Стефанович Цеханович

АО “Конструкторское бюро навигационных систем”

Email: ggsstt49@mail.ru
ORCID iD: 0009-0007-3748-7704

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

Россия, 121170, г. Москва, ул. Кульнева, 3

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

Дополнительные файлы