CC BY
41
В.А.Фурсов
ИДЕНТИФИКАЦИЯ ОПТИЧЕСКИХ ИСКАЖАЮЩИХ СИСТЕМ С ОТБОРОМ ИНФОРМАТИВНЫХ ФРАГМЕНТОВ ИЗОБРАЖЕНИЙ
ВВЕДЕНИЕ
Задача параметрической идентификации оптического канала состоит в определении параметров модели оптической искажающей системы по измерениям полей яркости неискаженного (входного) и искаженного (выходного) изображений. Обычно идентификация осуществляется эпизодически в ситуациях, когда наряду с прошедшим через оптическую систему выходным изображением появляется возможность фиксировать входное (неискаженное) изображение. Если характеристики искажающей системы в течение некоторого достаточно длительного промежутка времени после этого остаются неизменными, результаты идентификации могут использоваться для коррекции последующих искаженных изображений. Важнейшими требованиями, предъявляемыми к алгоритмам идентификации, являются точность и надежность определения характеристик.
Источником погрешностей идентификации искажающей системы являются ошибки в исходных данных, связанные как с погрешностями измерений полей яркости входного и выходного изображений, так и с ограниченностью размерности модели (ошибки аппроксимации). Погрешность в каждом отдельном эпизоде идентификации зависит от чувствительности к этим ошибкам. В частности, если задача плохо обуслоыена, то ошибки определения характеристик мот быть чрезвычайно высоки даже при достаточно точных измерениях полей яркости. Поэтому прсдста&тястся важным построение эффективных процедур анализа информативности данных с целью оценки возможности определения характеристик искажающей системы с требуемой точностью.
1. ПОСТАНОВКА ЗАДАЧИ
Модель искажающей оптической системы строится в виде разностного уравнения: фильтра с конечной импульсной характеристикой (КИХ-фильтра) априори заданного порядка. Если КИХ-фильтр обладает опорной областью О: (к,.к,) е Э, то связь выходных
отсчетов у(п],п2) с входными описывается соотношением [1]:
У(",.»2)= X Еь(к|.к2)х(п1-к1,п2-к2)^(п1,п2), (1)
где Цк^кД х(п - к(,п2-к2) - значения импульсной характеристики и входных
отсчетов в точках опорной области О соответственно , а ^|п1,п2| - ошибка, включающая
погрешности измерений точек поля яркости и ограничения на порядок модели.
Разностное уравнение (1) можно переписать в компактном векторном виде. Пусть ш - номер опорной области, в которой произведены наблюдения, тогда, введя обозначения:
Уш = У (" I»" г)» ^т = " 1' " 2 ) " *** скаЛЯРов- и
Хт =[ Х(П1'П2).....Х(П,-к,'П2"к2)'-Г'
ст =[ь(01>02),---5 - для векторов размерности N х 1 , где N - число
различных отсчетов в опорной области Э: (к,,к2) е О, в соответствии с (1) имеем:
ут = х!с + £ .
' т т ^т
Для М опорных областей получаем векторно-матричное соотношение:
(2)
У =
V '-Л V
= с +
.Ум. .^М -
= Бс + Е.
(3)
где Б - матрица М х N, а У, Е - векторы М х 1.
Задача идентификации заключается в определении оценки с вектора параметров с по наблюдениям Б, У при наличии ошибок в исходных данных - Е. В частности, оценка метода наименьших квадратов (МНК) получается в результате решения уравнения |2)
(4)
5т5с = Бт У, при этом ошибка идентификации Дс = [8т8]_18тЕ.
(5)
В соответствии с (5) ошибка Дс зависит от величины и направления вектора ошибок - Е, а также от свойств (обусловленности) матричного коэффициента усиления А = 5Т5.
Обычно, из соображений экономии памяти, идентификацию характеристик и обработку осуществляют последовательно, выделяя фрагменты изображений небольшой размерности. Интуитивно ясно, что высокая чувствительность к ошибкам в исходных данных, вызываемая плохой обусловленностью матрицы А, будет иметь место на фрагментах, не содержащих информативных элементов изображений, например, ровный "серый" фон. Отбраковка таких фрагментов путем предварительного визуального просмотра требует определенного опыта и затрат времени. Поэтому представляется важным построение простых формальных процедур анализа обусловленности задачи с
целью определения фрагментов изображений, на которых возможно определение характеристик искажающей системы с требуемой точностью.
Причиной плохой обусловленности матрицы А = является
мультиколлинеарность |2], понимаемая как сопряженность переменных, я&ляющихся компонентами векторов = матрицы Б = [в,^!...^]. Характеристикой
обусловленности являются следующие меры [2]: определитель матрицы А - <1е1:(А), минимальное собственное значение - ¿-„^(А) и спектральное число обусловленности к(А) = Хтах(А)/Хт(п( А), где - минимальное и максимальное собственные
значения матрицы А.
С использованием указанных мер можно прогнозировать погрешность решения задачи идентификации, если есть надежные предположения относительно величины нормы вектора ошибок в исходных данных. В частности, оценки сверху и снизу для квадрата эвклидовой нормы вектора даются неравенствами |3|:
Мг^"1,1Мг. о
Известна также оценка сверху для относительных ошибок в решениях |4]
|Дс|| К (А)
<-~~7~7-к + ч), (8)
1 - К(А)Ьд у ь А/
. 5А || §Ь ,
где а = || д.и' > "ь = Т^Т " °тн°сительнь1е ошибки задания матрицы А и вектора Ь:
X т
= А - А*, 6Ь = Ь - Ь*; а А* = Б* Б" и Ь*=5* У* - их точные значения. Из приведенных соотношений видно, что решающее значение при формировании ошибок идентификации имеют собственные значения матрицы А = 5Г5. В частности, если хотя бы одно из них близко к нулю, то соответствующие коэффициенты усиления ошибок становятся сколь угодно большими.
Попытка прогнозировать погрешность идентификации по соотношениям (6)-(8) обречена на неудачу. Во-первых, определение Хтт(А) и к(А) по вычислительной
сложности не уступает, а часто превосходит сложность собственно определения вектора с . Но даже если пойти на эти вычислительные затраты, то нужный эффект вряд ли будет Достигнут. Если задача действительно плохо обусловлена, то вычисление достаточно точных значений Хтш(А) и к(А) становится трудноразрешимой проблемой. В настоящей
работе предлагаются достаточно простые и эффективные (по сравнению с известными правилами Гершгорина и Адамара) методы обнаружения плохой обусловленности непосредственно по элементам матрицы А, с использованием которых далее строится процедура отбора информативных фрагментов изображения.
2. ПОСТРОЕНИЕ ОЦЕНОК ОБУСЛОВЛЕННОСТИ
Имеет место следующее утверждение. Пусть - семейство симметричных
вещественных ограниченных неотрицательно-определенных (типа А = БГ5) матриц порядка Ы: Л = {а: БрА* = с!к, к = 1,у, у< ы}, где с1к - заданные (в допустимой области)
константы.
Тогда максимальное и минимальное собственные значения для любой матрицы А из семейства А е удовлетворяют системе
Ут.Хк +\к =сГ, к = 1Я (9)
' I 1 V к
¡-1
где т., i = 1.N - 1 - целые положительные числа, удовлетворяющие условию
v-i
Zmi = N-l. (10)
¡-i
Сформулированное утверждение является результатом решения следующей задачи: найти экстремумы функции F{X) = X. для произвольного ie[l,N] при v ограничениях
вида
n
= k = T~v. (11)
¡«i
Отмстим, что ограничения (11) непосредственно вытекают из следующих свойств матриц
n
рассматриваемого класса: Ак) = Хк(А) и = SpA, а числа dk=SpAk легко
¡=i
вычисляются по элементам матрицы А = (а.. ) для любого к.
Из равенств (9), (10) вытекают следующие частные случаи. Если v = N. то, в силу (10) все m. = 1, i = l,N и система уравнений (9) может быть представлена в
единственном возможном варианте. n
к = Щ (12)
ы
Если вдобавок все числа Х.(А) различны, решение системы (12) существует и дает
решение полной проблемы собственных значений. При v=2 экстремальные собственные значения удовлетворяют системе
(N-lU, +XN =dr
(N-i)x;+x2N=d2, (13)
которая имеет аналитическое решение. В частности, имеют место следующие оценки для собственных значений матриц из множества
Л, = (а: SpA = d, = Const, SpA2 =d2 =Const):
maxX(A) = -Ml + V(Nß-l)(N-l)|, (14)
AcJ3 N 1 >
min MA) = ^fl - V(Nß - l)(N -1)1, (15)
Ле<4; N 1 J
d2 SpA2 1 1
где ß = -f = ---y, — sßs-. (16)
d2 (SpA)2 N N -1
Оценка (14) соответствует случаю, когда в (13) (N-1) совпадающих собственных значений X, меньше XN, являющегося в данном случае максимальным. Во втором случае - (15) собственное значение X , имеющее кратность (N-1), больше XN, которое теперь, в
свою очередь, явдяется минимально возможным.
С использованием опенок (14), (15) можно построить также верхнюю границу для числа обусловленности:
Равенство в (17) достигается лишь в случае р=Г<Л при этом К(А)=1. Во всех остальных случаях имеет место строгое неравенство, т.к. оценю! сверху (14) и снизу (15) соответствуют двум различным решениям (13).
При возрастании N и плохой обусловленности задачи оценки (14), (15), (17) могут не существовать в силу указанных в (16) ограничений на величину р. В контексте поставленной здесь задачи это не является существенным недостатком. Поскольку целью построения оценок является установление факта плохой обусловленности задачи, это как раз будет сигнализировать о недостаточной информативности исходных данных и необходимости перехода к формированию новой выборки наблюдений или уменьшения размерности модели.
3. ПОСТРОЕНИЕ ПРОЦЕДУРЫ ОТБОРА ИНФОРМАТИВНЫХ ФРАГМЕНТОВ
Подставляя оценки (14), (15), (17) в (6)-(8), можно построить соответствующие им оценки сверху и снизу для ошибок определения искомых параметров. Однако для построения решающего правила типа "разрешения" можно предложить более простой путь. Из (14), (15), (17) видно, что оценки для мер обусловленности являются
ей величину:
монотонными функциями величины Р = с17/с12. Удобно вместо р рассматривать обратную
2
р БрА 2.аи 2Л
Назовем характеристику Ф(А) показателем диагонального преобладания.
Непосредственно из равенств (18) легко устанавливается связь показателя Диагонального преобладания с известными мерами обусловленности. В частности, имеют место следующие соответствия:
если Ф(А) —► 1, то X ->0, К(А)->оо, dct(A)-»0,
(19)
/1 ( 1 I14
же Ф( А) -> N, то *т|П(А)->—, К( А) -> 1, аег(А) ] . (20)
с с ли
где А = А/БрА. Зависимости обратного свойства также возможны, но не обязательно.
Например, для особенной матрицы А, у которой одно собственное значение равно нулю, а N-1 собственных значений отличны от нуля и равны между собой, Ф(А) = N - 1, хотя А
вырождена. Вместе с тем выполнение неравенства
Ф( А) > N - 1 (21)
гарантирует невырожденность матрицы А. Более того, хорошая обусловленность, обычно, имеет место даже при незначительном превышении величины N-1.
Показатель диагонального преобладания имеет следующие важные свойства. Он инвариантен к умножению матрицы А (и Б) на скаляр. Однако если один из векторов-столбцов матрицы Б умножить на большое число, то Ф(А)->1, т.е показатель Ф(А)
чувствителен к большим различиям в масштабах измерений. В этом смысле он является интегральной характеристикой обусловленности. Однако это свойство часто служит причиной того, что "почти ортогональность" векторов е., 1 = I,N матрицы 5 может быть
замаскирована эффектом разницы в их нормах. Для установления степени сопряженности векторов V в данном случае показатель диагонального преобладания следует вычислять
для соответствующей матрицы сопряженности Я с элементами г = ^8. .я.)/
р. |. Ясно.
что если для всех \ * ], 1,3 = 1, N, г = 0, то достигается максимальное значение показателя Ф(А)=М. Если же все г = 1, ¡,] = то Ф(А) = 1 (матрица сопряженности Я
и соответствующая ей матрица А вырождены).
На основании свойств (19)-(21) предложена процедура отбора наиболее информативных фрагментов изображения, показанная на рис.1. Здесь показатель диагонального преобладания используется как самостоятельная мера обусловленности. В частности, величина показателя Ф(А) вычисляется для всех последовательно
просматриваемых фрагментов изображения. Сохраняются лишь те данные, для которых величина Ф(А) превышает ее значения на предыдущих фрагментах. После того как просмотрен последний фрагмент, проверяются условия типа Ф>Ф где Ф,
цт цт
допустимая величина показателя диагонального преобладания, Ф(А) и (или) Ф(Ю в данной задаче.
Экспериментально установлено, что в большинстве случаев целесообразно проверять два условия: первое - Ф(К)>М-1+е, где е - заданное малое положительное число, и второе - Ф( А) больше (1.1-5-1.2). Первое условие гарантирует невырожденность задачи, а второе контролирует недопустимое ухудшение обусловленности из-за большой разницы величин норм векторов, составляющих матрицу Б. Если требуемое значение показателя диагонального преобладания не достигается на данном изображении, следует принять меры к улучшению обусловленности, в частности, можно, например, подобрать подходящую опорную область, приводящую к снижению размерности задачи.
Рис 1. Общая схема алгоритма оценивания с отбором информативных фрагментов изображений.
4. ПРИДАНИЕ УСТОЙЧИВОСТИ К АНОМАЛЬНЫМ ОШИБКАМ
В соответствии со схемой рис.1 после выбора матрицы А и вектора Ь, соответствующих набору данных на фрагменте с достаточно хорошей обусловленностью, осуществляется переход к решению системы уравнений вида (4). Однако погрешность (5) решения задачи может быть достаточно большой даже при хорошей обусловленности задачи, если в исходных данных содержатся выделяющиеся (аномальные) ошибки. Более тогс^аномальные ошибки могут замаскировать факт плохой обусловленности, так как они обладают сильными регуляризующими свойствами. Поэтому показанная на рис.1 схема должна применяться в сочетании с мерами подавления аномальных ошибок.
Алгоритм, устойчивый к аномальным ошибкам, строится в виде процедуры робастного Ц-оценивания типа итерационного МНК [2]. В общем случае эта процедура
выглядит следующим образом [2,3). Вначале ищется обычная оценка МНК как решение Уравнения (4). Затем вычисляются невязки
е = у -хтс, т-ЦМ (22)
т т т ' '
и строится диагональная (весовая) матрица
порядка М. каждый элемент которой яапястся функцией соответствующей невязки: g = g(e )= е , т = 1,М, где 0^у<2.
°т с\ ш/ т
Далее строится новая система уравнений вида
5тОБс = 5тОУ, (24)
решение которой дает следующее приближение. С использованием вновь полученной оценки снова вычисляются невязки (22), весовая матрица (23) и строится система уравнений (24). Итерационный процесс останавливается, если
||ск-сЛ*8 или к>кЦт)
где с ск , - оценки, полученные на соседних итерациях, 5 - заданное (достаточно малое) положительное число, а кЦт - заданное допустимое число итераций.
В данном случае предложена трехшаговая процедура с изменяющимися от итерации к итерации значениями параметра у (у= 2,1,0). В частности, после получения оценок
МНК (1-й шаг, у=2) осуществляются еще две итерации с использованием следующих процедур формирования элементов матрицы О:
2-й шаг при у = 1 йт =
, (25)
е + ш I 0
3-й шаг при у=0 к = --. (26)
т О
Здесь е0 - малое положительное число, позволяющее избежать деления на нуль. Оно
вычисляется после каждой итерации в зависимости от величины нормы текущего вектора
м
невязок. Константа g0 выбирается из условия нормировки (3]: В =М, которое
Ш=1
реализуется в два этапа. Вначале вычисляются весовые коэффициенты по формулам
г м V1
т
«т
(25), (26) при Е0 = 1. Затем определяется новый нормирующий множитель = N и все коэффициенты умножаются на него: = .
5. РЕЗУЛЬТАТЫ ЭКСПЕРИМЕНТАЛЬНОЙ ПРОВЕРКИ
Решалась задача определения параметров импульсной характеристики по тестовым изображениям размером 128 х 128 |5|. Входное (неискаженное) изображение
формировалось в виде прямоугольников на постоянном фоне (рис.2). Размеры сторон прямоугольников и их яркости задавались в виде случайных равномерно распределенных чисел в диапазонах 1-И0 и 04-255 соответственно. Общее количество элементов (прямоугольников) задавалось равным 40. Они располагались на поле изображения случайным образом: по горизонтали в соответствии с равномерным законом распределения и с возрастающей плотностью вниз по вертикали с ростом порядкового номера фрагмента.
Выходное изображение формировалось "размазыванием" входного с использованием разностного уравнения вида:
у(п,,п2) = -5>(к)х(к),
(27)
к=0
где х(к) = х(п1-к,п2) + х(п1н-к,п2) + х(п ,п,-к) + х(п],п +к) .
Значения импульсной характеристики искажающего фильтра задавались равными
Ь(0)—0.5. 1)(1)=0.3, 11(2)=0.2. Уравнение (27) соответствует опорной области, показанной
на рис.3. В данном случае для нас имеет значение лишь то, что импульсный отклик в направлении осей N , 1М., опорной области одинаков, хотя частотная характеристика и не
обладает круговой симметрией.
N ч
N
Рис. 2 Рис. 3
Тестовое изображние Опорная область
Обрабатывавшиеся в каждом эпизоде идентификации фрагменты изображения имели вид полосок 5x128. При этом матрица Я для каждого фрагмента имела размерность 25x3. а соответствующая ей матрица А - 3x3. Искомый вектор оценок параметров с = [ с{0), с(1), с(2)]Т = [¿(О), Ь(1), Ь(2)1 определялся с использованием
описанного выше итерационного МЫ К.
Проведено два вида экспериментов: при отсутствии ошибок в исходных данных и при наличии импульсных помех. В первом эксперименте выявлялись величины показателей Ф( А) и Ф(Ю, при которых существенными оказываются даже ошибки
округления, а также степень влияния на точность идентификации робастной трехшаговои процедуры при отсутствии помех. Во втором эксперименте искаженное изображение дополнительно подвергаюсь аддитивным импульсным помехам с вероятностью 0,04. Цель этого эксперимента - иллюстрация эффективности трехшаговои робастной процедуры по сравнению с обычным МНК при появлении в исходных данных аномальных ошибок.
Результаты экспериментов приведены в таблице I. Здесь в первых двух (после номера фрагмента) столбцах приведены показатели диагонального преобладания ф(а| и вычисленные для всех фрагментов изображения, а в следующих - относительные среднеквадратическис ошибки оценивания на каждой итерации в первом и втором эксперименте соответственно (величины, соответствующие данным первою эксперимента, следует умножать на 10'5).
Таблица 1.
N фрагмента Ф(А) ф(Ю Относительные среднеквадратические ошибки (СКО)
Без помех (х10~5) При импульсных помехах
МНК 2 шаг 3 шаг МНК 2 шаг 3 шаг
1 1.000 1.000 _ _ _ - - -
2 1.005 1.735 _ _ _ - - -
3 1.027 1.036 0.10 0.30 1.40 0.000001 0.000003 0.000001
4 1.014 1.760 _ - - - - -
5 1.014 1.647 - _ - - - -
6 1.000 1.000 _ _ _ 0.027200 0.000250 0.000038
7 1.087 1.771 _ - 0.093880 0.002260 0.000037
8 1.000 1.000 _ _ _ _ _ -
9 1.093 2.136 0.20 0.90 0.20 21.99110 1.216270 0.000889
10 1.012 1.013 0.30 13.9 1.10 1.177410 0.034750 0.000056
11 1.220 2.053 0.30 0.10 0.00 0.267070 0.008150 0.000020
12 1.040 1.041 0.40 1.40 0.30 0.000004 0.000140 0.000003
13 1.050 2.742 0.10 0.10 0.20 6.370910 0.350150 0.000660
14 1.069 1.076 0.10 0.40 0.10 0.000001 0.000000 0.000001
15 1.024 1.023 4.30 38.2 0.20 0.592780 0.400750 0.174790
16 1.015 1.780 _ _ _ _
17 1.034 1.091 0.30 0.50 0.80 2.568100 0.430770 0.002190
18 1.022 1.051 0.50 0.50 0.90 0.737610 0.022180 0.000051
19 1.251 1.679 0.30 3.80 1.80 0.000003 0.000078 0.000018
20 1.093 1.263 4.80 455. 21.9 0.000000 0.000000 0.000003
21 1.134 1.147 0.20 1.10 0.40 2.244870 0.419580 0.038400
22 1.132 1.132 0.40 0.80 0.10 0.221670 0.049650 0.000480
23 1.032 1.043 0.40 1.60 0.60 0.000005 0.000013 0.000003
24 1.064 1.067 0.20 0.80 0.20 0.157280 0.000111 0.000016
25 1.000 1.000 - - - - -
Незаполненные графы в таблице соответствуют случаям аварийного завершения программы вследствие плохой обусловленности задачи. Отметим, что это произошло как раз на фрагментах, содержащих небольшое число (либо не содержащих вовсе) информативных элементов (прямоугольников). При этом показатель Ф(А) (а в некоторых случаях и Ф(Д)). как и следовало ожидать, оказался близким к единице. Интересно, что в двух случаях (фрагменты 6 и 7) после добавления в исходные данные импульсных помех не произошло аварийного завершения программы, несмотря на плохую обусловленность задачи. Это явлилось следствием того, что аномальные ошибки увеличили "угол" между вектором У и пространством матрицы Б.
ЗАКЛЮЧЕНИЕ
Данные экспериментов подтверждают факт тесной связи показателя диагонального преобладания с обусловленностью задачи. Поскольку тестовое изображение формировалось с использованием случайного набора случайно расположенных элементов, предложенные процедуры, по-видимому, будут работоспособны для достаточно широкого класса изображений.
Работа выполнена в лаборатории математических методов обработки изображений ИСОИ РАН.
Литература
1. Даджион Д., Мерсеро Р. Цифровая обработка многомерных сигналов. Пер. с англ. М.: Мир. 1988, 488 с.
2. Демиденко Е.З. Линейная и нелинейная регрессии.- М.: "Финансы и статистика", 1981, 303 с.
3. Фурсов В.А. Анализ точности и построение алгоритмов идентификации по малому числу наблюдений.- Изв. АН СССР. Техн. кибернетика, N6, 1991, с.130-135.
4. Ланкастер П. Теория матриц.- М.: Наука, 1978, 280 с.
5. Fursov V.A. Identification of optical distorting systems with selecting image informative fragments. - Workshop on Digital Image Processing and Computer Graphics. Proceedings SPIE. - 1994.- vol. 2363.
, r-
Вниманию читателей !
Международный центр научной н технической информации
принимает заказы на международный журнал
"ПРОБЛЕМЫ МАШИНОСТРОЕНИЯ И АВТОМАТИЗАЦИИ"
Главный редактор академик К.В.Фролов
Периодичность - 6 номеров в год Журнал распространяется только по заказам
Цена одного номера 3 USD Цена одного годового комплекта 15 USD С Оплата в рублях по курсу ЦБ РФ на момент расчетов)
Заказы принимаются по адресу: Россия, 125252. Москва, ул.Куусинена, 216, МЦНТИ, СОПИ Телефакс: (095) 943-00-89 Справки по телефону: 198-72-10
