Ме́тод Га́усса в небесной механике и астродинамике используется для первоначального определения параметров орбиты небесного тела по трём наблюдениям.
На практике для увеличения точности используется больше наблюдений, но в теории достаточно трёх. Кроме небесных координат объекта, необходимыми сведениями являются моменты наблюдений и земные координаты пунктов наблюдения.
В 1801 году была открыта Церера , но в течение некоторого времени её наблюдения были затруднены из-за близости к Солнцу, после чего было трудно снова найти её на небе. Карл Фридрих Гаусс поставил себе задачу определения её орбиты по имевшимся наблюдениям, за счёт чего и приобрёл мировую известность[1] . Однако описанный ниже метод годится только для определения орбит с фокусом в теле, с которого ведутся наблюдения, так что задача Гаусса была сложнее.
Геоцентрическая (θ) и геодезическая (ϕ) широта Вектор положения наблюдателя (в экваториальной системе координат ) можно вычислить, зная широту места наблюдения и местное звёздное время :
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 = 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}} — моменты наблюдений.
Вычисляются векторные произведения :
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}} } Вычисляются смешанные произведения :
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}} } Вычисляются позиционные коэффициенты:
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}} } Вычисляется модуль вектора положения наблюдателя в момент второго наблюдения:
R 2 2 = R 2 ⋅ R 2 {\displaystyle {R_{2}}^{2}=\mathbf {R_{2}} \cdot \mathbf {R_{2}} } Вычисляются коэффициенты полинома для поиска расстояния:
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 } — гравитационный параметр тела, вокруг которого происходит вращение.
Ищутся решения уравнения:
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}} — расстояние до объекта в момент второго наблюдения.
У кубического уравнения может быть до трёх действительных корней. В случае, если их больше одного, необходимо проверить каждый из них.
Вычисляются расстояния от точек наблюдения до объекта в каждый из моментов наблюдений:
ρ 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]} Вычисляются позиционные вектора объекта (в экваториальной системе координат ):
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}} } Вычисляются коэффициенты Лагранжа . Из-за этого пункта определение орбит становится неточным:
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}} Вычисляется вектор скорости объекта в момент второго наблюдения (в экваториальной системе координат):
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)} Теперь известно положение и скорость объекта в один момент времени. Значит, возможно определить параметры орбиты[2] .
Der, Gim J.. «New Angles-only Algorithms for Initial Orbit Determination.» Advanced Maui Optical and Space Surveillance Technologies Conference. (2012). Print (англ.)