静止軌道（GEO）の通信衛星は、見かけ上「地球から見て止まっている」優等生ですが、内部では絶え間ない姿勢調整が続いています。打ち上げ直後にきれいに地球指向していたアンテナは、数日放っておくとビームが衛星サービスエリアからずれ始めます。原因は重力でも磁場でもなく、太陽光がパドルやバスに当たって発生するわずか数μNmのトルク——太陽光圧（Solar Radiation Pressure, SRP）トルクです。
地表に立つ人間にとって、太陽光が「押す力」を持つというのは奇妙に聞こえます。手のひらに当たる真夏の太陽光ですら、押される感覚は皆無です。実際、地球付近での太陽光圧は $4.56\,\mu\mathrm{N/m^2}$ と極めて小さく、$10\,\mathrm{m}^2$ のパドルにかかる力でも $45\,\mu\mathrm{N}$、これは1円玉の重さの 1/20 にも満たない量です。ところが GEO 衛星は 15年間 真空中でこの力を浴び続け、しかも重心（CoM）と光圧中心（CoP）のオフセットによって発生するトルクは、リアクションホイールのモメンタムを少しずつ累積させ、最終的に飽和（saturation）させてしまいます。飽和したホイールは姿勢制御の能力を失い、ミッション継続にはスラスタや磁気トルカによる モメンタムアンローディング が必須となります。
この見えない外乱を理解することは、現代の通信衛星設計の中心課題のひとつです。地球観測衛星のポインティング精度、宇宙望遠鏡の長時間積分露光、月・惑星探査機の数年スパンのクルーズフェーズ——いずれも、SRPトルクの定量化と長期予算管理ができなければ成り立ちません。本記事では、SRPトルクの物理から重心-光圧中心オフセットの幾何、季節変動、ホイール飽和、アンローディング戦略までを順を追って解説し、Pythonで15年間のミッションプロファイルをシミュレーションします。
なお、太陽光圧そのものを「推進力」として積極的に利用するソーラーセイルの話題は別記事に譲り、本記事では一貫して「避けがたい外乱」としての SRP に焦点を絞ります。
本記事の内容
- 太陽光圧の物理（光子モーメンタム流束、反射係数、温度依存）
- 重心 $\bm{r}_{cm}$ と光圧中心 $\bm{r}_{cp}$ のオフセットから生じるSRPトルク $\bm{\tau} = (\bm{r}_{cp}-\bm{r}_{cm}) \times \bm{F}_{srp}$
- パドル展開時の非対称・アンテナ展開後の影によるオフセット増大
- 季節変動（地球公転による太陽方向の年周変化）
- リアクションホイールのモメンタム蓄積と飽和
- 磁気トルカ（LEO）／スラスタ（GEO）／重力傾度トルクによるアンローディング
- SRPによる軌道力学への二次効果（離心率の長期成長）
- Python による15年運用シミュレーションと姿勢誤差予算
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
「光が押す」とはどういうことか
真空中の太陽光は、波としての電磁波であると同時に、運動量 $p = E/c$ を持つ光子の流れでもあります。衛星表面が日射を吸収すると、光子の運動量がそのまま衛星に渡り、その方向に微小な力が発生します。さらに表面が鏡のように反射すれば、光子は来た方向に跳ね返り、反作用としてさらに大きな運動量が衛星に伝わります。「太陽光が押す」とは、毎秒地球軌道上の $1\,\mathrm{m}^2$ あたり $4.56 \times 10^{17}$ 個もの光子が連続的に衝突して、その運動量を表面に渡し続けている状態のことです。
地球軌道上での太陽光圧 $P_{\mathrm{srp}}$ は、太陽放射照度 $S_0 \approx 1361\,\mathrm{W/m^2}$（太陽定数）から次のように見積もれます。光子1個あたりのモーメンタムは $E/c$ なので、単位面積・単位時間あたりに到達する運動量、すなわち光圧は
$$ P_{\mathrm{srp}} = \frac{S_0}{c} = \frac{1361}{2.998 \times 10^8} \approx 4.56 \times 10^{-6}\,\mathrm{N/m^2} $$
となります。これは「1m² あたり 4.56μN の押す力」を意味します。たかが μN と思うかもしれませんが、これが衛星の数 m² の表面に、止まることなく、15年間ずっと作用し続けるという点が重要です。総モーメンタムは時間積分で効いてくるので、長期運用ではこの「ちりも積もれば」の効果が支配的になります。
太陽からの距離が $r$ AU の場所では、放射照度は $1/r^2$ で減るため、光圧も
$$ P_{\mathrm{srp}}(r) = P_{\mathrm{srp},0} \cdot \frac{1}{r^2} $$
となります。火星軌道（$r \approx 1.52$ AU）では地球の $1/2.3$、木星軌道（$r \approx 5.2$ AU）では $1/27$ に減衰します。逆に水星近傍では $r \approx 0.39$ AU で約 6.6 倍にもなり、ベピコロンボなどの探査機では SRP 外乱が桁違いに厳しくなります。本記事では一貫して地球公転軌道近傍（$r \simeq 1$ AU）を想定します。
ここまで「光圧」というスカラー的な圧力を見てきましたが、衛星に作用するのは方向性を持った力、さらにそこから幾何学的に派生するトルクです。次のセクションで、衛星表面に作用する SRP 力を厳密に書き下します。
反射係数と表面要素にかかる力
一般的な表面のSRP力モデル
衛星表面は単純な完全吸収体でも完全鏡面反射体でもなく、両者の混合に加えて拡散反射（Lambert 反射）まで含む複雑な光学特性を持ちます。微小表面要素 $dA$、法線ベクトル $\hat{\bm{n}}$、太陽から表面への単位ベクトル $\hat{\bm{s}}$（太陽 → 衛星向き）を考えます。太陽光が表面に当たる入射角を $\theta$ とすると、$\cos\theta = -\hat{\bm{s}}\cdot\hat{\bm{n}}$ です（$\hat{\bm{s}}$ が来る向き、$\hat{\bm{n}}$ が出る向きなので符号を反転して内積を取る）。
工学的によく使われる CR モデル（Coefficient of Reflectivity） に基づくと、表面要素 $dA$ に作用する力 $d\bm{F}$ は次のように書けます。
$$ d\bm{F} = -P_{\mathrm{srp}} \cos\theta \left\{ (1 – \rho_s)\hat{\bm{s}} + 2\left(\rho_s \cos\theta + \frac{1}{3}\rho_d\right)\hat{\bm{n}} \right\} dA $$
ここで $\rho_s$ は鏡面反射率（specular reflectivity）、$\rho_d$ は拡散反射率（diffuse reflectivity）、$1 – \rho_s – \rho_d$ が吸収率です。第1項 $(1-\rho_s)\hat{\bm{s}}$ は吸収と拡散反射に伴う「太陽方向へ押される」成分、第2項 $\hat{\bm{n}}$ 方向の項は法線方向の反射圧成分です。$1/3$ という因子は、Lambert 拡散面から半球状に再放射される光子のモーメンタム成分の平均値から来ます。
完全吸収体（$\rho_s = \rho_d = 0$）ならば
$$ d\bm{F}_{\mathrm{abs}} = -P_{\mathrm{srp}} \cos\theta\, \hat{\bm{s}}\, dA $$
で、力は単に「太陽から衛星へ」の方向です。完全鏡面反射体（$\rho_s = 1, \rho_d = 0$）ならば
$$ d\bm{F}_{\mathrm{spec}} = -2 P_{\mathrm{srp}} \cos^2\theta\, \hat{\bm{n}}\, dA $$
で、入射角が浅いほど（$\cos\theta \to 1$）力は2倍に増えますが、入射角が90°に近づくと急速にゼロに減衰します。これがソーラーセイルの「セイルを傾けて推力方向を制御する」物理的根拠ですが、衛星バスでは「斜めから日が当たると思いがけない方向にトルクが発生する」厄介な性質として現れます。
反射係数の温度依存
衛星熱制御材料（MLI、テフロン、銀蒸着 Kapton など）の反射係数は、温度や経年劣化で大きく変動します。実例として銀蒸着 Teflon は新品時に $\rho_s \simeq 0.85$ あったものが、紫外線・原子状酸素・宇宙放射線によって 10 年後には $0.5$ 以下まで劣化することが知られています。これは、SRP トルクの長期的にゆっくり変化する成分として姿勢制御系に効いてきます。設計時には材料の「end-of-life（EOL）光学特性」を用いて最悪ケースを評価するのが標準的です。
温度依存については、表面温度 $T$ が上昇すると分子レベルでの光散乱モードが変わり、$\rho_d$ がわずかに増加する一方 $\rho_s$ が減少する傾向があります。GEO の食（イクリプス）期前後では衛星表面温度が $-100^\circ\mathrm{C}$ から $+80^\circ\mathrm{C}$ まで激変するため、SRP トルクも食前後で系統的な変動を見せます。
ここまで「面要素に作用する力」を導いたので、次は衛星全体を構成する複数面の積分から、衛星に作用するトルクを構成します。
重心-光圧中心オフセットとトルクの正体
全体力とトルクの定義
衛星表面 $\Sigma$ 全体での総力 $\bm{F}_{\mathrm{srp}}$ と、重心まわりのトルク $\bm{\tau}_{\mathrm{srp}}$ は、表面要素の力 $d\bm{F}(\bm{r})$ を積分して
$$ \bm{F}_{\mathrm{srp}} = \int_\Sigma d\bm{F}(\bm{r}), \quad \bm{\tau}_{\mathrm{srp}} = \int_\Sigma (\bm{r} – \bm{r}_{cm}) \times d\bm{F}(\bm{r}) $$
で定義されます。$\bm{r}_{cm}$ は衛星の重心位置（質量中心）です。トルクが「ベクトルの外積」の形で現れるのは、力の作用点が重心からどれだけずれているかで回転モーメントが決まるからです。
ここで、SRP 力に対する等価な作用点として 光圧中心（Center of Pressure, CoP） $\bm{r}_{cp}$ を定義します。これは、全表面に分布した SRP 力 $d\bm{F}(\bm{r})$ を 1 点に集約したときに、合力 $\bm{F}_{\mathrm{srp}}$ と全トルク $\bm{\tau}_{\mathrm{srp}}$ を両立する点として、
$$ \bm{r}_{cp} \times \bm{F}_{\mathrm{srp}} = \int_\Sigma \bm{r} \times d\bm{F}(\bm{r}) $$
を満たすように定義されます（厳密には CoP は一意に決まらない場合があり、力の方向によって変動する「仮想的な作用点」になります）。この定義を使うと、重心まわりのトルクは
$$ \bm{\tau}_{\mathrm{srp}} = (\bm{r}_{cp} – \bm{r}_{cm}) \times \bm{F}_{\mathrm{srp}} $$
という非常にコンパクトな形にまとまります。$\bm{r}_{cp}$ と $\bm{r}_{cm}$ がぴったり一致すればトルクはゼロですが、現実の衛星でそれが起こることはまずありません。
なぜオフセットが避けられないのか
衛星設計者は重心と光圧中心を可能な限り一致させたいと願いますが、いくつかの避けがたい理由で必ずずれが発生します。
- 太陽電池パドルの非対称展開: 大型 GEO 通信衛星では片側のパドルが数 m〜十数 m に達します。両側パドルの寸法・反射率を完全に同一にすることは製造上不可能で、$\rho_s$ の左右差わずか 1% で数 cm の CoP オフセットが生じます。
- アンテナ展開後の影: 展開型大型反射鏡が衛星バスの一部に影を落とすと、その面の SRP 力が消えて CoP が影と反対側にずれます。
- 燃料消費による CoM 移動: 化学スラスタの推進薬を消費するにつれて重心が移動します。15年運用では数 cm 以上の CoM 移動が普通です。
- 熱変形: 太陽指向側と影側で温度差 $\Delta T \sim 100\,\mathrm{K}$ が生じ、CFRP 構造でも数 mm の熱変形が起こります。これが CoP の位置を周期的に変えます。
- MLI のしわや経年たるみ: 多層断熱材は完全な平面ではなく、微小なしわや経年劣化で局所的に反射特性が変化し、CoP を移動させます。
実機では、これらの効果を統合した結果、CoP オフセット $|\bm{r}_{cp} – \bm{r}_{cm}|$ はおよそ $5\,\mathrm{cm}$ から $30\,\mathrm{cm}$ 程度になることが多いです。
トルクの大きさの感覚
具体的な数値を入れてみます。GEO 通信衛星のバス断面積 $A = 15\,\mathrm{m}^2$、平均反射係数 $C_R \simeq 1.5$、CoP-CoM オフセット $|\Delta\bm{r}| = 0.1\,\mathrm{m}$ とすると、$\bm{F}_{\mathrm{srp}}$ の大きさは
$$ |\bm{F}_{\mathrm{srp}}| \approx C_R P_{\mathrm{srp}} A = 1.5 \times 4.56\times 10^{-6} \times 15 \approx 1.0 \times 10^{-4}\,\mathrm{N} $$
すなわち $100\,\mu\mathrm{N}$ 程度です。これに $10\,\mathrm{cm}$ のレバーアームをかけると
$$ |\bm{\tau}_{\mathrm{srp}}| \approx 1.0 \times 10^{-4} \times 0.1 = 1.0 \times 10^{-5}\,\mathrm{N\cdot m} = 10\,\mu\mathrm{N\cdot m} $$
となります。わずか 10μN·m。これがリアクションホイールが日々受け続ける外乱の典型値です。
ここまでで「瞬間の」SRP トルクの大きさが見えました。しかしこれだけでは姿勢制御系の真の難しさは見えません。地球の公転に伴って太陽方向が年周変化することで、SRP トルクは長期的に積分されて累積する——次のセクションでその構造を解きほぐします。
季節変動と長期累積トルク
太陽方向の年周変動
地球が太陽の周りを公転する以上、衛星から見た太陽方向 $\hat{\bm{s}}$ は1年で360°回転します。慣性座標系（J2000系）で太陽方向は、近似的に
$$ \hat{\bm{s}}(t) = \begin{pmatrix} \cos\lambda_\odot(t) \\ \sin\lambda_\odot(t)\cos\varepsilon \\ \sin\lambda_\odot(t)\sin\varepsilon \end{pmatrix} $$
と書けます。ここで $\lambda_\odot$ は黄経（年周で 0→360°）、$\varepsilon \approx 23.44^\circ$ は地軸傾斜です。
GEO 衛星では「地球指向（Z軸を地球向き）」と「太陽指向（パドル法線を太陽向き）」を両立させる必要があり、衛星本体は地球まわりに 1 日周期で回転、パドルは衛星本体に対し 1 日周期で逆回転して常に太陽を向くという二重回転構造を取ります。この座標系の中で SRP トルクを評価すると、トルクの1日成分（衛星座標系での1日周期）と年周成分（年周期）が重なった複雑な時系列になります。
累積モーメンタムの式
リアクションホイールに蓄えられるモーメンタム $\bm{h}(t)$ は、外乱トルクを時間積分したものです。
$$ \bm{h}(t) = \bm{h}_0 + \int_0^t \bm{\tau}_{\mathrm{srp}}(t’)\, dt’ $$
ここで重要なのは、$\bm{\tau}_{\mathrm{srp}}$ が周期的成分と DC（直流）成分に分けられることです。周期成分は積分しても増えずキャンセルしますが、DC 成分は積分するごとに線形に増えていきます。年間平均トルク $\bar{\bm{\tau}}$ の大きさが $1\,\mu\mathrm{N\cdot m}$ 程度残ると、
$$ |\Delta \bm{h}|_{\mathrm{year}} = \bar{\tau} \times T_{\mathrm{year}} = 10^{-6} \times 3.16\times 10^7 \approx 31.6\,\mathrm{N\cdot m\cdot s} $$
となり、1年あたり 30 N·m·s 級のモーメンタム蓄積が発生します。一般的なリアクションホイールの最大モーメンタム容量は数十〜数百 N·m·s なので、数年〜十年でホイールが飽和する計算です。
これが「GEO 衛星では数日単位で姿勢が崩れるが、その背後に15年スケールのモーメンタム累積という見えない時計が動いている」という、本記事冒頭で述べた現象の正体です。
ここで一つの自然な疑問が浮かびます——どうやってこの累積モーメンタムを「降ろす（アンロード）」のか？ 次のセクションで、軌道高度に応じた3つの代表的なアンローディング戦略を見ていきます。
モメンタムアンローディング戦略
LEO: 磁気トルカ
低軌道（LEO）衛星では、地球磁場が十分強いため磁気トルカ（magnetorquer） が主役になります。衛星に搭載した電磁石コイルが磁気モーメント $\bm{m}$ を生成すると、地球磁場 $\bm{B}$ との相互作用で
$$ \bm{\tau}_{\mathrm{mtq}} = \bm{m} \times \bm{B} $$
のトルクが発生します。LEO 高度では $|\bm{B}| \sim 25\,\mu\mathrm{T}$、典型的なトルカが $\bm{m} \sim 10\,\mathrm{A\cdot m^2}$ を出せるので、最大 $\bm{\tau} \sim 250\,\mu\mathrm{N\cdot m}$ が得られます。
ただし磁気トルカには制約があります。$\bm{B}$ に平行な成分のトルクは作れないため、3軸すべての制御は不可能で、衛星が軌道を1周する間に少しずつ各軸へモーメンタムを降ろします。標準的なアンローディング則は B-dot 則やそのモーメンタム版で、
$$ \bm{m} = -k\, \bm{h}_{\mathrm{rw}} \times \bm{B} / |\bm{B}|^2 $$
のように設計します。$k$ は正のゲイン。これにより $\bm{\tau} = \bm{m} \times \bm{B}$ が概ね $\bm{h}_{\mathrm{rw}}$ と反対向きを持ち、ホイールモーメンタムを少しずつ放出できます。
GEO: スラスタ
GEO 高度では地球磁場が $|\bm{B}| \sim 0.1\,\mu\mathrm{T}$ と LEO の 1/250 まで弱まり、磁気トルカは事実上使えません。代わりに化学スラスタまたは電気推進スラスタでアンローディングを行います。スラスタを CoM からのレバーアーム $\bm{r}_{\mathrm{thr}}$ だけずらした位置で噴射すると、
$$ \bm{\tau}_{\mathrm{thr}} = \bm{r}_{\mathrm{thr}} \times \bm{F}_{\mathrm{thr}} $$
のトルクが得られます。$\bm{F}_{\mathrm{thr}} \sim 10\,\mathrm{mN}$（電気推進）から $10\,\mathrm{N}$（化学）、レバーアーム $\sim 1\,\mathrm{m}$ で、トルクは $10\,\mathrm{mN\cdot m}$ から $10\,\mathrm{N\cdot m}$ と SRP 外乱より桁違いに大きく取れます。
GEO 通信衛星では、南北位置維持機動（NSSK: North-South Station Keeping） のスラスタ噴射と同時にアンローディングを行うのが標準的です。NSSK は黄道面と赤道面のずれ（傾斜角ドリフト $\sim 0.85^\circ/\mathrm{year}$）を補正するため、年に複数回スラスタを焚きますが、その噴射方向と作用点を巧妙に設計することで、軌道制御とモーメンタムダンプを一石二鳥でこなします。
重力傾度トルクとの併用
低軌道では重力傾度トルクもアンローディング補助に使えます。地球からの重力勾配により、衛星の主慣性軸が地球-衛星ベクトルに整列しようとする復元トルクが発生します。これを利用して受動的に長軸を地球に向け、SRP の累積を相殺する設計（GOCE 衛星など）もあります。
アンローディング制御則のまとめ
アンローディングは「ホイールモーメンタムをゼロ近傍に保つ」サブ制御として動作します。代表的な PI（比例-積分）型のアンロード則は
$$ \bm{u}_{\mathrm{unload}}(t) = K_p \bm{h}_{\mathrm{rw}}(t) + K_i \int_0^t \bm{h}_{\mathrm{rw}}(t’)\, dt’ $$
で、$\bm{u}_{\mathrm{unload}}$ がトルカ／スラスタへの指令値となります。比例項で「現在の蓄積モーメンタム」、積分項で「長期バイアス」を打ち消す構造です。
ここまでで「累積モーメンタムをどう降ろすか」のメカニズムが分かりました。次に、これらが軌道力学にも影響することを見ます。SRP は姿勢だけでなく軌道にもじわじわと作用しているのです。
SRPによる軌道力学への二次効果
離心率の長期成長
SRP の力は重心に作用する成分（並進力）と、重心まわりのトルク（回転モーメント）に分解できます。後者が姿勢を擾乱するのは前章までで述べたとおりですが、前者は衛星の軌道そのものをじわじわ歪めます。
GEO 衛星に対する SRP 並進力は、太陽方向に向かう $100\,\mu\mathrm{N}$ オーダーの定常的な力です。衛星質量 $m = 3000\,\mathrm{kg}$ とすると加速度は
$$ a_{\mathrm{srp}} = F/m \approx 10^{-4}/3000 \approx 3.3\times 10^{-8}\,\mathrm{m/s^2} $$
となり、地球重力の $0.225\,\mathrm{m/s^2}$（GEO 高度）に比べ約 $10^{-7}$ の摂動です。一見無視できる大きさですが、これが年単位で積分されると軌道要素に明確な変動を生みます。
特に顕著なのは離心率 $e$ の長期成長です。SRP の方向は年周で回転するため、Gauss の摂動方程式を年周平均すると、離心率ベクトル $\bm{e}$ が円を描いてゆっくり成長する離心率ベクトルの長周期項が現れます。年周平均された SRP 加速度を $a_{\mathrm{srp}}$、太陽の角速度を $n_\odot = 2\pi/(1\,\mathrm{year})$ とすると、離心率増加率は
$$ \dot{e}_{\mathrm{srp}} \sim \frac{3}{2}\frac{a_{\mathrm{srp}}}{v_{\mathrm{geo}}}\sin(\Omega_\odot t) $$
オーダーで、GEO 衛星では年間で $\Delta e \sim 10^{-4}$ の振幅を持ちます。これは赤道面内で衛星が東西に $\pm 40\,\mathrm{km}$ ふらつく量に相当し、通信衛星のサービスエリアからずれる原因のひとつとなります。東西位置維持（E-W Station Keeping） がスラスタによって年に複数回必要となる理由のかなりの部分は、この SRP 起因の離心率ドリフトです。
軌道-姿勢連成
さらに厄介なのは、SRP トルクと SRP 並進力が同じ太陽方向に依存するため、両者が完全に独立に扱えないことです。CoP オフセットが大きい衛星では、姿勢誤差が並進力の方向をわずかに変え、それが軌道に二次的影響を及ぼし、軌道誤差がさらに太陽指向誤差に跳ね返る——という連成ループが形成されます。詳細は軌道-姿勢統合シミュレータが必要ですが、本記事のシミュレーションでも単純化した形で扱います。
理論パートはここまでです。次は Python で 15 年運用シミュレーションを実装し、ここまで議論してきたモーメンタム累積・アンローディング・姿勢誤差予算を具体的に可視化します。
Pythonでの実装
SRPトルクモデル
まず単純化した衛星形状（直方体バス + 2枚のパドル）の SRP トルクを計算する関数を作ります。後で時系列シミュレーションに組み込むため、関数化しておきます。
import numpy as np
# 物理定数
SRP_AU = 4.56e-6 # N/m^2 地球軌道での太陽光圧
YEAR = 365.25 * 86400.0
DAY = 86400.0
OMEGA_ORBIT = 2*np.pi / YEAR
EPS_AXIAL = np.deg2rad(23.44) # 黄道傾斜角
def sun_unit_vector(t):
"""慣性系での太陽 → 衛星の単位ベクトル（時刻 t [s]）"""
lam = OMEGA_ORBIT * t
s = np.array([np.cos(lam),
np.sin(lam)*np.cos(EPS_AXIAL),
np.sin(lam)*np.sin(EPS_AXIAL)])
return -s / np.linalg.norm(s) # 衛星 → 太陽の符号を反転
def srp_force_torque(s_hat, surfaces, r_cm):
"""衛星表面の集合に対する総力 F とトルク tau を返す
s_hat: 太陽 → 衛星方向単位ベクトル
surfaces: list of dict(area, n_hat[面法線], r[面中心位置], rho_s, rho_d)
r_cm: 重心位置 (3-vector)
"""
F = np.zeros(3)
tau = np.zeros(3)
for surf in surfaces:
n = surf['n_hat']
cos_theta = -np.dot(s_hat, n)
if cos_theta <= 0: # 裏面は寄与しない（影の単純化）
continue
rho_s, rho_d = surf['rho_s'], surf['rho_d']
dF = -SRP_AU * cos_theta * (
(1 - rho_s)*s_hat
+ 2*(rho_s*cos_theta + rho_d/3.0)*n) * surf['area']
F += dF
tau += np.cross(surf['r'] - r_cm, dF)
return F, tau
このコードは前述の表面要素 SRP モデルを素直に実装したものです。cos_theta <= 0 の条件は、太陽光が裏側から当たる面（自己遮蔽）を簡易的に除外しています。実機では遮蔽ジオメトリは複雑ですが、本記事の範囲ではこの単純化で十分です。r_cm を引数で受けるのは、後に重心が燃料消費でずれていく効果を扱えるようにするためです。
GEO通信衛星モデル
次に、典型的な GEO 通信衛星の幾何を定義します。バス（直方体）の 6 面と、左右の太陽電池パドル（薄板）を表現します。
def make_geo_satellite(cm_offset=np.array([0.05, 0.02, 0.0])):
"""GEO 通信衛星（バス2.5×2.5×2.5m, パドル左右各 4×10m）
cm_offset: 公称重心からのずれ [m]（CoPオフセットの主因）
"""
surfaces = []
# バス6面
bus_size = 2.5
a_bus = bus_size**2
for axis, sign in [(0,+1),(0,-1),(1,+1),(1,-1),(2,+1),(2,-1)]:
n = np.zeros(3); n[axis] = sign
r = (bus_size/2.0)*n # 面中心
surfaces.append(dict(area=a_bus, n_hat=n, r=r,
rho_s=0.3, rho_d=0.2))
# 太陽電池パドル左右
paddle_area = 4.0 * 10.0
paddle_x = bus_size/2.0 + 5.0
# 左パドル（太陽追尾を Y軸法線として静的近似）
surfaces.append(dict(area=paddle_area, n_hat=np.array([0,+1,0]),
r=np.array([+paddle_x, 0, 0]),
rho_s=0.08, rho_d=0.15))
# 右パドル
surfaces.append(dict(area=paddle_area, n_hat=np.array([0,+1,0]),
r=np.array([-paddle_x, 0, 0]),
rho_s=0.08, rho_d=0.15))
return surfaces, cm_offset
ここでパドルの法線を +Y 方向に固定したのは「太陽指向制御によりパドル面が常に太陽の方を向いている」という GEO 衛星の運用を簡略化して再現するためです。実際にはパドルが衛星本体に対し1日周期で回転しますが、平均的な SRP トルクの大きさを見るにはこの近似で十分です。cm_offset をデフォルトで $(5\,\mathrm{cm}, 2\,\mathrm{cm}, 0)$ にしているのは、現実の衛星で典型的に観測される重心-光圧中心オフセットを再現するためです。
瞬間トルクの可視化
太陽方向を $0\,^\circ$ から $360^\circ$ まで回したとき、SRP トルクの3軸成分がどう変化するかをプロットします。
import matplotlib.pyplot as plt
surfaces, r_cm = make_geo_satellite()
angles = np.linspace(0, 2*np.pi, 361)
F_all = np.zeros((len(angles), 3))
tau_all = np.zeros((len(angles), 3))
for i, a in enumerate(angles):
s_hat = np.array([np.cos(a), np.sin(a), 0.0])
F_all[i], tau_all[i] = srp_force_torque(s_hat, surfaces, r_cm)
fig, axes = plt.subplots(2, 1, figsize=(10, 7), sharex=True)
labels = ['x', 'y', 'z']
for k in range(3):
axes[0].plot(np.rad2deg(angles), F_all[:,k]*1e6, label=f'F_{labels[k]}')
axes[1].plot(np.rad2deg(angles), tau_all[:,k]*1e6,
label=fr'$\tau_{labels[k]}$')
axes[0].set_ylabel('Force [μN]')
axes[1].set_ylabel('Torque [μN·m]')
axes[1].set_xlabel('Sun direction angle [deg]')
axes[0].set_title('SRP force and torque vs sun direction (GEO bus + panels)')
for ax in axes:
ax.grid(alpha=0.3); ax.legend(loc='upper right')
plt.tight_layout()
plt.savefig('srp_torque_vs_sun.png', dpi=150, bbox_inches='tight')
plt.show()
上のグラフから、SRP 力とトルクの両方が太陽方向の周期関数として大きく変動することが読み取れます。力 $\bm{F}$ は太陽方向に約 $100\,\mu\mathrm{N}$ で押されるベクトルとして $\sin/\cos$ 形で振動します。一方トルク $\bm{\tau}$ は CoP-CoM オフセットの方向に依存して、X軸まわり・Y軸まわりに $5\,\mu\mathrm{N\cdot m}$ 程度の振幅で振動します。注目すべきは、Y軸成分にゼロ点をまたぐ振動があるのに対し、Z軸成分には軌道平均で DC バイアスが現れる場合があるという点です。この DC 成分こそが、長期モーメンタム累積の真の駆動力です。
15年間のモーメンタム蓄積シミュレーション
次に、太陽方向が1年で回転する効果を含めて、15年間の累積モーメンタム $\bm{h}(t)$ を時間発展させます。
import numpy as np
import matplotlib.pyplot as plt
surfaces, r_cm = make_geo_satellite(cm_offset=np.array([0.08, 0.03, 0.0]))
dt = DAY * 0.5 # 12時間刻み
T_end = 15 * YEAR
ts = np.arange(0, T_end, dt)
h = np.zeros(3) # 累積モーメンタム
h_hist = np.zeros((len(ts), 3))
tau_hist = np.zeros((len(ts), 3))
for i, t in enumerate(ts):
s_hat = sun_unit_vector(t)
_, tau = srp_force_torque(s_hat, surfaces, r_cm)
h = h + tau * dt
h_hist[i] = h
tau_hist[i] = tau
# 飽和ライン
H_SAT = 50.0 # N·m·s（典型的なリアクションホイール容量）
fig, axes = plt.subplots(2, 1, figsize=(11, 7), sharex=True)
years = ts / YEAR
for k, lab in enumerate(['x','y','z']):
axes[0].plot(years, tau_hist[:,k]*1e6, label=fr'$\tau_{lab}$', lw=0.8)
axes[1].plot(years, h_hist[:,k], label=fr'$h_{lab}$')
axes[0].set_ylabel('SRP torque [μN·m]')
axes[1].set_ylabel('Accumulated momentum [N·m·s]')
axes[1].axhline(+H_SAT, color='k', ls='--', alpha=0.5,
label=f'±{H_SAT} N·m·s (saturation)')
axes[1].axhline(-H_SAT, color='k', ls='--', alpha=0.5)
axes[1].set_xlabel('Time [years]')
axes[0].set_title('SRP torque (top) and accumulated reaction wheel momentum (bottom)')
for ax in axes:
ax.grid(alpha=0.3); ax.legend(loc='upper right', fontsize=9)
plt.tight_layout()
plt.savefig('momentum_15yr_no_unload.png', dpi=150, bbox_inches='tight')
plt.show()
print(f"15年後の |h| = {np.linalg.norm(h):.1f} N·m·s")
print(f"飽和到達時刻（最初に |h|>50 を超える）: ",
end='')
exceeded = np.argmax(np.linalg.norm(h_hist, axis=1) > H_SAT)
if exceeded > 0:
print(f"{ts[exceeded]/YEAR:.2f} 年")
else:
print("飽和せず")
このシミュレーションの結果、トルク $\bm{\tau}$ は周期的に振動するものの、軌道平均で残るバイアス成分により累積モーメンタム $\bm{h}$ は時間に対し線形に増大することが見て取れます。デフォルトのパラメータ（CoM オフセット 8cm）では、おおむね 2〜3年で 50 N·m·s の飽和ラインを超えます。GEO 通信衛星の運用寿命15年に対し、アンローディングが無ければ早期にホイールが使い物にならなくなることがこれで定量的に分かります。1年あたりの累積モーメンタムが冒頭で見積もった 30 N·m·s と概ね一致していることも確認できます。
アンローディング制御の効果
次に、磁気トルカ／スラスタによるアンローディングを組み込み、ホイールが飽和しないことを確認します。今回は GEO を想定し「スラスタによる定期アンロード」を簡易モデル化します。具体的には、毎月1回 $\bm{h}$ を $\bm{h}_{\mathrm{target}} = \bm{0}$ に近づける離散的なリセットを入れます。
import numpy as np
import matplotlib.pyplot as plt
surfaces, r_cm = make_geo_satellite(cm_offset=np.array([0.08, 0.03, 0.0]))
dt = DAY * 0.5
T_end = 15 * YEAR
ts = np.arange(0, T_end, dt)
h = np.zeros(3)
h_hist = np.zeros((len(ts), 3))
unload_events = [] # アンロード実施時刻のログ
UNLOAD_INTERVAL = 30 * DAY # 30日に1回アンロード
UNLOAD_EFFICIENCY = 0.95 # 1回で 95% のモーメンタムを降ろす
next_unload = UNLOAD_INTERVAL
delta_v_total = 0.0 # アンロードに使った推進剤の効果（ΔV相当）
ISP, M_SC, R_ARM = 220.0, 3000.0, 1.0
G0 = 9.80665
for i, t in enumerate(ts):
s_hat = sun_unit_vector(t)
_, tau = srp_force_torque(s_hat, surfaces, r_cm)
h = h + tau * dt
# アンロード
if t >= next_unload:
# ΔH を吐き出すのに必要な推進剤量を簡易積算
dH = np.linalg.norm(h) * UNLOAD_EFFICIENCY
impulse = dH / R_ARM # 推力×時間 [N·s]
delta_v_total += impulse / M_SC
h = h * (1 - UNLOAD_EFFICIENCY)
unload_events.append(t)
next_unload += UNLOAD_INTERVAL
h_hist[i] = h
mass_used = M_SC * (1 - np.exp(-delta_v_total / (ISP*G0)))
fig, ax = plt.subplots(figsize=(11, 4.5))
years = ts / YEAR
for k, lab in enumerate(['x','y','z']):
ax.plot(years, h_hist[:,k], label=fr'$h_{lab}$', lw=0.9)
ax.axhline(+50, color='k', ls='--', alpha=0.4)
ax.axhline(-50, color='k', ls='--', alpha=0.4)
ax.set_xlabel('Time [years]'); ax.set_ylabel('Wheel momentum [N·m·s]')
ax.set_title('Reaction wheel momentum WITH monthly thruster unloading')
ax.grid(alpha=0.3); ax.legend()
plt.tight_layout()
plt.savefig('momentum_15yr_with_unload.png', dpi=150, bbox_inches='tight')
plt.show()
print(f"15年間でのアンロード回数: {len(unload_events)} 回")
print(f"累積 ΔV（簡易見積もり）: {delta_v_total:.3f} m/s")
print(f"消費推進薬: {mass_used:.2f} kg")
このプロットでは、各軸のモーメンタムが月単位の小さな鋸歯状（sawtooth）パターンで蓄積・放出を繰り返しているのが見えます。最大値は概ね $\pm 5\,\mathrm{N\cdot m\cdot s}$ 程度に抑えられ、リアクションホイールの飽和ライン $\pm 50\,\mathrm{N\cdot m\cdot s}$ には到達しません。15年間で約 180 回のアンロード機動が必要となり、累積ΔV は数 m/s、推進薬消費は数 kg オーダーになります。これがアンローディングのコストであり、衛星の燃料予算で必ず確保しておかなければならない量です。
反射係数劣化の長期影響
最後に、太陽電池パドルの反射係数が経年劣化する効果を組み込み、SRP トルクの長期変動を見ます。これは「BOL（Begin Of Life）と EOL（End Of Life）でアンロード頻度をどれだけ変える必要があるか」という設計判断に直結します。
import numpy as np
import matplotlib.pyplot as plt
def make_aged_satellite(t, t_total=15*YEAR):
"""時刻 t [s] における劣化を反映した衛星モデル"""
aging = t / t_total # 0 → 1
# MLI と銀テフロンの典型的劣化: ρ_s が指数的に減少
rho_s_bus = 0.30 * np.exp(-1.5 * aging)
rho_s_pad = 0.08 * np.exp(-1.0 * aging)
rho_d_bus = 0.20 + 0.10 * aging # 拡散が増える
rho_d_pad = 0.15 + 0.05 * aging
# CoM が燃料消費で移動
cm_offset = np.array([0.05 + 0.05*aging,
0.02 + 0.02*aging,
0.0])
surfaces = []
bus_size = 2.5
a_bus = bus_size**2
for axis, sign in [(0,+1),(0,-1),(1,+1),(1,-1),(2,+1),(2,-1)]:
n = np.zeros(3); n[axis] = sign
r = (bus_size/2.0)*n
surfaces.append(dict(area=a_bus, n_hat=n, r=r,
rho_s=rho_s_bus, rho_d=rho_d_bus))
for sgn in [+1, -1]:
surfaces.append(dict(area=40.0,
n_hat=np.array([0,+1,0]),
r=np.array([sgn*(bus_size/2+5), 0, 0]),
rho_s=rho_s_pad, rho_d=rho_d_pad))
return surfaces, cm_offset
# 各年で 30 日平均トルクを評価
year_marks = np.arange(0, 15.5, 0.5)
mean_torque_mag = []
for yr in year_marks:
t0 = yr * YEAR
surfaces, r_cm = make_aged_satellite(t0)
samples = np.linspace(t0, t0 + 30*DAY, 60)
taus = []
for ts_ in samples:
s = sun_unit_vector(ts_)
_, tau = srp_force_torque(s, surfaces, r_cm)
taus.append(np.linalg.norm(tau))
mean_torque_mag.append(np.mean(taus))
plt.figure(figsize=(10, 4.5))
plt.plot(year_marks, np.array(mean_torque_mag)*1e6, 'o-')
plt.xlabel('Mission age [years]')
plt.ylabel(r'Mean |$\tau_{SRP}$| over 30 days [μN·m]')
plt.title('Aging effect: SRP torque magnitude vs mission age')
plt.grid(alpha=0.3)
plt.tight_layout()
plt.savefig('srp_torque_aging.png', dpi=150, bbox_inches='tight')
plt.show()
print(f"BOL 平均トルク = {mean_torque_mag[0]*1e6:.2f} μN·m")
print(f"EOL 平均トルク = {mean_torque_mag[-1]*1e6:.2f} μN·m")
print(f"EOL/BOL 比 = {mean_torque_mag[-1]/mean_torque_mag[0]:.2f}")
このグラフから、ミッション後期（EOL 付近）の SRP トルクが初期（BOL）に比べ 約1.5〜2倍に増えることが分かります。原因は、反射率の低下で力の方向が変わることと、燃料消費による CoM 移動でレバーアームが伸びることの組み合わせです。これは、BOL と同じ運用方針では 15 年目にはホイール飽和まで耐えられない可能性を示唆しており、運用側では「年が進むごとにアンロード周期を短くする」アダプティブ運用が必要になります。実機ではこの「年に応じたアンロード計画」を、衛星オペレータが運用計画に組み込みます。
姿勢誤差予算の積み上げ
ここまでの結果を統合して、15 年運用での姿勢誤差予算を作ります。SRP トルクが姿勢制御サイクル内でホイールへの指令にどう跳ね返り、最終的に何度の指向誤差を生むかを概算します。
import numpy as np
I_sc = np.array([2000., 2500., 1800.]) # 慣性モーメント [kg·m^2]
tau_srp_max = 12e-6 # 最大SRPトルク [N·m]
T_cycle = 100.0 # 姿勢制御サイクル [s]
omega_bw = 0.01 # 制御帯域 [Hz]
# 高周波外乱に対する姿勢誤差: θ_err ~ τ / (I·(2πf_bw)^2)
theta_hi = tau_srp_max / (np.min(I_sc) * (2*np.pi*omega_bw)**2)
# 低周波 DC 外乱に対する誤差: PI制御で残留誤差は積分項で抑えるが
# 帯域内では θ_err ~ τ / (I·omega_bw·K_i) 程度
theta_low = tau_srp_max / (np.min(I_sc) * (2*np.pi*omega_bw)**2) * 0.1
# パドル熱変形による光学誤差: 推定 0.01 deg
theta_thermal = np.deg2rad(0.01)
# トータル誤差予算（RSS）
theta_total = np.sqrt(theta_hi**2 + theta_low**2 + theta_thermal**2)
print(f"高周波SRP誤差 = {np.rad2deg(theta_hi)*3600:.1f} arcsec")
print(f"低周波SRP誤差 = {np.rad2deg(theta_low)*3600:.1f} arcsec")
print(f"パドル熱変形誤差 = {np.rad2deg(theta_thermal)*3600:.1f} arcsec")
print(f"合計姿勢誤差(RSS) = {np.rad2deg(theta_total)*3600:.1f} arcsec")
出力から、SRP 外乱由来の指向誤差はおおむね 十数〜数十秒角（arcsec） のオーダーで、放送通信衛星のビーム指向要求（典型的に $0.05^\circ \simeq 180\,\mathrm{arcsec}$）を満たすのに対し、高精度地球観測や天文ミッションで要求される $1\,\mathrm{arcsec}$ クラスの指向にはそれ単独で予算超過となることが分かります。地球観測高精度衛星では、SRP モデルの精緻化（CoP の事前同定、表面光学特性の地上試験データ反映）と高帯域制御の組み合わせが必須となる理由がここにあります。
まとめ
本記事では、GEO 通信衛星の姿勢制御に「見えない敵」として立ちはだかる太陽光圧トルクを、物理から長期運用まで一貫して見ました。
- 物理の正体: 太陽光は地球軌道上で $P_{\mathrm{srp}} = S_0/c = 4.56\,\mu\mathrm{N/m^2}$ の圧力を持つ。これが衛星表面に作用し、反射係数 $\rho_s, \rho_d$ に応じた力 $d\bm{F} = -P_{\mathrm{srp}}\cos\theta[(1-\rho_s)\hat{\bm{s}} + 2(\rho_s\cos\theta + \rho_d/3)\hat{\bm{n}}]dA$ を生む。
- トルクの正体: 重心 $\bm{r}_{cm}$ と光圧中心 $\bm{r}_{cp}$ のオフセットから $\bm{\tau} = (\bm{r}_{cp}-\bm{r}_{cm}) \times \bm{F}_{\mathrm{srp}}$。実機では避けがたく数 cm〜数十 cm のオフセットが生じ、結果として $1\text{–}20\,\mu\mathrm{N\cdot m}$ のトルクが定常的に作用する。
- 累積モーメンタム: 軌道平均で残る DC 成分が時間積分され、年間 $\sim 30\,\mathrm{N\cdot m\cdot s}$ のモーメンタムをリアクションホイールに蓄積する。15年運用ではアンロードが無ければホイールが飽和。
- アンローディング: LEO は磁気トルカ（$\bm{\tau} = \bm{m}\times\bm{B}$）、GEO はスラスタ（$\bm{\tau} = \bm{r}_{\mathrm{thr}}\times\bm{F}_{\mathrm{thr}}$）が主役。GEO では NSSK と同時に行うのが標準。
- 軌道への二次効果: SRP の並進力は離心率を年周で振動させ、振幅 $\Delta e \sim 10^{-4}$、東西振れ $\pm 40\,\mathrm{km}$。これも東西位置維持機動の理由のひとつ。
- 長期劣化: 反射係数の経年劣化と CoM 移動により、EOL の SRP トルクは BOL の約1.5〜2倍。運用周期も年とともに短縮する必要がある。
数値シミュレーションでは、15 年運用におけるホイールモーメンタム累積、月1回スラスタアンロードの効果、反射係数劣化の長期影響、最終姿勢誤差予算（数十 arcsec）を定量化しました。SRP は「ちりも積もれば」の典型例であり、設計時点で誤差予算と燃料予算の両方を見積もる必要があります。
次のステップとして、以下の記事も参考にしてください。