CC BY
133
УДК 528.852
А.А. Бучнев, В.П. Пяткин
СГГ А, Новосибирск
КОНТРОЛИРУЕМАЯ И НЕКОНТРОЛИРУЕМАЯ КЛАССИФИКАЦИЯ МНОГОСПЕКТРАЛЬНЫХ ДАННЫХ ДИСТАНЦИОННОГО ЗОНДИРОВАНИЯ ЗЕМЛИ
Эффективность дистанционных исследований Земли из космоса во многом определяется используемыми методами тематической обработки данных дистанционного зондирования Земли (ДДЗЗ). При этом центральным при тематической обработке, безусловно, является модуль распознавания [1].
В данной работе представлены алгоритмы контролируемой и неконтролируемой классификации многоспектральных ДДЗЗ, включенные в модуль распознавания программного комплекса обработки ДДЗЗ, разработанного в ИВМиМГ СО РАН совместно с НИЦ “Планета” Роскомгидромета РФ [2].
1. КОНТРОЛИРУЕМАЯ КЛАССИФИКАЦИЯ
МНОГОСПЕКТРАЛЬНЫХ ДДЗЗ
Система контролируемой классификации (классификации с обучением) ДДЗЗ в программном комплексе состоит из семи классификаторов (один поэлементный классификатор и шесть объектных), основанных на использовании байесовской стратегии максимального правдоподобия, и двух объектных классификаторов, основанных на минимуме расстояния [3].
1.1. Поэлементная классификация.
Под элементом здесь понимается N -мерный вектор измерений (признаков) х = (х,•••,)Т , где N - число спектральных диапазонов. Предполагается, что векторы х имеют в классе щ нормальное распределение N(т, В) со средним т и ковариационной матрицей В . В этом случае байесовская стратегия максимального правдоподобия для поэлементного классификатора формулируется следующим образом.
Пусть О = (щ,...,щ) - конечное множество классов; р(щ) - априорная
вероятность класса щ. Тогда дискриминантная функция класса щ имеет вид:
gi (х) = !п( р(щ)) - 051п(1 В 1) - 0-5( х - т)Т В-1 (х - т)- С1)
Классическое решающее правило для классификатора принимает следующий вид: вектор х заносится в класс щ, если gl (х) > g (х) для всех у * I.
Для проверки гипотезы о нормальности обозначим через Т основанное
2
на распределении % пороговое значение для векторов класса щ :
Т = 1п р(щ ) - 0.5 А( N, О) - 0.51п| В |, (2)
где А(N, О) - критическое значение уровня Q распределения %2. Пусть I - переменная, значение которой зависит от параметра классификатора №г:
— да, если гкт = 1,
Т, і “ если гкт = 2,
т шіп Т, і=1 і ’ если гкт = 3,
т шах Т, і і =1 если гкт = 4,
т )/т і =1 если гкт = 5.
(3)
Тогда решающее правило для классификатора с учетом (3) принимает следующий вид [1]: вектор х заносится в класс а1, если gi (х) > gJ (х) для
всех і ф і и g¡ (х) > г . В противном случае вектор заносится в класс отклоненных векторов (класс с номером т +1).
Значения А(р, О) для размерности вектора р < 30 находятся из статистических таблиц. Например, для р = 4 и уровня О = 0.05 (т.е. 5% векторов может быть отклонено) А=9.488. При р > 30 для нахождения значений А( р,О) используется аппроксимация:
Ж1-0(Р) - Р
1 -
2
9 р
+ и
1-е
_2_
9 р
(4)
где и
1-е
значение стандартизованной нормальной величины для
вероятности 1 - Є (в частности, для уровня Є = 0.05 значение и095 * 1.645).
1.2. Объектная классификация.
Под объектом мы понимаем блок смежных векторов квадратной или крестообразной формы с длиной образующей ^ Поскольку физические размеры реально сканируемых пространственных объектов, как правило, больше разрешения съемочных систем, между векторами данных существуют взаимосвязи. Использование информации подобного рода дает возможность повысить точность классификации, если пытаться распознавать одновременно группу смежных векторов - объект в приведенном выше смысле. Решение об отнесении центрального элемента объекта к тому или иному классу принимается на основе результата классификации всего объекта.
Покажем размещение векторов в объектах для F=3:
для квадрата.
(5)
Классификаторы ОШОС, ОМЕАМС, ОМЕАКС, OMARKC работают с блоком в форме перекрестия, а классификаторы OINDS, ОМЕАМЗ,
^ х5 Л ^ х X х3 ^
X = Х4 х! х2 - для перекрестия, X = х 4 5х х 6
V х3 V х7 х8 х9)
г
<
3
ОМЕА^, OMARKS работают с квадратным блоком. Приведем алгоритмы работы этих классификаторов.
Пусть снова ^ = (щ,...,щ) - конечное множество классов, р(щ) -
априорная вероятность класса щ, X = (х1,..., хк) - объект, состоящий из к N -мерных векторов х1 = (х[,. х‘ы ).
ОГЫОС и OINDS основаны на предположении, что векторы внутри блока независимы. Будем рассматривать векторы, составляющие объект X, как один вектор размерности кЫ. Тогда дискриминантная функция класса щ имеет вид:
к
ё,(х) = ¥р(щ))- 0.5к 1п(| В1 |)- 0.5^(х1 - т )ГВ;1(У - т). (5а)
1=1
Решающее правило для данных классификаторов принимает следующий вид: центральный элемент объекта X заносится в класс щ , если ё (X) > ё (X) для всех у ф , и ё, (X) > I. В противном случае центральный элемент объекта заносится в класс отклоненных векторов. Здесь Г определяется по (2) и (3) с А = A(kN, О).
ОМЕАМС и OMEANIS классифицируют среднее блока в предположении, что векторы внутри блока независимы. Классифицируется вектор х, равный среднему по всем векторам объекта X:
1 к
х = - 2 х1 . (6)
к 1=1
Дискриминантная функция класса щ имеет вид:
ё,(X) = 1п(р(щ))-0.51п(| В |)-0.5к(х -т)ГВг-1(х -т). (6а)
Решающее правило для данных классификаторов имеет следующий вид: центральный элемент объекта X заносится в класс щ , если ё (X) > ё} (X)
для всех у ф , и ё, (X) > tl. В противном случае центральный элемент объекта
заносится в класс отклоненных векторов. Здесь ti определяется по (2) и (3).
OMEANC и ОМЕА^ классифицируют средний вектор (6) блока в предположении, что векторы внутри блока независимы и ковариационные матрицы равны единичной. Фактически это объектные классификаторы, решающие правила которых основаны на минимуме расстояния до центра класса. Дискриминантная функция класса щ имеет вид:
ё!(X) =1п(р(щ)) - 0.5(х - т)Г(х - т,). (6б)
Решающее правило для данных классификаторов имеет следующий вид: центральный элемент объекта X заносится в класс щ , если ё, (X) > ё} (X)
для всех у ф , и ё, (X) > tl. В противном случае центральный элемент объекта заносится в класс отклоненных векторов. Здесь ? = - А, где А>0 - число,
задаваемое пользователем.
OMARKC и ОМАКО объединяют пространственные характеристики объекта на основе модели каузального марковского случайного поля первого
и третьего порядка соответственно. Предполагается, что автокорреляционная функция элементов объекта разделима по вертикали и горизонтали:
согг (х., хы) = рАИ • р]н , (6в)
где х - вектор признаков на пересечении I -ой строки и & -ого столбца; Р - коэффициент пространственной корреляции между двумя соседними векторами в горизонтальном направлении;
Р - коэффициент пространственной корреляции между двумя соседними векторами в вертикальном направлении.
Пусть I - номер класса, 1,3 - размеры образа. Тогда:
1 М
р„=^ !р;„ (7)
N :=1
где
1 I 3-1
,,,_п I кх,:- <)(</.1 - ю
р'„ = 1 (3 1) '•1 ^---------------------, а 1 < : < N. (7а)
Соответственно
1 I 3-1
т: = ---£ Х!х: ,
/г I (3 -1) £ £ 3 ’
1 I 3
тГ = 77^7££х',, (7б)
I (3 - 1) /=1 3=2
1 I 3-1
^=„ , л 1 ££(х: - т:)2 >
I (3 - 1) - 1 /=1 з=1
1 I 3 Уаг: =-------1----£ £ (хк - т: )2.
Гг Г/ Г 1 ¿—¡¿—1у
I (3 - 1) - 1 /=1 з=2
Аналогичные формулы получаются для ри.
Для классификатора OMARKC (форма объекта - перекрестие) имеем решающие функции вида (см. (5))
5
с,(х) = g.(х1) - о.5£ g;(х. |х1). (7в)
3=2
Здесь gl (х1) - дискриминантная функция (1) поэлементного
классификатора, а в качестве примера для функций g!(ху.|х1) распишем функцию g!( х2|х1) (остальные получаются аналогично):
g! (хг\ х1) = N 1п(1 - р1) + 1п| Б\ + -^т( Ж 2 - 2р„Д 1 + РЖ),
1 - Ры
(8)
где
Ж’к1 =(х:- т)Тв:\х - т). (9)
Центральный элемент х1 объекта X относится в класс щ , если О(X) > О (X) для всех & ф , и О(X) > Тгс, где Т' - пороговое значение для отклоненных векторов класса щ :
ТС = 1п р(щ) - 0.5( А( N1,0) + (2[Ь/2]+ 1)1п|В,| + N[1 /2]1п((1 -р )(1 -рЫ))). (10)
Здесь А(N1,0) - критическое значение уровня 0 распределения х2, а Ь - количество векторов в объекте.
Для классификатора OMARKS (форма объекта - квадрат) имеем следующие решающие функции (см. (5)):
О, (Х) = gf (х^х2 , х4 , х5 ) + gг2 (х^ х3 , х5 , хб ) + gг2 (х^ х5 , х7 , х8 ) + g^ (х5 |хб , х8 , х9 ) +
+ g! (хз | хб ) + g! (хб | х9 ) + g! (х7\ х8 ) + g) (х^ х9 ) + gi (х9 ), (1 0а)
где g¡ (х9) - дискриминантная функция (1) поэлементного
классификатора, g!(х |х;) - функции вида (8), а в качестве примера для функций gl2 распишем функцию gг2(х^, х4, х5) (остальные получаются аналогично):
1
(1 -Р2)(1 -Р2)
ё] (х1 | х2 , х4 , х5 ) = -0. 5(1п| Д | + N 1п((1 - р;2 )(1 - р1щ )) + -------------------~----------— (Ж 1 - 2р; Д!2 -
- 2рЖ + 2рк,Р«Ж\5 + РЫЖ22 + 2рк,рмЖ24 - 2рЫрЖ5 + р2Ж44 - 2рИ,р1Ж45 + Рй Д5)) .
(11)
Центральный элемент х5 объекта X относится в класс щ , если О (X) > О (X) для всех & ф, и О (X) > Ть. , где Ть. - пороговое значение для отклоненных векторов класса щ :
т; = Тс -0.5(Ь0-5 -1№ 1п((1 -Р,2,)(1 -р’)) + 1п|В,|), (11а)
(Тгс берется из (10)).
Когда размер объекта равен 1, все объектные классификаторы становятся поэлементными классификаторами (1), основанными на стратегии максимального правдоподобия, за исключением OMEANC и ОМБА№, которые становятся поэлементными классификаторами, основанными на минимуме расстояния.
1.3. Обучение и работа классификаторов.
Необходимые для построения дискриминантных функций классов оценки статистических характеристик векторов, средних ковариационных матриц, коэффициентов пространственной корреляции между значениями координат соседних векторов в горизонтальном и вертикальном направлениях - определяются на основе векторов из обучающих выборок (полей). Кроме обучающих для каждого класса может быть задан набор контрольных полей.
Оценки коэффициентов пространственной корреляции ры и ри в /-ом классе вычисляются как средние по количеству обучающих полей:
Р; = (£ А,)/ 5, (11б)
I = 1
где 5 - количество обучающих полей в классе, р1ы - коэффициент пространственной корреляции в горизонтальном направлении для обучающего поля с номером I, вычисленный согласно (7).
Все классификаторы могут использоваться в двух режимах - тестовом и рабочем. По результатам работы классификаторов в тестовом режиме над векторами обучающих и контрольных полей рассчитываются матрица ошибок и оценки вероятностей правильной классификации. Известно, что эти оценки для векторов из обучающих полей являются, в среднем, оптимистическими, а для векторов из контрольных полей также в среднем являются пессимистическими. Анализируя эти данные, можно оценить (проконтролировать) качество обучения.
Как отмечалось выше, возможны ситуации, при которых нарушается условие о случайности поля измерений, следствием чего являются нулевые дисперсии в некоторых каналах. Тогда формулы типа (1) становятся неприменимыми, т.к. ковариационные матрицы В вырожденные. Для исправления подобных ситуаций в системе классификации предусмотрена функция добавления к спектральным каналам с нулевой дисперсией гауссовского шума с нулевой средней и единичной дисперсией.
Результатом работы классификаторов в рабочем режиме является одноканальное (байтовое) изображение, значениями пикселов которого являются номера классов. Это изображение окрашивается в предопределенные цвета, которые в интерактивном режиме могут быть заменены на цвета, определяемые пользователем. Кроме того, к этому изображению можно применить функцию постклассификации для удаления изолированных пикселов.
Система контролируемой классификации имеет следующие характеристики: число обучающих образов - до 9, число классов - до 15, число обучающих и контрольных полей в классе - до 10, размер каждого поля - до 50x50 векторов, размер объекта - от 1x1 до 11x11, размерность векторов данных не ограничивается.
На рис. 1 представлена тематическая карта ледовой обстановки в восточном секторе Арктики, полученная с использованием контролируемой классификации. В качестве распознаваемых изображений использовались мозаики радиолокационных и радиометрических снимков с ИСЗ «Океан-О1». Для выбора тестовых участков использовалось цветосинтезированное изображение.
Рис. 1
2. НЕКОНТРОЛИРУЕМАЯ КЛАССИФИКАЦИЯ ДДЗЗ
Неконтролируемая классификация (кластерный анализ) в программном комплексе представлен двумя алгоритмами - методом К -средних и методом анализа мод многомерной гистограммы [3].
2.1. Метод К - средних.
Первый подход основан на итеративной процедуре отнесения векторов признаков классам по критерию минимума расстояния от вектора до центра класса. Оптимальным считается такое разбиение входных векторов на кластеры, при котором внутриклассовый разброс не может быть уменьшен при переносе какого-либо вектора из одного кластера в другой.
Алгоритм состоит в выполнении следующих шагов:
1. На основе заданного соотношения а чистых и смешанных векторов производится разделение векторов на чистые и смешанные. Под смешанными мы понимаем векторы, компоненты которых либо формируются за счет попадания в поле зрения съемочной аппаратуры нескольких объектов, либо искажены влиянием фона. С этой целью вначале для исходного набора векторов измерений с помощью многомерного оператора Робертса рассчитывается градиентное изображение и одновременно строится
гистограмма градиентов. Исходя из заданного а , по гистограмме определяется порог, разделяющий векторы на смешанные и чистые.
2. Объединение чистых векторов в связные компоненты. На этом этапе все чистые векторы объединяются в связные компоненты, которые последовательно нумеруются. Соответствующий алгоритм, идейно близкий к алгоритму заполнения областей с произвольной границей по критерию связности [3], может выделять и нумеровать одновременно любое количество многосвязных областей без ограничений на их форму и ширину контуров. Для каждой связной компоненты вычисляется вектор средних значений.
3. Итеративная кластеризация векторов средних значений. Начальные центры кластеров определяются по следующей схеме. В качестве первых двух центров берется пара векторов, наиболее далеких друг от друга. Затем вся выборка делится на кластеры по критерию близости к выбранным центрам. В каждом кластере отыскивается вектор, наиболее далекий от центра. Для всех таких векторов рассчитывается суммарное расстояние до центров. В качестве нового центра берется вектор, для которого суммарное расстояние максимально, и процедура распределения векторов по кластерам повторяется.
4. Распределение связных компонент по кластерам. На этом этапе связные компоненты получают новые номера. Новый номер присваивается компоненте в соответствии с номером кластера, в который попал вектор средних значений этой компоненты.
5. Кластеризация смешанных векторов. На завершающем этапе производится неитеративная кластеризация смешанных векторов по принципу минимума расстояния до центров кластеров С,...,Ск. Смешанный вектор г будет отнесен к ближайшему кластеру щ, если
ІІС -г|| < 0.5тах||С -С\|, і,I =1,...,к, і ф I.
Рис. 3 иллюстрирует результат работы кластеризации методом К -средних значений изображения, представленного на рис. 2. Рис. 3 демонстрирует результат работы алгоритма для следующих входных данных: количество выделяемых кластеров 16, соотношение количества смешанных и чистых векторов а = 0.35.
Рис. 2
Рис. 3
2.2. Метод анализа мод многомерной гистограммы.
В основе второго подхода лежит предположение, что исходные данные являются выборкой из многомодового закона распределения, причем векторы, отвечающие отдельной моде, образуют кластер. Таким образом, задача сводится к анализу мод многомерных гистограмм. В программный комплекс включена реализация метода, описание основных шагов которого приведено в [3].
Гистограмма генерируется последовательным просмотром векторов данных и сравнением каждого вектора с текущим списком векторов. При этом либо изменяется соответствующее значение частоты, либо вектор добавляется в список. Для вычисления адресов векторов в списке используется хэш-кодирование. Первым шагом модального анализа является поиск ближайших соседей данного вектора списка среди других векторов списка. По определению вектор х есть ближайший сосед вектора у, если |хг - у.| < 1 для г = 1,_, N. Каждый из возможных ближайших
соседей данного вектора х может быть получен из него прибавлением вектора сдвига, компоненты которого принимают значения из множества {-1,0,1} . Алгоритмически г -ый вектор сдвига г = 1,2,..^ -1 можно получить, уменьшив на 1 каждый из коэффициентов представления числа г в троичной системе исчисления. Поскольку в реальной гистограмме присутствуют далеко не все ближайшие соседи, то для эффективного их поиска векторы предварительно упорядочиваются в многомерные бинарные деревья. В этом случае время поиска всех ближайших соседей данного вектора становится пропорциональным числу реально существующих соседей. При построении дерева векторы х рассматриваются как N -мерные ключи. Вначале рассчитываются
дисперсии по всем координатам векторов, и определяется координата ] , имеющая максимальную дисперсию. Медианное значение выборки по этой координате используется в качестве ключа для разделения множества векторов на два подмножества: в одно подмножество помещаются векторы, значение которых по координате ] меньше порогового значения, а в другое векторы, у которых значение координаты превосходит порог. Каждое из полученных подмножеств делится далее аналогичным образом.
Далее проводится локализация мод гистограммы. Вначале каждому вектору на основе анализа его ближайших соседей ставится в соответствие градиент. Вектору приписывается номер вектора с максимальным значением градиента. Если градиент меньше нуля, то это означает, что координаты вектора являются координатами локального максимума и вектору приписывается его собственный номер. В итоге каждой моде гистограммы сопоставляется ориентированный граф, корень которого соответствует точке моды. Если количество получаемых кластеров (количество локальных максимумов гистограммы) больше заданного порога, то проводится сглаживание гистограммы. Сглаживание осуществляется либо путем замены частоты к( х) вектора х на среднее значение частот его ближайших соседей, либо путем уменьшения “разрешения” векторов данных, т.е. делением компонент векторов на 2.
На завершающем этапе выполняется раскраска ориентированного графа одним цветом, т.е. всем вершинам графа присваивается значение, которое присвоено его корню. На рис. 4 приведен результат кластеризации исходного изображения (см. рис. 2) описанным выше методом. Выделено 13 кластеров при заданном максимальном количестве 16.
Существенные различия в рисунках 3 и 4 говорят о том, что надежность результатов кластеризации часто можно оценить лишь сравнением нескольких вариантов обработки.
Работа выполнена частично при финансовой поддержке Российского фонда фундаментальных исследований (проект № 07-07-00085).
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Природа Земли из космоса // А.П.Тищенко, В.В. Асмус, В.П. Пяткин и др. Л.: Гидрометеоиздат, 1984. - 151 с.
2. Асмус В.В., Бучнев А.А., Пяткин В.П. Программный комплекс для обработки данных дистанционного зондирования Земли // Труды XXXII Международной конференции “Информационные технологии в науке, образовании, телекоммуникации и бизнесе IT+SE,2005”, 20-30 мая 2005 г., Украина, Крым, Ялта-Гурзуф. - с. 229-232.
3. Асмус В.В. Программно-аппаратный комплекс обработки спутниковых данных и его применение для задач гидрометеорологии и мониторинга природной среды
// Диссертация (научный доклад) на соискание ученой степени доктора физикоматематических наук. На правах рукописи. Москва - 2002. - 75с.
© А.А. Бучнев, В.П. Пяткин, 2007
