Способ определения полноты древостоев

Изобретение относится к лесному хозяйству, предназначено для дистанционного мониторинга лесов и обработки изображений лесных массивов. Способ определения полноты древостоев заключается в получении изображения лесов, содержащих пробные площадки в виде цифровых матриц зависимости функции яркости изображения - сигнала I(х, у) от пространственных координат путем съемки фотокамерой, установленной на орбитальном комплексе наблюдения, расчете пространственного спектра, определении площади рельефа древесного полога S_p, получении эталонного ряда зависимостей известной относительной полноты пробных площадок от площади их рельефов, интерполировании смежных значений эталонного ряда по соотношению S_p/S_геом таксируемого участка, где S_геом - геометрическая площадь таксируемого участка. Способ позволяет повысить точность и статистическую устойчивость результата вычисления запаса лесного массива. 7 ил.

 

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

Одним из важнейших таксационных показателей насаждений при расчете запаса древостоев является полнота. Различают полноту абсолютную и относительную [см., например, Справочник, "Общесоюзные нормативы для таксации лесов", Колос, Москва, 1992 г., стр.116].

Полнота абсолютная - это сумма площадей сечений на высоте 1,3 м деревьев, элемента леса (яруса, древостоя) на единице площади [м²/га]. Определяется по данным пересчетов или путем закладки реласкопических пробных площадок и измерений полнометром Биттерлиха или призмой Анучина [см., например, Н.П.Анучин, Лесная таксация, учебник, 5^е изд., Лесная промышленность, 1982 г., Метод Биттерлиха, стр.314-318 - аналог]. В способе-аналоге на таксируемой площади выборочно закладываются 4-5 круговых площадок на 1 га, и из центра круговой площадки путем визирования прицелом реласкопа подсчитывают количество деревьев в полосе шкалы прибора, а затем количество подсчитанных через реласкоп деревьев умножают на соответствующий множитель шкалы прибора.

Точность способа (10...20%) зависит от количества заложенных круговых площадок и их размещения по площади участка. Кроме того, способ применим только при натурных измерениях, которым свойственны большая трудоемкость, документальная невоспроизводимость, недоступность отдаленных и сложных по рельефу районов. Для главнейших древесных пород суммы площадей поперечных сечений деревьев в полных нормальных насаждениях установлены опытным путем и представлены массовыми стандартными таблицами [см. там же, Справочник, таблица 34 "Стандартные значения сумм площадей сечений нормальных древостоев основных лесообразующих пород по классам бонитета", стр.126-131 - аналог].

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

При наличии данных об абсолютной полноте определяется по соотношению сумм площадей сечений на 1 га таксируемого (ΣG_T) и нормального древостоя той же породы и типа леса (ΣG_H):

П=ΣG_T/ΣG_H,

где ΣG_H берутся из стандартных таблиц сумм площадей сечений.

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

Недостатками способа визуальной оценки являются:

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

- отсутствие инструментария, количественного измерения какой-либо функции и средств измерений этой функции;

- субъективность способа и большой диапазон ошибок.

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

Ближайшим аналогом к заявляемому техническому решению является "Способ вычисления запаса лесных массивов" Патент RU №2.242.867, A 01 G, 23/00; G 03 B, 37/00, 2004 г.

В способе ближайшего аналога получают изображение лесов, содержащих пробные площадки, в виде цифровой матрицы зависимости яркости I(х, у) от пространственных координат, рассчитывают пространственный спектр Фурье матрицы и находят среднюю частоту F_ср, диаметр кроны среднего дерева Д_ср=1/F_cp, площадь рельефа древесного полога S_p, в виде мозаики аппроксимирующих рельеф треугольников, вычисляют полноту насаждения П через отношение площади рельефа пробной площадки S_рп к удельной площади рельефа древесного полога, домноженного на полноту пробной площадки П_п

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

определяют прикрепляющую точку огивы анализируемого насаждения, вычисляют запас по массовым таблицам как сумму пяти классов ступеней толщины Лорея V=ΣN_i·H_i·g_i·f_i, где S₀ - площадь изображения обрабатываемого массива, га.

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

- площадь рельефа древесного полога (размерность м²) не может быть вычислена по безразмерной функции яркости I(х, у), которая обычно квантуется в шкале 0...256 уровней, без введения метрического коэффициента глубины (размерность м/шаг квантования);

- отсутствует эталонная функция зависимости относительной полноты древостоя от площади рельефа древесного полога: П=П(S_р).

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

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

определяют площадь рельефа древесного полога, получают эталонный ряд зависимостей известной относительной полноты пробных площадок от вычисленных значений площади их рельефов, определяют относительную полноту таксируемого массива по соотношению S_p/S_геом интерполированием смежных значений эталонного ряда, где:

S_геом, S_p - геометрическая площадь таксируемого участка и площадь рельефа его древесного полога;

Н - средняя высота древостоя участка, вычисляемая по регрессионным зависимостям (для деревьев первой и второй групп по Анучину Н.П.)

; .

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

фиг.1 - исходное изображение лесного массива (распечатка);

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

фиг.3 - зависимость среднеквадратического отклонения яркости σ от изрезанности рельефа древесного полога;

фиг.4 - изменение огибающей пространственного спектра в зависимости от изрезанности рельефа древесного полога;

фиг.5 - зависимость средней высоты древостоя (деревьев первой и второй групп по Анучину Н.П.) от пространственного спектра H=H(1/F_cp);

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

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

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

Определение полноты, как отношение сумм площадей сечений деревьев применимо лишь для перечислительной таксации. При глазомерной таксации найти соотношение сумм площадей сечений деревьев, определяющих полноту насаждения, невозможно. Сумма площадей сечений всех деревьев, образующих насаждение, составляет по отношению к занимаемой насаждением территории весьма малую величину 0,002...0,004. Такие величины на глаз практически неуловимы. В наиболее распространенных насаждениях горизонтальные проекции крон деревьев составляют 0,4...0,8 от занимаемой насаждением площади. Такие величины уже легче установить на глаз. Благодаря ежегодному приросту на десятки см побегов и ветвей, образующих крону деревьев, существенно изменяются как проекция крон, так и площадь рельефа древесного полога, образуемая совокупностью крон. Следует ожидать, что метод расчета относительной полноты древостоев путем отслеживания динамики изменения площади рельефа древесного полога, как это иллюстрируется фиг 2, обладает более высокой чувствительностью, чем известные аналоги.

Из математики известно [см., например, Н.С.Пискунов "Дифференциальное и интегральное исчисление для ВУЗов", учебник, 5^е изд., том.2, §7 Вычисление площади поверхности, стр.73-74, Наука, Москва, 1964 г.], что площадь рельефа вычисляют, как поверхностный интеграл функции h(x, у), заданной в области Ф:

Аналогом поверхностной функции h(x, у) является функция яркости изображения I(х, у). Поскольку функция яркости I(х, у) обычно квантуется в шкале 0...256 уровней квантования, то чтобы получить размерность площади S в [м²], функция яркости должна иметь размерность высоты рельефа, [м], т.е. шкалу яркости необходимо прокалибровать. Для этого вычисляют среднеквадратическое значение яркости σ. Зависимость среднеквадратического отклонения яркости от изрезанности рельефа древесного полога иллюстрируется фиг.3. Двойную величину (2σ) отождествляют с высотой рельефа древесного полога Н. Тогда калибровочный коэффициент глубины изображения Δh рассчитывают как: Δh=Н/2σ. Абсолютное значение высоты Н определяют в зависимости от диаметра кроны среднего дерева и расчетной функции пространственного спектра матрицы |m×n| изображения [см., например, С.В.Белов, Н.Д.Дмитриев, А.Е.Колосова "Аэрофотосъемка и авиация в лесном хозяйстве", Учебное пособие, Всесоюзный заочный ЛТИ, 1962 г., стр.145, табл.15, а также Анучин Н.П. "Лесная таксация", 1982 г., стр.250 и Патент RU №2.133.565, 1999 г.]. Эта зависимость для деревьев первой и второй групп имеет вид:

первая группа - сосна, лиственница, береза, осина, ольха ();

вторая группа - ель, пихта, кедр, ясень, бук, дуб ().

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

Следующей задачей является непосредственный расчет площади рельефа. Как следует из аналитического выражения поверхностного интеграла, поскольку частные производные ∂h/∂х, ∂h/∂у неизвестны, а сама функция задана матрицей |m×n| дискретных отсчетов, интеграл может быть рассчитан только численным методом.

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

Из математики известно [см., например, Г.Корн, Т.Корн "Справочник по математике для научных работников и инженеров", перевод с англ., Наука, Физматгиз, М., 1970 г., стр.144], что наибольшую точность вычислений обеспечивает аппроксимация функций бесконечным тригонометрическим рядом, коэффициенты разложения которого суть коэффициенты Фурье. В настоящее время существуют пакеты программ специализированного математического обеспечения, как то: ER MAPPER 5.0 ["Пакет программ для обработки изображений в науках о Земле", GENASYS Inc.San Diego, USA, p.283-294; MATHCAD 6.0 PLUS издание 2^е стереотипное, М.: Информ. издат. Дом ФИЛИНЪ, 1997 г., стр.412], реализующих алгоритмы быстрого Фурье - преобразования (БПФ) от цифровых массивов информации. Процедура вычислений реализуется следующей последовательностью операций. Массив I(х, у) вводят в ПЭВМ. Стандартной программой БПФ (см. выше ER MAPPER или MATHCAD) вычисляют двумерный пространственный спектр Фурье обрабатываемого участка (фрагмента):

Затем обратным Фурье-преобразованием представляют рельеф I(х, у) бесконечным тригонометрическим рядом:

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

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

Непосредственно, аналитически, данные интегральные суммы вычислить затруднительно. По своему физическому смыслу, если I(х, у) есть сигнал, то частные производные [∂I/∂_х]², [∂I/∂y]² (скорость в квадрате) есть скорость флуктуации мощности процесса в функции пространственных координат. Известно [см., например, Заездный A.M., "Основы расчетов по статистической радиотехнике", М., Сов. Радио, 1969 г., стр.91-94], что скорость флуктуации мощности в функции пространственных координат отражает автокорреляционная функция сигнала b(х, у). Значение автокорреляционной функции в нуле b(х=0, у=0) равно средней мощности процесса, а σ² - есть мощность переменной составляющей. После преобразования получено, что площадь рельефа древесного полога представляется функцией:

где σ²(H) - дисперсия высоты рельефа древесного полога (прокалиброванная функция глубины изображения).

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

Относительная полнота, П	1	0,9	0,8	0,7	0,6	0,5	0,4	0,3
Относительная площадь рельефа Sp/S0 геометрич.	1,1	1,25	1,4	1,56	1,72	1,9	2,1	2,3

Используя эталонный ряд, относительную полноту таксируемого массива определяют по расчетной величине S_р/S_{геометрич} интерполированием (редуцированием).

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

Заявляемый способ может быть реализован на базе устройства по схеме фиг.7. Функциональная схема устройства фиг.7 содержит орбитальный комплекс наблюдения 1, типа Международной космической станции (МКС) с установленной на ее борту поворотной платформой 2 (±15° от надира) и цифровой фотокамерой 3 (типа KODAK DCS760). Съемка запланированных участков лесов 4, включение фотокамеры 3 и выставка поворотной платформы 2 осуществляет бортовой комплекс управления 5 (БКУ) по командам, передаваемым из центра управления полетом 6 (ЦУП) по радиолинии управления 7. Информацию изображений лесных массивов 4 записывают на бортовой видеомагнитофон 8 (типа "Нива") и в сеансах видимости МКС с наземных пунктов сбрасывает по автономному каналу передачи данных 9 на наземные пункты приема информации 10 (ППИ), где осуществляют запись массивов информации на видеомагнитофон 11 (типа "Арктур").

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

Обрабатывалось изображение лесного массива Щелковского лесхоза [снимок системы "Ресурс", заказ Госцентра "Природа" №11/93-42], содержащего пробные площадки, таксационные характеристики которых ежегодно возобновляются студентами Лесного факультета МГУЛ при производственной практике. Изображение вводилось в ПЭВМ путем сканирования снимка сканером типа "Panasonic" с разрешением 1024 точки на дюйм, масштаб снимка соответствовал 1:550. Уровень яркости изображения I(х, у) находится в интервале I_max=152, I_min=12, среднее значение =46, среднеквадратическое отклонение яркости σ=62.

Стандартной процедурой специализированного программного обеспечения MATH САД вычислялся пространственный спектр участка изображения. Расчетные значения огибающих пространственных спектров различных участков (1, 2) иллюстрируются фиг.4. Средняя частота пространственного спектра первого участка составила (фиг.4) F_cp=0,32.

Откуда диаметр кроны среднего дерева Д_ср=1/F_cp=3,1 м.

Регрессионная зависимость высоты древостоя от диаметра кроны среднего дерева (средней частоты пространственного спектра) иллюстрируется графиками фиг.5. Пространственное разрешение обрабатываемого участка составило: Δх=Δу=1,4 м/пиксель, Δh=Н/2σ [м/шаг квантования]=0,25 м. После определения масштабных значений х, у, h расчет площади рельефа осуществляют численным методом по следующей программе:

Текст программы вычисления площади рельефа древесного полога.

program Podschet;

uses crt;

const maxx=500;

Pixelx=3470;

Pixely=2510;

Lx=29.38; {cm}

Ly=21.25; {cm)

mashtab=550; {in 1 sm 550 m}

Kmem=100;

var S:real;

ss:string;

Buf1, Buf2, Buf3: array [1..maxx] of real;

f, ff:text;

z:byte;

i, j, k, kolx, koly, kolM, zmax, zmin: integer;

dx, dy, dh, dpx, dpy, Ax, Ay, z1, mm, sigm, sigm2:real;

function Sq (i:integer):real;

function Geron (x1, x2, x3, y1, y2, y3, z1, z2, z3: real): real;

var p, a12, a23, a31: real;

begin

a12:=sqrt(sqr(x1-x2)+sqr(y1-y2)+sqr(z1-z2));

a23:=sqrt(sqr(x2-x3)+sqr(y2-y3)+sqr(z2-z3));

a31:=sqrt(sqr(x3-x1)+sqr(y3-y1)+sqr(z3-z1));

p:=(a12+a23+a31)/2;

Geron:=sqrt(p*(p-a12)*(p-a23)*(p-a31));

end;

begin

Sq:=Geron(0,dx,0,0,0,dy,Buf1[i]*dh,Buf1[i+1]*dh,Buf2[i]*dh)

+Geron(0,dx,dx,dy,0,dy,Buf2[i]*dh,Buf1[i+1]*dh,Buf2[i+1]*dh);

end;

begin {main};

clrscr kolM:=0;

repeat

writein('Введите имя файла');

readin(ss);

assign(f,ss);

{$I-} reset(f); {$I+}

i:=IOresult;

until i=0;

assign(ff,'Rez.txt');

rewrite(ff);

koly:=0;

mm:=0; zmax:=0; zmin:=255;

while not eof(f) do

begin koly:=koly+1;

kolx:=0;

while not eoln(f) do

begin

read(f,z); kolx:=kolx+1; mm:=mm+z;

if z<zmin then zmin:=z;

if z>zmax then zmax:=z;

end readln(f)

end;

close(f);

mm:=mm/kolx/koly;

writeln('Среднее значение z', mm:12:7);

writeln('Максимальное и минимальное значения', zmax:4, zmin:4);

writeln('kolx=',kolx:10);

writeln('koly=',koly:10);

writeln(ff,'Среднее значение z', mm:12:7);

writeln(ff,'Максимальное и минимальное значения', zmax:4, zmin:4);

writeln(ff, 'kolx=', kolx:10);

writeln(ff, 'koly=', koly:10);

if (kolx>maxx) or (koly>maxx) then

begin

writeln('Данные не умещаются в массив');

halt

end;

dx:=Lx/(Pixelx-1)*mashtab;

dy:=Ly/(Pixely-1)*mashtab;

Ax:=Lx/(Pixelx-1)*(kolx-1)*mashtab;

Ay:=Ly/(Pixely-1)*(koly-1)*mashtab;

reset(f);

sigm2:=0;

while not eof(f) do

begin

while not eoln(f) do

begin

read(f,z); sigm2:=sigm2+sqr(z-mm)/(kolx*koly)

end;

readln(f);

end;

close(f);

sigm:=sqrt(sigm2);

writeln('Дисперсия', sigm:12:7);

writeln(ff, 'Дисперсия', sigm:12:7);

dh:=26.0/(2*sigm);

dh:=1;

reset(f);

for i:=1 to kolx do

begin

read(f,z); z1:=z*dh; Buf1[i]:=z1;

end;

readln(f);

for i:=1 to kolx do

begin

read(f,z); z1:=z*dh; Buf2[i]:=z1;

end;

readln(f);

for i:=1 to kolx do

begin

read(f,z); z1:=z*dh; Buf3[i]:=z1;

end;

readln(f);

for i:=1 to kolx-1 do S:=S+sq(i);

while not eof(f) do

begin

Buf1:=Buf2; Buf2:=Buf3;

for i:=1 to kolx do

begin

read(f,z); z1:=z*dh; Buf3[i]:=z1;

end;

readln(f);

for i:=1 to kolx-1 do S:=S+sq(i);

end;

Buf1:=Buf2; Buf2:=Buf3;

for i:=1 to kolx-1 do S:=S+sq(i);

close(f);

writeln('Имя файла-',ss);

writeln('dx=',dx:10:5,' dy=',dy:10:5,' dh=',dh:10:5);

writeln('Строк-', koly:5,' Столбцов-', kolx);

writeln('Ширина области', Ax:15:7);

writeln('Длина области', Ау:15:7);

writeln('Проекция площади', Ах*Ау:15:3);

writeln('Площадь поверхности', S:15:3);

writeln('Отношение площадей', S/(Ах*Ау):10:7);

writeln(ff,'Имя файла-',ss);

writeln(ff,'dx=',dx:10:5, 'dy=', dy:10:5,' dh=', dh:10:5);

writeln(ff,'Строк-',koly:5,' Столбцов-',kolx);

writeln(ff,'Ширина области',Ax:15:7);

writeln(ff,'Длина области',Ау:15:7);

writeln(ff,'Проекция площади', Ax*Ay:15:3);

writeln(ff,'Площадь поверхности', S:15:3);

writeln(ff,'Отношение площадей', S/(Ax*Ay):10:7);

close(ff);

end.{main}

Расчетная площадь рельефа древесного полога таксируемого массива составила S_p=5091452,129 м²; геометрическая площадь участка S₀=2951895,366 м²; отношение S_p/S₀=1,75. Относительная полнота (интерполирование смежных значений эталонного ряда) П=0,58. Использование заявляемой технологии позволяет проводить измерения относительной полноты древостоев с точностью единиц процентов, а сам метод может рассматриваться как метрологический.

Способ определения полноты древостоев, при котором получают изображение лесов, содержащих пробные площадки в виде цифровой матрицы |m×n| зависимости функции яркости изображения - сигнала I(х, у) от пространственных координат путем съемки фотокамерой, установленной на орбитальном комплексе наблюдения, рассчитывают пространственный спектр Фурье матрицы изображения, находят среднюю частоту F_cp. и диаметр кроны среднего дерева D_cp.=1/F_cp., отличающийся тем, что последовательно обрабатывают участки изображения, вычисляют среднеквадратическое отклонение σ сигнала матрицы изображения обрабатываемого участка, калибруют шкалу яркости изображения для чего вводят калибровочный коэффициент шкалы функции яркости Δh[м/шаг квантования]=Н/2σ, определяют площадь рельефа древесного полога по формуле

где х, у - пространственные координаты;

σ²(Н) - прокалиброванная функция глубины изображения,

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

Н - средняя высота древостоя участка, вычисляемая по регрессионным зависимостям.