Известия РАН. Серия физическая, 2021, T. 85, № 6, стр. 823-828

ВВЕДЕНИЕ

Развитие доступного метода измерения сдвиговой упругости скелетных мышц является актуальной и перспективной задачей акустики применительно к медицине [1, 2]. Доплеровская эластография используется для измерения упругости печени и почек [3]. Метод заключается в возбуждении в мягкой ткани низкочастотных сдвиговых волн на ее поверхности и последующей регистрации скорости сдвиговых волн с помощью ультразвуковых зондирующих импульсов. Однако большое затухание сдвиговых волн накладывает ограничения на точность измерений упругости глубоколежащих органов. Диагностика упругости скелетных мышц, расположенных вблизи поверхности, может быть проведена методом доплеровской эластографии достаточно точно, поскольку вблизи поверхности отношение сигнал/шум становится выше, чем в глубоколежащих органах.

ОПРЕДЕЛЕНИЕ УПРУГОСТИ МЯГКОЙ ТКАНИ ДОПЛЕРОВСКИМ МЕТОДОМ

Если в ткани сдвиговая упругость преобладает над сдвиговой вязкостью, то фазовая скорость сдвиговой волны c_t зависит только от модуля сдвига μ и плотности ρ ткани [4]. Тогда, измерив фазовую скорость сдвиговой волны, можно вычислить модуль сдвига мягкой ткани по формуле:

Измерение фазовой скорости сдвиговой волны производится с помощью ультразвуковых зондирующих импульсов, излучаемых пьезопреобразователем. Отражаясь от рассеивателей в мягкой ткани, УЗ зондирующие импульсы приобретают доплеровское смещение частоты, которое содержит в себе информацию о фазе колебаний рассеивателей. Метод обработки зондирующих импульсов, предложенный J. Sato [5] и позднее модифицированный Е.М. Тиманиным [3] с использованием современных электронных компонент, позволяет получить распределение фазы низкочастотных колебаний в мягкой ткани и определить фазовую скорость сдвиговой волны c_t.

Пусть излучаемые УЗ зондирующие импульсы с амплитудой s₀ и частотой ω₀ распространяются вдоль оси z: $s(t) = {{s}_{0}}\sin {{{{\omega }}}_{0}}t.$ Если проекция колебательной скорости рассеивателя ${{\upsilon }_{z}}\left( t \right)$ на ось z на трассе распространения зондирующего импульса меняется по закону:

где υ₀, Ω, и Φ – амплитуда, частота и фаза колебательной скорости рассеивателя соответственно, то отраженный импульс приобретет доплеровский сдвиг частоты: где c – скорость звука в мягкой ткани. Отраженный импульс, принятый УЗ пьезопреобразователем, будет иметь вид: где φ₀ – постоянный набег фазы на трассе распространения УЗ зондирующего импульса, $m = \frac{{2{{{{\omega }}}_{0}}{{\upsilon }_{0}}}}{{c\Omega }}.$ Из (4) следует, что фаза отраженного УЗ импульса определяется колебательной скоростью рассеивателя. Амплитуда и фаза колебательной скорости рассеивателя может быть вычислена путем квадратурной обработки [6] отраженного импульса. Для этого импульс ${{s}_{r}}(t)$ перемножается с прямой и квадратурной компонентами зондирующего импульса. После фильтрации удвоенной частоты 2ω₀ выделяется два квадратурных доплеровских сигнала: синусный Y_s и косинусный Y_c, из которых рассчитывается временной профиль проекции колебательной скорости рассеивателя ${{\upsilon }_{z}}\left( t \right)$ по формуле:

Используя рассчитанные профили ${{\upsilon }_{z}}\left( t \right)$ для двух рассеивателей, расстояние между которыми не превышает длины сдвиговой волны на частоте Ω, можно определить скорость распространения квазиплоской сдвиговой волны в области рассеивателей:

где Δr – расстояние между рассеивателями, α – угол между осью распространения УЗ импульсов и волновым вектором сдвиговой волны, Δψ – разность фаз колебаний рассеивателей. По измеренной фазовой скорости сдвиговой волны по формуле (1) вычисляется модуль сдвига среды.

Длина сдвиговых волн в мягких тканях на исследуемых частотах (150–250 Гц) варьируется в пределах нескольких см. На указанных частотах приповерхностная область мышечной ткани находится в ближнем поле низкочастотного излучателя, имеющего сложную пространственную структуру. Обычно используемый подход [7], основанный на расчете в приближении дальнего поля, может давать значительные погрешности в определении сдвиговой упругости ткани. Для корректного расчета сдвига фаз между рассеивателями необходимо учитывать структуру поля колебательной скорости вблизи низкочастотного излучателя.

РАСЧЕТ И ИЗМЕРЕНИЕ СМЕЩЕНИЙ В БЛИЖНЕМ ПОЛЕ НИЗКОЧАСТОТНОГО ИЗЛУЧАТЕЛЯ

Измерения смещений, создаваемых низкочастотным излучателем, проводились в желатиновом фантоме [8] (№ 1) 1 цилиндрической формы с радиусом 50 мм и высотой 60 мм (рис. 1). Модуль сдвига фантома μ = 11.85 ± 0.18 кПа был измерен статическим методом вдавливания шарика. В фантом был вморожен одноосный миниатюрный акселерометр 2 Brüel & Kjær 4374. Масса акселерометра 0.65 г, диаметр 5 мм, высота 6.3 мм. Эти размеры существенно меньше длины сдвиговой волны (23 мм на частоте 150 Гц), что позволяет проводить локальные измерения поля. Акселерометр располагался на оси цилиндрического фантома на расстоянии 20 мм от нижней поверхности.

В качестве низкочастотного излучателя 3 использовался брусок квадратного сечения со стороной 8 мм и длиной 50 мм, который совершал колебания вдоль оси z под действием вибратора 4 Brüel & Kjær 4810 на частоте 150 Гц. Подвижная часть вибратора, контрольный одноосный акселерометр 5 Brüel & Kjær 8305 и брусок жестко соединялись друг с другом с помощью резьбового соединения. Фантом закреплялся в подвижном держателе на высоте, позволяющей установить под ним низкочастотный излучатель. Акселерометры 2 и 5 подключались к двухканальному осциллографу через усилитель заряда и позволяли измерять амплитуду колебаний излучателя на поверхности фантома и области среды, где находился акселерометр 2. Амплитуда вертикальных колебаний излучателя составляла 0.5 мкм. В начальном положении x₀ излучатель располагался под акселерометром на расстоянии 20 мм от него. В процессе измерений излучатель смещался в горизонтальном направлении x на расстояние 20 мм с шагом 2.5 мм.

Численное моделирование проведено для двумерного случая в прямоугольнике со сторонами 100 × 60 мм (рис. 2а) методом конечных элементов (МКЭ) в программном пакете COMSOL Multiphysics. Модуль сдвига в области прямоугольника μ = 12 кПа, коэффициент Пуассона σ = 0.4999, плотность ρ = 1000 кг/м³. Величина коэффициента динамической вязкости составляла η = 2 Па · с, что характерно для желатиновых фантомов [9]. Границы AB, CD – жесткие стенки; BC, AF, ED – свободные границы. На отрезке FE находится поршневой излучатель бесконечной длины и шириной 8 мм, колеблющийся вдоль оси x с амплитудой 0.5 мкм на частоте 150 Гц. Прямоугольный контур PQRS размерами 5 × 6.3 мм изображает акселерометр. Нижняя сторона контура, как и в лабораторном эксперименте, расположена на высоте 20 мм от границы AD. Для учета конечных размеров акселерометра проводится усреднение поля смещений по нижней стороне PS контура PQRS. При моделировании излучатель смещался вдоль оси x на расстояние 20 мм с шагом 1 мм. В начальном положении x₀ центры излучателя и контура PQRS находились на одной вертикальной прямой. На рис. 2а цветом представлено мгновенное распределение поля вертикальных смещений u_z, рассчитанное для начального положения x₀ излучателя. На расстоянии порядка полутора длин волн сдвиговой волны (23 мм) от излучателя прослеживается регулярная структура поля с периодом повторения 23 мм. Вблизи границ AB, CD и BC она нарушается из-за появления отраженных волн и образования стоячей волны с периодичностью в два раза меньшей. Распределение амплитуды вертикальных смещений акселерометра на высоте 20 мм в зависимости от горизонтального смещения x, рассчитанное с помощью МКЭ, представлено на рис. 2б сплошной линией. Измеренное распределение амплитуды представлено точками, которые для наглядности соединены штриховой линией. Распределения нормированы на собственные максимальные значения. Видно качественное сходство расчетного и экспериментального распределения. Оба распределения имеют локальный минимум вблизи x = 5 мм, локальный максимум вблизи x = 8 мм и схожее увеличение амплитуды при x > 15 мм.

ИЗМЕРЕНИЕ УПРУГОСТИ ФАНТОМА ДОПЛЕРОВСКИМ МЕТОДОМ

В данной серии измерений использовался желатиновый фантом № 2 с модулем сдвига µ = = 14.03 ± 0.24 кПа, в который были внедрены два сферических пластиковых рассеивателя 6 диаметром 1 мм (рис. 1). Нижний рассеиватель был расположен на расстоянии 15 мм от основания фантома, верхний находился на 4 мм выше первого. Низкочастотным излучателем 3 служил полый цилиндр высотой 20 мм с внешним диаметром 20 мм и толщиной стенок 2 мм, что позволяло создавать максимальную амплитуду вертикальных смещений на оси цилиндра. На подвижном держателе над верхним основанием фантома закреплялся ультразвуковой пьезопреобразователь (УЗП) 7. В качестве согласующей среды между фантомом и УЗП использовался слой воды толщиной 4 мм. УЗП имел диаметр 10 мм и резонансную частоту 2.2 МГц. Центр УЗП и рассеиватели располагались на оси цилиндрического излучателя. В качестве источника сигналов для возбуждения УЗП и вибратора использовался генератор сигналов 8 RIGOL DG1062Z с двумя независимыми каналами. Импульсы, отраженные от рассеивателей и принятые УЗП, подавались на цифровой осциллограф 9 АКИП-4111/1, подключенный к компьютеру 10 и синхронизированный с генератором сигналов. Для защиты входа осциллографа во время генерации зондирующего импульса использовался диодный ограничитель амплитуды 11.

Вибратор возбуждал гармонические колебания цилиндрического излучателя вдоль вертикальной оси на частотах в диапазоне 130–160 Гц с амплитудой смещения 0.5 мкм. Возникающие в поле сдвиговой волны колебания сферических рассеивателей регистрировались с помощью зондирующих импульсов с частотой 2.2 МГц и длительностью 3.2 мкс. Период повторения импульсов 0.7 мс. Такие параметры зондирующих импульсов позволяли проводить измерения на глубине до 54 см с пространственным разрешением 2.4 мм. Импульсы, принятые УЗП во временном окне длительностью 70 мкс оцифровывались с частотой дискретизации 62.5 МГц и записывались в память компьютера. За одно измерение производилась запись 500 таких реализаций, что было ограничено памятью осциллографа. Время записи составляло 0.35 с, что позволяло записать от 45 до 56 периодов колебаний рассеивателей. Квадратурная обработка отраженных импульсов производилась в программе, разработанной в среде LabView. Алгоритм программы был проверен на модельных сигналах. Моделировались импульсы, отраженные от рассеивателей, колеблющихся под действием сдвиговой волны с частотой 100 Гц в среде с модулем сдвига 10 кПа. Положение рассеивателей задавалось в соответствие с рис. 1. Обрабатывались реализации общей длительностью 0.1 с, оцифрованные с частотой 25 МГц. Значение модуля сдвига, полученное из модельных сигналов, совпало с заданным значением 10 кПа с погрешностью не хуже 0.1%.

На рис. 3 символами представлены результаты измерения модуля сдвига желатинового фантома № 2 в диапазоне частот 130–160 Гц. Значение модуля сдвига 14.03 кПа, полученное из статических измерений показано сплошной линией. Значения μ, измеренные на частотах 130, 135, 155 и 160 Гц соответствуют статическому модулю в пределах погрешности измерений, не превышающей 5%. Однако на трех частотах значения модуля сдвига имеют существенный разброс. Скорее всего, такие особенности связаны со структурой ближнего поля и наличием переотражений сдвиговой волны в области рассеивателей.

Чтобы выяснить влияние структуры ближнего поля на измерения, было проведено моделирование экспериментальных условий, но с использованием более простого излучателя в виде бруска вместо полого цилиндра. Для того, чтобы минимизировать влияние переотражений, моделирование провели для частоты 250 Гц. Направление волнового вектора плоской сдвиговой волны получим как перпендикуляр к волновому фронту, проведенному из центра излучателя. На рис. 4 приведена зависимость фазы волны от расстояния l, отсчитываемого вдоль отрезка, совпадающего по направлению с построенным волновым вектором. Угол наклона отрезка к вертикальной оси α = 40°. Видно, что фаза меняется нелинейно, особенно вблизи границы, где имеются переотраженные волны, а также в начале отрезка, вблизи излучателя, где имеются волны утечки, связанные с распространением поверхностных волн. В диапазоне расстояний l = 15–40 мм, отсчитанных вдоль волнового вектора, фаза меняется линейно, поэтому в этой области можно проводить вычисления с использованием приближения плоских волн. Таким образом, можно выбирать оптимальную область измерений с учетом геометрии излучателя и частоты низкочастотных вибраций.

ЗАКЛЮЧЕНИЕ

Рассмотрены особенности применения УЗ доплеровского метода для измерения упругости скелетных мышц, залегающих на небольшой глубине (2–4 см) от поверхности кожи. Разработан и реализован алгоритм вычисления сдвигового модуля по измерению фазового сдвига двух рассеивателей, двигающихся в поле сдвиговой волны. Экспериментальная реализация достаточно проста и выполнена с использованием стандартной измерительной аппаратуры. Показано, что в отличие от измерений упругости глубоколежащих органов, измерения упругости мышц можно проводить за доли секунды при времени накопления, не превышающем 0.3 с, что существенно при проведении диагностических процедур. По нашим оценкам, при диагностике мышц, расположенных на глубине до 3 см, полезный сигнал может превышать в 2 раза уровень шума при колебаниях вибратора с амплитудой ускорения порядка 0.5 м/с², что сопоставимо с амплитудой обычных вибромассажеров. Небольшая глубина зондирования мышц (3 см) позволяет использовать вибрации с частотами 200–300 Гц, что обеспечивает пространственное разрешение предлагаемого метода на уровне 1–2 мм. Показано, что существуют определенные трудности в получении высокой степени точности измерений, обусловленные сложной структурой ближнего поля вибратора и отличием фазы сдвиговой волны в области рассеивателей от волны с плоским фронтом. Для корректного учета структуры ближнего поля предлагается проводить вычисления его структуры для конкретных геометрий вибраторов и корректировать вычисления фазовой скорости сдвиговых волн. В настоящее время ведутся работы по реализации конкретных схем расчетов для планируемых вибраторов.

Исследование выполнено при поддержке Российского научного фонда (проект № 19-72-00086) и Фонда развития теоретической физики “Базис”.