СПОСОБ ТАКСАЦИИ НАСАЖДЕНИЙ Российский патент 2002 года по МПК G03B37/00 A01G23/00 G06K9/36 

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

За период становления лесной науки разработано множество методов натурной таксации насаждений путем непосредственных измерений параметров деревьев на пробных площадках и последующего аналитического расчета таксационных характеристик по статистическим зависимостям (см. , например, Н.П. Анучин "Лесная таксация", учебник, 5-е изд., М., "Лесная промышленность", 1982 г., стр. 250, стр. 344-348 - аналоги).

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

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

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

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

Известен способ оценки запаса насаждений путем получения изображения лесного массива, преобразования его в дискретную матрицу отсчетов яркости, нахождения всех локальных максимумов матрицы, расчета среднего расстояния между максимумами, выделения контуров вокруг каждого максимума методами пространственного дифференцирования и оценки запаса по аналитическим зависимостям между рассчитанными числовыми характеристиками матрицы и элементами дерева среднего класса Лорея (см., например, патент РФ 2133565, кл. А 01 G 23/00, 23/02, 1999 г. - аналог). Способ-аналог может быть использован только при наличии высокодетальных снимков с разрешением менее 1 м, когда крона одиночного дерева отображается на снимке несколькими десятками пикселов. При меньшем числе пикселов (отсчетов) методы пространственного дифференцирования не работоспособны. Получить снимки с разрешением менее 1 м затруднительно, как по режимным ограничениям, так и по техническим возможностям космических систем.

Скрытые периодичности чередования крон древесного полога могут быть выявлены по снимкам среднего разрешения (более 1 м) путем вычисления спектра Фурье матрицы изображения.

Ближайшим аналогом по технической сущности с заявляемым способом является "Способ определения запаса насаждений" (патент РФ 2080051, кл. А 01 G 23/00, 1997 г.).

В способе ближайшего аналога получают изображение лесного массива, преобразуют функцию яркости изображения I(x, y) в матрицу |m×n| дискретных отсчетов, вычисляют характеристики электрического сигнала матрицы: огибающую пространственного спектра, математическое ожидание, среднеквадратическое отклонение σ, а запас насаждения рассчитывают по статистическим зависимостям между огибающей пространственного спектра, полнотой насаждения Р в функции от σ и параметрами отдельного дерева: высотой Н и сечением ствола g.

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

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

Решение поставленной задачи обеспечивается тем, что в способе таксации насаждений, при котором получают изображение лесного массива, преобразуют функцию яркости изображения I(x,y) в цифровую матрицу дискретных отсчетов размерностью |m×n| элементов, рассчитывают числовые характеристики матрицы и огибающую пространственного спектра, вычисляют таксационные параметры по статистическим зависимостям между числовыми характеристиками матрицы и элементами насаждений, рассчитывают математическое ожидание пространственного спектра (F_cp), осуществляют высокочастотную фильтрацию пространственного спектра с частотой среза, равной его математическому ожиданию, проводят обратное Фурье-преобразование отфильтрованного спектра, определяют полноту насаждения P как отношение количества выделенных в результате обратного Фурье-преобразования пикселей (n_c) к их общему количеству в матрице, находят прикрепляющую точку огивы насаждения к его изображению расчетом среднего количества деревьев в насаждении где D_cp=I/F_cp - средний размер крон древесного полога, S - площадь анализируемого участка насаждения, Н= 7•D_cp ^1,2 - средняя высота деревьев, а запас насаждения вычисляют по массовым таблицам, связывающим видовые числа f_i, высоту h_i, ступени толщины (d_1,3i), площади сечений и количество деревьев в насаждении N_i данной ступени толщины из соотношения:

Изобретение поясняется чертежами, где:
фиг.1 - исходное изображение анализируемого лесного массива;
фиг. 2 - визуализированный пространственный спектр Фурье сигнала матрицы изображения на комплексной плоскости;
фиг.3 - огибающая пространственного спектра Фурье:
а) распечатка с дисплея в логарифмическом масштабе программного обеспечения ПЭВМ,
б) реальная огибающая в масштабе исходного изображения;
фиг. 4 - маска оптимального высокочастотного фильтра на фоне спектра Фурье;
фиг. 5 - визуализированное изображение обратного Фурье-преобразования отфильтрованного спектра;
фиг. 6 - огива закономерности распределения деревьев в насаждении по естественным ступеням толщины;
фиг.7 - функциональная схема устройства, реализующего способ.

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

Техническая сущность изобретения заключается в следующем. Многие таксационные показатели, определяемые при оценке насаждений, не получают своего непосредственного отображения на снимках. Скрытые закономерности в строении насаждения и отдельных деревьев могут быть выявлены на основе использования нескольких независимых признаков изображения, таких как цвет, тон, текстура, топология и др. Чем больше независимых признаков изображения использовано при расчете, тем выше точность метода. Известны статистические зависимости между параметрами деревьев и элементами насаждения. Так, например, полнота насаждений связана с сомкнутостью полога, а диаметр крон - с высотой древостоя. Также установлено (см. , например, Н.П. Анучин "Лесная таксация", учебник, 5-е изд., М., "Лесная промышленность", с. 271-273), что распределение количества деревьев в насаждении по естественным ступеням толщины не зависит ни от породы, ни от бонитета, ни от полноты насаждения. Общей для всех насаждений и статистически устойчивой характеристикой распределения количества деревьев по естественным ступеням толщины является так называемая огива насаждения. Таким образом, исходными независимыми элементами насаждения являются: количество деревьев в насаждении N и распределение размеров крон деревьев по диаметрам D, м. Остальные таксационные параметры могут быть вычислены по аналитическим зависимостям, которые отражают статистические закономерности в виде функций регрессии этих параметров от независимых N и D. Для расчета независимых элементов N и D в заявляемом способе используют два независимых признака изображения: тон и текстуру. Текстура изображения содержит информацию о степени изрезанности древесного полога, образуемой кронами деревьев. Крона отдельного дерева при съемке сверху представляется некоторой колоколообразной фигурой. Вершина кроны освещена ярче, чем промежутки между кронами деревьев. Поэтому, текстура изображения повторяет периодичности чередования крон деревьев. Скрытая информация о периодичности чередования и размерах крон деревьев может быть выявлена путем вычисления двумерного спектра Фурье от анализируемой матрицы изображения (см., например, "Анализ пространственных частот" в книге Р. Дуда, П. Харт "Распознавание образов и анализ сцен", перевод с англ., Мир, Москва, 1976 г. , стр. 319-321).

Двумерный пространственный спектр L(F_x, F_y) от функции яркости изображения I(x, y) рассчитывают преобразованием Фурье:

Представленное аналитическое выражение является стандартной математической операцией, реализуемой алгоритмами быстрого Фурье-преобразования (БПФ), входящих в комплект специализированного программного обеспечения: "Пакет программ для обработки изображений в науках о Земле" (см., например, "ER MAPPER Reference", Earth Resource Mapping Pty Lid, Western Australia, 6005, 1998 г., р. 295). На фиг.1 представлено исходное изображение лесного массива, масштаб снимка М 1:10000. На фиг.2 представлен визуализированный двумерный пространственный спектр Фурье на комплексной плоскости после программной реализации операции двойного интегрирования I(x, y). На фиг.3 представлена огибающая пространственного спектра: а) распечатка с дисплея в логарифмическом масштабе программного обеспечения ER MAPPER, б) реальная функция в масштабе исходного изображения. Полученная функция по физическому смыслу представляет собой величину, обратную линейным размерам объектов подстилающей поверхности, отождествляемых с распределением диаметров крон деревьев D. Расчет огибающей реализуется стандартной программой, входящей в комплект программного обеспечения ER MAPPER 5,0 (см. там же, стр. 58, 286). Амплитуда А огибающей в каждой точке F_i определяет удельный вес соответствующей пространственной гармоники в общем спектре. Математическое ожидание (среднюю частоту) пространственного спектра находят из условия:

При этом расчетный средний диаметр крон анализируемого участка насаждения соответствует:

Другим независимым признаком изображения является тон. При преобразовании изображения в цифровую матрицу дискретных отсчетов функция яркости I(x, y) обычно квантуется в стандартной шкале от 0 до 256 уровней. Для расчета количества деревьев необходимо превратить исходное изображение в контурный рисунок, т.е. выделить контуры всех крон. Контур - это край участка изображения, где велик градиент изменения функции яркости I(x, y). Если разрешение снимка невелико и число пикселов недостаточно для расчета градиента, то эффективным методом выделения контурного рисунка является высокочастотная фильтрация (см., например, "Пространственная фильтрация" в книге Р. Дуда, П. Харт "Распознавание образов и анализ сцен", с. 329).

Чтобы эффективно выделить контуры крон на анализируемом изображении, необходимо оптимальным образом рассчитать передаточную характеристику К(F_x,F_y) высокочастотного фильтра. Неидентифицируемые разрывы между кронами, протяженные прогалины, опушки, соответствуют низким частотам пространственного спектра изображения и расположены на комплексной плоскости в начале координат. Следовательно, для повышения точности оценок необходимо срезать низкие частоты (НЧ), и подчеркнуть высокие частоты (ВЧ). Границей раздела (частотой среза) между НЧ и ВЧ служит рассчитанная выше характеристика - математическое ожидание спектра F_cp.

На фиг.4 представлена визуализированная маска высокочастотного фильтра с частотой среза F_cp. Высокочастотная фильтрация заключается в перемножении исходного пространственного спектра L(F_x,F_y) на передаточную функцию высокочастотного фильтра К(F_x, F_y), в результате чего получают отфильтрованный спектр L₀(F_x,F_y):
L₀(F_x,F_y)=L(F_x,F_y)•K(F_x,F_y).

Высокочастотная цифровая фильтрация реализуется программным расчетом на ПЭВМ с использованием стандартной процедуры (см., там же ER MAPPER 5,0 стр. 299).

После оптимальной фильтрации сигнала матрицы изображения осуществляют восстановление функции яркости изображения I_о(x, y) обратным Фурье-преобразованием:
I₀(x,y) = Ф^-1[L_o(F_x,F_y)].
Данная операция также реализуется стандартной программной процедурой, см. ER MAPPER 5,0, стр. 299.

Визуализированное изображение после операции обратного Фурье-преобразования представлено распечаткой с ПЭВМ на фиг.5. Как следует из фиг.5, визуально, операция ВЧ фильтрации состоит в ранжировании пикселей по яркости. Светлые участки изображения соответствуют площади, занимаемой кронами деревьев, темные участки - прогалинам. Отношение числа пикселей n_с, принадлежащих к светлым участкам, к общему числу пикселей |m×n| представляет собой характеристику насаждения - проекцию древесного полога на поверхность Земли или полноту (Р):

Располагая численным значением полноты, определяют среднее количество деревьев N_cp на анализируемом изображении через отношение площади участка S, домноженной на полноту P, к площади кроны среднего дерева:

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

Зависимость высоты дерева Н от диаметра кроны D представляется степенной функцией вида (см. аналог, патент РФ 2133565, 1999 г.)
H,м≅7•D^1,2(м).
Между высотой дерева Н и его толщиной d_1,3 существует зависимость, представленная в массовых таблицах А.В. Тюрина, Д.И. Товстолеса, В.К. Захарова (см. , например, "Общесоюзные нормативы для таксации лесов", Справочник, М. , "Колос", 1992 г., стр. 67-73). Кроме того, существуют таблицы всеобщих видовых чисел (f) М.Е. Ткаченко (см., например, Н.П. Анучин "Лесная таксация", 5-е издание, М., "Лесная промышленность", 1982 г., стр. 148).

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

где f_i - видовое число, g_i - площади поперечных сечений,
N_i, h_i - количество и высота деревьев i-й ступени толщины на анализируемом изображении.

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

Заявляемый способ может быть реализован на базе устройства по схеме фиг. 7. Функциональная схема устройства фиг.7 содержит орбитальный комплекс наблюдения 1 типа космического аппарата "Ресурс-Ф" с установленной на его борту фотокамерой 2 типа КФА-1000, осуществляющей съемку участков 3 подстилающей поверхности 4 по командам из Центра управления полетом 5, передаваемым по радиолинии управления 6. Заполненные кассеты с фотоматериалами съемки отстреливаются с борта космического аппарата по командам из ЦУПа и в виде капсул 7 доставляются в лабораторию обработки фотоматериалов 8. При обработке фотоматериалов осуществляется их предварительная геометрическая и радиометрическая коррекция. Материалы космической съемки в виде откорректированной серии снимков 9 поступают заказчику, в Центр тематической обработки 10. Посредством сканера 11 высокого разрешения типа "Panasonic" аналоговые снимки преобразуются в цифровые матрицы дискретных отсчетов и вводятся в ПЭВМ 12, включающей стандартный набор элементов в составе: процессора 13, оперативного запоминающего устройства 14, винчестера 15, клавиатуры 16, дисплея 17 и принтера 18. В качестве ПЭВМ 12 тематической обработки снимков 9 использовалась станция CVN со специальным программным обеспечением ER MAPPER 5.0, предварительно записанным на винчестер 15. Любое специализированное программное обеспечение в качестве единицы расстояния на плоскости изображения оперирует размерами одного элемента изображения - пикселя. Единицей пространственной частоты при этом считается величина, обратная элементу изображения
Следовательно, чтобы оперировать реальными линейными размерами объектов на изображении, необходимо связать масштаб снимка с пространственным разрешением одного пикселя.

Масштаб исходного снимка соответствовал M 1:10000. Для ввода в ПЭВМ снимок был отсканирован с разрешением 400 точек на дюйм. Стандартная матрица изображения имела размерность 512•512 отсчетов. Пространственное разрешение одного пикселя составило 0,64 м, а площадь анализируемого изображения S = |330•330|, м ≈11 га. Программное обеспечение ER MAPPER 5,0 выдает функцию огибающей пространственного спектра в логарифмическом, по основанию 2, масштабе. На фиг. 3а иллюстрируется распечатка этой функции с экрана дисплея. На фиг.3б представлена огибающая пространственного спектра, пересчитанная из логарифмического масштаба в масштаб реального изображения. Математическое ожидание пространственного спектра исходного снимка составило F_ср= 0,26, а средний диаметр крон древесного полога D_cp=I/F_cp=3,8 м. Средняя высота Н=7•D_cp ^1,2=34 м.

После высокочастотной фильтрации с частотой F_cp=0,26 и обратного Фурье-преобразования, произошло ранжирование пикселей изображения по яркости. Восстановленное после ВЧ фильтрации изображение иллюстрировано на фиг.5. Численный расчет количества светлых пикселей (n_c) к общему их числу в матрице дает значение полноты:

Количество деревьев в прикрепляющей точке огивы насаждения:

что составляет 66,3% от общего количества деревьев на 1 га.

Расчетный диаметр ствола среднего дерева в прикрепляющей точке огивы насаждения равен d_1,3=46 см. По массовым таблицам находят диаметры стволов и видовые числа деревьев насаждения всех ступеней толщины, а по огиве - их количество N_i в каждой ступени. Запас насаждения рассчитывают численными методами с использованием массовых таблиц. Расчетный запас насаждения для рассмотренного примера реализации составил
V = Σg_i•h_i•f_i•N_i = 320 м³/га.
Эффективность заявляемого способа оценивалась путем сравнения полученных таксационных характеристик насаждения с результатами наземной перечислительной таксации. Перечислительная таксация была выполнена на момент получения исходного изображения. Для анализируемого участка таксационные характеристики натурной таксации соответствовали: Р=0,4, d_1,3=48 см, Н=33 м, V=340 м³/га, что не выходит за пределы допуска (см., например, "Точность определения таксационных показателей" в справочнике "Общесоюзные нормативы для таксации лесов", М., "Колос", 1992 г., с. 121).

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

Изобретение относится к лесному хозяйству, методам дистанционного решения задач лесохозяйственного назначения. При реализации заявленного способа получают изображение лесного массива, преобразуют функцию яркости изображения I(x,y) в цифровую матрицу дискретных отсчетов размерностью |m×n| элементов, рассчитывают числовые характеристики матрицы и огибающую пространственного спектра, вычисляют таксационные параметры по статистическим зависимостям между числовыми характеристиками матрицы и элементами насаждений. Способ позволяет вычислить несколько независимых элементов насаждений по нескольким независимым признакам изображения. Технический результат - повышение точности и устойчивости результатов оценок расчетных параметров. 7 ил., 1 табл.

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

где D_ср= I/F_ср- средний размер крон древесного полога, S - площадь анализируемого участка насаждения, Н= 7•D_cр ^1,2 - средняя высота деревьев, а запас насаждения вычисляют по массовым таблицам, связывающим видовые числа f_i, высоту h_i, ступени толщины (d_1,3i), площади сечений и количество деревьев в насаждении N_i данной ступени толщины из соотношения V = Σq_i•h_i•f_i•N_i.