Rotor rub-impact identification method and system based on cyclic spectrum intrinsic orthogonal decomposition

Abstract

The invention relates to the field of rotary mechanical equipment vibration signal processing and rotor fault diagnosis, and provides a rotor rub-impact identification method and system based on cyclic spectrum eigen orthogonal decomposition, and the method comprises the steps: continuously collecting vibration signals of a continuously operating gas turbine and a combined cycle unit support bearing at equal intervals and equal length, and carrying out the preprocessing; a cyclic spectrum correlation mapping graph of the preprocessed vibration signals of the equal-length supporting bearing is obtained, a two-dimensional cyclic spectrum correlation mapping graph is obtained through amplitude normalization, and discrete estimation is carried out; constructing a snapshot matrix and performing intrinsic orthogonal decomposition; and determining the number of reconstructed modes according to the change condition of the modal energy ratio of each order, reconstructing an intrinsic orthogonal decomposition result, carrying out full-band integration along the spectral frequency direction to obtain a fault mode demodulation spectrum, and analyzing the correlation between the fault mode demodulation spectrum and the two-dimensional cyclic spectrum coherence mapping graph to determine the fault occurrence time. According to the method, the fault mode can be accurately extracted, the fault occurrence time can be identified, and the industrial field requirement of unit vibration monitoring in a continuous monitoring scene is met.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of vibration signal processing and rotor fault diagnosis of rotating machinery equipment, and particularly relates to a rotor rub-impact identification method and system based on cyclic spectrum eigen-orthogonal decomposition. BACKGROUND

[0002] In order to realize health management of core rotating equipment such as rotors in gas turbines and combined cycle units, ensure safe and stable long-term operation of the units, avoid unplanned shutdown of the units, improve economic benefits of operation of the units, and realize intelligent operation and maintenance of the units, the existing technology usually uses vibration acquisition equipment to continuously monitor vibration (including relative vibration of the rotor and vibration of the bearing seat) of the units in a continuous operation scenario, analyzes vibration data at different time periods, extracts weak fault features, identifies fault types, and determines fault occurrence time, which is crucial for realizing health management of the units.

[0003] Cyclic spectrum analysis can clearly show resonance information frequency bands caused by faults and discrete cyclic frequency structures caused by periodic motion mechanisms, and is an important tool for representing fault information of rotating parts such as rotors. Therefore, in the continuous monitoring scenario, how to analyze cyclic spectra of vibration signals at different time periods of the units, find fault modes from complex vibration data streams, retain main information of the fault modes, and determine fault occurrence time at the same time has become a problem that needs to be solved in state monitoring and fault diagnosis of gas turbines and combined cycle units in the power production industry. SUMMARY

[0004] The present application aims to at least solve one of the problems in the prior art, and provides a rotor rub-impact identification method and system based on cyclic spectrum eigen-orthogonal decomposition.

[0005] In one aspect of the present application, a rotor rub-impact identification method based on cyclic spectrum eigen-orthogonal decomposition is provided, comprising the following steps:

[0006] Step S1: continuously collecting vibration signals of support bearings of a continuously operated gas turbine and combined cycle unit at equal intervals and equal lengths, and pre-processing the vibration signals of the support bearings;

[0007] Step S2: obtaining a cyclic spectrum correlation mapping of the pre-processed vibration signals of the support bearings of equal lengths, performing amplitude normalization on the cyclic spectrum correlation mapping to obtain a two-dimensional cyclic spectrum coherence mapping, and performing discrete estimation on the two-dimensional cyclic spectrum coherence mapping to obtain discrete cyclic spectrum coherence mappings at different time periods;

[0008] Step S3: constructing a snapshot matrix according to the discrete cyclic spectrum coherence mappings at different time periods, and performing eigen-orthogonal decomposition on the snapshot matrix to determine proportions of energy of each order mode;

[0009] Step S4: determining the number of reconstructed modes according to the change of the proportion of the energy of each order mode, reconstructing the eigenvalue orthogonal decomposition result of the snapshot matrix according to the number of reconstructed modes, obtaining a cyclic spectrum coherence mapping of the fault mode, and integrating the cyclic spectrum coherence mapping of the fault mode in the spectral frequency direction to obtain a demodulation spectrum corresponding to the fault mode;

[0010] Step S5: analyzing the correlation between the demodulation spectrum corresponding to the fault mode and the two-dimensional cyclic spectrum coherence mapping to determine the occurrence time of the fault.

[0011] Optionally, the step S2 comprises the following steps:

[0012] Step S21: obtaining a cyclic spectrum correlation mapping CSC of the preprocessed equal-length support bearing vibration signal by using two-dimensional Fourier transform according to the following formula: x (α,f):

[0013] CSC x (α,f)=∫∫R x (t,τ)e -2πj(αt+fτ) dtdτ=∫∫E{x(t)x * (t+τ)}e -2πj(αt+fτ) dtdτ;

[0014] Wherein, x(t) represents the support bearing vibration signal at t time, * represents the conjugate operation, x * (t+τ) represents the conjugate operation on the support bearing vibration signal x(t+τ) at t+τ time, R x (t,τ) represents the two-dimensional time-varying autocorrelation function of the support bearing vibration signal, τ represents the time delay, E represents the mathematical expectation, e represents the natural constant, j represents the imaginary unit; α is the cyclic frequency, corresponding to the frequency domain mapping of time t; f is the spectral frequency, corresponding to the frequency domain mapping of time delay τ;

[0015] Step S22: amplitude normalizing the cyclic spectrum correlation mapping CSC x (α,f) according to the following formula to obtain a two-dimensional cyclic spectrum coherence mapping CSCoh x (α,f):

[0016]

[0017] Wherein, CSC x (0,f) represents the cyclic energy distribution of the support bearing vibration signal in the (0,f) plane, and CSC x (0,f+α) represents the cyclic energy distribution of the support bearing signal in the (0,f+α) plane.

[0018] Step S23: Apply the following formula to the two-dimensional cyclic spectrum coherence mapping map CSCoh x Discrete estimation of (α,f) yields:

[0019] S M×N (m,n)=CSCoh x (α m ,f n );

[0020] Among them, CSCoh x (α m ,f n ) represents the frequency α in the m-th cycle. m and the nth spectral frequency f n The cyclic spectral coherence value at point S M×N (m,n) represents a discrete cyclic spectral coherence mapping S of size M×N. M×N The elements in the m-th row and n-th column of the array are m = 1, 2, ..., M, n = 1, 2, ..., N, where M represents the number of cyclic frequency points and N represents the number of spectral frequency points.

[0021] Optionally, step S3 includes the following steps:

[0022] Step S31: Based on the discrete cyclic spectrum coherence mapping diagram at different times Construct a snapshot matrix X of size (M×N)×K, and represent the snapshot matrix X as follows:

[0023]

[0024] in, S represents the discrete cyclic spectrum coherence mapping at time 1, time 2, ..., time K, respectively. 1 (1,1), ..., S 1 (1, N) represent the elements in the first row and first column, ..., Nth column of the discrete cyclic spectral coherence map at time 1, respectively. 1 (2,1), ...,S 1 (2, N) represent the elements in the 1st, ..., Nth columns of the 2nd row of the discrete cyclic spectral coherence map at time 1, respectively, and S 1 (M,1), ...,S 1 (M, N) represent the elements in the M-th row and 1st column, ..., N-th column of the discrete cyclic spectral coherence map at time 1, respectively. 2 (1,1), ..., S 2 (1, N) represent the elements in the first row and first column, ..., Nth column of the discrete cyclic spectral coherence map at time 2, respectively, and S 2 (2,1), ...,S 2(2,1), …, S 2 (M,1), …, S 2 (M,N) represent the (M,1), …, S K (1,1), …, S K (1,N) represent the (1,1), …, S K (2,1), …, S K (2,N) represent the (2,1), …, S K (M,1), …, S K (M,N) represent the (M,1), …, S

[0025] Step S32: perform eigen-orthogonal decomposition on the snapshot matrix X based on singular value decomposition to obtain a decomposition modal matrix U and a modal energy value Λ according to the following formula:

[0026] X = U x Λ x V T ;

[0027] wherein the decomposition modal matrix U has a size of (M x N) x K and u i represents the i th modal obtained by eigen-orthogonal decomposition, and represents a real number column vector space with a dimension of M x N, i = 1, 2, …, K; Λ = diag(σ 1, σ 2, …, σ K ) is a diagonal matrix, representing modal energy; σ 1, σ 2, …, σ K represent the first-order modal energy, the second-order modal energy, …, and the K th-order modal energy, respectively, and σ 1 > σ 2 > … > σ K ; V represents a time coefficient, and T represents matrix transposition.

[0028] Step S33: calculate the specific proportion of each order of modal energy according to the following formula:

[0029]

[0030] wherein, represents the proportion of the i th-order modal energy σ i .

[0031] Optionally, the step S4 comprises the following steps:

[0032] Step S41: calculate a modal energy difference vector dσ, and the Jth element dσ J is represented as:

[0033] dσ J = σ J+1 - σ J , J = 1, 2, …, K-1;

[0034] wherein σ J represents a proportion of the Jth order modal energy, and σ J+1 represents a proportion of the (J+1)th order modal energy.

[0035] Step S42: determine the position of the largest modal energy change according to the difference value corresponding to each element in the modal energy difference vector dσ, and determine the reconstruction modal number Q according to the following formula:

[0036] Q = argmax{dσ};

[0037] wherein argmax{dσ} represents an index value corresponding to the maximum value of the elements in the modal energy difference vector dσ.

[0038] Step S43: reconstruct the decomposition modal matrix U according to the reconstruction modal number Q, to obtain a reconstruction vector convert the reconstruction vector into a two-dimensional matrix Y M×N as a cyclic spectral coherence mapping diagram of the fault modal, wherein:

[0039]

[0040] wherein, respectively represent the 1st, Mx1, M+1, Mx2, Mx(N-1)+1, and MxN elements in the reconstruction vector .

[0041] Step S44: perform full-band integration on the cyclic spectral coherence mapping diagram Y M×N of the fault modal in the spectral frequency direction according to the following formula, to obtain a demodulation spectrum corresponding to the fault modal:

[0042]

[0043] wherein FDS(α m ) represents a fault demodulation spectrum amplitude at the mth cyclic frequency α m , and Y(m, n) represents the mth row and nth column element in the cyclic spectral coherence mapping diagram Y M×N of the fault modal.

[0044] Optionally, the step S4 further comprises:

[0045] The demodulation spectrum corresponding to the fault mode is normalized to 0-1 according to the following formula:

[0046]

[0047] Wherein, NFDS(alpha m ) represents the normalized fault demodulation spectrum, min{FDS(alpha m )} represents the minimum value of FDS(alpha m ), and max{FDS(alpha m )} represents the maximum value of FDS(alpha m ).

[0048] Another aspect of the present application provides a rotor rub-impact identification system based on cyclic spectrum eigen-orthogonal decomposition, comprising:

[0049] The preprocessing module is configured to collect vibration signals of a support bearing of a continuously running gas turbine and combined cycle unit at equal intervals and equal lengths, and to preprocess the vibration signals of the support bearing;

[0050] The discrete estimation module is configured to obtain a cyclic spectrum correlation map of the preprocessed vibration signals of the support bearing of equal lengths, to perform amplitude normalization on the cyclic spectrum correlation map to obtain a two-dimensional cyclic spectrum coherence map, and to perform discrete estimation on the two-dimensional cyclic spectrum coherence map to obtain discrete cyclic spectrum coherence maps at different times;

[0051] The decomposition module is configured to construct a snapshot matrix according to the discrete cyclic spectrum coherence maps at different times, and to perform eigen-orthogonal decomposition on the snapshot matrix to determine the proportion of energy of each order mode;

[0052] The reconstruction module is configured to determine the number of reconstructed modes according to the proportion of energy of each order mode, to reconstruct the eigen-orthogonal decomposition result of the snapshot matrix according to the number of reconstructed modes, to obtain a cyclic spectrum coherence map of a fault mode, and to perform full-band integration on the cyclic spectrum coherence map of the fault mode in the spectral frequency direction to obtain a demodulation spectrum corresponding to the fault mode;

[0053] The analysis module is configured to analyze the correlation between the demodulation spectrum corresponding to the fault mode and the two-dimensional cyclic spectrum coherence map to determine the occurrence time of the fault.

[0054] Optionally, the discrete estimation module is specifically configured to:

[0055] The cyclic spectrum correlation map CSC(alpha,f) of the preprocessed vibration signals of the support bearing of equal lengths is obtained by two-dimensional Fourier transform according to the following formula: x

[0056] CSC​x (a,f) = ∫∫R x (t,τ)e -2πj(αt+fτ) dtdτ = ∫∫E{x(t)x * (t+τ)}e -2πj(αt+fτ) dtdτ;

[0057] wherein x(t) represents a support bearing vibration signal at t time, * represents a conjugate operation, x * (t+τ) represents a conjugate operation on the support bearing vibration signal x(t+τ) at t+τ time, R x (t,τ) represents a two-dimensional time-varying autocorrelation function of the support bearing vibration signal, τ represents a time delay, E represents a mathematical expectation, e represents a natural constant, j represents an imaginary unit; α is a cyclic frequency, corresponding to the frequency domain mapping of time t; f is a spectral frequency, corresponding to the frequency domain mapping of time delay τ;

[0058] The two-dimensional cyclic spectrum coherence mapping CSCoh x (α,f) is obtained by amplitude normalization of the cyclic spectrum correlation mapping CSC x (α,f) according to the following formula:

[0059]

[0060] wherein CSC x (0,f) represents a cyclic energy distribution of the support bearing vibration signal in the (0,f) plane, CSC x (0,f+α) represents a cyclic energy distribution of the support bearing vibration signal in the (0,f+α) plane;

[0061] The two-dimensional cyclic spectrum coherence mapping CSCoh x (α,f) is discretely estimated according to the following formula, and the following is obtained:

[0062] S M×N (m,n) = CSCoh x (α m ,d n );

[0063] wherein CSCoh x (α m ,f n ) represents a cyclic spectrum coherence value at the mth cyclic frequency α m and the nth spectral frequency f n , S M×n (m,n) represents an element in the mth row and the nth column of the discrete cyclic spectrum coherence mapping S M×N with a size of M×N, m=1,2,…,M, n=1,2,…,N, M represents a cyclic frequency point number, and N represents a spectral frequency point number.

[0064] Optionally, the decomposing module is specifically used for:

[0065] discrete cyclic spectrum coherence mapping diagrams at different time instants constructing a snapshot matrix X with a size of (M*N) * K, and representing the snapshot matrix X as:

[0066]

[0067] wherein, respectively represent the first column, the second column, …, and the Nth column element in the first row of the discrete cyclic spectrum coherence mapping diagram at the first time instant, S 1 respectively represent the first column, the second column, …, and the Nth column element in the first row of the discrete cyclic spectrum coherence mapping diagram at the first time instant, S 1 respectively represent the first column, the second column, …, and the Nth column element in the first row of the discrete cyclic spectrum coherence mapping diagram at the first time instant, S 1 respectively represent the first column, the second column, …, and the Nth column element in the first row of the discrete cyclic spectrum coherence mapping diagram at the first time instant, S 1 respectively represent the first column, the second column, …, and the Nth column element in the first row of the discrete cyclic spectrum coherence mapping diagram at the first time instant, S 1 respectively represent the first column, the second column, …, and the Nth column element in the first row of the discrete cyclic spectrum coherence mapping diagram at the first time instant, S 1 respectively represent the first column, the second column, …, and the Nth column element in the first row of the discrete cyclic spectrum coherence mapping diagram at the first time instant, S 2 respectively represent the first column, the second column, …, and the Nth column element in the first row of the discrete cyclic spectrum coherence mapping diagram at the first time instant, S 2 respectively represent the first column, the second column, …, and the Nth column element in the first row of the discrete cyclic spectrum coherence mapping diagram at the first time instant, S 2 respectively represent the first column, the second column, …, and the Nth column element in the first row of the discrete cyclic spectrum coherence mapping diagram at the first time instant, S 2 respectively represent the first column, the second column, …, and the Nth column element in the first row of the discrete cyclic spectrum coherence mapping diagram at the first time instant, S 2 respectively represent the first column, the second column, …, and the Nth column element in the first row of the discrete cyclic spectrum coherence mapping diagram at the first time instant, S 2 respectively represent the first column, the second column, …, and the Nth column element in the first row of the discrete cyclic spectrum coherence mapping diagram at the first time instant, S K respectively represent the first column, the second column, …, and the Nth column element in the first row of the discrete cyclic spectrum coherence mapping diagram at the first time instant, S K respectively represent the first column, the second column, …, and the Nth column element in the first row of the discrete cyclic spectrum coherence mapping diagram at the first time instant, S K respectively represent the first column, the second column, …, and the Nth column element in the first row of the discrete cyclic spectrum coherence mapping diagram at the first time instant, S K respectively represent the first column, the second column, …, and the Nth column element in the first row of the discrete cyclic spectrum coherence mapping diagram at the first time instant, S K respectively represent the first column, the second column, …, and the Nth column element in the first row of the discrete cyclic spectrum coherence mapping diagram at the first time instant, S K respectively represent the first column, the second column, …, and the Nth column element in the first row of the discrete cyclic spectrum coherence mapping diagram at the first time instant, S

[0068] The intrinsic orthogonal decomposition of the snapshot matrix X is performed based on singular value decomposition according to the following formula to obtain a decomposition modal matrix U and a modal energy value Λ:

[0069] X = U * Λ * V T ;

[0070] wherein the decomposition modal matrix U has a size of (M * N) * K and u i represents the i th modal obtained by the intrinsic orthogonal decomposition, and represents a real number column vector space with a dimension of M * N, i = 1, 2, …, K; Λ = diag (σ 1, σ 2, …, σ K ) is a diagonal matrix, representing modal energy; σ 1, σ 2, …, σ K respectively represent the 1 st order modal energy, the 2 nd order modal energy, …, the K th order modal energy, and σ 1 > σ 2 > … > σ K ; V represents a time coefficient, and T represents a matrix transpose;

[0071] According to the following formula, the specific proportion of each order modal energy is calculated:

[0072]

[0073] wherein, represents the proportion of the i th order modal energy σ i .

[0074] Optionally, the reconstruction module is specifically configured to:

[0075] calculate a modal energy difference vector dσ, and the J th element dσ J in the modal energy difference vector dσ is represented as:

[0076] dσ J = σ J+1 - σ J , J = 1, 2, …, K-1;

[0077] wherein, σ J represents the proportion of the J th order modal energy, and σ J+1 represents the proportion of the J+1 th order modal energy;

[0078] According to the difference value corresponding to each element in the modal energy difference vector dσ, the position of the largest modal energy change is determined, and the reconstruction modal number Q is determined according to the following formula:

[0079] Q = argmax {dσ};

[0080] wherein, argmax {dσ} represents an index value corresponding to the maximum value of the elements in the modal energy difference vector dσ;

[0081] According to the number Q of reconstruction modalities, the decomposition modality matrix U is reconstructed to obtain a reconstruction vector The reconstruction vector Y is converted into a two-dimensional matrix Y M×N The cyclic spectrum coherence mapping of the fault modality, wherein:

[0082]

[0083] wherein, respectively represent the 1st, Mx1, M+1, Mx2, Mx(N-1)+1, MxN elements in the reconstruction vector Y

[0084] The cyclic spectrum coherence mapping Y of the fault modality is integrated along the spectral frequency direction M×N to obtain a demodulation spectrum corresponding to the fault modality,

[0085]

[0086] wherein, FDS(a m ) represents the fault demodulation spectrum amplitude at the mth cyclic frequency a m , and Y(m, n) represents the mth row and n th column element in the cyclic spectrum coherence mapping Y of the fault modality. M×N

[0087] Optionally, the reconstruction module is further used for:

[0088] According to the following formula, the demodulation spectrum corresponding to the fault modality is normalized to 0-1:

[0089]

[0090] wherein, NFDS(a m ) represents the normalized fault demodulation spectrum, min{FDS(a m )} represents the minimum value of FDS(a m ), and max{FDS(a m )} represents the maximum value of FDS(a m ).

[0091] ​​​The application can extract the fault mode from the continuously collected vibration signal by adopting the cyclic spectrum analysis, equally spaced and equally long continuous acquisition of the vibration signal of the support bearing of the continuous operation unit, and the modal decomposition of the cyclic spectrum by the intrinsic orthogonal decomposition, and reconstructing the fault mode according to the distribution of the modal energy; the application can identify the fault occurrence time in the continuous monitoring scene of the gas turbine and the combined cycle unit by analyzing the correlation between the fault mode obtained by the intrinsic orthogonal decomposition reconstruction and the cyclic spectrum of the vibration signal obtained by the continuous monitoring; the application can be verified by the identification and occurrence time determination of the rotor rub-impact fault in the continuous monitoring scene, fully proving the unit state monitoring and fault diagnosis capability of the application in the continuous monitoring scene, and fully meeting the industrial field demand of the rotor vibration monitoring of the gas turbine and the combined cycle unit in the continuous monitoring scene. BRIEF DESCRIPTION OF DRAWINGS

[0092] One or more embodiments are illustrated by way of example in the figures that are part of this disclosure and which are illustrative, but not restrictive, of the embodiments, wherein elements having the same reference number designates like elements throughout the various figures, and the figures are not necessarily to scale.

[0093] Figure 1 A flow chart of a rotor rub-impact identification method based on the intrinsic orthogonal decomposition of the cyclic spectrum is provided for an embodiment of the application;

[0094] Figure 2 A structural schematic diagram of a test bench is provided for another embodiment of the application;

[0095] Figure 3 A schematic diagram of equally spaced and equally long continuously collected vibration signals of the support bearing is provided for another embodiment of the application;

[0096] Figure 4 A schematic diagram of equally long vibration signals collected at different times is provided for another embodiment of the application;

[0097] Figure 5 A two-dimensional cyclic spectrum coherence mapping diagram of vibration signals at different times is provided for another embodiment of the application;

[0098] Figure 6 A schematic diagram of all modes obtained based on the intrinsic orthogonal decomposition of the cyclic spectrum is provided for another embodiment of the application;

[0099] Figure 7 A distribution diagram of the proportion of all modal energies obtained based on the intrinsic orthogonal decomposition of the cyclic spectrum is provided for another embodiment of the application;

[0100] Figure 8A cyclic spectrum coherence mapping diagram of a reconstructed fault mode provided by another embodiment of the present application;

[0101] Figure 9 A schematic diagram of a demodulation spectrum corresponding to a fault mode obtained by integrating a fault mode in a spectral frequency direction provided by another embodiment of the present application;

[0102] Figure 10 A schematic diagram of a correlation coefficient variation of a fault mode and a two-dimensional cyclic spectrum coherence mapping diagram at different times provided by another embodiment of the present application. DETAILED DESCRIPTION

[0103] In order to make the objectives, technical solutions and advantages of the embodiments of the present application clearer, the embodiments of the present application will be described in detail below with reference to the drawings. However, those skilled in the art can understand that, in the embodiments of the present application, many technical details are proposed in order to make the readers better understand the present application. However, the technical solutions claimed by the present application can be implemented even without these technical details and various changes and modifications based on the following embodiments. The division of the following embodiments is for the convenience of description, and should not constitute any limitation on the specific implementation of the present application, and the embodiments can be combined and referenced with each other without contradiction.

[0104] One embodiment of the present application relates to a rotor rub-impact identification method based on cyclic spectrum eigen-orthogonal decomposition, and a flowchart thereof is shown as Figure 1 The method can identify a rotor rub-impact fault and its occurrence time from continuously monitored vibration signals of a support bearing of a gas turbine and a combined cycle unit.

[0105] Step S1: The vibration signals of the support bearing of the continuously running gas turbine and combined cycle unit are continuously collected at equal intervals and equal lengths, and the vibration signals of the support bearing are preprocessed.

[0106] Specifically, step S1 can continuously collect the vibration signals of the support bearing of the continuously running gas turbine and combined cycle unit at equal intervals and equal lengths by using an acceleration sensor, and perform preprocessing operations such as high-frequency sampling, mean removal, and pre-whitening on the vibration signals of the support bearing, so as to ensure that the subsequent analysis is not disturbed by the trend items and harmonic components generated by the normal vibration of the rotor.

[0107] Step S2: Obtain a cyclic spectrum correlation mapping diagram of the preprocessed vibration signals of the support bearing at equal lengths, perform amplitude normalization on the cyclic spectrum correlation mapping diagram, obtain a two-dimensional cyclic spectrum coherence mapping diagram, and perform discrete estimation on the two-dimensional cyclic spectrum coherence mapping diagram to obtain discrete cyclic spectrum coherence mapping diagrams at different times.

[0108] Exemplarily, step S2 includes the following steps S21 to S23.

[0109] Step S21: According to the following formula, a two-dimensional Fourier transform is used to obtain a cyclic spectrum correlation mapping Csc of the preprocessed equal-length support bearing vibration signal x (α,f):

[0110] Csc x (α,f) = ∫∫R x (t,τ)e -2πj(αt+fτ) j(αf+2πft) dtdτ = ∫∫E{x(t)x * (t+τ)}e -2πj(αt+fτ) j(αf+2πft) dtdτ;

[0111] Wherein, x(t) represents the support bearing vibration signal at t time, * represents the conjugate operation, x * (t+τ) represents the conjugate operation of the support bearing vibration signal x(t+τ) at t+τ time, R x (t,τ) represents the two-dimensional time-varying autocorrelation function of the support bearing vibration signal, τ represents time delay, E represents mathematical expectation, e represents natural constant, j represents imaginary unit; α is the cyclic frequency, corresponding to the frequency domain mapping of time t; f is the spectral frequency, corresponding to the frequency domain mapping of time delay τ.

[0112] Step S22: According to the following formula, the amplitude of the cyclic spectrum correlation mapping Csc x (α,f) is normalized to obtain a two-dimensional cyclic spectrum coherence mapping Cscoh x (α,f):

[0113]

[0114] Wherein, Csc x (0,f) represents the cyclic energy distribution of the support bearing vibration signal in the (0,f) plane, Csc x (0,f+α) represents the cyclic energy distribution of the support bearing vibration signal in the (0,f+α) plane.

[0115] Step S23: According to the following formula, the two-dimensional cyclic spectrum coherence mapping Cscoh x (α,f) is discretely estimated to obtain:

[0116] S M×N (m,n) = Cscoh x (α m ,f n );

[0117] Wherein, Cscoh x (α m ,f n ) represents the cyclic spectrum coherence of the support bearing vibration signal in the mth cyclic frequency α mand the n th spectral frequency f n the cyclic spectral coherence value at the (m, n) th element in the discrete cyclic spectral coherence map S M×N (m, n) represents the m th row, n th column element in the discrete cyclic spectral coherence map S M×N (m, n) represents the m th row, n th column element in the discrete cyclic spectral coherence map S

[0118] Step S3: Constructing a snapshot matrix according to the discrete cyclic spectral coherence maps at different time instants, and performing eigen-orthogonal decomposition on the snapshot matrix to determine the proportion of energy of each order mode.

[0119] For example, step S3 includes the following steps S31 to S33.

[0120] Step S31: Constructing a snapshot matrix X with a size of (M x N) x K according to the discrete cyclic spectral coherence maps at different time instants Step S31: Constructing a snapshot matrix X with a size of (M x N) x K according to the discrete cyclic spectral coherence maps at different time instants

[0121]

[0122] wherein, respectively represent the discrete cyclic spectral coherence maps at the 1 st time instant, the 2 nd time instant, …, and the K th time instant, respectively. 1 (1, 1), …, S 1 (1, N) represent the 1 st column, …, and the N th column elements in the 1 st row of the discrete cyclic spectral coherence map at the 1 st time instant, respectively. 1 (2, 1), …, S 1 (2, N) represent the 1 st column, …, and the N th column elements in the 2 nd row of the discrete cyclic spectral coherence map at the 1 st time instant, respectively. 1 (M, 1), …, S 1 (M, N) represent the 1 st column, …, and the N th column elements in the M th row of the discrete cyclic spectral coherence map at the 1 st time instant, respectively. 2 (1, 1), …, S 2 (1, N) represent the 1 st column, …, and the N th column elements in the 1 st row of the discrete cyclic spectral coherence map at the 2 nd time instant, respectively. 2 (2, 1), …, S 2 (2, N) represent the 1 st column, …, and the N th column elements in the 2 nd row of the discrete cyclic spectral coherence map at the 2 nd time instant, respectively. 2 (M, 1), …, S 2 (M, N) represent the 1 st column, …, and the N th column elements in the M th row of the discrete cyclic spectral coherence map at the 2 nd time instant, respectively. K (1, 1), …, S K(1,1), …, S K (2,1), …, S K (2,1), …, S K (2,1), …, S K (2,1), …, S

[0123] Step S32: Perform eigen-orthogonal decomposition on the snapshot matrix X based on singular value decomposition to obtain a decomposition modal matrix U and a modal energy value Λ according to the following formula:

[0124] X = U x Λ x V T ;

[0125] wherein the size of the decomposition modal matrix U is (M x N) x K and u i represents the i th modal obtained by eigen-orthogonal decomposition, and represents a real number column vector space with a dimension of M x N, i = 1, 2, …, K; Λ = diag(σ 1, σ 2, …, σ K ) is a diagonal matrix, representing modal energy; σ 1, σ 2, …, σ K represent the first-order modal energy, the second-order modal energy, …, and the K th-order modal energy, respectively, and σ 1> σ 2> …> σ K ; V represents a time coefficient, and T represents matrix transposition.

[0126] Step S33: Calculate the specific proportion of each order modal energy according to the following formula:

[0127]

[0128] wherein, represents the proportion of the i th-order modal energy σ i .

[0129] Step S4: Determine the number of reconstructed modes according to the proportion change of each order modal energy, and reconstruct the eigen-orthogonal decomposition result of the snapshot matrix according to the number of reconstructed modes to obtain a cyclic spectral coherence mapping diagram of a fault mode, and integrate the cyclic spectral coherence mapping diagram of the fault mode in the spectral frequency direction to obtain a demodulation spectrum corresponding to the fault mode, so as to identify a rotor fault according to the demodulation spectrum corresponding to the fault mode.

[0130] For example, step S4 includes the following steps S41 to S44.

[0131] Step S41: analyze the change of the modal energy proportion: calculate the modal energy difference vector dσ, and the Jth element dσ J is expressed as:

[0132] dσ J = σ J+1 - σ J , J = 1, 2, …, K-1;

[0133] wherein σ J represents the proportion of the Jth modal energy, and σ J+1 represents the proportion of the (J+1)th modal energy.

[0134] Step S42: determine the number of reconstructed modes: determine the position of the largest modal energy change according to the difference value corresponding to each element in the modal energy difference vector dσ, and determine the number of reconstructed modes Q according to the following formula:

[0135] Q = argmax {dσ};

[0136] wherein argmax {dσ} represents the index value corresponding to the maximum value of the elements in the modal energy difference vector dσ.

[0137] Step S43: reconstruct the decomposition modal matrix U according to the number of reconstructed modes Q to obtain the reconstructed vector convert the reconstructed vector into a two-dimensional matrix Y M×N , which is a cyclic spectral coherence mapping of the fault mode, wherein:

[0138]

[0139] wherein, respectively represent the 1st, Mx1, M+1, Mx2, Mx(N-1)+1, and MxN elements in the reconstructed vector .

[0140] Step S44: perform full-band integration on the cyclic spectral coherence mapping Y M×N of the fault mode in the spectral frequency direction according to the following formula to obtain the demodulation spectrum corresponding to the fault mode:

[0141]

[0142] wherein FDS (α m ) represents the fault demodulation spectrum amplitude at the mth cyclic frequency α m , and Y (m, n) represents the element in the mth row and the nth column in the cyclic spectral coherence mapping Y M×N of the fault mode.

[0143] For example, step S4 further includes: performing 0-1 normalization on the demodulation spectrum corresponding to the fault mode according to the following formula:

[0144]

[0145] Among them, NFDS(α) m ) represents the normalized fault demodulation spectrum, min{FDS(α) m )} represents FDS(α) m The minimum value of ), max{FDS(α) m )} represents FDS(α) m The maximum value of ).

[0146] Step S5: Analyze the correlation between the demodulation spectrum corresponding to the fault mode and the coherence mapping diagram of the two-dimensional cyclic spectrum to determine the time of fault occurrence.

[0147] Specifically, step S5 can use correlation coefficients, Euclidean distances, etc. to analyze the correlation between the demodulation spectrum corresponding to the fault mode and the coherence mapping diagram of the two-dimensional cyclic spectrum, thereby determining the time of fault occurrence.

[0148] In one specific embodiment, the rotor rubbing identification method based on cyclic spectrum intrinsic orthogonal decomposition provided by the above embodiments of the present invention is applied to the rubbing fault diagnosis of tie rod rotors for verification.

[0149] Take a test bench for a certain tie-rod rotor system as an example. Figure 2 As shown, the test bench comprises a drive motor 1, a support bearing 2, a bladeless tie rod disc 4, tie rod bolts 5, a bladed tie rod disc 6, a composite bearing 7, several couplings 8, a speed and torque measurement module 9, a magnetic powder brake 10, and a feed screw 11. The drive motor uses a rated voltage of 380V, a rated power of 37kW, a rated speed of 3000rpm, and a rated torque of 119N·m. The support bearing 2 and the composite bearing 7 bear the radial and axial loads, respectively. The rotor consists of two parts: the bladeless tie rod disc 4 and the bladed tie rod disc 6. Both are axially connected to the tie rod bolts 5, forming the entire rotor system. They are connected to the drive motor 1 and the magnetic powder brake 10 via several couplings 8. The magnetic powder brake provides the necessary torque load to the entire device. An acceleration sensor 3 is located on top of the support bearing 2 and is used to continuously collect vibration signals during the experiment. By analyzing the vibration signals, rotor rubbing fault diagnosis is achieved, and the time of fault occurrence is determined.

[0150] In one specific embodiment, to simulate the rubbing caused by insufficient radial clearance between the blades and the casing, soft copper nails were screwed into the radial position of the wheel disc. During the experiment, the radial clearance between the two could be adjusted by laterally moving the casing with a screw, simulating the mutual friction between the rotor blades and the casing. Therefore, during the lateral movement of the screw, the copper nails gradually contacted the casing, and at a certain moment, the copper nails and the casing began to rub against each other. During the vibration data acquisition process of the bearing housing, the acceleration signal acquisition frequency was set to 10240Hz, and vibration data were acquired at equal intervals for a duration of 1 second. During the experiment, the drive motor speed was set to 3000rpm, and the load was set to 20% of the full load of the magnetic powder brake. The copper nails rubbed against the casing once for every revolution of the rotor; therefore, the characteristic frequency of this type of fault was 50Hz. Continuous data acquisition was performed by laterally moving the screw at equal intervals and with equal lengths. Figure 3 The acceleration signal shown is obtained by preprocessing equal-length time-domain signals collected at different times, including mean removal and pre-whitening. Figure 4 The vibration signals shown indicate that it is difficult to detect them from... Figure 4 Impact components were found in the time-domain signal.

[0151] To illustrate the superiority of the rotor rubbing identification method based on cyclic spectrum eigenorthogonal decomposition provided in the above embodiments of the present invention, firstly, the two-dimensional cyclic spectrum coherence matrix of the whitened vibration signal is calculated to obtain... Figure 5 The two-dimensional cyclic spectrum coherence mapping diagram is shown. Next, using the rotor rubbing identification method based on cyclic spectrum eigenorthogonal decomposition provided in the above embodiments of the present invention, the decomposed modes are as follows: Figure 6 As shown, the modal energy distribution is as follows: Figure 7 As shown. According to the rotor rubbing identification method based on cyclic spectrum intrinsic orthogonal decomposition provided by the above embodiments of the present invention, reconstruction is performed based on the mode with the largest modal energy change. The cyclic spectrum coherence mapping diagram of the fault mode is obtained by reconstructing using the first decomposed mode, as shown in the figure. Figure 8 As shown, along the spectral frequency direction... Figure 8 The cyclic coherence map of the fault mode shown is integrated across the entire frequency band to obtain the demodulation spectrum corresponding to the fault mode, as shown below. Figure 9 As shown, from Figure 9 As can be seen, the characteristic frequencies of 50Hz and its harmonics are clearly identified in the weighted demodulation spectrum. Next, the fault modes and... Figure 4 The correlation relationship of the coherent mapping of each mode spectrum of the vibration signal is analyzed as follows: Figure 10 As shown, from Figure 10As can be seen from the mode numbered #15, there is a strong correlation with the fault mode, and the correlation coefficient is close to 1.0, indicating that the fault has occurred at this moment. This embodiment not only verifies the effectiveness of the rotor rub-impact identification method based on cyclic spectrum eigen-orthogonal decomposition provided in the above-mentioned embodiments of the application in extracting rotor rub-impact fault features, but also shows the advantage of the rotor rub-impact identification method based on cyclic spectrum eigen-orthogonal decomposition in determining the fault occurrence time.

[0152] The rotor rub-impact identification method based on cyclic spectrum eigen-orthogonal decomposition provided in the embodiments of the application uses cyclic spectrum analysis of vibration signals collected at different times, and uses eigen-orthogonal decomposition to decompose the cyclic spectrum mapping diagrams at different times, reconstructs the fault mode according to the distribution of the mode energy change, identifies the fault type according to the reconstructed mode, and further analyzes the correlation between the fault mode and the cyclic spectrum at different times, so as to identify the fault occurrence time in the continuous monitoring of the gas turbine and combined cycle unit. The rotor rub-impact identification method based on cyclic spectrum eigen-orthogonal decomposition provided in the embodiments of the application not only can accurately extract the fault mode for fault diagnosis of the unit, but also can identify the fault occurrence time of the unit in continuous monitoring, and can meet the industrial field demand for rotor vibration monitoring of the gas turbine and combined cycle unit in the continuous monitoring scene.

[0153] Another embodiment of the application relates to a rotor rub-impact identification system based on cyclic spectrum eigen-orthogonal decomposition, which comprises a preprocessing module, a discrete estimation module, a decomposition module, a reconstruction module and an analysis module.

[0154] The preprocessing module is used for continuously collecting support bearing vibration signals of the continuously running gas turbine and combined cycle unit at equal intervals and equal lengths, and pre-processing the support bearing vibration signals.

[0155] The discrete estimation module is used for obtaining a cyclic spectrum correlation mapping diagram of the pre-processed support bearing vibration signals of equal length, performing amplitude normalization on the cyclic spectrum correlation mapping diagram, obtaining a two-dimensional cyclic spectrum coherence mapping diagram, and performing discrete estimation on the two-dimensional cyclic spectrum coherence mapping diagram to obtain discrete cyclic spectrum coherence mapping diagrams at different times.

[0156] The decomposition module is used for constructing a snapshot matrix according to the discrete cyclic spectrum coherence mapping diagrams at different times, and performing eigen-orthogonal decomposition on the snapshot matrix to determine the proportion of the energy of each order mode.

[0157] The reconstruction module is configured to determine the number of reconstruction modes according to the proportion change of the energy of each order mode, reconstruct the eigenvalue orthogonal decomposition result of the snapshot matrix according to the number of reconstruction modes, and obtain a cyclic spectrum coherence mapping of the fault mode, and integrate the cyclic spectrum coherence mapping of the fault mode in the spectral frequency direction to obtain a demodulation spectrum corresponding to the fault mode.

[0158] The analysis module is configured to analyze the correlation between the demodulation spectrum corresponding to the fault mode and the two-dimensional cyclic spectrum coherence mapping, and determine the occurrence time of the fault.

[0159] For example, the discrete estimation module is specifically configured to:

[0160] According to the following formula, the cyclic spectrum correlation mapping CSC of the preprocessed equal-length support bearing vibration signal is obtained by using two-dimensional Fourier transform x (alpha, f):

[0161] CSC x (alpha, f) = integral integral R x (t, tau) e -2πj(αt+fτ) dtdtau = integral integral E{x(t)x * (t+tau)} e -2πj(αt+fτ) dtdtau;

[0162] Wherein, x(t) represents the support bearing vibration signal at t time, * represents the conjugate operation, x * (t+tau) represents the conjugate operation on the support bearing vibration signal x(t+tau) at t+tau time, R x (t, v) represents the two-dimensional time-varying autocorrelation function of the support bearing vibration signal, tau represents the time delay, E represents the mathematical expectation, e represents the natural constant, j represents the imaginary unit; Alpha is the cyclic frequency, which corresponds to the frequency domain mapping of time t; F is the spectral frequency, which corresponds to the frequency domain mapping of time delay tau;

[0163] According to the following formula, the cyclic spectrum correlation mapping CSC x (alpha, f) is amplitude normalized to obtain a two-dimensional cyclic spectrum coherence mapping CSCoh x (alpha, f):

[0164]

[0165] Wherein, CSC x (0, f) represents the cyclic energy distribution of the support bearing vibration signal in the (0, f) plane, CSC x (0, f+alpha) represents the cyclic energy distribution of the support bearing vibration signal in the (0, f+alpha) plane.

[0166] According to the following formula, the two-dimensional cyclic spectrum coherence mapping CSCoh xDiscrete estimation of (α,f) yields:

[0167] S M×N (m,n)=CSCoh x (α m ,f n );

[0168] Among them, CSCoh x (α m ,f n ) represents the frequency α in the m-th cycle. m and the nth spectral frequency f n The cyclic spectral coherence value at point S M×N (m,n) represents a discrete cyclic spectral coherence mapping S of size M×N. M×N The elements in the m-th row and n-th column of the array are m = 1, 2, ..., M, n = 1, 2, ..., N, where M represents the number of cyclic frequency points and N represents the number of spectral frequency points.

[0169] For example, the decomposition module is specifically used for:

[0170] Based on the discrete cyclic spectrum coherence mapping at different times Construct a snapshot matrix X of size (M×N)×K, and represent the snapshot matrix X as follows:

[0171]

[0172] in, S represents the discrete cyclic spectrum coherence mapping at time 1, time 2, ..., time K, respectively. 1 (1,1), ..., S 1 (1, N) represent the elements in the first row and first column, ..., Nth column of the discrete cyclic spectral coherence map at time 1, respectively. 1 (2,1), ...,S 1 (2, N) represent the elements in the 1st, ..., Nth columns of the 2nd row of the discrete cyclic spectral coherence map at time 1, respectively, and S 1 (M,1), ...,S 1 (M, N) represent the elements in the M-th row and 1st column, ..., n-th column of the discrete cyclic spectral coherence map at time 1, respectively. 2 (1,1), ..., S 2 (1, N) represent the elements in the first row and first column, ..., Nth column of the discrete cyclic spectral coherence map at time 2, respectively, and S 2 (2,1), ...,S 2 (2, N) represent the elements in the 1st, ..., Nth columns of the 2nd row of the discrete cyclic spectral coherence map at the 2nd time, respectively, and S 2(M,1),..., S 2 (M,N) represent the first column,..., the Nth column element in the Mth row of the discrete cyclic spectrum coherence mapping diagram at the second moment, respectively, S K (1,1),..., S K (1,N) represent the first column,..., the Nth column element in the first row of the discrete cyclic spectrum coherence mapping diagram at the Kth moment, respectively, S K (2,1),..., S K (2,N) represent the first column,..., the Nth column element in the second row of the discrete cyclic spectrum coherence mapping diagram at the Kth moment, respectively, S K (M,1),..., S K (M,N) represent the first column,..., the Nth column element in the Mth row of the discrete cyclic spectrum coherence mapping diagram at the Kth moment, respectively;

[0173] According to the following formula, the intrinsic orthogonal decomposition of the snapshot matrix X is carried out based on singular value decomposition, to obtain the decomposition modal matrix U and the modal energy value Λ:

[0174] X=U×Λ×V T ;

[0175] Wherein, the size of the decomposition modal matrix U is (M*N) * K and u i represents the i th modal obtained by intrinsic orthogonal decomposition, and represents a real number column vector space with a dimension of M*N, i=1, 2,..., K; Λ=diag(σ1,σ2,...,σ K ) is a diagonal matrix, representing the modal energy; σ1,σ2,...,σ K represent the first order modal energy, the second order modal energy,..., the Kth order modal energy, respectively, and σ1>σ2>...>σ K ; V represents the time coefficient, and T represents the matrix transpose;

[0176] According to the following formula, the specific proportion of each order modal energy is calculated:

[0177]

[0178] Wherein, represents the proportion of the i th order modal energy σ i .

[0179] Exemplarily, the reconstruction module is specifically used for:

[0180] Calculating the modal energy difference vector dσ, and the Jth element dσ J in the modal energy difference vector dσ is represented as:

[0181] dσ J =σJ+1 -σ J ,J=1,2,…,K-1;

[0182] wherein, σ J represents the proportion of the Jth modal energy, σ J+1 represents the proportion of the (J+1)th modal energy.

[0183] According to the difference value corresponding to each element in the modal energy difference vector dσ, the position of the largest modal energy change is determined, and the reconstruction modal number Q is determined according to the following formula:

[0184] Q = argmax{dσ};

[0185] wherein, argmax{dσ} represents the index value corresponding to the maximum value of the elements in the modal energy difference vector dσ.

[0186] According to the reconstruction modal number Q, the decomposition modal matrix U is reconstructed to obtain the reconstruction vector The reconstruction vector is converted into a two-dimensional matrix Y M×N as a cyclic spectrum coherence mapping diagram of the fault modal, wherein:

[0187]

[0188]

[0189] wherein, respectively represent the 1st, M×1, M+1, M×2, M×(N-1)+1, M×N elements in the reconstruction vector .

[0190] The cyclic spectrum coherence mapping diagram Y M×N of the fault modal is integrated along the spectral frequency direction to obtain the demodulation spectrum corresponding to the fault modal,

[0191]

[0192] wherein, FDS(α m ) represents the fault demodulation spectrum amplitude at the mth cyclic frequency α m , and Y(m,n) represents the mth row and nth column element in the cyclic spectrum coherence mapping diagram Y M×N of the fault modal.

[0193] Exemplarily, the reconstruction module is further used for:

[0194] According to the following formula, the demodulation spectrum corresponding to the fault modal is subjected to 0-1 normalization processing:

[0195]

[0196] wherein, NFDS (a m ) represents the normalized fault demodulation spectrum, min{FDS (a m )} represents the minimum value of FDS (a m ), and max{FDS (a m )} represents the maximum value of FDS (a m ).

[0197] The specific implementation method of the rotor rub-impact identification system based on cyclic spectrum eigen-orthogonal decomposition provided by the embodiment of the application can be referred to the rotor rub-impact identification method based on cyclic spectrum eigen-orthogonal decomposition provided by the embodiment of the application, which will not be described here again.

[0198] The rotor rub-impact identification system based on cyclic spectrum eigen-orthogonal decomposition provided by the embodiment of the application, compared with the prior art, adopts cyclic spectrum analysis of vibration signals collected at different times, and adopts eigen-orthogonal decomposition to decompose the cyclic spectrum mapping diagrams at different times, reconstructs the fault mode according to the distribution of the modal energy change, identifies the fault type according to the reconstructed mode, and further analyzes the correlation between the fault mode and the cyclic spectrum at different times, so as to identify the fault occurrence time in the continuous monitoring of the gas turbine and the combined cycle unit. The rotor rub-impact identification system based on cyclic spectrum eigen-orthogonal decomposition provided by the embodiment of the application can not only accurately extract the fault mode for the fault diagnosis of the unit, but also can identify the occurrence time of the unit fault in the continuous monitoring, so as to meet the industrial field demand of the rotor vibration monitoring of the gas turbine and the combined cycle unit in the continuous monitoring scene.

[0199] Those skilled in the art can understand that the above-mentioned embodiments are specific embodiments for implementing the present application, and in actual application, various changes can be made in form and details without departing from the spirit and scope of the present application.

Claims

1. A rotor rub identification method based on cyclic spectral eigenorthogonal decomposition, characterized in that, The method comprises the following steps: Step S1: continuously collecting vibration signals of support bearings of a continuously running gas turbine and a combined cycle unit at equal intervals and equal lengths, and preprocessing the vibration signals of the support bearings; Step S2: obtaining a cyclic spectrum correlation mapping of the preprocessed vibration signals of the support bearings of equal lengths, performing amplitude normalization on the cyclic spectrum correlation mapping to obtain a two-dimensional cyclic spectrum coherence mapping, and performing discrete estimation on the two-dimensional cyclic spectrum coherence mapping to obtain discrete cyclic spectrum coherence mappings at different time instants; Step S3: constructing a snapshot matrix according to the discrete cyclic spectrum coherence mappings at different time instants, and performing eigen-orthogonal decomposition on the snapshot matrix to determine proportions of energy of modes of different orders; Step S4: determining a number of reconstructed modes according to changes in the proportions of energy of the modes of different orders, reconstructing the eigen-orthogonal decomposition result of the snapshot matrix according to the number of reconstructed modes, obtaining a cyclic spectrum coherence mapping of a fault mode, and performing integration on the cyclic spectrum coherence mapping of the fault mode in a spectral frequency direction to obtain a demodulation spectrum corresponding to the fault mode; Step S5: analyzing a correlation between the demodulation spectrum corresponding to the fault mode and the two-dimensional cyclic spectrum coherence mapping to determine a time of occurrence of a fault.

2. The cyclic spectral-based eigenvalue orthogonal decomposition-based rubbing identification method according to claim 1, characterized in that, The step S2 comprises the following steps: Step S21: According to the following formula, a cyclic spectrum correlation map CSC of the pre-processed equal-length support bearing vibration signal is obtained by using two-dimensional Fourier transform x (a,f): CSC x (α,f) = ∫∫R x (t,τ)e -2πj(αt+fτ) dtdτ = ∫∫E{x(t)x * (t+τ)}e -2πj(αt+fτ) dtdτ; wherein x(t) represents the support bearing vibration signal at time t, * represents a conjugate operation, x * (t+τ) represents a conjugate operation on the support bearing vibration signal x(t+τ) at time t+τ, R x (t,τ) represents a two-dimensional time-varying autocorrelation function of the support bearing vibration signal, τ represents a time delay, E represents a mathematical expectation, e represents a natural constant, j represents an imaginary unit; α is a cyclic frequency, corresponding to a frequency domain mapping of time t; f is a spectral frequency, corresponding to a frequency domain mapping of time delay τ; Step S22: Apply the following formula to the cyclic spectrum correlation mapping diagram CSC. x Amplitude normalization is performed on (α,f) to obtain the two-dimensional cyclic spectrum coherence mapping CSCOh. x (α,f): wherein CSC x (0,f) represents the cyclic energy distribution of the support bearing vibration signal in the (0,f) plane, CSC x (0,f+α) represents the cyclic energy distribution of the support bearing vibration signal in the (0,f+α) plane; Step S23: Discrete estimate of the two-dimensional cyclic spectrum coherence map CSCoh(α, f) from x (α,f) is obtained as S M×N (m,n) = CSCoh x (α m ,f n ) ; wherein, CSCoh x (α m ,f n ) represents a cyclic spectral coherence value at the mth cyclic frequency α m and the nth spectral frequency f n , S M×N (m,n) represents the mth row, mth column element in the discrete cyclic spectral coherence map S M×N of size M x N, m = 1, 2, …, M, n = 1, 2, …, N, M represents the number of cyclic frequency points, and N represents the number of spectral frequency points.

3. The cyclic spectral-based eigenvalue orthogonal decomposition-based rubbing identification method according to claim 2, characterized in that, The step S3 comprises the following steps: Step S31: According to the discrete cyclic spectrum coherence mapping at different times A snapshot matrix X of size (M x N) x K is constructed, and the snapshot matrix X is represented as: wherein, S (1,1), S (1,2),..., S (1,N) respectively represent the first column, the second column,..., the Nth column element in the first row of the discrete cyclic spectrum coherence mapping diagram at the first time point, 1 1 S (2,1), S (2,2),..., S (2,N) respectively represent the first column, the second column,..., the Nth column element in the second row of the discrete cyclic spectrum coherence mapping diagram at the first time point, 1 1 S (2,1), S (2,2),..., S (2,N) respectively represent the first column, the second column,..., the Nth column element in the second row of the discrete cyclic spectrum coherence mapping diagram at the first time point, 1 1 S (M,1), S (M,2),..., S (M,N) respectively represent the first column, the second column,..., the Nth column element in the Mth row of the discrete cyclic spectrum coherence mapping diagram at the first time point, 2 2 S (1,1), S (1,2),..., S (1,N) respectively represent the first column, the second column,..., the Nth column element in the first row of the discrete cyclic spectrum coherence mapping diagram at the second time point, 2 2 S (2,1), S (2,2),..., S (2,N) respectively represent the first column, the second column,..., the Nth column element in the second row of the discrete cyclic spectrum coherence mapping diagram at the second time point, 2 2 S (M,1), S (M,2),..., S (M,N) respectively represent the first column, the second column,..., the Nth column element in the Mth row of the discrete cyclic spectrum coherence mapping diagram at the second time point, K K S (1,1), S (1,2),..., S (1,N) respectively represent the first column, the second column,..., the Nth column element in the first row of the discrete cyclic spectrum coherence mapping diagram at the Kth time point, K K S (2,1), S (2,2),..., S (2,N) respectively represent the first column, the second column,..., the Nth column element in the second row of the discrete cyclic spectrum coherence mapping diagram at the Kth time point, K K S (M,1), S (M,2),..., S (M,N) respectively represent the first column, the second column,..., the Nth column element in the Mth row of the discrete cyclic spectrum coherence mapping diagram at the Kth time point,​​​​​​​​​ Step S32: performing eigen-orthogonal decomposition on the snapshot matrix X based on singular value decomposition according to the following formula to obtain a decomposition mode matrix U and a mode energy value Λ: X = U x Λ x V T ; wherein the size of the decomposition modal matrix U is (M×N)×K and u i represents the i-th modal obtained by the eigen-orthogonal decomposition and represents a real number column vector space with a dimension of M×N, i=1, 2, …, K; Λ=diag(σ1, σ2, …, σ K ) is a diagonal matrix, representing modal energy; σ1, σ2, …, σ K respectively represent the 1st-order modal energy, the 2nd-order modal energy, …, the Kth-order modal energy, and σ1>σ2>…>σ K ; V represents a time coefficient, and T represents matrix transposition; Step S33: calculating specific proportions of energy of modes of different orders according to the following formula: wherein, represents the proportion of the i-th order modal energy σ i of the total energy.

4. The cyclic spectral-based eigenvalue orthogonal decomposition-based rubbing identification method according to claim 3, characterized in that, The step S4 comprises the following steps: Step S41: Calculate the modal energy difference vector dσ, and the Jth element dσ J is expressed as: dσ J = σ J+1 - σ J J = 1, 2, …, K - 1; wherein σ J represents the proportion of the Jth modal energy, σ J+1 represents the proportion of the (J+1)th modal energy; Step S42: determining a position of a largest change in mode energy according to a difference value corresponding to each element in a mode energy difference vector dσ, and determining a number of reconstructed modes Q according to the following formula: Q=argmax{dσ}; wherein, argmax{dσ} represents an index value corresponding to a maximum value of elements in the mode energy difference vector dσ; Step S43: reconstructing the decomposition modal matrix U according to the reconstruction modal number Q to obtain a reconstruction vector The reconstruction vector is converted into a two-dimensional matrix Y M×N The cyclic spectrum coherence mapping as a fault mode, wherein: wherein, respectively represent the 1st, Mx1, M+1, Mx2, Mx(N-1)+1, MxN elements in the reconstructed vector ​ Step S44: According to the following formula, the cyclic spectrum coherence map Y of the fault mode is mapped in the spectral frequency direction M×N Full-band integration is performed to obtain the demodulation spectrum corresponding to the fault mode: where FDS(α m ) represents the fault demodulation spectrum amplitude at the mth cycle frequency α m , and Y(m, n) represents the element in the mth row and the nth column of the cyclic spectrum coherence map Y M×N of the fault mode.

5. The cyclic spectral-based eigenvalue decomposition-based rubbing identification method according to claim 4, characterized in that, The step S4 further comprises: performing 0-1 normalization processing on the demodulation spectrum corresponding to the fault mode according to the following formula: where NFDS(a m ) represents the normalized fault demodulation spectrum, min{FDS(a m )} represents the minimum value of FDS(a m ), and max{FDS(a m )} represents the maximum value of FDS(a m ).

6. A rotor rub identification system based on cyclic spectral eigen-decomposition, characterized in that, It comprises: a preprocessing module configured to continuously collect vibration signals of support bearings of a continuously running gas turbine and a combined cycle unit at equal intervals and equal lengths, and to preprocess the vibration signals of the support bearings; a discrete estimation module configured to obtain a cyclic spectrum correlation mapping of the preprocessed vibration signals of the support bearings of equal lengths, to perform amplitude normalization on the cyclic spectrum correlation mapping to obtain a two-dimensional cyclic spectrum coherence mapping, and to perform discrete estimation on the two-dimensional cyclic spectrum coherence mapping to obtain discrete cyclic spectrum coherence mappings at different time instants; a decomposition module configured to construct a snapshot matrix according to the discrete cyclic spectrum coherence mappings at different time instants, and to perform eigen-orthogonal decomposition on the snapshot matrix to determine proportions of energy of modes of different orders; The reconstruction module is configured to determine a reconstruction modal number according to a change of a proportion of each order modal energy, reconstruct an eigenvalue orthogonal decomposition result of the snapshot matrix according to the reconstruction modal number, and obtain a cyclic spectral coherence mapping of a fault modal, and integrate the cyclic spectral coherence mapping of the fault modal in a full-band direction along a spectral frequency direction to obtain a demodulation spectrum corresponding to the fault modal. The analysis module is configured to analyze a correlation between the demodulation spectrum corresponding to the fault modal and the two-dimensional cyclic spectral coherence mapping, and determine a fault occurrence time.

7. The cyclic spectral-based eigenvalue orthogonal decomposition based rubbing identification system of claim 6, wherein, The decomposition module is specifically configured to: The cyclic spectrum correlation map CSC of the pre-processed equal-length support bearing vibration signal is obtained by using two-dimensional Fourier transform according to the following formula x (α, f): CSC x (α,f) = ∫∫R x (t,τ)e -2πj(αt+fτ) dtdτ = ∫∫E{x(t)x * (t+τ)}e -2πj(αt+fτ) dtdτ; wherein x(t) represents the support bearing vibration signal at time t, * represents a conjugate operation, x * (t+τ) represents a conjugate operation on the support bearing vibration signal x(t+τ) at time t+τ, R x (t,τ) represents a two-dimensional time-varying autocorrelation function of the support bearing vibration signal, τ represents a time delay, E represents a mathematical expectation, e represents a natural constant, j represents an imaginary unit; α is a cyclic frequency, corresponding to a frequency domain mapping of time t; f is a spectral frequency, corresponding to a frequency domain mapping of time delay τ; According to the following formula, the cyclic spectrum correlation mapping (CSC) is... x Amplitude normalization is performed on (α,f) to obtain the two-dimensional cyclic spectrum coherence mapping (CSCoh). x (α,f): wherein CSC x (0,f) represents the cyclic energy distribution of the support bearing vibration signal in the (0,f) plane, CSC x (0,f+α) represents the cyclic energy distribution of the support bearing vibration signal in the (0,f+α) plane; The two-dimensional cyclic spectral coherence map, CSCoh (a,f) is estimated by discretization according to the following formula: x (α,f) is obtained by discretization. S M×N (m,n) = CSCoh x (α m ,f n ) ; wherein, CSCoh x (α m ,f n ) represents a cyclic spectral coherence value at the mth cyclic frequency α m and the nth spectral frequency f n , S M×n (m,n) represents an element in the mth row and the nth column of a discrete cyclic spectral coherence map S M×N of size M x N, m = 1, 2, …, M, n = 1, 2, …, N, M represents the number of cyclic frequency points, and N represents the number of spectral frequency points.

8. The cyclic spectral-based eigenvalue orthogonal decomposition based rubbing identification system of claim 7, wherein, The decomposition module is specifically configured to: Discrete cyclic spectral coherence maps at different times A snapshot matrix X of size (M x N) x K is constructed, and the snapshot matrix X is represented as: wherein, S (1), S (2), …, S (K) respectively represent the discrete cyclic spectrum coherence mapping diagram at the 1st time, the 2nd time, …, the Kth time, 1 (1,1), …, S 1 (1,N) respectively represent the 1st column, …, the Nth column element in the 1st row of the discrete cyclic spectrum coherence mapping diagram at the 1st time, S 1 (2,1), …, S 1 (2,N) respectively represent the 1st column, …, the Nth column element in the 2nd row of the discrete cyclic spectrum coherence mapping diagram at the 1st time, S 1 (M,1), …, S 1 (M,N) respectively represent the 1st column, …, the Nth column element in the Mth row of the discrete cyclic spectrum coherence mapping diagram at the 1st time, S 2 (1,1), …, S 2 (1,N) respectively represent the 1st column, …, the Nth column element in the 1st row of the discrete cyclic spectrum coherence mapping diagram at the 2nd time, S 2 (2,1), …, S 2 (2,N) respectively represent the 1st column, …, the Nth column element in the 2nd row of the discrete cyclic spectrum coherence mapping diagram at the 2nd time, S 2 (M,1), …, S 2 (M,N) respectively represent the 1st column, …, the Nth column element in the Mth row of the discrete cyclic spectrum coherence mapping diagram at the 2nd time, S K (1,1), …, S K (1,N) respectively represent the 1st column, …, the Nth column element in the 1st row of the discrete cyclic spectrum coherence mapping diagram at the Kth time, S K (2,1), …, S K (2,N) respectively represent the 1st column, …, the Nth column element in the 2nd row of the discrete cyclic spectrum coherence mapping diagram at the Kth time, S K (M,1), …, S K (M,N) respectively represent the 1st column, …, the Nth column element in the Mth row of the discrete cyclic spectrum coherence mapping diagram at the Kth time. The snapshot matrix X is subjected to eigenvalue orthogonal decomposition based on singular value decomposition according to the following formula to obtain a decomposition modal matrix U and a modal energy value Λ: X = U x Λ x V T ; wherein the size of the decomposition modal matrix U is (M×N)×K and u i represents the i-th modal obtained by the eigen-orthogonal decomposition and represents a real number column vector space with a dimension of M×N, i=1, 2, …, K; Λ=diag(σ1, σ2, …, σ K ) is a diagonal matrix, representing modal energy; σ1, σ2, …, σ K respectively represent the 1st-order modal energy, the 2nd-order modal energy, …, the Kth-order modal energy, and σ1>σ2>…>σ K ; V represents a time coefficient, and T represents a matrix transpose; The specific proportion of each order modal energy is calculated according to the following formula: wherein, represents the i-th order modal energy σ i of the proportion.

9. The cyclic spectral-based eigenvalue orthogonal decomposition based rubbing identification system of claim 8, wherein, The reconstruction module is specifically configured to: a modal energy difference vector dσ, the jth element dσjof the modal energy difference vector dσ J is represented as: dσ J = σ J+1 - σ J J = 1, 2, …, K - 1; wherein σ J represents the proportion of the Jth modal energy, σ J+1 represents the proportion of the (J+1)th modal energy; The reconstruction module is specifically configured to: The position of the largest modal energy change is determined according to a differential value corresponding to each element in the modal energy difference vector dσ, and the reconstruction modal number Q is determined according to the following formula: Q = argmax {dσ}; According to the number of reconstruction modalities Q, the decomposition modality matrix U is reconstructed to obtain a reconstruction vector The reconstruction vector is converted into a two-dimensional matrix Y M×N The cyclic spectrum coherence mapping diagram as the fault modality, wherein: wherein, respectively represent the 1st, Mx1, M+1, Mx2, Mx(N-1)+1, MxN elements in the reconstructed vector ​ Cyclic spectrum coherence mapping of the fault mode in the spectral frequency direction Y M×N Full-band integration is performed to obtain the demodulation spectrum corresponding to the fault mode, where FDS(α m ) represents the fault demodulation spectrum amplitude at the mth cycle frequency α m , and Y(m, n) represents the element in the mth row and the nth column of the cyclic spectrum coherence map Y M×N of the fault mode.

10. The cyclic spectral-based eigenvalue orthogonal decomposition based rubbing identification system of claim 9, wherein, Wherein, argmax {dσ} represents an index value corresponding to a maximum value of the elements in the modal energy difference vector dσ. The reconstruction module is further configured to: The demodulation spectrum corresponding to the fault modal is subjected to 0-1 normalization processing according to the following formula: where NFDS(a m ) represents the normalized fault demodulation spectrum, min{FDS(a m )} represents the minimum value of FDS(a m ), and max{FDS(a m )} represents the maximum value of FDS(a m ).

Images