Как я искал боровики и дошёл до спутниковой классификации леса

Его самостоятельность подтверждается молекулярным филогенетическим исследованием Beugelsdijk et al. (2008): анализ участков ДНК (дезоксирибонуклеиновой кислоты) — ITS (Internal Transcribed Spacer, внутреннего транскрибируемого спейсера рибосомной РНК) и GAPDH (Glyceraldehyde-3-Phosphate Dehydrogenase, гена фермента глицеральдегид-3-фосфатдегидрогеназы) — надёжно выделил этот вид в отдельную эволюционную ветвь внутри секции Boletus.

Для моей задачи особенно полезна экологическая характеристика: вид образует эктомикоризу преимущественно с сосной обыкновенной (Pinus sylvestris) и приурочен к сухим, кислым, дренированным песчаным и супесчаным почвам. Среди морфологических маркеров — мясистая бугорчатая шляпка тёмного винно-красно-коричневого оттенка и выраженная рельефная сеточка в верхней части массивной ножки.

При этом исследование предостерегает от поверхностных выводов:

1. Связь с сосной не абсолютна: зафиксированы находки микоризы в еловых насаждениях.

2. Белый гриб под сосной — не обязательно B. pinophilus. Авторы описали образцы, собранные на песчаных почвах под соснами, которые фенотипически и по биотопу напоминали сосновый боровик, но по ДНК оказались экологической формой обычного белого гриба — Boletus edulis (в микологических статьях после первого упоминания род сокращается: B. pinophilus, B. edulis).

Извлечение почв: динамический веб-харвестинг

На портале soil-db.ru оцифрованная институтом им. Докучаева карта Фридланда не отдаётся готовым файлом. Она работает в веб-интерфейсе через фреймворк Livewire: при панорамировании окна браузер отправляет метод loadFeatures, а сервер шлёт геометрию полигонов в событии livewire-map:add-features.

Я написал сборщик на базе протокола Chrome DevTools Protocol (CDP), который прошёл всю Ленинградскую область сеткой из 55 перекрывающихся окон зума 9, вытянул векторные контуры и слил их в единый GeoJSON (RFC 7946) — формат хранения географических объектов на базе синтаксиса JSON.

Снапшот состояния, скопированный из активного окна, содержал кэш loadedIds — сервер считал эти участки уже загруженными и наотрез отказывался отдавать центральное окно Петербурга и Карельского перешейка. Пришлось сбрасывать идентификаторы и догружать недостающие 60 полигонов отдельным проходом. Итог: 364 полигона, 25 почвенных разностей. Главный фаворит — иллювиально-железистые подзолы (), визитная карточка сосновых боров на песках.

Маска леса: оконный доступ к COG

Для лесного покрова я взял продукт CGLS-LC100 (Buchhorn et al., 2020) за 2019 год. Это глобальный растр с шагом 100 м. Качать весь файл (около гигабайта) не потребовалось: GeoTIFF был оформлен по стандарту Cloud-Optimized GeoTIFF (COG).

С помощью библиотеки GDAL и её виртуальной файловой системы /vsicurl в связке с Python-библиотекой rasterio можно читать растр напрямую из удалённого хранилища через HTTP-запросы диапазонов байтов (range-requests). Скрипт вырезал ровно географическую рамку Ленобласти.

Далее пиксели класса 1 (вечнозелёный хвойный лес) были векторизованы функцией rasterio.features.shapes: получилось 3044 контура площадью более 1 км².

Ориентироваться только на теги OSM (leaf_type=needleleaved — тег, означающий иглолистную хвою вообще, а не конкретно сосну) не позволяла фрагментарность покрытия.

1. Экспертная оценка почв: типам почв Фридланда присвоили рейтинг пригодности (). В отбор пошли почвы от (прежде всего те самые иллювиально-железистые подзолы).

2. Растеризация: вектор почв перевели в растровую сетку, идентичную сетке CGLS-LC100.

3. Логическое пересечение: бинарная маска «хвойный лес почва ».

4. Полигонизация и отсев: преобразовали маску обратно в вектор через rasterio.features.shapes, выкинув мелкий мусор площадью до 0,3 км². Осталось 3265 зон общей площадью 10 240 км².

5. Зональная статистика подложки: для каждого полигона рассчитали зональное среднее доли песка в подпочве на глубине 60 см по глобальной сетке OpenLandMap / SoilGrids (Hengl et al., 2017) (разрешение 250 м).

6. Разделение на сосну и ель: если средняя песчаность превышала 60%, полигон получал флаг pine=true (условный сосняк); если ниже — отправлялся в категорию «прочий хвойный» (чаще сырые ельники).

Итог: 1048 «сосновых» зон (4348 км²) и 2217 прочих хвойных.

Флаг pine=true был смелым авансом: он означал не то, что модель распознала сосну, а лишь то, что хвойный лес стоит на относительно лёгкой песчаной почве. Пересечение разномасштабных данных (почвы 1:2,5 млн, песок 250 м, лес 100 м) не делает границы точнее их самого грубого компонента. Звёзды означали приоритет проверки ногами, а не плотность боровиков на квадратный метр.

Драматичный переезд с Leaflet на OpenLayers

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

• Стандартные веб-карты (OSM, Google) используют проекцию EPSG:3857 (WGS 84 / Pseudo-Mercator, или Web Mercator).

• Спутниковые тайлы Яндекса исторически опираются на эллипсоидальную проекцию Меркатора — EPSG:3395 (WGS 84 / World Mercator).

Если натянуть векторный оверлей из одной проекции на растр другой, на широте Санкт-Петербурга (60° с. ш.) возникает физический сдвиг геометрии на местности примерно в 18 километров. Библиотека MapLibre GL JS отпала сразу — она жёстко привязана к Web Mercator. Единственным движком, способным динамически трансформировать растровые и векторные слои в произвольные проекции через WebGL (браузерный интерфейс аппаратного ускорения графики) и Canvas прямо на мобильном процессоре, оказался OpenLayers. Проверил прототип на смартфоне: подложка Яндекса встала без перекосов, контуры легли пиксель в пиксель.

Векторные тайлы и сжатие в PMTiles

Отдавать мобильному браузеру «сырые» файлы GeoJSON лесов и болот весом по 35 МБ через мобильную сеть на лесной опушке — верная смерть. Данные требовалось тайлировать — разрезать на стандартные пространственные квадраты разного масштаба.

Конвейер генерации оформили в воспроизводимый bash-скрипт:

1. Утилита mapshaper переводила координаты в локально настроенную азимутальную равновеликую проекцию Ламберта LAEA (Lambert Azimuthal Equal Area, сохраняющую площади объектов) с центральной точкой 60° с. ш., 31° в. д. В плоских метрах пороги упрощения и фильтры площади работают геометрически строго, в отличие от искажающихся угловых градусов WGS84.

2. Команда -clean удаляла топологические наложения и осколки растеризации.

3. Упрощение геометрии алгоритмами Дугласа–Пекера (1973) / Висвалингама–Уайатта (1993) с флагом keep-shapes удаляло лишние вершины, не допуская распада мелких замкнутых полигонов. Контуры болот, например, «похудели» с 778 до 333 тысяч вершин.

4. Утилита tippecanoe паковала полигоны в векторные тайлы формата Mapbox Vector Tile (MVT) внутри единого контейнера PMTiles.

Формат PMTiles устроен гениально: это единый файл, содержащий пирамиду тайлов со встроенным индексным заголовком. Браузер через OpenLayers шлёт к серверу запросы диапазонов байтов (HTTP range-requests), скачивая только нужные несколько килобайт на текущий экран, а сервер отвечает кодом 206 Partial Content.

Архивы легли в облачное объектное хранилище Cloudflare R2, а перед ними встал Cloudflare Workers — бессерверный обработчик на JavaScript, который проверяет авторизацию и раздаёт байты с минимальной задержкой.

Болота: метод отрицания

Чтобы не тонуть в прямом смысле слова, через Overpass API (язык запросов к базе OpenStreetMap) мы выкачали все полигоны болот Ленобласти (natural=wetland). Утилита osmtogeojson собрала их в мультиполигоны с внутренними «дырами». Этот слой стал отрисовываться поверх кандидатов ярко-синим и оливковым цветом. Логика: сосновый бор часто тянется узкой песчаной гривой между непроходимыми верховыми топями. Видишь сухое ребро среди синевы — туда и надо идти.

Геология: Девон против ледника

Песок и песчаник — принципиально разные геологические тела. Песок — рыхлый, несцементированный обломочный материал. Песчаник — монолитная осадочная горная порода, где зёрна кварца намертво связаны природным минеральным цементом. Классический пример — знаменитые девонские обрывы реки Оредеж у Белогорки и Ям-Тёсова: этим породам свыше 380 миллионов лет, когда здесь плескалось тёплое мелководное море и плавали панцирные рыбы.

Но современный сосновый бор растёт не на девонских плитах. Его среда — это четвертичный чехол, сформированный последними оледенениями:

1. Морена (): несортированная смесь валунов, суглинков и гравия, которую ледник грёб перед собой. Это субстрат тяжёлый, плотный, на нём чаще поднимается тёмный ельник.

2. Флювиогляциальные отложения (): пески и гравийники, перемытые и отсортированные бешеными потоками талых ледниковых вод. Именно они слагают зандровые поля, озы и камовые холмы — сухие, прекрасно дренированные возвышенности.

3. Озёрно-ледниковые отложения (): накапливались в приледниковых подпруженных озёрах. Они состоят как из песков, так и из тончайших ленточных глин. Сам по себе индекс на карте вовсе не гарантирует сухой грунт под ногами.

4. Эоловые пески (): древние песчаные дюны, навеянные ветром в суровом перигляциальном климате.

Провал «в лоб»

Возникло ощущение, что решение лежит на поверхности: хвойные деревья вечнозелёные, а лиственные сбрасывают листву на зиму. Считаем нормализованный разностный вегетационный индекс NDVI (Rouse et al., 1974; Tucker, 1979):

Логика: сравним летний NDVI и весенний безлистный NDVI спутника Sentinel-2 (Drusch et al., 2012), найдём порог «просадки» зелёности и отделим хвойные от лиственных.

Мы проверили эту гипотезу на пикселях Ленобласти в Google Earth Engine (GEE) (Gorelick et al., 2017). Полный провал:

1. Бесснежная весенняя фенология: гистограмма сезонной разности NDVI внутри леса оказалась строго одногорбой. Никакого разделения на два колокола нет — единым порогом отделить хвою от листвы невозможно. К тому же озимые поля и луговины весной выдавали точно такую же «зелёность», как зрелый сосняк.

2. Зимний «снегометод»: идея искать вечнозелёные кроны над белым снегом разбилась о физику высоких широт. На 60° с. ш. солнце зимой едва поднимается на 6–15° над горизонтом. Кроны отбрасывают гигантские тени: 97% лесных пикселей сколлапсировали в физический минимум индекса.

Казалось, спутники бессильны. Но ресёрч научной литературы (Udali et al., 2021; Blickensdörfer et al., 2024) подсказал: проблема не в спутнике Sentinel-2, а в попытке решить многомерную задачу простым ручным порогом. Нужны размеченные обучающие выборки и машинное обучение.

Я снова обратился к OSM, но теперь не за готовыми сосняками, а за обучающими полигонами с тегом leaf_type. После отсева мелких скверов и дачных посадок (оставили строго контуры площадью от 3 гектаров) набралось 3 769 полигонов: 1 407 чисто хвойных, 451 лиственный и 1 911 смешанных. Смешанные массивы временно отложили, чтобы обучить бинарный классификатор.

Контуры ужали внутрь на 15 метров буфером отсечения граничных эффектов и насэмплировали обучающие точки с шагом 30 метров (76 760 точек). Для каждой точки в GEE собрали вектор признаков:

• Оптика Sentinel-2: летний и весенний NDVI, их разность (фенологический сдвиг), коротковолновые инфракрасные каналы SWIR (Short-Wave Infrared) B11 и B12, чувствительные к влажности хвои и геометрии полога, индекс NDII;

• Радар Sentinel-1: обратное рассеяние в поляризациях VV (вертикальная передача и приём) и VH (вертикальная передача, горизонтальный приём) летом и зимой;

• Рельеф Copernicus GLO-30: цифровая модель высот DEM с разрешением 30 метров, из которой вычислялись абсолютная высота и уклон поверхности.

Облачные запросы падали по таймауту пяти минут, поэтому выборку выгрузили пакетным CSV-файлом, а обучение запустили локально в Python на базе библиотеки scikit-learn. Модель — Random Forest (ансамбль решающих деревьев).

Честная валидация: как не обмануть себя

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

Мы применили пространственно-блочную кросс-валидацию: всю Ленинградскую область разбили на квадраты сетки 0,25°. Модель училась на одних блоках, а тестировалась на 30% совершенно независимых полигонов, которых никогда не видела в глаза.

Результат на 29 190 проверочных точках:

• Overall Accuracy (общая точность): 0,887;

• Каппа Коэна (): 0,772 (отличный показатель с поправкой на случайное угадывание);

• Полнота по классам: лиственный — 0,895, хвойный — 0,878.

Абляционный анализ: что реально работает?

Мы поочерёдно отключали группы признаков, чтобы понять их реальный вклад:

Выводы оказались поразительными:

1. Рельеф — ключевой драйвер: добавление высоты и уклона дало колоссальный скачок точности (+3,6 процентных пункта). Хвойные и мелколиственные леса в тайге чётко привязаны к геоморфологии.

2. Радар оказался лишним: 6 тяжелейших каналов Sentinel-1 дали прирост всего в 0,4 п.п. Для экономии памяти и времени радар был исключён из финального продакшена.

3. Анализ важности признаков по перестановкам (permutation importance) показал: лидирует летний NDVI (0,1635), за ним — уклон рельефа (0,0282) и канал SWIR B11 (0,0277). А та самая сезонная просадка NDVI, которая с треском провалилась в виде одиночного порога, внутри многомерной модели дала весомый вклад (0,0145).

Производственный релиз: компромиссы и оптимизация

Попытка рассчитать всю Ленобласть в 10-метровом разрешении в GEE зависла на 15,5 часов и была снята. На сетке 20 метров число пикселей сократилось в 4 раза, и расчёт уложился в 25 минут.

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

Итоговый слой получил название «Хвойный 20 м» (пометка , бета-версия). При нарезке растровых тайлов файл весил 466 МиБ, превышая лимит загрузки Cloudflare R2 в 300 МиБ. Мы квантовали альфа-канал прозрачности (задающий вероятность хвойного) с помощью Pillow до 16 фиксированных градаций: вес упал до 221 МиБ без визуальной деградации.

┌────────────────────────────────────────────────────────────────────────┐
│                   ИТОГОВАЯ ТОЧНОСТЬ КЛАССИФИКАТОРА                     │
│                                                                        │
│   МАТРИЦА ОШИБОК (ТОЧЕК В ТЕСТЕ: 29 190)                               │
│                         Назвала Лиственным      Назвала Хвойным        │
│   На деле Лиственный:         14 578                 1 706             │
│   На деле Хвойный:             1 580                11 326             │
│                                                                        │
│   • Overall Accuracy (Общая точность): 88.7%                           │
│   • Каппа Коэна (согласие сверх случайного): 0.772                     │
│   • Полнота детекции хвои (Recall): 87.8%                              │
└────────────────────────────────────────────────────────────────────────┘

Честный отрицательный результат: сосна против ели

А можно ли выделить саму сосну? Мы проверили это математически.

Среди хвойных точек попытались обучить модель предсказывать песчаный субстрат по спектру кроны:

• По оптике показатель ROC AUC (Receiver Operating Characteristic — Area Under the Curve, площадь под кривой ошибок) составил ровно 0,500 — эквивалент подбрасывания монетки.

• Радар поднял AUC лишь до 0,537.

Спектрально кроны сосны и ели на песчаных и суглинистых почвах в среднем диапазоне Sentinel-2 неотличимы. А обучающих меток конкретных пород в OSM катастрофически мало: сосны набралось около 150 выделов, а ели — всего 28.

Выделить породу сосны в чистом виде из бесплатного космоса и выборки OSM пока не получилось. Надо либо доразмечать датасет в ручную либо взять область шире. Но наш рабочий суррогат - это пересечение ML-слоя хвойного леса 20 м с песчаными контурами ГГК-1000 и сухим рельефом — на практике прекрасно локализует боровые биотопы.

Проблема лишайникового покрова

Беломошник - геоботанический тип леса, а не ландшафтный класс. Напочвенный лишайниковый покров со спутников среднего разрешения сквозь сомкнутый полог крон не считывается. Попытки найти повыдельные таксационные карты в открытых источниках (ФГИС ЛК, региональные лесохозяйственные регламенты) упёрлись в то, что по Ленинградской области геометрия лесоустройства в открытом доступе отсутствует. Поэтому пересечение факторов остаётся наиболее честным приближением.

Карты, геоданные и пространственные системы

• OSM (OpenStreetMap) — краудсорсинговый картографический проект с открытыми геоданными. Источник разметки лесов и контуров болот.

• CGLS (Copernicus Global Land Service) — служба глобального мониторинга земной поверхности в рамках европейской космической программы Copernicus.

• LC100 (Land Cover 100 m) — растровый продукт земного покрова с шагом пикселя 100 м (в проекте взят слой Forest-Type за 2019 год).

• ESA WorldCover — глобальная высокодетальная карта типов поверхности с разрешением 10 м от Европейского космического агентства (ESA). В проекте использован класс Tree cover версии v200.

• Copernicus DEM GLO-30 — глобальная цифровая модель высот (DEM — Digital Elevation Model) с шагом около 30 метров, источник данных о высоте и углах наклона склонов.

• GDAL (Geospatial Data Abstraction Library) — базовая библиотека для трансляции, трансформации и чтения растровых и векторных геоформатов.

• GeoTIFF — стандартизированный формат растровых графических файлов TIFF, содержащий встроенные географические метаданные о координатной привязке.

• COG (Cloud-Optimized GeoTIFF) — оптимизированный для облачных сред GeoTIFF с внутренней тайловой структурой, позволяющий скачивать фрагменты растра HTTP-запросами диапазонов байтов без полной загрузки файла.

• GeoJSON — открытый формат описания географических объектов (точек, ломаных, полигонов и их свойств), основанный на синтаксисе JSON.

• EPSG — реестр геодезических параметров Европейской группы нефтеразведки (European Petroleum Survey Group). Числовые коды обозначают конкретные координатные системы.

• EPSG:3857 (Web Mercator / Pseudo-Mercator) — общепринятая в веб-картографии цилиндрическая проекция на сфероиде (OSM, Google Maps).

• EPSG:3395 (World Mercator) — эллипсоидальная проекция Меркатора, применяемая в тайловых сервисах Яндекса.

• LAEA (Lambert Azimuthal Equal Area) — азимутальная равновеликая проекция Ламберта. Сохраняет истинные площади контуров, устраняя масштабные искажения при расчётах в метрах.

• WGS84 (World Geodetic System 1984) — общеземная эллипсоидальная геодезическая система, глобальный стандарт спутниковой навигации GPS.

• ГГК-1000/3 — Государственная геологическая карта масштаба 1:1 000 000 третьего поколения (ВСЕГЕИ).

• — геологический индекс флювиогляциальных (водно-ледниковых) отложений (пески, гравийники зандров и озов).

• — геологический индекс озёрно-ледниковых отложений (пески, супеси, ленточные глины приледниковых озёр).

• — геологический индекс эоловых отложений (дюнные перевеянные пески).

• — геологический индекс ледниковых моренных отложений (несортированные суглинки, валунники).

Спутниковое зондирование и машинное обучение

• NDVI (Normalized Difference Vegetation Index) — нормализованный относительный вегетационный индекс, количественная мера фотосинтетической активности растительности.

• SCL (Scene Classification Layer) — аппаратно-программный растровый слой классификации сцены спутника Sentinel-2, позволяющий автоматически отфильтровывать облака, тени крон, снежный покров и битые пиксели.

• SWIR (Short-Wave Infrared) — коротковолновый инфракрасный спектральный диапазон, отражение в котором определяется структурой крон и содержанием влаги в хвое/листве.

• B11 / B12 — каналы 11 и 12 сенсора MSI спутника Sentinel-2, соответствующие полосам SWIR.

• VV / VH — поляризационные моды радара с синтезированной апертурой (SAR) спутника Sentinel-1 (VV — вертикальное излучение и приём, VH — вертикальное излучение и горизонтальный приём).

• Random Forest (Случайный лес) — ансамблевый алгоритм классификации на базе комитетного голосования множества независимых решающих деревьев.

• ROC AUC (Receiver Operating Characteristic — Area Under the Curve) — метрика площади под кривой ошибок, отражающая интегральную способность бинарного классификатора разделять классы независимо от выбранного порога (0,5 — уровень случайного шума, 1,0 — идеальное разделение).

• Каппа Коэна () — статистическая метрика согласованности предсказаний модели с эталонной разметкой за вычетом вероятности случайного совпадения классов.

• Overall Accuracy (Общая точность) — отношение числа безошибочно классифицированных тестовых точек к общему размеру проверочной выборки.

Компьютерное зрение и векторизация

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

• HSV (Hue, Saturation, Value) — цилиндрическое представление цвета по тону, насыщенности и яркости.

• RANSAC / USAC / MAGSAC++ — семейство робастных алгоритмов оценки математических моделей на зашумлённых массивах данных с высокой долей грубых ошибок (выбросов).

Сетевой стек и облачная инфраструктура

• API (Application Programming Interface) — программный интерфейс для взаимодействия сетевых сервисов и передачи структурированных данных.

• 206 Partial Content — статусный код протокола HTTP, возвращающий запрошенный по смещению диапазон байтов файла.

• KV (Key–Value) — распределённое бессерверное хранилище типа «ключ — значение» платформы Cloudflare Workers.

• R2 — объектное S3-совместимое облачное хранилище Cloudflare с нулевой стоимостью исходящего трафика.

• WebGL (Web Graphics Library) — браузерный API для высокопроизводительного рендеринга 2D- и 3D-графики с использованием аппаратных ресурсов GPU.

• PMTiles — открытый облачно-ориентированный однофайловый архив для картографических векторных и растровых тайлов.

Биология, сокращения и единицы измерений

• ITS (Internal Transcribed Spacer) — внутренний транскрибируемый спейсер ДНК, универсальный молекулярный штрихкод в систематике грибов.

• B. pinophilus / B. edulis — общепринятое в биологической номенклатуре сокращение названий видов Boletus pinophilus (белый гриб сосновый) и Boletus edulis (белый гриб сетчатый/обыкновенный).

• ООПТ — особо охраняемая природная территория.

• РСФСР — Российская Советская Федеративная Социалистическая Республика.