CC BY
204
УДК 519.254
М.Ю. Катаев, А.А. Бекеров, А.К. Лукьянов
Методика выравнивания временных рядов вегетационного индекса NDVI, полученных по данным спектрорадиометра MODIS
Приводится описание методики, позволяющей выравнивать временные ряды вегетационного индекса NDVI (Normalized Difference Vegetation Index) для пространственного разрешения 250 м и ежедневных измерений. Эта проблема возникает при необходимости проводить анализ динамики вегетационного индекса при решении конкретных задач сельского и лесного хозяйства, МЧС и др. Методика базируется на коррекции значений вегетационного индекса в зависимости от параметров облачной маски и величины температуры. Приводится обработка реальных данных спутниковых измерений в районе г. Томска.
Ключевые слова: спектрорадиометр MODIS, временные ряды, вегетационный индекс, методика выравнивания.
doi: 10.21293/1818-0442-2016-19-1-35-39
В настоящее время спутниковые данные наблюдения являются важнейшим глобальным и периодическим источником информации, необходимой для изучения атмосферы и поверхности Земли. С каждым годом увеличивается число спутниковых систем, обладающих разнообразными характеристиками (пространственными, временными, спектральными), что порождает значительный рост получаемых данных и приводит к необходимости развивать соответствующие методы хранения, обработки и анализа, которые должны обеспечивать эффективную и быструю работу со спутниковыми данными. Особое внимание уделяется развитию программно-информационной составляющей, позволяющей пользователям искать, получать, обрабатывать и визуализировать спутниковую информацию.
Решение задачи оценки состояния окружающей среды, претерпевающей изменения в зависимости от влияния естественных и антропогенных факторов, невозможно без наличия достоверной информации об ее параметрах. В настоящее время не существует наземной системы мониторинга, которая обеспечила бы регулярное, с широким пространственным охватом получение параметров окружающей среды. Единственно возможной является космическая система мониторинга окружающей среды, которая обеспечивает большую обзорность, оперативность и регулярность получения информации. Получаемые из космоса данные необходимо накапливать, обрабатывать и анализировать. В результате обработки появляется возможность получать информацию об основных параметрах атмосферы и поверхности Земли. Изучая пространственно-временное поведение параметров, описывающих состояние поверхности Земли (температура поверхности, коэффициент отражения, вегетационный индекс и др.), можно извлекать массу связанной с этими параметрами информации. Такой информацией, например, может быть изменение, которое связано с отличием типичного поведения параметра во времени или по пространству (примером является пожар).
Получаемые с помощью разнообразных спутниковых систем мониторинговые данные о поверх-
ности Земли связаны с процессом непрерывного сбора информации для конкретной территории. Существующие спутниковые системы, такие как MODIS [http://modis.gsfc.nasa.gov] и Landsat [http://landsat.gsfc.nasa.gov], уже более десяти лет регулярно поставляют спутниковую информацию, обработка которой дает возможность оценивать параметры атмосферы (аэрозоль, облачность, осадки, влажность и др.) и поверхности Земли (вегетационные индексы, отражательную способность, температуру и др.).
Важнейшей информацией, которая позволяет исследовать происходящие на поверхности Земли процессы, являются вегетационные индексы [1-3]. Классификация величин индексов и их взаимосвязей позволяет четко выделять объекты на Земле, которые относятся к разным типам поверхности. Классификация типов поверхности, поиск изменений необходимы в различных технологических, научных или экологических приложениях, связанных с контролем состояния лесных массивов, пахотных земель, строительных территорий, оценкой площади выгоревшего леса и др. Однако решение этих задач опирается на один важный минус спутниковых оптических измерений в видимой и инфракрасной областях спектра, который связан с наличием мешающего фактора облачности и свойств прозрачности атмосферы. При этом отметим, что в таком случае сами характеристики поверхности чаще всего не меняются, а меняется значение измеренных сигналов и значит, величина вегетационного индекса, что подлежит коррекции.
Данная статья является развитием работ [4, 5] и содержит методику, которая предназначена для выравнивания рядов наблюдений вегетационного индекса, полученных по данным спектрорадиометра MODIS.
Постановка задачи
Измеряемое на борту спутника отраженное солнечное и уходящее инфракрасное излучение зависит от свойств взаимодействия с объектами земной поверхности и атмосферы [1, 5]. С формализованной точки зрения регистрируются мгновенные
значения яркости I(x, у, X, t) для матрицы, состоящей из (i, j) строк и столбцов (имеющих пространственные координаты широта (у) - долгота (x)). Каждое такое изображение получено для набора спектральных каналов c длиной волны X (зависит от типа спутникового прибора) при заданных значениях ф зенитного и азимутального 6 углов наблюдения, ф8 зенитного и азимутального 6s угла Солнца, конкретного состояния атмосферы, которое обычно определяется ее прозрачностью, коэффициентами отражения типов поверхности и рельефом. Прозрачность атмосферы определяется наличием аэрозольных составляющих, приземной дымки, перистой облачности и облаков.
Нами в работе используются данные измерений спутникового прибора MODIS (Moderate Resolution Imaging Spectroradiometer. Прибор MODIS - сканирующий спектрорадиометр среднего разрешения), предназначенный для глобального мониторинга из космоса характеристик океана, суши и атмосферы. MODIS состоит из двух приборов, один из которых (MODIS-N) снимает в надир, а ось съемки другого (MODIS-T) может быть отклонена. 36 спектральных каналов MODIS охватывают диапазон с длинами волн от 0,4 до 14,4 мкм. Съемка в двух первых спектральных каналах (620-670 и 841-876 нм) ведется с пространственным разрешением 250 м, в пяти спектральных каналах видимого и ближнего инфракрасного диапазона - с разрешением 500 м, а в остальных (диапазон от 0,4 до 14,4 мкм) - 1000 м. Радиометрическое разрешение исходных снимков составляет 12 бит. Траектория движения носителя (высота орбиты 705 км) и угол обзора системы позволяют охватывать измерениями полосу обзора шириной 2330 км. Это позволяет прибору MODIS за сутки получать изображение практически всей поверхности Земли, за исключением отдельных участков в районе полюсов.
Наличие такого набора спутниковых данных является одним из важнейших факторов, определяющих последующие методы и алгоритмы обработки космических изображений. Обработка возможна при использовании двух принципиально разных подходов: первый предполагает использование изображений, полученных в некий фиксированный момент времени, а второй основан на составлении рядов разновременных изображений. Для первого случая применяются методы визуальной интерпретации или автоматической классификации, а для второго используются математические подходы обработки и анализа временных рядов.
Временные ряды измерений I(X, t, x, у) для каждого пространственно локализованного пикселя с координатами (x, у) и спектрального канала X дают возможность определять динамические закономерности поведения параметров атмосферы или поверхности Земли. Если остановиться на задачах, связанных с изучением поведения параметров поверхности Земли, то в этом случае удобно использовать вегетационные индексы NDVI. Данный индекс
удобен для идентификации типа отражающей поверхности. NDVI (Normalized Difference Vegetation Index) - нормализованный относительный индекс растительности, количественный показатель фото-синтетически активной биомассы. Один из самых распространенных и используемых индексов для решения задач, использующих количественные оценки растительного покрова. Идея расчета NDVI основана на спектральной зависимости коэффициента отражения растений, которая в красной области спектра содержит максимум поглощения солнечной радиации хлорофиллом растений, а в инфракрасной области содержит область максимального отражения типом и листовой структурой растений. Отношение спутниковых измерений, полученных в этих двух спектральных областях, позволяет четко отделять растительные объекты от прочих. Рассчитывается вегетационный индекс с использованием данных по отражению на двух каналах: K 1 - данные измерений в первом (0,620-0,670 нм) и K 2 - (0,841-0,876 нм) каналах прибора MODIS по формуле [6]
NDVI = (Ki - K2)/(Ki + K2), (1)
где K 1, K 2 - измеренные значения прибором MODIS в соответствующих каналах.
Зимой значения NDVI минимальны, а летом они максимальны (рис. 1). Вегетационный индекс обладает естественной межгодовой периодичностью, а значит, временной ряд содержит составляющие: трендовую, сезонную и случайную. Вариации значений вегетационного индекса в каждый конкретный день года могут быть значительными, поскольку они обусловливаются многими естественными причинами, такими как облачность, дождь, ветер, температура, прозрачность атмосферы и др. На рис. 1 приведены значения вегетационного индекса для пикселя в районе г. Томска для измерений MODIS в течение двух лет.
Вегетационный индекс NDVI 0.6 -
~1—•—Г-1—Г-1—г 0 60 120 180 240 300 360 «О 480 540 Г>00 660 720
Номер дня
Рис. 1. Значения временного хода (два года) вегетационного индекса NDVI в районе г. Томска
Наличие шумов во временных рядах NDVI, вызванных влиянием облачности и различиями в геометрических условиях освещения и наблюдения поверхности, требует проведения предварительной
обработки спутниковых данных - выравнивания. Необходимость этого шага связана с тем, что статистические характеристики временных рядов значительно отличаются от гладких функций, какими и должны быть функции МБУ1, связанные с фенологическим циклом растений. Поэтому решение практических задач, например сельского хозяйства, где необходимо знать дату начала фенологического цикла и другие показатели, определяемые с точностью определения параметров годового цикла изменений кривой МЭУ1 (годовой максимум, начало и завершение фенологического цикла и др.). Для задач лесного хозяйства возникают те же проблемы, которые связаны с параметрами кривой МЭУ1, когда необходимо оценить изменение состояния лесного массива, при помощи статистических характеристик ряда: тренд, вид и параметры кривой МЭУ1 (максимум, ширина и др.).
Предлагаемая методика
Выравнивание выполняется для устранения влияния атмосферы за счет рассеяния и поглощения газами, ослабления излучения аэрозолем и облачности. Изменение измеренных сигналов в дни, когда наблюдается облачность, можно учесть с помощью расчета облачной маски. Расчет облачной маски основан на использовании спектральных каналов ближней ИК-области спектра [8] (данные измерений в третьем (459-479 нм) и шестом (1628-1652 нм) каналах прибора МОБ18 используются для расчета нормализованного разностного индекса снега
N081 = (К4 - Кб)/(К4 + Кб), (2)
где К4, К6 - измеренные значения прибором МОБ18 в соответствующих каналах).
Классификация участков изображения на типы: облачность, снег и чистая поверхность - проводится на основе пороговых условий: 1) «снег», если К4 > 0,1 и N081 > 0,4; 2) «облачность», если К4 > 0,1 и -0,2 < N081 < 0,4; 3) «полупрозрачная облачность», если К4 > 0,1 и -0,2 > N081 > -0,5; 4) «чистая поверхность» во всех остальных случаях). Затем выделяются участки изображения, покрытые облачностью, и формируется для каждого пикселя изображения признак разделения суша / облачная поверхность.
Методика выравнивания может быть представлена последовательностью действий:
1. По данным спектральных каналов 1 и 2 МОБ18 рассчитывается
2. Первый этап расчёта маски нормализованных значений N0^ состоит в том, что по данным спектральных каналов 4 и 6 МОБ18 определяется влияние облачности (пикселю присваивается значение 1, где нет облаков и дымки, в противном случае - 0).
3. Второй этап учета расчёта маски нормализованных значений К0У1 основан на том факте, что для бесснежного периода года К0У1 плавно растёт, а затем плавно снижается, следовательно, значения, которые намного меньше, чем предыдущие и последующие, неверны и должны быть исключены. Поэтому для всех элементов, у которых
маска равна 1, находятся максимальные значения К0У1, у которых маска также равна 1, справа и слева от текущего элемента на ширину окна (оптимальная ширина около 5). Если оба максимума больше 0,1 и оба максимума больше, чем значение текущего элемента, умноженное на коэффициент (из практических соображений взят равным 1,5), то маска текущего элемента становится равной 0.
4. Третий этап также основан на том факте, что К0У1 должен изменяться плавно, и состоит в удалении слишком сильных выбросов значений К0У1 как в большую, так и в меньшую сторону, от соседних. Если значение К0У1 с маской, равной 1, слишком сильно отличается от таких же стоящих перед ним и после него, а сами эти значения отличаются друг от друга слабо, то маска текущего элемента становится равной 0.
5. После этого для элементов с маской, равной 0, рассчитывается значение из тех, у которых маска равна 1, с помощью линейной интерполяции.
6. Проводится сглаживание методом скользящего среднего (окно 11 точек).
Далее, сопоставив значения ряда за один
год (для примера) с параметрами облачной маски согласно предлагаемой методике, получим кривую, которая является уже существенно более гладкой (рис. 2, а - сплошная кривая, по сравнению с исходным рядом рис. 2, б).
Однако точность определения параметров облачной маски не идеальна и возникают выбросы, но уже другой величины, которые легко исправить обычным сглаживанием или подгонкой некоторой известной функции [9-10].
Сравнение результатов, представленных на рис. 2, а и б для сглаженных значений (см. на рисунках пунктирные линии) и подгонки известного вида функции (для примера выбрана функция Гаусса -см. на рисунках точки), показывает, что сглаженные значения временного ряда N0^ существенно отличаются друг от друга по форме и амплитуде, что однозначно при определении тренда и параметров ряда приведет к значительным погрешностям.
Результаты
На рис. 3 представлены результаты обработки данных вегетационного индекса К0У1 для разнообразных типов поверхности (река, сосновый лес и город), а также их сглаженных значений. Из рис. 3 видно, что использование ежедневных значений измерений спутниковым прибором МОБ18 позволяет получать физически адекватную модель фенологического поведения вегетационного индекса [10] согласно спектральным коэффициентам отражения различных типов поверхности. Хорошо разделимы по форме и амплитуде выбранные для анализа (в данной статье) типы поверхности. Таким образом, разработанная нами методика является средством оценки временного ряда вегетационного индекса (с существенно отличающимися коэффициентами отражения) и значит может быть использована для решения практических задач.
Вегетационный индекс NDVI 0,5 0,4 0:3 0,2 0,1 0.0 -0:1
Номер дня
п—1—i—'—i—1—i—■—i—1—i—1—i—■—i—1—i—■—i—1—i—1—г 30 60 90 120 150 180 210 240 270 300 330 ЗбС а
Вегетационный индекс NDVI
0,6 т
0,5 -0,40,3 -
—|—1—|—I—|—I—|—I—]—I—|—1—|—I—|—I—г
120 150 ISO 210 240 270 .100 330 360 б Номер дня
Рис. 2. Выравнивание временного хода значений вегетационного индекса NDVI в районе г. Томска методом скользящего среднего (пунктир) и вписывания функции Гаусса (точки) для а - учета облачной маски и б - исходного ряда
Вегетационный индекс NDVI
•0,1
120 150 180 210 240 270 ЗЙО 330 3 Номер дня
Рис. 3. Значения временного ряда вегетационного индекса NDVI для разнообразных типов поверхности (река, сосновый лес и город), а также их сглаженных значений
Отметим, что данная методика пригодна для массовых вычислений за счет высокой скорости и точности учета разнообразных изменений значений вегетационного индекса. Учитывая, что значения вегетационного индекса зависят от температуры и влажности, данная методика подлежит уточнению, что будет сделано в следующих работах.
Заключение
В статье приведена методика выравнивания временного ряда вегетационного индекса NDVI, которая основана на нескольких шагах, связанных с расчетом непосредственно вегетационного индекса, облачной маски, двухэтапного учета облачной маски и сглаживания методом скользящего среднего. Предлагаемая методика апробирована для реальных спутниковых данных спектрорадиометра MODIS для территории, расположенной в окрестностях г. Томска. На данной территории были выбраны три пикселя с однородной структурой типов поверхности: река, город и лес. Методика показала высокую скорость и точность восстановления формы временного ряда вегетационного индекса и устойчивость к значительным отклонениям вегетационного индекса (выбросам).
Литература
1. Белов В.В. От физических основ, теории и моделирования к тематической обработке спутниковых изображений. - Томск: Изд-во Ин-та оптики атмосферы СО РАН, 2005. - 265 с.
2. Сухих В.И. Аэрокосмические методы в лесном хозяйстве и ландшафтном строительстве. - Йошкар-Ола, 2005. - 390 с.
3. Шовенгердт Р.А. Дистанционное зондирование. Модели и методы обработки изображений. - М.: Техносфера, 2010. - 582 с.
4. Катаев М.Ю. Обнаружение экологических изменений природной среды по данным спутниковых измерений / М.Ю. Катаев, А.А. Бекеров // Оптика атмосферы и океана. -2014. - Т. 27, № 7. - С. 652-656.
5. Катаев М.Ю. Интернет-информационная система накопления, обработки и анализа спутниковых данных MODIS / М.Ю. Катаев, А.А. Бекеров, А.К. Лукьянов // Доклады ТУСУРа. - 2015. - № 1(35). - С. 93-99.
6. Черепанов А.С. Вегетационные индексы / А.С. Черепанов // Геоматика. - 2011. - № 2. - С. 98-102.
7. 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. -№ 113. - P. 248-258.
8. Hall D.K. Development of methods for mapping global snow cover using moderate resolution imaging spectroradi-ometer data / D.K. Hall, G.A. Riggs, V.V. Salomonson // Remote Sensing of Environment. - 1995. - № 54. -P. 127-140.
9. Jonsson P. TIMESAT - a program for analyzing time-series of satellite sensor data / P. Jonsson, L. Eklundh // Computers and Geosciences. - 2004. - № 30. - P. 833-845.
10. Zhang X. Monitoring vegetation phenology using MODIS / X. Zhang, M. Friedl, C. Schaaf // Remote Sensing of Environment. - 2003. - № 84. - P. 471-475.
Катаев Михаил Юрьевич
Д-р техн. наук, профессор каф. автоматизированных
систем управления (АСУ) ТУСУРа,
профессор Юргинского технологического института
(филиала) Национального исследовательского
Томского политехнического университета
Тел.: (382-2) 70-15-36, 8-960-975-27-85
Эл. почта: [email protected]
Бекеров Артур Александрович
Аспирант Института мониторинга климато-экологических
систем СО РАН, Томск
Тел.: (382-2) 70-15-36
Эл. почта: [email protected]
Лукьянов Андрей Кириллович
Ассистент каф. АСУ ТУСУРа
Тел.: (382-2) 70-15-36
Эл. почта: [email protected]
Kataev M.Yu., Bekerov A.A., Lukyanov A.K. Alignment method of the time series for NDVI vegetation index derived from spectroradiometer MODIS
The article provides a description of the methods, allowing to fit the time series of the vegetation index NDVI (Normalized Difference Vegetation Index) for a spatial resolution of 250 m and daily measurements. This problem occurs when it is necessary to analyze the dynamics of the vegetation index to solve specific problems of agriculture and forestry, and could be useful for the Ministry of emergency situations, etc. The method is based on correction values of the vegetation index, depending on the parameters cloud mask and the values of temperature. The processing of real data obtained by satellite measurements near the city of Tomsk is described. Keywords: MODIS, time series, vegetation index, methods of alignment.
