УТОЧНЕНИЕ  ЗНАЧЕНИЙ  НОРМАЛИЗОВАННОГО  ВЕГЕТАТИВНОГО  ИНДЕКСА  (NDVI),  МЕТОДОМ  НАЛОЖЕНИЯ  ТРАНСПИРАЦИОННОЙ  МАСКИ

Жолобов  Денис  Алексеевич

канд.  техн.  наук,  декан  факультета  «Математики  и  информационных  технологий»  доцент  Астраханского  Государственного  Университета,  РФ,  г.  Астрахань

E-mail: 

Баев  Александр  Викторович

магистрант  Астраханского  Государственного  Университета,  РФ,  г.  Астрахань

E-mail: 

 

RECTIFICATION  OF  THE  NORMALIZED  VEGETATION  INDEX  (NDVI),  A  METHOD  OF  IMPOSING  TRANSPIRATION  MASK

Denis  Jolobov

candidate  of  Science,  Dean  of  the  Faculty  of  Mathematics  and  Information  Technologies  assistant  professor  of  Astrakhan  State  University,  Russia,  Astrakhan

Aleksandr  Baev

undergraduate  of  Astrakhan  State  University,  Russia,  Astrakhan

 

АННОТАЦИЯ

Целью  данной  работы  является  теоретическое  обоснование  применения  транспирационной  маски  для  уточнения  вегетативного  индекса  на  территориях  пустынь  и  полупустынь.  Был  использован  метод  нормализации  и  анализа  геоданных,  рассчитаны  сопутствующие  индексы.  Результатом  работы  стало  создание  эмпирического  метода  наложения  транспирационной  маски  на  вегетативный  индекс  (NDVI).

ABSTRACT

The  purpose  of  this  work  —  theoretical  substantiation  of  application  transpiration  masks  for  rectification  of  the  vegetation  index  in  desert  and  near  desert  territory.  We  used  the  method  of  normalization  and  analysis  of  geodata,  calculated  associated  indexes.  The  result  was  the  creation  of  the  empirical  method  overlay  transpiration  mask  to  vegetation  index  (NDVI).

 

Ключевые  слова:  нормализация  геоданных;  вегетативный  индекс.

Keywords:  normalization  of  geodata;  vegetative  index.

 

Введение

 

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

Описание  района  исследования

Основное  исследование  проводится  в  районе  Богдинско-Баскунчаского  заповедника,  включающего  в  себя  наивысшую  точку  Северного  Прикаспия  —  гору  Большое  Богдо  (абсолютная  отметка  +149,6  м;  протяженность  —  около  5  км)  и  уникальное  солёное  озеро  Баскунчак  (абсолютная  отметка  —21,3  м).  Остальная  территория  заповедника  представляет  собой  слегка  всхолмленную  равнину  (со  средней  абсолютной  отметкой  +15  -  +20  м)  с  общим  уклоном  поверхности  к  центру,  котловине  озера.  В  рельефе  местности  выделяются  поднятия  Куба-Тау  (+37  м)  и  Вак-Тау  (+22,4  м),  понижение  —  урочище  Кривая  лощина  (-15—9  м),  расположенные  восточнее  озера  Баскунчак.  Кроме  того,  рельеф  местности  сильно  осложнен  эрозионными  и  карстовыми  формами:  размеры  отдельных  котловин  достигают  100  м  в  диаметре  и  25  м  глубины.  Эрозионные  формы  широко  распространены  на  западном  берегу  озера  Баскунчак  и  представлены  многочисленными  оврагами  и  балками  длиной,  как  правило,  не  более  2  км.  Также  в  районе  имеется  несколько  понижений  —  лиманов.

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

Среди  современных  источников  геоданных  проще  всего  использовать  источник  свободных  геоданных  USGS(Геологической  Службы  Соединённых  Штатов)  на  сайте:  http://www.usgs.gov/  и  http://glovis.usgs.gov/  .  Данные  с  этих  ресурсов  представлены  в  формате  GeoTIFF  в  виде  непрерывных  наборов  сцен  для  различных  районов  мира,  в  том  числе  и  для  исследуемой  местности.  Представленные  на  ресурсах  геоданные  соответствуют  уровню  обработки  LG1:  «сырые  геоданные»  без  радиологической  и  атмосферной  нормализации.  Для  самостоятельного  проведения  нормализации  и  атмосферной  коррекции  —  доступны  файлы  настройки  спутниковых  сенсоров  с  набором  параметров  калибровки  спутниковых  датчиков:  минимумы  и  максимумы  пиксельной  яркости  на  изображении;  минимумы  и  максимумы  излучения  на  датчиках;  минимумы  и  максимумы  отражения  излучения  от  поверхности  земли  (Landsat  8)  и  многое  другое.  Кроме  данных  настройки  спутниковых  датчиков  в  архиве  присутствуют  несколько  файлов  в  формате  GeoTIFF,  распределённых  по  номерам  каналов  —  количество  и  состав  которых  для  разных  спутников  (Landsat  5,7,8)  различный  (см  таб.  1).

Как  видно  из  таблицы,  все  шесть  используемых  в  работе  «спектров»  находятся  у  разных  спутников  в  сходном  диапазоне  частот.  Поэтому,  для  удобства  использования  терминологии,  вместо  указания  спутниковых  каналов  в  дальнейшем  будем  использовать  самоназвание  спектров:  BLUE  —  синий,  GREEN  —  зелёный,  RED  —  красный,  NIR  —  ближний  инфракрасный,  SWIR1  —  средний  инфракрасный  1,  SWIR2  —  средний  инфракрасный  2  —  вне  зависимости  от  спутника  и  применять  к  ним  одинаковый  набор  вычислений.

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

Таблица  1.

Спутники  и  используемые  в  работе  каналы

 

Радиологическая  коррекция

Радиологическая  коррекция  является  первым  этапом  нормализации  сырых  геоданных  и  представляет  из  себя  математическую  операцию  перевода  значения  яркости  пикселей  геоснимка  в  значения  радиации  поступившей  на  датчики  спутника.  Для  такого  перевода  в  комплекте  данных  Landsat  присутствует  файл  коррекции  *_MTL.txt,  предельные  значения  из  которого  и  используются  на  этом  этапе  обработки  геоснимка.

Для  проведения  работы  использовалась  следующая  формула  перевода  яркости  в  значение  top  of  atmosphere  radiance  (TOA  radiance)  Thome  et  al.,  1994  [18],  Lu  et  al  2002  [22]:

 

L λ  =  DNcal  ×  Gainλ  +  Baisλ

 

где:  Lλ=  Спектральная  радиация,  пришедшая  на  сенсор  спутника 

DNcal =  Значения  яркости  пикселя  сырого  геоснимка

Gainλ =  Радиометрическое  усиление

Baisλ =  Радиометрическое  смещение

 

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

Атмосферная  коррекция

Следующим  этапом  нормализации  геоданных  является  уменьшение  влияния  атмосферы  на  снимок  и  перевод  значений  радиации,  дошедшей  до  сенсоров  спутника  (TOA  radiance),  в  значения  реально  отражённого  от  земли  спектрального  излучения  солнечного  света. 

Влияние  атмосферы  на  геоснимок  проявляется  в  целом  ряде  факторов:  угол  падения  и  отражения  солнечных  лучей,  прозрачность  атмосферы,  газовый  фактор  и  дымка  (Рис.  1).

 

Рисунок  1.  Факторы,  влияющие  на  попадание  отраженной  солнечной  радиации  на  сенсоры  спутника

 

Для  дальнейших  исследований  необходимо  провести  оптическую  коррекцию  (нормализацию)  данных  геоснимка  методом  Dark  Object  Subtraction  (DOS),  впервые  представленным  Chavez  (1996)  [23].  Суть  метода  состоит  в  нахождении  яркости  однопроцентно  тёмного  объекта  геоснимка  с  последующей  коррекцией  минимума  значений  каждого  пикселя  изображения,  относительно  спектральной  яркости  найденного  объекта.

Есть  два  основных  способа  поиска  1  %  темного  объекта  (Dark  Object)  для  метода  DOS:

1.  Эмпирический  метод  —  подразумевает  поиск  значений  в  ручном  режиме,  например,  с  использованием  инструмента  гистограмма  в  QGIS,  где  изменяя  нижний  порог  яркости  гистограммы  постепенно  находим  примерное  значение  яркости  искомого  тёмного  объекта.

2.  Вычислительный  метод  подразумевает,  что  суммарная  яркость  (от  0  до  n)  однопроцентно  тёмного  объекта  будет  соответствовать  0,01  %  от  суммарной  яркости  всех  пикселей  геоснимка  Sobrino  et  al.,  2004  [24]. 

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

После  определения  яркости  Dark  Object  (в  дальнейших  вычислениях  будем  обозначать  его  как  DNmin),  производим  атмосферную  коррекцию  по  методу  DOS  в  несколько  этапов:

1)  Вычисляем  значение  радиации,  соответствующее  значению  яркости  1  %  темного  объекта  (расчёт  производится  по  аналогии  с  TOA  radiance)

 

L λmin  =  DNmin  ×  Gainλ  +  Baisλ

 

где:  Lλmin=  Спектральная  радиация  для  1  %  тёмного  объекта 

DNmin =  Значения  яркости  пикселя  1  %  тёмного  объекта

Gainλ =  Радиометрическое  усиление

Baisλ =  Радиометрическое  смещение

 

 

2)  Рассчитываем  коэффициент  влияния  угла  падения  и  отражения  солнечных  лучей  для  1  %  тёмного  объекта

 

0.01  *  cos² θ  *  T_Z  *  E₀

L_1%  =  ---------------------------

π  *  d²

 

где:  L_1%  =  Коэффициент  влияния  угла  падения  и  отражения  солнечных  лучей  для  1%  тёмного  объекта

d  =  Расстояние  от  солнца  до  земли  в  астрономических  единицах  в  конкретный  день  съёмки  сцены  на  конкретной  местности  [29]

E₀   =  Коэффициент  солнечного  внеатмосферного  спектрального  излучения  (явно  представлен  как  табличные  данные  и  учитывается  при  калибровке  датчиков  Landsat  5  и  7,  для  Landsat  8  дополнительно  вычисляется)

θ   =  Зенитное  расстояние  для  солнца  в  радианах

T_Z   =  Мера  прохождения  излучения  от  солнца  до  земли,  в  методе  DOS2,  принимается  равным  cosθ

3)  Вычисляем  значение  атмосферной  дымки  (hazing)

 

Lλhaze  =  Lλmin  -  L_1%

 

где:  Lλhaze  =  значение  атмосферной  дымки  (hazing) 

L_1%   =  коэффициент  влияния  угла  падения  и  отражения  солнечных  лучей  для  1%  тёмного  объекта

Lλmin =  Спектральная  радиация  для  1  %  тёмного  объекта 

4)  Рассчитываем  атмосферно  скорректированные  значения  отражённой  солнечной  радиации  по  формуле:

 

π  *  (Lλ  —  Lλhaze)  *  d²

ρλ  =  ----------------------------------------

E₀  *  cosθ  *  T_Z

 

где:  ρλ  =  Атмоферно  скорректированные  значения  отражённой  солнечной  радиации

Lλ  =  Значения  радиации  пришедшей  на  сенсор  спутника

Lλhaze   =  Значение  атмосферной  дымки  (hazing)

d   =  Расстояние  от  солнца  до  земли  в  астрономических  единицах  в  конкретный  день  съёмки  сцены  на  конкретной  местности  [29]

E₀   =  Коэффициент  солнечного  внеатмосферного  спектрального  излучения  (  явно  представлен  как  табличные  данные  и  учитывается  при  калибровке  датчиков  Landsat  5  и  7  ,  для  Landsat  8  дополнительно  вычисляется  )

θ   =  Зенитное  расстояние  для  солнца  в  радианах

T_Z   =  Мера  прохождения  излучения  от  солнца  до  земли,  в  методе  DOS2,  принимается  равным  cosθ

 

Особенности  работы  с  растрами  спутника  Landsat  8

Часть  использованных  при  исследовании  сцен  были  сняты  спутником  Landsat  8  датчик  Oli.  Возникла  необходимость  их  сравнения  со  сценами  спутников  Lansat  5  (датчик  TM)  и  Landsat  8  (датчик  ETM+).  В  данном  контексте  известна  проблема  применения  стандартных  методов  атмосферной  коррекции  для  нового  спутника:  дело  в  том,  что  калибровка  датчиков  Oli  Landsat  8  производится  без  учёта  значений  коэффициента  солнечного  внеатмосферного  спектрального  излучения  Eo  (Sobrino  et  al.,  2004)  [24]  или,  в  других  источниках,  Esun.  Вместо  использования  данного  коэффициента,  в  файл  корректировки  *_MTL.txt,  были  добавлены  некоторые  новые  спектральные  параметры:  REFLECTANCE_MULT_BAND  —  усиление  значения  отражения  и  REFLECTANCE_ADD_BAND  —  смещение  значения  отражения  для  каждого  из  спектральных  датчиков.  В  результате,  по  задумке  авторов  изменений,  расчет  TOA  reflectance  для  Landsat  8  должен  приобретать  следующий  вид:

 

ρλ'   =  MρQcal+  Aρ

 

где:  ρλ'  =  Значение  верхнеатмосфеного  планетарного  отражения  радиации  (TOA  reflectance),  без  учёта  коррекции  по  углу  падения  и  отражения  солнечных  лучей.

Mρ  =  каналоспецифичный  мультипликативный  расчётный  фактор  для  (REFLECTANCE_MULT_BAND_x,  где  x  это  номер  канала)  —  усиление  значения  отражения

Aρ  =  каналоспецифичный  мультипликативный  расчётный  фактор  для  (REFLECTANCE_ADD_BAND_x,где  x  это  номер  канала)  —  смещение  значения  отражения 

Qcal   =  Значения  яркости  пикселя  сырого  геоснимка  (DN)

 

А  коррекция  TOA  reflectance  с  учётом  угла  падения  и  отражения  солнечных  лучей  выглядит  так:

 

 

где:  ρλ  =  Значение  верхнеатмосфеного  планетарного  отражения  радиации  (TOA  reflectance),  c  учётом  коррекции  по  углу  падения  и  отражения  солнечных  лучей.

θSE   =  Солнечный  элеватор.  Доступен  в  файле  *_MTL.txt  в  параметре  (SUN_ELEVATION).

θSZ   =  Солнечный  зенит;  θSZ=  90°  -  θSE

 

Атмосферная  коррекция  по  методу  DOS,  на  данный  момент,  не  реализована  для  нового  метода  обработки  геоснимков  Landsat  8.  Коме  того  сообщество  разработчиков  свободной  GIS  GRASS  указывает  на  одинаковые  значения  REFLECTANCE_MULT_BAND  и  REFLECTANCE_ADD_BAND  для  всех  каналов  снимка  [25],  что  не  может  быть  в  реальности.  Данная  группа  разработчиков  в  своём  модуле  нормализации  и  атмосферной  коррекции  i.landsat.toar  применяет  к  георастрам,  снятым  при  помощи  датчиков  Oli  Landsat  8,  те  же  математические  методы,  что  для  TM  Landsat  5  и  ETM+  Landsat  7,  а  недостающий  коэффициент  Eo  (Esun)  рассчитывает  по  формуле:

 

(PI  *  d^2)  *  RADIANCE_MAXIMUM

Esun  =  ---------------------------------------------------

REFLECTANCE_MAXIMUM

 

где:  Esun  =  (E₀)  вычисленный  коэффициент  солнечного  внеатмосферного  спектрального  излучения

d   =  расстояние  от  солнца  до  земли  в  астрономических  единицах  в  конкретный  день  съёмки  сцены  на  конкретной  местности  (координаты)  [27]

RADIANCE_MAXIMUM  =  каналоспецифичный  мультипликативный  расчётный  фактор  для  (RADIANCE_MAXIMUM_x,  где  x  это  номер  канала)  —  максимально  возможное  значение  поступающей  на  сенсор  радиации

REFLECTANCE_MAXIMUM  =  каналоспецифичный  мультипликативный  расчётный  фактор  для  (REFLECTANCE_MAXIMUM_x,  где  x  это  номер  канала)  —  максимальное  значение  радиации,  отражённой  от  поверхности  земли

 

Ещё  одним  нововведением  от  Landsat  8  стало  уменьшение  чувствительности  каналов  RED,  NIR  и  SWIR1  относительно  Landsat  7  и  5,  что  приводит  к  изменению  значений,  вычисляемых  из  этих  спектров,  индексов. 

Данную  проблему  пытался  решить  Neil  Flood,  2014  [26]  введением  дополнительных,  рассчитанных  им  эмпирическим  путём,  коэффициентов  для  каждого  канала  геоснимка.  В  результате,  коррекция  значений  TOA  reflectance  Landsat  8  в  значения  для  того  же  региона  и  тех  же  спектров  Landsat  7,  приобретает  следующий  вид:

 

  ρETM+  =  c0  +  c1  *  ρOLI

 

где:  ρETM  =  результат  конвертирования  значения  TOA  reflectance  из  Landsat  8  в  Landsat  7

ρOLI  =  вычисленное  для  Landsat8  значение  TOA  reflectance

c0  =  смещение  значения  отражённой  радиации  между  датчиками  OLI  и  ETM+

c1  =  усиление  значения  отражённой  радиации  между  датчиками  OLI  и  ETM+

 

Значения  c0  и  с1  были  сведены  автором  в  корректировочные  таблицы,  что  позволяет  применять  их  без  повторения  вычислений  представленных  Neil  Flood,  2014  [26]  в  своей  статье  «Continuity  of  Reflectance  Data  between  Landsat-7  ETM+  and  Landsat-8  OLI,  for  Both  Top-of-Atmosphere  and  Surface  Reflectance:  A  Study  in  the  Australian  Landscape»

Итог  нормализации  всего  множества  геоснимков  для  трёх  типов  спутников  Landsat:  5,7,8 

Несмотря  на  сложность  проведения  единообразной  обработки  трёх  типов  геоснимков  —  удалось  получить  однородные  нормализованные  сцены  исследуемой  области.  Все  вычисления  производились  с  использованием  языка  программирования  python.  Для  этого  была  написана  собственная  программа  нормализации  геоданных  с  использованием  интерфейса  GDAL(библиотеки  для  работы  с  растровыми  географическими  форматами),  тесно  связанного  с  расширением  numpy,  добавляющим  в  язык  программирования  поддержку  больших  многомерных  массивов  и  матриц,  а  также  систему  низкоуровневых  математических  функций  для  операций  с  этими  массивами.  Программа  состоит  из  следующих  структурных  элементов:

1.  Класса  преобразования  растра  в  массив.

2.  Класса  преобразования  вычисленного  массива  в  растр.

3.  Класса  сбора  данных  коррекции  при  помощи  парсера  *_MTL.txt  файла  и  данных  выгрузки  расстояний  от  земли  до  солнца  в  исследуемый  период  [29]

4.  Выполняемого  скрипта  самих  математических  вычислений  растров  для  проведения  процесса  нормализации.

Исходные  коды  программы  доступны  для  загрузки  по  адресу  [1],  кроме  того  исходные  GEOTIFF  и  нормализованные  спектральные  растры  доступны  по  адресу  [4]. 

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

Вегетативные  индексы

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

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

Расчет  большей  части  вегетационных  индексов  базируется  на  двух,  не  зависящих  от  прочих  факторов,  участках  кривой  спектральной  отражательной  способности  растений.  На  красную  зону  спектра  RED  (620—750  нм)  приходится  максимум  поглощения  солнечной  радиации  хлорофиллом,  а  на  ближнюю  инфракрасную  зону  NIR  (750—1300  нм)  —  максимальное  отражение  энергии  клеточной  структурой  листа.  То  есть  высокая  фотосинтетическая  активность,  связанная,  как  правило,  с  большой  прикорневой  фитомассой  растительности,  ведет  к  более  низким  значениям  коэффициентов  отражения  в  красной  зоне  спектра  и  большим  значениям  в  ближней  инфракрасной.  В  результате  математическая  разница  или  частное  от  деления  этих  показателей  позволяет  четко  отделять  растительность  от  прочих  природных  объектов  Jordan  (1969)  [85]. 

Наиболее  классическим  вариантом  вегетативного  индекса  на  сегодняшний  момент  является  нормализованный  разностный  ВИ  (Normalized  Difference  VI,  NDVI)  Rouse  et  al.  (1973)  [9],  рассчитываемый  по  формуле: 

 

NDVI   =  (NIR  —  RED)/(NIR+RED)

 

Индекс  может  принимать  значения  от  -1  ,  до  +1,  причём  любые  значения  ниже  нуля  не  имеют  логического  объяснения  как  проявление  активности  фотосинтезирующих  растений,  и  фактически  могут  игнорироваться  в  исследованиях  растительности  —  так  как  не  несут  смысловой  нагрузки.  Кроме  NDVI  существует  ряд  простых  ВИ:  инфракрасный  ВИ  (IPVI)  Crippen  (1990)  [11],  разностный  ВИ  (DVI)  Richardson  and  Everitt  (1992)  [9],  перпендикулярный  ВИ  (PVI)  Richardson  and  Wiegand  (1977)[14]  и  т.  д.,  рассмотрение  которых  не  планируется  в  рамках  данной  статьи.  Ограничимся  только  указанием  на  их  существование.

Теория  расчёта  вегетативных  индексов  исходит  из  того  что  в  местностях  с  низкой  плотностью  растительности  (менее  50  %),  к  котором  относится  и  большая  часть  площади  Богдинско  —  Баскунчаского  заповедника,  в  процессе  вычислении  ВИ  проявляется  влияние  «почвенного  шума»  —  состояния,  когда  значение  ВИ  искажается  посторонним  фактором  (почвы),  не  связанным  с  процессом  фотосинтеза.  Поэтому  было  введено  понятие:  «почвенная  линия»,  которая  имеет  один  наклон  в  пространстве  графика  RED-NIR  (рис.  2).

 

Рисунок  2.  График  почвенной  линии  относительно  точки  ВИ

 

Однако  часто  почвы,  на  одной  и  той  же  местности,  имеют  серьёзные  различия  —  что  можно  сказать  и  об  исследуемом  районе,  ввиду  наличия  большого  количества  типов  рельефа.  В  результате  почвенные  линии  имеют  разные  углы  наклона  на  одном  и  том  же  снимке.  Для  преодоления  влияния  «почвенного  шума»  была  разработана  группа  относительных  индексов  менее  чувствительных  к  плотности  растительного  покрова,  чем  NDVI.  Наиболее  активно  используемым,  из  почвенных  ВИ,  можно  считать  модифицированный  почвенный  ВИ-2  (Modified  Soil  Adjusted  VI,  MSAVI)  Qi  et  al.  (1994)  [18],  рассчитываемый  по  формуле:

 

SAVI   =  ((NIR-RED)/(NIR+RED-L))*(1+L)

 

Где:

 

L=1-(2  *  NIR  +  1  —  SQR(2  *  NIR  +  1)^2  —  8  *  (NIR  *  RED))/2

 

Индекс  обладает  тем  же  диапазоном  параметров,  что  и  классический  NDVI  (от  -1  до  +1)  и  той  же  логикой  использования:  отличие  заключается  в  уточнении  влияния  почвенной  линии  L  на  значение  NDVI. 

Помимо  «почвенного  шума»,  на  вычисление  ВИ  может  иметь  большое  влияние  состояние  атмосферы:  иногда  сильно  изменяемое  на  протяжении  одной  сцены,  особенно  на  территории  с  высотным  рельефом.  Атмосферное  влияние  изменяет  количество  света,  попадающее  на  сенсоры  спутника,  и  может  вызвать  ошибки  в  вычислении  индексов.  Несмотря  на  то,  что  подобные  атмосферные  феномены  не  типичны  для  исследуемого  района,  рассмотрим  индексы,  пытающиеся  решить  озвученную  проблему  без  применения  специальной  атмосферной  коррекции.  Эти  индексы  достигают  уменьшения  чувствительности  к  влиянию  атмосферы  ценой  уменьшения  динамического  диапазона.  В  целом,  они  менее  чувствительны  к  изменению  растительного  покрова  чем  NDVI,  но  если  растительность  невысока  (что  как  раз  и  наблюдается  в  районе  исследования),  они  подвержены  сильному  влиянию  «почвенного  шума». 

Среди  атмосферных  индексов  наиболее  известен  индекс  глобального  мониторинга  окружающей  среды  (Global  Environmental  Monitoring  Index,  GEMI),  разработанный  Pinty  and  Verstraete(1991)  [19]  и  вычисляемый  по  формуле:

 

GEMI  =  E*(1-0.25*E)-(RED-0.125)/(1-RED)

 

Где:

 

E  =  (2*(NIR^2-RED^2)+1.5*NIR+0.5*RED)/(NIR+RED+0.5)

 

Индекс  принимает  значения  от  0  до  +1.

Индекс  был  выведен  опытным  путём,  и  нечувствительность  к  атмосферному  влиянию  тоже  была  определена  чисто  эмпирически.  Существует  критика  использования  этого  метода:  Qi  et  al.  (1994)  [18]  показал  случаи  явных  ошибок  GEMI  из-за  почвенного  шума  при  низком  растительном  покрове,  что  критично  для  нашего  исследования.

Как  видно  из  уже  приведённых  примеров  ВИ,  любые  корректировки  расчёта  классического  NDVI  (почвенные  или  атмосферные)  приводят  к  увеличению  чувствительности  к  влиянию  помех  в  другой  мере  вычислений:  корректировки  почвенной  линии  L  —  увеличивают  атмосферные  помехи,  а  атмосферная  корректировка  —  влияние  «почвенного  шума»  на  территориях  с  низкой  плотностью  растительности.  Нахождение  сбалансированного  решения  относится  к  применимости  того  или  иного  вегетативного  индекса  к  конкретной  территории.

Водный  индекс

Помимо  ВИ,  ещё  один  тип  индексов  имеет  непосредственную  связь  с  жизнедеятельностью  растений  на  поверхности  земли,  водные  индексы  Water  Index  (WI).  Данная  связь  определяется  процессом  транспирации  наземных  растений.  В  исследуемом  районе,  транспирация  оказывает  гораздо  большее  постоянное  влияние  на  содержание  водяного  пара  в  приземном  слое  атмосферного  воздуха,  чем  другие  факторы  окружающей  среды.

Расчет  водного  индекса  базируется  на  двух  не  зависящих  от  прочих  факторов  участках  кривой  спектральной  поглощающей  способности  водяных  паров  в  приземном  слое  атмосферного  воздуха.  На  среднюю  инфракрасную  зону  SWIR1  (1500—2000  нм)  приходится  максимум  поглощения  солнечной  радиации  водяными  парами,  находящимися  в  нижних  слоях  атмосферы  у  самой  поверхности  земли,  а  на  ближнюю  инфракрасную  зону  NIR  (750—1300  нм)  —  максимальное  отражение  энергии  всей  поверхностью  земли  за  исключением  водоёмов.  В  результате  вычитания,  или  взвешенного  вычитания,  этих  двух  спектров  получаем  картину  приземной  концентрации  водяных  паров.  Расчёт  (Normalized  Difference  Water  Index,  NDWI)  Gao  (1996)  [27]  можно  выразить  формулой:

 

NDVI  =  (SWIR1  —  NIR)/(SWIR1  +  NIR)

 

Индекс  может  принимать  значения  от  -1  до  +1.  Причём  любые  значения  меньше  нуля  или  равные  нулю  могут  использоваться  как  маска  идентификации  водоёмов  и  других  объектов  с  открытой  водной  поверхностью.  Это  вызвано  тем,  что  отрицательные  и  нулевые  значения  NDWI  могут  появляться  только  при  значениях  NIR  =>  SWIR1,  что  характерно  для  поверхности  водоёмов,  активно  поглощающих  весь  спектр  инфракрасного  излучения.  Соответственно,  все  положительные  значения  индекса  можно  отнести  к  показателю  испарения  воды  с  поверхности  земли,  немалая  доля  которой,  в  засушливых  районах,  приходится  на  транспирацию  наземных  растений.

Анализ  соотношения  значений  вегетативного  и  водного  индексов. 

Рассмотрим  исследуемый  регион,  в  самый  сложный  для  вычисления  вегетативного  индекса,  период:  июнь-август.  В  это  время,  и  без  того  скудная  растительность  полупустыни  подвергается  еще  большей  деградации  и,  как  следствие,  «почвенный  шум»  значительно  изменяет  (смазывает)  значения  NDVI.  В  качестве  тестового  примера  выберем  сцену  LT51700262011152KHC01  со  спутника  Landsat  5  за  01.06.11  и  рассчитаем  для  неё  NDVI  (Рис.  3)

 

Рисунок  3.  Вегетативный  индекс  NDVI  сцены  LT51700262011152KHC01  для  исследуемой  области

 

Оценивая  получившуюся  карту  ВИ,  выделим  зоны  с  явно  неверными  значениями  вегетативного  индекса  (рис.  3):

Зона  №  1а.  Выделена  область  с  недавно  прошедшем  степным  пожаром,  её  яркость  сопоставима  с  рядом  расположенной  областью  пожаром  не  затронутой  (зона  1б)

Зона  №  2.  Здесь  слились  по  значениям  яркости  крупные  степные  балки,  с  искусственными  посадками  деревьев,  и  и  выше  расположенная  равнина.  Это  невозможно  в  середине  и  конце  вегетативного  периода.

Зона  №  3.  Наблюдается  однотипность  картины  ВИ  окрестностей  горы  Б.  Богдо,  что  не  может  быть  типично  для  исследуемой  территории.  Из-за  сложного  рельефа  местности  растительный  покров  имеет  ярко  выраженную  мозаичную  структуру.

Зона  №  4а.  Карстовые  формы  рельефа,  лишённые  древесной  растительности,  на  данном  изображении  сопоставимы  по  яркости  крупным  степным  балкам  с  искусственными  посадками  деревьев  (зона  4б).

Попытки  пересчитать  значения  вегетативного  индекса  с  использованием  методов  почвенных  ВИ  (SAVI  и  SARVI)  не  произвели  существенного  исправления  картины.  В  результате  возникла  идея  уточнить  NDVI  при  помощи  значения  транспирации  растений  как  корректирующего  фактора.  Для  этого  рассчитаем  еще  один  индекс,  непосредственно  связанный  с  интенсивностью  транспирации  —  водный  индекс. 

Для  составления  карты  водного  индекса  будут  использоваться  только  значения  NDWI>0,  так  как  нас  интересует  испарение  воды  c  поверхности  земли  (см.  глава  Водный  индекс),  приравняв  NDWI<0  к  0,  получим  следующую  формулу:

 

NDWIground  =  (NDWI  <  0)?  0,  NDWI

 

где:  NDWIground  —  значение  водного  индекса  в  околоземном  слое  почвы,  здесь  примерно  равный,  значению  транспирации  растений 

 

Для  краткости  NDWIground  в  дальнейшем  обозначим  как  NDWI.

Если  рассмотреть  изображения  обоих  рассчитанных  индексов  в  сравнении  (Рис.  4),  становится  хорошо  заметно,  что  в  большей  части  изображения  значения  водного  индекса  (NDWI)  инвертируют  значения  вегетативного  индекса  (NDVI).  Этому  явлению  есть  простое  объяснение: 

На  площадях  с  наибольшим  уровнем  ВИ  (а  значит,  теоретически,  значительной  прикорневой  фитомассой  растений)  —  будет  более  высокий  уровень  водяных  паров  за  счёт  транспирации,  что,  в  свою  очередь,  приводит  к  падению  водного  индекса  при  поглощении  отражённого  ИК  излучения  водяными  парами.

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

 

Рисунок  4.  Сравнение  изображение  карты  вегетативного  индекса  (слева)  и  водного  индекса  (справа)

 

Получившиеся  картина  инверсии,  независимо  вычисленных  индексов,  приводит  к  идее  учёта  влияния  значения  транспирации  (через  водный  индекс)  на  величину  ВИ.  Для  этого  попробуем  использовать  метод  вычитания  из  «позитива»  NDVI  «негатива»  NDWI  и  проанализируем  полученную  карту  разности  индексов  (рис.  5)

 

Рисунок  5.  Результат  вычитания  водного  индекса  из  вегетативного  индекса  (слева),  выделение  отрицательных  областей  разности  NDVI  и  NDWI  (справа)

 

Как  видно  на  рисунке  8  (слева),  разность  вычисленных  индексов  напоминает  NDVI  —  но  в  сильно  контрастном  исполнении.  Смазанность  фона  ВИ  теперь  отсутствует,  но  зато  в  резком  контрасте  проявляются  области  с  отрицательными  значениями  (рис.  5  справа).  С  точки  зрения  логики  NDVI,  получившаяся  карта  разности  не  может  отражать  проявление  фотосинтетической  активности  растений. 

Отрицательные  значение  разности  индексов  являются  проявлением  соотношения  NDWI  >  NDVI,  предположительно,  в  местностях,  с  наиболее  деградировавшей  в  летний  период  вегетации  растительностью.  Связано  это  с  тем,  что  при  почти  полном  отсутствии  фотосинтезирующей  биомассы,  значения  ВИ  будут  низкими  (даже  с  учётом  «почвенного  шума»),  а  водный  индекс  наоборот  высоким  —  из-за  низкого  содержания  водяного  пара.  Экспериментально  установлено,  что  для  исследуемой  области  подобная  картина  типична  для  середины  вегетативного  периода  (июнь,  июль,  август)  и  не  типична  для  начала  и  конца  периода  вегетации.  Соответственно,  использование  NDWI,  в  чистом  виде,  в  качестве  транспирационной  маски  вегетативного  индекса,  не  предоставляется  возможным. 

Для  превращения  водного  индекса  в  корректирующую  транспирационную  маску  для  NDVI,  можно  попытаться  применить  эмпирические  методы  с  последующей  проверкой  результатов  в  полевых  условиях.  Используем  следующую  логику:  если  значения  водного  индекса  велики  —  попробуем  их  «приглушить».  Наиболее  оптимально  в  качестве  подобного  «глушителя»  показала  себя  сумма  NDVI  и  NDWI,  а  произведение  данной  суммы  на  значение  водного  индекса  даёт  приемлемую  картину  транспирационной  маски  в  нужном  диапазоне  значений  (рис.  6): 

 

Tmask   =  NDWI  *  (NDVI  +  NDWI)

 

где:  Tmask  =  вычисленная  транспирационная  маска

 

Рисунок  6.  Сравнение  суммы  NDVI  и  NDWI  (слева)  с  вычисленной  транспирационной  маской  (справа)

 

Используем  разность  вегетативного  индекса  и  транспирационной  маски  для  вычисления  скорректированного  «контрастного»  транспирационного  ВИ:

 

TNDVI  =  NDVI  —  Tmask

 

где:  TNDVI  =  «контрастный»  транспирационный  ВИ

Tmask   =  вычисленная  транспирационная  маска

 

Проанализируем  получившуюся  карту  транспирационного  ВИ  (рис.  7)  относительно  замечаний  к  карте  NDVI,  описанных  в  начале  этой  главы:

Зона  №  1а.  Выделенная  область  с  недавно  прошедшим  степным  пожаром  теперь  имеет  намного  меньшую  яркость  относительно  рядом  расположенной  области,  не  повредившейся  выгоранию  (зона  1б)

Зона  №  2.  Яркость  ВИ  крупных  степных  балок,  с  искусственными  посадками  деревьев  теперь  намного  выше,  чем  яркости  рядом  расположенной  равнины.

Зона  №  3.  Однотипности  картины  ВИ  в  окрестностях  горы  Б.  Богдо  теперь  не  наблюдается,  хорошо  заметна  мозаичная  структура  растительности.

Зона  №  4а.  Карстовые  формы  рельефа,  лишённые  древесной  растительности,  теперь  имеют  наполовину  меньшую  яркость  относительно  крупных  степных  балок  с  искусственными  посадками  деревьев  (зона  4б).

 

Рисунок  7.  Транспирационный  вегетативный  индекс  ТNDVI  сцены  LT51700262011152KHC01  для  исследуемой  области

 

Заключение

Сравнение  изображения  всей  карты  исследуемой  области  в  режиме  NDVI  и  TNDVI  (рис.  8)  обращает  на  себя  внимание  более  низкая  яркость  и  высокая  контрастность  транспирационного  ВИ  относительно  NDVI.  Кроме  того,  отмеченные  в  главе  «Анализ  соотношения  значений  вегетативного  и  водного  индексов.»  неточности  вычисления  NDVI,  вызванные  сильным  «почвенным  шумом»,  были  успешно  нивелированы  применением  к  ВИ  транспирационной  маски.  Потому,  чисто  теоретически,  карта  TNDVI  больше  соответствует  исследуемой  области.  В  любом  случае,  требуются  полевые  исследования  для  подтверждения  или  опровержения  выдвинутой  теории  уточнения  ВИ  посредством  применения  транспирационной  маски. 

 

Рисунок  8  Сравнение  карты  NDVI  (слева)  и  транспирационного  ВИ  TNDVI  (справа)

 

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

1.Исходные  тексты  для  инструментария  норализации,  атмосферной  коррекции  и  калькуляции  растров  Landsat  5,7,8  [Электронный  ресурс]  —  Режим  доступа.  —  URL:  https://github.com/oldbay/raster_tools  (дата  обращения:  01.05.2015). 

2.Кренке  А.Н.,  Пузаченко  Ю.Г.  «Построение  карты  ландшафтного  покрова  на  основе  дистанционной  информации.»  Экологическое  планирование  и  управление  2008  №  2. 

3.Майорова  В.И.,  Банников  А.М.,  Гришко  Д.А.,  Жаренов  И.С.,  Леонов  В.В.,  Топорков  А.Г.,  Харлан  А.А.  «Контроль  состояния  сельскохозяйственных  полей  на  основе  прогнозирования  динамики  индекса  NDVI  по  данным  космической  мультиспектральной  и  гиперспектральной  съёмки.»  Наука  и  образование  2013  №  7.

4.Набор  нормализованных  спектральных  растров  и  вычисленных  индексов  для  текущей  статьи  [Электронный  ресурс]  —  Режим  доступа.  —  URL:  https://github.com/oldbay/paper_examples  (дата  обращения:  01.05.2015).

5.«Неформальное  сообщество  специалистов  в  области  ГИС  и  ДЗЗ»  [Электронный  ресурс]  —  Режим  доступа.  —  URL:  http://gis-lab.info/  (дата  обращения:  01.05.2015).

6.ALI  sensors,  Gyanesh  Chander,  Brian  L.  Markham  b,  Dennis  L.  Helder  «Summary  of  current  radiometric  calibration  coefficients  for  Landsat  MSS,  TM,  ETM+,and  EO-1»  Remote  Sensing  of  Environment  Volume  113,  Issue  5,  15  May  2009,  Pages  893—903

7.Meghan  Graham  MacLean,  Alexis  M.  Rudko,  Dr.  Russell  G.  Congalton  «Multi-temporal  image  analysis  of  the  coastal  watershed,  nh»  ASPRS  2010  Annual  Conference  San  Diego,  California  April  26-30,  2010

8.Jordan  C.F.  (1969)  "Derivation  of  leaf  area  index  from  quality  of  light  on  the  forest  floor,"  Ecology,  vol.  50,  —  pp.  663—666. 

9.Rouse  J.W.,  Haas  R.H.,  Schell  J.A.,  and  Deering  D.W.  (1973)  "Monitoring  vegetation  systems  in  the  great  plains  with  ERTS,"  _Third  ERTS  Symposium,  NASA  SP-351,  vol.  1,  —  pp.  309—317.

10.Kriegler  F.J.,  Malila  W.A.,  Nalepka  R.F.,  and  Richardson  W.  (1969)  "Preprocessing  transformations  and  their  effects  on  multispectral  recognition,  in  _Proceedings  of  the  Sixth  International  Symposium  on  Remote  Sensing  of  Environment,  University  of  Michigan,  Ann  Arbor,  MI,  —  pp.  97—131. 

11.Crippen  R.E.  (1990)  "Calculating  the  Vegetation  Index  Faster,"  Remote  Sensing  of  Environment,  vol  34.,  —  pp.  71—73.

12.Richardson  A.J.  and  Everitt,  J.  H.  (1992)  "Using  spectra  vegetation  indices  to  estimate  rangeland  productivity,  Geocarto  International,  vol.  1,  —  pp.  63—69. 

13.Lillesand  T.M.  and  Kiefer  R.W.  (1987)  Remote  Sensing  and  Image  Interpretation,  2nd  edition,  John  Wiley  and  Sons,  New  York,  Chichester,  Brisbane,  Toronto,  Singapore,  721  p. 

14.Richardson  A.J.  and  Wiegand  C.L.  (1977)  "Distinguishing  vegetation  from  soil  background  information,"  Photogrammetric  Engineering  and  Remote  Sensing,  vol.  43,  —  pp.  1541—1552.

15.Clevers  J.G.P.W.  (1988)  "The  derivation  of  a  simplified  reflectance  model  for  the  estimation  of  leaf  area  index,  Remote  Sensing  of  Environment,  vol  35.,  —  pp.  53—70.

16.Huete  A.R.  (1988)  "A  Soil-Adjusted  Vegetation  Index  (SAVI),"  Remote  Sensing  of  Environment,  vol.  25,  —  pp.  295—309. 

17.Baret  F.,  Guyot  G.,  and  Major  D.  (1989)  "TSAVI:  A  vegetation  index  which  minimizes  soil  brightness  effects  on  LAI  or  APAR  estimation,"  in  12th  Canadian  Symposium  on  Remote  Sensing  and  IGARSS  1990,  Vancouver,  Canada,  July`10—14.

18.Qi  J.,  Chehbouni  A.,  Huete  A.R.,  and  Kerr  Y.H.  (1994)  "Modified  Soil  Adjusted  Vegetation  Index  (MSAVI),"  Remote  Sensing  of  Environment,  vol.  48,  —  pp.  119—126.

19.Pinty  B.  and  Verstraete  M.M.  (1991)  "GEMI:  A  Non-Linear  Index  to  Monitor  Global  Vegetation  from  Satellites,"  Vegetation,  vol.  101,  15—20.

20.Kaufman  Y.J.,  Tanre  D.  (1992)  "Atmospherically  resistant  vegetation  index  (ARVI)  for  EOS-MODIS,  in  Proc.  IEEE  Int.  Geosci.  and  Remote  Sensing  Symp.  '92_,  IEEE,  New  York,  261—270.

21.Thome  K.J.,  Biggar  S.F.,  Gellman  D.L.,  Slater  P.N.  (1994)  —  Absolute-radiometric  calibration  of  Landsat-5  Thematic  Mapper  and  the  proposed  calibration  of  the  Advanced  Spaceborne  Thermal  Emission  and  Reflection  Radiometer.  Paper  presented  at  the  Geoscience  and  Remote  Sensing  Symposium,  1994.  IGARSS’94.  Surface  and  Atmospheric  Remote  Sensing:  Technologies,  Data  Analysis  and  Interpretation. 

22.Lu  D.,  Mausel  P.,  Brondizio  E.,  Moran  E.  (2002)  —  Assessment  of  atmospheric  correction  methods  for  Landsat  TM  data  applicable  to  Amazon  basin  LBA  research.  International  Journal  of  Remote  Sensing,  23  (13):  2651—2671. 

23.Chavez  P.S.  (1996).  Image-based  atmospheric  correction—revisited  and  improved.  Photogrammetric  Engineering  and  Remote  Sensing,  62(9),  1025—1036.

24.Sobrino  J.A.,  Jiménez-Muñoz  J.C.,  Paolini  L.  (2004)  —  Land  surface  temperature  retrieval  from  LANDSAT  TM  5.  Remote  sensing  of  Environment,  90  (4):  434—440.

25.«Calculates  top-of-atmosphere  radiance  or  reflectance  and  temperature  for  Landsat  MSS/TM/ETM+/OLI»  [Electronic  resource]  —  URL:http://grass.osgeo.org/grass65/manuals/i.landsat.toar.html  (accessed:  01.05.2015).

26.Neil  Flood  (2014)  Continuity  of  Reflectance  Data  between  Landsat-7  ETM+  and  Landsat-8  OLI,  for  Both  Top-of-Atmosphere  and  Surface  Reflectance:  A  Study  in  the  Australian  Landscape.  Remote  Sens.  2014,  6,  7952—7970

27.Gao  B.C.,  1996,  (NDWI—a  normalized  difference  water  index  for  remote  sensing  of  vegetation  liquid  water  from  space.  Remote  Sensing  of  Environment,  58,  —  pp.  257—26.

28.HANQIU  XU  (2006)  Modification  of  normalised  difference  water  index  (NDWI)  to  enhance  open  water  features  in  remotely  sensed  imagery  International  Journal  of  Remote  Sensing  Vol.  27,  №  14,  20  July  2006,  3025—3033

29.       «HORIZONS  Web-Interface»  [Electronic  resource]  —  URL:  http://ssd.jpl.nasa.gov/horizons.cgi  (accessed:  01.05.2015).