CN106943158A - 一种频谱参数实时计算方法 - Google Patents

Abstract

本发明公开了一种频谱参数实时计算方法，包括循环进行的步骤：取出一段最大流速曲线数据；对此数据进行预处理后；根据预处理后的数据求解两个滑动平均值MAbeat与MApeak，进一步求出阈值THR1；若有连续N个最大流速点对应的THR1大于MAbeat值，则认为此区域存在最高流速点，是感兴趣区域；统计每个感兴趣区域的长度，在每个感兴趣区域的长度范围内搜索，得到取出的这段数据中所有的最高流速点。本发明搜索阈值设置为最大流速曲线的滑动平均值，是根据最大流速曲线实时计算得到的，是完全自适应的。

Description

一种频谱参数实时计算方法

技术领域

本发明涉及一种计算方法，具体是一种频谱参数实时计算方法。

背景技术

超声多普勒技术被广泛用于人体血流的无损检测和测量，连续波多普勒(CW)和脉冲波多普勒(PW)技术均属于频谱多普勒技术，即对多普勒血流信号进行频谱分析，获得其频谱分布，从而可以估计出血管内血流速度的分布情况。

如图7所示为人体颈动脉多普勒信号声谱图。由该声谱图上提取最大频率(也称最大流速)，得到最大流速曲线，即声谱图的包络(上图曲线描绘了一个周期的最大流速曲线)，从中估计一些重要的临床参数，如心率、最高流速、平均流速等。

图7所示，曲线包络称为最大流速曲线，一个心动周期的最大流速曲线的峰值称为最高流速点(图示PS点)。

相近的技术方案CN 101301212B；对上述已有方案描述如下：

一种实时估计多普勒参数的方法及装置，用于超声诊断系统中，对运动组织或者血流声谱图的流速曲线自动进行多普勒参数实时计算处理。所述方法包括循环进行的步骤：用长度预先设置的数据缓冲区依次从所述流速曲线上取出一段数据，进行当前准心动周期的估计；确定最高流速的当前搜索阈值；根据所述阈值和准心动周期，搜索当前最高流速。

其中，阈值的设置，是初始设置或重新设置为当前预定时间长度T的一段流速曲线所对应的流速最大值和平均值的平均值；其中T的选取确保该曲线至少包含一个心动周期，T选取为2s；波峰搜索过程：先将每个写入缓冲区的最大流速值与搜索阈值比较，仅当连续有N2个最大流速值都大于该阈值时，以该N2个点中的第一个点对应的时刻为起点S，确定一预定长度范围内的时间段为所述波峰搜索期，本方案中N2对应时间长度为0.03秒，预定长度单位的时间段在本实施例中取为准心动周期的三分之一；进而，在此时间范围内搜索出最大曲线峰值，作为所述最高流速的估计值。再根据最高流速所在的位置实时估计包括当前心动周期、最低流速或心率在内的其他参数。

上述为确保搜索阈值的实时性和自适应性，使用刚刚计算出的当前心动周期内的最高流速和平均流速更新搜索阈值。在某些情况下，受外部因素印象，最大流速曲线的幅度可能发生变化。例如，患者身体移动可能导致多普勒声谱图及最大流速曲线全部变小，则以之前更新的阈值搜索最高流速，会出现错误，所以方案中还加了对搜索阈值有效与否的判断机制。

已有的方案中的搜索阈值以及最高流速的搜索时间范围都是根据前一个周期的参数计算得到的，在出现患者移动等特殊情况时，会出现错误，自适应差。

对搜索阈值有效性的判断与更新机制使运算复杂。

发明内容

本发明的目的在于提供一种频谱参数实时计算方法，以解决上述背景技术中提出的问题。

为实现上述目的，本发明提供如下技术方案：

一种频谱参数实时计算方法，包括循环进行的步骤：取出一段最大流速曲线数据；对此数据进行预处理后；根据预处理后的数据求解两个滑动平均值MAbeat与MApeak，进一步求出阈值THR1；若有连续N个最大流速点对应的THR1大于MAbeat值，则认为此区域存在最高流速点，是感兴趣区域；统计每个感兴趣区域的长度，在每个感兴趣区域的长度范围内搜索，得到取出的这段数据中所有的最高流速点。

作为本发明再进一步的方案：所述取出的这段数据包含多个心动周期的最大流速曲线，能够搜索出多个最高流速点。

与现有技术相比，本发明的有益效果是：本发明搜索阈值设置为最大流速曲线的滑动平均值，是根据最大流速曲线实时计算得到的，是完全自适应的；判断的方法不是用最大流速曲线值直接与阈值比较，而是用一段范围的不同时间长度内的平均值进行比较，当由于患者移动等其他情况影响，导致最大流速曲线的幅度整体变大或变小时，求出的滑动平均值也随着变化，仍能准确判断存在峰值点的感兴趣区域，避免了由于干扰带来的阈值无效问题，也不需要进行阈值有效与否的判断操作。算法具有更好的自适应性，鲁棒性更高，更简便，并且在每一段存在最高流速点(PS点)的感兴趣区域内搜索时，搜索的长度也是自适应的，具有更高效率。

附图说明

图1为频谱参数实时计算方法中PW实时描迹算法流程图。

图2为频谱参数实时计算方法中实时显示最近2个周期的最大流速曲线图。

图3为频谱参数实时计算方法中流速曲线及其二阶导数图。

图4为频谱参数实时计算方法中实施例1的结构框图。

图5为频谱参数实时计算方法中实施例2的结构框图。

图6为频谱参数实时计算方法中第一个周期的曲线图。

图7为现有技术中实时显示最近2个周期的最大流速曲线图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

请参阅图1～7，本发明实施例中，一种频谱参数实时计算方法，包括循环进行的步骤：取出一段最大流速曲线数据；对此数据进行预处理后；根据预处理后的数据求解两个滑动平均值MAbeat与MApeak，进一步求出阈值THR1；如果有连续N个最大流速点对应的THR1大于MAbeat值，则认为此区域存在最高流速点(PS点)，是感兴趣区域；统计每个感兴趣区域的长度，在每个感兴趣区域的长度范围内搜索，可以得到取出的这段数据中所有的最高流速点(取出的数据可能包含多个心动周期的最大流速曲线，能够搜索出多个最高流速点)。

本发明中，设置了两个阈值参数，THR1与THR2；

THR1＝MAbeat+beta*mean；

THR2＝W1；

对不同的最大流速点，MAbeat的取值不同，MAbeat求解方法为：以当前的最大流速点的时间位置为中间位置，取一段长度为W2的时间段内的最大流速点的平均值作为对应当前最大流速点的MAbeat；THR1是与当前点直接相关的，自适应强；beta是一个加权系数，取值在0～0.1之间；mean代表取出的所有数据的平均值。

MApeak的求解方法：以当前的最大流速点的时间位置为中间位置，取一段长度为W1的时间段内的最大流速点的平均值作为对应于当前最大流速点的MApeak值；

如果一个最大流速点对应的MApeak的数值大于THR1的数值，则认为它在感兴趣区域内；

其中MAbeat与MApeak相比，对应的时间窗长度更大；

搜索阈值设置为最大流速曲线的滑动平均值，是根据最大流速曲线实时计算得到的，是完全自适应的；判断的方法不是用最大流速曲线值直接与阈值比较，而是用一段范围的不同时间长度内的平均值进行比较，当由于患者移动等其他情况影响，导致最大流速曲线的幅度整体变大或变小时，求出的滑动平均值也随着变化，仍能准确判断存在峰值点的感兴趣区域，避免了由于干扰带来的阈值无效问题，也不需要进行阈值有效与否的判断操作。算法具有更好的自适应性，鲁棒性更高，更简便。并且在每一段存在最高流速点(PS点)的感兴趣区域内搜索时，搜索的长度也是自适应的，具有更高效率。

实施例1：

本发明方法用于超声诊断系统中，对运动组织或者血流声谱图的流速曲线自动进行多普勒参数实时计算处理，包括循环进行的步骤：

1.根据当前参数设置的需要实时描迹的最大流速曲线的周期数，设置缓冲区长度，从最大流速曲线取出一段数据；

图1所示，设置的实时显示最近1个周期的最大流速曲线；此时缓冲区长度对应为(1+N_period_rich)个周期；设N_period_rich取值1；即取出2个周期的包络数据；

图2所示，设置的实时显示最近2个周期的最大流速曲线；此时缓冲区长度对应为(2+N_period_rich)个周期，即取出3个周期的包络数据；

其中N_period_rich取1，2或者3。

最终存储在缓冲区中的最大流速曲线的点数设为N；

2.对取出的数据进行带通滤波处理，通带范围0.5～8Hz

下面步骤分为两部分，一部分用来求最高流速点，即最大流速曲线的峰值点(PS点)；一部分用来求ED点；

3.求PS点步骤如下：

(1)对经过带通滤波处理的最大流速曲线进行预处理，突出峰值点；

(2)求两个滑动平均值MAbeat，MApeak的数值，和两个阈值THR1，THR2；假设上一步处理后的曲线为y(n)，1<＝n<＝N，N代表取出数据的总长度；对应MAbeat的窗函数长度为W2,

THR1(n)＝MAbeat(n)+beta*Mean

其中，W2是预先设置的滑动平均窗的时间长度，取值范围在545ms～694ms；beta是预先设置的加权值，取值在0～0.1范围内；

其中，W1是预先设置的滑动平均窗长度，取值范围在54ms～120ms内；

THR2＝W1

(3)根据两个滑动平均值以及两个阈值判断感兴趣区域(存在PS点的区域)，并且记录每个感兴趣区域的时间长度Tlength(N_blocks)对取出的N个数据点，依次判断THR1(n)与MApeak(n)的大小；如果有连续N1(其中，N1>N_W1)个点对应的THR1的数值大于MAbeat的数值，则认为这个区域内存在PS点；其中N_W1是对应一段长度为W1的时间段内的最大流速数据点数。同时统计每个感兴趣区域内连续满足MApeak>THR1条件的点的数目，即每个感兴趣区域的时间长度Tlength；N_blocks代表在存储在缓冲区中的最大流速曲线中存在PS点的区域的个数，即PS点的个数。

(4)对每个感兴趣区域内搜索最大值，即该区域对应的心动周期内的最高流速点(PS点)，每个感兴趣区域的搜索时间范围即是上一步求出的对应该区域的Tlength。

至此，找到了取出的这段数据中所有的最高流速点(PS点)。

4.求ED点步骤如下：

(1)对步骤2中经过带通滤波处理的曲线数据首先求出其二阶导数曲线，二阶导数曲线幅度最大的点，对应最大流速曲线上的ED点；如图3所示，实线代表最大流速曲线，虚线代表最大流速曲线的二阶导数曲线；“o”代表最高流速点(PS点)；“*”代表ED点，可以看到ED点对应二阶导数曲线上的最大值点；

(2)对二阶导数曲线进行预处理，突出二阶导数曲线的峰值；

(3)求两个滑动平均值MAbeat，MApeak和两个阈值THR1、THR2，设二阶导数曲线为y1(n)，1<＝n<＝N，N代表缓冲区存储的数据点数

THR1(i)＝MAbeat(i)+beta*mean

其中，W2是预先设置的滑动平均窗的长度，这里W2取值范围在100ms～200ms；beta是预先设置的加权值，取值范围在0～0.1；

其中，W1是预先设置的滑动平均窗的时间长度，这里W1取值范围在1000ms～1.25s；在求PS点与求ED点的过程中用到的W1,W2，beta的取值范围不同。

THR2＝W1；

下面采用与求PS点过程相同的处理，根据两个滑动平均值以及两个阈值进行判断，找到取出的N个最大流速点数据中存在的所有的ED点。

进一步根据PS点和ED点进行其他参数的计算。

实施例2：

采用本发明的方法的一般用于超声频谱多普勒系统中；

一个典型的超声频谱多普勒系统的装置结构如图4所示，是一个超声频谱多普勒系统框图，分为发射超声波部分和接收超声波部分；本发明提出的方法用于对接收的超声波处理过程中的“参数实时计算”模块。

实施例3：

如图5所示，包括：

控制单元：用于用户设置需要实时描迹的多普勒频谱最大流速曲线的周期数目；这只是否需要描迹平均流速曲线；

数据存储单元：根据控制单元计算并设置缓冲区大小；实时读入并存储一定长度的最大流速曲线数据；

数据预处理单元：用于对取出的最大流速曲线进行预处理操作，消除干扰；

最高流速检测单元：用于实时计算两个滑动平均值和两个阈值，并且通过这几个值判断，计算峰值搜索范围大小；搜索得到所有的最高流速点(PS点)；

ED点检测单元：用于对数据进行进一步处理；得到最大流速曲线的二阶导数曲线，根据二阶导数曲线计算两个滑动平均值和两个阈值，以及计算每个区域的峰值搜索范围大小；搜索得到二阶导数曲线的峰值，对应的就是最大流速曲线上ED点的位置；

参数计算单元：应用搜索到的最高流速点(PS点)和ED点的坐标，计算相关参数；

参数输出单元：将计算得到的参数输出到显示屏进行实时显示。

本发明中用来搜索最高流速点(PS点)的阈值是根据当前最大流速曲线的数值实时更新计算的，是完全自适应的，避免了特殊情况下频谱幅度整体移动导致的阈值无效，无法找到PS点的情况；也免去了对搜索阈值有效性判断的单元；阈值随着最大流速点自适应的变化，提高了鲁棒性，和算法的稳定性，并且简化了算法。

在每个区域内搜索最高流速点(PS点)的搜索时间长度与该区域内满足MApeak>THR1的点数相关，是自适应变化的，使搜索更有效，也减小了搜索范围。

如图6所示，第一个周期的曲线，由于特殊情况，频谱曲线幅度值整体偏小，但是本发明提出的方法仍然准确估搜索到PS点和ED点的位置.

在求PS点和求ED点的预处理过程，可以包含，将待处理曲线x(n)小于零的部分置零，再对曲线求它的n次方(n>1)，或者求它的对数等任何用来突出峰值点的操作；例如，设预处理之后的曲线为

x_after(n)，1<＝n<＝N

或者

或者，

x_after(n)＝10^x(n)0≤n≤N。

求最大流速曲线的二阶导数的过程，可以首先对一阶导数进行带通滤波，以消除干扰，同样在求出二阶导数之后也可进行一定的带通滤波处理

对于本领域技术人员而言，显然本发明不限于上述示范性实施例的细节，而且在不背离本发明的精神或基本特征的情况下，能够以其他的具体形式实现本发明。因此，无论从哪一点来看，均应将实施例看作是示范性的，而且是非限制性的，本发明的范围由所附权利要求而不是上述说明限定，因此旨在将落在权利要求的等同要件的含义和范围内的所有变化囊括在本发明内。不应将权利要求中的任何附图标记视为限制所涉及的权利要求。

此外，应当理解，虽然本说明书按照实施方式加以描述，但并非每个实施方式仅包含一个独立的技术方案，说明书的这种叙述方式仅仅是为清楚起见，本领域技术人员应当将说明书作为一个整体，各实施例中的技术方案也可以经适当组合，形成本领域技术人员可以理解的其他实施方式。

Claims (2)

1.一种频谱参数实时计算方法，其特征在于，包括循环进行的步骤：取出一段最大流速曲线数据；对此数据进行预处理后；根据预处理后的数据求解两个滑动平均值MAbeat与MApeak，进一步求出阈值THR1；若有连续N个最大流速点对应的THR1大于MAbeat值，则认为此区域存在最高流速点，是感兴趣区域；统计每个感兴趣区域的长度，在每个感兴趣区域的长度范围内搜索，得到取出的这段数据中所有的最高流速点。

2.根据权利要求1所述的频谱参数实时计算方法，其特征在于，所述取出的这段数据包含多个心动周期的最大流速曲线，能够搜索出多个最高流速点。