Spatial anisotropy assessment of the forest vegetation heterogeneity at various azimuth angles of the radar polarimetric sensing

Full Text

ВВЕДЕНИЕ

В последние десятилетия отмечается тенденция увеличения количества спутников с оптической и радиолокационной аппаратурой для решения широкого спектра задач космического мониторинга (Shang et al., 2017; Бондур 2010, 2015; Бондур и др., 2015, 2016; Бондур, Савин, 1992; Бондур, Чимитдоржиев, 2008а, б). Важной областью применения космических средств является картирование земных покровов, том числе лесной растительности, подверженных воздействию антропогенных (Бондур, 2010; Бондур, Воробьев, 2015; Бондур и др., 2016) и других деструктивных факторов, например, таких как природные пожары (Бондур, 2011, 2015; Бондур и др., 2016; Бондур, Гордо, 2018) и их последствия (Бондур, 2015; Бондур, Гинзбург, 2015; Lappalainen et al., 2016). Перспективы использования радиолокационных изображений для решения задач исследования лесной растительности обусловлены высокой проникающей способностью радиолокационных эхо-сигналов в исследуемую среду и возможностью съемки в любое время суток и практически при любой погоде.

Как известно, текстура радиолокационных изображений в силу когерентности сигналов и используемых длин волн является уникальным признаком для оценки состояния и динамики различных типов земной поверхности (Бондур, Старченков, 2001; Бондур, Чимитдоржиев, 2008а). Одним из эффективных методов текстурного анализа является фрактальный подход, количественно характеризующий объекты со сложной пространственной структурой (Mandelbrot, 1982). В последние несколько лет интенсивные исследования были направлены на разработку эмпирических методов для оценки фрактальных размерностей спутниковых изображений как природных, так и искусственных поверхностей, созданных человеком. Например, в работе (Di Martino et al., 2018) рассмотрена роль пространственного разрешения радаров спутников Cosmo/SkyMed при построении карт фрактальной размерности природных сред и показано, что использование данных c разрешением на местности 3–5 м является более информативным. Там же отмечено «не фрактальное поведение» урбанизированных территорий в отличие от природных сред. Оценка особенностей текстуры, без фильтрации спекл-шума, позволила выполнить оценку пространственных изменений (Wang et al., 2016), сегментацию пикселов (Shang et al., 2017) и выявить разрушенные здания (Sun et al., 2016).

Наряду с развитием технологий текстурного анализа радарных данных развиваются другие методы радиолокационного зондирования, такие как интерферометрия, поляриметрия и томография. Значимые результаты получены при комплексном использовании этих методов, например, в поляриметрической томографии (Aghababaee, Sahebi, 2018) и поляриметрическом текстурном анализе (Weissgerber et al., 2016). Анализ текстуры путем фрактального анализа используется в интерферометрии для оценки сдвигов и последующего совмещения интерферометрических изображений (Danudirdjo, Hirose, 2013, 2015). По аналогии с указанными выше методами комплексирования можно предположить, что поляриметрический анализ текстуры радиолокационных изображений позволит получить дополнительную информацию об объектах исследования, в том числе лесной растительности.

Одним из методов поляризационного анализа является графическое представление значений коэффициента обратного радарного рассеяния для всех возможных состояний поляризационного эллипса, реализуемое в виде поляризационной сигнатуры (Van Zyl et al., 1987). Последняя является трехмерным графическим отображением поляриметрической информации об обратном радарном рассеянии (яркости поляриметрических радарных изображений): независимые переменные, представляющие собой угол эллиптичности и угол наклона поляризационного эллипса, откладываются по осям абсцисс и ординат, а аппликатой является интенсивность обратного радарного эхосигнала. Такой подход основан на усреднении величин обратного радарного рассеяния для некоторой совокупности соседних пикселов, что неполно представляет сложные явления, сопровождающие процессы рассеяния и переотражения радиоволн в неоднородных средах. Для оценки степени неоднородности лесной растительности следует оценить пространственные вариации радарного эхосигнала.

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

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

Фрактальный характер природных объектов может проявляться в измерениях различных параметров, описывающих их состояние. Рассмотрим в качестве примера стандартные профили фрагментов радиолокационных изображений, соответствующих лесным массивам и безлесным участкам. На рис. 1 представлены горизонтальные профили радиолокационных сигналов, соответствующих изображению, полученному с помощью аппаратуры спутника ALOS PALSAR‑1 на различных поляризациях.

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

 

Рис. 1. Профили радиолокационного эхо-сигнала на HH и HV поляризациях для РСА ALOS PALSAR‑1.

 

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

Радиолокационное изображение или его фрагмент (см., например, рис. 2а), можно представить в трехмерном виде (рис. 2б), где значения аппликаты будут характеризоваться величинами яркости (величинами ЭПР) радарного изображения.

В настоящее время разработано множество различных методов и алгоритмов расчета фрактальных размерностей, которые рассмотрены, например, в работах (Дмитриев и др., 2016; Чимитдоржиев и др., 2007; Dmitriev et al., 2016).

Рассмотрим в качестве сигнатуры радарных изображений лакунарность , которая также, как и фрактальная размерность введена Мандельбротом и является важным дополнением при анализе текстуры на основе фрактального подхода (Mandelbrot, 1982; Allain, Cloitre, 1991). Некоторые объекты могут иметь одинаковую фрактальную размерность, при этом обладая различной текстурой. Различить такие объекты и помогает лакунарность (Mandelbrot, 1982). Если рассматривать радарное изображение в трехмерном виде, то лакунарность связана с пространственным распределением размеров промежутков или понижений (впадин). Текстуры с низкой лакунарностью считаются однородными, поскольку все размеры промежутков одинаковы; и наоборот, текстуры с высокой лакунарностью неоднородны. Следует отметить, что текстуры однородные на небольших фрагментах, могут быть неоднородными на более крупных площадях и наоборот. Поэтому Λ можно считать масштабно-зависимой мерой текстуры.

Общий подход к оценке лакунарности состоит в оценке распределения масс заданного множества. Так, в работе (Allain, Cloitre, 1991) был предложен метод скользящих кубов, заключающийся в следующем. Использовался куб радиуса r, скользящий по сетке, покрывающей множество. Пусть n(M, r) число скользящих кубов радиуса r и массы M. Определим функцию вероятности распределения масс Q(M, r) как частное от деления величины n(M, r) на общее число кубов. Тогда лакунарность на масштабе r определяется следующим выражением (Allain, Cloitre, 1991):

$\Lambda \left(r\right)=\frac{{\displaystyle {\sum }_M{M}^{2}Q(M,r)}}{{\left[{\displaystyle {\sum }_{M}MQ(M,r)}\right]}^2}.$  (1)

 

Рис. 2. а — фрагмент радарного изображения в дециметровом диапазоне длин волн; б — его трехмерное представление.

 

Построение изображений лакунарности производилось в соответствии с алгоритмом, описанным в работе (Dong, 2000). Данный алгоритм основан на разностном методе подсчета кубов. Куб размером r × r × r (r = 2, 3, 4, 5, …) помещается в левом верхнем углу выделенного участка (окна) размером W × W на исходном изображении. W — нечетное число, удовлетворяющее условию r < W.

 

В зависимости от значений пикселов в скользящем окне размера r × r может потребоваться колонна высотой более одного куба для того, чтобы покрыть поверхность, составленную из пикселов изображения. Присвоим номера 1, 2, 3, … кубам начиная с нижнего. Пусть для каждого скользящего окна размера r × r минимальное и максимальное значения пикселов попадают в куб под номером u и v соответственно. Тогда относительная высота колонны:

 

$n_r\left(i,j\right)=v-u-1\text{\hspace{0.17em}},$ (2)

где i и j — координаты пиксела.

Когда скользящее окно движется в окне W×W, то имеем

$M_r=\sum _{i,j}{n}_r\left(i,j\right)$ . (3)

Для вычисления лакунарности в окне W × W масса М в выражении (1) заменяется массой Mr из (3). Далее, значение лакунарности присваивается центральному пикселу данного окна, а само окно движется через все изображение. Блок-схема построения поляризационной сигнатуры лакунарности представлена на рис. 3.

 

Рис. 3. Блок-схема построения сигнатуры лакунарности.

 

На первом шаге построения лакунарности диапазон изменения углов ориентации ψ и эллиптичности χ поляризационного эллипса разбивается на равные промежутки

ψ_i = ψ_{i −1} + Δψ (i = 0, 1, …, N; ψ0 = 0°; ψ_N  = 180°) и

χ_j = χ_{j −1} + Δχ (j = 0, 1, …, M; χ0 = –45°; χ_M = +45°). (4)

Здесь Δψ и Δχ фиксированные шаги изменения углов.

Для каждой комбинации углов (ψ_i, χ_j) в соответствии с работой (Van Zyl et al., 1987) вычисляется со- или кросс-поляризационная ЭПР. Затем, рассчитанные ЭПР преобразуются в изображения лакунарности с помощью алгоритма, описанного выше.

На полученных изображениях определяется область, для которой нужно построить сигнатуру. Значения, попадающие в данную область, усредняются для получения результирующей лакунарности. Таким образом, имеем одно значение Λ_ij для каждой пары углов (ψ_i, χ_j). Результирующий набор из L = (N + 1) · (M + 1) точек (ψ_i, χ_j, Λ_ij) образует поляризационную сигнатуру лакунарности (Dmitriev et al., 2018).

Для построения сигнатуры лакунарности разработано программное обеспечение на языке программирования С++, которое доступно на GitHub (PolSARSignatures, 2019). Оно основано на Orfeo ToolBox — свободном программном обеспечении для обработки данных дистанционного зондирования (Orfeo ToolBox, 2019).

ОПИСАНИЕ ТЕСТОВОГО ПОЛИГОНА И ЭКСПЕРИМЕНТАЛЬНЫХ ДАННЫХ

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

В силу высокой проникающей способности радиолокационного эхо-сигнала L‑диапазона для исследования были использованы данные космических поляриметрических РСА: спутника SIR-C, выполнявшего съемку в октябре 1994 г., и спутника ALOS, выполнявшего съемку в 2006–2009 гг. (табл. 1). Необходимо отметить отсутствие полностью поляриметрических радиолокационных изображений современного спутникового радара L‑диапазона ALOS‑2 PALSAR‑2 для исследуемого тестового участка.

Изображения, полученные РСА SIR-C формировались одновременно в двух частотных диапазонах (L и C с длинами волн λ = 24 см и λ = 5.6 см соответственно) на восходящей орбите с углом наклонения 62.6° и углом падения волны 23.9°. При этом зондирование выполнялось по направлению на юго-юго-восток (стрелка 1 на рис. 4).

Угол наклонения орбиты ALOS PALSAR‑1 составляет 98.1°, угол падения радиоволн — 24°; при зондировании на восходящей орбите направление зондирования близко к восточному (стрелка 2), а на нисходящей орбите направление зондирования изменяется на западное (стрелка 3).

С целью устранения спекл-шума для исходных радиолокационных данных проводилось некогерентное накопление. В результате размер пиксела на изображениях составил 24×24м.

 

Рис. 3. Блок-схема построения сигнатуры лакунарности.

 

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

SIR-C	ALOS PALSAR‑1 восходящая орбита	ALOS PALSAR‑1 нисходящая орбита
09.10.1994	28.06.2006	30.05.2006
10.10.1994	13.08.2006	15.07.2006
	28.09.2006	30.08.2006
	13.11.2006	15.10.2006
	31.03.2007
	16.05.2007
	16.11.2006
	02.04.2008
	05.04.2009

 

Рис. 5. Результаты наземных измерений азимутального распределения ветвей деревьев.

 

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

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

На пробных выделах было проведено таксационное описание и измерено азимутальное распределение неоднородностей в виде ветвей. На тестовой выборке для каждого из 100 деревьев определялся азимут направления, соответствующего максимальному количеству ветвей. Результаты этих измерений показали, что у 60–75% деревьев наибольшая густота ветвей наблюдается в диапазоне азимутальных углов от юго-востока до юго-запада (см. рис. 5).

Была выполнена также визуальная оценка наклона ветвей деревьев: распределение в промежутке углов 0°÷180° от горизонтали. В средней части стволов преимущественные наклоны крупных ветвей составляли 10°÷50°. Вообще, структура наклонов ветвей в вертикальной плоскости достаточно разнообразна и зависит от локальной плотности деревьев.

РЕЗУЛЬТАТЫ ИССЛЕДОВАНИЙ И ИХ АНАЛИЗ

На рис. 6 представлены классические сополяризованные сигнатуры (а, б, в), характеризующие усредненные, в определенном «окне», величины обратного радарного рассеяния.

Также построены сополяризованные сигнатуры:

• на основе второго момента (г, д, е);
• с использованием фрактального анализа — фрактальной размерности (ж, з, и) (Dmitriev et al., 2016) и лакунарности (к, л, м) (Dmitriev et al., 2018).

Они описывают пространственные вариации яркости радиолокационных изображений, по сути пространственные вариации сигналов обратного радарного рассеяния.

На рис. 6 по оси абсцисс отложен угол ориентации эллипса поляризации, по оси ординат — угол эллиптичности, по оси аппликат для классической сигнатуры — усредненная нормированная величина обратного радарного рассеяния, для фрактальной — фрактальная размерность, а также соответственно нормированная величина второго момента и величина лакунарности. Представленные сигнатуры получены усреднением данных раздельно для трех направлений радиолокационного зондирования.

 

Рис. 6. Сополяризованные сигнатуры: а–в — классические, г–е — сигнатуры второго момента, ж–и — фрактальные, к–м — лакунарности для различных орбит.

 

Традиционные поляризационные сигнатуры (см. рис. 6 а–в), полученные для трех различных направлений зондирования сходны по форме и соответственно одинаково характеризуют усредненное поляриметрическое радарное рассеяние тестовым участком леса. Сигнатуры второго момента (см. рис. 6 г–е) аналогичны традиционным и также сходны для трех случаев зондирования. Представленные традиционные поляризационные и сигнатуры второго момента являются хорошо изученными сигнатурами, характерными для двухгранного уголкового отражателя (VanZyl et al., 1987).

Рассмотрим фрактальную сигнатуру в зависимости от угла наклона поляризации, полагая что зависимость от угла эллиптичности симметрична относительно линейной поляризации (0° по оси ординат). Как отмечалось ранее в работе (Чимитдоржиев Т.Н. и др., 2007) фрактальная сигнатура, полученная по данным SIR-C (рис. 6ж), имеет симметричное распределение относительно угла наклона поляризации 90°, который соответствует вертикальной согласованной поляризации.

Значения фрактальной размерности более 2.85 (Бондур, Чимитдоржиев, 2008а) соответствуют значительным пространственным флуктуациям радиолокационного эхо-сигнала. При этом эти флуктуации коррелируют с пространственным распределением неоднородностей в виде стволов и крупных ветвей деревьев (Бондур, Чимитдоржиев, 2008а). Определенная симметричность сигнатуры относительно угла ориентации эллипса 90° (вертикальная поляризация), по-видимому, связана с аналогичным распределением ветвей с западной и восточной стороны деревьев при зондировании в южном направлении (в случае SIR-C).

Для радиолокационных изображений, полученных с борта спутника ALOS PALSAR‑1 при зондировании на восходящей орбите, т. е. в направлении на восток, получена фрактальная размерность более 2.97 в диапазоне изменения углов наклона эллипса поляризации 90° ÷ 170° (см. рис. 6з), т. е. наблюдается увеличение пространственных вариаций интенсивности обратного радарного рассеяния. При данной геометрии зондирования (см. стрелка 2 на рис. 4), когда нулевой угол ориентации эллипса поляризации совпадает с направлением на юг, подобное угловое распределение свидетельствует о том, что с южной стороны деревьев ветви образуют более густой слой отражения, а с северной стороны ветви располагаются реже, т. е. наблюдается анизотропия пространственного распределения ветвей.

Напротив, на нисходящих орбитах при зондировании в сторону запада (нулевой угол ориентации эллипса поляризации совпадает с направлением на север) увеличение пространственных вариаций наблюдается при углах ориентации эллипса поляризации 20° ÷ 90° (рис. 6и). Это также может быть обусловлено большей густотой ветвей, растущих в южную сторону.

 

Рис. 7. Гистограммы распределения третьего параметра Стокса для тестового участка леса: а — восходящая орбита (дата съемки: 28.09.2006), б — нисходящая орбита (дата съемки: 30.08.2006).

 

Сравним полученные выше результаты с использованием данных спутника ALOS PALSAR, с одним из традиционных способов оценки состояния поляризации. Рассмотрим гистограммы распределения значений 3-го параметра Стокса по изображениям тестового участка (см. рис. 7). Третий параметр Стокса S₂ характеризует соотношение компонент, линейно поляризованных под углами (с азимутами) +45° или –45°. Существуют следующие три случая для величин S₂:

1) S2 > 0, то волна поляризована преимущественно с азимутом +45°,

2) S2 < 0, то волна поляризована преимущественно с азимутом –45°,

3) S2 = 0, то в волне не преобладает ни одна из этих двух поляризаций.

Расчет изображений параметров Стокса выполнен с использованием программы PolSARpro version 5.0.4 при размере окна 5 × 5.

Гистограмма для восходящей орбиты, приведенная на рис. 7, демонстрирует некоторое преобладание положительных значений 3-го параметра Стокса, т. е. компонента линейно поляризованного обратного радарного рассеяния с углом наклона +45 градусов больше чем компонента с углом наклона поляризации –45 градусов.

При геометрии зондирования характерной для восходящей орбиты (стрелка 2 на рис. 4) это может свидетельствовать о большей деполяризации и большем количестве элементарных рассеивателей, расположенных под углом 45° с южной стороны деревьев.

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

Сигнатура лакунарности, представленная на рис. 6 к–м, также демонстрирует некоторую асимметричность относительно угла наклона 90 градусов, что особенно заметно для случая восходящей орбиты. Это подтверждает описанный выше эффект азимутальной анизотропии неоднородностей леса. Помимо этого, площадные максимумы красного цвета (см. рис. 6) могут характеризовать так называемую «степень заполняемости» пространства и соответствовать углам наклона ветвей деревьев.

ЗАКЛЮЧЕНИЕ

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

Предложенная новая поляризационная сигнатура для учета пространственных вариаций обратного радарного рассеяния — лакунарность, в отличие от фрактальной сигнатуры, характеризует степень неоднородности пространственных вариаций регистрируемых сигналов. Поляриметрическая сигнатура лакунарности позволяет определить “степень заполняемости” пространства, которая используется при анализе лесной растительности. Данная характеристика является дополнительной к сигнатуре фрактальной размерности. В целом такие сигнатуры более полно описывают пространственные вариации обратного радарного рассеяния лесной растительности при различных состояниях поляризационного эллипса.

Данные ALOS PALSAR-1 и ALOS‑2 PALSAR‑2 предоставлены японским аэрокосмическим агентством JAXA в рамках проекта ALOS‑2 RA‑6 (PI 3402 и PI 3092).

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

Работа выполнена в рамках апробации методов и технологий, разрабатываемых по государственному заданию ИФМ СО РАН, а также при финансовой поддержке Минобрнауки России (уникальный идентификатор проекта RFMEFI58317X0061).

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 of the Siberian Branch of the Russian Academy of Sciences

Email: tchimit@gmail.com
Russian Federation, Ulan-Ude

A. V. Dmitriev

Institute of Physical Materials Science of the Siberian Branch of the Russian Academy of Sciences

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

P. N. Dagurov

Institute of Physical Materials Science of the Siberian Branch of the Russian Academy of Sciences

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

References

Supplementary files