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

Эта страница опубликована в основном списке статей сайта
по адресу http://gis-lab.info/qa/atcorr-dos.html

В статье рассматривается атмосферная коррекция спутниковых снимков по методу 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λ - усиление значения радиации относительно яркости 1% тёмного объекта;

Baisλ - смещение значения радиации относительно яркости 1% тёмного объекта.

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, необходимость их сравнения со сценами спутников Landsat 5 (датчик TM) и Landsat 7 (датчик 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 — это номер канала) — максимальное значение радиации, отражённой от поверхности земли.

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

Данную проблему пытался решить Neil Flood, 2014 [9] введением дополнительных коэффициентов, рассчитанных им эмпирическим путём, для каждого канала геоснимка. В результате, коррекция значений TOA reflectance Landsat 8 в значения для того же региона и тех же спектров Landsat 7 приобретает вид формулы (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, добавляющим в Python поддержку больших многомерных массивов и матриц, а также систему низкоуровневых математических функций для операций с ними. Программа состоит из следующих структурных элементов:

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