Робастные оценки характеристик древостоев по изображениям видимого диапазона Текст научной статьи по специальности «Электротехника, электронная техника, информационные технологии»

ИНСТИТУТ СИСТЕМНЫХ ИССЛЕДОВАНИЙ ЛЕСА

РОБАСТНЫЕ ОЦЕНКИ ХАРАКТЕРИСТИК ДРЕВОСТОЕВ ПО ИЗОБРАЖЕНИЯМ ВИДИМОГО ДИАПАЗОНА

В.Ф. ДАВЫДОВ, академик РАЕН, проф. каф. БЖДМГУЛ, канд. техн. наук,

Ю.П. БАТЫРЕВ, ст. научн. сотр. ИсИЛМГУЛ, канд. техн. наук

Лесной массив (древостой) принято характеризовать набором таксационных характеристик, называемых элементами леса: запас, продуктивность, бонитет, состав, полнота и др. [1, 2]. Товарной характеристикой леса при сдаче его в аренду служит запас V (м3/га). В связи с новой арендной практикой ведения лесного хозяйства возникает необходимость оперативной дистанционной оценки запаса древостоев. Существуют различные методы наземной оценки запаса [1], среди которых известен метод среднего модельного дерева. Запас V (м3/га) равен произведению четырех параметров: N - среднего количества модельных деревьев в лесном массиве, Нср - средней высоты (м) модельного дерева, g - площади сечения (м2) модельного дерева на уровне 1,3 м., уср - среднего сбегового числа. С математической точки зрения, чтобы задача дистанционного определения запаса древос-тоев имела однозначное решение, необходимо вычислить как минимум четыре независимых параметра изображения лесного массива.

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

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

batyrev@mgul.ac.ru

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

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

Отображаемым на космическом снимке лесов объектом является древесный полог. Распределение яркостей в сечении древесного полога ортогональной плоскостью иллюстрируется рис. 1. Независимыми элементами древесного полога на изображении видимого диапазона являются:

N - число локальных максимумов на изображении, отождествляемое с количеством деревьев;

a - расстояние между локальными максимумами, отождествляемое с расстоянием между деревьями;

Sp - площадь рельефа древесного полога, содержащая информацию о сомкнутости крон, полноте древостоя;

Д - размах крон деревьев в насаждении, характеризуемый огибающей пространственного спектра функции I(x,y).

При программной обработке функции I(x,y) яркости древесного полога представляется возможным вычисление перечисленных 4-х независимых параметров и однозначной оценки таксационных характеристик по изображению видимого диапазона.

Вычисление количества локальных

максимумов на изображении

Поиск локальных максимумов двумерной функции яркости I(x,y) реализуется алгоритмом, иллюстрируемым рис. 2. Как известно, в точке экстремума производная функции равна нулю. Для нахождения про-

56

ЛЕСНОЙ ВЕСТНИК 7/2011

ИНСТИТУТ СИСТЕМНЫХ ИССЛЕДОВАНИИ ЛЕСА

Рис. 1. Сечение изображения древесного полога ортогональной плоскостью

изводной последовательно от начала массива разбивают матрицу изображения на окна размером 13^31 смежных элемента и рассчитывают производные по координатам x, у относительно центрального пикселя а/ = i (i+j) -1 (i -1, j) dx 2dx

< л = I(i, j +1) -1(i, j -1) ■ dy 2dj

В цифровой матрице приращения заменяют конечными разностями. Приращения по x, у тождественно равны одному пикселю. Последовательно анализируя окна конечных разностей, находят точки, для которых конечные разности стремятся к нулю. Подсчитывают количество таких точек и отождествляют их с числом деревьев N. Рассмотренный алгоритм реализуется специализированной программой, запатентованной авторским патентом № 2359229, 2007 г. [7].

Вычисление среднего расстояния между локальными максимумами

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

расстояние между локальными максимумами (а). С методической точки зрения учет двух параметров обеспечит большую точность. Для расчета среднего расстояния между локальными максимумами создают окно матрицы [m0 х n0] с центром хтекущ. Алгоритм расчета иллюстрируется рис. 3.

На оцифрованном изображении вокруг каждого пикселя располагается 8 смежных пикселей. Расстояние между локальными максимумами вычисляют по 8 направлениям. Поскольку среднее расстояние между деревьями находится между радиусами а3 и а4 [1], то достаточно ограничиться поиском пяти максимумов в «плавающем» окне. Расстояния между локальными максимумами находят по теореме Пифагора, катетами треугольников в которых являются перпендикуляры между строками (столбцами). Результат расчета усредняют по сумме составляющих. Количество деревьев в насаждении N2 рассчитывают из соотношения

N = S.

a

где S0 - геометрическая площадь обрабатываемого участка насаждения. Вычисление среднего расстояния между локальными максимумами реализуется специализированной программой, запатентованной авторским патентом [7]. Среднее число деревьев в насаждении N вычисляют как среднегеометрическое

N = V N • N ■

ЛЕСНОЙ ВЕСТНИК 7/2011

57

ИНСТИТУТ СИСТЕМНЫХ ИССЛЕДОВАНИИ ЛЕСА

столбец j - 1 столбец j столбец j + 1

строка i - 1 i - 1 j - 1 i - 1 j i - 1 j + 1

строка i i j - 1 i j i j+ 1

строка i + 1 i + 1 j - 1 i + 1 j i + 1 j + 1

к

я

<а

Ю

Я

о

н

о

окно 3 х 3

пиксел

строка m

Рис. 2. Алгоритм нахождения локальных максимумов

Определение полноты насаждения по его изображению

По определению [1] полнота насаждения характеризует эффективность использования им занимаемого пространства. Чем

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

58

ЛЕСНОЙ ВЕСТНИК 7/2011

ИНСТИТУТ СИСТЕМНЫХ ИССЛЕДОВАНИИ ЛЕСА

теграл функции I(x,y), заданной в некоторой области матрицей [m х n] элементов [3]

^=Я

1 +

дТ

дх

-|2

+

д!_

ду

-2

dx • dy.

Вычислительная процедура реализуется специализированной математической программой, запатентованной авторами патентом Ru № 2294622, 2007 г. [8, 9].

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

SP/S0 1,1 1,25 1,4 1,56 1,72 1,9 2,1 2,3

П 1 0,9 0,8 0,7 0,6 0,5 0,4 0,3

Используя эталонный ряд, относительную полноту таксируемого массива определяют по расчетному отношению Sp к геометрической площади S0.

Вычисление диаметров крон деревьев по изображению древесного полога

В системе древостоя вследствие взаимодействия особей между собой, конкуренции за свет, воду и питательные вещества пространственная структура подчиняется определенным закономерностям. Эффект взаимодействия деревьев и степень напряженности отношений проявляется в самой структуре древесного ценоза [2].

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

Разреженность или густота полога, диаметры крон, их геометрическая повторяемость участвуют в формировании изображения и его крупных дискретных контрастных элементов, что отражается на качестве текстуры или «шероховатости» изображения. Скрытые периодичности чередования крон

Т а б л и ц а 1 размеры крон основных древесных пород

Возраст, лет Порода Высота, м Диа- метр ствола, см Диаметр кроны, м

средн., м макс, м

С 13,4 10,7 1,4 2,4

Е 11,8 9,8 2,3 4,2

Б 15,9 10,8 1,9 4,0

ОС 15,0 16,4 2,3 4,3

С 18,3 15,9 1,7 3,0

Е 16,6 15,5 3,0 4,8

60 Б 19,3 15,7 2,4 4,5

ОС 19,9 22,1 3,1 4,8

ОЛ 22 24 3,0 4,0

Б 21 15 2,9 5,8

90 ОС 24 22 2,9 4,4

С 29 30 2,8 6,8

С 23,6 24,8 2,7 4,5

Е 23,4 25,2 4,3 6

110 Е 27 34 4,3 8,0

Примечание: На основе анализа обширных экспериментальных данных (Анучин Н.П. Лесная таксация. - М: Лесная промышленность, 1982. - с.250) главнейшие древесные породы разделены на 2-е группы. Первая группа - сосна, лиственница, береза, осина, ольха; вторая группа - ель, пихта, кедр, ясень, бук, дуб. По данным таблицы статистическая зависимость высоты деревьев от диаметра для 1-ой группы Н1 = 7-Д1-2 , для 2-ой группы Н2 = 5-Ди, g =120-Д°’8 (см2)

Рис. 3. Определение среднего расстояния между деревьями

ЛЕСНОИ ВЕСТНИК 7/2011

59

ИНСТИТУТ СИСТЕМНЫХ ИССЛЕДОВАНИИ ЛЕСА

АЧХ

0,2 0,4 0,6 0,8 1 1,2 1,4 1,6 1,8 2

F, 1/ м

Рис. 4. Амплитудно-частотные характеристики пространственных спектров функций яркости I(x,y)

деревьев могут быть выявлены путем вычисления пространственного спектра Фурье анализируемого участка. «Шероховатость» или изрезанность полога, его архитектура определяют частотно-спектральный образ растительности. По определению [4] пространственный спектр Фурье вычисляют как двумерное Фурье-преобразование от функции яркости I(x,y) участка изображения

m n

G( Fx, Fr) = ff I (x, y) x

Г 00 П

x exp I - j2n(Fx ■ x + Fy ■ y) I dx ■ dy .

Вычисление двумерного пространственного спектра Фурье входит в комплект специализированного программного обеспечения ER MAPPER 5.0 «Пакет программ для обработки изображений в науках о Земле».

АЧХ пространственных спектров двух изображений I(x,y) иллюстрируется рис. 4. Технология дистанционного расчета характеристик древостоев запатентована авторскими патентами [7-11].

По АЧХ вычисляют диаметр кроны Дср среднего (модельного) дерева: Дср = 1/F^. Среднюю частоту пространственного спектра находят из соотношения

Fcp Fmax

J A(f)df = f A(f)df .

F

cp

Средневзвешенная частота F делит площадь под кривой АЧХ пополам.

В свою очередь, диаметр крон деревьев связан регрессионными зависимостями

с высотой древостоя и ступенями толщины (табл. 1). Полученные числовые характеристики древесных пологов позволяют рассчитать запас насаждения либо аналитически [1], либо по табличным данным [2].

Библиографический список

1. Анучин, Н.П. Лесная таксация / Н.П. Анучин. - М: Лесная пром-сть, 1982. - С. 198-227.

2. Общесоюзные нормативы для таксации лесов. Справочник / М.: Колос, 1992. - С. 116-118

3. Пискунов, Н.С. Дифференциальное и интегральное исчисления для ВУЗов / Н.С. Пискунов: учебник, Т.2. - М.: Наука, 1964. - С. 73-74.

4. Белов, С.В. Аэрофотосъемка и авиация в лесном хозяйстве: уч. пособие. / Белов С.В., Дмитриев И.Д., Колосова А.Е. - Л: ВЗЛТИ, 1962.

5. Дуда, Р Распознавание образов и анализ сцен. / Дуда Р., Харт П. - М.: Мир, 1976.

6. Измерение продуктивности лесов и характеристик лесопроизрастания по изображениям видимого, ИК и СВЧ диапазона. / Отчет по НИР «УраганЛес». - М.: МГУЛ, 2005. - 59 с.

7. Способ определения количества деревьев в лесном массиве. / Давыдов В.Ф., Корольков А.В. - Патент РФ № 2359229, 2009.

8. Способ определения площади рельефа. / Давыдов В.Ф., Корольков А.В., Чернобровина О.К. - Патент РФ № 2253089, 2005.

9. Способ определения полноты древостоев. / Давыдов В.Ф., Корольков А.В., Новиков Е.К. - Патент РФ № 2294622, 2007.

10. Способ вычисления запаса лесных массивов. / Давыдов В.Ф., Корольков А.В., Шалаев В.С. и др. - Патент РФ № 2242867, 2004.

11. Способ определения состава насаждений. / Давыдов В.Ф., Корольков А.В., Бондур В.Г. и др. - Патент РФ № 2371910, 2009.

60

ЛЕСНОЙ ВЕСТНИК 7/2011