测量颅内椎-基底动脉血流动力学特性方法

摘要

本发明是一种颅内椎-基底动脉血流动力学分析方法。根据颅内椎-基底动脉血管的临床影像图像，建立解剖模型，本发明设计了颅内椎-基底动脉物理模型，并由此建立了相应的控制方程。通过方程，可以由测得颅内段椎-基底动脉段的几何参数、管腔直径和入口处的边界条件，求出椎-基底动脉的入口处的血流速度，求出各血管段的血流动力学特性参数。根据上述方法设计了由核磁共振技术和流体动力学分析软件组成的分析仪器。本发明对颅内椎-基底动脉血流动力学研究及临床指导都有十分重要意义。

说明书

技术领域

本发明属医疗技术领域，是一种测量椎-基底动脉血流动力学特性方法。 

背景技术

脑血管疾病是严重威胁人类健康的主要疾病之一，大量研究表明，脑血管疾病的发生、发展与脑循环系统血流动力学参数的异常改变密切相关。随着生物医学图像处理技术的飞速发展，以及大容量、高速度、多功能计算机软件和数值计算理论的不断完善，结合这些技术，利用计算流体力学的方法对血液血流动力学的参数进行数值模拟，研究脑血管系统的血流动力学机理，脑血管疾病发生与发展的过程，并对其结果进行处理分析。 

发明内容

本发明提供一种模拟颅内椎-基底动脉血流动力学特性方法，包括：选取正常人脑血管图像，利用核磁共振扫描技术测量得到整个脑血管图像，在磁共振ADW 4.2工作站中对整个脑血管图像进行处理，得到颅内椎-基底动脉血管成像。通过后处理技术精确获取椎-基底动脉详细的空间几何参数(即X，Y，Z坐标)，在左右椎动脉中对称选取垂直于血管截面的层面作为模型的入口边界条件，应用磁共振Fast PC CINE的方法测量此处的血流速度，采用计算机辅助设计(Computer Aided Design，CAD)的方法，通过磁共振所能提供的血管详细坐标以及入口处的边界条件，将数据记录后采用UG软件进行图像三维重建技术处理后，得到椎-基底动脉的物理模型，由此计算出此模型的血流动力学特性数据。 

参数设置 

采用0.8s为一个完整的心动周期。在0.8s的周期中，选取0.1s的时间作为初始时刻，左右椎动脉的流速经测算得到，通过数值公式计算，得到初始时刻0.1s时血液在血管中流动及壁面各种受力的数值结果。 

数值公式： 

控制方程 

对于血液在椎-基底动脉中流动及力学特性数值研究的基本思路是：首先血液在血管中流动的控制方程组，对控制方程组进行离散，由此获得差分方程并讨论方程的收敛性，结合血液在椎-基底动脉中流动的初始条件及边界条件，继而获得血液在椎-基底动脉中流动及力学特性。 

为了建立血液在椎-基底动脉中流动的数学模型，首先作如下假设： 

1)血液是均质各向同性不可压缩non-Newtonian流体； 

2)血液的流动满足Navier-Stokes方程； 

3)血管壁为刚性无滑移、不可渗透边界； 

4)除汇合及分叉处，血管为等截面圆管。 

基于以上的假设条件，血液在椎-基底动脉中流动过程可用下述控制方程组表示： 

连续性方程： 

$\frac{\partial \rho }{\partial t}+\mathrm{div}\left(\mathrm{\rho U}\right)=0$

能量方程： 

$\frac{\partial \left(\mathrm{\rho T}\right)}{\partial t}+\mathrm{div}\left(\mathrm{\rho UT}\right)=\mathrm{div}\left(\frac{\lambda }{c_p}\mathrm{grad}T\right)$

动量方程： 

$\frac{\partial \left(\mathrm{\rho u}\right)}{\partial t}+\mathrm{div}\left(\mathrm{\rho uU}\right)=\mathrm{div}\left(\eta \mathrm{grad}u\right)+S_u-\frac{\partial p}{\partial x}$

$\frac{\partial \left(\mathrm{\rho v}\right)}{\partial t}+\mathrm{div}\left(\mathrm{\rho vU}\right)=\mathrm{div}\left(\eta \mathrm{grad}v\right)+S_v-\frac{\partial p}{\partial y}$

$\frac{\partial \left(\mathrm{\rho w}\right)}{\partial t}+\mathrm{div}\left(\mathrm{\rho wU}\right)=\mathrm{div}\left(\eta \mathrm{grad}w\right)+S_w-\frac{\partial p}{\partial z}$

边界条件： 

T_f|_x＝0＝T_w，v_f|_x＝0＝v_w

初始条件： 

T_s|_t＝0＝T₀，u＝v＝w＝0 

湍流模型 

血管内的流动为层流和湍流混合并存，因此基于k-ω湍流模型的SST(The Shear Stress Transport)模型来描述层流与湍流并存的流动。 

具体计算公式如下： 

其中μ₀＝0.0035Pa，n＝0.6， 为剪切率。 

控制方程的离散 

采用了一种基于有限元的有限体积法，即CFX方法，完成方程的离散化。 

将血液在血管中流动的质量方程、动量方程、能量方程及反应方程进行离散，获得差分方程，结合血液在血管中流动过程的实际边界值条件进行计算。由于计算几何模型较复杂，故划分为非结构化网格，考虑到壁面边界层效应，近壁处的网格比较密集。网格的划分采用ICEM-CFD-12.0进行。 

上文所述的CFX方法的主要内容为：根据网格划分控制容积，并在一个控制容积内，应用高斯度定理对微分方程进行积分，使微分方程完成初步的离散，变为在各个积分点的 值以及 的导数值的方程，通过引入形状函数，使各积分点上 值以及 对应的梯度值用节点上的 值表示出来，从而完成方程的离散化。对于每个节点控制体积的选定，是将每个三角形单元的形心与其三个边的中点相连，因此围绕每一个内节点就会形成一个多边形，这种多边形就是这个节点的控制容积，如图所示：(图6、7)。 

在笛卡尔坐标下的通用微分方程为： 

$\frac{\partial }{\partial t}\left(\mathrm{\rho \phi }\right)+\frac{\partial }{{\partial x}_j}\left({\mathrm{\rho U}}_j\phi \right)=\frac{\partial }{{\partial x}_j}\left({\Gamma }_{\mathrm{eff}}\left(\frac{\partial \phi }{{\partial x}_j}\right)\right)+S_{\phi }$

对微分方程进行高斯积分： 

$\frac{d}{\mathrm{dt}}\underset{V}{\int }\mathrm{\rho \phi dV}+\underset{S}{\int }\rho U_j\phi {\mathrm{dn}}_j=\underset{S}{\int }{\Gamma }_{\mathrm{eff}}\left(\frac{\partial \phi }{{\partial x}_j}\right){\mathrm{dn}}_j+\underset{V}{\int }{S}_{\phi }\mathrm{dV}$

整理得： 

$\mathrm{\rho V}\left(\frac{\phi -{\phi }^0}{\mathrm{\Delta t}}\right)+\underset{\mathrm{ip}}{\Sigma }{\stackrel{·}{m}}_{\mathrm{ip}}{\phi }_{\mathrm{ip}}=\underset{\mathrm{ip}}{\Sigma }{\left({\Gamma }_{\mathrm{eff}}\frac{\partial \phi }{{\partial x}_j}{\mathrm{\Delta n}}_j\right)}_{\mathrm{ip}}+{\bar{S}}_{\phi }V$

从以上公式中可以看出求解方程需要知道各积分点上φ值以及φ的对应梯度值，这里引入有限元方法中 常用的形状函数，在CFX中均采用线性的形状函数(Ni)。 

$\phi =\underset{i=1}{\overset{N_{\mathrm{node}}}{\Sigma }}{N}_i{\phi }_i$

从而扩散项中的φ_ip的梯度项可以表示为： 

$\frac{\partial \phi }{\partial x}{|}_{\mathrm{ip}}=\underset{n}{\Sigma }\frac{{\partial N}_n}{\partial x}{|}_{\mathrm{ip}}{\phi }_n$

对于对流项φ_ip的处理，则采用如下的一个通用格式 

其中： 

φ_up为迎风节点的φ值；β可视为一个通用因子，通过定义不同的β值，完成不同的差分格式；▽φ为因变量φ的梯度值； 为从迎风节点到积分点的方向矢量。由以上方式的处理，各积分点上 值以及 的对应梯度值就可以用节点上的 值表示出来，从而完成方程的离散化。 

本发明的优点为，通过磁共振提供的血管详细坐标以及入口处的边界条件，将数据记录处理后得到椎-基底动脉的物理模型，由此模型进一步计算出血流动力学特性数据。 

附图说明

图1椎-基底动脉磁共振血管图像，为人体真实血管图像 

图2椎-基底动脉物理模型，是采用计算机辅助设计(Computer Aided Design，CAD)的方法，通过获取椎-基底动脉的几何参数(X，Y，Z坐标)，利用数字图像处理的三维重建技术得到解剖真实的人体椎-基底动脉的几何模型 

图3椎-基底动脉血液流线图 

图4椎-基底动脉血管外壁压力分布图 

图5椎-基底动脉血管外壁剪切力分布图 

图3-5都是表示椎-基底动脉血流动力学结果的图形 

图6有限控制容积的划分 

图7单元积分点分布示意图 

具体实施方式

1模型原始图像获取参数 

常规MRI 

采用美国GE公司Signa Excite HD 3.0T高场强MR扫描仪，应用标准头部8通道线圈。检查前告知受检者磁共振检查的注意事项，给受检者佩戴耳塞用来减小噪音。摆正受检者体位。以前后连合连线为轴位图像扫描基线。MRI研究包括轴位及矢状位T1-FLAIR、轴位T2-FSE序列。扫描参数：轴位T1-FLAIR：TR 3196ms，TE 7.9ms，TI 960ms；矢状位T1-FLAIR：TR 1750ms，TE 8.9ms，TI 920ms；轴位T2-FSE：PROP，TR5100ms，TE118ms；激励次数均为1次，层厚5.0mm，间隔1.0mm，矩阵512×256，视野240mm×240mm。 

MRA 

采用3D TOF法，扫描参数：TR 20ms，TE 3.2ms，层厚1.0mm。 

Fast Cine PC 

用3D TOF MRA血管图像进行定位，扫描层面与颅内段左右椎动脉入口处垂直。采用无相位卷折技术、呼吸补偿及流动补偿，使用外周脉搏门控，以便获得单位时间内的相位改变曲线(外周门控的获得以受检者的右手中指指腹采集而得)，设定30个时相。流动编码方向为Slice，与血液流动方向相反。扫描参数：TR/TE＝22/10ms，层厚(slice thickness)＝4mm，间隔(spacing)＝0mm，视野＝14×14，矩阵＝256×128，激励次数(NEX)＝1，翻转角(flip angle)＝20°，带宽(bandwidth)＝31.25。扫描时间20s左右(扫描时间随心率不同而有所差异)。 

2图像后处理及分析 

3D TOF MRA 

观察颅脑MRI及MRA图像，在ADW 4.2工作站用最大密度投影(maximum intensity projection，MIP)对MRA源图像进行重组，分别得到前后循环的图像。记录血管形态，测量与Fast Cine PC相对应的左右椎动脉管腔直径。 

Fast Cine PC 

Fast Cine PC序列可以得到相位图和幅度图两种图像。在ADW 4.2工作站对图像进行后处理。将图像放大4倍，窗宽窗位调至视野显示最清晰，在量值图上勾画出血管的截面，以3D TOF MRA MIP重组后图像测得的直径为参考，使所测管径横截面积与实际管腔面积最接近，再将此感兴趣区复制到相位图上，在相位图上，即可得到感兴趣区一个心动周期内30个时相的流速。测定感兴趣区的横截面积，进而得到管腔流速及单位时间内通过管腔的流量。 

3计算机辅助设计的多参数配置 

系统组成 

通常以具有图形功能的交互计算机系统为基础，主要设备有：计算机主机，图形显示终端，图形输入板，绘图仪，扫描仪，打印机，磁带机，以及各类软件。 

工程工作站 

一般指具有超级小型机功能和三维图形处理能力的一种单用户交互式计算机系统。它有较强的计算能力，用规范的图形软件，有高分辨率的显示终端，可以联在资源共享的局域网上工作，已形成最流行的cad系统。 

图形输入输出设备 

除了计算机主机和一般的外围设备外，计算机辅助设计主要使用图形输入输出设备。交互图形系统对cad尤为重要。图形输入设备的一般作用是把平面上点的坐标送入计算机。常见的输入设备有键盘、光笔、触摸屏、操纵杆、跟踪球、鼠标器、图形输入板和数字化仪。图形输出设备分为软拷贝和硬拷贝两大类。软拷贝设备指各种图形显示设备，是人机交互必不可少的；硬拷贝设备常用作图形显示的附属设备，它把屏幕上的图像复印出来，以便保存。常用的图形显示有三种：有向束显示、存储管显示和光栅扫描显示。有向束显示应用最早，为了使图像清晰，电子束必须不断重画图形，故又称刷新显示，它易于擦除和修改图形，适于作交互图形的手段。存储管显示保存图像而不必刷新，故能显示大量数据，且价格较低。光栅扫描系统能提供彩色图像，图像信息可存放在所谓帧缓冲存储器里，图像的分辨率较高。 

CAD软件 

除计算机本身的软件如操作系统、编译程序外，cad主要使用交互式图形显示软件、cad应用软件和数据管理软件3类软件。 

交互式图形显示软件用于图形显示的开窗、剪辑、观看，图形的变换、修改，以及相应的人机交互。cad应用软件提供几何造型、特征计算、绘图等功能，以完成面向各专业领域的各种专门设计。构造应用软件的四个要素是：算法、数据结构、用户界面和数据管理。数据管理软件用于存储、检索和处理大量数据，包括文字和图形信息。 

本发明通过CAD软件通过获取椎-基底动脉的几何参数(X，Y，Z坐标)，以及入口处椎动脉的流速，利用数字图像处理技术得到椎-基底动脉的几何模型。 

基本技术： 

主要包括交互技术、图形变换技术、曲面造型和实体造型技术等。 

在计算机辅助设计中，交互技术是必不可少的。交互式cad系统，指用户在使用计算机系统进行设计时，人和机器可以及时地交换信息。采用交互式系统，人们可以边构思、边打样、边修改，随时可从图形终端屏幕上看到每一步操作的显示结果，非常直观。 

图形变换的主要功能是把用户坐标系和图形输出设备的坐标系联系起来；对图形作平移、旋转、缩放、透视变换；通过矩阵运算来实现图形变换。 

4模型的边界条件参数设置 

1)血液是均质各向同性不可压缩non-Newtonian流体； 

2)血液的流动满足Navier-Stokes方程； 

3)血管壁为刚性无滑移、不可渗透边界； 

4)除汇合及分叉处，血管为等截面圆管。 

5)两条椎动脉入口处给定血液的速度vg0、va0和温度Tg0、Ta0，血液的速度采用正常形态血管椎动脉进口速度的测试值 

6)出口给定压力条件，出口相对压力设置参考人体正常脑部血压取pout＝10526Pa。 

7)血管外壁外面定温度条件，参考正常人体温度值血管外壁外面温度取37.5℃，初始的速度场u＝v＝w＝0，初始和血液入口温度T0＝36.5℃，初始压力场取prel＝10526Pa。 

8)血管壁面采用非滑移界面。对于收敛残差的设置，这里取均方根残差RMS＝10-4。