たけのこブログ

凡人が頑張って背伸びするブログ

ウェーブレット多重解像度解析(MRA)をSAR画像に適応してみた(python使用)

前回までのあらすじ

前回は、mueller行列に平均化処理をかけてノイズを除去することによって、zoomの背景画面にも使えるような合成開口レーダ(以下、SAR)画像の可視化に成功しました。

yukr.hatenablog.com

f:id:YuKR:20200330212017p:plain

きちんと前処理したSAR画像はとても綺麗ですが、いくつか弱点があります。例えば、以下のようなことが考えられてしまいます。

比較的平らな土地ってどこに分布してるの?(あるいは湖とか)

写真が諏訪湖周辺なので森林が多いこともあるけど、民家がある市街地とかってどこらへんにあるのかいまいち分からない(緑一色だから森に見えちゃう)

今回は上記の問題を解決するために、SAR画像にウェーブレット多重解像度解析(MRA)という処理をかけて低周波成分と高周波成分を取り出して、平らな土地や湖、および住宅地の分布を可視化したいと思います。

提案手法としては、以下の論文を参考にして実装しました。ミュラー平均化行列を施して抽出したSAR画像は前回の記事で生成しているので、今回はその時のSAR画像(諏訪湖周辺)をそのまま使用します(実は上下反転のままなんですが笑)。今回は諏訪湖の写真も表示します。

ci.nii.ac.jp

f:id:YuKR:20200424012250p:plain

今までのあらすじ

ちなみに、今まで取り組んできたSAR画像の生データの取り出し方や散乱特性の可視化方法などについては、以下の過去の記事を参考にしてください。

yukr.hatenablog.com

yukr.hatenablog.com

yukr.hatenablog.com

そもそも多重解像度解析とは?

多重解像度解析については真面目に解説するとそれだけでブログが一記事書けてしまうので、詳しい説明は以下を参考にしていただければと思います。超簡単に言ってしまうと、異なる解像度(周波数帯域)に分解して、SAR画像の場合は散乱特性の性質を同時に可視化できるということです。 主な使用用途としては、エッジ強調やノイズ除去に利用されることが多いです。

qiita.com

qiita.com

実際のソースコード

SAR画像のフルポラリメトリ解析に基づいて、全ての偏波成分(HH偏波、HV偏波、VV偏波)に関して多重解像度解析を行い、それぞれの結果をRGB形式で可視化することにしました。

解析に利用したソースコードは以下の通りです。

from PIL import Image
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from PIL import ImageOps
# ウェーブレット変換を行うためにPyWaveletモジュールを使いました
import pywt

def value(data):
    data = 10*np.log10(data.real**2 + data.imag**2)
    return data

def wave_plot(HH_data,HV_data,VV_data):
    data = np.stack([HH_data,HV_data,VV_data],axis=2)
    amax = data.max()
    amin = data.min()
    #scale = 255.0/(amax-amin)
    #img = data*scale
    # ウェーブレット変換の時はスケーリング処理がいらなかったです
    img = data
    img = img.astype(np.uint8)
    img = Image.fromarray(img).convert('RGB')
    print(np.max(img),np.min(img))
    plt.figure(figsize=(10,10))
    img.save('wavelet.jpg', quality=95)
    plt.imshow(img)

imgHH = value(Shh)
imgHV = value(Shv)
imgVV = value(Svv)

wp = pywt.WaveletPacket2D(data=imgHH, wavelet='db1', mode='symmetric')
hh = wp['aaa'].data # MRAを3回施した時のLL(低周波成分)を取得しています
# MRAを1回施した時のHH(高周波成分)をみたい場合はhh = wp['d'].dataで取得できます

wp = pywt.WaveletPacket2D(data=imgHV, wavelet='db1', mode='symmetric')
hv = wp['aaa'].data

wp = pywt.WaveletPacket2D(data=imgVV, wavelet='db1', mode='symmetric')
vv = wp['aaa'].data

wave_plot(hh,hv,vv)
print(np.max(hh),np.min(hh))

結果

元画像はこちらになります(再喝)

f:id:YuKR:20200330212017p:plain

では、実際に多重解像度解析を行なった結果と比較してみます。

MRAを3回施した時のLL(低周波成分)

三段階の多重解像度解析を施した後の低周波成分を取り出した結果がこちらになります。赤色の部分がHH偏波に対して多重解像度解析を施した低周波成分の結果になっていて、主に平坦な土地や湖などが抽出できていることが確認できました

元々、HH偏波のデータは比較的に地表被覆物による受信電力の差が明確です。この受信電力のみを用いて画像分類を用いてもそれなりの精度が確認できるほどなので、今回は低周波成分を取り出すことで、はっきりとした地表被覆の認識が可能になりました。

f:id:YuKR:20200424005751j:plain

MRAを1回施した時のHH(高周波成分)

f:id:YuKR:20200424010008p:plain

こちらはかなり見辛いのですが、よく見ると画像内の白いもやもやの部分が見えるかと思います。これが住宅地のような建物が確認できる地域を示していて、諏訪湖の周辺などを白いもやが覆っていることからそれが確認できるかと思います。少し見辛い結果になってしまいましたが、多重解像度解析によって建物のある地域の特定に成功しました

まとめ

今回はSAR画像に多重解像度解析を施すことによって、平坦な土地や建物が確認できる地域の可視化に成功しました。SAR画像はいろんな信号処理や画像処理の手法を適用できるので、いろんな解析を試してみたい方にはおすすめなので、興味がある方はぜひTellusにアカウントを登録して、専用の開発環境で衛星画像の解析をしていただけたらと思います!

実は先月から本当に色々なことがありまして、収入の事情で学生起業して個別事業主になったり、(全く別の分野で)一般科学誌向けにジャーナルを書いていたり、新規の研究のリサーチ(これもブログで触れてきた分野と全く異なる畑)もしていて、それに拍車をかけて請け負ってる仕事もあるので全然書く暇がないのですが、なんとか時間を作って解析をして行きたいと思いますのでよろしくお願いします(起業などに関してのエピソードは、また別記事でお話できたらと考えております)

次は差分干渉をするか、本格的に機械学習を試すか、またはSAR画像以外のデータを試して統合してみるか...汗 まだ未定ですが、何かしら面白い解析ができたらと思います。