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

Привет, Вы узнаете о том , что такое атмосферная коррекция, Разберем основные их виды и особенности использования. Еще будет много подробных примеров и описаний. Для того чтобы лучше понимать что такое атмосферная коррекция, метод dos , настоятельно рекомендую прочитать все из категории Цифровая обработка изображений.

В статье рассматривается атмосферная коррекция спутниковых снимков по методу DOS (Dark Object Subtraction). В конце статьи представлен готовый программный инструментарий для выполнения атмосферной коррекции.

Содержание 1 Подготовка геоданных 2 Конвертация значений яркости исходного GeoTIFF в показатели излучения на сенсоре 3 Атмосферная коррекция 4 Особенности работы с растрами Landsat 8 5 Программная реализация нормализации геоснимков для спутников Landsat 5,7,8 6 Список использованных источников

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

Как источник геоданных использовались данные USGS (Геологической Службы Соединенных Штатов) , доступные для свободного скачивания . Данные с этих ресурсов представлены в формате 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+ в показатели излучения на сенсоре', потому коснемся только практической реализации процесса.

Конвертирование является первым этапом нормализации сырых геоданных и представляет из себя математическую операцию перевода значения яркости пикселей геоснимка в значения радиации поступившей на датчики спутника. Об этом говорит сайт https://intellect.icu . Для такого перевода в комплекте данных 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 , Lu et al 2002 .

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

где:

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

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

Gainλ — усиление значения радиации относительно яркости исходного геоснимка;

Baisλ — смещение значения радиации относительно яркости исходного геоснимка.

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

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

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

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

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

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

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

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

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

В данной работе успешно применялся метод (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 - расстояние от солнца до земли в астрономических единицах в конкретный день съемки сцены на конкретной местности ;

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 — расстояние от Солнца до Земли в астрономических единицах в конкретный день съемки сцены на конкретной местности ;

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) или, как он указывается в других источниках, 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 для всех каналов снимка , чего не может быть в реальности. Данная группа разработчиков в своем модуле нормализации и атмосферной коррекции i.landsat.toar применяет к георастрам, снятым при помощи датчиков Oli Landsat 8, те же математические методы, что для TM Landsat 5 и ETM+ Landsat 7, а недостающий коэффициент Eo (Esun) рассчитывает по формуле (10).

(10)

где:

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

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

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

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

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

Данную проблему пытался решить Neil Flood, 2014 введением дополнительных коэффициентов, рассчитанных им эмпирическим путем, для каждого канала геоснимка. В результате, коррекция значений 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 .

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

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

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

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

Исследование, описанное в статье про атмосферная коррекция, подчеркивает ее значимость в современном мире. Надеюсь, что теперь ты понял что такое атмосферная коррекция, метод dos и для чего все это нужно, а если не понял, или есть замечания, то не стесняйся, пиши или спрашивай в комментариях, с удовольствием отвечу. Для того чтобы глубже понять настоятельно рекомендую изучить всю информацию из категории Цифровая обработка изображений