RU2416192C2

RU2416192C2 - Способ определения экологического состояния лесов

Info

Publication number: RU2416192C2
Application number: RU2009115787/21A
Authority: RU
Inventors: Валерий Григорьевич Бондур (RU); Валерий Григорьевич Бондур; Владимир Евгеньевич Воробьев (RU); Владимир Евгеньевич Воробьев; Елена Валентиновна Черепанова (RU); Елена Валентиновна Черепанова; Вячеслав Федорович Давыдов (RU); Вячеслав Федорович Давыдов; Евгений Геннадиевич Комаров (RU); Евгений Геннадиевич Комаров; Вера Алексеевна Фролова (RU); Вера Алексеевна Фролова
Original assignee: Государственное учреждение "Научный центр проблем аэрокосмического мониторинга" - ЦПАМ "АЭРОКОСМОС"
Priority date: 2009-04-27
Filing date: 2009-04-27
Publication date: 2011-04-20
Also published as: RU2009115787A

Abstract

Способ определения экологического состояния лесов включает зондирование с аэрокосмического носителя, измерение коэффициента спектральной яркости (КСЯ) древесного полога, вычисление по измеренным значениям КСЯ индексов жизненности (g) и красного поражения (r), получение двухмерных изображений поля яркости I(х, у) древесного полога в зеленой (G) и красной (R) полосах видимого диапазона, оценку состояния лесных массивов по измеренным значениям КСЯ и I(х, у). Осуществляют синхронное измерение КСЯ и I(х, у) одних и тех же участков соосно установленными на подвижной телеуправляемой платформе спектрометром и цифровыми камерами G и R диапазонов, поля зрения которых совмещены. Вычисляют вегетационный индекс (NDVI) насаждения анализируемого участка в ближнем ИК-диапазоне. Формируют результирующую матрицу (О) изображения путем попиксельного сложения матриц изображений (G) и (R). Программной обработкой сигнала матриц изображений вычисляют среднюю частоту пространственных спектров результирующей матрицы и матрицы (R) соответственно F_cp0, F_cpR, а также площади рельефов древесных пологов S_p0 и S_pR. Рассчитывают характеристики древостоев участка: полноту П₀, П_R и диаметр кроны среднего дерева D₀=1/F_cp0, D_R=1/F_cpR. Определяют весовой коэффициент α через отношения α=(П_R/П₀)(D_R/D₀)^1,2 коэффициента запаса ослабленных деревьев П_RD_R ^1,2 к общему коэффициентному запаса участка П₀D₀ ^1,2. Степень ослабленности Q древостоя участка определяют по регрессионной зависимости, прокалиброванной по измерениям эталонных участков, вида:

Q=0,6(gNDVI)^-1(r(1-NDVI))^1/3(П_R/П₀)(D_R/D₀)^1,2

Способ обеспечивает оперативность и точность оценки экологического состояния лесных массивов. 6 ил., 3 табл.

Description

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

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

Известны методы оценки лесопатологического и санитарного состояния насаждений при наземной (натурной) таксации лесов, такие как: визуальный (глазомерный); выборочно-измерительный; перечислительный метод [см., например, Справочник «Общесоюзные нормативы для таксации лесов», М.: Космос, 1992 г., стр.180-185, Таблицы 60-63, «Шкала оценки категорий состояния деревьев» и «Характеристика категорий состояния насаждений» - аналог].

В способе-аналоге визуальными признаками лесопатологии являются изменение окраски хвои, листьев от зеленого до бледно-зеленого и желтого, изреженность крон, суховершинность, усыхание ветвей. По визуальным признакам лесные массивы подразделяют на категории состояния: от I категории до VI категории.

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

При этом степень ослабленности насаждения оценивают величиной произведения категории состояния (от 1 до 6) на процентное соотношение количества деревьев в насаждении по категориям состояния (% N/N_Σ);

При средневзвешенной величине не более 1,5 насаждение считается здоровым, до 2,5 - ослабленным, до 3,5 - сильно ослабленным и до 4,5 - усыхающим.

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

Известен дистанционный «Способ ранней лесопатологической диагностики», Патент RU №2189732, A01G 23/00, 2002 г. - аналог.

В способе-аналоге получают оцифрованные значения функций яркости I(x,y) в виде матриц дискретных отсчетов (изображений) размерностью (m×n) элементов в зеленой G и красной R зонах спектра, вычисляют попиксельные отношения изображений, полученных в зонах G и R, составляют результирующую матрицу из этих отношений, на результирующем изображении методами пространственного дифференцирования выделяют контуры, рассчитывают функцию фрактальной размерности изображений внутри контуров и по численному значению фрактальной размерности, положению контуров и их площади судят о причинах, координатах и размерах выявленной патологии. Недостатками аналога являются:

- несоответствие получаемых количественных результатов установленным нормативным критериям оценки;

- отсутствие калибровочных характеристик для пересчета расчетных соотношений в нормативные показатели;

- статистическая неустойчивость результатов, поскольку не все признаки лесопатологии учитываются при расчетных оценках.

Ближайшим аналогом является «Способ оценки состояния лесов», патент RU №2038001, кл. A01G 23/00, 1995 г. - аналог.

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

g=G/(B+G+R)

и красного поражения

r=R/(B+G+R),

вычисляют регрессию хроматических коэффициентов

калибруют ее по измерениям контрольных площадок с известными категориями состояния растений на них и оценивают состояние лесного массива в баллах по соотношению текущих расчетных значений g, r. Известный способ имеет следующие недостатки:

- невысокое пространственное разрешение спектрометрических средств (от сотен метров), не позволяющее обнаруживать участки поражения в пределах одного га;

- статическая неустойчивость результатов, поскольку оценка проводится по единственному измеряемому параметру-коэффициенту спектральной яркости.

Задача, решаемая заявленным способом, состоит в реализации дистанционной оценки экологического состояния лесов по установленным нормативным критериям; обеспечении статической устойчивости результатов путем одновременного синхронного измерения нескольких признаков лесопатологии в видимом и ближнем ИК-диапазонах спектра; совместной обработке спектрометрических измерений и изображений, получаемых в зеленой G и красной R полосах видимого диапазона.

Технический результат достигается тем, что в способе определения экологического состояния лесов, включающем их зондирование с аэрокосмического носителя, измерение коэффициента спектральной яркости (КСЯ) древесного полога, вычисление по измеренным значениям КСЯ, индексов жизненности (g) и красного поражения (r), получение двухмерных изображений поля яркости I(x,y) древесного полога в зеленой (G) и красной (R) полосе видимого диапазона, оценку состояния лесных массивов по измеренным значениям КСЯ и I(x,y), дополнительно осуществляют синхронное измерение КСЯ и I(x,y) одних и тех же участков соосно установленными на подвижной телеуправляемой платформе спектрометром и цифровыми камерами G и R диапазонов, поля зрения которых совмещены, вычисляют вегетационный индекс (NDVI) насаждения анализируемого участка в ближнем ИК-диапазоне, формируют результирующую матрицу изображения путем попиксельного сложения матриц изображений (G) и (R), программной обработкой сигнала матриц изображений вычисляют среднюю частоту пространственных спектров результирующей матрицы и матрицы (R) соответственно F_cp0, F_cpR, а также площади рельефов древесных пологов S_p0 и S_pR, рассчитывают характеристики древостоев участка, полноту П₀, П_R и диаметр кроны среднего дерева D₀=1/F_cp0, D_R=1/F_cpR, определяют весовой коэффициент α через отношения α=(П_R/П₀)(D_R/D₀)^1,2 коэффициента запаса ослабленных деревьев П_RD_R ^1,2 к общему коэффициентному запаса участка П₀D₀ ^1,2, а степень ослабленности Q древостоя участка определяют по регрессионной зависимости, прокалиброванной по измерениям эталонных участков, вида:

Q=0,6(gNDVI)^-1(r(1-NDVI))^1/3(П_R/П₀)(D_R/D₀)^1,2

Изобретение поясняется чертежами, где:

Фиг.1 - коэффициенты спектральной яркости: 1 - смешанный лес, 2 - смешанный лес в стадии дигрессии, 3 - ветви сосны, 4 - почва, древесно-подзолистая, 5 - почва, серозем;

Фиг.2 - зависимость относительного объема фитомассы насаждений от NDVI;

Фиг.3 - функция регрессии, категорий состояния лесного массива;

Фиг.4 - зависимость полноты древостоя от отношения площади рельефа древесного полога к геометрической площади участка (S_p/S₀);

Фиг.5 - огибающие пространственного спектра изображений результирующего -1 и полученного в красном диапазоне - 2;

Фиг.6 - функциональная схема устройства, реализующая способ.

Техническая сущность способа состоит в следующем. Для достижения статистической устойчивости, достоверности и точности определения экологического состояния лесного массива необходимо измерить несколько биометрических параметров:

- параметр, характеризующий всю цветовую гамму хвои, листов;

- параметр, характеризующий общий объем фитомассы насаждения;

- параметр, характеризующий структуру (запас, полнота, количество деревьев) и соотношение здоровой и отмирающей растительности внутри общей фитомассы.

В видимом диапазоне длин волн определяющее влияние на спектр отражения и пропускания излучений оказывают хлорофилл и каротины. Хлорофилл имеет главный максимум поглощения в красной зоне спектра на длинах волн от 666 до 689 нм, а также ряд полос поглощения в синей зоне с максимумом поглощения 440-460 нм. Каротин наиболее интенсивно поглощает излучение с длинами волн 455 и 485 нм. Лишь небольшой минимум поглощения пигментами листьев вблизи 550 нм придает листьям (хвое) зеленую окраску. В районе красной границы видимого диапазона с ростом длины волны λ поглощение хлорофиллом и каротином резко падает, в связи с чем в ближнем инфракрасном диапазоне (БИК) отражательная способность листьев (хвои) достигает максимума 40% [см., например, Л.П.Чапурский. Отражательные свойства природных объектов в диапазоне 400-2500 нм. Ч.I, Мин. Обороны СССР, 1966 г., стр.40-51].

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

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

Как следует из графиков Фиг.1, обнаженные ветви сосны, почвы не имеют минимумов поглощения в красной области. Селектируемыми признаками экологического состояния леса по значениям КСЯ являются признаки цветности: зелености, желтизны, а также крутизна графика КСЯ в БИК. Количественной мерой дегрессии лесов по спектрометрическим измерениям (по прототипу) служат хроматические коэффициенты жизненности (зелености) g=ρ₅₅₀/(ρ₅₅₀+ρ₆₅₀) и красного поражения r=ρ₆₅₀/(ρ₅₅₀+ρ₆₅₀). Перечисленные коэффициенты охватывают все стадии дигрессии растительных сообществ. Количественное их изменение в зависимости от степени повреждения происходит монотонно. Признаком, не учтенным в аналогах и прототипе, является крутизна КСЯ в БИК. Для оценки этой крутизны используют так называемый нормализованный дифференциальный вегетационный индекс NDVI [см., например, Сборник «Аэрокосмические методы и геоинформационные системы в лесоведении и лесном хозяйстве», материалы вторго Всероссийского совещания, М.: из-во МГУЛ, 1998 г., стр.119-122, статья Жирин В.М. «Приближенная оценка фитомассы лесного покрова с использованием вегетационного индекса»].

Количественно NDVI характеризует контраст зеленой вегетирующей растительности с другими природными образованиями, равный:

NDVI=(ρ₉₀₀-ρ₆₇₀)/(ρ₉₀₀+ρ₆₇₀)

где ρ₆₇₀ - значения КСЯ в спектральном интервале 0,58-0,68 мкм; ρ₉₀₀ - значения КСЯ в спектральном интервале 0,73-1,1 мкм.

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

Таблица 1
Распределение значений фитомассы лесных земель и соответствующих значений NDVI
Интервалы июльских значений	Насаждения лиственницы	Естест-венные редкие лист- венницы	Заросли ерника и др.	Болота, водные поверхности	Участки горной тундры	Голые камени-стые россыпи	Всего г/га (площадь, %)
0.11-0.285	1.26(4.2)	0.97 (7)	0.26 (8.5)	0.02 (0.4)	3.26(25.1)	2.36 (54.8)	8.13(100)
0.286-0.335	1.82(6.1)	1.32(9.5)	0.29 (9.6)	0.07(1.3)	4.78 (36.8)	1.58(36.7)	9.86(100)
0.356-0.39	3.95 (13.2)	1.78(12.8)	0.44(14.5)	0.2 (3.6)	3.34 (25.7)	0.48(11.2)	14.72(100)
0.391-0.425	3.5(11.7)	2.5 (18.0)	0.34(11.4)	0.28 (5.0)	3.85 (29.6)	1.04 (24.3)	11.51(100)
0.426-0.495	2.75 (9.2)	2.4(17.3)	0.43(14.3)	0.2 (3.6)	4.73 (36.4)	0.82(19.2)	11.33(100)
0.496-0.53	1.7(5.7)	1.35(9.7)	0.5(16.5)	0	7.29(56.1)	0.52(12.0)	11.36(100)
Среднее значение фитомассы т/га типа	29.9	13.9	3.0	5.5	13.0	4.3

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

Чем больше объем фитомассы и индекс зелености, тем более здоровым является насаждение, критерием чего служит произведение (gNDVI). Наоборот, критерием ослабленности является произведение индекса красного поражения r на величину (1-NDVI). Поскольку нормативным показателем категории состояния является число (от 1 до 6), а само значение состояния леса не может измениться дискретно, то функция, характеризующая категорию состояния, должна быть монотонной и непрерывной.

В качестве монотонных и непрерывных функций регрессии математическим институтом им. Стеклова рекомендованы степенные функции. По измерениям эталонных участков и массиву данных [график Фиг.1 и таблица 1] получена калиброванная функция регрессии категорий состояния следующего вида:

K=0,6(gNDVI)^-1(r(1-NDVI))^1/3

Результаты расчетов представлены в примере реализации.

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

Вначале формируют результирующее изображение поля яркости I(x,y) путем попиксельного сложения матриц, полученных в полосах G и R. Тогда в качестве весовых коэффициентов оценки экологического состояния можно использовать отношение морфометрических параметров увядающей растительности (обработкой снимка полосы R) к морфометрическим параметрам результирующего изображения.

Параметрами, характеризующими морфологию лесного участка в целом вычисляемым программной обработкой снимков, являются полнота древостоя П и размеры крон древесного полога D. Полноту насаждения вычисляют через отношение площади рельефа древесного полога анализируемого участка S_p к его геометрической площади S₀ [см. Патент RU №2294622, 2007 г. «Способ определения полноты древостоев»]. График зависимости функции полноты древостоев П от соотношения S_p/S₀ иллюстрируется Фиг.4.

Количество деревьев в лесном массиве и распределение деревьев по диаметрам крон находят по огибающей пространственного спектра изображения поля яркости I(x,y) путем расчета его двумерного Фурье-преобразования [см. Патент RU №2242867, 2004 г. «Способ вычисления запаса лесных массивов»].

Огибающие пространственных спектров результирующего изображения и изображения R иллюстрируются графиком Фиг.5.

Уменьшение охвоенности, суховершинность, сухокронность уменьшают размеры диаметра крон, что сказывается в смещении огибающей пространственного спектра в сторону высоких частот. Диаметр кроны среднего дерева связан с пространственной частотой: D_cp=1/F_cp. От диаметра кроны зависят высота среднего дерева H_cp=7D_ср ^1,2 и ступень толщины.

Количество деревьев в насаждении зависит от полноты. Таким образом, запас насаждения в целом пропорционален двум параметрам, определяемым по изображению: полноте П и диаметру кроны D. Соответственно, весовой коэффициент, который согласно аналогу рассчитывают как отношение запаса ослабленных деревьев к общему запасу лесного массива, может быть представлен как: α=(П_R/П₀)(D_R/D₀)^1,2

Пример реализации способа.

Заявленный способ может быть реализован по схеме, представленной на Фиг.6.

Функциональная схема устройства содержит космический носитель 1 типа научного модуля 77КМЛ, стыкуемого с МКС, с телеуправляемой платформой 2, на которой соосно установлены средства измерений в составе цифрового спектрометрического модуля 3 типа МС-03, цифровая камера 4 с фильтром для съемки в зеленой полосе и камера 5 с фильтром для съемки в красной полосе видимого диапазона типа ДХС-1800. Включение средств измерений над заданным участком лесного массива 6 осуществляет бортовой комплекс управления 7 (PC.МКС) по суточной программе или разовым командам, передаваемым из центра управления полетом 8 по радиолинии управления 9. Отснятые кадры участков леса вместе со служебной информацией (время съемки, координаты, угол визирования платформы) записывают в буферное запоминающее устройство 10 и по каналу связи сбрасывают на наземный пункт приема информации 11. После предварительной обработки (выделение служебных признаков) информацию через ПУП перегоняют в Центр тематической обработки 12, в котором через адаптер 13 она вводится в ПЭВМ 14 для обработки и определения экологического состояния лесных массивов на отснятых кадрах. При обработке используют стандартную комплектацию ПЭВМ в составе процессора 15, оперативного запоминающего устройства 16, винчестера 17, дисплея 18, принтера 19, клавиатуры 20. Результаты обработки отображаются на сайте сети Интернет 21.

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

Таблица 2
Результаты измерений контрольных участков
Категория состояния леса на контрольном	g - индекс зелености	r - индекс желтизны	NDVI	1-NDVI
1	0,67	0,33	0,84	0,16
3	0,45	0,55	0,43	0,57
6	0,43	0,57	0,24	0,76

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

Таблица 3
Исходный массив измерений для получения регрессионной зависимости в виде степенных функций
Категория	g·NDVI	r(1-NDVI)
1	0.565	0.053
3	0.195	0.312
6	0.103	0.43

Записывают систему из трех уравнений с тремя неизвестными a, x, y:

1=a(0,565)^x/(0,053)^y

3=a(0,195)^x/(0,312)^y

6=a(0,103)^x/(0,43)^y

Логарифмированием система степенных уравнений сводится к линейным. Решение системы осуществляют по методу Крамера.

Получены следующие значения неизвестных: a=0,6, x=-1, y=1/3, что не противоречит логике, т.е. чем меньше индекс жизненности (gNDVI), тем более ослабленным является древостой.

Функция регрессии категорий состояния лесных массивов по спектрометрическим измерениям представляется в виде:

К=0,6 (gNDVI)^-1(r(1-NVDI))^1/3

В соответствии с «Общесоюзными нормативами для таксации лесов» [Справочник. М.: Колос, 1992 г., стр.184] степень ослабленности древостоя вычисляют как произведение категории состояния на отношение запаса ослабленных деревьев к общему запасу.

По снимкам цифровых камер в зеленом G и красном R диапазонах эталонного участка третьей категории состояния получены следующие морфометрические параметры древостоя (графики Фиг.4, Фиг.5).

Полнота древостоя П₀=0,6, П_R=0,5, диаметры крон среднего дерева D₀=1/F_cp0=1,66; D_R=1/F_cpR=1.2; (D_R/D₀)^1,2=0,69

Коэффициент α=0,5/0,6·0,69=0,58

По нормативному показателю состояния древостоя лесного массива Q=(3·0.58)=1.78 следует считать ослабленным.

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

Claims

Способ определения экологического состояния лесов, включающий их зондирование с аэрокосмического носителя, измерение коэффициента спектральной яркости (КСЯ) древесного полога, вычисление по измеренным значениям КСЯ индексов жизненности (g) и красного поражения (r), получение двухмерных изображений поля яркости I (х, у) древесного полога в зеленой (G) и красной (R) полосе видимого диапазона, оценку состояния лесных массивов по измеренным значениям КСЯ и I (х, у), отличающийся тем, что осуществляют синхронное измерение КСЯ и I (х, у) одних и тех же участков соосно установленными на подвижной телеуправляемой платформе спектрометром и цифровыми камерами G и R диапазонов, поля зрения которых совмещены, вычисляют вегетационный индекс (NDVI) насаждения анализируемого участка в ближнем ИК диапазоне, формируют результирующую матрицу (О) изображения путем попиксельного сложения матриц изображений (G) и (R), программной обработкой сигнала матриц изображений вычисляют среднюю частоту пространственных спектров результирующей матрицы и матрицы (R) соответственно F_cp0, F_cpR, а также площади рельефов древесных пологов S_p0 и S_pR, рассчитывают характеристики древостоев участка: полноту П₀, П_R и диаметр кроны среднего дерева D₀=1/F_cp0, D_R=1/F_cpR, определяют весовой коэффициент α через отношения α=(П_R/П₀)(D_R/D₀)^1,2 коэффициента запаса ослабленных деревьев П_RD_R ^1,2 к общему коэффициентному запаса участка П₀D₀ ^1,2, а степень ослабленности Q древостоя участка определяют по регрессионной зависимости, прокалиброванной по измерениям эталонных участков, вида:
Q=0,6(gNDVI)^-1(r(1-NDVI))^1/3(П_R/П₀)(D_R/D₀)^1,2.