Индексы вегетации NDVI и NDRE по снимкам Sentinel-2 и агродрона на ноутбуке: формулы, смещение L2A и фильтрация облаков

Индексы вегетации NDVI и NDRE по снимкам Sentinel-2 и агродрона на ноутбуке: формулы, смещение L2A и фильтрация облаков

Спутниковая группировка Copernicus Sentinel-2 даёт сельхозпредприятию бесплатную мультиспектральную съёмку с разрешением 10 и 20 метров на пиксель каждые 3–5 суток, а мультиспектральный агродрон детализирует опытные делянки и очаги сорняков до 2–5 см/пиксель. Локальный расчёт индексов на ноутбуке в QGIS + GDAL избавляет хозяйство от слепого доверия к «чёрному ящику» облачных подписок: агроном точно видит, не принята ли тень от кучевого облака за пятно азотного голодания пшеницы и не искажены ли значения обновлённой калибровкой спутника.

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

Спектральные каналы Sentinel-2 Level-2A: почему важен уровень атмосферной коррекции

Для сравнения состояния посевов между апрелем, маем и июнем скачивают продукт уровня Level-2A (BOA — Bottom of Atmosphere Surface Reflectance), в котором уже устранено рассеяние света атмосферой и аэрозольной дымкой (сырой уровень L1C отражает верхнюю границу атмосферы и плавает от даты к дате!).

Главная ловушка сырых файлов .SAFE с 2022–2026 годов: сдвиг -1000 (BOA_ADD_OFFSET)

Начиная с версии алгоритма обработки ESA Processing Baseline 04.00+, ко всем сырым 16-битным значениям пикселей (DN) в архивах IMG_DATA уровня Level-2A прибавлена константа +1000 (чтобы сохранить отрицательные шумы тёмных объектов в беззнаковом формате UInt16 без обрезки в ноль).

Что происходит, если агроном скачал сырые растры .jp2 из Copernicus Browser и просто вбил в «Калькулятор растров» старую школьную формулу (B08 - B04) / (B08 + B04) без вычитания тысячи? В числителе сдвиги +1000 взаимно уничтожаются, а в знаменателе складываются в лишние +2000, из-за чего все значения NDVI на весенней почве и слабых всходах искусственно занижаются на 15–35 %!

Правильный расчёт истинного коэффициента отражения поверхности (Surface Reflectance) в калькуляторе растров QGIS / gdal_calc:

# Истинное отражение каждого канала: SR = (DN - 1000.0) / 10000.0
# Корректная формула NDVI в QGIS Raster Calculator по сырым JP2 Baseline 04.00+:
( ( "B08@1" - 1000.0 ) - ( "B04@1" - 1000.0 ) ) / ( ( "B08@1" - 1000.0 ) + ( "B04@1" - 1000.0 ) )

(Всегда проверяйте инструментом «Определить объекты» один пиксель чистой воды или асфальта и XML-метаданные MTD_MSIL2A.xml: если вы скачали не сырой .SAFE, а уже гармонизированный Cloud-Optimized GeoTIFF из облачного каталога STAC, сдвиг мог быть вычтен поставщиком — тогда значения тёмной воды сразу лежат около 0,01..0,05, а не около 1050!)

Фильтрация теней от облаков и краевого эффекта лесополос до нарезки норм внесения

Полупрозрачная перистая дымка и особенно тень от кучевого облака резко роняют расчётный NDVI на совершенно здоровом участке поля на 0,10–0,18 пункта. Если пустить такой снимок напрямую в генератор карт-предписаний VRA в формате Shapefile или ISOXML, разбрасыватель высыплет максимальную дозу азота ровно по пятну вчерашней облачной тени, вызвав полегание хлебов!

Обязательная трёхступенчатая защита качества на ноутбуке:

1. Маскирование по растру SCL + визуальный контроль в естественных цветах (B04-B03-B02) и псевдоцветах SWIR (B12-B8A-B04): любые пиксели, попавшие под классы SCL 3 (тень облака), 8, 9, 10 (облачность и циррусы), переводятся в NoData, а не интерполируются как «низкая биомасса».
2. Внутренний буфер контура поля -10..-15 метров от края: крайний 10-метровый пиксель Sentinel-2 физически захватывает половину грунтовой дороги и крону лесополосы («смешанный пиксель»), завышая или занижая статистику всего поля. Вычтите внутренний буфер из эталонных контуров хозяйства (требования к самим контурам полей разобраны в статье про учёт границ в ЕФИС ЗСН и ФГИС «Зерно»).
3. Многолетняя композитная карта устойчивых зон плодородия: для базового внесения фосфора/калия и дифференцированного сева никогда не опирайтесь на один снимок одного года — усредните нормализованные растры 3–5 контрастных сезонов в фазу активной вегетации, чтобы отделить реальный рельеф и почвенные разности от разового огреха сеялки.

Сколько ресурсов ноутбука требуют космоснимки Sentinel-2 и мультиспектр БПЛА

Два сценария предъявляют совершенно разный счёт к железу ноутбука агронома:

• Спутниковая аналитика Sentinel-2 / Landsat на земельный банк 10–80 тыс. га: архив одной полной гранулы 100×100 км весит около 800–1100 МБ в сжатом JPEG2000. Распаковка JP2 на лету грузит ядра CPU и диск: скачивайте через STAC API (Cloud-Optimized GeoTIFF, COG) только прямоугольник (BBox) своих полей — тогда расчёт NDVI + NDRE + зональной статистики по 150 полям на ноутбуке с 6–8-ядерным CPU, 16–32 ГБ ОЗУ и быстрым NVMe SSD выполняется за 2–4 минуты.
• Фотограмметрическая сшивка мультиспектра с БПЛА (DJI Mavic 3 Multispectral, Geoscan и аналоги — про подключение базовых растровых слоёв в проект см. статью подложки и тайлы в QGIS): здесь за один вылет прилетают 2 000–6 000 16-битных TIFF-кадров в 4–5 узких каналах с панелью калибровки освещённости. Для сшивки ортофотоплана в Agisoft Metashape / Pix4Dfields / OpenDroneMap в полевой гостинице нужен производительный ноутбук уровня мобильной рабочей станции — 32–64 ГБ ОЗУ, NVMe SSD от 1–2 ТБ и дискретная видеокарта NVIDIA RTX с поддержкой CUDA/Vulkan с обслуженной системой отвода тепла (подробный разбор теплового режима под многочасовой нагрузкой есть в гиде по обслуживанию системы охлаждения мощного ноутбука).

Чек-лист расчёта карты вегетации перед выездом на осмотр полей

• [ ] Использован продукт Sentinel-2 Level-2A (или откалиброванный по солнечной панели мультиспектр БПЛА), проверено наличие/отсутствие сдвига -1000 в конкретном источнике данных.
• [ ] По слою SCL и визуальному осмотру полностью исключены полупрозрачная перистая облачность (Cirrus) и тени облаков.
• [ ] При сомкнутом пологе посевов (NDVI > 0,78) дополнительно рассчитан краснокраевой индекс NDRE по каналам B8A и B05.
• [ ] Из расчёта средних значений и гистограммы по полю вырезан краевой 10-метровый буфер дорог, столбов ЛЭП и лесополос.
• [ ] Намечены GPS-точки контрастных зон (минимум, среднее, максимум) для обязательного наземного осмотра с лопатой и листовой диагностикой до утверждения доз удобрений!

Частые вопросы

Означает ли низкий NDVI на пятне поля, что туда всегда нужно дать больше азота? Категорически нет! Космоснимок фиксирует лишь угнетённую надземную биомассу, но не называет её причину: в этой западине могло произойти вымокание корней, сильная кислотность почвы (pH < 4,8), поражение корневыми гнилями или очаг проволочника. Увеличивать норму удобрений на мёртвом вымокшем пятне убыточно — там применяют стратегию экономии («меньше в безнадёжную зону»), и ответ даёт только выезд агронома в точку по GPS!

Можно ли считать NDVI по фотографиям обычного бытового дрона с обычной RGB-камерой? Настоящий NDVI по обычной RGB-матрице посчитать физически невозможно, потому что у неё нет ближнего инфракрасного канала (NIR 840 нм) — там стоит ИК-отсекающий фильтр. По обычному цветному снимку считают только псевдовегетационные видимые индексы озеленения (VARI, ExG), которые годятся для подсчёта всходов подсолнечника/кукурузы и пятен сорняков по пашне, но не заменяют откалиброванный мультиспектр для расчёта доз азота.