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

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

一方で、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へお問い合わせください。

まずは、衛星データのダウンロードです。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データを使って、干渉解析という高度な処理が必要になってきますので、ここでは割愛します。

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

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

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

(1)しきい値の決定

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

(3)森林域の除外

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

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

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

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

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

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

(4)ポリゴン化

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

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

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

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

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

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

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

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

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

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

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

次に、再び「バッファ」ツールを開き、先ほどの実行結果を入力レイヤとし、「距離」には先ほど記入した数値に「マイナス」を付けて実行します（-10mや、-0.0001 degなど）。

このようなポリゴンが作成されました。拡大すると、ポリゴン内部の小さい穴が埋まっていたり、小さいポリゴンが消えたりしていることが分かると思います。

(6)微小ポリゴンの削除

最後に、小さいポリゴンを削除して終わりです。

レイヤリストに新しく追加されたポリゴンを右クリックで選択し、「編集モード切り替え」を押して編集モードにします。

手動で選択して削除する場合は、「シングルクリックによる地物選択」の「地物を選択」を選択した上で、消したいポリゴンを地図上で選択し、「選択物の削除」を押下します。

また、ある面積以下を一括で削除することもできます。

筆者の環境では各ポリゴンに面積の情報が付与されていなかったので、まず面積情報を付与します（面積情報が既にある場合はスキップ）。

レイヤリストからポリゴンを右クリックし「属性テーブルを開く」を押下します。属性テーブルが開いたら、上にある「フィールド計算機を開く」を押します。

出力する属性（フィールド）の名前を「area」など好きな文字とし、式に「$area」を入力します（画面中央の「ジオメトリ」＞「$area」を選択すると入力されます）。

「OK」を押すと、「area」という属性情報が追加されます。

次に、「値による地物選択」を選択します。

先ほど計算した「area」に対して、削除したい面積の条件を入力し、「地物を選択」を押下します。（筆者は「2000より小さい」という条件を入れました。）

条件に合ったポリゴンが選択されますので、「選択物の削除」を押します。選択されたポリゴンが削除されます。

編集したポリゴンは、忘れずに保存しましょう。「レイヤ編集内容を保存」を押し、編集モードを解除します。これで完成です。

「3. RGBカラー合成で浸水域を目立たせる方法」で作成したRGBカラー合成画像と重ねるとこのようになります。今回抽出した浸水域を黄色い枠で重ねています。RGBカラー合成図の赤い箇所が、抽出できています。

5. まとめ

本記事では、QGISを使って、ALOS-2のSAR画像から浸水域を抽出する方法についてまとめました。ご紹介したのは初歩的なものでしたが、SARの解析方法についてイメージを掴んでいただけたでしょうか？

Tellusで簡単にALOS-2等のSARデータがダウンロードできるので、是非試してみてください。

また、本記事でご紹介したTellusのQGIS付きの開発環境「Tellus Satellite Data Master with QGIS」は、現在Tellusのサイトからスペックと費用がご確認いただけます。お申込みはTellusへお問い合わせください。