三种电阻抗断层动态成像算法的仿真比较研究

【作 者】帅万钧，董秀珍，吕晓宁，晁勇

1 解放军总医院第一附属医院医学工程科，北京，100048

2 第四军医大学生物医学工程系，西安，710032

3 海军总医院全军航海航空医学专科中心，北京，100048

0 前言

电阻抗断层成像技术（Electrical Impedance Tomography，EIT）是新兴的功能成像技术，它根据生物体内不同组织在不同生理或病理条件下具有不同的电阻抗特性这一原理，通过给目标体内注入安全微弱的交流电信号，测量体表电极间电位差来重建体内的电阻抗分布或其变化的图像。与现有医学影像技术相比，该技术不使用核素或放射线，对人体无电离或辐射损伤，可多次测量、重复使用。而且，其设备结构简单、测量简便、成本低廉，是一种具有诱人应用前景的新型无损伤功能成像技术[1,2]。

EIT按成像目标划分，可分为静态成像和动态成像。静态成像是以电阻抗分布的绝对值为成像目标，而动态成像是以电阻抗分布的相对值（差异）为成像目标。一般而言，静态成像能够获得电阻抗的绝对分布信息，但对硬件采集系统要求苛刻，重建算法复杂，计算量大，目前仅限于实验室研究。动态成像以变化量为重建目标，使用测量数据的差值进行成像，减少了测量的系统误差，使该技术有望由实验室研究进入临床应用。

EIT成像算法一直以来就是研究的热点之一[3,4]。其动态成像算法主要有反投影算法(Back-Projection Algorithm)、牛顿一步重建法(Newton’s One-Step Error Reconstructor,NOSER)和灵敏度矩阵重建算法(Sensitivity Matrix Algorithm)。本文基于仿真研究，对此三种EIT动态成像算法的性能进行了比较，以探讨各自的优劣和适应性。

1 三种EIT动态成像算法

EIT所研究的是一个具有特殊边界条件的电场，其场域的数学描述可用一组偏微分方程来表达，EIT的逆问题求解可分为迭代数值计算和非迭代近似计算两大类[4]。非迭代线性近似方法常用于动态EIT成像，可表达为：

其中，△U为边界电压变化向量，△C为电导率变化向量，S为敏感矩阵，它表示△C与△U之间的线性敏感关系，由式（1）进一步可得：

式(2)表明可通过边界测量电压的变化量计算区域内扰动电导率的分布，此即为动态EIT重建的原理。其中敏感矩阵S是动态的，电极的位置和激励测量模式以及场域的边界外形都直接影响矩阵S求逆的稳定性。

1.1 反投影算法

等位线反投影算法思想源于X射线计算机断层成像（X-CT）的反投影技术[5]。其思路是：在电导率变化较小时，等电位线间的电导率的平均变化与相应等位线间的电位差变化成线性关系，并认为等位线分布在电导率变化前后不变。因而以相邻两条等位线为一个投影区域，将等位线间的电位差变化在此区域内投影，再将不同反投影区域内的反投影结果叠加起来，就可以得到电导率分布变化的结果。

该处理使得EIT重建不必对动态敏感矩阵S进行求逆运算，可通过近似线性处理直接得到反投影矩阵B，从而有△C=B·△U，使得计算量小，重建迅速。反投影矩阵B的计算方法如下：在将场域离散后得到的有限元模型上，假定各个有限单元的电导率均匀，利用有限元法计算EIT正问题，得到场域内所有节点的电位分布；然后根据每次激励下各个电极的电位和各个单元节点的电位大小，将各个单元划分到相应测量电极对的投影区域中。设有限元模型的单元数为M，EIT的电极对为N，则B为(M×N)矩阵，有△Ci=Bij·△Uj,(i=1,2,Λ,M; j=1,2,Λ,N)，Bij表示第j个边界电压变化投影到第i个有限单元时的系数。当该单元电位位于第j对测量电极电位之间时，令Bij=1，否则Bij=0。

1.2 牛顿一步重建法

NOSER算法是由美国Cheney教授于1990年首次提出的[6]。该算法实质是Newton法一步。Cheney认为，如果Newton法中设定的初始电导率足够接近于真实电导率分布，则Newton法中的迭代次数会大大减小，甚至只需要在初始电导率的基础上迭代一次就可以得到场域电导率分布的近似解。而且，Newton法之所以运算缓慢的原因在于每次迭代时都需要重新计算一次雅可比矩阵，如果只计算一步，运算量将大大降低。

NOSER算法确定初始电导率的策略是：在假设被测场域的电导率均匀分布的前提下，利用最小二乘法求得当前边界电压下所对应的场域电导率的分布。然后再计算雅可比矩阵，进行图像重建。

以此为基础，重庆大学罗辞勇博士进一步提出了FNOSER算法[7]，即快速(fast)NOSER算法。由于不需迭代，所以可事先设定初始电导率的分布，计算出雅可比矩阵。在实际检测中，直接读取该矩阵数据进行图像重建计算，从而节省大量运算时间。

1.3 灵敏度矩阵重建算法

灵敏度矩阵重建算法是由Murai和Kagawa于1985年提出来的[8]。其主要思想是基于Geselowitz的阻抗电场理论，利用有限元方法将场域等效为一电阻网络模型，根据这一网络模型建立总体传导阻抗变化与单个单元中电阻抗变化之间的对应关系，从而利用这一对应关系来修正阻抗分布。

Geselowitz的阻抗电场理论为：

其中σ为电导率初始分布，其相应的电位分布为φ(σ)，同时在电导率分布变化后的边界电位为ψφ(σ+△σ)，且△σ为电导率分布变化量。EIT成像的主要问题是通过测量边界电流驱动所导致的边界电压或边界节点间的传导阻抗Z来估计电导率分布σ+△σ，电位梯度 和 是电导率分布的函数。若对于已知的σ，则只需求出△σ即可得到新的电导率分布。进一步近似处理，有：

上式中的Sij表示第i对电极对于第j个单元的电导率灵敏度系数。S被称为灵敏度矩阵，对其求逆，即可得电导率变化值△σ。

2 仿真研究方法

所建立的场域有限元仿真模型如图1所示，在圆模型边界等间隔放置16个电极。为尽量模拟人体参数,模型半径设为14 cm，采用三角单元剖分，共有512个剖分单元。首先设置各个剖分单元电导率均匀分布，大小为0.2 S/m。从边界注入1 mA激励电流，通过数值计算得到扰动前的边界电极上的电压数据，以模拟实际检测中边界上所测量得到的电压数据Uo，然后设置扰动目标，其阻抗扰动量的大小按背景电导率的百分比计算，再通过数值计算获得存在扰动后的边界电极上的电压数据，记为U1。在实际应用中，U0和U1即为两个不同时刻检测的边界电压。对U0和U1按EIT算法进行图像重建，可获得到扰动前后所发生的阻抗分布变化图像。

图1 有限元仿真模型及电极分布Fig.1 Finite element Model and electrodes distribution

另外,为进一步讨论各算法的抗噪性能，在实际测量中用数据噪声进行了模拟。边界测量数据的加噪方法如下：

其中，Un为加噪前边界测量数据，U'n加噪后边界数据，Gn为高斯分布序列，其均值为0，方差为1，Kn为噪声幅度系数，a为噪声强度系数。利用计算机可以快速有效生成Gn，Kn计算方法如下：

其中，为边界数据的独立测量数。

3 结果及讨论

有限元仿真模型的参数如上所述，目标的扰动量设为10%，其位置和尺寸如图2第一行所示。图2和图3中“BP”、“NS”和“SE”分别表示该行图像为“反投影算法”、“NOSER”算法和“灵敏度矩阵算法”的成像。

3.1 图像分辨率的比较

从图2可见，三种算法对于径向不同位置的目标均能进行有效重建(前三列)；对于位置相距较远的多目标，其重建图像均能反映出目标的位置信息(后两列)；但对于相距较近的目标，重建图像仅能显示一个目标，无法分辨出实际目标个数(第四列)。相比较而言，NOSER算法的重建图像分辨率高，目标区域突出，灵敏度矩阵算法次之，但在二者图像的目标区域附近存在明显亮斑伪影。反投影算法的目标区域模糊，图像分辨率较低，但无亮斑伪影。总体上讲，在理想无噪声条件下，NOSER算法和灵敏度矩阵算法的图像分辨率较高；反投影算法的图像平滑，伪影相对较少。

3.2 抗噪性能的比较

图2 仿真模型上的反投影算法、NOSER算法和灵敏度矩阵算法的重建图像Fig.2 Images reconstructed by back-projection algorithm,Newton’s one-step error reconstructor algorithm and sensitivity matrix algorithm based on finite element model

前面的仿真研究中，边界电压数据是通过计算机模拟计算得到的，它们不含有任何的误差和噪声，是一种理想的测量数据。在实际应用中，测量数据不可能不耦合一定的噪声。因此，在一定噪声背景下研究算法的重建效果很有必要。以下按公式(6)为边界测量数据添加噪声,取噪声强度系数分别为0、0.01%、0.05%、0.1%和0.5%，其中0表示不添加噪声，如图3各列所示。

在图3中，NOSER算法(“NS”)的重建图像在噪声强度达到0.05%时，其图像分辨率已显著下降；噪声强度达到0.5%时，从重建图像上已无法分辨目标，重建失败。灵敏度矩阵算法(“SE”)受噪声影响更为严重，在0.1%的噪声强度下就无法对扰动目标成像。反投影算法(“BP”)在噪声强度小于0.1%时，扰动目标均可被有效重建成像，即使噪声强度达到0.5%，EIT成像才略受影响，但从重建图像上仍可分辨出扰动目标。因此，就抗噪性能而言，反投影算法最佳，NOSER算法次之，灵敏度矩阵算法最差。

图3 不同噪声条件下反投影、NOSER和灵敏度矩阵算法的重建成像Fig.3 Images reconstructed by back-projection algorithm,newton’s one-step error reconstructor algorithm and sensitivity matrix algorithm with different noise conditions

3.3 重建速度的比较

在现有条件(Intel Pentium 4 CPU 2.8 G，256 M内存)和剖分规模(512个有限单元)下，进行一次图像重建，反投影算法耗时不足1秒，而NOSER算法和灵敏度矩阵算法至少需要十几分钟、甚至几十分钟，主要原因在于它们均须进行一个大矩阵的正则化和求逆运算。

当然，可以按照快速(fast)NOSER算法思想[7]，事先以文件方式保存雅可比矩阵和灵敏度矩阵，重建时直接读取矩阵数据而不再计算，以提高运算效率。但NOSER算法和灵敏度矩阵算法对电导率分布初值要求较高，事先计算好的相应矩阵在实际应用中难以使用。而反投影算法对电导率分布初值要求较低，却可以事先按均匀分布计算好反投影矩阵，从而进一步减小应用中图像重建时间。

4 结论

以上仿真研究结果表明：反投影算法的重建图像虽然图像分辨率不高，但其图像平滑，无明显伪影，且该算法具有较强的抗噪性能和快速的运算能力，在临床检测的强噪声背景下，应优选此算法。NOSER算法的图像分辨率最好，且抗噪性能适中，在对运算速度要求不高时可作为首选。灵敏度矩阵算法是一种整体性能适中的算法，可作为EIT成像的备选算法。

[1]RH Bayford.Bioimpedance tomography (electrical impedance tomography)[J].Annu Rev Biomed Eng,2006,8: 63-91.

[2]T Muders,H Luepschen,C Putensen.Impedance tomography as a new monitoring technique[J].Curr Opin Crit Care,2010,16(3):269-275.

[3]M Bodenstein,M David,K Markstaller.Principles of electrical impedance tomography and its clinical application[J].Crti Care Med,2009,37(2): 713-724.

[4]刘国强.医学电磁成像[M].北京: 科学出版社,2006

[5]DC Barrber,AD Seagar.Fast reconstruction of resistance images[J].Clin Phys Physiol Meas,1987,8(Suppl A): 47-54.

[6]M Cheney,D Isaacson,JC Newell,et al.NOSER: An algorithm for solving the inverse conductivity problem[J].Int J Imag Syst Technol,1990,2: 66-75.

[7]罗辞勇.基于快速牛顿一步误差重构的电阻抗成像算法和实验研究[D].重庆大学,博士学位论文,2005

[8]T Murai and Y Kagawa.Electrical impedance computed tomography based on a finite-element model[J].IEEE Trans Biomed Eng,1985,32: 177-184.