Атмосферная коррекция по методу DOS: различия между версиями

Материал из GIS-Lab
Перейти к навигации Перейти к поиску
мНет описания правки
 
(не показано 14 промежуточных версий 2 участников)
Строка 1: Строка 1:
{{Статья|Черновик}}
{{Статья|Опубликована|atcorr-dos}}
==<b>Анотация</b>==
{{Аннотация|В статье рассматривается атмосферная коррекция спутниковых снимков по методу DOS (Dark Object Subtraction). В конце статьи представлен готовый программный инструментарий для выполнения атмосферной коррекции.}}


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


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


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


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


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


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


<br />


==<b>Конвертация значений яркости исходного GeoTIFF в показатели излучения на сенсоре</b>==
 
=='''Конвертация значений яркости исходного GeoTIFF в показатели излучения на сенсоре'''==


Теоретические основы процесса конвертирования подробно описаны в статье [http://gis-lab.info/qa/dn2radiance.html 'Конвертация данных TM, ETM+ в показатели излучения на сенсоре'], потому коснёмся только практической реализации процесса.
Теоретические основы процесса конвертирования подробно описаны в статье [http://gis-lab.info/qa/dn2radiance.html 'Конвертация данных TM, ETM+ в показатели излучения на сенсоре'], потому коснёмся только практической реализации процесса.
Строка 107: Строка 67:
Для проведения работы использовалась стандартная, описанная в документации NASA, формула (1) перевода яркости в значение top of atmosphere radiance (TOA radiance)  
Для проведения работы использовалась стандартная, описанная в документации NASA, формула (1) перевода яркости в значение top of atmosphere radiance (TOA radiance)  


<br />
 
::[[Файл:F_01.png]]    (1)
::[[Файл:F_01.png]]    (1)


<br />
 
где:
где:
::''Lλ'' - спектральная радиация, пришедшая на сенсор спутника;
::''Lλ'' спектральная радиация, пришедшая на сенсор спутника;
::''DNcal'' - значения яркости пикселя сырого геоснимка;
::''DNcal'' значения яркости пикселя сырого геоснимка;
::''Qcalmin'' - минимальное возможное значение пикселя геоснимка;
::''Qcalmin'' минимальное возможное значение пикселя геоснимка;
::''Qcalmax'' - максимальное возможное значение пикселя геоснимка;
::''Qcalmax'' максимальное возможное значение пикселя геоснимка;
::''LMINλ'' - минимальное значение спектральной радиации для конкретного датчика спутника для конкретного снимка;
::''LMINλ'' минимальное значение спектральной радиации для конкретного датчика спутника для конкретного снимка;
::''LMAXλ'' - максимальное значение спектральной радиации для конкретного датчика спутника для конкретного снимка.
::''LMAXλ'' максимальное значение спектральной радиации для конкретного датчика спутника для конкретного снимка.
 
 
Кроме того, для вычисления TOA radiance можно использовать упрощенную формулу (2) Thome et al., 1994 [3], Lu et al 2002 [4].


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


<br />
::''Lλ = DNcal × Gainλ + Baisλ''    (2)
::''Lλ = DNcal × Gainλ + Baisλ''    (2)


<br />
 
где:
где:
::''Lλ'' - спектральная радиация, пришедшая на сенсор спутника;
::''Lλ'' спектральная радиация, пришедшая на сенсор спутника;
::''DNcal'' - значения яркости пикселя сырого геоснимка;
::''DNcal'' значения яркости пикселя «сырого» геоснимка;
::''Gainλ'' - усиление значения радиации относительно яркости исходного геоснимка;
::''Gainλ'' усиление значения радиации относительно яркости исходного геоснимка;
::''Baisλ'' - смещение значения радиации относительно яркости исходного геоснимка.
::''Baisλ'' смещение значения радиации относительно яркости исходного геоснимка.
 
 
После обсчёта каждого пикселя геоснимка с использованием данной формулы получаем матрицу значений c плавающей точкой — нормализованные геоданные.  


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


<br />


==<b>Атмосферная коррекция</b>==
== Атмосферная коррекция ==


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


Материалы данного и последующих разделов, в некоторой степени, перекликаются со статьёй [http://gis-lab.info/qa/atcor.html "Атмосферная коррекция данных Landsat/ETM+ (COST метод)"], но в отличии от неё ориентированы на пакетную обработку ДДЗ из нескольких источников.
Материалы данного и последующих разделов в некоторой степени перекликаются со статьёй [http://gis-lab.info/qa/atcor.html "Атмосферная коррекция данных Landsat/ETM+ (COST метод)"], но в отличии от неё ориентированы на пакетную обработку ДДЗ из нескольких источников.


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


<br />
 
[[Файл:Dos_img_1.jpg]]
[[Файл:Dos_img_1.jpg]]
<br />
Рисунок 1 — Факторы, влияющие на попадание отраженной солнечной радиации на сенсоры спутника


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


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


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


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


<br />
''' 1. Вычисляем значение радиации, соответствующее значению яркости 1%-темного объекта (расчёт производится по аналогии с TOA radiance), по формуле (3). '''
 
 
::[[Файл:Dos_f_03.png]]    (3)
::[[Файл:Dos_f_03.png]]    (3)


<br />
 
где:
где:
::''Lλmin'' - спектральная радиация для 1% тёмного объекта;
::''Lλmin'' - спектральная радиация для 1% тёмного объекта;
Строка 175: Строка 136:
::''LMAXλ'' - максимальное значение спектральной радиации для конкретного датчика спутника для конкретного снимка.     
::''LMAXλ'' - максимальное значение спектральной радиации для конкретного датчика спутника для конкретного снимка.     


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


<br />
 
::''Lλmin = DNmin × Gainλ + Baisλ'' (4)
::''Lλmin = DNmin × Gainλ + Baisλ'' (4)


<br />
 
где:
где:
::''Lλmin'' - спектральная радиация для 1% тёмного объекта;
::''Lλmin'' - спектральная радиация для 1% тёмного объекта;
::''DNmin'' - значения яркости пикселя 1% тёмного объекта;
::''DNmin'' - значения яркости пикселя 1% тёмного объекта;
::''Gainλ'' - радиометрическое усиление;
::''Gainλ'' - усиление значения радиации относительно яркости 1% тёмного объекта;
::''Baisλ'' - радиометрическое смещение.
::''Baisλ'' - смещение значения радиации относительно яркости 1% тёмного объекта.
 


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


<br />
''' 2. Рассчитываем коэффициент влияния угла падения и отражения солнечных лучей для 1%-тёмного объекта по формуле (5). '''
 
::[[Файл:Dos_f_05.png]] (5)
::[[Файл:Dos_f_05.png]] (5)


<br />
 
где:
где:
::''L<sub>1%</sub>'' - коэффициент влияния угла падения и отражения солнечных лучей для 1% тёмного объекта;
::''L<sub>1%</sub>'' - коэффициент влияния угла падения и отражения солнечных лучей для 1% тёмного объекта;
Строка 202: Строка 163:
::''T<sub>Z</sub>'' - мера прохождения излучения от солнца до земли, в методе DOS2, принимается равным cosθ.
::''T<sub>Z</sub>'' - мера прохождения излучения от солнца до земли, в методе DOS2, принимается равным cosθ.


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


<br />
''' 3. Вычисляем значение атмосферной дымки (hazing) по формуле (6). '''
 
 
::''Lλhaze = Lλmin - L<sub>1%</sub>'' (6)
::''Lλhaze = Lλmin - L<sub>1%</sub>'' (6)


<br />
 
где:
где:
::''Lλhaze'' - значение атмосферной дымки (hazing);
::''Lλhaze'' значение атмосферной дымки (hazing);
::''L<sub>1%</sub>'' - коэффициент влияния угла падения и отражения солнечных лучей для 1% тёмного объекта;
::''L<sub>1%</sub>'' коэффициент влияния угла падения и отражения солнечных лучей для 1%-тёмного объекта;
::''Lλmin'' - спектральная радиация для 1% тёмного объекта.
::''Lλmin'' спектральная радиация для 1%-тёмного объекта.
 
 
''' 4. Рассчитываем атмосферно скорректированные значения отражённой солнечной радиации по формуле (7) '''


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


<br />
::[[Файл:Dos_f_07.png]] (7)
::[[Файл:Dos_f_07.png]] (7)


<br />
 
где:
где:
::''ρλ'' - атмоферно скорректированные значения отражённой солнечной радиации;
::''ρλ'' атмоферно скорректированные значения отражённой солнечной радиации;
::''Lλ'' - значения радиации пришедшей на сенсор спутника;
::''Lλ'' значения радиации, пришедшей на сенсор спутника;
::''Lλhaze'' - значение атмосферной дымки (hazing);
::''Lλhaze'' значение атмосферной дымки (hazing);
::''d'' - расстояние от солнца до земли в астрономических единицах в конкретный день съёмки сцены на конкретной местности [7];
::''d'' расстояние от Солнца до Земли в астрономических единицах в конкретный день съёмки сцены на конкретной местности [7];
::''E<sub>0</sub>'' - коэффициент солнечного внеатмосферного спектрального излучения (явно представлен как табличные данные и учитывается при калибровке датчиков Landsat 5 и 7 , для Landsat 8 дополнительно вычисляется);
::''E<sub>0</sub>'' коэффициент солнечного внеатмосферного спектрального излучения (явно представлен как табличные данные и учитывается при калибровке датчиков Landsat 5 и 7 , для Landsat 8 дополнительно вычисляется);
::''θ'' - зенитное расстояние для солнца в радианах;
::''θ'' зенитное расстояние для солнца в радианах;
::''T<sub>Z</sub>'' - мера прохождения излучения от солнца до земли, в методе DOS2, принимается равным cosθ.
::''T<sub>Z</sub>'' мера прохождения излучения от солнца до земли, в методе DOS2, принимается равным cosθ.
 
== Особенности работы с растрами Landsat 8 ==


<br />
Если часть сцен сняты спутником 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).
==<b>Особенности работы с растрами спутника Landsat 8</b>==


Если часть сцен сняты спутником 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).


<br />
::''ρλ' = MρQcal+ Aρ'' (8)
::''ρλ' = MρQcal+ Aρ'' (8)


<br />
 
где:  
где:  
::''ρλ''' - значение верхнеатмосфеного планетарного отражения радиации (TOA reflectance), без учёта коррекции по углу падения и отражения солнечных лучей;
::''ρλ''' значение верхнеатмосфеного планетарного отражения радиации (TOA reflectance), без учёта коррекции по углу падения и отражения солнечных лучей;
::''Mρ'' - каналоспецифичный мультипликативный расчётный фактор для (REFLECTANCE_MULT_BAND_x, где x это номер канала) - усиление значения отражения;
::''Mρ'' каналоспецифичный мультипликативный расчётный фактор для (REFLECTANCE_MULT_BAND_x, где x это номер канала) усиление значения отражения;
::''Aρ'' - каналоспецифичный мультипликативный расчётный фактор для (REFLECTANCE_ADD_BAND_x,где x это номер канала) - смещение значения отражения;
::''Aρ'' каналоспецифичный мультипликативный расчётный фактор для (REFLECTANCE_ADD_BAND_x, где x это номер канала) смещение значения отражения;
::''Qcal'' - значения яркости пикселя сырого геоснимка (DN).
::''Qcal'' значения яркости пикселя «сырого» геоснимка (DN).
 


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


<br />
 
::[[Файл:Dos_f_09.png]] (9)                               
::[[Файл:Dos_f_09.png]] (9)                               


<br />
 
где:  
где:  
::''ρλ'' - значение верхнеатмосфеного планетарного отражения радиации (TOA reflectance), c учётом коррекции по углу падения и отражения солнечных лучей;
::''ρλ'' значение верхнеатмосфеного планетарного отражения радиации (TOA reflectance), c учётом коррекции по углу падения и отражения солнечных лучей;
::''θSE'' - высота солнца над горизонтом. Доступна в файле _MTL.txt в параметре (SUN_ELEVATION);
::''θSE'' высота солнца над горизонтом. Доступна в файле _MTL.txt в параметре (SUN_ELEVATION);
::''θSZ'' - зенитное расстояние; θSZ = 90° - θSE.
::''θ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).


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


<br />
::[[Файл:Dos_f_10.png]] (10)
::[[Файл:Dos_f_10.png]] (10)


<br />
 
где:
где:
::''Esun'' - (E0) вычисленный коэффициент солнечного внеатмосферного спектрального излучения;
::''Esun'' (E0) вычисленный коэффициент солнечного внеатмосферного спектрального излучения;
::''d'' - расстояние от солнца до земли в астрономических единицах в конкретный день съёмки сцены на конкретной местности [7];
::''d'' расстояние от солнца до земли в астрономических единицах в конкретный день съёмки сцены на конкретной местности [7];
::''RADIANCE_MAXIMUM'' - каналоспецифичный мультипликативный расчётный фактор для (RADIANCE_MAXIMUM_x, где x это номер канала) - максимально возможное значение поступающей на сенсор радиации;
::''RADIANCE_MAXIMUM'' каналоспецифичный мультипликативный расчётный фактор для (RADIANCE_MAXIMUM_x, где x это номер канала) максимально возможное значение поступающей на сенсор радиации;
::''REFLECTANCE_MAXIMUM'' - каналоспецифичный мультипликативный расчётный фактор для (REFLECTANCE_MAXIMUM_x, где x это номер канала) - максимальное значение радиации, отражённой от поверхности земли.
::''REFLECTANCE_MAXIMUM'' каналоспецифичный мультипликативный расчётный фактор для (REFLECTANCE_MAXIMUM_x, где x это номер канала) максимальное значение радиации, отражённой от поверхности земли.
 
 
Ещё одной особенностью Landsat 8 стало уменьшение чувствительности каналов RED, NIR и SWIR1 относительно Landsat 7 и 5 , что приводит к изменению значений индексов, вычисляемых с использованием данных спектров.


<br />
Данную проблему пытался решить Neil Flood, 2014 [9] введением дополнительных коэффициентов, рассчитанных им эмпирическим путём, для каждого канала геоснимка. В результате, коррекция значений TOA reflectance Landsat 8 в значения для того же региона и тех же спектров Landsat 7 приобретает вид формулы (11).
Ещё одной особенностью Landsat8 стало уменьшение чувствительности каналов RED, NIR и SWIR1 относительно Landsat 7 и 5 , что приводит к изменению значений, вычисляемых с использованием данных спектров, индексов.


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


<br />
::''ρETM+ = c0 + c1 × ρOLI'' (11)
::''ρETM+ = c0 + c1 × ρOLI'' (11)


<br />
 
где:
где:
::''ρETM'' - результат конвертирования значения TOA reflectance из Landsat 8 в Landsat 7;
::''ρETM'' результат конвертирования значения TOA reflectance из Landsat 8 в Landsat 7;
::''ρOLI'' - вычисленное для Landsat8 значение TOA reflectance;
::''ρOLI'' вычисленное для Landsat8 значение TOA reflectance;
::''c0'' - смещение значения отражённой радиации между датчиками OLI и ETM+;
::''c0'' смещение значения отражённой радиации между датчиками OLI и ETM+;
::''c1'' - усиление значения отражённой радиации между датчиками OLI и ETM+.
::''c1'' усиление значения отражённой радиации между датчиками OLI и ETM+.
 


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


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


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


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


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


<br />
== Список использованных источников ==
==<b>Список использованных источников</b>==


# Официальный сайт Геологической Службы Соединённых Штатов - USGS [Электронный ресурс]. Режим доступа: http://www.usgs.gov.
# Официальный сайт Геологической Службы Соединённых Штатов - USGS [Электронный ресурс]. Режим доступа: http://www.usgs.gov.

Текущая версия от 23:05, 31 августа 2015

Эта страница опубликована в основном списке статей сайта
по адресу 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)


F 01.png (1)


где:

— спектральная радиация, пришедшая на сенсор спутника;
DNcal — значения яркости пикселя сырого геоснимка;
Qcalmin — минимальное возможное значение пикселя геоснимка;
Qcalmax — максимальное возможное значение пикселя геоснимка;
LMINλ — минимальное значение спектральной радиации для конкретного датчика спутника для конкретного снимка;
LMAXλ — максимальное значение спектральной радиации для конкретного датчика спутника для конкретного снимка.


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


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


где:

— спектральная радиация, пришедшая на сенсор спутника;
DNcal — значения яркости пикселя «сырого» геоснимка;
Gainλ — усиление значения радиации относительно яркости исходного геоснимка;
Baisλ — смещение значения радиации относительно яркости исходного геоснимка.


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


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

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

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

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


Dos img 1.jpg

Рис. 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).


Dos f 03.png (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).

Dos f 05.png (5)


где:

L1% - коэффициент влияния угла падения и отражения солнечных лучей для 1% тёмного объекта;
d - расстояние от солнца до земли в астрономических единицах в конкретный день съёмки сцены на конкретной местности [7];
E0 - коэффициент солнечного внеатмосферного спектрального излучения (явно представлен как табличные данные и учитывается при калибровке датчиков Landsat 5 и 7 , для Landsat 8 дополнительно вычисляется);
θ - зенитное расстояние для солнца в радианах;
TZ - мера прохождения излучения от солнца до земли, в методе DOS2, принимается равным cosθ.


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


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


где:

Lλhaze — значение атмосферной дымки (hazing);
L1% — коэффициент влияния угла падения и отражения солнечных лучей для 1%-тёмного объекта;
Lλmin — спектральная радиация для 1%-тёмного объекта.


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


Dos f 07.png (7)


где:

ρλ — атмоферно скорректированные значения отражённой солнечной радиации;
— значения радиации, пришедшей на сенсор спутника;
Lλhaze — значение атмосферной дымки (hazing);
d — расстояние от Солнца до Земли в астрономических единицах в конкретный день съёмки сцены на конкретной местности [7];
E0 — коэффициент солнечного внеатмосферного спектрального излучения (явно представлен как табличные данные и учитывается при калибровке датчиков Landsat 5 и 7 , для Landsat 8 дополнительно вычисляется);
θ — зенитное расстояние для солнца в радианах;
TZ — мера прохождения излучения от солнца до земли, в методе 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), без учёта коррекции по углу падения и отражения солнечных лучей;
— каналоспецифичный мультипликативный расчётный фактор для (REFLECTANCE_MULT_BAND_x, где x — это номер канала) — усиление значения отражения;
— каналоспецифичный мультипликативный расчётный фактор для (REFLECTANCE_ADD_BAND_x, где x — это номер канала) — смещение значения отражения;
Qcal — значения яркости пикселя «сырого» геоснимка (DN).


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


Dos f 09.png (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).


Dos f 10.png (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