Исследование Земли из Космоса, 2021, № 5, стр. 72-84

ВВЕДЕНИЕ

Развитие современных методов по дешифрированию данных спутниковой съемки позволяет получать информацию о состоянии лесов, их использовании, количественных и качественных показателях. Сочетание подходов по комплексному анализу оптической и радиолокационной съемки за более чем полувековую историю стало одним из эффективных способов получения информации о лесах на труднодоступных территориях. В отличие от оптической съемки, методов дешифрирования радиолокационной съемки можно считать их независимость от облачности атмосферы, погоды, освещенности. Несмотря на значительное развитие данного направления, до настоящего времени не существует эффективных методов определения таксационных показателей насаждений по спутниковым радарным съемкам. Данная тенденция объясняется отсутствием единых подходов, разнообразием применяемых методов анализа спутниковой съемки, региональной спецификой лесов, недостаточными объемами экспериментальных исследований. С целью практического внедрения различных алгоритмов таксации лесных насаждений по радарным данным, прежде всего, требуется определение оптимальных методов анализа информации по взаимосвязи показателей съемки с таксационными характеристиками насаждений и разработки на их принципах методов таксации лесов. В ранее проведенных исследованиях ИФМ СО РАН и компанией “СОВЗОНД” в районе дельты р. Селенги (Бурятия) по применению радиолокационной съемки для решения задач лесного хозяйства получены результаты, позволяющие определить высоту полога леса, полноту, покрытые и не покрытые лесной растительностью земли (Чимитдоржиев и др., 2014). В ходе проведенного эксперимента не удалось определить весь спектр таксационных характеристик насаждений, тем не менее, результаты исследований показали возможность применения данных радиолокационной съемки X диапазона для оценки высоты полога и сомкнутости леса. В последующем в 2016 г. в Германском центре авиации и космонавтики были проведены научно-исследовательские работы по разработке алгоритмов классификации покрытых и не покрытых лесной растительностью земель на основе спутниковых данных TanDEM-Х. В результате значительного объема проведенных исследований создана глобальная карта покрытых и не покрытых лесом территорий. Точность полученной карты оценивается более 90% (Martone et al., 2018). Значительные объемы исследований по оценки качественных и количественных характеристик лесов на основе радиолокационных данных проведены ООН на базе проекта National Forest Assessment – NFA (“Национальная оценка лесных ресурсов”) в рамках глобальной программы Sustainable Forest Managementin Changing Climate (“Стабильное лесоводство в изменяющемся климате”) (Koch, 2013). Работы показали перспективность использования радарных данных для определения запаса насаждений, лесопокрытой площади, сомкнутости.

Использование радиолокационных данных представляет значительный интерес также при изучении растительного покрова, ранее проведенные исследования в этом направлении показывают наличие взаимосвязи характеристик растительности с текстурными показателями радиолокационной съемки (Бондур, Чимитдоржиев, 2008), а также зависимости пространственной анизотропии неоднородностей лесной растительности при различных азимутальных углах радарного поляриметрического зондирования (Бондур, Чимитдоржиев, 2019). Учитывая перспективность применения данных радиолокационной съемки при различной облачности, специалистами ФГУП “Рослесинфорг”, “Центрлеспроект”, а также компании “Совзонд” для выявления вырубок проводился комбинированный анализ оптической съемки SPOT и радарного композита MTC. Далее из них составлялся мультивременной оптико-радарный композит, на котором новые вырубки автоматически выделялись розовым цветом и затем оконтуривались вручную. На основе работы были сделаны выводы, что методика выделения вырубок с помощью оптико-радарных мультивременных композитов аналогична выявлению изменений по оптическим мультивременным композитам. Радарные данные, за счет того, что при съемке можно не брать в расчет облачность и освещенность, позволяют гарантированно и оперативно получить необходимую информацию о вырубках к нужной дате с незначительной потерей точности картирования границ вырубки (Бахтинова и др., 2012). Значительный интерес представляет опыт реализации глобального проекта BIOMASAR по оценке биомассы бореальных лесов России, Северной Европы, Канады и Аляски, осуществленный с использованием нового алгоритма обработки данных, полученных с датчика Advanced Synthetic Aperture Radar (ASAR) спутника Европейского космического агентства (ESA) EnviSat. В ходе работ по осуществлению проекта разработаны методы дистанционного определения арктической лесной биомассы (New boreal forest…, 2010; Santoro et al., 2010). Проект BIOMASAR, подтвердил эффективность применяемого алгоритма обработки данных радарной съемки. Дистанционное определение объема роста запаса насаждений (GSV) стало возможным на уровне от 500 м³/га для бореальных лесов при использовании данных EnviSat. Для отработки алгоритма, получившего название ScanSAR, были заложены несколько пробных площадей в Скандинавии, Сибири и Канаде. Используя данный метод и архивы данных EnviSat, стало возможным издавать карты запасов насаждений ежегодно, при этом разрешение карт составляет 10 км и выше, а погрешность измерений не превышает 15% (New boreal forest…, 2010).

Развитие современных методов таксации в том числе с использованием данных спутниковой съемки имеет важное значение в стратегии планирования и ведения лесного хозяйства, а также в совершенствовании различных технологий лесоустройства (Рафаилов, Мартынюк, 2016; Моисеев и др., 2017).

Основной целью работ являлось проведение экспериментальных исследований в таежной зоне Западной Сибири, разработка алгоритмов определения таксационных характеристик лесов с использованием данных радарной съемки.

ОБЪЕКТЫ И МЕТОДЫ ИССЛЕДОВАНИЙ

Исследования проведены в ленточных борах Алтайского края. В связи с жесткими почвенно-климатическими условиями, основной лесообразующей древесной породой на территории района исследований является сосна обыкновенная, доля которой в составе лесов в среднем составляет 85% и меняется с 73% в северо-восточной части до 92% в юго-западной. В сложившихся лесорастительных условиях ленточных боров другие древесные породы имеют малое распространение. В пониженных элементах рельефа, где почва обладает большей влажностью, произрастают мягколиственные породы – береза и осина. Основной объем исследований проведен в Кулундинской, Барнаульской, Космолинской лентах.

При проведении исследований использовались актуальные данные лесоустройства, давность которого не превышала 5 лет. Для получения взаимосвязей таксационных характеристик насаждений с данными радиолокационной съемки были обработаны 204 космических снимков с космических сенсоров Sentinel-1 типа GRDH и SLC. Период охвата съемки – с 02.01.2018 по 28.12.2018 гг. (табл. 1). Обработка снимков проводилась с использованием программного обеспечения “SNAP”, а также отдельных пакетов программных комплексов Matlab, Rstudio и скриптов на языке программирования Python (Руководство пользователя…; Методические рекомендации…, 2020). Общий объем продуктов, полученных по результатам обработки и получения характеристик, используемых при лесотаксационном дешифрировании, составил более 2.45 Тб.

Для работы по оптической съемке использованы данные сенсора Sentinel-2 в количестве 10 снимков уровня обработки L1C зимнего периода. Период охвата съемки – с 01.12.2018–03.01.2020 гг. (табл. 2). Для снимков проводилась предварительная обработка с использованием алгоритмов DOS-коррекции с помощью модуля Semi-Automatic Classification программного комплекса QuantumGIS (QGIS). Общий объем обработанной информации по оптической съемке составил более 72 Гб.

Используемые в исследованиях показатели спутниковой съемки устанавливались опытным путем, на основе их статистических взаимосвязей с таксационными характеристиками насаждения. По степени достоверных взаимосвязей выделялись показатели, которые в дальнейшем были использованы для моделей зонирования территории района исследований по запасу и полноте насаждений. При анализе оптической съемки установлены взаимосвязи запаса и полноты насаждений с красным и инфракрасным каналом. Объединение этих каналов и, как следствие, увеличение статистической взаимосвязи возможнo при использовании индекса NDVI. При этом учитывались взаимосвязи данного индекса с запасом насаждений и полнотой в различные времена года. Лучшие результаты были достигнуты при использовании зимних снимков (Сидоренков и др., 2019).

Анализ радиолокационных данных осуществлялся по 22 показателям, на основе достоверных статистических зависимостей с таксационными характеристиками насаждений выделено 14 показателей, которые в дальнейшем и использовались для анализа данных. Перечень и расшифровка показателей радиолокационной съемки, используемых в работе, приведены в табл. 3.

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

В ходе полевых исследований заложено 134 пробных площади в различных насаждениях по породному составу, возрасту. Незначительное разнообразие по породному составу лесов Алтайского края позволило отработать алгоритмы определения таксационных характеристик лесных насаждений по данным радиолокационной съемки с сенсора Sentinel-1. Дифференциация насаждений осуществлялась с учетом ландшафтных особенностей местности, категорий труднодоступных участков. Первичный анализ проводился по материалам лесоустройства, определялись различия насаждений по запасам, полнотам, возрасту, планировался маршрут в связи с доступностью насаждений. Далее проводился анализ с использованием материалов оптической съемки в несколько этапов. На первом этапе осуществлялась дифференциация территории покрытых лесной растительностью земель на основе разделения по NDVI-индексу в зимний период покрытых оптическим съeмкам с сенсора Sentinel-2. Реализация данного алгоритма выполнялась через дерево решений в программном комплексе Envi 5.2. Результатом классификации являлась маска покрытых лесной растительностью земель. При этом из лесных земель исключались населенные пункты, дороги, лесозащитные полосы. Применение маски позволило выделить земли, покрытые лесной растительностью.

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

или где NDVI – (Normalized Difference Vegetation Index) – нормализованный относительный индекс растительности b8, b9 – восьмой и девятый каналы сенсора Sentinel-2

Коэффициент детерминации полученной модели R² = 0.6. Модель позволяет точно спрогнозировать запасы в средне и высокополнотных насаждениях в основном средневозрастных, спелых и перестойных классов возраста. Результаты анализа данных оптической съемки показывают, что сложность прогнозирования запаса в молодняках из-за их незначительной полноты, просветов полога.

Дифференциация района исследований по запасу насаждений по данным оптической съемки с сенсора Sentinel-2 позволила при проведении экспериментальных работ учитывать различие насаждений по данному признаку. В ходе экспериментальных работ заложено 134 пробных площади в насаждениях, различающихся по запасу, полноте, возрасту (рис. 2).

Результаты корреляционных взаимосвязей запаса, полноты, высоты насаждения с показателями радиолокационной съемки с сенсора Sentinel-1 (рис. 3, 4, 5) показывают сходство зависимостей запаса и полноты насаждений по месяцам года и по показателям. Данная тенденция прослеживается и в кластерном анализе данных, в котором показатели запаса насаждений и полноты отнесены к одной группе. Такая связь объясняется природой показателей, большая полнота насаждений определяет больший запас.

Результаты исследований показывают зависимости запаса и полноты насаждений от показателей GammaVH, GammaVV, GLCMMean, Sum. При этом большие коэффициенты корреляции характерны для периода января и февраля. Видимо, такая тенденция связана с наличием устойчивого снежного покрова, что способствует маскетированию нижних ярусов растительности, контрастности выделения верхнего яруса. Ранее проведенные исследования, а также данные литературных источников показывают, что сходные показатели наблюдаются при дешифрировании запаса и полноты насаждений по оптическим снимкам (Koch, 2013).

Оценка зависимостей высоты насаждений показывает устойчивую взаимосвязь в летний период с показателями GLCMCorrelation, GLCMMean, GammaVV (рис. 5). Значительное снижение корреляции для августа месяца скорее всего характерно из-за выпадения осадков в момент радиолокационной съемки.

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

Полученные результаты облака связей (распределения) запаса и полноты насаждений с показателями радиолокационной съемки за январь месяц показывают, что для разработки моделей прогноза предпочтительно использовать показатели удельной эффективной площади рассеяния (УЭПР) объекта в значении гамма-ноль на поляризации VH и поляризации VV (GammaVH, GammaVV), радарный индекс, представляющий сумму УЭПР объектов в значении гамма-ноль на VV и VH поляризациях (Sum), текстурный признак, суммарное среднее, являющееся произведением значения пикселя с частотой его появления с соседним пикселем (GLCMMean) (рис. 3, 4). Построение моделей прогноза высоты древостоев лучше осуществлять по показателям GLCMCorrelation, GLCMMean, GammaVV в летний период времени.

Анализ взаимосвязей таксационных характеристик насаждений на пробных площадях с показателями радиолокационной съемки с сенсора Sentinel-1 позволил установить, что оптимальный период получения достоверных зависимостей запаса, полноты насаждений является период зимы, высоты насаждений – период лета. Наблюдается взаимосвязь запаса, полноты и высоты насаждений с удельной эффективной площадью рассеяния (УЭПР) объекта в значении гамма-ноль на поляризации VV, VH (GammaVH, GammaVV), а также с текстурными признаками суммарным средним (GLCMMean) и дисперсией (GLCMVariance). Проведенный анализ данных по этим показателям позволяет провести зонирование модельной территории по запасу и полноте насаждений (рис. 6, 7).

Определение запаса насаждений на модельном участке проведено с использованием его зависимости с показателями УЭПР объекта в значении гамма-ноль на поляризации VV (GammaVH), а также с текстурным признаком – суммарным средним (GLCMMean) (рис. 6). Результаты дифференциации района исследований по показателям запаса насаждений имеют сходство с данными лесоустройства. Основное отличие наблюдается в молодняках, запас большинства из них превышает 10 м³ на га и в высокобонитетных спелых и перестойных насаждениях. Различие в молодняках объясняется давностью материалов лесоустройства (превышающим 5 лет) и появлением на этих участках самосева сосны или созданием лесных культур. Полученные результаты запаса по спелым и перестойным высокобонитетным насаждениям подтверждают правильность произведенных расчетов по материалам пробных площадей. В данных условиях произрастают насаждения 3–4 бонитета, которые в большинстве случаев достигают запаса 250–350 м³. Запасы 450 м³ не характерны для данной территории и скорее всего они связаны с неточностью данных лесоустройства.

Результаты зонирования территории модельного участка по полноте насаждений (рис. 8) сходны по своей сути с зонированием территории по запасу насаждений (рис. 6). Они на большинстве участков совпадают с данными лесоустройства, за исключением высокопродуктивных насаждений по данным лесоустройства с запасом более 400 м³ га. На нарушенных территориях, зафиксированных при материалах лесоустройства более 5 лет назад, появились молодняки в основном из самосева сосны. Полнота данных насаждений достигает 0.1.

Лучшие результаты по запасу насаждения были получены с применением модели полинома третьей степени. Уравнения запаса насаждений построены на основе его взаимосвязи с показателями удельной эффективной площади рассеяния поляризации VV (GammaVV) и значением пикселя с частотой его появления с соседним пикселем при вычислении МСВ (GLCMMean). Статический анализ показал возможность использования также логарифмического уравнения взаимосвязи. Эффективность зависимости оценивается коэффициентом корреляции прогнозных значений в интервале от 0.64 по логарифмической модели и 0.7 полиноминальной (рис. 7).

Полиноминальное уравнение связи запаса насаждений с показателями GammaVV, GLCMMean

Логарифмическое уравнение связи запаса насаждений с показателями GammaVV, GLCMMean:

где М – запас насаждения, GammaVV – УЭПР объекта в значении гамма-ноль на поляризации VV, GLCMMean – текстурный признак, являющийся произведением значения пикселя с частотой его появления с соседним пикселем.

Наибольшие показатели погрешности запаса наблюдаются в средне полнотных насаждениях в диапазоне от 100 до 300 м³ и могут в ряде случаев достигать 82 м³ на га. Тем не менее, полученные результаты исследований позволяют достаточно точно проводить оценку запаса низко и высоко продуктивных лесных насаждений. Возможно, значительные погрешности в среднепродуктивных насаждениях наблюдаются из-за особенностей фильтрации шумов при предварительных алгоритмах обработки снимка. При проведении дальнейших исследований в этом направлении работы будут продолжены с целью поиска оптимальных моделей предварительной обработки радиолокационных данных спутниковой съемки с сенсора Sentinel-1, позволяющих снизить погрешность определения запаса. Улучшение модели прогноза запаса возможно также путем создания модели выбросов и их корректировки. Проведенные исследования показывают также возможность определения высокопродуктивных запасов насаждений, превышающих 300 м³ на га.

Результаты статистического анализа взаимосвязи полноты насаждений с показателями GLCMMean, GammaVV, GammaVH показали, что лучшие результаты позволяет получить модель на основе связей полинома третьей степени с показателями GammaVV, GammaVH. Коэффициент корреляции прогнозных значений полноты по данной модели составляет 0.6 (рис. 9). Средняя ошибка по относительной полноте составляет около 0.2.

Уравнения связи полноты насаждений с показателями удельной эффективной площади рассеяния поляризации VV, VH (GammaVV, GammaVH):

где Р – полнота насаждения, GammaVV – УЭПР объекта в значении гамма-ноль на поляризации VV, GammaVР – УЭПР объекта в значении гамма-ноль на поляризации VH, GLCMMean – текстурный признак, являющийся произведением значения пикселя с частотой его появления с соседним пикселем.

ВЫВОДЫ

Полученные результаты прогнозирования запаса и полноты насаждений показывают возможность использовать данные Sentinel-1 для оценки этих показателей в лесах с третьим разрядом лесоустройства, то есть части эксплуатационных и резервных лесов. Достоверные взаимосвязи запаса и полноты насаждений наблюдаются с показателями удельной эффективной площади рассеяния (УЭПР) объекта в значении гамма-ноль на поляризации VH и поляризации VV (GammaVH, GammaVV), радарным индексом, представляющим сумму УЭПР объектов в значении гамма-ноль на VV и VH поляризациях (Sum), текстурным признаком − суммарное среднее, являющимся произведением значения пикселя с частотой его появления с соседним пикселем (GLCMMean). По статистическим показателям точности лучшей моделью является полиноминальное уравнение запаса (1) и полноты насаждений (2). Исследования показывают, что лучшим периодом для разработки моделей оценки запаса и полноты насаждения является январь и февраль месяцы. Видимо, такая тенденция связана с маскетированием снегом нижних ярусов растительности, что снижает шумы сигнала и как следствие приводит к большей точности уравнений взаимосвязи.

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