基于混合全變差正則化算法的電阻抗成像技術(shù)研究

左從磊， 李 靜

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

0 引 言

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

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

1 相關(guān)理論

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

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

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

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

(3)

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

2 算法推導(dǎo)

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

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

(4)

根據(jù)PDIPM理論，對偶(D)問題可表示為

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

(5)

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

原始問題是一個最小化問題，而對偶問題是一個最大化問題。通過消除其互補間隙GPD來求取可行域內(nèi)的最優(yōu)點

當且僅當式(8)成立，互補間隙GPD為0

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

(8)

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

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

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

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

gk+1=gk+λkδgk

(12)

為保證x的更新方向是對偶問題的下降方向，使用步長更新規(guī)則[8]

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

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

3 性能參數(shù)

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

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

(14)

算法的思路是平衡L1與L2罰項的阻尼作用，結(jié)合L2罰項的光滑和L1罰項的變分優(yōu)勢，重構(gòu)出分辨率高的圖像。算法在成像初期，L2罰項占據(jù)重構(gòu)的優(yōu)勢地位，而在成像后期L1罰項占據(jù)重構(gòu)的優(yōu)勢地位。在成像初期，考慮使逆問題的解快速收斂到真值的附近范圍;在成像后期，則考慮利用L1的變差特性增加邊界銳利度。一個重要的問題是正則化參數(shù)α的選擇，α過大或者過小都會對圖像的重構(gòu)結(jié)果產(chǎn)生影響。在實際成像中，α值具有較高的魯棒性[12]。基于先驗知識，α值在一定范圍內(nèi)的成像結(jié)果都有好的表現(xiàn)。本文HTV算法中α的較好范圍是10-4～10-6，仿真算例的重構(gòu)是在這個基礎(chǔ)上進行研究的。對于其他的算例，α的取值范圍會有所不同。權(quán)重參數(shù)η設(shè)定為[0,1]之間的等間距數(shù)，本文設(shè)置了20等分。選用三個指標以定量評判圖像重構(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的標準差，ρMN為協(xié)方差，表征N相對于M的非線性變化。SSIM(M,N)∈[0,1]，越接近1表明重構(gòu)質(zhì)量越好，其對于邊緣跳變有較好的區(qū)分能力。

4 仿真與討論

4.1 仿真平臺搭建

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

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

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

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

表1給出了不同算法重構(gòu)性能的定量對比。從表1和圖2結(jié)果來看，HTV算法的重構(gòu)結(jié)果在三個參數(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)尚好，但IRE值均在0.9以上，顯著高于HTV和TK法，SSIM值相比HTV法也非常小，究其原因是這兩種方法雖然可以一定程度上重構(gòu)分布范圍，但重構(gòu)電導(dǎo)率與真實值相差大。從總體上看，隨著目標區(qū)域的增大和形狀的復(fù)雜程度增加，重構(gòu)圖像相關(guān)系數(shù)都出現(xiàn)了一定的下降。

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

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

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

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

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

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

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

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

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

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

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

5 結(jié) 論

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