3.4. Анализ изображений для других образов-режимов [1981 Гренандер У. - Лекции по теории образцов: анализ образов]

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

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

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

(iii) Мы едва затронули задачу определения числа режимов. Если мы допустим, что оно может расти, то естественно ожидать лучшей аппроксимации, хотя численная сложность в конечном счете лишит результат привлекательности. Следовательно, необходимо тем или иным способом ввести в анализ изображений стремление к "простым" сегментациям.

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

Рассмотрим процесс выполнения программы ЭВМ при страничной организации памяти. Оперативная память разделяется на блоки, или страницы, имеющие обычно одинаковые размеры. В процессе выполнения программы происходят обращения к различным страницам оперативной памяти. Если перенумеровать страницы с помощью номеров 1, 2, 3, ..., n_p, то возникает цепочка целых чисел - временной образ, если интерпретировать время как число выполненных команд. Если страница, к которой производится обращение, уже отсутствует в оперативной памяти, то ее следует туда загрузить. Если же, с другой стороны, места в памяти недостаточно, то сначала необходимо удалить какую-то другую страницу. Интенсивность такого страничного обмена может легко достичь недопустимого уровня, и следует принять меры, препятствующие этому (см. Деннинг (1970)).

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

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

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

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

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

Это предположение подкрепляется одним эмпирическим исследованием (см. статью Фрайбергера, Гренандера и Сампсона (1975)). Данные относятся к машине IBM 360/67 с операционной системой CMS. Размер страницы памяти -4 килобайт, общее число страниц-64. Программа-монитор собирает данные об обращениях и регистрирует число обращений к каждой странице в пределах временного окна, в нашем случае -это тысяча команд. Размер этого окна достаточно мал для того, чтобы можно было подробно изучить особенности обращений. На самом деле можно сделать его больше, так как практически масштаб времени скорее характеризуется числом команд порядка 10-20 тыс.

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

Наибольший интерес для нас представляли большие программы, и мы остановились на компиляциях с Фортрана и ПЛ/1, а также на нескольких трансляциях с языка ассемблера. Данные записывались также на ленту для последующего анализа.

На рис. 3.4.1 представлены данные, характеризующие обращения к страницам памяти при компиляции программы с Фортрана. Более полно представляет картину рис. 3.4.2, на котором каждой горизонтальной линии соответствует одно временное окно. По горизонтали представлены различные страницы, которым в строке соответствует символ "х", если к данной странице в пределах данного окна производилось обращение. В противном случае соответствующий участок распечатки остается пустым. Отметим, что на рисунке представлена лишь часть процесса компиляции, а окна 63-260 исключены ради экономии места.

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

Нетрудно установить, почему возникают эти режимы. Это является следствием конструкции компилятора с Фортрана, и режимы соответствуют шести этапам работы компилятора, а именно вызову процедуры, разбору, размещению, объединению, генерации кода и выводу (см. IBM System 360 Operating System, FORTRAN IV (G) Compiler, Program Logic Manual).

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

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

Это частный случай Случая 3.4.3, т. 1 (образы-стационарные режимы), и значения α представляют здесь различные распределения вероятностей на множестве страниц. Введем эмпирический аналог, 64-мерный вектор f⁽ⁱ⁾ = (f⁽ⁱ⁾₁, f⁽ⁱ⁾₂,..., f⁽ⁱ⁾₆₄) где

f⁽ⁱ⁾_j - относительная частота обращений к странице с номером f на протяжении окон i, i + 1, i + 2, ..., i + N - 1. Здесь N - малое натуральное число, роль которого заключается в уменьшении Влияния случайных вариаций.

Далее из эвристических соображений получаем, что если окна i+ 1, i + 2, i + 2N - 1 приходятся на один и тот же режим, то норма

должна быть малой. Если она велика, то следует предполагать, что точка переключения лежит на дискретном отрезке [i, 1 + 2N - 1].

В соответствии с этими замечаниями наша эвристическая процедура должна сводиться к выборочному вычислению векторов частоты обращения к страницам памяти на N последовательных окнах и отысканию точек времени, в которых среднее квадратическое расстояние D между векторами частот принимает большие значения. Допустим теперь, что такое достаточно большое значение D обнаружено между векторами частот, вычисленными в окнах i - N, i - 1 и i, i + N- 1 соответственно. Попытаемся затем максимизировать значение D, так как это должно позволить более точно определить положение точки, разделяющей стационарные режимы. Векторы частот корректируются 2N раз, так чтобы значение D можно было вычислить между векторами с центром в каждом из 2N окон i - N, ..i - 1, i, i+ 1, ..i + N - 1. Если максимальное значение D достаточно велико, то окно, которому соответствует этот максимум, отождествляется с точкой переключения. Другими словами, этот алгоритм предполагает, что если расстояние между векторами частот обращения, вычисленными в последовательных группах из N окон, принимает максимальное значение в окне j, то окна j - N, ..., j - 1 и j, ..., j + N - 1 соответствуют совершенно различным режимам, имеющим различные векторы, характеризующие частоту обращения к страницам памяти. Это, однако, в определенном смысле идеализация, поскольку в принципе между режимами имеется переходный период, включающий по крайней мере несколько окон. Поэтому абсолютно точно локализовать точку переключения режимов невозможно.

Описанный только что алгоритм весьма подробно представлен блок-схемой, приведенной на рис. 3.4.3. Алгоритм запрограммирован на Фортране и реализован с помощью главной программы "РЕЖИМ" и подпрограммы "ПОИСК", которая вызывается для определения точки переключения путем максимизации расстояния D.

Очевидно, что выбор значения N - числа окон, на которых вычисляются векторы частоты, - имеет решающее значение для успешной работы алгоритма. Эта величина должна быть достаточно большой для того, чтобы могла быть получена хорошая оценка истинного характеристического вектора частот. С другой стороны, она должна быть меньше длины (выраженной в количестве окон) самого короткого стационарного режима из тех, что будут рас-смотрены. Важен также и выбор значения расстояния D, которое считается "достаточно" большим. На самом деле таких значений два. Первое из них DMAX1 - это величина, значение которой должно быть превышено для того, чтобы найти точное положение точки переключения режимов в последовательности окон. Второе значение DMAX2 - это величина, значение которой должно быть превышено максимальным значением D при поиске по последовательности в диапазоне 2Т окон (в подпрограмме "ПОИСК") окна, отождествляемого с точкой переключения режимов.

Выбор параметров TV, DMAX1, DMAX2 зависит от характера конкретного исследуемого класса программ. Алгоритм успешно работал при выделении режимов обращений к страницам для компилятора Фортрана со следующими параметрами: N = 20, DMAX1 = 0,99 и DMAX2 = 1,5. Подобный выбор N означает, что никакой диапазон, меньший 20 окон, или 20000 команд, не может быть идентифицирован как стационарный режим. Очевидно, что более короткие временные отрезки не имеют существенного практического значения и меньшие значения N затруднят точную оценку истинного характеристического вектора частот.

В правой части машинной распечатки, приведенной на рис. 3.4.2, указано значение расстояния D между последовательными векторами частот в диапазонах, включающих по 20 окон.

Точки переключения режимов для двух фортранных программ, определенные с помощью описанного алгоритма, показаны пунктирными линиями в тех случаях, когда эти точки не совпадают с точками переключений, определенными из интуитивных соображений (последние обозначены сплошными линиями). Максимумы функции расстояния, соответствующие этим разбиениям, указаны над или под пунктирными линиями. Минимальное значение D для любого окна, идентифицированного как точка переключения, равно 1,6771. Учитывая существование переходных периодов на стыках режимов, можно сказать, что данное разделение режимов, полученное с помощью описанного алгоритма, является настолько точным, насколько можно было рассчитывать. Отметим также, что для случаев компилирования программ с ПЛ/1 наилучшие результаты были получены при значениях DMAX1 = 1,5 и DMAX2 = 2,5. Эти более высокие пороговые значения оказались необходимыми в связи с большими размерами обрабатываемых множеств и значениями D.

Подробнее с изложенным материалом можно познакомиться по отчету Сампсона (1974), где, в частности, содержится критика выбора структуры образа.

Описанная выше процедура носит слишком частный характер: она не основана на явных утверждениях относительно структуры идеального образа и не учитывает деформаций. Для того чтобы устранить эти недостатки, мы провели еще одно исследование, с методологической точки зрения более общее, применительно к образам, похожим на рассмотренные в разд. 3.3.3, но деформированным из-за наличия шумов. Точнее, это означает, что данные х = I^D получаются из идеального изображения I, а время I = 1,2, ..., Т - дискретное.

Правые концевые точки режимов b_i (i = 1, ..., n - 1) соответствуют разрывам режимов, l_i - длины режимов и μ_i - средние значения режимов. При непрерывном времени это соответствует случаю L = D и q = 1 из предыдущего раздела.

Рассматриваемый здесь механизм деформаций представляет собой деформации в результате воздействия аддитивного шума, т. е. I^D = (х₁, х₂, . .., х_Т), где

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

где l_i - длина режима R_i, на множестве REG (множестве разбиений). Естественно, если разбиение reg тоньше разбиения reg', то q_x(reg) ≤ q_x(reg') и, избегая чрезмерной вычислительной сложности, мы приходим к тривиальному разбиению, где R_i = {i}. Если потребовать, чтобы n ≤ n_max, то численное решение (3.4.4) при больших Т вызовет непреодолимые вычислительные затруднения.

Альтернативный подход зиключается в том, чтобы, используя в качестве исходных диффузные временные образы, применить следующую бейесовскую процедуру. Допустим, что длины режимов l_i, средние μ_i, и шум e_t порождаются независимо, точки переключения образуют дискретный процесс восстановления и l_i имеет распределение F, т. е.

Средние значения режимов независимы и имеют одинаковые нормальные распределения

Без потери общности можно допустить, что λ = 0.

Эту ситуацию можно считать разновидностью случая 3.4.3 т. 1 с той разницей, что появляется дополнительная структура, представляющая свойство восстановления в точках переключения последовательных режимов.

если считать, что t = 0 - одна из точек восстановления.

При заданном разбиении вектор x = (x₁ ..., x_Т)' имеет нормальное условное распределение с нулевыми математическими ожиданиями и ковариациями

Обозначим через E_t матрицу размера l×l, все элементы которой равны единице, и через I_t единичную матрицу размера l×l. Тогда часть ковариационной матрицы, относящуюся к режиму длины l, можно записать так:

Отметим также, что полная ковариационная матрица х разделена на блоки, причем все блоки, находящиеся в окрестности главной диагонали, имеют вид (3.4.9). Из этого следует, что обратная матрица полной ковариационной матрицы точно так же разделена на блоки, причем блок (3.4.9) заменяется блоком следующего вида:

в последнем можно убедиться непосредственно, имея в виду, что E²_l = lE_l. Обозначим определитель M_l через D(l). Вычислив его, получаем

Тогда для функции условного распределения имеем

где х(R_i)- часть данных, относящихся к режиму {x_t, t∈R_i}.

Объединив (3.4.7) и (3.4.12), получаем функцию правдоподобия

Максимизация этого правдоподобия эквивалентна максимизации следующей величины:

Последнее эквивалентно максимизации величины

где H(l_i) = F(l_i) - F(l_i - 1), i =1, ..., n - 1, Н(l_n) = 1 - F (l_n - 1) и

(R_i) - среднее по i-му режиму.

При таком подходе числовая сложность n вводится в целевую функцию, так что разбиение, обеспечивающее максимизацию, уже не является наиболее мелким. Численная максимизация (3.4.15) все еще может требовать нереальных вычислительных затрат.

Для того чтобы обойти это препятствие, попробуем обратиться к эвристическим алгоритмам, в которых используется адаптация. Рассмотрим некоторое множество В операторов О, преобразующих пару из обрабатываемого в данный момент изображения I и I^d = (х₁, ..., x_Т) в некоторое новое изображение I':

Каждый оператор должен быть устроен так (эвристически), чтобы O(I, х) являлось в некотором смысле "лучшим" описанием х, чем I. Множество В обычно конечное, хотя каждый оператор может быть параметризован с помощью некоторого множества управляющих параметров.

Процесс поиска заключается в применении к некоторой начальной структуре I₀ последовательности преобразований

причем выбирать I₀ можно, исходя из любой смутной идеи, которая у нас появится относительно истинной природы данного изображения, либо совершенно произвольно (например, простейшая структура режима).

Тогда алгоритм восстановления Α можно рассматривать просто как некоторую последовательность операторов из множества В, т. е.

Определить такой алгоритм можно, задавая целиком последовательность или с помощью некоторой функции "преемника" J.

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

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

Теперь рассмотрим каждый оператор по отдельности.

где

_i - среднее по режиму R. Такая оценка приемлема только в тех случаях, когда разбиение reg близко истинному разбиению, и, следовательно, ей нельзя пользоваться на ранних этапах процесса поиска. Кроме того, влияние неточной оценки reg на оценку (SD1)² не совсем очевидно. Рассмотрим другую оценку:

здесь N - истинное число режимов, a μ_i - математические ожидания режимов. Следовательно, SD2² дает отклоненную оценку σ²; отклонение положительно и увеличивается при росте числа режимов и увеличении различия между средними соседних режимов.

Это отклонение в оценке SD2² объясняется членами (x_t+1 - x_t)², где t - истинные точки переключения; имеются возможности компенсировать это отклонение. Так, например, члены (x_t+1 - x_t)², имеющие чрезмерную величину, можно дезактивировать.

Можно предложить различные способы дезактивации. Пусть

Можно исключать верхнюю α-часть y_t или те значения y_t которые превышают некоторый фиксированный уровень отсечения, определяемый имеющимися данными.

Воспользуемся двухступенчатой процедурой. В первую очередь вычисляется среднее

затем те y_t значения которых превышают величину RLEVEL×ȳ, исключаются, и среднее, вычисленное по оставленным значениям y_t, используется в качестве оценки σ², т. е.

где скобки [ ] обозначают индикаторную функцию.

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

Отметим, что σ² в нашей модели данных - глобальная постоянная и получение ее как совокупной оценки (на основе SD² в противоположность SDI², где используются локальные признаки) отличается тем дополнительным преимуществом, что последующая организация локальных поисков становится независимой от этой операции.

STEST. Для оценки каждого из обнаруженных к этому этапу концов режимов (b_i) мы пользуемся критерием разделения

где

(оценка σ) - это рассмотренное выше среднее квадратическое отклонение (SD), а x_i - среднее по режиму R_i длины l_i; точка переключения b_i объявляется статистически значимой, если S(R_i R_i+1) ≥ CS. Это соответствует простому двухвыборочному критерию разности средних. Для совокупностей, имеющих нормальное распределение, целесообразно выбирать величину CS, равную приблизительно 2.

COMBINE. Оператор STEST применяется ко всем точкам переключения одновременно, что позволяет получить список статистически незначительных точек переключения. Затем эти точки исключаются. Отметим, что некоторые из остающихся после этого шага точек переключения могут быть незначительными, но они сохраняются в качестве рабочей гипотезы до тех пор, пока не будет проведена соответствующая проверка. Вся операция в целом названа COMBINE.

HTEST. Для проверки однородности каждого режима мы пользуемся простым критерием однородности Н (R_i), основанным на дисперсии, подсчитываемой для каждого режима.

Если режим R_i однороден, то при условии нормальности

CV соответствует числу средних квадратических ошибок относительно σ² = E[V_i]. Для совокупностей, имеющих нормальное распределение, целесообразно выбирать величину CV, равную приблизительно 1,5.

Отметим, что в предположении нормальности величина V_i действительно имеет χ²-распределение с l_i - 1 степенями свободы при условии однородности. Можно было бы вычислять уровень отклонения χ² или запоминать соответствующую таблицу (игнорируя влияние

²), но это было бы трудоемким, так как обычно такой тест приходится повторять многократно с различными l_i.

DIVIDE. Для определения правильной точки разделения режима, когда использование оператора HTEST позволило обнаружить неоднородность режима, мы выбираем подходящее число пробных точек разделения (NTRYDIV), скажем равномерно расположенных, вычисляем для каждой из этих точек величину "разности" между соответствующей левой и правой частями режима и выбираем в качестве новой точки переключения ту точку, в которой эта разность максимальна. Если длина рассматриваемого режима меньше, чем NTRYDIV+ 1, то используется максимально возможное число пробных точек.

Затем два полученных таким образом режима подвергаются изучению. Если длина левого режима I_l ≥ L и длина правого режима I_r ≥ L, то получаем новую точку переключения. Если I_l ≥ L и I_r ≥ L или I_l < L и L_r <1, то соответствующая концевая точка переносится в новую точку разделения (поскольку неоднородность может объясняться тем обстоятельством, что истинная точка переключения, которую предположительно представляет концевая точка, на самом деле расположена где-то внутри режима). Если I_l < L и I_r < L, то обе концевые точки переносятся в новую точку разделения (случай, когда обе концевые точки кажутся ложными).

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

Оператор HTEST применяется ко всем режимам, что позволяет получить список неоднородных режимов. Каждый из них затем делится или его концевая точка (точки) переносится, как было описано выше, причем все действия выполняются последовательно слева направо. Вся операция в целом определена как DIVIDE. Из рассмотренных до сих пор эта операция самая сложная.

Рис. 3.4.4. Пример применения оператора DIVIDE

На рис. 3.4.4 приведен пример применения операции DIVIDE при наличии сильной неоднородности. Оператор HTEST устанавливает однородность режима 4 и неоднородность девяти остальных режимов. Последующие действия включают следующие девять шагов:

Локальная коррекция точки переключения заключается в рассмотрении точек, образующих окрестность корректируемой точки (расположенные слева и справа от нее точки АМАХ -точки, в которых максимум достигается без нарушения с обеих сторон условия минимальной длины режима), и выборе в качестве скорректированной точки переключения точки, обеспечивающей максимальное разделение.

Коррекции выполняются во всех точках переключения последовательно слева направо; в целом весь процесс определен как ADJUST.

Чтобы ускорить процесс поиска оператор GUESS обращается к предварительному списку точек переключения (если они вообще есть), подготовленному оператором SDEVI, очищает его (слева направо) так, чтобы оказались исключенными режимы длины, меньшей L, и затем делит оставшиеся режимы на пробные режимы длины L каждый (если деление возможно, в противном случае в процедуру вносятся небольшие изменения). После этого вызывается оператор COMBINE для изучения полученной черновой структуры.

Можно применить и более тщательно проработанные гипотезы. Введем, например, для t = h, h + 1, ...,T - h локальные средние

Предполагается, что d_t,h принимает большие значения, когда t - истинная точка переключения. Таким образом,

"правдоподобное" множество кандидатов в приближенные точки переключения. Поскольку вероятно, что D_h состоит из связных кластеров, следует выделить эти кластеры и составить предварительный список точек переключения (скажем, всех точек, имеющих в соответствующем кластере максимальное значение d_t,h).

Машинные эксперименты показали, что при удачном выборе h и c_h данный метод работает вполне удовлетворительно при выделении ярко выраженных точек переключения. Мы, однако, не пользуемся им в наших дальнейших изысканиях, так как сознаем, что даже при очень "грубом" операторе GUESS, описанном выше, остальные адаптивные операторы, если применять их должным образом, обеспечат достаточно хорошие результаты.

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

При построении алгоритма эти операторы можно объединять различными способами. Целесообразным представляется следующий. В качестве исходной выбирается простая структура "1 режим" и к ней применяется оператор SDEVI, что обеспечивает оценку σ², необходимую для применения остальных операторов и, возможно, указывающую также на несколько точек переключения. Вторым, естественно, применяется оператор GUESS, а затем предварительное множество точек переключения корректируется. На этом предварительный этап обработки заканчивается. Дальнейшие усовершенствования заключаются в отыскании новых точек переключения или удалении явно ложных точек переключения с помощью оператора DIVIDE, коррекции новой структуры и последующей проверке точек переключения с помощью оператора COMBINE. Этот алгоритм можно записать как

здесь использованы следующие обозначения S = SDEVI, G = GUESS, А = ADJUST, D = DIVIDE и C = COMBINE.

Итак, функция "преемник J" задается как J: В² → В. Реально встречаются следующие сочетания операторов:

Если рассматривать три последовательно применяемых оператора DIVIDE-ADJUST-COMBINE как один, то наш алгоритм можно представить в виде Α = PII..., где P = PRELIM

(предварительная обработка), а I = DAC = IMPROVE (улучшение).

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

При составлении программы структура режима представляется матрицей, состоящей из трех строк и столбцов, каждый из которых характеризует некоторый режим. Первая строка содержит точки переключения [b₁, b₂, ..., b_n], вторая - сумму x_t по соответствующим режимам и третья-сумму квадратов x_t. На основе этой матрицы нетрудно вычислить критерий разделения S (/?,-, Ri + l) и критерий однородности И (/?,); функционирование алгоритма определяется заданием шести описанных выше параметров L, RLEVEL, CS, CV, NTRYDIV и AMAX.

(i) Точки переключения порождаются дискретным равновесным процессом восстановления, причем длины последовательных режимов определяются выражением l = 1 + B(s, p) , где B - бинормальная случайная величина.

(ii) Средние режимов порождаются независимо из распределения N(0, τ²).

(iii) Шумы порождаются независимо из распределения N (0, σ²).

Распечатка, приведенная на рис. 3.4.5, иллюстрирует процесс поиска на данных длины 200 со следующими параметрами: s = 72, р = 1/3 (математическое ожидание длины режима - 25), τ = 2 и σ = 1. Алгоритм имел следующие параметры: L = 10, NTRYDIV = 10, АМАХ = 5, CS = 1, 8, CV = 1, 5 и RLEVEL = 6, 6. Для контроля хода процесса восстановления величина средней остаточной дисперсии

определялась после каждого применения оператора.

В распечатке истинная структура режима приведена под словом ДАННЫЕ: в первом столбце содержатся точки переключения, а во втором -- соответствующие средние по режиму.

При поиске эвристические адаптивные операторы объединяются в два блока. Первый PRELIM = SDEVI - GUESS - ADJUST применяется в начале и только один раз; второй IMPROVE = DIVIDE - ADJUST - COMBINE применяется после блока PRELIM многократно.

На рис. 3.4.6 приведен график изменения величины средней остаточной дисперсии во времени в процессе этого поиска, а на рис. 3.4.7 представлены идеальное, деформированное и окончательное восстановленное изображения.

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

Рис. 3.4.6. Поведение средней остаточной дисперсии в процессе поиска

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

Рис. 3.4.7. Идеальные (*) деформированные (°) и восстановленные (⋅) изображения. Там где * на рис. отсутствует, имеет место * = ⋅ с точностью аппроксимации, использованной при воспроизведении

Практическая проверка описанного эвристического алгоритма проводилась на нескольких моделях, сходных с моделью 1, но отличающихся от нее некоторыми деталями.

Модель данных 2. Шум, порожденный независимо источником с равномерным распределением.

Модель данных 3. Шум, порожденный независимо источником с двумерным экспоненциальным распределением.

Модель данных 4. Средние по режиму образуют стационарный процесс в узком смысле, причем любые k последовательных средних μ_i,...,μ_i+k-1 совместно нормально распределены с нулевыми средними и трех диагональной ковариационной матрицей

Модель данных 5. Распределение длин режимов определяется смесью двух биномиальных распределений:

параметры s₁, р₁, s₂ и р₂ можно выбрать так, чтобы получить смесь длинных и коротких режимов.

Для каждой модели данных было сформировано по десять выборок, каждая длиной 200, однако в реальных условиях следует ожидать появления больших чисел. Во всех случаях соответствующие параметры выбирались так, чтобы математическое ожидание длины режима было равно 25, среднее квадратическое отклонение средних по режиму τ = 2 и среднее квадратическое отклонение шума σ = 1.

Для всех пяти моделей данных были приняты L = 10, NTRYDIV = 10, АМАХ = 5 и CS = 1, 8, а значения CV и RLEVEL выбирались, исходя из характера распределения шума. Для шума, распределенного нормально, CV = 1,5. Для равномерно распределенного шума при однородном режиме

то для равномерно распределенного шума мы выбрали CV = 1,5 √2/5 ≈ 1

При однородном режиме и шуме с двумерным экспоненциальным распределением

Выбор величины RLEVEL основывается на распределении

Таким образом, для нормально распределенного шума RLEVEL = 6,6, для равномерно распределенного шума RLEVEL = 4,86, для шума с двумерным экспоненциальным распределением RLEVEL=9.

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

где m_t-математическое ожидание, a

_t - оценка среднего значения.

Может оказаться полезным и критерий, зависящий лишь от истинного разбиения и его оценки. Пусть b₁, b₂, ..., b_N-1 - истинные точки переключения

₁,...,

_-1

- оценки точек переключения. Мы будем говорить, что

_i оценивает b_j, если

_i - ближайшая к b_j вычисленная точка, a b_j - ближайшая к

_i истинная

Пусть N_a - общее число таких ассоциированных пар. Введем

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

Пусть N_d - число различных истинных точек переключения в N_a ассоциированных парах; тогда число пропущенных точек переключения равно N - 1 - N_d. Введем

Величина D₂измеряет долю пропущенных истинных точек переключения.

Пусть

_d - число различных вычисленных точек переключения в N_a ассоциированных парах; тогда число избыточных точек переключения равно

- 1 -

_d. Введем

Величина D₃ измеряет долю избыточных точек переключения. Результаты эксперимента сведены в табл. 3.4.1. В таблице СТ

Таблица 3.4.1. Сводка результатов экспериментов

обозначает время вычислений в секундах, затраченное до завершения поиска, MRV - среднюю остаточную дисперсию и SD - оценку среднего квадратического отклонения шума. Строки Ave и s. d. содержат выборочные средние и средние квадратические отклонения соответственно, вычисленные по десяти выборкам.

В общем алгоритм хорошо работал на проверочных данных. Усредненная средняя квадратическая ошибка составила около 0,16 по сравнению с собственной дисперсией шума 1. Около 24 процентов истинных точек переключения оказались пропущенными, но по большей части они имели несущественный характер - небольшие различия средних у соседних режимов. Точность обнаруженных точек переключения составила около 2,3%, и, грубо говоря, на 200 точек данных приходилась одна избыточная точка переключения. Время вычислений составило около 6,4 с, что при работе с языком АПЛ приемлемо. Для практического использования этот алгоритм следует запрограммировать на более подходящем языке для того, чтобы уменьшить время работы центрального процессора. Окончательное значение средней остаточной дисперсии и оценка дисперсии шума очень близки к истинным величинам, как и следовало ожидать.

Из 50 рассмотренных выборок только в трех случаях поиск заканчивался циклами размера больше единицы (в двух случаях длина равна двум и в одном -трем).

Здесь следует заметить, что рассмотренные нами модели порождают стационарные случайные процессы {e_t}, {m_t} и {x_t}. Работая со спектральными плотностями, можно вывести не зависящий от времени линейный фильтр L который, будучи примененным к х, восстанавливает m наилучшим возможным в смысле наименьших квадратов способом. Таким образом, можно найти L минимизирующий L₂-ошибку, и вывести выражение для минимума ошибки при использовании линейно-квадратических методов. Для этого, естественно, необходимо знать параметры модели данных, и, таким образом, возникает новая задача, отличающаяся от той, которую мы рассматривали. Эту задачу мы изучаем в данном контексте, потому что она дает нам некий стандарт работоспособности, несмотря на то, что указанный метод требует наличия подробной априорной информации об образах в отличие от описанного выше.

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

(а) Точки переключения b_i образуют дискретный процесс восстановления следующего вида:

Будем считать, что данный процесс начался достаточно давно и потому можно пользоваться соответствующими асимптотическими результатами, или, что то же самое, рассматривать соответствующий равновесный процесс восстановления

(б) Средние режима _i образуют белый шум с нулевым математическим ожиданием и дисперсией τ².

(в) Шумы e_t образуют белый шум с нулевым математическим ожиданием и дисперсией σ₂.

Как и выше, m_t = μ_i, если t∈R_i и x_i = m_t + e_t Располагая выборкой x_t, t = 1, ..., T мы хотим воспользоваться теорией линейной фильтрации для восстановления m_t, t = 1,..., T. Пусть reg = реализация процесса переключения, тогда

На основании свойств равновесного процесса восстановления получаем, что

Спектральная плотность процесса m_t определяется как

Обозначив ковариационные функции процессов x_t, m_t и е_t через r_x(h), r_m(h) и r_e(h) соответственно, эти процессы можно записать в крамеровском представлении

a F_x, F_m и F_e - соответствующие спектральные плотности.

Чтобы построить не зависящий от времени линейный фильтр L который в применении к х восстанавливает m наилучшим возможным в смысле наименьших квадратов способом, рассмотрим

где γ - L₂(F_x) - функция с комплексными значениями.

Здесь dz_x = dz_m+ dz_e, dz_e ⊥ dz_m представлены в обычной символической записи.

Предполагая, что спектральные плотности равны f_m и f_e соответственно, получаем

Чтобы проиллюстрировать, как применяется (3.4.62), конкретизируем рассматриваемый процесс восстановления, потребовав, чтобы

длина режима I имела биномиальное распределение,

Тогда путем непосредственных вычислений получаем, что

Для того чтобы найти коэффициенты фильтра, представим функцию γ₀ ее разложением в ряд Фурье

где α = (τ/σ) - отношение сигнал/шум. Поскольку γ₀(λ) = γ₀(-λ), γ₀(λ + 2π) = γ₀(π) и целесообразно выбирать симметрический фильтр, то пусть

В качестве численной иллюстрации примем, что l_m = 1, s = 50, р = 0,4, μ = 21, α = 2, и вычислим υ(n), n = 1, ..., l_m + s- 1.

На рис. 3.4.8. приведены графики х, m и m^* - восстановленного сигнала, который как будто хорошо оценивает m, так что естественно воспользоваться критериями вида

и в тех случаях, когда D²_t,s ≥ const, мы предполагаем наличие точки переключения на отрезке [t, t + s].

Чтобы сопоставить качество работы эвристического адаптивного алгоритма и метода фильтрации, мы выбираем модель и параметры алгоритма в соответствии с центральным планом для 7 переменных с 22 центральными точками и 14 точками, расположенными на координатных осях; общее число наблюдений -100. Недостаток места не позволяет нам здесь вдаваться в подробности этого довольно сложного вычислительного эксперимента; читатель, интересующийся этими вещами, может обратиться к отчету Анга, (1974), в котором обсуждается примененная нами методика поверхности отклика. Преимущество эвристического адаптивного метода заключается в том, что он требует лишь небольших априорных сведений о структуре области и его реализация связана с умеренными затратами машинного времени. Ограничения метода фильтрации связаны с необходимостью иметь большее количество априорной информации, а также с тем обстоятельством, что этот метод использует только линейные операции в чисто нелинейной ситуации.

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

Рис. 3.4.8. Выборка идеальных (*) - m_t, искаженных (°) - х^* и восстановленных (⋅) - m_t сигналов

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

Первый аналогичен алгебре изображений, рассмотренной в разд. 3.3 (за исключением дискретности времени), но переход от одного режима к другому уже определяется процессом Бернулли; сравните также с допущением (3.4.63). Мы по-прежнему оперируем одним L-оператором, одним индексом класса образующих, и в каждой дискретной точке времени вероятность начала нового режима определяется некоторой величиной р. Под новым режимом в данном случае понимается изменение функции времени посредством подачи существенного импульса, возмущающего уравнение L - х = y + n. Символ у представляет здесь импульсы, большая часть которых равна нулю при малых р, а n - гауссов белый шум. Тогда справедлив следующий результат.

Теорема 3.4.1.Пусть образ-режим x_t порождается как решение уравнения

где n_t - гауссов шум N (0, σ²) и y_t независимо тождественно распределены, причем y_t = 0 с вероятностью р и y_t = N (А, τ²) с вероятностью 1 - р. Восстановление режимов по методу правдоподобия достигается выбором в качестве точек переключения всех точек времени t^* таких, что

Замечание. Получающийся алгоритм очень прост. Нужно просто вычислить (Lx_t)² при t = 1, 2, ..., n и найти те значения t^*, для которых справедливо неравенство (3.4.76). Число таких точек указывает нам (ориентировочно) на число режимов, и, если это необходимо, можно оценить также размер импульсов (см. приведенное ниже уравнение (3.4.78)).

Доказательство. Рассмотрим изображение, порожденное неизвестными режимами, на интервале (t_υ, t_υ+1) при 1 ≤ t_{υ ≤ n} для υ = 1, 2,..., l, t₀ = 1, t_l+1 = n. В каждой точке t_υ разностное уравнение (3.4.75) подвергается некоторым воздействиям y = N (А, τ²). Тогда функция совместного распределения для x₁, x₂,...,x_n, l, t₁, t₂,...,t_l, A₁, A₂,....,A_l принимает вид

здесь T₁ = {t₁, t₂,..., t_l}, a T₂ образуется остальными целыми числами от 1 до n. Оценки импульсов A_j по методу правдоподобия определяются с помощью следующего выражения:

Последнее выражение можно переписать в следующем виде:

Теперь остается лишь варьировать Т1У оно может быть любым подмножеством множества {1, 2, 3, ..., n}. Цель этих вариаций - увеличить L^*, или, что эквивалентно, сделать сумму в уравнении (3.4.81) по возможности минимальной. Очевидно, это можно обеспечить, проводя суммирование исключительно по отрицательным членам, и поэтому необходимо выбирать множество следующего вида:

На этом доказательство теоремы заканчивается.

Примечательно, насколько просто оказывается множество Т^*₁ при записи в такой форме с использованием оператора L. Теперь попытаемся проанализировать случай, когда число индексов класса образующих превышает единицу, так что имеется множество разностных операторов {L^α} и правила, указывающие, как осуществляются переходы α → α'(см. случай 3.4.4, т. 1, образы операторного режима), и модифицированные с учетом дискретности времени. Мы полагаем, что во всех точках времени переходы между режимами определяются как

Эти вероятности, естественно, могут выродиться в нуль или единицу-детерминированные переходы. Будем считать, что в пределах некоторого режима, имеющего индекс класса образующих α, L^αx равно белому шуму с параметрами N(0, σ²).

Следовательно, при порождении этих изображений отражается не только влияние белого шума, но также и влияние перескоков от одного индекса класса образующих к другому. Это показано схематически на примере, приведенном на рис. 3.4.9, где все стрелки, за исключением L² → L⁴ представляют нормальное поведение. С вероятностью р₂₄ система проскакивает L₃. Значения p_ij i = 1, 2, 3, 4, должны быть большими по сравнению с остальными с тем, чтобы обеспечивалось порождение режимов типа временных образов. В противном случае "режимы" будут столь коротки, что с ними можно уже не считаться.

Чтобы обеспечить бейесовское восстановление изображения, вычисляем совместную вероятность m

В уравнении (3.4.84) идеальное изображение I описывается с помощью α(t), t = 1, 2,...,Т, где α(T) постоянна на интервалах, соответствующих n режимам 1 = t₁ ≤ t < t₂, t₂ ≤ t < t₃,...., t_n ≤ t < T + 1 = t_n+1. Как и прежде, I^D = (х₁, х₂, ..., х_Т). Решая задачи такого типа, мы сталкиваемся с обычными затруднениями - для вычисления (3.4.84) необходимо знать α(0), а также ряд значений х₀, х_-1, х_-2>, ... в зависимости от максимального порядка операторов L. При больших Т знать некоторые начальные значения х уже не так важно; это означает, что до начала анализа какая-то малая часть деформированного изображения I^D исключается. Значение α(0) будет выбрано позже.

Перепишем уравнение (3.4.84) так, чтобы в более явном виде представить структуру режимов:

здесь l_υ - длина t_υ+1 - t_υ υ-го режима идеального изображения, β_υ - индекс класса образующих υ-го режима и

при q_ββ' = - ln p_ββ'. Уравнение (3.4.87) точно показывает вклад в ln l от каждого режима.

Чтобы максимизировать значение l(I^D, I) по I, необходимо отыскать некоторую функцию α(t), удовлетворяющую условию

Теперь задача приведена к форме, удобной для применения динамического программирования. Действительно, введем

Эта сумма принимает минимальное значение, когда функции α(⋅)

удовлетворяют условию α(s) = i. Тогда, воспользовавшись выражением (3.4.88), можно получить следующую рекурсию:

Это позволяет нам с помощью итерации дойти назад до значения s = 1 и получить последовательность α(⋅) являющуюся решением уравнения (3.4.88).

Теорема 3.4.2.Восстановление по методу правдоподобия для модели, предусматривающей перескоки с режима на режим, можно обеспечить при помощи последовательного решения методом динамического программирования уравнений (3.4.91) и (3.4.89).

Этот алгоритм еще не проверялся на числовых данных.