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

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

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

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

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

Если же главным при анализе космоснимков становится не форма, рисунок и взаимное расположение объектов, а их индивидуальные спектральные характеристики, то требуется дополнительная коррекция изображений. Кроме того, тепловой канал, содержащийся в материалах Landsat, начиная с № 3, позволяет достаточно точно определить температуру подстилающей поверхности. Для этого также потребуется произвести определённые вычисления с исходными значениями.

Пример, когда радиометрическая коррекция становится необходимой

Типичным примером такого анализа является вычисление вегетационных индексов по многозональным материалам, в первую очередь – нормализованного разностного вегетационного индекса NDVI. При интерпретации схем распределения NDVI крайне важными становятся не относительные изменения значений по одной сцене, а соотношения NDVI по одной сцене или разновременным сценам в абсолютном выражении.

Например, возьмём две сцены:

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

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

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

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

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

Дело в том, что в данном случае сцены были засняты разными спутниками, аппаратура которых была настроена по-разному. Следовательно, зарегистрированные ими изображения без дополнительной радиометрической коррекции не годятся для тонкого спектрального анализа.

Для чего может потребоваться вычисление температуры земной поверхности

Дистанционное определение температуры поверхности Земли находит применение при геотермическом картировании, которое используется для решения широкого круга фундаментальных и прикладных задач:

• анализа теплового потока Земли,
• исследования распространения вечномёрзлых пород,
• обнаружения природных пожаров,
• экологического мониторинга полигонов ТБО,
• обнаружения теплового загрязнения водоёмов сбросами промышленных вод.

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

Значения, содержащиеся в файлах каналов многозонального изображения, представляют собой безразмерную величину. Она пропорциональна интенсивности излучения, достигающего орбиты, на которой находится спутник Landsat. То, что мы получаем из файла – дискретное калиброванное значение пикселя $\ Q$, – номинировано в условных относительных единицах DN (Digital Numbers – числовые значения).

Требуется получить из $\ Q$ значения отражающей способности (альбедо) $\ \rho$ наземных объектов, видимых на космоснимке. Для материалов Landsat до 7 включительно, полученных из архива USGS, это производится следующим образом (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, взятых из архива USGS, корректирующие вычисления проводятся следующим образом:

$\rho ={\frac {2\cdot 10^{-5}~Q-0{,}1}{\sin \theta }},\qquad (3)$

$R=M_{R}~Q+A_{R}.\qquad (4)$

Калибровочные коэффициенты $\ M_{R}$ и $\ A_{R}$ можно найти в файле метаданных. Из списка параметров, приведённых на рис. 5, должны быть выбраны RADIANCE_MULT ($\ M_{R}$) и RADIANCE_ADD ($\ A_{R}$), соответствующие нужному каналу.

Рис. 5. Метаданные для вычисления интенсивности излучения по данным Landsat 8

После вычисления по представленным формулам значения альбедо $\ \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$, пришедшего на сенсор спутника и зарегистрированного соответствующим тепловым каналом. Величина $\ 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-7 для вычисления температуры, Вт/м²/мкм

Диапазон	Landsat 4	Landsat 5	Landsat 7
$\ K_{1}$	671.62	607.76	666.09
$\ K_{2}$	1284.30	1260.56	1282.71

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

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

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

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

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

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

См. также

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

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