CC BY
36
УДК 681.787:519.245
ДИНАМИЧЕСКОЕ ОЦЕНИВАНИЕ ПАРАМЕТРОВ ИНТЕРФЕРОМЕТРИЧЕСКИХ СИГНАЛОВ НА ОСНОВ Е ПОСЛЕДОВАТЕЛЬНОГО МЕТОДА МОНТЕ-КАРЛО1
М.А. Волынский3, И.П. Гурова, П.А. Ермолаев3, П.С. Скакова
а Университет ИТМО, Санкт-Петербург, Россия, [email protected]
Рассмотрен последовательный метод Монте-Карло применительно к задаче оценивания параметров интерферомет-рических сигналов, основанный на статистической аппроксимации апостериорной плотности вероятности распределения параметров. Приведено детальное описание алгоритма. Показана возможность использования минимума невязки предсказания и наблюдения в качестве критерия отбора элементов генерируемого на каждом шаге алгоритма множества векторов параметров. Проведен анализ влияния входных параметров на качество работы алгоритма. Получено, что среднее квадратичное отклонение ошибки оценки амплитуды типичных сигналов составило около 10% от максимального значения амплитуды. Показано, что ошибка оценки фазы имеет нормальное распределение. Проведен анализ характеристик работы алгоритма в зависимости от входных параметров, в частности, проведен анализ влияния количества отбираемых векторов параметров на результаты оценивания. На основании результатов моделирования для рассмотренного класса сигналов рекомендуется отбирать 30% от количества генерируемых векторов. Увеличение количества генерируемых векторов более 150 не влечет за собой значительного улучшения качества получаемых оценок. Последовательный метод Монте-Карло рекомендуется к использованию при динамической обработке интерферометрических сигналов для случаев, когда требуется повышенная устойчивость к нелинейностям изменения параметров сигнала и влиянию случайных помех.
Ключевые слова: интерферометрические сигналы, последовательный метод Монте-Карло.
DYNAMIC PARAMETERS ESTIMATION OF INTERFEROMETRIC SIGNALS BASED ON SEQUENTIAL MONTE CARLO METHOD1 M.A. Volynsky3, I.P. Gurov3, RA. Ermolaev3, P.S. Skakov3
а ITMO University, Saint Petersburg, Russia, [email protected]
The paper deals with sequential Monte Carlo method applied to problem of interferometric signals parameters estimation. The method is based on the statistical approximation of the posterior probability density distribution of parameters. Detailed description of the algorithm is given. The possibility of using the residual minimum between prediction and observation as a criterion for the selection of multitude elements generated at each algorithm step is shown. Analysis of input parameters influence on performance of the algorithm has been conducted. It was found that the standard deviation of the amplitude estimation error for typical signals is about 10% of the maximum amplitude value. The phase estimation error was shown to have a normal distribution. Analysis of the algorithm characteristics depending on input parameters is done. In particular, the influence analysis for a number of selected vectors of parameters on evaluation results is carried out. On the basis of simulation results for the considered class of signals, it is recommended to select 30% of the generated vectors number. The increase of the generated vectors number over 150 does not give significant improvement of the obtained estimates quality. The sequential Monte Carlo method is recommended for usage in dynamic processing of interferometric signals for the cases when high immunity is required to non-linear changes of signal parameters and influence of random noise. Keywords: interferometric signals, sequential Monte Carlo method.
Введение
Оптические бесконтактные высокоразрешающие методы исследования объектов широко востребованы в науке (медицина, биология, материаловедение) и научных приложениях (криминалистика, контроль качества производства) [1]. Благодаря высокой точности, наибольшее распространение получили интерферометрические методы. При этом получаемые данные представляют собой совокупность квазигармонических сигналов в смеси с помехами, содержащих информацию о свойствах объекта [2]. Для извлечения полезной информации применяются различные алгоритмы обработки, к которым предъявляются требования по быстродействию и устойчивости к различным типам сигналов и шумов [2-4].
Традиционно для анализа сигналов широко применяются алгоритмы, основанные на преобразовании Фурье, однако эти алгоритмы требуют наличия полной реализации сигнала, что ограничивает быстродействие. В рамках системного подхода к процессам формирования и анализа сигналов в интерферометре можно задать динамическую систему с помощью стохастических дифференциальных уравнений. Параметры такой системы могут быть оценены с помощью рекуррентных алгоритмов обработки данных [2-4]. Эти алгоритмы осуществляют последовательную обработку результатов измерений, используя априорную информацию о модели сигнала и статистических характеристиках шума.
1 Работа выполнена при финансовой поддержке Министерства образования и науки Российской Федерации и гранта Президента Российской Федерации для государственной поддержки молодых российских ученых - кандидатов наук № МК-1455.2014.8.
1 The work was done under financial support from the Russian Ministry of Education and Science and grant of the President
of the Russian Federation for state supporting of young Russian PhD scientists № МК-1455.2014.8.
Линейный фильтр Калмана [5] является оптимальным по критерию минимума среднего квадратичного отклонения (СКО) ошибки оценивания, однако этот алгоритм применим только для линейных систем с нормально распределенным аддитивным шумом. Формируемые в интерферометрических системах сигналы нелинейно зависят от параметров. Для обработки таких сигналов применяются алгоритмы нелинейной стохастической фильтрации. Примерами таких алгоритмов могут служить расширенный фильтр Калмана [6, 7] и нелинейный марковский фильтр [8, 9], основанные на аппроксимации нелинейных уравнений модели при помощи разложения в степенной ряд. Эти алгоритмы квазиоптимальны [2, 10] и неустойчивы к отклонениям начальных условий и шуму с распределением, отличным от нормального (см., например, [3]).
Для повышения устойчивости к различным типам нелинейности и распределениям шума в работе предлагается использовать последовательный метод Монте-Карло (Sequential Monte-Carlo method) [7, 11, 12], также упоминаемый в зарубежной литературе как фильтр частиц (Particle filter) [13], конденсационный алгоритм (The condensation algorithm) [14-15] или аппроксимация взаимодействующими частицами (Interacting particle approximation) [16]. Этот алгоритм базируется на статистической аппроксимации апостериорной плотности вероятности каждого параметра в векторе параметров системы на основании предыдущих наблюдений [7, 11-13].
Поскольку параметры интерферометрических сигналов нелинейно входят в модель наблюдения, а также сигналы подвержены влиянию случайных внешних воздействий, целесообразно исследовать возможность применения последовательного метода Монте-Карло в задаче динамического оценивания параметров интерферометрических сигналов, а также определить свойства алгоритма, включая его устойчивость к шумам и влияние начальных настроек на сходимость. Таким образом, в работе рассмотрено применение последовательного метода Монте-Карло для оценивания параметров интерферометрических сигналов с нелинейно изменяющейся амплитудой и фазой. Проведен анализ влияния входных параметров на качество работы алгоритма.
Последовательный метод Монте-Карло
Динамическая система, как известно, может быть задана в дискретной форме для каждого номера дискретного отчета к = 0..K-1 уравнениями системы и наблюдения:
в(к +1) = f (в(к), w(k )), (1)
s (к ) = h (в (к), n(k )), (2)
где в - вектор параметров системы с известной плотностью вероятности p(в(к)) на шаге к = 0, f() и h () - известные нелинейные дифференцируемые векторные функции, w и n - шумы системы и наблюдения, статистический закон распределения которых известен.
Работа алгоритма условно может быть разделена на четыре этапа: генерация случайного набора векторов параметров системы, предсказание возможных значений параметров на следующем шаге, отбор векторов, лучше всего удовлетворяющих поступившим наблюдениям, и коррекция плотности вероятности распределения параметров (рис. 1).
Рис. 1. Схема алгоритма динамического оценивания параметров последовательным методом Монте-Карло
Входными параметрами алгоритма является количество генерируемых случайных векторов М, пороговая вероятность отбора, статистические моменты априорной плотности вероятности распределения параметров. Рассмотрим каждый этап подробнее.
На первом этапе с учетом информации о распределении шумов и компонентов вектора параметров динамической системы генерируется множество в векторов, состоящее из N независимых векторов 0(£,/), где I = 0..М-1 - номер вектора 0 в множестве в. Количество генерируемых векторов выбирается с учетом
допустимых погрешностей оценки параметров и требованиям к быстродействию. На втором этапе в соответствии с уравнением (1) формируется множество предсказываемых значений вектора параметров:
0(k +1,i) = f (((k,i),w(k,i)).
На третьем этапе из элементов множества © выбираются векторы, лучше всего удовлетворяющие наблюдениям, зарегистрированным на текущем шаге. Этот выбор осуществляется на основе оценки вероятности совпадения каждого из векторов множества © с истинным вектором параметров на текущем шаге. Для этого с использованием уравнения (2) для каждого из векторов множества © вычисляется оценка наблюдения s(k). Описанную выше вероятность можно оценить как:
q(i) = p(k)=s(k)]|[0(k,i)=0(k)]) . (3)
Формула для вычисления q (i) зависит от характера функции распределения шума наблюдений. Например, если шум наблюдения n(k) аддитивен и распределен по нормальному закону с нулевым средним, уравнение (3) можно представить в виде (см., например, [7]):
q(i) = P(n(k)={s(k)-h(0(k,i))}) = ^ exp-2[s(k)-h(0(k,i))]T R-1[s(k)-h(0(k,i))] j, (4)
где m - количество элементов в векторе наблюдаемых значений; R - ковариационная матрица шума наблюдения; а - коэффициент нормировки. Выбор наиболее вероятных векторов осуществляется в соответствии с правилом
{V9(k,i)e©:q(i)<p], p e [0,1], (5)
где p - пороговое значение вероятности, определяющее минимальную условную вероятность совпадения вектора из множества © с истинным вектором параметров динамической системы.
На четвертом этапе проводится оценивание вектора параметров (например, он может быть вычислен как среднее значение выбранных векторов). На основании выбранных значений осуществляется коррекция плотности вероятности распределения параметров. Новое множество ©, которое используется для оценки параметров на следующем шаге, генерируется в соответствии со скорректированной плотно -стью вероятности распределения компонентов вектора параметров.
Оценивание параметров интерферометрических сигналов
Для реализации алгоритма необходима априорная информация о модели сигнала. Интерферомет-рический сигнал, как известно, можно представить дискретной последовательностью отсчетов
s(k) = B(k) + A(k) cos(0(k) + 5ф^)) + n(k) , (6)
где B(k) - фоновая составляющая; A(k) - амплитуда; n(k) - некоррелированный с сигналом белый шум с нулевым средним, распределенный по нормальному закону; Ф(к) - полная фаза сигнала; 5ф(k) - случайные флуктуации фазы. Полная фаза определяется как
Ф^) = ]Г 2nf (k')Дг ,
k'=0
где f (k') - частота; Дг - шаг дискретизации по независимой переменной.
С учетом модели интерферометрического сигнала (6) вектор параметров имеет вид
0 = (B,A, f,Ф)т ,
а векторные функции в уравнениях системы и наблюдения соответственно представляются в форме
h(0) = B + A cos Ф ,
f (0) = 0 + (0,0,0,27if Дг)T .
Плотность распределения вероятности каждого компонента вектора параметров на шаге k = 0 предполагается нормальной. Шумы системы и наблюдения считаются аддитивными, белыми и имеющими нормальное распределение с нулевым средним и некоторой дисперсией. Исходя из свойств типичных интерференционных картин и интерферометрических сигналов, можно заключить, что в процессе фильтрации закон распределения отклонений элементов вектора параметров останется нормальным, а параметры распределения (например, математическое ожидание и дисперсия) могут измениться.
Так как шум наблюдения в модели (6) можно считать распределенным по нормальному закону, для расчета вероятностей q (i) можно воспользоваться уравнением (4). Значения вероятностей q (i) обратны аргументу экспоненциальной функции, который прямо пропорционален квадрату евклидова расстояния между наблюдением s(k) и прогнозом сигнала s(k). Для уменьшения вычислительной сложности алго-
ритма при использовании данной модели целесообразно вместо максимизации функции (4) минимизировать невязку между наблюдением и прогнозом:
(7)
А =
h(ê(k ,i) )-s(k )
Использование выражения (7), однако, не позволяет вычислить вероятность q (г), необходимую для выбора векторов параметров по правилу (5). В этом случае можно осуществлять выбор М векторов (М < Ы), для которых невязка (7) минимальна.
Оценка вектора параметров на текущем шаге вычисляется как среднее значение выбранных в соответствии с правилом (5) элементов множества ©.
Рис. 2 иллюстрирует результаты динамического оценивания амплитуды и фазы интерферометри-ческого сигнала при помощи последовательного метода Монте-Карло. Параметры модели соответствуют реальным интерферометрическим сигналам, регистрируемым в системах корреляционной оптической когерентной томографии [17-21]. Фаза оцениваемого сигнала задавалась изменяющейся по нелинейному закону с отклонением ±2п рад на длине реализации. Количество отсчетов на периоде в модельном сигнале не превышает пяти. СКО аддитивного нормально распределенного белого шума составляет 10% от максимального значения амплитуды. Количество N генерируемых случайных векторов равняется 1000, количество М отбираемых - 300. Начальные дисперсии амплитуды и фазы соответственно равны 25 усл. ед. и п/4 рад.
к
и
S о
(D
S
и
(D
Se и
го
1
0,5 0 -0,5 -1
0
Номер дискретного отсчета k
а
«
« о <и g 25
H <D 20
У о S 15
О а п (0 а 10
о а « о <\> 5 0
ш H о я» н и s
w ю «
ч s
о S о
«
0 50 100 150 200 250 Номер дискретного отсчета k б
О "Л О "Л О "Л
^ ^ ^
<о <о <о <о <о <о ° Центральное значение интервала
ошибки оценки фазы, рад в
Рис. 2. Исходный сигнал и результат динамического оценивания амплитуды (а), ошибка оценки фазы (б) и гистограмма распределения ошибки оценки фазы на участке, содержащем полезную составляющую
сигнала (в)
Следует отметить, что поскольку последовательный метод Монте-Карло предполагает минимизацию невязки только для текущего дискретного отсчета к и не использует коэффициент усиления (как в алгоритмах калмановского типа), то исходная дисперсия шума сохраняются в процессе фильтрации. Указанный факт не всегда позволяет улучшить отношение сигнал/шум при использовании последовательного метода Монте-Карло, однако сохраняет устойчивость фильтра к нелинейностям изменения параметров сигнала. Так, для модельного сигнала с исходным СКО аддитивного шума 10% от максимального значения амплитуды, СКО ошибки оценки амплитуды составило 8% ее максимального значения. На рис. 2 показано, что на участке, содержащем полезную составляющую сигнала, ошибка фазы распределена по нормальному закону, а максимальное отклонение фазы равняется п/4 рад. В пределах всей реализации сигнала ошибка оценки фазы не превышает 3п/4 рад.
Анализ влияния входных параметров на результат работы алгоритма
Погрешность оценки параметров интерферометрического сигнала зависит от входных параметров алгоритма: плотности вероятности распределения параметров, количества генерируемых случайных век-
торов и количества отбираемых векторов, лучше всего удовлетворяющих текущим наблюдениям. В описанном выше примере плотность вероятности распределения отклонений параметров на каждом шаге является нормальной, и изменяются лишь параметры распределения. Ниже рассмотрено влияние соотношения количества генерируемых и отбираемых векторов на качество оценки параметров интерферо-метрических сигналов.
На рис. 3 представлены зависимости отношения сигнал/шум и СКО ошибки оценки амплитуды от количества отбираемых векторов. Количество генерируемых векторов фиксировано и равно 1000. Видно, что количество отбираемых векторов приблизительно одинаково для наиболее высокого отношения сигнал/шум и наименьшего значения ошибки амплитуды и составляет величину порядка 300. С увеличением этого параметра качество оценки падает.
На рис. 4 представлены зависимости отношения сигнал/шум восстановленного сигнала и СКО ошибки оценки амплитуды от количества генерируемых случайных векторов. Количество отбираемых векторов составляет 30% от количества генерируемых. Видно, что значительные изменения ошибки амплитуды и отношения сигнал/шум наблюдаются при значениях N меньше 150. Дальнейшее увеличение количества генерируемых векторов не влечет повышения качества оценивания.
s
S ч
К <и <ц ■
я g
О £ s w ю s
о
О «
О
30 25 20 15 10 5
0 200 400 600 800 Количество отбираемых векторов М
3
Ii
Т5 н
а и
и
S о
(D
S
и
(D
S
о и
10 8 6 4 2 0 -2
0 200 400 600 800 Количество выбираемых векторов М б
Рис. 3. Зависимости отношение сигнал/шум (а) и СКО оценки амплитуды (б) от количества выбираемых
случайных векторов параметров
м 12
о 10
о 8
а « и и
и ё tf 6
о в Ч
е и о и а И 4
И е В а н о о и и о 2
о о
и в
0
Г
к
о Я о
S
ю s S О
ч
о «
H
о
ч
!Г
ч а
25 20 15 10 5
200 400 600 800 Размер выборки N а
200 400 600 800 Размер выборки N б
Рис. 4. Зависимости отношение сигнал/шум (а) и СКО ошибки оценки амплитуды (б) от размера генерируемой выборки случайных векторов параметров
Заключение
В работе рассмотрено применение последовательного метода Монте-Карло для решения задачи динамического оценивания параметров интерферометрических сигналов.
Представленные результаты оценивания параметров модельных сигналов, характерных для систем корреляционной оптической когерентной томографии, показывают, что среднее квадратичное отклонение ошибки оценки амплитуды для типичных сигналов не превысило 8% от ее максимального значения при среднем квадратичном отклонении аддитивного шума 10%. В пределах участка, содержащего полезный сигнал, ошибка фазы распределена по нормальному закону. Ее максимальное отклонение не превысило п/4 рад на участке, содержащем полезную составляющую сигнала. Показано, что последовательный метод Монте-Карло является устойчивым к нелинейностям изменения параметров сигнала.
Проведен анализ влияния количества отбираемых векторов параметров на результаты оценивания. На основании результатов моделирования для рассмотренного класса сигналов рекомендуется отбирать 30% от количества генерируемых векторов. Увеличение количества генерируемых векторов более 150 не влечет за собой значительного улучшения качества получаемых оценок.
а
0
Литература
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. N 4. P. 514-521.
3. Van Kampen N. Stochastic Processes in Physics and Chemistry. North Holland, 1984. 464 p.
4. Степанов О.А. Основы теории оценивания с приложениями к задачам обработки навигационной информации. Ч. 1. Введение в теорию оценивания. СПб: Электроприбор, 2009. 496 с.
5. Kalman R.E. A new approach to linear filtering and prediction problems // Trans. ASME, J. Basic Eng. 1960. V. 82. P. 35-45.
6. Simon D. Using nonlinear Kalman filtering to estimate signals // Embedded Systems Design. 2006. V. 19. N 7. P. 38-53.
7. Simon D. Optimal state estimation: Kalman, H®, and Nonlinear Approaches. NY: John Wiley & Sons, Inc., 2006. 526 p.
8. Ярлыков М.С. Статистическая теория радионавигации. М: Радио и связь, 1985. 345 с.
9. Gurov I., Sheynihovich D. Interferometric data analysis based on Markov nonlinear filtering methodology // JOSA A. 2000. V. 17. N 1. P. 21-27.
10. 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.
11. Doucet A., de Freitas N., Gordon N. Sequential Monte Carlo methods in practice. NY: Springer-Verlag, 2001. 583 p.
12. Iba Y. Population Monte Carlo algorithms // Transactions of the Japanese Society for Artificial Intelligence. 2001. V. 16. N 2. P. 279-286.
13. Ristic B., Arulampalam S., Gordon N. Beyond the Kalman filter: Particle filters for tracking applications. Boston, Artech House, 2004. 318 p.
14. Isard M., Blake A. Contour tracking by stochastic propagation of conditional density // European Conference on Computer Vision. 1996. P. 343-356.
15. MacCormick J., Blake A. Probabilistic exclusion principle for tracking multiple objects // Proc. of IEEE International Conference on Computer Vision. 1999. P. 572-578.
16. Del Moral P. Measure valued processes and interacting particle systems. Application to nonlinear filtering problems // Annals of Applied Probability. 1998. V. 8. N 2. P. 438-495.
17. Гуров И.П. Оптическая когерентная томография: принципы, проблемы и перспективы / В сб.: Проблемы когерентной и нелинейной оптики / Под ред. И.П. Гуров, С.А. Козлов. СПб: СПбГУ ИТМО, 2004. С. 6-30.
18. Fercher A. Optical coherence tomography // Journal of Biomedical Optics. 1996. V. 1. N2. P. 157-173.
19. Гуров И.П., Жукова Е.В., Маргарянц Н.Б. Исследование внутренней микроструктуры материалов методом оптической когерентной микроскопии с перестраиваемой длиной // Научно-технический вестник информационных технологий, механики и оптики. 2012. № 3 (79). С. 40-45.
20. Гуров И.П., Жукова Е.В., Левшина А.В. Применение метода оптической когерентной томографии для изучения предметов искусства, выполненных в технике интарсии // Научно-технический вестник информационных технологий, механики и оптики. 2012. № 3 (79). С. 55-59.
21. Волынский М.А., Воробьева Е.А., Гуров И.П., Маргарянц Н.Б. Бесконтактный контроль микрообъектов методами интерферометрии малой когерентности и оптической когерентной томографии // Изв. вузов. Приборостроение. 2011. Т. 54. № 2. С. 75-82.
Волынский Максим - кандидат технических наук, доцент, Университет ИТМО, Санкт-
Александрович Петербург, Россия, [email protected]
Гуров Игорь Петрович - зав. кафедрой, доктор технических наук, профессор, Университет
ИТМО, Санкт-Петербург, Россия, [email protected]
Ермолаев Петр Андреевич - студент, Университет ИТМО, Санкт-Петербург, Россия,
Скаков Павел Сергеевич - ассистент, Университет ИТМО, Санкт-Петербург, Россия,
Maxim A. Volynsky - associate professor, PhD, ITMO University, Saint Petersburg, Russia,
Igor P. Gurov - Department head, D.Sc., Professor, ITMO University, Saint Petersburg,
Russia, [email protected]
Peter A. Ermolaev - student, ITMO University, Saint Petersburg, Russia,
Pavel S. Skakov - assistant, ITMO University, Saint Petersburg, Russia, [email protected]
Принято к печати 26.03.2014 Accepted 26.03.2014
