Using an unmanned aerial vehicle for soil moisture remote sensing by means of ultra-wideband electromagnetic impulses

Full Text

ВВЕДЕНИЕ

Внедрение роботизированных интеллектуальных технологий точного земледелия набрало инерцию поступательно развивающегося процесса. По прогнозам, в ближайшем будущем весь цикл полевых работ будет выполняться без участия человека с использованием беспилотных сельхозмашин и авиаплатформ [1]. До настоящего времени в производственный процесс ресурсосберегающего земледелия не внедрена технология дистанционного зондирования влажности почв (один из важнейших параметров, определяющий норму и сроки полива, оптимальную глубину и время сева, норму внесения удобрений и др.). Локальные измерения влажности почвы, осуществляемые в настоящее время путем отбора почвенных образцов (термостатно-весовой метод) или с использованием датчиков влажности почвы (методом рефлектометрических измерений во временной или частотной области, емкостным методом, нейтронным методом и т.д.), содержат ограниченную информацию о небольшой доли общей площади и затрудняют оценку пространственной неоднородности влажности почвы в пределах всего поля.

Применение беспилотных летательных аппаратов (БПЛА) для дистанционного зондирования влажности почв (а также состояния посевов) на больших площадях с высоким пространственным разрешением, по-видимому, является наиболее перспективной экономически обоснованной технологией, которая в ближайшем будущем войдет в роботизированные комплексы точного земледелия [1]. В последние десятилетия ведутся исследования по зондированию влажности почвы с использованием малых БПЛА методами микроволновой радиометрии в L-диапазоне частот (1.4 ГГц) [2, 3], пассивной радиолокации на основе приема отраженных от подстилающей поверхности сигналов глобальных навигационных спутниковых систем на частотах (1.2…1.6 ГГц) [4], сверхширокополосного (СШП) импульсного зондирования в надир (в диапазоне частот от сотен мегагерц до первых единиц гигагерц) [5]. Метод бокового обзора в БПЛА-рефлектометрах преимущественно используется для построения карт радиоизображения минных полей, а применительно к зондированию влажности почв широко не распространен. Регистрируемое в методе бокового обзора обратное (диффузное) рассеяние волны (на мелкомасштабных по сравнению с длиной волны неровностях поверхности почвы) в большей степени зависит от значения среднеквадратических отклонений высот неровностей поверхности почвы, чем от влажности ее поверхности [6], что ограничивает его применение. Являющийся развитием методик подповерхностного георадарного зондирования с приподнятой над поверхностью почвы антенной [7] метод зондирования влажности почвы в надир СШП-импульсами МГц-диапазона с борта БПЛА [5] в настоящее время представляется наиболее перспективным. (Компактных конструкций радиометров МГц-диапазона, используемых для измерения влажности почвы с борта малых БПЛА, в литературе нами не найдено.)

Надежность и воспроизводимость импульсных СШП-измерений влажности почвы обеспечивается проводимой внешней калибровкой радарной БПЛА-системы на различных высотах расположения фазового центра антенны над зондируемой поверхностью. В зависимости от метода калибровки радарной системы разделяют два способа восстановления влажности почвы:

а) как результат решения обратной задачи в частотной области в ходе минимизации нормы невязки между измеренной и рассчитанной функциями Грина задачи (диэлектрическое полупространство) [5];

б) на основе измерения модуля коэффициента отражения по максимуму огибающих зондирующих импульсов [8].

В первом численно-аналитическом методе [5] в частотной области выполняется полная калибровка передаточной характеристики антенно-фидерного тракта, представленного в виде четырехполюсника, элементы матрицы рассеяния которого находятся при различной высоте фазового центра антенны над эталонным отражателем (гладкая водная поверхность, металлический экран). В данном методе первичная измеряемая величина векторным анализатором цепей – это коэффициент отражения волны от зажимов антенны S₁₁(f), где f – частота. Для решения обратной задачи выбирается частотный диапазон, в котором наблюдаются наименьшие погрешности между моделью и измеренными амплитудными и фазовыми значениями функции Грина задачи. Во втором подходе [8] модуль коэффициента отражения измеряется как результат отношения максимумов огибающих импульсов, отраженных от подстилающей поверхности и калибровочного экрана (металл) либо от водной поверхности (с корректируемой поправкой на известный коэффициент отражения от воды).

Значительное влияние на надежность и воспроизводимость измерений влажности оказывает степень и характер мелкомасштабных неровностей поверхности почвенного покрова. Случайные высоты неровностей поверхности почвенного покрова существенно искажают спектр и временную форму зондирующего импульса, особенно это явление заметно при возрастании частоты в СВЧ-диапазоне (см. спектральное представление функций Грина в работе [9]). В МГц-диапазоне частот <1 ГГц амплитудный и фазовый спектр функции Грина имеют существенно меньшие флуктуации, обусловленные случайным рассеянием волн на неровностях поверхности почвы [10]. Отметим, что найденная таким образом максимальная частота, выше которой наблюдаются существенные случайные флуктуации амплитуды и фазы функции Грина задачи, может рассматриваться в качестве критерия гладкости зондируемой поверхности. При этом погрешность восстановления объемной влажности лежит в пределах от 0.04 см³/см³ [38] до 0.07 см³/см³ [9] без учета в алгоритмах инверсии высот неровностей поверхности почвы (от ~1 до ~5 см). Теоретически [11] и экспериментально (в лабораторных условиях на песке) [9] на лесных почвах с опадом [12], плоских участках строительных конструкций [13] было показано, что с разной степенью точности амплитуда и спектр (в диапазоне частот от сотен мегагерц до первых единиц гигагерц) отраженных СШП-импульсов могут быть описаны с использованием коэффициента отражения Френеля для гладкой поверхности с множителем в виде exp (–2k₀ ²σ_ш ²cos²θ^п) [13] (когерентная модель). Здесь θ^п — угол падения импульсной волны, k₀ = 2π/λ₀ = 2πf₀/c — волновое число в свободном пространстве, λ — длина волны, c — скорость света в вакууме, σ_ш — среднеквадратическое отклонение (СКО) высот неровностей зондируемой поверхности, f₀ – средняя частота зондирующего импульса. Стоит отметить, что, по-видимому, в результате аппроксимации, измеренной в широком диапазоне частот функции Грина (на одной какой-либо высоте) [5, 9–11] или усреднения максимумов огибающих зондирующих импульсов, измеренных на различных высотах [8], выполняются оценки когерентной составляющей функции Грина задачи и модуля коэффициента отражения соответственно. Однако до настоящего времени все еще недостаточно исследованы вопросы влияния мелкомасштабных неровностей поверхности сельскохозяйственных почв, находящихся под паром в естественных условиях периодической обработки, на флуктуации регистрируемых амплитуд зондирующих СШП-импульсов.

Следующим источником погрешности в методах дистанционного зондирования влажности почвы является неоднородность вертикального распределения влаги в поверхностном слое почвы. Сверхширокий спектр зондирующих импульсов принципиально может быть использован для прямого восстановления профилей влажности почвы. В связи с плохой обусловленностью данного класса обратных задач в большинстве случаев решение ищется на различных наборах параметрически заданных функций, описывающих глубинное распределение влаги в поверхностном слое почв [14]. Практическое решение подобных обратных задач получено для простых случаев, например, в работе [15] в качестве модельной среды использовался песчаник, комплексная диэлектрическая проницаемость (КДП) которого не обладает частотной дисперсией, а в работе [14] поверхность почвы была гладкой. Теоретические расчеты показывают [14], что в случае, если неоднородно увлажненная по глубине высыхающая почва интерпретируется в виде диэлектрически однородного полупространства, то восстановленные значения влажности совпадают со средней влажностью почвы, рассчитанной в слое различной толщины от ~2 см (L-диапазон) до ~7 см (P-диапазон); зондирование осуществлялось в надир, информативным признаком выступал модуль коэффициента отражения. В целом, как показано в экспериментах, использование L-диапазона частот в БПЛА-радиометрах ограничивает глубину зондирования поверхностного слоя почвы толщиной ~1.0…2.5 см [16]. Эксперименты показывают, что уменьшение частоты зондирования радиометра с 1.4 ГГц до 750 МГц, позволят увеличить глубину зондирования почвы (в зависимости от профиля влажности) до ~7…10.5 см [17]. Уменьшение рабочей частоты БПЛА-радиометров позволяет не только увеличить глубину зондирования, но и уменьшить рассеяние волны на неровностях поверхности почвы и элементах растительного покрова. В связи с этим может быть использован зарезервированный для радиоастрономических целей диапазон частот 406.1…410 МГц, также перспективным является развитие многочастотных [18] и импульсных автокорреляционных [19] СШП-методов радиометрического зондирования для платформ БПЛА. Эксперименты с СШП-импульсами МГц-диапазона частот от 200 МГц до 1.6 ГГц дают близкие с радиометрическим методом P-диапазона оценки глубины зондирования влажности почвы (содержание глинистой фракции 13…21 %) ~6 см [20]. Однако в настоящее время все еще остается до конца непонятной взаимосвязь между восстанавливаемыми (в приближении однородного полупространства) значениями влажности неоднородно увлажненной по глубине почвы из дистанционных измерений в МГц-диапазоне частот и средними значениями влажности, измеряемыми in-situ в поверхностном слое почвы различной толщины. В целом анализ показывает, что актуальность вопроса об оценке глубины зондирования влажности сельскохозяйственных почв дистанционными методами в МГц-диапазоне частот сохраняется.

Отметим также, что одним из факторов, влияющих на погрешность методик определения влажности почв [5, 8] является использование формулы Топпа [21], связывающей показатель преломления с объемной влажностью почвы. КДП почвы является входным параметром как в функции Грина задачи в методе, описанном в [5], так и в коэффициенте отражения в методах работ [8]. Подобная простая формула не учитывает частотную дисперсию КДП почв, а также была получена для спектра конкретно используемого импульса рефлектометра [21]. Частотная дисперсия КДП почв в МГц-диапазоне приводит к искажению частотного спектра и смещению средней частоты в низкочастотную область зондирующего СШП-импульса при его распространении в почве. Например, в работе [22] показано, что в модельных песчаных почвах при увеличении содержания бентонитовой глины с 0 % до 25 % средняя частота зондирующего импульса уменьшается примерно с 646…665 до 568…587 МГц.

Проведенный анализ показывает, что до настоящего времени возможности импульсного СШП-зондирования с БПЛА носителей влажности сельскохозяйственных почв в ходе их естественного севооборота, особенно на территории России, все еще слабо изучены. В данной статье исследованы возможности дистанционного зондирования СШП-импульсами с платформы БПЛА (квадрокоптер) влажности почв различной степени шероховатости, расположенных на полях ОПХ “Минино” ФИЦ КНЦ СО РАН, в районе пос. Минино Красноярского края.

1. ОРГАНИЗАЦИЯ ЭКСПЕРИМЕНТОВ

1.1. Местоположение и описание тестовых участков

На территории ОПХ “Минино” ФИЦ КНЦ СО РАН, расположенной в районе пос. Минино Емельяновского района, Красноярского края было выбрано три тестовых участка для проведения дистанционных рефлектометрических и контактных измерений влажности почвы, не покрытой растительным покровом (рис. 1). Почвы находились под паром в условиях естественного увлажнения и высыхания. В ходе экспериментов с 12 июня по 28 сентября 2022 г. на каждом из тестовых участков в различное время было подготовлено от одной до трех площадок с различной степенью шероховатости (от выровненной поверхности почвы граблями до естественной агропочвы после боронования). Размер каждой площадки составлял примерно 15 м × 15 м (рис. 2).

Рис. 1. Расположение тестовых участков в районе пос. Минино: координаты участка 1 — 56.0644 с.ш., 92.6967 в.д.; участка 2 — 56.0888 с.ш., 92.6660 в.д.; участка 3 — 56.0951 с.ш., 92.6654 в.д.

 

Рис. 2. Вид сверху на три подготовленные площадки с различной степенью шероховатости на тестовом участке 2.

 

Почвенный покров тестовых участков представлен выщелоченными и обыкновенными черноземными почвами; главная материнская порода — четвертичные лессовидные тяжелосуглинистые осадки, расположенные на древней денудационной равнине. Характеристики почвенных горизонтов тестовых участков в пахотном (А_пах) и подпахотном (АВ, В) слоях приведены в табл. 1.

 

Таблица 1. Характеристики почвенных горизонтов тестовых участков

Горизонт	D, см	H, %	ОВ, %	рНв.в, ед. pH	ρсух, г/см3	ИVКО, ммоль	MФ, %		MГФ**, %
							≤0.001*	≤0.01*
Участок 1
Апах	0—10	22.5	6.9	7.4	0.98	50	35.0	54.0	39.4
Апах	10—22	17.0	4.8	7.5	1.02	42	34.0	50.0	37.7
А	22—42	16.5	4.8	7.5	1.17	40	32.1	47.0	35.5
Участок 2
Апах	0—20	17.7	4.2	7.8	1.0	60	35	57	40.1
АB	20—30	17.1	4.0	8.0	1.1	46	39	61	44.1
B	30—45	15.7	3.2	7.8	1.1	44	38	55	41.9
Участок 3
Апах	0—20	22.3	7.2	7.9	1.1	60	28.1	53	33.8
АB	20—35	22.4	6.8	8.1	1.0	54	29.8	53	35.1
B	35—45	17.1	4.6	8.1	1.13	48	31.2	54	36.4

Примечания. D — глубина, H — содержание гумуса, ОВ — органическое вещество, рН_в.в — рН водной вытяжки, ρ_сух — плотность сухого сложения почвы, V_КО — емкость катионного обмена, M_Ф — весовое содержание фракции, M_ГФ — весовое содержание глинистой фракции. *Размер частиц приведен в мм. **Указана по стандарту USDA (U.S. Department of Agriculture), содержание физической глины по Качинскому соответствует фракции ≤ 0.01 мм.

 

Для каждой площадки зондирования были измерены следующие параметры шероховатости: среднеквадратическое отклонение и длина корреляции высот неровностей поверхности почвы. Параметры шероховатости рассчитывались по пяти различным профилям (длиной 3…5 м), отложенным на поверхности почвы с различным азимутом на основе цифровой модели рельефа (ЦМР), построенной для каждой из тестовых площадок синхронно с моментом зондирования БПЛА. ЦМР площадок была построена с использованием фотограмметрического метода (фотосъемка DJI AIR 2S с высоты 5 м под различными азимутальными углами). Локальные координаты задавались масштабной трехмерной линейкой (изготовленной из взаимно перпендикулярных деревянных брусков длиной 40 см и 1 м с сечением 19 мм × 40 мм), установленной на поверхности почвы.

В центре каждой площадки зондирования в пределах границ 3 м × 6 м в 30…80 точках с использованием датчика Decagon GS3, стержнями которого поверхность почвы протыкалась вертикально вниз до глубины 5 см, проводилось измерение влажности почвы. (Эффективная толщина слоя измерения влажности порядка 6…7 см.) Для датчиков GS3 использовалась заводская калибровка для минеральных почв. Также на каждой из площадок зондирования в одной точке измеряли влажность в поверхностном слое 6…7 см и по вертикальному профилю на глубине 3, 6, 9 и 15 см. Корреляционный анализ показал значимое расхождения между одиночными точечными и площадными измерениями влажности почвы в поверхностном слое 6…7 см, коэффициент детерминации R² = 0.29, СКО 3.5 % и систематическая ошибка +3.3 %. В результате влажность почвы, восстановленную по дистанционным импульсным СШП-измерениям, сопоставляли с площадными измерениями влажности почвы в поверхностном слое 6…7 см.

1.2. СШП БПЛА-рефлектометр

С использованием портативного векторного анализатора цепей (ВАЦ) CABAN R60 и ЛПА на диэлектрической подложке [23] был собран моностатический СШП-рефлектометр. ЛПА закреплена на БПЛА, так что максимум диаграммы направленности ориентирован в надир. SMA-разъем ЛПА соединен коаксиальным кабелем длиной 50 см с ВАЦ CABAN R60. СШП-рефлектометр был расположен на квадрокоптере (рама Tarot 650), снабженном лазерным дальномером LIDAR-Lite v3 (ошибка измерения расстояния ± 2.5 см на дистанциях до 40 м). Калибровочная плоскость лазерного дальномера практически совпадала с фазовым центром ЛПА. БПЛА управлялся полетным контроллером Pixhawk-4 и микрокомпьютером LattePanda V1. Разработанное программное обеспечение (на основе языка Python и библиотеки MAVSDK) управляло и синхронизировало движение БПЛА и рефлектометрические измерения ВАЦ CABAN R60. Управление ВАЦ CABAN R60 осуществлялось программно с использованием стандартных текстовых команд программируемых инструментов (SCPI) по протоколу TCP/IP (STREAM socket). Симулятор jMAVSim использовался для верификации алгоритмов работы всего комплекса в ходе выполнения полетных заданий. Внешний вид БПЛА рефлектометра приведен на рис. 3.

 

Рис. 3. Внешний вид БПЛА-рефлектометра.

 

Рефлектометрические измерения с борта БПЛА проводились в одной центральной точке на каждой из площадок зондирования при зависании БПЛА на девяти различных высотах над поверхностью почвы (примерно от 80…120 см до 6…7 м). Десятое измерение было калибровочным и выполнялось на высоте примерно 20 м.

2. МЕТОДИКА ИЗМЕРЕНИЙ КОЭФФИЦИЕНТА ОТРАЖЕНИЯ С ПЛАТФОРМЫ БПЛА

Конструкция ЛПА, состоящая из двухпроводной линии с симметричными вибраторами, не позволяет при ударном или сверхширокополосном непрерывном возбуждении сформировать излучаемый импульс, содержащий несколько колебаний поля [24]. В работах [23, 25] была предложена методика, позволяющая в присутствии конструкции БПЛА, сформировать ЛПА СШП-импульсы наносекундной длительности. В данной методике антенно-фидерный тракт БПЛА-рефлектометра, размещенного в воздушной среде над отражающим полупространством, представлен в виде четырехполюсника, элементы S-матрицы которого —

$r\left(f,d\right)=r_0\left(f\right)+G\left(f,d\right)Tr\left(f\right)$ (1)

– отыскиваются в ходе калибровки. Здесь r (f, d) — коэффициент отражения от входного соединения антенны (закрепленной на конструкции БПЛА), измеряемый ВАЦ CABAN R60, r₀ (f) – коэффициент отражения от входного соединения антенны при размещении БПЛА в неограниченном пространстве, Tr(f) – комплексная передаточная функция антенны, $G\left(f,d\right)\approx R\left(f\right)g\left(f,d\right)$ – функция Грина задачи, $g\left(f,d\right)=\text{exp}\left[4\mathrm{\pi }ifd\text{/}c\right]\text{/}\left(8\mathrm{\pi }d\right)$ – функция Грина в неограниченном пространстве, R(f) – коэффициент отражения волны, нормально падающей на границу раздела, d – расстояние от фазового центра антенны до границы воздух–полупространство. Эмпирические функции

$\mathbf{х}\left(f\right)=\left(\begin{array}{c}{r}_0\left(f\right)\\ Tr\left(f\right)\end{array}\right)$

могут быть найдены методом наименьших квадратов

$\underset{\mathrm{х}\left(f\right)}{\text{min}}{||\mathbf{Ах}\left(f\right)-\mathbf{В}||}_2$, x(f) = (A^TA)^–1A^TB

из переопределенной системы уравнений B = Ax, где

$\mathbf{В}=\left(\begin{array}{c}r\left(f,d_1\right)\\ ⋮\\ r\left(f,d_N\right)\end{array}\right),\mathbf{\text{А}}=\left(\begin{array}{cc}1& -g\left(f,d_1\right)\\ ⋮& ⋮\\ 1& -g\left(f,d_N\right)\end{array}\right)$,

которая сформирована на основе измерений r (f, d_p) на различных высотах d_p (p = 1,…,N) зависания БПЛА над плоской поверхностью с известными отражательными свойствами. В работе [23] приведены абсолютные значения и фазы r₀ (f), Tr(f), найденные описанным выше способом в диапазоне частот от f_min= 200 МГц до f_max= 1.3 ГГц с использованием N = 9 измерений r (f, d_p) на различных высотах (от 87 см до 5.7 м) зависания БПЛА над отражателем (6 м × 7 м), изготовленным из латунных сеток. При этом в уравнении (1) полагается $R\left(f\right)\equiv -1$.

C использованием найденных значений функций r₀ (f), Tr(f) и формулы (1) может быть рассчитан аналитический сигнал:

$\begin{array}{c}\dot{s}\left(t,d\right)=\frac{1}{s_K^{\text{огиб}}}{\int }_{f_{\text{min}}}^{f_{\text{max}}}df{K}_{\alpha }\left(f\right)\times \\ \times \frac{\left\{r\left(f,d\right)-r_0\left(f\right)\right\}}{Tr\left(f\right)}\text{exp}\left(-2\mathrm{\pi }ift\right)\end{array}$, (2)

временная форма $s\left(t,d\right)=2\text{Re\hspace{0.17em}}\dot{s}\left(t,d\right)$ и огибающая $s^{\text{огиб}}\left(t,d\right)=2\left|\dot{s}\left(t,d\right)\right|$ синтезированного СШП-импульса.

В формуле (2)

$s_K^{\text{огиб}}={\left|{\int }_{f_{\text{min}}}^{f_{\text{max}}}df{K}_{\alpha }\left(f\right)\text{exp}\left(-2\mathrm{\pi }ift\right)\right|}_{\text{max}}$

– нормировочная константа, необходимая для учета спектра оконной функции K_α(f). Синтезированный с использованием оконной функции Гаусса K_α(f) = exp(–0.5[(f — f₀)/α]²) импульс оказался длительностью ∆t_и = 1.8 нс (по уровню половинной амплитуды) с шириной спектра ∆f_и = 558 МГц (по уровню –10 дБ) от 456 до 1014 МГц и со средней частотой f₀ = 731 МГц. Ширина полосы пропускания оконной функции Гаусса K_α(f) задавалась порядка 1 ГГц (по уровню –6 дБ) с α = (4π)^–1. Форма синтезированного СШП-импульса и его спектр при отражении от латунной сетки показаны на рис. 4.

 

Рис. 4. Огибающая (1), временная форма (2) и нормированный спектр (3) синтезированного СШП-импульса. Расстояние между фазовым центром антенны и отражающим экраном (латунная сетка) примерно 87 см.

 

Предложенный способ калибровки позволяет проводить импульсные измерения модуля коэффициента отражения |R_поч(f₀)| электромагнитной волны от диэлектрического полупространства (на средней частоте f₀ зондирующего импульса) по формуле

$\left|R_{\text{поч}}\left(f_0\right)\right|=2d{s}_{\text{поч},\text{max}}\left(d\right)$, (3)

где $s_{\text{поч,\hspace{0.17em}max}}\left(d\right)={\left.s_{\text{поч}}^{\text{огиб}}\left(t,d\right)\right|}_{\text{max}}$ – максимум огибающей импульса, отраженного от диэлектрического полупространства (почва), при измерении на высоте d.

При различной высоте зависании БПЛА над латунным экраном модуль коэффициента отражения, измеренный согласно формуле (3):

$\left|R_{\text{мет}}\left(f_0\right)\right|=2d{s}_{\text{мет,max}} \left(d\right),$

близок к единице $\left|R_{\text{мет}}\left(f_0\right)\right|$ = 1.017±0.005 (относительно единицы, абсолютная погрешность в среднем составила около 2 %, а максимальная абсолютная погрешность — 11 % (рис. 5)).

Время задержки с коэффициентом детерминации 0.999 линейно пропорционально высоте d_лид, измеренной лазерным дальномером (рис. 5). Среднеквадратическое отклонение времени запаздывания отраженных импульсов от линейной зависимости составляет 0.15 нс (рис. 5, кривая 3), что эквивалентно СКО ≈ 2.3 см при определении импульсным методом высоты положения фазового центра антенны над металлическим экраном (подобная ошибка сопоставима с точностью измерения лазерным дальномером LIDAR-Lite v3, а также максимальным отклонением от плоскости латунных сеток, лежащих на имеющем неровности асфальте).

 

Рис. 5. Время распространения t_расп (1) и максимум огибающей s_{мет,мак}(d_лид) (2), синтезированного СШП-импульса, отраженного от металлического экрана (латунной сетки), при различной высоте d_лид положения фазового центра антенны над экраном. Аппроксимирующие зависимости и соответствующие R² и СКО: кривая 3 — t_расп = (0.18 ± 0.03)+(6.54 ± 0.01) d_лид, R² = 0.999 и СКО = 0.15 нс; кривая 4 — s_{мет, max}(d_лид) = (1.017 ± 4.910^–3)/(2d_лид), R² = 0.986 и СКО = 0.014 м^–1.

 

3. РЕЗУЛЬТАТЫ И ОБСУЖДЕНИЕ

А. С достаточной степенью уверенности можно утверждать, что для всех тестовых площадок, на которых в полевых условиях проводились рефлектометрические измерения с помощью БПЛА, характерен нормальный закон распределения высот мелкомасштабных неровностей поверхности почвы. На рис. 6а представлена нормированная частота p^изм(h_ш) распределения высот h_ш($\overrightarrow{\mathit{r}}$) неровностей поверхности почвы (измерения выполнены 21.09.2022 г. на тестовом участке 2, площадка 1), где $\overrightarrow{r}$ – радиус-вектор, пробегающий по поверхности почвы. Нормированная частота подобных распределений с высоким значением коэффициента детерминации R² = 0.981 описывается гауссовой функцией (см., например, на рис. 6а кривую p^теор(h_ш), аппроксимирующую измеренные значения p^изм(h_ш)). Для всех тестовых площадок в течение всех экспериментов коэффициент детерминации между измеренными p^изм(h_ш) и аппроксимированными гауссовой функцией p^теор(h_ш) распределений относительных частот появления высот неровностей поверхности почвы оказался равен R² = 0.964 (рис. 6б). Отметим, что для всех профилей была убрана “постоянная составляющая” (крупномасштабные неровности — рельеф) с использованием фильтрации на основе преобразования Фурье.

 

Рис. 6. Характеристики мелкомасштабных неровностей поверхности почвы тестовых участков: а) распределение высот неровностей h_ш, поверхности почвы p^изм(h_ш), измеренных фотограмметрическим методом 21.09.2022 на тестовом участке 2, площадка 1, при R² = 0.981; p^теор(h_ш) – соответствующая аппроксимация гауссовой функцией; б) корреляция между измеренными и аппроксимированными (гауссова функция) распределениями h_ш для всех тестовых участков в ходе всех экспериментов (1), линейная регрессия (2), при R² = 0.964; в, г) распределение измеренных фотограмметрическим методом СКО σ_ш и длины корреляции l_к высот мелкомасштабных неровностей поверхности почв в целом для всех тестовых участков в ходе всех экспериментов.

 

В ходе экспериментов СКО высот неровностей поверхности почвы σ_r тестовых площадок варьировалось от 0.39 см до 1.73 см (см. гистограмму на рис. 6в). Для большинства (порядка 70 %) профилей высот неровностей поверхности почвы на тестовых площадках характерна экспоненциальная автокорреляционная функция (для остальных видов зависимостей данный вопрос детально не изучался). Длина корреляции высот неровностей поверхности почвы l_к тестовых площадок варьировалась от 4.1 до 11.7 см (рис. 6г). Погрешность измерения σ_ш оценивалась как СКО σ_ш0 относительно высот неровностей вдоль двух горизонтальных профилей, проведенных по абсолютно плоским поверхностям (пластей) двух взаимно перпендикулярных брусков масштабной линейки. Эти СКО составляли σ_ш0 ≤ 3.5 мм для 79 % и σ_ш0 ≤ 5.5 мм для 89 % профилей.

Измеренные значения влажности почвы в пределах каждой тестовой площадки имели также близкий к нормальному закон распределения. На рис. 7а приведена гистограмма нормированных частот распределения p^изм(W) измеренных значений влажности почвы 21.09.2022 г. на тестовом участке 2 (площадка 1) и соответствующая аппроксимирующая гауссова функция p^теор(W). Для всех тестовых площадок в ходе всех экспериментов между измеренными p^изм(W) и аппроксимированными p^теор(W) гауссовой функцией распределениями коэффициент детерминации оказался равен R² = 0.895 (рис. 7б). Проведенный дополнительный анализ не обнаружил какой-либо связи между средними значениями и СКО влажности почвы в пределах каждого из участков. При этом СКО влажности почвы в пределах каждого из тестовых участков оказываются менее 3.5 % и 5.5 % соответственно для 60 % и 96 % случаев относительно среднего значения.

 

Рис. 7. Распределение значений объемной влажности почвы: а) измеренных 21.09.2022 p^изм(W) на тестовом участке 2, площадка 1 при R² = 0.869, p^теор(W) – соответствующая аппроксимация гауссовой функцией; б) корреляция между измеренными и аппроксимированными (гауссова функция) распределениями W для всех тестовых участков в ходе всех экспериментов (1), линейная регрессия (2) при R² = 0.895.

 

Б. В ходе экспериментов на каждой из тестовых площадок были проведены измерения временных форм и огибающих СШП-импульсов, отраженных от почвенного покрова при различной высоте зависания БПЛА над центром площадок зондирования. На рис. 8 представлены временные формы s(t, d) и огибающие s^огиб(t, d) СШП-импульсов, измеренные 11 сентября 2022 г. на тестовом участке 1 (площадка 1) при высоте полета БПЛА от 1.01 до 5.11 м.

 

Рис. 8. Временные формы s(t, d_лид) (сплошные линии) и огибающие s^огиб(t, d_лид) (штриховые линии) СШП-импульсов, измеренные 11 сентября 2022 г. на тестовом участке 1, площадка 1 при σ_ш = 0.5 см, l_к = 7.5 см, W = 25.5 % и различной высоте зависания БПЛА: d_лид= 1.01 (1), 1.63 (2), 2.33 (3), 3.20 (4), 4.17 (5) и 5.11 м (6).

 

Время задержки и затухание максимумов огибающих импульсов, изображенных на рис. 8, приведены на рис. 9а для тестового участка 1 (площадка 2) 17 сентября 2022 г., характеристика почвенного покрова – σ_ш = 1.7 см, l_к = 7.9 см, W = 9.3 %. Между высотой полета БПЛА d_лид измеренной лазерным дальномером, и временем (по максимуму огибающей) прихода отраженных от почвенного покрова СШП-импульсов имеется линейная зависимость (рис. 9а). Между калибровочной линейной зависимостью (см. рис. 5 и рис. 9а, линия 3: d_лид = –0.03+0.153 t_расп), полученной при регистрации отраженных импульсов от латунного экрана, и зависимостями (см. рис. 9а, 1 и 2), полученными при регистрации импульсов, отраженных от почвенного покрова различной степени шероховатости и влажности, наблюдаются высокие значения коэффициента детерминации, 0.995…0.999. (Построенные на основе всех экспериментальных данных подобные зависимости также характеризуются столь высокими значениями коэффициентов детерминации.) Вместе с тем, большие СКО высот неровностей поверхности почвы σ_ш = 1.7 см (см. рис. 9а, символы 2) в сравнении с σ_ш = 0.5 см (см. рис. 9а, символы 1) приводят к большим погрешностям в определении высоты БПЛА импульсным методом, СКО = 10.1 и 2.3 см соответственно, относительно значений лазерного дальномера. В целом имеется слабая линейная зависимость между дисперсией времен задержек, отраженных от почвы импульсов, и СКО высот неровностей поверхности почвы R² = 0.28 (анализ проведен для всех тестовых участков в ходе всех экспериментов, график не приводится).

 

Рис. 9. Время запаздывания (а) и максимумы огибающих (б) СШП-импульсов, отраженных от почвенного покрова: а – измерения выполнены 11 и 17 сентября 2022 г. на тестовом участке 1 соответственно на площадках 1 (1) и 2 (2); б — аппроксимации измеренных значений s_{поч, max}(d_лид) (символы 1 и 2) линейными зависимостями обозначены соответственно 4 и 5. Линейные зависимости 3 получены в ходе калибровки (см. разд. 3Б).

 

Между максимумами огибающих СШП-импульсов, отраженных от почвенного покрова, и двойным обратным расстоянием 1/2d_лид имеется линейная зависимость. В качестве примера на рис. 9б (1 и 2) показаны измеренные значения максимумов огибающих импульсов, отраженных от почвенного покрова тестового участка 1, площадки 1 и 2. Данные зависимости с коэффициентами детерминации R² = 0.989 и 0.928 и СКО, равными 6.6⋅10^–3 и 11.8⋅10^–3 1/м, описываются линейными зависимостями (см. рис. 9б, символы 1 и линия 4, символы 2 и линия 5). Линейная зависимость (см. рис. 9б, символ 3) s_{мет, max}(d_лид) = 1.017/2d_лид получена в ходе калибровки БПЛА-рефлектометра. Аналогичные линейные зависимости наблюдаются для всех подобных экспериментальных данных.

В соответствии с моделью (3) наклон линейных зависимостей (рис. 9б) определяет модуль коэффициента отражения. В качестве сопоставления на рис. 9б приведена калибровочная линия 3 (рис. 5), которая описывает значения максимумов огибающих импульсов s_{мет, max}(d_лид), отраженных от латунной сетки (тангенс угла наклона 1.017). Отметим, что значение модуля коэффициента отражения, определенное по наклону линий (рис. 9б) является оценкой когерентной компоненты, так как она получена в результате усреднения амплитуд отраженных от почвенного покрова импульсов при различной высоте полета БПЛА (кроме того, каждая точка, изображенная на рис. 9б, является результатом пяти измерений, которые проводились при зависании БПЛА на фиксированной высоте, сопровождавшихся небольшими вертикальными смещениями БПЛА в силу случайных факторов). В случае нормального закона распределения высот неровностей поверхности почвы (без учета вида автокорреляционной функции и параметра — длины корреляции) среднее значение (когерентная компонента) модуля коэффициента отражения может быть записана в виде [13]

$\begin{array}{c}\left|R^{\text{теор}}\left(f_0\right)\right|=\left|R_{\text{глад}}\left(f_0\right)\right|S_{\text{ш}}\left({\mathrm{\sigma }}_{\text{ш}}\right), \\ S_{\text{ш}}\left({\sigma }_{\text{ш}}\right)=\text{exp}\left[-2{\left(\frac{2\mathrm{\pi }{f}_0}{c}{\mathrm{\sigma }}_{\text{ш}}\right)}^2\right]\end{array}$, (4)

|R_глад(f₀)| – модуль коэффициента отражения Френеля от диэлектрического полупространства с гладкой границей для нормально падающей волны частотой f₀, S_ш(_ш) – фактор шероховатости.

В. В соответствии с методикой расчета модуля коэффициента отражения по наклону линейных зависимостей (см. рис. 9б и формулу (3)) были найдены модули коэффициентов отражения для всех тестовых площадок в ходе всех экспериментов. На рис. 10 приведены зависимости экспериментально измеренных модулей коэффициентов отражения $\left|R^{\text{изм}}\left(f_0\right)\right|$ (нормированных на фактор шероховатости S_ш(_ш)) от объемной влажности почвы, измеренной на тестовых участках. На рис. 10 отрезками отмечены интервалы СКО в определении коэффициентов отражения по линейным зависимостям (рис. 9б) и СКО объемной влажности почвы, измеренной на каждой из тестовых площадок. Наблюдаемая влажностная зависимость модуля коэффициента отражения (рис. 10), хорошо описывается следующей моделью:

$R^{\text{теор}}\left(f_0,W\right)|=\left|\frac{1-\sqrt{{\epsilon }_{\text{поч}}\left(f_0,{\mathrm{\rho }}_{\text{сух}},M_{\text{ГФ}},W\right)}}{1+\sqrt{{\epsilon }_{\text{поч}}\left(f_0,{\mathrm{\rho }}_{\text{сух}},M_{\text{ГФ}},W\right)}}\right|$, (5)

где ε_поч(f₀, ρ_сух, M_ГФ, W) – КДП почвенного покрова, ρ_сух – плотность сухого сложения, M_ГФ – фракция частиц в почве менее 2 мкм (весовое содержание физической глины по USDA, см. табл. 1). КДП почвенного покрова рассчитывали на основе диэлектрической модели Миронова [26]. Плотность сухого сложения задавалась равной 0.98 г/см³ среднему значению, измеренному в слое 0…6 см на всех тестовых полях. Фракция физической глины задавалась равной 37.8 % среднему значению, измерена (агрохимцентр ФГБУ ГЦАС “Красноярский”) в пахотном слое 0…20 см тестовых участков (см. табл. 1). Коэффициент детерминации и СКО между измеренными и рассчитанными значениями |R(f₀)| оказались равны R² = 0.693, СКО = 4.8⋅10^–3 (рис. 10).

 

Рис. 10. Зависимость нормированного на фактор шероховатости модуля коэффициента отражения, измеренного БПЛА-рефлектометром (1), в зависимости от объемной влажности почвы in-situ; модуль коэффициента отражения 2 рассчитан по формуле (5).

 

Разделив измеренные значения коэффициентов отражения (не нормированных на фактор шероховатости) на значения коэффициентов отражения, рассчитанные по модели (5), построена зависимость измеренного фактора шероховатости от СКО высот неровностей поверхности почвы (рис. 11).

 

Рис. 11. Зависимость фактора шероховатости S_ш(σ_ш), измеренного (1) на тестовых участках, от СКО высот неровностей поверхности почвы, σ_ш; функция распределения Гаусса (2) аппроксимирующая (1); фактор шероховатости (3), рассчитанный по модели (4).

 

Модель (4) фактора шероховатости S_ш(σ_ш) удовлетворительно описывает соответствующие измеренные значения (рис. 11, кривая 3). В силу того, что в формуле (4) в фактор шероховатости не входят свободные параметры для оптимизации, на рис. 11 для сравнения приведена также функция распределения Гаусса (кривая 2), аппроксимирующая экспериментальные зависимости (1) со СКО = 0.11. Видно, что зависимость, аппроксимирующая экспериментальные данные (рис 11, кривая 2), практически совпадает с моделью (4) фактора шероховатости (кривая 3), что дополнительно указывает на справедливость применения закономерности (4) для фактора шероховатости S_ш(σ_ш).

Г. Решим обратную задачу по нахождению объемной влажности почвы W в ходе минимизации функционала, представляющего собой норму невязки между измеренными $\left|R^{\text{изм}}\left(f_0\right)\right|$ и рассчитанными |$R^{\text{теор}}\left(f_0,{}_{\text{сух}},M_{\text{ГФ}},W\right)|$ значениями модулей коэффициентов отражения:

$W_j={\left.{\left[\left.\left|R_j^{\text{изм}}\left(f_0\right)\right.\right|-\left|\left.R_{}^{\text{теор}}\left(f_0,{\mathrm{\rho }}_{\text{сух}},M_{\text{ГФ}},W_j\right)\right|\right.\right]}^2\right|}_{\text{min}}$,(6)

где индекс j обозначает измерение коэффициента отражения на отдельной тестовой площадке.

Найдем объемные влажности двумя способами. В первом случае будем использовать измеренные значения $R_j^{\text{изм}}\left(f_0\right)$ без нормировки на фактор шероховатости, во втором случае с учетом нормировки (рис. 11, кривая 3). Найденные таким образом объемные влажности почвы W_вост относительно синхронно измеренных in-situ значений W_изм на тестовых площадках показаны на рис. 12. Найденные двумя способами объемные значения влажности почвы со СКО = 3.5…3.8 % совпадают с измеренными in-situ. При этом без учета фактора шероховатости найденные значения объемных влажностей почвы оказываются заниженными относительно измеренных не более чем на –2.1 % (в среднем на –1.2 %). Значения объемных влажностей почвы, найденных с учетом нормировки модуля коэффициента отражения, на фактор шероховатости практически не смещены (+0.2 %) относительно измеренных in-situ. Коэффициент детерминации между найденными и измеренными значениями объемной влажности почвы оказался равен 0.725 (с учетом S_ш(σ_ш)) и 0.698 (без учета S_ш(σ_ш)).

 

Рис. 12. Зависимость найденных в ходе решения обратной задачи и измеренных in-situ значений объемной влажности почвы: решение обратной задачи с корректировкой (1) и без корректировки (2) коэффициента отражения за счет фактора шероховатости поверхности почвы и соответствующие линии регрессии (3, 4).

 

ЗАКЛЮЧЕНИЕ

В результате долговременных систематически проведенных экспериментов по дистанционному зондированию почвенного покрова сельскохозяйственных полей, находящихся под паром в ходе естественных условий увлажнения и испарения, подтверждена практически значимая точность (3.5…3.8 %) измерения СШП-импульсами с борта БПЛА влажности почв. Поверхность почвенного покрова характеризовалась различными степенью и характером шероховатости в ходе периодического боронования и естественного разглаживания неровностей после дождей. Указанная погрешность измерения влажности почв находится в пределах естественных вариаций влажности почвы площадок зондирования.

Хорошее соответствие между значениями влажности почв, найденных из дистанционных наблюдений и измеренных контактным методом in-situ, позволяют судить о порядке толщины поверхностного слоя зондирования не менее чем 6…7 см (на частоте 731 МГц) в среднем для всех видов профилей влажности, степени и вида шероховатости поверхности почвенных покровов, складывавшихся на тестовых участках в ходе экспериментов. Данная оценка глубины зондирования близка к оценкам, полученным в работах [17] для сопоставимой частоты, и примерно в два раза больше (что соответствует кратности длин волн) по отношению к частоте зондирования 1.4 ГГц [16].

В диапазоне частот зондирующего импульса от 456 до 1014 МГц можно пренебречь влиянием диффузного рассеяния волн на случайных неровностях поверхности почвенного покрова тестовых участков (СКО высот неровностей меньше 2 см) в пределах указанной выше погрешности восстановления влажности почвы.

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

Авторы заявляют об отсутствии конфликта интересов.

ФИНАНСИРОВАНИЕ РАБОТЫ

Работа выполнена при финансовой поддержке Российского научного фонда (проект № 22-17-20042) и Красноярского краевого фонда науки (дог. № 42).

About the authors

K. V. Muzalevskiy

Kirensky Institute of Physics Federal Research Center KSC Siberian Branch Russian Academy of Sciences

Author for correspondence.
Email: rsdkm@ksc.krasn.ru
Russian Federation, Akademgorodok 50, bld. 38, Krasnoyarsk, 660036

S. V. Fomin

Kirensky Institute of Physics Federal Research Center KSC Siberian Branch Russian Academy of Sciences

Email: rsdkm@ksc.krasn.ru
Russian Federation, Akademgorodok 50, bld. 38, Krasnoyarsk, 660036

A. Yu. Karavayskiy

Kirensky Institute of Physics Federal Research Center KSC Siberian Branch Russian Academy of Sciences

Email: rsdkm@ksc.krasn.ru
Russian Federation, Akademgorodok 50, bld. 38, Krasnoyarsk, 660036

Z. Ruzicka

Kirensky Institute of Physics Federal Research Center KSC Siberian Branch Russian Academy of Sciences

Email: rsdkm@ksc.krasn.ru
Russian Federation, Akademgorodok 50, bld. 38, Krasnoyarsk, 660036

Yu. A. Leskova

Kirensky Institute of Physics Federal Research Center KSC Siberian Branch Russian Academy of Sciences

Email: rsdkm@ksc.krasn.ru
Russian Federation, Akademgorodok 50, bld. 38, Krasnoyarsk, 660036

A. G. Lipshin

Krasnoyarsk Research Institute of Agriculture Federal Research Center KSC Siberian Branch Russian Academy of Sciences

Email: rsdkm@ksc.krasn.ru
Russian Federation, Svobodny Prosp., 66, Krasnoyarsk

V. N. Romanov

Krasnoyarsk Research Institute of Agriculture Federal Research Center KSC Siberian Branch Russian Academy of Sciences

Email: rsdkm@ksc.krasn.ru
Russian Federation, Svobodny Prosp., 66, Krasnoyarsk

References

Supplementary files