CC BY
171
РАСПОЗНАВАНИЕ ПАХОТНЫХ ЗЕМЕЛЬ НА ОСНОВЕ МНОГОЛЕТНИХ СПУТНИКОВЫХ ДАННЫХ СПЕКТРОРАДИОМЕТРА MODIS И ЛОКАЛЬНО-АДАПТИВНОЙ КЛАССИФИКАЦИИ
Барталёв С.А., Егоров В.А., Лупян Е.А., Плотников Д.Е., Уваров И.А.
Учреждение Российской академии наук Институт космических исследований РАН
Аннотация
Разработан метод распознавания пахотных земель на основе многолетних временных рядов данных дистанционного зондирования Земли, получаемых спектрорадиометром MODIS со спутников Terra и Aqua. Метод предполагает вычисление по спутниковым данным набора признаков распознавания, учитывающих особенности сезонной и межгодовой динамики спектрально-отражательных характеристик используемых пахотных земель, отличающих их от других категорий сельскохозяйственных угодий и естественной растительности. Выявление пахотных земель выполняется с использованием алгоритма локально-адаптивной обучаемой классификации, учитывающей пространственную вариабельность значений признаков распознавания.
Использование метода позволило создать карту пахотных земель для всей территории России. Валидация полученной карты, выполненная на основе оптимума Парето с использованием опорных данных для тестового региона, позволила оценить точность метода и возможности его дальнейшего улучшения.
Ключевые слова: дистанционное зондирование, спутниковая спектрорадиометрия, признаки распознавания, обучаемая классификация, мониторинг сельскохозяйственных земель.
Введение
Оценка состояния и динамики земного покрова на основе дистанционного зондирования со спутников часто сопряжена с необходимостью распознавания объектов поверхности на основе данных измерений отражённого или собственного излучения. Для этого, как правило, используются методы обучаемой или не-обучаемой классификации в пространстве признаков, обеспечивающих разделимость типов распознаваемых объектов [1]. При этом стандартные методы классификации, разработанные применительно к данным непространственной природы, не позволяют в явном виде учитывать географическую изменчивость значений признаков распознавания, вызванную особенностями физико-химических, морфологических и структурных свойств объектов, часто подчиняющихся, в случае растительного покрова, закономерностям широтной зональности и высотной поясности.
В этой связи эффективность применения стандартных методов классификации для обработки данных дистанционного зондирования, как правило, снижается с увеличением размеров изучаемой территории. Так, будучи достаточно эффективными при обработке данных на относительно небольшие участки, указанные методы становятся практически не применимыми на глобальном уровне без принятия мер по снижению изменчивости значений признаков распознавания. К числу наиболее часто используемых мер такого рода относится стратификация территории с выделением однородных по некоторым критериям регионов (страт), отличающихся сниженной вариабельностью значений признаков распознавания. Несмотря на то, что данный подход позволяет повысить точность классификации данных спутниковых наблюдений, его недостатком является возможная рассогласованность результатов на границах смежных страт.
Локально-адаптивный алгоритм глобального картографирования LAGMA (Locally-Adaptive Glo-
bal Mapping Algorithm) [3] в силу наличия внутренне присущего ему механизма учёта географической изменчивости признаков распознавания свободен от указанных недостатков и может эффективно использоваться для классификации данных спутниковых наблюдений больших территорий.
Данные установленного на спутниках Terra и Aqua спектрорадиометра MODIS обеспечивают возможность картографирования растительного покрова на глобальном уровне и находят, в частности, широкое применение для мониторинга сельскохозяйственных земель [4,5].
Разработанный ранее в Институте космических исследований Российской академии наук (ИКИ РАН) метод выявления используемых пахотных земель на основе многолетних временных рядов данных MODIS нашёл применение при проведении Всероссийской сельскохозяйственной переписи 2006 года для получения независимой объективной информации [6]. Созданная на основе данного метода карта пахотных земель с пространственным разрешением 250 м вошла в состав информационного обеспечения Системы дистанционного мониторинга земель агропромышленного комплекса (СДМЗ АПК) [4].
Вместе с тем, дополнительные углубленные исследования показали наличие возможностей развития метода выявления пахотных земель по данным MODIS за счёт комбинированного использования набора информативных спектрально-временных признаков распознавания [7] и алгоритма локально-адаптивной классификации земного покрова [3].
Описанный в настоящей статье метод, основанный на использовании многолетних рядов данных MODIS и алгоритма LAGMA, обеспечивает возможность картографирования пахотных земель России с достаточной для решения многих практических задач точностью.
Алгоритм локально-адаптивной классификации Использование методов обучаемой классификации, позволяющих относить объекты классам некоторого множества, заданным своими описаниями в выбранном пространстве признаков, в общем случае предполагает следующую последовательность логических шагов [1]:
- построение статистических описаний (сигнатур) классов в пространстве признаков на основе репрезентативной выборки характерных эталонов;
- классификация объектов на основе значений их признаков и полученных описаний классов. Распознавание типов земного покрова по данным
спутниковых наблюдений основано, как правило, на использовании в качестве признаков значений спектральной яркости объектов (или полученных на их
основе различных производных признаков). Часто используемый статистический подход к классификации исходит из предположения нормального распределения значений признаков, оценки параметров которого формируют сигнатуры классов. Так, сигнатура каждого к-го класса включает в себя параметры вектора средних значений ик и ковариационной матрицы Еик признаков.
В основу алгоритма локально-адаптивной классификации положен подход, предполагающий на первом этапе формирование для рассматриваемой территории пространственного распределения сигнатур на основе репрезентативной совокупности эталонных объектов (обучающей выборки) с известной принадлежностью к одному из классов заданного множества. Блок-схема алгоритма приведена на рис. 1.
с
Начало
э
Выбрать текущий обрабатываемый узел О{р0-До)
Выбрать текущий обрабатываемый класс к
Вычислить Зк(РоЯо) Ск(РоШ Кк(РоЯо)
Выбрать текущий обрабатываемый пиксел С(х0;у0)
Определить ближайший к пикселу узел С&оЛо)
На основе и^{роЩо) вычислить р(Р(х-у)\щ) для всех классов
Классифи пиксел шровать р(Х;у)
( Конец )
Рис 1. Блок-схема алгоритма локально-
В качестве источников данных для формирования обучающей выборки часто выступают существующие тематические карты, данные наземных обследований, результаты визуальной интерпретации спутниковых изображений, знания экспертов.
Пространственное распределение локализованных сигнатур описывается их значениями в узлах G(p, q) регулярной прямоугольной сети с шагом d, где р и q - порядковые номера узлов по осям х и у соответственно (рис. 2). Обучение классификатора, включающее в себя два описанных ниже логических этапа, направлено на оценку в узлах G(p, q) параметров локализованных сигнатур ик (р, q) и Ек (р, q) для каждого к-го класса.
■адаптивной классификации LAGMA
Первоначально для каждого узла G(p, q) на основе расположенных в границах соответствующей клетки эталонных пикселов вычисляются величины: Б'к (р; д) - сумма значений /'-го признака к-го класса; С.' (р; д) - сумма произведений значений /-го иj-m признаков к-го класса;
N. (р, д) - количество эталонных пикселов к-го класса.
Указанные величины используются для оценки элементов Соу.. ' (р; д) ковариационной матрицы £к (р, д) на основе следующего выражения:
Со^' (р; д) = - Бк(р; д) ^(р; д).
кУИЧ' N. (р; д) N. (р; д) N. (р; д)
л й 2с1 „М ,4(1 5й х р=0 р=1 р=2 р=3 р=4
Рис 2. Локализация сигнатур классов в узлах регулярной сети
Оценка элементов Соук1 (р; q) ковариационной матрицы и средних значений признаков на основе (р; q) и Ык (р; q) позволяет получить параметры сигнатур ик (р; q) и Ек (р; q ).
Как видно из рис. 2, обучающая выборка эталонных пикселов, как правило, имеет пространственно неравномерное распределение. При этом в случае отсутствия эталонных пикселов к-го класса в окрестности узла G(р0; q0) вычисление параметров ик(р0;q0) и Ек(р0; q0) оказывается невозможным. Кроме того, при малых значениях Ык (р0;q0) возрастает влияние случайных ошибок в обучающей выборке. Для учёта этого фактора метод предусматривает задание порога репрезентативности Т, характеризующего минимально допустимое количество эталонных пикселов для оценки локальных сигнатур классов. Это приводит к появлению узлов, не обеспеченных значениями параметров сигнатур некоторых классов и, как следствие, к необходимости проведения второго этапа обучения классификатора.
На втором этапе для каждого не преодолевшего порог репрезентативности (Ык(р0, д0)<Т) узла 0(р0, q0) проводится обработка данных соседних клеток. Количество используемых соседних клеток зависит от величины Т, числа эталонных пикселов в клетке Nк(р0, д0) и соседних клетках Ык (р0 + Ар; q0 + Лq). Кроме того, методом предусмотрено задание параметров Lmin и Lnax, ограничивающих число используемых соседних клеток снизу и сверху соответственно.
Число используемых соседних клеток итеративно увеличивается, начиная с Lmin до величины, соответствующей достижению порога репрезентативности Т. Если порог Т не был преодолён при достижении числа ближайших клеток значения Lmax, то сигнатура узла G(р0; q0) для класса к считается несуществующей.
Расширение анализируемой для оценки параметров локальной сигнатуры класса окрестности осуществляется дискретно путём последовательного включения соседних клеток, находящихся на одинаковом удалении от узла G( р0; q0), и обобщения полученных на первом этапе обучения величин
( р0 +аР; qo +Аq):
ск1 (р0 + ар; qo + ад)
Ык (р0 + Ар; q0 + Лq). При этом новые характеристики узла Б'к*(р0; qo), С'к1>(р0; qo) и ык*(р0;д0)вы-числяются как суммы соответствующих величин, полученных для соседних узлов на первом этапе обучения.
По результатам второго этапа обучения, для класса к в узлах, для которых справедливо выражение Мк (р; q) > Т , определяются с использованием выражения аналогичного (1) параметры сигнатур ик*(р;q) и Ек*(р;q) для последующей обучаемой классификации.
В основу алгоритма локально-адаптивной классификации могут быть положены различные решающие правила, включая методы максимального правдоподобия, минимального расстояния или параллелепипеда [1]. К настоящему времени в составе разработанного программного комплекса LAGMA реализован алгоритм классификации на основе метода максимального правдоподобия.
В соответствии с решающим правилом максимума правдоподобия, пиксел Р(х; у) относится ко множеству ю, пикселов класса I в том случае, если для всех к=1, 2.. .т выполняется условие:
р(ю,)р(Р(х;у) | ю,) > р(ак)р(Р(х;у) | Юк), (2)
где р(ю,) и р(юк) - априорные вероятности классов I и к;
р(Р(х;у) | ю,) и р(Р(х;у)| юк) - плотности вероятности отнесения пиксела Р(х;у) ко множеству пикселов класса I и множеству пикселов класса к.
В свою очередь, плотность вероятности определяется по формуле:
р (Р(х; у)| Юк) =
ехр З- 2 (В(х; у) - Пк* (р; д))Т Ек * (р; д)-1 (В( х; у) - Пк* (р; д)) ^
п 1
' * X '
(3)
(2р)2 ек (р,д)
где В(х;у) - вектор признаков пиксела Р(х;у); п - число признаков.
При классификации используются параметры
к (р; д) и Ек*
ик (р; д) и Ек*(р; д) локализованных сигнатур, вы-
и
численных на втором этапе обучения. Для классификации пиксела P(x; y) используются сигнатуры ближайшего узла G(p; q), порядковые номера которого в строках и столбцах регулярной сетки с шагом d определяются по формулам p = x / d и q = y / d .
Наряду с сигнатурами, характеризующими локализованные значения признаков, при классификации используются априорные вероятности, полученные на основе данных об ареалах распространения различных типов земного покрова в пределах рассматриваемой территории. В решающем правиле (2) априорная вероятность класса используется в качестве коэффициента, определяемого в виде дробной величины в диапазоне от 0 до 1 для каждого класса и каждого пиксела территории. Априорные вероятности могут быть получены в результате обобщения обучающей выборки таким образом, что в окрестности эталонных пикселов заданного класса её значения максимальны и снижаются по эмпирически подобранным правилам по мере удаления.
Данные спектрорадиометра MODIS
Спектрорадиометр MODIS (Moderate Resolution Imaging Spectroradiometer) установлен на борту спутников Terra и Aqua. Ширина полосы захвата прибора составляет 2330 км при угле горизонтальной развёртки ± 55о. Прибор ведёт съёмку в 36 каналах видимого и инфракрасного (ИК) диапазонов длин волн с пространственным разрешением 250 м, 500 м и 1000 м в зависимости от канала.
Данные MODIS распространяются Геологической службой США (http://lpdaac.usgs.gov/main.asp). Для разработки описанных в настоящей работе методов была сформирована база данных MODIS, в которую вошли стандартные продукты M0D09 за период 2000-2010 годов. Набор данных включает измерения коэффициента спектральной яркости (КСЯ) в красном (620-670 нм) и ближнем ИК (841876 нм) каналах с пространственным разрешением 250 м. M0D09 содержат также измерения КСЯ в каналах голубого (459-479 нм) и среднего ИК (16281652 нм) диапазонов с пространственным разрешением 500 м. В дополнение к измерениям КСЯ, данные MOD09 содержат информацию об угловых характеристиках Солнца и прибора в момент съёмки. Данные MOD09 представлены в синусоидальной картографической проекции и разбиты на гранулы, регулярно покрывающие поверхность Земли.
Предварительная обработка данных MODIS
Предварительная обработка данных MODIS направлена на снижение влияния различных мешающих факторов и, следовательно, повышение эффективности тематической обработки данных. Согласно [8] предварительная обработка данных MODIS предусматривает выполнение следующих последовательных этапов:
1) маскирование пикселов по их значениям угловых характеристик освещения и наблюдения;
2) детектирование пикселов, подверженных влиянию снежного и облачного покровов;
3) детектирование пикселов, соответствующих участкам теней от облаков;
4) статистическая фильтрация временных рядов данных;
5) построение композитных изображений максимально свободных от влияния мешающих факторов.
Для выбора наблюдений, пригодных для дальнейшей обработки по геометрическим условиям освещения и наблюдения, используется совокупность пороговых критериев ПА > 400 и Б1А > 800, где ПА - зенитный угол наблюдения, 81А - зенитный угол Солнца.
Детектирование участков с наличием облачного и снежного покровов выполняется с использованием данных измерений КСЯ в голубом R3 (459-479 нм) и среднем ИК R6 (1628-1652 нм) каналах MODIS. В основе данного алгоритма детектирования лежит использование нормализованного разностного индекса снега NDSI, определяемого по формуле (3):
N081 = , (4)
R3 + R6
и априорных знаний об отражательных свойствах поверхности.
Примем на данном этапе, что каждый пиксель относится к одному из четырёх классов: облачность, полупрозрачная облачность, снег, свободная от влияния облачности и снега подстилающая поверхность (чистая поверхность).
Рассмотрим двумерное пространство признаков R3 и N081 (рис.3), которое разобьём на четыре области следующим образом:
«снег», если R3 > 0,05 и N081 > 0,1; «облачность», если R3 > 0,05 и -0,2 < N081 < 0,1; «полупрозрачная облачность», если R3 > 0,05 и - 0,35 < N081 < -0,2;
«чистая поверхность» во всех остальных случаях. Пиксели, окружающие классы «облачности» и «полупрозрачной облачности», также относятся к этим классам, если их значения R3 не меньше значений крайних пикселов, принадлежащих к классу.
Используя предполагаемое значение максимальной высоты облаков, на следующем шаге рассчитывается местоположение соответствующих им теней с использованием данных об углах наблюдения и положения Солнца (рис. 4).
Поскольку точные данные о высоте облаков отсутствуют, в качестве максимально возможной соответствующей величины было выбрано значение = 12 км, а алгоритм включил в себя построение зоны, заведомо содержащей подлежащие выделению пикселы. Если ввести прямоугольную декарто-ву систему координат с началом в данном соответствующем облаку пикселе с осью Ох, направленной на север, и осью Оу, направленной на восток, тогда радиус-вектор смещения тени на изображении относительно облака высотой Н задаётся в этой системе следующими координатами:
(5)
х = Н - ^(Р^(5))
у = Н (яп(у^(6) - sin(P)tg(5)) , где у - азимутальный угол наблюдения, 6 - зенитный угол наблюдения, р - азимутальный угол Солнца, д - зенитный угол Солнца.
ЫОБ! Чистая ^ поверхность Снег
Облачность ЯЗ
-0,2- Чистая поверхность -0,35- 0,05
Полупрозрачная облачность
Чистая поверхность
Рис. 3. Пороговое разделение пикселов «чистой поверхности», облачности и снежного покрова в двумерном пространстве R3-NDSI
Прибор
\ / Солнце
Рис. 4. Геометрическое расположение линии тени от облака на наблюдаемой спутниковым прибором поверхности Земли
Как правило, геометрическая линия тени, описываемая формулой (5), захватывает не только теневые пиксели, но и чистую поверхность, что обуславливает необходимость дополнительного анализа пикселов, отнесённых к классу геометрической тени. Так, при рассмотрении линии тени можно определить реальную границу затенённых участков по скачкообразному росту (рис. 5) отражательной способности Я2 в ближнем ИК канале (841-876 нм) MODIS.
Этап статистической фильтрации временных рядов направлен на выявление неверно выделенных теней, возникающих за счёт имеющего место ошибочного отнесения участков снега с недостаточно высоким значением индекса NDSI к классу облачности. Ложной тенью считается пиксел, не попадавший ни разу в класс «чистая поверхность» в течение двадцати дней (десять дней до даты начала интервала и десять дней после) при условии выпол-
нения следующего условия для пикселов геометрической тени:
R1(©*, t) > МК1 (©*, t) + стж (©*, t), (6)
где ст^(©*, t) - стандартное отклонение от среднего значения М^(©*, t).
Л в г
КСЯ
7 10 13 16 19 22 25 28 31 номер пиксела вдоль линии тени ---медианный фильтр -исходный ряд
Рис. 5. Профиль значений КСЯ в ближнем ИК канале MODIS по линии тени от облака. Зоны линии тени: А - зона остаточной облачности, Б - область тени, В - область избыточно выделенной тени
Фильтрация остаточного влияния мешающих факторов и детектирование сбойных пикселов осуществляется на основе анализа временных рядов данных MODIS за указанный выше двадцатидневный интервал времени. Отделение «зашумлённых» пикселей происходит по критерию превышения удвоенного стандартного отклонения ст^©*, t) от среднего значения Мм(©*, t):
^6(©*, t) - Мм(©*, t) | > 2 ^6(©*, t). (7)
В результате выполнения перечисленных выше шагов фильтрации данных MODIS формируется набор «масок» с указанием статуса каждого пиксела, используемых для последующего осреднения очищенных данных за заданные интервалы времени.
Предварительная обработка данных MODIS позволяет в значительной мере компенсировать влияние мешающих факторов и даёт возможность построения на основе ежедневных данных композитных изображений за различные промежутки времени с редуцированным влиянием облачного и снежного покровов, теней и аппаратурных шумов. В данной работе используются еженедельные композитные изображения для учёта сезонной динамики спектрально-отражательных характеристик земного покрова. На основе значений КСЯ в каналах композитных изображений проводится расчёт значений спектрального вегетационного индекса (ВИ). При решении задач мониторинга пахотных земель на больших территориях с меняющимися ярко-стными характеристиками почвенного покрова эффективным является использование почвенно-адаптивных ВИ. Перпендикулярный вегетационный индекс РУ1 в значительной мере независим от яркости почвенного покрова и тесно коррелирует с объёмом зелёной биомассы и концентрацией хлорофилла в растениях. Это делает возможным использование РУ1 для мониторинга сельскохозяйственной растительности и пахотных
земель больших территорий. Значение РУ1 рассчитывается в двумерном пространстве значений КСЯ красного и ближнего ИК диапазонов как евклидово расстояние от данной точки до линии почв. С использованием полученного в работе [7] уравнения для линии почв, РУ1 вычисляется по формуле
PVI(R1,R2) = -0,74^ + 0,67R2 -0,034, (8) где R1 и R2 соответствуют измерениям КСЯ в красном и ближнем ИК каналах. На основе рассчитанных значений РУ1 формируются временные серии данных с последующим сглаживанием временных рядов на основе скользящих статистических фильтров для снижения влияния случайного шума и заполнения пропусков данных.
Спектрально-динамические признаки и метод распознавания используемых пахотных земель Некоторые типы земного покрова, к числу которых относятся и используемые пахотные земли (ИПЗ), как правило, невозможно надёжно распознавать на основе одномоментных спутниковых данных или даже регулярных наблюдений в течение одного сезона. Так, опираясь на анализ только лишь сезонной динамики КСЯ, в зависимости от моментов наблюдений и текущего состояния пашни, участки ИПЗ часто можно ошибочно отнести к классу открытых почв или лугов (степей). Наличие многолетних рядов спутниковых наблюдений открывает возможности разработки новых информативных признаков распознавания ИПЗ и других типов земного покрова.
В данной работе используются признаки, построенные на основе многолетних рядов РУГ, учитывающие характерные отличия между естественной и сельскохозяйственной растительностью ввиду наличия на участках последней севооборота.
РУ1
0,35 0,300,25 0,200,150,10 0,05
0
.1
и.
Приведённые в качестве примера на рис. 6 временные ряды РУГ получены для двух участков пашни и естественной растительности, находящихся в близких почвенно-климатических условиях. Представленный рисунок показывает, что многолетние ряды РУГ рассматриваемых классов имеют характерные особенности, позволяющие разработать признаки распознавания ИПЗ. В частности, из рисунка видно, что культурная растительность начинает своё ежегодное развитие позже естественной (ввиду проведения посевных работ), интенсивнее развивается, достигая более высоких значений РУГ, а затем теряет хлорофилл и исчезает (в результате уборки урожая) раньше естественной растительности.
В таблице 1 приведены разработанные признаки распознавания пашни и даётся характеристика особенностей их использования. При этом признаки основной группы (первые три в таблице 1 ) обеспечивают хорошую локальную разделимость практически в любом регионе России и, следовательно, применимы для всей территории страны. Признаки распознавания вспомогательной группы (последние три в таблице 1) обладают хорошим уровнем разделимости лишь в отдельных регионах для земель с определёнными типами севооборота. Это ограничивает область их применения для картографирования пахотных земель, но позволяет использовать при решении вспомогательных задач, таких как формирование обучающей выборки или фильтрация ошибочно распознанных участков на этапе постклассификационной обработки результатов.
Таблица 2 даёт представление о разделимости участков пашни и других типов земель на основе признаков распознавания основной группы.
естественная
растительность пашня
01.01.2002 01.01.2003 01.01.2004 01.01.2005 01.01.2006 01.01.2007 Дата
Рис. 6. Примеры многолетних временных рядов PVI с характерными для пашни и естественной растительности межгодовыми и сезонными различиями. L1/2 иллюстрирует один из признаков распознавания ИПЗ (см. таблицу 1)
Таблица 1. Признаки распознавания используемых пахотных земель на основе многолетних рядов данных MODIS
Название Формула Описание Особенности
D K R C V
Индекс кратчайшего сезона вегетации Li/2 = т'пУ1 - 4), J=1..N PVI PVI(tL) = PVI(tF) = max , t, > t , t„ < t L max' F max Минимальный для ряда лет отрезок времени, когда в течение года значения РУ1 превышали половину его сезонного максимума. + - - + -
Индекс весеннего развития растительности MSI = min У PVI„ J=1..N ^ J iospw Многолетний минимум интеграла РУ1 в период 1 января - 15 июня каждого года наблюдений. - + + + -
Индекс сезонного снижения фитомассы N у ' PVIm'mosw NSMI = const J=N У У PVI. j=1 iesw Нормированная сумма многолетних сезонных минимумов РУ1 в период 15 мая - 15 сентября каждого года наблюдений. + + + - +
Индекс межгодовых различий динамики растительно сти K = = min (Cor(PVI(Year), PVI(Yeart))) i, Je1..N J Минимум всех возможных значений межгодовых корреляций временных рядов РУ1 + + + - -
Индекс межгодовой изменчивости фитомассы D = SD (У PVL) Year, Стандартное отклонение сумм накопленных за различные годы значений РУ1. - + + - -
Разностный индекс сезонного пика вегетации T = Med(PVIYearJ - PVIYear') . , , r ^ max mean ' J=1..N Многолетняя медиана разности максимального и среднего значения РУ1 + + + + -
Примечания:
1) В таблице представлены следующие обозначения особенностей признаков:
D - устойчивость к пропускам в данных наблюдений: (+) - устойчив и (-) - неустойчив. K - квантованность значений признака: (+) - слабая и (-) - сильная. R - динамический диапазон: (+) - широкий и (-) - узкий. C - зашумлённость: (+) - низкая и (-) - высокая.
V - вероятность перепутывания с участками изменений естественной растительности: (+) - низкая и (-) - высокая.
2) В формулах приняты следующие обозначения математических операций: Cor - корреляция; SD - стандартное отклонение; Med - медиана.
В основу индекса кратчайшего сезона вегетации Li/2 в качестве признака распознавания пахотных земель положено наличие, в силу проведения агротехнических мероприятий (обработка почвы перед севом и уборка урожая), тенденций к сокращению продолжительности вегетационного периода сельскохозяйственных культур по сравнению с естественной растительностью. Значение Li/2 рассчитывается как минимальный для ряда лет период, в котором значения PVI превышают порог, определяемый половиной максимального значения ВИ в текущем году (рис. 6).
Предпосылками для разработки и использования индекса весеннего развития растительности MSI являются характерные для естественной и сельскохозяйственной растительности различия времени начала роста. В отличие от естественной растительности, время появления всходов на пахотных землях непосредственно связано со сроками сева, определяемыми, в свою очередь, состоянием почвы и её готовностью к проведению посевных работ. Исклю-
чение составляют весенние всходы многолетних и озимых культур, которые уже имеют развитые с осени предыдущего года вегетативные органы. Таким образом, значения индекса MSI, определяемого как минимальная многолетняя сумма значений PVI в границах временного окна, охватывающего весенний период года, имеют, как правило, существенно большие значения на участках естественной растительности, чем на пахотных землях.
Временной интервал спутниковых наблюдений, используемых для вычисления индекса сезонного снижения фитомассы NSMI, соответствует периоду вегетации культурной растительности, но также включает в себя моменты ожидаемого минимального покрытия пашни растительностью и, следовательно, низких значений PVT. Это моменты осенней обработки почвы, скашивания многолетних трав или предпосевных работ. При этом предполагается, что в период вегетации естественная растительность накапливает существенный объём фитомассы.
Таблица 2. Оценка разделимости пашни и других типов земель на основе признаков распознавания, получаемых по многолетним данным MODIS
Название признака
Индекс кратчайшего сезона вегетации
Индекс весеннего развития растительности
Индекс сезонного снижения фитомассы
Фрагмент изображения
Показатели разделимости классов
Трансформированная дивергенция (0^2000)
Примечание: на двух первых (сверху) фрагментах изображения пахотные земли имеют тёмный оттенок, естественная растительность - светлый; на нижнем фрагменте - наоборот. В столбце «гистограмма» везде тёмный цвет соответствует классу «пашня», светлый - классу «естественная растительность».
В качестве индекса межгодовых различий динамики растительности К принимается значение минимального коэффициента линейной корреляции между временными рядами значений РУГ для всех возможных сочетаний пар лет заданного периода наблюдений. В основу разработки индекса К положено предположение о снижении его значений на пахотных землях (по сравнению с участками естественной растительности) за счёт наличия севооборота и связанных с ним больших межгодовых различий сезонной динамки РУГ.
Индекс межгодовой изменчивости фитомассы D, характеризуемый стандартным отклонением многолетних значений накопленной за год суммы РУГ, обеспечивает наиболее высокий уровень детектируе-мости пахотных земель при наличии в севообороте чистых паров. Вместе с тем, признак отличается значительными уровнями пространственной вариабельности и случайного шума. Кроме того, при использовании признака необходимо учитывать близость его значений для участков пахотных земель и некоторых
других типов объектов, таких как вырубки, гари и карьеры, с тем, чтобы минимизировать влияние этого фактора на точность распознавания.
Использование разностного индекса сезонного пика вегетации Т позволяет детектировать объекты, внутригодовой ход РУГ которых демонстрирует значительные отклонения максимальных значений летнего периода от среднегодовых значений. Очевидно, что в силу использования агротехники и ускоренного развития растительности на пахотных землях, для них характерны более существенные различия между максимальными и среднегодовыми значениями вегетационного индекса. Вычисление медианы указанных выше многолетних разностей РУГ позволяет исключить экстремальные значения, сопряжённые с более высоким уровнем случайного шума.
Используя критерий Колмогорова-Смирнова, установлено, что распределения значений признаков каждого из классов подчиняются закону нормального распределения с уровнем значимости а< 0,001. При использовании для распознавания метода мак-
симального правдоподобия можно оценить вероятность правильного детектирования класса как долю отнесённых к данному классу пикселов по отношению к числу пикселов класса в опорной выборке. На рис. 7 представлено пространственное распределение вероятности правильного детектирования пахотных земель на примере признака кратчайшего сезона вегетации.
Рис. 7. Пространственное распределение вероятности верного детектирования ИПЗ (%), основанное на оценке локальной разделимости по признаку L1/2
Из рисунка можно видеть, что значение вероятности демонстрирует заметную географическую вариабельность. При этом максимальная вероятность правильного детектирования относится к территориям, где заданные классы представлены примерно в равных пропорциях, а естественная растительность сильно отличается от пашни по средним значениям признака (например, восток европейской части России, юг Сибири, Дальний Восток). В то же время регионы с наибольшими или с наименьшими площадями пашни обладают более низкими значениями вероятности правильного распознавания в связи с сильным доминированием одного из классов. При этом вероятность верного распознавания не опускается ниже 80%, а для многих участков достигает 100%, что демонстрирует возможность использования признака кратчайшего сезона вегетации для выявления пахотных земель.
Для классификации был использован локально-адаптивный алгоритм и программный комплекс LAGMA. При этом ввиду особенностей разработанных признаков распознавания была выбрана схема независимой классификации с раздельным использованием каждого из признаков. В частности, индекс весеннего развития растительности не позволяет достоверно выявлять объекты класса ИПЗ, занятые на протяжении ряда лет многолетними травами, в то время как индекс сезонного снижения фитомас-сы позволяет это сделать. Для учёта этих особенностей результаты независимой классификации по каждому из этих признаков объединялись. Для исключения из результатов объединения не относящихся к сельскохозяйственным землям участков, что выражается в более высоком значении индекса кратчайшего сезона вегетации, применяется логическое умножение результата независимого выявления по
признаку L1/2 и результатов объединения признаков MSI и NSMI. Схему интеграции результатов независимой классификации по отдельным признакам можно представить в виде логического выражения AL = R(L1/2) П (R(MSI) U R(NSMI)), где AL - маска пахотных земель, R(X) - результат независимой классификации по признаку X, а L1/2, MSI и NSMI - соответствующие признаки из таблицы 1, П и U - операции логического умножения и сложения соответственно.
Результаты распознавания используемых пахотных земель
Поскольку признаки распознавания основаны на многолетних сериях данных, получаемые на их основе результаты характеризуют состояние пашни в соответствующие используемым спутниковым данным интервалы времени. При этом разработанный метод позволяет проводить ежегодное обновление карты пахотных земель, что даёт возможность выявлять изменения в режимах их использования. Это позволяет обнаруживать участки как введённой в сельскохозяйственный оборот, так и заброшенной пашни.
Результаты выявления используемых пахотных земель России на основе описанных спектрально-динамических признаков распознавания, соответствующих временному интервалу 2004-2009 гг., представлены на рис. 8.
На основе полученной маски пахотных земель были рассчитаны площади пашни в границах субъектов РФ. Сравнение результатов выявления используемых пахотных земель с данными Министерства сельского хозяйства РФ (www.gks.ru) представлено на рис.9.
Для валидации полученной карты и оценки уровня ошибок распознавания ИПЗ были использованы векторные карты границ полей, полученные по спутниковым данным высокого пространственного разрешения (5-10 м) на территорию Южного федерального округа. Фрагмент распознанных по данным MODIS участков ИПЗ в сравнении с векторными границами полей представлен на рис. 10.
Для количественной оценки точности распознавания ИПЗ на основе локально-адаптивной классификации и разработанных признаков был использован метод на основе оптимума Парето [9]. Метод определяет детальность тематической карты с помощью величин ошибок первого и второго рода, а также обеспечивает наглядное сопоставление результата распознавания с эталоном.
Точность распознавания ИПЗ всегда ограничена пространственным разрешением спутниковых данных. Идеальный продукт распознавания - лишённый погрешности, связанной с качеством спутниковых и обучающих данных, а также эффективностью метода классификации - будет обладать ошибками первого и второго рода, обусловленными пространственным разрешением.
Рис. 8. Карта пахотных земель России, полученная по данным МОББ анным Росстата, тыс.га
у=0,% я2=о,<, 1х+208
о
с^^
Яэг °
8эо
5000 4000 3000 2000 1000
О 1000 2000 3000 4000 5000 6000 По данным МОВ1Б, тыс.га
Рис 9. Сравнение площади ИПЗ по данным МОВ1Б с данными Минсельхоза РФ по посевным площадям.
Данные получены по 43 субъектам РФ, охватывающим более 97% площади пахотных земель России
При этом, согласно определению оптимума Па-рето, модификация такого продукта с целью снижения пропуска цели неизбежно приводит к повышению ложной тревоги и наоборот. В евклидовом пространстве величин ошибок первого и второго рода множество идеальных (Парето-оптимальных) продуктов описывается кривой Парето. Например, в работе [10] описана методика вычисления оптимума Парето на основе эталонных данных для заданной величины пространственного разрешения используемых спутниковых данных.
Рис. 10. Фрагмент маски ИПЗ (серый) с векторными границами полей (чёрный), используемыми для оценки точности классификации
Фактическая точность полученного продукта распознавания ИПЗ, определённая в виде величин ошибок первого и второго рода, может быть представлена в виде точки в указанном пространстве. Мерой качества распознавания является удалённость точки от линии оптимума, а взаимное расположение различных точек позволяет сравнивать качество различных продуктов.
В настоящей работе для оценки результатов распознавания пахотных земель при участии экспертов
был создан эталон для тестового участка на основе векторных границ сельскохозяйственных полей, оцифрованных по спутниковым данным высокого разрешения. На рис.11 представлены фрагменты, использованные для валидации продуктов на основе оптимума Парето. В настоящей работе расстояние до кривой измерялось по точке равенства ошибок первого и второго рода. Полученное значение евклидова расстояния 0,07 между точкой и оптимумом (рис.12) указывает на высокую точность распознавания пахотных земель на примере представленного фрагмента.
Рис. 11. Фрагменты полученной карты ИПЗ (а) и эталонной карты (б) для оценки точности классификации на основе оптимума Парето
Заключение Результаты экспериментальных исследований позволяют оценить эффективность разработанного метода распознавания пахотных земель на основе данных спектрорадиометра MODIS. Предложенный в работе оригинальный методический подход, предполагающий использование многолетних рядов спутниковых данных для разработки признаков распознавания, позволил обеспечить высокий уровень отделимости участков пашни от других категорий сельскохозяйственных земель и естественной растительности.
Разработанные признаки распознавания основаны на учёте сезонных и межгодовых вариаций спек-
трально-отражательных характеристик пахотных земель, отличительные черты которых формируются под влиянием проводимых на них агротехнических мероприятий и наличия севооборота.
Ложная тревога
(0,08; 0, Щ
/0,07
0 0,05 0,10 0,15 0,20
Пропуск цели
Рис. 12. Кривая оптимума Парето в пространстве ошибок первого и второго рода распознавания ИПЗ. Точки соответствуют минимально возможным при использовании данных MODIS значениям ошибок распознавания первого и второго рода (точка на кривой
выбрана при условии равенства ошибок обоих видов) и фактическим ошибкам распознавания ИПЗ (0,08;0,11)
При картографировании и мониторинге больших территорий использование локально-адаптивных алгоритмов классификации, в том числе для выявления пахотных земель, обосновано необходимостью учёта значительных пространственных вариаций используемых признаков распознавания. Использованный в работе алгоритм локально-адаптивной обучаемой классификации LAGMA продемонстрировал возможность надёжного выделения участков пахотных земель и построения соответствующей карты для всей территории России.
Высокий уровень согласованности полученных оценок площадей пахотных земель на уровне субъектов РФ с данными официальной статистики Мин-сельхоза РФ служит подтверждением эффективности метода. Сравнение ошибок распознания пахотных земель, оценённых на основе детальных опорных данных, и кривой оптимума Парето свидетельствует о том, что разработанный метод обеспечивает близкий к предельно достижимому уровню точности при заданной величине пространственного разрешения MODIS.
Литература
1. Cihlar, J. Classification by progressive generalization: a new automated methodology for remote-sensing multichannel data. / J. Cihlar, Q. Xiao, J. Beaubien, K. Fung, R.Latifovic // International Journal of Remote Sensing. -1998. - V.19. - P.2685-2704.
2. Дейвис, Ш.М. Дистанционное зондирование: количественный подход / Ш.М. Дейвис, Д.А. Ландгребе, Т.Л. Филипс; пер. с англ. - М.: Недра, 1983. - 415 с.
3. Уваров, И.А. Алгоритм и программный комплекс распознавания типов земного покрова на основе локаль-
но-адаптивной обучаемой классификации спутниковых изображений / И.А. Уваров, С.А. Барталёв // Современные проблемы дистанционного зондирования Земли из космоса: Физические основы, методы и технологии мониторинга окружающей среды, потенциально опасных явлений и объектов: сб. научных статей. - М: ООО «ДоМира», 2010. - Т.7, № 1. - С.353-365.
4. Лупян, Е.А. Технологии спутникового мониторинга в сельском хозяйстве России / Е.А. Лупян, С.А. Барталёв , И.Ю. Савин // Аэрокосмический курьер. - 2009.
- № 6. - С.47-49.
5. Барталёв, С.А. Классификация некоторых типов сельскохозяйственных посевов в южных регионах России по спутниковым данным MODIS / С.А. Барталёв, Е.А. Лупян, И.А. Нейштадт, И.Ю. Савин // Исследование Земли из космоса. - 2006. - №3. - С.68-75.
6. Барталёв, С.А. Метод выявления используемых пахотных земель по данным дистанционного зондирования со спутников / С.А. Барталёв, Е.А. Лупян, И.А. Нейштадт // Современные проблемы дистанционного зондирования Земли из космоса. Физические основы, методы и технологии мониторинга окружающей среды, потенциально опасных явлений и объектов: сб. научных статей. - М.: ООО «Азбука-2000». -2006. - В. 3. - Т.2. - С.271-280.
7. Плотников, Д.Е. Признаки распознавания пахотных земель на основе многолетних рядов данных спутникового спектрорадиометра MODIS / Д.Е. Плотников, С.А. Барталёв , Е.А. Лупян // Современные проблемы дистанционного зондирования Земли из космоса: Физические основы, методы и технологии мониторинга окружающей среды, потенциально опасных явлений и объектов: сб. науч. статей. - М: ООО «ДоМира», 2010.
- Т.7. № 1. - С.330-341.
8. Нейштадт, И.А. Построение безоблачных композитных спутниковых изображений MODIS для мониторинга растительности // Современные проблемы дистанционного зондирования Земли из космоса: сб. науч. статей. - М.: «Азбука-2000», 2006. - Т. 2. - С. 359-365.
9. Keeney, R.L. Decisions with multiple objectives: Preferences and value tradeoffs. - New York: Wiley, 1976.
10. Boschetti, L. Analysis of the conflict between omission and commission in low spatial resolution dichotomic thematic products: The Pareto Boundary / L. Boschetti, P.F. Stéphane, A.B. Pietro // Remote Sensing of Environment. - 2004 - Vol. 91. - P. 280-292.
References
1. Cihlar, J. Classification by progressive generalization: a new automated methodology for remote-sensing multichannel data. / Cihlar J., Xiao Q., Beaubien J., Fung K., Latifovic R. // International Journal of Remote Sensing. -1998. - V.19. - P.2685-2704.
2. Davis, S.M. Remote sensing: the quantitative approach / Davis S.M., Landgrebe D.A., Phillips T.L. // Moscow, Nedra, - 1983. - 415 p.
3. Uvarov, I.A. Algorithm and a software package for vegetation types locally-adaptive supervised classification / Uvarov I.A., Bartalev S.A. // Contemporary problems of Earth remote sensing: Physical basics, methods and technologies of environmental and hazardous phenomena monitoring. Scientific papers compilation. M. : - "Domira", - 2010. - V.7, - №1. - P.353-365.
4. Loupian, E.A. The remote sensing technologies for agricultural monitoring of Russia / Loupian E.A., Bartalev S.A., Savin I.U. // Aerospace Courier. - 2009. - №6. -P.47-49.
5. Bartalev, S.A. Several crops detection using MODIS in the South of Russia / Bartalev S.A., Loupian E.A., Neishtadt I.A., Savin I.U. // Issledovanie Zemli iz kosmosa. - 2006. - №3. - P.68-75.
6. Bartalev, S.A. Arable lands detection method based on remote sensing data. / Bartalev S.A., Loupian E.A., Neishtadt I.A. // Contemporary problems of Earth remote sensing: Physical basics, methods and technologies of environmental and hazardous phenomena monitoring. Scientific papers compilation. M.: - «Azbuka-2000», - 2006. -Issue 3. - V.2. - P.271-280.
7. Plotnikov, D.E. The recognition features to map arable lands based on multi-annual MODIS Earth observation data. / Plotnikov D.E., Bartalev S.A., Loupian E.A. // Contemporary problems of Earth remote sensing: Physical basics, methods and technologies of environmental and hazardous phenomena monitoring. Scientific papers compilation. M.: - "Domira", - 2010. - V.7, - №1. - P.330-341.
8. Neishtadt, I.A. Method of cloudless MODIS images creation for vegetation monitoring // Contemporary problems of Earth remote sensing. Scientific papers compilation. M.: - «Azbuka-2000», - 2006. - V.2. - P.359-365.
9. Keeney, R.L. Decisions with multiple objectives: Preferences and value tradeoffs / New York: - Wiley, - 1976.
10. Boschetti, L. Analysis of the conflict between omission and commission in low spatial resolution dichotomic thematic products: The Pareto Boundary / Boschetti L., Stéphane P.F., Pietro A.B. // Remote Sensing of Environment, - 2004 - Vol. 91. - P. 280-292.
RECOGNITION OF ARABLE LANDS USING MULTI-ANNUAL SATELLITE DATA FROM SPECTRORADIOMETER MODIS AND LOCALLY ADAPTIVE SUPERVISED CLASSIFICATION
S.A. Bartalev, V.A. Egorov, E.A. Loupian, D.E. Plotnikov, I.A. Uvarov Space Research Institute of the Russian Academy of Sciences
Abstract
The arable lands recognition method is developed based on multi-annual time-series of remote sensing data acquired by spectroradiometer MODIS on board of Terra and Aqua satellites. The method involves producing of the recognition features set, which exploits differences of seasonal and inter-annual changes of spectral reflectance for arable lands on one hand and other types of agricultural lands and natural vegetation on another hand. The arable lands recognition utilizes the locally-adaptive supervised classification algorithm, which accounts the spatial variability of the considered features for classes to be discriminated.
The developed method has been applied to produce the arable lands map for entire Russia. The arable lands map validation based on Pareto optimum approach and reference data has been performed for the test region in order to estimate the method's accuracy and potential for its further enhancement.
Key words: remote sensing, satellite spectroradiometer, supervised classification, recognition features, agricultural lands monitoring.
Сведения об авторах Барталёв Сергей Александрович, 1961 года рождения. В 1984 году окончил факультет прикладной космонавтики Московского института геодезии, аэрофотосъёмки и картографии по специальности «Исследование природных ресурсов аэрокосмическими методами». Доктор технических наук (2007) по специальности «Приборы и методы экспериментальной физики». Заведующий Лабораторией спутникового мониторинга наземных экосистем Института космических исследований Российской академии наук. Область научных интересов включает дистанционное зондирование, методы обработки изображений, технологии спутникового мониторинга экосистем. Автор более 200 научных публикаций. E-mail: b artalev@smis. iki. rssi. ru
Sergey Alexandrovich Bartalev (b. 1961) graduated (1984) from Moscow Institute of Geodesy, Aerophotosurvey and Cartography, Space Applications Faculty, majoring in Remote Sensing of Natural Resources. He has Doctor of Science diploma (2007) in experimental physics. He is working as head of Terrestrial Ecosystems Monitoring Laboratory at the Space Research Institute of Russian Academy of Sciences. His research interests include remote sensing, image processing methods, methods and technologies for monitoring of terrestrial ecosystems. He is author more then 200 scientific publications.
Егоров Вячеслав Александрович, 1979 года рождения. В 2001 году окончил Московский физико-технический институт (Государственный университет) в степени магистра по специальности «Прикладные математика и физика». Кандидат технических наук (2006) по специальности «Аэрокосмические исследования Земли, фотограмметрия». Работает в ИКИ РАН в отделе технологий спутникового мониторинга в должности научного сотрудника. Область научных интересов включает технологии дистанционного зондирования Земли, программирование, ГИС, методы оценки изменений в лесах, основанные на многолетних сериях данных спутниковых наблюдений. В списке публикаций В.А. Егорова 30 печатных работ. E-mail: egorov@d902. iki. rssi. ru
Viacheslav Alexandrovich Egorov (b. 1979) graduated (2001) from Moscow Institute of Physics and Technology (State University) (MIPT), majoring in Applied Mathematics and Physics. He has Philosophy doctor of Science diploma (2006) in Earth Aerospace research and photogrammetry. He is working as a researcher in the Terrestrial Ecosystems Satellite Monitoring laboratory (Space Monitoring Information Support department) of the RAS Space Research Institute. His current research interests include technologies for remote sensing, programming, GIS, methods of assessing changes in forests based on long-term series of satellite data. V.A. Egorov is co-author of 30 papers and scientific publications.
Лупян Евгений Аркадьевич, 1961 года рождения. Окончил Московский физико-технических институт в 1984 году. Доктор технических наук. В настоящее время заместитель директора ИКИ РАН, отвечает за направление «Исследование Земли из Космоса». Руководит отделом «Технологий спутникового мониторинга». Преподаёт в МФТИ и МГТУ им. Баумана. Более 25 лет работает в области создания методов систем и технологий, связанных с использованием данных дистанционного зондирования. Является участником многих проектов, посвящённых созданию и внедрению систем дистанционного мониторинга. Основной областью научных интересов в последние годы является создание методов, систем и технологий обработки данных дистанционных наблюдений, создание распределённых автоматизированных информационных систем сбора, обработки и распространения спутниковой информации. Автор более 200 научных работ.
E-mail: evgeny@d902. iki. rssi. ru
Evgeny Arkadievich Loupian (b.1961) graduated (1984) from Moscow Institute of Physics and Technology (State University). He has Doctor of Science diploma (1998) and works as Deputy Director of Space Research Institute of Russian Academy of Science, supervising Earth Remote Sensing research. He works at the development of methods and techniques of automated satellite data processing. He participates in many projects concerning satellite data accumulation, processing and distribution. He is author of more then 200 scientific publications.
Плотников Дмитрий Евгеньевич, 1984 года рождения. В 2007 года окончил Московский физико-технический институт (Государственный университет) в степени магистра по специальности «Прикладные математика и физика». Закончил аспирантуру Института космических исследований РАН по специальности «Приборы и методы экспериментальной физики». Работает в ИКИ РАН в отделе технологий спутникового мониторинга в должности младшего научного сотрудника. Область научных интересов включает дистанционное зондирование, программирование, ГИС, создание адаптивных алгоритмов классификации, использование временных серий, распознавание типов растительного покрова. В списке публикаций Д.Е. Плотникова 13 печатных работ. E-mail: [email protected]
Dmitry Evgenyevich Plotnikov (b.1984) graduated (2007) from Moscow Institute of Physics and Tecnology (State University) (MIPT), majoring in Applied Mathematics and Physics. He is working as a junior staff scientist in the Terrestrial Ecosystems Satellite Monitoring laboratory (Space Monitoring Information Support department) of the RAS Space Research Institute. His research interests include Earth remote sensing, GIS, programming, adaptive classification algorithms development, time series and vegetation recognition. He is currently co-author of 13 papers and scientific publications.
Уваров Иван Александрович, 1981 года рождения. В 2004 году окончил Московский государственный университет им. М.В. Ломоносова по специальности «География и картография». Работает младшим научным сотрудником в Лаборатории спутникового мониторинга наземных экосистем Института космических исследований Российской академии наук. Область научных интересов: дистанционное зондирование, ГИС, разработка программного обеспечения, методы обучаемой классификации, распознавание типов земного покрова. В списке публикаций И.А. Уварова 28 печатных работ. E-mail: uvarov@smis. iki. rssi. ru
Ivan Aleksandrovich Uvarov (b. 1981) graduated (2004) from Lomonosov Moscow State University majoring in geography and cartography. He is working as a junior staff scientist in the Terrestrial Ecosystems Satellite Monitoring laboratory of the RAS Space Research Institute. His area of interest is remote sensing, GIS, software development, methodology of supervised classification, land cover recognition. I.A. Uvarov is a co-author of 28 scientific publications.
Поступила в редакцию 11 октября 2010 г.
