Приборы и техника эксперимента, 2020, № 1, стр. 136-143

ВВЕДЕНИЕ

Работа посвящена задачам практической реализации биомедицинских диагностических приборов, использующих методы спектральной оптической когерентной томографии (о.к.т.). Спектральная о.к.т. [1] известна как метод неинвазивного исследaования приповерхностных внутренних структур биологических тканей (in vivo). Этот метод основан на зондировании исследуемого объекта маломощным низкокогерентным инфракрасным излучением. В данном методе осуществляется прием интерферометрическим способом инфракрасного излучения, рассеянного в обратном направлении внутренними неоднородностями исследуемой ткани. Последующая математическая обработка оптического спектра интерферометрического сигнала позволяет получить представление о пространственном распределении рассеивающих неоднородностей внутри исследуемого образца. Метод спектральной о.к.т. позволяет обеспечить сравнительно высокую скорость исследования (в данном случае >20 000 A-сканов/с) и высокое пространственное разрешение (~10 мкм).

В этой статье рассматриваются о.к.т.-системы, разработанные в ИПФ РАН. Они предназначены для эндоскопических исследований и для исследований приповерхностных тканей живого организма. В частности, в качестве иллюстраций в данной статье используются результаты о.к.т.-исследования барабанной полости уха человека.

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

СТРУКТУРА СПЕКТРАЛЬНОЙ О.К.Т.-СИСТЕМЫ

Данные о.к.т.-системы основаны на использовании одномодовой волоконной оптики и соответствующих оптических элементов. В основе схемы, представленной на рис. 1, лежит тандемное сочетание измерительного интерферометра Физо и компенсирующего интерферометра Майкельсона. Такая схема известна в литературе как common path OCT [2]. Излучение инфракрасного диапазона при распространении в исследуемой среде испытывает обратное рассеяние и возвращается в волоконный тракт устройства, где смешивается с опорным излучением, отраженным от торца оптического волокна в сканирующем зонде.

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

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

С целью исключения необходимости процедуры передискретизации регистрируемых отсчетов (для ликвидации их неэквидистантности по оптической частоте) в схеме рис. 1 применяются специальные призмы [3, 4]. Данный метод позволяет использовать быстрое преобразование Фурье (б.п.Ф.) при дальнейшей математической обработке.

Для регистрации спектральных компонент используется сенсор на основе линейки InGaAs-фотоэлементов.

Функции обработки и управления в схеме рис. 1 выполняет стандартный персональный компьютер (ПК). Для связи c ПК через стандартный интерфейс USB используется специальная схема сбора данных и управления [5, 6]. Для получения двумерного томографического среза процедура сбора данных синхронизирована с движением зеркала сканирующей системы.

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

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

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

где Z – оптическая разность хода опорной и рассеянной волны, c – скорость света, ω – круговая частота оптического излучения. В выражении (1) выделены автокорреляционная (${{S}_{{AC}}}({\omega })$) и кросс-корреляционная (${{S}_{{CC}}}({\omega })\cos (2{\omega }Z{\text{/}}c)$) составляющие [7].

Если в качестве исследуемого объекта представить только одну границу раздела двух сред с разными показателями преломления (идеализированный случай), то зарегистрированный в этом случае сигнал $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{I} ({\omega })$, будет иметь форму, представленную на рис. 2а. Для получения A-скана (элемента о.к.т.-изображения), соответствующего данному случаю (рис. 2б), из сигнала $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{I} ({\omega })$ следует исключить автокорреляционную составляющую ${{S}_{{AC}}}({\omega })$ и применить преобразование Фурье [7]. Следует отметить, что на рис. 2б в качестве значений F(Z) представлены значения модуля фурье-преобразования только для положительных значений Z.

Исключить автокорреляционную компоненту ${{S}_{{AC}}}({\omega })$ из $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{I} ({\omega })$ можно несколькими способами [7–11]. В основном автокорреляционная составляющая (показана стрелкой на рис. 2а) повторяет форму спектра источника. Но кроме этого она содержит когерентные помехи и результаты нежелательных переотражений в оптическом тракте о.к.т.-системы [7, 8]. Эти факторы являются причиной возникновения автокорреляционных артефактов, присущих спектральным о.к.т.-системам. Поэтому эффективное исключение автокорреляционной составляющей должно максимально подавлять автокорреляционные артефакты и минимально влиять на “полезный” сигнал.

Для решения этой и некоторых других задач в интерферометрической схеме используемой о.к.т.-системы (рис. 1) присутствует модулятор. С его помощью дополнительная оптическая длина пробега опорной волны для каждого следующего A-скана изменяется так, чтобы разность фаз между кросс-корреляционными составляющими соседних A‑сканов была π/2. Другими словами, регистрируемый сигнал в такой о.к.т.-системе имеет вид

где n – текущий номер А-скана. Результат применения данной аппаратной процедуры при регистрации изображения биологического объекта (барабанной перепонки) продемонстрирован на двумерной диаграмме рис. 3.

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

Яркость пикселя на диаграмме соответствует значению $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{I} ({\omega })$ в данной точке. Вертикальные столбцы данной диаграммы – это исходные данные A-сканов изображения. Их совокупность составляет B-скан, в представленном на рис. 3а случае число A-сканов – 512. Частота следования A‑сканов составляет ~20 000 Гц. Число элементов в линейке (число спектральных отсчетов) – 512.

Рис. 3б – это результат двумерного б.п.Ф. данных, представленных на рис. 3а. Здесь яркость пикселей соответствует логарифму модуля комплексного результата двумерного б.п.Ф. Область “полезного сигнала” на рис. 3б сосредоточена в левой верхней и правой нижней четверти диаграммы. В данном случае правая и нижняя области диаграммы – это области “отрицательных частот”. Симметрия областей полезного сигнала – это следствие эрмитовости фурье-образа. Таким образом, данные в правой нижней четверти диаграммы не содержат дополнительной информации и являются комплексно-сопряженными данным из левой верхней части. Область “полезного сигнала” смещена по горизонтали на четверть всего диапазона, так как каждый соседний А-скан регистрировался с возрастающим значением фазы интерференционного сигнала (за счет модуляции) на величину π/2, что соответствует λ/4. Такая модуляция производит эквивалентное смещение частоты сигнала из нулевого положения.

Для получения о.к.т.-изображения в данном случае можно совершить следующие преобразования. Исходные данные (рис. 3а) представим как ${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{I} }_{{w,x}}}$, где w – соответствует номеру спектрального отсчета, x – соответствует номеру A-скана в B-скане. Для каждой строки двумерного массива ${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{I} }_{{w,x}}}$ реализуем преобразование Гильберта следующим способом:

Результат преобразований (3)–(6) ${{I}_{{w,x}}}$ – это комплексный аналитический сигнал, в котором подавлена автокорреляционная составляющая.

Для синтеза о.к.т.-изображения достаточно получить фурье-образ столбцов двумерного массива ${{I}_{{w,x}}}$ т.е.

Таким образом, ${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{F} }_{{Z,x}}}$ – это двумерное о.к.т.-изображение (массив значений, представленный в децибелах).

Результат визуализации ${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{F} }_{{Z,x}}}$ представлен на рис. 4а. Для сравнения на рис. 4в представлен результат простой совокупности фурье-образов исходных A-сканов (при подстановке в (7) данных столбцов ${{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{I} }_{{w,x}}}$ вместо ${{I}_{{w,x}}}$). Для иллюстрации улучшения соотношения сигнал/шум на рис. 4б представлен результат преобразований (3), (4), (6), т.е. проведено только подавление автокорреляционной составляющей, преобразование (5) не применяется. Гистограмма, представленная на рис. 4г, свидетельствует об уменьшении фонового шума (зона I) примерно на 3 дБ при сохранении уровня “полезного сигнала” (зона II, обе зоны помечены на рис. 4а и 4б). Таким образом, применяемые преобразования (3)–(6) улучшают соотношение сигнал/шум и устраняют автокорреляционные артефакты. Необходимо отметить, что на рис. 4 визуализируются только “положительные” значения координаты Z (в терминах б.п.Ф.). Отображение “отрицательного” диапазона Z, соответствующего большим глубинам, для применяемой о.к.т. системы не целесообразно, так как глубже теряется фокусировка сканирующего пучка.

Следует отметить значение модуляции длины пробега опорной волны (2) для применяемых преобразований. Данный прием (в терминах двумерной диаграммы рис. 3б) позволяет переместить спектральную область “полезного сигнала” из зоны “постоянной составляющей”. В этом случае автокорреляционная и кросскорреляционная компоненты разнесены в фурье-пространстве, что дает возможность отфильтровать “полезный сигнал” c минимальными потерями. Кроме того, в результате такой модуляции длины пробега опорной волны диаграмма рис. 3б становится антисимметричной. За счет этого применяемое преобразование Гильберта компенсирует эрмитовость фурье-образа.

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

где N = 512 – число элементов в линейке (число спектральных отсчетов), Δω – ширина полосы регистрации спектрометра. В случае, когда применяются только преобразования (3), (4), (6) – рис. 5б, на изображении наблюдается перевернутое (“зеркальное” [7]) изображение объекта. При плавном приближении или удалении зонда относительно объекта основное и “зеркальное” изображения [7] (в случае рис. 5б) движутся в противоположных направлениях, что сбивает с толку неопытного пользователя. Применение преобразований (3)–(6) позволяет подавить “зеркальную” компоненту (рис. 5а), в результате чего на изображении виден только край барабанной перепонки в нижней части изображения, для которого условие дискретизации по разности хода не нарушено.

КОМПЕНСАЦИЯ ИСКАЖЕНИЙ, ВОЗНИКАЮЩИХ ПРИ БЫСТРОМ ПЕРЕМЕЩЕНИИ ЗОНДА

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

На рис. 6 представлена диаграмма двумерного б.п.Ф. исходных данных B-скана, зарегистрированного во время быстрого удаления зонда от барабанной перепонки. В данном случае области полезного сигнала переместились по горизонтали диаграммы за счет эффекта Доплера, обусловившего возникновение дополнительного фазового набега для кросскорреляционной составляющей сигнала (2). При этом в верхней половине диаграммы “полезный сигнал” оказался в области “отрицательных” пространственных частот. За счет этого в случае применения преобразования (3)–(6) о.к.т.-изображение будет выглядеть, как показано на рис. 7.

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

Для разрешения данной проблемы предлагается компенсировать влияние скорости перемещения зонда на смещение фазы спектральных компонент пространственных частот.

Для получения информации о текущей скорости перемещения зонда относительно исследуемого объекта можно воспользоваться частичной коррелированностью соседних A-сканов. Для используемой о.к.т.-системы характерно наличие частичного перекрытия соседних A-сканов, проиллюстрированное на рис. 8.

Это позволяет провести следующие манипуляции. Пусть ${{F}_{{Z,x}}}$ – это двумерный массив фурье-образов исходных A-сканов, т.е.

Тогда совокупную разность фаз между всеми соседними элементами двух соседних A-сканов

можно использовать для компенсации текущей скорости перемещения зонда ($F_{{Z,x + 1}}^{*}$ – комплексно-сопряженное с $F_{{Z,x + 1}}^{{}}$). Для этого следует вычислить массив кумулятивных фазовых множителей:

Коррекция исходных данных текущего B-скана будет выглядеть как

Эффективность данной процедуры определяется тем, что преобразование (8) является достаточно устойчивым к влиянию шумов и некоррелированных компонент соседних A-сканов. В качестве иллюстрации этого можно продемонстрировать результат двумерного б.п.Ф. от $I_{{w,x}}^{'}$ для случая, приведенного на рис. 6 и 7 (см. рис. 9).

В данном случае зона “полезного сигнала” всегда находится в одном месте диаграммы (рис. 9) – в зоне координаты X = 0, вне зависимости от скорости перемещения зонда.

Для получения о.к.т.-изображения в этом случае можно применить следующие манипуляции. Пусть $F_{{Z,X}}^{{''}}$ – результат двумерного б.п.Ф. комплексных скорректированных данных $I_{{w,x}}^{'}$ (рис. 9), тогда

Результат преобразований (9), (10) представлен на рис. 10.

Следует отметить, что преобразования (9), (10) в данном случае аналогичны преобразованию Гильберта, так как вследствие применения модуляции длины пробега опорной волны исходный зарегистрированный сигнал состоял из последовательности квадратурных спектральных компонент, которые были смещены по фазе процедурой коррекции.

ЗАКЛЮЧЕНИЕ

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

Структура описанных выше математических манипуляций позволяет использовать их в рамках метода асинхронных параллельных вычислений, описанных в [12], что обеспечивает высококачественную визуализацию исследуемых объектов в реальном времени.