Application of the Stacking-InSAR method for analyzing changes in forest canopy height

Full Text

ВВЕДЕНИЕ

Одним из эффективных методов мониторинга лесной среды является дифференциальная радиолокационная интерферометрия (ДРИ), которая позволяет оценивать состояние и изменения высоты/биомассы леса. В связи с климатическими изменениями и необходимостью проведения анализа объемов секвестрации углерода, такие оценки приобретают всё большую актуальность. Так, в работе (Solberg et. al., 2014) в рамках методологии REDD (Reducing Emissions from Deforestation and forest Degradation in developing countries ‒ сокращение выбросов в результате обезлесения и деградации) показано, что увеличение биомассы на 14.9 т/га соответствует увеличению на метр средней высоты леса. Такая высота измерялась с использованием метода ДРИ по данным бистатической пары радиолокаторов X-диапазона TerraSAR-X/TanDEM-X. Кроме того, в работе (Solberg et al., 2014) было предложено проводить оценки изменений высоты леса на основе сравнения глобальных цифровых моделей местности (ЦММ), полученных по данным радиолокационной интерферометрии C-диапазона в 2000 г. (SRTM) и X-диапазона TerraSAR-X/TanDEM-X. С учетом различий проникающей способности длин радиоволн X- и C-диапазонов и выявленных различных ошибок при построении ЦММ было заявлено о 50% точности при выявлении изменений высоты леса.

Результаты радиолокационного интерферометрического мониторинга высоты зрелых бореальных лесов Сибири также по данным бистатической радиолокационной системы X-диапазона TerraSAR-X/TanDEM-X представлены в работах (Zakharov et al., 2016; Чимитдоржиев и др., 2018). В этих работах предполагалось, что фазовые центры рассеяния электромагнитных волн X-диапазона в силу низкой проникающей способности располагались в верхнем ярусе сплошного лесного полога. Лесной полог ‒ совокупность крон сомкнувшихся деревьев, размещающихся в одном (для лесного подроста примерно равного возраста) или нескольких (для зрелого леса) ярусах. Соответственно, изменение высоты лесного полога характеризует динамику средней высоты леса. Подобное представление позволило предложить сравнительно простую методику интерферометрических площадных измерений высоты зрелого соснового леса на равнине (Чимитдоржиев и др., 2018). Суть методики заключалась в оценке разности высот расположения фазовых центров рассеяния для радиоволн X- и L-диапазонов длин волн. Было выявлено систематическое занижение действительной высоты зрелого леса: расхождение результатов радарной интерферометрии и подспутниковых измерений достигало 5.5 м при значениях полноты древостоя 0.5, 0.9 и 1, а при полноте древостоя от 0.6 до 0.8 это расхождение изменялось в диапазоне от 2 до 4 м. Занижение средней высоты деревьев, обнаруженное радарным интерферометрическим методом, при полноте леса 0.5 было интерпретировано как частичное рассеяние сигнала Х-диапазона от безлесного пространства между деревьями. Причиной занижения средней высоты деревьев при полноте 0.9 и 1 стала недостаточная проникающая способность волн L-диапазона под лесной полог зрелого леса с биомассой более чем 100 т/га (Чимитдоржиев и др., 2018). Методика ДРИ, помимо оценки состояния растительности в целом, позволяет также оценивать динамику зрелого и молодого леса. Для этого необходимы радиолокационные интерферометрические данные X- и L-диапазонов, а в настоящее время в свободном доступе имеются только данные радиолокаторов Sentinel-1A/B, функционирующих в C-диапазоне радиоволн. В связи с этим, основываясь на вышеизложенных результатах, для мониторинга изменений высоты лесного полога предлагается использовать подход, основанный на анализе долговременных временных рядов радарных интерферометрических данных C-диапазона, которые могут быть получены с помощью онлайн платформ обработки спутниковых изображений. Это позволит минимизировать возможные ошибки единичных измерений и оценить динамику растительности, например, при помощи метода Stacking-InSAR (Бондур и др., 2023; Bondur et al., 2024).

ОПИСАНИЕ ТЕСТОВОГО ПОЛИГОНА И ИСПОЛЬЗУЕМЫХ ДАННЫХ

В качестве тестового использовался полигон (см. рис. 1), который содержал молодой лес (подрост сосны одного возраста и одинаковой высоты), зрелый сосновый лес без изменений и лес с рубками прореживания, которые были произведены за рассматриваемый период времени. Подобная представительность лесного покрова позволяла получить относительные оценки для основных состояний хвойного леса. Его подробное описание приведено в работе (Дмитриев и др., 2022). Часть тестового полигона относилась к землям сельскохозяйственного назначения, однако в последние десятилетия 20 века данная местность не обрабатывалась. С начала 2000-х годов начался интенсивный рост соснового молодняка, состояние которого на лето 2007 г. показано на рис. 1. Красному цвету на рис. 1 соответствует спектральное отражение в 4 канале (NIR) мультиспектральной камеры AVNIR-2 спутника ALOS-1 (пространственное разрешение 10 м), зеленому цвету – спектральное отражение во 2 канале, а синему цвету – в 1 канале. Соответственно, подрост хвойных пород (молодой лес) выделяется более ярким красным цветом на фоне темно-красного зрелого соснового леса.

 

Рис. 1. Расположение тестового полигона и состояние лесного покрова на изображении с мультиспектральной камеры ALOS-1 AVNIR-2, полученном 16.07.2007 г. (каналы 4-2-1).

 

Как было показано ранее в работе (Дмитриев и др., 2022), залесение тестового полигона происходило неравномерно, что визуально можно наблюдать и на рис. 1. Первой начала интенсивно зарастать восточная часть полигона – контур розового цвета. Возраст этого молодого леса составлял 18–20 лет. В качестве контрольного участка для сравнения использовался фрагмент зрелого плотного леса (контур зелёного цвета на рис. 1), который оставался неизменным на всем протяжении периода наблюдений. Между этими участками расположен фрагмент леса с рубками прореживания, которые были выполнены в 2017–2018 гг. Временной ряд используемых радиолокационных данных охватывал этот период рубок. Для наглядности, далее по тексту количественные результаты представлялись на основе анализа профиля (см. рис. 1), проходящего от точки А до точки D. Соответственно зрелый плотный лес располагался на профиле от точки А до точки В, а лесной подрост на профиле от точки C до точки D. Между ними от точки В до точки С располагался лес с рубками прореживания.

Для обработки больших объемов радиолокационных данных, полученных со спутников Sentinel-1A/B, был использован функционал онлайн платформы Alaska Satellite Facility’s Hybrid Pluggable Processing Pipeline (HyP3) (Hogenson et al., 2020). Эта онлайн платформа предоставляет доступ к радиолокационным данным спутников Sentinel-1A/B, позволяет выбирать данные для радиолокационной интерферометрии и выполнять интерферометрическую обработку, начиная от совмещения до получения изображения развернутой фазы с минимальным пространственным разрешением 40 м.

Для расчетов были использованы все доступные данные (независимо от сезона) с 2014 по 2021 г.: 344 интерферограммы радиолокатора Sentinel-1B с временной базой 12, 24 и 36 дней на нисходящей относительной орбите 135 и согласованной вертикальной поляризации. Распределение интерферограмм по годам приведено в табл. 1.

 

Таблица 1. Количество развернутых интерферограмм

Год	2014	2015	2016	2017	2018	2019	2020	2021	Итого
Количество интерферограмм	15	15	30	79	83	42	29	51	344

 

ОПИСАНИЕ МЕТОДА РАСЧЕТА

Для расчета скорости смещений высоты фазовых центров рассеяния, фактически характеризующих динамику высоты лесного полога, вдоль линии обзора радиолокатора использовался метод Stacking-InSAR. Суть этого метода сводится к взвешенному суммированию всех N развернутых интерферометрических фаз φ_j с весами, пропорциональными их временной базе . Средняя скорость изменения фазы, описывающей смещения, записалась в виде (Zhang et al., 2021)

$\overline{\phi }=\sum _{j=1}^N{\phi }_j∆T_j/\sum _{j=1}^N∆T_j^2$.

Тогда средняя скорость деформаций вдоль линии обзора радиолокатора (англ. Line of Sight - LOS) выражается формулой

${\nu }_{деф}=-\lambda \overline{\phi }/4\mathrm{\pi }$,

где λ – длина волны радиолокатора.

С учетом пространственного разрешения полученных интерферометрических данных и самого метода суммирования взвешенных интерферометрических фаз были сделаны следующие предположения/допущения:

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

Ранее авторами в работе (Bondur et al., 2024) тропосферная задержка была скорректирована при помощи данных GACOS (Generic Atmospheric Correction Online Service for InSAR) (Yu et al., 2020), где использовались результаты модельных расчётов, в том числе, полученные с использованием многочисленных GPS-данных в районе исследования. В Байкальском регионе сеть GPS представлена единичными приборами, поэтому результаты, полученные с коррекцией по данным GACOS, следует сопоставить с результатами без такой коррекции.

РЕЗУЛЬТАТЫ ОБРАБОТКИ И МЕТОДИКА АНАЛИЗА

Первоначально по изображениям NDVI (с пороговым значением более 0.2), полученным по спутниковым изображениям Sentinel-2 для зимнего периода времени 2020 г., была сформирована маска вечнозелёного хвойного леса. Это позволило выделить на изображениях только объект исследования – сосновый лес, включая сосновый подрост различного возраста.

Далее по всему временному ряду данных интерферометрической фазы (см. табл. 1), с использованием метода Stacking-InSAR было получено поле скоростей изменения высоты фазовых центров рассеяния, которое фактически характеризовало динамику изменения высоты лесного полога за счёт линейного и радиального прироста стволов и ветвей деревьев. Результаты, полученные с коррекцией тропосферной задержки по данным GACOS (Yu et al., 2020), показали, что скорость прироста на участке лесного подроста составляла в среднем 1.6 см/год, а для участка зрелого леса 0.9 см/год. В то время, как без атмосферной коррекции скорость прироста составляла соответственно –0.3 см/год и –1.2 см/год. Это не согласовывалось с результатами ранее выполненных наземных полевых наблюдений. Таким образом, результаты без коррекции тропосферной задержки в данном случае являлись ошибочными. По этой причине и с учетом рекомендаций работы (Волкова и др., 2024) при расчетах с использованием развернутой интерферометрической фазы были использованы атмосферные поправки.

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

Для учёта перечисленных эффектов необходимо выбрать временные ряды данных только для бесснежного периода и при этом также исключать данные в сезонные периоды перехода от положительных температур к отрицательным и обратно. Для нивелирования временной декорреляции, вследствие изменения влажности или неучтенного влияния атмосферы в летний период времени использовались данные со средней когерентностью не менее 0.5 по всему тестовому полигону. В итоге, всем вышеперечисленным критериям с учетом ограничений по величине временной базовой линии (не более 36 дней) соответствовали 16 интерферограмм 2017 г., 17 интерферограмм 2018 г. и 15 за 2019 гг.

На рис. 2 представлены результаты, полученные методом Stacking-InSAR по этим временным рядам, а на рис. 3 соответствующие профили вдоль линии ABCD. Для удобства профили разделены на отрезки (см. рис. 1) с соответствующими линейными трендами для зрелого леса (АВ), леса с рубками прореживания (ВС) и молодого леса (CD). Представленные результаты показали, что максимальный прирост всех участков соснового леса наблюдался в вегетационный период 2018 г. и составлял величину порядка 20 см/год. При этом для зрелого леса без рубок разброс данных, полученных в 2018 г., варьировался в пределах 15÷22 см/год, а для леса с рубками отмечался разброс значений прироста от 10 до 22 см/год. Можно предположить, что смещение нижней границы значений прироста с 15 до 10 см/год было связано с рубками леса.

 

Рис. 2. Поле скоростей прироста сосны за вегетационный период: а – 2017 г., б – 2018 г., в – 2019 г.

 

Рис. 3. Скорости прироста сосны вдоль профиля ABCD за вегетационный период в разные годы.

 

В 2017 г. и 2019 г. для первых двух отрезков, характеризующих фрагменты тестового полигона для зрелого леса без рубки (см. участок AB на рис. 3) и с рубкой прореживания (см. участок BC на рис. 3), наблюдались сравнимые величины ‒ порядка 5±3 см/год.

Из литературы известно, что величина линейного и радиального прироста деревьев связана с количеством осадков. Так, например, в работе (Васильев и др., 2020) вычисление корреляции между временными рядами прироста древесины и количеством осадков проводилось, начиная с мая месяца года, предшествующего приросту, и заканчивалось в сентябре текущего года. При наземных исследованиях у образцов деревьев измерялись ширина годичных колец и ширина поздней древесины. Максимальная корреляционная связь была выявлена с суммой осадков для периода времени с мая по июль (r = 0.78). В данном случае количество осадков для тестового полигона в анализируемый период времени, полученным по данным (ERA5-Land Daily Aggregated, 2024), распределилось следующим образом: в 2016 г. – 3.99 см; в 2017 г. – 4.58 см; в 2018 г. – 5.62 см; в 2019 г. – 4.67 см. Различие между количеством осадков в 2018 г. и в 2017 г. составляло 1.04 см, а в 2018 по сравнению с 2019 г. – 0.95 см. Таким образом, в обоих случаях (2018 г. по отношению к 2017 и 2019 гг.) это различие составляло величину порядка 20% по сравнению с количеством осадков в 2017 г. Это вероятно и послужило причиной большего прироста лесной растительности в 2018 г.

Отсутствие существенных различий скорости прироста сосны, выявленных в 2017 и в 2019 годах для участков AB и BC, иллюстрирующихся на графике, представленном рис. 3, может быть связано с двумя факторами:

1. общий объем надземной биомассы леса после рубки прореживания остался в пределах известного порога насыщения 20 т/га для радиолокационных данных С-диапазона, выше которого оценка различий по биомассе затруднительна (Imhoff, 1995);
2. после прореживания увеличился прирост деревьев из-за увеличения пространства/площади.

Также не исключается совместное влияние этих двух причин. Это согласуется с данными, полученными в работе (Иванов и др., 2018).

Для участка CD (молодой лес) тренды скоростей прироста растительности в 2017 и в 2019 гг. заметно расходятся. В 2017 г. скорость поднятия лесного полога составляла 7 см/год по всему отрезку трассы на исследуемом полигоне. Далее интенсивный прирост растительности происходил в 2018 г., после которого в вегетационный период 2019 г. прирост существенно замедлился (2–3 см/год в северной части участка CD) вплоть до нулевых величин скорости в южной половине участка CD. Полевые наблюдения показали, что в южной части этого участка молодого леса (возрастом 18–20 лет) плотность растительности была существенно выше (0.36 дерева/м² при среднем диаметре стволов 11.4 см). При этом было отмечено частичное перекрытие крон деревьев. Такой эффект для участка CD (см. рис. 3) с учетом результатов работы (Иванов и др., 2018) можно объяснить тем, что молодой сосновый лес (с высокой плотностью) после интенсивного прироста в 2018 г. достиг определенного предела вследствие избыточной густоты – отсутствие свободного пространства для радиального прироста и соответственно линейного.

ЗАКЛЮЧЕНИЕ

Таким образом, показана принципиальная возможность радиолокационного интерферометрического мониторинга средней скорости изменения высоты фазовых центров рассеяния лесным пологом по временным рядам изображений интерферометрической фазы, полученной в результате использования подхода, основанного на использовании облачных расчетов. Полученные результаты согласуются с полевыми наблюдениями и могут быть использованы для расчета динамики надземной биомассы на основе зависимостей, полученных в работе (Solberg et al., 2014).

ИСТОЧНИК ФИНАНСИРОВАНИЯ

Исследования выполнены при поддержке гранта РНФ № 24-27-20044 (https://rscf.ru/project/24-27-20044/).

About the authors

V. G. Bondur

AEROCOSMOS Research Institute for Aerospace Monitoring

Author for correspondence.
Email: vgbondur@aerocosmos.info
Russian Federation, Moscow

T. N. Chimitdorzhiev

Institute of Physical Materials Science SB RAS

Email: vgbondur@aerocosmos.info
Russian Federation, Ulan-Ude

A. V. Dmitriev

Institute of Physical Materials Science SB RAS

Email: vgbondur@aerocosmos.info
Russian Federation, Ulan-Ude

Zh. D. Nomshiev

Institute of Physical Materials Science SB RAS

Email: vgbondur@aerocosmos.info
Russian Federation, Ulan-Ude

References

Supplementary files