Aero-engine performance detection method

Abstract

The invention relates to the technical field of engine testing, in particular to an aero-engine performance detection method, which comprises the following steps of: acquiring vibration signals mounted on an engine shell and blade tip timing data at the periphery of a gas compressor blade, and performing segmented processing on the continuous vibration signals according to an engine rotating speed period; and calculating time-frequency distribution of each signal segment to obtain a time-frequency energy matrix. According to the method, the vibration signals and the blade tip timing data are synchronously collected and segmented processing is carried out in combination with the rotating speed period, so that the frequency characteristics in different time periods can be completely extracted, and information missing caused by single time domain or frequency domain processing in a traditional method is avoided; furthermore, spectral density is calculated on a frequency domain and a cyclic frequency domain through a time frequency energy matrix, a two-dimensional spectrogram capable of reflecting energy distribution and cyclic characteristics at the same time is constructed, and full expression of characteristic information under multiple dimensions is guaranteed.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of engine testing, and particularly relates to an aero-engine performance detection method. BACKGROUND

[0002] As a key power device of an aircraft, the performance state of an aero-engine directly affects the operation safety of the aircraft and the safety of the crew on board, and therefore it is of great significance to detect the abnormality of the performance of the engine.

[0003] The prior art takes the overall vibration level or the numerical offset of a local parameter as the basis for judgment, resulting in a recognition blind area for the characteristic performance under a complex coupling state, and the threshold value of the overall parameter is easily affected by environmental noise and random disturbances to produce misjudgments, and the abnormal response of a local blade or a specific order is not sensitive enough to discover early local faults in time, for example, when the vibration amplitude of a certain order slightly increases, the overall vibration value may not significantly exceed the threshold value, thereby concealing potential hazards and delaying maintenance opportunities. Therefore, improvement is needed. SUMMARY

[0004] The purpose of the present application is to solve the problems existing in the prior art and to provide an aero-engine performance detection method.

[0005] In order to achieve the above-mentioned purpose, the present application adopts the following technical scheme, an aero-engine performance detection method, comprising the following steps:

[0006] Collecting vibration signals installed on the engine casing and blade tip timing data 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;

[0007] According to the time-frequency energy matrix, calculating the spectral density in the frequency domain and the cyclic frequency domain, establishing a two-dimensional cyclic spectrum graph, extracting the energy peak value and distribution along the target cyclic frequency line from the two-dimensional cyclic spectrum graph, and obtaining a modulation feature vector;

[0008] According to the blade tip timing data, extracting the first five engine order vibration amplitudes and resonance peak frequencies of each compressor blade, combining the parameters to form a multi-order blade vibration metric, fusing the modulation feature vector and the multi-order blade vibration metric, and obtaining an integrated blade dynamic state point;

[0009] According to the integrated blade dynamic state point, calculating the deviation from the preset health state mean vector, generating a blade state deviation vector, and using the blade state deviation vector and the health state covariance matrix to calculate an engine dynamic abnormality score.

[0010] Preferably, the time-frequency energy matrix obtaining step is:

[0011] According to the acceleration channel installed on the engine shell and the blade tip timing channel on the periphery of the compressor blade, the vibration signal timestamp sequence and the vibration signal amplitude sequence are read, the blade tip timing timestamp sequence and the blade tip timing identifier sequence are read, the time reference is unified and arranged in the collection order, and the vibration signal and the blade tip timing data are generated;

[0012] According to the vibration signal and the blade tip timing data, the adjacent same blade passing interval is calculated based on the blade tip timing timestamp sequence, the engine speed period sequence is generated, the vibration signal amplitude sequence is segmented based on the engine speed period sequence on the vibration signal timestamp sequence, and the signal segment sequence is formed;

[0013] 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 bin in 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 the frequency bin index to form a two-dimensional structure, and the time-frequency energy matrix is generated.

[0014] Preferably, the two-dimensional cyclic spectrum diagram obtaining step is:

[0015] According to the time-frequency energy matrix, the energy is accumulated in the time index direction according to the frequency bin index and the amplitude is normalized, the energy periodicity is counted in the frequency bin 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;

[0016] 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 bin 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 the two-dimensional cyclic spectrum diagram is generated.

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

[0018] 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 the peak amplitude distribution are counted, and the modulation feature vector is generated.

[0019] Preferably, the multi-order blade vibration metric obtaining step is:

[0020] According to the blade tip timing data, the passing time sequence is arranged according to the blade number, the adjacent passing time interval is calculated and the order index is mapped, the vibration amplitude of each blade is extracted along the order index and sorted according to the amplitude, the resonance peak frequency corresponding to the order is located, and the first five engine order vibration amplitudes and resonance peak frequencies are obtained;

[0021] According to the first five engine order vibration amplitude and resonance peak frequency, the unique index of blade number and order number is constructed, the vibration amplitude is normalized and the resonance peak frequency is standardized, the parameter sequence with fixed length is spliced in order and the blade number is marked, and the multi-order blade vibration metric is generated.

[0022] Preferably, the integrated blade dynamic state point acquisition step is:

[0023] According to the multi-order blade vibration metric, the corresponding parameters in the modulation feature vector are located according to the blade number and aligned according to the order index, the current working condition multi-dimensional feature point of a single blade is established and numbered and sorted, and the integrated blade dynamic state point is generated.

[0024] Preferably, the blade state deviation vector acquisition step is:

[0025] 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 dimension are checked for consistency, 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.

[0026] Preferably, the engine dynamic abnormality score acquisition step is:

[0027] According to the blade state deviation vector, the health state covariance matrix is read and checked for symmetry and positive definiteness, if not, only the diagonal elements are adjusted to positive definiteness with the minimum amplitude, the inverse is calculated and the square root of the diagonal elements is extracted as the standard deviation sequence, the health state covariance matrix inverse matrix and the standard deviation sequence are obtained.

[0028] According to the health state covariance matrix inverse matrix and the standard deviation sequence, the engine dynamic abnormality score is calculated.

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

[0030] In the present application, by synchronously collecting the vibration signal and the blade tip timing data and combining the rotation speed period for segmented processing, the frequency characteristics in different time periods can be completely extracted, avoiding the information missed in the traditional method due to single time domain or frequency domain processing. Further, the spectral density is calculated on the frequency domain and the cyclic frequency domain through the time-frequency energy matrix, a two-dimensional spectrum graph reflecting the energy distribution and the cyclic characteristics is constructed, the sufficient expression of the characteristic information in multiple dimensions is ensured, the modulation feature vector is obtained by extracting the energy peak value and the distribution of the target cyclic frequency line, the periodic vibration and the harmonic effect are captured, the multi-order blade vibration metric is formed by extracting the multi-order vibration amplitude and the resonance frequency of each blade using the blade tip timing data, and the modulation feature vector is fused, so that the dynamic state of a single blade can be quantified on multiple feature points. Finally, by comparing the deviation between the integrated blade dynamic state points and the health state mean value, and combining the comprehensive calculation of the health state covariance matrix, the abnormal score generated can quantify the degree of overall performance deviation, realize the identification of performance abnormalities under complex coupled working conditions, which not only improves the accuracy and robustness of abnormal detection, but also enhances the response capability to different levels of feature deviation, so that the potential faults of the engine can be identified and quantified at an early stage. BRIEF DESCRIPTION OF DRAWINGS

[0031] Figure 1 The figure is a schematic diagram of the steps of the present application. DETAILED DESCRIPTION

[0032] In order to make the purpose, technical scheme and advantages of the present application more clear, the present application is further described in detail below in combination with the drawings and examples. It should be understood that the specific examples described herein are only used to explain the present application and do not limit the present application.

[0033] Please refer to Figure 1 The present application provides a technical scheme, an aero-engine performance detection method, comprising the following steps:

[0034] Collecting the vibration signal installed on the engine shell and the blade tip timing data around the compressor blade, segmenting the continuous vibration signal according to the engine rotation speed period, calculating the time-frequency distribution of each signal segment, and obtaining the time-frequency energy matrix;

[0035] According to the time-frequency energy matrix, the spectral density on the frequency domain and the cyclic frequency domain is calculated, a two-dimensional cyclic spectrum graph is established, the energy peak value and the distribution along the target cyclic frequency line are extracted from the two-dimensional cyclic spectrum graph, and the modulation feature vector is obtained;

[0036] 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 to obtain an integrated blade dynamic state point;

[0037] 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.

[0038] The acquisition step of the time-frequency energy matrix is:

[0039] According to the acceleration channel installed on the engine shell and the blade tip timing channel around the compressor blade, 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 order, and the vibration signal and the blade tip timing data are generated;

[0040] According to the vibration signal and the blade tip timing data, the adjacent same blade passing interval is calculated according to the blade tip timing timestamp sequence, the engine speed period sequence is generated, the vibration signal amplitude sequence is segmented according to the engine speed period sequence, and the signal segment sequence is formed;

[0041] 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 bin in 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 the frequency bin index to form a two-dimensional structure, and the time-frequency energy matrix is generated.

[0042] Specifically, according to the acceleration channel installed on the engine shell and the blade tip timing channel around the compressor blade periphery, first, the time reference alignment of multi-source data is performed, specifically, the acquisition devices of the acceleration sensor and the blade tip timing sensor are synchronized by using a unified global positioning system clock signal, and the error of all time stamps is ensured to be within 1 microsecond, then the time stamp sequence and the corresponding amplitude sequence of the vibration signal, and the time stamp sequence and the blade identification sequence of the blade tip timing are read in parallel, and a preliminary data quality check is performed during the reading process, for example, checking whether the vibration signal amplitude exceeds the range of the sensor, setting a saturation threshold, which is determined according to the sensor specification, for example, for an acceleration sensor with a range of ±50g, the saturation threshold can be 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 of adjacent time stamps with the expected interval calculated in the last speed period, if the deviation exceeds the preset tolerance, the tolerance is set to 5% of the average speed period, for example, at 6000RPM, the period 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, the saturated vibration signal point is filled by linear interpolation of adjacent valid points, and the missing blade tip timing data point is estimated based on the linear extrapolation of the passing time of the same blade before and after, after alignment and verification, the vibration signal and the blade tip timing data containing the synchronous time stamp are generated.

[0043] According to the vibration signal and the blade tip timing data, the passing time of each blade is extracted from the blade tip timing timestamp sequence, the specific blade is identified through the blade tip timing identification sequence, the difference between the time when the blade passes through the same sensor twice is calculated, and this difference is the speed period of the engine at that time, in order to improve stability and accuracy, the speed periods calculated by all blades in a complete rotation period are averaged to obtain the average speed period of the current rotation, all average speed periods of the rotation are arranged in time sequence to form the engine speed period sequence, then the continuous vibration signal timestamp sequence is segmented using the engine speed period sequence, and the starting timestamp and the ending timestamp of each speed period are used as boundaries, the data of the corresponding time period is extracted from the vibration signal amplitude sequence, so that the long-term vibration signal is cut into a series of signal segments synchronized with the engine rotation period to form a signal segment sequence.

[0044] According to 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 selected based on the balance between the desired frequency resolution and 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 time window overlap ratio is set to an empirical value, which is 75% in this case. The Hanning window is selected as the window function. Then, the configured sliding time window is applied to each signal segment, and is shifted from the beginning to the end. For each intercepted signal segment in the window, a fast Fourier transform is performed to calculate the complex values of the frequency components. The energy value is obtained by taking the square of the modulus of the complex values, and the energy values are distributed to the corresponding frequency bins. At the same time, the time index corresponding to the center point of each time window is recorded. Finally, the energy values calculated by all time windows are arranged according to the time index and frequency bin index to form a two-dimensional data structure, where one dimension represents time and the other dimension represents frequency. The value of each element is the energy density of the frequency at that time. Thus, a time-frequency energy matrix is generated.

[0045] The acquisition steps of the two-dimensional cyclic spectrum diagram are as follows:

[0046] According to the time-frequency energy matrix, the energy is accumulated in the time index direction according to the frequency bin index and is amplitude-normalized. The energy periodicity is counted in the frequency bin 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.

[0047] 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 bin index is established. The energy density is filled according to the coordinate grid, and is interpolated at the boundary and smoothed at the discrete points. A two-dimensional cyclic spectrum diagram is generated.

[0048] Specifically, according to the time-frequency energy matrix, for calculating the frequency domain spectral density, first, sum all the energy values under each independent frequency bin index along the time index direction, the specific operation is, for each row in the matrix, add the numerical values of all columns in the row to obtain a scalar value representing the total energy of the frequency, after performing this operation on all rows, a one-dimensional array is obtained, the length of which is equal to the total number of frequency bins, and each element in the array corresponds to the total energy of a frequency bin, then the one-dimensional array is subjected to amplitude normalization processing, the maximum value of all elements in the array is calculated, and then each element in the array is divided by the maximum value, so that all energy values are scaled to the interval of 0 to 1, to obtain the frequency domain spectral density, then, for calculating the cyclic frequency domain spectral density, periodic analysis of the energy change with time in the time-frequency energy matrix is required, for each frequency bin index, that is, each row of the matrix, it is regarded as a time series, and the Fourier transform is performed on the time series to calculate the spectral components at different cyclic frequencies, the calculation formula of this process is Where S(a,f) is the cyclic spectral density function, a is the cyclic frequency, f is the carrier frequency (corresponding to the frequency bin 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 length of the observation time, in discrete implementation, this is achieved by applying the fast Fourier transform (FFT) to each row of data, and the result of the transformation is a complex number matrix, the amplitude part of which is the energy distribution at different cyclic frequencies, and the result is rearranged according to the cyclic frequency index to generate the cyclic frequency domain spectral density.

[0049] According to 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) of which represents the frequency grid index, ranging from the lowest frequency to the highest frequency, and the vertical axis (Y-axis) of which represents the cyclic frequency index, covering all calculated cyclic frequencies, then the energy density values in the cyclic frequency domain spectral density matrix are filled one by one into the corresponding coordinate points of the two-dimensional coordinate grid to form a discrete energy distribution graph, due to the discreteness of the calculation process, there may be data missing or discontinuity in the boundary area of the image, therefore, boundary interpolation processing is needed, a bilinear interpolation algorithm is adopted, for any coordinate point on the grid boundary without numerical value, the nearest four adjacent points with numerical value around it are found, and the interpolation value of the point is calculated by weighted average according to the distance from the point to the four adjacent points, the specific weight is inversely proportional to the distance, at the same time, in order to remove the random noise introduced by discrete calculation and make the features in the spectrum more smooth, a two-dimensional Gaussian smoothing filter is applied to the entire filled and interpolated grid data, the size of the Gaussian kernel is set to 5x5 pixels, and the standard deviation σ is set to 1.5, the size of the kernel and the standard deviation are selected according to the resolution and noise level of the spectrum by experience, for example, if the spectrum resolution is high, a smaller kernel can be selected, otherwise a larger kernel is selected, the standard deviation controls the degree of smoothing, the larger the value, the stronger the smoothing effect, the 5x5 Gaussian kernel is convolved with each point and its neighborhood in the coordinate grid, and the weighted average value calculated is used to replace the original value of the center point, after boundary interpolation and global smoothing, the final two-dimensional cyclic spectrum graph is generated.

[0050] The acquisition step of the modulation feature vector is:

[0051] According to the two-dimensional cyclic spectrum graph, 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, and the peak position and peak amplitude distribution are counted to generate the 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 for testing the performance of an aero-engine, characterized in that, Includes the following steps: 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. 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. 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. 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.

2. The method for testing the performance of an aero-engine according to claim 1, characterized in that, The steps for obtaining the time-frequency energy matrix are as follows: 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. 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. 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.

3. The method for testing the performance of an aero-engine according to claim 1, characterized in that, The steps for obtaining the two-dimensional cyclic spectrum are as follows: 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. 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.

4. The method for testing the performance of an aero-engine according to claim 1, characterized in that, The steps for obtaining the modulation feature vector are as follows: 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.

5. The method for testing the performance of an aero-engine according to claim 1, characterized in that, The steps for obtaining the multi-order blade vibration measurement are as follows: 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. 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.

6. The method for testing the performance of an aero-engine according to claim 1, characterized in that, The steps for obtaining the dynamic state points of the integrated blade are as follows: 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.

7. The method for testing the performance of an aero-engine according to claim 1, characterized in that, The steps for obtaining the blade state deviation vector are as follows: 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.

8. The method for testing the performance of an aero-engine according to claim 1, characterized in that, The steps for obtaining the engine dynamic anomaly score are as follows: Based on the blade state deviation vector, the health state covariance matrix is ​​read and its symmetry and positive definiteness are verified. If not, it is adjusted only by the minimum amplitude of the diagonal elements until it becomes 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.

9. The method for testing the performance of an aero-engine according to claim 8, characterized in that, The step of obtaining the engine dynamic anomaly score further includes: calculating the engine dynamic anomaly score based on the inverse matrix of the health state covariance matrix and the standard deviation sequence.

10. The method for testing the performance of an aero-engine according to claim 9, characterized in that, The formula for calculating the engine dynamic anomaly score is as follows: 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.

Images