前回、相関器出力を線形近似する『おおまかな導出』で測距分散の式を得ました。この回では推定理論の一般論であるクラメール・ラオ下限(CRLB)をコーシー・シュワルツの不等式から厳密に証明し、AWGN測距チャネルの実効帯域幅を使って遅延推定の理論限界を導出、前回の近似式との対応関係を確認します。ドップラー推定への双対的な拡張も扱います。
この回で学ぶこと
再生測距の回 で、私たちは相関器による遅延推定の分散を
σ τ 2 ≈ T c 2 2 T ( P R / N 0 ) \sigma_\tau^2 \approx \frac{T_c^{2}}{2\,T\,(P_R/N_0)} σ τ 2 ≈ 2 T ( P R / N 0 ) T c 2
という式で導きました。しかしその導出は、相関器出力を真の遅延近傍で線形近似し、PN符号の自己相関ピークの傾きを 1 / T c 1/T_c 1/ T c とみなす、という特定の受信機構造(相関器・線形判別器)を仮定した「おおまかな」計算 でした。あの時点では、「この式は本当に測距精度の限界を表しているのか、それとも別の受信機を使えばもっと良い精度が出せるのか」という問いには答えられていませんでした。
この回では、その問いに厳密に答えるための一般論――推定理論(estimation theory) ――を導入します。特定の受信機構造に依存せず、「どんな不偏推定量を使っても、これより良い精度は原理的に出せない」という理論的な下限 を与えるのが、この回の主役である**クラメール・ラオ下限(Cramér-Rao Lower Bound, CRLB)**です。ゴールは4つです。
フィッシャー情報量 を定義し、コーシー・シュワルツの不等式を使って、不偏推定量の分散が 1 / ( フィッシャー情報量 ) 1/(\text{フィッシャー情報量}) 1/ ( フィッシャー情報量 ) を下回れないことを証明 する。
AWGN測距チャネルにおける時間遅延(測距)推定のフィッシャー情報量を、受信信号の実効帯域幅 という概念を使って具体的に計算し、σ τ 2 ≥ 1 / ( 8 π 2 β 2 SNR ) \sigma_\tau^2 \geq 1/(8\pi^2\beta^2\,\text{SNR}) σ τ 2 ≥ 1/ ( 8 π 2 β 2 SNR ) という形の下限を導く。
この下限が、前回の「おおまかな」近似式 σ τ 2 ≈ T c 2 / ( 2 T ( P R / N 0 ) ) \sigma_\tau^2 \approx T_c^2/(2T(P_R/N_0)) σ τ 2 ≈ T c 2 / ( 2 T ( P R / N 0 )) とどう整合するのかを確認し、「チップレートが実効帯域幅の目安になっている」という対応関係を明らかにする。
ドップラー(周波数)推定にも同様のCRLBが存在し、今度は観測時間の長さ が精度を支配することを、時間・周波数の双対性から示す。
これは単なる数学的な厳密化ではありません。CRLBは実際のミッション設計において、「測距コードの帯域幅をどこまで広げれば要求精度に届くか」「ドップラー追尾の積分時間をどれだけ延ばせば重力科学実験に使えるデータになるか」を判断するための、理論的な性能限界の指標として使われています。
直感的な導入 — 「どんな定規を使っても超えられない目盛りの細かさ」
前回、測距精度がPN符号の自己相関ピークの鋭さで決まる様子を「距離という物差しの目盛りの細かさ」というアナロジーで説明しました。しかしそこで暗黙に使っていたのは、「相関器を使い、ピーク付近を線形近似する」という特定の測り方 でした。もし別の測り方(たとえば非線形な最尤推定器や、もっと複雑な統計処理)を使えば、もっと細かい目盛りが読めるのではないか——という疑問が当然浮かびます。
推定理論はこの疑問に、驚くほどきれいな答えを用意しています。答えは「測り方に依らない、データそのものが持っている情報量によって決まる、絶対的な下限が存在する 」というものです。受信信号 r ( t ) r(t) r ( t ) の中に、遅延 τ \tau τ についての情報がどれだけ含まれているかは、受信機の設計とは無関係にデータの統計的性質だけで決まります。この「データが持つ情報量」を定量化したものがフィッシャー情報量 であり、どんなに賢い推定器を設計しても、その分散はフィッシャー情報量の逆数を下回ることは決してできません。これがクラメール・ラオ下限の主張です。
言い換えると、前回導いた相関器の分散公式は「ある具体的な測り方の性能 」であり、これから導くCRLBは「どんな測り方をしても超えられない壁 」です。壁に対して実際の測り方がどれだけ近づけているか(効率, efficiency)を評価できて初めて、「この受信機は十分に良い設計なのか、まだ改善の余地があるのか」を判断できます。以下、この壁を数式で厳密に立てていきます。
フィッシャー情報量とクラメール・ラオの定理
設定
未知の決定論的パラメータ θ \theta θ を、確率密度 p ( x ; θ ) p(x;\theta) p ( x ; θ ) に従って生成される観測データ x x x から推定する問題を考えます。推定量(estimator)を θ ^ ( x ) \hat\theta(x) θ ^ ( x ) と書き、これが**不偏(unbiased)**であるとします。
E θ [ θ ^ ( x ) ] = ∫ θ ^ ( x ) p ( x ; θ ) d x = θ ( すべての θ について ) E_\theta[\hat\theta(x)] = \int \hat\theta(x)\, p(x;\theta)\, dx = \theta \qquad (\text{すべての } \theta \text{ について}) E θ [ θ ^ ( x )] = ∫ θ ^ ( x ) p ( x ; θ ) d x = θ ( すべての θ について )
フィッシャー情報量の定義
対数尤度 ln p ( x ; θ ) \ln p(x;\theta) ln p ( x ; θ ) の θ \theta θ に関する偏微分をスコア関数 と呼びます。フィッシャー情報量 は、スコア関数の分散(平均はゼロになることが以下で分かります)として定義されます。
I ( θ ) ≡ E [ ( ∂ ∂ θ ln p ( x ; θ ) ) 2 ] I(\theta) \equiv E\!\left[\left(\frac{\partial}{\partial\theta}\ln p(x;\theta)\right)^{2}\right] I ( θ ) ≡ E [ ( ∂ θ ∂ ln p ( x ; θ ) ) 2 ]
まず、スコア関数の期待値が恒等的にゼロであることを確認しておきます。確率密度の正規化条件 ∫ p ( x ; θ ) d x = 1 \int p(x;\theta)\,dx = 1 ∫ p ( x ; θ ) d x = 1 を θ \theta θ で微分すると(積分と微分の順序交換ができる、という正則条件のもとで)、
∫ ∂ p ( x ; θ ) ∂ θ d x = 0 ⟹ ∫ p ( x ; θ ) ∂ ln p ( x ; θ ) ∂ θ d x = 0 ⟹ E [ ∂ ln p ∂ θ ] = 0 \int \frac{\partial p(x;\theta)}{\partial\theta}\, dx = 0
\quad\Longrightarrow\quad
\int p(x;\theta)\,\frac{\partial \ln p(x;\theta)}{\partial\theta}\, dx = 0
\quad\Longrightarrow\quad
E\!\left[\frac{\partial \ln p}{\partial\theta}\right] = 0 ∫ ∂ θ ∂ p ( x ; θ ) d x = 0 ⟹ ∫ p ( x ; θ ) ∂ θ ∂ ln p ( x ; θ ) d x = 0 ⟹ E [ ∂ θ ∂ ln p ] = 0
を得ます(途中 ∂ p / ∂ θ = p ⋅ ∂ ln p / ∂ θ \partial p/\partial\theta = p\cdot \partial \ln p/\partial\theta ∂ p / ∂ θ = p ⋅ ∂ ln p / ∂ θ を使いました)。もう一度 θ \theta θ で微分すると、フィッシャー情報量のもう一つの等価な表現 が得られます。
∂ 2 p ∂ θ 2 = p [ ( ∂ ln p ∂ θ ) 2 + ∂ 2 ln p ∂ θ 2 ] \frac{\partial^2 p}{\partial\theta^2} = p\left[\left(\frac{\partial \ln p}{\partial\theta}\right)^{2} + \frac{\partial^2 \ln p}{\partial\theta^2}\right] ∂ θ 2 ∂ 2 p = p [ ( ∂ θ ∂ ln p ) 2 + ∂ θ 2 ∂ 2 ln p ]
この両辺を積分すると左辺は(正規化条件の2階微分より)ゼロになるので、
I ( θ ) = E [ ( ∂ ln p ∂ θ ) 2 ] = − E [ ∂ 2 ln p ∂ θ 2 ] I(\theta) = E\!\left[\left(\frac{\partial \ln p}{\partial\theta}\right)^{2}\right] = -E\!\left[\frac{\partial^2 \ln p}{\partial\theta^2}\right] I ( θ ) = E [ ( ∂ θ ∂ ln p ) 2 ] = − E [ ∂ θ 2 ∂ 2 ln p ]
という関係が成り立ちます。2階微分形は、対数尤度関数の「曲がり具合(尖り具合)」が情報量そのものであることを示しています。尤度がパラメータの変化に敏感に(鋭く)反応するほど、θ \theta θ についての情報がたくさん含まれているというわけです。以下ではこの2つの表現を状況に応じて使い分けます。
クラメール・ラオの定理の証明
定理 : 不偏推定量 θ ^ ( x ) \hat\theta(x) θ ^ ( x ) の分散は、次の下限を下回ることができない。
Var ( θ ^ ) ≥ 1 I ( θ ) \boxed{\;\text{Var}(\hat\theta) \geq \frac{1}{I(\theta)}\;} Var ( θ ^ ) ≥ I ( θ ) 1
証明 : 不偏性の条件 ∫ θ ^ ( x ) p ( x ; θ ) d x = θ \int \hat\theta(x)\,p(x;\theta)\,dx = \theta ∫ θ ^ ( x ) p ( x ; θ ) d x = θ を θ \theta θ で微分します。
∫ θ ^ ( x ) ∂ p ( x ; θ ) ∂ θ d x = 1 ⟹ ∫ θ ^ ( x ) p ( x ; θ ) ∂ ln p ∂ θ d x = 1 ⟹ E [ θ ^ ( x ) ⋅ ∂ ln p ∂ θ ] = 1 \int \hat\theta(x)\, \frac{\partial p(x;\theta)}{\partial\theta}\, dx = 1
\quad\Longrightarrow\quad
\int \hat\theta(x)\, p(x;\theta)\, \frac{\partial \ln p}{\partial\theta}\, dx = 1
\quad\Longrightarrow\quad
E\!\left[\hat\theta(x)\cdot\frac{\partial \ln p}{\partial\theta}\right] = 1 ∫ θ ^ ( x ) ∂ θ ∂ p ( x ; θ ) d x = 1 ⟹ ∫ θ ^ ( x ) p ( x ; θ ) ∂ θ ∂ ln p d x = 1 ⟹ E [ θ ^ ( x ) ⋅ ∂ θ ∂ ln p ] = 1
前節で示した通り E [ ∂ ln p / ∂ θ ] = 0 E[\partial \ln p/\partial\theta] = 0 E [ ∂ ln p / ∂ θ ] = 0 なので、両辺から θ ⋅ E [ ∂ ln p / ∂ θ ] = 0 \theta \cdot E[\partial \ln p/\partial\theta] = 0 θ ⋅ E [ ∂ ln p / ∂ θ ] = 0 を引いても等式は変わらず、
E [ ( θ ^ ( x ) − θ ) ⋅ ∂ ln p ∂ θ ] = 1 E\!\left[\big(\hat\theta(x) - \theta\big)\cdot\frac{\partial \ln p}{\partial\theta}\right] = 1 E [ ( θ ^ ( x ) − θ ) ⋅ ∂ θ ∂ ln p ] = 1
が成り立ちます。ここでコーシー・シュワルツの不等式 ( E [ U V ] ) 2 ≤ E [ U 2 ] E [ V 2 ] \big(E[UV]\big)^2 \leq E[U^2]\,E[V^2] ( E [ U V ] ) 2 ≤ E [ U 2 ] E [ V 2 ] を、U = θ ^ ( x ) − θ U = \hat\theta(x)-\theta U = θ ^ ( x ) − θ 、V = ∂ ln p / ∂ θ V = \partial \ln p/\partial\theta V = ∂ ln p / ∂ θ に適用すると、
1 = ( E [ ( θ ^ − θ ) V ] ) 2 ≤ E [ ( θ ^ − θ ) 2 ] ⋅ E [ V 2 ] = Var ( θ ^ ) ⋅ I ( θ ) 1 = \Big(E\big[(\hat\theta-\theta)V\big]\Big)^{2} \;\leq\; E\big[(\hat\theta-\theta)^2\big]\cdot E[V^2] \;=\; \text{Var}(\hat\theta)\cdot I(\theta) 1 = ( E [ ( θ ^ − θ ) V ] ) 2 ≤ E [ ( θ ^ − θ ) 2 ] ⋅ E [ V 2 ] = Var ( θ ^ ) ⋅ I ( θ )
両辺を I ( θ ) I(\theta) I ( θ ) で割れば、求める不等式 Var ( θ ^ ) ≥ 1 / I ( θ ) \text{Var}(\hat\theta) \geq 1/I(\theta) Var ( θ ^ ) ≥ 1/ I ( θ ) が得られます。■ \blacksquare ■
コーシー・シュワルツの不等式の等号成立条件(2つの確率変数が比例関係にあること、U = c ( θ ) V U = c(\theta)\,V U = c ( θ ) V )から、CRLBの等号が成立するのは ∂ ln p / ∂ θ = c ( θ ) ( θ ^ ( x ) − θ ) \partial \ln p/\partial\theta = c(\theta)\,(\hat\theta(x)-\theta) ∂ ln p / ∂ θ = c ( θ ) ( θ ^ ( x ) − θ ) という特別な関係が成り立つ場合に限られることも分かります。この条件を満たす推定量を**有効推定量(efficient estimator)**と呼び、後で見るように、AWGN通信路での最尤推定はこの条件を(高SNR極限で)満たします。
独立ガウス観測に対するフィッシャー情報量(補題)
以下の測距・ドップラーへの応用で繰り返し使う、簡単だが重要な補題を先に示しておきます。観測値 x 1 , … , x N x_1,\dots,x_N x 1 , … , x N が独立に、平均 μ i ( θ ) \mu_i(\theta) μ i ( θ ) ・分散 σ 2 \sigma^2 σ 2 (共通)の正規分布に従うとします。
ln p ( x ; θ ) = − 1 2 σ 2 ∑ i = 1 N ( x i − μ i ( θ ) ) 2 + const \ln p(x;\theta) = -\frac{1}{2\sigma^2}\sum_{i=1}^{N}\big(x_i - \mu_i(\theta)\big)^2 + \text{const} ln p ( x ; θ ) = − 2 σ 2 1 i = 1 ∑ N ( x i − μ i ( θ ) ) 2 + const
θ \theta θ で2階微分すると、
∂ 2 ln p ∂ θ 2 = 1 σ 2 ∑ i [ − μ i ′ ( θ ) 2 + ( x i − μ i ( θ ) ) μ i ′ ′ ( θ ) ] \frac{\partial^2 \ln p}{\partial\theta^2} = \frac{1}{\sigma^2}\sum_i \Big[-\mu_i'(\theta)^2 + \big(x_i-\mu_i(\theta)\big)\mu_i''(\theta)\Big] ∂ θ 2 ∂ 2 ln p = σ 2 1 i ∑ [ − μ i ′ ( θ ) 2 + ( x i − μ i ( θ ) ) μ i ′′ ( θ ) ]
期待値を取ると E [ x i − μ i ( θ ) ] = 0 E[x_i-\mu_i(\theta)]=0 E [ x i − μ i ( θ )] = 0 なので第2項は消え、
I ( θ ) = − E [ ∂ 2 ln p ∂ θ 2 ] = 1 σ 2 ∑ i = 1 N μ i ′ ( θ ) 2 \boxed{\;I(\theta) = -E\!\left[\frac{\partial^2\ln p}{\partial\theta^2}\right] = \frac{1}{\sigma^2}\sum_{i=1}^N \mu_i'(\theta)^2\;} I ( θ ) = − E [ ∂ θ 2 ∂ 2 ln p ] = σ 2 1 i = 1 ∑ N μ i ′ ( θ ) 2
という、平均がパラメータ θ \theta θ にどれだけ敏感に反応するか(μ i ′ ( θ ) \mu_i'(\theta) μ i ′ ( θ ) )の2乗和で情報量が決まる、というシンプルな結果が得られます。これを次節でそのまま使います。
AWGN測距チャネルにおけるCRLB: 実効帯域幅の登場
離散化モデルと連続極限
受信機が受け取る測距信号を、前回 と同様に
r ( t ) = s ( t − τ ) + n ( t ) r(t) = s(t-\tau) + n(t) r ( t ) = s ( t − τ ) + n ( t )
とします。s ( t ) s(t) s ( t ) は既知の(エネルギー E r ≡ ∫ s 2 ( t ) d t E_r \equiv \int s^2(t)\,dt E r ≡ ∫ s 2 ( t ) d t を持つ)受信波形、n ( t ) n(t) n ( t ) は片側スペクトル密度 N 0 N_0 N 0 (両側 N 0 / 2 N_0/2 N 0 /2 )の白色ガウス雑音、τ \tau τ が推定したい遅延です。
連続時間の雑音場を厳密に扱うには関数解析的な議論が必要になりますが、前回の離散チップモデルと同じ発想で、時間刻み Δ t \Delta t Δ t で受信波形をサンプリングして近似することにします。サンプル値は
r i = s ( i Δ t − τ ) + n i , n i ∼ N ( 0 , σ 2 ) , σ 2 = N 0 2 Δ t r_i = s(i\Delta t - \tau) + n_i, \qquad n_i \sim \mathcal{N}(0,\sigma^2), \qquad \sigma^2 = \frac{N_0}{2\Delta t} r i = s ( i Δ t − τ ) + n i , n i ∼ N ( 0 , σ 2 ) , σ 2 = 2Δ t N 0
という、前節の補題がそのまま使える独立ガウスサンプル列になります(雑音分散が Δ t \Delta t Δ t に反比例するのは、白色雑音を細かく刻むほど各サンプルの実効帯域が広がり分散が増えるためで、ranging-regen.mdのチップ雑音モデル σ n 2 = N 0 / ( 2 T c ) \sigma_n^2 = N_0/(2T_c) σ n 2 = N 0 / ( 2 T c ) の T c → Δ t T_c \to \Delta t T c → Δ t 版そのものです)。ここで μ i ( τ ) = s ( i Δ t − τ ) \mu_i(\tau) = s(i\Delta t-\tau) μ i ( τ ) = s ( i Δ t − τ ) なので、μ i ′ ( τ ) = − s ′ ( i Δ t − τ ) \mu_i'(\tau) = -s'(i\Delta t-\tau) μ i ′ ( τ ) = − s ′ ( i Δ t − τ ) 。補題を適用すると、
I ( τ ) = 1 σ 2 ∑ i s ′ ( i Δ t − τ ) 2 = 2 Δ t N 0 ∑ i s ′ ( i Δ t − τ ) 2 I(\tau) = \frac{1}{\sigma^2}\sum_i s'(i\Delta t-\tau)^2 = \frac{2\Delta t}{N_0}\sum_i s'(i\Delta t-\tau)^2 I ( τ ) = σ 2 1 i ∑ s ′ ( i Δ t − τ ) 2 = N 0 2Δ t i ∑ s ′ ( i Δ t − τ ) 2
Δ t → 0 \Delta t \to 0 Δ t → 0 の連続極限を取ると、和がリーマン和として積分に収束し、
I ( τ ) = 2 N 0 ∫ − ∞ ∞ [ s ′ ( t − τ ) ] 2 d t I(\tau) = \frac{2}{N_0}\int_{-\infty}^{\infty} \big[s'(t-\tau)\big]^2\, dt I ( τ ) = N 0 2 ∫ − ∞ ∞ [ s ′ ( t − τ ) ] 2 d t
という、サンプリングの刻み方に依存しない、遅延推定のフィッシャー情報量の式が得られます。s ′ s' s ′ の積分は τ \tau τ の値によらず一定なので、以下 τ = 0 \tau=0 τ = 0 として ∫ [ s ′ ( t ) ] 2 d t \int [s'(t)]^2 dt ∫ [ s ′ ( t ) ] 2 d t と書きます。
実効帯域幅による書き換え
この積分を周波数領域の言葉に翻訳します。s ( t ) s(t) s ( t ) のフーリエ変換を S ( f ) = ∫ s ( t ) e − j 2 π f t d t S(f) = \int s(t)e^{-j2\pi f t}dt S ( f ) = ∫ s ( t ) e − j 2 π f t d t とすると、微分の性質より s ′ ( t ) s'(t) s ′ ( t ) のフーリエ変換は j 2 π f S ( f ) j2\pi f\, S(f) j 2 π f S ( f ) です。パーセバルの定理を使うと、
∫ [ s ′ ( t ) ] 2 d t = ∫ ∣ 2 π f S ( f ) ∣ 2 d f = 4 π 2 ∫ f 2 ∣ S ( f ) ∣ 2 d f \int \big[s'(t)\big]^2\, dt = \int \big|2\pi f\, S(f)\big|^{2}\, df = 4\pi^2 \int f^{2}\,|S(f)|^{2}\, df ∫ [ s ′ ( t ) ] 2 d t = ∫ 2 π f S ( f ) 2 df = 4 π 2 ∫ f 2 ∣ S ( f ) ∣ 2 df
一方、信号のエネルギーも同じくパーセバルの定理から E r = ∫ s 2 ( t ) d t = ∫ ∣ S ( f ) ∣ 2 d f E_r = \int s^2(t)\,dt = \int |S(f)|^2\,df E r = ∫ s 2 ( t ) d t = ∫ ∣ S ( f ) ∣ 2 df です。そこで、信号のエネルギースペクトル密度 ∣ S ( f ) ∣ 2 |S(f)|^2 ∣ S ( f ) ∣ 2 の周波数に関する**2次モーメント(規格化した分散)**として、実効帯域幅(RMS帯域幅) β \beta β を次のように定義します。
β 2 ≡ ∫ f 2 ∣ S ( f ) ∣ 2 d f ∫ ∣ S ( f ) ∣ 2 d f = ∫ f 2 ∣ S ( f ) ∣ 2 d f E r \beta^{2} \equiv \frac{\displaystyle\int f^{2}\,|S(f)|^{2}\, df}{\displaystyle\int |S(f)|^{2}\, df} = \frac{\displaystyle\int f^{2}\,|S(f)|^{2}\, df}{E_r} β 2 ≡ ∫ ∣ S ( f ) ∣ 2 df ∫ f 2 ∣ S ( f ) ∣ 2 df = E r ∫ f 2 ∣ S ( f ) ∣ 2 df
β \beta β は「信号のエネルギーが周波数軸上でどれだけ広がっているか」を表す量で、単位はHzです。この定義を使うと ∫ [ s ′ ( t ) ] 2 d t = 4 π 2 β 2 E r \int [s'(t)]^2 dt = 4\pi^2\beta^2 E_r ∫ [ s ′ ( t ) ] 2 d t = 4 π 2 β 2 E r となるので、フィッシャー情報量は
I ( τ ) = 2 N 0 ⋅ 4 π 2 β 2 E r = 8 π 2 β 2 ( E r N 0 ) I(\tau) = \frac{2}{N_0}\cdot 4\pi^2\beta^2 E_r = 8\pi^2\beta^2\left(\frac{E_r}{N_0}\right) I ( τ ) = N 0 2 ⋅ 4 π 2 β 2 E r = 8 π 2 β 2 ( N 0 E r )
エネルギー対雑音密度比を SNR ≡ E r / N 0 \text{SNR} \equiv E_r/N_0 SNR ≡ E r / N 0 と定義すれば、
I ( τ ) = 8 π 2 β 2 ⋅ SNR \boxed{\;I(\tau) = 8\pi^2\beta^2\cdot \text{SNR}\;} I ( τ ) = 8 π 2 β 2 ⋅ SNR
というきれいな形にまとまります。クラメール・ラオの定理をそのまま適用すると、遅延推定の分散に対する理論限界が得られます。
σ τ 2 ≥ 1 8 π 2 β 2 SNR \boxed{\;\sigma_\tau^{2} \;\geq\; \frac{1}{8\pi^{2}\beta^{2}\,\text{SNR}}\;} σ τ 2 ≥ 8 π 2 β 2 SNR 1
この式が意味することは直感的です。信号のエネルギーが周波数軸上に広く分布している(実効帯域幅 β \beta β が大きい)ほど、遅延方向の「傾き」の情報が多く、精度の理論限界は良くなる(分散の下限が小さくなる)。 これは前回、自己相関ピークの傾き 1 / T c 1/T_c 1/ T c が精度を支配すると述べたことの、厳密な一般化になっています。
ranging-regen.mdの近似式との整合性
前回導いた式を思い出しましょう。積分時間 T T T 、測距SNR ( P R / N 0 ) (P_R/N_0) ( P R / N 0 ) のもとで、
σ τ 2 ≈ T c 2 2 T ( P R / N 0 ) \sigma_\tau^2 \approx \frac{T_c^2}{2\,T\,(P_R/N_0)} σ τ 2 ≈ 2 T ( P R / N 0 ) T c 2
でした。この式とCRLBを比較するために、まず両者の「SNR」の定義を揃えます。今回の SNR = E r / N 0 \text{SNR}=E_r/N_0 SNR = E r / N 0 は積分時間 T T T ぶんの受信エネルギーなので、E r = P R T E_r = P_R\, T E r = P R T とすれば SNR = T ( P R / N 0 ) \text{SNR} = T\,(P_R/N_0) SNR = T ( P R / N 0 ) となり、CRLBは
σ τ 2 ≥ 1 8 π 2 β 2 T ( P R / N 0 ) \sigma_\tau^2 \;\geq\; \frac{1}{8\pi^2\beta^2\, T\,(P_R/N_0)} σ τ 2 ≥ 8 π 2 β 2 T ( P R / N 0 ) 1
と書き直せます。前回の近似式・今回のCRLBはともに「T c 2 T_c^2 T c 2 または 1 / β 2 1/\beta^2 1/ β 2 」「1 / T 1/T 1/ T 」「1 / ( P R / N 0 ) 1/(P_R/N_0) 1/ ( P R / N 0 ) 」という同じ3つの依存性を持っており、これは偶然ではありません。両者を係数だけ比較すると、
T c 2 2 ⟷ 1 8 π 2 β 2 \frac{T_c^2}{2} \quad \longleftrightarrow \quad \frac{1}{8\pi^2\beta^2} 2 T c 2 ⟷ 8 π 2 β 2 1
という対応が見えてきます。もし実効帯域幅 β \beta β を、チップ周期 T c T_c T c を使って β ≈ 1 / ( 2 T c ) \beta \approx 1/(2T_c) β ≈ 1/ ( 2 T c ) という(オーダー評価としての)目安で置くと、
1 8 π 2 β 2 = 1 8 π 2 ⋅ 4 T c 2 = T c 2 2 π 2 ≈ T c 2 19.7 \frac{1}{8\pi^2\beta^2} = \frac{1}{8\pi^2}\cdot 4T_c^2 = \frac{T_c^2}{2\pi^2} \approx \frac{T_c^2}{19.7} 8 π 2 β 2 1 = 8 π 2 1 ⋅ 4 T c 2 = 2 π 2 T c 2 ≈ 19.7 T c 2
となり、前回の近似係数 T c 2 / 2 T_c^2/2 T c 2 /2 より小さい(=CRLBの方が良い、より小さい分散を許す)値になります。これは理にかなっています。CRLBはあらゆる不偏推定量の中での理論上の最良値であり、前回扱った「相関器出力を線形近似する」という具体的な(かつ最適とは限らない)受信機構造の分散は、CRLBを下回ることは決してできません。 実際、
σ τ , 相関器 2 ≈ T c 2 2 T ( P R / N 0 ) ≥ 1 8 π 2 β 2 T ( P R / N 0 ) ≈ σ τ , CRLB 2 \sigma_{\tau,\text{相関器}}^2 \approx \frac{T_c^2}{2T(P_R/N_0)} \;\geq\; \frac{1}{8\pi^2\beta^2\,T(P_R/N_0)} \approx \sigma_{\tau,\text{CRLB}}^2 σ τ , 相関器 2 ≈ 2 T ( P R / N 0 ) T c 2 ≥ 8 π 2 β 2 T ( P R / N 0 ) 1 ≈ σ τ , CRLB 2
という不等式が保たれており、両者の係数がおよそ1桁の範囲(約10倍)に収まっているのは、相関器ベースの推定が理論限界からそう遠くない、実用上「そこそこ効率の良い」推定方式であることを示しています。
ここで一点、数学的に注意すべき点があります。理想的な矩形パルス のチップ(ranging-regen.mdで自己相関を三角形近似したときの前提)は、波形の立ち上がり・立ち下がりが不連続なため、微分 s ′ ( t ) s'(t) s ′ ( t ) がデルタ関数的に鋭くなり、実効帯域幅 β \beta β は数学的には無限大に発散 してしまいます。これは「理想的に帯域無限のチップを使えば測距精度はいくらでも良くできる」という(非現実的な)結論を意味しており、CRLBの式が病的にゼロへ潰れてしまう例です。実際の測距波形は、送信機のフィルタリングやチャネルの帯域制限によって必ず有限の帯域幅に収まるため、実効帯域幅 β \beta β は送受信フィルタの設計によって決まる有限な値になります。チップレート R c = 1 / T c R_c=1/T_c R c = 1/ T c は、この現実的な実効帯域幅の目安(オーダー)として使える 、というのが両者の対応関係の正しい読み方です。CCSDS規格の測距波形設計(送受信フィルタの帯域指定)は、まさにこの「チップレートと実効帯域幅の対応」を踏まえて行われています。
ドップラー(周波数)推定のCRLB: 観測時間の双対性
同じ枠組みは、時間遅延だけでなくドップラーシフト(周波数)の推定にも適用できます。フーリエ変換において時間と周波数は対称な役割を持つ(時間・周波数の双対性 )ため、遅延推定の議論をそのまま「時間軸と周波数軸を入れ替えて」繰り返すことができます。
観測信号を、未知の周波数オフセット f d f_d f d (ドップラーシフト)を持つ複素基底帯域信号として
r ( t ) = s ( t ) e j 2 π f d t + n ( t ) , t ∈ [ − T 2 , T 2 ] r(t) = s(t)\,e^{j2\pi f_d t} + n(t), \qquad t \in \left[-\frac{T}{2}, \frac{T}{2}\right] r ( t ) = s ( t ) e j 2 π f d t + n ( t ) , t ∈ [ − 2 T , 2 T ]
とモデル化します(観測窓の中心を、信号エネルギーの「重心」に取っておくと、遅延推定のときの τ \tau τ 微分と対称な形で交差項がきれいに消えます)。遅延推定とまったく同じ手順(離散化 → ガウス補題 → 連続極限 → パーセバルの定理)を、τ \tau τ の代わりに f d f_d f d 、周波数領域の代わりに時間領域で繰り返すと、
I ( f d ) = 2 N 0 ∫ ( 2 π t ) 2 ∣ s ( t ) ∣ 2 d t = 8 π 2 T rms 2 ⋅ SNR I(f_d) = \frac{2}{N_0}\int (2\pi t)^2\, |s(t)|^2\, dt = 8\pi^2\, T_{\text{rms}}^2\cdot \text{SNR} I ( f d ) = N 0 2 ∫ ( 2 π t ) 2 ∣ s ( t ) ∣ 2 d t = 8 π 2 T rms 2 ⋅ SNR
という、遅延推定の式ときれいに対をなす結果が得られます。ここで T rms 2 T_{\text{rms}}^2 T rms 2 は、信号エネルギー ∣ s ( t ) ∣ 2 |s(t)|^2 ∣ s ( t ) ∣ 2 の**時間軸上の2次モーメント(実効的な観測時間の広がり)**です。
T rms 2 ≡ ∫ t 2 ∣ s ( t ) ∣ 2 d t ∫ ∣ s ( t ) ∣ 2 d t T_{\text{rms}}^{2} \equiv \frac{\int t^{2}\,|s(t)|^{2}\, dt}{\int |s(t)|^{2}\, dt} T rms 2 ≡ ∫ ∣ s ( t ) ∣ 2 d t ∫ t 2 ∣ s ( t ) ∣ 2 d t
したがって、ドップラー推定のCRLBは
σ f d 2 ≥ 1 8 π 2 T rms 2 SNR \boxed{\;\sigma_{f_d}^{2} \;\geq\; \frac{1}{8\pi^{2}\, T_{\text{rms}}^{2}\, \text{SNR}}\;} σ f d 2 ≥ 8 π 2 T rms 2 SNR 1
という、遅延推定のCRLBで β ↔ T rms \beta \leftrightarrow T_{\text{rms}} β ↔ T rms を入れ替えただけの形になります。これはまさに「遅延推定における実効帯域幅の役割を、ドップラー推定では実効観測時間が担う 」という美しい双対性を表しています。具体的に、長さ T T T の矩形観測窓(信号エネルギーが [ − T / 2 , T / 2 ] [-T/2, T/2] [ − T /2 , T /2 ] に一様に分布)の場合、一様分布の分散の公式より T rms 2 = T 2 / 12 T_{\text{rms}}^2 = T^2/12 T rms 2 = T 2 /12 となるので、
σ f d 2 ≥ 1 8 π 2 ( T 2 / 12 ) SNR = 3 2 π 2 T 2 SNR \sigma_{f_d}^{2} \geq \frac{1}{8\pi^2 (T^2/12)\,\text{SNR}} = \frac{3}{2\pi^2\, T^2\, \text{SNR}} σ f d 2 ≥ 8 π 2 ( T 2 /12 ) SNR 1 = 2 π 2 T 2 SNR 3
という具体的な下限が得られます。ここでの結論は明快です。観測(コヒーレント積分)時間 T T T を長く取るほど、ドップラー推定の分散は 1 / T 2 1/T^2 1/ T 2 で急速に改善します。 PLLの回 で扱ったループ帯域幅 B L B_L B L の設計は、雑音の取り込みと追従性のトレードオフというまた別の実用上の制約(探査機のダイナミクスへの追従)を扱うものでしたが、ここで得られたCRLBは、そうした実装上の制約を一切考えない「原理的にどこまで良くなれるか」という上限を与えるものである点に注意してください。実際の受信機(PLLやディジタル周波数推定器)の性能は、この理論限界に対してどれだけ近づけているかという効率 で評価されます。
なお、離散サンプルを使った古典的な周波数推定の文献(Rife & Boorstyn, 1974)では、サンプル数 N N N を使った 1 / N 3 1/N^3 1/ N 3 型の依存性が示されることがありますが、これは「SNRをサンプルあたりのエネルギーで定義するか、総エネルギーで定義するか」という規約の違いによるもので、本質的には今回導いた「実効観測時間が長いほど良い」という結論と矛盾しません。
実務での使われ方
CRLBは抽象的な数学的下限にとどまらず、実際の深宇宙ミッションの測距・追跡システム設計そのもの を左右する指標として使われています。
測距波形の帯域幅設計への直接的な指針 。CCSDS 401.0-B / 414.1-B が規定するPN測距コードのチップレートは、単に「速ければ速いほど良い」という漠然とした話ではなく、要求される測距精度からCRLBを逆算し、必要な実効帯域幅 β \beta β (したがってチップレート)を見積もる、というミッション設計プロセスの中で決められます。BepiColombo(ESA・JAXA共同の水星探査機)やJUICE(木星氷衛星探査機)のような、相対論的検証実験や重力科学を目的とする高精度ミッションでは、通常の航法用途よりも高いチップレートの測距モードが用意されており、これはCRLBが示す「帯域幅を広げるほど理論限界が良くなる」という関係を直接反映した設計判断です。
ドップラー追跡の積分時間設計 。カッシーニ探査機の太陽重力場・重力波探索実験や、木星探査機ジュノーの重力場マッピング(gravity science)といった実験では、コヒーレントなドップラー観測を数十分から数時間規模で継続することで、1 / T 2 1/T^2 1/ T 2 に従って周波数推定の理論限界を追い込み、探査機の軌道運動から惑星の重力ポテンシャルや相対論的効果をあぶり出しています。積分時間をどこまで延ばせば目標精度に届くかという見積もりの根拠に、今回導いたCRLBの形が使われます。
受信機の「効率」の評価基準 。DSN Telecommunications Link Design Handbook(810-005)のような設計文書では、実際に運用されている測距・ドップラー抽出アルゴリズムの達成精度を、CRLBに対する比(効率)として評価する考え方が用いられます。理論限界に対して実装がどれだけ近いかが分かれば、「アルゴリズムを改良する余地があるのか」「これ以上はハードウェアの帯域幅や積分時間を変えない限り改善できないのか」を切り分けられます。
リンク予算の中での位置づけ 。測距・ドップラーの要求精度(たとえば軌道決定に必要な精度)が先に与えられているとき、CRLBの式を逆に解けば「必要な実効帯域幅 β \beta β 」や「必要な積分時間 T T T 」が求まります。これはミッション設計の初期段階で、測距波形やアンテナ・送信機の仕様を決めるためのトップダウンな設計指標として機能します。
演習問題
クラメール・ラオの定理の証明で使ったコーシー・シュワルツの不等式の等号成立条件から、CRLBの等号が成立する「有効推定量」の条件式 ∂ ln p / ∂ θ = c ( θ ) ( θ ^ ( x ) − θ ) \partial \ln p/\partial\theta = c(\theta)(\hat\theta(x)-\theta) ∂ ln p / ∂ θ = c ( θ ) ( θ ^ ( x ) − θ ) が導かれることを確認してください。また、この節で示した独立ガウス観測に対するフィッシャー情報量の補題 I ( θ ) = ( 1 / σ 2 ) ∑ i μ i ′ ( θ ) 2 I(\theta) = (1/\sigma^2)\sum_i \mu_i'(\theta)^2 I ( θ ) = ( 1/ σ 2 ) ∑ i μ i ′ ( θ ) 2 を、最尤推定量 θ ^ M L \hat\theta_{ML} θ ^ M L がこの効率条件を満たすかどうかという観点から(定性的にでよいので)考察してください。
実効帯域幅 β = 1 \beta = 1 β = 1 MHz、SNR(= E r / N 0 =E_r/N_0 = E r / N 0 )が真数で 100 100 100 (すなわち20 dB)の測距リンクについて、遅延推定のCRLB σ τ 2 ≥ 1 / ( 8 π 2 β 2 SNR ) \sigma_\tau^2 \geq 1/(8\pi^2\beta^2\,\text{SNR}) σ τ 2 ≥ 1/ ( 8 π 2 β 2 SNR ) から σ τ \sigma_\tau σ τ の下限を求め、距離換算 σ R = ( c / 2 ) σ τ \sigma_R = (c/2)\sigma_\tau σ R = ( c /2 ) σ τ でメートル単位の下限を計算してください。
前回の演習問題3 で扱った、チップレート R c = 2 R_c = 2 R c = 2 Mchip/s(T c = 0.5 μ s T_c=0.5\ \mu\text{s} T c = 0.5 μ s )、積分時間 T = 10 T=10 T = 10 s、測距SNR ( P R / N 0 ) = 15 (P_R/N_0)=15 ( P R / N 0 ) = 15 dB-Hz のリンクを再び考えます。本文で示した目安 β ≈ 1 / ( 2 T c ) \beta \approx 1/(2T_c) β ≈ 1/ ( 2 T c ) を使ってCRLB σ τ 2 ≥ 1 / ( 8 π 2 β 2 T ( P R / N 0 ) ) \sigma_\tau^2 \geq 1/(8\pi^2\beta^2\, T(P_R/N_0)) σ τ 2 ≥ 1/ ( 8 π 2 β 2 T ( P R / N 0 )) を計算し、前回の演習で求めた相関器の分散 σ τ 2 ≈ T c 2 / ( 2 T ( P R / N 0 ) ) \sigma_\tau^2 \approx T_c^2/(2T(P_R/N_0)) σ τ 2 ≈ T c 2 / ( 2 T ( P R / N 0 )) の値と比較してください。CRLBの方が小さい(良い)値になっていることを確認し、その意味を本文の議論に沿って説明してください。
ドップラー推定のCRLB σ f d 2 ≥ 3 / ( 2 π 2 T 2 SNR ) \sigma_{f_d}^2 \geq 3/(2\pi^2 T^2\,\text{SNR}) σ f d 2 ≥ 3/ ( 2 π 2 T 2 SNR ) を使って、積分時間を T T T から 4 T 4T 4 T に延ばしたとき、周波数推定の標準偏差 σ f d \sigma_{f_d} σ f d が何分の1になるか求めてください。また、この「観測時間を延ばすほど精度が上がる」という性質と、PLLの回 で学んだ「ループ帯域幅 B L B_L B L を狭くするほど雑音耐性が上がるが探査機のダイナミクスへの追従が遅れる」というトレードオフとは、どのような関係にある(あるいはない)か、CRLBが「実装上の制約を含まない理論限界」であるという点を踏まえて自分の言葉で議論してください。
まとめと次回予告
この回では、これまで測距・ドップラー精度について「おおまかな導出」で済ませてきた分散の式を、推定理論の一般論から厳密に基礎づけました。フィッシャー情報量を定義し、コーシー・シュワルツの不等式からクラメール・ラオ下限 Var ( θ ^ ) ≥ 1 / I ( θ ) \text{Var}(\hat\theta)\geq 1/I(\theta) Var ( θ ^ ) ≥ 1/ I ( θ ) を証明し、AWGN測距チャネルの遅延推定に適用して σ τ 2 ≥ 1 / ( 8 π 2 β 2 SNR ) \sigma_\tau^2 \geq 1/(8\pi^2\beta^2\,\text{SNR}) σ τ 2 ≥ 1/ ( 8 π 2 β 2 SNR ) という理論限界を、実効帯域幅 β \beta β を軸に導出しました。この式は、前回の相関器ベースの近似式 σ τ 2 ≈ T c 2 / ( 2 T ( P R / N 0 ) ) \sigma_\tau^2 \approx T_c^2/(2T(P_R/N_0)) σ τ 2 ≈ T c 2 / ( 2 T ( P R / N 0 )) と同じ依存性の構造(チップ幅の2乗・積分時間の逆数・SNRの逆数)を持ち、チップレートが実効帯域幅の実用上の目安になっているという対応関係も確認しました。さらに時間・周波数の双対性から、ドップラー推定のCRLB σ f d 2 ≥ 3 / ( 2 π 2 T 2 SNR ) \sigma_{f_d}^2 \geq 3/(2\pi^2T^2\,\text{SNR}) σ f d 2 ≥ 3/ ( 2 π 2 T 2 SNR ) を導き、観測時間の長さが精度に効くことを数式で確認しました。
これで測距・追跡カテゴリの学習は一区切りです。トランスポンダから始まり、非再生測距・再生測距・トーン測距・DDOR・SBI、そして今回のCRLBによる理論的裏付けまで、探査機との距離・速度・角度をどう測るかという一連の技術を数式付きで追ってきました。次回からはシステム・運用 カテゴリに進み、まずはこれまで登場してきたDSN(Deep Space Network)そのものを主役に据えたDSN概論 を扱います。世界に散らばる巨大パラボラアンテナ局がどう配置され、どうスケジューリングされ、どうミッション間でリソースを分け合っているのか、通信理論の数式から少し離れて、運用というシステム全体の視点でDSNを見ていきます。
参考文献
H. L. Van Trees, Detection, Estimation, and Modulation Theory, Part I , Wiley
S. M. Kay, Fundamentals of Statistical Signal Processing: Estimation Theory , Prentice Hall
D. Rife, R. Boorstyn, “Single Tone Parameter Estimation from Discrete-Time Observations,” IEEE Transactions on Information Theory , Vol. 20, No. 5, 1974
CCSDS 414.1-B, Pseudo-Noise and Regenerative Ranging
CCSDS 401.0-B, Radio Frequency and Modulation Systems, Part 1: Earth Stations and Spacecraft
DSN Telecommunications Link Design Handbook, DSN No. 810-005(Ranging および Tracking Systems に関するモジュール)
J. H. Yuen (ed.), Deep Space Telecommunications Systems Engineering , JPL Publication 82-76