人体血管的应力测量装置、方法及存储介质

1.本技术涉及医学成像技术领域，特别涉及一种人体血管的应力测量装置、方法及存储介质。


背景技术：

2.相关研究表明血管壁的病变及其发展成心肌梗死、中风等多种心血管并发症的直接原因，多数血管壁的病变会导致动脉力学性质的改变，因此在体测量动脉力学性质对心血管疾病的早期筛查和诊断具有重大意义。
3.应力在血管的生长与重构中扮演着重要角色，因为血管中存在多种形式的应力，包括壁面切应力、壁内的环向/轴向正应力等，动脉中的细胞(如平滑肌细胞和内皮细胞)可以感受应力，并且通过改变内稳态环境对应力的改变做出应激性反应。因此对于动脉应力进行无创、在体测量，也是临床和健康领域迫切需要发展的一项关键技术。
4.相关技术中，对于血管力学性质的表征方面，主要是对动脉的弹性性质进行表征，而无法对血管的非线性特征和粘弹性特征进行表征，具有一定的局限性，并且现有分析算法的数据精度较低、效率低下，不利于临床使用；而对于血管应力的测量方面，主要是对血管环向应力的测量，受限于测量手段的有限性和测试条件的复杂性，导致血压的体侧量稳定性不高且操作复杂，不易于广泛使用。


技术实现要素：

5.本技术提供一种人体血管的应力测量装置、方法及存储介质，以解决相关技术中无法实现对血管非线性和粘弹性性质的系统表征导致具有一定局限性；对血管环向应力的测量稳定性不高、可靠性低且操作复杂等问题。
6.本技术第一方面实施例提供一种人体血管的应力测量装置，包括：超声测量模块，用于对目标位置的血管沿长轴切面或横轴切面进行超声成像，得到超声影像，通过心电信号触发的声辐射力对血管壁进行激励，使得所述血管壁激发出导波信号，采集信号得到超声成像序列；数据处理模块，用于识别所述超声影像得到血管半径和血管壁厚度，提取所述超声成像序列中导波信号的时空速度场，从所述时空速度场获取导波的实际频散曲线，基于所述血管半径、所述血管壁厚度和所述实际频散曲线进行血管材料参数及应力的反演，得到血管的实际应力和材料性质。
7.可选地，所述数据处理模块进一步用于：利用预设算法处理所述超声成像序列得到质点速度场；当沿血管长轴切面成像时，沿血管上壁中心线的位置提取出轴向导波的时空速度场；当沿血管横轴切面成像时，沿血管横截面的圆形中心线位置提取出环向导波的时空速度场。
8.可选地，所述数据处理模块进一步用于：对所述时空速度场进行加窗处理，得到处理后的时空速度场；对所述得到处理后的时空速度场进行二维傅里叶变换得到频谱图，基于所述频谱图的极值点确定导波的频散曲线。
9.可选地，所述数据处理模块进一步用于：根据所述血管半径和所述血管壁厚度确定血管几何；根据预设力学模型确定导波与所述血管几何之间的理论频散曲线；通过优化算法迭代求解所述实际频散曲线和所述理论频散曲线，将每次迭代求解得到参数的平均值作为最优血管材料参数，基于所述最优血管材料参数计算得到血管的实际应力和材料性质，其中，优化算法包括牛顿迭代法、模拟退火法、遗传算法。
10.可选地，所述理论频散曲线包括环向导波频散曲线和/或轴向导波频散曲线，其中，所述环向导波频散曲线的计算公式为：
11.c＝ω/re(k)，
12.所述轴向导波频散曲线的计算公式为：
[0013][0014]
其中，c表示相速度，ω表示角频率，k表示为角波数，re(k)表示取k的实部；ca(ω)表示轴向导波的相速度ca随频率ω变化的频散关系，r表示血管半径，n表示环向波数；
[0015]
ω与k关系由以下公式计算得到：
[0016][0017]
其中，s1和s2由以下四次方程求解得到：
[0018][0019]
其中，c
p2
＝κ
p
/ρf，κ
p
为水的体积模量，ρf为水的密度，t＝(αa+γ+αc)/3，i表示单位虚数，当求解轴向导波频散时取α＝αa，当求解环向导波频散时取α＝αc，h＝h/2表示半壁厚，ρ为血管壁密度，αa，αc，γ和β为血管的切线刚度参数，g和τ为血管的粘弹性参数。
[0020]
可选地，所述实际应力包括轴向应力和环向应力，其中，所述轴向应力σa的计算公式为：
[0021]
σa＝α
a-γ
[0022]
所述环向应力σc的计算公式为：
[0023]
σc＝α
c-γ
[0024]
其中，αa表示血管的轴向增量刚度(单位pa)，αc表示环向增量刚度(单位pa)，γ表示另一增量刚度系数(单位pa)。
[0025]
可选地，所述优化算法的目标函数为：
[0026][0027]
其中，表示实验测得的相速度，表示理论预测的相速度，α表示血管的增量刚度(对于应用于轴向或者环向导波时，则分别对应αa或αc)，γ为另一增量刚度参数，g为血管的粘弹性的松弛模量，τ为松弛特征时间，fi表示实验测得的频率，r表示血管半径，h表示血管壁厚，n表示实验测量数据点个数。
[0028]
可选地，所述超声测量模块包括：超声探头，用于对目标位置的血管沿长轴切面进行超声成像；超声主机，所述超声主机包含射频接收、射频发射端和聚焦声辐射力端，用于产生并发射声辐射力，对血管壁进行激励，并接收激发出的导波信号；超声系统，用于心电信号触发所述超声主机产生声辐射力。
[0029]
本技术第二方面实施例提供一种人体血管的应力测量方法，包括以下步骤：对目标位置的血管沿长轴切面或横轴切面进行超声成像，得到超声影像，通过心电信号触发的声辐射力对血管壁进行激励，使得所述血管壁激发出导波信号，采集信号得到超声成像序列；识别所述超声影像得到血管半径和血管壁厚度，提取所述超声成像序列中导波信号的时空速度场，从所述时空速度场获取导波的实际频散曲线，基于所述血管半径、所述血管壁厚度和所述实际频散曲线进行血管材料参数及应力的反演，得到血管的实际应力和材料性质。
[0030]
可选地，所述提取所述超声成像序列中导波信号的时空速度场，包括：利用预设算法处理所述超声成像序列得到质点速度场；当沿血管长轴切面成像时，沿血管上壁中心线的位置提取出轴向导波的时空速度场；当沿血管横轴切面成像时，沿血管横截面的圆形中心线位置提取出环向导波的时空速度场。
[0031]
可选地，所述从所述时空速度场获取导波的实际频散曲线，包括：对所述时空速度场进行加窗处理，得到处理后的时空速度场；对所述得到处理后的时空速度场进行二维傅里叶变换得到频谱图，基于所述频谱图的极值点确定导波的频散曲线。
[0032]
可选地，所述基于所述血管半径、所述血管壁厚度和所述实际频散曲线进行血管材料参数及应力的反演，得到血管的实际应力和材料性质，包括：根据所述血管半径和所述血管壁厚度确定血管几何；根据预设力学模型确定导波与所述血管几何之间的理论频散曲线；通过优化算法迭代求解所述实际频散曲线和所述理论频散曲线，将每次迭代求解得到参数的平均值作为最优血管材料参数，基于所述最优血管材料参数计算得到血管的实际应力和材料性质，其中，优化算法包括牛顿迭代法、模拟退火法、遗传算法。
[0033]
可选地，所述理论频散曲线包括环向导波频散曲线和/或轴向导波频散曲线，其中，所述环向导波频散曲线的计算公式为：
[0034]
c＝ω/re(k)，
[0035]
所述轴向导波频散曲线的计算公式为：
[0036][0037]
其中，c表示相速度，ω表示角频率，k表示为角波数，re(k)表示取k的实部；ca(ω)表示轴向导波的相速度ca随频率ω变化的频散关系，r表示血管半径，n表示环向波数；
[0038]
ω与k关系由以下公式计算得到：
[0039][0040]
其中，s1和s2由以下四次方程求解得到：
[0041][0042]
其中，c
p2
＝κ
p
/ρf，κ
p
为水的体积模量，ρf为水的密度，t＝(αa+γ+αc)/3，i表示单位虚数，当求解轴向导波频散时取α＝αa，当求解环向导波频散时取α＝αc，h＝h/2表示半壁厚，ρ为血管壁密度，αa，αc，γ和β为血管的切线刚度参数，g和τ为血管的粘弹性参数。
[0043]
可选地，所述实际应力包括轴向应力和环向应力，其中，所述轴向应力σa的计算公式为：
[0044]
σa＝α
a-γ
[0045]
所述环向应力σc的计算公式为：
[0046]
σc＝α
c-γ
[0047]
其中，αa表示血管的轴向增量刚度(单位pa)，αc表示环向增量刚度(单位pa)，γ表示另一增量刚度系数(单位pa)。
[0048]
可选地，所述优化算法的目标函数为：
[0049][0050]
其中，表示实验测得的相速度，表示理论预测的相速度，α表示血管的增量刚度(对于应用于轴向或者环向导波时，则分别对应αa或αc)，γ为另一增量刚度参数，g为血管的粘弹性的松弛模量，τ为松弛特征时间，fi表示实验测得的频率，r表示血管半径，h表示血管壁厚，n表示实验测量数据点个数。
[0051]
本技术第三方面实施例提供一种计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行，以用于实现如上述实施例所述的人体血管的应力测量方法。
[0052]
由此，本技术至少具有如下有益效果：
[0053]
本技术实施例基于超声方法和力学原理，设计了用于在体测量血管力学性质和应力的硬件装置系统，并发展了相应的数据处理方法，能够对人体的浅表动脉的力学性质及应力进行无创、实时的测量且分析效率较高；从表征的血管力学性质角度而言，通过动脉导波弹性成像方法对动脉的力学性质进行表征，可以表征血管的粘弹性、超弹性、各向异性等力学性质，有望为临床检测提供更多的定量化指导；从血管应力测量的角度而言，避免了血压测量的技术难题，而对血管轴向应力的测量可行性较高且可操作性较强。
[0054]
本技术附加的方面和优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本技术的实践了解到。
附图说明
[0055]
本技术上述的和/或附加的方面和优点从下面结合附图对实施例的描述中将变得明显和容易理解，其中：
[0056]
图1为根据本技术实施例提供的一种人体血管的应力测量装置的示意图；
[0057]
图2为根据本技术实施例的沿血管长轴切面进行超声测试的示例图；
[0058]
图3为根据本技术实施例的沿血管横轴切面进行超声测试的示例图；
[0059]
图4为根据本技术实施例的完整的超声成像序列示意图；
[0060]
图5为根据本技术实施例的当沿血管长轴切面成像时，沿血管上壁中心线的位置提取出时空速度场的示例图；
[0061]
图6为根据本技术实施例的当沿血管横轴切面成像时，沿血管上壁中心线的位置提取出时空速度场的示例图；
[0062]
图7为根据本技术实施例的对时空速度场进行加窗处理的示意图；
[0063]
图8为根据本技术实施例的由二维傅里叶变换得到频谱图；
[0064]
图9为根据本技术实施例的通过寻找极值点的方式得到频散曲线图；
[0065]
图10为根据本技术实施例的通过辉度图像对血管半径进行测量的示意图；
[0066]
图11为根据本技术实施例的通过辉度图像对血管下壁的壁厚进行测量的示意图；
[0067]
图12为根据本技术实施例的轴向导波反演方法的示意图；
[0068]
图13为根据本技术实施例的环向导波反演方法的示意图；
[0069]
图14为根据本技术实施例的颈动脉在舒张压和收缩压时对应的轴向导波频散曲
线；
[0070]
图15为根据本技术实施例的测量得到的颈动脉处的轴向应力；
[0071]
图16为根据本技术实施例的一种人体血管的应力测量方法的流程图。
具体实施方式
[0072]
下面详细描述本技术的实施例，所述实施例的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的，旨在用于解释本技术，而不能理解为对本技术的限制。
[0073]
近年来，心血管系统疾病已成为全球人类健康的头号杀手，居各种死因首位，相关研究表明，血管壁的病变是发展成心肌梗死、中风等多种心血管并发症的直接原因，多数血管壁的病变会导致动脉力学性质的改变，因此在体测量动脉力学性质对心血管疾病的早期筛查和诊断具有重大意义。
[0074]
应力在血管的生长与重构中扮演着重要角色，因为血管中存在多种形式的应力，包括壁面切应力、壁内的环向/轴向正应力等，动脉中的细胞(如平滑肌细胞和内皮细胞)可以感受应力，并且通过改变内稳态环境对应力的改变做出应激性反应。研究表明，壁面切应力的减小对于动脉斑块的形成有重要影响，血管环向应力的增大对动脉中弹力纤维的破损与断裂有显著关系，应力的改变也是斑块破裂、主动脉夹层等心血管疾病发生的重要因素。因此对于动脉应力进行无创、在体测量，也是临床和健康领域迫切需要发展的一项关键技术。
[0075]
动态弹性成像作为一项新兴的成像技术，在血管检测领域已经展现出一些技术优势，其中研究人员开发出的针对血管的弹性成像方法称为动脉导波弹性成像方法，该技术由聚焦声辐射力在血管壁中激励出弹性波，该弹性波以导波形式在血管壁内传播，再通过超声平面波快速成像技术对其进行成像，从而获取波场信息，结合波场信息即可对血管的力学性质进行表征，而从波场信息分析得到血管力学性质需要建立合适的力学模型，其中经典剪切波模型认为波速的平方和材料剪切模量满足线性关系，但是由于血管壁的有限厚度和其粘弹性的材料特性，使得其波场具有频散的特性，另一方面，由于其受到血压波动的影响，材料的非线性行为也会对波场产生重要影响。
[0076]
现有技术中未能对上述问题进行系统考虑，因此通过动脉导波成像技术分析得到的血管力学性质就存在不同程度的误差和不准确性。因此，本技术实施例提出了一种新的数据分析方法，从而对血管的力学性质进行了更准确的表征。
[0077]
另一方面，对于血管应力的测量，相关技术主要集中在对于血管环向应力的测量方法上，受限于测量手段的有限性和测试条件的复杂性(在体、无创)，目前并没有很好的对于血管轴向应力的测量方法。
[0078]
因此，本技术实施例利用血管导波成像技术，开发了一套同时对血管力学性质和血管应力(包括轴向和环向应力)进行测量的方法，一方面可以表征血管力学性质，另一方面也实现了血管轴向应力的测量。
[0079]
下面参考附图描述本技术实施例的人体血管的应力测量装置、方法及存储介质。具体而言，图1为本技术实施例所提供的一种人体血管的应力测量装置的示意图。
[0080]
如图1所示，该人体血管的应力测量装置10包括：超声测量模块100和数据处理模
块200。
[0081]
其中，超声测量模块100用于对目标位置的血管沿长轴切面或横轴切面进行超声成像，得到超声影像，通过心电信号触发的声辐射力对血管壁进行激励，使得血管壁激发出导波信号，采集信号得到超声成像序列；数据处理模块200用于识别超声影像得到血管半径和血管壁厚度，提取超声成像序列中导波信号的时空速度场，从时空速度场获取导波的实际频散曲线，基于血管半径、血管壁厚度和实际频散曲线进行血管材料参数及应力的反演，得到血管的实际应力和材料性质。
[0082]
其中，声辐射力可由产生特定声波的扬声器构成，会在特定位置产生压力密度的差异，即超声波在压力之下会出现形状的改变，并模拟出一种虚构的形状和力度，在此不做具体限定。
[0083]
可以理解的是，本技术实施例中的超声测量模块主要负责血管导波成像技术的可编程激励和定制化采集；而数据处理模块负责图像提取和数据分析，用于得到血管的力学性质及实际应力水平，能够对人体的浅表动脉的力学性质及应力进行无创、实时的测量以及数据分析，可行性较高且可操作性强。
[0084]
在本技术实施例中，超声测量模块包括：超声探头、超声主机和超声系统。
[0085]
其中，超声探头用于对目标位置的血管沿长轴切面进行超声成像；超声主机包含射频接收、射频发射端和聚焦声辐射力端，用于产生并发射声辐射力，对血管壁进行激励，并接收激发出的导波信号；超声系统用于心电信号触发超声主机产生声辐射力。
[0086]
其中，超声探头可选用普通线阵探头(探头中心频率5～15mhz)、或者由多个单阵元构成的点阵探头，在此不做具体限定。
[0087]
可以理解的是，本技术实施例中超声探头可根据不同视角测试长轴切面或横轴切面进行成像；超声主机用于产生并发射声辐射力，对血管壁进行激励，并接收激发出的导波信号；超声系统用于心电信号触发超声主机产生声辐射力；从而使得超声测量模块实现辉度模式成像、可编程声辐射力激励和高帧率采集功能。
[0088]
在本技术实施例中，数据处理模块进一步用于：利用预设算法处理超声成像序列得到质点速度场；当沿血管长轴切面成像时，沿血管上壁中心线的位置提取出轴向导波的时空速度场；当沿血管横轴切面成像时，沿血管横截面的圆形中心线位置提取出环向导波的时空速度场。
[0089]
其中，预设算法可以是loupas算法或是kasai算法等，可根据用户的实际需求进行选定，在此不做具体限定。
[0090]
可以理解的是，本技术实施例利用算法对超声成像序列进行处理得到质点速度场，对于平面波采集技术可以得到成像面内的质点速度场，当沿血管长轴切面成像时，沿血管上壁中心线的位置提取出轴向导波的时空速度场；当沿血管横轴切面成像时，沿血管横截面的圆形中心线位置提取出环向导波的时空速度场；对不同的成像采集有不同的提取方式从而提取出不同的时空速度场，以便于后续从时空速度场获取导波的实际频散曲线。
[0091]
在本技术实施例中，数据处理模块进一步用于：对时空速度场进行加窗处理，得到处理后的时空速度场；对得到处理后的时空速度场进行二维傅里叶变换得到频谱图，基于频谱图的极值点确定导波的频散曲线。
[0092]
其中，加窗处理可以对于无限长或持续时间很长的信号，由于存储和延时的问题，
计算处理时必须先对其进行截取，选择有限长的序列进行处理，在此不做具体限定。
[0093]
可以理解的是，本技术实施例对时空速度场进行加窗处理，得到处理后的时空速度场，对其进行二维傅里叶变换得到频谱图，并基于频谱图的极值点确定导波的频散曲线，反映出相速度随频率变化的关系，从而在时空速度场中获取导波的频散关系。
[0094]
在本技术实施例中，数据处理模块进一步用于：根据血管半径和血管壁厚度确定血管几何；根据预设力学模型确定导波与血管几何之间的理论频散曲线；通过优化算法迭代求解实际频散曲线和理论频散曲线，将每次迭代求解得到参数的平均值作为最优血管材料参数，基于最优血管材料参数计算得到血管的实际应力和材料性质，其中，优化算法包括但不限于牛顿迭代法、模拟退火法、遗传算法、粒子群优化算法等。
[0095]
需要说明的是，优化算法包括牛顿迭代法、模拟退火法、遗传算法、粒子群优化算法等，可以理解为优化算法并不仅限于上述算法，可根据实际需求进行选定，在此不做具体限定。
[0096]
其中，预设力学模型可以是用户事先设定的力学模型，可根据实际需求进行选定，在此不做具体限定。
[0097]
其中，理论频散曲线包括环向导波频散曲线和/或轴向导波频散曲线，其中，环向导波频散曲线的计算公式为：
[0098]
c＝ω/re(k)，
[0099]
轴向导波频散曲线的计算公式为：
[0100][0101]
其中，c表示相速度，ω表示角频率，k表示为角波数，re(k)表示取k的实部；ca(ω)表示轴向导波的相速度ca随频率ω变化的频散关系，r表示血管半径，n表示环向波数；
[0102]
ω与k关系由以下公式计算得到：
[0103][0104]
其中，s1和s2由以下四次方程求解得到：
[0105][0106]
其中，c
p2
＝κ
p
/
ρf，κ
p
为水的体积模量，ρf为水的密度，t＝(αa+γ+αc)/3，i表示单位虚数，当求解轴向导波频散时取α＝αa，当求解环向导波频散时取α＝αc，h＝h/2表示半壁厚，ρ为血管壁密度，αa，αc，γ和β为血管的切线刚度参数，g和τ为血管的粘弹性参数。
[0107]
其中，实际应力包括轴向应力和环向应力，其中，轴向应力σa的计算公式为：
[0108]
σa＝α
a-γ
[0109]
环向应力σc的计算公式为：
[0110]
σc＝α
c-γ
[0111]
其中，αa表示血管的轴向增量刚度(单位pa)，αc表示环向增量刚度(单位pa)，γ表示另一增量刚度系数(单位pa)。
[0112]
其中，优化算法的目标函数为：
[0113][0114]
其中，表示实验测得的相速度，表示理论预测的相速度，α表示血管的增量刚度(对于应用于轴向或者环向导波时，则分别对应αa或αc)，γ为另一增量刚度参数，g为血管的粘弹性的松弛模量，τ为松弛特征时间，fi表示实验测得的频率，r表示血管半径，h表示血管壁厚，n表示实验测量数据点个数。
[0115]
可以理解的是，本技术实施例根据血管半径和血管壁厚度确定血管几何，并根据力学模型确定导波与血管几何之间的理论频散曲线，将实验频散曲线和理论模型进行拟合，得到最优血管材料参数，基于最优血管材料参数计算血管的实际应力，通过数据处理模型分析得到的数据准确性高、效率快、稳定性强。
[0116]
根据本技术实施例提出的人体血管的应力测量装置，超声测量模块主要负责血管导波成像技术的可编程激励和定制化采集；而数据处理模块负责图像提取和数据分析，用于得到血管的力学性质及实际应力水平；使得从表征的血管力学性质角度而言，通过动脉导波弹性成像方法对动脉的力学性质进行表征，可以表征血管的粘弹性、超弹性、各向异性等力学性质，有望为临床检测提供更多的定量化指导；从血管应力测量的角度而言，避免了血压测量的技术难题，而对血管轴向应力的测量可行性较高且可操作性较强。由此，解决了相关技术中无法实现对血管非线性和粘弹性性质的系统表征导致具有一定局限性；对血管环向应力的测量稳定性不高、可靠性低且操作复杂等问题。
[0117]
下面将结合图2至图15对人体血管的应力测量装置进行详细阐述，主要从硬件系统、具体实施步骤和数据后处理方法三个方面进行描述，具体地：
[0118]
1、硬件系统包括：超声测量模块、数据处理模块和显示与输出模块。
[0119]
其中，超声测量模块由超声主机和超声探头构成。而超声主机包含rf(radio frequency，射频)接收/发射端、以及聚焦声辐射力端；超声探头可选用普通线阵探头(探头中心频率5～15mhz)、或者由多个单阵元构成的点阵探头；超声系统可实现辉度模式成像、
可编程声辐射力激励和高帧率(》5khz)采集功能。该超声设备需配备有心电采集模块，心电模块可以和声辐射力激励模块实现同步触发。
[0120]
其中，数据处理模块主要由计算机及搭载其上的软件构成，主要负责处理来自超声模块的图像数据，并对数据进行分析，从而得到血管力学性质和应力。
[0121]
其中，显示与输出模块主要由显示器构成，以屏显方式将测得的参数输出给用户，对于多次测量可以以实时更新的形式将数据输出给用户。
[0122]
2、具体实施方式是以受试者的测量颈总动脉为例，需要强调的是，本技术实施例可以对任何浅表动脉(如桡动脉、肱动脉等)进行测量；具体步骤如下：
[0123]
步骤1，用超声探头对血管沿长轴切面或横轴切面进行成像，超声探头轻放在皮肤表面，避免挤压血管。需要注意调节超声探头的视角，使得血管保持在视图的正中位置。其中图2(a)所示为沿血管长轴切面进行超声测试，图3(d)所示为沿血管横轴切面进行超声测试。在测试过程中，被测对象同时戴有心电监控装置，实时输出的心电信号用于声辐射力的触发。
[0124]
步骤2，通过心电信号触发声辐射力激励模块，从而在血管壁中激励出导波，而心电信号的触发阈值选取为心电波形的qrs峰，并设置特定的延迟时间以实现指定心动周期时刻(如收缩末期，舒张末期等)的超声信号采集。
[0125]
声辐射力的激励方式包括但不限于单一声辐射力(如图2(b)、3(e)所示)、移动声辐射力(多点声辐射力依次激励，图2(c)、3(f)所示)、同时聚焦多点声辐射力等，聚焦位置应当在血管壁或其附近位置处，使得能够激发出血管导波。
[0126]
步骤3，随后执行成像采集模式，实现血管导波信号的采集。
[0127]
具体实现成像采集的技术包括但不限于：平面波快速采集技术、a-line回波采集技术(至少需要两个采样点，即沿血管壁路径上的两点)等。其中，采样范围应至少包含血管的上壁和下壁；采样频率应当至少高于5khz；采集得到的原始数据可以是i/q数据或rf数据。
[0128]
在一个具体实施例中，采用平面波成像技术进行成像采集，时间为4ms，采样帧率为10khz，共采样40帧。一个完整的超声成像序列如图4所示，包括声辐射力阶段，空白阶段(由超声设备的实际性能而定，也可省去该阶段)，以及成像采集阶段。
[0129]
3、数据处理方法：
[0130]
(1)从超声系统采集得到i/q数据或rf数据，通过算法对原始数据处理，得到质点速度场v
x
(x,z,t)，具体的实现算法包括但不限于：loupas算法、kasai算法等。
[0131]
从质点速度场中提取血管导波信号，对于不同的成像采集技术可以有不同的提取方法，核心是至少获取到血管壁的两个质点位置处的时间-速度曲线。图5和图6展示了具体的实施例，对于平面波采集技术可以得到成像面内的质点速度场，当沿血管长轴切面成像时，沿血管上壁中心线的位置提取出时空速度场(图5所示)；当沿血管横轴切面成像时，沿血管横截面的圆形中心线位置提取出时空速度场(图6所示)。
[0132]
(2)从时空速度场中获取导波的频散关系，实现方法包括但不限于：二维傅里叶变换法，小波变换法等。
[0133]
作为一个具体实施例，首先对时空速度场进行加窗处理(图7)，之后由二维傅里叶变换得到频谱图(图8)，再通过寻找极值点的方式得到频散曲线(图9)，该频散曲线反映了
相速度随频率变化的关系。
[0134]
(3)对血管几何进行测量，包括血管半径r和壁厚h，需保证成像血管几何所处的时刻与声辐射力激励时刻基本一致(时间差小于10ms量级)；获取几何图像的方法包括但不限于常规b模式成像、快速成像技术等；测量几何的图像处理算法包括但不限于是阈值法、双峰拟合法、边缘检测法等。
[0135]
作为一个具体的实施例，选取快速成像得到的40帧中的一帧(本例中选取第一帧)，通过iq原始数据计算出血管的辉度图，如图10和图11所示。通过辉度图像对血管半径进行测量(图10)，以及对血管下壁的壁厚进行测量(图11)。
[0136]
下面将通过具体算法计算得出血管材料参数及应力的表征方法，具体如下：
[0137]
1、建立力学模型，得到血管导波频散曲线和血管材料参数、血压、几何等的理论关系，具体包括轴向导波的理论频散关系，以及环向导波的理论频散关系。
[0138]
理论模型将频散曲线和材料参数、几何进行了关联，因此为表征血管的材料参数，以及应力水平提供了理论指导。理论频散关系如下：
[0139][0140]
并且s1和s2由以下四次方程求解得到：
[0141][0142]
其中，c
p2
＝κ
p
/ρf，κ
p
为水的体积模量(2.2gpa)，ρf为水的密度(1000kg/m3)，t＝(αa+γ+αc)/3
[0143]
，其中i表示单位虚数，
[0144]
其中，当求解轴向导波频散时取α＝αa，当求解环向导波频散时取α＝αc，h＝h/2表示半壁厚，ρ为血管壁密度(取1000kg/m3)，αa，αc，γ和β为血管的切线刚度参数(量纲[pa])，g和τ为血管的粘弹性参数，分别表示松弛模量(量纲[1])和松弛特征时间(量纲[s])，ω为角频率(＝2πf，f为频率[hz])，k为角波数。
[0145]
由(1)式和(2)式可计算得到ω-k关系，再由下式计算得到相速度c
[0146]
c＝ω/re(k)
ꢀꢀꢀ
(3)
[0147]
(3)式即为所求的频散关系(即相速度c和频率f的关系)。
[0148]
对于计算血管轴向导波频散曲线，还需用到如下的修正公式：
[0149][0150]
其中ca(ω)表示轴向导波的相速度ca随频率ω变化的频散关系，c(ω)由(3)式决定，r为血管半径，n为环向波数，在此处取n＝2。
[0151]
综上，由(3)式可以计算理论的环向导波频散曲线，由(4)式可以计算理论的轴向导波频散曲线。
[0152]
2、应力的测量原理由如下的理论关系给出，轴向应力σa及环向应力σc和材料参数满足如下关系：
[0153]
σa＝α
a-γ
ꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0154]
σc＝α
c-γ
ꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(6)
[0155]
其中，αa表示血管的轴向增量刚度(单位pa)，αc表示环向增量刚度(单位pa)，γ表示另一增量刚度系数(单位pa)。
[0156]
3、建立反演方法，将实验频散曲线和理论模型进行拟合，得到材料参数。具体的，对于轴向导波，可以得到材料参数集(αa,γ,g,τ)；对于环向导波，可以得到材料参数集(αc,γ,g,τ)，再由材料参数可以计算得到应力。
[0157]
该反演方法包括但不限于：传统优化算法(如牛顿迭代法、最速下降法等)，智能优化算法(遗传算法、蚁群算法、粒子群算法、模拟退火法等)，机器学习方法(卷积神经网络等)。
[0158]
特别的，本技术实施例需要保护如下的一个具体的反演算法，即一套基于遗传算法的反演方法。对于轴向导波，其反演方法如图12所示。实验测得数据fi，ci(i＝1,2，
…
n,n表示测量的频散数据的数据点个数，一般取n＝40即可)，以及血管半径r和壁厚h。通过遗传算法迭代求解，得到最优的未知参数(αa,γ,g,τ)，为提高反演参数的稳定性，我们将该迭代求解过程重复m次(本实施例中取m＝3)，再对三次求解得到的参数集取平均值，从而得到最终的血管材料参数。再由(5)式即可计算得到血管的轴向应力σa。对于环向导波，其反演方法如图13所示，通过遗传算法迭代求解，得到血管环向的材料参数集(αc,γ,g,τ)，并由(6)式计算得到血管的环向应力σc。
[0159]
需要说明的是，在环向导波的实验数据中不需要血管半径r，因此仅就环向导波方法而言，在测量实验数据过程中还可以略去测量血管半径的步骤。
[0160]
在用遗传算法进行迭代求解时，需要定义目标函数，当其达到最小值时即解得最优参数解。优度目标函数定义为：
[0161][0162]
其中表示实验测得的相速度，表示理论预测的相速度，对于轴向导波则该理论解由(4)式决定，且α取为αa；对于环向导波则该理论解由(3)式决定，且α取为αc。n表示实验测量数据点个数(可取n＝40)。实验频散曲线的频段选取在500hz至1500hz之间，因
此n＝40表示在该频段内获取40个等分的频率-相速度数据点。
[0163]
4、作为一个具体的实施例，图14和图15展示了该方法的可行性和稳定性。图14为颈动脉在舒张压和收缩压时对应的轴向导波频散曲线，将实验曲线和理论模型进行拟合可得材料参数，图14亦画出了拟合的理论频散曲线。图15为测量得到的颈动脉处的轴向应力，可以观察到轴向应力在心动周期内存在明显的变化。
[0164]
需要说明的是，前文描述的轴向导波成像技术或环向导波技术，在具体的实施例中可以执行任一个而不必同时执行。具体来说，若执行的是轴向导波成像技术，则对血管轴向的力学性质以及血管轴向应力进行测量；若执行的是环向导波成像技术，则对血管环向的力学性质以及血管环向应力进行测量；若同时执行两个成像技术，则可以同时对血管双向的力学性质以及应力进行测量。
[0165]
其次参照附图描述根据本技术实施例提出的人体血管的应力测量方法。
[0166]
图16是本技术实施例的人体血管的应力测量方法的流程图。
[0167]
如图16所示，该人体血管的应力测量方法，包括以下步骤：
[0168]
在步骤s101中，对目标位置的血管沿长轴切面或横轴切面进行超声成像，得到超声影像，通过心电信号触发的声辐射力对血管壁进行激励，使得血管壁激发出导波信号，采集信号得到超声成像序列。
[0169]
可以理解的是，本技术实施例对目标位置的血管沿长轴切面或横轴切面进行超声成像得到超声影像，并通过心电信号触发声辐射力对血管壁进行激励，使血管壁激发出导波信号，采集信号得到超声成像序列，从表征的血管力学性质角度而言，通过动脉导波弹性成像方法对动脉的力学性质进行表征，可以表征血管的粘弹性、超弹性、各向异性等力学性质，有望为临床检测提供更多的定量化指导。
[0170]
在步骤s102中，识别超声影像得到血管半径和血管壁厚度，提取超声成像序列中导波信号的时空速度场，从时空速度场获取导波的实际频散曲线，基于血管半径、血管壁厚度和实际频散曲线进行血管材料参数及应力的反演，得到血管的实际应力和材料性质。
[0171]
可以理解的是，本技术实施例通过识别超声影像得到血管半径和血管壁厚度，提取超声成像序列中导波信号的时空速度场，并从中获取导波的实际频散曲线，基于血管半径、血管壁厚度和实际频散曲线进行血管材料参数及应力的反演，得到血管的实际应力和材料性质，从血管应力测量的角度而言，避免了血压测量的技术难题，而对血管轴向应力的测量可行性较高且可操作性较强。
[0172]
在本技术实施例中，提取超声成像序列中导波信号的时空速度场，包括：利用预设算法处理超声成像序列得到质点速度场；当沿血管长轴切面成像时，沿血管上壁中心线的位置提取出轴向导波的时空速度场；当沿血管横轴切面成像时，沿血管横截面的圆形中心线位置提取出环向导波的时空速度场。
[0173]
在本技术实施例中，从时空速度场获取导波的实际频散曲线，包括：对时空速度场进行加窗处理，得到处理后的时空速度场；对得到处理后的时空速度场进行二维傅里叶变换得到频谱图，基于频谱图的极值点确定导波的频散曲线。
[0174]
在本技术实施例中，基于血管半径、血管壁厚度和实际频散曲线进行血管材料参数及应力的反演，得到血管的实际应力和材料性质，包括：根据血管半径和血管壁厚度确定血管几何；根据预设力学模型确定导波与血管几何之间的理论频散曲线；通过优化算法迭
代求解实际频散曲线和理论频散曲线，将每次迭代求解得到参数的平均值作为最优血管材料参数，基于最优血管材料参数计算得到血管的实际应力和材料性质，其中，优化算法包括牛顿迭代法、模拟退火法、遗传算法。
[0175]
在本技术实施例中，理论频散曲线包括环向导波频散曲线和/或轴向导波频散曲线，其中，
[0176]
环向导波频散曲线的计算公式为：
[0177]
c＝ω/re(k)，
[0178]
轴向导波频散曲线的计算公式为：
[0179][0180]
其中，c表示相速度，ω表示角频率，k表示为角波数，re(k)表示取k的实部；ca(ω)表示轴向导波的相速度ca随频率ω变化的频散关系，r表示血管半径，n表示环向波数；
[0181]
ω与k关系由以下公式计算得到：
[0182][0183]
其中，s1和s2由以下四次方程求解得到：
[0184][0185]
其中，c
p2
＝κ
p
/ρf，κ
p
为水的体积模量，ρf为水的密度，t＝(αa+γ+αc)/3，i表示单位虚数，当求解轴向导波频散时取α＝αa，当求解环向导波频散时取α＝αc，h＝h/2表示半壁厚，ρ为血管壁密度，αa，αc，γ和β为血管的切线刚度参数，g和τ为血管的粘弹性参数。
[0186]
在本技术实施例中，实际应力包括轴向应力和环向应力，其中，
[0187]
轴向应力σa的计算公式为：
[0188]
σa＝α
a-γ
[0189]
环向应力σc的计算公式为：
[0190]
σc＝α
c-γ
[0191]
其中，αa表示血管的轴向增量刚度(单位pa)，αc表示环向增量刚度(单位pa)，γ表示另一增量刚度系数(单位pa)。
[0192]
在本技术实施例中，优化算法的目标函数为：
[0193][0194]
其中，表示实验测得的相速度，表示理论预测的相速度，α表示血管的增量刚度(对于应用于轴向或者环向导波时，则分别对应αa或αc)，γ为另一增量刚度参数，g为血管的粘弹性的松弛模量，τ为松弛特征时间，fi表示实验测得的频率，r表示血管半径，h表示血管壁厚，n表示实验测量数据点个数。
[0195]
需要说明的是，前述对人体血管的应力测量装置实施例的解释说明也适用于该实施例的人体血管的应力测量方法，此处不再赘述。
[0196]
根据本技术实施例提出的人体血管的应力测量方法，对目标位置的血管沿长轴切面或横轴切面进行超声成像得到超声影像，并通过心电信号触发声辐射力对血管壁进行激励，使血管壁激发出导波信号，采集信号得到超声成像序列，从表征的血管力学性质角度而言，通过动脉导波弹性成像方法对动脉的力学性质进行表征，可以表征血管的粘弹性、超弹性、各向异性等力学性质，有望为临床检测提供更多的定量化指导；通过识别超声影像得到血管半径和血管壁厚度，提取超声成像序列中导波信号的时空速度场，并从中获取导波的实际频散曲线，基于血管半径、血管壁厚度和实际频散曲线进行血管材料参数及应力的反演，得到血管的实际应力和材料性质，从血管应力测量的角度而言，避免了血压测量的技术难题，而对血管轴向应力的测量可行性较高且可操作性较强。由此，解决了相关技术中无法实现对血管非线性和粘弹性性质的系统表征导致具有一定局限性；对血管环向应力的测量稳定性不高、可靠性低且操作复杂等问题。
[0197]
本技术实施例还提供一种计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行时实现如上的人体血管的应力测量方法。
[0198]
在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本技术的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不是必须针对的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任一个或n个实施例或示例中以合适的方式结合。此外，在不相互矛盾的情况下，本领域的技术人员可以将本说明书中描述的不同实施例或示例以及不同实施例或示例的特征进行结合和组合。
[0199]
此外，术语“第一”、“第二”仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”的特征可以明示或者隐含地包括至少一个该特征。在本技术的描述中，“n个”的含义是至少两个，例如两个，三个等，除非另有明确具体的限定。
[0200]
流程图中或在此以其他方式描述的任何过程或方法描述可以被理解为，表示包括一个或更n个用于实现定制逻辑功能或过程的步骤的可执行指令的代码的模块、片段或部
分，并且本技术的优选实施方式的范围包括另外的实现，其中可以不按所示出或讨论的顺序，包括根据所涉及的功能按基本同时的方式或按相反的顺序，来执行功能，这应被本技术的实施例所属技术领域的技术人员所理解。
[0201]
应当理解，本技术的各部分可以用硬件、软件、固件或它们的组合来实现。在上述实施方式中，n个步骤或方法可以用存储在存储器中且由合适的指令执行系统执行的软件或固件来实现。如，如果用硬件来实现和在另一实施方式中一样，可用本领域公知的下列技术中的任一项或他们的组合来实现：具有用于对数据信号实现逻辑功能的逻辑门电路的离散逻辑电路，具有合适的组合逻辑门电路的专用集成电路，可编程门阵列，现场可编程门阵列等。
[0202]
本技术领域的普通技术人员可以理解实现上述实施例方法携带的全部或部分步骤是可以通过程序来指令相关的硬件完成，所述的程序可以存储于一种计算机可读存储介质中，该程序在执行时，包括方法实施例的步骤之一或其组合。
[0203]
尽管上面已经示出和描述了本技术的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本技术的限制，本领域的普通技术人员在本技术的范围内可以对上述实施例进行变化、修改、替换和变型。