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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Табл. 3.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$ – высота стояния Солнца над горизонтом в момент съёмки – может быть взят из файла с метаданными сцены (рис. 3.3). Значение это в файле приведено в градусах. Перед подстановкой в формулу его надо перевести в радианы.

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

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

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

где

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

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

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

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

Табл. 3.2. Соответствие между диапазонами электромагнитного излучения и номерами каналов спутников разных поколений миссии Landsat [12-16]

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

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

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

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

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

$T={\frac {K_{2}}{\ln {\left({\frac {K_{1}}{R}}+1\right)}}}-273{,}15\ (^{\circ }{\mbox{C}}),\qquad (3.5)$

Здесь $\ K_{1}$ и $\ K_{2}$ – калибровочные константы, значения которых можно взять из табл. 3.3.

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

Диапазон Landsat 4 Landsat 5 Landsat 7 K1 671.62 607.76 666.09 K2 1284.30 1260.56 1282.71

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

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

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

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

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

Значения температуры, определённые по каналам 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. 39. 2. Landsat 7 Science Data Users Handbook [Электронный ресурс] // NASA.GOV: сервер Национального управления США по воздухоплава-нию и исследованию космического пространства. URL: http://landsathandbook.gsfc.nasa.gov/pdfs/Landsat7_Handbook.pdf. (дата обращения: 03.09.2015). – P. 38. 3. EarthExplorer Help Documentation. [Электронный ресурс] // USGS.GOV: сервер Геологической службы США. URL: http://earthexplorer.usgs.gov/documents/helptutorial.pdf. (дата обращения: 05.08.2015). – 2013. – 73 pp. 4. An introduction and reference for MultiSpec. [Электронный ресурс] // PURDUE.EDU: сервер университета Пердью. URL: https://engineering.purdue.edu/~biehl/MultiSpec/MultiSpec_Intro_9_11.pdf. – 2011. – 189 pp. (дата обращения: 31.08.2015). 5. Tutorial: Combining separate image files into a single multispectral image file. [Электронный ресурс] // PURDUE.EDU: сервер университета Пердью. URL:https://engineering.purdue.edu/~biehl/MultiSpec/tutorials/MultiSpec_Tutorial_5.pdf. – 2013. – 10 pp. (дата обращения: 31.08.2015). 6. MultiSpec exercise: Image enhancement. [Электронный ресурс] // PURDUE.EDU: сервер университета Пердью. URL: https://engineering.purdue.edu/~biehl/MultiSpec/tutorials/MultiSpec_Exercise_2.pdf. – 2009. – 5 pp. (дата обращения: 31.08.2015). 7. MultiSpec exercise: Selecting areas and the coordinate view. [Электронный ресурс] // PURDUE.EDU: сервер университета Пердью. URL: https://engineering.purdue.edu/~biehl/MultiSpec/tutorials/MultiSpec_Exercise_7.pdf. – 2009. – 2 pp. (дата обращения: 31.08.2015). 8. MultiSpec exercise: Creating vegetation indices images. [Электронный ре-сурс] // PURDUE.EDU: сервер университета Пердью. URL: https://engineering.purdue.edu/~biehl/MultiSpec/tutorials/MultiSpec_Exercise_8.pdf. – 2010. – 7 pp. (дата обращения: 31.08.2015). 9. MultiSpec exercise: Overlay shape files on image window. [Электронный ре-сурс] // PURDUE.EDU: сервер университета Пердью. URL: https://engineering.purdue.edu/~biehl/MultiSpec/tutorials/MultiSpec_Exercise_6.pdf. – 2009. – 2 pp. (дата обращения: 31.08.2015). 10. MultiSpec exercise: Supervised classification. [Электронный ресурс] // PURDUE.EDU: сервер университета Пердью. URL: https://engineering.purdue.edu/~biehl/MultiSpec/tutorials/MultiSpec_Exercise_4.pdf. – 2009. – 7 pp. (дата обращения: 31.08.2015). 11. Landsat 7 Science Data Users Handbook [Электронный ресурс] // NASA.GOV: сервер Национального управления США по воздухоплава-нию и исследованию космического пространства. URL: http://landsathandbook.gsfc.nasa.gov/pdfs/Landsat7_Handbook.pdf. (дата обращения: 03.09.2015). – P. 117-120. 12. The Multispectral Scanner System // NASA.GOV: сервер Национального управления США по воздухоплаванию и исследованию космического пространства. URL: http://landsat.gsfc.nasa.gov/?p=3227. (дата обращения: 13.09.2015). 13. The Thematic Mapper // NASA.GOV: сервер Национального управления США по воздухоплаванию и исследованию космического пространства. URL: http://landsat.gsfc.nasa.gov/?p=3229. (дата обращения: 13.09.2015). 14. The Enhanced Thematic Mapper Plus // NASA.GOV: сервер Национального управления США по воздухоплаванию и исследованию космического пространства. URL: http://landsat.gsfc.nasa.gov/?p=3225. (дата обращения: 13.09.2015). 15. Operational Land Imager // NASA.GOV: сервер Национального управления США по воздухоплаванию и исследованию космического пространства. URL: http://landsat.gsfc.nasa.gov/?p=5107. (дата обращения: 13.09.2015). 16. Thermal Infrared Sensor // NASA.GOV: сервер Национального управления США по воздухоплаванию и исследованию космического пространства. URL: http://landsat.gsfc.nasa.gov/?p=5112. (дата обращения: 13.09.2015).