基于非精確Newton迭代正則化的EIT圖像重構(gòu)算法

趙一帆,王 靜

(黑龍江大學(xué) 數(shù)學(xué)科學(xué)學(xué)院,哈爾濱 150080)

0 引 言

電阻抗成像(Electrical impedance tomography,EIT)是一種新興的無損功能成像技術(shù),根據(jù)目標(biāo)體表面的邊界測(cè)量數(shù)據(jù)來估計(jì)目標(biāo)體內(nèi)部的電導(dǎo)率分布情況,具有非侵入性、設(shè)備輕便、成本低廉以及無損檢測(cè)等特點(diǎn),在醫(yī)療診斷、地理探測(cè)和工業(yè)檢測(cè)等領(lǐng)域有著廣泛的應(yīng)用前景。

圖像重構(gòu)是EIT的核心技術(shù),直接影響著成像的空間分辨率和實(shí)時(shí)性。EIT圖像重構(gòu)問題本身存在嚴(yán)重的不適定性,即觀測(cè)數(shù)據(jù)的微小變化會(huì)導(dǎo)致重構(gòu)參數(shù)發(fā)生較大的變化,而有限的測(cè)量數(shù)據(jù)以及噪聲的干擾都加劇了這種不適定性。為克服重構(gòu)問題的不適定性,往往需要借助不同的正則化技術(shù)來提高重構(gòu)圖像的質(zhì)量。國(guó)際上關(guān)于EIT圖像重構(gòu)方法的研究非常多,大致可分為3類：Tikhonov正則化為代表的直接方法、Landweber迭代為代表的迭代正則化方法,以及神經(jīng)網(wǎng)絡(luò)為代表的智能方法。Tikhonov正則化[1-2]是克服EIT重構(gòu)問題不適定性和病態(tài)性的有效方法,但正則化參數(shù)的選取較為復(fù)雜,直接影響成像效果,實(shí)際一般采用經(jīng)驗(yàn)值,具有一定的局限性。神經(jīng)網(wǎng)絡(luò)[3-5]是智能方法中研究較多的,也是當(dāng)前研究的熱點(diǎn),可直接建立邊界測(cè)量電壓數(shù)據(jù)與電導(dǎo)率分布之間的非線性映射關(guān)系,但其需要大量的樣本集來訓(xùn)練網(wǎng)絡(luò)模型,對(duì)訓(xùn)練樣本的完整性要求較高,且對(duì)未知模型無法預(yù)測(cè),在一定程度上也限制了其應(yīng)用。

迭代正則化方法中,最為經(jīng)典的Landweber迭代(Landweber iteration,LDI)[6-8]因其在實(shí)現(xiàn)圖像重構(gòu)的過程中只需用到數(shù)據(jù)擬合的梯度信息,易于實(shí)現(xiàn)且計(jì)算成本低,穩(wěn)定性好,也是目前圖像重構(gòu)中應(yīng)用最為成功的一種迭代方法。然而,靈敏度矩陣的病態(tài)性使得LDI迭代存在收斂速度較慢的不足,需要大量的迭代步來實(shí)現(xiàn)高質(zhì)量的成像,影響了其應(yīng)用,因此需要尋找更加快速的方法。近來,備受關(guān)注的非精確Newton迭代正則化方法[9-13]具有較快的收斂速度,是一種先線性化后正則化的方法。但需注意的是,線性化過程并不能消除問題本身的不適定性,仍需借助正則化策略處理局部線性化方程。該方法包括內(nèi)外兩層迭代：內(nèi)層迭代應(yīng)用迭代正則化方法求解局部線性化問題產(chǎn)生近似迭代序列,外層迭代為非精確Newton法用于更新迭代點(diǎn)列。穩(wěn)定性較好且格式簡(jiǎn)單的LDI迭代法[14]常常被用作為內(nèi)層迭代正則化策略,但鑒于LDI迭代法收斂速度較慢的缺陷,若將其作為內(nèi)層迭代策略,外層每一步迭代往往需要大量的內(nèi)層迭代次數(shù)來滿足停止準(zhǔn)則,造成內(nèi)層迭代計(jì)算量較大。2009年,韓波等首次提出將同倫攝動(dòng)迭代(Homotopy perturbation iteration,簡(jiǎn)稱HPI)用于求解非線性不適定反問題[15]。隨后,HPI迭代法被應(yīng)用于多個(gè)領(lǐng)域,如電阻抗成像[16]、地震勘探測(cè)井約束反演[17]等,數(shù)值結(jié)果表明：當(dāng)達(dá)到相似重構(gòu)精度時(shí),HPI迭代法僅需LDI迭代法大約一半的計(jì)算時(shí)間。近來,HPI迭代法作為一類加速策略得到了推廣應(yīng)用[18-19]?；诖?本文考慮采用HPI迭代作為內(nèi)層迭代并結(jié)合適當(dāng)?shù)牟介L(zhǎng)選取策略來加速內(nèi)層迭代速度,從而獲得一種更高效的非精確Newton迭代正則化方法。

1 EIT數(shù)學(xué)模型

目前廣泛采用的EIT數(shù)學(xué)模型為比較貼近實(shí)際的全電極數(shù)學(xué)模型[20]：

(1)

EIT問題的研究主要包括正問題和反問題。EIT正問題可歸結(jié)為求解邊值問題(1),即給定內(nèi)部電導(dǎo)率分布、注入電流和接觸阻抗,計(jì)算內(nèi)部電位u和電極電位U。正問題的求解可借助于有限元方法[21]對(duì)邊值問題(1)進(jìn)行離散化,可得到如下離散形式：

A(σ)b=f

(2)

其中,A(σ)為有限元離散矩陣,b為待定系數(shù) (用于計(jì)算u和U),f為包含電流激勵(lì)源的向量。

然而,在EIT實(shí)際應(yīng)用中,電導(dǎo)率分布σ是未知的,僅已知的是電極處的注入電流數(shù)據(jù)和測(cè)量電壓數(shù)據(jù),而測(cè)量電壓數(shù)據(jù)在數(shù)據(jù)采集過程中又不可避免地含有噪聲,這里記為Uδ。EIT觀測(cè)模型響應(yīng)可表示為如下的非線性算子方程

Uδ=F(σ)+e

(3)

這里,F表示模型空間到測(cè)量數(shù)據(jù)空間的投影,即求解正問題的過程,也稱之為正演算子,e表示測(cè)量誤差向量。因此,根據(jù)邊界電極處測(cè)量電壓數(shù)據(jù)Uδ反推內(nèi)部電導(dǎo)率分布σ的過程就稱之為EIT反問題,由于需要通過圖像的形式呈現(xiàn)出來,故也稱之為圖像重構(gòu)。

2 非精確Newton迭代正則化方法

2.1 LDI迭代法和HPI迭代法

EIT圖像重構(gòu)問題本身存在嚴(yán)重的不適定性,最小二乘法是處理該問題最為流行的方法,即極小化數(shù)據(jù)擬合項(xiàng)

(4)

最為經(jīng)典的LDI迭代格式為

σn+1=σn+μnF′(σn)*(Uδ-F(σn))

(5)

可視為(4)的最速下降法;HPI迭代法則是通過由同倫方法求解(4)的歐拉方程F′(σ)*(F(σ)-Uδ)=0得到的,其迭代格式為

σn+1=σn+μnF′(σn)*(2I-F′(σn)F′(σn)*)(Uδ-F(σn))

(6)

式中μn為適當(dāng)選取的步長(zhǎng)因子,這里采用比較適宜在EIT領(lǐng)域使用的變步長(zhǎng)因子選取準(zhǔn)則,滿足

(7)

2.2 INLDI迭代法和INHPI迭代法

鑒于非精確Newton迭代正則化方法具有較快的收斂速度,考慮將其用于處理EIT圖像重構(gòu)問題。該方法包括內(nèi)外兩層迭代：外層迭代為非精確Newton方法且以偏差原則作為停止準(zhǔn)則,內(nèi)層迭代這里考慮采用LDI迭代或HPI迭代并結(jié)合選取適當(dāng)?shù)牟介L(zhǎng)來加速。

首先構(gòu)造內(nèi)層迭代。假設(shè)已知當(dāng)前迭代步為σn,通過極小化問題

(8)

構(gòu)造相應(yīng)的迭代格式

σn,k+1=σn,k+μn,kdn,k

(9)

(10)

若選取

(11)

則為L(zhǎng)DI迭代格式;若選取

(12)

‖sn,k‖≤γ‖Uδ-F(σn)‖

(13)

‖Uδ-F(σn*)‖≤τδ≤‖Uδ-F(σn)‖,0≤n≤n*

(14)

作為停止準(zhǔn)則,其中n*=n*(δ,Uδ)為由(14)所確定的停止指標(biāo)。

算法1 EIT圖像重構(gòu)的非精確Newton迭代正則化方法Algorithm 1 Inexact Newton iterative regularization method for EIT image reconstruction

最后,以LDI迭代作為內(nèi)層迭代策略,引入基于LDI迭代的非精確Newton方法(Inexact Newton-LDI,INLDI),其格式為

(15)

以HPI迭代作為內(nèi)層迭代策略,引入基于HPI迭代的非精確Newton方法(Inexact Newton-HPI,INHPI),其格式為

(16)

以上求解EIT圖像重構(gòu)問題的非精確Newton迭代正則化方法的具體過程如算法1所示。

3 數(shù)值模擬

仿真模擬基于Matlab環(huán)境下Vauhkonen開發(fā)的EIDORS 2D軟件[22],建立16個(gè)電極的二維圓形EIT電極模型,采用傳統(tǒng)的相鄰電極電流激勵(lì)和相鄰電極電壓測(cè)量模式。設(shè)目標(biāo)區(qū)域Ω是坐標(biāo)系上規(guī)則的一個(gè)圓域,半徑為14,半徑無需標(biāo)注單位,僅作背景和目標(biāo)大小對(duì)比。為了避免反問題中的“Inverse Crime”,考慮了兩種不同規(guī)模的有限元網(wǎng)格剖分用于正反問題的求解,即用于正問題的網(wǎng)格由1 968個(gè)三角單元和1 049個(gè)結(jié)點(diǎn)構(gòu)成,如圖1(左一)所示,用于反問題的網(wǎng)格由492個(gè)三角單元和279個(gè)結(jié)點(diǎn)構(gòu)成,如圖1(左二)所示,其中網(wǎng)格周圍的綠色條紋表示理想電極的位置。重構(gòu)過程中,電極處的接觸阻抗假設(shè)已知,設(shè)置接觸阻抗大小為0.05。

圖1 有限元網(wǎng)格：細(xì)網(wǎng)格(左一)和粗網(wǎng)格(左二);三種測(cè)試模型：(a) 模型1;(b) 模型2;(c) 模型3Fig.1 FEM meshes：fine mesh (left one) and coarse mesh (left two),three test models：(a) model 1;(b) model 2;(3) model 3

為驗(yàn)證本文提出的INLDI迭代法和 INHPI迭代法在非線性EIT圖像重構(gòu)問題中的重構(gòu)性能,構(gòu)造了含有不同異物的三種測(cè)試模型作為研究對(duì)象,如圖1(a)～(c)所示,模型1包含1個(gè)圓形異物,模型2 包含2個(gè)方形異物,模型3包含5個(gè)不同異物,其中的標(biāo)尺顯示了各區(qū)域電導(dǎo)率σ的倒數(shù)(電阻率),藍(lán)色異物電阻率取值為1Ω·m,紅色背景電阻率取值為6Ω·m。為了便于對(duì)比分析,選取傳統(tǒng)的LDI迭代法、HPI迭代法與本文算法進(jìn)行比較。

真實(shí)測(cè)量數(shù)據(jù)由合成數(shù)據(jù)Usyn=F(σ*)(σ*為真實(shí)的電導(dǎo)率分布)計(jì)算得到,噪聲測(cè)量數(shù)據(jù)通過在合成數(shù)據(jù)上添加隨機(jī)噪聲產(chǎn)生,即Uδ=Usyn+δ·n,其中δ為噪聲水平,隨機(jī)變量n為服從[0,1]上的Gaussian分布。此外,為了比較EIT圖像重構(gòu)質(zhì)量,選擇較常用的相對(duì)誤差(RE)及相關(guān)系數(shù)(CC)作為衡量算法的客觀評(píng)價(jià)指標(biāo),其計(jì)算公式如下：

數(shù)值模擬的相關(guān)參數(shù)選取為τ=2.01,γ=0.8,初值σ0選取為背景值。由于過長(zhǎng)的計(jì)算時(shí)間并不符合EIT圖像重構(gòu)的需求,因此設(shè)置迭代的最大步數(shù)為nmax=10 000,kmax=50 000,即意味著即使未達(dá)到停止準(zhǔn)則,當(dāng)?shù)M(jìn)行到最大迭代步時(shí)停止迭代。

針對(duì)三種測(cè)試模型,為了測(cè)試四種方法的重構(gòu)效果,圖2描述了噪聲數(shù)據(jù)下相對(duì)誤差隨迭代時(shí)間的變化曲線,圖3描述了噪聲數(shù)據(jù)下相關(guān)系數(shù)隨迭代時(shí)間的變化曲線,用以分析算法的穩(wěn)定性和有效性。觀察圖2和圖3可以看到：四種方法可達(dá)到同等效果的RE值和CC值,說明具有相似的重構(gòu)精度和重構(gòu)質(zhì)量,但相比LDI和HPI迭代法,INLDI和INHPI這兩種迭代法較早地達(dá)到了穩(wěn)定狀態(tài),即RE曲線下降的更快、CC曲線上升的更快,這正體現(xiàn)了Newton迭代法具有收斂特性上的優(yōu)勢(shì),具有更快的收斂速度,而且也注意到,相比INLDI迭代法,INHPI迭代法具有更快的收斂速度,隨著噪聲水平的減少,INHPI迭代法在計(jì)算效率上的優(yōu)勢(shì)更加明顯,主要是內(nèi)層HPI迭代步引起的(具體見表1中的內(nèi)迭代步比較)。另一方面,從圖2和圖3中也可注意到：隨著噪聲水平的減少,四種方法獲得的RE值在變小而CC值在增大,說明重構(gòu)圖像的重構(gòu)精度和重構(gòu)質(zhì)量變得更好,但需要的計(jì)算時(shí)間卻大幅度增加了,這說明收斂速度與噪聲水平有關(guān)。

表1 LDI、HPI、INLDI以及INHPI迭代法重構(gòu)模型2的數(shù)值比較Table 1 Comparison by LDI,HPI,INLDI and INHPI for model 2

圖2 相對(duì)誤差曲線：(a) 模型1;(b)模型2;(c)模型3Fig.2 Relative error curves：(a) model 1;(b) model 2;(3) model 3

圖3 相關(guān)系數(shù)變化曲線：(a) 模型1;(b)模型;(c)模型3Fig.3 Correlation coefficient curves：(a) model 1;(b) model 2;(3) model 3

針對(duì)圖2和圖3中的結(jié)果,這里僅以模型2為例(見圖2(b)和圖3(b)),表1記錄了在不同噪聲水平下,四種方法的詳細(xì)數(shù)值結(jié)果比較,其中n*表示LDI迭代法與HPI迭代法的迭代終止步,nouter表示INHPI方法與INLDI方法的外層Newton迭代步,ninner表示INHPI方法與INLDI方法的內(nèi)層總迭代步,Time(s)表示計(jì)算所需時(shí)間。從表1可以看出,四種方法具有相似的重構(gòu)精度,但重構(gòu)效率明顯不同。相比LDI和HPI迭代法,INHPI和INLDI迭代法的重構(gòu)時(shí)間顯著減少,大大提高了重構(gòu)效率。另一方面,比較LDI和HPI迭代法,驗(yàn)證了HPI迭代法只需LDI大約一半的迭代步和計(jì)算時(shí)間;但比較INLDI和INHPI迭代法,INHPI迭代法總的內(nèi)層迭代次數(shù)比INLDI迭代法減少了約一半,高噪聲水平下的重構(gòu)效率雖有所提高,但并未像預(yù)期的那樣提高約一倍,主要原因在于需要調(diào)用正問題的外層迭代為收斂速度較快的Newton步(外層每步迭代耗時(shí)多,但需要的外層迭代步數(shù)少),內(nèi)層每步迭代耗時(shí)較少,迭代步數(shù)雖節(jié)約了一半左右,但在內(nèi)層迭代迭代步數(shù)基數(shù)較小的情況下,計(jì)算效率無法凸顯出來,但在低噪聲水平下,隨著需要的內(nèi)層迭代步數(shù)的劇增,相比INLDI迭代法,INHPI迭代法在計(jì)算效率上的優(yōu)勢(shì)就呈現(xiàn)出來了,改進(jìn)了成像效率。

最后,為了可視化重構(gòu)結(jié)果,鑒于四種方法具有相似的重構(gòu)精度和重構(gòu)質(zhì)量,針對(duì)三種測(cè)試模型,這里僅給出不同噪聲水平下INHPI迭代法的重構(gòu)圖像,如圖4所示。從圖4可以看出,在低噪聲水平下,重構(gòu)圖像均能夠較準(zhǔn)確地反映目標(biāo)體的大小、位置和性質(zhì),隨著噪聲水平的增加,重構(gòu)結(jié)果比較穩(wěn)定,說明對(duì)噪聲較為魯棒,不過隨著噪聲水平的增加,重構(gòu)圖像的分辨率有所降低。

圖4 基于不同噪聲數(shù)據(jù)的INHPI迭代法所得重構(gòu)圖像,從左到右依次為：真實(shí)模型,以及δ=5%,δ=2%,δ=0.5%,δ=0.05%噪聲數(shù)據(jù)下對(duì)應(yīng)的重構(gòu)圖像Fig.4 Reconstructed images for three test models by INHPI method under different noise levels.Reconstructions from left to right：true model,reconstructed images under δ=5%,δ=2%,δ=0.5%,δ=0.05%,respectively

4 結(jié) 論

基于LDI迭代法和HPI迭代法,提出了兩種非精確Newton迭代算法：INLDI迭代法和 INHPI迭代法,用于處理非線性EIT圖像重構(gòu)問題。其基本思想是應(yīng)用LDI迭代法或HPI迭代法逐步求解局部線性化的EIT觀測(cè)響應(yīng)模型來構(gòu)造內(nèi)層迭代序列,并結(jié)合適當(dāng)?shù)耐Ｖ箿?zhǔn)則終止內(nèi)層迭代,進(jìn)而將得到的近似重構(gòu)解作為新的外層迭代點(diǎn)并以偏差原則作為外層停止準(zhǔn)則。對(duì)LDI迭代法和HPI迭代法進(jìn)行了比較,并以重構(gòu)圖像的相對(duì)誤差和相關(guān)系數(shù)作為評(píng)價(jià)標(biāo)準(zhǔn)。結(jié)果表明,本文算法具有節(jié)約迭代步以及提高計(jì)算效率方面的優(yōu)越性,具有更快的收斂速度,可有效改進(jìn)成像效率。