♨️

宇宙からの地球観測12章干渉SAR

2024/11/22に公開

宇宙からの地球観測 
第12章は、 ついに干渉SAR (iSAR)です。 2回撮像するとどうして高さの詳細情報がとれるのか。数学が魔法のように感じる技術です。理解したいです。

問題12.1

Xバンド,Cバンド,LバンドのSARがある。波長はそれぞれ 3cm, 5.6cm, 24cm である。今,差分干渉図でアンラップをしない場合に取りうる,地殻変動量の範囲はいくらか?

基本用語の理解

干渉とは

本文より

干渉とは,ほんの少し異なる経路を通った2つの信号を足しあわせると,経路の長さの違いを検出し,それを可視化できることである。

SAR信号の複素数表現

本文の解説に下記の式12.1の説明がありますが、最初まったく理解できません。

式12.1

S_1 = A e^{-\frac{4\pi}{\lambda} r_r j}

ここで、

  • A は振幅
  • r_r は第1軌道からの距離 (下付きrはreferenceのr)
  • \lambda は波長
  • j は虚数単位(j^2 = -1

複素数の極形式

下記の解説を読むことでやっと理解できました。
https://manabitimes.jp/math/875#3

例として 複素数 z=\sqrt{3}+i は直交座標で原点からの距離が2で、X軸との角が \frac{\pi}{6}なので、
極形式では z=2(\cos{\frac{\pi}{6}}+i\sin{\frac{\pi}{6}})と表せる
オイラーの公式 e^{i\theta}=\cos{\theta}+i\sin{\theta} より
z=r(cos\theta+i\sin\theta)=re^{i\theta} と表せるので
\sqrt{3}+i=2e^{\frac{\pi}{6}i} と表すことができる

オイラーの公式、複素指数関数、eネイピア数

上記をさらに深く理解しようと思うと下記を理解する必要がでてきます。
何となくそうなるんだなという理解でまずは頭に置いておこうと思います。

https://manabitimes.jp/math/585

e^{i\theta}=\cos{\theta}+i\sin{\theta}

https://manabitimes.jp/math/714

e=\lim_{n \to \infty}\left (1+\dfrac{1}{n} \right)^n

e=2.7182818284,,,,

信号の指数形式への変換

複素数 z は次のように指数形式で表すことができます。

z = re^{j\theta}

ここで、

  • r = \sqrt{a^2 + b^2} は絶対値(振幅)
  • \theta = \tan^{-1}\left(\frac{b}{a}\right) は位相角

この形式を使うと、信号の振幅 A と位相 \theta を分けて考えることができ、計算が簡単になります。

電波の伝搬と位相

SAR(Synthetic Aperture Radar)で用いられる電波は、一定の波長 \lambda を持っています。電波は周期的な波として伝搬し、その位相は伝搬距離によって変わります。

  1. 波長と位相の関係:

    • 電波の1周期は波長 \lambda に相当します。
    • 1周期の位相は 2\pi ラジアンです。
  2. 距離と位相の変化:

    • 電波がある距離 r を伝搬するとき、その距離に対応する位相変化を考えます。

位相変化の計算

電波が伝搬する距離 r に対応する位相変化を計算するためのステップは以下の通りです。

  1. 距離を波長で割る:

    • 距離 r を波長 \lambda で割ると、電波が何周期分伝搬したかがわかります。
    \frac{r}{\lambda}
    • これは、電波が伝搬する距離が波長の何倍かを示します。
  2. 位相変化の計算:

    • 1周期の位相変化が 2\pi であることから、全位相変化は次のように計算できます。
    \text{位相変化} = \left( \frac{r}{\lambda} \right) \times 2\pi
  3. 往復距離の考慮:

    • SARでは、電波が送信されて地表で反射し、再び受信機に戻るため、実際の距離は往復距離になります。
    \text{往復距離} = 2r
    • したがって、位相変化も2倍になります。
    \text{位相変化} = \left( \frac{2r}{\lambda} \right) \times 2\pi = \frac{4\pi r}{\lambda}

式 (12.1) の位相の表現

改めて下記式 (12.1) では、SAR信号の第1画像の観測信号 S_1 が次のように表されている、指数部分が位相を表している。Aが振幅を表しているというのが理解できました。

S_1 = A e^{-\frac{4\pi}{\lambda} r_r j}

SAR信号の掛け算、その意味

SAR信号の掛け算は次のように行います。ここでは、第2信号 S_2 の共役複素数 S_2^* を使います。

S_1 S_2^* = A e^{-\frac{4\pi}{\lambda} r_r j} \cdot A e^{\frac{4\pi}{\lambda} r_s j}

共役複素数の使用理由

共役複素数 S_2^* を使う理由は、位相の差を計算するためです。具体的には、共役複素数を使うことで、位相が引き算の形になります。

  • 共役複素数は虚数部分の符号を反転させたものです。

    S_2^* = A e^{\frac{4\pi}{\lambda} r_s j}

位相差が抽出される理由

掛け算を展開すると、次のようになります。

S_1 S_2^* = A^2 e^{-\frac{4\pi}{\lambda} r_r j} \cdot e^{\frac{4\pi}{\lambda} r_s j}

指数の性質を使って整理すると、

S_1 S_2^* = A^2 e^{-\frac{4\pi}{\lambda} (r_r - r_s) j}

複素数の掛け算では、指数部が足し算されますが、共役複素数を使うと符号が反転するため、次のような位相差が抽出されます。

S_1 S_2^* = A^2 e^{j\phi}

このように、掛け算によって位相差 \phi が浮かび上がります。具体的には、次の式で位相差が抽出されることがわかります。

\phi = -\frac{4\pi}{\lambda} (r_r - r_s)

この位相差 \phi は、2つの信号の観測距離の差に基づいて計算されます。

SAR信号の掛け算

まず、2つの信号 S_1S_2 の積を考えます。ここでは S_2 の共役複素数 S_2^* を掛けています。
(本文式 12.7)

S_1 S_2^* = \sqrt{(ac + bd)^2 + (bc - ad)^2} e^{j\phi}

ここで、

  • a, bS_1 の実部と虚部
  • c, dS_2 の実部と虚部
  • \phi は位相差

1. 複素数の指数形式

式 12.7 は複素数の指数形式で表現されています。複素数の指数形式は、絶対値と位相角を用いて表現されます。

z = re^{j\theta}

ここで、

  • r は複素数の絶対値
  • \theta は位相角

2. 信号の積の実数部分と虚数部分

2つの信号の積を計算するとき、実数部分と虚数部分に分けることが重要です。

  • 実数部分 \Re(S_1 S_2^*):

    \Re(S_1 S_2^*) = ac + bd
  • 虚数部分 \Im(S_1 S_2^*):

    \Im(S_1 S_2^*) = bc - ad

3. 位相角の計算

位相角 \phi は次のように求められます。

\phi = \tan^{-1} \left( \frac{\Im(S_1 S_2^*)}{\Re(S_1 S_2^*)} \right)

これにより、2つの信号の積の位相差が抽出されます。

4. 位相差の物理的意味

SAR信号の位相差は、観測点からの距離の差に対応します。具体的には、次のように表されます。

\phi = -\frac{4\pi}{\lambda} (r_r - r_s)

ここで、

  • \lambda は波長
  • r_rr_s はそれぞれ第1軌道と第2軌道からの距離

位相差の数式解説

下記の図を用いて、位相差の式 12.12 を詳しく説明します。

位相差数式

(本文 式 12.12)

\phi = -\frac{4\pi}{\lambda} (r_r - r_s) = -\frac{4\pi}{\lambda} \left( \frac{B_{\text{perp}} z}{r_r \sin \theta} + dD + B_{\text{para}} \right)

この式は、2つのSARの軌道の違いが位相差 \phi に与える影響を示しています。各変数の意味を理解することで、この式がどのように構成されているかを理解できます。

各変数の意味

  • \lambda: 電波の波長
  • r_r: 第1軌道から観測点までの距離
  • r_s: 第2軌道から観測点までの距離
  • B_{\text{perp}}: 垂直基線長。SARの2つの軌道の間の垂直方向の距離。
  • B_{\text{para}}: 水平基線長。SARの2つの軌道の間の水平方向の距離。
  • \theta: SARの入射角。SARの視線が地面と成す角度。
  • z: 地面からの高さ。
  • dD: 地面の変動量。

図 の解説

図 は、2つのSARの軌道と地面の間の関係を示しています。

  • 第1軌道 (S1): 第1画像を取得するSARの軌道。
  • 第2軌道 (S2): 第2画像を取得するSARの軌道。

これらの軌道間の違いが基線長 (B_{\text{perp}}B_{\text{para}}) として示されています。

視線方向の地表変動量 dD

地表変動量 dD は次のように計算されます。

dD = -dz \cdot \cos \theta + dx \cdot \sin \theta

ここで、

  • dz は垂直方向の変動量
  • dx は水平方向の変動量

これらを組み合わせて、視線方向の変動量 dD を求めます。

式 12.12 の詳細

  1. 位相差の基本式:
\phi = -\frac{4\pi}{\lambda} (r_r - r_s)

これは、2つのSAR軌道からの距離の差に基づいて位相差を計算する基本式です。

  1. 距離の差の展開:
r_r - r_s = \frac{B_{\text{perp}} z}{r_r \sin \theta} + dD + B_{\text{para}}

ここで、

  • \frac{B_{\text{perp}} z}{r_r \sin \theta}: 垂直基線長 B_{\text{perp}} による距離の差
  • dD: 地面の変動による距離の差
  • B_{\text{para}}: 水平基線長 B_{\text{para}} による距離の差

これらを組み合わせて、全体の距離の差を求めます。

  1. 位相差の計算:

最終的に、これらの距離の差を位相差に変換します。

\phi = -\frac{4\pi}{\lambda} \left( \frac{B_{\text{perp}} z}{r_r \sin \theta} + dD + B_{\text{para}} \right)

視線方向の地表変動 dD

上記の式を変形することで、視線方向の地表変動 dD を計算するための式です。
(本文式 12.13)

dD = -\left( \frac{\lambda}{4\pi} \phi + \frac{B_{\text{perp}} z}{r_r \sin \theta} + B_{\text{para}} \right)

この式は、位相差と基線長の情報を組み合わせて、視線方向の地表変動 dD を計算します。具体的には、次の3つの項から構成されています。

  1. 位相差項: \frac{\lambda}{4\pi} \phi
  2. 垂直基線長項: \frac{B_{\text{perp}} z}{r_r \sin \theta}
  3. 水平基線長項: B_{\text{para}}

位相差の周期性に関する考察

上記12.13 で位相差に関する項 \frac{\lambda}{4\pi} \phi のみで考えると、
位相差 \phi が複数の2 \pi 周期を持つため、それに対応する地表変動 dD については最大で下記の変動をもつことになる。
(式 12.14)

dD = -\frac{\lambda}{2}

この式は、位相差が2 \pi の周期を持つため、地表変動も波長 \lambda の半分の周期で変動することを示しています。位相差 \phi が2 \pi に対応する変動量 dD を示します。

なぜ位相差が重要なのか

位相差は、SAR技術で地表の変動を検出するために重要な情報です。位相差を解析することで、以下のことがわかります。

  1. 地表の高さ変動:

    • 2つの異なる時点での観測距離の差に基づいて、地表の高さの変動を検出します。
  2. 移動量の計測:

    • 位相差から、地表の移動量(水平移動や垂直移動)を計測します。
  3. 高精度な地表マッピング:

    • 位相差情報をもとに、地表の高精度なマッピングを行います。

アンラップ

位相のアンラップとは?

もともと位相は -\pi から \pi のあいだの値としてでてくるように計算されます。それを積算した値に変換する処理がアンラップというのですね。 地殻の変動量についても位相の差が、1周、2周するについて 振幅分変化は積算されていくのがSAR干渉の縞模様ということですね!

回答12.1

回答自体はシンプルです
Xバンド,Cバンド,LバンドのSARがある。波長 3cm, 5.6cm, 24cm に対応する地殻変動量の範囲は
-\frac{\lambda}{4} から \frac{\lambda}{4}  なので
それぞれ -0.75〜0.75cm, -1.4〜1.4cm, -6〜6cm ですね。

問題12.2

r_r=700km, 0=45度,DEMの精度が10mの場合について,B_{perp} = 1 km, 0.1 km, 0.01 km, 0.001kmの4ケースについて,dDに及ぼす誤差成分を計算しなさい。

回答12.2

本文式12.13

dD = -\left( \frac{\lambda}{4\pi} \phi + \frac{B_{\text{perp}} z}{r_r \sin \theta} + B_{\text{para}} \right)

より、 与えられた条件が影響するのは 垂直基線長項: \frac{B_{\text{perp}} z}{r_r \sin \theta} のみです。

与えられた条件を代入します。

\sigma_{dD} = \left( \frac{B_{\text{perp}}}{700 \times 10^3 \times 0.707} \right) \times 10

B_{\text{perp}} : 誤差
1\times 10^3m : 0.0202m = 2.02 cm
1\times 10^2m : 2.02 mm
1\times 10^1m : 0.202 mm
1m : 0.0202 mm

問題 12.3

2つの人工衛星が地上の点をのぞき込む角度差により,周波数シフトが発生し,これにより2衛星データの干渉すべき帯域幅が狭められる。今,高度628km,衛星1の地表点への入射角を35度、衛星2はその内側1km を平行に飛行していた場合の,周波数シフト量はいくらになるか?
送信周波数は1.27 GHzとする。

周波数シフトと臨界垂直基線長(B_perp)

1. 周波数シフトとは?

周波数シフトとは、2つのSAR(合成開口レーダー)衛星が同じ地上の点を観測するとき、2つの衛星が異なる角度から地上を「のぞき込む」ことによって起こる現象。
このとき、信号の波長(つまり周波数)が角度の違いによってわずかに変化します。この「わずかな変化」が周波数シフトです。


2. なぜ周波数シフトが起こるの?

SAR1とSAR2が異なる角度(入射角)で地上の同じ点を観測すると、信号の「波長の見かけの長さ」が異なります。

  • 上の図で、SAR1とSAR2が地上で錯乱物にあたり作られる波長が異なる

波長のシフト量は以下の式で計算されます:
$$
\Delta \lambda = \frac{\lambda}{\tan \theta} \Delta \theta
$$
ここで、

  • \lambda:信号の波長
  • \theta:SAR1の入射角
  • \Delta \theta:SAR1とSAR2の入射角の差

周波数のシフト量に変換するには、次の関係を使います:
$$
\Delta f = \frac{f}{\tan \theta} \Delta \theta
$$

  • f:SAR信号の周波数(送信周波数)

3. 垂直基線長(B_perp)と干渉性の関係

SAR1とSAR2の間の「基線長(B_perp)」が大きくなると、2つのSAR衛星が観測する角度の差(\Delta \theta)が広がります。この結果:

  1. 周波数シフト量が大きくなる

    • 角度差が広がることで、観測波長のシフト量が大きくなります。
    • これにより、信号が占める帯域幅が広がります。
  2. 共通帯域幅が狭くなる

    • SAR1とSAR2が使う周波数帯域が一致する部分(共通帯域幅)が狭くなります。
    • 共通帯域幅が狭いと、干渉画像が正確に作れなくなる可能性があります。

4. 臨界基線長(Critical Baseline)

「臨界基線長」とは、SAR1とSAR2の基線長が長くなりすぎて共通帯域幅がなくなる直前の状態を指します。

  • 基線長が増えるほど干渉性が低下(0に近づく)し、最終的に共通帯域幅がなくなります。
  • 臨界基線長を超えると、SAR信号が重ならなくなり、干渉画像を生成できません。

臨界基線長の式:

B_{\text{crit}} = \frac{\lambda H}{\Delta f \cdot \tan \theta}

ここで:

  • \lambda:波長
  • H:衛星の高度
  • \Delta f:帯域幅
  • \theta:入射角

観察点:
以下の点でX-band, L-band それぞれの特性がでてくることがわかります

  • \Delta f(帯域幅)が狭くなると、臨界基線長 B_{\text{crit}} が小さくなる。
    • X-bandのほうが帯域幅は広げやすい
  • 波長がながいと臨界基線長は大きくなる
    • L-band 24 cm, X-band 3cm なので圧倒的にL-bandのほうが有利
  • 高度は高いほうが臨界基線長は大きくなる

SARの帯域幅とX-band/L-bandの違い

X-bandとL-bandでは、帯域幅が以下のように異なります。

周波数帯 周波数範囲 一般的な帯域幅 用途
X-band 約8~12GHz 100MHz~400MHz 高解像度画像、都市部や災害地域の観測
L-band 約1~2GHz 20MHz~80MHz 地殻変動、植生観測、長期間観測
  • X-band: 周波数が高く、帯域幅が広い。このため、細かい地形や構造物を高精細に観測可能。
  • L-band: 周波数が低く、帯域幅が狭い。しかし、森林や植生、地殻変動などを観測する際に信号減衰が少なく有利。

** 帯域幅が狭い場合の制約**

帯域幅が狭いと、SARの干渉性は保ちやすい一方で、以下のような制約が発生します:

  1. 基線長が短くなる:

    • SAR1とSAR2の間の垂直基線長 B_{\text{perp}} を短く保たないと、信号の重なり(共通帯域幅)が失われます。
    • 基線長が短いと、観測の「立体視効果」が弱くなり、地形起伏の情報(DEM生成など)の精度が低下します。
  2. 撮像条件が狭くなる:

    • SAR衛星の配置や飛行軌道に対して、厳密な条件(基線長の制約)が必要。
    • 広いエリアをカバーする撮像計画において柔軟性が失われます。

** 帯域幅を広げた場合のトレードオフ**

一方、帯域幅を広げるとどうなるでしょうか?

  • 利点:

    • 臨界基線長が長くなるため、SAR1とSAR2の位置関係に柔軟性が生まれます。
    • 広い立体視が可能になり、地形の起伏情報(高度分解能)が向上します。
  • 欠点:

    • 周波数シフトが大きくなり、干渉性が低下しやすい。
    • 雑音や位相ノイズの影響を受けやすく、精度が低下する可能性。

SARシステムの運用では、以下のように目的に応じたトレードオフが必要です:

  1. 精密な干渉画像が必要な場合(例:地殻変動の観測):

    • 帯域幅を狭くして、干渉性を優先。
    • 垂直基線長を短く抑え、柔軟性を犠牲にします。
  2. 広範囲の地形データが必要な場合(例:DEM生成):

    • 帯域幅を広げて、基線長の自由度を確保。
    • 干渉性の低下を補うため、複数回の観測やデータ補完を行う。

問題 12.3 回答

問題設定

  • 高度 H = 628 km
  • 衛星1の地表点への入射角 \theta = 35
  • 衛星2は衛星1の内側1kmを平行に飛行
  • 送信周波数 f = 1.27 GHz

問題

  • 周波数シフト量 \Delta f を求める

図 12.5 の式を用いた計算手順

  1. 入射角の変化 \Delta \theta の計算

    \Delta \theta は、衛星1と衛星2の入射角の違いを表します。内側1kmの飛行による入射角の変化を計算します。

    \Delta \theta = \theta - \tan^{-1} \left( \frac{H \tan \theta - 1.0 \, \text{km}}{H} \right)

    ここで、\theta = 35 度、H = 628 km を代入します。

    \Delta \theta = 35 - \tan^{-1} \left( \frac{628 \times \tan(35^\circ) - 1.0}{628} \right)
    \tan(35^\circ) \approx 0.7002

    よって、

    \Delta \theta = 35 - \tan^{-1} \left( \frac{628 \times 0.7002 - 1.0}{628} \right)
    \Delta \theta \approx 35 - \tan^{-1} \left( 0.6989 \right)
    \Delta \theta \approx 35 - 34.9387356 \approx 0.0612644 \, \text{度}

    角度をラジアンに変換します。

    0.0612644 \, \text{度} \approx 0.0612644 \times \frac{\pi}{180} \approx 0.001068 \, \text{ラジアン}
  2. 周波数シフト \Delta f の計算

    \Delta f = \frac{f}{\tan \theta} \Delta \theta

    ここで、\tan(35^\circ) \approx 0.7002 ですので、

    \Delta f \approx \frac{1.27 \times 10^9}{0.7002} \times 0.001068 \approx 1.93 \times 10^6 \, \text{Hz} = 1.93 \, \text{MHz}

衛星1と衛星2の間の周波数シフト量は約 1.93 \, \text{MHz} です。

まとめ

遠くの衛星2つの波の位相の差分を観測する。
その観測を実現するための条件、共通帯域幅 臨界基線長、 B_perp。
これらは理論的には数十年前に確立されていたもの。
遅ればせながら、これらの言葉の意味、深みを味あわせていただきました。

Discussion