Исследование возможности оценки возобновления лесной растительности после сплошных рубок по спутниковым данным LANDSAT (на примере Брянского полесья)

© 2019 г.

Е.И. Белова^*, Д.В. Ершов

Центр по проблемам экологии и продуктивности лесов РАН

Россия, 117997, г. Москва, ул. Профсоюзная, 84/32, стр. 14.

^*E-mail: belova@ifi.rssi.ru

Поступила в редакцию 23.09.2019

В статье рассмотрен метод дистанционной оценки восстановления лесных культур на территории Суземского лесничества Брянской области. Основной задачей является исследование возможности оценки возобновления леса на сплошных вырубках через посадку лесных культур по временному ряду вегетационных индексов NDVI и SWVI. Временной ряд построен по безоблачным летним композитным изображениям на основе данных Landsat 2002-2015 гг. В качестве опорных данных мы использовали лесоустроительные планы и таксационные описания 2003 года Холмченского участкового лесничества, которое входит в состав Суземского лесничества. В результате была получена линейная зависимость между значениями вегетационных индексов, их временной динамикой и таксационными характеристиками лесов (тип леса, породный состав, группы возраста). Это позволило с точностью 80,1% детектировать восстановленные после сплошных рубок леса Стеклянского, Краснослободского и Кокоревского участковых лесничеств, также входящих в состав Суземского лесничества Брянской области. Результаты работ полезны для практического применения в лесном хозяйстве, где остро стоит вопрос о мониторинге мероприятий по лесовосстановлению и лесоразведению.

Ключевые слова: лесовосстановление, временной ряд, вегетационные индексы, NDVI, SWVI.

Современная практика ведения лесного хозяйства основана на принципах устойчивого управления лесами, которые предполагают наличие оперативных актуальных знаний о состоянии леса. Основным источником данных о лесах служат материалы лесоустройства. На Европейской части Российской Федерации лесоустроительные работы проводятся с интервалом в 10-15 лет (Приказ Минприроды…, 2018). В свою очередь, при интенсивном ведении лесного хозяйства, где основной антропогенной нагрузкой на леса, как на природный ресурс, являются вырубки с целью заготовки древесины, такой интервал обновления материалов лесной таксации может быть слишком большим.

В соответствии с лесным кодексом РФ такие рубки допускаются только при условии последующего воспроизводства лесов на данных участках (Лесной кодекс…, 2006). Лесовосстановление осуществляется методами естественного, искусственного или комбинированного восстановления лесов (Лесной кодекс…, 2006). Согласно Правилам лесовосстановления, площади лесных участков, на которых проведено искусственное или комбинированное лесовосстановление, принято относить к землям, занятым лесными насаждениями, при достижении лесными растениями параметров главной лесной породы: возраст, количество деревьев, высота. Так, например, для сосны обыкновенной возраст составляет 9 лет (Приказ Минприроды…, 2016). При интервале лесоустроительных работ в 15 лет данные теряют свою актуальность. Поэтому для мониторинга мероприятий по лесовосстановлению, как естественному, так и через посадку лесных культур, а также для актуализации материалов лесной таксации и лесных карт необходим независимый метод оценки состояния древостоев на вырубках в период между проведением лесоустроительных работ.

Такой метод может быть реализован при использовании дистанционных данных, которые находятся в открытом доступе, как, например, данные Landsat (NASA, США) и другие.

В настоящее время среди огромного количества работ по изучению динамики растительности по космическим данным можно выделить несколько основных проектов: Global Forest Change (GFC) (Hansen 2013) и The Global Land Analysis and Discovery (GLAD) (Potapov 2015). Эти проекты используют для анализа изменения растительного покрова данные со спутников Landsat, которые характеризуются средним пространственным разрешением – 30 м. Также накоплен огромный временной архив данных (более 30 лет).

Результатами работ данных проектов являются спутниковые тематические продукты, отражающие прирост и сокращение леса. Точность результатов оценки, указанных выше характеристик, составила 87% для сокращения и 76% для прироста лесного покрова по GFC, и более 90% по GLAD продуктам.

В то же время у этих продуктов есть ряд недостатков:

1. в проекте GFC есть только данные по суммарному приросту по состоянию на 2012 г.;

2. проект GLAD реализован на территорию Восточной Европы;

3. не связаны с российским лесоустройством.

Первые исследования в этом направлении проводились авторами статьи для лесов Красноярского края в 2014 году (Белова и др., 2015). Эта статья является продолжением работы по исследованию возможности детектировать зарастание территорий вырубок молодняком на территории Красноярского края. Основной целью настоящей работы является исследование возможности дистанционной оценки возобновления леса на сплошных вырубках через посадку лесных культур и создание на базе полученных знаний автоматизированного метода. Результаты работ полезны для лесного хозяйства в части совершенствования методов государственного мониторинга воспроизводства лесов.

МАТЕРИАЛЫ И МЕТОДЫ

Регион исследования

В качестве тестового участка выбрана территория Суземского лесничества, расположенного на юго-западе Брянской области (рис. 1). Данная территория относится к району хвойно-широколиственных (смешанных) лесов европейской части РФ. В Суземском лесничестве, как и на всей Европейской части России, ведется активная лесохозяйственная деятельность. Работы по восстановлению проводят преимущественно через посадку лесных культур. По данным Росстата на 2012 г. лесовосстановление составило 3453 га, из них искусственное восстановление 3138 га. На вырубках в основном высаживаются культуры сосны и ели. В отличие от лесов Западной и Восточной Сибири, естественное восстановление растительности встречается редко. Основными лесообразующими породами тестового участка являются сосна, ель, береза и дуб, реже – ольха черная и осина.

Для исследования выбрано Холмченское участковое лесничество, входящее в состав Суземского лесничества.

 

Исходные данные

Для исследования использовали наземные опорные данные: лесоустроительные планы и таксационные описания выбранного участкового лесничества по состоянию на 2003 г.

Основные таксационные характеристики, которые использовались нами в ходе работ: площадь выдела, породный состав, возраст и тип леса.

В качестве дистанционных данных использовали временные серии летних безоблачных композитных спутниковых изображений. Композитные изображения были построены на основе данных Landsat-5/7/8 за 2002-2015 гг. (Белова и др., 2011). Временной ряд не является непрерывным. Летние композиты не были построены за 2003, 2004 и 2012 гг. из-за отсутствия необходимого количества спутниковых данных (изображений). Каждое композитное изображение состоит из трех спектральных каналов: красный (630-690 нм), ближний инфракрасный (770-900 нм) и средний инфракрасный (1550-1750 нм) с пространственным разрешением – 30 м.

Методика

По композитным изображениям для каждого года из временного ряда рассчитали вегетационные индексы NDVI (Normalized Difference Vegetation Index) и SWVI (Shortwave Vegetation Index) (1, 2):

NDVI=(NIR-RED)/(NIR+RED), (1)

SWVI=(NIR-SWIR)/(NIR+SWIR), (2)

где RED – красный спектральный канал (630-690 нм);

NIR – ближний инфракрасный спектральный канал (770-990 нм);

SWIR – средний инфракрасный спектральный канал (1550-1750 нм).

По лесоустроительному плану и таксационным описаниям на территорию Холмченского участкового лесничества были отобраны выделы с категориями земель: сомкнувшиеся лесные культуры разного возраста и состава, культуры в стадии восстановления (несомкнувшиеся на момент проведения таксации) разного состава и возраста, а также вырубки различной давности. По составу выделы с сомкнувшимися лесными культурами были поделены на 3 класса: ель, сосна, лиственные. К классу «лиственные» отнесли выделы, на которых произрастает дуб, осина или береза.

Для дальнейшей работы выделы были отфильтрованы по площади. Минимальная площадь была принята равной 1 га, что соответствует фрагменту изображения площадью 9 пикселов. Всего было выбрано 800 выделов. Границы выделов были векторизованы. Для исключения влияния пограничных пикселов, спектральная яркость которых формируется не только из яркостей растительности внутри выдела, но и из яркости растительности на его границе, внутри каждого выдела была построена буферная зона размером в один пиксел (рис. 2). Все пикселы каждого выдела за исключением буферной зоны рассматривались как один пространственный объект исследования.

 

Для каждого выдела рассчитали средние значения вегетационных индексов за каждый год (3):

𝑀_𝑥=(Σ₁ⁿx_i)/n, (3)

где x – значение вегетационного индекса i пиксела внутри выдела;

n – количество пикселов в выделе;

i – номер пиксела в выделе.

Далее для каждого выдела мы определили аккумулированные (накопленные) средние значения вегетационных индексов за каждый год (4), которые поместили в двумерное пространство SWVI—NDVI с целью получения множества аккумулированных кривых, характеризующих состояние растительности каждого выдела.

V_годi=Σⁱ₂₀₀₂(M_x)_i, i=2002, 2005, 2006…2015,, (4)

Где i – год;

M_x – среднее значение вегетационного индекса в выделе за i-год.

Например, для 2007 г. аккумулированные средние значения вегетационных индексов для выдела имеют следующий вид:

На рисунке 3 представлены аккумулированные кривые для девяти выделов.

 

РЕЗУЛЬТАТЫ И ИХ ОБСУЖДЕНИЕ

Форма связи аккумулированных значений вегетационных индексов лесных культур на вырубках носит линейный характер (рис. 3), а состояние растительности в выделе может характеризоваться значением угла наклона линейной функции: . Чем угол наклона больше, тем выше вероятность того, что в выделе либо отсутствует древесная растительность, либо она находится не в сомкнувшемся состоянии.

С помощью восстановленных связей можно рассчитать пороговые значения угла наклона для сомкнувшихся лесных культур (без разделения на породный состав). Для этого было решено уравнение вида:

Y = A × X + B, где B = 0,  (7)

В результате:

A = 1.829,  R² = 0.976

CKO = 0.126,

На графике (рис. 4) в виде точек отображены аккумулированные кривые для всех отобранных выделов с пороговыми линейными функциями:

Y = A × X,  (8)

Y = (A + 2CKO) × X,  (9)

Y = (A — 2CKO) × X,  (10)

Как видно из графика (рис. 4), большая часть точек, принадлежащих аккумулированным кривым выделов с сомкнувшимися лесными культурами, находится внутри интервала [A−2CKO;A+2CKO]. Таким образом, пользуясь данным интервалом – эталонным сектором, можно многократно определять состояние растительности на вырубках, где проводились мероприятия по лесовосстановлению, в любой момент времени, используя актуальные данные спутника Landsat.

 

Исследования также показали, что угол наклона линейной функции не зависит от породного состава насаждений, возраста насаждений и типа леса (рис. 5). Значения коэффициента линейной регрессии A для разных древесных пород близки, что не позволяет использовать предложенный метод для определения породы и возраста восстановленной лесной культуры.

Проверка полученных результатов

Апробацию метода проводили для тестовых выделов Стеклянского, Краснослободского и Кокоревского участковых лесничеств (входят в состав Суземского лесничества) с использованием пороговых значений пороговые значения коэффициента , полученных по данным анализа вегетационных индексов выделов Холмченского участкового лесничества.

На основе данных лесной таксации 2003 года отобрали 549 выделов площадью более 1 га с сомкнувшимися культурами разного возраста и состава, а также с несомкнувшимися культурами разного возраста и состава, сформировавшихся в результате вырубок разной давности. Как и для выделов Холмченского участкового лесничества, чтобы устранить влияние яркости пограничных пикселов выдела, внутри каждого выдела была построена буферная зона, исключенная из последующего анализа.

На территорию участковых лесничеств построили летние безоблачные композитные изображения за период 2002 – 2015 гг., по ним рассчитали вегетационные индексы NDVI и SWVI (1, 2) за указанный период.

По формуле (3) для каждого из 549 выделов были рассчитаны ежегодные средние значения вегетационных индексов. Затем для каждого выдела были построены аккумулированные кривые по формуле (4).

Из-за отсутствия доступных летних композитных изображений 2003 и 2004 гг., классификацию выделов проводили по спутниковым изображениям 2005 года, ближайшего по времени к материалам лесоустройства. По полученным пороговым значениям коэффициента линейной регрессии Холмченского участкового лесничества мы провели классификацию 549 выделов на классы: 1) сомкнувшиеся культуры разного породного состава и возраста; 2) несомкнувшиеся культуры и вырубки.

Результаты классификации мы сравнивали с таксационными описаниями этих участковых лесничеств. Результаты сравнения приведены в таблице 1.

Таблица 1 Проверка результатов классификации

Анализ ошибок показал, что точность классификации состояния лесовосстановления зависит от площади выдела (рис. 6), так как максимальное количество ошибок приходится на выделы площадью менее 2 га. Вероятно, это связано с пространственным разрешением прибора спутника Landsat (30 м или 0.09 га), который формирует изображение выделов менее 2 га со множеством смешанных» пикселей на границе лес – не лес. Возможно при использовании спутниковых изображений типа Sentinel-2 (10-20 м) позволит снизить величину ошибки для участков лесовосстановления площадью менее 2 га.

Сравнение полученных результатов с предыдущими работами

Аналогичные работы авторы статьи проводили на территории Усольского и Чунского лесничеств Красноярского края (Белова и др., 2015). В отличие от Суземского лесничества Брянской области в лесах Усольского и Чунского лесничеств Красноярского края преобладает естественное восстановление лесов. По данным Федеральной службы государственной статистики за 2012 г. искусственное восстановление на территории Красноярского края составило 7217 га, естественное – 49832 га.

 

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

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

 

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

ВЫВОДЫ

В ходе проделанной работы были изучены и проанализированы лесоустроительные материалы и временной ряд спутниковых безоблачных композитных изображений Landsat. Была получена зависимость между динамикой значений вегетационных индексов и состоянием растительности на сплошных рубках Европейской части России, которая характеризует зарастание территорий сплошных рубок при искусственном лесовосстановлении.

Основные выводы по проделанной работе:

1) данный алгоритм анализа восстановительной динамики древесной растительности позволяет характеризовать переход территории вырубки из непокрытой лесом площади в покрытую как в случае естественного восстановления, так и в случае искусственного метода восстановления.

2) полученные эталонные значения коэффициента линейной регрессионной модели можно многократно использовать при детектировании смыкания лесных культур в определенных лесорастительных условиях.

3) для выделов площадью менее 2 га необходимо использовать другие методы классификации или спутниковые данные более высокого пространственного разрешения, т.к. при применении модели для изображений Landsat в более чем 70% случаев метод классифицирует их ошибочно.

4) данный алгоритм не чувствителен к породному и возрастному составу насаждения. С его помощью можно фиксировать факт перехода вырубки из непокрытой в покрытую лесом территорию без возможности оценки породно-возрастной структуры древостоя.

БЛАГОДАРНОСТИ

Исследование выполнено за счет гранта Российского научного фонда (проект №19-77-30015). Подготовка и анализ тематических спутниковых изображений проводилась в рамках тематики научных исследований ЦЭПЛ РАН А18-118052590019-7).

СПИСОК ЛИТЕРАТУРЫ

Белова Е.И., Ершов Д.В. Опыт оценки естественного лесовосстановления на сплошных вырубках по временным рядам LANDSAT // Лесоведение. 2015. № 5. С. 339-345.

Белова Е.И. Ершов Д.В. Предварительная обработка временных серий изображений Landsat—TM/ETM+ при создании безоблачных композитных изображений местности // Современные проблемы дистанционного зондирования Земли из космоса. 2011. Т. 8. № 1. С. 73-82.

Приказ Минприроды России от 29.03.2018 № 122 «Об утверждении Лесоустроительной инструкции». URL: http://www.consultant.ru/document/cons_doc_LAW_296757/77e0bc4c1b732cffc4ef23c2113effb964769f33/ (дата обращения 28 августа 2019).

Лесной кодекс Российской Федерации от 04.12.2006 №200-ФЗ (ред. От 29.12.2017). URL: http://www.consultant.ru/ document/cons_doc_LAW_64299/710f2a2a776bb8cff9e204556177d7cb26f79dee/ (дата обращения 28 августа 2019).

Приказ Минприроды России от 29.06.2016 № 375 «Об утверждении Правил лесовосстановления». URL: http://www.consultant.ru/document/cons_doc_LAW_207285/4f0faeaf51f7cf71d1e928c355da9b3db3d3545d/ (дата обращения 28 августа 2019).

Федеральная служба государственной статистики 2018. URL: http://www.gks.ru/wps/wcm/connect/rosstat_main/rosstat/ru/statistics/environment/# (дата обращения 28 августа 2019).

Hansen M.C., Potapov P.V., Moore R., Turubanova S.A., Tyukavina A., Thau D., Stehman S.V., Goetz S.J., Loveland T R., Kommareddy A., Egorov A., Chini L., Justice C.O., Townshend J.R.G. High Resolution Global Maps of 21st Century Forest Cover Change // Science. 2013. Vol. 342. P. 850-853.

Potapov P.V., Turubanova S.A., Tyukavina A., Krylov A.M., McCarty J.L., Radeloff V.C., Hansen M.C. Eastern Europe’s forest cover dynamics from 1985 to 2012 quantified from the full Landsat archive // Remote Sensing of Environment. 2015. Vol. 159. P. 28-43.

Рецензент: к.б.н. Медведева М.А.