Digital mapping of organic carbon content and stocks in soils of Cis-Salair drained plain using the Google Earth Engine online platform and the Random Forest algorithm

Full Text

Введение

Тематические карты являются основным источником получения информации о пространственном распределении содержания и запасов почвенного органического углерода (ПОУ), которые помогают улучшить наши знания о плодородии почв, последствиях изменения климата, а также в почвенно-агроэкологическом мониторинге и моделировании. Цифровое почвенное картографирование (ЦПК) базируется на идеях В.В. Докучаева и Х. Дженни о связи свойств почв с факторами почвообразования. Геопространственные данные о факторах почвообразования в ЦПК принято называть предикторами или ковариатами окружающей среды (Hengl, 2007; McBratney et al., 2003), которые необходимы для геопространственного моделирования и создаются на основе данных дистанционного зондирования Земли, архивных источников и др. Установление основных факторов, влияющих на качественные и количественные характеристики почв, дает ценную информацию для понимания их пространственного распределения (Suleymanov et al., 2024).

Особое внимание в ЦПК уделяется изучению взаимосвязи между предикторами и различными свойствами почв, необходимой для обоснования выбора того или иного метода картографирования, которые в настоящее время отличаются большим разнообразием. К ним относятся обобщенные линейные модели, деревья классификации и регрессии, нейронные и нечеткие нейронные сети, геостатистика (Hengl, 2007). Алгоритм случайного леса является перспективным и более надежным с точки зрения уменьшения шума в данных и способен обрабатывать как количественные, так и качественные наборы геопространственных данных (Breiman, 2001; Gandhi, 2024). Потенциал алгоритма случайного леса в цифровом картографировании содержания и запасов ПОУ был продемонстрирован в следующих исследованиях (Cutler et al., 2007; Dharumarajan et al., 2017; Grimm et al., 2008; FAO and ITPS, 2020; Sreenivas et al., 2016; Suleymanov et al., 2024; Vågen et al., 2016).

Несмотря на широкое распространение различных алгоритмов машинного обучения, используемых в ЦПК, в Российской Федерации количество исследований в этом направлении недостаточно. В этом контексте цель исследования – провести цифровое картографирование содержания и запасов почвенного органического углерода (ПОУ), а также плотности сложения в пахотном слое (0–30 см) почв Предсалаирской дренированной равнины с использованием набора геопространственных данных, алгоритма случайного леса, реализованного на облачной онлайн-платформе Google Earth Engine (GEE), и педотрансферной функции. В настоящей работе проведено сравнение влияния предикторов на результаты картографирования содержания ПОУ, относящихся к 4-м факторам почвообразования (климат, рельеф, растительность, пространственное положение). Коэффициенты корреляции (R) между содержанием ПОУ и значениями предикторов учтены при формировании следующих наборов предикторов: 1) BIO11+RVI; 2) Longitude+CNBL; 3) SАТ10+CC+Texture; 4) 60 предикторов; 5) 42 (без кривизн рельефа, вегетационных индексов и предикторов с нулевыми значениями); 6) 37 (все с R > ±0.5); 7) 32 (все с R > ±0.3 без вегетационных индексов); 8) 27 (все с R > ±0.5 без вегетационных индексов); 9) 23 (без BIO1–19, кривизн рельефа, вегетационных индексов и предикторов с наличием нулевых значений в растровой карте).

Объекты и методы

Исследования проведены на территории Предсалаирской дренированной равнины (Предсалаирье) в лесостепной зоне юго-восточной части Западной Сибири (Новосибирская область, Тогучинский район, с. Усть-Каменка, 55°02´20˝ с.ш., 83°50´00˝ в.д.).

Участок пахотного угодья площадью 225 га и протяженностью 4 км занимает территорию водосборных бассейнов рек Ирба и Хайрузовка (рис. 1). На исследуемой территории распространены агрочерноземы, агротемно-серые и агросерые почвы (табл. 1).

Рис. 1. Карта высот территории исследования и схема отбора почвенных проб. Примечание для рис. 1, 3, 6: черными точками выделен обучающий набор данных (n = 42), красными – валидационный (n = 15).

Таблица 1. Преобладающие почвы на территории исследования

Согласно группировке почв по содержанию гумуса (Шишов и др., 2004), агрочерноземы в верхней части склона относятся к сильногумусированным (5–8%); агрочерноземы, агротемно-серые и агросерые почвы в средней и нижней частях склона – к среднегумусированным (3–5%). Почвы характеризуются слабокислой и близкой к нейтральной реакцией среды (рН водной суспензии) и тяжелосуглинистым гранулометрическим составом.

Отбор индивидуальных почвенных проб (n = 57) произведен буром из пахотного слоя (0–30 см) по нерегулярной сетке в начале июля 2013 года. Координаты точек опробования определены с помощью системы геопозиционирования (GPS, Garmin eTrex Vista). Почвенные пробы проанализированы на содержание ПОУ методом мокрого озоления по Тюрину (Аринушкина, 1970). Данные по содержанию ПОУ умножены на коэффициент 1.15, который, согласно исследованиям (FAO, 2021), позволяет учесть ошибку, связанную с неполным окислением органического углерода, определяемого методом Тюрина. Карта плотности сложения почв рассчитана на основе растровой карты содержания ПОУ с использованием педотрансферной функции (Abdelbaki, 2018):

𝐵𝐷 = 1.44 × exp(−0.03 × ПОУ),

где 𝐵𝐷 (г/см³) – плотность сложения почв; ПОУ (%) – содержание почвенного органического углерода.

Запасы ПОУ (ЗПОУ) рассчитаны по следующей формуле:

ЗПОУ = ПОУ (%) × 𝐵𝐷 (г/см³) × МС (см),

где ПОУ (%) – содержание почвенного органического углерода; 𝐵𝐷 (г/см³) – плотность сложения почв; МС – мощность слоя (см).

В качестве предикторов в модели случайного леса использовали: 19 биоклиматических переменных WorldClim (Fick, Hijmans, 2017); 5 климатических переменных, рассчитанных на основе WorldClim и данных и почвенно-климатического атласа (Сляднев, 1978); 8 вегетационных индексов, рассчитанных на основе Landsat 8 OLI (снимок от 7 июля 2013 г.); 26 морфометрических характеристик рельефа, рассчитанных на основе цифровой модели высот ALOS World 3D ver. 4.1; 2 переменные, характеризующие пространственное положение (табл. 2). В общем было использовано 60 предикторов, характеризующих такие факторы почвообразования как климат, растительность, рельеф и пространственное положение. Все растровые карты приведены к разрешению 30 × 30 м.

Таблица 2. Предикторы, используемые для моделирования содержания ПОУ

Обучающий набор данных (ОНД) состоял из 42 точек лабораторно-полевого обследования почв, валидационный (ВНД) – из 15 (рис. 1). Описательная статистика наборов данных по содержанию ПОУ приведена в таблице 3.

Таблица 3. Описательная статистика по содержанию ПОУ для общего, обучающего и валидационного наборов данных

Сокращения: SD – стандартное отклонение; CV – коэффициент вариации.

Оценку эффективности моделей проводили по следующим показателям: R² – коэффициент детерминации; RMSE (англ. Root Mean Square Error) – корень из среднеквадратической ошибки; MAPE (англ. Mean Absolute Percentage Error) – средняя абсолютная ошибка в процентах, MAE (англ. Mean Absolute Error) – средняя абсолютная ошибка.

Корень из среднеквадратической ошибки (RMSE, Root Mean Square Error) рассчитан следующим образом (Hengl, 2007):

$\text{RMSE}=\sqrt{\frac{1}{n}\sum _{j=0}^n{\left[\widehat{z}\left(s_j\right)-z\left(s_j\right)\right]}^2}$,

где n – количество наблюдений в анализируемом наборе данных; $\widehat{z}\left(s_j\right)$ – предсказанные значения; z(sj) – фактические значения в точках проверки. RMSE представлен в единицах измерения почвенного показателя. Значения RMSE, близкие к нулю, указывают на идеальное соответствие модели фактическим данным. Значение RMSE, меньше половины стандартного отклонения (SD) измеренных данных, указывает на высокую точность прогноза (Singh et al., 2004).

Средняя абсолютная ошибка в процентах (MAPE, Mean Absolute Percentage Error) рассчитана следующим образом (Афанасьев, Цыпин, 2008):

$\text{MAPE}(\%)=\frac{1}{n}\sum \frac{|y_t-\overline{y_t}|}{y_t}\times 100$,

где y_t – фактическое значение, $\overline{y_t}$– предсказанное значение, n – количество наблюдений в анализируемом наборе данных. $|y_t-\overline{y_t}|$ – показывает абсолютное значение (берется без отрицательного знака).

MAPE показывает отклонение предсказанного значения от фактического в среднем по выборке. Точность моделирования оценивается по значениям MAPE следующим образом (Афанасьев, Цыпин, 2008): <10% – высокая; 10–20% – хорошая; 20–50% – удовлетворительная; и >50% – неудовлетворительная точность.

Средняя абсолютная ошибка (МАЕ, Mean Absolute Error) рассчитана по следующей формуле:

$\text{MAE}=\frac{1}{n}\sum _{i=1}^n|y_t-\overline{y_t}|$,

где y_t – фактическое значение, $\overline{y_t}$ – предсказанное значение, n – количество наблюдений в анализируемом наборе данных. $|y_t-\overline{y_t}|$ – показывает абсолютное значение (берется без отрицательного знака). MAE представлена в единицах измерения исследуемого показателя.

Математические операции с геопривязанными растровыми картами проведены на облачной онлайн-платформе Google Earth Engine (Gorelick et al., 2017) и в программном обеспечении SAGA GIS (Conrad et al., 2015). Карты содержания ПОУ составлены с использованием алгоритма случайного леса (англ. Random Forest) на облачной онлайн-платформе Google Earth Engine (200 деревьев, значения гиперпараметров по умолчанию). Карта запасов ПОУ, а также карта плотности сложения почв составлены в SAGA GIS.

Результаты и обсуждение

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

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

Таблица 4. Коэффициенты корреляции (R) Спирмена между содержанием почвенного органического углерода (ПОУ) и значениями предикторов, характеризующих факторы почвообразования (p < 0.05)

Сокращения в обозначении показателей приведены в таблице 1.

Положительные корреляции заметной (R от 0.5 до 0.7) и высокой (R от 0.7 до 0.9) силы связи выявлены между содержанием ПОУ и следующими предикторами: BIO1, BIO12, BIO14, BIO16–BIO18, SАТ10, CH, AgPot, Elevation, Convexity, Texture, RSP, CND, CNBL, Longitude, Latitude, NDVI, CTVI, DVI, NRVI, RVI, SAVI, TTVI, TVI. Отрицательные корреляции заметной (R от -0.5 до -0.7) и высокой (R от -0.7 до -0.9) силы связи выявлены между содержанием ПОУ и следующими предикторами: BIO4–BIO11, BIO15, SradMean, CC, Evaporation, ValDepth. Регрессионный анализ методом пошагового введения переменных с использованием 35-ти предикторов показал, что наибольший коэффициент детерминации (R² = 0.75) получен при использовании в модели всего двух переменных (BIO11 и RVI). Остальные переменные исключены в связи с эффектом мультиколлинеарности (высокая сила связи между предикторами) или статистической незначимостью коэффициентов в уравнениях регрессии. После исключения переменных BIO11 и RVI из анализа с использованием аналогичного подхода выбраны следующие наборы предикторов: Longitude+CNBL (R²=0.65); SАТ10+CC+Texture (R² = 0.69).

При построении модели случайного леса нет ограничения на количество используемых предикторов, поэтому для сравнительного анализа подготовили еще 6 наборов предикторов: 60 предикторов; 42 (без кривизн рельефа, вегетационных индексов и предикторов с нулевыми значениями); 37 (все с R > ±0.5); 32 (все с R > ±0.3 без вегетационных индексов); 27 (все с R > ±0.5 без вегетационных индексов); 23 (без BIO1–19, кривизн рельефа, вегетационных индексов и предикторов с нулевыми значениями) (табл. 2, рис. 2).

Рис. 2. Показатель важности (IncNodePurity) отдельных переменных в модели случайного леса с наборами предикторов: а – BIO11+RVI; б – Longitude+CNBL; в – SАТ10+CC+Texture; г – 60 предикторов; д – 42 (без кривизн рельефа, вегетационных индексов и предикторов с нулевыми значениями); е – 37 (все с R > ±0.5); ж – 32 (все с R > ±0.3 без вегетационных индексов); з – 27 (все с R > ±0.5 без вегетационных индексов); и – 23 (без BIO1–19, кривизн рельефа, вегетационных индексов и предикторов с нулевыми значениями).

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

Результаты моделирования алгоритмом случайного леса с использованием различных наборов предикторов представлены на рис. 3.

Рис. 3. Карты пространственной изменчивости содержания ПОУ, составленные алгоритмом случайного леса на основе различных наборов предикторов: а – BIO11+RVI; б – Longitude+CNBL; в – SАТ10+CC+Texture; г – 60 предикторов; д – 42 (без кривизн рельефа, вегетационных индексов и предикторов с нулевыми значениями); е – 37 (все с R > ±0.5); ж – 32 (все с R > ±0.3 без вегетационных индексов); з – 27 (все с R > ±0.5 без вегетационных индексов); и – 23 (без BIO1–19, кривизн рельефа, вегетационных индексов и предикторов с наличием нулевых значений в растре).

Оценка точности моделирования содержания ПОУ показала, что у моделей с количеством предикторов 37 и 60 коэффициент детерминации (R² = 0.88) выше по сравнению с другими моделями, однако они не выбраны в качестве лучших в связи с более низкими коэффициентами детерминации, рассчитанными по валидационному набору данных (ВНД), и более высокими значениями RMSE (табл. 5, рис. 4).

Таблица 5. Показатели эффективности модели (ЭМ) случайного леса для цифрового картографирования содержания ПОУ

Сокращения: R² – коэффициент детерминации (рассчитано по ОНД); RMSE – корень из среднеквадратической ошибки; MAPE – средняя абсолютная ошибка в процентах; МАЕ – средняя абсолютная ошибка. RMSE, MAPE и MAE рассчитаны по ВНД.

Рис. 4. Зависимость между фактическим и предсказанным содержанием ПОУ в валидационном наборе данных (n = 15) для наборов предикторов: а – BIO11+RVI; б – Longitude+CNBL; в – SАТ10+CC+Texture; г – 60 предикторов; д – 42 (без кривизн рельефа, вегетационных индексов и предикторов с нулевыми значениями); е – 37 (все с R > ±0.5); ж – 32 (все с R > ±0.3 без вегетационных индексов); з – 27 (все с R > ±0.5 без вегетационных индексов); и – 23 (без BIO1–19, кривизн рельефа, вегетационных индексов и предикторов с нулевыми значениями).

Самый высокий R² (0.67) по ВНД отмечен у модели с 32 предикторами без вегетационных индексов (рис. 4ж). По всей видимости, использование вегетационных индексов приводит к снижению R² в моделях с 37 и 60 предикторами (рис. 4e, г).

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

Различие средних значений между предсказанным и фактическим содержанием ПОУ в обучающем и валидационном наборах данных статистически незначимо (рис. 5а, б).

Рис. 5. Сравнение средних значений ПОУ (фактических и предсказанных) по обучающему (а) и валидационному (б) наборам данных (карта содержания ПОУ на основе 32 предикторов).

Сравнение фактического содержания ПОУ с предсказанным моделью случайного леса на различных высотных ступенях показало удовлетворительное совпадение по средним значениям (рис. 6а, б). Предсказанные значения характеризовались меньшим диапазоном варьирования. В почвах первой высотной ступени содержание ПОУ было больше в 1.4–2.1 раза по сравнению с остальными (рис. 6а). Наибольшие значения содержания ПОУ выявлены в почвах верхней части склона с постепенным их снижением вниз по склону (рис. 6).

Рис. 6. Фактическое содержание ПОУ (а) и предсказанное моделью случайного леса (б) на различных высотных ступенях исследуемой территории (рис. 1). Примечание: объем выборки указан сверху над диаграммами размаха.

Карта содержания ПОУ, составленная с использованием 32 предикторов, с более высоким коэффициентом детерминации (0.67) по ВНД и более низкими значениями RMSE (0.72) и MAE (0.57) выбрана в качестве лучшей для дальнейшего моделирования запасов ПОУ (табл. 4, рис. 3ж, рис. 4ж). Согласно MAPE (22.1%), точность моделирования содержания ПОУ удовлетворительная (табл. 3).

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

Рис. 7. Карты пространственной изменчивости: а – плотности сложения почв в слое 0–30 см; б – запасов ПОУ в слое 0–30 см.

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

Таблица 6. Фактическая и предсказанная плотность сложения агротемно-серой типичной почвы разной степени смытости

RMSE составил 0.067 г/см³, а MAPE – 5%, что соответствует высокой точности моделирования плотности сложения почв. Использованная в настоящем исследовании педостранферная функция была разработана на основе географической базы данных обследований почв США (Abdelbaki, 2018). Однако расчет ошибки моделирования показал, что применение этой функции для почв исследуемой территории можно считать приемлемым.

На основе карт содержания ПОУ и плотности сложения почв проведено картографирование запасов ПОУ (рис. 7б), закономерности пространственного варьирования которого были аналогичны содержанию ПОУ (рис. 8а, рис. 6а). Значения по плотности сложения почв (предсказанные) варьировали в диапазоне от 1.22 до 1.35 г/см³ и увеличивались вниз по склону, т.е. имели обратный тренд распределения по сравнению с содержанием и запасами ПОУ (рис. 8б).

Рис. 8. Запасы ПОУ (а) и плотность сложения почв (б) на разных высотных ступенях исследуемой территории (см. рис. 1). Примечание: объем выборки указан сверху над диаграммами размаха.

Растровая карта исследуемой территории площадью 225 га состоит из 2768 пикселей размером 30 × 30 м с определенным значением запасов ПОУ в каждом элементе в расчете на гектар. Расчеты показали, что общие запасы ПОУ в пахотном слое (0–30 см) почв исследуемой территории площадью 225 га составили 28721.14 т или 28.72 кт.

Заключение

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

Содержание ПОУ в слое (0–30 см) варьировало от 1.3 до 6.1% по фактическим данным, а запасы ПОУ – от 84 до 203 т/га. Наибольшее содержание и запасы ПОУ выявлены в почвах верхней части склона (первая высотная ступень (280–300 м)), тогда как вниз по склону отмечено постепенное снижение значений этих показателей в 1.4–2.1 раза. Плотность сложения почв (предсказанные значения) варьировала в диапазоне от 1.20 до 1.36 г/см³ и увеличивалась вниз по склону, т.е. имела обратный тренд распределения по сравнению с содержанием и запасами ПОУ. Общие запасы ПОУ в пахотном слое (0–30 см) почв исследуемой территории площадью 225 га составили 28.7 кт.

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

Благодарности

Авторы выражают благодарность к.б.н. Савенкову О.А. и к.б.н. Смирновой Н.В. за помощь в проведении полевых работ, Черепахиной Л.Д. за определение содержания ПОУ в почвенных образцах.

Источник финансирования

Работа выполнена по государственному заданию ИПА СО РАН при финансовой поддержке Министерства науки и высшего образования Российской Федерации.

About the authors

N. V. Gopp

Author for correspondence.
Email: gopp@issa-siberia.ru
Russian Federation, Novosibirsk

T. V. Nechaeva

Email: gopp@issa-siberia.ru
Russian Federation, Novosibirsk

References

Supplementary files

Supplementary Files

Action

1. JATS XML

Download

2. Fig. 1. Elevation map of the study area and soil sampling scheme. Note for Figs. 1, 3, 6: the training data set (n = 42) is highlighted in black, and the validation data set (n = 15) in red.

Download (71KB)

Indexing metadata

3. Fig. 2. The importance index (IncNodePurity) of individual variables in the random forest model with predictor sets: a – BIO11+RVI; b – Longitude+CNBL; c – SАТ10+CC+Texture; g – 60 predictors; d – 42 (excluding relief curvatures, vegetation indices, and predictors with zero values); f – 37 (all with R > ±0.5); g – 32 (all with R > ±0.3 without vegetation indices); h – 27 (all with R > ±0.5 without vegetation indices). and – 23 (excluding BIO1–19, relief curvatures, vegetation indices, and predictors with zero values).

Download (151KB)

Indexing metadata

4. Fig. 3. Maps of spatial variability of SOC content compiled by the random forest algorithm based on different sets of predictors: a – BIO11+RVI; b – Longitude+CNBL; c – SАТ10+CC+Texture; g – 60 predictors; d – 42 (excluding relief curvatures, vegetation indices, and predictors with zero values); f – 37 (all with R > ±0.5); g – 32 (all with R > ±0.3 without vegetation indices); h – 27 (all with R > ±0.5 without vegetation indices). and – 23 (excluding BIO1–19, terrain curvatures, vegetation indices, and predictors with zero values in the raster).

Download (88KB)

Indexing metadata

5. Fig. 4. Relationship between actual and predicted SOC content in the validation dataset (n = 15) for predictor sets: a – BIO11+RVI; b – Longitude+CNBL; c – SАТ10+CC+Texture; g – 60 predictors; d – 42 (excluding terrain curvatures, vegetation indices, and predictors with zero values); f – 37 (all with R > ±0.5); g – 32 (all with R > ±0.3 without vegetation indices); h – 27 (all with R > ±0.5 without vegetation indices). and – 23 (excluding BIO1–19, terrain curvatures, vegetation indices, and predictors with zero values).

Download (64KB)

Indexing metadata

6. Fig. 5. Comparison of average SOC values (actual and predicted) for the training (a) and validation (b) datasets (SOC content map based on 32 predictors).

Download (30KB)

Indexing metadata

7. Fig. 6. Actual SOC content (a) and predicted by the random forest model (b) at different elevation levels of the study area (Fig. 1). Note: the sample size is indicated above the box plots.

Download (48KB)

Indexing metadata

8. Fig. 7. Spatial variability maps: a – soil bulk density in the 0–30 cm layer; b – SOC stocks in the 0–30 cm layer.

Download (40KB)

Indexing metadata