Обнаружение и мониторинг источников выбросов метана с помощью геопространственных возможностей Amazon SageMaker высокой частоты.

Выявление и контроль источников выделения метана с использованием высокочастотных геопространственных возможностей Amazon SageMaker.

Метан (CH4) – основной антропогенный парниковый газ, являющийся побочным продуктом нефте- и газодобычи, угледобывания, промышленного животноводства и утилизации отходов, среди прочих источников. Потенциал глобального потепления для CH4 составляет 86 раз больше, чем у CO2, и Межправительственная панель по изменению климата (МПИК) оценивает, что метан отвечает за 30 процентов наблюдаемого глобального потепления на сегодняшний день. Быстрое снижение выбросов CH4 в атмосферу является критическим аспектом в борьбе с изменением климата. В 2021 году ООН представила Глобальное обязательство по метану на Конференции по изменению климата (COP26) с целью принять “быстрые меры по метану для достижения будущего с температурой не выше 1,5 °C”. Обязательство подписали 150 стран, включая США и ЕС.

Раннее обнаружение и непрерывное мониторинг источников метана – ключевые аспекты значимых мер по метану и поэтому вызывают беспокойство как у политиков, так и у организаций. Как реализовать доступные и эффективные решения по обнаружению метана в масштабе – например, датчики метана на месте или установленные на самолеты спектрометры – представляют собой сложность, поскольку они часто неудобны или имеют высокую стоимость. С другой стороны, дистанционное зондирование, использующее спутники, может предоставить желаемые функциональные возможности по обнаружению на мировом уровне, с высокой частотой и стоимостью.

В этом блоге мы покажем вам, как вы можете использовать Спутниковые снимки Sentinel 2, хранящиеся в AWS Реестре открытых данных в сочетании с геопространственными возможностями Amazon SageMaker для обнаружения точечных источников выброса метана CH4 и их мониторинга со временем. Исходя из недавних исследований в сфере наблюдения Земли, вы узнаете, как реализовать собственный алгоритм обнаружения метана и использовать его для поиска и мониторинга утечек метана в различных местах по всему миру. В посте приведен сопроводительный код на GitHub, который предоставляет дополнительные технические подробности и помогает вам начать работу с вашим собственным решением по мониторингу метана.

Традиционно, выполнение сложных геопространственных анализов было сложной, трудоемкой и ресурсоемкой задачей. Геопространственные возможности Amazon SageMaker позволяют упростить процесс создания, обучения и развертывания моделей с использованием геопространственных данных для ученых-исследователей данных и инженеров машинного обучения. Используя возможности геопространственного модуля SageMaker, вы можете эффективно трансформировать или обогащать геопространственные наборы данных большого масштаба, ускорять построение моделей с помощью предварительно обученных моделей машинного обучения (ML) и исследовать прогнозы моделей и геопространственные данные на интерактивной карте с использованием графики, ускоренной 3D, и инструментов визуализации.

Дистанционное зондирование точечных источников метана с использованием мультиспектральных спутниковых снимков

Спутниковые методы обнаружения метана обычно основаны на уникальных характеристиках пропускания CH4. В видимом спектре метан имеет значение пропускания, равное или близкое к 1, что означает его невозможность обнаружения невооруженным глазом. Однако на определенных длинах волн метан поглощает свет (пропускание <1), что можно использовать для обнаружения. Для этого обычно выбираются длинные волны инфракрасного спектра (SWIR) (спектральный диапазон 1500-2500 нм), где метан наиболее обнаружим. Мультиспектральные спутниковые миссии, осуществляемые в гипер- и мультиспектральных диапазонах, охватывают эти диапазоны SWIR и, следовательно, представляют потенциальные инструменты обнаружения. Рисунок 1 показывает характеристики пропускания метана в спектре SWIR и покрытие SWIR различных кандидатов мультиспектральных спутниковых приборов (адаптировано из этого исследования).

Рисунок 1 - Трансмитантные характеристики метана в спектральном диапазоне SWIR и охват мультиспектральных миссий Sentinel-2

Рисунок 1 – Трансмитантные характеристики метана в спектральном диапазоне SWIR и охват мультиспектральных миссий Sentinel-2

Многие мультиспектральные спутниковые миссии ограничены либо низкой частотой повторных пролетов (например, гиперспектральная миссия PRISMA примерно раз в 16 дней), либо низким пространственным разрешением (например, Sentinel 5 с разрешением 7,5 км x 7,5 км). Дополнительная проблема – стоимость доступа к данным: некоторые специализированные констелляции работают как коммерческие миссии, что может затруднить доступ к данным о выбросах CH4 для исследователей, принимающих решения и других заинтересованных сторон из-за финансовых ограничений. Мультиспектральная миссия Sentinel-2 ЕКА, на основе которой разработано это решение, обеспечивает баланс между частотой повторных пролетов (примерно раз в 5 дней), пространственным разрешением (примерно 20 м) и открытым доступом (размещено в AWS Registry of Open Data).

Sentinel-2 имеет две полосы, охватывающие спектр SWIR (с разрешением 20 м): полосы 11 (центральная длина волны 1610 нм) и 12 (центральная длина волны 2190 нм). Обе полосы подходят для обнаружения метана, при этом полоса 12 обладает значительно большей чувствительностью к поглощению CH4 (см. рисунок 1). Интуитивно можно выделить два возможных подхода к использованию этих данных о светоотражении в спектре SWIR для обнаружения метана. Во-первых, можно сосредоточиться только на одной полосе SWIR (желательно на той, которая наиболее чувствительна к поглощению CH4) и вычислить разницу в светоотражении пикселя к пикселю между двумя различными спутниковыми пролетами. Второй подход – использовать данные из одного спутникового пролёта для обнаружения, используя две соседние спектральные полосы SWIR, которые имеют схожие свойства отражения поверхности и аэрозольного отражения, но различные характеристики поглощения метана.

Метод обнаружения, реализованный в этой статье, объединяет оба подхода. Мы опираемся на недавние исследования в области наблюдений Земли и вычисляем долевое изменение отражения верхней границы атмосферы (TOA) Δρ (то есть отражение, измеренное Sentinel-2 и включающее в себя вклад атмосферных аэрозолей и газов) между двумя спутниковыми пролётами и двумя SWIR полосами. Одна полоса базового пролёта, где метан отсутствует (база), и один полоса мониторингового пролета, где предполагается наличие активного источника метана (монитор). Математически это можно выразить следующим образом:

Уравнение 1Уравнение (1)

где ρ – отражение TOA, измеренное Sentinel-2, cmonitor и cbase вычисляются путём регрессии значений отражения TOA полосы 12 относительно значений полосы 11 по всей сцене (то есть ρb11 = c * ρb12). Дополнительные детали см. в исследовании по высокочастотному мониторингу аномальных источников метана с использованием мультиспектральных наблюдений спутника Sentinel-2.

Реализация алгоритма обнаружения метана с использованием геопространственных возможностей SageMaker

Для реализации алгоритма обнаружения метана мы используем геопространственную записную книжку SageMaker в рамках Amazon SageMaker Studio. Ядро геопространственной записной книжки предварительно оснащено необходимыми геопространственными библиотеками, такими как GDAL, GeoPandas, Shapely, xarray и Rasterio, что обеспечивает прямую визуализацию и обработку геопространственных данных в среде Python записной книжки. Инструкции по началу работы можно найти в руководстве по началу работы по использованию геопространственных возможностей SageMaker.

SageMaker предоставляет специализированный API, разработанный для упрощения получения спутниковых изображений через объединенный интерфейс с использованием вызова API SearchRasterDataCollection. SearchRasterDataCollection использует следующие входные параметры:

  • Arn: Амазонское имя ресурса (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), которые не включают коррекции атмосферы на уровне поверхности, которые делают изменения в аэрозольном составе и плотности (то есть утечки метана) незаметными.

Для определения потенциальных выбросов из конкретного источника мы нуждаемся в двух входных параметрах: координатах предполагаемого источника и отметке времени для мониторинга выбросов метана. Учитывая, что API SearchRasterDataCollection использует многоугольники или мульти-полигоны для определения интересующей области (AOI), наш подход заключается в расширении координат точки в прямоугольник, а затем использовании этого полигона для запроса изображений Sentinel-2 с помощью SearchRasterDateCollection.

В этом примере мы мониторим известную утечку метана, исходящую из нефтяного месторождения на Севере Африки. Этот пример является стандартным тестом из литературы по дистанционному зондированию и упоминается, например, в этом исследовании. Полностью выполняемый код предоставляется в репозитории GitHub amazon-sagemaker-examples. Здесь мы представляем только выбранные секции кода, которые являются основными строительными блоками для реализации решения по обнаружению метана с использованием возможностей 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    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 уровня 2    "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 и их соответствующих метаданных. В нем включены Cloud-Optimized GeoTIFFs (COG) для всех полос Sentinel-2 , а также изображения эскизов для быстрого предварительного просмотра изображений визуальных полос. Естественно, также возможен доступ к спутниковому изображению полного разрешения (построение RGB), показанному на Рисунке 2, который следует.

Рисунок 2Рисунок 2 – Спутниковое изображение (построение RGB) AOI

Как уже было подробно описано ранее, наш метод детекции основывается на дробных изменениях отражательной способности в спектральной области кратковолнового инфракрасного излучения над поверхностью атмосферы. Для того чтобы это работало, нахождение хороших базовых изображений является важным. Поиск хороших базовых изображений может стать утомительным процессом, требующим множество экспериментов. Однако, хорошая эвристика может сэкономить много времени при автоматизации этого процесса. Одной из эвристик, которая хорошо сработала в прошлых исследованиях, является следующая: в течение последних day_offset=n дней извлекайте все спутниковые изображения, удаляйте облачность и обрезайте изображение по границе AOI. Затем вычисляйте среднюю отражательную способность для полосы 12 в области AOI. Возвращайте идентификатор Sentinel-снимка с самой высокой средней отражательной способностью для полосы 12.

Данная логика реализована в следующем фрагменте кода. Ее основа основывается на том, что полоса 12 чувствительна к поглощению CH4 (см. Рисунок 1). Большее значение средней отражательной способности соответствует меньшему поглощению из источников, таких как выбросы метана, и поэтому предоставляет надежный показатель для сцены без выбросов.

def approximate_best_reference_date(lon, lat, date_to_monitor, distance_offset=1500, cloud_mask=True, day_offset=30):        #инициализируем AOI и другие параметры    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               #перебираем n=day_offset предыдущих дней    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')                   #получаем тайлы Sentinel-2 для текущей даты        s2_tiles_for_target_date = get_sentinel2_meta_data(target_date, aoi_geometry)                #перебираем доступные тайлы для текущей даты        for s2_tile_meta in s2_tiles_for_target_date:            s2_tile_id_to_test = s2_tile_meta['Id']            #извлекаем проверенные на облачность (при необходимости) данные по полосе L1C 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)            #вычисляем среднюю отражательную способность для полосы SWIR            mean_swir = target_band_data.sum() / target_band_data.count()                        #обеспечиваем, чтобы площадь области, видимой/не покрытой облаками, была достаточно большой            visible_area_ratio = target_band_data.count() / (target_band_data.shape[1] * target_band_data.shape[2])            if visible_area_ratio <= 0.7: #<-- обеспечиваем приемлемое покрытие облаками                continue                        #обновляем максимальный ref_s2_tile_id и ref_target_date при необходимости            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. Идентификаторы тайлов Sentinel-2 содержат информацию о миссии (Sentinel-2A/Sentinel-2B), уникальный номер тайла (например, 32SKA) и дату создания изображения, а также другую информацию, и уникально идентифицируют наблюдение (то есть сцену). В нашем примере процесс приближения указывает на 6 октября 2019 года (тайл 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 (Area of Interest) для заданной долготы и широты, даты мониторинга и базового тайла, получает необходимые спутниковые изображения и выполняет вычисление фракционного изменения отражаемости.

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)        #получить идентификатор тайла    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"]        #извлечь полосу 11 и 12 продукта Sentinel L1C для заданных тайлов S2    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

Запуск этого метода с ранее определенными параметрами дает фракционное изменение отражаемости TOA SWIR в виде xarray.DataArray. Мы можем провести первое визуальное исследование результатов, запустив простую команду plot() для данного массива данных. Наш метод показывает наличие метанового плумма в центре AOI, которая была неразличима на ранее видимом плане RGB.

Рисунок 3Рисунок 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

Рисунок 4 – RGB-изображение, доля изменения отражения в отражении TOA (спектр SWIR) и наложение метанового выброса для AOI

Проверка решения на реальных событиях выброса метана

В заключение, мы оцениваем наш метод на его способность правильно обнаруживать и уточнять утечки метана из различных источников и географических регионов. Сначала мы используем контролируемый эксперимент по выбросу метана, специально разработанный для валидации детектирования и количественной оценки выбросов точечного источника метана на суше на основе данных из космоса. В этом эксперименте 2021 года исследователи провели несколько выбросов метана в Эренберге, штат Аризона, в течение 19 дней. Запуск нашего метода детектирования для одного из проходов Sentinel-2 во время этого эксперимента дает следующий результат, показывающий метановый выброс:

Рисунок 5Рисунок 5 – Интенсивность метанового выброса для Контролируемого эксперимента по выбросу метана в Аризоне

Выброс, сгенерированный во время контролируемого выброса, четко определяется нашим методом детектирования. То же самое верно и для других известных реальных выбросов (на рисунке 6), таких как свалка в Восточной Азии (слева) или нефтегазовый объект в Северной Америке (справа).

Рисунок 6Рисунок 6 – Интенсивность метанового выброса для свалки в Восточной Азии (слева) и объекта нефтегазовой промышленности в Северной Америке (справа)

В итоге, наш метод может помочь определить выбросы метана как от контролируемого выброса, так и от различных реальных точечных источников по всему миру. Он работает лучше для точечных источников на суше с ограниченной окружающей растительностью. Он не работает для морских сцен из-за высокой поглощения (то есть низкой пропускной способности) спектра SWIR водой. Учитывая, что предлагаемый алгоритм детектирования полагается на вариации интенсивности метана, наш метод также требует предварительных наблюдений до выброса. Это может усложнить контроль за выбросами с постоянными уровнями выбросов.

Очистка

Чтобы избежать нежелательных расходов после завершения задачи по мониторингу метана, убедитесь, что вы завершили экземпляр SageMaker и удалите ненужные локальные файлы.

Заключение

Путем объединения геопространственных возможностей SageMaker с открытыми источниками геопространственных данных вы можете реализовать собственные сильно настраиваемые решения для удаленного мониторинга в масштабе. В этой статье акцент был сделан на обнаружении метана, которое является фокусной областью для правительств, НПО и других организаций, стремящихся обнаружить и в конечном итоге избежать вредных выбросов метана. Вы можете начать сегодня свой собственный путь в геопространственном анализе, создав Ноутбук с геопространственным ядром SageMaker и реализовав свое собственное решение обнаружения. Посетите репозиторий GitHub, чтобы начать создание собственного решения обнаружения метана на основе спутниковых данных. Также ознакомьтесь с репозиторием sagemaker-examples для дополнительных примеров и учебных пособий о том, как использовать геопространственные возможности SageMaker в других реальных приложениях дистанционного зондирования.