Акустический журнал, 2020, T. 66, № 2, стр. 206-212

ВВЕДЕНИЕ

При ультразвуковых исследованиях часто возникает необходимость формировать изображение области, которая отделена от ультразвуковой решетки набором промежуточных слоев с неизвестными толщинами и скоростями звука в них (рис. 1). Стандартный метод построения изображений основан на расчете времен распространения волн от передающего элемента решетки до визуализируемой точки (x₀, z₀) и от этой точки до приемного элемента с последующим суммированием принимаемых сигналов для всех пар приемо-передающих элементов после компенсации рассчитанных задержек [1–3].

Указанные задержки достаточно просто рассчитать, если решетка располагается непосредственно на границе области визуализации. Однако при наличии промежуточных слоев для каждого элемента решетки и точки фокусировки сначала необходимо найти точки пересечения лучами границ слоев [4–6]. Для одиночного слоя нахождение каждой такой точки пересечения требует решения системы уравнений, а для нескольких слоев возрастающий объем и сложность вычислений приводит к затруднительности практической реализации метода [7, 8]. Упростить решение данной проблемы можно путем разложения сигнала решетки в спектр плоских гармонических волн, однако в описанных алгоритмах учитывается в лучшем случае только один промежуточный слой [9–13]. Следует также отметить, что для реализации данного метода толщины слоев и скорости звука в них должны быть известны с достаточной точностью.

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

ТЕОРЕТИЧЕСКОЕ РАССМОТРЕНИЕ МЕТОДА

Пусть система ультразвуковой визуализации (рис. 1) построена на одномерной решетке элементов, излучение и прием ультразвуковых волн которыми описывается передаточными функциями $H\left( {{{k}_{1}},{\omega }} \right)$ и $H\left( {{{k}_{2}},{\omega }} \right)$ соответственно, где ω – частота, а k₁, k₂ – горизонтальные составляющие волновых векторов излучаемой и принимаемой плоских волн, соответственно. Тогда в линейном приближении выходной сигнал решетки как функция времени t и положений передающего и приемного элементов x₁ и x₂ может быть представлен в виде обратного фурье-преобразования [2, 15]:

где коэффициент $S\left( {{{k}_{1}},{{k}_{2}},{\omega }} \right)$ показывает связь между спектральными компонентами зондирующей и рассеянной волн, определенными в плоскости решетки. Предположим, что в этом сигнале могут быть выделены отклик ${{v}_{r}}\left( { - {{x}_{1}},{{x}_{2}},t} \right),$ отраженный от верхней границы визуализируемой области z = 0, и отклик ${{v}_{{{\text{in}}}}}\left( { - {{x}_{1}},{{x}_{2}},t} \right),$ образованный рассеянием на неоднородностях внутри этой области. При прохождении плоской волны с параметрами (k₁, ω) через набор слоев меняется ее амплитуда, и волна приобретает фазовую задержку где d_m – толщина m-ного слоя, k_zm1 – вертикальная составляющая волнового вектора: и C_m – скорость звука в слое. Предположим также, что k₁ не превосходит критических значений для материалов всех слоев и можно ограничиться рассмотрением распространения только продольных волн. Материалы слоев могут быть анизотропными, в этом случае скорости звука C_m зависят от направления распространения плоской волны. Соотношения, аналогичные (3), могут быть записаны для волны с параметрами $\left( {{{k}_{2}},{\omega }} \right),$ распространяющейся в обратном направлении от объекта к решетке. Тогда спектральный отклик, образованный рассеянием на неоднородностях внутри области визуализации, может быть представлен в виде: где T – амплитудный коэффициент, учитывающий прохождение через границы слоев, а функция рассеяния S₀ неоднородностей в слое отнесена к плоскости z = 0. Для идеального точечного отражателя с координатами $\left( {{{x}_{0}},{{z}_{0}}} \right)$ эта функция учитывает фазовые сдвиги при распространении волн из начала системы координат до этой точки и обратно:

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

где R – амплитудный коэффициент, учитывающий прохождение через слои и отражение от границы, δ – дельта функция.

На основе соотношений (4), (6) можно реализовать следующий метод формирования изображения. По спектру S_r сигнала v_r производится определение фазы произведения ${{H}^{2}}\left( {k,{\omega }} \right)\exp \left\{ {2i{\varphi }\left( k \right)} \right\}$ в полосе временных и пространственных частот, ограниченной передаточной функцией элементов решетки. Затем значения этой фазы для разных комбинаций вычитаются из фазы спектральной плотности S_in сигнала v_in, компенсируя тем самым фазовые искажения, вносимые слоями и элементами решетки. Полученный таким образом спектр может быть использован для формирования изображения в соответствии с известными методами, основанными на коррекции изменений пространственно-временного спектра, возникающими при его распространении вглубь визуализируемой области, и вычислении обратного фурье-преобразования [11–13].

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

где переменные s₁, s₂ имеют смысл проекций векторов медленности падающей и рассеянной плоских волн на ось x, соответственно. Подставляя (1) в (7) и учитывая известное свойство: для спектра откликов плоских волн можно получить выражение:

Учитывая, что элементы решетки являются узкими и обладают широкой диаграммой направленности, их передаточную функцию можно представить в виде произведения временного импульсного отклика p(t) и пространственного множителя H(s). Тогда спектр (9) можно представить в виде свертки по временному параметру τ:

где задает время распространения плоской волны через слои, а величины определяют задержку волн при распространении в области визуализации до отражателя и обратно в горизонтальном и вертикальном направлениях, соответственно. Компонента спектра плоских волн, отраженных от верхней границы, на основе (6) имеет вид: где $s = {{s}_{1}} = {{s}_{2}}.$

Для построения изображения исследуемой области сначала по рассчитанному спектру отклика границы W_r определяют задержку τ_L(s) в зависимости от параметра s. Далее в спектре отклика от внутренних неоднородностей W производят компенсацию задержек τ_L(s₁), τ_L(s₂) для падающей и отраженной плоских волн, соответственно. Для получения изображения точки с координатами (x, z) необходимо также компенсировать задержки τ_x(x), τ_z(z) и провести интегрирование по s₁, s₂:

Если имеется точечный отражатель, находящийся в точке (x₀, z₀), тогда формируемое изображение равно:

где A – амплитудная медленно меняющаяся функция и

Выражение (16) определяет по существу импульсный отклик рассматриваемой системы визуализации. На его основе можно получить оценку размеров импульсного отклика δx, δz, определяющих разрешающую способность в поперечном и продольном направлениях, соответственно.

Пусть временной сигнал p(t) имеет характерный период осцилляций T, а длительность, сравнимую с 2T. Тогда при z = z₀ подынтегральная функция p(τ–Δτ_x) имеет в области интегрирования временной линейный сдвиг, наклон которого пропорционален разности (x₀ – x). Результат интегрирования имеет максимум при x ≈ x₀ и становится малым, если на размере области интегрирования s_m временной сдвиг достигает половины периода T/2:

Значение s_m ограничено максимальным апертурным углом θ_m рассматриваемой системы формирования изображений ${{s}_{m}} \approx C_{L}^{{--1}}{\text{sin}}{{\theta }_{m}}$, откуда можно получить оценку

где λ = C_LT – характерная длина волны. Полученная оценка согласуется с известным значением разрешающей способности [17, 18].

Для оценки протяженности импульсного отклика в продольном направлении δz следует заметить, что в силу ограниченности длительности импульса p(τ) интеграл (16) при x = x₀ становится малым, если временной сдвиг Δτ_z достигает значения T/2. Поэтому вертикальный размер импульсного отклика можно оценить следующим образом:

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

Если точка, где расположен отражатель, имеет горизонтальную координату x₀, а фокусировка производится в точку x, последовательные слагаемые в дискретном представлении выражения (15) имеют относительные временные задержки вида (x₀ – x)Δs, где Δs – шаг дискретизации в спектральной области. Если сдвиг сравним или больше характерного временного периода сигнала T, то при суммировании множества таких откликов не происходит их полного взаимного подавления и имеет место генерация помехи. Поскольку разность координат |x₀ – x| ограничена горизонтальным размером области визуализации 2x_m, ограничение на интервал дискретизации Δs можно представить в виде:

Аналогичным образом для точек, разнесенных по вертикали на расстояние z_m, относительная задержка также должна быть меньше периода:

Используя параксиальное приближение и учитывая, что максимальное значение параметра s не превосходит $C_{L}^{{--1}}$, условие может быть представлено в виде:

Таким образом, чем больше размер области визуализации, тем более мелким должен быть шаг дискретизации спектра Δs.

ЭКСПЕРИМЕНТАЛЬНОЕ ИССЛЕДОВАНИЕ

В эксперименте использовалась плоская 64-х элементная решетка, изготовленная компанией Imasonic [20]. Ультразвуковые элементы решетки имели размеры 0.5 × 12 мм и были расположены с периодом 0.6 мм. Центральная частота и относительная полоса частотной характеристики решетки составляли f₀ = 5 МГц и 70%, соответственно. Возбуждение элементов решетки и прием отраженных сигналов осуществлялись с помощью разработанной 64-канальной электронной схемы, каналы которой могли работать на передачу и прием. Каждый канал был снабжен генератором импульсов с амплитудой 30 В и регулируемой длительностью, а в режиме приема посредством аналогового мультиплексора подключался к аналого-цифровому преобразователю с тактовой частотой 133 МГц и разрядностью 12 бит. Сигналы записывались независимо для всех возможных пар передающих и приемных каналов, характерное время записи полного набора 64 × 64 сигналов составляло 2…4 с. Записанные данные передавались в персональный компьютер для дальнейшей обработки посредством USB интерфейса.

Схема эксперимента представлена на рис. 2. Ультразвуковая решетка 1 была расположена на пластинке 2 из полистирола толщиной 25.4 мм. В качестве тестового образца использовался дюралюминиевый блок 5, отделенный от полистироловой пластинки стойками 4 с высотой 4.16 мм. Эта конструкция помещалась в резервуар 3 с водой, которая заполняла зазор между пластинками. В дюралюминии были выполнены 5 сквозных отверстий, расположенных параллельно поверхности блока и элементам ультразвуковой решетки. Центры отверстий располагались на глубинах 12, 10, 8, 6 и 4 мм с периодом по оси x 3 мм. Диаметр отверстий составлял 0.8 мм, что меньше длины волны ультразвука в дюралюминии 1.2 мм на центральной частоте. Толщина блока составляла 25 мм и существенно превосходила длину ультразвуковых элементов решетки, что позволяло избежать влияния боковых граней блока на распространение волн.

На рис. 3 представлен фрагмент пространственно-временного сигнала v_kj(t), принятого всеми элементами решетки при использовании элемента с номером k = 32 в качестве передающего. Уровень записанного знакопеременного сигнала представлен оттенками серого в соответствии со шкалой, показанной справа. В сигнале выделяется отклик v_r, отраженный от поверхности раздела вода–дюралюминий, и отклики v₁–v₅, отраженные от тестовых отверстий.

Спектр плоских волн W(s₁, s₂, τ), рассчитанный по измеренному пространственно–временного сигналу v_kj(t), показан на рис. 4 для значения s₂ = 0. В спектре наблюдается компонента W_r, даваемая откликом v_r от поверхности блока, и 5 компонент W₁–W₅, образованных отражениями от тестовых отверстий. С целью повышения наглядности представления данных на изображении амплитуда откликов от отверстий увеличена в 50 раз. Диапазон параметров |s₁|, |s₂|, для которых производился расчет спектра плоских волн, был ограничен величиной 0.13 мкс/мм, которое несколько меньше критического значения для алюминия 1/C_L ≈ 0.16 мкс/мм. Считая, что размер области визуализации 2x_m, z_m составляет 20 мм, оценки (22), (24) ограничивают шаг дискретизации величиной Δs = 0.01 мкс/мм. В расчетах он был выбран равным Δs = 0.13/20 = 0.0065 мкс/мм.

Компонента спектра W_r соответствует зеркальному отражению плоских волн от поверхности раздела, поэтому она имеет заметную величину только при s₁ ≈ s₂. На рис. 5 показан этот отклик, по которому согласно (14) была определена зависимость задержки плоских волн в слоях τ_L(s). Используя полученные зависимости, изображение блока с отверстиями было построено в соответствии с выражением (15). На полученном изображении (рис. 6) видно, что вертикальное и горизонтальное положения откликов корректно воспроизводят координаты тестовых отверстий. Оценка ширины протяженности импульсного отклика, произведенная по формулам (20) и (21), дает δx = 0.8 мм и δz = 0.6 мм соответственно, что согласуется с полученными результатами.

ЗАКЛЮЧЕНИЕ

Разработан метод формирования ультразвуковых изображений для конфигурации, в которой между визуализируемой областью и ультразвуковой решеткой располагается ряд слоев с неизвестными толщинами и значениями скоростей звука в них. Метод основан на разложении регистрируемого решеткой полного пространственно-временного сигнала в спектр плоских импульсных волн. На основе построенной теоретической модели показано, что задержки составляющих спектра, приобретаемые при прохождении через дополнительные слои, могут быть измерены по отклику, отраженному от границы области визуализации, и компенсированы при построении изображения. Показано также, что метод обеспечивает пространственную разрешающую способность, отвечающую классическому пределу, а максимальное значение шага дискретизации спектра, ограничиваемое условием отсутствия помех наложения, обратно пропорционально размеру области визуализации.

Для экспериментальной апробации метода использовался алюминиевый блок с тестовыми отверстиями, который исследовался через дополнительные слои из полистирола и воды. Экспериментальные результаты демонстрируют работоспособность метода и корректность полученных теоретических оценок.

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