CN1875887A - 声谱图包络曲线提取方法 - Google Patents

Abstract

一种声谱图包络曲线提取方法，用于频谱多普勒技术对血流速度的测量过程，包括步骤：射频超声回波信号经解调、滤波和模数转换后获得流体多普勒信号；对所述多普勒信号进行频谱分析获得对应功率谱P(f)；对一预定时刻的功率谱P(f)估计正向最大频率f_max+和反向最大频率f_max－；根据该两最大频率来确定噪声频率范围从而估计噪声平均功率E；利用所述噪声平均功率E来校正所述正向最大频率和反向最大频率。采用本发明方法，可以减小信噪比和带宽对包络曲线的影响，从而有利于血流参数的准确计算。

Description

声谱图包络曲线提取方法

技术领域  本发明涉及超声技术，特别涉及利用频谱多普勒技术测量流体流速的超声技术，尤其是涉及测量血流速度的声谱图包络曲线。

背景技术  超声技术被广泛用于医学成像和测量，其中利用多普勒技术来测量血流速度的研究或实践自1957年以来已大量地在物理、工程和临床展开。在连续波多普勒系统测量血流速度的过程中，超声诊断系统发射超声波到人体内某采样容积，并接收该采样容积内血细胞造成的散射回波。由于血细胞运动产生多普勒效应(散射回波与发射波具有一定的频率偏移，且该偏移与血细胞运动速度成正比)，所以通过测量频偏可以估计血管内血液流动状况。在脉冲波多普勒系统测量血流速度的过程中，系统是按一定时间间隔重复发射脉冲信号，而在连续两次发射之间的某一时刻接收散射回波，再通过测量不同散射回波与某一参考信号相位差随时间的变化率来实现流速测量。虽然使用脉冲波多普勒技术测量流速在原理上不依赖于多普勒效应，但最后测得的信号频偏与血细胞运动速度之间关系则与多普勒效应一致，因此，工程和临床上对脉冲波和连续波多普勒信号的处理方法一般相同。

射频超声散射回波信号经过正交解调、滤波后，多普勒信号频谱被从几兆赫兹搬迁到以0频率为中心，带宽约几千赫兹的音频范围内(在该范围内如此形成的信号通常称为音频多普勒信号，为了下面叙述方便，简称之为多普勒信号)。由于血液的黏性，血流速度在人体血管内具有一定分布，一般在近血管壁处接近0，而近血管中心处较大，因此多普勒信号频谱具有一定频谱宽度，其中最大频率正比于血管内最大血流速度。此外，由于血管内血流速度随心脏收缩和舒张而不断变化，检测到的多普勒信号频谱也随之不断变化。截取某一时刻的多普勒信号进行频谱分析，我们可以估计出该时刻血管内血流状况；若按一定时间间隔多次截取多普勒信号进行频谱计算，用灰度调制不同频率成分的功率大小，并按时间顺序显示，即可获得多普勒信号声谱图。

为了利用多普勒信号声谱图对血管疾病进行诊断，需要从声谱图上提取出一些参数，比如平均流速，最大流速，最小流速，S/D(收缩期、舒张末期最大流速比)，RI(阻尼指数)，PI(脉动指数)等。这些参数都可以基于声谱图包络曲线或平均频率曲线计算得到。所述声谱图包络曲线可以通过连接不同时刻多普勒信号频谱的最大频率来得到，该曲线幅度与血管内血流最大流速成正比。所述最大频率是估算平均频率和计算其它多普勒参数的基础，因此准确地估计最大频率在临床应用中具有重要意义。

传统估计最大频率的方法是：多普勒系统操作者判断显示中的声谱图满足要求特性后，冻结该声谱图，手工描出最大频率曲线，再由软件计算平均频率或其它相关参数。它具有明显的缺点：重复性差，估计精度低和无法实时估计。由于动脉血流搏动及噪声的影响，最大频率的估计一直是频谱多普勒技术应用中的一个难点。随着数字计算机技术的高速发展，研究者提出了许多估计最大频率的估计方法，包括百分比法(PM)、过阈值法(TCM)、改进的过阈值法(MTCM)、混合法(HM)、几何法(GM)、改进的几何法(MGM)和自适应阈值法(ATM)等。这些方法一般被用于基于快速傅里叶变换(FFT)的谱估计系统，但也可用于其它谱估计系统。在Evans等的“Doppler Ultrasound：Physic，Instrumentation and Signal Processing”(见2000年John Wiley & Sons第二版，)里介绍和比较了所述各种方法。百分比法计算量最小，但受信噪比和带宽影响最大。过阈值一类的方法由于随机噪声(其频谱幅度变化很大)的存在，若阈值设定过小就很容易检测到较大幅度的噪声分量而导致最大频率估计偏大，但若阈值设置过大又会导致最大频率估计偏小(详见1988年参考文献“Comparison of four digitalmaximum frequency estimators for Doppler ultrasound”，Utrasound Med Biol第14卷、第5册、第355-363页)。虽然Routh等人在参考文献US 5,287,753和“Evaluation of anautomated real-time spectral analysis technique”(详见Ultrasound Med Biol第一卷、第一册、第61-73页)，以及Mo在US 5,935,074中又提出一些已被用于现有商用超声成像系统的自适应估计阈值方法，但实际操作中往往为了增加最大频率估计的鲁棒性而不得不设置较大阈值。其中自适应阈值法是根据信号的信噪比来设定阈值的。由于信号和噪声在频谱上的重叠，单纯利用阈值区分信号和噪声会给信噪比的估计带来很大误差和方差，进而影响最大频率的估计精度和鲁棒性。几何法不存在阈值选取问题，可以实现最大频率的自动提取；但因它是以选取的功率谱峰值位置为参考点，而利用FFT估计频谱具有很大方差，峰值位置往往具有很大随机性，所以将影响最大频率的估计精度。改进的几何法直接将积分功率谱曲线减去连接该曲线最大点和最小点的直线，由此获得的曲线最大值位置即估计为最大频率。与其他方法相比，该法具有更好鲁棒性(详见参考文献“The performance of three maximumfrequency envelope detection algorithms for Doppler signals”，J Vasc Invest第一卷、第126-134页)；但实验显示当信号频谱较宽时，该法也存在明显的最大频率估计偏小问题。

图1是典型的多普勒信号处理流程示意框图：由超声探头接收到的回波信号经波束合成、解调、滤波、模数转换后获得数字多普勒信号，然后利用FFT等对多普勒信号进行频谱分析来获得声谱图，再在该声谱图上用改进的几何法、自适应阈值法等技术来估计所述多普勒信号的最大频率，从而提取出声谱图包络曲线；最后根据所述包络曲线上的特征点来计算表征所述多普勒信号的各个多普勒参数。

上述现有技术的主要不足在于：最大频率估计偏低，尤其频谱较宽时该偏差更加明显；另外，过阈值相关方法还存在阈值参数设置对频率估计影响大的缺点；使得多普勒参数的计算结果精度较差。

发明内容本发明要解决的技术问题是针对上述现有技术的不足，而提出一种方法，能准确提取多普勒信号声谱图的包络曲线，减小了信号信噪比和带宽的影响，从而有利于血流参数的准确计算。

为解决上述技术问题，本发明的基本构思为：采用初估和校正相结合的方法，先利用现有技术初步估计功率谱正向最大频率和反向最大频率，然后根据初估的频率范围(利用尽量多的噪声点)来确定准确的噪声平均功率；利用该噪声平均功率来校正所述功率谱曲线，并获得校正后较精确的最大值和最小值。此外，还可以利用平均流速的方向或正反向血流功率谱的积分，来确定最后的最大频率输出，使系统具有更好的鲁棒性。

作为实现本发明构思的技术方案是，提供一种声谱图包络曲线提取方法，用于频谱多普勒技术对血流速度的测量过程，包括步骤：

A.射频超声回波信号经过解调、滤波和模数转换后获得流体的多普勒信号；

B.对所述多普勒信号进行频谱分析，获得该多普勒信号随时间变化的功率谱P(f)；

C.对一预定时刻的功率谱P(f)估计所述多普勒信号对应的正向最大频率f_max+和反向最大频率f_max-；

尤其是，还包括步骤：

D.根据所述正向最大频率和反向最大频率来确定噪声的频率范围，并在该范围内估计噪声的平均功率E；

E.利用所述噪声平均功率E来校正所述正向最大频率和反向最大频率。

上述方案中，所述步骤C估计所述正向最大频率和反向最大频率采用的是改进的几何法。

上述方案中，所述步骤D中噪声频率范围是步骤C中所述最大频率的函数，表示为：[f₀，F(f_max-]+[F(f_max-)，f₁]，其中F(x)表示为变量x的函数，[f₀，f₁]为功率谱P(f)的频率范围。

上述方案中，所述步骤E包括操作：

a.根据当前和以往时刻的噪声平均功率E来设置一个阈值Th(E)；

b.将所述多普勒信号功率谱P(f)减去所述阈值，P(f)-Th(E)，并对该差值函数积分，寻找该积分后函数φ”(f)的最大值和最小值，该最大值和最小值位置分别对应校正后的正向最大频率F_max+和反向最大频率F_max-。

上述方案中，所述步骤E后还包括步骤：

F.根据流速方向来确定最大频率输出。

上述方案中，所述根据流速方向来确定最大频率输出，由所述反向最大频率F_max-和正向最大频率F_max+之间的功率谱成分来估计的平均频率F_mean，及该平均频率对应的平均流速方向来确定；或者依据正向和反向血流功率谱的积分值大小来确定。

采用上述技术方案，可以使最大频率的估计精度受信号带宽和信噪比影响小，特别是同时存在正向和反向血流时，能更容易准确判断和输出正向最大频率或反向最大频率。

附图说明  图1是典型的多普勒信号处理流程示意图

          图2是本发明最大频率和平均频率估计方法的一种实现流程图

          图3是改进的几何法提取多普勒信号最大频率原理图

          图4是对图3中最大频率估计的校正处理示意图；

          图5是声谱图最大频率估计结果示例图

          图6是获得的平均频率曲线示例图

具体实施方式  下面，结合附图所示之最佳实施例进一步阐述本发明。

图1多普勒信号处理流程中包括的声谱图包络曲线提取方法，一般包括步骤：

A.射频超声回波信号经过解调、滤波和模数转换后获得流体的多普勒信号；

B.对所述多普勒信号进行频谱分析，获得该多普勒信号随时间变化的功率谱P(f)；

C.对一预定时刻的功率谱P(f)估计所述多普勒信号对应的正向最大频率f_max+和反向最大频率f_max-；

进而系统可以根据不同时刻功率谱P(f)的所述最大频率来生成显示用的多普勒信号声谱图的两条包络曲线，或者系统直接基于上述方法步骤产生的各不同时刻所述最大频率进行血流参数计算。本发明方法对其中最大频率估计进行了一些改进，在步骤C之后还包括步骤：

D.根据所述正向最大频率和反向最大频率来确定噪声的频率范围，并在该范围内估计噪声平均功率E；

E.利用所述噪声平均功率E来校正所述正向最大频率和反向最大频率。

图2是本发明估计方法的一种实现流程图。所述步骤C中的正向最大频率f_max+和反向最大频率f_max-可以利用现有技术进行初估，所述流程图包括但不限于先用改进的几何法来初步估计所述正向最大频率和反向最大频率。同样，所述步骤B中的对所述多普勒信号进行频谱分析的方法包括但不限于短时的快速傅里叶变换法。因相关方法属于现有技术，不再赘述。

改进的几何法处理信号过程如图3所示。以图3(a)一预定时刻的功率谱函数P(f)为例，将该函数从频率f₀开始积分获得如图3(b)实曲线所示的

$\mathrm{\φ}\left(f\right)={\∫}_{f_0}^{f}P\left(f^{\′}\right){\mathrm{df}}^{\′}$

连接该曲线φ(f)最小点和最大点的直线(如图中虚线所示)可以表示为

$C\left(f\right){=\∫}_{f_0}^f{\mathrm{kdf}}^{\′},$ 其中斜率k＝[φ(f₁)-P(f₀)]/(f₁-f₀)，f₁为正向频率范围

则图3(b)中两线相减可以得到得到图3(c)所示的

                   φ′(f)＝φ(f)-C(f)该曲线φ’(f)最大值对应的位置即为估计的正向最大频率值f_max+，最小值对应的位置为估计的反向最大频率值f_max-。图3中两道竖线所对应位置即为本实施例中所述初步估计获得的正向最大频率和反向最大频率，分别等于49和10Hz。

大部分多普勒信号成分集中在所述估计的反向最大频率和正向最大频率之间，而噪声信号分布于整个频谱范围内，所以大于正向最大频率和小于反向最大频率的频率成分则基本上以噪声信号为主。正如改进的几何法已被介绍的不足之处，可以明显从图3看到正向最大频率的估计偏小。因此，本发明步骤D利用所述信号分布特性来更准确地估计噪声的平均功率，以便上述初步估计的结果进行校正。

所述噪声频率范围可以描述为所述最大频率的函数，表示为：

      [f₀，F(f_max-)]+[F(f_max-)，f₁]

其中F(x)表示为变量x的函数，[f₀，f₁为功率谱P(f)的频率范围为简单起见，可以取F(x)＝x，则噪声的平均功率可以表示为：

$E=\frac{{\∫}_{f_0}^{f_{\max -}}P\left(f\right)\mathrm{df}+{\∫}_{f_{\max +}}^{f_1}P\left(f\right)\mathrm{dt}}{f_{\max -}-f_0+f_1-f_{\max +}}$

此后，为了利用所述噪声平均功率E来校正正向最大频率和反向最大频率，先要根据当前和以往时刻的噪声平均功率E来设置一个阈值Th(E)。由于每个时刻估计出的噪声平均功率都有一定的误差，利用以前时刻估计的噪声平均功率对当前估计值进行平滑，从而减小噪声平均功率E的误差。

因此Th(E)可以表示为当前和以往时刻功率谱噪声平均功率的函数Th(E)＝F(E₁，E_t-1，E_t-2，...)，其中E_t为当前t时刻的噪声平均功率，E_t-i为以往t-i时刻的噪声平均功率。简单起见该函数可取为Th(E)＝a*E_t，即阈值仅为当前时刻t估计出的噪声平均功率的函数，其中a＞1为一可调节的常数。

接着，将所述多普勒信号功率谱P(f)减去所述阈值，P(f)-Th(E)，并对该差值函数积分：

${\mathrm{\φ}}^{\′\′}\left(f\right){=\∫}_{f_0}^f[P\left(f^{\′}\right)-\mathrm{Th}\left(E\right)]{\mathrm{df}}^{\′}$

$={\∫}_{f_0}^{f}P\left(f^{\′}\right){\mathrm{df}}^{\′}-{\∫}_{f_0}^f\mathrm{Th}\left(E\right){\mathrm{df}}^{\′}$

$=\mathrm{\φ}\left(f\right)-C^{\′}\left(f\right)$

可见该操作与改进的几何法非常相似，所不同的只是第二项积分式表示的直线具有不同的斜率。上式中直线C(f)(见图4(b))的斜率只与噪声的平均功率有关，而在改进的几何法中直线C(f)的斜率与信号加噪声的平均功率有关。随着调节系数a的增大，直线C(f)越来越接近直线C(f)，估计出的最大频率也逐渐接近改进的几何法的结果。

再对积分后函数φ”(f)寻找最大值和最小值，该最大值和最小值位置分别对应校正后的正向最大频率F_max+和反向最大频率F_max-。结果如图4(c)所示。竖点直线所对应的校正后的正向最大频率F_max+和反向最大频率F_max-分别等于70和7Hz。从图4(a)可见它们都具有更好的估计精度。

所述正向最大频率和反向最大频率得到校正后，可以更准确地用来确定多普勒血流信号的频率范围，可以取所述两条包络曲线之间(即反向最大频率F_max-和正向最大频率F_max+之间的频率范围)的功率谱成分来估计信号的平均频率F_mean，例如利用频率的一阶矩计算如下：

$F_{\mathrm{mean}}=\frac{{\∫}_{F\max -}^{F\max +}f*P\left(f\right)\mathrm{df}}{{\∫}_{F\max -}^{F\max +}P\left(f\right)\mathrm{df}}$

由于不受血流速度范围外噪声的影响，上述方法估计出的平均频率具有更好精确度。本实施例中，该平均频率为32.6Hz，图4(a)中竖虚线示范了该频率的对应位置。

本发明方法还包括根据流速方向来确定最大频率输出的步骤，与直接利用功率最大值频率点位置确定最大频率输出相比，具有更好的鲁棒性。因为正交解调后的两路多普勒信号可以用于双向血流的检测，若本实施例，小于0Hz的频率范围对应反向血流，大于0Hz的对应正向血流。可以利用所述反向最大频率F_max-和正向最大频率F_max+之间的功率谱成分来估计平均频率F_mean；再根据该平均频率对应的平均流速方向来确定最后输出的最大频率F_max。本实施例中，平均频率32.6Hz＞0Hz，因此判断此时血管内平均流速为正向，从而最后输出的最大频率F_max为70Hz，即选择正向最大频率输出。

上述输出最大频率还可以依据正反向血流功率谱的积分来确定。先分别对流体的正向和反向流向对应的功率谱进行积分；再根据积分值大小来确定：若正向血流功率谱积分值较大，则选择正向最大频率F_max+，否则选择反向最大频率F_max-。

对不同时刻功率谱进行最大频率和平均频率估计后，可以连接这些频率对应点拟合形成相应包络曲线，并进一步显示于声谱图上。由于噪声影响，由该方法提取出的最大频率曲线可能存在一些误差，可以通过线性或者非线性的数字滤波进一步平滑曲线，减小误差(因属于现有技术，不再赘述)。

利用上述技术对脉冲多普勒系统检测到的颈动脉正交多普勒信号进行声谱图包络提取实验，可以得到如图5所示的结果，最大频率曲线如图中黑色曲线所示，平均频率曲线如图6中黑色曲线所示，均很好地跟踪了声谱图频谱包络。此外，所述平均频率曲线也非常的平滑，受噪声影响小。

Claims (10)

1.一种声谱图包络曲线提取方法，用于频谱多普勒技术对血流速度的测量过程，包括步骤：

A.射频超声回波信号经过解调、滤波和模数转换后获得流体的多普勒信号；

B.对所述多普勒信号进行频谱分析，获得该多普勒信号随时间变化的功率谱P(f)；

C.对一预定时刻的功率谱P(f)估计所述多普勒信号对应的正向最大频率f_max+和反向最大频率f_max-；

其特征在于，还包括步骤：

D.根据所述正向最大频率和反向最大频率来确定噪声的频率范围，并在该范围内估计噪声平均功率E；

E.利用所述噪声平均功率E来校正所述正向最大频率和反向最大频率。

2.根据权利要求1所述的声谱图包络曲线提取方法，其特征在于：

所述步骤B中对所述多普勒信号进行频谱分析采用的是短时的快速傅里叶变换。

3.根据权利要求1所述的声谱图包络曲线提取方法，其特征在于，

所述步骤C估计所述正向最大频率和反向最大频率采用的是改进的几何法。

4.根据权利要求1所述的声谱图包络曲线提取方法，其特征在于：

所述步骤D中噪声频率范围是步骤C中所述最大频率的函数，表示为：[f₀，F(f_max-)]+[F(f_max-)，f₁]，其中F(x)表示为变量x的函数，[f₀，f₁]为功率谱P(f)的频率范围。

5.根据权利要求1所述的声谱图包络曲线提取方法，其特征在于，所述步骤E包括操作：

a.根据当前和以往时刻的噪声平均功率E来设置一个阈值Th(E)；

b.将所述多普勒信号功率谱P(f)减去所述阈值，P(f)-Th(E)，并对该差值函数积分，寻找该积分后函数φ”(f)的最大值和最小值，该最大值和最小值位置分别对应校正后的正向最大频率F_max+和反向最大频率F_max-。

6.根据权利要求5所述的声谱图包络曲线提取方法，其特征在于：

所述步骤a中阈值Th(E)可以表示为当前和以往时刻功率谱噪声平均功率的函数：Th(E)＝F(E_t，E_t-1，E_t-2...)，其中E_t为当前t时刻的噪声平均功率，E_t-i为以往t-i时刻的噪声平均功率。

7.根据权利要求1所述的声谱图包络曲线提取方法，其特征在于，所述步骤E后还包括步骤：

F.根据流速方向来确定最大频率输出。

8.根据权利要求7所述的声谱图包络曲线提取方法，其特征在于：

所述根据流速方向来确定最大频率输出，由所述反向最大频率F_max-和正向最大频率F_max+之间的功率谱成分来估计的平均频率F_mean，及该平均频率对应的平均流速方向来确定；或者依据正向和反向血流功率谱的积分值大小来确定。

9.根据权利要求7或8所述的声谱图包络曲线提取方法，其特征在于：

所述最大频率F_max还通过线性或者非线性的平滑滤波来降低误差。

10.根据权利要求7任一权利要求所述声谱图包络曲线提取方法，其特征在于：

所述步骤F还包括：将拟合不同时刻最大频率F_max得到的包络曲线同时显示于声谱图上。

Images