Aero-engine performance detection method

Abstract

The present application relates to the technical field of engine testing, in particular to an aero-engine performance detection method, comprising the following steps: collecting vibration signals installed on the engine shell and blade tip timing data around the compressor blade periphery, segmenting the continuous vibration signals according to the engine speed cycle, calculating the time-frequency distribution of each signal segment, and obtaining the time-frequency energy matrix. In the present application, the vibration signals and the blade tip timing data are synchronously collected and segmented by combining the speed cycle, so that the frequency characteristics in different time periods can be completely extracted, avoiding the missing information caused by single time domain or frequency domain processing in the traditional method. Further, the spectral density is calculated in the frequency domain and the cyclic frequency domain by the time-frequency energy matrix, a two-dimensional spectrum reflecting the energy distribution and the cyclic characteristics is constructed, and the sufficient expression of the characteristic information in multiple dimensions is ensured.

Description

Technical Field

[0001] This invention relates to the field of engine testing technology, and in particular to a method for testing the performance of an aero-engine. Background Technology

[0002] As a critical power unit for aircraft, the performance of aero engines directly affects the operational safety of aircraft and the safety of passengers and crew. Therefore, it is of great significance to detect abnormalities in engine performance.

[0003] Current technologies rely on overall vibration levels or local parameter deviations for judgment, leading to blind spots in identifying characteristic behaviors under complex coupling conditions. Furthermore, relying solely on overall parameter thresholds is susceptible to misjudgments due to environmental noise and random disturbances. Additionally, they are not sensitive enough to abnormal responses in local blades or specific orders, making it difficult to detect early local faults in a timely manner. For example, a slight increase in the vibration amplitude at a certain order may not significantly exceed the overall vibration threshold, thus masking potential problems and delaying maintenance. Therefore, improvements are needed. Summary of the Invention

[0004] The purpose of this invention is to address the shortcomings of existing technologies by proposing a method for testing the performance of aero-engines.

[0005] To achieve the above objectives, the present invention adopts the following technical solution: a method for testing the performance of an aero-engine, comprising the following steps:

[0006] Vibration signals from the engine housing and timing data from the blade tips around the compressor blades are collected. The continuous vibration signals are segmented according to the engine speed cycle, and the time-frequency distribution of each signal segment is calculated to obtain the time-frequency energy matrix.

[0007] Based on the time-frequency energy matrix, the spectral density in the frequency domain and the cyclic frequency domain is calculated, a two-dimensional cyclic spectrum is established, and the energy peaks and distributions along the target cyclic frequency line are extracted from the two-dimensional cyclic spectrum to obtain the modulation feature vector.

[0008] Based on the blade tip timing data, the first five engine order vibration amplitudes and resonant peak frequencies of each compressor blade are extracted, and the combined parameters are used to construct a multi-order blade vibration metric. The modulation feature vector is fused with the multi-order blade vibration metric to obtain the integrated blade dynamic state point.

[0009] Based on the integrated blade dynamic state points, the deviation between the blade state and the preset healthy state mean vector is calculated to generate a blade state deviation vector. The engine dynamic anomaly score is then calculated using the blade state deviation vector and the healthy state covariance matrix.

[0010] Preferably, the step of obtaining the time-frequency energy matrix is ​​as follows:

[0011] Based on the acceleration channel installed on the engine casing and the blade tip timing channel around the compressor blades, the vibration signal timestamp sequence and vibration signal amplitude sequence are read, the blade tip timing timestamp sequence and blade tip timing identifier sequence are read, the time base is unified and the data is organized according to the acquisition order to generate vibration signal and blade tip timing data.

[0012] Based on the vibration signal and blade tip timing data, the blade tip timing timestamp sequence is analyzed to calculate the passing interval of adjacent blades, generating an engine speed cycle sequence. The engine speed cycle sequence is used to locate segmentation points on the vibration signal timestamp sequence and slice the vibration signal amplitude sequence to form a signal segment sequence.

[0013] Based on the signal segment sequence, a fixed time window length and time window overlap ratio are set to perform sliding statistics on each signal segment. The energy of each frequency grid within each time window is counted and the time index corresponding to the center of the time window is recorded. The energy values ​​are arranged according to the time index and frequency grid index to form a two-dimensional structure, generating a time-frequency energy matrix.

[0014] Preferably, the step of obtaining the two-dimensional cyclic spectrum is as follows:

[0015] Based on the time-frequency energy matrix, the energy is accumulated in the time index direction according to the frequency grid index and the amplitude is normalized. The energy periodicity is statistically analyzed in the frequency grid direction according to the time index and rearranged according to the cyclic frequency index to generate the frequency domain spectral density and the cyclic frequency domain spectral density.

[0016] Based on the frequency domain spectral density and the cyclic frequency domain spectral density, a two-dimensional coordinate grid with cyclic frequency index and frequency grid index is established. The energy density is filled according to the coordinate grid, and interpolation is performed at the boundaries and smoothing is performed at discrete points to generate a two-dimensional cyclic spectrum.

[0017] Preferably, the step of obtaining the modulation feature vector is as follows:

[0018] Based on the two-dimensional cyclic spectrum, a set of target cyclic frequency lines is constructed according to the engine frequency index and integer multiple index. Equal step sampling and extreme value search are performed along each target cyclic frequency line, and the peak position and peak amplitude distribution are statistically analyzed to generate a modulation feature vector.

[0019] Preferably, the steps for obtaining the multi-order blade vibration measurement are as follows:

[0020] Based on the blade tip timing data, the passing time series is organized by blade number, the time interval between adjacent passing is calculated and mapped to the order index, the vibration amplitude of each blade is extracted along the order index and sorted by amplitude, the resonance peak frequency of the corresponding order is located, and the vibration amplitude and resonance peak frequency of the first five engine orders are obtained.

[0021] Based on the vibration amplitude and resonance peak frequency of the first five engine orders, a unique index of blade number and order number is constructed, vibration amplitude is normalized and resonance peak frequency is standardized, and the order is spliced ​​into a fixed-length parameter sequence and blade number is added to generate multi-order blade vibration measurement.

[0022] Preferably, the step of obtaining the dynamic state point of the integrated blade is as follows:

[0023] Based on the multi-order blade vibration measurement, the corresponding parameters are located in the modulation feature vector according to the blade number and aligned according to the order index. The current working condition multi-dimensional feature points of a single blade are established and numbered and sorted to generate the dynamic state points of the integrated blade.

[0024] Preferably, the step of obtaining the blade state deviation vector is as follows:

[0025] Based on the integrated blade dynamic state points, align them one by one with the dimensions and indices of the preset health state mean vector, verify that the numerical types and dimensions are consistent, remove missing and duplicate records, subtract each component while keeping the original index order unchanged, and obtain the blade state deviation vector.

[0026] Preferably, the step of obtaining the engine dynamic anomaly score is as follows:

[0027] Based on the blade state deviation vector, the health state covariance matrix is ​​read and its symmetry and positive definiteness are verified. If not satisfied, it is only adjusted up by the minimum amplitude of the diagonal elements until it is positive definite. The inverse is calculated and the square root of the diagonal elements is extracted as the standard deviation sequence to obtain the inverse matrix of the health state covariance matrix and the standard deviation sequence.

[0028] The engine dynamic anomaly score is calculated based on the inverse of the health status covariance matrix and the standard deviation sequence.

[0029] Compared with the prior art, the advantages and positive effects of the present invention are as follows:

[0030] In this invention, by synchronously acquiring vibration signals and blade tip timing data and performing segmented processing in conjunction with rotational speed cycles, the frequency characteristics within different time periods can be completely extracted, avoiding information missed due to single time-domain or frequency-domain processing in traditional methods. Furthermore, the spectral density is calculated in both the frequency and cyclic frequency domains using a time-frequency energy matrix, constructing a two-dimensional spectrum that simultaneously reflects energy distribution and cyclic characteristics. This ensures the full expression of feature information across multiple dimensions. The modulation feature vector is obtained by extracting the energy peak and distribution of the target cyclic frequency line, thus capturing periodic vibrations and harmonic effects. Simultaneously, the blade tip timing... Data extraction of multi-order vibration amplitude and resonant frequency of each blade forms a multi-order blade vibration metric, which is then fused with the modulation feature vector. This allows the dynamic state of a single blade to be quantified at multi-dimensional feature points. Finally, by comparing the deviation between the integrated blade dynamic state points and the mean of the healthy state, and combining this with the healthy state covariance matrix for comprehensive calculation, the generated anomaly score can quantify the degree of overall performance deviation, enabling the identification of performance anomalies under complex coupled operating conditions. This not only improves the accuracy and robustness of anomaly detection but also enhances the response capability to deviations at different levels of features, allowing potential engine faults to be identified and quantified at an early stage. Attached Figure Description

[0031] Figure 1 This is a schematic diagram of the steps of the present invention. Detailed Implementation

[0032] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0033] Please see Figure 1 This invention provides a technical solution, a method for testing the performance of an aero-engine, comprising the following steps:

[0034] Vibration signals from the engine housing and timing data from the blade tips around the compressor blades are collected. The continuous vibration signals are segmented according to the engine speed cycle, and the time-frequency distribution of each signal segment is calculated to obtain the time-frequency energy matrix.

[0035] Based on the time-frequency energy matrix, the spectral density in the frequency domain and the cyclic frequency domain is calculated, a two-dimensional cyclic spectrum is established, and the energy peaks and distributions along the target cyclic frequency line are extracted from the two-dimensional cyclic spectrum to obtain the modulation feature vector.

[0036] Based on the tip timing data, the vibration amplitude and resonant peak frequency of the first five engine orders of each compressor blade are extracted, and the combined parameters are used to construct a multi-order blade vibration metric. The modulation feature vector is fused with the multi-order blade vibration metric to obtain the integrated blade dynamic state point.

[0037] Based on the integrated blade dynamic state points, the deviation between the blade state and the preset healthy state mean vector is calculated to generate a blade state deviation vector. The engine dynamic anomaly score is then calculated using the blade state deviation vector and the healthy state covariance matrix.

[0038] The steps for obtaining the time-frequency energy matrix are as follows:

[0039] Based on the acceleration channel installed on the engine casing and the blade tip timing channel around the compressor blades, the vibration signal timestamp sequence and vibration signal amplitude sequence are read, the blade tip timing timestamp sequence and blade tip timing identifier sequence are read, the time base is unified and the data is organized according to the acquisition order to generate vibration signal and blade tip timing data.

[0040] Based on the vibration signal and blade tip timing data, the blade tip timing timestamp sequence is analyzed to calculate the passing interval of adjacent blades, generating an engine speed cycle sequence. The engine speed cycle sequence is used to locate segmentation points on the vibration signal timestamp sequence and slice the vibration signal amplitude sequence to form a signal segment sequence.

[0041] Based on the signal segment sequence, a fixed time window length and time window overlap ratio are set to perform sliding statistics on each signal segment. The energy of each frequency grid within each time window is counted and the time index corresponding to the center of the time window is recorded. The energy values ​​are arranged according to the time index and frequency grid index to form a two-dimensional structure, generating a time-frequency energy matrix.

[0042] Specifically, based on the acceleration channel installed on the engine casing and the blade tip timing channel around the compressor blades, the time reference of the multi-source data is first aligned. This involves using a unified global positioning system clock signal to synchronize the acquisition devices of the acceleration sensor and the blade tip timing sensor, ensuring that the error of all timestamps is within 1 microsecond. Then, the timestamp sequence and corresponding amplitude sequence of the vibration signal, as well as the timestamp sequence and blade identification sequence of the blade tip timing, are read in parallel. During the reading process, preliminary data quality checks are performed, such as checking whether the vibration signal amplitude exceeds the sensor's range. A saturation threshold is set, which is determined according to the sensor specifications. For example, for an acceleration sensor with a range of ±50g, the saturation threshold can be... The sampling point is set to 49g. Any sampling point exceeding this value is marked as invalid. At the same time, the continuity of the blade tip timing timestamp sequence is checked. By comparing the interval between adjacent timestamps with the expected interval calculated from the previous rotational speed cycle, if the deviation exceeds the preset tolerance, which is set to 5% of the average rotational speed cycle, for example, at 6000 RPM, the cycle is 10 milliseconds and the tolerance is 0.5 milliseconds, then the blade tip timing data point is marked as suspected missing. The marked data is processed. For example, saturated vibration signal points are filled by linear interpolation of the nearest valid points, and the missing blade tip timing data points are estimated based on the linear extrapolation of the passage time of the same blade before and after. After alignment and verification, vibration signals and blade tip timing data containing synchronization timestamps are generated.

[0043] Based on the vibration signal and blade tip timing data, the passage time of each blade is extracted from the blade tip timing timestamp sequence. Specific blades are identified by the blade tip timing identifier sequence, and the difference between the times when the blade passes the same sensor twice consecutively is calculated. This difference is one speed cycle of the engine at that moment. To improve stability and accuracy, the speed cycles calculated for all blades within a complete rotation cycle are averaged to obtain the average speed cycle of the current rotation. The average speed cycles of all rotations are arranged in chronological order to form an engine speed cycle sequence. Then, this engine speed cycle sequence is used to segment the continuous vibration signal timestamp sequence. Using the start and end timestamps of each speed cycle as boundaries, data for the corresponding time period is extracted from the vibration signal amplitude sequence, thereby dividing the long-term vibration signal into a series of signal segments synchronized with the engine rotation cycle, forming a signal segment sequence.

[0044] Based on the signal segment sequence, a short-time Fourier transform is applied to each signal segment to analyze its time-varying frequency characteristics. First, a fixed time window length and overlap ratio are set. The time window length is chosen to balance the desired frequency resolution with the time resolution, and its value is set to 1 / 16 of the average number of points in the signal segment. For example, if a signal segment contains 4096 sampling points, the time window length is set to 256 sampling points. The overlap ratio of the time window is set to an empirical value, here set to 75%. A Hanning window is used as the window function. Then, this configured sliding time window is applied to each signal segment. The process involves sliding the signal segment from beginning to end within each window. For each captured signal segment, a Fast Fourier Transform is performed to calculate the complex values ​​of each frequency component. The square of the modulus is then taken to obtain the energy value, which is assigned to the corresponding frequency grid. Simultaneously, the time index corresponding to the center point of each time window is recorded. Finally, the energy values ​​calculated from all time windows are arranged according to the time index and frequency grid index to construct a two-dimensional data structure, where one dimension represents time and the other represents frequency. The value of each element is the energy density of that frequency at that moment, thus generating a time-frequency-energy matrix.

[0045] The steps for obtaining a two-dimensional cyclic spectrum are as follows:

[0046] Based on the time-frequency energy matrix, the energy is accumulated in the time index direction according to the frequency grid index and the amplitude is normalized. The energy periodicity is statistically analyzed in the frequency grid direction according to the time index and rearranged according to the cyclic frequency index to generate the frequency domain spectral density and the cyclic frequency domain spectral density.

[0047] Based on the frequency domain spectral density and the cyclic frequency domain spectral density, a two-dimensional coordinate grid with cyclic frequency index and frequency grid index is established. The energy density is filled according to the coordinate grid, and the boundary is interpolated and smoothed at discrete points to generate a two-dimensional cyclic spectrum.

[0048] Specifically, based on the time-frequency energy matrix, to calculate the frequency domain spectral density, we first sum all energy values ​​under each independent frequency grid index along the time index direction. Specifically, for each row in the matrix, we sum the values ​​of all columns in that row to obtain a scalar value representing the total energy of that frequency. After performing this operation on all rows, we obtain a one-dimensional array whose length is equal to the total number of frequency grids. Each element in the array corresponds to the total energy of one frequency grid. Then, we perform amplitude normalization on this one-dimensional array, calculate the maximum value of all elements in the array, and then divide each element in the array by this maximum value, scaling all energy values ​​to the range of 0 to 1, thus obtaining the frequency domain spectral density. Next, to calculate the cyclic frequency domain spectral density, we need to perform periodic analysis on the energy changes over time in the time-frequency energy matrix. For each frequency grid index, i.e., each row of the matrix, it is treated as a time series. We perform a Fourier transform on this time series to calculate its spectral components at different cyclic frequencies. The formula for this process is as follows: Where S(α,f) is the cyclic spectral density function, α is the cyclic frequency, f is the carrier frequency (corresponding to the frequency grid index), E(t,f) is the energy value at time t and frequency f in the time-frequency energy matrix, and T is the total observation time. In the discrete implementation, this is done by applying a Fast Fourier Transform (FFT) to each row of data. The result of the transformation is a complex matrix, and its amplitude part is the energy distribution at different cyclic frequencies. This result is rearranged according to the cyclic frequency index to generate the cyclic frequency domain spectral density.

[0049] Based on the frequency domain spectral density and the cyclic frequency domain spectral density, a two-dimensional coordinate grid is first initialized. The horizontal axis (X-axis) represents the frequency grid index, ranging from the lowest to the highest frequency, and the vertical axis (Y-axis) represents the cyclic frequency index, covering all calculated cyclic frequencies. Then, the energy density values ​​from the cyclic frequency domain spectral density matrix are sequentially filled into the corresponding coordinate points of this two-dimensional coordinate grid, forming a discrete energy distribution map. Due to the discrete nature of the calculation process, data may be missing or discontinuous in the boundary regions of the image. Therefore, boundary interpolation processing is required. A bilinear interpolation algorithm is used. For any coordinate point on the grid boundary without a value, the four nearest neighbors with existing values ​​are found. A weighted average is calculated based on the distance from the point to these four neighbors. The interpolation at this point is weighted inversely to the distance. To remove random noise introduced by discrete calculations and smooth the features in the spectrum, a two-dimensional Gaussian smoothing filter is applied to the entire filled and interpolated grid data. The Gaussian kernel size is set to 5x5 pixels, and the standard deviation σ is set to 1.5. The kernel size and standard deviation are empirically selected based on the spectrum resolution and noise level; for example, a smaller kernel can be used if the spectrum resolution is high, and a larger kernel if the resolution is low. The standard deviation controls the degree of smoothing; a larger value results in a stronger smoothing effect. This 5x5 Gaussian kernel is convolved with each point in the coordinate grid and its neighborhood. The calculated weighted average value replaces the original value of the center point. After boundary interpolation and global smoothing, a visualized two-dimensional cyclic spectrum is finally generated.

[0050] The steps for obtaining the modulation feature vector are as follows:

[0051] Based on the two-dimensional cyclic spectrum, a set of target cyclic frequency lines is constructed according to the engine frequency index and integer multiple index. Equal step sampling and extreme value search are performed along each target cyclic frequency line, and the peak position and peak amplitude distribution are statistically analyzed to generate a modulation feature vector.

[0052] Specifically, based on the two-dimensional cyclic spectrum, the engine's average rotational frequency is first calculated from the tip timing data and mapped onto the cyclic frequency axis of the two-dimensional cyclic spectrum to obtain the rotational frequency index corresponding to the fundamental frequency. Then, integer multiples of this rotational frequency index (e.g., 1, 2, 3, 4, 5 times) are used as target cyclic frequencies, forming a set of target cyclic frequency lines containing five horizontal lines. For each target cyclic frequency line in this set, the energy values ​​corresponding to all frequency grid indices on that line are extracted to form a one-dimensional energy distribution curve. Next, extreme value search is performed on this one-dimensional curve to identify energy peaks. The condition for peak identification is that the energy value at a certain point is greater than the energy values ​​of its two adjacent points, and the energy value at that point must exceed a dynamically set threshold. The threshold is calculated based on the statistical characteristics of the current one-dimensional energy distribution curve. The specific calculation method is as follows: First, calculate the average and standard deviation of all energy values ​​on the curve. Then, set the threshold to the average plus 1.5 times the standard deviation. Only points that simultaneously satisfy the local maximum and are higher than the threshold are identified as valid peaks. For each identified valid peak, record its frequency grid index as the peak position and its normalized energy value as the peak amplitude. Finally, integrate the peak information extracted from all five target cyclic frequency lines and arrange them according to the order of the cyclic frequency lines (1x, 2x, etc.) and the magnitude of the peak amplitude. For example, extract the positions and amplitudes of the three peaks with the largest amplitudes for each line. Then, concatenate the 30 parameters (5 lines x 3 peaks x 2 parameters) of these 5 lines into a fixed-length one-dimensional array in a predetermined order to generate a modulation feature vector.

[0053] The steps for obtaining multi-order blade vibration measurements are as follows:

[0054] Based on the blade tip timing data, the passing time series is organized by blade number, the time interval between adjacent passing is calculated and mapped to the order index, the vibration amplitude of each blade is extracted along the order index and sorted by amplitude, the resonance peak frequency of the corresponding order is located, and the vibration amplitude and resonance peak frequency of the first five engine orders are obtained.

[0055] Based on the vibration amplitude and resonant peak frequency of the first five engine orders, a unique index of blade number and order number is constructed. Vibration amplitude is normalized and resonant peak frequency is standardized. The parameters are spliced ​​together in order of order into a fixed-length parameter sequence and a blade number identifier is added to generate multi-order blade vibration measurement.

[0056] Specifically, based on the blade tip timing data, an independent passage time series is first established for each blade. This series contains the precise timestamps of the blade passing the sensor within multiple consecutive engine rotation cycles. Next, to calculate the vibration characteristics of each blade, a baseline, vibration-free ideal passage time series is determined. This ideal series is calculated by taking the average speed cycle of the engine under stable operating conditions. Then, starting from the first passage time of the first blade, the ideal passage timestamps of all blades in all cycles are recursively derived according to the average speed cycle and blade geometric spacing. The actual passage time series of each blade is subtracted point-by-point from the corresponding ideal passage time series to obtain a time deviation series. This series reflects the early or delayed arrival of the blade tip due to vibration. Multiplying this time deviation series by the instantaneous tangential velocity of the blade tip converts it into a dynamic displacement vibration signal of the blade. Then, for each blade... The displacement vibration signal is subjected to Fast Fourier Transform (FFT) to obtain its vibration spectrum. From this spectrum, based on the engine rotation frequency, the frequency points corresponding to the engine order (1st harmonic, 2nd harmonic, 3rd harmonic, etc.) are identified, and the vibration amplitude at these frequency points is extracted. The vibration amplitudes of all orders are sorted in descending order, and the top five orders with the largest amplitudes are selected. For each of these five orders, a search window with a width of 5% of the theoretical frequency is set on the spectrum with its theoretical frequency as the center. For example, for an order with a theoretical value of 200Hz, the search range is 195Hz to 205Hz. The local maximum energy point is found within this window, and the precise frequency corresponding to this maximum point is the resonance peak frequency of that order. Finally, five data pairs are compiled for each blade, and each data pair contains an engine order vibration amplitude and a corresponding resonance peak frequency, thus obtaining the vibration amplitudes and resonance peak frequencies of the top five engine orders.

[0057] Based on the vibration amplitude and resonant peak frequency of the first five engine orders, a unique index is created for each order data of each blade, consisting of a blade number and an order number, such as "blade 08 - order 03". These data are then standardized. For vibration amplitude, a maximum-minimum normalization method is used. The maximum and minimum values ​​used for normalization are derived from the vibration amplitude envelope under design conditions obtained through statistical analysis of a large amount of historical data collected during long-term engine operation in a healthy state. For example, if historical data shows that the normal range of a certain order vibration amplitude is 0.2mm to 1.5mm, then the normalized value of a currently measured 0.8mm amplitude is (0.8-0.2) / (1.5-0.2), approximately 0.46. For the resonant peak frequency, the deviation from the theoretical frequency of that order is first calculated, and then... The deviation value is Z-score standardized using the formula (x-μ) / σ, where x is the frequency deviation of the current measurement, and μ and σ are the mean and standard deviation of the frequency deviation of this order in the health status historical database, respectively. For example, if the theoretical frequency of a certain order is 200Hz and the measured value is 201.5Hz, the deviation is +1.5Hz. If the mean of the deviation of this order in the health database is +0.5Hz and the standard deviation is 0.4Hz, the standardized value is (1.5-0.5) / 0.4, or 2.5. After normalization and standardization, the normalized vibration amplitude and standardized resonance peak frequency of the first five orders of each blade are spliced ​​together in order from low to high to form a parameter sequence of fixed length 10. Finally, the corresponding blade number is added before the sequence to generate a multi-order blade vibration metric.

[0058] The steps for obtaining the dynamic state points of the integrated blade are as follows:

[0059] Based on the multi-order blade vibration measurement, the corresponding parameters are located in the modulation feature vector according to the blade number and aligned according to the order index. The current working condition multi-dimensional feature points of a single blade are established and numbered and sorted to generate the dynamic state points of the integrated blade.

[0060] Specifically, based on the multi-order blade vibration measurement, the 10-dimensional parameter sequence corresponding to each blade is fused with the modulation feature vector representing the overall dynamic characteristics of the engine. First, the multi-order blade vibration measurement of each blade is matched with the same modulation feature vector by the blade number. This modulation feature vector is a parameter sequence of length 30 extracted from the two-dimensional cyclic spectrum, which reflects the overall vibration modulation information of the engine. The fusion method is feature splicing. Specifically, the 10-dimensional multi-order blade vibration measurement of each blade is taken as the first half and the 30-dimensional modulation feature vector is taken as the second half, and they are sequentially connected to generate a 40-dimensional feature vector for each blade. This 40-dimensional vector defines a point in the feature space, representing the complete dynamic state of the blade under the current operating condition. Finally, the 40-dimensional feature vectors of all blades are arranged according to the physical order of the blades on the compressor disk (e.g., from blade No. 1 to blade No. N) to form an N-row 40-column matrix. This set of matrices is the integrated blade dynamic state point.

[0061] The steps for obtaining the blade state deviation vector are as follows:

[0062] Based on the integrated blade dynamic state points, align them one by one with the dimensions and indices of the preset health state mean vector, verify that the numerical types and dimensions are consistent, remove missing and duplicate records, subtract each component while keeping the original index order unchanged, and obtain the blade state deviation vector.

[0063] Specifically, based on the integrated blade dynamic state points, a preset health state mean vector is first loaded. This vector is obtained through statistical analysis of at least 1000 hours of engine health operation history data. Specifically, under multiple standard test conditions, integrated blade dynamic state point data for all blades are collected. Then, the average value is calculated for each feature dimension (N×40 dimensions in total, where N is the total number of blades) in this massive dataset, forming a benchmark vector representing the center point of the ideal health state. Next, the currently collected integrated blade dynamic state points are aligned with this health state mean vector, comparing their dimensions and index order one by one to ensure a one-to-one correspondence of features. During this process, data quality checks are performed to verify that all values ​​are floating-point numbers and that all feature values ​​are dimensionless. The values ​​are standardized, and then the presence of missing values ​​is checked. Any missing values ​​are filled with the mean of the corresponding dimension in the healthy state mean vector. At the same time, the presence of abnormal duplicate records is detected by calculating the Euclidean distance between adjacent leaf state vectors. If two or more leaf state vectors are exactly the same, they are marked as data collection anomalies, and these duplicate records are replaced with the healthy state mean vector. After data preparation is completed, the 40-dimensional state vector of each current leaf is subtracted element-wise from the corresponding 40-dimensional component in the healthy state mean vector to generate a deviation vector representing the degree of deviation of each leaf from the healthy center. Finally, the deviation vectors of all leaves are concatenated end to end according to the physical order of the leaves to form a single, long-dimensional (N×40-dimensional) vector, which is the leaf state deviation vector.

[0064] The steps for obtaining the engine dynamic anomaly score are as follows:

[0065] Based on the blade state deviation vector, read the health state covariance matrix and verify its symmetry and positive definiteness. If it does not meet the requirements, adjust it only by the minimum amplitude of the diagonal elements until it is positive definite. Calculate the inverse and extract the square root of the diagonal elements as the standard deviation sequence to obtain the inverse matrix of the health state covariance matrix and the standard deviation sequence.

[0066] The engine dynamic anomaly score is calculated based on the inverse of the health status covariance matrix and the standard deviation sequence. The calculation formula is as follows:

[0067]

[0068] Among them, S RA Let P be the engine dynamic anomaly score, Δ be the blade state deviation vector, and P be the inverse of the healthy state covariance matrix. j Let σ be the j-th component of the blade state deviation vector. j w is the square root of the j-th element on the diagonal of the health status covariance matrix. j is the j-th element of the standardized deviation modulation weight, and m is the number of components in the blade state deviation vector.

[0069] Specifically, based on the leaf state deviation vector, the health state covariance matrix is ​​first retrieved from the historical health database. This matrix, also calculated based on the aforementioned 1000 hours of health operation data, describes the linear correlation between various feature dimensions under health conditions and their respective dispersion. After loading, the matrix is ​​immediately validated. The first step is to check its symmetry, i.e., whether the transpose of the matrix is ​​equal to itself. This is an inherent property of the covariance matrix; any asymmetry indicates data corruption or calculation errors. The second step is to verify its positive definiteness by attempting a Jolesi decomposition of the matrix. If the decomposition is successful, the matrix is ​​positive definite; if the decomposition fails, it indicates that the matrix is ​​non-positive definite, possibly due to collinearity in the data or numerical calculation errors. This leads to a regularization process that, when a non-positive definite condition is detected, initiates a regularization procedure by adding a tiny positive value to the main diagonal elements of the matrix. The initial size of this increment is set to one millionth of the smallest positive element on the diagonal. After the increment, the Jolesi decomposition test is performed again. If it still fails, the increment is doubled and the process is repeated until the matrix becomes positive definite. This process ensures the invertibility of the matrix. After confirming that the matrix is ​​positive definite, the Gauss-Jordand elimination method is used to calculate its inverse matrix. At the same time, all elements are extracted from the main diagonal of the regularized (if performed) health state covariance matrix, and the square root of each element is calculated to form a sequence containing the standard deviation of health state for each feature dimension. This yields the inverse matrix and standard deviation sequence of the health state covariance matrix.

[0070] formula: The above formula calculates a comprehensive engine dynamic anomaly score by combining the overall deviation of multidimensional features with the weighted individual deviation of key features. The first part of the formula... It is the Mahalanobis distance, which measures the distance from the current engine state point Δ to the center of the healthy state distribution and takes into account the correlation between features. It reflects the true degree of anomaly better than the Euclidean distance. The second part of the formula, [1+tanh(…)], is a dynamic modulation factor, which is determined based on the degree to which each feature deviates from its normal range (by…). (measure) and the importance of the feature (by weight w) j To achieve nonlinear amplification, the design employs a hyperbolic tangent function tanh to map the weighted average deviation to the (-1,1) interval. Adding 1 changes the modulation factor range to (0,2). This design allows the anomaly score to be amplified even when the overall Mahalanobis distance is small, even if certain key features deviate significantly, thereby improving the sensitivity to specific early failure modes.

[0071] Δ is the blade state deviation vector, representing the difference between the current engine state and the preset healthy state mean vector. This vector is obtained through the aforementioned step of "obtaining the blade state deviation vector based on the integrated blade dynamic state points...". Its dimension m is equal to the number of blades multiplied by the number of features of each blade. In this example, to simplify the calculation, four key feature dimensions of one blade are selected for demonstration, i.e., m = 4. A specific deviation vector is obtained through calculation.

[0072] P is the inverse of the health state covariance matrix, which reflects the decorrelation between features under healthy conditions. It is obtained through the aforementioned step of "obtaining the inverse of the health state covariance matrix and the standard deviation sequence based on the leaf state deviation vector...". For the four features in this example, the health state covariance matrix C calculated from the health database is... This matrix is ​​a symmetric positive definite matrix, and its inverse matrix P = C -1 The calculation yielded the following result:

[0073] Δ j Let Δ be the j-th component of the blade state deviation vector, and let Δ be the j-th element in the vector Δ. For example, Δ1 = 0.3, Δ2 = -0.5, Δ3 = 1.2, Δ4 = 0.1.

[0074] σ j Let $\frac{j}{j}$ be the square root of the $j$-th element on the diagonal of the health state covariance matrix, which is the standard deviation of the $j$-th feature under healthy conditions. This is obtained by taking the square root of the main diagonal elements of the health state covariance matrix $C$.

[0075] w j The j-th element of the standardized deviation modulated weights represents the importance of the j-th feature in assessing engine anomalies. This weight is determined based on the risk priority number corresponding to each feature in engine failure mode and effects analysis. The RPN integrates the probability, severity, and detectability of failure occurrence; a higher value indicates that the monitoring of the corresponding feature is more critical. The weight calculation formula is as follows: In this example, the four characteristics correspond to the first-order bending vibration amplitude, first-order torsional vibration frequency, aerodynamic load imbalance, and tip clearance variation of the blade, respectively. Based on expert evaluation and historical data analysis, their RPN values ​​are 80, 65, 90, and 40, respectively. The calculated weights are as follows:

[0076] m is the number of components in the blade state deviation vector, i.e., the total dimension of the feature. In this example, m = 4.

[0077] Calculations based on parameters:

[0078] First, calculate the Mahalanobis distance.

[0079]

[0080] Next, the modulation factor part [1+tanh(…)] is calculated:

[0081] Calculate the weighted terms

[0082]

[0083]

[0084] Calculate ∑w j and

[0085]

[0086] Calculate the modulation factor:

[0087]

[0088] Modulation factor = 1 + 0.485 = 1.485;

[0089] Finally, the engine dynamic anomaly score S is calculated. RA :

[0090] S RA =1.096·1.485≈1.627;

[0091] The results indicate that the calculated engine dynamic anomaly score is 1.627. This value combines the Mahalanobis distance between the current and healthy states in the multidimensional feature space with the weighted deviation of key features. Typically, a threshold system based on historical data and expert experience is preset to interpret this score; for example, a normal range is set to S. RA <1.5, the scope of concern is 1.5≤S RA <3.0, alarm range is S RA A score of ≥3.0 and the current score of 1.627 fall within the "range of concern," meaning that the engine's dynamic state has shown a slight abnormality. Although it has not yet reached the severity required to shut down the engine immediately, it has deviated from the normal health baseline. Operation and maintenance personnel need to focus on monitoring and further diagnostic analysis of the relevant characteristics, especially the third characteristic (aerodynamic load imbalance) which contributes the most to the score.

[0092] The above are merely preferred embodiments of the present invention and are not intended to limit the present invention in any other way. Any person skilled in the art may make changes or modifications to the above-disclosed technical content to create equivalent embodiments that can be applied to other fields. However, any simple modifications, equivalent changes, and modifications made to the above embodiments based on the technical essence of the present invention without departing from the scope of the present invention shall still fall within the protection scope of the present invention.

Claims

1. A method of detecting performance of an aeroengine, characterized in that, The method comprises the following steps: Collecting vibration signals installed on the engine shell and blade tip timing data on the periphery of the compressor blades, segmenting the continuous vibration signals according to the engine speed cycle, calculating the time-frequency distribution of each signal segment, and obtaining a time-frequency energy matrix; According to the time-frequency energy matrix, the spectral density in the frequency domain and the cyclic frequency domain is calculated, a two-dimensional cyclic spectrum diagram is established, the energy peak value and distribution along the target cyclic frequency line are extracted from the two-dimensional cyclic spectrum diagram, and a modulation feature vector is obtained; According to the blade tip timing data, the first five engine order vibration amplitudes and resonance peak frequencies of each compressor blade are extracted, the parameters are combined to form a multi-order blade vibration metric, the modulation feature vector and the multi-order blade vibration metric are fused, and an integrated blade dynamic state point is obtained; According to the integrated blade dynamic state point, the deviation between the preset health state mean vector is calculated, a blade state deviation vector is generated, and the engine dynamic abnormality score is calculated by using the blade state deviation vector and the health state covariance matrix; The acquisition step of the two-dimensional cyclic spectrum diagram is: According to the time-frequency energy matrix, the energy is accumulated in the time index direction according to the frequency grid index and is amplitude normalized, the energy periodicity is counted in the frequency grid direction according to the time index and is rearranged according to the cyclic frequency index, the frequency domain spectral density and the cyclic frequency domain spectral density are generated; According to the frequency domain spectral density and the cyclic frequency domain spectral density, a two-dimensional coordinate grid of the cyclic frequency index and the frequency grid index is established, the energy density is filled according to the coordinate grid, and the boundary is interpolated and the discrete points are smoothed, and a two-dimensional cyclic spectrum diagram is generated.

2. The method of claim 1, wherein, The acquisition step of the time-frequency energy matrix is: According to the acceleration channel installed on the engine shell and the blade tip timing channel on the periphery of the compressor blades, the vibration signal timestamp sequence and the vibration signal amplitude sequence are read, the blade tip timing timestamp sequence and the blade tip timing identification sequence are read, the time reference is unified and arranged according to the collection sequence, and the vibration signal and the blade tip timing data are generated; According to the vibration signal and the blade tip timing data, the adjacent same blade passing interval is calculated by analyzing the blade tip timing timestamp sequence, the engine speed cycle sequence is generated, the signal segment sequence is formed by positioning the division point on the vibration signal timestamp sequence according to the engine speed cycle sequence and slicing the vibration signal amplitude sequence; According to the signal segment sequence, a fixed time window length and a time window overlap ratio are set to slide and count each signal segment, the energy of each frequency grid in each time window is counted and the time index corresponding to the time window center is recorded, the energy values are arranged according to the time index and the frequency grid index to form a two-dimensional structure, and a time-frequency energy matrix is generated.

3. The method of claim 1, wherein, The acquisition step of the modulation feature vector is: According to the two-dimensional cyclic spectrum diagram, a target cyclic frequency line set is formed according to the engine rotation frequency index and the integer multiple index, equal-step sampling and extreme value searching are performed along each target cyclic frequency line, the peak position and peak amplitude distribution are counted, and a modulation feature vector is generated.

4. The method of claim 1, wherein, The acquisition step of the multi-order blade vibration metric is: According to the blade tip timing data, the through time sequence is arranged according to the blade number, the adjacent through time interval is calculated and the order index is mapped, the vibration amplitude of each blade is extracted along the order index and is sorted according to the amplitude, the corresponding resonance peak frequency of the order is located, and the first five engine order vibration amplitudes and resonance peak frequencies are obtained. According to the first five engine order vibration amplitudes and resonance peak frequencies, the unique index of the blade number and the order number is constructed, the vibration amplitude is normalized and the resonance peak frequency is standardized, the parameter sequence with a fixed length is spliced in the order sequence and the blade number identification is added, and the multi-order blade vibration metric is generated.

5. The method of claim 1, wherein, The acquisition step of the integrated blade dynamic state point is: According to the multi-order blade vibration metric, the corresponding parameters in the modulation feature vector are located according to the blade number and are aligned according to the order index, the current working condition multi-dimensional feature point of a single blade is established and is numbered and sorted, and the integrated blade dynamic state point is generated.

6. The method of claim 1, wherein, The acquisition step of the blade state deviation vector is: According to the integrated blade dynamic state point, the dimensions and indexes of the preset health state mean vector are aligned one by one, the numerical type and the dimension are checked, the missing and repeated records are cleared, the components are subtracted and the original index order is kept unchanged, and the blade state deviation vector is obtained.

7. The method of claim 1, wherein, The acquisition step of the engine dynamic abnormality score is: According to the blade state deviation vector, the health state covariance matrix is read and the symmetry and positive definiteness are checked, if not satisfied, only the diagonal element is adjusted to positive definite with the minimum amplitude, the inverse is calculated and the square root of the diagonal element is extracted as the standard deviation sequence, the health state covariance matrix inverse matrix and the standard deviation sequence are obtained.

8. The method of claim 7, wherein, The acquisition step of the engine dynamic abnormality score further includes: according to the health state covariance matrix inverse matrix and the standard deviation sequence, the engine dynamic abnormality score is calculated.

9. The method of claim 8, wherein, The calculation formula of the engine dynamic abnormality score is: ； in, This represents the engine dynamic anomaly score. This is the blade state deviation vector. Let the inverse of the health status covariance matrix be... The first of the blade state deviation vectors One portion, The diagonal of the health status covariance matrix The square root of each element. The first standardized deviation modulation weight One element, This represents the number of components in the blade state deviation vector.

Images