CN116453697B - 基于ffr拟合的冠脉狭窄血流动力学模拟方法及系统 - Google Patents

基于ffr拟合的冠脉狭窄血流动力学模拟方法及系统 Download PDF

Info

Publication number
CN116453697B
CN116453697B CN202211732598.4A CN202211732598A CN116453697B CN 116453697 B CN116453697 B CN 116453697B CN 202211732598 A CN202211732598 A CN 202211732598A CN 116453697 B CN116453697 B CN 116453697B
Authority
CN
China
Prior art keywords
pressure
coronary
coronary artery
stenosis
ffr
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202211732598.4A
Other languages
English (en)
Other versions
CN116453697A (zh
Inventor
赵英红
吴美静
陆晔
王慧慧
唐慧
唐璐
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xuzhou Medical University
Original Assignee
Xuzhou Medical University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Xuzhou Medical University filed Critical Xuzhou Medical University
Priority to CN202211732598.4A priority Critical patent/CN116453697B/zh
Publication of CN116453697A publication Critical patent/CN116453697A/zh
Application granted granted Critical
Publication of CN116453697B publication Critical patent/CN116453697B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G16INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16HHEALTHCARE INFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR THE HANDLING OR PROCESSING OF MEDICAL OR HEALTHCARE DATA
    • G16H50/00ICT specially adapted for medical diagnosis, medical simulation or medical data mining; ICT specially adapted for detecting, monitoring or modelling epidemics or pandemics
    • G16H50/50ICT specially adapted for medical diagnosis, medical simulation or medical data mining; ICT specially adapted for detecting, monitoring or modelling epidemics or pandemics for simulation or modelling of medical disorders
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • G06F30/23Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • G06F30/28Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/20Finite element generation, e.g. wire-frame surface description, tesselation
    • G06T17/205Re-meshing
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T19/00Manipulating 3D models or images for computer graphics
    • G06T19/20Editing of 3D images, e.g. changing shapes or colours, aligning objects or positioning parts
    • GPHYSICS
    • G16INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16HHEALTHCARE INFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR THE HANDLING OR PROCESSING OF MEDICAL OR HEALTHCARE DATA
    • G16H30/00ICT specially adapted for the handling or processing of medical images
    • G16H30/40ICT specially adapted for the handling or processing of medical images for processing medical images, e.g. editing
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00Details relating to CAD techniques
    • G06F2111/10Numerical modelling
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00Details relating to the application field
    • G06F2113/08Fluids
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14Force analysis or force optimisation, e.g. static or dynamic forces
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2210/00Indexing scheme for image generation or computer graphics
    • G06T2210/24Fluid dynamics
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2210/00Indexing scheme for image generation or computer graphics
    • G06T2210/41Medical
    • YGENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02TCLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Health & Medical Sciences (AREA)
  • Geometry (AREA)
  • Computer Hardware Design (AREA)
  • Medical Informatics (AREA)
  • General Engineering & Computer Science (AREA)
  • Public Health (AREA)
  • Epidemiology (AREA)
  • Primary Health Care (AREA)
  • Evolutionary Computation (AREA)
  • General Health & Medical Sciences (AREA)
  • Computer Graphics (AREA)
  • Software Systems (AREA)
  • Algebra (AREA)
  • Fluid Mechanics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computing Systems (AREA)
  • Radiology & Medical Imaging (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Pathology (AREA)
  • Databases & Information Systems (AREA)
  • Data Mining & Analysis (AREA)
  • Biomedical Technology (AREA)
  • Architecture (AREA)
  • Measuring Pulse, Heart Rate, Blood Pressure Or Blood Flow (AREA)

Abstract

本发明公开一种基于FFR拟合的冠脉狭窄血流动力学模拟方法,包括以下步骤:S1,建立冠脉狭窄个体化三维模型,获得网格化平滑冠脉狭窄三维模型;S2,对FFR实测数据进行拟合,量化心动周期内的时序血流压力输入;时序血流压力包括冠脉狭窄远端压力和主动脉压力;S3,添加时间轴进行血流模拟计算,量化计算血管内任意点及心动周期内任意时刻的血流动力学参数;采用不同狭窄病变程度的冠状动脉狭窄患者的影像数据建立个体化冠脉模型;融合临床FFR实测数据,对冠脉模型进行仿真计算,量化血流动力学实时参数分布,并对血流变异性提供实时的定量描述,分析血管的内皮细胞在不同的流体力学环境中得到的适应性反应。

Description

基于FFR拟合的冠脉狭窄血流动力学模拟方法及系统
技术领域
本发明公开一种基于FFR拟合的冠脉狭窄血流动力学模拟方法及系统,属于血流动力学数值模拟领域。
背景技术
现有技术对心脑血管内血液流动进行血流动力学数值分析,主要集中于冠脉狭窄血管的内皮细胞在不同流场和其他物理因素作用下、在不同作用时间内受到的破坏和改变,这些研究为探索该病症形成和发展的病理机制提供了丰富的理论依据,但是现有技术缺少基于冠状动脉狭窄的血液流动分析的血流动力学因素分析,模拟误差大。
现有技术CN202010033909公开一种基于零维血流动力学模型的FFR快速计算系统模型;所述系统包括:构建个性化零维血流动力学模型,狭窄阻力理论计算模型,确定狭窄阻力,FFR模型。静息状态下,通过冠脉后负荷和包括心脏在内的其他模型参数个性化零维模型。充血状态下,将阻力模型输出的狭窄阻力作为零维模型的输入,影响零维模型冠脉各分支的流量分配,再将零维模型输出的狭窄分支流量作为阻力模型的输入,重新计算狭窄阻力。两个模型如此反复迭代直到流量与阻力匹配,最终确定狭窄阻力,得到冠脉各分支压力,由充血状态下狭窄远端与主动脉根部平均压力的比值计算FFR。该系统可以快速准确地计算FFR,但是冠脉建模过程是根据CTA图像中每一支冠脉血管长度、血管直径d、横截面积A等冠脉的个性化解剖参数来建立物理血管模型,而非根据患者自身影像数据建立个体化三维模型,因此对于后续血流动力学参数的模拟效度会产生影响;模型血流运行采用电路元件(电容等)建立电学参数与血流动力学参数的简化的等效关系,与人体实际的临床数据拟合还存在一定差距,因此其结果降低了临床应用的参考价值;脉压计算采用的是人体正常冠脉流量分配作为计算条件,但FFR是对冠脉狭窄病变的评估,因此条件设置存在缺陷;现有技术CN202010033909虽然实现无创计算FFR,所得血流动力学结果仅限于血流压力计算,没有涉及其他重要血流动力学参数。
发明内容
本发明的目的是针对现有技术问题,本申请基于FFR拟合的冠脉狭窄血流动力学模拟方法及系统,采用不同狭窄病变程度的冠状动脉狭窄患者的影像数据建立个体化冠脉模型;融合临床FFR实测数据,对冠脉模型进行仿真计算,量化血流动力学实时参数分布;对血管实际所处的流体力学环境进行更全面的多维度血流动力学分析,并对血流变异性提供实时的定量描述,分析血管的内皮细胞在不同的流体力学环境中得到的适应性反应,为冠状动脉狭窄提供参考数据。
本申请技术方案如下:
一种基于FFR拟合的冠脉狭窄血流动力学模拟方法,具体包括以下步骤:
S1,建立冠脉狭窄个体化三维模型,基于三角网格顶点特征分解法对冠脉狭窄个体化三维模型进行平滑优化处理,采用曲面离散映射法对血管模型进行网格化处理,获得网格化平滑冠脉狭窄三维模型;
S2,获得临床FFR实测数据,对FFR实测数据进行拟合,加载到网格化平滑冠脉狭窄三维模型的入口边界,量化心动周期内的冠脉狭窄远端压力和主动脉压力输入;
S3,添加时间轴进行血流模拟计算,量化计算血管内任意点及心动周期内任意时刻的血流动力学参数。
步骤S1具体包括以下步骤:
步骤S101,获取影像数据
根据心电图(ECG)获取有效相对静止的心脏时,触发ECG前瞻性门控,进行双源冠状动脉计算机断层扫描血管造影(CCTA),扫描单个心动周期,基于快照冻结获得冠状动脉狭窄患者胸部增强CT医学DICOM格式图像;
进行双源冠状动脉计算机断层扫描血管造影(CCTA),扫描的感兴趣区域(ROI),冠状动脉CTA检查由于需要扫描不停运动的心脏,因此选择较高时间分辨力来“冻结”运动的心脏和冠状动脉,因为心脏运动是根据心电图(ECG)进行的有节律的重复运动,所以根据心电图(ECG)获取有效相对静止的心脏时进行扫描;应用自动触发ECG前瞻性门控扫描,扫描单个心动周期,同时,基于快照冻结获得胸部增强CT医学DICOM格式图像,快照冻结减少心脏运动伪影,提高心脏周期冠状动脉重建图像质量。
步骤S102,建立冠脉狭窄个体化三维模型;
根据人体图像特征,将冠状动脉狭窄患者胸部增强CT医学DICOM格式图像剔除肌肉、软组织、空气和钙化成分,保留冠脉血管;根据冠状动脉相连的主动脉确定冠状动脉位置;使用多阈值适应性算法分割出目标血管,利用图像灰度特征分别计算出主动脉、左右冠脉所属的初始HU值区间,通过多次筛选目标血管部位HU值区间,将CT医学DICOM格式图像生成若干个独立阈值区间,得到分割掩膜像(分割Masks),基于分割掩膜像建立图像三视图,图像三视图包括横断位图像、矢状位图像和冠状位图像,在横断位图像上勾勒冠状动脉边缘,结合冠状位图像和矢状位图像,完成冠状动脉的提取,生成冠脉狭窄个体化三维模型;
步骤S103,对冠脉狭窄个体化三维模型进行平滑处理,具体包括以下步骤:
将冠脉狭窄个体化三维模型表面划分为三角网格,三角网格用于模拟复杂物体的表面,三角网格包括顶点、边和面三要素;
将构建的冠脉狭窄个体化三维模型中三角网格的某一个顶点P为中心,取所有与其相邻的顶点P1…Pn-1和边组成的一阶邻域结构,n表示三角网格的顶点个数,将每个顶点都移动到相邻顶点的平均位置,在冠脉狭窄个体化三维模型的平滑处理中设置法线平滑和顶点拟合的迭代次数和平滑因子,获取平滑后的冠脉狭窄个体化三维模型,定义获取平滑后的冠脉狭窄个体化三维模型为平滑冠脉狭窄三维模型;
步骤S104,对平滑冠脉狭窄三维模型进行曲面网格划分:
对平滑冠脉狭窄三维模型的冠状动脉连续体进行离散处理,利用几何单元逼近,基于变形协调条件完成有限元网格划分,输出曲面网格划分后的平滑冠脉狭窄三维模型。
有限元网格划分采用曲面离散映射法,曲面离散映射法划分网格的过程具体包括以下步骤:根据平滑冠脉狭窄三维模型的边界函数,采用映射矩阵,将边界映射到二维参数空间,在二维参数空间中对所选定的冠状动脉分支区域进行网格划分,将二维参数反向映射到物理空间形成曲面网格。
步骤S2具体包括以下步骤:
步骤S201,获取临床患者FFR实测数据,FFR实测数据包括冠脉狭窄远端压力和主动脉压力的实测值,Pc-d表示第d个时刻的冠脉狭窄远端压力,Pa-d表示第d个时刻的主动脉压;
步骤S202,将实测冠脉狭窄远端压力和主动脉压进行曲线拟合,获得一个心动周期内不同入口边界条件下的流量曲线;
步骤S203,设置血液密度,根据冠脉狭窄程度设定各个患者的血液粘稠度;将FFR拟合后数据作为冠脉周期脉动血流压力载入血管入口边界,出口设置为各患者FFR中拟合后的冠脉狭窄远端压力的周期平均值,实现血流心动周期性实际压差;
实现血流心动周期性实际压差,最终能量收敛,满足质量与能量守恒。
步骤S202具体包括以下步骤:
S202-1,对冠脉狭窄远端压力进行拟合,获得拟合后的冠脉狭窄远端压力;
基于预测的FFR线性模型函数对冠脉狭窄远端压力进行拟合;以采样的血流时刻为横坐标,横坐标为d,d=1,2…D,d表示第d个采样时刻,D表示采样数量,以实测的冠脉狭窄远端压力为纵坐标,d为常数且d∈[1,D];拟合输出一个心动周期内冠脉狭窄远端压力曲线,动脉血流冠脉狭窄远端压力的预测值为:
f(pc_(d+1))=ω1Pc-12Pc-2+…ωdPc-d+b (1)
ωd表示基于FFR线性模型函数对冠脉狭窄远端压力进行拟合时多项式拟合的第d个系数,ωd属于常数;
式(1)为冠脉狭窄远端压力的FFR线性模型函数,f(pc_(d+1))表示第d+1个采样时刻冠脉狭窄远端压力的预测值;
冠脉狭窄远端压力的FFR线性模型函数输出预测的冠脉狭窄远端压力,预测的冠脉狭窄远端压力和冠脉狭窄远端压力的实测值构成损失函数,通过最小化损失函数求解ω和b并评估函数;Pc-d为一个心动周期内第d个采样时刻冠脉狭窄远端压力的实测值;
f(pc_(d+1))表示第d+1个采样时刻冠脉狭窄远端压力的预测值;
f(pc_(d+1))=ωTXPc+b (2)
式(2)为冠脉狭窄远端压力的FFR线性模型函数的向量表示,其中ω为向量系数;ω为ω1、ω2…ωd的向量表示,ωT表示ω的转置;
XPc是FFR线性模型函数的中Pc-1、Pc-2…Pc-d以行向量形式表示;
令xi=f(pc_d),yi=Pc-d,i=1,2…d…D;
f(pc_d)表示第d个采样时刻冠脉狭窄远端压力的预测值,Pc-d表示d时刻冠脉狭窄远端压力的实测值;式(3)为冠脉狭窄远端压力的损失函数,用于衡量预测值xi与实测值yi的差值,其中b*为b的偏置,ω*为ω向量的参数矩阵;
式(4)对两个参数ω和b求偏导,令式(4)等于0得到ω和b的最优解;b为函数FFR线性模型函数的参数,b*为参数b的偏置项;
i表示采样序号,i∈[1,D],D为采样数量;表示x1、x2…xD的均值;
式(5)和(6)分别为ω和b的最优解;
计算出ω和b的最优解后,代入公式(1),计算输出拟合后的冠脉狭窄远端压力;
S202-2,对主动脉压进行拟合,获得拟合后的主动脉压;
利用预测的FFR线性模型函数对测量好的主动脉压进行拟合;以采样的血流时刻为横坐标,横坐标为d,d=1,2…D,d表示第d个采样时刻,D表示采样数量,以实测的主动脉压为纵坐标,d为常数且d∈[1,D];拟合输出的是一个心动周期内主动脉压曲线,输出动脉血流实时主动脉压为:
f(Pa_(d+1))=ωa1Pa-1a2Pa-2+…ωadPa-d+ba (7)
式(7)中,ωad表示基于FFR线性模型函数对主动脉压进行拟合时多项式拟合的第d个系数,ωad属于常数,ba表示常数项;
式(7)为主动脉压的FFR线性模型函数,f(Pa_(d+1))表示第d+1个采样时刻主动脉压的预测值;
主动脉压的FFR线性模型函数输出预测的主动脉压,预测的主动脉压和主动脉压的实测值构成损失函数,通过最小化损失函数求解ωa和ba并评估函数;Pa-d为一个心动周期内第d个采样时刻主动脉压的实测值。
d表示采样时刻;d=1,2…D;
f(Pa_(d+1))=ωa TXPa+ba (8);
式(8)为主动脉压的FFR线性模型函数的向量表示,其中ωa为向量;ωa为ωa1、ωa2…ωad的向量表示,ωa T表示ωa的转置;
XPa是预测线性模型函数的中Pa-1、Pa-2…Pa-d以行向量形式表示;
令xa_i=f(Pa_d),ya_i=Pa-d,i=1,2…d…D;
f(Pa_d)表示第d个采样时刻主动脉压的预测值;Pa-d表示d时刻主动脉压的实测值;式(9)为损失函数,用于衡量预测值xi与实测值yi的差值,其中ba *为ba的偏置,ωa *为ωa向量的参数矩阵;
式(10)对两个参数ωa和ba求偏导,令式(10)等于0得到ωa和ba的最优解;ba为函数f(Pa_d)的参数,b* a为参数ba的偏置项;
i表示采样序号,i∈[1,D],D为采样数量;表示xa_1、xa_2…xa_D的均值;
式(11)、(12)为ωa和ba的最优解;计算出ωa和ba的最优解后,代入公式(7),计算输出拟合后的主动脉压;
步骤S3计算具体包括以下步骤:
步骤S301,在所述血流心动周期性实际压差条件下,通过单向流体渗流Stokes方程仿真冠脉血流流动,仿真计算结果包括血流动力学参量心动周期内全过程演变;血流动力学参量包括瞬时流量、血管层流断面平均流速、湍流强度、血管壁面剪切应力;
Stokes方程的作用是是确保血流动力学计算中遵循能量守恒,动量守恒,质量守恒。在恒定压力条件下,单向流体渗流Stokes方程为:
stokes方程是用于描述流体运动的方程,看作是流体运动的牛顿第二定律。该步骤中通过stokes方程求得u,根据u即得出仿真计算所需的血流流速以及时间。
其中P指的是血管内血流总体压力;μ为血流粘度系数;u为血流速度向量;F为体积力向量;ρ为固定常量血流密度,T为绝对温度293.15K;I表示单位矩阵、是梯度算子(在空间各方向上的全微分);/>是对血流速度求梯度,是指向速度标量场增长最快的方向;指/>即梯度算子点乘uu函数;/>表示各方向上的矢量和;
步骤S302,血流动力学参数包括瞬时流量、血管层流断面平均流速、湍流强度和血管壁面剪切应力;
在血管内取一个半径为r、长度为L的圆柱形流体段。
1)瞬时流量Q:
K为稠度系数,L为血液流过的流体段长度,n为流变指数,R为血管内径;
Δp=2Lτ/r,其中τ为血管壁面剪切应力,r为流体段半径,Δp为血流流过流体段压差;
τ由下方公式(18)得出;
2)血管层流断面平均流速V:
其中V为血管层流断面平均流速,μm为收缩期峰值时刻速度。
3)湍流强度B:
B=0.16*Re-1/8 (17)
B为湍流强度;Re为雷诺数;
4)WSS(血管壁面剪切应力)τ:
其中τ表示WSS(血管壁面剪切应力),R为血管内径,γ代表血管横截面上任意点到圆心的距离,当γ=R时计算WSS。
实时追踪冠脉任意点的血流参数变化,量化不同狭窄几何构型耦合作用下冠脉血流动力学动态特性演变规律。
一种基于FFR拟合的冠脉狭窄血流动力学模拟系统,具体包括模型建立单元、拟合单元和动力学参数计算单元;
模型建立单元建立冠脉狭窄个体化三维模型,基于三角网格顶点特征分解法对冠脉狭窄个体化三维模型进行平滑优化处理,采用曲面离散映射法对血管模型进行网格化处理,获得网格化平滑冠脉狭窄三维模型;
拟合单元融合临床FFR实测数据,对FFR实测数据进行拟合,加载到网格化平滑冠脉狭窄三维模型的入口边界,量化心动周期内的时序血流压力输入;时序血流压力包括冠脉狭窄远端压力和主动脉压力;
动力学参数计算单元添加时间轴进行血流模拟计算,量化计算血管内任意点及心动周期内任意时刻的血流动力学参数。
模型建立单元工作过程包括以下步骤:
步骤S101,获取影像数据
根据心电图(ECG)获取有效相对静止的心脏时,触发ECG前瞻性门控,进行双源冠状动脉计算机断层扫描血管造影(CCTA),扫描单个心动周期,基于快照冻结获得冠状动脉狭窄患者胸部增强CT医学DICOM格式图像;
进行双源冠状动脉计算机断层扫描血管造影(CCTA),扫描的感兴趣区域(ROI)为从气管隆突到心底,包括整个心脏,冠状动脉CTA检查由于需要扫描不停运动的心脏,因为心脏运动是根据心电图(ECG)进行的有节律的重复运动,所以根据心电图(ECG)获取有效相对静止的心脏时进行扫描;应用自动触发ECG前瞻性门控扫描,扫描单个心动周期,同时,基于快照冻结获得胸部增强CT医学DICOM格式图像,快照冻结减少心脏运动伪影,提高心脏周期冠状动脉重建图像质量。
步骤S102,建立冠脉狭窄个体化三维模型;
根据人体图像特征,将冠状动脉狭窄患者胸部增强CT医学DICOM格式图像剔除肌肉、软组织、空气、钙化成分,保留冠脉血管;根据冠状动脉相连的主动脉确定冠状动脉位置;使用多阈值适应性算法分割出目标血管,利用图像灰度特征分别计算出主动脉、左右冠脉所属的初始HU值区间范围,通过多次筛选目标血管部位HU值区间,将CT医学DICOM格式图像生成若干个独立阈值区间,得到分割掩膜像(分割Masks),基于分割掩膜像建立图像三视图,图像三视图包括横断位图像、矢状位图像和冠状位图像,在横断位图像上勾勒冠状动脉边缘,结合冠状位图像和矢状位图像,完成冠状动脉的提取,生成冠脉狭窄个体化三维模型;
步骤S103,对冠脉狭窄个体化三维模型进行平滑处理,具体包括以下步骤:
将冠脉狭窄个体化三维模型表面划分为三角网格,三角网格用于模拟复杂物体的表面,三角网格包括顶点、边和面三要素;
将构建的冠脉狭窄个体化三维模型中三角网格的某一个顶点P为中心,取所有与其相邻的顶点P1…Pn-1和边组成的一阶邻域结构,n表示三角网格的顶点个数,将每个顶点都移动到相邻顶点的平均位置,在冠脉狭窄个体化三维模型的平滑处理中设置法线平滑和顶点拟合的迭代次数和平滑因子,获取平滑后的冠脉狭窄个体化三维模型,定义获取平滑后的冠脉狭窄个体化三维模型为平滑冠脉狭窄三维模型;
步骤S104,对平滑冠脉狭窄三维模型进行曲面网格划分:
对平滑冠脉狭窄三维模型的冠状动脉连续体进行离散处理,利用几何单元逼近,基于变形协调条件完成有限元网格划分,输出曲面网格划分后的平滑冠脉狭窄三维模型;
有限元网格划分采用曲面离散映射法,曲面离散映射法划分网格的过程具体包括以下步骤:根据平滑冠脉狭窄三维模型的边界函数,采用映射矩阵,将边界映射到二维参数空间,在二维参数空间中对所选定的冠状动脉分支区域进行网格划分,将二维参数反向映射到物理空间形成曲面网格。
拟合单元工作过程具体包括以下步骤:
步骤S201,获取临床患者FFR实测数据,FFR实测数据包括冠脉狭窄远端压力和主动脉压力,记为Pc-d和Pa-d,Pc-d表示第d个时刻的冠脉狭窄远端压力,Pa-d表示第d个时刻的主动脉压;;
步骤S202,采用将实测冠脉狭窄远端压力和主动脉压进行曲线拟合,获得一个心动周期内不同入口边界条件下的流量曲线;
步骤S203,设置血液密度,根据冠脉狭窄程度设定各个患者的血液粘稠度;将患者FFR拟合后数据作为冠脉周期脉动血流压力载入血管入口边界,出口设置为各患者FFR中拟合后的冠脉狭窄远端压力和主动脉压的周期平均值,实现血流心动周期性实际压差;
实现血流心动周期性实际压差,最终能量收敛,满足质量与能量守恒。
FFR数据进行曲线拟合具体包括以下步骤:
S202-1,对冠脉狭窄远端压力进行拟合,获得拟合后的冠脉狭窄远端压力;
利用预测的FFR线性模型函数对测量好的冠脉狭窄远端压力Pc进行拟合;以采样的血流时刻为横坐标,横坐标为d,d=1,2…D,d表示第d个采样时刻,D表示采样数量,以实测的冠脉狭窄远端压力为纵坐标,d为常数且d∈[1,D];拟合输出一个心动周期内冠脉狭窄远端压力曲线,输出动脉血流冠脉狭窄远端压力的预测值为:
f(pc_(d+1))=ω1Pc-12Pc-2+…ωdPc-d+b (1)
ωd表示基于FFR线性模型函数对冠脉狭窄远端压力进行拟合时多项式拟合的第d个系数,ωd属于常数;
式(1)为冠脉狭窄远端压力的FFR线性模型函数,f(pc_(d+1))表示第d+1个采样时刻冠脉狭窄远端压力的预测值;
冠脉狭窄远端压力的FFR线性模型函数输出预测的冠脉狭窄远端压力,预测的冠脉狭窄远端压力和冠脉狭窄远端压力的实测值构成损失函数,通过最小化损失函数求解ω和b并评估函数;Pc-d为一个心动周期内第d个采样时刻冠脉狭窄远端压力的实测值;
f(pc_(d+1))第d+1个采样时刻冠脉狭窄远端压力的预测值;
f(pc_(d+1))=ωTXPc+b (2)
式(2)为冠脉狭窄远端压力的FFR线性模型函数的向量表示,其中ω为向量系数;ω为ω1、ω2…ωd的向量表示,ωT表示ω的转置;
XPc是FFR线性模型函数的中Pc-1、Pc-2…Pc-d以行向量形式表示;
令xi=f(pc_d),yi=Pc-d,i=1,2…d…D;
f(pc_d)表示第d个采样时刻冠脉狭窄远端压力的预测值,Pc-d表示d时刻冠脉狭窄远端压力的实测值;式(3)为冠脉狭窄远端压力的损失函数,用于衡量预测值f(xi)与实测值yi的差值,其中b*为b的偏置,ω*为ω向量的参数矩阵;
式(4)对两个参数ω和b求偏导,令式(4)等于0得到ω和b的最优解;b为函数FFR线性模型函数f(Pc-d)的参数,b*为参数b的偏置项;
i表示采样序号,i∈[1,D],D为采样数量;表示x1、x2…xD的均值;
式(5)和(6)为ω和b的最优解;
计算出ω和b的最优解后,代入公式(1),计算输出拟合后的冠脉狭窄远端压力;
S202-2,对主动脉压进行拟合,获得拟合后的主动脉压;
利用预测的FFR线性模型函数对测量好的主动脉压进行拟合;以采样的血流时刻为横坐标,横坐标为d,d=1,2…D,d表示第d个采样时刻,D表示采样数量,以实测的主动脉压为纵坐标,d为常数且d∈[1,D];拟合输出的是一个心动周期内主动脉压曲线,输出动脉血流实时主动脉压为:
f(Pa_(d+1))=ωa1Pa-1a2Pa-2+…ωadPa-d+ba (7)
式(7)中,ωad表示基于FFR线性模型函数对主动脉压进行拟合时多项式拟合的第d个系数,ωad属于常数,ba表示常数项;
式(7)为主动脉压的FFR线性模型函数,f(Pa_(d+1))表示第d个采样时刻主动脉压的预测值;
主动脉压的FFR线性模型函数输出预测的主动脉压,预测的主动脉压和主动脉压的实测值构成损失函数,通过最小化损失函数求解ωa和ba并评估函数;Pa-d为一个心动周期内第d个采样时刻主动脉压的实测值。
d表示采样时刻;d=1,2…D;
f(Pa_(d+1))=ωa TXPa+ba (8)
式(8)为主动脉压的FFR线性模型函数的向量表示,其中ωa为向量;ωa为ωa1、ωa2…ωad的向量表示,ωa T表示ωa的转置;
XPa是预测线性模型函数的未知数Pa-1、Pa-2…Pa-d以行向量形式表示;
令xa_i=f(Pa_d),ya_i=pa-d,i=1,2…d…D;
f(Pa_d)表示第d个采样时刻主动脉压的预测值;Pa-d表示d时刻主动脉压的实测值;式(9)为损失函数,用于衡量预测值xa_i与实测值yi的差值,其中ba *为ba的偏置,ωa *为ωa向量的参数矩阵;
式(10)对两个参数ωa和ba求偏导,令式(10)等于0得到ωa和ba的最优解;ba为函数f(Pa_d)的参数,b* a为参数ba的偏置项;
i表示采样序号,i∈[1,D],D为采样数量;
式(11)、(12)为ωa和ba的最优解;计算出ωa和ba的最优解后,代入公式(7),计算输出拟合后的主动脉压;
血流动力学参数计算单元工作过程具体包括以下步骤:
步骤S301,在所述血流心动周期性实际压差条件下,通过单向流体渗流Stokes方程仿真冠脉血流流动,仿真计算结果包括血流动力学参量心动周期内全过程演变;血流动力学参量包括瞬时流量、血管层流断面平均流速、湍流强度、血管壁面剪切应力;
Stokes方程的作用是是确保血流动力学计算中遵循能量守恒,动量守恒,质量守恒。在恒定压力条件下,单向流体渗流Stokes方程为:
stokes方程是用于描述流体运动的方程,看作是流体运动的牛顿第二定律。该步骤中通过stokes方程求得u,根据u即得出仿真计算所需的血流流速以及时间。
其中P指的是血管内血流总体压力;μ为血流粘度系数;u为血流速度向量;F为体积力向量;ρ为固定常量血流密度,公式(13)中T为绝对温度293.15K,I表示单位矩阵、是梯度算子(在空间各方向上的全微分);/>是对血流速度求梯度,是指向速度标量场增长最快的方向;/>指/>表示梯度算子点乘uu函数;/>表示各方向上的矢量和;
步骤S302,血流动力学参数包括瞬时流量、血管层流断面平均流速、湍流强度、血管壁面剪切应力;
在血管内取一个半径为r、长度为L的圆柱形流体段。
1)瞬时流量Q:
K为稠度系数、L为血液流过的流体段、n为流变指数、R为血管内径;
Δp=2Lτ/r,其中τ为血管壁面剪切应力,r为流体段半径,Δp为血流流过流体段压差;
τ可由下方公式(18)得出:
得出血管层流断面平均流速V:
其中V为血管层流断面平均流速,μm为收缩期峰值时刻速度。
2)湍流强度B:
B=0.16*Re-1/8 (17)
B湍流强度;Re为雷诺数;
3)WSS(血管壁面剪切应力)τ:
其中τ表示WSS(血管壁面剪切应力),R为血管内径,γ代表血管横截面上任意点到圆心的距离,当γ=R时计算WSS。
实时追踪冠脉任意点的血流参数变化,量化不同狭窄几何构型耦合作用下冠脉血流动力学动态特性演变规律。
相对于现有技术,本发明具有如下有益效果:
本申请公开了一种基于冠脉狭窄时间分辨率上的血流动力学数值模拟方法,实现了对于不同狭窄病变患者冠脉血管实现个体化精准建模;采用映射法形成冠脉模型的曲面网格划分;利用拟合患者实测FFR作为边界条件实现时间分辨率上的冠脉血流动力学数值模拟。
本申请对血管进行动态模拟,在一个心动周期内患者临床实测FFR数据的拟合曲线作为fluent软件中模型的出入口条件,对冠脉进行动态模拟;这样设置边界,使冠脉内血流的流态更接近于实际情况。
本申请在冠状动脉建模时,采用了分割算法,分割出的冠状动脉比由普遍采用的区域生长算法重建的模型更符合人体血管的特征。
本申请获得了冠脉出现狭窄时间分辨率上的血流动力学分布,得到了量化的,能够更直接的描述血流的流变模式。
本申请应用计算流体力学对冠状动脉患者进行数值模拟分析,以此获得冠状动脉的相关血流动力学参数,从而以无创的方式实现对冠脉狭窄血流动力学时间分辨率模拟。
附图说明
通过阅读参照以下附图所作的对非限制性实施例所作的详细描述,本申请的其它特征、目的和优点将会变得更明显:
附图1一种基于FFR拟合的冠脉狭窄血流动力学模拟方法流程图;
附图2本实施例患者FFR数据拟合结果---重度冠脉狭窄患者FFR拟合后数据曲线;
附图3本实施例患者FFR数据拟合结果---轻度冠脉狭窄患者FFR拟合后数据曲线。
附图4本实施例圆柱形流体段示意图。
具体实施方式
下面将结合本发明中的附图,对本发明的技术方案进行清楚、完整地描述,显然,所描述的实施例仅仅是本发明一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动条件下所获得的所有其它实施例,都属于本发明保护的范围。
如图1所示,一种基于FFR拟合的冠脉狭窄血流动力学模拟方法,具体包括以下步骤:
S1,建立冠脉狭窄个体化三维模型,基于三角网格顶点特征分解法对冠脉狭窄个体化三维模型进行平滑优化处理,采用曲面离散映射法对血管模型进行网格化处理,获得网格化平滑冠脉狭窄三维模型;
S2,融合临床FFR实测数据,对FFR实测数据进行拟合,加载到网格化平滑冠脉狭窄三维模型的入口边界,量化心动周期内的时序血流压力输入;时序血流压力包括冠脉狭窄远端压力和主动脉压力;
S3,添加时间轴进行血流模拟计算,量化计算血管内任意点及心动周期内任意时刻的血流动力学参数。
步骤S1具体包括以下步骤:
步骤S101,获取影像数据
根据心电图(ECG)获取有效相对静止的心脏时,触发ECG前瞻性门控,进行双源冠状动脉计算机断层扫描血管造影(CCTA),扫描单个心动周期,基于快照冻结获得冠状动脉狭窄患者胸部增强CT医学DICOM格式图像;
进行双源冠状动脉计算机断层扫描血管造影(CCTA),扫描的感兴趣区域(ROI)为从气管隆突到心底,包括整个心脏,冠状动脉CTA检查由于需要扫描不停运动的心脏,因此选择较高时间分辨力来“冻结”运动的心脏和冠状动脉,因为心脏运动是根据心电图(ECG)进行的有节律的重复运动,所以根据心电图(ECG)获取有效相对静止的心脏时进行扫描;应用自动触发ECG前瞻性门控扫描,扫描单个心动周期,同时,基于快照冻结获得胸部增强CT医学DICOM格式图像,快照冻结减少心脏运动伪影,提高心脏周期冠状动脉重建图像质量。
步骤S102,建立冠脉狭窄个体化三维模型;
将冠状动脉狭窄患者胸部增强CT医学DICOM格式图像导入Mimics软件(MIMICS是Materialise's interactive medical image control system),是Materialise公司发明的一种医学图像处理建模软件)中.根据人体图像特征(人体组织密度、图像HU值(HU:Hounsfiled Unit CT值)、各组织不同阈值),剔除肌肉、软组织、空气、钙化成分,保留冠脉血管;根据冠状动脉相连的主动脉确定冠状动脉位置;使用多阈值适应性算法分割出目标血管,利用图像灰度特征分别计算出主动脉、左右冠脉所属的初始HU值区间范围[A,B],通过多次筛选目标血管部位HU值区间,将CT医学DICOM格式图像生成若干个独立阈值区间,得到分割掩膜像(分割Masks),基于分割掩膜像建立图像三视图,图像三视图包括横断位图像、矢状位图像和冠状位图像,在横断位图像上勾勒冠状动脉边缘,结合冠状位图像和矢状位图像,完成冠状动脉的提取,生成冠脉狭窄个体化三维模型;
步骤S103,对冠脉狭窄个体化三维模型进行平滑处理,具体包括以下步骤:
将冠脉狭窄个体化三维模型表面划分为三角网格,三角网格用于模拟复杂物体的表面,三角网格包括顶点、边和面三要素;
将构建的冠脉狭窄个体化三维模型中三角网格的某一个顶点P为中心,取所有与其相邻的顶点P1…Pn-1和边组成的一阶邻域结构,n表示三角网格的顶点个数,将每个顶点都移动到相邻顶点的平均位置,在冠脉狭窄个体化三维模型的平滑处理中设置法线平滑和顶点拟合的迭代次数为4,平滑因子为0.4,获取平滑后的冠脉狭窄个体化三维模型,定义获取平滑后的冠脉狭窄个体化三维模型为平滑冠脉狭窄三维模型,平滑冠脉狭窄三维模型是光滑后的三维模型,应用于Fluent进行流体力学分析。
步骤S104,对平滑冠脉狭窄三维模型进行曲面网格划分:
对平滑冠脉狭窄三维模型的冠状动脉连续体进行离散处理,利用几何单元逼近,基于变形协调条件完成有限元网格划分,输出曲面网格划分后的平滑冠脉狭窄三维模型。
将平滑冠脉狭窄三维模型导入网格划分模块(ANSYS ICEM软件),选定所导入的平滑冠脉狭窄三维模型的冠状动脉分支区域,确定物理场,对所导入的平滑冠脉狭窄三维模型的冠状动脉连续体(“冠状动脉连续体”:连续体概念指的是描述具有力学性质的物体,本实施例指具有弹性的冠状动脉血管模型;)进行离散处理,利用几何单元逼近,基于变形协调条件完成有限元网格划分,输出曲面网格划分后的平滑冠脉狭窄三维模型。
有限元网格划分采用曲面离散映射法,曲面离散映射法划分网格的过程具体包括以下步骤:根据平滑冠脉狭窄三维模型的边界函数,采用映射矩阵,将边界映射到二维参数空间,在二维参数空间中对所选定的冠状动脉分支区域进行网格划分,将二维参数反向映射到物理空间形成曲面网格。
步骤S2具体包括以下步骤:
步骤S201,获取临床患者FFR实测数据,FFR实测数据包括冠脉狭窄远端压力和主动脉压力的实测值:将压力导丝传感器放置在冠脉狭窄处远端5cm处,冠脉注射200ug硝酸甘油,等血压恢复后,静脉泵入ATP/腺苷,移动指针向右到最大充血状态的波形平台期同时开始记录,分别采样点测量的冠脉狭窄远端压力和主动脉压,实测FFR数据包括不同标记时刻的冠脉狭窄远端压力和主动脉压的实测值,记为Pc-d和Pa-d,Pc-d表示第d个时刻的冠脉狭窄远端压力,Pa-d表示第d个时刻的主动脉压;在实验中不同时刻,0.0秒、0.2秒、0.4秒、0.6秒、0.8秒、1.0秒分别测量冠脉狭窄远端压力和主动脉压。
步骤S202,将实测冠脉狭窄远端压力和主动脉压进行曲线拟合,获得一个心动周期内不同入口边界条件下的流量曲线;
步骤S203,设置血液密度,将血液密度设置为1060kg/m3,根据冠脉狭窄程度设定各个患者的血液粘稠度(男性4.57±0.73mPa·s,女性3.90±0.51mPa·s);将患者FFR拟合后数据作为冠脉周期脉动血流压力载入血管入口边界,出口设置为各患者FFR中拟合后的冠脉狭窄远端压力的周期平均值,实现血流心动周期性实际压差;
实现血流心动周期性实际压差,最终能量收敛,满足质量与能量守恒。
步骤S2,获取临床FFR实测数据,并对FFR实测数据进行算法拟合,编程加载到曲面网格划分后的平滑冠脉狭窄三维模型的入口边界(入口边界包括流动变量和热变量在边界处的值),量化心动周期内的时序血流压力输入;时序血流压力包括冠脉狭窄远端压力和主动脉压力,冠脉狭窄远端压力和主动脉压力均属于动脉血管不同部位的血流压力,FFR是分别拟合冠脉狭窄远端压和主动脉压,分别输出的是拟合后冠脉狭窄远端压和主动脉压;步骤S202具体包括以下步骤:
S202-1,对冠脉狭窄远端压力进行拟合,获得拟合后的冠脉狭窄远端压力;
基于预测的FFR线性模型函数对冠脉狭窄远端压力进行拟合;以采样的血流时刻为横坐标,横坐标为d,d=1,2…D,d表示第d个采样时刻,D表示采样数量,以实测的冠脉狭窄远端压力为纵坐标,d为常数且d∈[1,D];拟合输出一个心动周期内冠脉狭窄远端压力曲线,输出动脉血流实时冠脉狭窄远端压力为:
f(Pc_(d+1))=ω1Pc-12Pc-2+…ωdPc-d+b (1)
ωd表示基于FFR线性模型函数对冠脉狭窄远端压力进行拟合时多项式拟合的第d个系数,ωd属于常数;
式(1)为冠脉狭窄远端压力的FFR线性模型函数,f(pc_(d+1))表示第d+1个采样时刻冠脉狭窄远端压力的预测值;
冠脉狭窄远端压力的FFR线性模型函数输出预测的冠脉狭窄远端压力,预测的冠脉狭窄远端压力和冠脉狭窄远端压力的实测值构成损失函数,通过最小化损失函数求解ω和b并评估函数;Pc-d为一个心动周期内第d个采样时刻冠脉狭窄远端压力的实测值;
f(pc_(d+1))表示第d+1个采样时刻冠脉狭窄远端压力的预测值;
f(pc_(d+1))=ωTXPc+b (2)
式(2)为冠脉狭窄远端压力的FFR线性模型函数的向量表示,其中ω为向量系数;ω为ω1、ω2…ωd的向量表示,ωT表示ω的转置;
XPc是FFR线性模型函数(预测线性模型函数)的中Pc-1、Pc-2…Pc-d以行向量形式表示;
令xi=f(pc_d),yi=Pc-d,i=1,2…d…D;
f(pc_d)表示第d个采样时刻冠脉狭窄远端压力的预测值,Pc-d表示d时刻冠脉狭窄远端压力的实测值;式(3)为冠脉狭窄远端压力的损失函数,用于衡量预测值f(xi)与实测值yi的差值,其中b*为b的偏置,ω*为ω向量的参数矩阵;
式(4)对两个参数ω和b求偏导,令式(4)等于0得到ω和b的最优解;b为函数FFR线性模型函数的参数,b*为参数b的偏置项;
i表示采样序号,i∈[1,D],D为采样数量;表示x1、x2…xD的均值;
式(5)和(6)为ω和b的最优解;
计算出ω和b的最优解后,代入公式(1),计算输出拟合后的冠脉狭窄远端压力;
S202-2,对主动脉压进行拟合,获得拟合后的主动脉压;
利用预测的FFR线性模型函数对测量好的主动脉压进行拟合;以采样的血流时刻为横坐标,横坐标为d,d=1,2…D,d表示第d个采样时刻,D表示采样数量,以实测的主动脉压为纵坐标,d为常数且d∈[1,D];拟合输出的是一个心动周期内主动脉压曲线,输出动脉血流实时主动脉压为:
f(Pa_(d+1))=ωa1Pa-1a2Pa-2+…ωadPa-d+ba (7)
式(7)中,ωad表示基于FFR线性模型函数对主动脉压进行拟合时多项式拟合的第d个系数,ωad属于常数,ba表示常数项;
式(7)为主动脉压的FFR线性模型函数,f(Pa_(d+1))表示第d+1个采样时刻主动脉压的预测值;
主动脉压的FFR线性模型函数输出预测的主动脉压,预测的主动脉压和主动脉压的实测值构成损失函数,通过最小化损失函数求解ωa和ba并评估函数;Pa-d为一个心动周期内第d个采样时刻(第0.00秒、0.02秒、0.04秒、0.06秒……)主动脉压的实测值。
d表示采样时刻;d=1,2…D;
f(Pa_(d+1))=ωa TXPa+ba (8);
式(8)为主动脉压的FFR线性模型函数的向量表示,其中ωa为向量;ωa为ωa1、ωa2…ωad的向量表示,ωa T表示ωa的转置;
XPa是预测线性模型函数的未知数Pa-1、Pa-2…Pa-d以行向量形式表示;
令xa_i=f(Pa_d),ya_i=Pa-d,i=1,2…d…D;
f(Pa_d)表示第d个采样时刻主动脉压的预测值;Pa-d表示d时刻主动脉压的实测值;式(9)为损失函数,用于衡量预测值xa_i与实测值ya_i的差值,其中ba *为ba的偏置,ωa *为ωa向量的参数矩阵;
式(10)对两个参数ωa和ba求偏导,令式(10)等于0得到ωa和ba的最优解;ba为函数f(Pa_d)的参数,b* a为参数ba的偏置项;
i表示采样序号,i∈[1,D],D为采样数量;表示xa_1、xa_2…xa_D的均值;
式(11)、(12)为ωa和ba的最优解;计算出ωa和ba的最优解后,代入公式(7),计算输出拟合后的主动脉压;
添加“Timeliner”时间轴将各患者FFR拟合后数据作为边界条件载入求解函数进行3D冠脉模型进行血流模拟计算;实现时间分辨率维度的冠脉血流演变全过程,量化计算血管内任意点及心动周期内任意时刻的血流动力学参数;步骤S3计算具体包括以下步骤:
S301,在所述血流心动周期性实际压差条件下,通过单向流体渗流Stokes方程仿真冠脉血流流动,仿真计算结果包括血流动力学参量心动周期内全过程演变;血流动力学参量包括瞬时流量、血管层流断面平均流速、湍流强度、血管壁面剪切应力;
Stokes方程的作用是是确保血流动力学计算中遵循能量守恒,动量守恒,质量守恒。在恒定压力条件下,单向流体渗流Stokes方程为:
stokes方程是用于描述流体运动的方程,看作是流体运动的牛顿第二定律。该步骤中通过stokes方程求得u,根据u即得出仿真计算所需的血流流速以及时间。
其中P指的是血管内血流总体压力;μ为血流粘度系数;u为血流速度向量;F为体积力向量;ρ为固定常量血流密度,T为绝对温度293.15K;I表示单位矩阵、是梯度算子(在空间各方向上的全微分);/>是对血流速度求梯度,是指向速度标量场增长最快的方向;/>指/>即梯度算子点乘uu函数;/>表示各方向上的矢量和;
根据牛顿粘性定律:血流粘度μ=τ·dv/dy,其中τ表示作用在流体平面上的剪应力,dv/dy是剪切应变率,斜率μ为血流粘度系数。
S302,血流动力学参数包括瞬时流量、血管层流断面平均流速、湍流强度和血管壁面剪切应力;
在血管内取一个半径为r(正常冠状动脉r为0.1cm)、长度为L的圆柱形流体段。
1)瞬时流量Q:
K为稠度系数,L为血液流过的流体段长度,n为流变指数,R为血管内径;
Δp=2Lτ/r,其中τ为血管壁面剪切应力,r为流体段半径,Δp为血流流过流体段压差;
τ由下方公式(18)得出;
2)血管层流断面平均流速V:
其中V为血管层流断面平均流速,μm为收缩期峰值时刻速度。
3)湍流强度B:
B=0.16*Re-1/8 (17)
B为湍流强度;Re为雷诺数;
4)WSS(血管壁面剪切应力)τ:
其中τ表示WSS(血管壁面剪切应力),R为血管内径,γ代表血管横截面上任意点到圆心的距离,当γ=R时计算WSS。
实时追踪冠脉任意点的血流参数变化,量化不同狭窄几何构型耦合作用下冠脉血流动力学动态特性演变规律。
一种基于FFR拟合的冠脉狭窄血流动力学模拟系统,具体包括模型建立单元、拟合单元和动力学参数计算单元;
模型建立单元建立冠脉狭窄个体化三维模型,基于三角网格顶点特征分解法对冠脉狭窄个体化三维模型进行平滑优化处理,采用曲面离散映射法对血管模型进行网格化处理,获得网格化平滑冠脉狭窄三维模型;
拟合单元融合临床FFR实测数据,对FFR实测数据进行拟合,加载到网格化平滑冠脉狭窄三维模型的入口边界,量化心动周期内的时序血流压力输入;时序血流压力包括冠脉狭窄远端压力和主动脉压力;
动力学参数计算单元添加时间轴进行血流模拟计算,量化计算血管内任意点及心动周期内任意时刻的血流动力学参数。
模型建立单元工作过程包括以下步骤:
步骤S101,获取影像数据
根据心电图(ECG)获取有效相对静止的心脏时,触发ECG前瞻性门控,进行双源冠状动脉计算机断层扫描血管造影(CCTA),扫描单个心动周期,基于快照冻结获得冠状动脉狭窄患者胸部增强CT医学DICOM格式图像;
进行双源冠状动脉计算机断层扫描血管造影(CCTA),扫描的感兴趣区域(ROI)为从气管隆突到心底,包括整个心脏,冠状动脉CTA检查由于需要扫描不停运动的心脏,因此选择较高时间分辨力来“冻结”运动的心脏和冠状动脉,因为心脏运动是根据心电图(ECG)进行的有节律的重复运动,所以根据心电图(ECG)获取有效相对静止的心脏时进行扫描;应用自动触发ECG前瞻性门控扫描,扫描单个心动周期,同时,基于快照冻结获得胸部增强CT医学DICOM格式图像,快照冻结减少心脏运动伪影,提高心脏周期冠状动脉重建图像质量。
步骤S102,建立冠脉狭窄个体化三维模型;
根据人体图像特征,将冠状动脉狭窄患者胸部增强CT医学DICOM格式图像剔除肌肉、软组织、空气、钙化成分,保留冠脉血管;根据冠状动脉相连的主动脉确定冠状动脉位置;使用多阈值适应性算法分割出目标血管,利用图像灰度特征分别计算出主动脉、左右冠脉所属的初始HU值区间范围[A,B],通过多次筛选目标血管部位HU值区间,将CT医学DICOM格式图像生成若干个独立阈值区间,得到分割掩膜像(分割Masks),基于分割掩膜像建立图像三视图,图像三视图包括横断位图像、矢状位图像和冠状位图像,在横断位图像上勾勒冠状动脉边缘,结合冠状位图像和矢状位图像,完成冠状动脉的提取,生成冠脉狭窄个体化三维模型;
步骤S103,对冠脉狭窄个体化三维模型进行平滑处理,具体包括以下步骤:
将冠脉狭窄个体化三维模型表面划分为三角网格,三角网格用于模拟复杂物体的表面,三角网格包括顶点、边和面三要素;
将构建的冠脉狭窄个体化三维模型中三角网格的某一个顶点P为中心,取所有与其相邻的顶点P1…Pn-1和边组成的一阶邻域结构,n表示三角网格的顶点个数,将每个顶点都移动到相邻顶点的平均位置,在冠脉狭窄个体化三维模型的平滑处理中设置法线平滑和顶点拟合的迭代次数和平滑因子,获取平滑后的冠脉狭窄个体化三维模型,定义获取平滑后的冠脉狭窄个体化三维模型为平滑冠脉狭窄三维模型;
步骤S104,对平滑冠脉狭窄三维模型进行曲面网格划分:
对平滑冠脉狭窄三维模型的冠状动脉连续体进行离散处理,利用几何单元逼近,基于变形协调条件完成有限元网格划分,输出曲面网格划分后的平滑冠脉狭窄三维模型;
有限元网格划分采用曲面离散映射法,曲面离散映射法划分网格的过程具体包括以下步骤:根据平滑冠脉狭窄三维模型的边界函数,采用映射矩阵,将边界映射到二维参数空间,在二维参数空间中对所选定的冠状动脉分支区域进行网格划分,将二维参数反向映射到物理空间形成曲面网格。
拟合单元工作过程具体包括以下步骤:
步骤S201,获取临床患者FFR实测数据,FFR实测数据包括冠脉狭窄远端压力和主动脉压力:将压力导丝传感器放置在冠脉狭窄处远端Qcm处,冠脉注射200ug硝酸甘油,等血压恢复后,静脉泵入ATP/腺苷,移动指针向右到最大充血状态的波形平台期同时开始记录,分别采样点测量的冠脉狭窄远端压力和主动脉压,实测FFR数据包括不同标记时刻的冠脉狭窄远端压力和主动脉压的实测值,记为Pc-d和Pa-d,Pc-d表示第d个时刻的冠脉狭窄远端压力,Pa-d表示第d个时刻的主动脉压;;
步骤S202,采用将实测冠脉狭窄远端压力和主动脉压进行曲线拟合,获得一个心动周期内不同入口边界条件下的流量曲线;
步骤S203,将血液密度设置为1060kg/m3,根据冠脉狭窄程度设定各个患者的血液粘稠度(男性4.57±0.73mPa·s,女性3.90±0.51mPa·s);将患者FFR拟合后数据作为冠脉周期脉动血流压力载入血管入口边界,出口设置为各患者FFR中拟合后的冠脉狭窄远端压力和主动脉压的周期平均值,实现血流心动周期性实际压差;
实现血流心动周期性实际压差,最终能量收敛,满足质量与能量守恒。
本实施例患者FFR数据拟合结果重度冠脉狭窄患者FFR拟合后数据曲线如图2所示;轻度冠脉狭窄患者FFR拟合后数据曲线如图3所示。
FFR数据进行曲线拟合具体包括以下步骤:
S202-1,对冠脉狭窄远端压力进行拟合,获得拟合后的冠脉狭窄远端压力;
利用预测的FFR线性模型函数对测量好的冠脉狭窄远端压力Pc进行拟合;以采样的血流时刻为横坐标,横坐标为d,d=1,2…D,d表示第d个采样时刻,D表示采样数量,以实测的冠脉狭窄远端压力为纵坐标,d为常数且d∈[1,D];拟合输出一个心动周期内冠脉狭窄远端压力曲线,输出动脉血流实时冠脉狭窄远端压力为:
f(pc_(d+1))=ω1Pc-12Pc-2+…ωdPc-d+b (1)
ωd表示基于FFR线性模型函数对冠脉狭窄远端压力进行拟合时多项式拟合的第d个系数,ωd属于常数;
式(1)为冠脉狭窄远端压力的FFR线性模型函数,f(pc_(d+1))表示第d+1个采样时刻冠脉狭窄远端压力的预测值;
冠脉狭窄远端压力的FFR线性模型函数输出预测的冠脉狭窄远端压力,预测的冠脉狭窄远端压力和冠脉狭窄远端压力的实测值构成损失函数,通过最小化损失函数求解ω和b并评估函数;Pc-d为一个心动周期内第d个采样时刻冠脉狭窄远端压力的实测值;
f(pc_(d+1))第d+1个采样时刻冠脉狭窄远端压力的预测值;
f(pc_(d+1))=ωTXPc+b (2)
式(2)为冠脉狭窄远端压力的FFR线性模型函数的向量表示,其中ω为向量系数;ω为ω1、ω2…ωd的向量表示,ωT表示ω的转置;
XPc是FFR线性模型函数(预测线性模型函数)的中Pc-1、Pc-2…Pc-d以行向量形式表示;
令xi=f(pc_d),yi=Pc-d,i=1,2…d…D;
/>
f(pc_d)表示第d个采样时刻冠脉狭窄远端压力的预测值,Pc-d表示d时刻冠脉狭窄远端压力的实测值;式(3)为冠脉狭窄远端压力的损失函数,用于衡量预测值f(xi)与实测值yi的差值,其中b*为b的偏置,ω*为ω向量的参数矩阵;
式(4)对两个参数ω和b求偏导,令式(4)等于0得到ω和b的最优解;b为函数FFR线性模型函数f(Pc-d)的参数,b*为参数b的偏置项;
i表示采样序号,i∈[1,D],D为采样数量;表示x1、x2…xD的均值;
式(5)和(6)为ω和b的最优解;
计算出ω和b的最优解后,代入公式(1),计算输出拟合后的冠脉狭窄远端压力;
S202-2,对主动脉压进行拟合,获得拟合后的主动脉压;
利用预测的FFR线性模型函数对测量好的主动脉压进行拟合;以采样的血流时刻为横坐标,横坐标为d,d=1,2…D,d表示第d个采样时刻,D表示采样数量,以实测的主动脉压为纵坐标,d为常数且d∈[1,D];拟合输出的是一个心动周期内主动脉压曲线,输出动脉血流实时主动脉压为:
f(Pa_(d+1))=ωa1Pa-1a2Pa-2+…ωadPa-d+ba (7)
式(7)中,ωad表示基于FFR线性模型函数对主动脉压进行拟合时多项式拟合的第d个系数,ωad属于常数,ba表示常数项;
式(7)为主动脉压的FFR线性模型函数,f(Pa_d)表示第d个采样时刻主动脉压的预测值;
主动脉压的FFR线性模型函数输出预测的主动脉压,预测的主动脉压和主动脉压的实测值构成损失函数,通过最小化损失函数求解ωa和ba并评估函数;Pa-d为一个心动周期内第d个采样时刻(第0.00秒、0.02秒、0.04秒、0.06秒……)主动脉压的实测值。
d表示采样时刻;d=1,2…D;
f(Pa_(d+1))=ωa TXPa+ba (8)
式(8)为主动脉压的FFR线性模型函数的向量表示,其中ωa为向量;ωa为ωa1、ωa2…ωad的向量表示,ωa T表示ωa的转置;
XPa是预测线性模型函数的未知数Pa-1、Pa-2…Pa-d以行向量形式表示;
令xa_i=f(Pa_d),ya_i=Pa-d,i=1,2…d…D;
f(Pa_d)表示第d个采样时刻主动脉压的预测值;Pa-d表示d时刻主动脉压的实测值;式(9)为损失函数,用于衡量预测值f(xi)与实测值yi的差值,其中ba *为ba的偏置,ωa *为ωa向量的参数矩阵;
式(10)对两个参数ωa和ba求偏导,令式(10)等于0得到ωa和ba的最优解;ba为函数f(Pa_d)的参数,b* a为参数ba的偏置项;
i表示采样序号,i∈[1,D],D为采样数量;
式(11)、(12)为ωa和ba的最优解;计算出ωa和ba的最优解后,代入公式(7),计算输出拟合后的主动脉压;
动力学参数计算单元工作过程具体包括以下步骤:
步骤S301,在所述血流心动周期性实际压差条件下,通过单向流体渗流Stokes方程仿真冠脉血流流动,仿真计算结果包括血流动力学参量心动周期内全过程演变;血流动力学参量包括瞬时流量、血管层流断面平均流速、湍流强度、血管壁面剪切应力;
Stokes方程的作用是是确保血流动力学计算中遵循能量守恒,动量守恒,质量守恒。在恒定压力条件下,单向流体渗流Stokes方程为:
stokes方程是用于描述流体运动的方程,看作是流体运动的牛顿第二定律。该步骤中通过stokes方程求得u,根据u即得出仿真计算所需的血流流速以及时间。
其中P指的是血管内血流总体压力;μ为血流粘度系数;u为血流速度向量;F为体积力向量;ρ为固定常量血流密度,T为绝对温度293.15K。I表示单位矩阵、是梯度算子(在空间各方向上的全微分);/>是对血流速度求梯度,是指向速度标量场增长最快的方向;指/>即梯度算子点乘uu函数;/>表示各方向上的矢量和;
根据牛顿粘性定律:血流粘度μ=τ·dv/dy,其中τ表示作用在流体平面上的剪应力,dv/dy是剪切应变率,斜率μ为血流粘度系数。
步骤S302,血流动力学参数包括瞬时流量、血管层流断面平均流速、湍流强度、血管壁面剪切应力;
如图4所示,在血管内取一个半径为r(正常冠状动脉r为0.1cm)、长度为L的圆柱形流体段。
1)瞬时流量Q:
K为稠度系数、L为血液流过的流体段、n为流变指数、R为血管内径;
Δp=2Lτ/r,其中τ为血管壁面剪切应力,r为流体段半径,Δp为血流流过流体段压差;
τ可由下方公式(18)得出:
得出血管层流断面平均流速V:
其中V为血管层流断面平均流速,μm为收缩期峰值时刻速度。
2)湍流强度B:
B=0.16*Re-1/8 (17)
B湍流强度;Re为雷诺数;
3)WSS(血管壁面剪切应力)τ:
其中τ表示WSS(血管壁面剪切应力),R为血管内径,γ代表血管横截面上任意点到圆心的距离,当γ=R时计算WSS。
实时追踪冠脉任意点的血流参数变化,量化不同狭窄几何构型耦合作用下冠脉血流动力学动态特性演变规律。
本申请做出了如下创新点:
(1)对血管进行动态模拟,在一个心动周期内患者临床实测FFR数据的拟合曲线作为fluent软件中模型的出入口条件,对冠脉进行动态模拟。这样设置边界,将会使冠脉内血流的流态更接近于实际情况。
(2)在冠状动脉建模时,采用了分割算法,分割出的冠状动脉比由普遍采用的区域生长算法重建的模型更符合人体血管的特征。
(3)实验获得了冠脉出现狭窄时间分辨率上的血流动力学分布,得到了量化的,能够更直接的描述血流的流变模式。
应用计算流体力学软件对冠状动脉患者进行数值模拟分析,以此获得冠状动脉的相关血流动力学参数,从而以无创的方式实现对患者的研究。相信此方法会在未来得到更多的应用。
在此处所提供的说明书中,说明了大量具体细节。然而,能够理解,本发明的实施例可以在没有这些具体细节的情况下被实践。在一些实例中,并未详细示出公知的方法、结构和技术,以便不模糊对本说明书的理解。
类似地,应当理解,为了精简本公开并帮助理解各个发明方面中的一个或多个,在上面对本发明的示例性实施例的描述中,本发明的各个特征有时被一起分组到单个实施例、图、或者对其的描述中。然而,并不应将该公开的方法解释成反映如下意图:即所要求保护的本发明要求比在每个权利要求中所明确记载的特征更多特征。更确切地说,如权利要求书所反映的那样,发明方面在于少于前面公开的单个实施例的所有特征。因此,遵循具体实施方式的权利要求书由此明确地并入该具体实施方式,其中每个权利要求本身都作为本发明的单独实施例。
本领域那些技术人员应当理解在本文所公开的示例中的设备的模块或单元或组间可以布置在如该实施例中所描述的设备中,或者可替换地可以定位在与该示例中的设备不同的一个或多个设备中。前述示例中的模块可以组合为一个模块或者此外可以分成多个子模块。
本领域那些技术人员可以理解,可以对实施例中的设备中的模块进行自适应性地改变并且把它们设置在与该实施例不同的一个或多个设备中。可以把实施例中的模块或单元或组间组合成一个模块或单元或组间,以及此外可以把它们分成多个子模块或子单元或子组间。除了这样的特征和/或过程或者单元中的至少一些是相互排斥之外,可以采用任何组合对本说明书(包括伴随的权利要求、摘要和附图)中公开的所有特征以及如此公开的任何方法或者设备的所有过程或单元进行组合。除非另外明确陈述,本说明书(包括伴随的权利要求、摘要和附图)中公开的每个特征可以由提供相同、等同或相似目的的替代特征来代替。
此外,本领域的技术人员能够理解,尽管在此所述的一些实施例包括其它实施例中所包括的某些特征而不是其它特征,但是不同实施例的特征的组合意味着处于本发明的范围之内并且形成不同的实施例。例如,在下面的权利要求书中,所要求保护的实施例的任意之一都可以以任意的组合方式来使用。
此外,所述实施例中的一些在此被描述成可以由计算机系统的处理器或者由执行所述功能的其它装置实施的方法或方法元素的组合。因此,具有用于实施所述方法或方法元素的必要指令的处理器形成用于实施该方法或方法元素的装置。此外,装置实施例的在此所述的元素是如下装置的例子:该装置用于实施由为了实施该发明的目的的元素所执行的功能。
这里描述的各种技术可结合硬件或软件,或者它们的组合一起实现。从而,本发明的方法和设备,或者本发明的方法和设备的某些方面或部分可采取嵌入有形媒介,例如软盘、CD-ROM、硬盘驱动器或者其它任意机器可读的存储介质中的程序代码(即指令)的形式,其中当程序被载入诸如计算机之类的机器,并被所述机器执行时,所述机器变成实践本发明的设备。
在程序代码在可编程计算机上执行的情况下,计算设备一般包括处理器、处理器可读的存储介质(包括易失性和非易失性存储器和/或存储元件),至少一个输入装置,和至少一个输出装置。其中,存储器被配置用于存储程序代码;处理器被配置用于根据该存储器中存储的所述程序代码中的指令,执行本发明的方法。
以示例而非限制的方式,计算机可读介质包括计算机存储介质和通信介质。计算机可读介质包括计算机存储介质和通信介质。计算机存储介质存储诸如计算机可读指令、数据结构、程序模块或其它数据等信息。通信介质一般以诸如载波或其它传输机制等已调制数据信号来体现计算机可读指令、数据结构、程序模块或其它数据,并且包括任何信息传递介质。以上的任一种的组合也包括在计算机可读介质的范围之内。
如在此所使用的那样,除非另行规定,使用序数词“第一”、“第二”、“第三”等等来描述普通对象仅仅表示涉及类似对象的不同实例,并且并不意图暗示这样被描述的对象必须具有时间上、空间上、排序方面或者以任意其它方式的给定顺序。
尽管根据有限数量的实施例描述了本发明,但是受益于上面的描述,本技术领域内的技术人员明白,在由此描述的本发明的范围内,可以设想其它实施例。此外,应当注意,本说明书中使用的语言主要是为了可读性和教导的目的而选择的,而不是为了解释或者限定本发明的主题而选择的。因此,在不偏离所附权利要求书的范围和精神的情况下,对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。对于本发明的范围,对本发明所做的公开是说明性的,而非限制性的,本发明的范围由所附权利要求书限定。

Claims (7)

1.一种基于FFR拟合的冠脉狭窄血流动力学模拟方法,其特征在于,具体包括以下步骤:
S1,建立冠脉狭窄个体化三维模型,基于三角网格顶点特征分解法对冠脉狭窄个体化三维模型进行平滑优化处理,采用曲面离散映射法对血管模型进行网格化处理,获得网格化平滑冠脉狭窄三维模型;
S2,获得临床FFR实测数据,对FFR实测数据进行拟合,加载到网格化平滑冠脉狭窄三维模型的入口边界,量化心动周期内的冠脉狭窄远端压力和主动脉压力输入;
S3,添加时间轴进行血流模拟计算,量化计算血管内任意点及心动周期内任意时刻的血流动力学参数;
步骤S2具体包括以下步骤:
步骤S201,获取临床患者FFR实测数据,FFR实测数据包括冠脉狭窄远端压力和主动脉压力的实测值,Pc-d表示第d个时刻的冠脉狭窄远端压力,Pa-d表示第d个时刻的主动脉压;
步骤S202,将实测冠脉狭窄远端压力和主动脉压进行曲线拟合,获得一个心动周期内不同入口边界条件下的流量曲线;
步骤S203,设置血液密度,根据冠脉狭窄程度设定各个患者的血液粘稠度;将FFR拟合后数据作为冠脉周期脉动血流压力载入血管入口边界,出口设置为各患者FFR中拟合后的冠脉狭窄远端压力的周期平均值,实现血流心动周期性实际压差;
步骤S202具体包括以下步骤:
S202-1,对冠脉狭窄远端压力进行拟合,获得拟合后的冠脉狭窄远端压力;
基于预测的FFR线性模型函数对冠脉狭窄远端压力进行拟合;以采样的血流时刻为横坐标,横坐标为d,d=1,2...D,d表示第d个采样时刻,D表示采样数量,以实测的冠脉狭窄远端压力为纵坐标;拟合输出一个心动周期内冠脉狭窄远端压力曲线,动脉血流冠脉狭窄远端压力的预测值为:
f(pc_(d+1))=ω1Pc-12Pc-2+…ωdPc-d+b (1)
ωd表示基于FFR线性模型函数对冠脉狭窄远端压力进行拟合时多项式拟合的第d个系数,ωd属于常数;
式(1)为冠脉狭窄远端压力的FFR线性模型函数,f(pc_(d+1))表示第d+1个采样时刻冠脉狭窄远端压力的预测值;
冠脉狭窄远端压力的FFR线性模型函数输出预测的冠脉狭窄远端压力,预测的冠脉狭窄远端压力和冠脉狭窄远端压力的实测值构成损失函数,通过最小化损失函数求解ω和b并评估函数;
f(pc_(d+1))=ωTXPc+b (2)
式(2)为冠脉狭窄远端压力的FFR线性模型函数的向量表示,其中ω为向量系数;ω为ω1、ω2…ωd的向量表示,ωT表示ω的转置;
XPc是FFR线性模型函数的中Pc-1、Pc-2…Pc-d以行向量形式表示;
令xi=f(pc_d),yi=Pc-d,i=1,2…d…D;
f(pc_d)表示第d个采样时刻冠脉狭窄远端压力的预测值;式(3)为冠脉狭窄远端压力的损失函数,用于衡量预测值xi与实测值yi的差值,其中b*为b的偏置,ω*为ω向量的参数矩阵;
式(4)对两个参数ω和b求偏导,令式(4)等于0得到ω和b的最优解;b为函数FFR线性模型函数的参数,b*为参数b的偏置项;
i表示采样序号,i∈[1,D];表示x1、x2…xD的均值;
式(5)和(6)分别为ω和b的最优解;
计算出ω和b的最优解后,代入公式(1),计算输出拟合后的冠脉狭窄远端压力;
S202-2,对主动脉压进行拟合,获得拟合后的主动脉压;
利用预测的FFR线性模型函数对测量好的主动脉压进行拟合;以采样的血流时刻为横坐标,横坐标为d,d=1,2…D,d表示第d个采样时刻,D表示采样数量,以实测的主动脉压为纵坐标;拟合输出的是一个心动周期内主动脉压曲线,输出动脉血流实时主动脉压为:
f(Pa_(d+1))=ωa1Pa-1a2Pa-2+…ωadPa-d+ba (7)
式(7)中,ωad表示基于FFR线性模型函数对主动脉压进行拟合时多项式拟合的第d个系数,ωad属于常数,ba表示常数项;
式(7)为主动脉压的FFR线性模型函数,f(Pa_(d+1))表示第d+1个采样时刻主动脉压的预测值;
主动脉压的FFR线性模型函数输出预测的主动脉压,预测的主动脉压和主动脉压的实测值构成损失函数,通过最小化损失函数求解ωa和ba并评估函数;
f(Pa_(d+1))=ωa TXPa+ba (8);
式(8)为主动脉压的FFR线性模型函数的向量表示,其中ωa为向量;ωa为ωa1、ωa2…ωad的向量表示,ωa T表示ωa的转置;
XPa是预测线性模型函数的中Pa-1、Pa-2…Pa-d以行向量形式表示;
令xa_i=f(Pa_d),ya_i=Pa-d,i=1,2...d...D;
f(Pa_d)表示第d个采样时刻主动脉压的预测值;式(9)为损失函数,用于衡量预测值xa_i与实测值ya_i的差值,其中ba *为ba的偏置,ωa *为ωa向量的参数矩阵;
式(10)对两个参数ωa和ba求偏导,令式(10)等于0得到ωa和ba的最优解;ba为函数f(Pa_d)的参数,b* a为参数ba的偏置项;
i表示采样序号,i∈[1,D],D为采样数量;
式(11)、(12)为ωa和ba的最优解;计算出ωa和ba的最优解后,代入公式(7),计算输出拟合后的主动脉压。
2.根据权利要求1所述的一种基于FFR拟合的冠脉狭窄血流动力学模拟方法,其特征在于,步骤S1具体包括以下步骤:
步骤S101,获取影像数据;
根据心电图获取有效相对静止的心脏时,基于快照冻结获得冠状动脉狭窄患者胸部增强CT医学DICOM格式图像;
步骤S102,建立冠脉狭窄个体化三维模型;
根据人体图像特征,将冠状动脉狭窄患者胸部增强CT医学DICOM格式图像剔除肌肉、软组织、空气和钙化成分,保留冠脉血管;根据冠状动脉相连的主动脉确定冠状动脉位置;使用多阈值适应性算法分割出目标血管,利用图像灰度特征分别计算出主动脉、左右冠脉所属的初始HU值区间,通过多次筛选目标血管部位HU值区间,将CT医学DICOM格式图像生成若干个独立阈值区间,得到分割掩膜像,基于分割掩膜像建立图像三视图,图像三视图包括横断位图像、矢状位图像和冠状位图像,在横断位图像上勾勒冠状动脉边缘,结合冠状位图像和矢状位图像,完成冠状动脉的提取,生成冠脉狭窄个体化三维模型;
步骤S103,对冠脉狭窄个体化三维模型进行平滑处理,具体包括以下步骤:
将冠脉狭窄个体化三维模型表面划分为三角网格,三角网格用于模拟复杂物体的表面,三角网格包括顶点、边和面三要素;
将构建的冠脉狭窄个体化三维模型中三角网格的某一个顶点P为中心,取所有与其相邻的顶点P1…Pn-1和边组成的一阶邻域结构,n表示三角网格的顶点个数,将每个顶点都移动到相邻顶点的平均位置,在冠脉狭窄个体化三维模型的平滑处理中设置法线平滑和顶点拟合的迭代次数和平滑因子,获取平滑后的冠脉狭窄个体化三维模型,定义获取平滑后的冠脉狭窄个体化三维模型为平滑冠脉狭窄三维模型;
步骤S104,对平滑冠脉狭窄三维模型进行曲面网格划分:
对平滑冠脉狭窄三维模型的冠状动脉连续体进行离散处理,利用几何单元逼近,基于变形协调条件完成有限元网格划分,输出曲面网格划分后的平滑冠脉狭窄三维模型。
3.根据权利要求2所述的一种基于FFR拟合的冠脉狭窄血流动力学模拟方法,其特征在于,
有限元网格划分采用曲面离散映射法,曲面离散映射法划分网格的过程具体包括以下步骤:根据平滑冠脉狭窄三维模型的边界函数,采用映射矩阵,将边界映射到二维参数空间,在二维参数空间中对所选定的冠状动脉分支区域进行网格划分,将二维参数反向映射到物理空间形成曲面网格。
4.根据权利要求1所述的一种基于FFR拟合的冠脉狭窄血流动力学模拟方法,其特征在于,
步骤S3计算具体包括以下步骤:
步骤S301,在所述血流心动周期性实际压差条件下,通过单向流体渗流Stokes方程仿真冠脉血流流动,仿真计算结果包括血流动力学参量心动周期内全过程演变;血流动力学参量包括瞬时流量、血管层流断面平均流速、湍流强度和血管壁面剪切应力;
单向流体渗流Stokes方程为:
通过stokes方程求得u;
其中P指血管内血流总体压力;μ为血流粘度系数;u为血流速度向量;F为体积力向量;ρ为固定常量血流密度,T为绝对温度293.15K;I表示单位矩阵,是梯度算子;/>是对血流速度求梯度,是指向速度标量场增长最快的方向;/>指/>表示梯度算子点乘uu函数;表示各方向上的矢量和;
步骤S302,血流动力学参数包括瞬时流量、血管层流断面平均流速、湍流强度和血管壁面剪切应力;
在血管内取一个半径为r、长度为L的圆柱形流体段;
1)瞬时流量Q:
K为稠度系数,L为血液流过的流体段长度,n为流变指数,R为血管内径;
Δp=2Lτ/r,其中τ为血管壁面剪切应力,r为流体段半径,Δp为血流流过流体段压差;
2)血管层流断面平均流速V:
其中V为血管层流断面平均流速,μm为收缩期峰值时刻速度;
3)湍流强度B:
B=0.16*Re-1/8 (17)
B为湍流强度;Re为雷诺数;
4)血管壁面剪切应力τ:
其中R为血管内径,γ代表血管横截面上任意点到圆心的距离,当γ=R时计算WSS。
5.一种基于FFR拟合的冠脉狭窄血流动力学模拟系统,其特征在于,具体包括模型建立单元、拟合单元和动力学参数计算单元;
模型建立单元建立冠脉狭窄个体化三维模型,基于三角网格顶点特征分解法对冠脉狭窄个体化三维模型进行平滑优化处理,采用曲面离散映射法对血管模型进行网格化处理,获得网格化平滑冠脉狭窄三维模型;
拟合单元融合临床FFR实测数据,对FFR实测数据进行拟合,加载到网格化平滑冠脉狭窄三维模型的入口边界,量化心动周期内的时序血流压力输入;
动力学参数计算单元添加时间轴进行血流模拟计算,量化计算血管内任意点及心动周期内任意时刻的血流动力学参数;
拟合单元工作过程具体包括以下步骤:
步骤S201,获取临床患者FFR实测数据,FFR实测数据包括冠脉狭窄远端压力和主动脉压力,记为Pc-d和Pa-d,Pc-d表示第d个时刻的冠脉狭窄远端压力的实测值,Pa-d表示第d个时刻的主动脉压的实测值;
步骤S202,采用将实测冠脉狭窄远端压力和主动脉压进行曲线拟合,获得一个心动周期内不同入口边界条件下的流量曲线;
步骤S203,设置血液密度,根据冠脉狭窄程度设定各个患者的血液粘稠度;将患者FFR拟合后数据作为冠脉周期脉动血流压力载入血管入口边界,出口设置为各患者FFR中拟合后的冠脉狭窄远端压力和主动脉压的周期平均值,实现血流心动周期性实际压差;
FFR数据进行曲线拟合具体包括以下步骤:
S202-1,对冠脉狭窄远端压力进行拟合,获得拟合后的冠脉狭窄远端压力;
利用预测的FFR线性模型函数对测量好的冠脉狭窄远端压力Pc进行拟合;以采样的血流时刻为横坐标,横坐标为d,d=1,2...D,d表示第d个采样时刻,D表示采样数量,以实测的冠脉狭窄远端压力为纵坐标;拟合输出一个心动周期内冠脉狭窄远端压力曲线,输出动脉血流的冠脉狭窄远端压力的预测值为:
f(pc_(d+1))=ω1Pc-12Pc-2+…ωdPc-d+b (1)
ωd表示基于FFR线性模型函数对冠脉狭窄远端压力进行拟合时多项式拟合的第d个系数,ωd属于常数;
式(1)为冠脉狭窄远端压力的FFR线性模型函数,f(pc_(d+1))表示第d+1个采样时刻冠脉狭窄远端压力的预测值;
冠脉狭窄远端压力的FFR线性模型函数输出预测的冠脉狭窄远端压力,预测的冠脉狭窄远端压力和冠脉狭窄远端压力的实测值构成损失函数,通过最小化损失函数求解ω和b并评估函数;Pc-d为一个心动周期内第d个采样时刻冠脉狭窄远端压力的实测值;
f(pc_(d+1))=ωTXPc+b (2)
式(2)为冠脉狭窄远端压力的FFR线性模型函数的向量表示,其中ω为向量系数;ω为ω1、ω2…ωd的向量表示,ωT表示ω的转置;
XPc是FFR线性模型函数的中Pc-1、Pc-2...Pc-d以行向量形式表示;
令xi=f(pc_d),yi=Pc-d,i=1,2...d...D;
f(pc_d)表示第d个采样时刻冠脉狭窄远端压力的预测值,Pc-d表示d时刻冠脉狭窄远端压力的实测值;式(3)为冠脉狭窄远端压力的损失函数,用于衡量预测值xi与实测值yi的差值,其中b*为b的偏置,ω*为ω向量的参数矩阵;
式(4)对两个参数ω和b求偏导,令式(4)等于0得到ω和b的最优解;b为函数FFR线性模型函数f(Pc_d)的参数,b*为参数b的偏置项;
i表示采样序号,i∈[1,D],D为采样数量;表示x1、x2…xD的均值;
式(5)和(6)为ω和b的最优解;
计算出ω和b的最优解后,代入公式(1),计算输出拟合后的冠脉狭窄远端压力;
S202-2,对主动脉压进行拟合,获得拟合后的主动脉压;
利用预测的FFR线性模型函数对测量好的主动脉压进行拟合;以采样的血流时刻为横坐标,横坐标为d,d=1,2...D,d表示第d个采样时刻,D表示采样数量,以实测的主动脉压为纵坐标;拟合输出的是一个心动周期内主动脉压曲线,输出动脉血流实时主动脉压为:
f(Pa_(d+1))=ωa1Pa-1a2Pa-2+…ωadPa-d+ba (7)
式(7)中,ωad表示基于FFR线性模型函数对主动脉压进行拟合时多项式拟合的第d个系数,ωad属于常数,ba表示常数项;
主动脉压的FFR线性模型函数输出预测的主动脉压,预测的主动脉压和主动脉压的实测值构成损失函数,通过最小化损失函数求解ωa和ba并评估函数;
d表示采样时刻;d=1,2...D;
f(Pa_(d+1))=ωa TXPa+ba (8)
式(8)为主动脉压的FFR线性模型函数的向量表示,其中ωa为向量;ωa为ωa1、ωa2…ωad的向量表示,ωa T表示ωa的转置;
XPa是预测线性模型函数的未知数Pa-1、Pa-2...Pa-d以行向量形式表示;
令xa_i=f(Pa_d),ya_i=Pa-d,i=1,2...d...D;
f(Pa_d)表示第d个采样时刻主动脉压的预测值;Pa-d表示d时刻主动脉压的实测值;式(9)为损失函数,用于衡量预测值xa_i与实测值ya_i的差值,其中ba *为ba的偏置,ωa *为ωa向量的参数矩阵;
式(10)对两个参数ωa和ba求偏导,令式(10)等于0得到ωa和ba的最优解;ba为函数f(Pa_d)的参数,b* a为参数ba的偏置项;
i表示采样序号,i∈[1,D],D为采样数量;表示xa_1、xa_2…xa_D的均值;
式(11)、(12)为ωa和ba的最优解;计算出ωa和ba的最优解后,代入公式(7),计算输出拟合后的主动脉压。
6.根据权利要求5所述的一种基于FFR拟合的冠脉狭窄血流动力学模拟系统,其特征在于,模型建立单元工作过程包括以下步骤:
步骤S101,获取影像数据
根据心电图(ECG)获取有效相对静止的心脏时,基于快照冻结获得冠状动脉狭窄患者胸部增强CT医学DICOM格式图像;
步骤S102,建立冠脉狭窄个体化三维模型;
根据人体图像特征,将冠状动脉狭窄患者胸部增强CT医学DICOM格式图像剔除肌肉、软组织、空气和钙化成分,保留冠脉血管;根据冠状动脉相连的主动脉确定冠状动脉位置;使用多阈值适应性算法分割出目标血管,利用图像灰度特征分别计算出主动脉、左右冠脉所属的初始HU值区间范围,通过多次筛选目标血管部位HU值区间,将CT医学DICOM格式图像生成若干个独立阈值区间,得到分割掩膜像,基于分割掩膜像建立图像三视图,图像三视图包括横断位图像、矢状位图像和冠状位图像,在横断位图像上勾勒冠状动脉边缘,结合冠状位图像和矢状位图像,完成冠状动脉的提取,生成冠脉狭窄个体化三维模型;
步骤S103,对冠脉狭窄个体化三维模型进行平滑处理,具体包括以下步骤:
将冠脉狭窄个体化三维模型表面划分为三角网格,三角网格用于模拟复杂物体的表面,三角网格包括顶点、边和面三要素;
将构建的冠脉狭窄个体化三维模型中三角网格的某一个顶点P为中心,取所有与其相邻的顶点P1…Pn-1和边组成的一阶邻域结构,n表示三角网格的顶点个数,将每个顶点都移动到相邻顶点的平均位置,在冠脉狭窄个体化三维模型的平滑处理中设置法线平滑和顶点拟合的迭代次数和平滑因子,获取平滑后的冠脉狭窄个体化三维模型,定义获取平滑后的冠脉狭窄个体化三维模型为平滑冠脉狭窄三维模型;
步骤S104,对平滑冠脉狭窄三维模型进行曲面网格划分:
对平滑冠脉狭窄三维模型的冠状动脉连续体进行离散处理,利用几何单元逼近,基于变形协调条件完成有限元网格划分,输出曲面网格划分后的平滑冠脉狭窄三维模型;
有限元网格划分采用曲面离散映射法,曲面离散映射法划分网格的过程具体包括以下步骤:根据平滑冠脉狭窄三维模型的边界函数,采用映射矩阵,将边界映射到二维参数空间,在二维参数空间中对所选定的冠状动脉分支区域进行网格划分,将二维参数反向映射到物理空间形成曲面网格。
7.根据权利要求6所述的一种基于FFR拟合的冠脉狭窄血流动力学模拟系统,其特征在于,动力学参数计算单元工作过程具体包括以下步骤:
步骤S301,在所述血流心动周期性实际压差条件下,通过单向流体渗流Stokes方程仿真冠脉血流流动,仿真计算结果包括血流动力学参量心动周期内全过程演变;血流动力学参量包括瞬时流量、血管层流断面平均流速、湍流强度和血管壁面剪切应力;
单向流体渗流Stokes方程为:
通过stokes方程求得;
其中P指的是血管内血流总体压力;μ为血流粘度系数;u为血流速度向量;F为体积力向量;ρ为固定常量血流密度,公式(13)中T为绝对温度293.15K,I表示单位矩阵、是梯度算子;/>是对血流速度求梯度,是指向速度标量场增长最快的方向;/>指/> 表示梯度算子点乘uu函数;/>将表示各方向上的矢量和;
步骤S302,血流动力学参数包括瞬时流量、血管层流断面平均流速、湍流强度、血管壁面剪切应力;
在血管内取一个半径为r、长度为L的圆柱形流体段;
1)瞬时流量Q:
K为稠度系数、L为血液流过的流体段、n为流变指数、R为血管内径;
Δp=2Lτ/r,其中τ为血管壁面剪切应力,r为流体段半径,Δp为血流流过流体段压差;
得出血管层流断面平均流速V:
其中V为血管层流断面平均流速,μm为收缩期峰值时刻速度;
2)湍流强度B:
B=0.16*Re-1/8 (17)
B湍流强度;Re为雷诺数;
3)血管壁面剪切应力τ:
R为血管内径,γ代表血管横截面上任意点到圆心的距离,当γ=R时计算血管壁面剪切应力。
CN202211732598.4A 2022-12-30 2022-12-30 基于ffr拟合的冠脉狭窄血流动力学模拟方法及系统 Active CN116453697B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211732598.4A CN116453697B (zh) 2022-12-30 2022-12-30 基于ffr拟合的冠脉狭窄血流动力学模拟方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211732598.4A CN116453697B (zh) 2022-12-30 2022-12-30 基于ffr拟合的冠脉狭窄血流动力学模拟方法及系统

Publications (2)

Publication Number Publication Date
CN116453697A CN116453697A (zh) 2023-07-18
CN116453697B true CN116453697B (zh) 2023-12-19

Family

ID=87129069

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211732598.4A Active CN116453697B (zh) 2022-12-30 2022-12-30 基于ffr拟合的冠脉狭窄血流动力学模拟方法及系统

Country Status (1)

Country Link
CN (1) CN116453697B (zh)

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105078440A (zh) * 2014-05-09 2015-11-25 西门子公司 无创计算冠状动脉狭窄的血液动力学指标的方法和系统
CN106539622A (zh) * 2017-01-28 2017-03-29 北京欣方悦医疗科技有限公司 基于血流动力学分析的冠脉虚拟支架植入方法和系统
CN108992057A (zh) * 2018-06-05 2018-12-14 杭州晟视科技有限公司 一种确定冠状动脉血流储备分数ffr的方法和装置
CN111241759A (zh) * 2020-01-13 2020-06-05 北京工业大学 一种基于零维血流动力学模型的ffr快速计算系统模型
CN112185551A (zh) * 2020-10-10 2021-01-05 北京工业大学 一种基于深度学习预测冠状动脉狭窄阻力的系统与方法
CN114530252A (zh) * 2022-03-03 2022-05-24 北京理工大学 冠状动脉血流动力学模拟仿真方法及装置
CN114677492A (zh) * 2022-04-13 2022-06-28 浙江大学 冠状动脉微循环血流动力学仿真分析方法及装置
CN114947909A (zh) * 2021-02-26 2022-08-30 复旦大学 基于狭窄病变前后的血液流量比值计算ffr的方法、系统

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10373700B2 (en) * 2012-03-13 2019-08-06 Siemens Healthcare Gmbh Non-invasive functional assessment of coronary artery stenosis including simulation of hyperemia by changing resting microvascular resistance

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105078440A (zh) * 2014-05-09 2015-11-25 西门子公司 无创计算冠状动脉狭窄的血液动力学指标的方法和系统
CN106539622A (zh) * 2017-01-28 2017-03-29 北京欣方悦医疗科技有限公司 基于血流动力学分析的冠脉虚拟支架植入方法和系统
CN108992057A (zh) * 2018-06-05 2018-12-14 杭州晟视科技有限公司 一种确定冠状动脉血流储备分数ffr的方法和装置
CN111241759A (zh) * 2020-01-13 2020-06-05 北京工业大学 一种基于零维血流动力学模型的ffr快速计算系统模型
CN112185551A (zh) * 2020-10-10 2021-01-05 北京工业大学 一种基于深度学习预测冠状动脉狭窄阻力的系统与方法
CN114947909A (zh) * 2021-02-26 2022-08-30 复旦大学 基于狭窄病变前后的血液流量比值计算ffr的方法、系统
CN114530252A (zh) * 2022-03-03 2022-05-24 北京理工大学 冠状动脉血流动力学模拟仿真方法及装置
CN114677492A (zh) * 2022-04-13 2022-06-28 浙江大学 冠状动脉微循环血流动力学仿真分析方法及装置

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Time-resolved simulation of blood flow through left anterior descending coronary artery: effect of varying extent of stenosis on hemodynamics;Zhao, YH et al.;《BMC CARDIOVASCULAR DISORDERS》;第第23卷卷(第第1期期);第1-13页 *
胡桃夹综合征血流动力学的血管建模研究;张世坤 等;《中国全科医学》;第3409-3413页 *

Also Published As

Publication number Publication date
CN116453697A (zh) 2023-07-18

Similar Documents

Publication Publication Date Title
US20210358635A1 (en) Systems and methods for image processing to determine blood flow
CN110853029B (zh) 用于基于医学图像自动预测血流特征的方法、系统和介质
CN111192316B (zh) 用于动脉分析和评估的深度学习
US10162932B2 (en) Method and system for multi-scale anatomical and functional modeling of coronary circulation
EP3160335B1 (en) Apparatus for determining a fractional flow reserve value
US10299862B2 (en) Three-dimensional quantitative heart hemodynamics in medical imaging
EP2942006B1 (en) Method and system for non-invasive computation of hemodynamic indices for coronary artery stenosis
CN108511075B (zh) 一种非侵入式获取血流储备分数的方法和系统
Xiong et al. Simulation of blood flow in deformable vessels using subject‐specific geometry and spatially varying wall properties
US11871995B2 (en) Patient-specific modeling of hemodynamic parameters in coronary arteries
Moosavi et al. Numerical simulation of blood flow in the left ventricle and aortic sinus using magnetic resonance imaging and computational fluid dynamics
Karabelas et al. Towards a computational framework for modeling the impact of aortic coarctations upon left ventricular load
EP3820357B1 (en) Patient-specific modeling of hemodynamic parameters in coronary arteries
Auricchio et al. A clinically applicable stochastic approach for noninvasive estimation of aortic stiffness using computed tomography data
Taylor et al. Patient-specific modeling of blood flow in the coronary arteries
CN116453697B (zh) 基于ffr拟合的冠脉狭窄血流动力学模拟方法及系统
Beier et al. Vascular hemodynamics with computational modeling and experimental studies
Yan et al. A multi-dimensional CFD framework for fast patient-specific fractional flow reserve prediction
Soni et al. A comprehensive review on CFD simulations of left ventricle hemodynamics: numerical methods, experimental validation techniques, and emerging trends
EP4322848A1 (en) Tracking segmental movement of the heart using tensors

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant