Tellus 2022/9/22

QGISを使って、SARから浸水域を抽出してみよう！意外と簡単？SARの解析方法を紹介

Tellusで提供しているQGIS付き開発環境を使って、PALSAR-2の浸水域抽出の方法をご紹介します。

日本は梅雨期の大雨や台風などによって水害が毎年のように発生します。最近では、パキスタンの豪雨の影響で3300万人以上が被災し、1300人以上が亡くなるという甚大な被害も発生しました。水害を未然に防ぐことも大変重要ですが、ひとたび被害が起こってしまった際に、少しでも早く被害の全容を把握し、適切な場所に適切な支援を行えるようにすることもまた重要です。

浸水域の全容把握は、一度に広い範囲の観測ができる人工衛星が得意とするところです。中でも水害には、「合成開口レーダ（SAR: Synthetic Aperture Radar）」の画像を使うのがおススメ。SARは、雲を透過するマイクロ波を使って観測を行うので、上空が厚い雨雲に覆われていても、それを透過して地上の様子を観測することができるからです。

SAR画像から浸水域を抽出した様子 Credit : Original data provided by JAXA

一方で、SARは、画像を見ても白黒でよく分からないし、観測原理を調べても難しい数式がたくさん出てくるしで、なんとなく上級者向けのイメージがありますよね。しかし、浸水域の抽出は、意外と簡単にできる方法があるんです。

この記事では、SARで浸水域を抽出するメカニズムと、QGISを使った簡単な解析方法についてご紹介します。簡単に解析できてしまうので、是非トライしてみてください！

1. SARで浸水域を抽出するメカニズム

SARは、マイクロ波を発射し、地表で跳ね返ってきたマイクロ波をとらえるセンサです。ざらざらした表面ほど多く電波が返ってきて白く見え、水面などつるつるした表面では電波が反射してしまうため黒く見えます。浸水して地表面が滑らかになった場所では、電波が衛星側に返ってこないため、受信信号が小さく（画像が暗く）なります。

SARのこの性質を利用して浸水した地域を推定することができます。本記事では、初歩的な２通りのやり方をご紹介します。

１つ目は、しきい値を決めて浸水域を抽出する方法です。発災後のSAR画像から、浸水と思われる箇所を見つけたら、その場所と同じか、より暗い領域を、浸水域として大まかに抽出します。

２つ目は、RGBカラー合成で浸水域を目立たせる方法です。RGBカラー合成を行い、災害が発生する前と後の画像を比較して、浸水した地域を強調させます。

今回は、災害前の画像を赤色に、災害後の画像を青色と緑色にして画像を合成します。災害後の青色と緑色の画像は、光の三原色から、重ねると水色にみえます。この処理によって、災害前の受信信号（衛星に返ってくる電波）が災害後と比べて大きい場合は赤色に、災害後の値の方が大きい場合は水色に、また受信信号に変化がない場合は白色や黒色に近い色で示される図を作成することができます。

浸水した場所は電波が反射して衛星側には返ってこなくなるため、災害前画像の方が受信信号が大きい箇所となり、このRGBカラー合成図では赤く強調して表示される、という仕組みです。

ビジネス事例 髙倉有希

災害時の衛星画像はどこを見に行けばいい？災害時の衛星画像や解析結果を共有する取り組みを紹介！

2. SAR画像のダウンロード

それでは早速、SAR画像から浸水域を抽出していきましょう。今回は、Tellusで販売している、TellusのデータにアクセスできるアドインがあらかじめインストールされたQGIS付きの開発環境「Tellus Satellite Data Master with QGIS」を使って、JAXAの「だいち2号」（ALOS-2）を対象に解析します。対象は「令和２年７月豪雨」の熊本県人吉市球磨川周辺とします。

※TellusのQGIS付きの開発環境「Tellus Satellite Data Master with QGIS」は、現在Tellusのサイトからスペックと費用がご確認いただけます。

まずは、衛星データのダウンロードです。TellusのQGIS付きの開発環境を使えば、QGIS画面上から欲しいデータが簡単に入手できるので、今回はこちらを利用します。なお、初期設定は省略していますので、Tellusから別途配布されるマニュアルを参照ください。今回使用したQGISのバージョンは3.22です。

では、QGISを立ち上げ、左側のメニューから「XYZ Tiles」を選択し、「OpenStreetMap」をダブルクリックします。

OpenStreetMapが表示されるので、衛星データをダウンロードしたいエリアに移動します。今回は、熊本県人吉市付近に移動します。

QGIS画面右上にあるTellusのプラグインアイコンをクリックします。検索条件は以下としました。

・対象データセット：【Tellus公式】PALSAR-2_L2.1
・日付で絞り込む：2020/07/01-2020/07/31
・AOI：「キャンバスに描画」をクリックし、人吉市を中心に描画

PALSAR-2は、ALOS-2のセンサ名です。またL2.1は、PALSAR-2データの処理レベルを示しています。L2.1は、数値標高データを用いて幾何補正（オルソ補正）が行われているデータで、扱いやすいため、初心者はまずこちらを利用するのが良いと思います。

検索結果は以下になります。4件の結果が表示されていますね。今回は「2020/7/4」の昼観測のデータをダウンロードしますので、該当のデータを選び、下の「ダウンロードするファイルを選択してください」をクリックして「ALOS2330172910-200704_webcog.tif」を選択します。保存先フォルダを指定して、ダウンロードを開始します。

マップ追加に関するダイアログが表示されるので「はい」を押します。すると、QGISの地図上に衛星データが表示されます。

とっても簡単に衛星画像をダウンロードすることができますね！TellusのQGIS付き開発環境が利用できれば、衛星データのダウンロードが非常に簡単にできて、とても便利だと思います。

3. RGBカラー合成で浸水域を目立たせる方法

まず先に、より簡単な方法である「RGBカラー合成で浸水域を目立たせる方法」を行います。

ここでは、発災前の画像も使いますので、先ほどと同じ手順で発災前の画像をダウンロードします。可能であれば、発災後画像と同じ観測条件（軌道方向、観測方向、観測モード、オフナディア角など）のものが望ましいです。今回は、「2016/4/16」の昼のデータをダウンロードします。（今回の事例は、浸水抽出にはあまり使用しないオフナディア角が大きい観測条件だったので、過去画像がかなり古くなっていました。）

次に、QGIS上部の「ラスタ」タブから「その他」＞「仮想ラスタの構築」を選択します。

著者の環境では、なぜかこの画面だけ英語なのですが、まず「Place each input file into a separate band」にチェックを入れ、その後、「Input layers」の右のボタンを選択します。

そして、発災前（2016/4/16）、発災後（2020/7/14）の順にチェックを入れ、下部の「実行」を押下します。

新しいレイヤが作成されました。

QGIS画面左側のレイヤリストにある、新しく作成されたレイヤ名を右クリックし、「プロパティ」を選択します。そして、バンドレンダリングの「青のバンド」で「バンド2」を選択し、「OK」を押します。

すると、発災前画像を赤、発災後画像を緑と青に割り当てて合成した、RGBカラー合成画像が作成できました。

このままでは見にくいので、先ほどのレイヤプロパティのレイヤレンダリングのうち、「輝度」や「コントラスト」等を上げることで、見やすくすることができます。

赤い部分（浸水域）が見やすくなりましたね。

国土地理院が、SNS画像やヘリが撮影した画像等を元に、浸水推定図を作成し公開していますので、今作成した画像に重ねてみましょう。

「令和2年7月豪雨に関する情報」は以下サイトから確認できます。
https://www.gsi.go.jp/BOUSAI/R2_kyusyu_heavyrain_jul.html

球磨川流域球磨川（2020年7月4日20時作成）の浸水推定図の浸水範囲の輪郭線（GeoJSON）をダウンロードし、先ほど作成したRGBカラー合成図の上に黄色枠として重ねると、このようになります。

どうでしょうか。RGBカラー合成画像で赤くなっている部分は、国土地理院が浸水域として示している箇所とだいたい合っていますね。

しかし一方で、RGBカラー合成図では、東側の市街地の浸水域は拾えていないことがわかります。市街地が完全に冠水していれば赤くなるのですが、住宅やビル等が完全に冠水することはほとんどありません。そのため、RGBカラー合成画像は、建物が密集しているような場所では、浸水していても赤くならないことが多いのです。

このように、SARでは、建物域の浸水域を抽出するのは原理的に簡単ではないのですが、解析方法に関する研究が進められており、いくつか方法があります。

例えば、浸水時は、水面と建物の壁で2回反射が起こり、むしろSARの受信信号が強くなることがあるため、それを抽出する方法（ここでは、水色になっている部分がそうです）。あるいは、干渉SARという解析手法で得られる「コヒーレンス（干渉性）」の減少から捉える方法。

後者の方が、より確実に抽出できるのですが、複素数データであるL1.1データを使って、干渉解析という高度な処理が必要になってきますので、ここでは割愛します。

衛星データ入門

干渉SAR（InSAR）とは-分かること、事例、仕組み、読み解き方-

4. しきい値を決めて浸水域を抽出する方法

次に、1つ目の方法である、「しきい値を決めて浸水域を抽出する方法」を行っていきましょう。全体的なフローは次の通りです。

(1)しきい値の決定
(2)しきい値以下の抽出
(3)森林域の除外
(4)ポリゴン化
(5)ポリゴンの平滑化（単純化）
(6)微小ポリゴンの削除

(1)しきい値の決定

浸水している可能性がある箇所にズームします。SAR画像は、表面の滑らかな箇所ほどSARの受信信号が小さく黒く表示されます。そのため、背景地図と見比べて、海・川・湖以外で黒い箇所があれば浸水している可能性があります。なお、SAR画像で白くなっている箇所は、反射が強いところであり、人工的な構造物が多い市街地などです。

左図：SAR画像、右図：OpenStreetMap
Credit: JAXA（左図）、OpenStreetMap contributors（右図）

①しきい値を決める簡単な方法

QGISの上方右側にある「地物情報表示」のアイコンをクリック
Credit : JAXA

浸水の可能性がある箇所の、SAR画像の画素値をチェックしましょう。QGISの「地物情報表示」のアイコンを押下してから、画像をクリックすることで、画素値が表示されます。暗いところを何か所かクリックすると、浸水域とそうでない場所のしきい値が見えてくると思います。（筆者はこの画像をみて、50くらいかなと思いました。）

一旦、えいやとしきい値を決めてみて、次に進みましょう。うまくいかなければ少ししきい値を変更して繰り返し、いい感じのしきい値を見つけてみてください。

なんて主観的なやり方…と感じるかもしれませんが、実際に浸水していることが報道されている箇所が1か所でも分かれば、そこと推定結果を見比べながらしきい値を調整することで、より確からしい全体的な浸水域を推定することもできます。

②より客観的にしきい値を決める方法

客観的にしきい値を決めたいという場合は、水域（海、川、湖など）のSARの画素値の平均などからしきい値を決めることができます。これは、自動的に解析したい場合や、現地の浸水域情報がなくしきい値を調整できない場合等に向いている手法です。

水域のSAR画像の値を知るために、土地被覆図を用意する必要があります。日本であればJAXAの「高解像度土地利用土地被覆図」が無料でダウンロードできますので、解析範囲と対応したデータをQGISに追加して表示します。
https://www.eorc.jaxa.jp/ALOS/jp/dataset/lulc_j.htm

左図：SAR画像、右図：土地被覆図
Credit: JAXA

JAXAの「高解像度土地利用土地被覆図」は、デジタル値が1の画素が水域 (Water bodies)を示します（上図では、濃い青色の箇所）。

QGIS上部の「ラスタ」タブから「ラスタ計算機」を呼び出し、式に以下を入力します。

if (”土地被覆データ” = 1, 1, 0/0)

出力レイヤ名を入力し「OK」を押下すると、水域が１、それ以外がnodataというデータになります。

次に、作成したラスタをポリゴンに変換します。「ラスタ」＞「変換」＞「ラスタのベクタ化」を選択し、実行します。

次に、「ラスタ」＞「抽出」＞「マスクレイヤで切り抜く」で、SAR画像を、先ほど作成した水域のポリゴンで切り抜きます。

これで、SAR画像のうち、水域部分が切り出されました。

レイヤリストに新しく追加されたラスタを右クリックし、「プロパティ」の「情報」を見ると、水域で切り抜かれたSAR画像の統計情報が表示されています。

「STATISTICS_MEAN=34.752987514579」とありますので、この画像の水域のSARの平均値は「約35」であることがわかりました。今回はしきい値を平均値である35として進めていきます（他にも、平均と標準偏差を組み合わせてしきい値を決定することもあります）。

(2)しきい値以下の抽出

では、SAR画像のしきい値以下の部分を抽出していきます。

ここからの処理は時間がかかるので、注目している領域があれば、先にそこだけを切り出してから先に進まれることをおすすめします（必須ではありません）。

切り出す場合は、「ラスタ」＞「抽出」＞「範囲を指定して切り抜く」で、入力レイヤにSARデータを指定して、切り抜く範囲を入力します。「実行」を押せば切り出しが実行されます。

クリップが終わったら、「ラスタ」＞「ラスタ計算機」に以下の式を入力し、しきい値（今回は35）以下を抽出します。

if (”クリップしたSAR画像のデータ” <= 35, 1, 0/0)

しきい値以下の領域が画素値１として出力されました。

(3)森林域の除外

浸水は一般に平坦な土地で生じるという性質を利用して、浸水の可能性が低い箇所を候補から除外していきます。ここでは、森林域を除外する方法を紹介します。他にも、数値標高モデル（DEM）を用いて、土地の傾斜角が一定以上の箇所を除外する方法もあります。

「(1)しきい値の決定」の②で使用した「高解像度土地利用土地被覆図」を使用します。（ここでも、注目している領域があれば、先にそこを切り出しておくのがおすすめです。）

「高解像度土地利用土地被覆図」では、デジタル値6〜9と11の画素が森林域を示すので、「ラスタ」＞「ラスタ計算機」で以下の式を入力し、「OK」を押下します（今回は、森林域以外を１として出力することに注意です）。

if (“土地被覆データ” = 6 OR “土地被覆データ” = 7 OR “土地被覆データ” = 8 OR “土地被覆データ” = 9 OR “土地被覆データ” = 11, 0/0, 1)

森林域以外の画素値が１で示されるデータができたら、次にそれをベクタ化するために、「ラスタ」＞「変換」＞「ラスタのベクタ化」を選択し、実行します。

森林域以外を示すポリゴンが作成されました。

そして、「ラスタ」＞「抽出」＞「マスクレイヤで切り抜く」で、(2)で作成したラスタを入力レイヤに、先ほど作った森林域以外のポリゴンをマスクレイヤとして実行すると、このように、森林域が除外されたデータが作成できます。

(4)ポリゴン化

先ほど作成した、森林域が除外されたラスタを、ポリゴンに変換します。

「ラスタ」＞「変換」＞「ラスタのベクタ化」を選択し、先ほど作成した画像を入力レイヤに指定して「実行」します。

ポリゴンが作成されました。

(5)ポリゴンの平滑化（単純化）

このままでは、ポリゴンの形状が複雑なので、形状を滑らかにしていきます。方法は色々ありますが、一例として「クロージング処理」での方法をご紹介します。クロージング処理は、今回のように2値化された画像に対して、膨張処理の後、収縮処理を行うことによって、ノイズを削除する方法です。

「ベクタ」＞「空間演算ツール」＞「バッファ」を選択します。

作成したポリゴンを入力レイヤに適用し、「距離」は10 m（「地理座標系」の場合は0.0001 deg）などとし、その他はデフォルトのまま「実行」します。

【参考】
バッファツールの「距離」の単位は、ポリゴンの座標系が「投影座標系」の場合メートルに、「地理座標系」の場合deg（度）になります。

今回は必須ではありませんが、メートルの方が馴染み深いので、使用するデータやプロジェクトの座標参照系（CRS）は、初めから「投影座標系」に変換しておいてもいいかもしれません。

投影座標系への変換候補としては、「WGS_1984_Web_Mercator_Auxiliary_Sphere」や「JGD2011 / UTM zone 〇〇N」等があります。

また、ここで使用する「バッファ」ツールでは、入力する「距離」が大きいほど、結果は単純になります。「距離」に入力する適切な値は、浸水域の大きさによって異なる可能性があるため、いくつか試して調整してみてください。