Определение динамических параметров группы объектов космического мусора в зоне ГСО по данным позиционных и фотометрических измерений

全文:

ВВЕДЕНИЕ

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

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

Важную роль в получении статистически значимых результатов играет отбраковка ошибочных или недостаточно точных наблюдений. Как показано в замечательной книге Шейнина “История теории ошибок” (Шейнин, 2007), проблема отбраковки измерений восходит к работам Гаусса и Лапласа и сколько-нибудь обоснованным является так называемое правило 3σ, где σ – среднеквадратичная ошибка, полученная в процессе улучшения орбиты по совокупности наблюдений. При этом, как показывает Шейнин, в каждой области науки, работающей с наблюдениями, существуют свои правила их отбраковки, полученные из практического опыта.

В астрономической практике, как отмечает Емельянов (2019) (стр. 280) параметр отбраковки выбирают в пределах от 3 до 6σ.

Опыт работы с наблюдениями объектов космического мусора (Авдюшев и др., 2022) показывает, что для этих объектов нужен более жесткий критерий, поскольку наблюдения содержат большое количество не устраняемых систематических ошибок.

В настоящей работе мы приводим результаты поиска такого критерия, а также описание алгоритма и программы машинной отбраковки и анализа процесса улучшения орбит объектов КМ, полученные при обработке наблюдений объектов ГСО, выполненных на оборудовании Центра коллективного пользования “Терскольская обсерватория” Института астрономии РАН. Кроме того, в работе представлены результаты совместного использования позиционных и фотометрических измерений для определения динамических параметров объектов.

УЛУЧШЕНИЕ ОРБИТ

В основе алгоритма лежит численная модель движения ИСЗ (Александрова и др. 2017), приспособленная для обработки позиционных наблюдений околоземных объектов и определения параметров движения по этим наблюдениям (Авдюшев и др., 2022).

Прогнозирование значений вектора положения спутника x, заданного в инерциальной геоцентрической системе координат, связанной с экватором J2000, на моменты наблюдений реализуется численным интегрированием системы дифференциальных уравнений движения 2-го порядка:

$\frac{{\mathrm{d}}^2\mathbf{x}}{\mathrm{d}{t}^2}=\mathbf{P}(t,\mathbf{x},\dot{\mathbf{x}},\gamma ),$ (1)

с заданными в эпоху t₀ начальными условиями $\mathbf{x}={\mathbf{x}}_0\left(t_0\right),\dot{\mathbf{x}}={\dot{\mathbf{x}}}_0\left(t_0\right)$.

Здесь P – равнодействующая сил, которая включает в себя (согласно IERS Conventions, 2010): притяжение Земли как протяженного тела (с учетом приливных деформаций) до гармоник 360-го порядка; притяжение Луны и Солнца, рассматриваемых как материальные точки; влияние светового давления в рамках сферической модели спутника.

Совместно с уравнениями движения интегрируются уравнения для вычисления изохронных производных

$\frac{{\mathrm{d}}^2}{\mathrm{d}{t}^2}\frac{\partial \mathbf{x}}{\partial \mathbf{q}}=\frac{\partial \mathbf{P}}{\partial \mathbf{x}}\frac{\partial \mathbf{x}}{\partial \mathbf{q}}+\frac{\partial \mathbf{P}}{\partial \dot{\mathbf{x}}}\frac{\partial \dot{\mathbf{x}}}{\partial \mathbf{q}}+\frac{\partial \mathbf{P}}{\partial \gamma }\frac{\partial \gamma }{\partial \mathbf{q}},$ (2)

где $\mathbf{q}={({\mathbf{x}}_0,{\dot{\mathbf{x}}}_0,\gamma )}^T$ – 7-мерный вектор определяемых параметров, в который, помимо оценок начальных координат и скоростей спутника, входит параметр $\gamma$, представляющий собой отношение площади миделевого сечения спутника A к его массе m.

Представление наблюдений и вычисление матрицы условных уравнений осуществляется традиционным способом (см., например, Авдюшев и др., 2022). Нахождения оцениваемых параметров производится с помощью итерационной процедуры Гаусса–Ньютона с обращением матрицы методом сингулярного разложения (SVD – Singular Value Decomposition):

${\mathbf{q}}_{k+1}={\mathbf{q}}_{\mathbf{k}}+\mathrm{\Delta }{\mathit{q}}_{k+1},$

$\mathrm{\Delta }{\mathbf{q}}_{k+1}=\mathrm{\Delta }{\mathbf{q}}_k-{\left({{{\mathbf{p}}^{\text{'}}}_{\mathbf{q}}}^T{{\mathbf{p}}^{\text{'}}}_{\mathbf{q}}\right)}^{-1}{{{\mathbf{p}}^{\text{'}}}_{\mathbf{q}}}^T({\mathbf{p}}^O-{\mathbf{p}}^C({\mathbf{q}}_k\left)\right).$ (3)

Здесь ${\mathbf{p}}^O$ и ${\mathbf{p}}^C={\mathbf{p}}^C\left(\mathbf{q}\right)$ – 2N-мерные векторы наблюдений и их модельных представлений наблюдений, соответственно, причем

${\mathbf{p}}^C={({\alpha }_1^C\cos {\delta }_1^C,{\delta }_1^C,\dots ,{\alpha }_N^C\cos {\delta }_N^C,{\delta }_N^C)}^T;$ (4)

N – число моментов наблюдений; ${\alpha }_i^C$ и ${\delta }_i^C$ $(i=\mathrm{1,}\dots ,N)$ – топоцентрические прямое восхождение и склонение спутника, отнесенные к экватору эпохи J2000; ${{\mathbf{p}}^{\text{'}}}_{\mathbf{q}}$ – матрица системы условных уравнений.

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

${\sigma }^2=\frac{{\left|\left|{\mathbf{p}}^O-{\mathbf{p}}^C\left({\displaystyle \hat{\mathbf{q}}}\right)\right|\right|}^2}{2N-7}$, (5)

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

$\mathbf{C}={\sigma }^2{\left({{{\mathbf{p}}^{\text{'}}}_{\mathbf{q}}}^T{{\mathbf{p}}^{\text{'}}}_{\mathbf{q}}\right)}^{-1}.$ (6)

По ней находятся ошибки определяемых параметров и число обусловленности Тодда, (Todd, 1949), представляющее собой отношение наибольшего собственного числа матрицы ${{{\mathbf{\text{p}}}^{\mathbf{\text{'}}}}_{\mathbf{q}}}^T{{\mathbf{p}}^{\text{'}}}_{\mathbf{q}}$ к наименьшему числу и характеризующее степень обусловленности задачи улучшения орбиты.

ОБРАБОТКА ИЗМЕРЕНИЙ БЛЕСКА ОБЪЕКТОВ

Наблюдения объектов ГСО представляются в виде табл. 1.

 

Таблица 1. Форма представления наблюдений

Эпоха		Прямое восхождение, ч., мин., с	Склонение, °, ΄, ″	Блеск, зв. величина
Дата	Время (UTC), ч., мин., с
10 09 2020	22 59 23.50	05 23 07.44	+02 28 45.27	17.1

 

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

Фотометрические измерения дают блеск спутника, измеренный в звездных величинах. Связь между звездными величинами спутника и Солнца задается формулой (McCue и др., 1971)

$m_{\text{sat}}=m_{\text{Sun}}-2.5\lg \frac{AkF\left(\mathrm{\psi }\right)}{{r^{\text{'}}}^2}$, (7)

где A – видимая наблюдателем поверхность аппарата, освещенная Солнцем; k – спектральный коэффициент отражения (может иметь значения от 0.1 до 1); $F\left(\mathrm{\psi }\right)$ – фазовая функция, $\mathrm{\psi }$ – угол в системе Солнце–космический объект–наблюдатель; ${\mathrm{m}}_{\text{Sun}}$ – звездная величина Солнца в видимом спектральном диапазоне (–26.74); $m_{\text{sat}}$ – звездная величина околоземного космического объекта (КО); $r^{\text{'}}$ – модуль топоцентрического вектора КО. Преобразуем формулу (7) и получим соотношение для определения A

$-0.4\left(m_{\text{sat}}-m_{\text{Sun}}\right)=\lg \frac{AkF\left(\mathrm{\psi }\right)}{{r^{\text{'}}}^2}$;

$\frac{AkF\left(\text{ψ}\right)}{{r^{\text{'}}}^2}={10}^{-0.4\left(m_{\text{sat}}-m_{\text{Sun}}\right)}$;

$A=\frac{{10}^{-0.4\left(m_{\text{sat}}-m_{\text{Sun}}\right)}{r^{\text{'}}}^2}{kF\left(\mathrm{\psi }\right)}$. (8)

Для случая отражения от диффузной сферы фазовая функция $F\left(\mathrm{\psi }\right)$ описывается как

$F\left(\mathrm{\psi }\right)=\frac{2}{3{\pi }^2}\left[(\pi -\mathrm{\psi })\mathrm{cos\psi }+\mathrm{sin\psi }\right]$, k = 0.1–0.18 (9)

Величину косинуса угла $\mathrm{\psi }$ в каждый момент наблюдения можно определить из численного моделирования движения

$\cos \left(\mathrm{\psi }\right)=\frac{{\mathbf{x}}_{\mathrm{E}}{\mathbf{x}}_{\mathrm{S}}}{{r^{\text{'}}}_{\mathrm{E}}{r^{\text{'}}}_{\mathrm{S}}}$, (10)

где ${\mathbf{x}}_{\mathrm{E}}$ и ${\mathbf{x}}_{\mathrm{S}}$ есть векторы положения наблюдателя и Солнца в системе координат, связанной со спутником (см. рис. 1).

$\begin{array}{l}{\mathbf{x}}_S=\left\{x_{\text{S}}-x_{\text{sat}};y_{\text{S}}-y_{\text{sat}};z_{\text{S}}-z_{\text{sat}}\right\},\\ {\mathbf{x}}_E=\left\{x_{\text{obs}}-x_{\text{sat}};y_{\text{obs}}-y_{\text{sat}};z_{\text{obs}}-z_{\text{sat}}\right\}.\end{array}$ (11)

 

Рис. 1. Конфигурация Солнце–спутник–наблюдатель.

 

По вычисленным мгновенным значениям А находится среднее значение величины A по всему используемому массиву наблюдений и вычисляется средняя оценка массы с использованием параметра A/m, найденного по позиционным наблюдениям.

ОТБРАКОВКА ИЗМЕРЕНИЙ. ОБОСНОВАНИЕ ВЫБРАННЫХ КРИТЕРИЕВ

Отбраковка измерений – самый ответственный этап в процессе улучшения орбит небесных тел, и околоземных объектов в том числе. В связи с этим был проведен численный эксперимент по выбору и обоснованию критериев отбраковки, результаты которого представлены в данной работе. Для эксперимента были выбраны объекты КМ, имеющие большое число наблюдений. Данные о наблюдениях объектов приведены в табл. 2.

 

Таблица 2. Описание использованных наблюдений объектов

№ объекта	Количество наблюдений	Период наблюдений, сут.	Интервал наблюдений внутри суток, ч.
10294	4700	06–23.12.2019	14:28–03:30
90008	1051	13–28.09.2021	16:41–02:01
90073	1749	06–23.12.2019	14:45–03:12

 

Номера объектов даны в соответствии с нумерацией в динамической базе данных космических объектов Института прикладной математики им. М.В. Келдыша РАН.

Рассматривались два правила отбраковки наблюдений:

Правило 1: |O–С|>$N{\sigma }^{\text{'}\text{'}}$; Правило 2: |O–С| >${\overline{M}}^{\text{'}\text{'}}$,

где N – целое число, $\overline{M}$ – пороговое значение O–C, по которому все наблюдения, удовлетворяющие Правилу 2, удаляются из процесса улучшения.

Поскольку традиционное Правило 1 плохо поддается строгой ранжировке, исследования проводились с Правилом 2. Полученные результаты легко согласуются с Правилом 1 по данным, приведенным на графиках. Точность наблюдений не ниже 0.01″, поэтому начальное значение $\overline{M}$ выбиралось равным 0.5″, а конечное – 3.5″.

Покажем, прежде всего, что даже самая жесткая отбраковка не позволяет полностью устранить систематику из представления наблюдений (табл. 3 и рис. 2). В табл. 3 приведены данные улучшения орбиты объекта 10294 с отбраковками при различных значениях $\overline{M}$. Здесь и далее в табл. 4 и табл. 5 σ, ″ – среднеквадратическая ошибка, полученная в процессе улучшения орбиты, “Σ ошибок в параметрах, км” – сумма модулей ошибок в прямоугольных координатах объекта, полученные в процессе улучшения орбиты.

 

Таблица 3. Данные о процессе улучшения орбиты объекта № 10294, полученные при полном наборе наблюдений и с использованием отбраковок

Способ отбраковки	Кол-во наблюд.	σ, ″	Число Тодда	Σ ошибок в параметрах, км	Ошибки параметров, x…, км, v..., км/с, A/m, м2/кг
без отбраковки	4700	1.426	2.216×105	0.02135	x ± 5.90×10–3; y ± 7.30×10–3; z ± 8.14×10–3; vx± 3.27×10–7; vy± 4.19×10–7; vz± 3.25×10–7; A/m ± 7.14×10–5
$\overline{M}=3^{\text{'}\text{'}}$	3856	0.8217	9.560×105	0.01897	x ± 4.70×10–3; y ± 8.59×10–3; z ± 5.68×10–3; vx± 3.99×10–7; vy± 6.93×10–7; vz± 1.99×10–7; A/m ± 1.33×10–4
$\overline{M}=2^{\text{'}\text{'}}$	3582	0.6587	2.080×106	0.02046	x ± 4.91×10–3; y ± 1.03×10–2; z ± 5.21×10–3; vx± 4.90×10–7; vy± 8.60×10–7; vz± 1.68×10–7; A/m ± 1.30×10–4
$\overline{M}=1^{\text{'}\text{'}}$	2666	0.4753	5.040×106	0.02271	x ± 5.79×10–3; y ± 1.20×10–2; z ± 4.88×10–3; vx± 6.11×10–7; vy± 1.08×10–6; vz± 1.49×10–7; A/m ± 2.87×10–4
$\overline{M}={0.5}^{\text{'}\text{'}}$	1365	0.2446	3.628×107	0.03178	x ± 6.61×10–3; y ± 1.94×10–2; z ± 5.70×10–3; vx ± 9.95×10–7; vy ± 1.60×10–6; vz ± 1.65×10–7; A/m ± 5.50×10–4

 

Рис. 2. Примеры представления наблюдений объекта № 10294 при различных уровнях их отбраковки: (а) – без отбраковки; (б) – $\overline{\mathit{M}}\mathit{=}{\mathbf{2}}^{\mathbf{\text{'}\text{'}}}$; (в) – $\overline{\mathit{M}}\mathbf{=}{\mathbf{1}}^{\mathbf{\text{'}\text{'}}}$; (г) – $\overline{\mathit{M}}\mathbf{=}\mathbf{0}\mathbf{.}{\mathbf{5}}^{\mathbf{\text{'}\text{'}}}$.

 

На рис. 2–4 на каждой панели сверху вниз представлены графики блеска и невязок для соответствующих отбраковок, перечисленных в табл. 3–5.

Далее аналогичные сведения приводятся по объекту № 90008 (табл. 4, рис. 3) и объекту 90073 (табл. 5, рис. 4).

 

Таблица 4. Данные о процессе улучшения орбиты объекта № 90008, полученные при полном наборе наблюдений и с использованием отбраковок

Способ отбраковки	Кол-во наблюд.	σ, ″	Число Тодда	Σ ошибок в параметрах, км	Ошибки параметров, x…, км; v..., км/с; A/m, м2/кг
без отбраковки	1051	1.308	2.989×105	0.03220	x ± 1.08×10–2; y ± 1.28×10–2; z ± 8.62×10–3; vx ± 1.42×10–6; vy ± 5.42×10–7; vz ± 1.00×10–6; A/m ± 3.39×10–4
$\overline{M}=3^{\text{'}\text{'}}$	933	0.8888	4.332×105	0.02375	x ± 7.89×10–3; y ± 9.90×10–3; z ± 5.95×10–3; vx ± 1.28×10–6; vy ± 5.78×10–7; vz ± 8.46×10–7; A/m ± 2.47×10–4
$\overline{M}=2^{\text{'}\text{'}}$	799	0.4402	1.641×106	0.01269	x ± 4.20×10–3; y ± 5.54×10–3; z ± 2.94×10–3; vx ± 1.24×10–6; vy ± 1.23×10–6; vz ± 7.07×10–7; A/m ± 1.39×10–4
$\overline{M}=1^{\text{'}\text{'}}$	746	0.3452	4.068×106	0.01143	x ± 4.12×10–3; y ± 4.94×10–3; z ± 2.36×10–3; vx ± 1.51×10–6; vy ± 1.63×10–6; vz ± 6.54×10–7; A/m ± 1.32×10–4
$\overline{M}={0.5}^{\text{'}\text{'}}$	319	0.2445	8.443×106	0.03107	x ± 1.23×10–2; y ± 1.57×10–2; z ± 3.06×10–3; vx ± 2.53×10–6; vy ± 8.42×10–7; vz ± 5.68×10–7; A/m ± 1.84×10–4

 

Рис. 3. Примеры представления наблюдений объекта № 90008 при различных уровнях их отбраковки: (а) – без отбраковки; (б) – $\overline{\mathit{M}}\mathit{=}{\mathbf{2}}^{\mathbf{\text{'}\text{'}}}$; (в) – $\overline{\mathit{M}}\mathbf{=}{\mathbf{1}}^{\mathbf{\text{'}\text{'}}}$; (г) – $\overline{\mathit{M}}\mathbf{=}\mathbf{0}\mathbf{.}{\mathbf{5}}^{\mathbf{\text{'}\text{'}}}$.

 

Таблица 5. Данные о процессе улучшения орбиты объекта № 90073, полученные при полном наборе наблюдений и с использованием отбраковок

Способ отбраковки	Кол-во наблюд.	σ, ″	Число Тодда	Σ ошибок в параметрах, км	Ошибки параметров, x…, км, v..., км/с, A/m, м2/кг
без отбраковки	1749	0.8391	7.590×105	0.01910	x ± 6.36×10–3; y ± 8.06×10–3; z ± 4.67×10–3; vx ± 4.60×10–7; vy ± 6.48×10–7; vz ± 4.64×10–7; A/m ± 9.62×10–5
$\overline{M}=3^{\text{'}\text{'}}$	1711	0.7449	8.719×105	0.01752	x ± 5.78×10–3; y ± 7.53×10–3; z ± 4.20×10–3; vx ± 4.42×10–7; vy ± 6.16×10–7; vz ± 4.22×10–7; A/m ± 8.90×10–5
$\overline{M}=2^{\text{'}\text{'}}$	1623	0.6250	1.209×106	0.01625	x ± 5.34×10–3; y ± 7.28×10–3; z ± 3.61×10–3; vx ± 4.30×10–7; vy ± 6.09×10–7; vz ± 3.79×10–7; A/m ± 8.52×10–5
$\overline{M}=1^{\text{'}\text{'}}$	1314	0.3711	6.890×106	0.02030	x ± 7.15×10–3; y ± 1.05×10–2; z ± 2.64×10–3; vx ± 6.17×10–7; vy ± 8.78×10–7; vz ± 2.80×10–7; A/m ± 1.16×10–4
$\overline{M}=0.5^{\text{'}\text{'}}$	881	0.2433	9.270×106	0.02006	x ± 8.11×10–3; y ± 9.66×10–3; z ± 2.27×10–3; vx ± 6.33×10–7; vy ± 8.30×10–7; vz ± 2.36×10–7; A/m ± 1.19×10–4

 

Рис. 4. Примеры представления наблюдений объекта № 90073 при различных уровнях их отбраковки: (а) – без отбраковки; (б) – $\overline{\mathit{M}}\mathit{=}{\mathbf{2}}^{\mathbf{\text{'}\text{'}}}$; (в) – $\overline{\mathit{M}}\mathbf{=}{\mathbf{1}}^{\mathbf{\text{'}\text{'}}}$; (г) – $\overline{\mathit{M}}\mathbf{=}\mathbf{0}\mathbf{.}{\mathbf{5}}^{\mathbf{\text{'}\text{'}}}$.

 

Далее следуют графики для объекта № 900073.

Как видно из табл. 3–5, при уменьшении порогового значения среднеквадратичная ошибка уменьшается. При этом обусловленность задачи может вести себя по-разному, а соответственно, и ошибки параметров: они ухудшаются, если ухудшается обусловленность задачи и уменьшаются, если порядок числа Тодда сохраняется. Если обратить внимание на результаты на рис. 2–4, можно заметить, что при усилении критерия отбраковки в представлении наблюдений остаются систематические ошибки. В связи с этим поиск оптимального критерия отбраковки проводился по точности результатов прогноза движения объектов, полученных при различных вариантах отбраковки.

Данные результатов улучшения и прогноз по всем трем объектам приведены на рис. 5. Сплошной темно-серой линией обозначены величины, полученные без отбраковки, а темно-серой точкой – вариант отбраковки, выбранный алгоритмом в качестве оптимального. На графиках сверху вниз:

• количество используемых наблюдений (после отбраковки);
• число обусловленности Тодда (логарифмическая шкала);
• сумма модулей ошибок определяемых координат объекта, обозначенная как “ошибка ${\sigma }_{\mathrm{x}}+{\sigma }_{\mathrm{y}}+{\sigma }_z,$ км”
• среднеквадратичная ошибка прогноза на 9-месячном интервале (логарифмическая шкала для объекта № 10294 и линейная для остальных объектов);
• среднеквадратичная ошибка на интервале улучшения орбиты.

 

Рис. 5. Данные для поиска оптимального критерия отбраковки на примере объектов с номерами: (а) – 10294; (б) – 90008; (в) – 90073.

 

Наибольший интерес представляют графики среднеквадратичной ошибки прогноза и суммы модулей среднеквадратичных ошибок определяемых координат. На всех графиках на рис. 5 видна зависимость между этими величинами. Если сумма ошибок координат меньше начального значения, то прогноз точнее. Особенно это наблюдается для случая с объектом № 10294, где σ прогноза может достигать наилучшего значения, равного 2.7″, и наихудшего – 625″.

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

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

Таким образом, в качестве оптимального значения числа $\overline{M}$ следует выбирать то, которое дает уменьшение среднеквадратичных ошибок определяемых параметров, при сохранении порядка числа Тодда, характеризующего обусловленность задачи улучшения орбиты.

ПРОГРАММНАЯ РЕАЛИЗАЦИЯ АЛГОРИТМА ОТБРАКОВКИ

В программном комплексе “Численная модель движения систем ИСЗ” реализованы два правила, с помощью которых пользователь может выбрать, каким способом будут удалены грубые наблюдения:

• Правило 1: |O–С| >$N{\sigma }^{\text{'}\text{'}}$.
• Правило 2: |O–С| >${\overline{M}}^{\text{'}\text{'}}$.

Пользователь выбирает интервал изменения параметров N, $\overline{M}$ и соответствующий для каждого из них шаг. Улучшение орбиты объекта выполняется при каждом варианте отбраковки итеративно до тех пор, пока не будет достигнута сходимость на первой итерации процедуры Гаусса–Ньютона (3). После выполнения всех вычислений результаты для каждой отбраковки сохраняются, и программа производит выбор оптимального варианта по следующим требованиям:

• выбираются все результаты со среднеквадратичной ошибкой, меньшей, чем в начальном варианте, полученном при полной выборке наблюдений;
• среди выбранных вариантов фиксируются те, у которых порядок числа обусловленности Тодда, либо совпадает с начальным вариантом, либо не превосходит его более чем на один порядок;
• из оставшихся вариантов отдается предпочтение тому, у которого сумма среднеквадратических ошибок координат объекта наименьшая.

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

НАИЛУЧШИЕ СИСТЕМЫ ПАРАМЕТРОВ

По всем группам наблюдений объектов (табл. 6) были получены системы параметров с помощью описанной выше программы.

 

Таблица 6. Описание использованных наблюдений 2022 г.

№ объекта	Кол-во наблюдений	Период наблюдений, сут.
10294	1031	05.01.2022 – 28.01.2022
	195	17.06.2022 – 02.07.2022
	349	25.10.2022 – 08.11.2022
90008	349	05.01.2022 – 17.01.2022
	219	17.06.2022 – 02.07.2022
	133	28.10.2022 – 08.11.2022
90073	288	05.01.2022 – 28.01.2022
	201	20.05.2022 – 02.07.2022
	135	01.11.2022 – 08.11.2022

 

Из систем параметров объектов, найденных по наблюдениям 2022 г., по принципам, изложенным выше, для каждого объекта были выбраны наиболее точные системы, приведенные ниже (табл. 7–9).

 

Таблица 7. Наилучшая система параметров 2022 г. для объекта № 10294

Январь 2022 г.
Кол-во наблюдений	A/m, м2/кг ± ср. кв. ошибка	Ср. кв. ошибка, ″	Эпоха
1005	0.01689 ± 0.000188	0.7323	17.01.2022 00:00:00.000 UT
Координаты и скорости		Ср. квадратичные ошибки
x, км	–40809.665401 ±	Δx, м	0.013099
y, км	–8365.372883	Δy, м	0.011203
z, км	–6643.812779	Δz, м	0.007419
vx, км/с	0.659010	Δvx, м/с	9.908506e-7
vy, км/с	–2.986419	Δvy, м/с	7.972262e-7
vz, км/с	–0.302336	Δvz, м/с	5.913179e-7
Элементы орбиты		Ср. квадратичные ошибки
a	42164.465418, км	Δa, м	0.013534
e	0.000890	Δe	1.698625e-7
i	10°41΄52.83″	Δi, ʺ	0.000012
Ω	313°59΄50.81″	ΔΩ, ʺ	0.000583
ω	115°23΄5.94″	Δω, ʺ	0.002668
М	122°34΄7.82″	ΔМ, ʺ	0.002583

 

Таблица 8. Наилучшая система параметров 2022 г. для объекта № 90008

Январь 2022 г.
Кол-во наблюдений	A/m, м2/кг ± ср. кв. ошибка	Ср. кв. ошибка, ″	Эпоха
244	0.00323 ± 0.000378	0.2265	09.01.2022 00:00:00.000 UT
Координаты и скорости		Ср. квадратичные ошибки
x, км	–39282.776554	Δx, м	0.006595
y, км	14388.501332	Δy, м	0.007737
z, км	–4004.365510	Δz, м	0.004707
vx, км/с	–1.010728	Δvx, м/с	4.759492e-7
vy, км/с	–2.888130	Δvy, м/с	7.737156e-7
vz, км/с	–0.389199	Δvz, м/с	2.540865e-7
Элементы орбиты		Ср. квадратичные ошибки
a	42158.507787, км	Δa, м	0.007078
e	0.003868	Δe	1.821231e-7
i	9°6΄26.46″	Δi, ʺ	0.000011
Ω	303°13΄15.27″	ΔΩ, ʺ	0.000012
ω	252°54΄52.43″	Δω, ʺ	0.000025
М	324°21΄21.18″	ΔМ, ʺ	0.000086

 

Таблица 9. Наилучшая система параметров 2022 г. для объекта № 90073

Январь 2022 год
Кол-во наблюдений	A/m, м2/кг ± ср. кв. ошибка	Ср. кв. ошибка, ″	Эпоха
135	0.08642 ± 0.002027	0.2897	09.01.2022 00:00:00.000 UT
Координаты и скорости		Ср. квадратичные ошибки
x, км	–41881.321996	Δx, м	0.019224
y, км	2632.634711	Δy, м	0.011702
z, км	5207.127258	Δz, м	0.010298
vx, км/с	–0.237058	Δvx, м/с	2.835499e-6
vy, км/с	–3.043755	Δvy, м/с	1.642194e-6
vz, км/с	–0.283246	Δvz, м/с	9.393271e-7
Элементы орбиты		Ср. квадратичные ошибки
a	42171.850937, км	Δa, м	0.019357
e	0.004340	Δe	1.7357999e-7
i	8°52΄11.95″	Δi, ʺ	0.000018
Ω	49°4΄10.77″	ΔΩ, ʺ	0.000052
ω	358°17΄50.28″	Δω, ʺ	0.000615
М	128°18΄55.53″	ΔМ, ʺ	0.000453

 

Для подтверждения правильности выбора наилучшей системы параметров 2022 г. для каждого объекта по всем системам параметров было проведено представление наблюдений 2023 г. (табл. 10–12).

 

Таблица 10. Данные о точности прогноза движения объекта № 10294

O–C	От января 2022 г.	От июня 2022 г.	От ноября 2022 г.
	Прогноз на июль 2023 г.
Δαʺ	56.7	102.4	85.5
Δδʺ	8.9	12.4	12.0
	Прогноз на ноябрь 2023 г.
Δαʺ	7.6	105.8	87.9
Δδʺ	2.3	16.0	14.8

 

Таблица 11. Данные о точности прогноза движения объекта № 90008

O–C	От января 2022 г.	От июня 2022 г.	От ноября 2022 г.
	Прогноз на июль 2023 г.
Δαʺ	142.3	120.0	118.5
Δδʺ	12.7	9.7	10.5
	Прогноз на ноябрь 2023 г.
Δαʺ	10.8	189.1	122.9
Δδʺ	3.5	19.4	14.5

 

Таблица 12. Данные о точности прогноза движения объекта № 90073

O–C	От января 2022 г.	От июня 2022 г.	От ноября 2022 г.
	Прогноз на июль 2023 г.
Δαʺ	1811.0	2762.5	1972.0
Δδʺ	271.8	412.0	293.4
	Прогноз на ноябрь 2023 г.
Δαʺ	791.1	1666.8	2577.4
Δδʺ	30.0	71.0	118.8

 

Таким образом, для всех трех объектов наилучшими следует считать системы параметров, полученные по наблюдения января 2022 г.

ИСПОЛЬЗОВАНИЕ ДАННЫХ ФОТОМЕТРИИ

Рассмотрим объект № 10294 по базе данных космических объектов Института прикладной математики им. М.В. Келдыша РАН. Этот объект является отработавшим итальянским спутником Sirio. Масса и парусность, которые спутник имел при запуске, известны: m = 218 кг, A/m = 0.009 м²/кг.

В табл. 13 представлены минимальное А_мин, максимальное А_макс и среднее А_ср значения площади миделевого сечения, параметр парусности объекта A/m, полученные по позиционным наблюдениям, среднеквадратическая ошибка $\sigma$ в конце итерационного процесса улучшения орбиты объекта и оценка его массы m.

 

Таблица 13. Определение параметров объекта № 10294

N	Период наблюдений	Амин, м2	Амакс, м2	Аср, м2	A/m, м2/кг	σ, ″	m, кг
1	03.03.2019 – 09.03.2019	0.2862	16.5093	0.8771	0.0142	0.3868	61.7750
2	12.09.2019 – 23.09.2019	0.6413	44.7173	2.5810	0.0093	0.2035	278.4369
3	06.12.2019 – 23.12.2019	0.4141	278.2284	3.7776	0.0158	0.9289	239.1645
4	22.06.2020 – 28.06.2020	0.5137	51.8611	2.1053	0.0110	0.2237	192.1120
5	10.09.2020 – 29.09.2020	0.5062	47.8325	2.5714	0.0104	0.4582	248.1700
6	12.09.2021 – 27.09.2021	0.6961	16.2754	1.8211	0.0120	0.3713	151.9454
7	11.10.2021 – 29.10.2021	0.7227	21.5851	1.8104	0.0176	0.4650	103.1416
8	03.11.2021 – 06.11.2021	1.0605	31.5174	1.9990	0.0261	0.5967	76.5682
9	17.06.2022 – 02.07.2022	0.4876	10.1689	1.6471	0.0135	0.3227	121.7657
10	25.10.2022 – 08.11.2022	0.6832	21.1805	1.5161	0.0190	0.3460	79.9016
11	15.09.2023 – 22.09.2023	0.3617	37.6496	2.6336	0.0168	0.8307	156.6146

 

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

Кроме того, на примере объекта № 10294 было рассмотрено влияние отбраковки наблюдений на процесс совместного использования позиционных и фотометрических измерений для определения динамических параметров объекта (табл. 14).

 

Таблица 14. Совместное использование позиционных и фотометрических измерений. Влияние отбраковки

N	Период наблюдений	$1^{\text{'}\text{'}}$		$2^{\text{'}\text{'}}$		$3^{\text{'}\text{'}}$
		A/m	m, кг	A/m	m, кг	A/m	m, кг
1	03.03.2019 – 09.03.2019	0.015	59	0.015	62	0.015	62
2	12.09.2019 – 23.09.2019	0.009	270	0.010	246	0.010	240
3	06.12.2019 – 23.12.2019	0.018	154	0.018	185	0.017	198
4	22.06.2020 – 28.06.2020	0.011	188	0.011	187	0.011	187
5	10.09.2020 – 29.09.2020	0.010	256	0.011	240	0.010	246
6	12.09.2021 – 27.09.2021	0.012	149	0.012	150	0.015	121
7	11.10.2021 – 29.10.2021	0.018	93	0.019	93	0.017	109
8	03.11.2021 – 06.11.2021	0.080	22	0.026	76	0.026	78
9	05.01.2022 – 28.01.2022	0.011	227	0.017	117	0.020	99
10	17.06.2022 – 02.07.2022	0.013	120	0.013	122	0.013	122
11	25.10.2022 – 08.11.2022	0.018	76	0.005	331	0.002	734
12	15.09.2023 – 22.09.2023	0.016	185	0.017	157	0.016	162
Значения при запуске		0.009	218	0.009	218	0.009	218
Средние		0.019	150	0.014	164	0.014	196
Средние без 1, 7, 8, 11		0.013	194	0.014	176	0.014	172

 

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

Проанализируем еще один объект, под номером 90073, который представляет собой небольшой фрагмент китайского метеорологического спутника Fengyun 2D, запущенного в декабре 2006 г. В работе (Bakhtigaraev и др., 2018) на основании исследования обширного наблюдательного материала получено, что максимальное значение парусности $\gamma$ этого объекта равно 0.18 м²/кг.

Как видно из представленных в табл. 15 данных, оценка средней массы спутника по фотометрическим измерениям колеблется около значения m = 9.156 кг. Оценка парусности объекта по всем имеющимся в нашем распоряжении массивам наблюдений дает среднее значение A/m = 0.1169 м²/кг.

 

Таблица 15. Определение параметров объекта № 90073

N	Период наблюдений	Амин, м2	Амакс, м2	Аср, м2	A/m, м2/кг	σ, ″	m, кг
1	3.09.2019 – 24.09.2019	0.0309	22.6135	0.6379	0.0925	0.5641	6.8995
2	06.12.2019 – 24.12.2019	5.4540×10–5	15.4045	0.4335	0.1001	0.2433	4.3291
3	22.06.2020 – 29.06.2020	0.2002	5.4631	1.1560	0.1295	0.2963	8.9285
4	11.09.2020 – 24.09.2020	0.0376	7.2564	0.8634	0.1230	0.3455	7.0197
5	12.09.2021 – 27.09.2021	0.0750	11.3801	1.2901	0.1418	0.5589	9.0979
6	11.10.2021 – 30.10.2021	0.0692	73.5874	1.3580	0.1146	1.0958	11.8517
7	03.11.2021 – 06.11.2021	0.1032	16.9116	1.1579	0.0983	0.5124	11.7816
8	05.01.2022 – 29.01.2022	0.1289	9.3395	1.5014	0.0864	0.2897	17.3731
9	17.06.2022 – 03.07.2022	0.0403	8.5674	1.3625	0.0809	0.2381	16.8338
10	28.10.2022 – 30.10.2022	0.2010	9.9854	1.1843	0.1488	0.2421	7.9572
11	01.11.2022 – 08.11.2022	0.1419	5.3395	0.8486	0.1388	0.2535	6.1133
12	17.09.2023 – 22.09.2023	0.0601	0.7655	0.2488	0.1479	0.9458	1.6819

 

Аналогичные оценки были получены по всем совместным наблюдениям объекта № 90008 (табл. 16). Этот объект, по предположениям наблюдателей, также является частью китайского метеорологического спутника Fengyun 2D и имеет форму диска (поскольку является крышкой).

 

Таблица 16. Определение параметров объекта № 90008

N	Период наблюдений	Амин, м2	Амакс, м2	Аср, м2	A/m, м2/кг	σ, ″	m, кг
1	25.09.2011 – 07.10.2011	0.2665	0.7309	0.4566	0.0130	0.2024	35.2042
2	30.08.2014 – 06.09.2014	0.3245	1.1734	0.5887	0.0157	0.2744	37.4446
3	05.03.2017 – 16.03.2017	0.3210	20.7162	1.2663	0.0308	0.3681	41.1658
4	10.09.2020 – 21.09.2020	0.3019	28.2536	2.7553	0.0118	0.4053	233.8612
5	30.05.2021 – 07.06.2021	0.3334	8.3031	1.3558	0.0108	0.3805	125.9946
6	12.09.2021 – 27.09.2021	0.3368	53.4874	5.6492	0.0126	0.3452	447.4261
7	05.01.2022 – 17.01.2022	0.2671	8.4292	0.5709	0.0032	0.2265	176.9545
8	17.06.2022 – 02.07.2022	0.2475	9.8093	0.7515	0.0104	0.2738	72.1387
9	28.10.2022 – 08.11.2022	0.2218	6.6798	0.8371	0.0184	0.2302	45.4413

 

Как показывают данные, приведенные в табл. 16, разброс всех определенных значений параметров и, в том числе, площади миделевого сечения спутника, очень велик, а соответственно велик и разброс масс. Определенные по данным фотометрии изменения площади миделевого сечения (рис. 6) подтверждают тот факт, что объект плоский, причем он наблюдается как бы сбоку, но при этом освещено ребро и основная плоскость тоже частично освещена. За все сеансы наблюдения максимальная освещенная плоская часть поверхности наблюдалась один раз в сентябре 2021 г., а минимальная – дважды в сентябре 2011 г. и в августе–сентябре 2014 г.

 

Рис. 6. Значения фазового угла и площади миделевого сечения объекта № 90008, полученные по измерениям блеска для всех сеансов наблюдений.

 

ЗАКЛЮЧЕНИЕ

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

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

Работа выполнена с использованием оборудования Центра коллективного пользования “Терскольская обсерватория”.

Работа выполнена в рамках государственного задания Министерства науки и высшего образования Российской Федерации (тема № FSWM-2020-0049).

作者简介

Т. Бордовицына

Томский государственный университет

编辑信件的主要联系方式.
Email: bordovitsyna@mail.ru
俄罗斯联邦, Томск

Н. Бахтигараев

Институт астрономии РАН

Email: nail@inasan.ru
俄罗斯联邦, Москва

П. Левкина

Институт астрономии РАН

Email: ayvazovskaya@inasan.ru
俄罗斯联邦, Москва

Н. Попандопуло

Томский государственный университет

Email: nikas.popandopulos@gmail.com
俄罗斯联邦, Томск

К. Салейко

Томский государственный университет

Email: ksenia.saleiko@yandex.ru
俄罗斯联邦, Томск

И. Томилова

Томский государственный университет

Email: irisha_tom@mail.ru
俄罗斯联邦, Томск

О. Новикова

Томский государственный университет

Email: freedgreed21@gmail.com
俄罗斯联邦, Томск

参考

补充文件