Метод подавления импульсных помех при обработке сигналов и изображений c использованием нелинейных фазовых фильтров

С.М. Широков, И.В. Григоров

МЕТОД ПОДАВЛЕНИЯ ИМПУЛЬСНЫХ ПОМЕХ ПРИ ОБРАБОТКЕ СИГНАЛОВ И ИЗОБРАЖЕНИЙ С ИСПОЛЬЗОВАНИЕМ НЕЛИНЕЙНЫХ ФАЗОВЫХ ФИЛЬТРОВ

ВВЕДЕНИЕ

Импульсные помехи (ИП) являются одним из основных факторов, ухудшающих качество изображений, и могут возникать как при передаче сигналов, несущих информацию об изображениях, по каналам связи ( в особенности коммутируемым проводным каналам низкого качества и радиоканалам), так и непосредственно в процессе их обработки (например, вследствие ошибок декодирования) [1-3]. Если ИП носят точечный характер, их влияние удается сравнительно просто свести к минимуму путем "исправления" пораженных помехой элементов изображения с помощью интерполяции по соседним элементам [1] или медианной фильтрации [3]. При обработке одномерных сигналов на практике широко используется ограничение, бланкирование и другие подобные нелинейные операции подавления ИП [4]. Однако все эти методы становятся малоэффективными, если помеха поражает не один, а достаточно много соседних элементов, т.е. становится соизмеримой по размерам с изображением (или по длительности - с сигналом).

Для подавления таких видов ИП в принципе можно использовать методы нелинейной фильтрации по минимуму среднего квадрата отклонения (СКО), основанные на представлении сигналов и изображений марковскими моделями в пространстве состояний [5,6]. Но возможности их практической реализации весьма ограничены - как из-за отсутствия точных решений у соответствующих уравнений для апостериорных распределений и необходимости использования многочисленных и не всегда оправданных приближений, так и потому, что вопрос о применимости марковских моделей для изображений часто остается открытым [6]. Кроме того, критерий минимума СКО далеко не всегда приемлем для оценки качества изображений и других видов непрерывных сообщений, особенно, если речь идет о субъективном восприятии.

В то же время значительный прогресс, достигнутый за последние годы в технике цифровой обработки сигналов, появление процессоров и других видов СБИС с высоким быстродействием значительно расширяют возможности реализации более сложных алгоритмов обработки, чем перечисленные выше простейшие нелинейные преобразования. Один из таких алгоритмов предложен авторами в [7] и включает в себя такие операции, которые позволяют сжать импульсы помехи до точечных практически без изменения полезного сигнала, что дает возможность затем подавить их одним из известных указанных выше методов. Цель данной работы - его обобщение на двумерные сигналы (изображения и поля).

1. МОДЕЛИ ИМПУЛЬСНЫХ ПОМЕХ

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

и (^ ¿2 ) = Е А* х

*

х Ч ( - ¿1* , ¿2 - ¿2* ; Т1к , Т2к ) ^ ((* ) где д (¿1 ; Т1, т2) - детерминированная функция ,

описывающая форму огибающей ИП, А , t , I , т ,

k ^ 2k ^

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

и (I, ©) = Ад (I, т) exp (ф*),

где I = (¿1 , ¿2) е Т - вектор координат из некоторой двумерной области Т, Тп = (т^ т2) и © = ( А, т1, т2, ф ) - векторы случайных параметров. В качестве простого приближения для формы ИП на изображениях можно использовать гауссовскую функцию

д ( т)=exp ^т - т

В каналах связи ИП обычно имеют осцилляции [8], для учета которых можно использовать модель вида

д ( т ) = exp

(

*„ жt

\

Sin

где г - одномерная временная координата, к^ - коэффициент формы, характеризующий затухание ИП. При кф>>0 ИП практически не имеет осцилля-ций, а при кф = 0 осцилляции принимают вид незатухающего гармонического колебания. Результаты статистического анализа реальных радиопомех [6,8] показывают, что для их амплитуд типично логнор-мальное распределение

1 ( 1п2 (А /

•(А) =

Аа\[2я

ехр

2ст2

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

Сехр[-(ги - тт)2];

т ^ т ^ т

w

Т ) =

2с

Т1 < Т2

0;

где С - нормирующая константа, определяемая гра-

2

ницами усечения, т, с -моменты исходной (неусеченной) гауссовской функции.

2. ВЫБОР ОПТИМАЛЬНОГО ПРЕДСТАВЛЕНИЯ СИГНАЛОВ И ПОМЕХ В ЗАДАЧАХ АМПЛИТУДНОЙ СЕЛЕКЦИИ

На вход приемного устройства поступает смесь сигнала 5

2

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

?о

остатка смеси (7) с подавленной ИП

2о

по селектированной функции Уо (&), а оператор W описывает линейный фильтр, обеспечивающий оптимальную оценку полезного сигнала € (1) в смеси (8).

¥ К г 0\" 1 г1 \У

ад

Рис. 1

Критерием оптимальности блоков нелинейной селекции и фильтра в целом при обработке временных сигналов обычно является минимум СКО

= (2о - го (

где

= (5 - € (2

€ = е-1КЕг.

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

д2 =

г, о

С)- )|2

где О - пространство реализаций.

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

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

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

Нетрудно показать, что при использовании ограничения указанное правило сводится к максимизации показателя селективности

7= Ри/Ро:

где Р и Ро - усредненные по реализациям пиковые мощности ИП и остаточной смеси вида (8) в области отображений (функций V (-Э) ).

При использовании бланкирования необходимо минимизировать среднюю ширину ИП в области отображений Зп .

Обоим правилам удовлетворяет отображение Е, обеспечивающее максимальное увеличение амплитуды и уменьшение длительности ИП или площади пораженного ИП участка изображения в смеси (7). Детальный анализ показывает, что такое отображение не может быть линейным и должно зависеть от амплитуды или мгновенной мощности 2(1) (в

2

е

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

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

д

г-+ аДУ + / (у) = 0,

дп

где а - постоянный коэффициент, Д - оператор Лапласа по координатам ^ и ¿2 , /- функция, определяющая вид нелинейности, У(п, г) - комплексная функция, значение которой при п= 0 и некотором п= I соответствует сигналам на входе и выходе блока Е:

у( о, г ) = г (г); у(/, г ) = V (г).

Уравнение (14) задает отображение функций, определенных в области Т , которое описывает преобразование сигнала г (г) в некотором двумерном пространственно распределенном нелинейном фильтре.

При этом обратное отображение задается сопряженным уравнением

дУ

где

т. е.

-г-+ аДУ + / (У) = 0,

у(0, г ) = ^ (г); У(1, г ) = ^ (г),

Е-1 = Е

Непосредственная реализация указанных операторов в аналоговой форме возможна для двумерных и одномерных полей оптического и некоторых других ( например, СВЧ ) диапазонов электромагнитных волн в нелинейных диэлектрических средах. При определенном выборе нелинейности, описываемой функцией /(У) в (14), можно за счет самовоздействия волн обеспечить селективную фокусировку участков поля с аномально большой интенсивностью, обусловленной наложением ИП, и тем самым - увеличение показателя селективности (13). Однако, возможности варьирования и оптимизации параметров нелинейного оператора в этом случае весьма ограничены. Поэтому основное применение рассматриваемый метод может найти при дискретной и цифровой обработке сигналов.

Для такой реализации оператор Е целесообразно представить, используя метод расщепления по физическим факторам [11,12], в виде произведения линейных (Ск) и нелинейных (Нк) операторов вида

Е = П С*Н* [У* ].

* =1

Каждый из линейных операторов Ск реализуется звеном с передаточной функцией

С*

которой соответствует импульсная характеристика

8* (г) = £0* ехР

га*

(¿2 + ) 2

где ю = (  , ) - вектор частот, п - число пар звеньев,

Дп* = П*+1 -п*,

-1

80* = (ДП* аг )2: 1

ак =-.

2Дп* а

Нетрудно заметить, что линейный интегральный оператор с ядром (21) представляет собой известное преобразование Френеля [2] и , таким образом, указанные линейные звенья - это фильтры Френеля. Нелинейные операторы Нк [Ук] в (19) соответствуют умножению функции Ук(т) = У(Пк,т) на зависящие от нее коэффициенты преобразования

Н* [У*]= ехр [/

Сопряженный оператор вследствие уже упоминавшегося свойства унитарности представляется в виде

Е* [У] = П Н* [У * ]С*.

Как видно из (20) и (23), рассматриваемые линейные и нелинейные звенья имеют единичные, не зависящие от частоты модули коэффициентов передачи и в совокупности образуют некоторый нелинейный фазовый фильтр с распределенными параметрами. Для его реализации в цифровой форме производится обычная дискретизация переменных ¿1, ¿2 и соответствующих им

частот ю1 , ю2 . Нелинейное преобразование (23) реализуется непосредственно в области определения сигнала, а преобразование Френеля - с помощью одного или двух двумерных БПФ [2,12]. При втором способе затраты времени больше, но зато выше точность.

*=1

3. ОПТИМИЗАЦИЯ ПАРАМЕТРОВ ФИЛЬТРА

Минимальное число пар звеньев, необходимое для реализации предъявленных к оператору Е требований п = 1. С увеличением п отображение Е становится менее критичным к форме ИП.

При п =1 результат отображения Е можно записать в явном виде:

'(t) = So Я exp { + r[2 (Г)]|

х 2(т) ёт1 йт2.

Для определения функции / в (25), доставляющей максимум показателя селективности (13),

необходимо учесть исходное требование к нелинейному оператору Н: селективное действие на ИП без существенного изменения остальных компонент смеси г(1), которые предполагаются слабыми по сравнению с ИП и (1). Это достигается при условии

f zo (t)]( <<f> (t)]( .

При этом знаменатель отношения (13) практически не зависит от выбора / и достаточно максимизировать пиковую мощность отклика (25) на ИП , т.е. при г (1) = V (1). Пик (25) в этом случае достигается при 1 = о и равен

)= So Ц exP

ia

(( + l (2

+ i

if[и (t)] >и (t) dt,dt2.

T

T

Для определения максимума пиковой мощности I v

ЦQ (t)dt,dt2 <jj(Q (t) (dt,dt2

которое переходит в равенство при Q(t) = | Q(t) | . В силу него для отдельной реализации ИП и (1) указанный максимум достигается (в предположении, что / - вещественная функция ) при

'

+ f [и (t )] + p(t ) = o

и равен

Um = max (v (0) (2 = (ft (2

JJ(u (t)(dtjdt2

где ф (t) = arg и (t).

Условие (29) есть условие согласования фильтра, имеющего импульсную характеристику (21) с ИП и (t) по фазе.

Выбор фазового фильтра обусловлен, как уже

отмечалось, требованием унитарности оператора Е.

2

В полностью согласованном фильтре величина и

ш

имела бы значение

um=(go (2

jj( (t)(2 dt1 dt2

= (So (2 Eu

где Е - энергия и (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) по параметрам Р после несложных преобразований получаем систему уравнений для оптимальных значений Р

j W

Г df [(и (ti, ©), ß(]|

ß

dt, dt, dt. dt© d © = o,

2

где

* = 1, 2,..., п;

а

ф(г, ©) = 1'2 ' + / (и, ©)+ф(г).

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

/ ((и (, в)=Ёв*И

тогда

д/ I I*

зрк 1 1

Подставляя (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 уравнения

Поскольку аналитическое вычисление таких интегралов громоздко или затруднено, для расчета коэффициентов (39), (40) целесообразно использовать численные методы.

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

4. РЕЗУЛЬТАТЫ СТАТИСТИЧЕСКОГО МОДЕЛИРОВАНИЯ

С целью проверки эффективности разработанного метода подавления ИП было проведено сравнительное статистическое моделирование, в процессе которого воспроизводились реализации сигнала и ИП, осуществлялось подавление ИП рассматриваемым и известным методами и определялись показатели качества обоих методов. При моделировании реализаций ИП задавались характеристики, указанные в п.1, для флуктуационной помехи была принята обычная модель в виде белого гауссовского шума, а для сигнала (изображения) - в виде гауссовского случайного поля с корреляционной функцией экспоненциального вида (воспроизводился фрагмент изображения размером 256x256 отсчетов). Параметры нелинейной фазовой функции, аппроксимированной полиномом второй степени вида (35), были рассчитаны путем численного решения системы уравнений (38) с коэффициентами, найденными по формулам (39) - (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.