Например, Бобцов

МЕТОД ДИНАМИЧЕСКОЙ ОБРАБОТКИ ДАННЫХ В СПЕКТРАЛЬНОЙ ОПТИЧЕСКОЙ КОГЕРЕНТНОЙ ТОМОГРАФИИ С КОМПЕНСАЦИЕЙ ВЛИЯНИЯ ДИСПЕРСИИ

МЕТОД ДИНАМИЧЕСКОЙ ОБРАБОТКИ ДАННЫХ …

Васильев Владимир Николаевич – Россия, Санкт-Петербург, Санкт-Петербургский национальный исследова-

тельский университет информационных технологий, механики и оптики,

доктор технических наук, профессор, ректор, vasilev@mail.ifmo.ru

Гуров Игорь Петрович

– Россия, Санкт-Петербург, Санкт-Петербургский национальный исследова-

тельский университет информационных технологий, механики и оптики,

доктор технических наук, профессор, зав. кафедрой, gurov@mail.ifmo.ru

Потапов Алексей Сергеевич

– Россия, Санкт-Петербург, Санкт-Петербургский национальный исследова-

тельский университет информационных технологий, механики и оптики,

профессор, доктор технических наук, pas.aicv@gmail.com

УДК 581.787:[519.6+517.443]
МЕТОД ДИНАМИЧЕСКОЙ ОБРАБОТКИ ДАННЫХ В СПЕКТРАЛЬНОЙ ОПТИЧЕСКОЙ КОГЕРЕНТНОЙ ТОМОГРАФИИ
С КОМПЕНСАЦИЕЙ ВЛИЯНИЯ ДИСПЕРСИИ
М.А. Волынский, И.П. Гуров
При формировании сигналов в спектральной оптической когерентной томографии при наличии дисперсии в среде спектральные интерференционные полосы испытывают частотную модуляцию из-за зависимости частоты полос от длины волны, что приводит к уширению спектра регистрируемого сигнала и ухудшению разрешающей способности спектрального интерферометра. В работе предложен метод динамической обработки данных спектральной интерференции на основе алгоритма дискретной линейной фильтрации Калмана с компенсацией влияния дисперсии в среде на разрешающую способность при исследовании оптически неоднородных частично прозрачных объектов. Алгоритм состоит в идентификации параметров (амплитуды и начальной фазы) гармонических составляющих интерференционного сигнала с фиксированным набором частот с помощью дискретного линейного фильтра Калмана. Использование информации о начальной фазе позволяет скомпенсировать влияние дисперсии и устранить нежелательные артефакты, что дает возможность повысить разрешающую способность спектральной оптической когерентной томографии. Представлены результаты обработки одномерных и двумерных сигналов в спектральной оптической когерентной томографии на примере исследования случайно-неоднородных рассеивающих сред в биомедицине. Ключевые слова: спектральная интерферометрия, оптическая когерентная томография, компенсация дисперсии в среде, преобразование Фурье, фильтр Калмана.
Введение
Методы бесконтактного контроля внутренней микроструктуры объектов необходимы для многих областей науки и техники. Одним из таких методов является оптическая когерентная томография (ОКТ), получившая распространение в целях неинвазивной диагностики объектов в биомедицине [1, 2]. Как известно, принцип ОКТ может быть реализован при использовании корреляционного или спектрального интерферометра [1–4]. В корреляционном интерферометре осуществляют перемещение оптической системы относительно исследуемого объекта. При этом интерференционные полосы малой когерентности формируются в пределах длины когерентности излучения при интерференции части измерительной волны, отраженной от поверхности непрозрачного объекта или от слоя частично прозрачного неоднородного объекта, находящихся от светоделителя на расстоянии, равном оптической длине пути опорной волны. В спектральном интерферометре оптическая длина пути опорной волны не равна оптической длине пути измерительной волны для всего диапазона высот рельефа непрозрачного объекта или глубины частично прозрачного объекта. На выходе интерферометра размещен спектральный прибор, позволяющий определить составляющую отраженной измерительной волны для различных длин волн. При этом, используя преобразование Фурье-спектра, зарегистрированного фотодетектором, можно определить расстояние и степень отражения от каждого слоя [4].
При наличии дисперсии в среде спектральные интерференционные полосы испытывают частотную модуляцию из-за зависимости частоты полос от длины волны, что приводит к уширению спектра регистрируемого сигнала и ухудшению разрешающей способности спектрального интерферометра.
Традиционный метод обработки данных в спектральной ОКТ основан на преобразовании Фурье (см., например, [4]). При этом для компенсации влияния дисперсии требуется точное определение амплитуд, частот и фаз информативных составляющих спектра сигнала, которое оказывается возможным только при регистрации полной реализации сигнала во всем диапазоне длин волн перед обработкой данных [5, 6], что снижает быстродействие системы. Повышение быстродействия возможно при получении динамических оценок параметров сигнала при рекуррентной обработке данных в системах спектральной ОКТ с перестраиваемой длиной волны источника излучения на основе фильтрации Калмана [7]. Алгоритм рекуррентной обработки данных в системах спектральной ОКТ предложен в [8, 9]. В настоящей работе рассматривается модифицированный метод динамической обработки сигналов спектральной интерференции на основе фильтра Калмана с возможностью повышения разрешающей способности за счет компенсации влияния дисперсии в среде.

50 Научно-технический вестник информационных технологий, механики и оптики,
2013, № 6 (88)

М.А. Волынский, И.П. Гуров

Основы метода спектральной оптической когерентной томографии

В спектральной интерферометрии регистрируют сигналы, пропорциональные значениям интен-

сивности [4], т.е.

L2
I ()  G() aR exp( j2r)   a(z) exp{ j2[r  n(z) z]}dz , 0

(1)

где G() – спектр источника излучении; aR – амплитуда опорной волны; j  1 ;   2 /  – волновое число, определяемое длиной волны ; 2r – оптическая длина пути опорной волны; a(z) – амплитуда пред-

метной волны, отраженной на глубине z в исследуемом образце в диапазоне глубин L; n(z) – изменение

показателя преломления по глубине среды. Полезная информация содержится в значениях a(z) , характе-

ризующих степень отражения предметной волны по глубине исследуемой среды; их необходимо опреде-

лить в результате обработки полученных значений I () .

В результате нормировки выражения (1) относительно спектра источника излучения интерферометрический сигнал в спектральной ОКТ определяется в области волновых чисел формулой (см., например, [5, 6])

2
S()  1 a(z) exp[ j2n0 z  jd ()]dz , 0

(2)

где n0  const – эффективный показатель преломления в среде; d – начальная фаза каждой из гармонических составляющих сигнала, обусловленная дисперсией в среде (для среды без дисперсии d  0 ).
Обозначим A() интеграл под модулем в выражении (2). Стандартный метод вычисления искомой

амплитуды отраженной волны a(z) состоит в применении обратного преобразования Фурье к выделен-

ной информативной составляющей A() сигнала S() :



a(z)   A() exp( j2n0 z)d .

(3)

0

Заметим, что преобразование (3) является непараметрическим, поскольку позволяет получать ре-

зультаты независимо от вида (параметров) преобразуемой функции A() .

Величину A() можно рассматривать как набор отдельных гармонических составляющих, опре-
деляемых параметрами – частотой, амплитудой и начальной фазой. При этом задача сводится к параметрической идентификации и может быть решена с помощью линейной фильтрации Калмана [8, 9]. Алгоритм дискретной фильтрации Калмана основан на рекуррентной процедуре предсказания значения сигнала на последующий шаг на основе информации, имеющейся на предыдущем шаге, с учетом известной параметрической модели сигнала и использовании ошибки предсказания для уточнения значений параметров. Применительно к спектральной интерферометрии это позволяет выбирать требуемое количество длин волн, обеспечивающее необходимую разрешающую способность и быстродействие.

Алгоритм обработки данных

Для описания алгоритма обработки данных перейдем к их дискретному представлению и введем

номер дискретного отсчета k = 0…К–1 в области волновых чисел, причем k  k , где  – шаг изменения волнового числа.

Пусть имеется последовательность отсчетов A(k)  A(k ) сигнала спектральной интерференции и
решается задача рекуррентной идентификации гармонической составляющей с некоторой фиксированной частотой f. В случае модели сигнала вида

H (k)  u cos(2fk  ) ,

(4)

где  – начальная фаза на частоте f, алгоритм рекуррентной идентификации наличия составляющей с

частотой f сводится к фильтрации параметров – амплитуды u и начальной фазы  гармонического сиг-

нала. Модель (4) может быть представлена в виде

H (k)  uс cos(2fk)  us sin(2fk) , а искомые амплитуда u и начальная фаза  определяются как

(5)

u  uc2  us2 ,   arctg(us / uc ) .

(6) (7)

Научно-технический вестник информационных технологий, механики и оптики, 2013, № 6 (88)

51

МЕТОД ДИНАМИЧЕСКОЙ ОБРАБОТКИ ДАННЫХ …

С учетом сказанного, поиск амплитуды и начальной фазы гармонической составляющей сигнала

для каждой частоты f сводится к оцениванию значений uс и us с уточнением оценки при поступлении

каждого отсчета сигнала k. Введем вектор параметров u  uс us T . Тогда выражение для априорной
оценки этого вектора параметров можно записать на k-ом шаге обработки в виде

u(k)  αu(k 1) ,

(8)

где вектор коэффициентов α задает линейный закон изменения амплитуд (при малых приращениях  можно считать, что   1 ). Апостериорная коррекция осуществляется на основе невязки предсказания и наблюдения:

u (k)  u(k)  P(k)[Hobs (k)  H (k)] .

(9)

Для вычислений по формуле (9) требуется определить коэффициент усиления

P(k)  RCT (k)[C(k)RCT (k)  Rn ]1 , где R – ковариационная матрица ошибки априорной оценки параметров; Rn – дисперсия шума наблюдения; С(k) – матрица перехода, элементы которой являются произ-

водными модели (5) по параметрам, C(k)  Hu (k) , H (k) – прогноз наблюдения, получаемый подстанов-
кой предсказанных параметров в модель (5). Если шаг дискретизации  отнормировать и принять равным единице, алгоритм фильтрации
(8)–(9) сведется к следующему рекуррентному выражению:

u(k)  u(k 1)  RCT (k)[C(k)RCT (k)  Rn ]1[Hobs (k)  H (k)] , где

(10)

C(k )



  

cos(2fk ) sin(2fk )

T  

,

(11)

H (k)  uc (k 1) cos(2fk)  us (k 1) sin(2fk) .

(12)

Следует отметить, что на каждом шаге фильтрации возможен пересчет ковариационной матрицы

R согласно выражению [5]

R(k)  [J  P(k)C(k)]R(k 1) ,

(13)

где J – единичная матрица. Использование (13) позволяет дополнительно повысить скорость сходимости алгоритма, однако увеличивает вычислительную сложность.
Из выражения (3) видно, что a(z) может быть представлена в виде набора δ-функций,

M 1
a(z)  aiδ(z  zi ) ,

(14)

i0

соответствующих гармоническим составляющим спектра и характеризующих положение слоев иссле-

дуемой среды на глубине zi, что эквивалентно представлению в частотной области (области длин волн) в виде M гармонических составляющих. С учетом выражения (14) возможна одновременная идентифика-

ция этих составляющих при использовании параллельной обработки (см. [9]).

В качестве примера работы алгоритма (10)–(13) на рис. 1, а, представлен исходный сигнал со сле-

дующими параметрами: нормированная частота f = 0,01, амплитуда 10 усл. ед., начальная фаза 2 рад,

аддитивный белый гауссовский шум с нулевым средним и дисперсией 1 усл. ед. На рис. 1, б, в, представ-

лены результаты оценивания параметров гармонической составляющей при использовании рассмотрен-

ного выше рекуррентного алгоритма. Видно, что после примерно 100 отсчетов (один период сигнала)

оценки выходят на квазистационарные значения, соответствующие истинным. Незначительные флуктуа-

ции оценок обусловлены наличием шума наблюдения с ненулевой дисперсией, что имеет место в реаль-

ных сигналах. Из гистограмм ошибок оценок амплитуды и фазы (рис. 1, д, е) видно, что они имеют од-

номодовый характер с нулевым математическим ожиданием и близки к гауссовской форме.

Поскольку для каждого сигнала оцениваемые параметры – константы, а их оценки принимают

квазистационарные значения в пределах одного периода сигнала, при реализации алгоритма (10)–(13)

используется критерий останова, предложенный в работе [8]. Указанный критерий основан на определе-

нии разности оценок амплитуды сигнала с заданной частотой на текущем k и предыдущем k – 1 отсчетах

и сравнении с некоторым пороговым значением  . Тогда критерий останова выражается в форме

u(k)  u(k 1)   .

(15)

Выбор порогового значения  может осуществляться различными способами и зависит от пред-
ставления исходного сигнала. Если значения сигнала принадлежат множеству действительных чисел, то в качестве  целесообразно выбирать среднее квадратичное отклонение шума наблюдения. В системах
ОКТ исходные данные, как правило, представляют собой изображения, представленные 256 или 4096 градациями серого для 8- и 12-битного представления соответственно. В этом случае в качестве  целе-

52 Научно-технический вестник информационных технологий, механики и оптики,
2013, № 6 (88)

М.А. Волынский, И.П. Гуров

сообразно выбирать минимально возможное изменение градаций серого. Особенности выбора критерия (15) рассмотрены в работах [8, 9].

Оценка амплитуды сигнала, усл. ед.

Значение сигнала,
усл. ед.

15

25 20

5 15

–5

10 5

–15 0

50 100 150 200

Номер дискретного отсчета

а

0 50 100 150 200 Номер дискретного отсчета б

Значение погрешности оценки сигнала, усл. ед.

Оценка начальной фазы сигнала, рад.

4
3 2
1
0 50 100 150 200 Номер дискретного отсчета в
20
15
10

1,5 1
0,5 0
–0,5 –1 0 50 100 150 200 Номер дискретного отсчета г
15
10

Вероятность, %

Вероятность, %

55

0 –2,8 –1,4

0

1,4 2,8

Погрешность оценки амплитуды,

0,01 усл. ед.

д

0

–18,7 –9,33 0

9,33 18,7

Погрешность оценки фазы,

0,0001 рад

е

Рис. 1. Исходный сигнал с частотой f = 0,01 (а), оценки амплитуды (б) и начальной фазы (в), погрешность предсказания значения сигнала (г) и гистограммы ошибок предсказания амплитуды (д) и начальной фазы (е) для 104 реализаций сигнала

Найденные в (6) и (7) значения амплитуд и начальных фаз могут быть использованы для компенсации влияния дисперсии в среде [5, 6]. Следует отметить, что начальные фазы для гармонических составляющих с амплитудой, близкой к нулю, не имеют физического смысла, поэтому целесообразно ввести пороговое значение оценки амплитуды (например, незначительно превышающее среднее квадратичное отклонение шума наблюдения) и считать начальные фазы гармонических составляющих с амплитудой, меньшей этого порогового значения, равными нулю. Значение начальной фазы для гармонических составляющих с ненулевой амплитудой может быть использовано для компенсации влияния дисперсии в среде с помощью алгоритма, предложенного в работе [5].
Как отмечается в [5], из-за влияния дисперсии происходит уширение полос в спектре интерферометрического сигнала, причем истинному положению полосы соответствует нулевая производная начальной фазы. С учетом этого можно ввести коэффициент пропорциональности на основе найденных значений начальной фазы на каждой частоте и корректировать амплитуду А-скана в соответствии с выражением

aˆ ( z )



1



1 

d

( z ) dz

 

a(z)

,

(16)

где aˆ(z) – скорректированная с учетом влияния дисперсии амплитуда.

В качестве примера на рис. 2, а, показан сигнал спектральной интерференции, полученный при исследовании двухслойного фантома из поливинилхлорида [10] с помощью спектрального оптического когерентного томографа Hyperion (ThorLabs, США) с источником излучения с центральной длиной волны 930 нм [11]. На рис. 2, б, показан результат обработки с компенсацией дисперсии и без компенсации. Видно, что наличие дисперсии в среде уширяет линии в спектре сигнала, а поправка (16) на основе учета значений начальной фазы позволяет скомпенсировать этот эффект, повышая, тем самым, разрешающую способность по глубине. Толщина верхнего слоя в фантоме при его изготовлении составляла около 100 мкм, что иллюстрируется на рис. 2, б.

Научно-технический вестник информационных технологий, механики и оптики, 2013, № 6 (88)

53

МЕТОД ДИНАМИЧЕСКОЙ ОБРАБОТКИ ДАННЫХ …

Значение сигнала, отн. ед.
Коэффициент отражение, отн. ед.

1 0,8
0,6 0,4
0,2 0 850 890 930 970 Длина волны, нм
а

1 0,5

0 100 200 Глубина, мкм б

300

Рис. 2. Сигнал спектральной интерференции (а) и полученный из него А-скан (б) с компенсацией дисперсии (сплошная линия) и без компенсации (пунктирная линия)

Обработка биомедицинских данных

На рис. 3 показан пример В-скана кожи на подушечке указательного пальца человека, полученный с помощью описанного в настоящей работе метода без учета начальных фаз (рис. 3, а, б) и с компенсацией дисперсии (рис. 3, в, г). Видно, что без учета начальных фаз на реконструируемых В-сканах присутствуют артефакты в виде вертикальных полос, вызванные уширением полосы частот, соответствующей каждой глубине объекта, вследствие влияния дисперсии в среде. На рис. 3, в, г, указанные артефакты отсутствуют.

аб
вг Рис. 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.
54 Научно-технический вестник информационных технологий, механики и оптики,
2013, № 6 (88)

М.А. Волынский, И.П. Гуров

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. 233981– 23413.
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).

Волынский Максим Александрович Гуров Игорь Петрович

– Россия, Санкт-Петербург, Санкт-Петербургский национальный иссле-

довательский университет информационных технологий, механики и

оптики,

кандидат

технических

наук,

доцент,

maxim.volynsky@gmail.com

– Россия, Санкт-Петербург, Санкт-Петербургский национальный иссле-

довательский университет информационных технологий, механики и

оптики, доктор технических наук, профессор, зав. кафедрой,

gurov@mail.ifmo.ru

Научно-технический вестник информационных технологий, механики и оптики, 2013, № 6 (88)

55