С.М. Широков, И.В. Григоров
МЕТОД ПОДАВЛЕНИЯ ИМПУЛЬСНЫХ ПОМЕХ ПРИ ОБРАБОТКЕ СИГНАЛОВ И ИЗОБРАЖЕНИЙ С ИСПОЛЬЗОВАНИЕМ НЕЛИНЕЙНЫХ ФАЗОВЫХ ФИЛЬТРОВ
ВВЕДЕНИЕ
Импульсные помехи (ИП) являются одним из основных факторов, ухудшающих качество изображений, и могут возникать как при передаче сигналов, несущих информацию об изображениях, по каналам связи ( в особенности коммутируемым проводным каналам низкого качества и радиоканалам), так и непосредственно в процессе их обработки (например, вследствие ошибок декодирования) [1-3]. Если ИП носят точечный характер, их влияние удается сравнительно просто свести к минимуму путем "исправления" пораженных помехой элементов изображения с помощью интерполяции по соседним элементам [1] или медианной фильтрации [3]. При обработке одномерных сигналов на практике широко используется ограничение, бланкирование и другие подобные нелинейные операции подавления ИП [4]. Однако все эти методы становятся малоэффективными, если помеха поражает не один, а достаточно много соседних элементов, т.е. становится соизмеримой по размерам с изображением (или по длительности - с сигналом).
Для подавления таких видов ИП в принципе можно использовать методы нелинейной фильтрации по минимуму среднего квадрата отклонения (СКО), основанные на представлении сигналов и изображений марковскими моделями в пространстве состояний [5,6]. Но возможности их практической реализации весьма ограничены - как из-за отсутствия точных решений у соответствующих уравнений для апостериорных распределений и необходимости использования многочисленных и не всегда оправданных приближений, так и потому, что вопрос о применимости марковских моделей для изображений часто остается открытым [6]. Кроме того, критерий минимума СКО далеко не всегда приемлем для оценки качества изображений и других видов непрерывных сообщений, особенно, если речь идет о субъективном восприятии.
В то же время значительный прогресс, достигнутый за последние годы в технике цифровой обработки сигналов, появление процессоров и других видов СБИС с высоким быстродействием значительно расширяют возможности реализации более сложных алгоритмов обработки, чем перечисленные выше простейшие нелинейные преобразования. Один из таких алгоритмов предложен авторами в [7] и включает в себя такие операции, которые позволяют сжать импульсы помехи до точечных практически без изменения полезного сигнала, что дает возможность затем подавить их одним из известных указанных выше методов. Цель данной работы - его обобщение на двумерные сигналы (изображения и поля).
1. МОДЕЛИ ИМПУЛЬСНЫХ ПОМЕХ
Импульсные помехи в каналах связи и на изображениях отличаются от других видов шумов тем, что имеют значительно меньшую длительность или размеры по сравнению с анализируемой областью одномерного или двумерного сигнала [1,8]. Как правило, их можно рассматривать как квазидетерминированные функции времени или пространственных координат, форма которых известна, а случайными являются только параметры. В зависимости от вида сигнала (модулированный сигнал в канале связи, комплексная амплитуда волнового поля, яркость точки двумерного изображения и т.п.) ИП может принимать значения на множествах комплексных, вещественных или положительных чисел. Используя, как наиболее общее, комплексное представление, можно записать выражение для ИП в виде
и (^ ¿2 ) = Е А* х
* (1)
х Ч ( - ¿1* , ¿2 - ¿2* ; Т1к , Т2к ) ^ ((* ) где д (¿1 ; Т1, т2) - детерминированная функция ,
описывающая форму огибающей ИП, А , t , I , т ,
k ^ 2k ^
т , - случайные параметры, имеющие смысл амплитуды, координат, размеров и фазы *-го импульса. Как правило, их можно считать независимыми с известными плотностями вероятностей. В рассматриваемых здесь задачах свойства последовательностей ИП в целом не существенны и алгоритмы фильтрации синтезируются с учетом характеристик отдельных импульсов. При этом начало координат можно связать с центром анализируемого импульса, так что модель одиночной ИП принимает вид
и (I, ©) = Ад (I, т) exp (ф*),
(2)
где I = (¿1 , ¿2) е Т - вектор координат из некоторой двумерной области Т, Тп = (т^ т2) и © = ( А, т1, т2, ф ) - векторы случайных параметров. В качестве простого приближения для формы ИП на изображениях можно использовать гауссовскую функцию
д ( т)=exp ^т - т
(3)
В каналах связи ИП обычно имеют осцилляции [8], для учета которых можно использовать модель вида
д ( т ) = exp
(
*„ жt
\
Sin ■
(4)
где г - одномерная временная координата, к^ - коэффициент формы, характеризующий затухание ИП. При кф>>0 ИП практически не имеет осцилля-ций, а при кф = 0 осцилляции принимают вид незатухающего гармонического колебания. Результаты статистического анализа реальных радиопомех [6,8] показывают, что для их амплитуд типично логнор-мальное распределение
1 ( 1п2 (А /
•(А) =
Аа\[2я
ехр
2ст2
(5)
где с, ц - параметры распределения. Закон распределения фаз импульсов часто принимается равномерным, а длительностей - усеченным нормальным
Сехр[-(ги - тт)2];
т ^ т ^ т
(6)
w
Т ) =
2с
Т1 < Т2
0;
где С - нормирующая константа, определяемая гра-
2
ницами усечения, т, с -моменты исходной (неусеченной) гауссовской функции.
2. ВЫБОР ОПТИМАЛЬНОГО ПРЕДСТАВЛЕНИЯ СИГНАЛОВ И ПОМЕХ В ЗАДАЧАХ АМПЛИТУДНОЙ СЕЛЕКЦИИ
На вход приемного устройства поступает смесь сигнала 5 (1), суммы сосредоточенной и флук-туационной помех п (1) с ИП и (1)
2 (1) = 5 (1) + П (1)+ и (1) . (7)
В наиболее общем виде рассматриваемый здесь класс нелинейных фильтров можно представить структурной схемой, показанной на рис. 1. Оператор Е задает отображение множества входных воздействий вида (7) на множество функций некоторой (в общем случае - новой) переменной , в котором обеспечена наилучшая селекция сигнала и импульсных помех, определяемая оператором К. Обратный оператор Е обеспечивает формирование оценки
?о (1 ) = 2 (1)- и (1)
остатка смеси (7) с подавленной ИП
2о (1 ) = 5 (1) + п (1) (8)
по селектированной функции Уо (&), а оператор W описывает линейный фильтр, обеспечивающий оптимальную оценку полезного сигнала € (1) в смеси (8).
¥ К г 0\" 1 г1 \У
ад
Рис. 1
Критерием оптимальности блоков нелинейной селекции и фильтра в целом при обработке временных сигналов обычно является минимум СКО
= (2о - го (
где
= (5 - € (2
€ = е-1КЕг.
(9)
(10)
(11)
В задачах обработки изображений, а также при приеме непрерывных сообщений, часто используется критерий приближения в равномерной метрике
д2 =
шах|5
г, о
С)- )|2
(12)
где О - пространство реализаций.
В отличие от задач оптимальной линейной фильтрации, где выбор базиса для представления сигналов не влияет на результат и определяется лишь удобствами реализации (например, применением БПФ), качество амплитудной селекции сигналов и ИП существенно зависит от такого выбора. Например, селекция с помощью ограничения или бланкирования ИП эффективна во временной или пространственной области, где сигнал и помеха различаются по амплитуде и длительности, и теряет смысл в спектральной области, где их спектры сходны.
Таким образом, возникает задача выбора оптимального отображения Е, при котором с учетом вида оператора селекции К обеспечивается наиболее эффективное подавление ИП по критерию (9).
Как и в других задачах оптимальной обработки сигналов, из общего критерия (9) можно вывести более простое и удобное для синтеза алгоритмов решающее правило. Оно зависит от выбранного способа селекции (оператора К).
Нетрудно показать, что при использовании ограничения указанное правило сводится к максимизации показателя селективности
7= Ри/Ро:
(13)
где Р и Ро - усредненные по реализациям пиковые мощности ИП и остаточной смеси вида (8) в области отображений (функций V (-Э) ).
При использовании бланкирования необходимо минимизировать среднюю ширину ИП в области отображений Зп .
Обоим правилам удовлетворяет отображение Е, обеспечивающее максимальное увеличение амплитуды и уменьшение длительности ИП или площади пораженного ИП участка изображения в смеси (7). Детальный анализ показывает, что такое отображение не может быть линейным и должно зависеть от амплитуды или мгновенной мощности 2(1) (в
2
е
противном случае, возможно пропорциональное увеличение пиков помехи и сигнала). Кроме того, должна быть простой и однозначной реализация обратного оператора Е . Этому требованию в наибольшей степени удовлетворяют унитарные (ортогональные) операторы.
Таким образом, отображение Е, удовлетворяющее поставленным требованиям, следует искать в классе операторов с унитарной нелинейностью. Один из классов таких операторов для функций непрерывных аргументов рассматривается в нелинейной квантовой механике [9], нелинейной оптике [10] и описывается эволюционными уравнениями шре-дингеровского типа
д
г-+ аДУ + / (у) = 0,
дп
(14)
где а - постоянный коэффициент, Д - оператор Лапласа по координатам ^ и ¿2 , /- функция, определяющая вид нелинейности, У(п, г) - комплексная функция, значение которой при п= 0 и некотором п= I соответствует сигналам на входе и выходе блока Е:
у( о, г ) = г (г); у(/, г ) = V (г).
(15)
Уравнение (14) задает отображение функций, определенных в области Т , которое описывает преобразование сигнала г (г) в некотором двумерном пространственно распределенном нелинейном фильтре.
При этом обратное отображение задается сопряженным уравнением
дУ
где
т. е.
-г-+ аДУ + / (У) = 0,
у(0, г ) = ^ (г); У(1, г ) = ^ (г),
Е-1 = Е
(16)
(17)
(18)
Непосредственная реализация указанных операторов в аналоговой форме возможна для двумерных и одномерных полей оптического и некоторых других ( например, СВЧ ) диапазонов электромагнитных волн в нелинейных диэлектрических средах. При определенном выборе нелинейности, описываемой функцией /(У) в (14), можно за счет самовоздействия волн обеспечить селективную фокусировку участков поля с аномально большой интенсивностью, обусловленной наложением ИП, и тем самым - увеличение показателя селективности (13). Однако, возможности варьирования и оптимизации параметров нелинейного оператора в этом случае весьма ограничены. Поэтому основное применение рассматриваемый метод может найти при дискретной и цифровой обработке сигналов.
Для такой реализации оператор Е целесообразно представить, используя метод расщепления по физическим факторам [11,12], в виде произведения линейных (Ск) и нелинейных (Нк) операторов вида
Е = П С*Н* [У* ]. (19)
* =1
Каждый из линейных операторов Ск реализуется звеном с передаточной функцией
С* (®) = ехр {-/'а(®12 + ю\ ) Дп* (20)
которой соответствует импульсная характеристика (функция Грина линейной части уравнения (14))
8* (г) = £0* ехР ■
га*
(¿2 + ) 2
(21)
где ю = ( ® , ) - вектор частот, п - число пар звеньев,
Дп* = П*+1 -п*,
-1
80* = (ДП* аг )2: 1
ак =-.
2Дп* а
(22)
Нетрудно заметить, что линейный интегральный оператор с ядром (21) представляет собой известное преобразование Френеля [2] и , таким образом, указанные линейные звенья - это фильтры Френеля. Нелинейные операторы Нк [Ук] в (19) соответствуют умножению функции Ук(т) = У(Пк,т) на зависящие от нее коэффициенты преобразования
Н* [У*]= ехр [/(У*)]. (23)
Сопряженный оператор вследствие уже упоминавшегося свойства унитарности представляется в виде
Е* [У] = П Н* [У * ]С*.
(24)
Как видно из (20) и (23), рассматриваемые линейные и нелинейные звенья имеют единичные, не зависящие от частоты модули коэффициентов передачи и в совокупности образуют некоторый нелинейный фазовый фильтр с распределенными параметрами. Для его реализации в цифровой форме производится обычная дискретизация переменных ¿1, ¿2 и соответствующих им
частот ю1 , ю2 . Нелинейное преобразование (23) реализуется непосредственно в области определения сигнала, а преобразование Френеля - с помощью одного или двух двумерных БПФ [2,12]. При втором способе затраты времени больше, но зато выше точность.
*=1
3. ОПТИМИЗАЦИЯ ПАРАМЕТРОВ ФИЛЬТРА
Минимальное число пар звеньев, необходимое для реализации предъявленных к оператору Е требований п = 1. С увеличением п отображение Е становится менее критичным к форме ИП.
При п =1 результат отображения Е можно записать в явном виде:
'(t) = So Я exp { + r[2 (Г)]|
(25)
х 2(т) ёт1 йт2.
Для определения функции / в (25), доставляющей максимум показателя селективности (13),
необходимо учесть исходное требование к нелинейному оператору Н: селективное действие на ИП без существенного изменения остальных компонент смеси г(1), которые предполагаются слабыми по сравнению с ИП и (1). Это достигается при условии
f zo (t)]( <<f> (t)]( .
(26)
При этом знаменатель отношения (13) практически не зависит от выбора / и достаточно максимизировать пиковую мощность отклика (25) на ИП , т.е. при г (1) = V (1). Пик (25) в этом случае достигается при 1 = о и равен
)= So Ц exP
ia
(( + l (2
+ i
if[и (t)] >и (t) dt,dt2.
(27)
T
T
Для определения максимума пиковой мощности I v(0) I используем известное неравенство для интеграла от любой комплексной функции Q(t)
ЦQ (t)dt,dt2 <jj(Q (t) (dt,dt2
(28)
которое переходит в равенство при Q(t) = | Q(t) | . В силу него для отдельной реализации ИП и (1) указанный максимум достигается (в предположении, что / - вещественная функция ) при
'(t2 + t2 )
+ f [и (t )] + p(t ) = o
и равен
Um = max (v (0) (2 = (ft (2
JJ(u (t)(dtjdt2
(29)
(3o)
где ф (t) = arg и (t).
Условие (29) есть условие согласования фильтра, имеющего импульсную характеристику (21) с ИП и (t) по фазе.
Выбор фазового фильтра обусловлен, как уже
отмечалось, требованием унитарности оператора Е.
2
В полностью согласованном фильтре величина и
ш
имела бы значение
um=(go (2
jj( (t)(2 dt1 dt2
= (So (2 Eu
(31)
где Е - энергия и (1). Ясно, что различие между (30) и (31) несущественно.
Таким образом, предлагаемое преобразование ИП основано на развитии известных методов сжатия частотно-модулированных импульсов в согласованных фильтрах, но в отличие от них здесь модуляция ставится в зависимость от интенсивности сигнала и фильтр является двумерным.
Как ясно из уравнения (29), оптимальная функция зависит от параметров импульса. Если они случайные, необходимо искать максимум пиковой мощности путем усреднения (27) с учетом их распределений:
P = |v (o) |2 = |So|2 j w (©)fjj exp {^ + if[ и (t, ©)]} и (t, ©)dti dt2 (2 d©,
(32)
где м (©) - плотность вероятности вектора параметров, Я - их область определения.
Рассмотрим параметрическую оптимизацию функционала (32). Пусть функция / определяется заданным заранее аналитическим выражением, содержащим вектор неизвестных параметров р,
/( и ) = /( и, в) .
Тогда точку экстремума функционала (32) можно найти обычным методом, по нулям производных. Приравнивая нулю производную выражения (32) по параметрам Р после несложных преобразований получаем систему уравнений для оптимальных значений Р
j W (©) Re jj jjexp{ [Ф (t1, ©)- Ф (t2, ©)]}и (t1, ©)и * (t2, ©)
Г df [(и (ti, ©), ß(]|
ß
dt, dt, dt. dt© d © = o,
(33)
2
где
* = 1, 2,..., п;
а (( + ¿22) (34)
ф(г, ©) = 1'2 ' + / (и, ©)+ф(г).
Способ решения указанных уравнений зависит от конкретного вида функций/ Пусть, например,
/ ((и (, в)=Ёв*И
тогда
д/ I I*
зрк 1 1
(35)
(36)
Подставляя (36) в (33) и линеаризуя полученные уравнения путем разложения экспоненты в ряд, после ряда преобразований получаем систему линейных уравнений вида
Ё а*в = Ь1-
или, в матричной форме
лр = Ь ,
где
а* = Л Л ехр{га (¿12 + ¿22 -т12 -т22 )}[г'/+*+2 -М,+*+2 }
Т Т
Ь1 = ЯеЦ Цехр(га(¿2 + -т^ -т^)}а,+2 ((, ¿2, т1, т2
Т Т
В выражения (39), (40) входят моментные функции
м+*+2 (1, г2) = | (и (, ©)(+* и (гр ©) и * (г2, ©) * (©) й©,
я
V,+*+2 (гР г2) = I (и (гр ©)(' (и (, ©)(* и (гр ©) и * (г2, ©) * (©) й©,
а,+2 (г1, г2) = | (и (г1, ©)( и (г1, ©)и * (г2, ©)* (©) й©,
При | А | ф 0 уравнения (38) имеют решение в = Л-1Ь. (44)
Поскольку аналитическое вычисление таких интегралов громоздко или затруднено, для расчета коэффициентов (39), (40) целесообразно использовать численные методы.
Алгоритм фильтрации и условия оптимизации фильтра для одномерных (временных) сигналов легко получаются из приведенных выше как частный случай.
4. РЕЗУЛЬТАТЫ СТАТИСТИЧЕСКОГО МОДЕЛИРОВАНИЯ
С целью проверки эффективности разработанного метода подавления ИП было проведено сравнительное статистическое моделирование, в процессе которого воспроизводились реализации сигнала и ИП, осуществлялось подавление ИП рассматриваемым и известным методами и определялись показатели качества обоих методов. При моделировании реализаций ИП задавались характеристики, указанные в п.1, для флуктуационной помехи была принята обычная модель в виде белого гауссовского шума, а для сигнала (изображения) - в виде гауссовского случайного поля с корреляционной функцией экспоненциального вида (воспроизводился фрагмент изображения размером 256x256 отсчетов). Параметры нелинейной фазовой функции, аппроксимированной полиномом второй степени вида (35), были рассчитаны путем численного решения системы уравнений (38) с коэффициентами, найденными по формулам (39) - (43). В качестве не-
(37)
(38)
(39)
(40)
(41)
(42)
(43)
линейного преобразования К ( в области образов Френеля ) выбрано, как наилучшее, бланкирование с экспериментально подобранным оптимальным порогом. Результаты подавления ИП оценивались по двум
2
критериям - СКО (10) д = в и квадрату модуля
2
максимального отклонения (12) дД = Д , а затем сравнивались с аналогичными показателями, полученными путем моделирования одного из наиболее эффективных известных методов подавления ИП на изображениях - медианной фильтрации [3] при тех же характеристиках сигнала и помех.
На рис.2 и 3 показаны полученные зависимости указанных показателей от среднего значения ширины
импульса помехи т по одной из координат для слу-
2
чая, когда дисперсия ширины составляла 10% от т© .
Зависимости для рассматриваемого метода показаны сплошными линиями, а для известного (медианного фильтра) - штриховыми. Размер окна при медианной фильтрации составлял 3x3 отсчета и его увеличение при выбранных параметрах ИП и сигнала не приводило к улучшению качества подавления ИП.
Из приведенных зависимостей видно, что для ИП малого размера по сравнению с рассматриваемой областью изображения, т.е. близких к точечным, предлагаемый метод, как и следовало ожидать, не дает выигрыша или даже имеет проигрыш (по второму критерию), так как нет резервов дополнительного сжатия ИП и более эффективной является обычная медианная фильтрация (или другой извест-
*=0
* =0
ный метод, например, интерполяция по соседним отсчетам). С увеличением размеров ИП все более заметным становится выигрыш, достигаемый предлагаемым методом по сравнению с известным, так как проявляются преимущества сжатия ИП. Например, при т& = 10 достигается выигрыш по СКО около 8 дБ, а по второму критерию - 6,5 дБ. Результаты, полученные при других значениях параметров ИП и сигнала, имеют аналогичный характер.
дех104
40 30 20 10
4 6
Рис.3
В заключение необходимо отметить, существуют резервы дальнейшего повышения эффективности подавления ИП предложенным методом за счет совершенствования вида аппроксимации нелинейной фазовой функции фильтра и увеличения числа его звеньев.
Для более обоснованных рекомендаций по выбору этих и других параметров необходимы дальнейшие эксперименты на реальных изображениях с использованием различных, в том числе субъективных критериев оценки их качества.
ЛИТЕРАТУРА
1. Ярославский Л.П. Введение в цифровую обработку изображений. - М. : Советское радио, 1979.
2. Ярославский Л.П. Цифровая обработка сигналов в оптике и голографии. - М. : Радио и связь, 1987.
3. Быстрые алгоритмы в цифровой обработке изображений / Под ред. Т.С. Хуанга : Пер. с англ. -М. : Радио и связь, 1984.
4. Финк Л.М. Теория передачи дискретных сообщений. - М. : Сов. радио, 1970.
5. Тихонов В.И., Кульман Н.К. Нелинейная фильтрация и квазикогерентный прием сигналов. -М. : Сов. радио, 1975.
6. Кловский Д. Д., Конторович В.Я., Широков С. М. Модели непрерывных каналов связи на основе СДУ. - М. : Радио и связь, 1984.
7. Широков С.М., Григоров И.В. Фильтрация сигналов на фоне импульсных помех с применением нелинейных ортогональных преобразований. Международная конференция и 50-я научная сессия РНТОРЭС им. А.С. Попова. Тезисы докладов, ч.2. -М. : 1995, с. 180.
8. Певницкий В.П., Полозок Ю.В. Статистические характеристики индустриальных радиопомех. -М.: Радио и связь, 1988.
9. Маслов В.П. Комплексные марковские цепи и континуальный интеграл Фейнмана. - М. : Наука, 1976.
10. Ахманов С.А., Выслоух В.А., Чиркин А.С. Оптика фемтосекундных лазерных импульсов. - М.: Наука, 1988.
11. Yevick D., Hermansson B. Solution analysis with the propagating beam method // Opt. Comm. -1983, v. 47, N2, P. 101 - 106.
12. Широков С.М. Различимость импульсов частично когерентного излучения в нелинейном оптическом канале. - Компьютерная оптика, 1993, вып. 13, c. 59 - 64.