CN107358156B - 基于希尔伯特-黄变换的超声组织定征的特征提取方法 - Google Patents
基于希尔伯特-黄变换的超声组织定征的特征提取方法 Download PDFInfo
- Publication number
- CN107358156B CN107358156B CN201710417075.3A CN201710417075A CN107358156B CN 107358156 B CN107358156 B CN 107358156B CN 201710417075 A CN201710417075 A CN 201710417075A CN 107358156 B CN107358156 B CN 107358156B
- Authority
- CN
- China
- Prior art keywords
- imf
- hilbert
- frequency
- spectrum
- feature
- 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
Links
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/08—Feature extraction
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/25—Fusion techniques
- G06F18/253—Fusion techniques of extracted features
Landscapes
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Data Mining & Analysis (AREA)
- Artificial Intelligence (AREA)
- Physics & Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Signal Processing (AREA)
- Life Sciences & Earth Sciences (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Evolutionary Biology (AREA)
- Evolutionary Computation (AREA)
- Complex Calculations (AREA)
Abstract
本发明公开了一种基于希尔伯特‑黄变换的超声组织定征的特征提取方法,包括以下步骤:采集组织样本的超声回波RF信号多帧,选择感兴趣区,构造超声RF时间序列;利用HHT算法,对超声RF时间序列进行特征提取,得到基于时间、频率和能量关系的样本特征向量;利用rankfeatures算法,进行特征筛选与特征融合,计算样本特征融合指数,把该特征融合指数作为定征组织的特征。本发明所提出的特征融合指数能够有效的反映不同类别组织样本之间的差异,适用于超声组织定征的特征提取领域。
Description
技术领域
本发明涉及超声组织定征技术领域,特别涉及一种基于希尔伯特-黄变换的超声组织定征的特征提取方法。
背景技术
由于超声波进入生物组织后与其相互作用的机理尚未十分明了,人们只能通过提取超声回波信息并作出解释来间接达到识别组织结构、成分、状态的目的,从而促使人们进行超声组织定征特征提取的研究。
已有的超声组织定征特征提取方法主要有三类:基于超声B型图、基于背散射声束RF信号和基于超声RF时间序列分析。基于超声B型图的研究中,主要是利用B型图像灰度分布的纹理特,该方法使用了超声图像灰度,因此易受超声诊断仪的、TGC的调整等成像参数等影响;基于背散射声束RF信号的方法,需要进行深度衰减补偿,而不同组织衰减能力存在差异,同时个体差异也会使得声波传播路径不同,从而严重影响该类方法的效果。而基于超声RF时间序列的方法,RF时间序列来源于同一位置同一深度不同帧的RF数据,在一定程度上降低了个体差异的影响。
传统的频谱分析方法大都是基于快速傅里叶变换(FFT)和小波变换(WT)。由于存在频谱泄露和栅栏效应,FFT分析非整数次谐波存在较大误差,通过加窗插值算法虽可以较好的消除频谱泄露和栅栏效应,但这些算法往往以降低频率分辨率为代价。WT中,由于高通和低通分解滤波器并非是完全理想的滤波器,各个频带之间存在严重的交叉现象,因此很难实现各个频带的严格划分。以上方法本质均是一种基于基函数展开的理论,信号分析结果很大程度上依赖于设计者的经验。
与传统的时频分析方法相比,希尔伯特-黄变换(Hilbert-Huang Transform,简称HHT)作为一种相对新兴的时频分析方法,在处理非线性非平稳信号中有着较为明显的优势,HHT可以根据信号本身的特点自适应的分解出一系列不同频率的分量,即固有模态函数(IMF),不用考虑基函数或者窗函数的选择等其他非自身因素的干扰,得到的结果能更好的反映信号本身的特点。
发明内容
本发明的目的在于克服现有技术的缺点与不足,提供一种基于希尔伯特-黄变换的超声组织定征的特征提取方法,所提出的特征融合指数能够有效的反映不同类别组织样本之间的差异,适用于超声组织定征的特征提取领域。
本发明的目的通过以下的技术方案实现:
一种基于希尔伯特-黄变换的可用于超声组织定征的特征提取方法,包括以下步骤:
S1:采集组织样本的超声回波RF信号多帧,选择ROI,构造超声RF时间序列;
S2:利用HHT算法,对超声RF时间序列进行特征提取,得到基于时间、频率和能量关系的样本特征向量;
S3:利用rankfeatures算法,进行特征筛选与特征融合,计算样本特征融合指数,把该特征融合指数作为定征组织的特征。
优选的,所述步骤S1中,构造超声RF时间序列的具体过程为:
S1-1:扫描组织,获取其超声回波RF信号多帧;
S1-2:解调某帧超声回波RF信号并显示B型图;
S1-3:在B型图上选取大小为M×N的ROI;
S1-4:对ROI内的每点取其前L帧RF信号构造出M×N个长度为L的超声RF时间序列。
优选的,所述步骤S2中,利用HHT算法,计算样本特征向量的方法为:
S2-1:取样本ROI内的一条时间序列x(t),进行集合经验模态分解,得到n个IMF分量;
S2-2:基于IMF分量提取时域特征:
IMF过零点数5个:IMF1-ZCs、IMF2-ZCs、IMF3-ZCs、IMF4-ZCs、IMF5-ZCs;
IMF方差5个:IMF1-Var、IMF2-Var、IMF3-Var、IMF4-Var、IMF5-Var;
IMF方差贡献率5个IMF1-VarR、IMF2-VarR、IMF3-VarR、IMF4-VarR、IMF5-VarR;
S2-3:对时间序列x(t)的n阶IMF,分别进行希尔伯特谱分析;
具体的,步骤S2-3包括:
S2-4:基于幅值函数和能量函数提取时域-能量特征:
5个IMF平均强度:IMF1-AvgA、IMF2-AvgA、IMF3-AvgA、IMF4-AvgA、IMF5-AvgA,
5个IMF能量:IMF1-Egy、IMF2-Egy、IMF3-Egy、IMF4-Egy、IMF5-Egy;
S2-5:基于瞬时频率提取时域-频率特征:
5个IMF最高频率:IMF1-MaxF、IMF2-MaxF、IMF3-MaxF、IMF4-MaxF、IMF5-MaxF;
S2-6:基于Hilbert边际谱提取能量-频率特征:
5个IMF平均中心频率:IMF1-MCF、IMF2-MCF、IMF3-MCF、IMF4-MCF、IMF5-MCF;1个x(t)平均中心频率Orig-MCF;
S2-7:基于Hilbert边际谱提取频域-能量特征:
x(t)的Hilbert边际谱熵Orig-EgyS;
x(t)的归一化Hilbert边际谱低频段能量、中低频段能量、中高频段能量、高频段能量Orig-MargS1、Orig-MargS2、Orig-MargS3、Orig-MargS4;
S2-8:基于Hilbert边际谱提取频域-能量曲线拟合特征:
x(t)的归一化Hilbert边际谱直线拟合斜率、截距:O-MLFSlope、O-MLFInterp;
一次多项式指数函数曲线拟合斜率、截距:O-MEFSlope、O-MEFInterp;
7个六阶曲线拟合系数:O-MSOFa0、O-MSOFa1、O-MSOFa2、O-MSOFa3、O-MSOFa4、O-MSOFa5、O-MSOFa6;
S2-9:基于Hilbert谱提取时域-频域-能量特征:
x(t)的Hilbert谱灰度统计直方图的统计特征——均值TFImgMean、方差TFImgSD、偏度TFImgSkew、峰度TFImgKurto;
S2-10:分别取样本ROI内的M×N条时间序列xi(t),其中i=0,1,2,...,M×N-1,重复执行步骤S2-1~S2-9,得到M×N条时间序列对应的M×N个特征向量,最后,再将这M×N个特征向量求平均得到得到一个特征向量,即样本特征向量。
优选的,步骤S2-1中得到n个IMF分量,n≥5。
具体的,步骤S2-1具体的过程包括:
步骤a:向x(t)中加入高斯白噪声序列,得到X(t)=x(t)+ω(t);
步骤c:如果h1(t)满足本征模函数IMF条件,那么h1(t)就是X(t)的第一IMF分量;否则重复步骤b,把h1(t)作为新的X(t),计算m11(t),再判断h11(t)=h1(t)-m11(t)是否满足IMF条件,如不满足,循环k次,直到h1k(t)=h1(k-1)(t)-m1k(t)满足IMF条件。记c1(t)=h1k(t)为X(t)的一阶IMF分量,此时,剩余分量r1(t)=X(t)-c1(t);
步骤d:将r1(t)作为新的X(t),重复步骤a~c,得到X(t)的第二个满足IMF条件的分量c2(t),循环执行步骤a~c,直到rn(t)成为一个单调函数,循环结束,从而得到信号X(t)的n(n≥5)阶IMF分量的集合;
步骤e:步骤a~d重复执行N次,每次添加不同的高斯白噪声序列ωi(t),分解得到N组IMF分量的集合,再对N组IMF分量中,对应阶数的分量求平均作为最终的IMF分量,ck(t)表示对x(t)经过EEMD分解后,得到的k阶IMF分量。
具体的,步骤S2-2包括:
步骤a:对于长度为L的IMF分量ci(t),统计其过零点数目的方法为:当1≤j≤L时,若c(j)*c(j-1)<0,则记为一个过零点,遍历[1,L]找出IMF分量ci(t)的所有过零点,即得到ci(t)的过零点数目,于是,可以得到前5阶IMF的过零点数IMF1-ZCs、IMF2-ZCs、IMF3-ZCs、IMF4-ZCs、IMF5-ZCs;
步骤b:按照公式计算前5阶IMF的方差IMF1-Var、IMF2-Var、IMF3-Var、IMF4-Var、IMF5-Var,再按照公式计算前5阶IMF的方差贡献率IMF1-VarR、IMF2-VarR、IMF3-VarR、IMF4-VarR、IMF5-VarR。
优选的,步骤S2-4中,对于时间序列x(t)的n阶IMF经过步骤S2-3的步骤a处理以后,可以得到n个长度为L的幅值函数ai(t),能量函数ei(t),(i=0,1,2,...,n-1),根据公式计算前5阶IMF的平均强度,再根据公式计算前5阶IMF的能量。
优选的,步骤S2-5中,对于时间序列x(t)的n阶IMF经过步骤S2-3的步骤a处理以后,可以得到n个长度为L的瞬时频率函数fi(t),(i=0,1,2,...,n-1),根据公式计算前5阶IMF的最高频率。
优选的,步骤S2-6中,对于时间序列x(t)的n阶IMF经过步骤S2-3的步骤a~c处理以后,可以得到n个包含m频率点数的Hilbert边际谱hi(f),以及x(t)的包含m个频率点数的Hilbert边际谱h(f),分别根据公式和其中hi(fj)表示第i个IMF分量的Hilbert边际谱的第j个频率点的幅值,fj表示第j个频率点,计算前5阶IMF的平均中心频率和x(t)的平均中心频率。
优选的,步骤S2-7中包括步骤:
优选的,步骤S2-8中采用最小二乘法对x(t)的归一化Hilbert边际谱ho(f)进行直线拟合得到O-MLFSlope和O-MLFInterp,进行一次多项式指数函数曲线拟合得到O-MEFSlope和O-MEFInterp,进行六阶曲线拟合得到O-MSOFa0、O-MSOFa1、O-MSOFa2、O-MSOFa3、O-MSOFa4、O-MSOFa5、O-MSOFa6。
优选的,步骤S2-9具体包括:
步骤a:对于时间序列x(t)的n阶IMF经过步骤S2-3的步骤a~c处理以后,可以得到x(t)的Hilbert幅值谱H(f,t),以时间t为横轴,频率f为纵轴,利用幅度控制灰度,得到时间序列x(t)的Hilbert谱时频灰度图像;
步骤b:再由x(t)的Hilbert谱时频灰度图像,以图像的灰度级作为横轴,图像上每一个灰度级出现的次数作为纵轴,得到时间序列x(t)的Hilbert谱时频灰度图像直方图hi,其中i=0,1,2,...,l-1。hi表示灰度直方图第i个灰度级像素出现的次数,l表示时频灰度图像的灰度级数;
优选的,步骤S3中,rankfeatures算法采用MATLAB软件的工具箱函数,该算法按照类别分离准则对关键特征进行排序,具体利用到函数
[IDX,Z]=rankfeatures(X,Group),
使用二元分类的独立性评价准则对特征进行排序,输入X是一个矩阵,每一列表示一个观测值向量而行数对应于原始的特征数;输入Group包含类标签,输出IDX是一个由X的行下标构成的列表,下标列表按照矩阵X中特征的显著性从高到低排列,输出Z是使用准则后,各个特征对应的权重。
优选的,步骤S3中利用rankfeatures算法计算样本特征融合指数,具体的步骤为:
S3-1:将所有样本特征向量构成的特征数据矩阵赋值给X,将所有样本标签构成的向量赋值给Group,输入到[IDX,Z]=rankfeatures(X,Group)函数进行处理,得到包含56个分量的特征权重向量Z;
S3-2:选择权重较大的m个特征进行加权求和,得到样本特征融合指数。
本发明与现有技术相比,具有如下优点和有益效果:
1、本发明采用希尔伯特-黄变换来处理超声RF时间序列,从时间、频率和能量关系的角度提取特征,发明所采取的HHT是更适合于非平稳、非线性信号分析的时频分析方法,提取的特征信息量更丰富。
2、本发明利用rankfeatures算法进行特征筛选,计算各个特征的权重,选择权重较大的特征进行加权求和,构造特征融合指数,本发明所提出的特征融合指数能够反映不同类别的组织样本之间的差异。
3、本发明基于超声RF时间序列,由于超声RF时间序列的信号取自同一位置,减小了超声波传播路径差异带来的误差,且不需要深度衰减补偿,从而可提高定征精度;对组织是否均匀并无特殊要求,应用范围更加广泛。
4、本发明所基于的超声RF时间序列可在常规B超检查中获取,不会产生额外的设备硬件开支。
附图说明
图1是实施例2方法的流程图;
图2是实施例2对照组和治疗组乳腺肿瘤区域病理图:图2(a)对照组;图2(b)治疗组;
图3是实施例2对照组和治疗组乳腺肿瘤区域B型图:图3(a)对照组;图3(b)治疗组;
图4是实施例2大鼠乳腺RF时间序列分布图。
具体实施方式
下面结合实施例及附图对本发明作进一步详细的描述,但本发明的实施方式不限于此。
实施例1
一种基于希尔伯特-黄变换的超声组织定征的特征提取方法,包括以下步骤:
S1:采集组织样本的超声回波RF信号多帧,选择感兴趣区域(Region OfInterest,简称ROI),构造超声RF时间序列;
S2:利用HHT算法,对超声RF时间序列进行特征提取,得到基于时间、频率和能量关系的样本特征向量;
S3:利用rankfeatures算法,进行特征筛选与特征融合,计算样本特征融合指数,把该特征融合指数作为定征组织的特征。
所述步骤S1中,构造超声RF时间序列的具体过程为:
S1-1:扫描组织,获取其超声回波RF信号多帧;
S1-2:解调某帧超声回波RF信号并显示B型图;
S1-3:在B型图上选取大小为M×N的ROI;
S1-4:对ROI内的每点取其前L帧RF信号构造出M×N个长度为L的超声RF时间序列。
所述步骤S2中,利用HHT算法,计算样本特征向量的方法为:
S2-1:取样本ROI内的一条时间序列x(t),进行集合经验模态分解(EnsembleEmpirical Mode Decomposition,简称EEMD),得到n(n≥5)个IMF分量,具体的过程为:
步骤a:向x(t)中加入高斯白噪声序列,得到X(t)=x(t)+ω(t);
步骤c:如果h1(t)满足本征模函数(IntrinsicModeFunction,简称IMF)条件,那么h1(t)就是X(t)的第一IMF分量;否则重复步骤b,把h1(t)作为新的X(t),计算m11(t),再判断h11(t)=h1(t)-m11(t)是否满足IMF条件,如不满足,循环k次,直到h1k(t)=h1(k-1)(t)-m1k(t)满足IMF条件。记c1(t)=h1k(t)为X(t)的一阶IMF分量,此时,剩余分量r1(t)=X(t)-c1(t);
步骤d:将r1(t)作为新的X(t),重复步骤a~c,得到X(t)的第二个满足IMF条件的分量c2(t),循环执行步骤a~c,直到rn(t)成为一个单调函数,循环结束,从而得到信号X(t)的n(n≥5)阶IMF分量的集合;
步骤e:步骤a~d重复执行N次,每次添加不同的高斯白噪声序列ωi(t),分解得到N组IMF分量的集合,再对N组IMF分量中,对应阶数的分量求平均作为最终的IMF分量,ck(t)表示对x(t)经过EEMD分解后,得到的k阶IMF分量。
S2-2:提取时域特征(基于IMF分量):IMF过零点数(5个)IMF1-ZCs、IMF2-ZCs、IMF3-ZCs、IMF4-ZCs、IMF5-ZCs,IMF方差(5个)IMF1-Var、IMF2-Var、IMF3-Var、IMF4-Var、IMF5-Var,IMF方差贡献率(5个)IMF1-VarR、IMF2-VarR、IMF3-VarR、IMF4-VarR、IMF5-VarR,具体过程为:
步骤a:对于长度为L的IMF分量ci(t),统计其过零点数目的方法为:当1≤j≤L时,若c(j)*c(j-1)<0,则记为一个过零点,遍历[1,L]找出IMF分量ci(t)的所有过零点,即得到ci(t)的过零点数目,于是,可以得到前5阶IMF的过零点数IMF1-ZCs、IMF2-ZCs、IMF3-ZCs、IMF4-ZCs、IMF5-ZCs;
步骤b:按照公式计算前5阶IMF的方差IMF1-Var、IMF2-Var、IMF3-Var、IMF4-Var、IMF5-Var,再按照公式计算前5阶IMF的方差贡献率IMF1-VarR、IMF2-VarR、IMF3-VarR、IMF4-VarR、IMF5-VarR;
S2-3:对时间序列x(t)的n阶IMF,分别进行希尔伯特谱分析(HilbertSpectrumAnalysis,简称HSA),具体过程如下:
步骤a:对一个IMF分量c(t)进行Hilbert变换式中P为柯西主值,然后求c(t)的解析信号表示成极坐标的形式z(t)=a(t)eiθ(t),其中幅值函数能量函数相位函数由相位函数可以计算出c(t)的瞬时频率
S2-4:提取时域-能量特征(基于幅值函数和能量函数):IMF平均强度(5个)IMF1-AvgA、IMF2-AvgA、IMF3-AvgA、IMF4-AvgA、IMF5-AvgA,IMF能量(5个)IMF1-Egy、IMF2-Egy、IMF3-Egy、IMF4-Egy、IMF5-Egy,具体方法为:
对于时间序列x(t)的n阶IMF经过步骤S2-3的步骤a处理以后,可以得到n个长度为L的幅值函数ai(t),能量函数ei(t),(i=0,1,2,...,n-1),根据公式计算前5阶IMF的平均强度,再根据公式计算前5阶IMF的能量。
S2-5:提取时域-频率特征(基于瞬时频率):IMF最高频率(5个)IMF1-MaxF、IMF2-MaxF、IMF3-MaxF、IMF4-MaxF、IMF5-MaxF,具体方法为:
S2-6:提取能量-频率特征(基于Hilbert边际谱):IMF平均中心频率(5个)IMF1-MCF、IMF2-MCF、IMF3-MCF、IMF4-MCF、IMF5-MCF,x(t)平均中心频率(1个)Orig-MCF,具体方法为:
对于时间序列x(t)的n阶IMF经过步骤S2-3的步骤a~c处理以后,可以得到n个包含m频率点数的Hilbert边际谱hi(f),以及x(t)的包含m个频率点数的Hilbert边际谱h(f),分别根据公式和其中hi(fj)表示第i个IMF分量的Hilbert边际谱的第j个频率点的幅值,fj表示第j个频率点,计算前5阶IMF的平均中心频率和x(t)的平均中心频率。
S2-7:提取频域-能量特征(基于Hilbert边际谱):x(t)的Hilbert边际谱熵(1个)Orig-EgyS,x(t)的归一化Hilbert边际谱低频段能量、中低频段能量、中高频段能量、高频段能量(4个)Orig-MargS1、Orig-MargS2、Orig-MargS3、Orig-MargS4,具体步骤如下:
S2-8:提取频域-能量曲线拟合特征(基于Hilbert边际谱):x(t)的归一化Hilbert边际谱直线拟合斜率、截距(2个)O-MLFSlope、O-MLFInterp,一次多项式指数函数曲线拟合斜率、截距(2个)O-MEFSlope、O-MEFInterp,六阶曲线拟合系数(7个)O-MSOFa0、O-MSOFa1、O-MSOFa2、O-MSOFa3、O-MSOFa4、O-MSOFa5、O-MSOFa6,具体方法为:
采用最小二乘法对x(t)的归一化Hilbert边际谱ho(f)进行直线拟合得到O-MLFSlope和O-MLFInterp,进行一次多项式指数函数曲线拟合得到O-MEFSlope和O-MEFInterp,进行六阶曲线拟合得到O-MSOFa0、O-MSOFa1、O-MSOFa2、O-MSOFa3、O-MSOFa4、O-MSOFa5、O-MSOFa6。
S2-9:提取时域-频域-能量特征(基于Hilbert谱):x(t)的Hilbert谱灰度统计直方图的统计特征——均值、方差、偏度、峰度(4个)TFImgMean、TFImgSD、TFImgSkew、TFImgKurto,具体步骤如下:
步骤a:对于时间序列x(t)的n阶IMF经过步骤S2-3的步骤a~c处理以后,可以得到x(t)的Hilbert幅值谱H(f,t),以时间t为横轴,频率f为纵轴,利用幅度控制灰度,得到时间序列x(t)的Hilbert谱时频灰度图像;
步骤b:再由x(t)的Hilbert谱时频灰度图像,以图像的灰度级作为横轴,图像上每一个灰度级出现的次数作为纵轴,得到时间序列x(t)的Hilbert谱时频灰度图像直方图hi,其中i=0,1,2,...,l-1。hi表示灰度直方图第i个灰度级像素出现的次数,l表示时频灰度图像的灰度级数;
S2-10:分别取样本ROI内的M×N条时间序列xi(t),其中i=0,1,2,...,M×N-1,重复执行步骤S2-1~S2-9,得到M×N条时间序列对应的M×N个特征向量,最后,再将这M×N个特征向量求平均得到得到一个特征向量,即样本特征向量。
所述步骤S3中,rankfeatures算法采用MATLAB软件的工具箱函数,该算法按照类别分离准则对关键特征进行排序,具体利用到函数
[IDX,Z]=rankfeatures(X,Group),
使用二元分类的独立性评价准则对特征进行排序,输入X是一个矩阵,每一列表示一个观测值向量而行数对应于原始的特征数;输入Group包含类标签,输出IDX是一个由X的行下标构成的列表,下标列表按照矩阵X中特征的显著性从高到低排列,输出Z是使用准则后,各个特征对应的权重。默认采用联合方差估计双样本t检验准则。
利用rankfeatures算法计算样本特征融合指数,具体的步骤为:
S3-1:将所有样本特征向量构成的特征数据矩阵赋值给X,将所有样本标签构成的向量赋值给Group,输入到[IDX,Z]=rankfeatures(X,Group)函数进行处理,得到包含56个分量的特征权重向量Z;
S3-2:选择权重较大的m(m<56)个特征进行加权求和,得到样本特征融合指数。
实施例2
本实施例除下述特征外其他结构同实施例1:
如图1所示,本实施例是一种基于希尔伯特-黄变换的超声组织定征的特征提取方法应用于乳腺癌组织化疗前后组织微结构变化的检测,包括如下步骤:
步骤一:构造超声RF时间序列
使用加拿大Ultrasonix公司生产的Sonix TOUCH及中心频率6.6MHz的宽频线阵超声探头扫描雌性BALB/C裸鼠的胸腔皮下组织,记录多帧超声回波RF信号。本实验采用27只4到6周龄的雌性BALB/C裸鼠作为实验对象。首先对人乳腺癌细胞MCF-7进行培养,再加入10%胎牛血清和100单位/ml请链霉素,置于温度为37℃、CO2浓度为5%的环境中,然后在无菌条件下将大约4*108个MCF-7细胞接种到裸鼠的胸腔皮下。在医学角度,这些裸鼠模型可当做无差别的样本进行研究分析。待肿瘤生长至直径约5mm时将乳腺癌动物模型分为治疗组和对照组,数量分别为14只和13只,治疗组裸鼠每天以3mg/kg的紫杉醇进行化疗,而对照组只给予生理盐水。图2是对照组和治疗组乳腺肿瘤区域病理图,可以观察到治疗组肿瘤细胞密度显著减少、细胞之间的缝隙增大等组织微结构的改变;
S101:扫描实验裸鼠乳腺肿瘤组织,获取其超声回波RF信号;
S102:选取第100帧超声回波RF信号解调显示B型图,图3是对照组和治疗组乳腺肿瘤区域B型图,对比两幅B型图,肿瘤区域的图像区别并不明显,无法识别组织微结构的改变;
S103:在B型图上选择20×70个点大小的ROI;
S104:对ROI内的每个点取其前256帧回波信号,即得到20×70个点的长度为256的超声RF时间序列,图4是超声RF时间序列分布图;
步骤二:利用HHT算法,计算实验裸鼠样本特征向量
按照计算样本特征向量的方法,对每个实验裸鼠样本ROI的时间序列进行特征提取,可以得到包含56个分量的样本特征向量;
步骤三:利用rankfeatures算法,计算裸鼠样本特征融合指数
S301:对27个裸鼠分别计算样本特征向量,并将27×56的特征数据矩阵的转置矩阵赋值给X,将27个样本标签构成的向量赋值给Group,利用公式[IDX,Z]=rankfeatures(X,Group),计算得到56个特征的权重向量;
S302:选择前16个权重较大的特征进行加权求和,经过计算,本实施例中,对照组14个样本的特征融合指数分别为:0.8334,1.2672,0.9598,16.4901,9.1725,1.0033,3.0090,1.2001,0.8029,1.6267,1.2665,2.0368,1.3963,1.3334;治疗组13个样本的特征融合指数分别为:1.6160,1.5381,0.7148,0.7618,0.7112,1.0661,0.8636,1.5183,0.9405,1.0340,0.8058,0.8384,0.6947。
对照组的特征融合指数的均值为3.0284,方差为19.6270;治疗组的特征融合指数的均值为1.0079,方差为0.1119。
由此可以看出,对照组特征融合指数的均值大于治疗组特征融合指数的均值,对照组特征融合指数的方差大于治疗组特征融合指数的方差,这一结果说明经过化疗的治疗组样本的特征融合指数小于对照组样本,并且波动较小。说明了特征融合指数可用于区分微结构发生改变的组织,因而可以作为组织定征的特征。
上述实施例为本发明较佳的实施方式,但本发明的实施方式并不受上述实施例的限制,其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化,均应为等效的置换方式,都包含在本发明的保护范围之内。
Claims (8)
1.基于希尔伯特-黄变换的超声组织定征的特征提取方法,其特征在于,包括以下步骤:
S1:采集组织样本的超声回波RF信号多帧,选择ROI,构造超声RF时间序列;所述构造超声RF时间序列的具体过程为:
S1-1:扫描组织,获取其超声回波RF信号多帧;
S1-2:解调某帧超声回波RF信号并显示B型图;
S1-3:在B型图上选取大小为M×N的ROI;
S1-4:对ROI内的每点取其前L帧RF信号构造出M×N个长度为L的超声RF时间序列;
S2:利用HHT算法,对超声RF时间序列进行特征提取,得到基于时间、频率和能量关系的样本特征向量;其中计算样本特征向量的具体过程为:
S2-1:取样本ROI内的一条时间序列x(t),进行集合经验模态分解,得到n个IMF分量;
S2-2:基于IMF分量提取时域特征:
IMF过零点数5个:IMF1-ZCs、IMF2-ZCs、IMF3-ZCs、IMF4-ZCs、IMF5-ZCs;
IMF方差5个:IMF1-Var、IMF2-Var、IMF3-Var、IMF4-Var、IMF5-Var;
IMF方差贡献率5个IMF1-VarR、IMF2-VarR、IMF3-VarR、IMF4-VarR、IMF5-VarR;
S2-3:对时间序列x(t)的n阶IMF,分别进行希尔伯特谱分析;
具体的,步骤S2-3包括:
步骤a:对一个IMF分量c(t)进行Hilbert变换式中P为柯西主值,然后求c(t)的解析信号表示成极坐标的形式z(t)=a(t)eiθ(t),其中幅值函数能量函数相位函数由相位函数可以计算出c(t)的瞬时频率
S2-4:基于幅值函数和能量函数提取时域-能量特征:
5个IMF平均强度:IMF1-AvgA、IMF2-AvgA、IMF3-AvgA、IMF4-AvgA、IMF5-AvgA;
5个IMF能量:IMF1-Egy、IMF2-Egy、IMF3-Egy、IMF4-Egy、IMF5-Egy;
S2-5:基于瞬时频率提取时域-频率特征:
5个IMF最高频率:IMF1-MaxF、IMF2-MaxF、IMF3-MaxF、IMF4-MaxF、IMF5-MaxF;
S2-6:基于Hilbert边际谱提取能量-频率特征:
5个IMF平均中心频率:IMF1-MCF、IMF2-MCF、IMF3-MCF、IMF4-MCF、IMF5-MCF;1个x(t)平均中心频率Orig-MCF;
S2-7:基于Hilbert边际谱提取频域-能量特征:
x(t)的Hilbert边际谱熵Orig-EgyS;
x(t)的归一化Hilbert边际谱低频段能量、中低频段能量、中高频段能量、高频段能量Orig-MargS1、Orig-MargS2、Orig-MargS3、Orig-MargS4;
S2-8:基于Hilbert边际谱提取频域-能量曲线拟合特征:
x(t)的归一化Hilbert边际谱直线拟合斜率、截距:O-MLFSlope、O-MLFInterp;
一次多项式指数函数曲线拟合斜率、截距:O-MEFSlope、O-MEFInterp;
7个六阶曲线拟合系数:O-MSOFa0、O-MSOFa1、O-MSOFa2、O-MSOFa3、O-MSOFa4、O-MSOFa5、O-MSOFa6;
S2-9:基于Hilbert谱提取时域-频域-能量特征:
x(t)的Hilbert谱灰度统计直方图的统计特征——均值TFImgMean、方差TFImgSD、偏度TFImgSkew、峰度TFImgKurto;
S2-10:分别取样本ROI内的M×N条时间序列xi(t),其中i=0,1,2,...,M×N-1,重复执行步骤S2-1~S2-9,得到M×N条时间序列对应的M×N个特征向量,最后,再将这M×N个特征向量求平均得到得到一个特征向量,即样本特征向量;
S3:利用rankfeatures算法,进行特征筛选与特征融合,计算样本特征融合指数,把该特征融合指数作为定征组织的特征。
2.根据权利要求1所述的特征提取方法,其特征在于,步骤S2-1中得到n个IMF分量,n≥5。
3.根据权利要求1所述的特征提取方法,其特征在于,步骤S2-1具体的过程包括:
步骤a:向x(t)中加入高斯白噪声序列,得到X(t)=x(t)+ω(t);
步骤c:如果h1(t)满足本征模函数IMF条件,那么h1(t)就是X(t)的第一IMF分量;否则重复步骤b,把h1(t)作为新的X(t),计算m11(t),再判断h11(t)=h1(t)-m11(t)是否满足IMF条件,如不满足,循环k次,直到h1k(t)=h1(k-1)(t)-m1k(t)满足IMF条件;记c1(t)=h1k(t)为X(t)的一阶IMF分量,此时,剩余分量r1(t)=X(t)-c1(t);
步骤d:将r1(t)作为新的X(t),重复步骤a~c,得到X(t)的第二个满足IMF条件的分量c2(t),循环执行步骤a~c,直到rn(t)成为一个单调函数,循环结束,从而得到信号X(t)的n(n≥5)阶IMF分量的集合;
4.根据权利要求1所述的特征提取方法,其特征在于,步骤S2-2包括:
步骤a:对于长度为L的IMF分量ci(t),统计其过零点数目的方法为:当1≤j≤L时,若c(j)*c(j-1)<0,则记为一个过零点,遍历[1,L]找出IMF分量ci(t)的所有过零点,即得到ci(t)的过零点数目,于是,可以得到前5阶IMF的过零点数IMF1-ZCs、IMF2-ZCs、IMF3-ZCs、IMF4-ZCs、IMF5-ZCs;
5.根据权利要求1所述的特征提取方法,其特征在于,步骤S2-4到步骤S2-8具体有:
步骤S2-4中,对于时间序列x(t)的n阶IMF经过步骤S2-3的步骤a处理以后,可以得到n个长度为L的幅值函数ai(t),能量函数ei(t),(i=0,1,2,...,n-1),根据公式计算前5阶IMF的平均强度,再根据公式计算前5阶IMF的能量;
步骤S2-5中,对于时间序列x(t)的n阶IMF经过步骤S2-3的步骤a处理以后,可以得到n个长度为L的瞬时频率函数fi(t),(i=0,1,2,...,n-1),根据公式计算前5阶IMF的最高频率;
步骤S2-6中,对于时间序列x(t)的n阶IMF经过步骤S2-3的步骤a~c处理以后,可以得到n个包含m频率点数的Hilbert边际谱hi(f),以及x(t)的包含m个频率点数的Hilbert边际谱h(f),分别根据公式和其中hi(fj)表示第i个IMF分量的Hilbert边际谱的第j个频率点的幅值,fj表示第j个频率点,计算前5阶IMF的平均中心频率和x(t)的平均中心频率;
步骤S2-7中包括步骤:
步骤S2-8中采用最小二乘法对x(t)的归一化Hilbert边际谱ho(f)进行直线拟合得到O-MLFSlope和O-MLFInterp,进行一次多项式指数函数曲线拟合得到O-MEFSlope和O-MEFInterp,进行六阶曲线拟合得到O-MSOFa0、O-MSOFa1、O-MSOFa2、O-MSOFa3、O-MSOFa4、O-MSOFa5、O-MSOFa6。
6.根据权利要求1所述的特征提取方法,其特征在于,步骤S2-9具体包括:
步骤a:对于时间序列x(t)的n阶IMF经过步骤S2-3的步骤a~c处理以后,可以得到x(t)的Hilbert幅值谱H(f,t),以时间t为横轴,频率f为纵轴,利用幅度控制灰度,得到时间序列x(t)的Hilbert谱时频灰度图像;
步骤b:再由x(t)的Hilbert谱时频灰度图像,以图像的灰度级作为横轴,图像上每一个灰度级出现的次数作为纵轴,得到时间序列x(t)的Hilbert谱时频灰度图像直方图hi,其中i=0,1,2,...,l-1,hi表示灰度直方图第i个灰度级像素出现的次数,l表示时频灰度图像的灰度级数;
7.根据权利要求1所述的特征提取方法,其特征在于,步骤S3中,rankfeatures算法采用MATLAB软件的工具箱函数,该算法按照类别分离准则对关键特征进行排序,具体利用到函数
[IDX,Z]=rankfeatures(X,Group),
使用二元分类的独立性评价准则对特征进行排序,输入X是一个矩阵,每一列表示一个观测值向量而行数对应于原始的特征数;输入Group包含类标签,输出IDX是一个由X的行下标构成的列表,下标列表按照矩阵X中特征的显著性从高到低排列,输出Z是使用准则后,各个特征对应的权重。
8.根据权利要求1所述的特征提取方法,其特征在于,步骤S3中利用rankfeatures算法计算样本特征融合指数,具体的步骤为:
S3-1:将所有样本特征向量构成的特征数据矩阵赋值给X,将所有样本标签构成的向量赋值给Group,输入到[IDX,Z]=rankfeatures(X,Group)函数进行处理,得到包含多个分量的特征权重向量Z;
S3-2:从得到的包含多个分量的特征权重向量Z中,选择权重大的前m个特征进行加权求和,得到样本特征融合指数。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201710417075.3A CN107358156B (zh) | 2017-06-06 | 2017-06-06 | 基于希尔伯特-黄变换的超声组织定征的特征提取方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201710417075.3A CN107358156B (zh) | 2017-06-06 | 2017-06-06 | 基于希尔伯特-黄变换的超声组织定征的特征提取方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN107358156A CN107358156A (zh) | 2017-11-17 |
CN107358156B true CN107358156B (zh) | 2020-05-19 |
Family
ID=60271035
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201710417075.3A Expired - Fee Related CN107358156B (zh) | 2017-06-06 | 2017-06-06 | 基于希尔伯特-黄变换的超声组织定征的特征提取方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN107358156B (zh) |
Families Citing this family (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN108614926A (zh) * | 2018-04-12 | 2018-10-02 | 西安交通大学 | 一种基于流形学习与希尔伯特-黄变换相结合的结构模态参数辨识方法 |
CN109214092A (zh) * | 2018-09-11 | 2019-01-15 | 吉林大学 | 基于希尔伯特黄变换的四分量钻孔应变数据异常提取方法 |
CN110032585B (zh) * | 2019-04-02 | 2021-11-30 | 北京科技大学 | 一种时间序列双层符号化方法及装置 |
CN113205000A (zh) * | 2021-04-07 | 2021-08-03 | 江苏师范大学 | 一种基于柯西Score极坐标图调制方式识别方法 |
CN113191232B (zh) * | 2021-04-21 | 2023-04-18 | 西安交通大学 | 基于多模态同源特征和XGBoost模型的电静液作动器故障识别方法 |
CN114509604B (zh) * | 2022-04-18 | 2022-09-02 | 国网江西省电力有限公司电力科学研究院 | 一种gis壳体暂态地电位升波形频谱分析方法及*** |
CN116028882B (zh) * | 2023-03-29 | 2023-06-02 | 深圳市傲天科技股份有限公司 | 用户标注和分类方法、装置、设备及存储介质 |
Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102855623A (zh) * | 2012-07-19 | 2013-01-02 | 哈尔滨工业大学 | 基于经验模态分解的心肌超声造影图像生理参数测量方法 |
CN105030279A (zh) * | 2015-06-24 | 2015-11-11 | 华南理工大学 | 一种基于超声射频时间序列的组织定征方法 |
-
2017
- 2017-06-06 CN CN201710417075.3A patent/CN107358156B/zh not_active Expired - Fee Related
Patent Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102855623A (zh) * | 2012-07-19 | 2013-01-02 | 哈尔滨工业大学 | 基于经验模态分解的心肌超声造影图像生理参数测量方法 |
CN105030279A (zh) * | 2015-06-24 | 2015-11-11 | 华南理工大学 | 一种基于超声射频时间序列的组织定征方法 |
Non-Patent Citations (3)
Title |
---|
Bouaguel Waad 等.A three-stage feature selection using quadratic programming for credit scoring.《Applied Artificial Intelligence》.2013,第27卷(第8期), * |
Nishant Uniyal 等.Ultrasound RF Time Series for Classification of Breast Lesions.《IEEE TRANSACTIONS ON MEDICAL IMAGING》.2015,第34卷(第2期), * |
R. Jurkonis 等.Algorithms and Results of Eye Tissues Differentiation Based on RF Ultrasound.《The Scientific World Journal》.2012,第2012卷(第2期), * |
Also Published As
Publication number | Publication date |
---|---|
CN107358156A (zh) | 2017-11-17 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN107358156B (zh) | 基于希尔伯特-黄变换的超声组织定征的特征提取方法 | |
EP3419516B1 (en) | Ultrasound blood flow imaging | |
Zhou et al. | High spatial–temporal resolution reconstruction of plane-wave ultrasound images with a multichannel multiscale convolutional neural network | |
JP5433097B2 (ja) | 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム | |
CN107346541B (zh) | 一种基于超声射频时间序列小波分析的组织定征方法 | |
EP2904975B1 (en) | Ultrasound observation device, operation method for ultrasound observation device, and operation program for ultrasound observation device | |
Besson et al. | A sparse reconstruction framework for Fourier-based plane-wave imaging | |
CN104240203A (zh) | 基于小波变换和快速双边滤波的医学超声图像去噪方法 | |
CN104411250B (zh) | 超声波观测装置、超声波观测装置的动作方法 | |
KR101610874B1 (ko) | 공간 일관성 기초 초음파 신호 처리 모듈 및 그에 의한 초음파 신호 처리 방법 | |
CN109003232B (zh) | 基于频域尺度平滑Shearlet的医学MRI图像去噪方法 | |
CN103381096A (zh) | 骨表微血管血流灌注分离检测与成像方法 | |
Pham et al. | Joint blind deconvolution and robust principal component analysis for blood flow estimation in medical ultrasound imaging | |
CN111248858B (zh) | 一种基于频域波数域的光声断层成像重建方法 | |
JPWO2012063930A1 (ja) | 超音波診断装置、超音波診断装置の作動方法および超音波診断装置の作動プログラム | |
CN110840484B (zh) | 自适应匹配最优声速的超声成像方法、装置及超声设备 | |
US9386965B2 (en) | Ultrasound machine providing composite image data | |
Chahuara et al. | Regularized framework for simultaneous estimation of ultrasonic attenuation and backscatter coefficients | |
Sultan et al. | Estimation of Cortical Bone Strength Using CNN-based Regression Model | |
Tsui et al. | An adaptive threshold filter for ultrasound signal rejection | |
Wang et al. | An easily-achieved time-domain beamformer for ultrafast ultrasound imaging based on compressive sensing | |
CN110869799A (zh) | 用于处理超声图像的方法和*** | |
Hansen et al. | In-vivo validation of fast spectral velocity estimation techniques | |
CN106037799B (zh) | 基于超声rf背散射信号时频分析的弹性参数成像方法 | |
Wang et al. | Gaussian wavelet based dynamic filtering (GWDF) method for medical ultrasound systems |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant | ||
CF01 | Termination of patent right due to non-payment of annual fee | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20200519 |