“It is not knowledge, but the act of learning, not possession, but the act of getting there, which grants the greatest enjoyment.”
– Carl Friedrich Gauss
这篇博客中主要的计算过程翻译自 Bachelor Thesis of Daniel Bed’at
s
ˇ
\check{\text{s}}
sˇ, Gauss’ calculation of Ceres’ orbit
Many thanks to him for allowing me to translate this work to Chinese. 感谢他同意我翻译他的论文产生了这篇博客。本文过长所以分为上下两篇。参考文献见下篇(Reference is at the end of second part).
高斯
大家可能都听过的一个故事是高斯7岁上小学的时候就发现了如何快速计算从1到100的和。
统计学的早期历史可以追溯到1795年,高斯18岁时,发明了最小二乘法和正态分布。后来1801年,Piazzi发现谷神星 Ceres 观测几十天然后跟丢后,当年高斯第一个根据有限的观测数据计算了谷神星的轨道并预测了Ceres的位置。1805年,Legendre勒让德发表的文章首次介绍了最小二乘法。高斯的文章发表于1809年,不过他声称自己1795年就知道怎么用最小二乘法。
很多科普的文章和视频介绍最小二乘法或正态分布时都会说这是高斯计算谷神星的轨道时用的方法。有的作者还会拿最小二乘法演示一下如何用多项式拟合一些数据。实际上高斯到底是解决了一个别人都没有解决的什么样的问题,他又是如何做到的呢?
1801年的1月1日,意大利天文学家 Giuseppe Piazzi 发现了一颗小行星。在观测了40天后到2月11日这颗小行星就因为天气和太阳的接近而无法追踪它的轨迹。他把它命名为Ceres Ferdinandea。人们面临的首要问题就是如何基于当时的少量观测数据来估计这颗小行星的轨道参数。
绘制出来这些点大概是下面这样。
开普勒三定律那时已经被天文学家们知晓,微积分也建立了。人们已经可以想象,自己是在绕着太阳做椭圆运动的地球上,观测绕着太阳做不同形状椭圆运动的行星,因此上面的观测数据的轨迹并不会让人感到太过奇怪。
但天文学家如果想要通过求解开普勒三定律中的椭圆方程来求行星的轨道,在仅有不到完整轨道1%的数据的情况下,这几乎是一个不可能的任务。24名天文学家组成了一个‘The Society for the Detection of a Missing World’,后来被亲切地称为‘Celestial Police’因为他们负有抓捕丢失的行星的责任。法国数学家Laplace认为这个任务是不可能完成的。
数学王子 (The Prince of Mathematics)
24岁的Gauss在基于开普勒三定律是正确的假设下花了三个月的时间解决了这个问题。他花了超过100小时进行大量计算没有出错。具体地来说,他需要基于已知的19个观测数据点(实际上他只用了3个)估计Ceres椭圆轨道的6个参数,如图所示
Source: Teets & Whitehead, The Discovery of Ceres: How Gauss Became Famous (1999)
Source: Math 221 Veronique Le Corvec、Jeffrey Donatelli、Jeffrey Hunt
分别解释下六个轨道参数。
确定行星的轨道平面相对地球轨道平面的位置需要两个参数:
- i:行星轨道平面法向量与地球轨道平面法向量(z轴)倾斜(inclination)的角度
-
Ω
\Omega
Ω: 平面相交线与固定的x轴之间的角度
行星轨道的形状和大小: - a: 椭圆轨道半长轴,基本确定轨道的大小
- e: 椭圆轨道偏心率,确定轨道的形状,用半短轴也可以确定轨道的形状(e=b/a)
然后有一个参数用来确定轨道的长轴在轨道平面内的角度: - ω \omega ω: 轨道长轴近日点与轨道平面交线(也叫line of nodes)的角度
上面五个参数已经可以完全确定轨道的位置。为了计算某个时刻行星的位置,除了以上这些参数,我们还需要知道行星在轨道上某个位置的时刻。
在200 Years of Least Squares Method中,使用的记号是:
- ω \omega ω arg. of perihelion 轨道的近日点相对于上升点连线的角度
- Ω \Omega Ω long. of ascend. node 上升点角度(轨道平面交线/点的角度)
- i inclination of orbit
- a semi-major axis 半长轴
- e eccentricity 偏心率
-
l
0
l_0
l0 mean heliocent. long.
J. Tennenbaum 和 B. Director 在 How Gauss Determined The Orbit of Ceres 中,详细地回顾了Gauss最早的工作 Summary Overview of the Method
Used To Determine the Orbits of the Two New Planets 中的方法。实际上在高斯最早的做法中,他只使用了三个数据点进行计算,并不是直接用最小二乘法基于大量的数据进行参数估计。当时矩阵计算的理论还处在发展之中,高斯消元法也是他发展的。高斯针对非线性的轨道参数计算采用了一些近似技巧,忽略了一些高阶小项。高斯后续的一系列改进的工作使得预测误差更低了,其中包括他后来所说的自己18岁时发明的最小二乘法。
这里参考之前的资料复现一下高斯是如何利用三个数据点计算轨道参数的。
主要参考Charles University 的 Daniel Bed’atˇs 的本科毕业论文 Causs’ calculation of Ceres’ orbit,还有前面提到的 J. Tennenbaum 和 B. Director 在 How Gauss Determined The Orbit of Ceres .
观测数据是
λ
\lambda
λ 和
β
\beta
β,而
ρ
\rho
ρ 是未知的距离。因此需要知道阶段A和B的公式。
Daniel的论文只有两章,第一章介绍椭圆轨道,主要讲轨道参数可以如何基于已有的变量计算出来,第二章讲高斯是如何建立方程来计算需要的参数的。
轨道参数如何计算
i 和 Ω \Omega Ω 的计算
假设我们已经得到两个时刻行星的位置向量
r
\bold{r}
r 和
r
′
′
\bold{r''}
r′′,分别在时刻
τ
\tau
τ 和
τ
′
′
\tau''
τ′′,以及
τ
<
τ
′
′
\tau<\tau''
τ<τ′′. 如何从观测数据得到
r
\bold{r}
r,
r
′
′
\bold{r}''
r′′后面叙述,目前假设我们已经知道了。
首先由于观测数据相隔时间很短,因此
r
\bold{r}
r
r
′
′
\bold{r}''
r′′ 的角度小于 180 度。外积
r
×
r
′
′
r \times r''
r×r′′ 给出了一个法向量
n
\bold{n}
n。记
e
3
=
(
0
,
0
,
1
)
T
e_3=(0,0,1)^T
e3=(0,0,1)T表示z轴的基向量,则
cos
i
=
n
⋅
e
∣
∣
n
∣
∣
∣
∣
e
∣
∣
=
n
3
∣
∣
n
∣
∣
\cos i=\frac{n\cdot e}{||n|| ||e||} =\frac{n_3}{||n||}
cosi=∣∣n∣∣∣∣e∣∣n⋅e=∣∣n∣∣n3
n
3
n_3
n3表示n的第三个坐标。因为i小于90度,所以通过arccos函数即可得到i。
为了求升交点的角度
Ω
\Omega
Ω,我们考虑法向量n在xy平面上的投影
(
n
1
,
n
2
,
0
)
T
(n_1,n_2,0)^T
(n1,n2,0)T. 因此,考虑
n
1
n_1
n1的正负,容易知道为正时投影向量顺时针旋转90度后就得到
Ω
\Omega
Ω,为负时需要旋转270度。
Ω
=
arctan
(
n
2
n
1
)
+
90
,
n
1
>
0
\Omega = \arctan \left( \frac{n_2}{n_1} \right) + 90, n_1 > 0
Ω=arctan(n1n2)+90,n1>0
Ω
=
arctan
(
n
2
n
1
)
+
270
,
n
1
<
0
\Omega = \arctan \left( \frac{n_2}{n_1} \right) + 270, n_1 < 0
Ω=arctan(n1n2)+270,n1<0
当 n 1 = 0 n_1=0 n1=0, 则 Ω = 180 \Omega=180 Ω=180。
计算轨道在平面内的角度 ω \omega ω
首先计算 r 和 升交点的角度
α
\alpha
α,然后利用之后得到的
θ
\theta
θ就可以得到
ω
\omega
ω。
ω
=
α
−
θ
ω
=
360
−
θ
+
α
\omega = \alpha-\theta \\ \omega=360-\theta+\alpha
ω=α−θω=360−θ+α
对于给定时间 t,必须找到角度 u 。
根据开普勒第二定律,“相同时间扫过的面积相同”。记轨道周期为P。
A
a
b
π
=
t
P
\frac{A}{ab\pi}=\frac{t}{P}
abπA=Pt
回顾一下椭圆的极坐标公式
r
(
θ
)
=
k
1
+
e
cos
θ
,
k
=
a
(
1
−
e
2
)
r(\theta)=\frac{k}{1+e\cos\theta} , k=a(1-e^2)
r(θ)=1+ecosθk,k=a(1−e2)
这里
θ
\theta
θ 也叫作位置
r
\bold{r}
r的 true anomaly。
b
2
=
a
2
−
d
2
4
=
a
2
−
4
a
2
e
2
/
4
=
a
2
(
1
−
e
2
)
=
a
k
b^2=a^2-\frac{d^2}{4}=a^2-4a^2e^2/4=a^2(1-e^2)=ak
b2=a2−4d2=a2−4a2e2/4=a2(1−e2)=ak
因此椭圆轨道的面积
π
a
b
=
π
a
3
/
2
k
\pi ab=\pi a^{3/2}\sqrt{k}
πab=πa3/2k.
记轨道周期为
t
p
t_p
tp,则
−
g
′
′
τ
′
′
−
τ
=
π
a
3
/
2
k
t
p
\frac{-g''}{\tau''-\tau}=\frac{\pi a^{3/2}\sqrt{k}}{t_p}
τ′′−τ−g′′=tpπa3/2k
另一方面,开普勒第三定律告诉我们
a
3
t
p
2
=
A
3
T
e
2
\frac{a^3}{t_p^2}=\frac{A^3}{T_e^2}
tp2a3=Te2A3
其中
T
e
T_e
Te 是地球周期,或者一年,A=1AU,是地球轨道的半长轴。上式两边开平方只取正数,则有
a
3
/
2
t
p
=
1
T
p
\frac{a^{3/2}}{t_p}=\frac{1}{T_p}
tpa3/2=Tp1
代入前面,
−
g
′
τ
′
′
−
τ
=
π
k
T
e
\frac{-g'}{\tau''-\tau}=\frac{\pi\sqrt{k}}{T_e}
τ′′−τ−g′=Teπk
将k独立出来则有
k
=
(
−
g
′
T
e
π
(
τ
′
′
−
τ
)
)
2
k=(\frac{-g'T_e}{\pi (\tau''-\tau)})^2
k=(π(τ′′−τ)−g′Te)2
下一步我们就需要估计
−
g
′
-g'
−g′,然后就可以得到 k 。
我们用极坐标下的面积积分来表示
−
g
′
=
1
2
∫
θ
θ
′
′
r
(
α
)
2
d
α
-g'=\frac{1}{2}\int_{\theta}^{\theta''}r(\alpha)^2d\alpha
−g′=21∫θθ′′r(α)2dα
其中
r
(
α
)
r(\alpha)
r(α) 表示角度
α
\alpha
α 对应的距离。
从这个积分我们得到一个面积的估计
−
g
′
≈
1
2
(
r
2
+
(
r
′
′
)
2
)
(
θ
′
′
−
θ
)
2
-g'\approx\frac{1}{2}\frac{(r^2+(r'')^2)(\theta''-\theta)}{2}
−g′≈212(r2+(r′′)2)(θ′′−θ)
其中
c
o
s
(
θ
−
θ
′
′
)
=
r
⋅
r
′
′
∣
∣
r
∣
∣
∣
∣
r
′
′
∣
∣
cos(\theta-\theta'')=\frac{r\cdot r''}{||r||||r''||}
cos(θ−θ′′)=∣∣r∣∣∣∣r′′∣∣r⋅r′′
这样我们就可以估计面积从而估计k。
计算a和e
我们先写下
r
=
k
1
+
e
cos
θ
and
r
′
′
=
k
1
+
e
cos
θ
′
′
r=\frac{k}{1+e\cos \theta} \text{ and } r''=\frac{k}{1+e\cos \theta''}
r=1+ecosθk and r′′=1+ecosθ′′k
重写成
e
cos
θ
=
k
r
−
1
and
e
cos
θ
′
′
=
k
r
′
′
−
1
e\cos \theta=\frac{k}{r}-1 \text{ and } e \cos\theta''=\frac{k}{r''}-1
ecosθ=rk−1 and ecosθ′′=r′′k−1
上面两式相除得到
k
r
′
′
−
1
k
r
−
1
=
cos
θ
′
′
cos
θ
=
cos
(
θ
+
(
θ
′
′
−
θ
)
)
cos
θ
=
cos
(
θ
′
′
−
θ
)
−
tan
θ
sin
(
θ
′
′
−
θ
)
\frac{\frac{k}{r''}-1}{\frac{k}{r}-1}=\frac{\cos \theta''}{\cos \theta}=\frac{\cos (\theta+(\theta''-\theta))}{\cos\theta}=\cos(\theta''-\theta)-\tan\theta\sin(\theta''-\theta)
rk−1r′′k−1=cosθcosθ′′=cosθcos(θ+(θ′′−θ))=cos(θ′′−θ)−tanθsin(θ′′−θ)
因为我们已经知道
θ
′
′
−
θ
\theta''-\theta
θ′′−θ,因此唯一不知道的就是
tan
θ
\tan\theta
tanθ。所以我们可以从已知的
k
,
r
,
r
′
′
,
sin
(
θ
′
′
−
θ
)
,
cos
(
θ
′
′
−
θ
)
k,r,r'',\sin(\theta''-\theta),\cos(\theta''-\theta)
k,r,r′′,sin(θ′′−θ),cos(θ′′−θ)得到
tan
θ
\tan\theta
tanθ。考虑角度的象限后可以得到
θ
\theta
θ。然后我们可以用前面的方法得到
ω
\omega
ω. 然后用极坐标公式可以得到 e,最后用k和e得到a。
最后一个变量 τ p \tau_p τp
这里选择一个特殊点的时刻作为
τ
p
\tau_p
τp,即轨道的perihelion(近日点)。我们需要引入一个新的定义和更多的几何结论。首先定义轨道上一点的 eccentric anomaly (离心偏角)记为E。
假设有一个半径为a的圆,圆心和椭圆中心重合。我们已经定义了true anomaly(位置的真实偏角)
θ
\theta
θ。P点到主轴的垂线与圆交于点Q,Q点与圆心的夹角就是E。
然后我们需要一些与E相关的结论来计算
τ
p
\tau_p
τp。
首先我们推导E和
θ
\theta
θ的关系。
引理 偏心率和
θ
\theta
θ满足
tan
θ
2
=
1
+
e
1
−
e
tan
E
2
\tan \frac{\theta}{2}=\sqrt{\frac{1+e}{1-e}}\tan \frac{E}{2}
tan2θ=1−e1+etan2E
证明:以
θ
\theta
θ在第一象限为例。其他象限可以利用椭圆的对称性证明。我们利用半角公式
tan
2
α
2
=
1
−
cos
α
1
+
cos
α
\tan^2 \frac{\alpha}{2}=\frac{1-\cos \alpha}{1+\cos \alpha}
tan22α=1+cosα1−cosα
根据图1.11,有以下等式
∣
C
R
∣
=
∣
C
F
∣
+
∣
F
R
∣
∣
C
R
∣
=
a
cos
E
∣
C
F
∣
=
a
e
∣
F
R
∣
=
r
cos
θ
a
cos
E
=
a
e
+
r
cos
θ
r
=
a
(
1
−
e
2
)
1
+
e
cos
θ
a
cos
E
=
a
e
+
a
(
1
−
e
2
)
1
+
e
cos
θ
cos
θ
cos
E
=
e
+
1
−
e
2
1
+
e
cos
θ
cos
θ
\begin{array}{l} | CR | = |CF| + |FR|\\ | CR | = a \cos E\\ | CF | = a e\\ | FR | = r \cos \theta\\ a \cos E = a e + r \cos \theta\\ r = \frac{a (1 - e^2) }{1 + e \cos \theta} \\ a \cos E = a e + \frac{a (1 - e^2) }{1 + e \cos \theta} \cos \theta\\ \cos E = e + \frac{1 - e^2}{1 + e \cos \theta} \cos \theta \end{array}
∣CR∣=∣CF∣+∣FR∣∣CR∣=acosE∣CF∣=ae∣FR∣=rcosθacosE=ae+rcosθr=1+ecosθa(1−e2)acosE=ae+1+ecosθa(1−e2)cosθcosE=e+1+ecosθ1−e2cosθ
代入半角公式则有
tan
2
E
2
=
1
−
(
e
+
1
−
e
2
1
+
e
cos
θ
cos
θ
)
1
+
(
e
+
1
−
e
2
1
+
e
cos
θ
cos
θ
)
=
(
1
−
e
)
(
1
+
e
cos
θ
)
−
(
1
−
e
2
)
cos
θ
(
1
+
e
)
(
1
+
e
cos
θ
)
+
(
1
−
e
2
)
cos
θ
=
(
1
−
e
)
(
1
+
e
cos
θ
−
(
1
+
e
)
cos
θ
)
(
1
+
e
)
(
1
+
e
cos
θ
+
(
1
−
e
)
cos
θ
)
=
1
−
e
1
+
e
⋅
1
−
cos
θ
1
+
cos
θ
=
1
−
e
1
+
e
⋅
tan
2
θ
2
\begin{array}{lll} \tan^2 \frac{E}{2} &= \frac{1 - \left( e + \frac{1 - e^2}{1 + e \cos \theta} \cos \theta \right)}{1 + \left( e + \frac{1 - e^2}{1 + e \cos \theta} \cos \theta \right)} & = \frac{ (1 - e) (1 + e \cos \theta) - (1 - e^2) \cos \theta}{(1 + e) (1 + e \cos \theta) + (1 - e^2) \cos \theta}\\ && \\ &= \frac{(1 - e) (1 + e \cos \theta - (1 + e) \cos \theta)}{(1 + e) (1 + e \cos \theta + (1 - e) \cos \theta)} & = \frac{1 - e}{1 + e} \cdot \frac{1 - \cos \theta}{1 + \cos \theta}\\ && \\ &=\frac{1-e}{1+e}\cdot \tan^2 \frac{\theta}{2} & \end{array}
tan22E=1+(e+1+ecosθ1−e2cosθ)1−(e+1+ecosθ1−e2cosθ)=(1+e)(1+ecosθ+(1−e)cosθ)(1−e)(1+ecosθ−(1+e)cosθ)=1+e1−e⋅tan22θ=(1+e)(1+ecosθ)+(1−e2)cosθ(1−e)(1+ecosθ)−(1−e2)cosθ=1+e1−e⋅1+cosθ1−cosθ
两角都位于第一象限,符号相同,开方后就得到我们的引理。
□
\square
□
因为前面我们已经确定了
θ
\theta
θ 和 e,所以我们根据这个引理就可以得到时刻
τ
\tau
τ 对应的位置
r
\bold{r}
r 对应的离心角 E 。下面的引理给出我们计算最后一个参数的方式:E 和
τ
p
\tau_p
τp 的关系。
引理 (Kepler’s equation)在时刻
τ
\tau
τ,我们有
E
−
e
sin
E
=
2
π
t
p
(
τ
−
τ
p
)
E-e\sin E=\frac{2\pi}{t_p}(\tau-\tau_p)
E−esinE=tp2π(τ−τp)
其中 E 是
τ
\tau
τ 时刻的偏心角弧度值,
t
p
t_p
tp 是行星运动的周期,
τ
p
\tau_p
τp 是行星最后一次处于升交点的时刻。
证明:令S表示行星在
τ
p
\tau_p
τp和
τ
\tau
τ之间扫过的面积。椭圆的面积是
π
a
b
\pi ab
πab。开普勒第二定律说
S
τ
−
τ
p
=
π
a
b
t
p
\frac{S}{\tau-\tau_p}=\frac{\pi ab}{t_p}
τ−τpS=tpπab
我们需要用E来表示S。下面的推导来自 Roy 2020 因为比较初级且有几何意义。我们考虑
θ
\theta
θ 小于180度的情况。椭圆的对称性可以给我们180度到360度之间的同样的结果。
考虑图1.12。我们暂时引入uv直角坐标系,来利用代数方程
u
2
a
2
+
v
2
b
2
=
1
\frac{u^2}{a^2}+\frac{v^2}{b^2}=1
a2u2+b2v2=1
外侧的圆的方程是
u
2
+
v
2
=
a
2
u^2+v^2=a^2
u2+v2=a2. 考虑上半平面,我们可以用u表达椭圆的纵坐标v
v
e
(
u
)
=
b
1
−
u
2
a
2
=
b
a
a
2
−
u
2
v_e(u)=b\sqrt{1-\frac{u^2}{a^2}}=\frac{b}{a}\sqrt{a^2-u^2}
ve(u)=b1−a2u2=aba2−u2
以及圆的纵坐标
v
c
(
u
)
=
a
2
−
u
2
v_c(u)=\sqrt{a^2-u^2}
vc(u)=a2−u2
对所有
u
∈
[
−
a
,
a
]
u\in [-a,a]
u∈[−a,a]。
因此对每个u,我们得到
v
e
(
u
)
=
b
a
v
c
(
u
)
v_e(u)=\frac{b}{a}v_c(u)
ve(u)=abvc(u)
因此从图1.12中我们可以观察到三个有用的关系
- ∣ R P ∣ = ( b / a ) ∣ R Q ∣ |RP|=(b/a)|RQ| ∣RP∣=(b/a)∣RQ∣。对应的 △ F R P \triangle FRP △FRP的面积是 △ F R Q \triangle FRQ △FRQ 面积的 b/a.
- 椭圆区域 R T P RTP RTP的面积是圆形中区域 R T Q RTQ RTQ 面积的 b/a.
- 椭圆区域 F T P FTP FTP的面积,也就是S,是圆形中区域 F T Q FTQ FTQ 面积的 b/a.
第2点可以由积分的线性性得到,我们把RTP的面积用积分表示然后就有
∫
R
1
a
v
e
(
u
)
d
u
=
∫
R
1
a
b
a
v
c
(
u
)
d
u
=
b
a
∫
R
1
a
v
c
(
u
)
d
u
\int_{R_1}^{a}v_e(u)du=\int_{R_1}^{a}\frac{b}{a}v_c(u)du=\frac{b}{a}\int_{R_1}^{a}v_c(u)du
∫R1ave(u)du=∫R1aabvc(u)du=ab∫R1avc(u)du
用第1和第2点就可以得到第3点结论。
现在为了计算S,我们只需要计算圆形中的面积 FTQ然后缩小 b/a 即可。
我们计算FTQ的方法是从扇形CTQ 中减去三角形 CFQ。扇形CTQ的面积可以用E占据圆形的比例表示,即
E
2
π
π
a
2
=
E
a
2
2
\frac{E}{2\pi}\pi a^2=\frac{Ea^2}{2}
2πEπa2=2Ea2
三角形CFQ的面积可以用底ae乘以高
a
sin
E
a\sin E
asinE得到,即
a
2
e
sin
E
2
\frac{a^2e\sin E}{2}
2a2esinE. 两者相减就得到FTQ的面积
a
2
(
E
−
a
sin
E
)
2
\frac{a^2(E-a\sin E)}{2}
2a2(E−asinE)
然后我们就得到
S
=
a
b
(
E
−
e
sin
E
)
2
S=\frac{ab(E-e\sin E)}{2}
S=2ab(E−esinE)
再由开普勒第二定律,
S
τ
−
τ
p
=
π
a
b
t
p
\frac{S}{\tau-\tau_p}=\frac{\pi ab}{t_p}
τ−τpS=tpπab
我们就得到
E
−
e
sin
E
2
(
τ
−
τ
p
)
=
π
t
p
\frac{E-e\sin E}{2(\tau-\tau_p)}=\frac{\pi}{t_p}
2(τ−τp)E−esinE=tpπ
这样就证明了引理。
高斯的方法
上一部分主要还是根据开普勒三定律和椭圆的性质可以得到的跟椭圆轨道参数有关的结论。下面介绍高斯是怎么基于观测数据中的 3 \bold{3} 3个点来计算出行星轨道的参数的。
首先在记号上,尽量沿用高斯使用的记号,用 τ , τ ′ , τ ′ ′ \tau,\tau',\tau'' τ,τ′,τ′′来表示Piazzi观测的三个数据点的时刻。各个参数针对 τ \tau τ表示,相应的也有 τ ′ , τ ′ ′ \tau',\tau'' τ′,τ′′的版本。
现在我们开始建立必要的坐标系并描述他们之间的关系。在前面我们已经定义了日心为原点的笛卡尔坐标系,Ceres的位置 r \bold{r} r表示为 ( x , y , z ) T (x,y,z)^T (x,y,z)T。这里我们把地球的日心坐标位置记为 ( X , Y , Z ) T (X,Y,Z)^T (X,Y,Z)T。
令 ξ η ζ \xi\eta\zeta ξηζ表示地心笛卡尔系统。这个坐标系的 ξ , η , ζ \xi, \eta, \zeta ξ,η,ζ-axis分别和日心坐标系的x,y,z-axis平行且同向。长度以天文单位表示。Ceres的位置记为 ( ξ , η , ζ ) T (\xi,\eta,\zeta)^T (ξ,η,ζ)T。
定义完笛卡尔坐标后,我们转到球坐标。记L为日心系地球的经度,B为日心系地球的纬度。这意味着L是日心出发的x轴沿逆时针方向到地球在xy平面的投影位置的角度。B是从日心观测的地球离xy平面的角度。最后记R为地球到太阳的距离。记 D = R cos B D=R\cos B D=RcosB 为R投影到xy平面上的长度。注意我们已经选择xy平面为地球的轨道平面,所以已经有 B=0,D=R。不过接下来我们仍然会使用B和D来保持记号的一致性。这里已经定义的参数在任何时刻都可以认为是已知的因为地球轨道的信息是相当完善的。
然后我们记
λ
\lambda
λ 和
β
\beta
β 分别为 Ceres 的地心系经度和纬度。经度
λ
\lambda
λ 是
ξ
\xi
ξ-axis 和 Ceres 在
ξ
η
\xi\eta
ξη 平面投影的逆时针角度。纬度
β
\beta
β 是 Ceres 到
ξ
η
\xi\eta
ξη 平面的角度。这两个角度是 Piazzi 在他的 Sicilian observatory 观测到的数据。最后,记
ρ
\rho
ρ 为地球到 Ceres 的距离,以及
δ
=
ρ
cos
β
\delta=\rho \cos \beta
δ=ρcosβ 为在
ξ
η
\xi\eta
ξη 平面上投影的长度。这两个距离是未知的,这里主要的目标就是求得他们。注意这里
λ
\lambda
λ 和
L
L
L 的值可以从 0 到 360 度。
β
\beta
β 可以从 -90 到 90 度变化。看图2.1。
再次回顾一下,
r
\bold{r}
r 表示 Ceres 的日心系位置向量,r 是其模长。
笛卡尔系和球坐标系之间的关系是直接的。
ξ
=
ρ
cos
λ
cos
β
η
=
ρ
sin
λ
cos
β
ζ
=
ρ
sin
β
\begin{array}{l} \xi= \rho \cos \lambda \cos \beta\\ \eta= \rho \sin \lambda \cos \beta\\ \zeta= \rho \sin \beta \end{array}
ξ=ρcosλcosβη=ρsinλcosβζ=ρsinβ
地球的两个坐标系之间的转换也是一样的。另外,两个笛卡尔系之间有简单的关系
ξ
=
x
−
X
η
=
y
−
Y
ζ
=
z
−
Z
\xi = x-X\\ \eta = y-Y\\ \zeta=z-Z
ξ=x−Xη=y−Yζ=z−Z
这个方法要求我们处理行星运动过程中扫过的面积。记
f
f
f 表示太阳和 Ceres 的时刻
τ
′
\tau'
τ′ 和
τ
′
′
\tau''
τ′′ 的两个位置之间形成的三角形的面积(看图2.2)。类似的,令
−
f
′
-f'
−f′ 表示
τ
\tau
τ 和
τ
′
′
\tau''
τ′′ 时刻的三角形的面积。令
f
′
′
f''
f′′ 表示时刻
τ
\tau
τ 和
τ
′
\tau'
τ′ 的三角形面积。这里选择定义
−
f
′
-f'
−f′ 是为了让
f
+
f
′
+
f
′
′
f+f'+f''
f+f′+f′′ 表示 Ceres 在三个时刻组成的小三角形的面积。为了下面用矩阵形式表示,我们记向量
f
=
(
f
,
f
′
,
f
′
′
)
T
\bold{f}=(f,f',f'')^T
f=(f,f′,f′′)T。
类似地,记地球的三个对应的面积为
F
,
F
′
,
F
′
′
F,F',F''
F,F′,F′′,以及
F
=
(
F
,
F
′
,
F
′
′
)
T
\bold{F}=(F,F',F'')^T
F=(F,F′,F′′)T。
然后我们考虑椭圆区域的面积而非三角形的面积。记 Ceres 对应的时间扫过的面积分别为
g
,
−
g
′
,
g
′
′
g,-g',g''
g,−g′,g′′(见图2.3)。
现在我们开始推导必要的方程。
根据开普勒第一定律,行星的轨道是一个椭圆而且太阳位于其中一个焦点。特别地,三个向量
r
,
r
′
,
r
′
′
\bold{r},\bold{r}',\bold{r}''
r,r′,r′′位于一个平面上,可以由它们中的任意两个确定该平面。因此存在两个常数
c
1
,
c
2
c1,c2
c1,c2 使得
r
=
c
1
r
′
+
c
2
r
′
′
(2.1)
\bold{r}=c_1\bold{r}'+c_2\bold{r}'' \tag{2.1}
r=c1r′+c2r′′(2.1)
我们用
r
′
\bold{r}'
r′ 和
r
′
′
\bold{r}''
r′′ 在两边做叉乘可以得到
r
×
r
′
=
c
2
(
r
′
′
×
r
′
)
and
r
×
r
′
′
=
c
1
(
r
′
×
r
′
′
)
(2.2)
\bold{r} \times\bold{r}' = c_2 (r'' \times r') \text{ and } \bold{r} \times\bold{r}'' = c_1 (r' \times r'') \tag{2.2}
r×r′=c2(r′′×r′) and r×r′′=c1(r′×r′′)(2.2)
叉乘有一个很好的性质:
∣
∣
a
×
b
∣
∣
||\bold{a}\times \bold{b}||
∣∣a×b∣∣ 等于向量
a
,
b
\bold{a},\bold{b}
a,b 张成的平行四边形的面积。因此,
2
f
2f
2f 是
r
′
,
r
′
′
\bold{r}',\bold{r}''
r′,r′′ 张成的面积。所以我们有这些关系
2
f
=
∣
∣
r
′
×
r
′
′
∣
∣
2f= ||\bold{r}'\times \bold{r}''||
2f=∣∣r′×r′′∣∣,
−
2
f
′
=
∣
∣
r
×
r
′
′
∣
∣
-2f'= ||\bold{r}\times \bold{r}''||
−2f′=∣∣r×r′′∣∣ ,
2
f
′
′
=
∣
∣
r
×
r
′
′
∣
∣
2f''=||\bold{r}\times\bold{r}''||
2f′′=∣∣r×r′′∣∣。
借助前面两式,我们可以得到
f
′
′
=
c
2
f
f''=c_2f
f′′=c2f 以及
−
f
′
=
c
1
f
-f'=c_1f
−f′=c1f。代入前面我们可以得到
f
r
=
−
f
′
r
′
−
f
′
′
r
′
′
也就是
f
r
+
f
′
r
′
+
f
′
′
r
′
′
=
0
(2.3)
f\bold{r}=-f'\bold{r}'-f''\bold{r}''\text{ 也就是 } f\bold{r}+f'\bold{r}'+f''\bold{r}''=0 \tag{2.3}
fr=−f′r′−f′′r′′ 也就是 fr+f′r′+f′′r′′=0(2.3)
这个式子可以用矩阵乘法来表达。相关的矩阵定义为
ϕ
=
(
x
x
′
x
′
′
y
y
′
y
′
′
z
z
′
z
′
′
)
=
(
r
∣
r
′
∣
r
′
′
)
and
Φ
=
(
X
X
′
X
′
′
Y
Y
′
Y
′
′
Z
Z
′
Z
′
′
)
\bold{\phi}= (\begin{array}{c} x \quad x' \quad x''\\ y \quad y' \quad y''\\ z \quad z' \quad z'' \end{array} ) = (\bold{r} | \bold{r}' | \bold{r}'') \text{ and } \bold{\Phi}= \left(\begin{array}{c} X \quad X' \quad X''\\ Y \quad Y' \quad Y''\\ Z \quad Z' \quad Z'' \end{array}\right)
ϕ=(xx′x′′yy′y′′zz′z′′)=(r∣r′∣r′′) and Φ=
XX′X′′YY′Y′′ZZ′Z′′
这样前式可以写成
ϕ
f
=
0
\phi \bold{f}=\bold{0}
ϕf=0。类似地我们可以得到另一个重要的方程
Φ
F
=
0
\Phi\bold{F}=\bold{0}
ΦF=0。从
ϕ
f
=
0
\phi\bold{f}=\bold{0}
ϕf=0 减去
Φ
F
\Phi\bold{F}
ΦF 再乘以
(
F
+
F
′
′
)
(F+F'')
(F+F′′) 得到
(
F
+
F
′
′
)
(
ϕ
f
−
Φ
f
)
=
(
F
+
F
′
′
)
(
−
Φ
f
)
(F+F'')(\phi\bold{f}-\Phi\bold{f})=(F+F'')(-\Phi\bold{f})
(F+F′′)(ϕf−Φf)=(F+F′′)(−Φf)
因此我们有
(
F
+
F
′
′
)
(
ϕ
−
Φ
)
f
=
−
Φ
(
F
+
F
′
′
)
f
(2.4)
(F+F'')(\phi-\Phi)\bold{f}=-\Phi(F+F'')\bold{f} \tag{2.4}
(F+F′′)(ϕ−Φ)f=−Φ(F+F′′)f(2.4)
因为
Φ
F
=
0
\Phi\bold{F}=\bold{0}
ΦF=0,所以也有
Φ
(
f
+
f
′
′
)
F
=
(
f
+
f
′
′
)
Φ
F
=
0
\Phi(f+f'')\bold{F}=(f+f'')\Phi\bold{F}=\bold{0}
Φ(f+f′′)F=(f+f′′)ΦF=0
因此我们可以把这一个零项加到 (2.4) 式的右边得到
(
F
+
F
′
′
)
(
ϕ
−
Φ
)
f
=
−
Φ
(
F
+
F
′
′
)
f
+
Φ
(
f
+
f
′
′
)
F
(F+F'')(\phi-\Phi)\bold{f}=-\Phi(F+F'')\bold{f}+\Phi(f+f'')\bold{F}
(F+F′′)(ϕ−Φ)f=−Φ(F+F′′)f+Φ(f+f′′)F
把
Φ
\Phi
Φ 因子提出就得到一个关键的方程
(
F
+
F
′
′
)
(
ϕ
−
Φ
)
f
=
Φ
(
(
f
+
f
′
′
)
F
−
(
F
+
F
′
′
)
f
)
(2.5)
(F+F'')(\phi-\Phi)\bold{f}=\Phi((f+f'')\bold{F}-(F+F'')\bold{f}) \tag{2.5}
(F+F′′)(ϕ−Φ)f=Φ((f+f′′)F−(F+F′′)f)(2.5)
下一步,我们把(2.5)式从笛卡尔坐标系转到球坐标系。可以定义以下有用的向量
w
=
(
cos
λ
,
sin
λ
,
tan
β
)
T
and
W
=
(
cos
L
,
sin
L
,
tan
β
)
T
\begin{array}{l} \bold{w} = (\cos \lambda, \sin \lambda, \tan \beta)^T \text{ and } \bold{W} = \end{array} (\cos L, \sin L, \tan \beta)^T
w=(cosλ,sinλ,tanβ)T and W=(cosL,sinL,tanβ)T
注意我们同时也定义了
w
′
,
w
′
′
\bold{w}',\bold{w}''
w′,w′′和
W
′
,
W
′
′
\bold{W}',\bold{W}''
W′,W′′。这里在记号上与Gauss原始的有所不同:Gauss使用的是
π
\pi
π 而不是
w
\bold{w}
w,而我们一直用来表示圆周的常数了。
注意这里
w
\bold{w}
w 和
W
\bold{W}
W 定义中的所有值都是已知的。而且注意
w
\bold{w}
w 乘以未知的
δ
=
ρ
cos
β
\delta=\rho\cos\beta
δ=ρcosβ 得到
δ
w
=
(
ρ
cos
λ
,
ρ
sin
λ
,
ρ
sin
β
)
T
=
(
ξ
,
η
,
ζ
)
T
\delta\bold{w}=(\rho\cos\lambda,\rho\sin\lambda,\rho\sin\beta)^T=(\xi,\eta,\zeta)^T
δw=(ρcosλ,ρsinλ,ρsinβ)T=(ξ,η,ζ)T
类似地,我们对地球的位置有
D
W
=
(
X
,
Y
,
Z
)
T
D\bold{W}=(X,Y,Z)^T
DW=(X,Y,Z)T。对于时刻
τ
′
,
τ
′
′
\tau',\tau''
τ′,τ′′ 也是一样。
再回忆地心和日心笛卡尔坐标系之间的简单关系
ξ
=
x
−
X
,
η
=
y
−
Y
,
ζ
=
z
−
Z
\xi=x-X,\eta=y-Y,\zeta=z-Z
ξ=x−X,η=y−Y,ζ=z−Z,以及对
τ
′
,
τ
′
′
\tau',\tau''
τ′,τ′′ 类似的关系。转化成矩阵表示就是
ϕ
−
Φ
=
(
ξ
ξ
′
ξ
′
′
η
η
′
η
′
′
ζ
ζ
′
ζ
′
′
)
\phi - \Phi = \left(\begin{array}{c} \xi \quad \xi' \quad \xi''\\ \eta \quad \eta' \quad \eta''\\ \zeta \quad \zeta' \quad \zeta'' \end{array}\right)
ϕ−Φ=
ξξ′ξ′′ηη′η′′ζζ′ζ′′
因此我们得到
ϕ
−
Φ
=
(
δ
w
∣
δ
′
w
′
∣
δ
′
′
w
′
′
)
\phi-\Phi=(\delta\bold{w}|\delta'\bold{w}'|\delta''\bold{w}'')
ϕ−Φ=(δw∣δ′w′∣δ′′w′′) 以及
Φ
=
(
D
W
∣
D
′
W
′
∣
D
′
′
W
′
′
)
\Phi=(D\bold{W}|D'\bold{W}'|D''\bold{W}'')
Φ=(DW∣D′W′∣D′′W′′)。利用这些关系以及三元恒等式
(
a
×
b
)
⋅
c
=
det
(
a
∣
b
∣
c
)
(\bold{a}\times\bold{b})\cdot\bold{c}=\det(\bold{a}|\bold{b}|\bold{c})
(a×b)⋅c=det(a∣b∣c),允许我们把式(2.5)转成三个新的方程,然后成为我们进行估计的开始,并且之后让我们可以估计未知数。我们也会利用一些行列式的性质,
det
(
a
∣
b
∣
b
)
=
0
\det(\bold{a}|\bold{b}|\bold{b})=0
det(a∣b∣b)=0 以及列的偶置换下行列式的值不变。
来推导第一个方程,我们在式(2.5)的两边和
(
w
′
′
×
w
)
(\bold{w}''\times\bold{w})
(w′′×w) 求内积。左边得到
(
w
′
′
×
w
)
⋅
(
F
+
F
′
′
)
(
ϕ
−
Φ
)
f
=
(
F
+
F
′
′
)
(
w
′
′
×
w
)
⋅
(
ϕ
−
Φ
)
f
=
(
F
+
F
′
′
)
det
(
w
′
′
∣
w
∣
(
ϕ
−
Φ
)
f
)
=
(
F
+
F
′
′
)
det
(
w
′
′
∣
w
∣
f
δ
w
+
f
′
δ
′
w
′
+
f
′
′
δ
′
′
w
′
′
)
=
(
F
+
F
′
′
)
(
det
(
w
′
′
∣
w
∣
f
δ
w
)
+
det
(
w
′
′
∣
w
∣
f
′
δ
′
w
′
)
+
det
(
w
′
′
∣
w
∣
f
′
′
δ
′
′
w
′
′
)
)
=
(
F
+
F
′
′
)
(
f
δ
det
(
w
′
′
∣
w
∣
w
)
+
f
′
δ
′
det
(
w
′
′
∣
w
∣
w
′
)
+
f
′
′
δ
′
′
det
(
w
′
′
∣
w
∣
w
′
′
)
)
=
(
F
+
F
′
′
)
f
′
δ
′
det
(
w
′
′
∣
w
∣
w
′
)
=
(
F
+
F
′
′
)
f
′
δ
′
det
(
w
∣
w
′
∣
w
′
′
)
\begin{array}{l} \begin{array}{c}\\ \end{array} (\bold{w}'' \times \bold{w}) \cdot (F + F'') (\phi - \Phi) \bold{f}\\ = (F + F'') (\bold{w}'' \times \bold{w}) \cdot (\phi - \Phi) \bold{f}\\ = (F + F'') \det (\bold{w}'' | \bold{w} | (\phi - \Phi) \bold{f})\\ = (F + F'') \det (\bold{w}'' | \bold{w} | f \delta \bold{w} + f' \delta' \bold{w}' + f'' \delta'' \bold{w}'')\\ = (F + F'') (\det (\bold{w}'' | \bold{w} | f \delta \bold{w}) + \det (\bold{w}'' | \bold{w} | f' \delta' \bold{w}') + \det (\bold{w}'' | \bold{w} | f'' \delta'' \bold{w}''))\\ = (F + F'') (f \delta \det (\bold{w}'' | \bold{w} | \bold{w}) + f' \delta' \det (\bold{w}'' | \bold{w} | \bold{w}') + f'' \delta'' \det (\bold{w}'' | \bold{w} | \bold{w}''))\\ = (F + F'') f' \delta' \det (\bold{w}'' | \bold{w} | \bold{w}')\\ = (F + F'') f' \delta' \det (\bold{w} | \bold{w}' | \bold{w}'') \end{array}
(w′′×w)⋅(F+F′′)(ϕ−Φ)f=(F+F′′)(w′′×w)⋅(ϕ−Φ)f=(F+F′′)det(w′′∣w∣(ϕ−Φ)f)=(F+F′′)det(w′′∣w∣fδw+f′δ′w′+f′′δ′′w′′)=(F+F′′)(det(w′′∣w∣fδw)+det(w′′∣w∣f′δ′w′)+det(w′′∣w∣f′′δ′′w′′))=(F+F′′)(fδdet(w′′∣w∣w)+f′δ′det(w′′∣w∣w′)+f′′δ′′det(w′′∣w∣w′′))=(F+F′′)f′δ′det(w′′∣w∣w′)=(F+F′′)f′δ′det(w∣w′∣w′′)
而等式右边是
(
w
′
′
×
w
)
⋅
Φ
(
(
f
+
f
′
′
)
F
−
(
F
+
F
′
′
)
f
)
=
(
f
+
f
′
′
)
(
w
′
′
×
w
)
⋅
Φ
F
−
(
F
+
F
′
′
)
(
w
′
′
×
w
)
⋅
Φ
f
=
(
f
+
f
′
′
)
det
(
w
′
′
∣
w
∣
Φ
F
)
−
(
F
+
F
′
′
)
det
(
w
′
′
∣
w
∣
Φ
f
)
=
(
f
+
f
′
′
)
det
(
w
′
′
∣
w
∣
F
D
W
+
F
′
D
′
W
′
+
F
′
′
D
′
′
W
′
′
)
−
(
F
+
F
′
′
)
det
(
w
′
′
∣
w
∣
f
D
W
+
f
′
D
′
W
′
+
f
′
′
D
′
′
W
′
′
)
=
(
f
+
f
′
′
)
(
F
D
det
(
w
′
′
∣
w
∣
W
)
+
F
′
D
′
det
(
w
′
′
∣
w
∣
W
′
)
+
F
′
′
D
′
′
det
(
w
′
′
∣
w
∣
W
′
′
)
)
−
(
F
+
F
′
′
)
(
f
D
det
(
w
′
′
∣
w
∣
W
)
+
f
′
D
′
det
(
w
′
′
∣
w
∣
W
′
)
+
f
′
′
D
′
′
det
(
w
′
′
∣
w
∣
W
′
′
)
)
=
(
f
′
′
F
−
f
F
′
′
)
(
D
det
(
w
′
′
∣
w
∣
W
)
−
D
′
′
det
(
w
′
′
∣
w
∣
W
′
′
)
)
+
(
(
f
+
f
′
′
)
F
′
−
(
F
+
F
′
′
)
f
′
)
D
′
det
(
w
′
′
∣
w
∣
W
′
)
=
(
f
′
′
F
−
f
F
′
′
)
(
D
det
(
w
∣
W
∣
w
′
′
)
−
D
′
′
det
(
w
∣
W
′
′
∣
w
′
′
)
)
+
(
(
f
+
f
′
′
)
F
′
−
(
F
+
F
′
′
)
f
′
)
D
′
det
(
w
∣
W
′
∣
w
′
′
)
\begin{array}{l} (\bold{w}'' \times \bold{w}) \cdot \Phi \left( (f + f'') \bold{F} - \left( F + F'' \right) \bold{f} \right)\\ = (f + f'') (\bold{w}'' \times \bold{w}) \cdot \Phi \bold{F} - (F + F'') (\bold{w}'' \times \bold{w}) \cdot \Phi \bold{f}\\ = (f + f'') \det (\bold{w}'' | \bold{w} | \Phi F) - (F + F'') \det (\bold{w}'' | \bold{w} | \Phi \bold{f})\\ = (f + f'') \det (\bold{w}'' | \bold{w} |FD \bold{W} + F' D' \bold{W}' + F'' D'' \bold{W}'') \\ \space - (F + F'') \det (\bold{w}'' | \bold{w} | fD \bold{W} + f' D' W' + f'' D'' W'')\\ = (f + f'') (FD\det (\bold{w}'' | \bold{w} | \bold{W}) + F' D' \det (\bold{w}'' | \bold{w} | \bold{W}') + F'' D'' \det (\bold{w}'' | \bold{w} | \bold{W}'')) \\ \space - (F + F'') (fD\det (\bold{w}'' | \bold{w} | \bold{W}) + f' D' \det (\bold{w}'' | \bold{w} | \bold{W}') + f'' D'' \det (\bold{w}'' | \bold{w} | \bold{W}''))\\ = (f'' F - f F'') (D\det (\bold{w}'' | \bold{w} | \bold{W}) - D'' \det (\bold{w}'' | \bold{w} | \bold{W}''))\\ \space + ((f + f'') F' - (F + F'') f') D' \det (\bold{w}'' | \bold{w} | \bold{W}')\\ = (f'' F - f F'') (D\det (\bold{w} \bold{| \bold{W} | \bold{w}''}) - D'' \det (\bold{w} | \bold{W''} | \bold{w}'')) \\ \space + ((f + f'') F' - (F + F'') f') D' \det (\bold{w} | \bold{W'} | \bold{w}'') \end{array}
(w′′×w)⋅Φ((f+f′′)F−(F+F′′)f)=(f+f′′)(w′′×w)⋅ΦF−(F+F′′)(w′′×w)⋅Φf=(f+f′′)det(w′′∣w∣ΦF)−(F+F′′)det(w′′∣w∣Φf)=(f+f′′)det(w′′∣w∣FDW+F′D′W′+F′′D′′W′′) −(F+F′′)det(w′′∣w∣fDW+f′D′W′+f′′D′′W′′)=(f+f′′)(FDdet(w′′∣w∣W)+F′D′det(w′′∣w∣W′)+F′′D′′det(w′′∣w∣W′′)) −(F+F′′)(fDdet(w′′∣w∣W)+f′D′det(w′′∣w∣W′)+f′′D′′det(w′′∣w∣W′′))=(f′′F−fF′′)(Ddet(w′′∣w∣W)−D′′det(w′′∣w∣W′′)) +((f+f′′)F′−(F+F′′)f′)D′det(w′′∣w∣W′)=(f′′F−fF′′)(Ddet(w∣W∣w′′)−D′′det(w∣W′′∣w′′)) +((f+f′′)F′−(F+F′′)f′)D′det(w∣W′∣w′′)
剩下的两个方程类似可从式(2.5)得到,通过两边和
(
W
′
×
w
)
(\bold{W}'\times\bold{w})
(W′×w) 以及
(
W
′
×
w
′
′
)
(\bold{W}'\times\bold{w}'')
(W′×w′′) 求内积得到。为了节约空间免去大量类似的计算,我们略过详细的推导,把结果写在这里。我们共可得到三个等式
(
F
+
F
′
′
)
f
′
δ
′
det
(
w
∣
w
′
∣
w
′
′
)
=
(
f
′
′
F
−
f
F
′
′
)
(
D
det
(
w
∣
W
∣
w
′
′
)
−
D
′
′
det
(
w
∣
W
′
′
∣
w
′
′
)
)
+
(
(
f
+
f
′
′
)
F
′
−
(
F
+
F
′
′
)
f
′
)
D
′
det
(
w
∣
W
′
∣
w
′
′
)
(2.6)
\begin{array}{l} (F + F'') f' \delta' \det (\bold{w} | \bold{w}' | \bold{w}'')\\ = (f'' F - f F'') (D\det (\bold{w} | \bold{W} | \bold{w}'') - D'' \det (\bold{w} | \bold{W}'' | \bold{w}''))\\ + ((f + f'') F' - (F + F'') f') D' \det (\bold{w} | \bold{W}' | \bold{w}'') \end{array} \tag{2.6}
(F+F′′)f′δ′det(w∣w′∣w′′)=(f′′F−fF′′)(Ddet(w∣W∣w′′)−D′′det(w∣W′′∣w′′))+((f+f′′)F′−(F+F′′)f′)D′det(w∣W′∣w′′)(2.6)
( F + F ′ ′ ) ( f ′ δ ′ det ( w ∣ w ′ ∣ W ′ ) + f ′ ′ δ ′ ′ det ( w ∣ w ′ ′ ∣ W ′ ) ) = ( f ′ ′ F − f F ′ ′ ) ( D det ( w ∣ W ∣ W ′ ) − D ′ ′ det ( w ∣ W ′ ′ ∣ W ′ ) ) (2.7) \begin{array}{c} (F + F'') (f' \delta' \det (w | w' | \bold{W}') + f'' \delta'' \det (\bold{w} | \bold{w}'' | \bold{W}'))\\ = (f'' F - f F'') (D \det (w | \bold{W} | \bold{W}') - D'' \det (w | \bold{W}'' | \bold{W}')) \end{array} \tag{2.7} (F+F′′)(f′δ′det(w∣w′∣W′)+f′′δ′′det(w∣w′′∣W′))=(f′′F−fF′′)(Ddet(w∣W∣W′)−D′′det(w∣W′′∣W′))(2.7)
( F + F ′ ′ ) ( f δ det ( w ′ ′ ∣ w ∣ W ′ ) + f ′ δ ′ det ( w ′ ′ ∣ w ′ ∣ W ′ ) ) = ( f ′ ′ F − f F ′ ′ ) ( D det ( w ∣ W ∣ W ′ ) − D ′ ′ det ( w ′ ′ ∣ W ′ ′ ∣ W ′ ) ) (2.8) \begin{array}{l} \begin{array}{c} (F + F'') (f \delta \det (\bold{w}'' | \bold{w} | \bold{W}') + f' \delta' \det (\bold{w}'' | \bold{w}' | \bold{W}')) \end{array}\\ = (f'' F - f F'') (D \det (\bold{w} | \bold{W} | \bold{W}') - D'' \det (\bold{w}'' | \bold{W}'' | \bold{W}')) \end{array} \tag{2.8} (F+F′′)(fδdet(w′′∣w∣W′)+f′δ′det(w′′∣w′∣W′))=(f′′F−fF′′)(Ddet(w∣W∣W′)−D′′det(w′′∣W′′∣W′))(2.8)
我们会用(2.7)和(2.8)来基于 δ ′ \delta' δ′ 来估计 δ ′ , δ ′ ′ \delta',\delta'' δ′,δ′′,然后用(2.6)来估计 δ \delta δ 自身。但在那之前,我们先来推导近似版本。
调整估计
在这部分,主要是引入一些(高阶小量)近似来化简前面的方程。这些结果在Gauss 1809 和 Teets & Whitehead 1999 中都有。这里作者自己推导了一遍。
由于Ceres在Piazzi的观测过程中只移动了很小的一个弧线(不到完整轨道的1%),Gauss注意到我们可以把时间差
τ
′
−
τ
,
τ
′
′
−
τ
′
,
τ
′
′
−
τ
\tau'-\tau,\tau''-\tau',\tau''-\tau
τ′−τ,τ′′−τ′,τ′′−τ视为无穷小量,然后用这一事实来舍去方程中更高阶的足够接近0的那些项。利用现代的极限的记号,我们可以分析这些项在时间差趋于0时的表现。特别地,我们考虑的一种假设的情形是
τ
,
τ
′
′
\tau,\tau''
τ,τ′′同时趋于
τ
′
\tau'
τ′。这里简单地记时间差
t
→
0
t \to 0
t→0.
下面我们考虑的变量记号并没有变,不过要注意现在要把变量看成时间的函数而不是静态的具体的值。下面这些式子是合理的
lim
t
→
0
r
=
lim
t
→
0
r
′
′
=
r
′
or
lim
t
→
0
D
=
lim
t
→
0
D
′
′
=
D
′
\lim_{t \rightarrow 0} \bold{r} = \lim_{t \rightarrow 0} \bold{r}'' = \bold{r'} \text{ or } \lim_{t \rightarrow 0} D = \lim_{t \rightarrow 0} D'' = D'
t→0limr=t→0limr′′=r′ or t→0limD=t→0limD′′=D′
接下来我们要证明多个引理,来估计上面方程中的各项的阶(order),记号为
o
o
o。当说一个函数 h(t)在t趋于0时的阶为n时,表示极限
h
(
t
)
t
n
\frac{h(t)}{t^n}
tnh(t)
有界。这时我们说 h(t) 是
o
(
t
n
)
o(t^n)
o(tn)的。这个定义等价于函数关于 t 的泰勒展开仅包含不小于 n 的阶数的项。这一部分我们反复用到的泰勒展开实际上只有sin和cos函数的
sin
α
=
∑
n
=
0
∞
(
−
1
)
n
α
2
n
+
1
(
2
n
+
1
)
!
=
α
−
α
3
6
+
⋯
cos
α
=
∑
n
=
0
∞
(
−
1
)
n
α
2
n
(
2
n
)
!
=
1
−
α
2
2
+
⋯
\sin \alpha=\sum_{n=0}^{\infty}(-1)^n\frac{\alpha^{2n+1}}{(2n+1)!}=\alpha-\frac{\alpha^3}{6}+\cdots\\ \cos\alpha =\sum_{n=0}^{\infty} (-1)^n\frac{\alpha^{2n}}{(2n)!}=1-\frac{\alpha^2}{2}+\cdots
sinα=n=0∑∞(−1)n(2n+1)!α2n+1=α−6α3+⋯cosα=n=0∑∞(−1)n(2n)!α2n=1−2α2+⋯
接下来我们说明一些引理来做支持我们做近似。
引理4 三角面积 f , − f ′ , f ′ ′ , F , − F ′ , F ′ ′ f,-f',f'',F,-F',F'' f,−f′,f′′,F,−F′,F′′ 都是1阶的。
引理5
f
′
′
F
−
f
F
′
′
f''F-fF''
f′′F−fF′′ 和
(
f
+
f
′
′
)
F
′
−
(
F
+
F
′
′
)
f
′
(f+f'')F'-(F+F'')f'
(f+f′′)F′−(F+F′′)f′ 都是4阶的。
这里用三角函数表示并代入泰勒展开可以验证2阶项被消掉了最小的就是4阶项。
引理6
det
(
w
∣
w
′
∣
w
′
′
)
\det(\bold{w}|\bold{w}'|\bold{w}'')
det(w∣w′∣w′′)是3阶的。
这里的证明用行列式根据最后一行的tan项作展开,然后用和角公式变形,再用sin的泰勒公式展开验证1阶项被消掉。
引理7
det
(
w
∣
w
′
∣
W
′
)
,
det
(
w
∣
w
′
′
∣
W
′
)
,
det
(
w
′
′
∣
w
∣
W
′
)
,
det
(
w
′
′
∣
w
′
∣
W
′
)
,
det
(
w
∣
W
′
∣
w
′
′
)
\det(\bold{w}|\bold{w}'|\bold{W}'), \det(\bold{w}|\bold{w}''|\bold{W}'), \det(\bold{w}''|\bold{w}|\bold{W}'), \det(\bold{w}''|\bold{w}'|\bold{W}'), \det(\bold{w}|\bold{W}'|\bold{w}'')
det(w∣w′∣W′),det(w∣w′′∣W′),det(w′′∣w∣W′),det(w′′∣w′∣W′),det(w∣W′∣w′′) 都是1阶的。
证明思路和引理6一样,先根据行列式对最后一行展开,并且利用
B
,
B
′
,
B
′
′
B,B',B''
B,B′,B′′都是0消去一项。剩下的项用泰勒公式展开后可以发现1次项消不掉,因此阶为1.
引理8 下面的三个式子
- D det ( w ∣ W ∣ w ′ ′ ) − D ′ ′ det ( w ∣ W ′ ′ ∣ w ′ ′ ) D\det(\bold{w}|\bold{W}|\bold{w}'')-D''\det(\bold{w}|\bold{W}''|\bold{w}'') Ddet(w∣W∣w′′)−D′′det(w∣W′′∣w′′)
- D det ( w ∣ W ∣ W ′ ) − D ′ ′ det ( w ∣ W ′ ′ ∣ W ′ ) D\det(\bold{w}|\bold{W}|\bold{W}')-D''\det(\bold{w}|\bold{W}''|\bold{W}') Ddet(w∣W∣W′)−D′′det(w∣W′′∣W′)
-
D
det
(
w
∣
′
′
W
∣
W
′
)
−
D
′
′
det
(
w
′
′
∣
W
′
′
∣
W
′
)
D\det(\bold{w}|''\bold{W}|\bold{W}')-D''\det(\bold{w}''|\bold{W}''|\bold{W}')
Ddet(w∣′′W∣W′)−D′′det(w′′∣W′′∣W′)
都是2阶的。
这里也利用 B , B ′ , B ′ ′ B,B',B'' B,B′,B′′为0,消去一些项,剩下的项合并并提取因式可以得到一些两个一阶三角函数的乘积,所以是2阶的。
几个结论
- (2.6)式的左侧 ( F + F ′ ′ ) f ′ δ ′ det ( w ∣ w ′ ∣ w ′ ′ ) (F + F'') f' \delta' \det (\bold{w} | \bold{w}' | \bold{w}'') (F+F′′)f′δ′det(w∣w′∣w′′)是 o ( t 5 ) o(t^5) o(t5),因为 F + F ′ ′ F+F'' F+F′′和 f ′ f' f′分别是1阶, δ ′ \delta' δ′是非零有界的, det ( w ∣ w ′ ∣ w ′ ′ ) \det(\bold{w}|\bold{w}'|\bold{w}'') det(w∣w′∣w′′)根据引理6是3阶的。右侧第一项 ( f ′ ′ F − f F ′ ′ ) ( D det ( w ∣ W ∣ w ′ ′ ) − D ′ ′ det ( w ∣ W ′ ′ ∣ w ′ ′ ) ) (f'' F - f F'') (D\det (\bold{w} | \bold{W} | \bold{w}'') - D'' \det (\bold{w} | \bold{W}'' | \bold{w}'')) (f′′F−fF′′)(Ddet(w∣W∣w′′)−D′′det(w∣W′′∣w′′)) 是6阶的,第二项 ( ( f + f ′ ′ ) F ′ − ( F + F ′ ′ ) f ′ ) D ′ det ( w ∣ W ′ ∣ w ′ ′ ) ((f + f'') F' - (F + F'') f') D' \det (\bold{w} | \bold{W}' | \bold{w}'') ((f+f′′)F′−(F+F′′)f′)D′det(w∣W′∣w′′)是5阶的。所以右边第1项我们可以视为约等于0而忽略。
- (2.7)式的左侧是3阶的,右边是6阶,可以近似为0.
- (2.8)式的左侧是3阶的,右边是6阶,可以近似为0.
初步的结果
在这部分,我们可以根据Gauss 1809和Teets & Whitehead 1999 的工作根据
δ
′
\delta'
δ′估计
δ
,
δ
′
′
\delta,\delta''
δ,δ′′。以此为目标,我们用上节结论中的第2和第3条,来估计式(2.7)和(2.8)。现在我们有
(
F
+
F
′
′
)
(
f
′
δ
′
det
(
w
∣
w
′
∣
W
′
)
+
f
′
′
δ
′
′
det
(
w
∣
w
′
′
∣
W
′
)
)
≈
0
(F + F'') (f' \delta' \det (w | w' | \bold{W}') + f'' \delta'' \det (\bold{w} | \bold{w}'' | \bold{W}')) \approx 0
(F+F′′)(f′δ′det(w∣w′∣W′)+f′′δ′′det(w∣w′′∣W′))≈0
以及
(
F
+
F
′
′
)
(
f
δ
det
(
w
′
′
∣
w
∣
W
′
)
+
f
′
δ
′
det
(
w
′
′
∣
w
′
∣
W
′
)
)
≈
0
(F + F'') (f \delta \det (\bold{w}'' | \bold{w} | \bold{W}') + f' \delta' \det (\bold{w}'' | \bold{w}' | \bold{W}')) \approx 0
(F+F′′)(fδdet(w′′∣w∣W′)+f′δ′det(w′′∣w′∣W′))≈0
我们除以
(
F
+
F
′
′
)
(F+F'')
(F+F′′)然后移项可以得到
f
′
′
δ
′
′
det
(
w
∣
w
′
′
∣
W
′
)
≈
−
f
′
δ
′
det
(
w
∣
w
′
∣
W
′
)
f'' \delta'' \det (\bold{w} | \bold{w}'' | \bold{W}') \approx -f' \delta' \det (w | w' | \bold{W}')
f′′δ′′det(w∣w′′∣W′)≈−f′δ′det(w∣w′∣W′)
以及
f
δ
det
(
w
′
′
∣
w
∣
W
′
)
≈
−
f
′
δ
′
det
(
w
′
′
∣
w
′
∣
W
′
)
f \delta \det (\bold{w}'' | \bold{w} | \bold{W}') \approx -f' \delta' \det (\bold{w}'' | \bold{w}' | \bold{W}')
fδdet(w′′∣w∣W′)≈−f′δ′det(w′′∣w′∣W′)
分离出 delta 的值就有
δ
′
′
≈
−
f
′
det
(
w
∣
w
′
∣
W
′
)
f
′
′
det
(
w
∣
w
′
′
∣
W
′
)
δ
′
and
δ
≈
−
f
′
det
(
w
′
′
∣
w
′
∣
W
′
)
f
det
(
w
′
′
∣
w
∣
W
′
)
δ
′
(2.9)
\delta'' \approx \frac{- f' \det (\bold{w} | \bold{w}' | \bold{W'})}{f'' \det (\bold{w | w'' | W'})} \delta' \space \text{and } \space \delta \approx \frac{- f' \det (\bold{w}'' | \bold{w}' | \bold{W}')}{f \det (\bold{w''} | \bold{w} | \bold{W'})} \delta' \tag{2.9}
δ′′≈f′′det(w∣w′′∣W′)−f′det(w∣w′∣W′)δ′ and δ≈fdet(w′′∣w∣W′)−f′det(w′′∣w′∣W′)δ′(2.9)
现在虽然我们有了估计的表达式,但右侧仍然有一些未知量
f
,
−
f
′
,
f
′
′
f,-f',f''
f,−f′,f′′。不过,他们只以比值的形式出现,这样可以让我们根据开普勒第二定律(恒定时间扫过面积相同)做估计。
g
τ
′
′
−
τ
′
=
−
g
′
τ
′
′
−
τ
=
g
′
′
τ
′
−
τ
′
\frac{g}{\tau''-\tau'}=\frac{-g'}{\tau''-\tau}=\frac{g''}{\tau'-\tau'}
τ′′−τ′g=τ′′−τ−g′=τ′−τ′g′′
也就是
g
′
′
−
g
′
⋅
τ
′
′
−
τ
τ
′
−
τ
=
g
−
g
′
⋅
τ
′
′
−
τ
τ
′
′
−
τ
′
=
1
\frac{g''}{-g'}\cdot \frac{\tau''-\tau}{\tau'-\tau}=\frac{g}{-g'}\cdot \frac{\tau''-\tau}{\tau''-\tau'}=1
−g′g′′⋅τ′−ττ′′−τ=−g′g⋅τ′′−τ′τ′′−τ=1
这样,我们在式(2.9)的左两边乘以1的式子,就可以得到
δ
′
′
≈
f
′
g
′
⋅
g
′
′
f
′
′
⋅
τ
′
′
−
τ
τ
′
−
τ
⋅
det
(
w
∣
w
′
∣
W
′
)
det
(
w
∣
w
′
′
∣
W
′
)
δ
′
\delta''\approx \frac{f'}{g'}\cdot\frac{g''}{f''}\cdot\frac{\tau''-\tau}{\tau'-\tau}\cdot\frac{\det (\bold{w}|\bold{w}'|\bold{W}')}{\det (\bold{w}|\bold{w}''|\bold{W}')}\delta'
δ′′≈g′f′⋅f′′g′′⋅τ′−ττ′′−τ⋅det(w∣w′′∣W′)det(w∣w′∣W′)δ′
以及
δ
≈
g
f
⋅
f
′
g
′
⋅
τ
′
′
−
τ
τ
′
′
−
τ
′
⋅
det
(
w
′
′
∣
w
′
∣
W
′
)
det
(
w
′
′
∣
w
∣
W
′
)
δ
′
\delta\approx \frac{g}{f}\cdot\frac{f'}{g'}\cdot\frac{\tau''-\tau}{\tau''-\tau'}\cdot\frac{\det (\bold{w}''|\bold{w}'|\bold{W}')}{\det (\bold{w}''|\bold{w}|\bold{W}')}\delta'
δ≈fg⋅g′f′⋅τ′′−τ′τ′′−τ⋅det(w′′∣w∣W′)det(w′′∣w′∣W′)δ′
我们一旦简单地近似这些比值
g
f
,
f
′
g
′
,
g
′
′
f
′
′
\frac{g}{f}, \frac{f'}{g'}, \frac{g''}{f''}
fg,g′f′,f′′g′′
为1,我们就基于
δ
′
\delta'
δ′得到
δ
,
δ
′
′
\delta,\delta''
δ,δ′′的简单估计。合理性在于当角度很小时,椭圆弧和弦几乎重合,所以面积差很小。另外,由于
g
f
>
1
,
f
′
g
′
<
1
,
g
′
′
f
′
′
>
1
,
\frac{g}{f}>1, \frac{f'}{g'}<1, \frac{g''}{f''}>1,
fg>1,g′f′<1,f′′g′′>1,
这些误差互相之间也会抵消一部分。后面我们可以看到,Gauss发明了一种聪明的方法来处理当误差变得太大时的情况去调整结果。
此外,Gauss 1809中也提出了一个稍微更好的近似方法。
记
p
,
p
′
,
p
′
′
p,p',p''
p,p′,p′′ 为
τ
,
τ
′
,
τ
′
′
\tau,\tau',\tau''
τ,τ′,τ′′ 三个时刻Ceres的位置。我们从一个直觉的近似开始,即三角形
p
p
′
p
′
′
pp'p''
pp′p′′ 的面积约为椭圆和弦
p
p
′
′
pp''
pp′′ 之间面积的 3/4,表示为
f
+
f
′
+
f
′
′
≈
3
4
(
f
′
−
g
′
)
or
f
′
−
g
′
≈
4
3
(
f
+
f
′
+
f
′
′
)
f+f'+f'' \approx\frac{3}{4}(f'-g') \text{ or } f'-g'\approx \frac{4}{3}(f+f'+f'')
f+f′+f′′≈43(f′−g′) or f′−g′≈34(f+f′+f′′)
给右式两边除以
−
f
′
-f'
−f′ 并加1后得到
g
′
f
′
≈
1
−
4
3
f
+
f
′
+
f
′
′
f
′
\frac{g'}{f'}\approx 1-\frac{4}{3}\frac{f+f'+f''}{f'}
f′g′≈1−34f′f+f′+f′′
右边的比值我们还是不知道,不过下一部分我们可以用已知的量和
r
′
r'
r′ 来估计(注意
r
′
r'
r′ 可以根据
δ
′
\delta'
δ′ 计算得到,且这里我们假设是已知的)。然后求倒数我们就得到上面需要估计的比值。
不管上面如何近似,我们已经获得 δ , δ ′ ′ \delta,\delta'' δ,δ′′ 的估计。乘以对应的向量 w , w ′ ′ \bold{w},\bold{w}'' w,w′′ 后,我们就得到 Ceres 的地心位置向量 δ w = ( ξ , η , ζ ) T \delta\bold{w}=(\xi,\eta,\zeta)^T δw=(ξ,η,ζ)T, δ ′ ′ w ′ ′ = ( ξ ′ ′ , η ′ ′ , ζ ′ ′ ) T \delta''\bold{w}''=(\xi'',\eta'',\zeta'')^T δ′′w′′=(ξ′′,η′′,ζ′′)T。然后我们加上地球的日心坐标 ( X , Y , Z ) , ( X ′ ′ , Y ′ ′ , Z ′ ′ ) (X,Y,Z), (X'',Y'',Z'') (X,Y,Z),(X′′,Y′′,Z′′) 就得到 Ceres 的日心坐标向量 r , r ′ ′ \bold{r},\bold{r}'' r,r′′。