СПОСОБ ДИСТАНЦИОННОГО ОПРЕДЕЛЕНИЯ ДЕГРАДАЦИИ ПОЧВЕННОГО ПОКРОВА
METHOD FOR REMOTE DETERMINATION OF DEGRADATION OF SOIL COVER

Изобретение относится к сельскому хозяйству и может найти применение при составлении государственного кадастра пахотных земель.
Почва - многокомпонентная среда. По определению основоположника почвоведения В.В. Докучаева [см, Н.Ф. Ганжара, «Почвоведение», учебник, межрегиональная ассоциация "Агрообразование", 2001 г., стр.5-9] почва - самостоятельное, естественноисторическое биокосное природное тело, представляющее собой открытую четырехфазную динамичную систему, включающую твердую фазу (полимерная органоминеральная система), жидкую фазу (вода, занимающая часть перового пространства с растворенными органическими и минеральными веществами), газовую фазу (почвенный воздух), живую фазу (почвенная биота, населяющие почву микроорганизмы).
По мере развития промышленности, прилегающие к промышленным зонам территории, подвергаются все усиливающемуся потоку техногенных выбросов и антропогенных нагрузок, приводящих к деградации почв.
Однако, на настоящий момент нет согласованного перечня как показателей деградации почв, так и методов их определения. Не установлены так же ПДК фоновых уровней загрязнителей.
Известен способ оценки устойчивости почв в условных баллах (см, «Почвенно-экологический мониторинг» учебное пособие под редакцией Д.С. Орлова, из-во МГУЛ, 1994, стр.77-76, таблица 18, «Экспертная оценка почвенного покрова по показателям, определяющим его интегральную устойчивость (баллы)» - аналог).
В способе-аналоге для оценки интегральной устойчивости определенной группы почв использованы 10 показателей: кислотность, увлажненность, интенсивность разложения растительных остатков, запас гумуса, рельеф, сельскохозяйственная освоенность, годичный прирост, насыщенность основания, почвообразующие породы, теплообеспеченность.
По сумме баллов выделены зоны: до 14 баллов - неустойчивая, от 18 до 20 - малоустойчивая, от 23 до 26 относительно устойчивая и устойчивая от 26 до 30 баллов.
Недостатками аналога следует считать:
- субъективность балльной экспертной оценки каждого из 10 показателей;
- неоперативность и трудоемкость наземного обследования больших площадей;
- недопустимо большое количество частных показателей, затрудняющих практическую применимость;
- отсутствие документально установленных границ зон.
Известна шкала ступеней (от 0 до 6) деградации почвы: [см, журнал, «Географический вестник», №2(17), 2011 г, статья «Экология и природопользование», стр.49-50] - аналог.
Таблица 1
Критерии для определения степени деградации почвы
Критерии
Степень деградации
0
1
2
3
4
5
Площадь обнаженного гумусового горизонта (А), %
0
<10
10-20
21-50
51-90
>90
Мощность абиотического наноса, см
0
<2
2-10
11-20
21-40
>40
Площадь обнаженной почвообразующей породы (С) или подстилающей породы (D), % от общей площади
0-2
3-5
6-10
11-15
16-25
>25
Уменьшение мощности почвенного профиля (А+В), % от исходного
0-1
1-3
3-25
26-50
51-75
>75
Уменьшение запасов гумуса в профиле почвы (А+В), % от исходного
<5
5-10
11-20
21-40
41-80
>80
Таблица 2
Дополнительные критерии определения степени деградации почв
Критерии
Степень деградации
0
1
2
3
4
5
Уменьшение содержания микроэлементов (Mn, Co, Mo, В, Cu, Fe) в % от средней степени обеспеченности
<5
5-10
11-20
21-40
41-80
>80
Уменьшение содержания подвижного фосфора, в % от средней степени обеспеченности
<5
5-10
11-20
21-40
41-80
>80
Уменьшение степени кислотности (pH сол.), в % от средней степени кислотности
<5
5-10
11-15
16-20
21-25
>25
Потери почвенной массы т/га/год
<2
2-5
6-25
26-100
100-200
>200
К недостаткам способа-аналога следует отнести:
- неопределенность расчета интегрального показателя при разных значениях частных составляющих;
- субъективность определения значений частных составляющих, неоперативность и большая трудоемкость их количественного расчета;
- отсутствие документальных измерений границ выделенных зон.
Ближайшим аналогом к заявленному техническому решению является способ дистанционных измерений коэффициента спектральной яркости почв и обнажений посредством спектрометров (типа СПИ-2, СПИ-74), установленных на самолетном носителе [см Л.И. Чапурский «Отражательные свойства природных объектов в диапазоне 400-2500 нм», часть 1, Министерство Обороны СССР, 1986 г, §5. Коэффициенты спектральной яркости (КСЯ) почв и обнажений, стр.20-39]. В способе ближайшего аналога измеряют КСЯ различных типов почв, вдоль трассы полета носителя, представляющих собой семейство практически не пересекающихся функций ρ(λ), ординаты которых увеличиваются с ростом длины волны (λ).
λ
гум
750
гум
<mrow><msub><mi>ρ</mi><mrow><mn>750</mn></mrow></msub><mo>=</mo><mi>A</mi><mo>×</mo><mi>B</mi><mo>×</mo><mrow><mo>(</mo><mrow><mn>1</mn><mo>+</mo><mfrac><mn>1</mn><mrow><mi>K</mi><mo>×</mo><msub><mi>H</mi><mrow><mi>г</mi><mi>у</mi><mi>м</mi></mrow></msub></mrow></mfrac></mrow><mo>)</mo></mrow></mrow>
Вариации КСЯ как по амплитуде, так и по спектру зависят от концентрации гумуса, минерализации почв и влажности почвенных покровов. В частности, получены эмпирические зависимости КСЯ (ρ=750) на длине волны λ=750 нм от концентрации гумуса в почве Н(%) вида
, (А=50; K=0,026; В=8,4) при коэффициенте корреляции между ρ и Н (%) порядка 0,9.
Недостатками способа ближайшего аналога являются:
- отсутствие установленной интегральной зависимости между деградацией почвенного покрова и измеренными значениями КСЯ;
2
- невысокое пространственное разрешение спектрометров типа СПИ, не позволяющее выделить на изображениях участки деградации почвенных покровов площадью менее ~2 км.
Задача, решаемая изобретением состоит в оперативном, количественном измерении степени деградации почвенных покровов, по дистанционным мультиспектральным изображениям, с выделением и оконтуриванием участков деградации, размерами 0,1 га, пригодных для составления кадастра пахотных земель.
Технический результат достигается тем, что способ дистанционного определения деградации почвенного покрова включает регистрацию полей яркости I(x,y) пахотных земель, содержащих тестовые эталонные участки многоканальным спектрометром, установленным на аэрокосмическом носителе в зональных диапазонах 450-515, 525-605, 630-690, 750-900 и 15 50-17 50 нм с одновременным получением цифровых изображений в каждом канале, расчет частных индексов состояния по комбинации зональных отношений сигналов I(x,y) в каналах для каждого пикселя изображения; индекса содержания гумуса [H, %], индекса засоленности NSI и влажностного индекса NDWI, выбор в качестве интегрального критерия деградации почв многопараметрической функции регрессии произведения перечисленных индексов, пронормированных относительно их значений для эталонных тестовых участков в виде степенных зависимостей:
<mrow><mi>D</mi><mo>=</mo><msup><mrow><mrow><mo>(</mo><mrow><mfrac><mrow><msub><mi>H</mi><mn>0</mn></msub></mrow><mi>H</mi></mfrac></mrow><mo>)</mo></mrow></mrow><mrow><mn>1,9</mn></mrow></msup><mo>⋅</mo><msup><mrow><mrow><mo>(</mo><mrow><mfrac><mrow><mi>N</mi><mi>S</mi><mi>I</mi></mrow><mrow><mi>N</mi><mi>S</mi><msub><mi>I</mi><mn>0</mn></msub></mrow></mfrac></mrow><mo>)</mo></mrow></mrow><mrow><mn>0,5</mn></mrow></msup><mo>⋅</mo><msup><mrow><mrow><mo>(</mo><mrow><mfrac><mrow><mi>N</mi><mi>D</mi><mi>W</mi><msub><mi>I</mi><mn>0</mn></msub></mrow><mrow><mi>N</mi><mi>D</mi><mi>W</mi><mi>I</mi></mrow></mfrac></mrow><mo>)</mo></mrow></mrow><mrow><mn>0,3</mn></mrow></msup></mrow>
пересчет математическими процедурами значений каждого пикселя изображения в масштаб интегрального показателя D, выделение контуров деградации на результирующей матрице и их визуализацию в виде обработанных снимков, где:
D - степень (категория) деградации 1, 2, 3, 4;
1,9; 0,5; 0,3 - показатели чувствительности D к частным индексам.
Изобретение поясняется чертежами, где:
фиг.1 - спектральные кривые гумусовых горизонтов;
фиг.2 - спектральные кривые засоленных горизонтов (солончаки);
фиг.3 - спектральная характеристика почвенного покрова в полосе 750 нм от концентрации гумуса;
фиг.4 - вид обработанного снимка с оконтуренными участками деградации;
фиг.5 - функциональная схема устройства, реализующая способ.
Техническая сущность изобретения состоит в следующем. Особенностью дистанционного зондирования Земли из космоса является зависимость сигнала измерений от характеристик среды (атмосферы, условий съемки: угла визирования на объект, высоты Солнца, времени суток и т.д.), что вызывает необходимость калибровки тракта измерений. Последнее в заявленном способе достигается путем одновременной съемки тестовых (эталонных) участков, по значениям сигналов от которых осуществляют в дальнейшем калибровку трактов измерений объектов подстилающей поверхности.
Из математики известно, что для однозначного решения системы уравнений, количество уравнений должно быть равным числу неизвестных.
Неизвестными в заявленном способе являются индексы - параметры: концентрация гумуса, засоленность, влагопотери. Следовательно, для однозначного их расчета необходимо иметь не менее трех независимых сигналов, что реализуется синхронными измерениями в перечисленных зонах спектрометра.
При количественном расчете самих индексов - параметров используют различные комбинации зональных отношений [см, например, Виноградов Б.В. «Аэрокосмический мониторинг экосистем», Наука, М, 1984 г]. Основное требование к зональным отношениям - адекватность измерений физическому процессу, высокий коэффициент корреляции между сигналом и расчетным параметром, устойчивость к внешним условиям съемки.
Концентрацию гумуса в почве (H) рассчитывают из соотношения [Виноградов Б.В.]:
<mrow><mi>H</mi><mo stretchy="false">[</mo><mi>%</mi><mo stretchy="false">]</mo><mo>=</mo><mfrac><mn>1</mn><mi>k</mi></mfrac><mo>⋅</mo><mi>ln</mi><mfrac><mrow><msub><mi>ρ</mi><mn>0</mn></msub><mo>−</mo><msub><mi>ρ</mi><mrow><mi>min</mi></mrow></msub></mrow><mrow><msub><mi>ρ</mi><mi>λ</mi></msub><mo>−</mo><msub><mi>ρ</mi><mrow><mi>min</mi></mrow></msub></mrow></mfrac></mrow>
,
где:
λ
ρ - коэффициент спектрального отражения почвы λ=750-900 нм,
0
ρ - тоже для безгумусовой почвообразующей породы,
min
ρ - тоже для многогумусовой почвы,
k - для почв средней полосы ~0,6.
<mrow><mi>N</mi><mi>S</mi><mi>I</mi><mo>=</mo><mfrac><mrow><mi>S</mi><msub><mi>I</mi><mrow><mi>max</mi></mrow></msub><mo>−</mo><mi>S</mi><mi>I</mi></mrow><mrow><mi>S</mi><msub><mi>I</mi><mrow><mi>max</mi></mrow></msub><mo>−</mo><mi>S</mi><msub><mi>I</mi><mrow><mi>min</mi></mrow></msub></mrow></mfrac></mrow>
<mrow><mi>S</mi><mi>I</mi><mo>=</mo><msqrt><mrow><msubsup><mi>ρ</mi><mn>1</mn><mn>2</mn></msubsup><mo>+</mo><msubsup><mi>ρ</mi><mn>3</mn><mn>2</mn></msubsup></mrow></msqrt></mrow>
Индекс засоленности
, где
,
1
3
ρ, ρ - КСЯ пикселя в первом и третьем каналах.
NIR
SWIR
Для определения нормализированного влажностного индекса (NDWI) используется измерения интенсивности в ближнем инфракрасном диапазоне ρ с характерной длинной волны λ=750-900 нм а также в коротковолновом инфракрасном диапазоне ρ с длинами волн λ=1550-1750 нм. При вычислении индекса NDWI используется следующая формула:
<mrow><mi>N</mi><mi>D</mi><mi>W</mi><mi>I</mi><mo>=</mo><mfrac><mrow><msub><mi>ρ</mi><mrow><mi>N</mi><mi>I</mi><mi>R</mi></mrow></msub><mo>−</mo><msub><mi>ρ</mi><mrow><mi>S</mi><mi>W</mi><mi>I</mi><mi>R</mi></mrow></msub></mrow><mrow><msub><mi>ρ</mi><mrow><mi>N</mi><mi>I</mi><mi>R</mi></mrow></msub><mo>+</mo><msub><mi>ρ</mi><mrow><mi>S</mi><mi>W</mi><mi>I</mi><mi>R</mi></mrow></msub></mrow></mfrac></mrow>
Для выявления скрытых закономерностей между измеряемыми индексами - параметрами деградации, и их вклада в результирующий (интегральный) показатель деградации использована многопараметрическая функция регрессии.
В соответствии с рекомендациями математического института им. Стеклова, в качестве монотонных(не разрывных) функций следует выбирать степенные функции вида:
<mrow><mi>D</mi><mo>=</mo><msup><mi>a</mi><mi>x</mi></msup><mo>⋅</mo><msup><mi>b</mi><mi>y</mi></msup><mo>⋅</mo><msup><mi>c</mi><mi>z</mi></msup></mrow>
,
где:
D - интегральный показатель деградации;
a, b, c - индексы-параметры, определяющие степень деградации;
x, y, z - показатели степени чувствительности единичного индекса-параметра в интегральном показателе.
Количественный расчет индексов-параметров и показателей степеней многопараметрической функции приведен в примере конкретной реализации.
Пример реализации способа
Заявленный способ может быть реализован по схеме фиг.5. Функциональная схема устройства содержит космическую систему LandSat (1) на космических аппаратах (2) которой установлены многоканальные радиометры (3), осуществляющие регистрацию отраженного от подстилающей поверхности солнечного потока в полосе сканирования (4), шириной 185 км с пространственным разрешением на пиксель ~30 м, в четырех зональных каналах.
Трассовую съемку запланированных регионов осуществляют по командам от бортового комплекса управления (5) на основе суточной программы управления полетом космического аппарата, передаваемой из центра управления системой (6), по заявкам потребителей.
Результаты измерений записываются в буферное запоминающее устройство (7) и в зонах радиовидимости космического аппарата с наземных пунктов приема (8) (г.Москва, г.Красногорск, г.Новосибирск) сбрасываются по каналу передачи данных (9).
После предварительной обработки информации по служебным признакам (номер витка, время съемки, координаты) на средствах (10) пункта приема информация попадает на сервер хранения данных (11).
Тематическую обработку изображений потребители осуществляют в центре обработки (12), где через устройство ввода и передачи (13) информация из сервера хранения поступает в электронно-вычислительную машину (14) со стандартным набором периферийных устройств: процессор (15), оперативное запоминающие устройство (16), накопитель на магнитных дисках (17), устройство отображения информации (18), печатающие устройство (19), клавиатура (20).
Обработанные снимки участков с расчетными значениями деградации почвенного покрова для пахотных земель региона помещают на сервер хранения (21) результатов обработки.
Состояние почвенных покровов оценивалось на примере участка в Рузском районе Московской области. Обработке подверглись 46 мультиспектральных снимков LandSat 4-5 ТМ с прилегающими окрестностями за период наблюдения с 2005 по 2011 годы. Пространственное разрешение снимков ~30 м/пиксель, сторона квадрата обработанного участка составляет 3030 м.
0
Состояние недеградированного пахотного слоя суглинистых почв Московской зоны оценивается следующими показателями: мощность гумусового горизонта 10-15 см, запасы гумуса порядка 50 т/га, процентное содержание гумуса Н (%) ~3%, индекс засоленности ~0,03, влажностный индекс ~0,24.
Усредненный данные результатов обработки сигналов на всем периоде наблюдений, представлены в таблице 3.
Таблица 3
Усредненные результаты обработки по всему периоду наблюдений
Категория состояния
Степень деградации
Потеря гумуса <maths num="8"><math display="block"><mrow><mfrac><mrow><msub><mi>H</mi><mn>0</mn></msub></mrow><mi>H</mi></mfrac></mrow></math><img file="00000008.tif" he="12" wi="8" img-content="undefined" img-format="tif" /></maths>
Индекс солености
Влагопотери <maths num="9"><math display="block"><mrow><mfrac><mrow><msub><mi>W</mi><mn>0</mn></msub></mrow><mi>W</mi></mfrac></mrow></math><img file="00000009.tif" he="12" wi="8" img-content="undefined" img-format="tif" /></maths>
1
норма
1,0
1,0
1,0
2
низкая
1,4
1,12
1,05
3
средняя
1,62
1,3
1,07
4
высокая
1,75
1,6
1,14
По данным таблицы 3, для получения многопараметрической функции регрессии, записывают следующую систему уравнений:
<mrow><mrow><mo>{</mo><mtable columnalign="left"><mtr><mtd><mi>D</mi><mo>=</mo><mn>2</mn><mo>=</mo><msup><mrow><mo>(</mo><mrow><mn>1,4</mn></mrow><mo>)</mo></mrow><mi>x</mi></msup><mo>×</mo><msup><mrow><mo>(</mo><mrow><mn>1,12</mn></mrow><mo>)</mo></mrow><mi>y</mi></msup><mo>×</mo><msup><mrow><mo>(</mo><mrow><mn>1,05</mn></mrow><mo>)</mo></mrow><mi>z</mi></msup></mtd></mtr><mtr><mtd><mi>D</mi><mo>=</mo><mn>3</mn><mo>=</mo><msup><mrow><mo>(</mo><mrow><mn>1,62</mn></mrow><mo>)</mo></mrow><mi>x</mi></msup><mo>×</mo><msup><mrow><mo>(</mo><mrow><mn>1,3</mn></mrow><mo>)</mo></mrow><mi>y</mi></msup><mo>×</mo><msup><mrow><mo>(</mo><mrow><mn>1,07</mn></mrow><mo>)</mo></mrow><mi>z</mi></msup></mtd></mtr><mtr><mtd><mi>D</mi><mo>=</mo><mn>4</mn><mo>=</mo><msup><mrow><mo>(</mo><mrow><mn>1,75</mn></mrow><mo>)</mo></mrow><mi>x</mi></msup><mo>×</mo><msup><mrow><mo>(</mo><mrow><mn>1,6</mn></mrow><mo>)</mo></mrow><mi>y</mi></msup><mo>×</mo><msup><mrow><mo>(</mo><mrow><mn>1,14</mn></mrow><mo>)</mo></mrow><mi>z</mi></msup></mtd></mtr></mtable></mrow></mrow>
После логарифмирования система степенных уравнений сводиться к линейной алгебраической. Решение системы осуществляют методами Крамера и Саррюса. Определитель системы:
<mrow><mi>Δ</mi><mo>=</mo><mrow><mo>|</mo><mrow><mtable><mtr><mtd><mrow><mi>lg</mi><mn>1,4</mn></mrow></mtd><mtd><mrow><mi>lg</mi><mn>1,2</mn></mrow></mtd><mtd><mrow><mi>lg</mi><mn>1,05</mn></mrow></mtd></mtr><mtr><mtd><mrow><mi>lg</mi><mn>1,62</mn></mrow></mtd><mtd><mrow><mi>lg</mi><mn>1,3</mn></mrow></mtd><mtd><mrow><mi>lg</mi><mn>1,07</mn></mrow></mtd></mtr><mtr><mtd><mrow><mi>lg</mi><mn>1,75</mn></mrow></mtd><mtd><mrow><mi>lg</mi><mn>1,6</mn></mrow></mtd><mtd><mrow><mi>lg</mi><mn>1,14</mn></mrow></mtd></mtr></mtable></mrow><mo>|</mo></mrow></mrow>
Алгебраические дополнения Δx, Δy, Δz получают путем замены коэффициентов при неизвестных свободными членами, т.е. столбцом из lg2, lg3, lg4. Получены следующие цифровые характеристики:
<mrow><mi>x</mi><mo>=</mo><mfrac><mrow><mi>Δ</mi><mi>x</mi></mrow><mi>Δ</mi></mfrac><mo>=</mo><mn>1,9</mn></mrow>
<mrow><mi>y</mi><mo>=</mo><mfrac><mrow><mi>Δ</mi><mi>y</mi></mrow><mi>Δ</mi></mfrac><mo>=</mo><mn>0,5</mn></mrow>
<mrow><mi>z</mi><mo>=</mo><mfrac><mrow><mi>Δ</mi><mi>z</mi></mrow><mi>Δ</mi></mfrac><mo>=</mo><mn>0,3</mn></mrow>
Для каждого пикселя изображения рассчитывают показатель D и формируют интегральную матрицу изображения D(x,y).
Психологически, восприятие образа человеком-оператором происходит на уровне контурного рисунка. Контурный рисунок получают путем вычисления градиента скалярной функции D(x,y) как:
<mrow><mi>g</mi><mi>r</mi><mi>a</mi><mi>d</mi><mi>D</mi><mrow><mo>(</mo><mrow><mi>x</mi><mo>,</mo><mi>y</mi></mrow><mo>)</mo></mrow><mo>=</mo><mfrac><mrow><mo>∂</mo><mi>D</mi></mrow><mrow><mo>∂</mo><mi>x</mi></mrow></mfrac><mo>⋅</mo><mi>i</mi><mo>+</mo><mfrac><mrow><mo>∂</mo><mi>D</mi></mrow><mrow><mo>∂</mo><mi>y</mi></mrow></mfrac><mo>⋅</mo><mi>j</mi></mrow>
[см. пример «Производная по направлению» в учебнике Н.С. Пискунов, «Дифференциальные и интегральные исчисления для ВТУЗов», 5-ое изд., т.1, Наука, М, 1964 г., стр.264-268]
Для получения контурного рисунка выбирают регулярный оператор с апертурой окна 2 на 2 элемента.
i, j
i, j+1
i+1, j
i+1, j+1
Элементы окна связанны по диагоналям (двум взаимно ортогональным направлениям) операцией вычитания. Известны поераторы Собеля, Робертса, Лапласа [см. например Дуда Р.О., Харт П.Е., «Распознавание образов и анализ сцен», перев. с англ., М, Мир, 1976 г, стр.287-288, §7.3 Пространственное дифференцирование].
В частности, оператор Робертса R(i,j) в каждой точке вычисляют как:
<mrow><mi>R</mi><mrow><mo>(</mo><mrow><mi>i</mi><mo>,</mo><mi>j</mi></mrow><mo>)</mo></mrow><mo>=</mo><mrow><mo>|</mo><mrow><mi>D</mi><mrow><mo>(</mo><mrow><mi>i</mi><mo>,</mo><mi>j</mi></mrow><mo>)</mo></mrow><mo>−</mo><mi>D</mi><mrow><mo>(</mo><mrow><mi>i</mi><mo>+</mo><mn>1,</mn><mi>j</mi><mo>+</mo><mn>1</mn></mrow><mo>)</mo></mrow></mrow><mo>|</mo></mrow><mo>−</mo><mrow><mo>|</mo><mrow><mi>D</mi><mrow><mo>(</mo><mrow><mi>i</mi><mo>+</mo><mn>1,</mn><mi>j</mi></mrow><mo>)</mo></mrow><mo>−</mo><mi>D</mi><mrow><mo>(</mo><mrow><mi>i</mi><mo>,</mo><mi>j</mi><mo>+</mo><mn>1</mn></mrow><mo>)</mo></mrow></mrow><mo>|</mo></mrow></mrow>
Выводят на экран точки, для которых величина R(i,j) больше порогового значения, установленного оператором.
Рассмотренный алгоритм реализуют следующей программой.
Программа выделения контуров на изображении
Результаты определения деградации почвенного покрова тестового участка в Московской области показаны в виде контуров деградации на карте схеме, иллюстрируются рисунком фиг.4. Области под номером (1) - водоемы, (2) - области с малой деградацией, (3) - области со средней деградацией, (4) - области с сильной деградацией.
Эффективность заявленного способа характеризуется оперативностью, достоверностью, масштабностью и документальностью получаемых результатов.