衛星画像による山火事の変化検出
— 2023年 ギリシャ・Evros 火災
Sentinel-2 の2時期差分・自動閾値・後処理・定量評価を、コア部分は自前実装で通した学習記録。 地理空間・リモートセンシングは未経験からの着手。
2023年8月にギリシャ北東部で発生した山火事を対象に、 発火前後の Sentinel-2 画像から焼失域を検出し、Copernicus EMS の公式被害評価と突き合わせた。 分光指数・自動閾値・モルフォロジ・評価指標はライブラリを呼ばず自前で実装し、 各段で確立ライブラリと数値を突き合わせている。
数値の読み方。 「IoU 0.87」は参照の定義を書かなければ意味を持たない。 同じ予測・同じ実装でも、参照を「延焼域全体 / 中被害以上 / 高被害以上」に変えるだけで IoU は 0.56 / 0.87 / 0.63 と動く。参照自身の Overall Accuracy は 96.33% だが、 OA と IoU は異なる指標なので、この値を IoU 差の測定限界には使わない。以降の数値は すべて参照定義とセットで示す。
第 1 層成果
第 2 層パイプラインと技術的な選択
各段で「何をしたか」より「なぜその選択をしたか」を残した。 地理空間データ特有の落とし穴は、ほとんどがこの選択の中にある。
B0 — 対象の選定(AOI と STAC 検索)
- Planetary Computer の STAC API で
sentinel-2-l2aを bbox・期間・雲量で検索。 - pre/post を同一 MGRS タイル(35TMF)から選んだ。同一タイルなら同一 UTM・同一原点なので、 後段の位置合わせで再投影の補間を挟まずに済む。これが B1 を無調整化した。
- L1C ではなく大気補正済みの L2A。時期間比較で大気差をノイズにしないため。
- 結果: pre 2023-08-05(雲 0.0%)、post 2023-09-12(雲 0.5%)。火災期間(8/19–8/29頃)を挟む。
B1 — 取得・整列(共通格子への co-registration)
- 各バンドは COG なので、110km のタイル全体ではなく AOI の窓だけを部分取得できる。
- AOI の経緯度 bbox を UTM へ変換し、10m 格子に snap して基準格子を作る
(EPSG:32635、1957×2462 px、約 19.6×24.6 km)。全バンドをここへ
reprojectする。 - リサンプリングをデータの性質で使い分ける。 反射率(連続値)は双一次、 SCL(シーン分類=カテゴリ)は最近傍。「4=植生」と「9=雲」を平均した 6.5 は存在しないため。
- L2A harmonize オフセット。 処理ベースライン 04.00 以降は
reflectance = (DN − 1000) / 10000。引き忘れると反射率が一律に底上げされる。 - SCL で雲・影・欠測を除外。今回は両時期とも有効画素 100%。
B2 — 分光指数
NDVI / NDWI / NBR はすべて正規化差 (a−b)/(a+b) という同じ型で、
バンドの割り当てが違うだけ。和で割ることで地形陰影・太陽高度の乗法的な違いが約分され、
分光的な形状だけが残る。一方 L2A のオフセットのような加法的なずれは分子で消えても
分母に残るため消えない(だから B1 で先に処理する)。
焼失に効くのは NBR。焼失は NIR 崩落と SWIR 上昇の両効きで分子を二重に押し下げる。 実測でも平均の低下幅は NDVI 0.26 に対し NBR 0.37 だった。
B3 — 変化検出(dNBR と Otsu 自動閾値)
dNBR = NBR_pre − NBR_post。指数を先に作ってから差分する (先にバンドを引くと、差の中に照度差が残ったままになる)。- 分類後比較ではなく画像差分。分類を経ると各時期の誤差が掛け算で積み上がる(0.9×0.9=0.81)。
- Otsu はクラス間分散を最大化する閾値を選ぶ。累積和2本で全閾値候補を
O(B)で一括計算し、 閾値だけでなく分散曲線そのものを配列で返す設計にした(安定度を見るため)。 - 結果: 閾値 +0.4469、変化 37.4% / 180.2 km²。
skimage.filters.threshold_otsuと 3.3e−08 で一致。
B4 — 後処理(空間的な事前知識の注入)
B3 のマスクは画素ごとに独立に作られている。焼跡は空間的にまとまるという知識を、
モルフォロジと最小マッピング単位で後から入れる。収縮・膨張・連結成分ラベリング
(ラン長 + Union-Find)は自前実装で、scipy.ndimage と厳密一致を確認した。
この段では設定を一つに決めなかった。後処理の強さは偽陽性の除去と真陽性の欠損の トレードオフだが、この時点では参照データが無く優劣を判定できない。勘で選ぶとその選択が 次段の評価に混入して検証できなくなるため、4段階すべてを次段に渡した。
B5 — 定量評価
参照は Copernicus EMS EMSN166(この火災の事後被害評価)。着手前に技術報告書を精読し、 以下を確認してから設計を決めた。
| 確認項目 | 内容 | 評価への影響 |
|---|---|---|
| post 画像 | Sentinel-2 2023-09-12 09:06 | 私たちと同一撮影。時刻ずれゼロだが循環性は最大 |
| 手法 | dNBR + 目視検証 | 同系統 = 整合性チェックであり独立検証ではない |
| MMU | 2500 m² | こちらの 1 ha より 4 倍細かい |
| 焼け残り | 2500 m² 未満は焼失域に含む | 内部で系統的な見逃し |
| 参照の精度 | Overall Accuracy 96.33% | 参照にも誤分類があることを示す補足情報 |
最良は s2_close × 中被害以上で IoU 0.8667
(Precision 0.9665 / Recall 0.8936)。
ただし s3_mmu との差は 0.0007 で、参照自身の精度よりはるかに小さい。
したがって「s2 が最良」とは書かず、「s2 と s3 は同等、どちらも後処理なしより良い」とする。
第 3 層学習の過程 — 外した予測と訂正
このケーススタディで最も学びが大きかったのは、うまくいった箇所ではなく 予測を立てて外した3回だった。いずれも外れ方が新しい理解を生んでいる。
① 二山性は差分に引き継がれる — 半分外れ
B2 で post の NBR が二山になったので、dNBR にも引き継がれると予測した。実際には 鋭い無変化ピーク + 平坦な肩という非対称な形になった。焼失強度が連続的に分布するためで、 Otsu は谷ではなく肩を切っていた。加えて無変化ピークは 0 ではなく +0.066 にあり、 「変化なし = dNBR 0」という素朴な前提が季節差で崩れることも分かった。
② 分散曲線が平坦なら閾値は不安定 — 誤り
閾値の安定度を分散曲線の平坦さで測ろうとしたが、合成した綺麗な二山のほうが平坦域は広かった (0.453 対 0.133)。谷が空なら、どこで切っても同じ結果になるので平坦なのだ。 正しい判定は平坦域に画素が何個あるかで、空なら安定・詰まっていれば不安定。 実データは 9.0%(43.2 km²)が詰まっており、閾値は一点ではなく 43 km² の幅を持っていた。
③ 2022年の再焼失域を系統的に見逃す — 外れ
参照の報告書に「2022年に焼けた区域が再び焼けた」とあり、既に低 NBR の場所は再焼失しても dNBR が伸びない(床効果)と予測した。実測では見逃しの 1.1% しか説明できず、 支配的だったのは閾値の差(83.7%)だった。 理屈が正しいことと、その効果が主要因であることは別である。
第 4 層技術者として得た知見
1. この手法が効いたのは理想条件だったから
IoU 0.87 は手法の一般的な性能ではなく、条件が揃った結果である。 両時期とも雲量ほぼゼロ、同一 MGRS タイルで位置ずれなし、単一事象、 そして参照が同じ post 画像・同系統の手法で作られている。
崩れる条件も具体的に見えた。季節差は無変化ピークを 0 から +0.066 ずらし、 分布を太らせて閾値を不安定にする。再焼失は床効果で検出感度を下げる。 雲があれば有効画素が減り、差分の母集団自体が欠ける。 そして参照定義が変われば同じ予測でも IoU は 0.56〜0.87 に動く。
根拠: B0 の雲量選定 / B3 の無変化ピーク +0.0656 と平坦域 43.2 km² / B5 の参照定義3通りの比較。
2. 検証における一致と不一致は対称ではない — 偽の一致のほうが危険
ライブラリとの照合が合わないことが4回あり、すべてバグではなく定義・規約の差だった (小領域除去の境界規約、Otsu のビン設計依存、モルフォロジの画像端の扱い、SWIR バンドの選択)。 不一致は調査を強制するので、結果的に理解が深まる。
危険なのは逆で、B3 で「自動閾値 +0.447 が公式基準 +0.44 とほぼ一致」と書いたとき、 一致に見えたので調査を止めた。B5 で測り直すと、両者は別のバンドで計算された別の量で、 一致は偶然に近かった。数値が近いことと、同じものを測っていることは別である。
根拠: B2〜B5 の照合4件 / B5 で参照クラス別に測った dNBR の系統的なずれ。
3. では、どういう一致なら信じてよいか
「一致は信じるな」では使えない。実用的な基準は 一致の情報量 = 満たされた制約の数 − 調整した自由度の数。
ずれの原因を「参照は SWIR に B11 を使ったのではないか」と考え、 符号を先に宣言してから(B12 のほうが焼失後の反射上昇が大きいので dNBR も大きくなるはず) 測定した。B11 で再計算すると、両端が閉じた2クラスで四分位範囲がまるごと公称帯に収まった (平均逸脱 +0.169 → +0.000)。制約は6個以上、調整した自由度はゼロ(バンドを差し替えただけ)、 そして逆に出れば仮説は死んでいた。B3 の「スカラー1個が近い」とは質が違う。
それでも断定はしない。B11 版は B12 版のほぼ 0.65 倍に比例するので、 同じ倍率を生む別処理でも同じ結果は作れる。技術報告書にバンドの明記もない。
根拠: src/exp_swir_band.py の実測(下図)。
4. 選択は、測定できる時点まで遅らせる
後処理の設定を B4 で一つに決めず、4段階すべてを評価段へ渡した。おかげで 「s2 と s3 の差は参照自身の精度より小さく、区別できない」という情報が得られた。 勘で片方を選んでいたら、この事実自体が観測されない。
選択を早く固定するほど、その選択の当否を後から検証できなくなる。 測定可能になるまで複数の候補を保持したまま運ぶコストは、多くの場合それに見合う。
根拠: B4 の4段階 × B5 の3参照定義 = 12通りの評価。
補足限界と、次にやること
- 参照は真値ではない。 Overall Accuracy 96.33%、しかも同じ post 画像・同系統の手法で 作られているため、共通の系統誤差は指標に現れない。これは独立検証ではなく整合性チェックである。
- AOI は火災全体の一部。 火災は 923 km² だが、評価したのは約 481 km² の窓の中(参照 AOI との積で 92.2%)。
- SWIR バンドの推定は未確定。 一致にもとづく推論で、文書による裏づけはない。
- 次の検証候補: 2022年火災以前の pre で再実行し、見逃しが誤検知に付け替わる量を測る / B11 版 dNBR でパイプライン全体を通し直す / SAR(log-ratio)への転移。