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

Эта страница является черновиком статьи.

В статье описывается упрощённая технология коррекции многозональных снимков 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. То, что мы получаем из файла – дискретное калиброванное значение пикселя $\ Q$, – номинировано в условных единицах DN (Digital Numbers – числовые значения). Требуется получить из $\ Q$ значения отражающей способности (альбедо) $\ \rho$ наземных объектов, видимых на космоснимке. Для материалов Landsat до 7 включительно, полученных из ГС США, это производится следующим образом (1):

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

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

Интенсивность излучения объекта, достигающего орбиты 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

Диапазон	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°С.

См. также

• Технические характеристики инструментов ДЗЗ и их носителей
• Получение бесплатных космических снимков Landsat TM,ETM+ через Glovis
• Работа с архивом материалов ДЗЗ через EarthExplorer
• Конвертация данных TM, ETM+ в показатели излучения на сенсоре
• Конвертация данных Landsat TM/ETM+ в значения температуры - Теория
• Определение расстояния от Земли до Солнца в момент получения снимка
• NDVI - теория и практика