技术领域
[0001] 本
发明属于磁共振成像计算机仿真领域。
背景技术
[0002] 心肌微循环是微动脉、毛细血管和小静脉之间的血液循环，是心肌细胞与血液进行物质交换的重要场所。心肌微循环灌注异常会引发相应心肌缺血等临床症状，严重威胁人类健康。体素内不相干运动磁共振成像(Intra-voxel Incoherent Motion Magnetic Resonance Imaging，IVIM MRI)技术是近年来新兴的一种微循环灌注成像手段，是在扩散加权磁共振成像技术
基础上发展而来的，其原理是利用流向随机的灌注血液在扩散敏感梯度作用下质子群
相位严重不相干的物理现象，生成多b值(衰减因子)磁共振衰减
信号，并最终从衰减信号中获得血液容积分数和血流速度等微循环灌注信息。由于IVIM MRI技术具有无需
显像剂以及
分辨率高等优点，目前国内外正在积极尝试将其应用于心肌微循环灌注的临床检测，其中，检测效果的评估尤为重要。
[0003] 成像技术检测效果评估方法包括活体实验定性评估和仿真实验定量评估。活体实验定性评估是在不明确被测活体组织真实生理参数的情况下，仅通过将活体组织成像结果与该组织普遍意义上的生理参数比对，粗略估计成像技术的检测效果。由于活体实验定性评估无法精确量化检测效果，所以需要与仿真实验定量评估相结合，才能获得全面、可靠的评估结论。传统的仿真实验定量评估先要制作能够模拟
生物组织生理特征的仿真实体，再通过统计实体成像结果与实体参数之间的误差，得到检测效果的定量评估结论。然而，由于心肌微循环的生理结构极为复杂，所以仿真实体的制作难度将明显增加。同时，制作理想的仿真实体需要精密的加工设备和复杂的制作工艺，这将显著增加评估成本和评估周期。基于上述原因，目前国内外仅能通过活体实验定性评估方法粗略估计心肌微循环灌注IVIM MRI的检测效果。由于缺乏精确量化的评估结论，严重阻碍了IVIM MRI技术在心肌微循环检测方面的临床应用。
发明内容
[0004] 本发明的目的是为了解决活体实验定性评估方法不能精确量化检测结果，以及传统的仿真实验定量评估方法评估难度大，评估成本高，评估周期长的问题，而提出了心肌微循环灌注体素内不相干运动磁共振成像仿真方法。
[0005] 本发明心肌微循环灌注体素内不相干运动磁共振成像仿真方法按下列步骤实现：
[0006] 一、采用网络避让方程(1)使将要生成的血管段朝向函数f1(x)为最小值的方向生长；
[0007]
[0008] 其网络避让方程(1)中x表示将要生成的血管段轴线的终点，xn表示已生成血管段的质心，N表示已生成血管的段数，βv表示衰减系数，Un表示避让权值，Un的表达式为：
[0009]
[0010] 其Un表达式(2)中Rn表示已生成的第n段血管的半径，Ln表示已生成的第n段血管的长度；
[0011] 采用边界避让方程(3)避免将要生成的血管段冲出组织边界，边界避让方程为：
[0012]
[0013] 其边界避让方程(3)中x表示将要生成的血管段轴线的终点，xb表示x在边界上的投影，βb表示衰减系数，Ub表示边界避让权值，Ub的表达式为：
[0014]
[0015] 其Ub表达式(4)中Un表示避让权值；
[0016] 采用最优分枝
角度方程(5)计算子段与母段的夹角；
[0017]
[0018] 其最优分枝角度方程(5)中θ1表示子段a与母段的夹角，r1表示子段a的半径，r2表示子段b的半径，r0表示母段的半径，基于θ1建立子段a的
流体分支约束方程(6)；
[0019]
[0020] 其流体分支约束方程(6)中Ul表示流体分支约束方程相对于网络避让方程(1)和边界避让方程(3)在最终的心肌微循环网络建模方程中的权重，β3表示抑制因子，v表示母段的方向矢量，xp表示将要生成的血管段轴线的起点；
[0021] 将网络避让方程(1)、边界避让方程(3)和流体分支约束方程(6)加和，建立心肌微循环网络建模方程(7)，心肌微循环网络建模方程(7)的表达式为：
[0022]
[0023] 二、模拟血流灌注，通过
水分子定向流动位移方程(8)计算血管内水分子定向流动的位移，水分子定向流动的位移方程为：
[0024]
[0025] 其水分子定向流动的位移方程(8)中Δp表示水分子流动位移， 表示血管轴线处的最大流速，r表示水分子到血管轴线的距离，R表示血管半径，τ表示单步行走时间；
[0026] 模拟血管外水分子的扩散运动，通过水分子扩散位移模值表达式(9)计算水分子扩散位移的大小，水分子扩散位移模值表达式(9)为：
[0027]
[0028] 其水分子扩散位移模值表达式(9)中Δx表示水分子单步行走的位移大小，m表示扩散空间维数，τ表示单步行走时间，D表示水分子自由扩散系数，q表示(0,1)之间的随机数；
[0029] 通过水分子扩散位移方向表达式(10)计算水分子的扩散位移方向，水分子扩散位移方向表达式(10)为：
[0030]
[0031] 其水分子扩散位移方向表达式(10)中 表示水分子扩散位移单位矢量在笛卡尔
坐标系x轴、y轴、z轴上的投影，θ表示扩散位移方向与z轴正方向的夹角，φ表示扩散位移在xOy平面上的投影与x轴正方向的夹角；
[0032] 三、基于扩散磁共振成像原理以及步骤二中水分子流动位移和水分子扩散位移，计算由水分子位移导致的相位离散，由水分子位移导致的相位离散公式(11)为：
[0033]
[0034] 其由水分子位移导致的相位离散公式(11)中 表示第i个水分子第j步行走导致的相位离散值，γ表示氢核的旋磁比， 表示扩散敏感梯度，δ表示扩散敏感梯度持续时间，表示第i个水分子第j步行走的位移；
[0035] 在双极梯度脉冲自旋回波序列作用下，设定双极梯度脉冲的时间间隔为Δ，水分子单步行走(包括流动和扩散)时间为τ，则水分子行走步数k＝Δ/τ(12)，在关闭双极梯度脉冲后，水分子i的位移导致的相位离散为Δ时间内水分子i每一步行走位移产生的相位离散之和，水分子i每一步行走位移产生的相位离散之和的表达式为 然后依据扩散磁共振成像原理，单位体素内所有水分子相位离散引起的磁共振信号衰减式(14)为：
[0036]
[0037] 其中S表示扩散敏感梯度作用下的磁共振信号，S0表示未施加扩散敏感梯度的磁共振信号，n表示单位体素内水分子的个数；
[0038] 由IVIM MRI中磁共振信号衰减表达式 和单位体素内所有水分子相位离散引起的磁共振信号衰减式(14)可得：
[0039]
[0040] 其中f表示血液容积分数，D表示扩散系数，D*表示伪扩散系数，b表示衰减因子，在双极梯度脉冲自旋回波序列中，b的表达式为：
[0041]
[0042] 四、基于由水分子位移导致的相位离散公式(11)、水分子i每一步行走位移产生的相位离散之和的表达式(13)和单位体素内所有水分子相位离散引起的磁共振信号衰减式(14)，计算不同 值对应的磁共振衰减信号S/S0，从而获得一系列S/S0的离散值，然后基于b的表达式(17)，计算不同 值对应的b值，从而获得一系列b的离散值，基于已获得的S/S0和b的离散值，最后对磁共振信号衰减表达式(15)进行非线性拟合，获得心肌微循环灌注模型的IVIM MRI仿真检测结果。
[0043] 本发明首先基于真实心脏微循环网络的几何和拓扑结构，利用网络避让
算法、边界避让算法及流体分支约束算法，结合计算机强大的运算能
力，建立虚拟心肌微循环网络模型；然后，基于心脏微循环灌注的
组织学知识，建立血管内血流灌注模型，基于蒙特卡洛法，建立血管外水分子扩散运动模型，进而建立心肌微循环灌注模型；再在计算机中模拟IVIMMRI机制，计算不同幅值的扩散敏感梯度对应的磁共振衰减信号；最后对衰减信号非线性拟合，获得心肌微循环灌注模型的IVIM MRI计算机仿真实验检测结果。并计算检测结果与心肌微循环灌注模型中各参数之间的误差，由此实现对IVIM MRI技术检测效果的定量评估。
[0044] 本发明心肌微循环灌注体素内不相干运动磁共振成像仿真方法包含以下有益效果：
[0045] (1)为心肌微循环灌注的IVIM MRI检测效果提供可靠的仿真定量评估结论，促进其临床应用；
[0046] (2)利用网络避让算法、边界避让算法及流体分支约束算法，结合计算机强大的运算能力，建立虚拟单位体素心肌微循环网络模型，无需考虑仿真实体的制作难度，并且能够最大程度地逼近心肌微循环网络真实的生理结构。
[0047] (3)在计算机中建立虚拟心肌微循环网络模型，无需精密的加工设备和复杂的制作工艺，相比传统的仿真实验评估方法，显著降低评估成本，缩短评估周期。
[0048] (4)由于建立心肌微循环灌注模型及模拟IVIM MRI机制均在计算机中进行，所以能够连续调整模型参数和成像参数( δ、Δ)，得到丰富的误差统计数据，相比只能获得少量误差统计数据的传统仿真实验定量评估方法，评估结果的准确性显著提高。
附图说明
[0049] 图1为心肌微循环灌注体素内不相干运动磁共振成像仿真方法的
流程图；
[0050] 图2为步骤一中网络避让算法的几何示意图；
[0051] 图3为将要生成的血管段轴线的终点在组织边界上的投影；
[0052] 图4为步骤一中流体分支约束算法的几何示意图，1—子段a，2—子段b，3—母段；
[0053] 图5为步骤二中心肌微循环网络中血流的Laminar流示意图；
[0054] 图6为步骤二中单个水分子在t1到t2时间段内扩散示意图；
[0055] 图7为
实施例一建立的500×500×500μm3体素内虚拟心肌微循环网络模型图；
[0056] 图8为实施例一在不同 值下对应的S/S0值和b值以及Levenberg–Marquardt的非线性拟合曲线。
具体实施方式
[0057] 具体实施方式一：本实施方式心肌微循环灌注体素内不相干运动磁共振成像仿真方法按下列步骤实现：
[0058] 一、采用网络避让方程(1)使将要生成的血管段朝向函数f1(x)为最小值的方向生长；
[0059]
[0060] 其网络避让方程(1)中x表示将要生成的血管段轴线的终点，xn表示已生成血管段的质心，N表示已生成血管的段数，βv表示衰减系数，Un表示避让权值，Un的表达式为：
[0061]
[0062] 其Un表达式(2)中Rn表示已生成的第n段血管的半径，Ln表示已生成的第n段血管的长度；
[0063] 采用边界避让方程(3)避免将要生成的血管段冲出组织边界，边界避让方程为：
[0064]
[0065] 其边界避让方程(3)中x表示将要生成的血管段轴线的终点，xb表示x在边界上的投影，βb表示衰减系数，Ub表示边界避让权值，Ub的表达式为：
[0066]
[0067] 其Ub表达式(4)中Un表示避让权值；
[0068] 采用最优分枝角度方程(5)计算子段与母段的夹角；
[0069]
[0070] 其最优分枝角度方程(5)中θ1表示子段a与母段的夹角，r1表示子段a的半径，r2表示子段b的半径，r0表示母段的半径，基于θ1建立子段a的流体分支约束方程(6)；
[0071]
[0072] 其流体分支约束方程(6)中Ul表示流体分支约束方程相对于网络避让方程(1)和边界避让方程(3)在最终的心肌微循环网络建模方程中的权重，β3表示抑制因子，v表示母段的方向矢量，xp表示将要生成的血管段轴线的起点；
[0073] 将网络避让方程(1)、边界避让方程(3)和流体分支约束方程(6)加和，建立心肌微循环网络建模方程(7)，心肌微循环网络建模方程(7)的表达式为：
[0074]
[0075] 二、模拟血流灌注，通过水分子定向流动位移方程(8)计算血管内水分子定向流动的位移，水分子定向流动的位移方程为：
[0076]
[0077] 其水分子定向流动的位移方程(8)中Δp表示水分子流动位移， 表示血管轴线处的最大流速，r表示水分子到血管轴线的距离，R表示血管半径，τ表示单步行走时间；
[0078] 模拟血管外水分子的扩散运动，通过水分子扩散位移模值表达式(9)计算水分子扩散位移的大小，水分子扩散位移模值表达式(9)为：
[0079]
[0080] 其水分子扩散位移模值表达式(9)中Δx表示水分子单步行走的位移大小，m表示扩散空间维数，τ表示单步行走时间，D表示水分子自由扩散系数，q表示(0,1)之间的随机数；
[0081] 通过水分子扩散位移方向表达式(10)计算水分子的扩散位移方向，水分子扩散位移方向表达式(10)为：
[0082]
[0083] 其水分子扩散位移方向表达式(10)中 表示水分子扩散位移单位矢量在笛卡尔坐标系x轴、y轴、z轴上的投影，θ表示扩散位移方向与z轴正方向的夹角，φ表示扩散位移在xOy平面上的投影与x轴正方向的夹角；
[0084] 三、基于扩散磁共振成像原理以及步骤二中水分子流动位移和水分子扩散位移，计算由水分子位移导致的相位离散，由水分子位移导致的相位离散公式(11)为：
[0085]
[0086] 其由水分子位移导致的相位离散公式(11)中 表示第i个水分子第j步行走导致的相位离散值，γ表示氢核的旋磁比， 表示扩散敏感梯度，δ表示扩散敏感梯度持续时间，表示第i个水分子第j步行走的位移；
[0087] 在双极梯度脉冲自旋回波序列作用下，设定双极梯度脉冲的时间间隔为Δ，水分子单步行走(包括流动和扩散)时间为τ，则水分子行走步数k＝Δ/τ(12)，在关闭双极梯度脉冲后，水分子i的位移导致的相位离散为Δ时间内水分子i每一步行走位移产生的相位离散之和，水分子i每一步行走位移产生的相位离散之和的表达式为 然后依据扩散磁共振成像原理，单位体素内所有水分子相位离散引起的磁共振信号衰减式(14)为：
[0088]
[0089] 其中S表示扩散敏感梯度作用下的磁共振信号，S0表示未施加扩散敏感梯度的磁共振信号，n表示单位体素内水分子的个数；
[0090] 由IVIM MRI中磁共振信号衰减表达式 和单位体素内所有水分子相位离散引起的磁共振信号衰减式(14)可得：
[0091]
[0092] 其中f表示血液容积分数，D表示扩散系数，D*表示伪扩散系数，b表示衰减因子，在双极梯度脉冲自旋回波序列中，b的表达式为：
[0093]
[0094] 四、基于由水分子位移导致的相位离散公式(11)、水分子i每一步行走位移产生的相位离散之和的表达式(13)和单位体素内所有水分子相位离散引起的磁共振信号衰减式(14)，计算不同 值对应的磁共振衰减信号S/S0，从而获得一系列S/S0的离散值，然后基于b的表达式(17)，计算不同 值对应的b值，从而获得一系列b的离散值，基于已获得的S/S0和b的离散值，最后对磁共振信号衰减表达式(15)进行非线性拟合，获得心肌微循环灌注模型的IVIM MRI仿真检测结果。
[0095] 本实施方式步骤一中利用网络避让算法，避免将要生成的血管段与已生成的血管相互交叠，同时促使将要生成的血管段向血管分布稀疏的空间延伸，网络避让算法中的βv表示衰减系数，降低βv可强化远处血管段对将要生成血管段的影响；边界避让算法中的βb也表示衰减系数，降低βb可强化远处边界对将要生成血管段的影响。步骤一建立虚拟心肌微循环网络模型时除了避免血管之间，血管与边界之间彼此交叠外，还要保证血管分支角度符合血液动力学。心肌微循环网络中血管通常以二叉树的形式产生新的分支，依据血液动力学，子段与母段的夹角应保证血液流经夹角时所受的粘性摩擦
应力最小。因此采用最优分枝角度方程。流体分支约束方程中β3表示抑制因子，β3越大，对非最优分支角度的抑制越大，反之亦然。
[0096] 本实施方式步骤二模拟血流灌注等同于模拟大量水分子的定向流动。心肌微循环网络中的血流以Laminar流为主，即水分子沿着与血管轴平行的方向作平滑直线运动，流速在血管轴线处最大，近壁处最小，血管内水分子的平均流速与最大流速之比等于0.5。
[0097] 具体实施方式二：本实施方式与具体实施方式一不同的是步骤一
中子段b与母段的夹角θ2按如下最优分枝角度方程计算：
[0098]
[0099] 具体实施方式三：本实施方式与具体实施方式一或二不同的是步骤四取15～20个不同的 值。
[0100] 具体实施方式四：本实施方式与具体实施方式一至三之一不同的是步骤四利用Levenberg–Marquardt(列文伯格–
马夸尔特)算法对磁共振信号衰减表达式(15)进行非线性拟合。
[0101] 实施例一：在500×500×500μm3大小的体素内建立心肌微循环网络模型。由于心肌微循环网络中小静脉不参与和心肌细胞之间的物质交换，并且小静脉中的血液流速很低，所以IVIM MRI技术对小静脉中的血流不敏感，因而本实施例只建立包含微动脉和毛细血管的心肌微循环网络模型。其中，微动脉和毛细血管的界定数据选用文献“Coronary microvascular dysfunction”；微动脉和毛细血管的半径、长度和段数数据选用文献“Morphometry of pig coronary arterial trees”和“Topology and dimensions of pig coronary capillary network”。在体素表面选取若干个随机点，作为微动脉和部分毛细血管进入体素的起点；
[0102] 一、采用网络避让方程(1)使将要生成的血管段朝向函数f1(x)为最小值的方向生长；
[0103]
[0104] 其网络避让方程(1)中x表示将要生成的血管段轴线的终点，xn表示已生成血管段的质心，N表示已生成血管的段数，βv表示衰减系数，Un表示避让权值，Un的表达式为：
[0105]
[0106] 其Un表达式(2)中Rn表示已生成的第n段血管的半径，Ln表示已生成的第n段血管的长度；
[0107] 采用边界避让方程(3)避免将要生成的血管段冲出组织边界，边界避让方程为：
[0108]
[0109] 其边界避让方程(3)中x表示将要生成的血管段轴线的终点，xb表示x在边界上的投影，βb表示衰减系数，Ub表示边界避让权值，Ub的表达式为：
[0110]
[0111] 其Ub表达式(4)中Un表示避让权值；
[0112] 采用最优分枝角度方程(5)计算子段与母段的夹角；
[0113]
[0114] 其最优分枝角度方程(5)中θ1表示子段a与母段的夹角，r1表示子段a的半径，r2表示子段b的半径，r0表示母段的半径，基于θ1建立子段a的流体分支约束方程(6)；
[0115]
[0116] 其流体分支约束方程(6)中Ul表示流体分支约束方程相对于网络避让方程(1)和边界避让方程(3)在最终的心肌微循环网络建模方程中的权重，β3表示抑制因子，v表示母段的方向矢量，xp表示将要生成的血管段轴线的起点；
[0117] 将网络避让方程(1)、边界避让方程(3)和流体分支约束方程(6)加和，建立心肌微循环网络建模方程(7)，心肌微循环网络建模方程(7)的表达式为：
[0118]
[0119] 式(7)的βv取2，βb取2，β3取1.5，Ul在生成微动脉血管段时等于Ub，在生成毛细血管段时等于40，循环执行式(7)建立心肌微循环网络模型，最终建成的心肌微循环网络模型中血液容积分数f(即体素内所有血管的容积和占体素容积的百分比)为14.45％；
[0120] 二、在已建立的心肌微循环网络模型中模拟血流灌注，通过水分子定向流动位移方程(8)计算血管内水分子定向流动的位移，水分子定向流动的位移方程为：
[0121]
[0122] 其水分子定向流动的位移方程(8)中Δp表示水分子流动位移， 表示血管轴线处的最大流速，r表示水分子到血管轴线的距离，R表示血管半径，τ表示单步行走时间，式(8)中， 取1.0mm/s，其对应的伪扩散系数D*为12.76×10-3mm2/s，τ取1.0ms；
[0123] 在已建立的心肌微循环网络模型中模拟血管外水分子的扩散运动，通过水分子扩散位移模值表达式(9)计算水分子扩散位移的大小，水分子扩散位移模值表达式(9)为：
[0124]
[0125] 其水分子扩散位移模值表达式(9)中Δx表示水分子单步行走的位移大小，m表示扩散空间维数，τ表示单步行走时间，D表示水分子自由扩散系数，q表示(0,1)之间的随机数，本实施例式(9)中，D取1.0×10-3mm2/s，m取3，τ取1.0ms；
[0126] 通过水分子扩散位移方向表达式(10)计算水分子的扩散位移方向，水分子扩散位移方向表达式(10)为：
[0127]
[0128] 其水分子扩散位移方向表达式(10)中 表示水分子扩散位移单位矢量在笛卡尔坐标系x轴、y轴、z轴上的投影，θ表示扩散位移方向与z轴正方向的夹角，φ表示扩散位移在xOy平面上的投影与x轴正方向的夹角；
[0129] 三、基于扩散磁共振成像原理以及步骤二中水分子流动位移和水分子扩散位移，计算由水分子位移导致的相位离散，由水分子位移导致的相位离散公式(11)为：
[0130]
[0131] 其由水分子位移导致的相位离散公式(11)中 表示第i个水分子第j步行走导致的相位离散值，γ表示氢核的旋磁比， 表示扩散敏感梯度，δ表示扩散敏感梯度持续时间，表示第i个水分子第j步行走的位移；
[0132] 在双极梯度脉冲自旋回波序列作用下，设定双极梯度脉冲的时间间隔为Δ，水分子单步行走(包括流动和扩散)时间为τ，则水分子行走步数k＝Δ/τ(12)，在关闭双极梯度脉冲后，水分子i的位移导致的相位离散为Δ时间内水分子i每一步行走位移产生的相位离散之和，水分子i每一步行走位移产生的相位离散之和的表达式为 然后依据扩散磁共振成像原理，单位体素内所有水分子相位离散引起的磁共振信号衰减式(14)为：
[0133]
[0134] 其中S表示扩散敏感梯度作用下的磁共振信号，S0表示未施加扩散敏感梯度的磁共振信号，n表示单位体素内水分子的个数；
[0135] 由IVIM MRI中磁共振信号衰减表达式 和单位体素内所有水分子相位离散引起的磁共振信号衰减式(14)可得：
[0136]
[0137] 其中f表示血液容积分数，D表示扩散系数，D*表示伪扩散系数，b表示衰减因子，在双极梯度脉冲自旋回波序列中，b的表达式为：
[0138]
[0139] 四、基于由水分子位移导致的相位离散公式(11)、水分子i每一步行走位移产生的相位离散之和的表达式(13)和单位体素内所有水分子相位离散引起的磁共振信号衰减式(14)，令双极梯度脉冲的时间间隔Δ等于50.0ms，扩散敏感梯度持续时间δ等于2.0ms，分别等于0、0.005、0.008、0.010、0.012、0.015、0.018、0.020、0.030、0.050、0.070、0.100、0.130、0.150、0.160、0.190、0.200、0.220、0.240(T/m)，计算磁共振衰减信号S/S0，并基于式(17)，计算b值，然后，基于获得的S/S0和b的离散值，利用Levenberg–Marquardt法对式(15)非线性拟合，最终获得心肌微循环灌注模型的IVIM MRI仿真检测结果，结果为：f＝
15.52％，D＝0.85×10-3mm2/s，D*＝9.82×10-3mm2/s。
[0140] 参照步骤一至三中，心肌微循环灌注模型的参数，f＝14.45％，D＝1.0×10-3mm2/s，D*＝12.76×10-3mm2/s，可得检测结果与心肌微循环灌注模型各参数之间的误差为：|Δf|＝1.07％，|D|＝0.15×10-3mm2/s，|ΔD*|＝2.94×10-3mm2/s。
[0141] 改变心肌微循环灌注模型参数 D以及成像参数Δ、δ和 的方向，能够获得更为丰富的误差统计数据，为心肌微循环灌注IVIM MRI的检测效果提供更加可靠的仿真定量评估结论。
