High resolution acoustic tomography based on backpropagation of waves

Full Text

Введение

Акустическая томография используется в медицине для диагностики и визуализации биологических органов [1–2], в неразрушающем контроле для обнаружения дефектов различных сплавов и структур [3–5], при изучении океана [6–7]. Методы акустической томографии также широко применяются в дефектоскопии при решении обратных задач – по измеренному сигналу, рассеянному неоднородностью, для определения ее формы и характеристик [8–11]. Обратные задачи не всегда имеют однозначное решение и крайне чувствительны к шумам измерений и несоответствиям математической модели и условий зондирования.

Метод обратного распространения волн в настоящее время все чаще используется при решении задач мониторинга неоднородных сред и в активной локации [12–13]. Этот метод позволяет учитывать фоновые неоднородности среды для визуализации источника. Также этот метод в литературе носит название The Time Reversal technique (TR). Метод TR можно использовать для обнаружения рассеивателей. Для достаточно малых рассеивателей можно использовать тот факт, что рассеиватели ведут себя как вторичные источники. В статье [14] в контексте неразрушающего контроля конструкций обсуждается вычислительное использование TR для визуализации точечных рассеивателей на основе обратного распространения поля из области измерений. Поскольку анализ зависит от времени, необходимо применять критерии остановки, чтобы остановить обратное моделирование, когда волны оптимально сфокусированы в месте расположения рассеивателя. Рассмотрены два критерия остановки, оба взяты из области обработки изображений: один основан на энтропии Шеннона, а другой – на норме ограниченной вариации. К аналогичному выводу пришли авторы работы [15], имеющие дело только с точечными рассеивателями, а предложенные ими алгоритмы остановки вообще не способны идентифицировать рассеиватели конечных размеров. В работе [16] авторы показали, что визуализация TR успешна только для эффективных точечных источников (или рассеивателей, которые являются вторичными источниками). Чтобы идентифицировать и охарактеризовать более крупные рассеиватели, в сочетании с TR используются дополнительные вычислительные инструменты.

Большинство методов томографии обратного рассеяния не учитывают фоновые препятствия и неоднородности среды [17–18]. Методы, учитывающие неоднородность среды, работают для частных случаев, например для плоскослоистых сред [19]. Для произвольных препятствий могут применяться приближенные методы расчета фазовых сдвигов волнового поля в неоднородных средах [20].

Наиболее перспективно развитие метода обратной временной миграции (Reverse Time Migration, RTM) [21–22]. Этот метод широко используется в сейсмологии и позволяет учитывать заранее известные фоновые преломляющие препятствия. Существующие методы RTM основаны на использовании короткоимпульсного зондирования при условии высокой локализации зондирующей волны в пространстве.

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

1. Схема измерений

Предлагается рассмотреть следующую схему измерений (рис. 1). Исследуемый объект погружается в емкость с водой на определенную глубину. На боковой стенке емкости с водой размещен ультразвуковой излучатель, который облучает исследуемый объект ультразвуковыми волнами. Волновое поле, рассеянное на объекте, регистрируется приемником в воде, размещенным на двухкоординатном сканере, перемещающемся над объектом контроля. Будем считать, что зондирующий сигнал является широкополосным. Сигнал может быть короткоимпульсным, шумоподобным, или любой другой формы, но иметь широкую полосу частот.

 

Рис. 1. Схема измерений

 

2. Метод решения обратной задачи

Задача состоит в том, чтобы визуализировать рассеивающий объект в среде по измеренному полю, прошедшему через фоновые преломляющие неоднородности. Предлагается применить метод обратной временной миграции (Reverse Time Migration – RTM), включающий обратное распространение волн (метод обратного времени) в сочетании с дополнительной фильтрацией поля прямой волны от излучателя. Для зондирования короткими импульсами метод RТМ предполагает, что произведение акустического давления поля волны обратного распространения на поле прямой волны позволит визуализировать рассеиватели в среде. Для случая произвольного зондирующего сигнала предлагается включить в алгоритм RTM дополнительную фильтрацию во временной области.

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

$\begin{array}{l}\tilde{P}\left(r\right)=\underset{x_m}{\int }\underset{y_m}{\int }\tilde{U}(x_m,y_m,\omega )\times \\ \times \frac{z}{2\pi r^2}i{k}_0\exp \left(i{k}_{0}r\right)d{x}_{m}d{y}_m,\end{array}$ (1)

где $r=\sqrt{{(x-x_m)}^2+{(y-y_m)}^2+z^2};$ $\mathbf{r}=(x,y,z)$ – точка наблюдения дифракционного поля; $(x_m,y_m)$ – эквивалентный точечный источник в плоскости $z=0$.

По формуле (1) можно восстановить дифракционное поле за экраном с отверстием, если измерить поле в отверстии $\tilde{U}(x_m,y_m,\omega ).$ В методе RTM необходимо изменить фронт волны на экране (зоне измерения), обратив время в обратную сторону. Рассмотрим обратное распространение поля от отверстия в экране (зона измерения) для восстановления поля источников в среде на основе формулы (1). Для этого предлагается инвертировать время, что эквивалентно замене мнимой единицы в (1) с $+i$ на $-i$. В результате после перехода во временную область получаем поле обратных волн:

$P\left(\mathbf{r},t\right)=\underset{-\infty }{\overset{\infty }{\int }}\underset{x_m}{\int }\underset{y_m}{\int }\tilde{U}(x_m,y_m,\omega )\frac{z}{2\pi r^2}\left(-i{k}_0\right)\times$ (2)

$\times \exp \left(-i{k}_{0}r\right)d{x}_{m}d{y}_m{e}^{-i\omega t}d\omega =$

$=\frac{z}{2\pi }\underset{x_m}{\int }\underset{y_m}{\int }\frac{1}{r^{2}c}\frac{\partial }{\partial t}U\left(x_m,y_m,\frac{r}{c}+t\right)d{x}_{m}d{y}_m.$

Предлагается метод RTM высокого разрешения на основе дополнительной фильтрации поля прямого распространения. Для случая произвольного зондирующего сигнала $S\left(t\right)$ линейную фильтрацию можно осуществлять следующим образом:

$B\left(\mathbf{r}\right)=\underset{t=0}{\overset{T}{\int }}P(\mathbf{r},t)\left|\mathbf{r}-{\mathbf{r}}_0\right|\cdot W\left(t-\frac{\left|\mathbf{r}-{\mathbf{r}}_0\right|}{c}\right)dt,$ (3)

где умножение на $\left|\mathbf{r}-{\mathbf{r}}_0\right|$ вводится для компенсации убывания поля зондирующего излучателя;$W\left(t\right)$  – импульсная реакция линейного фильтра, вычисляемая на основе обратной свертки с регуляризацией из зондирующего сигнала:

$W\left(t\right)=\underset{-\infty }{\overset{\infty }{\int }}\frac{{\tilde{S}}^{*}\left(\omega \right)}{\tilde{S}\left(\omega \right){\tilde{S}}^{*}\left(\omega \right)+\alpha }\exp \left(i\omega t\right)d\omega ,$ (4)

где $\tilde{S}\left(\omega \right)$ – спектр зондирующего сигнала; ${\tilde{S}}^{*}\left(\omega \right)$ – комплексно сопряженное значение $\tilde{S}\left(\omega \right)$; $\alpha$ – параметр регуляризации, значение которого подбирается для каждого типа зондирующих сигналов и определенного уровня шума. Желательно выбирать $\alpha >0$ в зависимости от дисперсии шума $\alpha ={\sigma }^2/\Delta \omega ,$ где ${\sigma }^2$ – дисперсия шума в измеряемом сигнале, $\Delta \omega$ – ширина полосы частот. Если ${\sigma }^2=0$, то в любом случае $\alpha$ следует выбирать больше 0, чтобы избежать деления на 0 в (4).

Функция $B\left(\mathbf{r}\right)$ является решением обратной задачи в приближении однократного рассеяния. По сути, (3) и (4) реализуют обратную свертку с регуляризацией сигнала и обеспечивают преобразование зондирующего сигнала в эквивалентный короткий дельта-подобный импульс. Классическая реализация RTM по существу выполняет согласованную фильтрацию сигнала,
которая максимизирует уровень интенсивности восстановленного изображения, но не обеспечивает наилучшего разрешения. В предлагаемом методе спектр зондирующего сигнала выравнивается по всем частотам, что позволяет значительно повысить разрешающую способность восстанавливаемого изображения. Таким образом, предлагаемый метод дает возможность повысить разрешающую способность по сравнению с классическим методом RTM.

3. Решение обратной задачи в неоднородной фоновой среде

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

Рассмотрим метод численного моделирования акустических процессов (решение прямой задачи), основанный на волновом уравнении для жидкой неоднородной среды:

$\Delta p-\frac{1}{c^2}\frac{{\partial }^{2}p}{\partial t^2}=S,$ (5)

где $p(x,y,z,t)$ – давление в среде; $c(x,y,z)$ – пространственное распределение скорости звука в среде; $S(x,y,z,t)$ – функция, описывающая источники. Численное решение этого уравнения осуществляется методом конечных разностей во временной области.

Предлагается восстанавливать волновое поле в среде (для решения обратной задачи) методом обратного распространения волн (техникой обращения времени). Этот метод инвариантен относительно замены времени, $t$ на $-t$ что позволяет сфокусировать волновое поле в окрестности локализованного источника колебаний, «переизлучая» принятые сигналы в обратном направлении. Согласно методу обратного распространения волн, в численной модели в области измерения задается множество источников, каждый из которых излучает сигнал, принятый в данной точке, но в обратном направлении. При распространении волн от таких источников они складываются в фазе в тех точках среды, где в определенный момент времени находятся рассеиватели. В области расположения источника сигнала волна обратного распространения сходится, что позволяет визуализировать изотопный источник в среде. Для проверки возможности восстановления поля обратного распространения было проведено численное моделирование. В качестве фоновой среды исследовалась вода (скорость звука 1440 м/с, плотность 1000 кг/м³). В качестве рассеивающих объектов рассматривались 4 точечных источника расположенных на разном расстоянии и высоте друг от друга (скорость звука 5500 м/с, плотность
4000 кг/м³). На верхней границе расчетной области располагались приемники, которые регистрировали дошедший до них сигнал. Рассеивающие объекты зондировались сферическим волновым пучком. Источник волнового пучка располагался сбоку и излучал волну в виде короткого импульса. При численном моделировании пространственный шаг составлял 0,25 мм, шаг по времени – 10 нс. В градациях серого на рис. 2, а показано смоделированное акустическое давление волны прямого распространения в момент времени после того, как произошло рассеяние. Прямоугольником обозначено препятствие плотностью 500 кг/м³ и скоростью звука 2880 м/с. Препятствие было введено для демонстрации возможности учета его влияния предлагаемым методом. На рис. 2, б показано поле обратных волн с учетом прямоугольного препятствия в момент прохождения волнами места расположения рассеивателей.

 

Рис. 2. Численное моделирование зондирования объектов через прямоугольную неоднородность (а – прямое распространение волны от источника и рассеянное на исследуемых объектах; б – обратное распространение поля от рассеивателей)

 

Трудность возникает при выборе момента времени, когда максимум поля обратного распространения будет точно сходиться в точке, где находятся рассеивающие объекты. Выбор этого момента времени связан с задержкой распространения прямого сигнала от зондирующего излучателя. Для восстановления томографического изображения рассеивающих объектов предлагается учитывать преобразованное по формуле (4) поле прямой волны от зондирующего излучателя. Первоначально считаем, что среда не содержит рассеивателей, и моделируем в ней распространение волны от источника. Далее проводится моделирование обратного распространения волн от приемников с обращением времени. Поле от обратной волны умножается на модифицированную прямую волну и интегрируется по времени, как в (3), для получения изображения рассеивателей. В области расположения рассеивателей прямая и обратная волны перекрываются, что при их наложении (путем умножения) позволяет визуализировать рассеивающий объект. На рис. 3, а показано восстановленное изображение рассеивающего объекта в виде точечных объектов по результатам численного моделирования методом RТМ. На рис. 3, б отражен результат восстановления изображения предложенным методом высокого разрешения.

 

Рис. 3. Восстановленное томографическое изображение рассеивающих объектов (а – результат, полученный методом RТМ; б – результат, полученный предложенным методом высокого разрешения)

 

Видно, что рассеивающие объекты визуализируются как методом RTM, так и предложенным методом высокого разрешения. Однако визуализированное изображение рассеивающих объектов предлагаемым методом имеет более высокое разрешение (рис. 3, б).

4. Экспериментальная проверка метода

Для экспериментальной проверки предложенного метода восстановления томографического изображения рассеивающего объекта в воде была собрана экспериментальная установка, состоящая из двухкоординатного сканера с установленным на нем ультразвуковым приемным датчиком. В емкость с водой помещались излучатель ультразвуковых сигналов с центральной частотой 200 кГц и рассеивающий объект в виде изогнутой синусоидальной поверхности из PLA пластика. Ультразвуковой приемник погружался в воду и перемещался по плоскости сканирования посредством сканера. В каждой точке пространственного сканирования с шагом 1,875 мм сигнал в приемном датчике измерялся в течение 10 мс. Сканирование осуществлялось сигналами с линейной частотной модуляцией в диапазоне частот от 20 до 350 кГц. Фотография исследуемого объекта представлена на рис. 4, а, а фотография экспериментальной установки – на рис. 4, б.

 

Рис. 4. Фотография тестового объекта (а) и экспериментальной установки (б)

 

В качестве источника ультразвуковых сигналов в воду был помещен пьезокерамический излучатель диаметром 10 мм. Форма зондирующего сигнала показана на рис. 5.

 

Рис. 5. Форма зондирующего сигнала (широкополосный сигнал в диапазоне частот от 20–350 кГц)

 

Приемник, перемещаясь в области 240 × 240 мм, регистрировал рассеянное от исследуемого объекта поле. Измеренные данные вносились в разработанную численную модель. При комбинировании решения с положительным и отрицательным ходом времени (3) было восстановлено томографическое изображение рассеивающей неоднородности. На рис. 6 показан результат восстановления томографического изображения рассеивающего объекта по экспериментальным данным.

 

Рис. 6. Восстановленное томографическое изображение рассеивающего объекта в виде криволинейной синусоидальной поверхности

 

На изображении визуализируется криволинейная синусоидальная поверхность. Искажение изображения вызвано неточным согласованием зондирующего сигнала с математически смоделированным. Тем не менее экспериментальные результаты подтверждают принципиальную применимость предложенного метода.

Заключение

Предложен метод акустической томографии высокого разрешения на основе обратной временной миграции, дополненного специальной предварительной фильтрацией сигнала поля прямого распространения. Метод учитывает влияние фоновых преломляющих неоднородностей среды. Поскольку рассматривается зондирование среды широкополосными сигналами, возможно получение трехмерных изображений рассеивающих объектов. Преимуществом предлагаемого метода является повышение разрешающей способности относительно метода RTM. Проведены экспериментальные исследования по визуализации объекта, помещенного в воду, показавшие применимость предлагаемого метода. Зондирование проводилось в диапазоне частот от 20 до 350 кГц. Предложенный метод улучшает разрешение томографических изображений и позволяет учитывать влияние препятствий.

Исследование выполнено при финансовой поддержке РФФИ в рамках проекта 20-32-90074.

About the authors

Dmitry Ya. Sukhanov

National Research Tomsk State University

Email: sdy@mail.tsu.ru
ORCID iD: 0000-0002-0805-4543

Doctor of Physical and Mathematical Sciences, professor of the Department of Radiophysics, National Research Tomsk State University, Tomsk, Russia.

Research interests: wave tomography, acoustics, digital signal processing.

Russian Federation, 36, Lenin Ave., Tomsk, 634050, Russia

Anzhela E. Kuzovova

National Research Tomsk State University

Author for correspondence.
Email: ang_kuz93@mail.ru

postgraduate student of the Department of Radiophysics, National Research Tomsk State University, Tomsk, Russia.

Research interests: acoustic tomography, acoustics, image restoration.

Russian Federation, 36, Lenin Ave., Tomsk, 634050, Russia

References

Supplementary files