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

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

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

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

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

С точки же зрения подготовки специалистов, основой закрепления теоретического курса «Проектирования РЭС» является лабораторный практикум. Поэтому разработка лабораторного программного комплекса для исследования эффективности пространственно-временной обработки изображений является важной и актуальной задачей.

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

1. Технико-экономическое обоснование

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

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

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

Применение персональных ЭВМ (электронных вычислительных машин) в лабораторном практикуме для исследования систем автоматического управления дает возможность выполнить:

анализ функционирования системы в различных режимах с учетом разброса параметров элементов и наличия дестабилизирующих факторов;

синтез структуры и принципиальной схемы устройств;

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

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

2. Теоретическая часть

2.1 Задача пространственно-временной обработки изображений при наличии шумов и помех

Пространственная модель изображения

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

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

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

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

Таким образом, линейная модель наблюдения изображения в условиях помех принимает вид [25]:

где — аддитивный шум.

Виды фильтров пространственной обработки изображений

В общем случае схему искажения и восстановления (фильтрации) изображений можно представить в виде, изображенном на рисунке 2.1.

Рисунок 2.1 — Схема искажения и восстановления изображений

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

Тогда ошибка оценивания в каждой точке изображения определяется выражением:

а средняя квадратичная ошибка (СКО) — через ее квадрат, то есть дисперсию ошибки:

.

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

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

Еще одним широко распространенным на практике критерием эффективности обработки изображений является связанный с рассмотренным критерий максимума пикового отношения сигнал-шум (ПОСШ) [25]:

.

Фильтры, основанные на порядковых статистиках

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

.

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

Необходимо отметить, что дисперсия медианы белого шума на выходе медианного фильтра определяется выражением [24]:

где — размер апертуры.

Равенство (2.1−6) показывает также, что для нормального белого шума дисперсия медианы приблизительно на больше, чем дисперсия для среднего. Следовательно, скользящее усреднение подавляет нормальный белый шум несколько лучше, чем медианный фильтр с такой же апертурой. Иначе говоря, чтобы медианный фильтр давал ту же дисперсию шума, что и скользящее усреднение в его апертуре должно быть на 57% больше точек.

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

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

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

Хотя медианный фильтр значительно более распространен в обработке изображений, чем остальные виды фильтров, основанные на порядковых статистиках, тем не менее, он не является единственным. Медиана представляет собой 50-й процентиль упорядоченного набора чисел, но, как следует из основ статистики, упорядочивание предоставляет много других возможностей. Например, использование 100-го процентиля приводит к так называемому фильтру максимума, который полезен при поиске на изображении наиболее ярких точек по отношению к окружающему фону. Процентиль 0 является фильтром минимума, используемым для поиска противоположных значений.

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

где — оператор псевдомедианной фильтрации, действующий следующим образом. Пусть — последовательность элементов, тогда псевдомедиана такой последовательности определяется так:

.

Медианные фильтры хорошо работают до тех пор, пока пространственная плотность импульсного шума невелика (эмпирическое правило — и не превышают 0.2). Модифицированная медианная фильтрация помогает справиться с импульсным шумом, вероятности которого превышают указанные значения. Дополнительное преимущество такого модифицированного медианного фильтра состоит в том, что такой фильтр «старается сохранить детали» в областях, искаженных не импульсным шумом. Обычный медианный фильтр таким свойством не обладает.

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

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

Ветвь А: если, перейти к ветви Б;

иначе результат равен .

Ветвь Б: если, результат равен ;

иначе результат равен ,

где — минимальное значение яркости в , — максимальное значение яркости в , — медиана значений яркости в , — значение яркости в точке .

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

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

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

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

· Фильтры, основанные на выборе максимального и минимального

значений.

Хотя медианные фильтры, безусловно, принадлежат к числу наиболее используемых в обработке изображений фильтров, основанных на порядковых статистиках, это отнюдь не единственный пример таких фильтров. Медиана представляет собой 50-ый процентиль упорядоченного набора чисел, однако использование иных статистических характеристик представляет много других возможностей. Например, использование 0-го процентиля приводит к фильтру, основанному на выборе минимального значения (минимальному фильтру), который задается выражением:

.

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

Использование же 100-го процентиля приводит к фильтру, основанному на выборе максимального значения (максимальному фильтру):

.

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

Весовая пространственная фильтрация

Пространственные методы фильтрации, рассмотренные ранее, в общем случае описываются выражением вида [22]:

где — искаженное изображение на входе фильтра, — восстановленное изображение на выходе фильтра, а — оператор фильтрации, определенный в некоторой окрестности, точки.

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

С целью повышения эффективности фильтрации, то есть для уменьшения квадрата СКО вида (2.1−3), используется весовая пространственная фильтрация.

Пространственная фильтрация изображения с помощью весового фильтра задается выражением общего вида:

где — весовая функция или функция окна, однозначно определенная в окрестности .

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

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

Временная фильтрация

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

Рисунок 2.2 — Схема временной фильтрации изображений

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

.

С учетом этого соответствующие выражения для фильтров, определяемых соотношениями (2.1−5) — (2.1−10), при временной фильтрации принимают следующий вид:

— среднеарифметический фильтр: ;

— среднегеометрический фильтр: ;

— среднегармонический фильтр: ;

— среднеконтргармонический фильтр:;

— медианный фильтр: ;

— псевдомедианный фильтр: ;

— минимальный фильтр: ;

— максимальный фильтр: .

Необходимо отметить, что, как и в случае пространственной фильтрации, для повышения эффективности временной фильтрации, то есть для уменьшения квадрата СКО вида (2.1−3), используется весовая временная фильтрация, которая задается выражением вида:

где — весовая функция, которая вместе с апертурой формирует маску.

Весовая временная фильтрация осуществляется таким же способом, что и пространственная.

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

2.2 Методы оптимизации при пространственно-временной обработке изображений

Задача весовой пространственно-временной фильтрации как нелинейная задача оптимизации

В общем случае при весовой пространственно-временной фильтрации выражение квадрата СКО ошибки приближения восстановленного и исходного изображений принимает вид:

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

Таким образом, при весовой пространственно-временной фильтрации изображений возникает задача нахождения такого оптимального вектора весовых коэффициентов, который обеспечивал бы минимум квадрата СКО (2.2−1), то есть минимум нелинейной функции :

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

Двухэтапный алгоритм весовой фильтрации

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

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

Необходимо отметить еще одну трудность, сопутствующую задаче вида (2.2−2). Имеется в виду разнородность компонент вектора весовых коэффициентов. Это может привести к появлению многомерных «оврагов» у целевой функции, что препятствует высокой скорости сходимости методов спуска при оптимизации в совокупном пространстве. Сходимость процесса повышается при циклическом переборе спусков в пространстве однородных весовых коэффициентов.

В связи с особенностями целевой функции (многоэкстремальность, не выпуклость, неаналитичность, плохая обусловленность) при решении задачи (2.2−2) необходимо использовать алгоритм, состоящий из двух этапов — глобального и локального [19,20]. На первом этапе с использованием одного из методов глобальной оптимизации находится грубое приближение оптимального вектора весовых коэффициентов (область глобального экстремума). А на втором этапе компоненты вектора определяются с заданной точностью одним из локальных методов. Очевидно, что время поиска и точность решения целиком определяются эффективностью используемых алгоритмов оптимизации.

Глобальный случайный поиск

Исторически первым и довольно долгое время единственным широко используемым методом поиска глобального экстремума являлся метод, названный впоследствии «мультистарт». Этот метод состоит в многократном (последовательном или параллельном) отыскании локальных минимумов из различных начальных точек, чаще всего случайных (такой вид алгоритма носит название «случайный мультистарт»).

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

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

Выбирается некоторое распределение вероятностей на W и полагается .

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

Согласно определенному правилу конструируется распределение вероятностей на W.

Проверяется условие останова, и если оно не выполняется, то осуществляется переход к шагу 2 с заменой .

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

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

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

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

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

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

Классификация методов глобальной оптимизации

Использование двухэтапного алгоритма весовой фильтрации предполагает на первом этапе применение процедуры, обеспечивающей поиск точки (области) глобального экстремума —. Математическая формулировка этого факта имеет следующую форму [8]:

при всех .

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

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

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

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

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

Алгоритмы глобальной оптимизации в настоящее время принято разделять на два класса [11,13,17]: использующие модель глобального поведения функции (рациональные) и не использующие таковой (эвристические).

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

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

Основные идеи построения рациональных и эвристических алгоритмов определяют следующую классификацию [11,17]:

Методы поиска глобального экстремума, основанные на использовании алгоритмов локальной оптимизации.

Алгоритмы минимизации функций с ограниченной скоростью изменения.

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

Глобальный случайный поиск.

Методы декомпозиции.

Анализ большого количества работ [8,11,17], посвященных методам глобальной оптимизации показывает, что с точки зрения объема вычислительных затрат, точности нахождения глобального экстремума и возможности программной реализации наиболее предпочтительными являются алгоритмы глобального случайного поиска и методы декомпозиции.

Метод динамического программирования

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

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

k-м шаге решения для всех и фиксированных значений параметрического вектора (вектора, полученного на предыдущем этапе) рассчитывается функция:

где ,…, — так называемые функции оптимального соответствия, построенные на () — м шаге.

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

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

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

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

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

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

применение интерполяции с использованием полинома Лагранжа при нахождении глобального экстремума функции оптимального соответствия;

использование процедуры адаптивного сдвига сеточных узлов на каждой итерации поиска.

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

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

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

Классификация методов локальной оптимизации

Методы поиска локального экстремума принято ранжировать следующим образом [12−18]:

1. Метод Ньютона и его модификации.

2. Квазиньютоновские методы.

3. Градиентные методы.

4. Методы прямого поиска.

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

Большинство методов минимизации являются итерационными по своей природе. Это означает, что, начиная с заданной точки, называемой начальным приближением, в дальнейшем получается последовательность точек, которая, в принципе, должна сходится к точке минимума целевой функции. Вычисление называется k-й итерацией, а точка — k-м приближением. В результате основное уравнение для k-й итерации записывается в виде:

,

где — скалярная величина (длина шага), обеспечивающая неравенство , — направление поиска, правила вычисления которого «определяют лицо алгоритма» минимизации и его эффективность.

Так, в методе Ньютона направление поиска определяется системой линейных алгебраических уравнений:

где — матрица Гессе, матрица вторых частных производных, а — вектор градиента функции .

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

В квазиньютоновских методах направление поиска вычисляется по формуле [14,17]:

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

где — поправочная матрица различных рангов. Существуют различные формулы пересчета приближений обращенной матрицы Гессе, использующие поправки разных рангов [14, 17]:

Симметричная формула ранга один (SR1-формула):

где, а — соответствующее приращение градиента.

Формула Бройдена-Флетчера-Гольдфарба-Шанно (BFGS-формула):

где .

Формула Дэвидона-Флетчера-Пауэлла (DFP-формула):

.

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

которое решается с использованием разложения матрицы по методу Холесского .

Каждой квазиньютоновской формуле для пересчета оценки обратной матрицы Гессе, соответствует (в том смысле, что из следует, где — единичная матрица) преобразованная формула пересчета матрицы [14, 15, 18]:

Преобразованная симметричная формула ранга один (SR1-формула):

где .

Преобразованная BFGS-формула:

.

Преобразованная DFP-формула:

где .

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

Имея разложение Холесского матриц, можно получать оценки их чисел обусловленности. Так справедливо неравенство, где — максимальный и минимальный диагональные элементы фактора. На промежуточных операциях такие оценки позволяют выявить ситуации, когда матрица обусловлена так плохо, что вычисляемый вектор из-за ошибок округления не будет направлением спуска. При этом факторы разложения можно подправить так, что обусловленность реализующей матрицы станет приемлемой [14, 15].

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

Используя в квазиньютоновском поиске факторизацию Холесского, легко обеспечить положительную определенность матрицы [14, 17]. Для этого на очередном шаге вычислений коэффициентов множителей Холесского соответствующий диагональный элемент матрицы увеличивается настолько, чтобы окончательное произведение было положительно определенным и хорошо обусловленным. Когда этот элемент выбран, очередной столбец множителя однозначно определяется требованием совпадения матриц и. Таким образом, в результате получается представление:

где — неотрицательная диагональная матрица, равная нулю, если матрица существенно положительно определена.

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

ожно показать [12, 14], что направление для очередной итерации в методе сопряженных градиентов при минимизации квадратичных функций вычисляется как линейная комбинация градиента в текущей точке и предыдущего направления :

где вычисляются по следующим теоретически эквивалентным формулам [12]:

Формула Флетчера-Ривса: .

Формула Полака-Рибьера: .

Теоретически этот метод конечен [12]: при точных вычислениях он должен найти решение не более чем за n шагов.

Метод сопряженных градиентов легко обобщается на задачи минимизации нелинейных функций общего вида. Для этого надо заменить используемую в нем формулу расчета какой-нибудь процедурой одномерного поиска и уточнить, будет ли направление всегда вычисляться по формуле (2.2−17) или допускаются периодические отступления. Последние принято называть рестартами или восстановлениями [11, 14]. Можно, например, принимая во внимание n-шаговую теоретическую сходимость в квадратичном случае, через каждые n шагов отказываться от направления (15) и брать в качестве очередного антиградиент. Соответствующий метод называют традиционным.

Традиционный метод сопряженных градиентов сходится в тех же предположениях, что и метод наискорейшего спуска [11,15]. Для достаточно широкого класса функций теоретическая скорость его сходимости при точной одномерной минимизации оказывается n-шаговой сверхлинейной. Однако, как показывает практический опыт, из-за ошибок округления реальная скорость сходимости метода сопряженных градиентов почти всегда линейна независимо от того, применяются рестарты или нет. Поэтому применение метода сопряженных градиентов следует считать успешным, когда удовлетворительное решение удалось получить менее чем за 5n шагов.

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

Конечно-разностная аппроксимация производных

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

Простейшая конечно-разностная аппроксимация градиента получается на основе разложения функции в ряд Тейлора [8,11]:

где — интервал конечной разности, связанный с j-й компонентой аргумента; - единичный вектор, j-я компонента которого равна единице; - ошибка порядка, связанная с отбрасываемыми членами ряда Тейлора и называемая ошибкой отбрасывания. Выражение (2.2−20) называется односторонней конечно-разностной аппроксимацией.

Если в формуле (2.2−20) выполняется еще один шаг тейлоровского разложения, то центральная конечно-разностная аппроксимация градиента имеет вид:

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

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

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

Неточность конечно-разностных аппроксимаций вида (2.2−20) — (2.2−22) связана с ошибками трех типов [8,11,15]:

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

Ошибка потери значащих цифр, связанная с разностью двух близких значений функции. Если — относительная точность вычисления (в лучшем случае, где m — число двоичных цифр, с которыми ведутся вычисления). Абсолютная ошибка вычисления имеет величину, тогда максимальная ошибка потери значащих цифр для аппроксимаций вида (2.2−20) и (2.2−21) равна .

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

Таким образом, ошибка численного определения градиента в (2.2−20) и (2.2−21) определяется ошибками отбрасывания и потери значащих цифр. Но если первая из них пропорциональна конечно-разностному интервалу, то вторая — обратной ему величины. Поэтому последствия варьирования для них противоположны. На практике конечно-разностный интервал выбирается из соображений минимизации вычислимой оценки суммарной ошибки.

Результаты многочисленных исследований [8−11,14,15] показали эффективность использования относительного конечно-разностного интервала, одинакового для всех компонент. При этом оценка абсолютного конечно-разностного интервала определяется выражением:

,.

В настоящее время принята следующая стратегия вычисления градиента. Сначала градиент надо вычислять по менее трудоемкой формуле (2.2−20), а затем, когда последовательность точек попадает в область, где норма градиента мала, надо переключиться на формулу (2.2−21). Переход осуществляется тогда, когда на какой-либо итерации использование алгоритмом односторонней конечно-разностной аппроксимации градиента не приводит к достаточному убыванию функции, но критерии останова еще не выполняются. В этом случае алгоритм продолжает свою работу с использованием центральной аппроксимации градиента до тех пор, пока не выполнятся или критерии останова, или условие обратного перехода. Условие перехода определяется следующим неравенством:

.

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

Методы условной оптимизации

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

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

.

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

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

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

где r — положительное число, а — непрерывные на интервале дифференцируемые функции, причем при. Вспомогательная функция имеет вид:

.

Она определена внутри W и неограниченно возрастает, если для какого-нибудь i. Наиболее часто используются следующие барьерные функции [11,17]:

логарифмическая функция ;

обратная функция .

Чем меньше значение r, тем меньше влияние второго слагаемого в выражении для. Это обеспечивает сходимость при точек к и как следствие существенное убывание значения и совпадение его в пределе с .

Модельная схема оптимизации в этом случае формулируется так [17]:

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

.

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

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

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