基于非精確Newton迭代正則化的EIT圖像重構算法

趙一帆,王 靜

(黑龍江大學 數學科學學院,哈爾濱 150080)

0 引 言

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

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

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

1 EIT數學模型

目前廣泛采用的EIT數學模型為比較貼近實際的全電極數學模型[20]：

(1)

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

A(σ)b=f

(2)

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

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

Uδ=F(σ)+e

(3)

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

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

2.1 LDI迭代法和HPI迭代法

EIT圖像重構問題本身存在嚴重的不適定性,最小二乘法是處理該問題最為流行的方法,即極小化數據擬合項

(4)

最為經典的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為適當選取的步長因子,這里采用比較適宜在EIT領域使用的變步長因子選取準則,滿足

(7)

2.2 INLDI迭代法和INHPI迭代法

鑒于非精確Newton迭代正則化方法具有較快的收斂速度,考慮將其用于處理EIT圖像重構問題。該方法包括內外兩層迭代：外層迭代為非精確Newton方法且以偏差原則作為停止準則,內層迭代這里考慮采用LDI迭代或HPI迭代并結合選取適當的步長來加速。

首先構造內層迭代。假設已知當前迭代步為σn,通過極小化問題

(8)

構造相應的迭代格式

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

(9)

(10)

若選取

(11)

則為LDI迭代格式;若選取

(12)

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

(13)

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

(14)

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

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

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

(15)

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

(16)

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

3 數值模擬

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

圖1 有限元網格：細網格(左一)和粗網格(左二);三種測試模型：(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

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

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

數值模擬的相關參數選取為τ=2.01,γ=0.8,初值σ0選取為背景值。由于過長的計算時間并不符合EIT圖像重構的需求,因此設置迭代的最大步數為nmax=10 000,kmax=50 000,即意味著即使未達到停止準則,當迭代進行到最大迭代步時停止迭代。

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

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

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

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

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

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

圖4 基于不同噪聲數據的INHPI迭代法所得重構圖像,從左到右依次為：真實模型,以及δ=5%,δ=2%,δ=0.5%,δ=0.05%噪聲數據下對應的重構圖像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 結 論

基于LDI迭代法和HPI迭代法,提出了兩種非精確Newton迭代算法：INLDI迭代法和 INHPI迭代法,用于處理非線性EIT圖像重構問題。其基本思想是應用LDI迭代法或HPI迭代法逐步求解局部線性化的EIT觀測響應模型來構造內層迭代序列,并結合適當的停止準則終止內層迭代,進而將得到的近似重構解作為新的外層迭代點并以偏差原則作為外層停止準則。對LDI迭代法和HPI迭代法進行了比較,并以重構圖像的相對誤差和相關系數作為評價標準。結果表明,本文算法具有節約迭代步以及提高計算效率方面的優越性,具有更快的收斂速度,可有效改進成像效率。