CC BY
48
Васильев Владимир Николаевич - Россия, Санкт-Петербург, Санкт-Петербургский национальный исследовательский университет информационных технологий, механики и оптики, доктор технических наук, профессор, ректор, [email protected] Гуров Игорь Петрович - Россия, Санкт-Петербург, Санкт-Петербургский национальный исследова-
тельский университет информационных технологий, механики и оптики, доктор технических наук, профессор, зав. кафедрой, [email protected] Потапов Алексей Сергеевич - Россия, Санкт-Петербург, Санкт-Петербургский национальный исследова-
тельский университет информационных технологий, механики и оптики, профессор, доктор технических наук, [email protected]
УДК 581.787: [519.6+517.443]
МЕТОД ДИНАМИЧЕСКОЙ ОБРАБОТКИ ДАННЫХ В СПЕКТРАЛЬНОЙ ОПТИЧЕСКОЙ КОГЕРЕНТНОЙ ТОМОГРАФИИ С КОМПЕНСАЦИЕЙ ВЛИЯНИЯ ДИСПЕРСИИ М.А. Волынский, И.П. Гуров
При формировании сигналов в спектральной оптической когерентной томографии при наличии дисперсии в среде спектральные интерференционные полосы испытывают частотную модуляцию из-за зависимости частоты полос от длины волны, что приводит к уширению спектра регистрируемого сигнала и ухудшению разрешающей способности спектрального интерферометра. В работе предложен метод динамической обработки данных спектральной интерференции на основе алгоритма дискретной линейной фильтрации Калмана с компенсацией влияния дисперсии в среде на разрешающую способность при исследовании оптически неоднородных частично прозрачных объектов. Алгоритм состоит в идентификации параметров (амплитуды и начальной фазы) гармонических составляющих интерференционного сигнала с фиксированным набором частот с помощью дискретного линейного фильтра Калмана. Использование информации о начальной фазе позволяет скомпенсировать влияние дисперсии и устранить нежелательные артефакты, что дает возможность повысить разрешающую способность спектральной оптической когерентной томографии. Представлены результаты обработки одномерных и двумерных сигналов в спектральной оптической когерентной томографии на примере исследования случайно-неоднородных рассеивающих сред в биомедицине. Ключевые слова: спектральная интерферометрия, оптическая когерентная томография, компенсация дисперсии в среде, преобразование Фурье, фильтр Калмана.
Введение
Методы бесконтактного контроля внутренней микроструктуры объектов необходимы для многих областей науки и техники. Одним из таких методов является оптическая когерентная томография (ОКТ), получившая распространение в целях неинвазивной диагностики объектов в биомедицине [1, 2]. Как известно, принцип ОКТ может быть реализован при использовании корреляционного или спектрального интерферометра [1-4]. В корреляционном интерферометре осуществляют перемещение оптической системы относительно исследуемого объекта. При этом интерференционные полосы малой когерентности формируются в пределах длины когерентности излучения при интерференции части измерительной волны, отраженной от поверхности непрозрачного объекта или от слоя частично прозрачного неоднородного объекта, находящихся от светоделителя на расстоянии, равном оптической длине пути опорной волны. В спектральном интерферометре оптическая длина пути опорной волны не равна оптической длине пути измерительной волны для всего диапазона высот рельефа непрозрачного объекта или глубины частично прозрачного объекта. На выходе интерферометра размещен спектральный прибор, позволяющий определить составляющую отраженной измерительной волны для различных длин волн. При этом, используя преобразование Фурье-спектра, зарегистрированного фотодетектором, можно определить расстояние и степень отражения от каждого слоя [4].
При наличии дисперсии в среде спектральные интерференционные полосы испытывают частотную модуляцию из-за зависимости частоты полос от длины волны, что приводит к уширению спектра регистрируемого сигнала и ухудшению разрешающей способности спектрального интерферометра.
Традиционный метод обработки данных в спектральной ОКТ основан на преобразовании Фурье (см., например, [4]). При этом для компенсации влияния дисперсии требуется точное определение амплитуд, частот и фаз информативных составляющих спектра сигнала, которое оказывается возможным только при регистрации полной реализации сигнала во всем диапазоне длин волн перед обработкой данных [5, 6], что снижает быстродействие системы. Повышение быстродействия возможно при получении динамических оценок параметров сигнала при рекуррентной обработке данных в системах спектральной ОКТ с перестраиваемой длиной волны источника излучения на основе фильтрации Калмана [7]. Алгоритм рекуррентной обработки данных в системах спектральной ОКТ предложен в [8, 9]. В настоящей работе рассматривается модифицированный метод динамической обработки сигналов спектральной интерференции на основе фильтра Калмана с возможностью повышения разрешающей способности за счет компенсации влияния дисперсии в среде.
Основы метода спектральной оптической когерентной томографии
В спектральной интерферометрии регистрируют сигналы, пропорциональные значениям интен-
сивности [4], т.е. I (а) = G(ct)
aR exp(j2аг) + Ja(z)exp{j2a[r + n(z)z]}dz
(1)
где О(ст) - спектр источника излучении; ак - амплитуда опорной волны; ] = >/-1; ст = 2 л / X - волновое число, определяемое длиной волны X; 2г - оптическая длина пути опорной волны; а(г) - амплитуда предметной волны, отраженной на глубине г в исследуемом образце в диапазоне глубин Ь; п(г) - изменение показателя преломления по глубине среды. Полезная информация содержится в значениях а(г), характеризующих степень отражения предметной волны по глубине исследуемой среды; их необходимо определить в результате обработки полученных значений I (ст).
В результате нормировки выражения (1) относительно спектра источника излучения интерферо-метрический сигнал в спектральной ОКТ определяется в области волновых чисел формулой (см., например, [5, 6])
S (а) =
1 + J a( z) exp[ j 2an0 z + (a)]dz
(2)
где n0 = const - эффективный показатель преломления в среде; <d - начальная фаза каждой из гармонических составляющих сигнала, обусловленная дисперсией в среде (для среды без дисперсии <d = 0).
Обозначим А (а) интеграл под модулем в выражении (2). Стандартный метод вычисления искомой амплитуды отраженной волны a( z) состоит в применении обратного преобразования Фурье к выделенной информативной составляющей A (а) сигнала S(а):
да
a(z) =JА(а)ехр(—/2ап0z)dа . (3)
0
Заметим, что преобразование (3) является непараметрическим, поскольку позволяет получать результаты независимо от вида (параметров) преобразуемой функции А (а).
Величину А (а) можно рассматривать как набор отдельных гармонических составляющих, определяемых параметрами - частотой, амплитудой и начальной фазой. При этом задача сводится к параметрической идентификации и может быть решена с помощью линейной фильтрации Калмана [8, 9]. Алгоритм дискретной фильтрации Калмана основан на рекуррентной процедуре предсказания значения сигнала на последующий шаг на основе информации, имеющейся на предыдущем шаге, с учетом известной параметрической модели сигнала и использовании ошибки предсказания для уточнения значений параметров. Применительно к спектральной интерферометрии это позволяет выбирать требуемое количество длин волн, обеспечивающее необходимую разрешающую способность и быстродействие.
Алгоритм обработки данных
Для описания алгоритма обработки данных перейдем к их дискретному представлению и введем номер дискретного отсчета к = 0.. ,Л"-1 в области волновых чисел, причем ак = к Да, где Да - шаг изменения волнового числа.
Пусть имеется последовательность отсчетов А (к) = А(ак) сигнала спектральной интерференции и решается задача рекуррентной идентификации гармонической составляющей с некоторой фиксированной частотой f В случае модели сигнала вида
H(к) = u cos(2fДа - <), (4)
где < - начальная фаза на частоте f алгоритм рекуррентной идентификации наличия составляющей с частотой f сводится к фильтрации параметров - амплитуды u и начальной фазы < гармонического сигнала. Модель (4) может быть представлена в виде
H(к) = uc cos(2fДа) + us sin(2fДа), (5)
а искомые амплитуда u и начальная фаза < определяются как
u = Vu2 + u2 , (6)
< = arctgu / uc). (7)
2
2
С учетом сказанного, поиск амплитуды и начальной фазы гармонической составляющей сигнала для каждой частоты / сводится к оцениванию значений ис и их с уточнением оценки при поступлении
каждого отсчета сигнала к. Введем вектор параметров и = (ис ил. . Тогда выражение для априорной оценки этого вектора параметров можно записать на к-ом шаге обработки в виде
и(к) = аи(к -1), (8)
где вектор коэффициентов а задает линейный закон изменения амплитуд (при малых приращениях Дст можно считать, что а = 1). Апостериорная коррекция осуществляется на основе невязки предсказания и наблюдения:
и (к) = и(к) + Р(к)[ ИоЫ (к) - Н (к)]. (9)
Для вычислений по формуле (9) требуется определить коэффициент усиления Р(к) = ЯСТ (к)[С(к)ЯСТ (к) + ]-1, где Я - ковариационная матрица ошибки априорной оценки параметров; - дисперсия шума наблюдения; С(к) - матрица перехода, элементы которой являются производными модели (5) по параметрам, С(к) = Н'и (к), Н(к) - прогноз наблюдения, получаемый подстановкой предсказанных параметров в модель (5).
Если шаг дискретизации Дст отнормировать и принять равным единице, алгоритм фильтрации (8)-(9) сведется к следующему рекуррентному выражению:
и(к) = и (к -1) + ЯСТ (к)[С(к)ЯСт (к) + К„\1[НоЫ (к) - Н (к)], (10)
где
= Г ссад У, (11)
^ ¡¡¡ш(27г/к))
Н(к) = ис (к -1) со8(2л/к) + и (к -1) 8т(2яД). (12)
Следует отметить, что на каждом шаге фильтрации возможен пересчет ковариационной матрицы Я согласно выражению [5]
И(к) = [I - Р(к )С(к )]И (к -1), (13)
где I - единичная матрица. Использование (13) позволяет дополнительно повысить скорость сходимости алгоритма, однако увеличивает вычислительную сложность.
Из выражения (3) видно, что а(г) может быть представлена в виде набора 5-функций,
М-1
а(г) = £ а,5(г - г,.), (14)
,=0
соответствующих гармоническим составляющим спектра и характеризующих положение слоев исследуемой среды на глубине г, что эквивалентно представлению в частотной области (области длин волн) в виде М гармонических составляющих. С учетом выражения (14) возможна одновременная идентификация этих составляющих при использовании параллельной обработки (см. [9]).
В качестве примера работы алгоритма (10)-(13) на рис. 1, а, представлен исходный сигнал со следующими параметрами: нормированная частота f = 0,01, амплитуда 10 усл. ед., начальная фаза 2 рад, аддитивный белый гауссовский шум с нулевым средним и дисперсией 1 усл. ед. На рис. 1, б, в, представлены результаты оценивания параметров гармонической составляющей при использовании рассмотренного выше рекуррентного алгоритма. Видно, что после примерно 100 отсчетов (один период сигнала) оценки выходят на квазистационарные значения, соответствующие истинным. Незначительные флуктуации оценок обусловлены наличием шума наблюдения с ненулевой дисперсией, что имеет место в реальных сигналах. Из гистограмм ошибок оценок амплитуды и фазы (рис. 1, д, е) видно, что они имеют од-номодовый характер с нулевым математическим ожиданием и близки к гауссовской форме.
Поскольку для каждого сигнала оцениваемые параметры - константы, а их оценки принимают квазистационарные значения в пределах одного периода сигнала, при реализации алгоритма (10)-(13) используется критерий останова, предложенный в работе [8]. Указанный критерий основан на определении разности оценок амплитуды сигнала с заданной частотой на текущем к и предыдущем к - 1 отсчетах и сравнении с некоторым пороговым значением е. Тогда критерий останова выражается в форме
и (к) - и (к -1) <8 . (15)
Выбор порогового значения е может осуществляться различными способами и зависит от представления исходного сигнала. Если значения сигнала принадлежат множеству действительных чисел, то в качестве е целесообразно выбирать среднее квадратичное отклонение шума наблюдения. В системах ОКТ исходные данные, как правило, представляют собой изображения, представленные 256 или 4096 градациями серого для 8- и 12-битного представления соответственно. В этом случае в качестве е целе-
сообразно выбирать минимально возможное изменение градаций серого. Особенности выбора критерия (15) рассмотрены в работах [8, 9].
15
5
ан
г
и цд
с ш -5
е и d
н -15
е ч
а
н
го
£ ■ 25 £ ^ 20 Н 5 15 ^ 10
0 50 100 150 200 Номер дискретного отсчета
а
!
и
о Я
О
«
о и
Л
4
Л £
5
и»
И Л
<D ет Я 5
о
0х
л" н о о и
о а
(D
m
50 100 150 200 Номер дискретного отсчета в
к н о о
X 0
(U &
о с
(U К X <D
£ X со
1,
5
50 100 150 200 Номер дискретного отсчета б
5
0,5 0
-0,5 -1
0 50 100 150 200 Номер дискретного отсчета г
-2,8 -1,4 0 1,4 2
15
гтт
-2,8 -1,4 0 1,4 2,8 Погрешность оценки амплитуды, х0,01 усл. ед.
д
-18,7 -9,33 0 9,33 Погрешность оценки фазы, х0,0001 рад е
18,7
Рис. 1. Исходный сигнал с частотой 1 = 0,01 (а), оценки амплитуды (б) и начальной фазы (в), погрешность предсказания значения сигнала (г) и гистограммы ошибок предсказания амплитуды (д) и начальной фазы (е) для 104 реализаций сигнала
Найденные в (6) и (7) значения амплитуд и начальных фаз могут быть использованы для компенсации влияния дисперсии в среде [5, 6]. Следует отметить, что начальные фазы для гармонических составляющих с амплитудой, близкой к нулю, не имеют физического смысла, поэтому целесообразно ввести пороговое значение оценки амплитуды (например, незначительно превышающее среднее квадратичное отклонение шума наблюдения) и считать начальные фазы гармонических составляющих с амплитудой, меньшей этого порогового значения, равными нулю. Значение начальной фазы для гармонических составляющих с ненулевой амплитудой может быть использовано для компенсации влияния дисперсии в среде с помощью алгоритма, предложенного в работе [5].
Как отмечается в [5], из-за влияния дисперсии происходит уширение полос в спектре интерферо-метрического сигнала, причем истинному положению полосы соответствует нулевая производная начальной фазы. С учетом этого можно ввести коэффициент пропорциональности на основе найденных значений начальной фазы на каждой частоте и корректировать амплитуду А-скана в соответствии с выражением
1 ё ф( г)
а( z) = 11 -
а( z);
(16)
Дст dz
где а( z) - скорректированная с учетом влияния дисперсии амплитуда.
В качестве примера на рис. 2, а, показан сигнал спектральной интерференции, полученный при исследовании двухслойного фантома из поливинилхлорида [10] с помощью спектрального оптического когерентного томографа Hyperion (ThorLabs, США) с источником излучения с центральной длиной волны 930 нм [11]. На рис. 2, б, показан результат обработки с компенсацией дисперсии и без компенсации. Видно, что наличие дисперсии в среде уширяет линии в спектре сигнала, а поправка (16) на основе учета значений начальной фазы позволяет скомпенсировать этот эффект, повышая, тем самым, разрешающую способность по глубине. Толщина верхнего слоя в фантоме при его изготовлении составляла около 100 мкм, что иллюстрируется на рис. 2, б.
0
а
0
и
и
S о о S И о
S
И
СО
1
0,6 0,4 0,2 0
850
890 930 Длина волны, а
970
15
(D
S
я s -е -е
m
â
Ч
<D
M н о
(U
s и
<4
0,5
о
1 1. ) \
i 1 № il 11
100
Глубина, мкм б
200
300
Рис. 2. Сигнал спектральной интерференции (а) и полученный из него А-скан (б) с компенсацией дисперсии (сплошная линия) и без компенсации (пунктирная линия)
Обработка биомедицинских данных
На рис. 3 показан пример В-скана кожи на подушечке указательного пальца человека, полученный с помощью описанного в настоящей работе метода без учета начальных фаз (рис. 3, а, б) и с компенсацией дисперсии (рис. 3, в, г). Видно, что без учета начальных фаз на реконструируемых В-сканах присутствуют артефакты в виде вертикальных полос, вызванные уширением полосы частот, соответствующей каждой глубине объекта, вследствие влияния дисперсии в среде. На рис. 3, в, г, указанные артефакты отсутствуют.
1
0
нм
в г
Рис. 3. B-скан кожи на подушечке указательного пальца человека, восстановленный без учета влияния дисперсии (а)-(б) и с компенсацией дисперсии (в)-(г)
Заключение
В системах спектральной оптической когерентной томографии с перестраиваемой длиной волны актуальна возможность компенсации влияния дисперсии среды на получаемые изображения. Для повышения быстродействия и достоверности результатов в системах спектральной оптической когерентной томографии предложен алгоритм обработки данных на основе линейного фильтра Калмана, позволяющий проводить идентификацию частот из заданного набора в сигнале и определять их параметры (амплитуду и начальную фазу). Использование информации о начальной фазе позволяет скомпенсировать влияние дисперсии и устранить нежелательные артефакты, возникающие при реконструкции В-сканов, что позволяет повысить разрешающую способность спектральной оптической когерентной томографии.
Апробация алгоритма при анализе биомедицинских данных демонстрирует адекватность предложенного метода и его применимость в системах спектральной оптической когерентной томографии с использованием современных источников излучения и повышенным быстродействием.
Работа выполнена при финансовой поддержке Министерства образования и науки Российской Федерации.
Литература
1. Tomlins P.H., Wang R.K. Theory, developments and applications of optical coherence tomography // J.
Phys. D: Appl. Phys. - 2005. - V. 38. - P. 2519-2535.
2. Drexler W., Fujimoto J.G. eds. Optical Coherence Tomography. Technology and Applications. - BerlinHeidelberg: Springer-Verlag, 2008. - 1346 p.
3. Волынский М.А., Воробьева Е.А., Гуров И.П., Маргарянц Н.Б. Бесконтактный контроль микрообъектов методами интерферометрии малой когерентности и оптической когерентной томографии // Изв. вузов. Приборостроение. - 2011. - Т. 54. - № 2. - С. 75-82.
4. Häusler G., Lindner M.V. «Coherence radar» and «Spectral radar» - new tools for dermatological diagnostics // J. Biomed. Opt. - 1998. - V. 3. - P. 21-31.
5. Hillmann D., Bonin T., Lührs C., Franke G., Hagen-Eggert M., Koch P., Hüttmann G. Common approach for compensation of axial motion artifacts in swept-source OCT and dispersion in Fourier-domain OCT // Opt. Expr. - 2012. - V. 20. - № 6. - P. 6761-6776.
6. Lippok N., Coen S., Nielsen P., Vanholsbeeck F. Dispersion compensation in Fourier domain optical coherence tomography using the fractional Fourier transform // Opt. Expr. - 2012. - V. 20. - № 21. - P. 23398123413.
7. Kalman R.E. A new approach to linear filtering and prediction problems // Trans. ASME, J. Basic Eng. -1960. - V. 82. - P. 35-45.
8. Волынский М.А., Гуров И.П. Рекуррентная обработка данных в спектральной оптической когерентной томографии на основе фильтрации Калмана // Научно-технический вестник информационных технологий, механики и оптики. - 2013. - № 2 (84). - С. 40-44.
9. Gurov I., Volynsky M. Recurrence signal processing in Fourier-domain optical coherence tomography based on linear Kalman filtering // Proc. SPIE. - 2013. - V. 8792. - P. 879203-1-879203-6.
10. Быков А.В., Волков М.В., Волынский М.А., Гуров И.П., Киннунен М., Маргарянц Н.Б., Попов А.П. Изготовление тканеимитирующих фантомов и капилляров и их исследование методом оптической когерентной томографии // Научно-технический вестник информационных технологий, механики и оптики. - 2013. - № 2 (84). - С. 54-59.
11. Спецификация прибора Hyperion на сайте компании ThorLabs [Электронный ресурс]. - Режим доступа: http://www.thorlabs.com/newgrouppage9.cfm?objectgroup_id=5701, свободный. Яз. англ. (дата обращения 11.09.2013).
Волынский Максим Александрович
Гуров Игорь Петрович
Россия, Санкт-Петербург, Санкт-Петербургский национальный исследовательский университет информационных технологий, механики и оптики, кандидат технических наук, доцент,
Россия, Санкт-Петербург, Санкт-Петербургский национальный исследовательский университет информационных технологий, механики и оптики, доктор технических наук, профессор, зав. кафедрой, [email protected]
