Методика сглаживания временного ряда вегетационного индекса NDVI Текст научной статьи по специальности «Математика»

УДК 528:581.5

М.Ю. Катаев, А.А. Бекеров, Д.В. Медвецкий

Методика сглаживания временного ряда вегетационного индекса NDVI

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

Ключевые слова: вегетационный индекс NDVI, временной ряд, сглаживание. doi: 10.21293/1818-0442-2017-20-1-85-88

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

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

Построение модели временного ряда вегетационного индекса NDVI в настоящее время связано с применением алгоритмов сглаживания различными методами. Эти методы позволяют из сильноизрезанной структуры временного ряда NDVI получить гладкую форму, удобную для анализа. Однако получаемая форма теряет, как правило, особенности изменения NDVI в различные фазы фенологического периода за счет сглаживания. Нами в данной статье приводится простая методика, удобная для автома-

тизации, позволяющая построить гладкую форму временного ряда NDVI в течение года и в основной части оставить особенности фенологического развития растительности.

Постановка задачи

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

Одной из проблем сглаживания временных рядов вегетационного индекса [5-9] является выявление ложных провалов и пиков, которые связаны с иными процессами, чем влияние облачного покрова. Чаще всего ложные провалы и пики возникают вследствие изменения таких метеорологических параметров, как температура (воздуха и поверхности) и влажность. Устранение изрезанности временного ряда эффективно устраняется временным сглаживанием, когда строят 8- или 16-дневные композиты (например, M0D09Q1 и MOD13Q1 для спутникового радиометра Terra MODIS). Результатом построения композитов являются временные ряды, вид которых показан на рис. 1.

На рис. 1 показаны временные ряды вегетационного индекса NDVI, для одного из исследуемых пикселей в районе г. Томска. Видно, что исходный временной ряд NDVI (однодневный), показанный на рис. 1, а, является весьма изрезанным, а композитный, приведенный на рис. 1, б, уже не имеет множества провалов и пиков. Композитный временной ряд вегетационного индекса является существенно ме-

нее изрезанным, однако потерявшим особенности, присущие однодневным измерениям, если бы каждое измерение было получено в безоблачном случае.

30 60 30 120 150

210 210 270

Номер И ОДСрСЕШЯ

а

0-9 -

> -

и

11,7 -

и

и 3

32 0,5-

Г

0.4 -

О

1

и <и-

ш

0,1 -

0,0 -

од-

-|—'—г

г 4

—I—1—I—I—I—I—|—1—I—I—|—1—I' I I

1(1 12 14 16 1Я 20 22 24

Рис.

] 1омер измерения

б

1. Временной ряд в течение нескольких лет: а - исходный временной ряд и б - 16-дневный композитный

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

В статье [10] приведены результаты сглаживания методом скользящего среднего и вписывания функции Гаусса во временной ряд МЭУ1, которые говорят, что полученные результаты не позволяют применять их на практике ввиду существенного отличия от исходного ряда. Поэтому вышеприведенные результаты заставили нас разработать методику сглаживания временного ряда МЭУ1, эффективную по скоростным, точностным характеристикам и удобную для автоматизации.

Описание методики сглаживания

Идею сглаживания временного ряда МЭУ1 проиллюстрируем рис. 2, на котором приведена сезонная динамика значений вегетационного индекса МЭУ1 (январь).

На рис. 2 видно, что в зимний период вегетационный индекс МБУ1 имеет низкое значение и сильно меняется в диапазоне [-0,06, 0] и [0, +0,06]. Так как реальное значение МЭУ1 для данного пикселя определяется отражательными свойствами снега и зданий (территория г. Томска), то можно предположить, что временные вариации МЭУ1 не должны быть такими значительными и лежать в положительной области.

0,07 0,06 0,05 0,04 0,03 0,02 0,01 0,00 -0,01 -0,02 -0,03 -0,04 -0,05 -0,06 -0,07 -0,08 -0,09 -0,10

1 ' I ' 1 ' I ' 1 ' I ' Г ' 1 ■ Г ' 1 ■ Г ' 1 К ¡0 12 14 16 1$ 20 22 24 26 28 30

Помер измерения

Рис. 2. Первые 30 дней (январь) вегетационного индекса (см. рис. 1, а)

Будем считать, что в течение года измеряется N = 365 значений вегетационного индекса МЭУ1 -£(/) с фиксированным шагом ( = 1, ..., Щ. Для каждого временного шага получается матрица значений, покрывающая область наблюдения (район г. Томска). Для примера рассмотрим временные изменения МЭУ1 в одной и той же точке (пикселе), для которой типы поверхности не изменялись за общее время измерений.

Методика сглаживания связана с идеей медианной фильтрации шума, а именно, применения алгоритма скользящей медианы. Вычисление значения скользящей медианы в t связано с анализом значений временного ряда в интервале [ - к, t + к], где (2к + 1) - ширина окна (или скользящая медиана). Основным преимуществом медианного фильтра перед многими другими сглаживающими фильтрами является устойчивость к наличию во временном ряде пропусков и выбросов. Одним из минусов медианного фильтра является невозможность вычислений, если точки ( + к) или ( - к) попадают на начальные или конечные точки временного ряда. Решение этой проблемы заключается в выборе на этих временных интервалах фильтров меньшего порядка (минимальная длина к = 3). В текущем окне ищется максимальное значение и им заменяется значение в точке t, далее окно смещается на один шаг, и операция повторяется снова. Так как во временном ряде МБУ1 могут быть провалы (см. рис. 1) за счет облачности до 10-12 дней, максимальная длина фильт-

ра нами выбрана равной к = 15. В этом случае после прохождения фильтром всего временного ряда получается новый массив значений Sn(t), который сглаживается далее простым сглаживающим фильтром (скользящее среднее), формула которого выглядит следующим образом:

Sm(t) = Sn(t + 1)/4 + Sn(t)/2 + Sn(t - 1)/4, где Sm(t) - сглаженный временной ряд вегетационного индекса.

Полученные результаты

Предложенная нами методика была апробирована при сглаживании данных, полученных со спутникового прибора МОБШ за период 10 лет и территории в районе г. Томска размером (500*1000 пикселей). На рис. 3 показан исходный и сглаженный, по предлагаемой методике ряд значений МЭУ1, для одного из пикселей исследуемой территории.

U

Ц>-|

0.9 -

> о.к -

а

£ (1J -

о

к

о |),й -

0,5-

0,4-

Ch

ч 03-

Si

и D IV-

е

0,1 -

од -

-0,1 -

Рис.

0 30 60 90 124 150 180 210 240 270 300 330 360 Помер измерения

3. Исходный и сглаженный с учетом предлагаемой методики ряд NDVI

Из рис. 3 видно, что предлагаемая методика временного ряда МЭУ1 позволяет получить достаточно гладкую форму, удобную для решения большинства основных практических задач. Отметим, что в основу гладкой формы ряда МЭУ1 легли все исходные значения ряда (ежедневные), что при сравнении с 16-дневными композитами носит более объективный характер. Кроме того, методика позволяет быстро обрабатывать большие объемы информации. Можно отметить, что основной проблемой получаемого сглаженного ряда являются выбросы, которые приводят к появлению локальных максимумов. В следующих статьях будет рассмотрена закономерность данных проявлений, когда каждая временная точка будет связываться с вариациями температуры и влажности.

Заключение

В данной статье рассмотрена простая вычислительная методика сглаживания временных данных вегетационного индекса МЭУ1. Основой предлагаемого подхода является использование алгоритма

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

Литература

1. Михайлов Н.Н. Использование временных рядов вегетационного индекса NDVI для мониторинга растительного покрова степной зоны Западной Сибири / Н.Н. Михайлов, Л.А. Михайлова, Н.Ф. Харламова и др. // Научные ведомости БелГУ. Сер. Естественные науки. -2010. - №15, вып. 12. - С. 25-33.

2. Спивак Л.Ф. Детектирование аномальных значений временных рядов вегетационных индексов // Современные проблемы дистанционного зондирования Земли из космоса. - 2014. - Т. 11, № 3. - С. 193-199.

3. Лисецкий Ф.Н. Использование космического мониторинга для изучения элементов водного баланса в целях адаптивного землеустройства агроландшафтов / Ф.Н. Лисецкий, Т.Н. Ковалева // Научные ведомости Бел-ГУ. Сер. Естественные науки. - 2011. - № 21, вып. 17. -С. 108-118.

4. Майорова В.И. Контроль состояния сельскохозяйственных полей на основе прогнозирования динамики индекса NDVI по данным космической мультиспектраль-ной и гиперспектральной съёмки / В.И. Майорова,

A.М. Банников, Д.А. Гришко и др. // Наука и образование: научное издание МГТУ им. Н.Э. Баумана. - 2013. -

B. 07. - С. 199-228.

5. Плотников Д.Е. Восстановление временных рядов данных дистанционных измерений методом полиномиальной аппроксимации в скользящем окне переменного размера / Д.Е. Плотников, Т. С. Миклашевич, С. А. Барта-лев // Современные проблемы дистанционного зондирования Земли из космоса. - 2014. - Т. 11, № 2. - С. 103-110.

6. Тронин А. А. Анализ длинных рядов вегетационного индекса территории РФ и регионов / А.А. Тронин, А.В. Киселев // Современные проблемы дистанционного зондирования Земли из космоса. - 2012. - Т. 9, № 1. -

C. 108-113.

7. Bradley B.A. A curve fitting procedure to derive inter-annual phenologies from time series of noisy satellite NDVI data / B.A. Bradley, R.W. Jacob, J.F. Hermance, J.F. Mustard // Remote Sensing of Environment. - 2007. - Vol. 106. -PP. 137-145.

8. Hird J.N. Noise reduction of NDVI time series: An empirical comparison of selected techniques / J.N. Hird, G.J. McDermid // Remote Sensing of Environment. - 2009. -Vol. 113. - PP. 248-258.

9. Jonsson P. TIMESAT - a program for analyzing time-series of satellite sensor data / P. Jonsson, L. Eklundh // Computers & Geosciences. - 2004. - Vol. 30. - PP. 833-845.

10. Катаев М.Ю. Методика выравнивания временных рядов вегетационного индекса NDVI, полученных по данным спектрорадиометра MODIS / М.Ю. Катаев, А.А. Бекеров, А.К. Лукьянов // Доклады ТУСУРа. - 2016. - № 1. -С. 35-39.

Медвецкий Дмитрий Владимирович

Магистрант каф. АСУ ТУСУРа

Тел.: (382-2) 70-15-36

Эл. почта: prodemon2012@gmail.com

Катаев Михаил Юрьевич

Д-р техн. наук, профессор каф. автоматизированных

систем управления (АСУ) ТУСУРа,

профессор Юргинского технологического института,

филиала Национального исследовательского

Томского политехнического университета

Тел.: (382-2) 70-15-36, 8-960-975-27-85

Эл. почта: kataev.m@sibmail.com

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

Аспирант Института мониторинга климатоэкологических

систем СО РАН, Томск

Тел.: (382-2) 70-15-36

Эл. почта: artur@bekerov.ru

Kataev M.Yu., Bekerov A.A., Medveckiy D.V. Method to smooth time series of NDVI

The method to smooth a time series of the Normalized Difference Vegetation Index (NDVI) for the period of time of one year is investigated. The main objective of the technique is to automate and to increase the speed and the efficiency in the presence of outliers. The authors describe the method and the results of its application to the smoothing of the NDVI derived from satellite radiometer MODIS.

Keywords: vegetation index, NDVI, time series, smoothing.