Нейросеть для классификации спутниковых снимков с помощью Tensorflow на Python

Автор оригинала: Pratyush Tripathy
  • Перевод
  • Tutorial


Это пошаговая инструкция по классификации мультиспектральных снимков со спутника Landsat 5. Сегодня в ряде сфер глубокое обучение доминирует как инструмент для решения сложных проблем, в том числе геопространственных. Надеюсь, вы знакомы с датасетами спутниковых снимков, в частности, Landsat 5 TM. Если вы немного разбираетесь в работе алгоритмов машинного обучения, то это поможет вам быстро освоить это руководство. А для тех, кто не разбирается, будет достаточным знать, что, по сути, машинное обучение заключается в установлении взаимосвязей между несколькими характеристиками (набором признаков Х) объекта с другим его свойством (значением или меткой, — целевой переменной Y). Мы подаём на вход модели много объектов, для которых известны признаки и значение целевого показателя/класса объекта (размеченные данные) и обучаем ее так, чтобы она могла спрогнозировать значение целевой переменной Y для новых данных (неразмеченных).

В чём заключается основная проблема со спутниковыми снимками?

Два и более класса объектов (например, застройки, пустыри и котлованы) на спутниковых снимках могут иметь одинаковые спектральные характеристикизначения, поэтому в последние лет двадцать их классификация представляет собой трудную задачу.

Из-за этого, возможно использовать классические модели машинного обучения с учителем и без, но их качество будет далеко от идеала. Они всегда обладают одними и теми же недостатками. Рассмотрим пример:



Если в качестве классификатора использовать вертикальную линию и двигать её вдоль оси Х, то классифицировать изображения домов будет непросто. Данные распределены так, что разделить их на классы с помощью одной вертикальной линии невозможно (в таких случаях говорят, что «объекты разных классов не являются линейно разделимыми»). Но это не означает, что дома вообще нельзя классифицировать!



Давайте с помощью красной линии разделим два класса. В данном случае классификатор определил большинство домов, однако один дом не был отнесен к своему классу, а ещё три дерева были отнесены к «домам» ошибочно. Чтобы не пропустить ни один дом, можно использовать классификатор в виде синей линии. Тогда будут охвачены все дома, то есть мы говорим, что у метрики recall (полнота) высокое значение. Однако не все классифицированные значения оказались домами, то есть мы одновременно получили низкое значение метрики precision (точность). Если же мы воспользуемся зелёной линией, то все изображения, классифицированные как дома, действительно будут домами, то есть классификатор покажет высокую точность. В этом случае полнота будет меньше, поскольку три дома окажутся неучтёнными. В большинстве случаев нам приходится находить компромисс между точностью и полнотой.

Эта проблема домов и деревьев аналогична проблеме с застройкой, пустырями и котлованами. Приоритет метрик классификации спутниковых снимков может меняться в зависимости от задачи. Например, если вам нужно удостовериться, что все застроенные территории без исключений классифицированы как застройка, и вы готовы смириться с наличием пикселей других классов с аналогичными сигнатурами, которые тоже будут классифицированы как застройка, тогда вам потребуется модель с высокой полнотой. А если вам важнее классифицировать именно застройку, без добавления пикселей других классов, и вы готовы отказаться от классификации смешанных территорий, тогда выбирайте классификатор с высокой точностью. В случае с домами и деревьями обычная модель будет использовать красную линию, выдерживая баланс между точностью и полнотой.

Используемые данные


В качестве признаков мы будем использовать значения шести диапазонов (band 2 — band 7) изображения из Landsat 5 TM, и попытаемся спрогнозировать двоичный класс застройки. Для обучения и тестирования будут использованы многоспектральные данные (изображения и слой с бинарным классом застройки) с Landsat 5 за 2011 года для Бангалора. А для прогнозирования будут использованы многоспектральные данные Landsat 5, полученные в 2005-м в Хайдарабаде.
Поскольку мы используем для обучения размеченные данные, это называется обучение с учителем.



Многоспектральные обучающие данные и соответствующий двоичный слой с застройкой.

Для создания нейросети воспользуемся Python—библиотекой Google Tensorflow. Также нам потребуются эти библиотеки:

  1. pyrsgis — для чтения и записи GeoTIFF.
  2. scikit-learn — для предобработки данных и оценки точности.
  3. numpy — для базовых операций с массивами.

А теперь без дальнейших проволочек давайте писать код.

Положите все три файла в директорию, пропишите в скрипте путь и названия входных файлов, а затем прочитайте файлы GeoTIFF.

import os
from pyrsgis import raster

os.chdir("E:\\yourDirectoryName")
mxBangalore = 'l5_Bangalore2011_raw.tif'
builtupBangalore = 'l5_Bangalore2011_builtup.tif'
mxHyderabad = 'l5_Hyderabad2011_raw.tif'

# Read the rasters as array
ds1, featuresBangalore = raster.read(mxBangalore, bands='all')
ds2, labelBangalore = raster.read(builtupBangalore, bands=1)
ds3, featuresHyderabad = raster.read(mxHyderabad, bands='all')

Модуль raster из пакета pyrsgis считывает геолокационные данные GeoTIFF и значения цифровых номеров (DN) в виде отдельных NumPy-массивов. Если вас интересуют подробности, почитайте здесь.

Теперь выведем на экран размер прочитанных данных.

print("Bangalore multispectral image shape: ", featuresBangalore.shape)
print("Bangalore binary built-up image shape: ", labelBangalore.shape)
print("Hyderabad multispectral image shape: ", featuresHyderabad.shape)

Результат:

Bangalore multispectral image shape: 6, 2054, 2044
Bangalore binary built-up image shape: 2054, 2044
Hyderabad multispectral image shape: 6, 1318, 1056

Как видите, в изображениях Бангалора такое же количество строк и столбцов, как и в бинарном слое (соответствующем застройке). Количество слоёв в многоспектральных изображениях в Бангалоре и Хайдарабаде тоже совпадает. Модель будет учиться решать, какие пиксели относятся к застройке, а какие нет, на основе соответствующих значений по всем 6-ти спектрам. Поэтому многоспектральные изображения должны иметь одинаковое количество признаков (диапазонов), перечисленных в одном и том же порядке.

Теперь превратим массивы в двумерные, где каждая строка представляет отдельный пиксель, потому что это нужно для работы большинства алгоритмов машинного обучения. Сделаем мы это с помощью модуля convert из пакета pyrsgis.


Схема реструктуризации данных.

from pyrsgis.convert import changeDimension

featuresBangalore = changeDimension(featuresBangalore)
labelBangalore = changeDimension (labelBangalore)
featuresHyderabad = changeDimension(featuresHyderabad)
nBands = featuresBangalore.shape[1]
labelBangalore = (labelBangalore == 1).astype(int)

print("Bangalore multispectral image shape: ", featuresBangalore.shape)
print("Bangalore binary built-up image shape: ", labelBangalore.shape)
print("Hyderabad multispectral image shape: ", featuresHyderabad.shape)

Результат:

Bangalore multispectral image shape: 4198376, 6
Bangalore binary built-up image shape: 4198376
Hyderabad multispectral image shape: 1391808, 6

В седьмой строке мы извлекли все пиксели со значением 1. Это помогает избежать проблем с пикселями без информации (NoData), которые часто имеют крайне высокие или низкие значения.
Теперь разделим данные на обучающую и валидационную выборки. Это нужно для того, чтобы модель не видела тестовых данных и работала так же хорошо с новой информацией. Иначе модель переобучится и будет хорошо работать только на обучающих данных.

from sklearn.model_selection import train_test_split

xTrain, xTest, yTrain, yTest = train_test_split(featuresBangalore, labelBangalore, test_size=0.4, random_state=42)

print(xTrain.shape)
print(yTrain.shape)

print(xTest.shape)
print(yTest.shape)

Результат:

(2519025, 6)
(2519025,)
(1679351, 6)
(1679351,)
test_size=0.4

означает, что данные разделены на обучающие и валидационные в соотношении 60/40.
Многим алгоритмам машинного обучения, включая нейросети, нужны нормализованные данные. Это означает, что они должны быть распределены в рамках заданного диапазона (в данном случае от 0 до 1). Поэтому, чтобы выполнить это требование, мы нормализуем признаки. Сделать это можно с помощью извлечения минимального значения с последующим делением на разброс (разницу между максимальным и минимальным значениями). Поскольку датасет Landsat является восьмибитным, минимальное и максимальное значение будут равны 0 и 255 (2 = 256 значений).

Обратите внимание, что для нормализации всегда лучше вычислять минимальные и максимальные значения на основе данных. Для упрощения задачи мы будем придерживаться восьмибитного диапазона по умолчанию.

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



# Normalise the data
xTrain = xTrain / 255.0
xTest = xTest / 255.0
featuresHyderabad = featuresHyderabad / 255.0

# Reshape the data
xTrain = xTrain.reshape((xTrain.shape[0], 1, xTrain.shape[1]))
xTest = xTest.reshape((xTest.shape[0], 1, xTest.shape[1]))
featuresHyderabad = featuresHyderabad.reshape((featuresHyderabad.shape[0], 1, featuresHyderabad.shape[1]))

# Print the shape of reshaped data
print(xTrain.shape, xTest.shape, featuresHyderabad.shape)

Результат:

(2519025, 1, 6) (1679351, 1, 6) (1391808, 1, 6)

Всё готово, давайте соберём нашу модель с помощью keras. Для начала воспользуемся sequential-моделью, добавляя слои один за другим. У нас будет один входной слой с количеством узлов, равным числу диапазонов (nBands) — в нашем случае их 6. Также будем использовать один скрытый слой с 14 узлами и функцией активации ReLu. Последний слой состоит из двух узлов для определения двоичного класса застройки с функцией активации softmax, которая подходит для вывода категоризированного результата. Подробнее о функциях активации написано здесь.

from tensorflow import keras

# Define the parameters of the model
model = keras.Sequential([
    keras.layers.Flatten(input_shape=(1, nBands)),
    keras.layers.Dense(14, activation='relu'),
    keras.layers.Dense(2, activation='softmax')])

# Define the accuracy metrics and parameters
model.compile(optimizer="adam", loss="sparse_categorical_crossentropy", metrics=["accuracy"])

# Run the model
model.fit(xTrain, yTrain, epochs=2)


Архитектура нейросети

Как упоминалось в строке 10, мы задаем adam в качестве оптимизатора модели (есть и несколько других). В данном случае будем использовать в качестве функции потерь перекрестную энтропию (en. categorical-sparse-crossentropy — подробнее о ней написано здесь). Для оценки качества работы модели будем использовать метрику accuracy.

Наконец, запустим обучение нашей модели в течение двух эпох (или итераций) на xTrain и yTrain. Это займёт какое-то время, в зависимости от размера данных и вычислительной мощности. Вот что вы увидите после компиляции:



Давайте спрогнозируем значения для валидационных данных, которые мы храним отдельно, и рассчитаем различные метрики точности.

from sklearn.metrics import confusion_matrix, precision_score, recall_score

# Predict for test data 
yTestPredicted = model.predict(xTest)
yTestPredicted = yTestPredicted[:,1]

# Calculate and display the error metrics
yTestPredicted = (yTestPredicted>0.5).astype(int)
cMatrix = confusion_matrix(yTest, yTestPredicted)
pScore = precision_score(yTest, yTestPredicted)
rScore = recall_score(yTest, yTestPredicted)

print("Confusion matrix: for 14 nodes\n", cMatrix)
print("\nP-Score: %.3f, R-Score: %.3f" % (pScore, rScore))

Функция softmax генерирует отдельные столбцы для значений вероятности для каждого класса. Мы используем только значения для первого класса («есть застройка»), как видно по шестой строке приведенного выше кода. Оценивать работу моделей геопространственного анализа не так просто, в отличие от других классических задач машинного обучения. Будет нечестно полагаться на обобщенную суммарную ошибку. Ключом к успешной модели является пространственное расположение. Таким образом, матрица ошибок (confusion matrix), точность и полнота могут дать более корректное представление о качестве работы модели.


Так в консоли отображаются матрица ошибок, точность и полнота.

Как видно из confusion matrix, есть тысячи пикселей, которые относятся к застройке, но классифицированы иначе, и наоборот. Однако их доля от общего объёма данных не слишком велика. Точность и полнота на тестовых данных превысили порог в 0,8.

Вы можете потратить больше времени и выполнить несколько итераций для поиска оптимального количества скрытых слоев, количества узлов в каждом скрытом слое, а также количества эпох для достижения желаемой точности. По мере необходимости, в качестве признаков можно использовать индексы дистанционного зондирования вроде NDBI или NDWI. При достижении желаемой точности используйте модель для прогнозирования застройки на основе новых данных и экспортируйте результат в GeoTIFF. Для подобных задач можно использовать аналогичную модель с небольшими изменениями.

predicted = model.predict(feature2005)
predicted = predicted[:,1]

#Export raster
prediction = np.reshape(predicted, (ds.RasterYSize, ds.RasterXSize))
outFile = 'Hyderabad_2011_BuiltupNN_predicted.tif'
raster.export(prediction, ds3, filename=outFile, dtype='float')

Обратите внимание, что мы экспортируем GeoTIFF со спрогнозированными значениями вероятности, а не с их бинаризированной по порогу версией. Позднее в ГИС-среде мы можем задать пороговое значение слоя типа float, как показано на рисунке ниже.


Слой застройки в Хайдарабаде, спрогнозированный моделью на основе многоспектральных данных.

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

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

Использованные данные и весь код лежат тут.
  • +24
  • 4,9k
  • 8
Инфосистемы Джет
147,81
Системный интегратор
Поделиться публикацией

Комментарии 8

    +1
    Следует отметить, что обученные на Хайдерабаде модели могут знатно фейлить на Петербурге в силу другого альбедо у типичных крыш в этом регионе (а также в силу другого биоразнообразия сельхозкультур, для NDVI).

    Бэнды в разных спутниках могут не полностью совпадать (например, Band 4 у Landsat не соответствует SWIR).

    Ну и интересно, где можно хоть какие-то размеченные данные нарыть :)
      +1
      Действительно, распределение признаков (в нашем случае, интенсивность в разных диапазонах на уровне пикселя снимка) для «застройки» могут отличаться не только из-за географического расположения и особенностей оборудования, но и даже в рамках одной локации, например, в разное время года.

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

      С сервиса Glovis можно выгрузить различные снимки (по разным типам систем) и слои по нужной территории
      +1
      классификация с помощью Tensorflow
      from tensorflow import keras

      По моему это не совсем корректно
        0
        Ваша версия?
          +1

          классификация с помощью Keras

        +1
        Ещё один этап предварительной обработки заключается в преобразовании матриц признаков из двумерных в трёхмерные, чтобы модель воспринимала каждую строку как отдельный пиксель (отдельный объект обучения).

        Сдается мне что в вашем примере с полносвязной сетью это лишнее действие, т.к. в такой сети всеравно каждый нейрон скрытого слоя соединен со всеми входами.
          0
          последний слой нужен для перевода выхода нейросети в бинарный вектор, соответствующий скору-вероятности (для каждого пикселя) для каждого из двух классов застройки (есть/нет)
            0
            Действительно, в полносвязной сети нейрон следующего слоя принимает все выходы («видит») с предыдущего слоя. В статье речь идет о подготовки векторов-признаков, то есть из исходного трехмерного изображения, получается несколько объектов обучения (равное числу пикселей), где признаками как раз и выступают интенсивности по разным каналам. Лучше всего это объясняет «Схема реструктуризации данных».

          Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.

          Самое читаемое