CN107480637A - 基于心音特征的心衰分期方法 - Google Patents

基于心音特征的心衰分期方法 Download PDF

Info

Publication number
CN107480637A
CN107480637A CN201710698340.XA CN201710698340A CN107480637A CN 107480637 A CN107480637 A CN 107480637A CN 201710698340 A CN201710698340 A CN 201710698340A CN 107480637 A CN107480637 A CN 107480637A
Authority
CN
China
Prior art keywords
heart failure
stages
signal
heart
phases
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.)
Granted
Application number
CN201710698340.XA
Other languages
English (en)
Other versions
CN107480637B (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.)
Chongqing University
Original Assignee
Chongqing University
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 Chongqing University filed Critical Chongqing University
Priority to CN201710698340.XA priority Critical patent/CN107480637B/zh
Publication of CN107480637A publication Critical patent/CN107480637A/zh
Application granted granted Critical
Publication of CN107480637B publication Critical patent/CN107480637B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/08Feature extraction
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B7/00Instruments for auscultation
    • A61B7/02Stethoscopes
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/0002Inspection of images, e.g. flaw detection
    • G06T7/0012Biomedical image inspection
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06VIMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V40/00Recognition of biometric, human-related or animal-related patterns in image or video data
    • G06V40/10Human or animal bodies, e.g. vehicle occupants or pedestrians; Body parts, e.g. hands
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/30Subject of image; Context of image processing
    • G06T2207/30004Biomedical image processing
    • G06T2207/30048Heart; Cardiac
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06VIMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V40/00Recognition of biometric, human-related or animal-related patterns in image or video data
    • G06V40/10Human or animal bodies, e.g. vehicle occupants or pedestrians; Body parts, e.g. hands
    • G06V40/15Biometric patterns based on physiological signals, e.g. heartbeat, blood flow

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • General Physics & Mathematics (AREA)
  • General Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Public Health (AREA)
  • Biomedical Technology (AREA)
  • General Engineering & Computer Science (AREA)
  • Signal Processing (AREA)
  • Molecular Biology (AREA)
  • Surgery (AREA)
  • Animal Behavior & Ethology (AREA)
  • Artificial Intelligence (AREA)
  • Acoustics & Sound (AREA)
  • Veterinary Medicine (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Radiology & Medical Imaging (AREA)
  • Quality & Reliability (AREA)
  • Human Computer Interaction (AREA)
  • Multimedia (AREA)
  • Measuring Pulse, Heart Rate, Blood Pressure Or Blood Flow (AREA)

Abstract

本发明提供了一种基于心音特征的心衰分期方法,其利用心音信号能够实时反映心脏生理和病理参数的特性,提取其相关的时频特征值并对其求熵构成特征向量,结合用于心衰分期识别的支持向量机分类模型进行心衰分期的辅助识别。由于心音信号属于无创信号,使得该方法方便快捷、成本低廉,能够为临床上心衰分期提供有效的参考信息;同时香农熵是从全局意义上评价信号指标,这样可以使特征向量更加稳定、有效,因此得到的心衰分期参考结果能够更好的结合其它心衰相关信息,帮助提高心衰分期判断的准确性。

Description

基于心音特征的心衰分期方法
技术领域
本发明涉及心音信号分析处理技术领域,具体涉及一种基于心音特征的心衰分期方法。
背景技术
心力衰竭是一种非常严重的心血管疾病,严重威胁着人类健康。在我国,心力衰竭的患病率低于西方国家,但却是我国唯一一种发病率、患病率和死亡率均逐年升高的心血管疾病,其病死率很高。但相关研究表明:若慢性心衰在早期得到及时的治疗,其病死率可以大大降低。所以,对患者心衰进行正确分期对于降低心衰死亡率具有非常重要的临床意义。根据美国心脏病协会,慢性心衰分为A、B、C、D四个时期,A、B期可以得到有效控制,C、D期较为严重,但通过治疗可以减轻症状。
目前临床上主要以病史、体格检查的结果为基础,结合超声心动图、胸部X线、心电图等辅助检查,从而完成对慢性心衰的分期判断。而使用较多的方法主要有:
(1)超声心动图。能够显示各心腔大小变化和心瓣膜结构,通过计算左室射血分数(LVEF)、舒张早期与晚期心室充盈速度最大值的比值(E/A)评估心脏的心功能状态并对慢性心衰进行分期。准确性较高,但价格昂贵,且部分指标的敏感性较差,无法对LVEF正常的心力衰竭患者的心脏收缩能力进行评估。
(2)NYHA心功能分级。NYHA心功能分级是根据心衰患者的活动耐量及症状分为四级:I级为日常活动无心力衰竭症状;II级为日常活动出现心力衰竭症状(呼吸困难、乏力);III级为低于正常活动出现心力衰竭症状;IV级为休息时出现心力衰竭症状。此方法简单易行,可用于慢性心衰的初步检查。但结果较为主观,评价准确性取决于患者的自我描述,缺乏客观性,同时还存在个体差异等缺陷。此外,该方法不能明确整个心衰过程,具有一定的局限性。
(3)机器学习法。通过选择病人的病史、性别、年龄等特征变量,并通过权值估计筛选出最佳的几个变量作为人工智能分类器的输入变量对心衰的严重程度进行评估和分期,也可以利用心电图中的心率变异性、每搏输出量等指标作为分类器的输入向量进行自动分期识别。该方法与之前的方法相比更加便捷、稳定。但是性别、年龄等变量的特征性不够充分,容易受个体差异的影响;心率变异性和每搏输出量等指标获取较为困难。
综上所述,现有的心力衰竭分期方法存在着分期便捷性不佳、稳定性不足或者准确率不高等弊端,在某些特定的情况下无法达到理想的状态,而心音作为一种心脏内的机械振动信号,获取便捷、无创,能实时有效的反映心脏的生理和病理信息,当心衰尚未发展到足以产生临床及病理改变之前,心音就会产生一定的变化,从而反映慢性心衰的病变情况。因此,提取心音的相关特征值可用于进行心衰分期的辅助识别。此外,利用熵值分析法处理信号特征值可以从全局意义上获取信号信息,定量描述信号的分布特性,使信号的识别更加准确;但目前心音仅用于心脏内异常心音的识别分类,而未用于心力衰竭的心衰分期辅助判断。
发明内容
针对现有技术中存在的不足,本发明的目的在于提供一种基于心音特征的心衰分期方法,为心衰分期识别提供一种无创、便捷、稳定性强的心衰分期辅助判断手段,从而能够用以结合其它心衰相关信息帮助提高心衰分期判断的准确性。
为实现上述目的,本发明采用了如下技术方案:
基于心音特征的心衰分期方法,用于对待测心音信号进行特征向量的提取,并根据所述特征向量识别得到待测心音信号的心衰分期参考结果;该方法包括如下步骤:
1)对待测心音信号X(n)进行预处理,得到处理后的信号XT(k),n=1,2,…,N,N表示待测心音信号X(n)的采样点总个数,k=1,2,…,K,K表示处理后的信号XT(k)的采样点总个数;
2)对预处理后的信号XT(k)进行小波包分解,得到子频带小波包系数;分别对各子频带小波包系数进行重构,得到重构后的各个子频带信号Dg(k),g∈{1,2,…,m},m为子频带信号总个数;
3)将所述m个子频带信号Dg(k)构成m×K的时间幅值矩阵,然后根据奇异值分解理论,对时间幅值矩阵进行奇异值分解得到若干个奇异值并求熵,得到相应的心音奇异熵值;利用能量谱公式分别计算m个子频带信号Dg(k)的能量谱以及重构信号XT(k)的总能量谱并求熵,得到相应的心音能量熵值;利用welch功率谱估计法分别估算m个子频带信号Dg(k)的功率谱并求熵,得到m个对应的心音功率谱熵值;
4)从得到的m个心音功率谱熵值中择优选取a个心音功率谱熵值,0<a≤m;
5)将选取的a个心音功率谱熵值连同所述心音奇异熵值和心音能量熵值的集合作为待测心音信号的特征向量
6)获取训练得到的用于心衰分期识别的支持向量机分类模型,对待测心音信号的特征向量进行分类识别,判定待测心音信号所属的心衰分期,得到心衰分期参考结果。
上述基于心音特征的心衰分期方法中,进一步,所述步骤1)中心音信号预处理的具体步骤为:
1.1)对待测心音信号X(n)进行多倍重采样得到重采样信号,再对重采样信号进行滤波得到滤波信号;
1.2)从所述滤波信号中截取一段至少包含两个心音周期的滤波信号段X(k)并进行小波分解,得到分解的各层细节系数di(k)和一个近似系数a(k),其中,k表示截取的滤波信号段中的采样点位置,k=1,2,…,K,K表示截取的滤波信号段X(k)的采样点总个数,i∈{1,2,…,I},I表示对截取的滤波信号段X(k)进行小波分解的总层数;
1.3)对各层细节系数di(k)进行阈值去噪处理,得到新的细节系数,再结合所述近似系数a(k)进行信号重构得到预处理后的信号XT(k)。
上述基于心音特征的心衰分期方法中,进一步,所述步骤1.3)中,对第i层细节系数di(k)进行阈值去噪处理所设定的去噪阈值δi为:
其中,σi为小波分解的第i层细节系数di(k)的噪声标准差,且符号|·|为绝对值运算符,median(|di(k)|)表示第i层细节系数di(k)绝对值的中值。
上述基于心音特征的心衰分期方法中,进一步,所述步骤2)具体为:
2.1)对预处理后的信号XT(k)进行b层小波包分解,得到2b个子频带小波包系数;
2.2)分别对2b个子频带小波包系数进行信号重构,得到m=2b个子频带信号Dg(k),g∈{1,2,…,m}。
上述基于心音特征的心衰分期方法中,进一步,所述步骤3)中获取心音奇异熵值的方式具体为:
3.1.1)将m个子频带信号Dg(k)构成m×K的时间幅值矩阵Dm×K
3.1.2)根据信号的奇异值分解理论,将时间幅值矩阵Dm×K分解为如下形式:
Dm×K=Um×l·Λl×l·Vl×K
所述Um×l和Vl×K分别为对时间幅值矩阵Dm×K进行奇异值分解的m×l和l×K的酉矩阵,Λl×l为l×l的对角矩阵,且该对角矩阵Λl×l中的l个主对角线元素λr即为时间幅值矩阵Dm×K的各奇异值,r∈{1,2,…,l};
3.1.3)计算心音奇异熵值WS
其中,Δpv表示第v阶增量小波奇异值,v∈{1,2,…,L},L为对角矩阵Λl×l中非零主对角线元素的个数,且:
λv表示对角矩阵Λl×l中的第v个非零主对角线元素;
所述步骤3)中获取心音能量熵值的方式具体为:
3.2.1)利用能量谱公式分别计算m个子频带信号Dg(k)的能量谱Eg
3.2.2)重构计算信号XT(k)的总能量谱并求熵,得到相应的心音能量熵值WE
所述步骤3)中获取心音功率谱熵值的方式具体为:
3.3.1)针对第g个子频带信号Dg(k),g∈{1,2,…,m},将其拆分为J个信号段dg,c(s),c∈{1,2,…,J},且允许相邻两个信号段之间存在重叠,每个信号段dg,c(s)包含M个采样点,s∈{1,2,…,M},从而按照Bartlett法分别估算每个子频带信号段dg,c(s)的功率谱Pg,c
其中,d2(s)为对子频带信号进行分段拆分的窗函数,e为自然常数;U为功率谱归一化因子,且
由此,分别估算得到各个子频带信号Dj(k)拆分的各个子频带信号段的功率谱;
3.3.2)针对第g个子频带信号Dg(k),g∈{1,2,…,m},计算其对应的心音功率谱熵值Wg
由此,分别计算得到各个子频带信号Dg(k)对应的功率谱心音功率谱熵值。
上述基于心音特征的心衰分期方法中,进一步,所述步骤4)中择优选取a个心音功率谱熵值的具体方式为:使用ReliefF方法的权重值更新方式对所得到的m个心音功率谱熵值的权值向量W(τ)进行多次迭代,选取经迭代后的权值向量W(τ)中较大权重值对应的a个心音功率谱熵值,0<a≤m。
上述基于心音特征的心衰分期方法中,进一步,所述步骤6)中用于心衰分期识别的支持向量机分类模型的训练方法为:
A)预先获取心衰分期分别属于A期、B期、C期、D期的若干个已知心音信号;
B)设定当前分期识别类别为心衰A期;
C)将获取的各已知心音信号随机划分作为训练样本和测试样本;分别获取各个训练样本的特征向量并记录每个训练样本各自所属的心衰分期标记得到训练样本集p∈{1,2,…,P},P表示训练样本的总个数;分别获取各个测试样本的特征向量构成测试样本集并记录每个测试样本各自所属的心衰分期标记q∈{1,2,…,Q},Q表示测试样本的总个数;其中,对于所属心衰分期为当前分期识别类别的训练样本和测试样本,其心衰分期标记的取值为1,对于所属心衰分期不为当前分期识别类别的训练样本和测试样本,其心衰分期标记的取值为-1;
D)利用训练样本集对一个支持向量机分类器进行训练;
E)利用测试样本集对训练得到的支持向量机分类器进行分类识别测试,分别得到所述支持向量机分类器对各个测试样本的心衰分期测试结果,并根据记录的各测试样本所属的心衰分期标记,统计心衰分期正确的测试识别率;若测试识别率未达到预设的预期目标识别率,则返回执行步骤D)进行重复训练;若测试识别率达到预设的预期目标识别率,则判定支持向量机分类器训练完成;
F)将训练完成的支持向量机分类器作为心衰A期的支持向量机分类器;
G)分别设定当前分期识别类别为心衰B期、C期和D期并执行步骤C)~E),从而分别得到心衰B期、C期和D期的支持向量机分类器;
H)将得到的心衰A期、B期、C期和D期的支持向量机分类器的集合作为训练得到的用于心衰分期识别的支持向量机分类模型。
上述基于心音特征的心衰分期方法中,进一步,所述训练得到的用于心衰分期识别的支持向量机分类模型中包含用于识别心音信号属于心衰A期、B期、C期和D期的四个支持向量机分类器;所述步骤6)中对待测心音信号的特征向量进行分类识别的具体方式为:
6.1)利用心衰A期的支持向量机分类器对待测心音信号进行分类识别,若识别判定待测心音信号属于心衰A期,则跳转执行步骤6.5);否则,执行步骤6.2);
6.2)利用心衰B期的支持向量机分类器对待测心音信号进行分类识别,若识别判定待测心音信号属于心衰B期,则跳转执行步骤6.5);否则,执行步骤6.3);
6.3)利用心衰C期的支持向量机分类器对待测心音信号进行分类识别,若识别判定待测心音信号属于心衰C期,则跳转执行步骤6.5);否则,执行步骤6.4);
6.4)利用心衰D期的支持向量机分类器对待测心音信号进行分类识别,若识别判定待测心音信号属于心衰D期,则执行步骤6.5);否则,判定未得到待测心音信号的心衰分期识别结果,然后执行步骤6.5);
6.5)根据识别结果确定待测心音信号的心衰分期参考结果。
相比于现有技术,本发明具有如下有益效果:
1、本发明基于心音特征的心衰分期方法,利用心音信号能够实时反映心脏生理和病理参数的特性,提取其相关的时频特征值构成特征向量,结合用于心衰分期识别的支持向量机分类模型进行心衰分期的辅助识别;由于心音信号属于无创信号,使得该方法方便快捷、成本低廉,能够为心衰分期临床上提供有效的参考信息。
2、本发明基于心音特征的心衰分期方法中,在提取在心音信号的相关特征值之后,使用香农信息熵的方法对特征值进行求熵处理;由于香农熵是从全局意义上评价信号指标,这样可以使特征向量更加稳定、有效。
3、利用本发明基于心音特征的心衰分期方法对心音信号进行处理识别得到的心衰分期参考结果,能够更好的结合其它心衰相关信息,帮助提高心衰分期判断的准确性。
附图说明
图1为本发明基于心音特征的心衰分期方法的流程图。
图2为心衰A期的心音信号示例图。
图3为心衰B期的心音信号示例图。
图4为心衰C期的心音信号示例图。
图5为心衰D期的心音信号示例图。
具体实施方式
本发明提供一种基于心音特征的心衰分期方法,为心衰分期识别提供一种无创、便捷、稳定性强的心衰分期辅助判断手段,能够用以结合其它心衰相关信息帮助提高心衰分期判断的准确性。
本发明基于心音特征的心衰分期方法流程如图1所示。首先获取通过心音传感器以设定的采样频率采样得到的作为待测对象的心音信号,作为待测心音信号,然后采用本发明方法对待测心音信号进行特征向量的提取,并根据所述特征向量识别得到待测心音信号的心衰分期参考结果;其流程具体包括如下步骤:
1)对待测心音信号X(n)进行预处理,得到处理后的信号XT(k),n=1,2,…,N,N表示待测心音信号X(n)的采样点总个数,k=1,2,…,K,K表示处理后的信号XT(k)的采样点总个数。
该步骤中对心音信号进行预处理的具体方式为:
1.1)先对待测心音信号X(n)进行多倍重采样得到重采样信号,再对重采样信号进行滤波得到滤波信号。
例如,对于原始采样频率为11025Hz的待测心音信号X(n)进行5倍重采样,得到2205Hz的重采样信号,然后再用高通滤波器和陷波滤波器对重采样信号进行滤波得到滤波信号。
1.2)为了减少数据处理量,从所述滤波信号中截取一段至少包含两个心音周期的滤波信号段X(k)并进行小波分解,得到分解的各层细节系数di(k)和一个近似系数a(k),其中,k表示截取的滤波信号段中的采样点位置,k=1,2,…,K,K表示截取的滤波信号段X(k)的采样点总个数,i∈{1,2,…,I},I表示对截取的滤波信号段X(k)进行小波分解的总层数。之所以截取的滤波信号段X(k)需要至少包含两个心音周期,是为了保证滤波信号段X(k)内包含有较完整的特征信息。
例如,对于2205Hz的滤波信号,其中每个心音周期大概包含1800个左右的采样点,从中截取至少包含两个心音周期的滤波信号段X(k)大概至少包含4000个左右的采样点;然后,可以使用coif3小波基对滤波信号段X(k)进行5层离散小波分解,得到5层细节系数di(k)和一个近似系数a(k),此时K=5,即i∈{1,2,3,4,5}。
1.3)对各层细节系数di(k)进行阈值去噪处理,得到新的细节系数,再结合所述近似系数a(k)进行信号重构得到预处理后的信号XT(k)。
其中,对任意的第i层细节系数di(k)进行阈值去噪处理所设定的去噪阈值δi可以选择为:
由于滤波信号段X(k)中的噪声,随着小波分解层数的增加,其噪声系数会逐渐减小,所以在全局去噪阈值中添加分母ln(1+i)。其中,σi为小波分解的第i层细节系数di(k)的噪声标准差,且符号|·|为绝对值运算符,median(|di(k)|)表示第i层细节系数di(k)绝对值的中值。
2)对预处理后的信号XT(k)进行小波包分解,得到子频带小波包系数;分别对各子频带小波包系数进行重构,得到重构后的各个子频带信号Dg(k),g∈{1,2,…,m},m为子频带信号总个数。
该步骤具体为:
2.1)对预处理后的信号XT(k)进行b层小波包分解,得到2b个子频带小波包系数。例如进行4层小波包分解,则得到16个子频带小波包系数
2.2)分别对2b个子频带小波包系数进行信号重构,得到m=2b个子频带信号Dg(k),g∈{1,2,…,m}。例如,对16个子频带小波包系数进行信号重构,即得到16个子频带信号Dj(k),此时m=16,即j∈{1,2,…,16}。
3)将所述m个子频带信号Dg(k)构成m×K的时间幅值矩阵,然后根据奇异值分解理论,对时间幅值矩阵进行奇异值分解得到若干个奇异值并求熵,得到相应的心音奇异熵值;利用能量谱公式分别计算m个子频带信号Dg(k)的能量谱以及重构信号XT(k)的总能量谱并求熵,得到相应的心音能量熵值;利用welch功率谱估计法分别估算m个子频带信号Dg(k)的功率谱并求熵,得到m个对应的心音功率谱熵值。
该步骤中,需要分别获取到心音奇异熵值、心音能量熵值和心音功率谱熵值。其中:
3.1、获取心音奇异熵值的方式具体为:
3.1.1)将m个子频带信号Dg(k)构成m×K的时间幅值矩阵Dm×K
3.1.2)根据信号的奇异值分解理论,将时间幅值矩阵Dm×K分解为如下形式:
Dm×K=Um×l·Λl×l·Vl×K
所述Um×l和Vl×K分别为对时间幅值矩阵Dm×K进行奇异值分解的m×l和l×K的酉矩阵,Λl×l为l×l的对角矩阵,且该对角矩阵Λl×l中的l个主对角线元素λr即为时间幅值矩阵Dm×K的各奇异值,r∈{1,2,…,l}。
3.1.3)为了定量描述信号的频率成分及分布特性,定义心音奇异熵值WS的计算式为:
其中,Δpv表示第v阶增量小波奇异值,v∈{1,2,…,L},L为对角矩阵Λl×l中非零主对角线元素的个数,且:
λv表示对角矩阵Λl×l中的第v个非零主对角线元素。
3.2、获取心音能量熵值的方式具体为:
3.2.1)利用能量谱公式分别计算m个子频带信号Dg(k)的能量谱Eg
3.2.2)信号的总能量等于各尺度下分解信号能量之和,因此,重构计算信号XT(k)的总能量谱并求熵,得到相应的心音能量熵值WE
3.3、由于心音信号是一种非线性非平稳的随机性信号,无法用解析信号精确地表示,其功率谱不能直接求得,所以可以采用Welch法来估计通过小波包变换的各尺度分解信号的功率谱。获取心音功率谱熵值的方式具体为:
3.3.1)针对第g个子频带信号Dg(k),g∈{1,2,…,m},将其拆分为J个信号段dg,c(s),c∈{1,2,…,J},且允许相邻两个信号段之间存在重叠,每个信号段dg,c(s)包含M个采样点,s∈{1,2,…,M},从而按照Bartlett法分别估算每个子频带信号段dg,c(s)的功率谱Pg,c
其中,d2(s)为对子频带信号进行分段拆分的窗函数,e为自然常数;U为功率谱归一化因子,且
由此,分别估算得到各个子频带信号Dg(k)拆分的各个子频带信号段的功率谱。
3.3.2)针对第g个子频带信号Dg(k),g∈{1,2,…,m},计算其对应的心音功率谱熵值Wg
由此,分别计算得到各个子频带信号Dg(k)对应的功率谱心音功率谱熵值。例如,若存在16个子频带信号Dg(k),则相应计算得到16个功率谱心音功率谱熵值。
4)从得到的m个心音功率谱熵值中择优选取a个心音功率谱熵值,0<a≤m。
若在步骤2)中对预处理后的信号XT(k)进行小波包分解的层数较多,则步骤3)中处理得到的功率谱心音功率谱熵值的个数也就相应较多。因此,若为了减少心音功率谱熵值作为心音信号的分期识别特征进行数据处理的数据量,那么可以从得到的m个心音功率谱熵值中择优选取a个心音功率谱熵值。至于择优选择的具体方式,可以根据应用情况采用不同的择优方案。例如,可以直接选择得到的全部功率谱心音功率谱熵值,即选择的心音功率谱熵值个数a=m;或者,可以使用ReliefF方法的权重值更新方式对所得到的m个心音功率谱熵值的权值向量W(τ)进行多次迭代,选取经迭代后的权值向量W(τ)中较大权重值对应的a个心音功率谱熵值,0<a≤m。使用ReliefF方法的权重值更新方式进行权值向量迭代的具体方式属于现有技术,并不是本发明的技术创新点,因此这里不再展开论述。
5)将选取的a个心音功率谱熵值连同所述心音奇异熵值和心音能量熵值的集合作为待测心音信号的特征向量
例如,若择优选取了6个心音功率谱熵值(即a=6),分别记为WP1、WP2、WP3、WP4、WP5和WP6,则连同所得到的心音奇异熵值WS和心音能量熵值WE,构成的心音信号的特征向量则为由此,便可以利用得到待测心音信号的特征向量,用以进行心衰分期识别。
6)获取训练得到的用于心衰分期识别的支持向量机分类模型,对待测心音信号的特征向量进行分类识别,判定待测心音信号所属的心衰分期,得到心衰分期参考结果。
该步骤中,需要获取到预先训练得到的用于心衰分期识别的支持向量机分类模型,用于对待测心音信号进行心衰分期参考结果的分类识别。预先训练得到的用于心衰分期识别的支持向量机分类模型,可以经过训练得到后加以存储,以备使用时调用获取。
具体而言,用于心衰分期识别的支持向量机分类模型的训练方法为:
A)预先获取心衰分期分别属于A期、B期、C期、D期的若干个已知心音信号。
在训练支持向量机分类模型的过程中,预先获取的已知心音信号用以划分作为训练样本和测试样本使用,训练样本用于实施训练,而测试样本则用于测试训练效果是否达到预期。训练样本和测试样本都需要包含A期、B期、C期、D期的已知心音信号。这些已知心音信号可以通过采用心音信号采集设备分别对已知为心衰A期、心衰B期、心衰C期和心衰D期的不同心衰患者进行心音信号采集而获得,针对每一种心衰分期采集的已知心音信号都需要采集若干个,且针对每一种心衰分期采集的已知心音信号数量越多,则越有利于在训练中使用。采集到的心衰A期、心衰B期、心衰C期、心衰D期的心音信号示例分别如图2、图3、图4、图5所示,可以看到,不同心衰分期的心音信号有所不同。
B)设定当前分期识别类别为心衰A期。
本发明采用了支持向量机分类模型,而单一的支持向量机分类器仅能够实现二分类(即判断是、否属于某一类别),因此本发明采用的支持向量机分类模型中需要使用多个支持向量机分类器来组合构成,所以需要针对每个不同的心衰分期分别进行训练。这里先设定当前分期识别类别为心衰A期,即表示先设定针对心衰A期进行支持向量机分类器的训练。
C)将获取的各已知心音信号随机划分作为训练样本和测试样本;分别获取各个训练样本的特征向量并记录每个训练样本各自所属的心衰分期标记得到训练样本集p∈{1,2,…,P},P表示训练样本的总个数;分别获取各个测试样本的特征向量构成测试样本集并记录每个测试样本各自所属的心衰分期标记q∈{1,2,…,Q},Q表示测试样本的总个数;其中,对于所属心衰分期为当前分期识别类别的训练样本和测试样本,其心衰分期标记的取值为1,对于所属心衰分期不为当前分期识别类别的训练样本和测试样本,其心衰分期标记的取值为-1。
该步骤中,提取每个训练样本和每个测试样本的特征向量的方法,可以参照前述步骤1)~5)对待测心音信号的特征向量提取方法来进行,将针对每个训练样本或测试样本提取选择的a个心音功率谱熵值连同提取的心音奇异熵值和心音能量熵值的集合作为相应训练样本或测试样本的特征向量。
D)利用训练样本集对一个支持向量机分类器进行训练。
这里,利用已知类别并带有分类标记的训练样本训练SVM分类器,目的是使得SVM分类器能够完成针对当前分期识别类别的心音信号分类识别;例如,当前分期识别类别为心衰A期,则训练样本集中仅属于心衰A期的已知心音信号训练样本的心衰分期标记取值为1,其它分期类别的训练样本的心衰分期标记取值为-1,由此经过训练,使得当前训练的SVM分类器能够针对心衰A期的心音信号和其它心衰分期的心音信号加以区分识别。而利用已知类别并带有分类标记的训练样本训练SVM分类器的具体方法流程属于现有技术,并不是本发明的技术创新点,因此这里不再展开论述。
E)利用测试样本集对训练得到的支持向量机分类器进行分类识别测试,分别得到所述支持向量机分类器对各个测试样本的心衰分期测试结果,并根据记录的各测试样本所属的心衰分期标记,统计心衰分期正确的测试识别率;若测试识别率未达到预设的预期目标识别率,则返回执行步骤D)进行重复训练;若测试识别率达到预设的预期目标识别率,则判定支持向量机分类器训练完成。
测试样本则用于测试对支持向量机分类器的训练识别效果是否达到预期目标识别率,若未达到,则重复训练,直至达到预期目标识别率为止。
F)将训练完成的支持向量机分类器作为心衰A期的支持向量机分类器。
在上述步骤中,在设定当前分期识别类别为心衰A期的情况下,经过训练完成而得到的支持向量机分类器则作为心衰A期的支持向量机分类器。
G)分别设定当前分期识别类别为心衰B期、C期和D期并执行步骤C)~E),从而分别得到心衰B期、C期和D期的支持向量机分类器。
分别设定当前分期识别类别为B期、C期和D期,并执行步骤C)~E)进行训练完成后得到的支持向量机分类器,则分别为心衰B期、心衰C期和心衰D期的支持向量机分类器。
H)将得到的心衰A期、B期、C期和D期的支持向量机分类器的集合作为训练得到的用于心衰分期识别的支持向量机分类模型。
由此,便训练得到了用于心衰分期识别的支持向量机分类模型。利用训练得到的用于心衰分期识别的支持向量机分类模型中对待测心音信号的特征向量进行分类识别的具体方式为:
6.1)利用心衰A期的支持向量机分类器对待测心音信号进行分类识别,若识别判定待测心音信号属于心衰A期,则跳转执行步骤6.5);否则,执行步骤6.2);
6.2)利用心衰B期的支持向量机分类器对待测心音信号进行分类识别,若识别判定待测心音信号属于心衰B期,则跳转执行步骤6.5);否则,执行步骤6.3);
6.3)利用心衰C期的支持向量机分类器对待测心音信号进行分类识别,若识别判定待测心音信号属于心衰C期,则跳转执行步骤6.5);否则,执行步骤6.4);
6.4)利用心衰D期的支持向量机分类器对待测心音信号进行分类识别,若识别判定待测心音信号属于心衰D期,则执行步骤6.5);否则,判定未得到待测心音信号的心衰分期识别结果,然后执行步骤6.5);
6.5)根据识别结果确定待测心音信号的心衰分期参考结果。
综上所述,本发明基于心音特征的心衰分期方法,利用心音信号能够实时反映心脏生理和病理参数的特性,提取其相关的时频特征值构成特征向量,结合用于心衰分期识别的支持向量机分类模型进行心衰分期的辅助识别;由于心音信号属于无创信号,使得该方法方便快捷、成本低廉,能够为心衰分期临床上提供有效的参考信息。本发明基于心音特征的心衰分期方法中,在提取在心音信号的相关特征值之后,使用香农信息熵的方法对特征值进行求熵处理;由于香农熵是从全局意义上评价信号指标,这样可以使特征向量更加稳定、有效。因此,利用本发明基于心音特征的心衰分期方法对心音信号进行处理识别得到的心衰分期参考结果,能够更好的结合其它心衰相关信息,帮助提高心衰分期判断的准确性。
最后说明的是,以上实施例仅用以说明本发明的技术方案而非限制,尽管参照较佳实施例对本发明进行了详细说明,本领域的普通技术人员应当理解,可以对本发明的技术方案进行修改或者等同替换,而不脱离本发明技术方案的宗旨和范围,其均应涵盖在本发明的权利要求范围当中。

Claims (8)

1.基于心音特征的心衰分期方法,其特征在于,用于对待测心音信号进行特征向量的提取,并根据所述特征向量识别得到待测心音信号的心衰分期参考结果;该方法包括如下步骤:
1)对待测心音信号X(n)进行预处理,得到处理后的信号XT(k),n=1,2,…,N,N表示待测心音信号X(n)的采样点总个数,k=1,2,…,K,K表示处理后的信号XT(k)的采样点总个数;
2)对预处理后的信号XT(k)进行小波包分解,得到子频带小波包系数;分别对各子频带小波包系数进行重构,得到重构后的各个子频带信号Dg(k),g∈{1,2,…,m},m为子频带信号总个数;
3)将所述m个子频带信号Dg(k)构成m×K的时间幅值矩阵,然后根据奇异值分解理论,对时间幅值矩阵进行奇异值分解得到若干个奇异值并求熵,得到相应的心音奇异熵值;利用能量谱公式分别计算m个子频带信号Dg(k)的能量谱以及重构信号XT(k)的总能量谱并求熵,得到相应的心音能量熵值;利用welch功率谱估计法分别估算m个子频带信号Dg(k)的功率谱并求熵,得到m个对应的心音功率谱熵值;
4)从得到的m个心音功率谱熵值中择优选取a个心音功率谱熵值,0<a≤m;
5)将选取的a个心音功率谱熵值连同所述心音奇异熵值和心音能量熵值的集合作为待测心音信号的特征向量
6)获取训练得到的用于心衰分期识别的支持向量机分类模型,对待测心音信号的特征向量进行分类识别,判定待测心音信号所属的心衰分期,得到心衰分期参考结果。
2.根据权利要求1所述基于心音特征的心衰分期方法,其特征在于,所述步骤1)中心音信号预处理的具体步骤为:
1.1)对待测心音信号X(n)进行多倍重采样得到重采样信号,再对重采样信号进行滤波得到滤波信号;
1.2)从所述滤波信号中截取一段至少包含两个心音周期的滤波信号段X(k)并进行小波分解,得到分解的各层细节系数di(k)和一个近似系数a(k),其中,k表示截取的滤波信号段中的采样点位置,k=1,2,…,K,K表示截取的滤波信号段X(k)的采样点总个数,i∈{1,2,…,I},I表示对截取的滤波信号段X(k)进行小波分解的总层数;
1.3)对各层细节系数di(k)进行阈值去噪处理,得到新的细节系数,再结合所述近似系数a(k)进行信号重构得到预处理后的信号XT(k)。
3.根据权利要求2所述基于心音特征的心衰分期方法,其特征在于,所述步骤1.3)中,对第i层细节系数di(k)进行阈值去噪处理所设定的去噪阈值δi为:
其中,σi为小波分解的第i层细节系数di(k)的噪声标准差,且符号|·|为绝对值运算符,median(|di(k)|)表示第i层细节系数di(k)绝对值的中值。
4.根据权利要求1所述基于心音特征的心衰分期方法,其特征在于,所述步骤2)具体为:
2.1)对预处理后的信号XT(k)进行b层小波包分解,得到2b个子频带小波包系数;
2.2)分别对2b个子频带小波包系数进行信号重构,得到m=2b个子频带信号Dg(k),g∈{1,2,…,m}。
5.根据权利要求1所述基于心音特征的心衰分期方法,其特征在于所述,所述步骤3)中获取心音奇异熵值的方式具体为:
3.1.1)将m个子频带信号Dg(k)构成m×K的时间幅值矩阵Dm×K
3.1.2)根据信号的奇异值分解理论,将时间幅值矩阵Dm×K分解为如下形式:
Dm×K=Um×l·Λl×l·Vl×K
所述Um×l和Vl×K分别为对时间幅值矩阵Dm×K进行奇异值分解的m×l和l×K的酉矩阵,Λl×l为l×l的对角矩阵,且该对角矩阵Λl×l中的l个主对角线元素λr即为时间幅值矩阵Dm×K的各奇异值,r∈{1,2,…,l};
3.1.3)计算心音奇异熵值WS
其中,Δpv表示第v阶增量小波奇异值,v∈{1,2,…,L},L为对角矩阵Λl×l中非零主对角线元素的个数,且:
λv表示对角矩阵Λl×l中的第v个非零主对角线元素;
所述步骤3)中获取心音能量熵值的方式具体为:
3.2.1)利用能量谱公式分别计算m个子频带信号Dg(k)的能量谱Eg
3.2.2)重构计算信号XT(k)的总能量谱并求熵,得到相应的心音能量熵值WE
所述步骤3)中获取心音功率谱熵值的方式具体为:
3.3.1)针对第g个子频带信号Dg(k),g∈{1,2,…,m},将其拆分为J个信号段dg,c(s),c∈{1,2,…,J},且允许相邻两个信号段之间存在重叠,每个信号段dg,c(s)包含M个采样点,s∈{1,2,…,M},从而按照Bartlett法分别估算每个子频带信号段dg,c(s)的功率谱Pg,c
其中,d2(s)为对子频带信号进行分段拆分的窗函数,e为自然常数;U为功率谱归一化因子,且
由此,分别估算得到各个子频带信号Dg(k)拆分的各个子频带信号段的功率谱;
3.3.2)针对第g个子频带信号Dg(k),g∈{1,2,…,m},计算其对应的心音功率谱熵值Wg
由此,分别计算得到各个子频带信号Dg(k)对应的功率谱心音功率谱熵值。
6.根据权利要求1所述基于心音特征的心衰分期方法,其特征在于,所述步骤4)中择优选取a个心音功率谱熵值的具体方式为:使用ReliefF方法的权重值更新方式对所得到的m个心音功率谱熵值的权值向量W(τ)进行多次迭代,选取经迭代后的权值向量W(τ)中较大权重值对应的a个心音功率谱熵值,0<a≤m。
7.根据权利要求1所述基于心音特征的心衰分期方法,其特征在于,所述步骤6)中用于心衰分期识别的支持向量机分类模型的训练方法为:
A)预先获取心衰分期分别属于A期、B期、C期、D期的若干个已知心音信号;
B)设定当前分期识别类别为心衰A期;
C)将获取的各已知心音信号随机划分作为训练样本和测试样本;分别获取各个训练样本的特征向量并记录每个训练样本各自所属的心衰分期标记得到训练样本集P表示训练样本的总个数;分别获取各个测试样本的特征向量构成测试样本集并记录每个测试样本各自所属的心衰分期标记 Q表示测试样本的总个数;其中,对于所属心衰分期为当前分期识别类别的训练样本和测试样本,其心衰分期标记的取值为1,对于所属心衰分期不为当前分期识别类别的训练样本和测试样本,其心衰分期标记的取值为-1;
D)利用训练样本集对一个支持向量机分类器进行训练;
E)利用测试样本集对训练得到的支持向量机分类器进行分类识别测试,分别得到所述支持向量机分类器对各个测试样本的心衰分期测试结果,并根据记录的各测试样本所属的心衰分期标记,统计心衰分期正确的测试识别率;若测试识别率未达到预设的预期目标识别率,则返回执行步骤D)进行重复训练;若测试识别率达到预设的预期目标识别率,则判定支持向量机分类器训练完成;
F)将训练完成的支持向量机分类器作为心衰A期的支持向量机分类器;
G)分别设定当前分期识别类别为心衰B期、C期和D期并执行步骤C)~E),从而分别得到心衰B期、C期和D期的支持向量机分类器;
H)将得到的心衰A期、B期、C期和D期的支持向量机分类器的集合作为训练得到的用于心衰分期识别的支持向量机分类模型。
8.根据权利要求1所述基于心音特征的心衰分期方法,其特征在于,所述训练得到的用于心衰分期识别的支持向量机分类模型中包含用于识别心音信号属于心衰A期、B期、C期和D期的四个支持向量机分类器;所述步骤6)中对待测心音信号的特征向量进行分类识别的具体方式为:
6.1)利用心衰A期的支持向量机分类器对待测心音信号进行分类识别,若识别判定待测心音信号属于心衰A期,则跳转执行步骤6.5);否则,执行步骤6.2);
6.2)利用心衰B期的支持向量机分类器对待测心音信号进行分类识别,若识别判定待测心音信号属于心衰B期,则跳转执行步骤6.5);否则,执行步骤6.3);
6.3)利用心衰C期的支持向量机分类器对待测心音信号进行分类识别,若识别判定待测心音信号属于心衰C期,则跳转执行步骤6.5);否则,执行步骤6.4);
6.4)利用心衰D期的支持向量机分类器对待测心音信号进行分类识别,若识别判定待测心音信号属于心衰D期,则执行步骤6.5);否则,判定未得到待测心音信号的心衰分期识别结果,然后执行步骤6.5);
6.5)根据识别结果确定待测心音信号的心衰分期参考结果。
CN201710698340.XA 2017-08-15 2017-08-15 基于心音特征的心衰分期方法 Active CN107480637B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710698340.XA CN107480637B (zh) 2017-08-15 2017-08-15 基于心音特征的心衰分期方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710698340.XA CN107480637B (zh) 2017-08-15 2017-08-15 基于心音特征的心衰分期方法

Publications (2)

Publication Number Publication Date
CN107480637A true CN107480637A (zh) 2017-12-15
CN107480637B CN107480637B (zh) 2019-08-30

Family

ID=60600684

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710698340.XA Active CN107480637B (zh) 2017-08-15 2017-08-15 基于心音特征的心衰分期方法

Country Status (1)

Country Link
CN (1) CN107480637B (zh)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109567856A (zh) * 2018-11-23 2019-04-05 四川长虹电器股份有限公司 一种基于智能平台的听诊位置智能提醒方法
CN109934243A (zh) * 2017-12-19 2019-06-25 中国科学院深圳先进技术研究院 Ecg数据分类方法、装置、电子设备及系统
CN110175996A (zh) * 2019-05-29 2019-08-27 中国人民解放军陆军军医大学第二附属医院 一种基于cmr双向循环网络的心衰自动分级方法
CN116525104A (zh) * 2023-06-26 2023-08-01 中国人民解放军总医院 一种心源性休克快速分期分诊系统、设备及可存储介质
CN116807486A (zh) * 2023-07-05 2023-09-29 淮安市楚州中医院 一种急性心衰的数字模型构件系统及方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1850007A (zh) * 2006-05-16 2006-10-25 清华大学深圳研究生院 基于心音分析的心脏病自动分类系统及其心音分段方法
CN101518439A (zh) * 2009-03-24 2009-09-02 重庆大学 一种基于心音和心电的心脏功能检测系统
CN104778342A (zh) * 2014-01-14 2015-07-15 济南大学 一种基于小波奇异熵的心音特征提取方法
CN105488520A (zh) * 2015-11-24 2016-04-13 昆明理工大学 基于多分辨奇异谱熵和svm的泄漏声发射信号识别方法
CN106226635A (zh) * 2016-07-14 2016-12-14 国网福建晋江市供电有限公司 一种配电网馈线故障类型识别方法及装置

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1850007A (zh) * 2006-05-16 2006-10-25 清华大学深圳研究生院 基于心音分析的心脏病自动分类系统及其心音分段方法
CN101518439A (zh) * 2009-03-24 2009-09-02 重庆大学 一种基于心音和心电的心脏功能检测系统
CN104778342A (zh) * 2014-01-14 2015-07-15 济南大学 一种基于小波奇异熵的心音特征提取方法
CN105488520A (zh) * 2015-11-24 2016-04-13 昆明理工大学 基于多分辨奇异谱熵和svm的泄漏声发射信号识别方法
CN106226635A (zh) * 2016-07-14 2016-12-14 国网福建晋江市供电有限公司 一种配电网馈线故障类型识别方法及装置

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
D. BOUTANA ETC: "The empirical mode decomposition and the SVD for abnormal heart sound signals detection and time-frequency analysis", 《2ND IET INTERNATIONAL CONFERENCE ON INTELLIGENT SIGNAL PROCESSING 2015 (ISP)》 *
赵志宏: "基于振动信号的机械故障特征提取与诊断研究", 《中国博士学位论文全文数据库 工程科技Ⅱ辑》 *

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109934243A (zh) * 2017-12-19 2019-06-25 中国科学院深圳先进技术研究院 Ecg数据分类方法、装置、电子设备及系统
CN109567856A (zh) * 2018-11-23 2019-04-05 四川长虹电器股份有限公司 一种基于智能平台的听诊位置智能提醒方法
CN110175996A (zh) * 2019-05-29 2019-08-27 中国人民解放军陆军军医大学第二附属医院 一种基于cmr双向循环网络的心衰自动分级方法
CN116525104A (zh) * 2023-06-26 2023-08-01 中国人民解放军总医院 一种心源性休克快速分期分诊系统、设备及可存储介质
CN116525104B (zh) * 2023-06-26 2023-11-17 中国人民解放军总医院 一种心源性休克快速分期分诊系统、设备及可存储介质
CN116807486A (zh) * 2023-07-05 2023-09-29 淮安市楚州中医院 一种急性心衰的数字模型构件系统及方法
CN116807486B (zh) * 2023-07-05 2024-04-26 淮安市楚州中医院 一种急性心衰的数字模型构件系统及方法

Also Published As

Publication number Publication date
CN107480637B (zh) 2019-08-30

Similar Documents

Publication Publication Date Title
CN107480637B (zh) 基于心音特征的心衰分期方法
CN107951485B (zh) 基于人工智能自学习的动态心电图分析方法和装置
CN108714026B (zh) 基于深度卷积神经网络和在线决策融合的细粒度心电信号分类方法
CN104523266B (zh) 一种心电信号自动分类方法
CN107822622A (zh) 基于深度卷积神经网络的心电图诊断方法和系统
CN106108889B (zh) 基于深度学习算法的心电图分类方法
CN107832737A (zh) 基于人工智能的心电图干扰识别方法
CN108030488A (zh) 基于卷积神经网络的心律失常的检测系统
CN109961017A (zh) 一种基于卷积循环神经网络的心音信号分类方法
CN107714023A (zh) 基于人工智能自学习的静态心电图分析方法和装置
CN109376751A (zh) 一种基于卷积神经网络的人脑功能网络分类方法
CN107811626A (zh) 一种基于一维卷积神经网络和s变换的心律失常分类方法
CN109907752A (zh) 一种去除运动伪影干扰与心电特征检测的心电诊断与监护方法及系统
CN107958213A (zh) 一种基于脑机接口辅助康复医疗的共空间模式与深度学习方法
CN106980746A (zh) 一种基于时序分析的通用无创血糖预测方法
CN111449644A (zh) 一种基于时频变换和数据增强技术的生物电信号分类方法
CN106096632A (zh) 基于深度学习和mri图像的心室功能指标预测方法
CN108090509A (zh) 一种数据长度自适应的心电图分类方法
CN111759345B (zh) 基于卷积神经网络的心脏瓣膜异常分析方法、系统和装置
CN108304887A (zh) 基于少数类样本合成的朴素贝叶斯数据处理系统及方法
CN110148466A (zh) 一种基于迁移学习的心冲击信号房颤计算机辅助诊断方法
CN109512423A (zh) 一种基于确定学习与深度学习的心肌缺血危险分层方法
CN109222963A (zh) 一种基于卷积神经网络的心电异常识别分类方法
CN110313894A (zh) 基于卷积神经网络的心率失常分类算法
CN109948396A (zh) 一种心拍分类方法、心拍分类装置及电子设备

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