甲烷（CH4）是主要的人为温室气体，它是石油和天然气开采、煤矿、大规模养殖和废物处理等许多来源的副产品。甲烷的全球变暖潜力是二氧化碳的86倍（CH4的全球变暖潜力是二氧化碳的86倍），气候变化政府间专门委员会（IPCC）估计，从迄今观测到的全球变暖中，甲烷占30％的责任。迅速减少甲烷泄漏到大气中是对抗气候变化的关键组成部分。 2021年，联合国在气候变化大会（COP26）上推出了全球甲烷承诺，旨在“采取迅速行动控制甲烷，实现1.5摄氏度的未来”。该承诺有150个签署者，包括美国和欧盟。
早期发现和持续监测甲烷源是采取有意义的甲烷行动的关键组成部分，因此它正在成为决策者和组织关注的问题。实施规模经济、有效的甲烷检测解决方案（例如现场甲烷探测器或飞机搭载的光谱仪）具有挑战性，因为它们通常不切实际或过于昂贵。另一方面，使用卫星的遥感技术可以提供利益相关者所需的全球规模、高频率、经济高效的检测功能。
在本博客文章中，我们将向您展示如何使用存储在AWS开放数据注册表上的Sentinel 2卫星图像结合Amazon SageMaker地理空间能力来检测甲烷排放的点源并监测其变化。借鉴地球观测文献的最新发现，您将了解如何实施自定义甲烷检测算法并使用该算法来检测和监测全球各地各种地点的甲烷泄漏。本文附有GitHub上的附带代码，提供附加的技术细节，并帮助您开始使用自己的甲烷监测解决方案。
传统上，进行复杂的地理空间分析是一项困难、耗时和资源密集的工作。Amazon SageMaker地理空间能力使得数据科学家和机器学习工程师能够更轻松地构建、训练和部署使用地理空间数据的模型。使用SageMaker地理空间能力，您可以高效地转换或丰富大规模地理空间数据集，加速使用预训练的机器学习（ML）模型构建模型，并使用3D加速图形和内置可视化工具在交互式地图上探索模型预测和地理空间数据。
使用多光谱卫星图像远程监测甲烷点源
基于卫星的甲烷感测方法通常依赖于甲烷的独特透过特性。在可见光谱中，甲烷的透过值等于或接近于1，意味着肉眼无法检测到。然而，在某些波长范围内，甲烷确实吸收光（透过率<1），这一特性可用于检测目的。通常选择短波红外（SWIR）光谱（1500-2500 nm光谱范围），在该范围内甲烷最容易被检测到。超光谱和多光谱卫星任务（即具有在电磁谱范围内由多个波长范围（波段）捕获图像数据的光学仪器）覆盖这些SWIR范围，因此具有潜在的检测功能。图1绘制了甲烷在SWIR光谱中的透过特性以及各种候选多光谱卫星仪器的SWIR覆盖范围（改编自this研究）。
图1 – 甲烷在SWIR光谱中的透过特性和Sentinel-2多光谱任务的覆盖范围
许多多光谱卫星任务要么受限于低重访频率（例如PRISMA Hyperspectral，约为16天），要么受限于低空间分辨率（例如Sentinel 5，为7.5 km x 7.5 km）。获取数据的成本是另一个挑战：一些专用星座作为商业任务运营，可能会由于财务限制而使甲烷排放的洞察力对研究人员、决策者和其他相关方不太容易获得。本解决方案基于ESA的Sentinel-2多光谱任务，平衡了重访频率（约为5天）、空间分辨率（约为20米）和开放访问（托管在AWS开放数据注册表上）的合适性。
Sentinel-2有两个覆盖SWIR光谱的波段（分辨率为20米）：波段-11（1610 nm中心波长）和波段-12（2190 nm中心波长）。这两个波段都适用于甲烷检测，而波段-12对CH4吸收具有更高的敏感性（参见图1）。直观上，对于使用此SWIR反射数据进行甲烷检测，有两种可能的方法。首先，可以专注于单个SWIR波段（理想情况下是对CH4吸收最敏感的波段），并计算两个不同卫星通道之间反射率的像素级差异。或者，可以使用一次单一卫星通道的数据进行检测，该卫星通道具有相似的地表和气溶胶反射特性，但具有不同的甲烷吸收特性的两个相邻光谱SWIR波段。
我们在这篇博文中实施的检测方法结合了这两种方法。我们利用最近从地球观测文献中获得的结果，并在两次卫星通道和两个SWIR波段之间计算大气顶部反射率Δρ（即Sentinel-2测量的反射率，包括大气气溶胶和气体的贡献）的分数变化，其中一次是不存在甲烷的基准通道（base），一次是怀疑存在活跃甲烷点源的监控通道（monitor）。数学上可表示为：
方程式（1）
其中ρ是由Sentinel-2测得的大气顶部反射率，cmonitor和cbase是通过回归波段-12的大气顶部反射率值与波段-11的反射率值得到的整个场景上的值（即ρ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表示经过处理为Level 2A（表面反射率，大气校正）的Sentinel-2图像。为了进行甲烷检测，我们将使用大气层顶（TOA）反射率数据（Level 1C），它不包括会使气溶胶组成和密度变化（即甲烷泄漏）无法检测到的地表大气层校正。
为了识别特定点源的潜在排放物，我们需要两个输入参数：疑似点源的坐标和指定的甲烷排放监测时间戳。由于SearchRasterDataCollection API使用多边形或多多边形来定义感兴趣区域（AOI），我们的方法涉及首先将点坐标扩展为边界框，然后使用该多边形使用SearchRasterDateCollection查询Sentinel-2图像。
在此示例中，我们监测来自非洲北部油田的已知甲烷泄漏。这是遥感文献中的标准验证案例，例如在此研究中引用。可在amazon-sagemaker-examples GitHub存储库上提供完整的可执行代码基础。在这里，我们仅强调代表使用SageMaker地理空间功能实施甲烷检测解决方案的关键构建块的选定代码部分。有关更多详细信息，请参阅存储库。
我们首先通过初始化示例案例的坐标和目标监测日期来开始。
#非洲北部油田的坐标和日期#参考：https://doi.org/10.5194/amt-14-2771-2021point_longitude = 5.9053point_latitude = 31.6585target_date = '2019-11-20'#每个方向上的边界框大小，以米为单位distance_offset_meters = 1500以下代码片段为给定的点坐标生成一个边界框，然后根据边界框和指定的监测日期搜索可用的Sentinel-2图像：
def bbox_around_point(lon, lat, distance_offset_meters): #从https://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html获取的赤道半径（km） 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项目及其对应的元数据列表。这些元数据包括所有优化云地理TIFF(COG)的Sentinel-2波段，以及缩略图图像，用于快速预览图像的可见光波段。当然，也可以访问完整分辨率的卫星图像（RGB绘图），如下所示的图2。
图2 – AOI的卫星图像（RGB绘图）
如前所述，我们的检测方法依赖于大气顶层（TOA）SWIR反射中的分数变化。为了使其有效，寻找一个好的基准线至关重要。找到一个好的基准线可能很快变得繁琐，需要大量的尝试和错误。然而，良好的启发式算法可以在自动化搜索过程中发挥重要作用。过去的一个day_offset=n天内，检索所有卫星图像，去除任何云层并裁剪到感兴趣区域（AOI）。然后计算AOI内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] #重塑形状 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"] #检索给定S2瓦片的Sentinel 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的原始卫星图像与掩蔽喷漏合成为单一复合图像来生成喷漏叠加效果。
图4- RGB图像，TOA反射率的分数变化（SWIR频谱）和AOI的甲烷喷漏叠加
使用真实世界甲烷排放事件进行方案验证
作为最后一步，我们评估我们的方法在从各种来源和地理区域正确检测和确定甲烷泄漏的能力。首先，我们使用一个专门设计用于验证陆地上甲烷点源检测和定量的受控甲烷释放实验。在这个2021年的实验中，研究人员在亚利桑那州埃伦伯格进行了一个19天的甲烷释放实验。在该实验期间的一个Sentinel-2过程中运行我们的检测方法，会产生以下结果，显示了甲烷喷漏：
图5- 亚利桑那州受控释放实验的甲烷喷漏强度
我们的检测方法明确识别出了受控释放期间产生的喷漏。对于其他已知的真实世界泄漏情况（见下图6），如东亚的垃圾填埋场（左）或北美的油气设施（右），也是如此。
图6- 东亚垃圾填埋场（左）和北美油气田（右）的甲烷喷漏强度
总之，我们的方法可以帮助确定全球各地受控释放和各种真实世界点源的甲烷排放。这对于陆地上周围植被有限的点源效果最佳。对于离岸场景由于水的SWIR谱吸收较高（即透射较低），该方法无效。鉴于我们的检测算法依赖于甲烷强度的变化，我们的方法也需要事前泄漏观察结果。这使得对于具有恒定排放率的泄漏的监测变得具有挑战性。
清理
在完成甲烷监测作业后为避免产生不必要的费用，请确保终止SageMaker实例并删除任何不需要的本地文件。
结论
通过将SageMaker的地理空间能力与开放地理空间数据源相结合，您可以在规模上实施自己高度定制的远程监测解决方案。本博文侧重于甲烷检测，这是政府、非政府组织和其他寻求检测和最终避免有害甲烷排放的机构的焦点领域。您可以立即开始自己进入地理空间分析之旅，通过启动带有SageMaker地理空间内核的笔记本并实现自己的检测方案。请查看GitHub存储库，开始构建基于卫星的甲烷检测解决方案。还请查看sagemaker-examples存储库，了解更多关于如何在其他真实远程感知应用中使用SageMaker地理空间能力的示例和教程。