CN112515704B - 一种基于超声的血管硬度测量方法 - Google Patents

Abstract

本发明公开了一种基于超声的血管硬度测量方法，步骤包括：选择取样门并矫正基线；解调得到IQ信号；划分成为M个子门；计算每一个子门中得被求和的时间系列的频谱；确定频带有效区域并计算血流速度；计算频带范围的平均速度；计算血流量速度‑时间曲线；计算位移‑时间矩阵；计算血管横截面积‑时间曲线；划分心动周期并计算舒张末期血管直径；计算脉冲波速度；计算血管壁厚度；计算血管的弹性模量；本发明使用多取样门技术提高了计算血流量速度的精确度，避免了传统使用M模式数据计算血管直径的同步困难，克服了传统弹性成像只能得到不同区域相对软硬变化和对组织运动敏感等不适合血管硬度测量的问题，实现了局部血管硬度的精确测量。

Description

一种基于超声的血管硬度测量方法

技术领域

本发明属于血管硬度测量领域，特别涉及一种基于超声的血管硬度测量方法。

背景技术

目前，心血管疾病的死亡率长期在城乡居民总死亡中高居首位，心血管疾病的及时防治成为当前政府迫在眉睫的任务。当前中国的心血管疾病患病率以及死亡率都仍处于上升阶段，研究报告推断现阶段我国现有心血管病患者人数超过了2.9亿，其中高血压患者约为2.45亿，脑卒患者约为1300万,冠心病患者约为1100万,肺原性心脏病患者约为500万,心力衰竭患者约为450万,风湿性心脏病患者约为250万,先天性心脏病患者约为200万。国家统计数据显示：心血管疾病死亡率高于肿瘤及其他疾病，位居所有死亡原因的首位，其中每5例死亡中就有 2例死于心血管病。如何及时的发现和诊断心血管疾病，如何能够及时快速的诊断心血管疾病的工作已经迫在眉睫。在所有的医学诊断中，超声医学诊断对于心血管疾病诊断有着极其明显的优势。超声医学诊断***的血流检查，尤其是超声诊断***的超声脉冲波多普勒成像技术，可以非常有效及时的对各种血流或血管疾病做出准确的早期分析。鉴于超声设备的无伤害和实时的特性，使得超声医学检查能够成为心血管疾病筛查的第一步。由此可见超声成像技术的发展尤其是超声多普勒技术的发展对心血管疾病的预防和诊断有着重要的意义。

传统医疗超声***使用超声回波的信号观测组织结构，可以检测血栓、心肌等组织的病变和发展情况，另一方面，基于超声多普勒效应的血流成像模式也可以用伪彩的方式显示血管内血流的分布情况，但这些无法反映血管的硬化情况。血管的软硬度与诸多心血管疾病如高血压、糖尿病等的关系非常密切，已经成为评价心血管疾病的一项重要指标，例如动脉硬化是动脉的一种非炎症性病变，可使动脉管壁增厚、***，失去弹性、管腔狭窄，伴随着高血压以及受累器官缺血等病症，对于早期的动脉硬化患者，大多数患者几乎无任何临床症状。对于中期的动脉硬化患者，大多数患者都或多或少有心悸、胸痛、胸闷、头痛、头晕、四肢凉麻、四肢酸懒、跛行、视力降低、记忆力下降、失眠多梦等临床症状。因此通过检测血管硬度早期发现和干预血管硬化的进程，对于预防和治疗心血管疾病具有十分重要的意义。

近年来，以超声弹性成像、声辐射力成像以及剪切波速成像等显示组织弹性技术的发展，可以定性或定量的反映检测区域软硬度，其中超声弹性成像和声辐射力成像可以得到定性的相对应变比反映不同区域的相对软硬度，而剪切波速成像通过测定剪切波在水平方向的传播速度，进而计算得到剪切模量和杨氏模量，实现定量的硬度测量，但剪切波波速测定结果容易受组织运动的影响，一般临床上主要用于肝脏纤维化的分级诊断。颈动脉的弹性成像主要用于测量血管斑块在血压作用下的运动和形变，观测斑块的弹性分布情况。而针对血管的硬度测量一般基于脉冲波速度来计算。

临床上通常采集基于两点测量法测量血管的脉冲波速度。首先，采集两处血管(一般是颈动脉与股动脉，另有颈动脉与桡动脉、肱动脉与踝动脉等)的压力曲线(也有采集血流曲线、血管运动曲线等)。然后，计算两处血管压力曲线之间的时间延迟，即脉搏在血管中传播的时间差。两测量点间血管的长度通过体表测量得到，脉冲波速度由血管长度除以时间延迟得到。传统测量方法通常是测量血管***上相距较远的两点间的脉冲波传播速度，并且假设两测量点间的血管走向为直线，得到的是血管硬度的整体平均值或全局值。实际上，血管硬度在血管***内并不是均匀分布的；另外，一些血管病变(如动脉粥样硬化斑块、腹主动脉瘤等)也会显著地改变局部血管的硬度。因此，局部血管硬度的测量更能反映动脉壁的病理变化特征。

发明内容

针对现有技术中的上述不足，本发明提供的一种基于超声的血管硬度测量方法通过多采样门血流量速度估计、自适应位移估计以及基于频域方式的血管直径估计，得到更为精确的局部血管硬度的测量结果，并结合一些边缘检测方法，实现血管壁厚度的自动计算，让***更便捷的完成血管硬度测量。

为了达到上述发明目的，本发明采用的技术方案为：

一种基于超声的血管硬度测量方法，包括以下步骤：

S1、选择指定的取样门，确定取样门的大小和位置；

S2、采集步骤S1中取样门的射频信号，经过正交解调得到IQ信号；

S3、将取样门划分为M个子采样门；

S4、按照指定的频谱计算方法，计算步骤S3中每一个子采样门中的被求和的时间系列的频谱；

S5、对各个子采样门得到的频谱进行上下包络检测，确定有效频带区域，在有效频带区域内计算血流平均速度；

S6、根据步骤S5中得到的血流平均速度，计算血流量速度-时间曲线；

S7、根据步骤S2得到的IQ信号，采用自适应频率估计的自相关方法计算位移-时间矩阵；

S8、根据步骤S2中IQ信号的包络计算初始血管直径，再根据血管直径得到血管横截面 -时间曲线；

S9、根据步骤S8得到的血管横截面划分心动周期，并计算舒张末期血管直径；

S10、根据步骤S8得到的血管横截面计算脉冲波速度；

S11、根据IQ信号的包络采用搜索的方法提取血管壁边缘并计算血管壁厚度；

S12、根据所得舒张末期血管直径、脉冲波速度和血管壁厚度，计算血管的弹性模量。

本发明的有益效果为：使用多取样门技术提高了计算血流量速度的精确度，同时获得血管直径，进而得到血管横截面-时间曲线，避免了传统使用M模式数据计算血管直径的同步困难，由于D和M扫描不能同时进行，所以血管横截面-时间曲线与血流量速度-时间曲线不能完全同步，本发明专利多取样门技术很好的避免了这个问题。在血管运动估计中引入了自适应频率衰减估计，提高了位移估计的准确度。设计了一种基于频率的心动周期划分算法，提高了心动周期自动划分的可靠性。克服了传统弹性成像只能得到不同区域相对软硬变化和对组织运动敏感等不适合血管硬度测量的问题，实现了局部血管硬度的精确测量。

进一步地，步骤S5具体过程为：用能量积分的方式计算上下包络，计算上包络时，从下往上在速度方向做累加，得到累积曲线，累积曲线切线趋近于0处则为上包络位置；计算下包络时，从上往下在速度方向做累加，得到累积曲线，累积曲线切线趋近于0处则为下包络位置；

计算有效频带范围的平均速度作为一个子取样门的速度，即血流平均速度，计算公式如下：

其中，

表示有效频带范围的血流平均速度，up表示上包络位置，down表示下包络位置， v_i表示第i个频谱位置对应的速度，p_i表示频谱能量值。

进一步地，步骤S6中血流量速度-时间曲线的计算公式表示为：

其中，

表示血流量速度-时间曲线，k表示子采样门序号，r_k表示第k个子采样门上的血管半径，N′+N-1是总的子采样门个数，N和N′分别表示血管中心向上和向下两个方向的子采样门个数，α表示采样门与血管壁的夹角，t_j表示第j个时间，

表示第j个时间上处于第k个子采样门上的血管位置的速度；本公式为用线积分求一个圆周面的血流量，得到血流量速度-时间曲线。

该进一步方案的有益效果为：传统单采样门在测量血流量时会忽略流速的空间信息的不足。在血管中的不同位置处，血流速度的分布实际上是不同的。计算PWV需要的血流量速度是血管横截面处的血流速度在横截面的线积分。如果血管发生了病变(例如血管斑块)或硬度发生变化，血管内不同位置到的血流速度会发生相应的变化，如果用传统的一个平均速度去计算整个血管的血流量速度，误差会很大

进一步地，步骤S7中位移-时间矩阵计算公式包括：

其中：

上述公式中Δz表示血管上下壁的位移，u表示窗口半径，l表示沿采样门方向轴的位置，c表示超声波的声度，R(0,1)_l,t表示第t个时间沿采样门方向轴向第l个位置的时间偏移1的自相关估计值，h表示位移计算窗口内的偏移游标，Conj表示取复数的共轭，R_phase(0,1)_l,t表示 R(0,1)_l,t的弧度相位角，Imag表示取复数虚部，Real表示取复数实部，atan表示求反正切值，

表示自适应估计的频率，f_d表示信号解调频率，f_s表示信号采样频率。

该进一步方案的有益效果为：由于使用了足够大的采样门，血管壁的回波信号也包含其中，因此可以得到血管上下壁的运动信息，避免了传统上用其它单独扫描方式来计算血管直径与血流流量的同步困难问题。并且自适应频率估计克服了衰减等带来的信号频率改变对位移估计结果的影响。

进一步地，步骤S8中血管横截面积-时间曲线的计算过程为：

得到整个取样门区域的位移-时间矩阵后，选择若干条IQ信号的包络，从中心往两侧搜索，先寻找两个方向上的最大值，再在反方向一定范围内寻找，找到初始血管壁的位置，得到多个直径，取平均作为初始直径，计算公式表示为：

其中，D₀表示血管壁初始直径，L₀表示初始上壁位置，Y₀表示初始下壁位置，W表示采样门宽度，S表示整个采样门内点数，α表示采样门与血管壁的夹角。

进一步地，步骤S9具体包括以下分步骤：

S91、将血流量速度-时间曲线进行傅立叶变换，取正频率方向幅值最大的频率的倒数，作为初始周期大小T0；

S92、对血流量速度-时间曲线每隔T0搜索一个最大值，对每个最大值位置的间隔先做1×3 的中值滤波，再进行平均处理后得到最终的周期大小T，周期T的倒数为心率；

S93、在血流量速度-时间曲线上以各个最大值为起点，在周期T的范围内，分别向前后采用梯度下降法搜索，找到最优的一个周期划分，每个周期长度均值为T，各个搜索梯度均方差最小；

S94、根据血管横截面-时间曲线中峰值后的波谷，取其中5到10个平均点，即为舒张末期血管直径。

进一步地，步骤S10中脉冲波速度计算过程如下：

将各个周期对应的血流量速度波形和血管横截面波形按照分割好的周期进行匹配，得到血流量速度-血管横截面环图，脉冲波速度为各个周期环图线性部分的直线斜率，脉冲波速度计算公式如下：

其中，V_pw表示脉冲波速度，Vol表示血流量速度波形，A表示血管横截面波形；

将各个周期的脉冲波速度V_pw先进行排序，去除指定比例的结果，再对剩下的数据进行平均处理，得到最终的结果。

该进一步方案的有益效果为：通过排序去除了多个心动周期数据中偏差过大的，再经过均值处理使脉冲波速度技术更为稳定可靠。

进一步地，步骤S11包括以下分步骤：

S111、首先通过从中心位置向上搜索最大值得到目标初始血管上壁中心，或向下搜索最大值得到目标初始血管下壁中心，搜索半径不超过正常血管半径加血管壁厚度乘以可调系数；

S112、沿上下两个方向做梯度搜索，找到一定搜索范围内梯度大于设定阈值，作为血管壁上下边缘；血管上下壁的阈值不同，血管上壁阈值范围更广；

S113、取不同时间的包络重复步骤S111和S112，得到同一位置不同时刻多个血管壁上下边缘；

S114、计算步骤S113中每一个时刻得到的血管壁厚度，首先滤除偏大过大的值，然后进行平均处理，得到血管厚度；

进一步地，步骤S12中血管的弹性模量的计算公式表示为：

其中，E表示血管的弹性模量，D表示舒张末期血管的直径，V_pw表示脉冲波速度，H表示血管壁的厚度，ρ表示血液密度。

附图说明

图1为本发明基于超声的血管硬度测量方法流程图；

图2为本发明中上下包络线示意图；

图3为本发明中多采样门血流量速度示意图；

图4为本发明中血流量速度-时间曲线示意图；

图5为本发明中血管横截面-时间曲线示意图；

图6为本发明中单个周期血流量速度-时间曲线图；

图7为本发明中单个周期血管横截面-时间曲线示意图；

图8为本发明中血管上壁和下壁厚度示意图；

图9为本发明中血流量速度-血管横截面环图。

具体实施方式

下面对本发明发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

如图1所示，一种基于超声的血管硬度测量方法，包括以下步骤：

S1、用户选择范围大到可以覆盖整个血管的取样门；同时用户选择取样门内感兴趣区域以显示参考多普勒频谱，并校正基线到适当位置，以确定所选择的区域为包含血流且频谱不发生混叠；

S2、确定取样门大小和位置后，采集取样门的射频信号，经过正交解调得到IQ信号(DIQ)；

S3、并将取样门(IQ信号)划分为M个子门；

S4、按照不同的频谱方法，例如FFT、Capon和APES等，计算每一个子门中的被求和的时间系列的频谱；各个子门的频谱计算会在时间方向每移动D个采样点做一次，也就是说频谱间隔的时间Δt＝D/PRF，PRF为脉冲重复频率。

S5、对各个子门得到的频谱进行上下包络检测确定有效频带区域，在有效频带区域内计算血流速度，有效频带如图2所示。上下包络采用能量积分的方式求得，上包络计算是从下往上在速度方向做累加，得到累积曲线，累积曲线切线趋近于0处则为上包络位置；下包络计算是从上往下在速度方向做累加，得到累积曲线，累积曲线切线趋近于0处则为下包络位置。一般取切线斜率小于0.05；

计算有效频带范围的平均速度作为一个子取样门的速度，即血流平均速度；

本发明实施例中，步骤S5中的计算公式为：

其中，

表示有效频带范围的血流平均速度，up表示上包络位置，down表示下包络位置， v_i表示第i个频谱位置对应的速度，p_i表示频谱能量值。

S6、根据步骤S5中得到的血流平均速度，计算血流量速度-时间曲线；

本发明实施例中，步骤S6中血流量速度-时间曲线的计算公式表示为：

其中，

表示血流量速度-时间曲线，k表示子采样门序号，r_k表示第k个子采样门上的血管半径，N′+N-1是总的子采样门个数，N和N′分别表示血管中心向上和向下两个方向的子采样门个数，α表示采样门与血管壁的夹角，t_j表示第j个时间，

表示第j个时间上处于第k个子采样门上的血管位置的速度；默认血管是圆周的，用线积分求一个圆周面的血流量，图3显示了多取样门的面积分示意图，再除以频谱对应的时间间隔，最后得到如图4 所示的血流量速度-时间曲线。

S7、根据步骤S2得到的IQ信号，采用自适应频率估计的自相关方法计算位移-时间矩阵；

本发明实施例中，步骤S7中位移-时间矩阵计算公式包括：

其中：

上述公式中Δz表示血管上下壁的位移，u表示窗口半径，l表示沿采样门方向轴的位置，c表示声音速度，R表示复数，R(0,1)_l,t表示第t个时间沿采样门方向轴向第l个位置的时间偏移1 的自相关估计值，h表示位移计算窗口内的偏移游标，Conj表示取复数的共轭，R_phase(0,1)_l,t表示R(0,1)_l,t的相位角，Imag表示取复数虚部，Real表示取复数实部，atan表示求反正切值，

表示自适应估计的频率，f_d表示信号解调频率，f_s表示信号采样频率。

位移矩阵计算是基于整个大的取样门来进行的，因为大的取样门覆盖了血管上下壁，因此就可以得到血管上下壁的位移，本实施例考虑衰减等带来的信号频率改变，因此采用自适应频率估计的自相关方法进行位移估计。

S8、根据步骤S2中IQ信号的包络计算初始血管直径，再根据血管直径得到血管横截面 -时间曲线；

本发明实施例中，步骤S8中血管横截面积-时间曲线的计算过程为：

得到整个取样门区域的位移-时间矩阵后，选择若干条IQ信号的包络，从中心往两侧搜索，先寻找两个方向上的最大值，再在反方向一定范围内寻找，找到初始血管壁的位置，得到多个直径，取平均作为初始直径，计算公式表示为：

其中，D₀表示血管壁初始直径，L₀表示初始上壁位置，Y₀表示初始下壁位置，W表示采样门宽度，S表示整个采样门内点数，α表示采样门与血管壁的夹角。

下一时刻的直径D₁，通过初始上下壁位移叠加对应的位移值，得到当前时刻的上下壁位置，代入上述公式即可。圆周的横截面积等于πD²/4，血管横截面积-时间曲线如图5所示。

S8、根据步骤S8得到的血管横截面划分心动周期，并计算舒张末期血管直径；

在计算过程中需要采集多个心动周期的数据，得到如图4和图5所示的多个周期的血流量速度与血管横截面-时间曲线，计算脉冲波速度时需要将其划分成单个周期，这里采用频率域的方法辅助划分。

本发明实施例中，步骤S9包括以下分步骤：

S91、将血流量速度-时间曲线进行傅立叶变换，取正频率方向幅值最大的频率的倒数，作为初始周期大小T0；

S92、对血流量速度-时间曲线每隔T0搜索一个最大值，对每个最大值位置的间隔先做1×3 的中值滤波，例如血流量速度-时间曲线包含8个周期就有8个最大值，对这8个位置的间隔先做1×3的中值滤波；之后再进行平均处理后得到最终的周期大小T，周期T的倒数为心率；

S93、在血流量速度-时间曲线上以各个最大值为起点，在周期T的范围内，分别向前后采用梯度下降法搜索，找到最优的一个周期划分，每个周期长度均值为T，各个搜索梯度均方差最小；因为在本发明专利的方法里面血流量速度计算和血管横截面是基于同一DIQ数据，因此它们是时间同步的。如图6和图7所示分别为划分后单个周期血流量速度-时间曲线和血管横截面-时间曲线；

S94、最后计算舒张末期血管直径，舒张末期是血管横截面-时间曲线中峰值后的波谷，如图7所示的尾部，取5到10个平均点。

S10、计算脉冲波速度；

本发明实施例中，步骤S10中脉冲波速度计算过程如下：

脉冲波速度需要计算出无反射期时的线性增长部分，将各个周期对应的血流量速度波形和血管横截面波形按照分割好的周期进行匹配，得到如图9所示的血流量速度-血管横截面环图，脉冲波速度为各个周期环图线性部分的直线斜率，脉冲波速度计算公式如下：

其中，V_pw表示脉冲波速度，Vol表示血流量速度波形，A表示血管横截面波形；

由于需要计算的是无反射期时的线性增长部分，所以对各个周期环图线性部分求直线斜率即脉冲波速度V_pw。将各个周期的脉冲波速度V_pw先进行排序，去除20％～40％比例的结果，再对剩下的数据进行平均处理，得到最终的结果。例如10个周期得到10个脉冲波速度，排序后去掉前后各2个，剩下6个平均得到最终的结果。

对于每一个环图，求解直线斜率采用逐点增加线性回归的拟合方法，即选取初始若干个点进行线性回归拟合，如果拟合误差大于设定阈值，则向前增加一个点，再进行拟合，直到达到迭代次数上限，则停止拟合。将所有小于设定阈值情况下的斜率排序，取其中一定范围 (80％～100％)的斜率值平均，得到这个环图的拟合的最终斜率。

这里的初始点数参考步骤S9计算的心率来设定，例如根据正常人在心率为75时候的平均快速射血期为0.09s，如果频谱间隔时间为0.0021s，换算到点数42个点左右。这里选取 40％～60％的点作为初始点集合，例如20个点；然后进行迭代线性拟合，直到满足条件。拟合误差阈值归一化范围0.7～0.95。迭代次数可以由平均快速射血期对应的点数减去初始点数，例如这里42-20＝22。

S11、根据IQ信号的包络采用搜索的方法提取血管壁边缘并计算血管壁厚度；

本发明实施例中，步骤S11包括以下分步骤：

S111、首先通过从中心位置向上搜索最大值得到目标初始血管上壁中心，或向下搜索最大值得到目标初始血管下壁中心，因为血管上壁和下壁厚度相同所以向下搜索血管上壁也一样，这里以上壁为例，为了减少初始搜索和误差，这里的搜索半径不超过(正常血管半径加血管壁厚度)乘以可调系数，以颈动脉为例，正常半径大约3mm，血管壁小于1.0mm，可调系数范围为0.8～1.5；

S112、沿上下两个方向做梯度搜索，找到一定搜索范围内梯度大于设定阈值，作为血管壁上下边缘；血管上下壁的阈值不同，血管上壁阈值范围更广，因为血管上壁向下临近血管，边缘梯度会相对大一些(阈值范围，梯度归一化后0.8～1.0)，而血管上壁向上一般临近组织，边缘梯度会相对小一些(阈值范围，梯度归一化后0.6～1.0)；

S113、取不同时间的包络重复步骤S111和S112，得到同一位置不同时刻多个血管壁上下边缘，血管上壁和下壁厚度示意图如图8所示；

S114、计算步骤S113中每一个时刻得到的血管壁厚度，首先滤除偏大过大的值(例如采用1×5的中值滤波等技术)，然后进行平均处理，得到血管厚度。

还可以采用另外两种方法：

方法一、基于B模式图像的测量，用户在进入血管硬度测量前，在B图像上完成血管壁的厚度计算，因为B图像上的厚度测量为常规测量，可以基于活动轮廓边缘检测进行测量，也可以基于深度学习模型；也可以基于深度学习模型例如FAST-RCNN等，相对基于B图像的测量结果更为精确，但血管壁厚度是作为血管硬度测量的前置输入条件。

方法二、基于DIQ的包络采用基于深度学***均值转换为距离值。

S12、根据所得到的脉冲波速度、舒张末期血管直径和血管壁厚度，计算血管的弹性模量。

本发明实施例中，步骤S12中血管的弹性模量的计算公式表示为：

其中，E表示血管的弹性模量，D表示舒张末期血管的直径，V_pw表示脉冲波速度，H表示血管壁的厚度，ρ表示血液密度。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。

Claims (8)

1.一种基于超声的血管硬度测量方法，其特征在于，包括以下步骤：

S1、选择指定的取样门，确定取样门的大小和位置；

S2、采集步骤S1中取样门的射频信号，经过正交解调得到IQ信号；

S3、将取样门划分为M个子采样门；

S4、按照指定的频谱计算方法，计算步骤S3中每一个子采样门中的被求和的时间系列的频谱；

S5、对各个子采样门得到的频谱进行上下包络检测，确定有效频带区域，在有效频带区域内计算血流平均速度；

S6、根据步骤S5中得到的血流平均速度，计算血流量速度-时间曲线；

S7、根据步骤S2得到的IQ信号，采用自适应频率估计的自相关方法计算位移-时间矩阵，其中位移-时间矩阵计算公式表示为：

其中：

上述公式中Δz表示血管上下壁的位移，u表示窗口半径，l表示沿采样门方向轴的位置，c表示超声波的声度，R(0,1)_l,t表示第t个时间沿采样门方向轴向第l个位置的时间偏移1的自相关估计值，h表示位移计算窗口内的偏移游标，Conj表示取复数的共轭，R_phase(0,1)_l,t表示R(0,1)_l,t的弧度相位角，Imag表示取复数虚部，Real表示取复数实部，atan表示求反正切值，

表示自适应估计的频率，f_d表示信号解调频率，f_s表示信号采样频率；

S8、根据步骤S2中IQ信号的包络计算初始血管直径，再根据血管直径得到血管横截面-时间曲线；

S9、根据步骤S8得到的血管横截面划分心动周期并计算舒张末期血管直径；

S10、根据步骤S8得到的血管横截面计算脉冲波速度；

S11、根据IQ信号的包络采用搜索的方法提取血管壁边缘并计算血管壁厚度；

S12、根据所得舒张末期血管直径、脉冲波速度和血管壁厚度，计算血管的弹性模量。

2.根据权利要求1所述的一种基于超声的血管硬度测量方法，其特征在于，所述步骤S5具体过程为：用能量积分的方式分别计算上下包络，计算上包络时，从下往上在速度方向做累加，得到累积曲线，累积曲线切线趋近于0处则为上包络位置；计算下包络时，从上往下在速度方向做累加，得到累积曲线，累积曲线切线趋近于0处则为下包络位置；

计算有效频带范围的平均速度作为一个子取样门的速度，即血流平均速度，计算公式如下：

其中，

表示有效频带范围的血流平均速度，up表示上包络位置，down表示下包络位置，v_i表示第i个频谱位置对应的速度，p_i表示第i个频谱位置对应的频谱能量值。

3.根据权利要求2所述的一种基于超声的血管硬度测量方法，其特征在于，所述步骤S6的具体过程包括：首先假定所测血管为圆周，之后用线积分计算一个圆周面的血流量，表示为：

其中，

表示血流量速度-时间曲线，k表示子采样门序号，r_k表示第k个子采样门上的血管半径，N′+N-1是总的子采样门个数，N和N′分别表示血管中心向上和向下两个方向的子采样门个数，α表示采样门与血管壁的夹角，

表示第j个时间上处于第k个子采样门上的血管位置的速度。

4.根据权利要求1所述的一种基于超声的血管硬度测量方法，其特征在于，所述步骤S8中血管横截面积-时间曲线的计算过程为：

得到整个取样门区域的位移-时间矩阵后，选择若干条IQ信号的包络，从中心往两侧搜索，先寻找两个方向上的最大值，再在反方向指定范围内寻找，找到初始血管壁的位置，得到多个血管壁直径，取平均作为初始血管壁直径，公式表示为：

其中，D₀表示血管壁初始直径，L₀表示初始上壁位置，Y₀表示初始下壁位置，W表示采样门宽度，S表示整个采样门内点数，α表示采样门与血管壁的夹角；对于下一时刻的血管壁直径D₁，通过初始上下壁位移叠加对应的位移值，得到当前时刻的上下壁位置L₁和Y₁，将当前时刻的上下壁位置L₁和Y₁代入上述公式即可得出下一时刻的血管壁直径D₁，以此类推，找到梯度改变最大的位置，即认为是血管壁；

则血管的横截面积表示为：πD²/4，以此得到血管横截面积-时间曲线。

5.根据权利要求4所述的一种基于超声的血管硬度测量方法，其特征在于，所述步骤S9包括以下分步骤：

S91、将血流量速度-时间曲线进行傅立叶变换，取正频率方向幅值最大的频率的倒数，作为初始周期大小T0；

S92、对血流量速度-时间曲线每隔T0搜索一个最大值，对每个最大值位置的间隔先做1×3的中值滤波，再进行平均处理后得到最终的周期大小T，周期T的倒数为心率；

S93、在血流量速度-时间曲线上以各个最大值为起点，在周期T的范围内，分别向前后采用梯度下降法搜索，找到最优的一个周期划分，每个周期长度均值为T，各个搜索梯度均方差最小；

S94、根据血管横截面-时间曲线中峰值后的波谷，取其中5到10个平均点，即为舒张末期血管直径。

6.根据权利要求5所述的一种基于超声的血管硬度测量方法，其特征在于，所述步骤S10中脉冲波速度计算过程如下：

将各个周期对应的血流量速度波形和血管横截面波形按照分割好的周期进行匹配，得到血流量速度-血管横截面环图，脉冲波速度为各个周期环图线性部分的直线斜率，脉冲波速度表示为：

其中，V_pw表示脉冲波速度，Vol表示血流量速度波形，A表示血管横截面波形；

将各个周期的脉冲波速度V_pw先进行排序，去除指定比例的结果，再对剩下的数据进行平均处理，得到最终的结果。

7.根据权利要求6所述的一种基于超声的血管硬度测量方法，其特征在于，所述步骤S11包括以下分步骤：

S111、首先通过从中心位置向上搜索最大值得到目标初始血管上壁中心，或向下搜索最大值得到目标初始血管下壁中心，搜索半径不超过正常血管半径加血管壁厚度乘以可调系数；

S112、沿上下两个方向做梯度搜索，找到一定搜索范围内梯度大于设定阈值，作为血管壁上下边缘；血管上下壁的阈值不同，血管上壁阈值范围更广；

S113、取不同时间的包络重复步骤S111和S112，得到同一位置不同时刻多个血管壁上下边缘；

S114、计算步骤S113中每一个时刻得到的血管壁厚度，首先滤除偏大过大的值，然后进行平均处理，得到血管壁厚度。

8.根据权利要求7所述的一种基于超声的血管硬度测量方法，其特征在于，所述步骤S12中血管的弹性模量的计算公式表示为：

其中，E表示血管的弹性模量，D表示舒张末期血管的直径，V_pw表示脉冲波速度，H表示血管壁的厚度，ρ表示血液密度。

Images