Численное моделирование сдвиговых волн в мягких биологических тканях в среде K-Wave

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

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

Распространение волн в упругих средах подразумевает расчет поля в каждой точке пространства, что означает работу с матрицами. Для этих целей удобно использовать пакет MATLAB, оптимизированный для работы с матрицами большого размера, а для удобства работы с пространством toolbox k-Wave, созданный для расчета распространения волн в упругих средах [1].

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

Численная модель набора скриптов k-wave подразумевает работу с пространством как инициализацию матрицы заданного размера, каждому элементу которой присваиваются необходимые физические параметры, необходимые для расчета поля в точке. В пространстве k-wave данная матрица имеет название kgrid.

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

radius = 6e-3;      % [m]

x_pos = 27.5e-3;    % [m]

y_pos = 45.5e-3;    % [m]

scattering_region = makeBall(Nx_tot, Ny_tot, Nz_tot, round(x_pos/dx), round(y_pos/dx), Nz_tot/2, round(radius/dx));

sound_speed_map(scattering_region == 1) = scattering_c0(scattering_region == 1);

density_map(scattering_region == 1) = scattering_rho0(scattering_region == 1);

Источник излучения в данной модели – точка. При необходимости можно группировать несколько точек в один излучатель.

Рисунок 1. Излучение из двух точечных источников

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

Рисунок 2. Фокусировка излучателей и распространение в плоскости

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

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

На рисунке 3 изображен процесс генерации сдвиговой волны. В отличии от предыдущих, эта симуляция происходит в трехмерном пространстве, и демонстрирует срезы во всех плоскостях, а также разделяет Normal и Shear для наглядности.

Рисунок 3 Генерация сдвиговой волны во времени

Для сравнения численного моделирования с физическим экспериментом использовались результаты, полученные в лаборатории «МедЛаб» Нижегородского госуниверситета с помощью исследовательской системы с открытой архитектурой V-1 Electronics фирмы «Verasonics» [2]. Данная система представляет собой универсальный ультразвуковой диагностический прибор, предназначенный для макетирования и отладки различных алгоритмов медицинской акустики.

Рисунок 4. Фокусировка в системе Verasonics(слева) и расчет в k-wave (справа)

Рисунок 5. Генерация сдвиговой волны в системе Verasonics(слева) и расчет в k-wave(справа)

При проведении исследований на акустической системе Verasonics использовался стандартный многоэлементный датчик L7–4. Рабочая частота была выбрана равной 5 МГц и использовалась стандартная амплитуда диагностических импульсов 50 В.

По результатам численного моделирование получено хорошее сопоставление с физическим моделированием, проведенным на акустической системе Verasonics.