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

2009年　8月JournalofBiomedicalEngineeringAugust　2009

基于GMRES和Tikhonov正则化的生物

电阻抗图像重建算法3

王化祥　范文茹Δ　胡　理

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

　　摘要　电阻抗层析成像(Electricalimpedancetomography,EIT)是利用被测物体场内部电导率分布不均匀性,通过边界注入电流,测量边界电压变化,重构被测场内电导率分布图像。由于EIT测量数据有限,场域存在严重的非线性,导致问题的欠定性。我们介绍了一种新的组合算法,利用GMRES算法生成Krylov子空间,并结合Tik2

honov正则化方法进行图像重建。该算法不仅改善了实时性,而且提高了成像质量及鲁棒性。

关键词　电阻抗层析成像　Tikhonov正则化　Krylov子空间　组合算法中图分类号　O441.4　　文献标识码　A　　文章编号　100125515(2009)0420701205

AHybridReconstructionMethodinElectricalImpedanceTomography

BasedonGMRESandTikhonovRegularization

WangHuaxiang　FanWenru　HuLi

(SchoolofElectricalEngineeringandAutomation,TianjinUniversity,Tianjin300072,China)

　　Abstract　Electricalimpedancetomography(EIT)isatechniqueforreconstructingtheconductivitydistribution

ofmeasuredfieldowingtoitscharacteristicsofbeingnon2homogeneous,ofinjectingcurrentattheboundaryofthemeasuredsubject,ofmeasuringthecorrespondingchangesinvoltage,andofreconstructingtheimageofthesubjectconsequently.However,thelimitedmeasurementdataofEIT,andtheseriousnonlinearityofthefieldresultinill2posedproblem,andtheresolutionofreconstructedimageispoor.Tosolvetheproblem,anewhybridalgorithmishereinproposed.ThemethodcombinesthecharacteristicsofKrylovsubspaceandTikhonovregularization,whichcanimprovetherealtimeperformance,thequalityandrobustnessofreconstructedimage.

Keywords　Electricalimpedancetomography(EIT)　　Tikhonovregularization　　Krylovsubspace　　Hy2bridalgorithm

1　引　言

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

3国家自然科学基金资助项目(60532020,50337020,

60472077,60301008);国家科技部支撑计划项目资助(2006BAI03A00)

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

EIT的图像重建实质上是利用边界测量电压值(V)重构物场内部电导率分布(σ)[2]。英国的Bar2ber和Brown[3]最早提出利用X2射线断层成像技术

中的反投影法求解场域内σ分布。尽管反投影法具有较好的鲁棒性,但是存在分辨率低,成像误差较大等缺点。目前常用的求解逆问题的算法可分为迭代成像算法和直接成像算法。迭代成像算法有Land2weber迭代法[4]、共轭梯度法[5]等,但存在计算量较

Δ通讯作者。E2mail:wenrufan@hotmail.com

大、收敛速度较慢等问题;直接成像算法有Tik2

702　　　　　　　　　　　　　　　　　　　生物医学工程学杂志　　　　　　　　　　　　　　　　　　　第26卷honov正则化算法[6]、截断奇异值分解(TSVD)[5]Rk×k,

h11h21

h12h22h32

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

由于逆问题求解的难点在于其病态性及不确定性,如何引入有效的算法提高重建计算的速度、精度和稳定性是算法研究的重点。传统的Tikhonov正则化算法,由于对解的滤波效应,使重构图像趋于平滑,从而丧失了应有的对比度和锐度。GMRES(GeneralizedMinimalResidualAlgorithm)算法[7]是正则化方法中的一种投影方法,尽管根据有限的子空间维度,可以获得正则化效果且计算速度快,但是难以选取合适的维度。因此,本文引入了一种混合算法GMRES2Tikhonov,它具有收敛速度快、精度高、鲁棒性好等特点。

Hk=

……ωω

0

……ω

hk,k-1

h1kh2k

0……hkk

(4)

…0

…

经过k步迭代,可得

JWk=WkHk+hk+1,kwk+1ek

R

(5)

式中:Wk=[w1,w2,…,wk]为Krylov子空间κk(J,δU)的标准正交矩阵,‖wi‖=1,i=1,2,…,k;ei=

(0,…,0,1i,0,…,0)

T

∈Rk+1。当hk+1,k=0时,迭代

过程结束。2　算　法

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

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

(k)(k)

σδσ-δU‖σ‖Jδ-δU‖=min‖J∈span2,δδσ

δU=U′(σσ=Jδσ0)δ

(1)(w1,w2,…,wk)(6)

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

界电压变化;J为Jacobian矩阵[9]。2.2　Tikhonov正则化方法

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

()

σ(0)式中δ,σk(k≥1)为线性系统的近似解,初始解δ

=0。

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

(k)()()

δσ　　=Wkyk,s.t.yk=

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

222

δσ-δU‖δσ-δσ3)‖　min{‖JL(2+λ‖2

σ

(2)

argmin‖(JWk)y-δU‖2,y

y

i

(k)

∈Rk(7)

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

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

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

33

　　经过Arnoldi过程的k步迭代,‖(JWk)y-δU‖2可转化为

‖(JWk)y-δU‖JWky-βv1‖2=‖2=

‖Wk+1[H󰁵ky-βe1]‖2=

‖H󰁵ky-βe1‖2

　　则

y

(k)

题,

min

J

λL

[10]

δσ-

δU0

(3)

(8)(9)

　　通过最小化GCV和鲁棒性。

函数,获得λ的最优值。实

H󰁵k∈R

=argmin‖H󰁵ky-βe1‖2

T

验表明,此正则化参数选取方法具有较高的精确度2.3　GMRES正则化算法

EIT逆问题经过线性化处理,变为典型的非对

式中:e1=(1,0,…,0)

(k+1)×k

∈Rk+1;Wk+1为正交矩阵;

()

,通过ArnoldiUpperHessenberg算

法快速求解获得。因此,yk为原问题在子空间内的

σ(k)=Wky(k)可以得到原问题解。解,通过δ

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

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

GMRES是基于Arnoldi方法[11]产生Krylov

子空间的基向量,并产生Hessenberg矩阵Hk∈

703第4期　　　　　　　　　王化祥等。　基于GMRES和Tikhonov正则化的生物电阻抗图像重建算法　　　　　　　　　Tikhonov正则化方法是通过增加罚函数滤除yλ直接获得,而λ的最优值可通过GCV方法求得。

(k)

系统中无效信息,而GMRES算法是基于将解空间投影到特殊的子空间,以滤除无效信息。二者有机

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

组合算法表示为

(k)(k)δσλ=Wkyλ

(k)

3　实验结果

3.1　仿真实验

本文正问题采用COMSOLMultiphysics@和

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

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

(10)

yλ

222

(JWk)y-δ=argmin{‖U‖Ry‖2+λ‖2}

y

(11)

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

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

yλ

(k)

=argmin{‖H󰁵ky-βe1‖+λ‖Ry‖}(12)y

2

2

2

22

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

σ3=0,得另一方面,式(2)取δ

(k)(k)(k)T

δσ‖Lλ‖RWkWkyλ‖Ryλ‖2=‖2=‖2,

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

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

r=

3

i=1N3

∑(σ󰁡)(σ-σ󰁡3)i-σi

N

N

(15)

i=1

3∑(σ󰁡)2∑(σ󰁡3)2i-σi-σ

i=1

WkWk=1

k

δσ‖Jλ-δU‖󰁵ky-βe1‖2=‖H2

()

T

(13)(14)

(k)

δσλ　　由上式可见,解泛数‖L‖2和残差泛数

(k)i

δσ‖Jλ-δU‖2可由投影空间的Tikhonov正则化解

σ为实际电导率分布;σ为电导率计算值;σ式中:󰁡3

为实际电导率分布平均值;σ󰁡为电导率计算平均值。分别对共轭梯度算法(CG)、Landweber迭代、GMRES以及组合算法进行相关计算,如图2所示。结果表明组合算法的稳定性高于单投影算法(GMRES),其收敛速度优于CG和Landweber迭代算法。

图1　图像重建效果图

Fig1　Imagereconstructionwithdifferentalgorithms

704　　　　　　　　　　　　　　　　　　　生物医学工程学杂志　　　　　　　　　　　　　　　　　　　第26卷σ-σ3‖‖　　同时,相对成像误差用于衡量成像精

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

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

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

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

(其中,GMRES和组合算法的子空间维数为60,Tikhonov正则化参数由GCV方法给出)Table2　Comparisonofthreedifferentalgorithmsintermsofcompu2

tationtime

(Thedimensionofsubspaceis60andtheregularizationparametersareautomaticallychosenusingGCVmethod)

图像重建时间(s)

Tikhonov316个单元格812个单元格2828个单元格GMRES2.8712363.26148413.849组合算法

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

Fig2　Comparisonofcorrelationcoefficientfordifferentalgorithms

(3objects)

3.2　人体肺部呼吸实验

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

Table1　Comparisonofsomedifferentalgorithmsintermsofrelative

imageerror(60iterations)

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

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

表3　各组织电导率分布

Table3　Conductivitydistributionofdifferenttissues

组织

电导率σ(1/Ωm)

心脏

0.67肺

0.04220.138

脊椎

0.006

脂肪

0.037

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

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

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

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

Fig3　Theimagereconstructionofhumanthoraxduringtherespirationprocess

(a)2(e)inspirationprocess;(f)2(i)expirationprocess

705第4期　　　　　　　　　王化祥等。　基于GMRES和Tikhonov正则化的生物电阻抗图像重建算法　　　　　　　　　rithmsforelectricalcapacitancetomography[J].MeasSciTechnol,2003,14:R12R13.

[6]

VAUHKONENM,VADASZD,KARJALAINENPA,etal.Tikhonovregularizationandpriorinformationinelectricalimpedancetomography[J].MedicalImaging,IEEETransac2tionson1998,17(2):2852293.

[7]

SAADYandSCHULTZMH.GMRES:Ageneralizedmini2malresidualalgorithmforsolvingnonsymmetriclinearsys2tems[J].SocietyforIndustryandAppliedMathematics,1986,7(3):8562869.

[8]

VAUHKONENM.Electricalimpedancetomographyandpri2orinformation[D].(PhD)Finland:UniversityofKuopio,1997:502.

[9]

GESELOWITZDB.Anapplicationofelectrocaidiographicleadtheorytoimpedanceplethysmography[J].IEEETransBionedEng,1971,18(1):38241.

[10]HANSENPC.Discreteinverseproblemsinsightandalgo2rithms.Atutorialwithmatlabexercises,ininformaticsandmathematicalmodelling[D].Denmark:TechnicalUniversityofDenmark,Denmark.2005:61272.

[11]

JACOBSENM.24.

[12]

HUL.Researchonimagereconstructionalgorithmbasedon

Modularregularizationalgorithms[D].

(PhD)Denmark:TechnicalUniversityofDenmark,2004:72

4　总　结

上述组合算法将Tikhonov正则化应用于Kry2

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

参

[1]

考文献

COHEN2BACRIEC,GOUSSARDYandGUARDOR.Reg2ularizedreconstructioninelectricalimpedancetomographyu2singavarianceuniformizationconstraint[J].MedicalIma2ging,IEEETransactionson1997,16(5):5622571.[2][3][4]

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

232.

BARBERDCandBROWNBH.Appliedpotentialtomo2graphy[J].JPhysE:SciInstrum,1984;17:7232733.WANGHX,WANGCandYINWL.Apre2iterationmeth2odfortheinverseprobleminelectricalimpedancetomography[J].InstrumentationandMeasurement.IEEETransactionson,2004,53(4):109321096.

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

sity,2007:60271.

(收稿:2007208201　　修回:2008203220)

[5]YANGWQandPENGLH.Imagereconstructionalgo2