脑组织牵拉变形矫正的XFEM 建模及其验证

王巍伟 李 萍 章琛曦 宋志坚*1(复旦大学数字医学研究中心，上海 200032)

2(上海市医学图像处理与计算机辅助手术重点实验室，上海 200032)

引言

以CT、MRI 等医学图像为基础的神经外科手术导航系统(image guided navigation system，IGNS)，对于辅助医生准确识别与切除脑部病变组织、避免损伤重要神经或功能区等具有重要的临床意义［1］。但是，在神经外科手术中，由于开颅后颅内压的改变、脑脊液的流出、牵拉、肿瘤切除等因素的影响，脑组织会产生较为严重的变形。这导致由术前数据建立的图像引导系统与实际术中脑组织情况产生较大偏差，大幅降低了IGNS 的定位精度。脑组织变形的产生可分为3 种类型:第1 种是硬脑膜打开后产生，第2 种是手术牵拉产生，第3 种是病灶切除产生。牵拉变形是外科医生为寻找到达大脑皮层下目标(如肿瘤)的路径，用脑压板拨开脑组织导致的脑组织变形［2］。目前，国内外的研究主要集中在第1 种脑组织变形上，较少考虑牵拉和肿瘤切除引起的第2、3 种脑组织变形，而这两种脑组织变形正是降低IGNS 术中导航精度的主要因素。矫正手术牵拉和切除引起脑组织变形，主要在于都会导致脑组织拓扑结构改变，且矫正牵拉变形是矫正切除变形的基础，因此矫正脑组织牵拉变形显得尤为重要。

目前主要有两类矫正脑组织变形的方法:术中影像和生物力学模型矫正。虽然术中影像如术中CT(intraoperative CT，iCT)和术中MR(intraoperative MR，iMR)能够较准确地矫正脑组织变形，但其缺点是成本高、操作不方便、容易引起感染且延长了手术时间，因此限制了术中影像矫正的临床应用。生物力学模型的矫正方法通过建立生物力学模型，然后利用有限元方法(finite element method，FEM)的求解来矫正脑组织变形，从而克服了术中影像方法的缺点。Ferrant 使用线弹性生物力学模型建模，通过删除牵拉路径上所有网格，实现模型的更新［3］。Miga 应用固结理论，建立线性多孔弹性生物力学模型，沿牵拉路径切割网格，产生两倍的网格节点，实现模型的更新［4］。动物实验表明，该模型可矫正平均66%的脑组织位移误差，降低80%的定位误差［5］。Sun 同样使用线性多孔弹性生物力学模型，通过FEM 计算，模拟牵拉脑组织引起的脑组织变形［6］。该方法利用立体摄像机获取脑压板的运动位移场，并以此为边界条件对模型进行求解。此方法获取大约75%的脑皮层运动。

以上研究者都利用FEM 来处理脑组织拓扑结构的改变，概括起来有两种处理方式:完全重新网格化或者将拓扑结构改变区域局部网格化［7］。不管FEM 采用哪种处理方法，都需要改变网格原有的拓扑结构，会破坏原有网格的连续性，而扩展有限元(extended finite element method，XFEM)［8］则克服了这个缺点。它处理脑组织变形首先不考虑拓扑结构变化产生的不连续性，而是直接将脑组织网格化，再根据裂缝的具体形状，沿着裂缝经过的网格节点加上新的虚自由度。这样处理脑组织牵拉产生的裂缝，能在不破坏原有网格连续性的基础上准确地表达裂缝产生的拓扑不连续性。Vigneron 首次提出将XFEM 应用于线弹性生物力学模型中，并利用基于iMR 图像的非刚性配准技术，获取牵拉前后脑组织的位移，以此作为边界条件求解模型［2，9-10］。此方法使用Canny 算子对配准前后的图像进行边缘检测，结果表明配准后图像的Hausdorff 距离减小。但该方法仅得到定性结果，没有给出XFEM 在矫正牵拉时引起脑组织变形有效性的定量性结果。

本研究提出，将XFEM 应用于线弹性生物力学模型，用来矫正牵拉变形的系统框架。选取脑组织体模作为试验对象，用iCT 图像作为评估的金标准;利用植入的不锈钢标记物，定量分析了XFEM 在矫正脑组织牵拉变形的有效性。同时，归纳比较不同尺寸网格对于XFEM 矫正牵拉变形的影响。

1 牵拉变形矫正系统框架

基于XFEM 矫正脑组织牵拉变形的系统框架如图1 所示。步骤1 和步骤2 为打开硬脑膜后的脑变形状态(以下称为术前状态)，得到硬膜打开后已矫正脑组织变形的脑组织网格;步骤3 通过iCT/iMR 扫描得到手术中牵拉产生的脑组织变形图像;步骤4 将牵拉变形图像通过点配准与术前图像配准，以确定在术前图像坐标系中脑压板初始位置和脑压板牵拉位移;步骤5 将获得的脑压板初始位置和脑压板牵拉位移引入到基于术前图像的XFEM线弹性模型中，作为边界条件;步骤6 根据步骤5 求得的边界条件和已知的零位移边界条件求解XFEM线弹性模型方程，据此，得到牵拉变形网格;步骤7更新脑组织术前图像，得到矫正脑组织牵拉变形后的图像;步骤8 将更新后的图像显示给用户。

1.1 XFEM 原理

FEM 无法很好地处理拓扑结构改变的不连续情况，这是由于决定网格节点生物力学属性的形函数(nodal shape function，NSF)是连续函数。XFEM扩展了FEM 的功能，通过增加裂缝处网格节点的虚自由度来更好地模拟裂缝等不连续情况，克服了FEM 处理裂缝时需要重新进行网格化的缺点。XFEM 可表达为

式中，uXFEM为XFEM 位移场。

式(1)等号右侧第1 部分是FEM 网格节点位移场，I 表示FEM 网格节点集合，φi是FEM 的NSF，ui是网格节点FEM 自由度。与FEM 不同，XFEM 在网格点集J 和K 两个集合上新增加了虚自由度，用这些虚自由度来表述裂缝的不连续性。点集J 和K的关系可表达为

图1 应用XFEM 建模矫正牵拉变形的系统框架Fig．1 XFEM framework for retraction correction

如图2 空心圆点所示，当裂缝将网格单元的所有连线都完全分割时，此类受影响网格节点组成点集J 。式(1)等号右侧第2 部分新增加Heaviside 函数来表达裂缝对点集J 产生的不连续性。单位阶跃函数H 用公式表示为

式中，x 为网格节点，x*为裂缝上与距离x 点最近点，en为x*处裂缝法向量，aj表示点集J 的虚自由度。

如图2 实心方点所示，当裂缝末端进入某网格单元但没有将该单元节点所有连线完全分割时，此类受影响网格节点组成受裂尖影响的点集K 。

式(1)等号右侧第3 部分表达了点集K 的不连续性，包括半径和角度的计算，Clk 为点集K 的新增虚自由度，Fl为点集K 的新增不连续函数，可表示为

式中，r 和θ 分别为点信K 中某点x 与裂尖形成的法向量和切向量的夹角。

1.2 脑组织分割与网格化

使用图像分割算法将术前采集图像中的颅骨和皮肤等组织分割出来，仅保留包含脑组织部分的图像，并完成对脑组织的三维重建。在此基础上，利用类八叉树算法将脑组织的三维几何模型划分为仅由均匀的六面体单元组成的网格。

图2 XFEM 原理示意Fig．2 The schematic diagram of XFEM

1.3 XFEM 应用于线性弹性模型

线弹性模型是应用在脑组织变形矫正中最简单有效的生物力学模型［11］，它将脑组织看作是一个各向同性的线性弹性物体，用一系列准静态偏微分方程［12］表示。牵拉变形矫正的核心问题就是使用XFEM 方法求解线性方程，即

式中，K 是网格节点XFEM 刚度矩阵，P 是网格节点受力向量，a 是网格节点位移向量。

K 由两个参数决定:杨氏模量E 和泊松比v，选取杨氏模量E = 3 000 Pa 和泊松比v = 0．45［13］。

1.4 边界条件获取与应用

为了求解式(5)中a 的唯一解，需要边界条件约束，即已知部分网格节点的位移，包括受裂缝影响网格节点位移和零位移网格节点。受裂缝影响网格节点包括脑组织牵拉的初始位置，即切开脑组织位置，可以从术前图像获取;还有这些网格节点受裂缝影响产生的位移，可以通过对比牵拉前后扫描影像上的脑压板坐标变化。零位移网格节点为脑干部位网格节点，可以将两种边界条件引入式(5)进行迭代运算，求得牵拉变形后网格。

1.5 牵拉变形结果更新

用XFEM 计算得到的牵拉变形网格更新术前图像获得变形后的图像，以此改进了反插值算法［12］。原有算法只能处理网格节点的位移信息，能正确更新部分脑组织(非牵拉区域)图像，无法正确更新牵拉区域图像。改进的算法能够正确处理牵拉变形网格中的裂缝信息，所以能更新包括牵拉区域在内的所有牵拉变形后的脑组织图像。

从变形后的六面体单元出发，寻找出单元内所有像素点;根据单元内8 个节点位移场计算出这些像素点的变形位移场;再根据位移场，找到像素点在未变形前图像中的位置。利用八邻域内三线性插值获得像素点的灰度值，将此作为变形后该点的灰度值。经过对所有变形后的单元进行上述处理，从而获得变形后的三维数据场，并对其进行三维显示。

2 试验验证

2.1 脑组织体模概述

为验证所提出的系统框架，选用脑组织体模，并采用iCT 图像作为验证金标准。脑组织体模如图3 所示，重1 ～2 kg，由6%水凝胶PVA-C 材料制备，外形接近大脑形态，杨氏模量E ≅3 000 Pa，变形在20%内，应力应变曲线接近线性。

图3 脑组织体模Fig．3 The brain phantom

将脑组织体模固定于有机玻璃容器中，在容器外侧不同高度和底部表面粘贴共10 个2 mm 的塑料标记物，用于配准获取边界条件。在双侧体模的额叶、顶叶和枕叶部位，分别植入4、13 和7 枚直径从1 ～2 mm 不等的不锈钢标记物。其中，顶叶为牵拉区域，13 枚标志物都位于较浅的裂缝周围组织中，额叶和枕叶除了标记物23 和24，其他均为浅部标记物。这些标记物随体模的运动而运动，因此可以对比这些标记物在体模内牵拉变形前后的位移变化，对比XFEM 框架矫正前后图像中这些标记物的位移变化，评估模型的精度。

2.2 图像分割与网格化

在术前图像的感兴趣区域设定种子点，通过区域增长算法寻求其最大连通域，并手动填充未选中的感兴趣区域;再使用腐蚀算法处理，将脑组织与有机玻璃容器之间残存的细微连接去除;最后，重建获得仅包含脑组织的三维几何模型。为便于XFEM 建模，因此对三维模型表面进行平滑，以此来简化模型表面。然后进行六面网格化。

2.3 边界条件获取与应用

选择螺旋模式，0.5 mm 层厚，无间隔扫描的CT(Siemens SOMATOM Definition AS)序列。首先，人工在体模上切开一道裂缝，将此设定为脑压板初始位置，进行扫描作为术前图像;然后，插入脑压板分别向右侧牵拉5.7 mm，向左侧牵拉4.8 mm，总计10.5 mm，进行术中扫描。两次扫描都重建为(512×512 × 421)像素的断层图像。最后，利用容器上的标记物，将术中的图像配准到术前图像坐标系中，获得从术中图像变换到术前图像的变换关系，找到术中图像中的脑压板在术前空间中的坐标，比较变形前后图像中的脑压板坐标变化，获得脑压板的牵拉位移。此外，零位移节点从脑干部位网格表面手工选取。

在体模的几何模型进行六面体网格化的过程中，为研究不同尺寸网格对XFEM 计算矫正牵拉变形的影响，采用5 种尺寸网格来比较矫正牵拉变形的结果，分别为2.5 mm (70 370个节点，65 044个单元)、3.5 mm (28 899 个节点，26 098个单元)、4.75 mm (13 221 个 节 点，11 648 个 单 元)，5.0 mm(11 614个节点，10 183个单元)和7.5 mm (4 373个节点，3 689个单元)。

建立XFEM 框架的线性方程组，利用工作站(Intel Xeon E5540 处理器，12GB 内存)进行迭代运算，求得网格节点的位移和与裂缝相关信息。

2.4 图像更新与验证结果评估

XFEM 计算得到变形网格后，用改进的反插值算法更新术前图像，得到牵拉变形后图像。

模型的预测误差由XFEM 框架矫正的图像与实际变形图像中不锈钢标记物坐标差计算得出，模型的矫正精度由矫正后图像中不锈钢标记物的位移与真实位移的百分比得出，有

预测误差 =‖B - C‖2

矫正精度 =(‖B - A‖2- ‖C - A‖2)× 100%式中，A 为术前图像，B 为矫正变形后图像，C 为实际术中图像。

3 试验结果

使用XFEM 矫正10.5 mm 的牵拉变形，植入体模的24 个不锈钢标记物，实际变形大小为0 ～6 mm，平均(2.3 ±1.4)mm。

在利用基于XFEM 计算得到的变形后网格对术前图像进行更新时，改进的反插值算法能够识别节点位移信息和裂缝信息。如图4 所示，分别将改进的反插值方法和传统的反插值法更新得到的图像与实际牵拉图像进行对比;(a)是利用传统的反插值算法更新得到的图像，明显看出牵拉裂缝处的形态无法准确更新，(b)是改进了反插值算法更新得到的图像，可看出已准确地更新了牵拉裂缝图像;(c)是实际牵拉图像。可以看出，(a)与(c)在图像裂缝处的形态相差较大，而(b)与(c)在裂缝处的形态吻合较好。

图4 术中牵拉图像的更新结果。(a)传统反插值法更新图像的结果;(b)改进的反插值法更新图像结果;(c)实际牵拉图像Fig．4 The images updated by two algorithms． (a)The warped images updated by traditional backinterpolation algorithm; (b) The warped images updated by modified back- interpolation algorithm;(c)Real retraction images

表1 列出采用XFEM 框架矫正相同牵拉变形时网格尺寸对模型矫正精度的影响，网格尺寸分别为2.5、3.5、4.75、5.0、7.5 mm。尺寸最小的2.5 mm网格无法计算出牵拉变形的结果;3.5 mm 网格的预测误差为0.4 mm，矫正精度为(87.1 ± 16.4%);4.75 mm 网格的预测误差为0.4 mm，矫正精度为(85.8 ± 16.2)%;5.0 mm 网格的预测误差为0.3 mm，矫正精度为(87.5 ±16.0)%;尺寸最大的7.5 mm 网格的预测误差为0.5 mm，矫正精度为(83.1±20.3)%。从预测上看，5.0 mm 网格的预测误差最小，7.5 mm 网格的预测误差最大，平均预测误差为0.4 mm。从矫正精度上看，5.0 mm 网格的矫正精度最高，7.5 mm 网格的矫正精度最低，平均矫正精度为85.9%。从计算时间上看，3.5 mm 网格的计算时间为4 440 s，远远长于其他几种网格的计算时间;7.5 mm 网格的计算时间最短，仅仅60 s。

表1 不同尺寸网格对XFEM 矫正牵拉变形结果的影响Tab． 1 The correction accuracy and calculating time of different size meshes

图5 是尺寸为5.0 mm 的网格矫正牵拉变形后图像中不锈钢标记物的位移和实际术中图像中不锈钢标记物的位移对比统计结果。在图5(a)中，菱形表示尺寸为5.0 mm 的网格矫正牵拉后图像中不锈钢标记物位移，正方形表示实际术中图像中不锈钢标记物位移，三角形表示对应每一个不锈钢标记物的预测误差。图5(b)为尺寸为5.0 mm 网格中每一个不锈钢标记物的矫正精度，所有的均大于50%;标记物5 ～18 正好位于顶叶牵拉区域，其矫正精度处于波动中;最大位移出现在标记物13，零位移出现在枕叶标记物24 和25。

图6 为5.0 mm 网格建立的XFEM 框架矫正牵拉变形的三维结果。图6(a)是术前扫描得到的网格，包含了裂缝位置，即脑压板牵拉的初始位置。图6(b)是XFEM 框架矫正牵拉变形网格的云纹图，红色表示变形最大位置，蓝色代表无变形位置;可以看出，受到牵拉影响，网格裂缝周围出现的大小不同的位移，最大位移出现在裂缝周围。图6(c)是XFEM 矫正牵拉变形更新得到的图像三维重建结果。

图5 尺寸5.0 mm 的网格矫正牵拉变形后图像中不锈钢标记物位移与实际术中图像中标记物位移的统计结果。(a)5.0 mm 网格矫正变形后图像中不锈钢标记物位移，实际术中图像中不锈钢标记物位移和预测误差;(b)网格中每一个不锈钢标记物的矫正精度Fig．5 The statistical results of metal markers in the updated images with mesh size 5.0 mm comparing with markers in intraoperative images．(a)The displacement of metal markers in the updated images and intraoperative images，and the correction errors of each marker;(b)The correction accuracy of each metal marker

从脑压板末端形态与受裂缝影响网格节点分布关系，探讨不同网格尺寸对XFEM 矫正牵拉变形的影响。如图7 所示，黑色半圆形曲线表示脑压板末端的形态，黑色网格为局部网格，深黑色点为插入脑压板后受其产生裂缝影响的网格节点。图7(a)是3.5 mm 网格，节点的分布完全能够反映脑压板末端形态;图7(b)是4.75 mm 网格，节点的分布大致能够反映出脑压板末端形态;图7(c)是5.0 mm 网格，节点的分布大致能够反映脑压板末端形态;图7(d)是7.5 mm 网格，节点的分布无法正确反映脑压板末端形态。

4 讨论

笔者提出了一种将扩展有限元(XFEM)应用于线弹性生物力学模型来矫正牵拉变形的系统框架，并利用脑组织体模，实际验证了该牵拉变形矫正框架的有效性。验证结果表明，该系统框架能够有效地矫正术中的脑组织牵拉变形，提高IGNS 术中导航的精确性及可靠性。

图6 XFEM 矫正框架的工作流程。(a)术前网格;(b)矫正牵拉后变形网格;(c)矫正牵拉变形后图像三维重建结果Fig．6 Illustration of the workflow of our XFEM framework． (a)Pre-deformed mesh;(b)The deformed mesh updated by our framework;(c)The three-dimensional brain phantom visualized by ray casting algorithm

在利用XFEM 模拟计算出的变形网格更新图像时，改进的反插值算法与原算法相比，能准确地更新受牵拉部位的脑组织图像，尽管图像的外轮廓较为平滑，但并不影响牵拉变形矫正结果评估的客观性。这是因为脑压板是对脑组织内部施加力，而脑组织表面的位移是皮层下运动传导的结果，外轮廓平滑不会影响脑组织内部的位移。此外，模型矫正精度的评估主要依靠植入脑组织内部的不锈钢标记物，脑组织表面平滑不会影响皮层下标记物的位置。

验证试验中植入脑组织体模的不锈钢标记物主要位于较浅脑组织，尤其是顶叶部位。由于受牵拉的脑组织一般较浅，因此可以最大程度地捕获到脑组织牵拉对周围组织的影响。Vigneron 并没有采用植入标记物方法，而是通过矫正前后图像的Hausdorff 距离来评估矫正结果，这不能定量评价矫正结果的好坏，同时评价结果受到非刚性配准精度的影响。

图7 选取不同尺寸网格，利用XFEM 框架矫正牵拉变形，受裂缝影响网格节点分布与脑压板末端形态之间的关系。(a)3.5mm 网格;(b)4.75mm 网格;(c)5.0mm 网格;(d)7.5 mm 网格Fig． 7 The relationship between nodes fully intersected by the crack and the end of the retractors in different sizes of meshes． (a)3.5mm mesh; (b ) 4.75mm mesh; (c)5.0mm mesh;(d)7.5mm mesh

在网格化的过程中，网格尺寸的选取应该考虑到牵拉变形的大小。在试验中，单侧脑压板牵拉位移平均为5.2 mm，通过5 种不同尺寸网格模型矫正结果的对比，发现与牵拉位移最接近的5.0 mm 网格矫正精度最高。小于牵拉位移的3.5 mm 网格的矫正精度与5.0 mm 网格的矫正精度相近，而计算时间远长于5.0 mm 网格的结果;选用远小于牵拉位移的2.5 mm 网格，无法计算出牵拉变形结果;大于牵拉位移的7.5 mm 网格的矫正精度最低。因此，网格尺寸相对于牵拉位移过小时，较多的节点和单元成倍地延长了计算时间，同时线性方程组可能无法迭代达到收敛。在网格相对于牵拉位移过大时，虽然系统计算负担减小，但会导致网格节点刚度矩阵K 相对于网格节点受力矩阵P 过大;在求解线性方程组中网格节点位移向量a 时，分母远大于分子，导致计算得到的网格节点位移肯定小于真实位移。

网格尺寸的选择，不仅与牵拉变形的大小有关，而且与该尺寸网格中受脑压板影响网格节点形态和脑压板末端的形态一致性相关。在选取不同尺寸网格矫正牵拉变形时，虽然网格尺寸对矫正精度影响较小，但是随着网格尺寸增大，模型的矫正精度有下降趋势，初步认为这可能与插入脑压板后受裂缝影响的网格节点分布有关。网格中受裂缝影响网格节点的分布越接近脑压板末端形态，对这些节点施加力就越接近真实脑压板末端作用于脑组织的力。但是，当网格尺寸变大时，网格中受裂缝影响节点的分布越来越不能反映脑压板末端形态，对节点施加的力就越不能真实地反映脑压板对脑组织的作用力。

在试验中，使用脑组织体模来验证线弹性XFEM 框架矫正脑组织牵拉变形的有效性，但这与开颅手术中脑组织的牵拉变形情况还是存在差距。后续工作准备将矫正系统框架应用于动物试验和临床试验，以进一步验证该矫正系统的可靠性和实用性。此外，XFEM 框架使用iCT 来获得术中脑组织牵拉变形图像，虽然得到的边界条件非常准确，但是iCT 存在辐射且费用高，因此后续将利用简单可靠的术中设备(如三维激光扫描仪或术中超声)来代替iCT/iMR，在不降低本系统框架矫正精度的前提下，增加临床的实用性。

［1］ Zhang C，Wang M，Song Z. A brain-deformation framework based on a linear elastic model and evaluation using clinical data［J］. IEEE Trans Biomed Eng，2011，58(1):191 -199.

［2］ Vigneron LM，Duflot MP，Robe PA，et al. 2D XFEM-based modeling of retraction and successive resections for preoperative image update［J］. Comput Aided Surg，2009，14(1 -3):1 -20.

［3］ Ferrant M，Nabavi A，Macq B，et al. Registration of 3D intraoperative MR images of the brain using a finite-element biomechanical model［J］. IEEE Trans Med Imaging，2001，20(12):1384 -1397.

［4］ Miga MI，Roberts DW，Kennedy FE，et al. Modeling of retraction and resection for intraoperative updating of images［J］.Neurosurgery，2001，49(1):75 -84，84 -85.

［5］ Lamprich BK，Miga MI. Analysis of model-updated MR images to correct for brain deformation due to tissue retraction［C］ //Medical Imaging 2003:Visualization，Image-Guided Procedures，and Display. San Diego:SPIE，2003:552 -560.

［6］ Sun H，Kennedy F，Carlson E，et al. Modeling of Brain Tissue Retraction Using Intraoperative Data［C］ //Barillot C，Haynor D，Hellier P，eds. Medical Image Computing and Computer-Assisted Intervention – MICCAI 2004. Berlin Heidelberg:Springer，2004:3217，225 -233.

［7］ Duflot M，Nguyen-Dang H. A meshless method with enriched weight functions for fatigue crack growth［J］. International Journal for Numerical Methods in Engineering，2004，59(14):1945 -1961.

［8］ Moës N. Dolbow J，Belytschko T. A finite element method for crack growth without remeshing［J］. International Journal for Numerical Methods in Engineering，1999，46(1):131 -150.

［9］ Vigneron LM，Noels L，Warfield SK，et al. Serial FEM/XFEMbased update of preoperative brain Images using intraoperative MRI［J］. Int J Biomed Imaging，2012，2012:872783.

［10］ Vigneron LM，Warfield SK，Robe PA，et al. 3D XFEM-based modeling of retraction for preoperative image update ［J］.Comput Aided Surg，2011，16(3):121 -134.

［11］ Wittek A，Hawkins T，Miller K. On the unimportance of constitutive models in computing brain deformation for imageguided surgery［J］. Biomech Model Mechanobiol，2009，8(1):77 -84.

［12］ Liu Y，Song Z. A robust brain deformation framework based on a finite element model in IGNS［J］. Int J Med Robot，2008，4(2):146 -157.

［13］ Platenik LA，Miga MI，Roberts DW，et al. In vivo quantification of retraction deformation modeling for updated image-guidance during neurosurgery［J］. IEEE Trans Biomed Eng，2002，49(8):823 -835.