重力科学とドップラー計測 — 惑星重力場を電波で読み解く

この回で学ぶこと

トランスポンダの回で、探査機と地上局のコヒーレントターンアラウンドがμm/sオーダーという驚くべき精度でドップラー速度を測定できることを見ました。そしてそれまでの回では、この精密な速度計測を専ら航法——すなわち「探査機が今どこにいて、どの速度で飛んでいるか」を知るための道具として扱ってきました。前回は、レンジ・レンジレートの観測量から探査機の軌道状態を推定する観測方程式の枠組み(線形化・最小二乗フィッティング)を学びました。

この回では、視点を180度転換します。同じコヒーレントドップラーという道具を使って、今度は「探査機」ではなく「探査機が飛んでいる相手の天体」——木星や土星のような惑星——の重力場そのものを知るために使うのです。これを重力科学 (gravity science) と呼びます。

一見、探査機の軌道決定と惑星の重力場推定は別の問題に見えるかもしれません。しかし実際には表裏一体です。探査機の軌道は重力場によって決まるのですから、逆に「探査機の軌道(のごく僅かなズレ)を精密に測れば、そこから重力場を復元できる」はずです。この回では、惑星重力ポテンシャルを球面調和関数で展開する方法、その不均一性が探査機の速度に及ぼす摂動、そしてドップラー残差から重力場係数を最小二乗法で逆算する枠組みを数式で追います。最後に、ジュノー探査機による木星重力場マッピングと、カッシーニ探査機を使った一般相対性理論の検証実験という、この技術がもたらした2つの代表的な科学的成果を見ていきます。

直感的な全体像

もし惑星が完全に均質な球体(質量が中心対称に分布している)だったなら、その重力場は理想的な点質量と同じで、探査機はどこを通ってもきれいなケプラー軌道(あるいは双曲線軌道)を描きます。このときのドップラーシフトの時間変化は、二体問題の運動方程式から寸分違わず予測でき、地上局が測るドップラー曲線は理論曲線とぴったり重なるはずです。

しかし実際の惑星は完全な球体ではありません。自転による偏平化(赤道方向に膨らむ)、内部の質量分布のムラ(マントルの対流、コアの偏心、木星なら深い大気循環に伴う質量の偏り)など、様々な理由で重力場は「理想的な点質量からのズレ」を持っています。探査機がそのような不均一な重力場の中を通過すると、場所によって受ける重力の強さや向きがわずかに変わり、軌道速度が理想的な二体問題からわずかにズレます。

このズレは、地上局が観測するドップラーシフトの中に理論値からの微小な揺らぎとして現れます。これをドップラー残差 (Doppler residual)——観測値から、単純な二体問題モデルによる予測値を差し引いたもの——と呼びます。ドップラー残差の時間変化のパターンを詳しく解析すれば、逆に「探査機をこのように揺らした重力場は、どんな形をしていたはずか」を推定できます。これはちょうど、水面に浮かべたコルクの揺れ方を精密に観測することで、水面下に隠れた岩の形を推定するようなものです。探査機は惑星の重力場という「見えない地形」を、自分の軌道の揺れを通じて我々に伝えてくれる、極めて精密なプローブなのです。

数式による定式化

重力ポテンシャルの球面調和展開

まず、不均一な重力場をどう数式で表現するかから始めます。惑星の重力ポテンシャル $U(r,\theta,\phi)$($r$: 惑星中心からの距離、$\theta$: 余緯度、$\phi$: 経度)は、惑星の質量分布がどれほど複雑であっても、球面調和関数 (spherical harmonics) という基底関数の重ね合わせとして厳密に展開できます。回転対称性が高い惑星(木星・土星など、経度方向にほぼ一様)を考える場合、経度に依存しない**帯調和項(zonal harmonics)**だけを残した簡略形が有用です。

$$U(r,\theta) = \frac{GM}{r}\left[1 - \sum_{n=2}^{\infty} J_n \left(\frac{R}{r}\right)^n P_n(\cos\theta)\right]$$

ここで、$GM$ は惑星の質量パラメータ、$R$ は基準半径(通常は赤道半径)、$P_n(\cos\theta)$ は $n$ 次のルジャンドル多項式、そして $J_n$ が帯調和係数と呼ばれる無次元の重力場係数です。$n=0$ の項(単なる $GM/r$、点質量ポテンシャル)は展開の外に置かれ、$J_n$ たちが「点質量からのズレ」の大きさを表します。

最も支配的なのは $n=2$ の項、係数 $J_2$ です。これは天体の扁平率(自転による赤道方向の膨らみ)に対応する項で、$P_2(\cos\theta) = \tfrac{1}{2}(3\cos^2\theta - 1)$ という形をしています。木星・土星のように自転が速いガス惑星では、$J_2$ は他の高次項に比べて桁違いに大きく、重力場の非一様性のほとんどを支配します。より高次の $J_4, J_6, \ldots$(偶数次のみが赤道対称な扁平構造に対応)は、内部質量分布のより細かい形状——木星の場合は深い大気循環に伴う質量の再配分——を反映します。実際には経度方向の非対称性を捉えるセクトラル項・テッセラル項($C_{nm}, S_{nm}$、$m\neq0$)も一般には必要ですが、この回では主に帯調和項による定式化に絞ります。

摂動加速度とドップラー残差への伝搬

探査機に働く重力加速度は、ポテンシャルの勾配 $\vec{a} = \nabla U$ で与えられます。これを二体問題の中心力項($n=0$)と、それ以外の摂動項に分けて書くと、

$$\vec{a}(\vec r) = -\frac{GM}{r^2}\hat r + \vec{a}_{pert}(\vec r; J_2, J_4, \ldots)$$

となります。摂動加速度 $\vec{a}_{pert}$ の大きさは、$J_2$ 項だけを取り出すと、赤道面通過時でおおむね

$$a_{J_2} \sim \frac{3}{2}\, J_2 \left(\frac{R}{r}\right)^2 \frac{GM}{r^2}$$

という程度のオーダーになります(厳密な方向・大きさは緯度 $\theta$ に依存しますが、ここではスケール感を掴むための近似式として扱います)。この摂動加速度が探査機の軌道速度にもたらす変化 $\delta \vec v(t)$ は、運動方程式 $\dot{\vec v} = \vec a$ を時間積分すれば得られ、フライバイのように接近通過が短時間で終わる場合、通過時間 $\tau$ のオーダーで

$$\delta v \sim a_{J_2}\, \tau$$

という大まかな見積もりができます(厳密には数値積分、あるいは摂動論の変分方程式で求めます)。

そして最後の1ステップが、トランスポンダの回で導いたコヒーレント2-wayドップラーの式そのものです。視線方向速度の摂動 $\delta v_r$(視線方向成分)は、ダウンリンク周波数 $f_c$ に対して

$$\delta f = \frac{2\,\delta v_r}{c}\, f_c$$

というドップラー周波数の揺らぎとして現れます。まとめると、重力場の非一様性($J_n$)→ 摂動加速度 → 速度の揺らぎ → ドップラー残差、という一直線の因果の鎖が数式でつながっていることが分かります。逆にこの鎖を遡れば、観測されたドップラー残差から重力場係数 $J_n$ を推定できるはずだ、というのがこの後の逆問題の考え方です。

ドップラー残差の定義と観測方程式

実務では、まず現時点で最良と考えられる重力場モデル(暫定的な $J_2^{(0)}, J_4^{(0)}, \ldots$ の値と軌道初期状態)を使って、探査機の軌道と、そこから予測されるドップラーシフトの時系列 $f_{model}(t)$ を計算します(これを参照軌道 (reference trajectory) と呼びます)。実際に地上局で観測されたドップラー周波数 $f_{obs}(t)$ との差が、ドップラー残差です。

$$\Delta f(t) \equiv f_{obs}(t) - f_{model}(t)$$

もし参照軌道の重力場モデルが完全に正しければ、$\Delta f(t)$ は観測雑音だけを反映したランダムな揺らぎになるはずです。しかし実際の重力場係数 $J_n$ の真値が、モデルで使った暫定値 $J_n^{(0)}$ からわずかにズレていれば、$\Delta f(t)$ には雑音に埋もれない系統的なパターンが残ります。これが重力場を「読み取る」手がかりです。

前回、レンジ・レンジレートの観測方程式を線形化・最小二乗フィッティングする枠組みを学びました。ここでも全く同じ枠組みが使えます。推定したい未知パラメータをベクトル $\vec x = (J_2, J_4, \ldots, \vec r_0, \vec v_0, \ldots)^T$(重力場係数に加えて、軌道初期状態そのものも通常同時推定します)とし、暫定値 $\vec x_0$ のまわりで観測量を線形化すると、

$$\Delta f_i \approx \sum_j \frac{\partial f_{model}}{\partial x_j}\bigg|_{\vec x_0} \Delta x_j + \varepsilon_i \equiv \sum_j H_{ij}\, \Delta x_j + \varepsilon_i$$

という観測方程式が得られます。ここで感度行列(design matrix)の要素 $H_{ij} = \partial f_{model}/\partial x_j$ は「重力場係数 $x_j$ が単位量だけ変化したら、時刻 $t_i$ のドップラーがどれだけ変化するか」という偏微分で、実務では軌道の運動方程式に付随する**変分方程式(状態遷移行列)**を数値積分することで計算されます。

観測を多数の時刻 $i=1,\ldots,M$ にわたって集めれば、行列形式で

$$\Delta \vec f = H\, \Delta \vec x + \vec\varepsilon$$

と書け、前回導いた重み付き最小二乗解と全く同じ形の正規方程式(normal equation)

$$\widehat{\Delta \vec x} = \left(H^T W H\right)^{-1} H^T W\, \Delta \vec f$$

を解くことで、重力場係数の補正量 $\widehat{\Delta \vec x}$、したがって改良された推定値 $\vec x_0 + \widehat{\Delta \vec x}$ が得られます。重み行列 $W$ には、個々のドップラー観測の精度(トランスポンダの回で導出した速度推定誤差 $\sigma_v$、あるいはそれに対応する周波数誤差 $\sigma_f$ の逆分散)がそのまま入ります。つまり、PLLのループSNRから始まった精度の連鎖が、最終的に「木星の重力場係数がどこまで正確に決まるか」というサイエンスの精度に直結しているのです。

数値でスケール感をつかむ

木星のジュノー探査機を例に、大まかな数値でこの連鎖を追ってみましょう。木星の $GM \approx 1.267\times10^{17}\ \text{m}^3/\text{s}^2$、赤道半径 $R \approx 7.149\times10^{7}\ \text{m}$、$J_2 \approx 1.470\times10^{-2}$ です。ジュノーの近木点(perijove)通過時、木星中心からの距離をおよそ $r \approx 1.06\,R \approx 7.58\times10^{7}\ \text{m}$ とすると、

$$\frac{GM}{r^2} \approx \frac{1.267\times10^{17}}{(7.58\times10^{7})^2} \approx 22.1\ \text{m/s}^2$$

が中心重力の加速度、そして $J_2$ 項による摂動加速度は

$$a_{J_2} \sim \frac{3}{2}\times(1.470\times10^{-2})\times\left(\frac{7.149\times10^7}{7.58\times10^7}\right)^2\times 22.1 \approx 0.43\ \text{m/s}^2$$

となり、中心重力の約2%に相当します。この摂動が仮に $\tau\approx60\ \text{s}$ 程度の時間作用し続けたとすると、$\delta v \sim 0.43\times60 \approx 26\ \text{m/s}$ もの速度変化がもたらされ、Xバンド($f_c=8.4\ \text{GHz}$)のドップラーシフトに換算すると

$$\delta f = \frac{2\times26}{3\times10^8}\times8.4\times10^9 \approx 1.5\times10^{3}\ \text{Hz}$$

というkHzオーダーの大きな信号になります。この規模の信号は、$J_2$ の値をあらかじめ暫定的にでも良く見積もっておかなければ、そもそも探査機を見失いかねないほど大きいことが分かります——だからこそ、$J_2$ は最初から参照軌道モデルに組み込んでおく必要があるのです。

サイエンスとして意味を持つのは、この暫定 $J_2^{(0)}$ からのさらに小さな補正量 $\Delta J_2$ です。感度 $\partial(\delta f)/\partial J_2 \sim \delta f / J_2 \approx 1.5\times10^3/1.47\times10^{-2} \approx 1.0\times10^{5}\ \text{Hz}$(単位 $J_2$ あたり)という比例関係を使うと、トランスポンダの回で導いたドップラー残差の測定精度 $\sigma_f \sim 10^{-4}\text{–}10^{-2}\ \text{Hz}$ 程度(積分時間・リンク条件による)から逆算すると、

$$\Delta J_2 \sim \frac{\sigma_f}{1.0\times10^{5}} \sim 10^{-9}\text{–}10^{-7}$$

という桁の分解能が単一パスからでも見込めることになります。実際にジュノーは30回を超える近木点フライバイを重ね、統計的に組み合わせることで、木星の $J_2$ を $J_2 = (14696.572\pm0.014)\times10^{-6}$(Iess et al., 2018, Nature)という、相対精度で $10^{-9}$ を切る驚異的な精度まで決定しています。

一般相対性理論の検証への応用: シャピロ遅延とドップラー

重力場の不均一性だけでなく、重力そのものが電波の伝搬に与える影響——一般相対性理論的効果——もコヒーレントドップラーで検出できます。太陽の重力場のように強い(といっても弱い場近似の範囲内ですが)重力ポテンシャルの近くを電波が通過すると、時空の歪みによって光路の伝搬時間がわずかに余分にかかります。これをシャピロ遅延 (Shapiro delay) と呼び、太陽近傍を通る地球〜探査機間の電波の余分な遅延は、近似的に

$$\Delta t_{Shapiro} \approx (1+\gamma)\,\frac{2GM_\odot}{c^3}\ln\!\left(\frac{r_E + r_S + d}{r_E + r_S - d}\right)$$

と表せます。ここで $r_E, r_S$ は太陽から見た地球・探査機の距離、$d$ は地球〜探査機間の距離、そして $\gamma$ がPPN(Parametrized Post-Newtonian)パラメータと呼ばれる、一般相対性理論を検証するための無次元パラメータです。一般相対性理論が正しければ $\gamma=1$ と予言され、$\gamma$ が1からズレていれば、それは一般相対性理論からの逸脱を意味します。

探査機が太陽の向こう側を通過する**太陽合(solar conjunction)**の前後では、電波の視線が太陽近傍を掠める幾何になり、$\Delta t_{Shapiro}(t)$ が時間とともに急激に変化します。この時間変化率がそのままコヒーレント2-wayドップラーの残差として現れます——通常のドップラー残差が「速度」由来だったのに対し、こちらは「時空の歪みによる余分な光路長の変化率」由来だという点が異なりますが、地上局が測るのは同じ $\Delta f$ という物理量です。太陽への最接近前後で符号が反転する鋭い山型のドップラー残差シグネチャが、シャピロ遅延の指紋です。

カッシーニ探査機は2002年の太陽合の際、X/Kaバンドの複数リンクを同時にコヒーレントに運用してこの残差を精密に測定し、

$$\gamma - 1 = (2.1\pm2.3)\times10^{-5}$$

という、当時としては一般相対性理論の最も精密な検証結果を得ました(Bertotti, Iess, Tortora, 2003, Nature)。この精度は、まさにこの回で見てきた重力場係数の推定と全く同じ最小二乗法の枠組み(観測方程式の未知パラメータに $\gamma$ を加えて解く)によって達成されています。

実務での使われ方

• ジュノー探査機(NASA/JPL、2016年〜木星周回)は、Ka-band Translator (KaT) と呼ばれる超高精度のKaバンドコヒーレントリンクを搭載し、近木点フライバイのたびにドップラー残差を取得して木星の重力場係数を推定しています。特に注目されたのは、$J_2, J_4$ に加えて奇数次の帯調和項($J_3, J_5, J_7,\ldots$、通常は南北対称な惑星ではゼロに近いはず)が有意にゼロでないことが発見された点で、これは木星の深部大気循環(赤道ジェット)が南北非対称であり、少なくとも深さ3000km程度まで及んでいることを示す証拠となりました(Iess et al., 2018, Nature)。
• **カッシーニ探査機(NASA/ESA/ASI、2004–2017年土星周回)**は、土星本体および衛星(特にエンケラドゥスやタイタン)の重力場マッピングにコヒーレントドップラーを使用しました。エンケラドゥスの複数回の近接フライバイによるドップラー計測からは、氷の地殻の下に液体の海が存在することを示す重力異常が検出され、後の間欠泉(プルーム)発見と合わせて、この衛星の海洋の存在を確定づける証拠の一つとなりました。
• **GRAIL(NASA、2011–2012年月周回)**は、2機の探査機が互いにKaバンドで衛星間ドップラー測距(Ka-Band Ranging System)を行うことで、月の重力場を従来の地球からの片方向トラッキングでは不可能だった空間分解能でマッピングしました。この手法は地球のGRACEミッションで培われた衛星間ドップラー技術の応用であり、月の裏側のクレーター構造や地殻の厚さの推定に貢献しています。
• 一般相対性理論の検証は、カッシーニのシャピロ遅延実験(PPNパラメータ $\gamma$)以外にも、太陽系内の惑星探査機のドップラー・レンジデータを用いた重力定数の時間変化の上限設定や、強い等価原理の検証など、複数の基礎物理実験に応用されています。これらはすべて、DSN局の高安定な周波数基準(水素メーザー)とコヒーレントトランスポンダが可能にする、ppm(百万分率)を遥かに超えるppb〜ppt級の精密計測に支えられています。

演習問題

1. 土星の $GM \approx 3.793\times10^{16}\ \text{m}^3/\text{s}^2$、赤道半径 $R\approx 6.027\times10^{7}\ \text{m}$、$J_2\approx1.629\times10^{-2}$ とします。カッシーニがある近接フライバイで土星中心から $r\approx1.1R$ の距離を通過したとして、本文の近似式 $a_{J_2}\sim\frac{3}{2}J_2(R/r)^2(GM/r^2)$ を使い、$J_2$ 項による摂動加速度の大きさを見積もってください。中心重力加速度 $GM/r^2$ に対する比(%)も求めてください。

2. 観測方程式 $\Delta f_i \approx H_i\,\Delta J_2 + \varepsilon_i$(未知パラメータが $J_2$ 一つだけの単純化されたケース)を考えます。3時刻の観測で感度が $H_1=8.0\times10^4,\ H_2=9.5\times10^4,\ H_3=7.2\times10^4$ Hz(単位 $J_2$ あたり)、対応するドップラー残差が $\Delta f_1=0.05,\ \Delta f_2=0.09,\ \Delta f_3=0.03$ Hz、すべての観測の重みが等しい($W=I$)とします。スカラー版の正規方程式 $\widehat{\Delta J_2} = \dfrac{\sum_i H_i \Delta f_i}{\sum_i H_i^2}$ を使って $\widehat{\Delta J_2}$ を求めてください。

3. カッシーニのシャピロ遅延実験では $\gamma-1=(2.1\pm2.3)\times10^{-5}$ という結果が得られました。この不確かさのオーダーが、トランスポンダの回で導いたドップラー速度計測精度(μm/s〜mm/sオーダー)や、太陽合の前後で電波が太陽近傍を通過する期間の長さと、どのように関係していると考えられるか、自分の言葉で説明してください(定性的な考察で構いません)。

4. 探査機の軌道決定(航法)と、この回で扱った重力場推定(重力科学)は、どちらも「観測方程式を線形化し、最小二乗法で未知パラメータを推定する」という全く同じ数学的枠組みを使います。両者の違いは、未知パラメータベクトル $\vec x$ に何を含めるかという一点に集約されると言えます。この観点から、なぜジュノーは1回のフライバイだけでなく30回以上のフライバイを重ねてようやく高精度な木星重力場を決定できたのか、観測方程式の行数(観測数 $M$)と未知パラメータ数の関係を踏まえて考察してください。

まとめと次回予告

コヒーレントドップラーは、探査機自身の軌道を知るための「航法」の道具であると同時に、視点を転換すれば、探査機が通過する天体の重力場そのものを読み取る「重力科学」の道具にもなります。重力ポテンシャルを球面調和関数($J_2, J_4,\ldots$)で展開し、その不均一性が引き起こす軌道摂動がドップラー残差として観測されること、そしてその残差を最小二乗法で逆問題として解くことで重力場係数を推定する枠組みを見ました。ジュノーによる木星深部大気循環の発見や、カッシーニによる一般相対性理論の精密検証は、いずれもこの枠組みの上に成り立つ成果です。

ここまで、地上局・探査機双方の周波数基準がどれほど正確かが、最終的な速度計測精度、ひいては重力場係数の決定精度に直結することを繰り返し見てきました。次回は、この精度の源泉である**原子時計とUSO(超高安定発振器)**そのものに焦点を当て、水素メーザーやルビジウム発振器がどのようにして $10^{-15}$ 台という驚異的な周波数安定度を実現しているのか、その物理と工学に踏み込んでいきます。

参考文献

• L. Iess et al., “Measurement of Jupiter’s asymmetric gravity field,” Nature 555, 220–222 (2018)
• B. Bertotti, L. Iess, P. Tortora, “A test of general relativity using radio links with the Cassini spacecraft,” Nature 425, 374–376 (2003)
• W. M. Kaula, Theory of Satellite Geodesy, Blaisdell (1966; Dover復刊)
• B. D. Tapley, B. E. Schutz, G. H. Born, Statistical Orbit Determination, Elsevier Academic Press
• M. T. Zuber et al., “Gravity Field of the Moon from the Gravity Recovery and Interior Laboratory (GRAIL) Mission,” Science 339, 668–671 (2013)
• C. M. Will, Theory and Experiment in Gravitational Physics, 2nd ed., Cambridge University Press
• DSN Telecommunications Link Design Handbook, DSN No. 810-005