CC BY
0
Вычислительные технологии, 2022, том 27, № 2, с. 105-121. © ФИЦ ИВТ, 2022 ISSN 1560-7534
Computational Technologies, 2022, vol. 27, no. 2, pp. 105-121. © FRC ICT, 2022 elSSN 2313-691X
ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ
D01:10.25743/ICT.2022.27.2.009
Оптико-микроволновая диагностика восстановления леса после пожаров
A.B. Дмитриев*, Т.Н. Чимитдоржиев, П.Н. Дагуров
Институт физического материаловедения СО РАН, 670047, Улан-Удэ, Россия
*
Поступила 21 января 2022 г., принята в печать Ц- февраля 2022 г.
Рассмотрена возможность обнаружения лесного подроста с помощью поляриметрических разложений и анализа временных рядов радиолокационных изображений ALOS-2 PALSAR-2 и Sentinel-1. Представленные результаты показывают принципиальную возможность уверенно определять рост лесных насаждений.
Ключевые слова: лесовосстановление, поляриметрические разложения, анализ временных рядов, ALOS-2 PALSAR-2, Sentinel-1/2.
Цитирование: Дмитриев A.B., Чимитдоржиев Т.Н., Дагуров П.Н. Оптико-микроволновая диагностика восстановления леса после пожаров. Вычислительные технологии. 2022; 27(2):105—121. D01:10.25743/ICT.2022.27.2.009.
Введение
Проблема масштабной оценки процессов восстановления леса после сплошных вырубок и лесных пожаров является актуальной в связи с климатическими изменениями и соответствующим стремлением мирового сообщества к дальнейшему низкоуглеродному развитию. Разработаны различные модели, прогнозирующие положительный или отрицательный углеродный баланс лесов к 2100 г. [1, 2]. Создаются глобальные карты, характеризующие распределение биомассы леса на основе данных дистанционного зондирования Земли (ДЗЗ) [3, 4]. Выполняются различные оценки поглощающей способности леса. Показано отсутствие значимых изменений для углеродного стока до и после прореживания леса — после выборочных рубок [5]. Эти исследования еще раз подтвердили теорию о существовании порога насыщения, после которого поглощение углерода лесными экосистемами изменяется незначительно. В [6] установлено, что сток углерода в зрелых тропических и бореальных лесах Сибири составляет в среднем 0.85 Пг/год (1 Пг = 1015 г = 1012 кг), а для лесного подроста 1.30 Пг/год. В [7] на примере Латинской Америки также подтверждается значимая роль лесного подроста: молодой лес после сплошных вырубок позволяет активно накапливать углерод в надземной биомассе, тем самым частично уравновешивая выбросы углерода в результате сжигания ископаемого топлива, деградации лесов и других антропогенных источников.
Подобные результаты исследований свидетельствуют о необходимости выполнения масштабной дистанционной оценки площадей лесовосстаповлепия после сплошных рубок и сильных лесных пожаров. Это в свою очередь предполагает разработку новых
и развитие существующих методов дистанционного зондирования процессов лееовое-становления. Одним из эффективных способов дистанционной оценки биофизических параметров растительных ресурсов является комплексное использование оптико-микроволновых данных [8-10],
В настоящей работе рассматривается возможность выявления подроста лесных насаждений при помощи анализа мультивременных данных радиолокатора с синтезированной апертурой (РСА) АЬОБ-2 I А1.ЯАН-2 с полной и двойной поляризацией (НН и НУ), Для сравнения также представлены результаты, полученные по временным рядам данных РСА С-диапазона Яеш то!-1 с двойной поляризацией (АА' и УН),
1. Методика исследований
Исследования основаны на сравнении результатов выделения подроста лесных насаждений при помощи нескольких методов,
• Поляриметрические преобразования, выделяющие основные физические механизмы рассеяния, например: объемное от лесной среды, поверхностное от шероховатой поверхности, уголковое от естественных и искусственных отражателей. Поскольку после лесных пожаров лесная растительность полностью уничтожается, предполагается наличие двухекачкового рассеяния (например, отражение от земной поверхности и остатков стволов деревьев) и поверхностного. После самосева и роста новых деревьев предполагается увеличение компоненты объемного рассеяния и уменьшения вкладов поверхностного и уголкового рассеяния,
• Анализ изменения интенсивности обратного радарного рассеяния и интерферо-метрической когерентности на различных поляриметрических комбинациях и в различных диапазонах. Ожидается, что рост лесных насаждений приведет к росту амплитуды радарного сигнала вследствие увеличения эффективной "шероховатости" отражающей поверхности, В то же время должно происходить уменьшение величины интерферометричеекой когерентности,
• Анализ временных рядов нормализованного индекса растительности (NDVI), полученных при помощи облачной платформы Google Earth Engine [11, 12],
1.1. Поляриметрические разложения
Поляризация, наряду с длиной волны, является ключевой характеристикой радиолокационного сигнала, распространяющегося в растительной среде и рассеивающегося на отдельных ее элементах. Большинство радаров спроектированы на передачу горизонтально (Н) или вертикально (V) поляризованных электромагнитных волн. Аналогично антенна радара может также принимать горизонтально или вертикально поляризованный отраженный сигнал. Некоторые радары могут принимать обе поляризации одновременно, В общем виде поле, рассеянное средой Ег, связано с падающим полем Ег соотношением
Акг
-[5 ]E
г
или в матричной форме
Егн Erv
Акг
Shh Shv Svh Svv
Eh ev
r
Здесь матрица ] называется матрицей рассеяния, а ее элементы = —
комплексные коэффициенты рассеяния при излучении па ^'-поляризации и приеме па г-поляризации, Множитель егкг/г, где к = 2п/\ — волновое число, определяет набег фазы и ослабление при распространении отраженной сферической волны от цели до радара, В радарной поляриметрии данная абсолютная фаза, как правило, не учитывается. Отметим, что при моностатическом зондировании выполняется условие 8ну = а матрица [5] становится симметричной,
С помощью матрицы рассеяния можно полностью описать детерминированные (точечные) рассеиватели, которые изменяют поляризацию падающей волны. Однако, для распределенных рассенвателей сделать невозможно. Поэтому такие рассеиватели зачастую представляют в виде большого числа случайно распределенных точечных целей. Тогда матрица рассеяния состоит из когерентной суммы всех матриц рассеяния данных элементарных рассенвателей. Таким образом, для анализа поляризационного состояния отраженного радарного сигнала необходимы статистические методы.
Наиболее общим подходом для полного описания распределенных рассенвателей является матрица когерентности [Т] размером 3x3, Она определяется как внешнее произведение трехмерных векторов рассеяния кр, записанных в базисе Паули
Здесь угловые скобки (•) обозначают операцию пространственного усреднения (некогерентного накопления, тиШ-1оо]«§), а знак "*" — комплексное сопряжение. По построению матрица [Т] — эрмитова и имеет положительные собственные значения. Ее вне-диагональные элементы показывают степень коррелированное™ элементов вектора кр. Основной целью поляриметрических разложений является представление отраженного от распределенной цели радиолокационного сигнала в виде вкладов элементарных механизмов рассеяния. Поляриметрические разложения широко применяются для интерпретации, классификации и сегментации радиолокационных данных. Как правило, они подразделяются на два класса: разложения, основанные на анализе собственных значений и собственных векторов матрицы когерентности [Т], и разложения, основанные на различных моделях зондируемой среды. Приведем кратко основные элементы теории,
1.1.1. Разложение, основанное на анализе собственных значений и собственных векторов (разложение Клауда^Потье)
Как было отмечено выше, матрица когерентности [Т] по построению самосопряженная (см, уравнение (2)), Следовательно, она может быть приведена к диагональному виду
¿яя + ¿уу
и имеет вид
[Т] = [ и ][Л][ и *т ],
где
Л 0 0 eii ei2 ei3
Л = 0 А2 0 , и = e2i е22 е23
0 0 A3 б3! е32 е33
Л — диагональная матрица неотрицательных собственных значений матрицы [Т], 0 ^ А3 ^ А2 ^ Ai, a [U] = [ei e2 e3] — унитарная матрица ортонормированных собственных векторов. Каждый из этих векторов может быть параметризован пятью углами [13, 14]
cos aj exp(i5ij) sin aj cos fy exp(íí2j-) sin aj sin fyj exp(i53j)
которые характеризуют различные механизмы рассеяния.
Идея, лежащая в основе данного разложения, называемого также декомпозицией Клауда-Потье, связана с представлением матрицы когерентности [Т] в виде некогерентной суммы трех независимых (т.е. ортогональных) матриц [Ti]:
[Т] = [U][Л][^*т] = [Ti] + [Т2] + [Тз] = Ai (eiei*T) + А2 (е2в2*т) + A3 (e^) . (3)
В качестве доминирующего механизма рассеяния в модели принимается механизм, соответствующий собственному вектору с максимальным собственным значением Ai; остальные два вектора описывают второй и третий по частоте встречаемости механизмы рассеяния. Из собственных значений матрицы когерентности можно получить два важных статистических параметра: энтропию Н и анизотропию А. Энтропия описывает степень случайности рассеяния
Н = Pi lo§3Pi, Pi
Ai
i=i
£ Ai
(4)
i=i
где Рг — вероятность каждого из трех механизмов. Если только одно из собственных значений матрицы когерентности отлично от нуля, с вероятностью единица наблюдается единственный механизм рассеяния (детерминистический случай) и Н = 0, Если, с другой стороны, выделенного механизма рассеяния нет и А1 = А2 = Аз, то Н = 1. Однако оба этих крайних случая в реальности не встречаются, обычно наблюдается промежуточное положение вещей. Приближением к Н =1 является рассеяние на облаке хаотически ориентированных диполей, которым моделируют кроны плотных лесов. Высокие значения энтропии указывают, как правило, на наличие растительности.
Анизотропия рассеяния определяется как нормализованная разность второго и третьего собственных значений
^ = А2 — Аз А2 + Аз
Она также изменяется в диапазоне от 0 до 1 и выражает отношение между вторым и третьим механизмами рассеяния. Высокие значения анизотропии, как правило, указывают на преобладание второго по частоте встречаемости механизма рассеяния. Анизотропия выступает в качестве дополнительного к энтропии параметра и позволяет интерпретировать механизмы рассеяния.
Физическая и геометрическая интерпретация механизмов рассеяния, представленных матрицами когерентности [Т^, [Т2] и [Т3] (3), может быть получена из соответствующих собственных векторов. Для каждого собственного вектора можно определить угол рассеяния aj = агеео8(|еу |), который лежит в диапазоне от нуля до 90° и связан с соответствующим механизмом рассеяния. Средневзвешенная величина этих углов дает еще один параметр в декомпозиции Клауда-Потье
а = р1«1 +р2а2 + Рз®з. (5)
Угол а принимает значения от нуля (соответствует механизму рассеяния от поверхно-
°
°
1.1.2. Модельное поляриметрическое разложение
Модельные разложения ставят своей целью представить рассеянное поле в виде суммы предопределенных (известных) компонентов рассеяния. Одним из первых таких модельных разложений было трехкомпонентное разложение Фримана- Дурдена [15]. Оно основано на трехкомпонентной модели мощности матрицы рассеяния. При этом общая мощность обратного рассеяния Р предполагается состоящей из суммы трех независимых слагаемых: мощности однократного (поверхностного) рассеяния, мощности эхо-сигнала, претерпевшего двукратное переотражение структурами, аналогичными двугранному уголковому отражателю, и мощности объемного (многократного) рассеяния:
Р = |^яя|2 + 2 |^Яу |2 + |2 = Р3 + Ра + ,
ветственно. Однократное поверхностное рассеяние — типичный механизм рассеяния от почвы или воды. Двойное переотражение встречается при наличии вертикально ориентированных целей на горизонтальной поверхности, оно часто встречается на застроенных территориях (отражение от земли и стен), а также в лесах (отражение от стволов и почвы). Объемное рассеяние встречается там, где имеется слой случайно ориентированных элементарных отражателей — например, ветвей деревьев, В модели используется корреляция между элементами матрицы рассеяния (1), описываемая ковариационной матрицей
(^яя1 ) ( л/25ЯЯ5ЯУ) (3нн3уу) [С]= (л/23ну3*нн) (2 |5Яу|) (л/28НуЗ*уу) . (6)
_ (ЗууЗ*нн) (л/25уу5яу) (^уу| ) _
Модель предполагает статистическую некоррелированность согласованной и перекрестной компонент матрицы рассеяния (1): ( ЗннЗ*Ну) ~ (ЗууЗ*Ну) ~ 0, а кроме того, ее симметричность Зну = Зун- Используется также ряд допущений, связанных с каждым из механизмов рассеяния в отдельности — таких, как равномерное распределение ориентации элементарных отражателей диполей в объемном слое, применимость брэг-говской модели к поверхностной составляющей рассеяния, диагональноеть матрицы рассеяния, характеризующей двойное переотражение. После параметризации элементов матрицы рассеяния для каждого из трех механизмов, использующих названные ограничения, авторы разложения получают набор уравнений
(\Shh |2> = \b\2fs + и2 fä + л, (\Shv I2) = f-,
(\Svv\2> = fs + fä + fv, {Shh S*v ) = bfs + afä + у,
где fä и fv — вклады поверхностного, двукратного (уголкового) и объемного рассеяния соответственно; а и b — степень корреляции горизонтальной и вертикальной компонент для двукратного и поверхностного рассеяния соответственно.
Таким образом, имеются четыре уравнения с пятью неизвестными. Последнее допущение вводится для того, чтобы получить единственное решение системы: в зависимости от знака действительной части компоненты (ShhS*v) доминирующим механизмом рассеяния объявляется либо поверхностное (и тогда a = —1), либо двукратное (и тогда b =1). Наилучшим образом модель работает в условиях, когда вклад либо поверхностного, либо двукратного переотражения невелик. Следует отметить, что второе часто имеет место для естественных покровов в L-диапазоне, И, наконец, после того как все неизвестные найдены, окончательное разложение мощности Р на три составляющие Ps, Pä и Pv производят по формулам
Ps = fs (1 + \b\2), Pä = fä(1 + H2), Pv = 8^.
1.2. Амплитудные радиолокационные и оптические данные
Широко используемым подходом к задаче определения параметров земного покрова является измерение и анализ коэффициента обратного рассеяния, характеризующего интенсивность радиолокационного сигнала. Применительно к проблеме лесовосстанов-ления анализ вариаций интенсивности обратного рассеяния может дать информацию о росте растительности. Так, по мере роста лесного подроста увеличиваются его биомасса, количество ветвей, усложняется геометрическая структура, В целом это приводит к увеличению обратного рассеяния и заметным эффектам деполяризации сигнала.
Для корректной интерпретации отраженных радиолокационных сигналов необходимо устранение эффектов, обусловленных особенностями наклонного зондирования космическими радиолокаторами. Для этого помимо стандартной геометрической коррекции полученных изображений с учетом рельефа местности (т, е, ортотранеформиро-вания) необходимо проводить и радиометрическую коррекцию (нормализацию). Радиометрическая коррекция с учетом рельефа (RTC - Radiometric Terrain Correction) [16] позволяет учесть "истинную" площадь отражающей поверхности, формирующей каждый отдельный пиксел изображения и тем самым нормализовать обратное рассеяние, В данной работе используются временные ряды интенсивности, выраженные величиной 7° (gamma nought). Для ее нахождения исходные комплексные данные подвергались следующим этапам обработки:
• применение точных орбитальных данных,
• устранение теплового шума,
• радиометрическая калибровка,
• некогерентное накопление (multi-looking),
• радиометрическая коррекция с учетом рельефа местности,
• геометрическая коррекция с учетом рельефа местности,
• перевод полученных значений в децибелы.
Все этапы обработки проводились в программном обеспечении ESA SNAP 8,0 [17],
2. Описание тестового полигона и используемые данные
Тестовый полигон расположен в центральной части Республики Бурятия на территории Курбинского лесничества (рис. 1), входящего в состав Байкальского горного лесного района. Это гористая, труднодоступная местность, покрытая различными типами лесной растительности. Сильный пожар летом 2007 г. уничтожил лес на значительной части исследуемой территории. В рамках тестового полигона были выбраны четыре различных участка для детального исследования, которые обозначены цифрами на врезке рис. 1. Цифрой 1 обозначен участок общей площадью 12 га, на котором после пожара начался процесс естественного лесовосстановления. Здесь лесной подрост представляет собой самосев молодых насаждений сосны высотой от 0.7 до 1.5 м (рис. 2). Цифрой 2 обозначен каменистый безлесный участок. Он был частично покрыт кедровым лесом, полностью уничтоженным во время сильного пожара. Участок 3 — это пойма небольшой реки, заросшая кустарником. В качестве контрольного выбран участок 4 преимущественно соснового леса, не затронутого пожаром.
Для исследования выбраны данные японского РСА АЬОЭ-2 РАЬ8А11-2 (табл. 1), работающего в Ь-диапазоне (длина волны 23.5 см). В настоящее время это единственно работающий РСА в данном диапазоне. Всего на исследуемую территорию доступно только два полностью поляриметрических изображения (режим съемки НВС^, разрешение 6 м). Также были отобраны семь сцен на двойной поляризации (режим съемки
Таблица 1. Поляриметрические данные ALOS-2 PALSAR-2 для двух режимов съемки Table 1. ALOS-2 PALSAR-2 Polarimetrie data for two imaging modes
Рис. 1. Общий вид исследуемых участков Fig. 1. Set up view
Рис. 2. Вид тестового участка 1 Fig. 2. Figure of test site 1
Режим съемки
HBQ (полная поляризация) FBD (двойная поляризация)
ALOS2018861050-140928 ALOS2062331050-150719
Идентификатор сцены
ALOS2277611050-190714 ALOS2331431050-200712
FBD, разрешение 10 м) е 2014 по 2020 г. Отметим, что последние шесть цифр в идентификаторе сцены (см, табл. 1) означают дату съемки в формате ГГММДД,
Дополнительно использовались изображения спутникового радиолокатора с синтезированной апертурой Европейского космического агентства Sentinel-IB, работающего в ('-диапазоне (длина волны 5,6 см). Большими достоинствами изображений, полученных данным спутником, являются довольно малая временная база, составляющая для одного спутника 12 дней, и свободный доступ к данным, что позволяет формировать .мулы нвременные ряды. При использовании двух спутников (Sentinel-IA и Sentinel-IB) временной интервал между съемками составляет шесть дней, однако на территории Бурятии доступен только 12-дневный интервал. Для анализа использовались снимки, полученные в интерферометричееком режиме IW (interferometrie wide) и представленные в виде SLC-продуктов (Single Look Complex), Пространственное разрешение изображений составляло 15 м в поляризационных режимах VV и VII. т, е, излучение на вертикальной поляризации, а прием на основной вертикальной поляризации и на кросс-компоненте, Всего обработано 111 изображений, полученных в период с 25 февраля 2017 г, по 28 февраля 2021 г,
3. Результаты
Рассмотрим результаты поляриметрических разложений и анализа обратного рассеяния отдельно для каждого набора данных,
3.1. РСА ALOS-2 PALSAR-2
3.1.1. Полностью поляриметрические данные
В табл. 2 приведены результаты разложения Фрнмана-Дурдена (см, подразд, 1,1,2) для каждого участка, выбранного для исследования, В процентах представлены средние значения на выбранном участке. Первое значение в ячейке относится к изображению, полученному 7 мая 2017 г., второе — 25 августа 2019 г. Значения в табл. 2 показывают увеличение компоненты объемного рассеяния на участках 1-3, Максимальное увеличение составляет 20% для участков 1 (лесной подрост) и 2 (безлесный каменистый участок), В первом случае лесной подрост вырос на 0,5-1,0 м за два года с соответствующим увеличением общего количества элементарных рассенвателей в виде стволов и ветвей деревьев. Увеличение биомассы леса вызывает эффект "экранирования", приводящий к уменьшению компоненты поверхностного рассеяния на те же 20%, Рост объемной составляющей рассеяния на участке 3 в основном связан с увеличением плотности кустарников. Для незатронутого пожаром леса изменения компонент рассеяния не являются существенными.
Таблица 2. Результаты разложения Фримана-Дурдена, % Table 2. Results of Freeman - Durden decomposition, %
N- участка 1 2 3 4
Ps 49.1/29 52/33.4 66/49.8 7.2/8.8
Pd 5.9/5.3 5.4/4.5 5/7.9 7.6/6.2
Pv 44.9/65.7 42.6/62 29/42.2 85.1/85
3.1.2. Данные с двойной поляризацией
Из-за небольшого количества полностью поляриметрических данных (всего две сцены), к сожалению, невозможно проанализировать динамику лесовосстановления. Однако имеются семь изображений, сделанных в режиме ИЗБ (поляризация НН/НУ), которые можно использовать для анализа (см. табл. 1). Для двойных поляриметрических данных можно выполнить только разложение Клауда-Потье [13] (см. подразд. 1.1.1). В этом случае матрица когерентности [Т] имеет размер 2x2 и, соответственно, можно выделить только два механизма рассеяния. На рис. 3 показана динамика изменения энтропии Н (4) и угла а (5) за период с 2014 по 2020 г. Как и ранее, полученные значения усреднялись в пределах каждого участка, для того чтобы получить для него единственное значение анализируемого параметра. Из рис. 3 следует, что ни энтропия, ни параметр а существенно не меняются на протяжении семи лет и не могут служить индикаторами процессов лесовосстановления.
3.1.3. Временные ряды интенсивности обратного рассеяния
Поскольку поляриметрическое разложение данных с двойной поляризацией не выявило значительных изменений на тестовых участках, было решено проанализировать вариации интенсивности рассеянного сигнала. На рис. 4 показаны временные зависимости радиометрически откалиброванной и нормализованной амплитуды для различных состояний поляризации. Анализ графиков амплитуды показывает, что отличить лесной подрост от других участков практически невозможно. Так, например, корреляция между кривой для участка 1 (лесной подрост) и кривой 2 (безлесный каменистый участок) равна 0.92 для поляризации НН и 0.96 для НУ. Единственный участок, который позволяют идентифицировать амплитудные кривые, — это зрелый лес (участок 4). За счет рассеяния на крупных ветвях деревьев уровень сигнала на кросс-поляризации в среднем на 4.5 дБ выше, чем на других участках.
Рис. 3. Временные зависимости энтропии (а) и угла а (б) Fig. 3. Time dependencies of entropy (a) and a angle (6)
а б
HH HV
20.07.201Д 20.07.2015 19.07.2016 19.07.2017 19.07.2018 19.07.2019 18.07.2020 20.07.2014 20.07.2015 19.07.2016 19.07.2017 19.07.2018 19.07.2019 18.07.2020
• 1 - «- -9 4
Рис. 4. Временные зависимости радиометрически откалиброванной (а) и нормализованной (б) амплитуды
Fig. 4. Time dependencies of radiometrically calibrated (a) and normalized (6) amplitude
Таким образом, данные Ь-диапазона на двойной поляризации не пригодны для мониторинга процесса лесовосстановления. Это связано главным образом с высокой проникающей способностью волн Ь-диапазона. Невысокий подрастающий лес на временном интервале в 7 лет, даже несмотря на свою относительную густоту, не вносит заметного затухания в отраженный радарный сигнал, как и не приводит к видимым эффектам деполяризации.
В связи с невозможностью отследить динамику лесовосстановления в L-диапазоне с помощью данных с двойной поляризацией было решено использовать свободно доступные данные европейского PC A Sentinel-IB, работающего в С-диапазоне (длина волны 5.6 см). Периодичность радиолокационной съемки Sentinel-IB для территории Сибири составляет 12 дней. Регулярное обследование участков проводится с 2017 г. Радар Sentinel-1 излучает волны с вертикальной поляризацией (VV и VH) в отличие от ALOS-2 PALSAR-2, который в основном работает при горизонтальной поляризации излучения
3.2.1. Разложение Клауда-Потье
На рис. 5 показаны рассчитанные временные ряды энтропии и угла а для каждого участка. В данном случае поляриметрическое разложение Н/а позволяет отделить участки с растительностью от безлесного участка по величине энтропии и механизму рассеяния. Однако, как и ранее, сами участки с растительностью не разделяются между собой. Также следует отметить выраженную сезонную зависимость как энтропии, так и угла а.
3.2.2. Временные ряды интенсивности обратного рассеяния
Рассмотрим временные ряды радиометрически откалиброванной и нормализованной амплитуды обратного рассеяния. Вследствие малой длины волны обратное отражение в С-диапазоне очень чувствительно к изменениям подстилающей среды, например к изменению влажности почвы, процессам замерзания/оттаивания [18] и т.п. Все эти факторы сильно влияют на величину амплитуды, что особенно заметно в летний и осенне-весенний периоды (см. сезонную зависимость на рис. 5), и выражается в большой дисперсии результатов измерений. Поэтому для оценки процессов лесовосстановления было решено рассмотреть более стабильные зимние данные, когда согласно рис. 5, а энтропия минимальна. Всего использовано 38 зимних изображений, взятых в период с начала ноября одного года и по конец февраля года, следующего за ним.
3.2. РСА Sentinel-1В
(НН и НУ).
а
б
0.5
0.4 -1-1-1-1-1-
15.01.2017 12.09.2017 10.05.2018 05.01.2019 02.09.2019 29.04.2020 —•—1 - ♦ - 2 3 4
15
15.01.2017 12.09.2017 10.05.2018 05.01.2019 02.09.2019 29.04.2020
—•—1 - ♦ - 2 3
Рис. 5. Временные зависимости энтропии (а) и угла а (б) для данных Sentinel-IB Fig. 5. Time dependencies of entropy (a) and a angle (6) for Sentinel-IB data
04.02.2017 -9
i
a
w
05.01.2019
05.12.2020 04.02.2017 -1- -ig
.....mL.*....
........
•vv
......J* I
£-23 <
Л.Н
5
VH
05.01.2019
..........................................-Jr:.............................*..............
—--......... ^Jftt........... .......
......Ш- .......... '
•::::............................. fr:
Рис. 6. Временные зависимости радиометрически откалиброванной (а) и нормализованной (б) амплитуды в зимние периоды времени
Fig. 6. Time dependencies of radiometrically calibrated (a) and normalized (6) amplitude during winter
Полученные результаты представлены на рис. 6. Анализ кривых на рис. 6, а для согласованной вертикальной поляризации показывает, что наблюдается устойчивый рост отраженного сигнала для участков 1 и 3. Сравнение средней величины амплитуды за зимний период 2017-2018 гг. с аналогичным периодом 2020-2021 гг. показывает ее рост на 11% для участка 1 и на 12% для участка 3. Также следует отметить практически параллельный ход линий тренда для этих участков, что говорит о схожих процессах рассеяния. Отраженный сигнал возрастает и для участков 2 и 4, но в количественном выражении он существенно меньше (менее 2 % для участка 2 и 6 % для участка 4).
Зависимости для кросс-поляризации (рис. 6, б) показывают аналогичный тренд для участков 1 и 3, но слабее выражены количественно (менее 4% для обоих участков). Отдельно следует отметить участок 4, который демонстрирует стабильность средних отражательных характеристик на протяжении четырех лет (изменения -0.6%). Заметно меняется и поведение обратного радарного рассеяния на участке 2, где видно небольшое уменьшение амплитуды рассеяния. Весьма интересно сравнение участков 2 и 3 в пределах одного зимнего сезона. Как отчетливо видно для сезонов 2017/2018, 2018/2019 и 2020/2021, радиолокационный сигнал последовательно возрастает в течение зимы. Очевидно, что это связано с ростом снежного покрова (который может превышать 1 м в данной местности) и увеличением компоненты объемного рассеяния в общем отраженном сигнале. Однако, если участок 2 — каменистая безлесная территория, то участок 3 — заросли кустарника в пойме реки. Отсюда можно сделать вывод, что в зимний период, когда кустарник сбрасывает листву, а почва промерзает, рассеяние на участке с лиственным лесным подростом становится аналогичным рассеянию на безлесном участке, т. е. лиственный лесной подрост становится "радиопрозрачным". Хвойный лесной подрост на участке 1 не демонстрирует подобного поведения. Это может служить косвенным индикатором для разделения лесного подроста по типу растительности.
Таким образом, в С-диапазоне на основании анализа амплитудных временных рядов можно достоверно обнаружить лесной подрост как хвойных, так и лиственных пород и отличить его от взрослого леса и безлесных участков. По косвенным признакам лесной подрост можно разделить и по породному составу. Следует отметить, что все эти выводы справедливы именно для молодого подроста. Если сохранятся тренды, показанные на рис. 6, то через 5-6 лет, когда хвойный подрост достигнет высоты 2.5-3 м, возникнет эффект насыщения и подрост невозможно будет отличить от взрослого леса. Однако в таком случае возросшая биомасса леса должна оказывать влияние на отражение уже в более длинноволновом L-диапазоне.
3.3. Sentinel-2A/B
В качестве дополнительного шага проанализирован нормализованный разностный индекс растительности (NDVI). Расчет NDVI выполнен с помощью облачной платформы Google Earth Engine (GEE) [11, 12]. Earth Engine — это облачная платформа для геопространственного анализа данных в планетарных масштабах. Она позволяет использовать огромные вычислительные мощности компании Google для изучения самых разнообразных проблем: потерь лесов, засухи, стихийных бедствий, эпидемий, продовольственной безопасности, управления водными ресурсами, изменения климата и защиты окружающей среды. С технической точки зрения GEE — это многопетабайтный каталог данных, интегрированный с высокопроизводительным кластером серверов для параллельных вычислений. Доступ к системе и управление осуществляются через интерфейс прикладного программирования (API). Пользователь создает сценарии обработки данных в интерактивной среде разработки на языке программирования JavaScript и обеспечивает оперативное создание прототипов и визуализацию результатов "на лету".
Каталог Earth Engine содержит множество общедоступных наборов геопространственных данных, в том числе данных европейских мультиспектральных сенсоров Sentinel-2A/B с уровнем обработки L1C и L2A. В видимом (каналы В2 — 0.490 мкм, ВЗ — 0.560 мкм, В4 — 0.665 мкм) и ближнем инфракрасном диапазоне (канал В8 — 0.842 мкм) данные имеют пространственное разрешение 10 м. Индекс растительности, согласно обзору [19], для данных Sentinel-2 вычисляется с помощью выражения
NDVI =
В8-В4 В8 + В 4'
В общей сложности для расчета временных рядов (рис. 7) использовано бо-
лее ста безоблачных изображений, полученных в бесснежные периоды времени. Из рис. 7 следует, что не может идентифицировать лесной подрост, но может помочь
в определении его породного состава. Например, кустарники на участке 3 имеют более высокое значение ]^ГОУ1, чем молодняк на участке 1. Скалистая безлесная местность (участок 2) имеет менее 0.3 в летний период. Это означает, что увеличение объ-
емного рассеяния, как было показано поляриметрическим разложением (см. табл. 2), по существу не связано с растительностью. Индекс растительности на участке 1 практически совпадает со значениям на участке 4, где расположен сосновый лес. Таким образом, можно сделать вывод, что лесной подрост на первом участке содержит преимущественно молодые сосны.
Рис. 7. Временные зависимости NDVI Fig. 7. Time dependencies of NDVI
Заключение
Рассмотрена возможность выявления подроста лесных насаждений при помощи анализа временных рядов радиолокационных изображений ALOS-2 PALSAR-2 и Sentinel-1. Показано, что в L-диапазоне (ALOS-2 PALSAR-2) с помощью модельных разложений полностью поляриметрических данных можно уверенно идентифицировать рост лесных насаждений. Однако данные с двойной поляризацией уже неспособны отделить лесной подрост от других типов растительности. Данные, полученные в ('-д иапазоне (Sentinel-1) на согласованной вертикальной поляризации, наоборот, позволяют определить лесной подрост из анализа временных рядов амплитуды радарного рассеяния, взятых в зимний период времени. Привлечение данных оптических сенсоров, в частности Sentinel-2A/B, позволяет разделить выявленный подрост по породному составу и устранить ошибки классификации на безлесных участках, которые показывают рост объемного рассеяния на модельных поляриметрических разложениях.
Благодарности. Авторы признательны японскому аэрокосмическому агентству JAXA за данные ALOS-2 PALSAR-2, предоставленные в соответствии с проектом EO-RA3 (PI: ER3A2N129), Работа выполнена в рамках государственного задания Института физического материаловедения СО РАН,
Список литературы
fl] Bellassen V., Luyssaert S. Carbon sequestration: Managing forests in uncertain times. Nature. 2014; 506(7487):153-155. D01:10.1038/506153a.
[2] Jiang M., Medlyn B.E., Drake J.E., Duursma R.A., Anderson I.C., Barton C.V.M., Boer M.M., Carrillo Y., Castañeda-Gómez L., Collins L., Crous K.Y., De Kauwe M.G., dos Santos B.M., Emmerson K.M., Facey S.L., Gherlenda A.N., Gimeno Т.Е., Hasegawa S., Johnson S.N., Kánnaste A., Macdonald C.A., Mahmud K., Moore B.D., Nazaries L., Neilson E.H.J., Nielsen U.N., Niinemets U., Noh N.J., Ochoa-Hueso R., Pathare V.S., Pendall E., Pihlblad J., Piñeiro J., Powell J.R., Power S.A., Reich P.B., Renchon A.A., Riegler M., Rinnan R., Rymer P.D., Salomon R.L., Singh B.K., Smith В., Tjoelker M.G., Walker J.K.M., Wujeska-Klause A., Yang J., Zaehle S., Ellsworth D.S. The fate of carbon in a mature forest under carbon dioxide enrichment. Nature. 2020; 580(7802):227—231. DOI:10.1038/s41586-020-2128-9.
[3] Global biomass. Available at: https://www.eo4sd-forest.info/global-biomass (accessed January 17, 2022).
[4] FAO. 2020. Global forest resources assessment 2020. Rome. Available at: https://www.fao. org/forest-resources-assessment/fra-2020/maps/en (accessed January 17, 2022).
[5] Pinnington E.M., Casella E., Dance S.L., Lawless A.S., Morison J.I.L., Nichols N.K., Wilkinson M., Quaife T.L. Understanding the effect of disturbance from selective felling on the carbon dynamics of a managed woodland by combining observations with model predictions. Journal of Geophysical Research: Biogeosciences. 2017; 122(4) :886-902. D01:10.1002/2017JG003760.
[61 Pugh T.A.M., Lindeskog M., Smith В., Poulter В., Arneth A., Haverd V., Calle L. Role of forest regrowth in global carbon sink dynamics. Proceedings of the National Academy of Sciences. 2019; 116(10):4382-4387. DOLIO. 1073/pnas. 1810512116.
[7] Chazdon R.L., Broadbent E.N., Rozendaal D.M.A., Bongers F., Zambrano A.M.A., Aide T.M., Balvanera P., Becknell J.M., Boukili V., Brancalion P.H.S.,
Craven D., Almeida-Cortez J.S., Cabral G.A.L., Jong Ben de, Denslow J.S., Dent D.H., DeWalt S.J., Dupuy J.M., Durán S.M., Espirito-Santo M.M., Fandino M.C., César R.G., Hall J.S., Hernández-Stefanoni J.L., Jakovac C.C., Junqueira А.В., Kennard D., Letcher S.G., Lohbeck M., Martinez-Ramos M., Massoca P., Meave J.A., Mesquita R., Mora F., Muñoz R., Muscarella R., Nunes Y.R.F., Ochoa-Gaona S., Orihuela-Belmonte E., Peña-Claros M., Pérez-García E.A., Piotto D., Powers J.S., Rodriguez-Velazquez J., Romero-Pérez I.E., Ruiz J., Saldarriaga J.G., Sanchez-Azofeifa A., Schwartz N.B., Steininger M.K., Swenson N.G., Uriarte M., Breugel M., Wal H., Veloso M.D.M., Vester H., Vieira I.C.G., Bentos T.V., Williamson G.B., Poorter L. Carbon sequestration potential of second-growth forest regeneration in the Latin American tropics. Science Advances. 2016. DOI: 10.1126/sciadv. 1501639.
[81 Lehmann E.A., Caccetta P., Lowell K., Mitchell A., Zhou Z.-S., Held A., Milne Т., Tapley I. SAR and optical remote sensing: Assessment of complementarity and interoperability in the context of a large-scale operational forest monitoring system. Remote Sensing of Environment. 2015; 156:335-348. D01:10.1016/j.rse.2014.09.034.
[9] Чимитдоржиев Т.Н., Дмитриев А.В., Кирбижекова И.И., Шерхоева А.А., Бал-тухаев А.К., Дагуров П.Н. Дистанционные оптико-микроволновые измерения параметров леса: современное состояние исследований и экспериментальная оценка возможностей. Современные проблемы дистанционного зондирования Земли из космоса. 2018. 15(4):9-24. DOI:10.21046/2070-7401-2018-15-4-9-24.
[10] Чимитдоржиев Т.Н., Дмитриев А.В., Дагуров П.Н. Технология совместного анализа временных рядов изображений интерферометрической когерентности Sentinel-1 и вегетационного индекса по данным Sentinel-2 для мониторинга сельскохозяйственных полей. Современные проблемы дистанционного зондирования Земли из космоса. 2020; 17(4):61-72. DOLIO.21046/2070-7401-2020-17-4-61-72.
[11] Gorelick N., Hancher M., Dixon M., Ilyushchenko S., Thau D., Moore R. Google Earth Engine: Planetarv-scale geospatial analysis for everyone. Remote Sensing of Environment. 2017; (202):18-27. D01:10.1016/j.rse.2017.06.031.
[12] Google Earth Engine. Available at: https://earthengine.google.com (accessed January 17, 2022).
[13] Cloude S. Polarisation: Applications in remote sensing. OUP Oxford; 2009: 472.
[14] Lee J.-S., Pottier E. Polarimetric radar imaging: Prom basics to applications. CRC Press: Boca Raton; 2009: 422.
[15] Freeman A., Durden S.L. A three-component scattering model for polarimetric SAR data. IEEE Transactions on Geoscience and Remote Sensing. 1998; 36(3):963-973. DOLIO.1109/36.673687.
[16] Small D. Flattening gamma: Radiometric terrain correction for SAR imagery. IEEE Transactions on Geoscience and Remote Sensing. 2011; 49(8):3081-3093. DOI: 10.1109/TGRS.2011.2120616.
[17] SNAP — ESA Sentinel Application Platform v8.0. Available at: http://step.esa.int (accessed January 17, 2022).
[18] Dagurov P.N., Chimitdorzhiev T.N., Dmitriev A.V., Dobrynin S.I. Estimation of snow water equivalent from L-band radar interferometrv: simulation and experiment. International Journal of Remote Sensing. 2020; 41(24) :9328-9359. DOLIO.1080/01431161.2020.1798551.
[19] Sentinel-2 Level-2A algorithm overview. Available at: https://sentinel.esa.int/web/ sentinel/technical-guides/sentinel-2-msi/level-2a/algorithm (accessed January 17, 2022)._
Вычислительные технологии, 2022, том 27, № 2, с. 105-121. © ФИЦ ИВТ, 2022 ISSN 1560-7534
Computational Technologies, 2022, vol. 27, no. 2, pp. 105-121. © FRC ICT, 2022 elSSN 2313-691X
INFORMATION TECHNOLOGIES
DOI: 10.25743/ICT.2022.27.2.009 Optics and microwave detection of forest restoration after fires
Dmitriev Aleksey V.*, Chimitdorzhiev Tumen N., Dagurov Pavel N.
Institute of Physical Materials Science SB RAS, 670047, Ulan-Ude, Russia
* Corresponding author: Dmitriev Aleksev V., e-mail: davQipms .bscnet .ru Received January 21, 2022, accepted February 14, 2022.
Abstract
The problem of large-scale assessment of forest restoration after artificial deforestation and wildfires is relevant in connection with climate change and the corresponding desire of the world community for further low-carbon development. One of the promising methods aimed at solving this problem is remote sensing of the Earth using space-based synthetic aperture radars (SAR). The paper proposes a comprehensive approach to assessing the dynamics of forest plantings growth using time series analysis of ALOS-2 PALSAR-2 and Sentinel-1 radars data, as well as Sentinel-2A/B optical sensors data. It is shown that with the help of model-based decomposition (Freeman-Durden decomposition) of fully polarimetric data in L-band (ALOS-2 PALSAR-2), the increasing of volume scattering component and the corresponding decrease in the surface component can confidently identify the growth of young forest. However, the data of ALOS-2 PALSAR-2 with dual polarization are not able to separate forest undergrowth from other types of vegetation over a seven-year observation interval. This is also true for the C-band. Thus, the polarimetric Cloud-Pottier decomposition of Sentinel-1 data allowing only separation for the areas with vegetation from the treeless ones. Time series analysis of radiometricallv corrected radar backscattering at vertical co-polarization in this band, imaged in the winter period of time, allows reliable determining of the dynamics for the growth of forest plantations. The use of freely available Sentinel-2A/B multispectral sensors data makes it possible to further divide the identified undergrowth by species composition and exclude classification errors of radar data in treeless areas, which show an increasing of volume backscattering component on model-based polarimetric decompositions.
Keywords: reforestation, polarimetric decompositions, time series analysis, ALOS-2 PALSAR-2, Sentinel-1/2.
Citation: Dmitriev A.V., Chimitdorzhiev T.N., Dagurov P.N. Optics and microwave detection of forest restoration after fires. Computational Technologies. 2022; 27(2):105—121. DOI: 10.25743/ICT.2022.27.2.009. (In Russ.)
Acknowledgements. The authors are grateful to the Japanese Aerospace Agency (JAXA) for the ALOS-2 PULSAR-2 data provided in accordance with the EO-RA3 project (PI: ER3A2N129). The research was carried out within the framework of the state task of Institute of Physical Materials Science SB RAS.
References
1. Bellassen V., Luyssaert S. Carbon sequestration: Managing forests in uncertain times. Nature. 2014; 506(7487): 153-155. D01:10.1038/506153a.
2. Jiang M., Medlyn B.E., Drake J.E., Duursma R.A., Anderson I.C., Barton C.V.M., Boer M.M., Carrillo Y., Castañeda-Gómez L., Collins L., Crous K.Y., De Kauwe M.G., dos Santos B.M., Emmerson K.M., Facey S.L., Gherlenda A.N., Gimeno T.E.,
120
A. B. , [Miirpiici;. T. H. 1 IIIMIII: H>I»KIICI;. II. H. , buy¡nm
Hasegawa S., Johnson S.N., Kannaste A., Macdonald C.A., Mahmud K., Moore B.D., Nazaries L., Neilson E.H.J., Nielsen U.N., Niinemets Ü., Noh N.J., Ochoa-Hueso R., Pathare V.S., Pendall E., Pihlblad J., Piñeiro J., Powell J.R., Power S.A., Reich P.B., Renchon A.A., Riegler M., Rinnan R., Rymer P.D., Salomón R.L., Singh B.K., Smith B., Tjoelker M.G., Walker J.K.M., Wujeska-Klause A., Yang J., Zaehle S., Ellsworth D.S. The fate of carbon in a mature forest under carbon dioxide enrichment. Nature. 2020; 580(7802):227-231. DC)I:10.1038/s41586-020-2128-9.
3. Global biomass. Available at: https://www. eo4sd-forest. info/global-biomass (accessed January 17, 2022).
4. FAO. 2020. Global forest resources assessment 2020. Rome. Available at: https://www.fao.org/ forest-resources-assessment/fra-2020/maps/en (accessed January 17, 2022).
5. Pinnington E.M., Casella E., Dance S.L., Lawless A.S., Morison JJ.L., Nichols N.K., Wilkinson M., Quaife T.L. Understanding the effect of disturbance from selective felling on the carbon dynamics of a managed woodland by combining observations with model predictions. Journal of Geophysical Research: Biogeosciences. 2017; 122(4):886-902. D01:10.1002/2017JG003760.
6. Pugh T.A.M., Lindeskog M., Smith B., Poulter B., Arneth A., Haverd V., Calle L. Role of forest regrowth in global carbon sink dynamics. Proceedings of the National Academy of Sciences. 2019; 116(10):4382-4387. D01:10.1073/pnas.l810512116.
7. Chazdon R.L., Broadbent E.N., Rozendaal D.M.A., Bongers F., Zambrano A.M.A., Aide T.M., Balvanera P., Becknell J.M., Boukili V., Brancalion P.H.S., Craven D., Almeida-Cortez J.S., Cabral G.A.L., Jong Ben de, Denslow J.S., Dent D.H., DeWalt S.J., Dupuy J.M., Durán S.M., Espirito-Santo M.M., Fandino M.C., César R.G., Hall J.S., Hernández-Stefanoni J.L., Jakovac C.C., Junqueira A.B., Kennard D., Letcher S.G., Lohbeck M., Martínez-Ramos M., Massoca P., Meave J.A., Mosquita R., Mora F., Muñoz R., Muscarella R., Nunes Y.R.F., Ochoa-Gaona S., Orihuela-Belmonte E., Peña-Claros M., Pérez-García E.A., Piotto D., Powers J.S., Rodríguez-Velazquez J., Romero-Pérez I.E., Ruiz J., Saldarriaga J.G., Sanchez-Azofeifa A., Schwartz N.B., Steininger M.K., Swenson N.G., Uriarte M., Breugel M., Wal H., Veloso M.D.M., Vester IT., Vieira I.C.G., Bentos T.V., Williamson G.B., Poorter L. Carbon sequestration potential of second-growth forest regeneration in the Latin American tropics. Science Advances. 2016. DOI: 10.1126/sciadv. 1501639.
8. Lehmann E.A., Caccetta P., Lowell K., Mitchell A., Zhou Z.-S., Held A., Milne T., Tapley I. SAR and optical remote sensing: Assessment of complementarity and interoperability in the context of a large-scale operational forest monitoring system. Remote Sensing of Environment. 2015; 156:335-348. D01:10.1016/j.rse.2014.09.034.
9. Chimitdorzhiev T.N., Dmitriev A.V., Kirbizhekova I.I., Sherhoeva A.A., Baltukhaev A.K., Dagurov P.N. Remote optical-microwave measurements of forest parameters: Modern state of research and experimental assessment of potentials. Sovremennye Problemy Distantsionnogo Zondirovaniya Zemli iz Kosmosa. 2018; 15(I):9 21. DOI: 10.21046/2070-7401-2018-15-4-9-24.
10. Chimitdorzhiev T.N., Dmitriev A.V., Dagurov P.N. Technology of joint analysis of Sentinel-1 interferometric coherence time series and vegetation index based on Sentinel-2 data for monitoring agricultural fields. Sovremennye Problemy Distantsionnogo Zondirovaniya Zemli iz Kosmosa. 2020; 17(4):61-72. DOI: 10.21046/2070-7401-2020-17-4-61-72.
11. Gorelick N., Hancher M., Dixon M., Ilyushchenko S., Thau D., Moore R. Google Earth Engine: Planetary-scale geospatial analysis for everyone. Remote Sensing of Environment. 2017; (202):18-27. D01:10.1016/j.rse.2017.06.031.
12. Google Earth Engine. Available at: https://earthengine.google.com (accessed January 17, 2022).
13. Cloude S. Polarisation: Applications in remote sensing. OUP Oxford; 2009: 472.
14. Lee J.-S., Pottier E. Polarimetric radar imaging: From basics to applications. CRC Press: Boca Raton; 2009: 422.
15. Freeman A., Durden S.L. A three-component scattering model for polarimetric SAR data. IEEE Transactions on Geoscience and Remote Sensing. 1998; 36(3):963-973. DOLIO.1109/36.673687.
16. Small D. Flattening gamma: Radiometric terrain correction for SAR imagery. IEEE Transactions on Geoscience and Remote Sensing. 2011; 49(8):3081-3093. DOL10.1109/TGRS.2011.2120616.
17. SNAP — ESA Sentinel Application Platform v8.0. Available at: http://step.esa.int (accessed January 17, 2022).
18. Dagurov P.N., Chimitdorzhiev T.N., Dmitriev A.V., Dobrynin S.I. Estimation of snow water equivalent from L-band radar interferometry: simulation and experiment. International Journal of Remote Sensing. 2020; 41(24):9328-9359. D01:10.1080/01431161.2020.1798551.
19. Sentinel-2 Level-2A algorithm overview. Available at: https://sentinel.esa.int/web/sentinel/ technical-guides/sentinel-2-msi/level-2a/algorithm (accessed January 17, 2022).
