CC BY
41
УДК 681.787:[519.245+519.6]
АНАЛИЗ ВЫЧИСЛИТЕЛЬНОЙ СЛОЖНОСТИ РЕКУРРЕНТНЫХ АЛГОРИТМОВ ОБРАБОТКИ ДАННЫХ В ОПТИЧЕСКОЙ КОГЕРЕНТНОЙ
ТОМОГРАФИИ М.А. Волынскийа, И.П. Гурова, П.А. Ермолаева, П.С. Скакова а Университет ИТМО, Санкт-Петербург, 197101, Российская Федерация, maxim.volynsky@gmail.com Аннотация. Рассмотрены основные принципы представления сигналов в оптической когерентной томографии с использованием формализма теории динамических систем; проведен сравнительный анализ вычислительной сложности алгоритмов динамического оценивания параметров сигналов в оптической когерентной томографии, таких как расширенный фильтр Калмана и последовательный метод Монте-Карло. Показано, что вычислительная сложность обработки одного отсчета сигнала при помощи расширенного фильтра Калмана полиномиально возрастает в зависимости от размера вектора параметров и вектора наблюдения, а сложность обработки отсчета сигнала последовательным методом Монте-Карло линейно зависит как от размеров вектора параметров и вектора наблюдения, так и от количества генерируемых случайных векторов. Приведены экспериментальные результаты оценивания времени обработки тестового сигнала при использовании каждого из алгоритмов. Показано, что время обработки сигнала, содержащего 500 дискретных отсчетов, при помощи расширенного фильтра Калмана в случае простейшей модели скалярного сигнала составляет примерно 0,1 с и возрастает при усложнении модели в несколько раз. Время обработки аналогичного сигнала при помощи последовательного метода Монте-Карло с использованием аналогичной простейшей модели и при фиксированном количестве генерируемых векторов составляет 0,7 с и при усложнении модели возрастает незначительно, примерно в 1,5 раза. Полученные результаты могут быть использованы при оценке ожидаемого времени обработки данных с помощью рекуррентных алгоритмов динамического оценивания параметров в системах оптической когерентной томографии.
Ключевые слова: оптическая когерентная томография, обработка интерферометрических сигналов, вычислительная сложность, расширенный фильтр Калмана, последовательный метод Монте-Карло.
Благодарности. Работа выполнена при финансовой поддержке Министерства образования и науки Российской Федерации.
^MPUTATIONAL COMPLEXITY ANALYSIS OF RECURRENT DATA PROCESSING ALGORITHMS IN OPTICAL COHERENCE TOMOGRAPHY M.A. Volynsky% I.P. Gurov% Р.А. Ermolaev% P.S. Skakov:i
а ITMO University, Saint Petersburg, 197101, Russian Federation, maxim.volynsky@gmail.com
Abstract. The paper deals with the basic principles of signals representation in optical coherence tomography with the usage of dynamic systems theory formalism. Computational complexity of algorithms for dynamic estimation of signals parameters is analyzed, such as extended Kalman filter and sequential Monte-Carlo method. It is shown that processing time of one discrete-time sample of the signal by extended Kalman filter increases polynomially with sizes of parameters vector and observation vector. Processing time of one discrete-time sample of the signal by sequential Monte-Carlo method depends linearly both on sizes of parameters vector and observation vector, and on the number of generating random vectors. Experimental results of processing time measurement by each algorithm are described. It is shown that processing time of the signal containing 500 discrete-time samples by extended Kalman filter in the case of the simplest model is approximately equal to 0.1 seconds and increases several times with complication of the model. Processing time of the same signal by sequential Monte-Carlo methods with fixed number of generated random vectors is equal to 0.7 seconds and slightly increases with complication of the model, approximately by 1.5 times. Obtained results may be used for estimation of expected data processing time by recurrent dynamic estimation algorithms in optical coherence tomography systems. Keywords: optical coherence tomography, interferometric signals processing, computational complexity, extended Kalman filter, sequential Monte-Carlo method.
Acknowledgements. The work has been carried out with the financial support of the Ministry of Education and Science of the Russian Federation.
Введение
Оптическая когерентная томография (ОКТ) является одним из перспективных методов исследования микроструктуры объектов с высоким разрешением и широко применяется в различных областях медицины и биологии [1-6]. Системы ОКТ подразделяются по типу регистрируемых сигналов на корреляционные [2] и спектральные [3]. В первом случае на выходе системы регистрируется интерферометриче-ский сигнал малой когерентности, а во втором случае - формируются полосы равного хроматического порядка. В этих сигналах содержится информация о внутренней микроструктуре объекта, для извлечения которой требуются специальные алгоритмы обработки сигналов [7].
Критическими требованиями к таким алгоритмам, помимо качества обработки, являются высокое быстродействие и малое количество требуемой памяти. Эти характеристики могут быть оценены в рамках теории вычислительной сложности [8, 9].
Традиционно для обработки сигналов в ОКТ используются алгоритмы на основе преобразования Фурье, однако их применение в ряде случаев сопряжено с недостатками, например, в системах, работающих в реальном времени, так как применение дискретного преобразования Фурье требует наличия полной реализации обрабатываемого сигнала.
Преимуществом рекуррентных алгоритмов динамического оценивания [7, 10], основанных на формализме стохастических дифференциальных уравнений [7], является использование априорной информации о процессе формирования сигнала и статистических характеристиках шума.
Наиболее распространенным алгоритмом динамического оценивания параметров является линейный фильтр Калмана [11, 12]. Этот алгоритм оптимален с точки зрения минимума средней квадратичной ошибки оценки. В системах ОКТ наблюдаемое значение сигнала зависит от параметров сигнала нелинейно, что ограничивает возможность использования алгоритма линейной фильтрации для оценивания параметров таких систем. Эта проблема может быть решена при помощи расширенного фильтра Калмана (РФК) [12, 13], в котором используется линеаризация уравнений динамической системы и (или) наблюдения путем разложения их в ряд Тейлора по параметрам.
Альтернативным подходом к оцениванию параметров нелинейных динамических систем является последовательный метод Монте-Карло (ПММК) [12, 14-17], работа которого основана на статистической аппроксимации плотности вероятности распределения параметров. ПММК является более перспективным методом обработки по сравнению с РФК, так как в условиях априорной неопределенности модели сигнала предоставляет больше возможностей по адаптации алгоритма путем изменения моделируемых плотностей вероятностей распределения параметров и более устойчив к неточности задания начальных условий [18].
В настоящей работе проведен сравнительный анализ вычислительной сложности РФК и ПММК на примере задачи обработки сигналов в корреляционной ОКТ. Даны рекомендации по использованию рассматриваемых алгоритмов.
Модель интерферометрического сигнала
Сигнал в корреляционной ОКТ может быть представлен параметрически в виде дискретной последовательности отсчетов [7]:
s(k) = B(k) + A(k) cos^(k) + 5ф(£)) + n(k), (1)
где k = 0..K-1 - номер дискретного отсчета; B(k) - фоновая составляющая; A(k)- амплитуда; n(k) - не коррелированный с сигналом белый шум, распределенный по нормальному закону с нулевым средним; 5ф^) - случайные флуктуации фазы; <^(k) - полная фаза сигнала, определяемая как
Ф^) = £ 2nf (k')Az ,
k '=0
где f (k') - частота; Az - шаг дискретизации.
Набор параметров в уравнении (1) можно записать в виде вектора параметров
0 = (B, A, f, Ф)т . (2)
Для применения рекуррентных алгоритмов динамического оценивания параметров требуется записать модель формирования интерферометрического сигнала (1), используя формализм теории динамических систем. Динамическая система задается в общем виде уравнениями системы и уравнением наблюдения:
0(k ) = f (0(k -1), w( k ) ), (3)
s(k ) = h (0(k ), n(k ) ), (4)
где f (0) и h(0) - известные нелинейные векторные функции, w и n - случайные шумы системы и наблюдения, а вектор 0 имеет известную плотность вероятности распределения параметров на шаге k = 0. С учетом (1) и (2) для случая неизменных фоновой составляющей, амплитуды и частоты в пределах шага дискретизации векторные функции в уравнениях (3) и (4) можно записать как
h(0) = B + A cos Ф ,
f (0) = 0 + (0,0,0,2rcf Az)T .
Динамическая обработка сигналов в ОКТ заключается в получении оценки вектора параметров (2) для каждого отсчета сигнала (1).
В зависимости от решаемой задачи и способа формирования сигналов для их описания могут быть использованы другие модели, в том числе модели с векторным наблюдением (когда результат вычисления функции h(0) является вектором). Примером может служить оценивание параметров сигнала по нескольким предшествующим отсчетам по последовательности видеокадров, что позволяет повысить точность оценки частоты сигнала.
Краткое описание алгоритмов
РФК инициализируется заданием нелинейных функций h(0) и f (0), начальных значений параметров на шаге k = 0 и ковариационных матриц шума системы и шума наблюдения. Начальные значения
параметров и ковариационные матрицы могут быть найдены путем предварительного анализа интерфе-рометрического сигнала.
Работа фильтра включает два этапа - предсказание и коррекцию по результатам наблюдения. На первом этапе происходит предсказание значений ковариационной матрицы ошибок и экстраполяция значений вектора параметров с учетом функции f (0):
0(к) = f (0(к -1)).
На этапе коррекции осуществляется корректировка предсказанного значения в(к) с учетом наблюдений на текущем шаге при помощи уравнения
в(к) = в(к)+P(k мк) - ы;в(к))],
где s(k) - вектор наблюдения на текущем шаге; Г (к) - матричный коэффициент усиления фильтра, определяющий вклад невязки наблюдения и предсказания в оценку вектора параметров. Коэффициент усиления вычисляется как
Р(к) = Я р. Н (к )т (Н (к р. Н(к )т + Я я )-1, где Ярг - предсказание ковариационной матрицы ошибок; - ковариационная матрица шума наблюдения, а Н(к) - матрица первых производных функции ^0) по компонентам вектора параметров. Подробнее особенности динамического оценивания параметров при помощи РФК и его модификаций рассмотрены в работах [7, 12, 19].
Инициализация алгоритма, реализующего ПММК, осуществляется не только заданием функций ^0), f (0) и начальных значений параметров на шаге к = 0, но и заданием количества генерируемых случайных векторов, правила отбора векторов (например, пороговой вероятности) и начальной функции плотности распределения параметров р(0(О)). В простейшем случае задаются моменты нормального распределения для каждого параметра.
Обработка одного отсчета сигнала при помощи ПММК может быть разделена на четыре этапа. На первом этапе в соответствии с известной плотностью вероятности распределения параметров происходит генерация набора из N случайных векторов параметров. На втором этапе для каждого сгенерированного вектора осуществляется экстраполяция значений параметров в соответствии с функцией f (0). На третьем этапе осуществляется отбор векторов параметров, которые наилучшим образом удовлетворяют наблюдениям. Критерием отбора может служить, например, пороговое значение некоторой меры совпадения [10, 14, 18]. На четвертом этапе работы ПММК оценивается искомый вектор параметров (например, как среднее значений отобранных векторов) и корректируется функция плотности вероятности распределения параметров в соответствии с распределением отобранных векторов. Подробнее процесс обработки квазигармонических сигналов при помощи ПММК рассмотрен в работе [18].
Теоретическая оценка вычислительной сложности алгоритмов
Мерой вычислительной сложности считается количество элементарных операций, необходимых для выполнения алгоритма на ЭВМ. При оценке сложности пользуются асимптотическими оценками, определяющими порядок роста количества элементарных операций в зависимости от размера входных данных. Вычислительная сложность алгоритмов динамического оценивания зависит от количества отсчетов К в обрабатываемом сигнале, размеров вектора параметров и вектора наблюдения. Обозначим размеры вектора параметров и вектора наблюдения как р х 1 и q х 1 соответственно. Тогда ковариационная матрица шума системы имеет размер р х р, шума наблюдения - q х q, ошибок - р х р. Размеры матриц производных функций ^0) и f (0) равны соответственно р х р и q х р, а размер коэффициента усиления
-р х q [12].
При обработке сигналов при помощи РФК наибольших временных затрат требуют расчеты, связанные с умножением и обращением матриц [12]. Простейший алгоритм умножения матриц требует 0(т3) операций, где т - размер умножаемых матриц, следовательно, сложность уравнений РФК, в которых используется умножение матриц, может быть оценена как 0(р3 + q3). Существуют более эффективные алгоритмы, применяемые для умножения больших матриц, однако их использование для матриц малого размера нецелесообразно.
Известно, что обращение матрицы эквивалентно умножению двух матриц [8, 9] и также может быть выполнено в простейшем случае за 0(т3) операций. Из этого следует, что вычислительная сложность расчета уравнений, содержащих обращение матриц, может быть оценена аналогично как 0(р3 + q3).
Так как функции f (0) и ^0) известны и не меняются в процессе динамического оценивания, уравнения для вычисления значений матриц их производных могут быть найдены априорно и заданы при инициализации фильтра. Сложность вычисления элементарных функций на ЭВМ оценивается как 0(1), так как зависит от разрядности используемой архитектуры и требуемой точности, а не от размера вход-
ных данных. Так, сложность вычисления функций f (0) и Ь(9) может быть оценена соответственно как O(p2) и O(pq), а сложность расчета их производных - как O(p3) и O(p2q).
Общая вычислительная сложность РФК применительно к обработке сигнала из K отсчетов составляет O(K(p3 + q3)). Время обработки отсчета сигнала увеличивается полиномиально в зависимости от размеров векторов параметров и наблюдения. В работе [20] описаны методы оптимизации процесса динамического оценивания параметров при помощи РФК.
На вычислительную сложность ПММК влияют не только размеры векторов параметров и наблюдения, но и количество генерируемых случайных векторов N. Сложность генерации одного случайного числа составляет O(1), следовательно, сложность первого этапа работы алгоритма - O(pN). Так как вычисление функций h(0) и f(0) осуществляется для каждого из сгенерированных векторов, сложность второго этапа оценивается как O(N(p2 + pq)).
Сложность третьего этапа зависит от метода отбора векторов. Когда критерием отбора является порог меры совпадения [10, 14, 18], сложность этапа равна O(N). Альтернативным подходом является отбор фиксированного количества векторов, лучше всего удовлетворяющих наблюдению, что требует выполнения сортировки случайных векторов. В этом случае при использовании методов быстрой сортировки сложность третьего этапа может быть оценена как O(NlogN).
Плотность распределения параметров корректируется на четвертом этапе независимо для каждого параметра и имеет сложность O(p). Общая вычислительная сложность ПММК в случае отбора векторов по критерию порога меры совпадения может быть оценена как O(N(p2 + pq)), а в случае отбора фиксированного количества векторов - как O(N(p2 + pq + logN)).
Экспериментальное сравнение времени работы алгоритмов
Экспериментальная оценка времени работы алгоритмов проводилась с использованием процессора Intel Core i7-4500U с номинальной тактовой частотой 1,8 ГГц. Тестовый сигнал содержал K = 500 отсчетов. В таблице представлено экспериментальное время работы алгоритмов в зависимости от размеров векторов параметров и наблюдения. Количество генерируемых векторов в ПММК N = 150.
Размер вектора параметров Размер вектора наблюдения Время работы РФК, с Время работы ПММК, с
2 1 0,1 0,7
3 1 0,2 0,7
4 1 0,2 0,7
3 8 0,3 0,8
4 8 0,4 0,9
5 8 0,5 1,0
6 8 0,8 1,0
Таблица. Экспериментальные результаты оценки времени работы РФК и ПММК в зависимости от размеров векторов параметров и наблюдения в используемой модели
Видно, что время работы РФК растет быстрее, чем время работы ПММК с увеличением размеров векторов параметров и наблюдения. Так, при больших размерах этих векторов обработка сигнала при помощи ПММК требует незначительно больше времени, чем при помощи РФК. На рисунке представлен экспериментальный график зависимости времени работы ПММК от количества генерируемых случайных векторов. Видно, что эта зависимость близка к линейной.
л
4 л и
и
5 о
5
н о ю л а ю о
(D
а m
100 200 300 400
Количество генерируемых векторов N
500
Рисунок. Зависимость времени работы ПММК от количества генерируемых векторов параметров N
4
3
2
1
0
При создании систем ОКТ необходимо учитывать и объемы памяти, используемой алгоритмами обработки данных, однако этот вопрос требует отдельного рассмотрения, выходящего за рамки данной работы.
Заключение
В работе проведен анализ вычислительной сложности расширенного фильтра Калмана и последовательного метода Монте-Карло применительно к обработке данных в оптической когерентной томографии. Показано, что количество элементарных операций, требуемых для работы алгоритмов, растет линейно в зависимости от количества отсчетов в обрабатываемом сигнале.
Так как в процессе работы расширенного фильтра Калмана используются операции умножения и обращения матриц, вычислительная сложность обработки одного отсчета сигнала полиномиально зависит от размеров вектора параметров и вектора наблюдения. Время обработки отсчета сигнала последовательным методом Монте-Карло линейно зависит как от размеров вектора параметров и вектора наблюдения, так и от количества генерируемых случайных векторов, которое обычно значительно превышает размеры этих векторов и составляет несколько сотен, вследствие чего время работы последовательного метода Монте-Карло выше, чем время работы расширенного фильтра Калмана.
При малых размерах вектора параметров и вектора наблюдения и при наличии достаточно точной априорной информации о начальных значениях параметров и характеристиках шума рекомендуется использовать расширенный фильтр Калмана, так как он обеспечивает более высокое быстродействие. Последовательный метод Монте-Карло допускает существенную априорную неопределенность модели сигнала и обеспечивает большую устойчивость к шуму и неизвестным начальным значениям параметров, однако требует более длительного времени работы, чем расширенный фильтр Калмана при малых размерах вектора параметров и вектора наблюдения. При больших размерах вектора параметров системы и вектора наблюдения рекомендуется использовать последовательный метод Монте-Карло, так как это требует незначительно большего времени, чем использование расширенного фильтра Калмана.
Литература
1. Deck L., de Groot P. High-speed non-contact profiler based on scanning white light interferometry // Applied Optics. 1994. V. 33. P. 7334-7338.
2. Гуров И.П. Оптическая когерентная томография: принципы, проблемы и перспективы. В кн.: Проблемы когерентной и нелинейной оптики / Под ред. И.П. Гурова, С. А. Козлова. СПб.: СПбГУ ИТМО,
2004. С. 6-30.
3. Fercher A. Optical coherence tomography // Journal of Biomedical Optics. 1996. V. 1. N 2. P. 157-173.
4. Волынский М.А., Воробьева Е.А., Гуров И.П., Маргарянц Н.Б. Бесконтактный контроль микрообъектов методами интерферометрии малой когерентности и оптической когерентной томографии // Изв. вузов. Приборостроение. 2011. Т. 54. № 2. С. 75-82.
5. Wyant J.C. Interferometric optical metrology: basic principles and new systems // Laser Focus with Fiveroptic Technology. 1982. V. 18. N 5. P. 65-71.
6. Alarousu E., Krehut L., Prykari T., Myllyla R. Study on the use of optical coherence tomography in measurements of paper properties // Measurement Science and Technology. 2005. V. 16. N 5. P. 1131-1137.
7. 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.
8. Грин Д., Кнут Д. Математические методы анализа алгоритмов. М.: Мир, 1987. 120 с.
9. Кормен Т., Лейзерсон Ч., Ривест Р., Штайн К. Алгоритмы: построение и анализ. 2-е изд. М.: Вильямс,
2005. 1296 с.
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. Kalman R.E. A new approach to linear filtering and prediction problems // Trans. ASME, J. Basic Eng. 1960. V. 82. P. 35-45.
12. Simon D. Optimal state estimation: Kalman, H®, and nonlinear approaches. NY: John Wiley & Sons, Inc.,
2006. 526 p.
13. Simon D. Using nonlinear Kalman filtering to estimate signals // Embedded Systems Design. 2006. V. 19. N 7. P. 38-53.
14. Doucet A., de Freitas N., Gordon N. Sequential Monte Carlo methods in practice. NY: Springer-Verlag, 2001. 583 p.
15. Gordon N.J., Salmond D.J., Smith A.F.M. Novel approach to nonlinear/non-Gaussian Bayesian state estimation // IEE Proceedings Part F: Radar and Signal Processing. 1993. V. 140. N 2. P. 107-113.
16. Ristic B., Arulampalam S., Gordon N. Beyond the Kalman filter: particle filters for tracking applications. Boston: Artech House, 2004. 318 p.
17. Gilks W., Berzuini C. Following a moving target - Monte Carlo inference for dynamic Bayesian models // Journal of the Royal Statistical Society. Series B: Statistical Maethodology. 2001. V. 63. N 1. P. 127-146.
18. Волынский М.А., Гуров И.П., Ермолаев П.А., Скаков П.С. Динамическое оценивание параметров ин-терферометрических сигналов на основе последовательного метода Монте-Карло // Научно-технический вестник информационных технологий, механики и оптики. 2014. № 3 (91). С. 18-23.
19. Ермолаев П. А. Динамическое оценивание параметров интерферометрических сигналов методом расширенной фильтрации Калмана второго порядка // Научно-технический вестник информационных технологий, механики и оптики. 2014. № 2 (90). С. 17-22.
20. Гуров И.П., Захаров А.С., Таратин М.А. Анализ и оптимизация вычислительного процесса нелинейной дискретной фильтрации Калмана // Изв. вузов. Приборостроение. 2004. Т. 47. № 8. С. 42-48.
Волынский Максим Александрович Гуров Игорь Петрович
Ермолаев Петр Андреевич Скаков Павел Сергеевич
кандидат технических наук, доцент, Университет ИТМО, Санкт-Петербург, 197101, Российская Федерация, maxim.volynsky@gmail.com доктор технических наук, профессор, заведующий кафедрой, Университет ИТМО, Санкт-Петербург, 197101, Российская Федерация, gurov@mail.ifmo.ru
студент, Университет ИТМО, Санкт-Петербург, 197101, Российская Федерация, Petr-ermolaev@hotmail.com
ассистент, Университет ИТМО, Санкт-Петербург, 197101, Российская Федерация, pavelsx@gmail.com
Maxim A. Volynsky Igor P. Gurov Peter A. Ermolaev Pavel S. Skakov
PhD, Associate professor, ITMO University, Saint Petersburg, 197101,
Russian Federation, maxim.volynsky@gmail.com
D.Sc., Professor, Department head, ITMO University, Saint Petersburg,
197101, Russian Federation, gurov@mail.ifmo.ru
student, ITMO University, Saint Petersburg, 197101, Russian Federation,
Petr-ermolaev@hotmail.com
assistant, ITMO University, Saint Petersburg, 197101, Russian Federation, pavelsx@gmail.com
Принято к печати 10.10.14 Accepted 10.10.14
