衛星画像で軽石検出に挑戦第2弾〜機械学習で検出〜

はじめに

こんにちは、UMITRONの宮山（@jkatagi）です。先日Google Earth Engine（以下GEE）を用いて軽石の検出にトライしてみました。予想より多くの方に閲覧していただき、嬉しい限りです。

上記記事ではかなり雑に閾値で区切って軽石の領域を検出しましたが、今回は別の手法として機械学習の手法である教師あり分類で検出にトライしてみようと思います。

なお今回も前回同様にSentinel-2を用いますが、JAXAの衛星であるGCOM-C/SGLIでも軽石を確認できるとのことです。GCOM-C/SGLIは空間解像度（1ピクセルあたりの大きさ）が250 mとSentinel-2の10 mに比べて荒いですが、その分広域を高頻度で観測できます。従っていつ到達したかを検知するには、そのような衛星を用いた方が望ましいかもしれません。

衛星「しきさい」(GCOM-C)による軽石観測画像について

教師あり分類とは

教師あり分類は機械学習の手法の一種で、人間が判読した教師情報をもとにアルゴリズムを用いて分類するものです。今回はこの教師あり分類のアルゴリズムの一つであるランダムフォレストを用いました（ランダムフォレストについては検索すると詳しい説明が出てくるので割愛します）。

分類条件

今回は以下のような条件で分類することにしました。

・あらかじめ陸域をマスク（除外）する
・クラスは4つ（海域・雲・雲影・軽石）
・教師情報は各クラス60点（以上）
・アルゴリズムはランダムフォレスト

あらかじめ陸域をマスクすることによって前回の記事のような陸域を軽石判定することがなくなります。またクラスはパッと思いついた上記4クラスにしました。教師情報の数は一旦60点と適当に決め、分類後足りないようであれば適宜増やす方針としました。アルゴリズムをランダムフォレストにしたのは、GEEに実装されているからです。

陸域をマスクする

陸域のマスクのための情報としてはMODIS（NASAのTerra/Aqua衛星搭載のセンサー）とSRTM（レーダーのミッション）から作成された MOD44W.005というプロダクトを用いました：

こちらは解像度250 mと少々荒いですが、ざっくり陸域を削れればいいと判断したため選択しました（本当はもう少し細かい Natural Earthというのがあるのですが、GeoJSONのimportに手こずったために断念しました）。なお、この後の画像の陸域のふちが若干斜めになっていますが、これはreprojectionの処理を省いたためです（繰り返しになりますがざっくりと削りたかったため省略）。

教師情報を取得する

次に教師情報を取得していきます。取り方としてはSentinel-2の画像を見ながら、ここは海域、ここは雲、ここは軽石といったようにポチポチしていきます。ツールはマーカー（Geometry）でFeatureCollectionとして保存しました。

水色のマーカーが海域、黄土色が軽石、白色が雲、濃い茶色が雲影です。

分類結果

その後、学習をし分類した結果が以下です：

分類後のマップでは軽石を目立たせるために赤色にしています（その他は青が海域、白が雲、黒が雲影）。これをみるとしっかりと軽石の領域が抽出できている一方で（※教師情報を取ったところ以外を眺めて判断）、雲の周りが軽石と誤分類されていることがわかります。これは以下二つの要因が考えられます。一つ目は単に教師情報が少ないことです。二つ目は軽石のスペクトル情報が雲または雲付近の海域と似ている可能性です。一つ目は単に教師情報を増やす（軽石/雲・軽石/海域の判断に迷う地点の教師を増やすとより良さそう）ことで解決しますが、二つ目はそれでは難しいと予想されます。そのため二つ目の対策として、モルフォロジーによる領域拡張を検討しました。

↑分類結果を透かしながら教師情報を増やす作業

最終的には以下のようなマップができました：

↓がモルフォロジー適用前で、比べてみると雲の周りの誤検出がだいぶ削減できました。

その一方で軽石領域まで雲判定されることも多くなってしまったため、一長一短かなと思います。

↑モルフォロジー適用後

↑モルフォロジー適用前

なお、GEEのすごいところはズームアウトするとその領域も同様の分類処理をして表示してくれるところです。なお沖縄周辺のみで学習したため、精度は落ちると思われます（左上は真っ赤になってしまっていますね…）。四国沖でも軽石の確認がされたとのことで全国での確認をしたいのですが、まだ検討が必要そうです。

なお今回作成したコードは以下のようになっています。教師情報はGUI上で取得したため、下記コードをそのまま走らせても走りませんのでその点ご留意ください。

function maskS2clouds(image, oceanImage) {
 var qa = image.select('QA60');
 var blue = image.select('B2');
 var green = image.select('B3');
 var red = image.select('B4');
 var b9 = image.select('B9');
 

 // Bits 10 and 11 are clouds and cirrus, respectively.
 var cloudBitMask = 1 << 10;
 var cirrusBitMask = 1 << 11;

 // Both flags should be set to zero, indicating clear conditions.
 var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
     .and(qa.bitwiseAnd(cirrusBitMask).eq(0));

 return image.updateMask(mask).divide(10000);
}

var sentinelImageCollection = ee.ImageCollection('COPERNICUS/S2_SR');
var oceanImage = ee.Image("MODIS/MOD44W/MOD44W_005_2000_02_24");
var oceanMask = oceanImage.select('water_mask').eq(1); // 1 is water

var startDate = '2021-10-25';
var endDate = '2021-10-30';

// roi is region of interest
var roi = ee.Geometry.Rectangle(127.61094820401469, 26.04176343366477, 128.3604519941186, 27.038911596761107)

var bands = ["B1", "B2", "B3", "B4", "B5", "B6", "B7", "B8", "B8A", "B9", "WVP"];

var sentinelImage = sentinelImageCollection
                 .filterDate(startDate, endDate)
                 // Pre-filter to get less cloudy granules.
                 .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE',20)) // 20
                 .map(maskS2clouds);

var maskedSentinelImage = sentinelImage.map(function(image) { 
 return image.updateMask(oceanMask)})
 .median();
 
// classifiy
var traningPoints = water.merge(cloudShade).merge(cloud).merge(pumice);
var traning = maskedSentinelImage.select(bands).sampleRegions({
 collection: traningPoints,
 properties: ['value'],
 scale: 10
})
var classifier = ee.Classifier.smileRandomForest(5).train({
 features: traning,
 classProperty: 'value',
 inputProperties: bands,
})
var classified = maskedSentinelImage.select(bands).classify(classifier);

// remove cloud
var cloudValue = 0;
var cloud = classified
 .select("classification")
 .eq(cloudValue);
var cloudMask = cloud.focal_max(1);
var maskedClassified = classified.where(cloudMask, cloudValue);


// viz parameter settings
var sentinelVizParams = {
 min: 0.0,
 max: 0.3,
 bands: ['B4', 'B3', 'B2'],
}
var classifiedVizParams = {
 min: 0,
 max: 3,
 // cloud, water, cloud shade, pumice
 palette: ['white', 'blue', 'black', 'red'],
}

// okinawa
Map.setCenter(128.0768285144796,26.703266227408598, 12);
Map.addLayer(maskedSentinelImage, sentinelVizParams, "sentinel2");
Map.addLayer(classified, classifiedVizParams, "clasified");
Map.addLayer(maskedClassified, classifiedVizParams, "maskedClassified");

おわりに

今回は教師あり分類によって軽石検出に挑戦しました。教師の取り方や使うバンドの選択などによって精度は変化すると思います。また光学画像は雲の影響が大きいため、光学以外の活用も期待されます。実際にJAXAでは合成開口レーダによって軽石の検出を試みているようです。