МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ
федеральное государственное автономное образовательное учреждение
высшего образования
Северный (Арктический) федеральный университет имени М.В. Ломоносова»
Кафедра информационных систем и технологий
(наименование выпускающей кафедры)
Карпов Александр Анатольевич
(фамилия, имя, отчество студента)
Институт
ИМИКТ
курс
группа
2
371466
09.04.02 Информационные системы и технологии
(код и наименование направления подготовки/специальности)
ВЫПУСКНАЯ КВАЛИФИКАЦИОННАЯ РАБОТА
(МАГИСТЕРСКАЯ ДИССЕРТАЦИЯ)
Разработка методики выявления и оценки площадей, пройденных лесными
пожарами, с использованием данных дистанционного зондирования
(наименование выпускной работы)
(обозначение)
Утверждена приказом №
__________
от «____» ___________ 20___ г.
Руководитель работы
Консультанты
Р.А. Алешко
Нормоконтроль
Рецензент
Зав. кафедрой
Р.А. Алешко
Н.В. Коновалова
Е.А. Деменкова
(подпись)
(дата)
(инициалы, фамилия)
Постановление Государственной экзаменационной комиссии от
Признать, что студент
______________
А. А. Карпов
(инициалы, фамилия)
выполнил и защитил выпускную работу с оценкой
Председатель ГЭК
Секретарь ГЭК
Архангельск 2016
МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ
федеральное государственное автономное образовательное учреждение
высшего профессионального образования
Северный (Арктический) федеральный университет имени М.В. Ломоносова»
Высшая школа информационных технологий и автоматизированных систем
(наименование института)
Кафедра информационных систем и технологий
(наименование выпускающей кафедры)
ЗАДАНИЕ ДЛЯ ПОДГОТОВКИ ВЫПУСКНОЙ
КВАЛИФИКАЦИОННОЙ РАБОТЫ
Студент 2 курса
Карпову Александру Анатольевичу
(фамилия, имя, отчество)
09.04.02 Информационные системы и технологии
(код и наименование направления подготовки/специальности)
Тема ВКР как Разработка методики выявления и оценки площадей,
пройденных лесными пожарами, с использованием данных дистанционного
зондирования жили были дпед жили были дпед жили были дпед жили были
дпед
Утверждено протоколом
заседания кафедры
№ _______ от «_______» _____________ 201
Срок сдачи выпускником законченной работы «____» __________ 201
Исходные данные к работе
г.
г.
Техническое задание
Основные разделы работы Анализ предметной области; получение исходных
данных и выбор инструментальных средств; анализ спутниковых снимков;
использование методики определения гарей на основе Google Earth Engine _а
____________________________________________________________________
____________________________________________________________________
____________________________________________________________________
Перечень подлежащих разработке вопросов Анализ предметной области;
получение спутниковых данных; первоначальный анализ снимков для
выявления гарей; создание методики выявления гарей на платформе Google
Earth Engine________________________________________________________
2
Перечень обязательных приложений к работе______________________________
_____________________________________________________________________
_____________________________________________________________________
_____________________________________________________________________
Перечень графического материала _______________________________________
_____________________________________________________________________
_____________________________________________________________________
_____________________________________________________________________
_____________________________________________________________________
База проведения исследований и внедрения результатов поиска______________
_____________________________________________________________________
_____________________________________________________________________
Консультанты по работе
по разделу __________________ ______ ________________________________________
должность, инициалы, фамилия
по разделу __________________ ______ ________________________________________
должность, инициалы, фамилия
по разделу __________________ ______ ________________________________________
должность, инициалы, фамилия
Дата выдачи задания « 6 » октября 2014 г.
Руководитель ВКР__________________
Р.А. Алешко
подпись
Задание принял к исполнению « 6 » октября 2014 г.
Подпись студента ________________________________
3
РЕФЕРАТ
Карпов Александр Анатольевич. Выпускная квалификационная работа
«Разработка методики выявления и оценки площадей, пройденных лесными
пожарами, с использованием данных дистанционного зондирования».
Руководитель выпускной квалификационной работы – Алешко Роман
Александрович.
Магистрантская работа. Пояснительная записка объемом 103 страниц,
содержит 50 рисунков, 12 таблиц, 22 источника, 2 приложения.
Ключевые слова: дистанционное зондирование, Landsat, лесной пожар, гарь,
Google Earth Engine.
Целью работы является разработка методики для выявления и оценки
площадей пройденными лесными пожарами.
Реализация проекта позволит произвести независимую оценку выгоревших
территорий и определить площадь ежегодно выгорающих территорий, использую
архивные данные. В качестве платформы для обработки снимков используется
Google Earth Engine.
В первой главе проведен анализ актуальности данной тематики и дан
краткий обзор существующих систем имеющих схожие функции.
Во второй и третей главе даны характеристики данных требуемых для
решения задачи и выбраны программные средства для проведения анализа и
создания автоматических методов.
В четвертой главе подробно рассмотрены создания методов и слоев
требуемых для решения задачи. Здесь также представлены полеченные результаты
и проведён их анализ, также рассмотрены перспективы развития решения
проблемы с использованием новых спутников
20 июня 2016
4
ОГЛАВЛЕНИЕ
ОПРЕДЕЛЕНИЕ, ОБОЗНАЧЕНИЕ, СОКРАЩЕНИЕ
7
ВВЕДЕНИЕ
8
1АНАЛИЗ ПРЕДМЕТНОЙ ОБЛАСТИ
9
1.1 Обоснования актуальности выбранной тематики
9
1.2 Обзор существующих систем
10
1.2.1 Global Fire Information Management System
1.2.1.1 Алгоритм обнаружения активных пожаров
10
11
MOD14/MYD14
1.2.1.2 Алгоритм детектирования сгоревших территорий MCD45
12
1.2.2 The European Forest Fire Information System
14
1.2.3 ИСМД-РосЛесХоз
15
1.2.4 ScanEx Fire Monitoring Service
18
2 ПОЛУЧЕНИЕ ИСХОДНЫХ ДАННЫХ И ВЫБОР
21
ИНСТРУМЕНТАЛЬНЫХ СРЕДСТВ
2.1 Получение исходных данных
21
2.2 Выбор инструментальных средств
22
3 АНАЛИЗ СПУТНИКОВЫХ СНИМКОВ
24
3.1 Предварительная обработка снимков
24
3.2 Анализ спектральных характеристик гарей
27
3.3 Классификация анализируемого изображения
35
3.4 Использование спектральных индексов для выявления гарей и
38
определения их площади.
3.5 Создания маски облачности и маски водных объектов.
44
3.6 Расчет точности масок гарей, построенных с помощью различных
47
спектральных индексов
4 ИСПОЛЬЗОВАНИЕ МЕТОДИКИ ОПРЕДЕЛЕНИЕ ГАРЕЙ НА ОСНОВЕ
51
СЕРВИСА GOOGLE EARTH ENGINE
4.1 Краткое описание Google Earth Engine
51
4.2 Создание безоблачного композита
52
4.2.1 Создание композита
52
5
4.2.2 Создания композита методом усреднения значения пикселей
54
4.2.3 Создание безоблачного композита, метод удаления облаков
55
4.2.4 Создание безоблачного композита с использованием нескольких
57
методов
4.3 Создание водной маски
60
4.4 Методы выявления выгоревших территорий
61
4.4.1 Метод выявления гарей по после-пожарному композиту.
61
4.4.2 Метод выявления гарей по до-пожарному и после-пожарному
64
композиту.
4.4.2.1 Описание метода
64
4.4.2.2 Описание метода поиска изменений Single-Channel Change
66
4.4.2.3 Полученные результаты
69
Detection
4.5 Анализ полученных результатов
72
4.6 Перспективы использования спутников Sentinel
78
ЗАКЛЮЧЕНИЕ
80
СПИСОК ИСПОЛЬЗОВАННЫХ ИСТОЧНИКОВ
81
ПРИЛОЖЕНИЕ А
85
ПРИЛОЖЕНИЕ Б
99
6
ОПРЕДЕЛНИЯ, ОБОЗНАЧЕНИЯ СОКРАЩЕНИЯ
ДДЗ – Данные дистанционного зондирования;
ГИС – Географическая информационная система;
FIRMS (англ. The Fire Information for Resource Management System) –
Система противопожарного и ресурсного менеджмента;
MODIS (англ. Moderate-resolution Imaging Spectroradiometer) –
спектрорадимотетр установленный на спутниках Terra и Aqua;
MOD14 – продукт MODIS, содержащий информацию о аномальных
температурах со спутника Aqua;
MYD14 – продукт MODIS, содержащий информацию о аномальных
температурах со спутника Terra;
MCD45– продукт MODIS, содержащий информацию о сгоревших площадях;
EFFIS (англ. The European Forest Fire Information System) – Европейская
противопожарная лесная информационная система;
NDVI – (англ. Normalized Difference Vegetation Index) – нормализированный
разностный вегетационный индекс;
SWVI – (англ. Short Wave Vegetation Index) – коротковолновый
вегетационный индекс)
NBR – (англ. Normalized Burned Ratio) – нормализированный индекс гари;
BAI – (англ. Burned Area Index); – индекс гари;
MIRBI – (англ. Mid-Infrared Burned Index) – индекс гари в среднем
инфракрасном диапазоне;
NIR – значение отражения в ближнем инфракрасном диапазоне спектра;
SWIRS (англ. Short Short Wave Infrared Short reflectance) – значение
отражения в ближнем коротковолновом инфракрасном диапазоне спектра.
SWIRL (англ. Short Wave Infrared Long reflectance) – значение отражения в
дальнем коротковолновом инфракрасном диапазоне спектра.
WRI – (англ. Water Ratio Index) Значения индекса WRI равный 1 или более
означает, что данный пиксель занят водным объектом.
API
–
(англ.
Application
Programming
программирования приложений);
7
Interface)
–
интерфейс
ВВЕДЕНИЕ
Проблема учета выгоревших лесных территорий сегодня не является такой
болезненной, как оперативное обнаружение очагов пожаров или выявление
незаконных рубок. Решение данной проблемы направлено на долгосрочную
перспективу, куда входит создание единой системы для эффективного управления
лесными ресурсами и оценки их запасов. Вторым фактором не приоритетности
данной задачи может быть занижение официальной статистики по выгоревшим
территориям. Когда данная методика могла бы сделать эти данными прозрачными.
Имея векторные контура гарей, можно отслеживать восстановление лесов на
местах пожарищ. Наличие такой информации, дает возможность проводить анализ
экосистем и делать прогнозы по их развитию. Это показывает высокую культуру
лесопользования, когда мы не только используем природные ресурсы, но
задумываемся о их восстановлении.
Данная проблема наиболее остро проявилась в 2010 и 2011 году. Из-за
аномально жаркого и засушливого лета на территории России произошло
множество лесных пожаров, и выгорели огромные лесные территории.
На сегодняшний день даже имеющиеся данные могут быть недоступны.
Например, в лесничествах всегда создавались схемы гарей. Сейчас эти данные
собираются с помощью приборов GPS, но в министерство природы и в Единый
Противопожарный Центр поступают данные только о выгоревшей площади.
Векторные данные не собираются.
Имея данные со спутника Landsat 5 TM возможно провести ретроспективный
анализ, и выявит гари с 1984 года. Это создаст данные для будущей системы, где
будет отслеживаться восстановление лесных ресурсов.
Создание методики для определения гарей это кирпичик в создание единой
системы природопользования.
8
1 АНАЛИЗ ПРЕДМЕТНОЙ ОБЛАСТИ
1.1 Обоснования актуальности выбранной тематики
Каждый год по завершении пожароопасного сезона в лесах региона
органами государственной власти и экологическими организациями проводится
учет выгоревших лесных территорий. Это необходимо для оценки последствий
природных пожаров, обновления информации о лесах, анализа эффективности
системы охраны лесов, выработки мер по её совершенствованию и планирования
природоохранных мероприятий.
Часто, размер ущерба, подсчитанный государственными структурами и
общественными экологическими и научными организациями, различается в
несколько раз. Такие расхождения вызывают споры о достоверности полученной
информации.
Необходима
единая
методика
оценки
площади
выгоревших
территорий, основанная на использовании дистанционных методов.
По данным Центра по проблемам экологии и продуктивности лесов РАН,
площадь гари на территории лесного фонда России в 5 раз превышает площадь
вырубки лесов, а размеры ежегодного ущерба от лесных пожаров соизмеримы с
величиной доходов от лесного хозяйства, а в отдельные годы значительно
превышают его. Исходя из этих данных, следует, что оценка площади и
местоположения территорий, пройденных лесными пожарами, является очень
актуальной задачей. Решение этой задачи непосредственно на местности может
привести к значительным временным и трудовым затратам. Использование данных
спутниковой съемки значительно упрощает процесс оценки гарей и, в ряде случаев,
значительно превосходит по точности оценку на местности [1, 2]
9
1.2 Обзор существующих систем
1.2.1 Global Fire Information Management System
FIRMS (The Fire Information for Resource Management System) была разработана
в университете Мэриленд с поддержкой Национального агентства США по
аэронавтике
и
космического
пространства
(NASA)
и
международной
продовольственной и сельскохозяйственной организацией ООН (UN Food and
Agriculture Organization). Данная система отражает информацию об активных
пожарах лишь с небольшой задержкой
по времени. Информация обновляется
каждые 8 часов.
Оперативная версия FIRMS известна, как Global Fire Information Management
System (GFIMS). На сайте данной системы можно найти информацию о пожарах и
архивные данные по ним. Существует большая потребность в данном ресурсе, в
особенности в условиях недостаточно отлаженной работы по мониторингу
пожаров работников служб, отвечающих за их обнаружение и тушение, в том
числе и в России.
К преимуществам использования информационной системы FIRM можно
отнести обзорность (данные предоставляются на весь мир, по России скачиваются
одним файлом), регулярность получения данных (несколько раз в день), точность
привязки на местности, независимость предоставляемой информации, легкость
использования пользователей сети Интернет, доступ к склейкам исходных снимков
на многие территории в удобном синтезе каналов.
Несомненно, есть и ограничения. Они в первую очередь связаны с низким
разрешением исходных снимков, автоматическими алгоритмами обработки и
задержкой предоставления получаемой информации, не позволяющей отслеживать
пожары в режиме реального времени. Система не позволяет отличить пожар от
любых других источников теплового излучения (на предприятиях, территориях
нефтедобычи и т.д.). Оперативные снимки MODIS, использующиеся для
мониторинга,
не
позволяют
детектировать
кратковременные, небольшие по площади пожары.
10
слабые,
низкотемпературные,
Рисунок 1.1 - Интерфейс GFIMS
Таким образом, система дает вполне качественную информацию о верховых
и сильных низовых пожарах. Однако для мониторинга некоторых торфяных и
травяных пожаров она не всегда удобна. [3]
Точки
аномальных
температур
возможно
скачать
с
EartData
(https://earthdata.nasa.gov/earth-observation-data/near-real-time/firms) сайта в формате
shp, kml, txt. Доступны как свежие данные, так и архивные данные начиная с 2006
года.
1.2.1.1Алгоритм обнаружения активных пожаров MOD14/MYD14
Для мониторинга активных пожаров используется стандартный продукт
MODIS Land MOD14/MYD14 (Fire and Thermal Anomalies). Данный продукт дает
информацию об очагах пожаров, т.е. о локациях, горящих в данный момент
(рисунок 1.2). Принцип детектирования пожаров основан на их сильном излучении
в среднем инфракрасном диапазоне.
Данная серия продуктов производятся на основе данных получаемых в
диапазонах 4 микрометра (каналы MODIS номер 21 и 22) и 11 микрометров (канал
MODIS номер 31) с разрешением 1 км. Для отделения облачности используется
вспомогательный
диапазон.
Стратегия
11
определения
пожара
основана
на
абсолютном
детектировании
пожара
(если
его
сила
достаточна
для
детектирования) и на детектировании относительно фоновых значений (для учета
изменчивости температуры поверхности и отражения солнечного света). Для
определения
ложных
срабатываний,
таких
как
солнечные
блики
или
незамаскированные участки береговой линии, а также для маскирования облаков
применяются дополнительные каналы 1 и 2 (250 метровые огрубленные до 1 км
разрешения, диапазоны 0.65 и 0.68 микрометров), также 500 метровый 7 канал
(диапазон 2.1 микрометра) и канал MODIS номер 32 (разрешение 1 км, диапазон 12
микрометров).
Рисунок 1.2 - Пример сгоревшей территории (черная область) и очагов MOD14A1
(красные залитые полигоны)
Terra MODIS получает данные дважды в день, также как и Aqua MODIS. Эти
четыре ежедневных наблюдения камерой MODIS помогают осуществлять
глобальный мониторинг процессов горения и помогают изучать влияние пожаров
на экосистемы, атмосферу и климат [4]
1.2.1.2 Алгоритм детектирования сгоревших территорий MCD45
Не вдаваясь в технические подробности, опишем основную идею,
используемую для выделения сгоревших площадей. Специально разработанный
для MCD45 алгоритм обнаружения сгоревших участков основан на анализе
временных серий ежедневных данных об отражательной способности поверхности.
Алгоритм использует двулучевую функцию отражательной способности (BRDF) и
12
позволяет детектировать только свежие территории, пройденные огнем, исключая
площади, сгоревшие ранее (например, в прошлом году или в прошлом месяце).
Основная идея заключается в ретроспективном анализе коэффициентов отражения
каждого пиксела и предсказании его последующего значения. В случае, если
предсказанное значение отличается от непосредственно наблюдаемого на
некоторую пороговую величину, то анализируемый пиксел считается пикселомкандидатом на отнесение к «сгоревшим». Для окончательного отнесения пиксела в
разряд «сгоревших» требуется, чтобы он прошел тест на временную устойчивость
(выполняется на основе данных об его отражательной способности в последующие
дни).
Продукт MCD45 доступен в двух форматах: HDF-EOS (официальный,
MCD45A1)
и
GeoTIFF
(подготовкой
данной
модификации
занимается
Мэрилендский университет).
MCD45A1 – ежемесячный продукт третьего уровня, представляющий собой
грид-сетку 500 метрового разрешения, находящийся в синусоидальной проекции и
содержащий попиксельную информацию о сгоревших площадях, начиная c 2000
года. Как и в случае прочих продуктов MODIS, имя файла MCD45A1 содержит
информацию о пространственном и временном охвате продукта, а также номер его
версии. [5]
Рисунок 1.3- Концептуальная схема обнаружения сгоревших площадей на основе
анализа коэффициента отражательной способности (ρ)
13
1.2.2 The European Forest Fire Information System
EFFIS (The European Forest Fire Information System) является комплексной
системой полного цикла от предупреждения о лесном пожаре и подготовке к его
ликвидации до выявления площади выгоревшей лесной территории и оценки
ущерба нанесенным пожаром. В данную систему входят 30 стран Европы и
Средиземноморского региона.
Danger forecast – прогноз опасности, forest fire events – лесные пожары, fire
detection – выявление пожаров, burnt area maps – картографирование выгоревших
территорий, land cover damage assessment – оценка пострадавших территорий,
emission assessment – оценка выброса CO2, Potential soil erosion estimate – оценка
эрозии почв, vegetation regeneration – восстановление растительности
Рисунок 1.4 - Цикл мониторинга пожаров в системе EFFIS
На рисунке 1.4 показаны различные модули EFFIS: прогнозирование
пожарной опасности, обнаружение пожаров, оценка выгоревших площадей, оценка
выбросов в атмосферу, оценка восстановления погоревших территорий.
Прогнозирование
пожарной
опасности
рассчитывается
по
двум
метеорологический, прогностических моделям French Météo-France и Deutsche
Wetter Diens и предоставляется прогноз на неделю вперед. Эти данные
14
Добавлено примечание ([УзМ1]): ??
используются для вычисления общего европейского индекса пожарной опасности,
основанного на Canadian Fire Weather Index (FWI).
Обнаружение лесных пожаров осуществляется с помощью данных,
полученных спектрометром MODIS, расположенных на спутниках Terra и Aqua.
Данные, полученные с MODIS обновляются ежедневно и доступны на веб-сервере
данной системы для пользователей в режиме реального времени.
EFFIS не стремится заменить системы мониторинга пожаров стран
европейского союза. Система разрабатывается главным образом для согласования
действий стран участников в случае крупного пожара. [6]
1.2.3 ИСМД-РосЛесХоз
«ИСДМ-Рослесхоз» является распределенной информационной системой,
состоящей нескольких блоков и подсистем, которые физически располагаются в
центральном (г. Пушкино Московской области) и региональных информационных
узлах (г. Москва, Новосибирск, Хабаровск, г. Красноярск. Имеется также
резервный узел в ИКИ РАН, на котором также тестируются модернизированные
элементы системы).
Главной особенностью «ИСДМ-Рослесхоз» является комплексный анализ
информации, связанной с лесными пожарами, которая поступает из разных
источников (метеорологическая информация, данные о результатах наземного и
авиационного мониторинга, поступающие от региональных лесопожарных служб и
данные космического мониторинга) [7]. При этом основой системы мониторинга
являются
данные
дистанционного
зондирования
Земли.
Это
позволяет
формировать однородную, независимую от человеческого фактора информацию о
лесных пожарах. К уникальной особенности системы можно отнести наличие
своеобразной обратной связи, когда в «ИСДМ-Рослесхоз» поступает информация о
подтверждении или опровержении сведений о загораниях, зарегистрированных из
космоса.
Официально утверждены ряд методик применения данных «ИСДМРослесхоз», в первую очередь для контроля достоверности сведений о пожарной
15
опасности в лесах и лесных пожаров [8]. Таким образом, «ИСДМ-Рослесхоз»
интегрирована в организационную систему охраны лесов от пожаров России.
Рисунок 1.5 - Структура ИСДМ-Рослесхоз
Вся информация (в том числе и обработанные ДЗЗ, поступающие в узлы
ИСДМ-Рослесхоз) структурируется, каталогизируется и архивируется в банк
данных и доступна из геоинформационной системы, интегрированной в webинтерфейс (рисунок 1.5). Выходные продукты в виде табличных отчетов или
географически привязанных композитных изображений в ГИС-формате с
атрибутивной информацией в любое время доступны зарегистрированным
пользователям, например, руководителю команды тушения лесного пожара в
любой точке Российской Федерации через сеть Интернет.
Основная информация, используемая системой:
1) Ежедневные информационные оперативные продукты, полученные на основе
данных спутников серии NOAA (прибор AVHRR) и спутников TERRA, AQUA
(прибор MODIS).
2) Ежедекадные информационные продукты для оценки последствий действия
пожаров, полученные на основе данных спутников TERRA, AQUA (прибор
MODIS) и SPOT (прибор VEGETATION).
16
3) Информационные продукты 1-2 раза в год на основе данных высокого и
среднего разрешения спутников SPOT, LandSat, Ресурс-ДК. [7]
Рисунок 1.6 - Интерфейс системы ИСДМ Рослесхоз
Возможности системы:
1) Прогнозирование и мониторинг пожарной опасности.
2) Детектирование и мониторинг лесных пожаров в динамике.
3) Оценка пройденной огнём площади.
4) Предварительная оценка повреждений насаждениям от пожаров (в том
числе выявление погибших насаждений).
5)
Сопоставление
данных
наземных,
авиационных
и
космических
ГИС-интерфейсе комплексной
информации
наблюдений, включающая обратную связь с регионами.
6)
Интеграция
в
одном
(топоосновы, ДЗЗ и атрибутивных данных) с целью поддержки управленческих
решений
в
области
мониторинга
лесопожарной
ситуации,
контроля
за
переданными субъектам полномочиями и оценки эффективности использования
субвенций.
7) Формирование данных по динамике изменений лесного фонда, не
связанной с воздействием лесных пожаров, для использования в качестве исходных
данных для систем мониторинга лесопользования, санитарного состояния лесов и
т.д. [9]
Технические характеристики системы:
17
1) Детектирование загораний на площади от 0,1 до 50 га (в зависимости от
метеоусловий и интенсивности горения);
2) Оперативность регистрации загораний в зависимости от частоты пролёта
спутников в различных регионах;
3) Точность определение координат в зависимости от широты точки
загорания;
4) Оценка пройденной огнём площади от 200 до 1000 га - с точностью до
30%, - более 1000 га с точностью до 5%. [9]
1.2.4 ScanEx Fire Monitoring Service
ScanEx Fire Monitoring Service является продуктом компании ScanEx
специализирующийся на выявлении лесных пожаров и оценке выгоревшей
территории.
Веб-сервис ScanEx Fire Monitoring Service предназначен для оперативного
получения спутниковых данных о местоположении очагов пожаров на территории
России и оценки площади гари и нанесенного ущерба. В веб-сервисе SFMS
отображается информация за последние четверо суток с ежедневным обновлением
данных по обстановке.
Очаги пожаров детектируются по данным спектрорадиометра MODIS,
который является ключевым инструментом на борту американских спутников Terra
и Aqua. Чувствительность приборов позволяет детектировать лесные и степные
пожары площадью от 1 га и более. В результате обработки геопривязанных
изображений можно локализовать положение горящих участков в пределах
области размером 1 х 1 км. Сервис SFMS обеспечивает возможность просмотра в
программе Google Earth растровых изображений, синтезированных в натуральных
цветах из данных оптических каналов MODIS, что позволяет оценить положение
облачного покрова и возможность обнаружения пожаров в интересующих районах
мониторинга. [10]
18
Рисунок 1.7– Сеть приемных станций компании ScanEx
Для расширения возможностей мониторинга в проекте используются
детальные мультиспектральные снимки со спутников SPOT 4 (разрешение 20
м/пиксел) и Landsat-5 (30 м/пиксел), которые позволяют наблюдать последствия
пожаров и картировать территории выгоревших участков местности. Детальная
информация на сервисе SFMS отображается в виде квиклуков всех пролетов SPOT
4 и Landsat-5, принимаемых станциями «УниСкан» сети ИТЦ «СКАНЭКС».
Пользователь может выбрать и заказать средствами сервиса малооблачные
детальные снимки районов мониторинга, на которых ранее датчиками MODIS
были обнаружены очаги пожаров.
Данные сервиса SFMS могут использоваться не только для обнаружения и
мониторинга пожаров, но и для наблюдения за другими «горячими» объектами,
например, факелами на нефтяных и газовых месторождениях.
Для доступа к сервису используется популярный веб-геосервис Google Earth.
Подключиться к информации о пожарах может каждый желающий, добавив в
программу Google Earth ссылку http://catalog.scanex.ru/sfms.kml. При обновлении
данной ссылки в интерфейс Google Earth загружается самая последняя информация
19
— данные о положении и свойствах детектированных очагов пожаров, обзорные
снимки МОDIS и квиклуки более детальной съемки SPOT 4 и Landsat-5. Данные
сгруппированы в разделы по четырем последним дням наблюдений, что позволяет
отслеживать как оперативную обстановку за текущий день, так и проводить
ретроспективный анализ изменения ситуации за трое предшествующих суток.
Кроме того пользователи могут добавлять собственные данные, интегрируя их с
уже имеющейся информацией о пожарах для совместного анализа и принятия
решений.
Доступ к данным MODIS и квиклукам SPOT 4, Landsat-5 на сервисе SFMS
бесплатный
возможность,
и
открытый.
настраивать
Для
зарегистрированных
пороговые
значения
пользователей
программы
есть
детектирования
«горячих» точек и по специальным ценам осуществлять в онлайн режиме заказ
снимков SPOT 4 и Landsat-5. [10]
Рисунок 1.8 - интерфейс сервиса SFMS
20
ПОЛУЧЕНИЕ
2
ИСХОДНЫХ
ДАННЫХ
И
ВЫБОР
ИНСТРУМЕНТАЛЬНЫХ СРЕДСТВ
2.1 Получение исходных данных
1). Изображения со спутников среднего разрешения Landsat 5 TM, Landsat 7
ETM+, Landsat 8 OLI являются оптимальными для решения магистерской задачи.
Снимки имеют каналы с видимым и инфракрасным спектром длин волн, что
является важным для выявления выгоревших территорий. Преимуществом
является большая полоса захвата равная 185 км, т.к дает возможность сделать
полное покрытие на Архангельскую область, используя всего несколько снимков.
Характеристики сенсоров, установленных на эти спутники, приводится в
приложении.
1) Данные снимки есть возможность скачать с сайта U. S Geological Survey.
(http://glovis.usgs.gov/). Периодичность съемки спутников Landsat 16 дней. Первые
снимки в архиве датируются с 1984 года, что дает возможность провести анализ
снимков, по которым уже имеется информация о пожарах. Большой архив снимков
дает возможность отследить динамику пожаров на определённые территории.
Выбор снимков производится по критериям, указанным в таблице 2.1.
Таблица 2.1 – Спецификация для выбора снимков
Критерий
Требуемый уровень
Разрешение снимка
Среднее (30 м на пиксель)
Дата съемки
Вегетационный период (желательно конец
августа)
Процент облачности
Не более 10 %
Уровень обработки
С произведенной геометрической коррекцией
Формат полученных данных
GeoTiff
21
Добавлено примечание ([УзМ2]): 2) ??
Рисунок 2.1 –Точки с аномальными температурами на территории
Архангельской области за 2011 год
2) Точки аномальных температур на территорию Архангельской области,
выявленные спектрометром MODIS, установленным на спутниках Terra и Aqua.
Данную информацию в векторном формате данных shp можно получить с сайта
FIRMS
(Fire
Information
for
Resource
Management
System)
(https://firms.modaps.eosdis.nasa.gov/download/request.php)
3) Снимки (SPOT, Worldview, QuickBird) высокого разрешения для
верификации полученных результатов.
4) Для сравнения полученных данных будут использоваться данные с сайтов
ИСМД-РОСЛЕСХОЗ
(https://nffc.aviales.ru/main_pages/open_map.shtml
и
Космоснимки-пожары (http://fires.ru/). Данные сайты содержат информацию по
площади выгоревших территорий и о датах пожаров.
22
2.2 Выбор инструментальных средств.
Для
выполнения
магистерской
диссертации
в
качестве
программного
обеспечения будет использоваться QGIS, как геоинформационная система. Данный
программный продукт имеют свободную лицензию.
Для расширения возможностей QGIS использовались следующие модули:
1) OpenLayers Plugin – модуль дает возможность использовать в проекте
карты и спутниковые снимки Google, Yandex, Bing, OpenStreetMap
2) Semi-Automatic Classification Plugin – модуль имеет возможность
скачивать снимки Landsat с координат заданных в проекте, производить
предварительную обработку и классификацию снимка.
QGIS имеет встроенный язык python, который дает возможности писать
скрипты для автоматизации шаблонных действий и создавать собственные модули,
используя встроенные команды.
ScanEx Image Processor будет использоваться для выполнения анализа
изображений и поиска изменений в снимках и получения результатов.
Для создания покрытия на всю Архангельскую область и проведения анализа
снимков в глобальном масштабе будет использован сервис Google Earth Engine.
23
3 АНАЛИЗ СПУТНИКОВЫХ СНИМКОВ
3.1 Предварительная обработка снимков
Для начала работы с космическими снимками следует произвести
предварительную
обработку,
которая
заключается
в
геометрической,
радиометрической коррекции, привязки снимка по опорным точкам к местности и
цветовая коррекция используемых снимков.
Снимки, скачанные с сайта Geological Survey Unite State, уже имеют данный
уровень обработки, поэтому можно перейти к следующему этапу.
На
данном
этапе
обработки
используется
модуль
Semi-Automatic
Classification Plugin для перевода значения пикселя на снимке в коэффициент
отражения
солнечной
энергии
и
выполнение
атмосферной
коррекции.
(рисунок 3.1).
Рисунок 3.1 - окно предварительной обработки плагина Semi-Automatic
Classification
Отражательная способность – отношение количества поступившего на
объект света и количество света им отраженного.
24
В
результате
операции
на
выходе
получается
канал,
содержащий
абсолютные безразмерные значения отражения. Преимущество таких данных в
том, что они не зависят от времени съемки, а зависят только от самого объекта. [11]
Расчет отражения солнечной энергии (TOA - Top of Atmosphere Reflectance)
осуществляется по формуле:
𝑄𝑝 =
π Lλ d2
(3.1)
Esunλ cosθ
где Lλ – спектральная плотность, количество приходящего на сенсор излучения, в
W/(m2 * ster *μm).;
d – расстояние от Земли до Солнца в астрономических единицах на дату съемки;
Esunλ – средняя солнечная внеатмосферная энергетическая освещенность является
постоянной для каждого канала (таблица 3.1), W/(m2 * µm);
сosθ – зенитный угол Солнца в момент съемки, содержится в MTL файле значение
SUN_ELEVATION [12]
Таблица 3.1 - Значения Esunλ для различных каналов
Канал
1
2
3
4
5
7
Landsat 5
1983
1769
1536
1031
220
83.44
Landsat 7
1997
1812
1533
1039
230.8
84.90
Следующим этапом нормализации данных является уменьшение влияния на
снимок атмосферы и перевод значений радиации, дошедшей до сенсоров спутника
(TOA radiance), в значения реально отражённого от земли спектрального излучения
солнечного света.
Влияние атмосферы на снимок проявляется в целом ряде факторов: угол
падения и отражения солнечных лучей, прозрачность атмосферы, газовый фактор и
дымка (рисунок 3.2).
25
Рисунок 3.2 - Факторы, влияющие на попадание отраженной солнечной
радиации на сенсоры спутника.
Dark Object Subtraction (DOS) – атмосферная коррекция, основанная на
анализе изображения. Данная методика впервые была предложена Chavez в
1996 г. [13] Суть метода состоит в нахождении яркости однопроцентного тёмного
объекта снимка с последующей коррекцией минимума значений каждого пикселя
изображения относительно спектральной яркости найденного объекта. Следует
отметить, что точность данных методик ниже, чем коррекции, основанной на
физических параметрах, но при отсутствии атмосферных измерений данная
методика может улучшить результат. [14]
𝐿𝑝 = 𝐿𝑚𝑖𝑛 − 𝐿𝐷𝑂1%
(3.2)
где Lmin – Спектральная плотность соответствует значению DN (Digital
Numbers), что от количества всех пикселей менее или равно 0,01% [15];
Ldo1% - спектральная плотность абсолютно черного тела, которая по оценке имеет
значения отражения 0,01.
Формула для снимков Landsat:
26
𝐿𝑚𝑖𝑛 = 𝑀𝐿 ∙ 𝐷𝑁𝑚𝑖𝑛 + 𝐴𝐿
(3.3)
Спектральная плотность абсолютно черного тела дана [15]:
𝐿𝐷𝑂1% =
𝐿𝑝 =
0,01∙[(𝐸𝑆𝑈𝑁𝜆∙𝑐𝑜𝑠𝜃𝑆 ∙𝑇𝑍 )+𝐸𝑑𝑜𝑤𝑛 ]∙𝑇𝑉
𝜋∙𝑑2
𝑀𝐿 ∙𝐷𝑁𝑚𝑖𝑛 +𝐴𝐿 −0.01∙[(𝐸𝑆𝑈𝑁𝜆 ∙𝑐𝑜𝑠𝜃𝑆 ∙𝑇𝑍) +𝐸𝑑𝑜𝑤𝑛 ]∙𝑇𝑉
𝜋∙𝑑 2
(3.4)
(3.5)
Существует несколько DOS методик (DOS1, DOS2, DOS3, DOS4),
основанные на различных значениях Tv, Tz и Edown. Наиболее простой методикой
является DOS1, где данные коэффициенты имеют следующие значение:
Tv = 1;
Tz = 1;
Edown = 0.
Формула спектральной плотности приводится к следующему виду:
𝐿𝑝 =
𝑀𝐿 ∙𝐷𝑁𝑚𝑖𝑛 +𝐴𝐿 −0.01∙𝐸𝑆𝑈𝑁𝜆 ∙𝑐𝑜𝑠𝜃𝑆
(3.6)
𝜋∙𝑑2
Значение отраженной энергии в этом случае равно:
𝑄𝑝 =
π Lλ−𝐿𝑞 d2
(3.7)
Esunλ cosθ
3.2 Анализ спектральных характеристик гарей
Первоначальный анализ спектральных характеристик гарей произведен по
снимку LT51750162011187KIS01 спутника Landsat 5 TM. Дата снимка 06.07.2011.
Местоположения снимка граница Архангельской области и республики Коми.
На снимке четко видны шлейфы дыма от лесных пожаров и свежие гари.
Помимо гарей 2011 года на снимки присутствую гари различных лет. Данная
информация отображена на рисунке 3.3.
Анализ производился с помощью модуля для QGIS Semi-Automatic
Classification по шести из семи каналов Landsat 5 TM. Шестой канал, который
отображает температуру поверхности земли, был исключен из классификации, т. к.
имеет больше величины значений и затрудняет построение графика.
27
Для проведения визуального анализа с целью поиска пострадавших от
пожаров территорий следует использовать следующие комбинации каналов 7-4-2
для Landsat 5,7 и 7-5-3 для Landsat 8. Эта комбинация дает изображение близкое к
естественным цветам, но в тоже время позволяет анализировать состояние
атмосферы и дым. Сухостойная растительность выглядит оранжевой, водаголубой. Песок, почва и минералы могут быть представлены очень большим
числом цветов и оттенков. Сгоревшие территории будут выглядеть ярко красными.
Эта комбинация используется для изучения динамики пожаров и пост-пожарного
анализа территории. Оливково-зеленый цвет характерен для лесных массивов и
более темный цвет является индикатором примеси хвойных пород.[16]
Рисунок 3.3 - Снимок LT51750162011187KIS0 с точками аномальных
температур и датами их обнаружения
28
Модуль расширяет возможности QGIS для работы со спутниковыми
снимками и позволяет производить их классификацию. Для этого требуется
отметить объекты на снимке при помощи полигонов с указанием класса и
подкласса. Все объекты были разделены на 5 классов: облако, растительность
(хвойные, лиственные, трава), не-растительность (почва, болота, вырубка), вода,
гарь (3 гари 2010 года и 1 гарь 2011 года). Всего создано 28 обучающих полигонов.
На
рисунке
3.4
приведен
пример
интерфейса
модуля
для
проведения
классификации.
Рисунок 3.4 - Создание обучающих полигонов для классификации.
По заданным обучающим полигонам модуль Semi-Automatic Classification
позволяет строить график отражения солнечной энергии в различных каналах.
Из графика на рисунке 3.5 видно, что наиболее близкими спектральными
характеристиками обладают объекты, отмеченные как хвойные, вырубки и болото.
Диапазоны значений отражения солнечной энергии для данных объектов
представлены на рисунке 3.6.
29
Рисунок 3.5 - Спектральный характеристики различных объектов
Рисунок 3.6 - Диапазоны значений отражения солнечной энергии от объектов
Из графика можно сделать вывод, что гари наиболее четко отделяются от
других объектов на снимке в 4 и 7 канале, а в 1,2 и 3 каналах имеют узкий диапазон
значений. Гарь 2011 года имеет более низкие значения отражения в 5 и 7 канале в
сравнении с гарями 2010 года.
Для проверки этого факта по снимку были дополнительно выделены
территории пострадавшие от огня в различное время и построен график
спектральных свойств гарей разных лет. Было отобрано 13 гарей: 2011г. – 2 шт.,
2010 г. – 4 шт., 2006 г. – 2 шт., 2004 г. – 1 шт., 2003 г. – 1 шт., 2001 г. – 1 шт.
30
По графику на рисунке 3.7 видно, что гари 2011 года имееют меньшие
значения отражения во всех каналах, наибольшая разница наблюдается в 4, 5 и 7
каналах. По таблице значений отражения гарей и отклонению в различных каналах
(рисунок
3.8)
полученной
из
модуля
Semi-Automatic
Classification
были
подобранны пороги для создания маски гарей. Были выбраны самые минимальные
и макимальные значения по каждому каналу и занесены в таблицу 3.2
Рисунок 3.7 - Спектральный характеристики гарей различных лет
Рисунок 3.8 - Значения отражения гарей 2011 и 2010 г. для различных каналов и
значения отклонения
31
Таблица 3.2 – Пороговые значения для выделения гарей
Минимальное значение
Максимальное значение
отражения
отражения
Канал 1
0,089
0,100
Канал 2
0,071
0,084
Канал 3
0,052
0,073
Канал 4
0,104
0,191
Канал 5
0,080
0,176
Канал 7
0,054
0,138
Наименование канала
Используя диапазоны из таблицы 3.2, была создана маска гарей.
Формула для растрового калькулятора:
0.089<=B1 && B1<=0.100 && 0.071<=B2 && B2<=0.084 && 0.052<=B3 &&
B3<=0.073 && 0.104<=B4 && B4<=0.191 && 0.080<=B5 && B5<=0.176 &&
0.054<=B7 && B7<=0.138 ? 1:0
Рисунок 3.9 – Маска территорий пострадавших от пожаров
Результат представлен на рисунке 3.8. Визуально видно, что маска
захватывает объекты, не являющиеся гарями, и содержит много шумов, так же не
вся площадь гари входит в маску.
32
Следовательно, нужно уменьшить диапазоны используемых значений для
получения маски лучшего качества. Из графика, представленного на рисунке 3.6
видно, что гари хорошо различимы от других объектов на снимке в 4 и 7 канале,
поэтому ориентируясь по графику, уменьшим диапазоны значений для этих
каналов. В 4 канале уменьшим максимальный порог до 0,160, а в 7 канале
увеличим минимальный порог до значения отражения 0,100.
Второй вариант формулы для растрового калькулятора:
0.089<=B1 && B1<=0.100 && 0.071<=B2 && B2<=0.084 && 0.052<=B3 &&
B3<=0.073 && 0.104<=B4 && B4<=0.160 && 0.080<=B5 && B5<=0.176 &&
0.100<=B7 && B7<=0.138 ? 1:0
Результат, полученный по второму варианту формулы представлен на
рисунке 3.10. Из данного рисунка видно, что шумов стало намного меньше, но
опять маска не покрывает полную площадь гари.
Рисунок 3.10 - Маска территорий пострадавших от пожаров
Для
улучшения качества создаваемой
маски был проведен анализ
спектральных характеристик гарей по дополнительным полигонам и проведен
анализ по графику спектральных кривых для получения более узких диапазонов
для уменьшения количества шумов. Данные представлены в таблице 3.3.
33
Таблица 3.3 - Пороговые значения для выделения гарей с учетом анализа
дополнительных полигонов
Минимальное значение
Максимальное значение
отражения
отражения
Канал 1
0,089
0,100
Канал 2
0,071
0,095
Канал 3
0,00
0,085
Канал 4
0,072
0,191
Канал 5
0,080
0,176
Канал 7
0,054
0,160
Наименование канала
Третий вариант формулы для растрового калькулятора:
0.089<=B1 && B1<=0.100 && 0.071<=B2 && B2<=0.095 && 0<=B3 &&
B3<=0.085 && 0.072<=B4 && B4<=0.191 && 0.080<=B5 && B5<=0.176 &&
0.054<=B7 && B7<=0.160 ? 1:0
Полученный результат представлен на рисунке 3.11. В данном варианте
маска полностью покрывает площади, пострадавшие от огня, но минусом является,
как и в первом варианте, что объекты, не относящиеся к гарям, попали в маску.
Далее был проведен анализ по другим снимкам, который дал худший
результат.
Следовательно, можно сделать вывод, что не существует постоянных
порогов для создания маски и значения отражения солнечной энергии от
поверхности Земли зависят от атмосферных условий и различаются даже для
снимков, сделанных с небольшим временным промежутком.
Причиной получения некачественного результата может быть подбор
пороговых значений вручную. Поэтому следует создать автоматический сбор
параметров с гарей по минимальным и максимальным порогам отражательных
характеристик. Данный алгоритм будет реализован в главе 4.
34
Рисунок 3.11 - Маска территорий пострадавших от пожаров
3.3 Классификация анализируемого изображения
Для снимка LT51750162011187KIS01 была произведена классификация,
используя тренировочные полигоны, с помощью которых мы собирали статистику
для построения спектральных кривых в первоначальном анализе изображения.
Классификация производилась с использованием плагина Semi-Automatic
Classification. Для классификации обычно выделяют четыре группы объектов:
застройка, растительность, почвы, вода. В нашем случае полигоны были созданы с
учетом этих групп кроме застройки, так как на снимке данная группа объектов не
была обнаружена. Дополнительный класс был создан для дыма от пожаров и
облаков, которые не входят в данную классификационную группу. В каждом
классе содержится подкласс. Например, в группу лес вошли подклассы хвойных,
лиственных и смешанных лесов. В группу почвы вошли болота, вырубки, гари и
открытые почвы. Данная структура классов входящих в макрокласс отображена на
рисунке 3.12.
35
Рисунок 3.12 - Схема классов входящих в макрокласс
Данный плагин позволяет использовать три методики для выполнения
классификации: метод наименьших расстояний, максимального сходства и
спектральной плотности.
В данной работе использовался алгоритм максимального сходства. Результат
классификации представлен на рисунке 3.13. Создав большое количество
тренировочных полигонов и достаточное количество подклассов можно с большой
точностью классифицировать изображение. Проблема заключается в том, что
спектральные характеристики, полученные с тренировочных полигонов, будут
подходить для классификации только данного изображения. Применить те же
полигоны для классификации снимка, полученного через год примерно в это же
время, не получится, так как результат классификации будет отличаться.
Следовательно, еще раз можно сделать вывод, что не существует постоянных
спектральных характеристик для выявления определенных объектов и создания
масок по постоянным пороговым значениям.
36
Рисунок 3.13 - Классификация снимка
Неточности классификации: тени на снимке классифицированы, как вода;
болота и дым от пожаров, так же попали в один класс; старая гарь и болота на
снимке классифицированы, как вырубка. Проводить улучшение классификации не
имеет смысла, так как данный метод не является универсальным.
.
37
3.4 Использование спектральных индексов для выявления гарей и
определения их площади.
Характерным
спектральная
признаком
отражательная
растительности
способность,
и
ее
состояния
характеризующаяся
является
большими
различиями в отражении излучения разных длин волн. Знания о связи структуры и
состояния растительности с ее спектрально отражательными способностями
позволяют
использовать
космические
снимки
для
картографирования
и
идентификации типов растительности и их стрессового состояния. Для работы со
спектральной информацией часто прибегают к созданию так называемых
«индексных»
изображений.
На
основе
комбинации
значений
яркости
в
определенных каналах, информативных для выделения исследуемого объекта, и
расчета
по
этим
значениям
«спектрального
индекса»
объекта
строится
изображение, соответствующее значению индекса в каждом пикселе, что и
позволяет выделить исследуемый объект или оценить его состояние. Спектральные
индексы, используемые для изучения и оценки состояния растительности,
получили общепринятое название вегетационных индексов [17].
Главным
достоинством
спектральных
индексов
является
их
универсальность. Мы можем использовать формулу спектрального индекса для
различных спутников, различных по времени снимков. Но если мы сравниваем
индексы различных годов, например, построенных по до-пожарному и после
пожарному снимку, то для получения достоверного результат требуется
использовать снимки желательно с одного спутника и снимки, полученные
примерно в одинаковое время года.
Для выявления выгоревших территорий на изучаемой области были
рассчитаны индексы NDVI, SWVI, NBR, BAI, MIRBI для каждого снимка. Для
определения контура и расчета площади гари была вычислена разность каждого
индекса dNDVI, dSWVI, dNBR, dBAI между снимком до пожара и после пожара.
NDVI - Normalized Difference Vegetation Index (Нормализированный
Разностный Вегетационный Индекс). Расчет NDVI базируется на двух наиболее
стабильных спектральных диапазонах – красном и ближнем инфракрасном. В
красной области спектра лежит максимум поглощения
38
солнечной радиации
хлорофиллом,
а
в
ближней
инфракрасной
области
находиться
область
максимального отражения солнечной радиации. Отношение этих показателей друг
к другу позволяет четко отделять и анализировать растительные объекты от прочих
природных объектов.
NIR−RED
NDVI =
(3.8)
NIR+RED
где NIR – значение отражения ближнем инфракрасном диапазоне спектра (4 канал
для Landsat 5 TM);
RED - значение отражения в видимом красном диапазоне спектра (3 канал для
Landsat 5 TM).
SWVI – Short Wave Vegetation Index (Коротковолновый вегетационный
индекс).
Средняя
область
инфракрасного
спектра
отражает
изменения
влагосодержания растения, в то время как в ближней области инфракрасного
спектра находится максимум отражения солнечного света для растений.
Совместное
использование
двух
диапазонов
повышает
точность
оценки
влагосодержания растения и компенсирует влияние оказываемой структурой листа
растения.
SWVI =
NIR−SWIR
(3.9)
NIR+SWIR
где NIR – значение отражения в ближнем инфракрасном диапазоне спектра; (4
канал для Landsat 5 TM);
SWIR - значение отражения в дальнем коротковолновом инфракрасном диапазоне
спектра. (5 канал для Landsat 5 TM).(7)
NBR – Normalized Burned Ratio (Нормализированный Индекс Гари).
Ближний инфракрасный спектральный диапазон чувствителен к структуре клеток
растительности, в то время как средний инфракрасный диапазон чувствителен к
влажности растений и имеет тенденцию к увеличению на открытых участках и
гарях.
39
NBR =
NIR−SWIR
(3.10)
NIR+SWIR
где NIR – значение отражения в ближнем инфракрасном диапазоне спектра; (4
канал для Landsat 5 TM);
SWIR - значение отражения в дальнем коротковолновом инфракрасном диапазоне
спектра. (7 канал для Landsat 5 TM).(8)
BAI – Индекс Гари (Burned Area Index).
BAI = (0.1+NIR)2
1
(3.11)
+(0.06+Red)2
где NIR – значение отражения солнечной энергии в ближнем инфракрасном
диапазоне спектра; (4 канал для Landsat 5 TM);
RED – значение отражение солнечной энергии в ближнем инфракрасном диапазоне
спектра; (3 канал для Landsat 5 TM);
MIRBI – Индекс гари в среднем инфракрасном диапазоне (Mid-Infrared
Burned Index)
MIRBI = 10 ∙ SWIRL − 9.8 ∙ SWIRS + 2.0
(3.12)
где SWIRS (Short Short Wave Infrared Short reflectance) - значение отражения в
ближнем коротковолновом инфракрасном диапазоне спектра. (7 канал для Landsat
5 TM);
SWIRL (Short Wave Infrared Long reflectance) - значение отражения в дальнем
коротковолновом инфракрасном диапазоне спектра. (5 канал для Landsat 5 TM).
Далее произведём выявление произошедших изменений между двумя
снимками используя разность индексов и приведем примеры полученных образцов
для dNDVI, dSWVI, dNBR, dBAI, dMIRBI.
Выявление изменений производим с помощью программы ScanEx Image
Processor с использование функции Single-Channel Change Detection. Как указано на
рисунке 3.14 задаем опорным растром (reference raster) до-пожарный индекс и
после пожарный индекс задаем сравниваемым растром (comparing raster), методы
выявления изменений указываем вычитание (subtraction) и ставим галочку
40
напротив создать побитовую маску, в которую войдут пикселе с вероятностью
изменений более 95% (Create binary raster with following probability threshold).
Рисунок 3.14 - Расчет dNDVI
На рисунке 3.16 показаны побитовые маски выявленных изменений с
помощью различных индексов. Анализ с использованием спектральных индексов
проводился
на
снимках
LT51820152010201KIS01
(2010
г.)
и
LT51820152011236KIS01 (2011 г.). Примерно в центре снимка расположен город
Онега. На рисунке 3.15 красным квадратом отмечено место расположения гари для
которой были созданы маски используя различные спектральные индексы
(рисунок 3.16).
41
Рисунок 3.15 - Снимок LT51820152011236KIS01 в естественных цветах
42
Рисунок 3.16 – Гарь в каналах 7-4-2 и маски, созданные с помощью различных
спектральных индексов
43
3.5 Создания маски облачности и маски водных объектов.
На обоих снимках присутствует облачность и территория занятая морем,
реками и озерами (водные объекты). Следовательно, требуется удалить территории
занятые облачностью и водными объектами из результата выявленных изменений.
Для этого требуется создать маску облачности для каждого снимка и водную маску
по одному из снимков и объединить все маски в одну. Объединенную маску нужно
будет вычесть из маски полученных гарей, чтобы отсечь ненужные территории из
дальнейшего анализа.
Маску
облачности
можно
построить
тремя
методами:
1)Наиболее простой метод подобрать порог по 3 каналу Landsat для
облачности и создать маску с помощью растровый калькулятора для данного
диапазона значений.
Рисунок 3.17 - Создание маски облачности по одному снимку методом подбора
порога для облачности
Как видно из рисунка 3.17 эмпирическим путем был подобран порог для
облачности от 0,1180 до 0,1295 по 3 каналу.
Формула для растрового калькулятора:
44
0.1180<B3<0.1295?0:1
2) Более сложный метод, который требует дополнительного безоблачного
снимка. Вычитаем из снимка, для которого строим маску безоблачный снимок и
если разность больше 0,03, то пиксель относится к облачности.
Формула для растрового калькулятора:
(B3_1-B3_2)>0.03?1:0
Рисунок 3.18 - Создание маски облачности по двум снимкам
Как видно из рисунка 3.18 помимо облачности маска захватывает дымку и
контур облака захватывается полностью. Использую данную маску можно более
качество отсечь территории занятые дымкой и облачностью.
3) С помощью программы grass, входящий в комплекс QGIS и i.landsat.acca.
Данная команда автоматическую оценку облачности по алгоритму ACCA.
Для построения маски водных объектов можно воспользоваться двумя
методами:
1) Наиболее простой метод, как и для маски облачности, подобрать порог по
3 каналу Landsat для водных объектов и создать маску с помощью растровый
калькулятора для данного диапазона значений.
45
2) Для создания водной маски использовать индекс WRI – Water Ratio Index.
Значения индекса WRI равный 1 или более означает, что данный пиксель занят
водным объектом.
𝑊𝑅𝐼 =
GREEN+RED
(3.13)
NIR+SWIR
Где GREEN – зеленый спектральный диапазон (2 канал Landsat 5 TM);
RED – красный спектральный диапазон (3 канал Landsat 5 TM);
NIR – близкий инфракрасный диапазон (4 канал Landsat 5 TM);
SWIR – средний инфракрасный диапазон (5 канал Landsat 5 TM)
Для создания маски по индексу WRI вводим в растровый калькулятор
следующую формулу: (b2+b3)/(b4+b5)>= 1 ? 1:0 [18]
На рисунке 3.19 представлена маска водных объектов, рассчитанная по
снимку LT51820152011236KIS01, черным выделены водные объекты.
Рисунок 3.19 - Водная маска, построенная с помощью индекса WRI
46
3.6 Расчет точности масок гарей, построенных с помощью различных
спектральных индексов.
Матрица ошибок представляет собой инструмент, использующий кросстабуляцию для показа того, как соотносятся значения совпадающих классов,
полученные из различных источников. В качестве источников могут выступать,
например, проверяемый растр (тематическая классификация) и опорный более
точный источник данных (растр или набор полевых данных в виде точек). При
интерпретации результатов обычно полагается, что проверяемый результат
потенциально является неточным, а проверочный растр хорошо отражает реальную
ситуацию. В противном случае, если проверочный растр также несовершенен,
нельзя говорить об «ошибке», а следует говорить о «разнице» между двумя
наборами данных. Для построения матрицы могут использоваться все ячейки
растра (пиксели) или выборка ячеек, расположенных случайно, стратифицировано
случайно или согласно какому-либо другому распределению.
Рисунок 3.20 - Маски гарей построенные с помощью различных индексов
47
По одной из осей матрицы записываются названия классов легенды
классификации проверяемого набора данных, по второй - классы легенды данных,
используемых для проверки. [19]
По снимку RapidEye с разрешением 5 м на пиксель вручную была создана
маска пострадавшей от пожара территории. Создан вектор по контурам гари и
затем произведена растеризация вектора. Данная маска использовалась как
эталонная для расчета ошибки масок полученных с помощью разных спектральных
индексов.
Рисунок 3.21 - Расчет точности маски NBR
Расчет ошибки матрицы производился в QGIS через плагин Semi-Automatic
Classification
(рисунок
3.21).
Точность
спектрального индекса NBR равна 86,58 %.
48
маски
построенной
с
помощью
Таблица 3.4 - Точности масок, построенных с помощью спектральных индексов.
Маска, рассчитанная по индексу
Точность, %
NDVI
93,88
NBR
86,58
MIRBI
59,14
SWVI
29,25
BAI
0,70
Полученная точность маски построенной с помощью индекса NDVI больше
точности маски построенной с помощью NBR. Но как видно из рисунка 3.22 маска
построенная по NDVI имеет больше ошибочных пикселе, чем построанная по NBR.
Если задать эталланой маску по NBR и сравнить с NDVI, то получается ошибка в
93%. Следовательно для получения контура гари более точным спектральным
индексом будет NBR.
Рисунок 3.22 - Сравнения масок построенных с помощью NDVI и NBR
49
Рисунок 3.23 - Сравнения масок построенных по индексам NBR и NDVI
(черное – соответствие, красное – несоответствие)
50
4 ИСПОЛЬЗОВАНИЕ МЕТОДИКИ ОПРЕДЕЛЕНИЕ ГАРЕЙ НА ОСНОВЕ
СЕРВИСА GOOGLE EARTH ENGINE
4.1 Краткое описание Google Earth Engine
Google Earth Engine – облачные сервис являющийся платформой для
обработки геопространтсвенных данных в масштабе Земли. Сервис включает в
себя
большой
объем
бесплатных
спутниковых
снимков
и
обеспечивает
высокопроизводительные, параллельные вычисления для проведения анализа
снимков, имеет API для JavaScript и Python.
API поддерживает проведение сложного геопространственного анализа,
включающего картографическую алгебру, работу с массивами, обработку
изображений, классификацию, выявления изменений, проведение временного
анализа, используя серию спутниковых снимков, конвертацию растра в вектор,
сбор статистики и построение графиков.
Преимуществами
данного
приложения
являются
высокая
скорость
обработки снимков, возможность обрабатывать снимки не скачивая их на
персональный компьютер, что сокращает время, потраченное на обработку данных
в разы, включает готовые алгоритмы, позволяющие создавать безоблачное
покрытие, используя снимки определённого временного диапазона.
Интерфейс веб-версии IDE для Earth Engine представлен на рисунке 4.1.
Earth Engine Code Editor использует синтаксис языка JavaScript для написания
скриптов для обработки данных.
Для выполнения магистерской работы Google Engine будет использоваться
для составления покрытия на территорию Архангельской области из снимков с
наименьшей облачностью и далее для расчета спектрального индекса NBR для
данного покрытия. Сравнивая индексы NBR различных покрытий, будет
выполнено выявление произошедших изменений. Для отсечений территории
занятой облаками и воды будет построена общая маска, включающая в себя
облачность для каждого покрытия, и маску воды. Данная маска требуется для
удаления облаков и водной территории из побитового слоя выявленных изменений,
так как это может привести к ошибке.
51
Zoom – увеличение, Geometry Tools – инструменты геометрии, Asset
Manager – менеджер загрузок, Script manager – менеджер скрипов, API
documentation – ИПП документация, Search for data –поиск данных, Imports –
импорт, Get a link(URL) to the script – получить ссылку на скрипт, Run the scrip –
запустить скрипт, Help button – помощь, Task manager – менеджер задания, Console
output – строка вывода, Inspect location, pixel values, object added to the map –
проверка положения, значения пикселя, объекты добавленные на карту, Layer
manager – менеджер слоев
Рисунок 4.1 - Интерес Earth Engine Code Editor
4.2 Создание безоблачного композита
4.2.1 Создание композита
Следующим важными шагом для выполнения проекта является создание
безоблачного композита на период после пожарного сезона (с 1 августа по 1
сентября). В качестве примера создадим композит на Архангельский регион.
Следующий скрипт создает коллекцию снимков Landsat 5 TM в значениях
отражения солнечной энергии от земной поверхности и фильтрует снимки по
52
выбранному периоду времени и границам Архангельской области. Далее создается
композит из всех снимков и вырезается по границам Архангельской области.
Композит отображается на экран в естественных цветах (рисунок 4.2). [20]
Скрипт для создания композита на Архангельскую область:
var Arhkangelsk_region =
ee.FeatureCollection('ft:1KJkBqYUi8_J2_qnl4rsGdxagmK1a28HQz6J4wzaJ');
var L5 = ee.ImageCollection('LANDSAT/LT5_L1T_TOA')
.filterBounds(Arhkangelsk_region)
.filterDate('2011-8-1', '2011-9-1')
.mosaic()
.clip(Arhkangelsk_region);
Map.addLayer(L5, {bands: ['B3', 'B2', 'B1'], max: 0.3}, 'composite 2011');
Map.setCenter(42, 63, 6);
Рисунок 4.2 - Композит, созданный из снимков Landsat 5 TM за период с 1 августа
по 1 сентября 2011 года
53
4.2.2 Создания композита методом усреднения значения пикселей
Следующий
накладывающихся
метод
создает
снимков.
композит
Метод
усредняя
значения
пикселей
ee.Algorithms.Landsat.simpleComposite()
используется для создания безоблачных композитов. [21]
Данный метод дал, на первый взгляд, неплохой результат, создав покрытие с
меньшей облачностью (рисунок 4.3). Но при увеличении масштаба можно увидеть
следы на месте удаленных облаков, так как удаление облаков производилось
посредством усреднения значения пикселя. Использование такого композита
может привести к неправильным результатам выявлением произошедших
изменений.
Рисунок 4.3 - создание безоблачного композита, используя метод
ee.Algorithms.Landsat.simpleComposite, справа - следы от облаков при
использовании метода за небольшой период времени
54
Пример скрипта:
var Arhkangelsk_region =
ee.FeatureCollection('ft:1KJkBqYUi8_J2_qnl4rsGdxagmK1a28HQz6J4wzaJ');
var L5 = ee.ImageCollection('LANDSAT/LT5_L1T');
var composite_2011 = ee.Algorithms.Landsat.simpleComposite({
collection: L5.filterDate('2011-8-1', '2011-9-1'),
asFloat: true});
composite_2011 = composite_2011.clip(Arhkangelsk_region)
Map.addLayer(composite_2011, {bands: ['B3', 'B2', 'B1'], max: 0.3}, 'composite
2011');
Map.setCenter(42, 63, 5);
4.2.3 Создание безоблачного композита, метод удаления облаков
В создании безоблачных композитов отлично показал себя метод удаления
облаков из снимков. Данный метод реализован на стандартной функции Google
Earth Engine ee.Algorithms.Landsat.simpleCloudScore(). Суть данного метода
заключается в удалении облаков из всех снимков, содержащихся в коллекции.
После выполнения удаления облачности из снимков создается композит и места
удаленных облаков заполняются снимками, которые не содержат облачности в
данном месте.
Данный метод имеет свои преимущества и недостатки. Преимуществом
является, то что пиксели не усредняются по значениям. К недостаткам можно
отнести неполной удаление облаков (рисунок 4.4) при установке фильтра
облачности на значения менее 20%, а при установке фильтра на 0% (полное
удаление облачности) удаляются большая часть композита и территория, имеющая
спектральные характеристики, которые можно отнести к минералам (песок). Также
при удалении облаков остаются контуры облака.
55
Рисунок 4.4 - создание безоблачного композита, используя метод
ee.Algorithms.Landsat.simpleCloudScore, справа - не все облака были удалены с
помощью данного метода
Пример скрипта:
var Arhkangelsk_region =
ee.FeatureCollection('ft:1KJkBqYUi8_J2_qnl4rsGdxagmK1a28HQz6J4wzaJ');
var L5 = ee.ImageCollection('LANDSAT/LT5_L1T_TOA')
.filterBounds(Arhkangelsk_region)
.filterDate('2011-8-1', '2011-9-1');
var fDeleteClouds = function(image) {
var cloud = ee.Algorithms.Landsat.simpleCloudScore(image);
var score = cloud.select(['cloud']).lte(20);
56
return image.updateMask(score);
};
var composite_2011 = L5.map(fDeleteClouds)
.mosaic()
.clip(Arhkangelsk_region);
Map.addLayer(composite_2011, {bands: ['B3', 'B2', 'B1'], max: 0.3}, 'composite
2011');
Map.setCenter(42, 63, 5);
4.2.4 Создание безоблачного композита с использованием нескольких
методов
Для создания более эффективного фильтра, следует использовать несколько
методов:
1) Метод ee.Algorithms.Landsat.simpleCloudScore(), который описан в главе
4.3.3
2) Для удаления остатков облачности и контуров облаков каждый пиксель
композита будет сравниваться с безоблачным композитом, созданным с помощью
стандартного метода ee.Algorithms.Landsat.simpleComposite(). Так как будет
использоваться период съемки за 10 лет для создания безоблачного композита, то в
данном случае данный метод создает качественный композит без следов
полностью неудаленных облаков.
Есть небольшой минус данного метода, т.к. удаляется территории, которые
по 3 каналу Landsat схожи с облачной (минералы), но территории пострадавшие от
огня не трогаются.
57
Рисунок 4.5 - Метод создания безоблачного композита, испозуя два метода
удаления облаков
Пример скрипта:
var Arhkangelsk_region =
ee.FeatureCollection('ft:1KJkBqYUi8_J2_qnl4rsGdxagmK1a28HQz6J4wzaJ');
var L5 = ee.ImageCollection('LANDSAT/LT5_L1T');
var free_clouds = ee.Algorithms.Landsat.simpleComposite({
collection: L5.filterDate('2000-1-1', '2010-9-1'),
asFloat: true});
var collection_2011 = ee.ImageCollection('LANDSAT/LT5_L1T_TOA')
.filterBounds(Arhkangelsk_region)
.filterDate('2011-8-1', '2011-9-1')
.sort('CLOUD_COVER',false);
//Функуция удаления облаков
var fDeleteClouds = function(image) {
var cloud = ee.Algorithms.Landsat.simpleCloudScore(image);
58
var score = cloud.select(['cloud']).lte(20);
return image.updateMask(score);
};
var composite_2011 = collection_2011.map(fDeleteClouds)
.mosaic()
.clip(Arhkangelsk_region);
var clouds_2011 = composite_2011.select('B3').subtract(free_clouds.select('B3'));
var cloud_mask_2011 = clouds_2011.expression('float(b("B3"))>0.03 ? 1:0');
var other_mask_1 =
cloud_mask_2011.expression('b("constant")>0?0:1').clip(Arhkangelsk_region);
var composite_2011 = composite_2011.updateMask(other_mask_1);
Map.addLayer(composite_2011, {bands: ['B3', 'B2', 'B1'], max: 0.3}, 'composite
2011');
Map.setCenter(42, 63, 6);
Рисунок 4.6 - Покрытие, созданное с использованием нескольких методов
59
4.3 Создание водной маски
Создание водной маски подробно описано в главе 3.5. Для создания водной
маски будет использоваться индекс WRI. Результат представлен на рисунке 4.7.
Водная маска создается по безоблачному композиту по снимкам с 2000 по 2010
год. Целью создания данной маски является отсечение бликов на воде, которые
могут попасть в композит рассчитанного индекса NBR и остаться при создании
маски гарей как шумы.
Рисунок 4.7 - Пример созданной водной маски и вычитание водной маски из
композита и композита индекса NBR.
60
4.4 Методы выявления выгоревших территорий
4.4.1 Метод выявления гарей по после-пожарному композиту.
Данный метод является наиболее простым и легко настраиваемым методом.
Его можно применить и для поиска других объектов на снимке, которые имеют
уникальные спектральные характеристики и легко отличаемы от иных объектов.
Суть метода заключается в создании полигонов для сбора статистики по
минимальному и максимальному значению пикселя в каждом канале снимка. Сбор
значений производится в шести каналах Landsat 5 TM с 1 по 5 и 7 канал с
исключением 6 термального канала. Данные значения используются как
минимальные и максимальные пороги для создания растровой маски гарей,
которая затем векторизуется и скачивается. Пороги, полученные по послепожарному композиту 2011 года приведены в таблице 4.1.
Таблица 4.1 - Пороги значений для определения гарей по композиту 2011
Канал
Минимальные значения
Максимальные значения
B1
0.07783
0.10242
B2
0.05162
0.07934
B3
0.03964
0.07635
B4
0.04287
0.12631
B5
0.04278
0.19052
B7
0.05289
0.14547
Преимущество данного метода заключается в простоте и скорости
выполнения скрипта. Из визуальной оценки работы скрипта можно сделать вывод,
что данный метод отлично выделяет гари. После получения вектора гарей можно
отсеять площади менее 3000 кв. м, чтобы уменьшить ошибку и отсеять шумы.
Алгоритм работы данного алгоритма и его результат представлен на рисунке 4.8.
Недостатком данного метода является выделение, как гарей на исследуемого
года, так и гарей прошлого и позапрошлого года. Есть возможность отделить гари
на исследуемый года, используя точки аномальных температур. Эта информация
61
является продуктом спектрометра MODIS, установленном на спутниках Terra и
Aqua, и содержится в отдельном наборе данных FIRMS-T21 в сервисе Google Earth
Engine. Но полученный результат будет не точным, т.к. пожар может быть закрыт
облаками и в этом случае информация о точках аномальных температур будет не
будет собрана.
Рисунок 4.8 - Метод выделения гари по после-пожарному композиту
Рисунок 4.9 - Полученные результаты методом анализа после-пожарного
композита.
62
Для получения площадей гарей полученная маска была векторизована в
QGIS для расчета площади пострадавших территорий. Полученные данные
представлены на рисунке 4.9.
Таблица 4.2 - Выявленная площадь гарей по композиту 2011 года по лесничествам
Лесничества
Площади, га
Архангельское
2460,53
Березниковское
172,89
Вельское
152,27
Верхнетоемское
6,82
Вилегодское
6,96
Выйское
248,70
Емецкое
195,71
Карпогорское
171,08
Котласское
92,50
Красноборское
2319,49
Лешуконское
8084,82
Мезенское
2379,89
Няндомское
0,65
Обозерское
7797,71
Онежское
9552,80
Пинежское
0,79
Приозерное
101,97
Пуксоозерское
8,27
Северодвинское
8626,11
Сурское
890,01
Холмогорское
4460,60
Шенкурское
641,90
Яренцкое
3304,02
по всем лесничествам:
51676,60
63
Рисунок 4.10 - Лесничества Архангельской области в градиентной заливке, где
зелёный цвет соответствует меньшей площади гари, красный максимальной.
Полученные данные за 2011 год
4.4.2 Метод выявления гарей по до-пожарному и после-пожарному
композиту.
4.4.2.1 Описание метода
Данный используется до-пожарный и после-пожарный композит для поиска
гарей. Используя данный метод можно с высокой точностью выявить гари на
исследуемый год, не учитывая гари за прошлые годы.
Для поиска гарей было создано два композита, используя снимки на период
августа за 2010 и 2011 год. Для облегчения анализа правильности работы
64
алгоритма из композитов были вырезаны только территории с точками аномальных
температур. Для этого был загружен набор данных FIRMS и выделен канал T21.
Далее, была произведена векторизация растра точек аномальных температур и
построен буфер дистанцией в 3 км, чтобы избежать обрезки гари.
Для каждого композита был рассчитан индекс NBR и проведен поиск
изменений методом single channel change detection, используя индексные
изображения. Для построения маски по значениям вероятностных изменений был
задан порог 95%. После получения маски для удаления шумов из маски были
удалены пикселы, которые имели значение NBR в после-пожарном композите
менее 0,2. После выполнения всех вышеописанных операций была произведена
векторизация маски и из вектора были удалены значения менее 3000 кв. м. для
отсеивания шумов. Схема работы данного алгоритма приведена на рисунке 4.10.
Рисунок 4.11 - Методика выявления гари по до-пожарному и после-пожарному
композиту
65
4.4.2.2 Описание метода поиска изменений Single-Channel Change Detection
Стандартными методами обнаружения изменений являются вычитание
яркостей одного снимка из яркостей другого. К положительным сторонам такого
простого способа поиска изменений, как вычисление разности, относится простота,
устойчивость и понятная интерпретация. К недостаткам - предположение о
линейной пропорциональности изменения яркостей искомым изменениям для всех
объектов на снимке.
Для успешной работы алгоритмов необходимо, чтобы обрабатываемые
растры удовлетворяли следующим требованиям:
1) Растровые каналы должны быть хорошо совмещенными.
2) Рекомендуется использовать изображения, полученные в близкие по
сезону даты
3) Обрабатываемые изображения должны иметь один базис (иметь
одинаковый размер в пикселях, пространственное разрешение и координату
верхнего левого угла)
Для уменьшения зависимости результатов от качества фотометрической
коррекции строится линейная регрессия анализируемого растра на опорный (по
данным в пределах заданной области анализа). Результатом выполнения операции
является растровый слой, значения яркости которого, приведенные в диапазон
значений 0-1.
Поиск изменений методом вычитания:
R = A-Ba
(4.1)
где R–промежуточный результат,
A–референсный растр,
Ba–результат регрессии анализируемого растра на референсный
Промежуточный результат используется для вычисления распределения, на
основании которого вычисляется вероятность отклонений. Для значений R,
соответствующих медиане распределения, вероятность изменений считается
нулевой, для значений, соответствующих вероятности 0.05 и 0.95 вероятность
66
изменений считается равной 90%, и т.д. Полученные значений вероятности
используются для формирования результирующих растров.
Рисунок 4.12 - демонстрация реализация алгоритма single channel change detection.
Создание маски облака по 3 каналу Landsat 5 TM.
Математическое описание алгоритма выявления изменений:
𝑥𝑖 – яркость i-го пикселя изображения A.
𝑦𝑖 – яркость i-го пикселя изображения B.
Пусть
𝑥𝑖 = 𝑎 ∙ 𝑦𝑖 + 𝑏 + 𝐸𝑖
67
Где a,b – коэффициенты линейной регрессии
Ei-случайная ошибка
Найдем a и b:
𝐸𝑖 = 𝑥𝑖 − 𝑎 ∙ 𝑦𝑖 − 𝑏
𝐸𝑖2 = (𝑥𝑖 − 𝑎 ∙ 𝑦𝑖 − 𝑏)2
Условие для нахождения a и b.(сумма квадратов ошибок должна быть
минимальной)
2
𝑆 = ∑𝑁
𝑖=1 𝐸𝑖 = 𝑚𝑖𝑛
𝑑𝑆
{𝑑𝑎
𝑑𝑆
𝑑𝑏
=0
=0
{
∑(−2) ∙ (𝑥𝑖 − 𝑎 ∙ 𝑦𝑖 − 𝑏) ∙ 𝑦𝑖 = 0
∑(2) ∙ (𝑥𝑖 − 𝑎 ∙ 𝑦𝑖 − 𝑏) ∙ 𝑦𝑖 = 0
Из второго уровня:
∑ 𝑥𝑖 − 𝑎 ∙ ∑ 𝑦𝑖 − ∑ 𝑏𝑖 = 0
∑ 𝑏 = 𝑁𝑏,
Где N – число пикселей
𝑏=
1
𝑁|
∙ (∑ 𝑥𝑖 − 𝑎 ∙ ∑ 𝑦𝑖 )
Подставим b в первое уравнение, заменим i на индекс j, чтобы избежать
путаницы:
1
∑(𝑥𝑖 − 𝑎 ∙ 𝑦𝑖 − ∙ (∑ 𝑥𝑗 − 𝑎 ∙ ∑ 𝑦𝑗 )) ∙ 𝑦𝑖 = 0
𝑁
1
1
𝑁
𝑁
∑ 𝑥𝑖 ∙ 𝑦𝑖 − 𝑎 ∙ ∑ 𝑦𝑖2 − ∙ ∑ 𝑥𝑗 ∙ ∑ 𝑦𝑖 + ∙ 𝑎 ∙ ∑ 𝑦𝑗 ∙ ∑ 𝑦𝑖 = 0
1
1
𝑁
𝑁
∑ 𝑥𝑖 ∙ 𝑦𝑖 − ∙ ∑ 𝑥𝑗 ∙ ∑ 𝑦𝑖 = 𝑎 ∙ ∑ 𝑦𝑖2 − ∙ 𝑎 ∙ ∑ 𝑦𝑗 ∙ ∑ 𝑦𝑖 = 0
𝑎=
∑ 𝑥𝑖 ∙ 𝑦𝑖 −
∑ 𝑦𝑖2 −
1
∙ ∑ 𝑥𝑗 ∙ ∑ 𝑦𝑖
𝑁
1
∙ ∑ 𝑦𝑗 ∙ ∑ 𝑦𝑖
𝑁
Вернемся к индексу i:
𝑎=
𝑏=
∑ 𝑥𝑖 ∙ 𝑦𝑖 −
∑ 𝑦𝑖2 −
1
𝑁
1
∙ ∑ 𝑥𝑖 ∙ ∑ 𝑦𝑖
𝑁
1
∙ ∑ 𝑦𝑖 ∙ ∑ 𝑦𝑖
𝑁
∙ (∑ 𝑥𝑖 − 𝑎 ∙ ∑ 𝑦𝑖 )
Где ∑ −это ∑𝑁
𝑖=1
Ba-изображение, что есть результат регрессии
68
Zi- яркость пикселя изображения Ba
𝑍𝑖 = 𝑎 ∙ 𝑦𝑖 + 𝑏
R – промежуточный результат (матрица)
ri – элемент матрицы R, 𝑟𝑖 𝜖[−255; 255] или 𝑟𝑖 𝜖[−1; 1]
в зависимости от минимального и максимального значения яркости
𝑟𝑖 = 𝑥𝑖 − 𝑧𝑖
из матрицы R можно построить ненормированное распределение
Рисунок 4.13 - ненормированное распределение
nv- число элементов матрицы R имеющих значение v
𝑣=255
∑𝑣=1
𝑣=−1 𝑛𝑣 = 𝑁 или ∑𝑣=−255 𝑛𝑣 = 𝑁
N – площадь фигуры у нормированного распределения равна 1
Площадь части фигуры:
𝑊(𝑉) = ∑𝑣𝑣=−1 𝑛𝑣
Вероятность в точке V:
𝑁
2
2
𝑁
𝑃(𝑉) = |𝑊(𝑉) − | ∙
Связь матриц P и R:
𝑟
𝑁
2
2
𝑁
𝑖
𝑝𝑖 = |∑𝑣=−1
𝑛𝑣 − | ∙
4.4.2.3 Полученные результаты
Используя разработанный метод обнаружения пострадавших от пожаров
территорий, были получены результаты, представленные в таблице 4.3. Анализ
проводился по снимкам Landsat 5 TM для 2007, 2010 и 2011 года. В анализе не
участвовали снимки 2000-2006 и 2008 и 2009 годов, так как на сервисе Google Earth
Engine нет снимков за эти годы. Для получения результатов рассчитывался индекс
NBR для до-пожарного и после-пожарного композита, потом эти слои скачивались
и обрабатывались в ScanEx Image Processor. Проводился поиск изменений через
функцию «Single channel change detection», полученная маска пропускалась через
растровый калькулятор, где фильтровался полученный результат. Пиксели в маске,
69
построенной с помощью «Change detection» не должны превышать значения 0,2 в
после пожарном растре NBR.
Таблица 4.3 - Площади гарей на территории Архангельской области за 2007, 2010,
2011 гг.
Год
2011
2010
2007
Площади, га
33540
5227
252
Рисунок 4.14 - Лесничества Архангельской области в градиентной заливке, где
зелёный цвет соответствует меньшей площади гари, красный максимальной.
Полученные данные за 2011 год по 2 методу
70
Визуальный анализ показал правдоподобность данных за 2011 год.
Полученные данные за 2010 и 2007 год имеют худшее качество. Возможная
причина такого результата некачественные, засвеченные снимки и недостаточное
покрытие за 2009-2010 и 2006-2007 год (покрытие строится из снимков, сделанных
в период августа). Фильтрация пикселей по значению менее 0,2 по послепожарному NBR растру для выделения гарей, так же была настроена по 2011 году.
Возможно, Landsat 5 TM мог иметь другие настройки в эти года и, следовательно,
требуется брать иной порог фильтрации. Улучшить качество покрытия так же
возможно расширив временной диапазон для снимков, например, брать снимки не
только за август, но и за сентябрь.
Таблица 4.4 - Площади гарей по лесничествам на территории Архангельской
области за 2011 гг.
Лесничество
Площади, га
Архангельское
1169,08
Березниковское
111,92
Вельское
0,51
Верхнетоемское
28,98
Выйское
513,55
Емецкое
126,08
Каргопольское
2,11
Котласское
327,98
Красноборское
3017,96
Лешуконское
5224,37
Мезенское
1105,61
Обозерское
6072,15
Онежское
6931,19
Пинежское
0,58
Плесецкое
1,60
Приозерное
28,40
Северодвинское
3694,77
Сурское
58,76
Устьянское
6,61
Холмогорское
2567,52
Шенкурское
141,41
Яренцкое
300,68
По всем лесничествам
31431,82
71
4.5 Анализ полученных результатов
Все полученные результаты оценки площадей гарей занесены в таблицу 4.6.
Благодаря этому можно сравнить официальные данные, полученные от Единого
Противопожарного Центра с продуктом MODIS MCD45 и двумя методами
дешифровки гарей которые были разработаны в магистерской диссертации.
Официальные данные имеют большие числа выгоревших территорий, чем
остальные. Это можно объяснить несколькими факторами:
1) Не все гари могли быть зафиксированы, т.к. композиты были созданы с
использованием снимков за период августа. По причине меньшей
облачности в это время. Лесные пожары могут продолжаться до октября.
2) Не все пожары можно зафиксировать, используя снимки Landsat. Низовые
пожары как правили, плохо выделяются на снимках. Пример представлен на
рисунке 4.15, где должна быть гарь на 523 га. Первый снимок за период
августа. На снимке виден активный пожар с характерным шлейфом дыма.
Далее, чтобы квартировать гарь на данном участке был взят снимок за
сентябрь, но на нем видно лишь небольшие полосы с характерными для гари
оттенками цветов. На следующем снимке красными точками отмечена маска
гари, созданная по 1 методу. Последний рисунок – это схема гари
предоставленная Вельским лесничеством.
Рисунок 4.15 - Не распознанная крупная гарь
72
3) Второй метод (использование до-пожарного и после-пожарного композита)
показал меньшие площади из-за отсутствия полного покрытия. Например,
центр композита Архангельской области на 2010 г. полностью был покрыт
облаками. Отсутствие данных 2010 года делает невозможным выявление
гарей в 2011 году по этому методу.
В сравнении с MCD45 данные, полученные первым и вторым методом за 2011
год больше в несколько раз. Хотя данные MCD45 получены с мая по октябрь. Это
объясняется склонностью данного алгоритма недооценивать выгоревшую площадь
и разрешения пикселя в этом наборе данных равно 500 м, что так же плохо влияет
на точность результатов.
В таблице 4.6 и далее метод с использованием только после-пожарного
композита будет называться 1 методом, а метод с использованием до-пожарного и
после-пожарного композита – 2 методом.
В таблице 4.5 были собраны данные по лесничествам за 2011 год полученные из
официальных источников и по 2-м методам.
90000
80000
70000
60000
Официальные данные
50000
1 метод
2 метод
40000
MCD45
30000
20000
10000
0
2011
2010
2009
2007
Рисунок 4.16 - Площади выгоревших территорий по годам
73
Таблица 4.5 - Сравнение полученных результатов на 2011 год с официальными
данными по лесничествам.
Площадь, га
Наименование
1 метод по точкам
лесничества
Официальная
аномальных
2 метода
температур
1. Архангельское
2388,28
2460,53
1169,08
2. Березниковкое
403,91
172,89
111,92
3. Вельское
1176,93
152,27
0,51
4. В-Тоемское
159,47
6,82
28,98
5. Вилегодское
35,02
6,96
6. Выйское
930,93
248,70
513,55
7. Емецкое
307,56
195,71
126,08
8. Каргопольское
2,11
9. Карпогорское
315,1
171,08
10. Кенозерское
11.Коношское
153,28
12.Котласское
207,255
92,50
327,98
13.Красноборское
4348,3
2319,49
3017,96
14.Лешуконское
9657,5
8084,82
5224,37
15.Мезенское
1240,11
2379,89
1105,61
16.Ненецкое
17.Няндомское
19,3
0,653
18.Обозерское
14528,3
7797,71
6072,15
19.Онежское
13932,9
9552,80
6931,19
20.Пинежское
153,85
0,79
0,58
21.Плесецкое
83,4
1,6
22.Приозерное
11517,1
101,97
28,4
23.Пуксоозерское
108,1
8,27
24.С-Двинское
6021,89
8626,11
3694,77
25.Соловецкое
2,1
890,01
26.Сурское
1747,91
58,76
27.Устьянское
21,62
6,61
4460,60
28.Холмогорское
5599,23
2567,52
641,90
29.Шенкурское
590,17
141,41
3304,02
30.Яренское
3960,96
300,68
31.Сийский
4,91
лесопарк
Всего
79615,39
51676,61
31431,82
Для сравнения результатов контуров выгоревших территорий были
запрошены схемы гарей у нескольких лесничеств. Удалось получить схемы
Вельского лесничества. Следовательно, результаты работы 1 методики можно
74
сравнить с полевыми данными. Результаты 2 методики сравнить не получится так
как из-за облачности в композите 2010 года на Архангельскую область
отсутствуют данные на эту территорию.
Таблица 4.6 - Сводная таблица выгоревших площадей из различных источников
Год
Площадь по
Площадь Площадь
Площадь,
Площадь
официальным по
по 1
подтвержденная по 2
данным
продукту методу, га
по аномальным методу, га
лесничеств,
MODIS
точкам, га
га
MCD45,
га
2011
79615
24339
72613
52996
33540
2010
14209
2095
6982
5888
5227
2009
183
327
1772
2
2007
929
10699
5
252
Таблица 4.7 – Официальные данные, полученные с Вельского лесничества,
Вельского участкового лесничества.
№ квартала
Официальная площадь, га Площадь по 1 методу, га
86
7,6
5,3
109
24,0
21,2
118
23,0
21,1
124
11,0
8,9
75
Рисунок 4.17 - Сравнения материалов по результатам пожаров 2011 года Вельского
лесничества с результатом выделения гарей по после-пожарному композиту
76
Рисунок 4.18 - Сравнения результатов полученных контуров с данными MCD 45
(MODIS). Средний рисунок – первый метод. Нижний – второй метод
77
4.6 Перспективы использования спутников Sentinel
Европейское Космическое Агентство (ESA) создает новое семейство
спутников Sentinels для развития программы «Copernicus». Каждая группа
спутников из семейства Sentinel состоит из двух аппаратов, что увеличивает
покрытие получаемых данных со спутников за счет более частой периодичности
съемки и увеличивает достоверность данных.
Данные спутники будут оснащены радарными и мультиспектральным
оборудованием для мониторинга поверхности суши, океана и атмосферы.
Спутники семейства Sentinels:
- Миссия Sentinel-1 будет представлять собой группировку из двух радарных
спутников на полярной орбите, оснащенных радаром с синтезированной апертурой
(SAR) для съемок в С-диапазоне. Съемки радарных спутников Sentinel-1 не будут
зависеть от погоды и времени суток. Миссия будет обеспечивать многие сервисы
Copernicus,
например,
мониторинг
покрытых
льдом
арктических
морей,
картографирование ледовых полей, мониторинг нефтяных разливов и обнаружение
кораблей с целью обеспечения безопасности, мониторинг подвижек земной
поверхности, картографирование лесов, внутренних вод и почв, поддержка
гуманитарных операций и управления кризисными ситуациями. Sentinel-1A,
запущен в апреле 2014 г., Sentinel-1B запущен в апреле 2016 г. Характеристике
спутников приводятся в приложении.
- Миссия Sentinel-2 оснащена мультиспектральным оборудованием высокого
разрешения для мониторинга поверхности Земли: мониторинг растительности,
почв, внутренних водных путей и прибрежных зон. Sentinel-2 так же предоставляет
информацию
для
экстренных
служб.
Спутники
находятся
на
полярной
орбите. Sentinel-2A был запущен в июне 2015 г., Sentinel-2B будет запущен во
второй половине 2016 г. Характеристика спутника приводится в приложении.
- Миссия Sentinel-3 имеет широкое предназначение для измерения морского
рельефа, измерения поверхности земли и моря, определения цвета земной и
морской поверхности с высокой точностью и надежностью. Миссия будет
поддерживать
прогнозирования
состояния
78
океана,
а
также
мониторинг
окружающей среды и климата. Sentinel-3A был запущен в феврале 2016 г. Sentinel3B планируется запустить в 2017 г.
- Миссия Sentinel-4 – атмосферный мониторинг с курсом на Meteosat Third
Generation-Sounder (MTG-S) спутника на геостационарной орбите.
- Миссия Sentinel-5 мониторинг атмосферы на полярной орбите на с помощью
MetOp Second Generation.
- Миссия Sentinel-5 Precursor атмосферный мониторинг, уменьшить разрыв данных
между Envisat и запущенного Sentinel-5.
- Миссия Sentinel-6 измерение глубины моря, изучение океанографии и климата.
[22]
Данные получаемые со спутников Sentinel-1 и Sentinel-2 являются
перспективными для выявления лесных территорий, пострадавших от огня.
Спутник Sentinel-2 имеет схожие диапазоны с Landsat 8 OLI, но обладает лучшим
пространственным разрешением в каналах с визуальным диапазоном 10 метров на
пиксель и имеет меньшую периодичность съемки. Sentinel-1 предоставляет данные
в радиометрическом диапазоне, что делает возможным выявление гарей без учета
облачности.
Данные, полученные с уже запущенных спутников находятся в открытом
доступе
и
могут
быть
скачаны
с
сайта
Copernicus
programme
(https://scihub.copernicus.eu). Снимки, сделанные спутниками доступны на сайте
после первоначальной обработки через 24часа.
79
ЗАКЛЮЧЕНИЕ
В ходе работы над выпускной квалификационной работой были созданы 2
методики выявления пострадавших от огня территорий:
1) Определение гарей только по после-пожарному композиты, используя
пороговые значения отражения солнечной энергии в различных каналах.
2) Определение гарей по до-пожарному и после-пожарному композиту,
выявляя изменения по спектральному индексу NBR
Для второго метода был выбран спектральный индекс NBR, т.к. данный
индекс более точно определяет площадь выгоревшей территории. Для этого была
произведена верификация полученных результатов по спутниковым изображениям
высокого разрешения.
Результаты были верифицированы со схемами гарей, полученных от
Вельского лесничества. Контуры гарей совпадают с данными схемами. Площадь
гарей, полученных 1 методом, отличается от данных лесничества в меньшую
сторону. Не совпадения результатов возможно как из-за неточности измерений
лесничеств, так и невысокого разрешения Landsat 5 TM (30 метров – 1 пиксель).
Верифицировать 2 метод по тем же схема не удалось из-за облачности в композите
2010 года на данной территории.
Были получены данные по выгоревшим территориям на Архангельскую
область за 2007, 2009, 2010 и 2011 годы.
Данные методики автоматически производят создание векторного слоя по
контурам гарей. Методики можно применить к любому другому региону России
или мира.
Перспективы развития - использование данных спутников Sentinel 1 и
Sentinel 2. Sentinel 1 использует радиодиапазон, что является преимуществом для
составления покрытия, так как можно не учитывать облачность. Sentinel 2 имеет
большее разрешения, большую полосу захвата и периодичность съемки в 5 дней,
что позволит создать полностью безоблачный композит за месяц.
Использование платформы
Google Earth Engine является
огромным
преимуществом, т.к. позволяет обрабатывать большие объемы данных с
применением облачных вычислений за короткое время.
80
СПИСОК ИСПОЛЬЗОВАННЫХ ИСТОЧНИКОВ
1) Барталев С.А., Егоров В.А. Детектирование сгоревших территорий РФ в 2010:
результаты
ИКИ
РАН,
электронный
источник,
режим
доступа:
http://gis-
lab.info/qa/fires-iki.html
2) Барталев С.А., Егоров В.А., Лупян Е.А. Адаптивный пороговый алгоритм
детектирования повреждений растительности пожарами на основе многолетней
статистической «нормы» сезонной динамики коротковолнового вегетационного
индекса // Седьмая всероссийская открытая ежегодная конференция "Современные
проблемы дистанционного зондирования Земли из космоса". Москва. ИКИ РАН. 1620 ноября 2009. Сборник тезисов конференции, 2009. С.272
3) GIS-LAB [Электронный ресурс] / Режим доступа свободный: http://gislab.info/qa/firms.html (Дата обращения: 10.06.2016)
4) GIS-LAB [Электронный ресурс] / Режим доступа свободный: http://gislab.info/qa/mod14a1.html (Дата обращения: 10.06.2016)
5) GIS-LAB [Электронный ресурс] / Режим доступа свободный: http://gislab.info/qa/mcd45.html(Дата обращения: 10.06.2016)
6) Jesús San-Miguel-Ayanz, Comprehensive Monitoring of Wildfires in Europe:The
European Forest Fire Information System (EFFIS), статья в журнале
7) Е. А. Лупян, С. А. Барталев, Д. В. Ершов, Р. В. Котельников, И. В. Балашов, М.
А. Бурцев, В. А. Егоров, В. Ю. Ефремов, В. О. Жарко, К. А. Ковганко, П. А.
Колбудаев, Ю. С. Крашенинникова, А. А. Прошин, А. А. Мазуров, И. А. Уваров, Ф.
В. Стыценко, И. Г. Сычугов, Е. В. Флитман, С. А. Хвостиков, П. П. Шуляк.
Организация работы со спутниковыми данными в информационной системе
дистанционного мониторинга лесных пожаров Федерального агентства лесного
хозяйства (ИСДМ-Рослесхоз) (рус.) // Современные проблемы дистанционного
зондирования Земли из космоса. — 2015. Т. 12 № 5. С. 222-250.
8) Приказ Министерства природных ресурсов и экологии Российской Федерации от
05.12.2014 № 540 «Об утверждении Методики осуществления оперативного
контроля за достоверностью сведений о пожарной опасности в лесах и лесных
пожарах»
81
9) Официальный сайт ИСДМ-Рослесхоз [электронный ресурс] / Режим доступа
свободный:
http://www.aviales.ru/default.aspx?textpage=25
(Дата
обращения:
10.06.2016)
10) Н. В. Евтушенко, В. В. Морозов Геоинформационный сервис спутникового
мониторинга пожаров // Экология производства. - 2010. - N
11) СканЭкс инженерно-технический центр руководство пользователя ScanEx
Image Processor v 4.0, Москва 2013
12) GIS-LAB [Электронный ресурс] / Режим доступа свободный: http://gislab.info/docs/revised_landsat5_tm_radiometric_calibration_procedures_and_postcalibrati
on_dynamic_ranges.pdf (Дата обращения: 10.06.2016)
13) - Chavez, P. S. (1996). Image-based atmospheric correction—revisited and improved.
Photogrammetric Engineering and Remote Sensing, 62(9), 1025 – 1036.
14) Semi-Automatic Classification [Электронный ресурс] / Режим доступа
свободный: http://fromgistors.blogspot.com/p/theinterface-2.html (Дата обращения:
10.06.2016)
15) Direction Magazine [Электронный ресурс] / Режим доступа свободный:
http://www.directionsmag.com/entry/land-cover-classification-of-cropland-a-tutorialusing-the-semi-automa/376137 (Дата обращения: 10.06.2016)
16) GIS-LAB [Электронный ресурс] / Режим доступа свободный: http://gislab.info/qa/landsat-bandcomb.html (21.08.2016)
17) . Зятькова Л.К., Елепов Б.С. У истоков аэрокосмического мониторинга
природной среды («Космос» – программе «Сибирь»): монография.– Новосибирск:
СГГА, 2007. – 380)
18) GIS-LAB [Электронный ресурс] / Режим доступа свободный: http://gislab.info/qa/grass7-landsat8
processing.html#.D0.A1.D1.80.D0.B0.D0.B2.D0.BD.D
0.B5.D0.BD.D0.B8.D0.B5_.D1.80.D0.B5.D0.B7.D1.83.D0.BB.D1.8C.D1.82.D0.B0.D1.
82.D0.BE.D0.B2_.D1.81_.D1.80.D0.B0.D1.81.D1.87.D0.B5.D1.82.D0.B0.D0.BC.D0.B
8_.D0.B2_ENVI_.D0.BF.D0.BE.D1.81.D0.BB.D0.B5_FLAASH (20.06.2015)
19) GIS-LAB [Электронный ресурс] / Режим доступа свободный: http://gislab.info/qa/error-matrix.html (20.06.2015)
20) Google Earth Engine [Электронный ресурс] / Режим доступа свободный:
https://developers.google.com/earth-engine/ic_filtering (Дата обращения: 10.06.2016)
82
21) Google Earth Engine [Электронный ресурс] / Режим доступа свободный:
https://developers.google.com/earth-engine/landsat (Дата обращения: 10.06.2016)
22)Copernicus program [Электронный ресурс] / Режим доступа свободный:
http://www.esa.int/Our_Activities/Observing_the_Earth/Copernicus/Overview4
обращения: 10.06.2016)
83
(Дата
ПРИЛОЖЕНИЕ A
Исходный код 1 метод:
var reg = 29;
var year = 2011;
var Russia =
ee.FeatureCollection('ft:1z7BRawAQRGo2kznrwVBCD74UuPJIICrvZ
PEbGIfA');
//Фильтруем регион по номеру региона
var region = Russia.filterMetadata('name', 'equals', reg);
var kvartalka =
ee.FeatureCollection('ft:1pXC2DUx7z67VRAebCqbdMJzA3y8AD3e5X8nvRBJ');
var MCD45 = ee.ImageCollection('MODIS/051/MCD45A1')
.filterBounds(region)
.filterDate('2011-05-01', '2011-10-01');
var HotSpotsKML =
ee.FeatureCollection('ft:1i1_y8xKqQqNCZj4jXtw9hNNRbaGymmiU2
x1NejNw');
var HotSpots = HotSpotsKML.filterMetadata('YEAR', 'equals',
year);
//Создание безоблачного композита
var L5 = ee.ImageCollection('LANDSAT/LT5_L1T');
var free_clouds = ee.Algorithms.Landsat.simpleComposite({
collection: L5.filterDate('2000-7-1', (year-1)+'-9-1'),
asFloat: true});
//Создание после-пожарного композита
var collection_2 =
ee.ImageCollection('LANDSAT/LT5_L1T_TOA')
.filterBounds(region)
85
.filterDate(year+'-08-01', year+'-09-01')
.sort('CLOUD_COVER',false);
//Функуция удаления облаков
var fDeleteClouds = function(image) {
var cloud =
ee.Algorithms.Landsat.simpleCloudScore(image);
var score =
cloud.select(['cloud']).lte(20);
return image.updateMask(score);
};
//Удаление облаков из композитов
var composite_2 = collection_2.map(fDeleteClouds)
.mosaic()
.clip(region);
var clouds_2 =
composite_2.select('B3').subtract(free_clouds.select('B3'))
;
var cloud_mask_2 = clouds_2.expression('float(b("B3"))>0.03
? 1:0').clip(region);
var min = composite_2.reduceRegions({
collection: geometry,
reducer: ee.Reducer.min(),
scale: 30,
});
var max = composite_2.reduceRegions({
collection: geometry,
reducer: ee.Reducer.max(),
scale: 30,
});
86
var MinFeature =
ee.Feature(min.first()).select(composite_2.bandNames());
var MinDictionary = MinFeature.toDictionary();
var minB1 = MinDictionary.get('B1');
var minB2 = MinDictionary.get('B2');
var minB3 = MinDictionary.get('B3');
var minB4 = MinDictionary.get('B4');
var minB5 = MinDictionary.get('B5');
var minB7 = MinDictionary.get('B7');
var minB1 = minB1.getInfo();
var minB2 = minB2.getInfo();
var minB3 = minB3.getInfo();
var minB4 = minB4.getInfo();
var minB5 = minB5.getInfo();
var minB7 = minB7.getInfo();
var MaxFeature =
ee.Feature(max.first()).select(composite_2.bandNames());
var MaxDictionary = MaxFeature.toDictionary();
var maxB1 = MaxDictionary.get('B1');
var maxB2 = MaxDictionary.get('B2');
var maxB3 = MaxDictionary.get('B3');
var maxB4 = MaxDictionary.get('B4');
var maxB5 = MaxDictionary.get('B5');
var maxB7 = MaxDictionary.get('B7');
var maxB1 = maxB1.getInfo();
var maxB2 = maxB2.getInfo();
var maxB3 = maxB3.getInfo();
var maxB4 = maxB4.getInfo();
87
var maxB5 = maxB5.getInfo();
var maxB7 = maxB7.getInfo();
print('B1: '+minB1+' - '+maxB1);
print('B2: '+minB2+' - '+maxB2);
print('B3: '+minB3+' - '+maxB3);
print('B4: '+minB4+' - '+maxB4);
print('B5: '+minB5+' - '+maxB5);
print('B7: '+minB7+' - '+maxB7);
var mask =
composite_2.expression('float(b("B1"))>'+minB1+'&&
float(b("B1"))<'+maxB1+
'&& float(b("B2"))>'+minB2+'&& float(b("B2"))<'+maxB2+
'&& float(b("B3"))>'+minB3+'&& float(b("B3"))<'+maxB3+
'&& float(b("B4"))>'+minB4+'&& float(b("B4"))<'+maxB4+
'&& float(b("B5"))>'+minB5+'&& float(b("B5"))<'+maxB5+
'&& float(b("B7"))>'+minB7+'&& float(b("B7"))<'+maxB7+' ?
1:0');
var mask = mask.updateMask(mask.neq(0));
//Создание вектора
var vector = mask.addBands(composite_2).reduceToVectors({
geometry: geometry3,
crs: composite_2.projection(),
scale: 30,
geometryType: 'polygon',
eightConnected: false,
labelProperty: 'zone',
reducer: ee.Reducer.mean(),
maxPixels: 5000000000
88
});
var options = {
title: 'Landsat 5 TM',
fontSize: 14,
hAxis: {title: 'TOA'},
vAxis: {title: 'Pixels'},
series: {
0: {color: 'blue'},
1: {color: 'green'},
2: {color: 'red'},
3: {color: 'magenta'},
4: {color: 'olive'},
5: {color: 'yellow'},
7: {color: 'maroon'}}
};
var histogram = Chart.image.histogram(composite_2,
geometry, 30)
.setSeriesNames(['B1', 'B2', 'B3', 'B4', 'B5', 'B7'])
.setOptions(options);
print(histogram);
Map.addLayer(composite_2, {bands: ['B7', 'B4', 'B2'],min:0,
max: 0.4}, 'collection 2');
Map.addLayer(HotSpots,{palette: 'ff0000'},'HotSpots');
Map.addLayer(mask,{palette: 'red'},'mask');
var burned = ee.Image(0).updateMask(0).paint(vector,
'000000',2);
89
var kvartalka29 =
ee.Image(0).updateMask(0).paint(kvartalka, '000000',2);
Map.addLayer(kvartalka29,{},'kvartalka');
Map.addLayer(MCD45,{},'MCD45');
Map.addLayer(burned,{palette: 'ff0000'},'burned');
Map.addLayer(MCD45,{},'MCD45');
Export.image.toDrive({
image: mask,
description: 'mask',
scale: 30,
region: geometry2,
maxPixels: 5000000000
});
90
Исходный код для метода 2:
//Выбор года для определения территорий пострадавших от
пожаров
var year = 2011;
//Выбираем вариант 1 - Выбор региона России, 2 - выбор
лесничества Архангельской области
//reg - номер региона, NAME_LES - названия лесничества
var variant = 1;
var reg = 29;
var NAME_LES = 'Мезенское';
//Значение NBR менее которого считается гарью
var BurnedValue = 0.3;
//Выбор года для определения территорий пострадавших от
пожаров
if (variant==1){
//Загружаем границы регионов России
var Russia =
ee.FeatureCollection('ft:1z7BRawAQRGo2kznrwVBCD74UuPJIICrvZ
PEbGIfA');
//Фильтруем регионг по номеру региона
var region = Russia.filterMetadata('name', 'equals', reg);
} else if(variant==2){
//Загружаем границы регионов России
var Arh_les = ee.FeatureCollection('ft:1Z55isShK22pFEL0i6O0NeyLJnj6B5ckrDQ6CnNi');
//Фильтруем регионг по номеру региона
var region = Arh_les.filterMetadata('NAME_LES', 'equals',
NAME_LES);
}
91
//KML файл с Hot Spots на Архангельскую область за 20062015 годы
var HotSpotsKML =
ee.FeatureCollection('ft:1i1_y8xKqQqNCZj4jXtw9hNNRbaGymmiU2
x1NejNw');
//Фильтруем Hot spots по году и по выбранному для анализа
региону
var HotSpots = HotSpotsKML.filterMetadata('YEAR', 'equals',
year).filterBounds(region);
//Каналы используемые в обработке снимков
var B1 = 'B1';
var B2 = 'B2';
var B3 = 'B3';
var B4 = 'B4';
var B5 = 'B5';
var B7 = 'B7';
//Условие выбора спутника для анализа
if(year==2012){
var satellite = 'Landsat 7 ETM+';
var sat = 'LE7_L1T_TOA';
} else if(year>=2013){
var satellite = 'Landsat 8 OLI';
var sat = 'LC8_L1T_TOA';
var B1 = 'B2'; //т. к. Landsat 8 отличается по каналам
var B2 = 'B3';
var B3 = 'B4';
var B4 = 'B5';
var B5 = 'B6';
92
var B7 = 'B7';
} else if (year>=2001){
var satellite = 'Landsat 5 TM';
var sat = 'LT5_L1T_TOA';
} else {
var satellite = 'Нет данных FIRMS';
var sat = 'LT5_L1T_TOA';
}
print(satellite);
print(sat);
//Даты для FIRMS период определения гарей
var StarFire = year + '-05-01';
var EndFire = year + '-10-01';
//Даты для создания безоблачного композита
var Start = '2000-7-1';
var End = '2010-9-1';
//Даты для создания до-пожарного композита
var StartUntilFire = (year-1) + '-08-01';
var EndUntilFire = (year-1) + '-09-01';
//Даты для создания после- пожарного композита
var StartAfterFire = year + '-08-01';
var EndAfterFire = year + '-09-01';
/////////////////////////FIRMS/////////////////////////////
///////
//Архив FIRMS
var FIRMS = ee.ImageCollection('FIRMS')
.filterBounds(region)
93
.filterDate(StarFire, EndFire);
var T21 = FIRMS.mosaic()
.clip(region)
.select('T21');
//Создание маски для векторизации
var zones = T21.gt(0);
zones = zones.updateMask(zones.neq(0));
//Создание вектора
var vectors = zones.addBands(T21).reduceToVectors({
geometry: region,
crs: T21.projection(),
scale: 500,
geometryType: 'polygon',
eightConnected: false,
labelProperty: 'zone',
reducer: ee.Reducer.mean()
});
//Функция для создания буфера для каждого полигона
var buffer = function(feature) {
return feature.buffer(3000);
};
//Расширение границ vectors на 1000 м
var bounds = vectors.map(buffer);
///////////////////////MODIS
MCD45////////////////////////////////
var MCD45 = ee.ImageCollection('MODIS/051/MCD45A1')
.filterBounds(geometry)
94
.filterDate(StarFire, EndFire);
///////////////////////Landsat/////////////////////////////
///////
//Создание безоблачного композита
var L5 = ee.ImageCollection('LANDSAT/LT5_L1T');
var free_clouds = ee.Algorithms.Landsat.simpleComposite({
collection: L5.filterDate(Start, End),
asFloat: true});
//Создание до-пожарного композита
var collection_1 = ee.ImageCollection('LANDSAT/'+sat)
.filterBounds(region)
.filterDate(StartUntilFire, EndUntilFire)
.sort('CLOUD_COVER',false);
//Создание после-пожарного композита
var collection_2 = ee.ImageCollection('LANDSAT/'+sat)
.filterBounds(region)
.filterDate(StartAfterFire, EndAfterFire)
.sort('CLOUD_COVER',false);
//Функуция удаления облаков
var fDeleteClouds = function(image) {
var cloud =
ee.Algorithms.Landsat.simpleCloudScore(image);
var score =
cloud.select(['cloud']).lte(20);
return image.updateMask(score);
};
//Удаление облаков из композитов
var composite_1 = collection_1.map(fDeleteClouds)
.mosaic()
95
.clip(bounds);
var composite_2 = collection_2.map(fDeleteClouds)
.mosaic()
.clip(bounds);
//Создание масок облачности и водной маски
var clouds_1 =
composite_1.select(B3).subtract(free_clouds.select(B3));
var cloud_mask_1 =
clouds_1.expression('float(b("'+B3+'"))>0.03 ? 1:0');
var clouds_2 =
composite_2.select(B3).subtract(free_clouds.select(B3));
var cloud_mask_2 =
clouds_2.expression('float(b("'+B3+'"))>0.03 ?
1:0').clip(bounds);
var water_mask =
free_clouds.expression('((float(b("'+B2+'"))+float(b("'+B3+
'")))/(float(b("'+B4+'"))+float(b("'+B5+'"))))>=1?1:0');
//Объединение всех масок в одну
var mask = cloud_mask_1.add(cloud_mask_2).add(water_mask);
var mask =
mask.expression('b("constant")>0?1:0').clip(bounds);
//Вычисляем индексы NBR
var NBR_1 = composite_1.expression('float(b("'+B4+'") b("'+B7+'")) / (b("'+B4+'") + b("'+B7+'"))');
var NBR_2 = composite_2.expression('float(b("'+B4+'") b("'+B7+'")) / (b("'+B4+'") + b("'+B7+'"))');
//Убираем из NBR_1 и NBR_2 облачность
var other_mask_1 =
cloud_mask_1.expression('b("constant")>0?0:1').clip(bounds)
;
var NBR_1 = NBR_1.updateMask(other_mask_1);
96
//var other_mask_2 =
mask.expression('b("constant")>0?0:1').clip(bounds);
var NBR_2 = NBR_2.updateMask(mask.eq(0));
//Поиск изменений
var mask_NBR_1 =
NBR_1.expression('b("B4")<'+BurnedValue+'?1:0').clip(bounds
);
var mask_NBR_2 =
NBR_2.expression('b("B4")<'+BurnedValue+'?1:0').clip(bounds
);
var bmask =
mask_NBR_2.expression('mask_NBR_2>mask_NBR_1?1:0', {
'mask_NBR_1' : mask_NBR_1.select('constant'),
'mask_NBR_2' : mask_NBR_2.select('constant'),
}).clip(bounds);
//Создание маски гарей
var burned_mask = NBR_2.lt(BurnedValue);
burned_mask = burned_mask.updateMask(burned_mask.eq(1));
//Создание вектора по маске
var burned_vectors =
burned_mask.addBands(NBR_2).reduceToVectors({
geometry: geometry,
//geometry-region
crs: NBR_2.projection(),
scale: 30,
geometryType: 'polygon',
eightConnected: false,
labelProperty: 'zone',
reducer: ee.Reducer.mean(),
maxPixels: 5000000000
97
});
//Отображение полученных данных
var region = ee.Image(0).updateMask(0).paint(region,
'000000',2);
Map.addLayer(region,{palette:'#000000'},'Boundary
region',0);
var display = ee.Image(0).updateMask(0).paint(vectors,
'000000',2);
Map.addLayer(display,{palette: '#4B0082'},'FIRMS'+(year),
0);
Map.addLayer(bmask,{},'bmask '+year,0);
Map.addLayer(HotSpots,{color:'ff0000'},'hot spots
'+year,0);
Map.addLayer(free_clouds, {bands: [B3, B2, B1],min:0.0,
max: 0.3, gamma: [0.95, 1.1, 1]}, 'free_clouds',0);
Map.addLayer(composite_1, {bands: [B7, B4, B2],min:0, max:
0.3}, 'collection '+(year-1),0);
Map.addLayer(composite_2, {bands: [B7, B4, B2],min:0, max:
0.3}, 'collection '+year,1);
Map.addLayer(NBR_1,{},'NBR '+(year-1),0);
Map.addLayer(NBR_2,{},'NBR '+(year),0);
//Map.addLayer(mask,{},'mask');cloud_mask_2 !!! Пока не
разобрался с change detection !!!
Map.addLayer(water_mask,{},'Water mask',0);
//Map.addLayer(cloud_mask_2,{},'cloud_mask '+(year),0);
var burned =
ee.Image(0).updateMask(0).paint(burned_vectors,
'000000',2);
Map.addLayer(burned,{palette: 'ff0000'},'burned
'+(year),0);
Map.addLayer(MCD45,{},'MCD '+(year),0);
98
//Экспорт вектора гарей
Export.table.toDrive({
collection: burned_vectors,
description:'burned',
fileFormat: 'KML'
});
//Экспорт слоев NBR
Export.image.toDrive({
image: NBR_1,
description: 'NBR_'+(year-1),
scale: 30,
region: geometry,
maxPixels: 2697651327
});
Export.image.toDrive({
image: NBR_2,
description: 'NBR_'+year,
scale: 30,
region: geometry,
maxPixels: 2697651327
});
99
ПРИЛОЖЕНИЕ Б
1
2
3
4
5
6
Видимый
(красный)
NIR (ближний
инфракрасный)
Видимый
(синий)
Видимый
(зеленый)
NIR (ближний
инфракрасный)
MIR (средний
инфракрасный)
0.841 – 0.876
0.459 – 0.479
0.545 – 0.565
500
1.230 – 1.250
1.628 – 1.652
8
0.405 – 0.420
2300
1-2
дня
12
0.438 – 0.448
(синий)
0.483 – 0.493
10
11
Видимый
0.526 – 0.536
12
(зеленый)
0.546 – 0.556
13 h
13 i
Видимый
14 h
(красный)
14 i
15
16
17
NIR (ближний
инфракрасный)
1000
0.662 – 0.672
0.673 – 0.683
0.743 – 0.753
0.862 – 0.877
2300
0.890 – 0.920
100
1-2
дня
12
1000
(м.)
разрешение
венное
Пространст
(бит)
разрешение
ическое
Радиометр
съемки
(км.)
250
2.105 – 2.155
Видимый
Период
обзора
0.620 – 0.670
7
9
полосы
Ширина
(мкм.)
диапазон
ый
Спектральн
ла
Номеркана
Таблица 1 - Основные технические характеристики MODIS
18
0.931 – 0.941
19
0.915 – 0.965
20
3.660 – 3.840
21
3.929 – 3.989
22
3.929 – 3.989
23
4.020 – 4.080
24
4.433 – 4.498
25
4.482 – 4.549
26
1.360 – 1.390
27
6.535 – 6.895
28
7.175 – 7.475
29
8.400 – 8.700
30
9.580 – 9.880
TIR (тепловой
31
32
33
34
35
36
инфракрасный)
10.780–
11.280
11.770–
12.270
13.185–
13.485
13.485–
13.785
13.785–
14.085
14.085–
14.385
101
(м.)
разрешение
венное
Пространст
(бит)
разрешение
ическое
Радиометр
съемки
(км.)
Период
обзора
полосы
Ширина
(мкм.)
диапазон
ый
Спектральн
ла
Номеркана
Продолжение таблицы 1
Таблица 2 - Основные характеристики сканера TM (Landsat 5)
Номера
каналов
1
2
3
4
5
6
7
Спектральный диапазон Ширина
полосы
обзора
(км)
Видимый
0.45 (синий)
0.515
Видимый
0.525 (зеленый)
0.605
Видимый
0.63 (красный)
0.690
NIR (ближний 0.75 185
инфракрасный) 0.90
NIR (ближний 1.55 инфракрасный) 1.75
TIR (тепловой 10.40 инфракрасный) 12.5
MWIR
2.09 (дальний ИК)
2.35
Период
съемки
РадиометПространственное
рическое
разрешение
разрешение м
(бит)
30
16 дней
8
120
30
Таблица 3 - Характеристики аппаратуры спутника Sentinel-1
Спектральный диапазон
С-диапазон
Периодичность съемки,
1-3
сутки
Режим
Номинальное
Ширина полосы
пространственное
съемки, км
Поляризация
разрешение, м
Stripmap (SM; single-look)
5x5
80
Interferometric Wide Swath
5x20
240
(IWS; single-look)
Extra Wide Swath (EWS;
400
single-look)
Interferometric Wide Swath
25x80
240
20х5
20х20
(IWS; 3 looks)
Wave mode (WM; singlelook)
102
Двойная (по
выбору —
HH/HV или
VV/VH)
Таблица 4 - Характеристики аппаратуры спутника Sentinel-2
Центр диапазона канала,
µm
Разрешение, м
Band 1 - Coastal aerosol
0.443
60
Band 2 - Blue
0.490
10
Band 3 - Green
0.560
10
Band 4 - Red
0.665
10
Band 5 - Vegetation Red Edge
0.705
20
Band 6 - Vegetation Red Edge
0.740
20
Band 7 - Vegetation Red Edge
0.783
20
Band 8 - NIR
0.842
10
Band 8A - Vegetation Red Edge
0.865
20
Band 9 - Water vapour
0.945
60
Band 10 - SWIR - Cirrus
1.375
60
Band 11 - SWIR
1.610
20
Band 12 - SWIR
2.190
20
Band 2 - Blue
0.490
10
Каналы
103
Отзывы:
Авторизуйтесь, чтобы оставить отзыв