Использование спутниковых данных высокого пространственного разрешения Sentinel-2 для выявления вырубок

В ИКИ РАН за отчётный период разработана технология выявления вырубок на основе спутниковых данных высокого пространственного разрешения Sentinel-2, полученных в условиях наличия снежного покрова на земной поверхности. Приборы серии Sentinel-2 позволяют регулярно получать изображения с периодичностью каждые 2–3 дня в средних широтах.

Каждый спутник Sentinel-2 оснащен мультиспектральным прибором MSI с рядом каналов в видимом/ближнем инфракрасном диапазоне (VNIR) и коротковолновом инфракрасном спектральном диапазоне (SWIR) (таблица 1).

Таблица 1 — Основные спектральные каналы прибора Sentinel-2

Канал	Длина волны, нм	Пространственное разрешение, м
Band 2 – Blue	492,4	10
Band 3 – Green	559,8	10
Band 4 – Red	664,6	10
Band 8 – NIR	832,8	10
Band 11 – SWIR	1613,7	20
Band 12 – SWIR	2202,4	20

 

Архив данных пополняется с 2019 г. продуктами уровня Level-2A с учтённым влиянием атмосферы, характеризующими отражательную способность земной поверхности. Эти данные могут быть получены через портал Copernicus Data Space Ecosystem (https://dataspace.copernicus.eu).

Специфика данных высокого разрешения позволяет выявлять чёткие границы большинства вырубок, которые при более низком пространственном разрешении размываются и становятся менее заметными. Использование композитных изображений снежного покрова снижает вероятность обнаружения, поскольку смешивание в случайных пропорциях яркостных характеристик здорового леса и открытой снежной поверхности, соответствующей вырубке, приводит к значительной неопределённости яркостных характеристик вырубки. Таким образом, для выявления вырубок используется временной ряд наблюдений, соответствующий снежному периоду и предварительно очищенный от влияния факторов облачности и сопутствующих теней. Логическая схема метода детектирования рубок приведена на рисунке 1.

 Рисунок 1 — Логическая схема метода детектирования рубок

Этапы детектирования вырубок включают предварительную обработку спутниковых данных, направленную на очищение набора исходных данных от влияния облачности, теней, шумов, а также на выделение территорий, не покрытых снегом; анализ двухлетнего временного ряда динамики значений набора спектральных характеристик; определение даты детектирования изменения состояния лесного покрова; а также пространственную коррекцию границ предполагаемых вырубок на основе анализа межгодовой динамики исследуемых значений.
Детектирование вырубок и зимних изменений основано на анализе временных рядов зимних наблюдений. В качестве исходных данных используются наблюдения, полученные в каналах RED, BLUE, NIR, SWIR, информативных при анализе как снежного, так и растительного покрова. Изменения значений в этих каналах, а также в индексах, построенных на их основе, позволяют предположить изменения класса и состояния растительного покрова и определить временной интервал изменений.

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

Данные проходят предварительную обработку, направленную в первую очередь на разделение облачного и снежного покрова, а также фильтрацию поверхности, не прокрытой снегом. Особенно важно исключить возможность перепутывания с рубками изменений в лесах, связанных с началом ранней вегетации. Поэтому, помимо индекса снега NDSI, при анализе временных рядов используется вегетационный индекс NDVI. При появлении рубки значение NDVI уменьшается, поскольку вырубка представляет собой более яркий объект по сравнению с лесом, и при анализе динамики индекса наблюдается его снижение в момент изменения состояния лесного покрова.
Для фильтрации облаков и тёмных объектов используется набор эмпирически полученных пороговых значений:

R_SWIR > 0,012, R_SWIR < 0,2, R_BLUE > 0,02.

Кроме того, на основе исходных данных и пространственного гистограммного анализа формируются маски бесснежных территорий, а дополнительная фильтрация территории, не покрытой снегом, осуществляется с использованием следующих условий:

R_RED + R_NIR + R_BLUE + R_SWIR > 0,15, NDVI > –0,5648NDSI + 0,092, NDSI_R + NDSI_B < 0.

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

Для фильтрации облачного покрова, который часто смешивается со снежным покровом, и может существенно исказить данные, после применения порогового подхода используется анализ двухлетнего временного ряда, за текущий и предшествующий годы. Величина среднего значения и стандартного отклонения в двухлетнем временном ряду позволяют выполнять статистическую фильтрацию выбросов в каналах RED, BLUE, NIR, SWIR. Пороги по числу допустимых стандартных отклонений подбираются адаптивно таким образом, чтобы суммарно не исключать более 20 % наблюдений по всем исследуемым характеристикам и не удалять из рассмотрения изменения лесного покрова, произошедшие в конце периода наблюдений.

Для временного ряда выполняется оценка межгодовой изменчивости статистических характеристик. Для этого для текущего и предыдущего года определяются среднее значение MEAN и величина стандартного отклонения DISP. Появление рубки предполагает повышение как средней яркости наблюдений, так и, соответственно, возрастание дисперсии. Сравнение этих величин позволяет дать предварительную оценку вероятности изменения состояния растительного покрова на основе среднего значения NDVI и яркости в NIR канале и стандартного отклонения в NIR канале:

(meanNIR_c – meanNIR_r) > 0, (dispNIR_c – DispNIR_r) > 0,
                             (meanNDVI_c – meanNDVI_r) < 0.                                (1)

Ближний инфракрасный канал отличается высоким уровнем сигнала по сравнению с красным и голубым каналами, при этом в нем отсутствует эффект поглощения для зеленых объектов, таких как хвойные леса. Изменения в NIR канале затем используются в качестве основы для пространственной коррекции границ вырубок.

Вегетационный индекс NDVI используется для фильтрации изменений, связанных с началом вегетационного сезона. При появлении рубки среднее значение NDVI не должно увеличиваться.

При выполнении условий (1) определяется число наблюдений, в которых значение яркости R_c превышает среднее значение MEAN_r за предыдущий год не менее чем на четыре стандартных отклонения:

R_c > MEAN_r + 4DISP_r).

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

Если суммарное число изменений Ch, т.е. количество дат, для которых в каналах RED, NIR и BLUE выполняется условие

Ch = Ch_RED + Ch_NIR + Ch_BLUE > 1,

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

Данные Sentinel-2 имеют нерегулярное временное разрешение. Кроме того, значительный объём данных отфильтровывается из-за облачного покрова и других мешающих факторов, в связи с чем точное определение даты осуществления рубки маловероятно. Алгоритм определяет дату детектирования рубки и предполагаемый временной диапазон возникновения изменения.

Для определения даты изменения состояния растительного покрова выполняется расчёт среднего значения и стандартного отклонения дисперсии до момента t — DispRED_R (с учётом данных за предыдущий год) — и после момента t — DispRED_L — для каждого наблюдения текущего года. Суммарная дисперсия до и после рассматриваемой даты позволяет оценить относительный перепад значений, например, в красном канале:

pRED = (RED_t – RED_t–1)/(DispRED_R + DispRED_L).

В момент времени t фиксируется изменение, если суммарный перепад значений в 3 каналах удовлетворяет условию

pRED + pNIR + pBLUE > 1,

а относительный перепад значений pR в одном или нескольких каналах является максимальным в рассматриваемом временном ряду. Кроме того, необходимым условием определения временного диапазон возникновения изменения [t-1; t] является увеличение дисперсии после даты t во всех трех каналах:

DispBLUER ≤ DispBLUEL, DispNIRR ≤ DispNIRL, DispREDR≤ DispREDL,

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

Даты изменения для каналов NIR, RED, BLUE подбираются независимо, после чего путем объединения полученных временных интервалов определяется предполагаемый диапазон возникновения вырубки.

На рисунке 2 представлен пример анализа двухлетнего набора наблюдений. Оранжевыми стрелками обозначены точки, отфильтрованные в ходе статистической фильтрации. В приведённом примере в обоих случаях фильтрация обусловлена скачком значений в SWIR-канале, который позволяет разделять облачный и снежный покров. Синими стрелками указана дата максимальной разделимости группы точек «до» и «после» по приведённому набору признаков.

Несмотря на многоуровневую фильтрацию, реальные данные характеризуются высоким разбросом значений, а значения КСЯ в каналах имеют тенденцию к увеличению к концу снежного сезона, что объясняется не только вероятным таянием снега, но и изменением положения Солнца. По этой причине набор условий для определения вырубок позволяет детектировать не всю площадь вырубки, и требуется дополнительный пространственный анализ для уточнения её границ. Для этого разработан пространственный фильтр, позволяющий присоединять точки в случае, если перепад среднего значения яркости в NIR-канале и величины стандартного отклонения вне границ вырубки (out) превышает соответствующие значения внутри области изменения (in) а отклонения межгодовых статистических характеристик находится в диапазоне 20 %:

meanNIR_in – meanNIR_out > –0,2meanNIR_in, dispNIR_in – dispNIR_out > –0,2dispNIR_in.

Рисунок 2 — Анализ двухлетнего временного ряда значений в каналах RED, BLUE, NIR, SWIR и индекса NDVI

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

На рисунке 3 приведены два примера применения пространственного фильтра для расширения территории лестной вырубки на основе анализа межгодовой динамики статистических величин КСЯ в NIR-канале.

      
а                                                        б                                                           в
                 
г                                                 д                                                      е
Рисунок 3 — Уточнение границ вырубок на основе анализа межгодовой динамики КСЯ: а, г — вырубки до применения пространственного фильтра; б, д — рубки после применения пространственного фильтра; в, е — RGB-синтез (2023–2022–2022) NIR-канала при 30°

На рисунках 4–7 показан пример работы алгоритма детектирования рубок в районе Ангары, где на момент наблюдения происходят активные вырубки леса. На рисунках 6 и 27 представлен увеличенный фрагмент, выделенный на рисунках 4 и 5 красной рамкой. На рисунках 5 и 7 результат детектирования вырубок приведен в датах детектирования изменения, в числовом диапазоне от 35 до 150 дней года.

      
 

Рисунок 4 — RGB-синтез (2023–2022–2022) NIR-канала при 30°; координаты центра изображения: 58°29ʹ02,86ʺ с.ш., 101°31ʹ06,93ʺ в.д.	Рисунок 5 — Пространственное распределение дат вырубок и зимних изменений 2023 г.

    

Рисунок 6 — RGB-синтез (2023–2022–2022) NIR канала при 30°; координаты центра изображения: 58°22ʹ13,03ʺ с.ш., 102°02ʹ19,01ʺ в.д.	Рисунок 7 — Пространственное распределение дат вырубок и зимних изменений 2023 г.

 

Границы временного диапазона возникновения рубки записываются как даты DATA_1 и DATA_2. Гистограммы распределения значений для тестовой территории (рисунок 4) приведены на рисунке 8. Распределение значений DATA_1 во временном диапазоне располагается левее распределение значений DATA_2. Гистограмма распределения временной неопределённости даты возникновения вырубки приведена на рисунке 9. Формы и пики гистограмм отражают наличие или отсутствие данных, что ограничивает точность определения даты вырубки. В частности, 50 % вырубок определяются с точностью до 5 дней, 75 % — до 17 дней, а 90 % — до 25 дней.

На изображениях видно, что метод позволяет детектировать как вырубки, хорошо различимые по композитным изображениям, так и едва заметные, возникающие поздно и оказывающие слабое влияние при построении композитов. Особенно это проявляется на вырубках за 2024 г.

      

Рисунок 8 — Гистограмма распределения границ временного интервала: до и после вырубки	Рисунок 9 — Гистограмма распределения временной неопределённости при детектировании вырубок

На рисунках 10–13 приведены примеры детектирования рубок за разные годы. На левых изображениях представлены зимние композитные изображения в RGB синтезе (2022–2023–2024) для NIR-канала при 30°. На правых изображениях можно видеть наложенные сверху рубки, детектированные за 2022 г. (красный цвет), 2023 г. (зелёный), 2024 г. (синий).

 

Рисунок 10 — RGB-синтез (2022–2023–2024) NIR канала при 30°; координаты центра изображения: 56°08ʹ10,19ʺ с.ш., 99°45ʹ59,22ʺ в.д.	Рисунок 11 — Рубки за 2022 г. (красный цвет), 2023 г. (зелёный), 2024 г. (синий)

Рисунок 12 — RGB-синтез (2022–2023–2024) NIR канала при 30°; координаты центра изображения: 56°27ʹ15,75ʺ с.ш., 99°13ʹ26,64ʺ в.д.	Рисунок 13 — Рубки за 2022 г. (красный цвет), 2023 г. (зелёный), 2024 г. (синий)