Патент на изобретение №2371910

(54) СПОСОБ ОПРЕДЕЛЕНИЯ СОСТАВА НАСАЖДЕНИЙ

(57) Реферат:

Способ включает получение изображения лесных массивов в виде цифровой матрицы из | m×n | элементов зависимости яркости I(х, у) от пространственных координат, расчет пространственного спектра Фурье. Кроме того, осуществляют нахождение средней частоты F_ср и диаметра кроны среднего дерева Д_ср=1/F_ср, определение высоты среднего дерева Н_ср(Д_ср). Также проводят расчет площади рельефа древесного полога S_р и определение полноты насаждения П, вычисление среднего количества деревьев на участке N_ср и нахождение прикрепляющей точки огивы насаждения, расчет запаса V по массовым таблицам при исходных значениях Д_ср, П, N_ср, Н_ср. Дополнительно синхронно осуществляют съемку участка видеоспектрометром на спектральных линиях цветности крон основных древесных пород с получением двумерных изображений в каждом спектральном канале. Вычисляют полноту древостоя (П_i) в каждом спектральном канале через параметры матрицы изображений

и запас отдельной породы V_i=П_iV, по полученной пропорции составляют формулу состава насаждения, где M_i, _i – математическое ожидание и среднеквадратическое отклонение сигнала матрицы изображения i-го спектрального канала; n – количество спектральных каналов (пород). Способ позволит обеспечить дистанционность определения состава насаждения. 6 ил.

Изобретение относится к лесному хозяйству, в частности к лесоустроительным работам с использованием космической съемки для расчета таксационных характеристик. Программная обработка изображений лесных массивов позволяет вычислить (дистанционно) ряд таксационных показателей: полноту насаждений, запас, прирост, бонитет [см., например, Патенты RU 2133565, 1999 г., 2242867, 2004 г, 2277325, 2006 г., 2294622, 2007 г.]. Однако методы дистанционного определения состава насаждений пока отсутствуют.

Известен способ определения состава насаждений путем натурной таксации [см., например, Н.П.Анучин, Лесная таксация, учебник, 5-е издание, Лесная промышленность, М., стр.203-206, § 43 Состав насаждений – аналог.]

В способе – аналоге смешанное насаждение в целом принимается за единицу, а участие в нем каждой породы выражается в десятых долях от единицы. Древесная порода, представленная наибольшей долей, называется преобладающей и ставится в формуле состава на первом месте. Конечной целью таксации леса является установление количества древесины из отдельных пород. Древесина учитывается не по количеству стволов, а по объему, т.е. запасу отдельных пород в виде формулы запаса, типа 6С3Е1Б, что означает 0,6 – запас сосны (С), 0,3 – запас ели (Е), 0,1 – запас березы (Б). Аналогично: Ос – осина, Ол – ольха, Ли – липа, Д – дуб.

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

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

Ближайшим аналогом к заявляемому техническому решению является «Способ вычисления запаса лесных массивов», Патент RU 2242867, 2004 г., A01G, 23/00; G03B 37/00. Способ ближайшего аналога включает операции получения изображения лесов, содержащих пробные площадки, в виде цифровой матрицы зависимости яркости I(х, у) от пространственных координат, расчет пространственного спектра Фурье матрицы, нахождение средней частоты F_cp и диаметра кроны среднего дерева

Д_ср=1/F_ср, а по массовым таблицам высоты и ступени толщины среднего дерева, программное вычисление площади рельефа лесного полога S_рп в виде суммы мозаики аппроксимирующих треугольников и полноты насаждения Р через отношение площади рельефа пробной площадки S_прп к удельной площади рельефа лесного полога S_рп/S₀, домноженного на полноту насаждения пробной площадки Р₀, вычисление среднего количества деревьев на участке N_cp=4S₀Р/ Д_ср² и определение прикрепляющей точки огивы анализируемого насаждения, вычисление запаса по массовым таблицам как суммы 5 классов ступеней толщины Лорея:

где S₀ – геометрическая площадь изображения обрабатываемого участка;

N_i – количество деревьев i-го класса ступени толщины;

Н_i – высота среднего дерева i-го класса ступени толщины, м;

g_i – площадь сечения среднего дерева i-го класса ступени толщины, м²;

f_i – видовое число среднего дерева i-го класса ступени толщины.

Недостатками ближайшего аналога являются:

– не определяется таксационный показатель как состав насаждения;

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

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

Техническое решение задачи состоит в том, что в способе определения состава насаждений, включающем получение изображения лесных массивов в виде цифровой матрицы из | m×n | элементов зависимости яркости I(х, у) от пространственных координат, расчет пространственного спектра Фурье, нахождение средней частоты F_cp и диаметра кроны среднего дерева Д_ср=1/F_cp, определение высоты среднего дерева Н_ср(Д_ср), расчет площади рельефа древесного полога S_p и определение полноты насаждения П, вычисление среднего количества деревьев на участке N_cp и нахождение прикрепляющей точки огивы насаждения, расчет запаса V по массовым таблицам при исходных значениях Д_ср, П, N_cp, Н_ср, дополнительно синхронно осуществляют съемку участка видеоспектрометром на спектральных линиях цветности крон основных древесных пород с получением двумерных изображений в каждом спектральном канале, вычисляют полноту древостоя (П_i) в каждом спектральном канале через параметры матрицы изображений:

и запас отдельной породы V_i=П_iV, по полученной пропорции составляют формулу состава насаждения,

где М_i, _i – математическое ожидание и среднеквадратическое отклонение сигнала матрицы изображения i-го спектрального канала;

n – количество спектральных каналов (пород).

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

фиг.1 – спектральные линии цветности крон основных древесных пород;

фиг.2 – мозаика изображений спектральных каналов (распечатка с ПЭВМ);

фиг.3 – числовые характеристики сигналов (М, ) матриц для расчета полноты древостоев отдельных пород;

фиг.4 – огибающая пространственного спектра Фурье сигнала изображения;

фиг.5 – зависимость полноты насаждения от соотношения площади рельефа S_p древесного полога и геометрической площади участка S₀;

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

Техническая сущность способа состоит в следующем. Отражательные свойства растительных сообществ характеризуются коэффициентом спектральной яркости (КСЯ) [см., например, Л.И.Чапурский «Отражательные свойства природных объектов в диапазоне 4002500 нм, ч.1, Мин. Обороны СССР, 1986 г., стр.40-47, КСЯ растительного покрова]. Спектры отражения крон деревьев формируются совокупными эффектами отражения, поглощения и пропускания лучистой энергии отдельными листьями, ветвями. Определяющее влияние на ход кривых спектрального отражения в видимом диапазоне оказывают хлорофилл и каратиноиды. До 90% лучистой энергии поглощается древесным пологом в процессе фотосинтеза и лишь небольшой максимум отражения в полосе 450550 нм придает растительности зеленую окраску. Внутри зеленой полосы листья, хвоя отличаются по цветности: темно-зеленые (ель, ольха), светло-зеленые (сосна, осина, береза). Цветность основных древесных пород и спектральные линии максимума их отражения в зеленой полосе иллюстрируются фиг.1.

В настоящее время в спектрометрии используются технические средства нового поколения, так называемые приборы «химического» зрения. Последние позволяют проводить измерения на отдельных спектральных линиях (шириной от 0,1 нм до 10 нм) с одновременным получением двумерных изображений в каждом спектральном канале [см., например, «Малый космический аппарат «Вулкан -Астрогон» с гиперспектрометром высокого разрешения». Инженерная записка, РАКА, ФГУП НИИЭМ, 2002 г., стр.9-10, структурная схема камеры «Астрогон-1»]. Таким образом, предоставляется возможность ранжировать кроны древесных пород по цветности и величине КСЯ внутри зеленой полосы видимого диапазона.

Изображения древесного полога, получаемые на спектральных полосах цветности крон различных пород, иллюстрируются мозаикой фрагментов фиг.2. Присутствие в отраженном от древесного полога световом потоке определенной полосы цветности свидетельствует о наличии в насаждениях той или иной породы. За пределами полосы цветности изображение зондируемого лесного массива в спектральном канале измерений имеет малую яркость (темный фон). Долю участия породы в древесном пологе характеризует полнота. По определению [см. аналог, Н.П.Анучин, учебник, стр.234-240, § 49, Полнота насаждений] полнота насаждения характеризует степень использования насаждением занимаемого пространства. Рельеф древесного полога характеризуется функцией зависимости яркости изображения I(х, у) от пространственных координат. Вершины крон освещены лучше и отражают падающий световой поток почти «зеркально». В промежутках между кронами деревьев свет отражается диффузно, поэтому их яркость на изображении меньше. В целом, рельеф древесного полога представляется мозаикой колоколообразных форм, составляющих насаждение древесных пород. Поскольку функция I(х, у) задана цифровой последовательностью дискретных отсчетов, представляется возможным вычислить объемы этих колоколообразных форм по каждой породе. Изображения спектральных каналов имеют единый масштаб и единую площадь участка (одинаковое число элементов | m×n |), поэтому объем занимаемого пространства каждой породой будет пропорционален сумме дискретных отсчетов яркости изображений спектральных каналов. Последнему условию удовлетворяет такая интегральная характеристика сигнала матрицы | m×n |, как математическое ожидание М_i. Другой характеристикой полноты древостоя является плотность стоящих деревьев, размах и взаимное перекрытие крон, т.е. степень изрезанности древесного полога. Данный признак древостоя пропорционален такой интегральной характеристике сигнала изображения, как дисперсия (мощность переменной составляющей). Чем плотнее расположены кроны, тем меньше изрезанность и меньше дисперсия сигнала. Расчет полноты (П_i) древостоя отдельной породы, в соответствующем породе спектральном канале, иллюстрируется фиг.3. Параметром, учитывающем КСЯ породы и степень изрезанности отдельного локуса этой породы в общем лесном массиве, служит отношение математического ожидания М_i к среднеквадратическому отклонению _i сигнала матрицы. Поскольку полнота древостоя занимает интервал |01|, то полученные частные значения П_i надо нормировать. Расчетная полнота П_i древостоя отдельной породы (в i-м канале) определяется как:

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

Алгоритм расчета включает определение размаха яркости I=I_max-I_min в каждом канале, расчет среднего значения яркости М_i и среднеквадратического отклонения _i в каждом канале.

Текст программы определения полноты древостоя в спектральных каналах цветности крон

program spectr;

const NN=10; {Количество матриц}

n=10; m=10;

type iso=array[1..n,1..m] of integer;

mass=array[1..NN] of real;

var sig, sr, rez:mass;

x:isо;

i,j,k:integer; S:real;

function srednee(var x:iso):real;

var i,j:integer; s:real;

begin

S:=0;

for i:=1 to n do

for j:=1 to m do

S:=S+x[i,j] / (n*m);

srednee:=S

end;

function srednee 2(var x:iso):real;

var i,j:integer; s:real;

begin

S:=0;

for i:=l to n do

for j:=1 to m do

S:=S+sqr(x[i,j]) / (n*m);

srednee2:=S

end;

begin

for k:=1 to NN do

begin

{Формирование данных}

for i:=1 to n do

begin

for j:=1 to m do

read(x[i, j]);

readln

end;

{Вычисление sig}

sr[k]:=srednee(x);

sig[k]:=sqrt(srednee2(x)-sqr(sr[k]))

end;

S:=0;

for k:=1 to NN do S:=S+sr[k] /sig[k];

for k:=1 to NN do rez[k]:=sr[k] /sig[k]/S;

for k:=1 to NN do write(rez[k]:10:5)

end.

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

Заявляемый способ может быть реализован на базе схемы, показанной на фиг.6. Схема содержит орбитальный комплекс наблюдения 1 типа космического аппарата (КА) «Ресурс» с установленной на его борту оптико-электронной камерой 2 (типа Геотон-1) и спектрометра 3 (типа Астрогон-1). Съемка запланированных участков местности и включение каналов 2, 3 осуществляет бортовой комплекс управления (БКУ) 4 по командам, передаваемым из Центра управления полетом (ЦУП) 5 по радиолинии управления 6. Информацию изображений лесных массивов 7 записывают на бортовое запоминающее устройство 8 и в сеансах видимости КА с наземных пунктов сбрасывают посредством телеметрической системы 9 по автономному каналу связи 10 на наземные пункты приема информации 11, где осуществляется запись информации на запоминающее устройство 12.

Предварительную обработку информации и формирование матрицы измерений осуществляют в Картографическом центре 13 Министерства Природных Ресурсов. Скомпонованные массивы матриц по запросам потребителей передаются в Региональные центры кадастрового учета лесов 14, где ведется архив 15 насаждений региона на базе стримеров (типа FT-120). Программный расчет площади рельефа древесного полога участков осуществляют на ПЭВМ 16 в стандартном наборе элементов: процессора 17, ОЗУ 18, винчестера 19, дисплея 20, принтера 21, клавиатуры 22. Расчетные значения состава насаждений участков помещают в базу региональных данных с выводом на сайт сети «Интернет» 23. Программу вычисления полноты отдельных пород записывают на винчестер 19.

Многие таксационные показатели, определяемые при оценке насаждений, не получают своего непосредственного отображения на снимках. Скрытые закономерности в строении насаждения и отдельных деревьев могут быть выявлены на основе использования нескольких независимых признаков изображения, таких как цвет, тон, текстура, топология и др. Чем больше независимых признаков изображения использовано при расчете, тем выше точность метода. Текстура изображения содержит информацию о степени изрезанности древесного полога, образуемой кронами деревьев. Крона отдельного дерева при съемке сверху представляется некоторой колоколообразной фигурой. Вершина кроны освещена ярче, чем промежутки между кронами деревьев. Поэтому, текстура изображения повторяет периодичности чередования крон деревьев. Скрытая информация о периодичности чередования и размерах крон деревьев может быть выявлена путем вычисления двумерного спектра Фурье от анализируемой матрицы изображений | m×n | (см., например, «Анализ пространственных частот» в книге Р.Дуда, П.Харт «Распознавание образов и анализ сцен», перевод с англ., Мир, Москва, 1976 г., стр.319-321). Двумерный пространственный спектр G (F_x, F_y) oт функции яркости изображения I(х, у) рассчитывают преобразованием Фурье:

Представленное аналитическое выражение является стандартной математической операцией, реализуемой алгоритмами быстрого Фурье-преобразования (БПФ), входящих в комплект специализированного программного обеспечения: «Пакет программ для обработки изображений в науках о Земле» (см., например, «ER MAPPER Reference», Earth Resource Mapping Pty Ltd, Western Australia, 6005, 1998 г., р.295).

На фиг.4 представлена огибающая пространственного спектра матрицы видеоизображения | m×n | элементов, полученная фотокамерой ДХС во всей полосе видимого диапазона. Полученная функция по физическому смыслу представляет собой величину, обратную линейным размерам объектов подстилающей поверхности, отождествляемых с распределением диаметров крон деревьев D. Расчет огибающей реализуется стандартной программой, входящей в комплект программного обеспечения ER MAPPER 5,0 (см. там же, стр.58, 286). Амплитуда А огибающей в каждой точке F_i определяет удельный вес соответствующей пространственной гармоники в общем спектре. Математическое ожидание (среднюю частоту) пространственного спектра находят из условия:

т.е. F_cp делит площадь под огибающей пополам (см. Патент RU 2183847, 2002 г., стр.7, 8).

При этом расчетный средний диаметр крон анализируемого участка насаждения соответствует: D_cp=1/F_cp (м). Из графика фиг.4 средняя частота F_ср=0,4(1/м), диаметр кроны среднего дерева D_cp=2,5 м. Существует среднестатистическая зависимость между диаметром кроны дерева, его высотой и площадью сечения ствола дерева (см. Патент RU 2133565, 1999 г., таблица 2). На основе анализа обширных экспериментальных данных главнейшие древесные породы разделены на две группы. Первая группа: сосна, лиственница, береза, осина, ольха; вторая группа: ель, пихта, кедр, ясень, бук, кедр. Среднестатистическая зависимость высоты деревьев от диаметра кроны для первой группы H_I7·D^1,2, для второй H_II5·D^1,1, площадь сечения g=120·D^0,8(см²).

Одним из методов исчисления запаса насаждений является способ средней модели [см. аналог Анучин Н.П., Лесная таксация, учебник. Лесная промышленность, 1992 г., стр.296. Способ средней модели]. Запас исчисляется как V, м³=H_cp·g_cp·N_cp·f_cp, где Н_ср – высота модельного дерева, g_cp – площадь сечения ствола, N_cp – среднее число модельных деревьев на участке, f_cp – среднее видовое число модельного дерева.

Среднее число модельных деревьев на участке определяют из соотношения:

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

По определению, полнота насаждения – это степень использования им занимаемого пространства. Характеристикой полноты на изображении является степень изрезанности древесного полога, или отношение площади рельефа древесного полога S_p к геометрической площади участка S₀. Площадь рельефа древесного полога S_p анализируемого участка вычисляют программным методом [текст программы приведен в ближайшем аналоге. Патент RU 2242867, 2004 г., с.7]. Зависимость полноты древостоя от соотношения S_p/S₀ иллюстрируется графиком фиг.5 [см. Патент RU 2294622, 2007 г., Способ определения полноты древостоев]. Для анализируемого участка расчетное соотношение S_p/S₀=l,4, полнота 0,5. Таким образом, среднее число модельных деревьев на одном гектаре (100 м × 100 м) составляет Среднее видовое число модельного дерева [см. аналог, Анучин Н.П., Лесная таксация, стр.105] составляет 0,4. При расчетных данных анализируемого участка H_cp=7·D^1,2_ср=21 м, g_cp=120·D^0,8_ср[см²]=0,252·10^-1 м², N_cp=1020 шт., f_cp=0,4, запас древостоя составит: 216 м³/га.

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

1-й спектральный канал, ель: M₁=176, ₁=22,

2-й спектральный канал, сосна: М₂=112, ₂=37,

3-й спектральный канал, береза: М₃=152, ₃=31.

Откуда полнота древостоев в спектральных каналах, соответственно, составила П₁=0,5, П₂=0,197, П₃=0,31. Запас пород, соответственно, ели 108 м³/га, березы 66 м³/га, сосны 42 м³/га. Формула состава насаждения: 5Е3Б2С.

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

Формула изобретения

Способ определения состава насаждений, включающий получение изображения лесных массивов в виде цифровой матрицы из |m × n| элементов зависимости яркости I(х, у) от пространственных координат, расчет пространственного спектра Фурье, нахождение средней частоты F_cp и диаметра кроны среднего дерева Д_ср=1/F_cp, определение высоты среднего дерева Н_ср(Д_ср), расчет площади рельефа древесного полога S_p и определение полноты насаждения П, вычисление среднего количества деревьев на участке N_cp и нахождение прикрепляющей точки огивы насаждения, расчет запаса V по массовым таблицам при исходных значениях Д_ср, П, N_cp, Н_ср, отличающийся тем, что дополнительно синхронно осуществляют съемку участка видеоспектрометром на спектральных линиях цветности крон основных древесных пород с получением двумерных изображений в каждом спектральном канале, вычисляют полноту древостоя (П_i) в каждом спектральном канале через параметры матрицы изображений

и запас отдельной породы V_i=П_iV, по полученной пропорции составляют формулу состава насаждения,
где M_i, _i – математическое ожидание и среднеквадратическое отклонение сигнала матрицы изображения i-го спектрального канала;
n – количество спектральных каналов (пород).

РИСУНКИ