基于人体结构先验信息的胸部电阻抗成像方法

王 琦 陈晓静 汪剑鸣* 李秀艳 段晓杰 王化祥

1(天津工业大学电子与信息工程学院，天津 300387)2(天津市光电检测技术与系统重点实验室，天津 300387)3(天津大学电气与自动化工程学院，天津 300072)

引言

电阻抗成像(electrical impedance tomography,EIT)是临床使用的一种新的成像工具，它的基本原理是：通过配置于人体体表的一组电极阵列，施加一定频率的低伏交变安全电流；再通过扫描电极测量电压数据，提取与人体生理、病理状态相关组织或器官的电阻抗信息；经数据采集单元并行处理后，送至重构计算机；运用求解被测物场电磁问题的图像重建算法，重建被测组织和器官的二维/三维图像[1-2]。由于EIT具有非侵入性、实时性、设备简单、功能性成像的优势，因此在医学、地质勘探和工业过程成像领域获得广泛应用[3]。

在当代医学中，人体肺部 EIT 检测，对于发现某些肺部疾病具有十分重要的科学研究意义和临床应用价值[4-5]。EIT 不仅可以检测肺通气能力，而且可提供与人体胸腔生理和病理变化直接相关的功能图像。该技术比应用肺活量计和流量计更方便，更适合于长期的临床监护，成为 ICU 病房中辅助呼吸机进行呼吸监测的有效手段。此外，由于EIT 具有无辐射、实时性强及便携性等优点，因此在肺损伤病人早期诊断和床旁实时监测方面具有不可比拟的优势[6]。

EIT问题包括正问题和逆问题。正问题是由生物模型内部的阻抗分布及边界激励信号，求物体内部或表面的电压和电流分布；逆问题是通过表面的电压、电流分布及边界激励信号，求生物体模型内部的阻抗分布。图像重建过程属于逆问题，它是EIT的核心技术，图像重建模型中的灵敏度矩阵构建依赖于正问题求解。通常采用有限元剖分方法求解电阻抗成像正、逆问题[7]，其中正问题通过网格剖分获得场域节点电势，逆问题通过网格剖分获得生物体模型内部电导率分布[8]。网格剖分需要已知重建对象的边界形状，才能得到更加逼近生物体真实电阻抗分布的剖分模型，进而提高重建图像的质量[9-10]。建立不同边界形状的胸部模型会影响成像的精度，因此构造符合人体真实边界形状的胸部模型是求解EIT问题的前提。目前国内外对于胸部电阻抗的研究多用圆形代替胸部结构，Saka提出，在实际EIT应用中，根据生物体边界形状建立椭圆几何模型，比圆形几何模型更能准确分析EIT问题[11]。严佩敏等用椭圆形模型代替胸部结构进行电阻抗成像[12]，但是由于胸部轮廓具有特异性，用近似人体胸部形状形成的统一模型[13-14]会引入测量误差以及成像误差。因此，获得正确反映人体胸部真实结构的电阻抗成像方法仍是一个需要探究的问题。

本研究提出一种基于人体结构先验信息的胸部电阻抗成像方法。通过评价指标——肺部区域比例(lung regional ratio，LRR)，评价该方法、基于椭圆形模型成像方法、基于圆形模型成像方法的成像效果，并用仿真实验进行了验证。结果表明：基于人体结构先验信息的胸部电阻抗成像方法能获得与人体结构更为接近的肺部区域比例，有效降低传统模型成像与人体真实胸部结构成像的误差，提高成像效果，更加真实地反映人体胸部结构。

1 方法

1.1 正问题建模方法

基于先验信息的正问题建模包括两部分。第一部分是对屏息时获得的胸部CT图像进行图像处理，得到人体胸、肺部形状和轮廓相对准确的信息，尤其是边界信息，为构建较准确的胸、肺部模型奠定基础，图1即为电极所在中心平面对应的胸部CT图像。第二部分是将图像处理后提取的轮廓边界信息转化为正问题模型。

1.1.1图像处理

为了准确地从CT图像中提取有用信息，需要对原始CT图像进行处理，以突出有用的图像信息。首先对图像进行预处理，包括去噪、滤波等；其次通过图像分割提取出感兴趣区域，该效果会直接影响到重建后模型的准确性；最后对图像分割后的图形进行轮廓提取，作为仿真软件模型的场域轮廓线。具体的图像处理流程如图2所示。

图1 胸部CT图像Fig.1 Chest CT tomography

图2 图像处理流程
Fig.2 Image processing flow chart

1)图像预处理。为了有效消除边缘噪声，提高轮廓线提取精度，本研究采用基于偏微分方程(partial differential equation′s，PDE′s)的图像方法。目前，以Perona&Malik模型[15]为代表的各向异性平滑方法已在边缘检测、图像增强、图像分割以及目标识别等领域得到广泛应用。本研究采用改进的Perona&Malik模型[15]，设u0(i,j)为原始CT灰度图像，引入时间变量t∈[0,T]，在求解偏微分方程的迭代过程中，将得到一组逐渐平滑的图像u(i,j,t)，即

(1)

式中：Gσ是高斯平滑模板；|u|是图像梯度模；扩散系数c(·)是关于|u|的非负函数，用于控制扩散速度。

2)阈值分割。为了得到肺部及胸部区域的轮廓线，需要对图像预处理后的CT图像进行阈值分割。本研究采用一种基于图像形态学中腐蚀、膨胀的阈值分割算法[16]。CT图像的分割步骤是：首先对预处理后的CT图像进行二值化处理，然后去除检查床提取胸部轮廓，再利用去除检查床的图像减去胸廓图像，即可得到肺部轮廓图像。最后通过对肺部轮廓和胸部轮廓的融合，得到肺部及胸部边缘。具体的阈值分割流程如图3所示。

图3 阈值分割流程Fig.3 Threshold segmentation flow chart

图4 胸部及肺部轮廓曲线Fig.4 Chest and lung contours

3)提取轮廓线。为了提高轮廓提取的精准度，本研究采用基于二值图像的轮廓提取算法[17]。设e(i,j)为二值化后CT图像的像素点，E(i,j)为应用规则处理后的CT图像像素点，判断当前像素是否为边界像素的规则如下：

(2)

E(i,j)=0时，当前像素不是边界像素，删除；E(i,j)=1时，当前像素是边界像素，予以保留。应用此规则对二值化CT图像的每个像素进行处理，保留在图像中的像素为图像边界，即可完成轮廓提取，最终得到较为平滑的胸部及肺部轮廓曲线如图4所示。

1.1.2构建模型

图5 胸部仿真模型Fig.5 Chest simulation model

基于本文第1.1.1节获得的胸部及肺部轮廓图像，构建胸部仿真正问题模型。将胸部及肺部轮廓曲线转化为CAD矢量图，运用有限元仿真软件Comsol，分别建立基于人体结构先验信息的胸部仿真模型、基于椭圆形模型的胸部仿真模型和基于圆形模型的胸部仿真模型，进行EIT正问题计算，从单元划分到计算都可以在短时间内得到计算结果数据。需要指出的是，对于圆形仿真模型，其直径是CT图片中最高点与最低点坐标之差；对于椭圆形仿真模型，其长轴是CT图片中X轴最大值与最小值坐标之差，短轴是CT图片中最高点与最低点坐标之差。3种重构模型的电导率是一致的，正问题初始电阻抗分布：胸腔是均匀分布的，不考虑内部器官组织。设置胸腔组织电导率为0.024 S/m，设置肺部的电导率为0.1 S/m[18]。胸腔组织电导率为0.024 S/m，脂肪电导率为0.024 S/m[18]仿真实验采用的是相邻驱动模式，对于有N个电极的EIT系统，相邻模式可以采集到N×(N-3)边界电压。本研究采用16个电极，胸部仿真模型如图5所示。

1.2 逆问题剖分方法

逆问题剖分要解决两个问题：首先，要有一个有效的图像边界，该边界可以通过正问题模型获取；其次，在该边界的内部进行逆问题剖分。为了避免逆问题风暴(inverse crime)[19]，正逆问题采用不同的网格剖分，本研究的正问题采用三角剖分，逆问题采用矩形剖分，即每个剖分单元为矩形。

1.2.1构建n×n的矩形网格

(3)

(4)

图6 n×n矩形网格Fig.6 n×nrectangular grid

图7 胸部模型上下两部分边界坐标点Fig.7 Chest model of the upper and lower parts of the boundary coordinates

1.2.2边界线性插值坐标

由本文第1.2.1节可知，剖分网格可以把模型沿X轴方向做n等分，获得剖分点横坐标集合A={x1,x2,x3…xn}，则模型上下两部分边界剖分点坐标集合分别为B={(x1,yU1),(x2,yUM2),…，(xn,yUn)}和C={(x1,yD1),(x2,yD2),…，(xn,yDn)}，其中yU1,yU2,…，yUn和yD1,yD2,…，yDn分别是上下边界剖分点的纵坐标。边界剖分点纵坐标可根据已知的边界剖分点横坐标A及边界点坐标(xu,yu)∈Rk1×2和(xd,yd)∈Rk2×2通过插值方法获得。本研究采用线性插值的方法，上下两部分边界剖分点纵坐标均需要通过插值方法得到。以获得上边界剖分点纵坐标为例：设(xui,yui)和(xu(i+1),yu(i+1))是两个相邻的边界点，若两边界点间隔s个边界剖分点，则每个边界剖分点的纵坐标如下：

yU(m+l)=yui+k(x(m+l)-xui)

(5)

最终得到上下两部分边界剖分点坐标集合B和C，线性插值后的坐标如图8红点所示。模型内部的坐标可通过模型边界坐标获得，具体方法见本文第1.2.3节。

图8 线性插值后的坐标点Fig.8 Coordinates after linear interpolation

线性插值后的边界坐标点围成的轮廓与原图像轮廓的对比如图9所示。其中第一行是三个CT样本，第二行中黑色线代表CT图像边界轮廓，红色点表示线性插值后的坐标点。由图9可知，线性插值后的边界坐标点围成的轮廓与原轮廓基本一致。

图9 线性插值后的边界坐标点围成的轮廓与原轮廓对比Fig.9 Comparison between the original contour and the contour enclosed by the boundary points after linear interpolation

图10 基于人体结构先验信息的胸部电阻抗成像方法的有效网格筛选方法Fig.10 An effective mesh selection for chest electrical impedance tomography based on prior information of human body

1.2.3筛选有效网格

将包含在模型边界点坐标范围内的网格定义为有效网格,定义变量H(Xi,Yi)为网格的有效值，有

(6)

若H(Xi,Yi)的值为1，说明该矩形网格在模型边界内部记为有效网格；若H(Xi,Yi)的值为0，则其为无效网格。

通过对图9中的样本1进行边界获取及剖分，然后筛选得到有效网格，如图10所示。其中，灰色网格为有效网格，白色网格为无效网格。将有效网格的坐标数据、编号以及数量保存下来，用于后期的图像重建。为了更加直观地显示网格剖分过程，绘制流程如图11所示。

图11 基于人体结构先验信息的胸部电阻抗成像方法的网格剖分流程Fig.11 Mesh generation flow chart of chest electrical impedance tomography based on prior information of human body

网格精细程度将决定图像分辨率，即n值越大，图像分辨率越高。本研究采用n=32的分辨率坐标值，共1 024个矩形剖分单元。图12(a)为剖分后绘制的矩形网格线，其中每个小矩形就是一个剖分单元，它比较直观地显示剖分效果，且能较直观地描绘胸部模型的形状。为了对比剖分效果，64×64的矩形网格线如图12(b)所示，128×128的矩形网格线如图12(c)所示。由上述剖分网格线可知，在有限元剖分时，随着网格精细程度的提高，图像轮廓更加接近真实场域形状。

图12 不同精细程度的剖分网格。(a)32×32；(b)64×64；(c)128×128矩形剖分网格Fig.12 Different square mesh.(a)32×32；(b)64×64；(c)128×128

1.3 评价方法

为评估不同方法构成的胸部模型与人体胸部真实结构的差异，引入一种基于重建图像结果的评价指标——肺部区域比例(lung regional ratio，LRR)，其定义如下：

(7)

式中，TV表示肺部通气面积，SV表示胸部截面除肺部以外区域的面积[20]。

若用F∈Rl×1、E∈Rl×1分别表示肺部区域像素点的集合、胸部截面除肺部以外区域像素点的集合(l表示有效的剖分点数)，则集合F、E可分别表示为

(8)

(9)

式中：t(x,y)是像素点(x,y)的像素值；T为肺部区域阻抗阈值，表示肺部区域最低阻抗值。

那么，TV和SV可分别表示为

(10)

(11)

本研究采用最大类间方差算法[21]确定区分肺部区域和非肺部区域的阈值。设f(x,y)为重建图像的位置(x,y)处的灰度值，具有L个灰度级(本研究中L选为20，即将归一化的电导率范围平均分为20份)，灰度级用变量b表示。设灰度级为bi的所有像素个数为fi，其中i=0,1,…,L-1，且b0

(12)

选择第t个灰度级作为图像像素的阈值，用于区分胸腔C0和肺部组织C1，此时胸腔C0的灰度级包括b0,b1，…，bt-1，肺部组织C1的灰度级包括bt,bt+1，…，bL-1，则胸腔与肺部组织的分割阈值bt可表示为

bt=arg.max{∂2(t)} (0≤t

(13)

式中，∂2(t)是类间方差，可表示为

∂2(t)=ω0(λ0-λ)2+ω1(λ1-λ)2

(14)

分别计算不同成像方法重建图像的LRR，与基于CT图像计算的人体真实LRR进行比较，并计算相对误差，用来评价本研究提出方法的有效性。相对误差的计算公式为

(15)

式中,l为成像方法计算得出的LRR，L为基于CT图像计算的人体真实LRR。

笔者共研究30个肺部健康的人体CT图像样本，样本取自5个年龄段的人群，分别为：20～30，30～40，40～50，50～60，60～70岁；每个年龄段共6个样本，其中男性、女性各3个。每个样本分别采用基于人体结构先验信息的胸部电阻抗成像方法、基于椭圆形模型成像方法、基于圆形模型成像方法进行成像。在重建图像中，重构断面为电极所在的中心平面，重建的灰度值根据电导率的真实值归一化至[1,3]范围内。由于动态成像方法具有相对简单、测量误差小、成像质量更高的优点[22]，所以本研究采用基于共轭梯度的动态成像方法来求解电阻抗成像的逆问题。逆问题的剖分及成像采用Matlab2012a软件，所有的结果在处理器为Intel(R)Core(TM)i5-3 230 M CPU、内存为4GB的计算机上获得。

1.4 显著性检验

对研究的30个样本进行显著性检验，根据实验的特点使用配对样本t检验的方法[23]，分别对3种成像方法计算的LRR值以及3种成像方法与人体真实LRR值的相对误差做显著性检验。以3种成像方法算出的LRR值为例，显著性检验具体步骤如下：

1) 建立假设检验，确定检验水准。H0表示基于成像方法算出的LRR与真实值不存在显著性差异，H1表示基于成像方法算出的LRR与真实值存在显著性差异。双侧检验，检验水准α=0.05。

2) 计算统计量为

(16)

3) 确定P值，做出统计推断。运用统计软件确定P值：若P>α，则H0假设成立；若P<α，则H1假设成立。

2 结果

用基于不同模型的成像方法，对样本的图像进行重建比较。以图9中的3个样本为例，3种成像方法重建图像结果如图13所示，其中分别用黑色实线、红色实线标出了肺部区域轮廓、胸部区域轮廓。不同成像方法的统计学分析如表1所示。

图13 不同成像方法重建图像结果Fig.13 Reconstruction of image results with different imaging methods(From the left to the right: sample 1～sample 3)

从图13的结果可看出，边界形状对成像效果有一定影响。对比3种方法重建图像的结果，发现存在明显的差异，其中基于人体先验信息的胸部电阻抗成像方法可获得更加接近人体胸部的边界形状。

由表1可见，对于3种基于不同模型方法成像结果的LRR均与基于CT图像计算的人体真实LRR存在差异，但是基于人体结构先验信息的胸部电阻抗成像方法得到的LRR更加接近真实值。经计算可知：基于人体结构先验信息的胸部电阻抗成像方法与人体真实LRR的相对误差最小。通过显著性分析发现：基于人体结构先验信息的胸部电阻抗成像方法得到的LRR与基于CT图像计算的人体真实LRR相比，P>0.05，大于显著性水平；基于椭圆形模型成像方法、基于圆形模型成像方法与基于CT图像计算LRR相比，P<0.05；基于椭圆形模型成像方法、基于圆形模型成像方法与基于人体结构先验信息的胸部电阻抗成像方法的相对误差相比，P<0.05，均小于显著性水平。由此说明，基于人体结构先验信息的成像方法能有效降低成像误差，提高成像效果。

表1 不同成像方法计算LRR统计学分析(均值±标准差)Tab.1 Statistical analysis of LRR by different imaging method (mean±SD)

注：与LRR真实值相比,a:P>0.05, b:P<0.05; 与基于人体结构先验信息成像方法的相对误差相比，c:P<0.05

Note: Compared with real LRR, a:P>0.05, b:P<0.05; Compared with the relative error of imaging method based on priori information of human body structure, c:P<0.05

3 讨论

EIT在人体胸部应用中，边界形状的不确定性是一个难点，由于胸部轮廓具有特异性，建立不同边界形状会影响性能和重建效果[24-25]，所以获取真实的人体胸部边界形状是求解EIT问题的前提。

在胸部电阻抗成像问题中，建立符合人体真实胸部及肺部边界形状的模型并进行有效剖分，对成像效果有着积极作用。传统圆形模型、椭圆形模型成像方法均不能真实反映人体胸部结构，成像误差较大。所以，如何正确获取人体胸部及肺部边界形状来提高成像效果是一个难点。本研究提出的基于人体结构先验信息的胸部电阻抗成像方法，通过对CT图片进行图像处理提取胸部及肺部轮廓，为正问题和逆问题提供图像边界的先验信息，基于边界先验信息进行正问题胸部建模，避免了统一模型引入的测量误差。同时，基于边界先验信息提出一种有效的图像逆问题剖分方法，在减小误差的同时提高成像效果，从而使重建图像形状更接近真实情况。

研究结果反映出所提出的基于人体结构先验信息的胸部电阻抗成像方法的有效性，与传统基于椭圆形模型成像方法、基于圆形模型成像方法相比，该方法在成像效果上更加接近人体胸部边界形状；不同成像方法的评价指标——肺部区域比例(LRR)的计算结果表明，该方法可获得与人体真实结构更接近的肺部区域比例。所提出方法的LRR与基于CT图像计算的真实LRR之间的相对误差(3.71%±1.77%)显著小于基于椭圆形模型(10.29%±3.30%)和基于圆形模型(12.74%±2.87%)这两种成像方法。该方法显著降低了传统模型成像方法与人体真实胸部结构的成像误差，成像效果得到明显提高。由于本研究提出方法的标准差与真实LRR相对误差的标准差最小，所以该成像方法更稳定。分别就这三种成像结果的LRR值及其与真实值的相对误差进行统计学分析，结果表明：基于人体结构先验信息的胸部电阻抗成像方法与真实的LRR值不存在显著性差异，基于椭圆形模型成像方法、基于圆形模型成像方法与真实的LRR值存在显著性差异，基于椭圆形模型成像方法、基于圆形模型成像方法与基于人体结构先验信息的胸部电阻抗成像方法的相对误差存在显著性差异。这是由于胸部轮廓具有特异性，用圆形模型或椭圆形模型代替人体胸部模型会引入测量误差以及成像误差。基于人体结构先验信息的胸部电阻抗成像方法可有效解决胸部EIT图像重建的边界轮廓特异性问题，提高成像质量。

在本研究中，利用屏息状态下的CT图像获取胸部及肺部的边界信息。由于呼吸过程中肺部区域、胸部区域会随着呼吸而变化，所以无法获取确定的轮廓边界信息。在下一步工作中，将考虑呼吸过程对胸廓形状的影响，融入形变先验信息，并且进一步提高剖分精度，改进成像质量。

4 结论

本研究提出一种基于人体结构先验信息的胸部电阻抗成像方法，通过现有图像处理方法获得胸部、肺部轮廓，为正问题和逆问题提供图像边界的先验信息，同时基于边界先验信息提出一种有效的图像逆问题剖分方法，从而使重建图像形状更接近人体胸部真实情况，有效降低用圆形、椭圆形代替胸部结构或用近似人体胸部形状形成统一模型引入的误差，改善成像效果。