Методы лесотаксационного дешифрирования радиолокационной спутниковой съемки в горных лесах Южной Сибири

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

Важность решаемых вопросов отражена в основных стратегических и отраслевых документах в сфере развития лесного хозяйства России [1, 2]. Проведенные ранее исследования показывают, что применение методов дешифрирования данных спутниковой съемки для горных лесов имеет значительное различие с аналогичными работами в условиях равнинных лесов [3–5].

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

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

Приведенные методы лесотаксационного дешифрирования спутниковой съемки обоснованы их верификацией по таксационным данным сети из 169 пробных площадей.

Цель исследования – изучить возможность использования радиолокационных данных в решении задач лесотаксационного дешифрирования по определению количественных и качественных характеристик насаждений для эксплуатационных и резервных горных лесов.

Объекты и методы исследований

Наземные исследования для сбора экспериментального материала выполнены в Южно-Сибирской горной зоне: в Алтае-Саянском горно-таежном, Алтае-Саянском горно-лесостепном и Байкальском горном лесных районах. Исследование охватило все выделенные В.Н. Смагиным с соавт. [9] в горах Южной Сибири лесорастительные области: Алтае-Саянскую горную, Центральноазиатскую котловинно-горную, Восточнотувинско-Южноза-байкальскую горную и Прибайкальскую горную.

В Алтае-Саянской горной лесорастительной области пробные площади с отбором образцов были заложены в Северной Алтае-Саянской горной лесорастительной провинции пихтовых и кедровых лесов (Джебашско-Амыльском округе черневых и горно-таежных пихтовых и кедровых лесов и Осевом Западносибирском округе горно-таежных и подгольцово-таежных кедровых лесов), в Восточнокузнецко-Минусинской котловинно-горной лесорастительной провинции лиственничных и сосновых лесов (Кузнецко-Ала-таусском округе подтаежно-лесостепных лиственничных, горно-таежных темнохвойных лесов) и Восточносаянской горной лесорастительной провинции кедровых лесов (Бирюсинско-Китой-ском округе подтаежных сосновых и горно-таежных кедровых лесов).

Для проведения исследований заложено 169 пробных площадей, по которым получены таксационные данные, дано описание типов леса, лесорастительных условий, живого напочвенного покрова и почв (рис. 1). Результаты проведенных за 5 лет наземных исследований позволили сформировать сеть пробных площадей, включающую

Ка(па<**ыа шлалкм >фобиых luomiri ■ Крас1ка|к-»»м крм. Рсса)4лмж Tuu. Pkiij&wk Xaucwa. P

Рис. 1. Карта пробных площадей, заложенных для проведения исследований

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

Пробные площади закладывали по общепринятым в лесоустройстве методам [10] – круговыми площадками постоянного радиуса. Для повышения точности данных радиус площадок в приспевающих, спелых и перестойных насаждениях составлял 20 м, а в средневозрастных – 16 м. В молодняках учёт осуществлялся площадками небольшого радиуса – 1,83 м по диагоналям или параллельным линиям в пределах пробной площади 50*50 м. На пробной площади выполнялась таксация древостоя, описание живого напочвенного покрова и закладка почвенного разреза. Для каждой пробной площади составлена карточка таксационно-дешифровочного участка. При таксации деревьев верхнего яруса, помимо учета всех таксационных показателей, дополнительно определяли их видимость и жизнеспособность. Возраст насаждений оценивался по годичным кольцам на кернах.

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

Следующим шагом предварительной обработки данных РЛИ была радиометрическая калибровка – процесс перевода «сырых» значений пикселей (DN, digital number) к значениям УЭПР. В зависимости от проекции линейного разрешения РЛИ по дальности, существует три варианта нормировки эффективной площади рассеивания (ЭПР) геометрической площадью:

• а)    бета-ноль – ЭПР объекта, нормированная геометрической площадью, у которой линейный размер поперек линии пути находится в плоскости наклонной дальности;

• б)    сигма-ноль – ЭПР объекта, нормированная геометрической площадью, у которой линейный размер поперек линии пути находится в плоскости наземной дальности;

• в)    гамма-ноль – ЭПР объекта, нормированная геометрической площадью, у которой линейный размер поперек линии пути находится в плоскости, перпендикулярной к линии визирования.

С учетом расположения территории исследования ниже 60° с. ш. радиометрическая калибровка осуществлялась в 2 этапа: пересчет DN значения бета-ноль и нормирование до значений гамма-ноль с коррекцией уклонов по цифровой модели рельефа (ЦМР). В последующем осуществлялась координатная привязка радиометрически откалиброванного изображения (ортотрансформирование). В процессе данной процедуры устранялись искажения, связанные с рельефом местности, с помощью ЦМР. Географические координаты записывались в систему WGS84; использовалась SRTM-модель рельефа с шагом пикселя 30 м.

На следующем шаге выполнено подавление спекл-шума – снижение влияния спекл-шума, возникающего на РЛИ вследствие когерентного наложения откликов от нескольких объектов в элементе разрешения изображения. В работе использовался метод некогерентного накопления, заключающийся в усреднении изображения в окне N*M пикселей с последующим укрупнением элементов изображения. Для процедуры выбран размер «окна», равный 2*2, из-за небольшого размера пробных площадей – до 1 га.

Вторым шагом в процедуре предварительной подготовки данных РЛИ было создание и применение маски затенений и переналожений вершин и склонов с целью исключения территорий, на которых возможно получение ошибочных данных для анализа. В горной местности в зависимости от углов визирования для объектов разной высоты и наклона рельефа могут наблюдаться разные смещения склонов гор по наклонной дальности: растяжение, сжатие склонов, переналожение вершин и оснований склонов, радиолокационные тени. Первичная корректировка геометрических искажений РЛИ из-за эффектов сжатия и растяжения склонов осуществляется на этапе ортотрансформирования путем подбора интерполяторов ЦМР высокого порядка – бикубического и на основе sinc-функции. Однако в зонах переналоже-ний и затенений данная коррекция невозможна. Эта ситуация требует исключения из обработки и анализа этих зон путем маскирования. Маска затенений и переналожений строилась на основе данных геореференцированного РЛИ с Sentinel-1, снятого в интерферометрическом широкозахватном режиме (Interferometric Wide Swath, IW), уровня обработки GRDH (Ground Range Detected High Resolution – амплитудное РЛИ в проекции наклонной дальности c высоким уровнем детальности) в модуле SAR-Simulation Terrain Correction программного комплекса SNAP.

Подготовка данных радарной съемки по приведенному алгоритму позволила получить первичные данные для лесотаксационного дешифрирования. Ранее проведенные исследования показывают, что использование только данных спутниковой съемки, без связей с особенностями произрастания насаждений, не позволяет получить количественные и качественные характеристики лесов [8]. Наши исследования позволили разработать методы интеграции факторов, влияющих на продуктивность лесов, в методы лесотаксационного дешифрирования.

Обсуждение

Район проведения исследований характеризуется высоким разнообразием природных факторов, оказывающих влияние на породный состав и продуктивность лесов, процессы лесооб-разования. В частности, район включает зоны тайги, лесостепи, степи и полупустыни, сухой степи. Разнообразие природно-климатических условий осложняется циркуляцией воздушных потоков, в том числе определяющих специфику формирующихся высотно-поясных комплексов. В районе исследований, согласно районированию горных лесов по В.Н. Смагину, выделены 4 лесорастительные области, 14 провинций, 29 округов.

С продвижением в глубь материка в районе исследований уменьшается годовая сумма осадков, происходит резкое колебание температуры воздуха и ее понижение. С продвижением на юг усиливается аридизация климата и снижается количество осадков. Согласно ранее проведенным исследованиям, оптимальные условия произрастания для основных лесообразующих пород наблюдаются на большей части Алтае-Са-янской горной области. Доминирующей формой рельефа является среднегорный эрозионный (800–2 000 м над ур. моря), большая часть которого занята лесами [9]. На высотах 1 500–2 000 м часто встречается рельеф древних поверхностей выравнивания, представленный широкими плато, покрытыми лесом или россыпями камней. Поверхности выравнивания заняты в основном подгольцово-таежными, подгольцово-субаль-пийскими, реже горно-таежными природными комплексами с преобладанием в составе пихты, кедра, иногда лиственницы и ели.

Особенности рельефа и высотная зональность оказывают значительное влияние на природно-климатические факторы и, как следствие, на растительность. Среди элементов рельефа основными факторами, влияющими на лесную растительность, являются высота местности над уровнем моря, экспозиция и крутизна склона [11, 12]. Разнообразие рельефа в горной местности приводит к значительному различию природных комплексов. При этом высота границы леса зависит от широты расположения местности, специфики рельефа.

В среднем с возрастанием абсолютной высоты на 100 м происходит снижение температуры на 0,5 °С, но на северных склонах падение температуры может быть ниже, чем на южных [13].

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

Все перечисленные процессы имеют существенную неоднородность, так как зависят от экспозиции и крутизны склона. Экспозиция склона определяет уровень радиации, распределение осадков, в том числе снега, и скорость ветра [14, 15]. Все эти факторы оказывают влияние на распределение лесообразующих пород в зависимости от экспозиции склона и высотной зональности. Приход солнечной радиации на склоны северной и южной экспозиции сильно различается. В пределах одного климатического пояса в условиях сложного горного рельефа склоны разной экспозиции имеют разный климат, что приводит к различиям и в лесообразовательном процессе. Наиболее высокие температуры отмечаются на юго-западных склонах и, как следствие, для этих склонов характерна пирогенная динамика формирования лесов. Северные склоны прогреваются хуже, на них в основном преобладают темнохвойные породы. Роль экспозиции склонов резко возрастает в районах с сухим климатом, где различия растительных сообществ на горных склонах выражены наиболее отчетливо [16, 17].

Крутизна склона влияет на возможность произрастания растений, образование почвы и ее прогрев. Каждый 1° южного склона в сторону увеличения перемещает условия произрастания на 10° географической широты к югу и, наоборот, каждый градус северного склона – к северу. Даже незначительный склон оказывает влияние на лесную растительность [18–20]. От крутизны склонов зависит степень эрозии и, соответственно, возможность почвообразовательного процесса. При высоких показателях крутизны склона, более 70°, скалы становятся практически недоступными для поселения лесной растительности: она селится в выступах, промежутках между скалами, где есть возможность для формирования почвообразовательного процесса.

Проведённые под руководством профессора В.Н. Смагина работы выявили значительные различия в породном составе и растительности горных лесов юга Сибири в зависимости от широты, долготы, высотной зональности, экспозиции склона, особенностей климата [9].

Обобщение ранее проведенных исследований и анализ результатов экспериментальных работ показывают возможность создания математических моделей роста леса с учетом особенностей рельефа. Их эффективность зависит от специфики воздействия на произрастание лесной растительности региональных факторов, в частности, высоты над уровнем моря, экспозиции и крутизны склона. Перечисленные показатели можно получить из цифровых моделей рельефа, в том числе SRTM (NASA Shuttle Radar Topography Mission) [21, 22]. Оптимальным вариантом для целей лесотаксационного дешифрирования является использование моделей рельефа, полученных по данным лидарной аэрофотосъёмки, что исключает искажение данных в области тени при радиолокационной съемке.

Для решения задач лесотаксационного дешифрирования и разработки математических моделей необходим критерий, которым можно охарактеризовать комплекс изменяющихся условий. На основе работ В.Н. Смагина и проведенных нами исследований таким критерием, комплексно учитывающим воздействие природных факторов, определена верхняя граница леса. Для создания карты границ леса использовались данные со спутника Метеор-М, аппарата комплекса многозональной спутниковой съёмки (КМСС) с полосой захвата 900 км, разрешением 60 м. Анализ спутниковых данных позволил получить уравнение связи верхней границы леса с высотой над уровнем моря, экспозицией и крутизной склона, широтой и долготой. Математическая модель позволяет осуществлять прогноз высоты верхней границы леса для Южно-Сибирской горной зоны (рис. 2).

Зависимость можно выразить уравнением полиномиальной регрессии (1):

ВГЛ = -164240,689044 + 1867,06010105 хХ

+ 3103,25029464xY + 215,646402577хК +

2109,23354366хЭ - 34,9410542175 xXxY —

2,46454303617хХхК +           (1)

0,230633676384xYxK — 31,0209047004хХхЭ

+ 14,0276065139xYxЭ —

79,8877612692хКхЭ + 0xXxYxК +

0xXxYxЭ + 0,949230352913хХхКхЭ — 0,127414214933 xYxКxС + 0xXxYxКxС, где:

ВГЛ – высота верхней границы леса, м;

• X – восточная долгота, градус;

• Y – северная широта, градус;

• К – крутизна склона, градус;

Э – экспозиция склона.

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

Получена модель прогноза класса бонитета насаждений в зависимости от природных факторов. Результаты работы модели представлены

Рис. 2. Фрагмент карты верхней границы леса

на рис. 3. Модель можно описать уравнением факторной регрессии (2):

B = -87,9154224999 + 0,937624218007хЭ + 5,89803386852хК + 0,0914915031603хВНМ

+ 0,0526150661517хВГЛ -

0,41589335311хЭхК —

0,00274025455906хЭхВНМ —

0,00555276512657хКхВНМ —

0,000460769802699хЭхВГЛ—        (2)

0,00359820294032хКхВГЛ -

5,07543622073e - 005хВНМхВГЛ +

0,000428870647473хЭхКхВНМ +

0,000247069532734хЭхКхВГЛ +

1,25719843157e - 006хЭхВНМх ВГЛ3,2077470722е - 006хКхВНМхВГЛ -2,32337283732e - 007хЭхКхВНМхВГЛ, где:

ВГЛ – высота верхней границы леса, м;

К – крутизна склона, градус;

Э – экспозиция склона;

ВНМ – высота над уровнем моря, м.

Оценка наблюдаемых и прогнозных значений показывает высокий уровень корреляции модели (R=0,96, F=12, p=0). На рис. 4 показано,

Ill I g r M * Ы M

Рис. 3. Дешифрирование класса бонитета насаждения на основании верхней границы леса (применение маски границы леса)

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

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

Результаты

С учетом результатов ранее проведенных исследований [3, 4, 8, 23, 24] были определены показатели радиолокационной съемки, имеющие связь с таксационными характеристиками насаждений. Разработка моделей определения запаса, относительной полноты, среднего возраста насаждений осуществлялась с использованием следующих показателей радиолокационной спутниковой съемки:

^ Gamma_VH - УЭПР объекта в значении гамма-ноль на поляризации VH.

^ Gamma_VV - УЭПР объекта в значении гамма-ноль на поляризации VV.

^ RVI - модификация радарного вегетационного индекса для двух поляризаций VV и VH. Представляет собой отношение четырех УЭПР объектов в значении гамма-ноль на поляризации VH к сумме УЭПР объектов в значении гамма-ноль на VV и VH поляризациях.

^ GLCMMean - суммарное среднее. Текстурный признак. Является произведением значения пикселя с частотой его появления с соседним пикселем.

/ GLCMCorrelation - корреляция. Текстурный признак. Мера линейной связи уровней яркостей пар пикселей.

В уравнениях зависимостей применяются обозначения:

B – класс бонитета насаждений;

ВНМ – высота над уровнем моря, м;

ВГЛ – высота верхней границы леса, м;

Э – экспозиция склона;

К – крутизна склона, градус.

В моделях оценка экспозиции осуществлялась согласно классификации: 1 – Север, 2 – Северо-Восток, 3 – Восток, 4 – Юго-Восток, 5 – Юг, 6 – Юго-Запад, 7 – Запад, 8 – Северо-Запад.

Перечисленные показатели используются во всех формулах, представленных ниже.

Реализация расчетов предусматривает использование информации о зонировании территории района исследований по классу бонитета лесных насаждений и показателю верхней границы леса. Зонирование территории по классам бонитета лесных насаждений можно осуществить с использованием материалов лесоустройства, при их отсутствии или неактуальности – по данным цифровых моделей рельефа.

Апробация разработанных методов лесотаксационного дешифрирования проведена на территории Саяно-Шушенского заповедника с определением показателей запаса и полноты насаждений.

Результаты исследований показывают возможность определения запаса насаждений с высокой достоверностью (3). Коэффициент корреляции между наблюдаемыми и прогнозными значениями составляет R = 0,74 (F = 18, p = 0) (рис. 5).

Связь запаса насаждений (М) с показателями радиолокационной съемки характеризуется полиномиальной зависимостью, которую можно представить формулой (3):

М = 583,639211956 – 3938,63692387xGamma_VH + 58308,9862794xGamma_VH~2 + 644,947333587xGamma_VV -1044,7974875 Gamma_VVA2

+ 216,809476031XRVI -

79,0660294597xRVIA2 +

11,0477411441xGLCMMean -

0,0704335968988 GLCMMeanA2 - 195,947 (3) 842693 ^GLCMCorrelation -

61,1594919669xGLCMCorrelation^2 -

198,2247664xB + 16,4666729983 BA2 -

0,0325109965151хВНМ +

2,59467437601е-005хВНМА2 -0,0461101777147ХВГЛ + 5,20766016225e 006 ВГЛА2 -34,0988187075ХЭ + 3,93852859431 ЭА2 + 2,64464290178ХК - 0,0949154300019 КА2.

Рис. 4. Оценка эффективности модели по прогнозу класса бонитета насаждения

Наблюдаемые значения запаса, м3

Рис. 5. Оценка эффективности модели по прогнозу запаса насаждения

Результаты работы модели по формуле представлены на рис. 6.

Полнота насаждения (4) является наиболее трудно определяемым показателем. Уровень

Рис. 6. Дешифрирование запаса насаждений по данным радиолокационной съемки

Рис. 7. Дешифрирование полноты насаждений по данным радиолокационной съемки

корреляции между прогнозными и экспериментальными значениями не превышает 0,5 (F = 5, p = 0). Это может быть связано с тем, что большинство исследований проведено в сложных двухъярусных насаждениях. Отражение сигнала S-диапазона от верхнего полога насаждений, слабое его проникновение под полог ограничивает определение показателя относительной полноты насаждений (Р):

P = 1,43606484464 –

8,93928619974xGamma_VH +

114,169968515xGamma_VH~2 +

2,56959386264xGamma_VV -

6,17665299064*Gamma_W~2 +

0,117400468196xRVI +

0.144535500605 RVK2 +

0,0160772011437xGLCMMean -

0,000117121519459xGLCMMean^2 -

1,51816304489xGLCMCorrelation +     (4)

0,622054320851 Gl .CMCoiTelation^ -

0,192946906605xB +

0,0316702329601 B"2 -

0,00029731827993xBHM +

5,22537524734e - 008xBHM~2 -

0,000311058319558xBrn +

1,55849566358e - 007xBEn~2 -

0,0168759474307xЭ +

0,00148781112477 1K2 +

0,0115133245654 x K0,000419238141324 x K^2.

Результаты работы модели по формуле представлены на рис. 7.

Выводы

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

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