Cпутниковый снимок в ложных цветах (зелёный, красный, ближний инфракрасный) с пространственным разрешением 3 метра и наложенной маской зданий из OpenStreetMap (спутниковая группировка PlanetScope)

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

В этой статье я постараюсь дать примерное представление о космической съёмке Земли, расскажу о трудностях, с которыми можно столкнуться, начиная работу со спутниковыми снимками: предварительная обработка, алгоритмы для анализа и библиотеки Python для работы со спутниковыми снимками и геоданными. Так что все, кому интересна область компьютерного зрения, добро пожаловать под кат!

Итак, начнём с областей спектра, в которых спутники выполняют съёмку, и рассмотрим характеристики съёмочной аппаратуры некоторых спутников.

Спектральные сигнатуры, атмосферные окна и спутники

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


Отражательная способность сухой травы, асфальта и свежой базальтовой лавы

Как и поверхность Земли, газы в атмосфере также имеют уникальные спектральные сигнатуры. При этом не каждое излучение проходит через атмосферу Земли. Спектральные диапазоны, проходящие через атмосферу, называют «атмосферными окнами», а спутниковые сенсоры настраиваются для измерения в этих окнах.


Атмосферные окна

Рассмотрим спутник с одним из самых высоких пространственных разрешений — WorldView-3 (Пространственное разрешение — величина, характеризующая размер наименьших объектов, различимых на изображении; надир — направление, совпадающее с направлением действия силы гравитации в данной точке).

Характеристики съёмочный аппаратуры WorldView-3

Спектральный диапазон	Наименование	Диапазон, нм	Пространственное разрешение в надире, м
Панхроматический (1 канал; охватывает видимую часть спектра)		450-800	0,31
Мультиспектральный (8 каналов)	Прибрежный	400-450	1,24
	Синий	450-510
	Зелёный	510-580
	Желтый	585-625
	Красный	630-690
	Красный край	705-745
	Ближний инфракрасный	770-895
	Ближний инфракрасный	860-1040
Многополосный в коротком инфракрасном диапазоне (8 каналов)	SWIR-1	1195-1225	3,70
	SWIR-2	1550-1590
	SWIR-3	1640-1680
	SWIR-4	1710-1750
	SWIR-5	2145-2185
	SWIR-6	2185-2225
	SWIR-7	2235-2285
	SWIR-8	2295-2365

Помимо приведённых каналов WorldView-3 имеет ещё 12 каналов, предназначенных специально для атмосферной коррекции, — CAVIS (Clouds, Aerosols, Vapors, Ice, and Snow) с разрешением 30 м в надире и длинами волн от 0,4 до 2,2 мкм.


Пример снимка с WorldView-2; панхроматический канал


Пример снимков с различным пространственным разрешением

Другие интересные спутники — SkySat-1 и его близнец SkySat-2. Интересны они тем, что умеют снимать видео продолжительностью до 90 секунд  над одной территорией с частотой 30 кадров/сек и пространственным разрешением 1,1 м.


Видеосъёмка со спутников SkySat

Пространственное разрешение панхроматического канала — 0,9 м, мультиспектральных каналов (синий, зелёный, красный, ближний инфракрасный) — 2 м.

Ещё некоторые примеры:

1. Группировка спутников PlanetScope выполняет съёмку с пространственным разрешением 3 м в красном, зелёном, синем и ближнем инфракрасном диапазоне;
   
2. Группировка спутников RapidEye выполняет съёмку с пространственным разрешением 5 м в красном, крайнем красном, зелёном, синем и ближнем инфракрасном диапазоне;
   
3. Серия российских спутников «Ресурс-П» выполняет съёмку с пространственным разрешением 0,7-1 м в панхроматическом канале, 3-4 м в мультиспектральных каналах (8 каналов).
   

Гиперспектральные сенсоры в отличие от мультиспектральных делят спектр на множество узких диапазонов (примерно 100-200 каналов) и имеют пространственное разрешение уже другого порядка: 10-80 м.

Гиперспектральные снимки доступны не так широко, как мультиспектральные. Космических аппаратов, на борту которых установлены гиперспектральные сенсоры, немного. Среди них Hyperion на борту спутника NASA EO-1 (выведен из эксплуатации), CHRIS на борту спутника PROBA, принадлежащего Европейскому космическому агентству, FTHSI на борту спутника MightySatII исследовательской лаборатории военно-воздушных сил США, ГСА (гиперспектральная аппаратура) на российских космических аппаратах «Ресурс-П».


Обработанный гиперспектральный снимок с бортового сенсора CASI 1500; до 228 каналов; спектральный диапазон 0,4 – 1 нм


Обработанный снимок с космического сенсора EO-1; 220 каналов; спектральный диапазон 0,4 – 2,5 нм

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

Предварительная обработка снимков

Перед тем, как приступить к анализу снимков, необходимо выполнить их предварительную обработку:

1. Геопривязка;
   
2. Ортокоррекция;
   
3. Радиометрическая коррекция и калибровка;
   
4. Атмосферная коррекция.
   

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

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

Приблизительная геопривязка вычисляется по исходному положению спутника на орбите и геометрии изображения. Уточнение геопривязки выполняется по наземным точкам привязки (Ground Control Points – GCP). Эти контрольные точки ищутся на карте и на снимке, а зная их координаты в разных системах координат, можно найти параметры преобразования (конформное, аффинное, перспективное либо полиномиальное) из одной системы координат в другую. Поиск GСP выполняется при помощи GPS-съёмки [ист. 1 с. 230, ист. 2] либо при помощи сопоставления двух снимков, на одном из которых точно известны координаты GCP, по ключевым точкам.

Ортокоррекция изображения – процесс геометрической коррекции изображений, при котором устраняются перспективные искажения, развороты, искажения вызванные дисторсией объектива и другие. Изображение при этом приводится к плановой проекции, то есть такой, при которой каждая точка местности наблюдается строго вертикально, в надир.

Так как спутники осуществляют съемку с очень большой высоты (сотни километров), то при съёмке в надире искажения должны быть минимальными. Но космический аппарат не может всё время снимать в надире, иначе пришлось бы очень долго ждать момента,  когда он пройдет над заданной точкой. Для устранения этого недостатка спутник «доворачивают», и большинство кадров получаются перспективными. Следует заметить, что углы съемки могут достигать 45 градусов, и при большой высоте это приводит к значительным искажениям.

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

Модель камеры спутника представляется в виде обобщенных аппроксимирующих функций (рациональных полиномов — RPC коэффициентов), а высотные данные могут быть получены в результате наземных измерений, при помощи горизонталей с топографической карты, стереосъемки, по радарным данным или из общедоступных грубых цифровых моделей рельефа: SRTM (разрешение 30-90 м) и ASTER GDEM (разрешение (15-90 м).

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

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


Удаление сбойных пикселей и вертикальных полос

Радиометрическая коррекция выполняется двумя методами:

1. C использованием известных параметров и настроек съемочного прибора (корректировочных таблиц);
   
2. Cтатистически.
   

В первом случае необходимые корректировочные параметры определяются для съемочного прибора на основе длительных наземных и полетных испытаний. Коррекция статистическим методом выполняется путем выявления дефекта и его характеристик непосредственно из самого изображения, подлежащего коррекции.

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

$$

где $$ — энергетическая яркость для спектральной зоны $$;
$$ — сырые значения яркости;
$$ — калибровочный коэффициент;
$$ — калибровочная константа.

Электромагнитное излучение спутника, прежде чем  будет зафиксировано датчиком, пройдет через атмосферу Земли дважды. Имеется два главных эффекта влияния атмосферы — рассеяние и поглощение. Рассеяние происходит в случае, когда имеющиеся в атмосфере частицы и молекулы газа взаимодействуют с электромагнитным излучением, отклоняя его от первоначального пути. При поглощении часть энергии излучения преобразуется во внутреннюю энергию поглощающих молекул, в результате чего происходит нагревание атмосферы. Влияние рассеяния и поглощения на электромагнитное излучение меняется при переходе от одной части спектра к другой.


Факторы, влияющие на попадание отраженной солнечной радиации на сенсоры спутника

Существуют различные алгоритмы выполнения атмосферной коррекции (например метод DOS — Dark Object Subtraction). Входными параметрами для моделей служат: геометрия расположения Солнца и датчика, атмосферная модель для газообразных компонентов, модель аэрозоля (тип и концентрация), оптическая толщина атмосферы, коэффициент поверхностного отражения и спектральные каналы.
Для атмосферной коррекции можно также применять алгоритм удаления дымки с изображения — Single Image Haze Removal Using Dark Channel Prior (реализация).


Пример работы Single Image Haze Removal Using Dark Channel Prior

Индексные изображения

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

Примеры индексных изображений

Название индекса	Формула	Применение
Индекс содержания оксида железа	Красный/синий	Для выявления содержания окислов железа
Индекс содержания глинистых минералов	Отношение значений яркости в пределах среднего инфракрасного канала (CИК). CИК1 / CИК2, где CИК1 это диапазон от 1,55 до 1,75 мкм, CИК2 это диапазон от 2,08 до 2,35 мкм	Для выявления содержания глинистых минералов
Индекс содержания железистых минералов	Отношение значения яркости в среднем инфракрасном (СИК1; от 1,55 до 1,75 мкм) канале к значению яркости в ближнем инфракрасном канале (БИК). СИК1 / БИК	Для выявления содержания железистых минералов
Индекс красноцветности (RI)	Основан на различии отражательной способности красноцветных минералов в красном (К) и зеленом (З) диапазонах. RI = (К — З) / (К + З)	Для выявления содержания оксида железа в почве
Нормализованный дифференциальный индекс снега(NDSI)	NDSI – это относительная величина, характеризуемая различием отражательной способности снега в красном (К) и коротковолновом инфракрасном (КИК) диапазонах. NDSI=(К — КИК) / (К + КИК)	Используется для выделения территорий, покрытых снегом. Для снега NDSI > 0,4
Водный индекс (WI)	WI=0,90 мкм / 0,97 мкм	Применяется для определения содержания воды в растительности по гиперспектральным снимкам
Нормализованный дифференциальный вегетационный индекс (NDVI)	Хлорофилл листьев растений отражает излучение в ближнем инфракрасном (БИК) диапазоне электромагнитного спектра и поглощает в красном (К). Отношение значений яркости в этих двух каналах позволяет четко отделять и анализировать растительные от прочих природных объектов. NDVI=(БИК — К) / (БИК + К)	Показывает наличие и состояние растительности. Значения NDVI варьируют в пределах от -1 до 1; Густая растительность: 0,7; Разряженная растительность: 0,5; Открытая почва: 0,025; Облака: 0; Снег и лед: -0,05; Вода: -0,25; Искусственные материалы (бетон, асфальт): -0,5

Работа со спутниковыми снимками на Python

Одним из форматов, в которых принято хранить спутниковые снимки, является GeoTiff (ограничимся только им). Для работы с GeoTiff на Python можно воспользоваться библиотекой gdal либо rasterio.

Для установки gdal и rasterio лучше воспользоваться conda:

conda install -c conda-forge gdal
conda install -c conda-forge rasterio

Другие библиотеки для работы со спутниковыми снимками легко ставятся через pip.

Чтение GeoTiff через gdal:

from osgeo import gdal
import numpy as np

src_ds = gdal.Open('image.tif')
img = src_ds.ReadAsArray()
#height, width, band
img = np.rollaxis(img, 0, 3)

width = src_ds.RasterXSize
height = src_ds.RasterYSize

gt = src_ds.GetGeoTransform()
#геокоординаты крайних точек снимка
minx = gt[0]
miny = gt[3] + width*gt[4] + height*gt[5] 
maxx = gt[0] + width*gt[1] + height*gt[2]
maxy = gt[3]

Систем координат у спутниковых снимков довольно много. Их можно разделить на две группы: географические системы координат (ГСК) и плоские системы координат (ПСК).

В ГСК единицы измерения угловые и координаты представлены десятичными градусами. Наиболее известная ГСК — WGS 84 (EPSG:4326).

В ПСК единицы измерения линейные и координаты могут быть выражены как метры, футы, километры и т.д., поэтому их можно линейно интерполировать. Для перехода из ГСК в ПСК применяют картографические проекции. Одна из наиболее известных проекций — проекция Меркатора.

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

Скрипт для растеризации полигонов:

import numpy as np
import rasterio
#fiona умеет работать с разными форматами хранения полигонов, например shp и geojson
import fiona
import cv2

from shapely.geometry import shape
from shapely.ops import cascaded_union
from shapely.ops import transform
from shapely.geometry import MultiPolygon

def rasterize_polygons(img_mask, polygons):
    if not polygons:
        return img_mask
    
    if polygons.geom_type == 'Polygon':
        polygons = MultiPolygon([polygons])
    
    int_coords = lambda x: np.array(x).round().astype(np.int32)
    exteriors = [int_coords(poly.exterior.coords) for poly in polygons]
    interiors = [int_coords(pi.coords) for poly in polygons
                 for pi in poly.interiors]
    
    cv2.fillPoly(img_mask, exteriors, 255)
    cv2.fillPoly(img_mask, interiors, 0)
    return img_mask


def get_polygons(shapefile):
    
    geoms = [feature["geometry"] for feature in shapefile]
    polygons = list()
    for g in geoms:
            s = shape(g)
            if s.geom_type == 'Polygon':
                if s.is_valid:
                    polygons.append(s)
                else:
                    #устраняем самопересечения полигона
                    polygons.append(s.buffer(0))
            elif s.geom_type == 'MultiPolygon':
                for p in s:
                    if p.is_valid:
                        polygons.append(p)
                    else:
                        #устраняем самопересечения полигона
                        polygons.append(p.buffer(0))

    mpolygon = cascaded_union(polygons)
    
    return mpolygon

#читаем снимок и geojson файл с полигонами
src = rasterio.open('image.tif')
shapefile = fiona.open('buildings.geojson', "r")

#геокоординаты снимка (системы координат снимка и полигонов должны совпадать)
left = src.bounds.left
right = src.bounds.right
bottom = src.bounds.bottom
top = src.bounds.top

height,width = src.shape

#извлекаем мультиполигоны
mpolygon = get_polygons(shapefile)

#так как система координат плоская, то можем получить координаты полигонов уже на самом изображении
mpolygon = transform(lambda x, y, z=None: (width * (x - left) / float(right - left),                                                       height - height * (y - bottom) / float(top - bottom)), mpolygon)

#рисуем маску
real_mask = np.zeros((height,width), np.uint8)
real_mask = rasterize_polygons(real_mask, mpolygon)

Во время дешифрования снимков может понадобиться операция пересечения полигонов (например, автоматически разметив здания по карте, хочется также автоматически убрать разметку зданий в тех местах снимка, где есть облака). Для этого есть алгоритм Уайлера-Атертона, но он работает только с полигонами без самопересечений. Для устранения самопересечений нужно проверить пересечение всех ребер с другими ребрами полигона и добавить новые вершины. Эти вершины разобьют соответствующие им ребра на части. В библиотеке shapely есть метод для устранения самопересечений — buffer(0).

Для перевода из ГСК в ПСК можно воспользоваться библиотекой PyProj (либо сделать это в rasterio):

#переводим координаты из географической системы координат epsg:4326 в координаты плоской системы координат epsg:32637 (ПСК для Москвы) 
def wgs84_to_32637(lon, lat):
    from pyproj import Proj, transform
    
    inProj = Proj(init='epsg:4326')
    outProj = Proj(init='epsg:32637')
    
    x,y = transform(inProj,outProj,lon,lat)
    
    return x,y

Метод главных компонент

Если снимок содержит более трех спектральных каналов, можно создать цветное изображение из трех главных компонент, уменьшив тем самым объем данных без заметной потери информации.

Такое преобразование также проводят для серии разновременных снимков, приведенных в единую систему координат, для выявления динамики, которая ярко проявляется в одной или двух компонентах.

Скрипт для сжатия 4-х канального изображения до 3-х канального:

from osgeo import gdal
from sklearn.decomposition import IncrementalPCA
from sklearn import preprocessing
import numpy as np
import matplotlib.pyplot as plt

#читаем снимок, меняем порядок на height, width, band и вырезаем часть изображения, чтобы PCA быстрее справился с данными
src_ds = gdal.Open('0.tif')
img = src_ds.ReadAsArray()
img = np.rollaxis(img, 0, 3)[300:1000, 700:1700, :]
height, width, bands = img.shape

#готовим данные для PCA
data = img.reshape((height * width, 4)).astype(np.float)
num_data, dim = data.shape
pca = IncrementalPCA(n_components = 3, whiten = True)

#cтандартизация данных (достаточно и центрирования)
x = preprocessing.scale(data)

#сжимаем каналы и приводим к виду удобному для визуализации
y = pca.fit_transform(x)
y = y.reshape((height, width, 3))

#нормализуем яркость пикселей
for idx in range(0, 3):
    band_max = y[:, :, idx].max()
    y[:, :, idx] = np.around(y[:, :, idx] * 255.0 / band_max).astype(np.uint8)

Снимок с группы спутников PlanetScope (red, green, blue без цветовой коррекции)


Снимок с группы спутников PlanetScope (green, red, near infrared)


Снимок, полученный при помощи метода главных компонент

Метод спектрального разделения (Spectral Unmixing)

Метод спектрального разделения применяют для распознавания на снимках объектов, размер которых значительно меньше размера пикселя.

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


Смешанная спектральная кривая


Дешифрованный снимок

Сегментация спутниковых снимков

На данный момент state-of-the-art результаты в задачах бинарной сегментации изображений показывают модификации U-Net модели.


Архитектура U-Net модели (размер изображения на выходе меньше размера входного изображения; делается это из-за того, что сеть хуже предсказывает на краях изображения)

Автор U-Net разработал архитектуру на основе другой модели — Fully convolutional network (FCN), особенностью которой является наличие только свёрточных слов (не считая max-pooling).
U-Net отличается от FCN тем, что добавлены слои, в которых max-pooling заменён на up-convolution. Таким образом, новые слои постепенно увеличивают выходное разрешение. Также признаки с энкодер-части комбинируются с признаками из декодер-части, чтобы модель могла делать более точные предсказания за счёт дополнительной информации.

Модель, в которой отсутствует проброс признаков из энкодер-части в декодер-часть, называется SegNet и на практике показывает результаты хуже U-Net.


Max-pooling


Up-convolution

Отсутствие в U-Net, Segnet и FCN слоёв, привязанных к размеру изображения, позволяет подавать на вход одной и той же сети изображения разного размера (размер снимка должен быть кратным количеству фильтров в первом свёрточном слое).

В keras это реализуется так:

inputs = Input((channel_number, None, None))

Обучение, как и предсказание, можно проводить либо на фрагментах изображения (кропах), либо на всём изображении целиком, если позволяет память GPU. При этом в первом случае:
1) больше размер батчей, что хорошо скажется на точности модели, если данные зашумлены и неоднородны;
2) меньше риск переобучения, т.к. данных гораздо больше, чем при обучении на изображениях полного размера.

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

U-Net  –  простая и мощная архитектура для задачи бинарной сегментации, на github можно найти не одну реализацию для любого DL фреймворка, но для сегментации большого количества классов эта архитектура проигрывает другим архитектурам, например, PSP-Net. Здесь можно прочитать интересный обзор по архитектурам для семантической сегментации изображений.

Определение высоты зданий

Высоту зданий можно определять по их теням. Для этого необходимо знать: размер пикселя в метрах, длину тени в пикселях и sun (solar) elevation angle (угол солнца над горизонтом).


Геометрия солнца, спутника и здания

Вся сложность задачи в том, чтобы как можно точнее сегментировать тень здания и определить длину тени в пикселях. Проблем добавляет ещё и наличие облаков на снимках.

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



Пример работы алгоритма определения высоты зданий по тени

На снимках с группы спутников PlanetScope с пространственным разрешением 3 м ошибка определения высоты зданий по MAPE (mean absolute percentage error) составила ~ 30%. Всего было рассмотрено 40 зданий и один снимок. Однако на субметровых снимках исследователи получили ошибку всего 4-5%.

Заключение

Дистанционное зондирование Земли даёт много возможностей для проведения аналитики. Например, такие компании, как Orbital Insight, Spaceknow, Remote Sensing Metrics, OmniEarth и DataKind, на основе съёмки со спутников выполняют мониторинг производства и потребления в США, анализ урбанизации, трафика, растительности, экономики и т.д. При этом сами снимки становятся всё более доступными. Например, снимки со спутников Landsat-8 и Sentinel-2 с пространственным разрешением больше 10м находятся в открытом доступе и постоянно обновляются.

В России компании Совзонд, СканЭкс, Ракурс, Гео-Альянс и Северная Географическая Компания также проводят геоаналитику по спутниковым снимкам и являются официальными дистрибьюторами компаний-операторов спутников ДЗЗ: ОАО «Российские космические системы» (Россия), DigitalGlobe (США), Planet (США), Airbus Defence and Space (Франция-Германия) и др.

P.S. Вчера мы запустили онлайн-соревнование по спутниковым снимкам, состоящее из двух задач:

1. Сегментация зданий и автомобилей и подсчёт автомобилей;
   
2. Определение высот зданий.
   

Так что все, кто кому интересна область компьютерного зрения и дистанционного зондирования Земли, присоединяйтесь!

Мы планируем провести ещё одно соревнование по снимкам Роскосмоса, Airbus Defence and Space и PlanetScope.