CC BY
115
ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА
2012 Математика и механика № 4(20)
УДК 519.6
Е.С. Шерина, А.В. Старченко
ЧИСЛЕННЫЙ МЕТОД РЕКОНСТРУКЦИИ РАСПРЕДЕЛЕНИЯ ЭЛЕКТРИЧЕСКОГО ИМПЕДАНСА ВНУТРИ БИОЛОГИЧЕСКИХ ОБЪЕКТОВ ПО ИЗМЕРЕНИЯМ ТОКА НА ГРАНИЦЕ1
Электроимпедансная томография (ЭИТ) применяется для восстановления неизвестного распределения электрической проводимости внутри объектов живой природы, обладающих неоднородной структурой, по измерениям силы и напряжения электрического тока на границе. В данной работе предложен численный метод решения обратных задач ЭИТ, опирающийся на использование дифференциальной эволюции.
Ключевые слова: электроимпедансная томография (ЭИТ), неструктурированная треугольная сетка, разностная схема, коэффициентная обратная задача, генетические алгоритмы оптимизации, метод дифференциальной эволюции.
В настоящее время в медицине и промышленности для компьютерной диагностики применяется способ получения томографического изображения, называемый электроимпедансной томографией (ЭИТ) или томографией приложенных потенциалов [1-5]. Он основан на использовании электрического тока в качестве средства, зондирующего исследуемый объект. Цель ЭИТ состоит в том, чтобы восстанавливать неизвестное распределение электрической проводимости внутри объекта с неоднородной структурой, используя электрические измерения только на границе объекта. Слабый электрический ток высокой частоты подводится к объекту исследования посредством системы токопроводящих электродов, прикрепляемых к его поверхности. На них выполняются измерения электрического напряжения или тока с целью получения исходной совокупности данных для реконструкции. С физической точки зрения, ЭИТ как диагностический метод с некоторой погрешностью воспроизводит распределение электрической проводимости или электрического сопротивления (импеданса) внутри объекта. Поскольку различные вещества характеризуются своими значениями электрической проводимости (сопротивления), знание характера распределения электрической проводимости позволяет установить внутреннюю структуру объекта, что может быть использовано в медицине и промышленности [3]. Однако восстановление распределения скалярной функции проводимости внутри области по измерениям электрического тока на границе является некорректной задачей и требует применения специальных методов реконструкции, включающих решение обратной задачи [3, 4].
В данной работе предложен метод решения задач ЭИТ для исследования внутренней структуры биологических объектов, основанный на алгоритме дифференциальной эволюции [6, 7]. Использованный стохастический метод - дифференци-
1 Работа выполнена по заданию Министерства образования и науки РФ № 8.4859.2011 при финансовой поддержке ФЦП «Научные и научно-педагогические кадры инновационной России на 2009-2013» (госконтракт № 14.B37.21.0667), РФФИ (грант №12-01-31050 мол_а).
альная эволюция - относится к классу эвристических эволюционных алгоритмов, применяемых для решения задач оптимизации. Подобные алгоритмы итеративно моделируют эволюционный процесс, в котором в качестве популяции рассматривается множество допустимых решений задачи. В процессе моделирования учитываются основные механизмы теории биологической эволюции, такие, как мутация, скрещивание и селекция. Применение генетических алгоритмов представляет интерес для решения обратных задач, поскольку они позволяют обойти главную трудность задач этого класса, а именно, характерную некорректность их постановки. В случае задач ЭИТ их некорректный характер вызван нарушением устойчивости процесса получения решения, которое возникает из-за недостаточной точности и количества входных данных.
В предложенном методе алгоритм оперирует популяцией векторов электрической проводимости, среди которых может оказаться возможное решение задачи. Поиск «оптимального» решения осуществляется путём циклического изменения текущей популяции по законам эволюции и оценки значения целевой функции для каждого вектора популяции. По завершению работы алгоритма выбирается представитель популяции, для которого целевая функция принимает наименьшее значение, и данный вектор считается приближенным решением задачи. Метод дифференциальной эволюции был выбран в качестве способа минимизации в силу его простоты, высокой скорости работы и прямолинейной стратегии поиска. Принимая во внимание результаты, полученные в работе [6], можно сделать вывод, что метод дифференциальной эволюции в большинстве случаев превосходит все другие методы оптимизации по требуемому количеству оценок функции, необходимых для определения минимума целевой функции с заданной точностью. Алгоритм легко использовать, поскольку он прост и имеет небольшое число управляющих параметров, которые можно выбрать из заранее определенных числовых интервалов.
Постановка задачи
Разработанный метод решения обратных задач ЭИТ реализован на примере двумерной области D = D и Г (рис. 1), которая представляет собой схематическое изображение некоторого биологического объекта, имеющего неоднородную структуру. Границы неоднородностей внутри области заранее известны. В каждой подобласти коэффициент электрической проводимости с имеет некоторое постоянное значение, поэтому ст(х, у) можно рассматривать как кусочно-постоянную
функцию от координат (х,у) е D . Необходимо определить значения коэффициента с внутри области, если априори известны измерения электрического напряжения и силы тока на границе Г.
Чтобы получить исходные данные для поставленной задачи, через электропроводящий объект, который находится в воздушной среде, пропускается электрический ток небольшой силы 1-5 мА высокой частоты 100-150 кГц. Ток инжектируется через расположенные локально на его поверхности электроды, на которых выполняются измерения разности электрических потенциалов [3, 4]. Между электродами граница Г контактирует с воздухом. Использована модель ЭИТ с идеально проводящими электродами, называемая в научной литературе shunt model [3, 4, 8], и реализована одна из схем организации измерений, в которой электроды включаются попарно и электрическое напряжение определяется на ак-
тивной паре. Обращаясь к физической постановке задачи, предполагаем, что учитываемые моделью процессы описываются электрическим и магнитным полями, которые удовлетворяют стационарным уравнениям Максвелла. В случае применения ЭИТ к объектам биологической природы интенсивность магнитного поля становится пренебрежимо малой из-за особенности в проводимости биологических тканей [3, 4], и поэтому магнитная компонента не будет рассматриваться. При таких условиях вектор напряженности электрического поля равен градиенту электрического потенциала ф, распределение которого внутри объекта исследования будет удовлетворять уравнению эллиптического типа [9]. В силу отсутствия внутри объекта источников электрического поля сумма втекающих /ш и вытекающих /ои4 токов равняется нулю согласно закону сохранения заряда, /1п + /ои = 0. Также предполагаем, что на границе, где инжектируется ток, известно значение потенциала электрического тока фт.
С учётом описанной выше физической постановки математическая модель ЭИТ представляет собой краевую задачу для уравнения эллиптического типа в частных производных [3, 4]:
Рис. 1. Пример модели биологического объекта Б с границей Г, на которую прикреплены два электрода с площадью поверхности £1п = £ои4. В объект инжектируется электрический ток величиной /1п
(1)
(2)
Ф(ХУ) = Фт , (ХУ)еГт;
(3)
о (Х У Ки Х,У) = ^о* = ^ ( х У) 6 Гои
дп
(4)
где о - коэффициент электрической проводимости, ф - потенциал электрического поля, п - единичный вектор внешней нормали, Jn - плотность тока на границе, ^1п,
- площадь контактной поверхности 1п- и оШ-электродов соответственно. Уравнение (1) показывает зависимость скалярных функций электрического потенциала ф(х,у) и проводимости о(х,у) внутри объекта. Краевое условие (2) выражает контакт тела с воздухом, (3) и (4) - контакт тела с электродами.
При некотором априори заданном распределении значений функции проводимости о и значениям электрического потенциала ф1п и электрического тока /1п на границе из краевой задачи (1) - (4) можно найти неизвестное распределение потенциалов ф электрического поля внутри области Б . В общем случае эта задача не может быть решена аналитически [3, 4].
На основе постановки (1) - (4) можно сформулировать задачу по восстановлению неизвестного распределения электрической проводимости о(х,у) в области Б, если заданы значения силы инжектируемого тока /1п,р (р = 1,...,Ь) и измерены значения разности электрического потенциала Лфртешшгеа на границах контактов тела с электродами, т.е.:
ф(х У ) = фт,р , (^ У) 6 Гт,р ; (5)
ф(х,у) = фш,р +Дфтеакигей, (х,У)б Го*,р, р = 1,...,Ь, (6)
где Ь - количество использованных конфигураций электрического тока и активных пар электродов.
Переходя к принятой терминологии [1-5], задача (1) - (4) с функцией электрического потенциала ф(х,у) в качестве неизвестного будет в дальнейшем называться прямой задачей ЭИТ. Задача поиска неизвестной функции проводимости о(х,у) в области Б , которая удовлетворяет уравнению (1) в Б и граничным условиям (2) - (6) на Г, - это обратная задача ЭИТ.
В данной работе для решения обратной задачи ЭИТ предложено использовать метод, осуществляющий перебор допустимых распределений функции проводимости о(х,у), для каждого из которых по решению прямой задачи вычисляется значения электрического потенциала внутри и на границе области. Основываясь на результатах решения совокупности Ь прямых задач с заданным распределением проводимости о(х,у), выполняется поиск «оптимального» оор1;(х,у), удовлетворяющего условию минимума функционала
/(ст) = £ (Дфр (а)-Дфтеакигей)2, (7)
р=1
где Лфр(о) = фо1й,р(о) - ф1П,р(о) - разность электрических потенциалов, вычисленных при р-ой конфигурации электродов, заданном значении силы электрического тока /1п,р и при заданном распределении о(х,у) в Б , Лфртеашге<1 - измеренная разность потенциалов при тех же параметрах тока и конфигурации активной пары электродов. Следует отметить, что в процесс минимизации вовлечены все возможные пары электродов, поскольку при подобном подходе немаловажной деталью является общее количество выполненных измерений. На каждой паре электродов через объект инжектируется М образцов тока, следовательно, максимально возможное число измерений равно Ь = MxQ(Q - 1)/2, где Q - общее количество электродов. Большой набор входных данных обеспечивает разрешимость обратной задачи ЭИТ и допускает восстановление соответствующего числа параметров электрической проводимости о(х,у) [3, 4].
Численный метод решения дифференциальной задачи
Для практической реализации математической модели и разработки численного метода решения прямой задачи ЭИТ была построена неструктурированная треугольная сетка О} (например, расчётная сетка на рис. 2) с помощью сеточного генератора GMSH [10] или GAMBIT [11]. Контур области D и неоднородности внутри неё задавались набором из 40 точек. Полученная сетка строилась как сетка Делоне, для этого использовались соответствующие встроенные возможности рассмотренных построителей сеток.
Рис. 2. Расчётная сетка
Следует сделать несколько замечаний, касающихся расчётных сеток для ЭИТ. Оптимальные сетки для решения задач ЭИТ должны обладать хорошим качеством (с геометрической и вычислительной точек зрения), чтобы представлять потенциал электрического поля с достаточной точностью для последующего определения электрической проводимости. Это означает, что при моделировании нужно адекватно представлять форму поверхности визуализируемой области и геометрию электродов. Кроме того, сетка должна обеспечивать большую плотность узлов в областях с большими градиентами электрического поля, чтобы уменьшить влияние численных ошибок, в частности, вблизи поверхностей электродов [4]. Для этого необходимо сгущать сетку (измельчать её ячейки) вблизи электродов. Однако часто подобное построение сетки вблизи электродов не увеличивает точность решения и иногда может привести к большим вычислительным затратам. При решении обратной задачи интерес представляют изменения проводимости внутри области, которые в любом случае будут найдены в масштабе не меньшем чем размеры электродов, поэтому параметризация проводимости будет неизбежно грубее, чем аппроксимация производных для потенциала [4]. Кроме того, сетка должна учитывать подобласти, в которых проводимость известна или является константой. Ещё один аспект заключается в том, что напряженность электрического поля изменяется в зависимости от используемых образцов силы тока, потому на практике следует выбирать сетки, подходящие для всех конфигураций тока, следовательно, в некоторых случаях сетки могут отличаться излишне хорошим качеством вдали от активных электродов.
Использованные программы-генераторы сеток GAMBIT и GMSH позволяют осуществлять проверку качества построения по различным параметрам («скошен-
ность» ячеек сетки и соотношение сторон треугольных ячеек). Улучшение сетки можно проводить вручную или используя операцию сглаживания. Сглаживание сетки подразумевает корректировку расположения узлов построенной сетки с целью улучшения единообразия расстояния между узлами по всей поверхности.
Обобщая вышесказанное, можно отметить, что неструктурированные расчётные сетки являются оптимальным классом сеток для двумерной задачи ЭИТ, так как они гибко отображают геометрию объекта и позволяют сгущать узлы сетки в областях больших градиентов. При генерации расчётных сеток для ЭИТ предпочтение отдается сеткам с «правильными» ячейками, т.е. ячейками, элементы которых близки к правильным треугольникам.
Разностный аналог дифференциальной задачи (1) - (4) построен на основе конечно-разностного метода конечных объёмов (МКО) [12-14], который обеспечивает выполнение интегрального закона сохранения для каждого конечного объёма и для всей области в целом. В качестве конечных объёмов выбраны треугольные элементы сетки, что позволяют разделить границы областей с разными физическими свойствами, а именно, характеризующиеся различными значениями электрической проводимости с(х,у). Приближенные значения функции потенциала электрического поля ф(х,у) и функции электрической проводимости с(х,у) ищутся в центрах масс треугольных ячеек (рис. 3).
Распределение электрической проводимости, появляющееся в различных приложениях метода ЭИТ, обычно является кусочно-постоянным. Такой случай подходит, например, для медицинской ЭИТ, так как различные ткани в организме описываются разными значениями функции проводимости с разрывами на границах внутренних органов. В данной работе рассмотрен случай кусочно-постоянной функции электрической проводимости с внутри области Б . При дискретизации краевой задачи использованы условия сопряжения (граничные условия 4-го рода)
для обеспечения идеального контакта двух сред, имеющих разные проводимости. При использовании формулы Грина при интегрировании (1) по конечному объему АЛ^Л^Лз (рис. 3), формулы средних прямоугольников для приближенного вычисления криволинейного интеграла и условий сопряжения на границах ячеек с различной электрической проводимостью, получена следующая разностная схема [14]:
* *
СТМ1 ФМ1 -фр (А/ ) , СТМ2 Фм2-фр (А/
* * о/С , сМ'- 42' + * * т/о , сМ'-
СТР +СТМ1 2 (З12 + З12 ) СТР +СТМ2 2 (З23 + З23 )
*
стМз фМз -фР , Ч2 п/п\ /оч
+——^ / 3 о* ч (А/з:) = 0,(р)еОА; (8)
аР +аМ3 2 (1 + З31 )
Фм,. -Фр = 0, (М , )6(И )а^ (9)
I Ап
Фм, = Фт +-------11 ч7 , (Мг ) 6 (И )т , (10)
2стм, (А/у )2
Фм, - Фт =-----°Ц) ^ 2 , (Мг ) 6 ( И Ь , (11)
стм,. (А/у)
где фР ~ ф(Р), фМ1 ~ ф(М1), фМ2 = ф(М2), фМ3 ~ ф(М3) - значения сеточной функции
электрической проводимости, А/ц - длина ребра Л1Л1, Апу - расстояние от точки Р
до ребра ЛЛу, Ап1 - расстояние от точки М, до ребра Л Лу, Зу = Яалл/, Яу = $алам
(рис. 3). Кроме того, введены следующие обозначения:
* ст(Р) * ст(М,)
°р =, ^м, = -Ц^ , 1 = 1, 2, 3.
Апу ' Апу
При аппроксимации граничных условий (2) - (4) (формулы (9) - (11)) использовались фиктивные конечные объемы уИ, которые симметрично достраивались к треугольным ячейкам сетки, две вершины которых лежат на границе Г.
В ходе решения задачи теоретически были исследованы характеристики полученной разностной схемы. Погрешность аппроксимации разностной задачи (8) -(11) с учётом граничных условий имеет первый порядок точности О(Дп), Ап = шах(Ап у, Ап* }. С помощью принципа максимума и мажоранты Гершгорина
была доказана устойчивость разностной схемы (8) - (11), что в совокупности обеспечивает сходимость решения разностной задачи (8) - (11) к точному решению дифференциальной задачи (1) - (4) с первым порядком точности.
Записывая полученную разностную схему для каждой ячейки сетки на некоторой триангуляции области Б , рассмотренный подход к решению прямой задачи ЭИТ приводит к большой системе линейных алгебраических уравнений (СЛАУ). Поскольку использована треугольная неструктурированная сетка, то матрица системы является разреженной: в каждой её строке содержится по четыре ненулевых элемента. В расположении элементов в строке и матрице в целом не наблюдается чёткой структуры, кроме того, что один элемент строки всегда находится на главной диагонали матрицы. Причём для внутренних узлов сетки этот коэффициент равен сумме остальных ненулевых положительных элементов, стоящих в данной строке, взятых с обратным знаком. Для решения составленной несимметричной
СЛАУ использован стабилизированный метод бисопряжённых градиентов Б1-Св8ТЛБ [15]. Тестовые расчёты задачи (1) - (4) выполнены для области Б (рис. 1), на которой вводились различные расчётные сетки.
Алгоритм реконструкции
Обратная задача ЭИТ является некорректно поставленной задачей, поэтому для восстановления требуемого количества неизвестных коэффициентов электрической проводимости с нужно использовать информацию большого набора независимых измерений. Согласно методике ЭИТ совокупность измерений электрических напряжений на границе объекта накапливается путём многократного пропускания электрического тока через него. Схема измерений учитывает все возможные пары электродов, при этом на каждой паре перебирается несколько образцов инжектируемых токов, т.е. используются токи различной силы. Собранные данные о разности потенциалов на электродах обеспечивают обратную задачу ЭИТ необходимыми дополнительными условиями. Как правило, достаточное количество априорной информации сужает класс допустимых решений задачи, в результате чего она становится условно корректной [3]. Учитывая вышесказанное, обратная задача ЭИТ сводится к краевой задаче (1) - (4) с дополнительными данными в качестве совокупности измерений электрических напряжений на границе объекта (5) - (6) для подбора с°р4(х,у), обеспечивающего минимум функционала (7) с заданной точностью.
В работе предложена итерационная процедура для решения рассматриваемой коэффициентной обратной задачи, в основу которой положен алгоритм дифференциальной эволюции [6, 7]. Эвристический алгоритм глобальной оптимизации на непрерывных пространствах, известный как метод дифференциальной эволюции, был предложен Р. Сторном и К. Прайсом [6] в качестве нового эволюционного подхода к минимизации различных нелинейных и недифференцируемых функций из непрерывного пространства. Согласно данному алгоритму, оценка параметра электрической проводимости с осуществляется путём перебора различных вариантов решений прямой задачи (1) - (4). Подробнее, дифференциальная эволюция оперирует популяцией постоянного размера №, состоящей из вещественнозначных векторов, которые являются возможными вариантами решения задачи. Они выбираются из множества допустимых решений Е с Яп . Для анализа вариантов алгоритм использует различия между индивидами популяции. Поиск «оптимального» решения осуществляется посредством циклического изменения текущей популяции по эволюционным правилам, таким как мутация, скрещивание и селекция. Внутренним циклом алгоритма анализируются все векторы текущей популяции, принадлежащие поколению О, они по очереди выбираются вектором-мишенью. На каждой итерации формируется новый вектор по вектору-мишени и произвольно выбранным векторам рассматриваемого поколения, который впоследствии служит образцом для оценки вектора-мишени. Сравнивая вектор-мишень с вектором-образцом по значению целевой функции, определяется, какой из них перейдёт в новую популяцию (поколение О + 1). Итерационный процесс алгоритма завершается по одному из двух возможных исходов: будет найден индивид популяции, для которого целевая функция принимает значение меньшее некоторого фиксированного уровня, или число поколений достигнет заданного максимального количества Ошах. В последнем случае из полученной популяции
выбирается вектор, которому соответствует наименьшее значение целевой функции, и он считается приближенным решением задачи.
В задаче ЭИТ целевая функция для каждого вектора популяции представляется в следующем виде:
где с1,0 - вектор распределения электрической проводимости внутри объекта, , = 1,...,ЫР, ЫР- размер популяции, О- индекс поколения, которому принадлежит популяция. Афг(сг’°) - значения разности электрических потенциалов, вычисленные на границе области при решении прямых задач ЭИТ с пробными значениями с. Требуется найти глобальный условный минимум целевой функции / (с) на множестве Е.
Основные этапы алгоритма дифференциальной эволюции для решения обратной задачи ЭИТ [6, 7]:
1) Инициализация. Поколение О = 0. На данном этапе с помощью датчика случайных чисел генерируется начальная популяция, состоящая из № п-мерных векторов вида
где п - количество неизвестных параметров электрической проводимости внутри области (например, для рис. 1 п = 2); гаМк[0,1] обозначает равномерно распределенное случайное число из диапазона [0,1].
Для вычисления целевой функции векторов популяции решаются прямые задачи со сформированными значениями с1,0, і = 1,...,ЛР, для всех конфигураций инжектируемых токов на различных парах электродов.
2) Мутация. Цикл ] = 1,...ДР. В текущей популяции последовательно перебираются все её представители с целью проверки их конкурентоспособности по сравнению с другими векторами, т.е. все элементы один раз выполняют функцию вектора-мишени с1 = <С’°. Для каждого проверяемого вектора с4 формируется локальный вектор се согласно правилу
телей популяции, причём ни один из них не равен вектору-мишени с, и все они попарно различны, а,Ь,с 6 {1,2,...,ЫР}.
3) Скрещивание. Компоненты мутировавшего вектора се смешиваются с компонентами вектора-мишени с‘, чтобы получить вектор-образец с5. Он формируется по правилу
где СЯ - константа скрещивания, ик 6 [0,1] - равномерно распределенное случайное число, гпЬг (у) 6 {1,2,..., п} - произвольно выбранный индекс, который гаран-
Векторы са,а, сь,а, сс,а случайным образом выбираются из оставшихся представи-
аек, если ик < СЯ или к = гпЬг (}), ак, если ик > СЯ и к Ф гпЬг(}), к = 1,...,п,
тирует, что вектор с5 всегда будет получать, по крайней мере, один компонент вектора се.
4) Селекция. На данном этапе решается совокупность прямых задач с параметрами вектора с и находятся значения Афр(с5) = фои4,Р(с5) - фшр(сО, р = 1,.,Ь. Если вектор-образец с5 приводит к меньшему значению целевой функции /(с), чем вектор-мишень с‘, то вектор-образец заменяет проверяемый вектор с‘ в следующем поколении; иначе, сохраняется вектор с(.
5) О = О + 1 и выполняется проверка условий. Итерационный процесс завершается, если значение целевой функции / (сЬе!4,0) от лучшего вектора популяции сЬе8‘°, принадлежащего поколению О, меньше фиксированного ограничения е, е > 0, или число поколений достигло заданного максимального количества Отах, иначе алгоритм переходит к пункту 2.
Обобщая вышесказанное, работа описанного алгоритма в применении к задаче ЭИТ опирается на два действия: решение большого количества прямых задач ЭИТ на каждой итерации и минимизации квадратичного функционала / от разности между вычисленным и измеренным электрическим напряжением на границе объекта.
Отметим, что минимизация методом дифференциальной эволюции регулируется несколькими параметрами, а именно, ЫР, СЯ и ЕО. Они влияют на скорость сходимости и работоспособность процесса поиска. Параметр ЫР выбирается от 5п до 10п [6]. Чтобы выполнялось условие пункта 2, ЫР должен быть равен, по крайней мере, 4. Константа скрещивания СЯ 6 [0,1] определяется пользователем. Большие значения СЯ (СЯ = 0,9 или СЯ = 1,0) обычно ускоряют сходимость алгоритма. Параметр ЕО 6 (0, 2], Е = 0,5 является хорошим начальным выбором. Если популяция сходится преждевременно, то ЕО и/или ЫР следует увеличивать. Значения Е меньшие 0,4, а также большие 1 лишь изредка оказываются эффективными [6, 7].
Результаты расчётов
Для проверки работоспособности описанного алгоритма проведены численные эксперименты на модели, включающей объект с одной областью неоднородности и четырьмя электродами (2 = 4), прикреплёнными на его поверхность (рис. 2). Всего рассмотрено 4 токовые конфигурации (М = 4) и сформировано семейство
измерений {р, Фттей, Фтигригей}, Р = 1,. .,Ь (Ь = 24), для каждой пары электродов, для чего были решены прямые задачи со следующими значениями электрической проводимости: с1 = 202,4, с2 = 0,5 (внутреннее включение). Расчёты проводились с двойной точностью. Результаты решения прямых задач с разным расположением электродов на границе для случая кусочно-однородной среды приведены на рис. 4.
Из рисунка видно, что разность потенциалов электрического поля напрямую зависит от пути прохождения электрического тока, его силы и значения электрической проводимости внутренних составляющих объекта. В силу того, что разнообразные вещества характеризуются различными значениями электрической проводимости, анализируя её распределение внутри объекта, можно судить о его структурных особенностях и составе. Благодаря существенному различию электрической проводимости биологических тканей, ЭИТ позволяет реконструировать контрастные изображения.
Феи -Ф«, = 0,0155
Рис. 4. Примеры распределения значений электрического потенциала в области Б (рис. 2) при инжектировании тока через различные пары электродов, /;п = 1,5 мА
При решении обратной задачи в соответствии с рекомендациями статьи [6, 7] размер популяций ЫР выбран равным от 5п до 10п, где п = 2 - количество неизвестных параметров электрической проводимости, константа скрещивания СЯ е [0, 1], константа мутации Еа е (0,2]. Начальное распределение электрической проводимости задавалось по формуле
о/0 = 0,001 + гаМк[0, 1] 1000, к = 1,2, I = 1,...,ЫР.
При таких условиях после инициализации дифференциальной эволюции целевая функция имеет значение порядка 1.
При решении рассматриваемой в данной работе обратной задачи методом дифференциальной эволюции применялся также другой вариант операции мутации в п.2:
ое = 0Ье8‘,е + ¥0(оъ’° - ос,°), (13)
в котором использовался лучший на каждой итерации вектор значений электрической проводимости аЬе84,°. Расчеты, проведенные для двух вариантов операции мутации (формулы (12) и (13)), продемонстрировали существенное преимущество (13), поскольку этот вариант значительно ускоряет итерационный процесс поиска с заданной точностью.
Для представленных на рис. 2 области и покрывающей ее сетки в 130 треугольных ячеек были проведены расчеты с различными значениями параметров
метода дифференциальной эволюции (ИР, СЯ, Ро, Отах = 500 или 1000, е > 0 - допустимый минимум целевой функции, при достижении которого поиск прекращался). Цель вычислительного эксперимента заключалась в определении оптимальных значений параметров метода, обеспечивающих быструю сходимость процесса поиска к глобальному минимуму целевой функции и хорошую точность определения истинных значений электрической проводимости. Результаты выполненных расчетов приведены в таблице. Оценка быстроты и качества поиска была получена на основе осреднения результатов 10-ти запусков программы решения обратной задачи при одних и тех же значениях параметров метода.
Расчёты параметра электрической проводимости
№ е Кол-во итераций (медиана) ИР Ро СЯ Вычисленная электрическая проводимость Вычисленная электрическая проводимость 02 Значение целевой функции .Г
1 10-6 10 20 0,9 0,4 182,728 13,774 5,8-10-7
2 10-7 18 20 0,9 0,4 198,927 2,687 4,75-10-8
3 10-8 34 20 0,9 0,4 200,959 1,112 3,89-10-9
4 10-9 39 20 0,9 0,4 201,965 0,646 3,88-10-10
5 10-ю 48 20 0,9 0,4 202,364 0,538 5,71-10-11
6 10-11 66 20 0,9 0,4 202,188 0,554 7,17-10-12
7 10-12 500 20 0,9 0,4 202,312 0,497 3,57-10-12
8 10-11 187 20 0,9 0,1 202,338 0,505 4,25-10-12
9 10-10 36 20 0,7 0,4 202,203 0,567 4,57-10-11
10 10-11 56 20 0,7 0,4 202,338 0,51 6,58-10-12
11 10-12 500 20 0,7 0,4 202,37 0,506 1,95-10-12
12 10-10 33 20 гаш1о(0,2] 0,4 202,13 0,569 5,35-10-11
13 10-11 47 20 гаш1о(0,2] 0,4 202,163 0,574 5,4910-12
14 10-12 229 20 гаиёо(0,2] 0,4 202,305 0,544 9,17-10-13
15 10-10 22 20 0,7 0,9 202,241 0,543 6,33-10-11
16 10-11 32 20 0,7 0,9 202,292 0,561 8,58-10-12
17 10-12 500 20 0,7 0,9 202,141 0,587 5-10-12
18 10-10 500 4 0,7 0,4 206,853 114,9 1,37-10-5
19 10-10 98 10 0,7 0,4 197,149 13,694 2,52-10-7
20 10-10 36 30 0,7 0,4 202,134 0,592 4,56-10-11
В таблице варианты № 1-7 представляют результаты исследования влияния точности определения минимума целевой функции е. Видно, что при недостаточно малом е метод не успевает сойтись (по точности воспроизведения истинных значений электрической проводимости), так как значения целевой функции достигают заданного уровня раньше, чем найдено приемлемое решение. Кроме того, следует отметить, что при е = 10-12 увеличение максимального числа итераций метода дифференциальной эволюции с 500 до 1000 не вносит существенного улучшения в качество решения обратной задачи.
Варианты расчетов № 5-7, 9-14 характеризуют влияние выбора параметра Ро, который определяет темп изменения пробных значений электрической проводимости на этапе мутации (п.2 метода дифференциальной эволюции). Из таблицы видно, что при ИР = 20, Отах = 500, СЯ = 0,4 наименьшие затраты на получение решения обратной задачи достигаются в том случае, когда параметр Ро есть случайная величина, равномерно распределенная в диапазоне (0,2] и заново генери-
руемая с помощью датчика псевдослучайных чисел на каждой итерации метода
G < Gmax.
Варианты расчетов № 6, 8, 9-11 и 15-17 (см. таблицу) направлены на исследование влияния значения константы скрещивания 0 < CR < 1 на количество итераций метода и точность решения обратной задачи. Анализирую представленные в таблице числовые значения, можно отметить, что большие значения параметра CR, например CR = 0,9, ускоряют сходимость метода. Наоборот, значения CR, близкие к 0, замедляют процесс поиска, но в итоге приводят к более точным результатам определения истинных значений сь с2. Наиболее оптимальной по скорости сходимости метода и точности решения обратной задачи является комбинация Fg = 0,7 и CR = 0,4 (вариант № 10).
Варианты расчетов № 9, 18-20 показывают влияние размера выбранной популяции NP при Fg = 0,7, CR = 0,4, е = 10-10 на поведение метода дифференциальной эволюции в процессе поиска минимума целевой функции. Значение NP менялось в диапазоне 2n + 15n. Представленные в таблице результаты указывают на то, что при небольшом размере популяции (NP < 5n) метод не позволяет с требуемым качеством найти значения электрической проводимости. При NP > 10n метод сходится практически за одно и то же количество итераций с близкой точностью.
Заключение
В данной работе предложен численный метод решения задач ЭИТ на основе алгоритма дифференциальной эволюции, который был реализован для оценки значений электрической проводимости на двумерной модели области, включающей одну неоднородность. В результате исследования прямой задачи ЭИТ установлено, что электрические свойства сред различной природы (например, биологических тканей), путь прохождения электрического тока и его сила напрямую влияют на значение и распределение разности электрических потенциалов внутри изучаемого объекта. Проведенные по предложенному методу решения обратной задачи ЭИТ численные эксперименты продемонстрировали работоспособность метода на крупной расчётной сетке. Два параметра электрической проводимости восстанавливались с использованием стохастического алгоритма дифференциальной эволюции, были опробованы различные вариации управляющих параметров метода. После серии тестовых запусков наилучшие результаты показала комбинация: CR = 0,4, NP = 10n и Fg = 0,7.
Подход к реконструкции ЭИТ, включающий, в частности, многократное решение прямых задач с различными параметрами, и структура алгоритма дифференциальной эволюции позволяют говорить о возможности эффективного распараллеливания на многопроцессорных системах. Параллельная реализация позволит применять метод на более детальных сетках с меньшими временными затратами на вычисления и решать задачу с большим числом неизвестных параметров электрической проводимости.
ЛИТЕРАТУРА
1. Electrical impedance tomography [Электронный ресурс]. Электрон. дан. URL: http:// www.eit.org.uk/ (дата обращения: 15.02.2011).
2. Электроимпедансная томография (ЭИТ) [Электронный ресурс]. Электрон. дан. URL: http://www.cplire.ru/mac/etomo/index.html (дата обращения: 15.02.2011).
3. Borcea L. Electrical impedance tomography // Inverse Problems. 2002. № 18. P. 99-136.
4. Lionheart W., Polydorides N., Borsic A. Electrical Impedance Tomography: Methods, History and Applications. Manchester, 2004. P. 62.
5. Пеккер Я.С., Бразовский К.С. Электроимпедансная томография. Томск: Изд-во НТЛ, 2004. 192 с.
6. Storn R. Differential evolution - a simple and efficient heuristic for global optimization over continuous spaces // J. Global Optimization. 1997. No. 11. P. 341-356.
7. Li Y., Xu G., Guo L., Wang L. Resistivity parameters estimation based on 2D real head model using improved differential evolution algorithm // Engineering in Medicine and Biology Society, 28th Annual International Conference of the IEEE. New York, 2006. P. 6720-6723.
8. Isaacson D. Problems in impedance imaging // Lecture Notes in Physics. 1993. V. 422/1993. P. 62-70.
9. Седов Л.И. Механика сплошной среды. М.: Наука, 1983. Т. 1. 528 с.
10. Gmsh: a three-dimensional finite element mesh generator with built-in pre- and postprocessing facilities [Электронный ресурс]. Электрон. дан. URL: http://geuz.org/gmsh/ (дата обращения: 15.02.2011).
11. ANSYS - Simulation Driven Product Development [Электронный ресурс]. Электрон. дан. URL: http://www.ansys.com/ (дата обращения: 15.02.2011).
12. Патанкар С. Решение задач теплообмена и динамики жидкости: пер. с англ. В.Д. Виленский. М.: Энергоатомиздат, 1984. 124 с.
13. Елизарова Т.Г. Квазигазодинамические уравнения и методы расчета вязких течений. М.: Научный Мир, 2007. 350 с.
14. Шерина Е.С. Численный метод решения прямой задачи электроимпедансной томографии // Перспективы развития фундаментальных наук: труды VI Международной конференции студентов и молодых ученых. Томск, 2009.
15. Van der Vorst H. Bi-CGSTAB: A fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems // SIAM J. Sci. Statist. Comput. 1992. V. 13. P. 631644.
Статья поступила 14.04.2012 г.
Sherina E.S., Starchenko A.V. NUMERICAL METHOD FOR RECONSTRUCTING THE ELECTRICAL IMPEDANCE DISTRIBUTION IN BIOLOGICAL OBJECTS USING CURRENT MEASUREMENTS AT THE BOUNDARY. Using boundary measurements of electric current and voltage, electrical impedance tomography (EIT) is applied for reconstructing the unknown electrical conductivity distribution within living objects with a heterogeneous structure. This paper describes the numerical approach to solve EIT inverse problems in terms of the Differential Evolution algorithm.
Keywords: electrical impedance tomography (EIT), unstructured triangular mesh, difference scheme, coefficient inverse problem, genetic optimization algorithms, Differential Evolution.
SHERINA Ekaterina Sergeevna (Tomsk State University)
E-mail: [email protected]
STARCHENKO Alexander Vasil’evich (Tomsk State University)
E-mail: [email protected]
