На сайте Directionsmag.com размещена статья, в которой Лука Конгедо (Luca Congedo) представляет свое пособие для выполнения классификаций сельхозугодий. В частности, он использует снимки Landsat на территорию штата Канзас, неподалеку от города Улисс, используя новую версию плагина (2.3.2) полуавтоматический классификации в QGIS.

В данной статье представлено пособие для выполнения классификаций сельхозугодий. В частности, задачей является классификация снимков Landsat на территорию американского штата Канзас, недалеко от города Улисс, используя новую версию плагина (2.3.2) Полуавтоматический классификации в QGIS.

Перед тем как начать изучать пособие, пожалуйста, посмотрите следующее видео по теме, которое очерчивает область исследования и обеспечивает очень полезную информацию для интерпретации изображений в ненатуральной цветопередаче (видео предоставлено Европейским космическим агентством/ESA). Кроме того, с кратким описанием того, что мы собираемся классифицировать, можно ознакомиться на сайте ESA.

Как показано в предыдущем видео, территория исследования в основном покрыта пахотными землями, и мы собираемся классифицировать снимки Landsat 8, полученные в апреле 2013 года (доступны на сайте Геологической службы США - USGS), используя полуавтоматический подход (то есть все изображение классифицируется автоматически на основе нескольких образцов "обучающих (ROI)" пикселей, которые используются в качестве эталонов).

Перед тем как начать процесс классификации, мы собираемся преобразовать значения исходной яркости (DN) в значения коэффициента отражения (reflectance). Коэффициент отражения определяется как соотношение между отраженной и падающей энергией на поверхность (здесь доступна дополнительная информация о преобразовании исходных значений в коэффициенты отражения для снимков Landsat).

Атмосферная коррекция значений коэффициента отражения полезна во время предварительного этапа обработки, потому что электромагнитная энергия, которая измеряется со спутника, зависит от атмосферных эффектов (например поглощения или рассеяния). Таким образом, преобразование изображений в значtния отражательной яркости позволяет сравнивать в том числе и мозаики различных спутниковых снимков (например, Landsat 5 и Landsat 8), и улучшает результаты классификации.

Плагин полуавтоматической классификации предоставляет возможности для простой коррекции с помощью метода DOS1 (краткое объяснение о преобразования изображений в отражение и атмосферная коррекция DOS1 проиллюстрированы в желтой рамке в конце этого поста). Стоит отметить, что GRASS GIS также имеет инструмент для расчета отражения на верхней границе атмосферы (см. i.landsat.toar) и других DOS методов.

Во-первых, нужно загрузить снимки. Изображение представляет собой подмножество всех сцен, в том числе каналов Landsat (каждый канал представлен в виде 16-битного растра):
• 2 - Синий;
• 3 - Зеленый;
• 4 - Красный;
• 5 - Ближний инфракрасный;
• 6 - Коротковолновый инфракрасный 1;
• 7 - Коротковолновый инфракрасный 2.

Следующий алгоритм иллюстрирует основные этапы полуавтоматической классификации (см. рис. 1)

1. Преобразование растровых каналов из DN в значения отражательной яркости

Этот этап является частью предварительной обработки изображения, которая улучшает дальнейшие возможности использования снимков за счет преобразования значений DN (то есть исходной яркости) в физические значения отражения на верхней границе атмосферы (TOA). Кроме того, мы собираемся применить простую атмосферную коррекцию с помощью метода DOS1 (вычитание темного объекта 1).
В качестве примера приведем следующие изображения - на них отображены три различных спектральных диаграммы объектов (воды, растительности, и голой почвы), рассчитанные по исходным изображением (DN), коэффициенту отражения TOA и исправленным алгоритмом DOS1 изображениям. Как вы можете видеть, коррекция DOS1 изменяет значения отражательной способности, особенно сильно в синем канале. Обратите внимание, что номер канала на изображениях показан в порядковой нумерации для конкретного снимка, так, например: канал 1 на диаграмме является каналом 2 на Landsat 8 (то есть это синий канал).

Рис. 2: График для исходного изображения

Рис. 3: График по коэффициенту отражения TOA

Рис. 4: График по исправленному алгоритмом DOS1 изображению

- Откройте QGIS и запустите плагин полуавтоматической классификации, в главном меню выберите вкладку Pre processing > Landsat;
- Выберите каталог, содержащий каналы снимка Landsat (а также метафайл MTL.txt), и выберите выходной каталог, куда будут сохраняться преобразованные каналы;
- Проверьте, установлен ли флажок напротив знаечния "Применить DOS1 атмосферную коррекцию", и нажмите "Выполнить преобразование" для преобразования каналов Landsat в коэффициенты отражения;
После завершения процесса, преобразованные каналы будут загружены в QGIS (см. рис. 5).

2. Определение входящих значений для классификации

Сейчас мы создадим цветной композит изображения. А именно, композит RGB = 543 (что является эквивалентом RGB = 432 для Landsat 7) полезен для интерпретации пахотных земель (как показано в предыдущем видео ЕКА), вследствие того, что здоровая растительность отражает большую часть падающего света в ближней инфракрасной зоне.
Для того, чтобы создать цветовой композит в QGIS, мы создаем виртуальный растр. Затем определяем канал, который мы собираемся, чтобы классифицировать. Наконец, мы создаем тренировочный шейп-файл, в котором будут храниться результаты трансформирования (т.е. участки, представляющие для нас интерес), определяющие классы растительного покрова.

В меню Растр необходимо выбрать Разное> Создать виртуальный растр (каталог); нажмите кнопку Выбрать... и выберите каналы 3, 4, и 5; определите выходной файл (например rgb.vrt); отметьте галочкой Отдельно (каналы будут сохранены по отдельности) и нажмите OK (см. рис. 6);

В плагине dock ROI нажмите на кнопку Канал, расположенную рядом с кнопкой Выберите изображение, тогда откроется вкладка с группой каналов; нажмите кнопку Выбрать все, а затем Добавить установочные растры (выберите каналы по названиям, отсортированным по алфавиту в порядке возрастания, используя ↑ и ↓ стрелки) — см. рис. 7;

Для того чтобы создать учебный шейп (полигональный слой, имеющий несколько полей атрибутов, которые будут хранить значения ROI), в меню dock ROI нажмите кнопку Создать шейп и выберите, куда сохранить его (например ROI.shp) — см. рис. 8.

3. Создание ROI

Нам нужно создать несколько ROI, учитывая спектральную изменчивость классов земного покрова. Важно, чтобы эти эталоны представляли из себя однородные участки изображения, поэтому мы будем создавать их путем наращивания территорий вокруг пиксела (то есть сегментация изображения, за счет группировки однородных пикселей).

Ниже приведены классы почвенно-растительного покрова, которые мы собираемся определять на снимке:
- сельхозозяйственные земли (например, поля с зеленой растительностью);
- низкорослая растительность (например, поля без зеленой растительности или кустарники);
- застроенные территории (например, искусственные участки, здания и асфальт);
- фермы;
- голая почва;
- вода.

Как вы можете видеть на композитном изображении (см. рис. 9), цвет пахотных земель может варьироваться от темно-красного до светло-красного цвета. Мы должны создать несколько эталонов для каждого из этих цветов.

Для создания ROI нажмите кнопку "+", расположенную рядом с кнопкой "Создать ROI", а затем нажмите на любой пиксель изображения; приблизьтесь к этому участку карты и нажмите на красный пиксель в круге; через несколько секунд появится ROI многоугольник (полупрозрачный оранжевый многоугольник);

В поле определение для ROI введите краткое описание (например, сельскохозяйственное поле; заметьте, что это описание не будет использоваться в процессе классификации, но полезно для различения последующего использования ROI); поле ID (то есть идентификатор) используется в качестве источника для классификации почвенно-растительного покрова, поэтому важно, чтобы каждая категория имела уникальное значение ID (ROI, имеющие одинаковый идентификатор оцениваются как единый ROI). Так как мы должны будем набрать более одного ROI в каждом классе, то установите ID=11, первая цифра будет указывать класс, а вторая номер ROI (другие вариации одной и той же культуры будут иметь значения от 11 до 19, таким же образом, мы устанавливаем ID от 21 до 29 для низкорослой растительности и так далее);

Для того, чтобы сохранить ROI нажмите кнопку Сохранить ROI в шейп-файл;
В главном интерфейсе плагина выберите вкладку инструменты ROI > Спектральный образ, которая отображает наборы набранных ROI, и выберите пункт сельхоз, если поставить плажок напротив σ, то будет показано стандартное отклонение каждого ROI (см. рис. 10).
Повторите эти действия для каждого класса (для того, чтобы получить необходимый ROI, измените значения Минимального размера ROI и радиуса поиска); после сбора нескольких ROI, полезно визуализировать их спектральные характеристики для того, чтобы избежать набора спектрально похожих элементов. Ниже приведено несколько примеров спектральных образов для различных классов почвенно-растительного покрова (рис. 11-15).
Пользователи могут скачать итоговый шейп ROI (в котором автор сохранил 14 образцов).

4. Классификация территории исследования

При выполнении классификации бывает полезно выполнять ее предварительный просмотр, чтобы оценить качество (например, общей проблемой является то, что открытая почва классифицируется как застроенные территории и наоборот). Если результат предварительного просмотра не очень хорош, то мы должны вернуться к фазе создания ROI, удалить некачественные и/или создать новые.
Классификация может быть выполнена с использованием одного из нескольких алгоритмов, в этот раз мы собираемся использовать классификацию Спектральных углов.
В меню классификации, в графе Предварительный просмотр классификации установите размер = 500 (это сторона изображения предварительного просмотра классификации в пикселях); нажмите кнопку +, а затем нажмите на изображение; через несколько секунд предварительный просмотр классификации будет загружен (с цветами по умолчанию);
Для лучшего понимания результатов интерпретации классов, полезно назначить цвета в изображении классификации для каждого класса.
Кроме того, мы можем сохранить стиль визуализации в QGIS, и затем в любой момент воспользоваться этим файлом в качестве стиля по умолчанию для каждой классификации или ее предварительного просмотра.
В панели Слои (см. рис. 16), щелкните левой кнопкой мыши на классификации предварительного просмотра и выберите Свойства > Стиль; выберите тип Псевдоцвета по одному каналу и измените цвета и подписи классов, в соответствии с учебным шейп-файлом, а затем нажмите кнопку Сохранить стиль, чтобы сохранить изменения. Файл назовите, например classification.qml;
В графе стиль классификации нажмите на кнопку Выбрать QML, выберите файл classification.qml, следующая классификация или ее пред.просмотр будут загружены уже в соответствии с этим стилем;
Когда мы будем удовлетворены результатом предварительного просмотра, можно будет выполнить окончательную классификацию (выходным файлом классификации является растровый файл в TIF-формате).
Нажмите кнопку Выполнить классификацию и выберите, куда сохранить выходной файл (например classification.tif) — см. рис. 17.
На сайте можно скачать окончательные результаты классификации почвенно-растительного покрова.

5. Расчет точности классификации

Это полезно (и необходимо), чтобы оценить точность классификации, для того чтобы понять ее качество и выявить ошибки на карте.
Пожалуйста, обратите внимание, что эта классификация выполнена только для демонстрационных целей, так как точная классификация почвенно-растительного покрова потребует дополнительных данных и полевого обследования.
Однако для того, чтобы оценить точность классификации, мы можем сравнить итоговую классификацию земельного покрова с теми ROI, которые мы создали (теоретически, пиксели изображения, принадлежащие к определенному ROI должны быть классифицированы как соответствующий ему класс).

Выберите вкладку Постобработка>Точность в главном интерфейсе плагина (см. рис. 18); выберите classification.tif, затем выберите классификацию для оценки и выберите шейп ROI, нажмите кнопку Рассчитать матрицу ошибок. После этого матрица будет отображена, вы сможете сохранить ее, нажав на кнопку сохранить матрицу ошибок в файл;

Как правило, точность классификации считается хорошей, если значения этой матрицы больше 80%. В этой классификации, как вы можете увидеть, общая точность классификации составляет около 91,4% (матрица сообщений об ошибках приведена ниже). Кроме того, полезно рассмотреть ошибки для одиноких классов (то есть количество недобранных пикселей и взятых излишне - омиссия и комиссия) — см. рис. 19.

6. Расчет статистики классификации

Классы итоговой классификации содержат больше основных классов почвенно-растительного покрова, чем мы определили на шаге 3 (это произошло потому, что мы назначили различные идентификаторы для ROI, принадлежащих к одному и тому же классу).

Для того чтобы оценить сельхозяйственные угодья в этом районе, мы перекодируем растр классификации в меньшее количество классов, в соответствии с основными классами земельного покрова.
Из Processing Toolbox перейдите к SAGA > Инструменты сетки > переклассифицировать значения реклассификации сетки; в окне инструмента, во вкладке Сетка selectclassification.tif; выберите [2] простую таблицу — см. рис. 20;

Под таблицей нажмите на кнопку ..., чтобы открыть окно полной таблицы; справочная таблица имеет 3 колонки: минимум, максимум, и новое. Столбцы минимум и максимум определяют диапазон переклассификационных значений, в соответствии с приведенной ниже таблицей (см. таблицу 1).
Нажмите "OK" и после через несколько секунд будут загружены результаты переклассифицированных значений.
Теперь мы можем создать отчет о классификации, который покажет проценты и площадь классов земельного покрова.
Выберите вкладку Запустить обработку > Отчет о классификации в главном интерфейсе плагина; выберите Переклассифицировать, затем классификация и нажмите Рассчитать отчет о классификации; через несколько секунд будет выведен отчет (см. таблицу 2 ниже).

Эти результаты показывают, что в исследуемом районе (в момент получения изображений) насчитывается около 738 кв. км сельскохозяйственных угодий (класс 1, около 10% от исследуемой области) и 6231 кв. км низкорослой растительности (класс 2, около 87%). Кроме того, можно заметить, что около 14 кв. км занимают фермы (класс 4) и около 18% застроено (класс 3).
Конечно, эти цифры являются лишь демонстрационными результатами; для более точной оценки территории (и, возможно идентификации типов культур), необходимы данные полевых обследований, которое помогут улучшить качество набора ROI и правильнее оценивать точность классификации.
Однако целью этого урока является намерение показать, как проводить классификацию почвенно-растительного покрова и то, как это может быть выполнено просто и быстро с помощью программ с открытым исходным кодом, которые полезны сразу для нескольких видов деятельности, таких как мониторинг растительного покрова или точное земледелие.
Ниже приведено видео этого урока.

Преобразование изображений Landsat в значения отражения и атмосферная коррекция DOS1

Исходное отражение (radiance) является "потоком энергии, излучаемым единицей земной поверхности в заданном направлении", "Исходное отражение - это то, что измеряется на датчике и в некоторой степени зависит от коэффициента отражения" (НАСА, 2011, р. 47).
Спектральное отражение на датчике сенсора (Lλ) измеряется в [Вт/(кв.м*стерридиан*мкм)] и для снимков Landsat он определяется следущим документом: (https://landsat.usgs.gov/Landsat8_Using_Product.php):
Lλ = ML * Qcal + AL
где:
- ML = зависящий от канала коэффициент пересчета, для снимков Landsat определен в метаданных (RADIANCE_MULT_BAND_x, где х номер зоны)
- AL = зависящий от канала добавочный коэффициент пересчета, для снимков Landsat определен в метаданных (RADIANCE_ADD_BAND_x, где х номер зоны)
- Qcal = Квантованные и калиброванные стандартные значения пикселей (DN)
"Для относительно четких сцен со спутника Landsat, сокращение различий между сценами может быть достигнуто путем нормализации солнечного излучения путем преобразования спектральных значений, в значения планетарного отражения или альбедо. Это объединяет поверхностное и атмосферное отражение Земли и вычисляется по следующей формуле" (НАСА, 2011, стр. 119):
ρp = (π * Lλ * d^2)/ (ESUNλ *cosθs)
где:
- pp = безразмерное значение планетарного отражения, это "отношение того, как отражает Земля по сравнению с общим значением отражения энергии" (НАСА, 2011, стр. 47)
- Lλ = Спектральное отражение на датчике (отражение на спутнике)
- d = расстояние от Земли до Солнца в астрономических единицах (поставляется с метаданными Landsat 8, и также можно получить из файла Excel по ссылке http://landsathandbook.gsfc.nasa.gov/excel_docs/d.xls)
- ESUNλ = среднее солнечное экзо-атмосферное излучение,
- Θs = зенитный угол Солнца в градусах, который равен θs = 90 ° - θe, где θe - высота Солнца
Стоит отметить, что снимки Landsat 8 сопровождаются коэффициентами пересчета для каждого из каналов, которые позволяют производить прямой пересчет исходныз значений DN значения отражения TOA. Тем не менее, влияние атмосферы (т.е. изменение отражательной способности, которое изменяется с длиной волны) следует рассматривать в аспекте измерения значений отражения на поверхности земли. Согласно исследованиям Морана и др.. (1992), отражательное значение земли (ρ) это:
ρ = [π * (Lλ - Lp) * d^2]/ [Tv * ((ESUNλ * cosθs * Tz) + Edown)]
где:
- Lp, отражение расстояния
- Tv, прозрачность атмосферы в направлении съемки
- Tz, прозрачность атмосферы в направлении освещения
- Edown, нисходящая рассеяная освещенность
Поэтому нам необходимо несколько атмосферных измерений для расчета ρ (исправлений на основе физики атмосферы). В качестве альтернативы, можно использовать метод расчета этих параметров для изображения, без натурных измерений во время получения изображения.
Метод вычитания темного объекта (DOS) - это группа атмосферных поправок для космических снимков.
Чавес (1996) объясняет, что «основное предположение состоит в том, что в пределах изображения некоторые пиксели находятся в полной тени и их излучение, полученное на спутнике, обусловлено атмосферным рассеянием (рассеиванием расстояния). Это предположение рассматривается в сочетании с тем, что очень немногие цели на поверхности Земли являются абсолютно черными, так, предполагается, что один процент минимального отражения лучше, чем ноль процентов". Стоит отметить, что точность этого метода, как правило, ниже, чем поправки, сделанные по данным физических наблюдений, но они бывают очень полезны, когда никакие атмосферные измерения не доступны, поскольку они могут улучшить оценку поверхностного отражения.
Отражение на расстоянии дается согласно работе (Sobrino и др., 2004):
Lp = Lmin – LDO1%
где:
- Lmin = "отражение, которое соответствует цифровым значениям камеры, для которых сумма всех пикселей с такими значениями меньше или равна 0,01% всех пикселей снимка" (Sobrino и др. др.., 2004, с. 437), поэтому этому значению присваивается исходное цифровое значение счетчика (DNmin)
- LDO1% = отражение темного объекта, предполагается, что оно имеет значение коэффициента отражения 0,01
Для снимков Landsat:
Lmin = ML * DNmin + AL
The radiance of Dark Object is given by (Sobrino, et al., 2004):
LDO1% = 0,01* [(ESUNλ * cosθs * Tz) + Edown] * Tv / (π * d^2)
Здесь отражение расстояния это:
Lp = ML* DNmin + AL – 0,01* [(ESUNλ * cosθs * Tz) + Edown] * Tv / (π * d^2)

Lp = ML* DNmin + AL – 0,01* ESUNλ * cosθs / (π * d^2)
Есть несколько методов DOS (например DOS1, DOS2, DOS3, DOS4), основанные на различных исходных предположениях о Tv, TZ, и Edown.
Простейший метод - это DOS1, где сделаны следующие допущения (Моран и др., 1992):
- ТВ = 1
- Tz = 1
- Edown = 0
Поэтому отражение расстояния:
Lp = ML * DNmin + AL - 0,01 * ESUNλ * cosθs / (π * D ^ 2)
И в результате отражение на земной поверхности определяется по формуле:
ρ = [π * (Lλ - Lp) * д ^ 2] / (ESUNλ * cosθs)
Esun [Вт / (м2 * мкм)] значения для датчиков Landsat приведены в следующей таблице (таблица 3).
 * согласно Chander & Markham (2003)
** согласно Finn и др.(2012)

Для Landsat 8, Esun может быть рассчитана, как описано в документе: http://grass.osgeo.org/grass65/manuals/i.landsat.toar.html
Esun = (π * d2) * RADIANCE_MAXIMUM / REFLECTANCE_MAXIMUM
где RADIANCE_MAXIMUM и REFLECTANCE_MAXIMUM предоставляются с метаданными снимков.