はじめに

SAR業界から遠のいてしまいましたが、博士課程まで勉強させてもらっていたので、いろいろと共有したいことはあります。とにかくやっている人が少ない分野なのですが、ここ数年の宇宙分野の機運の高まりによって、SAR業界も日の目を見ることが近そうな状況です。

それに伴って勉強する人も増えている（という勝手な希望）と思うので、そうした方々に少しでも報いるべく、ここでは、合成開口レーダの後方散乱について、知っていることをつらつらと書いていこうと思います。

後方散乱係数って何よ

モノスタティックレーダシステム（送受信アンテナが同一位置にある（とみなせる））を考えます。

アンテナが電波を送信した際、何らかの観測対象（以下ターゲット）で散乱したとき、送信アンテナ側に返ってくる散乱波を後方散乱と言います。逆に電波の送信方向にいっちゃうのを前方散乱と言います。レーダの受信信号なのに後方散乱という呼び方に私は最初、違和感満点でしたが、電波自身の目線で考えるから前方後方となるわけですね。

これを衛星搭載SAR（航空機でもいいのだが）の文脈でとらえますと、斜め下方向に照射した送信マイクロ波は、数100㎞を伝播し地面にぶつかります。光の反射と同じで大部分の電波は鉛直方向にのみ方向を変えてそのまま再び空中へ散乱していきます。しかし、地面には凹凸もありますし、そこには建物や木々などの地物もあるので、一部の電波は後方散乱します。モノスタティックのシステムでは、この後方散乱を受信することでターゲットまでの距離（さらには物性や速度など）を推定することができるというわけです。

さて、後方散乱の定性的な理解についてはこの通りですが、これはサイエンスの話なので、何らかの定量化が必要です。そこで、後方散乱係数、後方散乱断面積といったワードが必要になります。これらが何なのかを先に言うと、後方散乱係数は、散乱係数の後方散乱成分、後方散乱断面積は散乱断面積の後方散乱成分ということになります。

・・・めちゃくちゃ意味のない説明ですね。それでは散乱係数や、散乱断面積って何よ？って話になりますね。

散乱断面積って？

例えば、電波を送信し何らかのターゲットからの後方散乱を受信アンテナで受信して、その強度として1という値を得たとします。これをそのターゲットの固有の散乱の強さを表す値として考えていいでしょうか。

そんなことはないですね。この値には、そもそもどのくらいの電界の電波を送信したのか、とか受信アンテナのゲインとか、ターゲットまでの距離なんかの情報がコミコミになっています。つまり、ターゲットの物性とは無関係の情報によって汚れてしまっています。このような値ですと、例えば違うレーダシステムで計測した場合に、同じターゲットであっても異なる値が計測されてしまい、比較が不可能となります。

そこで後方散乱断面積の出番です。後方散乱断面積は、計測値から不要な情報を取り除くというか補正することで、シンプルにターゲットの物性としての後方散乱の強さを表します。後方散乱断面積を使えば、異なるセンサ間であっても、値のフェアな比較ができます（入射角とかが違うとそれも考えなくてはいけませんが、それは後述）

具体的に式で示したいところですが、いくらでもわかりやすく書いてくださっている先達がいらっしゃるのでここではリンクを載せておきます。次の資料の「有効反射断面積」となっているパラメータが、後方散乱の文脈で言うところの後方散乱断面積となります。 physics.thick.jp

散乱係数って？

それでは散乱係数って何でしょうか。散乱断面積では散乱におけるエネルギーに主眼が置かれていました。散乱係数はこれの振幅と位相での表現だと思ってもらって大丈夫です。つまり、空気とターゲットとの境界面において、入射電界と散乱電界の比率で表します。振幅と位相を表すため、通常複素数を使って $s=Ae^{j\phi}$ のように表します。

規格化散乱断面積？

話はこれで終われば素直なのですが、例えば国産航空機SARであるPiSAR-L2の次のリンクを見ますと、「規格化後方散乱断面積」という単語が現れます（正規化、とする訳もあります。英語だとNormalized radar cross section）。

https://www.eorc.jaxa.jp/ALOS/Pi-SAR-L2/cal_valt.html

要約すると、生データに対して補正を施すことで規格化後方散乱断面積に変換できます、という内容です。補正については、また後程触れるとして、規格化、とはいったいなんのことでしょう。何らかの基準値で規格化するのだということがわかりますが、いったい何の値を使うのでしょうか。

ここからはうろ覚えでこんな感じじゃなかったかな程度の知識なので、もし誤っていることが分かった方にはぜひ訂正のコメントを頂ければと思います。

規格化は導体球（金属球）の散乱断面積を基準として行います。

次の新潟大学山口先生の資料が非常にわかりやすいです。

http://www.wave.ie.niigata-u.ac.jp/yamaguchi/education/information/sensing_system.pdf

5ページの1.3レーダ散乱断面積（RCS)から少し読み進めますと、図1.7で正規化後方散乱断面積を求めていますが、これも金属球の散乱断面積を基準値としています。理論値がわかっている、波長依存性がない、散乱が等方的、偏波依存性がないということで基準として適しているのだと思います。（１点注意事項ですが、波長依存性がない、と書きましたが、 $\pi a^ 2$ という値は、波長に対して十分大きい径の球の散乱断面積なので、本当は波長依存性はあります。あるのですが、 $a$ が定数で波長を小さくしていった場合、値が収束します。この性質はほかのシンプルなターゲットでは見られません。なので波長依存性がない、と書きました。）

なるほど、金属球の散乱断面積の理論値で規格化するのか!

・・・、径の $a$ はどう決めるの？と思いますよね。ここがうろ覚えなのですが、取得した衛星データの解像度（1画素の面積）と等価になるように $a$ を決めるのではなかったでしょうか。たぶん。つまりデータの解像度と同じ投影面積を持った金属球が存在していたとすると、その規格化後方散乱断面積は1となります。規格化の説明のために金属球を引き合いに出しましたが、単に解像度で割るでもいいです。規格化した後の値の大きさの具体的なイメージを持つために、金属球を考えるとイメージが少し湧くかなと思います。

衛星搭載SAR生データの補正（ラジオメトリック補正）

さて、ここまでで理論、というか座学はおしまいです。この説では実際に生データを読みだしてきて、後方散乱係数、後方散乱断面積に補正してみたいと思います。これはラジオメトリック補正と呼ばれる処理です。

今回はこのデータを使って、確認していこうと思います。

JAXA謹製のサンプルデータ

ftp://ftp.eorc.jaxa.jp/pub/ALOS-2/1501sample/020_tokyo/0000008660_001001_ALOS2014410740-140829.zip

使う言語はMATLABです。生データの読み方については、手前味噌で恐縮ですが、この記事にて記しています。 qiita.com

生データの確認

まずは読み込んだ生のI,Qを値から振幅を計算し、表示してみます。

% imageDataに複素数データを格納したとする
img = abs(imageData(1:2:end,1:2:end)); % 絶対値の計算（重いので間引いた）
img = imresize(img, 0.25); % 重いのでさらに縮小
imshow(img, []); % 表示

実行してみるとわかるかと思いますが、真っ暗になります。これは、ダイナミックレンジが広く、かつ、値の大きな建物などの地物は全画素の中で少ないからです。表示する値の範囲の調整が必要です。

imtool(img);

でヒストグラムを見つつ、白黒のレンジをGUIで変えることができます。だいたい下図のようにしきい値を決めます。

f:id:deokuradar:20200130165238j:plain — imtoolでのグレースケールのしきい値決めの状況

出てくる画像はこんな感じになります。

f:id:deokuradar:20200130165345j:plain — 表示されるSAR画像

ちなみに最小値は

min(img(:))

ans =
  -4.6345e+05

あれ？！マイナス？おかしくない？これはどうやらimresizeがデフォルトでcubicで内挿を加えるのが原因みたいです。bilinearに指定し直すと

img = imresize(img, 0.25, 'Method', 'bilinear');
min(img(:))

ans =
   3.8519e+04

ということで無事に正の値が出てきました。

ということで最小値でも38,000という値になっていて、規格化されていないことは明白な状況です。

ラジオメトリック補正式の確認

ラジオメトリック補正式については、PALSAR-2の公式ドキュメントに記載があります。

PALSAR-2プロダクトフォーマット説明書

こちらの73ページに次のように記されています。

校正係数（CF）

レベル1.1 ：σ0= 10*log10<I2+Q2> + CF -32.0

レベル1.5/3.1：σ0= 10*log10 + CF

本式は、該当するピクセルの後方散乱係数がアンサンブル平均<>で求まること、つまり、求めたい点のまわりについての平均処理で求まることを表す。ここで、I, Qはレベル1.1の、DNはレベル1.5/3.1のピクセル値である。

ということです。このあたり後方散乱断面積と後方散乱係数、規格化後方散乱断面積の用語がごっちゃになってるんじゃないのと思います。実際私も後方散乱断面積の意味で後方散乱係数って使ったりしてましたし、まぁ大きく違うものを指しているわけでもないし、そこまで気にしなくてもいいのかもしれません。この記事では説明の都合上、上述のように係数は振幅と位相を含む複素数の形で表される波ベースのもの、断面積はスカラで表されるエネルギーベースのものを考えています。

CFの値自体はまさにこの説明が書いてあるフィールドを読めばいいです。

fid = fopen('LED-ALOS2014410740-140829-UBSL1.1__A', 'rb', 'ieee-be');
offset = 720+4096+4680+16384;
fseek(fid, offset+20, 'bof');
char(fread(fid, [1,16], 'char'))
fclose(fid);

ans =
    '     -83.0000000'

ラジオメトリック補正

それではこの式を使ってラジオメトリック補正をしてみます。簡単ですね。

規格化後方散乱断面積のラジオメトリック補正

sigma0 = 10*log10(img.^2) - 83 - 32.0;

補正後の値

imtoolで補正後の値を見てみましょう。こちらがヒストグラムを見ながら適当に決めたしきい値。

f:id:deokuradar:20200131071830j:plain — imtoolでの後方散乱断面積のしきい値決めの状況

こちらがその後方散乱断面積の画像

f:id:deokuradar:20200131071855j:plain — 表示される後方散乱断面積画像

ちなみにこれは「規格化」後方散乱断面積なので単位は無次元になります。もともとは散乱断面積なので[m2]です。

値に注目すると大体-20から+5くらいになっていますね。経験上ですが、このレンジになっていればまぁまず正しくラジオメトリック補正できていると言えます。

少し細かく見ていきましょう。

水面上

f:id:deokuradar:20200131072722j:plain — 水面上の規格化後方散乱断面積

値は-16くらいですね。鏡面反射と呼びますが、水面はなめらかで後方散乱させる地物もないので、前方散乱がかなり強くなり、後方散乱は微弱です。入射角や波の立ち方に大きく依存するのですが、30度後半から50度くらいのレンジでは暗く映っていて-15以下くらいが経験上は適正値です。

木

f:id:deokuradar:20200131072955j:plain — 木の規格化後方散乱断面積

値は-7くらいですね。枝葉で構成されるや樹冠部分でのランダムな散乱（体積散乱）が生じています。今回のLバンドだと、それらを透過して、地面からも散乱があります（表面散乱）ので、これらのミックスによりまぁまぁの後方散乱があります。こちらも波長とか樹種とか地面状況に依存してきますが、おおむね-5から-12とかです。あくまで経験上の感覚値ですが。

建物

f:id:deokuradar:20200131073134j:plain — 建物の規格化後方散乱断面積

10くらいですね。地面と建物の二面構造が強い後方散乱を生み出します。なので値も高くなり、5以上が感覚値です。建物の向きによってかなり違います。

後方散乱係数のラジオメトリック補正は、I,Qの段階でのスケーリングしてあげて、-83-32の項がなくてもいいようにしてあげればいいので

$\displaystyle{ \begin{aligned} 10\log_{10}((kI)^2+(kQ)^2) &= 10\log_{10}(I^2+Q^2) -83 - 32 \\ \end{aligned} }$

これを $k$ について解けばいいです。数式打つのがしんどいのでここまでにします。結果についても振幅が変わるだけなので特に表示まではしません。

TerraSAR-Xの場合は？

以上がPALSAR-2の場合でした。TerraSARシリーズの場合はどうでしょうか。図まで出そうかと思いましたが、力尽きてしまったので、ラジオメトリック補正式のリンクを載せておきます。ちなみにラジオメトリック校正は、radiometric calibrationですので、ググるときはこれで。

https://spacedata.copernicus.eu/documents/12833/14537/TerraSAR-X_RadiometricCalculations

基本的に一緒です。校正係数もメタデータに記載がありますので、そこを読んでくればOK。sigma naughtとか beta naughtとかは入射角とか地形の影響とかを考慮しているものになりますが、今回は割愛します。

Sentinel-1A/1Bの場合は？

フリーで使えるのがうれしいセンチネルシリーズです。同じくリンクをば。

https://sentinel.esa.int/web/sentinel/radiometric-calibration-of-level-1-products

冷静に考えるとセンチネルのラジオメトリック補正はやったことがないので校正係数とかよくわかりません。ググった感じだとラジオメトリック補正されたプロダクトになっているのかな？？また詳細がわかったら追記するかもしれません。

注意点

後方散乱係数の計算についていくつか注意点を書きます。

偏波SARにおいて

偏波間の関係が問題になるケース、例えば偏波間の位相差とか相関係数とかの計算においては校正前のままで特に問題はありません。位相差や相関係数の計算は基本複素数同士の位相差のばらつきを計算するのであって、絶対値には左右されません。次の式からもその点は明らかです。 $s _ A,s _ B$ は任意偏波の散乱係数、 $\sum$ 記号は空間平均を表します。それぞれの偏波成分を $a,b$ でスケーリングしてあげても、位相差、相関係数は変わらない、という式となります。

偏波間の位相差

$\displaystyle{ \begin{aligned} \angle{\theta_ {AB}} &= \angle \sum {s_ {A}{s_B}^*} \\ &= \angle ab\sum {s_ {A}{s_B}^*} \\ &= \angle \sum {(as_ {A})({bs_B}^*)} \end{aligned} }$

偏波相関

$\displaystyle{ \begin{aligned} r_{AB} &= \frac{|\sum {s_ {A}{s_B}^*}|}{\sqrt{\sum {s_ {A}{s_A}^*}}\sqrt{\sum {s_ {B}{s_B}^*}}} \\ &= \frac{ab|\sum {s_ {A}{s_B}^*}|}{ab\sqrt{\sum {s_ {A}{s_A}^*}}\sqrt{\sum {s_ {B}{s_B}^*}}} \\ &= \frac{|\sum {(as_ {A})({bs_B}^*)}|}{\sqrt{\sum {(as_ {A})({as_A}^*)}}\sqrt{\sum {(bs_ {B})({bs_B}^*)}}} \\ \end{aligned} }$

ただし、散乱電力分解など絶対値が問題となるケースでは、やはり後方散乱係数への変換が必要です。

干渉SARにおいて

干渉SARでは、2時期間の計測値の位相差やコヒーレンスが問題となります。つまり、上述の偏波間の位相差や相関係数と数式で言えば同じであり、やはり後方散乱係数への変換の有無は結果に影響しません。

入射角などの観測条件において

実はSARが使いづらい元凶ともいえる問題がこれなのですが、入射角や計測に利用する波長が異なる場合、後方散乱係数の単純比較は意味がありません。一般的に地物には散乱の入射角依存性がありますし、波長依存性があります。

ですので、異なる入射角や計測波長で同一地物を観測して、異なる後方散乱係数を得たとして、それが計測間での地物の変化によるものなのか、単に計測条件の違いによるものなのかの判断はデータのみでは不可能です。もちろん、グランドトゥルースとして事前もしくは事後に地物の状況がわかっていれば判断はつきます。もし、地物の物性（形状や誘電率）が変わっていなくて、計測条件によって後方散乱係数が変化している場合は、その変化分は計測条件の違い、ということになりますし、地物の物性が変化している場合は、得られた後方散乱係数の違い（もしくは同じであっても）は、地物の変化による成分と計測条件の違いによる成分の和となっています。

ややこしいですよね？なので基本的に時期間での後方散乱係数の違いを観察したい場合、計測波長や入射角などがそろった同一衛星同一軌道のデータで比較します。

ってかこんなこと知らなくても専用ソフトで読みだせば自動でラジオメトリック補正してくれるし

その通りです。ただ使う分にはそれでいいかと思います。しかし、ソフトが裏側で何をしているかについての理解はいろんな場面であった方がいいと思います。少なくともSARを専門で行う技術者にとっては必須知識と思います。

おわりに

SARの後方散乱係数について、知っていることを記録しました。読んだ方の何かの参考になれば幸いと思いつつ、間違ってたらすみません。ここまで書いておいてなんですが、SARの真骨頂は位相情報なのでそのあたりもいずれ書きたいと思います。

ひろい集め学習

レーダ画像解析技術者が技術全般についてやったこと学んだことを書いておく

SAR（合成開口レーダ）の後方散乱について