Картографирование древостоев с использованием
космических снимков

Работа профинансирована "Программой фундаментальных исследований Республики Коми", 2000 г.

Введение
Объекты исследования 
Методика работы
Подготовка данных
Построение регрессии 
Расчет регрессии для темнохвойных пород
Проверка результатов расчета регрессии для темнохвойных
Расчет регрессии для светлохвойных пород 
Расчет регрессии для мелколиственных пород 
Проверка регрессии по формуле (4)
Заключение 
Список литературы 

Введение

В прошлые годы, использование вычислительной техники позволило значительно увеличить эффективность работы органов государственного управления за счет автоматизации документооборота и автоматизированной обработки информации. В то же время, такие мощные средства автоматизации как использование систем информационного обеспечения принятия решений практически не применялись из-за сравнительной дороговизны исходной информации для них (цифровые модели местности и их «семантическое», смысловое наполнение) и необходимости ее быстрого и независимого получения в достаточно больших объемах.

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

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

 Решению этих проблем могут способствовать три фактора появившихся в последнее время. Во первых, это широкое внедрение геоинформационных систем (ГИС), позволяющих проводить сложный анализ с использованием цифровой картографической информации и, во вторых, возросшее качество и доступность данных дистанционного зондирования Земли (ДДЗЗ) или космических снимков, которые могут быть с успехом использованы для оперативного и независимого получения обобщенной информации о территории. Третьим фактором, является возможность современных ГИС работать c данными дистанционного зондирования, проводить их обработку и анализ совместно с картографической и семантической (смысловой) информацией. Таким образом, совместное использование ГИС и ДДЗЗ должно позволить решить многие из поставленных задач.

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

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

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

Исходя из статистической значимости построенных уравнений регрессии появляется возможность определить степень пригодности отечественных космических природоресурсных систем, основанных на использовании камеры МСУ-Э (спутники «Ресурс-О1», «Океан», модуль «Природа») для решения задач упрощенной детальной инвентаризации лесов по ряду параметров. 

В качестве «наземных наблюдений» были использованы данные государственного лесоустройства. В практике работ по дешифрированию космических снимков спутника серии РЕСУРС-О1 данные лесоустройства используются впервые, поэтому одной из целей данной работы было выяснение приделов их пригодности и задач, в которых материалы лесной таксации способны полностью или частично «подменить» данные наземных наблюдений.

Вычисления проводились с использованием геоинформационных систем ARC\INFO 7.2.1.и ArcView 3.2, а так же пакета для статистического анализа SPSS 8.0.

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

Объекты исследования

В качестве объектов исследования выбраны насаждения Кедвинского лесничества Ухтинского лесхоза. Лесоустройство проведено в 1996 г. (картографические материалы изданы в 1997 г.) Северным лесоустроительным предприятием, расположенным в г. Вологда. Данные лесоустройства представлены в виде цифровой карты базового масштаба 1:25000, в условной системе координат.

Для каждого таксационного выдела, показанного на цифровой повыдельной карте в атрибутивной таблице представлены таксационные данные в формате кодировки, принятом для АСУ-ЛР [1] . Наиболее важные для дальнейшей работы виды таксационных характеристик приведены в таблице 1.

Таблица 1.  Наиболее важные атрибутивные данные, описывающие лесоустроительные выдела, которые были использованы в работе.

Лесоустроительная цифровая карта была преобразована в систему координат Гаусса-Крюгера 10-ой зоны, для ее совмещения с космическими снимками.

Лесоустроительная карта была разделена на два не перекрывающихся участка, один из которых был использован в качестве основного (далее – Основной участок или участок №1) и содержал 748 лесных выделов, другой использовался в качестве дополнительного участка (далее – Дополнительный участок или участок №2) для уточнения отдельных положений и выводов и содержал 5470 выделов. Карта схема расположения участков приведена в Приложении 1.

Участок №1 был таксирован методом инвентаризации лесов при повторном лесоустройстве (III разряд). Участок №2 таксирован в основном методом инвентаризации лесов при повторном лесоустройстве (III разряд), а частично методом камерального дешифрирования аэрофотоснимков (III разряд). [2]

Для проверки корректности построения регрессии были использованы цифровые лесоустроительные карты Сюзъюнского лесничества. Участок, по которому проводилась проверка регрессии находится в западной части лесничества и включает в себя 4500 выделов. Участок таксирован методом инвентаризации лесов при повторном лесоустройстве (III разряд).

Методика работы

Рисунок 1. Спектр растительности. Пояснения к рисунку даны в тексте.

Можно видеть, что полоса поглощения хлорофилла (зона 2) попадает в 3-ий спектральный канал камеры МСУ-Э, а Зона 1 во 2-ой канал.

Действительная ситуация, при работе с данными дистанционного зондирования значительно сложнее. Дело в том, что из-за низкой разрешающей способности камеры МСУ-Э (около 40 метров), реально оказывается наблюдаемым не спектр отдельных деревьев, а интегральный спектр некоторого участка местности, сравнимого по размерам с разрешающей способностью снимка. Таким образом реально возможно изучение только этого «интегрального» спектра целого растительного сообщества. «Интегральный» спектр растительности в целом обладает теми же характеристиками, что и спектр отдельного дерева. Следовательно, при разработке метода автоматического дешифрирования мы можем пользоваться теми же характерными особенностями спектра.

Введем следующие основные обозначения:

COUNT10, COUNT20, COUNT200 – число в лесной формуле, соответствующее содержанию сосны, ели (ели, кедра и лиственницы) и мелколиственных пород,

RELECT10, REFLECT20, REFLECT200 – значения спектральных яркостей измеренные по космическому снимку в 1-ом, 2-ом и 3-ем каналах соответственно.

Решение поставленной задачи состоит в том, чтобы установить математико-статистические зависимости линейной регрессии вида:

COUNTn = K1n* REFLECT1 + K2n * REFLECT2 + K3n * REFLECT3 (1)

При этом содержание древесных пород каждого из 3-х видов определяется независимо от содержания другого вида. Это уравнение определяет некоторую «плоскость» в пространстве факторов (спектральных каналов) и носит название «спектральный индекс».

Зная уравнения вида (1), построенные на основе изучения данных опытного участка для каждого типа древесных пород, в дальнейшем можно применять их к калиброванным космическим снимкам на другие участки местности и, тем самым строить тематические карты (в GRID виде) [3]

Для выполнения работы были использованы два снимка полученные спутником Ресурс-О1 №3, (камера МСУ-Э) в 1996 и 1998 годах. Сведения о снимках приведены в таблице 2.

Таблица 2. Данные о космических снимках.

Космические снимки были трансформированы в систему координат Гаусса-Крюгера, 10 зона, по трансформационным точкам, полученным с цифровой карты масштаба 1:200000. Достигнута удовлетворительная степень совмещения по твердым контурам (дорогам, просекам) с лесоустроительной картой. Снимок №1 оттрансформирован несколько хуже, чем снимок №2. В первом случае ошибки несовпадения контуров достигают 50-80 метров, во втором случае 30-60 метров.

Подготовка данных

Среди основных древесных пород, произрастающих на территории Республики Коми наибольшее распространение и хозяйственное значение имеют: ель, сосна и береза. На территории опытных участков так же в незначительных количествах произрастают: кедр (сосна кедровая), пихта, лиственница. По данным лесоустройства, примеси указанных пород до 20% от общего состава по запасу имеются только на 5% всех участков. Имеются 6 (менее 1%) участков, на которых в качестве основной породы указана лиственница. В виду того, что указанные породы не имеют важного хозяйственного значения, было принято решение о группировке всех встречающихся видов пород по их спектральным характеристикам на: темнохвойные, светлохвойные и мелколиственные породы.

В группу темнохвойных пород входят: ель и кедр; в группу светлохвойных: сосна и лиственница; в группу мелколиственных: береза, осина и ивы. Регрессионные зависимости и определялись для указанных групп.

Исходными данными для решения поставленных задач является набор отношений вида:

NV, WLK1, KPLK1, … WLK5, KPLK5, COUNT10, COUNT20, COUNT200, REFLECT1, REFLECT2, REFLECT3, D1, D2, D3  (2)

Где:

1. NV – номер выдела
2. WLK1, KPLK1, … WLK5, KPLK5 – коэффициенты лесных формул для 1-ого яруса и типы древесных пород
3. COUNT10 - суммарный коэффициент для сосны и лиственницы (светлохвойные)
4. COUNT20 – суммарный коэффициент для ели и кедра (темнохвойные)
5. COUNT200 – суммарный коэффициент для всех мелколиственных пород
6. REFLECT1 – среднее значение спектральной яркости для 1-ого спектрального канала
7. REFLECT2 - среднее значение спектральной яркости для 2-ого спектрального канала
8. REFLECT3 - среднее значение спектральной яркости для 3-его спектрального канала
9. D1, D2, D3 – дисперсии величин спектральной яркости для 1..3 спектральных каналов

Величины параметров WLKn, KPLKn, совпадают с соответствующими по обозначению параметрами из АСУ-ЛР (см. таблицу 1) . Величины COUNTn формируются из данных АСУ-ЛР методом суммирования коэффициентов лесных формул (WLK) для соответствующих типов лесных пород.

Величины REFLECTn и Dn были получены путем автоматизированной совместной обработки космических снимков и цифровой повыдельной карты. Обработка проводилась следующим образом:

1. Лесоустроительная повыдельная карта и космические снимки были загружены в ArcView.
2. Космические снимки были конвертированы в GRID файлы.
3. При помощи ESRI Spatial Analyst была выполнена операция Summarize Zones для всех выделов. В качестве поля, определяющего зоны для которых рассчитывались статистические характеристики использовалось поле NV (номер квартала) таблицы лесоустроительных данных. В результате выполнения данной операции была получена таблица вида: NV, REFLECTn, D для каждого канала (для каждого снимка три таблицы).
4. Из таблиц были удалены лишние поля, затем три таблицы были сведены в одну при помощи средств ГИС ArcView.
5. Путем выполения операции JOIN полученная единая таблица была объединена с таблицей лесоустроительных данных. В результате и была получена таблица вида (2), которая является исходной для расчета статистических зависимостей.

Полученная таблица вида (2) была преобразована в формат статистического пакета SPSS для дальнейшей обработки.

Затем было произведено «отсеивание» ошибочных данных, которые могли появиться в результате неточного совмещения цифровой лесоустроительной карты и космических снимков, а также особенностей формирования повыдельной карты. Рассмотрим рисунок 2. На рисунке изображен фрагмент цифровой повыдельной карты и 3-ий канал космического снимка. Можно видеть яркую белую полосу, проходящую через изображение квартала №234022 – это лесная дорога шириной около 10 м. Такая ситуация может значительно исказить рассчитанные спектральные характеристики.

Для того, чтобы выяснить каким образом мы должны скорректировать данную ошибку построим гистограмму для вариации спектральной яркости каждого выдела в 3-ем спектральном канале (рисунок 3), в котором наиболее хорошо различимы лесопокрытые и нелесопокрытые территории. Путем изучения нескольких участков по которым проходят лесные дороги было установлено, что вариация на всех них превышает 15%. На рисунке 3 можно видеть небольшое увеличение количества участков с вариацией более 16%. Таким образом из исходных данных были удалены все участки вариация значений спектральной яркости в которых превышает 16%. Количество таких участков составило около 2% от общего их числа.

Заключительным этапом подготовки данных для вычислений лесной формулы явилась их группировка («агрегирование») по признаку «коэффициент лесной формулы по данным лесоустройства», т.е. по значениям COUNT20, COUNT10 и COUNT200. Необходимость группировки вызвана следующим. Рассмотрим рисунок 4, на котором приведена гистограмма распределения коэффициентов лесной формулы для темнохвойных пород. Видно, что это распределение крайне неравномерно.

Рисунок 4. Гистограмма для распределения коэффициента лесной формулы для темнохвойных пород.

Следовательно, данные нуждаются в «выравнивании» по их значимости для корректного построения регрессии. Выравнивание заключалось в агрегировании по коэффициентам лесной формулы, с вычислением средних арифметических значений для REFLECT1, REFLECT2 и REFLECT3. В дальнейшем, если не указано иначе, мы будем использовать обозначения REFLECTn для обозначения выровненных значений. Как показали исследования, условные распределения REFLECTn при фиксированных значениях COUNTn близки к нормальным (см рисунок 5), поэтому среднее арифметическое даст хорошую оценку.

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

Построение регрессии

Математико-статистическая модель имеет вид:

(3)

где e – остатки, случайная составляющая с независимыми значениями.

Предсказанные значения будем обозначать как PREDn.

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

PRED10 +PRED20 +PRED200 = 10 (4)

Расчет регрессии для темнохвойных пород.

На рисунке 6 показано поле значений для REFLECT1, REFLECT2 и REFLECT3 при зависимой переменной COUNT20. Вертикальными линями показан 95% доверительный интервал для математического ожидания.

Строим регрессию с независимой переменной REFLECT1 для участка №1. График для нормализованных отклонений приведен на рисунке 7.

Для проверки корректности построения регрессии проверяем результат на основании данных участка №2. При анализе результатов необходимо иметь ввиду, что снимки для участка №1 и участка №2 не были калиброваны (проблему калибровки предполагается решить в дальнейшем), по этой причине коэффициенты регрессии не будут совпадать. Однако можно получить определенные выводы относительно ошибки регрессии и значимости коэффициентов.

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

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

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

• Построены цифровые карты содержания темнохвойных пород, на основании найденных уравнений регрессии в виде GRID с шагом сетки 30 метров. В Приложении 3 приведена карта для участка №1.
• При помощи ГИС ArcView и расширения Spatial Analyst было вычислено среднее значение лесной формулы для лесных выделов, полученное по данным лесоустройства.
• Вычислены отклонения предсказанных значений от значений, полученных по лесоустройству ( E = PRED20 - COUNT20 ).

В результате проверки были получены данные о распределении ошибки (E) и ее основные статистические характеристики.

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

Проверка по исходным данным.

Проверка по данным Сюзьюнского лесничества. 

Расчет регрессии для светлохвойных пород.

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

Согласно экспериментальным данным основным спектральным каналом для определения % содержания светлохвойных должен быть 1-ый спектральный канал, однако, из рассмотрения графиков видно, что зависимость не выражена. В то же время в спектральных каналах 2 и 3 наблюдается слабая зависимость. Однако из экспериментальных исследований известно, что в данных спектральных каналах нет сколько-нибудь ярко выраженных особенностей  спектра, присущих древесным породам. Эта зависимость является следствием влияния какого-то другого фактора. Сравним графики для 2-ого и 3-его спектральных каналов приведенные рис 13 и рис 6. Можно видеть, что для темнохвойных пород «наклон» кривой противоположен по сравнению с аналогичной кривой для светлохвойных пород [6] . Возможно, это является следствием более низкой величины плотности биомассы крон (единица массы на единицу площади) для светлохвойных пород, по сравнению с темнохвойными. Так или иначе эта особенность может служить критерием для выделения светлохвойных пород, при условии что лесорастительные условия на всех «тестируемых» участках сходны.

Очевидно, что однозначно решить вопрос о выборе независимых переменных не удается, по этой причине было проведено три вида расчетов. В первый и второй расчет были включены все данные, не сгруппированные по значениям COUNT10. Эти вычисления использовались в основном для проверки значимости переменных, введенных в уравнения. Первый расчет был выполнен для значений COUNT10>4 (т.е. там, где светлохвойные можно считать основной породой). Во второй расчет были включены все данные. Необходимо заметить, что результаты второго расчета для Участка №2 показали, что R=0.03 и, следовательно результаты регрессии низко значимы, по этой причине они не были включены в итоговую таблицу. И, наконец, третий расчет был сделан по сгруппированным (агрегированным) данным для обоих участков. Расчет регрессии для Участка №2 оказался низко значимым (результаты не включены в итоговую таблицу). В третьем расчете была сделана попытка исключить переменную REFLECT3, потому что, результаты расчетов №1 и 2 показали ее низкую значимость.

Результаты расчета регрессии (сводная таблица).

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

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

Из сравнения рис 13 и 14 можно видеть, что наблюдается хорошее повторение зависимости REFLECT2(COUNT10) и REFLECT3(COUNT10). В регрессию была включена только переменная REFLECT3, т.к. значения в этом канале показывают наиболее гладкий и монотонный ход.

Результаты расчета регрессии для участка №1.

Результаты расчета регрессии для участка №2.

Расчет регрессии для мелколиственных пород

На рисунке 17 показано поле значений для REFLECT1, REFLECT2 и REFLECT3 при зависимой переменной COUNT200 (Участок №1). Вертикальными линями показан 95% доверительный интервал для математического ожидания.

Сильный скачок кривой и увеличение дисперсии для COUNT200 = 6 является следствием неблагоприятного сочетания ряда факторов в исходных данных. Он объясняется малым количеством выделов с COUNT200 = 6, их малыми размерами и неточным их совмещением с космическим снимком.

Результаты построения регрессии для COUNT200 c независимой переменной REFLECT1, участок №1.

Результаты построения регрессии для COUNT200 c независимой переменной REFLECT1, участок №2.

Проверка регрессии по формуле (4).

Как уже было сказано, для проверки корректности построения регрессионных уравнений можно использовать формулу (4). Кроме того, важно проверить результаты построения регрессии вычислив коэффициенты лесной формулы в среднем по выделу и сравнить их c данными лесоустройства. Выполнение этих видов проверок проводилось следующим образом. Для снимка №1 было выполнено вычисление коэффициентов лесных формул для темнохвойных, светлохвойных и мелколиственных пород, согласно найденным уравнениям регрессии. Для светлохвойных пород использовались результаты расчета №1 для 1-ого участка. Вычисленные коэффициенты для каждого пиксела космического снимка были представлены в GRID форме. Обозначим эти GRID карты как GRID10 для светлохвойных пород, GRID20 для темнохвойных и GRID200 для мелколиственных.

Выполнение проверки по формуле (4) заключается в построении производной GRID карты, показывающей отклонения d для каждого пиксела. Исключим из этого уравнения PRED10, т.к. регрессия для него оказалась малозначимой.

PRED20 + PRED200 –10 = d (6)

Гистограмма для d, для участка №1 приведена на рисунке 21. Можно видеть, что большинство ошибок лежит в пределах -1-0, а ошибок меньших -1 нет вообще. Соответствующий картографический материал, показывающий пространственное распределение ошибок приведен в Приложении 3. Из его рассмотрения видно, что большинство «ошибочных» пикселов сосредоточено на границах с нелесными землями.

Ниже представлены результаты аналогичного вычисления расхождений Е для мелколиственных пород (участок №1).

Заключение

В процессе выполнения работы были проведены вычисления, позволяющие построить зависимости между % содержанием темнохвойных, светлохвойных и мелколиственных древесных пород на лесных выделах и интегральными спектральными характеристиками этих выделов полученными по данным дистанционного зондирования (камера МСУ-Э). 

Для оценки коэффициентов лесной формулы использовался аппарат линейной регрессии. Средняя ошибка оценки составила:

• для темнохвойных пород 1.45 при высокой значимости регрессии,
• для светлохвойных 1.49, при очень низкой значимости,
• для мелколиственных 1.66, при высокой значимости. 

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

Анализ графиков стандартизованных отклонений показывает, что линейная регрессия в целом достаточно хорошо описывает исходные данные, но для более точного описания необходимо использовать кривые более высокого порядка (3-его).

В дальнейших исследованиях необходимо выполнить следующее:

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

Список литературы

1. Автоматизированная система управления лесными ресурсами. Описание комплекта задач 0273329.М4181.005 3В. КарНИИЛП, 1988 г.
2. Проект организации и ведения лесного хозяйства Ухтинского лесхоза. Т1 Объяснительная записка, г. Вологда, Северо-западное лесоустроительное предприятие. Лесоустройство 1996 г
3. В.А. Колемаев и др., Теория вероятностей и математическая статистика, Москва 1991 г.
4. Справочник по прикладной статистике (ред. Э.Ллойд и У.Линдерман), пер с англ. Москва, 1989 г.
5. А.Н. Поляков, Практикум по лесной таксации и лесоустройству, Москва, 1998 г.