Rotor damping ratio identification method based on key phase sensing measurement

Through the rotor damping ratio identification method based on bond sensing measurement, the Hampel filtering algorithm and vectorized calculation optimized signal processing is used to solve the real-time and accuracy problems of modal identification in rotary power machinery, and the efficient and stable monitoring and fault warning of the rotor system are achieved.

CN120408559AActive Publication Date: 2025-08-01BEIJING UNIV OF CHEM TECH +1
View PDF 4 Cites 0 Cited by

Patent Information

Application Number
CN202510746635.4
Authority / Receiving Office
CN · China
Patent Type
Applications(China)
Current Assignee / Owner
Filing Date
2025-06-05
Publication Date
2025-08-01
Estimated Expiration
2045-06-05

AI Technical Summary

Technical Problem

There are data processing methods in existing monitoring systems that are difficult to achieve real-time analysis and early warning response in rotary power machinery. Multi-source signal feature extraction is susceptible to background noise interference, resulting in the accumulation of modal parameter identification errors, affecting the accuracy of fault warning.

Method used

The rotor damping ratio identification method based on bond-phase sensing measurement is adopted. Time domain signals are collected in real time through bond-phase sensors, harmonic interference is suppressed using improved Hampel filtering algorithm, and signal processing flow is optimized through vectorized calculations, rotor damping ratio identification model is constructed, harmonic interference components are eliminated, and identification efficiency and accuracy are improved.

Benefits of technology

It realizes efficient, stable and accurate modal identification of the rotor system, improves the timeliness and accuracy of real-time monitoring and fault warning, and is suitable for online monitoring and fault diagnosis.

✦ Generated by Eureka AI based on patent content.

Smart Images

  • Figure CN120408559A_ABST
    Figure CN120408559A_ABST
Patent Text Reader

Abstract

The invention discloses a rotor damping ratio identification method based on key phase sensing measurement, and belongs to the related technical field of rotating machinery. The method comprises the following steps: carrying out band-pass filtering on a key phase signal; after the auto-power spectral density is calculated, vectorized Hampel filtering is implemented, and a window matrix is constructed to detect and replace an abnormal value formed by harmonic waves; reconstructing a complex frequency spectrum and carrying out inverse transformation to generate a time domain pulse attenuation signal; a reverse autoregression model is adopted to extract a feature root, and stability parameters are automatically obtained in combination with a clustering algorithm. According to the method, the modal parameters of the rotor are analyzed through the key-phase signals, and the application of the key-phase signals in dynamic analysis of the rotor is expanded; in the signal analysis process, a power spectral density function is regarded as a special time sequence, an excellent harmonic interference effect can be realized by adopting vectorization Hampel filtering, and a circulating sliding window is replaced by vectorization operation, so that the calculation efficiency is remarkably improved, and a reliable technical means is provided for real-time monitoring of the rotating machinery.
Need to check novelty before this filing date? Find Prior Art

Description

Technical Field

[0001] The design of the present invention belongs to the technical field related to rotating machinery, and specifically relates to a method for identifying the rotor damping ratio based on key phase sensing measurement. Background Art

[0002] As the core power device in industrial fields such as petrochemical, power, and metallurgy, the operating stability of rotating power machinery directly affects the reliability of the industrial production system. With the evolution of modern industrial equipment towards high power density, supercritical parameters, and extreme flow rates, the accurate identification of the dynamic characteristics of rotor systems faces severe challenges. Engineering practice shows that constructing an early warning mechanism by real-time monitoring of key parameters such as vibration spectrum, pressure pulsation, and temperature field distribution is an effective means to achieve equipment status prediction. However, there are deficiencies in existing monitoring systems. On the one hand, data processing methods are difficult to achieve real-time analysis and early warning response of massive monitoring data, resulting in a significant lag in fault early warning. On the other hand, during the extraction of multi-source signal features, it is easily interfered by background noise, causing the cumulative amplification of modal parameter identification errors, which seriously affects the accuracy of early warning. Among them, modal analysis, as a key technology, can effectively identify the dynamic characteristics of structures and has extremely important engineering value in many aspects such as the optimal design of structures, fault diagnosis, damage identification, and health monitoring. However, traditional OMA modal analysis mostly performs modal identification under high signal-to-noise ratios in the laboratory. For equipment such as engines, turbines, and compressors, during the actual engineering operation process, the system structure is often excited by non-white noise, such as harmonic excitation, which poses challenges to traditional operational modal parameter identification methods.

[0003] Therefore, optimizing the data processing flow to reduce the data volume, solving the harmonic interference problem, and improving the identification efficiency and accuracy will contribute to achieving efficient, stable, and accurate modal identification of rotating power machinery, which has great engineering significance for the stable operation of equipment. Summary of the Invention

[0004] The technical objective of the present invention is to propose a method for identifying the rotor damping ratio based on keyphasor sensing measurement. This method combines the industrial data processing in the condition monitoring of rotating machinery with the dynamic mechanism, identifies the modal parameters of the rotor system through the time-domain signals collected in real time by the keyphasor sensor, uses an improved Hampel filtering algorithm to suppress harmonic interference, and optimizes the signal processing flow based on the vectorized calculation method. Among them, a rotor damping ratio identification model based on the analysis of keyphasor sensor signals is constructed, which expands the application scope of keyphasor signals in the identification of rotor dynamic parameters; the Hampel filter is introduced to suppress non-Gaussian noise in the original signal, and the harmonic interference components are effectively removed by dynamically adjusting the sliding window threshold, solving the problem of modal aliasing caused by periodic process disturbances in traditional industrial data processing; the matrix vectorized calculation method is used to reconstruct the filter algorithm in parallel, and the sliding window filtering algorithm is reconstructed into a two-dimensional matrix block operation, optimizing the industrial data processing flow and improving the overall efficiency. This method effectively solves the problems of insufficient real-time performance and limited accuracy of modal parameter identification existing in the existing monitoring technology for rotor damping ratio identification.

[0005] To achieve the above objective, the technical solution adopted by the present invention is a method for identifying the rotor damping ratio based on keyphasor sensing measurement, including:

[0006] S1. Obtain the keyphasor signal of the keyphasor sensor on the vibration measurement plane of the rotating machinery, perform a fast Fourier transform (FFT) on the keyphasor signal, calculate the modulus square of the complex result in the frequency domain, and obtain the auto-power spectral density function of the keyphasor signal characterizing the vibration characteristics of the rotor system through single-sided spectrum normalization.

[0007] S2. Perform Hampel filtering on the auto-power spectral density function of the keyphasor signal. First, construct a dynamic sliding window on the auto-power spectral density function of the keyphasor signal according to the frequency band distribution characteristics of harmonic interference in industrial data; the width of the dynamic sliding window is dynamically adjusted according to the rotor speed, and the width of the dynamic sliding window is inversely proportional to the speed; secondly, use the vectorized processing method to convert the spectral data in the dynamic sliding window into a two-dimensional matrix block, and calculate the median and median absolute deviation (MAD) of each window based on the statistical mechanism; finally, perform filtering on the power spectral density function according to the preset threshold to remove harmonic interference and obtain the Hampel-filtered power spectral density function.

[0008] S3. Perform an inverse Fourier transform on the Hampel-filtered power spectral density function, and normalize the time-domain decay signal after the inverse Fourier transform to ensure the numerical stability of subsequent parameter identification, and obtain the Hampel-filtered time-domain pulse decay signal.

[0009] S4. Construct an autoregressive model in reverse based on the time-domain pulse decay signal after Hampel filtering, and calculate the corresponding stability parameters in each order of the rotor system, including damping ratio and natural frequency; according to the statistical distribution characteristics of industrial data, use the 3σ criterion to eliminate abnormal data from the stability parameters, and use the mean value of normal operating condition data as the identification result of the dynamic characteristics of the rotor system.

[0010] Further, the S1 includes the following steps:

[0011] S1.1. Synchronously collect the key-phase signal of the measurement plane where the key-phase sensor of the rotating machinery is located at a sampling frequency not lower than 2 kHz; perform preliminary filtering on the collected key-phase signal, and use a band-pass filter to remove high-frequency noise and part of the low-frequency interference. The cut-off frequency of the band-pass filter is set according to the spectral distribution characteristics of the key-phase signal.

[0012] S1.2. Perform a fast Fourier transform (FFT) on the key-phase signal, calculate the modulus square of the frequency-domain complex result, and obtain the auto-power spectral density function of the key-phase signal characterizing the vibration characteristics of the rotor system through single-sided spectrum normalization.

[0013] Further, the S2 includes the following steps:

[0014] S2.1. Supplement a specific number of data points (usually zero-value filling or mirror-symmetric extension) at the start and end positions of the auto-power spectral density function of the key-phase signal. Form a dynamic sliding window centered on each original data point (non-supplemented data point), and each dynamic sliding window is represented as each row of the matrix to form a window matrix of the auto-power spectral density function of the key-phase signal.

[0015] S2.2. Obtain the median of each dynamic sliding window by sorting the power spectral density window matrix of the key-phase signal and taking the middle value, and calculate the absolute difference between each element in the window and the window median. Take the median of these absolute differences as the median absolute deviation (MAD).

[0016] S2.3. According to a preset threshold, compare the difference between each element in the sliding window and the median with the product of the absolute median deviation (MAD) and the threshold to determine outliers in the power spectral density window matrix of the key-phase signal. The determined outlier points are replaced with the median of the corresponding window, so as to batch eliminate the harmonic interference in the key-phase signal and obtain the power spectral density function without harmonics.

[0017] Further, the S4 includes the following steps:

[0018] S4.1. Perform inverse autoregressive modal parameter identification on the time-domain pulse attenuation signals respectively to obtain the characteristic roots of the discrete models of each order, and calculate the natural frequencies and damping ratios corresponding to each order in the measured keyphasor signals based on the characteristic roots.

[0019] S4.2. Use the 3σ clustering algorithm to perform automatic clustering operations on the natural frequencies and damping ratios of different orders in the keyphasor signals. Through clustering analysis, similar dynamic characteristic parameters are grouped into one category, thereby realizing the efficient identification and classification of the dynamic characteristics of the rotor system. Take the mean value of the normal condition data as the identification result of the dynamic characteristics of the rotor system to ensure the reliability and accuracy of the result.

[0020] The present invention has the following advantages:

[0021] The rotor stability parameter identification method based on anti-harmonic interference of keyphasor signals proposed by the present invention has remarkable innovation and practicability. This technology expands the use of keyphasor signals in rotor dynamics analysis. Compared with the existing technology where keyphasor signals are only used for speed extraction and phase-based synchronous sampling, this method constructs a power spectrum function of keyphasor signals and then analyzes the modal parameters of the rotor. In the signal processing link, this method introduces the Hampel filtering algorithm. The signal processed by the Hampel filter greatly restores the modal characteristics of the rotor system in the power spectrum, which simplifies the subsequent modal parameter extraction process. Compared with traditional signal processing methods, the computational complexity of the Hampel filtering algorithm is at a relatively low level. This advantage makes it particularly suitable for real-time or near-real-time data processing requirements. In key fields such as online monitoring and fault diagnosis in industrial production, the timeliness of data processing is crucial. Moreover, the present invention further vectorizes the Hampel filtering method. Through vectorization, the calculation speed of this method is significantly improved. In the online monitoring scenario, a fast calculation speed means that subtle changes in the operating state of the rotor system can be captured more timely, thereby providing more efficient and accurate operating data for the system. This has inestimable value for realizing real-time monitoring and fault warning of the rotor system. Description of the Drawings

[0022] The following further elaborates on the present invention in detail in conjunction with the drawings and specific implementation methods. The drawings described herein are used to provide a further understanding of the present invention, form a part of this application, and do not constitute an improper limitation to the present invention.

[0023] Figure 1 is the implementation flowchart of the embodiment of this application;

[0024] Figure 2 is the amplitude-frequency curve diagram of the analytical signal power spectral density function in the embodiment of this application;

[0025] Figure 3 Schematic diagram of the vectorized Hampel filtering method in the embodiment of the present application;

[0026] Figure 4 Amplitude-frequency curve diagram after removing harmonic interference from the power spectral density function in the embodiment of the present application;

[0027] Figure 5 Time-domain pulse decay signal curve diagram before filtering in the embodiment of the present application;

[0028] Figure 6 Time-domain pulse decay signal curve diagram after filtering in the embodiment of the present application;

[0029] Figure 7 Schematic diagram of the two-dimensional steady-state diagram of the fitting order and natural frequency in the embodiment of the present application;

[0030] Figure 8 Schematic diagram of the two-dimensional steady-state diagram of the fitting order and damping ratio in the embodiment of the present application;

[0031] Figure 9 Automatic clustering result diagram of the steady-state diagram data corresponding to the attenuation signal in the embodiment of the present application.[[ID=2 June 26]] Detailed implementation manners

[0032] The technical solutions of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. The description of the exemplary embodiments is merely illustrative and in no way limits the present disclosure and its application or use. The present disclosure can be implemented in many different forms and is not limited to the embodiments described herein. These embodiments are provided to make the present disclosure thorough and complete, and to fully convey the scope of the present disclosure to those skilled in the art.

[0033] The embodiments of the present application will be described in detail below. Please refer to Figure 1 , the present invention discloses a rotor damping ratio identification method based on key-phase sensing measurement, including:

[0034] S1. Obtain the key-phase signal of the key-phase sensor on the vibration measurement plane of the rotating machinery, and synchronously collect the signal at a sampling frequency not lower than 2 kHz; perform a fast Fourier transform (FFT) on the key-phase signal, calculate the modulus square of the complex result in the frequency domain, and obtain the auto-power spectral density function characterizing the vibration characteristics of the rotor system through single-sided spectrum normalization.

[0035] Obtaining the key-phase signal of the rotor: A small raised metal sheet or a groove is attached to the shaft surface as the key-phase, and an eddy current displacement sensor is used to measure the pulse signal at the key-phase position.

[0036] S1.1 Synchronously collect the keyphasor signal of the measurement plane where the keyphasor sensor of the rotating machinery is located at a sampling frequency not lower than 2 kHz; perform preliminary filtering on the collected keyphasor signal, and use a band-pass filter to remove high-frequency noise and some low-frequency interference. The cut-off frequencies of the filter are set according to the spectral distribution characteristics of the keyphasor signal.

[0037] Input the collected keyphasor signal z(t) into the band-pass filter, and its passband range is determined according to the frequency characteristics of the keyphasor signal. The passband range is [0.8 f0, 1.2(M + 1) f0] (f0 is the real-time rotational frequency, M is the order of the harmonic to be concerned about, and the default value of M is 5), to filter out asynchronous vibration interference;

[0038] Example: If the rotational frequency f0 = 50 Hz is measured according to the keyphasor signal of the rotating machinery and the order of the harmonic to be concerned about M = 5, then set the passband of the band-pass filter to 40 Hz to 360 Hz to ensure that the effective signal is retained while filtering out high-frequency and low-frequency interference outside this range.

[0039] S1.2 Perform a fast Fourier transform (FFT) on the keyphasor signal, calculate the square of the modulus of the complex result in the frequency domain, and obtain the auto-power spectral density function characterizing the vibration characteristics of the rotor system through single-sided spectrum normalization.

[0040] Perform a fast Fourier transform (FFT) on the keyphasor signal z(t) to obtain Z(w), and then calculate the square of the modulus of the complex amplitude after the Fourier transform to obtain the auto-power spectral function of the keyphasor signal. Specifically, the expression of the power spectral density function is:

[0041]

[0042] Among them, G(w) represents the power spectral density at the angular frequency w, w is the angular frequency; Z(w) is the Fourier transform of the keyphasor signal z(t), that is , t is the time, j is the imaginary unit (a commonly used symbol in engineering), and e is the natural constant in mathematics (the base of the natural logarithm, approximately equal to 2.71828); is the square of the modulus of the complex amplitude after the Fourier transform.

[0043] In an embodiment, the keyphasor signal of a rotor system with a flexible foundation at a rotational speed of 2000 r / min and a sampling rate fs = 20 kHz is analyzed. First, perform a power spectral density analysis on the keyphasor signal z(t) of the rotor system with a flexible foundation, and obtain the power spectral density function G(w) diagram of the keyphasor signal as shown in Figure 2 shown. From Figure 2It can be seen that in addition to the rotation frequency (f0 = 33.3 Hz) representing the normal operation of the rotor system, there are also a large number of multiple frequencies and fractional frequencies. Near these frequency components, there are also a large number of spike signals, which are distributed at different frequency points, namely the so-called harmonic interference.

[0044] S2. Perform Hampel filtering on the self-power spectral density function of the key phase signal. First, according to the frequency band distribution characteristics of harmonic interference in industrial data, construct a dynamic sliding window on the self-power spectral density function of the key phase signal. The window width is dynamically adjusted according to the rotor speed (the window width is inversely proportional to the speed). Secondly, adopt a vectorization processing method to convert the spectral data within the sliding window into a two-dimensional matrix block, and calculate the median and absolute median difference (MAD) of each window based on statistical mechanisms. Finally, perform filtering on the power spectral density function according to a preset threshold to remove harmonic interference and obtain the self-power spectral density function after Hampel filtering.

[0045] In one embodiment, please refer to Figure 3 , the vectorized Hampel filtering method includes:

[0046] S2.1 Supplement a specific number of data points (usually zero-value filling or mirror symmetry extension) at the start and end positions of the self-power spectral density function of the key phase signal. Form a dynamic sliding window centered on each original data point (non-supplemented data point). Each dynamic sliding window is represented as each row of the matrix, and a window matrix of the self-power spectral density function of the key phase signal is formed.

[0047] Set the dynamic sliding window length L = 2k + 1, where k is the half-length of the window. L is dynamically adjusted according to the rotor speed. When L < 3, it is forced to be set to L = 3 (i.e., the window half-length k = 1) to avoid filtering failure due to too small a window.

[0048] Preset the window length L1. When L1 is an even number, the window length is L = L1 - 1; when L1 is an odd number, the window length is L = L1. Specifically, the preset window length is expressed as:

[0049]

[0050] where, f s is the sampling frequency (Hz) of the key phase signal of the vibration measurement plane key phase sensor of the rotating machinery, which is determined by the data acquisition system hardware; f r is the rotor fundamental frequency (Hz), which is obtained by calculating the pulse interval of the key phase signal; ⌊⋅⌋ represents the floor operation symbol.

[0051] Pad k (half window length) zeros before and after the power spectral density function of the key phase signal G(w) to ensure sufficient data at the boundaries during sliding. If the length of the original signal is n, the length of the extended signal is n + 2k. Specifically, the extended power spectral density function of the key phase signal is expressed as:

[0052]

[0053] where is the power spectral density function of the key phase signal after signal extension; G1, G2, …, G n are the elements in G(w), and n is the signal length of the power spectral density function of the key phase signal G(w).

[0054] Taking k adjacent points before and after the element G i in G(w) as the center, a dynamic sliding window Z of a fixed length is formed i ={G i-k , G i-k+1 , G i-k+2 ,…, G i , …, G i+k}, where i = 1, 2, 3, ……, n; A series of dynamic sliding windows are generated by traversing the extended power spectral density function of the key phase signal. Each sliding window is represented as a row of the window matrix, and finally a power spectral density window matrix with the shape of n×L is constructed, specifically expressed as:

[0055]

[0056] S2.2 For the window matrix of the power spectral density of the key phase signal, obtain the median of each dynamic sliding window by sorting and taking the middle value, and calculate the absolute difference between each element in the window and the median of the window, and take the median of these absolute differences as the median absolute deviation (MAD).

[0057] Calculate the median within each dynamic sliding window. After sorting the data within each dynamic sliding window, the (k + 1)-th position of each window is the median, specifically expressed as:

[0058]

[0059] where m is the median vector of the rows in the power signal window matrix; is the (k + 1)-th element after sorting the elements in the i-th dynamic sliding window Z i , where Z i is the i-th row of the Z G window matrix.

[0060] For each dynamic sliding window, calculate the absolute difference between each element and the median value in the window, and take the median of these absolute differences as the median absolute deviation (MAD). The specific expression is as follows:

[0061]

[0062] where Z G ij is the element in the i-th row and j-th column of the power spectral density window matrix Z G of the key phase signal; is the row MAD vector, , is the median absolute deviation of the i-th dynamic sliding window; m i is the i-th element of the median vector of the power signal window matrix.

[0063] S2.3 According to a preset threshold, compare the difference between each element in the sliding window and the median value with the product of the absolute median difference (MAD) and the threshold, and determine the outliers in the power spectral density window matrix of the key phase signal. The determined outliers are replaced with the median value of the corresponding window, so as to eliminate the harmonic interference in the batch of key phase signals and obtain the power spectral density function without harmonics.

[0064] Compare the product of the calculated median absolute value deviation and the set threshold (usually a multiple based on MAD) with each element of the power spectral density window matrix. If the absolute value of the difference between the element and the median value in the window exceeds the product of MAD and the threshold, mark the element as an outlier, and then batch replace the outliers through boolean indexing to obtain the power spectral density function G d (w). Specifically, the outlier detection is expressed as:

[0065] If the current point Z G ij satisfies: , then determine that Z G ij is an outlier.

[0066] where t0 is the outlier detection threshold, usually taken as 3.

[0067] Regarding the discrete amplitude signals in the power spectral density function as special "time series", the steep spikes formed by the harmonics in the key phase signal can just be regarded as "outliers" because they are far from the group. Using the Hampel filter to detect and remove the outliers can achieve an excellent effect of removing harmonic interference.

[0068] For example, for Figure 2 the power spectral density diagram G(w) of the key phase signal shown, the sampling rate fs = 20 kHz, and the rotor fundamental frequency f r=f0 = 33.3Hz, it is calculated that the signal window length L = 299, then the window half-length k = 149; the vectorized Hampel filter is used to process the power spectral density function of the key phase signal. When setting parameters, the threshold t0 is set to 3. After being processed by this filter, the result as shown in Figure 4 is obtained. By comparing the images before and after processing, it can be found that Figure 2 the harmonic interference that originally existed has disappeared, and the entire spectrum has become smoother, being more able to accurately reflect the true characteristics of the compressor vibration signal.

[0069] S3. Perform the inverse Fourier transform on the power spectral density function after Hampel filtering, and normalize the complex signal after the inverse transform to ensure the numerical stability of subsequent parameter identification, obtaining the time-domain pulse decay signal after Hampel filtering.

[0070] After completing the filtering process of the power spectral density function, in order to further restore the complete characteristics of the signal, it is necessary to supplement zeros in the G d (w) function to construct the complex frequency spectrum Z d (w).

[0071] The inverse Fourier transform is used to restore the complex frequency spectrum in the frequency domain to the time-domain pulse decay signal. Specifically, the time-domain pulse decay signal z1(t) is expressed as:

[0072]

[0073] After obtaining the time-domain signal z1(t), in order to ensure the stability of numerical calculation and make the dimensions of the signal consistent, the signal is processed by maximum absolute value normalization. The specific normalization is expressed as:

[0074]

[0075] For example, for the power spectral density functions shown in Figure 2 and Figure 4 , the power spectral density function G(w) without filtering and the power spectral density function G d (w) after filtering are subjected to the inverse Fourier transform, and then maximum amplitude normalization processing is performed respectively, and then the time-domain pulse decay signal corresponding to the power spectral density function before filtering shown in Figure 5 and the time-domain pulse decay signal corresponding to the power spectral density function after filtering shown in Figure 6 can be obtained. By comparing the time decay signals in Figure 5 and Figure 6 , for the time decay signal without filtering, the Figure 5 time-domain signal is covered by the rotating frequency and is approximately a cosine signal; for the Figure 6The time-domain signal effectively eliminates the rotational frequency and harmonic interference and can better represent the original attenuation characteristics.

[0076] S4. Based on the time-domain pulse attenuation signal after the Hampel filtering, construct an autoregressive model in reverse, and calculate the corresponding stability parameters (damping ratio and natural frequency) in each order of the rotor system; according to the statistical distribution characteristics of industrial data, use the 3σ criterion to eliminate abnormal data from the stability parameters, and use the mean value of the normal operating condition data as the identification result of the dynamic characteristics of the rotor system.

[0077] S4.1 Perform autoregressive modal parameter identification in reverse on the filtered time-domain pulse attenuation signal respectively to obtain the characteristic roots of the discrete models in each order; and calculate the natural frequency and damping ratio corresponding to each order in the measured key phase signal based on the characteristic roots.

[0078] Specifically, adopt a single-input single-output autoregressive discrete model for the S(t) time-domain pulse attenuation signal:

[0079] [[ID=,13]]

[0080] where S is the time-domain attenuation signal after the Hampel filtering; N is the signal sequence of the discretized time-domain attenuation signal; p is the order of the autoregressive model in reverse; b i is the system characteristic parameter; W is a white noise sequence with a mean of 0, and the sampling interval is ;

[0081] The sequence of the time-domain attenuation signal S after the Hampel filtering is 1, 2, …, N - 1. According to the above autoregressive model in reverse, construct a linear equation system:

[0082]

[0083] where S[p] abbreviates the above formula as:

[0084]

[0085] In the formula .

[0086] Perform a z-transform on the above formula to obtain:

[0087] [[ID=,42]]

[0088] In the formula , z is the discrete transformation factor; S(z) is the z-transform of the output signal S [N] ; W(z) is the z-transform of the output signal W [N] ;

[0089] In the actual signal, there is a certain amount of noise interference signal. Therefore, SVD is used to decompose the noise of A:

[0090]

[0091] where A is the original matrix of size N - p×p; U is an orthogonal matrix of size N - p×N - p (left singular value vector matrix), and U r is the direction of the retained true signal, and U r+ is the direction of high - order noise; V is the transpose of an orthogonal matrix of size p×p (right singular vector matrix), is the transpose of the right singular vector matrix retaining the first r principal components, is the transpose of the right singular vector matrix of the high - order noise part after decomposition; S is a diagonal matrix of size N - p×p (singular value matrix), is the retained part after SVD decomposition, with the order of r×r, is the high - order noise part; the minimum noise reduction order r needs to satisfy the following conditions:

[0092]

[0093] Therefore, the denoised A matrix is expressed as:

[0094]

[0095] Based on the above - mentioned linear equations, the parameter vector b parameter estimation of the model can be obtained:

[0096]

[0097] For the characteristic equation of the system , the characteristic roots corresponding to each order of the system modes can be obtained ;

[0098]

[0099] Furthermore, the natural frequency, damping ratio, and logarithmic decrement of the rotor system can be obtained

[0100]

[0101] S4.2 Use the clustering algorithm to perform an automatic clustering operation on the natural frequencies and damping ratios of different orders in the key - phase signal. By clustering analysis, similar dynamic characteristic parameters are grouped into one category, so as to realize the efficient identification and classification of the dynamic characteristics of the rotor system. Using the mean value of the normal - condition data as the identification result of the dynamic characteristics of the rotor system ensures the reliability and accuracy of the result.

[0102] After obtaining the natural frequency and damping ratio processed from the time-domain pulse attenuation signal, the automatic statistical clustering algorithm is used to further eliminate outliers for automatic monitoring and display.

[0103] Specifically, use the clustering method to cluster the identification results of the inverse autoregressive parameters of different orders

[0104]

[0105] where is the natural frequency, is the damping ratio; and are the average values of the natural frequency and damping ratio; and are the standard deviations.

[0106] After multiple iterations and eliminating outliers, calculate the mean values of the natural frequency and damping ratio after statistical clustering ( , ), and the estimated value of the system stability identification result can be obtained.

[0107] For example, perform damping ratio and natural frequency identification on the Figure 6 time-domain attenuation signal of the key phase signal, use the inverse autoregressive model to identify the stability parameters, and set the order p = 300 of the inverse autoregressive model; calculate the natural frequency and damping ratio of the rotor system with a flexible foundation, as shown in Figure 7 , 8 , where the abscissa is the order and the ordinates are the damping ratio and natural frequency respectively, and use the clustering method to obtain the final result, as shown in Figure 9 , where the center of the ellipse represents the final parameter estimation value ( ) = (49.28 Hz, 0.75%), and this value is the stability identification result of the rotor system based on the key phase signal.

[0108] The above is the specific implementation manner of the present invention, but the protection scope of the present invention is not limited thereto. Those skilled in the art of this technology can make appropriate adjustments without departing from the principle of the present invention, and these adjustments should be included in the protection scope of the present invention.

Claims

1. A method for identifying the rotor damping ratio based on keyphasor sensing measurement, characterized in that The steps include: S1. Obtaining a key phase signal of a rotating machinery vibration measuring plane key phase sensor; performing a fast Fourier transform (FFT) on the key phase signal, calculating the square of the modulus of the complex number result in the frequency domain, and obtaining an autopower spectrum density function of the key phase signal that characterizes the vibration characteristics of the rotor system through unilateral spectrum normalization; S2. Performing Hampel filtering on the power spectrum density function of the key phase signal: first, constructing a dynamic sliding window on the power spectrum density function of the key phase signal according to the frequency band distribution characteristics of the harmonic interference in the industrial data; the width of the dynamic sliding window is dynamically adjusted according to the rotor speed, wherein the width of the dynamic sliding window is inversely proportional to the speed; secondly, using a vectorized processing method, the spectrum data in the dynamic sliding window is converted into a two-dimensional matrix block, and the median and absolute median difference (MAD) of each window are calculated based on a statistical mechanism; finally, filtering the power spectrum density function according to a preset threshold to remove harmonic interference and obtain a power spectrum density function after Hampel filtering; S3, performing an inverse Fourier transform on the power spectral density function after the Hampel filter, and normalizing the time domain attenuation signal after the inverse Fourier transform to ensure the numerical stability of subsequent parameter identification, and obtaining a time domain pulse attenuation signal after the Hampel filter; S4. Construct an inverse autoregressive model based on the time-domain pulse attenuation signal after Hampel filtering, and calculate the corresponding stability parameters of the rotor system in each order, including the damping ratio and the natural frequency; according to the statistical distribution characteristics of industrial data, use the 3σ criterion to eliminate abnormal data of the stability parameters, and use the mean of the normal operating data as the dynamic characteristics identification result of the rotor system.

2. The rotor damping ratio identification method based on keyphasor sensing measurement according to claim 1, characterized in that The S1 comprises the following steps: S1.

1. Synchronously collect the key phase signal of the measuring plane where the rotating mechanical key phase sensor is located at a sampling frequency of not less than 2 kHz; perform preliminary filtering on the collected key phase signal, using a bandpass filter to remove high-frequency noise and some low-frequency interference. The cutoff frequency of the bandpass filter is set according to the spectral distribution characteristics of the key phase signal; S1.2, performing fast Fourier transform (FFT) on the key phase signal, calculating the square of the modulus of the complex number result in the frequency domain, and obtaining the key phase signal auto-power spectrum density function that characterizes the vibration characteristics of the rotor system through unilateral spectrum normalization.

3. A rotor damping ratio identification method based on keyphasor sensing measurement according to claim 1, characterized in that, The S2 comprises the following steps: S2.1, adding a specific number of data points at the starting and ending positions of the bond phase signal self-power spectrum density function, forming a dynamic sliding window with each original data point as the center, each dynamic sliding window is represented as each row of the matrix, and constitutes a window matrix of the bond phase signal self-power spectrum density function; S2.2, sorting the power spectrum density window matrix of the key phase signal and taking the middle value to obtain the median of each dynamic sliding window, and calculating the absolute difference between each element in the window and the middle value in the window, and taking the median of these absolute differences as the median absolute deviation (MAD); S2.

3. Determine outliers in the power spectral density window matrix of the key phase signal by comparing the differences between the elements in the sliding window and the median with the product of the absolute median difference (MAD) and the threshold according to the preset threshold; replace the determined outliers with the median of the corresponding window, thereby batch eliminating the harmonic interference in the key phase signal and obtaining the power spectral density function without harmonics.

4. A method for identifying the rotor damping ratio based on keyphasor sensing measurement according to claim 1, characterized in that, The said S4 includes the following steps: S4.

1. Respectively perform inverse autoregressive modal parameter identification on the time-domain pulse decay signal to obtain the characteristic roots of the discrete models of each order, and calculate the natural frequencies and damping ratios corresponding to each order in the measured key phase signal based on the characteristic roots; S4.

2. Use the 3σ clustering algorithm to automatically cluster the natural frequencies and damping ratios of different orders in the key phase signal; through cluster analysis, group similar dynamic characteristic parameters into one category, thereby realizing the efficient identification and classification of the dynamic characteristics of the rotor system.

Citation Information

Patent Citations

  • Rotary machine stability identification method and device, computer equipment and storage medium

    CN111289275A

  • Distributed spaceborne SAR synchronous phase estimation method

    CN118584446A

  • Classification apparatus, classification method, and classification program

    JP2022186422A

  • Method of Processing a Noisy Sound Signal and Device for Implementing Said Method

    US20070255535A1