CC BY
244
М.Ю. Катаев, А.А. Бекеров, П.В. Шалда. Анализ трендов временных рядов вегетационного индекса NDVI
81
УДК 528:581.5
М.Ю. Катаев, А.А. Бекеров, П.В. Шалда
Анализ трендов временных рядов вегетационного индекса NDVI
Рассматривается методика анализа трендов временных рядов вегетационного индекса NDVI (Normalized Difference Vegetation Index) за промежуток времени в несколько лет. Методика основывается на построении долговременного (несколько лет), годового и полугодового линейных трендов и последующего анализа их коэффициентов. Приводятся результаты оценки параметров трендов вегетационного индекса NDVI полученных по данным спутникового радиометра MODIS.
Ключевые слова: вегетационный индекс NDVI, тренды временных рядов. doi: 10.21293/1818-0442-2017-20-1-81-84
Задача выявления изменений в общем случае является областью хорошо изученной в статистике [1], обработке сигналов [2] и теории управления [3]. Тем не менее, большинство методов из данных областей недостаточно хорошо показывают себя при работе с массивными многовременными пространственно-временными наборами многоспектральных данных дистанционного зондирования Земли (ДЗЗ). Данный факт обусловлен ограничениями, такими как вычислительная сложность, и неспособностью использовать в полной мере эффект сезонности и пространственно-временной избыточности, которые присущи данным дистанционного зондирования Земли.
Предлагаемая методика поиска изменений, происходящих на поверхности Земли, основана на последовательном анализе коэффициентов линейных трендов различной длительности: несколько лет, годовой и полугодовой. Широкое распространение на практике имеет подход использования единого тренда по многолетним данным или в рамках одного года. Однако обнаружить таким способом можно лишь существенные изменения типов поверхности. Эта работа является продолжением работ [4, 5].
Постановка задачи
Данные дистанционного зондирования Земли для случая многоспектральных данных состоят из измеренного отраженного солнечного излучения от поверхности Земли или излучения типами поверхности в ИК-каналах. Естественно, что измеренные данные зависят от времени года, состояния атмосферы и поверхности Земли. Измеренные данные являются функциями спектрального коэффициента отражения, которые связаны с типами поверхности. Естественно, что измеряемые сигналы, таким образом, зависят от наличия того или иного типа поверхности в пятне обзора спутникового прибора. На основе нескольких спектральных каналов можно построить алгоритмы (расчет индексов) выделения типов поверхности. Одним из таких индексов является нормализованный вегетационный индекс МЭУ1.
Общий вид таких наборов спутниковых данных представлен на рис. 1, а, где х и у - размер спутникового изображения (число пикселей), - длина волны (/ = 1, ... Ы, N - число спектральных каналов)
и / - время измерения. На рис. 2, б представлен вид зависимости нормализованного разностного вегетационного индекса от времени, который отражает биомассу растительного покрова в определенной области по шкале от -1 до 1.
NDVI
MIl^lilJUHJyiyiLAl1
-jiftudiu^^
б t, дни
Рис. 1. Набор многоспектральных данных ДЗЗ, полученных в различное время: а - общий вид и б - их преобразование в нормализованный разностный вегетационный индекс
В настоящее время существуют данные с различными характеристиками (различное число спектральных каналов, пространственное и временное разрешение). В данной работе использовались данные, полученные со спутника TERRA от прибора MODIS. Данные этого спутника поступают для обработки ежедневно. Размер сцены одного снимка составляет 2340 на 2030 км, из которого вырезается исследуемая территория (район с центром в г. Том -ске, размером 500*1000 пикселей). Данные предоставляются в форматах HDF, который нами преобразуется в GeoTIFF, так как этот формат больше подходит для проведения анализа временных рядов.
y
t
x
n
а
Основной задачей данной работы является анализ временных рядов данных вегетационного индекса NDVI [4, 5], который рассчитывается на основе данных спутникового прибора MODIS. Проводимый анализ связан с оценкой трендов временных рядов NDVI, с целью обнаружения долговременных изменений (несколько лет), в течение года и полугодовых. Знание такой информации позволяет оценивать текущее состояние территории пикселя, опираясь на историю предыдущих состояний.
Анализ временных рядов
В работе используются данные обработки уровня L2 в виде продукта MOD09Q1. По данным этих двух каналов рассчитывался вегетационный индекс нормированной разности - Normalized Difference Vegetation Index (NDVI). Этот индекс является показателем, связанным с количеством растительности на исследуемой территории. Временные ряды представляют собой последовательность пикселей в различные моменты времени. Различие заключается в том, что во время наличия облачности значения вегетационного индекса существенно отличаются в меньшую сторону по отношению к моментам времени, когда в момент измерений облачности нет. Для получения равномерной по времени сетки значений вегетационного индекса, нами проводится для ранее полученных значений NDVI процедура нормализации [6]. Это делается для привлечения стандартных методов анализа временных рядов. Нами анализ временных рядов, главным образом, используется для определения изменений процесса развития вегетационного индекса во времени.
Существует два подхода, которые являются наиболее общими для анализа временных рядов, а именно: долговременный анализ рядов на основе исторических данных (ранее полученных) и анализ последовательности рядов в режиме реального времени. Для обнаружения изменений для первого подхода существуют следующие наиболее популярные методы: множественный регрессионный анализ и динамическая байесовская сеть.
На рис. 2 приведен исходный временной ряд вегетационного индекса NDVI за период времени 2009-2016 гг. и его нормализованное значение (темная гладкая кривая). Исходный ряд NDVI, как видно, является сильно изрезанным, что связано с тем, что в момент съемки спектрорадиометром исследуемой территории над ней находилась облачность. Хорошо видно, что нормализация приводит к существенно более гладкой форме временного ряда вегетационного индекса, что важно для проведения анализа и устойчивости получаемых результатов.
2009 2010 2011 2012 2013 2014 2015
Рис. 2. Временной ряд вегетационного индекса NDVI за 2009-2016 гг.
Исследование трендов временных рядов
Временной ряд вегетационного индекса можно разложить на отдельные составляющие: тренд, сезонная составляющая и случайные остатки (аномалии). Таким образом, значения вегетационного индекса N01 в момент времени t представляются как £(/) = &г(0 + + АS(t), (1)
где Str(t), Sseas(t), АS(t) - величина тренда N0'VI, сезонная составляющая значений N01 и отклонения значений вегетационного индекса (аномалия); t - номер дня измерений.
Разложение временного ряда N01 на компоненты осуществляется поэтапно. Вначале выделяется тренд в линейной форме: Str(t) = а + Ъ4, где а -средняя скорость изменения индекса; Ъ - константа. Для оценивания параметров линейного тренда используется метод наименьших квадратов. Далее из исходного ряда удаляется величина тренда и для остатка определяется сезонный ход. В данной работе проводится анализ линейного тренда для всего времени наблюдений (несколько лет) - долговременный тренд, годовой и полугодовой тренд.
Выявление изменений на временных данных с помощью линейных трендов
Одними из наиболее популярных методов выявления изменений на временных сериях являются методы с применением различной формы трендов (полиномиальные, нелинейные и др.). При исследовании временных серий с использованием данных дистанционного зондирования Земли с применением трендов нами предлагается выделить следующие подходы к оценке тренда по временным данным: полугодовые, годовые и долговременные.
На рис. 3, а-в представлены виды рассматриваемых трендов. Для их построения был применен метод линейной регрессии с использованием метода наименьших квадратов и в качестве неизвестных оценены значения коэффициентов модели а и Ъ.
i
2009 2010 2011 2012 2013 2014 2015
а
1
2009 2010 2011 2012 2013 2014 2015
б
1
в
Рис.3. Долговременный тренд (а), годовой тренд (в) и полугодовой тренд (в)
В случае с долговременным трендом значения коэффициентов а и Ъ имеют значения 0,302 и -0,03,
М.Ю. Катаев, А.А. Бекеров, П.В. Шалда. Анализ трендов временных рядов вегетационного индекса МБУ!
83
что свидетельствует об убывающем тренде (отрицательное значение коэффициента Ъ).
При использовании годового тренда, когда линейные тренды оцениваются за каждый год отдельно, появляется возможность более детально оценить изменения, произошедшие в каком-либо году временной серии. В табл. 1 представлены величины а и Ъ всех линейных функций годовых промежутков тренда.
Таблица 1 Значения коэффициентов a и b линейного годового тренда
Год a b
2009 0,21338 0,0343
2010 0,3384 -0,0217
2011 0,3304 -0,0293
2012 0,3052 -0,0423
2013 0,28978 -0,0437
2014 0,2702 -0,0538
2015 0,20163 0,0004
Из табл. 1 видно, что коэффициент Ъ на протяжении промежутка времени, начиная с 2010 г. и заканчивая 2014 г., имел тенденцию к убыванию, что хорошо видно на рис. 3, б. Эти значения хорошо коррелируют с тем фактом, что отрицательное изменение вегетационного индекса совпадает с началом строительства на исследуемой территории. Стройка связана с уменьшением растительности на территории, что ведет к убыванию индекса N01.
Еще одним вариантом построения тенденции временной серии является тренд с полугодовыми промежутками времени, который позволяет еще более детально оценить изменения, протекающие на изучаемой территории (см. рис. 3, в). В табл. 2 представлены величины а и Ъ всех линейных функций полугодовых промежутков тренда.
Таблица 2 Значения чисел а и Ь линейных функций полугодовых промежутков тренда_
Год Полугодие a b
2009 1 -0,117 0,004
2 0,546 -0,003
2010 1 0,009 0,003
2 0,474 -0,002
2011 1 0,0996 0,003
2 0,556 -0,003
2012 1 0,011 0,003
2 0,524 -0,003
2013 1 -0,009 0,002
2 0,371 -0,002
2014 1 0,0936 0,001
2 0,349 -0,002
2015 1 0,005 0,002
2 0,359 -0,002
Из рис. 3, в видно, что при построении данного типа тренда образуются треугольники различных типов относительно положения вершины, амплитуды и симметрии сторон. Наблюдая типы треугольников за весь период наблюдений, можно отметить,
что по ним можно оценить основную тенденцию происходящих изменений. Беря в расчёт 2009 г. как опорный, можно сделать вывод, что с каждым последующим годом углы, противолежащие основанию, имеют тенденцию к увеличению, что в свою очередь свидетельствует об уменьшении растительной биомассы в данной исследуемой области. Таким образом, результаты изучения трендов различной длительности позволяет взглянуть на состояние растительности в исследуемом пикселе более детально.
Заключение
С увеличением доступности данных дистанционного зондирования для широкого круга исследователей методы их обработки и анализа стали стремительно развиваться. Длительный период измерений позволяет решать такую задачу, как выявление изменений. Обнаружение и классификация изменений типов поверхности Земли является очень важной задачей, т.к. в целом может оценивать как гло -бальные, так и пространственно локальные изменения для определенного участка земной поверхности. В предлагаемой вниманию статье предлагаются методика и результаты оценки изменений на основе анализа коэффициентов линейных трендов временных рядов вегетационного индекса NDVI. На основе обработки большого массива данных для территории, где проводились строительные работы, продемонстрирована эффективность предлагаемого подхода. Визуальный анализ и величины коэффициентов линейного тренда отчетливо показывают момент наступления изменения (2010 г.), после чего наблюдается тенденция к спаду вегетационного индекса, что свидетельствует о значительном уменьшении растительной биомассы. Через некоторое время, когда строительные работы завершены, начинается озеленение территории, что и наблюдается нами в 2015 г., когда происходит рост (пусть и незначительный) тренда.
В данной статье были рассмотрены некоторые простые методы анализа временных серий для данных дистанционного зондирования Земли. При этом в самих алгоритмах были использованы данные, которые прошли обработку на сглаживание и удаления шума. Представленные методы, несмотря на простоту, продемонстрировали состоятельность и наглядность. Стоит отметить, что при работе с временными сериями, которые являются данными достаточно массивными, например, в данных примерах использовался набор данных, состоящий из двух с лишним тысяч снимков, очень важно, чтобы используемый алгоритм обладал хорошим быстродействием и простотой, что минимизирует вычислительные ошибки.
Работа выполнена в Центре космического мониторинга Земли ТУСУРа.
Литература
1. Inclan C. Use of cumulative sums of squares for retrospective detection of changes of variance / C. Inclan, G.C. Tiao // Journal of the American Statistical Association. -1994. - № 89т(427). - РР. 913-923,
2. Gustafsson F. Adaptive Filtering and Change Detection. - John Wiley & Sons, 2000. - 233 р.
3. Lai T.L. Sequential changepoint detection in quality control and dynamical systems // Journal of the Royal Statistical Society. Series B (Methodological). - 1995. - № 57 (4). -PP. 613-658.
4. Катаев М.Ю. Обнаружение экологических изменений природной среды по данным спутниковых измерений / М.Ю. Катаев, А.А. Бекеров // Оптика атмосферы и океана. - 2014. - № 7. - С. 652-656.
5. Катаев М.Ю. Интернет-информационная система накопления, обработки и анализа спутниковых данных MODIS / М.Ю. Катаев, А.А. Бекеров, А.К. Лукьянов // Доклады ТУСУРа. - 2015. - № 1. - С. 93-99.
6. Катаев М.Ю. Методика выравнивания временных рядов вегетационного индекса NDVI, полученных по данным спектрорадиометра MODIS / М.Ю. Катаев, А.А. Бекеров, А.К. Лукьянов // Доклады ТУСУРа. - 2016. - № 1. -С. 35-39.
Катаев Михаил Юрьевич
Д-р техн. наук, профессор каф. автоматизированных систем управления (АСУ) ТУСУРа, профессор Юргинского технологического института,
филиала Национального исследовательского Томского политехнического университета Тел.: (382-2) 70-15-36, 8-960-975-27-85 Эл. почта: kataev.m@sibmail.com
Бекеров Артур Александрович
Аспирант Института мониторинга климатоэкологических
систем СО РАН, Томск
Тел.: (382-2) 70-15-36
Эл. почта: artur@bekerov.ru
Шалда Павел Валерьевич
Магистрант каф. АСУ ТУСУРа
Тел.: (382-2) 70-15-36
Эл. почта: pavel.shalda93@gmail.com
Kataev M.Yu., Bekerov A.A., Shalda P.V. Trend analysis for time series of NDVI
The trend analysis of time series of the Normalized Difference Vegetation Index (NDVI) for the period of few years is reviewed. The method is based on the establishing long-term trend (several years), annual and semi-annual linear trends and further assessment of their ratios. The results of the parameter estimates for NDVI trends derived from satellite radiometer MODIS are described.
Keywords: vegetation index, NDVI, trends, time series.
