Астрономический вестник, 2022, T. 56, № 3, стр. 206-216

Определение предварительной орбиты в некомпланарном случае

В. Б. Кузнецов *

Институт прикладной астрономии РАН
Санкт-Петербург, Россия

* E-mail: vb.kuznetsov@iaaras.ru

Поступила в редакцию 25.10.2021
После доработки 02.12.2021
Принята к публикации 16.12.2021

Полный текст (PDF)

Аннотация

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

Ключевые слова: задача двух тел, определение орбиты, метод Нелдера–Мида, комета Борисова

ВВЕДЕНИЕ

Проблема определения предварительной орбиты, в ее классической постановке ограниченной кеплеровской задачей двух тел, остается актуальной уже более трехсот лет. Методы, используемые для ее решения, вместе с ограничениями им присущими были описаны в предыдущей работе (Кузнецов, 2019). Там же, для определения предварительной орбиты по трем наблюдениям, в моменты времени t1, t2 и t3, был рассмотрен универсальный метод, основанный на трансцендентных уравнениях, который позволяет находить решение без каких‑либо ограничений на интервалы времени между наблюдениями и эксцентриситет искомой орбиты. Обобщим этот метод на случай n ≥ 3 наблюдений в моменты времени ti (i = 1, 2, …, n). Обозначим через ri – радиус-векторы искомой орбиты в соответствующие моменты времени. Тогда уравнения для гелиоцентрических радиус-векторов будут иметь следующий вид:

(1)
${{{\mathbf{r}}}_{i}} = {{{\mathbf{e}}}_{i}}{{\rho }_{i}} - {{{\mathbf{R}}}_{i}},$
где ei – единичные векторы наблюденного направления на объект, ρi – топоцентрические расстояния до объекта, Ri – векторы положения Солнца относительно топоцентра. Некомпланарным назовем случай, когда множества векторов ei и Ri не принадлежат одной плоскости – плоскости эклиптики (Кузнецов, 2019). Для него определение орбиты сводится к решению системы из n –1 уравнения относительно двух топоцентрических расстояний до объекта.

В силу того, что орбита в задаче двух тел плоская и векторы ri связаны условием компланарности, мы можем составить смешанные произведения, равные нулю, из любых их троек. Таким образом возможно выразить все топоцентрические расстояния через два выбранных (обычно это ρ1 и ρn).

Параметр орбиты p определяется из следующего соотношения (Stumpff, 1959):

(2)
$p = \frac{{ \pm {{r}_{1}}\left| {{{{\mathbf{r}}}_{{\left[ {(n + 1)/2} \right]}}} \times {{{\mathbf{r}}}_{n}}} \right| \mp {{r}_{{\left[ {(n + 1)/2} \right]}}}\left| {{{{\mathbf{r}}}_{1}} \times {{{\mathbf{r}}}_{n}}} \right| \pm {{r}_{n}}\left| {{{{\mathbf{r}}}_{1}} \times {{{\mathbf{r}}}_{{\left[ {(n + 1)/2} \right]}}}} \right|}}{{ \pm \left| {{{{\mathbf{r}}}_{{\left[ {(n + 1)/2} \right]}}} \times {{{\mathbf{r}}}_{n}}} \right| \mp \left| {{{{\mathbf{r}}}_{1}} \times {{{\mathbf{r}}}_{n}}} \right| \pm \left| {{{{\mathbf{r}}}_{1}} \times {{{\mathbf{r}}}_{{\left[ {(n + 1)/2} \right]}}}} \right|}} \geqslant 0,$
где ${{r}_{i}} = \left| {{{{\mathbf{r}}}_{i}}} \right|$, верхние и нижние знаки ± и $ \mp $ в (2) соответствуют положительным и отрицательным знакам синуса угла между векторами в последующем векторном произведении, а индекс [(n + 1)/2] представляет собой антье (целую часть) числа (n + 1)/2. Для (2) числитель и знаменатель не могут иметь разные знаки.

Для нахождения неизвестных топоцентрических расстояний мы можем построить систему из n – 1 трансцендентного уравнения Шефера для определения орбиты по векторам ri, ri +1 и интервалу времени [ti, ti +1] (Шефер, 2010; Кузнецов, 2019):

(3)
$\begin{gathered} {{f}_{i}} = \pm \left[ {\sqrt {{{r}_{i}}{{r}_{{i + 1}}} + {{{\mathbf{r}}}_{i}}{{{\mathbf{r}}}_{{i + 1}}}} \pm \frac{1}{{\sqrt 8 }}X\left( {{{x}_{{(i + 1)i}}}} \right)} \right. \times \\ \left. {\frac{{}}{{}} \times \,\,\left( {{{r}_{i}} + \sqrt {2\left( {{{r}_{i}}{{r}_{{i + 1}}} + {{{\mathbf{r}}}_{i}}{{{\mathbf{r}}}_{{i + 1}}}} \right)} \left( {2{{x}_{{(i + 1)i}}} - 1} \right) + {{r}_{{i + 1}}}} \right)} \right], \\ \sqrt {{{r}_{i}} + \sqrt {2\left( {{{r}_{i}}{{r}_{{i + 1}}} + {{{\mathbf{r}}}_{i}}{{{\mathbf{r}}}_{{i + 1}}}} \right)} \left( {2{{x}_{{(i + 1)i}}} - 1} \right) + {{r}_{{i + 1}}}} - \\ - \,\,k\left( {{{t}_{{i + 1}}} - {{t}_{i}}} \right) - kL\left( {{{\rho }_{i}} - {{\rho }_{{i + 1}}}} \right) = 0, \\ \end{gathered} $
где $i = \overline {1,n - 1} $, индекс (i + 1)i состоит из двух чисел i + 1 и i, k – постоянная Гаусса или ее аналоги, L = 0.00576832 сут/а. е. – постоянная аберрации (или ее аналоги), X(x(i+ 1)i) – гипергеометрическая функция: X(x) = 4/3 F(1, 3, 5/2; x):

(4)
$\begin{gathered} X\left( x \right) = \frac{1}{{1 - x}}\left( {1 + \frac{2}{3}\frac{1}{{\sqrt {1 - x} \left( {\sqrt {1 - x} + 1} \right)}}} \right. \times \\ \times \,\,\left. {\left( {1 - \frac{1}{5}\frac{{\sqrt {1 - x} - 1}}{{\sqrt {1 - x} + 1}}F\left( {\frac{1}{2},1,\frac{7}{2};\frac{{\sqrt {1 - x} - 1}}{{\sqrt {1 - x} + 1}}} \right)} \right)} \right). \\ \end{gathered} $

Переменные x(i+ 1)i выражаются через (1) и (2) следующим образом:

(5)
${{x}_{{(i + 1){\kern 1pt} i}}} = \frac{1}{2} \pm \frac{{{{r}_{i}}{{r}_{{i + 1}}} - {{{\mathbf{r}}}_{i}}{{{\mathbf{r}}}_{{i + 1}}} - p\left( {{{r}_{i}} + {{r}_{{i + 1}}}} \right)}}{{2p\sqrt {2\left( {{{r}_{i}}{{r}_{{i + 1}}} + {{{\mathbf{r}}}_{i}}{{{\mathbf{r}}}_{{i + 1}}}} \right)} }},$
где ±, как и ранее, соответствует положительным и отрицательным знакам синуса угла между векторами. Таким образом, задачу определения некомпланарной орбиты можно свести к решению системы (3) относительно двух топоцентрических расстояний.

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

БЕЗРАЗМЕРНЫЕ ПЕРЕМЕННЫЕ НА КРУГЕ И КВАДРАТЕ

В работе (Кузнецов, 2019) были показаны недостатки использования топоцентрических расстояний в качестве переменных: область поиска решений в некомпланарном случае представляет собой открытый квадрант. Для того, чтобы этого избежать, предлагалось использовать нормированные безразмерные переменные, составляющие единичный вектор N = (Nx, Ny, Nz), представляющий собой вектор нормали к плоскости искомой орбиты. Топоцентрические расстояния ρi (i = 1, 2, …, n) выражаются через N следующим образом (Sarnecki, 1997):

(6)
${{\rho }_{i}} = \frac{{{\mathbf{N}}{{{\mathbf{R}}}_{i}}}}{{{\mathbf{N}}{{{\mathbf{e}}}_{i}}}} \geqslant 0.$

Условие единичности N может быть определено так:

(7)
$N_{x}^{2} + N_{y}^{2} + N_{z}^{2} = 1.$

После подстановки (6) во все уравнения, содержащие ρi, получаем аналог (3) для fi (i = 1, 2, …, n – 1) с ограничением (7) и безразмерными неизвестными {Nx, Ny, Nz}:

(8)
${{f}_{i}}({\mathbf{N}}) = 0.$

Для уравнений (8) область решения ограничивается поверхностью единичной сферы. Точки (при n = 3) или области (при n > 3) пересечения fi являются возможными направлениями N. Так как кривые (3) симметричны относительно плоскости эклиптики, то можно выбрать одну из полусфер с прямым или обратным орбитальным движением. Тогда условие (7) можно выразить как:

(9)
${{N}_{z}} = \sqrt {1 - N_{x}^{2} - N_{y}^{2}} $
и область решения единичной сферы отобразить на единичную область круга эклиптики. Следовательно, мы возвращаемся к двумерному случаю с неизвестными Nx и Ny.

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

(10)
$\left. \begin{gathered} {{N}_{x}} = {{N}_{{xs}}}\sqrt {1 - 0.5N_{{ys}}^{2}} , \hfill \\ {{N}_{y}} = {{N}_{{ys}}}\sqrt {1 - 0.5N_{{xs}}^{2}} , \hfill \\ \end{gathered} \right\}$
где Nxs и Nys – координаты, соответствующие квадрату.

ПОРЯДОК ПОЛОЖЕНИЙ ОБЪЕКТА НА ОРБИТЕ

При решении системы (8), как правило, не учитывается порядок точек пересечения лучей зрения (наблюдений) с орбитой. Поэтому возможно появление решений, для которых порядок соответствующих точек орбиты не соответствует хронологическому. Для их контроля можно воспользоваться следующими условиями (для n наблюдений):

(11)
$0.99999 \leqslant \left| {\frac{{\Delta {{\theta }_{{n1}}}}}{{\sum\limits_{l = 1}^{n - 1} {\Delta {{\theta }_{{(l + 1)l}}}} }}} \right| \leqslant 1.00001,$
где Δθij величины углов между ri и rj.

ОСОБЫЕ ТОЧКИ В ОБЛАСТИ ПОИСКА

Для безразмерных переменных {Nxs, Nys} существуют особые точки, в которых ρi (i = 1, 2, …, n) из (6) становятся неопределенными. В зависимости от того, обращается ли числитель или знаменатель в ноль, значение соответствующего топоцентрического расстояния может стать либо бесконечно малым, либо бесконечно большим. К тому же положительность расстояний требует одного знака для числителя и знаменателя. 2n уравнений, соответствующих равным нулю числителям и знаменателям (6), задают сетку из 2n кривых, попарное пересечение которых задает границы областей, в которых может быть определено искомое решение. Разность координат точек пересечения этих кривых позволяет определить наименьший масштаб триангуляции. Уравнения в “квадратных” координатах будут иметь следующий вид (i, j = 1, 2, …, n):

(12)
$\left. \begin{gathered} {{f}_{{{\text{числ}}\_i}}}\, = \,{{X}_{i}}{{N}_{{xs}}}\sqrt {1\, - \,0.5N_{{ys}}^{2}} \, + \,{{Y}_{i}}{{N}_{{ys}}}\sqrt {1\, - \,0.5N_{{xs}}^{2}} \, + \\ + \,{{Z}_{i}}\sqrt {1\, - \,N_{{xs}}^{2}\, - \,N_{{ys}}^{2}\, + \,N_{{xs}}^{2}N_{{ys}}^{2}} \, = \,0, \\ {{f}_{{{\text{знам}}\_j}}}\, = \,{{e}_{{jx}}}{{N}_{{xs}}}\sqrt {1\, - \,0.5N_{{ys}}^{2}} \, + \,{{e}_{{jy}}}{{N}_{{ys}}}\sqrt {1\, - \,0.5N_{{xs}}^{2}} \, + \\ + \,{{e}_{{jz}}}\sqrt {1\, - \,N_{{xs}}^{2}\, - \,N_{{ys}}^{2}\, + \,N_{{xs}}^{2}N_{{ys}}^{2}} \, = \,0, \\ \end{gathered} \right\}$
где Ri = (Xi, Yi, Zi) – эклиптические векторы положения Солнца относительно топоцентра, ej = (ejx, ejy, ejz) – единичные эклиптические векторы наблюденного направления на объект.

Система (12) нелинейная и для поиска ее решения удобно применить метод продолжения решения по параметру с наилучшей параметризацией (Шалашилин, Кузнецов, 1999). Задачу можно свести к нормальной системе из трех обыкновенных дифференциальных уравнений (ОДУ) в виде (i, j = 1, 2, …, n) (Кузнецов, 2021):

(13)
$\left. \begin{gathered} \frac{{d{{N}_{{xs}}}}}{{d\tilde {s}}} = {{\left( {{{{\mathbf{R}}}_{i}} \times {{{\mathbf{e}}}_{j}}} \right)}_{x}}\sqrt {1 - 0.5N_{{xs}}^{2}} + \\ + \,\,\frac{{0.5{{{\left( {{{{\mathbf{R}}}_{i}} \times {{{\mathbf{e}}}_{j}}} \right)}}_{y}}{{N}_{{xs}}}{{N}_{{ys}}}}}{{\sqrt {1 - 0.5N_{{ys}}^{2}} }}, \\ \frac{{d{{N}_{{ys}}}}}{{d\tilde {s}}} = {{\left( {{{{\mathbf{R}}}_{i}} \times {{{\mathbf{e}}}_{j}}} \right)}_{y}}\sqrt {1 - 0.5N_{{ys}}^{2}} + \\ + \,\,\frac{{0.5{{{\left( {{{{\mathbf{R}}}_{i}} \times {{{\mathbf{e}}}_{j}}} \right)}}_{x}}{{N}_{{xs}}}{{N}_{{ys}}}}}{{\sqrt {1 - 0.5N_{{xs}}^{2}} }}, \\ \frac{{d\mu }}{{d\tilde {s}}} = \\ = \frac{{\left( {1 - 0.5\left( {N_{{xs}}^{2} + N_{{ys}}^{2}} \right)} \right)\left[ {\left( {{{{\mathbf{R}}}_{i}} \times {{{\mathbf{e}}}_{j}}} \right){\mathbf{N}}} \right]}}{{\sqrt {\left( {1 - 0.5N_{{xs}}^{2}} \right)\left( {1 - 0.5N_{{ys}}^{2}} \right)\left( {N_{{xs}}^{2} - 1} \right)\left( {N_{{ys}}^{2} - 1} \right)} }}, \\ \end{gathered} \right\}$
где µ ∈ [0, 1] – параметр гомотопии, $\tilde {s}$ – длина дуги кривой множества решений. В третьем уравнении в квадратных скобках представлено смешанное произведение векторов в эклиптической системе координат, а в первых двух уравнениях нижние индексы у векторного произведения обозначают соответствующие компоненты. Следует интегрировать систему (13) в направлении роста параметра $\tilde {s}$, пока не будет достигнуто значение µ = 0. Решение систем (13) даст нам n2 точек, минимальные и максимальные координаты которых позволят ограничить на нашем квадрате прямоугольную область с наиболее тонкой структурой области решений и, возможно, наиболее плотным расположением искомых решений системы (8).

ТРИАНГУЛЯЦИЯ ОБЛАСТИ ПОИСКА РЕШЕНИЯ

Решение системы (8) предлагается искать в виде минимума целевой функции

(14)
${{f}_{{goal}}}\left( {\mathbf{N}} \right) = \sum\limits_{i = 1}^{n - 1} {f_{i}^{2}\left( {\mathbf{N}} \right)} ,$
где n – 1 – число уравнений системы (8), а n будет, соответственно, равно числу наблюдений. Причем поиск основывается на использовании триангуляции, т.е. на разбиении области возможных решений на непересекающиеся треугольники (Самотохин, Хуторовский, 2014). Эти треугольники необходимы как начальное приближение для переменных симплексов – подвижных треугольников, по которым выполняется поиск минимума целевой функции методом Нелдера–Мида (Химмельблау, 1975).

Вся максимально возможная область решений представляет собой квадрат с полустороной единичной длины. В качестве базовой триангуляции разобьем этот квадрат на 1024 “малых квадрата”, каждый с длиной стороны, равной 0.0625. Затем каждый из этих “малых квадратов” разобьем на четыре равных треугольника с общей вершиной в центре (Кузнецов, 2021). Этими равнобедренными треугольниками, с длиной основания 0.0625, мы заполняем всю область, для которой хотя бы одна из вершин треугольника удовлетворяет условию положительности n уравнений (6). Исключением является “область особых точек” (описанная в предыдущем разделе). Она увеличивается во все стороны до совпадения с границами ближайшего “малого квадрата” и получает вид прямоугольника с длинами сторон кратными 0.0625. Далее мы разбиваем эту область на квадраты в два раза меньшего размера, т.е. с длиной стороны 0.03125. Затем мы разбиваем полученные “малые квадраты” на треугольники с общей вершиной в центре. Для них также выполняется проверка того, что хотя бы для одной из вершин все ρi > 0 (i = 1, 2, …, n) и верно условие (11).

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

Если левая вершина не удовлетворяет условиям, то мы переносим ее на середину основания (рис. 1). Так же поступаем и с правой вершиной, если для нее не выполнены условия.

Рис. 1.

Перенос левой вершины треугольника.

В случае средней вершины опускаем ее на половину высоты (рис. 2).

Рис. 2.

Перенос средней вершины треугольника.

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

Рис. 3.

Перенос левой и средней вершин треугольника.

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

Рис. 4.

Перенос левой и правой вершин треугольника.

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

Таким образом мы заполняем всю область возможных решений треугольниками указанных видов. Они упорядочены в последовательность по мере построчного заполнения с левого нижнего угла “большого квадрата” до правого верхнего угла. За “большими треугольниками” в последовательности по тому же правилу идут “малые”.

ОПРЕДЕЛЕНИЕ РАНГОВ ТРЕУГОЛЬНИКОВ

Полученное множество из нескольких сотен треугольников, с точки зрения перспективы поиска решений, можно разделить на несколько классов. В работе (Самотохин, Хуторовский, 2014) вводится понятие ранга треугольника для случая n = 3, т.е. двух уравнений системы (8). Обобщим их подход для n > 3. Предлагается выполнить линейную интерполяцию возможного решения по значениям функций fi и fj (i, j = 1, 2,…, n–1) в вершинах треугольника, обозначим их как fi(m) (m = 1, 2, 3). Тогда координаты точки линейной интерполяции ($N_{{xs(0)}}^{{ij}},N_{{ys(0)}}^{{ij}}$) можно выразить так:

(15)
$\left. \begin{gathered} N_{{xs(0)}}^{{ij}} = \\ = {{N}_{{xs(1)}}}\, + \,{{\xi }_{{ij}}}\left( {{{N}_{{xs(2)}}}\, - \,{{N}_{{xs(1)}}}} \right)\, + \,{{\zeta }_{{ij}}}\left( {{{N}_{{xs(3)}}}\, - \,{{N}_{{xs(1)}}}} \right), \\ N_{{ys(0)}}^{{ij}} = \\ = {{N}_{{ys(1)}}}\, + \,{{\xi }_{{ij}}}\left( {{{N}_{{ys(2)}}}\, - \,{{N}_{{ys(1)}}}} \right)\, + \,{{\zeta }_{{ij}}}\left( {{{N}_{{ys(3)}}}\, - \,{{N}_{{ys(1)}}}} \right), \\ \end{gathered} \right\}$
где (Nxs(m), Nys(m)) (m = 1, 2, 3) – координаты соответствующих вершин,

(16)
$\left. \begin{gathered} {{\xi }_{{ij}}} = \frac{{{{f}_{{j(1)}}}{{f}_{{i(3)}}} - {{f}_{{i(1)}}}{{f}_{{j(3)}}}}}{{{{f}_{{i(2)}}}{{f}_{{j(3)}}} - {{f}_{{i(1)}}}{{f}_{{j(3)}}} - {{f}_{{i(2)}}}{{f}_{{j(1)}}} - {{f}_{{i(3)}}}{{f}_{{j(2)}}} + {{f}_{{i(3)}}}{{f}_{{j(1)}}} + {{f}_{{i(1)}}}{{f}_{{j(2)}}}}}, \hfill \\ {{\zeta }_{{ij}}} = \frac{{{{f}_{{i(1)}}}{{f}_{{j(2)}}} - {{f}_{{j(1)}}}{{f}_{{i(2)}}}}}{{{{f}_{{i(2)}}}{{f}_{{j(3)}}} - {{f}_{{i(1)}}}{{f}_{{j(3)}}} - {{f}_{{i(2)}}}{{f}_{{j(1)}}} - {{f}_{{i(3)}}}{{f}_{{j(2)}}} + {{f}_{{i(3)}}}{{f}_{{j(1)}}} + {{f}_{{i(1)}}}{{f}_{{j(2)}}}}}. \hfill \\ \end{gathered} \right\}$

Точка ($N_{{xs(0)}}^{{ij}},N_{{ys(0)}}^{{ij}}$) будет внутри треугольника при выполнении следующих условий:

(17)
${{\xi }_{{ij}}} > 0,\,\,\,\,{{\zeta }_{{ij}}} > 0,\,\,\,\,{{\xi }_{{ij}}} + {{\zeta }_{{ij}}} < 1.$

Указанная интерполяция является обобщением метода хорд для решения нелинейного уравнения от одной переменной. При этом необходимо учесть, что подобная линеаризация при выполнении (17) не гарантирует существование реального решения внутри треугольника. Здесь условия fi = 0 и fj = 0 (i, j = 1, 2, …, n–1) соответствуют двум прямым линиям – линиям нулевого уровня, пересекающимся в точке ($N_{{xs(0)}}^{{ij}},N_{{ys(0)}}^{{ij}}$).

Для каждой пары функций {fi, fj} (i, j = 1, 2, …, n – 1) опишем возможные подранги треугольников (Самотохин, Хуторовский, 2014):

• 0 – ни одна из линий нулевого уровня не пересекает ни одной из сторон треугольника (значения fi и fj для всех вершин имеют один знак);

• 1 – только одна из линий нулевого уровня имеет пересечения со сторонами треугольника (значения либо fi, либо fj имеют разные знаки для разных вершин);

• 2 – обе линии нулевого уровня пересекают стороны треугольника, но точка ($N_{{xs(0)}}^{{ij}},N_{{ys(0)}}^{{ij}}$) находится вне треугольника (значения fi и fj имеют разные знаки для разных вершин, но условия (17) не выполняются);

• 3 – точка ($N_{{xs(0)}}^{{ij}},N_{{ys(0)}}^{{ij}}$) лежит внутри треугольника, но одна из сторон треугольника не пересекается ни одной из линий нулевого уровня (условия (17) выполняются, но суммы знаков произведений значений fi и fj для вершин каждой из сторон нигде не равны нулю);

• 4 – точка ($N_{{xs(0)}}^{{ij}},N_{{ys(0)}}^{{ij}}$) лежит внутри треугольника и все три стороны треугольника пересекаются линиями нулевого уровня (условия (17) выполняются и есть стороны, для которых суммы знаков произведений значений fi и fj для их вершин равны нулю).

Ранг треугольника определяется по максимальному значению подрангов. Подобное ранжирование удобно для уменьшения объема вычислений. При наличии треугольников с рангом 4 и 3 поиск решения следует начать с них. Затем, при отсутствии результата или для проверки его полноты, следует перейти к треугольникам с рангом 2. Обычно этого достаточно для нахождения всех приемлемых решений. К треугольникам с рангом 1 следует обращаться только в случае отсутствия решений, для того, чтобы в этом окончательно убедиться.

МЕТОД НЕЛДЕРА–МИДА

После того, как выполнено ранжирование треугольников, можно переходить к поиску решений. Для этого хорошо подходит метод Нелдера–Мида (Nelder, Mead, 1965; Кузнецов, 2021). Описание алгоритма и текст программы на языке Фортран представлены в статьях (Потемкин, 2005; Арушанян, 2006).

ПОИСК РЕШЕНИЙ И ПРЕДСТАВЛЕНИЕ НАБЛЮДЕНИЙ

Как было описано в разделе “Определение рангов…” поиск решений осуществляется по всем треугольникам с рангами 4, 3 и 2. Полученные решения необходимо попарно сравнить и отбросить кратные. Условием кратности будет расстояние между точками меньше заданного малого числа. Затем мы переводим координаты N в квадрате на координаты на круге (10) и от эклиптических координат переходим к экваториальным ${{{\mathbf{N}}}^{{eqv}}} = \left\{ {N_{x}^{{eqv}},N_{y}^{{eqv}},N_{z}^{{eqv}}} \right\}$. Подставляя Neqv в (6), получим значения n топоцентрических расстояний ρi (i = 1, 2, …, n). После их подстановки в (1) определим векторы гелиоцентрического положения для всех моментов времени. С этими данными можно вычислить орбиты для всех найденных решений. Представление наблюдений наилучшим образом (значения “O–C”) позволит выбрать искомое решение.

ПРИМЕРЫ

В качестве трех численных примеров рассмотрим определение орбиты межзвездной кометы 2I Борисова на короткой дуге (в уравнениях (2), (3) и (5) берутся только верхние знаки).

Определение орбиты по трем наблюдениям

В первом примере, как и ранее, рассмотрим самый простой вариант – определение орбиты по трем наблюдениям.

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

Построим графики уравнений (3) в “квадратных” координатах по осям Nxs и Nys. Закрасим серым цветом области с отрицательными топоцентрическими расстояниями и параметром орбиты p, а также там, где (11) не выполняется, соответственно, области возможных решений будут окрашены белым (рис. 5, рис. 6).

Рис. 5.

Графики уравнений (8) для кометы Борисова. Общий вид.

Рис. 6.

Окрестности точек пересечения кривых (12) для кометы Борисова.

На рис. 5 в левой части графика расположена большая область белого цвета. На рис. 6 представлена ее правая часть. Найдем координаты девяти точек пересечения (12) из решения системы (13) (табл. 2).

Таблица 1.  

Наблюдения кометы Борисова

t (UT)
(год, месяц, день)
α (2000)
(час, мин, с)
δ (2000)
(град, мин, с)
Обсерватория
2019 09 08.630 642 08 44 37.105 +30 57 54.54 Mauna Kea
2019 09 28.234 820 09 21 15.673 +24 15 06.14 Optical Ground Station, Tenerife
2019 10 18.14 757 09 57 58.12 +15 25 08.0 University of Szeged, Piszkesteto Stn. (Konkoly)
2019 09 18.146 976 09 02 26.422 +27 55 25.06 Observatori Astronomic del Montsec
Таблица 2.  

Значения Nij (i, j = 1, 2, 3), полученные из уравнений (12) для кометы Борисова

j\i 1 2 3
1 (–0.0836, –0.3138) (0.0200, –0.2477) (0.0907, –0.2003)
2 (–0.0718, –0.2711) (0.0152, –0.1893) (0.0640, –0.1418)
3 (–0.0376, –0.1439) (0.0063, –0.0791) (0.0238, –0.0530)

Все девять точек расположены в прямоугольнике [–0.0836, 0.0907] × [–0.3138, –0.0530]. Увеличим границы до величин, кратных 0.0625: [–0.1250, 0.1250] × [–0.3750, 0.1250] – получим прямоугольник с соотношением сторон 2 : 3. Его удобно разбить на 8 × 12 “малых” квадратов и 384 “малых” треугольника, соответственно. Доля “больших” квадратов составит 1000, а треугольников – 4000. После проверки на покрытие областей возможных решений останется 1021 треугольник: 937 “больших” и 84 “малых” (рис. 7, рис. 8). Триангуляция по правому краю рис. 7 соответствует левой границе области решений вследствие проекции трехмерной сферы на плоскость.

Рис. 7.

Триангуляция областей возможных решений (8) для кометы Борисова. Общий вид.

Рис. 8.

Триангуляция малых областей возможных решений (8) для кометы Борисова.

Ранжирование треугольников дает: 1 – третьего ранга, 21 – второго, 83 – первого и 916 – нулевого.

Поиск минимумов целевой функции (14) производился методом Нелдера–Мида при следующих заданных параметрах: α = 1, β = 0.5, γ = 2, ε = 10–16.

Все треугольники третьего и второго ранга дали решения.

Реальных решений, конечно же, меньше 22 и необходимо отбросить все дубликаты, оставив в каждом случае лишь одно с наименьшим значением целевой функции. В качестве минимального расстояния между различными решениями была принята величина 10–5. После его проверки были получены три набора решений. Минимальные значения целевая функция принимает в решениях, представленных в табл. 3.

Таблица 3.  

Решения для кометы Борисова

Nxs Nys fgoal
1 –0.85067 –0.25668 2.60 × 10–17
2 –0.58011 –0.47090 2.5 × 10–17
3 –0.00008 0.00012 8.80 × 10–17

Из табл. 3 легко увидеть, что все три минимума являются глобальными и незначительно отличаются от нуля.

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

Таблица 4.  

Значения ρ1, ρ2 и ρ3, полученные из (6) для различных орбит кометы Борисова

Орбита ρ1 (а. е.) ρ2 (а. е.) ρ3 (а. е.)
“1” 1.82702 1.64577 1.39277
“2” 3.50614 3.07836 2.67465
“3” 0.00048 0.00053 0.00056

Как видно из табл. 4, имеются три решения, причем одно из них близко к орбите Земли. Элементы орбиты этих решений представлены в табл. 5, где в качестве эталона приведены элементы орбиты кометы Борисова на эпоху 2020 05 31.0, полученные с учетом всех возмущений в Minor Planet Center (MPC, 2020).

Таблица 5.  

Элементы решения орбит для кометы Борисова

Орбита T0 /М0 (угл. град) a
(а. е.)
e i
(угл. град)
ω
(угл. град)
Ω
(угл. град)
“1” 197.847 0.7856 0.616 59.464 341.862 283.772
“2” 2019 12 08.59 –0.853 3.350 44.063 209.153 308.136
“3” 261.368 0.999 0.017 0.008 250.832 214.262
MPC 2019 12 08.55 –0.851 3.357 44.053 209.127 308.149

Здесь: T0 – момент прохождения перигелия (год, месяц, день). У нас имеется эллиптическая орбита “1”, гиперболическая “2” и орбита, соответствующая орбите Земли “3”. Элементы второго решения хорошо совпадают с элементами эталонной орбиты.

Теперь рассмотрим представление наблюдений орбитами, как показано в табл. 6.

Таблица 6.  

Представление наблюдений кометы Борисова полученными орбитами

Орбита t1 t2 t3
Δα
(угл. с)
Δδ
(угл. с)
Δα
(угл. с)
Δδ
(угл.с)
Δα
(угл.с)
Δδ
(угл.с)
“1” 0.0 –4.6 × 10–10 0.035 –0.009 4.4 × 10–10 3.0 × 10–10
“2” 7.9 × 10–11 –1.1 × 10–10 –0.001 0.003 1.8 × 10–10 –1.1 × 10–10
“3” 1.6 × 10–7 4.7 × 10–8 13.520 16.060 1.2 × 10–7 5.5 × 10–10

Значения (“O–C”), для момента t2 в табл. 6, приемлемы для всех трех решений. Особенно это касается решений “1” и “2”. Для выбора между ними необходимо привлечь наблюдение для момента t4. Значения (“O–C”) в момент t4 показаны в табл. 7.

Таблица 7. 

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

Орбита t4
Δα (угл. с) Δδ (угл. с)
“1” 139.75 305.40
“2” 2.19 1.47
“3” –3.6 × 104 –2.8 × 104

Орбита “2” лучше других представляет четвертое наблюдение и может считаться искомой.

Определение орбиты по четырем наблюдениям

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

Рассмотрим задачу определения орбиты по четырем наблюдениям. Так наблюдения, представленные в табл. 1, частично изменят порядок: четвертое станет вторым, а второе и третье соответственно третьим и четвертым.

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

После проверки на покрытие областей возможных решений, число треугольников увеличится до 1071, за счет роста числа “больших” до 987, при неизменных 84 “малых”. Изменения произойдут и в ранжировании треугольников: 3 – третьего ранга, 47 – второго ранга, 80 – первого и 941 – нулевого.

Поиск минимумов целевой функции (14) производился методом Нелдера–Мида с теми же параметрами, что и в предыдущем примере. После проверки минимальных расстояний между решениями были получены три набора. Минимальные значения целевая функция принимает в решениях, представленных в табл. 8.

Таблица 8.  

Решения для кометы Борисова

Nxs N fgoal
1 –0.85684 –0.21368 3.22× 10–1
2 –0.57989 –0.47090 1.68 × 10–6
3 –0.00013 0.00021 3.87 × 10–3

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

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

Таблица 9.  

Значения ρ1, ρ2, ρ3 и ρ4, полученные из (6) для различных орбит комет Борисова

Орбита ρ1 (а. е.) ρ2 (а. е.) ρ3 (а. е.) ρ4 (а. е.)
“1” 1.73609 1.66563 1.56634 1.32386
“2” 3.50775 3.29747 3.07972 2.67592
“3” 0.00081 0.00073 0.00084 0.00079

Как видно из табл. 9, для трех решений, по сравнению с предыдущим примером (табл. 4), больше всего изменились топоцентрические расстояния для первого решения. Элементы орбиты этих решений представлены в табл. 10. Сравнение ее с табл. 5 демонстрирует, что и здесь сильнее всего изменились элементы орбиты “1” и то, что орбита “2” стала ближе к решению MPC.

Таблица 10.  

Элементы решения орбит для кометы Борисова

Орбита T0 / М0
(угл. град)
a
(а. е.)
e i
(угл. град)
ω
(угл. град)
Ω
(угл. град)
“1” 172.455 0.7257 0.659 59.757 344.962 281.349
“2” 2019 12 08.56 –0.851 3.357 44.052 209.133 308.149
“3” 249.197 0.998 0.017 0.014 254.427 212.879
MPC 2019 12 08.55 –0.851 3.357 44.053 209.127 308.149

Теперь рассмотрим представление наблюдений орбитами, как показано в табл. 11. Если сравнивать представленные в ней результаты с табл. 6, то t4 соответствует t3, а значения для t3 нужно сравнивать с данными из табл. 7.

Таблица 11.  

Представление наблюдений кометы Борисова полученными орбитами

Орбита t1 t2 t3 t4
Δα
(угл. с)
Δδ
(угл. с)
Δα
(угл. с)
Δδ
(угл. с)
Δα
(угл. с)
Δδ
(угл. с)
Δα
(угл. с)
Δδ
(угл. с)
“1” 0.0 0.0 –35.76 61.81 –228.0 –322.7 1.8 × 10–10 1.2 × 10–10
“2” 1.6 × 10–10 6.9 × 10–11 7.8 × 10–6 –7.2 × 10–4 –2.81 –1.40 3.5 × 10–10 –2.4 × 10–10
“3” 7.4 × 10–8 8.5 × 10–8 –4.1 × 103 –1.5 × 104 1.7 × 104 110.7 1.4 × 10–7 –2.8 × 10–8

Значения (“O–C”), для моментов t2 и t3 в табл. 11, приемлемы только для второго решения. Третье решение, несмотря на то, что значение целевой функции для него меньше, явно уступает в представлении этих наблюдений первому.

Определение орбиты по пяти наблюдениям

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

Таблица 12.  

Пятое наблюдение кометы Борисова

t (UT)
(год, месяц, день)
α (2000)
(час, мин, с)
δ (2000)
(град, мин, с)
Обсерватория
2019 10 08.18 294 09 39 41.06 +20 07 17.3 G. Pascoli Observatory, Castelvecchio Pascoli

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

Область возможных решений оказалась близкой к области в первом примере, как и ее триангуляция: 1021 треугольник: 937 “больших” и 84 “малых”. Изменения произошли в ранжировании треугольников: 1 – четвертого ранга, 5 – третьего, 53 – второго, 84 – первого и 878 – нулевого.

Поиск минимумов целевой функции (14) проводился методом Нелдера–Мида при тех же параметрах, что и ранее.

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

Таблица 13.  

Решения для кометы Борисова

Nxs Nys fgoal
1 –0.84139 –0.23417 0.99580
2 –0.58006 –0.47090 9.36 × 10–7
3 –5.96 × 10–5 9.16 × 10–5 3.84 × 10–3

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

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

Таблица 14.  

Значения ρ1, ρ2, ρ3, ρ4 и ρ5, полученные из (6) для различных орбит комет Борисова

Орбита ρ1 (а. е.) ρ2 (а. е.) ρ3 (а. е.) ρ4 (а. е.) ρ5 (а. е.)
“1” 1.76527 1.69786 1.60078 1.48674 1.36104
“2” 3.50661 3.29627 3.07862 2.87126 2.67491
“3” 0.00034 0.00026 0.00039 0.00032 0.00038

Элементы орбиты этих решений представлены в табл. 15.

Таблица 15.  

Элементы решения орбит для кометы Борисова

Орбита T0/М0
(угл. град)
a (а. е.) e i
(угл. град)
ω
(угл. град)
Ω
(угл. град)
“1” 188.372 0.7532 0.629 58.304 344.678 282.780
“2” 2019 12 08.58 –0.853 3.351 44.061 209.145 308.139
“3” 261.820 0.999 0.017 0.006 251.579 213.058
MPC 2019 12 08.55 –0.851 3.357 44.053 209.127 308.149

Теперь рассмотрим представление наблюдений орбитами, как показано в табл. 16. Здесь сохраняется отношение точностей представления центральных наблюдений по решениям, не совпадающее с отношением значений целевой функции для случаев “1” и “3”.

Таблица 16.  

Представление наблюдений кометы Борисова полученными орбитами

t1 t2 t3 t4 t5
Δα Δδ Δα Δδ Δα Δδ Δα Δδ Δα Δδ
1 0.0 –4 × 10–11 103 347 –34.8 66.7 –157 –251 9 × 10–11 –9 × 10–11
2 0.0 1 × 10–10 2.18 1.57 –0.008 0.21 1.49 0.65 9 × 10–11 –8 × 10–11
3 8 × 10–8 8 × 10–7 7 × 104 4 × 104 3 × 103 3 × 103 8 × 103 666 1 × 10–7 –4 × 10–6

Примечание: единицы измерений во всех колонках – “угл. с”.

Как и ранее, значения (O–C), для моментов t2, t3 и t4 в табл. 16 приемлемы только для второго решения. Дальнейшее увеличение числа наблюдений внутри интервала для определения орбиты не изменит настоящей картины, а только лишь увеличит объем вычислений.

ЗАКЛЮЧЕНИЕ

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

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

  1. Apyшaнян O.Б. http://num-anal.srcc.msu.ru/lib_na/cat/mn/mnb6r.htm, 2006.

  2. Кузнецов В.Б. К вопросу об определении предварительной орбиты небесного тела // Астрон. вестн. 2019. Т. 53. № 6. С. 456–466. (Kuznetsov V.B. Parabolic orbit determination. Camparison of the Olbers method and algebraic equations // Sol. Syst. Res. 2019. V. 53. № 6. P. 462–472.)

  3. Кузнецов В.Б. Определение параболической орбиты. Поиск решения в методе алгебраических уравнений // Астрон. вестн. 2021. Т. 55. № 3. С. 1–9. (Kuznetsov V.B. Defenition of a parabolic orbit: the search of a solution with algebraic equations // Sol. Syst. Res. 2021. V. 55. № 4. P. 359–367.)

  4. Пoтeмкин B. http://fortran-90.pvbk.spb.ru/min.html#FM28, 2005.

  5. Самотохин А.С., Хуторовский З.Н. Метод первоначального определения параметров околоземных орбит по трем угловым измерениям // Препр. ИПМ им. М.В. Келдыша. 2014. № 44. 31 с. URL: http://library.keldysh.ru/preprint.asp?id=2014-44

  6. Химмельблау Д. Прикладное нелинейное программирование. М.: Мир, 1975. 536 с.

  7. Шалашилин В.И., Кузнецов Е.Б. Метод продолжения решения по параметру и наилучшая параметризация (в прикладной математике и механике). М.: Эдиториал УРСС, 1999. 224 с.

  8. Шефер В.А. Новый метод определения орбиты по двум векторам положения, основанный на решении уравнений Гаусса // Астрон. вестн. 2010. Т. 44. № 3. С. 273–288. (Shefer V.A. New method of orbit determination from two position vectors based on solving Gauss’s equations // Sol. Syst. Res. 2010. V. 44. № 3. P. 273–288.)

  9. Fong C. Analytical methods for squaring the disc // Seoul ICM 2014.

  10. Nelder J.A., Mead R. A simplex method for function minimization // Computer J. 1965. V. 7 (4). P. 308–313.

  11. The International Astronomical Union. Minor Planet Center, 2020, https://minorplanetcenter.net/db_search/show_object?utf8=%E2%9C%93&object_id=2I

  12. Sarnecki A.J. A projective approach to orbit determination from three sight–lines // Celest. Mech. and Dyn. Astron. 1997. V. 66. P. 425–451.

  13. Stumpff K. Himmelsmechanik. Bd. 1. Berlin, 1959. 682 p.

Дополнительные материалы отсутствуют.