Алгоритм расчета степени деградации почвы

В работе предлагается метод оценки степени деградации почв, базирующийся на вегетационном и других индексах подстилающей поверхности. Метод подтвержден на основе наземных наблюдений. Особое внимание уделяется исследованию деградации почвы, вызванной замусориванием территорий. Разработана методика анализа состояния объектов захоронения отходов по данным космических изображений. Получены карты деградации почв для области, включающей Московскую область и также соседние (Тверскую, Ярославскую, Ивановскую, Владимирскую, Рязанскую, Тульскую, Калужскую и Смоленскую области). Приведены результаты работы алгоритма на примере фрагмента этой карты.

Ключевые слова: дистанционное зондирование, линейная регрессия, деградация почвы, индексы подстилающей поверхности, степень деградации почвы.

A new method for assessment of soil degradation is proposed. The method is based on vegetation and other surface indexes. The method confirmed by ground observations. The special attention is given to study of soil degradation caused by anthropogenic pollution of territories. Developed method for assessment of condition of waste sites based on remote sensing data. Computed soil degradation maps for the following regions: Moscow, Tver, Yaroslavl, Ivanovo, Vladimir, Ryazan, Tula, Kaluga and Smolensk.

Keywords: remote sensing, linear regression, soil degradation, underling surface indexes, extent of soil degradation.

Введение

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

1. Заболачивание – избыточное увлажнение почвы, характеризующееся повышенной кислотностью и низкой плодородностью почвы, выходом на поверхность стоячих или проточных вод, но без постоянного слоя воды на поверхности;
2. Технологическая деградация – в результате длительного использования почвы под жилые или промышленные застройки, разработки полезных ископаемых;
3. Замусоривание – покрытие поверхности почвы слоем мусора, препятствующим нормальному функционированию почвы;
4. Засоление – накопление в верхних горизонтах почвы легкорастворимых солей, препятствующих развитию большинства растений;
5. Опустынивание – характеризуется иссушением земли, увяданием растительности, снижением связности почвы.

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

Для оценки деградации в работе используются n мультиспектральных изображений серий Landsat 4-5 TM: .

Изображения расположены в хронологическом порядке съемки. Каждое изображение состоит из 7 каналов, а каждый канал представлен в виде вектора-столбца длины m=s₁s₂, где s₁ и s₂ – число строк и столбцов матрицы канала изображения. Т.о. изображение можно представить в виде матрицы: .

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

Перечислим используемые индексы состояния растительности [2]: индекс условий роста VCI; индекс солености SI; индекс тепла NT; индекс влаги NW; индекс света NL.

Индекс VCI рассчитывается по формуле:

, (1)

где q­­₃ - коэффициент спектральной яркости пикселя на третьем канале, q₄ – коэффициент спектральной яркости пикселя на четвертом канале.

Индекс SI рассчитывается по формуле:

, (2)

где q₁ – коэффициент спектральной яркости пикселя на первом канале.

Индекс NT рассчитывается по формуле:

, (3)

где q₆ – коэффициент спектральной яркости пикселя на шестом (тепловом) канале, T – нормированная температура поверхности.

Индекс NW рассчитывается по формуле:

, (4)

где W – количество осадков в области пикселя в день съемки (в мм РТ. С.).

Индекс NL рассчитывается по формуле:

, (5)

где l – продолжительность светового дня (в часах), b – бал облачности.

Если T определяется для каждого пикселя, то W, l и b постоянно для некоторой области и определяется по истории метеоданных. На метеостанциях, распределенных по территории РФ, измеряются погодные параметры с некоторой частотой и вводятся в базу метеоданных. Каждая метеостанция контролирует свою область, для каждого пикселя которой каждый погодный параметр считается постоянным. Например, продолжительность солнечного сияния L зависит от l и b. При этом продолжительность светового дня l является в общем случае функцией от времени и от географических координат (широты и долготы). Однако можно для каждого региона получить свою функцию l(t). Так, для Московского региона самый длинный день t₁ 21 июня и его продолжительность составляет 17 часов и 33 минуты. Самый короткий день t₂ 22 декабря и его продолжительность составляет 4 часа и 12 минут. Продолжительность светового дня L0 для Московского региона меняется во времени по гармоническому закону, изображенному на рис. 1. L0 выражена в долях от продолжительности суток (24 часа).

Рисунок 1 – Продолжительность светового дня от времени года для Московского региона

Бал облачности b – концентрация облаков на небе, равно как и количество осадков W, находится в базах метеоданных для соответствующей станции (на один регион станций может быть несколько). Как и L бал облачности b лежит в диапазоне от 0 до 1. Параметр b редко присутствует в базах метеоданных для российских регионов и является платным ресурсом. Поэтому можно положить для всех снимков территории РФ b=1.

Приведенные индексы состояния растительности определяются для каждой даты. Например, индексы NT, NW и NL задают, насколько хорошо растительный покров реагирует на тепло, влагу и свет соответственно в данный момент времени. Т.о. мы имеем дело с системой управления, входными воздействиями которой являются 3 необходимых вида питания растительности (тепло, влага и свет), а реакциями на входные воздействия – качество растительного покрова. Одному входному воздействию соответствует одно значение индекса, временному ряду входных воздействий – соответственно временной ряд индексов. Значения всех индексов лежит в диапазоне от 0 до 1, и чем больше значение индекса, тем хуже качество растительности.

По временному ряду индексов состояния растительности определяются степени деградаций почвы (СДП). Для этого проводятся процедуры темпоральной обработки. Допустим, для некоторого пикселя поверхности земли мы получили n временных отчетов значений индекса солености SI. Решаем задачу линейной регрессии и строим линейный тренд для этого временного ряда. Полученное уравнение линейного тренда: , где k – коэффициент наклона, b – свободный член. Число, пропорциональное k – значение степени деградации по индексу солености для данного пикселя. Пикселей исследуемой территории много – столько, сколько произведение числа строк на число столбцов матрицы изображения (m пикселей). Поэтому целесообразно провести процедуру множественной линейной регрессии и найти значения степеней деградации для всех пикселей изображения одновременно.

Для оценки СДП используются различные индексы состояния растительности [3]:

1. Индекс солености SI (получаем СДП по индексу солености RSI);
2. Индекс тепла NT (получаем СДП по индексу тепла RNT);
3. Индекс влаги NW (получаем СДП по индексу влаги RNW);
4. Индекс света NL (получаем СДП по индексу света RNL).

СДП рассчитывается по формуле:

.

СДП лежит в диапазоне от -1 до 1. Если СДП положительна, то присутствует деградация, причем чем больше СДП, тем больше деградация почвы и растительности. Если СДП отрицательна, то деградации нет, происходит восстановление почвы и растительности по данному индексу. Если RNT₁<RNT₂<0, то растительность с RNT₁ лучше реагирует на тепло, чем растительность с RNT₂, что говорит о лучшем состоянии первой.

Значения постоянных a и b для разных степеней деградации почвы показаны в таблице 1. СДП означает скорость деградации, изменение воспроизводства растительного покрова во времени. В таблице 2 показаны критерии степени деградации почвы [4].

Значения СДП отображаются в виде бинарных изображений СДП. Пиксели этих изображений принимают два значения: 0 или 1. Например, на изображении слабой деградации значение пикселя равно 1, если СДП этого пикселя лежит в диапазоне (0…0.2], и 0 – в противном случае.

Изображения СДП не всегда удобно воспринимать. Поэтому вводятся функциональная зависимость площади деградации от степени деградации S(R). Каждая точка с координатами [r s] функции S(R) задает площадь деградации s поверхности земли, имеющую степень деградации не меньше r. Т.о. функция S(R) монотонно убывает и имеет область определения 0..1.

Таблица 1

Значения коэффициентов a и b

Таблица 2

Критерии степеней деградации почвы

Блок-схема алгоритма расчета СДП показана на рис. 2. Стадии алгоритма:

1. Предобработка;
2. Получение индексов состояния растительности и заполнения временного хранилища данными индексов SI, NT, NW и NL;
3. Получение СДП по индексам солености, тепла, влаги и света по данным временного хранилища и заполнение ими хранилища данных;
4. Вывод результатов работы алгоритма на экран.

Входными данными алгоритма являются мультиспектральные изображения серий Landsat 4-5 TM пространственного разрешения 30 м с 7 спектральными каналами. Снимки берутся малооблачные за долгосрочный период наблюдений – от 10 лет. Выходной информацией являются изображения степеней деградации почвы, зависимостей площади деградации от степени деградации в графической и табличной форме.

Шаг 1. Ввод данных. На специализированных сайтах свободного доступа по скачиванию космических изображений настраиваются параметры скачивания космических изображений:

1. Контрольная область: задается многоугольником вводом географических координат вершин с клавиатуры или нажатиями кнопки мыши на карте Земли;

2. Число скачиваемых изображений: чем точнее должен быть расчет, тем больше предельное число скачиваемых изображений;

3. Облачность, т.е. концентрация облаков на снимке: вводится, как можно меньшей;
4. База данных космических изображений: Landsat 4-5 TM;
5. Время суток съемки: дневные снимки.

Рисунок 2 – Блок-схема алгоритма расчета степени деградации почвы

Шаг 2. Получение снимков и метаданных. С помощью программы автоматического скачивания космических изображений осуществляется процедура скачивания изображений на компьютер или сервер. Метаданные обычно присутствуют в папках с изображениями. В случае отсутствия метаданных с сайта скачивается файл расширения *.xls (файл Microsoft Excel) с основными метаданными всех скаченных снимков.

Шаг 3. Геопривязка. В программе ENVI+IDL осуществляется процедура географической привязки всех скаченных изображений либо в интерфейсном режиме, либо программированием. Получаем многоканальное изображение I, вектор номеров снимков S и вектор номеров каналов C:

,

где I_i – вектор-столбец канала C_i мультиспектрального изображения S_i. Номера мультиспектральных каналов C_i лежат в диапазоне от 1 до 7 (для снимков серий Landsat 4 TM и Landsat 5 TM).

Шаг 4. Получение маски облачности. Маска облачности получается для каждого мультиспектрального изображения с помощью алгоритма автоматической оценки облачности ACCA. На выходе получаем бинарное изображение со значениями: 0, если облачности нет и 1, если облачность есть.

Шаг 5. Получение участков изображений и временных данных. Режем геопривязанное изображение I на квадраты заданной стороны (например, 100 пикселей). Получаем участки изображений с p слоями. В каждом участке удаляем слои, для которых в маске облачности существует хотя бы одно значение, равное 1. Также удаляем слои, в которых существует хотя бы одно значение коэффициента спектральной яркости пикселя, равное 0.

По метаданным и векторам S и C строим данные слоев для каждого участка:

,

где s и c – вектор-столбцы номеров снимков и каналов мультиспектральных изображений, t1 и t2 – вектор-столбцы годов и дней годов съемки. Элементы векторов s, c и t расположены в соответствии с порядком слоев участка. Сохраняем данные слоев в текстовых файлах.

Шаг 6. Получение метеоданных. Метеоданные можно получить на специализированных сайтах архивов погоды. Обычно они сохраняются в виде матриц, столбцы которых – требуемые для задачи параметры погоды и временные отчеты. Нам нужны данные о количестве осадков и балах облачности. Сохраняем эти данные для метеостанций, контролирующих область исследования.

Шаг 7. Получение данных о количестве осадков и балах облачности. Зная метеостанцию, контролирующую участок, временные ряды количеств осадков и балов облачности, полученные на шаге 6, а также данные слоев участков, для каждого участка получаем векторы:

,

где W_i – количество осадков, b_i – бал облачности для слоя i. Можно положить b_i=1.

Шаг 8. Расчет индексов состояния растительности. Для всех участков по формуле (1) для всех пар 3-их и 4-ых каналов (хранятся в векторах s и c) рассчитываем VCI. Для всех участков по формуле (2) для всех пар 1-ых и 3-их каналов рассчитываем SI. Для всех участков, зная VCI, по формуле (3) для всех 6-ых каналов рассчитываем NT. Для всех участков, зная VCI и VW, по формуле (4) для всех 6-ых каналов рассчитываем NW. Для всех участков, зная VCI и L, по формуле (5) для всех 6-ых каналов рассчитываем NL. L находим по формуле (5), зная Vb и зависимость продолжительность светового дня от времени года для данной контрольной области.

Шаг 9. Расчет степеней деградации почвы. Для каждого участка, зная вектор t и временной ряд индексов солености, находим СДП по индексу солености RSI. Для каждого участка, зная вектор t и временной ряд индексов тепла, находим СДП по индексу тепла RNT. Для каждого участка, зная вектор t и временные ряды индексов влаги и света, находим СДП по индексам влаги и света RNW и RNL.

Шаг 10. Вывод данных. Выводим 2 вида данных:

1. Изображения слабой, средней и сильной степеней деградации почвы для контрольной области, сложенной из участков;
2. Функциональные зависимости площадей деградаций от степеней деградаций для контрольной области.

Алгоритм расчета СДП был проведен для территории Московского региона и окружающих его регионов (контрольная область). Входными данными алгоритма являются 1033 космических мультиспектральных изображения серий Landsat 4-5 TM за период 1983-2012 годы. Выходной информацией являются:

1. (а)	(б)	(в)
   Рисунок 3 – а) Изображение степени деградации по индексу тепла в общем случае; б) Изображение участков сильной деградации по индексу тепла; в) Изображение участков среднего восстановления по индексу тепла

   Карты деградаций почвы по индексам солености, тепла, влаги и света контрольной области с разрешением 30 м;
2. Зависимости площади деградации почвы от степеней деградаций в табличной и графической форме;
3. Изображения слабой, средней и сильной деградаций.

Покажем результаты работы алгоритма на фрагменте контрольной области. В качестве фрагмента рассмотрим полигон ТБО Кучино Московского региона (Балашихинский район) и его ближайшую окрестность.

На рис. 3а показано изображение СДП по индексу тепла. Изображения участков сильной деградации почвы – рис. 3б, среднего восстановления почвы – рис. 3в. На рис. 3б и 3в области деградации (восстановления) отмечены черным. Видно, что почва сильно деградаровала на участке складирования, причем на тех картах складирования участка, на которых происходит в настоящий момент захоронение отходов. Среднее восстановление почвы наблюдается в основном вблизи водоема (р. Пехорка). Сильного восстановления не наблюдается ни для одной СДП. По изображениям всех СДП, то можно заключить:

1. Наибольшей деградации по всем СДП подвержен участок складирования;
2. Отмечена засолоненность и заболачиваемость участка складирования и его окрестности, в особенности прилегающего поля;
3. Лес имеет свойство экранировать влияние загрязнения на почву, тогда как поле более чувствительно к деградации;
4. Если деградация наблюдается по одному индексу, то по остальным в большинстве случаев она также наблюдается.

Рисунок 4 – Значения площади деградации S почвы, имеющую степень деградации не меньше R: (1) – по индексу солености, (2) – по индексу тепла, (3) – по индексу влаги, (4) – по индексу света

На рис. 4 показаны графики зависимости площади деградации территории в гектарах от степени деградации: 1 – для индекса солености, 2 – для индекса тепла, 3 – для индекса влаги, 4 – для индекса света. Чем быстрее убывает характеристика, тем лучше состояние почвы. Однако сравнивать разные характеристики (например, по индексу солености и тепла) нельзя, поскольку каждая из них оказывает свое влияние на состояние растительности. Из графиков видно, что площадь сильной деградации больше 20 га.

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

ЗАКЛЮЧЕНИЕ

1. В работе показана возможность оценки степени деградации почвы по космическим снимкам на примере снимков Landsat 4-5 TM.

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

3. Исследование деградируемости почвы позволяет строить прогнозы состояния растительности и дает возможность принять необходимые меры по улучшению или не ухудшению плодородности земли и ее воспроизводства.

1. Методические рекомендации по выявлению деградированных и загрязненных земель: Прил к письму Ком. РФ по земел. Ресурсам и зелеустройству от27.03.1995 №3-15/582;

2. Crippen R.E. (1990) “Calculationg the vagitation index faster”, Ramote sensing of Environment, vol 34., pp 71-73;

3. Шовенгердт Р.А., Дистанционное зондирование. Модели и методы обработки изображений, Москва: Техносфера, 2010. – 560 с., 32 с. Цв. Вкл. ISBN 978-5-94836-244-1;
4. Известия ЮФУ. Технические науки. №6 (131). 2012 г., стр. 131-134, В.Г. Бондур, А.Б. Мурынин, А.А. Рихтер, М.А. Шахраманьян, «РАЗРАБОТКА АЛГОРИТМА ОЦЕНКИ СТЕПЕНИ ДЕГРАДАЦИИ ПОЧВЫ ПО МУЛЬТИСПЕКТРАЛЬНЫМ ИЗОБРАЖЕНИЯМ»;
5. Щербенко Е.В., Мониторинг засухи по данным космических съемок, 4 Всероссийская открытая конференция «Современные проблемы дистанционного зондирования Земли из космоса», Москва, ИКИ РАН, 13-17 ноября 2006 г., Сборник тезисов конференции.