CC BY
115
УДК 504.064.37
М.Ю. Катаев, С.Г. Катаев, А.А. Бекеров
Методика поиска изменений из анализа спутниковых данных спектрорадиометра MODIS
Статья посвящена методам нахождения изменений на поверхности Земли из анализа спутниковых данных спектрорадиометра MODIS. Приведен обзор типичных подходов к решению задачи поиска изменений и предлагается новая методика. Обсуждаются результаты применения предлагаемой методики к обработке данных спутниковых измерений. Ключевые слова: дистанционные спутниковые методы, спектрорадиометр MODIS, спутниковое изображение, априорные данные, методы обработки и анализа.
В настоящее время спутниковые данные наблюдения являются важнейшим глобальным и периодическим источником информации, необходимой для изучения атмосферы и поверхности Земли. С каждым годом увеличивается число спутниковых систем, обладающих разнообразными характеристиками (пространственными, временными, спектральными), что порождает значительный рост получаемых данных и приводит к необходимости развивать соответствующие методы хранения, обработки и анализа, которые должны обеспечивать эффективную и быструю работу со спутниковыми данными. Особое внимание уделяется развитию программно-информационной составляющей, позволяющей пользователям искать, получать, обрабатывать и визуализировать спутниковую информацию.
Получаемые с помощью разнообразных спутниковых систем мониторинговые данные связаны с процессом непрерывного сбора информации об интересующих параметрах исследуемого объекта на конкретной территории для определения тенденций их изменения. Существующие спутниковые системы, такие как MODIS [http://modis.gsfc.nasa.gov] и Landsat [http://landsat.gsfc.nasa.gov], уже более десяти лет регулярно поставляют спутниковую информацию, обработка которой дает возможность оценивать параметры атмосферы (аэрозоль, облачность, осадки, влажность и др.) и поверхности Земли (вегетационные индексы, отражательную способность, температуру и др.).
Важнейшей информацией, которая позволяет исследовать происходящие на поверхности Земли процессы, являются вегетационные индексы [1—3]. Классификация величин индексов и их взаимосвязей позволяет четко выделять объекты на Земле, которые относятся к разным типам поверхности. Изучение временных и пространственных изменений типов поверхности относится к задачам, которые в мировой литературе называются «поиск изменений или change detection». Поиск изменений необходим в различных приложениях, связанных с контролем состояния лесных массивов, пахотных земель, строительных территорий, оценкой площади выгоревшего леса и др. Решение конкретной задачи поиска изменений представляет собой достаточно сложный процесс, поскольку при этом приходится оперировать большими массивами не только спутниковых, но и иных измерений, что требует своеобразной организации вычислительных процедур в информационном комплексе.
Данная статья является развитием работ [4, 5] по разработке интернет-информационной системы, предназначенной для накопления, обработки и анализа данных измерений, полученных спек-трорадиометром MODIS.
Обзор подходов поиска изменений. На рис. 1 показана схема формирования отраженного и рассеянного солнечного излучения, которое является основой спутниковых сигналов в видимой и ближней ИК-областях спектра (0,3-3 мкм). Именно в этой области спектра расположены основные спектральные каналы многих спутников, предназначенных для исследования поверхности Земли. Относительно пикселя, куда направлено поле зрения спутникового прибора, спутник и Солнце располагаются под разными углами (cpsat, psun), что требует учета их величин при предварительной подготовке данных для обработки. Основные вклады в измеряемый сигнал дают отраженное от поверхности Земли солнечное излучение (I1) и излучение, которое формируется из рассеянного в воздухе солнечного излучения на частицах аэрозоля и газов, переотраженного излучения от расположенной вне пикселя площадки Т, а также рассеянного в воздухе излучения (I2). Таким образом,
поступающее на спутниковый прибор излучение является суммой и несет в себе информацию, как это видно из рис. 1, о параметрах атмосферы и поверхности Земли:
I(X) = Ix(X) + I2(X) = 4m(X) cos(9sun) р(Х) Taat(X) Tsun(^) + /2(Г), (1)
где X - спектральный канал спутникового прибора; Isun(X) - солнечное излучение на границе атмосферы; р(Х) - коэффициент спектрального отражения поверхности; /..: (л). / (л) — спектральное пропускание атмосферы на пути от Солнца до поверхности Земли и до спутника.
Из уравнения (1) видно, что основным элементом, который связан со свойствами поверхности является коэффициент спектрального отражения р(Х). Величина этого коэффициента зависит от типа поверхности и ее состояния (мокрый-сухой, вода-снег-лед и др.). Заметим, что исследование свойств поверхности требует предварительной обработки поступающей спутниковой информации и эта процедура называется атмосферной коррекцией. Без проведения такой операции получаемые характеристики поверхности земли будут зависеть от условий прохождения солнечного излучения в атмосфере.
На рис. 2 приведены основные подходы к обнаружению изменений на поверхности Земли. Рассмотрим вкратце основное их содержание [6-13].
Т Пиксел Т
Рис. 1. Схема формирования спутникового сигнала отраженного солнечного излучения
(--N
Пространство
признаков ч_
(--(--^ (--(--^ (--^
Для одного Линейные Классификация Для Временные
пикселя преобразования одного пикселя ряды --/ ч_/ ч_/ ч_/ ч___/
Рис. 2. Обзор методов основных подходов к обнаружению изменений на поверхности Земли
Двухвременные изображения позволяют построить алгоритмы обнаружения изменений, основанные на сравнении двух изображений Д^) и Д?2) с разными датами (^ и ?2). На выходе алгоритма получается, как правило, маска изменений С(г, у) = /(?1)//(?2) (или С(г, у) = Д^) — Д?2)), которая представляет собой матрицу, состоящую из пикселей, в которой С(г, у) = 1, если состояние поверхности Земли в пикселе изменилось, и С(г, у) = 0 в противном случае. Развитие этого подхода связано с тем, что вводят некоторые уровни изменений и тогда при определении маски С, происходит сравнение получаемых величин с порогом. Особенностью этих методов является то, что изменения имеют знак измерений (положительный или отрицательный), и это привносит сложность при разработке алгоритмов. Подход, основанный на учете пороговых соотношений, может быть реализован с помощью иного варианта, когда на основе изменений строят распределение значений разности двух изображений, представленное на рис. 3. В этом случае порогом будут являться границы крыльев полученной гистограммы, где сконцентрированы максимальные положительные и отрицательные изменения, а в центре расположены типичные значения, в которых изменения минимальны. Так как вариации состояния атмосферы от измерения к измерению могут быть значительными, для выравнивания величин измеряемых сигналов, значения одного из сигналов нормируют на среднее ц и дисперсию ст другого сигнала. Тогда пересчитанное значение (1п), например, второго сигнала по уровню первого будет таким: /и(?2) = ст1(/(?2) — ц2)/ст2 + ц1. В этом случае элементы матрицы С рассчитываются так: С(г, у) = Д^) - /и(?2). Эти методы весьма эффективны для исследования ло-
кальных территории, но их основной недостаток заключается в том, что они крайне чувствительны к шуму измерений и неточности атмосферной коррекции.
Рис. 3. Распределение значений разности двух изображений
Методы обнаружения на основе линейных преобразований являются гораздо более надежными, чем пиксель-ориентированные подходы, однако их эффективность достигается для осредненных пространственных областей. Именно поэтому они менее восприимчивы к шуму и вариациям пороговых значений. Эти методы основаны на факте, что за относительно небольшой промежуток времени значения измеряемых сигналов, формирующих изображения, меняются незначительно.
Это позволяет записать линейную зависимость между измеряемыми изображениями в виде Д^) = А — В-/и(?2), где А, В - искомые коэффициенты. Тогда поиск изменений связан с фиксацией отклонения коэффициентов А, В от некоторых средних значений, полученных за определенный промежуток времени. Методы поиска изменений, основанные на подходах классификации, в значительной мере относятся к классу вычислительно трудоемких, так как приходится выполнять массовые матричные вычисления, например, как в подходе главных компонент при расчете собственных векторов и значений или кластеризации. К этой категории можно отнести методы, основанные на подходах искусственного интеллекта: деревья решений, нейронные сети, опорные векторы и др. Сложность в этих подходах заключается в выборе центров классов и их количества, формировании соответствующей по точности обучающей выборке.
При наличии спутниковых сигналов для определенной территории в течение длительного времени (многовременные изображения) можно строить алгоритмы поиска изменений характеристик типов земной поверхности с учетом этой «исторической» информации.
Нами рассмотрены многие из указанных подходов и разработан собственный метод, который представлен ниже по тексту.
Описание метода. Основа метода оценки изменений состояния исследуемой территории состоит в использовании классического подхода разложения временного ряда КБУ1 (рис. 4) на составляющие: трендовую, сезонную и случайную (аномалии) [14]. На рис. 4 приведены значения вегетационного индекса для одного из пикселей измерений, выполненных спутниковым прибором МОБК, в районе г. Томска. Зимой значения КБУ1 минимальны, а летом они максимальны. Вариации значений вегетационного индекса связаны с естественными причинами, такими как дождь, ветер, прозрачность атмосферы и др.
Принцип разложения временных рядов вегетационного индекса на составляющие основывается на предпосылке, что значения КБУ1, порожденные процессами меньшего масштаба, колеблются около величин, сформированных процессами большего масштаба. Так, средние месячные значения КБУТ в годовом цикле совершают колебания около долговременного тренда, аномалии средних ме-
—.—|—.—|—,—,—|—.—|—,—|—,—|—.—|—.—|-
50 100 150 200 250 300 350 400 450 500 550 600 650 700
Время, день
Рис. 4. Значения вегетационного индекса NDVI за два года (2005 и 2006 гг.)
сячных значений колеблются около годового хода. Таким образом, средние месячные значения вегетационного индекса NDVI в конкретный месяц t представляются как
N(t) = Ntr(t) + Ns(t) + AN(t), (2)
где Ntr(t) + Ns(t) + AN(t) - величина долговременного тренда NDVI, среднее значение NDVI конкретного месяца т и отклонение средних месячных значений конкретных месяцев (аномалия); t — условный порядковый номер месяца (t = 1, ..., n), начиная с первого элемента временного ряда.
Разложение ряда на компоненты осуществляется поэтапно. Вначале выделяется долговременный линейный тренд: Nr (t) = at + ß, где a — средняя скорость изменения индекса; ß - константа. Для оценивания параметров линейного тренда используется метод наименьших квадратов. Далее из исходного ряда удаляется величина тренда и для остатка определяется годовой ход. Годовой ход можно определять разными способами. Одним из них является стандартный расчет средних многолетних значений конкретного месяца. Другой способ определения годового хода заключается в описании его гармонической функцией (первый член ряда Фурье):
Ns (t)= a • sin ^(t - 1)j+b • cost -1)j , (3)
где постоянные для временного ряда коэффициенты a и b рассчитываются как
a = -2 £ N (t) • sin (|(t-1)j; b = 2 £n (t) • cos^(t-1)).
Для определения амплитуды годового хода используется выражение
A = y¡a2 + b2 . (4)
Используя амплитуду годового хода А и величину фазового сдвига ф
Ф = arctg(a/b),
NDVI представим в виде
Ns (t) = Acosjj6(t -1) -ф^ .
Аномалии (остаток) температуры получаются как разность между фактическим значением температуры конкретного месяца и суммой долговременного тренда и сезонной составляющих:
AN (t) = N (t) - [ Ntr (t)+Ns (t)].
Таким образом, определив для каждого пикселя измеренного изображения спутниковым прибором основные параметры детерминированных составляющих временного ряда (среднее значение, амплитуда и фаза временного ряда) и зная остаток AN(t), мы имеем возможность сравнивать эти величины для различных временных отрезков. Понятно, что чем больше отличий в этих значениях за определенный промежуток времени, тем более значимые изменения произошли на этой территории (данном пикселе).
Практическая апробация метода. Для практической апробации предлагаемого подхода, нами сформирована выборка спутниковых данных спектрорадиометра MODIS со средним пространственным разрешением за десять лет (2005-2014 гг). Данными являются ежедневные продукты MODIS: атмосферно скорректированный продукт глобального дневного поверхностного отражения (Surface Reflectance) 250 м MOD09GA/MOD09GQ для спутника Terra и MYD09GA/MYD09GQ для спутника AQUA. Измерения заимствованы с ресурсов NASA:
[https://lpdaac.usgs.gov/data_access/daac2disk] или [https://lpdaac.usgs.gov/data_access/data_pool]. Спутниковые данные записаны в формате HDF (иерархический формат данных) [http://www.hdfgroup.org], который представляет собой контейнер растровых слоев и метаинформа-ции, доступных для просмотра через программу HDFView [www.hdfgroup.org/products/java/hdfview].
На рис. 5 приведены величины отклонений среднего значения для территории одного из районов г. Томска (а - промежуток времени 2005-2014 гг. и б - 2005-2010 гг.).
На рис. 5, а, б приведена территории (40x40 пикселей) постройки одного из микрорайонов г. Томска (по осям обозначены номера пикселей, по оси абсцисс - номера пикселей по долготе, а ординат - по широте). Видно, что на период начала строительства (см. рис. 5, а), изменения природной среды были минимальными (серый цвет), и лишь только первый высотный дом и инфраструктура стройки отразилась в изменениях (темно-серый цвет). Через четыре года (см. рис. 5, б) в мик-
рорайоне было построено около 10 домов различной высотности, и изначальное состояние природной среды существенно изменилось, что четко фиксирует предлагаемая нами методика. Обнаружено, что естественные изменения укладываются в диапазон значений в абсолютном выражении [0-6], а [7-10] при появлении строительных или иных объектов (дороги, свалки и др.).
а б
Рис. 5. Величины отклонений среднего значения для территории одного из районов г. Томска: а - промежуток времени 2005-2010 гг. и б - 2005-2014 гг.)
Таким образом, разработанный нами метод может служить средством оценки изменений на территории (изменения, связанные с заменой одного типа поверхности на другой с существенно отличающимися коэффициентами отражения).
Заключение. В статье приведен метод поиска изменений, основанный на вычислении параметров, характеризующих основные составляющие временного ряда среднемесячных значений вегетационного индекса NDVI за отдельные промежутки времени. Сравнение этих параметров позволяет определить не только величину отклонения, но и направление (знак). Предлагаемый метод апробирован для спутниковых данных спектрорадиометра MODIS (для расчета NDVI использованы первый и второй спектральные каналы с пространственным разрешением 250 м на пиксель) для территории, расположенной в окрестностях г. Томска. Методика показала высокую скорость и точность определения знака и величины отклонений и оказалась устойчивой к случайным погрешностям и отклонениям вегетационного индекса.
Литература
1. Белов В.В. От физических основ, теории и моделирования к тематической обработке спутниковых изображений. - Томск : Изд-во Ин-та оптики атмосферы СО РАН, 2005. - 265 с.
2. Сухих В.И. Аэрокосмические методы в лесном хозяйстве и ландшафтном строительстве. -Йошкар-Ола. 2005. - 390 с.
3. Шовенгердт Р.А. Дистанционное зондирование. Модели и методы обработки изображений / Р. А. Шовенгердт. - М.: Техносфера, 2010. - 582 с.
4. Катаев М.Ю. Обнаружение экологических изменений природной среды по данным спутниковых измерений / М.Ю. Катаев, А.А. Бекеров // Оптика атмосферы и океана. - 2014. - Т. 27, № 7. -С. 652-656.
5. Катаев М.Ю. Интернет-информационная система накопления, обработки и анализа спутниковых данных MODIS / М.Ю. Катаев, А.А. Бекеров, А.К. Лукьянов // Доклады ТУСУРа. - 2015. -№ 1(35). - С. 93-99.
6. Gombay E. Change detection in autoregressive time series // Journal of Multivariate Analysis. -2008. - № 99(3). - P. 451-464.
7. Hawkins D.M. Fitting multiple change-point models to data // Computational Statistics & Data Analysis. - 2001. -№ 37(3). - P. 323-341.
8. Kasetkasem T. An image change detection algorithm based on Markov random field models / T. Kasetkasem, P. Varshney // IEEE Transactions on Geoscience and Remote Sensing. - 2002. - № 40(8). -P. 1815-1823.
9. Liu Y. Analysis of four change detection algorithms in bi-temporal space with a case study / Y. Liu, S. Nishiyama, T. Yano // International Journal of Remote Sensing. - 2004. - № 25. - P. 2121-2139.
10. Lu D. Change detection techniques / D. Lu, P. Mausel, E. Brond'izio, E. Moran // International Journal of Remote Sensing. - 2003. - № 25(12). - P. 2365-2401.
11. Lunetta R.S. Land-cover change detection using multi-temporal MODIS NDVI data / R.S. Lunetta, J.F. Knight, J. Ediriwickrema, J.G. Lyon // Remote Sensing of Environment. - 2006. -№ 105(2). - P. 142-154.
12. Radke R. Image change detection algorithms: a systematic survey / R. Radke, S. Andra, O. Al-Kofahi, B. Roysam // IEEE Transactions on Image Processing. - 2005. - № 14(3). - P. 294-307.
13. Melgani F. Unsupervised change-detection methods for remote-sensing images / F. Melgani, G. Moser, S.B. Serpico // Optical Engineering. - 2002. - № 41(12). - P. 3288-3297
14. Кусков А.И. Проблемы исследования геофизических полей /А.И. Кусков, С.Г. Катаев // Вестник ТГПУ. - 2000. - Вып. 2 (18). - С. 21-27.
Катаев Михаил Юрьевич
Д-р техн. наук, профессор каф. автоматизированных систем управления (АСУ) ТУСУРа, профессор Юргинского технологического института (филиала) Национального исследовательского Томского политехнического университета Тел.: 8-960-975-27-85, +7 (382-2) 70-15-36 Эл. почта: kataev.m@sibmail.com
Катаев Сергей Григорьевич
Д-р техн. наук, профессор каф. общей физики Томского государственного педагогического университета
Тел.: +7 (382-2) 52-17-51
Эл. почта: sgkataev@sibmail.com
Бекеров Артур Александрович
Аспирант Института мониторинга климатоэкологических систем СО РАН, Томск Тел.: +7 (382-2) 70-15-36 Эл. почта: artur@bekerov.ru
Kataev M.Yu., Kataev S.G., Bekerov A.A.
Technique for change detection using the analysis of satellite data radiometer MODIS
The article presents method of the change detection on the Earth's surface using of the analysis of satellite data radiometer MODIS. The authors provide an overview of typical approaches to the problem of change detection and propose a new approach. The results of applying the proposed method to process data of satellite measurements are described.
Keywords: remote satellite methods, MODIS, satellite image, apriori data, methods of processing and analysis.
