Aero-engine strong non-stationary cross frequency signal time-frequency analysis method and system

Abstract

The invention provides an aero-engine strong non-stationary cross frequency signal time-frequency analysis method and system, and relates to the technical field of aero-engine signal processing and fault diagnosis, and the method comprises the steps: obtaining a vibration signal, and obtaining initial time-frequency distribution; performing adaptive ridge line extraction on the initial time-frequency distribution, and determining a time-frequency block range of an intersection region; carrying out local sub-blocking on the time frequency distribution in the intersection area to obtain optimized time frequency distribution; performing component reconstruction on the vibration signal by adopting Vold-Kalman order filtering to obtain final time frequency distribution; according to the method, the ridge line fracture time-frequency block is determined through bidirectional retrieval repair, the intersection area and the hierarchical cascade reconstruction component are optimized through comprehensive evaluation indexes, the ridge line continuity and instantaneous frequency extraction precision of the intersection area are remarkably improved, and high-resolution representation and accurate separation of strong non-stationary intersection frequency signals are achieved.

Description

Technical Field

[0001] This invention relates to the field of aero-engine signal processing and fault diagnosis technology, specifically to a time-frequency analysis method and system for strong non-stationary cross-frequency signals of aero-engines. Background Technology

[0002] As the core power unit of an aircraft, the operating status of the aero-engine directly affects flight safety. The compressor, a critical component of the aero-engine, generates vibration signals with strong non-stationary characteristics under complex operating conditions. Particularly during the operation of a dual-rotor system, the frequency components of the high- and low-pressure rotors often intersect in the time-frequency domain, forming typical crossover frequency signals. Accurate analysis of these signals is crucial for early fault diagnosis and anomaly identification.

[0003] Traditional time-frequency analysis methods face numerous challenges when processing cross-frequency signals. While the short-time Fourier transform is widely used, it suffers from severe energy dispersion in the time-frequency crossover region, making it difficult to clearly distinguish neighboring frequency components. Wavelet transform, limited by fixed time-frequency resolution, struggles to simultaneously provide a precise characterization of both the time and frequency domains. Synchronous compression transform enhances energy concentration through frequency domain redistribution, but the extracted time-frequency ridges in the crossover region often exhibit breaks, leading to inaccurate instantaneous frequency estimation and impacting subsequent feature extraction and fault diagnosis.

[0004] In recent years, methods such as local maximum synchronous compression transform have improved time-frequency resolution to some extent. However, when processing cross-frequency signals from dual rotors of aero-engines, problems such as insufficient ridge continuity in the cross-region and energy leakage near the instantaneous frequency trajectory still exist. The root cause of these problems is that existing methods lack a dedicated processing mechanism for the cross-region, and cannot effectively cope with the energy dispersion and ridge breakage caused by time-frequency cross-intersection.

[0005] In actual operation, aero-engine vibration signals exhibit significant multi-component and strongly non-stationary characteristics. As a crucial tool for analyzing non-stationary signals, time-frequency analysis (TFA) has become a core method for fault diagnosis and performance evaluation of complex rotating machinery. TFA can jointly characterize the energy distribution and instantaneous frequency features of a signal in both time and frequency dimensions. However, the high and low voltage rotational frequencies and their harmonics of a dual-rotor system often approach, overlap, or even actually cross during operation, resulting in severe overlap of local frequency band energy and characteristic confusion. Traditional TFA generally exhibits energy dispersion and ambiguous characterization in this "crossing neighborhood," making it difficult to meet the needs of accurate diagnosis.

[0006] Existing Time-Frequency Analysis (TFA) methods mainly fall into two categories: one is adaptive analysis based on multimodal separation, such as Empirical Mode Decomposition (EMD), Local Mean Decomposition (LMD), and Variational Nonlinear Chirp Mode Decomposition (VNCMD), which decompose complex signals into physically meaningful intrinsic mode functions and then extract instantaneous features. However, these methods are susceptible to mode aliasing, resulting in spurious components. The other category is analysis methods based on basis function transformation, such as Velocity Synchronous Linear Chirplet Transform (VSLCT) and Scaled Basis Chirplet Transform (SBCT). Although these methods improve the concentration of time-frequency energy to some extent, they still struggle to maintain high resolution in frequency crossover regions, leading to unclear instantaneous frequency extraction.

[0007] To further enhance time-frequency focusing capabilities, researchers have proposed various post-processing redistribution strategies, such as the redistribution method (RM) and Synchronous Compression Transform (SST) and its improved versions, including Multiple Synchronous Compression (MSST), Local Maximum SST (LMSST), Synchronous Redistribution Transform (SRT), and Local Maximum Synchronous Chirplet Transform (LMSSCT). Although these methods have made progress in energy redistribution and ridge extraction, most rely on the accuracy of the initial time-frequency representation and require that the signal components have been well separated in the time-frequency domain. Therefore, it is still difficult to obtain a complete characterization under frequency crossover conditions.

[0008] In recent years, research on cross-frequency signals has continued to advance, with methods such as Generalized Linear Chirplet Transform (GLCT), Adaptive Linear Chirplet Transform (ALCT), Generalized Modulation Basis Transform (GCBT), and Adaptive Linear Chirplet Synchronous Extraction Transform (ALCSET) emerging one after another. However, when dealing with cross-regions, they still face challenges such as insufficient time-frequency accuracy, computational complexity, or errors in selecting the optimal rotation angle. Entropy-matched Synchronous Compressed Chirplet Transform (EMSSCT) searches for the optimal Chirplet rate through Rényi entropy and combines it with a synchronous compression strategy, showing outstanding performance in global energy accumulation of strong non-stationary signals. However, it is highly dependent on the initial EMCT results and is prone to frequency redistribution misjudgment at cross-regions, leading to ridge breakage and feature loss, which limits its application in the complex operating conditions of aero-engines.

[0009] Chinese patent document CN114674410A discloses a method for instantaneous frequency estimation of underwater acoustic signals with time-varying component numbers. It adopts extended synchronous compression transform combined with automatic ridge extraction technology and uses ridge fusion standard to handle the time-varying component number problem, which can overcome noise interference. However, this method lacks a dedicated ridge repair mechanism for the time-frequency crossover region. When processing the crossover frequency signal of dual rotors of aero-engines, its ability to solve the problem of ridge breakage in the crossover region is limited.

[0010] Chinese patent document CN120336835A discloses a time-frequency analysis method for compressor dynamic pressure signals. By constructing a dual instantaneous frequency estimation operator and iteratively performing multiple local maximum synchronous compression transformations, the time-frequency energy concentration is improved and energy leakage is eliminated. However, this method does not specifically address the cross-region and lacks precise positioning and local optimization mechanisms for ridge fracture boundaries. It still has shortcomings in the multi-component separation and reconstruction of strong non-stationary cross-frequency signals. Summary of the Invention

[0011] The purpose of this invention is to provide a time-frequency analysis method and system for strong non-stationary cross-frequency signals of aero-engines. While maintaining the global energy concentration advantage of EMSSCT, it overcomes the limitations of traditional methods in terms of energy aliasing, ridge line breakage, and parameter sensitivity in the cross-frequency region. It achieves high-precision characterization and reconstruction of strong non-stationary cross-frequency signals, significantly improves the accuracy and stability of fault feature extraction under harmonic cross-frequency conditions of multi-rotor systems of aero-engines, and provides reliable data support for early fault diagnosis and in-service health monitoring.

[0012] To achieve the above objectives, the present invention provides the following technical solution: A time-frequency analysis method for strong non-stationary cross-frequency signals from aero-engines includes the following steps: S1: Acquire the vibration signal of the aero-engine, and perform linear frequency modulated wavelet transform on the vibration signal to obtain the initial time-frequency distribution; S2: Adaptive ridge extraction is performed on the initial time-frequency distribution to obtain an initial ridge set; for ridges with breaks in the initial ridge set, the break boundaries are located by bidirectional retrieval to determine the time-frequency block range of the intersection region; S3: Within the time-frequency block range of the intersection region, the time-frequency distribution under different transformation parameters is locally re-blocked, and an ideal time-frequency block is selected from the candidate time-frequency blocks based on the preset evaluation index to replace the time-frequency block at the corresponding position in the original time-frequency distribution, thereby obtaining the optimized time-frequency distribution; S4: Extract the complete ridge line based on the optimized time-frequency distribution, and reconstruct the vibration signal according to the instantaneous frequency information of the complete ridge line to obtain each component signal; S5: Perform local synchronous compression transformation on each component signal, superimpose the transformation results of each component to obtain the final time-frequency distribution, and identify the abnormal state of the aero-engine based on the final time-frequency distribution.

[0013] Further: The step of adaptive ridge extraction of the initial time-frequency distribution in S2 includes: S201: Identify local maxima on each time slice of the initial time-frequency distribution to form a candidate extremum matrix; S202: Perform envelope estimation on the candidate extreme value matrix to suppress noise spurious peaks; S203: A recursive strategy is used to correlate candidate peaks over time to track and form the initial ridge set.

[0014] Furthermore: In step S203, which employs a recursive strategy to correlate candidate peaks over time, a scoring function with a penalty term is used to select matching points. The expression for the scoring function is:

[0015] In the formula, At any moment Target frequency at the previous moment The optimal candidate frequency set obtained To predict the step size based on the frequency obtained from the target motion model, It is a weighting factor. It represents the target transition probability, and m is the penalty order. The linear frequency modulation transformation result in time index ,frequency The time-frequency coefficient at that location.

[0016] Furthermore: the steps in S2 for locating the fracture boundary through bidirectional retrieval include: S211: Construct an elimination function to remove the time-frequency energy corresponding to the complete components in the initial ridge set, and obtain the time-frequency distribution of the residual signal; S212: Perform ridge tracking on the time-frequency distribution of the residual signal along the positive time axis, and record the first time point when it cannot be extended as the left endpoint; S213: Perform ridge tracking on the time-frequency distribution of the residual signal in reverse along the time axis, and record the time point at which the first extension is impossible as the right endpoint; S214: Define the time-frequency block boundary where the ridge line is missing based on the left endpoint and the right endpoint.

[0017] Furthermore: the elimination function The expression is:

[0018] In the formula, The instantaneous frequency trajectory of the target. The neighborhood width, For frequency.

[0019] Further: The preset evaluation index is a linear weighted sum of the Rényi entropy and the reciprocal of the maximum amplitude. The step in S3 of selecting the ideal time-frequency block from the candidate time-frequency blocks based on the preset evaluation index includes: S31: Calculate the Rényi entropy and maximum instantaneous amplitude for each candidate time-frequency block; S32: Perform a reciprocal transformation on the instantaneous maximum amplitude to obtain the amplitude evaluation value; S33: Construct a comprehensive evaluation index, the expression of which is:

[0020] in For the first A comprehensive index of candidate time-frequency blocks, For the first Rényi entropy values ​​of candidate time-frequency blocks For the first Amplitude evaluation quantity, These are the weighting parameters for the Rényi entropy value and the magnitude evaluation value, respectively. The candidate time-frequency block with the smallest comprehensive evaluation index is selected as the ideal time-frequency block.

[0021] Furthermore: In step S4, the vibration signal is reconstructed based on the instantaneous frequency information of the complete ridge line. The Volk-Kalman order filtering method is used to separate the component signals by applying differential smoothing constraints on the amplitude trajectory and minimizing the reconstruction error.

[0022] Furthermore: In the Volk-Kalman order filtering method, the amplitude change is modeled as a low-order polynomial function, a structural equation in the state space is constructed, and the optimal amplitude estimate is solved by least squares optimization with a regularization term. The cost function expression of the least squares optimization is:

[0023] In the formula, It is a weighting factor used to balance the reconstruction error term. With smoothing regularization term This adjusts the equivalent bandwidth of VKF. This is the reconstruction error vector between the observed signal and the carrier model output. This is the magnitude deviation vector under the differential smoothing constraint; Let be the vector of component magnitudes to be estimated. The optimal solution is found here. The carrier matrix is ​​constructed from the instantaneous frequency of this component. For the observation matrix, For the observed signal vector, This represents the transpose operator.

[0024] Further: In step S5, where each component signal undergoes local synchronous compression transformation, the energy of each component signal is finely aggregated in the frequency direction to obtain a component-level time-frequency representation. The formula for this component-level time-frequency representation is:

[0025] in, This refers to the high-precision time-domain components output by VKF. This is a local maximum synchronous compression transform. This is a component-level time-frequency representation.

[0026] A time-frequency analysis system for strongly non-stationary cross-frequency signals of aero-engines based on any of the methods described above, characterized in that it includes: The initial time-frequency distribution module is used to acquire the vibration signal of the aero-engine, and perform linear frequency modulated wavelet transform on the vibration signal to obtain the initial time-frequency distribution. The fracture localization module is used to adaptively extract ridge lines from the initial time-frequency distribution to obtain an initial ridge line set; for ridge lines in the initial ridge line set that have fractures, the fracture boundary is located by bidirectional retrieval to determine the time-frequency block range of the intersection region; The local optimization module is used to locally subdivide the time-frequency distribution under different transformation parameters within the time-frequency block range of the intersection region, select the ideal time-frequency block from the candidate time-frequency blocks based on the preset evaluation index, and replace the time-frequency block at the corresponding position in the original time-frequency distribution to obtain the optimized time-frequency distribution. The component reconstruction module is used to extract the complete ridge line based on the optimized time-frequency distribution, and to reconstruct the vibration signal into components based on the instantaneous frequency information of the complete ridge line to obtain each component signal. An anomaly identification module is used to perform local synchronous compression transformation on each component signal, superimpose the transformation results of each component to obtain the final time-frequency distribution, and identify the abnormal state of the aero-engine based on the final time-frequency distribution.

[0027] Compared with the prior art, the present invention has the following advantages: I. This invention addresses the technical bottlenecks of traditional EMSSCT methods in handling strong non-stationary cross signals, such as ambiguous energy distribution and broken ridge lines in the cross region. It proposes a cross region localization and local representation optimization framework based on an improved adaptive ridge line extraction algorithm. By accurately locking the cross block boundary through bidirectional forward and reverse ridge line tracking, and introducing a novel TF block evaluation index jointly constructed from Rényi entropy and the inverse of amplitude, robust selection of ideal sub-blocks and energy consistency optimization are achieved. This effectively improves the time-frequency resolution and instantaneous frequency continuity of the cross region, solves the problem of inaccurate instantaneous frequency estimation caused by broken ridge lines in the cross region in existing methods, improves ridge line continuity and instantaneous frequency extraction accuracy, and significantly enhances the ability to fully represent the cross frequency trajectory.

[0028] II. This invention deeply integrates the Improved Entropy Matched Chirplet Transform (IEMCT) with Vold-Kalman Order Filtering (VKF) and Local Maximum Synchronous Compression Transform (LMSST) to form a complete IEMSSCT analysis framework. Based on the locally optimized global time-frequency representation, VKF is used to achieve high-precision amplitude and phase estimation and separation of multi-component signals. Then, LMSST is used to complete energy re-aggregation and fine reconstruction. This not only ensures the energy concentration and feature integrity of strongly non-stationary crossover signals, but also maintains high computational efficiency, breaking through the technical limitation of existing methods that make it difficult to balance separation accuracy and real-time performance in multi-component crossover scenarios. Attached Figure Description

[0029] Figure 1 A flowchart of a time-frequency analysis method for strong non-stationary cross-frequency signals of aero-engines provided by the present invention; Figure 2 A schematic diagram of the structure of a time-frequency analysis system for strong non-stationary cross-frequency signals of aero-engines provided by the present invention; Figure 3 (a) in the figure is the time-domain waveform of the experimental signal acquired by the dual-rotor fault simulation test bench for aero-engines; Figure 3 (b) in the figure is the spectrum of the experimental signals acquired by the dual rotor fault simulation test bench for aero-engines; Figure 3 (c) in the figure is a speed trend graph of the experimental signals collected by the dual rotor fault simulation test bench for aero-engines; Figure 3 (d) in the figure is the STFT result diagram of the experimental signals acquired by the aero-engine dual-rotor fault simulation test bench; Figure 4 (a) in the figure is the local optimization result of EMCT; Figure 4 (b) in the figure is the local optimization result of IEMCT; Figure 5 (a) in the image is a magnified view of the result corresponding to the IEMSSCT method; Figure 5(b) in the image is a magnified view of the intersection region corresponding to the EMSSCT method; Figure 5 (c) in the image is a magnified view of the intersection region corresponding to the ALCT method; Figure 5 (d) in the image is a magnified view of the intersection region corresponding to the ALCSET method; Figure 5 (e) in the image is a magnified view of the intersection region corresponding to the IEMSSCT method; Figure 5 (f) in the image is a magnified view of the non-intersecting region corresponding to the EMSSCT method; Figure 5 (g) in the image is a magnified view of the non-intersecting region corresponding to the ALCT method; Figure 5 (h) in the image is a magnified view of the non-intersecting region corresponding to the ALCSET method; Figure 5 In the image, (i) is a magnified view of the non-intersecting region corresponding to the IEMSSCT method; Figure 6 (a) in the figure is the IF trajectory extraction result of the EMSSCT method; Figure 6 (b) in the figure is the result of IF trajectory extraction using the ALCT method; Figure 6 (c) in the figure is the result of IF trajectory extraction by the ALCSET method; Figure 6 (d) in the figure is the IF trajectory extraction result of the IEMSSCT method; Figure 7 (a) in the figure is the reconstruction result of the EMSSCT method; Figure 7 (b) in the figure is the reconstruction result of the ALCT method; Figure 7 (c) in the figure is the reconstruction result of the ALCSET method; Figure 7 (d) in the figure is the reconstruction result of the IEMSSCT method. Detailed Implementation

[0030] The technical solution of the present invention will now be clearly and completely described with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0031] Example 1 like Figure 1 As shown, this invention provides a time-frequency analysis method for strong non-stationary crossover frequency signals in aero-engines, comprising the following steps: S1: Acquire the vibration signal of the aero-engine, and perform linear frequency modulated wavelet transform on the vibration signal to obtain the initial time-frequency distribution; S2: Adaptive ridge extraction is performed on the initial time-frequency distribution to obtain an initial ridge set; for ridges with breaks in the initial ridge set, the break boundaries are located by bidirectional retrieval to determine the time-frequency block range of the intersection region; S3: Within the time-frequency block range of the intersection region, the time-frequency distribution under different transformation parameters is locally re-blocked, and an ideal time-frequency block is selected from the candidate time-frequency blocks based on the preset evaluation index to replace the time-frequency block at the corresponding position in the original time-frequency distribution, thereby obtaining the optimized time-frequency distribution; S4: Extract the complete ridge line based on the optimized time-frequency distribution, and reconstruct the vibration signal according to the instantaneous frequency information of the complete ridge line to obtain each component signal; S5: Perform local synchronous compression transformation on each component signal, superimpose the transformation results of each component to obtain the final time-frequency distribution, and identify the abnormal state of the aero-engine based on the final time-frequency distribution.

[0032] According to a specific implementation of this embodiment, the step of adaptive ridge extraction of the initial time-frequency distribution in S2 includes: S201: Identify local maxima on each time slice of the initial time-frequency distribution to form a candidate extremum matrix; S202: Perform envelope estimation on the candidate extreme value matrix to suppress noise spurious peaks; S203: A recursive strategy is used to correlate candidate peaks over time to track and form the initial ridge set.

[0033] According to a specific implementation of this embodiment, in step S203, which uses a recursive strategy to perform data association on candidate peaks in the time dimension, a scoring function with a penalty term is used to select matching points. The expression of the scoring function is:

[0034] In the formula, At any moment Target frequency at the previous moment The optimal candidate frequency set obtained To predict the step size based on the frequency obtained from the target motion model, It is a weighting factor. It represents the target transition probability, and m is the penalty order. The linear frequency modulation transformation result in time index ,frequency The time-frequency coefficient at that location.

[0035] According to a specific implementation of this embodiment, the step of locating the fracture boundary by bidirectional retrieval in S2 includes: S211: Construct an elimination function to remove the time-frequency energy corresponding to the complete components in the initial ridge set, and obtain the time-frequency distribution of the residual signal; S212: Perform ridge tracking on the time-frequency distribution of the residual signal along the positive time axis, and record the first time point when it cannot be extended as the left endpoint; S213: Perform ridge tracking on the time-frequency distribution of the residual signal in reverse along the time axis, and record the time point at which the first extension is impossible as the right endpoint; S214: Define the time-frequency block boundary where the ridge line is missing based on the left endpoint and the right endpoint.

[0036] According to a specific implementation of this embodiment, the elimination function The expression is:

[0037] In the formula, The instantaneous frequency trajectory of the target. The width of the neighborhood.

[0038] According to a specific implementation of this embodiment, the preset evaluation index is a linear weighted sum of the Rényi entropy value and the reciprocal of the maximum amplitude value. The step in S3 of selecting the ideal time-frequency block from the candidate time-frequency blocks based on the preset evaluation index includes: S31: Calculate the Rényi entropy and maximum instantaneous amplitude for each candidate time-frequency block; S32: Perform a reciprocal transformation on the instantaneous maximum amplitude to obtain the amplitude evaluation value; S33: Construct a comprehensive evaluation index, the expression of which is:

[0039] in For the first A comprehensive index of candidate time-frequency blocks, For the first Rényi entropy values ​​of candidate time-frequency blocks For the first Amplitude evaluation quantity, These are the weighting parameters for the Rényi entropy value and the magnitude evaluation value, respectively. The candidate time-frequency block with the smallest comprehensive evaluation index is selected as the ideal time-frequency block.

[0040] According to a specific implementation of this embodiment, in the step of reconstructing the vibration signal based on the instantaneous frequency information of the complete ridge line in S4, the Volk-Kalman order filtering method is adopted. By applying differential smoothing constraints on the amplitude trajectory and minimizing the reconstruction error, the separation of each component signal is achieved.

[0041] According to a specific implementation of this embodiment, in the Vold-Kalman order filtering method, the amplitude change is modeled as a low-order polynomial function, a structure equation in the state space is constructed, and the optimal amplitude estimate is solved by least squares optimization with a regularization term. The cost function expression of the least squares optimization is as follows:

[0042] In the formula, It is a weighting factor used to balance the reconstruction error term. With smoothing regularization term This adjusts the equivalent bandwidth of VKF. This is the reconstruction error vector between the observed signal and the carrier model output. This is the magnitude deviation vector under the differential smoothing constraint; Let be the vector of component magnitudes to be estimated. The optimal solution is found here. The carrier matrix is ​​constructed from the instantaneous frequency of this component. For the observation matrix, For the observed signal vector, This represents the transpose operator.

[0043] According to a specific implementation of this embodiment, in step S5, where each component signal is subjected to local synchronous compression transformation, the energy of each component signal is finely aggregated in the frequency direction to obtain a component-level time-frequency representation. The formula for the component-level time-frequency representation is as follows:

[0044] in, This refers to the high-precision time-domain components output by VKF. This is a local maximum synchronous compression transform. This is a component-level time-frequency representation.

[0045] Example 2 like Figure 2 As shown, the present invention also provides a time-frequency analysis system for strongly non-stationary cross-frequency signals of aero-engines, comprising: The initial time-frequency distribution module is used to acquire the vibration signal of the aero-engine, perform linear frequency modulated wavelet transform on the vibration signal, and obtain the initial time-frequency distribution. The fracture localization module is used to adaptively extract ridges from the initial time-frequency distribution to obtain an initial ridge set. For ridges with fractures in the initial ridge set, the fracture boundary is located through bidirectional retrieval to determine the time-frequency block range of the intersection area. The local optimization module is used to locally subdivide the time-frequency distribution under different transformation parameters within the time-frequency block range of the intersection region. Based on the preset evaluation index, it selects the ideal time-frequency block from the candidate time-frequency blocks and replaces the time-frequency block at the corresponding position in the original time-frequency distribution to obtain the optimized time-frequency distribution. The component reconstruction module is used to extract the complete ridge line based on the optimized time-frequency distribution, and to reconstruct the vibration signal into components based on the instantaneous frequency information of the complete ridge line to obtain each component signal. The anomaly identification module performs local synchronous compression transformation on each component signal, superimposes the transformation results of each component to obtain the final time-frequency distribution, and identifies the abnormal state of the aero-engine based on the final time-frequency distribution. The functions and method steps of each module of the device correspond to each other, realizing high-precision time-frequency analysis and anomaly identification of the strong non-stationary cross-frequency signal of the aero-engine.

[0046] Example 3 like Figures 3-7 As shown in the figure, this embodiment provides a method for analyzing strongly non-stationary cross-frequency signals of aero-engines based on the Improved Entropy Matched Synchronous Compressed Chirplet Transform (IEMSSCT). The specific steps are as follows: In step S1, it is first necessary to monitor the basic information of the equipment operation. Specifically, considering the operating speed and physical characteristics of the equipment itself, the layout and linear monitoring range of the vibration sensors for monitoring are determined. During signal acquisition, the setting of the sampling frequency is crucial. According to the sampling theorem, to avoid frequency aliasing, the sampling frequency should be set to at least 2.56 times higher than the highest detection frequency. Specifically, the sampling rate of the vibration signal is set to no less than 20000Hz, and the maximum measurement range of vibration acceleration is 500g.

[0047] In step S2, local maxima are identified in each time slice to form a candidate extremum matrix. The envelope representation is then obtained by using a sequential statistical filter to suppress noise spurious peaks.

[0048] In the formula, Therefore The length of the sliding window centered on it. It is usually determined by the minimum distance between adjacent maxima:

[0049] The aforementioned envelope and candidate peak set serve as the observation input for the tracking phase. During the tracking process, the method employs a recursive strategy of "prior target → candidate matching": for Target frequency at time The candidate set at the current time The scoring function with a penalty term is used to select the most likely match, thereby achieving data association in the time dimension:

[0050] In the formula, It is a weighting factor. It is the target transition probability. It is the penalty order. Based on the probability density assumption, the target transition follows a Gaussian distribution, i.e.

[0051] In the formula, To represent a mean is The standard deviation is The Gaussian probability density function. It is the derivative of the time-frequency ridge at the target.

[0052] Finally, the mathematical expression for the TF ridge line extracted by this method is as follows:

[0053] However, the aforementioned improved adaptive time-frequency ridge extraction algorithm is based on the premise that the initial EMCT result itself is relatively complete. To address potential breakage, this invention first constructs a function to remove the complete components obtained from the initial EMCT result, thereby acquiring the residual signal. The function constructor is as follows:

[0054] In the formula, The neighborhood width is used to control the range of elimination in the frequency direction. A further endpoint localization strategy based on sequential and reverse time-axis retrieval is proposed: For the residual signal, sequential ridge tracking is performed along the time axis, and the time point at which further extension fails is recorded as the candidate left endpoint for that fragment; similarly, reverse tracking is performed on the residual signal, and the time point at which further extension fails is recorded as the candidate right endpoint. This two-way retrieval aims to accurately define the boundaries of TF blocks with missing ridges by identifying "forward interruptions" and "reverse interruptions" respectively, providing a clear domain for subsequent local refinement and energy redistribution.

[0055] After determining the ridge line missing boundary, the intersection region is confined to several small TF blocks. The sequential and reverse retrieval based on the time axis can effectively locate the left and right endpoints of the region, providing precise triggering conditions for subsequent implementation of encrypted time-frequency grid and local redistribution.

[0056] In step S3, the local representation optimization strategy operates within the located cross-window, aiming to automatically identify and replace the "ideal TF block" from the candidate sub-regions to obtain a local time-frequency representation with an energy distribution closer to the true value and continuous ridge lines. The overall process can be divided into four stages: mapping and extraction, local re-blocking, comprehensive evaluation and selection, and local replacement and splicing.

[0057] Mapping and candidate sub-block extraction: coordinates of the intersection point , Mapping to discrete TF block indices. Let the time starting point and frequency starting point of the global time-frequency representation be respectively... , The block lengths in the time and frequency directions are respectively , The mapping relationship is as follows:

[0058] Therefore, an index set is constructed:

[0059] Extract sets from the EMCT time-frequency plots corresponding to different discrete CR parameters. All corresponding candidate sub-regions will be used as the objects of subsequent evaluation.

[0060] Local subdivision: Following the approach of EMSSCT, each candidate sub-region is processed according to an encryption factor. , Local re-blocking is performed to improve the discriminability of nearest neighbor IFs. After block division, a set of candidate cells is obtained, and the minimum Rényi entropy value of each candidate cell block can be calculated and selected. As the final result, however, selecting the TF block solely based on the minimum Rényi entropy can lead to a misjudgment of "minimum entropy but severely insufficient energy"; therefore, it is necessary to introduce other information as a basis for judgment.

[0061] Construction and numericalization of comprehensive evaluation indicators: In constructing the comprehensive evaluation index, this invention employs a linear weighted form of two types of original quantities: Rényi entropy and the inverse of amplitude. For the candidate sub-block set within the cross-domain... Let the Rényi entropy of each candidate sub-block be:

[0062] And record the maximum instantaneous amplitude as:

[0063] To unify the use of amplitude and entropy information for performance evaluation, a reciprocal transformation is applied to the amplitude in a monotonic direction that maintains the principle of "the smaller the better":

[0064] in It is a very small constant to ensure the numerical stability of the calculation process.

[0065] Based on this, a linearly weighted composite index is defined. :

[0066] in These are the weighting parameters for the Rényi entropy and the reciprocal of the maximum amplitude, respectively. The final "ideal TF block" is determined by the following formula:

[0067] Replacement of ideal TF blocks and global representation update: After selecting the ideal sub-block within each intersection region, its corresponding local TFR data is replaced in the corresponding position of the original global TFR to obtain a new global time-frequency representation TFR3. After the replacement is completed, standard ridge extraction operations are performed on TFR3 to obtain a more complete and continuous ridge trajectory IFk, thus providing a stable input for subsequent high-precision IF estimation and component reconstruction.

[0068] In step S4, VKF, by applying differential smoothing constraints on the amplitude trajectory and aiming to minimize reconstruction error, can achieve robust separation of harmonic components of any time-varying IF while maintaining phase consistency. It is particularly suitable for energy distribution and amplitude recovery in cases of crossed or near-near IFs. Based on Fourier transform, the vibration signal measured in practical engineering can be expressed as:

[0069] In the formula, For the instantaneous phase of the signal, This represents the error term containing noise. For mechanical vibration signals acquired in actual engineering projects, their amplitude changes are usually relatively smooth, so low-order polynomial functions can be used to model them. Based on this, the VKF process model can be defined as:

[0070] In the formula, Represents the difference operator. Indicates the difference order, representing the order of the magnitude sequence. Order difference operator, additional terms This represents process noise, indicating non-target signals or interference components. Assume the difference order... Considering all samples, the structural equations in state space can be constructed in matrix form as follows:

[0071] Therefore, the above formula can be expressed as:

[0072] In the formula, This represents the carrier matrix. The carrier matrix of each component can be based on its instantaneous frequency. Perform the calculation.

[0073] To solve the above equation, the problem can be transformed into a least-squares optimization problem with a regularization term. By constructing a cost function and finding its minimum value, the optimal amplitude estimate can be obtained, thereby achieving signal component reconstruction. The least-squares calculation formula is as follows:

[0074] In the formula, It is a weighting factor used to balance and To control VKF bandwidth.

[0075] Therefore, the amplitude matrix of the component with the largest amplitude. It can be calculated as follows:

[0076] Based on the above steps, the component signal with the largest amplitude can be calculated as follows:

[0077] Based on the aforementioned analysis of the VKF method, combined with the extracted instantaneous frequency ridges... This information can be further used to separate and reconstruct the original signal with high precision. The reconstructed signal can accurately recover its amplitude variation trend while maintaining the integrity of the original signal's phase information.

[0078] To obtain high-precision time-domain components of VKF output Subsequently, the present invention performs LMSST on each component to finely focus the energy of the component in the frequency direction, thereby obtaining a component-level time-frequency representation.

[0079] After completing all the above steps, the LMSST characterization results of each component are superimposed on the same time-frequency plane to construct the final global time-frequency representation of IEMSSCT.

[0080] To better illustrate the implementation effect of this invention, a dual-rotor fault simulation test bench for an aero-engine is used as an example. This test bench mainly consists of two frequency converters, an electric spindle for the outer rotor, an outer rotor coupling, bearing housings, the outer rotor, a No. 1 intermediate bearing, a No. 2 bearing housing (containing a bearing casing), bearing housings, the inner rotor, an inner rotor coupling, and an electric spindle for the inner rotor. This test bench can simulate the operating conditions of a dual-rotor aero-engine, mainly by adjusting the frequencies of the two frequency converters to change the speeds of the inner and outer rotors, thereby simulating the operation of the high- and low-pressure rotors of an aero-engine.

[0081] During the experiment, two eddy current sensors were used to acquire the key phase signals of the dual rotors, obtaining the actual rotational speeds of the high- and low-pressure rotors of the aero-engine. Two eddy current sensors (Bently 3300 XL, 3.94 mV / μm) were used to acquire the rotational speed signals of the two rotors, and an accelerometer (BK4511, 1 mV / (m / s²)) was used to collect the rotor vibration signals. The data acquisition equipment was a 16-channel LMS SCADAS, specifically designed for vibration and acoustic testing.

[0082] The experimental signal was sampled at a frequency of 51200 Hz for 10 seconds. The invention was used to process acceleration vibration signals, and the effectiveness was verified using key phase signals.

[0083] like Figure 3 As shown, (a) is the time-domain waveform of the rotor harmonic vibration signal; (b) is the spectrum analysis obtained by applying a fast Fourier transform to the signal; (c) is the speed variation trend of the high-voltage and low-voltage rotors, where N1 represents the rotational frequency of the high-voltage rotor and N2 represents the rotational frequency of the low-voltage rotor.

[0084] To analyze the vibration signal, the original signal was downsampled to 64Hz. The time-frequency result generated by the Short Time Fourier Transform (STFT) is as follows: Figure 3 As shown in (d), the acquired signal after the test setup is completed contains three main components: the rotational frequency S1 of the high-voltage rotor, the rotational frequency S2 of the low-voltage rotor, and its second harmonic frequency S3. Furthermore, there is a time-frequency overlap between the high-voltage rotor's rotational frequency S1 and the second harmonic component S3 of the low-voltage rotor's rotational frequency. Therefore, the vibration signal obtained from the experiment is a signal with strong time-varying characteristics and frequency overlap.

[0085] First, we analyze the local optimization performance of this invention. For example... Figure 4 As shown, by locating the intersection region and performing secondary block division, the improved EMCT (IEMCT) results obtained have a time-frequency characterization of energy concentration at the local intersection position, providing a good data foundation for subsequent signal separation and reconstruction.

[0086] The time-frequency results generated by IEMSSCT are as follows Figure 5 As shown in (a), IEMSSCT clearly characterizes the main components of the experimental signal, and the time-frequency crossover phenomenon between components S1 and S3 is also clearly expressed. Figure 3 Compared to the STFT results in (d), the IEMSSCT results show clearer signal characteristics, which helps to accurately extract signal features.

[0087] To verify the effectiveness of the proposed invention, time-frequency representations obtained by the EMSSCT, ALCT, and ALCSET methods are demonstrated, such as... Figure 5 As shown in (bi), the above methods can effectively represent the time-frequency characteristics of the S1 component in the cross-region of the experimental signal. Although the above methods can redistribute the diffuse energy around the instantaneous frequency trajectory of the S1 component, the TFR of the S2 and S3 components of the signal cannot accurately characterize the instantaneous features of each component. In contrast, IEMSSCT provides a complete and clear time-frequency representation.

[0088] The IF trajectories of the main components of the signal extracted from the results of EMSSCT, ALCT, ALCSET, and IEMSSCT are as follows: Figure 6 As shown, EMSSCT, ALCT, and ALCSET all exhibit ridge line breakage, resulting in incomplete extraction. EMSSCT, however, extracts TF ridges that are complete and highly consistent with the actual IF trajectory. This greatly facilitates subsequent analysis of the principal components, demonstrating that under real-world conditions, the ridge line extraction in this invention is complete and highly accurate.

[0089] Signal reconstruction is further achieved by utilizing the IF trajectory extracted by IEMSSCT. Figure 7 It can be seen that the reconstructed signal retains the amplitude and frequency characteristics of the original signal, fully demonstrating the superior ability of IEMSSCT in time-frequency crossover signal reconstruction.

[0090] By introducing key technologies such as improved adaptive ridge extraction, local characterization optimization, and component-level reconstruction, high-resolution characterization and accurate separation of vibration features in intersection regions are achieved. This method fully leverages the advantages of locally refined time-frequency grids and comprehensive evaluation indices, significantly improving ridge continuity and instantaneous frequency extraction accuracy in intersection regions while ensuring global energy concentration. Through the collaborative work of IEMCT, Volk-Kalman order filtering (VKF), and Local Maximum Synchronous Compression Transform (LMSST), a complete multi-component high-precision time-frequency analysis framework is constructed, enabling anomaly identification and early fault diagnosis of strong non-stationary signals under complex aero-engine operating conditions, providing reliable data support and engineering application assurance for flight safety.

[0091] The above embodiments are only for illustrating the technical concept and features of the present invention, and are intended to enable those skilled in the art to understand the content of the present invention and implement it accordingly. They should not be construed as limiting the scope of protection of the present invention. All equivalent transformations or modifications made in accordance with the spirit and essence of the present invention should be covered within the scope of protection of the present invention.

Claims

1. A time-frequency analysis method for strongly non-stationary cross-frequency signals of aero-engines, characterized in that: Includes the following steps: S1: Acquire the vibration signal of the aero-engine, and perform linear frequency modulated wavelet transform on the vibration signal to obtain the initial time-frequency distribution; S2: Adaptive ridge extraction is performed on the initial time-frequency distribution to obtain an initial ridge set; for ridges with breaks in the initial ridge set, the break boundaries are located by bidirectional retrieval to determine the time-frequency block range of the intersection region; S3: Within the time-frequency block range of the intersection region, the time-frequency distribution under different transformation parameters is locally re-blocked, and an ideal time-frequency block is selected from the candidate time-frequency blocks based on the preset evaluation index to replace the time-frequency block at the corresponding position in the original time-frequency distribution, thereby obtaining the optimized time-frequency distribution; S4: Extract the complete ridge line based on the optimized time-frequency distribution, and reconstruct the vibration signal according to the instantaneous frequency information of the complete ridge line to obtain each component signal; S5: Perform local synchronous compression transformation on each component signal, superimpose the transformation results of each component to obtain the final time-frequency distribution, and identify the abnormal state of the aero-engine based on the final time-frequency distribution.

2. The time-frequency analysis method for aero-engine strongly non-stationary cross-frequency signals according to claim 1, characterized in that: The step of adaptive ridge extraction of the initial time-frequency distribution in S2 includes: S201: Identify local maxima on each time slice of the initial time-frequency distribution to form a candidate extremum matrix; S202: Perform envelope estimation on the candidate extreme value matrix to suppress noise spurious peaks; S203: A recursive strategy is used to correlate candidate peaks over time to track and form the initial ridge set.

3. The time-frequency analysis method for aero-engine strongly non-stationary cross-frequency signals according to claim 2, characterized in that: In step S203, which employs a recursive strategy to correlate candidate peaks over time, a scoring function with a penalty term is used to select matching points. The expression for the scoring function is: In the formula, At any moment Target frequency at the previous moment The optimal candidate frequency set obtained To predict the step size based on the frequency obtained from the target motion model, It is a weighting factor. It represents the target transition probability, and m is the penalty order. The linear frequency modulation transformation result in time index ,frequency The time-frequency coefficient at that location.

4. The time-frequency analysis method for aero-engine strongly non-stationary cross-frequency signals according to claim 1, characterized in that: The steps for locating the fracture boundary through bidirectional retrieval in S2 include: S211: Construct an elimination function to remove the time-frequency energy corresponding to the complete components in the initial ridge set, and obtain the time-frequency distribution of the residual signal; S212: Perform ridge tracking on the time-frequency distribution of the residual signal along the positive time axis, and record the first time point when it cannot be extended as the left endpoint; S213: Perform ridge tracking on the time-frequency distribution of the residual signal in reverse along the time axis, and record the time point at which the first extension is impossible as the right endpoint; S214: Define the time-frequency block boundary where the ridge line is missing based on the left endpoint and the right endpoint.

5. The time-frequency analysis method for aero-engine strongly non-stationary cross-frequency signals according to claim 4, characterized in that: The elimination function The expression is: In the formula, The instantaneous frequency trajectory of the target. The neighborhood width, For frequency.

6. The time-frequency analysis method for aero-engine strongly non-stationary cross-frequency signals according to claim 1, characterized in that: The preset evaluation index is a linear weighted sum of the Rényi entropy and the reciprocal of the maximum amplitude. The step in S3 of selecting the ideal time-frequency block from the candidate time-frequency blocks based on the preset evaluation index includes: S31: Calculate the Rényi entropy and maximum instantaneous amplitude for each candidate time-frequency block; S32: Perform a reciprocal transformation on the instantaneous maximum amplitude to obtain the amplitude evaluation value; S33: Construct a comprehensive evaluation index, the expression of which is: in For the first A comprehensive index of candidate time-frequency blocks, For the first Rényi entropy values ​​of candidate time-frequency blocks For the first Amplitude evaluation quantity, These are the weighting parameters for the Rényi entropy value and the magnitude evaluation value, respectively. The candidate time-frequency block with the smallest comprehensive evaluation index is selected as the ideal time-frequency block.

7. The time-frequency analysis method for aero-engine strongly non-stationary cross-frequency signals according to claim 1, characterized in that: In step S4, which involves reconstructing the vibration signal based on the instantaneous frequency information of the complete ridge, the Volk-Kalman order filtering method is used. By applying differential smoothing constraints on the amplitude trajectory and minimizing the reconstruction error, the separation of each component signal is achieved.

8. The time-frequency analysis method for aero-engine strongly non-stationary cross-frequency signals according to claim 7, characterized in that: In the Volk-Kalman order filtering method, the amplitude change is modeled as a low-order polynomial function, constructing a structure equation in the state space. The optimal amplitude estimate is solved by least squares optimization with a regularization term. The cost function expression of the least squares optimization is as follows: In the formula, It is a weighting factor used to balance the reconstruction error term. With smoothing regularization term This adjusts the equivalent bandwidth of VKF. This is the reconstruction error vector between the observed signal and the carrier model output. This is the magnitude deviation vector under the differential smoothing constraint; Let be the vector of component magnitudes to be estimated. The optimal solution is found here. The carrier matrix is ​​constructed from the instantaneous frequency of this component. For the observation matrix, For the observed signal vector, This represents the transpose operator.

9. The time-frequency analysis method for aero-engine strongly non-stationary cross-frequency signals according to claim 1, characterized in that: In step S5, where each component signal undergoes local synchronous compression transformation, the energy of each component signal is finely aggregated in the frequency direction to obtain a component-level time-frequency representation. The formula for this component-level time-frequency representation is as follows: in, The time-domain component of the VKF output. This is a local maximum synchronous compression transform. This is a component-level time-frequency representation.

10. A time-frequency analysis system for strongly non-stationary cross-frequency signals of aero-engines using the method of any one of claims 1-9, characterized in that, include: The initial time-frequency distribution module is used to acquire the vibration signal of the aero-engine, and perform linear frequency modulated wavelet transform on the vibration signal to obtain the initial time-frequency distribution. The fracture localization module is used to adaptively extract ridge lines from the initial time-frequency distribution to obtain an initial ridge line set; for ridge lines in the initial ridge line set that have fractures, the fracture boundary is located by bidirectional retrieval to determine the time-frequency block range of the intersection region; The local optimization module is used to locally subdivide the time-frequency distribution under different transformation parameters within the time-frequency block range of the intersection region, select the ideal time-frequency block from the candidate time-frequency blocks based on the preset evaluation index, and replace the time-frequency block at the corresponding position in the original time-frequency distribution to obtain the optimized time-frequency distribution. The component reconstruction module is used to extract the complete ridge line based on the optimized time-frequency distribution, and to reconstruct the vibration signal into components based on the instantaneous frequency information of the complete ridge line to obtain each component signal. An anomaly identification module is used to perform local synchronous compression transformation on each component signal, superimpose the transformation results of each component to obtain the final time-frequency distribution, and identify the abnormal state of the aero-engine based on the final time-frequency distribution.

Images