CN107036808B - Gearbox of wind turbine combined failure diagnostic method based on support vector machines probability Estimation - Google Patents

Gearbox of wind turbine combined failure diagnostic method based on support vector machines probability Estimation Download PDF

Info

Publication number
CN107036808B
CN107036808B CN201710232999.6A CN201710232999A CN107036808B CN 107036808 B CN107036808 B CN 107036808B CN 201710232999 A CN201710232999 A CN 201710232999A CN 107036808 B CN107036808 B CN 107036808B
Authority
CN
China
Prior art keywords
sample
sensor
probability
support vector
vector machines
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.)
Active
Application number
CN201710232999.6A
Other languages
Chinese (zh)
Other versions
CN107036808A (en
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.)
Zhejiang University ZJU
Original Assignee
Zhejiang University ZJU
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 Zhejiang University ZJU filed Critical Zhejiang University ZJU
Priority to CN201710232999.6A priority Critical patent/CN107036808B/en
Publication of CN107036808A publication Critical patent/CN107036808A/en
Application granted granted Critical
Publication of CN107036808B publication Critical patent/CN107036808B/en
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

Abstract

The invention discloses a kind of gearbox of wind turbine combined failure diagnostic method based on support vector machines probability Estimation.The invention mainly comprises following two modules: first module is vibration signal processing and fault signature extraction algorithm based on overall experience mode decomposition (EEMD);Second module is the training pattern method and probability Estimation algorithm of support vector machines.The present invention is respectively established for different sensors, and when analysis is integrated the result of all classifiers, improves the accuracy of diagnosis.The present invention gives detailed algorithm description to emulate signal as test object, and passes through experimental verification validity of the algorithm in terms of gear-box combined failure diagnosis.

Description

Gearbox of wind turbine combined failure diagnosis based on support vector machines probability Estimation Method
Technical field
The invention belongs to fault signature extraction and method for diagnosing faults field, more particularly to one kind are general based on support vector machines The fault signature for wind-driven generator group wheel box of rate estimation extracts and method for diagnosing faults.
Background technique
In the past few years, with the fast development of Wind Power Generation Industry, a large amount of Wind turbines are constantly built by deployment, wind Electric single-machine capacity is also being continuously increased, and the accident because of caused by wind power generating set failure happens occasionally, and causes huge warp Ji loss, hinders the development speed of Wind Power Generation Industry.Gear-box is located in engine rooms of wind power generators, is that wind-driven generator transmission is dynamic The main component of power and the important pivot of connection main shaft and generator, its operation can be interfered by many factors, such as wind speed wave Dynamic and load dynamic change.Simultaneously as environment locating for Wind turbines is generally relatively more severe, fluctuations in wind speed and load variation Acutely and frequently, therefore gearbox parts easily Aging Damage during long-play, all kinds of failures are generated.Therefore, it grinds Study carefully a set of feasible Fault Diagnosis of Gear Case method to be of great significance for the stable operation of Wind turbines.
In order to monitor the operating status of gear case of blower, it generally will be installed several vibrating sensors outside gear-box, pass The vibration signal measured is sent to host computer by sensor, is analyzed using certain diagnostic method vibration signal and is speculated at this time The operating status of gear-box.With the continuous accumulation of gear case of blower vibration data, several machine learning algorithms are also introduced into To Fault Diagnosis of Gear Case field.Support vector machines is a kind of machine learning algorithm proposed on the basis of Statistical Learning Theory, Suitable for the troubleshooting issue under Small Sample Size, it is widely used in mechanical breakdown field.The essence of algorithm of support vector machine Disaggregated model exactly is established on the basis of existing gear-box vibration data, when new sample data temporarily, to enter data into Disaggregated model, to judge mechanical operating status.However, the support vector machines that scholars propose at present still has several offices It is sex-limited: firstly, the output result of existing algorithm of support vector machine is only category label, when support vector machines is judged by accident When, which can not provide more information with regard to entirely ineffective.Secondly, current algorithm of support vector machine solution is most It is single failure problem, i.e., the failure that sample occurs is preferably at most one, less to combined failure diagnosis research.Finally, true Gearbox fault signal has the characteristics that non-stationary, non-linear, existing method showed when extracting the feature of such signal it is poor, Influence the final accuracy rate of diagnosis of support vector machines.
The present invention is directed to the shortcoming of existing method, provide complete set based on support vector machines probability Estimation Method for diagnosing faults, algorithm mainly include following two module: first module is based on overall experience mode decomposition (EEMD) Vibration signal processing and fault signature extraction algorithm;Second module is that the training pattern method of support vector machines and probability are estimated Calculating method.The novelty of the diagnostic method is mainly reflected in the following aspects: being good at processing using overall experience mode decomposition The characteristics of non-linear, non-stationary signal, the vibration signal there are combined failure is decoupled, the feature of different faults is separated. Meanwhile the energy ratio of faults energy and the approximate entropy of faults confusion degree are chosen when extracting fault signature, optimization event Hinder Feature Selection.In addition, this method is respectively established for different sensors in Training Support Vector Machines model, point The result of these classifiers is integrated when analysis, improves the accuracy of diagnosis.
Summary of the invention
In order to improve existing method, the diagnosis of gear-box combined failure is able to solve the purpose of the present invention is to provide one kind and is asked The method of topic, and improve the accuracy of diagnosis.
The purpose of the invention is achieved by the following technical solution, comprising the following steps:
Step (1), from the vibration signal extracted in gear-box monitoring data library under different working condition as sample.Assuming that M sensor is installed on gear-box, according to the corresponding sample of each sensor of the sequential processes of sensor, initializes m=1, To m-th of sensor, the fault signature of its respective sample is sought in the steps below.
Step (2) carries out EEMD decomposition to the sample point x (t) in sample, and EEMD decomposable process is as follows:
(2-1) initializes the number of iterations N=1, and predetermined maximum number of iterations is K;
The white Gaussian noise w (t) of computer random generation is added as new x (t) in (2-2) in x (t), calculates x (t) Upper lower envelope average value m1
(2-3) calculates x (t) and m1Difference h10, judge h10It whether is intrinsic mode function, if it is, enabling eigen mode Measure c1=h10, and it is transferred to (2-5), if it is not, then being transferred to (2-4);
(2-4) calculates k h repeatedly1(k-1)-m1k=h1k, until h1kMeet intrinsic mode function condition, enables c1=h1k;Its In, m1kFor h1(k-1)Upper lower envelope average value;
(2-5) separates c from x (t)1, enable r1=x (t)-c1, by r1Regard new x (t) as, constantly repeats (2-2)-(2- 4) it, is often repeated once, so that it may obtain new intrinsic modulus, i.e. r1-c2=r2,...,rn-1-cn=rn, stop standard until meeting Then;
The intrinsic modulus and residual volume that (2-6) record epicycle is decomposed, it is assumed that isolate q intrinsic modulus, then enable c1N =c1,c2N=c2,...,cqN=cq,rqN=rq, and the number of iterations N is added 1, if N is less than predetermined maximum number of iterations K, turn Enter (2-2), otherwise calculating final intrinsic modulus and residual volume, calculation formula isWith
Step (3) calculates the energy ratio and approximate entropy of each intrinsic modulus and residual volume, constitutes m-th of sensor and corresponds to sample This feature vector.
Step (4) using m-th sensor respective sample feature vector construction optimal decision function to training sample into Row classification, categorised decision function are as follows:
F (x)=sgn (wTx+b)
Wherein w and b is respectively the weight coefficient vector sum deviation of the categorised decision function, and x is sampling feature vectors, this Two parameters can determine by solving objective function, objective function are as follows:
s.t.yi(wTx+b)-1+ξi≥0
Wherein Q indicates the number of samples in m-th of sensor, ξiFor slack variable, C is penalty coefficient, yiFor sample mark Label, the classifier constructed at this time are denoted as SVMm.
Probability between the class of step (5) calculating sampleIts Middle i, j indicate i-th kind of operating status, jth kind operating status, and f is the value that sample point inputs that decision function obtains, and A and B are to pass through Following objective function is minimized to obtain
Sample point x (t) belongs to the Probability p of i-th kind of operating status under the model for step (6) estimationi=P (y=i | x), Calculation method is to optimize following objective function
Step (7), to next sensor repeat step (2)-step (6), even m adds one, establish it is next support to Amount machine classifier;Until all the sensors analysis finishes, indicates that fault diagnosis model is established and complete;
Test sample input model is carried out fault diagnosis by step (8).Specifically, when test sample data are acquired from multiple Different sensors, then need to belong in sample the vibration data separation of different sensors, and sequentially inputs and had built up The output probability of all models is finally averaging, obtains final result by the supporting vector machine model of respective sensor.
For the present invention compared with existing method for diagnosing faults, beneficial effect includes: of the invention by EEMD algorithm introducing gear-box The characteristics of fault signature of vibration signal extracts, and agrees with gear-box vibration signal time-varying, non-linear, non-stationary, and can be effective Solve the phenomenon that combined failure is mutually coupled;The sample fault signature that the present invention extracts is energy ratio and approximate entropy, is demonstrate,proved through experiment Real diagnosis with higher;The present invention belongs to the probability of different conditions using support vector machines output sample, is able to solve The shortcomings that traditional support vector machine model can not provide fault message when judging by accident;The method of the present invention is for different sensors point Do not establish corresponding diagnostic model, when diagnosis is again integrated the result of each diagnostic model, effectively improves the spirit of algorithm Activity and stability.
Detailed description of the invention
Fig. 1 algorithm frame schematic diagram;
Fig. 2 simulation sample point time domain waveform;
Intrinsic modulus and residual time domain waveform of 1 vibration signal of Fig. 3 simulation sample point sensor after EEMD algorithm process Figure.
Specific embodiment
Specific implementation method of the invention is described further with reference to the accompanying drawing:
Fig. 1 is total algorithm block schematic illustration of the invention, as shown in the figure, it is assumed that sample signal is acquired from gear-box M sensor, then sample x can be decomposed into x according to sensor1,x2,...,xM, vibration performance f is extracted respectively1,f2,...,fM, and Establish different supporting vector machine models.When new sample input, the diagnostic result of comprehensive all supporting vector machine models is obtained Last conclusion out.
In order to verify effectiveness of the invention, several groups emulation signal is constructed by Matlab to simulate gear-box vibration Signal.Emulation experiment setting is as follows: number of sensors M=2, can use sensor 1 and sensor 2, while also illustrating that and needing to instruct Practice two supporting vector machine models.Fig. 2 is a simulation sample point time domain waveform randomly selected, it can be seen that it has two-way Signal respectively indicates the signal that two sensors measure, x1For the signal that sensor 1 is surveyed, x2The signal surveyed for sensor 2.Failure Type set is two kinds, is indicated with failure 1 and failure 2.Simulation sample can be divided into four classes according to respective state: normal sample, 1 sample of failure, 2 sample of failure and combined failure (failure 1 and failure 2 exist simultaneously) sample.Training sample be normal sample, Three kinds of 1 sample of failure and 2 sample of failure, combined failure sample is as test sample.Through discussed above it is found that proposed by the present invention Supporting vector machine model is generated by normal sample and single failure sample training, but can be used to diagnose combined failure sample.
Specific step is as follows for execution:
Step (1) generates several normal samples, 1 sample of failure and 2 sample of failure as training sample using Matlab. As it is assumed that being equipped with 2 sensors on gear-box, then all there is two-way sensor signal in all sample points, by the suitable of sensor Each corresponding sample of sensor of sequence processing, initializes m=1.
Step (2) corresponds to the sample point x (t) in sample to m-th of sensor and carries out EEMD decomposition, the road Tu3Wei Mou sensing The intrinsic modulus and surplus time domain waveform that device signal is decomposed through EEMD.The specific decomposable process of EEMD is as follows:
(2-1) initializes the number of iterations N=1, and predetermined maximum number of iterations is K;
The white Gaussian noise w (t) of computer random generation is added as new x (t) in (2-2) in x (t), calculates x (t) Upper lower envelope average value m1
(2-3) calculates x (t) and m1Difference h10, judge h10It whether is intrinsic mode function, if it is, enabling eigen mode Measure c1=h10, and it is transferred to (2-5), if it is not, then being transferred to (2-4);
(2-4) calculates k h repeatedly1(k-1)-m1k=h1k, until h1kMeet intrinsic mode function condition, enables c1=h1k.Its In, m1kFor h1(k-1)Upper lower envelope average value;
(2-5) separates c from x (t)1, enable r1=x (t)-c1, by r1Regard new x (t) as, constantly repeats (2-2)-(2- 4) it, is often repeated once, so that it may obtain new intrinsic modulus, i.e. r1-c2=r2,...,rn-1-cn=rn, stop standard until meeting Then;
The intrinsic modulus and residual volume that (2-6) record epicycle is decomposed, it is assumed that isolate q intrinsic modulus, then enable c1N =c1,c2N=c2,...,cqN=cq,rqN=rq, and the number of iterations N is added 1, if N is less than predetermined maximum number of iterations K, turn Enter (2-2), otherwise calculating final intrinsic modulus and residual volume, calculation formula isWith
Step (3) calculates the energy ratio and approximate entropy of each intrinsic modulus and residual volume, constitutes m-th of sensor and corresponds to sample This feature vector.
Step (4) classifies to training sample using the feature vector construction optimal decision function of m sample class, classifies Decision function are as follows:
F (x)=sgn (wTx+b)
Wherein w and b is respectively the weight coefficient vector sum deviation of the categorised decision function, and x is sampling feature vectors, this Two parameters can determine by solving objective function, objective function are as follows:
s.t.yi(wTx+b)-1+ξi≥0
Wherein Q indicates the number of samples in m-th of sensor, ξiFor slack variable, C is penalty coefficient, yiFor sample mark Label, the classifier constructed at this time are denoted as SVMm.
Probability between the class of step (5) calculating sampleIts Middle i, j indicate i-th kind of operating status, jth kind operating status, and f is the value that sample point inputs that decision function obtains, and A and B are to pass through Following objective function is minimized to obtain
Sample point x (t) belongs to the Probability p of i-th kind of operating status under the model for step (6) estimationi=P (y=i | x), Calculation method is to optimize following objective function
Step (7), parameter m value add 1, to next sensor repeat step (2)-step (6), establish it is next support to Amount machine classifier;Until all the sensors analysis finishes, indicate that fault diagnosis model foundation finishes.
Combined failure sample input model is carried out fault diagnosis by step (8).Specifically, by being acquired in combined failure sample SVM1 is inputted from the vibration signal of sensor 1, inputs SVM2 for being acquired in combined failure sample from the vibration signal of sensor 2, Probability Estimation is respectively obtained, finally resulting probability is averaging, obtains final sumbission.
Table 1 is the diagnostic result of 5 combined failure samples.From table it can be seen that, either SVM1 or SVM2, all without Method individually detects combined failure sample while there are two types of failures for tool, and after their probability Estimation is integrated, it finds that therefore The probability that barrier 1 and failure 2 occur is closer to, thus it can be concluded that sample has the conclusion there are two types of failure simultaneously, demonstrate this The validity of invention algorithm.
Several combined failure sample diagnostic results of Table I
It should be understood that above-described embodiment is preferred embodiments of the present invention, but embodiments of the present invention are not It is restricted to the described embodiments.For those of ordinary skills, it can be modified or changed according to the above description, and All these modifications and variations should all belong to the protection domain of appended claims of the present invention.

Claims (2)

1. a kind of gearbox of wind turbine combined failure diagnostic method based on support vector machines probability Estimation, it is characterised in that: Corresponding supporting vector machine model is established for each sensor, the output result of model is certain sample point category under the model The output result of all models is averaging, is somebody's turn to do to each operating status by the corresponding probability under different operating statuses Sample point belongs to the true probability of the operating status, specifically:
Step (1), from the vibration signal extracted in gear-box monitoring data library under different working condition as sample;Assuming that in tooth M sensor is installed on roller box, according to the corresponding sample of each sensor of the sequential processes of sensor, m=1 is initialized, to m A sensor seeks the fault signature of its respective sample in the steps below;
Step (2) corresponds to the sample point x (t) in sample to m-th of sensor and carries out EEMD decomposition, obtains final intrinsic modulus And residual volume;
Step (3) calculates the energy ratio and approximate entropy of each intrinsic modulus and residual volume, constitutes m-th of sensor and corresponds to sample Feature vector;
Step (4) divides training sample using the feature vector construction optimal decision function of m-th of sensor respective sample Class, categorised decision function are as follows:
F (x)=sgn (wTx+b)
Wherein w and b is respectively the weight coefficient vector sum deviation of the categorised decision function, and x is sampling feature vectors, the two Parameter can determine by solving objective function, objective function are as follows:
s.t.yi(wTx+b)-1+ξi≥0
Wherein Q indicates the number of samples in m-th of sensor, ξiFor slack variable, C is penalty coefficient, yiFor sample label, this When the classifier that constructs be denoted as SVMm;
Probability between the class of step (5) calculating sampleWherein i, j Indicate i-th kind of operating status, jth kind operating status, f is the value that sample point inputs that decision function obtains, and A and B are to pass through minimum Change following objective function to obtain
Sample point x (t) belongs to the Probability p of i-th kind of operating status under the model for step (6) estimationi=P (y=i | x), calculating side Method is to optimize following objective function
Step (7) repeats step (2)-step (6) to next sensor and establishes next support vector machines even m adds one Classifier;Until all the sensors analysis finishes, indicates that fault diagnosis model is established and complete;
Test sample input model is carried out fault diagnosis by step (8);Specifically, when test sample data are acquired from multiple and different Sensor then needs that the vibration data separation of different sensors will be belonged in sample, and sequentially inputs the correspondence having had built up The output probability of all models is finally averaging, obtains final result by the supporting vector machine model of sensor.
2. the gearbox of wind turbine combined failure diagnosis side according to claim 1 based on support vector machines probability Estimation Method, it is characterised in that: EEMD decomposable process is as follows in step 2):
(2-1) initializes the number of iterations N=1, and predetermined maximum number of iterations is K;
The white Gaussian noise w (t) of computer random generation is added as new x (t) in (2-2) in x (t), calculates the upper of x (t) Lower envelope average value m1
(2-3) calculates x (t) and m1Difference h10, judge h10It whether is intrinsic mode function, if it is, enabling intrinsic modulus c1 =h10, and it is transferred to (2-5), if it is not, then being transferred to (2-4);
(2-4) calculates k h repeatedly1(k-1)-m1k=h1k, until h1kMeet intrinsic mode function condition, enables c1=h1k, wherein m1k For h1(k-1)Upper lower envelope average value;
(2-5) separates c from x (t)1, enable r1=x (t)-c1, by r1Regard new x (t) as, constantly repeats (2-2)-(2-4), it is every heavy It is multiple primary, so that it may to obtain new intrinsic modulus, i.e. r1-c2=r2,...,rn-1-cn=rn, until meeting stopping criterion;
The intrinsic modulus and residual volume that (2-6) record epicycle is decomposed, it is assumed that isolate q intrinsic modulus, then enable c1N=c1, c2N=c2,...,cqN=cq,rqN=rq, and the number of iterations N is added 1, if N is less than predetermined maximum number of iterations K, it is transferred to (2- 2), otherwise calculating final intrinsic modulus and residual volume, calculation formula isWith
CN201710232999.6A 2017-04-11 2017-04-11 Gearbox of wind turbine combined failure diagnostic method based on support vector machines probability Estimation Active CN107036808B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710232999.6A CN107036808B (en) 2017-04-11 2017-04-11 Gearbox of wind turbine combined failure diagnostic method based on support vector machines probability Estimation

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710232999.6A CN107036808B (en) 2017-04-11 2017-04-11 Gearbox of wind turbine combined failure diagnostic method based on support vector machines probability Estimation

Publications (2)

Publication Number Publication Date
CN107036808A CN107036808A (en) 2017-08-11
CN107036808B true CN107036808B (en) 2019-04-19

Family

ID=59535911

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710232999.6A Active CN107036808B (en) 2017-04-11 2017-04-11 Gearbox of wind turbine combined failure diagnostic method based on support vector machines probability Estimation

Country Status (1)

Country Link
CN (1) CN107036808B (en)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107830996B (en) * 2017-10-10 2020-11-03 南京航空航天大学 Fault diagnosis method for aircraft control surface system
CN108507783B (en) * 2018-03-14 2020-06-23 湖南大学 Rotary machine composite fault diagnosis method based on group decomposition
CN108872781B (en) * 2018-05-08 2020-11-06 广东昊阳电力建设有限公司 Analysis method and device based on intelligent inspection of electric power facilities

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR101249576B1 (en) * 2010-09-13 2013-04-01 한국수력원자력 주식회사 Rotating Machinery Fault Diagnostic Method and System Using Support Vector Machines
CN102944416B (en) * 2012-12-06 2015-04-01 南京匹瑞电气科技有限公司 Multi-sensor signal fusion technology-based fault diagnosis method for wind turbine blades

Also Published As

Publication number Publication date
CN107036808A (en) 2017-08-11

Similar Documents

Publication Publication Date Title
CN107036808B (en) Gearbox of wind turbine combined failure diagnostic method based on support vector machines probability Estimation
CN105973594B (en) A kind of rolling bearing fault Forecasting Methodology based on continuous depth confidence network
Han et al. Intelligent fault diagnosis method for rotating machinery via dictionary learning and sparse representation-based classification
CN102944769A (en) Fault diagnosis method of power transformer based on extreme learning machine
CN104915747A (en) Electricity generation performance evaluation method of generator set and equipment thereof
CN103064008A (en) Nolinear analog circuit soft fault diagnostic method based on Hilbert-huang transform
CN105065212A (en) Checking method and system of wind generation sets of wind power plant
CN104300557B (en) Power system active power oscillations type identification method based on support vector machine
CN108732528A (en) A kind of digitalized electrical energy meter method for diagnosing faults based on depth confidence network
CN109738776A (en) Fan converter open-circuit fault recognition methods based on LSTM
CN109613428A (en) It is a kind of can be as system and its application in motor device fault detection method
CN109033612A (en) A kind of Diagnosis Method of Transformer Faults based on vibration noise and BP neural network
Reddy et al. Application of artificial neural networks for identification of unbalance and looseness in rotor bearing systems
CN103335840A (en) Intelligent diagnosis method for faults of mining drilling machine gearbox
CN105157821A (en) Rotary machinery vibration fault diagnosis and quantitative analysis method
CN110287552A (en) Based on the motor bearings fault diagnosis method and system for improving random forests algorithm
CN110160789A (en) A kind of Wind turbines Method for Bearing Fault Diagnosis based on GA-ENN
Altaf et al. Fault diagnosis in a distributed motor network using Artificial Neural Network
Long et al. A CBA-KELM-based recognition method for fault diagnosis of wind turbines with time-domain analysis and multisensor data fusion
Cao et al. Fault diagnosis of wind turbine gearbox based on deep bi-directional long short-term memory under time-varying non-stationary operating conditions
Cao et al. Intelligent fault diagnosis of wind turbine gearbox based on Long short-term memory networks
CN105116323B (en) A kind of electrical fault detection method based on RBF
CN108548669A (en) A kind of engineer equipment driving unit fault diagnostic method and system
Durbhaka et al. Fault behaviour pattern analysis and recognition
CN105426665B (en) Method is determined based on the DYNAMIC RELIABILITY of status monitoring

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