Содержание к диссертации
Введение
1. Методы интерпретации данных дистанционного зондирования 9
1.1. Развитие методов дистанционного зондирования 9
1.2. Методы предварительной обработки спутниковых изображений 12
1.2.1. Радиометрическая коррекция 12
1.2.2. Расчет абсолютных значений яркости 13
1.2.3. Создание композитных изображений и масок 13
1.2.4. Геометрическая трансформация изображений 14
1.2.5. Радиометрическая коррекция влияния рельефа местности 16
1.3. Дополнительные признаки для оценки параметров растительности по данным дистанционного зондирования 16
1.3.1. Зональные отношения 17
1.3.2. Локальные пространственные характеристики 19
1.4. Интерпретация качественных и количественных характеристик 23
1.5. Методы классификации 23
1.5.1. Неуправляемая классификация 23
1.5.2. Классификация с обучением 24
1.5.3. Методы снижения размерности пространства признаков 29
1.5.4. Оценка качества классификации 33
1.6. Методы оценки количественных характеристик по данным дистанционного зондирования 34
1.6.1. Линейные регрессионные модели 34
1.6.2. Методы декомпозиции спектральных смесей 36
1.7. Обобщенный контекст сцены 39
1.8. Актуальные проблемы тематической интерпретации дистанционных данных 42
2. Методы сопряженного анализа данных спутниковой съемки и ГИС для оценки качественных и количественных характеристик лесных территорий 44
2.1. Принцип сопряженного анализа разнородных данных в средствах ГИС 44
2.2. Метод автоматизированного формирования эталонных выборок 46
2.2.1. Расчет статистических сигнатур выделов 46
2.2.2. Исключение пограничных зон выделов 48
2.2.3. Учет размера и формы выделов 49
2.2.4. Анализ степени пространственной однородности выделов 50
2.2.5. Стратификация выделов по степени однородности. 54
2.3. Метод непараметрической регрессии для оценки количественных характеристик растительности 55
2.4. Методика для классификации категорий земель и оценки запасов надземной фитомассы древесного яруса 64
2.5. Выводы к разделу 2 78
3. Состав и структура комплекса ГИС "Леса Средней Сибири" 80
3.1. Назначение, функции и структура банка данных многоуровневой ГИС "Леса Средней Сибири" 80
3.2. Структура лесоинвентаризационной ГИС локального масштабного уровня 82
3.2.1. Цифровая картографическая основа 82
3.2.2. База данных ГИС о лесном фонде 83
3.2.3. Материалы космической съемки 86
3.3. Аппаратурные, программные и функциональные средства комплекса ГИС 88
3.3.1. Аппаратурные средства для обработки геоданных 88
3.3.2. Программные средства 89
3.3.3. Функциональные возможности комплекса на базе PC Arclnfo и ArcView GIS 91
3.4. Разработанное программное обеспечение 93
4. Оценка породного состава и запаса фитомассы лесных насаждений по данным спутниковой съемки Landsat ТМ с использованием ГИС 96
4.1. Характеристика объекта исследования 96
4.2. Выбор данных спутниковой съемки 101
4.3. Предварительная обработка данных 103
4.4. Классификация категорий земель 115
4.5. Количественная оценка запаса надземной фитомассы по данным спутниковой съемки Landsat ТМ методами линейной и непараметрической регрессии 122
4.5.1. Линейные регрессионные модели оценки запаса 122
4.5.2. Оценка запаса древостоев по методу непараметрической регрессии 133
Основные научные результаты и выводы 138
Список использованной литературы 140
Приложение 151
- Локальные пространственные характеристики
- Методика для классификации категорий земель и оценки запасов надземной фитомассы древесного яруса
- Функциональные возможности комплекса на базе PC Arclnfo и ArcView GIS
- Линейные регрессионные модели оценки запаса
Локальные пространственные характеристики
В пространственных, геометрических и географических характеристиках изображений земной поверхности содержится очень большой объем информации. Опытный оператор-дешифровшик, вероятно, полагается на такую информацию (т. е. на форму, текстуру, ориентацию, контекст и другие характеристики) в большей степени, чем на спектральные признаки. Что касается машинных методов интерпретации и понимания изображений, то здесь пространственные аспекты еще не нашли должного учета. Можно, конечно, упомянуть множество «перспективных методов интерпретации» на основе пространственной информации, но есть все основания утверждать, что они не нашли широкого употребления в повседневной практике. Необходимость в развитии методов «пространственного мышления» для нужд дистанционного зондирования сейчас выдвигается на первое место [96].
Локальные пространственные признаки определяются относительно областей изображения небольшого размера, состоящих обычно из нескольких (возможно, до сотен) ближайших соседей данного элемента. Наиболее знакомый пример локального признака — это текстура объекта. Текстура представляет собой локальный пространственный признак изображения, который легко воспринимается зрительной системой человека. Хотя она давно используется как характерное свойство классов изображений земной поверхности, общего определения текстуры еще нет. Можно согласиться с тем, что текстура связана со степенью локальной изменчивости в изображении, однако разнообразие ее проявлений, видимо, очень велико.
Самым очевидным способом описания текстуры в некоторой области изображения представляется оценка соответствующего распределения пространственных частот [103] Расчет двумерного распределения пространственного энергетического спектра выполняется на основе преобразования Фурье [21] Другой близкий по идеологии, но более экономный с точки зрения вычислительных затрат прием заключается в определении локальной пространственной частоты пересечения изменяющегося уровня яркости с уровнем средней локальной яркости [80].
Более детальное описание текстуры можно получить путем расчета относительных частот, с которыми каждый уровень яркости встречается вблизи каждого другого уровня яркости. В результате расчета образуется набор так называемых «матриц пространственной смежности», которые, однако, непосредственно не используются [67]. По этим матрицам рассчитывается ряд статистических характеристик, служащих описанием текстуры. Они, например, несут в себе информацию о зернистости, контрастности и направленности текстуры. Экспериментально установлено, что такого рода признаки более эффективны при идентификации классов земной поверхности, чем пространственные спектры, однако связанные с ними вычислительные затраты тоже велики [103].
Многие современные исследования сосредоточены на понимании фундаментальной природы текстуры. Предложены различные модели, в которых изображение рассматривается как двумерный случайный процесс, состоящий из полезного сигнала и шума [15, 70], но практического применения такие модели пока не нашли.
Локальные пространственные признаки служат основой для разбиения (сегментации) сцены на однородные области или «объекты» разных типов. Среди таких признаков можно отметить локальную текстуру [55] и локальную однородность [73]. В работах [17, 97] разработана модель сцены, в которой классификация каждого элемента зависит от данных, содержащихся в четко определенной локальной окрестности (элемент, как правило, находится в центре). Показано, что точность классификации можно повысить, однако вычислительные затраты таковы, что практическая реализация этого подхода возможна лишь на вычислительной системе, обладающей высокоэффективной архитектурой [45, 93] .
При дистанционном изучении растительности неоднократно отмечалось перспективность использования текстурных признаков аэрокосмических изображений, которые имеют биологическую природу и отражают пространственную структуру полога древостоев или их сочетаний [33]. Например, в работе [3] предлагается использовать выборочные смешанные моменты совместного распределения вероятностей значений Vi=V(x1) и Vj=V(Xj) элементов изображения, находящихся на расстоянии от одного до трех шагов дискретного растра. Совместное распределение P(V(Xi), V(Xj)/r) оценивалось частотами nCkCi{r)lN(r) , где пс с (г) - число пар точек {xif Xj) с плотностями V(xi)=Ck и V(XJ) =С] ; {С1г С2, ., Сп) возможные значения плотностей V(x); п - число уровней квантования V(x); г - расстояние между отсчетами, г = 1, 2, 3, и Л/(г)-число пар элементов в окне фрагмента, имеющих взаимное расстояние г. Для расчета индексов изменчивости можно использовать дискретные (квантованные) значения яркостей в спектральных каналах спутникового изображения или составленное на их основе композитное изображение, значения пикселей которых дискретно меняется в интервале 0 - 255.
Число различных пикселей в окне (NDC), число пикселей отличных от центрального {CVN) и число различных пар пикселей (ВСМ) [85] также могут использоваться в качестве характеристик пространственной неоднородности участков дистанционного изображения. Как показали наши предыдущие исследования [29], вышеперечисленные пространственные признаки дают хорошие результаты при оценке продукционных характеристик древостоев по спутниковым изображениям со спутника «Ресурс-01» (МСУ-Э, разрешение 30x45м).
Методика для классификации категорий земель и оценки запасов надземной фитомассы древесного яруса
Для выполнения поставленных в работе задач классификации категорий земель и оценки запасов надземной фитомассы была выработана последовательность этапов обработки и анализа дистанционных данных (рис. 2.4) с использованием данных лесоинвентаризационной ГИС, реализованной на базе пакета ArcView GIS [59] с модулем растрово-векторного пространственного анализа Spatial Analyst, пакета IDRISI и дополнительных разработанных на языке Avenue процедур (см. раздел 3).
Предварительная обработка спутниковых изображений
Первичная радиометрическая коррекция
Первичная радиометрическая коррекция изображений включает в себя процедуру выравнивания яркости сцены в строках различных детекторов. Поканальная коррекция изображений производится по методу, рассмотренному в разделе 1.2.1 и реализованному в процедуре DESTRIPE пакета IDRISI.
Создание композитного изображения
Композитное изображение формируется для визуального анализа спутниковой сцены с целью дальнейшего выбора на нем точек привязки к топографической основе банка данных ГИС.
Для изображений спутника Lndsat ТМ композитное RGB-изображение рассчитывается, как линейная комбинация 5, 4 и 3 каналов (средний инфракрасный, ближний инфракрасный и красный спектральные диапазоны соответственно, см. таб. 4.4.). Формирование композитного осуществляется процедурой COMPOSIT пакета IDRISI.
Создания изображений-масок по методу разрядной декомпозиции яркостного сигнала При наличии в сцене спутникового изображения облаков, их теней и водных поверхностей необходимо выполнить исключение соответствующих им пикселей из дальнейшего анализа. Это достигается при помощи изображений-масок, в которых пикселям с нулевым значением соответствуют исключаемые области. Изображение маски формируется путем «отсечения» пикселей композитного изображения со слишком большими (облака, блики) и низкими (тени, водная поверхность) значениями яркости. Для этих целей был предложен метод разрядной декомпозиции яркостного сигнала [29]
Полученные слои содержат компоненты яркостного сигнала объектов изображения, причем каждый i-ый слой отражает величину изменения исходного сигнала из диапазона величиной в 22. Это позволяет разделять в двоичных слоях элементы с различной яркостью и частотой изменения сигнала.
Изображения-маски, соответствующие старшим разрядам значений яркости пикселя (обычно - разряды б или 7), в большинстве случаев отображают помехи от ярких объектов (облака блики, водной поверхности и т.п.), исключаемых из дальнейшего анализа.
Метод двоичной декомпозиции реализуется в процедуре BDC, разработанной на языке Avenue пакета ArcView GIS.
Совмещение дистанционных изображений с картами ГИС На этом этапе снимки трансформируются в картографическую проекцию для дальнейшего совмещения с векторными слоями банка данных ГИС о лесном фонде. В качестве базовой проекции для совместного анализа пространственных данных удобно использовать универсальную поперечную проекцию Меркатора (UTM) [76]. Параметры проекции выбираются в зависимости от номера территориальной зоны района исследований. Например, для стандартной картографической трапеции 0-4 6 с центральным меридианом 93 с.ш., номер зоны в проекции UTM равен 46.
На основании известных географических координат угловых точек спутниковых сцен вычисляются координаты всех остальных пикселей изображения, которые затем подвергаются преобразованию проекции с использованием процедуры PROJECT (IDRISI).
Для совмещения спутниковых изображений с топографической основой ГИС формируется сеть точек привязки (реперных точек), опознаваемых на композитном изображении и векторной топографической основе. Нанесение и редактирование точек привязки осуществляется при помощи средств визуального анализа ArcView GIS. Для геометрической трансформации растра используется процедура GRID.WARP (ArcView GIS, Spatial Analyst), выполняющая полиномиальное преобразование заданного порядка по методу, описанному в п. 1.2.4.
Выбор необходимого количества точек привязки зависит от степени используемого полиномиального преобразования и желаемой точности совмещения, которая контролируется при помощи таблицы неувязок точек привязки, и общей среднеквадратичной ошибки (RMSE). Ошибка RMSE в единицах измерения карты (обычно это метры) не должна превышать размер пикселя дистанционного изображения (30 м - для Landsat ТМ) В противном случае выполняется проверка точности выбора точек привязки, увеличивается их количество и повышается степень полиномиального преобразования.
Вторичная радиометрическая коррекция Обработка спутниковых изображений для местности с пересеченным рельефом должна включать радиометрическую коррекцию влияния топографических эффектов. Для этого использовался алгоритм коррекции по Ламберту, описанный в разделе 1.2.5. Обработка растрового изображения с целью коррекции осуществляется при помощи разработанной на Avenue процедуры RADNORM. Процедура использует растровые слои углов наклона и экспозиции местности. Входными параметрами процедуры являются зенитный угол Солнца и азимут направления на Солнце, определяемые из файла описания спутниковой сцены.
Расчет дополнительных признаков
Для оценки параметров растительности по данным дистанционного зондирования могут использоваться 4 группы признаков:
1) Значения яркостей пикселей в спектральных каналах изображения.
2) Значения зональных отношений (вегетационных индексов NDVI, NDVIc, см. п. 1.3.1)
3) Характеристики пространственной изменчивости значение индексов R, Er D, F, NDC, CVN, ВСМ (см. п. 1.3.2).
4) Дополнительные априорные контекстуальные признаки (значения высоты, экспозиции и угла наклона поверхности при наличии цифровой модели рельефа).
Расчет зональных отношений
На основании исходных изображений спектрозональных каналов: красного (R) , ближнего инфракрасного (IR) и среднего инфракрасного (MR) выполняется расчет вегетационных индексов с помощью процедуры IMIGE CALCULATOR (ArcView GIS, Spatial Analyst). В результате чего формируются два новых растровых изображения, значения пикселей которых соответствуют значениям NDVI и NDVIc. Рассчитанные вегетационные индексы затем используются в качестве дополнительных переменных в моделях оценки характеристик растительности по алгоритмам классификации, линейной и непараметрической регрессии.
Расчет индексов пространственной изменчивости
Индексы R, Н, D, F, NDC, CVN и ВСМ, рассмотренные в разделе 1.3.2, рассчитываются по ранее созданному композитному изображению или изображению одного из спектральных каналов. Для оценки изменчивости полога растительности может быть использован ближний инфракрасный канал (4-й канал для Landsat ТМ) Расчет индексов изменчивости осуществляется процедурой TEXTURE (IDRISI). Результатом расчетов на данном этапе являются изображения индексов изменчивости, также используемые в моделях оценки характеристик растительности.
Функциональные возможности комплекса на базе PC Arclnfo и ArcView GIS
Отображение географические данных на экране компьютера, включая возможность легко создавать карты по существующим источникам пространственных данных. Поддерживаются форматы шейп-файлов ARC/INFO, PC ARC/INFO, ArcCAD, AutoCAD (DXF, DWG), Intergraph (DGN). Импорт картографических данных из Maplnfo, Atlas GIS и ASCII.
Отображение табличных данных, включая возможность импортировать табличные данные и присоединять их к картографическим данным для пространственного анализа. Возможно прямое использование таблиц баз данных из ASCII, dBase, INFO, ACCESS, Oracle, FoxBase, SQL Server, Sybase, Paradox, DB2, Ingres и любых ODBC/SQL совместимых баз данных.
Создание и редактирование пространственных данных.
Определение атрибутов (описаний) любых объектов карты.
Отображение объектов темы различными цветами и символами согласно значениям их атрибутов.
Выбор объектов по значениям их атрибутов.
Визуализация табличных данных с использованием диаграммы, возможность отображения их на карте.
Возможность расчета статистики по любым атрибутам объектов.
Выбор объектов на основе анализа их пространственного расположения по отношению к другим объектам (пространственные запросы).
Создание компоновки карты и ее вывод на принтер или плоттер.
Создание компоновки карты и ее экспорт для использования в другой программе.
Возможность разработки программ (скриптов) на языке Avenue (ArcView GIS) или макросов на языке SML (PC Arclnfo) для реализации адаптированных приложений с использованием большого набора функций анализа пространственных геоданных в векторной и растровой (с использованием модуля Spatial Analyst) форме.
Дополнительные возможности растрово-векторного анализа на базе модуля Spatial Analyst для ArcView GIS
Модуль расширения Spatial Analyst предоставляет функциональные возможности пространственного анализа по растровым изображениям (грид-темам) и векторным картам ArcView GIS.
С помощью модуля расширения Spatial Analyst возможно выполнение следующих функций пространственного анализа: вычисление расстояний между пространственными объектами, определение пространственной близости, вычисление плотности объектов растровых и векторных карт, интерполяция поверхности, вычисление уклона, экспозиции, построение отмывки рельефа, построение изолиний, расчет статистики по пикселям и зонам растровых и векторных карт, расчет гистограмм и кросс-табуляций площадей, формирование запросов к растровой карте, выполнение логических и арифметических операций между элементами растровых карт (картографический калькулятор), расчет статистик соседства, переклассификация данных.
Для обработки растровых изображений спутниковой съемки также используется пакет IDRISI for Windows v. 2.
Возможности IDRISI включают: функции предварительной обработки растров (фильтрация, радиометрическая коррекция и калибровка); возможности расчета текстурных характеристик изображений; процедуры классификации изображений (методы параллелепипеда, минимального расстояния, максимального правдоподобия и др.); возможности расчета статистических характеристик, таблиц соответствия и Каппа-статистик, а также ряд других функций обработки и анализа растровых изображений.
Линейные регрессионные модели оценки запаса
На основании анализа работ, выполненных по данным спектрозональной аэрофотосъемки [10, 61], зависимость запаса от спектральной яркости в красном и инфракрасном диапазонах предполагалась обратно пропорциональной, причем, по мнению авторов [61], в линейных моделях оценки целесообразно использовать натуральный логарифм от значения запаса (In(Z)).
Выборки зависимой и независимых переменных для построения линейных регрессионных моделей формировались по эталонным выделам территории Предивинского лесничества (площадь 90189 га) Для классификации множеств эталонных выделов, было использовано два подхода, определяющих область применения линейных моделей. В первом случае выделы группировались по преобладающей породе (с коэффициентом состава более 7 единиц) и наиболее представленным возрастным стадиям (табл. 4.10). Расчет моделей выполнялся отдельно для каждой выделенной группы.
Во втором случае эталонные выборки выделов формировались с учетом однородности лесорастительных условий и возрастных стадий насаждений. В этом случае породный состав насаждений внутри группы уже не ограничивался одной преобладающей породой, а имел некоторый смешанный состав. Области однородных лесорастительных условий были выделены по комплексу характеристик, включающих сочетание типологических единиц растительности, почв и элементов рельефа, по методике Рыжковой В. А. [41] Для группирования выделов по типам лесорастительных условий использовалась карта восстановительных рядов растительности [41, 42], составленная Рыжковой В.А. на основе данных лесоустройства, цифровой модели рельефа и полевых исследований. Выбор восстановительных рядов, однородных по лесорастительным условиям, в качестве основы для формирования моделей оценки запасов осуществлялся как для повышения качества моделей (в силу однородности не учитываемых моделью природных факторов), так и с целью оценки запасов древесины на различных восстановительных стадиях. Характеристики восстановительных рядов растительности приведены в приложении 1.
Используемая карта восстановительных рядов была совмещена в ГИС с векторной картой лесотаксационных выделов для территории Предивинского лесничества Болылемуртинского лесхоза. Также как и в первом случае, для обеспечения статистической достоверности линейных моделей расчет проводили только для наиболее представленных на данной территории восстановительных рядов и возрастных стадий (табл. 4.11) .
Для предварительного анализа зависимости запаса древесины от дистанционных характеристик спутникового изображения использовался статистический аппарат простого линейного регрессионного анализа. Для оценки степени линейной связи между независимыми и зависимой переменными помимо коэффициента корреляции Пирсона (г) использовался квадрат коэффициента корреляции (R - коэффициент детерминации), характеризующий степень адекватности модели данным эталонной выборки.
Для повышения качества оценок запаса также были построены множественные линейные регрессионные модели, количество и состав переменных в которых определялся по методу максимального модифицированного R (см. п. 2.4.1). Верхняя граница уровня статистической независимости моделей (р) для простых линейных зависимостей принималась равной 5% (р 0.05), для множественных линейных моделей верхнее значение р уменьшалось в зависимости от состава и числа независимых переменных в модели [8 9].
Рассчитанные простые линейные модели оценки запаса для отдельных древесных пород в различном возрасте и возрастных стадий приведены в таблице 4.12, а для восстановительных рядов и возрастных стадий - в таблице 4.13. В таблицах приведены только те модели, для которых значение р 0.05. На рисунке 4.11 показаны графики статистически независимых моделей (р 0.05), для которых достигается максимальное значение R (выделены жирным шрифтом в таблицах 4.12 и 4.13).
Анализ полученных моделей показал, что все статистически независимые модели, рассчитанные для преобладающих пород, а также большинство моделей для восстановительных рядов подтверждают ожидаемую обратную связь спектральной яркости в каналах СІ, С2 и СЗ с запасом биомассы насаждений, что объясняется уменьшением отраженной солнечной радиации с ростом удельного по площади содержания хлорофилла, поглощающего излучение в соответствующих спектральных диапазонах. Связь запаса и отраженной радиации в ближнем инфракрасном диапазоне (С4) в большинстве случаев оказалась положительной, что, согласно ранее проделанным работам [32], объясняется влиянием внутренней структуры листьев на отражение солнечного излучения в этом диапазоне волн.
Во всех случаях для каналов среднего инфракрасного диапазона (С5 и С7) наблюдается отрицательная корреляция с запасом, вызванная поглощением волн указанного диапазона (см. табл. 4.3) содержащейся в листьях водой [32]. Вегетационные индексы {NDVI и NDVIc) во всех статистически независимых моделях, как и ожидалось, присутствуют с положительным знаком.
Статистически независимые модели с NDVIc в основном получены для выборок с преобладанием хвойных насаждений (пихты (табл. 4.12) и 13-го восстановительного ряда елово-пихтовых насаждений (табл. 4.13)). Однако в большинстве случаев статистически независимая корреляция с запасом наблюдается для средних инфракрасных каналов С5 и С7 (рис. 4.11 в, г) .
Сравнительный анализ таблиц 4.12 и 4.13 показал, что для восстановительных рядов корреляция (значения R ) дистанционных признаков и запаса в среднем достигает больших значений, чем для однородных по породе участков. Это можно объяснить тем, что для однородных по комплексу лесорастительных условий выделов влияние неучтенных факторов на дисперсию дистанционных признаков меньше, чем для выделов, однородных только по породному составу и возрасту насаждений.
Отдельно следует отметить тот факт, что в состав статистически независимых моделей, построенных для восстановительных рядов 3 и 2 (табл. 4.13) вошли модели, где независимой переменной является индекс пространственного разнообразия (Я) , рассчитанный для окна 3x3 пикселя. Причем, изменение запаса с ростом индекса Я (неоднородности) во всех приведенных моделях является обратно пропорциональной. Наиболее высокое значение R2 в моделях с участием Я достигается для восстановительного ряда 3 и возрастной стадии 60 - 79 лет (рис. 4.11д) Влияние степени однородности насаждений на удельный запас в восстановительных рядах прослеживается, скорее всего, из-за смешанного породного состава насаждений в рядах, а также пространственной мозаичностью насаждений со средней и низкой полнотой. Так, например, для восстановительного ряда 3 характерно участие в составе насаждений березы, осины и сосны с незначительным преобладанием березы в возрастной стадии 60-79 лет. Модель (рис. 4.11д) отражает снижение среднего запаса выдела с ростом индекса пространственной изменчивости (мозачности) Я, что объясняется снижением полноты и сомкнутости березовых насаждений первого яруса.
Простые линейные модели, хотя и позволяют проанализировать функциональную зависимость запаса насаждений от дистанционных признаков, все же не обеспечивают достаточной точности количественных оценок, о чем свидетельствуют невысокие значения R . С целью повышения точности оценок были рассмотрены множественные линейные регрессионные модели, рассчитанные по тем же выборкам для пород (таб. 4.14) и восстановительных рядов (таб. 4.15) .
