基于GMRES和Tikhonov正则化的生物电阻抗图像重建算法

生物医学工程学杂志

Journal of Biomedical E ngineering

V ol.26　N o.4

August　2009

基于G MRES和Tikhonov正则化的生物

电阻抗图像重建算法3

王化祥　范文茹Δ　胡　理

(天津大学电气与自动化工程学院,天津300072)

　　摘要　电阻抗层析成像(Electrical impedance tomography,EIT)是利用被测物体场内部电导率分布不均匀性,通过边界注入电流,测量边界电压变化,重构被测场内电导率分布图像。由于EIT测量数据有限,场域存在严重的非线性,导致问题的欠定性。我们介绍了一种新的组合算法,利用GMRES算法生成Krylov子空间,并结合Tik2 honov正则化方法进行图像重建。该算法不仅改善了实时性,而且提高了成像质量及鲁棒性。

关键词　电阻抗层析成像　Tikhonov正则化　Krylov子空间　组合算法

中图分类号　O441.4　　文献标识码　A　　文章编号　100125515(2009)0420701205

A H ybrid R econstruction Method in Electrical Impedance Tomography

B ased on G MRES and Tikhonov R egularization

W ang H uaxiang　F an Wenru　H u Li

(School of Elect rical Engineering and A utomation,Tianj in Universit y,Tianj in300072,China)

　　Abstract　Electrical impedance tomography(EIT)is a technique for reconstructing the conductivity distribution of measured field owing to its characteristics of being non2homogeneous,of injecting current at the boundary of the measured subject,of measuring the corresponding changes in voltage,and of reconstructing the image of the subject consequently.However,the limited measurement data of EIT,and the serious nonlinearity of the field result in ill2 posed problem,and the resolution of reconstructed image is poor.To solve the problem,a new hybrid algorithm is herein proposed.The method combines the characteristics of Krylov subspace and Tikhonov regularization,which can improve the real time performance,the quality and robustness of reconstructed image.

K ey w ords　Electrical impedance tomography(EIT)　　Tikhonov regularization　　Krylov subspace　　Hy2 brid algorithm

1　引　言

电阻抗层析成像(Elect rical impedance tomo2 grap hy,EIT)是近年来新兴的一种医学成像技术,它利用人体器官阻抗特性的不同,以及生理活动或病变导致的人体阻抗变化这一特点,将阵列电极配置于人体体表,通过对电极加载不同频率的交流激励信号,检测测量电极的电压信号,经过数据处理提取其中的电特性(电导率)信息,并通过图像重建算

3国家自然科学基金资助项目(60532020,50337020, 60472077,60301008);国家科技部支撑计划项目资助(2006BAI03A00)

Δ通讯作者。E2mail:wenrufan@hot mail.com 法,显示各组织器官截面电导率分布。与传统的成像技术相比,具有非侵入性、无损伤、可连续监测、价格低廉等优点。研究表明,EIT对于监测肺功能、检测乳腺癌、研究空腹状态等具有实用价值[1]。

EIT的图像重建实质上是利用边界测量电压值(V)重构物场内部电导率分布(σ)[2]。英国的Bar2 ber和Brown[3]最早提出利用X2射线断层成像技术中的反投影法求解场域内σ分布。尽管反投影法具有较好的鲁棒性,但是存在分辨率低,成像误差较大等缺点。目前常用的求解逆问题的算法可分为迭代成像算法和直接成像算法。迭代成像算法有Land2 weber迭代法[4]、共轭梯度法[5]等,但存在计算量较大、收敛速度较慢等问题;直接成像算法有Tik2　　　　　　　　　　　　　　　　　　　生物医学工程学杂志　　　　　　　　　　　　　　　　　　　第26卷

honov正则化算法[6]、截断奇异值分解(TSVD)[5]

等,如何选取合适的正则化参数成为问题的关键。

由于逆问题求解的难点在于其病态性及不确定

性,如何引入有效的算法提高重建计算的速度、精度

和稳定性是算法研究的重点。传统的Tikhonov正

则化算法,由于对解的滤波效应,使重构图像趋于平

滑,从而丧失了应有的对比度和锐度。GMRES

(Generalized Minimal Residual Algorit hm)算法[7]

是正则化方法中的一种投影方法,尽管根据有限的

子空间维度,可以获得正则化效果且计算速度快,但

是难以选取合适的维度。因此,本文引入了一种混

合算法GMRES2Tikhonov,它具有收敛速度快、精

度高、鲁棒性好等特点。

2　算　法

2.1　线性化物理模型

EIT线性化物理模型[8]可写为:

δU=U′(σ

)δσ=Jδσ(1)

式中:δσ为电导率变化;δU为电导率变化引起的边

界电压变化;J为J acobian矩阵[9]。

2.2　Tikhonov正则化方法

Tikhonov正则化方法是基于最小二乘原理,通

过加入罚函数,将病态性问题转化为良态问题,即

　min

σ

{‖Jδσ-δU‖22+λ2‖L(δσ-δσ3)‖22(2)

式中:δσ3为设定的初始值,一般取δσ3=0;L为正

则化矩阵,通常取L=I;λ为正则化参数,λ>0。

Tikhonov正则化问题可以转化为最小二乘问

题,

min

J

λL

δσ-

δU

(3)

　　通过最小化GCV[10]函数,获得λ的最优值。实验表明,此正则化参数选取方法具有较高的精确度和鲁棒性。

2.3　G MRES正则化算法

EIT逆问题经过线性化处理,变为典型的非对称线性问题。共轭梯度、Lanczos等传统的迭代算法,主要是解决无限维对称问题。因此,基于非对称问题,作者引入广义最小残差算法,它基于Krylov 子空间[11]产生的投影算法,求解速度快。

GMRES是基于Arnoldi方法[11]产生Krylov 子空间的基向量,并产生Hessenberg矩阵H k∈

R k×k,

H k=

h11h12……h

1k

h21h22 (2)

0h32ω…

…ωω…

0…0h k,k-1h kk

(4)

经过k步迭代,可得

J W k=W k H k+h k+1,k w k+1e R k(5)式中:W k=[w1,w2,…,w k]为Krylov子空间κk(J,

δU)的标准正交矩阵,‖w

i‖=1,i=1,2,…,k;e i= (0,…,0,1i,0,…,0)T∈R k+1。当h k+1,k=0时,迭代过程结束。

基于Krylov子空间求解方程(1)的近似解,可转化为最小二乘问题,

‖Jδσ(k)-δU‖=min

δσ

‖Jδσ-δU‖2,δσ(k)∈s p an (w1,w2,…,w k)(6)式中,δσ(k)(k≥1)为线性系统的近似解,初始解δσ(0) =0。

而近似解δσ(k)可写为

　　δσ(k)=W k y(k),s.t.y(k)=

arg min

y

‖(J W k)y-δ

i

U‖2,y(k)∈R k(7)　　经过Arnoldi过程的k步迭代,‖(J W k)y-δU‖2可转化为

‖(J W k)y-δU‖2=‖J W k y-βv1‖2=

‖W k+1[ H k y-βe1]‖2=

‖ H k y-βe1‖2(8)　　则

y(k)=argmin‖ H k y-βe1‖2(9)式中:e1=(1,0,…,0)T∈R k+1;W k+1为正交矩阵;

H k∈R(k+1)×k,通过Arnoldi Upper Hessenberg算法快速求解获得。因此,y(k)为原问题在子空间内的解,通过δσ(k)=W k y(k)可以得到原问题解。

投影法求解最小二乘问题,可有效减少条件数,增强了抗噪能力,而且保留了原系统的有效信息,减少了计算时间。但对Krylov子空间,关键是如何选取合适的维数k。维数选择小,成像效果不一定收敛到最佳值;维数选择大,重构成像发散。

2.4　G MRES2Tikhonov组合算法

207

第4期　　　　　　　　　王化祥等。　基于GMRES 和Tikhonov 正则化的生物电阻抗图像重建算法　　　　　　　　　

Tikhonov 正则化方法是通过增加罚函数滤除系统中无效信息,而GMRES 算法是基于将解空间投影到特殊的子空间,以滤除无效信息。二者有机

组合,将原问题投影到有效子空间中,通过增加罚项对子空间中的问题进行求解,并将求得的解映射回原空间,达到双重正则化的效果。

组合算法表示为

δσ(k )λ=W k y λ

(k )(10)

y λ

(k )

=arg min y

{‖(J W k )y -δU ‖22+λ2‖R y ‖2

2}

(11)

式中,R =L ×W k 为Krylov 子空间的正则化矩阵。

将式(8)代入式(11),得

y λ

(k )

=arg min y

{‖ H k y -βe 1‖2

2

+λ2

‖R y ‖22

}(12)

由于 H k 为上三角阵,上式时间复杂度为O (k ),而单独采用Tikho nov 正则化方法时间复杂度为O (min (m ,n )m n )(m :矩阵的行数,n :矩阵的列数)。可见组合算法有效减少了运算时间。

另一方面,式(2)取δ

σ3=0,得‖L δσλ(k )‖2=‖R W T

k W k y λ(k )‖2=‖R y λ(k )‖2,

W T

k W k =1

(13)‖J δσλ(k

)

-δU ‖2=‖ H k y -βe 1‖2

(14)

　　由上式可见,解泛数‖L δσλ(k )

‖2和残差泛数

‖J δσλ(k )-δU ‖i

2可由投影空间的Tikhonov 正则化解

y λ(k )

直接获得,而λ的最优值可通过GCV 方法求得。

3　实验结果

3.1　仿真实验

本文正问题采用COMSOL Multiphysics @和

Matlab @软件求解,COMSOL 具有建模方便及求解精确等优点。采用COMSOL 自动生成剖分网格,边界自动细分。逆问题采用方形网格进行求解计算。

本实验采用16电极模型,电流激励模式。采用两种不同的电导率分布,考虑实际测量中存在随机误差,因此加入±011%的高斯白噪声,分别用Landweber 迭代算法、TSVD 算法及Tikhonov 2GMRES 组合算法进行图像重建,如图1所示。从

成像效果可见组合算法优于前两种算法。

　　实验通过COMSOL 和Matlab 采用有限元方法进行正问题求解。由于仿真中实际电导率分布已知,引入相关系数对重建图像质量进行评判:

r =

∑N

i =1(σi - σ)(σi

3

- σ3)∑N

i =1

(σi - σ)2∑N

i =1

(σ3i - σ

3)2(15)

式中:σ3

为实际电导率分布;σ为电导率计算值; σ3

为实际电导率分布平均值; σ为电导率计算平均值。分别对共轭梯度算

法(C G )、Landweber 迭代、GMRES 以及组合算法进行相关计算,如图2所示。结果表明组合算法的稳定性高于单投影算法(GMRES ),其收敛速度优于C G 和Landweber 迭代算法。

图1　图像重建效果图

Fig 1　Im age reconstruction with different algorithms

307

　　　　　　　　　　　　　　　　　　　生物医学工程学杂志　　　　　　　　　　　　　　　　　　　第26卷

　　同时,相对成像误差‖σ-σ3‖‖σ3‖用于衡量成像精

度。对4种不同电导率分布分别采用几种常用方法

和组合算法进行比较,如表1所示。可见组合算法具有较高的精确度和稳定性

。

图2　不同算法的相关系数比较(3个物体)

Fig 2　Comp arison of correlation coeff icient for different algorithms

(3objects)

表1　几种算法对不同电导率分布成像的相对误差比较(迭代60次)

T able 1　Comp arison of some different algorithms in terms of relative

im age error (60

iterations)

此外,如上文所述, H k 为上三角结构,求解式(12)的

时间复杂度为O(k )。因此对大规模问题求解时,投影算法有效地减少了运算时间,从而提高了成像速度。对不同数量的网格剖分使用三种成像算法的计算时间进行比较(计算机配置:Intel (R)4,CPU 3106GH z ,内存512M),如表2所示。

表2　不同算法图像重建实时性比较

(其中,GMRES 和组合算法的子空间维数为60,Tikhonov 正则化参数由GCV 方法给出)T able 2　Comp arison of three different algorithms in terms of compu 2

tation time

(The dimension of subspace is 60and t he regularization parameters are automatically chosen using GCV met hod )

图像重建时间(s )

Tikhonov

GMRES 组合算法

316个单元格　0.493179

2.871236 2.595958812个单元格10.573979

3.261484 2.7006952828个单元格

580.165768

13.849

7.4581

3.2　人体肺部呼吸实验

基于人体肺部组织结构特性和电导率分布特性

(见表3)等先验知识,利用COMSOL 建立肺部模型,进行正向问题计算[12],生成J acobian 矩阵。

表3　各组织电导率分布

T able 3　Conductivity distribution of different

tissues

组织

心脏

肺

脊椎

脂肪

电导率σ(1/Ωm )

0.670.04220.138

0.006

0.037

　　人体实验中,将16个氯化银电极等间隔粘贴于

人体胸部,注入幅值为1mA 频率为500k Hz 的正弦电流,采用相邻激励、测量模式,将呼气末作为参考数据,采用动态成像法重建每个采样时刻的图像。采用组合算法对呼吸运动进行过程成像,如图3所示。

图3　肺部呼吸运动图像重建(a )～(e )吸气过程;(f )～(i )呼气过程

Fig 3　The im age reconstruction of hum an thorax during the respiration process

(a )2(e )inspiration process ;(f )2(i )expiration process

407

4　总　结

上述组合算法将Tikhonov正则化应用于Kry2 lov子空间以解决投影问题,由于奇异值向量在Krylov子空间中的前k个向量包含了J acobian矩阵的最优信息,从而有效地减少了计算量。由上述实验证明,组合算法不仅可以提高成像质量和精度;并且当选择合适的子空间维数时,组合算法的实时性和鲁棒性明显低于Tikhonov正则化法和GMRES算法。

参考文献

[1]CO H EN2BACRIE C,GOUSSARD Y and GUARDO R.Reg2

ularized reconstruction in electrical impedance tomography u2

sing a variance uniformization constraint[J].Medical Ima2

ging,IEEE Transactions on1997,16(5):5622571.

[2]毕德显.电磁场理论[M].北京:电子工业出版社,1985:652

232.

[3]BARBER D C and BROWN B H.Applied potential tomo2

graphy[J].J Phys E:Sci Instrum,1984;17:7232733. [4]WAN G H X,WAN G C and YIN W L.A pre2iteration met h2

od for t he inverse problem in elect rical impedance tomography

[J].Instrumentation and Measurement.IEEE Transactions

on,2004,53(4):109321096.

[5]YAN G W Q and PEN G L H.Image reconstruction algo2

rit hms for electrical capacitance tomography[J].Meas Sci

Technol,2003,14:R12R13.

[6]VAU H KON EN M,VADASZ D,KARJ ALAIN EN P A,et

al.Tikhonov regularization and prior information in electrical

impedance tomography[J].Medical Imaging,IEEE Transac2

tions on1998,17(2):2852293.

[7]SAAD Y and SCHUL TZ M H.GMRES:A generalized mini2

mal residual algorit hm for solving nonsymmetric linear sys2

tems[J].Society for Industry and Applied Mat hematics,

1986,7(3):8562869.

[8]VAU H KON EN M.Electrical impedance tomography and pri2

or information[D].(PhD)Finland:University of Kuopio,

1997:502.

[9]GESELOWITZ D B.An application of electrocaidiographic

lead t heory to impedance plet hysmography[J].IEEE Trans

Bioned Eng,1971,18(1):38241.

[10]HANSEN P C.Discrete inverse problems insight and algo2

rit hms.A tutorial wit h matlab exercises,in informatics and

mat hematical modelling[D].Denmark:Technical University

of Denmark,Denmark.2005:61272.

[11]J ACOBSEN M.Modular regularization algorit hms[D].

(PhD)Denmark:Technical University of Denmark,2004:72

24.

[12]HU L.Research on image reconst ruction algorit hm based on

‘soft2field’characteristics[D].(MSc)China:Tianjin Univer2

sity,2007:60271.

(收稿:2007208201　　修回:2008203220)

507