Numerical method for reconstructing the electrical impedance distribution in biological objects using current measurements at the boundary.pdf В настоящее время в медицине и промышленности для компьютерной диагностики применяется способ получения томографического изображения, называемый электроимпедансной томографией (ЭИТ) или томографией приложенных потенциалов [1-5]. Он основан на использовании электрического тока в качестве средства, зондирующего исследуемый объект. Цель ЭИТ состоит в том, чтобы восстанавливать неизвестное распределение электрической проводимости внутри объекта с неоднородной структурой, используя электрические измерения только на границе объекта. Слабый электрический ток высокой частоты подводится к объекту исследования посредством системы токопроводящих электродов, прикрепляемых к его поверхности. На них выполняются измерения электрического напряжения или тока с целью получения исходной совокупности данных для реконструкции. С физической точки зрения, ЭИТ как диагностический метод с некоторой погрешностью воспроизводит распределение электрической проводимости или электрического сопротивления (импеданса) внутри объекта. Поскольку различные вещества характеризуются своими значениями электрической проводимости (сопротивления), знание характера распределения электрической проводимости позволяет установить внутреннюю структуру объекта, что может быть использовано в медицине и промышленности [3]. Однако восстановление распределения скалярной функции проводимости внутри области по измерениям электрического тока на границе является некорректной задачей и требует применения специальных методов реконструкции, включающих решение обратной задачи [3, 4]. В данной работе предложен метод решения задач ЭИТ для исследования внутренней структуры биологических объектов, основанный на алгоритме дифференциальной эволюции [6, 7]. Использованный стохастический метод - дифференциальная эволюция - относится к классу эвристических эволюционных алгоритмов, применяемых для решения задач оптимизации. Подобные алгоритмы итеративно моделируют эволюционный процесс, в котором в качестве популяции рассматривается множество допустимых решений задачи. В процессе моделирования учитываются основные механизмы теории биологической эволюции, такие, как мутация, скрещивание и селекция. Применение генетических алгоритмов представляет интерес для решения обратных задач, поскольку они позволяют обойти главную трудность задач этого класса, а именно, характерную некорректность их постановки. В случае задач ЭИТ их некорректный характер вызван нарушением устойчивости процесса получения решения, которое возникает из-за недостаточной точности и количества входных данных. В предложенном методе алгоритм оперирует популяцией векторов электрической проводимости, среди которых может оказаться возможное решение задачи. Поиск «оптимального» решения осуществляется путём циклического изменения текущей популяции по законам эволюции и оценки значения целевой функции для каждого вектора популяции. По завершению работы алгоритма выбирается представитель популяции, для которого целевая функция принимает наименьшее значение, и данный вектор считается приближенным решением задачи. Метод дифференциальной эволюции был выбран в качестве способа минимизации в силу его простоты, высокой скорости работы и прямолинейной стратегии поиска. Принимая во внимание результаты, полученные в работе [6], можно сделать вывод, что метод дифференциальной эволюции в большинстве случаев превосходит все другие методы оптимизации по требуемому количеству оценок функции, необходимых для определения минимума целевой функции с заданной точностью. Алгоритм легко использовать, поскольку он прост и имеет небольшое число управляющих параметров, которые можно выбрать из заранее определенных числовых интервалов. Постановка задачи Разработанный метод решения обратных задач ЭИТ реализован на примере двумерной области D = D и Г (рис. 1), которая представляет собой схематическое изображение некоторого биологического объекта, имеющего неоднородную структуру. Границы неоднородностей внутри области заранее известны. В каждой подобласти коэффициент электрической проводимости с имеет некоторое постоянное значение, поэтому ст(x, y) можно рассматривать как кусочно-постоянную функцию от координат (x, y) е D . Необходимо определить значения коэффициента с внутри области, если априори известны измерения электрического напряжения и силы тока на границе Г. Чтобы получить исходные данные для поставленной задачи, через электропроводящий объект, который находится в воздушной среде, пропускается электрический ток небольшой силы 1-5 мА высокой частоты 100-150 кГц. Ток инжектируется через расположенные локально на его поверхности электроды, на которых выполняются измерения разности электрических потенциалов [3, 4]. Между электродами граница Г контактирует с воздухом. Использована модель ЭИТ с идеально проводящими электродами, называемая в научной литературе shunt model [3, 4, 8], и реализована одна из схем организации измерений, в которой электроды включаются попарно и электрическое напряжение определяется на активной паре. Обращаясь к физической постановке задачи, предполагаем, что учитываемые моделью процессы описываются электрическим и магнитным полями, которые удовлетворяют стационарным уравнениям Максвелла. В случае применения ЭИТ к объектам биологической природы интенсивность магнитного поля становится пренебрежимо малой из-за особенности в проводимости биологических тканей [3, 4], и поэтому магнитная компонента не будет рассматриваться. При таких условиях вектор напряженности электрического поля равен градиенту электрического потенциала ф, распределение которого внутри объекта исследования будет удовлетворять уравнению эллиптического типа [9]. В силу отсутствия внутри объекта источников электрического поля сумма втекающих Im и вытекающих Iout токов равняется нулю согласно закону сохранения заряда, Iin + Iout = 0. Также предполагаем, что на границе, где инжектируется ток, известно значение потенциала электрического тока фт. Рис. 1. Пример модели биологического объекта D с границей Г, на которую прикреплены два электрода с площадью поверхности Sin = Sout. В объект инжектируется электрический ток величиной Iin in С учётом описанной выше физической постановки математическая модель ЭИТ представляет собой краевую задачу для уравнения эллиптического типа в частных производных [3, 4]: (1) ^Z^XA = 0, (х,у) е Г (2) (3) ф(х у ^ф^ (х у )еГт; дф(х,у) О (X у) Sout-I-= -^out = huU ^ у ) е Г, (4) где о - коэффициент электрической проводимости, ф - потенциал электрического поля, n - единичный вектор внешней нормали, Jn - плотность тока на границе, Sin, Sout - площадь контактной поверхности in- и out-электродов соответственно. Уравнение (1) показывает зависимость скалярных функций электрического потенциала ф(х,у) и проводимости c(x,y) внутри объекта. Краевое условие (2) выражает контакт тела с воздухом, (3) и (4) - контакт тела с электродами. При некотором априори заданном распределении значений функции проводимости с и значениям электрического потенциала фш и электрического тока Iin на границе из краевой задачи (1) - (4) можно найти неизвестное распределение потенциалов ф электрического поля внутри области D . В общем случае эта задача не может быть решена аналитически [3, 4]. На основе постановки (1) - (4) можно сформулировать задачу по восстановлению неизвестного распределения электрической проводимости с(х,у) в области D, если заданы значения силы инжектируемого тока Iin,p (p = 1,...,L) и измерены значения разности электрического потенциала Лфртешшге 0, или число поколений достигло заданного максимального количества Gmax, иначе алгоритм переходит к пункту 2. Обобщая вышесказанное, работа описанного алгоритма в применении к задаче ЭИТ опирается на два действия: решение большого количества прямых задач ЭИТ на каждой итерации и минимизации квадратичного функционала f от разности между вычисленным и измеренным электрическим напряжением на границе объекта. Отметим, что минимизация методом дифференциальной эволюции регулируется несколькими параметрами, а именно, NP, CR и FG. Они влияют на скорость сходимости и работоспособность процесса поиска. Параметр NP выбирается от 5п до 10п [6]. Чтобы выполнялось условие пункта 2, NP должен быть равен, по крайней мере, 4. Константа скрещивания CR е [0,1] определяется пользователем. Большие значения CR (CR = 0,9 или CR = 1,0) обычно ускоряют сходимость алгоритма. Параметр FG е (0, 2], F = 0,5 является хорошим начальным выбором. Если популяция сходится преждевременно, то FG и/или NP следует увеличивать. Значения F меньшие 0,4, а также большие 1 лишь изредка оказываются эффективными [6, 7]. Результаты расчётов Для проверки работоспособности описанного алгоритма проведены численные эксперименты на модели, включающей объект с одной областью неоднородности и четырьмя электродами (Q = 4), прикреплёнными на его поверхность (рис. 2). Всего рассмотрено 4 токовые конфигурации (M = 4) и сформировано семейство измерений {lin,p, Ф™ей, Ф™ей}, p = 1,...,L (L = 24), для каждой пары электродов, для чего были решены прямые задачи со следующими значениями электрической проводимости: о1 = 202,4, о2 = 0,5 (внутреннее включение). Расчёты проводились с двойной точностью. Результаты решения прямых задач с разным расположением электродов на границе для случая кусочно-однородной среды приведены на рис. 4. Из рисунка видно, что разность потенциалов электрического поля напрямую зависит от пути прохождения электрического тока, его силы и значения электрической проводимости внутренних составляющих объекта. В силу того, что разнообразные вещества характеризуются различными значениями электрической проводимости, анализируя её распределение внутри объекта, можно судить о его структурных особенностях и составе. Благодаря существенному различию электрической проводимости биологических тканей, ЭИТ позволяет реконструировать контрастные изображения. При решении обратной задачи в соответствии с рекомендациями статьи [6, 7] размер популяций NP выбран равным от 5n до 10n, где n = 2 - количество неизвестных параметров электрической проводимости, константа скрещивания CR е [0, 1], константа мутации FG е (0,2]. Начальное распределение электрической проводимости задавалось по формуле О^0 = 0,001 + randk[0, 1] 1000, k = 1,2, i = 1,...,NP. При таких условиях после инициализации дифференциальной эволюции целевая функция имеет значение порядка 1. При решении рассматриваемой в данной работе обратной задачи методом дифференциальной эволюции применялся также другой вариант операции мутации в п.2: О = 0bestG + FG(ob,G - gc,°), (13) в котором использовался лучший на каждой итерации вектор значений электрической проводимости obest,G. Расчеты, проведенные для двух вариантов операции мутации (формулы (12) и (13)), продемонстрировали существенное преимущество (13), поскольку этот вариант значительно ускоряет итерационный процесс поиска с заданной точностью. Для представленных на рис. 2 области и покрывающей ее сетки в 130 треугольных ячеек были проведены расчеты с различными значениями параметров метода дифференциальной эволюции (NP, CR, FG, Gmax = 500 или 1000, е > 0 - допустимый минимум целевой функции, при достижении которого поиск прекращался). Цель вычислительного эксперимента заключалась в определении оптимальных значений параметров метода, обеспечивающих быструю сходимость процесса поиска к глобальному минимуму целевой функции и хорошую точность определения истинных значений электрической проводимости. Результаты выполненных расчетов приведены в таблице. Оценка быстроты и качества поиска была получена на основе осреднения результатов 10-ти запусков программы решения обратной задачи при одних и тех же значениях параметров метода. out 9out -Ф„ = 0,0135 9out -Ф„ = 0,0155 lout Фout -Ф,п = 0,022 Рис. 4. Примеры распределения значений электрического потенциала в области D (рис. 2) при инжектировании тока через различные пары электродов, /in = 1,5 мА Фои -Ф„ = 0,01 Расчёты параметра электрической проводимости № е Кол-во итераций (медиана) NP Fg CR Вычисленная электрическая проводимость Вычисленная электрическая проводимость Значение целевой функции f 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-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 randG(0,2] 0,4 202,13 0,569 5,35-10-11 13 10-11 47 20 randG(0,2] 0,4 202,163 0,574 5,49-10-12 14 10-12 229 20 randG(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 характеризуют влияние выбора параметра FG, который определяет темп изменения пробных значений электрической проводимости на этапе мутации (п.2 метода дифференциальной эволюции). Из таблицы видно, что при NP = 20, Gmax = 500, CR = 0,4 наименьшие затраты на получение решения обратной задачи достигаются в том случае, когда параметр FG есть случайная величина, равномерно распределенная в диапазоне (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. Подход к реконструкции ЭИТ, включающий, в частности, многократное решение прямых задач с различными параметрами, и структура алгоритма дифференциальной эволюции позволяют говорить о возможности эффективного распараллеливания на многопроцессорных системах. Параллельная реализация позволит применять метод на более детальных сетках с меньшими временными затратами на вычисления и решать задачу с большим числом неизвестных параметров электрической проводимости.

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