Применение метода регуляризации Тихонова для численного решения обратной задачи электрокардиографии Текст научной статьи по специальности «Математика»

УДК 517.958, 519.632

А.М. Денисов, Е.В. Захаров, А.В. Калинин, В.В. Калинин

ПРИМЕНЕНИЕ МЕТОДА РЕГУЛЯРИЗАЦИИ ТИХОНОВА ДЛЯ ЧИСЛЕННОГО РЕШЕНИЯ ОБРАТНОЙ ЗАДАЧИ ЭЛЕКТРОКАРДИОГРАФИИ1

(кафедра математической физики факультета ВМиК, НЦССХ им.А.Н. Бакулева,

e-mail: [email protected])

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

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

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

При формулировке обратной задачи электрокардиографии используется следующая биофизическая модель [3, 4]. Грудная клетка рассматривается как проводник второго рода с постоянным коэффициентом электропроводности, занимающий ограниченную область пространства и окруженный диэлектрической средой — воздухом. Электрическое поле сердца описывается моделями электродинамики стационарных токов, т. е. его потенциал удовлетворяет уравнению Пуассона. Источники электрического поля расположены в толще сердечной мышцы. В области пространства, ограниченного внешней поверхностью сердца и поверхностью тела, источники поля отсутствуют и потенциал поля удовлетворяет уравнению Лапласа. На поверхности грудной клетки, контактирующей с воздухом (граница "проводник-диэлектрик"), нормальная составляющая напряженности электрического поля сердца равна нулю. На этой же части поверхности известен потенциал электрического поля сердца, полученный в результате экспериментальных измерений. Требуется найти потенциал электрического поля на внешней поверхности сердца.

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

2. Математическая постановка задачи. Пусть <> — область в пространстве НЛ. ограниченная снаружи замкнутой поверхностью Si, а изнутри замкнутой поверхностью ¿>2- Поверхности Si и S2

1 Работа выполнена при поддержке Российского фонда фундаментальных исследований, проект № 08-01-00314.

предполагаются достаточно гладкими. Требуется найти функцию и(х,у,г), такую, что

(ж,у,.г)еП, (1)

1, (2) (ж,у,г)€бг1. (3)

Эта задача называется задачей Коши для уравнения Лапласа. Задача Коши для уравнения Лапласа является классическим примером некорректно поставленной задачи. Исследованию ее единственности и условной устойчивости посвящено большое число работ (см., например, [5] и цитированную там литературу). Для численного решения этой задачи был разработан ряд регуляризирующих алгоритмов, основанных на различных принципах [6-9].

В данной работе рассматривается численный алгоритм решения задачи Коши для уравнения Лапласа (1)-(3) в трехмерной области, основанный на методе регуляризации Тихонова [10].

Аи(х, у, г) = О,

и(х,у,г) = 11(х,у,г), ди

— (ж, у, г) = 0,

3. Метод регуляризации Тихонова для решения обратной задачи электрокардиографии. Рассмотрим смешанную краевую задачу

Аи(х, у, ¿) = 0, (ж, у, ¿) € О, (4)

и(х,у,г) = «(ж, у, г), (ж, у, г) € (5) ди

— (ж,у,г) = 0, (ж,у,г)€бг1. (6)

Обозначим решение этой задачи для заданной функции г>(ж, у, г), (ж, у, ¿) € <?2) через и(х, у, г; у). Задача (4)-(6) определяет оператор А, отображающий потенциал V на поверхности Б2 в потенциал на поверхности

Ау = и(х, у, г; г>), (ж, у, г) € <51-

Тогда задача (1)-(3) может быть сформулирована как задача решения операторного уравнения первого рода:

Ау = 11(х,у,г), (ж, у, г) € ¿>1. (7)

Применим для ее решения метод регуляризации Тихонова [10]. Пусть для точных значений и(х,у,г), (ж, у, ¿) € ¿»1, существует точное решение уравнения (7) д(х, у, г), (ж, у, ¿) € <?2) но С/(ж, у, ¿) неизвестна, а задано ее приближение [/¿(ж, у, г), (ж, у, г) € ¿»1, такое, что ЦС/д — ¿^Ц^^ ^ Требуется, зная [/¿(х, у, г) и величину погрешности 8, построить приближенное решение «¿(ж, у, г). Рассмотрим функционал

М>] = - + а , а > 0. (8)

Приближенное решение уа определяется как элемент, реализующий минимум функционала Ма[и], в котором параметр регуляризации а должным образом зависит от величины погрешности 8, т.е. а = ск(<5). Приближенное решение уа может быть найдено из уравнения

ау + А*Ау = А*Щ, (9)

представляющего собой необходимое условие минимума функционала (8). Оператор А* отображает функцию д(х,у,г), заданную на поверхности Бг, в функцию, заданную на поверхности ¿>2) и определяется следующей смешанной задачей:

Дго(ж, у, г) = 0, (ж,у, г) £ !1, (10)

и)(х,у,г)=0, (ж, у, г) € ¿»25 (11)

ди;

— (х,у,г)=д(х,у,г), (ж, у, г) € Яь (12)

Обозначив решение этой задачи для заданной функции д через получим [11, 12], что

Таким образом, операторы А и А*, входящие в уравнение (9), определяются смешанными задачами (4)-(6) и (10)—(12) соответственно.

Значение параметра регуляризации а = а(5) может быть найдено из принципа невязки [13]

\\Ауа-и8\\Ы81) = 6.

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

решения, вместо функционала (8) можно ввести функционал

М(

>;£] = \\Ау - и6\\12{31) + а\\у - у\\12{32) , а > 0.

Необходимое условие минимума, из которого находится приближенное решение, имеет в этом случае вид

а{у-Ъ) + А*АУ = А*и6. (13)

Из изложенного выше следует, что задача построения приближенного решения сводится к задаче решения уравнения (9) или уравнения (13). Перейдем к методам и алгоритмам численного решения этих уравнений и краевых задач (4)-(6), (10)—(12).

4. Метод граничных интегральных уравнений. Как было показано в предыдущих разделах, основные вычислительные процедуры сводятся к нахождению решений уравнений Эйлера (9) или (13). Для этого необходимо решать смешанные краевые задачи (4)-(6) и (10)—(12). Среди методов численного решения задач данного типа в вычислительной практике обычно применяются метод конечных элементов (дискретизация всей ЗБ-области) и метод граничных интегральных уравнений (см., например, [14, 15]).

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

Для перехода от краевых задач к граничным интегральным уравнениям обычно применяют два подхода. В первом решение ищется в виде потенциала простого или двойного слоя, во втором неизвестные краевые условия находятся как решения интегральных уравнений, непосредственно следующих из основной формулы Грина (третьей формулы Грина). Смешанные краевые задачи типа (4)-(6) и (10)—(12) удобнее решать вторым методом.

Рассмотрим область О (рис. 1) с достаточно гладкой границей <90 = = и ¿?2. Потенциал и в О (включая границу) удовлетворяет уравнению Лапласа. Для точек границы <90 с учетом фундаментального решения уравнения Лапласа можно записать интегральное представление решения:

2тт-и(М) =

или

2тг • и(М) +

1

Г

где М, Р Е 50, г

евклидово расстояние между точками МиР, 1 /г

Рис. 1. Область

(14)

фундаментальное решение

уравнения Лапласа в й3, Пр — единичный вектор внутренней нормали в точке Р, д(Р) =

Разобьем поверхность <ЭП на граничные элементы dsf. S = ds\ U ds2 U ... U dsn. Введем систему из n линейно независимых базисных элементов cpi, ср2, ■ ■ ■, <рп (характеристических функций), определенных следующим образом:

'<Pi(s) = 1, s G dsi, jpi(s) = 0, s £ dsi.

Значение потенциала и его нормальной производной представляются в виде разложения по системе базисных функций cpi (кусочно-постоянная аппроксимация):

г= 1

п

i(s) = fa' cPi(s)-'

(15)

1=1

где щ — значение «(«), а Д — значение в центре тяжести г-го граничного элемента. Подставляя (15) в (14), получим

Ни =

где элементы матриц Н и О определяются так:

е(1/г„)

5,

(16)

<-окЫ1м±й

J апр ■

hij —

I

dil/rj,)

дп

^ds

гфз-,

2тг, i = //'.

9ij =

ds.

Поскольку граница 5 разбита на части 51 и 52 (рис. 1), на каждой из которых задаются граничные условия, то уравнение (16) удобнее представить в виде

г* 1 Г„

(17)

'ЯЦ Н\2 U\ G н G12 Ql

JI-JI IITJ и2_ р21 G 22 Я2.

Н = [///,-/]• О = [//*•/]•

где индекс к обозначает номер поверхности с фиксированной точкой Л /. а индекс I — номер поверхности с переменной точкой Р.

При решении смешанных краевых задач (4)-(6) и (10)—(12) заданы граничные условия (С^,^). Таким образом, систему (17) с учетом известных условий можно переписать в виде

(18)

Из представления (18) можно получить в явном виде матричные уравнения для решения смешанной краевой задачи, т. е. нахождения и д2- Поскольку смешанная краевая задача является корректно поставленной по Адамару, то матрицы данных уравнений являются достаточно хорошо обусловленными и относительно легко обращаются.

~Нц Н\2 U\ GH G12 'Qi

JI-JI IITJ р2\ G 22 Я2 „

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

• задавалась расчетная область П;

• задавалась система источников электрического поля;

• решалась задача вычисления поля заданных источников на поверхностях 52 и ¿>! по известным соотношениям теории потенциала;

• на значение потенциала на поверхности 1 накладывался гауссовский шум величиной от 1 до 15%;

• решалась обратная задача вычисления потенциала на поверхности Б2 по зашумленным данным;

• полученное решение сравнивалось с эталонным.

В качестве расчетной области использовалась геометрическая конфигурация, приведенная на рис. 1. Внутри поверхности $2 задавались источники дипольного и квадрупольного типов, а также случайно распределенные в пространстве единичные заряды разных знаков. Система зарядов формировалась таким образом, чтобы конфигурация электрического поля на поверхности $2 моделировала электрическое поле сердца. Триангуляция поверхности проводилась на основе алгоритмов, предложенных в работе [17]. Число граничных элементов дискретизации поверхности составляло порядка 600-1500 треугольников.

Далее итерационным способом решались уравнения (9) и (13). На каждом шаге итерации методом граничных интегральных уравнений решались прямые смешанные краевые задачи (4)-(6) и (10)—(12) и вычислялись значения операторов А и А*. Параметр регуляризации а подбирался по принципу невязки.

Для всех конфигураций заданных источников метод устойчиво сходился к решению с уровнем погрешности от 2 до 15% при уровне шума 1%. Точность решения в рамках используемых моделей источников не зависела от характера восстанавливаемого потенциала (с учетом ограничений метода граничных элементов на близость источников к границе).

Рис. 2. Конфигурация из трех распределенных по объему источников поля: заданная конфигурация (а); восстановленная конфигурация (б)

Качество решения можно оценить из рис. 2, 3. На этих рисунках изображено распределение потенциала поля на верхней части внутреннего эллипсоида (см. поверхность $2 на рис. 1) как функция и = /(ж, у); ж, у принадлежат верхней части $2.

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

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

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

-50

а

б

-1 -5

0

4

3

2

5

1

а

б

Рис. 3. Конфигурация из трех источников, сгруппированных в верхней части внутренней поверхности: заданная конфигурация (а); восстановленная конфигурация (б)

1. Барр Р., Спэк М. Решения обратной задачи, выраженные непосредственно в форме потенциала // Теоретические основы электрокардиологии. М.: Медицина, 1979.

2. Бокерия J1.A., Шакин В. В., Мир с кий Г. В., П о л я к о в а И. П. Численные методы электрофизиологической оценки состояния сердца. М.: ВЦ АН СССР, 1987.

3. MacLeod R. S., Brooks D.H. Recent progress in inverse problems in electrocardiology // IEEE Eng. in Med. Bio. Mag. 1998. 17. N 1. P. 73-83.

4. Ramanathan C., Ghanem R. N., Jia P., Kyungmoo R., Rudy Y. Noninvasive electrocardiographic imaging for cardiac electrophysiology and arrhythmia // Nature Medicine. 2004. 10. P. 422-428.

5. Лаврентьев M. M., Романов В. Г., Шишатский С. П. Некорректные задачи математической физики и анализа. М.: Наука, 1980.

6. Лат тес Р., Лионе Ж.-Л. Метод квазиобращения и его приложения. М.: Мир, 1970.

7. Козлов В. А., Мазь я В. Г., Фомин А. В. Об одном итерационном методе решения задачи Коши для эллиптических уравнений // ЖВМиМФ. 1991. 31. № 1. С. 64-74.

8. Гласко В. Б. Обратные задачи математической физики. М.: Изд-во МГУ, 1984.

9. Вабищевич П.Н., Гласко В.Б., Криксин Ю.А. О решении одной задачи Адамара с помощью ре-гуляризирующего по А.Н. Тихонову алгоритма // ЖВМиМФ. 1979. 19. № 6. С. 1462-1473.

10. Тихонов А.Н., Арсенин А. Я. Методы решения некорректных задач. М.: Наука, 1986.

11. Kabanikhin S.I.,Karchevsky A.L. Optimizational method for solving the Cauchy problem for an elliptic equation //J. Inverse Ill-Posed Problems. 1995. 3. P. 21-46.

12. Dihn Nho Hao, Lesnic D. The Cauchy problem for Laplace's equation via the conjugate gradient method // IMA J. of Appl. Math. 2000. 65. P. 199-217.

13. Морозов В. А. О принципе невязки при решении операторных уравнений методом регуляризации // ЖВМиМФ. 1968. 8. № 2. С. 295-309.

14. Бреббия К., Теллес Ж., Вроубел Л. Методы граничных элементов. М.: Мир, 1987.

15. Р ozr ikidis С. A practical guide to boundary elements methods with the software library BEMLIB. Chapman & Hall/CRC, 2002.

16. Тихонов A.H., Самарский А. А. Уравнения математической физики. M.: Наука, 2004.

17. Гольник Э.Р., Вдовиченко А. А., Успехов А. А. Построение и применение препроцессора генерации, управления качеством и оптимизации сеток триангуляции контактных систем / / Информационные технологии. № 4. Воронеж: Изд-во ВГТУ, 2004.

СПИСОК ЛИТЕРАТУРЫ

Поступила в редакцию 04.09.07