ねこふくろうのメモ置き場𝕏
記事一覧に戻る

天体観測で地球の公転軌道の離心率をもとめる 2/2~フィッティング編~

公開日:2026/08/28

導入

前の記事では天体観測で公転の移動角度を計測する方法を考えました。
今回はその観測結果から離心率を求めるためのフィッティング手法を解説します。
具体的には「移動経路の面積の公式の導出」「観測データの作成」「フィッティング」「フィッシャー情報行列による誤差推定」の4つのパートに分けて説明する。

今回使用したコードはgithubで公開しています。

本文

移動経路の面積を求める公式の導出

まず前の記事で導出を後回しにした、近日点からの移動角度θ\thetaに対する公転軌道の移動経路の面積S(θ,α)S(\theta, \alpha)を導出する。
ここでα\alphaは短径と長径の比を表す。
fig-1: 公転軌道の位置関係と角度
fig-1: 公転軌道の位置関係と角度
導出方針は「Oを中心とするAからBまでの弧の面積（Sarc(θ,α)S_{arc}(\theta, \alpha)）」から「OPBの3点を結ぶ三角形の面積（Stri(θ,α)S_{tri}(\theta, \alpha)）」を引いて求める。
そのために上の図から簡単に分かる式を記載しておく。
公転軌道の楕円は式(1)で表せる。
x2α2+y2=1(1)\frac{x^2}{\alpha^2} + y^2 = 1 \tag {1}
また、PBを結ぶ直線は式(2)で表せる。
y=xtan⁡(θ−π2)−1−α2(2)y=x\tan(\theta- \frac{\pi}{2})-\sqrt{1-\alpha^2} \tag {2}
∠AOB\angle AOBを求めるためにBの座標を求める必要があるため、次の項で求める。

Bの座標を求める

式(1)に式(2)を代入すると
x2α2+(xtan⁡(θ−π2)−1−α2)2=1(3)\frac{x^2}{\alpha^2} + \left(x\tan\left(\theta-\frac{\pi}{2}\right) - \sqrt{1-\alpha^2}\right)^2 = 1 \tag{3}
となり、式(3)をxで解くと式(4)を得る。
x=α2(1−α2tan⁡(θ−π2)±1+tan⁡2(θ−π2))α2tan⁡2(θ−π2)+1(4)x= \frac{\alpha^{2} \left(\sqrt{1 - \alpha^{2}} \tan{\left(\theta-\frac{\pi}{2}\right)} \pm \sqrt{1+\tan^2\left(\theta-\frac{\pi}{2}\right)}\right)}{\alpha^{2} \tan^{2}{\left(\theta-\frac{\pi}{2}\right)} + 1} \tag{4}
となる。
0≤θ≤π0\leq \theta \leq \piの時0≤x0\leq xであり、π≤θ≤2π\pi \leq \theta \leq 2\piの時x≤0x \leq 0となるためθ\thetaの値によって式(4)のプラスを採用するかマイナスを採用するかが変わるため、どちらを採用するべきかを考える必要がある。
まず、当然分母は正であるため分子のみに着目すればよいことが分かる。
また、tan⁡(θ−π2)≤0\tan\left(\theta-\frac{\pi}{2}\right) \leq 0の時は1−α2tan⁡(θ−π2)≤0<1+tan⁡2(θ−π2)\sqrt{1 - \alpha^{2}} \tan{\left(\theta-\frac{\pi}{2}\right)} \leq 0 < \sqrt{1+\tan^2\left(\theta-\frac{\pi}{2}\right)}である。
なおかつ、tan⁡(θ−π2)>0\tan\left(\theta-\frac{\pi}{2}\right) > 0の時は1−α2tan⁡(θ−π2)≤tan⁡(θ−π2)<1+tan⁡2(θ−π2)\sqrt{1 - \alpha^{2}} \tan{\left(\theta-\frac{\pi}{2}\right)} \leq \tan\left(\theta-\frac{\pi}{2}\right) < \sqrt{1+\tan^2\left(\theta-\frac{\pi}{2}\right)}である。
以上より常に1−α2tan⁡(θ−π2)≤1+tan⁡2(θ−π2)\sqrt{1 - \alpha^{2}} \tan{\left(\theta-\frac{\pi}{2}\right)} \leq \sqrt{1+\tan^2\left(\theta-\frac{\pi}{2}\right)}である。よってxxは
x={α2(1−α2tan⁡(θ−π2)+1+tan⁡2(θ−π2))α2tan⁡2(θ−π2)+1(0≤θ≤π)α2(1−α2tan⁡(θ−π2)−1+tan⁡2(θ−π2))α2tan⁡2(θ−π2)+1(π<θ≤2π)(5) x= \begin{cases} \frac{\alpha^{2} \left(\sqrt{1 - \alpha^{2}} \tan{\left(\theta-\frac{\pi}{2}\right)} + \sqrt{1+\tan^2\left(\theta-\frac{\pi}{2}\right)}\right)}{\alpha^{2} \tan^{2}{\left(\theta-\frac{\pi}{2}\right)} + 1} & (0 \le \theta \le \pi) \\ \frac{\alpha^{2} \left(\sqrt{1 - \alpha^{2}} \tan{\left(\theta-\frac{\pi}{2}\right)} - \sqrt{1+\tan^2\left(\theta-\frac{\pi}{2}\right)}\right)}{\alpha^{2} \tan^{2}{\left(\theta-\frac{\pi}{2}\right)} + 1} & (\pi < \theta \le 2\pi) \\ \end{cases} \tag{5}
で与えられる。yyは上で求めたxxを代入すればよい。

Oを中心とするAからBまでの弧の面積（Sarc(θ,α)S_{arc}(\theta, \alpha)）を求める

∠AOB\angle AOBはarctan⁡(y/x)\arctan(y/x)で表せる。これより楕円の弧の面積は式(6)で表せる。
Sarc(θ,α)=α2arctan⁡(yx)(6)S_{arc}(\theta, \alpha) = \frac{\alpha}{2} \arctan{\left( \frac{y}{x} \right)} \tag{6}
xxとyyは前項で求めたとおりθ\thetaとα\alphaにのみ依存するため注意すること。
また、tan⁡(θ)=tan⁡(θ+π)\tan(\theta) = \tan(\theta+\pi) であるため実装の際は工夫する必要がある。
方法は複数あり、本筋ではないため説明は省く。

恒星と中点と惑星の３点を結ぶ３角形の面積を求める

三角形の面積も図から分かるとおり
Stri(θ,α)=x1−α22(7)S_{tri}(\theta, \alpha) = \frac{x\sqrt{1-\alpha^2}}{2} \tag{7}
と書ける。
このとき、(π<θ≤2π)(\pi < \theta \le 2\pi)ではStri(θ,α)S_{tri}(\theta, \alpha)は負の値になる。
ただし、fig-2で分かる通り(π<θ≤2π)(\pi < \theta \le 2\pi)では弧の面積と三角形の面積を足す必要がある。
そのため面積を加算するか減算するかはStri(θ,α)S_{tri}(\theta, \alpha)の正負で吸収されるため場合分けは不要である。
fig-2: 公転軌道の位置関係と角度(π<θ≦2π)
fig-2: 公転軌道の位置関係と角度(π<θ≦2π)

求めた２つの面積から移動経路の面積を求める

最初に述べた通り移動経路の面積は楕円の弧の面積から三角形の面積を引けば良いため式(8)となる。
S(θ,α)=Sarc(θ,α)−Stri(θ,α)(8)S(\theta, \alpha) = S_{arc}(\theta, \alpha) - S_{tri}(\theta, \alpha) \tag{8}
以上により式(8)で公転軌道の移動面積を求めることができるようになった。
これを利用し数日後の地球の位置などを推定することができる。
以下ではそれを利用し観測データのフィッティングを行う。

観測データを作成する

Skyfieldライブラリを使用して2025年の地球の公転データから観測角度を計算しそれにノイズを加えた。
観測期間は以下の2つの期間に分けて観測した
・2025/01/15 〜 2025/06/16（近日点付近から半年分）
・2025/07/10 〜 2025/12/27（遠日点付近から半年分）
分けた理由はy/xが無限に発散することを避けるためである。

ノイズパラメータ

各観測日の地球から太陽方向ベクトルを計算し、0日目からの角度変化を公転移動角度YtY_tとして取得する。
実際の観測では計測誤差があるため、真値に正規分布ノイズを加えた：
Yt=Yttrue+w,w∼N(0,σ2)(9)Y_t = Y_t^{\text{true}} + w, \quad w \sim N(0, \sigma^2) \tag{9}
今回はσ2=0.01\sigma^2 = 0.01を使用した。
本来であれば、前の記事で述べたように観測地点を変える必要があるためそのたびにノイズが重なっていく。
今回は問題の簡単化のために一定にした。

図示して見る

ノイズを加えた観測データと真値の差をプロットするとfig-3とfig-4のようになる。
fig-3: 誤差のプロット
fig-3: 誤差のプロット
fig-4: 観測データと真値
fig-4: 観測データと真値
fig-4では誤差が小さくて真値のラインしか見えないので見る必要はない。
ただ、ほぼ直線に見えることだけ把握しておいてほしい。
これは地球の公転軌道がほぼ真円に近いことを表す。

フィッティングする

最小二乗法でパラメータα^,θ^0,x^P\hat{\alpha}, \hat{\theta}_0, \hat{x}_Pを推定する。
xPx_{P}は遠日点を通過する前に観測した最後の角度の真値である。
xPx_{P}も予測する値として導入する理由は正確な面積速度を予測するためxPx_{P}もノイズを除去した値を推定する値に入れて目的関数に組み込む必要がある。
目的関数は式(10)とする。
q(α,θ0,xP)=∑t=1T(Yt−xt(α,θ0,xP))2(10)q(\alpha, \theta_0, x_P) = \sum_{t=1}^{T} \bigl(Y_t - x_t(\alpha, \theta_0, x_P)\bigr)^2 \tag {10}
ここでxtx_tはt毎に式(11)を数値解析的に解いた値である。
S(θ0+xP)−S(θ0)TP−S(θ0+xt)−S(θ0)Tt=0(11)\frac{S(\theta_0 + x_P) - S(\theta_0)} {T_P} - \frac{S(\theta_0 + x_t) - S(\theta_0)} {T_t} = 0 \tag {11}
式(11)は前の記事で紹介したケプラーの第二法則を利用したものである
TTは観測開始時刻からの経過時間を表す。

予測結果

最適化の結果は以下の通りであった：
α^≈0.99984896,θ^0≈0.25076657,x^N≈4.18790703\hat{\alpha} \approx 0.99984896, \quad \hat{\theta}_0 \approx 0.25076657, \quad \hat{x}_N \approx 4.18790703

fig-5は観測値と予測値から真値を引いたものである。
fig-5から予測値のほうが良い結果を得られていることが分かる。
これによりフィッティングが正常に機能していることが確認できた。
fig-5: 観測データと真値
fig-5: 観測データと真値
次はフィッシャー情報行列を計算することで数値的に予測値の信頼度を測ってみる。

誤差推定【フィッシャー情報行列】

フィッシャー情報行列とは

フィッシャー情報行列は簡単に説明すると「データがパラメータについてどれだけ正確な情報を持っているかを表すもの」である。
詳しい説明は別の資料を参考にしてもらいたい。
今回重要なのは次の項の式(12)と式(13)のように推定したパラメータの分散を推定することができることである。

QQのフィッシャー情報行列を求める

推定量の分散の下限はクラメール・ラオの不等式により
Var(Q^)≥[I(Q)]−1(12)\mathrm{Var}(\hat{Q}) \geq [I(Q)]^{-1} \tag {12}
で与えられる。ここでI(θ)I(\theta)はフィッシャー情報行列である。
フィッシャー情報行列は以下の式で表せる。
I(Q)=E ⁣x∣Q[∂2log⁡p(x∣Q)QQ⊤](13)I(Q)=\mathbb{E}\!_{x\mid Q}\left[\frac{\partial^2 \log p(x\mid Q)} {QQ^\top} \right] \tag {13}
p(x∣Q)=∏n=1N12πσ2exp⁡(−(xt−yt)22σ2)(14)p(x\mid Q) = \prod_{n=1}^N \frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(x_t-y_t)^2}{2\sigma^2}\right) \tag {14}
ただしQ=(α,θ0,xP)⊤Q=(\alpha,\theta_0, x_P)^\topである。
p(x∣Q)p(x\mid Q)を整理すると以下の式に変形できる。
p(x∣Q)=(12πσ2)Nexp⁡(−∑t=1N(xt−yt)22σ2)(14)p(x\mid Q) = \left(\frac{1}{\sqrt{2\pi\sigma^2}}\right)^N \exp\left(-\sum_{t=1}^N\frac{(x_t-y_t)^2}{2\sigma^2}\right) \tag {14}
よってlog⁡p(x∣Q)\log p(x\mid Q)は以下のようになる
log⁡p(x∣Q)=Nlog⁡(12πσ2)−12σ2∑t=1N(xt−yt)2(15)\log p(x\mid Q) = N \log \left(\frac{1}{\sqrt{2\pi\sigma^2}}\right) - \frac{1}{{2\sigma^2}} \sum_{t=1}^N{(x_t-y_t)^2} \tag {15}
よってフィッシャー情報行列を求めるためにr,s∈{θ0,α,xP}r, s \in \{ \theta_0, \alpha, x_P \}として
∂2log⁡p(x∣Q)∂r∂s=12σ2∑t=1N∂2(xt−yt)2∂r∂s(16 )\frac{\partial^2 \log p(x\mid Q)} {\partial r \partial s} = \frac{1}{{2\sigma^2}} \sum_{t=1}^N{\frac {\partial^2 (x_t-y_t)^2}{\partial r \partial s}} \tag {16 }
を計算する必要がある。
式を見れば分かるとおり、xtx_tさえrrで微分できれば連鎖律により解析的に計算できる。
xtx_tは式(11)を満たす変数xtx_tであり、陰関数定理を適用できるため∂xt/∂r{\partial x_t}/{\partial r}を解析的に求めることができる。以下のような式で計算できる。
∂xt∂r=−∂(S(θ0+xP)−S(θ0)TP−S(θ0+xt)−S(θ0)Tt)∂r(∂(S(θ0+xP)−S(θ0)TP−S(θ0+xt)−S(θ0)Tt)∂xt)−1(17)\frac{\partial x_t}{\partial r} = - \frac {\partial \left(\frac{S(\theta_0 + x_P) - S(\theta_0)} {T_P} - \frac{S(\theta_0 + x_t) - S(\theta_0)} {T_t}\right)} {\partial r} \left( \frac {\partial \left(\frac{S(\theta_0 + x_P) - S(\theta_0)} {T_P} - \frac{S(\theta_0 + x_t) - S(\theta_0)} {T_t}\right)} {\partial x_t} \right)^{-1} \tag {17}
ただし、本来は式(11)のヤコビ行列の行列式が0でないことを確かめなければならないがフィッティングの結果が安定しているためヤコビ行列の行列式は0でないと推測できるため計算は省く。

※1　実際に導出しようとすると長くなるので計算は省く。
※2　陰関数定理を説明しようとすると長くなるため詳細が分からない方は別の資料を参照してもらいたい。簡単に説明するとg(x, y)=0という関数に対してy=f(x)という式が存在するかどうかを判定する定理である。

これにより式(13)は解析的に解くことができ、その逆行列I−1I^{-1}の対角成分が各パラメータの推定誤差の下限を与える：
I−1≈(9.25×10−7−7.77×10−4−1.29×10−5−7.77×10−44.92−0.176−1.29×10−5−0.1761.02×10−2)I^{-1} \approx \begin{pmatrix} 9.25 \times 10^{-7} & -7.77 \times 10^{-4} & -1.29 \times 10^{-5} \\ -7.77 \times 10^{-4} & 4.92 & -0.176 \\ -1.29 \times 10^{-5} & -0.176 & 1.02 \times 10^{-2} \end{pmatrix}

この結果からα\alphaは十分に良い結果を得たと言える一方でθ0\theta_0の分散はだいぶ大きいことが分かる。
そのためほとんど信頼できない結果であると言ってよい。
xPx_Pの分散は観測データの作成の際に加えたノイズの分散とほとんど変わらない（むしろ悪くなっている）。
そのため、推定値と観測値でほとんど信頼性が変わらないことが分かる。

これによりα\alphaを求めるために作成した関数達はα\alphaを求める目的では効果的であることが分かる。
一方で地球の現在地（これも気になるデータではある）を推定するには不十分であることも分かる。

今回のゴール設定がα\alphaから離心率をもとめることであったため、目的を達成するには十分である。
計算すると離心率は0.0173となる。一般的に地球の公転軌道の離心率は0.0167と紹介されるため悪くない値である。
ただ、観測誤差をだいぶ小さくとっているため実際に行うともっと悪い結果がでるはずであるため、観測データを増やすことで最頻値を求めるなどの工夫が必要である。
今回の方法の場合、観測を開始する時刻は何時でも良いため1日に1回の観測ではなく何度も観測し、それでフィッティングすれば観測精度を上げることができる。

最後に

以下は感想なので読み飛ばしてもらって結構です。
今回は公転軌道の面積公式の導出からフィッティング・誤差推定までの一連の流れをまとめました。
フィッティングの方法を求めるまでに紆余曲折がありました。
最初はパラメータのフィッティングをカルマンフィルタで実施しようと考えていました。
理由は観測結果が時系列であるため、カルマンフィルタにより徐々にパラメータを真値に近づけていこうと考えていたためです。
ただ、いろいろ難しく考えた結果、単純な方法が最も良いというありがちな結論になったのは若干モヤモヤが残っています。
せっかくカルマンフィルタについて学習したのでどこかでそれを利用した記事を書きたいと考えています（ネタはありませんが...）。
他にもxtx_tを計算するためにどうするかという方法もありました。
最初の最初はxtx_tを解析的に求め、どれくらいの観測精度が必要かを直接計算しようと考えていましたが、解析的にはどうにもなりそうにないことが分かり諦めました。
今回は個人で観測を行うには地球一周旅行を行う必要がでましたが、今後同じ地点から地球の離心率を求める方法が思いつけばそれを紹介したいとおもいます（地軸の傾きを求めれば可能かも？）。
ここまで読んでいただきありがとうございました。
ぜひ次の記事も読んでもらえますと幸いです。

参考文献

片山徹『非線形カルマンフィルタ』朝倉書店、2011年、ISBN 978-4-254-20148-2。


記事一覧に戻る