动态感应电流电阻抗断层成像的算法仿真
关键词: 感应电流;电阻抗;计算机模拟;算法;体层摄影术;牛顿迭代法
摘 要:目的 在仿真模型的基础上，研究一种用感应电流激励的动态电阻抗断层成像算法—牛顿迭代法的特性及其对独立测量数的依赖性. 方法 利用有限元法对成像区域进行离散，然后再用牛顿迭代法针对不同的线圈数进行求解、成像，以资比较. 结果 对同一目标，分别就不同的线圈数得到了成像结果，表明牛顿迭代法对电导率扰动的定位是基本准确的，成像误差随线圈数的增加而减小. 结论 用牛顿迭代法解动态感应电流电阻抗断层成像的逆问题是可行的，但在独立测量数小于剖分单元数的情况下，迭代过程不会准确的收敛于实际的电导率分布，而是一种近似;在独立测量数大于剖分单元数的情况下，迭代过程可以收敛于实际的电导率分布，从而得到高质量的重构图像.
Algorithm of dynamic induced current electrical impedance tomography based on simulation phan-tom
MA Jian-Ying，DONG Xiu-Zhen，QIN Ming-Xin，LIU Rui-Gang，YOU Fu-Sheng，XIANG Hai-Yan Department of Medical Electronic Engineering，Fac-ulty of Biomedical Engineering，Fourth Military Medical University，Xi'an710033，China
Keywords:induced current;electric impedance;computer simulation;algorithms;tomography;Newton it-eration method
Abstract:AIM To study the property of the algorithm of dynamic induced current electrical impedance tomography-Newton iteration method and its dependency on independent measurement.METHODS FEM method was used to dis-cretize the imaging area，and then Newton iteration method was used to solve the problem and reconstruct images with different numbers of coils for comparision.RESULTS With different coils，different images were obtained.The results indicated that Newton iteration method could be used to lo-cate the conductivity perturbation accurately，and that when the number of independent measurement increased，the imag-ing error would decrease.CONCLUSION Newton iteration method is a feasible method for the inverse problem of dy-namic induced current electrical impedance tomography.But，when the number of independent measurement is less than the number of FEM element，the iterative process will not converge accurately at the practical conductivity and the re-sult is just an approximate.When the number of independent measurements is larger than the number of FEM element，the iterative process can converge accurately at the practical conductivity，by which high quality images will be obtained.
0 引言
感应电流电阻抗断层成像(induced current elec-tric impedance tomography，ICEIT)是电阻抗断层成像(EIT)技术的一个新的分支，该技术利用载流线圈在成像区域内产生变化的磁场，由电磁感应原理在成像区域内产生感应电场，继而通过区域周围的电极测量区域边界上的电位，经一定的重构算法来得到成像区域内电阻抗分布的图像.ICEIT按成像目标的不同可分为静态和动态两种，静态ICEIT以成像区域内的电阻抗的绝对分布为成像目标，动态ICEIT则以成像区域内的电阻抗变化量的分布为成像目标.到目前为止，只有ICEIT的动态算法可见诸报道.我们假设成像区域为一圆形区域，测量电极等间隔分布在区域的边界上.激励电流为频率是50kHz的正弦电流，大小是500mA.激励线圈的圆心等间隔分布在以成像区域的圆心为圆心的一个圆周上.成像算法需要两组测量数据之差，一组是电导率无扰动时的边界电压，即背景的边界电压，此时电导率为均匀分布，另一组为电导率发生扰动时的边界电压.独立测量数等于电极数减一再乘以线圈数.对给定的有限元剖分模型，设成像区域内的每个单元内的电导率为常数[1] .
动态感应电流电阻抗成像模型如Fig1所示.
图1 动态感应电流电阻抗成像模型示意 略
1 正向问题
ICEIT的正向问题是指已知成像区域电导率的分布和激励线圈的电流，求成像区域的电位分布.由麦克斯韦方程组出发，可得ICEIT问题的微分方程[1] :·[(σ+jωε)ψ]=-jωA·(σ+jωε)2 A=μ0 (σ+jωε)(-ψ-jωA) (1)其中，ψ为成像区域的标量电位，σ为成像区域的电导率，ω是激励电流的正弦角频率，ε为真空介电常数，A是成像区域的矢量磁位，μ0 是真空磁导率，为矢量微分算子，j为虚数单位.边界条件是[1] :ψn=jωA·n (2)由于电导率发生扰动时的矢量磁位和没有扰动时的矢量磁位相比，差别很小，因此，可假定二者相等，这样可使电位的计算得以简化.又由于在50kHz电流的激励下，成像区域内的位移电流和传导电流相比可忽略不记，可假设位移电流为零.由这两个假设，并将电位的实虚部分离，可将微分方程组(1)简化为[1] :·(σψ)=-ωA·σψn=-ωA·n (3)和微分方程(3)相对应的能量泛函是[2] :
f(ψ)=12 ∫Ω [σ ψ 2 ]dΩ+ ∫Ω σωA·ψdS=min (4)Fig2是用于有限元法的剖分模型，剖分规模为4层和6层，节点数分别为81和169，单元数分别为128和288.用有限元法将成像区域离散，计算(4)的左侧并对节点电压矢量求导，可得关于节点电压的线 性方程组[1] :s(σ)v=b(σ) (5)其中，v为节点电压矢量，σ是单元电导率组成的矢量，S(σ)为系数矩阵，b(σ)是常数矢量，他们均是σ的函数.解此方程组，即可得节点电压值v.
图2 用于有限元法的剖分模型 略
2 牛顿迭代法
ICEIT的逆问题是指已知激励线圈的电流和成像区域的边界电压，求成像区域电导率的分布.由于边界电压的测量值中，既有标量电位的成分，又有感应电动势的成分[1] ，由于电导率发生扰动前后，感应电动势变化很小，经过两组测量值相减，可抵消感应电动势的成分，这样，仅能得到标量电位的差值，所以，仅能用于动态算法[1，2] .由方程(5)可得:v=S(σ)-1 b(σ).边界电压差表示为:g=c(s(σ)-1 b(σ)-s(σ0 )-1b(σ0 )) (6)其中g为相邻电极的标量电位差矢量，c为关联矩阵，σ0 为已知的背景电导率.对于n个激励线圈，可得n个形如(6)的方程组，将这些方程组并到一块，即可得逆问题的方程组，可表示为:G=S(σ)-S(σ0 ) (7)这是关于电导率σ的非线性方程组.其中未知数的个数就等于剖分单元的个数.把(7)改写为:F(σ)=S(σ)-S(σ0 )-G=0，用牛顿迭代法解非线性方程组的步骤如下:(1)令i=0，设σi =σ
0 .(2)在σi 处F(σ)=0将用它的切平面F(σ)-F(σi )=F(σ)σ(σ-σi )代替，解之可得:σ-σ
i =△σ i =pinv(F(σ)σ)·(F(σ)-F(σi ))其中，pinv(X)是求广义逆的函数[3] .(3)令σi+1 =σi +△σ i ，i=i+1，重复步骤(2)和(3).
在一定的判据下，迭代若干次数后，即可得方程(7)的在一定误差范围内的近似解.
3 成像结果
Fig3是针对不同的线圈数得到的成像结果.其中，激励线圈的半径为25cm，成像区域的半径为15cm，线圈中心距成像区域的圆心9cm，激励电流为500mA，剖分规模为4层，节点数81，单元数128，电极数为16，背景电导率为100Ωm，扰动也为100Ωm.成像误差定义为:error=1NΣN i=1 (σi -σi' )2其中，N为剖分单元数，σ
i 为第i个单元的实际电导率，σi' 为第i个单元的计算电导率.在计算的过程中，系数矩阵的条件数有随线圈数的增加而增加的趋势，经过适当地选取校正因子，可将系数矩阵的条件数控制在108 数量级以内[4，5] ，这时对成像结果的影响基本可忽略.
Fig4是另一组成像结果，其中剖分规模为6层，电极数16，节点数为169，单元数位288，其他设置与Fig1相同.
图3 - 图4 略
以上试验结果表明，牛顿迭代法对扰动区域的定位是准确的，但当独立测量数小于剖分单元个数的个数时，所得扰动区域的面积明显大于扰动区域的实际面积，随着线圈的增加，成像误差将逐渐减小，当独立测量数大于剖分单元个数的个数时，所得扰动区域的面积明显的减小，逼近实际的面积.成像结果的精度和线圈数并不是直线关系.
4 讨论
用牛顿迭代法解ICEIT的逆问题是可行的，当独立测量数大于剖分单元个数时，迭代过程可收敛于电导率的真实分布，这时将得到高质量的图像，但是，由于受实际系统条件的限制，如测量系统的测量精度等，电极数和线圈数不能无限的增加，到一定的程度，将因为系统误差的增大而使成像质量恶化.这使有限元的剖分规模受到一定的限制.另外，在某些情况下，即使独立测量数很大，牛顿迭代法也可能不收敛，或者收敛速度很慢[6] ，如何合理地处理这些问题，还需作进一步的工作.
参考文献:
[1]Gencer NG，Kuzuglu M，Ider YZ.Electrical impedance tomog-raphy using induced current [J].IEEE Trans Med Imaging，1994;13(2):338-350.
[2]Ruan WX，Guatdo R，Adler A.Experimental evaluation of iter-ative reconstruction methods for induced current electrical impedance tomography [J].IEEE Trans Med Imaging，1996;15(2):180-187.
[3]Tang M，Dong X，Qin M，Fu F，Shi X，You F.Electrical impedance tomography reconstruction algorithm based on gener-al inversion theory and finite element method [J].Med Bio Eng Comput，1998;36(4):395-398.
[4]Liu RG，Dong XZ，Qin MX，Tang MX，You FS，Shi XT.A regularized method in reconstructed algorithm of electrical impedance tomography [J].Di-si Junyi Daxue Xuebao(J Fourth Mil Med Univ)，2000;21(1):107-109.
[5]Hua P，Woo EJ，Webster JG，Tompkins WG.Iterative recon-struction methods using regularization optional current patterns in electrical impedance tomography [J].IEEE Trans Med Imag，1991;10(4):621-628.
[6]Liu RG，Dong XZ，Qin MX，You FS，Shi XT，Fu F，Tang MX.Primary images of static electrical impedance tomography based on physical phantom [J].Di-si Junyi Daxue Xuebao(J Fourth Mil Med Univ)，2001;22(4):297-300.
基金项目: 中国-以色列国际合作项目(99M-0021426)
通讯作者: 董秀珍.Tel.(029)3374848Ext.408 Email.bmee@fmmu.edu
作者简介: 马建英(1970-)，男(汉族)，陕西省蒲城县人.硕士生(导师董秀珍，秦明新).Tel.(029)3374848 Email.bmee@fmmu.edu.cn
第四军医大学生物医学工程系医学电子工程教研室，陕西西安710033
编辑 何扬举, 百拇医药(马建英 董秀珍 秦明新 刘锐岗 尤富生 向海燕)
摘 要:目的 在仿真模型的基础上，研究一种用感应电流激励的动态电阻抗断层成像算法—牛顿迭代法的特性及其对独立测量数的依赖性. 方法 利用有限元法对成像区域进行离散，然后再用牛顿迭代法针对不同的线圈数进行求解、成像，以资比较. 结果 对同一目标，分别就不同的线圈数得到了成像结果，表明牛顿迭代法对电导率扰动的定位是基本准确的，成像误差随线圈数的增加而减小. 结论 用牛顿迭代法解动态感应电流电阻抗断层成像的逆问题是可行的，但在独立测量数小于剖分单元数的情况下，迭代过程不会准确的收敛于实际的电导率分布，而是一种近似;在独立测量数大于剖分单元数的情况下，迭代过程可以收敛于实际的电导率分布，从而得到高质量的重构图像.
Algorithm of dynamic induced current electrical impedance tomography based on simulation phan-tom
MA Jian-Ying，DONG Xiu-Zhen，QIN Ming-Xin，LIU Rui-Gang，YOU Fu-Sheng，XIANG Hai-Yan Department of Medical Electronic Engineering，Fac-ulty of Biomedical Engineering，Fourth Military Medical University，Xi'an710033，China
Keywords:induced current;electric impedance;computer simulation;algorithms;tomography;Newton it-eration method
Abstract:AIM To study the property of the algorithm of dynamic induced current electrical impedance tomography-Newton iteration method and its dependency on independent measurement.METHODS FEM method was used to dis-cretize the imaging area，and then Newton iteration method was used to solve the problem and reconstruct images with different numbers of coils for comparision.RESULTS With different coils，different images were obtained.The results indicated that Newton iteration method could be used to lo-cate the conductivity perturbation accurately，and that when the number of independent measurement increased，the imag-ing error would decrease.CONCLUSION Newton iteration method is a feasible method for the inverse problem of dy-namic induced current electrical impedance tomography.But，when the number of independent measurement is less than the number of FEM element，the iterative process will not converge accurately at the practical conductivity and the re-sult is just an approximate.When the number of independent measurements is larger than the number of FEM element，the iterative process can converge accurately at the practical conductivity，by which high quality images will be obtained.
0 引言
感应电流电阻抗断层成像(induced current elec-tric impedance tomography，ICEIT)是电阻抗断层成像(EIT)技术的一个新的分支，该技术利用载流线圈在成像区域内产生变化的磁场，由电磁感应原理在成像区域内产生感应电场，继而通过区域周围的电极测量区域边界上的电位，经一定的重构算法来得到成像区域内电阻抗分布的图像.ICEIT按成像目标的不同可分为静态和动态两种，静态ICEIT以成像区域内的电阻抗的绝对分布为成像目标，动态ICEIT则以成像区域内的电阻抗变化量的分布为成像目标.到目前为止，只有ICEIT的动态算法可见诸报道.我们假设成像区域为一圆形区域，测量电极等间隔分布在区域的边界上.激励电流为频率是50kHz的正弦电流，大小是500mA.激励线圈的圆心等间隔分布在以成像区域的圆心为圆心的一个圆周上.成像算法需要两组测量数据之差，一组是电导率无扰动时的边界电压，即背景的边界电压，此时电导率为均匀分布，另一组为电导率发生扰动时的边界电压.独立测量数等于电极数减一再乘以线圈数.对给定的有限元剖分模型，设成像区域内的每个单元内的电导率为常数[1] .
动态感应电流电阻抗成像模型如Fig1所示.
图1 动态感应电流电阻抗成像模型示意 略
1 正向问题
ICEIT的正向问题是指已知成像区域电导率的分布和激励线圈的电流，求成像区域的电位分布.由麦克斯韦方程组出发，可得ICEIT问题的微分方程[1] :·[(σ+jωε)ψ]=-jωA·(σ+jωε)2 A=μ0 (σ+jωε)(-ψ-jωA) (1)其中，ψ为成像区域的标量电位，σ为成像区域的电导率，ω是激励电流的正弦角频率，ε为真空介电常数，A是成像区域的矢量磁位，μ0 是真空磁导率，为矢量微分算子，j为虚数单位.边界条件是[1] :ψn=jωA·n (2)由于电导率发生扰动时的矢量磁位和没有扰动时的矢量磁位相比，差别很小，因此，可假定二者相等，这样可使电位的计算得以简化.又由于在50kHz电流的激励下，成像区域内的位移电流和传导电流相比可忽略不记，可假设位移电流为零.由这两个假设，并将电位的实虚部分离，可将微分方程组(1)简化为[1] :·(σψ)=-ωA·σψn=-ωA·n (3)和微分方程(3)相对应的能量泛函是[2] :
f(ψ)=12 ∫Ω [σ ψ 2 ]dΩ+ ∫Ω σωA·ψdS=min (4)Fig2是用于有限元法的剖分模型，剖分规模为4层和6层，节点数分别为81和169，单元数分别为128和288.用有限元法将成像区域离散，计算(4)的左侧并对节点电压矢量求导，可得关于节点电压的线 性方程组[1] :s(σ)v=b(σ) (5)其中，v为节点电压矢量，σ是单元电导率组成的矢量，S(σ)为系数矩阵，b(σ)是常数矢量，他们均是σ的函数.解此方程组，即可得节点电压值v.
图2 用于有限元法的剖分模型 略
2 牛顿迭代法
ICEIT的逆问题是指已知激励线圈的电流和成像区域的边界电压，求成像区域电导率的分布.由于边界电压的测量值中，既有标量电位的成分，又有感应电动势的成分[1] ，由于电导率发生扰动前后，感应电动势变化很小，经过两组测量值相减，可抵消感应电动势的成分，这样，仅能得到标量电位的差值，所以，仅能用于动态算法[1，2] .由方程(5)可得:v=S(σ)-1 b(σ).边界电压差表示为:g=c(s(σ)-1 b(σ)-s(σ0 )-1b(σ0 )) (6)其中g为相邻电极的标量电位差矢量，c为关联矩阵，σ0 为已知的背景电导率.对于n个激励线圈，可得n个形如(6)的方程组，将这些方程组并到一块，即可得逆问题的方程组，可表示为:G=S(σ)-S(σ0 ) (7)这是关于电导率σ的非线性方程组.其中未知数的个数就等于剖分单元的个数.把(7)改写为:F(σ)=S(σ)-S(σ0 )-G=0，用牛顿迭代法解非线性方程组的步骤如下:(1)令i=0，设σi =σ
0 .(2)在σi 处F(σ)=0将用它的切平面F(σ)-F(σi )=F(σ)σ(σ-σi )代替，解之可得:σ-σ
i =△σ i =pinv(F(σ)σ)·(F(σ)-F(σi ))其中，pinv(X)是求广义逆的函数[3] .(3)令σi+1 =σi +△σ i ，i=i+1，重复步骤(2)和(3).
在一定的判据下，迭代若干次数后，即可得方程(7)的在一定误差范围内的近似解.
3 成像结果
Fig3是针对不同的线圈数得到的成像结果.其中，激励线圈的半径为25cm，成像区域的半径为15cm，线圈中心距成像区域的圆心9cm，激励电流为500mA，剖分规模为4层，节点数81，单元数128，电极数为16，背景电导率为100Ωm，扰动也为100Ωm.成像误差定义为:error=1NΣN i=1 (σi -σi' )2其中，N为剖分单元数，σ
i 为第i个单元的实际电导率，σi' 为第i个单元的计算电导率.在计算的过程中，系数矩阵的条件数有随线圈数的增加而增加的趋势，经过适当地选取校正因子，可将系数矩阵的条件数控制在108 数量级以内[4，5] ，这时对成像结果的影响基本可忽略.
Fig4是另一组成像结果，其中剖分规模为6层，电极数16，节点数为169，单元数位288，其他设置与Fig1相同.
图3 - 图4 略
以上试验结果表明，牛顿迭代法对扰动区域的定位是准确的，但当独立测量数小于剖分单元个数的个数时，所得扰动区域的面积明显大于扰动区域的实际面积，随着线圈的增加，成像误差将逐渐减小，当独立测量数大于剖分单元个数的个数时，所得扰动区域的面积明显的减小，逼近实际的面积.成像结果的精度和线圈数并不是直线关系.
4 讨论
用牛顿迭代法解ICEIT的逆问题是可行的，当独立测量数大于剖分单元个数时，迭代过程可收敛于电导率的真实分布，这时将得到高质量的图像，但是，由于受实际系统条件的限制，如测量系统的测量精度等，电极数和线圈数不能无限的增加，到一定的程度，将因为系统误差的增大而使成像质量恶化.这使有限元的剖分规模受到一定的限制.另外，在某些情况下，即使独立测量数很大，牛顿迭代法也可能不收敛，或者收敛速度很慢[6] ，如何合理地处理这些问题，还需作进一步的工作.
参考文献:
[1]Gencer NG，Kuzuglu M，Ider YZ.Electrical impedance tomog-raphy using induced current [J].IEEE Trans Med Imaging，1994;13(2):338-350.
[2]Ruan WX，Guatdo R，Adler A.Experimental evaluation of iter-ative reconstruction methods for induced current electrical impedance tomography [J].IEEE Trans Med Imaging，1996;15(2):180-187.
[3]Tang M，Dong X，Qin M，Fu F，Shi X，You F.Electrical impedance tomography reconstruction algorithm based on gener-al inversion theory and finite element method [J].Med Bio Eng Comput，1998;36(4):395-398.
[4]Liu RG，Dong XZ，Qin MX，Tang MX，You FS，Shi XT.A regularized method in reconstructed algorithm of electrical impedance tomography [J].Di-si Junyi Daxue Xuebao(J Fourth Mil Med Univ)，2000;21(1):107-109.
[5]Hua P，Woo EJ，Webster JG，Tompkins WG.Iterative recon-struction methods using regularization optional current patterns in electrical impedance tomography [J].IEEE Trans Med Imag，1991;10(4):621-628.
[6]Liu RG，Dong XZ，Qin MX，You FS，Shi XT，Fu F，Tang MX.Primary images of static electrical impedance tomography based on physical phantom [J].Di-si Junyi Daxue Xuebao(J Fourth Mil Med Univ)，2001;22(4):297-300.
基金项目: 中国-以色列国际合作项目(99M-0021426)
通讯作者: 董秀珍.Tel.(029)3374848Ext.408 Email.bmee@fmmu.edu
作者简介: 马建英(1970-)，男(汉族)，陕西省蒲城县人.硕士生(导师董秀珍，秦明新).Tel.(029)3374848 Email.bmee@fmmu.edu.cn
第四军医大学生物医学工程系医学电子工程教研室，陕西西安710033
编辑 何扬举, 百拇医药(马建英 董秀珍 秦明新 刘锐岗 尤富生 向海燕)