АКУСТИКА ЖИВЫХ СИСТЕМ. ^^^^^^^^^^^^ БИОМЕДИЦИНСКАЯ АКУСТИКА
УДК 534.2:534.7
ВОССТАНОВЛЕНИЕ ПРОСТРАНСТВЕННЫХ РАСПРЕДЕЛЕНИЙ СКОРОСТИ ЗВУКА И ПОГЛОЩЕНИЯ В МЯГКИХ БИОТКАНЯХ ПО МОДЕЛЬНЫМ ДАННЫМ УЛЬТРАЗВУКОВОГО ТОМОГРАФИРОВАНИЯ
© 2014 г. В. А. Буров, Д. И. Зотов, О. Д. Румянцева
Московский государственный университет им. М.В. Ломоносова Физический факультет, 119991, ГСП-1, Москва, Ленинские горы Тел.: (495) 939-3081; факс: (495) 932-8820 E-mail: burov@phys.msu.ru Поступила в редакцию 20.01.2014 г.
Для восстановления пространственных распределений акустических характеристик мягких биологических тканей — скорости звука и коэффициента поглощения — используется двухшаговый алгоритм. Знание этих распределений актуально для целей ранней диагностики доброкачественных и злокачественных новообразований в биотканях, прежде всего, — в молочной железе. На первом шаге оцениваются крупномасштабные распределения, которые на втором шаге уточняются с высоким разрешением. Представляются результаты восстановления на основе модельных исходных данных. Иллюстрируется принципиальная необходимость предварительного восстановления крупномасштабных распределений и последующего учета их на втором шаге. Привлечение для обработки технологии CUDA позволяет получать итоговые изображения форматом 1024 х 1024 отсчетов всего за несколько минут.
Ключевые слова: акустическая томография, пространственные распределения скорости звука и коэффициента поглощения, крупномасштабные неоднородности, тонкая структура.
DOI: 10.7868/S0320791914040029
ВВЕДЕНИЕ
Работа направлена на развитие методов ультразвуковой томографии, предназначенных для восстановления пространственных распределений акустических характеристик мягких биологических тканей — скорости звука и коэффициента поглощения в них. Знание этих характеристик актуально для целей ранней диагностики доброкачественных и злокачественных новообразований в биотканях и, прежде всего, — в молочной железе. Для получения экспериментальных данных от исследуемого объекта, т.е. рассеивателя (биоткани, в медицинских приложениях), используется двумерная круговая томографическая схема, в которой приемоизлучающие преобразователи располагаются на окружности, охватывающей томографируемый объект со всех сторон.
Поскольку в общем случае объект сильно искажает зондирующее его поле, то при обработке экспериментальных данных необходимо принимать во внимание процессы многократного рассеяния волн внутри рассеивателя. Как следствие, алгоритмы, основанные на лучевом [1], первом борновском [2—8] или рытовском [2—6] приближениях, имеют весьма ограниченную область применимости в прикладных системах. Нелинейность зависимости экспериментальных данных (рассеянного поля) от характеристик рассеивателя, т.е.
нелинейность соотношений, описывающих процессы многократного рассеяния, означает необходимость использования более сложных процедур обработки данных. Наибольшей относительной простотой обладает рассматриваемый ниже двухшаговый алгоритм, позволяющий в ряде случаев получить приемлемые результаты после двух этапов, существенно различных по своему характеру. Строго говоря, двухшаговый алгоритм не является итерационным, так как он с самого начала предполагает выполнение всего одного цикла, в то время как итерационные процедуры предполагают неопределенное количество циклов постепенного приближения к решению. Своим возникновением двухшаговый подход обязан специфике задач медицинской томографии, при которой неоднородности в фазовой скорости звука и большие размеры областей с такими неоднородностя-ми вызывают сильное искажение волновых фронтов в ткани. Это, в свою очередь, приводит к невозможности прямого применения линеаризованного борновского приближения. Тем не менее, двухшаговый подход, используя пространственные особенности проблемы и за счет этого как бы предельно упрощая многократное повторение итерационного цикла, с одной стороны, содержит много общего с "классическими" итерационными методами, а с другой стороны, ис-
пользует на своем втором шаге линеаризованное приближение, однако на неоднородном фоне.
Предположим, что рассматривается достаточно типичная ситуация, когда нужно восстановить рассеиватель, состоящий из мелких характерных деталей объекта (тонкой структуры), которые присутствуют на фоне крупномасштабных неизвестных неоднородностей фазовой скорости и коэффициента поглощения. Мелкие детали имеют линейный размер от нескольких десятых долей длины волны, лежащей в миллиметровом диапазоне, до нескольких длин волн; размер крупномасштабных неоднородностей составляет несколько длин волн и более. Полная оценка всего рассеива-теля, содержащего как крупные фрагменты, так и мелкие детали с "борновским" контрастом, может быть осуществлена на основе идеи использования борновского алгоритма, но с учетом неоднородности фона [9, 10], предварительно оцениваемого на основе экспериментальных данных. В [11] работоспособность двухшагового алгоритма была подтверждена численным моделированием. Обсуждаемые ниже численные примеры по восстановлению крупномасштабных неоднород-ностей и тонкой структуры объекта служат дальнейшей иллюстрацией возможностей двухшагового алгоритма. Используемая двумерная томографическая схема соответствует конструкции опытного образца линейного ультразвукового медицинского томографа [12], действующего по принципу линейного (т.е. без изменения частоты) рассеяния волн. Квазиточечные приемоизлучаю-щие преобразователи расположены на антенном кольце радиуса Я0; томографируемый объект находится внутри этого кольца. Каждый из преобразователей может выступать в роли как излучателя, так и приемника.
ДВУХШАГОВЫЙ АЛГОРИТМ
ВОССТАНОВЛЕНИЯ РАСПРЕДЕЛЕНИЙ СКОРОСТИ ЗВУКА И ПОГЛОЩЕНИЯ
На первом шаге предполагается восстановление распределения по объекту крупномасштабных неоднородностей — областей со значительными отклонениями по скорости звука с(г) или амплитудному коэффициенту поглощения а( г, ю) от постоянных значений с0 и а0, которые предварительно определяются в фоновой среде без объекта, например, по методике [13]. Для восстановления используется процедура, которая далее условно будет называться "квазирентгеновским" алгоритмом [11, 14, 15], поскольку восстановление осуществляется в приближении прямолинейного распространения акустических сигналов через объект, т.е. эффекты рефракции не учитываются. Ниже кратко упоминаются только основные моменты этой процедуры.
Экспериментальные данные Э К8 измеряются для различных пар "излучатель (индекс ^—при-
емник (индекс К)"; геометрические положения излучателя и приемника в полярной системе координат с началом в центре антенного кольца задаются
радиус-векторами г5 = {я0,ф^} и гя = {я0,фЯ} соответственно. Восстановление крупномасштабных распределений медленности 1/ с0(г) и амплитудного коэффициента поглощения а 0(г, ю) происходит по одному и тому же алгоритму времяпролетного типа на основе экспериментальных данных Э . Для оценки пространственного распределения медленности полагается Э = tRS - tЯS, где tRS и ! — времена распространения сигналов от излучателя S до приемника Я в присутствии объекта и без него соответственно. При этом в каждой фиксированной
точке г восстанавливается функция 3(г) = —----.
С0(г) С0
Для оценки пространственного распределения поглощения полагается Э ^ = - ), где
Л^, А^ — эффективные амплитуды сигналов, приходящих на приемник К от излучателя S, в присутствии объекта и без него. Значения Эхарактеризуют степень дополнительного затухания волны относительно ослабления за счет расходимости волнового фронта; в этом случае 3(г) = а0(г,ю) - а0(ю). Если сигнал при обработке не раскладывается на монохроматические составляющие, то вместо частотно-зависимых коэффициентов поглощения а 0(г, ю) и а 0(ю) будут фигурировать их значения, усредненные по ю в рабочем частотном диапазоне с шириной Дю, т.е.
-—-Аш -Аш -7Га® -А®
а0(г) и а0 ; тогда ^(г) = а0(г) - а0 .
Восстановление распределения 3(г) в каждой фиксированной точке г = {х, у} (в декартовых координатах) по известным значениям Э = Э^, ф!) осуществляется из соотношения
2п
3(г) = [Я(г, фS) 4п2 •>
1
где
g(г, Фs) = Иш
81,82 —>+0
тт2/ 0Ч
н (г, фs)
ё ф S,
(1)
А 2
' р?1 -25- 2п ^
I +.
0 р"с1+252
81и
0
' Фя
0
Фs
ео8
0
'фя
0
хЭ^, ФЯ )ёфЯ
+ ! + ^ | (ео8 Y)Э(фS, фП
Н 2(г, ) = Я2
-у
X 008 ф5 + у 81П
0\ Ф5 ).
г = {х, у};
У = у(г, ) = агс1ё
X 81П ф5 - у 008 ф5
У е
п п
I. 2,2J
ехс1
; Фя
Я0 - (х 008 ф^ + у 81п ф^ ) (ф£ + 2у + я) , фЯхс1 е [0,2п).
Смысл величин, входящих в (1), (2), заключается в следующем. Величина Н — расстояние между излучателем Б и точкой г; у — угол между направлениями из точки Б на начало координат О и
на точку г. Особый угол фяс1 соответствует такому
лехс1
положению приемника Я , при котором хорда,
о т>ехс1
соединяющая точки Б и Я , проходит через точку г. При расчете функции ^(г, ф^.), согласно (2), выделяется малая окрестность (не обязательно
\ Г ехс1 ехс1 по "1 г
симметричная) I фЯ - 25^ фЯ + 2о2 I особого угла, поскольку подынтегральное выражение в точ-
0 ехс1 г г^
ке ф Я = ф Я имеет особенность. Однако результат вычислений функции g(г, ) не зависит от выбираемых малых значений 5! > 0 и 52 > 0.
Смысл процедуры второго шага описан в [11, 16]. В основе лежит известный принцип восстановления слабых рассеивателей (т.е. рассеивате-лей, слабо искажающих падающее на них поле), однако на уже известном неоднородном фоне. Такой фон задается найденными из первого шага крупномасштабными распределениями С0(г) и
-Аш
а0(г) . Ниже будут отмечены основные моменты этой процедуры и особенности ее технической реализации в упомянутом опытном образце ультразвукового томографа.
В двумерном случае падающее на рассеиватель поле и0(г, г5,0 представляет собой цилиндрическую волну с амплитудным спектром А(ю). Это падающее поле при |г5 - г| > 1.5^ шах (Кшах — максимальная рабочая длина волны) имеет вид
«0(г, г5,0 =
1
| А(ю)ехрI -Ш + I—|г - г5| Iйю, г е %
V (г) =
Г С2 С0
2
С0
С0 (г) С (г),
| Б(ю, ¿0)йю + с^Тр^г) х
X V
21 1 ] гБ(ю,
л/Р(г)
7)) V
2
йю - 2с0/ х
(3)
ю
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.