解析ノートブック 2020/2/26

【コード付き】復興状況を人工衛星からモニタリング！北海道地震後の夜間光から確認してみた

本記事では、広域にわたる月次データが安定して提供されてきたVIIRS夜間光画像データを使って被災地の長期モニタリングの例を示します。夜間光データからは、人の社会経済活動が分かるとされています。

近年、日本やアジアでは規模の大きい災害が多発しています。最近では、日本国内ですと2018年西日本大豪雨、2019年首都圏台風・豪雨で甚大な被害を受けました。

日本国内だけでなく、アジアでも多数の災害に見舞われ、2013年フィリピン・ハイヤン台風、2018年インド・ケララ州台風等、大きな被害をもたらす災害が発生しました。このような自然災害について、気候変動により想定外の災害が多くなると言われています。
災害の発生をコントールすることが極めて難しい一方で、災害が起こったときの対応の改善を目的とした国際的な枠組み、「仙台防災枠組2015-2030」が2015年に合意されました。

この枠組みでは宇宙技術や地理空間情報技術の利活用が課題に挙げられており、Tellusが提供する衛星データ分析プラットフォームが貢献する領域です。

災害の被害を議論するときの基本的な考え方として、災害が発生した位置や範囲だけでなく、その領域にどの程度の社会経済が影響を受けるかを計量する必要があります。具体的には、どれだけ大規模な台風が海上に発生しても陸に上がるまでに弱くなっていれば被害は比較的小さいと言えます。

逆に、地震の直接的な影響がなくても、電力の供給元である発電所が地震によって稼働停止すれば、大変な影響を被ることになります。
さらに、被害の位置・領域を把握するだけでなく、復旧・復興の進捗をモニタリングする必要もあります。復旧・復興は、被害の把握に比べて長期的になされるため、長期的な観測が求められます。地上踏査だけでは、長期にわたる広域観測の費用が甚大な規模に及びますが、地球観測衛星は情報の精細さは落ちるものの、地上踏査に比べてきわめて低い費用で実現できます。

本記事では、広域にわたる月次データが安定して提供されてきたVIIRS夜間光画像データを使って被災地の長期モニタリングの例を示します。

1. 夜間光画像から得られる社会経済の情報

地球観測衛星による夜間光画像は、もともとDMSPという衛星から軍事活動の検出を目的に観測されてきましたが、建物等から漏れる人工光が社会経済の状況を表すという着想から、観測から得られる夜間高強度を使った社会経済指標の推定に役立てるための研究が重ねられてきました。

具体的には、電力消費量、GDP、人口密度が挙げられます。2012年にはDMSPに加えて、VIIRSという衛星による夜間光画像データが利用できるようになりました。DMSPは空間分解能が1 kmであるのに対し、VIIRSは0.5 kmに向上したうえ、データレンジが拡がり情報量も向上しました。

夜間光画像は雲の影響を強くうけるため、ある箇所を毎日通過するとしても、必ずしも地上の様子が毎日観測されるとは限りません。このような制約に対し、NOAA NGDCはDMSPやVIIRSによる夜間光画像を使いやすい形式で配信するための研究開発を積み重ねてきました。

具体的には、画素毎に1年間や1ヶ月間のうち雲の影響が無い日の観測の平均値などを計算して、1年間や1ヶ月間の観測値として集約することで、年別・月別の雲無し夜間光画像を作成しています。なお、2019年10月にNOAA NGDCはデータ更新を終了し、現在ではColorado School of Minesがデータを処理・配信しています。

夜間光画像からGDPや人口等の社会経済情報を得ることができれば、時間と費用がかかる統計調査の代替として、年別・月別で社会経済の変化を分析することができます。こうして得られるデータは、災害の復旧・復興モニタリングにとても有用です。

というのも、災害からの復旧・復興は長い期間を要し、調査の費用は被災後の地域にとって重荷になることが想像されます。被災地が広域にわたる場合は費用が更に倍増するだけでなく、都市部以外での復旧・復興の遅れを見逃すおそれもあります。夜間光画像では地上の詳しい様子を知るには限界がありますが、広域動態の全体像を長期にわたってモニタリングするにはうってつけです。

以下、まだ検証が必要な手順ではありますが、2018年北海道胆振東部地震を対象とした夜間光画像を使った復旧・復興モニタリングの例を示します。十分な検証がなされていない方法につき、結果はあくまで試行結果としてご覧ください。

2. 夜間光画像を使った復旧・復興モニタリングのデータ処理

災害直後は被害のために社会経済活動が停滞することが多くあります。その状況は、夜間光画像では夜間光強度の減少として表れます。逆に、災害後に復旧の進行につれて社会経済が活発になると、災害直後に比べて減少の度合いが小さくなる、あるいは増加に転ずることが考えられます。この仮説に基づき、月別の夜間光画像を使って、月毎の変化の分析を試みます。

2.1. 範囲と時期の設定

北海道胆振東部地震は2018年9月6日に起こりました。震災前後の比較のため、2018年8月から変化を分析します。実際の状況については例えばこちらのリンクをご覧ください。

今回のスクリプトでは画像を処理するのにシェルスクリプト（bash）を通してGDALのユーティリティプログラムを利用します。自身の環境にGDALのユーティリティプログラムをインストールしていない、もしくはスクリプトが動かない場合には、UbuntuやDebian等のAPT系の環境で、下記コマンドを実行してください。

sudo apt update; sudo apt install gdal-bin

Anacondaでデータ処理するには下記コマンドでインストールしてください。

conda install -c conda-forge gdal

以下のセルでは、緯度経度の範囲と、データ分析の時期始めを指定します。緯度経度の正確な座標値はGoogle Earth等を使って確認します。

XMAX=142.0
XMIN=140.9
YMAX=43.39
YMIN=42.72

Y_BEGIN=2018
M_BEGIN=8

2.2. データのダウンロード

VIIRS夜間光画像の月別雲無しコンポジット画像のファイル名・パスは、年月を含む体系的な命名法なので、繰り返し処理により自動化します。以下の例では上記に示した時期始めから5ヶ月間のデータをダウンロードします。

mkdir -p src
cd src
for i in $(seq 0 4); do
    TARGET=$(date +%Y%m --date="$Y_BEGIN/$M_BEGIN/01 $i month")
    BASE=https://eogdata.mines.edu/wwwdata/viirs_products/dnb_composites/v10/
    wget -nd -nc -r $BASE/$TARGET/vcmcfg/ --accept-regex=".*SVDNB_npp_${TARGET}01-${TARGET}.*_75N060E_.*.tgz"

    rm -f index.html

 # ファイルサイズに注意。展開後のTIFはファイルサイズが2GBを超えるものもある。
    TGZ=*SVDNB_npp_${TARGET}01-${TARGET}*.tgz
    tar xaf $TGZ $(tar taf $TGZ | grep "SVDNB_npp_.*_vcmcfg_v10_.*.avg_rade9h.tif")
done
cd ..

3. データの処理

衛星リモートセンシングの変化検出の手法のうち、最もシンプルな方法として、差分画像に対しμ±2σ（μ:平均値; σ:標準偏差）を上限・下限として、それらを超える画素値を変化とみなす方法があります。正規分布においては、μ±2σの区間に全データの95%が含まれ、上限・下限を超えるデータは2.5%ずつになります。2.5%は全体の傾向に比べて異常に大きい・小さいことから、何かしらの大きな変化が地上にあったとみなします。これにより対象領域の全体的なバイアスを除くことができます。

今回の試行では、震災が起こってから5ヶ月間にわたって月毎の変化の過程を分析します。8月から9月には震災により停電を始めとした社会経済への影響があったことから、夜間光強度は減少すると考えられますが、9月以降は復旧・復興の進行に従って、増加に転ずると考えられます。

以下のセルでは、対象範囲で切り抜いてから、1ヶ月毎の差分を計算し、それぞれの差分についてμ±2σを閾値として変化を検出します。

# 指定範囲で切り抜き
mkdir -p subset
for TIF in src/*.tif; do
    gdal_translate -projwin $XMIN $YMAX $XMAX $YMIN $TIF subset/subset-$(basename $TIF)
done

mkdir -p diff change
for i in $(seq 0 3); do
    # 1ヶ月ずらして YYYYMM をセット
    TARGET_0=$(date +%Y%m --date="$Y_BEGIN/$M_BEGIN/01 $i month")
    TARGET_1=$(date +%Y%m --date="$Y_BEGIN/$M_BEGIN/01 $(expr $i + 1) month")
    
    # 連続する2ヶ月間の差分を計算
    DIFF=diff/diff-${TARGET_1}_${TARGET_0}.tif
    gdal_calc.py --calc="A-B" --outfile=$DIFF --type=Float32 -A subset/*_${TARGET_1}01-*.tif -B subset/*_${TARGET_0}01-*.tif 

    # 差分の統計値を計算し、μ±2σを計算
    eval $(gdalinfo $DIFF -stats | grep STATISTICS)
    UPPER_TH=$(perl -e "print $STATISTICS_MEAN + 2 * $STATISTICS_STDDEV")
    LOWER_TH=$(perl -e "print $STATISTICS_MEAN - 2 * $STATISTICS_STDDEV")
    
    # diff > μ+2σ は 3, diff < μ-2σ は 1,  μ-2σ < diff  $UPPER_TH, 3, where(A > $LOWER_TH, 2, 1))" -A $DIFF --outfile=$CHANGE_RAS --type=Byte
    gdal_sieve.py -st 3 $CHANGE_RAS # 極小ポリゴン（3ピクセル未満）の除外
    gdal_polygonize.py $CHANGE_RAS -f GeoJSON $CHANGE_VEC
done

はじめに夜間光画像を切り抜いた結果をQGISで表示します（図1）。QGISをダウンロードしていない方は、無料なのでこちらからダウンロードのうえ、データを確認してみてください。

札幌市市街地を中心に高い輝度値が拡がる様子がわかります。8月の観測では対象地域北部が全てゼロになっていますが、データの抜け落ちと思われます。12月には輝度値の高い範囲が他の月に比べて拡がっています。これは、クリスマスをはじめとした休暇シーズンによる観光旅行者の増加にともなう経済活動の増加が要因と考えられます。

図1：夜間光画像を月別に表示した結果。左上：2018年8月、右上：同年9月、左中：同年10月、右中：同年11月、左下：同年12月

それぞれの月の夜間光画像を見比べるだけでは時系列の変化をとらえることは難しいので、変化検出の結果をQGISで表示します（図2）。8月から9月にかけて減少変化の領域が札幌市を中心に周辺の都市部で見られます。北東部の岩見沢市周辺の都市部で夜間光強度の増加変化が見られますが、2018年8月のデータ抜け落ちにより輝度値がゼロになってしまったことが原因です。

9月から10月にかけて夜間光強度の減少変化が続いていますが、領域は小さくなっていることから復旧の成果が表れていることが想像されます。10月から11月にかけては大部分で夜間光強度の増加変化が見られ、11月から12月にかけては全ての領域で増加変化が見られました。この変化の動態は復旧・復興の過程に即しているように思われます。

以上の結果から、夜間光強度の変化の観測が復旧・復興の進行具合をモニタリングするのに役立ちそうという期待はありますが、夜間光画像は季節性がきわめて強い時系列データのため、季節性を調整する必要があります。特に、12月はクリスマス等のために都市部で際立って夜間光強度が高くなる傾向があり、復旧・復興の進行状況にかかわらず前月に比べて増加する可能性はあります。

まだまだ検証と改良が必要な手法ではありますが、これが上手くできれば、大規模災害の後に長期的な復旧・復興の進捗を効果的に管理し、支援の見落としや投入の過度な集中を避けることができます。この手法の研究開発にご興味のある方はお気軽にご連絡・ご相談ください。