пїњ

–ј«Ќќ—“Ќџ≈ —’≈ћџ Ќј ќ—Ќќ¬≈ ћ≈“ќƒј  ќЌ≈„Ќџ’ ќЅЏ®ћќ¬ ƒЋя «јƒј„» ЁЋ≈ “–ќ»ћѕ≈ƒјЌ—Ќќ… “ќћќ√–ј‘»»

¬≈—“Ќ»  “ќћ— ќ√ќ √ќ—”ƒј–—“¬≈ЌЌќ√ќ ”Ќ»¬≈–—»“≈“ј
2014 ћатематика и механика є 3(29)
”ƒ  519.6
≈.—. Ўерина, ј.¬. —тарченко
–ј«Ќќ—“Ќџ≈ —’≈ћџ Ќј ќ—Ќќ¬≈ ћ≈“ќƒј  ќЌ≈„Ќџ’ ќЅЏ®ћќ¬ ƒЋя «јƒј„» ЁЋ≈ “–ќ»ћѕ≈ƒјЌ—Ќќ… “ќћќ√–ј‘»»1
ѕолучены разностные схемы дл€ решени€ задачи электроимпедансной томографии. ¬ывод схем осуществл€етс€ с помощью метода конечных объЄмов на неструктурированных сетках. „исленное сравнение дл€ конечных объЄмов разной формы выполнено на двух тестовых задачах с аналитическим решением. ѕолученные результаты также сравниваютс€ с решением по методу конечных элементов.
 лючевые слова: электроимпедансна€ томографи€, метод конечных объЄмов, метод конечных элементов, разностные схемы, неструктурированные сетки.
1. ¬ведение
Ёлектроимпедансна€ томографи€ (Ё»“) - это быстро развивающийс€ неинвазивный метод визуализации со спектром применени€, включающим медицинскую томографию, неразрушающий контроль материалов и контроль промышленных процессов [1-5]. ћетод Ё»“ восстанавливает неизвестное распределение или оценивает несколько параметров электрической проводимости (или полного импеданса) внутри объекта, если измерено электрическое напр€жение или сила тока на его границе. ѕо реконструкции даЄтс€ оценка внутренней структуры объекта. Ё»“ позвол€ет наблюдать динамические процессы и работать со статическими изображени€ми.
ќтмечают несколько преимуществ Ё»“ перед другими видами томографии: компактный размер оборудовани€, простота устройства, невысока€ стоимость и быстрый сбор данных.  роме того, в медицинских приложени€х ценитс€ безопасность дл€ биологических тканей. Ќаправление Ё»“ представл€етс€ перспективным в разных приложени€х, но его широкое внедрение ограничено невысоким качеством получаемых изображений. Ќаравне с необходимостью усовершенствовани€ аппаратных средств томографии следует удел€ть внимание разработке математической теории Ё»“ и новых численных методов решени€ такой задачи.
¬осстановление проводимости с помощью Ё»“ €вл€етс€ нелинейной и некорректной обратной коэффициентной задачей. —ложность задачи вызвана нарушением услови€ устойчивости в классическом определении корректности по јдама-ру [2].  ак следствие, реконструкци€ статических изображений Ё»“ очень чувствительна к ошибкам измерений прибора и погрешности аппроксимации дифференциальной задачи разностной, а также качеству физического моделировани€ рассматриваемого объекта. ƒанное исследование направлено на изучение особенностей численного моделировани€ пр€мой задачи Ё»“, €вл€ющейс€ важной компонентой обратной задачи. ќтдельное внимание уделено вопросу уменьшени€ погрешности аппроксимации, котора€ по€вл€етс€ в ходе дискретизации задачи. ¬се
1 –абота выполнена при финансовой поддержке ћинобрнауки –‘ в рамках государственного задани€ є 2014/223 (код проекта 2382).
результаты получены дл€ двумерной статической задачи и при необходимости могут быть перенесены на трЄхмерный случай или использоватьс€ в динамической реконструкции.
¬ Ё»“ на поверхность объекта прикрепл€етс€ набор электродов (рис. 1, а), затем выполн€етс€ сери€ измерений, в которых через объект пропускаетс€ электрический ток небольшой силы и высокой частоты. ќбычно, сила тока выбираетс€ в интервале от 1 до 5 мј с частотой от 10 до 150 к√ц. —хема проведени€ эксперимента - схема подведени€ тока и сбора данных - определ€етс€ исследователем. ћежду электродами граница объекта контактирует с воздухом. –езультирующее напр€жение измер€етс€ на электродах. ¬с€ сери€ данных Ђэлектрический ток -напр€жениеї используетс€ дл€ реконструкции неизвестного распределени€ электрической проводимости внутри объекта.
¬ данной статье использована модель Ё»“ с идеально провод€щими электродами, котора€ называетс€ в научной литературе gap model [6].
ѕредполагаетс€, что учитываемые моделью физические процессы описываютс€ электрическим и магнитным пол€ми, которые удовлетвор€ют стационарным уравнени€м ћаксвелла. ¬ данной работе рассматриваетс€ случай применени€ Ё»“ к объектам биологической природы. “кани провод€т электричество, потому что они содержат ионы, которые перенос€т электрический зар€д. ѕроводимость биологических тканей обладает определенной спецификой и зависит от частоты тока.  огда электромагнитные волны проход€т через биологический объект, интенсивность магнитного пол€ становитс€ пренебрежимо малой [1, 2], поэтому магнитна€ компонента в уравнени€х ћаксвелла не рассматриваетс€. ¬ этом случае система уравнений ћаксвелла может быть с использованием закона ќма сведена к уравнению эллиптического типа в частных производных. ¬ силу отсутстви€ внутри объекта источников электрического пол€, а также с учЄтом 1-го правила  ирхгофа дл€ электрической цепи, сумма втекающих и вытекающих токов равн€етс€ нулю - выполн€етс€ закон сохранени€ зар€да. “акже предполагаетс€, что на электроде, где инжектируетс€ ток, известно значение потенциала электрического тока.
— учЄтом прин€тых условий математическа€ модель Ё»“ дл€ двумерного случа€ (рис. 1, а) может быть представлена краевой задачей [2, 7]:
2. ѕостановка задачи Ё»“
2.1. ‘изическа€ постановка задачи
2.2. ћатематическа€ постановка задачи
V-(cVcp) = 0, (х, y )eQ ;
(1)
(2)
dn k=1
(3)
где о(х,у) - электрическа€ проводимость, ф(х,у) - потенциал электрического пол€, п - единичный вектор внешней нормали, Jk - плотность электрического тока на
k-м электроде, ek - k-й электрод, L - число электродов. ‘ункци€ электрического потенциала pe H1 (Q), плотность тока j =ст- Vp е L2 (Q).
ѕо прин€той в литературе терминологии задача (1) - (3) с функцией электрического потенциала ф(х,у) в качестве неизвестного будет в дальнейшем называтьс€ пр€мой задачей Ё»“. ќбратной задачей Ё»“ именуетс€ задача (1) - (3) по поиску неизвестной функции проводимости о(х,у) в области Q, дополненна€ данными измерений электрического напр€жени€ на электродах.
3. ѕостроение численных схем дл€ задачи Ё»“
„увствительность реконструкции Ё»“ к ошибкам измерений и погрешности аппроксимации дифференциальных операторов ограничивает круг подход€щих численных методов, которыми можно с хорошим качеством решать обратную задачу.  ак правило, дифференциальные задачи привод€т к дискретному аналогу каким-либо сеточным методом. —реди них можно выделить следующие группы: конечно-разностные (ћ –), конечно-элементные (ћ Ё), гранично-элементные (ћ√Ё), конечно-объЄмные (ћ ќ) и другие. «адачу Ё»“ изначально ассоциировали с электродинамикой и моделированием цепей, поэтому был опробован, а впоследствии и наиболее распространилс€ подход по ћ Ё. ћ Ё обладает гибкостью, когда необходимо быстро улучшить качество аппроксимации без сложных изменений алгоритма. Ќо в [8] отмечаетс€, что ћ Ё не удовлетвор€ет условию непрерывности электрического тока, тем самым не соблюдаетс€ закон сохранени€ на элементе сетки и области в целом. — ћ Ё ищетс€ решение не конкретной дифференциальной задачи, а обобщЄнное решение вариационной задачи, выведенной из исходной и записанной в интегральной форме. ¬ данной работе привод€тс€ некоторые результаты построени€ конечно-объЄмных схем дл€ пр€мой задачи Ё»“, которые сравниваютс€ с решением ћ Ё по нескольким показател€м.
3.1. ƒискретизаци€ расчЄтной области
Ќа начальном этапе решени€ задачи выполн€етс€ дискретизаци€ расчЄтной области rah, и определ€етс€ набор сеточных функций электрического потенциала U и проводимости ch. —оздание расчЄтных сеток - отдельна€ область исследовани€. ¬ насто€щее врем€ доступны свободно распростран€емые и коммерческие программные продукты дл€ генерации сеток (GMSH, GAMBIT, QMG, Netgen, Tgrid, FEMLAB и др.).
¬ данном исследовании численные эксперименты проводились на треугольных неструктурированных сетках (M €чеек, N узлов), построенных в программах GAMBIT и GMSH. ѕример расчЄтной сетки приведЄн на рис. 1, б. «аметим, что сетка детализована вдоль границы области, потому что градиент потенциала быстро измен€етс€ вблизи контактов электродов. ƒл€ оценки качества построенных сеток с треугольными €чейками обычно вводитс€ какой-нибудь числовой показатель. „тобы осуществл€ть контроль качества построенных сеток была выбрана группа параметров: коэффициент деформации элемента сетки, угловой коэффициент, соотношение сторон треугольного элемента, соотношение радиуса вписанной к радиусу описанной окружности дл€ элемента сетки.  ак правило, точность расчЄтов снижаетс€, если элементы треугольной сетки сильно отличаютс€ от правильных треугольников. —ледует заметить, что решающим фактором при выборе сетки будет не только набор показателей качества и размер сетки, но и еЄ способ-
ность обеспечивать работоспособность конкретного примен€емого численного метода в целом. ƒолжна быть обеспечена вычислительна€ экономичность расчЄтов и достаточна€ точность приближенного решени€ полученных сеточных уравнений, поэтому, по-видимому, следует говорить о том, насколько хорошо сетка подходит дл€ примен€емого численного метода.
Д е5 ^
С'6^'"

( \е
о ≤
™≥к /™≥5
Ч--'^™14
–ис. 1. ƒвумерна€ модель области реконструкции ќ (а), на границе сё обозначены электроды вк, к = 1,...,16; расчЄтна€ сетка со сгущением к электродам (б), размер сетки -ћ = 4168 элементов, N = 2181 узел; расчЄтна€ сетка дл€ тестовой задачи с точечными электродами 1 и 13 (в), размер сетки - ћ = 136 элементов, N = 81 узел [8]
10
0
3.2. ¬арианты конечных объЄмов
 огда дифференциальна€ задача с помощью сеточного метода преобразуетс€ в разностную, то выбираетс€ способ соответстви€, по которому неизвестные значени€ сеточной функции св€зываютс€ с расчЄтной сеткой. ‘ункции электрического потенциала ф(х,х) и проводимости а(х,_у) замен€ютс€ дискретными аналогами, т.е. сеточными функци€ми ык и соответственно, ык е иъ,ак е!к, ÷ и - гильбертовы пространства. Ќа триангул€ции области обычно рассматриваютс€ два варианта - значение сеточных функций определ€етс€ в некоторой точке внутри €чейки сетки (рис. 2, а) или в узле сетки (рис. 2, б, в).  онечный объЄм в обоих случа€х строитс€ вокруг точки, с которой св€зываютс€ значени€ сеточной функции. ¬ насто€щей работе обсуждаетс€ реализаци€ ћ ќ на трЄх конечных объЄмах: треугольный конечный объЄм (рис. 2, а), барицентрический конечный объЄм (рис. 2, б) и конечный объЄм - €чейка ƒирихле - ¬ороного (рис. 2, в).
–ис. 2. ‘рагменты расчЄтной сетки: а - треугольный конечный объЄм с центром в точке –0; б - барицентрический конечный объЄм, построенный вокруг вершины –0; в - конечный объЄм - €чейка ƒирихле - ¬ороного, - построенный вокруг вершины –0
Ѕарицентрическа€ €чейка строитс€ вокруг вершины –0 по точкам пересечени€ медиан и точкам середин сторон треугольников. ¬ €чейке ƒирихле - ¬ороного вместо барицентров выбираютс€ точки пересечени€ срединных перпендикул€ров. “реугольные конечные объЄмы соответствуют €чейкам сетки, в них положение точки –0 задаетс€ в центре масс, точке пересечени€ срединных перпендикул€ров или биссектрис.
3.3.¬ывод численных схем
—огласно принципу ћ ќ дл€ каждого конечного объЄма записываетс€ ин-тегро-интерпол€ционный баланс. ”равнение (1) интегрируетс€ по конечному объЄму ќт,
[ ”-(а”ф)0 = 0. (4)
ѕт
ѕо формуле √рина уравнение (4) преобразуетс€ к равенству нулю криволинейного интеграла по границе от плотности тока
аЧ<» = 0. (5)
до дп
д0т
ƒальнейша€ оценка интеграла в левой части (5) зависит от формы конечного объЄма и способа аппроксимации градиента потенциала в подынтегральном выражении. ‘ункци€ проводимости задаетс€ кусочно-посто€нной по треугольным элементам сетки, 0 < с < о(х,у) < —, с и — - константы.
Ќа треугольных конечных объЄмах (рис. 2, а) интеграл в правой части (5) вычисл€етс€ по сторонам текущего треугольника ќт по интерпол€ционной формуле, в данном случае - средних пр€моугольников. ќпробовано несколько вариантов расположени€ точки –0 внутри конечного объЄма, с которой ассоциировались сеточные функции проводимости ск и электрического потенциала ик. ƒискретные значени€ определ€лись в точке пересечени€ медиан, срединных перпендикул€ров или биссектрис треугольника. ѕроизводные аппроксимировались в серединах сторон конечного объЄма, причЄм оценка выполн€лась с учЄтом условий сопр€жени€ на границе двух смежных €чеек, например дл€ грани у (рис. 2, а):
а^”0 ’^‘ ~
дп )  V дп
(‘^ =(‘)^
чтобы схема стала чувствительной к изменению свойств токопровод€щей среды. ѕотоки через границу конечного объЄма могут быть аппроксимированы разными соотношени€ми. ќдин из предложенных подходов рассматривает аппроксимацию
(дф^–∞ и^ ир Ћ дфЋр ћр и^
вида I Ч I Ђ Ч1------ и I Ч I и Ч1-----L, где Ћп–|) - длина отрезка перпенди-
V дп )^ ƒп–() V дп )щ ƒп–1
кул€ра из точки –0 до стороны у, јп–; - от точки –1 до стороны //', а и–0 = ф(х–0,у–0). ѕреобразовани€ с использованием условий сопр€жени€ привод€т к системе сеточных уравнений:
ј–г (и–1 - и–0 )ƒ1у + ј–2 (и–2 - и–0 )ƒ1]к + ј–3 (и–3 - и–0 )ƒ1к = 0, (6)
где
ќ" р ќ р
р0 р
јм Ч Ч
„ ќр јпр + ќр ƒпр
„ Ч1,2,3,
0 „ „ о
дл€ текущего конечного объЄма; /1у обозначена длина стороны у. «аметим, что при таком выборе конечного объЄма граничные услови€ Ќеймана (2) - (3) записываютс€ без погрешности аппроксимации.
Ќа барицентрических конечных объЄмах [9] (рис. 2, б) функци€ потенциала
приближаетс€ на множестве базисных функций {“п (х,у1, где N - количество
узлов сетки, полиномом
'(х у ) Ч ип “ п (^ ” )
(7)
в котором ип - значение сеточной функции и в п-й вершине сетки.
 ажда€ базисна€ функци€, как и в ћ Ё, задаетс€ локально на треугольной €чейке сетки по узловым значени€м, причем еЄ значение не равно нулю только в одной вершине треугольника. –ассмотрим треугольник /–0–т–т+1, т = 1, дл€ каждой вершины запишетс€ сво€ базисна€ функци€:
1
'(х у) Ч - 1
“(т) р0
(х,”) Ч
1
2^Д
1
у
”р
1 т
”рт+
“(т) I
2SД
лр о X
”–о

”рт+

(т) р т+1
(x, ”) Ч
2Sm
1 хр р0 ”р0 1 Ч 1 хр р0 ”р0
1 хр т ”рт 1 хр т ”рт
1 х ” 1 хрт+1 ”рт+1
(8)
где Sm - площадь треугольника ƒ–0–т–т+1, т = 1, 2,...ћ–, ћ– - число треугольников вокруг точки –0. ƒл€ одного треугольника базисные функции в сумме дают 1. “огда распределение потенциала на каждом треугольнике интерполируетс€ линейной функцией, принимающей значени€ потенциала в его вершинах. ƒл€ ј–0–т–т+1 в выражении (7) останутс€ только 3 ненулевые базисные функции,
ик (х,”) Ч ир “рт) (х,”) + и
р “р) (X,”) + ир
т
“рт) (х,у).  роме того, градиент по-
тенциала принимает посто€нное значение на /р0ртрт+1. ≈сли подставить градиент в интеграл по контуру (5) и проинтегрировать по участкам границы многоугольного объЄма, которые проход€т в др0ртрт+1, то будет получен вклад в разностную формулу от этого треугольного элемента с посто€нной проводимостью ст. —хема дл€ конечного объЄма вокруг узла р0 (рис. 2, б) запишетс€ в виде
ћр ќ г
^ ((”рт - (т+1)2 + (хр
' 4S
тЧ1 т
- х-,
)2)+
)(хрт+1 - ’рт ))
+ (9)
+ирт ((”рт+1 - ”р0)(”рт - ”рт+1)+(хр0 - х
+ирт+1 ((”р0 - ”рт )(”рт - ”рД +1 ) + (^т ~ хр0 )(^т+1 ~ ^т )» Ч 0,
где суммирование выполн€етс€ по всем треугольным элементам сетки с общей вершиной р0, причем, когда значение индекса т в (9) больше ћр, то т = 1.
пЧ1
Ќа объЄмах ƒирихле - ¬ороного (рис. 2, в) используетс€ другой способ аппроксимации производной в (5). ‘ункци€ потенциала задаетс€ линейной функци-
√радиент линейного потенциала принимает посто€нное значение на ј–0–т–т+1. “огда подстановка градиента и значени€ вектора внешней нормали к текущему отрезку границы конечного объЄма ƒирихле в (5) даЄт вклад в коэффициенты дл€ вершины –0 от ј–0–т–т+1. —хема дл€ конечного объЄма (рис. 2, в) вокруг узла –0 имеет вид
где через —т обозначена точка пересечени€ срединных перпендикул€ров в ј–0–т–т+1; ят, ят+1 - середины сторон –0–т и –0–т+1 соответственно. ѕод |-| подразумеваетс€ длина отрезка.
»нтегральные соотношени€ (5) дл€ граничных узлов на барицентрических объЄмах и объЄмах ƒирихле - ¬ороного аппроксимируютс€ с использованием значений производных от функции ык(х,у), вз€тых из краевых условий (2), (3).  аждое соотношение вносит вклад только в соответствующую компоненту правой части —Ћј”. ¬ схеме дл€ треугольных конечных объЄмов краевые услови€ (2), (3) подставл€ютс€ в криволинейный интеграл дл€ граничных €чеек.
ѕоскольку в сеточных построител€х нельз€ гарантировать Ђправильностьї треугольных €чеек, ћ ќ на €чейках ƒирихле - ¬ороного может давать неоднозначные результаты. ≈сли точка пересечени€ срединных перпендикул€ров одного из треугольников не лежит внутри, то в коэффициентах дл€ вершин накапливаетс€ погрешность.
–ешение задачи Ќеймана дл€ эллиптического уравнени€ (1) - (3) определ€етс€ с точностью до некоторой аддитивной посто€нной [10], при практической реализации у решаемых —Ћј” отсутствует диагональное преобладание. ќднако в измерительной системе Ё»“ всегда есть электрод заземлени€, который также считаетс€ точкой отсчЄта распределени€ электрического потенциала в объекте. Ётот факт можно использовать дл€ вычислительных целей и задавать на одном электроде условие ƒирихле.
«апись соотношени€ (6), (9) или (10) на всей расчЄтной сетке приводит к разреженной симметричной системе линейных алгебраических уравнений (—Ћј”) с диагональным преобладанием относительно неизвестных сеточных значений потенциала при использовании услови€ ƒирихле на одном из электродов. „исла обусловленности систем завис€т от размера сетки, пор€дка нумерации узлов, свойств треугольных элементов сетки, распределени€ проводимости и т.п. Ќа тес-
ей ык(х, у) = ат(х-х–0) + №т(уЧу–0)+и–0 по т-му треугольному элементу сетки:
ат = ({и–тЧи–о)(у–т+\Чу–о) Ч (и–т+1Чи–о)(у–тЧу–о))/(28т),
№т = ((и–т+1 и–0)(х–т ’–0) Ч (и–т и–0)(х–т+1 xP0))/(2Sm), т = ((х–т-х–0)(у–т+1-у–0) Ч (х–т+1Чх–0)(у–тЧу–о))/2-
\/√
(10)
+
3.4. »сследование свойств численных схем
товой сетке (рис. 1, в) число обусловленности равно 133 дл€ конечно-объЄмных схем (9) и (10) и конечно-элементной системы, 1596 - дл€ конечно-объЄмной схемы (6) на однородной задаче (о(х,у) = 1 —м/м в ÷); 1428 и 7647 - на неоднородной задаче (о(х,у) = 5 —м/м в круге радиуса г0 = 5 см и о(х,у) = 1 —м/м вне круга) соответственно. Ќа тестовой сетке (рис. 1, б) получено число обусловленности 2,4104 и 6,1104 на однородной задаче, 3,9104 и 2,6105 - на неоднородной задаче соответственно.
ѕостроенные конечно-объЄмные схемы локально консервативны, так как на любой грани, раздел€ющей соседние конечные объЄмы, выход€щий поток равен вход€щему потоку. Ћокальна€ консервативность построенных разностных схем позвол€ет говорить и об их интегральной консервативности: выполнении интегрального закона сохранени€ тока дл€ всей области ^. »менно выполнение этого услови€ в разностном виде обеспечивает разрешимость полученных систем сеточных уравнений [11]. ”стойчивость разностных схем (6), (9) и (10) доказываетс€ с использованием мажоранты √ершгорина и принципа максимума [12].
3.5. ¬ыбор численного метода решени€ —Ћј”
ƒл€ решени€ полученной —Ћј” были опробованы различные итерационные методы: якоби, √аусса - «ейдел€, верхней релаксации (80я), методы подпространств  рылова Ѕёв и Ѕ1—в8“ЋЅ. ¬ качестве критери€ завершени€ итерационного процесса использовались услови€ малости нев€зки и ошибки (<е = 10-10). ¬ычислительную устойчивость и высокую скорость сходимости на ћ ќ- и ћ Ё-системах сеточных уравнений продемонстрировал 80я со значением параметра релаксации ю = 1,9. Ётот метод выбран дл€ дальнейших вычислений. «аметим, что дл€ ћ ќ- и ћ Ё-схем нет необходимости хранить матрицы коэффициентов целиком. ћетоды решени€ —Ћј” были модифицированы дл€ работы со сжатыми данными, в которых записаны только ненулевые компоненты матриц.
4. “естирование на задаче с аналитическим решением
ѕроведена сери€ экспериментов на задаче с аналитическим решением, в которых оценивалась точность полученных разностных методов. –ассматривалась двумерна€ область ќ радиуса я с круговым включением радиуса г0 в центре (рис. 1, а). јналитическое решение дл€ модели Ё»“ (1) - (3) с неоднородностью в центре круга записываетс€ в пол€рных координатах через р€д ‘урье [13]:
и (г,9) =
ад 1
я 2
п=1 п
(
ад 1
я 2 -
п=1 п
јп I Ч
я
I
я
(ап соЇ(п9) + №п Їш(п9)), 0 < г < г0
-Ѕј -я
(11)
( соЇ(п9) + №п Ї≥п(п9)), г0 < г < я;
—п = -
1
јп =-(1 + *)—п, Ѕп = јп -1;
\2п ' п 9
1+*-(1 - *) € '
1 2п 1 2п
ап = Ч ≤ ](9)соЇ(n9)d9, №п =Ч ≤ ](9)Ї≥п(n9)d9,
ѕ 0 ѕ 0
где ƒ0) - распределение плотности тока по поверхности круга, * = о1/о2, о1 - элек-
трическа€ проводимость включени€, с2 - фонова€ электрическа€ проводимость круга.
ќтносительна€ погрешность вычислений определена через аналитическое решение по формуле
и Ч и
an max an mm
h h где и num - электрическим потенциал, найденный численным методом, и an - аналитическое решение (11), uan_max, uan_min - максимальное и минимальное значение аналитического потенциала в области Q соответственно.
“ест 1. ѕерва€ группа тестов выполнена на сетке (рис. 1, в), представленной в статье G. Dong и др. [8], где авторы провели сравнение ћ Ё- с ћ ќ-подходом, разработанным на барицентрических конечных объЄмах. јвторы не упоминают использование базисных функций (8), но говор€т о линейной интерпол€ции потенциалов по треугольным элементам сетки. ¬ывод формул конечно-объЄмной схемы и детали аналитического решени€ дл€ тестов в [8], в основном, упущены. ќтдельным пунктом выделено различие обработки граничных условий в ћ ќ- и ћ Ё-подходах, когда сила тока задаетс€ на точечных электродах. „исленные эксперименты выполнены на однородном круге с посто€нной проводимостью. ќтносительна€ погрешность находитс€ в пределах 0,3 % дл€ ћ ќ и 0,5 % дл€ ћ Ё в большинстве граничных и внутренних узлов. “очки 1 и 13, где задавались граничные услови€, были исключены из сравнени€ со ссылкой на сингул€рность аналитического решени€ в этих точках и бесконечное значение потенциала в них. ћаксимальна€ относительна€ ошибка вблизи этих точек возрастает до 2,56 % дл€ ћ ќ и 7,4 % дл€ ћ Ё [8].
¬ насто€щей работе приведены результаты сравнени€ ћ ќ с трем€ вариантами выбора конечного объЄма (рис. 2) и ћ Ё на однородном объекте. “есты выполнены дл€ круга радиуса R = 13 см. ƒл€ однородного случа€ электрическа€ проводимость задавалась с = 1 —м/м. Ёлектрический ток инжектируетс€ на паре узлов 1-13 (рис. 1, в), плотность тока J = 1 ј/м2, J13 = -1 ј/м2. ¬ычисленна€ функци€ электрического потенциала сравниваетс€ с аналитическим решением (11) в граничных узлах 1-24 и промежуточных узлах вдоль оси ќх между точками 1 и 13.
 онечно-объЄмные схемы (9) и (10) дают схожее приближенное решение, которое проиллюстрировано графиками (рис. 3). Ќа графиках приводитс€ аналитическое и численное решение в граничных узлах (рис. 3, а), а также вдоль сечени€ области осью ќх (рис. 3, б). ќтносительна€ погрешность методов находитс€ в пределах 0,25 % в граничных точках и 1 % - во внутренних точках вдали от электродов. јппроксимаци€ граничных условий приводит к погрешности до 2,75 % в узлах, относ€щихс€ к электродам. ѕодобные результаты объ€сн€ютс€ тем, что распределение потенциала на каждом элементе сетки использует линейную интерпол€цию потенциала в узлах. ћатрицы коэффициентов схожи на одинаковой сетке, —Ћј” отличаютс€ правой частью, поскольку граничные услови€ обрабатываютс€ по-разному.
“акже были проведены расчЄты ћ Ё дл€ однородного круга. ќтносительна€ погрешность метода сохран€етс€ на том же уровне. ¬ариационна€ формулировка ћ Ё дл€ задачи Ќеймана (1) - (3) получена по [14, 15] на конечномерном подпространстве Uh е H1 (Q), где приближенное решение uh может быть представ-
лено в виде линейной комбинации кусочно-линейных базисных функций (7). –еша€ составленные сеточные уравнение дл€ метода √алеркина относительно иь, находим приближение электрического потенциала в узлах расчЄтной сетки.
ѕриближЄнное решение по ћ ќ на треугольных конечных объЄмах получено в центрах масс треугольных €чеек. ƒл€ сравнени€ с аналитическим решением на границе области численное решение интерполировалось из центров треугольников в граничные точки по направлению перпендикул€ра с использованием разложени€ в р€д “ейлора. Ќа графиках (рис. 3) метод по схеме (6) демонстрирует относительную погрешность до 0,75 % в граничных точках, где нет контакта с электродами, и 5,06 %о - на электродах. Ёлектрический потенциал вдоль сечени€ 1-13 интерполируетс€ по значени€м, которые вычислены в паре смежных объЄмов, имеющих общую сторону на оси ќх. ќтносительна€ погрешность вдоль сечени€ составл€ет 0,17 %о во внутренних узлах, а максимальное значение ошибки 1,49 % находитс€ вблизи электродов 1 и 13.
о4 »
ю
14^
о
и
л ч 2-
о
ќ
и
н
ќ
”гол, рад
------јналитическое решение
Х Х Х ћ ќ, барицентрические €чейки ќќќ ћ ќ, треугольные €чейки
6
<й и ю
¬ 4-1-
х, см
2-
о
о
и
н
о
3 4 5
”гол, рад

~\
-15 -10 -5 0 5 10 15
х, см
–ис. 3. —равнение приближенного ћ ќ-решени€ с аналитическим на круге дл€ случа€ однородного объекта с точечными электродами: а, б - распределение приближенного и аналитического решени€ по внешней границе круга и по сечению у = 0; в, г - относительна€ погрешность вычислений
6
в
г
0
6
7
“ест 2. ¬ другой серии экспериментов конечно-объЄмные разностные схемы (6), (9) и (10) протестированы на модели с электродами размером 2 см (рис. 1, а) дл€ случа€ однородного и неоднородного объекта. “есты выполнены дл€ круга радиуса я = 13 см с внутренним включением размером г0 = 5 см.  огда рассматривалс€ однородный случай, электрическа€ проводимость задавалась как в тесте
1. ¬ экспериментах с неоднородным объектом проводимость включени€ и фона задавалась а1 = 5 —м/м и а2 = 1 —м/м соответственно. –асчЄты выполнены на сетке 4186 элемент (рис. 1, б). ѕлотность тока задана на паре диаметрально размещЄнных электродов в\ и е9, = 1 ј/м2, «13 = -1 ј/м2.
Ќа графиках (рис. 4) аналитическое решение задачи сравниваетс€ с численным ћ ќ-решением дл€ схемы на барицентрических объЄмах. ѕриближенное решение по схемам (9) и (10) незначительно отличаетс€ друг от друга и от ћ Ё решени€ в граничных узлах, расположенных вдали от активной пары электродов, что может быть проиллюстрировано графиком (рис. 4, а). ќтносительна€ погрешность в граничных узлах находитс€ на уровне 0,01 % с отклонением до 0,34 % в узлах, относ€щихс€ к электродам. ¬доль сечени€ объекта осью ќх аналитическое решение сравниваетс€ с вычисленной функцией потенциала, котора€ интерполируетс€ в 100 промежуточных точках (рис. 4, б). ќтносительна€ ошибка достигает
0,8 % внутри объекта и увеличиваетс€ до 4,5 % на границе, что можно объ€снить неточностью интерпол€ции.
–езультаты дл€ конечно-объЄмной разностной схемы (6) показали, что в граничных точках относительна€ погрешность возрастает до 0,9 % в местах контакта активных электродов и держитс€ на уровне 0,2-0,4 % на участках, граничащих с воздухом. Ќа поперечном сечении круга относительна€ ошибка - до 1,5 % во внутренних точках и 5 % вблизи электродов.
”гол, рад х, см
–ис. 4. —равнение приближенного ћ ќ-решени€ с аналитическим решением на круге дл€ случа€ однородного и неоднородного объекта при пропускании тока через электроды е1 - е9: а - распределение потенциала по внешней границе круга; б - распределение потенциала вдоль сечени€ у = 0
ѕри проведении расчЄтов дл€ пр€мой задачи Ё»“ с неоднородным кругом видно, что профиль функции потенциала измен€етс€, когда электрический ток
проходит через неоднородность. ћатематическа€ модель (1) - (3) пригодна дл€ изучени€ изотропных объектов и ограничиваетс€ случаем идеального контакта электродов.
“ест 3. ¬ыполнены расчЄты на группе измельчающихс€ сеток одинакового качества на задаче с неоднородной проводимостью. Ѕолее детальные сетки получены разделением треугольных элементов на четыре части. „тобы сохранить качество сетки, в исходном треугольнике соедин€лись середины сторон, и образовывалось 4 новых треугольных элемента. –езультаты дл€ конечно-объЄмных схем и метода конечных элементов приведены в таблице и говор€т о численной сходимости методов на последовательности сеток с уменьшающимс€ шагом.
ќтносительна€ ошибка конечно-объемных разностных схем и метода конечных элементов
Ќомер схемы* –азмер сетки, €чеек ќтносительна€ ошибка, %
ѕо всей области ѕо граничным узлам
ћаксимум —реднее ћаксимум —реднее
1 790 3,202 0,256 3,202 0,300
3160 2,382 0,211 2,382 0,292
12640 1,795 0,159 1,795 0,230
2 790 1,156 0,060 1,156 0,101
3160 0,347 0,059 0,347 0,067
12640 0,292 0,064 0,292 0,064
3 790 1,157 0,059 1,157 0,101
3160 0,345 0,058 0,345 0,066
12640 0,290 0,063 0,290 0,063
4 790 1,156 0,060 1,156 0,101
3160 0,347 0,059 0,347 0,067
12640 0,292 0,064 0,292 0,064
* ќбозначени€: конечно-объЄмна€ схема дл€ треугольных конечных объЄмов - 1, дл€ барицентрических конечных объЄмов - 2, дл€ конечных объЄмов ƒирихле - ¬ороного - 3, метод конечных элементов - 4.
6. «аключение
ƒл€ численного решени€ пр€мой задачи электроимпедансной томографии на неструктурированных сетках построены конечно-объЄмные аппроксимации. ¬ качестве конечных объЄмов выбирались барицентрические €чейки, €чейки ƒирихле
- ¬ороного, треугольные конечные объЄмы. »сследование аппроксимационных свойств полученных разностных схем затруднено вследствие нерегул€рной структуры вычислительной сетки. ”стойчивость разностных схем доказывалась с использованием принципа максимума и мажоранты √ершгорина. —ходимость провер€лась численно на последовательности сеток с кратно уменьшающимис€ размерами €чеек. Ќа основе сравнительного анализа получаемых численных решений по построенным разностным схемам с точным решением модельной задачи дл€ однородного и неоднородного дисков и с численным решением на основе метода конечных элементов установлено, что высокую точность расчЄтов (сравнимую с ћ Ё) обеспечивают разностные схемы на барицентрических €чейках, €чейках ƒирихле - ¬ороного. Ќесколько худшие результаты, тем не менее с высокой точностью, показывает разностна€ схема на треугольных конечных объЄ-
мах, котора€ позвол€ет без погрешности аппроксимации учесть граничные услови€ Ќеймана.
ƒл€ решени€ пр€мой задачи Ё»“ разработан комплекс программ на €зыке C/C++.
Ћ»“≈–ј“”–ј
1. ѕеккер я.—., Ѕразовский  .—., ”сов ¬.ё. и др. Ёлектроимпедансна€ томографи€. “омск: »зд-во Ќ“Ћ, 2004. 192 с.
2. Lionheart W., Polydorides N., Borsic A. Electrical Impedance Tomography: Methods, History and Applications. Manchester, 2004. 62 p.
3. Brown B.H. Electrical impedance tomography (EIT): a review // J. Med. Eng. Technol. 2003. V. 27. No. 3. P. 97-108. doi: 10.1080/0309190021000059687.
4. Holder D.S. Medical impedance tomography and process impedance tomography: a brief review // Meas. Sci. Technol. 2001. V. 12. P. 991-996.
5. Barber D.C., Brown B.H. Applied Potential Tomography // J. Phys. E: Sci. Instrum. 1984. V. 17. No. 9. P. 723-733. doi:10.1088/0022-3735/17/9/002.
6. Somersalo E., Cheney M., Isaacson D. Existence and uniqueness for electrode models for electric current computed tomography // SIAM J. Appl. Math. 1992. V. 52. No. 4. P. 10231040.
7. Ўерина ≈.—., —тарченко ј.¬. „исленный метод реконструкции распределени€ электрического импеданса внутри биологических объектов по измерени€м тока на границе // ¬естник “омского государственного университета. ћатематика и механика. 2012. є 4(20). —. 36-49.
8. Dong G., Zou J., Bayford R.H., et al. The comparison between FVM and FEM for EIT forward problem // IEEE Trans. Magnetics. 2005. V. 41. No. 5. P. 1468-1471. doi: 10.1109/ TMAG.2005.844558.
9. Li J., Yuan Y. Numerical simulation and analysis of generalized difference method on triangular networks for electrical impedance tomography // App. Math. Modelling. 2009. V. 33. No. 5. P. 2175-2186.
10. “ихоновј.Ќ., —амарский ј.ј. ”равнени€ математической физики. ћ.: Ќаука, 1977. 735 с.
11. ћарчук√.». ћетоды вычислительной математики. ћ.: Ќаука, 1989. 608 с.
12.  рылов ¬.»., Ѕобков ¬.¬., ѕ ћонастырный.». ¬ычислительные методы. ћ.: Ќаука, 1977. 399 с.
13. Isaacson D. Distinguishability of conductivities by electric current computed tomography // IEEE Trans. Medical Imaging. 1986. V. 5. No. 2. P. 91-95.
14. Vauhkonen M. Electrical impedance tomography and prior information. PhD. Kuopio, 1997. 110 p.
15. –о€к ћ.Ё., —оловейчик ё.√., Ўурина Ё.ѕ. —еточные методы решени€ краевых задач математической физики. Ќовосибирск: »зд-во Ќ√“”, 1998. 120 с.
—тать€ поступила 14.04.2014 г.
Sherina E.S., Starchenko A.V. FINITE VOLUME SCHEMES FOR THE ELECTRICAL IMPEDANCE TOMOGRAPHY PROBLEM. In electrical impedance tomography, the currents are applied on electrodes placed on the surface of an object. The electrical conductivity is reconstructed in the interior of the object using voltage measurements on its surface. Knowing the conductivity distribution provides information about the internal object's structure which can be useful for medical and industry applications. Instability of the EIT problem causes difficulties challenging a successful reconstruction. Since a static EIT imaging is sensitive to the measurement noise and approximation errors, this paper addresses the problem of reducing the latter. The finite volume method is presented for solving the EIT forward problem, which is a significant part of the inverse problem. Finite volume schemes were obtained for unstructured grids. The schemes were derived for three kinds of finite volumes, which can be considered relative to triangulation
of the domain. The approaches are based on vertex-centered and cell-centered discretization, where the numerical solution is associated with mesh vertices or mesh elements, respectively. In the first case, a finite volume approximation was introduced on barycentric volumes and Dirichlet
- Voronoi volumes on the assumption of a linear distribution of electric potential. Triangular finite volumes were utilized for approximation based on cell-centered discretization. Both cases suggested a piecewise constant conductivity distribution over grid cells. Numerical comparison for the obtained finite volume schemes was carried out on a test problem that can be solved analytically. The results were compared to a solution by the finite element method.
Keywords: electrical impedance tomography, finite volume method, finite element method, difference schemes, unstructured meshes
Sherina Ekaterina Sergeevna (M. Sc., Tomsk State University, Tomsk, Russian Federation) E-mail: sherina@math.tsu.ru
Starchenko Aleksandr Vasilyevish (Doctor of Physics and Mathematics, Prof., Tomsk State University, Tomsk, Russian Federation)
E-mail: starch@math.tsu.ru
REFERENCES
1. Pekker Ya.S., Brazovskiy K.S.,. Usov V.Yu, Plotnikov M.P., Umanskiy O.S. Elektroim-pedansnaya tomografiya. Tomsk, NTL Publ., 2004. 192 p. (in Russian)
2. Lionheart W., Polydorides N., Borsic A. Electrical Impedance Tomography: Methods, History and Applications. Manchester, 2004. 62 p.
3. Brown B.H. Electrical impedance tomography (EIT): a review (2003) J. Med. Eng. Technol. V. 27. No 3, pp. 97-108. doi: 10.1080/0309190021000059687.
4. Holder D.S. Medical impedance tomography and process impedance tomography: a brief review (2001) Meas. Sci. Technol. V. 12, pp. 991-996.
5. Barber D.C., Brown B.H. Applied Potential Tomography (1984) J. Phys. E: Sci. Instrum. V. 17. No. 9, pp. 723-733. doi:10.1088/0022-3735/17/9/002.
6. Somersalo E., Cheney M., Isaacson D. Existence and uniqueness for electrode models for electric current computed tomography (1992) SIAM J. Appl. Math. V. 52. No. 4, pp. 10231040.
7. Sherina E.S., Starchenko A.V. Chislennyy metod rekonstruktsii raspredeleniya elek-tricheskogo impedansa vnutri biologicheskikh ob"ektov po izmereniyam toka na granitse (2012) Vestnik Tomskogo gosudarstvennogo universiteta. Matematika i mekhanika. No. 4(20), pp. 36-49. (in Russian)
8. Dong G., Zou J., Bayford R.H., Xinshan M., Shangkai G., Weili Y., Manling G. The comparison between FVM and FEM for EIT forward problem (2005) IEEE Trans. Magnetics. V. 41. No. 5, pp. 1468-1471. doi: 10.1109/TMAG.2005.844558.
9. Li J., Yuan Y. Numerical simulation and analysis of generalized difference method on triangular networks for electrical impedance tomography (2009) App. Math. Modelling. V. 33. No. 5, pp. 2175-2186.
10. Tikhonov A.N., Samarskiy A.A. Uravneniya matematicheskoy fiziki. Moscow, Nauka Publ., 1977. 735 p.
11. Marchuk G.I. Metody vychislitel'noy matematiki. Moscow, Nauka Publ., 1989. 608 p. (in Russian)
12. Krylov V.I., Bobkov V.V., Monastyrnyy P.I. Vychislitel'nye metody. Moscow, Nauka Publ., 1977. 399 p.
13. Isaacson D. Distinguishability of Conductivities by Electric Current Computed Tomography (1986) IEEE Trans. Medical Imaging. V. 5. No. 2, pp. 91-95.
14. Vauhkonen M. Electrical impedance tomography and prior information. PhD. Kuopio, 1997. 110 p.
15. Royak M.E., Soloveychik Yu.G., Shurina E.P. Setochnye metody resheniya kraevykh zadach matematicheskoy fiziki. Novosibirsk: NGTU Publ., 1998. 120 p. (in Russian)

пїњ