CN102944360A - Imbalanced signal extracting method of dual-rotor system with slight speed difference - Google Patents

Abstract

The invention provides an imbalanced signal extracting method of a dual-rotor system with slight speed difference, and aims to solve the problem of difficultly in extracting imbalanced signals of the system in the prior art. With regard to un-damped real-value sinusoid signals, an improved Prony method (IPM for short) can be adopted to rather precisely estimate the frequency of each sinusoid signal. For the characteristics of beat vibration signals in the dual-rotor system with slight speed difference4, a band-pass tracking filtering method is adopted to reduce noise and then samples are extracted from filtered data, the IPM is adopted to realize frequency identification of two adjacent frequency signals; and then a least square method is used for identifying the amplitude and the phase of the beat vibration signals.

Description

A kind of unbalanced signal extracting method of Dual Rotor System with Little Rotation Speed Difference

Technical field

The invention belongs to the mechanical measurement technique field, particularly a kind of unbalanced signal extracting method of Dual Rotor System with Little Rotation Speed Difference.

Background technology

Dual Rotor System with Little Rotation Speed Difference is widely used in the industries such as oil, chemical industry, pharmacy, two synthetic signals that shake of clapping of power frequency component that frequency is close in the vibration signal of this type systematic, and this extracts to unbalanced signal and brings larger difficulty.

Unbalanced signal extracting method for this type of rotor-support-foundation system has two classes at present, one class is not understand shooting method, clap the peak and clap paddy from time domain waveform analysis, then according to the bat principle of shaking, amplitude through simply calculating two rotor oscillations and phase place, the method needs a bat above Measuring Time of cycle of shaking at least, and needs could determine by complicated algorithm which rotor two amplitudes respectively belong to, and the imbalance of one-shot measurement goes heavily rate to only have about 70%.Equations of The Second Kind is to separate shooting method, isolates two power frequency components by signal processing method, thereby extracts separately amplitude and phase place.Traditional solution shooting method generally adopts relevant function method, and the method need to be done digital integration in the cycle is shaken in whole bat could obtain degree of precision, the same shortcoming that needs Measuring Time long that exists, and accurately obtain to clap also difficult of the cycle of shaking..

It is poor that above-mentioned two methods all need to measure accurately the speed of two rotors, less fast difference measurements error will produce larger impact to the signal extraction precision, and internal rotor generally is enclosed in the outer rotor and speed probe can't be installed in the real system, adds the small speed of actual embedded measurement system Measurement accuracy poor also being difficult for and realizes.Although adopt fft analysis can obtain signal frequency from vibration signal, the data sample in needing for a long time just can obtain higher frequency resolution, and the method real-time is too poor and be not suitable for engineering and use.

Summary of the invention

The present invention proposes a kind of unbalanced signal extracting method of Dual Rotor System with Little Rotation Speed Difference, the problem that is difficult to extract with the unbalanced signal that solves in the prior art for this system.

The Prony method can be used than the short data sample and realize the high precision spectrum estimation, be widely used in Fault Signal Analyses in HV Transmission identification, for the real-valued sinusoidal signal of non-decay, adopt improved Prony method (Improved Prony Method is called for short IPM) can estimate more accurately the frequency of each sinusoidal signal.The present invention is directed to bat in the poor birotor of the dead slow speed signal characteristic that shakes, adopt first the logical tracking filter method noise reduction of band, then sample drawn from filtered data, adopt the IPM method to realize the frequency identification of two near by frequency signals, adopt at last least square method to identify amplitude and the phase place of clapping the signal that shakes.

Technical scheme of the present invention is that a kind of unbalanced signal extracting method of Dual Rotor System with Little Rotation Speed Difference may further comprise the steps:

A1 makes that the sampling period is T, and the vibration signal take T as the cycle to described Dual Rotor System with Little Rotation Speed Difference carries out equal interval sampling, obtains one group of data sequence { D (kT) }, k=1, and 2 ... K, K are a positive integer;

A2, getting centre frequency is f ₀, quality factor q＞2 adopt the digital band pass tracking filter that { D (kT) } carried out filtering, obtain { y (kT) }, k=1, and 2 ... K, wherein, the pulsed transfer function of digital band pass tracking filter is

$H\left(z\right)=k\frac{1-z^{-2}}{1+a_1{z}^{-1}+a_2{z}^{-2}},$ Each coefficient in the formula

$k=\frac{\tan \left(\frac{f_{0}T}{2}\right)}{Q+\tan \left(\frac{f_{0}T}{2}\right)+Q{\tan }^2\left(\frac{f_{0}T}{2}\right)}$

$a_1=\frac{2Q{\tan }^2\left(\frac{f_{0}T}{2}\right)-2Q}{Q+\tan \left(\frac{f_{0}T}{2}\right)+Q{\tan }^2\left(\frac{f_{0}T}{2}\right)}$

$a_2=\frac{Q-\tan \left(\frac{f_{0}T}{2}\right)+Q{\tan }^2\left(\frac{f_{0}T}{2}\right)}{Q+\tan \left(\frac{f_{0}T}{2}\right)+Q{\tan }^2\left(\frac{f_{0}T}{2}\right)};$

A3 with a default rule, makes data to { y (kT) } and extracts, the acquisition data sequence x[n] }, n=0 in the formula, 1 ... N-1 is time series, and N is a positive integer;

A4 establishes { x[n] } and contains consisting of of noisy sinusoidal signal

x[n]＝a ₁cos(ω ₁n+φ ₁)+a ₂cos(ω ₂n+φ ₂)+e(n)(1)

N=0 in the formula, 1 ... N-1 is time series, a ₁, a ₂Be unbalanced signal amplitude, ω ₁, ω ₂Be unbalanced signal frequency and φ ₁, φ ₂Be the unbalanced signal phase place,

The capable difference equation matrix of structure N-4 is

X[n]＝[x[n]+x[n+4],x[n+1]+x[n+3],x[n+2]]（2）

The Prony method shows, the homogeneous solution of above-mentioned equation has been specified e ^{I ω}In 4 rank multinomial coefficients, because sampled data all is real number, e ^{I ω}In 4 rank polynomial expressions can be rewritten into formula about cos (2 ω), thereby can solve with the Chebyshev polynomial expression.

By matrix X structure variance matrix G,

$G=\frac{1}{M}\left[\begin{array}{ccc}\mathrm{\Σ}{X}_{n,1}^2& \mathrm{\Σ}{X}_{n,1}{X}_{n,2}& \sqrt{2}\mathrm{\Σ}{X}_{n,1}{X}_{n,3}\\ \mathrm{\Σ}{X}_{n,1}{X}_{n,2}& \mathrm{\Σ}{X}_{n,2}^2& \sqrt{2}\mathrm{\Σ}{X}_{n,2}{X}_{n,3}\\ \sqrt{2}\mathrm{\Σ}{X}_{n,1}{X}_{n,3}& \sqrt{2}\mathrm{\Σ}{X}_{n,2}{X}_{n,3}& 2\mathrm{\Σ}{X}_{n,3}\end{array}\right]---\left(3\right)$

M=N-4 in the formula, all summations are limited to up and down from 0 to M-1, and the minimal eigenvalue of variance matrix G is by its secular equation polynomial expression

c ₁λ ³+c ₂λ ²+c ₃λ+c ₄＝0（4）

Coefficient ask for, this coefficient c is

$c=\left[\begin{array}{c}1\\ -G_{\mathrm{1,1}}-G_{\mathrm{2,2}}-G_{\mathrm{3,3}}\\ G_{\mathrm{1,1}}{G}_{\mathrm{2,2}}+G_{\mathrm{1,1}}{G}_{\mathrm{3,3}}+G_{\mathrm{2,2}}{G}_{\mathrm{3,3}}-G_{\mathrm{1,2}}^2-G_{\mathrm{1,3}}^2-G_{\mathrm{2,3}}^2\\ -G_{\mathrm{1,1}}{G}_{\mathrm{2,2}}{G}_{\mathrm{3,3}}-2G_{\mathrm{1,2}}{G}_{\mathrm{1,3}}{G}_{\mathrm{2,3}}+G_{\mathrm{1,1}}{G}_{\mathrm{2,3}}^2+G_{\mathrm{2,2}}{G}_{\mathrm{1,3}}^2+G_{\mathrm{3,3}}{G}_{\mathrm{1,2}}^2\end{array}\right]---\left(5\right)$

Then minimal eigenvalue is ${\mathrm{\λ}}_0=-c_2/3-(S+T)/2+i\sqrt{3}(S-T)/2---\left(6\right)$

Wherein S, T computing formula are,

$S=\sqrt[3]{R+\sqrt{Q^3+R^2}}$

$T=\sqrt[3]{R-\sqrt{Q^3+R^2}}$

$Q=(3c_3-c_2^2)/9$

$R=(9c_2{c}_3-27c_4-2c_2^3)/54---\left(7\right)$

Minimal eigenvalue λ ₀The characteristic of correspondence vector is

$v=\left[\begin{array}{c}-{\mathrm{\λ}}_0^2+(G_{\mathrm{2,2}}+G_{\mathrm{3,3}}){\mathrm{\λ}}_0-G_{\mathrm{2,2}}{G}_{\mathrm{3,3}}+G_{\mathrm{2,3}}^2\\ G_{\mathrm{1,2}}{G}_{\mathrm{3,3}}-G_{\mathrm{1,3}}{G}_{\mathrm{2,3}}-G_{\mathrm{1,2}}{\mathrm{\λ}}_0\\ \sqrt{2}(G_{\mathrm{1,3}}{G}_{\mathrm{2,2}}-G_{\mathrm{1,2}}{G}_{\mathrm{2,3}}-G_{\mathrm{1,3}}{\mathrm{\λ}}_0)\end{array}\right]---\left(8\right)$

Obtain the minimal eigenvalue λ of G matrix by matrix decomposition ₀With its characteristic of correspondence vector

Then the total least square solution of matrix X is

$v=\frac{1}{\sqrt{2}{v}_{{\mathrm{\λ}}_0}[K+1]}\left[v_{{\mathrm{\λ}}_0}\right[1\left]v_{{\mathrm{\λ}}_0}\right[2]\·\·\·v_{{\mathrm{\λ}}_0}[K\left]\sqrt{2}{v}_{{\mathrm{\λ}}_0}\right[K+1\left]\right]---\left(9\right)$

Known cos (2 θ)=2 (cos (θ)) by the Chebyshev polynomial expression ²-1(10)

According to the Prony Method And Principle, the frequency of any in homogeneous solution v and 2 unknown signalings all satisfies equation

2cos(2ω _k)v ₁+2cos(ω _k)v ₂+v ₃＝0（11）

Then the frequency estimation of unknown signaling is

${\widehat{\mathrm{\ω}}}_1={\cos }^{-1}\left(\frac{v_2+\sqrt{v_2^2-4v_1{v}_3+8v_1^2}}{4v_1}\right)$

${\widehat{\mathrm{\ω}}}_2={\cos }^{-1}\left(\frac{v_2-\sqrt{v_2^2-4v_1{v}_3+8v_1^2}}{4v_1}\right)---\left(12\right)$

Obtain the frequencies omega of unbalanced signal in the Dual Rotor System with Little Rotation Speed Difference vibration signal ₁And ω ₂

A5 utilizes least square method to estimate ω ₁And ω ₂Amplitude and phase place, formula (1) can be write as

x[n]＝a ₁cos(ω ₁n)cos(φ ₁)-a ₁sin(ω ₁n)sin(φ ₁)

+a ₂cos(ω ₂n)cos(φ ₂)-a ₂sin(ω ₂n)sin(φ ₂)+e(n)（13）

If y ₁=a ₁Cos (φ ₁), y ₂=a ₁Sin (φ ₁), y ₃=a ₂Cos (φ ₂), y ₄=a ₂Sin (φ ₂), treat that estimated parameter is y ₁, y ₂, y ₃, y ₄, according to each time observed reading x[n] (n=1,2 ... N), can draw the system of linear equations (14) of Ay=x, ask its least square solution can pick out each parameter

$\left[\begin{array}{cccc}\cos \left({\mathrm{\ω}}_1{t}_1\right)& -\sin \left({\mathrm{\ω}}_1{t}_1\right)& \cos \left({\mathrm{\ω}}_2{t}_1\right)& -\sin \left({\mathrm{\ω}}_2{t}_1\right)\\ \cos \left({\mathrm{\ω}}_1{t}_2\right)& -\sin \left({\mathrm{\ω}}_1{t}_2\right)& \cos \left({\mathrm{\ω}}_2{t}_2\right)& -\sin \left({\mathrm{\ω}}_2{t}_2\right)\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \cos \left({\mathrm{\ω}}_1{t}_N\right)& -\sin \left({\mathrm{\ω}}_1{t}_N\right)& \cos \left({\mathrm{\ω}}_2{t}_N\right)& -\sin \left({\mathrm{\ω}}_2{t}_N\right)\end{array}\right]\left[\begin{array}{c}{y}_1\\ y_2\\ y_3\\ y_4\end{array}\right]=\left[\begin{array}{c}x\left[1\right]\\ x\left[2\right]\\ \·\\ \·\\ \·\\ x\left[N\right]\end{array}\right]---\left(14\right)$

The frequencies omega of unbalanced signal then ₁And ω ₂The estimated value of amplitude and phase place is

$a_1=\sqrt{y_1^2+y_2^2},$ $a_2=\sqrt{y_3^2+y_4^2}$

φ ₁＝arctan(y ₂/y ₁),φ ₂＝arctan(y ₄/y ₃)（15）。

The present invention is directed to bat in the poor birotor of the dead slow speed signal characteristic that shakes, adopt first the logical tracking filter method noise reduction of band, then sample drawn from filtered data, adopt the IPM method to realize the frequency identification of two near by frequency signals, at last before the filtering data as sample (because of wave filter influential to unbalanced signal), adopt amplitude and the phase place of two unbalanced signals of least square method identification.

Description of drawings

Fig. 1 is that the present invention is used for the process flow diagram that faint unbalanced signal extracts.

Embodiment

The present invention is a kind of unbalanced signal extracting method for Dual Rotor System with Little Rotation Speed Difference, comprises following content:

Digital band pass tracking filter link, the pulsed transfer function of wave filter is

$H\left(z\right)=k\frac{1-z^{-2}}{1+a_1{z}^{-1}+a_2{z}^{-2}}$

Each coefficient is as follows in the formula:

$k=\frac{\tan \left(\frac{f_{0}T}{2}\right)}{Q+\tan \left(\frac{f_{0}T}{2}\right)+Q{\tan }^2\left(\frac{f_{0}T}{2}\right)}$

$a_1=\frac{2Q{\tan }^2\left(\frac{f_{0}T}{2}\right)-2Q}{Q+\tan \left(\frac{f_{0}T}{2}\right)+Q{\tan }^2\left(\frac{f_{0}T}{2}\right)}$

$a_2=\frac{Q-\tan \left(\frac{f_{0}T}{2}\right)+Q{\tan }^2\left(\frac{f_{0}T}{2}\right)}{Q+\tan \left(\frac{f_{0}T}{2}\right)+Q{\tan }^2\left(\frac{f_{0}T}{2}\right)}$

Wherein Q is the wave filter quality factor, f ₀Centered by frequency, T is the sampling period.

This wave filter is unattenuated in centre frequency place amplitude, without phase shift, but the Q value is when larger, near the phase shift variations the centre frequency is violent, and transient process is long.For this reason, the Q value should not be selected excessive, is more than 2 times of inner and outer rotors difference on the frequency but generally should make filter passbands wide.

Parameter identification link based on two signals of IPM:

If contain consisting of of a noisy K sinusoidal signal

x[n]＝a ₁cos(ω ₁n+φ ₁)+a ₂cos(ω ₂n+φ ₂)+e(n)(1)

N=0 in the formula, 1 ... N-1 is time series, amplitude a ₁, a ₂, frequencies omega ₁, ω ₂And phase ₁, φ ₂Be unknown parameter.

The capable difference equation matrix of structure N-4 is as follows:

X[n]＝[x[n]+x[n+4]，x[n+1]+x[n+3]，x[n+2]]（2）

Known by the Prony Method And Principle, the homogeneous solution of above-mentioned DIFFERENCE EQUATIONS is e ^{I ω}In the polynomial coefficient in 4 rank, when sampled data is real number, about e ^{I ω}4 rank polynomial expressions can use cos (2 ω) to substitute, thereby can utilize the Chebyshev polynomial solving.

For asking the total least square solution of equation (2), first by matrix X structure modified variance matrix G,

$G=\frac{1}{M}\left[\begin{array}{ccc}\mathrm{\Σ}{X}_{n,1}^2& \mathrm{\Σ}{X}_{n,1}{X}_{n,2}& \sqrt{2}\mathrm{\Σ}{X}_{n,1}{X}_{n,3}\\ \mathrm{\Σ}{X}_{n,1}{X}_{n,2}& \mathrm{\Σ}{X}_{n,2}^2& \sqrt{2}\mathrm{\Σ}{X}_{n,2}{X}_{n,3}\\ \sqrt{2}\mathrm{\Σ}{X}_{n,1}{X}_{n,3}& \sqrt{2}\mathrm{\Σ}{X}_{n,2}{X}_{n,3}& 2\mathrm{\Σ}{X}_{n,3}\end{array}\right]---\left(3\right)$

M=N-4 in the formula, all summations are limited to from 0 to M-1 up and down.Because the G matrix is symmetrical matrix, only need the upper triangle of calculating or lower triangle to get final product.

The minimal eigenvalue of G can pass through its secular equation polynomial expression

c ₁λ ³+c ₂λ ²+c ₃λ+c ₄＝0（4）

Coefficient ask for.Coefficient c is

$c=\left[\begin{array}{c}1\\ -G_{\mathrm{1,1}}-G_{\mathrm{2,2}}-G_{\mathrm{3,3}}\\ G_{\mathrm{1,1}}{G}_{\mathrm{2,2}}+G_{\mathrm{1,1}}{G}_{\mathrm{3,3}}+G_{\mathrm{2,2}}{G}_{\mathrm{3,3}}-G_{\mathrm{1,2}}^2-G_{\mathrm{1,3}}^2-G_{\mathrm{2,3}}^2\\ -G_{\mathrm{1,1}}{G}_{\mathrm{2,2}}{G}_{\mathrm{3,3}}-2G_{\mathrm{1,2}}{G}_{\mathrm{1,3}}{G}_{\mathrm{2,3}}+G_{\mathrm{1,1}}{G}_{\mathrm{2,3}}^2+G_{\mathrm{2,2}}{G}_{\mathrm{1,3}}^2+G_{\mathrm{3,3}}{G}_{\mathrm{1,2}}^2\end{array}\right]---\left(5\right)$

Then minimal eigenvalue is

${\mathrm{\λ}}_0=-c_2/3-(S+T)/2+i\sqrt{3}(S-T)/2---\left(6\right)$

Wherein S, T computing formula are

$S=\sqrt[3]{R+\sqrt{Q^3+R^2}}$

$T=\sqrt[3]{R-\sqrt{Q^3+R^2}}$

$Q=(3c_3-c_2^2)/9$

$R=(9c_2{c}_3-27c_4-2c_2^3)/54---\left(7\right)$

λ ₀The characteristic of correspondence vector is

$v=\left[\begin{array}{c}-{\mathrm{\λ}}_0^2+(G_{\mathrm{2,2}}+G_{\mathrm{3,3}}){\mathrm{\λ}}_0-G_{\mathrm{2,2}}{G}_{\mathrm{3,3}}+G_{\mathrm{2,3}}^2\\ G_{\mathrm{1,2}}{G}_{\mathrm{3,3}}-G_{\mathrm{1,3}}{G}_{\mathrm{2,3}}-G_{\mathrm{1,2}}{\mathrm{\λ}}_0\\ \sqrt{2}(G_{\mathrm{1,3}}{G}_{\mathrm{2,2}}-G_{\mathrm{1,2}}{G}_{\mathrm{2,3}}-G_{\mathrm{1,3}}{\mathrm{\λ}}_0)\end{array}\right]---\left(8\right)$

Obtain the minimal eigenvalue λ of G matrix by matrix decomposition ₀With its characteristic of correspondence vector

Then the total least square solution of X is

$v=\frac{1}{\sqrt{2}{v}_{{\mathrm{\λ}}_0}[K+1]}\left[v_{{\mathrm{\λ}}_0}\right[1\left]v_{{\mathrm{\λ}}_0}\right[2]\·\·\·v_{{\mathrm{\λ}}_0}[K\left]\sqrt{2}{v}_{{\mathrm{\λ}}_0}\right[K+1\left]\right]---\left(9\right)$

Known by the Chebyshev polynomial expression,

cos(2ω)＝2(cos(ω)) ²-1（10）

The frequency of any in homogeneous solution v and 2 unknown signalings all satisfies equation

2cos(2ω _k)v ₁+2cos(ω _k)v ₂+v ₃＝0（11）

Then the frequency estimation of unknown signaling is

${\widehat{\mathrm{\ω}}}_1={\cos }^{-1}\left(\frac{v_2+\sqrt{v_2^2-4v_1{v}_3+8v_1^2}}{4v_1}\right)$

${\widehat{\mathrm{\ω}}}_2={\cos }^{-1}\left(\frac{v_2-\sqrt{v_2^2-4v_1{v}_3+8v_1^2}}{4v_1}\right)---\left(12\right)$

After signal frequency estimates, utilize least square method (LSM) can estimate amplitude and phase place.Formula (1) can be write as

x[n]＝a ₁cos(ω ₁n)cos(φ ₁)-a ₁sin(ω ₁n)sin(φ ₁)

+a ₂cos(ω ₂n)cos(φ ₂)-a ₂sin(ω ₂n)sin(φ ₂)+e(n)（13）

If y ₁=a ₁Cos (φ ₁), y ₂=a ₁Sin (φ ₁), y ₃=a ₂Cos (φ ₂), y ₄=a ₂Sin (φ ₂), treat that estimated parameter is y ₁, y ₂, y ₃, y ₄, according to each time observed reading x[n] (n=1,2 ... N), can draw the system of linear equations (14) of Ay=x, ask its least square solution can pick out each parameter

$\left[\begin{array}{cccc}\cos \left({\mathrm{\ω}}_1{t}_1\right)& -\sin \left({\mathrm{\ω}}_1{t}_1\right)& \cos \left({\mathrm{\ω}}_2{t}_1\right)& -\sin \left({\mathrm{\ω}}_2{t}_1\right)\\ \cos \left({\mathrm{\ω}}_1{t}_2\right)& -\sin \left({\mathrm{\ω}}_1{t}_2\right)& \cos \left({\mathrm{\ω}}_2{t}_2\right)& -\sin \left({\mathrm{\ω}}_2{t}_2\right)\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \cos \left({\mathrm{\ω}}_1{t}_N\right)& -\sin \left({\mathrm{\ω}}_1{t}_N\right)& \cos \left({\mathrm{\ω}}_2{t}_N\right)& -\sin \left({\mathrm{\ω}}_2{t}_N\right)\end{array}\right]\left[\begin{array}{c}{y}_1\\ y_2\\ y_3\\ y_4\end{array}\right]=\left[\begin{array}{c}x\left[1\right]\\ x\left[2\right]\\ \·\\ \·\\ \·\\ x\left[N\right]\end{array}\right]---\left(14\right)$

Then the estimated value of amplitude and phase place is

$a_1=\sqrt{y_1^2+y_2^2},$ $a_2=\sqrt{y_3^2+y_4^2}$

φ ₁＝arctan(y ₂/y ₁),φ ₂＝arctan(y ₄/y ₃) （15）

The below extracts with certain Dual Rotor System with Little Rotation Speed Difference (the inner and outer rotors frequency is respectively 28.283Hz and 28.1Hz) unbalanced signal specific implementation method of the present invention is described as an example.

Process flow diagram divided for 5 steps finished as shown in Figure 1, and the concrete operations in per step are as follows:

1, getting the sampling period is 1000Hz, and equal interval sampling obtains 1000 data;

2, get quality factor q=20, centre frequency is 28.1Hz, adopts digital band pass tracking filter algorithm to carry out filtering operation;

3, from the 401st of data after the filtering, per 16 data extract one, extract altogether 20 data as sample;

4, with formula (2)～(12) in the data substitution above-mentioned IP M method, obtain the frequencies omega of two unbalanced signals ₁, ω ₂

5, with ω ₁, ω ₂As given value, adopt 1000 data in (1) as sample, according to formula (14) structure system of linear equations, obtain y by the least square solution of system of equations ₁~ y ₄, obtain at last amplitude and the phase place of two signals according to formula (15).

The constructive simulation signal is as follows:

x(t)＝a ₁cos(2πf ₁t+φ ₁)+a ₂cos(2πf ₂t+φ ₂)+e(n)（16）

F wherein ₁=28.1, f ₂=28.283(be rotating speed be respectively 1566 and 1577rpm), a ₁=1, a ₂=0.5, φ ₁=30, φ ₂=230, e (n) is white noise signal.

When signal to noise ratio snr=50, shown in 10 measurement result following tables 1.

Table 1 is 10 measurement results when SNR=50

Sequence number	f 1	f 2	a 1	a 2	ω 1	φ 2
							1	26.10042594	26.28008422	1.0026	0.5069	30.3252	50.7706
2	26.09998165	26.2807692	1.0004	0.5051	30.3741	50.8167
							3	26.09995846	26.28186484	1.002	0.5032	30.1599	50.2225
4	26.10016344	26.28014628	1.006	0.5084	30.3457	50.405
							5	26.10058556	26.2820633	1.0049	0.5045	30.0185	49.8685
6	26.10121111	26.27879443	1.0137	0.5157	30.3361	50.1702
							7	26.10032799	26.28130369	1.0049	0.5057	30.1575	50.1219
8	26.10048734	26.28268253	1.0029	0.5025	29.9647	49.8353
							9	26.10057793	26.28099092	1.0064	0.5075	30.1613	50.0876
10	26.10033701	26.28045194	1.0062	0.5078	30.2561	50.2589

Claims (1)

1. the unbalanced signal extracting method of a Dual Rotor System with Little Rotation Speed Difference is characterized in that, may further comprise the steps:

A1 makes that the sampling period is T, and the vibration signal take T as the cycle to described Dual Rotor System with Little Rotation Speed Difference carries out equal interval sampling, obtains one group of data sequence { D (kT) }, k=1, and 2 ... K, K are a positive integer;

A2, getting centre frequency is f ₀, quality factor q＞2 adopt the digital band pass tracking filter that { D (kT) } carried out filtering, obtain { y (kT) }, k=1, and 2 ... K, wherein, the pulsed transfer function of digital band pass tracking filter is

$H\left(z\right)=k\frac{1-z^{-2}}{1+a_1{z}^{-1}+a_2{z}^{-2}},$ Each coefficient in the formula

$k=\frac{\tan \left(\frac{f_{0}T}{2}\right)}{Q+\tan \left(\frac{f_{0}T}{2}\right)+Q{\tan }^2\left(\frac{f_{0}T}{2}\right)}$

$a_1=\frac{2Q{\tan }^2\left(\frac{f_{0}T}{2}\right)-2Q}{Q+\tan \left(\frac{f_{0}T}{2}\right)+Q{\tan }^2\left(\frac{f_{0}T}{2}\right)}$

$a_2=\frac{Q-\tan \left(\frac{f_{0}T}{2}\right)+Q{\tan }^2\left(\frac{f_{0}T}{2}\right)}{Q+\tan \left(\frac{f_{0}T}{2}\right)+Q{\tan }^2\left(\frac{f_{0}T}{2}\right)};$

A3 with a default rule, makes data to { y (kT) } and extracts, the acquisition data sequence x[n] }, n=0 in the formula, 1 ... N-1 is time series, and N is a positive integer;

A4 establishes { x[n] } and contains consisting of of noisy sinusoidal signal

x[n]＝a ₁cos(ω ₁n+φ ₁)+a ₂cos(ω ₂n+φ ₂)+e(n)(1)

N=0 in the formula, 1 ... N-1 is time series, a ₁, a ₂Be unbalanced signal amplitude, ω ₁, ω ₂Be unbalanced signal frequency and φ ₁, φ ₂Be the unbalanced signal phase place,

The capable difference equation matrix of structure N-4 is

X[n]＝[x[n]+x[n+4],x[n+1]+x[n+3],x[n+2]]（2）

The Prony method shows, the homogeneous solution of above-mentioned equation has been specified e ^{I ω}In 4 rank multinomial coefficients, because sampled data all is real number, e ^{I ω}In 4 rank polynomial expressions can be rewritten into formula about cos (2 ω), thereby can solve with the Chebyshev polynomial expression,

By matrix X structure variance matrix G,

$G=\frac{1}{M}\left[\begin{array}{ccc}\mathrm{\Σ}{X}_{n,1}^2& \mathrm{\Σ}{X}_{n,1}{X}_{n,2}& \sqrt{2}\mathrm{\Σ}{X}_{n,1}{X}_{n,3}\\ \mathrm{\Σ}{X}_{n,1}{X}_{n,2}& \mathrm{\Σ}{X}_{n,2}^2& \sqrt{2}\mathrm{\Σ}{X}_{n,2}{X}_{n,3}\\ \sqrt{2}\mathrm{\Σ}{X}_{n,1}{X}_{n,3}& \sqrt{2}\mathrm{\Σ}{X}_{n,2}{X}_{n,3}& 2\mathrm{\Σ}{X}_{n,3}\end{array}\right]---\left(3\right)$

M=N-4 in the formula, all summations are limited to up and down from 0 to M-1, and the minimal eigenvalue of variance matrix G is by its secular equation polynomial expression

c ₁λ ³+c ₂λ ²+c ₃λ+c ₄＝0（4）

Coefficient ask for, this coefficient c is

$c=\left[\begin{array}{c}1\\ -G_{\mathrm{1,1}}-G_{\mathrm{2,2}}-G_{\mathrm{3,3}}\\ G_{\mathrm{1,1}}{G}_{\mathrm{2,2}}+G_{\mathrm{1,1}}{G}_{\mathrm{3,3}}+G_{\mathrm{2,2}}{G}_{\mathrm{3,3}}-G_{\mathrm{1,2}}^2-G_{\mathrm{1,3}}^2-G_{\mathrm{2,3}}^2\\ -G_{\mathrm{1,1}}{G}_{\mathrm{2,2}}{G}_{\mathrm{3,3}}-2G_{\mathrm{1,2}}{G}_{\mathrm{1,3}}{G}_{\mathrm{2,3}}+G_{\mathrm{1,1}}{G}_{\mathrm{2,3}}^2+G_{\mathrm{2,2}}{G}_{\mathrm{1,3}}^2+G_{\mathrm{3,3}}{G}_{\mathrm{1,2}}^2\end{array}\right]---\left(5\right)$

Then minimal eigenvalue is ${\mathrm{\λ}}_0=-c_2/3-(S+T)/2+i\sqrt{3}(S-T)/2---\left(6\right)$

Wherein S, T computing formula are,

$S=\sqrt[3]{R+\sqrt{Q^3+R^2}}$

$T=\sqrt[3]{R-\sqrt{Q^3+R^2}}$

$Q=(3c_3-c_2^2)/9$

$R=(9c_2{c}_3-27c_4-2c_2^3)/54---\left(7\right)$

Minimal eigenvalue λ ₀The characteristic of correspondence vector is

$v=\left[\begin{array}{c}-{\mathrm{\λ}}_0^2+(G_{\mathrm{2,2}}+G_{\mathrm{3,3}}){\mathrm{\λ}}_0-G_{\mathrm{2,2}}{G}_{\mathrm{3,3}}+G_{\mathrm{2,3}}^2\\ G_{\mathrm{1,2}}{G}_{\mathrm{3,3}}-G_{\mathrm{1,3}}{G}_{\mathrm{2,3}}-G_{\mathrm{1,2}}{\mathrm{\λ}}_0\\ \sqrt{2}(G_{\mathrm{1,3}}{G}_{\mathrm{2,2}}-G_{\mathrm{1,2}}{G}_{\mathrm{2,3}}-G_{\mathrm{1,3}}{\mathrm{\λ}}_0)\end{array}\right]---\left(8\right)$

Obtain the minimal eigenvalue λ of G matrix by matrix decomposition ₀With its characteristic of correspondence vector Then the total least square solution of matrix X is

$v=\frac{1}{\sqrt{2}{v}_{{\mathrm{\λ}}_0}[K+1]}\left[v_{{\mathrm{\λ}}_0}\right[1\left]v_{{\mathrm{\λ}}_0}\right[2]\·\·\·v_{{\mathrm{\λ}}_0}[K\left]\sqrt{2}{v}_{{\mathrm{\λ}}_0}\right[K+1\left]\right]---\left(9\right)$

Known cos (2 ω)=2 (cos (ω)) by the Chebyshev polynomial expression ²-1(10)

The frequency of any in homogeneous solution v and 2 unknown signalings all satisfies equation

2cos(2ω _k)v ₁+2cos(ω _k)v ₂+v ₃＝0（11）

Then the frequency estimation of unknown signaling is

${\widehat{\mathrm{\ω}}}_1={\cos }^{-1}\left(\frac{v_2+\sqrt{v_2^2-4v_1{v}_3+8v_1^2}}{4v_1}\right)$

${\widehat{\mathrm{\ω}}}_2={\cos }^{-1}\left(\frac{v_2-\sqrt{v_2^2-4v_1{v}_3+8v_1^2}}{4v_1}\right)---\left(12\right)$

Obtain the frequencies omega of unbalanced signal in the Dual Rotor System with Little Rotation Speed Difference vibration signal ₁And ω ₂

A5 utilizes least square method to estimate ω ₁And ω ₂Amplitude and phase place, formula (1) can be write as

x[n]＝a ₁cos(ω ₁n)cos(φ ₁)-a ₁sin(ω ₁n)sin(φ ₁)

+a ₂cos(ω ₂n)cos(φ ₂)-a ₂sin(ω ₂n)sin(φ ₂)+e(n)（13）

If y ₁=a ₁Cos (φ ₁), y ₂=a ₁Sin (φ ₁), y ₃=a ₂Cos (φ ₂), y ₄=a ₂Sin (φ ₂), treat that estimated parameter is y ₁, y ₂, y ₃, y ₄, according to each time observed reading x[n] (n=1,2 ... N), can draw the system of linear equations (14) of Ay=x, ask its least square solution can pick out each parameter,

$\left[\begin{array}{cccc}\cos \left({\mathrm{\ω}}_1{t}_1\right)& -\sin \left({\mathrm{\ω}}_1{t}_1\right)& \cos \left({\mathrm{\ω}}_2{t}_1\right)& -\sin \left({\mathrm{\ω}}_2{t}_1\right)\\ \cos \left({\mathrm{\ω}}_1{t}_2\right)& -\sin \left({\mathrm{\ω}}_1{t}_2\right)& \cos \left({\mathrm{\ω}}_2{t}_2\right)& -\sin \left({\mathrm{\ω}}_2{t}_2\right)\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \cos \left({\mathrm{\ω}}_1{t}_N\right)& -\sin \left({\mathrm{\ω}}_1{t}_N\right)& \cos \left({\mathrm{\ω}}_2{t}_N\right)& -\sin \left({\mathrm{\ω}}_2{t}_N\right)\end{array}\right]\left[\begin{array}{c}{y}_1\\ y_2\\ y_3\\ y_4\end{array}\right]=\left[\begin{array}{c}x\left[1\right]\\ x\left[2\right]\\ \·\\ \·\\ \·\\ x\left[N\right]\end{array}\right]---\left(14\right)$

The frequencies omega of unbalanced signal then ₁And ω ₂The estimated value of amplitude and phase place is

$a_1=\sqrt{y_1^2+y_2^2},$ $a_2=\sqrt{y_3^2+y_4^2}$

φ ₁＝arctan(y ₂/y ₁),φ ₂＝arctan(y ₄/y ₃)（15）。

Images