Разработка и верификация технологии прогнозирования коэффициентов спектральной яркости подстилающей поверхности по данным MODIS Terra/Aqua

УДК 528.88

Разработка и верификация технологии прогнозирования коэффициентов спектральной яркости подстилающей поверхности по данным MODIS — Terra/Aqua *

Д.А. Прокопов1, А.П. Жуков2

1 Управление по обеспечению мероприятий в области гражданской обороны, чрезвычайных ситуаций и пожарной безопасности

в Алтайском крае

2 Алтайский государственный университет

Development and Verification of the Technology for Land Surface Reflectances Prediction from MODIS — Terra/Aqua Data

D.A. Prokopov1, A.P. Zhukov2

1 Administration for Provision of Civil Defense, Emergency Situation, and Fire Safety Measures in Altai Region

2 Altai State University

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

Целью данной работы является разработка и верификация технологии построения ожидаемых (прогностических) КСЯ подстилающей поверхности по данным спектрорадиометра MODIS, установленного на платформах Terra и Aqua. Обсуждается алгоритм расчета прогностических КСЯ, основанный на использовании модели двунаправленного коэффициента отражения оперативного алгоритма MODIS. Представлены результаты сравнения прогностических и наблюдаемых данных по КСЯ в первых семи каналах MODIS. На основании анализа полученных данных делается вывод о пригодности прогностических КСЯ каналов 1, 3, 4 MODIS для поиска пикселей, в которых произошли изменения свойств подстилающей поверхности.

Ключевые слова: Западная Сибирь, изменение свойств подстилающей поверхности, коэффициент спектральной яркости, двунаправленный коэффициент отражения, MODIS, Terra, Aqua.

DOI 10.14258/izvasu(2014)1.2-37

Changes in properties of land surface are caused by manmade and natural phenomena including emergencies. Fast detection of the changes is possible with the near-real time processing of satellite data only. Identification of the territory, where the changes occurred, is performed by comparison of expected and observed surface reflectances. Obtaining the reliable estimates of the expected reflectances in accordance with previous satellite observations of the region land surface for future identification is a key element for the task.

The goal of this research is development and verification of the technology for expected (forecast) reflectances of the land surface from the data of MODIS spectroradiometer, installed on Terra and Aqua platforms. An algorithm for predicted reflectancies calculation based on bidirectional reflectance distribution function model from MODIS operational algorithm is discussed. A comparison of predicted and observed reflectances in first seven MODIS channels is presented. Our analysis shows that predicted reflectances from 1, 3, 4 MODIS channels are useful for the detection of pixels, where land surface properties were changed.

Key words: Western Siberia, changing the properties of the land surface, reflectance, bidirectional reflectance distribution function, MODIS, Terra, Aqua.

* Работа выполнена при частичной финансовой поддержке гранта РФФИ № 12-07-00545 и ИП 69 СО РАН.

Введение. Установление характера изменений подстилающей поверхности (ПП) Земли, обусловленных как процессами глобального потепления, естественными вариациями климата и экстремальными природными явлениями, так и антропогенным воздействием человека является одной из актуальных задач наук о Земле. Необходимые для решения этой задачи в масштабах региона объективные данные по коэффициентам спектральной яркости [1-3] подстилающей поверхности (ПП) с требуемым временным разрешением могут быть получены лишь с использованием результатов спутниковых наблюдений. Обнаружение изменений свойств подстилающей поверхности, обусловленных техногенными и природными явлениями, носящими характер чрезвычайной ситуации, при оперативном мониторинге ПП возможно лишь при использовании данных, получаемых наземной станцией региона с космической платформы в режиме реального времени.

Установление зон территории, в которых произошли изменения, осуществляется путем сопоставления ожидаемых на момент съемки коэффициентов спектральной яркости (КСЯ) поверхности с наблюдаемыми [4]. Построение достоверных оценок ожидаемых КСЯ подстилающей поверхности для будущих известных условий съемки по данным предыдущих мониторинговых спутниковых наблюдений территории региона является ключевым элементом для данной задачи.

Целью данной работы является разработка и верификация технологии построения ожидаемых (прогностических) КСЯ подстилающей поверхности по данным спектрорадиометра MODIS (MODerate resolution Imaging Spectroradiom-eter) [5], установленного на платформах Terra [6] и Aqua [7]. Обсуждается алгоритм расчета прогностических КСЯ, основанный на использовании модели двунаправленного коэффициента отражения [1,2,8] оперативного алгоритма MODIS [9]. Представлены результаты сравнения прогностических и наблюдаемых данных по КСЯ в первых семи каналах MODIS. На основании анализа полученных в работе данных делается вывод о пригодности прогностических КСЯ каналов 1, 3, 4 MODIS для поиска пикселей, в которых произошли изменения свойств подстилающей поверхности.

1. Восстановление параметров модели двунаправленного коэффициента отражения по данным MODIS. Период повторения орбит спутников Terra [6] и Aqua [7] составляет 16 дней. В каждый 16-дневный цикл в некоторой зоне подстилающей поверхности региона наблюдается MODIS при различных углах освещения (в, y>j) и наблюдения (ß, ) (рис. 1). В силу этого откорректированные на атмосферные эффекты данные MODIS [10] в каждом цикле наблюдений

могут быть использованы для восстановления параметров модели двунаправленного коэффициента отражения (ДКО).

Рис. 1. Угловые переменные освещения и наблюдения: в — зенитный угол Солнца; ё — зенитный угол спутникового прибора; рг — азимутальный угол Солнца; срг — азимутальный угол прибора

Процедура восстановления ДКО базируется на утверждении, что линейная комбинация индикатрисы оптико-геометрической модели Квео и расчетные данные в рамках теории переноса Кр дополняющие друг друга при описании прохождения квантов солнечного излучения через однородные и неоднородные растительные покровы, позволяют удовлетворительно моделировать ДКО подстилающей поверхности (см. обсуждение этого вопроса в [3,11]).

В данной работе использовалась модель ДКО операционного алгоритма МОБК [9]. В этой модели ДКО подстилающей поверхности для канала Л МОБК представляется в виде

Д8 (Л,М,<) = /1(Л)К;8о +

+ /2 (Л)^р (М, <) +

+ /з

где

1,

tp

(тг/2 - g) cos g + sing cos в + sin ß

Kgeo = O(0, ß, Ф) - sec в - sec ß+

1

+ — (1 + cosg) sec в sec ß,

О = max[0, —{t — sin t cos t) (sec в + sec ß)], n

Рис. 2. Блок-схема программного комплекса вычисления КСЯ по модели двунаправленного коэффициента отражения

cos t = max[-1,

где

. ,oV^2 + (tg0tgflsiny)2

mm{2-----, 1) ,

v sec в + sec $ n

D = max[0, \Jtg2 в + tg2 & - 2 tg в tg & cos ip\,

cos £ = cos в cos $ + sin в sin $ sin

£ = max(-1, min(arccos £, 1)).

Коэффициенты разложения fk (Л) являются параметрами, которые необходимо определить по результатам спутниковых наблюдений. Если ps (Л, в;, $i ,ipi), l = 1,...,m — коэффициенты спектральной яркости подстилающей поверхности, полученные по результатам m спутниковых измерений в течение 16-дневного цикла, то минимизируя значение квадратичной ошибки

О _ 1 (pa(A,gt,flt,yt)-i?a(A,gt,i9t,yt))3

d ^ шм

i=i

можно определить Д(Л). В этом выражении с! = (т — 3) — число степеней свободы; шл,\ — вес, присвоенный измерению ря (Л, в\ ,$\, ^).

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

n í m

ЕЕ

i=1 \l=1

Kk,i

WA,l Ki,i

fi =

к

—, (2) ^ * fiJA Г

l= 1

ШЛА

где Kk,i = K(в;,$;,<fi),k = 1(iso), 2(tp), 3(geo). Из (2) находим

nm

fk = E E i=i \i=i

K

-i

k,l

EKi,j РГ

j=1

E;

K

-i

k,i

= шл,l Ki,i

— обратная матрица системы (2).

Последовательность обработки откорректированных на атмосферные эффекты данных МОБК при получении коэффициентов разложения и альбедо обсуждается в [12,13].

2. Прогнозирование коэффициентов спектральной яркости по данным ЫСО]^.

Вычислительные эксперименты по прогнозированию КСЯ проведены с использованием разработанного программного комплекса, блок-схема которого приведена на рисунке 2. Этот комплекс позволяет вычислять КСЯ по модели ДКО (1) на синусоидальной сетке МОБК [14, 15] с пространственным разрешением 0.5; 1 и 5.6 км.

Для верификации разработанного программного комплекса были проведены сравнения результатов его работы с данными МОБК, представленными на сайте [16]. При вычислениях КСЯ в этой тестовой серии экспериментов были использованы продукты МСБ43С1 и МСБ43С2 с данными о коэффициентах разложения Д (Л) на сетке 5.6 км. Расчеты проведены для надир-ного варианта наблюдения ($ = 0). Полученные по модели ДКО (1) коэффициенты спектральной яркости сравнивались с данными продукта МСБ43С4 [16], который содержит КСЯ именно для надирной геометрии наблюдения. Распределение разностей КСЯ, вычисленных по модели ДКО и полученных из продукта МСБ43С4 для первого канала МОБК, приведено на рисунке 3. Видно, что модельные значения КСЯ практически совпадают с измеренными: абсолютное значение различия не превышает 0.0001. Для остальных спектральных каналов получены аналогичные результаты.

m

Рис. 3. Распределение разностей КСЯ, вычисленных по модели ДКО, и полученных из продукта МСБ43С4 для первого канала МОБТЯ

е 40

Рис. 4. Распределение разностей между модельными КСЯ Ямод, вычисленными по (1) с использованием Д для 16-дневного периода (18.06.2011 — 03.07.2011), и измеренными спектрорадиометром МОБТЯ Яизм для 03.07.2011

0

0.0002

0.0001

0

0.0001

0.0002

—изм —мод

0

0.1

0.08

0.06

0.04

0.02

0

0.02

0.04

0.06

0.08

0.1

- изм - - мод

е 40

Рис. 5. Распределение разностей между модельными КСЯ Ямод, вычисленными по (1) с использованием Д для 16-дневного периода (18.06.2011 — 03.07.2011), и измеренными спектрорадиометром МОБК Яизм для 04.07.2011

0

0.1

0.08

0.06

0.04

0.02

0

0.02

0.04

0.06

0.08

0.1

-изм —мод

Таблица 1

Среднеквадратичные отклонения разностей прогностических и измеренных КСЯ

Дата восста- с 24.05.12 с 24.05.12 с 12.07.10 с 12.07.10 с 15.04.10 с 20.07.11 с 18.06.11 с 18.06.11

новления по 08.06.12 по 08.06.12 по 27.07.10 по 27.07.10 по 30.04.10 по 04.08.11 по 03.07.11 по 03.07.11

Дата наблю- 08.06.12 09.06.12 27.07.10 28.07.10 01.05.10 05.08.11 04.07.11 04.07.11

дения

Платформа/ прибор Тегга/ МСЖК Тегга/ МСЖК Тегга/ МСШК Тегга/ МСШК Тегга/ МСШК Тегга/ МСЖК Адиа/ МСЖК Тегга/ МСЖК

Тайл Ь23У03 Ь23У03 Ь23У03 Ь23У03 Ь23У03 Ь22У03 Ь23У03 Ь23у03

Канал 1 0.008 0.010 0.010 0.016 0.048 0.028 0.012 0.012

Канал 2 0.017 0.022 0.021 0.029 0.040 0.035 0.024 0.024

Канал 3 0.008 0.009 0.008 0.016 0.052 0.024 0.010 0.012

Канал 4 0.008 0.010 0.009 0.016 0.050 0.027 0.011 0.012

Канал 5 0.017 0.034 0.030 0.030 0.049 0.067 0.068 0.022

Канал 6 0.017 0.020 0.017 0.022 0.029 0.039 — 0.019

Канал 7 0.011 0.016 0.015 0.017 0.025 0.034 0.018 0.014

Таким образом, в случае отсутствия изменений свойств ПП в рассматриваемом пикселе региона отличие КСЯ, восстановленных по модели ДКО, от измеренных не должно превышать суммарную погрешность модели и MODIS. Для того чтобы провести процедуру сравнения корректно, необходимо исключить облачные, теневые и другие пиксели, которые могли бы внести погрешность в результаты сравнения. Для исключения этих факторов были использованы данные продуктов MOD09GA (для спутника Terra) и MYD09GA (для спутника Aqua), полученных в [16]. В этих продуктах содержится информация о состоянии облачного покрова, наличии теней и других объектах наблюдения, КСЯ которых мог бы негативно сказаться на точности проводимых сопоставлений.

Для проверки достоверности данных, получаемых в рамках развиваемого подхода по прогнозированию КСЯ ПП, была проведена новая серия вычислительных экспериментов. В этих экспериментах были использованы данные MCD43A1 [16] с коэффициентами fk(Л) на сетке 500 м. Для задания углов наблюдения и освещения используются M0D09GA/MYD09GA.

В первом эксперименте сравнивались модельные (прогностические) коэффициенты спектральной яркости Дмод, вычисленные по (1) с использованием коэффициентов разложения f для 16-дневного периода (18.07.2011 (00:00 UTC) -03.07.2011 (23:99 UTC)), с измеренными спек-трорадиометром MODIS спектральными яркостями Дизм в последний день этого периода (т.е. 03.07.2011) для территории Алтайского края и соседних регионов (тайл h23v03 [14,15]). При этом были исключены недостоверные пиксели, о которых упоминалось выше.

На рисунке 4 приведены распределения разностей Дизм — Дмод для 100000 пикселей для каждого канала. Видно, что для всех каналов MODIS значения Дизм — Дмод не выходят за пределы интервала [—0.06; 0.06]. Вместе с тем необходимо отметить, что более 90% разностей для каналов 1, 3, 4 находятся в интервале [—0.015; 0.005].

Во втором эксперименте проверялась возможность использования восстановленных за 16-дневный период коэффициентов спектральной яр-

кости для прогнозирования КСЯ следующего 17-го дня. В этих расчетах геометрия наблюдения и измеряемые КСЯ соответствовали следующему дню — 04.07.2011. Результаты сопоставлений приведены на рисунке 5.

Аналогичные эксперименты были проведены и для других циклов измерений MODIS в 20102012 гг. Информация о некоторых экспериментах и полученных в них среднеквадратичных отклонениях для Дизм — Дмод представлены в таблице 1.

3. Заключение. Оперативное обнаружение изменений свойств подстилающей поверхности региона возможно лишь при использовании данных спутниковых наблюдений, получаемых наземной станцией с космической платформы в режиме реального времени. Установление зон, в которых произошли изменения свойств подстилающей поверхности, осуществляется путем сопоставления ожидаемых на момент съемки коэффициентов спектральной яркости поверхности с измеренными спутниковым прибором. Построение достоверных оценок ожидаемых КСЯ подстилающей поверхности для будущих известных условий съемки по данным предыдущих мониторинговых спутниковых наблюдений территории региона является ключевым элементом для данной задачи.

В данной статье представлены результаты разработки и верификации технологии построения ожидаемых КСЯ подстилающей поверхности по данным спектрорадиометра MODIS, установленного на платформах Terra и Aqua. Обсуждаются алгоритм расчета прогностических КСЯ, основанный на использовании модели двунаправленного коэффициента отражения оперативного алгоритма MODIS. Анализируются результаты серии вычислительных экспериментов по верификации разработанного программного комплекса и проверки достоверности данных, получаемых в рамках развитого подхода прогнозирования КСЯ ПП.

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

Библиографический список

1. Nicodemus F.E., Richmond J.C., Hsia J.J. et al. Geometrical Considerations and Nomenclature for Reflectance. — Washington, DC: National Bureau of Standards, US Department of Commerce, 1977. — № 160.

2. Schaepman-Strub G., Schaepman M.E., Painter T.H. et al. Reflectance Quantities in Optical Remote Sensing — Definitions and Case Studies // Remote Sens. Environ. — 2006. — Vol. 103.

3. Лагутин А.А., Никулин Ю.А., Шмаков И.А. и др. Восстановление характеристик подстилающей поверхности Сибирского региона по данным спектрорадиометра MODIS // Вычислительные технологии. — 2006. — Т. 11.

4. Roy D.P., Jin Y., Lewis P.E., Justice C.O. Prototyping a Global Algorithm for Systematic Fire Affected Area Mapping Using MODIS Time Series Data // Remote Sens. Environ. — 2005. — Vol. 97.

5. Salomonson V.V., Barnes W.L., Maymon P.W. et al. MODIS: Advanced Facility Instrument for Studies of the Earth as a System // IEEE Trans. Geosci. Remote Sens. — 1989. — Vol. 27, №2.

6. Kaufman Y.J., Herring D.D., Ranson K.J., Collatz G.J. Earth Observing System AMI Mission to Earth // IEEE Trans. Geosci. Remote Sens. — 1998. — Vol. 36, №4.

7. Parkinson C.L. Aqua: An Earth-observing Satellite Mission to Examine Water and Other Climate Variables // IEEE Trans. Geosci. Remote Sens. — 2003. — Vol. 41, №2.

8. Тимофеев Ю.М., Васильев А.В. Основы теоретической атмосферной оптики. — Спб., 2007.

9. Schaaf C.B., Gao F., Strahler A.H. et al. First Operational BRDF, Albedo Nadir Reflectance Products from MODIS // Remote Sens. Environ. — 2002. — Vol. 83.

10. Vermote E.F., Saleous N.S., Justice C.O. Atmospheric Correction of MODIS Data in the

Visible to Middle Infrared: First Results // Remote Sens. Environ. - 2002. - Vol. 35.

11. Roujean J.-L., Leroy M., Deschamps P.Y. A Bidirectional Reflectance Model of the Earth's Surface for the Correction of Remote Sensing Data //J. Geophys. Res. - 1992. - Vol. 97.

12. Лагутин А.А., Шмаков И.А., Никулин Ю.А. и др. Альбедо подстилающей поверхности по данным спектрорадиометра MODIS/ Terra // Известия Алт. гос. ун-та. — 2008. — №1(57).

13. Шмаков И.А., Лагутин А.А., Хворо-ва Л.А. Математические технологии восстановления альбедо и NDVI по данным спектрорадиометра MODIS // Фундаментальные и прикладные исследования в математической экологии и агроэкологии: международная школа-семинар. — Барнаул, 2012.

14. National Aeronautics and Space Administration [Electronic resource]. URL: http://modis-land.gsfc.nasa.gov/MODLAND_grid. html.

15. Wolfe R.E., Nishihama M., Fleiga A.J. et al. Achieving sub-pixel geolocation accuracy in support of MODIS land science // Remote Sens. Environ. -2002. - Vol. 83.

16. NASA Land Data Products and Services [Electronic resource]. URL: https://lpdaac. usgs.gov/products/modis_products _table/.