基于Deflation技術(shù)的預(yù)調(diào)制Restarted GMRES算法*

陳 鋒，王化祥

(天津大學(xué)電氣與自動化工程學(xué)院，天津300072)

電學(xué)層析成像(Electrical Tomography，ET)技術(shù)[1]是近幾年快速發(fā)展起來的一種檢測技術(shù)，相比于其他的CT技術(shù)而言，其具有非侵入性、便攜性、無輻射以及價(jià)格低廉等優(yōu)勢，被廣泛應(yīng)用于工業(yè)和醫(yī)學(xué)等領(lǐng)域，例如工業(yè)過程中的多相流檢測以及人體病變組織監(jiān)護(hù)等。而電阻層析成像(Electrical Resistance Tomography，ERT)、電容層析成像(Electrical Capacitance Tomography，ECT)以及電磁層析成像(Electrical Magnetic Tomography，EMT)[2－4]是三種比較常見的ET技術(shù)。

從模塊結(jié)構(gòu)上看，ET系統(tǒng)主要包括傳感器模塊、數(shù)據(jù)采集模塊和圖像重建模塊所構(gòu)成[5]，對被測對象給出激勵信號(通常是電流/電壓信號)后，通過前端傳感器模塊以及硬件電路對測得的電信號進(jìn)行采集，同時(shí)根據(jù)上位機(jī)的要求選取激勵和測量電極，測取電壓/電容值，進(jìn)行初步處理后傳輸?shù)胶蠖恕？杉?xì)分為給出激勵信號的信號源、負(fù)責(zé)采集敏感場中的電信號的數(shù)據(jù)采集以及控制選通電極部分以及負(fù)責(zé)傳輸信號的傳輸接口。圖像重建與分析模塊是通過編寫上位機(jī)的軟件系統(tǒng)來實(shí)現(xiàn)的，主要負(fù)責(zé)控制邊界數(shù)據(jù)來完成數(shù)據(jù)采集模塊的測量過程，并通過編寫設(shè)計(jì)好的通訊模塊(通常采用是USB通訊)來接收數(shù)據(jù)采集模塊傳輸過來的數(shù)據(jù)，最后根據(jù)事先編寫好的成像算法來完成圖像重建工作。

圖1 ET系統(tǒng)示意圖

電學(xué)層析成像的圖像重建需要對逆問題進(jìn)行求解，而求解過程中存在著非線性、欠定性以及病態(tài)性嚴(yán)重等難題，使得圖像重建可能不收斂，或者即使收斂，但獲得的圖像分辨率較低。本文所研究的電學(xué)層析成像過程中涉及到非對稱矩陣的線性方程組的求解過程，而GMRES算法通常被廣泛應(yīng)用于求解線性方程和最小二乘問題，是Krylov子空間迭代法中的一種方法，可以用于求解Krylov子空間中的近似解，同時(shí)，預(yù)條件技術(shù)使其成為求解大型非對稱線性方程組的一種有效方法[6]。

1 GMRES算法

廣義最小殘差GMRES(Generalized Minimal Residual)算法是由Martin Schultz和Yousef Saad提出的一種迭代算法[7]，主要應(yīng)用于一些系數(shù)矩陣非對稱的大型稀疏線性方程組。可以被應(yīng)用于電學(xué)層析成像過程中涉及到非對稱矩陣的線性方程組的求解過程。

通常，在電阻抗層析成像過程中，需要求解以下線性方程組:

而靈敏度矩陣S通常不是方陣，以此無法用g=S－1z來求灰度值向量g。此時(shí)將預(yù)調(diào)制器矩陣M引入進(jìn)來，先假定將M作用于左邊，即將M－1左乘到方程兩邊，則所要求解的EIT方程組變?yōu)?

因此，GMRES在預(yù)調(diào)制下的n步迭代可以描述如下:Preconditioned GMRES(m)[8]

(1)首先，先確定一個初始值g0以及m的值;

(2)接著計(jì)算:

其中，j=1，2，…，m，這里將 Hm定義為(m+1)×m 的上Heissenberg陣，其非零元素為系數(shù)hij

(3)近似解的形成

找到向量ym并將其進(jìn)行極小化:

式中:e1=[1，0，…，0]T∈Rm+1，y∈Rm

式中:Vm是由正交的基向量{v1，v2，…，vm}構(gòu)成的一個n×m矩陣。

由于Vm完全正交，隨著Krylov子空間增大，傳統(tǒng)的完全GMRES算法所需的計(jì)算時(shí)間明顯變長以及占用的空間也有著顯著增加。而Restarted GMRES算法可以將Krylov子空間的規(guī)模限定在一個比完全GMRES算法中所使用的規(guī)模小得多的一個固定值上。但是，Restarted GMRES算法的收斂速率不僅僅取決于特征值的分布，同時(shí)也取決于Krylov子空間的規(guī)模，其收斂速度比完全GMRES算法慢，并且重開啟過程會丟失一些信息，諸如最小的Ritz值[9]。因此，為了加速收斂速度，可將最小特征值對應(yīng)的近似特征向量添加到Krylov子空間中，這樣可以非常有效地估算出特征譜中的最小特征向量。這里也可以將Amoldi向量保存下來以便下一次循環(huán)使用，這樣可以有效地消除Restarted GMRES算法中病態(tài)性的影響[10]。

2 線性逼近

假設(shè)在EIT系統(tǒng)中，電導(dǎo)率的分布上的變化很小，因此可以通過線性化方法來求解逆問題，通過使用離散有限元法，忽略高階項(xiàng)，將EIT近似成如下形式:

其中，δσ∈Rn×1是電導(dǎo)率的變化(n 是網(wǎng)格數(shù))，δU∈Rm×1是邊界電壓的擾動(m 測量的電壓值)，J∈Rm×n是Jacobi矩陣，通過Geselowitz靈敏度定理，可以得到Jacobi矩陣如下:

其中，u(Id)為第d個驅(qū)動模型的電勢，u(Im)為第m個驅(qū)動模型的假定電壓，則圖像重建過程就等價(jià)于求解最小二乘問題:

3 基于Deflation技術(shù)的預(yù)調(diào)制Restarted GMRES算法

基于Deflation技術(shù)的預(yù)調(diào)制Restarted GMRES算法(Restarted GMRES Preconditioned by Deflation Technique Algorithm)，是GMRES的一種改進(jìn)，可以很好地提高速度和精度，下文描述為DEFLGMRES(m，l)。

考慮到矩陣J不是一個方陣，GMRES算法可以用于求解式(14):

式(15)中給定一個初值x0，采用完全GMRES算法從仿射子空間x0+Kk尋求一個近似解。其中:

通過施加的Galerkin條件可以得到:

采用 Amoldi方法，設(shè) v1=r0/‖r0‖，β=‖r0‖，當(dāng)進(jìn)行到第k步時(shí)，可以得到:

以上兩式中的 Vk=[v1，v2，…，vk]是基于 Krylov子空間 Kk的正交集，是 Heissenberg，其非零元素Hij是在Amoldi方法中定義的，e1=(1，0，0，…，0)T∈Rn。

因此，可以得到近似解如下:

通過式(19)和式(20)可以得到余差范數(shù)如下:

根據(jù)式(21)和式(22)，完全GMRES算法可以描述如下:

其中

用NNZ來表示JTJ中的非零元素個數(shù)，則k步迭代的完全GMRES算法需要消耗的時(shí)間為2k2n+2kNNZ，所占用的空間總數(shù)為(k+3)n+k2/2[11]。很顯然，完全GMRES算法隨著迭代步驟k的增加，所耗費(fèi)的代價(jià)也更高。而Restarted GMRES算法每m步迭代后將會重頭開始迭代，這樣可以很有效地降低計(jì)算速度和對所需要占用空間的需求，而這里的m在大多數(shù)情況下遠(yuǎn)小于總的迭代次數(shù)n。

以下是Restarted GMRES算法(GMRES(m))的步驟:

(1)設(shè)置容差值ε和x0;

(2)應(yīng)用Amoldi方法分別計(jì)算Vm、β以及;

(3)通過求解式(23)和式(24)來得到Xm;

(5)設(shè)x0=xm，則返回第2步。

4 預(yù)調(diào)制器的結(jié)構(gòu)

JTJ的最小特征值會降低Restarted GMRES算法的收斂速度，因此使用預(yù)調(diào)節(jié)器來移除JTJ的最小特征值，從而提高Restarted GMRES算法的收斂速度。

將JTJ定義為A，假定P是r最小特征根所對應(yīng)的不變子空間，U是P的一個正交基。

引理 1[12]如果 T=UTAU，M=In+U(1/|λn|TIr)UT，則 M 是非奇異的，且 M－1=In+U(|λn|T－1－Ir)UT。AM－1的特征值為 λr+1，λr+2，…，λn，|λn|，|λn|為最后一個r值的多重根。

證明設(shè)Z=[U，W]為Rn的一個正交基，設(shè):

式中T=UTAU是將A限制在P空間內(nèi)是將A限制在W空間內(nèi)。

重寫

并設(shè):

其中，In－r為單位矩陣，T是非奇異的，是可逆的，表示如下:

因此可得:

因此可得 AM－1的特征值為。

引理1證明了預(yù)控制器M可將最小的特征值λ1，λ2，…，λr用|λn|替換。

5 實(shí)驗(yàn)仿真與結(jié)果分析

本論文采用了兩種模型，首先將仿真中背景設(shè)置為礦化水(σwater=1 ms/cm)，而場內(nèi)的物體設(shè)置為油(σoil=0.002 ms/cm)，采用16電極激勵形式。同時(shí)，為了和實(shí)際測量系統(tǒng)的噪聲水平一致，將±1%的高斯隨機(jī)噪聲加入到仿真電壓中。

對 Full-GMRES，GMRES(m)和 DEFLGMRES(m，l)進(jìn)行定量分析，EIT中的相對成像誤差如下:

設(shè)置參數(shù)迭代次數(shù)m=10，l=3，圖2在模型情況下仿真，結(jié)果顯示，DEFLGMRES(m，l)算法在提高收斂速度的同時(shí)成像誤差低，穩(wěn)定性好。

圖 2Full-GMRES，GMRES(m)和DEFLGMRES(m，l)相對成像誤差

分別采用 Full-GMRES，GMRES(m)和 DEFLGMRES(m，l)算法，對圖像進(jìn)行重建，如圖3所示。結(jié)果顯示，開始Full-GMRES算法獲得的圖像最佳，但是隨著迭代次數(shù)的增加，圖像質(zhì)量下降;而GMRES(m)和DEFLGMRES(m，l)算法隨著迭代次數(shù)的增加，所成圖像質(zhì)量有明顯提高，更重要的是，隨著迭代次數(shù)的增加，DEFLGMRES(m，l)算法的成像質(zhì)量優(yōu)于GMRES(m)算法。

圖3 圖像重建

表1分別顯示了 Full-GMRES，GMRES(m)和DEFLGMRES(m，l)算法在到達(dá)各自最低成像誤差時(shí)所需要的計(jì)算時(shí)間。結(jié)果顯示，DEFLGMRES(m，l)算法成像誤差小，且收斂速度快。

表1 Full-GMRES，GMRES(m)和 DEFLGMRES(m，l)算法計(jì)算時(shí)間比較

由上述的建模和仿真以及實(shí)驗(yàn)數(shù)據(jù)可以發(fā)現(xiàn)，DEFLGMRES算法在成像圖像重建上，提高了計(jì)算速度和成像分辨率，同時(shí)也節(jié)省了硬件資源的空間，這意味著它適合于在線操作，例如多相流檢測，而且今后可以將此算法由2D成像領(lǐng)域推廣到3D范圍，這樣，其速度和占用空間小的優(yōu)勢可以得到更好地發(fā)揮。

[1]王超，王化祥.醫(yī)學(xué)電阻抗成像系統(tǒng)電極結(jié)構(gòu)優(yōu)化設(shè)計(jì)[J].第四軍醫(yī)大學(xué)學(xué)報(bào)，2001，22(1):78－79.

[2]王研.電阻抗成像電極系統(tǒng)優(yōu)化設(shè)計(jì)仿真研究[J].中華醫(yī)學(xué)研究雜志，2005，5(8):759－761.

[3]董峰，鄧湘，徐立軍，等.過程層析成像技術(shù)綜述[J].儀器儀表用戶，2001，8(1):69－7.

[4]唐磊.電學(xué)成像系統(tǒng)圖像重建算法研究及可視化軟件設(shè)計(jì)[D].天津:天津大學(xué)，2008.

[5]何泳成.電學(xué)層析成像重建算法研究及軟件系統(tǒng)設(shè)計(jì)[D].天津:天津大學(xué)，2010.

[6]梅丹.基于解空間分解的GMRES算法及其在圖像處理中的應(yīng)用[J].計(jì)算機(jī)與數(shù)字工程，2009，37(12)，139.

[7]黃云清，VAN der Vorst V A.Some Observations on the Convergence Behavior of GMRES[J].湘潭大學(xué)學(xué)報(bào)，1989，1(5):16－22.

[8]熊新平.并行預(yù)條件GMRES方法[J].計(jì)算機(jī)工程與設(shè)計(jì)，1994，4(3):25－32.

[9]J Erhel K B，Pohl B.Restarted GMRES Preconditioned by Deflation[J].Journal of Computational and Applied Mathematics，1996，69:303－318.

[10]Morgan R B.A Restarted GMRES Method Augmented with Eigenvectors，SIAM J MATRIX ANAL APPL，1995，6.1154－71.

[11]Saad Y.Iterative Methods for Sparse Linear Systems[M].Society for Industrial and Applied Mathematics，2003.

[12]Erhel J K B，B Pohl.Restarted GMRES Preconditioned by Deflation[J].Journal of Computational and Applied Mathematics，1996，69:303－18.