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

陳 鋒，王化祥

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

電學層析成像(Electrical Tomography，ET)技術(shù)[1]是近幾年快速發(fā)展起來的一種檢測技術(shù)，相比于其他的CT技術(shù)而言，其具有非侵入性、便攜性、無輻射以及價格低廉等優(yōu)勢，被廣泛應(yīng)用于工業(yè)和醫(yī)學等領(lǐng)域，例如工業(yè)過程中的多相流檢測以及人體病變組織監(jiān)護等。而電阻層析成像(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ù)上位機的要求選取激勵和測量電極，測取電壓/電容值，進行初步處理后傳輸?shù)胶蠖恕？杉毞譃榻o出激勵信號的信號源、負責采集敏感場中的電信號的數(shù)據(jù)采集以及控制選通電極部分以及負責傳輸信號的傳輸接口。圖像重建與分析模塊是通過編寫上位機的軟件系統(tǒng)來實現(xiàn)的，主要負責控制邊界數(shù)據(jù)來完成數(shù)據(jù)采集模塊的測量過程，并通過編寫設(shè)計好的通訊模塊(通常采用是USB通訊)來接收數(shù)據(jù)采集模塊傳輸過來的數(shù)據(jù)，最后根據(jù)事先編寫好的成像算法來完成圖像重建工作。

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

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

1 GMRES算法

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

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

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

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

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

(2)接著計算:

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

(3)近似解的形成

找到向量ym并將其進行極小化:

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

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

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

2 線性逼近

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

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

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

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

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

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

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

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

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

以上兩式中的 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算法需要消耗的時間為2k2n+2kNNZ，所占用的空間總數(shù)為(k+3)n+k2/2[11]。很顯然，完全GMRES算法隨著迭代步驟k的增加，所耗費的代價也更高。而Restarted GMRES算法每m步迭代后將會重頭開始迭代，這樣可以很有效地降低計算速度和對所需要占用空間的需求，而這里的m在大多數(shù)情況下遠小于總的迭代次數(shù)n。

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

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

(2)應(yīng)用Amoldi方法分別計算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 實驗仿真與結(jié)果分析

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

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

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

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

分別采用 Full-GMRES，GMRES(m)和 DEFLGMRES(m，l)算法，對圖像進行重建，如圖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)算法在到達各自最低成像誤差時所需要的計算時間。結(jié)果顯示，DEFLGMRES(m，l)算法成像誤差小，且收斂速度快。

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

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

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

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

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

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

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

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

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

[8]熊新平.并行預(yù)條件GMRES方法[J].計算機工程與設(shè)計，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.