CC BY
38
УДК 535.41: 681.7.014.3
АНАЛИЗ ПАРАМЕТРОВ ВОЛНОВОГО ПАКЕТА НА ОСНОВЕ ДИНАМИЧЕСКОГО ОЦЕНИВАНИЯ ФАЗЫ СПЕКТРАЛЬНЫХ
СОСТАВЛЯЮЩИХ Е.А. Воробьева, И.П. Гуров
Рассмотрен метод динамического анализа параметров волнового пакета, состоящий в выделении и оценивании амплитуд и фаз отдельных монохроматических составляющих волнового пакета. Показано, что определение положения максимума интенсивности волнового пакета можно осуществить по критерию равенства нулю начальных фаз монохроматических составляющих, оцениваемых при помощи дискретного нелинейного фильтра Калмана. Ключевые слова: волновой пакет, спектральные составляющие, начальная фаза, фильтр Калмана.
Введение
При исследовании различных объектов оптическими методами используют разнообразные источники излучения. В ряде случаев излучение представляет собой совокупность монохроматических электромагнитных колебаний с различными частотами, образующих при распространении в пространстве волновой пакет. В зависимости от соотношения амплитуд и фаз составляющих волнового пакета возможно формирование распределений интенсивности различного вида. В частности, в области совпадения или малого отличия фаз формируется распределение интенсивности в форме импульса, положение максимума которого соответствует нулевому значению фазы всех монохроматических составляющих волнового пакета.
Свойства волнового пакета широко используются в оптической когерентной томографии (ОКТ) [1], поскольку идентификация положения максимума распределения интенсивности позволяет определить условие равенства оптической длины пути опорной и измерительной волн, тем самым выделить часть излучения, отраженного от отдельного слоя исследуемого частично прозрачного объекта, и в результате получить послойную томограмму внутренней микроструктуры объекта.
Определить положение максимума интенсивности можно различными методами, например, в результате детектирования огибающей распределения интенсивности [2] или нелинейной фильтрации распределения интенсивности [3]. Однако возможен и иной подход, предлагаемый в настоящей работе и состоящий в выделении и оценивании амплитуд и фаз отдельных монохроматических составляющих волнового пакета. При этом оказывается возможным определение положения максимума интенсивности по критерию равенства нулю начальных фаз монохроматических составляющих.
В работе рассматривается возможность фильтрации набора монохроматических составляющих в динамическом режиме с автоматическим определением положения максимума интенсивности импульса, представляющего волновой пакет.
Многочастотная модель интерферометрического сигнала для волнового пакета в системах ОКТ
Интерферометрический сигнал, характеризующий волновой пакет, можно представить в виде суммы гармонических составляющих с амплитудами a. и частотами f
s( f, г ) = ц/ (f, z) = a i cos (2nfz ) , (1)
где ц - коэффициент преобразования интенсивности /(f, z) в сигнал; z - независимая переменная. В случае огибающей гауссовой формы для интенсивности волнового пакета преобразование Фурье G(f) огибающей также имеет гауссову форму:
at = G(f) = exp(-n2 (f0 ±Sf )2 a2), (2)
где f0 - среднее значение частоты, f0 = 1/ (X), (X} - среднее значение длины волны излучения; 5f - отклонения частоты составляющих волнового пакета от среднего значения частоты.
При регистрации составляющие (1) на всех частотах суммируются, и для наблюдения доступен суммарный сигнал вида
s( z) = a,. cos(2nf z).
Поскольку в системах ОКТ разность хода опорной и измерительной волн меняется, в дискретном случае z = кAz , где Az - шаг дискретизации, регистрируемая последовательность отсчетов имеет вид
s(k) = ^a... cos(2nf к Az).
Для определения значения параметра a в уравнении (2) рассмотрим спектр волнового пакета. Ширина спектра сигнала на уровне 1/e
fc = 1/ (na). (3)
Обозначим количество периодов косинусоиды, укладывающихся в интервале ст , как Q. С учетом (3) можно записать Q = гост / Х .
С другой стороны, интервал 2ст есть ширина спектра на уровне 1/е и примерно соответствует длине когерентности излучения 1с = (Х)2 / ДХ, где ДХ - ширина спектра излучения. Нетрудно показать, что
ДХ / (Х) = Д/ / /0, откуда ДХ = (Х)Д/ / /0 = (Х)2 ДХ. С учетом (3) имеем Д/ = 2/с = 2/ гост и получаем
ДХ = 2( Х^2 / гост, откуда
ст=2(^2=^(Х.. (4)
го ДХ го ДХ
Поскольку ^Х) в (4) определяется источником излучения, т.е. задается в известной степени произвольно, можно отнормировать (4), приняв первый сомножитель равным единице. Тогда параметр ст* = (Х^ / ДХ будет соответствовать относительной ширине спектра, т.е. числу, которое показывает, во
сколько раз средняя длина волны больше, чем ширина спектра. Для типичных источников излучения, используемых с системах ОКТ, ст * находится в интервале 10-50.
Восстановление параметров волнового пакета методом дискретной нелинейной фильтрации Калмана
Параметры интерферометрического сигнала с учетом отражения оптического излучения в образце можно описать в виде вектора
е = (Ф, (к), б,, (к), / (к) )т (5)
где Ф,(к) и е,(к) - полная и начальная фазы, соответствующие частоте /(к), к=0, 1, ..., К - номер дискретного отчета в области независимой переменной; К - общее количество отсчетов; , = 0, ±1,..., ±(Ь -1) / 2, (2Ь+1) - количество частот в волновом пакете.
Оценку вектора параметров (5) в динамическом режиме можно получить с использованием фильтра Калмана [4]. При этом проводится обработка последовательности отсчетов, определяемой заданной моделью, которая выражается уравнением наблюдения, включающим влияние шума наблюдения. Эволюция вектора параметров характеризуется уравнением системы с матрицей перехода от предыдущего отсчета (состояния системы) к последующему отсчету.
В интерферометрической системе уравнение наблюдения может быть записано как
s(k) = Xа соз(Ф,. (к,е,,/(к))) + п(к), где п(к) - последовательность отсчетов шума наблюдения.
Матрица перехода А(к) для вектора параметров в уравнении системы имеет размерность (6Ь + 3) х (6Ь + 3) и определяется как
(1 0 2%Дг ... 0 ^
А(к) =
\
0 10 0
0 0 1 0
... ... 0
0 0 0 0 1
Алгоритм дискретной нелинейной фильтрации Калмана, подробно описанный, например, в работе [5], с использованием рассмотренной выше многочастотной модели позволяет восстановить параметры, т.е. амплитуду и фазу составляющих интерферометрического сигнала. Однако при использовании модели с достаточно большим количеством частот Ь быстродействие алгоритма обработки заметно снижается. Для получения результатов, приведенных в работе, была использована модель сигнала с 30 частотами.
На рис. 1, а, представлен исходный смоделированный сигнал. Огибающая сигнала имеет два равных по амплитуде максимума в точках, которым соответствуют дискретные отсчеты сигнала 2980 и 6770. Несущая частота интерферометрического сигнала /0 = 0,015 , ст* = 15 .
Дискретный набор частот определяется как / = /0 ± ,Д/ , где Д/ = 3 • 10-5. В данном случае максимальное отклонение по частоте от несущей составляет Д/ • (Ь /2) = 0,0005. Таким образом, набор частот попадает в диапазон (/-,/+), где граничные частоты равны 0,0145 и 0,0155 соответственно. На
рис. 1, б, представлен пример спектра огибающей сигнала. Пунктирными линиями обозначен диапазон, соответствующий интервалу , /с ). На рис. 1, б, также отмечены первые 9 частот, которые использованы в модели интерферометрического сигнала в алгоритме фильтрации Калмана.
л ч А
Й и 0
2 « I н
« о «
о
и
<и
к
2000 4000 6000 8000 Номер дискретного отсчета
10000
0,013
0,0145 0,015 0,0155 Частота б
0,017
Рис. 1. Смоделированный сигнал (а) и спектр огибающей сигнала (б)
Покажем, что алгоритм дискретной нелинейной фильтрации Калмана позволяет восстановить параметры волнового пакета. На рис. 2, а, представлены восстановленные значения частот для трех сигналов, которым соответствует центральная частота (сплошная линия) и две боковые частоты /+ = / + А/ и /_ = /0 - А/ (пунктирные линии). Как видно на рис. 2, а, значение частот мало корректируется в процессе фильтрации, поскольку было задано корректное начальное значение.
„ 7
й н о н о Й
0,01504
0,01502
0,015
0,01498
« й й о К л
н
й ^
й К
2000 4000 6000 8000 Номер дискретного отсчета
10000
2000 4000 6000 8000 Номер дискретного отсчета
б
10000
Рис. 2. Восстановленные значения частот (а) и начальных фаз для трех составляющих
волнового пакета (б)
Наибольший интерес представляют восстановленные значения начальных фаз составляющих сигнала, поскольку на основе этих значений можно восстановить положения максимума огибающей сигнала, соответствующие нулевым фазам всех составляющих. На рис. 2, б, представлены восстановленные значения начальных фаз для сигналов, которым соответствуют частоты /0, /+ и /_. Жирной линией показана начальная фаза сигнала с центральной частотой /0, тонкими линями - начальные фазы сигналов с боковыми частотами /+ и /_. Как видно на рис. 2, б, значение начальной фазы остается постоянным на интервале после появления очередного максимума огибающей и до появления следующего максимума. Важно отметить, что на представленном графике значение начальной фазы для сигнала, соответствующего центральной частоте, превышает значение 2л. Это связано с особенностями реализации алгоритма восстановления параметров сигнала. В данном случае начальная фаза сигнала есть остаток от деления на 2л восстановленного значения.
Для вычисления положений максимума огибающей рассмотрим три частоты /0, /+ и /_. Значение сигнала можно представить в виде
) = а0 соБ(е0 + 2л/0к) + а1 соб^ + 2л/+ к) + а2 С0Б(е2 + 2л/_к) где е0, е1 и е2 - начальные фазы сигналов с частотами /0, /+ и /_ соответственно. Значения а0, а1 и а2 вычисляются по формуле (2).
При формировании волнового пакета положению максимума огибающей соответствует дискретный отсчет сигнала, в котором начальные фазы составляющих сигнала, формирующих волновой пакет, совпадают и равны нулю. Можно записать условия максимума огибающей волнового пакета
е0 + 2л/0к0 = е1 + 2л/+ к0 + 2лп1 = е2 + 2л/_к0 + 2лп2, (6)
0
а
0
0
а
где п1 и п2 - целые числа; к0 - положение максимума огибающей волнового пакета. В выражении (6) 2пп1 и 2пп2 - это фазы, которые «набегут» за счет разностей между /0 и /+ , /0 и /_ соответственно. Поскольку значение Д/ задано, то можно записать
2^ = 2пп2 = 2%Д/к0,
откуда
ко = ^ = П-. (7)
0 Д/ Д/
На рис. 3, а, представлены восстановленные значения положения максимумов огибающей волнового пакета. Пунктирной линией обозначены истинные значения положения максимумов. Видно, что положение первого максимума алгоритм восстанавливает достаточно точно. Однако положение второго максимума определено не совсем точно, поскольку алгоритм учитывает влияние предыдущего максимума огибающей. Для устранения этого недостатка можно использовать переключаемую модель сигнала [3].
1
0,8
т о
5!
¡Г
к «
о <и
9 1°,6
Тл Я
«0,4
и о
0,2
2000 4000 6000 8000 10000 Номер дискретного отсчета
2000 4000 6000 8000 Номер дискретного отсчета
10000
а б
Рис. 3. Восстановленные значения положений максимумов огибающей сигнала, полученные с использованием уравнения (7) (а) и критерия равенства полной фазы сигналов 2п (б)
Другой способ оценки положения максимумов огибающей основан на анализе полных фаз составляющих сигнала, формирующих волновой пакет. Положению максимумов огибающей соответствуют дискретные отсчеты, для которых полная фаза сигналов кратна 2п. На рис. 3, б, представлены оценки максимума огибающей, полученные с использованием описанного критерия. Для трех сигналов, которым соответствуют частоты /0, /+и/_, были вычислены полные фазы. На графике дискретные отсчеты, для которых полные фазы трех сигналов кратны 2п, отмечены дельта-функциями. Алгоритм корректно оценил положение существующих максимумов огибающей волнового пакета, однако при этом возникают ложные максимумы. Ложный максимум в начале дискретной последовательности обусловлен начальной подстройкой дискретного нелинейного алгоритма фильтрации Калмана. Использование только трех частот для определения положения максимумов является причиной ложных результатов в остальных отсчетах дискретной последовательности. При увеличении количества частот надежность результатов заметно возрастает, однако это происходит за счет повышения вычислительной сложности.
1
к 0,6
о
^ 0,4
<Ц
О 0,2
А
и
0 2000 4000 6000 8000 10000 0 0,005 0,01 0,015 0,02 0,025 0,03 Номер дискретного отсчета Частота
а б
Рис. 4. Двухчастотный сигнал (а) и его спектр (б)
Рассмотрим волновой пакет с изменяющейся несущей частотой (рис. 4). Положению максимума огибающей волнового пакета соответствует дискретный отсчет 4729. Несущая частота волнового пакета до отсчета 4729 равна 0,015, после - 0,011. Рис. 4, б, иллюстрирует спектр сигнала.
0
0
Алгоритм дискретной нелинейной фильтрации Калмана с использованием многочастотной модели сигнала позволяет восстановить центральную частоту волнового пакета. На рис. 5 представлена оценка восстановленной центральной частоты волнового пакета. На интервале отсчетов от 0 до 4729 оценка значения частоты постепенно приближается к заданному значению 0,015. После отсчета 4729 алгоритм фиксирует изменение частоты. Восстановленное значение частоты постепенно корректируется до истинного значения 0,011.
0,016
Я 0,015
ч
£ 0,014
" 0,013
га '
Н
н 0,012
ей
0,01
0 2000 4000 6000 8000 10000 Номер дискретного отсчета Рис. 5. Восстановленные значения центральной частоты
Таким образом, алгоритм дискретной нелинейной фильтрации Калмана с использованием многочастотной модели позволяет восстановить параметры интерферометрических сигналов, формирующих волновой пакет, в том числе при локальном изменении среднего значения частоты сигнала.
Заключение
Проведенное рассмотрение особенностей формирования и анализа интерферометрических сигналов в системах оптической когерентной томографии на основе многочастотной модели сигнала показывает возможность получения динамических оценок параметров волнового пакета, в том числе при изменениях значения средней частоты сигнала волнового пакета. Алгоритм нелинейной фильтрации Калмана с использованием предложенной многочастотной модели интерферометрических сигналов позволил получить оценки начальной фазы частотных составляющих волнового пакета и определить положение максимумов огибающей по критерию равенства нулю начальных фаз составляющих и критерию равенства полной фазы целому числу фазовых циклов. При реализации рассмотренного алгоритма необходимо учитывать возрастание вычислительной сложности при увеличении числа учитываемых в модели сигнала частотных составляющих. Обработка большего числа составляющих обеспечивает повышение точности определения положения максимума огибающей волнового пакета.
Работа выполнена при финансовой поддержке Министерства образования и науки Российской Федерации.
Литература
1. Гуров И.П. Оптическая когерентная томография: принципы проблемы и перспективы // Проблемы когерентной и нелинейной оптики / Под ред. И.П. Гурова и С.А. Козлова. - СПб: СПбГУ ИТМО, 2004. - С. 6-30.
2. Воробьева Е.А., Гуров И.П., Петерсон М.В. Исследование метода обработки сигналов в системах оптической когерентной томографии с повышенным быстродействием // Изв. вузов. Приборостроение. - 2010. - № 3. - С. 65-73.
3. Gurov I., Volynsky M., Zakharov A. Evaluation of multilayer tissue in optical coherence tomography by the extended Kalman filtering method // Proc. SPIE. - 2007. - V. 6734. - P. 67341.
4. Kalman R.E. A new approach to linear filtering and prediction problems // Transactions of the ASME. -Journal of Basic Engineering, Series D. - 1960. - V. 82. - P. 34-45.
5. Гуров И.П., Захаров А.С. Анализ характеристик интерференционных полос методом нелинейной фильтрации Калмана // Оптика и спектроскопия. - 2004. - Т. 96. - № 2. - С. 210-216.
Воробьева Елена Александровна - ООО «Моторола-Мобилити», инженер-программист,
Гуров Игорь Петрович - Санкт-Петербургский национальный исследовательский университет
информационных технологий, механики и оптики, доктор технических наук, профессор, зав. кафедрой, [email protected]
