Атмосферная коррекция по методу DOS

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

Анотация

Недавно была защищена магистерская диссертация по ГИС тематике, решил выдержки из данной работы опубликовать в виде цикла статей. В первой статье рассмотрим атмосферную коррекцию по методу DOS (Dark Object Subtraction). В конце статьи представлен готовый программный инструментарий для выполения атмосферной коррекции.

Подготовка геоданных

Среди источников геоданных использовались данные USGS (Геологической Службы Соединённых Штатов) [1] доступных для свободного скачивания [2]. Данные с этих ресурсов представлены в формате GeoTIFF в виде непрерывных наборов сцен для различных районов мира. Представленные на ресурсах геоданные соответствуют уровню обработки LG1: «сырые геоданные» не конвертированные в значения излучения на сенсоре. Для самостоятельного проведения конвертации и атмосферной коррекции - доступны файлы настройки спутниковых сенсоров с набором ряда параметров: минимумы и максимумы пиксельной яркости на изображении; минимумы и максимумы излучения на датчиках; минимумы и максимумы отражения излучения от поверхности земли (Landsat 8) и многое другое. Кроме данных настройки спутниковых датчиков в архиве присутствуют несколько файлов в формате GeoTIFF, распределённых по номерам каналов, количество и состав которых для разных спутников (Landsat 5,7,8) различный (Таблица 1).

Таблица 1 — Спутники и используемые в работе каналы

Спутник и датчик	Порядковый номер канала (band)	Название спектра	Длинны волн (нм)
Landsat 5 TM	1	Blue - Синий	450–520
	2	Green - Зелёный	520–600
	3	Red - Красный	630–690
	4	NIR — Ближний ИК	760–900
	5	SWIR1 — Средний ИК 1	1550–1750
	7	SWIR2 — Средний ИК 2	2080–2350
Landsat 7 ETM+	1	Blue - Синий	450–520
	2	Green - Зелёный	520–600
	3	Red - Красный	630–690
	4	NIR — Ближний ИК	770–900
	5	SWIR1 — Средний ИК 1	1550–1750
	7	SWIR2 — Средний ИК 2	2080–2350
Lansat 8 Oli	2	Blue - Синий	450—515
	3	Green - Зелёный	525—600
	4	Red - Красный	630—680
	5	NIR — Ближний ИК	845—885
	6	SWIR1 — Средний ИК 1	1560—1660
	7	SWIR2 — Средний ИК 2	2100—2300

Как видно из таблицы, шесть наиболее используемых «спектров» находятся у разных спутников в сходном диапазоне частот. Поэтому, для удобства использования терминологии, вместо указания спутниковых каналов в дальнейшем будем использовать самоназвание спектров: BLUE - синий, GREEN - зелёный, RED - красный, NIR - ближний инфракрасный, SWIR1 - средний инфракрасный 1, SWIR2 - средний инфракрасный 2 - вне зависимости от спутника и применять к ним одинаковый набор вычислений.

Как уже упоминалось, геоданные по всем каналам предоставлены уровнем обработки LG1 - в сыром виде, это означает что предоставленные канальные GeoTIFF не более чем спозиционированные на местности яркостные «фотографии», которые, в таком виде, не могут быть использованы для ГИС исследований. Поэтому данные уровня обработки LG1 необходимо нормализовать - то есть провести конвертацию значений яркости в показатели излучения на сенсоре и атмосферную коррекцию.

Конвертация значений яркости исходного GeoTIFF в показатели излучения на сенсоре

Теоретические основы процесса конвертирования подробно описаны в статье 'Конвертация данных TM, ETM+ в показатели излучения на сенсоре', потому коснёмся только практической реализации процесса.

Конвертирование является первым этапом нормализации сырых геоданных и представляет из себя математическую операцию перевода значения яркости пикселей геоснимка в значения радиации поступившей на датчики спутника. Для такого перевода в комплекте данных Landsat присутствует файл коррекции _MTL.txt, предельные значения из которого и используются на этом этапе обработки геоснимка.

Для проведения работы использовалась стандартная, описанная в документации NASA, формула (1) перевода яркости в значение top of atmosphere radiance (TOA radiance)

(1)

где:

Lλ - спектральная радиация, пришедшая на сенсор спутника;

DNcal - значения яркости пикселя сырого геоснимка;

Qcalmin - минимальное возможное значение пикселя геоснимка;

Qcalmax - максимальное возможное значение пикселя геоснимка;

LMINλ - минимальное значение спектральной радиации для конкретного датчика спутника для конкретного снимка;

LMAXλ - максимальное значение спектральной радиации для конкретного датчика спутника для конкретного снимка.

Кроме того для вычисления TOA radiance можно использовать упрощенную формулу (2) Thome et al., 1994 [3], Lu et al 2002 [4].

Lλ = DNcal × Gainλ + Baisλ (2)

где:

Lλ - спектральная радиация, пришедшая на сенсор спутника;

DNcal - значения яркости пикселя сырого геоснимка;

Gainλ - радиометрическое усиление;

Baisλ - радиометрическое смещение.

После обсчёта каждого пикселя геоснимка с использованием данной формулы, получаем матрицу значений c плавающей точкой - нормализованных геоданных.

Атмосферная коррекция

Следующим этапом нормализации геоданных является уменьшение влияния атмосферы на снимок и перевод значений радиации, дошедшей до сенсоров спутника (TOA radiance), в значения реально отражённого от земли спектрального излучения солнечного света.

Материалы данного и последующих разделов, в некоторой степени, перекликаются со статьёй "Атмосферная коррекция данных Landsat/ETM+ (COST метод)", но в отличии от неё ориентированы на пакетную обработку ДДЗ из нескольких источников.

Влияние атмосферы на геоснимок проявляется в целом ряде факторов: угол падения и отражения солнечных лучей, прозрачность атмосферы, газовый фактор и дымка (Рисунок 1).

Рисунок 1 — Факторы, влияющие на попадание отраженной солнечной радиации на сенсоры спутника

Для дальнейших исследований необходимо провести оптическую коррекцию (нормализацию) данных геоснимка методом Dark Object Subtraction (DOS), впервые представленным Chavez (1996) [5]. Суть метода состоит в нахождении яркости однопроцентно тёмного объекта геоснимка с последующей коррекцией минимума значений каждого пикселя изображения, относительно спектральной яркости найденного объекта.

Есть два основных способа поиска 1% темного объекта (Dark Object) для метода DOS:

1. эмпирический метод - подразумевает поиск значений в ручном режиме, например, с использованием инструмента гистограмма в QGIS, где изменяя нижний порог яркости гистограммы постепенно находим примерное значение яркости искомого тёмного объекта;
2. вычислительный метод подразумевает, что суммарная яркость (от 0 до n) однопроцентно тёмного объекта будет соответствовать 0.01% от суммарной яркости всех пикселей геоснимка Sobrino et al., 2004 [6].

В данной работе успешно применяли метод (2), хорошо показавший себя при обработке большого количества геоснимков исследуемого района.

После определения яркости Dark Object (в дальнейших вычислениях будем обозначать его как DNmin), производим атмосферную коррекцию по методу DOS в несколько этапов:

1.Вычисляем значение радиации, соответствующее значению яркости 1% темного объекта (расчёт производится по аналогии с TOA radiance), по формуле (3).

(3)

где:

Lλmin - спектральная радиация для 1% тёмного объекта;

DNmin - значения яркости пикселя 1% тёмного объекта;

Qcalmin - минимальное возможное значение пикселя геоснимка;

Qcalmax - максимальное возможное значение пикселя геоснимка;

LMINλ - минимальное значение спектральной радиации для конкретного датчика спутника для конкретного снимка;

LMAXλ - максимальное значение спектральной радиации для конкретного датчика спутника для конкретного снимка.

Или с использованием упрощённой формы (4).

Lλmin = DNmin × Gainλ + Baisλ (4)

где:

Lλmin - спектральная радиация для 1% тёмного объекта;

DNmin - значения яркости пикселя 1% тёмного объекта;

Gainλ - радиометрическое усиление;

Baisλ - радиометрическое смещение.

2.Рассчитываем коэффициент влияния угла падения и отражения солнечных лучей для 1% тёмного объекта по формуле (5).

(5)

где:

L_1% - коэффициент влияния угла падения и отражения солнечных лучей для 1% тёмного объекта;

d - расстояние от солнца до земли в астрономических единицах в конкретный день съёмки сцены на конкретной местности [7];

E₀ - коэффициент солнечного внеатмосферного спектрального излучения (явно представлен как табличные данные и учитывается при калибровке датчиков Landsat 5 и 7 , для Landsat 8 дополнительно вычисляется);

θ - зенитное расстояние для солнца в радианах;

T_Z - мера прохождения излучения от солнца до земли, в методе DOS2, принимается равным cosθ.

3.Вычисляем значение атмосферной дымки (hazing) по формуле (6).

Lλhaze = Lλmin - L_1% (6)

где:

Lλhaze - значение атмосферной дымки (hazing);

L_1% - коэффициент влияния угла падения и отражения солнечных лучей для 1% тёмного объекта;

Lλmin - спектральная радиация для 1% тёмного объекта.

4.Рассчитываем атмосферно скорректированные значения отражённой солнечной радиации по формуле (7).

(7)

где:

ρλ - атмоферно скорректированные значения отражённой солнечной радиации;

Lλ - значения радиации пришедшей на сенсор спутника;

Lλhaze - значение атмосферной дымки (hazing);

d - расстояние от солнца до земли в астрономических единицах в конкретный день съёмки сцены на конкретной местности [7];

E₀ - коэффициент солнечного внеатмосферного спектрального излучения (явно представлен как табличные данные и учитывается при калибровке датчиков Landsat 5 и 7 , для Landsat 8 дополнительно вычисляется);

θ - зенитное расстояние для солнца в радианах;

T_Z - мера прохождения излучения от солнца до земли, в методе DOS2, принимается равным cosθ.

Особенности работы с растрами спутника Landsat 8

Если часть сцен сняты спутником Landsat 8 датчик Oli, необходимость их сравнения со сценами спутников lansat5 (датчик TM) и landsat8 (датчик ETM+). В данном контексте известна проблема применения стандартных методов атмосферной коррекции для нового спутника: дело в том, что калибровка датчиков Oli Landsat 8 производится без учёта значений коэффициента солнечного внеатмосферного спектрального излучения Eo (Sobrino et al., 2004) [6] или ,в других источниках, Esun. Вместо использования данного коэффициента, в файл корректировки _MTL.txt, были добавлены некоторые новые спектральные параметры: REFLECTANCE_MULT_BAND - усиление значения отражения и REFLECTANCE_ADD_BAND - смещение значения отражения для каждого из спектральных датчиков. В результате, по задумке авторов изменений, расчет TOA reflectance для Landsat 8 должен производится по формуле (8).

ρλ' = MρQcal+ Aρ (8)

где:

ρλ' - значение верхнеатмосфеного планетарного отражения радиации (TOA reflectance), без учёта коррекции по углу падения и отражения солнечных лучей;

Mρ - каналоспецифичный мультипликативный расчётный фактор для (REFLECTANCE_MULT_BAND_x, где x это номер канала) - усиление значения отражения;

Aρ - каналоспецифичный мультипликативный расчётный фактор для (REFLECTANCE_ADD_BAND_x,где x это номер канала) - смещение значения отражения;

Qcal - значения яркости пикселя сырого геоснимка (DN).

А коррекция TOA reflectance с учётом угла падения и отражения солнечных лучей вычисляется по формуле (9).

(9)

где:

ρλ - значение верхнеатмосфеного планетарного отражения радиации (TOA reflectance), c учётом коррекции по углу падения и отражения солнечных лучей;

θSE - высота солнца над горизонтом. Доступна в файле _MTL.txt в параметре (SUN_ELEVATION);

θSZ - зенитное расстояние; θSZ = 90° - θSE.

Сообщество разработчиков свободной GIS GRASS указывает на одинаковые значения REFLECTANCE_MULT_BAND и REFLECTANCE_ADD_BAND для всех каналов снимка [8], что не может быть в реальности. Данная группа разработчиков в своём модуле нормализации и атмосферной коррекции i.landsat.toar применяет к георастрам, снятым при помощи датчиков Oli Landsat 8, те же математические методы, что для TM Landsat 5 и ETM+ Landsat 7, а недостающий коэффициент Eo (Esun) рассчитывает по формуле (10).

(10)

где:

Esun - (E0) вычисленный коэффициент солнечного внеатмосферного спектрального излучения;

d - расстояние от солнца до земли в астрономических единицах в конкретный день съёмки сцены на конкретной местности [7];

RADIANCE_MAXIMUM - каналоспецифичный мультипликативный расчётный фактор для (RADIANCE_MAXIMUM_x, где x это номер канала) - максимально возможное значение поступающей на сенсор радиации;

REFLECTANCE_MAXIMUM - каналоспецифичный мультипликативный расчётный фактор для (REFLECTANCE_MAXIMUM_x, где x это номер канала) - максимальное значение радиации, отражённой от поверхности земли.

Ещё одной особенностью Landsat8 стало уменьшение чувствительности каналов RED, NIR и SWIR1 относительно Landsat 7 и 5 , что приводит к изменению значений, вычисляемых с использованием данных спектров, индексов.

Данную проблему пытался решить Neil Flood, 2014 [9] введением дополнительных, рассчитанных им эмпирическим путём, коэффициентов для каждого канала геоснимка. В результате, коррекция значений TOA reflectance Landsat8 в значения для того же региона и тех же спектров Landsat7, приобретает вид формулы (11).

ρETM+ = c0 + c1 × ρOLI (11)

где:

ρETM - результат конвертирования значения TOA reflectance из Landsat 8 в Landsat 7;

ρOLI - вычисленное для Landsat8 значение TOA reflectance;

c0 - смещение значения отражённой радиации между датчиками OLI и ETM+;

c1 - усиление значения отражённой радиации между датчиками OLI и ETM+.

Значения c0 и с1 формулы (11) были сведены автором в корректировочные таблицы, что позволяет применять их без повторения вычислений представленных в статье Neil Flood, 2014 [9].

Программная реализация нормализации геоснимков для трёх типов спутников Landsat: 5,7,8

Все вычисления производились с использованием языка программирования python. Для этого была написана собственная программа нормализации геоданных с использованием интерфейса GDAL (библиотеки для работы с растровыми географическими форматами), тесно связанного с расширением numpy, добавляющим в язык программирования поддержку больших многомерных массивов и матриц, а также систему низкоуровневых математических функций для операций с этими массивами. Программа состоит из следующих структурных элементов:

1. Класса преобразования растра в массив.
2. Класса преобразования вычисленного массива в растр.
3. Класса сбора данных коррекции при помощи парсера _MTL.txt файла и данных выгрузки расстояний от земли до солнца [7].
4. Выполняемого скрипта самих математических вычислений растров для проведения процесса нормализации.

Исходные коды программы доступны для загрузки по адресу: https://github.com/oldbay/raster_tools, кроме того исходные GEOTIFF и нормализованные спектральные растры доступны по адресу: https://github.com/oldbay/paper_examples

Список использованных источников

1. Официальный сайт Геологической Службы Соединённых Штатов - USGS [Электронный ресурс]. Режим доступа: http://www.usgs.gov.
2. Интреактивный WEB-интерфейс получения данных дистанционного зондирования с ресурса Геологической Службы Соединённых Штатов [Электронный ресурс]. Режим доступа: http://glovis.usgs.gov.
3. Thome K.J., Biggar S.F., Gellman D.L., Slater P.N. (1994) - Absolute-radiometric calibration of Landsat-5 Thematic Mapper and the proposed calibration of the Advanced Spaceborne Thermal Emission and Reflection Radiometer. Paper presented at the Geoscience and Remote Sensing Symposium, 1994. IGARSS’94. Surface and Atmospheric Remote Sensing: Technologies, Data Analysis and Interpretation.
4. Lu D., Mausel P., Brondizio E., Moran E. (2002) - Assessment of atmospheric correction methods for Landsat TM data applicable to Amazon basin LBA research. International Journal of Remote Sensing, 23 (13): 2651-2671.
5. Chavez, P. S. (1996). Image-based atmospheric correction—revisited and improved. Photogrammetric Engineering and Remote Sensing, 62(9), 1025 – 1036.
6. Sobrino J. A., Jiménez-Muñoz J.C., Paolini L. (2004) - Land surface temperature retrieval from LANDSAT TM 5. Remote sensing of Environment, 90 (4): 434-440.
7. «HORIZONS Web-Interface» [Electronic resource]. Access mode: http://ssd.jpl.nasa.gov/horizons.cgi.
8. «Calculates top-of-atmosphere radiance or reflectance and temperature for Landsat MSS/TM/ETM+/OLI» [Electronic resource]. Access mode: http://grass.osgeo.org/grass65/manuals/i.landsat.toar.html.
9. Neil Flood (2014) Continuity of Reflectance Data between Landsat-7 ETM+ and Landsat-8 OLI, for Both Top-of-Atmosphere and Surface Reflectance: A Study in the Australian Landscape. Remote Sens. 2014, 6, 7952-7970