CC BY
33
Радіоелектроніка біомедичних технологій
РАДІОЕЛЕКТРОНІКА БІОМЕДИЧНИХ ТЕХНОЛОГІЙ
UDC 621.372.061
FEATURES OF SOLVING THE ELECTRICAL IMPEDANCE TOMOGRAPHY INVERSE PROBLEM BY ZONES CONDUCTIVITIES
METHOD
I. Sushko, postgraduate student,
O. Rybin, Doctor of Science (Technics), professor,
National Technical University of Ukraine
“Kyiv Polytechnic Institute ”, Kyiv, Ukraine
ОСОБЛИВОСТІ РОЗВ’ЯЗАННЯ ЗВОРОТНОЇ ЗАДАЧІ ІМПЕДАНСНОЇ ТОМОГРАФІЇ МЕТОДОМ ЗОН ПРОВІДНОСТІ
Сушко І.О.,аспірантка, Рибін О.І,. д.т.н., проф.,
Національний технічний університет України "Київський політехнічний інститут", м. Київ, Україна
Introduction
The solution of Electical Impedance Tomography (EIT) inverse problem allows to operate with derivatives of transfer resistances on surface zones conductivities by zones conductivities method [1, 2]. These zones consist of great number of square simple [3] finite elements; values of derivatives on zones conductivities are commensurate with values of transfer resistances - that is not true for derivatives on square simple finite elements (FE) conductivities [5, 6]. The phantom model consists of small size FE (number of elements is 800 - 4000 usually). It requires to make and to inverse the derivatives matrix with order of 800 - 4000 at each iteration. We can conclude that derivatives matrix is always ill-conditioned considering very low sensitivity of transfer resistance to conductivity change of one small finite element [5].
The classical algorithm for solving the inverse problem
The classical algorithm for solving the EIT inverse problem [1] is as follows:
1. To divide (digitize) the study object section to FE (pixels) replacing finite elements with their equivalent schemes. The number of FE М is determined the number of independent transfer resistances (section outline voltages at different source connections) necessary measurements.
2. To carry out М measurements of independent transfer resistances at different current source connections.
106
Вісник Національного технічного університету України "КПІ" Серія — Радіотехніка. Радіоапаратобудування.-2012.-№51
Радіоелектроніка біомедичних технологій
3. To solve the direct problem with known approximate distribution of some FE surface conductivities ck . To take values of all surface conductivities equal without priori information about image.
The system of phantom equilibrium equations should be formulated for solving the EIT inverse problem and solve it concerning the voltages on measuring electrodes. The combinations of inverse conductive matrix elements give required transfer resistances.
YxU = I; Y 1 xI = ZxI = U
4. To compare measured and calculated values of transfer resistances and to find errors AZ j or AU .
5. To estimate numerical value of discrepancy 5 using chosen discrepancy norm (Chebyshev, Hemming or Euclid).
6. The computations must be stopped and calculated values Ck should be assigned to corresponding FE if 5 < 8гр. The corresponding color or “gray” level
should be assigned to each conductive element ck and displayed as image element. Here 5гр — chosen boundary value of discrepancy, under which the image is considered satisfactory.
If 5 > 5гр then go to 7 item of the algorithm.
7. To calculate derivatives dZj / dak on all current surface conductivities
for all M transfer resistances Z j using current inverse matrix Z .
8. To formulate the system of equations
-(dZj /d<7k)x A<jk = AZj (la)
or ak=Aak+ak (lb)
where Ack — current values corrections of surface conductivities decreasing the discrepancy norm.
9. To solve the equation (1) and to find unknown increments of conductivities
-[dUj /dak] 1 xAU = Ac .
10. To correct the values of all computing surface conductivities.
<Уь
A<jk + ak
where Gk — precise current value of conductivity.
11. To go to item 3 of the algorithm.
(2)
Вісник Національного технічного університету України "КПІ" Серія — Радіотехніка. Радіоапаратобудування.-2012.-№51
107
Радіоелектроніка біомедичних технологій
There are array of “narrow” places in this algorithm, which require discussions about next modernization.
Since the structure of phantom electrical scheme is rarefied as you know, the solving methods of linear equation systems with rarefied coefficient matrix [9]
can be used successfully for the calculating of the current column of voltages U
although the order of matrix Z N-1000... 4000 with М-1000... 4000.
However, the derivatives computations by the known formulas require sufficiently large costs of machine time for all M2 = (1000)2...(4000)2 elements in expression (1). But the biggest problem is the multiple solution of correction equation (2) (iterative procedure).
The matrix [5U j / dak ] in (2) has the order M, but is not related with any structure and is not rarefied. The derivatives matrix inversion with high order gives the large error by reason of ill-conditionality and its big size. The reduction of derivatives matrix size is possibly using zones conductivities method [1, 2, 5]. However, there are situations in which the derivatives matrix from transfer resistances on surface zones conductivities is ill-conditioned using this method. The order of such matrix according to zones configuration is often equal 6.30. The derivatives matrix ill-conditionality for zones is determined by low sensitivity of the specific values of transfer resistances to large changes of the phantom inhomogeneity conductivity, especially if such conductivity is greater than “background” conductivity.
Consider the example to illustrate the emergence of ill-conditioned derivatives matrix. The homogenous phantom (pic.1) consists of FE with surface conductivity of each element <r = 1.
(2) (1)F
(3)
(4)
(5)
(6)
(0)
(15)
(3).
Ш
(4)
(7)
(8)
(1)!° (0)
(9) (15)
(14)
(13.
(12)
(11)
(10)
(14.
it
4
(13.
(12,
(5)
47) (8)
(9)
W10)
(11)
І
Pic.1. Phantom with 14 conductive zones
Pic.2. Phantom with 6 conductive zones
108
Вісник Національного технічного університету України "КПІ" Серія — Радіотехніка. Радіоапаратобудування.-2012.-№51
Радіоелектроніка біомедичних технологій
(3)
(4)
(1)
(5)
'Ш
(оЩ\
(15^^^ШШэ) (15Щ
(14^^^Ш?Т10) (14)
(Із) ■■ (11) (
(21
(3)
(4)
(7) (1)\
(13)
п
m
V-
(12.
(5)
(6)
(7)
(8)
(9)
(11)
(10)
Pic.3. Phantom with inhomogeneity Pic.4. Phantom with 7 conductive zones
The electrode numeration is given in brackets. Herewith the initial phantom can be divided to 14 zones [1, 2]. Zones I.. .IV (pic. 1) are combined in one (zone 1 - pic.2), further there are zones II - V, zones X...XIV are combined in zone VI for simplicity. The current source (/=1) is connected between nodes 0.8.
The derivatives matrix of phantom pic.2 has the form (3)
0,032186 -0,10484 -0,13998 -0,1320 -0,067327 -0,013301
-0,042626 -0,12226 -0,14889 -0,16738 -0,10731 -0,039082
-0,106697 -0,12933 -0,1610 -0,21398 -0,12325 -0,060073
-0,060073 -0,12325 -0,21398 -0,1610 -0,12933 -0,106697
-0,039082 -0,10731 -0,16738 -0,14889 -0,12226 -0,042626
-0,013301 -0,067327 -0,1320 -0,13998 -0,10484 0,032186
Consider the simplest case as calculating example with surface conductive of
“measured” phantom a = 0,5. Then the differences AU between the “measured” voltages (transfer resistances) in phantom nodes 1, 4, 7, 9, 12, 15 and computing values in the same nodes (with surface conductivities a = 1) are Au = 0,42526, Au4 = 0,62754, Au7 = 0,79461, Aug = 0,79461, Au12 = 0,62754,
Au15 = 0,42526 respectively.
It is necessary to inverse the derivatives matrix (3) for calculating the column of surface conductivities increments for zones I.VI A a.
The inverse matrix (dZj / dak) has the following form:
(dZ. / da )—1 = [aaa-a-a- c]T, where C — row with next values c = \ab c — a — b — c], a = -2,2641 • 1013, b = 6,7819 -1013, c = —2,7084 -1013; T— transponse sign.
To make sure that the obtained matrix is really inverse to derivatives matrix (2) (although with significant errors in this example) multiply both matrix.
Вісник Національного технічного університету України "КПІ" Серія — Радіотехніка. Радіоапаратобудування.-2012.-№51
109
Радіоелектроніка біомедичних технологій
The main diagonal of the product (instead of ones) has values:
Diag = [1,0005; 0,99951; 1,001; 1,0005; 1,0005; 1,0005],
the elements of resulting matrix with numbers 1,3; 1,4; 1,6 are equal to zero, every other elements are non-zero (instead of zero) with values within the limits of 6,104-10-5 — 1,9531-10 3.
So the inverse matrix computational error is not much less than 2 -10-3. The results of low-order matrix inversion don’t allow to calculate zones conductive increments with given deviations between “measured” voltages and computing
values. The multiplication of derivatives matrix and column Aak with all values equal to -1 is exactly equal to AU. So the solution is Aa = [-1 -1 -1 -1 -1 - 1f, where — transponse sign.
It can be proved by number substituting. The solution (1b) in the form of (2) considering the presented values of inverse matrix [dU; / dak ] 1 gives the column Aak = [000000]T, which in fact does not provide even approximate equality (1b).
Herewith matrix determinant [dZnep / dat ] is non-zero, but it is very “small”
value [6], resulting in “big” values of inverse matrix elements [dZnep / dat]-1. So
the equation (2) is “machine infinity multiplied by machine zero” indeterminacy. The disclosure of this indeterminacy is possibly solving the equation (2) by methods for solving the ill-conditioned systems of equations [7, 8]. This solution is not complicated based on the fact that order of system is small (6...30). Another possible way to increase the calculation accuracy is the choice of current source connection to zones phantom.
Asymmetric source connection
The calculations were carried out for source positions with constant phantom position (pic.1). It was current source connections to the pair of nodes 1 — 9; 2 — 10; 3 — 11; 14 — 6; 15 — 7. The biggest derivatives matrix determinant
is on source connection 4 - 12. Consider the phantom example (pic.1) with source connection to 4 - 12 nodes (maximum asymmetry and voltages sensitivity to zones conductivities changes). The surface conductivity of each zone is a = 1 as before. So the derivatives matrix from voltages on electrodes (from transfer resistances) on zones surface conductivities becomes asymmetric (instead of (2)) and its determinant becomes maximum. The product of derivatives matrix with 14 order and inverse matrix gives ones in the main diagonal and maximum non-zero element out of main diagonal is equal to 2,4869 -10-14. If use homogenous phantom with surface conductivity a = 1 as initial approximation and required phantom (also homogenous) has surface conductivity a = 1/2 than voltages increments on electrodes (number of electrodes is started from 0) take the following values:
110
Вісник Національного технічного університету України "КПІ" Серія — Радіотехніка. Радіоапаратобудування.-2012.-№51
Радіоелектроніка біомедичних технологій
(AU)Т = [0,6352; 0,59699; 0,53557; 0,45665; 0,44431; 0,51905; 0,5825; 0,63552; 0,68853; 0,75199; 0,82673; 0,81439; 0,73547; 0,67405] (4)
The product of column with voltages increments and inverse derivatives matrix (2) gives the column of values in this case
Aa = [-1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1]T, (5)
where T— transponse sign. The received values are accurate and entirely consistent with values had to be obtained (with good conditioned derivatives matrix) for phantom pic.2 with 6 zones. Herewith the maximum difference between results obtained by substitution (5) to (1b) and initial value AU to (4) is equal 8,8818-10-16. The discrepancy is very small (instead of symmetric source connection to the phantom) and the results of derivatives matrix inversion allow to organize the iterative procedure.
Obtain the following zones specific conductivities with current source connection to 4.. .12 nodes and zones phantom position pic.1 imposed to the phantom with inhomogeneity pic.3 (a = 10):
ax = 1,0002 ; a2 = 0,9956; a3 = 0,9978; a4 = 1,0564; a5= 1,0196;
a= 1,2491; a= 1,1381; a= 1,5971; a= 1,8264; a10 = 2,3954;
an= 2,8988; a12 = 3,8421; a13 = 5,0670; a14 = 10,8375.
The summary non-normalized discrepancy calculated as difference between voltages “measured” for phantom pic.3 and computing values for conductive zones in iterative procedure is equal to 0,1212. The completely satisfactory results were obtained also for source connection to another opposite about phantom pic.1.pairs of nodes.
The regularization procedure for derivatives matrix inversion
of zones phantom
The most common way for ill-conditioned systems of equations (1) is regularization method by A. Tyhonov [10] using for solving the EIT problems. The process of iterative problem solving with regularization (on the basis of Newton-Raphson classical method) is as follows.
1. To make the matrix [dU^ / dak ]T;
2. To find the product of matrixes [dUJ /dak]T x[dUJ /dak] = Y , the obtained matrix is symmetric;
3. To make the diagonal matrix D y from the main diagonal of matrix Y.
4. To make equation choosing “small” positive number p<1 from equation
(1b) = = _ _ _
(Y + P x Dy ) x [Aak ] = -[0U7- / da ]T x AU (6)
Вісник Національного технічного університету України "КПІ" Серія — Радіотехніка. Радіоапаратобудування.-2012.-№51
111
Радіоелектроніка біомедичних технологій
that coincides with (1b) for p = 0
5. To solve (6)
[Aak ] = -(V + p x Dv )-1 x \dUj / dak ] x AU (7)
Consider phantom with zones pic.4 with imposed inhomogeneity pic.3 with surface conductivity cr = 10 on the background cr = 1 to illustrate abilities of method. The parameter P = 0,35/3. The computational results (1.. .5 iterations) are given in table 1. The phantom pic.3 is stationary and imposed phantom with zones turns for each source connection.
Table 1. Computing values of zones conductivities with source connection to opposite electrode pairs___________________________________________________
Source connection to electrode pairs
№ zone 0..8 1..9 2..10 3..11 4..12 5..13 6..14 7..15
1 1,2768 1,7557 1 1 1 1 2,4380 8,4564
2 1,5880 1,3305 1 1 1 1 1,2141 2,0282
3 1,1014 1,0311 1 1 1 1 1,002 0,9971
4 0,9662 1,0019 1 1 1 1 1,0027 1,0034
5 1,2528 1,0633 1 1 1 1 1,0677 1,0624
6 7,3275 1,5208 1 1 1 1 1,4937 1,1703
7 39,931 2,9746 1 1 1 1 1,9304 0,6805
discrepancy 0,0591 0,0606 0,1232 0,1029 0,0905 0,1068 0,0824 0,0456
8..0 9..1 10..2 11..3 12..4 13..5 14..6 15..7
1 9,8885 9,0161 5,4137 1 1 1,3525 2,6154 1,6271
2 2,0855 1,8450 1,6158 1 1 2,6135 1,2036 1,2948
3 1,0424 1,0834 1,0169 1 1 1,7990 1,0495 1,0818
4 1,0168 1,0295 1,0045 1 1 1,2837 1,0015 1,0216
5 1,0621 1,0645 0,9297 1 1 1,6579 1,0221 1,0780
6 1,4118 1,6174 1,8088 1 1 2,8508 4,0250 2,0290
7 -0,6185 0,8847 1,7347 1 1 0,4720 3,1252 2,9988
discrepancy 0,03575 0,0394 0,0617 0,4616 0,1352 0,0960 0,0840 0,0612
The results of zones surface co nductivities rec onstruction are conformed to distribution pic.3. This is easily seen considering zones with current source connection in positions 0.8 and 8.0. Some of connections are insensitive to changes in surface conductivities of inhomogeneity phantom, but it does not affect on result of summing the phantoms.
Conclusions
1. Zones conductivities method allows to carry out the image reconstruction in improved conditions (in terms of complexity and accuracy of solving problems).
2. It is necessary to use regularization methods for solving ill-conditioned systems of equations or proposed in this article asymmetric source connection to zones phantom to eliminate ill-conditionality of derivatives matrix.
3. The regularization method by A. Tyhonov allows to solve (as calculating results show) ill-conditioned matrixes like derivatives matrixes from transfer re-
112
Вісник Національного технічного університету України "КПІ" Серія — Радіотехніка. Радіоапаратобудування.-2012.-№51
Радіоелектроніка біомедичних технологій
sistances (nodes voltages) on FE surface conductivities. The additional phantom structuring by zones enables to decrease radically the order of derivatives matrix. It significantly increases the speed of EIT inverse problem solving. It demonstrates the perspective application of zones conductivities method with regularization method by A. Tyhonov.
References
1. Рибіна І. О. Розв’язання зворотної задачі імпедансної томографії методами зон провідностей та зворотної проекції / І. О. Рибіна, О. І. Рибін, О. Б. Шарпан // Вісник НТУУ «КПІ». Сер. Радіотехніка. Радіоапаратобудування.— 2011. — № 45.— С. 5 — 18.
2. Рибіна І. О. Метод променів провідностей та моделювання фантома в імпеданс-ній томографії / І. О. Рибіна // Вісник ЖДТУ. — 2010. — т.8. — 4. — С. 21 — 28.
3. Рибіна І. О. Моделювання кінцевого елемента в імпедансній томографії /
О. І. Рибіна, Є. В. Гайдаєнко // Вісник НТУУ «КПІ». Сер. Радіотехніка. Радіоапаратобудування. — 2010. — № 41. — С. 19 — 24.
4. Рибіна І. О. Обчислення похідних від передаточного опору по поверхневій провідності кінцевих елементів при розв’язанні зворотної задачі імпедансної томографії методом зон провідностей / І. О. Рибіна, О. І. Рибін, О. Б. Шарпан // Вісник НТУУ «КПІ». Сер. Радіотехніка. Радіоапаратобудування. — 2011. — № 44. — С. 5 — 21.
5. Сушко И. А. Сравнение классического метода решения обратной задачи с методом «зон» проводимости / И. А. Сушко, А. И. Рыбин // Вісник НТУУ «КПІ». Сер. Радіотехніка. Радіоапаратобудування. — 2012. — № 48. — С. 19 —24.
6. Сушко І. О. Потенційна чутливість імпедансної томографії / І. О. Сушко,
Є. В. Гайдаєнко, О. А. Якубенко // Вісник НТУУ «КПІ». Сер. Радіотехніка. Радіоапаратобудування. — 2012. — № 49. — С. 19 — 24.
7. Бахвалов Н. С. Численные методы. — М. : Наука, 1973. — 631с.
8. Марчук Г. И. Методы вычислительной математики. — М. : Наука, 1977. — 331с.
9. Сушко І. О. Алгоритм розв’язання прямої задачі імпедансної томографії методом модифікацій / І. О. Сушко // Вісник НТУУ «КПІ» Сер. Радіотехніка. Радіоапаратобудування. — 2011. — № 47. — С.165 — 175.
10. Тихонов А. Н. Методы решения некорректных задач / А. Н. Тихонов,
В. Я. Арсенин. — М. : Наука, 1979.
References
1. Rybina I. O., Rybin O. I., Sharpan O. B. Rozvyazannia zvorotnoi zadachi impedansnoi tomografii metodami zon provodimosti ta zvorotnoi proektsii. Visnyk NTUU “KPI”. Ser. Ra-diotehnika. Radioaparatobuduvannia, 2011, no. 45, pp. 5-18.
2. Rybina I. O. Metod promeniv providnosti ta modeliuvannia fantoma v impedansnii tomografii. Visnyk ZhDTU, 2010, vol. 8, pp.21-28.
3. Rybina I. O., Gaidaienko E. V. Modeliuvannia kintsevogo elementa v impedansnii tomografii. Visnyk NTUU “KPI”. Ser. Radiotehnika. Radioaparatobuduvannia, 2010, no. 41, pp. 19-24.
4. Rybina I. O., Rybin O. I., Sharpan O. B. Obchyslennia pohidnyh vid peredatochnogo oporu po poverhnevii providnosti kintsevyh elementiv pry rozvyazanni zvorotnoi zadachi impedansnoi tomografii metodom zon providnosti. Visnyk NTUU “KPI”. Ser. Radiotehnika. Radioaparatobuduvannia, 2011, no. 44, pp. 5-21.
5. Sushko I. A., Rybin A. I. Sravnenie klassicheskogo metoda resheniia obratnoi zadachi s metodom “zon” provodimosti. Visnyk NTUU “KPI”. Ser. Radiotehnika. Radioaparatobuduvannia, 2012, no. 48, pp.19-24.
Вісник Національного технічного університету України "КПІ" Серія — Радіотехніка. Радіоапаратобудування.-2012.-№51
113
Радіоелектроніка біомедичних технологій
6. Sushko I. O., Gaidaienko E. V., Jakubenko O. A. Potenciina chutlyvist impedansnoi tomografii. Visnyk NTUU “KPI”. Ser. Radiotehnika. Radioaparatobuduvannia, 2012, no. 50, pp. 19-24.
7. Bahvalov N. S. Chislennyie metody. Moscow, Nauka, 1973, 631p.
8. Marchuk G. I. Metody vychislitelnoi matematiki. Moscow, Nauka, 1977, 331p.
9. Sushko I. O. Algoritm rozvyazannia priamoi zadachi impedansnoi tomografii metodom modyfikacii. Visnyk NTUU “KPI”. Ser. Radiotehnika. Radioaparatobuduvannia, 2011, no. 47, pp.165-175.
10. Tihonov A. N., Arsenin V. Ja. Metody resheniia nekorrektnykh zadach. Moscow, Nauka, 1979.
Сушко І.О., Рибін О.І. Особливості розв’язання зворотної задачі імпедансної томографії методом зон провідності. Проведено порівняння сходимості ітераційної процедури розв ’язання зворотної задачі імпедансної томографії класичним методом Ньютона-Рафсона, методом регуляризації за А.Н. Тихоновим та методом асиметричного підключення джерела струму до зонного фантома. Показано погану зумовленість матриці похідних, що приводить до дуже великих похибок. При довільному розташуванні зонного фантома відносно джерела струму регуляризація за Тихоновим дає задовільні результати. Обрані для розв ’язання зворотної задачі методи легко програмуються і зручні на практиці завдяки їх використанню разом із методом зон провідності.
Ключові слова: імпедансна томографія, метод зон провідності, зворотна задача, погано зумовлена матриця похідних, фантом, регуляризація.
Сушко И.А., Рыбин А.И. Особенности решения обратной задачи импедансной томографии методом зон проводимости. Проведено сравнение сходимости итерационной процедуры решения обратной задачи импедансной томографии классическим методом Ньютона - Рафсона, методом регуляризации по А.Н. Тихонову и методом асимметричного подключения источника тока к зонному фантому. Показано плохую обусловленность матрицы производных, что приводит к непозволительно высоким погрешностям. При любом положении зонного фантома относительно источника тока регуляризация по Тихонову даёт удовлетворительный результат. Избранные для решения обратной задачи методы легко программируются и удобны на практике благодаря их использованию совместно с методом зон проводимости.
Ключевые слова: импедансная томография, метод зон проводимости, обратная задача, плохо обусловленная матрица производных, фантом, регуляризация.
Sushko I., Rybin O. Features of solving the Electrical Impedance Tomography inverse problem by zones conductivities method. The comparison of the iterative process convergence for the EIT inverse problem solving by classical Newton-Raphson method, regularization method by A. Tyhonov and asymmetric current source connection to zones phantom is carried out. Ill-conditionality of the derivatives matrix is shown. It leads to very big errors. The regularization method by A. Tyhonov gives satisfactory results at any position of zones phantom relatively current source connection. Chosen methods for the inverse problem are programmed easy and are convenient in practice, because of their usage with zones conductivities method.
Keywords: Electrical Impedance Tomography, zones conductivities method, inverse problem, ill-conditioned derivatives matrix, phantom, regularization.
114
Вісник Національного технічного університету України "КПІ" Серія — Радіотехніка. Радіоапаратобудування.-2012.-№51
