CN103654776A - 融入颅骨电阻率非均匀分布信息的电阻抗断层成像方法 - Google Patents

Abstract

本发明涉及一种融入颅骨电阻率非均匀分布信息的电阻抗断层成像方法，利用被测对象头部CT数据提供的解剖结构信息以及颅骨板障层厚度同颅骨电阻率值的统计参数模型建立包含颅骨电阻率非均匀分布的头部二维有限元模型，并基于此模型实现电阻抗断层图像重建。该方法将颅骨组织电阻率的真实分布信息快速、自动地融入成像算法，修正了图像重构的结果，改善了成像目标的定位精度、图像空间分辨率，提高了头部电阻抗断层成像的图像质量，改进了头部电阻抗成像实际使用中的易用性。

Description

融入颅骨电阻率非均匀分布信息的电阻抗断层成像方法

技术领域

本发明属于电阻抗断层成像技术领域，涉及一种通过电阻抗断层成像技术来获取被对象体内阻抗分布图像的方法，特别涉及一种融入颅骨电阻率非均匀分布信息的电阻抗断层成像方法。

背景技术

电阻抗断层成像技术（Electrical Impedance Tomography，EIT）是继形态、结构成像之后新一代的无损伤功能成像技术。其基本原理是根据人体内不同组织在不同生理、病理状态下具有不同电导率（电阻率），通过在人体表面施加安全驱动电流（电压），在体外测量响应电压（电流）信号进而重建出反映人体内部电导率分布或变化的图像。中国专利申请（申请号：99115855.5），题为《一种电阻抗断层成像方法》，中国专利申请（申请号：03134598.0），题为《一种用于床旁图像监护的电阻抗断层成像方法》，对此技术方案进行了详细描述，本申请以此作为相关参考文献。

根据成像目标不同，EIT可分为动态成像和静态成像两种方式。动态EIT是对阻抗变化的相对值成像，可以在测量过程中抑制共同噪声的同时，利用人体阻抗的变化直接反映出人体生理状态或病理过程这一特点，为临床提供具有高时间分辨率的图像监护，这是X射线计算机断层成像和核磁共振成像无法比拟的优点，因此电阻抗动态图像监护技术的研究具有重要的应用价值和意义。目前，EIT技术已应用于腹腔出血图像监护、肺通气图像监护、脑功能图像监护等方面的研究。而静态成像，能够重构出绝对电阻抗分布，更具有临床应用价值，但由于EIT逆问题的病态性问题，少许的重构输入误差则会引起很大的成像误差甚至导致成像失败，所以静态重构方法仍在理论研究阶段。

在先前的研究中，EIT成像方法通常使用圆域作为成像场域，并将初始电阻率分布近似为均匀分布，这种假设可以控制成像场域有限元模型剖分规模、简化计算，并减少逆问题病态性对重构的影响。但随着EIT应用于人体测量，发现缺失被测对象内部解剖结构和阻抗分布信息，会造成较大的模型误差，以至重构图像与真实情况不符，无法满足应用需求。因此，有EIT系统通过医学结构成像技术获取人体解剖结构信息融入成像算法，如中国专利申请（申请号：200910022777.7），题为《一种结构信息融合的电阻抗断层成像方法》所述，将头部各组织层解剖结构融合入成像方法，修正EIT重构计算误差。但在此方法中，并未考虑颅骨层电阻率分布的特点对EIT成像的影响。

因颅骨组织学结构的特异性，其阻抗特性具有以下特点：颅骨电阻率明显高于其它头部组织，对流入颅腔的电流具有很大的阻碍作用；六类颅骨结构（标准三层骨、准三层骨、标准密致骨、准密致骨、齿状骨缝、鳞状骨缝）的电阻率差异显著；标准三层骨和准三层骨的外侧密致骨（皮质骨）同中间板障层（髓质骨）电阻率差异显著，其电阻率由板障层的厚度决定，而板障层厚度分布不均一，不同部位、不同个体均不一样；板障层厚度同颅骨电阻率之间存在显著的统计关系（Tang,C.,et al.,Correlation between structure and resistivity variations of the livehuman skull.IEEE Trans.Biomed.Eng,2008.55（9）:p.2286-2292.）（Tang,C.,et al.,Modeling the frequency dependence of theelectrical properties of the live human skull.Physiol Meas,2009.30（12）:p.1293-301.）。以上特点决定了颅骨电阻率的非均匀分布。如果在EIT成像中，将颅骨近似为电阻率均匀分布的一层结构，那么重构计算中融入的颅腔内电流密度分布形式同被测对象的真实形式将完全不同，造成EIT图像空间分辨率退化，成像目标定位不准确，严重制约了EIT技术的临床应用。

综上所述，在头部EIT成像中，为提高成像目标的定位精度和空间分辨率、改善重构图像的质量，需要一种能够融入被测对象个体特异的颅骨电阻率非均匀分布信息的电阻抗成像方法。

发明内容

本发明的目的在于，提供一种融入颅骨电阻率非均匀分布信息的电阻抗断层成像方法，以克服本邻域图像重构过程中因颅骨电阻率分布信息的缺失造成的模型误差，提高成像定位精度和空间分辨率、改善重构图像质量。

为实现上述任务，本发明采取如下的技术解决方案：

一种融入颅骨电阻率非均匀分布信息的电阻抗断层成像方法，其特征在于，按如下步骤进行：

步骤1，获取被测对象头部电阻抗成像测量电极所在层面的轴位CT扫描数据；

步骤2，根据CT扫描数据自动划分头部的解剖结构，取得头部组织的分层结构信息；

步骤3，对颅骨层CT数据进行处理，获取密质骨和板障层分布信息；

步骤4，将头部组织的分层结构信息作为限定条件构建头部逼真有限元模型；

步骤5，配置颅骨电阻率非均匀分布以及头部其它组织的电阻率分布信息；

步骤6，由EIT测量系统采集头部实时阻抗测量数据，重构出融入颅骨电阻率非均匀分布信息的EIT图像并显示。

所述的获取测量电极所在层面的轴位CT扫描数据是通过在被测对象头部表面粘贴CT成像标记物确定头部轴位CT序列数据中EIT测量电极所在的层面。

所述的自动提取头部组织的分层结构信息，其具体步骤如下：

1）对CT图像采用直方图约束的自适应多阈值模糊C均值聚类分割算法提取颅骨层结构信息；

2）对CT图像中颅骨内侧区域采用非局部均值修正的自适应模糊C聚类分割算法提取脑室结构信息；

3）对颅骨内侧边界采用二值位图操作，自动绘制脑脊液层结构信息。

所述的获取颅骨层密质骨和板障层分布信息是指使用非局部均值修正的自适应模糊C聚类分割算法对颅骨层进行处理，提取密质骨、板障层和骨缝的结构分布信息。

所述的构建头部逼真有限元模型是指利用头部组织的分层信息、基于质量及尺寸可控的自适应限定Delaunay三角网格剖分技术构建用于头部EIT重构计算的有限元模型。

所述的配置头部组织电阻率分布信息，具体步骤包括：

①利用颅骨板障层厚度同颅骨电阻率值的统计参数模型，对颅骨层单元的电阻率属性进行赋值，实现颅骨电阻率非均匀分布信息的融入；

②根据人体头部组织电阻抗频谱特性测量结果，对头部其它组织层单元的电阻率属性进行赋值。

所述的由EIT系统采集被测对象头部实时阻抗测量数据，重构出融合头部先验信息的EIT图像并显示，具体方法为：利用头部EIT测量系统实时采集的数据和基于标准Tikhonov正则化方法的最小二乘EIT重构算法进行图像重构计算，将计算得到的电阻率分布利用云图方式显示。

本发明的融入颅骨电阻率非均匀分布信息的电阻抗断层成像方法，能够应用于电阻抗成像系统的临床脑部功能监测、头部物理模型实验、计算机仿真中。实现了被测对象的具有颅骨电阻率非均匀分布和真实解剖结构的头部有限元模型的构建、融合先验信息的电阻抗断层图像重建。所带来的优点是：通过融入颅骨电阻率非均匀分布和头部组织空间结构信息，改善了EIT成像目标的定位精度、图像空间分辨率，提高了图像质量，以更准确地反应被测对象脑部组织的电阻抗变化。

附图说明

图1是本发明融入颅骨电阻率非均匀分布信息的电阻抗断层成像方法的原理框图；

图2是图像分割聚类中心自动选取的流程图；

图3是提取头皮和颅骨组织分布边界信息的流程图；

图4是提取脑脊液层和脑室组织分布边界信息的流程图；

图5是得到的头部各层组织分布边界信息；

图6是构建头部逼真有限元模型的流程图；

图7是配置颅骨电阻率非均匀分布信息的流程图；

图8是融入颅骨电阻率非均匀分布的头部逼真有限元模型；

图9是本发明实施例的实验物理模型；

图10是融入颅骨电阻率非均匀分布信息的电阻抗断层成像结果；

以下结合附图和实施例对本发明作进一步的详细说明。

具体实施方式

本发明的融入颅骨电阻率非均匀分布信息的电阻抗断层成像方法，具体实施步骤为：

步骤1，获取被测对象头部电阻抗成像测量电极所在层面的轴位CT扫描数据；

步骤2，根据CT扫描数据自动划分头部的解剖结构，取得头部组织的分层结构信息；

步骤3，对颅骨层CT数据进行处理，获取密质骨和板障层分布信息；

步骤4，将头部组织的分层结构信息作为限定条件构建头部逼真有限元模型；

步骤5，配置颅骨电阻率非均匀分布以及头部其它组织的电阻率分布信息；

步骤6，由EIT系统采集被测对象头部实时阻抗测量数据，重构出融入颅骨电阻率非均匀分布信息的EIT图像并显示。

步骤1中，获取被测对象头部电阻抗成像测量电极所在层面的轴位CT扫描数据，其操作方法是：使用CT成像标记物标定被测对象头部EIT电极位置并行头部轴位CT扫描，将EIT系统连接入CT工作站使用的PACS系统，获取被测对象DICOM格式的头部轴位CT扫描数据；选取CT成像标记物即电极标记位置层面的CT图像。

步骤2中，提取头部解剖结构先验信息，其具体步骤如下：

①对CT图像采用直方图约束的多阈值模糊C均值（Fuzzy C-Means，FCM）聚类分割算法分割出头皮层和颅骨层，其方法是：

a、按照以下公式⑴统计CT图像的灰度直方图：

$H\left(k\right)=\frac{1}{\mathrm{MN}}\underset{i=0}{\overset{M-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{N-1}{\mathrm{\Σ}}}\mathrm{\δ}(f(i,j)-k)\left\{\begin{array}{c}k\∈\{\mathrm{0,1,2},...,G-1\}\\ \mathrm{\δ}\left(0\right)=1\\ \mathrm{\δ}(k\≠0)=0\end{array}\right.---\left(1\right)$

式中：M、N为图像像素尺寸，G为图像中最大灰度级。

将灰度直方图按以下公式⑵转换为势直方图：

$P_H\left(k\right)=\frac{1}{P_{\mathrm{MAX}}}\underset{g=0}{\overset{G-1}{\mathrm{\Σ}}}\frac{H\left(g\right)}{1+\mathrm{\α}{(g-k)}^2}---\left(2\right)$

式中：P_MAX＝max(P_H(g),gò{0,1,2,…,G-1})，α为趋近系数。

由图像势直方图提取最大峰值对应的灰度值和所有局部峰值中对应的最大灰度值作为直方图约束的多阈值FCM分割算法的初始分类中心ν。

b、按以下公式（3）、（4）计算约束系数ω和隶属度μ

${\mathrm{\&omega;}}_{\mathrm{ik}}=\left\{\begin{array}{cc}\frac{k-{\mathrm{\&upsi;}}_{i-1}}{{\mathrm{\&upsi;}}_i-{\mathrm{\&upsi;}}_{i-1}},& {\mathrm{\&upsi;}}_{i-1}\&le;k\&le;{\mathrm{\&upsi;}}_{i}i=\mathrm{2,3},...,C\\ \frac{{\mathrm{\&upsi;}}_{i+1}-k}{{\mathrm{\&upsi;}}_{i+1}-{\mathrm{\&upsi;}}_i},& {\mathrm{\&upsi;}}_i\&le;k<{\mathrm{\&upsi;}}_{i+1}i=\mathrm{1,2},...,C-1\\ 1,& {\mathrm{\&upsi;}}_0\&le;k<{\mathrm{\&upsi;}}_1\mathrm{or}{\mathrm{\&upsi;}}_C\&le;k<{\mathrm{\&upsi;}}_{C+1}\\ 0,& \mathrm{otherwise}\end{array}\right.---\left(3\right)$

${\mathrm{\&mu;}}_{\mathrm{ik}}=\frac{{\left[\frac{{\mathrm{\&omega;}}_{\mathrm{ik}}}{{(k-{\mathrm{\&upsi;}}_i)}^2}\right]}^{\frac{1}{m-1}}}{{\mathrm{\&Sigma;}}_{p=1}^{G-1}{\left[\frac{{\mathrm{\&omega;}}_{\mathrm{pk}}}{{(k-{\mathrm{\&upsi;}}_p)}^2}\right]}^{\frac{1}{m-1}}}\left\{\begin{array}{c}k=\mathrm{0,1},...,G-1\\ i=\mathrm{1,2},...,C\end{array}\right.---\left(4\right)$

式中：m为模糊系数。

c、将约束系数ω和隶属度μ代入公式（5）计算新的聚类中心，

${\mathrm{\&upsi;}}_i=\frac{{\mathrm{\&Sigma;}}_{k=0}^{G-1}{\mathrm{\&mu;}}_{\mathrm{ik}}^m\mathrm{kH}\left(k\right)}{{\mathrm{\&Sigma;}}_{k=0}^{G-1}{\mathrm{\&mu;}}_{\mathrm{ik}}^{m}H\left(k\right)}i=\mathrm{1,2},...,C---\left(5\right)$

如果同前一聚类中心相比，满足条件υ^(n-1)-υ⁽ⁿ⁾＜ε（n表示计算次数、ε表示计算停止阈值），则停止计算，否则返回步骤b。

d、按照隶属度μ对图像进行多阈值划分，取得头皮层、颅骨层、脑实质区域的二值分割模板。

②对CT图像中颅骨内侧区域采用非局部均值修正的模糊C均值聚类分割算法分割出脑室区域，其方法是：

a、将颅骨内侧区域定义为待分割区域I，对I中第j个像数选取7×7领域窗口N_j，按下述公式（6）计算中心像素点同窗口内其它像素点灰度值间的局部约束系数ω_l；

${\mathrm{\&omega;}}_1(x_k,x_j)=e^{-\frac{{|x_k-x_j|}^2}{{\mathrm{\&sigma;}}^2}}---\left(6\right)$

式中，σ²为窗口内像素点灰度值的统计方差。

b、对I中第j个像数选取21×21搜索窗口S_j，按以下公式（7）计算中心像素点同窗口内其它像素点灰度值间的非局部约束系数ω_nl；

${\mathrm{\&omega;}}_{n1}(x_k,x_j)=\frac{1}{Z\left(x_j\right)}U(x_k,x_j)---\left(7\right)$

式中，U表示以像素点K与j为中心的两个领域窗间高斯核函数，Z表示以搜索窗内所有像素点作为K时函数U计算值的和。

c、按以下公式（8）、（9）分别计算局部距离函数

和非局部距离函数 $d_{n1}^2;$

$d_1^2(x_j,{\mathrm{\&upsi;}}_i)=\frac{{\mathrm{\&Sigma;}}_{k\&Element;N_j}{\mathrm{\&omega;}}_1(x_k,x_j)d^2(x_k,{\mathrm{\&upsi;}}_i)}{{\mathrm{\&Sigma;}}_{k\&Element;N_j}{\mathrm{\&omega;}}_1(x_k,x_j)}---\left(8\right)$

式中，ν为初始聚类中心。

$d_{n1}^2(x_k,{\mathrm{\&upsi;}}_i)=\underset{k\&Element;S_j}{\mathrm{\&Sigma;}}{\mathrm{\&omega;}}_{n1}(x_k,x_j)d^2(x_k,{\mathrm{\&upsi;}}_i),\left\{\begin{array}{c}0\&le;{\mathrm{\&omega;}}_{n1}(x_k,x_j)\&le;1\\ \underset{k\&Element;S_j}{\mathrm{\&Sigma;}}{\mathrm{\&omega;}}_{n1}(x_k,x_j)=1\end{array}\right.---\left(9\right)$

d、按以下公式（10）计算非局部均值修正FCM算法的距离函数D²；

$D^2(x_j,{\mathrm{\&upsi;}}_i)=(1-\mathrm{\&lambda;})d_1^2(x_j,{\mathrm{\&upsi;}}_i)+\mathrm{\&lambda;}{d}_{n1}^2(x_j,{\mathrm{\&upsi;}}_i)---\left(10\right)$

式中，λ为权重因子，取值（0，1）。

e、按以下公式（11）、（12）进行迭代运算，直到满足预设的终止条件为止；

${\mathrm{\&mu;}}_{\mathrm{ij}}={\left(\underset{k=1}{\overset{c}{\mathrm{\&Sigma;}}}{\left(\frac{D(x_j,{\mathrm{\&upsi;}}_i)}{D(x_j,{\mathrm{\&upsi;}}_k)}\right)}^{2/(m-1)}\right)}^{-1}---\left(11\right)$

${\mathrm{\&upsi;}}_i=\frac{{\mathrm{\&Sigma;}}_{j=1}^n{\mathrm{\&mu;}}_{\mathrm{ij}}^m{x}_j}{{\mathrm{\&Sigma;}}_{j=1}^n{\mathrm{\&mu;}}_{\mathrm{ij}}^m}---\left(12\right)$

式中，m表示模糊因子；

f、按照隶属度μ对图像进行多阈值划分，取得脑实质层和脑室层的二值分割模板。

③自动绘制脑脊液层边界，其方法为：将获取的脑实质层二值分割模板外侧进行二值图像腐蚀操作，利用7×7操作矩阵，将脑实质层外侧向内缩进4个像素点，同颅骨层内侧形成约2mm宽的间隔，将此间隔定义为脑脊液层二值分割模板。

④取得各组织层边界，其方法为：遍历各组织层二值分割模板选取任一边界点作为种子点，利用生长算法由种子点向八个方向搜索边界点，直到返回种子点为止，保存各组织层边界信息。

步骤3中，提取颅骨层密致骨和板障层的分布信息，其具体步骤如下：

根据步骤2的结果提取颅骨层的CT数据，按骨窗进行窗宽窗位变化。利用步骤2所述直方图约束的多阈值FCM算法，自动获取颅骨内、外板密质骨以及中间板障层的空间分布信息，建立颅骨密致骨层的二值分割模板。

步骤4中，将各层边界离散化，构建头部二维有限元模型并配置电阻率分布先验信息，其具体步骤如下：

①依据各层组织边界曲线的曲率定义边界采样点，其方法为：头皮外侧边界依据电极位置定义16个采样点，并按电极顺序编号，电极点间安装一定间隔距离平均定义采样点；其余各层组织边界按同一间隔距离平均定义采样点；遍历颅骨内侧、脑脊液层内侧、脑室外侧的边界点，将曲率不存在的点定义为采样点，即边界上一阶导数不存在且不趋于无穷的点。

②构建头部二维有限元模型，其方法为：

a、将边界离散化中获取的离散点定义为限定端点，各层采样点间的线段定义为限定线段，利用自适应限定Delaunay三角网格剖分算法进行初始剖分。具体实现中，首先采用标准Delaunay三角网格剖分算法根据限定点进行网格构建。因为标准方法不考虑限定条件，在得到的网格中部分限定线段会被算法删除，因此使用剔除三角单元的方法加入被删除的限定线段，恢复网格的限定条件。

b、利用三角单元径（外接圆半径）边（单元最短边长）比和面积作为网格质量和尺寸的控制因子，对初始剖分网格进行质量及尺寸可控的Delaunay三角网格优化，形成头部三角形网格。具体实现中，优化不满足质量和尺寸要求的三角单元，采用插入偏心点的方法进行。首先计算被优化三角单元的外接圆圆心，如果此圆心与被优化单元最短边所构成新三角形的径边比符合质量要求，则该点定义为偏心点；否则，在此圆心点与被优化单元最短边中点的连线上寻找一点，使该点与最短边构成新三角形的径边比符合质量要求，则该点定义为偏心点。将偏心点插入网格，按Delaunay准则更新网格，直到所有网格单元均符合质量要求和尺寸要求。

步骤5中，配置颅骨电阻率非均匀分布以及头部其它组织的电阻率分布信息，其具体步骤为：

①根据网格中各层组织边界的限定线段信息，对网格进行组织区域划分，标记网格的组织属性；

②配置颅骨电阻率非均匀分布，其方法为：

a、取颅骨层质心为原点，建立坐标系。

b、将颅骨二值分割模板绕原点进行逆时针旋转，间隔角度为1°。每旋转一次计算颅骨区域与X轴重叠的像素个数，这个值表示颅骨在此角度的像数厚度。同理，将密质骨二值分割模板绕原点逆时针旋转，间隔角度为1°。每旋转一次计算密质骨区域与X轴重叠的像数个数，表示密质骨在此角度的像数厚度。则颅骨PTD可由下式（13）计算：

由于骨缝结构的存在，会形成密致骨二值分割模板中内、外侧密致骨均不连续的区域，或单侧密致骨不连续区域，使密致骨在空间分布上不连续。对于骨缝区域的PTD均按80%计算，对应齿状骨缝的电阻率值（电阻率最低的颅骨组织结构）。

c、采用10°为间隔角度，从前述360个原始PTD采样值中抽取36个作为计算采样值。

为了保证骨缝电阻率信息的还原，每个计算采样值在其前后5个原始PTD采样值中搜索骨缝采样点，如果存在则用骨缝所在角度对应的原始PTD采样值作为该角度的计算采样值。这样就得到了间隔角度为10°的36个颅骨PTD计算采样值，形成36个电导率离散分布数据。

d、根据颅骨电阻率同PTD对应的统计参数模型计算颅骨电阻率值。

Y_{颅骨电阻率}＝25900.5-284.2·PTD           （14）

e、利用这36个离散电阻率分布值和对应的旋转角度为颅骨区域的三角形单元赋予电导率值；方法为：

遍历颅骨区域内的所有三角形单元，计算每个三角形单元的内切圆圆心，求出该点与颅骨质心连线与X轴的夹度∠S（范围为[0°,360°）），利用36个PTD计算采样值对应的离散电阻率分布值对其进行线性插值，得到角∠S对应的电导率值。

③按以下规则对其它组织层的网格电阻率属性赋值：头皮层：2.27Ω·m，脑脊液层：0.56Ω·m，脑实质层：6.67Ω·m，脑室层：0.56Ω·m，形成用于EIT动态成像计算的头部二维有限元模型。

所述步骤6中，由EIT系统采集被测对象头部实时阻抗测量数据，重构出融合头部先验信息的EIT图像并显示，其具体步骤如下：

①得到被测对象头部实时阻抗测量数据，其方法为：使用具有16电极和USB通讯接口的EIT测量系统采集激励频率50KHz、激励电流1mA时的边界电压数据，包括每一个电极注入激励电流时其它邻近电极间测量得到的电压。

②建立头部实时阻抗测量动态数据，其方法为：选取被测对象脑部生理或病理状态发生变化之前的边界电压数据作为前景数据V_f、选取状态发生之后的边界电压数据作为背景数据V_b，则EIT测量动态数据ΔV＝V_b-V_f，其中ΔV表示由于生理或病理变化造成的脑部组织阻抗改变所引起的边界电压的变化。当然，这里被测对象也可以是物理模型、仿真模型，状态变化则是物理模型内部出现扰动目标或仿真模型内部设置了电阻率变化区域。

③构建融合被测对象头部解剖结构和电阻率分布先验信息的EIT重构算法，其方法为：按公式（15）计算融合头部组织解剖结构和电阻率分布的雅克比矩阵，其表示边界电压变化同模型单元电阻率变化之间的关系；

$J=\frac{\&PartialD;U_{\mathrm{ij}}\left(\mathrm{\&rho;}\right)}{\&PartialD;{\mathrm{\&rho;}}_k}\left\{\begin{array}{c}i,j=\mathrm{1,2},...,N\\ k=\mathrm{1,2},...,M\end{array}\right.---\left(15\right)$

式中，M、N表示头部有限元模型的单元数和节点数，U_ij为第i个电极对激励下、第j个电极对的测量电压，P=（P₁，P₂，…，P_M）^T为单元的电阻率属性。

这样，在头部有限元模型中融入的被测对象头部解剖结构和组织电阻率分布信息形成的约束条件便转化成重构算法使用的雅克比矩阵。

基于标准Tikhonov正则化方法的最小二乘EIT重构算法，利用公式（16）重构，式中Δρ为重构计算得到的被测对象头部生理或病理状态发生变化引起的脑部组织电阻率变化量、R为J^TJ矩阵主对角线元素、λ为Tikhonov正则化参数，其中λ的选取可以使用L型曲线法。

Δρ＝[J^TJ+λR]^-1J^TΔV             （16）

④显示EIT图像，其方法为：将具有同一节点的单元电阻率的均值定义为该节点电阻率属性，在相邻两个节点间按一定步长进行电阻率属性的插值，根据插值结果利用云图方式显示。

以下是发明人给出的具体实施例。

参见附图1，图1给出了一种融入颅骨电阻率非均匀分布信息的电阻抗断层成像方法的原理框图。首先于被测对象头部眉弓以上2厘米、枕骨粗隆以上1厘米处的EIT电极分布层面粘贴CT成像标记物，其中CT成像标记物尺寸为粘贴面直径1厘米的圆柱形；由CT扫描设备行间隔小于5毫米的头部轴位扫描，扫描数据经CT工作站连入的PACS系统读入电阻抗断层成像软件系统，筛选CT成像标记物清晰的CT数据进行均值滤波预处理；根据选取的CT扫描数据构建融入颅骨电阻率非均匀分布信息的头部逼真有限元模型；在被测对象头部按标记位置粘贴测量电极并接入电阻抗测量系统，由软件系统配置激励和测量模式开始测量，将得到的数据通过USB接口传入计算机；根据建立的被测对象头部逼真模型和成像算法计算头部组织电阻率变化量，按云图方式绘制电阻率变化量的图像并显示输出。

参见附图2的图像分割聚类中心自动选取的流程图。为了实现自动分割、提高算法的运算速度，这里依据图像灰度直方图自动选取分类数和初始聚类中心。将预处理后的CT图像按脑部观察窗调整窗宽窗位，对非空气域的图像，即灰度值30～254的范围按公式⑴统计灰度直方图，其中趋近系数取0.1，由此取得的灰度直方图曲线过于杂乱难以分析，于是按公式⑵绘制图像灰度势函数直方图。根据曲线微分的零点，自动判断出代表头皮组织和颅骨组织的两个聚类中心。

参见附图3的提取头皮和颅骨组织分布边界信息的流程图。将前步选取的聚类中心代入公式（3）和（4），计算约束系数和隶属度。由公式（5）更新聚类中心值并按ε=0.001定义结束条件。如果更新前后的聚类中心满足结束条件则停止计算，否则再次迭代计算直到满足结束条件为止。计算得到的隶属度表示图像灰度值属于每个分类的比重，将具有相同灰度值的像素点划分到所占比重最大的分类中，属于该分类的像素点用1表示，否则用0表示，形成二值分割模板并提取边界。

参见附图4的提取脑脊液层和脑室组织分布边界信息的流程图。脑室组织分布的边界信息由非局部均值修正的模糊C均值聚类算法进行图像分割提取。图像分割区域局限在颅内内侧的脑组织区域，包括灰度较深的脑室区域和灰度较浅的脑实质区域。仍然依据直方图进行分类类别和初始聚类中心的自动选取，并在分割区域内按公式（6）和（7）计算局部和非局部约束函数。将结果代入公式（8）、（9）和（10）集散距离函数的值，其中权重因子λ=0.2。由公式（11）和（12）更新隶属度和聚类中心，其中模糊因子m=2，，这里的隶属度表示图像像素点属于每个分类的比重。仍然以ε=0.001为结束条件，计算完成后按隶属度分配像素点形成二值分割模板并提取边界。脑脊液层组织分布的边界信息由二值图像图形腐蚀操作提取。该实施例中选取被测对象的头部CT扫描数据的Pixel Spacing属性为0.5859mm，所以沿颅骨内侧边界向内缩进4个像素点，形成间隔2.3436mm的脑脊液层内侧边界。各层组织的边界组合形成被测对象头部组织分布的边界信息，如附图5所示。

参见附图6的构建头部逼真有限元模型的流程图。本实施例中按曲率变化关系确定了502个限定点和NN限定线段，其中曲率不存在的点位于颅骨内侧骨突的顶点与对应的脑脊液层上的点。由自适应限定Delaunay三角网格剖分及其优化算法得到的头部三角形网格包含1685单元和NN节点。

参见附图7的配置颅骨电阻率非均匀分布信息的流程图。标记单元组织属性后，以颅骨质心为原点建立坐标系（如图5所示）并按公式（13）计算颅骨的PTD采样值。根据统计参数模型（14）对颅骨层单元电阻率属性进行插值配置；根据文献记载的组织电阻率数据对其它组织层单元进行配置，并在电极节点定义注入激励电流的强度，形成被测对象的融入颅骨电阻率非均匀分布信息的头部逼真有限元模型，如附图8所示。

根据得到的有限元模型，按公式（15）计算雅克比矩阵，并使用附图9所示的物理模型生成EIT边界测量数据。按公式（16）计算头部组织电阻率变化量并按云图方式建立图像，如附图10所示。

以上对本发明及其实施例的描述，并不局限于此，因此发明并不限于具体实施方式中所述的实施例，凡是由本领域技术人员根据发明的计算方案得出的等效替换和增加，同样属于本发明保护的范围。

Claims (7)

1.一种融入颅骨电阻率非均匀分布信息的电阻抗断层成像方法，其特征在于，按如下步骤进行：

步骤1，获取被测对象头部电阻抗成像测量电极所在层面的轴位CT扫描数据；

步骤2，根据CT扫描数据自动划分头部的解剖结构，取得头部组织的分层结构信息；

步骤3，对颅骨层CT数据进行处理，获取密质骨和板障层分布信息；

步骤4，将头部组织的分层结构信息作为限定条件构建头部逼真有限元模型；

步骤5，配置颅骨电阻率非均匀分布以及头部其它组织的电阻率分布信息；

步骤6，由EIT测量系统采集头部实时阻抗测量数据，重构出融入颅骨电阻率非均匀分布信息的EIT图像并显示。

2.如权利要求1所述的方法，其特征在于，所述的获取测量电极所在层面的轴位CT扫描数据是通过在被测对象头部表面粘贴CT成像标记物确定头部轴位CT序列数据中EIT测量电极所在的层面。

3.如权利要求1所述的方法，其特征在于：所述的自动提取头部组织的分层结构信息，其具体步骤如下：

1）对CT图像采用直方图约束的自适应多阈值模糊C均值聚类分割算法提取颅骨层结构信息；

2）对CT图像中颅骨内侧区域采用非局部均值修正的自适应模糊C聚类分割算法提取脑室结构信息；

3）对颅骨内侧边界采用二值位图操作，自动绘制脑脊液层结构信息。

4.如权利要求1所述的方法，其特征在于，所述的获取颅骨层密质骨和板障层分布信息是指使用非局部均值修正的自适应模糊C聚类分割算法对颅骨层进行处理，提取密质骨、板障层和骨缝的结构分布信息。

5.如权利要求1所述的方法，其特征在于：所述的构建头部逼真有限元模型是指利用头部组织的分层信息、基于质量及尺寸可控的自适应限定Delaunay三角网格剖分技术构建用于头部EIT重构计算的有限元模型。

6.如权利要求1所述的方法，其特征在于，所述的配置头部组织电阻率分布信息，具体步骤包括：

①利用颅骨板障层厚度同颅骨电阻率值的统计参数模型，对颅骨层单元的电阻率属性进行赋值，实现颅骨电阻率非均匀分布信息的融入；

②根据人体头部组织电阻抗频谱特性测量结果，对头部其它组织层单元的电阻率属性进行赋值。

7.如权利要求1所述的方法，其特征在于，所述的由EIT系统采集被测对象头部实时阻抗测量数据，重构出融合头部先验信息的EIT图像并显示，具体方法为：利用头部EIT测量系统实时采集的数据和基于标准Tikhonov正则化方法的最小二乘EIT重构算法进行图像重构计算，将计算得到的电阻率分布利用云图方式显示。

Images