CC BY
337
УТОЧНЕНИЕ ЗНАЧЕНИЙ НОРМАЛИЗОВАННОГО ВЕГЕТАТИВНОГО ИНДЕКСА (NDVI), МЕТОДОМ НАЛОЖЕНИЯ ТРАНСПИРАЦИОННОЙ
МАСКИ
Жолобов Денис Алексеевич
канд. техн. наук, декан факультета «Математики и информационных
технологий» доцент Астраханского Государственного Университета, РФ,
г. Астрахань E-mail: denis.jolobov@aspu. ru Баев Александр Викторович магистрант Астраханского Государственного Университета, РФ,
г. Астрахань E-mail: [email protected]
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.
Спутники и используемые в работе каналы
Порядковый
Спутник и датчик номер канала Название спектра
(band)
Landsat 5 TM 1 Blue — Синий
Длинны волн (нм)
450—520
i) DociFreezer
created by free version of
Спутник и датчик
Landsat 7 ETM+
Lansat 8 Oli
Порядковый номер канала (band)
2
3
4
5 7 1 2
3
4
5 7 2
3
4
5
6 7
Название спектра
Green — Зелёный Red — Красный NIR — Ближний ИК SWIR1 — Средний ИК 1 SWIR2 — Средний ИК 2 Blue — Синий Green — Зелёный Red — Красный NIR — Ближний ИК SWIR1 — Средний ИК 1 SWIR2 — Средний ИК 2 Blue — Синий Green — Зелёный Red — Красный NIR — Ближний ИК SWIR1 — Средний ИК 1 SWIR2 — Средний ИК 2
Длинны волн (нм)
520—600 630—690 760—900 1550—1750 2080—2350 450—520 520—600 630—690 770—900 1550—1750 2080—2350 450—515 525—600 630—680 845—885 1560—1660 2100—2300
Радиологическая коррекция
Радиологическая коррекция является первым этапом нормализации сырых геоданных и представляет из себя математическую операцию перевода значения яркости пикселей геоснимка в значения радиации поступившей на датчики спутника. Для такого перевода в комплекте данных Landsat присутствует файл коррекции *_MTL.txt, предельные значения из которого и используются на этом этапе обработки геоснимка.
Для проведения работы использовалась следующая формула перевода яркости в значение top of atmosphere radiance (TOA radiance) Thome et al., 1994 [18], Lu et al 2002 [22]:
LX = DNcal x GainX + BaisX
где: LX= Спектральная радиация, пришедшая на сенсор спутника
ВЫса1= Значения яркости пикселя сырого геоснимка Gainl= Радиометрическое усиление BaisX= Радиометрическое смещение
После обсчёта каждого пикселя геоснимка с использованием данной формулы, получаем матрицу значений 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)
LXmin = DNmin x GainX + BaisX
где: LXmin= Спектральная радиация для 1 % тёмного объекта DNmin= Значения яркости пикселя 1 % тёмного объекта GainX= Радиометрическое усиление BaisX= Радиометрическое смещение
2) Рассчитываем коэффициент влияния угла падения и отражения солнечных лучей для 1 % тёмного объекта
0.01 * cos2в * Tz * Eo
Li% =---------------------------
п * d2
где: L1% = Коэффициент влияния угла падения и отражения солнечных лучей для 1% тёмного объекта
d = Расстояние от солнца до земли в астрономических единицах в конкретный день съёмки сцены на конкретной местности [29]
E0 = Коэффициент солнечного внеатмосферного спектрального излучения (явно представлен как табличные данные и учитывается при калибровке датчиков Landsat 5 и 7, для Landsat 8 дополнительно вычисляется) в = Зенитное расстояние для солнца в радианах
TZ = Мера прохождения излучения от солнца до земли, в методе DOS2, принимается равным cosG
3) Вычисляем значение атмосферной дымки (hazing)
LXhaze = LXmin - L1%
где: LXhaze = значение атмосферной дымки (hazing)
L1% = коэффициент влияния угла падения и отражения солнечных лучей для 1% тёмного объекта
LXmin= Спектральная радиация для 1 % тёмного объекта
4) Рассчитываем атмосферно скорректированные значения отражённой солнечной радиации по формуле:
п * (LX — LXhaze) * d2
pX =----------------------------------------
Eo * cose * Tz
где: pX = Атмоферно скорректированные значения отражённой солнечной радиации
LX = Значения радиации пришедшей на сенсор спутника LXhaze = Значение атмосферной дымки (hazing)
d = Расстояние от солнца до земли в астрономических единицах в конкретный день съёмки сцены на конкретной местности [29]
E0 = Коэффициент солнечного внеатмосферного спектрального излучения ( явно представлен как табличные данные и учитывается при калибровке датчиков Landsat 5 и 7 , для Landsat 8 дополнительно вычисляется ) в = Зенитное расстояние для солнца в радианах
TZ = Мера прохождения излучения от солнца до земли, в методе DOS2, принимается равным cosG
Особенности работы с растрами спутника 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 должен приобретать следующий вид:
pXf = MpQcal+ Ap
где: pXf = Значение верхнеатмосфеного планетарного отражения радиации (TOA reflectance), без учёта коррекции по углу падения и отражения солнечных лучей.
Mp = каналоспецифичный мультипликативный расчётный фактор для (REFLECTANCE_MULT_BAND_x, где x это номер канала) — усиление значения отражения
Ap = каналоспецифичный мультипликативный расчётный фактор для (REFLECTANCE_ADD_BAND_x,где x это номер канала) — смещение значения отражения
Qcal = Значения яркости пикселя сырого геоснимка (DN)
А коррекция TOA reflectance с учётом угла падения и отражения солнечных лучей выглядит так:
где: рХ = Значение верхнеатмосфеного планетарного отражения радиации (TOA reflectance), c учётом коррекции по углу падения и отражения солнечных лучей.
OSE = Солнечный элеватор. Доступен в файле *_MTL.txt в параметре (SUN_ELEVATION).
OSZ = Солнечный зенит; OSZ= 90° - OSE
Атмосферная коррекция по методу 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) рассчитывает по формуле:
рХ1
PÁ =
cos (ASE) sin(eSE)
(PI * dA2) * RADIANCE_MAXIMUM
Ю DociFreezer
created by free version of
Esun =---------------------------------------------------
REFLECTANCE_MAXIMUM
где: Esun = (Eo) вычисленный коэффициент солнечного внеатмосферного спектрального излучения
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, приобретает следующий вид:
pETM+ = c0 + c1 * pOLI
где: pETM = результат конвертирования значения TOA reflectance из Landsat 8 в Landsat 7
pOLI = вычисленное для Landsat8 значение TOA reflectance c0 = смещение значения отражённой радиации между датчиками OLI и ETM+
c1 = усиление значения отражённой радиации между датчиками OLI и ETM+
Значения c0 и с1 были сведены автором в корректировочные таблицы, что позволяет применять их без повторения вычислений представленных Neil Flood, 2014 [26] в своей статье «Continuity of Reflectance Dat a 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. Наиболее активно используемым, из почвенных ВИ, можно считать модифицированный
сп Z
RED
почвенный ВИ-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)A2 — 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)
Где:
DociFreezer
created by free version of
E = (2*(NIRA2-REDA2)+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, так как нас интересует испарение воды с поверхности земли (см. глава Водный индекс), приравняв NDWI<0 к 0, получим следующую формулу:
ШМ1%гоипа = (ЫБМ1 < 0)? 0, ЫБМ1
где: ЫБМ^гои^ — значение водного индекса в околоземном слое почвы, здесь примерно равный, значению транспирации растений
Для краткости NDWIground в дальнейшем обозначим как NDWI.
Если рассмотреть изображения обоих рассчитанных индексов в сравнении (Рис. 4), становится хорошо заметно, что в большей части изображения значения водного индекса (NDWI) инвертируют значения вегетативного индекса (NDVI). Этому явлению есть простое объяснение:
На площадях с наибольшим уровнем ВИ (а значит, теоретически, значительной прикорневой фитомассой растений) — будет более высокий уровень водяных паров за счёт транспирации, что, в свою очередь, приводит к падению водного индекса при поглощении отражённого ИК излучения водяными парами.
На площадях, где завышенные значения ВИ генерируются не биомассой растений, а высоким уровнем «шума» открытой почвы, нет высокого уровня транспирации и, следовательно, значения водного индекса так же высоки.
Рисунок 4. Сравнение изображение карты вегетативного индекса (слева) и
водного индекса (справа)
Получившиеся картина инверсии, независимо вычисленных индексов, приводит к идее учёта влияния значения транспирации (через водный индекс) на величину ВИ. Для этого попробуем использовать метод вычитания из «позитива» NDVI «негатива» NDWI и проанализируем полученную карту разности индексов (рис. 5)
Рисунок 5. Результат вычитания водного индекса из вегетативного индекса (слева), выделение отрицательных областей разности NDVI и NDWI
(справа)
Как видно на рисунке 8 (слева), разность вычисленных индексов напоминает NDVI — но в сильно контрастном исполнении. Смазанность фона ВИ теперь отсутствует, но зато в резком контрасте проявляются области с отрицательными значениями (рис. 5 справа). С точки зрения логики NDVI, получившаяся карта разности не может отражать проявление фотосинтетической активности растений.
Отрицательные значение разности индексов являются проявлением соотношения NDWI > предположительно, в местностях, с наиболее
деградировавшей в летний период вегетации растительностью. Связано это с тем, что при почти полном отсутствии фотосинтезирующей биомассы, значения ВИ будут низкими (даже с учётом «почвенного шума»), а водный индекс наоборот высоким — из-за низкого содержания водяного пара. Экспериментально установлено, что для исследуемой области подобная картина типична для середины вегетативного периода (июнь, июль, август) и не типична для начала и конца периода вегетации. Соответственно, использование NDWI, в чистом виде, в качестве транспирационной маски вегетативного индекса, не предоставляется возможным.
Для превращения водного индекса в корректирующую транспирационную маску для NDVI, можно попытаться применить эмпирические методы с последующей проверкой результатов в полевых условиях. Используем следующую логику: если значения водного индекса велики — попробуем их «приглушить». Наиболее оптимально в качестве подобного «глушителя» показала себя сумма NDVI и NDWI, а произведение данной суммы на значение водного индекса даёт приемлемую картину транспирационной маски в нужном диапазоне значений (рис. 6):
Ттазк = ЫБМ1 * (N07! + №М) где: Ттазк = вычисленная транспирационная маска
Рисунок 6. Сравнение суммы NDVI и NDWI (слева) с вычисленной транспирационной маской (справа)
Используем разность вегетативного индекса и транспирационной маски для вычисления скорректированного «контрастного» транспирационного ВИ:
ШВУ! = N071 — Ттазк
где: TNВУI = «контрастный» транспирационный ВИ Ттазк = вычисленная транспирационная маска
Проанализируем получившуюся карту транспирационного ВИ (рис. 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 «Multitemporal 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,
DociFreezer
created by free version of
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).
