Исследование Земли из Космоса, 2022, № 1, стр. 68-79

Радарные данные Sentinel 1

В работе использованы находящиеся в открытом доступе Sentinel 1 радарные данные С-диапазона IW (interferometric wide swath) моды с поляризацией VV и VH и пространственным разрешением 10 м и временным разрешением 12 дней. Угол обзора составляет 42.66°–42.81°. Количество рассмотренных сеансов съемки равно 18 с 14.10.2019 по 5.5.2020, что определялось временем нахождения температурных датчиков в почве. Работа с изображениями S1 осуществлялась с помощью программы SNAP (https://sentinel.esa.int/web/sentinel/toolboxes/sentinel-1). Предварительная обработка данных включала выделение фрагмента с исследуемой областью и радиометрическую калибровку.

На рис. 2 приведены графики изменения коэффициента обратного рассеяния (КОР) для обеих поляризаций за период 14.10.2019 по 5.5.2020 для 6 тестовых площадок. Отметим особенности поведения графиков: для кросс-поляризации характерна дифференциация значений КОР для разных тестовых площадок, а именно, в порядке убывания значений КОР: лес, целина, серые лесные почвы, черноземы. Для согласованной поляризации максимумы КОР (лес) и минимумы КОР (черноземы) сохраняются (за исключением 2–3 сеансов из 18), однако четкой дифференциации значений КОР между серыми лесными почвами и целиной нет.

Корреляция радарных данных с температурой воздуха, температурой поверхности почвы и температурой почвы на глубине 5 см представлена на примере тестовых площадок №№ 2 и 3 (табл. 2, размер выборки N = 18).

Получилось, что для целины (площадка № 3) корреляция выше для всех трех случаев. Построим регрессионное соотношение для случая с наибольшим коэффициентом корреляции Спирмена ${{\rho }_{S}}$ = 0.81 (целина) между КОР кросс-поляризации и температурой воздуха (рис. 3, а) и температурой почвы на глубине 5 см (рис. 3, б).

Мультиспектральные данные Sentinel 2

В работе использованы мультиспектральные данные Sentinel 2 с высоким временным, пространственным и спектральным разрешением. Периодичность мультиспектральной съемки каждым спутником составляет 10 дней, а при работе двух аппаратов – 5 дней. В мультиспектральной камере 13 каналов с разным пространственным разрешением от 10 до 60 м. Использовались данные съемочной системы S2 с уровнем обработки L2A в виде альбедо на нижней границе атмосферы (с атмосферной коррекцией). Работа с изображениями S2 осуществлялась программой SNAP. Использованы безоблачные сеансы съемки за апрель–май 2019 и 2020 гг.

Выбор исходных изображений S2

Для спутниковой оценки содержания органического вещества и глины в почве накладываются ограничения на выбор мультиспектральных изображений, связанные с выделением на изображениях сухих почв без растительности (оголенных почв). В работе (Dematte et al., 2018) приведены такие условия, которые для каналов S2 представлены следующим образом (Castaldi et al., 2019): 1) нулевая облачность в районе исследования, 2) значение вегетационного индекса NDVI = = (В8 – В4)/(В8 + В4) < 0.35 для исключения зеленой растительности, 3) разность коэффициентов отражения (КО) между каналами В3 и В2 и каналами В4 и В3 должна быть больше 0 (использование данных фильтров улучшает выделение почвы (Dematte et al., 2018)), 4) значение NBR = = (В11 – В12)/(В11 + В12) должно быть NBR < = 0.05, что позволяет выделить на изображении пикселы с сухой оголенной почвой. Почвенная влага усиливает поглощение света, и КО резко снижается. Спектральные каналы В11 и В12 сильно коррелируют с влажностью почвы (Musick, Pelletier, 1988), и их разность позволяет различать спектры сухой, влажной почвы, а также спектры, связанные с растительностью. Выбор порога для NBR сильно влияет на количество пикселов на изображении, удовлетворяющих условиям сухой почвы. Увеличение порога NBR до 0.15 приводит к уменьшению числа “нужных” пикселов в 2 раза (Castaldi et al., 2019), т.е. увеличение NBR приводит к ухудшению точности моделей для определения параметров почвы по спутниковым данным.

В результате для исследуемой территории были выбраны оптические снимки S2 L2A за 23.4.2019 и 22.4.2020, для которых выполнено условие отсутствия облачности, NDVI меняется от 0.17 до 0.24, разность КО между каналами В3 и В2 и каналами В4 и В3 больше 0, значения NBR несколько превышают порог 0.05, а именно, NBR меняется от 0.1 до 0.158. Приведенные значения выполняются для всех площадок, кроме площадки № 6 (лес), для которой значение NDVI > 0.5 и NBR > 0.17.

Оценка содержания гумуса в верхнем слое почвы по данным S2

Потери органического углерода почвы (SOC- soil organic carbon) являются одной из главных причин деградации пахотных земель. Таким образом, пространственный и временной мониторинг SOC является чрезвычайно важной задачей, решению которой посвящены многочисленные работы (Орлов и др., 2001; Караванова, Орлов. 1996; Castaldi et al., 2019; Gholizadeh et al., 2018 и многие другие). Данными авторами получены регрессионные модели с экспоненциальной связью содержания SOC и коэффициента отражения (КО), и разной степенью отрицательной корреляции между содержанием SOC и значениями КО в зависимости от спектрального канала.

Для создания регрессионной модели связи КО с Н необходимы как спутниковые, так и наземные данные, причем количество образцов почв должно быть не меньше 20 для определения корреляции. Кроме того, исследуемые тестовые участки должны иметь близкий гранулометрический состав. Дело в том, что чем более мелкие частицы почвы, тем больше КО от этих почв (Караванова, 2003). То есть модели для определения Н являются локальными, и использование имеющихся в литературе моделей требует корректировки для исследуемой территории.

В работах (Castaldi et al., 2019; Gholizadeh et al., 2018) показано, что наилучшую корреляцию между содержанием гумуса в почве и КО имеют спектральные каналы S2 В4–В6 и В11, В12. Для качественной оценки изменений содержания гумуса в почвах исследуемых тестовых участков приведем графики разности значений КО спектральных каналов S2 B4–B6, B11, B12 для сеансов съемки 23.4.2019 и 22.4.2020 (рис. 4) для каждой тестовой площадки, и оценим изменения КО за год.

Учитывая тот факт, что КО данных каналов имеет отрицательную корреляцию с содержанием гумуса в почве, оценим по знаку разности в какую сторону произошло изменение за год содержания гумуса. Разность значений КО положительна для площадки № 1 (чернозем, несмытая пашня) для всех каналов, кроме В12. Это свидетельствует о том, что значения КО для площадки № 1 увеличились в 2020 г. в сравнении с 2019 г., т.е. содержание Н уменьшилось. Для чернозема смытого (площадка № 2) есть уменьшение разности для четырех каналов и увеличение разности для одного. Для площадки № 3 (целина) разность КО для всех каналов отрицательная, что может свидетельствовать об увеличении содержания Н. Аналогичная ситуация и для площадки № 4 (серые лесные почвы несмытые). Для слабосмытых серых почв (площадка № 5) для трех каналов есть положительная разность, для двух – отрицательная.

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

Покровским в 1927 г. впервые было предложено уравнение экспоненты, определяющее связь содержания гумуса и значений КО (по Орлов и др., 2001):

где ${{\rho }_{{750}}}$ – КО на длине волны 750 нм, ${{\rho }_{{750,~h}}}$ – КО многогумусной почвы, ${{\rho }_{{750,~0}}} = \left( {{{\rho }_{{750,h}}} + А} \right)$ – КО безгумусной почвы, H – содержание гумуса, k – коэффициент, определяющий крутизну экспоненциального графика.

При использовании экспоненциальной зависимости КО от содержания гумуса в сухой оголенной почве следует учитывать три момента: 1) при небольшом содержании гумуса (до 3%) его определение наиболее точное, но разброс значений КО, соответствующих данному содержанию гумуса, велик, 2) при значениях содержания гумуса 6–7% и более определение содержания гумуса по КО неточное, влияние оказывают при одинаковом содержании гумуса влажность почвы, условия обработки (Караванова, 2003), 3) решающим фактором является выбор дат спутниковой съемки, позволяющих выявить пикселы с сухой оголенной почвой.

В настоящей работе выяснить корреляционную связь КО с содержанием гумуса на 5 тестовых участках Новосибирской области не представляется возможным из-за недостаточного числа наземных точек измерений. Судя по имеющимся наземным измерениям содержания гумуса в почвах тестовых участков (табл. 1), предполагаем, что кривые зависимости КО от Н (гумус) представлены выполаживающими участками, где точность определения содержания Н от КО невелика.

В данной работе для определения содержания Н по спутниковым данным используется экспоненциальная модель с параметрами для черноземов и серых лесных почв, представленная в (Караванова, Орлов, 1996).

Отметим, что для соответствия этой модели надо взять данные оптики для длины волны 750 нм. Для S2 – это канал В6 с длиной волны 740 нм, шириной полосы 15 нм и пространственным разрешением 20 м. Для использования параметров уравнения экспоненты по (Караванова, Орлов, 1996) в локальных условиях исследуемых тестовых участков необходима корректировка этих параметров. Одним из возможных вариантов для оценки содержания гумуса в почвах тестовых участков является использование для черноземов следующих параметров: ${{\rho }_{{750,~h}}}$ = 8.0, А = 29.1 и k = 0.1256. Для серых лесных почв ${{\rho }_{{750,~h}}}$ = 8.5, А = 40.5 и k = 0.28. В табл. 3 приведены вычисленные значения содержания гумуса Н для пяти тестовых площадок при данных коэффициентах уравнения экспоненты.

Отличие полученных значений процентного содержания гумуса Н по уравнению (табл. 3) от значений, полученных в лабораторных условиях (табл. 1), составляет по абсолютной величине не более 1.5% для 2019 г. Источники погрешностей лежат в первую очередь в недостаточном количестве наземных данных. Для каждого типа почв (черноземы, серые лесные почвы) необходимо порядка 20 и более наземных измерений, одна часть которых используется для получения параметров экспоненты и вторая часть для валидации полученного уравнения. Чем больше таких наземных данных, тем надежнее формула. Вторым моментом, приводящим к погрешностям, является выбор исходного оптического снимка, для которого наиболее важным параметром является значение спектрального индекса NBR. Увеличение значения NBR > 0.05 значительно сокращает число пикселов, удовлетворяющих условию сухой оголенной почвы. И тем не менее некоторые выводы о количественном содержании гумуса в почве сделать можно, прежде всего, об изменении значений Н за год (2019–2020 гг.). Наибольшее отличие получено для площадки № 1 – уменьшение содержания гумуса за год на 2.4%. Для остальных площадок изменения за год незначительные от 0.3 до 0.9% (см. табл. 3).

Оценка содержания глины в верхнем слое почвы по оптическим данным спутника S2

Многочисленные работы посвящены оценке компонентов текстуры почвы с использованием дистанционных оптических данных (Украинский, Чепелев, 2011; Shabou et al., 2015; Bousbih et al., 2019; Vaudouret al., 2019; Gholizadehet al., 2018 и другие). Так в работе (Gholizadehet al., 2018) показано, что наилучшей корреляцией с содержанием глины в почве обладает коэффициент отражения канала В7 спутника S2, а также спектральные индексы V = В8/В4 (Vegetation index), SAVI = 1.5 × (В8 – В4)/(В8 – В4 + 0.5) (Soil Adjusted Vegetation Index) и другие. Отмечена слабая корреляция S2 с илом и песком. Авторы (Bousbih et al., 2019) получили, что каналы В11 (SWIR1) и В12 (SWIR2) S2 наиболее чувствительны к изменениям содержания глины в почве, и обладают отрицательной корреляцией с этим содержанием. Причем, для использования этих каналов условием для получения количественной оценки содержания глины является сухая оголенная почва. В работе (Hengl, 2007) введен спектральный индекс глины CI = B11/B12 (Clay Index), который отличается сильной отрицательной корреляцией с содержанием глины в почве. Следует отметить, что количественная оценка содержания глины в почве проводится локально для исследуемой территории, и точность оценки напрямую связана с числом тестовых измерений образцов почвы.

На рис. 5, а приведены графики значений CI = = B11/B12 для двух сеансов съемки S2 23.04.2019 и 22.04.2020, что позволяет сделать качественное сравнение изменений за год в содержании глины в верхнем слое почвы исследуемых тестовых площадок № 1–№ 5. Изменения за год по содержанию глины не коснулись площадки № 5, для площадок № 2–№ 4 есть небольшое увеличение содержания глины, и небольшое уменьшение для площадки № 1.

Количественной оценкой процентного содержания глины в почве может стать регрессионная модель (Castaldi et al., 2019) экспоненциальной связи содержания в почве органического углерода с КО. Для случая с глиной в качестве переменной в экспоненте используем индекс глины CI = = В11/В12, имеющий отрицательную корреляцию с содержанием глины в почве. Формула для количественной оценки процентного содержания глины в почве такая: Clay(%) = 802 × exp(–2.69CI) – для черноземов, Clay(%) = 5123.6 × exp(–4.29CI) – для серых лесных почв. Часть наземных измерений содержания глины в почве тестовых участков Новосибирской области была использована для получения параметров экспоненты, а оставшаяся часть – для валидации. Полученные формулы носят локальный характер, и точность полученных коэффициентов в них напрямую зависит от числа наземных измерений.

Проверка на наличие корреляции CI с элементами текстуры почвы показала, что отрицательная корреляция CI с илом и физической глиной отсутствует. В формулах Clay – это частицы с диаметром больше 0.001 мм и меньше 0.01 мм.

На рис. 5, б приведены графики значений процентного содержания глины в почвах пяти тестовых участков (измеренные значения) и подсчитанные по формулам выше для 2020 г.

Качественная оценка содержания глины в почвах тестовых участков по значению индекса глины не противоречит количественной оценке по регрессионной модели.

Оценка влажности в верхнем слое почвы по оптическим и радарным данным S1 и S2

Рассмотрим возможность оценки влажности почвы по оптическим данным, радарным данным, а также по совместному использованию оптических и радарных данных.

Оптические данные. В работе (Burapapol, Nagasawa, 2016) используется индекс NDDI (normalized difference drought index) для оценки влажности почвы:

где ${\text{NDWI}} = \frac{{{\text{B8A}} - {\text{B}}11}}{{{\text{B8A}} + {\text{B}}11}}$ или ${\text{NDWI}} = \frac{{{\text{B8A}} - {\text{B}}12}}{{{\text{B8A}} + {\text{B}}12}}$ (normalized difference wet index), В8А, В11 и В12 – спектральные каналы S2. Авторами показано, что бóльшие значения NDDI соответствуют меньшим значениям влажности почвы. В данной работе для вычисления NDWI использовался канал В12. График NDDI для исследуемых площадок с черноземами (площадки № 1, № 2, № 3 (целина)) и серыми лесными почвами (№ 4, № 5, № 6 (лес)) показан на рис. 6 для сеансов оптической съемки за апрель (20, 27 и 30) и май (10, 12 и 17) 2020 г.

Выявилась сильная дифференциация значений NDDI для разных площадок для сеанса съемки 20.4.2020 и группировка тестовых участков по близким значениям влажности почвы. Самой большой влажностью обладает почва площадки № 6 (лес). Самой малой влажностью отличаются площадки №№ 4, 5 (серые лесные почвы, пашня) и № 3 (черноземы, целина). Между этими крайними значениями группируются площадки №№ 1 и 2 (черноземы несмытые и смытые, пашня). При переходе от апреля к маю дифференциация площадок по значениям NDDE значительно уменьшается, приводя к близости значений влажности почвы для всех тестовых площадок.

Радарные данные. В работе (N. Baghdadi et al., 2016) предложена полу- эмпирическая модель для определения влажности оголенной почвы, основанная на модели Dubois (Dubois et al., 1995). Параметры модели получены на базе большого числа измерений для сенсоров Х, С и L диапазонов длин волн с углами обзора 20°–45°. Авторы (Baghdadi et al., 2016) демонстрируют более точную оценку влажности почвы. Приведем формулы этой модели для поляризаций VV и VH, что соответствует поляризациям моды IW GRD Sentinel 1:

где $\sigma _{{{\text{VV}}}}^{0}$ и $\sigma _{{{\text{VH}}}}^{0}$ – коэффициенты обратного рассеяния для поляризаций VV и VH, θ – угол обзора, m_v – объемная влажность почвы в %, f(h) – функция, зависящая от шероховатости поверхности почвы. Напомним, что для применимости модели Dubois (Dubois et al., 1995) необходимо выполнение ряда условий, а именно, шероховатость поверхности kh ≤ 2.5, m_v ≤ 35%, θ ≥ 30°, $k = {{2\pi } \mathord{\left/ {\vphantom {{2\pi } \lambda }} \right. \kern-0em} \lambda }$.

Из формул (3), (4) для определения m_v получаем следующее выражение:

где $A = {{10}^{{1.15~}}}~{{\left( {\cos \theta } \right)}^{{0.6794}}}$, $B = {{{{{(\sigma _{{{\text{VV}}}}^{0})}}^{{0.44}}}} \mathord{\left/ {\vphantom {{{{{(\sigma _{{{\text{VV}}}}^{0})}}^{{0.44}}}} {{{{(\sigma _{{{\text{VH}}}}^{0})}}^{{0.71}}}}}} \right. \kern-0em} {{{{(\sigma _{{{\text{VH}}}}^{0})}}^{{0.71}}}}}$, C = 0.00429 ctgθ, (6)

Для угла обзора 43° значение А = 11.4204 (lg10А = = 1.05768), C = 0.0046.

Из полученных численных значений влажности почвы можно использовать только значения ≤35%.

Совместное использование радарных и оптических данных. Авторами (Y. Bao et al., 2018) предложен новый метод восстановления поверхностной влажности почвы SSM (surface soil moisture) в условиях частичной растительности на основе совместного использования радарных данных Sentinel 1 и оптических данных Landsat 8 OLI и модели водяного облака (water cloud model). Авторами показано, что 1) для поверхности, покрытой растительностью, NDWI = (B8A – B11)/(B8A + B11) наиболее подходит для устранения влияния растительности при оценке влажности почвы, 2) данные Sentinel 1 VV поляризации более подходят для восстановления SSM в сравнении с HV поляризацией в том числе из-за большей точности, 3) модель позволяет получать SSM с высокой точностью, о чем свидетельствует коэффициент корреляции R = 0.911 со среднеквадратичной ошибкой 0.053 см³/см³ между измеренной и полученной по модели SSM. С учетом положений 1)–3) формула для определения объемной влажности почвы SSM в слое 0–5 см такова:

где SSM измеряется в см³/см³, $\sigma _{{{\text{VV}}}}^{0}$ в дБ. Следует отметить (Bao et al., 2018), что модель (7) игнорирует влияние шероховатости поверхности и типа растительности.

В данной работе для вычисления NDWI в (7) использовался канал В12, как и выше при вычислении NDDI. Для вычисления SSM использовались радарные данные за 23 апреля и 5 мая 2020 г., а оптические, соответственно, за 22 апреля и 10 мая 2020 г. В последнем случае из-за отсутствия безоблачных сеансов с 30 апреля по 10 мая. На рис. 7 показаны графики изменения SSM (см³/см³) для двух соседних дат радарной съемки 23 апреля и 5 мая 2020 г. Отметим отрицательное значение SSM для площадки № 1, которое является следствием ограниченной точности модели (Bao et al., 2018).

Полученные оценки влажности верхнего слоя почвы как качественные по оптическим данным (рис. 6), так и количественные по радарным данным (учитывая значения, не превышающие 35%) и совместном использовании радарных и оптических данных (рис. 7) показали их непротиворечивость, общие тренды и возможность сравнить степень влажности верхнего слоя почвы для разных участков дистанционно при отсутствии наземных измерений.