Задача (формулировка с математической моделью и критериями).
Дано: начальные осреднённые навигационные элементы (или позиции и скорости) двух околоземных объектов $O_1,O_2$ в момент $t_0$, их неопределности (ковариационные матрицы) $C_{0,1},C_{0,2}$, физические параметры (масса $m$, площадь поперечного сечения $A$, коэффициент лобового сопротивления $C_d$, параметры теплового ответа и спина, задающие эффект Ярковского), модель плотности верхних слоёв атмосферы $\rho (h)$. Требуется: 1) критерии долгосрочной стабильности орбит; 2) оценка вероятности столкновения за заданный интервал $T$, учитывая атмосферное сопротивление и эффект Ярковского.
Модель движения (суммарный режим):
$$\ddot{\mathbf{r}}=-\mu \frac{\mathbf{r}}{r^3}+{\mathbf{a}}_{\mathrm{pert}}+{\mathbf{a}}_{\mathrm{drag}}+{\mathbf{a}}_{\mathrm{Y}},$$ где $\mu$ — гравитационный параметр Земли, ${\mathbf{a}}_{\mathrm{pert}}$ — прочие гравитационные возмущения (включая нецентральность, Луна, Солнце), ${\mathbf{a}}_{\mathrm{drag}}$ — атмосферное сопротивление, ${\mathbf{a}}_{\mathrm{Y}}$ — вектор ускорения от эффекта Ярковского.
Модель атмосферного сопротивления (в аэродинамической форме):
$${\mathbf{a}}_{\mathrm{drag}}=-\frac{1}{2}{C}_d\frac{A}{m}\rho (h)v_{\mathrm{rel}}{\mathbf{v}}_{\mathrm{rel}},$$ где ${\mathbf{v}}_{\mathrm{rel}}$ — скорость относительная атмосфере, $h=r-R_{\oplus }$. (Для малых тел в ВБО $\rho (h)$ экспоненциальна и даётся моделью; $\rho$ и $C_d$ — с погрешностями.)
Упрощённая параметризация эффекта Ярковского (в орбитальной трёхосевой системе: радиальная $\hat{\mathbf{r}}$, тангенциальная $\hat{\mathbf{t}}$, нормальная $\hat{\mathbf{h}}$):
$${\mathbf{a}}_{\mathrm{Y}}=A_1\hat{\mathbf{r}}+A_2\hat{\mathbf{t}}+A_3\hat{\mathbf{h}},$$ где $A_i$ — параметры, зависящие от тепловых свойств, размера, спина; в задаче они трактуются как неизвестные величины с априорным распределением (оценки и дисперсии).
Связь возмущений с изменением орбитальных элементов: используйте уравнения Гаусса. В частности, компонентное выражение для полуоси (пример):
$$\dot{a}=\frac{2}{n\sqrt{1-e^2}}\left(e\sin fR+\frac{p}{r}T\right),$$ где $R,T,N$ — проекции дополнительного ускорения на радиальную, тангенциальную и нормальную оси, $n=\sqrt{\mu /a^3}$, $p=a(1-e^2)$, $f$ — истинная аномалия. Аналогично получаются $\dot{e},\dot{i},\dot{\Omega },\dot{\omega }$ (см. справочную формулу Гаусса для полного набора).
Критерии долгосрочной устойчивости (рекомендуемые количественные тесты):
1. Сходимость/линейная устойчивость:
- Посчитать наибольший Ляпуновский показатель ${\lambda }_{\max }$ для динамической пары (или одиночных орбит) по фазовым траекториям; требование устойчивости: характеристическое время разлога $T_L=1/{\lambda }_{\max }$ должно быть значительно больше целевого горизонта $T$ (например $T_L\gg T$).
2. Сохранение топологии орбит:
- Минимальное расстояние пересечения орбит (MOID) должно оставаться больше гарантийного запаса $\Delta$ на весь интервал: ${\min }_{t\in [0,T]}\mathrm{MOID}(t)>\Delta$ (обычно $\Delta$ = сумма радиусов + запас).
3. Секулярные дрейфы:
- Оценить среднюю скорость изменения полуоси $\overline{da/dt}$ из усреднённых уравнений (учитывая $A_2$ для Ярковского и среднее значение тангенциальной составляющей сопротивления). Требование стабильности: ∣Δa∣=∣∫0Tda/dt‾ dt∣<ϵa|\Delta a| = \left|\int_0^T \overline{da/dt}\,dt\right| < \epsilon_a∣Δa∣= ∫0T da/dt dt <ϵa (порог ${ϵ}_a$ задаётся).
4. Резонансные и переходные условия:
- Проверить пересечение резонансных зон с Землёй/планетами; если орбита пересекает зону резонанса, устойчивость снижается резко.
Оценка вероятности столкновения (алгоритм):
1. Задайте статистику входных неопределённостей: начальные ковариации $C_{0,i}$ и априорные распределения для $A_1,A_2,A_3$ и параметров атмосферы.
2. Стратегии расчёта:
- Монте‑Карло (надёжно, но дорого): сэмплировать наборы (начальные состояния + параметры Ярковского + параметры атмосферы), численно интегрировать уравнения движения для двух тел на $[t_0,t_0+T]$, фиксировать события близкого сближения или столкновения $(\Vert {\mathbf{r}}_1-{\mathbf{r}}_2\Vert <R_1+R_2)$. Оценка вероятности: частота столкновений в выборке.
- Линейная аппроксимация около ожидаемого времени сближения (метод Ошика/Чесли): найти моменты потенциального узла (времена пересечения узлов орбит), вычислить среднее относительное положение ${\mathbf{\mu }}_{\perp }$ в трансверсальной плоскости и ковариацию ${\Sigma }_{\perp }$. Тогда вероятность столкновения при нормальном приближении:
$$P_{\mathrm{coll}}={\iint }_{|{\mathbf{r}}_{\perp }|<R}\mathcal{N}({\mathbf{r}}_{\perp };{\mathbf{\mu }}_{\perp },{\Sigma }_{\perp })d^2{r}_{\perp },$$ где $R=R_1+R_2$. Интеграл есть двумерная нормальная вероятность попадания в круг; вычисляется численно или через преобразование в канонический вид.
3. Учёт неневязок:
- При больших нелинейностях/длинных горизонтах избегать только-линейных приближений; предпочесть методику частичной линейности с реанализом (detect potential encounters, затем локальный МС вокруг этих моментов).
4. Практические критерии достоверности:
- Если для малого числа выборок вероятность $P_{\mathrm{coll}}$ превосходит заданный порог, требуется уточнённый анализ с увеличением числа итераций и/или наблюдений для уменьшения ковариации.
Рекомендованная пошаговая процедура для решения задачи на практике:
1. Сформировать модель: задать ${\mathbf{a}}_{\mathrm{drag}}$ и стохастическую модель ${\mathbf{a}}_{\mathrm{Y}}$ (априорные распределения).
2. Продвинутый численный интегратор (адаптивный, учёт атмосферы и негравитаций) для ансамблевой интеграции.
3. Вычислить MOID(t), найти потенциальные моменты близких сближений; для них применить линейную оценку вероятности и/или локальный МС.
4. Оценить ${\lambda }_{\max }$ и секулярные $\Delta a,\Delta e$ за интервал $T$ — дать вердикт устойчивости.
5. Отчёт: дать $P_{\mathrm{coll}}$ с доверительными интервалами и условиями чувствительности (на какие параметры результат сильнее зависит: $A_2$, площадь/масса, плотность атмосферы).
Замечания по параметрам и масштабу:
- Эффект Ярковского даёт медленный секулярный дрейф $a$ (важен для десятилетних и более горизонтов) и сильно зависит от спина и размеров (для тел 10–100 m возможно существенен).
- Атмосферный драг значим для глубоко перигейных орбит и при низких высотах (менее ~300 km); для больших высот эффект мал, но непредсказуемые всплески плотности могут давать скачки.
- Для коротких горизонтов (годы) и небольших неопределённостей чаще достаточно локальной линейной оценки; для десятилетий/веков предпочтителен ансамбль/сэмплинг.
Итог: сформулируйте модель уравнения движения с ${\mathbf{a}}_{\mathrm{drag}}$ и стохастической ${\mathbf{a}}_{\mathrm{Y}}$, используйте уравнения Гаусса для перехода к элементам и вычисления секулярных дрейфов; критерии устойчивости — сравнительный анализ Ляпунова‑времени, сохранение MOID и малость секулярных изменений $|\Delta a|,|\Delta e|$; вероятность столкновения вычисляется методом ансамблей (MC) или через линейную нормальную аппроксимацию относительной ковариации в момент потенциального сближения (формула интеграла двумерной нормали над кругом приведена выше).