基于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.