使用 Amazon SageMaker 地理空间功能检测和高频监测甲烷排放点源

甲烷（CH4）是一种主要的人为温室气体，是石油和天然气开采、煤炭开采、大规模动物养殖和废物处理等来源的副产品。甲 烷的全球变暖潜能值是二氧化碳 的86倍， 政府间气候变化专门委员会（IPCC）估计， 迄今为止 ， 甲烷是观测到的全球变暖的30％的原因 。迅速减少向大气中泄漏的甲烷是应对气候变化的关键组成部分。2021年，联合国在气候变化会议（COP26） 上推出 了《 全球甲烷承诺 》，目标是 “对甲烷采取快速行动，使1.5摄氏度的未来触手可及”。该承诺有 150个签署国 ， 包括美国和欧盟。

尽早发现和持续监测甲烷源是对甲烷采取有意义行动的关键组成部分，因此已成为决策者和组织都关注的问题。大规模实施经济实惠、有效的甲烷检测解决方案，例如现场甲烷探测器或 机载光谱仪 ，具有挑战性，因为它们通常不切实际或昂贵得令人望而却步。另一方面，使用卫星进行遥感可以提供利益相关者所期望的全球规模、高频率和具有成本效益的探测功能。

在这篇博客文章中，我们将向您展示如何使用 托管在 亚马逊云科技 开放数据 注册表上的 S entinel 2 卫星图像 与 A mazon S ageMaker 地理空间 功能相结合，来检测甲烷排放的点源并随时间推移对其进行监控。根据 地球观测文献 的最新发现 ， 您将学习如何实现自定义的甲烷探测算法，并使用它来检测和监测来自全球不同地点的甲烷泄漏。这篇文章包含 GitHub 上的 附带代码 ，该代码提供了更多技术细节，可帮助您开始使用自己的甲烷监测解决方案。

传统上，进行复杂的地理空间分析是一项困难、耗时和资源密集型的工作。 Ama@@ zon SageMaker 地理空间 功能使数据科学家和机器学习工程师能够更轻松地使用地理空间数据构建、训练和部署模型。使用 SageMaker 地理空间功能，您可以高效地转换或丰富大规模地理空间数据集，使用预训练的机器学习 (ML) 模型加快模型构建，并使用三维加速图形和内置可视化工具在交互式地图上浏览模型预测和地理空间数据。

使用多光谱卫星图像对甲烷点源进行遥感

基于卫星的甲烷传感方法通常依赖于甲烷的独特透射率特性。在可见光谱中，CH4 的透射率值等于或接近 1，这意味着肉眼无法检测到。但是，在某些波长下，甲烷确实会吸收光线（透射率<1），这一特性可以用于探测目的。为此，通常选择短波长红外 (SWIR) 光谱（1500—2500 nm 光谱范围），这是最容易检测到甲烷的地方。超光谱和多光谱卫星任务（即使用光学仪器捕获电磁频谱多个波长范围（波段）内的图像数据的任务）涵盖这些短波红外范围，因此是潜在的探测仪器。 图 1 绘制了 SWIR 光谱中甲烷的透射率特性以及各种候选多光谱卫星仪器（改编自 本研究）的 SWIR 覆盖范围。

图 1 — SWIR 光谱和 Sentinel-2 多光谱任务覆盖范围中甲烷的透射率特性

许多多谱段卫星任务要么受到重访频率低（例如， PRISMA高光谱 约为16天）或空间分辨 率低（例如，7.5 km x 7.5 km 的 Sentin el 5 ）的限制。访问数据的成本是另一个挑战：一些专用星座作为商业任务运行，由于资金限制，可能会使研究人员、决策者和其他相关方不太容易获得甲烷排放见解。该解决方案以欧空局的 Sentinel-2多光谱任务 为基础，在重访率（大约 5 天）、空间分辨率（大约 20 m）和开放获取（托管在 亚马逊云科技 开放数据注册表 上）之间取得了适当的平衡。

Sentinel-2 有两个波段覆盖短波红外光谱 （分辨率为 20 m）：波段 11（中心波长 1610 nm）和波段 12（中心波长 2190 nm）。这两个波段都适用于甲烷检测，而波段12对甲烷吸收的灵敏度要高得多（见图 1）。直观地讲，使用这种 SWIR 反射率数据进行甲烷探测有两种可能的方法。首先，你可以只关注单个短波红外波段（最好是对甲烷吸收最敏感的波段），然后计算两个不同卫星通道的反射率逐像素差异。或者，您可以使用来自单个卫星通道的数据进行检测，方法是使用两个相邻的光谱 SWIR 波段，这两个波段具有相似的表面和气溶胶反射特性，但具有不同的甲烷吸收特性。

我们在本篇博文中实现的检测方法结合了这两种方法。我们借鉴了 地球观测文献 的最新发现 ，计算了两次卫星通道 和两个 SWIR 波段之间大气顶部 (TOA) 反射率 Δ（即 Sentinel-2 测量的反射率，包括来自大气气溶胶和气体的贡献）的分数变化；一次不存在甲烷的基线通道（基线）和怀疑存在活性甲烷点源的一次监测通道（监测）。从数学上讲，这可以表示如下：

方程 (1)

_{其中 ² 是 Sentinel-2 测量的 TOA 反射率，c ^监视器 和 c ^base 是通过将整个场景中波段 12 的 TOA 反射值与波段 11 的 TOA 反射值进行回归计算得出的（即 q b11 = c * β b12）。} 更多详情，请参阅这份关于通过 多光谱Sentinel-2卫星观测对异常甲烷点源进行 高频监测 的研究。

使用 SageMaker 地理空间功能实现甲烷检测算法

为了实现甲烷检测算法，我们在亚马逊 SageMaker Studio 中使用了 SageMaker 地理空间笔记本。地理空间笔记本内核预装了基本的地理空间库，例如 GDAL 、G eop andas、 Shapely 、 xar ray 和 Rasterio ，允许在 Python 笔记本环境中直接可视化和处理地理 空间数据。请参阅 入门指南，了解 如何开始 使用 SageMaker 地理空间功能。

SageMaker 提供了专门构建的 API ， 旨在使用 SearchRasterDataCollection API 调用，便于通过统一界面检索卫星图像。 searchRasterDataCollection 依赖于以下输入参数：

• Arn ：查询的栅格数据集合的亚马逊资源名称 (ARN)
• areaoFinterest ：代表搜索查询感兴趣区域的多边形对象（采用 GeoJSON 格式）
• timeRangeFilter ：定义感兴趣的时间范围，表示为 {startTime:，endTime:}
• Property Filters ：还可以纳入补充属性过滤器，例如最大可接受云量的规范

此方法支持从各种栅格数据源进行查询，可以通过调用 ListrasterDataCollecti ons 来浏览这些数据源。我们的甲烷探测实现使用 Sentinel-2 卫星图像，可以使用以下 ARN 进行全球引用：arn: aws: sagemaker-geospatial: us-w est-2:378778860802 ： raster-data-collection/public/nmqj48dcu3g7ayw8。

此 ARN 代表 Sentinel-2 图像，该图像已处理到 2A 级（表面反射率，经过大气校正）。出于甲烷探测的目的，我们将使用大气顶部 (TOA) 反射率数据（1C 级），其中不包括表面水平的大气校正，这些校正会使气溶胶成分和密度的变化（即甲烷泄漏）无法检测到。

为了识别来自特定点源的潜在排放，我们需要两个输入参数：可疑点源的坐标和甲烷排放监测的指定时间戳。鉴于 searchRasterDataCollection API 使用多边形或多面体来定义感兴趣区域 (AOI)，我们的方法包括先将点坐标扩展为边界框，然后使用该多边形使用 SearchRasterDateCollection 查询 Sentinel-2 图像。

在此示例中，我们监测了源自北非油田的已知甲烷泄漏事件。这是遥感文献中的标准验证案例，例如 在本 研究中被引用。 亚马逊 sagemaker 示例 GitHub 存储 库中提供了完全可执行的代码库。 在这里，我们仅重点介绍代表使用 SageMaker 地理空间功能实现甲烷检测解决方案的关键组成部分的部分代码。有关更多详细信息，请参见存储库。

我们首先初始化示例案例的坐标和目标监控日期。

#coordinates and date for North Africa oil field
#see here for reference: https://doi.org/10.5194/amt-14-2771-2021
point_longitude = 5.9053
point_latitude = 31.6585
target_date = '2019-11-20'
#size of bounding box in each direction around point
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],
    ])

#generate bounding box and extract polygon coordinates
aoi_geometry = bbox_around_point(point_longitude, point_latitude, distance_offset_meters)
aoi_polygon_coordinates = geometry.mapping(aoi_geometry)['coordinates']

#set search parameters
search_params = {
    "Arn": "arn:aws:sagemaker-geospatial:us-west-2:378778860802:raster-data-collection/public/nmqj48dcu3g7ayw8", # Sentinel-2 L2 data
    "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))
        }
    },
}
#query raster data using SageMaker geospatial capabilities
sentinel2_items = geospatial_client.search_raster_data_collection(**search_params)

响应包含匹配的 Sentinel-2 项目及其相应元数据的列表。其中包括 适用于所有 Sentinel -2 波段的云优化 GeoTiffs (COG) ，以及用于快速预览 图像 视觉波段 的缩略图 。当然，也可以访问全分辨率卫星图像（RGB 图），如下图 2 所示。

图 2 — AOI 的卫星图像（RGB 图）

如前所述，我们的探测方法依赖于大气顶部 (TOA) SWIR 反射率的部分变化。要使之奏效，确定良好的基准至关重要。寻找良好的基准可能很快就会成为一个繁琐的过程，需要进行大量的反复试验。但是，良好的启发式方法可以大大促进搜索过程的自动化。以下是对过去调查过的案例非常有效的搜索启发式方法：在过去的 day_offset=n 天内，检索所有卫星图像，移除所有云层，然后将图像剪辑到 AOI 范围内。然后计算 AOI 上的平均波段 12 反射率。返回波段 12 中平均反射率最高的图像的 Sentinel 图块 ID。

此逻辑在以下代码摘录中实现。其理由取决于波段12对甲烷吸收高度敏感的事实（见图 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):
    
    #get flattened arrays for regression
    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)
        
    #fit linear model using least squares regression
    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):
    
    #get correction factors
    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)
    
    #get corrected fractional reflectance change
    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):
    
    #get bounding box
    aoi_geometry = bbox_around_point(lon, lat, distance_offset)
    
    #get S2 metadata
    s2_meta_monitor = get_sentinel2_meta_data(date_monitor, aoi_geometry)
    
    #get tile 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"]
    
    #retrieve band 11 and 12 of the Sentinel L1C product for the given S2 tiles
    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)
    
    #compute corrected fractional reflectance change
    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

使用我们之前确定的参数运行此方法会得出 XArray.DataArray 形式的 SWIR TOA 反射率的部分变化。 我们可以通过在此数据数组上运行简单的 plo t () 调用来对结果进行首次目视检查。我们的方法揭示了AOI中心存在甲烷羽流，而在之前看到的RGB图中这是无法检测到的。

图 3 — TOA 反射率（SWIR 频谱）的分数反射率变化

作为最后一步，我们提取已识别的甲烷羽流，并将其叠加在原始的 RGB 卫星图像上，以提供重要的地理背景。这是通过阈值来实现的，其实现方式如下所示：

def get_plume_mask(change_in_reflectance_tif, threshold_value):
    cr_masked = change_in_reflectance_tif.copy()
    #set values above threshold to nan
    cr_masked[cr_masked > treshold_value] = np.nan 
    #apply mask on nan values
    plume_tif = np.ma.array(cr_masked, mask=cr_masked==np.nan) 
    
    return plume_tif

就我们的案例而言，阈值为 -0.02 的反射率小数变化会产生不错的结果，但这可能会因场景而异，因此您必须针对您的特定用例进行校准。下面的图 4 说明了如何通过将 AOI 的原始卫星图像与被掩盖的羽流合并为显示甲烷羽流在其地理环境中的单个合成图像来生成羽流叠加层。

图 4 — RGB 图像、TOA 反射率（SWIR 光谱）的分数反射率变化以及 AOI 的甲烷羽流叠加

使用真实世界的甲烷排放事件进行解决方案验证

作为最后一步，我们将评估我们的方法是否能够正确检测和查明来自不同来源和地区的甲烷泄漏情况。首先，我们使用专为 验证陆上甲烷排放的天基点源探测和量化而设计的甲 烷释放控制实验。 在2021年的这项实验中，研究人员在19天内在亚利桑那州的埃伦伯格进行了几次甲烷释放。在该实验期间，对其中一条 Sentinel-2 通道运行我们的探测方法，得出以下显示甲烷羽流的结果：

图 5 — 亚利桑那州控制释放实验的甲烷羽流强度

我们的探测方法可以清楚地识别控制释放期间产生的羽流。来自东亚垃圾填埋场（左）或北美石油和天然气设施（右）等来源的其他已知现实世界泄漏事件（如下图6所示）也是如此。

图 6 — 东亚垃圾填埋场（左）和北美油气田（右）的甲烷羽流强度

总而言之，我们的方法可以帮助识别来自控制释放和全球各种现实世界点源的甲烷排放。这最适合周围植被有限的陆上点源。它不适用于海上场景，因为水对短 波红外光谱的吸收率很高（即透过率低） 。 鉴于拟议的检测算法依赖于甲烷强度的变化，我们的方法还需要在泄漏前进行观察。这可能会使以恒定的排放速率监测泄漏变得困难。

清理

为避免甲烷监测任务完成后产生不必要的费用，请确保终止 SageMaker 实例并删除所有不需要的本地文件。

结论

通过将 SageMaker 地理空间功能与开放的地理空间数据源相结合，您可以大规模实施自己的高度定制的远程监控解决方案。这篇博文重点关注甲烷检测，这是寻求检测并最终避免有害甲烷排放的政府、非政府组织和其他组织的重点领域。通过使用 SageMaker 地理空间内核启动笔记本并实现自己的检测解决方案，您可以立即开始自己的地理空间分析之旅。查看 GitHub 存储库 ，开始构建自己的基于卫星的甲烷探测解决方案。另请查看 sagemaker 示例 库，以获取有关如何在其他现实世界遥感应用程序中使用 SageMaker 地理空间功能的更多示例和教程。

作者简介

卡斯滕·施罗尔博士 是 AW S 的解决方案架构师。他支持客户利用数据和技术来推动其IT基础设施的可持续性，并构建云原生数据驱动的解决方案，从而在各自的垂直领域实现可持续运营。Karsten 在应用机器学习和运营管理领域攻读博士学位后加入 亚马逊云科技。他对以技术为基础的社会挑战解决方案充满热情，并且喜欢深入研究构成这些解决方案基础的方法和应用架构。

Janosch Woschitz 是 亚马逊云科技 的高级解决方案架构师，专门研究地理空间 AI/ML。凭借超过 15 年的经验，他支持全球客户利用 AI 和 ML 开发利用地理空间数据的创新解决方案。他的专业知识涵盖机器学习、数据工程和可扩展的分布式系统，再加上强大的软件工程背景和自动驾驶等复杂领域的行业专业知识。