Коррекция материалов Landsat

Материал из GIS-Lab
Перейти к: навигация, поиск
Эта страница является черновиком статьи.


В статье описывается упрощённая технология коррекции многозональных снимков Landsat: радиометрическая за неоднородность режима съёмки, а также пересчёт значений теплового канала в температуру

Содержание

Зачем нужна коррекция

Материалы Landsat в том виде, как они предоставляются пользователю, уже в значительной мере скорректированы, например, за кривизну Земли, особенности рельефа снимаемой территории, масштабирование в пределах возможных значений регистрируемых величин и т.п. Однако все эти коррекции направлены главным образом на обеспечение надёжного визуального дешифрирования наземных объектов. Если же главным при анализе космоснимков становится не форма, рисунок и взаимное расположение объектов, а их индивидуальные спектральные характеристики, то требуется дополнительная коррекция изображений. Типичным примером такого анализа является вычисление вегетационных индексов по многозональным материалам, в первую очередь, – нормализованного разностного вегетационного индекса NDVI. При интерпретации схем распределения NDVI крайне важным становится не относительные изменения значений по одной сцене. Более важны соотношения NDVI по одной или разновременным сценам в абсолютном выражении. Например, возьмём две сцены:

  • LE71870132004212EDC01 от 2004-07-30, зарегистрированную Landsat 7.
  • LT51860132004213KIS00 от 2004-07-31, зарегистрированную Landsat 5.

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

Рис. 1. Сравнение двух схем NDVI, построенных по сценам, зарегистрированными разными спутниками Landsat с промежутком в 1 сутки

Сравнивая эти схемы можно решить, что северная тайга в последний день июля 2004 года испытала быстрый расцвет! Судя по значениям NDVI, которые в среднем изменились на 18%, растительность повсеместно из категории угнетённой разреженной перешла в разряд пышной сомкнутой. Этот курьёзный случай может быть объяснён только одним образом – нельзя сравнивать несравнимое. Прежде чем проводить вычисления NDVI, надо было выполнить коррекцию значений каналов каждой сцены в соответствии с калибровочными коэффициентами каждой из них. Дело в том, что в данном случае сцены были засняты разными спутниками, аппаратура которых была настроена по-разному. Следовательно, зарегистрированные ими изображения без дополнительной радиометрической коррекции не годятся для тонкого спектрального анализа.

Радиометрическая коррекция

Значения, содержащиеся в файлах каналов многозонального изображения, представляют собой безразмерную величину. Она пропорциональна интенсивности излучения, достигающего орбиты, на которой летает спутник Landsat[1]. То, что мы получаем из файла – дискретное калиброванное значение пикселя  \ Q , – номинировано в условных относительных единицах DN (Digital Numbers – числовые значения). Требуется получить из  \ Q значения отражающей способности (альбедо)  \ \rho наземных объектов, видимых на космоснимке. Для материалов Landsat до 7 включительно, полученных из архива ГС США, это производится следующим образом (1):

 \rho = \frac {\pi R ~ d^2} {E \cos \theta}, \qquad (1)

где  \ R – интенсивность излучения от объекта, достигшего орбиты Landsat,  \ d – расстояние между Землёй и Солнцем,  \ E – коэффициент светимости для каждого канала,  \ \theta – высота стояния Солнца над горизонтом в момент съёмки.

Интенсивность излучения объекта, достигающего орбиты  \ R вычисляется по формуле (2):

 R = \left ( R_{\mbox{max}} - R_{\mbox{min}} \right ) \frac {Q - 1} {254} + R_{\mbox{min}}, \qquad (2)

Калибровочные коэффициенты  \ R_{\mbox{min}} и  \ R_{\mbox{max}} должны быть взяты из текстового файла метаданных с именем «*_mtl.txt», который поставляется вместе со сценой в одном архивном файле.

Если открыть этот файл для просмотра в текстовом редакторе, то там нужно найти строки, содержащие параметры, как показано на рис. 2. Из всего набора параметров следует выбрать коэффициенты, соответствующие определённому каналу.

Рис. 2. Содержимое файла с метаданными сцены Landsat 7, необходимыми для получения интенсивности излучения на сенсоре

Значение интенсивности излучения на сенсоре  \ R имеет размерность Вт/м²/ср/мкм, т.е. мощность излучения, падающая на единицу площади земной поверхности сквозь телесный угол в 1 стерадиан, взятого относительно единицы длины волны излучения.

Расстояние между Землёй и Солнцем  \ d может быть вычислено по приближённой формуле:

 d = 1 - 0{,}01668 \cos \left ( i \frac {2 \pi} {365} \right )

где  \ i – порядковый номер в году дня получения изображения. Его можно определить по дате регистрации сцены, которая также находится в файле метаданных (рис. 3):

Рис. 3. Метаданные сцены: дата регистрации (DATE_ACQUIRED)

Если для вычислений использовать MS Excel, то в нём для получения значения  \ i удобно использовать формулу =Ref-ДАТА(ГОД(Ref);1;1), где Ref – ссылка на ячейку, где записана дата регистрации.

Величина  \ d измеряется в астрономических единицах (а.е.).

Значения светимости  \ E для каждого канала могут быть взяты из следующей таблицы (табл. 1):

Табл. 1. Коэффициенты светимости для каналов Landsat, Вт/м²/мкм

Каналы 1 2 3 4 5 7 8
Landsat 4 1957 1825 1557 1033 214.9 80.72
Landsat 5 1957 1826 1554 1036 215.0 80.67
Landsat 7 1970 1842 1547 1044 225.7 82.06 1369

И последний, необходимый для вычисления альбедо, параметр  \ \theta – высота стояния Солнца над горизонтом в момент съёмки – может быть взят из файла с метаданными сцены (рис. 4). Значение это в файле приведено в градусах. Перед подстановкой в формулу его надо перевести в радианы.

Рис. 4. Метаданные сцены: высота Солнца (SUN _ELEVATION)

Для материалов Landsat 8, взятых из архива ГС США, расчётные формулы имеют гораздо более простой вид:

 \rho = \frac {R} { \sin \theta}, \qquad (3)

где

 R = 2 \cdot 10^{-5} Q - 0{,}1. \qquad (4)

После вычисления по представленным формулам значения альбедо  \ \rho по разным каналам могут быть использованы, например, для определения величины NDVI:

 \mbox{NDVI} = \frac {\rho_{_{\mbox{NIR}}} - \rho_{_{\mbox{RED}}}} {\rho_{_{\mbox{NIR}}} + \rho_{_{\mbox{RED}}}}.

Здесь NIR – номер канала ближнего инфракрасного диапазона, RED – номер канала красного диапазона. Определить какие номера каналов соответствует эти диапазонам можно с помощью табл. 2.

Табл. 2. Соответствие между диапазонами электромагнитного излучения и номерами каналов спутников разных поколений миссии Landsat[2][3][4][5][6]

Диапазон Landsat 1 Landsat 2 Landsat 3 Landsat 4 Landsat 5 Landsat 7 Landsat 8
RED 5 5 5 3 3 3 4
NIR 6 6 6 4 4 4 5
Тепловой 8 6 6 61, 62 10, 11

Определение температуры земной поверхности

По значениям тепловых каналов (табл. 2) можно определить температуру подстилающей поверхности. Теоретически точность оценки температуры около 0,5°С, однако дымка в атмосфере занижает значения на несколько градусов.

Исходными данными для определения температуры служат значения интенсивности излучения  \ R , пришедшего на сенсор спутника и зарегистрированного соответствующим тепловым каналом, в соответствии с методикой, описанной в предыдущем разделе (формулы 2, 4).

Для Landsat 4-5 температура  \ T вычисляется по следующей формуле:

 T = \frac {K_2} { \ln{ \left( \frac {K_1} {R} + 1 \right) }} - 273{,}15 \ (^{\circ} \mbox{C}), \qquad (5)

Здесь  \ K_1 и  \ K_2 – калибровочные константы, значения которых можно взять из табл. 3.

Табл. 3. Калибровочные константы для вычисления температуры, Вт/м²/мкм

Диапазон Landsat 4 Landsat 5 Landsat 7
 \ K_1 671.62 607.76 666.09
 \ K_2 1284.30 1260.56 1282.71

Пример результатов определения значений температуры подстилающей поверхности по данным теплового канала показан на рис. 5.

Рис. 5. Тепловая аномалия города Воронеж в марте 1985 г. Значения шкалы – °С

Материалы Landsat 7 содержит два тепловых канала: 61 и 62, которые отличаются настройками усиления сигнала. Для вычисления температуры можно использовать в формуле (5) любой из них или оба – они дают почти одинаковые значения (±0,1°С).

Почти также и с помощью той же формулы (5) производится вычисление температуры по данным Landsat 8, за исключением того, что значения констант  \ K_1 и  \ K_2 надо брать отдельно для каждого теплового канала данной сцены из её файла метаданных (рис. 6).

Рис. 6. Метаданные сцены Landsat 8: калибровочные константы для вычисления температуры

Значения температуры, определённые по каналам 10 и 11 (они различаются охватываемыми интервалами длин волн теплового диапазона) Landsat 8, отличаются друг от друга на 1,5-3°С. В общем случае можно осреднить их.

Литература

  1. Landsat 7 Science Data Users Handbook [Электронный ресурс] // NASA.GOV: сервер Национального управления США по воздухоплаванию и исследованию космического пространства. URL: http://landsathandbook.gsfc.nasa.gov/pdfs/Landsat7_Handbook.pdf. (дата обращения: 03.09.2015). – P. 117-120.
  2. The Multispectral Scanner System // NASA.GOV: сервер Национального управления США по воздухоплаванию и исследованию космического пространства. URL: http://landsat.gsfc.nasa.gov/?p=3227. (дата обращения: 13.09.2015).
  3. The Thematic Mapper // NASA.GOV: сервер Национального управления США по воздухоплаванию и исследованию космического пространства. URL: http://landsat.gsfc.nasa.gov/?p=3229. (дата обращения: 13.09.2015).
  4. The Enhanced Thematic Mapper Plus // NASA.GOV: сервер Национального управления США по воздухоплаванию и исследованию космического пространства. URL: http://landsat.gsfc.nasa.gov/?p=3225. (дата обращения: 13.09.2015).
  5. Operational Land Imager // NASA.GOV: сервер Национального управления США по воздухоплаванию и исследованию космического пространства. URL: http://landsat.gsfc.nasa.gov/?p=5107. (дата обращения: 13.09.2015).
  6. Thermal Infrared Sensor // NASA.GOV: сервер Национального управления США по воздухоплаванию и исследованию космического пространства. URL: http://landsat.gsfc.nasa.gov/?p=5112. (дата обращения: 13.09.2015).

См. также

Константин Силкин 22:34, 24 сентября 2015 (MSD)

Персональные инструменты
Пространства имён

Варианты
Действия
Статьи
Спецпроекты
Инструменты