УДК 681.327.12
Ю. С. Бехтин
Рязанский государственный радиотехнический университет
Д. В. Титов
Курский государственный технический университет
ВЕЙВЛЕТ-ОБРАБОТКА ИК-ИЗОБРАЖЕНИЙ ДЛЯ КОМПЕНСАЦИИ ДРЕЙФА ВОЛЬТОВОЙ ЧУВСТВИТЕЛЬНОСТИ ЭЛЕМЕНТОВ ФОТОЭЛЕКТРОННЫХ МОДУЛЕЙ
Дрейф вольтовой чувствительности элементов матриц фотоэлектронных модулей рассматривается как мультипликативный шум на изображениях. Целью обработки является текстурная сегментация изображения на уровне вейвлет-преобразования, чтобы в зависимости от типа сегмента применять различные схемы формирования оценок вейвлет-коэффициентов.
Введение. Некоторые отечественные фотоэлектронные модули (ФЭМ) строятся на матрицах инфракрасного (ИК) диапазона, которые имеют достаточно большой (до 30%) разброс вольтовой чувствительности их элементов [1]. Неоднородность профиля чувствительности значительно искажает оригинал изображения, формируемого на выходе ФЭМ с помощью электронного коммутатора приводя, иногда к почти полной неузнаваемости объектов оператором. Чаще всего разработчики ФЭМ предлагают для компенсации разброса чувствительности (или геометрического шума чувствительности) использовать предварительную калибровку по „низкой" и „высокой" температурам [1], которая заключается в оценке коэффициентов передачи электронного канала для каждого элемента матрицы. К сожалению, на практике наблюдается дрейф коэффициентов передачи во время работы ФЭМ, что делает предварительную калибровку неэффективной. Эта проблема является особенно острой для КРТ-фотоприемников (кадмий — ртуть — теллур), где нет охлаждения ФЭМ жидким азотом, температура среды и ФЭМ изменяется с большой частотой и в широких пределах и где прерывания рабочего режима на калибровку производятся каждые 20-30 мин. Очевидно, чтобы увеличить интервал между калибровками, необходимо принимать дополнительные меры по удержанию коэффициентов передачи фотоэлементов в некотором допустимом диапазоне. Одним из вариантов является цифровая обработка изображений, если геометрический шум чувствительности рассматривать как мультипликативный шум с единичным средним, воздействующим на оригинальное изображение.
Поскольку размер кадра на выходе ФЭМ относительно большой, то для подавления мультипликативных помех в реальном времени желательно иметь как можно меньший объем передаваемой информации. Одним из инструментов, где и фильтрация, и сжатие данных производятся одновременно, является двумерное вейвлет-преобразование на основе сепарабель-ных базисов [2].
Статистическая модель геометрического шума чувствительности. В целях выявления вида закона распределения коэффициентов передачи фотоэлементов для большого числа матриц ФЭМ были проанализированы результаты оценок их параметров. На рис. 1 приведена гистограмма коэффициентов передачи фотоэлементов КРТ-матрицы ФЭМ, изготовленного одним из предприятий Москвы. Как видно из рис. 1 матрица содержит дефектные и условно-годные элементы, сигналы с которых заменяются на другие, вычисленные по какому-либо алгоритму компенсации дефектов [1], чаще всего на сигналы соседних, годных элементов.
Убирая из рассмотрения дефектные и условно-годные фотоэлементы, путем усреднения по 30 изделиям была получена типовая гистограмма коэффициентов передачи фотоэлементов матриц ФЭМ. На основании типовой гистограммы была выдвинута гипотеза о гауссовском законе распределения коэффициентов передачи с единичным средним. Гипотеза проверялась по критерию X с альтернативными гипотезами о виде закона распределения: треугольное, двойное экспоненциальное (Лапласа), логнормальное, Стьюдента. Проверка по критерию % подтвердила гауссов характер плотности вероятности коэффициентов передачи годных фотоэлементов матриц ФЭМ.
250
Рис. 1
Модель изображения на уровне вейвлет-преобразования. Дрейф коэффициентов передачи приводит к мультипликативной модели
I = Хг, (1)
где Г — наблюдаемое изображение, X — неизвестный оригинал, 2 — мультипликативный га-уссовский шум с единичным средним. Для оценки поведения трех составляющих модели (1) применяют коэффициенты вариации, которые являются индикаторами однородности изображения [2]. Коэффициент вариации вычисляется в пределах некоторого окна вокруг наблюдаемого пикселя. Поскольку переменные Xи 2 в (1) независимы, то можно записать:
Му = Мх Мг = Мх , (2)
Су = Оу / Мг , Сх = Ох / Мх , С2 = О2 1Мг ,
где С(...) — коэффициенты вариации, вычисляемые с помощью оценок математических ожиданий му, Мх, Мг и дисперсий о2 , о\ , а\ в пределах заданных окон. Обычно окно для оценки коэффициента вариации С2 имеет меньший размер по сравнению с окном для оценки коэффициента вариации Су [2, 3].
С учетом того, что CY = (M [X 2]M [Z2] -цХ)/V-X, значение квадрата коэффициента вариации для истинной интенсивности изображения имеет вид:
C2 = Ч^Ч (3)
x C2 +1
Сравнение оценок коэффициентов вариации CY и C£z позволяет различать три ситуации, возникающие при анализе зашумленного изображения, и, следовательно, проводить его обработку по следующему принципу:
f (y), если Cz < Cy < Cmax(Cx * 0) £ =\ y, если Cy > Cmax(Cx > 0) (4)
, если (§Y < Cz (CX = 0) ,
где Cmax = max{Cz}. Если CY > Cmax (CX > 0), то данная ситуация соответствует попаданию в участок либо с крайне разнородной текстурой, либо ФЭМ с точечными источниками или перепадами яркостей при переходе окном границ объекта; в этом случае необходимо сохранить наблюдаемое значение сигнала без изменений, поскольку влияние шума незначительно. При (§Y < C£z (CX = 0) наблюдается однородная текстура; очевидно, что усредняющий фильтр может восстановить оригинал. Наконец, ситуация C£z < CY < Cmax(CX * 0) является неопределенной с точки зрения проверки на однородность текстуры, следовательно, для получения оценки оригинала должен быть применен какой-либо алгоритм фильтрации f(y).
Поведение вейвлет-коэффициентов при действии шума. Если применить схему вейв-лет-декомпозиции изображения Малла (Mallat) [2], где разложение идет по ветви аппроксимации, то можно определить оператор W[ 1 ], формирующий горизонтальные, вертикальные и диагональные коэффициенты на каждом уровне j. Поскольку вейвлет-преобразование является результатом последовательных сверток, и если после децимации размер окрестности точки (a, b) на каждом масштабе не изменяется, то
W[ 1] = G[i] П H[i], Wy (a, b) = [W [l]Y]a,b, (5)
i=1
где
WY (a, b) = XX HlkGY (a + k, b +1); (6)
k I
здесь k и l определяют величину окрестности вокруг наблюдаемого пикселя, при этом
Hk = hi; G = g[1-1V-1. (7)
Тогда вейвлет-декомпозиция зашумленного изображения с мультипликативной моделью вида (1) может быть представлена следующим образом:
WY = W[ 1 ]Y = W[ 1 ]xz = W[ 1 ]X + W[ 1]X(z -1) = W[ 1]X + W[ 1]xzc = WX + Ws, (8)
где wx = w[1]x, we= w[1](xzc) — центрированные и некоррелированные случайные процессы с нулевыми математическими ожиданиями.
Одним из авторов данной статьи в [3] показано, что на уровне вейвлет-преобразования
С 2
т.. — Си
r2 _ (Q)
Сх _ j+с|)' (Q)
где KlJ] )"' CWy _°Wy / и cwz _ VK2jl cz • Сравнивая соотношения (4) и (Q), можно
n
k l
сделать вывод о возможности использования коэффициентов вариаций Сщ и Сщ для классификации текстуры изображения (сегментации).
Поиск оценок коэффициентов вариации CW и CWz происходит в пределах окрестностей
(окон) каждого вейвлет-коэффициента высокочастотных субполос и соответствующего ему пикселя из низкочастотной области (аппроксимации) на текущем уровне декомпозиции через
оценки <jЩ (a, b) _ — £ w^(k, l) и $Y (a, b) _ — £ Y(k, l), где D — область соседних точек, оп-
Y D (k,l)eD D (k,l)eD
ределяемых окном. Оценки (( можно найти аналогичным образом с меньшим окном или
Wz
использовать оценки коэффициентов вариации Сz , вычисленные по массиву коэффициентов
аппроксимации соответствующего уровня.
Следуя логике соотношения (4), определим возможные ситуации при сравнении оценок коэффициентов вариации вейвлет-коэффициентов:
1) если <Cw2 < CWY < CW max , то WY _ f (wY ) , CWmax _ max {CWZ } ;
2) если CwY > Сщ max , то WX _ WY ;
3) если Сщг < CWz, то wx _ 0 •
Восстановленное (отфильтрованное) изображение X получается в результате обратного вейвлет-преобразования над коэффициентами wx .
Поиск оценок вейвлет-коэффициентов при неоднородной текстуре. Пусть каким-либо образом определены априорные плотности вероятности PW■ (wX) и PW (wY) для случайных процессов WX и WY . Тогда для обеспечения минимально возможного значения байесовского риска следует вычислять такое значение оценки вейвлет-коэффициента Wx , при котором апостериорная условная плотность вероятности PW. \W (wX | wY) принимает максимальное значение:
PWX\WY (wX 1WY) ^max. (10)
WX
Апостериорная условная плотность вероятности определяется по соотношению Байеса [6], которая с учетом (14) принимает вид
п /14 PWY\WX
(WY \ WX ) PWX (WX ) PWа\WX (WS \ WX ) PWX (WX ) пп
PW \W (wX | WY) _—^^-X-_—^^-X-. (11)
Щх|ЩЛ Pwy (Wy) PWY (WY)
Оценка вейвлет-коэффициента М'х , при которой достигается максимум апостериорной условной плотности вероятности, вычисляется с использованием логарифмирования соотношения (11):
d
dwX
(ln(Pws\Wx (wE \ ^х)) + ln(Pwx (WX))) \WX_WX _ 0. (12)
Таким образом, необходимо определить априорную условную плотность вероятности для вейвлет-коэффициентов шума Рш__^ (w„ | wX) и априорную плотность вероятности
Р\¥х (Wx), описывающую поведение неискаженных вейвлет-коэффициентов. Для этого используем обобщенное распределение Гаусса, широко применяемое в последнее время.
Обобщенное распределение Гаусса случайной величины X с нулевым математическим ожиданием имеет вид [4]:
/х(х) = ЖР,с)ехр{-((р,а)|х|)} , -о< х <+оо, р>0, (13)
где с2 — дисперсия, р — параметр, определяющий эксцесс для кривой распределения, при этом
Я(Р, а) =
P-G(P, а)
G (Р, а) =
Г(3/ Р) Г(1/ Р)
1/2
(14)
2Г(1/ Р)
Частными случаями (14) являются распределение Лапласа (р=1) и распределение Гаусса (р =2). Таким образом, чтобы применить выражение (13) для поиска оценки W€x по оценке (12), необходимо найти оценки дисперсий и параметров р и р_. для аппроксимации соответствующих плотностей вероятности вейвлет-коэффициентов (wx) и ^ (| Wx).
Если оценки дисперсий и параметров рх и рн найдены, то оценка М'х согласно (13) определяется через численное решение уравнения (приводится без учета знака при вейвлет-коэффициенте wx )
рвОЕ |w7 - wx |Р2_1 +РxGx¡Wx ^-1 = 0. (15)
Для общего случая оценки дисперсий и параметров р x и р~ могут быть найдены через соответствующие центральные моменты выборочных данных.
Параметр обобщенного распределения Гаусса р вычисляется через численное решение уравнения [4]:
Г(3/ р) . (16)
л/Г(1/ Р) Г (5/ р)
Для вычисления второго и четвертого моментов процессов X и 5 используем тот факт, что оценки локальных дисперсий с^ найдены в процессе анализа текстуры, при этом
Ц2,^ =с1¥г , = -1 X М'^(к'1).
(k ,l )eD
Результаты моделирования. Для проверки разработанного алгоритма в среде Matlab была сформирована последовательность кадров, где геометрический шум чувствительности моделировался с помощью датчика нормально распределенных случайных чисел с единичным средним и дисперсией, которая увеличивалась от кадра к кадру по линейному закону. В эксперименте использовалось быстрое трехуровневое вейвлет-преобразование, тип вейвле-та — CDF 9.7.
Рис. 2
Рис. 3
На рис. 2 и 3 представлены некоторые зашумленные и обработанные кадры из видеопоследовательности соответственно. Изменения средней квадратической ошибки (СКО) для искаженных и восстановленных кадров показаны на рис. 4 где показаны изменения СКО по кадрам (в увеличенном временном масштабе): без обработки (кривая 1), с обработкой предложенным алгоритмом (кривая 2).
ат
Рис. 4
Из рис. 3 и 4 видно, что предложенный алгоритм вейвлет-обработки изображений матричных ФЭМ обеспечивает относительно высокое качество восстановления оригинала, что позволяет увеличить интервал между калибровками ФЭМ до нескольких часов.
список литературы
1. Кругликов С. В. Методы и средства подавления структурных помех многоэлементных фотоприемников // Аналитический обзор № 4628 за 1970—1987 гг. М., 1989.
2. Гонсалес Р., Вудс Р. Цифровая обработка изображений. М.: Техносфера, 2006.
3. БехтинЮ. С. Алгоритм вейвлет-фильтрации зашумленных изображений // Вестник РГРТА. № 15. 2004.
4. Multiresolution MAP Despeckling of SAR Images Based on Locally Adaptive Generalized Gaussian pdf Modeling / F. Argenti, T. Bianchi, L. Alparone // IEEE Trans. on Image Processing. 2006. Vol. 15, No. 11. P. 3385—3399.
Рекомендована кафедрой Поступила в редакцию
вычислительной техники 01.09.07