CN105928701A - 基于相关分析的emd过程中有效imf的判定方法 - Google Patents

基于相关分析的emd过程中有效imf的判定方法 Download PDF

Info

Publication number
CN105928701A
CN105928701A CN201610279173.0A CN201610279173A CN105928701A CN 105928701 A CN105928701 A CN 105928701A CN 201610279173 A CN201610279173 A CN 201610279173A CN 105928701 A CN105928701 A CN 105928701A
Authority
CN
China
Prior art keywords
imf
component
correlation coefficient
emd
pseudo
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
CN201610279173.0A
Other languages
English (en)
Other versions
CN105928701B (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.)
Xi'an Triumph Electronic Technology Co., Ltd.
Shijiazhuang Tiedao University
Original Assignee
Shijiazhuang Tiedao 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 Shijiazhuang Tiedao University filed Critical Shijiazhuang Tiedao University
Priority to CN201610279173.0A priority Critical patent/CN105928701B/zh
Publication of CN105928701A publication Critical patent/CN105928701A/zh
Application granted granted Critical
Publication of CN105928701B publication Critical patent/CN105928701B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • 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
    • G01M13/02Gearings; Transmission mechanisms
    • G01M13/028Acoustic or vibration analysis

Landscapes

  • Physics & Mathematics (AREA)
  • Acoustics & Sound (AREA)
  • General Physics & Mathematics (AREA)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)

Abstract

本发明公开了一种基于相关分析的EMD过程中有效IMF的判定方法,涉及专门适用于特定应用的数字计算或数据处理的设备或方法技术领域。所述方法包括如下步骤:采集故障振动信号,将该振动信号作为原始信号;将原始信号使用EMD方法进行分解得到若干个IMF分量;计算每一个IMF分量与原始信号的相关系数,画出相关系数分布图;利用噪声判断准则确定噪声IMF分量;利用伪分量判断准则确定伪IMF分量;将判定为噪声的IMF分量和伪IMF分量去除,得到有效的IMF分量。所述方法有效的解决了EMD过程产生的IMF分量中的噪声分量和伪分量的判定问题,进而得到有效分量,提高EMD的分析效果。

Description

基于相关分析的EMD过程中有效IMF的判定方法
技术领域
本发明涉及专门适用于特定应用的数字计算或数据处理的设备或方法技术领域,尤其涉及一种基于相关分析的EMD过程中有效IMF的判定方法。
背景技术
经验模态分解(Empirical Mode Decomposition,EMD)方法由Huang在1998年提出,至今已在多个领域得到了十分广泛的应用。信号经过EMD的“筛分”过程,将信号进行分解得到若干个不同频率的本征模态分量(IMF),频率从高到低排列,其频率成分和带宽随信号的不同而不同,体现了该方法的自适应性,因此EMD方法可以作为一个自适应滤波器。由于EMD方法的这一特点,我们可以根据所要研究信号的特点和要求,选择性地将适当的IMF相加,进行信号重构,把我们感兴趣的某一频率范围内的信号特征凸显出来,从而达到滤波器滤波的效果。信号经EMD分解后得到有限个频率从高到低的IMF分量,其中分解得到的前几个阶数较小的IMF对应于信号的高频成分,通常认为它们主要包含了原信号中的噪声成分;分解得到的后几个阶数较大的IMF对应于信号的低频成分,通常认为它们包含的噪声成分很少。
EMD去噪方法的主要思想是把原信号经过层层分解成为不同特征时间尺度的IMF分量,有选择性地将某些IMF分量重新组合起来,得到一个新的信号,就可以组成低通滤波器、高通滤波器、带通滤波器或带阻滤波器,从而达到去噪的目的。大多数基于EMD的去噪方法都是将IMF分量中阶次较低的高频分量作为噪声直接去除,但是具体应该去除几个、去除哪几个分量目前还没有特定的准则,很多情况下都要依靠人的经验,不能排除人的主观因素对其造成的影响。
发明内容
本发明所要解决的技术问题是提供一种基于相关分析的EMD过程中有效IMF的判定方法,所述方法有效的解决了EMD过程产生的IMF分量中的噪声分量和伪分量的判定问题,进而得到有效分量,提高EMD的分析效果。
为解决上述技术问题,本发明所采取的技术方案是:一种基于相关分析的EMD过程中有效IMF的判定方法,其特征在于所述方法包括如下步骤:
采集故障振动信号,将该振动信号作为原始信号;
将原始信号使用EMD方法进行分解得到若干个IMF分量;
计算每一个IMF分量与原始信号的相关系数,画出相关系数分布图;
根据相关系数分布图并利用噪声判断准则确定噪声IMF分量;
根据相关系数分布图并利用伪分量判断准则确定伪IMF分量;
将判定为噪声的IMF分量和伪IMF分量去除,得到有效的IMF分量。
进一步的技术方案在于:采用压电加速度传感器采集故障振动信号。
进一步的技术方案在于:假设原始信号经EMD方法分解得到N个频率从高到低排列的IMF分量,噪声成分主要集中在信号的高频段,即EMD方法分解得到的阶数较小的IMF当中,因此存在一个数字k,使得IMF1~IMFk是以噪声为主导模态的IMF分量;伪分量主要集中在阶数较大的IMF分量中,因此存在一个数字h,使得从第h个IMF往后的所有分量IMFh~IMFN都是伪IMF分量;剩下的IMFk+1~IMFh-1是以有用信号为主导模态的分量。
进一步的技术方案在于:所述的计算每一个IMF分量与原始信号的相关系数的方法包括如下步骤:
相关系数ρxy的计算公式为:
ρ x y = Σ n = 0 ∞ x ( n ) y ( n ) [ Σ n = 0 ∞ x 2 ( n ) Σ n = 0 ∞ y 2 ( n ) ] 1 / 2
式中x(n)表示EMD分解前的原始信号,y(n)表示IMF分量。
进一步的技术方案在于:所述的根据相关系数分布图并利用噪声判断准则确定噪声IMF分量的方法包括如下步骤:
计算各个IMF分量与原信号的相关系数,画出相关系数分布图,在所述相关系数分布图中找到第一个相关系数发生突变的点K,该点对应的IMF分量为IMFk,将IMF1~IMFk分量视为噪声IMF分量。
进一步的技术方案在于:将各个IMF分量与原信号的相关系数点依次连接构成相关系数曲线,曲线方向发生反转的相关系数对应的点为相关系数发生突变的点。
进一步的技术方案在于:所述的根据相关系数分布图并利用伪分量判断准则确定伪IMF分量的方法包括如下步骤:
计算各个IMF分量与原信号的相关系数,画出相关系数分布图,从图中找到第一个相关系数小于A的点h,该点对应的IMF分量为IMFh,将IMFh+1~IMFN分量视为伪分量。
进一步的技术方案在于:所述A的值为0.1。
采用上述技术方案所产生的有益效果在于:假设原始信号经EMD方法分解得到N个频率从高到低排列的IMF分量,噪声成分主要集中在信号的高频段,即EMD方法分解得到的阶数较小的IMF当中,因此存在一个数字k,使得IMF1~IMFk是以噪声为主导模态的IMF分量;伪分量主要集中在阶数较大的IMF分量中,因此存在一个数字h,使得从第h个IMF往后的所有分量IMFh~IMFN都是伪IMF分量;剩下的IMFk+1~IMFh-1是以有用信号为主导模态的分量。因此,所述方法有效的解决了EMD过程产生的IMF分量中的噪声分量和伪分量的判定问题,进而得到有效分量,提高EMD的分析效果。
附图说明
图1是本发明所述方法的流程图;
图2齿轮磨损故障振动信号时域图;
图3a-3b齿轮磨损故障振动信号EMD分解得到的若干个IMF图;
图4是若干个IMF分量的相关系数分布图;
图5重构信号的时域图。
具体实施方式
下面结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例仅仅是本发明的一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。
在下面的描述中阐述了很多具体细节以便于充分理解本发明,但是本发明还可以采用其他不同于在此描述的其它方式来实施,本领域技术人员可以在不违背本发明内涵的情况下做类似推广,因此本发明不受下面公开的具体实施例的限制。
如图1所示,本发明公开了一种基于相关分析的EMD过程中有效IMF的判定方法,包括如下步骤:
采用压电加速度传感器采集故障振动信号,将该振动信号作为原始信号;
将原始信号使用EMD方法进行分解得到若干个IMF分量。
计算每一个IMF分量与原始信号的相关系数,画出相关系数分布图,相关系数ρxy的计算公式为:
ρ x y = Σ n = 0 ∞ x ( n ) y ( n ) [ Σ n = 0 ∞ x 2 ( n ) Σ n = 0 ∞ y 2 ( n ) ] 1 / 2
式中x(n)表示EMD分解前的原始信号,y(n)表示IMF分量。
根据相关系数分布图并利用噪声判断准则确定噪声IMF分量,确定噪声IMF分量的方法如下:计算各个IMF分量与原信号的相关系数,画出相关系数分布图,在所述相关系数分布图中找到第一个相关系数发生突变的点K,该点对应的IMF分量为IMFk,将IMF1~IMFk分量视为噪声IMF分量。将各个IMF分量与原信号的相关系数点依次连接构成相关系数曲线,曲线方向发生反转的相关系数对应的点为相关系数发生突变的点。
根据相关系数分布图并利用伪分量判断准则确定伪IMF分量,确定伪IMF分量的方法如下:计算各个IMF分量与原信号的相关系数,画出相关系数分布图,从图中找到第一个相关系数小于0.1的点h,该点对应的IMF分量为IMFh,将IMFh+1~IMFN分量视为伪分量。
将判定为噪声的IMF分量和伪IMF分量去除,得到有效的IMF分量。
设原始信号经EMD方法分解得到N个频率从高到低排列的IMF分量,噪声成分主要集中在信号的高频段,即EMD方法分解得到的阶数较小的IMF当中,因此存在一个数字k,使得IMF1~IMFk是以噪声为主导模态的IMF分量;伪分量主要集中在阶数较大的IMF分量中,因此存在一个数字h,使得从第h个IMF往后的所有分量IMFh~IMFN都是伪IMF分量;剩下的IMFk+1~IMFh-1是以有用信号为主导模态的分量。因此,所述方法有效的解决了EMD过程产生的IMF分量中的噪声分量和伪分量的判定问题,进而得到有效分量,提高EMD的分析效果。
实验分析:
采集齿轮磨损故障振动信号,信号时域波形如图2所示。将故障信号进行EMD分解,得到12个IMF和一个趋势项,趋势项对本发明没有影响,故将其忽略不计,各个IMF的时域图如图3a-3b所示。
计算每一个IMF与齿轮磨损故障振动信号的相关系数,如下表所示:
表1 齿轮磨损信号EMD分解得到的IMFs的相关系数
画出各个IMF的相关系数分布图,如图4所示。相关系数分布突变的第一个分量是IMF3,根据噪声分量判断准则,噪声和有效分量的分界点是k=3,则IMF1~IMF3就是以噪声为主导模态的分量;根据伪分量判断准则,相关系数小于0.1的第一个分量是IMF8,有效分量和伪分量的分界点是h=8,则IMF8~IMF12就是伪分量。将以噪声为主导模态的分量和伪分量全部去除,剩下的所有分量即为有效的IMF,将其相加得到重构信号。重构信号时域图如图5所示。与图3相比,峰值之间具有明显的周期性,计算得出周期为16.67Hz,与齿轮磨损故障频率16.97Hz很接近。因此,当本方法应用到齿轮故障判断中,可以判断出该齿轮是否发生故障。

Claims (8)

1.一种基于相关分析的EMD过程中有效IMF的判定方法,其特征在于所述方法包括如下步骤:
采集故障振动信号,将该振动信号作为原始信号;
将原始信号使用EMD方法进行分解得到若干个IMF分量;
计算每一个IMF分量与原始信号的相关系数,画出相关系数分布图;
根据相关系数分布图并利用噪声判断准则确定噪声IMF分量;
根据相关系数分布图并利用伪分量判断准则确定伪IMF分量;
将判定为噪声的IMF分量和伪IMF分量去除,得到有效的IMF分量。
2.如权利要求1所述的基于相关分析的EMD过程中有效IMF的判定方法,其特征在于:
采用压电加速度传感器采集故障振动信号。
3.如权利要求1所述的基于相关分析的EMD过程中有效IMF的判定方法,其特征在于:
假设原始信号经EMD方法分解得到N个频率从高到低排列的IMF分量,噪声成分主要集中在信号的高频段,即EMD方法分解得到的阶数较小的IMF当中,因此存在一个数字k,使得IMF1~IMFk是以噪声为主导模态的IMF分量;伪分量主要集中在阶数较大的IMF分量中,因此存在一个数字h,使得从第h个IMF往后的所有分量IMFh~IMFN都是伪IMF分量;剩下的IMFk+1~IMFh-1是以有用信号为主导模态的分量。
4.如权利要求1所述的基于相关分析的EMD过程中有效IMF的判定方法,其特征在于,所述的计算每一个IMF分量与原始信号的相关系数的方法包括如下步骤:
相关系数ρxy的计算公式为:
ρ x y = Σ n = 0 ∞ x ( n ) y ( n ) [ Σ n = 0 ∞ x 2 ( n ) Σ n = 0 ∞ y 2 ( n ) ] 1 / 2
式中x(n)表示EMD分解前的原始信号,y(n)表示IMF分量。
5.如权利要求1所述的基于相关分析的EMD过程中有效IMF的判定方法,其特征在于,所述的根据相关系数分布图并利用噪声判断准则确定噪声IMF分量的方法包括如下步骤:
计算各个IMF分量与原信号的相关系数,画出相关系数分布图,在所述相关系数分布图中找到第一个相关系数发生突变的点K,该点对应的IMF分量为IMFk,将IMF1~IMFk分量视为噪声IMF分量。
6.如权利要求5所述的基于相关分析的EMD过程中有效IMF的判定方法,其特征在于:将各个IMF分量与原信号的相关系数点依次连接构成相关系数曲线,曲线方向发生反转的相关系数对应的点为相关系数发生突变的点。
7.如权利要求1所述的基于相关分析的EMD过程中有效IMF的判定方法,其特征在于,所述的根据相关系数分布图并利用伪分量判断准则确定伪IMF分量的方法包括如下步骤:
计算各个IMF分量与原信号的相关系数,画出相关系数分布图,从图中找到第一个相关系数小于A的点h,该点对应的IMF分量为IMFh,将IMFh+1~IMFN分量视为伪分量。
8.如权利要求7所述的基于相关分析的EMD过程中有效IMF的判定方法,其特征在于,所述A的值为0.1。
CN201610279173.0A 2016-04-29 2016-04-29 基于相关分析的emd过程中有效imf的判定方法 Active CN105928701B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201610279173.0A CN105928701B (zh) 2016-04-29 2016-04-29 基于相关分析的emd过程中有效imf的判定方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201610279173.0A CN105928701B (zh) 2016-04-29 2016-04-29 基于相关分析的emd过程中有效imf的判定方法

Publications (2)

Publication Number Publication Date
CN105928701A true CN105928701A (zh) 2016-09-07
CN105928701B CN105928701B (zh) 2017-07-18

Family

ID=56837666

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201610279173.0A Active CN105928701B (zh) 2016-04-29 2016-04-29 基于相关分析的emd过程中有效imf的判定方法

Country Status (1)

Country Link
CN (1) CN105928701B (zh)

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106568589A (zh) * 2016-11-04 2017-04-19 东南大学 一种基于经验小波变换碰摩声发射消噪方法
CN107392867A (zh) * 2017-07-14 2017-11-24 中国科学院遥感与数字地球研究所 基于emd的sar图像自适应辐射均衡方法
CN108387373A (zh) * 2017-12-06 2018-08-10 上海电力学院 基于相关系数改进变分模态分解的滚动轴承故障诊断方法
CN108427067A (zh) * 2018-06-12 2018-08-21 国网江苏省电力有限公司宜兴市供电分公司 一种开关柜局部放电故障检测方法、装置及系统
CN112378414A (zh) * 2020-11-20 2021-02-19 深圳信息职业技术学院 一种基于pm2.5健康出行的路径规划装置及方法
CN113616213A (zh) * 2021-07-29 2021-11-09 山东大学 一种基于bp神经网络及改进的emd方法的心电信号去噪方法、设备及存储介质
CN114112013A (zh) * 2021-11-04 2022-03-01 北京建筑大学 古建筑安全性测定方法和装置、电子设备和存储介质
CN116341993A (zh) * 2023-05-29 2023-06-27 无锡兴达泡塑新材料股份有限公司 一种用于聚苯乙烯生产过程中状态监测方法及系统

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102620928A (zh) * 2012-03-02 2012-08-01 燕山大学 基于小波半软阈值和emd的风电齿轮箱故障诊断方法
CN104006961A (zh) * 2014-04-29 2014-08-27 北京工业大学 基于经验模态分解与倒频谱的摆线锥齿轮故障诊断方法
US20150066390A1 (en) * 2013-08-30 2015-03-05 National Central University Error measuring method of gear
CN104748961A (zh) * 2015-03-30 2015-07-01 中国矿业大学 基于svd分解降噪和相关性eemd熵特征的齿轮故障诊断方法
CN105046025A (zh) * 2015-08-31 2015-11-11 浙江大学 一种核磁共振多相流测量中各相分离的方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102620928A (zh) * 2012-03-02 2012-08-01 燕山大学 基于小波半软阈值和emd的风电齿轮箱故障诊断方法
US20150066390A1 (en) * 2013-08-30 2015-03-05 National Central University Error measuring method of gear
CN104006961A (zh) * 2014-04-29 2014-08-27 北京工业大学 基于经验模态分解与倒频谱的摆线锥齿轮故障诊断方法
CN104748961A (zh) * 2015-03-30 2015-07-01 中国矿业大学 基于svd分解降噪和相关性eemd熵特征的齿轮故障诊断方法
CN105046025A (zh) * 2015-08-31 2015-11-11 浙江大学 一种核磁共振多相流测量中各相分离的方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
林丽等人: "基于相关系数的EMD 改进算法", 《计算机与数学工程》 *

Cited By (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106568589A (zh) * 2016-11-04 2017-04-19 东南大学 一种基于经验小波变换碰摩声发射消噪方法
CN107392867A (zh) * 2017-07-14 2017-11-24 中国科学院遥感与数字地球研究所 基于emd的sar图像自适应辐射均衡方法
CN108387373A (zh) * 2017-12-06 2018-08-10 上海电力学院 基于相关系数改进变分模态分解的滚动轴承故障诊断方法
CN108427067A (zh) * 2018-06-12 2018-08-21 国网江苏省电力有限公司宜兴市供电分公司 一种开关柜局部放电故障检测方法、装置及系统
CN112378414A (zh) * 2020-11-20 2021-02-19 深圳信息职业技术学院 一种基于pm2.5健康出行的路径规划装置及方法
CN113616213A (zh) * 2021-07-29 2021-11-09 山东大学 一种基于bp神经网络及改进的emd方法的心电信号去噪方法、设备及存储介质
CN113616213B (zh) * 2021-07-29 2023-02-28 山东大学 一种基于bp神经网络及改进的emd方法的心电信号去噪方法、设备及存储介质
CN114112013A (zh) * 2021-11-04 2022-03-01 北京建筑大学 古建筑安全性测定方法和装置、电子设备和存储介质
CN116341993A (zh) * 2023-05-29 2023-06-27 无锡兴达泡塑新材料股份有限公司 一种用于聚苯乙烯生产过程中状态监测方法及系统
CN116341993B (zh) * 2023-05-29 2023-07-25 无锡兴达泡塑新材料股份有限公司 一种用于聚苯乙烯生产过程中状态监测方法及系统

Also Published As

Publication number Publication date
CN105928701B (zh) 2017-07-18

Similar Documents

Publication Publication Date Title
CN105928701A (zh) 基于相关分析的emd过程中有效imf的判定方法
CN108845250B (zh) 基于振动信号特征提取的有载分接开关故障识别方法
CN108801630B (zh) 单通道盲源分离的齿轮故障诊断方法
Dong et al. Non-iterative denoising algorithm for mechanical vibration signal using spectral graph wavelet transform and detrended fluctuation analysis
CN102697495B (zh) 基于总体平均经验模式分解的二代小波肌电信号消噪方法
CN103961092B (zh) 基于自适应阈值处理的脑电信号去噪方法
CN103984866A (zh) 一种基于局域均值分解的信号去噪方法
CN105411565A (zh) 基于广义尺度小波熵的心率变异性特征分类方法
CN105320969A (zh) 基于多尺度Renyi熵的心率变异性特征分类方法
CN111007316B (zh) 基于fft与dwt的混合谐波检测改进方法
CN105677035A (zh) 基于eemd和小波阈值的运动想象脑电信号消噪方法
CN101259016A (zh) 实时自动检测癫痫特征波的方法
CN109871733A (zh) 一种自适应海杂波信号去噪方法
CN104391336B (zh) 用于地球天然脉冲电磁场数据处理的时频谱分析方法
CN102855408A (zh) 基于ica的改进emd过程中imf判定方法
CN112957056B (zh) 利用协同网络的肌肉疲劳等级特征的提取方法及系统
CN110151175A (zh) 基于ceemd与改进小波阈值的表面肌电信号消噪方法
CN103761424A (zh) 基于二代小波和ica的肌电信号降噪与去混迭方法
CN105573104A (zh) 基于改进emd的手表检测降噪方法
CN114358093A (zh) 一种电力设备中局部放电的检测方法及设备
CN115730199A (zh) 一种滚动轴承振动信号降噪和故障特征提取方法和系统
CN106706122A (zh) 基于相关系数和emd滤波特性的碰摩声发射信号降噪方法
CN109901224A (zh) 一种地震资料低频信号保护压制噪声方法
CN116610907B (zh) 基于变分模态分解的齿轮振动信号特征提取方法
CN108020761B (zh) 一种局部放电去噪方法

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
CB03 Change of inventor or designer information
CB03 Change of inventor or designer information

Inventor after: Hao Rujiang

Inventor after: Wang Shengjiang

Inventor after: Li Daiyong

Inventor after: Shen Yingming

Inventor after: Niu Zhilei

Inventor after: Yang Hongna

Inventor after: Li Hui

Inventor before: Hao Rujiang

Inventor before: Wang Shengjiang

Inventor before: Niu Zhilei

Inventor before: Yang Hongna

Inventor before: Li Hui

TA01 Transfer of patent application right
TA01 Transfer of patent application right

Effective date of registration: 20170621

Address after: 050043 Hebei province Shijiazhuang North Second Ring Road No. 17

Applicant after: Shijiazhuang Tiedao University

Applicant after: Xi'an Triumph Electronic Technology Co., Ltd.

Address before: 050011 Hebei province Shijiazhuang North Second Ring Road No. 17

Applicant before: Shijiazhuang Tiedao University

GR01 Patent grant
GR01 Patent grant