Как с помощью трансферного обучения обнаружить вулканы на Гавайях

Геопространственная сегментация изображений с использованием топографических данных и трансферного обучения

Data Science — это не только данные о клиентах. К старту нашего флагманского курса рассмотрим пример геопространственной семантической сегментации, где с помощью данных цифровой модели рельефа отобразим шлаковые конусы на Гавайях.

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

Разберёмся сначала, что такое «шлаковые конусы», и как с помощью QGIS быстро создать набор данных с размеченными данными. Затем проведём обработку данных и отбор/обучение модели, а после протестируем обученную модель на новых данных.

Что такое «шлаковые конусы»?

«Шлаковые конусы — это простейший тип вулкана. Они формируются частицами и сгустками застывшей лавы, выброшенными из одного жерла. <…> Большинство шлаковых конусов венчает чашеобразный кра́тер». Геологическая служба США.

Вместо тысячи слов:

Итак, нам нужны конусообразные объекты с крутыми склонами и круглым кратером посредине. Посмотрим, как они выглядят в наборе данных.

Загрузка и разметка данных

Для набора данных используем предоставленные USGS¹ данные ЦМР — топографические данные с крупнейшего Гавайского острова, где видна морфология шлаковых конусов. Нанесём на карту шлаковые конусы спящего вулкана Мауна-Кеа, чтобы построить модель отображения конусов соседнего вулкана Хуалалаи.

Для разметки шлаковых конусов используем QGIS — отличное приложение для работы с геопространственными данными. Визуализируем шлаковые конусы, загружая данные в формате .geotiff и получая на выходе контурные горизонтали с интервалами 10 метров (подсказка: ищите небольшие плотно расположенные концентрические окружности):

Это шлаковые конусы на основе уже определённых пространственных критериев. Нанесём их на карту, чтобы создать размеченные маски. Для простоты будем считать шлаковые конусы круглыми. Это значит, что метками будут круговые многоугольники. Создадим шейп-файл этих фигур и начнём загружать его метки и данные ЦМР geotiff:

#Loading a geotiff shapefile
def readimage():
 print("reading image…")
 with rasterio.open(path_image, "r") as ds:
 arr = ds.read()
 arr = np.transpose(arr, (1, 2, 0))
 #Clip negative values to 0 
 arr = arr.clip(0)
 print(arr.shape)
 return arr
#Loading the shapefile mask
def readmask(link):
print("reading mask…")
geo = gpd.read_file(link)
with rasterio.open(path_image) as src:
raster = src.read()
geo = geo.to_crs(src.crs)
out_image, out_transform = rasterio.mask.mask(src, geo.geometry, filled=True)
masks = out_image[0,:,:]
#Set the mask labels to 1 and the rest to 0
masks[np.where(masks<=0)] = 0
masks[np.where(masks>0)] = 1
masks = masks.astype(np.int8)
masks = np.expand_dims(masks, axis=2)
return masks

Предварительная обработка данных

В массиве ЦМР у нас весь остров. Уберём всё вокруг размеченной маски Мауна-Кеа, а этот массив из 2D переделаем в 3D с тремя каналами. Зачем? Для совместимости с моделью сегментации, которая использует предварительно обученный энкодер. Этому энкодеру обычно требуется трёхканальное входное изображение, например, цветное:

#Crop both the image and masks to the same Mauna Kea area
image = image[ymin:ymax, xmin:xmax,:]
masks = masks[ymin:ymax, xmin:xmax,:]

#Stack the image array and normalise
image = np.dstack((image, image, image))
image = (image — image.min())/(image.max() — image.min())
 
original_size = image.shape

Получится форма (4000, 6500, 3). Визуализируем изображение и маску:

Дублировав одно и то же изображение трижды, применим фильтры. Чтобы шлаковые конусы выделялись на фоне окружения, увеличим контраст. И добавим фильтр Собеля, чтобы чётче обозначались их круговые формы:

#Contrast enhancing
image_eq = exposure.equalize_adapthist(image, clip_limit=0.05)
#Sobel filter and normalising
image_sobel = filters.sobel(np.squeeze(image))
image_sobel = (image_sobel — image_sobel.min())/(image_sobel.max() — image_sobel.min())
#concatenate standard image, equalised and sobel together
images = np.dstack((image[:,:,0], image_sobel[:,:,0], image_sobel[:,:,0]))

Из этого большого массива мы использовали три канала. Разделим массив на мелкие массивы, которые отправим в модель. Вот результат разделения:

Количество образцов, высота, ширина, количество каналов:

#Making image tiles
size = 224
step = int(size/2)
patch_arr = skimage.util.view_as_windows(image, (size, size, layer.shape[2]), step = step)
output = patch_arr.reshape((-1,) + (size, size, layer.shape[2]))

Процесс повторяется для входного изображения и масок, чтобы они дополняли друг друга. Итог — такие входные данные:

Какие интересные изображения! Итоговая форма массива будет такой: (1938, 224, 224, 3). Примеров для создания модели не так много, ещё и поэтому нужна предварительно обученная модель. Последний этап — разделение данных на обучающую и контрольную выборки:

x_train, x_val, y_train, y_val = train_test_split(images, masks, test_size=0.2, shuffle=True, random_state=123)
print(x_train.shape)
print(x_val.shape)
print(y_train.shape)
print(y_val.shape)
y_train = tf.cast(y_train, tf.float32)
y_val = tf.cast(y_val, tf.float32)

Построение и обучение модели

Воспользуемся, по-видимому, лучшим вариантом архитектуры сегментации — моделью UNET — с предварительно обученным энкодером из Keras. Наш энкодер — InceptionResNetV2, на наборе данных imagenet у него хорошая производительность, хотя подойдёт любая предварительно обученная модель:

input_shape = (size, size, 3)
inception = InceptionResNetV2(include_top = False, weights = "imagenet", input_tensor = input_shape)
inception.summary()
layer_names = [layer.name for layer in model.layers]

Нужно добавить соединения быстрого доступа между энкодером и декодером. Но сначала создадим декодер UNET, загрузив предварительно обученный энкодер с входными размерами и получив список имён слоёв для него.

Подходящий декодер будет создан после установки include_top в False. С помощью декодера выходные данные энкодера возвращаются к исходным размерам входного изображения, а мы получаем соответствующую маску сегментации:

x = Conv2DTranspose(num_filters, (2,2), strides=2, padding="same")(inputs)
x = Conv2D(size, 2, padding="same", dilation_rate = 1, kernel_initializer = "he_normal")(x)
x = BatchNormalization()(x)
x = Dropout(0.2)(x)
x = Activation("LeakyReLU")(x)

Складывая эти блоки один на другой от большего к меньшему фильтру, получим декодер. Соединения быстрого доступа добавляются к слоям активации, где размер входного слоя — 224, 112, 56 и 28. В зависимости от используемой модели размер слоя может отличаться. Тогда потребуется заполнение:

skip_connection_2 = inception.get_layer(index = 3).output 
skip_connection_2 = ZeroPadding2D(( (1, 0), (1, 0))  
(skip_connection_2)

Посмотрите: четвёртый слой (индекс 3) надо заполнить до нужной формы. Благодаря соединениям быстрого доступа энкодер и декодер связываются конкатенацией выходных данных указанных уровней активации и транспонированных блоков декодера:

x = Concatenate()([x, skip_connection_2])

Создав декодер, завершаем модель конечным выходным слоем, используя сигмоидную функцию активации, ведь наша задача — бинарная сегментация, то есть определение, со шлаковым конусом мы имеем дело или нет:

outputs = Conv2D(1, (1,1), padding="same", activation="sigmoid")(last_layer)

При компиляции модели в функции потерь должна учитываться несбалансированность классов. Количество шлаковых конусов относительно фоновой области небольшое, поэтому воспользуемся Tversky Loss. 

Есть много разных функций потерь. Чтобы измерить процент перекрывающихся пикселей между входной маской и маской сегментации, возьмём коэффициент сходства Жаккара:

model.compile(tf.keras.optimizers.Adam(5e-5), 
loss= tversky_loss, metrics=[jaccard_coefficient],)
reduce_lr = ReduceLROnPlateau(monitor="loss", factor=0.2, patience=15, verbose=1, mode="max", min_lr=1e-6)
early = tf.keras.callbacks.EarlyStopping(monitor="val_loss", patience=10)
history = model.fit(x_train, y_train, epochs=100, batch_size = 32, validation_data = (x_val,y_val), callbacks = [reduce_lr, early])

При низкой скорости обучения 1e-5 и ReduceLROnPlateau в течение 100 эпох моделью достигнуто чуть более 80% IoU валидации (пересечения по объединению). Не так уж плохо. 

Улучшить показатель способна более совершенная стратегия разметки, без упрощений в виде круговых фигур в качестве меток, или другая модель энкодера. Можно поэкспериментировать и с замораживанием разных слоёв энкодера.

Тестирование модели на другом вулкане

Посмотрим, как хорошо модель сегментации распознаёт шлаковые конусы вулкана Хуалалаи:

Проходим те же этапы предварительной обработки входных данных, чтобы сделать их совместимыми с моделью. Снова убираем всё лишнее, применяем фильтры, разделяем массив на тайлы и получаем форму (420, 224, 224, 3).

Увеличивать размер образцов здесь не нужно, поэтому обойдёмся без получения плиток изображений. Вот какие здесь входные тайлы:

После использования модели для формирования прогнозов:

y_pred = model.predict(tiled_Hualalai_image)

Посмотрите на сегментированные тайлы, прежде чем снова объединять тайлы изображений, чтобы сделать маску сегментации с размерами исходного изображения:

Прогнозы в норме, но у модели проблемы с сегментированием небольших шлаковых конусов. Они могут быть связаны с размером шлаковых конусов в обучающей выборке Мауна-Кеа, которые обычно больше.

Посмотрим на изображение сегментации, наложенное на исходное входное изображение. Чтобы обнаружить шлаковые конусы было легче, наложим спрогнозированные маски сегментации на изображение с фильтром Собеля.

Шлаковые конусы изображены как маленькие точки или кольца:

Изображение выглядит намного лучше:

• Расположение масок и шлаковых конусов вдоль вулкана в целом совпадает (северо-запад и юго-восток).

• Сегментированы шлаковые конусы различных размеров, а не только большие.

• Все шлаковые конусы по форме в целом соответствуют маскам сегментации.

Заключение

Мы разобрали пример использования сегментации изображений с геопространственными наборами данных для идентификации шлаковых конусов на Гавайях. Преобразуя маски сегментации из входных данных ЦМР обратно в исходные координаты, можно сделать из них наборы данных для дальнейшего анализа (например, анализа распределения размеров шлаковых конусов или изменения их плотности). Спасибо за внимание!

Продолжить изучение Data Science или Python вы сможете на наших курсах:

• Профессия Data Scientist

• Профессия Fullstack-разработчик на Python

Узнайте подробности здесь.