Кластерный анализ и классификация с обучением многоспектральных данных дистанционного зондирования Земли Текст научной статьи по специальности «Математика»

УДК 004.4:528.9

Кластерный анализ и классификация с обучением многоспектральных данных дистанционного зондирования Земли

В.В. Асмус", А.А. Бучнев6, В.П. Пяткин6*

a Научно-исследовательский центр "Планета", РОСКОМГИДРОМЕТ 123242 Россия, Москва, Большой Предтеченский пер., 7 b Институт вычислительной математики и математической геофизики СО РАН

630090 Россия, Новосибирск, пр. Лаврентьева, 61

Получена 15.11.2008, окончательный вариант 10.12.2008, принята к печати 10.03.2009 Рассматриваются вопросы, связанные с проблемой выбора адекватных алгоритмов распознавания многоспектральных данных дистанционного зондирования Земли. Представлена система контролируемой классификации, основанная на стратегии максимального правдоподобия для нормально распределенных векторов признаков. Описывается система кластерного анализа, включающая алгоритм К-средних и метод анализа мод многомерной гистограммы.

Ключевые слова: дистанционное зондирование Земли, распознавание данных, контролируемая классификация, неконтролируемая классификация, кластерный анализ, решающее правило, обучение классификатора, метод K-средних, многомерная гистограмма.

1. Введение

Эффективность дистанционных исследований Земли из космоса во многом определяется используемыми методами тематической обработки данных дистанционного зондирования Земли (ДДЗЗ). При этом центральной при тематической обработке, безусловно, является система классификации. В модуль распознавания программного комплекса обработки ДДЗЗ, разработанного в ИВМиМГ СО РАН совместно с НИЦ "Планета" Роскомгидромета РФ, включены алгоритмы контролируемой и неконтролируемой классификации многоспектральных ДДЗЗ [1,2]. Центральный вопрос интерпретации ДДЗЗ (вопрос повышения качества дешифрирования) непосредственно связан с проблемой выбора адекватного алгоритма распознавания. Возникающие при этом трудности обусловлены следующими причинами:

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

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

* Corresponding author E-mail address: [email protected]

1 © Siberian Federal University. All rights reserved

данных и данных, предъявляемых на распознавание ("загрязнение" выборок смешанными векторами измерений, т.е. векторами, которые образуются при попадании в элемент разрешения съемочной системы нескольких природных объектов, неточное соответствие обучающих данных, получаемых с помощью кластеризации, истинным тематическим классам, помехи аппаратуры, влияние атмосферных условий и т.п.) [2].

Таким образом, современный опыт автоматизированного распознавания ДДЗЗ показывает: заранее практически невозможно установить, какой алгоритм будет лучше с точки зрения соотношения точности классификации и стоимости. Поэтому в распознающую систему целесообразно закладывать несколько алгоритмов и выбор оптимального алгоритма проводить эмпирически на этапе обучения по результатам классификации тестовых данных. Выбранный алгоритм используется затем для распознавания всего набора векторов измерений.

2. Контролируемая классификация

Разработанная в ИВМиМГ СО РАН совместно с НИЦ "Планета" система контролируемой классификации (классификации с обучением) ДДЗЗ в программном комплексе состоит из семи классификаторов (один поэлементный классификатор и шесть объектных), основанных на использовании байесовской стратегии максимального правдоподобия, и двух объектных классификаторов, основанных на минимуме расстояния. Эта система является частью программного комплекса по обработке ДДЗЗ [3].

2.1. Поэлементная классификация. Под элементом здесь понимается Ж-мерный вектор измерений (признаков) х = (х\,... )Т, где N — число спектральных диапазонов. Предполагается, что векторы х имеют в классе ш^ нормальное распределение N(ш^, В¿) со средним и ковариационной матрицей В^. В этом случае байесовская стратегия максимального правдоподобия для поэлементного классификатора формулируется следующим образом [4-6].

Пусть И = (ш1,... ,шт) — конечное множество классов, р(ш^) — априорная вероятность класса ш^. Тогда дискриминантная функция класса ш^ имеет вид

Классическое решающее правило для классификатора принимает следующий вид: вектор х заносится в класс ш^, если дг(х) > д^(х) для всех ] = г.

Для класса ш^ следующим образом определим параметр Т:

д^(х) = 1п(р(ш<)) - 0.51п(|В<|) - 0.5(х - ш^)ТВ-1(х -

(1)

Т = 1пр(ш1) - 0.5Л(Ж^) - 0.51п ВЛ,

(2)

где Л(Ж^) — критическое значение уровня Q распределения х2. Пусть ^ — переменная, значение которой зависит от параметра классификатора ЬНт:

2

и

г=1

т ГТ1

тах ТI,

тт Т,

Т

-Ж,

г=1

т

если Шг = 1, если Шг = 2,

если Шг = 3,

если Шг = 4,

(3)

Т /ш, если ЬНг = 5.

Тогда решающее правило для классификатора с учетом (3) принимает следующий вид [2]: вектор х заносится в класс шг, если дг(х) > д^ (х) для всех ] = г и дг(х) > Ьг (при Ьг = = Тг это означает следующее ограничение расстояния Махалонобиса до центра класса: (х—тг)ТБ-1(х—тг) < Л(М, Q)). В противном случае вектор заносится в класс отклоненных векторов (класс с номером т +1).

Значения Л(Ы, д) для размерности вектора N < 30 находятся из статистических таблиц. Например, для N = 4 и уровня Q = 0.05 (т. е. 5% векторов может быть отклонено) Л = 9.488. При N > 30 для нахождения значений Л(М^) используется аппроксимация

¿-¿ю«N ^ - ^+, (4)

где — значение стандартизованной нормальной величины для вероятности 1 — Q (в

частности, для уровня Q = 0.05 значение мо.95 « 1.645).

2.2. Объектная классификация. Под объектом мы понимаем блок смежных векторов квадратной или крестообразной формы. Поскольку физические размеры реально сканируемых пространственных объектов, как правило, больше разрешения съемочных систем, между векторами данных существуют взаимосвязи. Использование информации подобного рода дает возможность повысить точность классификации, если пытаться распознавать одновременно группу смежных векторов - объект в приведенном выше смысле. Рассмотрим вектор (объект) X(х\,..., хь)Т, состоящий из смежных N—мерных векторов хг, г = 1, . . . , Ь (например, в окрестности 3*3, 5*5,... элементов; мы работаем с объектами двух видов -квадратными либо крестообразными). Решение об отнесении центрального элемента объекта тому или иному классу принимается на основе результата классификации всего объекта.

Такой подход порождает целое семейство решающих правил. Во-первых, это использование принципа голосования, т.е. независимая классификация элементов объекта и отнесение центрального элемента к тому классу, которому было отнесено большинство элементов объекта. Во-вторых, это применение текстурных операторов (простейший пример - описание объекта через вектор средних составляющих его элементов) с последующим отнесением центрального элемента классу, к которому был отнесен параметр, характеризующий X. В-третьих, описание объекта случайным марковским полем, т. е. р(Х\шг) = р(х\\х2,... ,хь; ).. .р(хь\^г). В этом случае модель выглядит следующим образом. Пусть вектор х имеет в классе нормальное распределение N(тг, Бг) со средним тг и ковариационной матрицей Бг. Тогда вектор также нормально распределен в классе со средним Мг размерности МЬ и ковариационной матрицей Кг размерности NЬx NЬ. Оценка этой матрицы при больших значениях МЬ (требуется очень большое количество обучающих данных), а также ее обращение на практике трудно реализуемо. Поэтому введем упрощающие структурные предположения. Если считать, что корреляция между элементами объекта во всех зонах съемки одинакова, то ковариационную матрицу Кг можно представить в

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

ь

матрицу Бг. Если Кг является единичной, то р(Х\шг) = П р(хг\шг) и мы имеем известное

г=1

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

Приведем алгоритмы работы некоторых объектных классификаторов. Предположим, что векторы внутри блока независимы. Будем рассматривать векторы, составляющие объект X, как один вектор размерности МЬ. Тогда дискриминантная функция класса шг имеет вид

ь

дг(Х) = 1п(р(шг)) - 0.5Ь 1п( | В г |) - 0.5^(хг - ш^)Т В-1 (х1 - ш>).

1=1

Решающее правило для данных классификаторов принимает следующий вид: центральный элемент объекта X заносится в класс шг, если дг(Х) > д^ (X) для всех ] = г и дг(Х) > Ьг. В противном случае центральный элемент объекта заносится в класс отклоненных векторов. Здесь определяется по (2) и (3) с Л = Л(ЬЖ, Q).

Снова считаем, что векторы внутри блока независимы. Классифицируется вектор х,

равный среднему по всем векторам объекта X:

1 Ь

х = 1^хг. (5)

г=1

Дискриминантная функция класса шг имеет вид

д.¿(X) = 1п(р(шг)) - 0.511п(1Вг1) - 0.5Ь(х - ш^)ТВ-1(х - шг).

Решающее правило для данных классификаторов таково: центральный элемент объекта X заносится в класс шг, если дг(X) > д^(X) для всех ] = г и дг(X) > ¿г. В противном случае центральный элемент объекта заносится в класс отклоненных векторов. Здесь Ьг определяется по (2) и (3).

Классифицируется средний вектор (5) блока в предположении, что векторы внутри блока независимы и ковариационные матрицы равны единичной. Фактически это объектные классификаторы, решающие правила которых основаны на минимуме евклидова расстояния до центра класса. Дискриминантная функция класса шг имеет вид

дг(X) = 1п(р(шг)) - 0.5(х - шг)ТВ-1(х - шг).

Решающее правило для данных классификаторов: центральный элемент объекта X заносится в класс шг, если дг^) > д^(X) для всех ] = г и дг^) > ¿г. В противном случае центральный элемент объекта заносится в класс отклоненных векторов. Здесь Ьг = -Л, где Л > 0 - число, задаваемое пользователем.

Система классификации содержит также объектные классификаторы, основанные на модели каузального марковского случайного поля первого и третьего порядка.

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

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

среднем пессимистическими. Анализируя эти данные, можно оценить (проконтролировать) качество обучения.

Как отмечалось выше, возможны ситуации, при которых нарушается условие о случайности поля измерений, следствием чего выступают нулевые дисперсии в некоторых каналах. Тогда формулы типа (1) становятся неприменимыми, т. к. ковариационные матрицы B вырожденные. Для исправления подобных ситуаций в системе классификации предусмотрена функция добавления к спектральным каналам с нулевой дисперсией гауссовского шума с нулевым средним и единичной дисперсией.

Результатом работы классификаторов в рабочем режиме служит одноканальное (байтовое) изображение, значениями пикселов которого являются номера классов. Это изображение окрашивается в предопределенные цвета, которые в интерактивном режиме могут быть заменены на цвета, определяемые пользователем. Кроме того, к этому изображению можно применить функцию редактирования, которая определяется как уточнение карты классификации на основе учета контекста без изменения перечня ранее выделенных классов. Эта функция может работать в двух режимах: в режиме Vote, при котором центральный пиксел окрестности 3*3 заменяется модой гистограммы окрестности, и в режиме Alísame, при котором центральный пиксел такой же окрестности меняется только тогда, когда все окружающие его пикселы имеют одинаковое значение.

Система контролируемой классификации имеет следующие характеристики: число обучающих образов - до 9, число классов - до 15, число обучающих и контрольных полей в классе - до 10, размер каждого поля - до 50*50 векторов, размер объекта - от 1*1 до 11*11, размерность векторов данных не ограничивается. На рис. 1 представлена тематическая карта ледовой обстановки в восточном секторе Арктики, полученная с использованием контролируемой классификации. В качестве распознаваемых изображений использовались мозаики радиолокационных и радиометрических снимков с ИСЗ "Океан-О1". Для выбора тестовых участков использовалось цветосинтезированное изображение.

3. Неконтролируемая классификация

Неконтролируемая классификация (кластерный анализ) в программном комплексе представлена двумя алгоритмами - методом K-средних и методом анализа мод многомерной гистограммы [3].

3.1. Метод K-средних. Этот подход основан на итеративной процедуре отнесения векторов признаков классам по критерию минимума расстояния от вектора до центра класса. Оптимальным считается такое разбиение входных векторов на кластеры, при котором внутриклассовый разброс не может быть уменьшен при переносе какого-либо вектора из одного кластера в другой.

Алгоритм состоит в выполнении следующих шагов:

1. На основе заданного соотношения а чистых и смешанных векторов производится разделение векторов на чистые и смешанные. Под смешанными мы понимаем векторы, компоненты которых либо формируются за счет попадания в поле зрения съемочной аппаратуры нескольких объектов, либо искажены влиянием фона. С этой целью вначале для исходного набора векторов измерений рассчитывается градиентное изображение и одновременно строится гистограмма градиентов. Исходя из заданного а, по гистограмме определяется порог, разделяющий векторы на смешанные и чистые.

2. Объединение чистых векторов в связные компоненты. На этом этапе все чистые век-

торы объединяются в связные компоненты, которые последовательно нумеруются. Соответствующий алгоритм, идейно близкий к алгоритму заполнения областей с произвольной границей по критерию связности [3], может выделять и нумеровать одновременно любое количество многосвязных областей без ограничений на их форму и ширину контуров. Для каждой связной компоненты вычисляется вектор средних.

Рис. 1

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

4. Распределение связных компонент по кластерам. На этом этапе связные компоненты получают новые номера. Новый номер присваивается компоненте в соответствии с номером кластера, в который попал вектор средних этой компоненты.

5. Кластеризация смешанных векторов. На завершающем этапе производится неитеративная кластеризация смешанных векторов по принципу минимума расстояния до центров кластеров С\,... . Смешанный вектор г будет отнесен к ближайшему кластеру если

\\Ci — z\\ < 0.5 max\\Ci — Ci jj, i,l = l,...,k, k = l.

Рис. 2 Рис. 3

Рис. 3 иллюстрирует результат работы кластеризации методом ^-средних изображения, представленного на рис. 2. Рис. 3 демонстрирует результат работы алгоритма для следующих входных данных: количество выделяемых кластеров 15, соотношение количества смешанных и чистых векторов а = 0.35.

3.2. Метод анализа мод многомерной гистограммы. В основе второго подхода лежит предположение, что исходные данные являются выборкой из многомодового закона распределения, причем векторы, отвечающие отдельной моде, образуют кластер. Таким образом, задача сводится к анализу мод многомерных гистограмм. В программный комплекс включена реализация метода, описание основных шагов которого приведено в [3].

Гистограмма генерируется последовательным просмотром векторов данных и сравнением каждого вектора с текущим списком векторов. При этом либо изменяется соответствующее значение частоты, либо вектор добавляется в список. Для вычисления адресов векторов в списке используется хэш-кодирование. Первым шагом модального анализа является поиск ближайших соседей данного вектора списка среди других векторов списка. По определению вектор x есть ближайший сосед вектора y, если \xi — yi\ < l для i = l,..., N. Каждый из возможных ближайших соседей данного вектора x может быть получен из него прибавлением вектора сдвига, компоненты которого принимают значения из множества { — 1, 0,1}. Алгоритмически i-й вектор сдвига, i = l, 2,..., 3N — l, можно получить, уменьшив на 1 каждый из коэффициентов представления числа i в троичной системе счисления. Поскольку в реальной гистограмме присутствуют далеко не все ближайшие соседи, то для эффективного их поиска векторы предварительно упорядочиваются в многомерные бинарные деревья. В этом случае время поиска всех ближайших соседей данного вектора становится пропорциональным числу реально существующих соседей. При построении дерева векторы x рассматриваются как N-мерные ключи. Вначале рассчитываются дисперсии по всем координатам векторов и определяется координата j, имеющая максимальную дисперсию. Медианное значение выборки по этой координате используется в качестве ключа для разделения множества векторов на два подмножества: в одно подмножество помещаются векторы, значение которых по координате j меньше порогового значения, а в другое — векторы, у которых значение координаты превосходит порог. Каждое из полученных подмножеств делится далее аналогичным образом.

Далее проводится локализация мод гистограммы. Вначале каждому вектору на основе анализа его ближайших соседей ставится в соответствие градиент. Вектору приписывается номер вектора с максимальным значением градиента. Если градиент меньше нуля, то это означает, что координаты вектора являются координатами локального максимума и вектору приписывается его собственный номер. В итоге каждой моде гистограммы сопоставляется ориентированный граф, корень которого соответствует точке моды. Если количество получаемых кластеров (количество локальных максимумов гистограммы) больше заданного порога, то проводится сглаживание гистограммы. Сглаживание осуществляется либо путем замены частоты Н(х) вектора х на среднее значение частот его ближайших соседей, либо путем уменьшения "разрешения" векторов данных, т. е. делением компонент векторов на 2.

На завершающем этапе выполняется раскраска ориентированного графа одним цветом, т.е. всем вершинам графа присваивается значение, которое присвоено его корню.

Рис. 4

На рис. 4 приведен результат кластеризации исходного изображения (рис. 2) описанным методом. Выделено 15 кластеров. Сравнивая рис. 3 и 4, можно заметить существенные различия: на рис. 3 лучше разделены сельскохозяйственные угодья, а на рис. 4 лучше "проработаны" водные поверхности (включая зоны паводкового затопления). Эти различия говорят о том, что надежность результатов кластеризации часто можно оценить лишь сравнением нескольких вариантов обработки.

4. Заключение

Практика решения конкретных прикладных задач дистанционного зондирования Земли с использованием предлагаемого подхода подтверждает его эффективность [2].

Работа выполнена частично при финансовой поддержке РФФИ (грант №07-07-00085).

Список литературы

[1] Remote Sensing: The Quantitative Approach // Edited by P.H. Swain and S.M. Davis. USA, McGraw-Hill, Inc., 1978. - 396 p.

[2] Асмус В. В. Программно-аппаратный комплекс обработки спутниковых данных и его применение для задач гидрометеорологии и мониторинга природной среды. Дис. д-ра физ.-математ. наук. На правах рукописи. - М., 2002. - 75 с.

[3] Асмус В. В., Бучнев А. А., Пяткин В. П. Программный комплекс для обработки данных дистанционного зондирования Земли. Труды XXXII Международной конференции "Информационные технологии в науке, образовании, телекоммуникации и бизнесе, IT+SE'2005", приложение к журналу "Открытое образование", 20-30 мая 2005, Украина, Крым, Ялта-Гурзуф. С. 229-232.

[4] Дж. Ту, Р. Гонсалес. Принципы распознавания образов. - М.: Мир, 1978. - 411 с.

[5] Гонсалес Р., Вудс Р. Цифровая обработка изображений. - М.: Техносфера, 2005, - 1072 с.

[6] Marques de Sa J.P. Pattern Recognition: Concepts, Methods and Applications. SpringerVerlag, Berlin, Heidelberg, 2001, - 318 р.

The Cluster Analysis and Classification with Training of Multispectral Data of Earth Remote Sensing

Vasily V. Asmus", Alexey A. Buchnev6 and Valery P. Pyuatkin6

a Research center "Planeta" 7 Bolshoy Predtechensky, Moscow, 123242 Russia b Institute of Computational Mathematics and Mathematical Geophysics SB RAS

6 pr. Ak. Lavrentjeva, Novosibirsk, 630090 Russia

It is obtained questions, connected with the problem of choosing appropriate algorithms of recognition of multispectral data of Earth remote sensing. It is submitted the system of supervised classification, based on a strategy of maximum probability for vectors of indications having the normal distribution. It is described the system of cluster analysis, including an algorithm for K-means method and analyzing method of mode of multidimensional histogram.

Key words: Earth remote sensing, data recognition, supervised classification, unsupervised classification, cluster analysis, decision rule, classifier training, K-means method, multidimensional histogram.