СПОСОБ ОПРЕДЕЛЕНИЯ ПРИРОСТА ЗАПАСА НАСАЖДЕНИЙ Российский патент 2006 года по МПК A01G23/00 

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

Текущий прирост по запасу является одним из главных показателей продуктивности древостоев, позволяющий на ограниченном по объему опытном материале, используя типовые и стандартизованные модели, составлять в интересах хозяйственного использования региональные таблицы хода роста [см., например, Справочник "Общесоюзные нормативы для таксации лесов", М., Колос, стр.299-338, табл.130-140, стр.346].

Известны методы определения текущего прироста путем повторных измерений на постоянных или временных пробных площадках [см. там же "Справочник", стр.355-350 - аналог]. В способе-аналоге подбирают модельные деревья по ступеням толщины, по двухметровым отрезкам срубленных модельных деревьев, замеряют прирост, по диаметру аппроксимируют модельные деревья таксационными показателями, близкими к средним по высоте, диаметру, форме ствола, а процент текущего прироста вычисляют по формуле

где M₁, М₂...M_n - запасы отдельных ступеней толщины;

P₁, Р₂...Р_n - проценты текущего прироста модельных деревьев по ступеням толщины;

М - запас растущей части древостоя.

Способ трудоемок, требует длительного периода наблюдений, не обеспечивает достаточную точность определения текущего прироста.

Для решения практических задач, не требующих высокой точности, используют упрощенные методы расчетов. Известны расчетные формулы Н.П.Анучина, В.В.Загреева, М.Л.Дворецкого и др., а также таблицы текущего прироста древостоев по запасу. Во всех формулах предлагается на одних и тех же деревьях, отобранных методом случайной выборки, проводить измерение средних таксационных характеристик различных элементов леса (не менее 0,2 состава) [см., там же, Справочник, «Общесоюзные нормы для таксации лесов», стр.356-357, табл.144 - аналоги]. Недостатками аналогов являются большая трудоемкость и неоперативность работ, ограниченный объем выборки измерений.

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

Массив пикселей содержит информацию как о параметрах отдельных деревьев, так и о статических характеристиках элементов леса в целом. Ближайшим по технической сущности является «Способ таксации насаждений», патент RU №2183847, А 01 G 23/00, 2002 г.

В способе ближайшего аналога получают изображение лесного массива, преобразуют функцию яркости изображения I (х, у) в цифровую матрицу дискретных отчетов размерностью |m×n| элементов, вычисляют огибающую пространственного спектра и его математическое ожидание (F_cp), осуществляют высокочастотную фильтрацию пространственного спектра с частотой среза, равной его математическому ожиданию, проводят обратное Фурье-преобразование отфильтрованного спектра, определяют полноту насаждения (П) как отношение количества выделенных в результате обратного Фурье-преобразования пикселей (n_c) к общему количеству в матрице, находят прикрепляющую точку огивы насаждения к его изображению расчетом среднего количества деревьев в насаждении

где D_cp - средний диаметр крон древесного полога, S - площадь анализируемого участка насаждения, - средняя высота деревьев, а запас вычисляют по массовым таблицам, связывающим видовые числа f_i, высоту h₁, ступени толщины (d_1,3i) площади сечений и количество деревьев в насаждении N_i данной ступени толщины из соотношения

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

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

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

Задача, решаемая заявляемым способом, состоит в повышении точности, оперативности, статической устойчивости результата путем определения нескольких интегральных характеристик древостоя по его изображению. Поставленная задача решается тем, что в способе определения прироста запаса насаждений, включающем получение изображений лесных массивов в виде зависимости яркости I (х, у) от пространственных координат, преобразование функции яркости изображения в цифровую матрицу дискретных отчетов размерностью |m×n| элементов, программный расчет пространственного спектра изображения анализируемой матрицы и определение прикрепляющей точки D_cp=1/F_cp огивы насаждения анализируемого участка, расчет запаса насаждений по числовым характеристикам анализируемой матрицы и массовым таблицам, дополнительно снимки лесных массивов, содержащих пробные площадки, осуществляют на временном лаге в несколько лет в одни и те же вегетационные периоды года, программным расчетом матрицы вычисляют площадь рельефа древесного полога S_p анализируемого участка, а прирост запаса насаждений Р_м рассчитывают по формуле

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

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

D_cp - средний диаметр кроны древесного полога;

F_cp - средняя частота пространственного спектра матриц изображений.

Изобретение поясняется чертежами, где фиг.1 - исходное изображение участка лесного массива (распечатка с ПЭВМ);

фиг.2 - сечение древесного полога вертикальной плоскости;

а) рельеф исходного и б) прирастающего полога;

фиг.3 - семейство огибающих пространственного спектра матриц с различными значениями Fcp;

фиг.4 - огива насаждения и ее прикрепляющая точка;

фиг.5 - вид автокорреляционной функции изображения: а) одномерной б) двухмерной автокорреляционной функции в аксонометрии;

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

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

Изображение лесного массива представляется функцией дискретных отчетов яркости изображения I (х, у) от пространственных координат. Падающий световой поток по разному отражается от изрезанного древесного полога (фиг.2). Вершины крон деревьев освещены лучше, в то время как в промежутках между кронами отражение носит диффузный характер. Следовательно, текстура изображения содержит скрытую информацию о геометрии древостоя, его пространственной структуре, т.е. закономерность чередования светлых и темных пикселей изображения зависит от рельефа древесного полога. За счет ежегодного роста ветвей и молодых побегов рельеф полога изменяется быстрей, чем толщина стволов. Следует ожидать, что метод оценки прироста по динамике изменения диаметров крон и площади рельефа обладает более высокой чувствительностью, чем способы-аналоги. Задача состоит в том, чтобы адекватно связать интегральные показатели матрицы изображения как площадь рельефа полога и огибающую пространственного спектра со средними таксационными параметрами древостоя. Известно, что распределение количества деревьев в насаждении по естественным ступеням толщины не зависит ни от типа породы, ни от полноты, ни от бонитета и характеризуется так называемой огивой насаждения (см., например, Н.П.Анучин «Лесная таксация», учебник, 5-е изд., М.: Лесная промышленность, с.271-273). Прикрепляющая точка огивы насаждения делит ступени толщины деревьев в насаждении в пропорции 60% тонких и 40% толстых . Чтобы вычислить прирост Р_м в процентах, нет необходимости рассчитывать детальные таксационные характеристики. Для статической устойчивости результата оценок следует оперировать средними показателями. Из огивы насаждения и аналитических соотношений способа-аналога следует, что текущий прирост можно вычислить как

где D_1cp, D_2cp - средние диаметры крон исходного и прирастающего насаждения, равны соотношению D_1cp=1/F_cp1, D_2cp=1/F_cp2;

N_1cp, N_2ср - среднее количество деревьев данной ступени толщины в насаждении.

Средний диаметр крон анализируемого участка рассчитывают путем вычисления спектра Фурье матрицы |m×n| функции яркости I (х, у) стандартной процедурой специализированного программного обеспечения (см. Пакет программ для обработки изображений в науках о Земле, «ER MAPPER Reference», Earth Resource Mapping Pty Lid, Western Australia, 6005, 1998, p.295, см. также способ ближайшего аналога патент RU №2183847, 2002 г.). На фиг.3 представлено семейство огибающих пространственного спектра а) исходного, б) прирастающего насаждения со средними значениями пространственных частот F_cp1 и F_cp2. Полученная функция по физическому смыслу представляет собой величину, обратную линейным размерам объектов подстилающей поверхности, отождествляемой с распределением крон деревьев по диаметру Д. Второй сомножитель (N_cp) формулы процента прироста свяжем с площадью рельефа (S_p) древесного полога.

По определению полнота древостоя характеризует степень использования им занимаемого пространства. Относительная полнота, выражаемая в долях от полноты (П) нормального древостоя, определяется по сомкнутости полога и плотности стояния деревьев [см., например, Справочник, «Общесоюзные нормативы для таксации лесов», М.: Колос, стр.116]. Нормальным, имеющим полноту П=1 считается участок, в котором нет ни одного недостающего и ни одного лишнего дерева. При известной полноте пробной площадки (П_п) полноту участка определяют через отношение суммы площадей проекции крон к линейной площади таксируемого насаждения

или через отношение площадей рельефа их пологов: анализируемого участка и пробной площадки S_pn

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

При известном аналитическом выражении функции I (х, у) площадь рельефа вычисляют как поверхностной интеграл от этой функции [см., например, Н.С.Пискунов «Дифференциальное и интегральное исчисление для втузов», учебник, 5-е изд., том 2. М.: Наука, 1964 г., стр.73-74]

Поскольку сама функция I (х, у) задана массивом дискретных отсчетов, для вычисления двойного интеграла необходимо иметь аналитическое выражение аппроксимирующей функции. Наибольшую точность вычислений обеспечивает аппроксимация функций бесконечным тригонометрическим рядом, коэффициенты разложения которого суть коэффициенты Фурье (см, например, Г.Корн, Т.Корн, «Справочник по математике для научных работников и инженеров» перев. с англ., М.: Наука, Физматгиз, 1970, стр.144]. Процедура вычисления площади рельефа древесного полога реализуется следующей последовательностью операций. Изображение заданного участка посредством сканера вводят в ПЭВМ. Стандартной программой быстрого Фурье образования (БПФ) вычисляют двухмерный пространственный спектр Фурье обрабатываемого участка

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

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

Непосредственно аналитически данные, интегральные суммы вычислить затруднительно. По физическому смыслу, если I (х, у) есть амплитуда сигнала, то частные производные есть скорость флюктуации мощности процесса по пространственным координатам. Известно [см., например, А.М.Заездный «Основы расчетов по статической радиотехнике», М.: Сов. Радио, 1969 г., стр.91-94], что скорость флюктуации мощности в функции пространственных координат отражает автокорреляционная функция сигнала b(х, у). Автокорреляционная функция сигнала вычисляется как обратное Фурье-преобразование от энергетического спектра стандартной программой процедурой специализированного программного обеспечения [см., там же ER MAPPER]. Представляя пространственный спектр канонической функцией вида получено, что огибающие автокорреляционной функции, характеризующие скорость флюктуации мощности процесса, имеют вид

На фиг.5а представлен вид автокорреляционной функции b(х) по одной из координат. Значение автокорреляционной функции в нуле равно средней мощности процесса, т.е. постоянной составляющей плюс переменной составляющей (дисперсии) σ². На фиг.5б представлен вид двухмерной функции b(х, у) в аксонометрии.

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

Приведенные выше функциональные соотношения реализуются следующей специализированной программой для ПЭВМ.

Текст программы

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

Заявляемый способ может быть реализован на базе устройства по схеме фиг.6. Функциональная схема устройства фиг.6 содержит орбитальный комплекс наблюдения 1 типа космического аппарата «Ресурс-Ф» с установленной на его борту фотокамерой 2 типа КФА-1.000, осуществляющей съемку участков 3 подстилающей поверхности 4 по командам из Центра управления полетом 5, передаваемым по радиолинии управления 6. Заполненные кассеты с фотоматериалом съемки отстреливается с борта космического аппарата по командам из ЦУПа и в виде капсул 7 доставляется в лабораторию обработки фотоматериалов 8 Федерального Научного Центра «Природа». При обработке фотоматериалов осуществляется их предварительная геометрическая и фотометрическая коррекция. Материалы космической съемки в виде откорректированной серии снимков 9 поступают заказчику в Центр тематической обработки 10. Посредством сканера 11 высокого разрешения типа Panasonic аналоговые снимки преобразуются в цифровые матрицы |m×n| дискретных отсчетов I (х, у) и вводятся в ПЭВМ 12, включающей набор следующих элементов: процессор 13, оперативное запоминающее устройство 14, винчестер 15, дисплей 16, принтер 17, клавиатура 18. В качестве ПЭВМ 12 тематической обработки снимков использовалась станция SYN со специальным программным обеспечением ER MAPPER 5.0, предварительно записанным на винчестер 15. Любое специализированное программное обеспечение в качестве единицы расстояния на плоскости изображения оперирует размерами одного элемента изображения - пикселя. Единицей пространственной частоты при этом считается величина, обратная элементу изображения F=1/элемент.

Следовательно, чтобы оперировать реальными линейными размерами объектов на изображении, необходимо связать масштаб снимка с пространственным разрешением одного пикселя. При обработке использован снимок Щелковского лесхоза масштаба 1:10000. Для ввода в ПЭВМ снимок был отсканирован с разрешением 600 точек на дюйм, пространственное разрешение пикселя ˜0,50 м.

Стандартная матрица изображения |m×n| составила |512×512| элементов. Площадь анализируемого участка S=|256×256| м² или ˜6,55 га. Математическое ожидание частоты пространственного спектра составило F_cp1=0,28, а средний диаметр D_1cp=3,5 м. Уровень яркости изображения находился в интервале I_max=135 I_min=23, среднее значение яркости 64, дисперсия по дисперсия по у Площадь рельефа древесного полога S_p1=32,2 га.

Повторная съемка проведена через 5 лет. Средняя частота пространственного спектра составила F_cp2=0,26, D_cp1=3,8, а площадь рельефа того же участка S_p2=26,3 га. Оттуда прирост по запасу составил

или 5,6% в год.

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

Способ может быть использован в лесном хозяйстве для дистанционного решения лесохозяйственных задач на основе программной обработки изображений лесов. Способ заключается в том, что получают изображение лесных массивов в виде зависимости функции яркости изображения I (х, у) от пространственных координат. Преобразуют функцию яркости изображения в цифровую матрицу дискретных отсчетов размерностью |m×n| элементов. Осуществляют программный расчет пространственного спектра изображения анализируемой матрицы и определяют прикрепляющую точку D_cp=1/F_cp огивы насаждения анализируемого участка. Осуществляют расчет запаса насаждений по числовым характеристикам анализируемой матрицы и массовым таблицам. Изображения лесных массивов, содержащих пробные площадки, осуществляют на временном лаге в несколько лет в одни и те же вегетационные периоды года. Программным расчетом матрицы вычисляют площадь рельефа древесного полога S_p анализируемого участка, а прирост запаса насаждений Р_м рассчитывают по формуле

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

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

D_cp - средний диаметр кроны древесного полога;

F_cp - средняя частота пространственного спектра матриц изображений.

Это позволит повысить точность, оперативность, достоверность, статистическую устойчивость результата. 6 ил.

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

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

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

D_cp - средний диаметр кроны древесного полога;

F_cp - средняя частота пространственного спектра матриц изображений.