SARの勉強３：スペックルその１

正直言ってよくわかっていない

スペックルノイズ、あるというのは認識できるんですが、どうやって発生するかというのがわからない。チュートリアルをやってみたところ、自信がどんどんなくなりました。

developers.google.com

チュートリアル自体は素晴らしい

ノイズの発生するメカニズムを数式を使って簡潔かつ華麗に解説されています。フランクフルト空港の実例を使って、どうしてノイズがある統計分布に従うかというのを説明されています。ただ最大の問題は、私の波やら電磁波に関する知識が乏しすぎて理解が追い付かないこと。最初の前提条件に平均と分散がわかればOKと書いていますが、そんなことはなかった。

結構頑張って勉強したんだけど大きくわからないのが２か所

式1.7で分母がルートNである理由：何らかの数で割らないと大変なことになるから調整のために割るんだろうというのはわかるけど、どうしてルートN？（Nじゃないの？）
式1.15の式変形。微分の絶対値をかけたりするところの進め方が付いていけない。

とはいえ、受け入れることにしました。ちゃんとガンマ分布になるし説得力があります。全体の流れとしては、

1.あるピクセルの拡散係数には真の値に加えてノイズがのっかっている。観測値はN個のサブピクセル（的なもの）に含まれる波の和で表す。

2.サブピクセルの波は実数部虚数部それぞれ独立している。中心極限定理から両者とも平均値は正規分布に従う。で、ピクセルの波の二乗和、自由度2のカイ二乗分布に従う。

3.それをうまいこと式変形すると、観測値の確率分布を指数分布で表すことができる（=SLCの状態）。

4.Sentinel-1はマルチルック処理をしているから、3を複数回しているということ。なので、GRDは指数分布を拡張した確率分布である、ガンマ分布で近似できる。

というような感じと思われます。

自分でもやってみたい

お手本コードはGoogle Colabで開けるのでその通りできたけど、他の場所でも試してみます。広そうな空港ということで、千歳空港。

滑走路、滑走路わきの草地、周りの平地林でVVとVHの値がどう分布しているか確認。

f:id:camarao:20210917122501p:plain — 滑走路（赤）、草地（オレンジ）、平地林（黄緑）のピクセルで確認

お手本にそってグラフを書く

基本的にお手本のとおりだけど、対象地を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') #ヒストグラムの描写