CC BY
69
Радіотехнічні кола та сигнали
РАДІОТЕХНІЧНІ КОЛА ТА СИГНАЛИ
УДК 621.372.061
РЕШЕНИЕ ПРЯМОЙ ЗАДАЧИ ИМПЕДАНСНОЙ ТОМОГРАФИИ МЕТОДАМИ ТЕОРИИ ЦЕПЕЙ
И.А. Рыбина, магистрант
Национальный технический университет Украины "Киевский политехнический институт", г. Киев, Украина
Томографические методы исследований (неразрушающего контроля) последние 30 лет (с момента появления промышленных рентгеновских томографов) находят все более широкое распространение в радиотехнике, медицине, механике и т.д. [1—7]. Среди методов томографических исследований импедансная томография отличается относительной простотой технического оборудования и значительной сложностью требуемого для получения образов математического аппарата.
Математическую задачу импедансной томографии можно разбить на две основных — задачу анализа (вычисление напряжений или передаточных проводимостей по обводу томографического сечения - фантома — при известных поверхностных или объемных проводимостях тканей) и задачу синтеза (восстановление по внешним измерениям - проекциям — внутреннего распределения поверхностных или объемных проводимостей фантома).
Сложность математической задачи как анализа, так и синтеза заключается (в отличие от аналогичных задач рентгеновской, гамма, ядерного магнитного резонанса, ультразвуковой томографии) в том, что токи между электродами (сигналы - носители информации) измерительного от источника тока протекают по всей плоскости (и объему) фантома, а не по прямой линии между электродами. В остальных вышеуказанных томографических системах сигнал — носитель информации распространяется прямолинейно («линейная» томография) и на него не влияют ткани, расположенные в фантоме вне этой прямой линии.
Вследствие указанных причин реконструкция образа (синтез) представляет собой сложную задачу (в отличие от «линейной» томографии) и решается итерационными методами.
В силу итерационного характера обратной задачи импедансной томографии особое значение приобретает способ решения прямой задачи, от точности и быстродействия решения которой зависит успех итерационной процедуры [8].
4
Вісник Національного технічного університету України "КПІ" Серія — Радіотехніка. Радіоапаратобудування.-2010.-№43
Радіотехнічні кола та сигнали
Для решения прямой задачи, как правило, используют дискретизацию фантома (метод конечных элементов [8, 9]), при которой вся плоскость (объем) фантома разбивают на ячейки с постоянной в их пределах поверхностной (объемной) удельной проводимостью (первая дискретизация).
Одним из направлений решения задачи анализа является использование методов теории цепей. При этом модель фантома представлена линейной пассивной электрической цепью, которая на частотах тока источника f < 100 кГц имеет резистивный характер, а на частотах f> 500 кГц содержит комплексные сопротивления. При этом каждый конечный элемент следует заменить его эквивалентной электрической схемой, причем такие схемы подключаются одна к другой уже своими узлами (вторая дискретизация).
Для определенности в дальнейшем рассмотрим особенности предлагаемого алгоритма решения прямой задачи на примере плоского фантома.
Для этого выберем в качестве конечного элемента квадрат со стороной а, электрической моделью которого является четырехугольник с двумя диагональными ветвями, полученный из четырех равнобедренных прямоугольных треугольников с гипотенузами, равными а. Размеры всех четырехугольников одинаковы, а веса проводимостей резистивных ветвей пропорциональны значениям удельных поверхностных проводимостей конечных элементов.
При таком представлении каждый конечный элемент является пикселем изображения (уровень серого или цвет соответствует значению проводимости) внутренней структуры фантома, а задача анализа заключается в замене таких пикселей их электрическими эквивалентными схемами с дальнейшим составлением системы уравнений узловых напряжений фантома
_ YV = I, _ (1)
где Y - матрица проводимостей порядка N; I - столбец независимых источников тока /о размера /V/1, причем, если источник подключен к первому и общему узлам эквивалентной схемы, то столбец 7 имеет лишь один
ненулевой элемент /0 в первой строке; U - вектор-столбец искомых узловых напряжений размера N*1, причем для получения проекции (напряжений на измерительных электродах) необходимыми являются лишь 6, 14 или 30 (в зависимости от количества электродов - 8, 16 или32) на внешнем обводе фантома.
При этом визуальное качество получаемого в результате реконструкции изображения зависит от числа M пикселей (конечных элементов) фантома. В зависимости от числа пикселей различают «грубую» и «тонкую» модели фантома. Так, для квадратного фантома «грубая» модель содержит
Вісник Національного технічного університету України "КПІ" Серія — Радіотехніка. Радіоапаратобудування.-2010.-№43
5
Радіотехнічні кола та сигнали
M=16*16=256, а «тонкая» - M=32*32=1024 пикселя. Ясно, что как для «грубой», так и для «тонкой» моделей визуальное качество образа является неудовлетворительным (каждая строка фантома содержит ha на эквивалентной лишь 16 или 32 пикселя соответственно). При этом число независимых узлов N эквивалентной схемы составляет для квадратного фантома для «грубой» модели N=288, а для «тонкой» N=1198, т.е. таковым будет порядок системы уравнений (1).
Естественное желание повысить визуальное качество образа приводит к увеличению числа пикселей. Так при формате изображения M=256x256=65536 порядок N системы (1) составит N=65048, а для M=1024x1024=1048576 порядок N =1050624. В случае круглого фантома соотношение между количеством пикселей фантома и порядком системы (1) поменяется, однако для M>250 порядок системы уравнений N мало отличается от числа пикселей. Таким образом, для решения прямой задачи необходимо решить систему уравнений равновесия фантома высокого порядка. Однако наличие жесткой (заданной формой конечных элементов) структуры фантома позволяет значительно упростить вычисления.
Рассмотрим в качестве примера фрагмент фантома (рис. 1 а). Из рис. 1 а видно, что в системе уравнений (1) строка а матрицы Y содержит лишь 9 ненулевых элементов не зависимо от порядка матрицы N (на рис. 1 а показаны те диагональные элементы ba, da, fa, ha схеме конечного элемента рис. 1 б, которые учитываются при составлении a-то уравнения в (1)). Таким образом, матрица проводимостей Y будет сильно разреженной, поскольку в каждой ее строке (столбце) количество ненулевых элементов не превышает девяти.
Рис.1
Решение систем линейных уравнений высокого порядка с разреженными матрицами коэффициентов является одной из актуальных задач теории цепей [11—14]. Следует выделить два подхода к решению разрежен-
6
Вісник Національного технічного університету України "КПІ" Серія — Радіотехніка. Радіоапаратобудування.-2010.-№43
Радіотехнічні кола та сигнали
ной системы уравнений: приведение разреженной матрицы к ленточному виду [13] и метод модификаций [14—17].
Приведение матрицы Y к ленточному (в данном случае — девятидиагональному) виду позволяет радикально сократить число операций при вычислении решения U. Однако, при решении обратной задачи большое значение имеет быстрое и экономное вычисление производных от узловых напряжений (передаточных сопротивлений) по поверхностным проводимостям конечных элементов, что легче всего [18—21] осуществить, исходя из решения уравнения (1) путем обращения матрицы Y,
U=Y~4 = Z7 (2)
Поэтому для решения прямой задачи в дальнейшем будем использовать метод модификаций для обращения матрицы проводимостей [15, 16]. Идея метода заключается не в обычной процедуре - составление системы уравнений равновесия цепи (1) и ее решение (2) известными методами, а в получении решения (матрицы Z в (2)) непосредственно, используя информацию о структуре анализируемой цепи.
Для решения прямой задачи методом модификаций достаточно иметь обратные матрицы элементарных блоков, каковыми являются обратные матрицы для конечных элементов рис. 1 б, полученные любым способом (в том числе и методом модификаций). Обращенная матрица узловых проводимостей для электрической эквивалентной схемы конечного элемента имеет вид [10]
1И =f-1 = ^
где kt = 4/аа ; а - сторона квадрата; ot - удельная поверхностная проводимость конечного элемента.
При построении фантома задача анализа заключается в последовательном объединении эквивалентных схем соответствующих конечных элементов, с обратными матрицами, взвешенными коэффициентами кь полученными при первой дискретизации.
Для пояснения процедуры объединения подсхем рассмотрим схему рис. 2 а, на которой изображен базовый квадратный конечный элемент (один из узлов которого соединен с общим узлом - «землей») и объединяе-
6 — V 2 4 —V 2
3
4-V2 3 2
3
3
3 —2V2 3
4—V2 3
2 ' 3
4-V2
3
6-V2 3 J
(3)
Вісник Національного технічного університету України "КПІ" Серія — Радіотехніка. Радіоапаратобудування.-2010.-№43
7
Радіотехнічні кола та сигнали
мым с ним «незаземленным» конечным элементом. Для «незаземленного» конечного элемента обратная матрица не существует, однако это неудобство можно обойти, если воспользоваться алгоритмом «выращивания» незаземленной подсхемы [22].
Рис.2
Для объединения двух подсхем А и В, из которых подсхема В не соединена с общим узлом (рис.2 а) достаточно выбрать первый из объединяемых узлов (на рис. 2 а это узел №6) и, «заземлив» его, вычислить обратную матрицу подсхемы В. Составить матрицу Ztq, образованную матрицами ZA и ZB, расположенными в ее диагонали
Г^ц Z12 213
Z21 22 2 2гг
22 2 2гг
0 0 0
0 0 0
.0 0 0
0 0 о-
0 0 0
0 0 0
Z2’2' 22', ^2'5
2,2' 2„ 2,5
2з 2' 2, , 2 -3 -3 J
(4)
В полученной матрице 2^о в столбцы нулей (слева от матрицы ZA ) подставить столбец этой матрицы, номер которого соответствует номеру узла объединяемого с условно заземленным узлом матрицы ZB (в рассматриваемом примере это столбец №3). В строки с нулями матрицы ZT0 подставить строку матрицы ZA, номер которой соответствует номеру объединяемого узла этой матрицы. К каждому элементу матрицы ZB прибавить элемент zu матрицы ZA, индексы і которого совпадают с номером объединяемого узла.
8
Вісник Національного технічного університету України "КПІ" Серія — Радіотехніка. Радіоапаратобудування.-2010.-№43
Радіотехнічні кола та сигнали
В рассматриваемом примере в соответствии с вышеописанным правилом из матрицы Z10 (4) получим
Г*и *12
Z21 *22
*31 *3 2
*31 *3 2
*3 1 *3 2
-*3 1 *32
^13 *23 *3 3 *3 3 *3 3 *33
13 *13 *13
23 *23 *23
33 *3 3 *3 3
*2 '2 ' "Т *3 3 *2Г4 + Z33 *2 5 *33
*42 г + *33 *44 + Z33 *45 + *3 3
*52г Т *зз *54 Т *33 *55 + *33
(5)
По матрице ZY1 получим далее обратную матрицу для схемы рис. 2 а после выращивания в ней бесконечной проводимости между узлами №2 и №2 . Для матрицы ZTl такое «выращивание» выполняется по формуле [15-18]
Zz2=Zzi-f%zo%oz’
(6)
ГДQ^Z22-Zir-Zr2+Zrr,L=[(Z12- Zir) (Z22-Z2r)- (Z-j2- г5Ґ)У
foz = t^21_ *2ri) (Z22 ~ *2гз) "
понирования.
Все значения ZtJ берутся из матрицы ZT1
(_Z2-~ Z2B)]; — знак транс-
Описанный выше алгоритм позволяет путем «выращиваний» получить обратную матрицу одной «линейки» конечных элементов на рис. 2,б.
При вычислении обратной матрицы всего фантома в целом необходимо объединять отдельные линейки (рис. 2,б), для чего следует использовать операции алгоритма приведенного выше. Так, для схемы рис.2,6 обратная матрица с условно «заземленным» узлом Л®8 существует и объединение матриц двух «линеек» осуществляется аналогично объединению узлов №6 и №3 для подсхем рис.2 а. Далее, «выращивая» связи между остальными узлами со штрихами и соответствующими им узлами с номерами без штрихов по формуле (6), получим обратную матрицу фрагмента фантома, состоящего из двух «линеек» и т.д. В случае фантома круглой или овальной формы длина отдельных «линеек» будет меняться в соответствии со структурой, полученной в результате первой дискретизации.
Поскольку при изменении узлов подключения источника тока нумерация внешних узлов фантома меняется (общий узел, к которому подключен источник, перемещается вдоль обвода), интерес представляет переход от обратной матрицы для схемы с «заземленным» узлом а к обратной матрице той же схемы с заземленным узлом в. Такой переход легко осуществить, воспользовавшись алгоритмом [22], который проиллюстрируем на
Вісник Національного технічного університету України "КПІ" Серія — Радіотехніка. Радіоапаратобудування.-2010.-№43
9
Радіотехнічні кола та сигнали
матрице третьего порядка. Пусть найдена исходная обратная матрица Ё(-а^с заземленным узлом а
(7)
На базе этой матрицы построим такую же матрицу, однако с учетом «заземленного» узла и записи элементов матрицы (7), которые имеют смысл входных и передаточных сопротивлений относительно входов a, b, c, d
^(0+0X0+0) ^(ь+0X0+0) ^(c+0X0+0) X"(t£+0X0+0)
^(о+оХо+ь)
Z(b+o}(o+b)
^(c+oXo+b)
^(d+oXo+b)
Z(o +0Х0+С) ^(b + oXo+c) ^(c+oXo+c) Z(d+0X0+c)
^(o+oXo + d) ^(b+oXo+d)
^(.c+o'Ko+d)
^(d+oXo+d)-
(8)
При «заземлении» узла d заменим в матрице (8) все индексы а =0 на d. Получим матрицу Z(-d^ с заземленным узлом d и «оторванным» от «земли» узлом а
’^(o+dXo+d) Z(b+d) (o+d) %{c+d) (o+d) -^(d+d)(0+d)
^(o+d) (b + d) Z(b + d)(b + d) ^(c+d) (b + d) Z(d+d)(b+d)
Z(0 + d)(c+d)
Z(b + d)(c + d) ^(c+d) (c+d) j^(d+d)(c + d)
Z(0+d)(d+d) ' Z(b + d)(d+d) ^(c + dXd+d) Z(d+d)(d+d)-
(9)
Раскрывая скобки в индексах (9), получим
Kd) _
обратную матрицу схемы с «заземленным» узлом d и «оторванным» от «земли» узлом а, записанную через элементы матрицы той же схемы, но с «заземленным» узлом а и «незаземленным» узлом b.
При использовании предлагаемого алгоритма решения прямой задачи методом модификаций (обращения матрицы проводимостей системы уравнений равновесия фантома) операции над нулевыми элементами матрицы проводимости (как и при решении прямой задачи методом Гаусса после приведения исходной матрицы проводимостей к ленточному виду) не выполняются. Однако при этом вычисляется обратная матрица, которая в
10
Вісник Національного технічного університету України "КПІ" Серія — Радіотехніка. Радіоапаратобудування.-2010.-№43
Радіотехнічні кола та сигнали
случае фантома, используемого для моделирования в данной работе, всегда будет полностью заполненной. В отличие от решения системы уравнений методом Г аусса информация, которая содержится в обратной матрице системы уравнений равновесия фантома, позволяет значительно упростить алгоритмы вычисления производных от функций цепи (передаточных сопротивлений) по параметрам (удельным поверхностным проводимостям).
Следует также отметить, что метод модификаций позволяет находить обратные матрицы для фантома в случае изменения (например, уточнения) проводимостей отдельных его областей. Для этого достаточно сформировать обратную матрицу для этой отдельной области (подсхемы), удельная проводимость конечных элементов которой в сумме с удельной проводимостью той же области на исходном фантоме даст требуемое новое значение. Для коррекции достаточно «вырастить» короткие замыкания между соответствующими узлами фантома и корректирующей подсхемы.
Аналогично можно сформировать фантом и его обратную матрицу, построив сначала фантом с одинаковой проводимостью всех конечных элементов, а потом скорректировать указанным способом удельные проводимости отдельных областей.
Выводы
1. Предложен метод прямого формирования обратной матрицы фантома обобщенным методом модификаций, для реализации которого не требуется построения матрицы проводимости, а требуется лишь информация о структуре фантома (способе соединения конечных элементов) и удельных проводимостях отдельных конечных элементов.
2. Метод позволяет корректировать значения элементов обратной матрицы фантома при коррекции проводимостей отдельных его областей без нового построения обратной матрицы.
3. Результаты решения прямой задачи содержат избыточную информацию (для вычисления проекций — узловых напряжений — необходимы лишь п < 64 элемента обратной матрицы из тысяч - миллионов при матрицах Y порядковN= 32... 1000.
4. Избыточная для решения прямой задачи информация в виде неиспользуемых значений обратной матрицы является необходимой при вычислении производных от передаточных сопротивлений удельным проводимостям конечных элементов.
5. В работе рассмотрен пример плоского фантома. Поскольку токи между электродами текут по всему объему, при измерениях на физических объектах на результаты вычислений будут содержать дополнительные проекции на плоскость фантома.
6. При рассмотрении объемного фантома математический аппарат решения прямой задачи останется неизменным: увеличится только порядок
Вісник Національного технічного університету України "КПІ" Серія — Радіотехніка. Радіоапаратобудування.-2010.-№43
11
Радіотехнічні кола та сигнали
обратной матрицы и модель конечного элемента.
Литература
1. Cormack A.M. Early two-dimensional reconstruction and recent topics stemming from it // Nobel Lectures in Physiology or Medicine 1971—1980. World Scientific Publishing Co.—1992.—P/551—563.
2. Hounsfield G.M. Computed Medical Imaging // Nobel Lectures in Physiology or Medicine 1971—1980. World Scientific Publishing Co.—1992.—P.568—586.
3. Lauterbur P.C. All science is interdisciplinary — from magnetic moments to umole-cules to men // Lex Prix Nobel. The Nobel Prizes 2003. Nobel Foundation.—2004.—P.245— 251.
4. Mansfield P. Snap-short MRI // Lex Prix Nobel. The Nobel Prizes 2003/ Nobel Foundation.—2004.—P.245—251.
5. Brown B.H., Barber D.C. Electrical Impedance Tomography // Clinical Physics and Physiological Measurement.—1992.—v.13.—Sappl. A, 207 p.
6. Физика визуализации изображений в медицине. Под ред. С.Уэбба.—М.: Мир, 1991, т. 1, 2.—408 с.
7. Электроимпедансная томография / Я.С. Пеккер, К.С. Бразовский, В.Ю. Усов, М.П. Плотников, О.С. Уманский.— Томск: ООО «Издательство научно-технической литературы»,2004.—190с.
8. Murray T., Kagawa Y. Electrical Impedance Computed Tomography Based on a Finite Elements Model // IEEE Trans. On Biomed. Eng.—1985.—v.32.—P.177—184.
9. Сильвестр П., Феррари З. Метод конечных элементов для радиоинженеров и инженеров-электриков.— М.: Мир, 1986.—229с.
10. Рибіна І.О., Гайдаєнко Є.В. Моделювання кінцевого елемента в імпеданс ній томографії // Вісник НТУУ «КПІ». Сер. Радіотехніка. Радіоапаратобудування.—2010.— №4.—С.19—24.
11. Марчук Г.И. Методы вычислительной математики.— М.: Наука, 1980.—535с.
12. Бахвалов Н.С. Численные методы.— М.: Наука, 1973.—631 с.
13. Глориозов Е.А., Ссорин В.Г., Сыпчук П.П. Введение в автоматизацию схемотехнического проектирования.— М.: Сов. Радио, 1976.—224 с.
14. Дмитришин Р.В., Шаповалов Ю.И. Диакоптический алгоритм анализа сложных линейных цепей на ЭВМ // Автоматизация проектирования в электронике.— Киев: Техніка, 1975.—Вып.12.—С.42—46.
15. Рыбин А.И. Решение задач моделирования обращением матрицы методов взаимных производных // Радиоэлектроника.—1978.—№6.—С.35—47. (Изв. высш. учеб. заведений).
16. Рыбин А.И. Численно-символьный метод анализа электрических цепей обобщенным методом модификаций. // Праці Інституту електродинаміки НАН України: Сб. наукових праць.—2002.—№1(1).—С.26—30.
17. Основи теорії кіл: Підручник для ВНЗ. Ч.2 /Ю.О. Коваль, Л.В. Гринченко, І.О. Милютченко, О.І. Рибін / За заг. ред. В.М. Шокала та В.І. Правди.—Харків: ХНУРЕ: Колегіум, 2006.—668 с.
18. Трохименко ЯК., Каширский И.С., Рыбин А.И. Статистический анализ линейных электронных цепей постоянного тока // Радиоэлектроника.-1974.-№6.-С.69-73. (Изв. высш. учеб. заведений).
19. Трохименко Я.К., Каширский И.С., Рыбин А.И. Вероятностный аналіз линейных электронных цепей переменного тока // Радиоэлектроника.—1975.—№6.—С.35—40. (Изв. высш. учеб. заведений).
12
Вісник Національного технічного університету України "КПІ" Серія — Радіотехніка. Радіоапаратобудування.-2010.-№43
