CC BY
52
НЕЛИНЕЙНАЯ СТОХАСТИЧЕСКАЯ ФИЛЬТРАЦИЯ СИГНАЛОВ В ИНТЕРФЕРОМЕТРАХ С ЧАСТИЧНО КОГЕРЕНТНЫМ ОСВЕЩЕНИЕМ
М.А. Волынский, И.П. Гуров, А.С. Захаров
Предложены модели интерференционных сигналов малой когерентности, получаемых при анализе многослойных сред. Исследованы характеристики разрешающей способности метода нелинейной дискретной фильтрации Калмана применительно к интерференционным сигналам в оптической когерентной томографии. Проведено сравнение результатов, полученных с использованием различных моделей сигналов.
Введение
Бесконтактные методы неразрушающего контроля объектов имеют большое значение при решении многих научно-технических задач. Методы оптической когерентной томографии (ОКТ) обеспечивают наиболее высокое разрешение по сравнению с другими известными методами томографии и получили широкое распространение в биомедицинских исследованиях, исследованиях микроструктуры различных материалов и в других областях.
Метод ОКТ основан на формировании интерферометрических сигналов малой когерентности, которые содержат полезную информацию о свойствах исследуемых случайно-неоднородных сред. Для извлечения полезной информации требуются быстродействующие помехоустойчивые методы и алгоритмы обработки сигналов ОКТ.
Обработка сигналов основывается на использовании априорной информации о свойствах полезной составляющей сигнала и помех [1]. Эту информацию можно представить в виде математической модели, учитывающей особенности процесса формирования сигналов. При описании и анализе интерферометрических сигналов в ОКТ особый интерес представляет использование стохастических моделей и методов стохастической дискретной фильтрации параметров сигналов ОКТ. В частности, для решения поставленных задач оптимальным методом фильтрации (квазиоптимальным в нелинейном случае) является метод Калмана-Бьюси [2, 3].
В работе рассматриваются модели интерференционных сигналов малой когерентности, получаемые при исследованиях многослойных сред, и результаты обработки ин-терферометрических сигналов с применением расширенного фильтра Калмана.
Стохастические модели исследуемой среды и интерферометрических сигналов
При рассмотрении двумерного сечения неоднородной среды каждая граница двух слоев среды может быть представлена в виде функции 2. (х). Изменение рельефа каждого слоя можно представить как решение нестационарного уравнения Ланжевена [4, 5]:
й2 (X)
= -«2. (X) + Ж (X), (1)
ах
где ж (х) - формирующий шум, / - номер слоя исследуемой структуры (слои считаются изотропными). Постоянная а задает динамику случайного изменения высоты рельефа слоя. Функция ж(х) такова, что {'ж.(х)^ = 0 и (ж.(х)жк(х= (0.12)8.к8(х - х'), где О. - спектральная плотность .-ой компоненты формирующего шума, 3]к - символ Кронекера, 5(х - х') - дельта-функция Дирака.
Для каждой точки поверхности проведено моделирование интерференционного сигнала, получаемого в соответствующей точке плоскости наблюдения при изменении положения опорного отражателя в интерферометре малой когерентности в форме
5(X, 2) = s0 (X, 2) + ^ (X, 2) С08(Ф(2) + р(2)) + п(X, 2) , (2)
где 50 (х, 2) - фоновая составляющая сигнала, (х, 2) - огибающая, 2 - координата по глубине слоя, Ф (2) - фаза сигнала, 5р( 2) - равномерно распределенный фазовый шум, вызванный неравномерностью движения опорного отражателя, п (х, 2) - случайный белый гауссовский шум наблюдения [5]. Фаза сигнала равна
Ф (2) = 2п/ 2 +р, (3)
где /0 - несущая частота, которой соответствует центральная длина волны = с / /0, р - начальная фаза. Огибающую целесообразно определить в виде суммы гауссовых кривых
(2 - (х))2
,(х, z) = Х A exP
2
О
(4)
где Aj - постоянные, а - параметр, зависящий от свойств источника излучения и имеющий смысл длины когерентности [4, 6]. Коэффициенты Aj учитывают эффекты, связанные с отражением от исследуемого слоя и от предыдущих слоев, а также поглощение в слоях.
Полагая процесс излучения стационарным и эргодическим, можно использовать теорему Винера-Хинчина, согласно которой корреляционная функция является преобразованием Фурье спектральной плотности мощности процесса. В свою очередь, преобразованием Фурье функции гауссовой формы, как известно, является функция такой же формы. В работе [7] показано, что в случае источника с гауссовым спектром гауссова огибающая интерференционного сигнала имеет полуширину, примерно равную длине когерентности излучения источника.
При дискретном перемещении опорного отражателя «текущая» координата z = kAz, где к - номер дискретного отсчета сигнала, Az - шаг смещения опорного отражателя. При этом модель сигнала определяется формулой
s (х, к) = s0( х, z) + sm (х, к )cos(2^ fk Az + р + 5р( z)) + n (х, к), (5)
где дискретная последовательность значений амплитуды сигнала имеет вид
2
X х к) = Х A exp
(к Az - zt (х))2
2
О
(6)
Для исследования характеристик метода дискретной фильтрации Калмана были смоделированы интерференционные сигналы двух типов: с одним максимумом и сигналы, полученные при отражении измерительной волны от двуслойной среды как простейшего примера многослойной, т.е. сигналы с двумя истинными максимумами. Ко всем сигналам добавлялись фоновая составляющая и случайные флуктуации фазы, вызываемые вибрациями опорного отражателя при перемещении и погрешностью перемещений.
При моделировании сигналов использовались следующие параметры: центральная длина волны 0,8 мкм, шаг смещения опорного отражателя = 0,016 мкм, шаг дискретизации в боковом направлении Ах = 1 мкм, длина когерентности излучения источника 6 мкм.
Примеры сигналов представлены на рис. 1.
ч и
к н о
ей Ч ей X и К о
Ч К X <и
й X СП
2,5
1,5
0,5
-0,5
ч и
к н о
ч й X и
к
о ч к
X
<и р
й X СП
8
16 24
32
40 48
56
Длина пути измерительной волны, мкм
а)
2,5 2 1,5 1
0,5 0
-0,5
8
16 24
32
40
48
56
Длина пути измерительной волны, мкм
б)
Рис. 1. Примеры исходных интерференционных сигналов с фоновой составляющей (а - с одним максимумом; б - с двумя максимумами)
Пример отклонений фазы от линейного закона, распределенных по равновероятному закону, показан на рис. 2. При моделировании флуктуации фазы составляли от 0 до 0,5 рад.
0 8 16 24 32 40 48 56 Длина пути измерительной волны, мкм
Рис. 2. Флуктуации фазы сигнала
При обработке сигналов необходимо устранение шумов для корректного определения положений максимумов огибающей сигналов. Ниже рассматриваются особенно-
2
1
0
0
0
сти использования нелинейной фильтрации Калмана для обработки интерферометриче-ских сигналов.
Алгоритм дискретной нелинейной фильтрации Калмана
Дискретный линейный фильтр Калмана рассмотрен, например, в работах [2-4, 8]. Приведем основные соотношения, определяющие алгоритм нелинейной обработки ин-терферометрических сигналов в рассматриваемом случае.
Для дискретного нелинейного фильтра Калмана уравнение наблюдения имеет
вид:
= Ь(0(к)) + п(к), (7)
где И(0(к)) - полезная составляющая векторного сигнала, п(к) - шум наблюдения. Уравнение системы определяется в форме
0(к+1) = 0(к) + Щк)) & + w(k), (8)
где 0(к) - предыдущее состояние системы, w(k) - шум системы, Д2 - шаг дискретизации. Функция Г(0(к)) описывает эволюцию системы. Алгоритм нелинейной фильтрации рассмотрен более подробно в [4, 9].
Максимум огибающей интерференционного сигнала можно выделить, включив в вектор параметров амплитуду и фазу сигнала. При этом вектор параметров примет вид: 0 = [50, 5т, Ф]Т. (9)
Функция наблюдения представляется как
И(0) = 50 + Эт С0Б(Ф). (10)
Матрица перехода В для вектора параметров (9) в простейшем случае определяется как
В
' 1 0 0Л
0 1 0 0 0 1
(11)
Vй и V
В качестве компонента вектора параметров возможно использование непосредственно положения максимума огибающей. Функция 2(х) была включена в вектор параметров, например, в работе [10], при этом функция наблюдения имела вид:
и(0) = 50 + А ехр[-(кД2- 2(х))2/а2] сов(ф), (12)
где к - номер дискретного отсчета, Д2 - шаг дискретизации, однако использование такого вектора параметров не было достаточно исследовано. Ниже приводятся результаты использования вектора параметров вида
0 = [50, Ф, А, а, 2 ]Т. (13)
Отметим, что параметры А и а при фильтрации практически не корректируются ввиду того, что в реальных интерферометрических системах они квазиконстантны.
Аналогично (11), матрица перехода В для вектора параметров (13) в простейшем случае определяется как
0 0 0 0Л
В
0 10 0 0 0 0 10 0 0 0 0 1 0
V0 0 0 0 1 у
(14)
Из (14) видно, что значение фазы сигнала на шаге к+1 предсказывается как значение фазы на шаге к, при этом поправка на 2п/0Ак не производится.
Ниже представлены результаты восстановления параметров интерференционных сигналов с применением моделей (9)—(11) и (12)-(14), соответственно.
Анализ полученных результатов
На рис. 3 представлены исходный сигнал без фоновой составляющей и результат восстановления огибающей сигнала с применением модели (9)—(11).
0 8 16 24 32 40 48 56 Длина пути измерительной волны, мкм
Рис. 3. Интерференционный сигнал в условиях шума и оценка огибающей
Среднее квадратическое отклонение (СКО) оценки положения максимума составило 0,5 мкм.
Рассмотрим результаты, получаемые при использовании вектора параметров, включающего непосредственно координату максимума огибающей, т. е. координату по глубине исследуемого слоя.
На рис. 4 показана оценка глубины расположения исследуемого слоя, полученная при обработке сигнала рис. 3 с максимумом в точке 32 мкм.
й
№ е
я 3
н е 2 а
я ю
о и
л и
о о
п а
а
к му
н м
е я и о
О к
а
ма
40
30
й 20
10
8 16 24 32 40 48 56 Длина пути измерительной волны, мкм
0
0
Рис. 4. Оценка положения максимума огибающей интерференционного сигнала
Видно, что преимущество использованного подхода заключается в возможности предсказания положения максимума огибающей раньше его фактического появления (для рассматриваемого примера в точке 20 мкм), которое обеспечивается использованием заданной модели гауссовской огибающей и априорной информации о длине когерентности излучения источника.
Из рис. 4 видно, что оценка положения максимума имеет участки с относительно стабильным отклонением оценки от среднего значения на участке. Установив некоторый критерий, согласно которому такие «полочки» можно считать областью между границами отражающих слоев образца, можно выделить координаты границ слоев среды. Следует отметить, что при таком рассмотрении число слоев динамически опреде-
ляется в процессе фильтрации и не требует учета в модели. В эксперименте для идентификации границ отражающих слоев использовались следующие критерии: минимальная длина «полочки» и максимально возможное отклонение от среднего значения оценки в рамках указанного квазиоднородного участка. Выбор конкретных значений этих параметров определяется требованиями к разрешающей способности метода, длиной когерентности излучения источника и априорной информацией (при ее наличии) о структуре исследуемой среды.
На рис. 5 показана оценка глубины расположения исследуемого слоя, полученная при обработке сигнала с двумя максимумами в точках 30 и 40 мкм соответственно. При минимальной длине квазистабильного участка 6 мкм и максимально допустимом отклонении от среднего значения оценки в рамках одного участка, равном 2 мкм, ошибка определения глубины слоев составила 0,94 мкм и 1,30 мкм для первого и второго слоев, соответственно.
« и
§ 5 = §
|Ц ей
« Ю
8 3
ч о с
ей И К <и
а О
50 40 30 20 10 0
0 8 16 24 32 40 48 56 Длина пути измерительной волны, мкм
Рис. 5. Оценка положений максимумов огибающей интерференционного сигнала
Отметим, что для сигнала с одним максимумом (рис. 4) ошибка оценки слоя составила примерно 0,7 мкм.
При описанном подходе СКО оценки толщины слоев исследуемой структуры, усредненное по всем точкам, составило 3,3% в то время как при использовании модели (9) - (11) СКО оценки той же структуры составило 2,7%.
«
к й-к 3
! 8 « ю о к
ч о с
ей
и
К <и
а О
и о
ей
^
К о И ей
60 50 40
!30 20
10
0
>
8 16 24 32 40 48 56
Длина пути измерительной волны, мкм
Рис. 6. Оценка положения максимума огибающей интерференционного сигнала с использованием различных начальных условий для положения максимума
На рис. 6 приведены оценки максимума огибающей сигнала, полученные с использованием различной априорной информации о положении максимума для случаев,
0
когда в начале фильтрации априорное положение максимума предполагалось равным 0, 24, 32, 40, 48 и 53 мкм.
Из рис. 6 видно, что быстрота обнаружения максимума зависит от выбора начальных условий. Если заданное начальное значение положения максимума огибающей не превосходит истинного значения положения максимума, то коррекция оценки положения максимума начинается с первых точек сигнала. Если априорная оценка положения максимума превышает истинное значение, то коррекция оценки начинается с тех точек, где ожидается появление полезного сигнала. При этом на первых точках коррекция не происходит, что выражается в наличии плоских начальных участков кривых (рис. 6).
Заключение
В реальных интерферометрических системах движение опорного отражателя вносит фазовые отклонения, коррекция которых невозможна при линейной фильтрации. При использовании метода дискретной нелинейной фильтрации Калмана одновременная коррекция фазы и амплитуды сигнала позволяет повысить эффективность фильтрации.
Проведенные исследования показали, что с точки зрения наименьшей ошибки, предпочтительно использовать в качестве восстанавливаемых параметров фазу и амплитуду интерференционного сигнала. Преимущество использования в векторе параметров положения максимума огибающей сигнала состоит в возможности упреждающей оценки положения максимума и отсутствии в модели требований к априорной информации о количестве слоев исследуемой среды. Критерии, устанавливаемые для обнаружения границ слоев, вносят ограничения на разрешающую способность метода.
Использование метода дискретной нелинейной фильтрации Калмана дает более точные результаты при меньшей вычислительной сложности по сравнению с другими методами и обеспечивает возможность динамической обработки интерферометриче-ских сигналов.
Литература
1. Васильев В.Н., Гуров И.П. Компьютерная обработка сигналов в приложении к ин-терферометрическим системам. СПб: БХВ-Санкт-Петербург, 1998.
2. Kalman R.E. A new approach to linear filtering and prediction problems //Trans. ASME, J. Basic Eng. 1960. V. 82. P. 35-45.
3. Балакришнан А.В. Теория фильтрации Калмана. М.: Мир, 1988.
4. Захаров А.С. Нелинейный анализ стохастических параметров интерференционных систем. Диссертация кандидата технических наук. СПб: СПбГУ ИТМО. 2005.
5. Tuzlukov V.P. Signal Detection Theory. Boston: Birkhauser, 2001.
6. Коломийцов Ю.В. Интерферометры. Основы инженерной теории, применение. Л.: Машиностроение, 1976.
7. Гуров И.П. Оптическая когерентная томография: принципы, проблемы и перспективы. //В кн.: Проблемы когерентной и нелинейной оптики /Под ред. И. П. Гурова и С.А. Козлова, СПб: СПбГУ ИТМО, 2004. С. 6-30.
8. Справочник по прикладной статистике / под ред. Э. Ллойда, У. Лидермана. М.: Финансы и статистика. 1989. Т. 2. С. 421-470.
9. Alarousu E., Gurov I., Hast J., Myllyla R., Zakharov A. Optical coherence tomography of multilayer tissue based on the dynamical stochastic fringe processing. //Proc. SPIE. 2003. V. 5149. P. 13-20.
10. Gurov I., Sheynihovich D. Interferometric data analysis based on Markov non-linear filtering methodology //J. Opt. Soc. Am. A. 2000. V.17. P. 21-27.
