宙畑では、衛星観測技術である合成開口レーダ(Synthetic Aperture Radar: 以下SAR)についていくつか解析事例を紹介しています。

その中でも、「偏波（へんぱ）」と呼ばれるSARの観測方法を使って都市や自然物、湖などの構造物や自然物の推定した記事があります。

SAR画像の偏波を使った構造物推定（偏波の成分分解）
https://sorabatake.jp/13778/

SAR画像を使って、構造物と自然物を識別する！四成分散乱分解
https://sorabatake.jp/14537/

これらの記事では一時期のSAR画像の偏波情報(HH,HV,VV偏波)を使用して構造物や自然物を可視化していて、時系列の変化を見てはいません。

時期の異なる複数のSAR画像を用いて、時系列解析を試みた事例としては差分干渉がありますが、こちらは単一偏波の情報のみを扱った解析になります。また、差分干渉自体が時期の異なる複数の画像間での地盤変化の可視化に利用される手法であり、構造物の推定に用いられるものではありません。

そこで、本記事では複数枚のSAR画像を使った時系列間の違いの検証や、複数枚のSAR画像の偏波情報を用いた構造物、自然物の推定を試みたいと思います。

具体的な内容は以下の通りです。

・二枚のSARの強度画像の位置合わせを行い、差分検出(Change Detection)をしてみる
・差分検出での位置合わせに用いたパラメータを使ってSARの生データを位置合わせして最大相関化アルゴリズムによるコヒーレンス最適化(Coherence Optimization)を行い、得られたコヒーレンス情報を可視化してみる

また、今回の解析はTellusで提供されている一般的な仮想環境で実行することが可能ですので、ローカルで色々と開発環境をセットする必要もなくデータ処理することが可能です。興味がある方はぜひTellusの開発環境に登録してSARデータを触ってみてください。

☆Tellusの登録はこちらから
https://www.tellusxdp.com/ja/developer/

使用した仮想環境(さくらのクラウド )

CPU：4Core
メモリ：8GB
ディスク：SSD 100GB

差分検出(Change Detection)

必要なフルポラリメトリックSAR画像のペアを取得する

差分検出とは単純な時系列間のSAR画像を比較するための手法で、SARの強度画像の差分を取ることで、時系列間の間にどのような変化が起きたかを可視化する手法です。

まずは二枚のSAR画像がないと検証できないので、以下のようなソースコードで同じ撮像地域のSAR画像を探していきました。

今回用いるSAR画像は、4種類の偏波情報を有している「フルポラリメトリック」、略して「フルポラ」の画像です。具体的には’polarisations’というパラメータが’HH+HV+VV+VH’となっている画像です。

import os, requests
import geocoder # こちらはTellusの仮想環境にデフォルトで入っていないので、!pip install geocoderでgeocoderをインストールしてください

# Fields
BASE_API_URL = "https://file.tellusxdp.com/api/v1/origin/search/" # https://www.tellusxdp.com/docs/api-reference/palsar2-files-v1.html#/
ACCESS_TOKEN = "YOUR-TOKEN"
HEADERS = {"Authorization": "Bearer " + ACCESS_TOKEN}
TARGET_PLACE="tokyo"
SAVE_DIRECTORY="./SAR_{}/".format(TARGET_PLACE)

# Functions
def rect_vs_point(ax, ay, aw, ah, bx, by):
    return 1 if bx > ax and bx < aw and by > ay and by < ah else 0
def get_scene_list_slc(_get_params={}):
    query = "palsar2-l11"
    r = requests.get(BASE_API_URL + query, _get_params, headers=HEADERS)
    if not r.status_code == requests.codes.ok:
        r.raise_for_status()
    return r.json()

def get_slc_scenes(_target_json, _get_params={}):
    # get file list
    r = requests.get(_target_json["publish_link"], _get_params, headers=HEADERS)
    if not r.status_code == requests.codes.ok:
        r.raise_for_status()
    file_list = r.json()['files']
    dataset_id = _target_json['dataset_id'] # folder name
    # make dir
    if os.path.exists(SAVE_DIRECTORY + dataset_id) == False:
        os.makedirs(SAVE_DIRECTORY + dataset_id)
    # downloading
    print("[Start downloading]", dataset_id)
    for _tmp in file_list:
        r = requests.get(_tmp['url'], headers=HEADERS, stream=True)
        if not r.status_code == requests.codes.ok:
            r.raise_for_status()
        with open(os.path.join(SAVE_DIRECTORY, dataset_id, _tmp['file_name']), "wb") as f:
            f.write(r.content)
        print("  - [Done]", _tmp['file_name'])
    print("finished") 
    return

# Entry point
def main():
    # extract slc list around the address
    gc = geocoder.osm(TARGET_PLACE, timeout=5.0) # get latlon
    slc_list_json = get_scene_list_slc({'after':'2019-1-01T00:00:00', 'polarisations':'HH+HV+VV+VH', 'page_size':'1000'})  # get single look complex list
    target_places_json = [_ for _ in slc_list_json['items'] if rect_vs_point(_['bbox'][1], _['bbox'][0], _['bbox'][3], _['bbox'][2], gc.latlng[0], gc.latlng[1])] # lot_min, lat_min, lot_max...
    target_ids = [_['dataset_id'] for _ in target_places_json]
    print("[Matched SLCs]", target_ids)
    for target_id in target_ids:
        target_json = [_ for _ in slc_list_json['items'] if _['dataset_id'] == target_id][0]
        # download the target file
        get_slc_scenes(target_json)
        #break # if remove, download all data
    
if __name__=="__main__":
       main()

暫く待つと、以下のような出力結果が返ってきます。

[Matched SLCs] ['ALOS2250690700-190113', 'ALOS2250912890-190115', 'ALOS2251130730-190116', 'ALOS2265402890-190423', 'ALOS2267250700-190505', 'ALOS2267472890-190507']
[Start downloading] ALOS2250690700-190113
  - [Done] BRS-HH-ALOS2250690700-190113-HBQR1.1__A.jpg
  - [Done] BRS-HV-ALOS2250690700-190113-HBQR1.1__A.jpg
  - [Done] BRS-VH-ALOS2250690700-190113-HBQR1.1__A.jpg
  - [Done] BRS-VV-ALOS2250690700-190113-HBQR1.1__A.jpg
  - [Done] IMG-HH-ALOS2250690700-190113-HBQR1.1__A
  - [Done] IMG-HV-ALOS2250690700-190113-HBQR1.1__A
  - [Done] IMG-VH-ALOS2250690700-190113-HBQR1.1__A
  - [Done] IMG-VV-ALOS2250690700-190113-HBQR1.1__A
  - [Done] LED-ALOS2250690700-190113-HBQR1.1__A
  - [Done] TRL-ALOS2250690700-190113-HBQR1.1__A
  - [Done] VOL-ALOS2250690700-190113-HBQR1.1__A
  - [Done] summary.txt
finished
(以下略)

ここで、[Matched SLCs]に書かれているALOS〜から続く配列に注目してください。こちらが今回の処理でTARGET_PLACE=tokyoの周辺で取得できたフルポラSAR画像の情報が記録されています。ハイフンのあとは撮像された日時を示していますが、その前の4桁(例えば、ALOS2250690700-190113 にある0700の部分)が撮影範囲の中心地域を示しています。この中心地域がなるべく同じSAR画像同士で解析するのが望ましいので、今回は以下の二つのSAR画像で解析を行いました。

target_id1 = ‘ALOS2267472890-190507′
target_id2 = ‘ALOS2250912890-190115′

PALSAR-2の観測パラメータなどの詳細は公式サイトに掲載されているのでご覧ください。
定常運用期間の観測実績について

では、使用するフルポラリメトリックSAR画像のペアが見つかったので、画像を切り出していきます。今回は千葉県野田市周辺を使うことにしました。

import glob
import numpy as np
import re
import struct

target_id1 = 'SAR_tokyo/ALOS2267472890-190507'
target_id2 = 'SAR_tokyo/ALOS2250912890-190115'

class Get_image():
    def __init__(self, dataset_id=''):
        self.dataset_id = dataset_id
    
    def SAR_data(self):
        names = ['HH','HV','VV']
        lists = glob.glob(self.dataset_id + '/*')
        CEOS_lists = sorted({k for k in lists if re.search('IMG-(HH|HV|VV)', str(k))})
        i=0
        for filename in CEOS_lists:
            fp = open(os.path.join('/home','jovyan','work',filename),mode='rb')
            col = 2000
            row = 2000
            cr = np.array([3000,3000+col,8000,8000+row],dtype="i4");
            
            fp.seek(236)
            nline = int(fp.read(8))
            fp.seek(248)
            ncell = int(fp.read(8))
            nrec = 544 + ncell*8

            nline = cr[3]-cr[2]
            fp.seek(720)
            fp.read(int((nrec/4)*(cr[2])*4))
            data = struct.unpack(">%s"%(int((nrec*nline)/4))+"f",fp.read(int(nrec*nline)))
            data = np.array(data).reshape(-1,int(nrec/4))
            data = data[:,int(544/4):int(nrec/4)]
            slc = data[:,::2] + 1j*data[:,1::2]
            slc = slc[:,cr[0]:cr[1]]
            print(np.shape(slc))
            np.save(os.path.join('/home','jovyan','work',self.dataset_id,names[i]),slc)
            i+=1


SAR_get = Get_image(target_id1)
SAR_get.SAR_data()
print('GET SAR SLC DATA AS NUMPY!!')

SAR_get = Get_image(target_id2)
SAR_get.SAR_data()
print('GET SAR SLC DATA AS NUMPY!!')

ポラリメトリック補正と8bitの正規化をかけて、コントラストの調整などを行うと以下のような形でPALSAR-2の疑似カラー画像を取得できます。

ちなみに、対象地域をGoogleマップで確認すると以下の地域になります。真ん中にある江戸川を境にして左側が埼玉県越谷市、右側が千葉県野田市になります。

差分検出と位置合わせの重要性

次に先ほど取得したSAR画像のペアを使用して差分検出を行なっていくのですが、まずは位置合わせを行わずにSARの強度画像による差分検出の結果をお見せしたいと思います(ソースコードは位置合わせのバージョンも含めて後述しています)。

この画像は単純に強度画像同士の引き算をしているだけなので、全く同じ画像であれば真っ黒（ピクセルの値がゼロ）になるものです。

幸い二つのSAR画像はかなり同じ地点を撮像していたみたいなので、無事に差分が取れていると思いますが、よく見ていただけると分かるように白い小さな細いスジのようなものが確認できるかと思います。これは位置合わせが正確でないことを示していて、ほとんど同じ場所を撮像しているのにも関わらず、差分検出を行なってみると僅かに歪みが発生していることが確認できるかと思います。この歪みに関しては画像同士を正確に位置合わせすることによって軽減させることができます。

今回はSAR画像でよく使われる回転不変位相限定相関法(RIPOC: Rotation Invariant Phase Only Correlation)という手法で位置合わせを行いました。RIPOCに関しては独自で実装しようと思ったのですが、以下のサイトにopencvと組み合わせて簡単に実装されていたので参考にさせて頂きました。

RIPOC – 二枚の画像間の傾きと縮尺を求める – kakasi’s blog

target_id1 = 'SAR_tokyo/ALOS2267472890-190507'
target_id2 = 'SAR_tokyo/ALOS2250912890-190115'


# ここで呼び出してるライブラリは後のコヒーレンス最適化の処理でも使用します
from scipy.optimize import leastsq
import scipy, scipy.fftpack
import cv2
import matplotlib.pyplot as plt
from numpy import inf
import numpy as np

def ripoc(a, b, m = None):
    g_a = np.asarray(cv2.cvtColor(a, cv2.COLOR_BGR2GRAY), 'float')
    g_b = np.asarray(cv2.cvtColor(b, cv2.COLOR_BGR2GRAY), 'float')

    h, w = g_a.shape
    hy = np.hanning(h)
    hx = np.hanning(w)
    hw = hy.reshape(h, 1)*hx
    
    f_a = np.fft.fftshift(np.log(np.abs(np.fft.fft2(g_a*hw))))
    f_b = np.fft.fftshift(np.log(np.abs(np.fft.fft2(g_b*hw))))

    if not m:
        l = np.sqrt(w*w + h*h)
        m = l/np.log(l)

    center = (w/2, h/2)
    flags = cv2.INTER_LANCZOS4 + cv2.WARP_POLAR_LOG
    p_a = cv2.warpPolar(f_a, (w, h), center, m, flags)
    p_b = cv2.warpPolar(f_b, (w, h), center, m, flags)
    (x, y), e = cv2.phaseCorrelate(p_a, p_b, hw)

    angle = y*360/h
    scale = (np.e)**(x/m)
    M = cv2.getRotationMatrix2D(center, angle, scale)
    t_b = cv2.warpAffine((g_b), M, (w, h))
    (x, y), e = cv2.phaseCorrelate(g_a, t_b)

    return x, y, angle, scale

def scatter_normalize(data):
    data = 10*np.log10(data) -83.0 -32.0
    data = np.array(255*(data-np.amin(data))/(np.amax(data)-np.amin(data)),dtype="uint8")
    return data

def make_image(HH_path,HV_path,VV_path,flag):
    Shh = scatter_normalize(np.load(HH_path))
    Shv = scatter_normalize(np.load(HV_path))
    Svv = scatter_normalize(np.load(VV_path))
    img=np.zeros((2000, 2000, 3), np.uint8)
    img[:,:,0]=Shh
    img[:,:,1]=Shv
    img[:,:,2]=Svv
    return cv2.flip(img, flag)

def registration(img1,img2,flag):
    # 位置合わせをしなかったとき
    if flag==0:
        sigma = np.array(img1,dtype=float)-np.array(img2,dtype=float)
        img = np.array(255*(sigma-np.amin(sigma))/(np.amax(sigma)-np.amin(sigma)),dtype="uint8")
        plt.imshow(img)
        plt.imsave('raw_change.png', img)
    # 位置合わせを行うとき
    if flag==1:
        x, y, angle, scale = ripoc(img1, img2)
        print(x, y, angle, scale)

        h, w, ch = img1.shape
        M = cv2.getRotationMatrix2D((w/2,h/2), angle, scale)
        M[0][2] -= x
        M[1][2] -= y
        dst = cv2.warpAffine(img2, M, (w, h))
        
        sigma = np.array(img1,dtype=float)-np.array(dst,dtype=float)
        img = np.array(255*(sigma-np.amin(sigma))/(np.amax(sigma)-np.amin(sigma)),dtype="uint8")
        plt.imshow(img)
        plt.imsave('ripoc_change.png', img)
        return M, h, w
        

img1 = make_image('/home/jovyan/work/'+ target_id1 + '/HH.npy',
            '/home/jovyan/work/'+ target_id1 + '/HV.npy',
            '/home/jovyan/work/'+ target_id1 + '/VV.npy',
            1)
img2 = make_image('/home/jovyan/work/'+ target_id2 + '/HH.npy',
            '/home/jovyan/work/'+ target_id2 + '/HV.npy',
            '/home/jovyan/work/'+ target_id2 + '/VV.npy',
            1)

# 位置合わせしない時の差分検出
registration(img1,img2,0)
# RIPOCで位置合わせした時の差分検出
M, h, w = registration(img1,img2,1)

上の図が差分検出の結果になりますが、位置合わせを行った後の方が余計な白いスジがなくなり、正確な差分検出ができていることが確認できると思います。

画像の位置合わせの変化量ですが、今回のケースではx方向のズレが-0.6671478707661436ピクセル、y方向のズレが1.140806499636824ピクセルでした。些細な変化だと思いますが、例えば差分干渉解析を行う場合には顕著なノイズを生じてしまう恐れがあるので、SAR画像解析においては重要な処理になります。

では、実際に差分検出で得られた結果を見ていきましょう。今回は冬と春での差分検出の結果になりますので主に森林地域では差分が確認できますが、右下の部分に他とは明らかに違う変化が生まれていることが確認できるかと思います(下図の赤枠で囲まれている部分)。

こちらに関してGoogleマップで確認してみたところ、以下の地域での変化であることが推測できました。

こちらについてですが、どうやら千葉県流山市の開発事業による工事が行われており、その時の工事状況による変化であることが考えられました。

今回は位置合わせを行ってSAR画像同士を差分検出しただけの単純な処理ですが、これだけでも多くの情報を取得することができます。対象地域の細かい解析を行う前に変化検出を行って変化のある地域を特定してから、その地域の偏波シグネチャを取得して散乱特性の違いを把握してグランドトゥルースする手法は今でも行われてるSAR画像の代表的な解析方法の一つだと思います。

最大相関化法を用いてコヒーレンスを最適化し、コヒーレンス情報を可視化してみる

では、先ほどの位置合わせで行ったパラメータ(M, h, w)を使用して複素数の情報を保持している生のSARデータの位置合わせを行い、コヒーレンスを最適化していきたいと思います。

コヒーレンス（coherence）とは、干渉性とも訳され、電波がどれだけ同じ形を保っているかを示す言葉です。

今回は二時期の観測に対しての「コヒーレンス」なので、二時期の観測データがどれほど似ているかと意味しています。

コヒーレンスは二つのSARデータの同じ位置の複素相関で表現することができます。SAR画像の複素相関は同じ座標位置の散乱成分とその複素共役を乗算することで求められます。
相関というくらいなので、コヒーレンスが高い地域は二つのSARデータで変化がなかった地域(主に人工物など)が該当します。一方でコヒーレンスが低い地域は二つのSARデータの時系列間で変化起きた地域(河川や表面散乱するような地域)が該当します。

前回の宙畑の記事では差分干渉と呼ばれる処理を行うことで、干渉性を評価して地形の変化情報を取り出していました。今回は土地被覆分類の一環として人工物の検出を行うためにコヒーレンス情報を取得したいと思います。

また、純粋に二つのSAR画像の複素相関を算出するのではなく、二つのSAR画像の相関を最大化する位相状態を推定する手法である相関最大化アルゴリズムを適用することによってコヒーレンスの最適化を行うことにしました。

では、まず始めに複素情報を所持してる生のSAR画像の位置合わせとその時の差分結果を表示してみます。先ほどのソースコードに続いて、以下のコードを実行します。

# SARの生データのレジストレーション
HH1 = np.load('/home/jovyan/work/'+ target_id1 + '/HH.npy')
HV1 = np.load('/home/jovyan/work/'+ target_id1 + '/HV.npy')
VV1 = np.load('/home/jovyan/work/'+ target_id1 + '/VV.npy')
HH2 = np.load('/home/jovyan/work/'+ target_id2 + '/HH.npy')
HV2 = np.load('/home/jovyan/work/'+ target_id2 + '/HV.npy')
VV2 = np.load('/home/jovyan/work/'+ target_id2 + '/VV.npy')

def affine(data, M, h, w):
    # 8bitに正規化しないとRIPOCが上手く行かないので、
    # 正規化した方で変換パラメータだけ取得して、生データに対して変換パラメータを適応させる
    result = np.zeros((np.shape(data)[0],np.shape(data)[1]), dtype=np.complex)
    result.real = cv2.warpAffine(data.real, M, (w, h))
    result.imag = cv2.warpAffine(data.imag, M, (w, h))
    return result

# slaveのSAR画像の位置合わせを行う
t_HH2 = affine(HH2, M, h, w)
t_HV2 = affine(HV2, M, h, w)
t_VV2 = affine(VV2, M, h, w)


#cv2.flip()のflagの設定について：
#flag = 0: 上下反転
#flag > 0: 左右反転
#flag < 0: 上下左右反転
def make_raw_image(HH,HV,VV,flag):
    Shh = 10*np.log10(HH) -83.0 -32.0
    Shv = 10*np.log10(HV) -83.0 -32.0
    Svv = 10*np.log10(VV) -83.0 -32.0
    img=np.zeros((2000, 2000, 3), np.uint8)
    img[:,:,0]=Shh
    img[:,:,1]=Shv
    img[:,:,2]=Svv
    return cv2.flip(img, flag)

master = make_raw_image(HH1,HV1,VV1,1)
slave = make_raw_image(t_HH2,t_HV2,t_VV2,1)

sigma = np.array(master,dtype=float)-np.array(slave,dtype=float)
img = np.array(255*(sigma-np.amin(sigma))/(np.amax(sigma)-np.amin(sigma)),dtype="uint8")
plt.imshow(img)
plt.imsave('change_SAR_complex.png', img)

実際に位置合わせを行った結果は以下のようになります。

ほぼ真っ黒な画像を得ることができていることから、正常に位置合わせができていることが確認できました。

では、次に相関最大化アルゴリズムによってコヒーレンスの最適化を行っていきます。元論文はこちらになります。
Polarimetric SAR interferometry – IEEE Journals & Magazine

こちらの手法についてですが、以下の手順で実装を行いました。

1. それぞれのSAR画像をパウリ分解する

2. パウリ分解して得られた偏波ベクトルの外積を用いて、三種類の共分散行列を生成する

3.干渉コヒーレンスが最大となるように、ラグランジュの未定乗数法を用いて最適化を行う
干渉コヒーレンスについては以下の式で表現することができます。

こちらの最適化の手順に関しては、具体的に以下の式を計算して最適化すると元論文などにも書かれています。

こちらの解法については機械学習などで知られている正準相関分析の解法と同じ手順で求めれば大丈夫です。具体的には、以下の行列の最大の特異値に対応する特異ベクトルを求めることによって達成できます。

4.最後に、得られた特異ベクトルの情報を使って以下の式でコヒーレンス情報を計算します

具体的なソースコードは以下のようになります。こちらの処理時間については、今回の解析対象地域であれば10分以内に処理が終わります。

def Complex_coherence(Shh1,Shv1,Svv1,Shh2,Shv2,Svv2,x,y):
    
    T11 = np.zeros((x,y,3,3), dtype=np.complex)
    T22 = np.zeros((x,y,3,3), dtype=np.complex)
    Omega12 = np.zeros((x,y,3,3), dtype=np.complex)
    
    for i in range(x):
        for j in range(y):
            k11 = Shh1[i,j]+Svv1[i,j]
            k12 = Svv1[i,j]-Shh1[i,j]
            k13 = 2*Shv1[i,j]
            k21 = Shh2[i,j]+Svv2[i,j]
            k22 = Svv2[i,j]-Shh2[i,j]
            k23 = 2*Shv2[i,j]
            
            k1 = np.array([k11, k12, k13]).reshape(3,1)
            k2 = np.array([k21, k22, k23]).reshape(3,1)
            
            T11[i,j,:,:] = np.outer(k1, np.conjugate(k1).T)
            T22[i,j,:,:] = np.outer(k2, np.conjugate(k2).T)
            Omega12[i,j,:,:] = np.outer(k1,np.conjugate(k2).T)
            
    result = np.zeros((x,y), dtype=np.complex)
    for i in range(x):
        for j in range(y):     
            Sigma = np.dot(np.dot(np.sqrt(T11[i,j,:,:]),Omega12[i,j,:,:]),np.sqrt(T22[i,j,:,:]))
            U, S, V = np.linalg.svd(Sigma, full_matrices=False)
            result[i,j] = np.dot(np.dot(np.conjugate(U[0]).T,Omega12[i,j,:,:]),V[0])
    return result

result = Complex_coherence(HH1,HV1,VV1,t_HH2,t_HV2,t_VV2,np.shape(HH1)[0],np.shape(HH1)[1])

では、実際のコヒーレンス情報(インターフェロメトリック位相)を可視化してみます。

# コヒーレンスの高い地域を可視化
coh = 10*np.log10(abs(result)) -83.0 -32.0
coh = np.angle(np.array(coh, dtype=np.complex64))
coh[coh == -inf] = 0
print(coh)
# 8bitで正規化して余計なノイズをなくす
coh = np.array(255*(coh-np.amin(coh))/(np.amax(coh)-np.amin(coh)),dtype="uint8")
print(coh)
img = cv2.flip(np.array(coh, dtype=np.float64), 1)
plt.imshow(np.array(img, dtype=np.float64))
plt.imsave('phase.jpg', img, cmap = "jet")

上図の結果で青く表示されている部分がコヒーレンスの高い地域(人工物など)になります。先ほどお見せした疑似カラー画像と重ね合わせた結果を下図に示していますが、HH偏波(赤色)が強い値を持つ建物などの人工物が青色で綺麗に表示できていることが確認できました。

このように、二つのフルポラリメトリSAR画像を組み合わせることで、コヒーレンスが強くなりやすい人工物のある地域だけを綺麗に可視化することができました。

この情報を使用することによって、大まかな都市域の面積を地域ごとに比較したり、時代ごとの人工物の変化を追跡することが可能です。また、都市域の大まかなラベリング情報としての使い道としても利用することができるので、機械学習を用いた土地被覆分類の特徴量としても利用することができると思います。

さらにコヒーレンス情報の使い道はこれだけではなく、差分干渉による地盤変動の解析にも使用されています。

また、先ほどの相関最大化アルゴリズムで得られたコヒーレンス情報と差分干渉の過程で取り出した高度情報を組み合わせることで、バイオマスの推定などにも応用することができます(多偏波干渉SARと言います)。

まとめと感想

今回は複数のフルポラリメトリSAR画像を組み合わせて変化検出を行い、相関最大化アルゴリズムを用いてコヒーレンス情報を取り出して人工物のある地域の可視化を行いました。

差分検出は位相情報を扱わず単純な強度画像の差分で表示するだけのシンプルな解析手法ですが、SAR画像を用いることで時系列間の本来気づきにくい変化を簡単に観測することができるのでオススメです。また、コヒーレンス情報も大まかに建物のある地域だけを把握したいときに活かすことができると思います。

今回ご紹介した手法の応用である「多偏波干渉SAR」も今後チャレンジしてみたいと思います。