автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Численные методы решения задачи электроимпедансной томографии в случае кусочно-постоянной проводимости

Московский государственный университет имени М.В.Ломоносова факультет вычислительной математики и кибернетики

На правах рукописи

005533291

Гаврилов Сергей Вадимович

Численные методы решения задачи электроимпедансной томографии в случае кусочно-постоянной проводимости

Специальность 05.13.18 - «Математическое моделирование, численные методы и комплексы программ»

АВТОРЕФЕРАТ

диссертации на соискание ученой степени кандидата физико-математических наук

1У СЕН ДШ

Москва - 2013

005533291

Работа выполнена на кафедре математической физики факультета вычислительной математики и кибернетики Московского государственного университета имени М.В.Ломоносова.

Научный доктор физико-математических наук, профессор,

руководитель: заведующий кафедрой математической физики факультета ВМК МГУ имени М.В.Ломоносова Денисов Александр Михайлович. Официальные доктор физико-математических наук, профессор оппоненты: кафедры математики физического факультета МГУ

имени М.В.Ломоносова Боголюбов Александр Николаевич;

доктор физико-математических наук, начальник научно-исследовательской лаборатории ФГУП «Всероссийский научно-исследовательский институт автоматики им. Н.Л.Духова» Сафронов Сергей Иванович.

Ведущая федеральное государственное бюджетное

организация: образовательное учреждение высшего

профессионального образования «Московский авиационный институт (национальный

исследовательский университет)».

Защита состоится « 9 » октября 2013 г. в 15 ч. 30 мин. на заседании диссертационного совета Д 501.001.43 при Московском государственном университете имени М.В.Ломоносова, расположенном по адресу: 119991, ГСП-1, Москва, Ленинские горы, МГУ, 2-й учебный корпус, факультет вычислительной математики и кибернетики, ауд. 685.

С диссертацией можно ознакомиться в фундаментальной библиотеке МГУ имени М.В.Ломоносова.

Автореферат разослан «_»_2013г.

Ученый секретарь диссертационного совета,

доктор физико-математических наук, о у

профессор Захаров Е.В.

Общая характеристика работы

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

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

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

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

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

4

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

Цель и задачи работы

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

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

3. Программная реализация предложенных методов и проведение вычислительных экспериментов с целью определяют их эффективности.

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

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

электроимпедансной томографии в диагностике онкологических и кардиологических заболеваний.

Апробация работы. Основные результаты работы докладывались на следующих конференциях и семинарах.

1. V международная конференция "Математические идеи ПЛ.Чебышева и их приложение к современным проблемам естествознания"(Обнинск, 14-18 мая 2011 года).

2. Международная конференция, посвященная 80-летию со дня рождения академика М.М.Лаврентьева "Обратные и некорректные задачи математической физики"(Новосибирск, 5-12 августа 2012 года).

3. Научная конференция "Тихоновские чтения"(Москва, МГУ имени М.В.Ломоносова, 29.10.2012-02.11.2012).

4. Научно-исследовательский семинар кафедры математической физики факультета ВМК МГУ имени М.В.Ломоносова.

5. 4-я Международная конференция «Функциональные пространства. Дифференциальные операторы. Общая топология. Проблемы математического образования.», посвящённая 90-летию со дня рождения члена-корреспондента РАН, академика Европейской академии наук Л.Д. Кудрявцева (Москва, РУДН, 25-29 марта 2013 года).

6. Научный семинар "Обратные задачи математической физики"под руководством профессоров А.Б.Бакушинского, А.В.Тихонравова и А.Г.Яголы в Научно-исследовательском вычислительном центре МГУ имени М.В.Ломоносова.

7. Семинар по математическим методам в естественных науках под руководством профессора А.Н. Боголюбова на физическом факультете МГУ имени М.В.Ломоносова.

8. Научный семинар "Моделирование и идентификация тепловых процессов в технике"под руководством чл.-корр. РАН, д.т.н., профессора О.М.Алифанова в федеральном государственном бюджетном образовательном учреждении высшего профессионального образования «Московский авиационный институт (национальный исследовательский университет)».

Публикации. Материалы диссертации опубликованы в 8 печатных работах, из них 4 статьи в рецензируемых журналах, рекомендованных ВАК [1, 2, 4, 7] статья [6] и 3 тезиса докладов на конференциях [3, 5, 8].

Структура диссертации. Диссертация состоит из введения, обзора литературы, 3 глав, заключения и библиографии. Общий объем диссертации 110 страниц, включая 20 рисунков. Библиография включает 94 наименований на 11 страницах.

Общая характеристика работы

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

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

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

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

Математическая постановка задачи такова. На плоскости рассматривается односвязная ограниченная область О, с границей Го. Область 0.1 — односвязная область с границей Гь такая, что е П. Границы Го и Г1 достаточно гладкие. Через По обозначается область

Рассматривается следующая задача. В предположении, что граница Го известна и со, о\ — заданные положительные постоянные, требуется определить неизвестную границу Г1 и функцию и{х,у), такие что: и € С(й),и(х,у) = щ{х,у), (х,у) е = 0,1), где и0 6 С2(П0) Л Сг1(Щ,'"1 е С2{Пг) Л С1^)-

(1)

(2)

(3)

(4)

(5)

&щ(х,у) = 0, (ж,1/)еПг,1 = 0,1;

Щ{х,у) = П1(х,у), (х,у) е Гь его—' (».УЗбГь

ио(х,у) = /(ж,у), (х,у)£Г0,

дио(х,у) . . , . _ дп

где /(х, у), д(х, у), (х, у) е Г0 - заданные функции.

Задача (1)-(5) описывает модель среды с кусочно-постоянной проводимостью.

В §1.1 разрабатывается численный метод решения поставленной задачи, основанный на интерпретации уравнений (1)-(5) как обратной задачи к задаче Дирихле (1)-(4), в которой требуется определить неизвестные Г\ и и(х,у) по дополнительному условию (5). Эта обратная задача в диссертационной работе называется обратной задачей Дирихле-Неймана.

В §1.1.1 выводится операторное уравнение для функции, определяющей неизвестную границу Гх. Для параметризации границы Г\ делается предположение, что класс неизвестных кривых Г\ таков, что известна точка М0, являющаяся общим центром звездности для кривых 1\ из этого класса, кроме того, кривые Гх задаются в полярной системе координат с центром в точке Мо функциями г(ф) : r(ip) е С2 [0,27г] и 11г11с2[0 2тг] — Со. гДе со - фиксированное число.

Решение задачи Дирихле (1)-(4) ищется в виде суммы двух потенциалов простого слоя

и(М) = [ Р) In (—) dip + (сто - ах) [ v{P) In (—) dip, (6) J \PmpJ J \PmpJ

Го Ti

где точка M с координатами (х, у) принадлежит области

С использованием введенного представления (6), уравнений (4), (3) и свойств потенциалов простого слоя выводится система интегральных уравнений для плотностей ц(Р) и v{P)

[ ц{Р) In (—) dip + (<т0 - ox) / v(P) In (—) dip = f{M), M 6 Го, J \PmpJ J \pMpJ

Го Г!

■kv(M) + , , [ v(P)ir- Ь (—) dlP+ (00 + 01) У дпт \pMpJ

Г° (8) +&Lpi f v(P)J- in (J-) dlP = 0, Me Гх.

Гх

где nm - внутренняя нормаль к кривой Гх в точке М.

С использованием условия (5), представления (6) и свойств потенциала

простого слоя строится уравнение

-мм)+j

Г° (9)

+(<70 - <гх) [ "(P)J~1п (■—) dlP = 9{М), м е Го, J dnm \pmpJ

Г1

где Пт - внутренняя нормаль к кривой Го в точке М.

Уравнения (7), (8), (9), записанные в полярной системе координат, определяют нелинейное операторное уравнение относительно неизвестной функции г(ф)

Ar = g (10)

Для вычисления функции (Аг)(ф) при заданной г(ф), определяющей границу Гх, нужно решить систему интегральных уравнений (7), (8) и определить плотности потенциалов v{P), а затем вычислить значение оператора, стоящего в левой части уравнения (9).

В §1.1.2 строится итерационный метод решения обратной задачи Дирихле-Неймана. В качестве начального приближения г0(ф) неизвестного контура Гх принимается окружность. Переход от приближения гп(ф), полученного на n-ом шаге итерационного процесса к приближению гп+\(ф) осуществляется по следующему алгоритму.

В окрестности функции гп(ф) уравнение (10) линеаризуется и

получается линейное уравнение для функции рп{"Ф), представляющей собой поправку к функции гп{ф)\

В[гп]рп = дп. (11)

Решением этого уравнение является функция рп(Ф), с которой вычисляется гп+1(ф) = гп{ф) + рп(ф), п = 0,1,2... . (12)

При численной реализации итерационного метода на отрезке [0,2ж\ вводится сетка ^ и сеточные аналоги ггп функций гп(ф). После замены интегралов, входящих в выражения, определяющие оператор В[гп), на квадратурные формулы задача решения уравнения (11) сводится к решению системы линейных алгебраических уравнений и для решения полученной системы линейных алгебраических уравнений применяется метод регуляризации Тихонова. Величина параметра регуляризации согласовывается с точностью задания исходной информации и с шагом итерационного процесса п.

В §1.2 разрабатывается численный метод решения двумерной задачи электроимпедансной томографии, основанный на интерпретации уравнений (1)-(5) как обратной задачи к задаче Неймана (1)-(3), (5), в которой требуется определить Г\ и и(х,у) по дополнительному условию (4). Эта обратная задача в диссертационной работе называется обратной задачей Неймана-Дирихле.

В §1.2.1 выводится операторное уравнение для функции, определяющей неизвестную границу Гь Для класса кривых, в котором ищется граница Гь аналогично §1.1.1 делается предположение о существовании единого центра звездности и граница Г1 параметризуется в

полярной системе координат функцией г{ф). Для решения задачи Неймана (1)-(3), (5) используется представление в виде суммы потенциалов простого слоя (6). Это представление позволяет построить систему интегральных уравнений для нахождения плотностей потенциалов при фиксированном контуре Гь в которую входят условия (3), (5) и дополнительное условие, задающее значение решения и(М) в некоторой фиксированной точке Мх б Г0.

Уравнения для определения плотностей потенциалов простого слоя в совокупности с условием (4), записанным с помощью представления (6) определяют нелинейное операторное уравнение относительно неизвестной функции г(^).

В §1.2.2 строится итерационный метод решения обратной задачи Неймана-Дирихле. Принципы построения итерационного метода аналогичны §1.1.2. В окрестности функции гп(ф) — приближения неизвестного контура, полученного на п-ом шаге итерационного процесса, нелинейное операторное уравнение для функции г(ф) линеаризуется и получается линейное уравнение для функции Рп(.Ф), представляющей собой поправку к функции гп(ф). Это уравнение решается с применением метода регуляризации Тихонова.

В §1.3 разрабатывается третий численный метод решения двумерной задачи электроимпедансной томографии для кусочно-постоянной проводимости в случае одного измерения на границе.

Метод основан на введении вспомогательного контура, относительно близкого к Го, и включает в себя два этапа. На первом решается задача Коши для уравнения Лапласа и находится функция

щ(М) = /2(М), М еТ2 (13)

на вспомогательном контуре Г2, таком что Г2 € П0 и Г2 содержит внутри себя контур IV На втором этапе с найденной функцией /2(М) решается задача (1)-(4), (13).

В §1.3.1 приводится численный метод решения задачи Коши для уравнения Лапласа для нахождения значения потенциала на вспомогательном контуре. Математическая постановка задачи такова. Рассматривается односвязная область Пг. ограниченная снаружи контуром Го и изнутри контуром Гг. Требуется определить значение функции щ(М),М £ Гг, если известно, что в области П2 функция и0(М) является решением задачи Коши для уравнения Лапласа

Ди0(М) = 0, М б П2, (14)

«о(М) = /(М), Мб Г0) (15)

^ = мег0. (16)

Строится численный метод решения задачи (14)-(16), основанный на методе граничных интегральных уравнений [13, 9] и методе регуляризации Тихонова [И].

В §1.3.2 разрабатывается численный метод определения границы Гх в задаче (1)-(4),(13). Метод основан на представлении решения краевой задачи в виде суммы потенциалов простого и двойного слоя, построении операторного уравнения для неизвестной границы и решении его итерационным методом, на каждом шаге которого полученное операторное уравнение линеаризуется и его решение сводится к решению системы линейных алгебраических уравнений, при котором применяется метод регуляризации Тихонова.

В §1.4 представлены сведения о программной реализации

Рисунок 1. Результаты вычислительного эксперимента для метода Дирихле-Неймана

предложенных численных методов и приведены результаты вычислительных экспериментов, проведенных с использованием разработанного комплекса программ в среде программирования МАТЬАВ. В §1.4 представлена схема проведения вычислений и полученные результаты для каждого из трех методов. На рис. 1 показаны результаты одного из вычислительных экспериментов, в котором для решения задачи элекгроимпедансной томографии использовалась программная реализация численного метода решения обратной задачи Дирихле-Неймана. В этом вычислительном эксперименте были выбраны сг0 = -3, аг = 1. Функция /(ф) имела вид

/(ф) = 50 ^ехр[—4эш2(ф/2)\ - ехр[-4со82(^/2)]^.

;

С заданными Г0, Гь сг0, <У\ и ¡(ф) решалась задача Дирихле (1)-(4)

I

и находилась функция д(ф). В эту функцию вносилась погрешность и получалась функция д5(ф), такая что ||д(ф) - дб{ф)\\ь2[о121г]/ЫШ121оМ = | 0,05. Функция дй(ф) использовалась в качестве исходных данных для численного решения обратной задачи. На рис. 1 приведены точный контур

|

Гх, начальное приближение Г? и результат решения обратной задачи rj4, полученный при выборе равномерных сеток на контуре Г0 в 100 узлов и на контуре Гх в 60 узлов, в результате 14-ти итераций с начальным приближением Г°. Критерием останова служило достижение уровня погрешности по невязке.

В рамках проведенных вычислительных экспериментов все три предложенных численных метода показали схожую эффективность.

Результаты первой главы были опубликованы в работах [1, 2, 3, 6].

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

Математическая постановка задачи такова. В пространстве рассматривается односвязная ограниченная область í"2 с границей Го- fii — односвязная область с границей Гх, такая, что Í2x € П. Поверхности Го и Гх достаточно гладкие. Через П0 обозначается область íl0 = П \

Рассматривается следующая краевая задача. Требуется определить функцию и{М), такую что: и G С(Щ,и(М) = щ(М),М £ = 0,1), где щ е C2(íí¿) П = 0,1),

Здесь 0о, ох заданные положительные постоянные, а f{M) известная функция, непрерывная и не постоянная на Го-

Задача электроимпедансной томографии формулируется следующим

А щ(М) = 0, Meüi, i = 0,1

(17)

(18)

м е Гх,

(19)

(20)

образом. Пусть в краевой задаче (17)-(20) поверхность Го, постоянные сто, сг\ и функция /(М) на Го заданы, а поверхность Г1 неизвестна. Требуется определить Гь если задана дополнительная информация о решении задачи (17М20):

ме Г„, (21)

где д(М) известная функция, непрерывная на Г0, а п - внутренняя нормаль к Го.

В §2.1 выводится операторное уравнение для функции, определяющей неизвестную поверхность Г1. Делается предположение, что класс неизвестных поверхностей Гх таков, что известна точка Мо, являющаяся общим центром звездности для поверхностей Гх из этого класса, кроме того, поверхности Гх задаются в сферической системе координат с центром в точке М0 функциями г(0,</з) : г (в, ¡р) £ С2{[0,7г] х [0,2т:]} и 1М1с2{[о,7г]х[о,2тг]} ^ Со, где со - фиксированное число.

Для решения краевой задачи (17)-(20) вводится представление в виде суммы двух потенциалов простого слоя

и(М) = [[ ц{Р)—йвр + (ао - <71) /У и(Р)—{[8р, Меп, (22) 3 ] Рмр Рмр

То Г!

где рмр ~ расстояние между точками М и Р.

С использованием представления (22), условий (19), (20) и свойств потенциала простого слоя получается система интегральных уравнений для плотностей р(Р), и{Р)- Условие (21), записанное с учетом представления (22), в совокупности с системой уравнений для плотностей потенциалов простого слоя определяет нелинейное операторное уравнение относительно неизвестной функции г(9,ф)

иг = д. (23)

В §2.2 строится итерационный метод решения операторного уравнения (23). В качестве начального приближения го(9,<р) для неизвестной поверхности принимается сфера. На каждом шаге итерационного метода уравнение (23) линеаризуется в окрестности функции г„(в,(р), получается линейное операторное уравнение относительно поправки к приближению r„(9, ip), которое решается с применением метода регуляризации Тихонова.

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

Численная реализация итерационного метода проводится аналогично

итерационному методу решения задачи Дирихле-Неймана в двумерном

случае. На рис. 2 приведены результаты одного из вычислительных

экспериментов, проведенных с использование разработанного комплекса

программ. Схема проведения вычислительного эксперимента была такова.

С заданными Г0, Гь сто, и f(6,tp) решалась задача Дирихле (17)-

(20) и находилась функция g(9,tp), представляющая собой значение „ ди(М)

нормальной производной ^ на поверхности Г0. В эту функцию вносилась погрешность и получалась функция gs(6,<p), такая что II9(6, <р) - ffi(6,Iy)IU2{[0,T]x[0,27r]}/||5(él,^)IU2{MX[0,2^} = 0,01. Функция gs(0, <р) использовалась в качестве исходных данных для численного решения обратной задачи. На рис. 2 приведены некоторые результаты решения обратной задачи. Расчеты проводились при выборе равномерных сеток на Г0 {(0<,Pj);i = 0,1,..,31;j = 0,1, ..,63} и на 1\ {{ви^)\г = 0,1,.., 17; j = 0,1,.., 35}. Критерием останова служило достижение уровня погрешности по невязке.

Результаты второй главы были опубликованы в работах [4, 5].

В третьей главе разрабатываются численные методы решения задачи

-40

-30

-20

-10

20

30

10

-50 -40 -30 -20 -10

0 10 20 30 40 50

Рисунок 2. Результаты вычислительного эксперимента. Кривые Fj, Г^ и Г]4 — сечения точной поверхности Ti, начального приближения и найденной в результате 14 итераций поверхности FJ4 соответственно, (а) — сечение плоскостью, содержащей главную ось эллипса и ось OZ; (б) — сечение плоскостью, проходящей через начало координат и перпендикулярной главной оси симметрии эллипса; (в) — сечение плоскостью, перпендикулярной оси OZ и проходящей через начало координат.

электроимпедансной томографии для кусочно-постоянной проводимости в случае нескольких измерений на внешней границе.

Математическая постановка задачи формулируется так. Рассматривается односвязная ограниченная область О на плоскости или в пространстве, имеющая границу Го- — односвязная область с границей Гх, такая, что fii <Е П. Границы Г0 и Ti достаточно гладкие. Через обозначается область

По = ÍÍ \ ПТ.

В предположении, что граница Го известна и а0, 0\ — заданные положительные постоянные, требуется определить неизвестную границу Гх

и функции и>(М): и> G С(П),и>(М) = и{(М),(М) е Q¿(i = 0,1), где

и{ е С2(П0) П СЧЩХ е C\íli) П с1(й[),

Au¡{M) = О, (М) €^,» = 0,1;

(24)

4(М) = и{{М), {М)£ТХ, (25)

ди{{М) ди{(М)

(М)еГь (26)

4(М) = Р(М), (М) е Г0, (27)

~ (М) е Г0, (28)

где 3 = 1,2,...,тп. Здесь пара функций Р(М) и д^{М) соответствует j-мy измерению потенциала и его нормальной производной на внешней границе, а га — число измерений.

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

В §3.1.1 строится численный метод решения поставленной задачи. При каждом фиксированном у, функция р{х, у) берется в качестве начальных данных задачи Дирихле (24)-(27), и получается обратная задача Дирихле-Неймана с дополнительным условием (28), определяемым функцией д*(х, у). Для этой задачи по принципам, изложенным в §1.1, строится операторное уравнение. Для его построения вводится представление функции и^(х,у) в виде суммы двух потенциалов простого слоя. С использованием условий (26), (27), (28) и свойств потенциала простого слоя получается операторное уравнение для функции г(ф), которая параметризует границу Гх в полярной системе координат,

АЬ = д]. (29)

После построениея операторов Л7 для всех 3 = 1,2, ...,т обратная задача

Рисунок 3. Результаты вычислительного эксперимента для задачи с 1 измерением (а) и 3-мя измерениями (б). Кривые Г\ и, Г? — точный контур Г1 и начальное приближение. Г" и г^7 — кривые, полученные в результате 14 итераций в вычислительном эксперименте с 1 измерением (а) и в результате 17 итераций в вычислительном эксперименте с 3 измерениями.

(24)-(28) сводится к решению операторного уравнения

Атг = дт (30)

где оператор Ат = (^А1,А2,..., отображает функцию г(-ф) в вектор-

функцию дт = (^дг,д2, ■

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

В §3.1.2 приводятся сведения о программной реализации предложенного итерационного метода и результаты проведенных вычислительных экспериментов. На рис. 3 представлено сравнение

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

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

Результаты третьей главы были опубликованы в работах [3, 7, 8].

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

Основные результаты

1. Разработаны и программно реализованы численные методы определения неизвестной границы неоднородности в двумерной задаче электроимпедансной томографии при кусочно-постоянной проводимости в случае одного измерения на внешней границе;

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

3. Для двумерной и трехмерной задач элекгроимпедансной томографии с несколькими измерениями на внешней границе предложены и про1раммно реализованы итерационные методы определения границы неоднородности.

Список публикаций

[1] Гаврилов C.B., Денисов A.M. Численный метод определения границы неоднородности в задаче Дирихле для уравнения Лапласа в кусочно-однородной среде//Ж. вычисл. матем. и матем. физ. 2010. т.50. № 8. с. 1462 - 1470.

[2] Гаврилов C.B., Денисов А.М. Численные методы определения границы неоднородности в краевой задаче для уравнения Лапласа в кусочно-однородной среде//Ж. вычисл. матем. и матем. физ. 2011. т.51. № 8. с.1476 - 89.

[3] Гаврилов C.B., Денисов A.M. Итерационные методы определения границы неоднородности в краевой задаче для уравнения Лапласа в кусочно-однородной среде/А^ международная конференция "Математические идеи ПЛ.Чебышева и их приложение к современным проблемам естествознания Обнинск, 14-18 мая 2011 г. Тезисы докладов. с.71.

[4] Гаврилов C.B., Денисов А.М. Итерационный метод решения трехмерной задачи электроимпедансной томографии в случае кусочно-постоянной проводимости и одного измерения на границе//Ж. вычисл. матем. и матем. физ. 2012. т.52. № 8. с. 1426-36.

[5] Гаврилов C.B., Денисов A.M. Итерационный метод решения задачи

электроимпедансной томографии в случае кусочно-постоянной проводимости//Международная конференция, посвященная 80-летию со дня рождения академика М.М.Лаврентьева "Обратные и некорректные задачи математической физики Новосибирск, 5-12 августа 2012 г. Тезисы докладов, с. 186.

[6] Гаврилов C.B. Численный метод решения задачи электроимпедансной томографии в случае кусочно-постоянной проводимости и одного измерения на границе/УПрикладн. матем. и информ. 2012. т.41. с.38-47.

[7] Гаврилов C.B. Итерационный метод решения трехмерной задачи электроимпедансной томографии в случае кусочно-постоянной проводимости и нескольких измерений на границе//Вычисл. методы и программ. 2013. т.14. с.26-30.

[8] Гаврилов C.B. Численный метод решения трехмерной задачи электроимпедансной томографии в случае кусочно-постоянной проводимости и нескольких измерений на границе//4-я Международная конференция «Функциональные пространства. Дифференциальные операторы. Общая топология. Проблемы математического образования.», посвящённая 90-летию со дня рождения члена-корреспондента РАН, академика Европейской академии наук Л.Д. Кудрявцева, Москва, РУДН, 25-29 марта 2013 г. Тезисы докладов. с.403-04.

Цитированная литература

[9] Денисов А.М.,Захаров Е.В.,Калинин A.B..Калинин В.В. Численные методы решения некоторых обратных задач электрофизиологии

сердца//Дифференц. ур-ния. 2009. т.45. № 7. с.1014-1022.

[10] Тихонов А.Н. Об устойчивости обратных задач//Докл. АН СССР. 1943. Т.39 № 5. С 195-98.

[11] Тихонов А.Н., Арсенин В.Я. Методы решения некорректных задач.//Москва: Наука. 1974.

[12] Choi М.Н., Као T.-J., Isaacson D., Saulnier G.J., Newell J.C. A simplified model of mammography geometry for breast cancer imaging with electrical impedance tomography//Proceedings of the 26th Annual International Conference of the IEEE EMBS, San Francisco, CA, USA, September 2004, pp.960-963.

[13] Dihn Nho Hao, Lesnic D. The Cauchy problem for Laplace's equation via the conjugate gradient method//IMA Journal of Applied Mathematics (2000) 65, 199-217.

[14] Saeed S. Babaeizadeh, Dana H., Brooks D.H. Electrical impedance tomography for piecewise constant domains using boundary element shape-based inverse solutions/ЯЕЕЕ Trans Med Imaging 26(5) pp.637-47 (2007)

[15] Saulnier G.J., Blue R.S., Newell J.C., Isaacson D., Edic P.M. Electrical impedance tomography/ЯЕЕЕ Signal Processing Magazine. 2001. vol. 18, № 6, pp. 31-43.

[16] Zou Y„ Guo Z. A review of electrical impedance techniques for breast cancer detection//Medical Engineering and Physics. 2003. vol. 25, № 2, pp. 19-90.

Напечатано с готового оригинал-макета

Издательство ООО "МАКС Пресс" Лицензия ИД N 00510 от 01.12.99 г. Подписано в печать 02.09.2013 г. Формат 60x90 1/16. Усл.печл. 1,0. Тираж 100 экз. Заказ 263. Тел./факс: (495) 939-3890,939-3891. 119992, ГСП-2, Москва, Ленинские горы, МГУ им. М.В. Ломоносова, 2-й учебный корпус, 527 к.

Московский государственный университет имени М.В.Ломоносова

На правах рукописи

04201361424

Гаврилов Сергей Вадимович

Численные методы решения задачи электроимпедансной томографии в случае кусочно-постоянной проводимости

Специальность 05.13.18 - «Математическое моделирование, численные

методы и комплексы программ»

Диссертация на соискание ученой степени кандидата физико-математических наук

Научный руководитель д.ф.-м.н., профессор Денисов Александр Михайлович

Москва - 2013

Содержание

Введение 4

Обзор литературы 9

Глава 1. Численные методы решения двумерной задачи электроимпедансной томографии для кусочно-постоянной проводимости в случае одного измерения на границе 17

1.1 Метод решения обратной задачи Дирихле-Неймана..... 18

1.1.1 Операторное уравнение для неизвестной границы

в случае обратной задачи Дирихле-Неймана..... 18

1.1.2 Итерационный метод решения обратной задачи Дирихле-Неймана.................... 23

1.2 Метод решения обратной задачи Неймана-Дирихле..... 28

1.2.1 Операторное уравнение для неизвестной границы

в случае обратной задачи Неймана-Дирихле..... 28

1.2.2 Итерационный метод решения обратной задачи Неймана-Дирихле.................... 32

1.3 Метод, основанный на введении вспомогательного контура . 37

1.3.1 Продолжение значений потенциала на вспомогательный контур........................... 37

1.3.2 Численный метод решения задачи с данными на вспомогательном контуре................ 41

1.4 Программная реализация численных методов и результаты вычислительных экспериментов................ 51

Глава 2. Численный метод решения трехмерной задачи электроимпедансной томографии для кусочно-постоянной проводимости в случае одного измерения на границе 60

2.1 Операторное уравнение для неизвестной поверхности . ... 61

2.2 Итерационный метод решения обратной задачи....... 66

2.3 Программная реализация численного метода и результаты вычислительных экспериментов................ 72

Глава 3. Численные методы решения задачи электроимпедансной томографии для кусочно-постоянной проводимости в случае нескольких измерений на границе 79

3.1 Двумерная задача......................... 79

3.1.1 Численный метод решения............... 80

3.1.2 Программная реализация численного метода и результаты вычислительных экспериментов ..... 84

3.2 Трехмерная задача........................ 86

3.2.1 Численный метод решения............... 87

3.2.2 Программная реализация численного метода и результаты вычислительных экспериментов ..... 92

Заключение 99

Литература 100

Введение

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

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

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

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

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

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

Цель и задачи работы

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

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

3. Программная реализация и проведение вычислительных экспериментов для предложенных методов.

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

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

Апробация работы. Основные результаты работы докладывались на следующих конференциях.

1. V международная конференция "Математические идеи П.Л.Чебышева и их приложение к современным проблемам естествознания"(Обнинск, 14-18 мая 2011 года).

2. Международная конференция, посвященная 80-летию со дня рождения академика М.М.Лаврентьева "Обратные и некорректные задачи математической физики"(Новосибирск, 5-12 августа 2012 года).

3. Научная конференция "Тихоновские чтения"(Москва, МГУ имени М.В.Ломоносова, 29.10.2012-02.11.2012).

4. 4-я Международная конференция «Функциональные пространства. Дифференциальные операторы. Общая топология. Проблемы математического образования.», посвящённая 90-летию со дня рождения члена-корреспондента РАН, академика Европейской академии наук Л.Д. Кудрявцева (Москва, РУДН, 25-29 марта 2013 года).

Публикации. Материалы диссертации опубликованы в 8 печатных работах, из них 4 статьи в рецензируемых журналах, рекомендованных ВАК [7, 8, 10, 13] статья [12] и 3 тезиса докладов на конференциях [9, 11, 14].

Структура диссертации. Диссертация состоит из введения, обзора литературы, 3 глав, заключения и библиографии. Общий объем диссертации 110 страниц, включая 20 рисунков. Библиография включает 94 наименований на 11 страницах.

Обзор литературы

Томографические методы в настоящее время широко применяются в различных областях науки, техники и медицины. Началом их возникновения можно считать открытие в 1895 году В.К.Рентгеном ионизирующего излучения, которое впоследствии было названо рентгеновским. Это открытие привлекло значительный интерес научного сообщества и уже с 1896 года стали проводиться исследования этого излучения и разрабатываться технологии его практического применения. В течение двух последующих десятилетий развитие рентгенологических методов привело к формированию нового направления в медицине -радиологии [25]. Помимо развития рентгенографии первая половина 20 века характеризуется появлением устройств, использующих ультразвуковые волны. В 1940 году был создан первый коммерческий аппарат для ультразвуковой дефектоскопии [4]. Первая половина 20 века примечательна еще одним важным открытием, давшим основу для возникновения магнитно-резонансной томографии: в 1938 г. И.Раби открыл явление ядерного магнитного резонанса в молекулярных пучках [1], а в 1946 г. Ф.Блох и Э.М.Парселл параллельно открыли ядерный магнитный резонанс в жидкостях и твердых телах [28].

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

проведены первые испытания компьютерного томографа [38]. Спустя год в 1973 г. П.Лотербур опубликовал первые снимки, полученные с помощью метода магнитно-резонансной томографии [74]. С тех пор появилось множество новых разновидностей томографических методов, существенный прогресс достигнут в совершенствовании технологий и повышении их эффективности, значительно расширилась область применения томографических методов. С современным состоянием томографии можно ознакомиться, например, по книгам [26, 59, 34, 88].

Одним из видов томографии, активно развивающихся в настоящее время, является электроимпедансная томография. Этот метод позволяет получать представление об электрических характеристиках вещества внутри исследуемого объекта по одновременным измерениям значений потенциала и электрического тока на его поверхности [43]. В качестве этих электрических характеристик выступают электрическая проводимость и диэлектрическая проницаемость [24]. Определение этих характеристик представляет большой практический интерес, поскольку знание распределения проводимости или диэлектрической проницаемости внутри объекта позволяет делать выводы о его строении.

Различные виды горных пород, минералов и полезных ископаемых характеризуются разной электропроводностью [70], что делает возможным использование электроимпедансной томографии в геофизических исследованиях. Примером служат обнаружение месторасположения полезных ископаемых [53, 81], исследования уровня загрязненности и поиск источников загрязнения в почве [82, 83, 51].

Различные виды биологических тканей или одинаковые ткани в разном физиологическом состоянии также имеют разные значения проводимости и диэлектрической проницаемости [35, 87].

Знание распределения этих электрических характеристик позволяет идентифицировать расположение тканей и получать медицинские снимки. Медицинскими приложениями электроимпедансной томографии являются мониторинг физиологической активности сердца и легких [77], обнаружение жидкости в брюшной полости [85], онкологическая диагностика [41], визуализация при проведении абляции злокачественных опухолей [65, 80, 93] и ряд других приложений [57, 75, 44, 58, 60].

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

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

Пусть П С К1, с1 > 2 ограниченная область с границей Г и во всех точках области О определена функция ^{х, и), представляющая собой значение комплексной проводимости: 7(2,0;) = <т(х) + ше(х), где и — частота электромагнитного поля. Комплексная проводимость или адмиттанс — физическая величина, объединяющая свойства электрической проводимости а и диэлектрической проницаемости е [5].

Для случая изотропной среды и электромагнитного поля фиксированной частоты при определенных ограничениях [48] из уравнений Максвелла можно получить уравнение

Ч(-у{х,и)Чф(х,ш)) = 0, я е Q, (1)

где ф(х,ш) — скалярный потенциал. Условие Дирихле для этого уравнения выглядит следующим образом

<j>{x,J) = V(x,u), хет. (2)

Условие Неймана выглядит так

7(ж,ц) 1 =/(з?,ц), хеТ. (3)

Здесь п — нормаль к границе Г, а 1(х,и) — функция, представляющая распределение электрического тока на внешней границе. Для уравнений (1)-(3) ставится следующая задача. Требуется определить неизвестную функцию 7(х,и>), исходя из знания оператора Дирихле-Неймана

АуУ^х.сй) = 7 (х,ш)———,

для всех функций V(x, и). Здесь ф(х, со) — решение краевой задачи (1),(2).

Впервые приведенная математическая задача была рассмотрена в работе А.Кальдерона [46]. Последние десятилетия велось активное теоретическое исследование обратных задач такого типа, для широкого класса функций доказана единственность решения в двумерном случае [76] и в трехмерном пространстве [90, 78, 79].

Некоторые результаты достигнуты в разработке численных методов

решения математической задачи, сформулированной А.Кальдероном. Первые алгоритмы были основаны на линеаризации уравнения для решения обратной задачи. При этом рассматривался статический случай {и = 0) и делалось предположение о близости проводимости к постоянной величине, в окрестности которой исходное уравнение линеаризовывалось. Методы основанные на этом принципе излагаются в работах [36, 39].

Другой подход к численному решению задачи Кальдерона заключался в последовательном восстановлении проводимости вглубь исследуемой области. Значение проводимости на внешней границе считалось известным, что является оправданным допущением, поскольку её измерение в большинстве приложений не представляет сложности. Методы с применением послойного восстановления проводимости рассматриваются в работах [50, 62, 89, 91].

Эффективный метод численного решения был разработан для задачи, в которой предполагалось, что проводимость исследуемой области отличается от постоянной величины лишь в конечном числе подобластей, в которых проводимость являлась непрерывной функцией, отличной от постоянной величины. Для этого случая был разработан алгоритм, который позволял определять для фиксированной точки области, принадлежит ли она одному из вкраплений, или находится в области, где проводимость постоянна. Критерий принадлежности был основан на свойствах оператора Неймана-Дирихле (Ла)-1, который использовался в качестве исходной информации. Методы, основанные на этом подходе рассматривались в работах [45, 67].

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

свод