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

УДК 519.248

Е.А. Барышева, А.К. Протасов, М.С. Шестакова

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

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

Введение

При разработке моделей дистанционного исследования Земли важную роль играет проблема поиска устойчивых взаимосвязей энергетических, спектральных, геометрических и биофизических параметров земных покровов при формировании образов растительных ценозов в условиях фенологических вариаций природных комплексов. Разработка физико-математических уравнений, описывающих связи спектральных яркостей и биометрических параметров растений на основе законов взаимодействия электромагнитного излучения с веществом, является очень сложной теоретической задачей, требующей знания большого количества параметров, характеризующих физическое состояние земных покровов. В то же время новые разделы кибернетики, основанные на современных методах математической статистики, связанные с непараметрическими методами извлечения информации из экспериментальных данных, позволяют решать широкий класс задач распознавания сложных образов, восстанавливать стохастические зависимости между переменными и оценивать качество моделей. С помощью регрессионного анализа на основе дистанционных и согласованных с ними наземных наблюдений достаточно просто можно находить устойчивые связи для системы «почва — растительность» и восстанавливать такие параметры, как высота и густота растительности, оценивать проективное покрытие почвы растительностью, влажность и состав почвы, определять объем фитомассы, биопродуктивность и другие характеристики растительности.

При построении решающих правил распознавания образов и обнаружения событий наиболее предпочтительным в силу высокой формализации является статистический байесов подход. Именно в рамках этого подхода при полной априорной информации получены оптимальные решающие правила из условий минимума функционала потерь, называемого риском. Информация о реальных ситуациях часто задана в виде выборки наблюдений, возможно, классифицированных учителем. Качество процедур классификации, использующих непараметрические оценки неизвестных условных функций плотности, естественно характеризовать оценкой средних потерь по выборке материала обучения с помощью функционала эмпирического риска. Эффективность решающих правил существенно возрастает, если система используемых признаков содержит контекстную информацию. В простейшем случае контекстная информация содержится в самом же изображении, а именно в окрестности анализируемой точки, выявляемой системой текстурных признаков. Так как текстура — пространственное локальное свойство изображения, измерение её характеристик должно быть ограничено областями относительной однородности, поэтому необходимо установить границы области квазиоднородности текстуры, например, с помощью метода автоматической сегментации. Следует заметить, что при всех положительных качествах текстурных признаков их основной недостаток заключается в том, что достаточно часто при анализе изображений в рамках скользящего окна будут попадаться не только однородные текстуры какого-либо одного класса, но и комбинации классов и границы разделов различных текстур. В этом случае значение текстурного признака несет в себе смешанные свойства объектов подстилающей поверхности Земли и, чтобы усилить роль центральных пикселов, относительно которых производятся те или иные оценки, следует расширить набор текстурных признаков добавлением спектральных компонентов для этого пиксела и разнообразных вегетационных признаков.

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

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

Алгоритм кластерного анализа на основе модификации ЕМ-подхода декомпозиции смешанного распределения

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

ЕМ-алгоритм предполагает организацию итерационного процесса, состоящую из последовательности следующих шагов: ОЦЕНИВАНИЕ, при этом оцениваются параметры распределений; МАКСИМИЗАЦИЯ, при этом вычисляется значение критерия качества и процедура повторяется до получения приемлемых значений критерия качества.

Рассмотрим вариант простого ЕМ-алгоритма кластерного анализа при следующих предположениях [1-3]:

1) исходные данные в виде смешанной выборки получены из генеральной совокупности с априорно известным количеством С классов и имеют вид векторных наблюдений (без указаний на принадлежность к классам) х^-.^хдг, хейл, N — объем смешанной выборки;

2) априорные вероятности классов = Р(сор , у = 1,...,С, не известны;

3) известен параметрический вид одномодальных условных функций плотности, являющихся моделями классов (кластеров) /(х | соу-, в -), у = 1,...,С ;

4) неизвестными параметрами в задаче также являются векторные параметры распределений 6,-, 0 е Ит одной и той же размерности для разных классов.

Функция плотности вероятности распределения отдельного наблюдения имеет вид смеси

Ах I в) = X /-

/=1 ¡=1

где 0 = (01,...,вс) — параметры соответствующих одномодальных функций плотности. Требование одномодальности обусловлено тем, что наша задача заключается в выявлении текстурно-однородных элементов многоспектральных изображений при недостатке априорной информации о выделяемых классах.

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

/V

^„.„Хдг |в;.) = П^х*1

Если все неизвестные векторные компоненты 0 удастся оценить по смешанной выборке х1,...,хдг, то мы сможем декомпозировать смесь на составные одномодальные распределения и выделить унимодальные составляющие, которые ответственны за классы однородности данных. Предполагается, что смесь (2) идентифицируема [1-3]. Это трудно формализуемое условие на содержательном уровне означает, что все локальные моды «достаточно ярко выражены» и тем или иным способом их можно выявить. Предполагая дифференцируе-мость ,...,Хдг) в (2), получим необходимые условия для оценивания 0, вычисляя градиент выражения

ь = 1>ЯхА|е)

1

следующим образом:

N

-V,

у=1

Если при I ф у элементы векторов 04 и 0у функционально независимы и ввести апостериор ные вероятности по формуле Байеса вида

о/ I 1 = \хкА)

то необходимое условие экстремума (3) по каждой из компонент б1 вектора 0г, / = 1,...,С будет иметь вид

= X I 1©1»в£) = 0, ¿ = 1 ,...,с.

Действительно, подставим (5) в (4), получим

V ь у

01 01

У=1

но

У01 {1п/(хА |

откуда, после подстановки (8) в (7), и получаем (6).

N

Теперь найдем необходимые условия для оценивания р1, дифференцируя I = £ 1п /(хА | в).

с

С учетом (5) Дх,, | 0) = и ограничений р1 = 1 , вводимых в функцию (3) с

Р(со;|х^,0г) ¿=1

помощью множителей Лагранжа, получим

— = I у

1

дь _ *

Р

/(х^, I СО;,0;) - ХЫ = 0 ,

¿Р

А = 1 Рь

Так как множитель X входит в произведение Хр> , то вводя переобозначение, получим

.V

I

к = 1

¿Р(о>, |хл>е,)-ЛГЛ=0

к=1

Для дальнейшего конкретизируем вид условных функций плотности, которые являются моделями кластеров.

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

Учитывая тот факт, что модели текстурно-однородных участков многоспектральных изображений описываются одномодальными распределениями, естественно, по крайней мере

на начальном этапе, обратиться к гауссовым аппроксимациям неизвестных функций плотности. Кроме того, большие объемы выборочных данных, появляющиеся в задачах структурирования аэрокосмической информации, не позволяют решать задачу декомпозиции смесей при большом числе неизвестных. Поэтому будем предполагать неизвестными лишь априорные вероятности кластеров р] = у?(со;), у = 1,...,С, и векторные математические ожидания 0^. = ц. 5 у = 1,...,С. Ковариационные матрицы будем считать приблизительно одинаковыми и «похожими» на ковариационную матрицу, оцениваемую по смешанной выборке х^.^х^ объема N :

N

^ N 1 ^

G = — X - £)(х* - , А = ^ X ,

И=\ к=1

причем для каждого кластера ковариационная матрица регулируется дисперсией кластера

G,=—G, у = 1,..., С,

В этом случае плотность смеси имеет вид

-1/2

m 1 -

х ехр

N С /е-11-1

В случае неизвестных математических ожиданий

л/2

-I

lgf(xft I (o¿,£¿) « -lg

(2n) f \G:

¡1/2

Дифференцируя квадратичную форму (14) по вектору щ , получим следующее выражение [4]:

V 1 gf(xk | соi,\ii) = Gi1(xk-\ii).

Подставив в уравнение (6), получим

X P(coi \ х^)^** - щ)] = О, i = 1 ,...,С

к = 1

Умножая на слева, находим оценку вектора математического ожидания

Йг

JV

fe=l_

N

k=l

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

N

Aifí + l]

fr=i

Л = 1

Аналогично из (10) получаем рекуррентное выражение для нахождения р1:

у=1

Теперь основной задачей использования выражений (18) и (19) является нахождение хороших приближений для параметров {ц} и {р}. Предварительно заметим, что в случае близости пар значений и цу их можно объединить по критерию

Теперь вернемся к задаче оценивания локальных экстремумов смешанного распределения в предположении, что моды смеси достаточно ярко выражены. Необходимо найти приближенные оценки локальных мод функции плотности /(х) , эти моды будут позиционировать экстремумы одномодальных распределений /(х|сог,04), £ = 1,...,С. Воспользуемся следующим конструктивным способом оценивания локальных мод многомерной многоэкстремальной функции плотности смеси (1) [5]. Определим к -мерный шар, содержащий небольшое количество тп выборочных точек, так что т «Ы . Затем центр шара помещаем последовательно в каждую из выборочных точек и подсчитываем число соседствующих с центральной точкой элементов выборки, попавших в шар. Если для некоторой исследуемой точки количество попавших в ее окрестность точек, определенную шаром, больше некоторого среднестатистического значения т или максимально, то в этой точке и будет локальный экстремум функции плотности, а сама точка будет приближенной оценкой локальной моды [5]. Точнее, локальную моду можно оценить, вычислив математическое ожидание точек, попавших в окрестность шара. Для целей поиска кластерообразующих выборочных данных (под кластерообразующими точками мы понимаем те из выборочных данных, окрестности которых содержат максимальное количество точек) можно воспользоваться методом декомпозиции гистограмм, который заключается в следующем. Для каждой выборочной точки из х^-.^х^ строится гистограмма распределения расстояний от этой точки, называемой опорной, до всех точек выборки. В первом приближении естественно воспользоваться простейшей метрикой Евклида. Перебирая последовательно все точки выборки, наблюдая распределение расстояний, выбираем ту из опорных кластерообразующих точек, гистограмма которой левоориентирована в большей степени. Это и означает, что выбранная точка служит оценкой локальной моды. Для построения байесова решающего правила обнаружения и выделения экстремальных (минимальных) значений расстояний на полученной гистограмме распределения расстояний до всех точек выборки необходимо прежде всего восстановить вероятностные модели ситуаций А0 — ядро кластера и Аг — не входящие в кластер точки и оценить их априорные вероятности. Для этого необходимо декомпозировать полученную гистограмму распределения расстояний на два распределения, одно из которых распределение расстояний до подвыборки кластера, а второе — распределение расстояний до остальных не принадлежащих кластеру точек выборки, то есть надо идентифицировать составляющие

где /0(х) _ распределение расстояний от некоторой центральной точки до точек кластера;

д (х) _ распределение расстояний до удаленных точек; Р,Я — априорные вероятности ситуаций А0 и А1 соответственно, причем Р + ф =1. Возникает задача оптимизации квадратичного критерия качества идентификации смеси (21) следующего вида:

где fi,xj) — гистограмма распределения всех расстояний от пробной точки до всех векторов выборки, а в — вектор неизвестных параметров, состоящий из компонента Р и параметров функций плотности /о(х) и Д(х), принадлежащих параметрическим семействам функций. Следует заметить, что задача восстановления компонентов смеси имеет решение лишь в случае её идентифицируемости. Это трудно формализуемое и проверяемое условие с геометрической точки зрения означает, что ^(х) и ^(х) должны иметь ярко выраженные моды. Поэтому доля или мера участков с экстремальными значениями градиентов должна быть достаточно высокой, чтобы функция плотности /0(х) могла проявить свою форму. Задача декомпозиции смешанного распределения ввиду высокой неопределенности не всегда имеет решение. Необходимо привлечь априорные данные о форме составляющих распределений. В связи с этим обратимся к следующему факту математической статистики, связанному с теорией экстремальных значений. Известно, что плотность распределения минимумов N независимых случайных величин в асимптотике растущего числа наблюдений N -> оо I типа (распределение минимальных значений Гумбеля [6]) имеет следующий вид:

следую ГЦ1

Ux) = f{x;\i,o) = -exp а

-{x-\i)-e

х-ц

-оо < х < +оо , -со < ц < +оо , а > 0, ГДе ц — параметр (мода) центра распределения; сг — масштаб распределения, причем стандартные оценочные математическое ожидание \х и дисперсия ст связаны с параметрами распределения /0(х) ц и су следующим образом: ц = Д-0,577а, о = 1,283а. В качестве нами было выбрано распределение Джонсона с параметрами с (нижняя граница х), X (размах выборки), г|, у (параметры формы). Учитывая, что часть параметров можно оценить по выборочным данным [6], фактически вектор неизвестных параметров имел лишь три компонента и 0 = (Р,х\,у)т , где Т — знак транспонирования. Задача оптимизации критерия (22) решалась с привлечением адаптивных методов поиска экстремума Я.З. Цыпкина [7]. После того как смесь идентифицирована, построим байесово решающее правило проверке кластер. Байесово решающее правило выявляет

ки двух гипотез:

Но

на выборке х1,...,х7

— кластер и Н1 кластерообразующее подмножество х},..., х^

и = агщта.х{Р^(х),Я^(х)},

Юл)

N, « N ,

где и — принимаемое решение или номер принимаемой гипотезы; и <= {0,1} и это подмножество в дальнейшем исключается из анализа. Повторяя этот процесс, находим оценки локальных мод по всей выборке. После того как «центры» предполагаемых кластерных претендентов выявлены, можно оценить обобщенную ковариационную матрицу по всей выборке. Следует отметить, что на этом этапе количество классов должно быть существенно больше того значения числа классов, которое должно быть задано как окончательное. Теперь, когда центры — математические ожидания кластеров — определены, воспользуемся гауссовой аппроксимацией с обобщенной ковариационной матрицей для описания вероятностных моделей кластеров и построим байесово решающее правило распознавания пикселов, не вошедших в материал, по которому формируется обучение. В этом случае байесово решающее правило последующего агрегирования выборочных данных в кластеры имеет вид

где

f(x \ со,) =

arg max P(a)!)f(x \ со,), 1=1,..„с

I ^

1,2

1=1

С — оценка ковариационной матрицы с использованием всей смешанной выборки; — собственные значения ковариационной матрицы, / = 1,...,к .

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

Следующим шагом работы алгоритма является укрупнение классов путем слияния близких кластеров. В качестве меры близости кластеров можно взять расстояние Бхатача-рия для равных ковариационных матриц

Ц,)Т -С'1

-(^-ц,), i ф у, U] = 1,...,с',

где С' > С — промежуточное количество кластеров.

Смысл следующего приближения заключается в том, чтобы уточнить моды классов и оценить ковариационные матрицы классов с учетом того, что эти выборки засорены. Так как выборочные данные кластеров изменились, вновь можно оценить математические ожидания классов и ковариационные матрицы.

На втором шаге работы алгоритма кластеры уже проявляют свою индивидуальность и классификацию данных естественно провести байесовым решающим правилом с гауссовыми функциями плотности, для которых подсчитаны ковариационные матрицы классов Gi,

I = 1,...,С'. Дальнейшее укрупнение классов естественно производить по критерию минимума расстояния Бхатачария с учетом различных ковариаций

1 Ч\ 1 / \Т Сг{ + (г- , N. 1

Полученные оценки математических ожиданий центров кластеров теперь могут быть уточнены по рекуррентному соотношению (18), а априорные вероятности — по (19).

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

Ох-Вт

тогда относительный критерий качества кластеризации имеет вид

Я = -77---,

Л \\*3(х,у)-Ъ(х,у)(1хау

с

где ${х,у) — исходное изображение, а р(х,у) = — объединение кластеров всего

¿=1

изображения; — «раскрашенный» средними значениями кластер с номером г = 1,...,С ;

С — количество выделенных кластеров; Б = £>х£>у — область интегрирования, представляющая собой весь анализируемый квадрат. Причем «раскрашивание» кластеров производится средними радиояркостями (или их модальными значениями), взятыми из исходных данных для каждого кластера. Работу алгоритма кластерного анализа проиллюстрируем следующим практическим примером.

Кластерный анализ геосистемы Большого Васюганского болота по результатам спутниковых наблюдений

Уникальность и научная ценность Большого Васюганского болота заключаются не только в его огромных размерах, влияющих на климат Западной Сибири, но и в сложной ландшафтной структуре этой гидрогеоморфологической системы. Оно находится на стыке двух ботанико-географических подзон — южно-таежных и подтаежных лесов. Южная часть системы Большого Васюганского болота в значительной мере низинная, в то время как северная — почти исключительно верховая. Средняя Иртышско-Обская подпровинция кедрово-березовых лесов характеризуется преобладанием березовых лесов (72%). Древний рельеф территории — выположенная равнина. В целом болота представляют собой хорошо выраженные природные образования, имеющие широкое площадное распространение и характеризующиеся рядом специфических черт, местоположением, наличием определенных типов растительности и ее особым пространственным расположением, наличием открытых водных поверхностей и т.д. Эти признаки, наряду с труднодоступностью и сложностями в организации наземных исследований болотных территорий, делают весьма эффективным изучение болот с помощью различных дистанционных методов анализа аэрокосмической информации. Ввиду сказанного возникает необходимость структурирования экосистем Большого Васюганского болота с целью выявления однородных природных комплексов по данным дистанционных наблюдений с привлечением статистических методов тематического компьютерного дешифрирования и ГИС. Решение этой задачи, в частности, позволит экстраполировать данные ключевых участков, полученные экспедициями томских университетов и в лабораториях, на другие малодоступные участки болота с подобными текстурными характеристиками, выявленными в результате кластерного анализа.

ТЕХНИЧЕСКИЕ НАУКИ

12-----------

Птя иллюстрации работы программы были использованы снимки территории Большого

характер рисунка, его положенИю в рельефе, характеру питания,

ретным рир д поэтому количество дешифрируемых контуров полей на рас-

ВЕЕЕ

озера Пешкова, на рис. ¿ . ^„„„^ г,Тхтттт, ппа класса —

кластерного анализа, представленного 20 классами, на рис. 4 выДелень лишь два класса

„„„, „ ,ñr;,pr с пастительностью». Масштаб снимков и карты 1.20U UUU.

вепховья речных долин); 4) более дренируемые участки экотона, а также экосистемы по водораздельным „Дренируемые: 1) заболоченные леса (с "i«»«^^™^ верхних частях склонов речных долин; 2) смешанные леса в средних и нижних частях верхних час. ь няиболее сухие возвышенные местоположения (высокие 2 и 3

ГрТы,"ад«! первог^и второго порядка в реки третьего порядка с темнохвой-нГ"—; 4) антропогенные экосистемы - приноселковые пашни, вырубки и гари.

Заключение

Изучение структуры болотных экосистем имеет несколько прикладных аспектов. Болота могу™тьсяУ1дикаторами морфогенетических типов рельефа и месторождении полез-ГьГископаеГых? неотектонических движений, биогеохимических характеристик ландшафта гидрологических особенностей водосборов отдельных территорий и т.д. Использование 1я'нных дистанционного зондирования и тематического распознавания с применением кон-?е"Той инТормаГи позволило в первом приближении выявить структуру однородных природных комплексов, соответствующих группам типов ландшафтных местностей.

Литература

1. Дуда Р. Распознавание образов и анализ сцен / Р. Дуда, П. Харт ; пер. с англ. - М. :

МИР^11й6ва^наА. Классификация многомерных данных / С.А. Айвазян, З.И. Бежаева,

п R Гтяповеоов. - М. : Статистика, 1974. - 238 с.

3 пГкладная статистика: классификация и снижение размерности : справ, изд. / под

ред. С.А. Айвазяна. - М. : Финансы и статистика, 1989. 607 с.

4. Pao С.Р. Линейные статистические методы и их применение / С.Р. Pao , пер. с англ.

- М. : Наука, 1968. - 548 с.

Рис. 1 _ фрагмент Васюганского болота у озера Пешкова

рис_ 2 — Фрагмент Васюганского болота

Рис. 3 — Результат кластерного анализа

Рис. 4 — Результат кластерного анализа

5. Жиглявский A.A. Математическая теория глобального случайного поиска / A.A. Жиглявский. - JI. : Изд-во Ленингр. ун-та, 1985. - 296 с.

6. Хан Г. Статистические модели в инженерных задачах / Г. Хан, С. Шапиро ; пер. с

англ. - М. : Мир, 1969. - 388 с.

8. Цыпкин Я.З. Адаптация и обучение в автоматических системах / Я.З. Цыпкин. - М.

: Наука, 1968. - 400 с.

Барышева Елена Александровна

Студентка кафедры автоматизированных систем управления ТУСУРа Телефон: (3822) 41 90 25 (дом.) Эл. почта: krosovka@nm.ru

Протасов Александр Константинович

Студент факультета управления качеством Западно-Сибирского филиала Российского гос. университета инновационных технологий и предпринимательства Телефон: (3822) 49 11 11, доп. 12 64 Эл. почта: prot@iao.ru

Шестакова Марианна Сергеевна

Аспирантка отделения послевузовского профессионального образования при ТУСУРе и ИОА СО РАН Телефон: (3822) 49 11 11, доп. 12 64 Эл. почта: prot@iao.ru.

Е.А. Barisheva, А.К. Protasov, M.S. Shestakova

A Cluster analysis for determining of structure of land cover geosystems the using results of space monitoring

The algorithm of cluster analysis is developed, which is approximated to the classical statistical interpreting of decomposition of the probability distribution mixtures restored using mixed selections. The algorithm performance is illustrated with solution of an applied problem of structuring the space shots of the Vasugan Bog.

УДК 621.371 С.Г. Госенченко

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

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

Высотный профиль индекса преломления воздуха строится (восстанавливается) путем расчета индекса преломления ТУ на определенных высотах по известной формуле Дебая [1]

где е — упругость водяного пара на соответствующей высоте; Р — атмосферное давление;

{ — температура воздуха.

Исходными для расчета в случае стандартного зондирования являются давление на «нулевой» высоте, температура воздуха и относительная влажность и, полученная по показаниям «сухого» и «смоченного» термометров аэрологического зонда и пересчитываемая в упругость водяного пара по формуле [2]