Организация региональных научно-исследовательских работ требует привлечения спутниковых данных, отражающих состояние местности и факторы развития природно-антропогенных систем. Спутниковые данные Landsat 8 OLI предоставляются пользователям уже после прохождения ими спектральной коррекции и получения географической привязки. Однако их применение в условиях расчлененного рельефа горно-таежной местности с крутыми (более 40 °) склонами затруднено ввиду неравномерной освещенности территории. Эти особенности изображения могут привести к неточностям в обработке изображений и дальнейшей интерпретации результатов.
Усилия, направленные на устранение этого дефекта – топографическую коррекцию, необходимо выполнять аналитику, работающему с изображениями, что требует внесения корректив в процесс исследования и приобретения соответствующих компетенций. Целью настоящей работы является адаптация алгоритма топографической коррекции спутниковых изображений SCS + C [1; 2] к использованию с данными Landsat 8 OLI и цифровой моделью рельефа (ЦМР) SRTM [3] и реализация его в виде автоматизированной программы Landsat Shadow Removal (LSR), написанной на языке Python с применением графического интерфейса QT.
Материалы и методы исследования
Применение методов топографической коррекции космических изображений обусловлено наличием неровностей земной поверхностей и несоответствия угла падения солнечных лучей и нормали, проведённой к земной поверхности (рис. 1).
Применяемый нами алгоритм SCS + C представляет собой доработанный метод SCS (sun-canopy-sensor) и отличается от него использованием C-параметра, позволяющего учитывать особенности рассеяния света и освещенности территории. Впервые данный алгоритм был предложен авторами [1]. Вычисление значений спектральной яркости каналов дистанционного изображения производится в соответствии с формулой (1):
(1)
где Ln – скорректированные значения спектральной яркости, L – исходные значения спектральной яркости, i – угол падения солнечных лучей, α – крутизна склона в градусах, θ – зенитный угол Солнца, С – параметр диффузной освещенности неба и вычисляемый как отношение (2):
(2)
где a – свободный член, b – угловой коэффициент уравнения линейной регрессии для зависимости L(cos(i)), сформулированной как (3):
L = a + b cos i. (3)
В соответствии с [2], работа алгоритма топографической коррекции SCS + C, в приложении к изображениям Landsat 8 OLI, включает стадии ввода изображений и их метаданных, предварительную обработку космических изображений и цифровой модели рельефа (digital elevation model, DEM), расчет C-параметра, применение топокоррекции, реализуемой для каждого канала в отдельности, и дальнейшую заверку результата (рис. 2).
Стоит отметить, что применение данных радиолокационной модели SRTM более релевантно применению космических данных, чем использование материалов топографических съёмок, так как учитывает высоту крон деревьев и кустарниковой растительности, т.е. отражающей поверхности, наблюдаемой на космической основе.
Рис. 1. Угловые величины, используемые в алгоритмах топографической коррекции (обозначения в тексте)
Рис. 2. Структура алгоритма SCS+C для реализации в виде программы ЭВМ (по [2], с дополнениями и изменениями)
Подготовленная на основе этого алгоритма программа Landsat Shadow Removal (LSR), представляемая с открытым исходным кодом на языке Python [4], использует в работе библиотеки numpy (работа с матрицами и многомерными данными), gdal/ogr (открытие geotiff и ESRI shape форматов), pandas (вывод статистической информации в виде таблиц MS Excel), scikit-image для трансформации SRTM-изображения к Landsat, а также библиотеки richdem и elevation (применяемые для расчета карты уклонов земной поверхности и азимута (экспозиции) склонов [5]).
Применение библиотек обработки геоданных gdal/ogr позволяет не только открывать популярные форматы файлов географических информационных систем, но и сохранять в этих форматах результаты обработки. Следует отметить, что использование программы подразумевает обязательное нахождение обрабатываемых данных в одной и той же метровой картографической проекции.
Визуальный интерфейс программы реализован с помощью библиотеки pyQT и обладает структурой, куда входят модули Input («Ввод»), Configure («Конфигурация») и Process («Обработка») (рис. 3).
С помощью модуля ввода, соответственно, выбирается каталог, содержащий космические снимки, каталог для вывода, векторный shape-файл, чей картографический охват (экстент) соответствует интересующей исследователя области. Shape-файл может быть сгенерирован в любом картографическом приложении. Если такой файл не был загружен пользователем до обработки, то используется охват первого из снимков в очереди.
Приложение использует каналы Landsat с разрешением 15 м, т.е. все имеющиеся каналы, исключая панхроматический канал 8, а также цифровую модель рельефа SRTM с разрешением 90 м, которая обрезается по охвату Landsat 8 OLI и чье разрешение приводится к Landsat за счет интерполяции.
После выбора каталога с исходными данными программа проверяет наличие там файлов и выбирает файлы с расширением имени tif. Файлы с таким расширением, находящиеся в одной папке, рассматриваются как относящиеся к одной спутниковой сцене и добавляются в очередь на обработку. Определение принадлежности файла к тому или иному спектральному каналу происходит на основе постфикса имени файла _b %N %., находящегося перед расширением tif, где %N % соответствует номеру спектрального канала.
Исходные данные должны быть в одной и той же координатной проекции и обладать существенным перекрытием или совпадением охвата [6].
Помимо обрезки изображений по рамке площади интереса (AOI) и применения топографической коррекции, программа может формировать изображения, применимые для различных целей. В модуле конфигурации можно выбрать рассчитываемые программой комбинации спектральных каналов (4-3-2, 7-4-2, 7-6-4, 6-5-2, 4-5-3, 6-4-2, 7-6-5) и спектральные индексы. Последние включают в себя нормализованный разностный индекс растительности (NDVI), окислы железа (IOA), гидроксильные изменения (HA), глинистые минералы (CM), кумулятивные изменения (CA) и главные компоненты изображения (PC).
Рис. 3. Интерфейс и модули приложения Landsat Shadow Removal (LSR)
Индекс кумулятивных изменений (CA – cumulative alterations) в соответствии с [7] вычисляется как композитное изображение RGB, в котором красному каналу соответствует HA, зелёному IOA, а голубому их полусумма. Этот продукт обработки спутниковых данных позволяет установить минералы – окислы железа и гидроокислы среди вторичных изменений субстрата.
Также в окне модуля конфигурации указывается азимут Солнца и угол его возвышения над горизонтом, что и необходимо для применения топографической коррекции освещенности территории. Ход выполнения программы, возникающие при этом уведомления, предупреждения и возможные ошибки выводятся в текстовом поле, находящемся на главном окне программы и в системной консоли.
Результаты исследования и их обсуждение
Для тестирования программы использованы спутниковые данные Landsat 8 OLI, сцена (105/029), дата съемки 27 октября 2015, отражающие природные и техногенные комплексы острова Кунашир (юг Большой Курильской гряды). Спутниковые изображения этой территории могут быть использованы для идентификации полей развития гидротермальных изменений, связанных с благородно- и редкометалльной минерализацией, температурных аномалий и молодых вулканических структур. Выбор сезона съемки для анализа территории острова Кунашир диктуется относительным уменьшением плотности растительного покрова, затушевывающего геологическое строение. Взаимная корреляция спектральных каналов Landsat 8 OLI для территории исследований приведена в таблице. Заметно, что между каналами с близкими длинами волн наблюдается более высокая корреляция, сохраняющаяся после применения алгоритма SCS + C.
В качестве подстилающих пород территории выступают неоген-четвертичные вулканические образования (андезиты, базальты, риолиты, дациты), а также осадочные и вулканогенно-осадочные образования, включающие морские террасы. На участках рельефа, находящегося выше снеговой линии, присутствуют скопления снега и льда [8].
Существенный перепад высот (до 600 м) вместе с относительно невысокими углами Солнца (31-32 °) на снимках приводит к появлению в изображении контрастных теней, проявленных как на спектрозональных изображениях, так и на отношениях каналов, образуя ложные аномалии, затрудняющие анализ.
Матрица корреляции использованных каналов Landsat 8 для сцены (105/029), дата съемки 27 октября 2015
band1 |
band2 |
band3 |
band4 |
band5 |
band6 |
band7 |
|
band1 |
1 |
0,998 |
0,975 |
0,970 |
0,523 |
0,155 |
0,187 |
band2 |
0,998 |
1 |
0,983 |
0,973 |
0,536 |
0,173 |
0,202 |
band3 |
0,975 |
0,983 |
1 |
0,991 |
0,660 |
0,320 |
0,342 |
band4 |
0,970 |
0,973 |
0,991 |
1 |
0,682 |
0,350 |
0,376 |
band5 |
0,523 |
0,536 |
0,660 |
0,682 |
1 |
0,857 |
0,841 |
band6 |
0,155 |
0,173 |
0,320 |
0,350 |
0,857 |
1 |
0,991 |
band7 |
0,187 |
0,202 |
0,342 |
0,376 |
0,841 |
0,991 |
1 |
Рис. 4. Сопоставление изображений до и после применения топокоррекции, а также их гистограммы. Обозначения: M – среднее значение для канала, S – стандартное отклонение для канала; цвета графиков соответствуют каналам RGB
Устранение теней и улучшение качества изображений является важной задачей, выполненной на этапе предварительной обработки материалов. Для эффективного применения программы LSR необходимо использовать метаданные снимка, включающие азимут Солнца и его положение над горизонтом (для условий съемки возвышение Солнца над горизонтом составляет 31,24 °, азимут светила 163,13 °). Сравнение фрагментов спектрозональных изображений (комбинированы каналы 7-4-2), полученных с помощью полнофункционального прототипа программы LSR с применением алгоритма SCS + C и без него, показано на рис. 4.
Сравнение результирующих изображений и их частотных кривых также свидетельствует о том, что изображение в общем стало светлее, стали менее контрастными или исчезли тени вдоль берега моря на северо-востоке участка. Левые пики, соответствующие недостаточно экспонированным участкам на кривых, сместились правее, правые пики, соответствующие ярким отражениям снега и льда, сохранили свое положение. Рельеф участков, занятых снегом и льдом, стал менее заметным.
В абсолютных значениях яркости для красного и зелёного каналов характерно увеличение среднего значения и уменьшение стандартного отклонения. Для голубого канала наблюдается обратная картина. Участок изображения, в пределах которого наиболее заметно осветление в условиях применения топографической коррекции при качественной оценке, демонстрирует усиление детализации изображения за счет проявления форм рельефа, ранее слабо различимых в затенённых участках.
Выводы
В результате выполнения исследования, направленного на адаптацию алгоритма SCS + C для обработки данных Landsat 8 OLI, получена альфа-версия описанного приложения Python, которое может быть выполнено в системах Windows и Unix. Программа обладает интуитивным графическим интерфейсом, облегчающим освоение ее использования новыми пользователями.
Алгоритмы, реализованные в программе LSR, были использованы для предварительной оценки территории острова Кунашир с помощью метода логистической регрессии (MaxEnt) [8]. Приложение топографической коррекции LSR было написано с использованием библиотеки gdal/ogr, хорошо зарекомендовавшей себя в подготовке свободного программного обеспечения для работы с пространственными данными и отличающейся высокой скоростью и надежностью. Открытость и кросс-платформенность исходного кода программы позволяют заинтересованным исследователям применять ее в своих проектах и расширять возможности. Проект также нуждается в тестировании на новых территориях с иными условиями рельефа и отражательными характеристиками подстилающего субстрата. Наибольшая отдача при использовании программы достигается применением современных географических информационных систем, таких как QGIS, задействованных для подготовки shape-файла области интереса и дальнейшего открытия и компоновки получаемых продуктов обработки спутниковых данных.