Динамическое Оценивание параметров интерферометрических сигналов методом расширенной фильтрации Калмана второго порядка Текст научной статьи по специальности «Электротехника, электронная техника, информационные технологии»

ФОТОНИКА И ОПТОИНФОРМАТИКА PHOTONICS AND OPTOINFORMATICS

УДК 681.787:[519.24+519.6]

ДИНАМИЧЕСКОЕ ОЦЕНИВАНИЕ ПАРАМЕТРОВ ИНТЕРФЕРОМЕТРИЧЕСКИХ СИГНАЛОВ МЕТОДОМ РАСШИРЕННОЙ ФИЛЬТРАЦИИ КАЛМАНА ВТО РОГО ПОРЯДКА1 П.А. Ермолаев"

a Санкт-Петербургский национальный исследовательский университет информационных технологий, механики и оптики (Университет ИТМО), Санкт-Петербург, Россия, [email protected]

К алгоритмам обработки данных в интерферометрических системах предъявляются требования высокого разрешения и быстродействия. В рекуррентных алгоритмах, основанных на параметрическом представлении сигналов, осуществляется последовательная обработка отсчетов сигнала. В ряде случаев рекуррентные алгоритмы позволяют повысить быстродействие и увеличить разрешение по сравнению с классическими методами обработки. Зависимость наблюдаемого интерферометрического сигнала от параметров его модели и стохастической природы формирования шумов в системе является в общем случае нелинейной. Для обработки таких сигналов целесообразно применение алгоритмов нелинейной стохастической фильтрации. Примером такого алгоритма может служить расширенный фильтр Кал-мана, основанный на линеаризации уравнений системы и наблюдения с использованием первых производных по вектору параметров. Для уменьшения ошибки аппроксимации нелинейных уравнений модели предложено использовать расширенный фильтр Калмана второго порядка, в котором при оценивании параметров дополнительно учитываются вторые производные уравнений модели по вектору параметров. Приведены примеры реализации алгоритма для различных наборов оцениваемых параметров. Предложенный алгоритм позволяет повысить качество обработки данных в интерферометрических системах, сигналы в которых формируются в соответствии с рассмотренными в работе моделями. Полученное среднеквадратическое отклонение оценки амплитуды от ее максимального значения не превысило 4%. Показано, что отношение сигнал/шум восстановленного сигнала увеличено по сравнению с исходным более чем на 60%.

Ключевые слова: расширенный фильтр Калмана второго порядка, анализ интерферометрических сигналов.

DYNAMIC ESTIMATION FOR PARAMETERS OF INTERFERENCE SIGNALS BY THE SECOND ORDER EXTENDED KALMAN FILTERING1

P.A. Ermolaev"

a Saint Petersburg National Research University of Information Technologies, Mechanics and Optics (ITMO University), Saint Petersburg, Russia, [email protected]

Data processing in the interferometer systems requires high-resolution and high-speed algorithms. Recurrence algorithms based on parametric representation of signals execute consequent processing of signal samples. In some cases recurrence algorithms make it possible to increase speed and quality of data processing as compared with classic processing methods. Dependence of the measured interferometer signal on parameters of its model and stochastic nature of noise formation in the system is, in general, nonlinear. The usage of nonlinear stochastic filtering algorithms is expedient for such signals processing. Extended Kalman filter with linearization of state and output equations by the first vector parameters derivatives is an example of these algorithms. To decrease approximation error of this method the second order extended Kalman filtering is suggested with additionally usage of the second vector parameters derivatives of model equations. Examples of algorithm implementation with the different sets of estimated parameters are described. The proposed algorithm gives the possibility to increase the quality of data processing in interferometer systems in which signals are forming according to considered models. Obtained standard deviation of estimated amplitude envelope does not exceed 4% of the maximum. It is shown that signal-to-noise ratio of reconstructed signal is increased by 60%. Keywords: the second order extended Kalman filter, interferometer data analysis.

Введение

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

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

1 Работа выполнена при государственной финансовой поддержке ведущих университетов Российской Федерации (субсидия 074-U01) и в рамках гранта Президента Российской Федерации для государственной поддержки молодых российских ученых - кандидатов наук № МК-1455.2014.8.

The work is done with the state financial support of the leading universities of the Russian Federation (subsidy 074-U01) and under a grant of the President of the Russian Federation for state support of young Russian PhD scientists № MK-1455.2014.8.

1

информации о свойствах объекта в ряде случаев требуются быстродействующие алгоритмы и методы обработки сигналов [2].

Алгоритмы обработки, основанные на преобразовании Фурье, не всегда удовлетворяют критериям быстродействия и требуют в качестве входных данных полную реализацию сигнала, что затрудняет возможность их использования в реальном времени. Существует класс рекуррентных алгоритмов, основанных на представлении состояния динамической системы с помощью стохастических дифференциальных уравнений [2, 3]. В этих алгоритмах используется априорная информация о модели формирования сигналов и статистических характеристиках шумов. Зависимость регистрируемого сигнала от параметров модели в общем случае является нелинейной, поэтому большой интерес для обработки данных в интерфе-рометрических системах представляют алгоритмы нелинейной фильтрации Калмана (см., например, [4, 5]). Наиболее распространенным среди таких алгоритмов является расширенный фильтр Калмана (РФК) [4, 5], в котором осуществляется линеаризация математической модели при помощи разложения в ряд Тейлора и отбрасывания составляющих, содержащих производные выше первого порядка.

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

Расширенный фильтр Калмана второго порядка

Как известно, линейный фильтр Калмана оптимален по критерию минимума среднеквадратической ошибки [7], однако неприменим для нелинейных систем. В этом случае применяются алгоритмы нелинейной фильтрации Калмана, являющиеся лишь квазиоптимальными [2, 8]. Одним из таких алгоритмов является РФК [5]. В РФК второго порядка [6] при аппроксимации нелинейных уравнений модели учитываются дополнительные составляющие ряда Тейлора, содержащие вторые производные уравнений системы и наблюдения по вектору параметров, что позволяет повысить точность аппроксимации этих уравнений.

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

0(к) = f (0(к -1)) + «т(к),

»(к) = Ь(0(к)) + п(к),

где 0(к) - вектор параметров системы; »(к) - результаты наблюдений; к = 0..К-1 - номер дискретного отсчета; f(•) и Ь() - известные нелинейные дифференцируемые векторные функции, описывающие соответственно эволюцию системы во времени и полезную составляющую наблюдаемого сигнала; «'(к) и п(к) - шум системы и шум наблюдения соответственно.

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

0(k ) = f (0(k-1)) + 2 £ ф Tr

2 1=1

R(k -1)

в( k-1)

30 2

где f. имеет смысл /-го уравнения в векторной функции системы; N - количество оцениваемых параметров; Тг[-] - оператор вычисления следа матрицы. Вектор ф,. представляется как вектор-столбец размера N значения всех элементов которого, кроме /-го, равны нулю, а /-й элемент равен единице; Я (к -1) -ковариационная матрица ошибок на предыдущем шаге.

Предсказание ковариационной матрицы ошибок имеет вид

Я(к) = Е(к)Я(к - 1)Е(к)т + Я„ , (2)

где Яп - ковариационная матрица шума системы, а Е(к) является первой производной уравнения системы по вектору параметров:

^ (0)

F(k)=-

50

(3)

0( к-1)

Коррекция вектора параметров с учетом наблюдаемого значения сигнала представлена как

0(к) = 0(к) + Р(к)[8(к) - И(0(к))] - п(к), (4)

где Р(к) является коэффициентом усиления фильтра:

Р(к) = Я(к)Н(к)т (и(к)Я(к)Н(к)т + Яя )-1. (5)

Коэффициент, учитывающий вклад второй производной функции наблюдения по вектору параметров, п(к), представляется в форме

1 M

n(k) = -P(k)£ ф,Tr [D¿ (k)Rpr (k)] .

(б)

В уравнениях (5)-(6) Яп - ковариационная матрица шума наблюдения; Н(к) - первая производная функции наблюдения по вектору параметров; Б,, (к) - матрицы, учитывающие вторые производные /-го компонента векторного уравнения наблюдения: 5Ь(0)

H(k ) =-

D, (k ) =

se

е( k )

s 2h,. (e)

se2

(7)

(8)

е( k )

Коррекция ковариационной матрицы ошибок осуществляется с учетом предварительно рассчитанного коэффициента усиления

Я (к) = (I - Р (к )Н (к ))К (к), (9)

где I - единичная матрица. Подробнее РФК второго порядка рассмотрен в [6].

Оценивание параметров интерферометрических сигналов

Для реализации алгоритма необходимо знать модель формирования наблюдений в интерферомет-рической системе. Интерферометрический сигнал можно представить как дискретную последовательность отсчетов в виде

ь(к) = В(к) + А(к)соБ(Ф(к) + 5ф(к)) + п(к), (10)

где В(к) - фоновая составляющая; А(к) - амплитуда; п(к) - белый гауссовский шум с нулевым средним, не коррелированный с сигналом; Ф (к) - полная фаза сигнала; 5ф(к) - случайные флуктуации фазы. Полная фаза сигнала имеет вид

Ф(к) = Х 2п/(к')Аг ,

к'=0

где /(к') - частота; А2 - шаг дискретизации, имеющий смысл изменения разности хода опорной и измерительной волны в исследуемой точке объекта при перемещении опорного отражателя в интерферометре. С учетом модели интерферометрического сигнала (10) вектор параметров и уравнения системы и наблюдения можно представить в виде

0 = (В, А, /, Ф)г,

И(0) = В + Асоб Ф,

f (0) = 0 + (0 0 0 2 п/А2 ) .

С учетом соотношений (3), (7) и (11)—(13) могут быть рассчитаны матрицы Е(к) и Н(к): (1 0 0 0^

(ll) (l2)

(13)

F(k ) =

sf (e)

se

e(k-l)

о 1 о о о о 1 о о о 2nAz 1

H(k ) =

sh(e)

se

= (l cos Ф о - A sin Ф).

е( k )

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

Векторная функции наблюдения (12) представляется единственным уравнением, следовательно, существует лишь одна матрица D(k):

0 ^ - sin Ф

D(k ) =

s2h.. (e)

se2

е( k )

Г о о о

о

о - sin Ф о - A cos Ф

i=l

Оценку вектора параметров (11) получают путем последовательного применения соотношений (1), (2), (4) и (9) к поочередно поступающим на вход алгоритма отсчетам исходного сигнала.

На рис. 1 представлены результаты оценивания огибающей и фазы модельного интерферометриче-ского сигнала. Отношения сигнал/шум исходного и восстановленного сигналов составляют соответственно 10,7 дБ и 19,1 дБ. Среднее квадратичное отклонение ошибки оценки амплитуды составляет 3,7% от ее максимального значения и имеет нормальное распределение. Максимальное отклонение оценки фазы не превышает п/10 рад. С уменьшением числа отсчетов, на которые приходится полезная составляющая сигнала, в 10 раз среднее квадратичное отклонение амплитуды увеличивается до 7%, максимальное отклонение фазы - до п/4, а отношение сигнал/шум восстановленного сигнала снижается до 15,2 дБ в зависимости от характеристик ковариационной матрицы ошибок. На рис. 1, в, видно, что наибольшие отклонения значения фазы имеют до появления в сигнале полезной составляющей.

200

S ^

и

и

S о

(D

S

и

(D

£ и

го

1

0,5

ч

" 0 «

о -0,5

-1

j г \

- Ji

I'll IP

J 1 ¡I

250 500 750 1000 Номер дискретного отсчета а

Я 0 5 3 0,5 0,25

д 0

Н-0,25

-0,5

«

о

н

ч о да

a a-

е

и

н е ч а -e

а

н

CO

150 100

250 500 750 1000 Номер дискретного отсчета б

(D

S

и

(D

и

8 W н О

0 250 500 750 1000 Номер дискретного отсчета в

Рис. 1. Результаты оценивания огибающей интерферометрического сигнала (а), изменения фазы (б) и график отклонения оценки фазы от линейного закона (в)

Оценивание положения максимумов огибающей сигналов

Для систем корреляционной оптической когерентной томографии, основанной на принципах интерферометрии малой когерентности и широко применяющейся в исследованиях биологических тканей и структуры многослойных покрытий [9-16], характерны сигналы, огибающая которых представляется суммой кривых гауссовой формы [10]. Положения максимумов огибающей такого сигнала соответствуют положениям границ слоев исследуемой среды [9, 10]. Амплитуду таких сигналов Л (к) целесообразно представить в виде

'-(к& - 2))2"

A(k ) = Х A (k)exp

(14)

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

Учитывая форму представления амплитуды (14), вектор параметров (11) можно дополнить параметром z , обозначающим положение текущего максимума огибающей:

0 = (B, A, f, Ф, z)T .

При этом параметром A является максимальное значение амплитуды сигнала (14). Тогда уравнения системы и наблюдения и матрица F(k) выражаются в следующей форме:

-(Mz - z)2\ _ 'cos Ф ,

h(ö) = B+A exp ^ f (0 ) = о + (0 0 0 2 nf Az 0 y

0

0

2

a

F(0(k )) =

5f (0)

50

0( k-1)

( 1 0 0 0 01

0 1 0 0 0

= 0 0 1 0 0

1) 0 0 2nAz 1 0

,0 0 0 0 1,

Значения вторых производных всех компонентов векторного уравнения системы по вектору параметров равны нулю, аналогично предыдущей рассмотренной модели. Матрицы Б(к) и Н(к) могут быть вычислены при помощи уравнений (7) и (8).

На рис. 2 представлены результаты оценивания положения максимума огибающей интерферомет-рического сигнала. Реальные максимумы огибающей находятся в дискретных отсчетах под номерами 400 и 600. Видно, что предсказание положения первого максимума осуществляется до его фактического появления в сигнале (в точке, в которой амплитуда составляет около 50% от максимального значения).

250 500 750 1000 Номер дискретного отсчета

0 250 500 750 1000 Номер дискретного отсчета б

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

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

При оценивании положения максимума огибающей не требуется наличие априорной информации о количестве и взаимном расположении границ слоев исследуемой среды, однако качество работы фильтра в значительной степени зависит от начальных условий: априорное значение параметра 2 должно быть меньше координаты расположения реальных максимумов, иначе они не будут учтены в процессе оценивания. Среднее квадратичное отклонение нормально распределенной ошибки восстановленной амплитуды составляет 3,4% от ее максимального значения. Соотношение сигнал/шум, восстановленное по результатам оценивания сигнала, было увеличено с 14,6 до 24,1 дБ по сравнению с соотношением сигнал/шум исходного сигнала.

Заключение

В работе рассмотрено решение задачи динамического оценивания параметров интерферометриче-ских сигналов в условиях воздействия шума при помощи расширенного фильтра Калмана второго порядка. Указаны особенности реализаций алгоритма для различных моделей сигнала, в том числе для модели с несколькими максимумами огибающей, характерной для сигналов в оптической когерентной томографии. Максимальное отклонение результата оценивания фазы от линейного закона не превышает п/10 рад при начальном отношении сигнал/шум порядка 10 дБ и нескольких десятках отсчетов на одном периоде сигнала. Наибольшие отклонения фаза имеет до появления полезной составляющей сигнала. Ошибки оценивания амплитуды с использованием обеих моделей не превышают 4% от ее максимального значения. Соотношение сигнал/шум восстановленного сигнала увеличено более чем на 60% по сравнению с исходным сигналом. Уменьшение числа отсчетов, на которые приходится полезная составляющая сигнала, ведет к увеличению ошибки оценивания амплитуды и фазы, понижению отношения сигнал/шум восстановленного сигнала и к отсутствию устойчивых значений оценки положений максимумов огибающих. Предложенный алгоритм удовлетворяет требованиям к быстродействию и качеству обработки данных в интерферометри-ческих системах, сигналы в которых формируются в соответствии с рассмотренными моделями.

Расширенный фильтр Калмана второго порядка является квазиоптимальным алгоритмом. Результаты его работы и сходимость в значительной степени зависят от априорной информации о модели фор-

0

а

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

Литература

1. Malacara D. Optical Shop Testing. NY: Wiley, 1978. 862 p.

2. Gurov I., Volynsky M. Interference fringe analysis based on recurrence computational algorithms // Optics and Lasers in Engineering. 2012. V. 50. № 4. P. 514-521.

3. Van Kampen N. Stochastic Processes in Physics and Chemistry. North Holland, 1984. 464 p.

4. Simon D. Using Nonlinear Kalman Filtering to Estimate Signals // Embedded Systems Design. 2006. V. 19. N 7. P. 38-53.

5. Волынский М.А., Гуров И.П., Захаров А.С. Динамический анализ сигналов в оптической когерентной томографии методом нелинейной фильтрации Калмана // Оптический журнал. 2008. T. 75. № 10. C. 89-94.

6. Simon D. Optimal state estimation: Kalman, H®, and Nonlinear Approaches. NY: John Wiley & Sons, Inc., 2006. 526 p.

7. Kalman R.E. A new approach to linear filtering and prediction problems // Trans. ASME, J. Basic Eng. 1960. V. 82. P. 35-45.

8. Gurov I., Ermolaeva E., Zakharov A. Analysis of low-coherence interference fringes by the Kalman filtering method // JOSA A. 2004. V. 21. № 2. P. 242-251.

9. Гуров И.П. Оптическая когерентная томография: принципы, проблемы и перспективы // Проблемы когерентной и нелинейной оптики / Под ред. И.П. Гурова, С.А. Козлова. СПб: СПбГУ ИТМО, 2004. С. 6-30.

10. Gurov I., Volynsky M., Zakharov A. Evaluation of multilayer tissues in optical coherence tomography by the extended Kalman filtering method // Proc. SPIE. 2007. V. 6734. art. N 67341P.

11. Dresel T., Häusler G., Ventzke H. Three-dimensional sensing of rough surfaces by coherence radar // Appl. Opt. 1992. V. 31. P. 919-925.

12. Deck L. de Groot P. High-speed non-contact profiler based on scanning white light interferometry // Appl. Opt. 1994. V. 33. P. 7334-7338.

13. Fercher A. Optical coherence tomography // J. Biomed. Opt. 1996. V. 1. N 2. P. 157-173.

14. Гуров И.П., Жукова Е.В., Маргарянц Н.Б. Исследование внутренней микроструктуры материалов методом оптической когерентной микроскопии с перестраиваемой длиной // Научно-технический вестник информационных технологий, механики и оптики. 2012. № 3 (79). С. 40-45.

15. Гуров И.П., Жукова Е.В., Левшина А.В. Применение метода оптической когерентной томографии для изучения предметов искусства, выполненных в технике интарсии // Научно-технический вестник информационных технологий, механики и оптики. 2012. № 3 (79). С. 55-59.

16. Волынский М.А., Воробьева Е.А., Гуров И.П, Маргарянц Н.Б. Бесконтактный контроль микрообъектов методами интерферометрии малой когерентности и оптической когерентной томографии // Изв. вузов. Приборостроение. 2011. Т. 54. № 2. С. 75-82.

17. Ярлыков М.С. Статистическая теория радионавигации. М: Радио и связь, 1985. 344 с.

18. Gurov I., Sheynihovich D. Interferometric data analysis based on Markov nonlinear filtering methodology // JOSA A. 2000. V. 17. N 1. P. 21-27.

Ермолаев Петр Андреевич -

Petr A. Ermolaev

Принято к печати 10.01.14 Accepted 10.01.14

студент, Санкт-Петербургский национальный исследовательский университет информационных технологий, механики и оптики (Университет ИТМО), Санкт-Петербург, Россия, [email protected]

student, Saint Petersburg National Research University of Information Technologies, Mechanics and Optics (ITMO University), Saint Petersburg, Russia, [email protected]