Исследование Земли из Космоса, 2023, № 5, стр. 3-12

ВВЕДЕНИЕ

Одним из эффективных инструментов дистанционных измерений биофизических параметров леса являются поляриметрические радиолокаторы с синтезированной апертурой (РСА) (Cloude, 2009; Lee, Pottier, 2009; Moreira et al., 2013; Richards, 2009), которые позволяют получить достаточно полную информацию о состоянии и изменениях лесной растительности (Bondur et al., 2022; Dobson et al., 1992; Koyama et al., 2022, 2019; Le Toan et al., 1992; Yu, Saatchi, 2016; Бондур и др., 2019; Дмитриев и др., 2022; Чимитдоржиев и др., 2018), в т.ч. для бореальных лесов Сибири (Bondur et al., 2022; Бондур и др., 2019; Чимитдоржиев и др., 2018). Для оценки леса более эффективно использовать радиолокационные данные L-диапазона, поскольку в большинстве случаев радиолокационный эхо-сигнал проникает под лесной полог и отражается обратно от подстилающей поверхности (Bondur et al., 2022; Koyama et al., 2022, 2019; Бондур и др., 2019; Чимитдоржиев и др., 2018). При этом происходит многократное рассеяние радарного эхо-сигнала от стволов и ветвей деревьев. В случае лесного подроста, временное увеличение обратного радарного рассеяния в направление на радар показывает увеличение размеров и густоты молодого леса (Bondur et al., 2022; Дмитриев и др., 2022). Однако достаточно быстро возникает эффект “насыщения”, выражающийся в том, что обратное рассеяние от лесной среды возрастает до некоторого предела, до определенной величины надземной биомассы леса, и в дальнейшем перестает расти. Например, для хвойного леса в P-диапазоне это 200 тонн/га (Dobson et al., 1992), в L-диапазоне 100 т/га (Dobson et al., 1992; Le Toan et al., 1992; Yu and Saatchi, 2016). Тем не менее рост деревьев (увеличение биомассы) может продолжаться, хотя уже и не такой интенсивный. Так, ранее авторами в (Дмитриев и др., 2022), на примере соснового подроста, было показано устойчивое возрастание величины обратного радиолокационного рассеяния неоднородностями леса на согласованной горизонтальной поляризации в L-диапазоне длин волн и последующее “насыщение”. В тоже время соответствующий тренд для кросс-поляризации не достиг уровня контрольного участка зрелого леса, что объясняется различием в механизмах рассеяния радарного сигнала на неоднородностях лесного полога (стволы и ветви) на различных поляризациях. Этот результат, подтверждающий классические зависимости величины эффективной поверхности рассеяния (ЭПР) от биомассы (Dobson et al., 1992; Le Toan et al., 1992; Yu and Saatchi, 2016), показывает необходимость совместного анализа поляриметрических данных.

К настоящему времени разработано множество методов и алгоритмов анализа полностью поляриметрических радарных изображений (Cloude and Pottier, 1997; Freeman, 2007; Freeman and Durden, 1998; Krogager et al., 1997; Touzi, 2007; Yamaguchi et al., 2005). Эти методы показали высокую эффективность при количественных измерениях биофизических параметров леса. К сожалению, полностью поляриметрический режим радиолокационного зондирования используется достаточно редко. Так, в случае космическим радиолокатором ALOS-2 PALSAR-2 (L-диапазон) глобальная съемка в полностью поляриметрическом режиме проводится один раз в 5 лет, что в лучшем случае дает только два изображения на текущий момент времени. Более распространен режим съемки с двойной поляризацией, как правило на согласованной горизонтальной (HH) или вертикальной (VV), и одним из соответствующих вариантов кросс-поляризации (HV или VH). Такой режим радиолокационного зондирования экономичнее с точки зрения объема данных и площади покрытия (полоса обзора шире). Однако существенно снижается информативность и количество поляриметрических декомпозиций, пригодных для анализа. В известной нам литературе, для радарных изображений с двойной поляризацией предложена только декомпозиция Клауд–Потье (Cloude, 2007), в литературе и программных продуктах именуемая как декомпозиция H-α. Однако, как было показано в работе (Ji, Wu, 2015) классификация природных объектов по физическим механизмам рассеяния малоэффективна для режимов съемки с двойной поляризацией (HH-HV или VV-VH).

Тем не менее, можно предположить, что с ростом лесной растительности параметры H и α разложения Клауда–Потье будут также последовательно расти вдоль кривой, ограничивающей область их допустимых значений. Это позволит обнаружить и оценить динамику роста лесного молодняка с помощью данных с двойной поляризацией. Данный подход был использован в работе (Guo et al., 2018), в которой декомпозиция Клауда–Потье данных РСА Sentinel-1 (C-диапазон, длина волны 5.6 см) с двойной поляризацией VV-VH была использована с целью извлечения дополнительного параметра, позволившего провести классификацию с обучением для различных сельскохозяйственных культур.

Целью данной работы является оценка возможностей декомпозиции Клауда–Потье спутниковых радиолокационных данных с двойной поляризацией HH и HV в L-диапазоне длин волн для мониторинга залесения/лесовосстановления.

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

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

Для исследования использованы временные ряды данных спутниковых радиолокаторов ALOS-1 PALSAR-1 (2007–2010 гг.) и ALOS-2 PALSAR-2 (2015–2020 гг.) на восходящей орбите в режиме FBD (Fine Beam Dual, разрешение 10 м/пиксель) на согласованной горизонтальной (HH) и кросс-поляризациях (HV), полученные в сезон положительных среднесуточных температур. Изображения получены в L-диапазоне (длина волны излучения 24 см). Углы падения составляли 37.2° и 33.2° для PALSAR-1 и PALSAR-2 соответственно. Список использованных изображений приведен в табл. 1. Последние цифры в идентификаторе сцены ALOS-2 PALSAR-2 обозначают дату съемки в формате: первые две цифры обозначают год, следующие две цифры месяц, последние – день месяца. Вследствие значительного влияния влажности земных покровов на обратное отражение радиолокаторов, табл. 1 дополнена данными об осадках в день съемки (ERA5 Daily Aggregates, 2023), полученными из каталога Google Earth Engine (Gorelick et al., 2017).

ВЫДЕЛЕНИЕ УЧАСТКОВ ЛЕСНОГО ПОДРОСТА ПО СТЕПЕНИ ПЛОТНОСТИ ДРЕВОСТОЯ

Как было показано ранее в (Дмитриев и др., 2022), залесение тестового полигона происходило неравномерно, что визуально можно наблюдать и на рис. 1. Первой интенсивно начала зарастать восточная часть – Участок 1, обозначенная синим полигоном на рис. 2, а. Затем – западная часть, на которой наблюдается различие в плотности молодняка, а также наличие не зарастающих мест. Поэтому, для уменьшения ошибок при поляриметрическом анализе лесного подроста, необходимо предварительно сегментировать тестовый полигон, по возможности исключая безлесные участки. Для этого было получено и проанализировано изображение, значения пикселов которого представляли собой минимальные значения интенсивности обратного радарного рассеяния на кросс-поляризации из всего временного ряда данных, представленных в табл. 1. Это позволило выделить в один класс пикселы с наиболее низкими величинами обратного радарного рассеяния – безлесные участки, которые затем были исключены из рассмотрения. В другой класс попадают фрагменты с наибольшими значениями обратного радарного рассеяния – молодой плотный лес, расположенный в западной части тестового полигона (Участок 2 – пурпурные полигоны). Также выделяется класс с промежуточными величинами минимальной интенсивности, представляющий собой лесной подрост, не формирующий сплошной полог, т.е. редколесье (Участок 3 – желтые полигоны). Результаты сегментации представлены на рис. 2, а. В качестве контрольного участка для сравнения также будем использовать фрагмент зрелого плотного леса, который оставался неизменным на всем протяжении периода наблюдений (Участок 4, синий контур).

Для валидации результатов сегментации был дополнительно проведен анализ изменения индекса NDVI за период между 2007 и 2020 гг. С помощью функционала Google Earth Engine по данным Landsat-5, -8 было получено изображение (см. рис. 2, б), представляющее собой разность усредненных 32-дневных значений NDVI, взятых в зимний сезон 2020 г. и аналогичный сезон 2007 г. Зимний сезон был выбран с целью оценки динамики роста только хвойных насаждений. В целом, результаты анализа вегетационного индекса согласуются с результатами анализа обратного радарного рассеяния.

ДЕКОМПОЗИЦИЯ КЛАУДА–ПОТЬЕ

Физической основой радиолокационной поляриметрии является зависимость поляризационного состояния радарного сигнала от физического механизма обратного рассеяния земными объектами. Как правило, поляриметрические данные представлены в виде многоканального изображения, в каждом канале которого находится один из элементов матрицы рассеяния. С помощью матрицы рассеяния можно полностью описать детерминированные (точечные) рассеиватели, которые изменяют поляризацию падающей волны. Однако для распределенных по площади рассеивателей этого сделать невозможно. Поэтому были разработаны различные методы ее разложения (декомпозиции), основной целью которых является представление отраженного от распределенной цели радиолокационного сигнала в виде вкладов элементарных механизмов рассеяния (Cloude, Pottier, 1997; Freeman, Durden, 1998; Krogager et al., 1997; Touzi, 2007; Yamaguchi et al., 2005). Например, такими механизмами являются однократное (от поверхности), объемное (от разнонаправленных неоднородностей растительности) и двухкратное (последовательно от вертикальной структуры затем от поверхности, а также в обратном порядке – от поверхности и затем от структуры, ориентированной вертикально).

Метод Клауда–Потье (Cloude, Pottier, 1997) является математической моделью, которая для данных с полной поляризацией выявляет основной механизм рассеяния из восьми возможных типов. Типы выделяются в соответствии со значениями параметрического угла α и энтропии H, получаемыми из собственных значений матрицы когерентности. Метод Клауда–Потье широко применяется для интерпретации, классификации и сегментации радиолокационных данных.

Применение метода Клауда–Потье для случая двойной поляризации впервые было дано в (Cloude, 2007), однако в данной работе не было приведено деление на зоны в плоскости H-α. Это было сделано позднее в (Ji, Wu, 2015), в которой показано расположение основных механизмов рассеяния в плоскости H-α (при минимальных величинах энтропии) для полной и двойной поляризации (рис. 3). Как можно заметить на рис. 3, для двойных поляризаций HH-HV и HV-VV основные механизмы рассеяния не различаются. Более того, в отличие от случая с полной поляризацией, когда границы между зонами на плоскости H-α точно разделяют механизмы рассеяния, в случае с двойной поляризацией зоны перекрываются между собой, а границы сильно размыты. Последнее является причиной того, что однозначная классификация по физическим механизмам рассеяния для случая двойной поляризации не возможна. Однако, как указывалось выше, возможна оценка динамики сельскохозяйственных посевов в C-диапазоне длин волн (Guo et al., 2018).

РЕЗУЛЬТАТЫ ДЕКОМПОЗИЦИИ КЛАУДА-ПОТЬЕ ДАННЫХ С ДВОЙНОЙ ПОЛЯРИЗАЦИЕЙ

Рассмотрим результаты анализа временного ряда данных радаров ALOS-1, -2 PALSAR-1, -2 с согласованной горизонтальной и кросс-поляризациями. Предварительная обработка (импорт, калибровка, фильтрация спекл-шума) и декомпозиция Клауда–Потье были выполнены в программном обеспечении SNAP (SNAP, 2023). Как и следовало ожидать с учетом работы (Ji, Wu, 2015), расположение результатов декомпозиции на плоскости H-α всех имеющихся данных (18 изображений, см. табл. 1) показало существенное взаимное перекрытие. Поэтому для предварительного анализа изменений на тестовых участках результаты применения декомпозиции представлены только за две даты: 11 июля 2007 г. (синий цвет на рис. 4, а–4 в) и 15 июня 2020 г. (оранжевый цвет на рис. 4, а–4 в). Границы зон в плоскости H-α представлены для случая двойной поляризации в соответствие c (Ji, Wu, 2015).

В целом, представленный на рис. 4, результат показывает отсутствие значительных изменений механизмов радарного рассеяния для всех трех участков, несмотря на значительное увеличение NDVI (см. рис. 2, б) с 2007 по 2020 г. Это связано в том числе с обнаруженным в работе (Ji, Wu, 2015) эффектом смешивания физических механизмов рассеяния и соответственно отсутствием четких границ в отличие от варианта декомпозиции полностью поляриметрических радиолокационных данных. Однако, на указанном временном интервале изменения поляриметрических характеристик сравнимы с флуктуациями облака точек тестовых участков. Поэтому проанализируем статистические характеристики тестовых участков для результатов, показанных на рис. 4. Для этого определим средние величины, стандартное отклонение S и вариативность $Var\left( x \right) = ({{S}_{x}}{\text{/}}\bar {x}) \times 100\% ~$ на участках 1–4 (см. табл. 2). Результаты, представленные в табл. 2 показывают, что энтропия на участках лесного подроста в среднем возросла на 0.08–0.11, на контрольном участке (зрелый лес) на 0.07. Величины разброса (стандартное отклонение и вариативность) наоборот снизились на 0.004–0.037 и на 2–6% соответственно, а на контрольном участке не претерпели изменений. Параметр α также демонстрирует увеличение средних значений на 3.11°–4.63° на участках лесного подроста и на 4.72° на участке зрелого леса.

Для проверки значимости изменений двух средних больших выборок с известными дисперсиями (т.е. изменений на тестовых участках с 2007 по 2020 г.), а также изменений дисперсии в математической статистике используются значения статистического критерия Z_расч и критерия Фишера F_расч, которые определяются по следующим формулам:

где $\bar {x}$ и $\bar {y}$ − выборочные средние значения, S_x и S_y – стандартные отклонения выборок, n_x и n_y – количество элементов в выборках.

Изменения считаются статистически значимыми при выполнении следующих условий (уровень значимости 0.05):

Как видно по табл. 2 значимость увеличения средних значений и уменьшение дисперсии энтропии Н наблюдается для всех участков за исключением контрольного Участка 4. На Участках 1 и 3 также наблюдается значимое снижение величины разброса и вариации параметра α (на 0.44–1.33 и 5–10% соответственно). На 2-м и 4-м участках видно слабое увеличение величины разброса и вариации на 0.13–0.30 и на 0–3% соответственно. Проверка роста дисперсии посредством снижения уровня значимости до 0.01 показала, что критические значения критерии Фишера увеличиваются до 1.06 и 1.19 для 2-го и 4-го участков соответственно, при этом условия значимости увеличения дисперсии на этих участках нарушается. Во остальных случаях значимость не изменилась.

Статистика, приведенная в табл. 2 позволяет утверждать, что достаточно равномерно по всем участкам увеличилась хаотичность (энтропия) обратного радарного рассеяния, что является признаком увеличения количества и размера неоднородностей (стволов и ветвей) соизмеримых с длиной волны радиолокационного эхо-сигнала. Стандартное отклонение для энтропии уменьшается для всех участков лесного подроста, что указывает на уменьшение разнообразия характеристик участков залесения. Так это, например, может показывать уменьшение количества и/или площади безлесных фрагментов и редколесья, т.е. участки к 2020 г. более равномерно заросли молодыми насаждениями. В целом такое сравнение достаточно информативно и демонстрирует возможности декомпозиции при рассмотрении усредненных величин.

Рассмотрим теперь усредненные в пределах каждого тестового участка результаты декомпозиции Клауда–Потье, предполагая, что это позволит избежать многочисленных наложений/перекрытий множества точек, полученных для каждой даты рассматриваемого временного ряда (см. табл. 1). Результаты декомпозиции представлены на рис. 5. На всех графиках красным ромбом выделено среднее (за 2007–2020 гг.) значение для контрольного Участка 4. Это значение является некоторым эталоном, к которому постепенно смещаются средние значения энтропии и угла α лесного подроста. Как можно заметить, усредненные значения для контрольного лесного участка расположены в зоне Z4 (шероховатая поверхность), что противоречит установленному разграничению зон (Ji, Wu, 2015) для декомопозиции с двойной поляризацией в плоскости H-α. Зрелый лес, как правило, объемно рассеивает радиолокационный эхо-сигнал. Соответственно, результаты декомпозиции для зрелого леса должны располагаться в зонах Z5 и Z8 (случайно ориентированные диполи), чего не наблюдается. Вместе с тем на рис. 5 отмечается определенная динамика при использовании средних величин: со временем отмечается переход из зоны Z1 (слабо шероховатая поверхность) в Z4 (шероховатая поверхность). Подобный результат получен в (Guo et al., 2018) для оценки роста сельскохозяйственных посевов по данным C-диапазона с двойной согласованной вертикальной и кросс-поляризациями. Таким образом можно утверждать, что наблюдаемая на рис. 5 динамика изменений параметров H-α показывает увеличение размеров элементов молодого подроста, соизмеримых с длиной волны L-диапазона, а также их количества.

ЗАКЛЮЧЕНИЕ

Анализ результатов поляриметрической декомпозиции Клауда–Потье радиолокационных данных L-диапазона с двойной поляризацией за 11 июля 2007 г. и 15 июня 2020 г. показал, что мониторинг залесения/лесовосстановления на основе классификации по физическим механизмам рассеяния на плоскости H-α невозможен. Поэтому была выполнена статистическая оценка изменений усредненных параметров разложения (энтропии и угла α) за указанные даты. Эти усредненные параметры показывают значимый рост (уровень значимости 0.05) за период наблюдения. Рассмотрение временного ряда таких средних значений для трех тестовых участков залесения на плоскости H-α выявило следующее: с одной стороны, возможна качественная оценка динамики для всех участков, с другой стороны, все три участка, которые относятся по своим отражательным характеристикам к различным классам залесения, демонстрируют сходную динамику по отношению к участку зрелого плотного леса, принятого в качестве контрольного. Вместе с тем следует отметить, что среднее значение для зрелого леса располагается в зоне, соответствующей физическому механизму рассеяния от шероховатой поверхности, вместо зоны с механизмом объемного, в результате многократных переотражений стволами и ветвями деревьев, рассеяния. Это несоответствие позволяет сделать вывод, что предложенное в работе (Ji, Wu, 2015) разделение на зоны различных физических механизмов рассеяния не корректно применительно к оценке лесной среды по данным L-диапазона длин волн. В целом, поляриметрическая декомпозиция Клауда–Потье для случая двойной поляризации совместно со статистическим анализом позволяет уверенно выявить динамику залесения/лесовосстановления.