太陽同期軌道の理論と設計 — J2摂動を利用した軌道設計をわかりやすく解説

Landsat衛星が撮影した地球の画像を見たことがあるでしょうか。何十年にもわたって同じ地域を同じ照明条件で撮影し続け、森林破壊の進行や都市の拡大、氷河の後退を記録しています。なぜこれほど一貫した画像を撮り続けられるのでしょうか? 答えは太陽同期軌道(Sun-Synchronous Orbit, SSO)にあります。

地球観測衛星にとって、太陽の当たる角度が撮影のたびに変わってしまうのは致命的です。影の向きが変わると、植生の反射スペクトルの比較ができなくなり、季節変化の長期モニタリングが不可能になります。太陽同期軌道は、地球の扁平さ(J2摂動)を巧みに利用して、衛星が常に同じ太陽角度のもとで地表を観測できるようにする軌道です。

太陽同期軌道を理解すると、以下のような応用が見えてきます。

  • 地球観測ミッション設計: Landsat、Sentinel-2、だいち(ALOS)のような衛星の軌道パラメータを設計できる
  • リモートセンシングの画像品質向上: 一定の照明条件を保つことで、マルチテンポラル解析(時系列比較)の精度を高められる
  • 回帰軌道の設計: 定期的に同じ地上軌跡を通る軌道を設計し、定点観測を実現できる

本記事の内容

  • 太陽同期軌道の直感的な理解と物理的背景
  • J2摂動による昇交点赤経の歳差の数学的導出
  • 太陽同期条件の導出と高度・傾斜角の関係
  • 楕円軌道への拡張
  • 回帰軌道(repeat ground track)の設計原理
  • 地方太陽時(LST)の意味と選択基準
  • Pythonによる太陽同期軌道の設計と可視化

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

太陽同期軌道とは — 概念の直感的理解

地球観測の基本的な悩み

地球を周回する衛星を想像してください。衛星は地球のまわりを約90分で1周しますが、その間に地球は太陽のまわりを少しだけ公転しています。1日あたり約 $1°$ です。もし衛星の軌道面が宇宙空間に対して固定されていたら(慣性空間で固定)、地球が太陽のまわりを回るにつれて、衛星が地表を通過する時刻が毎日少しずつずれていきます。

たとえば、最初は午前10時に赤道を通過していた衛星が、3か月後には午後1時に通過するようになります。太陽の高度角が変わるため、同じ場所でも影の長さや方向が大きく異なり、画像の比較が困難になります。

軌道面を回す発想

この問題を解決するには、衛星の軌道面を毎日少しずつ回転させて、地球の公転に追従させればよいのです。地球は太陽のまわりを365.25日で $360°$ 回るので、1日あたり

$$ \dot{\Omega}_{\text{sun}} = \frac{360°}{365.2422 \, \text{日}} \approx 0.9856°/\text{日} $$

だけ、軌道面を東向きに(赤経が増加する方向に)回転させれば、衛星と太陽の幾何学的関係が一年を通じて保たれます。

ここで $\Omega$ は昇交点赤経(Right Ascension of the Ascending Node, RAAN)と呼ばれる軌道要素で、軌道面が赤道面と交わる点(昇交点)の方向を春分点から測った角度です。

自然の力を利用する

「軌道面を毎日0.9856° 回転させる」と聞くと、スラスターで常時噴射が必要に思えるかもしれません。しかし、実際にはロケットエンジンを使う必要はありません。地球が完全な球ではなく、赤道方向にわずかに膨らんだ回転楕円体であることが、まさにこの回転を自然にもたらしてくれるのです。

地球の赤道膨らみによる重力の非対称性は、J2摂動と呼ばれます。この摂動は衛星の軌道面に対してトルクを与え、昇交点赤経 $\Omega$ を時間とともに変化させます。適切な軌道傾斜角と高度の組み合わせを選べば、このJ2摂動による $\Omega$ の変化率をちょうど $+0.9856°/\text{日}$ に設定することができます。これが太陽同期軌道の基本原理です。

太陽同期軌道の原理が直感的にわかったところで、次にこの原理を数学的に定式化するために、地球の公転と昇交点赤経の歳差について詳しく見ていきましょう。

地球の公転と昇交点赤経の歳差

慣性系から見た地球の公転

太陽同期条件を正確に理解するために、慣性座標系(春分点と恒星を基準とする座標系)で考えます。地球は太陽のまわりを1恒星年($T_{\text{sid}} = 365.2564$ 日)で公転していますが、太陽同期軌道の設計では1太陽年(回帰年: $T_{\text{trop}} = 365.2422$ 日)を使います。これは、私たちが気にするのは「太陽に対する衛星の向き」であり、春分点を基準とした太陽の見かけの動きが回帰年で一周するためです。

地球から見た太陽の黄経(ecliptic longitude)の変化率は

$$ \dot{\lambda}_{\text{sun}} = \frac{360°}{365.2422 \, \text{日}} \approx 0.9856°/\text{日} $$

です。太陽同期軌道では、衛星の昇交点赤経の歳差率がこの値に一致することを要求します。

昇交点赤経の意味の復習

昇交点赤経 $\Omega$ は、衛星が南半球側から北半球側へ赤道を横切る点(昇交点)の方向を、春分点 $\gamma$ から東向きに測った角度です。$\Omega$ が時間とともに変化するということは、軌道面が赤道面に対する傾斜角を保ったまま、地球の自転軸のまわりをコマの歳差のように回転することを意味します。

J2摂動がない場合(地球が完全な球の場合)、$\Omega$ は一定です。しかし実際の地球は赤道方向に膨らんでいるため、この膨らみが衛星にかける重力トルクによって $\Omega$ が時間とともに変化します。

ここで疑問が生まれます。J2摂動による $\Omega$ の変化率は具体的にどのような式で表されるのでしょうか? 次のセクションでこの重要な式を導出します。

J2摂動による昇交点赤経の変化率の導出

地球の重力ポテンシャル

地球が完全な球であれば、重力ポテンシャルは質点の場合と同じ $U = -\mu/r$ です。しかし実際の地球は赤道方向にわずかに膨らんでおり、この非球形性を帯調和関数(zonal harmonics)で表現します。主要な非球形項はJ2項で、地球の重力ポテンシャルは

$$ U = -\frac{\mu}{r}\left[1 – J_2 \left(\frac{R_E}{r}\right)^2 P_2(\sin\phi)\right] $$

と書けます。ここで $\mu = GM_{\oplus} = 3.986 \times 10^{14} \, \text{m}^3/\text{s}^2$ は地球の重力パラメータ、$R_E = 6{,}378.137 \, \text{km}$ は地球の赤道半径、$J_2 = 1.08263 \times 10^{-3}$ は地球のJ2係数、$P_2$ は2次のルジャンドル多項式、$\phi$ は地心緯度です。

2次のルジャンドル多項式は

$$ P_2(\sin\phi) = \frac{3\sin^2\phi – 1}{2} $$

です。$\sin\phi$ は衛星の軌道面と赤道面の幾何学的関係から、軌道要素を使って表現できます。

摂動方程式からの導出

J2摂動がRAANに与える影響を求めるには、ラグランジュの惑星方程式(Lagrange’s planetary equations)を用います。この方程式は、摂動ポテンシャル $R$ に対する軌道要素の時間変化率を与えます。昇交点赤経 $\Omega$ に関するラグランジュの惑星方程式は

$$ \frac{d\Omega}{dt} = \frac{1}{na^2\sqrt{1-e^2}\sin i} \frac{\partial R}{\partial i} $$

です。ここで $n$ は平均運動(mean motion)、$a$ は軌道長半径、$e$ は離心率、$i$ は軌道傾斜角です。

J2摂動ポテンシャルを軌道要素で表現し、1軌道周期にわたって平均化する操作(secular perturbation, 永年摂動)を行うと、短周期変動が消えて長期的なトレンドだけが残ります。

平均化された変化率

詳細な計算(摂動関数の軌道要素への展開と1周期平均)を経ると、J2摂動による昇交点赤経の永年変化率は次のように得られます。

$$ \dot{\Omega}_{J2} = -\frac{3}{2} n J_2 \left(\frac{R_E}{a}\right)^2 \frac{\cos i}{(1-e^2)^2} $$

この式の物理的意味を1つずつ確認しましょう。

  • $\cos i$ の項: 軌道傾斜角 $i$ が $90°$ より小さい(順行軌道)なら $\cos i > 0$ で $\dot{\Omega}_{J2} < 0$(西向きの歳差)。$i > 90°$(逆行軌道)なら $\cos i < 0$ で $\dot{\Omega}_{J2} > 0$(東向きの歳差)。太陽同期軌道は東向きの歳差が必要なので、$i > 90°$ の逆行軌道でなければなりません。
  • $(R_E/a)^2$ の項: 軌道高度が低いほど($a$ が小さいほど)歳差率の絶対値が大きくなります。地球の膨らみに近い方が、その影響を強く受けるためです。
  • $n$ の項: 平均運動 $n = \sqrt{\mu/a^3}$ なので、高度が低いほど周回が速く、J2摂動の影響を受ける頻度も高くなります。
  • $(1-e^2)^{-2}$ の項: 離心率が大きい(楕円的な)軌道ほど、近地点付近で地球に接近する時間が長くなるため、歳差率が大きくなります。

平均運動 $n$ を $\mu$ と $a$ で明示的に書くと

$$ \dot{\Omega}_{J2} = -\frac{3}{2} \sqrt{\frac{\mu}{a^3}} \, J_2 \left(\frac{R_E}{a}\right)^2 \frac{\cos i}{(1-e^2)^2} $$

整理すると

$$ \dot{\Omega}_{J2} = -\frac{3}{2} \frac{J_2 R_E^2 \sqrt{\mu}}{a^{7/2}} \frac{\cos i}{(1-e^2)^2} $$

となります。この式が太陽同期軌道の設計における最も重要な方程式です。

J2摂動によるRAANの変化率がわかったところで、いよいよこの式を使って太陽同期条件を具体的に導出しましょう。

太陽同期条件の導出

条件式

太陽同期軌道の定義は、「衛星の昇交点赤経の歳差率が、地球の太陽まわりの公転による太陽方向の変化率に等しい」ことです。数式で書くと

$$ \dot{\Omega}_{J2} = \dot{\Omega}_{\text{sun}} $$

です。右辺は地球の公転角速度で

$$ \dot{\Omega}_{\text{sun}} = \frac{2\pi}{T_{\text{trop}}} = \frac{2\pi}{365.2422 \times 86400} \approx 1.991 \times 10^{-7} \, \text{rad/s} $$

これを度に変換すると $\dot{\Omega}_{\text{sun}} \approx 0.9856°/\text{日}$ です。

J2摂動の式を代入すると、太陽同期条件は

$$ -\frac{3}{2} \frac{J_2 R_E^2 \sqrt{\mu}}{a^{7/2}} \frac{\cos i}{(1-e^2)^2} = \dot{\Omega}_{\text{sun}} $$

となります。ここで $\dot{\Omega}_{\text{sun}} > 0$(東向き)であり、左辺にマイナス符号があるため、$\cos i < 0$ でなければなりません。つまり $i > 90°$ です。これは、太陽同期軌道が必ず逆行軌道(retrograde orbit)であることを数学的に示しています。

$\cos i$ について解く

上の条件式を $\cos i$ について整理します。

$$ \cos i = -\frac{2 \dot{\Omega}_{\text{sun}} a^{7/2}(1-e^2)^2}{3 J_2 R_E^2 \sqrt{\mu}} $$

この式が太陽同期軌道の設計方程式です。軌道長半径 $a$ と離心率 $e$ を指定すれば、太陽同期条件を満たす軌道傾斜角 $i$ が一意に決まります。

物理的な制約

$\cos i$ は $-1$ 以上 $1$ 以下でなければならないため、上の式の右辺の絶対値が $1$ を超えてはいけません。これは $a$(つまり高度)に上限があることを意味します。高度が高すぎると、J2摂動の影響が弱くなりすぎて、必要な歳差率を得ることができなくなるのです。

具体的には、円軌道($e = 0$)で $i = 180°$($\cos i = -1$)のとき

$$ a_{\max} = \left(\frac{3 J_2 R_E^2 \sqrt{\mu}}{2 \dot{\Omega}_{\text{sun}}}\right)^{2/7} $$

を計算すると、$a_{\max} \approx 12{,}352 \, \text{km}$、つまり高度約 $5{,}974 \, \text{km}$ が理論上の上限となります。実用的な太陽同期軌道は高度600〜900 kmに集中しており、傾斜角は $96°$ 〜 $99°$ 程度です。

太陽同期条件の一般式が得られたので、次に最も基本的なケースとして円軌道の場合を詳しく見ていきましょう。

円軌道における高度と軌道傾斜角の関係

円軌道での簡略化

多くの太陽同期軌道衛星は離心率がほぼゼロの円軌道(near-circular orbit)を使用します。$e = 0$ とすると、設計方程式は大幅に簡略化されます。

$(1-e^2)^2 = 1$ なので

$$ \cos i = -\frac{2 \dot{\Omega}_{\text{sun}} a^{7/2}}{3 J_2 R_E^2 \sqrt{\mu}} $$

また、円軌道では $a = R_E + h$($h$ は高度)なので、高度 $h$ を指定すれば太陽同期傾斜角 $i$ が直ちに計算できます。

具体的な数値例

代表的な高度での太陽同期傾斜角を計算してみましょう。

高度 $h = 600 \, \text{km}$ の場合、$a = R_E + h = 6{,}978.137 \, \text{km} = 6.978 \times 10^6 \, \text{m}$ です。

まず右辺の各量を代入します。

$$ \frac{2 \dot{\Omega}_{\text{sun}}}{3 J_2 R_E^2 \sqrt{\mu}} = \frac{2 \times 1.991 \times 10^{-7}}{3 \times 1.08263 \times 10^{-3} \times (6.378 \times 10^6)^2 \times \sqrt{3.986 \times 10^{14}}} $$

分母を計算すると

$$ 3 \times 1.08263 \times 10^{-3} \times 4.068 \times 10^{13} \times 1.996 \times 10^7 \approx 2.638 \times 10^{18} $$

したがって

$$ \frac{2 \dot{\Omega}_{\text{sun}}}{3 J_2 R_E^2 \sqrt{\mu}} \approx \frac{3.982 \times 10^{-7}}{2.638 \times 10^{18}} \approx 1.510 \times 10^{-25} \, \text{m}^{-7/2} $$

$a^{7/2} = (6.978 \times 10^6)^{7/2}$ を計算して掛けると

$$ \cos i \approx -0.1047 $$

よって

$$ i = \arccos(-0.1047) \approx 96.01° $$

同様に、高度 $h = 800 \, \text{km}$ では $i \approx 98.6°$、$h = 500 \, \text{km}$ では $i \approx 97.4°$ となります。

高度と傾斜角の関係の特徴

計算結果から、以下の傾向が読み取れます。

  1. 高度が上がるほど傾斜角は大きくなる — 高度が上がると $a^{7/2}$ が大きくなるため、$|\cos i|$ が増加し、$i$ は $90°$ からさらに離れます。これは、高い軌道ではJ2摂動が弱いため、より大きな $|\cos i|$ で「テコの腕」を増やして歳差率を確保する必要があるためです。
  2. 傾斜角は常に $90°$ よりやや大きい — 実用的なSSO高度(400〜1500 km)では $96°$ 〜 $104°$ 程度です。
  3. $a^{7/2}$ の急峻な依存性 — $7/2 = 3.5$ 乗という強い依存性のため、高度がわずかに変わるだけで必要な傾斜角が大きく変化します。

このように、円軌道の場合は高度と傾斜角が1対1に対応します。しかし、実際のミッションでは円軌道とは限りません。次に、離心率がゼロでない楕円軌道に太陽同期条件を拡張してみましょう。

楕円軌道への拡張

一般的な太陽同期条件

楕円軌道では離心率 $e \neq 0$ となり、設計方程式は

$$ \cos i = -\frac{2 \dot{\Omega}_{\text{sun}} a^{7/2}(1-e^2)^2}{3 J_2 R_E^2 \sqrt{\mu}} $$

のままです。ここで注意すべき点は、$(1-e^2)^2$ の因子です。$e > 0$ のとき $1 – e^2 < 1$ なので $(1-e^2)^2 < 1$ となり、$|\cos i|$ が円軌道の場合より小さくなります。

これは直感的に理解できます。楕円軌道では衛星が近地点で地球に接近するため、J2摂動の効果がより強く蓄積されます。したがって、同じ長半径 $a$ と傾斜角 $i$ の円軌道よりも歳差率が大きくなります。太陽同期歳差率を達成するために必要な $|\cos i|$ は、その分だけ小さくてよい($i$ が $90°$ に近くてよい)のです。

凍結軌道との組み合わせ

楕円軌道の太陽同期衛星では、近地点引数 $\omega$ の変動が問題になります。J2摂動は $\omega$ も変化させ

$$ \dot{\omega}_{J2} = -\frac{3}{2} n J_2 \left(\frac{R_E}{a}\right)^2 \frac{5\cos^2 i – 1}{2(1-e^2)^2} $$

と表されます。$\omega$ が変化すると、近地点と遠地点の位置が軌道面内で回転し、地上軌跡のパターンが崩れます。

この問題を解決するために、凍結軌道(frozen orbit)という概念が使われます。$\dot{\omega}_{J2} = 0$ となる条件は

$$ 5\cos^2 i – 1 = 0 \quad \Rightarrow \quad \cos i = \pm \frac{1}{\sqrt{5}} \quad \Rightarrow \quad i \approx 63.4° \, \text{or} \, 116.6° $$

ですが、より一般的には高次のJ項(J3, J4など)を考慮した凍結軌道条件を使います。実際の地球観測衛星では、太陽同期条件と凍結軌道条件を同時に満たすように $a$, $e$, $i$ を設計します。たとえば、Landsat-7は高度約705 km、離心率 $e \approx 0.001$、傾斜角 $i \approx 98.2°$ の凍結太陽同期軌道を使用しています。

楕円軌道まで考えると設計の自由度が増えますが、地球観測衛星にはもう一つ重要な要件があります。それは、一定の周期で同じ地上軌跡を繰り返す回帰軌道です。次にこの概念を見ていきましょう。

回帰軌道(Repeat Ground Track)の設計

回帰軌道の必要性

地球観測衛星が同じ地点を定期的に観測するためには、衛星の地上軌跡が一定の日数の後に正確に繰り返される必要があります。これを回帰軌道(repeat ground track orbit)と呼びます。

地上軌跡が繰り返されるとはどういうことか、直感的に考えましょう。衛星は1日に約15周(高度600 km程度の場合)地球を回りますが、地球も自転しているため、連続する周回の地上軌跡は西にずれていきます。$D$ 日後に衛星がちょうど $N$ 周回して元の地上軌跡に戻るならば、$D$ 日を回帰日数、$N$ 周を回帰周回数と呼びます。

回帰条件の導出

衛星の軌道周期を $T_{\text{sat}}$ とします。衛星が1周する間に、地球は角度 $\omega_E T_{\text{sat}}$ だけ自転します。ここで $\omega_E = 7.2921 \times 10^{-5} \, \text{rad/s}$ は地球の自転角速度です。

しかし、太陽同期軌道ではRAANが $\dot{\Omega}_{\text{sun}}$ の割合で変化しているため、衛星が赤道を通過する地理経度の変化率は $\omega_E – \dot{\Omega}_{\text{sun}}$ ではなく、正確には衛星のノーダル周期(nodal period)$T_n$ と地球の自転の関係で決まります。

ノーダル周期とは、衛星が昇交点から次の昇交点まで移動する時間で、J2摂動による近地点引数 $\omega$ の変化も考慮した実効的な周期です。

$$ T_n = \frac{2\pi}{n + \dot{\omega}_{J2} + \dot{M}_{J2}} $$

ここで $\dot{M}_{J2}$ はJ2による平均近点角の永年変化率です。

回帰条件は、$D$ 日後に衛星が正確に同じ地上軌跡に戻ることを要求します。これは、$D$ 日間に衛星が描く地上軌跡のパターンが一巡することに相当し

$$ \frac{N}{D} = \frac{T_{\text{day,nodal}}}{T_n} $$

と書けます。ここで $T_{\text{day,nodal}} = 2\pi / (\omega_E – \dot{\Omega}_{\text{sun}})$ はノーダル日(太陽に対する地球の自転周期を、RAANの歳差で補正したもの)です。

$N$ と $D$ は互いに素な正の整数で、$N/D$ は有理数でなければなりません。

回帰軌道の設計手順

実際の設計は以下のように進めます。

  1. まず回帰日数 $D$ と回帰周回数 $N$ を決めます。たとえば $D = 16$ 日、$N = 233$ 周(Landsat)。
  2. $N/D$ の値から必要なノーダル周期 $T_n$ を計算します。
  3. $T_n$ から軌道長半径 $a$ を反復計算で求めます(J2摂動の影響があるため解析解はなく、反復法が必要)。
  4. 求めた $a$ と太陽同期条件から傾斜角 $i$ を決定します。

たとえばLandsatの場合、16日で233周回するので

$$ \frac{N}{D} = \frac{233}{16} = 14 + \frac{9}{16} = 14.5625 \, \text{周/日} $$

これから1日あたり約14.56周回、つまり1周あたり約98.88分(ノーダル周期)が必要です。この値から高度を逆算すると約705 kmが得られます。

代表的な回帰軌道パラメータ

衛星 回帰日数 $D$ 回帰周回数 $N$ 高度 [km] 傾斜角 [°]
Landsat-7/8/9 16 233 705 98.2
Sentinel-2 10 143 786 98.6
だいち2号(ALOS-2) 14 207 628 97.9
SPOT-6/7 26 385 694 98.2

回帰日数 $D$ が小さいほど同じ地点を頻繁に観測できますが、隣接する地上軌跡の間隔が広くなるトレードオフがあります。逆に $D$ が大きいと軌跡の間隔が密になりますが、再訪間隔が長くなります。

回帰軌道の設計は地上軌跡の空間分解能と時間分解能のバランスを決定します。しかし、地球観測にはもう一つ重要なパラメータがあります。それは、衛星が赤道を通過するときの現地時刻 — 地方太陽時です。

地方太陽時(LST)の意味と選び方

地方太陽時とは

地方太陽時(Local Solar Time, LST)は、ある地点において太陽の位置から定義される時刻です。太陽が南中(最も高い位置に来る)するのが正午(12:00 LST)です。

太陽同期軌道では、衛星が昇交点(赤道を南から北へ横切る点)を通過するときの地方太陽時をLTAN(Local Time of Ascending Node)と呼びます。太陽同期軌道の最大の特徴は、このLTANが一年を通じて(ほぼ)一定に保たれることです。

なぜLTANが重要か

LTANの値は、衛星が地表を観測するときの太陽の高度角と方位角を決定します。

  • LTAN = 10:30(午前): 適度な太陽高度で影が残り、地形の凹凸がわかりやすい。Landsat(10:00頃)やSentinel-2(10:30頃)はこのタイプです。光学観測に最適とされます。
  • LTAN = 06:00(夜明け/夕暮れ): 太陽がちょうど地平線上にあり、ターミネーター(昼夜境界線)付近を飛行します。影が非常に長く、微小な地形変化の検出に有利です。
  • LTAN = 12:00(正午): 太陽が最も高く、影が最小になります。海面温度の計測や大気観測に適しています。
  • LTAN = 18:00: 夕暮れの太陽同期軌道で、合成開口レーダー(SAR)衛星によく使われます。SARは自ら電波を照射するため太陽光は不要ですが、18:00軌道は地上の電力需要ピーク時にデータ配信が可能という運用上の利点があります。

LTANとRAANの関係

LTANは昇交点赤経 $\Omega$ と太陽の赤経 $\alpha_{\text{sun}}$ の差で決まります。

$$ \text{LTAN} = 12 \, \text{h} + \frac{\Omega – \alpha_{\text{sun}}}{15°/\text{h}} $$

ここで $15°/\text{h}$ は地球の自転角速度($360°/24 \, \text{h}$)です。太陽同期軌道では $\Omega – \alpha_{\text{sun}}$ が一定に保たれるため、LTANも一定になります。

LTANの選び方

LTANの選択は、ミッションの目的に直結する重要な設計判断です。

用途 推奨LTAN 理由
光学地球観測 10:00〜10:30 適度な太陽高度、影あり
高分解能光学 10:30〜11:00 十分な照度、影で3D情報
海面・大気温度 13:00〜14:00 日中最高温度に近い
SAR 06:00 or 18:00 太陽光不要、運用上の利点
気象観測 09:00, 14:00等 他衛星と補完するタイミング

LTANを決定すると、打ち上げ日時に応じた初期RAAN $\Omega_0$ が決まり、太陽同期軌道の全パラメータが確定します。

ここまでで太陽同期軌道の理論を一通り学びました。次は、これらの理論をPythonコードで実装し、実際に軌道パラメータを計算してみましょう。

Pythonによる太陽同期軌道の設計と可視化

基本定数と太陽同期傾斜角の計算

まず、高度を指定して太陽同期条件を満たす軌道傾斜角を計算する関数を実装します。円軌道と楕円軌道の両方に対応します。

import numpy as np
import matplotlib.pyplot as plt

# --- 基本定数 ---
mu = 3.986004418e14       # 地球の重力パラメータ [m^3/s^2]
R_E = 6378.137e3          # 地球の赤道半径 [m]
J2 = 1.08263e-3           # J2係数
T_trop = 365.2422 * 86400 # 1太陽年(回帰年) [s]
omega_sun = 2 * np.pi / T_trop  # 太陽同期歳差率 [rad/s]

def sso_inclination(h_km, e=0.0):
    """
    太陽同期軌道の傾斜角を計算する。
    h_km: 高度 [km] (円軌道の場合は軌道高度、楕円の場合は平均高度)
    e: 離心率
    戻り値: 傾斜角 [deg]
    """
    a = (R_E + h_km * 1e3)  # 軌道長半径 [m]
    cos_i = -2 * omega_sun * a**(7/2) * (1 - e**2)**2 / (3 * J2 * R_E**2 * np.sqrt(mu))
    if np.abs(cos_i) > 1:
        return np.nan  # 太陽同期条件を満たせない
    return np.degrees(np.arccos(cos_i))

この関数は、前のセクションで導出した設計方程式をそのまま実装したものです。高度と離心率を引数に取り、太陽同期条件を満たす傾斜角を返します。$|\cos i| > 1$ となる場合は太陽同期軌道が存在しないため nan を返します。

高度 vs 傾斜角の関係をプロット

高度200 kmから2000 kmまでの範囲で、太陽同期傾斜角がどのように変化するかを可視化します。

# 高度の範囲を設定
h_range = np.linspace(200, 2000, 500)  # [km]

# 円軌道(e=0)の場合
inc_circular = [sso_inclination(h, e=0.0) for h in h_range]

# 楕円軌道の場合(e=0.01, 0.05)
inc_e001 = [sso_inclination(h, e=0.01) for h in h_range]
inc_e005 = [sso_inclination(h, e=0.05) for h in h_range]

# プロット
fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(h_range, inc_circular, 'c-', linewidth=2, label='e = 0 (circular)')
ax.plot(h_range, inc_e001, '--', color='#4dabf7', linewidth=1.5, label='e = 0.01')
ax.plot(h_range, inc_e005, ':', color='#fab005', linewidth=1.5, label='e = 0.05')

# 代表的な衛星をプロット
satellites = {
    'Landsat-8': (705, 98.2),
    'Sentinel-2': (786, 98.6),
    'ALOS-2': (628, 97.9),
    'Terra': (705, 98.2),
}
for name, (h, i) in satellites.items():
    ax.plot(h, i, 'o', markersize=8, zorder=5)
    ax.annotate(name, (h, i), textcoords="offset points",
                xytext=(10, 5), fontsize=9)

ax.set_xlabel('Altitude [km]', fontsize=12)
ax.set_ylabel('Inclination [deg]', fontsize=12)
ax.set_title('Sun-Synchronous Orbit: Altitude vs Inclination', fontsize=14)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_xlim(200, 2000)
ax.set_ylim(95, 130)
plt.tight_layout()
plt.savefig('sso_altitude_inclination.png', dpi=150)
plt.show()

このグラフから、いくつかの重要な特徴が読み取れます。

  1. 高度と傾斜角はほぼ単調に増加する関係にある — 高度600 kmでは約 $97°$、800 kmでは約 $99°$ と、実用的な範囲では高度が上がるにつれて傾斜角も緩やかに増加します。これはJ2摂動の強さが高度とともに減少するため、より大きな $|\cos i|$ で補う必要があることと一致します。
  2. 離心率が大きくなると、同じ高度でも必要な傾斜角が小さくなる — $e = 0.05$ のカーブは円軌道のカーブよりも下側にあります。近地点で地球に接近するため、J2摂動の効果が増幅されることの反映です。
  3. 代表的な地球観測衛星は理論曲線上にきれいに乗っている — Landsat-8やSentinel-2のプロットが理論曲線にほぼ一致しており、計算の妥当性が確認できます。

RAAN変化率の計算と検証

次に、J2摂動によるRAANの変化率を計算し、太陽同期条件(0.9856 deg/day)と比較する関数を実装します。

def raan_rate(a_m, i_deg, e=0.0):
    """
    J2摂動による昇交点赤経の変化率を計算する。
    a_m: 軌道長半径 [m]
    i_deg: 傾斜角 [deg]
    e: 離心率
    戻り値: dΩ/dt [deg/day]
    """
    i_rad = np.radians(i_deg)
    n = np.sqrt(mu / a_m**3)  # 平均運動 [rad/s]
    dOmega_dt = -1.5 * n * J2 * (R_E / a_m)**2 * np.cos(i_rad) / (1 - e**2)**2
    return np.degrees(dOmega_dt) * 86400  # [deg/day]

# 検証: 高度700kmの太陽同期軌道
h_test = 700  # [km]
i_test = sso_inclination(h_test)
a_test = (R_E + h_test * 1e3)

print(f"高度: {h_test} km")
print(f"太陽同期傾斜角: {i_test:.4f} deg")
print(f"RAAN変化率: {raan_rate(a_test, i_test):.6f} deg/day")
print(f"太陽同期要求値: {360/365.2422:.6f} deg/day")

出力結果を見ると、RAAN変化率が太陽同期要求値(0.985647 deg/day)と非常によく一致していることが確認できます。これは、sso_inclination 関数が正しく太陽同期条件を実装していることの数値的な検証です。

回帰軌道の設計

回帰軌道のパラメータを計算する関数を実装します。回帰日数 $D$ と回帰周回数 $N$ を指定して、対応する軌道高度を反復法で求めます。

def nodal_period(a_m, i_deg, e=0.0):
    """
    ノーダル周期を計算する(J2摂動による補正を含む)。
    """
    i_rad = np.radians(i_deg)
    n = np.sqrt(mu / a_m**3)
    p = a_m * (1 - e**2)

    # J2による平均運動の補正
    delta_n = (3/2) * n * J2 * (R_E / p)**2 * (1 - 1.5 * np.sin(i_rad)**2)

    # J2による近地点引数の変化率
    omega_dot = -(3/2) * n * J2 * (R_E / p)**2 * (2.5 * np.sin(i_rad)**2 - 2) / (1 - e**2)

    # ノーダル周期
    n_eff = n + delta_n + omega_dot  # 簡略化(厳密にはMの変化率を使う)
    # より正確なノーダル周期の計算
    T_kep = 2 * np.pi / n
    T_nodal = T_kep * (1 - (3/2) * J2 * (R_E / p)**2 *
              (3.5 * np.sin(i_rad)**2 - 4 + 2.5 * (1 - e**2)))
    # 近似的なノーダル周期(一次補正)
    T_nodal = T_kep * (1 - 1.5 * J2 * (R_E / a_m)**2 *
              (3 - 2.5 * np.sin(i_rad)**2) / (1 - e**2)**2)
    return T_nodal

def repeat_ground_track(N_rev, D_days, e=0.0, tol=1e-6, max_iter=100):
    """
    回帰軌道の設計: N周回/D日の回帰軌道パラメータを計算する。
    N_rev: 回帰周回数
    D_days: 回帰日数
    e: 離心率
    戻り値: (高度[km], 傾斜角[deg], ノーダル周期[s])
    """
    omega_E = 7.2921159e-5  # 地球の自転角速度 [rad/s]

    # ノーダル日
    T_nodal_day = 2 * np.pi / (omega_E - omega_sun)

    # 必要なノーダル周期
    T_n_required = T_nodal_day * D_days / N_rev

    # 初期推定: ケプラー周期から高度を推定
    a_guess = (mu * (T_n_required / (2 * np.pi))**2)**(1/3)

    a = a_guess
    for _ in range(max_iter):
        h_km = (a - R_E) / 1e3
        i_deg = sso_inclination(h_km, e)
        if np.isnan(i_deg):
            return None, None, None

        T_n = nodal_period(a, i_deg, e)
        # ニュートン法的な補正
        da = (T_n_required - T_n) / T_n * a / 3  # T ∝ a^(3/2)から
        a = a + da
        if abs(T_n - T_n_required) / T_n_required < tol:
            break

    h_km = (a - R_E) / 1e3
    i_deg = sso_inclination(h_km, e)
    T_n = nodal_period(a, i_deg, e)

    return h_km, i_deg, T_n

この実装では、まずケプラー周期から軌道長半径の初期推定値を求め、次にJ2摂動を考慮したノーダル周期と太陽同期傾斜角を反復的に更新しています。回帰条件とSSO条件を同時に満たす解を見つけるのがポイントです。

# 代表的な回帰軌道を計算
missions = [
    ("Landsat-8/9", 233, 16),
    ("Sentinel-2", 143, 10),
    ("ALOS-2", 207, 14),
    ("SPOT-6/7", 385, 26),
]

print(f"{'衛星':15s} {'N/D':>8s} {'高度[km]':>10s} {'傾斜角[deg]':>12s} {'Tn[min]':>10s}")
print("-" * 60)
for name, N, D in missions:
    h, inc, Tn = repeat_ground_track(N, D)
    if h is not None:
        print(f"{name:15s} {N}/{D:>2d}     {h:8.2f}   {inc:10.4f}   {Tn/60:8.3f}")

出力結果から、計算されたパラメータが各衛星の公称値とよく一致していることがわかります。Landsat-8/9は高度約705 km、傾斜角約98.2度と計算され、これは実際の軌道パラメータと合致しています。回帰日数 $D$ が大きい SPOT-6/7(26日)はLandsatと似た高度ですが、地上軌跡の間隔がより密になるという違いがあります。

地上軌跡の可視化

太陽同期回帰軌道の地上軌跡を1回帰周期分プロットしてみましょう。ここでは簡略化のため、球面地球上の円軌道の地上軌跡を計算します。

def ground_track(h_km, i_deg, duration_days, dt=60):
    """
    地上軌跡を計算する(球面地球近似)。
    h_km: 高度 [km]
    i_deg: 傾斜角 [deg]
    duration_days: シミュレーション日数
    dt: 時間刻み [s]
    戻り値: (経度の配列, 緯度の配列) [deg]
    """
    omega_E = 7.2921159e-5  # 地球の自転角速度 [rad/s]
    a = R_E + h_km * 1e3
    n = np.sqrt(mu / a**3)  # 平均運動 [rad/s]
    i_rad = np.radians(i_deg)

    # RAAN変化率(太陽同期なので固定値)
    dOmega = omega_sun  # [rad/s]

    total_time = duration_days * 86400
    t = np.arange(0, total_time, dt)

    # 軌道面内の角度(真近点角の近似として平均近点角を使用)
    theta = n * t  # 近似: 円軌道ではM ≈ θ

    # RAAN
    Omega = dOmega * t  # 初期値Ω0=0

    # 赤道座標系での位置
    # 衛星の緯度
    lat = np.degrees(np.arcsin(np.sin(i_rad) * np.sin(theta)))

    # 衛星の経度(地球固定座標系)
    # 軌道面内の赤経方向成分
    alpha = np.arctan2(np.cos(i_rad) * np.sin(theta), np.cos(theta)) + Omega
    lon = np.degrees(alpha - omega_E * t)

    # 経度を -180〜180 に正規化
    lon = ((lon + 180) % 360) - 180

    return lon, lat

# Landsat-8相当の軌道(16日回帰)
h_landsat = 705  # [km]
i_landsat = 98.2  # [deg]

lon, lat = ground_track(h_landsat, i_landsat, duration_days=1, dt=30)

fig, ax = plt.subplots(figsize=(14, 7))
ax.scatter(lon, lat, c=np.linspace(0, 1, len(lon)), cmap='cool',
           s=0.5, alpha=0.8)
ax.set_xlabel('Longitude [deg]', fontsize=12)
ax.set_ylabel('Latitude [deg]', fontsize=12)
ax.set_title('Ground Track of SSO Satellite (h=705 km, i=98.2°, 1 day)',
             fontsize=14)
ax.set_xlim(-180, 180)
ax.set_ylim(-90, 90)
ax.grid(True, alpha=0.3)
ax.set_aspect('equal')
plt.tight_layout()
plt.savefig('sso_ground_track_1day.png', dpi=150)
plt.show()

このグラフから、太陽同期軌道の地上軌跡の特徴が明確にわかります。

  1. 衛星は南北にジグザグに振動しながら西から東へ移動する — 軌道傾斜角が約98°なので、緯度は約 $\pm 82°$ の範囲を振動します。南北の折り返し点の緯度は $|90° – i|$ ではなく、逆行軌道なので $|i – 90°| = 8.2°$ だけ赤道から離れた高緯度帯までカバーできます。
  2. 連続する地上軌跡は西にずれている — 地球の自転により、各周回の軌跡は約25°(高度705 kmの場合)西にシフトします。16日後にはこの隙間が完全に埋まり、同じパターンが繰り返されます。

16日分の回帰軌跡の可視化

# 16日分(1回帰周期)の地上軌跡
lon_16d, lat_16d = ground_track(h_landsat, i_landsat, duration_days=16, dt=30)

fig, ax = plt.subplots(figsize=(14, 7))
# 色を日ごとに変える
n_points_per_day = int(86400 / 30)
colors = plt.cm.viridis(np.linspace(0, 1, 16))

for day in range(16):
    start = day * n_points_per_day
    end = (day + 1) * n_points_per_day
    if end > len(lon_16d):
        end = len(lon_16d)
    ax.scatter(lon_16d[start:end], lat_16d[start:end],
               c=[colors[day]], s=0.3, alpha=0.6, label=f'Day {day+1}')

ax.set_xlabel('Longitude [deg]', fontsize=12)
ax.set_ylabel('Latitude [deg]', fontsize=12)
ax.set_title('Repeat Ground Track (h=705 km, 16-day cycle, Landsat-like)',
             fontsize=14)
ax.set_xlim(-180, 180)
ax.set_ylim(-90, 90)
ax.grid(True, alpha=0.3)
ax.set_aspect('equal')
plt.tight_layout()
plt.savefig('sso_ground_track_16days.png', dpi=150)
plt.show()

16日分の地上軌跡を重ねたグラフでは、地球全体がほぼ均等に覆われている様子が見て取れます。日ごとに色を変えてプロットしているため、各日の軌跡が少しずつ西にシフトしながら、16日後に初日の軌跡のすぐ隣に戻ってくるパターンが確認できます。これがまさに回帰軌道の設計意図であり、全球を均等にカバーしつつ、16日周期で同じ地点を再訪する仕組みです。

高度・傾斜角の感度解析

最後に、設計パラメータの感度を分析します。目標の太陽同期歳差率からのずれが、どの程度のRAANドリフトを引き起こすかを計算します。

# 傾斜角の微小な誤差がRAANドリフトに与える影響
h_nominal = 705  # [km]
i_nominal = sso_inclination(h_nominal)
a_nominal = R_E + h_nominal * 1e3

# 傾斜角の誤差範囲
di_range = np.linspace(-0.1, 0.1, 100)  # [deg]

# 各誤差でのRAANドリフト(太陽同期要求値との差)
raan_drift = []
for di in di_range:
    rate = raan_rate(a_nominal, i_nominal + di)
    drift = rate - 360 / 365.2422  # [deg/day]
    raan_drift.append(drift * 365.25)  # [deg/year]

fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(di_range, raan_drift, 'c-', linewidth=2)
ax.axhline(y=0, color='gray', linestyle='--', alpha=0.5)
ax.axvline(x=0, color='gray', linestyle='--', alpha=0.5)

ax.set_xlabel('Inclination Error [deg]', fontsize=12)
ax.set_ylabel('RAAN Drift [deg/year]', fontsize=12)
ax.set_title('Sensitivity: Inclination Error → RAAN Drift (h=705 km)', fontsize=14)
ax.grid(True, alpha=0.3)

# 典型的な制御精度の範囲を表示
ax.axvspan(-0.02, 0.02, alpha=0.15, color='cyan', label='Typical control accuracy')
ax.legend(fontsize=10)

plt.tight_layout()
plt.savefig('sso_sensitivity.png', dpi=150)
plt.show()

print(f"公称傾斜角: {i_nominal:.4f} deg")
print(f"傾斜角 0.01° の誤差によるRAANドリフト: "
      f"{np.interp(0.01, di_range, raan_drift):.3f} deg/year")
print(f"傾斜角 0.05° の誤差によるRAANドリフト: "
      f"{np.interp(0.05, di_range, raan_drift):.3f} deg/year")

感度解析のグラフから、太陽同期軌道の維持に必要な精度が明らかになります。

  1. 傾斜角 $0.01°$ の誤差で、RAANは年間約 $0.3°$ ドリフトする — これは比較的小さな値ですが、10年間では $3°$ に達し、LTANにして約12分のずれになります。
  2. 傾斜角 $0.05°$ の誤差では年間約 $1.5°$ のドリフト — これは無視できない量で、数年以内にミッション要求を逸脱する可能性があります。
  3. 感度はほぼ線形 — この高度範囲では、傾斜角誤差とRAANドリフトの関係がほぼ直線的です。

このため、実際の太陽同期軌道衛星は定期的な軌道制御(傾斜角調整マヌーバ)を行って、太陽同期条件からのドリフトを許容範囲内に収めています。

代表的な太陽同期軌道衛星

Landsat計画

Landsat計画は1972年に始まった地球観測プログラムで、現在はLandsat-8(2013年打ち上げ)とLandsat-9(2021年打ち上げ)が運用中です。両衛星は同じ軌道面内を8日ずらして飛行し、16日回帰軌道の全球カバレッジを8日ごとに完成させています。

  • 高度: 705 km
  • 傾斜角: 98.2°
  • 回帰周期: 16日(233周回)
  • LTAN: 約10:00

Landsatのデータは無料で公開されており、50年以上にわたる地球表面の変化記録として、気候変動研究や土地利用変化のモニタリングに不可欠な資源となっています。

Sentinel-2

ESA(欧州宇宙機関)のCopernicus計画で運用される光学観測衛星です。

  • 高度: 786 km
  • 傾斜角: 98.6°
  • 回帰周期: 10日(143周回)
  • LTAN: 約10:30

Landsatよりやや高い軌道で、回帰周期は10日と短く設定されています。Sentinel-2A と Sentinel-2B の2機体制で5日ごとの再訪が実現されています。

だいち2号(ALOS-2)

JAXAが運用するLバンド合成開口レーダー(SAR)衛星です。

  • 高度: 628 km
  • 傾斜角: 97.9°
  • 回帰周期: 14日(207周回)
  • LTAN: 約12:00(降交点)

SARは自ら電波を照射して観測するため、昼夜・天候に関係なくデータを取得できます。地殻変動の検出(InSAR)や森林バイオマスの推定など、光学衛星では困難な観測に威力を発揮します。

各衛星の軌道パラメータ比較

これらの衛星はいずれも太陽同期軌道の基本原理を共有しながら、ミッション目的に応じて高度・回帰周期・LTANを最適化しています。高い空間分解能が必要なら低軌道(短い周期)、広域を頻繁にカバーしたいなら短い回帰日数、特定の照明条件が必要ならLTANの調整 — というように、太陽同期軌道の設計は多くのトレードオフの結果として決定されます。

まとめ

本記事では、太陽同期軌道の理論と設計について、基本原理から数学的導出、Pythonによる実装まで解説しました。

  • 太陽同期軌道の原理: 地球のJ2摂動(赤道膨らみによる重力非対称性)を利用して、昇交点赤経を年間 $360°$(1日あたり $0.9856°$)東向きに歳差させることで、衛星と太陽の幾何学的関係を一年を通じて一定に保つ
  • 設計方程式: $\cos i = -\frac{2 \dot{\Omega}_{\text{sun}} a^{7/2}(1-e^2)^2}{3 J_2 R_E^2 \sqrt{\mu}}$ により、高度と離心率から太陽同期傾斜角が一意に決まる。実用的なSSO高度(400〜1500 km)では傾斜角は $96°$ 〜 $104°$ の逆行軌道となる
  • 回帰軌道: 回帰日数 $D$ と回帰周回数 $N$ の組み合わせで、定期的に同じ地上軌跡を通る軌道を設計できる。Landsatは16日233周回、Sentinel-2は10日143周回
  • 地方太陽時(LST): 衛星が赤道を通過する時刻を制御することで、観測に最適な太陽照明条件を確保する。光学観測では午前10時〜10時半が多用される
  • 感度と軌道維持: 傾斜角のわずかな誤差がRAANドリフトを引き起こすため、定期的な軌道制御が必要

太陽同期軌道は、J2摂動という「自然の力」を巧みに利用した軌道設計の傑作です。一見すると衛星運用の障害となるJ2摂動を、逆にミッション目的の達成に活用するという発想は、軌道力学の美しさを象徴しています。

次のステップとして、以下の記事も参考にしてください。