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

Полный текст

Введение В настоящее время при решении задач обработки и интерпретации экспери- ментальных данных часто возникает необходимость рассмотрения обратной за- дачи, заключающейся в восстановлении неизвестного входного воздействия по результатам регистрации откликов на выходе средств измерения. В большинстве случаев это задача компенсации искажающего действия аппаратной функции, обеспечивающая улучшение разрешающей способности различного рода изме- рительных приборов и систем [1, 2]. В случае, когда для обработки доступна только часть искаженного сигнала, без начальных условий, задача становится недоопределенной и, соответственно, некорректной [3, 4]. Примером таких задач являются задачи восстановления смазанных и расфо- кусированных изображений [5, 6, 7]. Так, в современных космических системах дистанционного зондирования Земли (КСДЗЗ) изображение формируется с по- мощью устройств с зарядовой связью, работающих в режиме временной задерж- ки и накопления оптического сигнала. Для правильной работы таких приборов необходимо, чтобы скорость космического аппарата была точно согласована с периодом опроса светочувствительной матрицы. На практике такое равенство может нарушаться из-за ошибки вычисления скорости спутника [8]. В результате изображение подстилающей поверхности оказывается смазанным вдоль траекто- рии движения летательного аппарата. Конструктивные особенности светочув- ствительных элементов позволяют получить параметры функции рассеяния точ- ки. Способы определения параметров смаза представлены в [9, 10]. Полученные в результате несоответствия скоростей искажения имеют одну пространствен- ную составляющую, что позволяет перейти от двумерной задачи к одномерной и существенно снизить объем вычисляемых данных. Решение двумерной задачи восстановления смазанного изображения представлено, например, в [11]. Постановка задачи Каждая строка полученного с помощью КСДЗЗ смазанного изображения может быть представлена как свертка значений строки исходного изображения х(m) с известной функцией рассеяния точки h0 N0 1 xCM (m)  h0 (i)x(m  i) , (1) i0 здесь N0 - величина весовой функции фильтра - представляет собой величину смаза. Значение h0 для всех строк одинаково. Задача реконструкции полученного с помощью (1) изображения сводится к нахождению функции x’(m), в некотором роде близкой к x(m), по имеющимся значениям xCM(m). Отсутствие начальных значений xCM(m) переводит данную за- дачу в класс обратных некорректно поставленных. Существующие в настоящее время известные подходы к синтезу оптималь- ных алгоритмов реконструкции изображений либо требуют для своей реализации априорной информации, которая не всегда доступна, либо сталкиваются с вы- числительными проблемами, связанными с некорректностью обратных задач и необходимостью использования регуляризующих процедур [12, 13]. В данной статье рассматривается подход, основанный на построении модели весовой функции обратного (восстанавливающего) фильтра (рис. 1). Построение аппроксимационных моделей на базисе стохастических функций при решении задач восстановления сигналов и изображений описано в [14]. В представляемом подходе в качестве критерия адекватности модели ис- пользуется минимум квадратической погрешности рассогласования весовой функции тракта «прямой - обратный фильтр». Синтез обратных КИХ-фильтров на основе критерия моментов представлен в [15]. Решением задачи восстановления изображения будем считать нахождение функции h(), представляющей собой весовую функцию обратного фильтра, поз- воляющего получить с помощью операции свертки оценку восстановленного изображения х’(m): N 1 x'(m)  h(i)xCM (m  i) . (2) i0 h h0 Прямой (искажающий) фильтр Обратный (восстанавлива- ющий) фильтр Исходный сигнал Искаженный сигнал Восстановленный сигнал Весовая функция тракта H Рис. 1. Тракт «прямой - обратный фильтр» Рассмотрим случай, когда прямой и обратный фильтры являются КИХ- фильтрами. Весовыми функциями этих фильтров, соответственно, будут h0(i), i=0,…,N0-1 и h(i), i=0,…,N-1. Тогда весовая функция тракта «прямой - обратный фильтр» будет равна i H (i)  h0 (v)h(i  v) . (3) v0 i  0,...,N  N0  2 В идеальном случае должно быть выполнено условие H (i)  (i) где δ(i) - символ Кронекера. i  0,...,N  N0  2 , (4) Это условие принципиально выполнено быть не может, потому что обрат- ный фильтр должен быть БИХ-фильтром. Поэтому всегда будет иметь место по- грешность, квадратичное значение которой будет равно N  N0 1    i 0 (H (i)  (i))2  (H (0) 1)2  N  N0 1 H (i)2 . (5) i 1 При заданном прямом фильтре оно будет зависеть от весовой функции h(i) обратного фильтра и от N. Связь погрешности восстановления сигнала с квадратичной погрешностью Неискаженный сигнал x(m) и его восстановленное значение x’(m) связаны между собой соотношением x'(m)  N  N0 2 H (i)x(m  i) . (6) i0 Погрешность восстановления равна (m)  x'(m)  x(m) , или с учетом (6) (m)  N  N0 2 N  N0 2 H (i)x(m  i)  x(m)  i0  H (i)x(m  i)  (H (0) 1)x(m)  i1 (7) где N  N0 2  i x(m  i) i0   H (0) 1 при i  0 .  i  H (i) при i  1 В соответствии с неравенством Коши - Буняковского из (7) находим, что или 2 (m)  N  N0 2  i 0 2 (i)  N  N0 2  x2 (m  i) , i 0   N  N0 2 N  N0 2 2 (m)   (H (0) 1)2   H 2 (i)  x2 (m  i) .  Принимая во внимание (5), получим i 1  i 0 2 (m)   N  N0 2  x2 (m  i) . (8) i 0 Отсюда следует, что значение ε должно быть как можно меньшим. Синтез обратного фильтра по минимуму квадратичной погрешности Найдем решение данной задачи в частотной области. Введем в рассмотрение частотные характеристики фильтров   W0 ( j)     W ( j)    N  N0 2 N0 1 h0 (i) exp( ji) i0 N0 1 h(i) exp( ji) i0 (9) H ( j)   H (i) exp( ji)  W0 ( j)W ( j) i0 С учетом этих соотношений квадратичная погрешность (5) примет вид     2  W ( j)W0 ( j) 1W ( j)W0 ( j) 1d    (10) Для обеспечения минимума этой погрешности значения весовой функции обратного фильтра должны исходить из условия d dh(k) Из (8) с учетом (9) находим  0, k  0,...,N 1 . (11)     d  N 1   2     2 h(v)  W0 ( j) cos((k  v))d  W0 ( j)exp( jk)d dh(k)   v0  2    2       (12) С учетом этого соотношения из (10) находим, что min  1 H (0) . (13) Справедливы соотношения:      2 N0 1  2  W0 ( j) cos(k)d  h0 (v)h0 (v  k)          vk  W0 ( j) exp( jk)d  h0 (k)  2      С учетом этих соотношений формула (12) примет вид d  N 1   2 h(v)E0 ( k  v )  h0 (k), где dh(k)    v0  N0 1 E0 (k)  h0 (v)h0 (v  k), k  0,...,N0 1 . (14) vk В итоге для определения значений h(k) необходимо решить систему уравне- ний: N 1 h(v)E0 ( k  v )  h0 (k), k  0,...,N 1 . (15) v0 Алгоритм решения этой системы уравнений таков:  g(0,0)  1  (0)  E0 (0)  1 k p1   g(k,0)   (k 1)  g(k 1, p 1)E0 ( p), k  1,...,N0  2  N 1  1 0 (16) g(k,0)   (k 1)  g(k 1, p 1)E0 ( p), k  N0 1,...,N 1  p1 g(k,v)  g(k 1,v 1)  g(k 1, k 1  v)g(k,0), v  1,...,k 1  g(k, k)  1 (k )  (k 1)(1  g 2 (k,0), k  1,...,N 1 N 1 g(q, k)g(q,0) h(k)  h0 (0)  qk (q) , k  0,...,N 1 (17) Вычисляя значение минимальной погрешности из (3), имеем: H (0)  h(0)h0 (0) . Подставив сюда h(0) из (17), получаем: 2 N 1 g 2 (q,0) H (0)  h0 (0)  q0 (q) . Тогда в соответствии с (13) минимальное значение погрешности будет равно 2 N 1 g 2 (q,0) min  1  h0 (0)  q0 Для справки приведем еще такие формулы: (q) . (18) N  N0 2 2 N 1N 1    k 0 H (k )   h(i)h( j)E0 ( i  j ) i0 j0  N  N0 2  d  H 2 (k )   k 0   dh(v) N 1 2 h( j)E0 ( v  j ) j0 (19)  N 1N 1    1  2H (0)   h(i)h( j)E0 ( i  j )  i0 j0    Синтез обратного фильтра по минимуму квадратичной погрешности при внесении ряда ограничений, наложенных на его весовую функцию Исходя из принципов физической реализуемости весовая функция синтезируемого фильтра должна удовлетворять условиям N 1 (i,v)h(i)  (v),v  0,..., p . (20) i0 Так как функции H(i) и h(i) однозначно взаимосвязаны, то условия, нало- женные на функцию H(i), могут быть сведены к условиям, наложенным на h(i). Например, если функция H(i) должна удовлетворять условиям N  N0 2 (i,v)H (i)  (v),v  0,..., p i0 то эти условия могут быть записаны в виде (20), где N0 1 , (21)  (i, v)   h0 (k) (i  k, v), v  0,..., p . (22) k 0 Задачу синтеза обратного фильтра по минимуму квадратичной погрешности с учетом ограничений физической реализуемости, наложенных на его весовую функцию, будем решать методом Лагранжа. Для этого рассмотрим функцию p N 1 f (h(0),...,h(N 1))    2 Av (i,v)h(i) , где Av - коэффициенты Лагранжа. v0 i0 Значения h(i) будем определять из условий минимума функции f(h(0),…,h(N-1)): df (h(0),...,h(N 1))  0, k  0,...,N 1. dh(k) В результате на основании рассмотренных ранее выводов для обеспечения этих условий получим следующую систему уравнений: N 1 N h(i)E0 ( k  i )  h0 (k)   Av(k,v),k  0,...,N 1. i0 v0 Решив ее, получим N 1 g(q,i)  p  h(i)    g(q,0)h0 (0)   Ak D(q, k) . (23) Здесь qi (q)    k 0  q D(q, k)  g(q,v)(v, k) ; (24) v0 g(q,k), Ψ(q) - рассмотренные ранее величины. Подставив h(i) из (23) в (20), получим следующую систему уравнений для определения значений коэффициентов Лагранжа: p N 1 D(q,v)D(q, k) N 1D(q,v)  Ak  (q)   (q) g(q,0)h0 (0)  (v),v  0,..., p . k 0 q0 q0 Для простоты дальнейших выкладок введем обозначения: N 1 D(q, v)D(q, k) f (v, k)   q0 N 1 D(q, v) (q) (25) F (v)   q0 (q) g(q,0)h0 (0)  (v) В этом случае система уравнений примет вид p  Ak f (v, k)  F (v) . (26) k 0 Алгоритм решения этой системы уравнений будет таким. Сначала по форму- лам (24) и (25) определим значения f(v,k), v=0,…,p; k=0,…,v и F(v), v=0,…,p. За- тем выполним вычисления по формулам:  k 1 E(v, k)  f (v, k)  C(k, q)E(v, q)   C(v, k)   E(v, k) E(k, k) q0 (27) k  0,...,v; v  0,..., p     v1 (v)  F (v)  C(v, k)(k)   v  0,..., p k 0 A  (k)  k E(k, k) p  AvC(v, k),k  p,...,0 vk 1 Подставив найденные значения Ak в (23), получим искомые значения весовой функции обратного фильтра. Далее по формуле (2) получим значение восстанов- ленного изображения х’(m). Апробация результатов Для апробации алгоритма было взято изображение, полученное в процессе дистанционного зондирования Земли. Благодаря конструктивным особенностям регистрирующего аппаратуры изображение имеет 1024 градации серого, что в 4 раза выше значений, принятых в распространенных форматах хранения графиче- ских файлов на персональных компьютерах. Из тестового изображения был взят тестовый фрагмент 480 на 285 пикселей, над которым и проводились экспери- менты (рис. 2, а). Каждая строка тестового изображения была обработана фильтром (1) с весо- вой функцией: 1 h0  N . 0 Таким образом, был выполнен смаз изображения вдоль горизонтальной оси на различное количество пикселей N0 = 3, 4, 5, …, 10. Пример полученного сма- занного изображения приведен на рис. 2, б. Рис. 2. Эталонное изображение Рис. 3. Изображение, смазанное на 10 пикселей Далее были предприняты попытки восстановления по алгоритмам (16) и (27). Для количественной оценки качества восстановления использовалась от- носительная среднеквадратическая погрешность (ОСП), вычисляемая по форму- ле ОСП  M 1  (x '( j)  x( j))2 , j 0  M 1 x( j)2 j 0 (28) где х(j) - значение пикселя строки эталонного изображения; х’(j) - значение пикселя строки восстановленного изображения; M - количество пикселей в строке. Данный способ оценки возможен благодаря наличию и эталонного, и вос- становленного изображения. Полученная величина показывает степень отклоне- ния результата от исходного изображения. Другие способы оценки качества вос- становления изображений приведены, например, в [16]. Значения ОСП смазанных изображений без восстановления приведены в табл. 1. Таблица 1 ОСП невосстановленных изображений N0 3 4 5 6 7 8 9 10 ОСП 0,0505 0,0565 0,0672 0,0743 0,0833 0,0913 0,0974 0,1041 Восстановление изображений по алгоритму (16) и (17) без наложенных ограничений. Анализ формул (14) - (19) показывает, что результат восстановле- ния зависит от значений весовых функций прямого и обратного фильтров: N и N0. Некоторые значения погрешности восстановления приведены в табл. 2. Таблица 2 Результаты восстановления по алгоритму (16), 17) N0 N ОСП N0 N ОСП N0 N ОСП 3 3 0,044 5 5 0,0617 7 7 0,062 3 5 0,043 5 9 0,0494 7 13 0,056 3 7 0,043 5 13 0,0602 7 19 0,060 4 4 0,061 6 6 0,0593 8 8 0,075 4 7 0,053 6 11 0,0543 8 15 0,066 4 10 0,060 6 16 0,0539 8 22 0,077 Графики зависимости ОСП от N0 показаны на рис. 4. При увеличении длины весовой функции обратного фильтра N процесс восстановления становится более неустойчивым (рис. 5), что не позволяет значительно увеличивать N. Рис. 4. Зависимость ОСП от N0 и N: 1 - значение ОСП при N0 = 8; 2 - значение ОСП при N0 = 4; 3 - значение ОСП при N0 = 5; 4 - значение ОСП при N0 = 7; 5 - значение ОСП при N0 = 6; 6 - значение ОСП при N0 = 3 Рис. 5. Фрагмент строки эталонного и восста- новленного изображений при N0 = 3: x(m), x’(m) - значения пикселей; m - номер пикселя; 1 - восстановленное изображение при N = 13; 2 - эталонное изображение Восстановление изображений по алгоритму (27) c наложенными ограниче- ниями. Для выполнения условий (21) и (22) были выбраны следующие соотно- шения: H (0)  h(0)h0 (0) 1 ; N  N0 1 H (i)iv  (v),v  0,..., p 1; i0  (i, v)   h (k)(i  k) . N0 1 v 0 k 0 Анализ алгоритма построения весовой функции обратного фильтра (20) - (27) показывает, что результат восстановления зависит от значений весовых функций прямого и обратного фильтров: N и N0, а также p из (22). При проведе- нии экспериментов с различными значениями N0, N и p наименьшие значения ОСП получались при N= 2*p-1. Зависимость ОСП от значений N и p показана на рис. 6. Рис. 6. Зависимость ОСП от N: 1 - N0 = 6; 2 - N0 = 5; 3 - N0 = 4; 4 - N0 = 3 Рис. 7. Зависимость ОСП от порядка модели с различными N0: 1 - N0 = 8; 2 - N0 = 5; 3 - N0 = 4; 4 - N0 = 7; 5 - N0 = 6; 6 - N0 = 3 Некоторые полученные результаты представлены в табл. 3. Таблица 3 Результаты восстановления по алгоритму (27) N0 P ОСП N0 P ОСП N0 P ОСП 3 2 0,034 5 2 0,044 7 2 0,055 3 3 0,044 5 3 0,046 7 3 0,068 3 4 0,054 5 4 0,065 7 4 0,067 4 2 0,029 6 2 0,046 8 2 0,063 4 3 0,042 6 3 0,062 8 3 0,074 4 4 0,065 6 4 0,069 8 4 0,068 Выводы На основе анализа известных подходов к синтезу оптимальных алгоритмов реконструкции сигналов, основанных на использовании регуляризующих проце- дур при решении некорректных обратных задач, и связанных с этим вычисли- тельных проблем предложен метод построения цифровых фильтров для решения обратных задач восстановления сигналов, временных рядов и изображений с использованием аппроксимационного подхода. Сформулирована постановка задачи восстановления смазанных и расфоку- сированных изображений в космических системах дистанционного зондирования Земли, показана возможность снижения объема вычислений при обработке дву- мерных массивов данных путем перехода к решению одномерной задачи. Разработаны и исследованы алгоритмы синтеза обратного фильтра для вос- становления сигналов с известной весовой функцией прямого фильтра на основе критерия минимума квадратической погрешности рассогласования весовой функции всего тракта «прямой - обратный фильтр» для реконструкции смазан- ных изображений. Проведена апробация алгоритмов на модельных примерах и при обработке реальных изображений, полученных при дистанционном зондировании Земли. Для количественной оценки качества восстановления использовалась относи- тельная среднеквадратическая мера различия эталонного и восстановленного сигналов (изображений). Приведенные результаты апробации (табл. 1, 2, 3) пока- зывают, что использование данного подхода позволяет уменьшить погрешность восстановления и аппроксимации, что дает преимущество при решении задач аппроксимации и восстановления данных.

Об авторах

Виталий Иванович Батищев

Самарский государственный технический университет

(д.т.н., проф.), профессор кафедры «Информационные технологии» Россия, 443100, г. Самара, ул. Молодогвардейская, 244

Игорь Иванович Волков

Самарский государственный технический университет

(к.т.н., доц.), доцент кафедры «Информационные технологии» Россия, 443100, г. Самара, ул. Молодогвардейская, 244

Алексей Георгиевич Золин

Самарский государственный технический университет

Email: zolin.a.g@gmail.com
(к.т.н., доц.), доцент кафедры «Информационные технологии» Россия, 443100, г. Самара, ул. Молодогвардейская, 244

Список литературы

Дополнительные файлы