CC BY
50
УДК 681.787:519.245
РЕКУРРЕНТНЫЙ АЛГОРИТМ ОБРАБОТКИ ИНТЕРФЕРОМЕТРИЧЕСКИХ СИГНАЛОВ НА ОСНОВЕ МУЛЬТИОБЛАЧНОЙ МОДЕЛИ
ПРЕДСКАЗАНИЯ М.А. Волынский3, И.П. Гурова, П.С. Скакова
а Университет ИТМО, Санкт-Петербург, Россия, [email protected]
Аннотация. Рассматривается модификация рекуррентного алгоритма обработки дискретной последовательности отсчетов интерферометрического сигнала, который основан на предсказании последующего отсчета сигнала при задании набора («облака») значений вектора параметров сигнала методом Монте-Карло, сравнении с измеренным значением сигнала и использованием полученной невязки для уточнения значений параметров сигнала на каждом шаге дискретизации. В предлагаемом модифицированном алгоритме использована концепция мультиоблачной модели предсказания, когда формируется набор нормально распределенных облаков с математическими ожиданиями, выбранными на основе критерия минимума невязки предсказания и наблюдения. Проведена экспериментальная апробация предлагаемого метода на примере оценивания начальной фазы полос в интерферометрии фазового сдвига, при этом дисперсия оценки сигнала, реконструированного по оцененной начальной фазе, относительно исходного сигнала не превышает 2% от максимального значения сигнала. Показано, что использование предложенного алгоритма позволяет избежать 2п-неоднозначности и обеспечить устойчивое восстановление фазы интерференционной картины сложного вида без использования априорной информации о распределении фазы интерференционных полос. Использование предлагаемого алгоритма применительно к оцениванию параметров интерферометрических сигналов позволяет повысить устойчивость фильтра к влиянию случайных помех и понизить требования к точности априорного задания параметров фильтрации по сравнению с обычной (однооблачной) реализацией последовательного метода Монте-Карло.
Ключевые слова: анализ интерферометрических сигналов, байесовский фильтр, последовательный метод Монте-Карло, мультиоблачное предсказание.
Благодарности. Работа выполнена при государственной финансовой поддержке ведущих университетов Российской Федерации (субсидия 074-U01) и в рамках гранта Президента Российской Федерации для государственной поддержки молодых российских ученых - кандидатов наук № МК-1455.2014.8.
THE RECURRENT ALGORITHM FOR INTERFEROMETRIC SIGNALS PROCESSING BASED ON MULTI-CLOUD PREDICTION MODEL M.A. Volynsky% I.P. Gurov% P.S. Skakov"
а1ТМО University, Saint Petersburg, Russia, [email protected]
Abstract. The paper deals with modification of the recurrent processing algorithm for discrete sequence of interferometric signal samples. The algorithm is based on subsequent reference signal prediction at specifying a set ("cloud") of values for signal parameters vector by Monte Carlo method, comparison with the measured signal value and usage of the residual for enhancing the values of signal parameters at each discretization step. The concept of multi-cloud prediction model is used in the proposed modified algorithm. A set of normally distributed clouds is created with expectation values selected on the base of criterion of minimum residual between prediction and observation values. Experimental testing of the proposed method applied to estimation of fringe initial phase in the phase shifting interferometry has been conducted. The estimate variance of the signal reconstructed according to estimated initial phase from initial signal does not exceed 2% of the maximum signal value. It has been shown that the proposed algorithm application makes it possible to avoid the 2n-ambiguity and ensure sustainable recovery of interference fringes phase of a complicated type without involving a priori information about interference fringe phase distribution. The usage of the proposed algorithm applied to estimation of interferometric signals parameters gives the possibility for improving the filter stability with respect to influence of random noise and decreasing requirements for accuracy of a priori filtration parameters setting as compared with conventional (single-cloud) implementation of the sequential Monte Carlo method.
Keywords: interferometric signals analysis, Bayesian filter, sequential Monte Carlo method, multi-cloud prediction. Acknowledgements. The work was done under government financial support for the leading universities of the Russian Federation (grant 074-U01) and within the Russian Federation President grant for the government support of young Russian PhD scientists № МК-1455.2014.8.
Введение
Интерферометрические методы измерений благодаря их высокой точности востребованы во многих областях науки и техники (см., например, [1, 2]). В реальных интерферометрических системах регистрируемые сигналы искажены случайными помехами различной природы. Повышение вычислительной мощности при использовании современных компьютерных технологий позволило реализовать ряд быстродействующих алгоритмов обработки интерферометрических сигналов и картин полос, устойчивых к случайным внешним воздействиям (см., например, [3, 4]), однако дальнейшее повышение устойчивости алгоритмов и минимизация погрешностей обработки данных в интерферометрических системах является актуальной задачей.
Алгоритмы, основанные на преобразовании Фурье, требуют наличия полной реализации сигнала, что ограничивает их быстродействие. Альтернативный подход состоит в использовании рекуррентных алгоритмов, что позволяет обрабатывать данные по мере их поступления. С учетом возможных случай-
ных помех при формировании и регистрации сигналов целесообразно использовать формализм стохастических дифференциальных уравнений для описания моделей сигналов и эволюции их параметров.
В последние годы активно исследовались рекуррентные алгоритмы обработки интерферометриче-ских сигналов на основе фильтра Калмана [4-8]. Результаты этих исследований позволили выявить ряд недостатков фильтра Калмана, главным из которых является неустойчивость фильтра в случае использования нелинейных математических моделей сигнала и (или) эволюции его параметров. Так, например, расширенный фильтр Калмана применительно к анализу интерферометрических сигналов предполагает локальную линеаризацию модели квазигармонического сигнала, что делает фильтр неоптимальным по критерию минимума средней квадратической ошибки оценки параметров сигнала. Особый интерес представляют алгоритмы, предполагающие замену использования нелинейных функциональных преобразований статистической аппроксимацией оцениваемых параметров. Одним из таких алгоритмов является сигма-точечный фильтр Калмана (Unscented Kaiman filter) [9], примененный в [8] для задачи обработки данных в оптической когерентной томографии. Сигма-точечный фильтр Калмана предполагает замену оцениваемых параметров на так называемые сигма-точки [9], вычисляемые на основе априорной информации о статистических моментах отклонений параметров, что позволяет отказаться от использования нелинейных функциональных преобразований с сохранением устойчивости алгоритма. Последовательный метод Монте-Карло [10, 11] основывается на более общем подходе по сравнению с сигма-точечным фильтром Калмана и состоит в статистической аппроксимации оцениваемых параметров с использованием множества («облака») значений параметров и соответствующего набора значений сигнала, что при обработке повышает устойчивость метода к случайным помехам благодаря выбору наиболее близких к истинным элементов облака на основе байесовского подхода.
При использовании общепринятой [10, 11] реализации последовательного метода Монте-Карло с однооблачным предсказанием в ряде случаев истинные значения параметров лежат близко к границам облака или за его пределами при прогнозировании их локализации в центре облака. Использование муль-тиоблачного предсказания позволяет более корректно оценивать параметры, истинные значения которых не попадают в центр конкретного облака. Таким образом, в настоящей работе предлагается модификация последовательного метода Монте-Карло путем использования мультиоблачной модели предсказания, что позволяет дополнительно повысить устойчивость фильтра к неточно заданным при инициализации алгоритма дисперсиям параметров.
Последовательный метод Монте-Карло с мультиоблачной моделью предсказания
Наблюдаемый сигнал может быть описан в дискретном времени t(k) = Mi, где к = 0, 1, ..., K-1 -номер дискретного отсчета, ht - шаг дискретизации, с помощью уравнения наблюдения
s (к ) = h(0(k )) + n(k ), (1)
где в - вектор параметров; h(-) - известная векторная функция; n - белый гауссовский шум наблюдения с нулевым средним. Задача оценивания элементов вектора в по набору наблюдений сигнала s сводится к оцениванию условного математического ожидания £[в | s]. Если функция h в уравнении (1) линейна, то задача решается с помощью линейного фильтра Калмана [10, 12]. В общем случае необходимо определить распределение вероятностей р(к) = р[в(к) | s(к)], после чего выбрать наиболее вероятное значение вектора параметров в .
Функцияр(к) может быть определена рекурсивно с помощью байесовской формулы [13]
р[в | s(к +1)] = рШ +1)|в]Р[в|*(к)] , (2)
J p[s( к+1)|в] р[в | s ( к )]^в
где р[в | s(к)] - предполагаемое распределение, спрогнозированное на шаге к,
p[s (к ) | в ] = [(2тс )к diag(Rn )Г0,5 exp{-0,5[s (к ) - h (в (к ))]Rn '[s (к ) - h (в (к ))]} , (3)
Rn - ковариационная матрица шума наблюдения; diag(Rn ) - ее диагональные элементы.
Фильтр, построенный на основе формул (1)-(3), может быть реализован с помощью последовательного метода Монте-Карло [10, 11]. Алгоритм включает следующие шаги: генерация набора из N векторов параметров с известным (прогнозируемым) законом распределения р[в | s(к)] и заданными статистическими моментами (математическим ожиданием и дисперсией); вычисление для каждого вектора весового коэффициента как величины, обратной невязке предсказания, с помощью функции h и наблюдения; отбор M векторов (M < N) с максимальными весами; их сдвиг в соответствии с известной моделью эволюции параметров; формирование нового набора из N векторов с математическим ожиданием, равным математическому ожиданию отобранных M векторов. Следует отметить, что рассматриваемый алгоритм аналогичен генерации случайных векторов с распределением (2) методом Монте-Карло. Детально реализация последовательного метода Монте-Карло применительно к обработке интерферометрических сигналов рассмотрена, например, в работе [14].
Использование одного облака векторов параметров для предсказания параметров как математического ожидания облака и генерация нового облака с вычисленным математическим ожиданием на следующем шаге фильтрации иногда не позволяют в достаточной степени минимизировать невязку предсказания и наблюдения на следующем шаге фильтрации, что особенно заметно в случае, если истинные значения параметров лежат близко к границам облака или за его пределами. Указанный недостаток однооб-лачного предсказания при обработке параметров интерферометрических сигналов проявляется, например, при оценивании фазы в точках, где значения фазы близки к 2п рад, и малые ошибки могут вызвать 2п-неустойчивость.
При использовании мультиоблачной модели предсказания предлагается заменить генерацию одного облака для каждого из параметров на генерацию М облаков с математическими ожиданиями, равными значениям параметров соответствующих М векторов. При этом предсказание значений параметров как математического ожидания полученного мультиоблака оказывается ближе к истинным значениям в случаях, когда однооблачный подход не позволяет корректно оценить значения параметров сигналов или вносит 2п-неустойчивость.
На рис. 1 иллюстрируется разница между однооблачной и мультиоблачной моделями предсказания. Видно, что математическое ожидание мультиоблака (серый круг на рис. 1, б) ближе к истинному значению параметров (черный круг на рис. 1), чем математическое ожидание, предсказанное в соответствие с традиционной интерпретацией алгоритма (прозрачный круг на рис. 1).
х сх х ^х)
- • По- Г ^ ■
а б
Рис. 1. Однооблачное (а) и мультиоблачное (б) предсказание. * - точки, выбранные по критерию наибольшей вероятности; о - однооблачное предсказание; • - мультиоблачное предсказание; • - наблюдение (истинное значение); эллипсами обозначены границы облаков (размеры полуосей
равны дисперсиям параметров)
Обработка интерферометрических сигналов на примере оценивания фазы
В простейшем случае скалярного наблюдения модель интерферометрического сигнала в уравнении (1) может быть задана как
И(в(к)) = Л (к) соб(Ф(к)), (4)
где Л(к) - амплитуда сигнала, Ф(к) - полная фаза, 0 = (Л,ф). В случае существенной нелинейности функции Ф(к) в выражении (4) интерес представляет способность фильтра правильно оценивать ее значение без дополнительной априорной информации о характере нелинейности. Для примера сгенерируем сигнал с квадратично возрастающей фазой Ф(к) (рис. 2) и зададим линейную модель эволюции параметров
0(к +1) = 0(к) + (0,Дф), (5)
где Дф - фиксированное приращение фазы. Для сигнала на рис. 1 закон изменения фазы Ф(к) = 0,62 к2, отношение сигнал/шум составляет 20 дБ.
1
ев
I 0,5
» Я
о о 0
Ш К
к е
и о
£ -0,5
га
И
м -1
0
50 100 150 200 Номер дискретного отсчета к
250
Рис. 2. Пример интерферометрического сигнала с квадратично возрастающей фазой и линейно
убывающей амплитудой
При использовании однооблачного подхода выбор параметра Дф в уравнении (5) существенно влияет на корректность оценки фазы: при выборе слишком малого значения Дф оценка фазы корректна лишь на первых отсчетах сигнала, где приращение фазы действительно невелико, при увеличении значения Дф оценка на последних отсчетах сигнала улучшается, однако появляются ошибки на 2п рад на пер-
вых отсчетах. На рис. 3 показаны результаты оценивания фазы сигнала, изображенного на рис. 2, при использовании однооблачной (рис. 3, а, б) и мультиоблачной (рис. 3, в) модели предсказания. Входные параметры алгоритма: генерируемое количество точек - 200, отбираемое количество точек - 10, среднее квадратичное отклонение амплитуды и фазы - 0,6 усл. ед. и 0,03 рад соответственно, приращение фазы -1 рад для рис. 3, а, в, и 1,2 рад для рис. 3, б. На рис. 3, а, б, пунктирными линиями иллюстрируются истинные значения фазы, сплошными линиями - оценки фазы с описанными ошибками. Использование мультиоблачного предсказания делает алгоритм робастным по отношению к выбору значения параметра Дф (рис. 3, в).
250 200
150
100
§ я
я О
250 200
150
й 100 я
0 50 100 150 200 250 Номер дискретного отсчета к
50
0 50 100 150 200 250 Номер дискретного отсчета к б
я О
50
0 50 100 150 200 250 Номер дискретного отсчета к
а О в
Рис. 3. Оценка фазы сигнала (рис. 2) с использованием однооблачного предсказания с малым (а) и большим (б) приращением фазы Дф (пунктирные линии - истинные значения фазы, сплошные линии -оценки фазы с описанными ошибками) и с использованием мультиоблачной модели предсказания (в)
Оценивание фазы сигнала с помощью последовательного метода Монте-Карло с мультиоблачной моделью предсказания может быть использовано при анализе экспериментальных картин интерференционных полос в интерферометрии фазового сдвига [4, 15, 16], поскольку фазовый сдвиг, вносимый, например, с помощью пьезоэлемента, в условиях реального эксперимента оказывается не всегда одинаковым и не идеально линейным.
а б
Рис. 4. Интерференционная картина с нормированной амплитудой (а) и оценка начальной фазы с использованием предложенного алгоритма (б). Размер исследуемой поверхности 3*3 мм2
На рис. 4 иллюстрируется результат работы алгоритма при восстановлении начальной фазы полос для серии из 50 интерференционных картин, зарегистрированных при исследовании гладкой металлической поверхности в интерферометре Майкельсона. На рис. 4, а, представлена интерференционная картина, завершающая серию. Дисперсия оценки сигнала, реконструированного по оцененной начальной фазе исходного сигнала, не превышает 2% от максимального значения сигнала. Следует отметить, что при использовании алгоритма с однооблачной моделью предсказания и теми же выходными параметрами дисперсия составляет 34%, что свидетельствует о расходимости фильтра для примерно трети точек изображения. Из рис. 4 видно, что использование предложенного алгоритма позволяет избежать 2п-неоднозначности и обеспечить устойчивое восстановление фазы интерференционной картины сложного вида без использования априорной информации о распределении фазы интерференционных полос.
Заключение
В работе предложен модифицированный последовательный метод Монте-Карло на основе мульти-облачной модели предсказания. Использование алгоритма применительно к оцениванию параметров ин-терферометрических сигналов позволяет повысить устойчивость фильтра к влиянию случайных помех и понизить требования к точности априорного задания параметров фильтрации по сравнению с обычной (однооблачной) реализацией последовательного метода Монте-Карло.
Экспериментальная апробация предложенной модификации алгоритма применительно к обработке интерферограмм фазового сдвига показала возможность устойчивого восстановления фазы в условиях искажений распределения фазы.
Литература
4.
5.
6.
7.
9.
Malacara D. Optical Shop Testing. 2nd ed. NY: Wiley, 1992. 792 p.
Hariharan P. Optical interferometry // Reports on Progress in Physics. 1991. V. 54. N 3. P. 339-390. Huntley J.M. Automated fringe pattern analysis in experimental mechanics: a review // Journal of Strain Analysis for Engineering Design. 1998. V. 33. N 2. P. 105-125.
Васильев В.Н., Гуров И.П. Компьютерная обработка сигналов в приложении к интерферометриче-ским системам. СПб: БХВ-Санкт-Петербург, 1998. 240 с.
Gurov I., Ermolaeva E., Zakharov A. Analysis of low-coherence interference fringes by the Kalman filtering method // Journal of the Optical Society of America A. 2004. V. 21. N 2. P. 242-251. Gurov I., Volynsky M. Interference fringe analysis based on recurrence computational algorithms // Optics and Lasers in Engineering. 2012. V. 50. N 4. P. 514-521.
Захаров А.С. Нелинейный анализ стохастических параметров интерференционных систем: дис. ... канд. техн. наук ... по спец. 05.13.01. защищена 20.12.05. СПб, 2005. 157 с.
Волынский М.А. Рекуррентные алгоритмы обработки данных в оптической когерентной томографии: дис. ... канд. техн. наук ... по спец. 05.13.01. защищена 20.12.11. СПб, 2011. 112 с. Wan E.A., van der Merwe R. The unscented Kalman filter / In: Kalman Filtering and Neural Networks. NY: John Wiley & Sons, 2001. P. 221-280.
10. Simon D. Optimal state estimation: Kalman, H®, and Nonlinear Approaches. NY: John Wiley & Sons, Inc., 2006. 526 p.
11. Doucet A., de Freitas N., Gordon N. Sequential Monte Carlo methods in practice. NY: Springer-Verlag, 2001. 583 p.
12. Перов А.И. Статистическая теория радиотехнических систем. М.: Радиотехника, 2003. 400 с.
13. 0ksendal B.K. Stochastic Differential Equations: An Introduction with Applications. 6th ed. Berlin: Springer, 2003. 379 p.
14. Волынский М.А., Гуров И.П., Ермолаев П.А., Скаков П.С. Динамическое оценивание параметров ин-терферометрических сигналов на основе последовательного метода Монте-Карло // Научно-технический вестник информационных технологий, механики и оптики. 2014. № 3 (91). С. 18-23.
15. Zhao W, Cao G. A real-time adaptive phase-shifting interferometry // Proc. of SPIE - The International Society for Optical Engineering. 2012. V. 8493. Art. N 849313.
16. Chen L.-C., Yeh S.-L., Tapilouw A.M., Chang J.-C. 3-D surface profilometry using simultaneous phase-shifting interferometry // Optics Communications. 2010. V. 283. N 18. P. 3376-3382.
Волынский Максим Александрович Гуров Игорь Петрович Скаков Павел Сергеевич
Maxim A. Volynsky Igor P. Gurov Pavel S. Skakov
кандидат технических наук, доцент, Университет ИТМО, Санкт-
Петербург, Россия, [email protected]
доктор технических наук, профессор, зав. кафедрой, Университет
ИТМО, Санкт-Петербург, Россия, [email protected]
ассистент, Университет ИТМО, Санкт-Петербург, Россия,
PhD, Associate professor, ITMO University, Saint Petersburg, Russia, [email protected]
D.Sc., Professor, Department head, ITMO University, Saint Petersburg, Russia, [email protected]
assistant, ITMO University, Saint Petersburg, Russia, [email protected]
Принято к печати 02.04.14 Accepted 02.04.14
