CN110207967A - 一种基于小波包能量特征与互相关的状态识别方法与系统 - Google Patents

一种基于小波包能量特征与互相关的状态识别方法与系统 Download PDF

Info

Publication number
CN110207967A
CN110207967A CN201910512776.4A CN201910512776A CN110207967A CN 110207967 A CN110207967 A CN 110207967A CN 201910512776 A CN201910512776 A CN 201910512776A CN 110207967 A CN110207967 A CN 110207967A
Authority
CN
China
Prior art keywords
cross
sample
correlation
signal
state
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
CN201910512776.4A
Other languages
English (en)
Other versions
CN110207967B (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.)
Dalian Maritime University
Original Assignee
Dalian Maritime 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 Dalian Maritime University filed Critical Dalian Maritime University
Priority to CN201910512776.4A priority Critical patent/CN110207967B/zh
Publication of CN110207967A publication Critical patent/CN110207967A/zh
Application granted granted Critical
Publication of CN110207967B publication Critical patent/CN110207967B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01HMEASUREMENT OF MECHANICAL VIBRATIONS OR ULTRASONIC, SONIC OR INFRASONIC WAVES
    • G01H17/00Measuring mechanical vibrations or ultrasonic, sonic or infrasonic waves, not provided for in the preceding groups
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01MTESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
    • G01M13/00Testing of machine parts
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01MTESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
    • G01M99/00Subject matter not provided for in other groups of this subclass

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Testing Of Devices, Machine Parts, Or Other Structures Thereof (AREA)

Abstract

本发明提供一种基于小波包能量特征与互相关的状态识别方法与系统。本发明包括:对采集到的样本信号截段进行小波包分解;计算小波包分解后各个子空间内能量值;将能量值按子空间顺序组成特征向量;基于各样本的特征向量构建样本库;将待检测元件特征向量与样本库中的样本特征向量做互相关分析,其中互相关系数最大时对应的样本状态即为所识别待测信号状态。本发明简单易行,通过软件编程将算法写入硬件平台中,即自动样本制作与状态识别,结果准确性高,识别过程所需时间短,与其他识别方法相比,自动化程度更高。

Description

一种基于小波包能量特征与互相关的状态识别方法与系统
技术领域
本发明涉及机械故障检测领域,尤其涉及一种基于小波包能量特征与互相关的状态识别方法与系统。
背景技术
大多数工程机械的执行元件是旋转机械。有些旋转机械长时间在重载条件下工作会发生机械磨损或疲劳导致机械故障,它的故障可能会导致巨大的经济损失甚至是严重的安全事故,因此对旋转机械进行故障识别的需求十分迫切。故障识别方法主要可分两部分:故障特征提取和故障特征识别。
故障特征提取分为基于物理模型和基于信号数据两大类:物理模型的方法需要相关的动态系统知识,这很难应用于复杂的实际环境中;相比之下,基于信号处理的方法因为不需要精确的物理系统模型已广泛应用于旋转机械的故障诊断中。在旋转机械故障识别领域,傅里叶变换,经验模态分解,小波包分析等是近年来最新的几种信号处理方法,但都有其局限性。其中傅里叶变换及其延伸的短时傅里叶变换很难处理机械故障信号中的信号突变部分;经验模态分解经常受边界效应和模态混淆的影响;小波包分析是小波变换的延伸,它同时对低频与高频成分进行分解。小波包能量法是将原信号分解后的每个子空间内能量值组合成向量作为故障特征,目前该方法广泛应用于旋转机械故障识别。
故障特征识别现阶段多采用支持向量机与神经网络等方法。其中支持向量机方法核函数的选取对识别结果有很大影响,当信号未知时可能无法选择适合的核函数导致结果偏差大;神经网络的弊端在于训练网络所需要的数据量非常大,而且算法上非常复杂难以实现。现有技术在原理的理解及实现上所需的知识成本与时间成本都较高,不利于实际应用的推广与使用。
于此同时,现有技术无法对小波包分解层数、进行状态识别的信号截取长度进行预先判断,只能人为手动进行特征向量的制作。为了解决传统故障特征识别方法的弊端,本发明介绍了一种基于小波包能量特征与互相关的状态识别方法,可以实现无需手动进行数据处理即可生成适合的故障特征向量与进行状态识别。
发明内容
根据上述提出的技术问题,而提供一种基于小波包能量特征与互相关的状态识别方法与系统。本发明采用的技术手段如下:
一种基于小波包能量特征与互相关的状态识别方法,包括如下步骤:
对采集到的样本信号截段进行小波包分解;
计算小波包分解后各个子空间内能量值;
将能量值按子空间顺序组成特征向量;
基于各样本的特征向量构建样本库,所述样本包括正常状态的正常元件样本和至少一种故障状态的故障元件样本;
将待检测元件特征向量与样本库中的样本特征向量做互相关分析,其中互相关系数最大时对应的样本状态即为所识别待测信号状态。
进一步地,特征向量表示如下:
其中:是元件运行状态为φ的信号截段经过小波包变换后在第j层第k个位置子空间的能量值,l为信号截段的采样点数,是第j层小波包变换的能量值按子空间顺序组成的特征向量。
进一步地,所述样本特征向量制作部分小波包层数j0选择由下式表示:
j0={j|max[D(j)]}
其中D(j)的每个元素代表小波包分解层数为j时,遍历元件所有运行状态时互相关系数差值的最小值,
特征向量TX(τ)对应的状态为X,τ代表时滞,D(j)的计算式如下所示:
其中:min(A)是变量A中的最小值,Δ为元件所有运行模式的集合。
进一步地,为了平衡识别准确率与识别时间之间的关系,合理的截段采样点数l0由下式确定:
其中f(l)是变量F(l)的拟合函数,F(l)的每个元素代表信号截段采样点数为l时,遍历元件所有运行状态时互相关系数差值的最小值,根据已确定的j0,F(l)由下式定义:
为了降低因波动导致的F(l)数值误差,采用二阶幂函数对其进行拟合处理,其基函数如下所示:
f(l)=a·lb+c。
一种基于上述小波包能量特征与互相关的状态识别方法的状态识别系统,包括:
采集单元,用于采集样本元件和待检测元件的声信号;
信号转换单元,其与采集单元相连,用于将采集到的声信号转换成数字信号;
信号处理单元,用于计算待检测元件特征向量,并将其与样本库中的样本特征向量做互相关分析,其中互相关系数最大时对应的样本状态即为所识别待测信号状态。
本发明具有以下优点:
1、简单易行。通过软件编程将算法写入硬件平台中,即自动样本制作与状态识别。2、结果准确性高。经多次测试识别准确率都在90%以上,效果显著。3、识别过程所需时间短。4.与其他识别方法相比,本方法可以实现信号截段长度的自动选取,自动化程度更高。
基于上述理由本发明可在机械故障检测领域广泛推广。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案,下面将对实施例或现有技术描述中所需要使用的附图做以简单地介绍,显而易见地,下面描述中的附图是本发明的一些实施例,对于本领域普通技术人员来讲,在不付出创造性劳动性的前提下,还可以根据这些附图获得其他的附图。
图1为本发明基于小波包能量特征与互相关的状态识别方法流程图。
图2为本发明实施例中故障元件实物图。
图3为本发明实施例中采集到的同步液压马达齿轮磨损的信号波形图。
图4为本发明实施例中采集到的同步液压马达齿轮生锈的信号波形图。
图5为本发明实施例中采集到的同步液压马达端盖磨损的信号波形图。
图6为本发明实施例中采集到的同步液压马达正常状态的信号波形图。
图7为本发明实施例中不同小波包分解层数时互相关系数的曲线图。
图8为本发明实施例中不同小波包分解层数时D(j)的曲线及j0的取值。
图9为本发明实施例中不同信号截段长度时的拟合曲线图。
图10为本发明实施例中不同信号截段长度时F(l)和f(l)的曲线。
图11为本发明实施例中不同信号截段长度时f(l)/l的曲线及l0的取值
图12为本发明实施例中状态识别处理系统框图。
图13为本发明实施例中改变外界条件时,本发明与支持向量机识别准确率对比图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚,下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例是本发明一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。
如图1所示,本发明公开了一种基于小波包能量特征与互相关的状态识别方法,包括如下步骤:
对采集到的样本信号截段进行小波包分解;
计算小波包分解后各个子空间内能量值;
将能量值按子空间顺序组成特征向量;
基于各样本的特征向量构建样本库,所述样本包括正常状态的正常元件样本和至少一种故障状态的故障元件样本;
将待检测元件特征向量与样本库中的样本特征向量做互相关分析,其中互相关系数最大时对应的样本状态即为所识别待测信号状态。
其中特征向量表示如下:
其中:是元件运行状态为φ的信号截段经过小波包变换后在第j层第k个位置子空间的能量值,l为信号截段的采样点数,是第j层小波包变换的能量值按子空间顺序组成的特征向量。
所述样本特征向量制作部分小波包层数j0选择由下式表示:
j0={j|max[D(j)]}
其中D(j)的每个元素代表小波包分解层数为j时,遍历元件所有运行状态时互相关系数差值的最小值,所述运行状态为故障类型和正常运行状态的统称。
特征向量TX(τ)对应的状态为X,τ代表时滞,D(j)的计算式如下所示:
其中:min(A)是变量A中的最小值,Δ为元件所有运行模式的集合。
随着信号截段采样点数参数l的增大,不同状态对应的特征向量差异性越大,这有利于准确度的提高但会使得识别时间增加。为了平衡识别准确率与识别时间之间的关系,合理的截段采样点数l0由下式确定:
其中f(l)是变量F(l)的拟合函数,F(l)的每个元素代表信号截段采样点数为l时,遍历元件所有运行状态时互相关系数差值的最小值,根据已确定的j0,F(l)由下式定义:
为了降低因波动导致的F(l)数值误差,采用二阶幂函数对其进行拟合处理,其基函数如下所示:
f(l)=a·lb+c。
本发明还公开了一种基于上述小波包能量特征与互相关的状态识别方法的状态识别系统,包括:
采集单元,用于采集样本元件和待检测元件的声信号;
信号转换单元,其与采集单元相连,用于将采集到的声信号转换成数字信号;
信号处理单元,用于计算待检测元件特征向量,并将其与样本库中的样本特征向量做互相关分析,其中互相关系数最大时对应的样本状态即为所识别待测信号状态。
实施例
本实施例以同步液压马达产生的故障信号为例,如图2所示,故障类型包括齿轮磨损、齿轮生锈和端盖磨损,基于小波包能量特征与互相关的状态识别方法步骤如下:
A、搭建同步液压马达信号采集系统,包括液压系统、采集卡、电脑等;
B、将算法编写程序导入单片机,完成信号处理系统内部结构的搭建;
C、将故障样本信号导入信号处理系统完成信号处理系统的初始化。
本发明步骤A中搭建声信号采集系统过程可描述为:
A1、将声传感器置于同步液压马达附近,其用于采集故障声信号;
A2、声传感器接采集卡连在电脑上,采集卡将声传感器接收的信号转换成数字信号并储存在电脑上。
本发明步骤C所述的完成信号处理系统的初始化过程可描述为:
利用小波包能量特征与互相关的识别方法中小波包层数与信号截取长度的选取算法,对二者参数进行计算并将其设置为状态识别的基础参数。
1.如图8所示,其中样本特征向量制作部分小波包层数j0选择由下式表示:
j0={j|max[D(j)]}
如图7所示,其中D(j)的每个元素代表小波包分解层数为j时,遍历元件所有运行状态时互相关系数差值的最小值。特征向量TX(τ)对应的状态为X,τ代表时滞。D(j)的计算式如下所示:
其中:min(A)是变量A中的最小值。Δ为元件所有运行模式的集合。
2.随着信号截段采样点数参数l的增大,不同状态对应的特征向量差异性越大,这有利于准确度的提高但会使得识别时间增加。如图11所示,为了平衡识别准确率与识别时间之间的关系,合理的截段采样点数l0由下式确定:
其中f(l)是变量F(l)的拟合函数,F(l)的每个元素代表信号截段采样点数为l时,遍历元件所有运行状态时互相关系数差值的最小值。其中相关系数差值的拟合曲线如图9所示。根据已确定的j0,F(l)由下式定义:
如图10所示,为了降低因波动导致的F(l)数值误差,对其进行拟合处理。经试验采用二阶幂函数拟合效果最好,其基函数如下所示:
f(l)=a·lb+c
图10中,竖坐标两个量分别是F和f(纵坐标量是互相关系数的差值,无量纲的标量,无单位)将算法编写的程序导入计算机中,利用采集卡将采集到的信号导入计算机中进行处理分析得到初始化参数。整个信号采集分析系统如图12所示。
对小波包能量互相关法的可行性进行验证,具体地,采集同步液压马达四种状态的声信号:同步液压马达齿轮磨损;同步液压马达齿轮生锈;同步液压马达端面磨损;同步液压马达正常状态,信号波形分别如图3-6所示。
同步液压马达及声信号采集处理参数如表1所示:
表1信号采集系统参数
根据本发明算法计算出适合的初始参数小波包合理层数j0与信号截取长度l0
为了验证识别准确率,在表1条件下对每种信号测试数据进行1000次识别并统计准确率。识别结果如表2所示:
表2小波包能量互相关法的识别结果
改变出口压力分别采用互相关与支持向量机方法(核函数为RBF)进行识别,统计二者准确率如图13所示。可见本发明方法在本条件下由于支持向量机法。
为了对比互相关法与SVM法(100组数据训练,核函数为RBF)的识别用时,在处理器型号为Intel Corei5-7400,RAM为8G的PC上利用matlab软件统计每种方法的计算运行时间。对长度为29,样本状态种类为4的1000组特征向量数据进行识别时间统计,1000次识别过程用时如表3所示:
表3互相关法与SVM法的识别用时
从表3中可以看出,互相关法识别用时大约为SVM法的四分之一,其识别速度比RBF内核的SVM法更快。
最后应说明的是:以上各实施例仅用以说明本发明的技术方案,而非对其限制;尽管参照前述各实施例对本发明进行了详细的说明,本领域的普通技术人员应当理解:其依然可以对前述各实施例所记载的技术方案进行修改,或者对其中部分或者全部技术特征进行等同替换;而这些修改或者替换,并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。

Claims (5)

1.一种基于小波包能量特征与互相关的状态识别方法,其特征在于,包括如下步骤:
对采集到的样本信号截段进行小波包分解;
计算小波包分解后各个子空间内能量值;
将能量值按子空间顺序组成特征向量;
基于各样本的特征向量构建样本库,所述样本包括正常状态的正常元件样本和至少一种故障状态的故障元件样本;
将待检测元件特征向量与样本库中的样本特征向量做互相关分析,其中互相关系数最大时对应的样本状态即为所识别待测信号状态。
2.根据权利要求1所述的基于小波包能量特征与互相关的状态识别方法,其特征在于,特征向量表示如下:
其中:Eφj,l(k)是元件运行状态为φ的信号截段经过小波包变换后在第j层第k个位置子空间的能量值,l为信号截段的采样点数,Tφj,l(k)是第j层小波包变换的能量值按子空间顺序组成的特征向量。
3.根据权利要求2所述的基于小波包能量特征与互相关的状态识别方法,其特征在于,所述样本特征向量制作部分小波包层数j0选择由下式表示:
j0={j|max[D(j)]}
其中D(j)的每个元素代表小波包分解层数为j时,遍历元件所有运行状态时互相关系数差值的最小值,
特征向量TX(τ)对应的状态为X,τ代表时滞,D(j)的计算式如下所示:
其中:min(A)是变量A中的最小值,Δ为元件所有运行模式的集合。
4.根据权利要求3所述的基于小波包能量特征与互相关的状态识别方法,其特征在于,为了平衡识别准确率与识别时间之间的关系,合理的截段采样点数l0由下式确定:
其中f(l)是变量F(l)的拟合函数,F(l)的每个元素代表信号截段采样点数为l时,遍历元件所有运行状态时互相关系数差值的最小值,根据已确定的j0,F(l)由下式定义:
为了降低因波动导致的F(l)数值误差,采用二阶幂函数对其进行拟合处理,其基函数如下所示:
f(l)=a·lb+c。
5.根据权利要求1~4任一项所述的基于小波包能量特征与互相关的状态识别方法的状态识别系统,包括:
采集单元,用于采集样本元件和待检测元件的声信号;
信号转换单元,其与采集单元相连,用于将采集到的声信号转换成数字信号;
信号处理单元,用于计算待检测元件特征向量,并将其与样本库中的样本特征向量做互相关分析,其中互相关系数最大时对应的样本状态即为所识别待测信号状态。
CN201910512776.4A 2019-06-13 2019-06-13 一种基于小波包能量特征与互相关的状态识别方法与系统 Expired - Fee Related CN110207967B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910512776.4A CN110207967B (zh) 2019-06-13 2019-06-13 一种基于小波包能量特征与互相关的状态识别方法与系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910512776.4A CN110207967B (zh) 2019-06-13 2019-06-13 一种基于小波包能量特征与互相关的状态识别方法与系统

Publications (2)

Publication Number Publication Date
CN110207967A true CN110207967A (zh) 2019-09-06
CN110207967B CN110207967B (zh) 2020-12-01

Family

ID=67792510

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910512776.4A Expired - Fee Related CN110207967B (zh) 2019-06-13 2019-06-13 一种基于小波包能量特征与互相关的状态识别方法与系统

Country Status (1)

Country Link
CN (1) CN110207967B (zh)

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102566555A (zh) * 2012-02-10 2012-07-11 安徽建筑工业学院 基于模式识别的白色家电工作状态监测方法
CN105022912A (zh) * 2015-05-28 2015-11-04 北京交通大学 基于小波主成分分析的滚动轴承故障预测方法
CN105424366A (zh) * 2015-12-14 2016-03-23 长春工业大学 基于eemd自适应消噪的轴承故障诊断方法
CN105740822A (zh) * 2016-02-01 2016-07-06 杭州杰牌传动科技有限公司 一种机械故障诊断优化方法及系统
CN108414219A (zh) * 2018-01-31 2018-08-17 杭州携测信息技术有限公司 一种齿轮电动机的齿轮故障诊断方法和系统
CN108444704A (zh) * 2018-03-30 2018-08-24 华中科技大学 一种滚动轴承早期故障诊断方法
CN108710889A (zh) * 2018-04-02 2018-10-26 天津大学 一种汽车发动机缺缸故障诊断方法
KR20190019713A (ko) * 2017-08-18 2019-02-27 인하대학교 산학협력단 무인 항공기 음향 식별을 위한 서포트 벡터 머신에 기반한 음향 특징 추출 및 분류 방법 그리고 시스템

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102566555A (zh) * 2012-02-10 2012-07-11 安徽建筑工业学院 基于模式识别的白色家电工作状态监测方法
CN105022912A (zh) * 2015-05-28 2015-11-04 北京交通大学 基于小波主成分分析的滚动轴承故障预测方法
CN105424366A (zh) * 2015-12-14 2016-03-23 长春工业大学 基于eemd自适应消噪的轴承故障诊断方法
CN105740822A (zh) * 2016-02-01 2016-07-06 杭州杰牌传动科技有限公司 一种机械故障诊断优化方法及系统
KR20190019713A (ko) * 2017-08-18 2019-02-27 인하대학교 산학협력단 무인 항공기 음향 식별을 위한 서포트 벡터 머신에 기반한 음향 특징 추출 및 분류 방법 그리고 시스템
CN108414219A (zh) * 2018-01-31 2018-08-17 杭州携测信息技术有限公司 一种齿轮电动机的齿轮故障诊断方法和系统
CN108444704A (zh) * 2018-03-30 2018-08-24 华中科技大学 一种滚动轴承早期故障诊断方法
CN108710889A (zh) * 2018-04-02 2018-10-26 天津大学 一种汽车发动机缺缸故障诊断方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
李刚,张宏亮,李欣,杜思琪: "基于EEMD和小波包的机车轴箱轴承故障诊断方法", 《兰州交通大学学报》 *

Also Published As

Publication number Publication date
CN110207967B (zh) 2020-12-01

Similar Documents

Publication Publication Date Title
Guo et al. Rolling bearing fault classification based on envelope spectrum and support vector machine
CN107013449B (zh) 基于深度学习的声音信号识别压缩机故障的方法及系统
CN111914883B (zh) 一种基于深度融合网络的主轴轴承状态评估方法及装置
CN108426713A (zh) 基于小波变换和深度学习的滚动轴承微弱故障诊断方法
CN110657985B (zh) 基于奇异值谱流形分析的齿轮箱故障诊断方法及系统
CN104697767B (zh) 一种基于振动分析的转子系统故障诊断方法及装置
CN103900824B (zh) 基于瞬时转速聚类分析的柴油机故障诊断方法
CN111504645B (zh) 一种基于频域多点峭度的滚动轴承故障诊断方法
CN110487547B (zh) 基于振动图和迁移学习的变工况下滚动轴承故障诊断方法
CN106441896A (zh) 滚动轴承故障模式识别及状态监测的特征向量提取方法
CN110261116A (zh) 一种轴承故障检测方法及装置
CN107963239A (zh) 一种基于音频的运载火箭故障检测装置及检测方法
CN112461934B (zh) 一种基于声发射的航空发动机叶片裂纹源定位方法
CN105678343A (zh) 基于自适应加权组稀疏表达的水电机组噪声异常诊断方法
CN108388692A (zh) 基于分层稀疏编码的滚动轴承故障特征提取方法
CN112304611A (zh) 一种基于深度学习的轴承故障的诊断方法
CN109605128B (zh) 一种基于功率谱熵差的铣削颤振在线检测方法
CN114755017B (zh) 一种跨域数据驱动无监督领域共享网络的变转速轴承故障诊断方法
CN111898443A (zh) 一种fdm型3d打印机送丝机构流量监测方法
CN108153987A (zh) 一种基于超限学习机的液压泵多故障诊断方法
CN111898644A (zh) 一种无故障样本下航天液体发动机健康状态智能识别方法
CN106708009A (zh) 一种基于支持向量机聚类的船舶动力定位测量系统多故障诊断方法
CN111665050A (zh) 一种基于聚类k-svd算法的滚动轴承故障诊断方法
CN114548295A (zh) 基于多尺度领域自适应网络的轴承故障分类系统及方法
CN114263621A (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
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20201201

CF01 Termination of patent right due to non-payment of annual fee