СПОСОБ ЭКОЛОГИЧЕСКОГО МОНИТОРИНГА ЛЕСОВ Российский патент 2010 года по МПК A01G23/00 

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

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

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

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

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

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

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

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

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

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

g=G/(B+Gh-R)

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

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

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

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

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

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

Технический результат достигается тем, что способ экологического мониторинга лесов, включающий дистанционную регистрацию полей яркости лесной растительности аэрокосмическими средствами, причем дистанционную регистрацию полей яркости лесной растительности осуществляют путем зондирования много- или гиперспектральным датчиком в зеленой G (450-550 нм), красной R (550-670 нм) и ближней инфракрасной БИК (670-950 нм) зонах спектра с одновременным получением цифровых изображений для каждой зоны, вычисляют математическое ожидание сигналов (M_G, M_R, M_БИК) в каждой зоне, формируют матрицу результирующего изображения путем попиксельного сложения изображений G, R, БИК, рассчитывают признаки лесопатологии в виде индекса жизненности g=M_G/(M_G+M_R), индекса поражения R=M_R/(M_G+M_R), нормированного дифференциального индекса продуцирующей фитомассы NDVI=(M_БИК-M_R)/(M_БИК+M_R), площади рельефов древесных пологов изображения R и результирующего изображения, соответственно S_pR, S_p0, среднюю частоту пространственных спектров изображения R и результирующего изображения 0, соответственно F_cpR, F_cp0, а степень ослабленности Q древостоя участка площадью So, определяют по калиброванной эталонной регрессионной зависимости вида:

Q≈0,6(NDVIg)^-1[r(1-NDVI)^1/3(П_R/П₀)(D_R/D₀)^1,2

где П_R, П₀ - расчетная полнота древостоев, вычисляемая через площади рельефов соответствующих матриц П_R=S_PR/S₀, По=S_р0/S₀;

D_R, D₀ - диаметры крон среднего дерева, равного соответственно D_R=1/F_cpR, D₀=1/F_cp0.

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

Фиг.1 - математическое ожидание матрицы цифрового изображения в диапазонах G, R, БИК; 1 - здоровый лес, 2 - смешанный лес в стадии дигрессии;

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

Фиг.3 - функция регрессии, критерий зависимости состояния лесного массива от соотношения индексов (g NDVT) и r(1-NDVI);

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

Фиг.5 - огибающие пространственных спектров изображений: результирующего F_cp0 и R диапазона F_cpR;

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

Техническая сущность способа состоит в следующем

Согласно аналогу [стр.182, табл.60 «Шкала оценки категорий состояния деревьев»] визуальными признаками ослабленности насаждений являются:

- изменение окраски хвои, листьев от насыщенного зеленого цвета до бледно-зеленого, желто-зеленого, желто-оранжевого;

- уменьшение длины хвои, размеров листьев, уменьшение облиственности, охвоенности, объема продуцирующей фитомассы;

- ажурность крон деревьев, укороченный прирост ветвей, разрывы в древесном пологе;

- сухокронность, суховершинность отдельных особей, изреженность древостоя в целом.

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

Изменение цветности насаждений изменяет соотношение между долей поглощаемой и отраженной солнечной радиацией в процессе фотосинтеза в спектральных полосах G и R. Количественной характеристикой эффективности фотосинтеза является расчет индексов зелености g и увядания r через средние значения матриц цифровых изображений: g=М_G/(М_G+М_R) и r=М_R/(М_R+М_R).

Изменение математического ожидания матриц изображений в зависимости от изменения цветности насаждений изображено на графике фиг.1 (1 - здоровый лес, 2 - смешанный лес в стадии дигрессии). Перечисленные коэффициенты охватывают все стадии дигрессии растительных сообществ, количественное их изменение в зависимости от степени угнетения происходит монотонно. Признаком, не учтенным в аналоге и прототипе, является крутизна графика функции при переходе из R полосы в БИК. Для оценки этой крутизны используют так называемый нормализованный дифференциальный вегетационный индекс NDVI [см., например, Сборник, «Аэрокосмические методы и геоинформационные системы в лесоведении и лесном хозяйстве», материалы второго Всероссийского совещания, из-во МГУЛ, М., 1998 г, стр.119-122, статья, Жирин В.М. «Приближенная оценка фитомассы лесного покрова с использованием вегетационного индекса».].

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

NDVI=(M_БИК-M_R)/(M_БИК+M_R).

Распределение значений фитомассы лесных земель (тонн/га абсолютно сухого вещества по июльским измерениям) представлены в таблице 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.

В заявленном способе при расчетах используют метод отношений средних значений цифровых изображений I(х,у), полученных в спектральных зонах G, R, БИК, синхронно, при одних и тех же внешних условиях съемки (высота Солнца, угол визирования, метеоусловия), что обеспечивает исключение систематических ошибок и статическую устойчивость результатов измерений. Количественное изменение категорий состояния насаждения от соотношения перечисленных индексов g, г, NDVI иллюстрируется графиком на Фиг.3. Получение синхронных измерений в зонах G, R, БИК может обеспечивать гиперспектрометр типа «Астрогон-1» [см., например, «Малый космический аппарат «Вулкан-Астрогон» с гиперспектрометром высокого разрешения», эскизный проект, РАКА, НИИЭМ, НТЦ «Реагент», 2002 г., стр.6-14, Гиперспектрометр высокого пространственного разрешения «Астрогон» для решения задач ДЗЗ с КА] либо иные оптико-электронные многоканальные системы.

Гиперспектрометр «Астрогон-1» имеет шесть спектральных диапазонов в интервале от 0,25 до 2,5 мкм со спектральным разрешением 1-3 нм в каждом диапазоне с пространственным разрешение до 1 м и оцифровкой амплитуды сигнала в шкале до 12 двоичных разрядов.

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

Поскольку снимок в каждой спектральной зоне G, R, БИК получают автономно, то массивы изображений следует считать независимыми. Тогда дисперсия результирующего изображения σ_E ²=σ_G ²+σ_R ²+σ_БИК ² существенно превосходит дисперсию слагаемых, чем достигается повышение контрастности результирующего изображения и как следствие точности расчетов.

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

О применимости программного вычисления диаметров крон древесных пологов [см., например, Патент RU №2242867, 2004 г. «Способ вычисления запаса лесных массивов»]. Огибающие пространственных спектров результирующего изображения и R изображения иллюстрируются графиками на фиг.5.

Уменьшение охвоенности, суховершинность, сухокронность уменьшают размер диаметра крон, что сказывается в смещении пространственного спектра в сторону высоких частот. Диаметр кроны среднего дерева связан с пространственной частотой D_cp=1/F_cp. От диаметра кроны зависит высота среднего дерева H_cp=7(D_cp)^1,2 и запас насаждения [м³/га].

Запас насаждения зависит также и от количества деревьев на участке, т.е. от полноты насаждения.

Полноту насаждения вычисляют программным методом, путем расчета площади рельефа древесного полога [см. Патент RU №92294622, 2007 г, «Способ определения полноты древостоев»].

График функции полноты древостоя П от соотношения площади рельефа древесного полога S_p к геометрической площади участка S₀ (равной произведению числа пикселей в матрице на пространственное разрешение одного пикселя Δ) приводится на Фиг.4. Вычисляют соответственно полноту древостоя матрицы результирующего изображения П_о и полноту древостоя матрицы R изображения П_R.

В соответствии с «Общесоюзными нормативами для таксации лесов» [аналог. Справочник, М., Колос, 1992 г, стр.184] степень экологической ослабленности древостоя Q вычисляют как произведение категории состояния (график на фиг.3.) на коэффициент α отношения запаса ослабленных деревьев П_RD_R ^1,2 к общему запасу участка П₀D₀ ^1,2, α=П_R/П₀(D_R/D₀)^1,2.

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

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

Заявленный способ может быть реализован на базе устройства по схеме Фиг.6. Функциональная схема устройства содержит космическую платформу наблюдения 1, типа лабораторного модуля 77 КМЛ, стыкуемого с международной космической станцией (МКС). На космической платформе устанавливается много- или гиперспектральная оптико-электронная камера высокого пространственного разрешения 2, например гиперспектрометр типа «Астрогон-1», осуществляющий съемку участков лесов 3 по командам от бортового комплекса управления (БКУ) 4, на основе программ, закладываемых в БКУ из центра управления полетом 5 по радиолинии командного управления 6.

Получаемые синхронно в трех спектральных зонах G, R, БИК изображения леса заносятся в буферное запоминающее устройство 7. По командам БКУ в зонах радиовидимости космической платформы и наземных пунктов приема информации ППИ 8 данные сбрасываются по мобильному каналу связи 9. После предварительной обработки кадров по служебным признакам (номер витка, время съемки, координаты участка) на ППИ 8 информация через ЦУП передается в центр тематической обработки 10, где через канал связи 11 направляется в ПЭВМ тематической обработки 12, в состав ПЭВМ тематической обработки входит: процессор 13, оперативное ЗУ 14, винчестер 15, дисплей 16, принтер 17, клавиатура 18. Результаты обработки экологического состояния лесных участков выводится на веб-сервер 19.

Исходными для экологического мониторинга являются оценки одного и того же участка лесного массива, одновременно получаемые в трех зонах спектра G, R, БИК.

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

Математическое ожидание сигналов матриц в шкале квантования 0…255 уровней участков 1,2 и 4 категорий состояния и расчетные значения индексов g, r, NDVI представлены в таблице 2.

Таблица 2 Категория состояния M_G М_R M_БИК g r NDVI 1-NDVI Эталон.площ. 1 38 22 196 0.7 0.3 0.8 0.2 2 43 47 190 0.51 0.49 0.6 0.4 4 50 77 182 0.4 0.6 0.41 0.59

Для получения функции регрессии состояния (k) в виде многопараметрической зависимости от лесопатологических индексов, математическим институтом им. Стеклова рекомендованы монотонные степенные функции типа k=а b^x c^y… Будем искать функцию регрессии в виде зависимости k=a(g NDVI)^x[r(1-NDVI]^y.

По данным таблицы 2 исходная система уравнений примет вид:

1=α(0,7 0,8)^x(0,3 0,2)

2=α(0,51 0,6)^x(0,49 0,4)

4=α(0,4 0,41)^x(0,6 0,59)

Логарифмированием система степенных уравнений сводится к линейной. Решение системы осуществлялось по правилу Крамера-Сарриуса. Получены следующие значения неизвестных: а=0,6; х=-1; у=0,3; что не противоречит логике, т.к. чем меньше индекс жизненности (g* NDVI), тем более ослабленным является насаждение (т.е. обратно пропорциональная зависимость). Таким образом, монотонная функция регрессии категорий состояния насаждений, прокалибрированная по изображениям эталонных участков, представляется в аналитическом виде (график на фиг.3) как

k=0,6/(g NDVI)^-1[r(1-NDVI]^0,3

Далее в качестве примера рассмотрим результаты программной обработки изображений участка 2-й категории, они иллюстрируются графиками Фиг.4, Фиг.5.

По графику Фиг.4 определяем среднюю частоту пространственных спектров результирующей матрицы F_ср0=0,75 и матрицы F_cpR=1,25.

Расчетная полнота насаждений, вычисленная по соответствующим матрицам П₀=0,6; П_R=0,5

Коэффициент α=П_R/П₀(D_R/D₀)^1,2=0,5/0,6(D_R/D₀)^1,2

По нормативному показателю состояние насаждений лесного массива Q=(3*0,45), что <1,5. Насаждение следует считать здоровым.

Заявленный способ может быть реализован на существующей технической базе Российского Сегмента Международной Космической Станции (PC МКС) после стыковки с ним лабораторного модуля 77 КМЛ в 2010 г, разработанного РКК «Энергия» им. С.П.Королева на базе отечественных космических аппаратов дистанционного зондирования Земли (ДЭЗ), типа «Ресурс-ДК», на базе иностранных систем ДЭЗ.

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

Способ экологического мониторинга лесов включает дистанционную регистрацию полей яркости лесной растительности аэрокосмическими средствами. Дистанционную регистрацию полей яркости лесной растительности осуществляют путем зондирования много- или гиперспектральным датчиком в зеленой G(450-550 нм), красной R(550-670 нм) и ближней инфракрасной БИК (670-950 нм) зонах спектра с одновременным получением цифровых изображений для каждой зоны. Вычисляют математическое ожидание сигналов (M_G, M_R, М_БИК) в каждой зоне, формируют матрицу результирующего изображения путем попиксельного сложения изображений G, R, БИК. Рассчитывают признаки лесопатологии в виде индекса жизненности g=M_G/(M_G+M_R), индекса поражения R=M_R/(M_G+M_R), нормированного дифференциального индекса продуцирующей фитомассы NDVI=(M_БИК-M_R)/(M_БИК+M_R), площади рельефов древесных пологов изображения R и результирующего изображения, соответственно S_pR, S_p0, среднюю частоту пространственных спектров изображения R и результирующего изображения 0, соответственно F_cpR, F_ср0. Степень ослабленности Q древостоя участка площадью S₀, определяют по калиброванной эталонной регрессионной зависимости вида:

Q≈0,6(NDVIg)^-1[r(1-NDVI)^1/3(П_R/П₀)(D_R/D₀)^1,2

где: П_R, П₀ - расчетная полнота древостоев, вычисляемая через площади рельефов соответствующих матриц П_R=S_рR/S₀, П₀=S_p0/S₀;

D_R, D₀ - диаметры крон среднего дерева, равного соответственно D_R=1/F_срR, D₀=1/F_ср0.

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

Способ экологического мониторинга лесов, включающий дистанционную регистрацию полей яркости лесной растительности аэрокосмическими средствами, причем дистанционную регистрацию полей яркости лесной растительности осуществляют путем зондирования много - или гиперспектральным датчиком в зеленой G (450-550 нм), красной R (550-670 нм) и ближней инфракрасной БИК (670-950 нм) зонах спектра с одновременным получением цифровых изображений для каждой зоны, вычисляют математическое ожидание сигналов (M _G , M _R , М _БИК ) в каждой зоне, формируют матрицу результирующего изображения путем попиксельного сложения изображений G, R, БИК, рассчитывают признаки лесопатологии в виде индекса жизненности g=M _G /(M _G +M _R ), индекса поражения R=M _R /(M _G +M _R ), нормированного дифференциального индекса продуцирующей фитомассы NDVI=(M _БИК -M _R )/(M _БИК +M _R ), площади рельефов древесных пологов изображения R и результирующего изображения, соответственно S _pR , S _p0 , среднюю частоту пространственных спектров изображения R и результирующего изображения 0, соответственно F _cpR , F _cp0 , а степень ослабленности Q древостоя участка площадью So, определяют по калиброванной эталонной регрессионной зависимости вида:
Q≈0,6(NDVIg) ^-1 [r(1-NDVI) ^1/3 (П _R /П ₀ )(D _R /D ₀ ) ^1/2 ,
где П _R , П ₀ - расчетная полнота древостоев, вычисляемая через площади рельефов соответствующих матриц П _R =S _рR /S ₀ , П ₀ =S _p0 /S ₀ ;
D _R , D ₀ - диаметры крон среднего дерева, равного соответственно D _R =1/F _cpR , D ₀ =1/F _cp0 .