Рекуррентная обработка данных в спектральной оптической когерентной томографии на основе фильтрации Калмана Текст научной статьи по специальности «Электротехника, электронная техника, информационные технологии»

2

УДК 681.787:[517.977+517.443]

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

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

Предложен метод динамической обработки данных в спектральной оптической когерентной томографии на основе фильтрации Калмана. Показана возможность адаптивного выбора необходимого количества спектральных отсчетов для обеспечения требуемого разрешения для конкретного исследуемого образца. Приведены примеры обработки экспериментальных данных, полученных при исследовании биологических объектов.

Ключевые слова: оптическая когерентная томография, спектральная интерферометрия, преобразование Фурье, фильтр Калмана.

Введение

Бесконтактный контроль внутренней микроструктуры объектов необходим для многих областей науки и современных технологий. Одно из перспективных направлений представляет оптическая когерентная томография (ОКТ), широко используемая в биомедицине [1-3]. Методы ОКТ основаны на принципах интерферометрии малой когерентности применительно к рассеивающим неоднородным средам и обеспечивают разрешение до долей микрометра при восстановлении томографических изображений, представляющих трехмерную внутреннюю микроструктуру неоднородных сред на глубине проникновения оптического излучения.

Как известно, методы ОКТ основываются на использовании корреляционного или спектрального интерферометра [1-5]. В корреляционном интерферометре осуществляют перемещение оптической системы относительно исследуемого объекта. При этом интерференционные полосы малой когерентности формируются в пределах длины когерентности излучения при интерференции части измерительной волны, отраженной от поверхности непрозрачного объекта или от слоя частично прозрачного неоднородного объекта, находящихся от светоделителя на расстоянии, равном оптической длине пути опорной волны.

В спектральном интерферометре оптическая длина пути опорной волны не равна оптической длине пути измерительной волны для всего диапазона высот рельефа непрозрачного объекта или исследуемой области частично прозрачного объекта. На выходе интерферометра размещен спектральный прибор, позволяющий определить составляющую отраженной измерительной волны для каждой из длин волн. При этом, используя преобразование Фурье-спектра, зарегистрированного приемником излучения, можно определить расстояние и степень отражения от каждого слоя [5]. Основное преимущество методов спектральной ОКТ состоит в отсутствии необходимости перемещения оптико-механических элементов, что обеспечивает более высокое быстродействие по сравнению с корреляционной ОКТ.

Вместо использования спектрального прибора, размещенного на выходе интерферометра, можно освещать интерферометр монохроматическим излучением с перестраиваемой последовательно во времени длиной волны. Преимущество метода спектральной ОКТ с перестраиваемой длиной волны состоит в возможности получения информации о микроструктуре объекта на всем наблюдаемом участке, а не только вдоль линии, соответствующей щели спектрального прибора.

В ОКТ с перестраиваемой длиной волны особенно важно обеспечить динамическую обработку получаемых данных, задавая шаг по длине волны и количество длин волн в зависимости от требований к разрешающей способности и быстродействию для конкретного исследуемого объекта.

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

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

ФОТОНИКА И ОПТОИНФОРМАТИКА

Формирование сигналов спектральной интерференции

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

сивности (см., например, [5])

I (k) = G(k)

aR exp(-j2kr) + Ja(z)exp{-j2k[r + ns(z)z]}dz

(1)

где G(к) - спектр источника излучения; aR - амплитуда опорной волны; j =V-1; к = 2 л / X - волновое число, определяемое длиной волны X; 2r - оптическая длина пути опорной волны; a(z) - амплитуда предметной волны, отраженной на глубине z в исследуемом образце в диапазоне глубин L; ns (z) - изменение показателя преломления по глубине среды. Полезная информация содержится в значениях a(z), характеризующих степень отражения предметной волны по глубине исследуемой среды, и которые необходимо определить в результате обработки полученных значений I (к).

В результате нормировки выражения (1) относительно спектра источника излучения и в предположении отсутствия дисперсии в среде (т.е. ns (z) = const) интерферометрический сигнал в спектральной ОКТ определяется в области волновых чисел формулой

S (к) =

1 + J a( z) exp(- j 2knsz )dz

(2)

Обозначим как A(k) интеграл под модулем в выражении (2). В результате разложения A(k) на вещественную и мнимую части, A(k) = x + jy, получим:

S(k) =| 1 + A(k) |2 =| 1 + x + jy |2 = 1 + (x2 + y2) + (x + jy) + (x - jy) = 1 +1 A(k) |2 + A(k) + A*(k). (3)

Стандартный метод вычисления искомой амплитуды предметной волны a(z) состоит в применении обратного преобразования Фурье к выделенной в выражении (3) информативной составляющей A(k) сигнала S(k), k > 0 :

да

a(z) =J A(k)exp(j2knsz)dk . (4)

0

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

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

Пусть имеется последовательность значений A(n) и решается задача рекуррентной идентификации некоторой гармонической составляющей с фиксированной частотой f содержащейся в этой последовательности. В рассматриваемом случае n обозначает номер в последовательности значений волнового числа k, задаваемых в спектральной ОКТ.

Искомая гармоническая составляющая определяется моделью H(n) = u cos(2nfnAk), где u - амплитуда; Ак - шаг дискретизации в области волновых чисел, соответствующий шагу изменения длины волны. Рекуррентный алгоритм идентификации гармонической составляющей сводится к предсказанию ее амплитуды u для шага n:

u (n) =au(n -1), (5)

где коэффициент a задает линейный закон изменения амплитуды между соседними отсчетами (при постоянной амплитуде следует считать, что a = 1 ), и ее апостериорной коррекции на основе невязки предсказания и наблюдения:

й(n) = u (n) + P(n) [A(n) - C(n)u (n)]. (6)

В формуле (6) коэффициент перехода C (n) является производной модели по искомому параметру, C(n) = H (n), коэффициент усиления P(n) = RprC(n)[C2(n)Rpr + Rn , где Rpr - дисперсия ошибки априорной оценки оцениваемого параметра; Rn - дисперсия шума наблюдения [6].

2

2

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

ч 1Ч КС(п)[А(п) - С(п)и(п -1)]

и(п) = и (п-1) + —---, (7)

С (п)Ярг + Я„ ' ^

где С(п) = со&(2п/п).

Значения Крг и Яп выбираются исходя из априорной информации о сигнале и шуме. Поскольку

отклонение начальной оценки искомого параметра от истинного значения амплитуды влияет лишь на скорость сходимости фильтра и не влияет на конечный результат, начальное значение оценки амплитуды может быть установлено в известной степени произвольно. Исследования показали, что наличие ненулевой начальной фазы также не влияет на результат оценивания амплитуды, а влияет только на скорость сходимости. Ввиду этого можно пренебречь дополнительной оценкой начальной фазы.

Примеры обработки сигналов в спектральной интерферометрии

Для проверки работы алгоритма (7) был использован модельный сигнал спектральной интерференции (рис. 1), соответствующий 500 значениям длины волны. Введем определение частоты как величины, обратной количеству отсчетов на одном периоде гармонической составляющей. Отметим, что при этом для данной последовательности минимальное значение частоты составляет 0,002, что соответствует одному периоду на всей последовательности. Сигнал на рис. 1 содержит две гармонические составляющие и описывается моделью

s(n) = со^к/п) + 0,5 ооб(2п/2п) + м>(п), (8)

где Л = 0,01; /2 = 0,1; ^(п) - белый гауссовский шум с дисперсией 0,01.

100 200 300 400 Номер дискретного отсчета а

100 200 300 400 500 Номер дискретного отсчета б

Рис. 1. Сигнал спектральной интерференции (а) и оценки амплитуд гармонических составляющих сигнала

с частотами 0,01, 0,1 и 0,05 (б)

& 1

£0,8

§ 0,6

* 0,4

<4

£ 0,2

<Ц

Я

О П

0,05

0,10 0,15 Частота и

0,19 0,24 0,29

0,09 0,095 0,1 0,105 Частота и б

0,11

Рис. 2. Результат идентификации гармонических составляющих с частотами из интервала [0,002, 0,3] (а) и его увеличенный фрагмент в окрестности частоты 0,1 (б)

При идентификации гармонических составляющих сигнала (8) была введена дополнительная (гипотетическая) составляющая с частотой 0,05, отсутствующая в этом сигнале, для определения помехоустойчивости алгоритма идентификации в условиях воздействия шума. Видно, что после обработки 50 отсчетов (что соответствует половине периода гармонической составляющей с частотой Л) оценки выходят на стационарные значения, равные 1; 0,5; 0 для частот 0,01; 0,1; 0,05 соответственно. Незначительные флуктуации оценок обусловлены наличием шума наблюдения с ненулевой дисперсией, что имеет место в реальных системах ОКТ. Ввиду указанного выше отсутствия влияния начальной фазы на результат идентификации гармонической составляющей анализируемый сигнал может иметь нецелое количество периодов на длине реализации. Как известно, если при преобразовании Фурье используется нецелое число периодов, то в таких случаях появляются артефакты (краевой эффект). При моделирова-

а

нии нами использовано целое число периодов для корректности сравнения характеристик предлагаемого метода с преобразованием Фурье.

На рис. 2 представлен спектр сигнала, показанного на рис. 1. Идентификация гармонических составляющих проводилась для частот в интервале [0,002, 0,3]. Увеличенный фрагмент спектра в окрестности частоты ^ = 0,1 демонстрирует разрешающую способность метода.

Из рис. 2 видно, что ширина спектра выделенной гармонической составляющей с частотой ^ = 0,1 на уровне 0,2 примерно равна удвоенному минимальному значению частоты (0,002), что соответствует разрешающей способности преобразования Фурье.

Выбор количества отсчетов

При практическом использовании систем спектральной ОКТ важно минимизировать количество отсчетов (длин волн) без снижения разрешения. При использовании дискретного преобразования Фурье при обработке сигналов в ОКТ требуемое количество отсчетов прямо пропорционально разрешающей способности по глубине. При исследовании произвольных объектов априорная информация о «необходимом» разрешении по глубине отсутствует. Быстрая сходимость метода рекуррентной обработки данных на основе линейного фильтра Калмана позволяет вводить критерий останова, т.е. ограничения количества отсчетов при достижении разрешения, достаточного для конкретного образца.

Можно ввести критерий останова на основе разности оценок амплитуды сигнала с заданной частотой на текущем п и предыдущем (п - 1) отсчетах, сравнивая эту разность с некоторым пороговым значением е . Тогда критерий останова выражается выполнением условия

й(п) - й(п -1) <е . (9)

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

На рис. 3 показан пример восстановления двумерного томографического изображения (В-скана) силиконовой модели кровеносного сосуда на 128, 596, 1023 шагах обработки и эволюция оценки одного из столбцов В-скана (А-скана) в процессе обработки. Критерий останова (9) выполняется на 312 отсчете, в то время как достижение аналогичного результата при использовании алгоритма обработки на основе быстрого преобразования Фурье требует 1024 отсчета. Таким образом, в данном примере требуемое количество отсчетов (длин волн) снижено примерно в три раза без снижения качества представления томограммы.

г--~— ~

100[ 200 300 400 5001 Номер дискретного отсчета, п

а б

Рис. 3. Оценка В-скана тест-объекта на 128, 596 и 1023 шагах обработки (а) и эволюция оценки А-скана при поступлении новых отсчетов сигнала (б). По вертикальной оси отложены инвертированные значения А-скана в зависимости от номера дискретного отсчета сигнала. Ширина В-сканов - 2 мм,

глубина (размер А-скана) - 1,5 мм

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

На рис. 4 показаны примеры В-скана личинки серой мясной мухи Sarcophagidae, полученные с помощью описанного в настоящей работе метода при использовании 98 отсчетов (длин волн) и с помощью метода на основе преобразования Фурье при 1024 отсчетах. Отличие соответствующих элементов изображений не превышает 12 градаций яркости из 255, т.е. погрешность составляет не более 5%. Отме-

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

а б

Рис. 4. B-сканы личинки серой мясной мухи Sarcophagidae, полученные с помощью описанного в настоящей работе метода (а) и с помощью метода на основе преобразования Фурье (б). Ширина В-сканов - 2 мм, глубина - 1 мм

Заключение

В системах спектральной оптической когерентной томографии с перестраиваемой длиной волны актуальна динамическая обработка данных, позволяющая сократить количество спектральных отсчетов (длин волн) при сохранении высокого разрешения. Традиционные методы обработки сигналов спектральной интерференции, основанные на преобразовании Фурье, не удовлетворяют этому требованию. Для повышения быстродействия систем спектральной оптической когерентной томографии в работе предложен алгоритм обработки данных на основе линейного фильтра Калмана, позволяющий проводить идентификацию гармонических составляющих с заданным набором частот и их вклад (амплитуду гармоник). Введенный критерий останова алгоритма позволяет использовать для восстановления полезной информации о внутренней микроструктуре исследуемого образца минимально необходимое количество спектральных отсчетов (длин волн) при обеспечении требуемого разрешения по глубине.

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

Описанный в работе алгоритм реализован в форме программного модуля [7]. Авторы выражают благодарность М.В. Волкову за полезное обсуждение особенностей программной реализации описанного в работе метода.

Работа выполнена при финансовой поддержке Министерства образования и науки Российской Федерации.

Литература

1. Tomlins P.H., Wang R.K. Theory, developments and applications of optical coherence tomography II J. Phys. D: Appl. Phys. - 2005. - V. 3S. - P. 2519-2535.

2. Optical Coherence Tomography. Technology and Applications I W. Drexler, J.G. Fujimoto, eds. - BerlinHeidelberg: Springer-Verlag, 200S. - 1346 p.

3. Гуров И.П., Волынский М.А., Жукова Е.В., Маргарянц Н.Б. Исследование растительных тканей методом оптической когерентной микроскопии II Научно-технический вестник информационных технологий, механики и оптики. - 2012. - № 5 (S1). - С. 42-47.

4. Волынский М.А., Воробьева Е.А., Гуров И.П., Маргарянц Н.Б. Бесконтактный контроль микрообъектов методами интерферометрии малой когерентности и оптической когерентной томографии II Изв. вузов. Приборостроение. - 2011. - Т. 54. - № 2. - С. 75-S2.

5. Häusler G., Lindner M.V. «Coherence radar» and «Spectral radar» - new tools for dermatological diagnostics II J. Biomed. Opt. - 199S. - V. 3. - P. 21-31.

6. Kalman R.E. A new approach to linear filtering and prediction problems II Trans. ASME, J. Basic Eng. -1960. - V. S2. - P. 35-45.

7. Волынский М. А. Программный модуль «Рекуррентная обработка данных в спектральной оптической когерентной томографии» II Свидетельство № 2012617SS1 от 31.0S.2012.

Волынский Максим Александрович - Санкт-Петербургский национальный исследовательский университет

информационных технологий, механики и оптики, кандидат технических наук, доцент, maxim.volynsky@gmail.com Гуров Игорь Петрович - Санкт-Петербургский национальный исследовательский университет

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