一种超声频谱完整心动周期标记的方法与流程

本发明涉及到医学领域的超声技术，特别涉及到一种超声频谱完整心动周期标记方法。

背景技术：

超声技术被广泛应用于医学成像和测量。在脉冲波多普勒系统测量血流速度的过程中，系统按照一定时间间隔重复发射脉冲信号，在连续两次发射之间的某一时刻接收散射回波，再通过解调、滤波后，多普勒信号频谱从几兆赫兹搬移到以零频为中心，带宽几千赫兹的音频范围内。对解调出的音频信号进行两路处理，一路声音播放，通过扬声器可以听到多普勒信号的声音，一路进行傅里叶变换显示为以时间为横轴，频率为纵轴的功率谱图。在临床上对血管疾病进行诊断时，经常会截取多普勒信号声谱图中的某一段进行分析，从中提取一些参数，比如平均流速，最大流速，最小流速，收缩期、舒张末期最大流速比，阻尼指数，脉动指数等。而计算所有这些参数都必需以心动周期内血流的峰值为基础，因此如何准确从包络曲线中估算出包络峰值对计算其他医学参数指数具有重要意义。

传统血流动力学参数标记方法是：操作者判断声谱图满足性能要求后，冻结该声谱图并手工描记出收缩期峰值流速、舒张期末峰流速等标记点，或者用户自行选择一正常的信噪比较强的心动周期，直接由系统显示出该区域最大包络值，简单搜索最大最小值即可计算，获得血流医学参数指标。随着计算机计算的高速发展，研究者提出了多种频谱参数自动计算的方法。专利US6050,948采用阈值处理的方法估算包络最大值，该方法计算简单，计算效率高，搜索也很简单，但是这种阈值法判断环节比较单一，在存在多峰状态时很容易引起误判；专利CN 200710074477.4中，可以有效避免单个心动周期内多个峰值引起的误判，同时计算简单，效率较高，但是对于杂波较多的情况，很容易引起误判；专利CN201210200603.7，效率相对较高、对正常流速包络曲线峰值能正常标记，同时能避免幅值较大杂波的干扰，但是该算法参数设置需要借助先验主观信息，这并不能满足所有特殊情况的峰值搜索情况。

图1展示了现有技术中一种多普勒信号处理流程，是目前常用的一种超声频谱参数实时自动计算方法，一般对于信噪比较低或者复杂的频谱分布情况提取较为困难，因此目前大多数相关的频谱参数实时自动计算方法对谱的质量要求相对较高，而如何有效提取复杂频谱的医学参数成为目前研究的方向。

技术实现要素：

针对现有技术的不足，本发明提出了一种超声频谱完整心动周期标记的方法，用于实时计算包络线每个心动周期内的峰值，即使对复杂频谱，也可准确计算包络线每个心动周期内的峰值，从而基于此可准确计算其他相关血流医学参数指标。

一种超声频谱完整心动周期标记的方法，其中，包含以下步骤：

1)预存储经过解调、滤波、DFT处理的t-Ts到t时刻时间内的声功率谱数据Spectral(i，j)以及对应的正向包络数据EnvlopP(i)、负向包络数据EnvlopN(i)；

其中：i表示时间方向谱线序号,t-Ts≤i≤t；

j表示频率方向谱线的点序号；

之后，进入步骤2)；

2)取经步骤1)存储的t-Ts时刻到t时刻的谱数据Spectral(i，j)以及对应的正向包络数据EnvlopP(i)、负向包络数据EnvlopN(i)，进行计算以确定当前声功率谱的状态分布，具体的计算确定步骤如下：

21)获取缓存区内所有的从t-Ts到t时刻时间内的声功率谱数据Spectral(i，j)，计算其平均值，再计算所有从t-Ts到t时刻时间内的声功率谱数据Spectral(i，j)的正向累加和、负向累加和，公式如下：

从t-Ts到t时刻时间内的声功率谱数据Spectral(i，j)的平均值为MeanSpectral(j)，具体计算公式如下：

$MeanSpectral\left(j\right)=\frac{\underset{i=t-T_s}{\overset{i=t}{\Σ}}Spectral(i,j)}{n};$

从t-Ts到t时刻时间内声功率谱数据Spectral(i，j)的正向累加和为SumSpectralP，具体计算公式如下：

$SumSpectralP=\underset{j=baseline}{\overset{j=n-1}{\Σ}}MeanSpectral\left(j\right);$

从t-Ts到t时刻时间内声功率谱数据Spectral(i，j)的负向累加和为SumSpectralN，具体计算公式如下：

$SumSpectralN=\underset{j=0}{\overset{j=basline-1}{\Σ}}MeanSpectral\left(j\right);$

其中:baseline表示量化后的基线位置，0≤baseline＜n；

n表示单根谱线采样点最大值；

Spectral(i，j)表示i时刻离散谱线数据；

22)计算缓存区内t-Ts到t时刻时间之间的正向包络数据EnvlopP(i)、负向包络数据EnvlopN(i)的正向累加和、负向累加和，公式如下：

所有谱线包络的正向累加和为：

$SumEnvlopP=\underset{i=t-Ts}{\overset{i=t}{\Σ}}\left(Envlop\right(i)-baseline+1);$

所有谱线包络的负向累加和为：

$SumEnvlopN=\underset{i=t-Ts}{\overset{i=t}{\Σ}}(baseline-Envlop(i\left)\right)$

其中：Envlop(i)表示从t-Ts时刻到t时刻之间谱线的包络值；

SumEnvlopP表示谱线包络正向累加和；

SumEnvlopN表示谱线包络负向累加和；

23)取步骤21)中频谱的正向累加和SumSpectralP、负向累加和SumSpectralN，取步骤22)中包络的正向累加和SumEnvlopP、负向累加和SumEnvlopN，分别比较判断频谱的正向累加和SumSpectralP与频谱的负向累加和SumSpectralN、包络的正向累加和SumEnvlopP与包络的负向累加和SumEnvlopN，以区分缓存区内t-Ts到t时刻之间正向包络EnvlopP(i)、负向包络EnvlopN(i)的分布状态，分布状态公式如下：

如果SumEnvlopP≥SumEnvlopN&SumSpectralP≥SumSpectralN，那么

Envlop(i)＝EnvlopP(i)，进入步骤3)；

如果SumEnvlopP＜SumEnvlopN&SumSpectralP＜SumSpectralN，那么

Envlop(i)＝n-EnvlopN(i)，进入步骤3)；

否则返回步骤1)，重新开始；

其中：&表示逻辑“与”运算；

n表示单根谱线采样点最大值，与步骤21)中一致；

Envlop(i)为确定状态的声功率谱包络，包括正向包络EnvlopP(i)或者负向包络EnvlopN(i)。

3)对步骤2)中已确定状态的声功率谱包络Envlop(i)进行降干扰处理，降干扰处理后的包络记为Envlop_F(i)，降干扰处理步骤如下：

31)首先对Envlop(i)去除毛刺干扰，去除毛刺干扰计算公式如下：

Envlop_M_F(i)＝Get2[Reshape[Envlop(i-1),Envlop(i),Envlop(i+1)]]；

Envlop_M_F(i)为去除毛刺干扰后的包络数据；

32)其次对Envlop_M_F(i)抑制高频干扰，公式如下：

Envlop_H_F(i)＝[Envlop_M_F(i-1)+Envlop_M_F(i)+Envlop_M_F(i+1)]/3；

Envlop_H_F(i)为抑制高频干扰后的包络数据；

33)最后对Envlop_H_F(i)凸出峰值，公式如下：

Envlop_F(i)＝Max[Envlop_H_F(i-1),Envlop_H_F(i),Envlop_H_F(i+1)]；

Envlop_F(i)为凸出峰值处理后的包络数据；

其中：Reshape[]函数表示对[]中三个数进行从小到大排序；

Get2[]函数表示取[]内序列的第2个值；

Max[]表示取[]中三个值的最大值；

4)分类确定步骤3)中降干扰后的包络数据Envlop_F(i)是否具有完整的波峰特征，其中分类确定步骤如下；

41)预设定包络判断阈值Envlop_sdth并存储；

42)将步骤3)中降干扰后的包络数据Envlop_F(i)进行N等分，对每一等分计算其标准差，标准差公式如下：

$Envlop\_sd\_x\left(m\right)=\sqrt{\frac{1}{M}\underset{k=1}{\overset{k=M}{\Σ}}{\left(Envlop\right(k)-u)}^2}$

其中：M为每等分的采样点数；

u为每一等分的平均值；

m＝1、2、3……N；

k＝1、2、3……M；

Envlop(k)表示对包络数据Envlop_F(i)进行N等分后的第k个等分序列表达式；

43)计算步骤42)中Envlop_sd_x(m)的算术平均数std_mean，同时在Envlop_sd_x(m)中搜索大于std_mean的所有标准差，依次存储到新的标准差序列Envlop_sd_y(q)中，并计算新标准差序列Envlop_sd_y(q)的算术平均数Envlop_sdm；

其中：q＝1、2、3……N；

44)判断预存储的阈值Envlop_sdth是否小于新标准差序列Envlop_sd_y(q)的算术平均数Envlop_sdm，如是，判定为本次存储的数据具有完整的波峰特征，则进入步骤5)，如否，判定为本次存储的数据不具有波峰特征，则返回步骤1)；

5)分类确定步骤33)中降干扰处理后的包络数据Envlop_F(i)的判断阈值Ath，其中判断阈值Ath的计算步骤如下：

51)取步骤33)中包络数据Envlop_F(i)，计算其算术平均数Envlop_F_mean；

52)对包络数据Envlop_F(i)进行L等分，对每一等分按顺序依次求平均值，将各平均值依次存入新的包络序列Envlop_F_x(l)，搜索获得新包络序列Envlop_F_x(l)中所有大于Envlop_F_mean的值，将搜索到的各值依次存入序列Envlop_F_y(p)，之后，计算序列Envlop_F_y(p)的算术平均数，序列Envlop_F_y(p)的算术平均数即为峰值搜索阈值Ath；

其中：l＝1、2、3……L；

P＝1、2、3……P，P<L；

P表示为Envlop_F_y(p)序列中值的个数；

6)计算包络数据的峰值位置并标记，具体计算步骤如下：

61)分别从步骤3)、步骤52)获取包络序列Envlop_F(i)、判断阈值Ath，搜索包络序列Envlop_F(i)中大于判断阈值Ath的包络区间，记为Envlop_I(r)；

其中：r＝1、2、3……L；

62)对包络区间Envlop_I(r)作一阶后向差分运算，获得新的序列从中查找最接近0的值并记录下此值的位置，此位置即为峰值位置；所述的一阶后向差分运算的公式如下：

$\&dtri;Envlop\_I\left(r\right)=Envlop\_I\left(r\right)-Envlop\_I(r-1)$

其中：表示前向差分算子；

7)记录相邻两个峰值的位置，所述两个峰值之间的包络区间即为一个完整的心动周期。

本发明提供的超声频谱完整心动周期标记的方法，与传统的超声频谱参数实时自动计算方法相比，无需外部预测或者估计过多的经验参数，在信噪比较低、混叠、多峰等复杂情况下也能准确搜索到包络峰值，更有利于最大值以及完整心动周期的估计，从而精确估算血流的各项参数指标；即使对复杂频谱，也可准确计算包络线每个心动周期内的峰值，从而基于此可准确计算其他相关血流医学参数指标。

附图说明

图1为现有技术中多普勒信号处理流程图；

图2为本发明一种超声频谱完整心动周期标记的方法的实施流程图；

图3为本发明方法分类判断波峰状态步骤的实施流程图；

图4为本发明方法分类搜索波峰阈值步骤的实施流程图；

图5为一个心动周期标识示意图；

图6为本发明超声频谱完整心动周期标记的方法一种具体实施例的最终效果图。

具体实施方式

本发明提供了一种超声频谱完整心动周期标记的方法，具体包含以下步骤：

1)预存储经过解调、滤波、DFT处理的t-Ts到t时刻时间内的声功率谱数据Spectral(i，j)以及对应的正向包络数据EnvlopP(i)、负向包络数据EnvlopN(i)；

其中：i表示时间方向谱线序号,t-Ts≤i≤t；

j表示频率方向谱线的点序号；

之后，进入步骤2)；

2)取经步骤1)存储的t-Ts时刻到t时刻的谱数据Spectral(i，j)以及对应的正向包络数据EnvlopP(i)、负向包络数据EnvlopN(i)，进行计算以确定当前声功率谱的状态分布，具体的计算确定步骤如下：

21)获取缓存区内所有的从t-Ts到t时刻时间内的声功率谱数据Spectral(i，j)，计算其平均值，再计算所有从t-Ts到t时刻时间内的声功率谱数据Spectral(i，j)的正向累加和、负向累加和，公式如下：

从t-Ts到t时刻时间内的声功率谱数据Spectral(i，j)的平均值为MeanSpectral(j)，具体计算公式如下：

$MeanSpectral\left(j\right)=\frac{\underset{i=t-T_s}{\overset{i=t}{\&Sigma;}}Spectral(i,j)}{n};$

从t-Ts到t时刻时间内声功率谱数据Spectral(i，j)的正向累加和为SumSpectralP，具体计算公式如下：

$SumSpectralP=\underset{j=basline}{\overset{j=n-1}{\&Sigma;}}MeanSpectral\left(j\right);$

从t-Ts到t时刻时间内声功率谱数据Spectral(i，j)的负向累加和为SumSpectralN，具体计算公式如下：

$SumSpectralN=\underset{j=0}{\overset{j=basline-1}{\&Sigma;}}MeanSpectral\left(j\right);$

其中:baseline表示量化后的基线位置，0≤baseline＜n；

n表示单根谱线采样点最大值；

Spectral(i，j)表示i时刻离散谱线数据；

22)计算缓存区内t-Ts到t时刻时间之间的正向包络数据EnvlopP(i)、负向包络数据EnvlopN(i)的正向累加和、负向累加和，公式如下：

所有谱线包络的正向累加和为：

$SumEnvlopP=\underset{i=t-Ts}{\overset{i=t}{\&Sigma;}}\left(Envlop\right(i)-baseline+1);$

所有谱线包络的负向累加和为：

$SumEnvlopN=\underset{i=t-Ts}{\overset{i=t}{\&Sigma;}}(baseline-Envlop(i\left)\right)$

其中：Envlop(i)表示从t-Ts时刻到t时刻之间谱线的包络值；

SumEnvlopP表示谱线包络正向累加和；

SumEnvlopN表示谱线包络负向累加和；

23)取步骤21)中频谱的正向累加和SumSpectralP、负向累加和SumSpectralN，取步骤22)中包络的正向累加和SumEnvlopP、负向累加和SumEnvlopN，分别比较判断频谱的正向累加和SumSpectralP与频谱的负向累加和SumSpectralN、包络的正向累加和SumEnvlopP与包络的负向累加和SumEnvlopN，以区分缓存区内t-Ts到t时刻之间正向包络EnvlopP(i)、负向包络EnvlopN(i)的分布状态，分布状态公式如下：

如果SumEnvlopP≥SumEnvlopN&SumSpectralP≥SumSpectralN，那么

Envlop(i)＝EnvlopP(i)，进入步骤3)；

如果SumEnvlopP＜SumEnvlopN&SumSpectralP＜SumSpectralN，那么

Envlop(i)＝n-EnvlopN(i)，进入步骤3)；

否则返回步骤1)，重新开始；

其中：&表示逻辑“与”运算；

n表示单根谱线采样点最大值，与步骤21)中一致；

Envlop(i)为确定状态的声功率谱包络，包括正向包络EnvlopP(i)或者负向包络EnvlopN(i)；

3)对步骤2)中已确定状态的声功率谱包络Envlop(i)进行降干扰处理，降干扰处理后的包络记为Envlop_F(i)，降干扰处理步骤如下：

31)首先对Envlop(i)去除毛刺干扰，去除毛刺干扰计算公式如下：

Envlop_M_F(i)＝Get2[Reshape[Envlop(i-1),Envlop(i),Envlop(i+1)]]；

Envlop_M_F(i)为去除毛刺干扰后的包络数据；

32)其次对Envlop_M_F(i)抑制高频干扰，公式如下：

Envlop_H_F(i)＝[Envlop_M_F(i-1)+Envlop_M_F(i)+Envlop_M_F(i+1)]/3；

Envlop_H_F(i)为抑制高频干扰后的包络数据；

33)最后对Envlop_H_F(i)凸出峰值，公式如下：

Envlop_F(i)＝Max[Envlop_H_F(i-1),Envlop_H_F(i),Envlop_H_F(i+1)]；

Envlop_F(i)为凸出峰值处理后的包络数据；

其中：Reshape[]函数表示对[]中三个数进行从小到大排序；

Get2[]函数表示取[]内序列的第2个值；

Max[]表示取[]中三个值的最大值；

4)分类确定步骤3)中降干扰后的包络数据Envlop_F(i)是否具有完整的波峰特征，其中分类确定步骤如下；

41)预设定包络判断阈值Envlop_sdth并存储；

42)将步骤3)中降干扰后的包络数据Envlop_F(i)进行N等分，对每一等分计算其标准差，标准差公式如下：

$Envlop\_sd\_x\left(m\right)=\sqrt{\frac{1}{M}\underset{k=1}{\overset{k=M}{\&Sigma;}}{\left(Envlop\right(k)-u)}^2}$

其中：M为每等分的采样点数；

u为每一等分的平均值；

m＝1、2、3……N；

k＝1、2、3……M；

Envlop(k)表示对包络数据Envlop_F(i)进行N等分后的第k个等分序列表达式；

43)计算步骤42)中Envlop_sd_x(m)的算术平均数std_mean，同时在Envlop_sd_x(m)中搜索大于std_mean的所有标准差，依次存储到新的标准差序列Envlop_sd_y(q)中，并计算新标准差序列Envlop_sd_y(q)的算术平均数Envlop_sdm；

其中：q＝1、2、3……N；

44)判断预存储的阈值Envlop_sdth是否小于新标准差序列Envlop_sd_y(q)的算术平均数Envlop_sdm，如是，判定为本次存储的数据具有完整的波峰特征，则进入步骤5)，如否，判定为本次存储的数据不具有波峰特征，则返回步骤1)；

5)分类确定步骤33)中降干扰处理后的包络数据Envlop_F(i)的判断阈值Ath，其中判断阈值Ath的计算步骤如下：

51)取步骤33)中包络数据Envlop_F(i)，计算其算术平均数Envlop_F_mean；

52)对包络数据Envlop_F(i)进行L等分，对每一等分按顺序依次求平均值，将各平均值依次存入新的包络序列Envlop_F_x(l)，搜索获得新包络序列Envlop_F_x(l)中所有大于Envlop_F_mean的值，将搜索到的各值依次存入序列Envlop_F_y(p)，之后，计算序列Envlop_F_y(p)的算术平均数，序列Envlop_F_y(p)的算术平均数即为峰值搜索阈值Ath；

其中：l＝1、2、3……L；

P＝1、2、3……P，P<L；

P表示为Envlop_F_y(p)序列中值的个数；

6)计算包络数据的峰值位置并标记，具体计算步骤如下：

61)分别从步骤3)、步骤52)获取包络序列Envlop_F(i)、判断阈值Ath，搜索包络序列Envlop_F(i)中大于判断阈值Ath的包络区间，记为Envlop_I(r)；

其中：r＝1、2、3……L；

62)对包络区间Envlop_I(r)作一阶后向差分运算，获得新的序列从中查找最接近0的值并记录下此值的位置，此位置即为峰值位置；所述的一阶后向差分运算的公式如下：

$\&dtri;Envlop\_I\left(r\right)=Envlop\_I\left(r\right)-Envlop\_I(r-1)$

其中：表示前向差分算子；

7)记录相邻两个峰值的位置，所述两个峰值之间的包络区间即为一个完整的心动周期。

图2显示了本发明的一种实施流程。

下面讲述了本发明方法的一个具体实施例，可以有效提取复杂频谱的医学参数：

1.设定长度为Ts的硬件缓存区存储声功率谱数据以及对应的包络数据，在计算PW参数时无论实时还是电影模式均从硬件缓存区中读取声功率谱数据以及对应的包络数据，当声功率谱数据刷新完预设定的时间范围Ts，则需重新存储新的声功率谱数据以及对应的包络数据到Ts秒的缓存区；

2.对Ts秒缓存区内所有声功率谱数据进行Ts维度上的平均后获得一根声功率平均谱线记为MeanSpectral，以基线为中心分别计算MeanSpectral正向累加和和负向累加和，分别记为MeanSpectralP、MeanSpectralN；计算Ts秒内所有谱线包络的正向累加和和负向累加和，分别记为SumEnvlopP、SumEnvlopN；如果SumEnvlopP大于或等于SumEnvlopN，同时MeanSpectralP大于或等于MeanSpectralN，取正向包络；如果SumEnvlopP小于SumEnvlopN并且同时MeanSpectralP小于MeanSpectralN，取负向包络；如果有其他情况，则需要重新缓存Ts秒的数据并重新计算判断阈值以确定包络使用情况，如图3所示；

3.对缓存区内声功率谱数据对应的包络数据进行m点中值滤波，m可取3、5、7……点，可去除孤立点的干扰；对中值滤波后的包络数据进行n点的平滑处理，n可以取3、5、7…..点，平滑包络数据上的凹凸缺陷；对平滑后的包络数据进行极值滤波，放大波峰、波谷的极值效应；

4.预设置包络判断阈值，记为Envlop_sdth，把缓存区内存储的包络数据N等分，对每一等分计算其标准差记为Envlop_sd_x(m)；对Envlop_sd_x(m)序列求平均，记为std_mean，把Envlop_sd_x(m)中所有大于std_mean的数据取出并求所取出的大于std_mean的数据的平均值，记为Envlop_sdm；判断如果Envlop_sdm大于等于Envlop_sdth，则判定为当前存储的数据包含心动周期的特征，继续后续运算，如果判断Envlop_sdm小于Envlop_sdth，则需重新刷新Ts秒数据，返回步骤3)重新计算；

5.对Envlop_F(i)进行L等分，每一等分求平均值后获得新的包络序列Envlop_F_x(l)，计算Envlop_F(i)序列的算术平均数记为Envlop_F_mean，搜索新包络序列Envlop_F_x(l)中所有大于Envlop_F_mean值的平均包络值，记为Envlop_F_y(p)，计算序列Envlop_F_y(p)的算术平均数，此即为峰值搜索阈值Ath，如图4所示；

6.搜索大于判断阈值Ath的包络区域，并记录其值和标签值，记为(p，Envlop_new(p))；确定包络区间函数Envlop_new(p)的起点S(s，Envlop_new(s))和终点E(e，Envlop_new(e))；确定包络区间函数Envlop_new(p)的差分函数Envlop_diff(p)；由E、S、Envlop_new、Envlop_diff确定当前包络区间的最大值的标签Max；

7.在实际病人血流检测中，需要去除一些常见的异常情况，如果最大值Max为E，则确定为单调上升区间，去除此点；如果最大值Max为S，则确定为单调下降区间，去除此点；

8.找出完整心动周期峰值后，需按照峰值位置查找波谷位置，临床应用的标记线为标记两个波谷位置，如图5所示；

9.在本发明方法实现的应用程序界面上设置自动计算打开和关闭开关，用户在实时状态下可以自行选择打开和关闭，大大方便了临床的应用，最终效果图如图6所示。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应所述以权利要求的保护范围为准。