Method for estimating rotational speed of multi-shaft rotating machinery based on displacement vibration signal

Abstract

This application belongs to the field of rotating machinery signal processing technology, and discloses a method for estimating the rotational speed of multi-axis rotating machinery based on displacement vibration signals. The method involves acquiring displacement vibration signals and performing short-time Fourier transform to obtain a time-frequency distribution matrix; selecting a safe harmonic set for each axis system based on the estimated rotational speed range; calculating the sharpness evaluation value of each safe harmonic and selecting the optimal harmonic as the target tracking harmonic; constructing a cost function that fuses amplitude and frequency mutation penalties; performing ridge tracking in the denoised local time-frequency matrix to obtain a preliminary estimated time-frequency curve; constructing a demodulation phase function based on this curve; performing generalized demodulation and narrowband filtering on the original signal to extract clean harmonic components; and finally reconstructing the final rotational speed curve of the corresponding axis system through instantaneous phase differentiation. This method achieves high-precision, anti-interference decoupling and estimation of the rotational speed of each axis system under strong interference conditions.

Description

Technical Field

[0001] This application belongs to the field of rotating machinery signal processing technology, specifically relating to a method for estimating the rotational speed of multi-axis rotating machinery based on displacement vibration signals. Background Technology

[0002] Rotating machinery, such as aero engines, gas turbines, and various multi-axis gearboxes, relies on rotational speed as a key parameter for characterizing operational status, fault diagnosis, and health management. Traditional speed measurement depends on direct contact or non-contact sensors, such as photoelectric encoders or magnetoelectric tachometers. However, in extreme service environments such as high temperature, high pressure, and high speed, sensors suffer from installation difficulties, susceptibility to damage, and reduced reliability. Furthermore, in compact or enclosed mechanical systems, space constraints may prevent the deployment of sensors. Therefore, tachometer-free methods for extracting rotational speed information from readily available vibration signals have become an important research direction in the field of rotating machinery condition monitoring.

[0003] Existing speed estimation methods without tachometers are typically based on time-frequency analysis of vibration signals. For a single-rotor system, its rotational frequency and its harmonics are relatively sparsely distributed in the frequency spectrum. By directly tracing the energy ridges in the time-frequency spectrum, the rotational speed can be estimated relatively accurately. These methods assume that there is a clear and singular correspondence between the dominant frequency components of the vibration signal and the rotor speed.

[0004] However, practical industrial rotating machinery often employs a multi-shaft design. For example, a dual-rotor aircraft engine comprises a high-speed rotating high-pressure rotor and a low-speed rotating low-pressure rotor. Each rotor, rotating independently, excites a vibration response containing harmonics of its rotational frequency. When multiple rotors operate simultaneously, their vibration signals linearly superimpose at the sensor measurement points. The rotational frequencies of different rotors and their integer harmonic components approach or even overlap each other in the frequency domain, forming complex harmonic cross-interference phenomena. Furthermore, background noise, structural resonance, and transient impact components further contaminate the spectrum. In this context, a single, clear, and continuous energy ridge corresponding to a specific rotor no longer exists in the time-frequency spectrum. Directly applying ridge-tracking methods for single-rotor systems easily locks onto the frequency components of interfering shafts or erroneous noise peaks, leading to severe deviations or even complete failure of the speed estimation results.

[0005] In summary, when processing vibration signals from multi-axis rotating machinery, existing technologies suffer from overlapping and intertwining of harmonic components generated by different rotors in the frequency domain, resulting in severe cross-interference and spectral aliasing. This makes it impossible to reliably separate and extract the rotational speed information of a single target shaft system. Summary of the Invention

[0006] To address the problem of inaccurate speed extraction in multi-axis rotating machinery due to the cross-interference of harmonics from different rotors in existing technologies, this application provides a speed estimation method for multi-axis rotating machinery based on displacement vibration signals. By integrating harmonic filtering, sharpness evaluation, robust ridge tracking, and generalized demodulation phase analysis, this method achieves high-precision, anti-interference decoupling and estimation of the speed of each shaft system under strong interference background, providing a reliable technical means for condition monitoring without tachometers.

[0007] To achieve the above technical objectives, this application specifically adopts the following technical solution:

[0008] In one aspect of this application, a method for estimating the rotational speed of multi-axis rotating machinery based on displacement vibration signals is provided, comprising the following steps:

[0009] S1. Signal Acquisition and Time-Frequency Transformation: The displacement vibration signal of the rotating machinery is acquired, and the displacement vibration signal is subjected to short-time Fourier transform to obtain a time-frequency distribution matrix. The time-frequency distribution matrix contains the distribution of the displacement vibration signal amplitude on the time-frequency plane.

[0010] S2. Cross-interference suppression and safety harmonic screening: Based on the estimated rotational speed range of each shaft system of the rotating machinery, calculate the theoretical harmonic frequency of each shaft system; by comparing the closeness of the theoretical harmonic frequencies of different shaft systems in the corresponding frequency band of the time-frequency distribution matrix, screen and eliminate harmonics with cross-interference, and determine the safety harmonic set of each shaft system.

[0011] S3. Optimal Harmonic Selection Based on Sharpness Evaluation: For each candidate harmonic in the safety harmonic set of each axis system, a local time-frequency matrix centered on the theoretical frequency of the candidate harmonic is extracted from the time-frequency distribution matrix. The sharpness evaluation value of each local time-frequency matrix is ​​calculated. The sharpness evaluation value comprehensively reflects the energy concentration and bandwidth of the local time-frequency matrix. The candidate harmonic with the highest sharpness evaluation value is selected as the target tracking harmonic of the corresponding axis system.

[0012] S4. Ridge tracking based on cost function: For the target tracking harmonics, the corresponding local time-frequency matrix is ​​extracted from the time-frequency distribution matrix and denoised; a cost function that integrates the normalized amplitude term and the frequency change penalty term is constructed. In the denoised local time-frequency matrix, ridge tracking is performed by minimizing the cost function to obtain the preliminary estimated frequency curve of the corresponding axis system.

[0013] S5. Generalized demodulation and speed reconstruction: Based on the preliminary estimated rotational frequency curve, a demodulation phase function is constructed, and the displacement vibration signal is subjected to generalized demodulation and narrowband filtering to extract the pure harmonic components corresponding to the target tracking harmonics; the instantaneous phase is extracted from the pure harmonic components, and the instantaneous phase is differentiated to reconstruct the final speed curve of the corresponding shaft system.

[0014] In one implementation, step S1 of performing a short-time Fourier transform on the displacement vibration signal includes: filling the beginning and end of the displacement vibration signal with zeros, windowing it using a rectangular window function, and performing a sliding transform with a preset overlap rate to obtain a complex-form time-frequency distribution matrix. :

[0015]

[0016] in, Indicates the time frame index. Indicates frequency index, The complex value at time frame m and frequency index k represents the amplitude and phase information of the signal at that time and frequency point; Index of time samples within the window; Denotes the step size, and satisfies , For window length, The number of overlapping points; For discrete displacement vibration signals after zero-value filling, in time and position Samples at the location; For length is The rectangular window function; This is the kernel function for the Discrete Fourier Transform, used to map time-domain signals to the frequency domain. It is an imaginary number;

[0017] By taking the modulus of the time-frequency distribution matrix, an amplitude spectrum matrix reflecting the time-frequency variation of the displacement vibration signal amplitude is obtained. .

[0018] In one implementation, step S2, screening and removing harmonics with cross-interference, includes: calculating the relative deviation ratio between the theoretical frequencies of any two harmonics in different axis systems; if the relative deviation ratio is less than a preset interference determination threshold, then the two corresponding harmonics are determined to have cross-interference and are removed.

[0019] In one implementation, in step S3, the sharpness evaluation value for:

[0020]

[0021] in, Let h be the spectral quality of the low-h-th order candidate harmonic of the i-th rotation axis. Weighted by frequency proximity; Energy concentration; For average normalized bandwidth, To prevent zero constant.

[0022] In one implementation, the energy concentration for:

[0023]

[0024] in, Let be the local time-frequency matrix corresponding to the i-th rotating axis and the h-th harmonic. For time frame indexing, This is a local frequency index.

[0025] In one implementation, step S4 includes: reconstructing the local time-frequency matrix corresponding to the target tracking harmonic into a time-domain narrowband signal, constructing the Hankel matrix of the time-domain narrowband signal and performing singular value decomposition, reconstructing the denoised matrix by retaining the main singular values ​​according to a preset energy threshold, and then converting the denoised matrix back to the time-frequency domain.

[0026] In one implementation, in step S4, the cost function for:

[0027]

[0028]

[0029]

[0030] in, Candidate frequency index for the current frequency point In time frame The normalized amplitude, Based on the current candidate frequency index Ridge frequency index from the previous moment Distance between Calculated distance penalty term, and These are the positive weighting coefficients.

[0031] In one implementation, the distance penalty term for:

[0032]

[0033] in, Bandwidth parameters for controlling penalty sensitivity.

[0034] In one implementation, the ridge tracking in step S4 employs a bidirectional tracking strategy, including forward tracking from the first time frame to the last time frame, and reverse tracking from the last time frame to the first time frame with path correction.

[0035] In one implementation, the demodulation phase function described in step S5 for:

[0036]

[0037] in, For the i-th axis at time Smooth estimation of the frequency of change, Let be the optimal harmonic order of the i-th rotating axis. For instantaneous frequency estimation of the target harmonic, The average frequency of the target harmonic.

[0038] In one implementation, in step S5, the final rotational speed curve of the corresponding shaft system is reconstructed. for:

[0039]

[0040] in, Let be the final estimated rotational speed of the i-th shaft at time t. Let be the instantaneous frequency of the i-th rotating axis at time t. Let be the optimal harmonic order of the i-th rotating axis.

[0041] The beneficial effects of this application are as follows:

[0042] 1) Effectively overcome harmonic cross-interference: By pre-calculating and eliminating theoretical harmonics with frequency overlap between different shaft systems, an independent set of safe harmonics for each shaft is constructed, which avoids tracking errors caused by spectral aliasing from the source and significantly improves the reliability of speed decoupling in multi-rotor systems.

[0043] 2) Achieve high-precision speed reconstruction: By introducing a clarity index that integrates energy concentration and bandwidth, the harmonic component with the best signal-to-noise ratio and the clearest time-frequency ridge can be automatically selected as the tracking object; combined with generalized demodulation technology, the phase information is directly extracted and differentiated, so that the final reconstructed speed curve can reflect transient micro-perturbations, achieving an accuracy that surpasses traditional time-frequency ridge tracking methods.

[0044] 3) Strong robustness and engineering applicability: Ridge tracking employs a cost function that includes amplitude attraction and frequency jump penalty, effectively suppressing frequency jumps caused by noise and transient impacts, ensuring the continuity of the estimated curve. The method in this application has a clear flow, does not rely on tachometer signals, and is suitable for online or offline analysis under complex operating conditions. Attached Figure Description

[0045] Figure 1 This is a flowchart illustrating the method for estimating the rotational speed of multi-axis rotating machinery based on displacement vibration signals, as described in this application.

[0046] Figure 2 These are discrete displacement vibration signals collected in the embodiments of this application;

[0047] Figure 3 The rotation phases of shaft 1 and shaft 2 are predicted based on vibration signals in this embodiment of the application; where hollow dots represent the rotation phases measured by the photoelectric speed sensor, and solid lines represent the rotation phases predicted in this embodiment.

[0048] Figure 4 The rotational speeds of shaft 1 and shaft 2 are predicted based on vibration signals in this embodiment of the application; where hollow dots represent the rotational speeds measured and calculated by the photoelectric speed sensor, and solid lines represent the rotational speeds predicted in this embodiment. Detailed Implementation

[0049] The technical solution of this application will be clearly and completely described below with reference to specific embodiments. However, those skilled in the art will understand that the embodiments described below are only some embodiments of this application, not all embodiments, and are only used to illustrate this application, and should not be regarded as limiting the scope of this application. Based on the embodiments in this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0050] To address the challenge of overlapping and difficult-to-separate rotor harmonics in vibration signals of multi-axis rotating machinery, this application does not blindly track harmonics across the entire time-frequency plane. Instead, based on prior system information, it actively identifies and eliminates interfering frequency components originating from other shaft systems, defining a "safe" harmonic selection range without cross-influence for each target rotor. Furthermore, by quantitatively evaluating the energy concentration and contour clarity of each candidate harmonic in the time-frequency domain, the component with the highest characterizing quality is selected for tracking. To cope with noise and transient changes in the actual signal, an optimization algorithm combining amplitude characteristics and path continuity is employed to ensure the stability of frequency ridge tracking. Finally, abandoning the approach of simply relying on ridge frequency values, generalized demodulation technology is used to convert the selected harmonic components into steady-state signals, and the rotational speed is reconstructed by extracting and differentiating their instantaneous phase. This achieves robust and high-precision decoupling and extraction of the rotational speed information of each independent rotor from a strongly interfering mixed signal.

[0051] Exemplary embodiments

[0052] This embodiment uses a dual-axis rotor test bench as a simulation object for multi-axis rotating machinery to collect and verify data.

[0053] The dual-axis rotor test bench adopts a coaxial structure, including an inner rotating shaft (shaft 1) and an outer rotating shaft (shaft 2). The inner rotating shaft passes through the center of the hollow outer rotating shaft, and the two shafts are connected and coupled to each other via an intermediate bearing. Each end of the inner rotating shaft is supported by two independent bearing seats, while one end of the outer rotating shaft is supported by the intermediate bearing, and the other end is supported by an independent bearing seat. The inner and outer rotating shafts are each driven by two independent drive motors to achieve independent or coupled rotation of the dual-rotor system.

[0054] In this embodiment, during data acquisition, the inner and outer rotating shafts are controlled to maintain a constant rotational speed, and the sampling duration is set to 3 seconds. Radial vibration signals are acquired using an eddy current displacement sensor. Displacement time-domain vibration signals are obtained without relying on a tachometer.

[0055] Let the collected discrete displacement vibration signal be... The sampling frequency is The number of sampling points is The displacement vibration signal is represented as a one-dimensional time series vector in the computer:

[0056]

[0057] in, Represents a discrete-time index. This time series serves as the input data for subsequent time-frequency analysis and rotational speed estimation.

[0058] Reference Figure 1 As shown, the specific process is as follows:

[0059] 1) Signal preprocessing and time-frequency transformation

[0060] The collected discrete displacement vibration signals ( Figure 2 A short-time Fourier transform (STFT) is performed to obtain the time-frequency distribution containing the time-varying characteristics of the rotational speed. Considering the boundary effects and frequency resolution requirements in practical signal processing, the specific processing steps in this embodiment are as follows:

[0061] a. Signal filling

[0062] To reduce the boundary truncation effect during subsequent time-frequency analysis, the displacement vibration signal is first analyzed. The beginning and end of the window are padded with zero values. The sliding window length is set to... Then the fill length Set to half the window length, i.e. , This represents the rounding function. The filled signal. Increasing the length ensures that the data at the boundaries of the original signal are also located in the center of the sliding window, thus allowing for complete analysis.

[0063] b. Sliding window settings

[0064] Based on the code logic, this embodiment selects a rectangular window function. As a truncation function. Although the rectangular window has a larger side lobe leakage, it has an advantage in main lobe width, which is beneficial for resolving adjacent frequency components in the frequency-dense region of the dual rotor. Window length The power of 2 was chosen to accelerate the calculation process of the Fast Fourier Transform (FFT).

[0065] Set the window overlap rate (In this embodiment) Then the number of overlapping points A higher overlap rate helps improve resolution on the time axis, thereby capturing instantaneous fluctuations in rotational speed.

[0066] c. Short-Time Fourier Transform (STFT) Calculation

[0067] Using the rectangular window function set above For the filled signal Windowing and discrete Fourier transform are performed to obtain the complex form of the time-frequency distribution matrix. The calculation formula is as follows:

[0068]

[0069] in, Indicates the time frame index. Indicates frequency index, The complex value at time frame m and frequency index k represents the amplitude and phase information of the signal at that time and frequency point; Index of time samples within the window; Denotes the step size, and satisfies ; For discrete displacement vibration signals after zero-value filling, in time and position Samples at the location; This is the kernel function for the Discrete Fourier Transform, used to map time-domain signals to the frequency domain. It is an imaginary number.

[0070] d. Time-frequency spectrum generation

[0071] The calculated complex time-frequency distribution matrix Taking the modulus, we obtain the amplitude spectrum matrix. The amplitude spectrum matrix not only reflects the distribution of signal energy in the time-frequency plane, but also includes the fundamental frequency and its harmonic components of each rotor in the multi-axis system.

[0072]

[0073] Received This serves as the foundational input data for subsequent steps involving harmonic crosstalk suppression, sharpness evaluation, and ridgeline tracing.

[0074] 2) Multi-axis harmonic cross-interference detection and safe harmonic screening

[0075] Because multi-axis rotating machinery contains multiple independently rotating rotors, and the displacement vibration signal is affected by the coupling effect of the vibrations of each rotor, the rotational frequencies of different rotors and their harmonic harmonics are prone to aliasing in the frequency spectrum. In order to avoid erroneously tracking the frequency components of interfering shafts in subsequent steps, harmonic orders that pose a risk of cross-interference are pre-emptively eliminated by calculating the theoretical frequency trajectory.

[0076] a. Constructing a set of theoretical harmonic frequencies

[0077] First, based on the preset operating conditions of the multi-axis rotating machinery, determine the nominal speed of each shaft system. Assume the shaft system includes... One rotating shaft (in this embodiment) ), No. The nominal rotational speed of each shaft is (Unit: RPM)

[0078] Therefore, the nominal fundamental frequency of the i-th rotating shaft is calculated. (Unit: Hz):

[0079]

[0080] Define the range of harmonic orders to be searched. For the i-th axis of rotation, define an initial set of harmonic orders. The set contains a series of integers from low to high order (e.g., taking the order). Based on the nominal fundamental frequency, calculate the theoretical frequency value of the k-th harmonic of the i-th rotating shaft. :

[0081] .

[0082] The theoretical frequencies of all axes and all orders constitute the initial set of candidate harmonics.

[0083] b. Harmonic Cross-Interference Determination

[0084] When determining the safe harmonic set for each rotating shaft, interference is judged based on the calculated theoretical frequencies of each candidate harmonic.

[0085] For any two different axes of rotation i and j ( ), and their respective arbitrary harmonic orders k and ( ), calculate the relative deviation ratio of the theoretical frequencies of the two. .

[0086] To accurately measure frequency similarity, the relative deviation ratio is defined as the ratio of the absolute value of the difference between two theoretical frequencies to the smaller of the two frequencies:

[0087] .

[0088] Set an interference detection threshold (In this embodiment, it is set) That is, the minimum allowable frequency interval is 0.5%.

[0089] Decision logic: If Then determine the k-th harmonic of the i-th rotating axis and the j-th rotating axis. The first harmonics pose a risk of cross-interference in the frequency domain.

[0090] Marking process: Once an interference risk is determined, the harmonic order k is marked as "unsafe", indicating that it is at risk of being interfered with and should be removed from subsequent processing; conversely, if the relative deviation ratio of a harmonic with any harmonic of all other axes is greater than the threshold, it is marked as "safe".

[0091] c. Generate a safe harmonic set

[0092] After removing the harmonic orders that have cross-interference, it is also necessary to remove extremely low-order harmonics with weak energy (such as the low-frequency noise region near the fundamental frequency) and extremely high-order harmonics that exceed the sampling frequency limit, based on the signal-to-noise ratio characteristics of the actual displacement vibration signal.

[0093] Finally, for the i-th axis of rotation, the selected retained orders constitute the safe harmonic set, which is denoted as […]. Any order in this set All meet the conditions of no cross-interference and being within the effective frequency band.

[0094]

[0095] in, and These are the upper and lower limits of the order set based on the sensor bandwidth. This indicates harmonics that are marked as safe. Output This will serve as a candidate set for selecting target tracking harmonics in subsequent steps.

[0096] 3) Selection of optimal harmonic components based on sharpness index

[0097] The safe harmonic set obtained after screening Although interference caused by multi-axis crossovers has been eliminated, the remaining harmonic components still exhibit significant differences in energy intensity, noise interference level, and bandwidth. To ensure the robustness of subsequent ridge tracking, this step introduces a comprehensive evaluation metric: sharpness. The harmonic with the clearest spectral characteristics and the most concentrated energy is selected from the safe harmonic set as the optimal tracking target.

[0098] a. Adaptive extraction of local time-frequency features

[0099] For each candidate safe harmonic order of the i-th axis of rotation First, in the generated amplitude spectrum matrix Extract the corresponding local time-frequency data.

[0100] Based on the theoretical center frequency of this harmonic As a baseline, set the initial search bandwidth. To correct the deviation between the actual rotational speed and the nominal rotational speed, a two-step search method is used to locate the actual frequency band:

[0101] The first step is to define the scope. Find the actual frequency position corresponding to the energy maxima in the amplitude spectrum matrix.

[0102] The second step involves using the initial frequency location as the new center frequency, and then extracting the bandwidth from the amplitude spectrum matrix again. Time-frequency data blocks.

[0103] The extracted data blocks ultimately constitute the local time-frequency matrix corresponding to the candidate harmonic, denoted as . Where m is the time frame index. This is a local frequency index. This matrix is ​​used for subsequent calculations of the sharpness rating for this harmonic.

[0104] b. Calculation of Feature Indicators

[0105] For the extracted local time-frequency matrix Calculate the following three key characteristic indicators respectively:

[0106] Indicator 1: Energy Concentration

[0107] Energy concentration reflects the proportion of the main harmonic peak energy to the total energy of a local frequency band. A higher proportion indicates a better signal-to-noise ratio. The calculation formula is:

[0108]

[0109] in, To prevent tiny positive numbers with a denominator of zero.

[0110] Indicator 2: Average Normalized Bandwidth

[0111] The sharpness of the spectral lines is evaluated using the full width at half maximum (FWHM) criterion. For each time frame m, the amplitude on both sides of the main peak decreases to the peak value. Calculate the instantaneous bandwidth of the time frame at the specified frequency point. The instantaneous bandwidth of all time frames is averaged and normalized to obtain the average normalized bandwidth:

[0112]

[0113] in, This represents the total number of frequency points in the local frequency band. This represents the total number of frames for the local time-frequency data block on the time axis. A smaller average normalized bandwidth value indicates sharper spectral lines and higher frequency resolution.

[0114] Indicator 3: Frequency proximity weight

[0115] To prevent incorrectly locking onto noise peaks of non-target axes, the deviation between the actual extracted center frequency and the theoretical center frequency is calculated. If the deviation is within acceptable limits, a weight is set. If the deviation is too large, it indicates that the component may be a false signal. Set... .

[0116] c. Sharpness scoring and optimal harmonic determination

[0117] Based on the calculated feature indicators, a comprehensive clarity evaluation function is constructed. This is used to quantify the tracking suitability of each candidate harmonic. Sharpness rating value. for:

[0118] .

[0119] in, Let be the spectral quality of the low h-th order candidate harmonic of the i-th axis of rotation. This function indicates that the more concentrated the energy and the narrower the spectral bandwidth (the sharper the spectral line), the higher the sharpness evaluation value.

[0120] Iterating through the safe harmonic set of the i-th axis of rotation For each order h, calculate its respective sharpness rating. The candidate harmonic with the highest evaluation value is selected as the optimal harmonic order for this rotating shaft. :

[0121]

[0122] At the same time, the local time-frequency matrix corresponding to the optimal harmonic is retained. Its frequency boundary index is used as input data for subsequent ridgeline tracing steps.

[0123] 4) SVD Denoising and Ridge Tracing Based on Cost Function

[0124] In actual operating conditions, displacement vibration signals contain background noise and transient impact components, which may cause the energy ridges in their time-frequency distribution to become broken or blurred. To obtain a continuous, smooth, and accurate instantaneous frequency trajectory, this step processes the local time-frequency matrix corresponding to the target harmonic. The processing consists of two main parts: first, the local time-frequency matrix is ​​denoised to enhance the effective components; then, a specific cost function is constructed, and the optimal frequency path is traced in the denoised time-frequency data by minimizing this function.

[0125] a. Local band enhancement based on SVD

[0126] First, the local time-frequency matrix corresponding to the extracted optimal harmonic is... Transform back to the time domain. Using the inverse short-time Fourier transform (ISTFT), reconstruct the complex matrix into a narrowband time-domain waveform signal, denoted as . .

[0127] To effectively separate noise and periodic components in this time-domain narrowband waveform signal, a method is constructed. Hankel matrix .set up Length is Take the embedding dimension as (usually taken) (left and right), then the Hankel matrix is ​​constructed as follows:

[0128]

[0129] in, For the matrix Perform Singular Value Decomposition (SVD):

[0130]

[0131] in, Let be a diagonal matrix composed of singular values, and ; It is a left singular vector matrix; This is a right singular vector matrix. Based on a preset energy threshold criterion (e.g., retaining the frontmost vectors with an energy percentage exceeding 80%)... (single values), reconstructed denoised Hankel matrix .

[0132]

[0133] in, They are matrices The List; will Restored to a one-dimensional time domain signal Then, it is subjected to STFT transformation again to obtain a clean local time-frequency matrix with significantly suppressed background noise, denoted as . This is used for subsequent ridgeline tracing calculations.

[0134] b. Cost function construction

[0135] In order to To accurately track frequency ridges, a path cost function combining amplitude advantage and continuity constraint is defined. .

[0136] For any candidate frequency index of time frame m Its cost function consists of two parts:

[0137] ① Normalized amplitude term ( ): It is expected that the ridge line is located at the point of strongest energy.

[0138] Let the maximum amplitude of the spectrum at the current time be... The minimum amplitude is Then the current candidate frequency index normalized amplitude for:

[0139]

[0140] The amplitude cost term is: ,in, This is a custom positive weighting coefficient, typically 1.

[0141] ② Distance penalty section :

[0142] To suppress irregular jumps in the time-frequency ridge line caused by noise and interference, and to ensure the temporal continuity of the estimated frequency trajectory, an optimization mechanism based on frequency mutation penalty is introduced in the ridge line tracing step.

[0143] Let the previous moment be The determined ridge frequency index is Current candidate frequency index The distance to it is The distance penalty term is defined using the Gaussian function form. Quantify the transition:

[0144]

[0145] in, Bandwidth parameters for controlling penalty sensitivity.

[0146] This distance penalty With a positive weighting coefficient Multiplication forms the distance penalty component in the total cost function. :

[0147]

[0148] in, The value is usually 3.

[0149] In summary, the total cost function for:

[0150] .

[0151] c. Two-way ridge tracking

[0152] To further improve the robustness of tracking, this embodiment adopts a bidirectional tracking strategy:

[0153] Forward tracing: from the start time Recursion to the termination time At each moment, iterate through all possible... Select to make The lowest frequency point is used as the ridge point at the current moment to update the path.

[0154] Reverse optimization: Starting from the endpoint of the forward tracing, from... Recursion Utilize the cost function again Perform path correction.

[0155] Finally, the precise frequency ridge index sequence of the i-th rotating shaft at the optimal harmonic is obtained. Map the index back to the physical frequency value, add the starting offset of the local frequency band, and divide by the optimal harmonic order. This allows us to obtain a preliminary estimate of the rotational frequency of the i-th shaft. .

[0156]

[0157] in, This is a conversion function from frequency index to Hertz.

[0158] 5) Generalized demodulation and phase differential speed reconstruction

[0159] Although step 4) obtains a preliminary estimate of the rotational speed curve through ridge tracing, this curve is based on a time-frequency grid derived from a short-time Fourier transform (STFT). Its frequency resolution is limited by the window length, and it exhibits a stepped discontinuity on the time axis. To obtain a high-precision, continuous rotational speed curve that reflects transient micro-fluctuations, this step utilizes generalized demodulation (GD) technology to frequency-shift the original displacement vibration signal, converting the time-varying frequency band signal into a narrowband signal with a stable center frequency. Then, through narrowband filtering and instantaneous phase differentiation, the final rotational speed curve is reconstructed.

[0160] a. Frequency ridge smoothing and demodulation operator construction

[0161] First, the initial estimated frequency sequence for the i-th output axis. Curve fitting and smoothing are performed. Polynomial fitting or spline smoothing algorithms are used to obtain a time-continuous smoothed frequency function. .

[0162] To convert the time-varying harmonic components corresponding to this shaft into signals with a constant center frequency, a phase demodulation operator needs to be constructed. The theoretical mean center frequency of this harmonic component is then calculated. :

[0163]

[0164] in, The total duration of the signal. The selected optimal harmonic order.

[0165] Next, the demodulation phase function used to compensate for frequency fluctuations is calculated. This function is defined as the time integral of the difference between the instantaneous harmonic frequency and the mean frequency:

[0166] .

[0167] b. Generalized demodulation transform and narrowband filtering

[0168] Using Hilbert transform to transform the original real-value displacement vibration signal Constructed as a complex analytic signal :

[0169]

[0170] in, is the Hilbert transform; j is the imaginary unit.

[0171] Using the constructed demodulation phase function For analytical signals Perform a generalized demodulation operation. This is achieved by introducing a reverse rotation factor. This cancels out the frequency fluctuations of the target harmonic components, making them appear in the transformed domain with a center frequency of The stable signal. The demodulated signal is denoted as... :

[0172] .

[0173] right The spectrum is obtained by performing a Fourier transform. At this point, a center frequency of [value] is designed. Bandpass filters with extremely narrow bandwidth (For example, a rectangular window filter or a Gaussian filter) to accurately filter out background noise and interference components from other axes:

[0174]

[0175] The filtered spectrum Perform an inverse Fourier transform and multiply by a positive rotation factor. The restoration yields a pure analytical signal in the time domain containing only the optimal harmonic component of the i-th rotation axis. :

[0176]

[0177] in, This represents the inverse Fourier transform.

[0178] c. Phase extraction and speed reconstruction

[0179] The pure analytical signal obtained from the i-th rotating axis Extracting instantaneous phase To solve the phase winding problem, the phase is unwound:

[0180]

[0181] in, This indicates the unwinding operation; This indicates the calculation of the phase angle of a complex number. The output of this embodiment is... like Figure 3 As shown.

[0182] By performing time differentiation on the continuous phase curve after unwinding, the precise instantaneous frequency of the optimal harmonic component is obtained. :

[0183]

[0184] Finally, based on the optimal harmonic order By inversely calculating the ratio of the rotational frequency to the rotational frequency of the i-th shaft, the final high-precision rotational speed of the i-th shaft can be obtained. (Unit: rad / s):

[0185]

[0186] The output of this embodiment like Figure 4 As shown in the figure. Through the above steps, this application successfully achieved precise decoupling and reconstruction of the rotational speeds of each rotor in a strong interference and multi-axis environment.

[0187] Although the embodiments of this application have been described above in conjunction with the accompanying drawings, this application is not limited to the specific embodiments and application fields described above. The specific embodiments described above are merely illustrative and instructive, not restrictive. Those skilled in the art can make many other forms based on the guidance of this specification and without departing from the scope of protection of the claims of this application, and all of these are within the scope of protection of this application.

Claims

1. A method for estimating the rotational speed of multi-axis rotating machinery based on displacement vibration signals, characterized in that, Includes the following steps: S1. Collect the displacement vibration signal of the rotating machinery, and perform a short-time Fourier transform on the displacement vibration signal to obtain a time-frequency distribution matrix. The time-frequency distribution matrix contains the distribution of the displacement vibration signal amplitude on the time-frequency plane. S2. Based on the estimated rotational speed range of each shaft system of the rotating machinery, calculate the theoretical harmonic frequency of each shaft system; by comparing the closeness of the theoretical harmonic frequencies of different shaft systems in the frequency band corresponding to the time-frequency distribution matrix, screen out and eliminate harmonics with cross-interference, and determine the safe harmonic set of each shaft system. S3. For each candidate harmonic in the safety harmonic set of each axis system, extract the local time-frequency matrix centered on the theoretical frequency of the candidate harmonic from the time-frequency distribution matrix, and calculate the sharpness evaluation value of each local time-frequency matrix. The sharpness evaluation value comprehensively reflects the energy concentration and bandwidth of the local time-frequency matrix. The candidate harmonic with the highest sharpness evaluation value is selected as the target tracking harmonic for the corresponding axis system; The clarity evaluation value for: in, The spectral quality of the low h-th order candidate harmonic of the i-th rotation axis; The frequency proximity weight represents the deviation between the actual extracted center frequency and the theoretical center frequency. Energy concentration; For average normalized bandwidth, To prevent zero constant; S4. For the target tracking harmonics, extract the corresponding local time-frequency matrix from the time-frequency distribution matrix and perform noise reduction processing; construct a cost function that integrates the normalized amplitude term and the frequency change penalty term, and perform ridge tracking by minimizing the cost function in the noise-reduced local time-frequency matrix to obtain the preliminary estimated frequency curve of the corresponding axis system. S5. Based on the preliminary estimated rotational frequency curve, construct a demodulation phase function, perform generalized demodulation and narrowband filtering on the displacement vibration signal, and extract the pure harmonic components corresponding to the target tracking harmonics; extract the instantaneous phase from the pure harmonic components, differentiate the instantaneous phase, and reconstruct the final rotational speed curve of the corresponding shaft system.

2. The method for estimating the rotational speed of multi-axis rotating machinery based on displacement vibration signals according to claim 1, characterized in that, Step S1, performing a short-time Fourier transform on the displacement vibration signal, includes: filling the beginning and end of the displacement vibration signal with zeros, windowing it using a rectangular window function, and performing a sliding transform with a preset overlap rate to obtain a complex time-frequency distribution matrix. : in, Indicates the time frame index. Indicates frequency index, The complex value at time frame m and frequency index k represents the amplitude and phase information of the signal at that time and frequency point; Index of time samples within the window; Denotes the step size, and satisfies , For window length, The number of overlapping points; For discrete displacement vibration signals after zero-value filling, in time and position Samples at the location; For length is The rectangular window function; This is the kernel function for the Discrete Fourier Transform, used to map time-domain signals to the frequency domain. It is an imaginary number; By taking the modulus of the time-frequency distribution matrix, an amplitude spectrum matrix reflecting the time-frequency variation of the displacement vibration signal is obtained. .

3. The method for estimating the rotational speed of multi-axis rotating machinery based on displacement vibration signals according to claim 1, characterized in that, The step S2 of screening and removing harmonics with cross-interference includes: calculating the relative deviation ratio between the theoretical frequencies of any two harmonics in different axis systems; if the relative deviation ratio is less than a preset interference judgment threshold, then the two corresponding harmonics are determined to have cross-interference and are removed.

4. The method for estimating the rotational speed of multi-axis rotating machinery based on displacement vibration signals according to claim 1, characterized in that, The energy concentration for: in, Let be the local time-frequency matrix corresponding to the i-th rotating axis and the h-th harmonic. For time frame indexing, This is a local frequency index.

5. The method for estimating the rotational speed of multi-axis rotating machinery based on displacement vibration signals according to claim 1, characterized in that, In step S4, the noise reduction process includes: reconstructing the local time-frequency matrix corresponding to the target tracking harmonic into a time-domain narrowband signal, constructing the Hankel matrix of the time-domain narrowband signal and performing singular value decomposition, reconstructing the noise-reduced matrix by retaining the main singular values ​​according to a preset energy threshold, and then converting the noise-reduced matrix back to the time-frequency domain.

6. The method for estimating the rotational speed of multi-axis rotating machinery based on displacement vibration signals according to claim 1, characterized in that, In step S4, the cost function for: in, Candidate frequency index for the current frequency point In time frame The normalized amplitude, Based on the current candidate frequency index Ridge frequency index from the previous moment Distance between Calculated distance penalty term, and These are the positive weighting coefficients.

7. The method for estimating the rotational speed of multi-axis rotating machinery based on displacement vibration signals according to claim 6, characterized in that, The distance penalty item for: in, Bandwidth parameters for controlling penalty sensitivity.

8. The method for estimating the rotational speed of multi-axis rotating machinery based on displacement vibration signals according to claim 1, characterized in that, The ridge tracking described in step S4 adopts a bidirectional tracking strategy, including forward tracking from the first time frame to the last time frame, and reverse tracking from the last time frame to the first time frame with path correction.

9. The method for estimating the rotational speed of multi-axis rotating machinery based on displacement vibration signals according to claim 1, characterized in that, The demodulation phase function described in step S5 for: in, For the i-th axis at time Smooth estimation of the frequency of change, Let be the optimal harmonic order of the i-th rotating axis. For instantaneous frequency estimation of the target harmonic, The average frequency of the target harmonic.

Images