CN103869162B  Dynamic signal phasor measurement method based on time domain quasisynchronization  Google Patents
Dynamic signal phasor measurement method based on time domain quasisynchronization Download PDFInfo
 Publication number
 CN103869162B CN103869162B CN201410078765.7A CN201410078765A CN103869162B CN 103869162 B CN103869162 B CN 103869162B CN 201410078765 A CN201410078765 A CN 201410078765A CN 103869162 B CN103869162 B CN 103869162B
 Authority
 CN
 China
 Prior art keywords
 quasi
 frequency
 time domain
 signal
 phasor measurement
 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
Links
 238000000691 measurement method Methods 0.000 title claims abstract description 17
 238000005070 sampling Methods 0.000 claims abstract description 53
 238000004422 calculation algorithm Methods 0.000 claims abstract description 19
 238000005259 measurement Methods 0.000 claims abstract description 15
 238000001228 spectrum Methods 0.000 claims abstract description 8
 230000001360 synchronised Effects 0.000 claims description 15
 238000000034 method Methods 0.000 claims description 9
 238000004450 types of analysis Methods 0.000 claims 1
 238000004458 analytical method Methods 0.000 abstract description 6
 230000000875 corresponding Effects 0.000 description 5
 238000005316 response function Methods 0.000 description 3
 238000010586 diagram Methods 0.000 description 2
 230000000694 effects Effects 0.000 description 2
 238000004364 calculation method Methods 0.000 description 1
 238000006243 chemical reaction Methods 0.000 description 1
 150000001875 compounds Chemical class 0.000 description 1
 238000010276 construction Methods 0.000 description 1
 238000009795 derivation Methods 0.000 description 1
 238000005516 engineering process Methods 0.000 description 1
 238000001914 filtration Methods 0.000 description 1
 239000004065 semiconductor Substances 0.000 description 1
 238000004088 simulation Methods 0.000 description 1
 230000003595 spectral Effects 0.000 description 1
 238000010183 spectrum analysis Methods 0.000 description 1
Abstract
The invention discloses a dynamic signal phasor measurement method based on time domain quasisynchronization. The dynamic signal phasor measurement method based on the time domain quasisynchronization comprises the following steps: estimating the fundamental wave frequency of a sampled signal by a time domain quasisynchronization sampling algorithm; carrying out time domain quasisynchronization on the sampled signal by the estimated value of the fundamental wave frequency; reconstructing a quasisynchronization sampling sequence by Newton interpolation; carrying out frequency domain analysis on the reconstructed quasisynchronization sampling sequence by FFT (fast Fourier transform) to obtain a dynamic signal phasor measurement result. According to the method, the influence on the measurement precision of the traditional phasor measurement method by frequency spectrum leakage caused by nonsynchronous sampling is avoided. The computation complexity of the dynamic signal phasor measurement method based on the time domain quasisynchronization is smaller than the computation complexity of the traditional phasor measurement method based on the discrete Fourier transform, and the dynamic signal phasor measurement method based on the time domain quasisynchronization is easy to realize in an embedded system.
Description
Technical field
The present invention relates to signal phasor measurement field, specifically one kind are based on time domain quasi synchronous Dynamic Signal phasor measurement
Method.
Background technology
With widely using of the nonlinearloads such as power electronic equipment, semiconductor device, power quality problem layer goes out not
Thoroughly.How measuring and to analyze the practical situation of power system in real time, thus improving the quality of power supply, having become as power train in recent years
Emphasis in system research field and focus.
At present, in power system, voltage and current is all the variations per hour doing sinusoidal variations in time, all can use phasor representation.
But, IEEE Standard 13441995 define phasor premise be under steady state conditions, a reactor, that is, the amplitude of signal, frequency and
Phase angle all keeps constant, and under rated frequency, the instantaneous measure of phase angle keeps constant with absolute reference time relatively.But
In practical application, when signal frequency deviates rated frequency, phase angle changes with frequency, thus introducing error it is difficult to obtain accurately
Value of calculation.
Existing phasor measurement algorithm mainly has cross zero detecting method, discrete Fourier transform (DFT) method etc..Cross zero detecting method
Be compare intuitively a kind of synchronous phasor measuring method it is only necessary to by the zero crossing moment of tested power frequency component with sometime mark
Standard compares and can draw phase angle difference；Cross zero detecting method principle is simple, it is easy to accomplish, but its precision is not high, and be easily subject to harmonic wave, make an uproar
The impact of sound component.When electrical network medium frequency no offsets, DFT algorithm can accurately measure amplitude and the phase place of signal, and its
Computational accuracy is not affected by Constant Direct Current component and integer harmonic component；But when mains frequency offsets, due to asynchronous
The spectrum leakage that sampling causes, the precision of phasor measurement can decline rapidly.
For this reason, a kind of Dynamic Signal phasor measurement method studying high accuracy is ensureing power system safety and stability operation
Aspect is significant.
Content of the invention
The invention provides a kind of be based on time domain quasi synchronous Dynamic Signal phasor measurement method, it is to avoid nonsynchronous sampling
The impact to existing phasor measurement method certainty of measurement for the spectrum leakage causing；Measured based on time domain quasi synchronous Dynamic Signal phase
Quantity algorithm computational complexity is less than the existing phasor measurement algorithm based on discrete Fourier transform, and quasi synchronous based on time domain
Dynamic Signal phasor measurement algorithm is easy to realize in embedded systems.
For solving abovementioned technical problem, solution proposed by the present invention is：Estimated using time domain quasisynchro sampling algorithm
The fundamental frequency of sampled signal, does time domain plesiochronousization using fundamental frequency estimated value to sampled signal, by Newton interpolating method
Reconstruct quasisynchro sampling sequence, carries out frequencydomain analysiss using FFT to the quasisynchro sampling sequence of interpolation reconstruction, obtains Dynamic Signal
Phasor measurement result.
Technical scheme is as follows：
One kind is based on time domain quasi synchronous Dynamic Signal phasor measurement method, comprises the following steps：
Step one：Electric power signal is sampled, individually sampled signal is saved as crude sampling sample after sampling, with
When to sampled signal add wave filter, to filter harmonic wave and noise jamming；
Step 2：Adopt quasisynchro sampling algorithm to estimate fundamental frequency the sampled signal after filtered process, obtain base
Wave frequency estimated value f_{g}；
Step 3：Using the frequency estimation f obtaining in step 2_{g}And sampling number N calculates in the signal period
Obtain quasisynchro sampling cycle λ, with λ as steplength, using Newton interpolating method to the discrete sequence in sample original in step one
Row do time domain plesiochronousization and process, and interpolation reconstruction obtains quasisynchro sampling sequence；
Step 4：To the quasisynchro sampling sequence obtaining in step 3, a signal period is intercepted using rectangular window, carries out
FFT spectrum is analyzed, and obtains the frequency domain information of signal, and calculates frequency, amplitude and the phase parameter of electric power signal, obtains dynamic
Signal phasor measurement result.
Described method, the selection rule of step one median filter is：
Choose triangular selfconvolution window digital band pass FIR filter, lower stopband edge frequency is 40Hz, lower passband edge frequency
For 46Hz, upper passband edge frequency is 54Hz, and upper stopband edge frequency is 60Hz, passband ripple 0.01, stopband ripple 0.1.Its
In, the method for designing of triangular selfconvolution window digital band pass FIR filter is：First, according to the finger to stopband attenuation and intermediate zone
Mark requires, and selects triangular selfconvolution window, and estimates length of window；Secondly, the frequency response function of construction ideal digital wave filter,
And ideal unitary impulse response is obtained according to ideal frequency response function；Finally, to impulse response function plus triangular selfconvolution window
Obtain design result.
Described method, in step 2, quasisynchro sampling algorithm parameter selection rule is：
Single iteration points D is 64, and iterationses P is 5.
Mode below by way of theoretical derivation illustrates to the technique effect that the present invention reaches.
In time domain, continuous harmonic signal may generally be expressed as following form
In formula, k is overtone order, represents fundamental wave during k=1；A_{k}For kth subharmonic amplitude；T is the time；F is fundamental signal
Frequency；θ_{k}Initial phase angle for kth subharmonic.
Ignore the quantization error in analogdigital conversion process, and the various random error in measurement process, using sampling frequency
Rate is f_{s}Analysis System for Power Quality obtain N_{c}Individual sample
In formula, n is more than or equal to 0 and to be less than or equal to N_{c} 1 integer.
In order to carry out fundamental frequency estimation, filtered using the sample that triangular selfconvolution window band filter obtains to sampling
Ripple is processed, and filters higher hamonic wave, obtains signal as follows
U (n)=A_{1}sin(2πnf/f_{s}+θ_{1}) in (3) formula, A_{1}For fundamental voltage amplitude；θ_{1}For fundamental wave initial phase angle.
Quasisynchro sampling algorithm recurrence formula is as follows
In formula：The subscript " 1 " of " X " represents the 1st quadrature computing；D be defined synchronous sampling algorithm single iteration points, D=
64；ρ_{i}For the corresponding weight coefficient of numerical quadrature formula, for complexification compound trapezoidal integeration, ρ_{0}=ρ_{D}=0.5, ρ_{1}=ρ_{2}=...=
ρ_{D1}=1.
Estimate to obtain phase difference θ in the Δ t time, then signal fundamental frequency using quasisynchro sampling algorithm recurrence formula
Estimated value f_{g}For
In formula, ω represents fundamental wave angular frequency, and Δ θ represents fundamental wave phase angle difference, and Δ t represents the corresponding time difference of Δ θ.
According to the quasisynchronous algorithm parameter of patent requirements of the present invention, estimated frequency error is less than 3 × 10^{10}Hz.Using (6)
Frequency estimation, the original sampled signal sample in step one is carried out with time domain plesiochronousization, interpolation reconstruction obtains plesiochronous
Sample sequence.Interpolation polynomial is
P (x)=u [x_{0}]+u[x_{0},x_{1}](xx_{0})+L+u[x_{0},x_{1},L,x_{m}](xx_{0})L(xx_{m1}) (7)
In formula, x is sampled signal time interval, and subscript m represents mth sampling time interval, and u [...] is difference coefficient, u [x_{0}]=u
(x_{0}), k jump business is
With rectangular window, is intercepted to interpolation reconstruction signal one signal period, carry out frequencydomain analysiss using FFT, in signal spectrum
Middle searching harmonic wave corresponding peak value spectral line, by analyzing frequency, phase place and the magnitude parameters of peak value and then acquisition electric power signal,
Complete the phasor measurement of Dynamic Signal.
In sum, of the present invention a kind of it is based on time domain quasi synchronous Dynamic Signal phasor measurement method, it is to avoid non
Under synchronized sampling, spectrum leakage and fence effect that FFT causes, and the method computational complexity is less than traditional windowed interpolation
Fft algorithm is it is easy to realize in embedded systems.
The invention will be further described below in conjunction with the accompanying drawings.
Brief description
Fig. 1 is the program flow diagram realized in the present invention based on time domain quasi synchronous Dynamic Signal phasor measurement；
Fig. 2 is quasisynchro sampling algorithm iteration schematic diagram in the present invention；
Fig. 3 is Newton interpolating method reconfiguration principle figure in the present invention.
Specific embodiment
The program circuit that the present invention realizes based on the quasi synchronous harmonic analysis and measurement method of time domain is as shown in Figure 1.
As shown in figure 1, the first step, using Analysis System for Power Quality, the signal of input is sampled, obtain N_{c}Individual sample
This, N_{c}For natural number, f_{s}For sample frequency.
Second step, with triangular selfconvolution window band filter to the N collecting_{c}Individual sample is filtered denoising.Three
The parameter of angle selfconvolution window band filter is：Triangular selfconvolution window band filter, lower stopband edge frequency is 40Hz, lower logical
Belt edge frequency is 46Hz, and upper passband edge frequency is 54Hz, and upper stopband edge frequency is 60Hz, passband ripple 0.01, stopband
Ripple 0.1.
3rd step, is iterated to the signal after filtering and noise reduction using quasisynchro sampling algorithm, estimates the base of measured signal
Wave frequency, obtains fundamental frequency f_{g}, wherein quasisynchro sampling single iteration points D is 64, and iterationses P is 5.Frequency Estimation
Comprise the following steps that：
In t_{1}The P*D+1 discrete sample values in moment are iterated, as shown in Fig. 2 the fundamental wave respectively obtaining P iteration is real
Portion X^{J} _{a}With fundamental wave imaginary part X^{J} _{b}, then t_{1}Fundamental phase θ in moment_{1}For
T can be obtained in the same manner_{2}Fundamental phase θ in moment_{2}, then fundamental frequency estimated value f_{g}For
In formula, ω represents fundamental wave angular frequency, and Δ θ represents fundamental wave phase angle difference, and Δ t represents the corresponding time difference of Δ θ.
4th step, using fundamental frequency estimated value f_{g}, using principle shown in Fig. 3, time domain is done to crude sampling sample plesiochronous
Change is processed, the signal sequence u under the conditions of interpolation reconstruction synchronized sampling_{i}(k).
5th step, with sequence u to interpolation reconstruction for the rectangular window_{i}K () is blocked, intercept the discrete sequence of a signal period
Row u_{i}(k_{i})(k_{i}=0,1 ..., N1).Carry out FFT spectrum analysis to intercepting sequence, obtain.
In formula, k_{0}=f_{g}N/f_{s}；W_{R}K () is rectangular window frequency spectrum.
6th step, to Y (k_{i}) be analyzed, in Y (k_{i}) the middle searching corresponding peak value of each harmonic, divided by analyzing peak value
Huo Qu not the frequency of each harmonic, amplitude and phase place.
So far, complete the measurement of Dynamic Signal phasor, the method for the present invention is applicable to electric power signal fundamental wave frequency deviation
The situation of no frequency deviation.
Method provided by the present invention is below applied to carry out emulation experiment, to verify the reliability of method provided by the present invention
Property.
Produce an emulation signal using MATLAB, signal model is as follows
Y (t)=220sin (2 π f_{0}t+θ_{1})+2.3936sin(6πf_{0}t+θ_{3})+1.3442sin(10πf_{0}t+θ_{5}) (12) formula
In, fundamental frequency f_{0}Value in 49.5～50.5Hz, steplength 0.1Hz；Fundamental voltage amplitude A_{1}For 220V, phase theta_{1}=π/3rad；3
Subharmonic and 5 order harmonic components amplitudes are respectively A_{3}=2.3936V, A_{5}=1.3442V, phase place is respectively θ_{3}=π/4rad, θ_{5}=
π/6rad.
Using sample frequency f_{s}The Analysis System for Power Quality collection N of=3200Hz_{c}=512 samples are analyzed.Emulation
Parameter value is as follows
Quasisynchro sampling algorithm：Single iteration points D is 64, and iterationses P is 5
Newton interpolating method：Interpolation is counted as 11 points
FFT frequencydomain analysiss：Using rectangular window, the fft analysis of N=64 points
As shown in table 1, in table 1, aEb represents a × 10 to simulation result^{b}.
Referring to table 1, when frequency fluctuation is ± 0.5Hz, electric power signal fundamental frequency estimates that absolute error is respectively less than
8.70E10, fundamental voltage amplitude absolute error is respectively less than 4.60E9, and fundamental phase absolute error is respectively less than 5.90E11；3 subharmonic
Amplitude Estimation error is respectively less than 1.10E08, and 3 subharmonic phase angle absolute errors are respectively less than 3.20E09；5 subharmonic Amplitude Estimation
Error is respectively less than 9.80E09, and 5 subharmonic phase angular estimation absolute errors are respectively less than 2.50E08.Therefore, demonstrate institute of the present invention
The feasibility of Dynamic Signal phasor measurement method providing and correctness.
Table 1 carries out the result of emulation experiment using method provided by the present invention
Claims (2)
1. a kind of time domain quasi synchronous Dynamic Signal phasor measurement method that is based on is it is characterised in that comprise the following steps：
Step one：Electric power signal is sampled, individually sampled signal is saved as crude sampling sample after sampling, simultaneously right
Sampled signal adds wave filter, to filter harmonic wave and noise jamming, then the sampled signal after filtered process is carried out step 2
Process；
Step 2：Adopt quasisynchro sampling algorithm to estimate fundamental frequency the sampled signal after filtered process, obtain fundamental wave frequency
Rate estimated value f_{g}；
Step 3：Using fundamental frequency estimated value f obtaining in step 2_{g}And sampling number N calculates in the signal period
To quasisynchro sampling cycle λ, with λ as steplength, using Newton interpolating method to the discrete serieses in sample original in step one
Do time domain plesiochronousization to process, interpolation reconstruction obtains quasisynchro sampling sequence；
Step 4：To the quasisynchro sampling sequence obtaining in step 3, a signal period is intercepted using rectangular window, carries out FFT
Spectrum analyses, obtain the frequency domain information of signal, and calculate frequency, amplitude and the phase parameter of electric power signal, are dynamically believed
Number phasor measurement result；
The selection rule of step one median filter is：
Choose triangular selfconvolution window band filter, lower stopband edge frequency is 40Hz, lower passband edge frequency is 46Hz, upper logical
Belt edge frequency is 54Hz, and upper stopband edge frequency is 60Hz, passband ripple 0.01, stopband ripple 0.1.
2. method according to claim 1 it is characterised in that in step 2 quasisynchronous algorithm parameter selection rule be：
Single iteration points D is 64, and iterationses P is 5.
Priority Applications (1)
Application Number  Priority Date  Filing Date  Title 

CN201410078765.7A CN103869162B (en)  20140305  20140305  Dynamic signal phasor measurement method based on time domain quasisynchronization 
Applications Claiming Priority (1)
Application Number  Priority Date  Filing Date  Title 

CN201410078765.7A CN103869162B (en)  20140305  20140305  Dynamic signal phasor measurement method based on time domain quasisynchronization 
Publications (2)
Publication Number  Publication Date 

CN103869162A CN103869162A (en)  20140618 
CN103869162B true CN103869162B (en)  20170208 
Family
ID=50907913
Family Applications (1)
Application Number  Title  Priority Date  Filing Date 

CN201410078765.7A Active CN103869162B (en)  20140305  20140305  Dynamic signal phasor measurement method based on time domain quasisynchronization 
Country Status (1)
Country  Link 

CN (1)  CN103869162B (en) 
Families Citing this family (16)
Publication number  Priority date  Publication date  Assignee  Title 

CN104181389B (en) *  20140702  20170215  中国农业大学  Phasor measurement method in electricpower system 
CN104122443B (en) *  20140804  20170215  国家电网公司  Adjacent harmonic and interharmonic separation and measurement method under IEC (international electrotechnical commission) framework 
CN104198872B (en) *  20140929  20171024  徐雪松  Online equipment for monitoring power quality and method 
CN104407197B (en) *  20141127  20170627  湖南大学  A kind of method of the signal phasor measurement based on trigonometric function iteration 
CN104459419A (en) *  20141229  20150325  国网四川省电力公司电力科学研究院  Measuring error compensation method and device for nonsynchronous periodic signal adoption 
CN104898041A (en) *  20150630  20150909  绵阳市维博电子有限责任公司  Phasesensitive rail signal detecting method and system 
CN105388361B (en) *  20151231  20180123  武汉大学  Twoway interpolation synchronizes the FFT electric harmonic detection methods of sample sequence 
CN106405229B (en) *  20160830  20190430  威胜集团有限公司  A kind of fundamental wave and harmonic wave electric energy gauging method 
CN107085144B (en) *  20170428  20190820  珠海泰芯半导体有限公司  A kind of method of rapid survey Harmonious Waves in Power Systems 
CN108362939B (en) *  20180131  20200623  成都泰格微波技术股份有限公司  Frequency domain parameter measuring method of linear frequency modulation signal 
JP6649419B2 (en) *  20180326  20200219  ファナック株式会社  Encoder signal processing device and encoder 
CN109507495B (en) *  20181017  20201215  华北水利水电大学  Variablewindowlength quasisynchronization gridconnected parameter measurement method 
CN109632070B (en) *  20190110  20201124  东北大学  Digital array time domain quasisynchronous calibration method based on Newton interpolation 
CN110068729B (en) *  20190422  20211105  南京磐能电力科技股份有限公司  Signal phasor calculation method 
CN110687377B (en) *  20191012  20210730  广东电网有限责任公司  Online monitoring data processing method and device for distributed energy system 
CN112067887B (en) *  20200909  20210827  山东大学  Method for calculating phase quantity under condition of sampling value loss based on filter orthogonal characteristic 
Citations (5)
Publication number  Priority date  Publication date  Assignee  Title 

US5400366A (en) *  19920709  19950321  Fujitsu Limited  Quasisynchronous detection and demodulation circuit and frequency discriminator used for the same 
EP1367794A2 (en) *  19980130  20031203  Matsushita Electric Industrial Co., Ltd.  Modulation method and radio communication system 
CN101915874A (en) *  20100720  20101215  北海市深蓝科技发展有限责任公司  Harmonic wave detection method based on Fourier transformation 
CN102788901A (en) *  20120814  20121121  上海电器科学研究院  High accuracy synchronous dynamic phasor measurement method 
CN103543335A (en) *  20131030  20140129  国家电网公司  Method for measuring synchronous phasor 

2014
 20140305 CN CN201410078765.7A patent/CN103869162B/en active Active
Patent Citations (5)
Publication number  Priority date  Publication date  Assignee  Title 

US5400366A (en) *  19920709  19950321  Fujitsu Limited  Quasisynchronous detection and demodulation circuit and frequency discriminator used for the same 
EP1367794A2 (en) *  19980130  20031203  Matsushita Electric Industrial Co., Ltd.  Modulation method and radio communication system 
CN101915874A (en) *  20100720  20101215  北海市深蓝科技发展有限责任公司  Harmonic wave detection method based on Fourier transformation 
CN102788901A (en) *  20120814  20121121  上海电器科学研究院  High accuracy synchronous dynamic phasor measurement method 
CN103543335A (en) *  20131030  20140129  国家电网公司  Method for measuring synchronous phasor 
NonPatent Citations (2)
Title 

基于准同步采样的电力系统谐波与间谐波在线检测方法研究;周峰;《中国博士学位论文全文数据库（电子期刊）》;20121015(第10期);论文第3章第32页最后第1段第33页第8段 * 
广域测量系统中相量算法的研究;姜桂秀;《中国优秀硕士学位论文全文数据库（电子期刊）》;20060715(第7期);摘要和论文第5章第25页第1段第30页第3段 * 
Also Published As
Publication number  Publication date 

CN103869162A (en)  20140618 
Similar Documents
Publication  Publication Date  Title 

CN103869162B (en)  Dynamic signal phasor measurement method based on time domain quasisynchronization  
CN103245832B (en)  Based on harmonic wave timefrequency characteristic method for parameter estimation and the analyser of quick Stransformation  
Yang et al.  A precise calculation of power system frequency  
CN102435844B (en)  Sinusoidal signal phasor calculating method being independent of frequency  
CN103454497B (en)  Based on the method for measuring phase difference improving windowed DFT  
CN107247182B (en)  Interharmonic component reduction method based on measured phasor data  
CN104049144A (en)  Synchronous phasor measurement implementing method with filteredout attenuation direct current components  
CN104635094A (en)  Method for improving PMU (power management unit) synchronous phasor measurement precision  
CN103257271A (en)  Device and method for detecting micro grid harmonic wave and interharmonics based on STM32F107VCT6  
CN103116064A (en)  Method and device for detecting voltage fluctuation and flicker based on energy operator and spectrum correction  
CN203287435U (en)  A micro electrical network harmonic wave and interharmonic wave test apparatus based on an STM32F107VCT6  
CN101915874A (en)  Harmonic wave detection method based on Fourier transformation  
CN103904693B (en)  Based on the synchronized method that frequency self adaptation Virtual shipyard is estimated  
CN104502698A (en)  Method and system for measuring frequency of electric power signal  
CN103543331B (en)  A kind of method calculating electric signal harmonic wave and mAcetyl chlorophosphonazo  
CN105487034A (en)  0.05level electronic transformer verification method and system  
CN103983849B (en)  A kind of Electric Power Harmonic Analysis method of realtime highprecision  
CN109507480A (en)  A kind of harmonic detection method and device of neighbouring fundamental wave/harmonic wave  
CN102313857A (en)  Method and device for analyzing fault recording data of power system  
CN105334381A (en)  Method and device for measuring AC active power  
KR100839436B1 (en)  The method of power frequency estimation using the difference between the gain and cosine and sine filter  
Kaiser et al.  Estimation of power systems amplitudes, frequencies, and phase characteristics using energy operators  
Petrovic  Frequency and parameter estimation of multisinusoidal signal  
CN105467209B (en)  A kind of new metal oxide arrester leakage current analysis method  
CN202119835U (en)  Unstable harmonic and interharmonic measuring instrument 
Legal Events
Date  Code  Title  Description 

C06  Publication  
PB01  Publication  
C10  Entry into substantive examination  
SE01  Entry into force of request for substantive examination  
C14  Grant of patent or utility model  
GR01  Patent grant 