Finite volume schemes for the electrical impedance tomography problem.pdf 1. Введение Электроимпедансная томография (ЭИТ) - это быстро развивающийся неинва-зивный метод визуализации со спектром применения, включающим медицинскую томографию, неразрушающий контроль материалов и контроль промышленных процессов [1-5]. Метод ЭИТ восстанавливает неизвестное распределение или оценивает несколько параметров электрической проводимости (или полного импеданса) внутри объекта, если измерено электрическое напряжение или сила тока на его границе. По реконструкции даётся оценка внутренней структуры объекта. ЭИТ позволяет наблюдать динамические процессы и работать со статическими изображениями. Отмечают несколько преимуществ ЭИТ перед другими видами томографии: компактный размер оборудования, простота устройства, невысокая стоимость и быстрый сбор данных. Кроме того, в медицинских приложениях ценится безопасность для биологических тканей. Направление ЭИТ представляется перспективным в разных приложениях, но его широкое внедрение ограничено невысоким качеством получаемых изображений. Наравне с необходимостью усовершенствования аппаратных средств томографии следует уделять внимание разработке математической теории ЭИТ и новых численных методов решения такой задачи. Восстановление проводимости с помощью ЭИТ является нелинейной и некорректной обратной коэффициентной задачей. Сложность задачи вызвана нарушением условия устойчивости в классическом определении корректности по Адама-ру [2]. Как следствие, реконструкция статических изображений ЭИТ очень чувствительна к ошибкам измерений прибора и погрешности аппроксимации дифференциальной задачи разностной, а также качеству физического моделирования рассматриваемого объекта. Данное исследование направлено на изучение особенностей численного моделирования прямой задачи ЭИТ, являющейся важной компонентой обратной задачи. Отдельное внимание уделено вопросу уменьшения погрешности аппроксимации, которая появляется в ходе дискретизации задачи. Все результаты получены для двумерной статической задачи и при необходимости могут быть перенесены на трёхмерный случай или использоваться в динамической реконструкции. 2. Постановка задачи ЭИТ В ЭИТ на поверхность объекта прикрепляется набор электродов (рис. 1, а), затем выполняется серия измерений, в которых через объект пропускается электрический ток небольшой силы и высокой частоты. Обычно, сила тока выбирается в интервале от 1 до 5 мА с частотой от 10 до 150 кГц. Схема проведения эксперимента - схема подведения тока и сбора данных - определяется исследователем. Между электродами граница объекта контактирует с воздухом. Результирующее напряжение измеряется на электродах. Вся серия данных «электрический ток -напряжение» используется для реконструкции неизвестного распределения электрической проводимости внутри объекта. 2.1. Физическая постановка задачи В данной статье использована модель ЭИТ с идеально проводящими электродами, которая называется в научной литературе gap model [6]. Предполагается, что учитываемые моделью физические процессы описываются электрическим и магнитным полями, которые удовлетворяют стационарным уравнениям Максвелла. В данной работе рассматривается случай применения ЭИТ к объектам биологической природы. Ткани проводят электричество, потому что они содержат ионы, которые переносят электрический заряд. Проводимость биологических тканей обладает определенной спецификой и зависит от частоты тока. Когда электромагнитные волны проходят через биологический объект, интенсивность магнитного поля становится пренебрежимо малой [1, 2], поэтому магнитная компонента в уравнениях Максвелла не рассматривается. В этом случае система уравнений Максвелла может быть с использованием закона Ома сведена к уравнению эллиптического типа в частных производных. В силу отсутствия внутри объекта источников электрического поля, а также с учётом 1-го правила Кирхгофа для электрической цепи, сумма втекающих и вытекающих токов равняется нулю - выполняется закон сохранения заряда. Также предполагается, что на электроде, где инжектируется ток, известно значение потенциала электрического тока. (1) (2) (3) 2.2. Математическая постановка задачи С учётом принятых условий математическая модель ЭИТ для двумерного случая (рис. 1, а) может быть представлена краевой задачей [2, 7]: V-(cVcp) = 0, (x,y)еП ; dp ст- = Jk, (xy)е ek, k = L ; дп ст-= 0, (x,y)е dQ/ u ek ; дп k=1 где o(x,y) - электрическая проводимость, ф(х,у) - потенциал электрического поля, п - единичный вектор внешней нормали, Jk - плотность электрического тока на k-м электроде, ek - k-й электрод, L - число электродов. Функция электрического потенциала фе H1 (Q), плотность тока j =ст • Уф е L2 (Q). По принятой в литературе терминологии задача (1) - (3) с функцией электрического потенциала 9(x,y) в качестве неизвестного будет в дальнейшем называться прямой задачей ЭИТ. Обратной задачей ЭИТ именуется задача (1) - (3) по поиску неизвестной функции проводимости o(x,y) в области дополненная данными измерений электрического напряжения на электродах. 3. Построение численных схем для задачи ЭИТ Чувствительность реконструкции ЭИТ к ошибкам измерений и погрешности аппроксимации дифференциальных операторов ограничивает круг подходящих численных методов, которыми можно с хорошим качеством решать обратную задачу. Как правило, дифференциальные задачи приводят к дискретному аналогу каким-либо сеточным методом. Среди них можно выделить следующие группы: конечно-разностные (МКР), конечно-элементные (МКЭ), гранично-элементные (МГЭ), конечно-объёмные (МКО) и другие. Задачу ЭИТ изначально ассоциировали с электродинамикой и моделированием цепей, поэтому был опробован, а впоследствии и наиболее распространился подход по МКЭ. МКЭ обладает гибкостью, когда необходимо быстро улучшить качество аппроксимации без сложных изменений алгоритма. Но в [8] отмечается, что МКЭ не удовлетворяет условию непрерывности электрического тока, тем самым не соблюдается закон сохранения на элементе сетки и области в целом. С МКЭ ищется решение не конкретной дифференциальной задачи, а обобщённое решение вариационной задачи, выведенной из исходной и записанной в интегральной форме. В данной работе приводятся некоторые результаты построения конечно-объёмных схем для прямой задачи ЭИТ, которые сравниваются с решением МКЭ по нескольким показателям. 3.1. Дискретизация расчётной области На начальном этапе решения задачи выполняется дискретизация расчётной области rah, и определяется набор сеточных функций электрического потенциала uh и проводимости ch Создание расчётных сеток - отдельная область исследования. В настоящее время доступны свободно распространяемые и коммерческие программные продукты для генерации сеток (GMSH, GAMBIT, QMG, Netgen, Tgrid, FEMLAB и др.). В данном исследовании численные эксперименты проводились на треугольных неструктурированных сетках (M ячеек, N узлов), построенных в программах GAMBIT и GMSH. Пример расчётной сетки приведён на рис. 1, б. Заметим, что сетка детализована вдоль границы области, потому что градиент потенциала быстро изменяется вблизи контактов электродов. Для оценки качества построенных сеток с треугольными ячейками обычно вводится какой-нибудь числовой показатель. Чтобы осуществлять контроль качества построенных сеток была выбрана группа параметров: коэффициент деформации элемента сетки, угловой коэффициент, соотношение сторон треугольного элемента, соотношение радиуса вписанной к радиусу описанной окружности для элемента сетки. Как правило, точность расчётов снижается, если элементы треугольной сетки сильно отличаются от правильных треугольников. Следует заметить, что решающим фактором при выборе сетки будет не только набор показателей качества и размер сетки, но и её способность обеспечивать работоспособность конкретного применяемого численного метода в целом. Должна быть обеспечена вычислительная экономичность расчётов и достаточная точность приближенного решения полученных сеточных уравнений, поэтому, по-видимому, следует говорить о том, насколько хорошо сетка подходит для применяемого численного метода. e5 7 1э\/ \/ V у\ 5 Рис. 1. Двумерная модель области реконструкции Q (а), на границе dQ обозначены электроды ek, к = 1,...,16; расчётная сетка со сгущением к электродам (б), размер сетки -M = 4168 элементов, N = 2181 узел; расчётная сетка для тестовой задачи с точечными электродами 1 и 13 (в), размер сетки - M = 136 элементов, N = 81 узел [8] 10 0 ---ev^X ^"Чез ( e \ Q J ev 3.2. Варианты конечных объёмов Когда дифференциальная задача с помощью сеточного метода преобразуется в разностную, то выбирается способ соответствия, по которому неизвестные значения сеточной функции связываются с расчётной сеткой. Функции электрического потенциала ф(х,х) и проводимости а(х,,у) заменяются дискретными аналогами, т.е. сеточными функциями uh и ah соответственно, uh е Uh,ah eSh, U и Zh - гильбертовы пространства. На триангуляции области обычно рассматриваются два варианта - значение сеточных функций определяется в некоторой точке внутри ячейки сетки (рис. 2, а) или в узле сетки (рис. 2, б, в). Конечный объём в обоих случаях строится вокруг точки, с которой связываются значения сеточной функции. В настоящей работе обсуждается реализация МКО на трёх конечных объёмах: треугольный конечный объём (рис. 2, а), барицентрический конечный объём (рис. 2, б) и конечный объём - ячейка Дирихле - Вороного (рис. 2, в). P2 P Рис. 2. Фрагменты расчётной сетки: а - треугольный конечный объём с центром в точке P0; б - барицентрический конечный объём, построенный вокруг вершины P0; в - конечный объём - ячейка Дирихле - Вороного, - построенный вокруг вершины P0 P Барицентрическая ячейка строится вокруг вершины P0 по точкам пересечения медиан и точкам середин сторон треугольников. В ячейке Дирихле - Вороного вместо барицентров выбираются точки пересечения срединных перпендикуляров. Треугольные конечные объёмы соответствуют ячейкам сетки, в них положение точки P0 задается в центре масс, точке пересечения срединных перпендикуляров или биссектрис. 3.3.Вывод численных схем Согласно принципу МКО для каждого конечного объёма записывается ин-тегро-интерполяционный баланс. Уравнение (1) интегрируется по конечному объёму Qm, | У-(стУф)^Q = 0. (4) Qm По формуле Грина уравнение (4) преобразуется к равенству нулю криволинейного интеграла по границе от плотности тока стдф dl = 0. (5) dQ dn dQm Дальнейшая оценка интеграла в левой части (5) зависит от формы конечного объёма и способа аппроксимации градиента потенциала в подынтегральном выражении. Функция проводимости задается кусочно-постоянной по треугольным элементам сетки, 0 < c < o(x,y) < C, c и C - константы. На треугольных конечных объёмах (рис. 2, а) интеграл в правой части (5) вычисляется по сторонам текущего треугольника Qm по интерполяционной формуле, в данном случае - средних прямоугольников. Опробовано несколько вариантов расположения точки P0 внутри конечного объёма, с которой ассоциировались сеточные функции проводимости ch и электрического потенциала uh. Дискретные значения определялись в точке пересечения медиан, срединных перпендикуляров или биссектрис треугольника. Производные аппроксимировались в серединах сторон конечного объёма, причём оценка выполнялась с учётом условий сопряжения на границе двух смежных ячеек, например для грани ij (рис. 2, а): стдфТ° г дф" dn )R v dn К =(ф)^ чтобы схема стала чувствительной к изменению свойств токопроводящей среды. Потоки через границу конечного объёма могут быть аппроксимированы разными соотношениями. Один из предложенных подходов рассматривает аппроксимацию (дфУ° uR - uP (дфУ* up - uR вида I - I « -1-- и I - I и -1-L, где AnP(1 - длина отрезка перпендиv dn )r1 Anp v dn )r1 AnP1 куляра из точки P0 до стороны ij, AnPi - от точки P1 до стороны ij, а uP0 s 9(xP0,yP0). Преобразования с использованием условий сопряжения приводят к системе сеточных уравнений: AP1 (uP1 - uP0 ) Alij + AP2 (uP2 - uP0 ) Aljk + AP3 (uP3 - uP0 ) Alkr = 0, (6) а р а р P0 Pq q = 1,2,3, где AM = q ар AnP +ар ДпР 0 q q 0 для текущего конечного объёма; Aly обозначена длина стороны ij. Заметим, что при таком выборе конечного объёма граничные условия Неймана (2) - (3) записываются без погрешности аппроксимации. На барицентрических конечных объёмах [9] (рис. 2, б) функция потенциала приближается на множестве базисных функций {Чn (x,y)}}= , где N - количество узлов сетки, полиномом (7) '(x,У) = XипЧn (x,У), n=1 в котором un - значение сеточной функции и в n-й вершине сетки. Каждая базисная функция, как и в МКЭ, задается локально на треугольной ячейке сетки по узловым значениям, причем её значение не равно нулю только в одной вершине треугольника. Рассмотрим треугольник AP0PmPm+1, m = 1, для каждой вершины запишется своя базисная функция: 1 yPo У yPm+ ЛР 0 x У yP '(*, у ) = - 1 Ч (m) Po Ч(m) i 2S„ 1 ypm+ (8) (m) p m+\ (x У) = Ч 1 (х,У) = m 2S„ 2S„ 1 Хр0 УР0 _ 1 = 2 1 Хр0 УР0 1 Pm yPm , Sm 1 Pm yPm 1 x У 1 xPm+1 yPm+1 где Sm - площадь треугольника AP0PmPm+1, m = 1, 2,...MP, MP - число треугольников вокруг точки P0. Для одного треугольника базисные функции в сумме дают 1. Тогда распределение потенциала на каждом треугольнике интерполируется линейной функцией, принимающей значения потенциала в его вершинах. Для AP0PmPm+1 в выражении (7) останутся только 3 ненулевые базисные функции, uh (x,У) = uP ЧPm) (x,y) + и P ЧP) (x,У) + Up m Ч(pm) (x,y). Кроме того, градиент по тенциала принимает постоянное значение на AP0PmPm+1. Если подставить градиент в интеграл по контуру (5) и проинтегрировать по участкам границы многоугольного объёма, которые проходят в AP0PmPm+1, то будет получен вклад в разностную формулу от этого треугольного элемента с постоянной проводимостью cm. Схема для конечного объёма вокруг узла P0 (рис. 2, б) запишется в виде Mp а X 1ST К ((yPm - yPm+1)2 + (xP - x-, ' 4S m=1 m + (9) +upm ((ypm+1 - Ур0)(ypm - ypm+1)+(xp0 - x +UPm+1 ((УР0 - ypm )(ypm ~ Урт +1 ) + (^m ~~ ^0 )(^m+1 - ^m )И = 0, где суммирование выполняется по всем треугольным элементам сетки с общей вершиной Р0, причем, когда значение индекса m в (9) больше Mp, то m = 1. )2)+ )(xPm+1 - xPm )) На объёмах Дирихле - Вороного (рис. 2, в) используется другой способ аппроксимации производной в (5). Функция потенциала задается линейной функцией uh(x, y) = am(x-xP0) + bm(y-yP0)+uP0 по m-му треугольному элементу сетки: am = ((uPm-uP0)(yPm+1-yP0) - (uPm+1-uP0)(УPm-УP0))/(2Sm), bm = ((uPm+1-uP0)(xPm-xP0) - (uPm-uP0)(xPm+1-xP0))/(2Sm), Sm = ((xPm xP0)(yPm+1 -yP0) - (xPm+1-xP0)(yPm-yP0))/2. Градиент линейного потенциала принимает постоянное значение на AP0PmPm+1. Тогда подстановка градиента и значения вектора внешней нормали к текущему отрезку границы конечного объёма Дирихле в (5) даёт вклад в коэффициенты для вершины P0 от AP0PmPm+1. Схема для конечного объёма (рис. 2, в) вокруг узла Р0 имеет вид ARmCm J(xP - xP )2 + (yP - yP )2 . V Pm P0 ' y Pm ■/p0/ (10) MP m-1 - uP P0 - up) \CmRm+11 -0, ( P+1 - xp0)2 +(ypm+1 - yp0)2 где через Cm обозначена точка пересечения срединных перпендикуляров в AP0PmPm+1; Rm, Rm+1 - середины сторон P0Pm и P0Pm+1 соответственно. Под |-| подразумевается длина отрезка. Интегральные соотношения (5) для граничных узлов на барицентрических объёмах и объёмах Дирихле - Вороного аппроксимируются с использованием значений производных от функции uh(x,y), взятых из краевых условий (2), (3). Каждое соотношение вносит вклад только в соответствующую компоненту правой части СЛАУ. В схеме для треугольных конечных объёмов краевые условия (2), (3) подставляются в криволинейный интеграл для граничных ячеек. Поскольку в сеточных построителях нельзя гарантировать «правильность» треугольных ячеек, МКО на ячейках Дирихле - Вороного может давать неоднозначные результаты. Если точка пересечения срединных перпендикуляров одного из треугольников не лежит внутри, то в коэффициентах для вершин накапливается погрешность. Решение задачи Неймана для эллиптического уравнения (1) - (3) определяется с точностью до некоторой аддитивной постоянной [10], при практической реализации у решаемых СЛАУ отсутствует диагональное преобладание. Однако в измерительной системе ЭИТ всегда есть электрод заземления, который также считается точкой отсчёта распределения электрического потенциала в объекте. Этот факт можно использовать для вычислительных целей и задавать на одном электроде условие Дирихле. 3.4. Исследование свойств численных схем Запись соотношения (6), (9) или (10) на всей расчётной сетке приводит к разреженной симметричной системе линейных алгебраических уравнений (СЛАУ) с диагональным преобладанием относительно неизвестных сеточных значений потенциала при использовании условия Дирихле на одном из электродов. Числа обусловленности систем зависят от размера сетки, порядка нумерации узлов, свойств треугольных элементов сетки, распределения проводимости и т.п. На тестовой сетке (рис. 1, в) число обусловленности равно 133 для конечно-объёмных схем (9) и (10) и конечно-элементной системы, 1596 - для конечно-объёмной схемы (6) на однородной задаче (o(x,y) = 1 См/м в Q); 1428 и 7647 - на неоднородной задаче (o(x,y) = 5 См/м в круге радиуса r0 = 5 см и o(x,y) = 1 См/м вне круга) соответственно. На тестовой сетке (рис. 1, б) получено число обусловленности 2,4104 и 6,1104 на однородной задаче, 3,9104 и 2,6105 - на неоднородной задаче соответственно. Построенные конечно-объёмные схемы локально консервативны, так как на любой грани, разделяющей соседние конечные объёмы, выходящий поток равен входящему потоку. Локальная консервативность построенных разностных схем позволяет говорить и об их интегральной консервативности: выполнении интегрального закона сохранения тока для всей области Именно выполнение этого условия в разностном виде обеспечивает разрешимость полученных систем сеточных уравнений [11]. Устойчивость разностных схем (6), (9) и (10) доказывается с использованием мажоранты Гершгорина и принципа максимума [12]. 3.5. Выбор численного метода решения СЛАУ Для решения полученной СЛАУ были опробованы различные итерационные методы: Якоби, Гаусса - Зейделя, верхней релаксации (SOR), методы подпространств Крылова BiCG и BiCGSTAB. В качестве критерия завершения итерационного процесса использовались условия малости невязки и ошибки (

Скачать электронную версию публикации