Ме́тод Га́усса
в
небесной механике
и
астродинамике
используется для первоначального определения
параметров орбиты
небесного тела по трём наблюдениям.
На практике для увеличения точности используется больше наблюдений, но в теории достаточно трёх. Кроме
небесных координат
объекта, необходимыми сведениями являются моменты наблюдений и земные координаты пунктов наблюдения.
История
В 1801 году была открыта
Церера
, но в течение некоторого времени её наблюдения были затруднены из-за близости к Солнцу, после чего было трудно снова найти её на небе.
Карл Фридрих Гаусс
поставил себе задачу определения её орбиты по имевшимся наблюдениям, за счёт чего и приобрёл мировую известность
. Однако описанный ниже метод годится только для определения орбит с фокусом в теле, с которого ведутся наблюдения, так что задача Гаусса была сложнее.
Вектор положения наблюдателя
Геоцентрическая (θ) и геодезическая (ϕ) широта
Вектор положения наблюдателя (в
экваториальной системе координат
) можно вычислить, зная
широту
места наблюдения и местное
звёздное время
:
R
n
=
[
R
e
1
−
(
2
f
−
f
2
)
sin
2
ϕ
n
+
H
n
]
cos
ϕ
n
(
cos
θ
n
I
^
+
sin
θ
n
J
^
)
+
[
R
e
(
1
−
f
)
2
1
−
(
2
f
−
f
2
)
sin
2
ϕ
n
+
H
n
]
sin
ϕ
n
K
^
{\displaystyle \mathbf {R_{n}} =\left[{R_{e} \over {\sqrt {1-(2f-f^{2})\sin ^{2}\phi _{n}}}}+H_{n}\right]\cos \phi _{n}(\cos \theta _{n}\ \mathbf {\hat {I}} +\sin \theta _{n}\ \mathbf {\hat {J}} )+\left[{R_{e}(1-f)^{2} \over {\sqrt {1-(2f-f^{2})\sin ^{2}\phi _{n}}}}+H_{n}\right]\sin \phi _{n}\ \mathbf {\hat {K}} }
или:
R
n
=
R
e
cos
ϕ
n
′
cos
θ
n
I
^
+
R
e
cos
ϕ
n
′
sin
θ
n
J
^
+
R
e
sin
ϕ
n
′
K
^
,
{\displaystyle \mathbf {R_{n}} =R_{e}\cos \phi '_{n}\cos \theta _{n}\ \mathbf {\hat {I}} +R_{e}\cos \phi '_{n}\sin \theta _{n}\ \mathbf {\hat {J}} +R_{e}\sin \phi '_{n}\ \mathbf {\hat {K}} ,}
где:
R
n
{\displaystyle \mathbf {R_{n}} }
— вектор положения наблюдателя;
R
e
{\displaystyle R_{e}}
— экваториальный радиус тела, на котором находится наблюдатель;
f
{\displaystyle f}
— сплюснутость тела у полюсов (например, для
Земли
— 0.003353);
ϕ
n
{\displaystyle \phi _{n}}
— геодезическая широта;
ϕ
n
′
{\displaystyle \phi '_{n}}
— геоцентрическая широта;
H
n
{\displaystyle H_{n}}
— высота;
θ
n
{\displaystyle \theta _{n}}
— местное звёздное время.
Вектор направления на объект
Вектор направления на объект может быть вычислен с помощью
склонения
и
прямого восхождения
:
ρ
^
n
=
cos
δ
n
cos
α
n
I
^
+
cos
δ
n
sin
α
n
J
^
+
sin
δ
n
K
^
{\displaystyle \mathbf {{\hat {\rho }}_{n}} =\cos \delta _{n}\cos \alpha _{n}\ \mathbf {\hat {I}} +\cos \delta _{n}\sin \alpha _{n}\ \mathbf {\hat {J}} +\sin \delta _{n}\ \mathbf {\hat {K}} }
,
где:
ρ
^
n
{\displaystyle \mathbf {{\hat {\rho }}_{n}} }
— единичный вектор направления на объект;
δ
n
{\displaystyle \delta _{n}}
— склонение;
α
n
{\displaystyle \alpha _{n}}
— прямое восхождение.
Определение орбиты
Далее нужно получить вектор расстояния до объекта, а не только единичный вектор направления на него.
Шаг 1
Вычисляются интервалы между наблюдениями:
τ
1
=
t
1
−
t
2
{\displaystyle \tau _{1}=t_{1}-t_{2}}
τ
3
=
t
3
−
t
2
{\displaystyle \tau _{3}=t_{3}-t_{2}}
τ
=
t
3
−
t
1
,
{\displaystyle \tau =t_{3}-t_{1},}
где
t
n
{\displaystyle t_{n}}
— моменты наблюдений.
Шаг 2
Вычисляются
векторные произведения
:
p
1
=
ρ
^
2
×
ρ
^
3
{\displaystyle \mathbf {p_{1}} =\mathbf {{\hat {\rho }}_{2}} \times \mathbf {{\hat {\rho }}_{3}} }
p
2
=
ρ
^
1
×
ρ
^
3
{\displaystyle \mathbf {p_{2}} =\mathbf {{\hat {\rho }}_{1}} \times \mathbf {{\hat {\rho }}_{3}} }
p
3
=
ρ
^
1
×
ρ
^
2
{\displaystyle \mathbf {p_{3}} =\mathbf {{\hat {\rho }}_{1}} \times \mathbf {{\hat {\rho }}_{2}} }
Шаг 3
Вычисляются
смешанные произведения
:
D
0
=
ρ
^
1
⋅
p
1
=
ρ
^
1
⋅
(
ρ
^
2
×
ρ
^
3
)
{\displaystyle D_{0}=\mathbf {{\hat {\rho }}_{1}} \cdot \mathbf {p_{1}} =\mathbf {{\hat {\rho }}_{1}} \cdot (\mathbf {{\hat {\rho }}_{2}} \times \mathbf {{\hat {\rho }}_{3}} )}
D
11
=
R
1
⋅
p
1
D
12
=
R
1
⋅
p
2
D
13
=
R
1
⋅
p
3
{\displaystyle D_{11}=\mathbf {R_{1}} \cdot \mathbf {p_{1}} \qquad D_{12}=\mathbf {R_{1}} \cdot \mathbf {p_{2}} \qquad D_{13}=\mathbf {R_{1}} \cdot \mathbf {p_{3}} }
D
21
=
R
2
⋅
p
1
D
22
=
R
2
⋅
p
2
D
23
=
R
2
⋅
p
3
{\displaystyle D_{21}=\mathbf {R_{2}} \cdot \mathbf {p_{1}} \qquad D_{22}=\mathbf {R_{2}} \cdot \mathbf {p_{2}} \qquad D_{23}=\mathbf {R_{2}} \cdot \mathbf {p_{3}} }
D
31
=
R
3
⋅
p
1
D
32
=
R
3
⋅
p
2
D
33
=
R
3
⋅
p
3
{\displaystyle D_{31}=\mathbf {R_{3}} \cdot \mathbf {p_{1}} \qquad D_{32}=\mathbf {R_{3}} \cdot \mathbf {p_{2}} \qquad D_{33}=\mathbf {R_{3}} \cdot \mathbf {p_{3}} }
Шаг 4
Вычисляются позиционные коэффициенты:
A
=
1
D
0
(
−
D
12
τ
3
τ
+
D
22
+
D
32
τ
1
τ
)
{\displaystyle A={\frac {1}{D_{0}}}\left(-D_{12}{\frac {\tau _{3}}{\tau }}+D_{22}+D_{32}{\frac {\tau _{1}}{\tau }}\right)}
B
=
1
6
D
0
[
D
12
(
τ
3
2
−
τ
2
)
τ
3
τ
+
D
32
(
τ
2
−
τ
1
2
)
τ
1
τ
]
{\displaystyle B={\frac {1}{6D_{0}}}\left[D_{12}\left(\tau _{3}^{2}-\tau ^{2}\right){\frac {\tau _{3}}{\tau }}+D_{32}\left(\tau ^{2}-\tau _{1}^{2}\right){\frac {\tau _{1}}{\tau }}\right]}
E
=
R
2
⋅
ρ
^
2
{\displaystyle E=\mathbf {R_{2}} \cdot \mathbf {{\hat {\rho }}_{2}} }
Шаг 5
Вычисляется модуль вектора положения наблюдателя в момент второго наблюдения:
R
2
2
=
R
2
⋅
R
2
{\displaystyle {R_{2}}^{2}=\mathbf {R_{2}} \cdot \mathbf {R_{2}} }
Шаг 6
Вычисляются коэффициенты полинома для поиска расстояния:
a
=
−
(
A
2
+
2
A
E
+
R
2
2
)
{\displaystyle a=-\left(A^{2}+2AE+{R_{2}}^{2}\right)}
b
=
−
2
μ
B
(
A
+
E
)
{\displaystyle b=-2\mu B(A+E)}
c
=
−
μ
2
B
2
,
{\displaystyle c=-\mu ^{2}B^{2},}
где
μ
{\displaystyle \mu }
—
гравитационный параметр
тела, вокруг которого происходит вращение.
Шаг 7
Ищутся решения уравнения:
r
2
8
+
a
r
2
6
+
b
r
2
3
+
c
=
0
,
{\displaystyle {r_{2}}^{8}+a{r_{2}}^{6}+b{r_{2}}^{3}+c=0,}
где
r
2
{\displaystyle r_{2}}
— расстояние до объекта в момент второго наблюдения.
У кубического уравнения может быть до трёх действительных корней. В случае, если их больше одного, необходимо проверить каждый из них.
Шаг 8
Вычисляются расстояния от точек наблюдения до объекта в каждый из моментов наблюдений:
ρ
1
=
1
D
0
[
6
(
D
31
τ
1
τ
3
+
D
21
τ
τ
3
)
r
2
3
+
μ
D
31
(
τ
2
−
τ
1
2
)
τ
1
τ
3
6
r
2
3
+
μ
(
τ
2
−
τ
3
2
)
−
D
11
]
{\displaystyle \rho _{1}={\frac {1}{D_{0}}}\left[{\frac {6\left(D_{31}{\dfrac {\tau _{1}}{\tau _{3}}}+D_{21}{\dfrac {\tau }{\tau _{3}}}\right){r_{2}}^{3}+\mu D_{31}\left(\tau ^{2}-{\tau _{1}}^{2}\right){\dfrac {\tau _{1}}{\tau _{3}}}}{6{r_{2}}^{3}+\mu \left(\tau ^{2}-{\tau _{3}}^{2}\right)}}-D_{11}\right]}
ρ
2
=
A
+
μ
B
r
2
3
{\displaystyle \rho _{2}=A+{\frac {\mu B}{{r_{2}}^{3}}}}
ρ
3
=
1
D
0
[
6
(
D
13
τ
3
τ
1
−
D
23
τ
τ
1
)
r
2
3
+
μ
D
13
(
τ
2
−
τ
3
2
)
τ
3
τ
1
6
r
2
3
+
μ
(
τ
2
−
τ
1
2
)
−
D
33
]
{\displaystyle \rho _{3}={\frac {1}{D_{0}}}\left[{\frac {6\left(D_{13}{\dfrac {\tau _{3}}{\tau _{1}}}-D_{23}{\dfrac {\tau }{\tau _{1}}}\right){r_{2}}^{3}+\mu D_{13}\left(\tau ^{2}-{\tau _{3}}^{2}\right){\dfrac {\tau _{3}}{\tau _{1}}}}{6{r_{2}}^{3}+\mu \left(\tau ^{2}-{\tau _{1}}^{2}\right)}}-D_{33}\right]}
Шаг 9
Вычисляются позиционные вектора объекта (в
экваториальной системе координат
):
r
1
=
R
1
+
ρ
1
ρ
^
1
{\displaystyle \mathbf {r_{1}} =\mathbf {R_{1}} +\rho _{1}\mathbf {{\hat {\rho }}_{1}} }
r
2
=
R
2
+
ρ
2
ρ
^
2
{\displaystyle \mathbf {r_{2}} =\mathbf {R_{2}} +\rho _{2}\mathbf {{\hat {\rho }}_{2}} }
r
3
=
R
3
+
ρ
3
ρ
^
3
{\displaystyle \mathbf {r_{3}} =\mathbf {R_{3}} +\rho _{3}\mathbf {{\hat {\rho }}_{3}} }
Шаг 10
Вычисляются
коэффициенты Лагранжа
. Из-за этого пункта определение орбит становится неточным:
f
1
≈
1
−
1
2
μ
r
2
3
τ
1
2
{\displaystyle f_{1}\approx 1-{\frac {1}{2}}{\frac {\mu }{{r_{2}}^{3}}}{\tau _{1}}^{2}}
f
3
≈
1
−
1
2
μ
r
2
3
τ
3
2
{\displaystyle f_{3}\approx 1-{\frac {1}{2}}{\frac {\mu }{{r_{2}}^{3}}}{\tau _{3}}^{2}}
g
1
≈
τ
1
−
1
6
μ
r
2
3
τ
1
3
{\displaystyle g_{1}\approx \tau _{1}-{\frac {1}{6}}{\frac {\mu }{{r_{2}}^{3}}}{\tau _{1}}^{3}}
g
3
≈
τ
3
−
1
6
μ
r
2
3
τ
3
3
{\displaystyle g_{3}\approx \tau _{3}-{\frac {1}{6}}{\frac {\mu }{{r_{2}}^{3}}}{\tau _{3}}^{3}}
Шаг 11
Вычисляется вектор скорости объекта в момент второго наблюдения (в экваториальной системе координат):
v
2
=
1
f
1
g
3
−
f
3
g
1
(
−
f
3
r
1
+
f
1
r
3
)
{\displaystyle \mathbf {v_{2}} ={\frac {1}{f_{1}g_{3}-f_{3}g_{1}}}\left(-f_{3}\mathbf {r_{1}} +f_{1}\mathbf {r_{3}} \right)}
Шаг 12
Теперь известно положение и скорость объекта в один момент времени. Значит, возможно определить параметры орбиты
.
Примечания
(неопр.)
. Дата обращения: 11 марта 2020.
15 мая 2012 года.
(неопр.)
. Дата обращения: 11 марта 2020.
10 ноября 2020 года.
Литература
Der, Gim J.. «New Angles-only Algorithms for Initial Orbit Determination.» Advanced Maui Optical and Space Surveillance Technologies Conference. (2012). Print
(англ.)