日本は、地震、台風、豪雨などといった自然災害が発生しやすい国です。近年、こういった災害による被害が広域化・同時多発化する傾向が見られ、今後も、南海トラフ地震のような巨大地震や豪雨災害など甚大な災害の発生が懸念されています。災害の初動対応では、被害状況をいち早く把握することが必要不可欠ですが、被害が広域化・同時多発化することで対応がより難しくなります。このような中、衛星画像を用いることでより早く、より効率的に概況調査を実施できると期待されています。特に、この章で扱う合成開口レーダ(SAR)は全天候型(観測時の天候に左右されにくい)であり、夜間でも観測可能であるため、災害発生時に迅速な対応ができます。
SARでは、受信信号の強度と位相を記録しています。各観測での位相差を計算し、地表面の変位を計測するのがこの章の前半で扱う差分干渉解析です。差分干渉解析によって、地震や火山活動に伴う地殻変動や地すべり、工事に伴う地盤沈下といった地表面変位の状況を把握することができます。このうち、この章では地震による地殻変動を対象とします。
一方、複数時期の観測データを利用して受信信号の強度変化を調べることで、地表面の様子が著しく変化した箇所を把握することができます。これがこの章の後半で扱う内容で、土砂災害や浸水箇所の概況調査が可能です。
この章ではJAXAの陸域観測技術衛星2号「だいち2号」(ALOS-2)に搭載された合成開口レーダPALSAR-2のプロダクトのうち、Level 1.1(L1.1)とLevel 2.1(L2.1)を使用します。L1.1は複素データ(実部、虚部のデータを格納している)のため、受信信号の強度に加え、位相情報も保持しています。差分干渉解析には位相情報が必要であるため、この章の前半ではL1.1を使用します。
L1.1から複素数の大きさを計算して実画像とした上で、実際の地表面上の位置に投影されたものがL1.5、そこからさらに地形による倒れこみを補正(オルソ補正)したものがL2.1です。この章の後半では受信信号の強度を利用した解析を行うため、そのまま地図投影ができるL2.1を使用します。ただし、複素データではないので位相情報が失われている、つまり、L1.5やL2.1を差分干渉解析に用いることはできないのでご注意ください。
TelluSARを用いた干渉SARの活用方法を学びます。干渉SARは地表面の変化を調べる解析手法で、主に地殻変動や地盤沈下、地すべりなどによる地表面変位の状況を把握するために使われています。干渉SARを初めて聞いたという方は、宙畑に干渉SARの概要や解析事例を紹介した記事があるので、そちらを参照してください。⇒「干渉SAR(InSAR)とは-分かること、事例、仕組み、読み解き方-」
さて、実際に干渉SARを行ってみましょう。使用する衛星データはALOS-2/PALSAR-2のL1.1(位相情報を含む)です。Tellusに搭載されているPALSAR-2オリジナルデータを用いて1つずつ処理を行う方法は、宙畑の「【コード付き】地盤の沈降が分かる干渉SAR解析をTellusでやってみた」で説明されています。が、複雑な処理なので、衛星データに慣れてない方にはハードルが高いかもしれません。そこで、TelluSARの差分干渉解析機能を使います。TelluSARではシーン(画像)を選んでしまえば自動で処理を実行してくれるため、誰でも簡単に干渉SAR解析を行うことができます。TelluSARの基本的な使用方法は宙畑の「ついに登場!公式解析ツール「TelluSAR(てるーさー)」 APIの使い方を徹底解説!」で紹介されているので参考にしつつ、ここではより実践的に、平成30年北海道胆振東部地震を対象にTelluSARの差分干渉解析機能を用いて差分干渉画像を作成し、地殻変動の状況を確認します。
平成30年9月6日未明に、北海道胆振東部地方中東部を震源とする地震が発生しました。最大震度は厚真町で観測された震度7で、厚真町を中心に斜面崩壊が多発、札幌市清田区など市街地での液状化現象や北海道全域での大規模停電が発生しました。
それでは、使用するシーン(画像)を選んで解析を始めましょう。利用できるシーンはTelluSARの検索機能で確認できますが、まずは、TelluSAR(API)を使うための事前準備をします。
必要なモジュール、ライブラリをインポートします。
import requests, datetime, json, pprint
from skimage import io
from io import BytesIO
マーケットトークンを発行します。 ※get_market_tokenを使用しています。
TOKEN = '自分のトークンをはる'
provider_id = 'tellus-product'
tool_label = 'tellusar-api-f'
# マーケットトークン有効期限(30分に設定)
# デフォルトでは5分、最長で60分まで設定可能
expires_at = (datetime.datetime.now(datetime.timezone(datetime.timedelta(hours=+9))) + datetime.timedelta(minutes=30)).isoformat()
def get_market_token(payload={}):
url = 'https://sdk.tellusxdp.com/api/manager/v1/auth/api_access_token/token'
headers = {
'Authorization': 'Bearer ' + TOKEN
}
r = requests.post(url, headers=headers, data=json.dumps(payload))
if r.status_code is not 200:
raise ValueError('status error({}).'.format(r.status_code))
return json.loads(r.content)
# マーケットトークンを発行する
ret = get_market_token({'provider_id': provider_id, 'tool_label': tool_label, 'expires_at':expires_at})
market_token = ret['token']
print(market_token)
差分干渉解析では2時期に撮影されたSAR画像のペアを用いて、その間に生じた地表面変位を解析するので、基準となるシーン(primary画像)とその後に観測されたシーン(secondary画像)が必要となります。今回は地震の影響を把握することが目的なので、primary画像は地震前から、secondary画像は地震後から選ぶことになります。TelluSARで利用できるシーンを検索する際にパラメータとして位置情報を与えておくと、その範囲を含むシーンを確認することができるので、震源の位置を与えておきます。
# 解析に用いるシーンが次の範囲を含む
# 平成30年北海道胆振東部地震の震源のおおよその位置
parameter1 = {
'left_bottom_lon' : 142.000,
'left_bottom_lat' : 42.690,
'right_top_lon' : 142.000,
'right_top_lat' : 42.690
}
解析に用いるシーンの位置条件を設定したら、その範囲を含むシーンの一覧を確認します。一覧ではシーンの観測範囲や偏波の情報も得られますが、情報量が膨大になるので、ここでは観測時刻とシーンIDのみを表示します。
※get_free_sceneを使用しています。
# 利用できるシーン一覧を取得する
def get_free_scene(market_token, payload={}):
url = 'https://tellusar.tellusxdp.com/api/v1/search'
headers = {
'Authorization': 'Bearer ' + market_token
}
r = requests.get(url, headers=headers, params=payload)
if r.status_code is not 200:
print(r.status_code)
return json.loads(r.content)
get_free_scene_result = get_free_scene(market_token, parameter1)
# シーンの枚数を表示
print(len(get_free_scene_result['data']['scenes']))
# シーンの観測時刻とシーンIDのみを表示
scene_list = []
for dict in get_free_scene_result['data']['scenes']:
newdict = {key:dict[key] for key in ['scene_id','observation_datetime']}
scene_list.append(newdict.copy())
pprint.pprint(scene_list)
震源の位置を含む観測シーンが73枚あることがわかります。まずは、primary画像を選んでみましょう。地震は2018年9月6日に発生しています。一覧の'observation_datetime'を見ていくと、地震の直前では、2018年8月23日に観測されたシーンがあるようなので、そのシーンを選ぶことにします。(時刻はUTC表記なので、日本時刻が知りたい場合は+9時間してください。)
# 2018年8月23日のシーンID
scene_id = 'ALOS2229472760-180823'
primary画像を選んだので、次はsecondary画像を選びます。secondary画像を選ぶ際に必要な条件は、以下の通りです。条件を満たしていないシーンを選ぶと差分干渉解析を実施することができません。
TelluSARでは基準となるシーンの'scene_id'を指定するとこれらの条件を満たすsecondary画像の候補を表示してくれます。
※get_scene_afterを使用しています。
# ペアの候補となるシーン一覧を取得する
# after_scene_id
def get_scene_after(scene_id, market_token, payload={}):
url = 'https://tellusar.tellusxdp.com/api/v1/search/{}/afters'.format(scene_id)
headers = {
'Authorization': 'Bearer ' + market_token
}
r = requests.get(url, headers=headers, params=payload)
if r.status_code is not 200:
print(r.status_code)
return json.loads(r.content)
get_scene_after_result = get_scene_after(scene_id, market_token, {})
print(len(get_scene_after_result['data']['scenes']))
scene_after_list = []
for dict in get_scene_after_result['data']['scenes']:
newdict = {key:dict[key] for key in ['scene_id','observation_datetime']}
scene_after_list.append(newdict.copy())
pprint.pprint(scene_after_list)
9枚のシーンが候補として表示されました。この中からsecondary画像を選ぶ時のポイントは「観測日」です。地表面の被覆状態の変化は干渉性の良し悪し(結果の良し悪し)に影響するため、primary画像に対して、観測日が近い、または、同一季節のシーンを用いることが望ましいです。地震が発生した日の正午頃に観測されたシーンがあるため、そのシーンを選ぶことにします。
今回は「観測日」に注目してシーンを選定しましたが、他にも、衛星の観測条件を考慮することがあります。差分干渉解析では、衛星が軌道上の同一の箇所に来た時に観測したシーンのペアを用いていますが、衛星の位置が完全に一致することはありません。このときの軌道間距離の衛星視線方向に直交する成分を垂直基線長(Bperp)と呼び、このBperpが大きいと干渉性が低下します(結果が悪くなります)。BperpはTellus上で確認することはできませんが、JAXAが運用しているAUIG2(ALOS及びALOS-2の観測データに関する各種オンラインサービスを提供するシステム)で調べることができます。
# 2枚目のscene_id
# 2018年9月6日(地震直後に観測)のシーンID
scene_after_id = 'ALOS2231542760-180906'
用いるシーンのペアを選ぶことができたので、差分干渉解析を実施し、差分干渉画像を作成します。先ほど紹介した宙畑の記事から、差分干渉解析のフローを以下に示しますが、TelluSARではこれらの処理を自動で実施するので、誰でも簡単に差分干渉画像を得ることができます。
干渉SARの処理フロー( https://sorabatake.jp/12465/ )
処理パラメータとしてルック数とフィルター回数を指定できます。SAR画像にはスペックルノイズと言われるごま塩状ノイズが生じますが、このノイズを低減する方法としてルック処理というものがあります。ルック処理はSARに特有な処理方法であり、厳密には少し違いますが、複数のピクセルの平均値を代表値としてダウンサンプリングする方法だと思ってください。平均値を計算する画素数を指定することができ、'nlook_rg'がレンジ方向(電波照射方向)、'nlook_az'がアジマス方向(衛星進行方向)のルック数を示しています。また、生成された干渉縞について、移動平均フィルターもノイズ低減手法として有効であり、'filter'でその回数を指定することができます。
ルック処理のイメージ(厳密にはもう少し複雑な処理をしています。)
ルック数の違いによる解像度の変化(ルック数のより大きい右側の画像の方が解像度が低く滑らかになります。Original data provided by JAXA)
できるだけノイズの影響を低減してきれいな差分干渉画像を得るため、ルック数とフィルター回数は指定できる最大値に設定しています。今回の解析対象は地震による地殻変動であり、空間スケールが比較的大きい事象です。この場合は、ノイズ低減処理(ルック処理・フィルター)を強めにかけても問題ありませんが、注目している事象が局所的な事象である場合、これらの処理を強めにかけると求めたい変位が埋もれて見えなくなってしまう可能性があります。そのような場合はルック数を下げるか、フィルター回数を下げて調整してください。
# 処理パラメータ
# ノイズをできるだけ低減するため、ルック数とフィルター回数を最大に設定
parameter2 = {
'before_scene_id': scene_id,
'after_scene_id': scene_after_id,
'polarisation': 'HH',
'nlook_rg': 5,
'nlook_az': 7,
'filter': 2,
}
処理パラメータを設定したら、処理を実行します。
※request_workを使用しています。
# 処理を実行する
# work_idが返却される
def request_work(market_token, payload={}):
url = 'https://tellusar.tellusxdp.com/api/v1/works'
headers = {
'Authorization': 'Bearer ' + market_token
}
r = requests.post(url, headers=headers, params=payload)
print(r.url)
if r.status_code is not 200:
print(r.status_code)
return json.loads(r.content)
work_result = request_work(market_token, parameter2)
print(work_result)
処理を開始すると処理Id(work_id)が返されます。既に処理された結果が存在する場合は'exist_flag': Trueとなります。処理時間は30分~数時間程度かかります。
処理ステータスの確認は、work_idを指定する方法と一覧で確認する方法があります。
※get_work, get_work_listを使用しています。
# work_idを指定して処理結果を取得する
def get_work(work_id, market_token, payload={}):
url = 'https://tellusar.tellusxdp.com/api/v1/works/{}'.format(work_id)
headers = {
'Authorization': 'Bearer ' + market_token
}
r = requests.get(url, headers=headers, params=payload)
print(r.url)
if r.status_code is not 200:
print(r.status_code)
return json.loads(r.content)
work_id = work_result['data']['work_id']
get_work_result = get_work(work_id, market_token, {})
print(get_work_result)
# 処理結果の一覧を取得する
def get_work_list(market_token, payload={}):
url = 'https://tellusar.tellusxdp.com/api/v1/works'
headers = {
'Authorization': 'Bearer ' + market_token
}
r = requests.get(url, headers=headers, params=payload)
print(r.url)
if r.status_code is not 200:
print(r.status_code)
return json.loads(r.content)
get_work_result_list = get_work_list(market_token, {})
print(get_work_result_list)
作成された差分干渉画像を確認してみます。タイル画像(png)を取得して表示する方法とGeoTIFFを表示する方法がありますが、タイル画像の表示は宙畑で紹介されているので、ここではGeoTIFFを表示します。
url = 'https://tellusar.tellusxdp.com/api/v1/works/{}/tifs/fringe_diff.tif'.format(work_id)
headers = {
'Authorization': 'Bearer ' + market_token
}
r = requests.get(url, headers=headers)
print(r.url)
if r.status_code is not 200:
print(r.status_code)
fringe=io.imread(BytesIO(r.content))
io.imshow(fringe)
差分干渉画像(Original data provided by JAXA)
差分干渉画像(GeoTIFF)をダウンロードします。位相に応じた色付けがされており、地理情報を持っているので、GIS上にそのまま表示することができます。
with open('fringe.tif', 'wb') as f:
for data in r.iter_content(1024):
f.write(data)
ダウンロードした差分干渉画像をGIS上に表示してみました。
画像の下2/3は海なので差分干渉解析結果が得られず砂嵐状になっています。上1/3は陸域で干渉性が比較的良く、震源の南東に目玉状の変状が見られています。
差分干渉解析結果から地表面変位量を読み解きます。変位を読む際は衛星側から読むのが一般的ですが、変状箇所が目玉状に集まっている場合はその外側から読んでいくこともあります。目玉状の変状の外側、陸域の大部分は水色を呈しています。目玉の外(A)から目玉の中心(B)に向かうにつれて、水色→黄色→ピンクのように色変化しています。この色変化を凡例に照らし合わせると、黄色→ピンクの色変化は-πを超えてしまいますが、これは、得られる位相が-π~+πに折り畳まれているためなので、凡例をもう1周分追加して読み取ってください。おおよそπ分強、衛星に近づいていると読み取れます。この位相量は利用するレーダ波の波長に対応しており、波長の半分が2πに相当します(衛星と地表面間の往復分なので、波長の半分となります)。ALOS-2の場合、波長は約23.6cmなので、2π=11.8cmとして換算でき、地点Aに対して地点Bは6cm強、衛星に近づいていると推定できます。このように、差分干渉解析で得られる地表面変位は相対変位であるので、変位量を読み取る際は基準となる点を設定する必要があります。
差分干渉解析で得られた変位量は衛星視線方向の変位量です。SARによる観測は、衛星から斜め下方向に電波を照射しているため、得られる変位量は上下成分と東西成分が合成されたものです。ここで、今回用いたシーンの観測条件を確認します。
※search_fileに変更を加えています。
def search_file_id(scene_id, params={}, next_url=''):
if len(next_url) > 0:
url = next_url
else:
url = 'https://file.tellusxdp.com/api/v1/origin/search/palsar2-l11/{}'.format(scene_id)
headers = {
'Authorization': 'Bearer ' + TOKEN
}
r = requests.get(url, params=params, headers=headers)
if not r.status_code == requests.codes.ok:
r.raise_for_status()
return r.json()
ret = search_file_id(scene_id, {})
pprint(ret)
衛星進行方向はディセンディング('orbit_direction': 'D')で、右側観測('observation_direction': 'R')なので、得られた差分干渉画像と衛星の位置関係は次の通りです。
得られた結果は東側からの観測結果であり、衛星に近づく変位が出ていた場合、隆起・東向き変位が生じている可能性があります。ここから隆起量あるいは東向き変位を求めたい場合はどちらか一方の変位がないと仮定をおく必要があります。このほか、反対側からの観測(アセンディング・右側観測)シーンのペアについて差分干渉解析を実施することで西側からの観測結果を得ることができ、東側からの観測結果と合わせることで、準上下方向と準東西方向に分離することができます。これを2.5次元解析と呼びます。同時性が保てないことが課題ではありますが、胆振東部地震の場合は運良く、primaryとsecondaryともに、ディセンディングの観測(今回用いたペア)の半日後の観測があったため、それらのシーン(アセンディングのペア)を用いた解析を合わせると、今回得られた変状は隆起を示しており、震源の南東側が6cm強の隆起を示していると判断できたようです。(国土地理院の解析例)
TelluSARの出力結果には、強度画像(primary及びsecondary)、差分干渉画像に加えて、コヒーレンス画像があります。コヒーレンス(干渉性)とは、位相情報の推定の良し悪しの指標であり、コヒーレンスが高いほど信頼性が高いと考えます。
コヒーレンス画像を見てみると、陸域が比較的高い値を示し(明るい)、海では低い値(暗い)となることがわかります。差分干渉画像で海は砂嵐状になっていましたが、このように位相がランダムとなる(位相の推定ができていない、信頼性が低い)箇所ではコヒーレンスが低い値となります。逆に言えば、コヒーレンスが低い箇所の位相はランダムなので地表面変位を読み解くことが難しく、差分干渉画像からコヒーレンスの低い箇所をマスクしておくと解釈しやすくなることがあります。
コヒーレンス画像(Original data provided by JAXA)
コヒーレンスが低くなる要因として地表面の状態の著しい変化があることが考えられます。この特徴を利用した被害箇所の解析方法として、この章の後半で扱う2時期のSAR強度画像を用いたRGB疑似カラー合成画像にコヒーレンスの情報を加えて変位箇所を可視化するMulti Temporal Coherence mapping(MTC)やコヒーレンスの変化に着目した3時期コヒーレンス解析といったものがあります。
これまでInSAR界隈では、基準となる画像を"master"、2枚目の画像を"slave"と呼ぶことが主流でした。おそらく、皆さんがInSARの解析事例をネット等で調べると、この用語を使っているものが多く出てくるのではないでしょうか。
昨今、"Black Lives Matter"が話題となっていますが、この流れがInSARの分野にまで波及し、"master" を "primary" に、 "slave" を "secondary" に言い換えるようにしましょうという動きが出てきています。本稿ではこの考えに賛同して、"primary"と"secondary"という用語で説明するようにしています。
ここからはSARの強度画像を用いた解析方法を学びます。強度画像は地表面上の散乱体の有無や凹凸を反映したものであり、散乱体が有る場合や凹凸が粗い場合により強い後方散乱が得られるため、大きい値をとります。土砂災害や浸水の被害範囲では地表面の様子が著しく変化するため強度画像に変化が生じ、災害発生時を挿む二時期の強度画像の差分をとることで被害範囲を可視化することが期待できます。
Tellusには、ALOS-2/PALSAR-2(L2.1)のオリジナルデータが搭載されているので、今回はそのデータを使います。オリジナルデータの取得方法は宙畑の「【ゼロからのTellusの使い方】提供元オリジナルデータを取得する」にて解説されているので参考にしつつ、平成30年7月豪雨時の岡山県倉敷市真備町の浸水を対象に、二時期の強度画像を用いた浸水範囲の推定を行います。
平成30年6月28日から7月8日にかけて、梅雨前線や台風の影響により、西日本を中心に全国的に広い範囲で豪雨が発生しました。西日本から東海地方の多くの地点で24、48、72時間雨量の観測史上最大値が更新され、西日本の広域にわたって土砂災害や洪水、浸水被害が発生しました。
Tellus上で、岡山県倉敷市真備町周辺のポリゴン(GeoJSON)を用意します。ポリゴンの作成方法は漁業分野を参照してください。
作成したポリゴン(GeoJSON)をTellus開発環境にアップロードした後、読み込んで範囲の情報を確認します。
import geopandas as gpd
aoi = gpd.read_file('aoi.geojson', encoding='utf-8')
bounds = aoi['geometry'].bounds
bounds
使用するシーンの条件を設定します。真備町の浸水被害は2018年7月7日に発生しているので、第一に、この直後に観測されたシーンが必要となります。次に、今回は二時期のデータを用いて浸水被害を推定するため、浸水前のシーンも必要となります。差分干渉解析のときと同様に、ここでも観測日が重要となります。時期が離れると人工改変等の影響が出てしまう、あるいは、季節変化の影響が出ることもあるので、できるだけ期間の短い、あるいは、季節を揃えたシーンを選択することが望ましいです。
上記の条件に合うシーンを検索します。期間を1年間として、関心領域を含むようにしたときのパラメータは次のように指定できます。
#使用するシーンの条件
parameter2 = {
'after' : '2017-06-01T00:00:00',
'before' : '2018-07-10T00:00:00',
'left_bottom_lon' : bounds['minx'],
'left_bottom_lat' : bounds['miny'],
'right_top_lon' : bounds['maxx'],
'right_top_lat' : bounds['maxy']
}
使用するシーンの条件を設定したら、Tellusに搭載されているPALSAR-2オリジナルデータを検索します。受信信号の位相情報は必要ないので、地図投影されたL2.1データを使用します。
※search_fileを使用しています。
import requests
from pprint import pprint
TOKEN = '自分のトークンをはる'
provider_id = 'tellus-product'
tool_label = 'tellusar-api-f'
def search_file(params={}, next_url=''):
if len(next_url) > 0:
url = next_url
else:
url = 'https://file.tellusxdp.com/api/v1/origin/search/palsar2-l21'
headers = {
'Authorization': 'Bearer ' + TOKEN
}
r = requests.get(url, params=params, headers=headers)
if not r.status_code == requests.codes.ok:
r.raise_for_status()
return r.json()
ret = search_file(parameter2)
pprint(ret['count'])
pprint(ret['items'])
比較する2枚の画像を選ぶ時の条件は差分干渉解析のときと同じです。ここでおさらいしておきましょう。
SAR観測は偏波によって感度が異なります(同じ対象でも受信信号の強度が変化する)。また、観測する地形と衛星の位置関係によって観測できない範囲が生じたり、歪みが生じます。観測条件を揃えない場合、これらの要因がノイズとして現れるようになるので、観測条件は同一である必要があります。
さて、設定した期間に2枚の観測シーンがあることが確認できました。運良く条件(観測条件が全て同一であること)にもあてはまっています。2枚目の方は、日本時刻で7月7日の24時ごろに観測されているので、災害直後の画像として適しています。1枚目の方は、2枚目の観測と同じ条件で観測されている、かつ、期間も3ヶ月前と比較的短いので、災害前のシーンとして採用することにします。
使用するシーンが決まったら、シーンId(dataset_id)を指定してオリジナルデータをダウンロードします。
※publish_file, download_fileを使用しています。
import requests
from pprint import pprint
def publish_files(dataset_id):
url = 'https://file.tellusxdp.com/api/v1/origin/publish/palsar2-l21/{}'.format(dataset_id)
headers = {
'Authorization': 'Bearer ' + TOKEN
}
r = requests.get(url, headers=headers)
if not r.status_code == requests.codes.ok:
r.raise_for_status()
return r.json()
published1 = publish_files('ALOS2210180680-180414')
published2 = publish_files('ALOS2222600680-180707')
pprint(published1)
pprint(published2)
def download_file(url):
headers = {
'Authorization': 'Bearer ' + TOKEN
}
r = requests.get(url, headers=headers)
if not r.status_code == requests.codes.ok:
r.raise_for_status()
return r.content
def save_file(data, name):
with open(name, "wb") as f:
f.write(data)
file1 = download_file(published1['files'][1]['url'])
save_file(file1, published1['files'][1]['file_name'])
file2 = download_file(published2['files'][1]['url'])
save_file(file2, published2['files'][1]['file_name'])
print('done')
これで、オリジナルデータを手に入れることができました。このオリジナルデータを用いて解析をしていきます。
※オリジナルデータのTellus外への持ち出しは禁止されているので気をつけてください!
2時期の強度画像の差分の表現方法は大きく分けて2つあります。方法①は2時期の強度画像にそれぞれ別の色付けをし、重ね合わせることで変化箇所を可視化する方法です。方法②は2時期の強度画像の画素値の差を直接計算する方法です。
ここでは、これら2つの方法について両方学べるようになっていますが、方法②はやや中級者向けの内容になるので、まずは方法①の内容をきちんとおさえておきましょう。
2時期の強度画像を用いて変化した箇所をわかりやすく表現する方法として、RGB疑似カラー合成画像を作成し、変化箇所を可視化します。災害前の画像をR(赤)に、災害後の画像をG(緑)とB(青)に割り当てることで、災害前後で変化が大きい箇所は赤色、または、青色(シアン色)に、変化が小さい箇所は白黒として表現できます。
オリジナルデータ(GeoTIFF)を読み込み、処理時間節約のために関心領域で画像をクリップ、RGB疑似カラー合成画像を作成し、浸水範囲を推定します。ここではGDALを使いますが、処理は基本的に基礎分野と同じなので、詳細は基礎分野を参照してください。
必要なモジュール、ライブラリをインポートします。
from osgeo import gdal, gdalconst
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
オリジナルデータを読み込み、座標系を確認します。座標系の説明も基礎分野を参照してください。
img_before_original = gdal.Open('IMG-HH-ALOS2210180680-180414-UBSR2.1GUA.tif', gdalconst.GA_ReadOnly)
img_after_original = gdal.Open('IMG-HH-ALOS2222600680-180707-UBSR2.1GUA.tif', gdalconst.GA_ReadOnly)
src_prj=img_before_origina.GetProjection()
src_info=img_before_origina.GetGeoTransform()
print(src_prj)
print(src_info)
座標系はUTMであることが確認できます。今回は関心領域の範囲を緯度経度で設定してあるので、解析する画像の座標系も緯度経度に変換して画像をクリップします。
# 画像の座標系を変換、関心領域で画像をクリップ
gdal.Warp('ALOS2_180414_clip.tif', img_before_original, dstSRS='EPSG:4326', outputBounds=(bounds['minx'], bounds['miny'], bounds['maxx'], bounds['maxy']))
gdal.Warp('ALOS2_180707_clip.tif', img_after_original, dstSRS='EPSG:4326', outputBounds=(bounds['minx'], bounds['miny'], bounds['maxx'], bounds['maxy']))
print("done")
クリップした画像を表示して確認します。その前に、画素値のヒストグラムを確認し、表示するレンジを決定します。ヒストグラムに偏りが見られる場合は表示するレンジを適切に設定しておかないと、表示された画像が暗すぎる、もしくは逆に、明るすぎて見にくいことがあります。
img_before = gdal.Open('ALOS2_180414_clip.tif', gdalconst.GA_ReadOnly)
img_after = gdal.Open('ALOS2_180707_clip.tif', gdalconst.GA_ReadOnly)
array_before = img_before.ReadAsArray()
array_after = img_after.ReadAsArray()
plt.figure(figsize=(20,5))
plt.subplot(1,2,1)
plt.hist(array_before.flatten(),bins=300,range=(1, 65545))
plt.title('before')
plt.subplot(1,2,2)
plt.hist(array_after.flatten(),bins=300,range=(1, 65545))
plt.title('after')
plt.show()
ヒストグラムを確認すると画素値はだいたい0-10000に集まっているので、表示するレンジを0-10000に設定します。
plt.figure(figsize=(10,10))
plt.subplot(1,2,1)
plt.imshow(array_before,vmax=10000,cmap='gray')
plt.title('before')
plt.subplot(1,2,2)
plt.imshow(array_after,vmax=10000,cmap='gray')
plt.title('after')
plt.show
Original data provided by JAXA
災害前後の画像を表示してみました。SARでは反射(後方散乱)した電波の強度を記録しているため、光学画像(カラー)とは違って白黒の画像になります。明るいところが強い反射が得られている箇所、暗いところが弱い反射かあるいは反射が得られていない箇所を示しています。災害前後を比較してみると、災害後(浸水後)の画像では暗い範囲が増えているのがわかります。浸水時は水面で電波が鏡面反射してしまうため後方散乱が得られにくく、災害前後で暗く変化した箇所は水面になった(浸水した)可能性があります。
レーダ波が水面で鏡面反射する様子
災害前後の2枚の画像からRGB疑似カラー合成画像(R:災害前、G:災害後、B:災害後)を作成したとき、土砂崩れ等で散乱体となる樹木が消失した場合や浸水によって受信信号が得られにくい(強度が小さくなる)ときは赤色に、反対に、土砂崩れ等の崩壊土や浸水が引いた後に漂流物が堆積した場合はその堆積物が散乱体となるので強度が大きくなり、青色(正しくは、緑+青なのでシアン色)を呈します。災害前後で変化が見られない箇所は白黒として表示されます。
RGBカラーなので、画像のデータ型はByte型にしておいた方が扱いやすいです。変換する際は、表示した時のレンジ0-10000をスケールパラメータに指定しておくと、先ほど表示した画像の濃淡のままの画像ができます。
# スケールを指定して画像をByte型に変換
gdal.Translate('ALOS2_180414_8bit.tif', img_before, outputType=gdal.GDT_Byte, scaleParams=[[0,10000]])
gdal.Translate('ALOS2_180707_8bit.tif', img_after, outputType=gdal.GDT_Byte, scaleParams=[[0,10000]])
print("done")
画像をByte型に変換したら、画像を結合してRGB疑似カラー合成画像を作成します。
img_before_8bit = gdal.Open('ALOS2_180414_8bit.tif', gdalconst.GA_ReadOnly)
img_after_8bit = gdal.Open('ALOS2_180707_8bit.tif', gdalconst.GA_ReadOnly)
array_before_8bit = img_before_8bit.ReadAsArray()
array_after_8bit = img_after_8bit.ReadAsArray()
Xsize=img_before_8bit.RasterXSize
Ysize=img_before_8bit.RasterYSize
dtype=gdal.GDT_Byte
band=3
outfile = 'RGBcomposite.tif'
outRGB= gdal.GetDriverByName('GTiff').Create(outfile, Xsize, Ysize, band, dtype)
outRGB.SetProjection(img_before_8bit.GetProjection())
outRGB.SetGeoTransform(img_before_8bit.GetGeoTransform())
outRGB.GetRasterBand(1).WriteArray(array_before_8bit)
outRGB.GetRasterBand(2).WriteArray(array_after_8bit)
outRGB.GetRasterBand(3).WriteArray(array_after_8bit)
outRGB.FlushCache()
plt.figure(figsize=(8,8))
image = mpimg.imread(outfile)
plt.imshow(image)
plt.title("RGB composite")
plt.show
Original data provided by JAXA
RGB擬似カラー合成画像を表示すると赤く色づけられた箇所が目立ちます。SAR画像に見られる変状が何によるものか判読しやすくするために、地理院タイルを重ねてみましょう。川が赤くなっているのは川の増水を反映していると考えられ、画像の中央あたりに見られる赤色の範囲が浸水をしている範囲と推定されます。
国土地理院が公開している推定浸水範囲(https://www.gsi.go.jp/BOUSAI/H30.taihuu7gou.html)と重ね合わせてみると赤色箇所と概ね一致しているように見られます。このように、RGB疑似カラー合成画像は変化箇所(2時期の画像で差分が大きい箇所)を色付けして強調するため、被害範囲を把握する際に便利な表現方法です。
2つ目の方法は2時期の画像の画素値から差分を計算し、閾値を決めて変化箇所と非変化箇所に分け、浸水範囲を抽出します。このように閾値に応じて分類する方法をレベルスライス法といいます。
GDALを用いて画像から画素値の配列を作成し、計算していきます。方法①でクリップした画像から始めますが、画素値を計算する場合は、方法①の可視化に比べてノイズの影響を受けやすくなるので、フィルターをかけてノイズを低減します。フィルターの種類はノイズが低減できれば特に決まりはありませんが、ここではメディアンフィルターを使うことにします。メディアンフィルターは外れ値のようなスパイク状のノイズを補正しやすいという特徴があります。
img_before = gdal.Open('ALOS2_180414_clip.tif', gdalconst.GA_ReadOnly)
img_after = gdal.Open('ALOS2_180707_clip.tif', gdalconst.GA_ReadOnly)
array_before = img_before.ReadAsArray()
array_after = img_after.ReadAsArray()
# フィルター処理
def median_filter(img_array, wsize):
d = int((wsize-1)/2)
flt = img_array.copy()
height = img_array.shape[0]
width = img_array.shape[1]
# ウィンドウ内のピクセル値の中央値を出力
# 画像端は元画像の値を保持
for y in range(d, height - d):
for x in range(d, width - d):
flt[y][x] = np.median(img_array[y-d:y+d+1, x-d:x+d+1])
return flt
# 3×3のメディアンフィルターを適用
array_before_flt = median_filter(array_before, wsize=3)
array_after_flt = median_filter(array_after, wsize=3)
# 災害前の画像をフィルター適用前後で表示
# 画像の一部分を表示
plt.figure(figsize=(10,10))
plt.subplot(1,2,1)
plt.imshow(array_before[1300:1800, 2000:3000],vmax=10000,cmap='gray')
plt.title('original')
plt.subplot(1,2,2)
plt.imshow(array_before_flt[1300:1800, 2000:3000],vmax=10000,cmap='gray')
plt.title('filtered')
plt.show
Original data provided by JAXA
フィルター適用後の画像(右図)のノイズが低減し、きれいになったことが確認できます。このフィルター後の画像の画素値を用いて、2時期の差分を計算します。画像のデータ型はUInt16ですが、差分を計算すると負の値となる可能性があるので、データ型をInt32に変換してから引き算します。その後、閾値を決めて浸水の可能性がある変化箇所とそうでない箇所に分類します。浸水の可能性がある場合は強度は小さくなっているはずなので、災害前後の画素値を比較すると災害後の方が画素値が小さく、差分は正の値になります。画素値は0-10000に固まっていたのでその10%として、災害前後で画素値に1000以上の差があれば変化箇所と分類することにします。
# 差分の計算と閾値による分類
def create_geotiff_from_array(outfile, array, img):
Xsize=img.RasterXSize
Ysize=img.RasterYSize
dtype=gdal.GDT_Byte
band=1
outTiff= gdal.GetDriverByName('GTiff').Create(outfile, Xsize, Ysize, band, dtype)
outTiff.SetProjection(img.GetProjection())
outTiff.SetGeoTransform(img.GetGeoTransform())
outTiff.GetRasterBand(1).WriteArray(array)
outTiff.GetRasterBand(1).SetNoDataValue(0) # 0をデータなしに設定
outTiff.FlushCache()
print("done")
# データ型変換(UInt16→Int32)
array_before_flt2 = array_before_flt.astype(np.int32)
array_after_flt2 = array_after_flt.astype(np.int32)
# 2時期の差分を計算
array_diff = array_before_flt2 - array_after_flt2
array_diff2 = array_diff.copy()
# 差分が1000以上の場合は変化あり(1)、1000未満の場合は変化なし(0)
array_diff2[array_diff2 < 1000] = 0
array_diff2[array_diff2 >= 1000] = 1
# 元の画像(クリップ後)に合わせてGeoTIFFの作成
create_geotiff_from_array('diff2.tif', array_diff2, img_before)
差分の閾値を1000として変化箇所と非変化箇所の分類をしました。変化箇所がある程度固まっている範囲もありますが、ごま塩状に散らばっている変化箇所も見られます。ここから浸水の可能性が高い箇所を抽出していきます。
第一に浸水の可能性がない範囲(海・河川といった元々水域や山地など)をマスクします。ここでは簡単に、洪水浸水想定区域を利用し、被害が想定される区域で実際に被害がありそうな範囲はどこか衛星データから解析するという視点で進めます。
洪水浸水想定区域は国土数値情報(https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-A31-v2_1.html)からダウンロードしたポリゴンデータ(Shapefile)を利用します。このポリゴンデータを画像化(ラスター化)し、配列に変換してマスクする範囲を指定します。ポリゴンデータのラスター化はgdal.Raterize()を使用します。
# ポリゴンのラスター化
# 元の画像から画素の解像度を取得
img_prj=img_before.GetGeoTransform()
xPixelSize = img_prj[1] # x方向の解像度
yPixelSize = -img_prj[5] # y方向の解像度
# ラスター化
mask=gdal.Rasterize('mask.tif', # 出力ラスター
'A31-12_33.shp', # 入力ポリゴン
outputSRS='EPSG:4326', # 出力ラスターの座標系
xRes=xPixelSize, yRes=yPixelSize, # 出力ラスターの解像度
outputBounds=[bounds['minx'], bounds['miny'], bounds['maxx'], bounds['maxy']], # 出力ラスターの範囲
burnValues=1.0, # 出力ラスターの画素値
outputType=gdal.GDT_Byte # 出力ラスターのデータ型
)
mask=None
print("done")
白い範囲(1)が洪水浸水想定区域です。その範囲にある変化箇所だけを残して、それ以外の黒い範囲(0)をマスクします。
# マスク処理
mask = gdal.Open('mask.tif', gdalconst.GA_ReadOnly)
array_diff_mask = array_diff2.copy()
array_mask = mask.ReadAsArray()
# マスク画像で画素値が'0'となる箇所をマスク
array_diff_mask[array_mask == 0] = 0
# 元の画像(クリップ後)に合わせてGeoTiffの作成
create_geotiff_from_array('diff_mask.tif', array_diff_mask, img_before)
マスク処理の後は、ごま塩状に存在する小さい変化箇所は信頼性に乏しくノイズである可能性が高いので、空間的にある程度の大きさをもつ変化箇所だけを抽出します。画像上でつながった画素の数を数える方法もありますが、ここでは面積を計算する方法を行ってみます。画素数よりも面積で表現した方が直感的にわかりやすくなりますが、面積を計算するために変化箇所画像からポリゴンを作成する必要があります。
画像のポリゴン化にはgdal.Polygonize()を使います。また、ここからOGRも登場します。OGRはGDALのベクター(ポリゴンデータなど)版と思ってください。
# ラスターのポリゴン化
from osgeo import ogr, osr
img_diff_mask = gdal.Open('diff_mask.tif', gdalconst.GA_ReadOnly)
# 座標系を取得
srs = osr.SpatialReference()
srs.ImportFromWkt(img_diff_mask.GetProjection())
# ポリゴン化するバンド及びマスクバンドを取得
imgband = img_diff_mask.GetRasterBand(1)
maskband = imgband.GetMaskBand()
# 出力先ポリゴンの用意
drv = ogr.GetDriverByName("ESRI Shapefile")
dst_ds = drv.CreateDataSource('diff_mask.shp')
dst_layer = dst_ds.CreateLayer('diff_mask', srs=srs)
# ポリゴン化
gdal.Polygonize(imgband, # ポリゴン化するバンド
maskband, # マスクバンド
dst_layer, # 出力先ポリゴンのレイヤー
-1, # 出力ポリゴンの属性フィールド(画素値が格納される)、-1とした場合は画素値は格納されない
["8CONNECTED=8"] # ある画素に隣接した画素を、4方向だけでなく、8方向で連結したポリゴンを作成
)
dst_ds = None
print("done")
ポリゴン化しましたが、座標系は緯度経度なので、このままでは面積の計算ができません。そこで、座標系をUTM(オリジナル画像の座標系EPSG:6690)に変換します。
# ベクターデータの座標変換
# ポリゴンデータの読み込み
s_ds = ogr.Open('diff_mask.shp')
s_lyr = s_ds.GetLayer()
s_srs = s_lyr.GetSpatialRef()
# 出力の座標系を指定
t_srs = osr.SpatialReference()
t_srs.ImportFromEPSG(6690)
# 座標系の変換方法を設定
coord_trans = osr.CoordinateTransformation(s_srs, t_srs)
# 出力先ポリゴンの用意
t_driver = ogr.GetDriverByName('ESRI Shapefile')
t_ds = t_driver.CreateDataSource('diff_mask_utm.shp')
t_lyr = t_ds.CreateLayer('diff_mask_utm',
srs=t_srs,
geom_type=ogr.wkbMultiPolygon # ベクターの形式をマルチポリゴンに設定
)
# 入力ポリゴンのレイヤー定義を取得
s_lyr_defn = s_lyr.GetLayerDefn()
# 出力ポリゴンに入力ポリゴンの属性フィールドをコピー
for i in range(0, s_lyr_defn.GetFieldCount()):
fieldDefn = s_lyr_defn.GetFieldDefn(i)
t_lyr.CreateField(fieldDefn)
# 出力ポリゴンのレイヤー定義を取得
t_lyr_defn = t_lyr.GetLayerDefn()
# 入力ポリゴンを出力の座標系に変換し、出力ポリゴンにコピー
for feature in s_lyr:
geom = feature.GetGeometryRef()
geom.Transform(coord_trans)
t_feature = ogr.Feature(t_lyr_defn)
t_feature.SetGeometry(geom)
# 出力ポリゴンに入力ポリゴンの属性値を格納
for i in range(0, t_lyr_defn.GetFieldCount()):
t_feature.SetField(t_lyr_defn.GetFieldDefn(i).GetNameRef(), feature.GetField(i))
t_lyr.CreateFeature(t_feature)
t_feature = None
feature = None
s_ds = None
t_ds = None
print("done")
座標系をUTMに変換したら、面積を計算します。ポリゴンの属性フィールドに'Area'を作成し、そこに面積を格納します。
# ポリゴンの面積を計算
# ポリゴンの読み込み、第2引数を'1'にすると上書きモード
poly = ogr.Open('diff_mask_utm.shp', 1)
poly_lyr = poly.GetLayer()
# 面積を格納する属性フィールドを作成
new_field = ogr.FieldDefn("Area", ogr.OFTInteger)
new_field.SetWidth(10)
poly_lyr.CreateField(new_field)
# 面積を取得し、属性フィールドに格納
for feature in poly_lyr:
geom = feature.GetGeometryRef()
area = round(geom.GetArea())
feature.SetField("Area", area)
poly_lyr.SetFeature(feature)
poly = None
print("done")
属性フィールドに'Area'が追加され、面積(m2)が格納されました。
ALOS-2の観測モードSM1の解像度は3mです。解像度の10倍程度あれば信頼性が高いと仮定すると、(3m×10)2 = 900m2となりますが、切りの良いところで1000m2としてポリゴンを抽出します。
# ポリゴンの面積で抽出
s_ds = ogr.Open('diff_mask_utm.shp', 0)
s_lyr = s_ds.GetLayer()
s_srs = s_lyr.GetSpatialRef()
# 面積が1000m2以上を選択
s_lyr.SetAttributeFilter("Area >= 1000")
# 選択したポリゴンを新しいポリゴンにコピー
t_driver = ogr.GetDriverByName('ESRI Shapefile')
t_ds = t_driver.CreateDataSource('diff_extract.shp')
t_lyr = t_ds.CreateLayer("diff_extract", srs=s_srs, geom_type=ogr.wkbMultiPolygon)
s_lyr_defn = s_lyr.GetLayerDefn()
for i in range(0, s_lyr_defn.GetFieldCount()):
fieldDefn = s_lyr_defn.GetFieldDefn(i)
t_lyr.CreateField(fieldDefn)
t_lyr_defn = t_lyr.GetLayerDefn()
for feature in s_lyr:
t_feature = ogr.Feature(t_lyr_defn)
geom = feature.GetGeometryRef()
t_feature.SetGeometry(geom)
for i in range(0, t_lyr_defn.GetFieldCount()):
t_feature.SetField(t_lyr_defn.GetFieldDefn(i).GetNameRef(), feature.GetField(i))
t_lyr.CreateFeature(t_feature)
s_ds = None
t_ds = None
print("done")
抽出結果を洪水浸水想定区域と重畳してみると、解析範囲内のほとんどの区域が浸水の影響を受けた可能性が見て取れますが、右下の区域では被害が少なかったと想定されます。
一方、方法①と同様に国土地理院の推定浸水範囲を重畳してみると今回の解析結果はやや過剰に抽出しており、改善の余地があるといえます。フィルターの種類・強さの変更、各閾値の調整、補助情報(マスク)の追加など色々と試してみてください。
真備町での浸水を対象として、被害範囲の推定を行ってみました。浸水の場合は土砂崩れ等と違い、災害発生後から時間が経つにつれて水が引いてしまうため、災害発生後すぐに観測された画像が必要になります。光学画像は天候による影響を受けるため、観測できるタイミングがあったとしても雲に覆われてしまって見えないということがあります。そこで役に立つのがSAR画像です。SARから照射されるレーダ波は雲を透過することができるため、平成30年7月豪雨時の悪天候であっても地表面の様子を捉えた画像を撮像することができました。また、使用した画像は夜間に観測されたものですが、SARは夜間でも観測できることが強みです。
今回用いた災害直後の画像はちょうど災害が発生した日に観測されたものでした。これは元々予定されていた観測ではありません。ALOS-2は災害時に計画を変更して被災地を観測することになっており、まさに平成30年7月豪雨の被害を把握するために観測され、今回はその画像を使用しました。このような緊急観測画像は災害発生直後に災害の概況把握をすることができる貴重な情報源として、実際に各関係機関に提供され、災害時の初動対応に役立てられています。各関係機関というのは日本だけではありません。宇宙技術によるアジア太平洋地域の災害監視を目的とした国際協力プロジェクト「センチネル・アジア」にもALOS-2の画像が提供され、アジア各国の災害時概況調査に貢献しています。
この章の前半で、平成30年北海道胆振東部地震を対象に行った差分干渉解析で使用した画像も地震発生日に観測された画像でしたが、これも緊急観測によるものです。今回は差分干渉解析によって地殻変動を確認してみましたが、同じタイミングで観測された画像ペア(災害前:ALOS2229472750-180823、災害後:ALOS2231542750-180906)で土砂災害の様子を確認することもできます。土砂崩れ等の場合は地表面の様子が著しく変化し、干渉性が低下しているため差分干渉解析を行うことはできませんが、そのような著しい変化が生じている場合は強度画像の差分をとることで被害状況を把握できます。興味があれば、実際に解析してみてください。以下のような図が作成でき、土砂崩壊が多数発生している様子(土砂崩壊箇所が赤色、崩壊土砂の堆積箇所が青色)が見て取れます。
この章ではTellusに搭載されたALOS-2/PALSAR-2プロダクトを使用しましたが、解析できる期間が限られています。その期間外で解析をしようとするとデータを購入する必要がありますが、SAR衛星には無償でデータ提供しているものもあります。ヨーロッパ宇宙機関(ESA)が運用中のSentinel-1は、Copernicus Open Access Hubから無償でデータをダウンロードすることができます。
災害直後の画像が無い場合もありますが、令和元年東日本台風(台風19号)の際はSentinel-1による観測が実施されています。使用するプロダクトタイプはGRD(ALOS-2のL1.5相当)です。GeoTIFF形式で取得できるので今回学んだ方法でRGB疑似カラー合成画像を作成してみると、埼玉県の川越市や坂戸市などの浸水状況を可視化することができます。
この章では、ALOS-2/PALSAR-2データを用いて、前半では位相情報を利用する差分干渉解析、後半では強度情報を利用して2時期の差分解析を学習しました。SAR画像は私たちが目で見ているようなイメージとは違うものなので理解するのが大変ですが、観測時の天候に左右されにくく、夜間でも観測可能なので、特に災害時に真価を発揮します。この教材で学習した方法をほかの災害についても試してみてください。SAR画像についての理解が深まると思いますし、今までとは違った(衛星を使った)視点から災害を見ることで、新しい気付きがあるかもしれません。ただし、災害の規模や様態によっては被害状況が必ずしも得られるわけではなく、衛星による調査はあくまでも概況調査であるということを忘れないでください。