CN1954235B - 用于超声波信号的局部频谱分析的改进方法和设备 - Google Patents

用于超声波信号的局部频谱分析的改进方法和设备 Download PDF

Info

Publication number
CN1954235B
CN1954235B CN200480035699XA CN200480035699A CN1954235B CN 1954235 B CN1954235 B CN 1954235B CN 200480035699X A CN200480035699X A CN 200480035699XA CN 200480035699 A CN200480035699 A CN 200480035699A CN 1954235 B CN1954235 B CN 1954235B
Authority
CN
China
Prior art keywords
coefficient
partial estimation
signal
frequency
matrix
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related
Application number
CN200480035699XA
Other languages
English (en)
Other versions
CN1954235A (zh
Inventor
伦纳尔多·玛索蒂
艾娜·比亚吉
西蒙纳·格兰奇
卢卡·布莱斯奇
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Esa Ott Co Ltd
Actis Active Sensors S R L
Esaote SpA
Actis Active Sensors Srl
Original Assignee
Esa Ott Co Ltd
Actis Active Sensors S R L
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Esa Ott Co Ltd, Actis Active Sensors S R L filed Critical Esa Ott Co Ltd
Publication of CN1954235A publication Critical patent/CN1954235A/zh
Application granted granted Critical
Publication of CN1954235B publication Critical patent/CN1954235B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/52Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S15/00
    • G01S7/52017Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S15/00 particularly adapted to short-range imaging
    • G01S7/52023Details of receivers
    • G01S7/52036Details of receivers using analysis of echo signal for target characterisation
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S15/00Systems using the reflection or reradiation of acoustic waves, e.g. sonar systems
    • G01S15/88Sonar systems specially adapted for specific applications
    • G01S15/89Sonar systems specially adapted for specific applications for mapping or imaging
    • G01S15/8906Short-range imaging systems; Acoustic microscope systems using pulse-echo techniques
    • G01S15/8977Short-range imaging systems; Acoustic microscope systems using pulse-echo techniques using special techniques for image reconstruction, e.g. FFT, geometrical transformations, spatial deconvolution, time deconvolution

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Acoustics & Sound (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • General Physics & Mathematics (AREA)
  • Ultra Sonic Daignosis Equipment (AREA)
  • Analysing Materials By The Use Of Radiation (AREA)
  • Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)
  • Investigating Or Analyzing Materials By The Use Of Ultrasonic Waves (AREA)

Abstract

本方法通过时间-频率变换之类的滤波而将射频信号分解到子频带中。从通过所述分解而采集的系数获取频谱系数矩阵,并且从该矩阵获取局部估计量,特别地,所述局部估计量是由内插多项式的系数构成的。局部估计量的统计分布在重叠在超声帧上的窗口中被评估。这些频谱系数分布直方图的构造提供了一个参数,该参数与局部估计量相结合,提供加权的局部估计量,所述加权局部估计量包含可以在对接受超声分析的器官中的特定结构进行识别的过程中使用的频谱信息。

Description

用于超声波信号的局部频谱分析的改进方法和设备
技术领域
本发明涉及一种用于处理超声波信号的方法,尤其涉及基于超声波方法来进行无侵害临床测试的领域。此外,本发明还涉及一种用于实施超声波信号处理方法的超声设备。
背景技术
超声扫描是一种用于对人类和动物身体的内部器官进行检查的方法,此外,它还是一种用于对包括非生物和无机物结构在内的其他类型的结构进行无侵害检查的方法。众所周知,该扫描是以产生一系列超声波信号(形成一个声束)为基础的,这些超声波信号将被引导到病人体内,并且将会激励一个超声波响应信号,该响应信号由存在于所观察的器官内部的散射体和/或反射体通过反射和/或散射声波而产生。发射和接收超声波信号的探测器将这个响应信号发送到接收和处理系统,此后,该信号将会变换成一个二维超声图像。
商用的超声扫描器提供的是与超声波信号幅度相关的图像,其中该超声波信号是由受到声束撞击的不同组织或界面反向传播的。
受到碰撞的不同组织或不同接触面是基于其回声特性而被检测到的,所述回声特性则由其声阻抗决定的,而所述声阻抗则又是密度与声速的乘积。
常规超声扫描器的目的是借助于时域中检测的幅度来重构和保持被检查组织的空间顺序,而这可以提供深度上的空间位置。针对信号所执行的调节处理旨在将信噪比增至最大,并在深度发生变化的时候使调查灵敏度保持恒定。
当在介质中传播声波并且声波遭遇到不均匀区域的时候,其中某些能量会朝着传感器反向传播,某些能量则会传播到组织内。
反向传播的声学信号包括各个反射体和/或散射体提供的时间重叠成分,由于分布的随机特性以及散射体和反射体的不同形状和机械特性,因此,结果得到的是一个形状极为复杂的信号。而传感器则会将所检查的介质反向传播的声学信号变换成被称为射频信号的电信号。
图1显示的是形成与生物组织有关的超声描迹的射频信号的部分。特别地,图1显示的是受到超声检查的器官的一部分返回的超声信号描迹的一部分。在这里,横坐标显示的是响应时间,也就是产生传感器得到的返回信号的散射体或反射体的深度位置。而在纵坐标上显示的则是信号幅度。对所检查的器官而言,其超声图像是从这种轨迹或描迹集合中重构的,其中所述轨迹或描迹集合将会形成信号帧。有时,这些描迹或轨迹也被称为扫描线。
第一种针对所述射频信号所执行的简单处理操作是借助于包络提取操作来提取信号幅度。在图2中,字母I表示的是对图1所示的射频信号执行包络操作所得到的结果。应该指出的是,包络信号跟随的是射频信号的正值波峰,由此提供了信号幅度的估计。
当前,在所有的超声成像技术中都使用了包络操作,该操作构成了单独描迹的A模式显示(其中A代表幅度)以及合成图像的B模式显示(其中B代表亮度)的基础,其中每一个单独的A模式描迹的幅度都是依照灰度或色标来进行编码的.
在现代的超声扫描器中,其中使用了超声波信号的包络,也就是说,检测到的信号是在A模式、B模式以及M模式描迹中描绘的。
对这种解决方案而言,包含在射频信号中的信息是有损耗的,也就是说,在接收传感器产生的形式中包含载波频率与所述接收传感器频带的中心频率近似相等的信号,其中该信号经过了具有相位阶跃(频率峰值)的调频(非对称频谱),并且在介质呈现出非线性效应(谐波频谱的产物)的情况下具有因为非线性失真所导致的形态。
图3显示的是图1的射频描迹在常规A模式显示中所具有的包络信号。应该指出的是,这其中不可挽回地损耗了关于原始信号频率和相位的信息。
特别地,如果通过包络提取操作而消除了射频:那么下列信息将会丢失:
●从超声和生物组织(局部声阻抗,该阻抗依取决于介质密度,所述密度又取决于入射超声波强度、对比剂中存在的气泡在谐波和次谐波上的次级辐射等等)之间的非线性交互中获取的信息;
●从线性交互中获取的信息,也就是因为选择性吸收、尤其是依赖于深度的选择性吸收以及轴向变化(scalloping)效应(在由于不连续性而被反向散射的信号中呈现的幅度频谱形状的缺口(indentation),其中所述不连续性包括空间中周期性分布的点类型的不连续性)所导致的信号频谱的修正。
因此,用于医疗诊断的常规超声扫描系统是以观察A模式描迹或B模式图像为基础的,但是该图像只显示了与所遭遇的不同组织结构相关联的回声信号的幅度。
由此存在这样一种情况,其中发现常规方法因为灵敏度低或是特异性差而无法获取可靠的病状诊断,由此将会受到限制。但是,射频信号不仅包含了与幅度有关的信息,而且还包含了频率和相位信息(该信息与回声因为散射体个体差异而延迟到达传感器的时延相关联)。
最新的信号处理方法倾向于提取用于表征信号的所有参数,尤其是通常所说的频谱内容。在从传感器返回回声信号的时候,直接施加于射频信号的频谱分析可以从回声信号的所有特征中提取信息。
用于“读取”频域中的信号的频谱分析方法可以提供特征值,以便对不同类型的信号加以区分,进而对与之交互的生物结构加以区分。
图4中的图示(A)再次显示了图1中的频谱,而图示(B)和(C)显示的则是与图4A中的描迹的两个部分有关的信号的频谱内容,由此显示了组织中的结构差异如何导致产生不同类型的超声波脉冲,进而产生相应频谱中的不同频率变化。图4B和4C所显示的信号的这两个部分的频谱内容是从0Hz延伸到12MHz的。应该指出的是,这两个频谱具有邻近4MHz(该信号涉及的是用标称频率为5MHz的超声波探测器实施的扫描操作)的大致相同的中心频率(可以发现峰值的频率),但却具有不同的频率变化,这是因为第一个图形具有成分(大于5MHz),这些成分在信号的其他部分中较少出现。
在实践中,由于传播路径存在差异,因此,从微观结构反向散射的超声波信号会在接收传感器表面以一种破坏性或建设性方式产生干涉,由此导致图像中出现“斑点”效应.对取决于干涉类型的回声信号幅度差异、即回声特性(echogenicity)同样会受到组织反射率的影响,而所述反射率是由微观结构成分的机械特性、密度以及弹性确定的;其中举例来说,含量很高的胶原质将会增大局部机械阻抗,由此将会提升回声特性.因此,以下特性有助于展现典型的超声波图像:确定斑点特性的微观结构的空间结块,以及用于确定回声特性的微观结构的机械特性.这样一来,受随机构成的各种形态因素的影响,空间上的平均幅度变化似乎并没有提供组织的单义特征.
如果可以在斑点上提取有关空间周期性(空间频率)的信息,那么可以获取与被检查结块的任何空间周期性有关的信息。
射频超声波信号在超声波信号的特定幅度和相位调制方面将会受到图像斑点所产生的信息的影响。
在时域中存在的信息在频域中被保留,并且斑点的影响将会转变成轴向变化效应,这种效应包含了依照图4D所示图示的幅度变化。实际上,以超声波束碰撞N个散射体或反射体(它们是结合超声波信号波长来区分的,但是这种区别并不是非常明确的)为例,对依照傅立叶变换所计算的接收回声信号而言,其功率谱密度|X(ω)|2可以表示如下:
| X ( ω ) | 2 = | S ( ω ) | 2 [ Σ i = 1 N r i 2 + 2 Σ i = 1 N - 1 Σ j = i + 1 N r i r j cos ω ( 2 ( l i - l j ) / c ) ]
其中:
|S(ω)|2:所传送的信号的功率谱密度
ri:第i个反射体的反射系数
li:第一个与第i个反射体之间的距离
c:超声速度
如果只有两个反射体并且这两个反射体之间的距离等于1,那么该表达式可以采用一种更简单的形式:
|X(ω)|2=|S(ω)|2[r1 2+r2 2+2r1r2 cos(2ωl/c)]
是由两个反射体反向传播的信号的峰值之间的时间距离。
从这个示意性公式中可以看出,频谱是用周期与反射体的距离成比例的余弦函数来加权的;这样则可以确定频谱功率密度的最大和最小值。这种效应被称为轴向变化,在图4D中为两个在超声波束传播方向上相互距离(用S0、S1、S2表示)可变的反射体示意性显示了这种效应。
为了区别两个无法以其他方式加以区分的反射体,也就是说,这两个反射体的距离小于所传送的脉冲的持续时间,我们可以在频谱域中评估最小值之间的距离。对多个反射体而言,初始信号频谱被周期不同的重叠振荡加权。此外,上述公式还包括确定反向散射信号幅度的反射体的反射系数。
从上述公式中可以看出,时域中存在的信息影响频谱形状,因此应该将关于时域幅度分布变化的精确描述提供给那些用于区分生物结构的部件。
如这些示范性考虑因素所示,射频信号的频谱分析是一种借助相关超声波信号来表征组织的强有力手段。
近年来,随着新型信号处理算法的开发,有可能依照反射信号中感应的频谱变化来调查组织与超声波束之间的交互,由此获取更多的组织信息。实际上,目前已经发现,病症结构可以依照其选择性频率过滤作用的差别来识别。
然而,迄今为止,已开发的方法尚未被证实有效,这其中的原因要么是因为这些方法不适合可靠地表征组织,要么是因为所用算法需要过多的计算能力,这样则无法以一种能够实时显示的方式来执行所述算法,由此不允许那些将数据输出速率作为基本需要的临床应用。
EP-A-1,341,003描述了一种用于对超声波信号进行频谱分析的全新方法,该方法克服了现有技术中的限制和缺陷,并且允许提取产生图像斑点的信息。在本说明书中,这个在先公开出版物的内容将被全部引入并构成说明书的一个构成部分。
发明内容
由EP-A-1,341,003描述的方法开始,本发明的目的是改进和扩展这类频谱分析所提供的可能性,从而获取与被检查结构有关的更重要和更精确的信息。
实质上,对于这里描述的发明而言,该发明的前提即为EP-A-1,341,003中描述的方法的基础。这里描述的发明允许借助以下处理来改进频谱描述:恰当选择母小波(wavelet)和/或选择将信号的整个谱带分解为数量可变并具有任意宽度和位置的子带,由此从被调查的采样中提取尽可能多的信息和/或更加有效使用局部频谱描述中的系数特性。
在实践中,依照第一个方面,本发明是以一种组合局部估计量的思想为基础的,在下文中也将这些局部估计量称为局部频谱参数,并且所述局部估计量可以借助于一种与EP-A-1,341,003中的描述相类似的过程,使用与所述局部估计量的各数值类别中的分布直方图形状相关的信息来获取。在实践中,根据从被检查的结构收集的射频信号帧,获取局部频谱参数或局部估计量的一个或多个帧或矩阵。随后,通过沿局部估计量的矩阵来移动具有预先确定的尺寸的窗口,可以对这些估计量进行统计处理,以便确定局部估计量在不同数值类别中的分布直方图。与它们的形状有关的参数则是从这些直方图中获取的,其中举例来说,该参数可以是标准偏差。然后,局部估计量连同分布直方图的形状参数将会组合在一起,其中所述组合可以是相乘、相除或是其他方式的处理。这样则可以获取一个新的局部加权估计量矩阵,以便允许将估计量的数值区间与所调查的采样的均匀区域相关联。然后,这些区域将会通过组织学分析来标识。原始超声帧的这些区域的显示是借助于重叠在常规B模式图像上的色码来进行的。
本发明的方法可以缩写成RULES方法(射频超声波局部估计量),因为该方法是以使用局部估计量为基础的,其中所述局部估计量包含了与被检查的结构所返回的超声波射频信号相关的频谱信息。
从下文中可以清楚了解,局部估计量是用恰当次数例如四次或五次内插多项式的系数来表示的,估计为先前细分数字化采样射频信号的不同频带所获取的系数。作为选择,局部估计量也可以包括内插多项式中不同次数的系数的组合。所述不同次数的系数的组合可以在用分布直方图的形状因数加权每一个系数之前执行,也可以在其后执行。
依照本发明的一个不同方面,本发明包括对内插多项式的系数进行处理,以便在没有使用分布直方图的形状参数来对系数进行加权的情况下,从所述多项式的不同次数的两个或多个系数的组合中获取局部估计量。此外,与现有技术中的方法相比,由于在这种情况下使用了内插多项式的若干个组合系数而不是单个系数,因此可以获取更多的信息。
对细分到频带、也就是分解到频带或子频带中的处理而言,该处理可以结合不同的准则来实现。此外,对局部估计量的分布直方图而言,其形状参数可以用多种方式进行选择。特别地,其中可以:
1.使用不同的母小波来适配所调查的组织的典型斑点;
2.使用不同的分解等级来获取不同的空间分辨率以及使用可以组合获取的整个频谱的不同细分;
3.借助于通过特定参数来对直方图进行针对性调查,使用内插多项式的系数组合,其中所述系数由涉及它们的统计分布特性的“参数”加权,并且举例来说,所述特定参数可以是:
-总体最多的类别
-标准偏差
-峰态(kurtosis)
-较高次数时刻
-对称索引
恰当加权的系数的组合为已处理的超声图像的每个窗口(可以设置该窗口的大小)创建一个新的矩阵,并且该矩阵将会包含这些组合。新的统计是在该矩阵上执行的,并且在这里将会选择与所要表征的各种结构相关联的类别。为了提升区别组织的能力,也就是识别“典型”类别的能力,在这里将会再次使用类别的组合分析。
特别地,在附图中显示了依照本发明的方法的有利实施例,并且在下文中将会参考某些实施例例示以及所执行的实验测试来对其进行更详细的描述。
依照一个不同方面,不同于执行局部估计量与所述估计量的分布直方图的形状参数的组合或加权操作,依照本发明,两个或多个局部估计量可以组合在一起,这种组合可以被用作组合局部估计量。例如可以产生若干个局部估计量矩阵,在每一个矩阵中,各系数由内插多项式中的不同次数的系数构成。然后,这些矩阵的相应位置中的系数可以被组合在一起,例如将所有矩阵或某些矩阵的系数相加或相乘。最后,借助这种方式,可以获取一个系数矩阵,其中每一个由先前获取的矩阵中所包含的局部估计量的函数组成。
附图说明
现在,我们将会从描述和附图中更清楚地理解本发明,其中所述描述和附图显示的是本发明的非限制性实际实施例。更具体的说:
图1至4A-4D显示的是超声描迹以及该描迹的两个部分中的频谱内容表示;
图5显示的是超声扫描设备的简化框图;
图6A显示的是由传感器发送到所检查的器官内部的超声波激励信号;
图6B至6D分别显示的是RF响应信号、用于过滤数字化RF信号的Daubechies 16小波以及通过分解成八个频带而获取的DWPT的系数;
图7A~7I示意性显示了依照本发明的方法而对超声波信号执行的处理序列;
图7J、7K示意性显示了依照本发明的处理过程的修改实施例;
图8A和8B分别显示的是颈动脉和前列腺的B模式超声波图像,其中图9A~17B中显示的结果是在该图像上获取的;
图9A~9C显示的是特别用于颈动脉检查的三个四次内插多项式变化,其中这三个多项式变化涉及的是脂质斑、血液以及钙化;
图10显示的是与颈动脉的脂质斑部位相关的四次内插多项式的各阶系数的分布直方图;
图11显示的是血液中内插多项式系数的分布直方图;
图12显示的是钙化区域中内插多项式系数的分布直方图;
图13A~13C显示的是通过与图9A~9C中的应用相同的应用而获取的五阶内插多项式变化;
图14、15和16与图10、11和12相似,显示的是脂质班、血液以及钙化的不同区域的五次内插多项式的不同系数的分布直方图;
图17A和17B显示的是用于前列腺分析的四阶内插多项式系数的直方图;
图18A显示的是在具有可变幅度以及低频区域细分增多的频带中,借助离散小波将经过数字化和采样的射频信号分解到六个频带或子频带中的级联滤波器的框图;
图18B示意性显示了六个分解频带;
图18C显示的是通过将图18A、18B中的分解处理应用于颈动脉的超声波响应信号而获取的四阶内插多项式变化;
图19、20和21显示的是受脂质斑和钙化现象影响的颈动脉的一部分的B模式超声图像,此外还分别显示了受到脂质斑、血液(游离血管)和钙化影响的三个部位中内插多项式变化以及内插多项式中不同次数的系数的分布直方图;
图22显示的是借助DWPT变换获取的分解树的结构;以及
图23A和23B显示的是局部估计量的分类过程的说明性图示。
具体实施方式
图5显示的是超声扫描设备的高度示意性框图。数字1表示的是用2概括性指示的扫描头所具有的传感器,其中该传感器发出超声波激励信号并且接收超声波响应信号。扫描头的模拟输出端3与采集卡5相连。后者将会产生经过采样和数字化的射频信号,该信号则构成了用7概括性表示的集成FIR滤波器组的输入。滤波器组7可以包括任何商用集成滤波器,例如Harris生产的43168电路、Gennum生产的GF191电路、Plessey生产的GEC16256电路或是Graychip生产的GC2011电路。FIR滤波器7与用9表示的数字信号处理器(DSP)相关联,该处理器充当整个处理单元的控制器。其中举例来说,所述微控制器9可以是Texas生产的TMS32031DSP,也可以是其他等价产品。
数字11表示的是与滤波器7以及控制器9相关联的存储器。存储器11与模数转换器13以及扫描头2的数字级相连。
这里描述的电路是已知的,并且在这里不必对其进行更详细的描述。实质上,该设备是以如下方式工作的。传感器1向包含所检查器官的身体的内部发出一系列具有如图6A所示形式的超声波激励信号,其中所述器官可以是前列腺、肝或其他部位。受到传感器1所发射的信号超声波撞击的组织将会返回一个散射或反射信号,并且该信号具有图1的信号格式(射频信号)。
所述射频信号由紧随在时间增益补偿(TGC)单元之后的超声波扫描器的接收电路获取,然后由采集卡5在例如40MHz的恰当采样频率上以12比特的分辨率(或是其他恰当的分辨率)进行数字化处理.在实践中,其中将会为单独的描迹执行采集处理,直至完成一个帧或是帧的一部分.为每一个单独描迹所获取的信号采样将被保存在采集卡的存储体中,直至采集到一个完整的帧或是帧的一部分为止;只有在这个时候,所保存的数据才会发送到包含了DSP 9以及FIR滤波器7的处理机卡,其中所述处理机卡将会依照本发明的方法来执行下文中详细描述的操作。在将数据传送到处理机卡的同时,采集卡将会开始获取下一个帧的描迹信号采样。而处理机卡12的处理时间则足够短,以便允许在下一个帧的采样信号采集时间以内结束处理。由此可以实现实时操作。
处理机卡12执行信号处理,尤其是使用集成滤波器7的时间-频率变换以及下文中更详细描述的其他操作,来获取信号频谱的表示特征。
图7A~7E以图表方式概述了针对输入信号所执行的操作,其中该操作被用于获取矩阵表示,所述矩阵包含与使用依照本发明的方法所确定的局部估计量的恰当组合相关的信息。为了对在具体实例中的不同处理阶段所执行的操作进行详细描述,现在将会参考EP-A-1,341,003中的描述,但是在这里不会对其进行详细描述。
对采集卡5所获取的经过采样和数字化的信号而言,针对该信号所进行的处理是从应用时间-频率变换以便将信号编码到子频带中开始的。本质上,对经过采集、采样和数字化的帧的各个轨迹而言,通过使用时间-频率变换、特别是借助离散小波分组分解而将信号分解到一系列的频带或子频带中,其中所述时间-频率变换诸如小波变换。滤波处理由滤波器组7执行。用以分解数字化信号的频带或子频带数量取决于所使用的连续滤波操作的数量。在下文的说明中,其中将会参考一种分解到八个频带或子频带的解码过程,并且在下文中,这些频带或子频带是用频带0、频带1、......、频带7表示的。然而,在下文中应该清楚了解,该值并不是约束;与之相反,对本发明的方法而言,其创新方面在于能够通过解码到具有不同宽度(也就是说,幅度是随频带而改变的)的大量频带或子频带而从射频超声波信号中提取更多用于表征组织的信息。
滤波系数可以通过使用Daubechies 16小波来定义。举例来说,图6C显示了Daubechies 16小波的变化,并且图6D显示了通过对图6B中的经过数字化和采样的RF信号进行滤波而获取的DWPT系数的变化。然而,所使用的时间-频率变换的类型有可能作为所分析的结构或组织类型的函数而改变。例如,在这里可以对Daubechies 16使用不同的小波。一般来说,如果要对颈动脉进行检查,那么可以使用双正交的3.7小波来表征脂质斑,在下文中将会参考某些实验数据来对此进行说明。
滤波器系数的数量一般很少,例如对Daubechies 16滤波器而言,其系数是32个,因此,这些系数很容易在商用电路中实现,所述电路通常提供系数数量有限的滤波器,如果该滤波器不是线性相位类型的滤波器,同时可以借助级联的两个部分来获取最多128个系数,那么每一个部分的系数不会超过64个,并且该滤波器是按4抽选的。对提供小波变换的滤波器而言,其系数以递归方式定义并且是已知的。例如,在这里可以使用US-A-6066098中定义的系数,其中该文献的内容将被引入到本说明书中。
图7A示意性显示了由超声扫描器的扫描头采集的帧.在下文中,这个帧称为“输入帧”,它包含了数量为n的多个轨迹.每一个轨迹都包括由传播超声激励信号的组织反射或反向传播的波前表示的模拟信号.在执行采样和数字化处理之后,这时将会得到包含了n个轨迹的帧(图7B),其中每一个轨迹都包括一系列数字化信号采样.这个帧则被称为“数字化帧”.
通过将时间-频率变换应用于数字化帧,也就是说,通过对包含形成数字化帧(在图7B中示意性显示)的系数集合的采样数字化信号进行过滤,获取频带或子频带中的分布:在实践中,对从1到n的各个轨迹而言,在每一个时刻,也就是对时间方向(也就是信号传播方向)上的每一个采样而言,在由图7C示意性显示并用频带0、频带1、......、频带7表示的每个频带中都有系数。
在所示实例中,其中有八个频带并且所有这些频带都具有相同的宽度。然而如上所述,有必要理解的是(在下文中将会参考某些实验数据来进行说明),用以分解信号的频带的数量也可以更多。在具有数量更多的解码子频带的情况下,可以从射频信号中获取与其频谱特性相关的更多信息,并且由此获取更多关于被调查结构的信息。此外,细分到子频带中也可以借助不恒定的频率间隔来实现,也就是在幅度可变以及具有任何形式的滤波器的频带或子频带中实现。分解类型是遵循实验标准选择的,其中所述实验标准作为受到超声调查的组织或其他结构的类型的函数。
在下文中,包含时间-频率变换系数的矩阵是作为“变换矩阵”表示的;所述矩阵的系数则被表示为“频谱系数”。这是三维矩阵,该矩阵的维度由轨迹数量、每个轨迹的采样数量、用以细分射频信号的频带或子频带数量给出的。
后续操作包括:为频谱系数矩阵的每一个点(由此为输入帧的各点,其中频谱系数矩阵是借助上文中说明的操作而从所述输入帧获取的)提取信号的局部估计量。
为此目的,对于每一个时刻、也就是沿着时间维度(也就是深度)的每一个点,沿着每一个单独的轨迹来获取不同频带0~7的系数,并且确定近似这些系数的变化的内插多项式。在图7D、7E中示意性显示了这个操作。更为特别的是,图7D以常规时刻pk为例显示了与轨迹1到轨迹n相关的不同子频带中的频谱系数。在图7E的第一和第二个图示中,不同频带0~7的所有系数在位置I,k和n,k都被显示成是相互对准的,其中PI1表示内插多项式,也就是对在轨迹1上在深度k处对准的频谱系数进行最优近似的多项式,而PIn表示与轨迹n上在深度k处对准的频谱系数相关的内插多项式。
这个多项式的一个或多个系数构成局部估计量,其中如下所述,所述局部估计量允许提取与射频信号频谱相关的必要信息。实际上,这个参数表示的是信号频谱特性的变化。这是因为依照定义,该参数涉及的是为频谱系数矩阵中的每个点近似每个频带中的频谱系数变化的多项式的属性,由此间接涉及输入帧的属性。另一方面,从内插多项式PI提取的局部频谱估计量提供了在输入帧的指定点的RF信号频谱特性变化的定性表示。
由于内插多项式是在输入超声帧的每一个点计算的,因此这个操作的结果——对内插多项式的每一个系数而言——将会是包含频谱估计量集合的矩阵。在下文中,该矩阵被称为“局部估计量矩阵”或“局部估计量帧”,并且在图7F中显示了用于内插多项式的单个系数的矩阵。在这里,局部估计量的各值是用aij表示的,其中i从1变化到n,n是矩阵中的轨迹数量;j则从1变化到m,m是每一个轨迹的系数数量.有必要了解的是,根据内插多项式的不同阶的各系数的意义,有可能为多项式中的每个系数创建局部估计量的矩阵.因此,对四阶内插多项式而言,其中有可能构造五个之多的局部估计量矩阵,并且每一个矩阵用于多项式的五个系数中的一个.可以使用更高次数的多项式来构造数量更多的矩阵.然而在单个超声波信号的频谱表征过程中,并非所有的多项式系数都是必须使用的.对这些系数而言,其意义重要与否是依照所调查的组织或结构类型而以实验方式确定的.在下文中将会对某些非限制性实例进行描述.
最小二乘法近似可以用于计算多项式的系数值。在给出了数量为n并且对应于横坐标x0、x1......、xn-1的纵坐标y0、y1、......、yn-1的情况下,这种方法可以用于确定次数m<n-1的多项式:
P(x)=a0+a1x+......+amxm
该多项式给出了被检测数据的最佳近似值。
考虑近似值与观测值之间差值的平方和S
S = Σ j = 0 n - 1 ( a 0 + a 1 x j + . . . . . + a m x j m - y j ) 2
这些系数必须以使S最小的方式来进行选择,以使多项式能够尽可能精密地估计数据。通过相对于系数导出S以及将结果均衡为零,在给出了以及
Figure G200480035699XD00143
的情况下,这时将会得到一个等式系统,其可以采用矩阵形成而被写为:
a=Q-1v
其中a是多项式P(x)的系数矢量。
对不同类型的组织以及不同类型的病症状态来说,可以观察到的是,内插多项式的系数值采取以统计方式确定的值。通过定义数值类别,可以进一步观察到统计分布具有依赖于组织类型及其可能具有的任何病症状态的可变形状,其中所述统计分布即为内插多项式的每一个系数的分布直方图。依照本发明的一个方面,该方法是以这个发现为基础的,并且该方法使用了直方图的形状参数来从超声波信号中提取数量更多的频谱信息。
为了理解在不同状态以及对于不同组织直方图形状的变化的重要性和程度,有必要对图8A至17B加以参考。
举例来说,在对颈动脉进行分析以便检查颈动脉斑是否存在或是检查其稠度的过程中发现:通过使用一种借助双正交的3.7小波而解码到十六个频带的处理以及通过使用四阶内插多项式,不但多项式的三阶系数具有用于表征被调查组织的有效分布,而且多项式的其他系数同样具有用于表征被调查组织的有效分布。此外还可以发现,通过分析分布直方图可以显示:分布直方图的不同形态取决于接受检查的不同类型的组织。
图8A显示的是用于获取下述实验结果的颈动脉的B模式超声图像。此外,在该图像上还标识了三个窗口,这三个窗口则是以脂质斑、游离的血管部位以及钙化为中心的,其中所述游离血管部位即为血液。局部估计量矩阵的系数是在这些组织结构互不相同的部位中计算的,由此显示这些估计量具有依照所调查的组织类型而改变的统计分布。
图9A、9B和9C描绘的是三个四阶内插多项式的变化,这些变化分别涉及的是斑块、血液和钙化。也就是说,在遵循上文所述的数字化射频信号的解码和处理方法的情况下,对涉及脂质斑的图像部分而言,对频谱系数(也就是变换矩阵的系数)进行最优近似的四阶多项式具有图9A中的变化,对游离血管部位而言,该多项式具有图9B中的变化,对存在钙化现象的部位而言,该多项式具有图9C中的变化。
图10、11和12分别显示的图9A、9B和9C中的内插多项式的五个系数在不同的参考类别中所具有的统计分布直方图,其中所述内插多项式是通过分析颈动脉斑的一部分的不同部位获得的.参考类别是在横坐标上表示的,总体频率则是在纵坐标上表示的.
对所考虑的五个系数(在这里用A0、A1、A2、A3、A4表示)中的每一个系数而言,其中构造了与三个不同的窗口相对应的三个直方图,这三个窗口分别位于脂质斑部位(图10)、游离血管部位(图11)以及遭受钙化的部位(图12),其中所述游离血管部位即为充满血液的部位。从图11中的直方图可以看出,在游离血管部位中,所有次数的系数分布都是非常窄的。与之相关的分布直方图则具有极其有限的分散。相反,与健康组织(血液,图11)中的相同系数相比,在受脂质斑和钙化影响的部位(图10和12)中,二阶系数(A2)特别具有完全不同的统计分布,此外,对脂质斑以及钙化而言,所述二阶系数是完全不同的。换句话说,通过对用于表示二阶系数统计分布的直方图的变化加以分析,可以区分所检查的部位包含的是游离血管(血液,图11)、脂质斑(图10)还是钙化(图12)。
这意味着脂质斑和钙化会在反向传播的信号中引入不同的频谱变动,这种变动可以使用内插多项式的二阶系数的分布直方图的定性评估来加以检测,其中所述内插多项式近似频谱系数的变化。
相似考虑可以应用于与三阶系数相关的直方图的变化,其中所述三阶系数在所考虑的三个窗口中具有不同的分布,其中这三个窗口分别是以脂质斑、血液以及钙化为中心的。涉及一阶系数分布的直方图仅仅存在程度较小的差别,但是这种差别仍旧是非常明显的。
通过使用一种解码到数量众多的频带的处理,可以使用次数大于四的内插多项式,例如五阶内插多项式,而在上述附图所指示的数据实例中,所述频带的数量则是十六个。这样一来,甚至那些与不同系数的统计分析相关联的更重要的实验结果也是可以获取的。图13A、13B和13C显示了五阶内插多项式的变化,该多项式近似颈动脉的频谱系数变化,其中所述变化是借助使用双正交3.7小波解码到16个子频带来获取的。
图14、15和16中显示的直方图描绘的是用于图8A中所示窗口的内插多项式的不同次数的系数(A0~A5)的统计分布,其中所述窗口分别是以脂质斑、血液以及受钙化影响的部位为中心的。
特别地,通过将图9A、9B、9C与13A、13B、13C相比较,可以在脂质斑和钙化部位观察到内插多项式的不同变化。内插多项式系数的直方图在不同部位中的不同变化与此是对应的,其中与涉及四阶多项式的直方图中发现的差别相比,不同类型的组织之间的统计分布差别甚至更大。
同样,所关注的结果是在关于前列腺的超声检查中获取的。在图8B中显示了下文所要阐明和论述的用于获取实验结果的B模式超声图像。此外,该图像还显示了以所检查的肿瘤组织部位为中心的窗口。
同样,对前列腺而言,使用四阶内插多项式以及将数字化射频信号解码到16个频带或子频带,记录内插多项式的一阶、二阶、三阶以及四阶系数的统计分布,其中这些分布在健康组织与肿瘤组织之间是存在差别的,对分布直方图的形状而言也是如此.图17A和17B分别为肿瘤组织和健康组织显示了用于一阶、二阶、三阶以及四阶阶数的系数a1和a2(图17A)以及系数a3和a4(图17B)的不同参考类别的总体频率直方图,此外,还包括不存在组织的情况(在水采样中测得的噪声).通过比较这些直方图可以看出,由于对所有的四个系数而言,不同系数的直方图形式会在健康组织与瘤组织之间发生变化,因此所有系数都提供了可以用于识别肿瘤组织的信息.
一方面,这些实验结果显示了内插多项式的不同系数的相关性(用于频谱分析),另一方面则显示了局部估计量的不同数值类别中分布直方图的形状因数在用于识别肿瘤组织或被检查器官中的特定特征中的重要性。
因此,依照第一个方面,根据本发明的方法,内插多项式的一个或多个系数的分布直方图的形状因数被用于获取可以在超声扫描器的正常B模式图像上重叠显示的信息。
依照一个不同的方面,本发明的方法其特征在于使用若干局部估计量,也就是使用同一内插多项式的不同次数的若干系数。
在下文中将参考图7G~7I来描述一种可能的处理方法。该实例考虑单个局部估计量矩阵,也就是内插多项式的单个系数的值。所考虑的局部估计量矩阵可以是与内插多项式的任何系数有关的矩阵,而所述内插多项式在调查被检查的组织类型的过程中是非常重要的。对所关注的每一个局部估计量矩阵而言,下文描述的程序是可以重复执行的,其中每一个矩阵都是通过收集那些在超声图像中通过内插多项式的特定系数而取的数值来获得的。此外,还可以以一个局部估计量矩阵开始,其中该矩阵包含了与同一多项式的不同次数的系数相关联的矩阵的组合。举例来说,在这里可以产生一个收集二阶系数的局部估计量矩阵,以及一个收集三阶系数的不同的局部估计量矩阵,然后则从这两个矩阵产生一个局部估计量矩阵,在该矩阵中,每一个系数都是所述前两个矩阵的相应系数之和。在一个不同的实施例中,其中也可以使用在两个或多个局部估计量矩阵中收集的局部估计量的乘积作为局部估计量。此外,使用内插多项式的所有系数同样是可行的。对四阶多项式而言,其中可以将五个系数组合在一起(将其相加或相乘,或是对其施加其他类型的处理),以便从五个矩阵中获取单个矩阵。
返回到所描述的实例,该实例使用了内插多项式的单个系数,并且由此使用单个矩阵,其中该实例是从图7F中的局部估计量矩阵开始的,而局部估计量的数值的分布直方图则是在不同的参考类别中确定的。
使一个具有恰当定义的大小I×J的窗口从所考虑的局部估计量矩阵上移过。该窗口位于第一位置(图7G),而落入该窗口内的局部估计量矩阵系数aij的分布直方图则是在这个位置中计算的。
如参考图10以及下文所述,相对于实验结果而言,对给定多项式的不同阶数而言,其系数的分布直方图将会根据窗口所包围的被检查图像的特定部分中发现的组织类型而发生变化,从而采用不同的形状。一旦将大小为I×J的调查窗口定位在局部估计量矩阵上的给定位置,那么将会完成分布直方图,并且计算所述直方图的一个形状因数,例如标准偏差σ。通过将所述窗口沿着局部估计量矩阵的行与列一次移动一步,来获取一系列与局部估计量矩阵的系数数量相等的标准偏差值。通过用σij来表示与使用位置ij的窗口计算的局部估计量矩阵系数的分布直方图有关的标准偏差值,可以构造出一个系数数量与局部估计量矩阵系数aij的数量相等的矩阵,但是该矩阵的值在考虑分布直方图的形状参数(在这个具体实例中则是“标准偏差σij”)的情况下进行了修改和处理。这个矩阵将会包括m×n个值bij
bij=f(aij,σij)
其中f是aij和σij的类函数。举例来说,局部估计量矩阵的每一个系数aij都可以与相应值σij相除,由此可以得到:
bij=aijij
然而有必要理解的是,将系数bij关联于系数aij以及相应分布直方图的形状因数σij的函数也可以不同于上述函数。例如,在这里可以执行乘法而不是除法,此外还可以执行这些数值的不同组合。
此外,在这里还可以为内插多项式中的不同次数的多个系数考虑相应的多个局部估计量矩阵,并且将这些矩阵组合在一起,使用相关分布直方图的形状参数来对其进行加权。
所用直方图的形状因数或形状参数也可以不同于标准偏差。举例来说,在这里可以使用峰态、对称索引、较高次数时刻或是用于表示对象数量最多的类别的参数。通常,局部估计量矩阵的系数可以通过任何在特定窗口中表征所述矩阵系数的分布直方图的形状因数来处理,其中举例来说,该窗口以局部估计量矩阵的每一个单独系数的位置为中心。
在下文中将大小为m×n、收集系数bij的矩阵称为加权局部估计量矩阵。这意味着所述矩阵的每个系数被相应窗口的分布直方图的形状因数加权,其中“加权”意味着以某种方式将该系数与形状参数相组合,所述形状参数可以是标准偏差或另一恰当参数。
有必要理解的是,对可用于超声波信号的频谱表征的内插多项式的不同次数而言,在这里可以为这些次数的某些或所有系数单独执行所描述的操作。对收集内插多项式的一般次数的系数的单个局部估计量矩阵来说,针对该矩阵所执行的操作对所有那些被考虑的次数而言可以是相同的,尽管也可以用不同的方式来对内插多项式的不同次数的系数进行处理。例如,通过将三阶系数与标准偏差参数相除,可以对三阶系数进行加权,而二阶系数既可以进行加权,也可以执行另一种与峰态之类的不同参数相结合的操作,此外还可以执行一种不同于与局部估计量以及形状参数相除的操作。
无论如何,一般概念是组合表示局部估计量的系数值和考虑特定窗口中所述系数的直方图的形状的系数。
图7J示意性指示了在使用所述内插多项式的所有系数而不是单个系数来获取更复杂和更精细的加权局部估计量集合时所执行的过程。由此,局部估计量矩阵是一个大小为m*n*k的三维矩阵,其中m是每一个轨迹上的采样数量,n是经过数字化和采样的超声帧中的轨迹数量,k则是内插多项式(次数为k-1的多项式)系数的数量。该矩阵的系数是用aij (1)、......、aij (k)表示的。
系数aij (1)、aij (2)、......、aij (k)将会相互组合,以便获得组合局部估计量矩阵的系数cij。这个矩阵是一个大小为n*m的二维矩阵。它的每一个系数都是如下给出的:
cij=f(aij (1),aij (2),......,aij (k))
其中f是类函数。举例来说,该函数可以是三维矩阵的位置ij的所有k个系数的总和或乘积。
在这个矩阵cij中,分布直方图的形状系数σij可以是标准偏差,并且该系数是借助上文中参考图7G、7H所描述的过程获取的。加权局部估计量矩阵同样是用bij表示的,由此,总体矩阵fij是从形状系数cij以及形状参数σij中获取的。
在概念方面,该过程与参考图7D~7H所描述的过程相似,但不用之处在于:包含源自内插多项式的所有系数而不是其中一个系数的信息的系数(加权局部估计量)包含在最终矩阵中.
仅仅使用内插多项式中的某些系数的中间解决方案同样是可以实现的。
在一个实施例变化中,其中首先可以通过形状系数σij来执行单个系数aij的加权,随后则可以仅仅将经过加权的局部估计量相互组合。实质上,在这种情况下,当确定了频谱系数矩阵之后,确定内插多项式,此外还会从这些多项式的系数中确定估计量aij (k)的矩阵。对每一个k值来说,确定作为直方图形状的统计分布,并且计算一个形状因数,例如标准偏差σij (k)。然后计算加权的局部估计量:
bij (k)=f(aij (k),σij (k))
最后,组合加权的局部估计量而得到:
cij=f’(bij (1),bij (2),......,bij (k))
其中f和f’是类函数。
依照本发明的一个不同方面,该过程并没有对系数或局部估计量进行加权。所述过程的这个不同实施例不允许使用分布直方图的形状中所包含的信息,但是仍旧允许使用内插多项式的若干个系数的信息内容。图7K示意性描绘了这种情况下的处理过程。其中相同的符号表示与图7J中的部件相同或等价的部件。在这里不再通过处理矩阵cij来获取形状因数σij,但是其系数(未经加权的局部估计量)会经过统计分析,以便直接获取在这里用fij’表示的总体矩阵。
在获取了加权局部估计量矩阵的情况下,这个矩阵的系数值将会依照数值类别来进行分类,其中所述数值类别是通过实验发现的与特定组织结构具有双单义关联的数值类别。换句话说,通过实验可以证实,某一种类型的组织(例如前列腺的瘤组织)将会改变超声波激励信号,由此某一个经过加权的局部估计量将会采取落入某个数值类别的值。然后,通过对超声扫描器进行编程,可以对输入帧执行上文描述的操作,从而获取加权局部估计量矩阵。之后,以这种方式获取的矩阵系数将会与瘤组织的数值类别特性进行比较,其中对那些数值落入瘤组织的特性类别内的加权局部估计量而言,由这些加权局部估计量矩阵的点形成的彩色图像将会重叠在黑白的B模式图像上。这个彩色指示为超声操作者提供了必要信息,由此显示了瘤组织的存在与否及其位置。
应该理解的是,在这里可以为不同类型的组织或结构执行分类。举例来说,如果希望显示健康组织的存在,那么通过识别某种用于表征健康组织存在的加权局部估计量的数值类别即可满足需要。同样,如果希望在相同的图像上同时或者按顺序显示两种不同类型的组织或结构(例如脂质斑和颈动脉中的钙化)的存在,那么只要执行以下操作即可满足需要:识别表征这两种结构类型的两种数值类别、然后使用这些类别的特征值来过滤局部估计量矩阵的系数,以及同时或是按顺序在B模式超声图像上显示结果,其中在必要以及优选的情况下,所述显示是用两种不同的颜色实现的。
对于为某个局部估计量表征特定类型组织的数值类别而言,要想以实验方式来识别该数值类别,可以通过检查实验采样来执行该处理,其中在所述实验采样中借助了直方图来分析所关注的结构。以前列腺为例,其中将会对健康组织、腺癌以及腺纤维瘤的均匀部位进行识别。然后,如上所述,作为响应而从实验采样扫描中获取的射频超声波信号将被处理,以便获取单个加权局部估计量矩阵。
B模式图像的点与矩阵的每个系数是对应的,这一点归因于用于获取局部估计量矩阵的方式。因此,由于某种结构(例如腺癌)所处的部位是已知的,在这些部位中确定局部估计量的统计分布,并且识别出具有最大总体频率的类别。就这个特定器官而言,该类别与腺癌的存在与否是相关联的。因此,在每一次对类似的器官执行超声分析的时候,其中都会对经过加权的局部估计量进行分类,并且那些落入腺癌的参考类别内的加权局部估计量将会通过在B模式图像上覆盖色标来高亮显示,其中在所述B模式图像中,对其中加权局部估计量数值落入腺癌特征类别内的图像区域而言,与之相对应的象素是彩色的。此外,相似的过程还被用于其他类型的组织(对前列腺而言则是健康组织、腺纤维瘤)。
更详细的处理设想的是:在对参考类别进行识别的阶段以及对接受超声调查的器官进行分析的阶段,作用于至少一对局部估计量矩阵之上。现在将借助图23A和23B来描述基于三个不同局部估计量的此类处理的实例。对实验采样而言,在这里将会对其再次进行超声帧采集处理,然后,该采样将被实施上文所述的处理,从而获取三个加权局部估计量矩阵。
第一个二维图形分别在横坐标和纵坐标上显示了两个加权局部估计量矩阵的所有系数值,例如通过内插多项式的零阶和一阶系数所获取的值。第二个图形显示的是借助零阶和二阶系数所获取的加权局部估计量。图23A在横坐标上显示了借助用于前列腺的多项式的一阶系数而获取的加权局部估计量矩阵的系数。在纵坐标上则显示了借助同一多项式的零阶系数所获取的加权局部估计量矩阵的系数。图23B在横坐标上显示了借助二阶系数获取的加权局部估计量,并且在纵坐标上显示了借助零阶系数所获取的加权局部估计量。
因此,由所使用的加权局部估计量中的某一个值以及其他值组成的一对数值(点的坐标)是与这两个图形中每一个图形的每一个点相关联的。
以这种方式获取的图形可以重叠在输入超声帧上。因此,输入超声帧的每一个点可以关联于那些与该点相对应的被检查加权局部估计量的数值类别对。
借助组织学分析,在超声帧上将会标识那些特性(健康、瘤等等)均匀的组织部位。在图23A和23B中,这些部位是用三个不同的灰度级区分的。在这两个图形中,由于加权局部估计量的数值对是与超声帧中的每一个点相对应的,因此可以将用于所述加权局部估计量的数值类别关联于所检查的组织中的每一个均匀部位,其中所述类别与这个特定类型的组织具有单义关联性。
为了以双单义方式来识别每一个部位,非常重要的是选择一个集合,该集合包含了加权局部估计量数值类别对或三元组,并且这些加权局部估计量并未重叠在与不同部位相关联的其他估计量上。其中举例来说,用于提取这些集合的操作可以通过逻辑“或”运算以及“与”运算来完成。
如所述,对于以这种方式确定的集合而言,其总体频率在常规的B模式图像上是用色标表示的,所述色标对每一个集合而言都是不同的,此外,所述色标直接与总体频率的强度成比例。
正如图17A、17B中以实验方式所指出的那样,对希望显示的部位(在本实例中则是瘤组织)而言,通常只要将其与那些从内插多项式的所有系数中获取的图像的总和相组合,就可以高亮显示该部位,在这种情况下,所述部位的高亮显示是通过借助逻辑“与”操作而与那些从内插多项式的所有系数中获取的数值集合进行组合来实现的.
对于未加权的局部估计量,例如通过将内插多项式的不同次数的系数相互组合所得到的局部估计量而言,用于在通过组织学分析而检测的组织特性与加权本地估计量的数值类别之间进行实验相关的上述过程同样是有效的。
对上文中论述以及图9A~17B中表示的实验结果而言,其中使用的是将RF信号分解到16个相同宽度的频带或子频带的处理。这种分解处理在前列腺检查以及颈动脉检查中提供了非常有用的结果。选择这两种类型的器官则是因为它们具有大小不同的组织结构。与颈动脉相比,前列腺的结构通常具有较大的尺寸。对于反射超声波激励信号的不同大小的结构而言,这些结构将会导致在返回的信号中出现不同的频率分布。而通过执行细分到16个相同幅度的频带的处理,可以获取从被检查的器官返回并且分辨率足够高的RF信号。
尽管如此,一般来说,以上描述可以了解,分解到频带中的选择(频带的数量和幅度)有可能依照所要分析的组织或结构类型而改变。在被分析的器官中,所研究的结构越小,分解到高频的分辨率也就越大,反之亦然。这一点可以通过下文中参考图18~21所论述的实验结果来加以确认,其中所述实验结果同样是通过处理颈动脉的某一部分的超声图像获取的。
图19、20和21显示的是所分析的颈动脉部分的B模式图像。使用一个离散小波将经过数字化和采样的RF信号解码到六个幅度可变的频带中。图18A显示的是描绘所使用的滤波处理的树形结构,图18B显示的则是六个分解频带。可以观察到,从低频到高频,频带的幅度增加。信号的采样频率则是用fc表示的。图18C显示了通过将前述分解处理应用于图19、20和21中的B模式图像而获取的内插多项式的变化。
图19的顶部显示的是所检查的颈动脉的B模式图像,其下方是四阶内插多项式的变化(与图18C相对应),再下来则是在以脂质斑为中心的窗口(在B模式图像上重叠显示)级计算的内插多项式的不同次数的系数A0-A4的分布直方图。具有相同图像布局的图20显示的是在以游离血管部位、即血液为中心的窗口中获取的内插多项式以及分布直方图。同样,在这个范例中,窗口也是在颈动脉的B模式图像上重叠显示的。最后,具有相同配置的图21显示的是在一个以B模式图像中被辨认出存在钙化现象的部分为中心的窗口中获取的结果。
通过比较图19、20和21可以看出,在分解到幅度从低频到高频逐渐递增的六个频带的情况下,在所观察的三个部位(脂质斑、血液以及钙化)中,内插多项式的变化并未显著改变。而这正是所用分解没有提供适合检查这类器官这一事实的第一个指示。出现这种情况是因为所识别的结构的尺寸很小,由此它们提供的频谱信息都位于响应频谱的高频部分,也就是说,在这种情况下,所使用的分解处理需要具有更大幅度的频带。与之相反,对分解到16个幅度相等的频带中的处理而言,其在相同器官上获取的结果(图9A~16)显示借助该分解处理所获取的分辨率足以提供有效结果。
从图19、20和21的分布直方图中可以注意到,实际上(虽然内插多项式的变化在脂质斑、血液和钙化这三个区域中基本上是相同的),不同系数的次数的分布直方图的形状是变化的.然而,在这三个被调查的部位中,相应系数的分布直方图之间的差别并不像使用基于解码到16个幅度相等的频带的信号处理所得到的差别那样显著和明显.最后,对借助两种解码到子频带的处理所获取的结果而言,对其结果所进行的比较显示:如果就频带数量及其幅度方面而选择了恰当的解码处理,那么则允许将本发明的分析方法适配成依照所调查的组织或结构类型来获取最优的结果.解码或分解类型可以用实验方式进行选择,此外还可以借助实验来选择最适合的内插多项式次数,其最重要的系数,以及所述系数的分布直方图的最相关形状系数.
图22例示了通过离散小波分组变换(DWPT)的射频信号频谱分解树。依照所执行的超声调查的类型以及所检查的结构的特征,该调查方法可以使用沿着图22中的树形结构的不同分解可能性中的任何一种。
很明显,附图仅仅显示的是本发明的实际实施例,在不脱离本发明的原理范围的情况下,这些实施例的形式和装置是可以发生变化的。附加权利要求中给出的任何参考数字的用途仅仅是为参考上文以及附图所进行的阅读提供帮助,它们并没有限制权利要求的范围。

Claims (27)

1.一种用于对从接受超声检查的结构返回的射频超声信号进行频谱分析的方法,该方法包括以下步骤:
a)将超声激励信号传送到被检查的所述结构的一部分;
b)接收来自所述结构的射频响应信号;
c)应用一系列滤波操作,以便获得将射频响应信号的频带分解成多个频带;
d)从所述滤波操作产生的系数计算局部估计量,所述局部估计量包含关于射频信号频谱的信息;
其中,所述局部估计量与表示所述局部估计量的统计分布形状的参数组合到超声图像的一部分中。
2.如权利要求1所述的方法,其中用于细分所述射频信号的频带覆盖射频信号的整个频带。
3.如权利要求1所述的方法,其中分解后的多个频带是具有不同宽度和位置的频带。
4.如权利要求1所述的方法,还包括以下步骤:
为超声输入帧产生经过采样和数字化的帧;
将所述经过采样和数字化的帧分解到所述频带中;
产生频谱系数矩阵,该矩阵包含所述滤波操作所产生的系数或从中推导的系数;
为经过采样和数字化的帧的至少一些点确定相应的内插多项式PI,这个多项式近似所述频谱系数在用以分解射频信号的不同频带中的变化;
为所述点从内插多项式的至少一个系数获取所述局部估计量,所述局部估计量构成局部估计量矩阵。
5.如权利要求4所述的方法,其中所述局部估计量的每一个由相应内插多项式的一个系数构成。
6.如权利要求4所述的方法,其中根据内插多项式的至少两个系数而为每一个点确定至少两个局部估计量,以便产生三维的局部估计量矩阵。
7.如权利要求4所述的方法,其中所述局部估计量的每一个由相应内插多项式的多个系数的组合来构成。
8.如权利要求4所述的方法,其中所述局部估计量的每一个与包含所述局部估计量的窗口内的所述局部估计量的分布直方图的形状系数相结合,以便获取加权的局部估计量。
9.如权利要求8所述的方法,还包括以下步骤:
确定比所述局部估计量矩阵尺寸更小的窗口中的所述局部估计量的统计分布;
为每一个所述窗口确定所述统计分布的形状参数特性;
为每一个窗口组合所述形状参数和相应局部估计量,以便获得加权的局部估计量。
10.如权利要求4所述的方法,其中使用相应内插多项式的不同系数而为经过采样和数字化处理的超声帧的同一个点获取的若干加权局部估计量相互组合。
11.如权利要求1所述的方法,其中所述滤波操作是通过使用时间-频率变换来获得的。
12.如权利要求11所述的方法,其中所述时间-频率变换是小波。
13.如权利要求11所述的方法,其中所述时间-频率变换是离散小波分组变换DWPT。
14.如权利要求8所述的方法,还包括以下阶段:确定加权局部估计量的统计分布,以及创建数值类别的集合,其中所述集合能够双单义地识别被采样的超声帧上的均匀部分.
15.如权利要求8所述的方法,其中使用所述加权局部估计量产生的彩色图像重叠在超声图像上。
16.如权利要求15所述的方法,其中所述彩色图像是选择落入与预定组织结构具有双单义关系的参考类别内的加权局部估计量来产生的。
17.一种用于对从接受超声检查的结构返回的射频超声波信号进行频谱分析的方法,该方法包括以下步骤:
a)向所要检查的所述结构的一部分传送超声波激励信号;
b)接收来自所述结构的输入射频响应信号;
c)为输入超声帧产生经过采样和数字化的帧;
d)将滤波系列应用于所述经过采样和数字化的帧,以便将射频响应信号的频带分解成多个频带;
e)产生频谱系数矩阵,该矩阵包含所述滤波操作所产生的系数或从中推导的系数;
f)为经过采样和数字化的帧的至少一些点确定相应的内插多项式PI,其中该多项式近似所述频谱系数在用以分解射频信号的各频带中的变化;
g)为所述点从内插多项式系数获取局部估计量,相互组合内插多项式的不同阶的至少两个系数。
18.如权利要求17所述的方法,其中用以细分所述射频信号的频带覆盖射频信号的整个频带。
19.如权利要求17所述的方法,其中分解后的多个频带是不同宽度和位置的频带。
20.如权利要求17所述的方法,其中所述滤波操作是使用时间-频率变换来获得的。
21.如权利要求20所述的方法,其中所述时间-频率变换是小波。
22.如权利要求20所述的方法,其中所述时间-频率变换是离散小波分组变换DWPT。
23.如权利要求17所述的方法,还包括以下阶段:确定局部估计量的统计分布,以及创建数值类别的集合,该集合能够以双单义方式来识别被采样的超声帧上的均匀部分。
24.如权利要求17所述的方法,其中使用所述局部估计量产生的彩色图像重叠在超声图像上。
25.如权利要求24所述的方法,其中所述彩色图像是选择落入与预定组织结构具有双单义关系的参考类别内的局部估计量来产生的。
26.一种超声设备,该设备包括超声探测器、用于采集和处理从接受超声检查的结构返回的射频信号的部件,其中:所述采集和处理部件经过编程来执行以下步骤:
a)将超声波激励信号传送到被检查的所述结构的一部分;
b)接收来自所述结构的射频响应信号;
c)应用一系列滤波操作,以便获得将射频响应信号的频带分解成多个频带;
d)从所述滤波操作产生的系数来计算局部估计量,其中所述局部估计量包含关于射频信号频谱的信息;
其中,所述局部估计量与表示所述局部估计量的统计分布形状的参数组合到超声图像的一部分中。
27.一种超声设备,该设备包括超声探测器、用于采集和处理从接受超声检查的结构返回的射频信号的部件,其中,所述采集和处理部件经过编程而执行以下步骤:
a)向所要检查的所述结构的一部分传送超声波激励信号;
b)接收来自所述结构的输入射频响应信号;
c)为输入超声帧产生经过采样和数字化的帧;
d)将滤波系列应用于所述经过采样和数字化的帧,以便获得将射频响应信号的频带分解成多个频带;
e)产生频谱系数矩阵,该矩阵包含所述滤波操作所产生的系数或从中推导的系数;
f)为经过采样和数字化的帧中的至少一些点确定相应的内插多项式PI,其中该多项式近似所述频谱系数在用以分解射频信号的频带中的变化;
g)为所述点从内插多项式系数获取局部估计量,相互组合内插多项式的不同阶的至少两个系数。
CN200480035699XA 2003-10-08 2004-10-05 用于超声波信号的局部频谱分析的改进方法和设备 Expired - Fee Related CN1954235B (zh)

Applications Claiming Priority (3)

Application Number Priority Date Filing Date Title
ITFI2003A000254 2003-10-08
IT000254A ITFI20030254A1 (it) 2003-10-08 2003-10-08 Metodo e dispositivo perfezionati per l'analisi spettrale
PCT/IT2004/000548 WO2005033738A1 (en) 2003-10-08 2004-10-05 Improved method and device for local spectral analysis of an ultrasonic signal

Publications (2)

Publication Number Publication Date
CN1954235A CN1954235A (zh) 2007-04-25
CN1954235B true CN1954235B (zh) 2010-05-05

Family

ID=34401303

Family Applications (1)

Application Number Title Priority Date Filing Date
CN200480035699XA Expired - Fee Related CN1954235B (zh) 2003-10-08 2004-10-05 用于超声波信号的局部频谱分析的改进方法和设备

Country Status (7)

Country Link
US (1) US7509861B2 (zh)
EP (1) EP1671156B1 (zh)
CN (1) CN1954235B (zh)
AT (1) ATE452348T1 (zh)
DE (1) DE602004024699D1 (zh)
IT (1) ITFI20030254A1 (zh)
WO (1) WO2005033738A1 (zh)

Families Citing this family (27)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
ITFI20020034A1 (it) * 2002-02-27 2003-08-27 Actis Active Sensors S R L Metodo e dispositivo per l'analisi spettrale di un segnale ecografico
DE102005040192A1 (de) * 2005-08-25 2007-03-01 Robert Bosch Gmbh Verfahren und Vorrichtung zur Erkennung von Quietschgeräuschen
KR20080064155A (ko) 2005-10-14 2008-07-08 어플라이드 리써치 어쏘시에이츠 뉴질랜드 리미티드 표면 특징을 모니터링하는 방법 및 장치
CN101356518B (zh) * 2005-11-04 2010-07-21 特克特朗尼克公司 用于多域标记的方法、系统和装置
ATE517584T1 (de) 2006-12-05 2011-08-15 Elesta S R L Ein satz optischer fasern zur perkutanen ablativen behandlung
US9549713B2 (en) 2008-04-24 2017-01-24 Boston Scientific Scimed, Inc. Methods, systems, and devices for tissue characterization and quantification using intravascular ultrasound signals
WO2009132188A1 (en) * 2008-04-24 2009-10-29 Boston Scientific Scimed, Inc. Methods, systems, and devices for tissue characterization by spectral similarity of intravascular ultrasound signals
CN101458331B (zh) * 2009-01-04 2011-09-07 中国人民解放军海军工程大学 一种多普勒声纳测试的声对接装置
JP4691732B1 (ja) * 2010-07-30 2011-06-01 国立大学法人 岡山大学 組織抽出システム
US9179844B2 (en) 2011-11-28 2015-11-10 Aranz Healthcare Limited Handheld skin measuring or monitoring device
CN102628941A (zh) * 2012-03-12 2012-08-08 陈明霞 一种宽频带超声波回波信号的重构方法
CN103020479A (zh) * 2012-12-28 2013-04-03 上海交通大学 基于非线性调频小波变换的信号瞬时频率估计方法
DE102013201975A1 (de) 2013-02-07 2014-08-07 Siemens Aktiengesellschaft Verfahren und Vorrichtung zur Verbesserung der SAFT-Analyse bei unregelmäßiger Messung
US10955273B2 (en) * 2013-04-04 2021-03-23 Texas Instruments Incorporated Extended range ADC flow meter
CN103735284B (zh) * 2013-12-26 2015-10-28 华南理工大学 基于线性扫描的三维超声弹性成像中rf信号估计方法
CN107567309A (zh) 2015-05-05 2018-01-09 波士顿科学国际有限公司 有设于超声成像系统换能器上的可膨胀材料的系统和方法
US10013527B2 (en) 2016-05-02 2018-07-03 Aranz Healthcare Limited Automatically assessing an anatomical surface feature and securely managing information related to the same
US11116407B2 (en) 2016-11-17 2021-09-14 Aranz Healthcare Limited Anatomical surface assessment methods, devices and systems
CN106771588A (zh) * 2016-12-15 2017-05-31 中国电子科技集团公司第四十研究所 一种时间窗口内时隙信号的频谱测试方法
JP2020508174A (ja) * 2017-02-27 2020-03-19 ラトガース,ザ ステート ユニバーシティ オブ ニュー ジャージー 肝臓および腎臓の癌の改善された診断のための計算超音波
US11903723B2 (en) 2017-04-04 2024-02-20 Aranz Healthcare Limited Anatomical surface assessment methods, devices and systems
US10670722B2 (en) * 2017-08-15 2020-06-02 Samsung Electronics Co., Ltd. Increase depth resolution and depth accuracy in ToF sensors by avoiding histogrammization
CN108535636A (zh) * 2018-05-16 2018-09-14 武汉大学 一种模拟电路基于参数随机分布邻近嵌入胜者为王的故障特征提取方法
CN112054798B (zh) * 2019-06-06 2023-09-19 中国科学院苏州纳米技术与纳米仿生研究所 信号采集电路及其采样频率调节方法、计算机存储介质
CN110333285B (zh) * 2019-07-04 2021-07-27 大连海洋大学 基于变分模态分解的超声兰姆波缺陷信号识别方法
CN112270282B (zh) * 2020-11-03 2021-12-10 华北电力大学 一种利用矩阵谱模的功率信号滤波方法和系统
CN113237951A (zh) * 2021-05-11 2021-08-10 重庆大学 一种基于形状上下文动态时间规整的金属板疲劳损伤超声导波检测方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1162250A (zh) * 1994-11-01 1997-10-15 舍林股份公司 超声处理以及执行这种处理的电路
US6066098A (en) * 1997-06-13 2000-05-23 Esaote S.P.A. Method for enhancing the diagnostic power of ultrasonic scanning systems by using real-time spectral maps, and device operating by said method
CN1440726A (zh) * 2002-10-01 2003-09-10 深圳迈瑞生物医疗电子股份有限公司 全数字超声频谱多普勒成像方法及装置

Family Cites Families (17)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
DE2702557C3 (de) 1977-01-22 1979-10-11 Labora Mannheim Gmbh Fuer Labortechnik, 6800 Mannheim Vorrichtung zum Bestimmen der Blutsenkungsgeschwindigkeit
FR2644915A1 (fr) * 1989-03-22 1990-09-28 Inst Nat Sante Rech Med Procede et dispositif d'analyse spectrale en temps reel de signaux instationnaires complexes
IT1233510B (it) 1989-04-05 1992-04-03 Diesse Diagnostica Apparecchio per la preparazione e la determinazione degli esami della velocita' di sedimentazione di liquidi organici ed altro
SE465140B (sv) 1989-12-13 1991-07-29 Tesi Ab Foerfarande och anordning foer att bestaemma blods saenkningsreaktion
US6506606B1 (en) 1995-06-06 2003-01-14 Brigham And Women's Hospital Method and apparatus for determining erythrocyte sedimentation rate and hematocrit
IT1286631B1 (it) 1996-05-16 1998-07-15 Diesse Diagnostica Apparecchio per la preparazione e la determinazione degli esami della velocita' di sedimentazione di liquidi organici ed altro
IT1286630B1 (it) 1996-05-16 1998-07-15 Diesse Diagnostica Una provetta per esami biologici di liquidi organici con apparecchi elettro-ottici
US5914272A (en) 1996-06-19 1999-06-22 Becton Dickinson And Company Test method for determining the erythrocyte sedimentation rate and a surfactant for use therein
AU4201597A (en) 1996-07-16 1998-02-09 Cdr S.R.L. Erythrosedimentation rate (esr) test apparatus
US6081612A (en) * 1997-02-28 2000-06-27 Electro Optical Sciences Inc. Systems and methods for the multispectral imaging and characterization of skin tissue
US5935074A (en) * 1997-10-06 1999-08-10 General Electric Company Method and apparatus for automatic tracing of Doppler time-velocity waveform envelope
GB9804560D0 (en) 1998-03-05 1998-04-29 Zynocyte Ltd Method of measuring erthrocyte sedimentation rate (esr),plasma viscosity and plasma fibrinogen of a blood sample
US6251077B1 (en) * 1999-08-13 2001-06-26 General Electric Company Method and apparatus for dynamic noise reduction for doppler audio output
US6523003B1 (en) * 2000-03-28 2003-02-18 Tellabs Operations, Inc. Spectrally interdependent gain adjustment techniques
DE10021286B4 (de) * 2000-05-02 2005-03-10 Kara Can Verfahren und Vorrichtung zur Kompression und/oder Dekompression von Daten
CA2463673C (en) 2001-10-19 2009-01-13 Monogen, Inc. Filtration system and method for obtaining a cytology layer
ITFI20020034A1 (it) 2002-02-27 2003-08-27 Actis Active Sensors S R L Metodo e dispositivo per l'analisi spettrale di un segnale ecografico

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1162250A (zh) * 1994-11-01 1997-10-15 舍林股份公司 超声处理以及执行这种处理的电路
US6066098A (en) * 1997-06-13 2000-05-23 Esaote S.P.A. Method for enhancing the diagnostic power of ultrasonic scanning systems by using real-time spectral maps, and device operating by said method
CN1440726A (zh) * 2002-10-01 2003-09-10 深圳迈瑞生物医疗电子股份有限公司 全数字超声频谱多普勒成像方法及装置

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
MASOTTI L ET AL.tissue differentiation based on radiofrequency echographicsignal local spectral content(rules:radiofrequency ultrasoniclocal estimator).2003 IEEE ULTRASONICS SYMPOSIUM PROCEEDINGS.HONOLULU,HAWAII2.2003,21030-1033. *
SUZUKI K ET AL.CEPSTRAL ANALYSIS OF ULTRASOUND IN CHRONICLIVER DISEASE-A PRELIMINARY STUDY IN THENON-INVASIVE EVALUATION OF STRUCTURAL CHANGE.FRONTIERS OF MEDICAL AND BIOLOGICAL ENGINEERING,VSP.ZEIST,NL3 4.1991,3(4),269-281.
SUZUKI K ET AL.CEPSTRAL ANALYSIS OF ULTRASOUND IN CHRONICLIVER DISEASE-A PRELIMINARY STUDY IN THENON-INVASIVE EVALUATION OF STRUCTURAL CHANGE.FRONTIERS OF MEDICAL AND BIOLOGICAL ENGINEERING,VSP.ZEIST,NL3 4.1991,3(4),269-281. *

Also Published As

Publication number Publication date
EP1671156A1 (en) 2006-06-21
ITFI20030254A1 (it) 2005-04-09
DE602004024699D1 (de) 2010-01-28
CN1954235A (zh) 2007-04-25
US20060277998A1 (en) 2006-12-14
WO2005033738A1 (en) 2005-04-14
US7509861B2 (en) 2009-03-31
EP1671156B1 (en) 2009-12-16
ATE452348T1 (de) 2010-01-15

Similar Documents

Publication Publication Date Title
CN1954235B (zh) 用于超声波信号的局部频谱分析的改进方法和设备
Cincotti et al. Frequency decomposition and compounding of ultrasound medical images with wavelet packets
Czerwinski et al. Line and boundary detection in speckle images
Andria et al. Linear filtering of 2-D wavelet coefficients for denoising ultrasound medical images
CN105793729B (zh) 使用多频率波形的超声图像形成和/或重建
EP1599122B1 (en) Non-invasive plaque characterization system
US8834376B2 (en) Dispersive ultrasound technology as a diagnostic device for traumatic brain injuries
Abbey et al. Observer efficiency in discrimination tasks simulating malignant and benign breast lesions imaged with ultrasound
Parker et al. Fine-tuning the H-scan for discriminating changes in tissue scatterers
US7815573B2 (en) Method and device for spectral analysis of an echographic signal
Jeong et al. Soft tissue differentiation using multiband signatures of high resolution ultrasonic transmission tomography
Perreault et al. Speckle simulation based on B-mode echographic image acquisition model
Granchi et al. Multidimensional spectral analysis of the ultrasonic radiofrequency signal for characterization of media
US20200225347A1 (en) Fine-tuning the h-scan for visualizing types of tissue scatterers
Kadah et al. Principal Component Analysis Based Hybrid Speckle Noise Reduction Technique for Medical Ultrasound Imaging
Kim et al. Multiband tissue differentiation in ultrasonic transmission tomography
AU2020103375A4 (en) Speckle Denoising System for Ultrasound Images with Framelet Transform and Gaussian Filter
EP1581118A1 (en) Detection of small-size defects in medical ultrasonic imaging
Zemp et al. Generalized NEQ for assessment of ultrasound image quality
Andria et al. Wavelet based methods for speckle reduction in ultrasound images
Saniie et al. Machine Learning and Modeling of Ultrasonic Signals for High-Fidelity Data Compression
Goudarzi Inverse Problem Formulation and Deep Learning Methods for Ultrasound Beamforming and Image Reconstruction
Bafaraj Speckle Reduction inMedical Ultrasound Imaging
Paridar et al. A Low-Complexity Compressive Sensing Algorithm for Point Target Detection Using Ultrafast Plane Wave Imaging
Al-Bayati Development of an algorithm for restoration and enhancement of ultrasound medical imaging

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C14 Grant of patent or utility model
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20100505

Termination date: 20141005

EXPY Termination of patent right or utility model