CN103149470B - Method of distinguishing transformer magnetizing rush current by transformer winding vibration - Google Patents

Abstract

The invention relates to a method distinguish magnetizing rush current of a power transformer and discloses the method and a calculating method of distinguishing the magnetizing rush current only by depending on a real-time feature of transformer winding vibration and the method is independent of transformer current. The method takes advantages of saturation of magnetic circuit of the transformer to cause change features of winding electrodynamic force. A winding vibration signal is collected and feature modal of the vibration signal is drawn. The magnetizing rush current is distinguished by relative relation of the vibration feature modal. The feature modal chooses 50 HZ power frequency, double frequency and treble frequency signal of the vibration signal. The power frequency modal is corresponding to direct magnetic bias of a magnetic circuit and the treble frequency modal is corresponding to second harmonic current of the magnetizing rush current. A data sample of a power frequency period of vibrating in short time and data formed by a plurality of periodic extension is adopted. The feature modal is drawn by a calculation method of a narrow-band digital filter bank sets so as to distinguish the magnetizing rush current quickly in a power frequency period. The method has the advantages of being good in antijamming capability and good in robustness.

Description

Utilize the method for transformer winding vibration identification transformer excitation flow

Technical field

The present invention relates to the excitation surge current of power transformer and the identification of saturation effects, belong to Power Transformer Condition monitoring and relay protection field, can be used for forming transformer differential protection new principle and criterion.The present invention is applicable to the excitation surge current identification of transformer of generating plant step-up transformer, electrical power trans mission/distribution system, both can be directly used in transformer main protection, and also can be used for running state of transformer on-line monitoring.

Background technology

Power transformer is the key equipment of electric system, realizes the energy transfer and conversion of different electric pressure.Realize voltage at Power Plant Side raise and be connected to power transmission network, realize energy in the transmission of different electric pressure and electricity consumption in grid side.The fault of power transformer not only can burn transformer, more can directly cause large-area power-cuts or electric power to be sent limited.Transformer adopts current differential protection as main protection, and the core technology of transformer differential protection be sensitive, reliably distinguish excitation surge current and short-circuit current.The method of existing identification excitation surge current is mainly identified by the difference current waveform of transformer primary side, secondary current, concrete principle and method mainly comprise a few class: the ratio 1) calculating second harmonic component in difference current, when ratio reaches certain numerical value, think generation excitation surge current; 2) size of the interval angle in recognition differential streaming current, thinks that current waveform exists the long period near zero-crossing point, when the interval angle of formation reaches certain numerical value, is judged to excitation surge current; 3) symmetry of the upper and lower half-wave of recognition differential streaming current waveform, think remove point periodic component difference current waveform upper and lower half-wave existence larger asymmetric time, be then judged to excitation surge current; 4) introduce the voltage calculating excitation reactance of former limit, secondary, think when comparatively macromutation occurs excitation reactance, be judged to excitation surge current.From difference current waveform character, identify that the principle of excitation surge current affects by the Transfer characteristic of primary current mutual inductor and secondary current mutual inductor completely.Existing principle fails very effectively to realize the identification to power transformer excitation surge current, still there is tranformer protection tripping and malfunction that identification error causes in engineering.

The basic reason that power transformer produces excitation surge current is its saturation effects.The manufacturing process of power transformer is very complicated, and large-scale power transformer has that magnetic flux density is increasing, transformer is more and more close to the technological trend that saturated critical point runs.The voltage cataclysm that the fault of electric system and operation cause all may cause the saturation effects of transformer, and then produces excitation surge current.The magnetic flux of static exciter branch road is made up of stable state magnetic flux and transient state magnetic flux, and stable state magnetic flux is relatively stable, and transient state magnetic flux is the key causing saturation effects.Saturation effects after total magnetic flux exceedes the saturation flux of transformer core, starts to occur excitation surge current.

While saturation effects produces excitation surge current, also cause Transformer Winding by the vibration under instantaneous electric power excitation.In the technology of existing identification excitation surge current, basket vibration feature has been left in the basket.The present invention proposes the basket vibration feature identification excitation surge current utilizing saturation effects to cause.The electric power of Transformer Winding is directly proportional to the electric current flow through in winding and magnetic flux density, in excitation surge current process, the maximum rapid increasing of electric current is to 6-8 times of load current, magnetic flux density significantly improves simultaneously, these two factors all cause the electrodynamic increase of winding, and Direct driver transformer winding vibration intensity increases suddenly.On the other hand, the winding of transformer and the manufacturing process of iron core make winding physical construction ideal, and compared with generator unit stator winding technique, usual transformer is not easy to produce winding and loosens, and this makes transformer winding vibration have better consistance.By detecting basket vibration identification saturation effects and excitation surge current, having physical concept and robustness more clearly relative to other by the recognition methods of excitation surge current Current Waveform Characteristics, is a kind of brand-new monitoring method.The core technology main points of the method are: 1) how to identify vibration performance rapidly; 2) noise how in filtering vibration signal, accomplishes reliable recognition.Excitation surge current is the same with short-circuit current, is the electromagnet phenomenon that Transformer occurs suddenly, needs to detect fast and accurately identifies that excitation surge current just can be applied to the locking of shoving of transformer differential protection.This point is different from the technology and calculating method of rotating machinery vibrating failure diagnosis, therefore needs the sudden change by identifying fast in vibration survey, utilizes limited vibrating data collection to extract characteristic component and identifies.

Summary of the invention

The present invention relates to the excitation surge current identification of power transformer, disclose a kind of independent of transformer current, the method for the real-time characteristic identification transformer excitation flow that only relies on transformer winding vibration and algorithm thereof.The method utilizes the change of the saturated winding electric dynamic characteristic caused of transformer magnetic circuit, from multiple station acquisition basket vibration signals of Transformer Winding, extracts the multi-modal feature of vibration signal, utilizes the relative scale of vibration performance mode to differentiate excitation surge current.Characteristic modes selects 50Hz mode, 100Hz mode and 150Hz mode in vibration signal, wherein the DC magnetic biasing in the corresponding magnetic circuit of 50Hz mode, the second harmonic current in the corresponding excitation surge current of 150Hz mode.In addition the present invention adopts the vibration short time-window data in a power frequency period to input as calculating, utilize the cooperation of the multicycle continuation of short time-window data and arrowband infinite impulse response filtering algorithm to realize arrowband modal characteristics to extract, thus identify excitation surge current fast within a power frequency period time, and there is desirable antijamming capability and robustness.The present invention directly reflects the saturation effects of power transformer and the bias direct current component of magnetic flux, has the sensitivity of desirable identification excitation surge current, real-time and reliability.

The vibration of Transformer mainly comprises transformer primary side, vice-side winding vibration, core vibration, cooling fan vibration.Wherein, basket vibration is directly proportional to the size of current in winding and magnetic flux density; Core vibration is relevant to transformer working voltage, and the vibration of no-load running is based on core vibration, and core vibration is little in Transformer change; Fan for cooling transformer and the vibration at work of other annexes less, and the characteristic frequency of vibration signal is in below 50Hz, is easy to distinguish.For excitation surge current, flow through larger current in winding, magnetic circuit occurs saturated, this main body that transformer under this operating mode is vibrated is basket vibration.Basket vibration through the support section of transformer oil and iron core from internal delivery to transformer shell, less in the vibration survey difference of the diverse location of device shell.

The driving source of transformer winding vibration is the electric power of winding suffered by magnetic field, and electric power is directly proportional to the momentary current in winding and magnetic induction density, as formula (1):

F＝nILB (1)

Wherein F is electric power, and n is the number of turn, and I is electric current, L be wire in the length perpendicular to magnetic direction, B is magnetic induction density.

Magnetic induction density when there is saturation effects is as formula (2):

$B\left(t\right)=B_0{e}^{-\frac{t}{{\mathrm{\τ}}_B}}+B_M\cos ({\mathrm{\ω}}_{0}t+c)---\left(2\right)$

Namely by transformer remanent magnetism B ₀the DC component exponentially decayed, and a peak value is B _minterchange steady-state component composition.τ _bbe the time constant of remanent magnetism decay, c is initial phase, ω ₀for the angular velocity that 50Hz is corresponding.

Electric current when there is saturation effects and excitation surge current are as formula (3):

$I\left(t\right)=I_0{e}^{-\frac{t}{{\mathrm{\τ}}_{I0}}}+I_M\cos ({\mathrm{\ω}}_{0}t+d)+I_{2M}{e}^{-\frac{t}{{\mathrm{\τ}}_{I2}}}\cos ({2\mathrm{\ω}}_{0}t+e)---\left(3\right)$

Namely be made up of the second harmonic current component of decaying DC component, power current component, decay.τ _i0the time constant of DC current decay, τ _i2be the time constant of second harmonic current decay, d is power frequency component initial phase, and e is second harmonic component initial phase, ω ₀for the angular velocity of 50Hz.

Formula (2), (3) are brought in (1), definition k ₀(t), k ₂t () represents the ratio of the DC component fundametal compoment relative to second harmonic component in electric current, k _bremanent magnetism attenuation components in (t) expression magnetic induction density and the ratio of AC compounent, that is:

$k_B\left(t\right)=\frac{B_0{e}^{-\frac{t}{{\mathrm{\τ}}_B}}}{B_M};$ $k_0\left(t\right)=\frac{I_0{e}^{-\frac{t}{{\mathrm{\τ}}_{I_0}}}}{I_M};$ $k_2\left(t\right)=\frac{I_{2M}{e}^{-\frac{t}{{\mathrm{\τ}}_{I_2}}}}{I_M}$

Then:

F∝B _M(k _B(t)+cos(ω ₀t+c))×I _M(k ₀(t)+cos(ω ₀t+d)+k ₂(t)cos(2ω ₀t+e)) (4)

F∝B _MI _M((k _B(t)k ₀(t)+0.5cosθ′)+

(k ₀(t)cos(ω ₀t+θ ₀)+k _B(t)cos(ω ₀t+θ ₁)+0.5k ₂(t)cos(ω ₀t+θ ₂))+ (5)

(k _B(t)k ₂(t)cos(2ω ₀t+θ ₃)+0.5cos(2ω ₀t+θ ₄))+

0.5k ₂(t)cos(3ω ₀t+θ ₅))

Namely DC component, fundametal compoment, two harmonics, frequency tripling component is comprised in electric power.Wherein two harmonics are main bodys, and first-harmonic and frequency tripling component are the characteristic signals because saturation effects causes.

The present invention reflects the characteristic modes of saturation effects by extracting in transformer winding vibration, comprise 50Hz mode, 100Hz mode, 150Hz mode, utilize characteristic modes ratio to identify excitation surge current, concrete steps and technical essential as follows:

(1) sensor configuration.Acceleration vibration transducer is installed respectively in the top and bottom of each phase iron core winding of power transformer, sensor signal output line adopts screen layer and ground connection, vibration signal accesses vibration measurement device through cable, and do continuous shaking sampled measurements, sampling rate is not less than 1.8kHz.For 1.8kHz, the every cycle of corresponding power frequency is sampled 36 points.

(2) timing 5ms calculates the amplitude of 2 times of power frequency components of a Three-Phase Transformer basket vibration signal, i.e. 100Hz component amplitude, using this moment up-to-date 36 sampled datas as input, is designated as { S ₀, S ₀..., S ₃₅, utilize complete cycle fourier algorithm, i.e. formula (6), (7), (8), calculate and obtain signal amplitude A, wherein N=2.

$a=\underset{i=0}{\overset{35}{\mathrm{\Σ}}}\left(S_i\sin (\frac{\mathrm{\π}}{18}\·i\·N)\right)---\left(6\right)$

$b=\underset{i=0}{\overset{35}{\mathrm{\Σ}}}\left(S_i\cos (\frac{\mathrm{\π}}{18}\·i\·N)\right)---\left(7\right)$

$A=\sqrt{a^2+b^2}---\left(8\right)$

Form the time series of amplitude A, be designated as { A (k-n), ..., A (k-3), A (k-2), A (k-1), A (k) }, differentiate adjacent 20ms in the 100Hz component amplitude time series whether amount of undergoing mutation, criterion is | A (k)-A (k-4) | and > Δ A, calculating the Sudden Changing Rate of the upper and lower vibration measuring point of A, B, C three-phase windings respectively, when there being a Sudden Changing Rate to exceed Sudden Changing Rate setting threshold value Δ A, then judging that Sudden Changing Rate starts.

(3) after step (2) identifies Sudden Changing Rate startup, the characteristic component in narrow band filter group algorithm accurate Calculation vibration signal is adopted.Choose 20ms and vibrate sampled value, after doing the periodic extension of multiple cycle, as the input that characteristic quantity calculates, by each vibration measuring point of A, B, C three-phase, respectively calculated rate be the mode 1 of 50Hz, the frequency mode 2 that is 100Hz, the frequency mode 3 that is 150Hz, concrete steps are as follows:

A () continuation 20ms sampled signal, to 36 vibration sampled value { S as input ₀, S ₁..., S ₃₅do the continuation in 25 cycles, be designated as { S ₀, S ₁..., S ₃₅, S ₃₆..., S _{36 × 25-1}, amount to 900 input values;

B (), by IIR filter group, extract 50Hz, 100Hz and 150Hz characteristic modes, bank of filters comprises arrowband bandpass filtering and narrow-band band-elimination filtering.50Hz characteristic modes is extracted, adopts the bandreject filtering of 50Hz bandpass filtering cascade 100Hz; 100Hz characteristic modes is extracted, adopts 100Hz bandpass filtering cascade 50Hz bandreject filtering, then cascade 150Hz bandreject filtering; 150Hz characteristic modes is extracted, adopts 150Hz bandpass filtering cascade 100Hz bandreject filtering, then cascade 200Hz bandreject filtering.Adopt above-mentioned cascading filter group to do filtering to 900 of continuation input values to calculate, last 36 that choose output export as filtering, are designated as { X ₈₆₄, X ₈₆₅..., X ₈₉₉.Formula (6), (7), (8) are adopted to calculate last cycle { X ₈₆₄, X ₈₆₅..., X ₈₉₉amplitude, get N and be respectively 1,2,3, obtain the amplitude of the 50Hz characteristic modes of this input signal, 100Hz characteristic modes, 150Hz characteristic modes, be designated as A respectively _m1, A _m2, A _m3.

(4) relative scale of characteristic modes amplitude is calculated according to these 3 relative ratios, utilize fuzzy membership method identification excitation surge current.For Characteristic Ratios its fuzzy membership λ is set ₁for: be 0 when this Characteristic Ratios is less than 0.08 degree of membership, it is 1, time between 0.08 to 0.15 that this Characteristic Ratios is greater than 0.15 degree of membership, from 0 to 1 linear value.For Characteristic Ratios its fuzzy membership λ is set ₃for: be 0 when this Characteristic Ratios is less than 0.15 degree of membership, it is 1, time between 0.15 to 0.25 that this Characteristic Ratios is greater than 0.25 degree of membership, from 0 to 1 linear value.For Characteristic Ratios its fuzzy membership λ is set ₁₃for: be 0 when this Characteristic Ratios is less than 0.2 degree of membership, it is 1, time between 0.2 to 0.4 that this Characteristic Ratios is greater than 0.4 degree of membership, from 0 to 1 linear value.

(5) fuzzy total degree of membership λ=0.4 λ is set ₁+ 0.4 λ ₃+ 0.2 λ ₁₃, when fuzzy total degree of membership λ > 0.75 judges excitation surge current occurs.As fuzzy total degree of membership 0.75 > λ > 0.3, can not determine and whether excitation surge current occurs, now time delay 20ms recalculates the vibration characteristic signals of next power frequency cycle, continues to identify excitation surge current.Not excitation surge current when λ < 0.3 is judged to.

Accompanying drawing explanation

Fig. 1 is the process flow diagram of the vibration performance identification excitation surge current utilizing Transformer Winding.First differentiate whether 100Hz characteristic modes vibration mutation amount occurs, if the amount of undergoing mutation, then start magnetizing inrush current discrimination algorithm; Get a cycle, namely 20ms vibrates sampled value, extracts characteristic modes and amplitude thereof, and calculates corresponding fuzzy membership; Carry out fuzzy diagnosis according to total fuzzy membership λ, be wherein confirmed to be after excitation surge current through fuzzy diagnosis, get back to through time delay involution the amplitude that timing calculates mode of oscillation 2; After fuzzy diagnosis really admits a fault excitation surge current, get back to timing through time delay involution equally and calculate mode of oscillation 2 amplitude.When fuzzy diagnosis is not confirmed whether as excitation surge current, next power frequency cycle vibration signal of sampling, continues to extract vibration performance mode and does fuzzy discrimination.

Fig. 2 is the processing flow chart extracting characteristic modes from the vibration signal of a cycle.The vibration signal illustrating 20ms is through periodic extension, and the bank of filters respectively by bandpass filtering and bandreject filtering cascade obtains the processing procedure of 50Hz, 100Hz and 150Hz vibration performance modal amplitudes.

Fig. 3 is the relative scale of characteristic modes amplitude corresponding fuzzy membership λ ₁, λ ₃and λ ₁₃value function.

Fig. 4 is the principle schematic of the saturated initiation excitation surge current of transformer magnetic circuit.The magnetic flux of static exciter branch road is made up of stable state magnetic flux and remanent magnetism magnetic flux, and stable state magnetic flux is relatively stable, and remanent magnetism magnetic flux is the key causing saturation effects, decays in time.Saturation effects after total magnetic flux exceedes the saturation flux of transformer core, starts to occur excitation surge current, and Transformer Winding produces abnormal vibrations simultaneously.

Embodiment

Install 6 piezoelectric type vibration acceleration transducers altogether at the support zone of transformer shell near the top and bottom of three-phase windings, arrange vibration survey unit on the spot, the electric signal that sensor exports is through shielded cable access vibration survey unit.Vibration survey unit, to 6 road vibration signal samplings, realizes vibration mutation amount and starts differentiation, vibration performance Frequency extraction and excitation surge current fuzzy diagnosis.When identifying excitation surge current; coordinated with protection equipment for transformer, Transformer's Condition Monitoring device or system by the switching value output signal of vibration survey unit, transmit characteristic by communication to protection equipment for transformer, Transformer's Condition Monitoring device or system simultaneously.

Sampling rate selects 1.8KHz, and namely 20ms samples 36 vibratory outputs.The 100Hz characteristic component amplitude that timing 5ms utilizes formula (6), (7), (8) calculate 6 road vibration signals, calculate 6 tunnel amplitude seasonal effect in time series to be separated by the variable quantity of 4 simultaneously, and compare with threshold value Δ A, Dang You mono-tunnel amplitude variable quantity exceedes Δ A and then judges that Sudden Changing Rate starts.

After differentiation Sudden Changing Rate starts, get 36 point sampling values from current vibration sample pointer position, and do 25 periodic extensions, obtain the long data sample of 900.Bank of filters modal calculation is as shown in Figure 2 done respectively to every road long data sample, obtains the characteristic modes amplitude A on each road _m1, A _m2, A _m3.Wherein in digital filter bank, the center passband bandwidth of narrow band filter is 5Hz, and stopband center frequency range is 5Hz.Adopt second order Butterworth bandpass filter to realize narrow-band filtering, the transport function of wave filter is as formula (9).

$\frac{Y}{X}=\frac{b_0{z}^0+b_1{z}^{-1}+b_2{z}^{-2}}{a_0{z}^0+a_1{z}^{-1}+a_2{z}^{-2}}---\left(9\right)$

Each filter coefficient is as follows:

1) bandpass filter of 50Hz characteristic modes is:

a[]＝{1，-1.9594，0.9897}，b[]＝{0.0864，0，-0.0864}

100Hz rejection filter is:

a[]＝{1，-1.873，0.9933}，b[]＝{0.9966，-1.873，0.9966}；

2) bandpass filter of 100Hz characteristic modes is:

a[]＝{1，-1.8699，0.9899}，b[]＝{0.1701，0，-0.1701}

50Hz rejection filter is:

a[]＝{1，-1.9628，0.9931}，b[]＝{0.9966，-1.9628，0.9966}

150Hz rejection filter is:

a[]＝{1，-1.7264，0.9935}，b[]＝{0.9968，-1.7264，0.9968}

3) bandpass filter of 150Hz characteristic modes is:

a[]＝{1，-1.7236，0.9903}，b[]＝{0.2488，0，-0.2488}

100Hz rejection filter is:

a[]＝{1，-1.873，0.9933}，b[]＝{0.9966，-1.873，0.9966}

200Hz rejection filter is:

a[]＝{1，-1.7762，0.9953}，b[]＝{0.2865，0，-0.2865}

Above-mentioned algorithm realization is for the narrow-band filtering of instantaneous 20ms sampled data, narrow band filter group is calculated to last 36 point data exported, utilizes formula (6), 50Hz, 100Hz, 150Hz characteristic component amplitude that (7), (8) calculate each road vibration signal.For the vibration measurement being placed on power transformer device shell, its key is the non-characteristic frequency signal under suppression electromagnetic interference environment, and this algorithm has desirable filter effect.

Utilize the characteristic modes amplitude A on 6 tunnels _m1, A _m2, A _m3relative scale calculate its degree of membership λ respectively ₁, λ ₃, λ ₁₃, calculate total degree of membership λ.Flow process with reference to accompanying drawing 1 differentiates whether excitation surge current occurs, identification structure for the top and bottom vibration transducer of every phase winding adopts "or" logic, as long as there is a measuring point generation excitation surge current, then judge this phase generation excitation surge current, and instantaneous closed vibration survey unit is to the output switch parameter node that excitation surge current should occur mutually.

Claims (2)

1. one kind utilizes the method for power transformer phase-splitting basket vibration characteristic modes identification transformer excitation flow, the method is by extracting the vibration performance mode of Transformer Winding, utilize the relative scale identification saturation effects between different characteristic mode, and then differentiate the generation of excitation surge current, specifically comprise step:

1) sensor configuration: install vibration transducer respectively in the top and bottom of each phase iron core winding of power transformer, vibration signal is positioned at the vibration survey unit of transformer in-situ through cable access, do continuous shaking sampled measurements, sample frequency is not less than 1.8KHz;

2) differentiation that starts of Sudden Changing Rate: select 1.8KHz for sampling rate, timing 5ms calculates 2 times of power frequency component amplitudes of a Three-Phase Transformer vibration signal, i.e. 100Hz component amplitude, using this moment up-to-date 36 sampled datas as input, is designated as { S ₀, S ₁..., S ₃₅; complete cycle fourier algorithm is utilized to calculate signal amplitude and formation time sequence; be designated as { A (k-n); ...; A (k-2); A (k-1); A (k) }; differentiate in the amplitude time series of 100Hz component, whether adjacent 20ms exists Sudden Changing Rate; criterion is | A (k)-A (k-4) | and > Δ A; calculating the Sudden Changing Rate of the upper and lower vibration measuring point of A, B, C three-phase respectively, when there being the Sudden Changing Rate of a measuring point to exceed Sudden Changing Rate setting threshold value Δ A, then judging that Sudden Changing Rate starts;

3) amplitude of characteristic modes is calculated: when step 2) identify after Sudden Changing Rate starts, choose current 20ms and vibrate sampled value, do multiple cycle continuation, calculate the long data sample of input as characteristic quantity; To the mode 3 respectively vibrating measuring point adopt narrow band digital filter group to extract respectively frequency is the mode 1 of 50Hz, frequency is 100Hz mode 2 mutually, frequency is 150Hz, get last cycle data that the filtering of narrow band digital filter group exports, for compute mode amplitude, computing method adopt complete cycle fourier algorithm, and the amplitude result of calculation of mode 1, mode 2, mode 3 is designated as A respectively _m1, A _m2, A _m3;

4) differentiation of excitation surge current: the relative scale calculating characteristic modes amplitude utilize fuzzy membership method identification excitation surge current; Calculate the fuzzy membership λ of 3 ratios respectively ₁, λ ₃and λ ₁₃, setting time between 0.08 to 0.15, λ ₁from 0 to 1 linear value; time between 0.15 to 0.25, λ ₃from 0 to 1 linear value; time between 0.2 to 0.4, λ ₁₃from 0 to 1 linear value; Fuzzy total degree of membership λ=0.4 λ ₁+ 0.4 λ ₃+ 0.2 λ ₁₃, when fuzzy total degree of membership λ > 0.75 judges excitation surge current occurs, through time delay involution, get back to step 2); As fuzzy total degree of membership 0.3 < λ < 0.75, then get back to step 3), reselect up-to-date 20ms and vibrate sampled data, continue to identify excitation surge current; Not excitation surge current when λ < 0.3 is judged to, get back to step 2 through time delay involution).

2. a kind of method utilizing power transformer phase-splitting basket vibration characteristic modes identification transformer excitation flow according to claim 1, the amplitude of wherein said calculating characteristic modes, it is characterized in that the input utilizing short time-window sampled data continuation acquisition long data sample to calculate as characteristic modes, the periodicity of continuation is not less than 25; The narrow band digital filter group of each characteristic modes is made up of narrow band filter and some rejection filter cascades, and wherein the bandpass center frequency range of narrow band filter is 5Hz, and the stopband center frequency range of rejection filter is 5Hz.