EP.36 で、衛星データが波長ごとの測定値であることを扱いました。この回はその応用として、 を扱います。
同じ場所を違う時点で観測し、差から変化を見つける。造成、伐採、被災、作付けの変化。原理は単純ですが、そのまま引き算すると誤検出だらけになります。
1. 何が誤検出を生むか
2時点の画像を引き算したとき、変化していないのに差が出る要因を先に把握します。
| 要因 | 何が起きるか | 対処 |
|---|---|---|
| 季節 | 植生の状態が違う | 同じ季節同士で比べる |
| 太陽の角度 | 影の出方が変わる | 近い時刻・時期を選ぶ |
| 雲・霞 | 明るさが変わる | で除外 |
| 観測角度 | 見え方が変わる | 同じ軌道の観測を選ぶ |
| 大気の状態 | 全体の値がずれる | 補正済みのデータを使う |
| 位置のずれ | 別の場所を比べてしまう | 位置合わせを確認 |
1行目が最も影響が大きい。春と秋を比べれば、何も変わっていなくても植生の値は大きく違います。「前年の同じ時期」と比べるのが基本になります。
6行目は見落とされやすいのですが、数十メートルずれているだけで、まったく別の場所を比較することになります。特に境界（建物の縁、道路沿い）で大きな差が出て、変化があったように見えます。提供されているデータは通常補正済みですが、確認する価値はあります。
2. 差を取る前に条件を揃える
実務の手順は、比較する2時点を選ぶところから始まります。ここで大半が決まります。
- 1同じ季節の、近い時期を選ぶ — 前年同月が基本
- 2雲の少ない観測を選ぶ — 雲被覆率で絞る
- 3同じ提供元・同じ処理段階を使う — 混ぜない
- 4同じ範囲に切り出す — 位置を揃える
- 5雲マスクで、雲の画素を両方から除外する — 片方だけでは不十分
5番目が重要です。片方に雲があれば、その画素は両方とも使えません。「時点Aは晴れているが時点Bは雲」という画素を比較すると、雲の明るさが変化として出ます。
import numpy as np
def change_index(before_nir, before_red, after_nir, after_red, before_cloud, after_cloud): """2時点の植生指標の差を出す。雲のある画素は除外する。
*_cloud: True が雲（提供元の雲マスクから作る） """ def ndvi(nir, red): denom = nir + red # 0除算を避ける（分母が0の画素は無効にする） return np.where(denom == 0, np.nan, (nir - red) / np.where(denom == 0, 1, denom))
v_before = ndvi(before_nir, before_red) v_after = ndvi(after_nir, after_red)
# どちらか一方でも雲なら、その画素は使わない invalid = before_cloud | after_cloud | np.isnan(v_before) | np.isnan(v_after)
diff = v_after - v_before diff[invalid] = np.nan
valid_rate = (~invalid).mean() return diff, valid_rate
# 使える画素の割合が低ければ、その比較は信頼できない# diff, rate = change_index(...)# if rate < 0.5:# print(f"使える画素が {rate:.0%} しかありません。別の日付を検討してください")使える画素の割合を返しているのが要点です。雲を除外した結果、半分も残らないなら、その比較から結論を出すのは危険です。判断できないことを判断できるようにしておきます。 のような指標を計算する前に、そもそも計算できる画素がどれだけあるかを確認する、という順序です。
3. どこからを「変化」とするか
差を取っても、どこからを変化とみなすかの閾値が要ります。ここは対象によって変わります。
import numpy as np
def suggest_threshold(diff, n_sigma=2.0): """差の分布から、変化とみなす閾値の目安を出す。
大半の画素は変化していないという前提で、 分布の中心から離れたものを変化候補とする。 """ valid = diff[~np.isnan(diff)] if valid.size == 0: return None
# 外れ値に強い指標を使う（一部の大きな変化に引きずられないように） median = np.median(valid) mad = np.median(np.abs(valid - median)) spread = mad * 1.4826
return { "中央値": median, "散らばり": spread, "減少の閾値": median - n_sigma * spread, "増加の閾値": median + n_sigma * spread, "変化候補の割合": float( np.mean(np.abs(valid - median) > n_sigma * spread) ), }
# 実データで確認するrng = np.random.default_rng(0)# 大半は変化なし（0付近）、5%が明確に減少（伐採を想定）diff = np.concatenate([ rng.normal(0, 0.03, 9500), rng.normal(-0.35, 0.05, 500),])for k, v in suggest_threshold(diff).items(): print(f" {k:14s}: {v:+.4f}" if isinstance(v, float) else f" {k}: {v}")大半の画素は変化していないという前提を使っています。分布の中心付近が「変化なし」で、そこから離れたものが変化候補。外れ値に強い指標を使うのは、一部の大きな変化に閾値が引きずられないためです。
ただし、この前提が成り立たない場合もあります。広範囲が一斉に変わる（災害、大規模な造成）ときは、「変化していない画素」が少数派になります。その場合は、別の場所を基準にするなどの工夫が要ります。
もう1つ、この方法は必ず過剰に検出します。分布の裾を切る以上、変化していない画素の一定割合が候補に入るためです。実際に試すと、本当の変化が5%のデータで、候補は8%前後になります。差の3%は誤検出です。EP.37 の用途（候補を絞る）では許容できますが、候補の件数をそのまま「変化した面積」として報告してはいけません。
4. 面としてまとめる
画素単位の変化は、ばらついて点在します。実務で意味があるのは、まとまった面としての変化です。
- 孤立した画素を除く — 1画素だけの変化はノイズの可能性が高い
- 隣接する画素をまとめる — 面として扱う
- 面積で絞る — 一定以上の広さのものだけを対象に
- 形を見る — 直線的な境界は人為的な変化を示唆
3番目が実務的な判断です。何平方メートル以上を対象とするかを決めると、検出結果の件数が現実的な数に収まります。決めないと、数万件の変化候補が出て確認できません。
1画素が10メートル四方なら、1画素は100平方メートルです。それより小さい変化は原理的に捉えられません。使うデータの解像度から、検出できる下限が決まるということを、最初に確認してください。
5. 検出したあとの確認
変化を検出しても、それが何かは分かりません。EP.36 で書いたとおり、衛星データは測定値であって、意味は付いていません。
| 検出された変化 | 考えられる原因 |
|---|---|
| 植生の値が大きく減少 | 伐採、造成、災害、収穫、季節のずれ |
| 植生の値が増加 | 植林、作付け、回復、季節のずれ |
| 水面らしき反応が拡大 | 浸水、貯水、新しい施設 |
| 明るさだけが変化 | 大気や観測条件の可能性 |
どの行にも「変化ではない可能性」が含まれている点に注意してください。特に最下行は、そもそも地表が変わっていない可能性が高い。
確認の手段としては、他のデータと突き合わせるのが有効です。建築確認の情報、土地利用の記録、報道。EP.35 で扱った地図データと重ねると、その場所に何があるかが分かります。ただし次回以降で扱うとおり、地図データにも網羅率の問題があるため、そちらも万能ではありません。複数の情報を重ねて確からしさを上げる、という進め方になります。
6. 何に使えるか
最後に、実務での使いどころを整理します。個別の判断材料ではなく、広域の把握に向きます。
| 用途 | 向き | 理由 |
|---|---|---|
| 広域の土地利用の変化を把握 | 向く | 現地調査では網羅できない |
| 災害後の被害範囲の概観 | 向く | 早期に広範囲を確認できる |
| 造成・開発の兆候を見つける | 向く | 定期的に確認できる |
| 個別の物件の状態を判断 | 向かない | 解像度が足りない |
| 変化の理由を特定 | 向かない | 他の情報が必要 |
| 厳密な面積の算定 | 向かない | 誤差が大きい |
上3つは「候補を絞る」用途です。全域を人が確認するのは不可能なので、衛星で候補を出して、そこだけ詳しく見る。この使い方なら、多少の誤検出があっても実務が回ります。見逃しは困るが、誤検出は確認すれば済むという非対称があるため、やや広めに検出する設定が実務的です。
逆に、下3つのように個別の判断や厳密な数値を出す用途には向きません。誤検出を許容できるかどうかが、使えるかどうかの分かれ目です。
そのまま引き算すると誤検出だらけになる。最も効くのは同じ季節同士で比べること。雲は両時点で除外し（片方だけでは不十分）、使える画素の割合を確認する。閾値は大半が変化していないという前提で分布から決める。面としてまとめて面積で絞らないと件数が多すぎる。そして変化が何かは分からないので、他のデータと突き合わせてください。
この記事の感想を教えてください
あなたの 1 クリックで、本当にこの記事は更新されます。「もっと詳しく」「続編希望」が一定数集まった記事は、 ふくふくが 実際に内容を拡充したり続編記事を公開 します。 送信したリアクションはお使いのブラウザに記録され、再カウントされません。