Synthesis of the reverse filter on the criterion of minimum quadratic error of the weight function mismatch of the path «direct - reverse filter»

Full Text

Введение В настоящее время при решении задач обработки и интерпретации экспери- ментальных данных часто возникает необходимость рассмотрения обратной за- дачи, заключающейся в восстановлении неизвестного входного воздействия по результатам регистрации откликов на выходе средств измерения. В большинстве случаев это задача компенсации искажающего действия аппаратной функции, обеспечивающая улучшение разрешающей способности различного рода изме- рительных приборов и систем [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) пока- зывают, что использование данного подхода позволяет уменьшить погрешность восстановления и аппроксимации, что дает преимущество при решении задач аппроксимации и восстановления данных.

About the authors

Vitaly I Batishchev

Samara State Technical University

(Dr. Sci. (Techn.)), Professor 244, Molodogvardeyskaya st., Samara, 443100, Russian Federation

Igor I Volkov

Samara State Technical University

(Ph.D. (Techn.)), Associate Professor 244, Molodogvardeyskaya st., Samara, 443100, Russian Federation

A G Zolin

Samara State Technical University

Email: zolin.a.g@gmail.com
(Ph.D. (Techn.)), Associate Professor 244, Molodogvardeyskaya st., Samara, 443100, Russian Federation

References

Supplementary files