摘要本文提出了一种新的电阻抗成像(EIT)图像重建算法—基于修正B_k的拟Newton类算法，并通过数值试验验证了该算法的有效性。新算法基于目标函数梯度向量的快速算法，利用基于修正B_k的拟Newton法来求解EIT逆问题。与著名的Newton-Raphson算法相比，该算法的突出优点是计算量小，数值稳定性好，具有整体收敛性。

A Quasi-Newton Algorithm for Electrical Impedance Tomography

—a Switch Algorithm Based on Modified B_k

Li Wenbo

(Institute of Space Medical Engineering, Beijing 100094)

Wang Boliang

(Department of Electronic Technology, National University of Defense Technology, ChangSha 410073)

Sun Guiju

(Institute of Space Medical Engineering, Beijing 100094)

Xiang Fei

(Institute of Space Medical Engineering, Beijing 100094)

Liu Xishun

(Department of Electronic Technology, National University of Defense Technology, ChangSha 410073)

Abstract

In this paper we propose a new image reconstruction algorithm for electrical impedance tomography-Quasi-Newton algorithm, based on modified B_k. Numerical tests were made to verify the effectiveness of this algorithm. Based on the fast algorithm for calculating gradient vectors of object function, this new algorithm makes use of Quasi-Newton algorithm based on modified B_k to solve the inverse problem of EIT. Comparing this new algorithm with the famous Newton-Raphson method, the outstandng advantages of this new algorithm are that the workload of computation can be reduced, and numerical stability and integral convergence can be obtained.

Key words:Electrical impedance tomography; Reconstruction algorithm; Quasi-Newten algorithm; Invere problem

0引言

电阻抗成像(electrical impedance tomography, EIT)作为一种功能性医学成像技术的研究始于70年代末。由于它具有设备简单、成本低廉、对人体无害、不要求特殊的工作环境等优点，使其成为近年来生物医学工程界的研究热点之一。

EIT的数学模型由下述椭圆形方程边值问题描述^［1］：

(1)

,(2)

(3)

其中Ω为物体所在的空间区域，为其边界，U为求解区域的电势分布，V、J分别为边界电压和边界电流密度分布，ρ为待求电阻率分布。与一般的偏微分方程的定解问题不同，这里的系数ρ和电势分布U都是未知的，而且我们要求解的正是物体的电阻率分布ρ，因此这是一个二阶椭圆形方程未知系数辨识的逆问题(inverse problem)。

对于区域Ω，在边界电流密度J给定的条件下，记电阻率分布到边界电压分布之间的映射为F(ρ)=V，其中ρ、V分别属于一适当的函数空间(离散化后为两个有限维空间)。EIT算法包括两部分：求解正问题和逆问题。在ρ设定的情况下，求F(ρ)的过程称为EIT的正问题。考虑到EIT问题的特点，一般用有限元方法求解EIT正问题，此时求解区域离散成许多小单元(三角形或四边形)，ρ以某种插值函数的形式表示。记F的离散形式为f，则f为两个有限维空间之间的映射：R^m→Rⁿ,f(ρ)=v。在ρ未知的情况下，利用已知的边界条件求ρ的问题称为EIT的逆问题。EIT的逆问题通常要用最优化方法求解。

实践中，为了获取尽可能多的关于电阻率分布ρ的信息，一般采用多次施加电流的方法，进行多次边界电压测量。设共有ρ次线性独立的电流模式注入，相应的正问题的解及边界电压测量值分别记为f_i(ρ)和v_i(i=1,2,…,p)，则EIT逆问题离散化后成为下面的非线性最小二乘问题：

minφ(ρ)，(4)

其中

(5)

对于EIT逆问题，现在已有许多重建算法，其中理论上较完善、重建效果较好的算法是Newton-Raphson算法及其改进形式^［2］，本质上是求解上述非线性最小二乘问题(4)、(5)的Gauss-Newton法和LM法。但是，Newton-Raphson算法的缺点也是明显的，主要是：计算量大，而且随空间分辨率的提高，计算量迅速增长；数值稳定性不够好；只有局部收敛性。针对Newton-Raphson算法的这些缺点，文献［3］导出了一种目标函数梯度向量的快速算法，在此基础上提出了一种基于修正H_k的拟Newton法(组合变尺度法)用于求解EIT逆问题。与Newton-Raphson算法相比，拟Newton法的性能明显要好：拟Newton法的计算量随空间分辨率提高而增长的速度比Newton-Raphson算法要慢得多；而且该算法具有整体收敛性，因此迭代初值的选取比较方便；拟Newton法的数值稳定性也明显好于Newton-Raphson算法，因此拟Newton法有更强的抗扰动能力。由于上述优点，使拟Newton类算法在求解EIT逆问题方面成为一类很有竞争力的算法。

1用基于修正B_k的拟Newton法求解EIT逆问题

1.1数值稳定性分析

我们知道，拟Newton法是求解一般无约束最优化问题的一种重要方法。它的基本思想是：利用初始给定的对称正定矩阵B_k和H_k作为对目标函数的Hesse矩阵或Hesse矩阵的逆矩阵的估计，选定迭代初值，利用目标函数的梯度向量和B_k或H_k矩阵计算出下降方向，进行迭代。在迭代过程中不断修正B_k和H_k，使算法在不用Hesse矩阵的情况下具有Newton法的某些优良性质，因此它是一类非常有效的算法。

常用的修正B_k的公式有两个，分别称为关于B_k的DFP公式和BFGS公式；常用的修正H_k的公式也有两个，分别称为关于H_k的DFP公式和BFGS公式^［4］。实际计算时，DFP公式和BFGS公式可单独使用，也可以结合起来构成所谓的开关算法(又称组合变尺度法)。一般认为，组合变尺度法的数值稳定性较普通的拟Newton算法更好。

文献［3］导出了求解EIT逆问题的基于修正H_k的拟Newton法，同样我们也可以用基于修正B_k的拟Newton法来求解EIT问题。而且，文献［5］中指出，基于修正B_k的拟Newton法比基于修正H_k的拟Newton法有更高的数值稳定性。为了以后讨论方便，我们先给出修正B_k的DFP公式，

(6)

和BFGS公式，

(7)

以及修正H_k的DFP公式，

(8)

和BFGS公式，

(9)

其中s_k=ρ^k+1-ρ^k，y_k=g_k+1-g_k；ρ^k为第k次迭代时的电阻率分布向量，g_k为第k次迭代时的梯度向量。比较上面四个公式可以看出，修正B_k需要计算y_ky_Tk/s_Tky_k；修正H_k需要计算s_ss_Tk/s_Tky_k。假定S_k、y_k均有计算误差ε，则计算y_ky_Tk/s_Tky_k和计算s_ks_Tk／S^T_ky_k的误差分别为O(［y_k²₂/s_Tky_k］·ε)和O(［‖S_k‖²₂/s_T