Как картографировать городскую запечатанность почв по космическим снимкам?

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

Автор: Юрий Дворников

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

Исследователи из научного центра SUNLAB и МГУ выбрали десять российских городов, расположенных в разных широтах на территории Европейской России, в Сибири и на Дальнем Востоке, охватывающих биоклиматические зоны от субарктической до степной (как показано на карте), для того чтобы определить, какой именно способ спектральной трансформации снимков лучше всего подходит для картографирования запечатанности почв.

Спутниковые снимки разных сезонов

Для территории каждого города были выбраны снимки (или дневные композиты), сделанные в период начала (GSO – growing season onset), середины (GSM – growing season middle) и конца (GSE – growing season end) вегетационного периода весной – осенью 2019/2020 г. Мы определили центры этих периодов по 16-дневному тренду данных MODIS NDVI 250m и трендам температуры воздуха (на высоте 2 м) по данным реанализа ERA5-Land. Начало и конец вегетационного периода определялись как среднее двух дат: 1) дата резкого увеличения/снижения значений NDVI; 2) дата пересечения порога температуры воздуха +5 °C. Периоды GSO и GSE были определены как одна неделя до и после соответствующих дат. GSM определялся как период плато значений MODIS NDVI.

Мы использовали данные ESA Sentinel-2 в виде значений спектрального отражения «на сенсоре» (TOA – top-of-atmosphere), а также атмосферно-скорректированные значения спектральной яркости (SR – surface reflectance) при помощи алгоритма sen2cor, доступный через платформу облачных вычислений Google Earth Engine. Изображения TOA использовались для получения значений спектральной трансформации TasseledCap.

Для оценки принадлежности пикселя к классу «непроницаемой поверхности» использовалось несколько подходов спектрального преобразования. Эти преобразования включали попиксельную и субпиксельную классификацию. Анализировали вегетационные индексы: модифицированный NDVI (mNDVI) и модифицированный SAVI (mSAVI). Для этих двух индексов мы протестировали последовательность порогов от 0.1 до 0.7 с шагом 0.01, всего 61. Кроме того, мы также протестировали набор опубликованных индексов застройки (NDBI²), улучшенный индекс застройки – IBI³, BCI⁴, улучшенный нормализованный индекс непроницаемой поверхности – ENDISI⁵, перпендикулярный индекс непроницаемой поверхности – PISI⁶ и индекс городского композита – UCI⁷ с пороговыми последовательностями от -0.5 до 0.4 с шагом 0.01, всего 91.

Для каждого порогового значения дополнительно производилась субпиксельная классификация с использованием созданных конечных членов семи классов: открытая почва, древесная растительность, кустарники, луга, искусственные поверхности с высоким, низким и средним альбедо. Таким образом, для каждого изображения/композита было получено 1336 оценок, то есть всего мы сделали 54 776 оценок для всех городов/сезонов. Мы верифицировали качество наших оценок на основе данных экспертного дешифрирования в пределах нескольких тестовых участков.

Реакция показателей на уровне отдельного города

Исследование показало, что IBI и PISI были чувствительны к порогу во всех случаях, а это означает, что индексы могут показать худшие результаты, если будет выбрано неподходящее пороговое значение. Лучшие метрики были получены с порогами  -0,09 для IBI и 0,04 для PISI. Небольшие сдвиги порогов PISI (ниже 0,02 и выше 0,05) приводят к быстрому снижению метрик коэффициента детерминации (R²) и среднеквадратического отклонения (RMSE). Наилучшие результаты получены на основе индексов ENDISI, UCI и mNDVI. При пороговых значениях -0,15, -0,47 и 0,35 они могут достичь максимального R², равного 0,73, 0,79 и 0,76 соответственно. Эти индексы, наряду с mSAVI являются менее чувствительными к пороговому значению. mNDVI и ENDISI показали стабильные результаты в широком диапазоне пороговых значений (0,2–0,45 для mNDVI и −0,25–0 для ENDISI).

При сравнении лучших показателей отдельно по всем городам явного «лучшего» сезона выявлено не было: в шести городах (Череповец, Москва, Пущино, Калуга, Курск, Ростов-на-Дону) лучшие показатели были получены при использовании двулетних медианных композитов (МC). Однако эти показатели также были очень близки к другим сезонам (в некоторых случаях GSM или GSE). Для Мурманска, Новосибирска, Курска и Находки самые высокие показатели получены для сцен Sentinel-2, сделанных в устоявшийся (GSM) вегетационный период. В среднем, R² для всех случаев превышал 0,85 (с максимальным значением 0,97 для Апатитов и минимальным 0,63 для Находки), а среднее значение RMSE составляло <9% (4,94 % для Апатитов – минимальное значение и 14,88 % для Москвы – максимальное значение). Все индексы, кроме BCI (биофизический композитный индекс), хотя бы один раз фигурировали в итоговой таблице лучших метрик, хотя в целом индексы застройки встречались чаще, чем вегетационные.

Выявление наиболее стабильных индексов

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

Для всех медианных композитов (MC) UCI (-0,49) и mNDVI (0,41) точно отражали запечатанность во всех городах. PISI смог показать аналогичную картину только для нескольких городов (Новосибирск, Пущино, Калуга, Находка), IBI – только для Череповца. В других городах PISI (0,01) и IBI (-0,14) имели артефакты в пределах классов открытой почвы и древесной растительности. За период GSM UCI (-0,49) и mNDVI (0,41) были лучшими для всех городов, кроме Москвы и Находки. Все четыре выбранных индекса имели близкие метрики при рассмотрении всех тестовых территорий (n=254) во исследуемых городах.

Среди глобальных наборов данных, ESA Worldcover (2020) и GISA-10m показали близкие характеристики по сравнению с UCI и mNDVI, тогда как DLR WSF2019 недооценил, а ESRI Land Cover значительно переоценил процент непроницаемых площадей на тестовых участках во всех десяти городах (n=254). Завышение показателей данных ESRI связано с неопределённостью внутригородских ландшафтов – городские газоны и зелёная инфраструктура нанесены на карту как непроницаемые. Эти завышенные оценки вносят значительный вклад в окончательные оценки запечатанности для отдельных городов. В среднем WSF2019 недооценил квазиистинную запечатанность на 22 %, ESRI Land Cover – переоценил на 164 %, а ESA Worlcover – переоценил всего на 3 %. Статистически значимых различий между оценками ESA и mNDVI/UCI выявлено не было. Набор данных GISA-10m имеет тенденцию переоценивать реальную картину, но незначительно.

Мелкоочаговые городские структуры (гаражи, жилые малоэтажные дома), очевидно, вносят большой вклад в общую ошибку запечатанности при разрешении данных 10 м. При картографировании с использованием трёх наиболее эффективных индексов (mNDVI, PISI, UCI), а также набора данных ESA Woldcover, самые высокие остатки при обнаружении запечатанности ожидаются в гаражах – для этого типа землепользования они статистически отличались от других типов землепользования во всех случаях, кроме mNDVI, который менее чувствителен к выбросам значений.

Итоги: влияние выбранных подходов и сезонов

Мы обнаружили, что среди рассмотренных индексов mNDVI (0,41) и UCI (-0,49), применённые к медианным композитам, являются наиболее подходящими для картографирования запечатанных поверхностей в городских ландшафтах, как показали метрики качества и визуальный осмотр полученных бинарных классификаций. Индексы PISI и IBI показали высокую точность определения запечатанности, но они имели большее количество артефактов, преимущественно в пределах сельскохозяйственных полей и густой древесной растительности. Тем не менее, порог PISI 0,01 показал хорошую картину запечатанности для нескольких городов (Новосибирск, Пущино, Калуга, Находка), IBI (-0,14) – для Череповца.

Тени в населённых пунктах от густой древесной растительности, сельскохозяйственные поля и водные объекты являются основными источниками ошибок при автоматизированном дешифрировании запечатанности с использованием оптических данных. Эти эффекты можно частично устранить, используя летние медианные композиты. Поскольку пиксели класса теней часто ошибочно классифицируются как водные объекты, мы применяли маску воды по данным OpenStreetMap.

Наши данные отчетливо показали преимущество данных GSM и MC над GSO и GSE при дешифрировании запечатанности. Вероятно, это связано с более контрастным спектральным откликом искусственного материала по сравнению с растительностью / открытой почвой: при картографировании запечатанности с помощью изображений, сделанных в течение установившегося вегетационного периода, существует меньше шансов, что пиксели открытой почвы будут ошибочно классифицированы как запечатанные и наоборот.

Исследование выявило, что наибольшая точность может быть достигнута путём применения пороговых значений mNDVI и UCI (0,41 и -0,49 соответственно) для летних медианных композитов данных Sentinel-2A/B (наивысшие значения R² > 0,82 и наименьшие значения RMSE < 10 %), если производить валидацию «по площади». В случаях, когда количество безоблачных снимков недостаточно, следует использовать данные за установившийся вегетационный период. Малые по размеру антропогенные объекты обладают наибольшей погрешностью при разрешении 10 м, но в целом точность их картографирования не превышает 20 %. Среди глобальных наборов данных только ESA Worldcover 10-m продемонстрировал сопоставимые показатели R² и RMSE. Более того, сравниваемые глобальные наборы данных показали значительные различия (до десятков %) между оценками запечатанной поверхности для выбранных десяти городов, что подчёркивает необходимость дальнейшей валидации этих данных. Наши результаты могут быть успешно использованы для картографирования годовой и даже сезонной динамики запечатанности городской среды.

Полный материал исследования опубликован по ссылке: https://www.sciencedirect.com/science/article/pii/S1569843223002947.

Дистанционный анализ и моделирование выполнены при поддержке проекта Российского научного фонда №19-77-300-12. Обработка данных и подготовка статьи выполнены при грантовой поддержке РУДН (грант №202194-2-000), а также в рамках государственного задания (№122111000095-8).

Ссылки

1. Karger, D.N., Conrad, O., Bohner, J., Kawohl, T., Kreft, H., Soria-Auza, R.W., Zimmermann, N.E., Linder, P., Kessler, M., 2017. Climatologies at high resolution for the Earth land surface areas. Sci. Data 4, 170122. 10.1038/sdata.2017.122. ↩︎
2. Zha, Y., Gao, J., Ni, S., 2003. Use of normalized difference built-up index in automatically mapping urban areas from TM imagery. Int. J. Remote Sens. 24, 583–594. https://doi.org/10.1080/01431160304987. ↩︎
3. Xu, H., 2008. A new index for delineating built-up land features in satellite imagery. Int. J. Remote Sens. 29, 4269–4276. https://doi.org/10.1080/01431160802039957. ↩︎
4. Deng, C., Wu, C., 2012. BCI: A biophysical composition index for remote sensing of urban environments. Remote Sens. Environ. 127, 247–259. https://doi.org/10.1016/j.rse.2012.09.009 ↩︎
5. Chen, J., Yang, K., Chen, S., Yang, C., Zhang, S., He, L., 2019. Enhanced normalized difference index for impervious surface area estimation at the plateau basin scale. J. Appl. Remote Sens. 13, 016502 https://doi.org/10.1117/1.jrs.13.016502. ↩︎
6. Tian, Y., Chen, H., Song, Q., Zheng, K., 2018. A novel index for impervious surface area mapping: Development and validation. Remote Sens. 10 https://doi.org/10.3390/rs10101521. ↩︎
7. Zhang, Y., Zhang, H., Lin, H., 2014. Improving the impervious surface estimation with combined use of optical and SAR remote sensing images. Remote Sens. Environ. 141, 155–167. https://doi.org/10.1016/j.rse.2013.10.028 ↩︎