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

Меня зовут Сергей Кукуруз, я руковожу ML‑проектами в центре технологий для общества Yandex Cloud. В этой статье расскажу, как мы совместно со студентами Школы анализа данных (ШАД), а также с движением добровольцев «СтопБорщевик» автоматизировали этот процесс с помощью машинного обучения. Мы подробно разберём технический пайплайн: от нормализации GeoTIFF‑файлов и извлечения признаков (индекс CIVE) до обучения модели в CatBoost.

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

Коротко: что за модель и зачем мы её подготовили

Наша ML‑модель — это часть сервиса для экологического мониторинга, который мы сделали в Yandex Cloud вместе со студентами ШАДа. Мы разработали его, чтобы помочь владельцам земельных участков (особенно крупных) контролировать распространение опасных растений, в том числе борщевика. 

Сервис анализирует спутниковые снимки и выдаёт изображение с размеченными зарослями борщевика и площадью заражения. Также можно выгрузить векторный файл (например, GeoJSON) для работы в ГИС‑программах. 

Итог работы сервиса: красным обозначены заросли опасного растения
Итог работы сервиса: красным обозначены заросли опасного растения

Мы обучили модель на 55 больших спутниковых снимках (это дало датасет в 10 тысяч изображений). Но прежде чем я расскажу о процессе, нужно дать контекст: в чём особенность и опасность борщевика, а также почему мы создали сервис именно таким.

Почему борщевик Сосновского — это проблема

Раньше это растение культивировали на корм для животных: оно вырастает до 3–5 метров в высоту, даёт много зелёной массы, неприхотливо и быстро размножается. Но довольно быстро обнаружился недостаток: при попадании сока на кожу появлялись ожоги. Кроме того, после откорма борщевиком молоко и мясо давали неприятный запах и привкус.

Заросли борщевика
Заросли борщевика

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

К 2017 году борщевик Сосновского был обнаружен в 18 регионах на площади около 170 тысяч гектаров. Чтобы легче представить — это квадрат с длиной стороны в 41 км. К 2022 году растение было зафиксировано уже в 40 регионах, но площадь заражения не оценивалась. 

Бороться с вредителем сложно из‑за нескольких факторов: 

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

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

  • У владельцев не всегда хватает ресурсов. Если участок небольшой, то обнаружить и ликвидировать растение — посильная задача. А вот администрациям парков и других крупных участков можно только посочувствовать. 

Высота борщевика Сосновского — более метра, однако могут встречаться и экземпляры высотой до 3 и даже 4–5 метров
Высота борщевика Сосновского — более метра, однако могут встречаться и экземпляры высотой до 3 и даже 4–5 метров

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

Мы подумали, что нужен инструмент, который будет: 

  • Простым. Чтобы им мог пользоваться любой человек, без технических навыков. 

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

  • Бесплатным. Чтобы им могли пользоваться все. Значит, оно должно быть недорогим в разработке и обслуживании. 

Он позволил бы быстро решать сразу несколько задач: 

  • Находить новые заражённые области. 

  • Оценивать результаты ликвидации.

  • Аккумулировать данные — для научных работ или принятия решений.

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

Сервис под капотом

Если нам нужно обучить модель корректно относить объект к одной из двух категорий «вредитель» / «не вредитель» — это задача пиксельной классификации (семантической сегментации). Для неё хорошо подходит алгоритм градиентного бустинга (почему — расскажу чуть позже). 

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

  1. Выделили по контуру области с борщевиком. 

  2. Считали вегетативный индекс и подбирали для него порог: если значение было выше, область закрашивалась, если ниже — нет. 

Пример такой разметки
Пример такой разметки

Поэтому нам и потребовались именно спутниковые снимки: они не только дают значения RGB‑каналов (по ним рассчитывается вегетативный индекс CIVE), но и обладают метаданными, которые позволяют привязать закрашенную область к определённым координатам.

Далее мы обучили модель с помощью библиотеки CatBoost. Ещё 2000 фрагментов спутниковых снимков использовали для тестирования. 

После обучения у нас получился такой пайплайн: 

  1. Чтение спутникового снимка. 

  2. Нормализация данных.

  3. Разбивка на тайлы.

  4. Извлечение признаков.

  5. Классификация с помощью CatBoost.

  6. Векторизация.

  7. Представление на карте.

Давайте разберём подробнее каждый шаг. 

Чтение спутникового снимка

Пользователь загрузил GeoTIFF‑файл. Далее мы открываем его с помощью библиотеки rasterio — она нужна для чтения, записи и обработки растровых геоданных. Проще говоря, она позволяет работать с географическими изображениями, где каждый пиксель привязан к координатам на Земле.

Мы читаем BGR‑каналы и конвертируем их в RGB. Одновременно сохраняем метаданные: аффинное преобразование (привязка к координатам), CRS (система координат) и размер изображения.

def read_image(path: Path) -> tuple[np.ndarray, Any, Any, Any]:
    # ...
    with rasterio.open(path) as src:
        bgr_data = src.read()
        return np.dstack([bgr_data[2], bgr_data[1], bgr_data[0]]), src.transform, src.crs, src.shape

Нормализация данных

Спутниковые снимки часто содержат аномально яркие или тёмные пиксели. Это могут быть: 

  • блики (засветки от облаков или воды);

  • мёртвые или горячие пиксели сенсора;

  • глубокие тени зданий или рельефа;

  • артефакты сжатия.

Из‑за выбросов вариативность значений может быть огромной, что способно сильно ухудшить качество классификации. Поэтому мы проводим нормализацию значений яркости каждого из трёх каналов (R, G, B) по отдельности.

Алгоритм такой: 

  • для каждого канала берётся 0,3-й перцентиль (min_value) и 99,7-й перцентиль (max_value); 

  • всё, что ниже 0,3%, обрезается до min_value;

  • всё, что выше 99,7%, обрезается до max_value; 

  • оставшийся диапазон линейно растягивается в [0:255].

def normalize_channel(channel: np.ndarray, percent: float = 0.3) -> np.ndarray:
    min_value, max_value = np.percentile(channel, percent), np.percentile(channel, 100 - percent)
    stretched = (channel - min_value) / (max_value - min_value)
    return (np.clip(stretched, 0, 1) * 255.0).astype(np.uint8)

Тут ещё стоит отметить, что на этапе обучения мы использовали инструмент RandomUnderSampler для борьбы с дисбалансом в распределении классов. Он берёт все объекты миноритарного класса, случайно удаляет часть объектов мажоритарного класса и делает классы примерно одинаковыми по размеру.

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

Разбивка на тайлы

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

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

def get_mask(image, height, width, hogweed_model, size=256):
    mask = np.zeros((height, width)).astype(np.uint8)
    for i in range(0, height, size):
        for j in range(0, width, size):
            # ... обработка краёв ...
            small_image = image[v:v + size, h:h + size, :]
            x = collect_features(small_image)
            positive_probs = hogweed_model.predict_proba(x)[:, 1]
            threshold = 0.95
            predictions = (positive_probs >= threshold).astype(int)
            small_mask = predictions.reshape((256, 256))
            # ...

Извлечение признаков

Этот этап ещё называется feature engineering. Для каждого тайла формируется таблица, где каждый пиксель — строка с 20 признаками:

  • 4 пиксельных признака: R‑, G‑, B‑значение конкретного пикселя плюс его значение вегетационного индекса CIVE.

  • 16 глобальных (контекстных) признаков тайла. Для каждого из трёх каналов вычисляется среднее значение по всему тайлу плюс: 

    • среднее (mean); 

    • стандартное отклонение (std);

    • медиана (median);

    • минимальное (min);

    • и максимальное (max) значение.

Вот как выглядит таблица:

#

Признак

Тип

Описание

1–3

R, G, B

Пиксельный

Значения каналов конкретного пикселя

4

CIVE

Пиксельный

Вегетационный индекс этого пикселя

5–9

R_mean, R_std, R_median, R_min, R_max

Глобальный

Статистики красного канала по тайлу 256 × 256

10–14

G_mean, G_std, G_median, G_min, G_max

Глобальный

Статистики зелёного канала по тайлу

15–19

B_mean, B_std, B_median, B_min, B_max

Глобальный

Статистики синего канала по тайлу

20

CIVE_mean

Глобальный

Средний CIVE по тайлу

Индекс CIVE (англ. Color Index of Vegetation Extraction) рассчитывается по формуле: 0,441  R − 0,881  G + 0,385 * B + 18,78745. Он хорошо выделяет зелёную растительность.

Глобальные признаки — это «контекст»: они одинаковы для всех пикселей в одном тайле и дают модели информацию об окружающей обстановке (лес, поле, город). 

Пиксельные признаки — конкретный цвет и «зелёность» данной точки. Комбинация позволяет бустингу отделить «зелёный пиксель в окружении зелени» (лес) от «зелёный пиксель определённого оттенка в контексте с определёнными статистиками» (борщевик).

def collect_features(image):
    global_features = spectral_features(image)
    cive, cive_mean = calculate_cive(image)
    global_features.append(cive_mean)
    cive_expanded = np.expand_dims(cive, axis=-1)
    pixel_features = np.concatenate([image, cive_expanded], axis=-1)
    pixel_features = pixel_features.reshape(-1, 4)
    combined_features = add_global_features_to_pixels(pixel_features, global_features)
    return combined_features

Классификация

Обученная с помощью CatBoost модель получает признаки и предсказывает вероятность борщевика для каждого пикселя. Мы применяли высокий порог в 0,95 — чтобы минимизировать ложные срабатывания. 

В результате мы получим бинарную маску: вектор из 0 и 1, который показывает, какие пиксели выбраны, а какие нет.

Векторизация 

На этом этапе мы конвертируем бинарную маску в полигоны с помощью уже знакомой вам библиотеки rasterio. Далее упаковываем их в GeoDataFrame с помощью библиотек shapely и geopandas. 

Затем с помощью операции spatial join мы определяем, в каком регионе находится каждый полигон борщевика, сопоставляя его с шейп‑файлом регионов РФ, после чего вычисляем площадь каждого полигона в квадратных метрах.

Что такое spatial join?

Это операция, которая соединяет два геопространственных набора данных по их взаимному расположению (пересечению, вложенности и так далее)

Вот код:

def get_features(mask: np.ndarray, transform: Any, connectivity=4) -> Generator[Polygon, None, None]:
    for geom, value in features.shapes(mask, transform=transform, connectivity=connectivity):
        if value:
            yield shape(geom)

def gdf_from_features(vec_features: Iterable[Polygon], feat_crs: Any) -> gpd.GeoDataFrame:
    return gpd.GeoDataFrame(vec_features, columns=["geometry"], crs=feat_crs)

Готово! После этого данные остаётся выгрузить на карту и представить пользователю. 

Детектированная область выделена жёлтым внизу
Детектированная область выделена жёлтым внизу

Для оценки точности алгоритмов объектного обнаружения используется метрика IOU (Intersection Over Union). Она измеряет степень перекрытия между предсказанной областью и эталонной разметкой. У нашего решения этот показатель 0,75, что довольно неплохо. 

Пара слов о градиентном бустинге

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

  • Размер датасета. Для обучения нейросети сегментации (U‑Net, DeepLab и пр.) нужны десятки тысяч размеченных снимков. У нас было только 55 спутниковых изображений, которые дали нам ~3,6 млн пикселей‑строк, но для табличного бустинга это комфортный объём.

  • Ручной feature engineering вместо эмбеддингов. Вместо того чтобы заставлять нейросеть самой учить представления, мы сформулировали задачу как табличную, выделив 20 осмысленных признаков на пиксель. Градиентный бустинг — это де‑факто стандарт для табличных данных, он стабильно побеждает нейросети на подобных задачах.

  • Возможность посмотреть feature importance. То есть понять, какие признаки важны. С нейросетью это сложнее.

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

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

Если что‑то упустил — напишите в комментариях, для чего вы бы применили классификацию объектов на спутниковых снимках?

Что в итоге

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

Сейчас на карте размечены очаги заражения общей площадью 421 га в 17 регионах европейской части страны

Сегодня территорию Москвы и Подмосковья уже удалось проанализировать полностью. Планируем, что к лету 2026 года мониторингом будет охвачено 100 тысяч км² территории в Тверской и Ярославской областях: здесь мы рассчитываем на вклад пользователей и волонтёров. 

Кроме того, сервис можно масштабировать и на другие опасные инвазивные растения: например, для эхиноцистиса лопастного и золотарников канадского и гигантского.

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