基于Deflation技術的預調制Restarted GMRES算法*

陳 鋒，王化祥

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

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

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

圖1 ET系統示意圖

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

1 GMRES算法

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

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

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

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

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

(2)接著計算:

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

(3)近似解的形成

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

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

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

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

2 線性逼近

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

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

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

3 基于Deflation技術的預調制Restarted GMRES算法

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

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

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

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

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

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

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

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

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

其中

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

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

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

(2)應用Amoldi方法分別計算Vm、β以及;

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

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

4 預調制器的結構

JTJ的最小特征值會降低Restarted GMRES算法的收斂速度，因此使用預調節器來移除JTJ的最小特征值，從而提高Restarted GMRES算法的收斂速度。

將JTJ定義為A，假定P是r最小特征根所對應的不變子空間，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值的多重根。

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

式中T=UTAU是將A限制在P空間內是將A限制在W空間內。

重寫

并設:

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

因此可得:

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

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

5 實驗仿真與結果分析

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

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

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

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

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

圖3 圖像重建

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

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

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

[1]王超，王化祥.醫學電阻抗成像系統電極結構優化設計[J].第四軍醫大學學報，2001，22(1):78－79.

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

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

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

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

[6]梅丹.基于解空間分解的GMRES算法及其在圖像處理中的應用[J].計算機與數字工程，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]熊新平.并行預條件GMRES方法[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.