左從磊， 李 靜

(中國(guó)科學(xué)技術(shù)大學(xué) 工程科學(xué)學(xué)院，安徽 合肥 230027)

0 引 言

電阻抗成像(electrical impedance tomography,EIT)技術(shù)作為一種新型的成像技術(shù)，具有無(wú)輻射、非侵入性、成本低等優(yōu)勢(shì)，在醫(yī)學(xué)成像、地質(zhì)勘察、流體成像、材料檢測(cè)等領(lǐng)域都有廣闊的應(yīng)用前景[1,2]。EIT的目標(biāo)是重構(gòu)檢測(cè)對(duì)象的電導(dǎo)率，一般在邊界放置若干個(gè)測(cè)量電極，并向成像物體中注入一定幅值的安全激勵(lì)電流，通過(guò)測(cè)量并處理邊界電壓來(lái)重構(gòu)電導(dǎo)率分布。由于從邊界電壓獲得的數(shù)據(jù)數(shù)目遠(yuǎn)小于未知的電導(dǎo)率分布，EIT重構(gòu)結(jié)果具有不唯一性和不適定性；另一方面，測(cè)量容易受到噪聲的影響，邊界電壓的微小擾動(dòng)會(huì)導(dǎo)致重構(gòu)過(guò)程發(fā)生很大變化，使得EIT成為一個(gè)嚴(yán)重病態(tài)的逆問(wèn)題。常用的L2范數(shù)正則化算法使逆問(wèn)題的病態(tài)性得到改善，但由于L2范數(shù)的處處連續(xù)性質(zhì)，使重構(gòu)圖像的邊界不理想，同時(shí)重構(gòu)圖像往往具有比較嚴(yán)重的偽跡和重影，這些都制約了L2范數(shù)正則化的運(yùn)用[3～5]。近年來(lái)，L1范數(shù)的使用被越來(lái)越多的研究者所考慮[6,7]。

本文基于原始對(duì)偶內(nèi)點(diǎn)方法[8]，通過(guò)組合L1和L2范數(shù)，實(shí)現(xiàn)混合全變差(hybrid total variation,HTV)正則化重構(gòu)算法，有效提高了圖像重構(gòu)質(zhì)量，并在此基礎(chǔ)上進(jìn)一步研究了算法的抗噪聲性能。

1 相關(guān)理論

在各項(xiàng)研究中，共軛梯度(conjugate gradient，CG)法[9]、范數(shù)正則化等方法被廣泛應(yīng)用于EIT逆問(wèn)題求解。其中,常用的L2范數(shù)正則化方法通過(guò)添加罰函數(shù)項(xiàng)對(duì)解起到一定的阻尼作用，使逆問(wèn)題解保持穩(wěn)定，可以表示為

式中l(wèi)j為相鄰邊長(zhǎng)度，gm(j)和gn(j)分別為相鄰單元元素值。

將L1和L2正則化組成混合罰函數(shù)項(xiàng)，以發(fā)揮L1和L2范數(shù)的優(yōu)勢(shì)。定義混合HTV正則化算法的目標(biāo)函數(shù)為

α2‖L2(g-g0)‖2}

(3)

式中Lj為變差矩陣，為式(2)的矩陣形式，L2為NOSER正則化矩陣，α1和α2分別為L(zhǎng)1范數(shù)和L2范數(shù)的正則化參數(shù)。

2 算法推導(dǎo)

式(3)是一個(gè)凸優(yōu)化問(wèn)題，由于L1范數(shù)的不可微分性質(zhì)，求解需要引入平滑參數(shù)。而傳統(tǒng)的算法對(duì)于平滑參數(shù)的選擇十分苛刻，其值過(guò)大或者過(guò)小對(duì)圖像的質(zhì)量都有很大影響[12]。Anderson K D等人[8]已經(jīng)證明原始對(duì)偶內(nèi)點(diǎn)法(primal-dual interior-point method,PDIPM)在求解L1范數(shù)上具有巨大優(yōu)勢(shì)。PDIPM算法通過(guò)引入對(duì)偶變量將不等式約束轉(zhuǎn)化為等式約束[7]，對(duì)平滑參數(shù)的選擇具有高的魯棒性。本文基于PDIPM推導(dǎo)了HTV算法的求解框架如下。把式(3)標(biāo)記為原始(P)問(wèn)題

α2‖L2(g-g0)‖2}

(4)

根據(jù)PDIPM理論，對(duì)偶(D)問(wèn)題可表示為

α2‖L2(g-g0)‖2]},‖x‖≤1

(5)

式中x為引入的對(duì)偶變量，L1為二維EIT有限元問(wèn)題的L1范數(shù)矩陣。為得到內(nèi)部最小化條件，式(5)對(duì)g求一階偏導(dǎo)數(shù)并令其值為0

原始問(wèn)題是一個(gè)最小化問(wèn)題，而對(duì)偶問(wèn)題是一個(gè)最大化問(wèn)題。通過(guò)消除其互補(bǔ)間隙GPD來(lái)求取可行域內(nèi)的最優(yōu)點(diǎn)

當(dāng)且僅當(dāng)式(8)成立，互補(bǔ)間隙GPD為0

Ljg-xi|Ljg|=0,?i,j

(8)

如果相鄰單元的元素值相同，Ljg為0，上述條件不可微分。引入平滑參數(shù)β，克服不可微的問(wèn)題

使用Newton方法進(jìn)行計(jì)算。式(6)和式(9)分別對(duì)原始變量g和對(duì)偶變量x求一階偏導(dǎo)數(shù)。為方便描述，式(6)記為Eq6，式(9)記為Eq9

最終得到HTV算法的PDIPM求解框架如式(11)

基于式(11)可以對(duì)原始變量和對(duì)偶變量進(jìn)行迭代更新求解。g的更新步長(zhǎng)使用優(yōu)化的線搜索方法，其中,k為迭代次數(shù)，λk為線搜索得到的步長(zhǎng)因子

gk+1=gk+λkδgk

(12)

為保證x的更新方向是對(duì)偶問(wèn)題的下降方向，使用步長(zhǎng)更新規(guī)則[8]

xk+1=xk+min(1,φ*)δxk

式中φ*為第k次迭代時(shí)xi與可行域邊界+1和-1之間的最小距離，N為x的維度。

3 性能參數(shù)

如式(3)所示，g的更新取決于三個(gè)部分：基于二次函數(shù)的殘差項(xiàng)、L1罰函數(shù)項(xiàng)以及L2罰函數(shù)項(xiàng)。對(duì)L1和L2罰項(xiàng)的系數(shù)做一定調(diào)整，可以改變罰項(xiàng)對(duì)逆問(wèn)題解的懲罰作用。通過(guò)引入權(quán)重系數(shù)對(duì)罰項(xiàng)系數(shù)做如下的變形操作

α(η|L1g|+(1-η)‖L2(g-g0)‖2) }

(14)

算法的思路是平衡L1與L2罰項(xiàng)的阻尼作用，結(jié)合L2罰項(xiàng)的光滑和L1罰項(xiàng)的變分優(yōu)勢(shì)，重構(gòu)出分辨率高的圖像。算法在成像初期，L2罰項(xiàng)占據(jù)重構(gòu)的優(yōu)勢(shì)地位，而在成像后期L1罰項(xiàng)占據(jù)重構(gòu)的優(yōu)勢(shì)地位。在成像初期，考慮使逆問(wèn)題的解快速收斂到真值的附近范圍;在成像后期，則考慮利用L1的變差特性增加邊界銳利度。一個(gè)重要的問(wèn)題是正則化參數(shù)α的選擇，α過(guò)大或者過(guò)小都會(huì)對(duì)圖像的重構(gòu)結(jié)果產(chǎn)生影響。在實(shí)際成像中，α值具有較高的魯棒性[12]?；谙闰?yàn)知識(shí)，α值在一定范圍內(nèi)的成像結(jié)果都有好的表現(xiàn)。本文HTV算法中α的較好范圍是10-4～10-6，仿真算例的重構(gòu)是在這個(gè)基礎(chǔ)上進(jìn)行研究的。對(duì)于其他的算例，α的取值范圍會(huì)有所不同。權(quán)重參數(shù)η設(shè)定為[0,1]之間的等間距數(shù)，本文設(shè)置了20等分。選用三個(gè)指標(biāo)以定量評(píng)判圖像重構(gòu)的質(zhì)量，一是圖像重構(gòu)誤差I(lǐng)RE，二是圖像相關(guān)系數(shù)ICC，三是結(jié)構(gòu)相似度(structural similarity,SSIM)。

1)圖像重構(gòu)誤差(image reconstruction error,IRE)

IRE=‖G-‖/‖G‖

(15)

式中G為設(shè)定電導(dǎo)率，為重構(gòu)電導(dǎo)率。IRE越小，表明重構(gòu)誤差越小。

2)圖像相關(guān)系數(shù)(image correlation coefficient,ICC)

3)SSIM

式中uM和uN分別為模型M和重構(gòu)結(jié)果N的均值，ρM和ρN分別為M和N的標(biāo)準(zhǔn)差，ρMN為協(xié)方差，表征N相對(duì)于M的非線性變化。SSIM(M,N)∈[0,1]，越接近1表明重構(gòu)質(zhì)量越好，其對(duì)于邊緣跳變有較好的區(qū)分能力。

4 仿真與討論

4.1 仿真平臺(tái)搭建

基于MATLAB平臺(tái)開(kāi)發(fā)二維圓域仿真算例，將背景區(qū)域電導(dǎo)率設(shè)定為1 S/m，目標(biāo)區(qū)域設(shè)定為2 S/m。為測(cè)試算法對(duì)不同目標(biāo)的重構(gòu)情況，建立了單目標(biāo)模型、雙目標(biāo)模型、多目標(biāo)模型，不同模型的目標(biāo)形狀和位置分布也不盡相同。為保證測(cè)量條件相同，均采用16電極模型，激勵(lì)電流幅值為1 mA，激勵(lì)模式和測(cè)量模式均為相鄰模式。通過(guò)二維仿真數(shù)值模擬，測(cè)試HTV算法在電導(dǎo)率分布重構(gòu)中的性能表現(xiàn)，并與TK法、CG法、NOSER法進(jìn)行了比較。

4.2 無(wú)噪聲情況下重構(gòu)成像

圖1對(duì)比了無(wú)噪聲時(shí)HTV算法以及幾種經(jīng)典的重構(gòu)算法的性能表現(xiàn)。HTV算法在單目標(biāo)、雙目標(biāo)、多目標(biāo)模型中都展現(xiàn)出更好的重構(gòu)效果。1)HTV算法的成像邊界最清晰；Tikhonov法、CG法和NOSER法的成像邊界模糊，證明L1范數(shù)在重構(gòu)空間階躍性變化時(shí)具有獨(dú)特的優(yōu)勢(shì)。2)HTV算法的形狀重構(gòu)最準(zhǔn)確;三種模型的目標(biāo)區(qū)域及形狀被都被較好地定位和重構(gòu)，相比之下其余方法形變較為嚴(yán)重。3)HTV算法的成像偽跡最少;TK法和CG法的鋸齒較嚴(yán)重，NOSER法的重構(gòu)目標(biāo)范圍過(guò)大且偽跡多。

圖1 無(wú)噪聲時(shí)不同算法的重構(gòu)圖像

表1給出了不同算法重構(gòu)性能的定量對(duì)比。從表1和圖2結(jié)果來(lái)看，HTV算法的重構(gòu)結(jié)果在三個(gè)參數(shù)上都表現(xiàn)最好，其次是TK算法。HTV算法三種模型的重構(gòu)圖像IRE值分別是TK法的約78 %，58 %，95 %；ICC值分別是TK法的約105 %，121 %，107 %；SSIM值分別是TK法的約110 %，128 %，113 %。而CG法和NOSER法重構(gòu)在ICC上表現(xiàn)尚好，但I(xiàn)RE值均在0.9以上，顯著高于HTV和TK法，SSIM值相比HTV法也非常小，究其原因是這兩種方法雖然可以一定程度上重構(gòu)分布范圍，但重構(gòu)電導(dǎo)率與真實(shí)值相差大。從總體上看，隨著目標(biāo)區(qū)域的增大和形狀的復(fù)雜程度增加，重構(gòu)圖像相關(guān)系數(shù)都出現(xiàn)了一定的下降。

表1 無(wú)噪聲時(shí)不同算法的重構(gòu)性能參數(shù)

圖2 無(wú)噪聲時(shí)不同算法的重構(gòu)性能參數(shù)柱狀圖

4.3 含噪聲情況下重構(gòu)成像

在實(shí)際應(yīng)用中，噪聲是不可避免的。為了測(cè)試算法的抗噪聲性能，在正問(wèn)題生成的測(cè)量電壓中加入了SNR為40 dB的隨機(jī)高斯噪聲。

圖3和表2給出了當(dāng)SNR為40 dB時(shí)算法的重構(gòu)結(jié)果對(duì)比。NOSER算法的抗噪性能最好，與理想無(wú)噪聲情況下相比重構(gòu)圖像變化最小，三種模型的重構(gòu)結(jié)果中ICC高于其余算法。而HTV算法因?yàn)榻Y(jié)合了NOSER正則化罰項(xiàng)，抗噪性能較好。相較于無(wú)噪聲情況，當(dāng)目標(biāo)靠近邊界時(shí)出現(xiàn)了一定的形變，但是HTV算法的IRE值在幾種算法中最小，抑制了偽跡且有清晰的邊界。在醫(yī)學(xué)成像領(lǐng)域，不同醫(yī)學(xué)組織之間的邊界是具有價(jià)值的圖像信息[6]。從這一點(diǎn)上而言，HTV算法以抗噪性能略微下降的代價(jià)獲得更清晰的圖像是值得的。Tikhonov法的重構(gòu)結(jié)果表現(xiàn)不佳，模型2和模型3的重構(gòu)目標(biāo)區(qū)域已經(jīng)很難分辨，圖像的偽跡和鋸齒比較嚴(yán)重。CG法表現(xiàn)不佳，其重構(gòu)電導(dǎo)率值與真實(shí)值差異大，相關(guān)系數(shù)和結(jié)構(gòu)相似度均較小。

圖3 含噪聲時(shí)不同算法的重構(gòu)圖像

表2 含噪聲時(shí)不同算法的重構(gòu)性能參數(shù)

4.4 真實(shí)肺部圖像重構(gòu)

人體在進(jìn)行肺部呼吸時(shí)，電導(dǎo)率會(huì)發(fā)生一定的變化。Andler A等人[13]采集了人體呼吸時(shí)肺部EIT數(shù)據(jù)。16個(gè)測(cè)量電極被等間距地放置在一個(gè)成年人的胸部周?chē)?，并以相鄰模式進(jìn)行電流激勵(lì)和電壓測(cè)量。幅值為10 mA的激勵(lì)電流通過(guò)相鄰的電極對(duì)注入，并在其余相鄰的電極之間測(cè)量出電壓差。采用Andler A等人采集的數(shù)據(jù)，應(yīng)用HTV算法和其余算法對(duì)肺部圖像進(jìn)行重構(gòu)。

圖4是不同算法的重構(gòu)圖像結(jié)果對(duì)比，結(jié)果表明HTV算法相較于其他方法成像質(zhì)量最高，重構(gòu)圖像邊界清晰、偽跡較少，形狀最接近真實(shí)肺部。因此HTV算法提高了肺部重構(gòu)圖像的重構(gòu)質(zhì)量，在實(shí)際應(yīng)用中有顯著優(yōu)勢(shì)和良好的應(yīng)用前景。

圖4 真實(shí)肺部圖像重構(gòu)圖像結(jié)果對(duì)比

5 結(jié) 論

本文創(chuàng)新使用L2范數(shù)的NOSER罰函數(shù)項(xiàng)和L1罰函數(shù)項(xiàng)的HTV正則化重構(gòu)算法，并基于原始對(duì)偶內(nèi)點(diǎn)法推導(dǎo)出算法的求解框架。通過(guò)建立不同大小和形狀的仿真測(cè)試模型，在EIT正問(wèn)題解的基礎(chǔ)上，使用HTV算法實(shí)現(xiàn)了圖像重建，并將其與經(jīng)典Tikhonov、共軛梯度、牛頓一步誤差重構(gòu)算法進(jìn)行了比較。結(jié)果表明：HTV算法在理想無(wú)噪聲情況下的重構(gòu)表現(xiàn)優(yōu)于其他算法，重構(gòu)圖像的邊界更為清晰且形狀最準(zhǔn)確。圖像的重構(gòu)誤差I(lǐng)RE最小，ICC和SSIM最高。在含噪聲情況下(SNR為40 dB)時(shí)HTV算法也有良好的抗噪性能，可以重構(gòu)出高質(zhì)量的圖像?；谌梭w肺部呼吸EIT采集數(shù)據(jù)進(jìn)行重構(gòu)，證明HTV算法重構(gòu)成像清晰，有較好的實(shí)用價(jià)值和應(yīng)用前景。