CC BY
71
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ
Таблица 2
Основные вегетационные индексы, применяемые на практике при обработке данных ДЗЗ
Вегетационный индекс Особенности
NDVI Нормализованный относительный индекс растительности NIR - RED NDVI = NIR + RED
RVI Относительный ВИ NIR RVI = RED
IPVI Инфракрасный ВИ NIR NDVI+1 IPVI= = NIR+RED 2
WDVI Взвешенный разностный ВИ WDVI = NIR - S-RED, S - наклон почвенной линии
SAVI Почвенный ВИ NIR - RED SAVI= (NIR + RED+L)(1+L) где L - корректирующий коэффициент
TSAVI Трансформированный почвенный ВИ s(NIR -sRED-a) TSA VI = (aNIR + RED - as+X (1+s2)) где a - координата пересечения почвенной линии с осью NIR, S - наклон почвенной линии, X - коэффициент коррекции, для уменьшения почвенного шума (часто принимают равным 0.08)
MSAVI Модифицированный почвенный ВИ NIR - RED MSAVI= (1+L) NIR + RED + L где, L = 1-2S-NDVIWDVI S - наклон почвенной линии.
MSAVI2 Модифицированный почвенный ВИ - 2 XS4VI2 NIR - RED 2 NIR+1-yj (2 NIR+1)2 -8( NIR - RED) (NIR+RED+L)(1+L) ’ L = 2
GEMI Индекс глобального мониторинга окружающей среды (RED-0 125) 2( NIR - RED) +1,5 NIR + 0,5RED GEMI=E(1-0,25E) ; E= 1-RED NIR + RED +0,5
ARVI ВИ, устойчивый к влиянию атмосферы NIR - Rb ARVI = NIR + Rb где, Rb = RED -A-(RED-BLUE). Как правило, А=1, но при малом покрытии растительности и неизвестном типе атмосферы А=0,5
SARVI Почвенный ВИ, устойчивый к влиянию атмосферы NIR - Rb SARVI= (NIR + Rb)(1+L)
130
ЛЕСНОЙ ВЕСТНИК 4/2012
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ
Окончание таблицы 2
Вегетационный индекс Особенности
GVI ВИ зелености GVI = -0.29MSS4 - 0.56MSS5 + 0.6MSS6 + 0.49MSS7 GVI = -0.2848TM1 - 0.2435TM2 - 0.5436TM3 + 0.7243TM4 + 0.0840TM5 - 0.1800TM7 где MSSn - в n канале сенсора MSS, аналогично для сенсора TM.
SVI Спектральный вегетационный индекс. Обобщает спектральный образ растительного полога, характеризуя состояние его зеленой биомассы.
SWVI Индекс вариации влажности почвы. Участвует в детектировании повреждений лесов пожарами по данным SPOT-Vegetation
EVI Расширенный вегетационный индекс. Усовершенствованный вариант NDVI, который рассчитывается по формуле: 2( NIR - RED) EVI = (NIR +C1 ■RED-C2 -BL UE+L)(1+L) где NIR, RED и BLUE - коэффициенты спектральной яркости в ближней инфракрасной (0,840-0,876 мкм), красной (0,620-0,670 мкм) и голубой (0,459-0,479 мкм) зонах; L - корректирующий коэффициент для учета характера подстилающей растительностью поверхности («1); С1 («6) и С2 («7,5) - коэффициенты, регулирующие степень использования голубой спектральной зоны в атмосферной коррекции красной спектральной зоны
DVI Разностный ВИ. DVI = NIR - RED
PVI Перпендикулярный ВИ PVI = sin(a)-NIR - cos(a)-RED a - угол между почвенной линией и осью NIR.
PRI Фотохимический индекс отражения. Индекс используется при оценке производительности растительности и выявлении стрессового состояния. Применяется для измерения изменений в каратоноидной пигментации (чувствителен к пигменту ксантофил) в живой листве. Каратоноидный пигмент может указывать на фотосинтетическую эффективность света или скорость поглощения углекислого газа листвой на единицу поглощенной энергии PRI = R531 — R570 R531 + R570 R - КСЯ соответствующих частот
TVI Трансформированный ВИ
MTVI Модифицированный трансформированный индекс MTVI(m) = L5[L2(Rto(S - RbanJ - 2.5(RM - RW2)V [(2Rband5 + 1)2 - (6RW5 - 5 Rba„J - 0.5] m -измерение мультиспектральным датчиком RapidEyeTM.
MCARI Модифицированный показатель поглощения хлорофилла MCARI(m) = [((Rba„M - Rba„J - °.2) (Rband fRband2)] ("Rband JRband 3) m -измерение мультиспектральным датчиком RapidEye™.
вегетационных индексов. Вегетационный индекс (ВИ) - это показатель, рассчитываемый в результате операций с разными спектральными диапазонами (каналами) и имеющий отношение к параметрам растительности в данном пикселе снимка. Эффективность ВИ определяется особенностями отражения электромагнитных волн. Индексы выведены, главным образом, эмпирически, применяются на
каждом конкретном участке с определенными особенностями. Существует большое число различных ВИ (Табл. 2), многие прикладные программы позволяют создавать удобные для пользователя собственные ВИ (ENVI, ERDAS IMG). Известно более 150 видов различных ВИ, однако из опыта мониторинга лесов центрального региона РФ число ключевых и самых распространенных ВИ сокращается до 5-10.
ЛЕСНОЙ ВЕСТНИК 4/2012
131
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ
Основное предположение по использованию ВИ состоит в том, что некоторые математические операции с разными каналами ДЗЗ могут дать полезную информацию о растительности. Это подтверждается множеством опытов. Второе предположение - это идея, что открытая почва на снимке будет формировать в спектральном пространстве прямую линию (почвенную линию). Почти все распространенные вегетационные индексы используют только соотношение красного и ближнего инфракрасного каналов, предполагая, что в ближней инфракрасной области лежит линия открытой почвы. Подразумевается, что эта линия означает нулевое количество растительности.
Большинство вегетационных индексов очень плохо работает для территорий с разреженным растительным покровом, спектр снимка в этом случае в основном зависит от почвы. Почвы могут различаться очень сильно по отражению, даже если для анализа используются очень широкие спектральные диапазоны.
Исходя из основного определения лесного выдела в лесоустроительных работах и учета методов ДЗЗ, предлагается создание специально адаптированных для космического мониторинга выделов на основе классификации полей значений ВИ лесного участка без участия эксперта в автоматическом режиме. Формирование дистанционно-ориентированных выделов (ДОВ) осуществляется алгоритмами кластеризации на заранее определенном участке леса, который классифицирован методами управляемой классификации с участием эксперта. Учитывая наибольшую популярность и проведя сравнительный анализ преимуществ и недостатков основных алгоритмов кластеризации (PAM, k-means, CURE, MST, Fuzzy C-means, HCM), выбрали алгоритм ISODATA (Iterative Self-Organizing Data Analysis Techniques). Его осуществление подразумевает первоначальное определение некоторого числа параметров при участии эксперта:
1. K - число возможных кластеров;
2. I - максимальное число итераций;
3. Qn - пороговый уровень для минимального числа пикселей в каждом кластере (необходимо для отбраковки кластеров);
4. Qs - пороговый уровень для среднеквадратического отклонения значений ВИ (необходимо для разделения пикселей);
5. Qc - пороговый уровень разницы (расстояния в пространстве значений) между двумя значениями пикселя поля ВИ (необходимо для слияния пикселей).
Алгоритм ISODATA в задаче формирования ДОВ выглядит следующим образом:
Шаг 1. Произвольно выбирается k (оно может и не равняться K), инициализируются кластерные центры: m m2,... тк из общего набора значений ВИ (х., i=1, 2,.,N}.
Шаг 2: Для каждого из N значений определяются наиболее корректные кластерные центры исходя из того, что
х е W , если Dl (х, m.)= min( DL (х, m.), i=1,.,k},
где W - кластер;
Dl - разность значений или расстояние в
пространстве значений).
То есть определяется близость значений к произвольно выбранному кластерному центру.
Шаг 3: Затем отбрасываются кластеры с количеством пикселей меньше Q т.е. если для какого-то j, R< Qk, то кластер W отбрасывается и к<—к-1, где R. - пиксели в окрестности кластера m., Qk - пороговое число пикселей k-го кластера.
Шаг 4: Выполняется преобразование кластерных центров согласно
mj = -1 Z х, (j=1,.,k);
Rj xW
Шаг 5: Вычисляется среднее расстояние Dj в пространстве значений кластера W с полученными кластерными центрами
Dj = Z Dl(Xmj), j=1,---,k)
Rj xW
Шаг 6: Рассчитывается общее среднее расстояние значений по всем кластерам 1 к
d = —ZR,dt.
K j=i
Шаг 7: Проверяется условие. Если k<K/2 (т.е. слишком мало кластеров), то выполняются операции расщепления полученных кластеров (Шаг 8), а если k>2K (слишком много кластеров), то операции слияния
132
ЛЕСНОЙ ВЕСТНИК 4/2012
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ
кластеров (Шаг 9). В иных случаях проверяется условие порога итераций (Шаг 14) кластеризации, если оно достигает I, то процесс прекращается, иначе все выполняется заново (Шаг 2).
Шаг 8: Операции расщепления кластеров проходят в несколько этапов.
Формируется вектор значений среднеквадратических отклонений о. = [ор ..., oj для каждого кластера
(j)
о
-1 X (х - m )2, j=1,^,k),
Rj
где о. - среднеквадратическое отклонение значений в кластере W .
Шаг 9: Затем выполняется поиск максимальной компоненты в векторе о. .
Шаг 10: Если для какого-либо о вы-
max
полняется условие
0max > Qs> Dj > D, R > 2Q^ то разделяется m• на два новых кластерных центра m + и
m~ , которые образуются прибавлением ±5 соответствующего значению отах кластера. Где 5 = а о , для произвольного а >0. Прежний кластерный центр удаляется и k<—k+1. Алгоритм переходит на начальный Шаг 2, иначе на Шаг 14.
Шаг 11: Это первый шаг по слиянию кластеров. Вычисляется попарно расстояние в пространстве значений (разница) D между двумя кластерными центрами
D = Dl (т ,, mi), для всех i ij. Сортируется k(k-1)/2 значений D . в порядке возрастания.
Шаг 12: Определяются не более P наименьших D . , которые меньше чем Qc и сортируются в порядке возрастания: D < Di2j2 <
~ ~ i
Шаг 13: Выполняется попарное объединение кластеров по кластерным центрам с условием о том, что если ни mh ни mh не используются в текущей итерации, то они сливаются в один кластерный центр
m =-----1--[R m. + R. m. ],
L ir ir JL JL
R + R,
lL JL
где L = 1, ..., P
После чего предыдущие варианты кластерных центров удаляются и k——k-1 и выполняется Шаг 2.
Шаг 14: Процесс завершается, если максимальное число итераций достигнуто I, в ином случае возвращаются к Шагу 2.
Алгоритм ISODATA, реализованный в ENVI, использует поле значений вегетационного индекса, для определения соответствующего кластера (класса) для каждого пикселя. Процесс начинается с назначения случайного (приближенного) среднего значения кластера (m) и повторяется до тех пор, пока это значение не достигнет величины среднего для каждого кластера исходных данных. Начальные средние значения кластеров распределяются равномерно вдоль центрального вектора пространства значений.
При первой итерации кластеризации поле значений ВИ1 равномерно разбивается на области, центром каждой из которых являются средние значения кластеров (mpm2, ... mn). Пиксели анализируются с левого верхнего угла изображения к нижнему правому. Вычисляется разница между значением пикселя ВИ1 и средним значением кластера. Пиксели назначаются в тот кластер, где разница минимальна (см. Шаг 2). Затем рассчитывают реальные средние значения ВИ1 в полученных кластерах, их средние значения меняются из-за попавших новых пикселей (m’pm’2, ... m’n). Выполняется вторая итерация, в процессе которой повторяют кластеризацию с новыми средними значениями и рассчитывают границы кластеров. После этого определяют новые средние значения и выполняют новую итерацию. В процессе второй итерации снова определяется минимальная разница между значениями ВИ1 пикселя и новыми средними значениями самих кластеров, по окончании пиксели перераспределяются. Перераспределения выполняются до тех пор, пока все пиксели поля значений ВИ1 с заданной вероятностью (порог сходимости) не попадут в свой кластер. В случае неспособности перераспределения пикселя в кластер алгоритм ограничивается заданным числом итераций I. Окончательно сформированные кластеры поля значений ВИ1 запоминаются для поиска пересечений с кластерами полей значений ВИ2, который формируется аналогично вышеописанной процедуре (рис.2). На этапе
ЛЕСНОИ ВЕСТНИК 4/2012
133
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ
ВИ1 ВИ2
о СО ДОВ^—0
Рис. 2. Этап формирования ДОВ. Пространственное пересечение кластеров ВИ1 и ВИ2
кластеризации отдельных ВИ указываются ограничения по минимальному числу пикселей в ДОВ (полученных классах). Следующий этап процесса формирования ДОВ подразумевает ограничение и по максимальному количеству пикселей, их образующих, исходя из пространственного разрешения снимка.
Формирование ДОВ может быть выполнено после поиска пространственных пересечений большего числа кластеров различных ВИ, однако в этом случае снижается возможность их пересечения для выделения общего участка. Основные этапы формирования ДОВ представлены на рис. 3. В случае отсутствия пространственных пересечений кластеров ДОВ формируется на основе приоритетного ВИ, который выбирается экспертом.
К наиболее простым в реализации вегетационным индексам можно отнести NDVI, EVI, ARVI ,IPVI, GEMI. Их расчет можно осу-
Рис. 3. Основные этапы формирования ДОВ
ществить в специализированном программном продукте ENVI. Ниже приведем пример формирования ДОВ на основе одного приоритетного индекса NDVI по спутниковому снимку Landsat ETM+ за июль 2010 г.
Согласно рис.3 первоначально определяется область исследования, т.е. территория леса. Шаг 1: Общий снимок подвергается управляемой или неуправляемой классификации с участием эксперта. По полученному классу лесов создаются векторные слои для ограничения от других классов. Выбирается конкретный интересующий лесной массив (рис. 4).
Шаг 2: Проводится вычисление необходимых вегетационных индексов. В нашем примере ограничимся расчетом наиболее популярного ВИ - NDVI, определяя его как приоритетный (рис. 5).
Шаг 3: Выполняется кластеризация полей значений ВИ алгоритмом ISODATA (рис. 6). При его реализации вводится ряд параметров, которые определяют итоговую картину сформированных ДОВ. На сегодня остается открытой задачей поиск оптимального варианта кластеризации и выработка критериев отбора наилучших ДОВ, поэтому в нашем случае этот процесс возложен на участие эксперта.
Шаг 4: Так как других ВИ не было, то сразу выполняется процесс векторизации классифицированных участков леса (рис.7). В итоге получаем повыделенное разбиение лесного массива, которое ориентировано именно на дистанционный спутниковый мониторинг с учетом однородности значений ВИ.
Модель формирования ДОВ встраивается в общий функционал специализированного программного продукта ENVI. Ее
134
ЛЕСНОЙ ВЕСТНИК 4/2012
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ
Рис. 4. Выделенный лесной массив
Рис. 5. Поле значений NDVI исследуемого участка леса
Рис. 6. Кластеризация полей ВИ
Рис. 7. Векторный слой ДОВ
реализация осуществляется в диалоговом исполнении с минимальным участием пользователя. Основные этапы выполняются в автоматическом режиме. Так как автоматический режим не исключает появления ошибок и браков, предложенный в настоящий момент вариант модели формирования ДОВ подразумевает тесное участие эксперта. Дальнейшие исследования будут направлены на отладку конкретных этапов алгоритма с минимизацией субъективного воздействия пользователя.
Концепция использования ДОВ может быть полезна при непрерывном многомесячном или многолетнем мониторинге лесов с выявлением количественных изменений значений вегетационных индексов, которые, в свою очередь, имеют известные регрессионные связи с биологическими параметрами леса. Интересен многолетний контроль изменения непосредственно самих границ ДОВ, сформированных с помощью спутниковых снимков определенной даты годового сезона.
Таким образом, метод мониторинга на основе ДОВ может дать объективную обобщенную информацию о состоянии участка леса и его изменении во времени, упрощая при этом процедуру поиска важных качественных и количественных признаков.
Библиографический список
1. Официальный сайт ФГУП «Рослесинфорг». [Электронный ресурс]: Сайт ФГУП «Рослесинфорг». Режим доступа http://www.roslesinforg.ru/activity/ forest.
2. Журнал Промышленник России №10 (121) 2010.
3. Сухих, В.И. Аэрокосмические методы в лесном хозяйстве и ландшафтном строительстве: учебник / В.И. Сухих. - Йошкар-Ола: МарГТУ, 2005. -392 с.
4. Барталев, С.А. Разработка методов оценки состояния и динамики лесов на основе данных спутниковых наблюдений: автореф. дисс. ... д-р техн. наук / С.А. Барталев. - М.: ИКИ РАН, 2007. - 48 стр.
5. Данилин, И.М. Высокие технологии XXI века для аэрокосмического мониторинга и таксации лесов / И.М. Данилин, Е.М. Медведев, Н.И. Абэ и др.
ЛЕСНОЙ ВЕСТНИК 4/2012
135
