探査機搭載時計の長期ドリフトモデリング — 二次多項式モデルと較正運用

この回で学ぶこと

原子時計とUSOの回では、発振器の周波数安定度をアラン分散で定量化しました。あの回で扱ったのは、積分時間 $\tau$ が秒〜時間オーダーの短期的な統計的揺らぎ——ホワイト周波数雑音、フリッカー雑音、ランダムウォーク——であり、いわば「発振器がどれだけ細かく震えているか」の話でした。また双方向時刻転送の回では、離れた2つの時計の瞬間のオフセット $\Delta T$ を測る方法を学びましたが、そこでは測定の数秒〜数分のあいだ「オフセットは一定」とみなせることを前提にしていました。

しかし探査機のミッションは数ヶ月、数年、ときに数十年続きます。このタイムスケールで探査機搭載のUSO(Ultra-Stable Oscillator)を眺めると、まったく別の顔が見えてきます。周波数が統計的に揺らぐのではなく、一方向にじわじわとずれ続けるのです。打ち上げ直後には公称値ぴったりだった発振周波数が、1年後には $10^{-9}$ のオーダーで系統的にシフトしている——これは雑音ではなく、水晶の経年変化(エージング)や放射線損傷の蓄積という物理過程が生む決定論的なドリフトです。

この回では、この長期ドリフトを扱うための標準的な道具である二次多項式時計モデル

$$x(t) = x_0 + y_0\, t + \frac{1}{2} D t^2 + \epsilon(t)$$

を導入し、(1) なぜこの形なのか、(2) 各パラメータを実測データからどう推定するのか、(3) 推定を怠ると軌道決定にどんな系統誤差が混入するのか、(4) 実際の深宇宙ミッションはどんな較正運用を回しているのか、を順に見ていきます。アラン分散が「揺らぎの統計」だったのに対し、今回は「トレンドの推定」——同じ発振器を、統計学でいう定常過程の分析から回帰分析の対象へと視点を切り替える回です。

直感的な全体像 — 「震え」と「ずれ」は別物

腕時計にたとえてみましょう。安物のクォーツ時計でも、1秒のあいだの刻みの揺らぎは人間には知覚できないほど小さい——これが短期安定度、アラン分散の世界です。ところが同じ時計を1ヶ月放置すると、確実に数秒〜数十秒進むか遅れるかしています。しかも多くの場合、「毎月だいたい同じ方向に、だいたい同じだけ」ずれます。これが長期ドリフトです。

重要なのは、この2つがまったく異なる対処を要求することです。

• 短期の揺らぎは確率的で予測不可能。統計的に特徴付けて(アラン分散)、平均化や積分時間の設計で影響を抑えるしかない。
• 長期のドリフトは決定論的で予測可能。パラメータを推定してモデル化すれば、補正によってほぼ取り除ける。

アラン分散の回で「線形ドリフト成分は差分を取ることで打ち消される」と述べたのを思い出してください。あれはアラン分散が短期安定度の指標として機能するための工夫でしたが、裏を返せば、アラン分散はドリフトの情報を意図的に捨てているということです。ドリフトはアラン分散にとっては邪魔者でしたが、数ヶ月にわたる探査機運用にとっては主役級の誤差要因であり、独立にモデル化して推定してやる必要があります。

そして、ここで双方向時刻転送が効いてきます。TWSTTや往復測距による時刻比較は「ある瞬間のオフセット $\Delta T$」を与える測定でした。この測定を数週間〜数ヶ月にわたって繰り返し、得られたオフセットの時系列に多項式をフィットする——これが長期ドリフト推定の基本戦略です。1回の測定では「いまのずれ」しか分かりませんが、時系列として眺めれば「ずれの速さ(周波数オフセット)」と「ずれの加速(ドリフト率)」まで読み取れるのです。

数式による定式化

二次多項式時計モデル

探査機の時計が表示する時刻を、理想的な基準時刻 $t$ の関数として $C(t) = t + x(t)$ と書き、時刻誤差 $x(t)$(単位: 秒)を次のようにモデル化します。これが時刻・周波数の分野で**標準的な時計モデル(clock model)**と呼ばれる形です。

$$x(t) = x_0 + y_0 (t - t_0) + \frac{1}{2} D (t - t_0)^2 + \epsilon(t)$$

各項の意味を確認しましょう($t_0$ はモデルの基準エポックです)。

• $x_0$ [s]: 時刻オフセット。エポック $t_0$ における時計のずれそのもの。双方向時刻転送で測った $\Delta T$ に対応します。
• $y_0$ [無次元]: 周波数オフセット。エポックにおける相対周波数偏差 $y(t_0)$。周波数が $y_0$ だけずれていれば、時刻誤差は毎秒 $y_0$ 秒ずつ線形に積み上がります。
• $D$ [1/s]: 周波数ドリフト率。相対周波数偏差が単位時間あたりに変化する割合で、この回の主役です。周波数が線形に変化すれば、その積分である時刻誤差は二次関数的に成長します。
• $\epsilon(t)$: 上の決定論的成分で説明しきれない確率的残差。ランダムウォーク周波数雑音などが積分されて時刻に現れた成分で、こちらの統計的性質を特徴付けるのがアラン分散の役割でした。

相対周波数偏差 $y(t)$ は時刻誤差の微分なので、このモデルのもとでは

$$y(t) = \frac{dx}{dt} = y_0 + D (t - t_0) + \dot{\epsilon}(t)$$

となります。つまりこのモデルは「周波数が時間に対して線形にドリフトする」という仮定と等価です。より高次のドリフト(周波数の二次変化)を入れることも形式上は可能ですが、後述するようにエージングの物理は時間とともに減速する(対数則に従う)ため、数ヶ月程度の運用区間を区切って二次モデルを当てはめ、区間ごとにパラメータを更新するのが実務の定石です。

この3パラメータ+確率成分という構成は特定機関のローカルルールではなく、事実上の業界標準です。GPSの航法メッセージが放送する衛星時計補正係数 $(a_{f0}, a_{f1}, a_{f2})$ はまさにこの $(x_0, y_0, D/2)$ に対応しますし、DSNの軌道決定ソフトウェア(ODP、およびその後継のMONTE)やGNSSの精密解析ソフトも、探査機・衛星の時計をこの多項式+確率過程の形で状態変数として保持しています。

数値で感覚を掴む — ドリフトの二次的蓄積

ドリフト項のおそろしさは、時刻誤差が時間の二乗で効いてくる点にあります。探査機USOとして典型的な、1日あたり $1 \times 10^{-11}$ のドリフト率を考えましょう。

$$D = \frac{1 \times 10^{-11}}{86400\ \text{s}} \approx 1.16 \times 10^{-16}\ \text{s}^{-1}$$

一見、無視できそうなほど小さい数字です。ところが較正なしで30日間($t = 2.59 \times 10^6$ s)放置すると、蓄積する時刻誤差は

$$\frac{1}{2} D t^2 = \frac{1}{2} \times (1.16 \times 10^{-16}) \times (2.59 \times 10^6)^2 \approx 3.9 \times 10^{-4}\ \text{s}$$

つまり約0.39ミリ秒。One-way計測でこの時刻誤差がそのまま距離誤差に化けると、$c \times 3.9\times10^{-4} \approx 1.2 \times 10^5$ m、100 km超に相当します。周波数側で見ても、30日後には $y = D t \approx 3 \times 10^{-10}$ となり、USOの短期安定度($\sigma_y \sim 10^{-13}$)より3桁も大きい。長期運用において、モデル化されないドリフトが誤差収支の中でいかに支配的になるかが分かります。

ドリフトの物理的起源

$D$ という1個のパラメータの背後には、複数の物理過程が畳み込まれています。

(1) 水晶のエージング(経年変化)。 水晶発振器の物理設計の回で見たように、USOの心臓部は恒温槽制御された水晶振動子(OCXO)です。水晶の共振周波数は振動子の質量と実効弾性に依存するため、(a) 電極や振動子表面への分子の吸着・脱離による質量変化、(b) 製造時に導入された結晶格子の応力の緩和、(c) 実装マウントの応力緩和、といった微視的過程が周波数をゆっくり動かします。エージング速度は稼働直後がもっとも速く、時間とともに対数関数的に減速する($\Delta f/f \propto \ln(1 + t/\tau_a)$ という経験則)ことが知られており、高品質OCXOで1日あたり $10^{-11}$〜$10^{-10}$、十分にエージングが進んだ宇宙用USOで $10^{-12}$〜$10^{-11}$/日が典型的なオーダーです。打ち上げ前に地上で長期間通電して「枯らす」のは、この初期の速いエージング領域を通過させておくためです。

(2) 放射線損傷の蓄積。 深宇宙では探査機は継続的に銀河宇宙線・太陽粒子にさらされ、水晶内部に格子欠陥と電離損傷が蓄積します。累積線量(TID)に対する周波数感度は、掃引処理(sweeping)された高純度水晶で低線量域において1 radあたり $10^{-12}$〜$10^{-11}$ 程度のオーダーであり、長期的には単調な周波数シフトとして、また太陽フレア時にはステップ状の変化として現れます。木星圏など高放射線環境のミッションでは、遮蔽設計とあわせてこの効果が時計モデルの更新頻度を左右します。

(3) 温度サイクル応力。 恒温槽があっても、探査機の姿勢変更・電力モード変更・日陰通過などによる外界温度の変動はわずかに槽内へ漏れ込みます。温度変化そのものによる短期的な周波数変動に加え、繰り返しの熱サイクルが実装部の応力状態を変え、エージング曲線を「リセット」あるいは屈曲させることがあります。実測されたドリフト率が姿勢運用の変更後に不連続に変わる、という現象の一因です。

これらはいずれも「単調・緩慢・準決定論的」という共通の性格を持つため、まとめて二次多項式の $D$ に押し込み、区間ごとに再推定する、というのがモデリング上の割り切りです。

最小二乗フィットによるパラメータ推定

では $(x_0, y_0, D)$ を実測からどう決めるか。素材は双方向時刻転送(探査機の場合は往復測距とテレメトリのSCLKタイムスタンプの組み合わせ、いわゆるSCLK-SCET相関)で得られる時刻オフセットの時系列です。時刻 $t_1, t_2, \dots, t_N$ における測定値を $\tilde{x}_i = x(t_i) + v_i$($v_i$ は測定雑音、分散 $\sigma^2$)とすると、観測方程式は

$$\tilde{x}_i = x_0 + y_0 (t_i - t_0) + \frac{1}{2} D (t_i - t_0)^2 + v_i$$

これは未知パラメータベクトル $\boldsymbol{\theta} = (x_0, y_0, D)^\mathsf{T}$ に対して線形なので、通常の最小二乗法がそのまま使えます。$\Delta t_i \equiv t_i - t_0$ として計画行列を

$$H = \begin{pmatrix} 1 & \Delta t_1 & \tfrac{1}{2}\Delta t_1^2 \\ 1 & \Delta t_2 & \tfrac{1}{2}\Delta t_2^2 \\ \vdots & \vdots & \vdots \\ 1 & \Delta t_N & \tfrac{1}{2}\Delta t_N^2 \end{pmatrix}$$

と置けば、正規方程式の解と推定誤差共分散は

$$\hat{\boldsymbol{\theta}} = (H^\mathsf{T} H)^{-1} H^\mathsf{T} \tilde{\boldsymbol{x}}, \qquad P = \sigma^2 (H^\mathsf{T} H)^{-1}$$

で与えられます。実務手順としては、(1) 数日〜数週間おきに時刻比較セッションを実施してオフセット系列を蓄積し、(2) この二次フィットを走らせ、(3) 残差 $\tilde{x}_i - H\hat{\boldsymbol{\theta}}$ を検査してモデルの妥当性(残差が白色的か、系統的な曲がりが残っていないか)を確認する、という流れになります。

推定精度について重要な性質を1つ。測定間隔を一定、観測区間の長さを $T$ とすると、$H^\mathsf{T}H$ の各要素は $T$ のべきでスケールするため、ドリフト率の推定分散はおおよそ

$$\mathrm{Var}(\hat{D}) \sim \frac{\sigma^2}{N} \cdot \frac{c_D}{T^4}$$

($c_D$ は定数)と、観測区間の4乗に反比例して改善します。二次の係数は「曲がり」を見て決めるため、長い区間ほど曲がりがはっきり見え、劇的に精度が上がるのです。ただし区間を延ばしすぎると、今度は $\epsilon(t)$ に含まれるランダムウォーク周波数雑音(アラン分散のログログプロットで傾き $+1/2$ の領域)の蓄積が測定雑音を上回り、さらにエージングの対数則からの乖離(二次モデル自体の破綻)も効き始めます。フィット区間の長さは、測定雑音・確率的雑音・モデル誤差の三つ巴のトレードオフで決まる——ここでもアラン分散のU字カーブが、区間設計の指針として顔を出します。

軌道決定への系統誤差混入と時計パラメータの同時推定

未較正のドリフトを放置すると何が起きるか。One-wayのドップラー・レンジング観測量には、幾何学的な情報(距離・視線速度)と時計誤差が加算的に混ざります。たとえばOne-wayレンジ観測量は近似的に

$$\rho_{\text{obs}}(t) = \rho_{\text{true}}(t) + c \left[ x_0 + y_0 (t-t_0) + \frac{1}{2} D (t-t_0)^2 \right] + (\text{雑音})$$

という形になります。厄介なのは、右辺第2項の各成分が軌道パラメータの誤差シグネチャと部分的に相関することです。定数項 $c\,x_0$ は距離バイアスと、線形項 $c\,y_0 t$ は視線速度の誤差と、二次項は視線加速度——すなわち探査機に働く力(重力・太陽輻射圧)のモデル誤差——とそれぞれ見分けがつきにくい。軌道決定フィルタは観測残差を最小にしようとして、時計のずれを軌道側のパラメータに誤って吸収し、軌道解が系統的に歪みます。重力場推定を行う科学ミッションでは、この汚染が重力場係数の誤差として科学成果に直結します。

標準的な対策は、カルマンフィルタによる軌道決定で学んだ状態空間の枠組みに、時計パラメータを状態ベクトルごと組み込んで同時推定することです。軌道の状態 $(\boldsymbol{r}, \boldsymbol{v})$ に時計状態 $(x, y, D)$ を連接し、時計部分のダイナミクスを

$$\frac{d}{dt}\begin{pmatrix} x \\ y \\ D \end{pmatrix} = \begin{pmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} x \\ y \\ D \end{pmatrix} + \boldsymbol{w}(t)$$

という線形系(いわゆる3状態時計モデル)で与えます。プロセス雑音 $\boldsymbol{w}(t)$ の強度は、その発振器のアラン分散から逆算して設定します——ホワイト周波数雑音は $y$ への、ランダムウォーク周波数雑音は $D$ 方向への雑音強度に対応し、短期統計(アラン分散)と長期モデル(多項式)がここで1つのフィルタに統合されるわけです。GPSの精密軌道・時計決定でも、DSNの深宇宙軌道決定でも、この「軌道と時計の同時推定」は標準装備であり、時計状態に適切なプロセス雑音を与えることで、幾何と時計の相関を観測の幾何学的多様性(複数局・長期弧・DDORなど独立観測型の混合)によって解きほぐしています。

実務での使われ方

• GPS航法メッセージの時計補正係数: 各GPS衛星は自身の時計モデル $(a_{f0}, a_{f1}, a_{f2})$——時刻オフセット・周波数オフセット・ドリフト率(の1/2)——を航法メッセージで放送しており、受信機は擬似距離計算の前にこの二次多項式で衛星時計誤差を補正します。係数は管制セグメント(Colorado Springsの主管制局)が全世界の監視局データへの推定で決定し、定期的にアップロードします。ルビジウム原子時計を積む現行世代のGPS衛星のドリフト率は $10^{-14}$〜$10^{-13}$/日のオーダーで、水晶USOより2〜3桁小さいものの、ゼロではないため二次項が現に運用されています。
• DSNのSCLK-SCET相関ファイル: 深宇宙探査機の搭載時計(SCLK)と地上時刻(SCET)の対応は、往復測距とテレメトリタイムスタンプから定期的に推定され、NAIF/SPICEのSCLKカーネルとして配布されます。カーネルの中身は区分的な線形(+レート)係数の列であり、時計モデルの更新履歴そのものです。較正セッションの頻度は要求時刻精度とUSOのドリフト率から逆算して計画され、ミッションによって週次〜月次程度で更新されます。
• カッシーニ・ニューホライズンズのUSO実測: 重力科学・電波掩蔽観測にUSOを用いたカッシーニでは、USOの周波数が地上較正データとして継続的に監視され、エージングによるドリフトが年オーダーで追跡されました。実測されるドリフト率はまさに $10^{-12}$〜$10^{-11}$/日のオーダーで、稼働年数とともに減速していく対数則的な振る舞いを示します。ニューホライズンズはOne-wayドップラー航法を採用したため、USOの時計モデル精度がそのまま航法精度に直結する設計であり、巡航中の定期的な時刻比較パスがミッション計画に組み込まれていました。
• 規格・ハンドブック: 周波数ドリフト(エージング)の定義と測定法はIEEE Std 1139に、宇宙用水晶発振器の環境感度(放射線・温度・加速度)の整理は米陸軍研究所のVigのチュートリアルにまとまっています。DSN側の周波数・時刻系の要求はDSN 810-005のFTSモジュールに規定されており、地上局側の水素メーザーにも(桁こそ小さいものの)同じ3状態モデルによる監視が適用されています。地上の時刻標準の世界でも、BIPMがTAI計算に用いる各機関時計の重み付けは、まさにこの種の時計モデルの予測残差に基づいています。

演習問題

1. あるUSOのドリフト率が $D = 5 \times 10^{-12}$/日($= 5.79 \times 10^{-17}\ \text{s}^{-1}$)、較正直後の時刻オフセット・周波数オフセットは完全にゼロに較正されたとする。(a) 60日後の相対周波数偏差 $y$、(b) 60日間で蓄積する時刻誤差 $\frac{1}{2}Dt^2$、(c) その時刻誤差がOne-wayレンジ観測に与える等価距離誤差 $c \cdot \frac{1}{2}Dt^2$ をそれぞれ求めよ。
2. 時刻比較を10日おきに5回($t_i = 0, 10, 20, 30, 40$ 日)実施し、オフセット測定値 $\tilde{x}_i = \{0.0,\ 12.1,\ 28.3,\ 48.2,\ 72.1\}$ マイクロ秒を得た。二次モデル $\tilde{x} = x_0 + y_0 t + \frac{1}{2}D t^2$ を最小二乗フィットし(手計算でもスプレッドシートでもよい)、$y_0$ と $D$ を推定せよ。単位換算($\mu$s/日 → 無次元、/日 → /s)に注意すること。
3. 本文中の推定分散のスケーリング $\mathrm{Var}(\hat{D}) \propto \sigma^2 T^{-4}/N$ を前提に、観測区間を2倍(測定回数も2倍)にするとドリフト率の推定精度(標準偏差)は何倍改善するか。また、区間を無制限に延ばしても精度改善が続かない理由を、ランダムウォーク周波数雑音とエージングの対数則の2つの観点から説明せよ。
4. One-wayドップラー観測だけで軌道決定を行う場合、時計の周波数オフセット $y_0$ の誤差シグネチャは軌道側のどの量の誤差と見分けにくいか。また、ドリフト率 $D$ の誤差シグネチャはどうか。それぞれ本文の観測方程式をもとに考察し、この縮退を緩和するためにTwo-way観測やDDORのような独立観測を混ぜることがなぜ有効かを説明せよ。

まとめと次回予告

この回では、アラン分散が捉える短期の統計的揺らぎとは別の顔——数ヶ月〜数年スケールの決定論的な長期ドリフト——を扱いました。時刻オフセット $x_0$・周波数オフセット $y_0$・ドリフト率 $D$ からなる二次多項式時計モデルは、GPSの放送時計補正からDSNのSCLK-SCET相関、軌道決定フィルタの3状態時計モデルまで、時刻を扱うあらゆる現場に共通する語彙です。ドリフトの物理的起源(水晶のエージング、放射線損傷、熱サイクル応力)を押さえたうえで、双方向時刻比較の時系列への最小二乗フィットという推定の実務手順、そして未較正ドリフトが軌道解へ系統誤差として漏れ込む構造と、状態ベクトル拡張による同時推定という対策を確認しました。「短期はアラン分散で統計的に、長期は多項式で決定論的に」という役割分担と、その両者がプロセス雑音の設定を介してカルマンフィルタの中で再会する、という全体像を掴んでもらえたなら、この回の目標は達成です。

次回は視点を変えて、1つのアンテナで複数の探査機を同時に追跡する技術を扱う予定です。火星周回機が何機も同じ視野に入る現代、DSNの貴重なアンテナ時間を複数ミッションで共有するMSPA(Multiple Spacecraft Per Aperture)のような運用がどんな信号処理と周波数計画で成立しているのか——本日の時計の話とも無縁ではない、周波数資源の共有の話です。

参考文献

• IEEE Std 1139-2008, IEEE Standard Definitions of Physical Quantities for Fundamental Frequency and Time Metrology — Random Instabilities
• J. R. Vig, “Introduction to Quartz Frequency Standards,” US Army Research Laboratory, SLCET-TR-92-1(エージング・放射線感度の章)
• IS-GPS-200, NAVSTAR GPS Space Segment/Navigation User Segment Interfaces(放送時計補正係数 $a_{f0}, a_{f1}, a_{f2}$ の規定)
• J. W. Chaffee, “Relating the Allan Variance to the Diffusion Coefficients of a Linear Stochastic Differential Equation Model for Precision Oscillators,” IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 34, no. 6 (1987)
• L. Galleani, “A Tutorial on the Two-State Model of the Atomic Clock Noise,” Metrologia, vol. 45, no. 6 (2008)
• DSN Telecommunications Link Design Handbook, DSN No. 810-005(Frequency and Timing Subsystem に関するモジュール)
• NAIF/SPICE Documentation, SCLK Required Reading(SCLKカーネルと時刻相関係数)
• T. D. Moyer, Formulation for Observed and Computed Values of Deep Space Network Data Types for Navigation, JPL Deep-Space Communications and Navigation Series