一种基于脑血流动力学的自动化脑动脉瘤破裂风险分析系统的制作方法

1.本发明属医学图像处理及应用领域，涉及医学影像后续处理分析技术，具体涉及一种基于脑血流动力学的自动化脑动脉瘤破裂风险分析系统，尤其是一种三维血流动力学自动化计算流程，其计算程序包括模型预处理模块、示踪颗粒模块、径迹统计模块、骨架提取模块、网格化分模块以及流体求解模块设计可自动执行输出血流场量结果的程序。可完成血流仿真计算，判断动脉瘤破裂风险，方便用物理仿真技术辅助手术方案，有助于临床应用。


背景技术：

2.研究报道，颅内动脉瘤发生的确切原因尚不十分清楚，有研究表明血流动力学(hemodynamics)、遗传学以及血管生物力学等多种因素均在颅内动脉瘤的发生中起重要作用。目前研究通常认为，血流动力学因素在动脉瘤的发生、发展以及破裂过程中起重要作用，从而使得调整动脉瘤的血流动力学成为颅内动脉瘤治疗的目标之一。
3.现有技术的血流动力学实验研究可以分为两类：动物在体实验研究和体外模型实验研究。研究实践显示，在体实验能够得到最为真实的数据，但体内实验存在周期较长，实验结果受个体差异影响大，长时间的控制和监测血流动力学参数极为困难等局限性；体外模型试验在某种程度上弥补了某些局限比如用彩色多普勒检测技术和磁共振血管成像技术等，但体外实验尚存在比如模型实验与实际流动条件很难做到近似，且从模型实验中所提取到的流动参数极其常规和有限等局限；因此，血流动力学仿真显得很有必要，其可在一定程度地弥补实验本身的不足。
4.cfd是通过计算机数值计算和图像显示，对包含有流体流动和热传导等相关物理现象所做的系统分析，其基本思想是将原来在时间域和空间域上联系的物理量的场，如速度场和压力场，用一系列有限个离散点上的变量值的集合来代替，通过一定的原则和方式建立起关于这些离散点上变量之间关系的代数方程组，然后求解代数方程组获得场变量的近似值。cfd可以看做是流动基本方程控制下对流动的数值模拟，通过这种数值模拟得到极其复杂问题的流场内各个位置上的基本物理量(如速度、压力、温度、浓度等)的分布，以及这些物理量随时间的变化情况。目前利用cfd软件仿真血流动力学主要分为以下流程：首先是利用医学图像处理软件例如mimics，将由临床采集的dicom格式图像转化，重建患者血管三维几何(stl格式)保存；然后把三维血管几何导入计算流体网格划分软件如ansys—icem生成数值模拟用计算网格；然后针对模型入口和出口需要设置边界条件，如：血压、血液流速、流量等；设置血液和血管壁的属性参数，如：密度、粘度、弹性等；最后开始计算到收敛后提取可视化参数；数值模拟研究的可靠性取决于上述步骤是否完全正确进行，欲获得可靠的数值模拟研究结果，研究者必须对每一步骤严格质量控制，争取每一步所得参数和患者实际情况一样。
5.临床实践中，对于曾破裂过的形态简单的颅内动脉瘤，可直接选择相应的手术治疗，但对于某些曾破裂过的复杂动脉瘤(包括多发动脉瘤、畸形动脉瘤)，仿真分析对于临床
诊疗可提供很大帮助，有助于医生更深入认识动脉瘤的血液动力学特征，帮助判断动脉瘤壁的最大受力位置和可能破裂的位置，医生参考其计算结果，结合附近血管的解剖结构，可选择合适的手术方式(介入栓塞、开颅动脉瘤夹闭)；如：如果数值模拟预测出囊状动脉瘤瘤底von mises应力明显增大，选择介入栓塞可能增加弹簧栓刺破动脉瘤壁导致颅内出血的风险，该种病例最好选择开颅夹闭术；当然，通过数值模拟研究结果选择何种手术方式为最佳以往研究尚未形成定论，需要进一步开展大量颅内动脉瘤的数值模拟研究，以及数值模拟研究的结果与手术方式选择和预后的关系的研究。
6.颅内动脉瘤数值计算属交叉学科的研究，血流计算基于cfd程序，血管壁的力学分析基于cae程序，整个分析过程涉及流体力学、数学、计算机软件、计算机仿真、有限元法、材料学、拓扑学、医学等，使用商业软件分析的流程也比较偏力学专业，目前较少有完全掌握上述研究方法的人员，所述方法尚未在临床工作中大量推广。仿真过程中最为重要的一环是区域离散过程，即需要对流体区域或者固体区域划分计算网格然后才能进行计算求解，网格质量直接决定了结果准确度。如果从stl模型出发，将网格制作过程和计算过程打包成一个和使用者交互最少的“黑盒子”，那么整个仿真过程则变为影像提取，模型制作，边界定义这三阶段，并且这三个阶段的设置均属于临床医生的专业范畴，仿真门槛将大为降低。
7.基于现有技术的现状，本技术的发明人拟提供一种基于脑血流动力学的自动化脑动脉瘤破裂风险分析系统，尤其是一种三维血流动力学自动化计算流程，所述分析系统可完成血流仿真计算，判断动脉瘤破裂风险，方便用物理仿真技术辅助手术方案，有助于临床应用。
8.与本发明有关的参考文献有：
9.[1]updegrove,a.,wilson,n.,merkow,j.,lan,h.,marsden,a.l.and shadden,s.c.,simvascular
‑
an open source pipeline for cardiovascular simulation,annals of biomedical engineering(2016).doi:10.1007/s10439
‑
016
‑
1762
‑8[0010]
[2]刘峰瑞,黄建平,黎忠,王纳秀等.圆柱容器内生成单粒径颗粒堆积的两阶段重力落球算法[j].核技术,2017。


技术实现要素：

[0011]
本发明的目的是基于现有技术的现状，提供一种基于脑血流动力学的自动化脑动脉瘤破裂风险分析系统，尤其是一种三维血流动力学自动化计算流程，本发明的风险分析系统可自动执行输出血流场量结果的程序，完成血流仿真计算，判断动脉瘤破裂风险，方便用物理仿真技术辅助手术方案，有助于临床应用。
[0012]
具体的，本发明提供了一种基于脑血流动力学的自动化脑动脉瘤破裂风险分析系统，其包括计算三维血流动力学的自动化流程，其特征在于，所述自动化流程中的血流自动化计算程序包括，模型预处理模块、示踪颗粒模块、径迹统计模块、骨架提取模块、网格化分模块以及流体求解模块；
[0013]
所述自动化流程只需使用者提供血管stl模型，并指定模型中血流的进出口面和相应的边界条件，其程序可自动执行输出血流场量结果。
[0014]
本发明中，所述的示踪颗粒模块中基于颗粒动力学的方法提取血管模型的拓扑学信息，避免传统基于拓扑学的骨架提取算法的算法不稳定的问题。
[0015]
本发明中，所述的径迹统计模块中基于蒙特卡洛思想对示踪颗粒进行注入，颗粒的初始位置在端面上随机分布，颗粒的初始速度也呈随机分布，不断重复注入过程增强颗粒对血管几何的捕获能力，当径迹类别稳定后结束注入过程。
[0016]
本发明中，所述的骨架提取模块中对径迹统计得到的初始骨架进行最大内切球“矫正”，矫正后的骨架提供血管瘤位置和体积信息。
[0017]
本发明中，所述的网格划分模块中用分割思想，将模型依据骨架分割为子模型，对每个子模型单独进行网格划分，子网格之间不存在拓扑约束。
[0018]
本发明中，所述的流体求解模块中的流体求解器针对不同的血流区域做耦合求解。
[0019]
本发明中，在和临床医生少量交互的情况下自动计算出血流的动力学参数，需要临床医生依据dicom影像提供stl格式的血管模型，程序对stl进行重新三角化，在模型文件外形不变的前提下尽可能得使得三角面片分布均匀，便于后续网格划分和颗粒
‑
血管壁的碰撞检测；对于模型中尖锐的区域，程序对模型做适当平滑接近血管真实的特征；基于vtk库提供临床医生一个模型交互界面，用于指定血管簇的进口和出口及相关的血流流动边界条件，这也是本发明中程序唯一需要和用户交互的环节。基于目前采用的是切割命名操作，需临床医生在进出口定性地做切面标记进出口位置，对于进口边界，临床医生可以在程序输入卡上指定进口面沿着外法线方向延长一段距离，为了保证流体流动充分发展，在血管簇的每个进出口端面，程序开启示踪颗粒注入操作，按照蒙特卡洛思想注入不同批次不同速度的颗粒群，颗粒轨迹稳定后对示踪颗粒轨迹做统计分析，标记血管分支和内部存在涡旋轨迹的区域，然后进一步基于示踪颗粒膨胀算法对统计得到的血管骨架做最大内切球矫正，矫正后的骨架和最大内切球半径分布可以供临床使用；骨架信息得到后，基于涡旋区域对模型做分割，对分割后的模型集合单独做带边界层的网格剖分，对于剖分不成功的区域，退化为不带边界层的网格，但作为网格质量下降的折衷处理方式，无边界层的网格尺寸设置偏小；网格划分需要用户提前提供进口血流条件，从而基于血流雷诺数来决定网格边界层的厚度等网格尺寸参数；最后综合临床医生指定的边界条件和已生成的分区网格进行流体计算，其具体执行需要调用针对多区网格的专用流体求解器；
[0020]
本发明中整个程序执行流程如图3所示；所述的核心算法分为四个模块：统计示踪颗粒轨迹/血管骨架提取/模型网格划分和流体求解器四大部分。
[0021]
本发明中，对于给定的stl模型，首先程序自动识别出血管簇所有的首末端面，将后面待处理的端面移除使得模型在端面出呈开放状态(后续注入示踪颗粒)，然后程序遍历每个血管簇首末端面，利用蒙特卡洛的思想，从端面处随机抛出一些球形颗粒，让这些颗粒在随机初始速度的作用下沿着端面入口向前运动，设置整个体系的外力场来驱动颗粒运动(作用类似流体流动方程中的压差)。此时这个体系仅存在颗粒
‑
颗粒间作用力和颗粒
‑
三角单元间的作用力(即颗粒
‑
血管壁作用)，追踪颗粒轨迹，用统计的方法可以归纳出颗粒群的拓扑特性，为后续骨架提取和网格划分提供基础信息。
[0022]
本发明提供了一种基于脑血流动力学的自动化脑动脉瘤破裂风险分析系统，尤其是一种三维血流动力学自动化计算流程，其计算程序包括模型预处理模块、示踪颗粒模块、径迹统计模块、骨架提取模块、网格化分模块以及流体求解模块设计可自动执行输出血流场量结果的程序。其主程序基于python语言，复杂功能模块基于c 封装为动态库供主程序
调用。所述流程只需使用者提供血管stl模型，并指定模型中血流的进出口面和相应的边界条件，即可自动执行输出血流场量结果。本发明通过与临床医生尽可能少的交互可完成血流仿真计算，判断动脉瘤破裂风险分别为，极高危，高危，中危，低危组，方便用物理仿真技术辅助手术方案，有助于临床应用。本发明极大地降低了血流仿真的技术门槛，可方便临床医生在术前仿真观察血流细节，为手术方案提供参考信息。
附图说明
[0023]
图1是离散单元法的执行流程。
[0024]
图2是流体计算的流程图。
[0025]
图3是主程序整体的执行流程。
具体实施方式
[0026]
实施例1
[0027]
计算颗粒在外力作用下的运动采取颗粒离散元法(dem，流程如图1所示)，将颗粒之间和颗粒
‑
面片的受力等效为弹簧
‑
阻尼模型，弹簧和阻尼分别在法向n和切向t施加；对于两个接触颗粒i和j，半径分别为ri和rj，颗粒i速度角速度的控制方程为：
[0028][0029][0030]
其中mi，ii分别代表颗粒i的质量和转动惯量，nc为与颗粒i相接触的邻居颗粒数目，程序采用简化的hertz
–
mindlin
–
deresiewicz接触模型计算颗粒间的接触力：
[0031][0032][0033]
其中k和γ代表硬度和阻尼系数；代表颗粒间的相对速度；μ代表摩擦系数。颗粒间的切向形变和法向形变分别为：
[0034][0035][0036]
切向和法向上的刚度系数k和阻尼系数γ分别为：
[0037][0038][0039]
[0040]
γ
n
＝γ
t
[0041]
其中r*,m*和e*分别为：
[0042][0043][0044][0045]
e,g和υ分别是杨氏模量,剪切模量与泊松比，ε是恢复系数，
[0046]
最大的时间步长由瑞利时间步长决定：
[0047][0048][0049]
颗粒的位置、速度等随时间演化，采用蛙跳算法(leap
‑
frog algorithm)迭代更新，该算法三阶精度且时间可逆，使用当前时间步的位置、前半个时间步的速度与当前时刻颗粒的受力进行位置与速度的更新，
[0050][0051][0052]
其中分别表示t时刻颗粒的速度，位置和受力，
[0053]
dem算法采用邻居列表的方式来搜索颗粒邻居，完成邻居搜索后紧接着进行严格的碰撞检测，依据颗粒和其邻居之间的相对位置来施加颗粒
‑
颗粒和颗粒
‑
血管壁相互作用力，最后利用牛顿定律，更新每个颗粒的位移，速度和角速度等动力学信息，流程如图1所示。
[0054]
示踪颗粒轨迹第一个特点是针对不同进出口的连通特性，类似微分方程中场量的散度定义，多条轨迹可能会在某一段区域交汇，也可能会在某个点分离，可以将某个截面上的颗粒质量流量来作为拓扑学中线条分与合的标志，质量流量增加的地方代表血流分支的合并点，反之亦然。示踪颗粒从某个端面进入血管网络，在颗粒
‑
颗粒和血管壁之间的作用下向前运动直至从某个端面流出血管网络，期间一条完整的颗粒轨迹可以表征一支血管通路，遍历所有端面即可得到血管通路集合。将血管通路合并可以形成一个初步的血管模型骨架；
[0055]
模型骨架应该能够以简洁的方式体现模型整体的拓扑结构,同时还能很好地保持模型的分枝结构。为了提高模型骨架提取的质量,需要保证模型骨架点能很好地集中在模型中轴位置。统计示踪颗粒轨迹得到的通路集合取决于蒙卡过程中的示踪方式，本身并不一定临近血管中心，所以需要基于三维最大内切球原则对其做一次“矫正”；具体方法是在
这些通路上间隔放置一定半径的虚拟颗粒，让颗粒半径逐步增加，将颗粒
‑
颗粒作用取消，单纯考虑颗粒
‑
壁面相互作用力，增加颗粒半径直至在血管壁作用下移动并达到稳定时即可得到矫正后的血管通路，此时血管通路上虚拟颗粒位置处最大内切球的半径信息均已得到，血管系统的拓扑信息可供临床应用，该算法思想类似“圆柱容器内生成单粒径颗粒堆积的两阶段重力落球算法”中的颗粒膨胀算法，通过颗粒动力学算法来避过纯拓扑学算法上的一些局限[2]，增加了算法的鲁棒性更适合嵌入到自动化流程中；
[0056]
示踪颗粒轨迹另一个特点是局域范围内的涡旋特性，可以看作血管内部入口出口重叠情况下的连通特性(如血管瘤内部的血流涡旋)，类似场量的旋度定义，示踪颗粒有时会从血管壁的某个位置流出然后再流入，其轨迹是闭合曲线，此时血流通量散度为0但旋度不为0，在涡旋的内部示踪颗粒质量流量非0，在上下游示踪颗粒质量流量不变，无法判断涡旋的存在，但此时经过矫正后的血管骨架在涡旋区域会退化为相对很短的分支线段，所以对这些非用户标记的端面，可以在分支线段代表的区域内进行示踪颗粒质量流量检测，从而判断出涡旋区域。
[0057]
实施例2
[0058]
基于涡旋区域往往是临床所关心的区域，通常也是影响网格质量的关键位置，传统对血管模型整体做网格划分，很难对所有区域都划分得当，经常需要反复折衷尝试多次才能得到一个比较完整的整体网格，需要大量的人机交互；本发明中程序按照涡旋区域将血管的stl模型进行分割，分割后的每个子集进行独立网格剖分，不同血管单元之间不要求网格一致性，由此，每个网格单元拓扑形状和血流特性相对单一，可以对每一个血管子集生成符合其流动特性的流体网格，是整个仿真流程自动化实现的关键；
[0059]
对血管子集生成网格子集后，流体求解器相应地需要支持不同子血管交接面的信息传递(传统的cfd典型算法piso/simple只针对单个流动区域)，本发明程序中流体求解器采用piso算法，分别对每个子网格中的流动进行迭代计算，同时在不同子血管交接面进行速度压力等场量的插值；基于计算流体力学(cfd)计算流体动力学信息，遵循的控制方程主要包括连续性方程：
[0060][0061]
对于密度为常数的流动上式简化为：
[0062][0063]
动量守恒方程：
[0064][0065]
左边第一项代表微团速度的变化率；左边第二项代表对流项；右边第一项代表压力梯度项；右边第二项代表作用在微团上的外部体积力；右边第三项代表微团应力；最后一项代表重力；对于各向同性流体，由流体本够关系有：
[0066]
[0067]
程序采取基于雷诺平均n
‑
s方程的湍流模型：
[0068][0069][0070]
在计算开始前还需要临床医生指定血管边界条件，这也是算法唯一需要和使用者交互的地方，使用者需要移动界面移动/缩放/旋转截面，将此矩形截面与stl模型做切割得到的切面作为血流进出口面，然后在保存窗口输入面类型(进口/出口),进口边界(血流流量/流速分布模型/拉伸长度)，出口边界(出口压力)等流动信息；流体求解器基于用户指定的边界和物性参数进行计算,对于每个子网个区域，首先给定的初始压力或上一迭代步压力求解动量方程，因为求得的速度变量并不一定满足连续性方程所以需要求解压力泊松方程求解得到压力，利用求得的压力修正速度，使之能够满足连续性方程。如果速度不满足动量方程则重复该循环，直至收敛为止，程序流程如图2所示。
[0071]
本发明可完成血流仿真计算，判断动脉瘤破裂风险分别为，极高危，高危，中危，低危组，方便用物理仿真技术辅助手术方案，有助于临床应用。