SARの勉強５：データを見てみる
アニメーションにするといいらしい
前処理について理解が少し進んだので、データを見てみる。
が、何をみてるかよくわからない。教材動画ではアニメーションにするとわかりやすいといっているので、やってみる。 コードはスライドが公開されているので、そこのリンクからコピペ＋場所やら範囲やらを変更したらできる。
スライド docs.google.com
先生のコードはこちら https://code.earthengine.google.com/3a7843212acec68df670188fdb29e814
アニメーション
岩木山
網走
尾瀬
那覇空港
山梨県全体
信楽周辺
八郎潟
一番面白かったのがこれ。田んぼに水が入るのが見て取れる。 八郎潟をもうちょっと見ることにした（次回）。
まとめ
アニメーションにしたらわかる、ということだったけど、わかるようなわからないような。
一つの画像をみるよりはわかる。ただ、そこがどういうところなのかを知らないとわからない。
GIFは簡単にできた。わかる/わからないは置いておいて、とりあえず作るといい。
SARの勉強４：Terrain Flattening その２
DEMの違い：SRTMとALOS
Google Earth Engineに入っている全球で使えるDEMは2種類。どちらも解像度30mだけど、SRTMはおよそ90m解像度のものを内挿してして30m1に、ALOSのは2.5m？だかの高精細なものをアップスケールして30mにしているという違いがある。この違いがTerrain Flatteningしたときはどう変わるか見てみる。
生駒
元のデータは2020年8月のDESCENDINGのVVデータ。違いは微妙だな。ALOSの方が落ち着いて見えるかな？程度。
黒部ダム
もっと急な場所を見てみたいので黒部ダム。元のデータは2020年8月のDESCENDINGのVVデータ。
薄緑色に抜けているのは山の陰になって値が取得できなかった場所。 これはALOSの勝ち。SRTMは細かい谷の陰にできるシワシワ感があるのに対し、ALOSは全体になじんだトーン。
というわけで、ALOSの方がよさそう。 一つ気になったのが、ALOSのDSMは都会の大きなビルがわかるくらい細かいから、そういう場所（ビルの側面とか）の斜度をもとに補正すると極端な値になったりしないのかな？ということ。ただ新宿でやってみてみたけど、（見た目では）変な感じになってなかったから心配しなくていいのかも。
すごい急なところはどうなる？
黒部をみて思った。なんだかんだ言って、補正にも限界はある。谷間はしんどいし、尾根の向こう側もデータが取りづらい。これはもうどうしようもないんだろうなとわかりました。
で、じゃあ、何度くらいのところまでなら補正でしのげるのか？というのが気になる。 その周辺の地形と斜面方位も関係するはずから、どうやるべと思ったんですが、とりあえずグラフをいろいろ書いてみます。 まずはお手本？というか、超平らな場所はどうなんだろうと
傾斜は違うけど、広範囲に土地被覆に変化がなさそうな場所ということで、原生的な？ところを選んでみる。ざっとみたらスペックルノイズが気になったから、２か月間の平均をとっているけど、それでもノイズは残っていそう（特にVV）。グラフをRで書いてるけど、y軸の最大を適当に制限しているのでそれ以上の外れ値は無視しています。
白神山地
奄美大島南部
御蔵島
印象とまとめ
特にVVはノイズが残っているみたいで点のばらつきが大きい。とはいえ、VVもVHも全体をとおしてみると30度くらいを境に、森林でその値はないんじゃない？という高い値が目立ってくる。というわけで30度を超える斜面は止めた方がよさそう。 で、こういう異常値の点の色は揃っていてAESCENDINGで青色、DESCENDINGで黄色から緑なので、衛星の進行方向によって値が厳しくなる斜面方位がある（＝陰になるとか）というのが推測できる。
一方20度くらいでもそれなりに高い値は出ているから、何度だったら大丈夫というのはちょっと言えない感じ。例えば傾斜が小さくても陰になっているとかあるだろうし。
まとめ
傾斜は小さいほうが良い。30度以上は止めること推奨。 何度だったら大丈夫かはわからない、（できたら）現地データとの比較で決めるとか。
データ取得のコードはこちら：https://code.earthengine.google.com/f6c94988e988cc1b90c7c382f2809970
SARの勉強４：Terrain Flattening その１
Radiometric Terrain Flatteningがよくわかってない
さて、次の補正ですが、何をどうする補正かよくわかっていません。そもそも日本語で何というのか、わかってません。超おおまかな理解としては、斜面だとマイクロ波の入射角が変わるから真の（水平面だった時の）値より大きくなったり小さくなったりする。それを衛星の角度のデータとDEMから補正しようということだと思います。
わかっていないから、まずはデータを見てみようと。そこで、以前からお世話になっているGoogle Earth EngineでSentinel-1の前処理をするスクリプトに登場いただき（下のリンク）、補正前後の様子を見てみる。
処理前後の比較
見た目
場所は生駒の遊園地のあたりとして画像を見てみる。
処理前は画像左半分は白っぽく、右半分は黒っぽいのに対し、処理後は画像左右の明暗の差は減って平坦な感じに見える。 この画像では場所の様子がわかりづらいのでGEEで見たほうが納得しやすい↓
比較用作図のコード：https://code.earthengine.google.com/83d0a3bd6a4b923f3e03320ebd70a230
散布図を描く
散布図を書いて左右の明暗が減ったか確認してみる。 同じく生駒で、2000年8月上旬に取得された画像を題材にみてみる。衛星の進行方向ASCENDINGとDESCENDINGのそれぞれをみる。データ取得まではGEEでやったけど、散布図はRで書いた。数点大きく外れる値があったけど散布図のy軸は適当に上限を決めて表示してます。
図から見て取れるように、処理前はAESCENDINGとDESCENDINGのVVの値が東経135.675と135.680の中間くらいで交差している。処理後は全体に混ざり合って交差は見えない。南北に走る生駒山地の主稜線が135.6775くらいにあって、そこを境に衛星から近い・遠いがでる（入射角がガラッと変わるから）せいで、VVの値が変わる。だけど、処理することでその影響が緩和されたということだろう。
ASCENDINGとDESCENDINGそれぞれで見てみるとこんな感じ。
VHも似た感じだった（グラフ以外省略）。
こうなるっていうのが少しわかった（続く）。
SARの勉強３：スペックル その３
複数画像ImageCollectionのmean()とmedian()はどうだろう？
SARの勉強シリーズ１回目でも触れたけど、論文によっては複数画像の平均値だったり、一筆の畑の平均値だったりを使うことでスペックルノイズに対処しているのがある。どのくらいノイズ除去効果があって、エッジは残せるのか？
中心極限定理だろ
スペックル１回目で出てきた中心極限定理。分布の形を問わず、平均μ、標準偏差σの集団からN個サンプリングするとその平均はμ、標準偏差σ/√Nになる。
となると、一定期間でピクセルの散乱強度係数が変わらないと仮定したとき、その期間に得られたN枚の画像の平均から得られるピクセルの値は平均μ、標準偏差はσ/√Nとなると考えられる。つまり、データのバラつき（=スペックルノイズに起因）は標準偏差で表すと1/√Nで減っていくはず。
空間で平均をとっても同じように使える。畑や田んぼ、同じような林など散乱特性が同一と扱っても問題ないような一定面積の場合でも平均の分布を取れば、平均は変わらず、バラつきは小さくなるはず。
というわけで、単純に平均をとるだけでフィルタリング効果的なの得られるんじゃない？と思われるわけです。
やってみる
処理に使うのは3枚（前後1枚ずつ：±12日）と8枚（だいたい前後1か月ずつ）として、平均と中央値を取ってみた。 まずは画像から。前回と場所は同じ。
ターミナル付近
東側のゴルフ場
滑走路南端
感想
いい感じでノイズは減っているけど細部も残っている。フィルタリングよりも断然細部がくっきり。
3より8の方が綺麗な画像。
平均と中央値はあまりかわらない？
というわけで、全体としていいんじゃないかと思いました。
平均とか分散とか
CVの棒グラフ
比較のため1枚を使ったフィルタリングの結果も載せた。３枚の平均でもImproved Sigma Leeフィルタと同程度くらいはノイズが下がってる。8枚だともっといい感じ。
ヒストグラム
平均3枚
平均8枚
中央値3枚
中央値8枚
中心極限定理はあってる？
今日のまとめ
平均値でもフィルタと同じようなノイズ削減効果が期待できる。
舗装された滑走路は平均をとるとノイズが増えた。ツルツルしたところはそうなるんだと推測。
3枚の平均でもフィルタと同等（詳しくはフィルタによるけど）までノイズが減る。枚数を増やせば増やすほどノイズは減る。
平均のフィルタに対する利点は２点：一つは細部がつぶれずに残ること。もう一つは処理時間が短いこと。フィルタはGoogle Earth EngineでのAssetsへの書き出しに２時間程度かかったけど、平均だと2～3分で終わった。
フィルタの平均に対する利点は１点、ただし決定的：１枚の画像だけで処理できること。もちろんSentinel-1は無料だし、GEEでサクサク処理できるから何枚使ってもいいんだけど、時間の経過で地表の状態が変わったら使えない。農地のように短時間で状態がガラッと変わる場所は、時間的な平均をとる処理は正しくない。
空間的な平均も時間的な平均と同じようにノイズ削減効果が期待できそう（やってないけど）。
平均と中央値では平均がいいと思う。似たような結果だったけど中心極限定理を踏まえてとなると平均の方が妥当に思う。
フィルタをかけて平均をとったり、平均をとってからフィルタをかけたりすると更にノイズは減るのだろうとは思う。まあ、なんだかんだいってそこまで必要な場合はそんなにないのかもしれない。自分の使いたいデータ分析の方法でどの程度までノイズを許容するかという話だろうから。
平均をとる枚数を増やせば増やすほどノイズは減るけど、減り方はだんだん鈍くなる（1/√枚数）。また、長期間の情報をつかうことになるから、地表の状態が同じという条件を満たすのが難しくなる。たとえば「夏」だと森林かわらない！とか自信を持って言える期間に限定するのがいいと思う。
というわけで、スペックル編はいったんここまで。 勉強になりました。
SARの勉強３：スペックル その２
フィルタリングの効果を見てみたい
ちょっと前に公開されたGoogle Earth Engineで使えるSentinel-1の前処理法を使ってフィルタリングする。以前の続きです。 前処理法はこちら https://github.com/adugnag/gee_s1_ard
フィルタリングの方法は5種類
スペックルノイズを除去する前処理法は色々あって、今でも新しい方法が提唱されたりしているようですが、ここで紹介されているのは５種類。それぞれ１枚の画像だけをつかうのか、複数枚で時間方向にもフィルタをかけるのか？で5x10種類の処理が可能。
- Boxcar フィルタ
- Gamma maximum a-posterior (MAP) フィルタ
- Lee フィルタ
- Improved Lee Sigma speckle フィルタ
- Refined Lee speckle フィルタ
それぞれの方法で一長一短あるようで、検索すると比較する論文とか出てきますが、まあ分析対象の場所に寄りけりということなんでしょう。ざっと見た感じ最近ではImproved Lee Sigma が多く使われてるよう。ちなみにこの前処理法の論文(Mullisa et al 2021)では自分で調べてねというようなことが書いてあります。
Users can evaluate the performance of speckle filters qualitatively by visually checking the speckle reduction in homogeneous regions and the preservation of subtle features, such as roads and point scatterers, in the filtered images or quantitatively using the equivalent number of looks (ENL) or coefficient of variation in homogeneous regions in the image.
ほかに決めること
フィルタリングの窓サイズ
GEEのee.Kernelを使うので正方形のn x nの窓ではなく中心のピクセルから直径nピクセルの円として指定（数字は奇数のみ）。とりあえず5としてみました。
時間方向にフィルタリングするとき使用する画像の数
8枚分を使うとしました。これは前後まとめて３か月分のデータくらいでやってみようという見積もり。
どうやって比べるか
引用した論文の通り、画像を見た眼でみるのと均質な場所での統計量をみます。地形補正(Terrain Flattening)済みの画像を基準として、処理したものと比較。
フィルタリングした
画像をみる
ターミナル付近
この辺。
まずは１枚を処理したやつ。R：VV、G：VH、B：VV/VHで表示。
時間方向にもフィルタリングすると
東側のゴルフ場
１枚を処理
時間方向にフィルタ
滑走路南端
ここだけ色の値を変えています（値が低くて真っ暗になるから）。バンドの割り振りは同じ。 まずは１枚。
取得時期の異なる複数枚でフィルタリング。
感想
フィルタをかけるとやっぱりぼやける。ぼやけた印象が強い順にBoxcar ＞ Lee ＞ Improved Lee Sigma＞GMA ＞＝Refined Lee ＞ 処理なし。
フィルタをかけるとノイズは減る。ノイズの少ない印象順にBoxcar ＞ Lee ＞ Improved Lee Sigma ＞ GMA ＞＝Refined Lee ＞ 処理なし。
複数枚を使った場合、ノイズは増えてぼやけた印象は減る。細部の様子がよりはっきり見える。
理由は分からないけどGMAは画像を書き出すと変に切り取られた。ちょっとずれた場所なので要注意。
平均とか分散とか
とりあえずまとめ
フィルタリングするとCV（変動係数＝標準偏差/平均）は大体下がる。下がり幅（=ノイズの減り）が大きいのはLeeやImproved Lee Sigma、反対にあまり下がらないのはRefined Lee。
複数枚を使うとCVは上がる。
元々CVの高い滑走路_VVはフィルタリングすると処理前よりCVが上がることがある。これは平均的な値をとる＝ノイズに引っ張られちゃうっていうことかな。
ヒストグラムも書いてみた。
まずは処理前。
まずはBoxcar。処理前の分布は薄い色で表示。
GMA
Lee
Improved Lee Sigma
Refined Lee
まあ、こんな感じでフィルタをかけると値のばらつきが減って山が尖るということなんだね。グラフだらけになるので複数枚のグラフは載せないけど、全体として複数枚の方が山は平たんになります。
いったんまとめ
フィルタリングするとノイズは減ってエッジも減る。
程度はあるけど、ノイズだけを減らしてエッジを残すというのはできない。どっちも一緒に減る。
なので、自分の対象とする場所の特性に合わせて丁度いいバランスのフィルタ（とその設定）を探すのがポイント。
例えば分類に使うのであれば、フィルタリングした後のヒストグラムの分布の重なり具合をみて調整するとかできそう。必要以上に強いノイズ除去はせずにとどめておくような加減が探れそう。
想像するに、いろんな細かい土地被覆が複雑に分布しているところはエッジ重視。反対に、スケールの大きな、一つの土地被覆パッチが大きいところではノイズを減らすのに重きを置くといい気がする。
林縁とか土地被覆の境目（特にそこに道路が通ってたりすると）、反射の割合が変わるから突然周りと値が全然違う場所ができることがある（衛星の角度による）。なので、農地や森林の代表データを取りたいというときはそういう端っこを避けて取った方が綺麗なデータがとれるはず。
今回は建物を入れるのを忘れてた。
Google Earth Engine上でフィルタリングした画像データをAssetsにエクスポートする処理がめっちゃ時間かかる。今回のサイズで２時間程度。Googleドライブへの書き出しだと2分くらいで終わるのに…多分コンピューターも頑張らないとできない処理なんだと思う。
これだけ色々やってようやく雰囲気がつかめてきた。今回は窓サイズを変えていないけど、そっちも調整できます。
コード置き場
前回に続いて慣れないGoogle Colabを使ってPythonを書きました。少しコードを書き残しておく。相変わらずよくわかんない。
注目するポイントの周りの画像を切り出して保存する処理。画像をどう取りに行けばわからなかったので、requestsでお茶を濁したけどもっといい方法があるんだと思う：
from IPython.display import Image, display_jpeg import requests #見たい場所の座標指定。.buffer()の引数mで範囲指定。 poi_1 = ee.Geometry.Point(141.6837, 42.7878).buffer(500) # near the terminal 500m poi_2 = ee.Geometry.Point(141.7278, 42.7741).buffer(500) # golf course poi_3 = ee.Geometry.Point(141.69298, 42.76525).buffer(500) # south end of a runway #for loop の仕込み用リスト pois = [poi_1, poi_2] pois_title = ["poi_1", "poi_2"] prep = [flattened, monoBoxcar, monoGamma, monoLee, monoLeeSigma, monoRefinedLee] prepro_title = ["flattenOnly", "monoBoxcar", "monoGamma", "monoLee", "monoLeeSigma", "monoRefinedLee"] #メインの関数 def toRgbImage(img, roi, file_name): rgbImage = ee.Image.rgb(img.select('VV'), img.select('VH'), img.select('VV').divide(img.select('VH'))) #バンドのRGBへの割り振り #画像が書き出されるurlを参照 url = rgbImage.getThumbUrl({ 'min': [0.01, 0.0032, 1.25],'max': [1.0, 0.31, 31.62], 'dimensions': [512,512], 'region': roi, 'format':'jpg'}) #poi_1,2の場合の色指定 # 'min': [0.01, 0.0032, 1.25],'max': [0.1, 0.01, 31.62], 'dimensions': [512,512], 'region': roi, 'format':'jpg'}) #poi_3の場合の色指定 thumb = Image(url = url) #urlからGoogle Colabに画像を持ってくる response = requests.get(url) outputImg = response.content with open(file_name, "wb") as aaa: aaa.write(outputImg) files.download(file_name) #自分のPCにダウンロード return thumb #ループで各データ、各地点の画像を取得する for i in range(len(prep)): for j in range(len(pois)): toRgbImage(prep[i], pois[j], prepro_title[i]+'_'+pois_title[j]+'.jpg')
重ね合わせたヒストグラムを書く：
from scipy.interpolate import make_interp_spline import matplotlib.pyplot as plt import numpy as np import pandas as pd %matplotlib inline #フィルタ処理済みの使うデータのリスト prep = [flattened, monoBoxcar, monoGamma, monoLee, monoLeeSigma, monoRefinedLee] prepro_title = ["flattenOnly", "monoBoxcar", "monoGamma", "monoLee", "monoLeeSigma", "monoRefinedLee"] #注目する土地被覆のポリゴン lc = [aoiGrass, aoiForest, aoiPavement] lc_title = ['GrassLand', 'Forest', 'Pavement'] #線の色 color = ['darkorange', 'gold','darkgreen', 'greenyellow', 'red', 'salmon'] # ヒストグラムを書く関数 def combHist(prepType, pol): y_df = pd.DataFrame() for i in range(len(lc)): hist = prepType.select(pol).reduceRegion( ee.Reducer.fixedHistogram(0, 0.55, 500),lc[i]).get(pol).getInfo() # fixedHistogram(min, max, steps) a1 = np.array(hist) x1 = a1[:, 0] # array of bucket edge positions y1 = a1[:, 1]/np.sum(a1[:, 1]) # normalized array of bucket contents y_df[lc_title[i]] = y1 hist2 = flattened.select(pol).reduceRegion( ee.Reducer.fixedHistogram(0, 0.55, 500),lc[i]).get(pol).getInfo() # fixedHistogram(min, max, steps) a2 = np.array(hist2) x2 = a2[:, 0] # array of bucket edge positions y2 = a2[:, 1]/np.sum(a2[:, 1]) # normalized array of bucket contents y_df[lc_title[i]+'_flattened'] = y2 y_df['x'] = x1 return y_df # リストにしたデータごとにループ処理でグラフを書く for k in range(len(prep)): plt.figure(figsize=(10,8)) y_df = combHist(prep[k], 'VV') for j in range(6): model=make_interp_spline(y_df['x'], y_df.iloc[:, [j]], 3) xs=np.linspace(0,0.5,100) ys=model(xs) plt.plot(xs, ys, c = color[j], label = y_df.columns[j] ) plt.title(prepro_title[k]+'_VV') plt.xlim(0, 0.5) plt.legend() plt.savefig(prepro_title[k]+'_VV.png') files.download(prepro_title[k]+'_VV.png') for k in range(len(prep)): plt.figure(figsize=(10,8)) y_df = combHist(prep[k], 'VH') for j in range(6): model=make_interp_spline(y_df['x'], y_df.iloc[:, [j]], 3) xs=np.linspace(0,0.12,200) ys=model(xs) plt.plot(xs, ys, c = color[j], label = y_df.columns[j]) plt.title(prepro_title[k]+'_VH' ) plt.xlim(0, 0.12) plt.legend() plt.savefig(prepro_title[k]+'_VH.png') files.download(prepro_title[k]+'_VH.png')
SARの勉強３：スペックル その１
正直言ってよくわかっていない
スペックルノイズ、あるというのは認識できるんですが、どうやって発生するかというのがわからない。 チュートリアルをやってみたところ、自信がどんどんなくなりました。
チュートリアル自体は素晴らしい
ノイズの発生するメカニズムを数式を使って簡潔かつ華麗に解説されています。フランクフルト空港の実例を使って、どうしてノイズがある統計分布に従うかというのを説明されています。ただ最大の問題は、私の波やら電磁波に関する知識が乏しすぎて理解が追い付かないこと。最初の前提条件に平均と分散がわかればOKと書いていますが、そんなことはなかった。
結構頑張って勉強したんだけど大きくわからないのが２か所
式1.7で分母がルートNである理由：何らかの数で割らないと大変なことになるから調整のために割るんだろうというのはわかるけど、どうしてルートN？（Nじゃないの？）
式1.15の式変形。微分の絶対値をかけたりするところの進め方が付いていけない。
とはいえ、受け入れることにしました。ちゃんとガンマ分布になるし説得力があります。 全体の流れとしては、
1.あるピクセルの拡散係数には真の値に加えてノイズがのっかっている。観測値はN個のサブピクセル（的なもの）に含まれる波の和で表す。
2.サブピクセルの波は実数部虚数部それぞれ独立している。中心極限定理から両者とも平均値は正規分布に従う。で、ピクセルの波の二乗和、自由度2のカイ二乗分布に従う。
3.それをうまいこと式変形すると、観測値の確率分布を指数分布で表すことができる（=SLCの状態）。
4.Sentinel-1はマルチルック処理をしているから、3を複数回しているということ。なので、GRDは指数分布を拡張した確率分布である、ガンマ分布で近似できる。
というような感じと思われます。
自分でもやってみたい
お手本コードはGoogle Colabで開けるのでその通りできたけど、他の場所でも試してみます。 広そうな空港ということで、千歳空港。
滑走路、滑走路わきの草地、周りの平地林でVVとVHの値がどう分布しているか確認。
お手本にそってグラフを書く
基本的にお手本のとおりだけど、対象地をGeoJSONで入力するのが面倒だったので、GEEのassetに保存したFeatureCollectionを呼び出して使いました。
import ee # Trigger the authentication flow. ee.Authenticate() # Initialize the library. ee.Initialize() import matplotlib.pyplot as plt import numpy as np from scipy.stats import norm, gamma, f, chi2 import IPython.display as disp %matplotlib inline aoi = ee.FeatureCollection('Assetの保存したところ/ファイル名').geometry(); grassland = ee.FeatureCollection('Assetの保存したところ/ファイル名'); forest = ee.FeatureCollection('Assetの保存したところ/ファイル名'); pavement = ee.FeatureCollection('Assetの保存したところ/ファイル名'); aoiGrass = grassland.geometry() aoiForest = forest.geometry() aoiPavement = pavement.geometry() s1_fl = ee.Image(ee.ImageCollection('COPERNICUS/S1_GRD_FLOAT') .filterBounds(aoi) .filterDate(ee.Date('2020-08-08'), ee.Date('2020-08-10')) .first() .clip(aoi)) def histFunction(aoi_sub, pol, title): hist = s1_fl.select(pol).reduceRegion( ee.Reducer.fixedHistogram(0, 0.5, 500),aoi_sub).get(pol).getInfo() # fixedHistogram(min, max, steps) mean = s1_fl.select(pol).reduceRegion( ee.Reducer.mean(), aoi_sub).get(pol).getInfo() variance = s1_fl.select(pol).reduceRegion( ee.Reducer.variance(), aoi_sub).get(pol).getInfo() beta = variance/mean alpha = mean/beta a = np.array(hist) x = a[:, 0] # array of bucket edge positions y = a[:, 1]/np.sum(a[:, 1]) # normalized array of bucket contents fig, ax = plt.subplots(1,1) plt.grid() ax.plot(x, y, '.', label='data') plt.plot(x, gamma.pdf(x, alpha, 0, beta)/1000, '-r', label='gamma') plt.legend() plt.title(title + " " + pol) plt.text(0.7, 0.3, "mean = {}".format(round(mean,4)), transform=ax.transAxes) plt.text(0.7, 0.2, "variance = {}".format(round(variance,5)), transform=ax.transAxes) return plt.show() histFunction(aoiGrass, 'VV', 'GrassLand') #ヒストグラムの描写
ヒストグラムを見る
大体あってる。すごい。ガンマ分布最高。
というわけで、次回以降はフィルタリングするとどうなるかを見てみたいと思います。
SARの勉強２：入射角の影響
SAR勉強シリーズを始める
ちょっと勉強したいと思ったら結構なボリュームを勉強しないと分からないことが片付かないSAR。ちょっとずつやれば、そのうち理解が深まるんじゃないかと期待してシリーズを始めることに。まず、前回ノイズとエラーの原因として挙げたことから見てみようかと。で、まあ簡単に見られそうな入射角から取り組んでみます。
問題のおさらい
Google Earth Engineで提供されているSARデータは、衛星から斜め下方に照射されるマイクロ波のうち、衛星に帰ってきたものの強度（とタイミング的なやつが本当はあるけど入っていない）を測った結果をきれいにして画像にしたもの。衛星からはある幅をもってマイクロ波を照射するため、衛星に近い側と遠い側では（地表面が平坦であるとき）入射角が変わる。波の反射する角度やら反射しやすさやらは入射角によって変化するので、地表の状態が同じであっても（つまり土地被覆がおなじでも）、SARの画像データの値が違うはず。それはどうなのということ。
実際に画像を見てみると、偏波VVのSentinel-1画像の海上では境目に線がはっきり見える。この線の東西で地表の状態が違うことは考えにくいから、入射角の影響だと見て取れる。ただ陸上はよくわかんないな。それでVHの方には見えない。
検証方法
土地被覆が同じで傾斜がない場所であれば（本来の）後方散乱強度は同じはず、というわけなんですが、検証できるのは広大な場所でないといけない。Sentinel-1の一枚の画像は数百キロくらいだから、取りうる角度を網羅しようと思ったら相当広くないとしんどい。というわけで、３か所を考えてみました。 それぞれの場所で、Sentinel-1画像をモザイクして得られる値（ここではσ0）と入射角の散布図を描きます1。その後、前回紹介したTerrain Flatteningの前処理（下記リンク）を施して同じような散布図を描いて、補正処理の効果をみます。散布図に線を引っ張ってみましたが検定はかけてません。
アマゾンの森林
アマゾン広大だし、平坦な場所多そうだし。
オーストラリア南部の荒野
検索したら世界一広い平地と出てきたので。多分砂漠だと思います。
五大湖のひとつ、スペリオール湖
横に広そうな湖ということで。
まとめ
というわけで、まとめると
- 入射角が大きいほど値は小さくなる。
- 入射角の影響はVVで一番大きく次にVH。VV/VHの比をとると影響はほとんどなさそう。
- 前処理（Terrain Flattening）をすると、VVでは入射角の影響が緩和される。VHは微妙。VV/VHは元々影響がなさそうなので前後での傾向に変化は少ない。
さらに実用を踏まえて付け加えると
- 一枚の画像で小面積（例えば市町村とか以下のサイズ）を対象とした分析をする場合、入射角はそれほど変わらないのであまり気にする必要はない。ただし、画像が複数枚のモザイクの場合、あるいは衛星の進行方向の違うデータ（ASCENDINGとDESCENDING）を混ぜて使う場合は別。
- そもそも森林、草地灌木（砂漠のイメージ）、水域でかなりσ0の値が違う。この土地被覆による違いと比べると入射角の影響による違いは随分小さい。なので分類のためにSentinel-1を使うときは、入射角の前処理でそんなにカリカリしなくても分類できちゃいそう。
- 分類ではなく、たとえばσ0と植被率を回帰するとか（できるか知らないけど）の場合は当然補正をすべし。
- ものすごく平坦な場所とかでない限り、どうせTerrain Flattening （地形補正）はする。その時に入射角の影響が緩和されるならよし。
今回のコード
前で紹介しているMullissa et al. 2021のコードを使わせてもらっています。GitHubからリンクを辿ると前処理用のコードがコピーされるので、それをrequireで読み込んでいます。上の8行目までを適当に変えれば、３箇所の結果の散布図が出てきます。
https://code.earthengine.google.com/8f48feec90c619deb614489765718072
-
なお普通にデータカタログに出てくるSentinel-1はlog値（デシベル）のやつなんですが、前処理のお手本コードでは小数点表記のものを使っていたのでそれに乗っかってみました。
この、FLOATというσ0のデータは検索しても出てこなくって、いったいどこから出てきてるのかよくわかりませんでした。FLOATのデータの説明も公式にちゃんとありました（一番下に" If you want to use the underlying collection with raw power values (which is updated faster), see COPERNICUS/S1_GRD_FLOAT.“との記述）。↩
