Алгоритм оценки долговременных вариаций MODIS NDVI

Важнейшей информацией, которая позволяет исследовать происходящие на поверхности Земли процессы, являются данные дистанционного зондирования Земли [1]. С их помощью решается обширный класс задач, связанных с оценкой изменения климата, мониторингом стихийных бедствий, изучением динамики геосистем и т.д. [2, 3]. При этом исследование вариаций характеристик растительного покрова приобретает особое значение на фоне происходящих климатических изменений [4]. Для отслеживания динамики растительного покрова широко используются вегетационные индексы, среди которых наиболее распространен нормализованный относительный вегетационный индекс NDVI (Normalized Difference Vegetation Index), основывающийся на отражении и поглощении растительностью электромагнитного излучения в красном и ближнем инфракрасном диапазонах соответственно.

В настоящее время широкое распространение получили съемочные системы TM, ETM+ и OLI (спутников серии Landsat), AVHHR (NOAA), MODIS (Terra и Aqua), данные которых находятся в свободном доступе. Снимки со спутников Landsat обладают более высоким пространственным разрешением по сравнению с AVHHR и MODIS, но имеют частоту съемки один раз в 16 сут и часто оказываются неприменимы для анализа из-за наличия облачности. Для задач непрерывного мониторинга растительного покрова оптимальным представляется использование данных MODIS высокого временного разрешения (частота съемки может достигать до четырех раз в сут) и более высоким пространственным разрешением по сравнению с AVHHR c широким охватом съемки. При этом необходима разработка алгоритмов для обработки большого объема данных спутниковой информации.

Одним из наиболее используемых средств анализа временных рядов вегетационных индексов является программа TIMESAT [5]. Данная программа позволяет аппроксимировать временные ряды с целью определения фенологических параметров растительности с использованием асимметричной гауссовой функции, двойной логистической функции, фильтра Савицкого-Голея. TIMESAT обладает развитым инструментарием для удаления пиков, сглаживания данных, использования информации о типах ландшафтов. Недостатки TIMESAT: необходимость конвертации данных в бинарный формат с использованием стороннего программного обеспечения, отсутствие возможности полиномиальной аппроксимации трендовой компоненты временного ряда и анализа вклада различных фаз вегетационного периода в общий тренд.

Цель данной работы является создание алгоритма оценки многолетних изменений состояния растительности с использованием временных рядов вегетационного индекса MODIS NDVI. В результате разработан алгоритм обработки рядов NDVI высокого временного разрешения, проведены расчеты трендов NDVI для обширной территории Западного Забайкалья, представлены примеры изменения растительного покрова с использованием снимков сверхвысокого пространственного разрешения.

Материалы и методы

В работе использованы данные спектрорадиометра MODIS, установленного на спутниках Terra и Aqua. Данные доступны для загрузки из портала NASA EOSDIS. Исходными данными послужили временные ряды NDVI тематического продукта MOD13Q1, который помимо нормализованного разностного вегетационного индекса NDVI содержит усовершенствованный вегетационный индекс EVI [3]. Продукт представлен растровыми 16-дневными композитами с пространственным разрешением 250 м в синусоидальной проекции. Для создания композита используется метод максимального значения MVC (Maximum Value Compositing). Композит формируется из максимальных значений NDVI за 16 дней, таким образом исключая влияние атмосферной составляющей. Кроме того, для сокращения оптической дистанции выбирают данные со спутника, находящегося ближе надиру, что также уменьшает атмосферное влияние.

На рис. 1. представлена блок-схема разработанного алгоритма.

Значения NDVI из тематического продукта MOD13Q1 записывают в куб данных размерностью 4800×4800×N, где первые две координаты являются пространственными, третья – временной, где N – количество композитов. Данные содержатся в иерархическом hdf-файле и представлены в виде NDVI×10000.

Рис. 1. Блок-схема алгоритма обработки данных

Fig. 1. Flowchart of the data processing algorithm

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

Для избавления от случайной составляющей временных рядов проводится их сглаживание с использованием фильтра Савицкого-Голея, который часто применяется для снижения уровня шумов [6]. Порядок полинома равен 2, ширина окна равна количеству композитов в году.

Чтобы исключить влияние на тренд значений вегетационного индекса открытой почвы и снега, временной ряд ограничивается значениями NDVI для периода вегетации [6]. За начало (SOS) и конец (EOS) вегетационного периода принято значение, когда NVDI становится соответственно выше и ниже порогового значения. В работах [4, 7] в качестве пороговых значений вегетационного периода приняты значения NDVI, равные 0,14 и 0,1. Поскольку для разных сенсоров пороговое значение NDVI может отличаться, то в алгоритме оно задается пользователем. Далее определяется длительность вегетационного периода (LOS), постоянная для каждого пикселя. За начало периода вегетации выбирается наиболее позднее значение, а за конец – наиболее раннее, т.е. длительность периода выбирается минимальной (рис. 2).

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

На рис. 3 представлен пример реализации алгоритма для одного временного ряда NDVI.

Для каждого временного ряда проводится оценка статистической значимости трендов с использованием критерия Фишера. В качестве нулевой гипотезы принимаем, что тренд равен 0. Значимым считаем тренд с 95%-ным уровнем (p < 0,05). Пользователь имеет возможность задавать уровень значимости.

Рис. 2. Схема определения вегетационного периода

Fig. 2. Scheme of the definition of growing season

Рис. 3. Пример обработки временного ряда NDVI. Черным цветом показаны исходные данные, красным – отфильтрованный ряд, синим – нелинейный тренд, коричневым – линейный тренд

Fig. 3. Example processing of NDVI time series. Black color shows the initial data, red – filtered series, blue – non-linear trend, brown – linear trend

Реализована возможность оценки вклада фаз вегетационного периода в долговременную динамику NDVI. Так, определяются тренды максимума вегетации POS, суммы значений в фазе роста от SOS до POS и в фазе окончания вегетации от POS до EOS.

Выходные данные реализации алгоритма представляются в растровом изображении tif-формата.

Результаты и обсуждение

Разработанный алгоритм апробирован для территории Западного Забайкалья. Входными данными послужили композиты NDVI за период с 2000 по 2016 г. В результате попиксельной обработки получено более 46 млн регрессионных уравнений и оценены их коэффициенты наклона (тренд).

Построена карта пространственного распределения трендов NDVI в Западном Забайкалье (рис. 4). Зоны с отрицательным трендом NDVI, характеризующим устойчивое снижение фотосинтетически активной биомассы, приурочены к степным растительным сообществам, которые в основном расположены в межгорных котловинах. Территории с наибольшей площадью отрицательных трендов NDVI располагаются в долинах рр. Селенга, Уда, Хилок и Баргузин. Зоны с положительным трендом вегетационного индекса NDVI приурочены к лесным ландшафтам.

В результате анализа наиболее высоких и низких значений линейных трендов NDVI выявлены области с наиболее интенсивными изменениями ландшафтов. В результате визуальной верификации трансформации ландшафтов с помощью спутниковых изображений сверхвысокого пространственного разрешения выявила, что изменения, в основном, обусловлены антропогенным влиянием. На рис. 5 представлена серия снимков сервиса Google Earth, на которых – 65 –

Рис. 4. Карта линейных трендов NDVI в Забайкалье с 2000 по 2016 гг. Зеленым и красным цветами обозначены статистически значимые (p < 0,05) положительные и отрицательные тренды соответственно

Fig. 4. Map of linear trends of NDVI in Transbaikalia from 2000 to 2016. Green and red colors indicate significant (p < 0.05) positive and negative trends

в)

Рис. 5. Изменение ландшафтов в долине р. Топка (Бичурский район, Республика Бурятия): а – 2005 г.; б – 2010 г.; в – 2015 г.; 1 – участок вырубки леса; 2 – участок лесовосстановления; 3 – разработка Окино-Ключевского угольного карьера

б)

Fig. 5. Changes of landscapes in the valley of the river Topka: a – 2005; b – 2010; c – 2015; 1 – the site of deforestation; 2 – the site of reforestation; 3 – development of the Okino-Klyuchevskiy coal mine

а)

б)

в)

Рис. 6. Тренды NDVI с 2000 по 2015 гг.: а – участок вырубки леса; б – участок лесовосстановления; в – разработка Окино-Ключевского угольного карьера

Fig. 6. NDVI Trends from 2000 to 2015: a – a site of deforestation; b – a site of reforestation; c – development of the Okino-Klyuchevskiy coal mine зафиксированы рубка и восстановление леса, а также разработка угольного карьера. Показано, что на ранее вырубленных участках леса процесс лесовосстановления до фоновых значений NDVI продолжался около 12 лет (с 2002 по 2014 гг.). На рис. 6 изображены соответствующие временные ряды NDVI.

В результате автоматизированной обработки временной серии композитов NDVI обнаружены 283 участка с экстремально положительными трендами NDVI (области лесовосстановления) общей площадью 145,3 км2, 470 участков с экстремально негативными трендами охватывают 919,2 км2. Разработанную методику оценки ландшафтных изменений необходимо уточнить с учетом региональных физико-географических особенностей для обширной территории Забайкалья.

На рис. 6 представлены графики NDVI, характерные для различных случаев изменения ландшафтов. Вырубка леса (рис. 6 а ) проявляется в резком снижении NDVI с 0,7 (хвойные леса) до 0,4 (степная растительность). Лесовосстановление характеризуется постепенным ростом NDVI (рис. 6 б ). Разработка угольного карьера проявляется в снижении NDVI до значений 0,15-0,2, характерных для открытых почв (рис. 6 в ).

Выводы

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

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

Исследование выполнено при финансовой поддержке РФФИ в рамках проекта № 17-05-01059 и в рамках государственного задания БИП СО РАН.