CC BY
78
К 110-летию СПбГУ ИТМО
УДК 621.391: 681.787
Е. А. Воробьева, И. П. Гуров, М. В. Петерсон
ИССЛЕДОВАНИЕ МЕТОДА ОБРАБОТКИ СИГНАЛОВ В СИСТЕМАХ ОПТИЧЕСКОЙ КОГЕРЕНТНОЙ ТОМОГРАФИИ С ПОВЫШЕННЫМ БЫСТРОДЕЙСТВИЕМ
Рассмотрен метод асинхронной амплитудной демодуляции сигналов в оптической когерентной томографии (ОКТ), основанный на выполнении последовательных операций квадратичного преобразования информативной составляющей интер-ферометрического сигнала и фильтрации нижних частот. Исследована погрешность восстановления огибающей сигнала при использовании субдискретизации с прореживанием последовательности отсчетов сигнала в условиях влияния аддитивного шума. Представлены результаты восстановления огибающей при использовании фильтра нижних частот Баттерворта для выделения полезной составляющей сигналов в ОКТ.
Ключевые слова: оптическая когерентная томография, субдискретизация сигналов, амплитудная демодуляция.
Введение. Исследования внутренней структуры различных объектов имеют важное значение для биомедицины, материаловедения и высоких технологий. Одним из переспектив-ных методов неразрушающей диагностики микроструктуры объектов является оптическая когерентная томография (ОКТ) [1—8], обеспечивающая высокую разрешающую способность при восстановлении томографических изображений, представляющих трехмерную внутреннюю микроструктуру неоднородных объектов и сред на глубине проникновения оптического излучения. В обзорных статьях [6—8] подробно описаны основные принципы ОКТ и результаты, достигнутые в последние годы.
Одно из важных направлений исследований и разработок — создание систем ОКТ с повышенным быстродействием. Ввиду высокой разрешающей способности ОКТ необходимо обрабатывать большие объемы информации. В последние годы ведутся разработки систем ОКТ, обеспечивающих в режиме реального времени наблюдение и анализ трехмерных микрообъектов и неоднородных сред, изменяющих свое состояние, что особенно актуально для биомедицинских исследований [7].
Известно (см., например, [9, 10]), что наиболее высокую разрешающую способность обеспечивают широкопольные системы ОКТ корреляционного типа. Однако быстродействие таких систем ограничивается характеристиками электронного тракта регистрации и обработки сигналов ОКТ. Повышение быстродействия достигается при использовании метода субдискретизации применительно к узкополосным сигналам [11], когда значение частоты дискретизации ОКТ-сигналов устанавливается ниже предела Найквиста для частоты интерференционных полос [12]. При этом возможно значительно уменьшить количество видеокадров,
регистрируемых в широкопольных системах ОКТ и соответственно снизить требования к системе регистрации и обработки информации.
В работе [13] детально проанализирован метод субдискретизации применительно к классической интерферометрии и продемонстрирована возможность превышения предела Найквиста за счет использования предварительной информации о свойствах интерферомет-рических сигналов. В работах [12, 14] показано, что метод субдискретизации при определенных условиях позволяет минимизировать погрешность обработки сигналов ОКТ до единиц процентов при выигрыше в быстродействии примерно на порядок.
Метод субдискретизации сигналов ОКТ, предложенный и исследованный в работах [12, 14], включает операцию фильтрации параметров сигналов ОКТ нелинейным фильтром Калмана, что обеспечивает возможность обработки сигналов с изменяющимися значениями несущей частоты интерференционных полос, обусловленными влиянием внешних дестабилизирующих факторов. Однако этот метод относится к классу методов синхронной демодуляции, поэтому при значительных изменениях или непредсказуемых скачках несущей частоты полос в ряде случаев не обеспечивается необходимая устойчивость обработки сигналов ОКТ.
Для решения поставленной задачи можно совместно использовать метод субдискретизации и асинхронной амплитудной демодуляции сигналов ОКТ. Проведенные исследования показали [15], что предпочтительным вариантом асинхронной амплитудной демодуляции является выполнение последовательных операций устранения фоновой составляющей сигнала ОКТ, возведения в квадрат информативной составляющей интерферометрического сигнала и фильтрации нижних частот.
В настоящей работе получены оценки выигрыша в быстродействии при использовании метода субдискретизации с прореживанием последовательности отсчетов и представлены результаты исследований точности метода [15] в условиях влияния помех.
Алгоритм обработки ОКТ-сигналов. Интерферометрический сигнал малой когерентности в ОКТ, сформированный в фиксированной точке (х, у) плоскости наблюдения, при изменении координаты г по глубине исследуемой среды после удаления фоновой составляющей можно представить в виде
г ) = А( г )соБ Ф( г), (1)
где А( г) и Ф( г) — огибающая и фаза полос соответственно. Огибающая полос содержит полезную информацию о степени отражения измерительной волны по глубине исследуемой среды. Фаза полос Ф( г) выражается как
Ф( г) = 2/ г+ф+5ф( г ), (2)
где /о — несущая частота полос, ф — начальная фаза в точке г = 0, которую можно исключить без ограничения общности рассмотрения, 5ф(г) учитывает случайные изменения фазы, обусловленные, в частности, влиянием внешних факторов. Полезная составляющая сигнала (1) искажается также влиянием аддитивного шума наблюдения п(г), при этом
г) = А( г) соБ Ф( г) + п( г). (3)
Необходимо получить оценки величины А(г) при обработке сигнала ^(г) в (3) с высокой разрешающей способностью, помехоустойчивостью и быстродействием.
Схема алгоритма обработки сигнала [15] приведена на рис. 1.
Сигнал ^(г) обрабатывается в полосно-пропускающем фильтре для устранения фоновой составляющей и высокочастотного шума. Полоса пропускания фильтра устанавливается достаточно широкой для передачи огибающей интерференционных полос с требуемой разрешающей способностью с учетом возможных отклонений несущей частоты и фазы полос. Полезная составляющая интерферометрического сигнала (1) возводится в квадрат:
52 (z) = 2 A2 (z)[1+cos 2Ф(z)] = I A2 (z)+2 A2 (z) cos 2Ф(z). (4)
Второе слагаемое в (4) с учетом (2) изменяется быстрее по сравнению с первым и устраняется при помощи фильтра нижних частот с частотой среза, соответствующей максимальной частоте в спектре квадрата огибающей. Извлечение корня квадратного из оставшегося первого слагаемого дает оценку огибающей интерференционных полос.
Рис. 1
Преимущество алгоритма (см. рис. 1) состоит в его простоте и отсутствии строгих априорных ограничений на значения параметров интерферометрического сигнала. Алгоритм обеспечивает получение точных результатов, если значение несущей частоты полос /0 (с учетом случайных отклонений фазы) является более высоким по сравнению с удвоенным значением частоты Найквиста /N = 2/max, т.е. /0 >4/max, для максимальной частоты /max в спектре огибающей.
При использовании рассматриваемого метода необходимо учитывать свойства алгоритма обработки ОКТ-сигналов при различных соотношениях несущей частоты интерферомет-рического сигнала и частоты дискретизации, включая область субдискретизации по отношению к критерию Найквиста.
Характеристики спектра при субдискретизации сигналов. При дискретизации сигнала (1) получают последовательность отсчетов
s(k) = A(k) cos Ф(к), (5)
где номер отсчета к соответствует дискретному значению координаты zk = кAz, к = 0,1,..., K,
Az — шаг дискретизации.
Рассмотрение процесса взятия отсчетов непрерывного сигнала зависит от способа определения этой операции. В теории дискретизации непрерывных сигналов обычно предполагается (см., например, [11, 16]), что процесс дискретизации состоит в умножении непрерывного сигнала на решетку дельта-функций с шагом Az . При этом, согласно известной теореме о свертке, спектр дискретизованного сигнала представляет собой набор спектров исходного непрерывного сигнала, разнесенных по частоте на величину /s = 1 / Az .
В работе [12] было отмечено, что формирование спектральных порядков при дискретизации имеет место, если последовательность отсчетов априори включает дискретные отсчеты с нулевыми значениями между дельта-функциями. Если субдискретизация проводится с пропуском всех прореженных отсчетов, то наблюдается перенос спектра исходного сигнала на более низкую несущую частоту без формирования других спектральных порядков. Указанное свойство объясняется тем, что при пропуске участков непрерывного сигнала между точками дискретизации масштаб времени задается частотой дискретизации, но не исходным сигналом.
На рис. 2, а, б приведен пример ОКТ-сигнала в виде последовательности отсчетов s(k) и показан модуль спектра сигнала. Сигнал содержит шесть отсчетов на периоде несущей частоты, нормированное значение шага дискретизации Az принято равным единице, при этом значение несущей частоты примерно 0,17. Получаемые значения модуля спектра определяются известным равенством Парсеваля (правилом сохранения энергии) для дискретного представления спектра. При моделировании к полезной составляющей сигнала (см. рис. 2, а) добавлен аддитивный шум n(k) с центрированным нормальным распределением и дисперсией
ап = 0,01, что соответствует среднему квадратическому значению (СКО) шума примерно 3 %
от максимального значения огибающей сигнала.
При прореживании последовательности отсчетов сигнала (взятии каждого пятого отсчета в рассматриваемом примере) и заполнении прореживаемых отсчетов нулевыми значениями (см. рис. 2, в) спектр содержит спектральные порядки (см. рис. 2, г). В случае взятия каждого пятого отсчета и пропуска других (рис. 2, д) спектр сигнала (рис. 2, е) имеет вид, подобный спектру исходного сигнала (рис. 2, б).
а)
«(к), о.е. 1
0
б)
-1
$(/
200 100
0 200 400 600 800 к
I, , , .
-0,4 -0,2
0,2 0,4 / о.е.
в)
s(k), о.е. 1
-1
г)
S/)
40 -I
0
200 400 600 800 к
д)
«(к), о.е. 1
-1
е) $(/
30 20 10 0
■м - ИМЯ км м«
-0,4 -0, 0 0,2 0,4 /, о.
-0,4
-0,2
0,2
0,4 /, о.е.
Рис. 2
Получаемое при субдискретизации новое значение несущей частоты /0 (см. рис. 2, г) определяется известным условием наложения (маскирования) частот (см., например, [17]), которое имеет вид
/0 = 2т/м ± /0, (6)
где т = 1, 2, ...
Обратимся к рис. 2, в котором /0 = 0,17, /« = 0,2 и /^ = 0,1. При субдискретизации (рис. 2, г) при т = 1 наблюдаются значения частоты 0,03 и 0,37; 0,23 и 0,43 (третий и пятый пик от нуля) соответствуют значениям т = 2 и 3, т.е. 2-2/^ -/0 = 2-2-0,1-0,17 = 0,23; 2-3/^ -/0 = 2-3-0,1-0,17 = 0,43 . Отметим, что при введенных выше допущениях значения час-
0
0
0
0
0
0
тоты выше 0,5 в спектре не имеют смысла, поскольку понятие (априори неизвестной) частоты не определено при взятии менее двух отсчетов на периоде.
Рассмотрим частотные соотношения при субдискретизации более подробно. Введем операцию дискретизации с использованием гипотетического опорного сигнала
( -) = ы^п^г) (7)
с частотой . Найдем произведение полезной составляющей сигнала (1) и опорного (7):
ад=А(.-)Со8(2я/0г)с°8(ад=2 д.-ксо^я«, +Л ).-]+С05[2Я(/0-/а).-]}. (8)
Пусть взятие отсчетов осуществляется в точках = к А-, соответствующих максимумам опорного сигнала, т.е. 2п/3-к = 2пк , откуда получаем -к = к // . При этом из (8) получим дискретную последовательность отсчетов в форме
%(к)=2 А(к ){^[2п(/о + / )к / / ]+соа[2я(/о - / )к / / ]} =
=1А (к ){соБ[2л;(/о / / )к+2пк]+^[2л(/о / / )к - 2пк ]} =
= А (к )^[2п(/о/ /8 )к ]. (9)
Режиму субдискретизации соответствует отношение /о / / > о, 5, когда на период несущей частоты /о в (2) приходится менее двух отсчетов. В общем случае при субдискретизации /о //У =т+8, /у = / /2, т > 1 — целое число, 8 — дробная часть, о < 8 < 1. При этом с учетом периодичности и четности функции косинуса значения ~ (к) в (9) одинаковы для частот /о и 2т/у ±/о в соответствии с соотношением (6).
Проведенное рассмотрение спектров при субдискретизации показывает, что во избежание проявления эффекта наложения частот предпочтительно использовать методику прореживания с пропуском промежуточных отсчетов (рис. 2, д, е), при которой снижается общее число отсчетов в выборке. Поскольку значение шага дискретизации А- сигнала (см. рис. 2, д), как и ранее, принято равным единице и при прореживаниии на каждый период исходного сигнала приходится меньшее число отсчетов (что эквивалентно уменьшению общей длительности реализации сигнала), преобразование спектра сигнала определяется известным свойством масштабирования (подобия) для преобразования Фурье
Е{^(ак А-)} = (1/ а)Б (/ / а), (1о)
где обозначает операцию преобразования Фурье, Б(/) — спектр исходного сигнала, а — постоянная (в рассматриваемом примере а = 5). Из (1о) следует, что при прореживании происходит уширение спектра, пропорциональное степени прореживания а, что иллюстрируется на рис. 2, е (ср. с рис. 2, б).
Исследование точности алгоритма обработки при субдискретизации ОКТ-сигналов. Полученные в результате субдискретизации значения сигнала в рассматриваемом алгоритме возводятся в квадрат, согласно (4). При этом имеет место преобразование спектра, которое иллюстрируется на рис. 3, б, г для исходной (рис. 3, а) и прореженной (рис. 3, в) последовательности соответственно.
Из выражения (4) следует, что полезная информация об огибающей сигнала А(к) в (5) содержится в составляющей спектра (см. рис. 3, г) в окрестности нулевой частоты. Для выделения этой составляющей следует использовать фильтр нижних частот с частотой среза, соответствующей максимальной частоте в указанной составляющей спектра. Правильный выбор частоты среза с учетом рассмотренного выше эффекта уширения спектра при субдискретизации позволяет минимизировать погрешности получаемых оценок огибающей.
Из (4) и (10) следует, что частота среза фильтра /с должна удовлетворять условию
/с = 2а/тах. (11)
Кроме того, во избежание наложения составляющих спектра (рис. 3, г) должно быть обеспечено значение несущей частоты интерференционных полос
/0 ^ 40/тах . (12)
Источниками погрешностей оценки огибающей являются наличие второго слагаемого в соотношении (4) и влияние шума п(г) в (3). Заметим, что квадратичное преобразование шума приводит к ненулевому среднему значению, которое не вносит значительных изменений в огибающую при достаточно высоких значениях отношения сигнала к шуму, поэтому в качестве основного источника погрешностей следует оценивать влияние второго слагаемого в (4).
а)
«(к), о.е.
б)
120 80 40 0
0
в)
^(к), о.е.
200 400 600 800 к
г)
1
-0,4 -0,2 0 0,2 0,4 /, о
0,8
0,4
S(f) 30 -Ъ
20
10
0
40
80
120
160
-0,4
-0,2
0,2
0,4 /, о.е.
Рис. 3
(13)
При рассмотрении удвоенного второго слагаемого в (4) с учетом (2)
(г ) = А2 (г ) сов(4п/0 г+25ф) нетрудно показать, что при флуктуациях фазы с нулевым средним значением ^бф = 0 сред-
(14)
нее значение
= А2( г)в соБ 4/г
где
р = ( со8(2бф)>. (15)
Отметим, что именно средние значения сигналов наблюдаются при формировании томограммы, состоящей из набора огибающих сигналов вида (5).
Известно (см., например, [18]), что при гауссовой плотности вероятности отклонений
фазы с дисперсией аф соотношение (15) может быть представлено следующим образом:
в = ехр(-аф /2).
(16)
0
0
к
Таким образом, с учетом (14) и (16) влияние неинформативной составляющей (13) снижается при увеличении флуктаций фазы интерференционных полос и способствует помехоустойчивости рассматриваемого алгоритма асинхронной амплитудной демодуляции.
Влияние случайных отклонений частоты 5/ можно учесть в выражении (13) в форме
s|( z ) = Л2( z )cos[4n(/ +5/ ) z ].
(17)
Известно (см., например, [19]), что мгновенные значения частоты определяются производной от значений фазы (отклонения фазы в реальной физической ОКТ-системе предполагаются дифференцируемыми), причем в случае гауссовой плотности вероятности отклонений фазы плотность вероятности производной также гауссова, и оба значения плотности вероятности статистически независимы в совпадающие моменты времени. Следовательно, случайные отклонения частоты в среднем снижают влияние неинформативной составляющей (13) аналогично флуктуациям фазы.
Проведенное рассмотрение показывает возможность устранения влияния составляющей
(13) при вычислении скользящего среднего с интервалом усреднения Ь«(2а/тах)-1 подобно
тому, как это осуществляется в радиометрии [2о]. Однако при значительных отклонениях несущей частоты интерференционных полос в ОКТ возможно частичное наложение составляющих спектра (см. рис. 3, г) и невыполнение условия (11), что вызывает возрастание погрешности оценки огибающей интерферометрического сигнала. Поэтому целесообразно использовать фильтр нижних частот с более высокой частотной селективностью по сравнению со скользящим средним.
В качестве фильтра нижних частот для выделения полезной составляющей спектра (рис. 3, г) использовался фильтр Баттерворта с передаточной характеристикой вида
\H (/)|2 = Г 1+(/ / / )
2M
(18)
где М — порядок фильтра.
На рис. 4 представлен пример восстановления огибающей ОКТ-сигнала (см. рис. 2, а) с использованием рассматриваемого метода с прореживанием отсчетов (а=5) при частоте среза фильтра, соответствующей условию (11), и величине отношения СКО шума к максимальному значению огибающей сигнала 4 %. Сравнение восстановленной огибающей ОКТ-сигнала (сплошная кривая) и исходной огибающей (пунктир) показывает, что метод характеризуется высокой точностью.
Исследования точности получаемых результатов проводились при моделировании ОКТ-сигналов с получением оценок величины СКО (5А) полученных значений огибающей А(к) от истинных значений в виде
-,1/2
40
80
Рис. 4
120
160
5A =
£ (() - A(k) )2 k 1 k=1
(19)
На рис. 5 показаны кривые зависимости погрешности восстановления огибающей от выбранного значения частоты среза фильтра нижних частот при степени прореживания а=5. Вычисления проводились при усреднении по ансамблю из 100 реализаций ОКТ-сигнала при
1
различных уровнях шума (сплошная кривая — сигнал без шума; пунктир — к сигналу добавлен шум с СКО 4 %; штрихпунктир — к сигналу добавлен шум с СКО 6 %).
Рис. 5
Из рисунка видно, что при малых значениях частоты среза фильтра нижних частот имеют место существенные погрешности ввиду ограничения ширины спектра полезной составляющей в окрестности нулевой частоты. При увеличении частоты среза погрешность снижается до пренебрежимо малых значений при отсутствии влияния шума. При дальнейшем расширении полосы пропускания фильтра погрешность вновь возрастает за счет влияния составляющей (13). Следует отметить, что увеличение уровня шума приводит к заметному возрастанию погрешности: ввиду уширения спектра полезной составляющей при прореживании (см. рис. 3, г) пропорционально увеличивается остаточная дисперсия шума в полезной полосе частот.
Заключение. Проведенные исследования точности метода асинхронной амплитудной демодуляции сигналов в системах ОКТ, основанного на прореживании последовательности отсчетов сигнала, квадратичном преобразовании информативной составляющей и фильтрации нижних частот, показывают возможность значительного повышения быстродействия получения оценок огибающей сигналов при умеренном возрастании погрешности (до нескольких процентов).
При использовании метода с прореживанием последовательности отсчетов возможен выбор частоты дискретизации, существенно более низкой по сравнению со значением, устанавливаемым критерием Найквиста по отношению к частоте интерференционных полос. Однако при этом следует учитывать эффект уширения спектра огибающей сигнала и, как следствие, возрастание влияния остаточной составляющей дисперсии шума.
Рассмотренный метод прост в реализации и характеризуется высокой помехоустойчивостью по отношению к случайным отклонениям фазы и частоты интерференционных полос при умеренных требованиях к точности задаваемого значения частоты среза фильтра. Метод представляет интерес для использования в системах ОКТ в режиме экспресс-контроля исследуемых объектов при повышенном быстродействии.
список литературы
1. Huang D., Swanson E. A., Lin C. P. et al. Optical coherence tomography // Science. 1991. Vol. 254. Iss. 5035.
P. 1178—1181.
2. Fercher A. F. Optical coherence tomography // J. Biomed. Opt. 1996. Vol. 1. P. 157—173.
3. Schmitt J. M. Optical coherence tomography: a review // IEEE J. Select Topics Quant. Electron. 1999. Vol. 5.
P. 1205—1215.
4. Optical low-coherence reflectometry and tomography / Ed. by B. R. Masters // SPIE Milstone Series. 2001.
Vol. MS165.
5. Handbook of optical coherence tomography / Ed. by B. E. Bouma and G. J. Tearney. NY: Marcel Dekker Inc., 2002.
6. Fercher A. F., Drexler W., Hitzenberger C. K., Lasser T. Optical coherence tomography — principles and applications // Rep. Prog. Phys. 2003. Vol. 66. P. 239—303.
7. Гуров И. П. Оптическая когерентная томография: принципы, проблемы и перспективы // Проблемы когерентной и нелинейной оптики. СПб: СПбГУ ИТМО, 2004. С. 6—30.
8. Tomlins P. H., WangR. K. Theory, developments and applications of optical coherence tomography // J. Phys. D: Appl. Phys. 2005. Vol. 38. P. 2519—2535.
9. Dubois A., Grieve K., Moneron G., Lecaque R., Vabre L., Boccara C. Ultrahigh-resolution full-field optical coherence tomography // Appl. Opt. 2004. Vol. 43. P. 2874—2883.
10. Васильев В. Н., Гуров И. П. Сравнительный анализ методов оптической когерентной томографии // Изв. вузов. Приборостроение. 2007. Т. 50, № 7. С. 30—40.
11. Макс Ж. Методы и техника обработки сигналов при физических измерениях. М.: Мир, 1983. С. 86—89.
12. Васильев В. Н., Гуров И. П., Захаров А. С., Таратин М. А. Высокопроизводительная обработка сигналов с узкополосным спектром на основе метода субдискретизации и нелинейной фильтрации Калмана // Изв. вузов. Приборостроение. 2006. Т. 49, № 8. С. 47—54.
13. Greivenkamp J. E. Sub-Nyquist interferometry // Appl. Opt. 1987. Vol. 26. P. 5245—5258.
14. Gurov I., Taratin M., Zakharov A. High-speed signal evaluation in optical coherence tomography based on sub-Nyquist sampling and Kalman filtering method // AIP Conf. Proc. 2006. Vol. 860. P. 146—150.
15. Alarousu E., Gurov I., Kalinina N., Karpets A., Margariants N., Myllyla R., Prykari T., Vorobeva E. Full-field high-resolving optical coherence tomography system for evaluating paper materials // Advanced Laser Technologies 2007. Proc. SPIE. 2008. Vol. 7022. P. 7022.
16. Цикин И. А. Дискретно-аналоговая обработка сигналов. М.: Радио и связь, 1982. С. 85—95.
17. Бендат Дж., Пирсол А. Измерение и анализ случайных процессов. М.: Мир, 1974. С. 261—264.
18. Васильев В. Н., Гуров И. П. Компьютерная обработка сигналов в приложении к интерферометрическим системам. СПб: БХВ Санкт-Петербург, 1998. С. 83.
19. Ахманов С. А., Дьяков Ю. Е., Чиркин А. С. Введение в статистическую радиофизику и оптику. М.: Наука, 1981. С. 122.
20. Тихонов В. И. Статистическая радиотехника. М.: Радио и связь, 1982. С. 358—360.
Сведения об авторах
Елена Александровна Воробьева — Санкт-Петербургский государственный университет информационных
технологий, механики и оптики, кафедра компьютерной фотоники; младший науч. сотрудник
Игорь Петрович Гуров — д-р техн. наук, профессор; Санкт-Петербургский государственный
университет информационных технологий, механики и оптики, кафедра компьютерной фотоники; E-mail: [email protected]
Максим Владимирович Петерсон — студент; Санкт-Петербургский государственный университет информационных технологий, механики и оптики, кафедра компьютерной фотоники
Рекомендована кафедрой Поступила в редакцию
компьютерной фотоники 25.12.09 г.
