アマゾンセイジメーカーの地理情報能力を使用したメタン排出ポイント源の検出と高周波監視

アマゾンセイジメーカーによる地理情報能力を活用したメタン排出ポイント源の検出と高周波監視

メタン（CH4）は、石油やガスの抽出、石炭採掘、大規模な動物飼育、廃棄物処理など、他のさまざまな源からの副産物である主要な人為的温室効果ガスです。CH4の地球温暖化潜在能はCO2の86倍であり、気候変動に関する政府間パネル（IPCC）によると、メタンはこれまでの観測された地球温暖化の30％を引き起こしていると推定されています。大気中へのCH4の漏洩を急速に減らすことは、気候変動との闘いにおいて重要な要素をなしています。2021年には、国連が気候変動会議（COP26）で「グローバルメタン宣言」を発表し、「1.5℃の未来を実現するためのメタンへの迅速な対策を行う」ことを目指しています。この宣言には、米国やEUなど150以上の署名国が参加しています。

メタン源の早期発見と継続的なモニタリングは、メタンに関する意味のある対策の重要な要素であり、政策立案者や組織の関心事となっています。現場で使用することができる手頃な価格で効果的なメタン検出ソリューション（例：現場用メタン検出器または機内に取り付けられた分光計）を実施することは、しばしば非現実的または高額です。一方、衛星を使用した遠隔センシングは、関係者が望むような地球規模、高頻度、費用対効果の高い検出機能を提供することができます。

このブログポストでは、AWS Registry of Open DataにホストされるSentinel 2衛星画像を使用し、Amazon SageMakerの地理空間能力と組み合わせて、CH4排出源のポイントを検出し、時間とともにモニタリングする方法を紹介します。地球観測文献の最新の研究成果を活用しながら、世界中のさまざまな場所からのメタン漏洩を検出およびモニタリングするためのカスタムメタン検出アルゴリズムを実装する方法を学びます。この投稿には、GitHubの関連するコードが含まれており、追加の技術的な詳細を提供し、独自のメタンモニタリングソリューションを始めるのに役立ちます。

従来、複雑な地理空間分析の実行は困難で、時間がかかり、リソースを消費する作業でした。Amazon SageMakerの地理空間能力は、データサイエンティストや機械学習エンジニアが地理空間データを使用してモデルを構築、トレーニング、展開することを容易にします。SageMakerの地理空間能力を使用すると、大規模な地理空間データセットの変換や拡張を効率的に行うことができます。事前に学習済みの機械学習（ML）モデルを使用してモデルの構築を加速し、3Dアクセラレーショングラフィックスと組み込みの可視化ツールを使用してモデルの予測と地理空間データをインタラクティブな地図上で探索することができます。

• ニューラルネットワークは人間のような言語一般化を達成する
• 「大規模な言語モデルが医療テキスト分析に与える影響」
• 「ReactJSとChatGPT：ウェブアプリに対話型AIを構築する方法」

多光谱衛星画像を使用したメタンポイント源の遠隔センシング

衛星を利用したメタンセンシング手法では、CH4の固有の透過特性に依存することが一般的です。可視スペクトルでは、CH4は透過値が1と等しいか非常に近いため、肉眼では検出できません。しかし、一部の特定の波長では、メタンは光を吸収します（透過値<1）。この特性は検出のために利用されます。そのため、短波長赤外（SWIR）スペクトル（1500〜2500 nmの波長範囲）が通常選ばれ、そこでCH4が最も検出されます。ハイパースペクトルおよび多光谱衛星ミッション（すなわち、電磁スペクトル全体で複数の波長範囲（バンド）の画像データをキャプチャする光学機器を持つもの）は、このSWIR範囲をカバーし、検出に有効な手段となります。図1はSWIRスペクトルでのメタンの透過特性と各候補の多光谱衛星機器のSWIR範囲を示しています（こちらの研究を参考に作成）。

図1- SWIRスペクトルにおけるメタンの透過特性とSentinel-2マルチスペクトルミッションのカバー範囲

多くの多スペクトル衛星ミッションは、再訪頻度が低い（例：PRISMAハイパースペクトルは約16日）または空間分解能が低い（例：Sentinel 5は7.5 km x 7.5 km）など、制約があります。データへのアクセスコストは追加の課題です：一部の専用コンステレーションは商業ミッションとして運用されており、資金の制約により、研究者、意思決定者、およびその他の関係者へのCH4排出に関する洞察が研究の公開性を十分に行き渡らなくしています。ESAのSentinel-2マルチスペクトルミッションは、再訪率（約5日）、空間分解能（約20 m）およびオープンアクセス（AWS Registry of Open Dataにホストされます）の適切なバランスを実現しています。

Sentinel-2にはSWIRスペクトルをカバーする2つのバンド（20 mの解像度）があります：バンド11（1610 nmの中心波長）とバンド12（2190 nmの中心波長）。両バンドはメタンの検出に適しており、バンド12はCH4吸収に対して非常に高い感度を持っています（図1参照）。直感的には、メタン検出のためにこのSWIR反射データを使用する場合、2つのアプローチが考えられます。1つ目は、異なる2つの衛星パス間でピクセルごとの反射率の差を計算するだけである（理想的には、CH4吸収に最も敏感なバンド）。または、表面およびエアロゾルの反射特性は類似していますが、メタン吸収特性が異なる2つの隣接するスペクトルSWIRバンドからのデータを使用して、検出のために単一の衛星パスのデータを使用することもできます。

このブログ記事では、両方のアプローチを組み合わせた検出方法を実装します。地球観測文献からの最近の研究成果を活用し、2つの衛星パスおよび2つのSWIRバンド間の大気上端（TOA）反射率Δρ（つまり、Sentinel-2によって測定される反射率（大気エアロゾルとガスの寄与を含む））の分数の変化を計算します。一方のベースラインパスではメタンが存在しない（base）、もう一方のモニタリングパスではアクティブなメタンポイントソースが疑われる場合（monitor）。数学的には、以下のように表現できます：

式（1）

ここで、ρはSentinel-2によって測定されるTOA反射率であり、cmonitorおよびcbaseはシーン全体でのバンド11のTOA反射率値とバンド12のTOA反射率値の回帰によって計算されます（つまり、ρb11 = c * ρb12）。詳細については、多スペクトルSentinel-2衛星観測による異常メタンポイントソースの高頻度モニタリングに関するこの研究を参照してください。

SageMakerジオスペーシャル機能を使用したメタン検出アルゴリズムの実装

メタン検出アルゴリズムを実装するために、Amazon SageMaker Studio内のSageMakerジオスペーシャルノートブックを使用します。ジオスペーシャルノートブックカーネルには、GDAL、GeoPandas、Shapely、xarray、およびRasterioなどの必要なジオスペーシャルライブラリが予め装備されており、Pythonノートブック環境でジオスペーシャルデータを直接視覚化および処理することができます。SageMakerジオスペーシャル機能の使用を開始する方法については、はじめにガイドを参照してください。

SageMakerはAPIを提供し、SearchRasterDataCollection API呼び出しを使用して、衛星画像の取得を容易にするために設計されたものです。 SearchRasterDataCollectionは、次の入力パラメータに依存します：

• Arn：クエリされたラスターデータコレクションのAmazonリソース名（ARN）
• AreaOfInterest：検索クエリの興味領域を表すポリゴンオブジェクト（GeoJSON形式）
• TimeRangeFilter：興味のある時間範囲を定義します。 {StartTime: <string>, EndTime: <string>}
• PropertyFilters：最大許容クラウドカバーなどの補足的なプロパティフィルタも組み込むことができます

このメソッドは、ListRasterDataCollectionsを呼び出すことで探索できるさまざまなラスターデータソースのクエリをサポートしています。当社のメタン検出実装では、Sentinel-2衛星画像が使用され、次のARNを使用してグローバルに参照できます：arn:aws:sagemaker-geospatial:us-west-2:378778860802:raster-data-collection/public/nmqj48dcu3g7ayw8。

このARNはSentinel-2のイメージを表しており、レベル2A（表面反射率、大気補正）の処理が行われています。メタン検出の目的では、大気上層部（TOA）の反射率データ（レベル1C）を使用しますが、これにはエアロゾル組成と密度の変化（つまり、メタン漏れ）を検出できない表面レベルの大気補正は含まれていません。

特定のポイントソースからの潜在的な排出物を特定するには、2つの入力パラメータが必要です：疑わしいポイントソースの座標とメタン排出物の監視のための指定されたタイムスタンプです。 SearchRasterDataCollection APIはポリゴンまたはマルチポリゴンを使用して興味領域（AOI）を定義するため、私たちのアプローチはまずポイント座標をバウンディングボックスに展開し、それからそのポリゴンを使用してSearchRasterDateCollectionを使用してSentinel-2画像をクエリすることです。

この例では、北アフリカの油田から起因する既知のメタン漏れを監視しています。これはリモートセンシング文献の標準的な検証ケースであり、たとえばこの研究で参照されています。実行可能なコードベースはamazon-sagemaker-examples GitHubリポジトリで提供されています。ここでは、SageMakerの地理情報能力を使用したメタン検出ソリューションの実装のキービルディングブロックを示す選択されたコードセクションのみを強調表示します。詳細についてはリポジトリを参照してください。

まず、例のケースの座標とターゲットの監視日を初期化します。

#北アフリカの油田の座標と日付#参考にはここをご覧ください：https://doi.org/10.5194/amt-14-2771-2021    point_longitude = 5.9053    point_latitude = 31.6585    target_date = '2019-11-20'#ポイントの周りのバウンディングボックスのサイズ    distance_offset_meters = 1500

次のコードスニペットは、指定されたポイント座標に対するバウンディングボックスを生成し、バウンディングボックスと指定された監視日に基づいて利用可能なSentinel-2イメージを検索します：

def bbox_around_point(lon, lat, distance_offset_meters):    #Equatorial radius (km) taken from https://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html    earth_radius_meters = 6378137    lat_offset = math.degrees(distance_offset_meters / earth_radius_meters)    lon_offset = math.degrees(distance_offset_meters / (earth_radius_meters * math.cos(math.radians(lat))))    return geometry.Polygon([        [lon - lon_offset, lat - lat_offset],        [lon - lon_offset, lat + lat_offset],        [lon + lon_offset, lat + lat_offset],        [lon + lon_offset, lat - lat_offset],        [lon - lon_offset, lat - lat_offset],    ])#バウンディングボックスを生成してポリゴン座標を抽出するaoi_geometry = bbox_around_point(point_longitude, point_latitude, distance_offset_meters)aoi_polygon_coordinates = geometry.mapping(aoi_geometry)['coordinates']#検索パラメータを設定するsearch_params = {    "Arn": "arn:aws:sagemaker-geospatial:us-west-2:378778860802:raster-data-collection/public/nmqj48dcu3g7ayw8", # Sentinel-2 L2データ    "RasterDataCollectionQuery": {        "AreaOfInterest": {            "AreaOfInterestGeometry": {                "PolygonGeometry": {                    "Coordinates": aoi_polygon_coordinates                }            }        },        "TimeRangeFilter": {            "StartTime": "{}T00:00:00Z".format(as_iso_date(target_date)),            "EndTime": "{}T23:59:59Z".format(as_iso_date(target_date))        }    },}#SageMaker地理情報能力を使用してラスターデータをクエリするsentinel2_items = geospatial_client.search_raster_data_collection(**search_params)

応答には、一致するSentinel-2アイテムとその対応するメタデータのリストが含まれています。これには、すべてのSentinel-2バンドのクラウド最適化ジオTIFF（COG）、およびイメージのビジュアルバンドのクイックプレビュー用のサムネイル画像が含まれます。もちろん、図2に示すように、フル解像度の衛星画像（RGBプロット）にもアクセスできます。

図2 – AOIの衛星画像（RGBプロット）

以前に詳細に説明したように、検出手法は大気上層（TOA）SWIR反射率の部分的な変化に依存しています。これを実現するには、適切なベースラインの同定が重要です。適切なベースラインを見つけることは、試行錯誤を伴う手間のかかるプロセスになることがよくありますが、良いヒューリスティクスを使用するとこの検索プロセスを自動化することができます。過去の調査事例でうまく機能した検索ヒューリスティクスは次のとおりです：過去のday_offset=n日間の衛星画像を取得し、クラウドを除去し、対象のAOIに画像をクリップします。そして、AOI全体でバンド12の反射率の平均値を計算します。バンド12の平均反射率が最も高い画像のSentinelタイルIDを返します。

この論理は以下のコード抜粋に実装されています。その根拠は、バンド12がCH4吸収に高い感度を持つことです（図1を参照）。平均反射率がより高い値は、メタン排出などのソースからの吸収が低いことを示し、排出物のないベースラインシーンを示す強力な指標となります。

def approximate_best_reference_date(lon, lat, date_to_monitor, distance_offset=1500, cloud_mask=True, day_offset=30):        #initialize AOI and other parameters    aoi_geometry = bbox_around_point(lon, lat, distance_offset)    BAND_12_SWIR22 = "B12"    max_mean_swir = None    ref_s2_tile_id = None    ref_target_date = date_to_monitor               #loop over n=day_offset previous days    for day_delta in range(-1 * day_offset, 0):        date_time_obj = datetime.strptime(date_to_monitor, '%Y-%m-%d')        target_date = (date_time_obj + timedelta(days=day_delta)).strftime('%Y-%m-%d')                   #get Sentinel-2 tiles for current date        s2_tiles_for_target_date = get_sentinel2_meta_data(target_date, aoi_geometry)                #loop over available tiles for current date        for s2_tile_meta in s2_tiles_for_target_date:            s2_tile_id_to_test = s2_tile_meta['Id']            #retrieve cloud-masked (optional) L1C band 12            target_band_data = get_s2l1c_band_data_xarray(s2_tile_id_to_test, BAND_12_SWIR22, clip_geometry=aoi_geometry, cloud_mask=cloud_mask)            #compute mean reflectance of SWIR band            mean_swir = target_band_data.sum() / target_band_data.count()                        #ensure the visible/non-clouded area is adequately large            visible_area_ratio = target_band_data.count() / (target_band_data.shape[1] * target_band_data.shape[2])            if visible_area_ratio <= 0.7: #<-- ensure acceptable cloud cover                continue                        #update maximum ref_s2_tile_id and ref_target_date if applicable            if max_mean_swir is None or mean_swir > max_mean_swir:                 max_mean_swir = mean_swir                ref_s2_tile_id = s2_tile_id_to_test                ref_target_date = target_date    return (ref_s2_tile_id, ref_target_date)

この方法を使用すると、適当なベースラインの日付と対応するSentinel-2タイルIDを概算することができます。Sentinel-2タイルIDには、ミッションID（Sentinel-2A / Sentinel-2B）、一意のタイル番号（例：32SKA）、および画像が撮影された日付などの情報が含まれており、観測を一意に識別します。この例では、概算プロセスにより、2019年10月6日（Sentinel-2タイル：S2B_32SKA_20191006_0_L2A）が最も適したベースライン候補として示されます。

次に、基準日とモニタリングしたい日の反射率の修正された変化率を計算することができます。修正係数 c（前述の式 1 を参照）は以下のコードで計算されます：

def compute_correction_factor(tif_y, tif_x):        # 回帰のためにフラット化した配列を取得するために最少二乗法を使用する   y = np.array(tif_y.values.flatten())    x = np.array(tif_x.values.flatten())    np.nan_to_num(y, copy=False)    np.nan_to_num(x, copy=False)            # 最少二乗法を使用して線形モデルにフィットする    x = x[:,np.newaxis] #reshape    c, _, _, _ = np.linalg.lstsq(x, y, rcond=None)    return c[0]

式 1 の完全な実装は以下のコードスニペットで与えられています：

def compute_corrected_fractional_reflectance_change(l1_b11_base, l1_b12_base, l1_b11_monitor, l1_b12_monitor):        # 修正係数を取得する    c_monitor = compute_correction_factor(tif_y=l1_b11_monitor, tif_x=l1_b12_monitor)    c_base = compute_correction_factor(tif_y=l1_b11_base, tif_x=l1_b12_base)        # 修正された分光反射率の変化を取得する    frac_change = ((c_monitor*l1_b12_monitor-l1_b11_monitor)/l1_b11_monitor)-((c_base*l1_b12_base-l1_b11_base)/l1_b11_base)    return frac_change

最後に、上記のメソッドをエンドツーエンドの手順に組み込んで、指定した経度と緯度、モニタリング日、基準タイルでの AOI を特定し、必要な衛星画像を取得し、分光反射率の変化を計算します。

def run_full_fractional_reflectance_change_routine(lon, lat, date_monitor, baseline_s2_tile_id, distance_offset=1500, cloud_mask=True):        # バウンディングボックスを取得する    aoi_geometry = bbox_around_point(lon, lat, distance_offset)        # S2 のメタデータを取得する    s2_meta_monitor = get_sentinel2_meta_data(date_monitor, aoi_geometry)        # タイル ID を取得する    grid_id = baseline_s2_tile_id.split("_")[1]    s2_tile_id_monitor = list(filter(lambda x: f"_{grid_id}_" in x["Id"], s2_meta_monitor))[0]["Id"]        # センチネル L1C プロダクトのバンド 11 と 12 を取得する    l1_swir16_b11_base = get_s2l1c_band_data_xarray(baseline_s2_tile_id, BAND_11_SWIR16, clip_geometry=aoi_geometry, cloud_mask=cloud_mask)    l1_swir22_b12_base = get_s2l1c_band_data_xarray(baseline_s2_tile_id, BAND_12_SWIR22, clip_geometry=aoi_geometry, cloud_mask=cloud_mask)    l1_swir16_b11_monitor = get_s2l1c_band_data_xarray(s2_tile_id_monitor, BAND_11_SWIR16, clip_geometry=aoi_geometry, cloud_mask=cloud_mask)    l1_swir22_b12_monitor = get_s2l1c_band_data_xarray(s2_tile_id_monitor, BAND_12_SWIR22, clip_geometry=aoi_geometry, cloud_mask=cloud_mask)        # 修正された分光反射率の変化を計算する    frac_change = compute_corrected_fractional_reflectance_change(        l1_swir16_b11_base,        l1_swir22_b12_base,        l1_swir16_b11_monitor,        l1_swir22_b12_monitor    )        return frac_change

前に決定したパラメータでこのメソッドを実行すると、SWIR TOA 反射率の分数変化が xarray.DataArray として得られます。このデータ配列に対して plot() を実行することで結果の視覚的なインスペクションが行えます。このメソッドは、以前に見た RGB プロットでは検出できなかった AOI の中心にメタンプルームの存在を示しています。

図3 – TOA 反射率の分数変化（SWIR スペクトル）

最後のステップとして、特定されたメタンプルームを抽出し、生のRGB衛星画像にオーバーレイして重要な地理的コンテキストを提供します。これは、次に示すように実装されるしきい値によって実現されます：

def get_plume_mask(change_in_reflectance_tif, threshold_value):    cr_masked = change_in_reflectance_tif.copy()    #しきい値を上回る値をnan値に設定する    cr_masked[cr_masked > treshold_value] = np.nan     #nan値にマスクを適用する    plume_tif = np.ma.array(cr_masked, mask=cr_masked==np.nan)         return plume_tif

この場合、反射率の変化が-0.02のしきい値を得ると、良好な結果が得られますが、シーンごとに異なる可能性があり、特定のユースケースに合わせてキャリブレーションする必要があります。次の図4は、AOI（Area of Interest）の生の衛星画像とマスクされたプルームを組み合わせて、地理的なコンテキストでメタンプルームを表示する単一の複合イメージを生成する方法を示しています。

図4 – RGB画像、TOA反射率（SWIRスペクトル）の分散反射率変化、およびAOIのメタンプルームオーバーレイ

実世界のメタン排出イベントによる解決策の検証

最後のステップとして、さまざまなソースと地理的な場所からのメタン漏洩を正確に検出し特定する能力について、メソッドを評価します。まず、メタン漏洩のスペースベースの点源検出と陸上メタン排出量の定量化の検証に特化した制御メタン放出実験を使用します。この2021年の実験では、研究者が19日間にわたってアリゾナ州エーレンバーグで複数のメタン放出を行いました。その実験の期間中にSentinel-2パスの1つで実行すると、メタンプルームを示す次の結果が得られます：

図5 – アリゾナ制御放出実験のメタンプルームの濃度

制御放出中に生成されたプルームは、私たちの検出方法によって明確に識別されます。同様の結果が、東アジアの埋立地（左）や北米の油田・ガス施設（右）などの他の既知の実世界の漏洩（以下の図6）でも見られます。

図6 – 東アジアの埋立地（左）および北米の油田・ガス施設（右）のメタンプルームの濃度

まとめると、私たちの方法は、制御放出および世界中のさまざまなポイントソースからのメタン排出を識別するのに役立ちます。これは、周囲の植生が制限されている陸上のポイントソースに最適です。水によるSWIRスペクトルの吸収（つまり透過率の低さ）のため、沖合のシーンでは機能しません。提案された検出アルゴリズムはメタン強度の変動に依存しているため、リークの定常排出率のモニタリングは困難になる場合があります。

クリーンアップ

メタンモニタリングのジョブが完了した後、不要な料金が発生しないように、SageMakerインスタンスを終了し、不要なローカルファイルを削除してください。

結論

オープンジオスペーシャルデータソースとSageMakerの地理空間能力を組み合わせることで、スケールに対応した独自の高度にカスタマイズされたリモートモニタリングソリューションを実装することができます。このブログ投稿では、政府、NGO、他の組織が有害なメタン排出を検出し、最終的に回避するためにメタン検出に焦点を当てています。SageMakerの地理空間カーネルを使用してノートブックを起動し、独自の検出ソリューションを実装することで、今日から地理空間解析の旅に出ることができます。自身の衛星ベースのメタン検出ソリューションを構築するためには、GitHubリポジトリを参照してください。また、他の現実世界のリモートセンシングアプリケーションでSageMakerの地理空間能力を使用する方法に関するさらなる例やチュートリアルについては、sagemaker-examplesリポジトリもご覧ください。

We will continue to update VoAGI; if you have any questions or suggestions, please contact us!