PALSAR2を使った多偏波SARの解析手法(ミュラー行列+マルチルック処理)

はじめに

お久しぶりです、新星竹です。久しぶりにSAR解析をしました。取り組んでいる研究の解析やらジャーナルの執筆やらで佳境に入っていて、尚且つ新しい委託のお仕事を探していたり決まった後にフルコミットしたりしていた関係で、SAR画像解析をする時間が中々作れませんでした汗

今までは、TellusのGPUコンピューティング環境のセットアップからSAR画像のバイナリデータを解析してRGBで散乱特性の可視化などを行なってきました。今回は、SAR画像が解析に使えるように前処理工程を実装します。せっかくフルポラリメトリ解析ができる環境がTellusにはあるので、今回は電波が入射した表面の偏波散乱特性を表現するミュラー行列をSAR画像のデータから算出して、ミュラー行列の特徴としてパワーに関する表現であることから加算平均が可能なので、マルチルック処理を並行して行なって後方散乱係数を取得します。マルチルック処理とは、コヒーレントな電波を用いて対象に電波を送信して対象物を計測する場合、原理的に生じる画素値のノイズ(スペックル)を低減させるために行われるものです。最終的に、散乱波のパワーに関するSARデータは、対象の単位面積当たりの散乱波のパワーである後方散乱係数によって記述されます。

この過程で、SAR画像の前処理工程は完了するので、そこからシグネチャを抽出して対象の散乱の振る舞い方を解析したり、様々な画像分類アルゴリズムを使用して土地の被覆分類などを行うことが可能になります！(つまり、いろんな解析ができます笑)

実装アルゴリズム

この解析は、画像をクロップして切り出さずに解析を行なった場合(つまり、TellusのAPIからHH,HV,VV成分をそのまま取得したデータで解析を行なった場合、かなりの解析時間がかかります(ミュラー行列を算出するだけで3時間、そこからマルチルック処理が90分くらいかかるので、全体で一つのSAR撮影画像を前処理するのに5時間もかかります笑)。なので、この処理は画像を切り抜いた状態で解析を行うか、サーバ上でジョブを投げて管理しないとTellusの開発環境上では実践しにくいので注意してください汗

また、globで取り出しているのは.npyとしてバイナリデータを処理したものを使っています。SARのバイナリデータからnumpyへの変換方法は前回の記事を参考にしてください。また、今回のマルチルック処理は、アジマス方向に4ルックのマルチルック処理を施しました(つまり、縦方向に1/4に圧縮されています)。

yukr.hatenablog.com

以下、ミュラー加算平均アルゴリズムの実装(python)

import requests
from PIL import Image
import urllib.request
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import codecs
import struct
from io import BytesIO
from PIL import ImageOps
import re
import time
import os

import glob
SAR_list = glob.glob('~/dataset/*.npy')

i=0
res = []
for data in SAR_list:
    SAR = np.load(data)
    print(SAR.shape[0],SAR.shape[1],SAR.dtype)
    res.append(SAR)
    i=i+1
res = np.array(res)


def multilook(im):
    original_width = im.shape[1]
    original_height = im.shape[0]
    width = original_width
    height = original_height // 4
    resized_image = np.zeros(shape=(height, width), dtype=np.float64)
    scale = 4

    for i in range(height):
        for j in range(width):
            temp = 0
            for x in range(scale):
                for y in range(scale):
                    temp += im[i*scale + x, j]
            resized_image[i, j] = temp/(scale*scale)

    return resized_image


def mueller(HH,HV,VV,x,y):
    
    m11 = np.zeros((x,y))
    m12 = np.zeros((x,y))
    m13 = np.zeros((x,y))
    m14 = np.zeros((x,y))
    m21 = np.zeros((x,y))
    m22 = np.zeros((x,y))
    m23 = np.zeros((x,y))
    m24 = np.zeros((x,y))
    m31 = np.zeros((x,y))
    m32 = np.zeros((x,y))
    m33 = np.zeros((x,y))
    m34 = np.zeros((x,y))
    m41 = np.zeros((x,y))
    m42 = np.zeros((x,y))
    m43 = np.zeros((x,y))
    m44 = np.zeros((x,y))
    
    for i in range(x):
        for j in range(y):
            a = HH[i,j]
            b = a.conjugate()
            c = HV[i,j]
            d = c.conjugate()
            e = VV[i,j]
            f = e.conjugate()

            m11[i,j] = (0.25 * (a*b + 2*c*d + e*f)).real
            m12[i,j] = (0.25 * (a*b - e*f)).real
            m13[i,j] = (0.5 * ((a*d).real + (c*f).real)).real
            m14[i,j] = (0.5 * ((c*b).imag + (e*d).imag)).real

            m21[i,j] = (0.25 * (a*b - e*f)).real
            m22[i,j] = (0.25 * (a*b - 2*c*d + e*f)).real
            m23[i,j] = (0.5 * ((a*d).real - (c*f).real)).real
            m24[i,j] = (0.5 * ((c*b).imag + (c*f).imag)).real
            
            m31[i,j] = (0.5 * ((a*d).real + (c*f).real)).real
            m32[i,j] = (0.5 * ((a*d).real - (c*f).real)).real
            m33[i,j] = (0.5 * ((a*f).real + c*d)).real
            m34[i,j] = (0.5 * (e*b).imag).real
            
            m41[i,j] = (0.5 * ((c*b).imag + (e*d).imag)).real
            m42[i,j] = (0.5 * ((c*b).imag + (c*f).imag)).real
            m43[i,j] = (0.5 * (e*b).imag)
            m44[i,j] = (0.5 * ((c*d) - (a*f).real)).real
 
    
    m11 = multilook(m11)
    m12 = multilook(m12)
    m13 = multilook(m13)
    m14 = multilook(m14)
    m21 = multilook(m21)
    m22 = multilook(m22)
    m23 = multilook(m23)
    m24 = multilook(m24)
    m31 = multilook(m31)
    m32 = multilook(m32)
    m33 = multilook(m33)
    m34 = multilook(m34)
    m41 = multilook(m41)
    m42 = multilook(m42)
    m43 = multilook(m43)
    m44 = multilook(m44)
    
    Shh_conj_Shh = 2 * m12 + 2 * m11 - m33 - m44
    Shv_conj_Shv = m33 + m44
    Svv_conj_Svv = 2 * m11 - 2 * m12 - m33 - m44
    Re_Svv_conj_Shh = m33 - m44
    Im_Svv_conj_Shh = 2 * m34
    Re_Shv_conj_Shh = m13 + m23
    Im_Shv_conj_Shh = m14+  m24
    
    Re_Shh = np.sqrt(Shh_conj_Shh )
    abs_Shv = np.sqrt(Shv_conj_Shv)
    abs_Svv = np.sqrt(Svv_conj_Svv)
    
    lines = Re_Shh.shape[0]
    pixels = Re_Shh.shape[1]
    
    Shh = np.zeros( (lines,pixels), dtype=np.complex64)
    Im_Shh = np.zeros( (lines,pixels), dtype=np.complex64)
    ###############
    Shh.real = Re_Shh
    Shh.imag = Im_Shh
    ###############
    
    abs_Shv_conj_Shh = np.sqrt(Re_Shv_conj_Shh * Re_Shv_conj_Shh + Im_Shv_conj_Shh * Im_Shv_conj_Shh)
    

    
    Re_Shv = abs_Shv * Re_Shv_conj_Shh / abs_Shv_conj_Shh
    Im_Shv = abs_Shv * Im_Shv_conj_Shh / abs_Shv_conj_Shh
    
    Shv = np.zeros( (lines,pixels), dtype=np.complex64)
    ###############
    Shv.real = Re_Shv
    Shv.imag = Im_Shv
    ###############
    
    abs_Svv_conj_Shh = np.sqrt(Re_Svv_conj_Shh * Re_Svv_conj_Shh + Im_Svv_conj_Shh * Im_Svv_conj_Shh)
    
    
    Re_Svv = abs_Svv * Re_Svv_conj_Shh / abs_Svv_conj_Shh
    Im_Svv = abs_Svv * Im_Svv_conj_Shh / abs_Svv_conj_Shh
    
    Svv = np.zeros( (lines,pixels), dtype=np.complex64)
    ###############
    Svv.real = Re_Svv
    Svv.imag = Im_Svv
    ###############
    
    return Shh, Shv, Svv
    


Shh, Shv, Svv = mueller(res[0],res[1],res[2],res[0].shape[0],res[0].shape[1])

np.save(os.path.join('/home/ssh_admin','Shh'),Shh)
np.save(os.path.join('/home/ssh_admin','Shv'),Shv)
np.save(os.path.join('/home/ssh_admin','Svv'),Svv)

    
def SAR_plot(HH_data,HV_data,VV_data):
    data = np.stack([HH_data,HV_data,VV_data],axis=2)
    #data = 10*np.log10(data.real**2 + data.imag**2)-83-32 #-83-32はいらないです。PALSAR2の説明書にはキャリブレーションの仕様で導入されていたんですが、入れると変になります汗
    data = 10*np.log10(data.real**2 + data.imag**2)
    amax = data.max()
    amin = data.min()
    scale = 255.0/(amax-amin)
    img = data*scale
    img = img.astype(np.uint8)
    img = Image.fromarray(img).convert('RGB')
    plt.figure(figsize=(10,10))
    img.save('SAR_01.jpg', quality=95)
    plt.imshow(img)
    
Shh = np.load('Shh.npy')
Shv = np.load('Shv.npy')
Svv = np.load('Svv.npy')

SAR_plot(Shh,Shv,Svv)