Research algorithm restoring signals in the basis of exponential functions

Full Text

Введение

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

В большинстве случаев это задача компенсации искажающего действия аппаратной функции, обеспечивающая улучшение разрешающей способности различного рода измерительных приборов и систем [1, 2]. В случае, когда для обработки доступна только часть искаженного сигнала, без начальных условий, задача становится недоопределенной и, соответственно, некорректно поставленной [3, 4].

Решение таких задач требует методов регуляризации, базирующихся на привлечении априорной информации о решении, которая может быть как качественной (неотрицательность и гладкость решения), так и количественной [5, 6, 7]. Общим фундаментальным свойством методов регуляризации является их ориентация на принципиально смещенные решения. Это же свойство является фундаментальным для аппроксимационных методов, имеющих явные перспективы в решении различного рода обратных задач [8, 9], в том числе и задач восстановления сигналов [10, 11, 12]. На основе обобщенных принципов получения информации об исследуемых объектах экспериментальным путем по результатам измерений аппроксимационный подход к восстановлению сигналов позволяет использовать аналитические модели функциональных характеристик этих объектов, выбранных на основе априорной информации с учетом целей получения того или иного результата.

Постановка задачи

В общем случае искаженный сигнал х_см(m) может быть представлен как свертка значений исходного сигнала х_исх(m) с известной весовой функцией прямого (искажающего) фильтра h₀:

$x_{см}(m)={\displaystyle \sum _{i=0}^{N_0-1}{h}_0(i)x_{исх}(m-i)},$ (1)

здесь N₀ $–$ величина весовой функции фильтра, представляет собой количество значений h₀.

Задача реконструкции полученного с помощью (1) сигнала сводится к нахождению функции x_вст(m), в некотором роде близкой к x_исх(m), по имеющимся значениям x_см(m). Данная задача представляет собой обратную задачу. Отсутствие начальных значений x_см(m) переводит данную задачу в класс некорректно поставленных.

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

$y=\exp (-a\cdot x).$

Применение других видов базисных функций в задачах восстановления сигналов (полиномиальных, тригонометрических) рассмотрено, например, в [13]. Построение аппроксимационных моделей на базисе стохастических функций при решении задач восстановления сигналов и изображений описано в [12].

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

Также на практике используются и другие критерии адекватности; так, например, синтез обратных фильтров на основе критерия моментов представлен в [13], а синтез нелинейных адаптивных фильтров для решения обратных задач восстановления сигналов $–$ в [14].

Решением задачи восстановления сигнала будем считать нахождение функции h(i), представляющей собой весовую функцию обратного фильтра, позволяющего получить с помощью операции свертки оценку восстановленного сигнала х_вст(m):

$x_{вст}(m)={\displaystyle \sum _{i=0}^{N-1}h(i)x_{см}(m-i)}.$ (2)

Представим оценку восстановленного сигнала в виде следующей модели:

$x_{вст}^{м}={\displaystyle \sum _{v=0}^p{C}_v\exp (-a\cdot m(v+1))},$ (3)

где p $–$ порядок модели;

C_v $–$ параметры модели;

а $–$ коэффициент, связанный с шагом дискретизации.

Подставив x^м_вст(m) в (2), получим выходной сигнал прямого фильтра в виде

$x_{см}^{м}={\displaystyle \sum _{v=0}^p{C}_v\mu (v)\exp (-am(v+1)),}$ (4)

где

$\mu (v)={\displaystyle \sum _{k=0}^{N_0-1}{h}_0(k)\exp (ak(v+1)).}$

Связь погрешности восстановления сигнала с квадратичной погрешностью

Значения Cv будем определять по выборке {y_i}, i=0, …, N-1 на основе обеспечения минимума квадратической погрешности:

$E={\displaystyle \sum _{i=m-N+1}^m{(x_{см}^{м}(i)-x_{см}(i))}^2={\displaystyle \sum _{i=0}^{N-1}{(x_{см}^{м}(m-i)-x_{см}(m-i))}^2}}.$ (5)

Значение погрешности Е зависит от величин С₀, …, С_р.

Для обеспечения минимума этого значения должно быть выполнено условие

$\frac{dE}{d{C}_k}=\mathrm{0,}k=\overrightarrow{\mathrm{0,}p}.$

С учетом (3) это условие примет вид

${\displaystyle \sum _{i=0}^{N-1}(x_{см}^{м}(m-i)-x_{см}}(m-i))\frac{d{x}_{см}^{м}(m-i)}{d{C}_k}=\mathrm{0,}k=\overrightarrow{\mathrm{0,}p}.$ (6)

$$ Подставив сюда x^м_см(m-i) из (4), получим следующее выражение:

$\begin{array}{c}{\displaystyle \sum _{v=0}^p{C}_v\mu (v){\displaystyle \sum _{i=0}^{N-1}\exp (-a(v+k+2)(m-i))=}}\\ ={\displaystyle \sum _{i=0}^{N-1}{x}_{см}(m-i)\exp (-a(k+i)(m-i))},k=\overrightarrow{\mathrm{0,}p.}\end{array}$ (7)

Левую и правую части системы (7) умножим на exp(am(k+1)) и приведем ее к виду

$\begin{array}{c}{\displaystyle \sum _{v=0}^p{C}_v\mu (v)\exp (-am(v+1))B(k+v)=}\\ ={\displaystyle \sum _{i=0}^{N-1}{x}_{см}\exp (ai(k+1)),k=\overrightarrow{\mathrm{0,}p,}}\end{array}$ (8)

где

$B(n)=\frac{1-\exp (aN(n+2))}{1-\exp (a(n+2))},n=\overrightarrow{\mathrm{0,2}p}.$ (9)

Введем следующее обозначение:

$A_v(m)=C_v\exp (-am(v+2)).$ (10)

С учетом этого модель (3) будет выглядеть так:

$x_{вст}^{м}(m)={\displaystyle \sum _{v=0}^p{A}_v(m),}$ (11)

а система уравнений (8) примет вид

$\begin{array}{c}{\displaystyle \sum _{v=0}^p\mu (v)A_v(m)B(k+v)=}\\ ={\displaystyle \sum _{i=0}^{N-1}{x}_{см}(m-i)\exp (ai(k+1)),k=\overrightarrow{\mathrm{0,}p}.}\end{array}$ (12)

 

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

Обозначим элементы матрицы, обратной матрице коэффициентов B(k+v) системы (12), как β(n₁,n₂). Тогда решение системы уравнений будет таким:

$A_v(m)=\frac{1}{\mu (v)}{\displaystyle \sum _{i=0}^{N-1}{x}_{см}(m-i){\displaystyle \sum _{k=0}^p\beta (v,k)\exp (ai(k+1)),v=\overrightarrow{\mathrm{0,}p}.}}$ (13)

Подставив A_v(m) из (13) в (10), получим формулу восстановления исходного сигнала:

$x_{вст}^{м}(m)={\displaystyle \sum _{i=0}^{N-1}h(i)x_{см}(m-i),}$ (14)

где

$h(i)={\displaystyle \sum _{v=0}^p\frac{1}{\mu (v)}{\displaystyle \sum _{k=0}^p\beta (v,k)\exp (ai(k+1)).}}$ (15)

Здесь h(i) $–$ весовая функция обратного (восстанавливающего) фильтра.

Представленный способ построения весовой функции обратного фильтра требует определить элементы обратной матрицы коэффициентов β(n₁,n₂). Для выполнения этой процедуры был разработан более простой, с точки зрения реализации и объемов вычислений, алгоритм. Найдем коэффициенты K(n₁,n₂), которые вычислим по следующей формуле:

$\left\{\begin{array}{c}\Psi (k,v)=B(k+v)-{\displaystyle \sum _{q=0}^{v-1}K(v,q)\Psi (k,q)}\\ K(k,v)=\frac{\Psi (k,v)}{\Psi (v,v)}\\ v=\overrightarrow{\mathrm{0,}k}\\ k=\overrightarrow{\mathrm{0,}p}\end{array}\right.$ (16)

Используя полученные коэффициенты, получим выражение для нахождения весовой функции обратного фильтра:

$\left\{\begin{array}{c}g(m,m)=1\\ g(m,q)=-{\displaystyle \sum _{k=q+1}^{m}g(m,k)K(k,q)}\\ q=\overrightarrow{m-\mathrm{1,0}}\\ h(i)={\displaystyle \sum _{k=0}^p\exp (ai(k+1)){\displaystyle \sum _{q=k}^p\frac{g(q,k)}{\Psi (q,q)}{\displaystyle \sum _{v=0}^q\frac{g(q,v)}{\mu (v)}}}}\end{array}\right.$ (17)

 

Апробация результатов

Одним классом задач восстановления сигналов являются задачи восстановления смазанных и расфокусированных изображений [15, 16, 17]. Так, в системах дистанционного зондирования Земли, используемых в том числе в космических исследованиях, изображение формируется с помощью устройств с зарядовой связью, работающих в режиме временной задержки и накопления оптического сигнала. Эксплуатация таких приборов требует, чтобы скорость космического аппарата была точно согласована с периодом опроса светочувствительной матрицы. На практике такое равенство может нарушаться из-за погрешностей вычисления скорости спутника [18]. В результате изображение подстилающей поверхности оказывается смазанным вдоль траектории движения летательного аппарата. Конструктивные особенности светочувствительных элементов позволяют получить параметры функции рассеяния точки, являющейся в нашем случае весовой функцией прямого фильтра. Способы определения параметров смаза представлены в [19, 20]. Полученные в результате несоответствия скоростей искажения имеют одну пространственную составляющую, что позволяет прейти от двумерной задачи к одномерной и существенно снизить объем вычисляемых данных. Решение двумерной задачи восстановления смазанного изображения представлено, например, в [21].

Апробация алгоритма проводилась на двух типах изображений: на изображениях, полученных в процессе дистанционного зондирования Земли, и на специально созданных контрастных изображениях для проверки поведения алгоритмов в условиях резкого перепада яркостей. В случае с первым типом изображений из-за конструктивных особенностей регистрирующего аппаратуры полученные данные имеют 1024 градации серого, что в 4 раза выше значений, принятых в распространенных форматах хранения графических файлов на персональных компьютерах. Для согласованности результатов измерений контрастные изображения были сформированы в том же формате. Пример тестового эталонного изображения, полученного в результате дистанционного зондирования Земли, показан на рис. 1, пример контрастного изображения приведен на рис. 2.

 

Рис. 1. Эталонное изображение участка Земли

 

Рис. 2. Контрастное изображение

 

Каждая строка тестового изображения была обработана по формуле (1) фильтрами с весовыми функциями:

$h_0(i)=\frac{1}{N_0},i=\overrightarrow{\mathrm{0,}{N}_0-1}$; (18)

$h_0(i)=\frac{ki+b}{N_0(b+\mathrm{0,5}k(N_0+1))-k},i=\overrightarrow{\mathrm{0,}{N}_0-1}$; (19)

$h_0(i)=\frac{-ki+b}{N_0(b-\mathrm{0,5}k(N_0+1))+k},i=\overrightarrow{\mathrm{0,}{N}_0-1}$. (20)

Схематичное изображение весовых функций из формул (18), (19), (20) приведено на рис. 3.

 

Рис. 3. Схематичное изображение весовых функций из формул (18), (19), (20)

 

Таким образом, было выполнено размытие изображения вдоль горизонтальной оси на различное количество пикселей N₀ = 3, 4, 5, …, 10. Для весовых функций (19) и (20) брались следующие значения: k=0.2, b = 1.0.

В ходе проведения экспериментов были предприняты попытки восстановления по алгоритмам (16), (17) и (14). Для количественной оценки качества восстановления использовалась относительная среднеквадратическая погрешность (ОСП), вычисляемая по формуле

$ОСП=\sqrt{\frac{{\displaystyle \sum _{j=0}^{M-1}{(x_{вст}^{м}(j)-x_{исх}(j))}^2}}{{\displaystyle \sum _{j=0}^{M-1}{x}_{ucх}{(j)}^2}}},$ (21)

где х_исх(j) $–$ значение пикселя строки эталонного изображения;

х^м_исх(j) $–$ значение пикселя строки восстановленного изображения;

M $–$ количество пикселей в строке.

Поскольку в процессе апробации одновременно есть доступ и к эталонному, и к восстановленному изображению, данный способ позволяет оценить качество исследуемых алгоритмов. Полученная величина показывает степень отклонения результата от исходного изображения. Другие способы оценки качества восстановления изображений приведены, например, в [22].

Значения ОСП смазанных изображений земной поверхности без восстановления приведены в табл. 1.

 

Таблица 1 ОСП невосстановленных изображений

N0 ОСП	3	4	5	6	7	8	9	10
h0 (18)	0,0328	0,0512	0,0694	0,0855	0,0983	0,1082	0,1155	0,121
h0 (19)	0,0334	0,0448	0,0584	0,0710	0,0813	0,0895	0,096	0,101
h0 (20)	0,0336	0,0587	0,0812	0,1007	0,1161	0,1276	0,1361	0,1421

 

Анализ представленного в работе алгоритма показывает, что результат восстановления зависит от значений весовых функций прямого и обратного фильтров N , N₀ и порядка модели P. Некоторые значения погрешности восстановления изображений земной поверхности с разными видами весовых функций прямого фильтра приведены в табл. 2, 3, 4. Минимальные значения погрешности выделены жирным шрифтом.

 

Таблица 2 Значение ОСП восстановления в случае весовой функции (18)

N0=5
N P	3	4	5	6	7	8	9	10	11	12
16	0.038	0.058	0.057	0.056	0.055	0.051	0.046	0.035	0.031	0.030
17	0.038	0.058	0.057	0.057	0.368	0.068	0.046	0.032	0.031	0.030
18	0.038	0.057	0.057	0.057	0.066	0.057	0.046	0.036	0.031	0.030
19	0.038	0.057	0.057	0.053	0.051	0.057	0.045	0.035	0.102	0.032
20	0.038	0.057	0.057	0.056	0.055	0.057	0.049	0.070	0.116	0.029
N0=6
N P	4	5	6	7	8	9	10	11	12	13
16	0.064	0.063	0.060	0.057	0.054	0.060	0.049	0.040	0.041	0.039
17	0.063	0.063	0.068	0.178	0.058	0.060	0.050	0.040	0.040	0.039
18	0.063	0.063	0.061	0.053	0.055	0.060	0.050	0.040	0.041	0.103
19	0.063	0.063	0.062	0.051	0.055	0.060	0.052	0.115	0.051	0.038
20	0.063	0.063	0.06	0.052	0.055	0.058	0.070	0.129	0.045	0.038
N0=7
N P	7	8	9	10	11	12	13	14	15	16
10	0.052	0.048	0.047	0.045	0.061	0.046	0.046	0.050	0.079	0.068
11	0.051	0.048	0.047	0.052	0.053	0.047	0.046	0.053	0.065	0.273
12	0.052	0.057	0.056	0.065	0.053	0.046	0.043	0.047	0.043	0.067

 

Таблица 4 Значение ОСП восстановления в случае весовой функции (20)

N0=5
N P	3	4	5	6	7	8	9	10	11	12
27	0.252	0.244	0.183	0.118	0.066	–	0.120	0.083	0.132	0.223
28	0.885	0.130	0.161	0.122	0.080	–	0.171	0.101	0.139	0.258
29	0.132	0.359	0.165	0.132	0.030	–	0.043	0.084	0.132	0.214
30	0.128	0.436	0.185	0.093	0.078	–	0.043	0.076	0.071	0.268
31	0.083	0.263	0.134	0.093	0.050	–	0.088	0.082	0.143	0.261
N0=6
N P	3	4	5	6	7	8	9	10	11	12
32	0.080	0.241	0.392	0.141	0.144	–	0.138	0.182	0.260	0.209
33	0.068	0.251	0.204	0.139	0.099	–	0.113	0.137	0.261	0.178
34	0.067	0.267	0.173	0.145	0.151	–	0.117	0.155	0.244	0.304
35	0.510	0.419	0.219	0.157	0.155	–	0.111	0.123	0.237	0.309
36	0.498	0.318	0.298	0.170	0.140	–	0.111	0.130	0.232	0.274
N0=7
N P	11	12	13	14	15	16	17	18	19	20
14	0.160	0.198	0.276	0.354	0.430	0.397	0.578	0.633	0.722	0.765
15	0.261	0.225	0.271	0.355	0.463	0.079	0.588	0.724	0.704	0.763
16	0.165	0.200	0.264	0.356	0.432	0.074	0.484	0.653	0.708	0.765
17	0.154	0.197	0.294	0.354	0.147	0.518	0.451	0.628	0.705	0.765
18	0.198	0.202	0.294	0.352	0.399	0.491	0.527	0.628	0.712	0.810

 

Графики зависимости ОСП от N при различных N₀ и Р в случае минимальной погрешности показаны на рис. 4, а график зависимости от Р $–$ на рис. 5. Эти графики приведены для случая весовой функции прямого фильтра (18).

 

Рис. 4. Зависимость ОСП от N при заданных N0 и Р: 1 N0 =3, P=14; 2 N0 =4, P=19; 3 N0 =5, P=20; 4 N0 =6, P=20; 5 N0 =7, P=12; 6 N0 =8, P=4; 7 N0 =9, P=2

 

Рис. 5. Зависимость ОСП от P при заданных N0 и N: 1 - N0 =3, N=7; 2 - N0 =4, N=9; 3 - N0 =5, N=12; 4 - N0 =6, N=13; 5 - N0 =7, N=13; 6 - N0 =8, N=8; 7 - N0 =9, N=4.

 

На рис. 6 представлен график зависимости ОСП от N при различных N₀ и Р в случае минимальной погрешности для случая весовой функции прямого фильтра (18).

 

Рис. 6. Зависимость ОСП от N: 1 N0=3, P=27; 2 N0=4, P=33; 3 N0=5, P=40; 4 N0=6, P=38;

 

Результаты восстановления контрастного изображения (см. рис. 2) для случая весовой функции прямого фильтра (18) представлены в табл. 5. Однако для получения представления работы алгоритма на контрастном изображении рассмотрим вид полученного восстановленного сигнала (рис. 7). Из рисунка видно, что большой вклад в погрешность вносят выбросы на границах перепада яркостей, но далее алгоритм формирует сигнал с низким значением ОСП. Погрешность восстановления в этом случае можно несколько уменьшить; учитывая, что максимальные и минимальные значения пикселя изображения не должны выходить из интервала [0; 1023], явно заменим выходящие за диапазон величины на 0 и 1023.

 

Рис. 7. Фрагмент восстановленной строки контрастного изображения

 

Таблица 5 Значение ОСП восстановления контрастного изображения

N0=5
N P	5	6	7	8	9	10	11	12	13	14
24	0.095	0.118	0.146	0.072	0.099	–	0.112	0.128	0.181	0.162
25	0.092	0.118	0.072	0.070	0.094	–	0.158	0.115	0.170	0.156
26	0.092	0.118	0.129	0.068	0.535	–	0.107	0.132	0.107	0.162
27	0.095	0.118	0.129	0.066	0.541	–	0.107	0.136	0.134	0.157
28	0.088	0.118	0.130	0.069	0.457	–	0.110	0.134	0.156	0.155
N0=6
N P	4	5	6	7	8	9	10	11	12	13
31	0.150	0.150	0.160	0.140	0.189	–	0.083	0.108	0.092	0.108
32	0.133	0.133	0.162	0.139	0.226	–	0.082	0.110	0.110	0.097
33	0.136	0.136	0.163	0.139	0.225	–	0.082	0.106	0.155	0.097
34	0.136	0.136	0.168	0.139	0.196	–	0.083	0.129	0.091	0.111
35	0.137	0.137	0.164	0.139	0.120	–	0.082	0.117	0.111	0.098
N0=7
N P	8	9	10	11	12	13	14	15	16	17
14	0.156	0.139	0.137	0.107	0.102	0.106	0.099	0.115	0.127	0.167
15	0.156	0.139	0.137	0.164	0.121	0.109	0.097	0.108	0.126	0.115
16	0.156	0.139	0.156	0.105	0.114	0.098	0.095	0.114	0.156	0.103
17	0.164	0.137	0.136	0.138	0.114	0.101	0.097	0.139	0.092	0.103
18	0.183	0.137	0.136	0.136	0.115	0.101	0.097	0.114	0.092	0.106

 

Выводы

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

Сформулирована постановка задачи восстановления сигналов, временных рядов и изображений в случае одномерной функции рассеяния точки. Представлен алгоритм, позволяющий снизить объем вычислений при нахождении значений весовой функции обратного фильтра.

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

Приведенные результаты апробации (см. табл. 2, 3, 4, 5) показывают, что использование данного подхода позволяет уменьшить погрешность восстановления и аппроксимации. Наибольшее преимущество при решении задач аппроксимации и восстановления данных представленный подход демонстрирует в случае весовой функции прямого фильтра типа (18). В настоящее время ведутся исследования, направленные на снижение погрешности восстановления сигналов, имеющих точки разрыва.

About the authors

Vitaly I. Batishchev

Samara State Technical University

Email: zolin.a.g@gmail.com

Dr. Sci. (Techn.), Professor

Russian Federation, 244, Molodogvardeyskaya st., Samara, 443100

Aleksey G. Zolin

Samara State Technical University

Author for correspondence.
Email: zolin.a.g@gmail.com

Ph.D. (Techn.), Associate Professor

Russian Federation, 244, Molodogvardeyskaya st., Samara, 443100

References

Supplementary files

There are no supplementary files to display.