Журнал вычислительной математики и математической физики, 2022, T. 62, № 5, стр. 854-860

Решение двумерной обратной задачи квазистатической эластографии с помощью метода малого параметра

А. С. Леонов^1, *, Н. Н. Нефедов^2, **, А. Н. Шаров^2, ***, А. Г. Ягола^2, ****

¹ НИЯУ “МИФИ”
115409 Москва, Каширское ш., 31, Россия

² МГУ, физический факультет, кафедра математики
119992 Москва, Ленинские горы, Россия

^* E-mail: asleonov@mephi.ru
^** E-mail: nefedov@phys.msu.ru
^*** E-mail: scharov.aleksandr@physics.msu.ru
^**** E-mail: yagola@physics.msu.ru

Поступила в редакцию 31.10.2021
После доработки 31.10.2021
Принята к публикации 14.01.2022

EDN: LGWCWE
DOI: 10.31857/S0044466922050076

Полный текст (PDF)

Аннотация

Рассматривается прямая и обратная задачи двумерной квазистатической эластографии в рамках модели деформации исследуемой биологической ткани как упругого тела, подвергаемого поверхностному сжатию. Возникающие смещения ткани в приближении плоского линейного деформированного состояния упругого тела описываются краевой задачей для уравнений в частных производных с коэффициентами, которые определяются модулем Юнга и постоянным коэффициентом Пуассона ткани. Эта краевая задача содержит малый параметр, что позволяет решить ее с помощью теории регулярных возмущений уравнений в частных производных. Исследуется процедура такого решения и при некоторых предположениях выписываются простые формулы для решения как прямой, так и обратной задач двумерной квазистатической эластографии. Представлены результаты численных экспериментов по решению прямой и обратной модельных задач по предлагаемым формулам. Полученные результаты достаточно хорошо отражают модельные решения. Расчеты по формулам требуют долей микросекунд на персональном компьютере средней производительности для достаточно мелких сеток, так что предлагаемый подход с использованием малого параметра может быть применен при онкологической диагностике в реальном времени. Библ. 16. Фиг. 2.

Ключевые слова: двумерная квазистатическая эластография, обратные задачи, метод малого параметра, регуляризация.

1. ВВЕДЕНИЕ

В последние десятилетия бурно развивается отрасль онкологической диагностики, называемая эластографией (см. [1]–[5]). Она объединяет ряд специфических методов, основанных на различиях в механических характеристиках здоровой и опухолевой ткани определенных типов. Все эти методы структурно состоят из следующих этапов: воздействие на исследуемую часть тела поверхностными силами; измерение (или вычисление) возникающих деформаций исследуемой биологической ткани; определение характеристик упругости ткани путем решения обратной задачи. Деформации (смещения) тканей определяются по данным ультразвуковых исследований или с помощью магнитно-резонансного метода. Принципиальным моментом в такой схеме является решение обратной математической задачи в рамках некоторой модели биологической ткани.

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

В теории эластографии большое распространение получила модель, основанная на уравнениях квазистатической теории упругости (см. [5]–[7]). В ней участок исследуемой ткани, характеризуемый распределениями механических модулей, представляется как линейно-упругое изотропное тело, подвергаемое малым поверхностным сжатиям. Решение обратной задачи для такой модели позволяет найти по известным смещениям ткани распределение этих модулей в рассматриваемой области и, тем самым, найти характерные онкологические включения с повышенным значением модулей. Квазистатическая модель адекватно описывает исследуемые ткани, но обратная задача поиска механических модулей для нее оказывается достаточно сложной. Многочисленные численные эксперименты показали, что решить такую обратную задачу в реальном времени невозможно (см. [6]–[10]).

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

2. ПРЯМАЯ И ОБРАТНАЯ ЗАДАЧИ ДВУМЕРНОЙ ЭЛАСТОГРАФИИ В ПРИБЛИЖЕНИИ ПЛОСКОГО ДЕФОРМИРОВАННОГО СОСТОЯНИЯ

Будем считать, что исследуемый участок двумерного сечения биологической ткани моделируется областью $\Omega = ( - \infty ,\infty ) \times (0,h) \subset R_{{xy}}^{2}$, а сама ткань характеризуется распределением модуля Юнга $E = E(x,y)$ и постоянным коэффициентом Пуассона $\nu = 0.495$. Этот слой подвергается на нижней поверхности воздействию направленной вертикально вверх силы, а верхняя поверхность слоя неподвижна. Такие предположения характерны для постановок задач двумерной эластографии (см., например, [11]). В модели двумерного плоского линейного деформированного состояния упругого тела Ω (см., например, [5]–[8]) связь горизонтальных и вертикальных смещений ткани $u(x,y),w(x,y)$ с распределением модуля Юнга и коэффициентом Пуассона описывается системой уравнений в частных производных вида

(1)

$\begin{array}{*{20}{c}} {(1 - \nu )\frac{\partial }{{\partial x}}\left( {E\frac{{\partial u}}{{\partial x}}} \right) + \frac{1}{2}(1 - 2\nu )\frac{\partial }{{\partial y}}\left( {E\frac{{\partial u}}{{\partial y}}} \right) + \nu \frac{\partial }{{\partial x}}\left( {E\frac{{\partial w}}{{\partial y}}} \right) + \frac{1}{2}(1 - 2\nu )\frac{\partial }{{\partial y}}\left( {E\frac{{\partial w}}{{\partial x}}} \right) = 0,} \\ {\frac{1}{2}(1 - 2\nu )\frac{\partial }{{\partial x}}\left( {E\frac{{\partial u}}{{\partial y}}} \right) + \nu \frac{\partial }{{\partial y}}\left( {E\frac{{\partial u}}{{\partial x}}} \right) + \frac{1}{2}(1 - 2\nu )\frac{\partial }{{\partial x}}\left( {E\frac{{\partial w}}{{\partial x}}} \right) + (1 - \nu )\frac{\partial }{{\partial y}}\left( {E\frac{{\partial w}}{{\partial y}}} \right) = 0,} \\ {(x,y) \in \Omega ,} \end{array}$

с граничными условиями на ${{\Gamma }_{1}} = \{ (x,y):x \in ( - \infty ,\infty ),y = 0\} $ (нижняя граница полосы Ω):

(2)

$\begin{array}{*{20}{c}} {\frac{1}{2}(1 - 2\nu )E(x,y)\frac{{\partial u}}{{\partial y}} + \frac{1}{2}(1 - 2\nu )E(x,y)\frac{{\partial w}}{{\partial x}} = 0,} \\ {\nu E(x,y)\frac{{\partial u}}{{\partial x}} + (1 - \nu )E(x,y)\frac{{\partial w}}{{\partial y}} = (1 + \nu )(1 - 2\nu ){{f}_{y}}(x),} \end{array}$

и на ${{\Gamma }_{2}} = \{ (x,y):x \in ( - \infty ,\infty ),y = h\} $: $u(x,y) = w(x,y) = 0$ (верхняя граница полосы Ω). Здесь ${{f}_{y}}$ есть вертикальная компонента давления на поверхность ${{\Gamma }_{1}}$. Определение смещений $u(x,y),w(x,y)$ по известным коэффициентам $E(x,y),{\kern 1pt} {\kern 1pt} \nu $ составляет прямую задачу двумерной эластографии. При условиях ${{f}_{y}}(x) \in {{C}^{1}}({{\Gamma }_{1}}),{\kern 1pt} {\kern 1pt} E(x,y) \in {{C}^{1}}(\bar {\Omega })$ и $0 < {{E}_{1}} \leqslant E(x,y) \leqslant {{E}_{2}}$ она однозначно разрешима на классе функций $u(x,y),{\kern 1pt} w(x,y) \in W_{2}^{1}(\Omega )$ (это следует, например, из [12, гл. 2, 11], где получены более общие результаты). Здесь ${{E}_{{1,2}}}$ – известные константы. Соответствующие обратные задачи эластографии можно ставить по-разному, но все они заключаются в нахождении распределения модуля Юнга $E(x,y)$ по известным смещениям или функционалам от них (см. [5]–[11]). Ниже это будет уточнено.

Величину $\varepsilon = \frac{1}{2} - \nu \approx 0.005$ ($2\varepsilon \approx 0.01$) можно считать малым параметром. Тогда прямую задачу (1), (2) можно переписать в форме:

$\begin{gathered} \frac{\partial }{{\partial x}}\left( {E\frac{{\partial u}}{{\partial x}} + E\frac{{\partial w}}{{\partial y}}} \right) = - 2\varepsilon \left[ {\frac{\partial }{{\partial y}}\left( {E\frac{{\partial u}}{{\partial y}} + E\frac{{\partial w}}{{\partial x}}} \right) + \frac{\partial }{{\partial x}}\left( {E\frac{{\partial u}}{{\partial x}} - E\frac{{\partial w}}{{\partial y}}} \right)} \right], \hfill \\ \frac{\partial }{{\partial y}}\left( {E\frac{{\partial u}}{{\partial x}} + E\frac{{\partial w}}{{\partial y}}} \right) = - 2\varepsilon \left[ {\frac{\partial }{{\partial x}}\left( {E\frac{{\partial u}}{{\partial y}} + E\frac{{\partial w}}{{\partial x}}} \right) - \frac{\partial }{{\partial y}}\left( {E\frac{{\partial u}}{{\partial x}} - E\frac{{\partial w}}{{\partial y}}} \right)} \right],{\kern 1pt} \quad (x,y) \in \Omega , \hfill \\ \end{gathered} $

с граничными условиями

$\begin{gathered} {{\left( {E(x,y)\frac{{\partial u}}{{\partial x}} + E(x,y)\frac{{\partial w}}{{\partial y}}} \right)}_{{y = 0}}} = 4\varepsilon (1 + \nu ){{f}_{y}} + 2\varepsilon {{\left( {E(x,y)\frac{{\partial u}}{{\partial x}} - E(x,y)\frac{{\partial w}}{{\partial y}}} \right)}_{{y = 0}}},\quad (x,y) \in {{\Gamma }_{1}}, \hfill \\ \varepsilon {{\left( {\frac{{\partial u}}{{\partial y}} + \frac{{\partial w}}{{\partial x}}} \right)}_{{y = 0}}} = 0,{\kern 1pt} \quad (x,y) \in {{\Gamma }_{1}};\quad {\kern 1pt} u(x,y) = w(x,y) = 0,\quad (x,y) \in {{\Gamma }_{2}}. \hfill \\ \end{gathered} $

В дальнейшем функцию давления, фигурирующую в первом граничном условии, будем обозначать как $F(x) = 4\varepsilon (1 + \nu ){{f}_{y}}(x)$.

Выше была отмечена одна из важнейших проблем в практических приложениях прямой и обратной задач эластографии – получить решение быстро и достаточно точно. Тогда соответствующий алгоритм можно использовать в реальной диагностике. Численное решение прямой задачи (1), (2), основанное на методе конечных элементов, с этой точки зрения сравнительно эффективно и требует несколько секунд на ПК средней производительности при достаточно подробных сетках переменных x, y. Однако все известные нам методы решения двумерных обратных задач эластографии требуют значительно большего времени (от десятков минут до часов в зависимости от размеров сеток). Ниже предлагается метод приближенного решения прямой задачи как системы уравнений в частных производных с малым параметром. Полученные простые аналитические формулы для решений можно использовать для “быстрого” приближенного решения обратной задачи.

3. ИССЛЕДОВАНИЕ ПРЯМОЙ ЗАДАЧИ: НЕДООПРЕДЕЛЕННОСТЬ ЗАДАЧ ДЛЯ ПРИБЛИЖЕНИЙ МАЛОГО ПАРАМЕТРА, СОВМЕСТНОСТИ ЭТИХ ЗАДАЧ

Пусть $U = (u(x,y),w(x,y))$, где $u(x,y),w(x,y) \in W_{2}^{1}(\Omega )$. Введем дифференциальные операторы, действующие из $W_{2}^{1}(\Omega ) \times W_{2}^{1}(\Omega )$ в ${{L}_{2}}(\Omega )$, предполагая, что коэффициент $E(x,y)$ это гарантирует:

$\begin{gathered} {{L}_{0}}(U) = \left( {\frac{\partial }{{\partial x}}\left( {E\frac{{\partial u}}{{\partial x}} + E\frac{{\partial w}}{{\partial y}}} \right),\quad \frac{\partial }{{\partial y}}\left( {E\frac{{\partial u}}{{\partial x}} + E\frac{{\partial w}}{{\partial y}}} \right)} \right), \\ {{L}_{1}}(U) = \left( {\frac{\partial }{{\partial y}}\left( {E\frac{{\partial u}}{{\partial y}} + E\frac{{\partial w}}{{\partial x}}} \right) + \frac{\partial }{{\partial x}}\left( {E\frac{{\partial u}}{{\partial x}} - E\frac{{\partial w}}{{\partial y}}} \right),\quad \frac{\partial }{{\partial x}}\left( {E\frac{{\partial u}}{{\partial y}} + E\frac{{\partial w}}{{\partial x}}} \right) - \frac{\partial }{{\partial y}}\left( {E\frac{{\partial u}}{{\partial x}} - E\frac{{\partial w}}{{\partial y}}} \right)} \right). \\ \end{gathered} $

Определим также операторы граничных условий:

$\begin{gathered} {{l}_{0}}(U) = \left( {{{{\left( {E(x,y)\frac{{\partial u}}{{\partial x}} + E(x,y)\frac{{\partial w}}{{\partial y}}} \right)}}_{{y = 0}}},\;0} \right), \\ {{l}_{1}}(U) = \left( {{{{\left( {E(x,y)\frac{{\partial u}}{{\partial x}} - E(x,y)\frac{{\partial w}}{{\partial y}}} \right)}}_{{y = 0}}},\frac{1}{2}{{{\left( {\frac{{\partial u}}{{\partial y}} + \frac{{\partial w}}{{\partial x}}} \right)}}_{{y = 0}}}} \right). \\ \end{gathered} $

Тогда прямую задачу можно записать в форме:

(3)

$\begin{array}{*{20}{l}} {{{L}_{0}}(U) = - 2\varepsilon {{L}_{1}}(U),\quad (x,y) \in \Omega ,} \\ {{{l}_{0}}(U) = 2\varepsilon {{l}_{1}}(U) + (F(x),\,{\kern 1pt} 0),\quad {\kern 1pt} x \in ( - \infty ,\infty )} \\ {U{{{\text{|}}}_{{y = h}}} = 0.} \end{array},$

Решая (3) формально по методу малого параметра для регулярных возмущений (см., например, [13], [14]), получаем $U(x,y) = \sum\nolimits_{n = 0}^\infty {{(2\varepsilon )}^{n}}{{U}_{n}}(x,y)$. Здесь векторные функции ${{U}_{n}}(x,y)$, ${\kern 1pt} n = 0,1,2...{\kern 1pt} ,$ суть решения задач

$\begin{array}{*{20}{l}} {{{L}_{0}}({{U}_{0}}) = 0,{\kern 1pt} \quad (x,y) \in \Omega ,} \\ {{{l}_{0}}({{U}_{0}}) = (F,0),\quad x \in ( - \infty ,\infty ),} \\ {{{U}_{0}}{{{\text{|}}}_{{y = h}}} = 0,} \end{array}\quad {\kern 1pt} \begin{array}{*{20}{l}} {{{L}_{0}}({{U}_{{n + 1}}}) = - {{L}_{1}}({{U}_{n}}),{\kern 1pt} \quad (x,y) \in \Omega ,} \\ {{{l}_{0}}({{U}_{{n + 1}}}) = {{l}_{1}}({{U}_{n}}),\quad x \in ( - \infty ,\infty ),} \\ {{{U}_{{n + 1}}}{{{\text{|}}}_{{y = h}}} = 0,} \end{array}\quad n = 0,1,...{\kern 1pt} ,$

и предполагается, что все эти задачи имеют единственные решения в классе $W = W_{2}^{1}(\Omega ) \times W_{2}^{1}(\Omega )$ с нормой $\left\| U \right\|_{W}^{2} = \left\| u \right\|_{{W_{2}^{1}}}^{2} + \left\| w \right\|_{{W_{2}^{1}}}^{2}$. Обозначим через $P(U)$ линейный оператор решения задач второго типа из этой группы. Тогда

(4)

$\begin{gathered} {{U}_{1}} = P({{U}_{0}}),\quad {{U}_{2}} = P({{U}_{1}}) = {{P}^{2}}({{U}_{0}}),{\kern 1pt} \quad ...,\quad {\kern 1pt} {{U}_{n}} = {{P}^{n}}({{U}_{0}}),\quad ...{\kern 1pt} \Rightarrow \\ U = {{U}_{0}} + \sum\limits_{n = 1}^\infty {{(2\varepsilon )}^{n}}{{P}^{n}}({{U}_{0}}) = \sum\limits_{n = 0}^\infty {{(2\varepsilon )}^{n}}{{P}^{n}}({{U}_{0}}). \\ \end{gathered} $

Ряд Неймана (4) сходится, если формально $2\varepsilon {\kern 1pt} ||{\kern 1pt} P({{U}_{0}}){\kern 1pt} {{||}_{W}}\; < 1$. В этом случае для приближенного решения задачи (3) надо решить задачу для ${{U}_{0}}$ и затем проделать несколько итераций по вычислению величин ${{U}_{n}} = P({{U}_{{n - 1}}}),$, $n = 1,2,...$.

Задача для ${{U}_{0}}$ имеет следующий вид:

(5)

$\begin{array}{*{20}{l}} {\frac{\partial }{{\partial x}}\left( {E\frac{{\partial u}}{{\partial x}} + E\frac{{\partial w}}{{\partial y}}} \right) = 0,\quad \frac{\partial }{{\partial y}}\left( {E\frac{{\partial u}}{{\partial x}} + E\frac{{\partial w}}{{\partial y}}} \right) = 0,{\kern 1pt} \quad (x,y) \in \Omega ,} \\ {{{{\left( {E(x,y)\frac{{\partial u}}{{\partial x}} + E(x,y)\frac{{\partial w}}{{\partial y}}} \right)}}_{{y = 0}}} = F(x),\quad (x,y) \in {{\Gamma }_{1}};} \\ {u(x,y) = w(x,y) = 0,{\kern 1pt} \quad (x,y) \in {{\Gamma }_{2}}.} \end{array}$

Из первых двух уравнений получим

(6)

$E(x,y)\frac{{\partial u}}{{\partial x}} + E(x,y)\frac{{\partial w}}{{\partial y}} = C = {\text{const}},\quad (x,y) \in \Omega ,$

а из граничного условия на ${{\Gamma }_{1}}$ тогда следует

${{\left( {E(x,y)\frac{{\partial u}}{{\partial x}} + E(x,y)\frac{{\partial w}}{{\partial y}}} \right)}_{{y = 0}}} = C = F(x)\quad {\kern 1pt} \forall x$.

Это значит, что задача для ${{U}_{0}}$ разрешима только для постоянных давлений: $F(x) = {{F}_{0}}$. Даже при таком предположении задача оказывается недоопределенной, так как из (6) имеем

$\begin{array}{*{20}{l}} {E(x,y)\frac{{\partial u}}{{\partial x}} + E(x,y)\frac{{\partial w}}{{\partial y}} = {{F}_{0}},\quad (x,y) \in \Omega ,} \\ {u(x,y) = w(x,y) = 0,{\kern 1pt} \quad (x,y) \in {{\Gamma }_{2}},} \end{array}$

и найти однозначно $u(x,y),w(x,y)$ нельзя. Отсюда следует, что для предлагаемой схемы решения краевой задачи (1), (2) в задачу для ${{U}_{0}}$ нужно вводить дополнительные предположения. Аналогичные проблемы возникают при решении задач для ${{U}_{n}},\;{\kern 1pt} n > 0$.

4. ДОПОЛНИТЕЛЬНЫЕ ПРЕДПОЛОЖЕНИЯ. АНАЛИТИЧЕСКОЕ РЕШЕНИЕ ЗАДАЧ ДЛЯ ПРИБЛИЖЕНИЙ МАЛОГО ПАРАМЕТРА. РЕШЕНИЕ ОБРАТНОЙ ЗАДАЧИ

Будем дополнительно предполагать, что горизонтальные смещения $u(x,y)$ “малы”: формально $u(x,y) = 0$, $(x,y) \in \Omega $. Тогда из (5) получим для нулевого приближения ${{w}_{0}}$:

$\begin{array}{*{20}{l}} {\frac{\partial }{{\partial x}}\left( {E\frac{{\partial {{w}_{0}}}}{{\partial y}}} \right) = 0,\quad \frac{\partial }{{\partial y}}\left( {E\frac{{\partial {{w}_{0}}}}{{\partial y}}} \right) = 0,\quad (x,y) \in \Omega ,} \\ {{{{\left( {E(x,y)\frac{{\partial {{w}_{0}}}}{{\partial y}}} \right)}}_{{y = 0}}} = F(x),\quad (x,y) \in {{\Gamma }_{1}};\quad {{w}_{0}}(x,y) = 0,{\kern 1pt} \quad (x,y) \in {{\Gamma }_{2}}.} \end{array}$

Отсюда, как и при получении равенства (6), выводим, что $F(x) = {{F}_{0}}$ и

(7)

$E(x,y)\frac{{\partial {{w}_{0}}}}{{\partial y}} = {{F}_{0}},\quad (x,y) \in \Omega ;\quad {{w}_{0}}(x,h) = 0.$

Поэтому нулевое приближение вертикального смещения есть

(8)

${{w}_{0}}(x,y) = {{F}_{0}}\int\limits_h^y {\frac{{d\eta }}{{E(x,\eta )}}} .$

Формулы (8) и (7) можно использовать при приближенном решении прямой и обратной задач эластографии в нулевом приближении.

Задача для следующего приближения ${{w}_{1}}(x,y)$ принимает вид

$\begin{array}{*{20}{l}} {{{L}_{0}}({{U}_{1}}) = - {{L}_{1}}({{U}_{0}}),\quad {{U}_{0}} = (0,{{w}_{0}}),\quad {{U}_{1}} = (0,{{w}_{1}}),} \\ {{{l}_{0}}({{U}_{1}}) = {{l}_{1}}({{U}_{0}}),\quad x \in {{\Gamma }_{1}},} \\ {{{U}_{1}}{{{\text{|}}}_{{y = h}}} = 0} \end{array}$

или

(9)

$\begin{array}{*{20}{l}} {\frac{\partial }{{\partial x}}\left( {E(x,y)\frac{{\partial {{w}_{1}}}}{{\partial y}}} \right) = - \frac{\partial }{{\partial y}}\left( {E(x,y)\frac{{\partial {{w}_{0}}}}{{\partial x}}} \right),{\kern 1pt} \quad \frac{\partial }{{\partial y}}\left( {E(x,y)\frac{{\partial {{w}_{1}}}}{{\partial y}}} \right) = - \frac{\partial }{{\partial x}}\left( {E(x,y)\frac{{\partial {{w}_{0}}}}{{\partial x}}} \right),{\kern 1pt} \quad (x,y) \in \Omega ,} \\ {{{{\left. {E(x,y)\frac{{\partial {{w}_{1}}}}{{\partial y}}} \right|}}_{{y = 0}}} = - E(x,y)\frac{{\partial {{w}_{0}}}}{{\partial y}}{{|}_{{y = 0}}} = - {{F}_{0}},\quad 0 = \frac{1}{2}{{{\left. {\frac{{\partial {{w}_{0}}}}{{\partial x}}} \right|}}_{{y = 0}}},\quad {{{\left. {{{w}_{1}}} \right|}}_{{y = h}}} = 0.} \end{array}$

Фигурирующие здесь уравнения в частных производных совместны не всегда. Например, для достаточно гладкого коэффициента $E(x,y)$, точнее при $E(x,y) \in {{C}^{3}}(\Omega )$, необходимо выполнение условия

(10)

$\frac{{{{\partial }^{2}}}}{{\partial {{x}^{2}}}}\left( {E(x,y)\frac{{\partial {{w}_{0}}}}{{\partial x}}} \right) = \frac{{{{\partial }^{2}}}}{{\partial {{y}^{2}}}}\left( {E(x,y)\frac{{\partial {{w}_{0}}}}{{\partial x}}} \right),\quad (x,y) \in \Omega .$

Равенство (10) есть ограничение на класс допустимых функций $E(x,y)$, которое достаточно трудно интерпретируется на практике. Сделаем более ясное дополнительное предположение, гарантирующее выполнение условия (10) – считаем, что модуль Юнга слабо зависит от x в следующем смысле:

(11)

$\frac{{\partial {{w}_{0}}}}{{\partial x}} = \frac{\partial }{{\partial x}}\left( {{{F}_{0}}\int\limits_h^y {\frac{{d\eta }}{{E(x,\eta )}}} } \right) \approx 0,{\kern 1pt} \quad (x,y) \in \Omega .$

Тогда задача (9) переходит в задачу

$\begin{array}{*{20}{l}} {\frac{\partial }{{\partial x}}\left( {E(x,y)\frac{{\partial {{w}_{1}}}}{{\partial y}}} \right) = 0,{\kern 1pt} \quad \frac{\partial }{{\partial y}}\left( {E(x,y)\frac{{\partial {{w}_{1}}}}{{\partial y}}} \right) = 0,\quad (x,y) \in \Omega ,} \\ {{{{\left. {E(x,y)\frac{{\partial {{w}_{1}}}}{{\partial y}}} \right|}}_{{y = 0}}} = - {{F}_{0}},\quad {{{\left. {{{w}_{1}}} \right|}}_{{y = h}}} = 0,} \end{array}$

схожую с задачей для ${{w}_{0}}$. Поэтому ${{w}_{1}} = - {{w}_{0}} = - {{F}_{0}}\int_h^y \frac{{d\eta }}{{E(x,\eta )}}$. Аналогичные рассуждения можно провести для всякого ${{w}_{n}},n > 1:{\kern 1pt} {\kern 1pt} {{w}_{n}} = - {{w}_{{n - 1}}}$. При этом используется одно и то же предположение (11) о модуле Юнга. Соответственно, ряд Неймана (4) будет иметь вид

(12)

$w = {{w}_{0}} + 2\varepsilon {{w}_{1}} + {{(2\varepsilon )}^{2}}{{w}_{2}} + ... + {{(2\varepsilon )}^{n}}{{w}_{n}} + ... = {{w}_{0}}\sum\limits_{n = 0}^\infty {{( - 1)}^{n}}{{(2\varepsilon )}^{n}} = \frac{{{{w}_{0}}}}{{1 + 2\varepsilon }} = \frac{{{{F}_{0}}}}{{(1 + 2\varepsilon )}}\int\limits_h^y {\frac{{d\eta }}{{E(x,\eta )}}} .$

Таким образом, при сделанных предположении малости смещений $u(x,y)$ и предположении (11) нулевое приближение рассматриваемой прямой задачи с точностью до нормирующего множителя $\frac{1}{{(1 + 2\varepsilon )}} \approx 0.99$ совпадает с ее точным решением. Это позволяет использовать уравнения (8), (12) для приближенного нахождения функции $E(x,y)$ при решении обратной задачи. Например, это можно сделать в форме

(13)

$E(x,y) = \frac{{{{F}_{0}}}}{{(1 + 2\varepsilon )}}{{\left( {\frac{{\partial w}}{{\partial y}}} \right)}^{{ - 1}}},{\kern 1pt} \quad (x,y) \in \Omega ,$

если вычислять частную производную с помощью регуляризованных методов, устойчивых к возмущениям данных $w(x,y)$ обратной задачи.

Отметим, что при эластографической диагностике часто важно знать не абсолютные значения модуля Юнга, а отношение $E(x,y){\text{/}}{{E}_{0}}$, где ${{E}_{0}}$ – характерная величина модуля Юнга здоровой ткани. Существенное превышение единицы для такого отношения есть признак, на который следует обратить внимание при диагностике.

5. ЧИСЛЕННОЕ РЕШЕНИЕ ПРЯМОЙ И ОБРАТНОЙ ЗАДАЧ

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

(14)

$\begin{gathered} E(x,y) = 32 + 100\exp \left( { - \frac{{{{{(x + 0.2)}}^{2}} + {{{(y + 0.01)}}^{2}}}}{{0.05}}} \right) + 152\exp \left( { - \frac{{{{{(x - 0.2)}}^{2}} + {{{(y + 0.15)}}^{2}}}}{{0.053}}} \right) + \\ + \;50\exp \left( { - \frac{{{{{(y - 2(x + 0.9))}}^{2}}}}{{0.085}}} \right). \\ \end{gathered} $

Для распределения (14) решим прямую задачу нахождения смещений $w(x,y)$ по заданному коэффициенту $E(x,y)$ в “перевернутой и сдвинутой” области конечного размера $\Omega = [ - 1,1] \times [ - 0.5,0.5]$. Полагаем, что к верхней границе этой области приложено давление ${{F}_{0}} = 1$, а левый и правый края области свободны. Решение соответствующей задачи (1), (2), полученное с помощью метода конечных элементов с высокой точностью (“точное решение”), представлено на фиг. 1а. Для сравнения найдем приближенное решение прямой задачи по формуле (12), которая получена для бесконечной полосы при указанных в разд. 4 дополнительных предположениях (фиг. 1б). Видно, что оно, по крайней мере качественно, отражает основные особенности точного решения $w(x,y)$ прямой задачи.

Фиг. 1.

(а) – “Точное” решение прямой задачи (1), (2) с помощью метода конечных элементов. (б) – Приближенное решение ${{w}_{0}}(x,y)$ прямой задачи по формуле (12).

Теперь решим обратную задачу: по данным $w(x,y)$, которые изображены на фиг. 1а, т.е. по решению точной прямой задачи (1), (2), восстановим модуль Юнга. Для этого применим приближенную формулу (13), используя в ней регуляризованную процедуру численного дифференцирования функции $w(x,y)$ (см., например, [15], [16]). Результат в сравнении с точным распределением модуля (14), нормированным на фоновое значение модуля Юнга ${{E}_{0}} = \min \{ E(x,y):{\kern 1pt} (x,y) \in \Omega \} $, показан на фиг. 2. Полученное приближенное решение обратной задачи с качественной точки зрения достаточно хорошо отражает структуру неоднородности модуля Юнга. С количественной точки зрения оно с удовлетворительной точностью представляет величины отношения $E(x,y){\text{/}}{{E}_{0}}$.

Фиг. 2.

Сравнение точного (а) и восстановленного (б) распределений модуля Юнга при приближенном решении обратной задачи по методу малого параметра. Модули нормированы на фоновое значение ${{E}_{{{\text{min}}}}}$.

Расчет по формуле (13) на персональном компьютере средней производительности занимает 0.5 мс для сеток размера 100 × 100. Скорость расчета и достаточная точность получаемого приближенного решения позволяют надеяться, что этот подход можно использовать в реальной онкологической диагностике. Можно также использовать данное приближенное решение в качестве начального приближения при решении обратной задачи двумерной квазистатической эластографии для исходной системы (1), (2) по методам из работ [8]–[10].

Список литературы

Gao L., Parker K., Lerner R. et al. Imaging of the elastic properties of tissue – a review // Ultrasound Med. Biol. 1996. V. 22. P. 959–977.
Ophir J., Alam S., Garra B. et al. Elastography: ultrasonic estimation and imaging of the elastic properties of tissues // Proc. Inst. Mech. Eng. Part H: J. Eng. Med. 1999. V. 213. P. 203–233.
Greenleaf J.F., Fatemi M., Insana M. Selected methods for imaging elastic properties of biological tissues // Annu. Rev. Biomed. Eng. 2003. V. 5. P. 57–78.
Parker K.J., Taylor L.S., Gracewski S. et al. A unified view of imaging the elastic properties of tissue // J. Acoust. Soc. Am. 2005. V. 117. P. 2705–2712.
Doyley M. Model-based elastography: a survey of approaches to the inverse elasticity problem // Phys Med Biol. 2012. V. 57. P. R35–R73.
Oberai A.A., Gokhale N.H., Feijoo G.R. Solution of inverse problems in elasticity imaging using the adjoint method // Inverse Probl. 2003. V. 19. P. 297–313.
Richards M., Barbone P., Oberai A. Quantitative three-dimensional elasticity imaging from quasi-static deformation: a phantom study // Phys. Med. Biol. 2009. V. 54. P. 757–779.
Leonov A.S., Sharov A.N., Yagola A.G. A posteriori error estimates for numerical solutions to inverse problems of elastography // Inverse Probl. Sci. Eng. 2017. V. 25. P. 114–128.
Leonov A.S., Sharov A.N., Yagola A.G. Solution of the inverse elastography problen for parametric classes of inclusions with a posteriori error estimate // J. Inverse Ill-Posed Probl. 2017. V. 26. P. 1–7.
Leonov A.S., Sharov A.N., Yagola A.G. Solution of the three-dimensional inverse elastography problem for parametric classes of inclusions // Inverse Probl. Sci. Eng. 2021. V. 29. N. 8. P. 1055–1069.
Rychagov M., Khaled W., Reichling S. et al. Numerical modeling and experimental investigation of biomedical elastographic problem by using plane strain state model // Fortsch. Der Akustik. 2003. V. 29. P. 586–589.
Ладыженская О.А. Краевые задачи математической физики. М.: Наука, 1973.
Тихонов А.Н., Васильева А.Б., Свешников А.Г. Дифференц. ур-ния. М.: Наука, 1980.
Треногин В.А. Функциональный анализ. М.: Наука, 1980.
Тихонов А.Н., Гончарский А.В., Степанов В.В., Ягола А.Г. Численные методы решения некорректных задач. М.: Наука, 1990.
Леонов А.С. Решение некорректно поставленных обратных задач. Очерк теории, практические алгоритмы и демонстрации в МАТЛАБ. М.: Либроком, 2009.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Журнал вычислительной математики и математической физики

Архивы выпусков Информация о журнале Отправить рукопись в журнал