A Spectrum Analysis Method for Undersampled Harmonic Signals

Abstract

This invention provides a spectral analysis method for undersampled harmonic signals, comprising the following steps: acquiring multi-channel leaf tip timing signals; constructing multiple snapshot matrices based on the multi-channel leaf tip timing signals; calculating a cross-correlation matrix; calculating the following power expression of the cross-correlation matrix; calculating the steering vector; and obtaining the function beamforming output. The technical solution of this invention can be widely applied in the field of harmonic signal processing technology to solve the problem of processing undersampled harmonic signals.

Description

Technical Field

[0001] This invention relates to a method for processing undersampled harmonic signals, particularly a method for processing undersampled signals when identifying the natural frequency of a rotating blade. Background Technology

[0002] Rotating blades are widely used in turbomachinery and are an important component of it. Currently, testing techniques for the vibration of rotating blades during operation can be broadly categorized into contact and non-contact methods. Contact testing typically involves attaching strain gauges to the blade surface and using slip rings or wireless transmission devices to introduce the test signal. However, these contact testing methods not only result in structural complexity and may affect the pre-designed flow field, but also make it difficult to monitor all blades over extended periods.

[0003] Blade tip timing (BTT) technology, as a non-contact testing method for rotating blades, has received widespread attention in recent years. Its advantage lies in its ability to perform long-term testing and monitoring of a row of rotating blades using a non-contact approach. The core idea of ​​blade tip timing technology is to calculate the difference between the theoretical and actual arrival times of the bond phase pulses for each blade, thereby calculating the blade tip displacement. When sampling the blade tip displacement of a rotating blade using blade tip timing, each blade is sampled only once per revolution on each blade tip timing sensor. Due to the limited space on the casing wall, it is difficult to install a large number of dedicated blade tip timing measurement points, resulting in the blade tip displacement measured by the blade tip timing signal generally being undersampled.

[0004] The natural frequency identification of rotating blades can be used to monitor blade condition; when a blade cracks, its natural frequency changes. However, due to the undersampling characteristics of the BTT signal, the spectrum obtained directly from the BTT signal is aliased, making it impossible to obtain accurate natural frequency identification results.

[0005] Therefore, suppressing spectral aliasing and obtaining accurate natural frequencies of blades has been extensively studied in BTT signal processing. The single-degree-of-freedom method can obtain the blade's natural frequencies using a single BTT sensor without post-processing, but this method must be performed under frequency sweep conditions. When identifying blade natural frequencies using MUSIC (Multiple Signal Classification) methods, the number of natural frequencies of the blade needs to be estimated in advance to obtain the noise subspace. In reality, this information is unknown. In MUSIC methods, if the estimated number of frequencies is too low, a certain natural frequency will be lost; if the estimated number of natural frequencies is too high, interference will not be completely removed, affecting the performance of the MUSIC algorithm in identifying natural frequencies of rotating blades. The actual tip displacement signal of a rotating blade is complex, containing not only synchronous vibration components related to the rotational frequency but also asynchronous vibration components related to the blade's natural frequencies, making it difficult to estimate the number of frequencies. Summary of the Invention

[0006] To address the aforementioned problems with existing non-contact testing techniques for rotating blades, this invention provides a spectral analysis method for undersampled harmonic signals.

[0007] The technical solution of the present invention is as follows:

[0008] A method for spectral analysis of undersampled harmonic signals, comprising the following steps:

[0009] S1. Acquire multi-channel leaf tip timing signals;

[0010] S2. Construct a multi-block snapshot matrix based on the multi-channel leaf tip timing signal:

[0011] Y = AX + W

[0012] in, Let K be the number of snapshots in the matrix. This is the transfer matrix from frequency to tip displacement. Let W be the pseudo-amplitude matrix of the frequency to be solved, and W be Gaussian white noise.

[0013] S3. Calculate the cross-correlation matrix C:

[0014]

[0015] S4. Calculate the following power expression for the cross-correlation matrix:

[0016]

[0017] Where U is a unitary matrix, σ is a singular value of the matrix, and ξ is an exponent greater than 1;

[0018] S5, Calculate frequency f p The guide vector s p :

[0019]

[0020] Wherein, d1…d M This is the delay factor;

[0021] S6, Function Beamforming Output FB(f p ):

[0022]

[0023] Optionally, the tip displacement signal y obtained in step S1 from the m-th sensor during the l-th revolution. m [l] is:

[0024]

[0025] Where T is the sampling period, P is the number of frequency components of the blade, and c p Let be the coefficient of the p-th frequency component. For the corresponding phase,

[0026]

[0027] f h d represents the highest frequency expected in the analysis. m This is the delay factor.

[0028] Optionally, step S2 further includes: dividing matrix A into blocks, with M rows forming one block, for a total of M-K+1 blocks, and the i-th block A{i} being:

[0029]

[0030] Optionally, matrix X is:

[0031]

[0032] Optionally, in step S4, ξ is greater than or equal to 100.

[0033] The technical effects of this invention are as follows:

[0034] The spectrum analysis method for undersampled harmonic signals of this invention is based on function beamforming. It introduces power operations into the measurement data and conventional beamforming to suppress spectral aliasing while preserving the energy of the target frequency components. Furthermore, this method eliminates the need to estimate the number of frequencies in the BTT signal, effectively avoiding the omission of blade natural frequencies or the introduction of additional noise. In summary, the technical solution of this invention achieves its objective.

[0035] The further effects of the above-mentioned alternative methods will be explained in detail below with reference to specific implementation methods. Attached Figure Description

[0036] Figure 1 A schematic diagram illustrating the principle of timed measurement of blade tip displacement.

[0037] Figure 2 Construct a graph for a multi-block snapshot matrix.

[0038] Figure 3 This is a flowchart of the method of the present invention.

[0039] Figure 4 This is a comparison chart of the recognition success rates of various methods in the simulation experiment.

[0040] Figure 5 This is a visualization of the simulation results of various methods under different signal-to-noise ratios.

[0041] Figure 6 The recognition success rate of function beamforming under different ξ values.

[0042] Figure 7 This is a visualization of the simulation results for different values ​​of ξ.

[0043] Figure 8 The average running time for beamforming of the function under different ξ values.

[0044] Figure 9 This is a diagram showing the overall layout of the experiment.

[0045] Figure 10 The images show the time-frequency diagram and spectrum of the laser displacement sensor signal under frequency sweeping conditions.

[0046] Figure 11 The spectrum of the laser displacement sensor at 78% of its maximum rotational speed.

[0047] Figure 12 The image shows the frequency identification results of the MUSIC algorithm under different numbers of harmonics in the experiment.

[0048] Figure 13 The image shows the identification results of conventional beamforming and beamforming algorithms with different ξ values ​​in the experiment. Detailed Implementation

[0049] The technical solution of the present invention will be described in detail below with reference to the accompanying drawings.

[0050] Tip timing indirectly measures tip displacement by comparing the theoretical and actual arrival times of the blade tip. For example... Figure 1As shown, using an OPR (Once per-revolution) sensor as a reference, its circumferential azimuth angle is assumed to be 0°. For simplicity, only the tip displacement of the b-th blade measured by the m-th sensor is considered. When the blade is not vibrating, the time it takes for blade b to rotate to sensor m can be determined as the theoretical arrival pulse of that blade. During the actual rotation of the blade, vibration will occur due to aerodynamic forces and other excitations, causing the actual arrival pulse of blade b measured by the m-th sensor to lead or lag. The tip displacement can be calculated by combining this with the blade disk rotation speed.

[0051] Assuming the sampling period is T, then the arrival time t of the m-th sensor is... m It can be represented as:

[0052]

[0053] In equation (1), θ m Let be the circumferential azimuth angle of the m-th sensor. The tip displacement signal acquired by the m-th sensor in the l-th revolution can be expressed as:

[0054]

[0055] In equation (2), l = 0, 1, ..., L-1 represents the number of revolutions the blade makes, P represents the number of frequency components of the blade, and c p Let be the coefficient of the p-th frequency component. Let be the corresponding phase.

[0056]

[0057] In equation (3), f h Let be the highest frequency expected to be analyzed. Further, t can be... m Represented as T0 and delay coefficient d m The product form:

[0058] t m =d m T0. (4)

[0059] Substituting formula (4) into formula (2), we get:

[0060]

[0061] Considering all M sensors, the construction method of the multi-block snapshot matrix is ​​as follows: Figure 2 As shown.

[0062] The tip displacement signal acquired by BTT can be represented in matrix form as follows:

[0063] Y = AX + W (6)

[0064] In formula (6) This is a multi-snapshot matrix of blade tip displacement signals acquired by BTT, where K is the number of snapshots and N is the number of blade rotations. Figure 2 The number of columns in a multi-block snapshot matrix. This is the transfer matrix from frequency to tip displacement. Let Y be the pseudo-amplitude of the frequency to be solved, and W be Gaussian white noise with the same dimension as Y. Divide matrix A into N-K+1 blocks, with each block consisting of M rows. The i-th block A{i} can be represented as:

[0065]

[0066] The pseudo-amplitude matrix X for each frequency can be represented as:

[0067]

[0068] Thus, the problem of identifying the natural frequency of the rotating blade using the BTT signal is transformed into the process of finding X given Y and A in equation (6) when noise W is present.

[0069] Conventional beamforming is a method widely used in radar and sound source identification. Its core idea is to focus the signal measured by the sensor onto each scan point of interest using a steering vector, obtaining the relative energy magnitude of each scan point to determine the target's direction or location. In this invention, beamforming is used in BTT signal processing, where each frequency point of interest is considered a scan point. A steering vector is constructed from the BTT test data Y to each frequency of interest based on matrix A, thereby obtaining the relative energy magnitude at each frequency and ultimately identifying the blade's natural frequency.

[0070] When using beamforming, the cross-correlation matrix C of the BTT signal measurement must first be calculated:

[0071]

[0072] For a single frequency f p In terms of its guiding vector s p It can be written as:

[0073]

[0074] Based on the steering vector and cross-correlation matrix, the frequency f can be obtained. p Beamforming output B(f) at the location p ):

[0075]

[0076] In equation (11), ||·||2 represents the 2-norm of the vector. The natural frequency of the blade can be determined by the output magnitude of the beamforming at different frequencies.

[0077] However, timing the blade tip displacement is a severely undersampled process due to the limited number of sensors that can be installed on the casing wall. When the sampling rate is insufficient, beamforming will produce severe frequency aliasing, affecting the identification of the natural frequencies of the rotating blades.

[0078] The point propagation function describes the contribution of a single frequency component to the beamforming output. Frequency component f p The energy-mean-inspired average tip displacement contribution induced by the sensor can be expressed as:

[0079]

[0080] In equation (12), a p Let f represent the p-th column of matrix A. Assuming the components at different frequency levels are uncorrelated, then in the measured cross-correlation matrix C, the frequency component f... p Amount of contribution for:

[0081]

[0082] For the entire matrix C, we have:

[0083]

[0084] Taking into account equations (12), (13), and (14), we can obtain:

[0085]

[0086] In formula (15):

[0087]

[0088] Then h(f) p ,f p′ This is called the point spread function (PSF), which characterizes the frequency component f. p′ The component at the point affects beamforming f p The contribution of the output. Ideally, the PSF should present a main lobe at the target frequency, while exhibiting side lobes due to energy attenuation at non-target frequencies. By optimizing the side lobe suppression capability of the PSF, the aliasing frequency and the target frequency can be effectively distinguished, thereby improving the robustness of intrinsic frequency identification.

[0089] Perform eigenvalue decomposition on the cross-correlation matrix C:

[0090]

[0091] In equation (17), U represents the unitary matrix, and σ represents the singular value of the correlation matrix. Introducing the exponent ξ, C is calculated. 1 / ξ as follows:

[0092]

[0093] According to formula (18), the output of beamforming optimization by introducing an exponential term is the following function: Beamforming output:

[0094]

[0095] In equation (19), ξ is an exponent greater than 1.

[0096] Assume that the signal contains only the frequency component f. p The following relationship holds:

[0097] σ2=σ3=…=σ M(N-K+1) =0, (20)

[0098]

[0099] In equation (21), tr(·) represents the trace of the matrix; in equation (22), u1 represents the first column of matrix U, i.e., the eigenvector corresponding to the first singular value. Combining equations (18) to (22), we can obtain:

[0100]

[0101] According to the Cauchy-Buniakowsky-Schwarz inequality, the point propagation function should satisfy:

[0102]

[0103] If and only if f p =f p′ The equality holds. From equation (24), it can be seen that at the target frequency, the value of the point propagation function is close to 1, while at non-target frequencies, the value of the point propagation function is between 0 and 1. From equation (23), it can be seen that when ξ is greater than 1, the output at non-target frequencies can be effectively suppressed, while the output at the target frequency remains basically unchanged. By suppressing the output at non-target frequencies, function beamforming can suppress frequency aliasing, effectively suppress the output at aliasing frequencies, and thus identify the blade's natural frequency.

[0104] A flowchart for identifying the natural frequency of a rotating blade using function beamforming is shown below. Figure 3 As shown. Figure 3As shown, the natural frequency of the blade is identified through function beamforming without requiring prior knowledge of the number of harmonics. The following section discusses... Figure 3 The steps of the spectrum analysis method for undersampled harmonic signals of the present invention will be described in detail.

[0105] S1. Acquire multi-channel leaf tip timing signal

[0106] In this step, existing technology is used to obtain multi-channel blade tip timing signals for rotating blades.

[0107] S2, Construct a multi-block snapshot matrix

[0108] In this step, the following is adopted: Figure 2 The method shown is used to construct a multi-block snapshot matrix.

[0109] S3. Calculate the cross-correlation matrix

[0110] In this step, the cross-correlation matrix C is calculated using the aforementioned equation (9).

[0111] S4. Calculate the power formula for the cross-correlation matrix.

[0112] In this step, the power form of the correlation matrix C is calculated using the aforementioned equation (18).

[0113] S5. Calculate the guide vector

[0114] In this step, the frequency f is calculated using the aforementioned equation (10). p The guide vector s p .

[0115] S6, Function Beamforming Output

[0116] In this step, the function beamforming output FB(f) is obtained using the aforementioned equation (19). p ).

[0117] The technical effects of the present invention are verified through numerical simulation experiments.

[0118] In the simulation experiment, the blade rotation speed was set to 78% of its maximum speed. Two natural frequencies were set for the blade, 45Hz and 220Hz, with amplitudes of 0.01mm and 0.02mm respectively, and initial phases of both were set to 0°. Six BTT sensors were installed circumferentially, with circumferential angles of 16°, 159°, 177°, 181°, 240°, and 276° respectively.

[0119] The frequency range identified by the algorithm is 0–300 Hz, with a frequency interval of 0.1 Hz. The proposed function beamforming method is compared with the block MUSIC method and conventional beamforming methods to verify the effectiveness and superiority of the method of the present invention. In the simulation, the number of blade rotations is set to 209, the number of blocks is 30, and the number of snapshots is 180, i.e., the dimension of the constructed multi-block snapshot matrix is ​​180×180, and a rectangular window function is used.

[0120] In the simulation experiment, the signal-to-noise ratio (SNR) range was set from -5dB to 50dB, with an SNR interval of 5dB. The function beamforming method ξ was set to 30. To evaluate the accuracy of different algorithms in identifying the natural frequency of the blade, 100 Monte Carlo simulations were performed at each SNR, and the recognition success rate χ% was defined as:

[0121]

[0122] In equation (25), successful recognition is defined as the frequency deviation corresponding to the two highest energy peaks in the output of each algorithm being within 1 Hz. The calculated recognition success rates of each method are as follows: Figure 4 As shown in the figure. Since the MUSIC method requires an estimate of the number of frequencies, MUSIC-1 in the figure caption represents an estimated frequency count of 1, and so on. Through... Figure 4 It is evident that the success rate of the MUSIC method depends on the estimated number of frequencies. When the estimated number of frequencies is 1, the MUSIC method fails to identify the two natural frequencies of the blade. For the MUSIC method, at low signal-to-noise ratios, inaccurate noise subspace estimation when selecting 10 frequency components for reconstruction results in a lower success rate compared to selecting only 2 or 5 frequencies. The function beamforming method, on the other hand, does not require prior knowledge of the number of frequencies and achieves a 100% success rate at a signal-to-noise ratio of 5 dB. Conventional beamforming methods, however, fail to identify the natural frequencies of the blade at all.

[0123] The recognition results of each method were visualized when the signal-to-noise ratio was selected as 5dB and 30dB, respectively. For the MUSIC method, only the results of selecting 1 frequency and 5 frequencies, i.e., MUSIC-1 and MUSIC-5, are displayed. The visualization results are as follows: Figure 5As shown in the visualization, the MUSIC-1 method, because it only retains one frequency component, cannot identify the 45Hz frequency component, resulting in a 0% recognition success rate. The conventional beamforming method, due to severe frequency aliasing, exhibits a significant increase in energy at 62Hz, making it unable to identify the 45Hz frequency component, thus also resulting in a 0% recognition success rate. At a signal-to-noise ratio of 5dB, the function beamforming result produces some sidelobes, but their amplitude is lower than the target frequency and does not affect the target frequency recognition effect. Therefore, it can be concluded that the function beamforming method can accurately identify the target frequency even at low signal-to-noise ratios without requiring prior information on the number of frequencies.

[0124] The performance of the function beamforming method depends to some extent on the setting of the power ξ. Therefore, the impact of different ξ settings on the frequency identification performance of the function beamforming method is further analyzed.

[0125] The signal-to-noise ratios were selected as 5dB and 30dB, and the ξ range was 10. 0 ~10 4 Twenty points are evenly distributed in the middle. Similarly, 100 Monte Carlo simulations are performed at each signal-to-noise ratio to calculate the success rate of frequency identification. The results are as follows: Figure 6 As shown, when ξ is 1, its output is equivalent to the traditional beamforming method. It can be observed that as ξ increases, the recognition success rate of the described function beamforming method increases, up to 100%. At a signal-to-noise ratio of 5dB, setting ξ to around 20 achieves a 100% recognition success rate. However, at higher signal-to-noise ratios, a smaller value for ξ achieves a 100% recognition success rate.

[0126] The frequency identification results are visualized under different ξ conditions. Figure 7 . Figure 7 The image shows the frequency identification results with ξ of 2, 5, 100, and 1000 at a signal-to-noise ratio of 5 dB. Figure 7 As shown, when ξ is small, the function beamforming method has a poor effect on suppressing aliasing. When ξ is 2, the target frequency component still cannot be identified. When ξ is 5, sidelobes of a magnitude comparable to the target frequency are generated, interfering with the frequency identification results. With further increases in ξ, when ξ is 100 and 10000, frequency aliasing is effectively suppressed, and the target frequency can be identified.

[0127] To evaluate the algorithm efficiency under different ξ settings, the average time consumed by 100 Monte Carlo simulations was calculated, and the results are as follows: Figure 8As shown, the signal-to-noise ratio (SNR) has no significant impact on the algorithm's runtime. Furthermore, the runtime of the algorithm does not exhibit a significant correlation with the value of ξ. Under different conditions, the described function beamforming method can complete one computation in an average of 0.25 seconds. The difference between the maximum and minimum average time required by the method of this invention does not exceed 0.15 seconds.

[0128] The accuracy and effectiveness of the method of this invention in identifying the natural frequency of rotating blades were verified using a rotating bladed disk experimental setup. The overall layout of the experimental setup is as follows: Figure 9 As shown in the figure. The circumferential arrangement of the BTT sensors is consistent with the numerical simulation, and the blade rotation speed is set to 78% of its maximum speed. A set of magnets is placed on the blade's rotation trajectory to provide an external excitation force to excite the blade's natural frequency components. In the experiment, the number of snapshots is set to 180, the number of blade rotations is 209 (corresponding to 209 rotations of the bladed disk), the number of blocks is 30, and the constructed multi-block snapshot matrix has a dimension of 180×180, using a rectangular window function consistent with the parameters in the numerical simulation.

[0129] Two pairs of blades and laser displacement sensors were installed on the impeller to ensure the dynamic balance of the overall rotor system. The laser displacement sensors were used to collect the blade displacement information during rotation, which was then compared with the blade's natural frequency identified via BTT signals. The sampling frequency was 3000 Hz. Data from only one blade was collected for analysis during the experiment.

[0130] To determine the natural frequency of the blade, a frequency sweep experiment was first conducted. The sweep conditions were as follows: initially maintaining 67% of the maximum rotational speed, then increasing the speed to 83% of the maximum speed after signal acquisition began, and then maintaining 83% of the maximum speed. The time spectrum and frequency spectrum of the signal acquired by the laser displacement sensor during the above process are shown below. Figure 10 As shown, a frequency resolution of 0.1 Hz was chosen when calculating the spectrum. The time-frequency graph shows two consistent frequency components, which are the two natural frequencies of the blade. Two distinct resonance peaks are visible in the spectrum, with frequencies of approximately 45 Hz and 219 Hz, representing the approximate frequency range of the two natural frequencies of the rotating blade.

[0131] exist Figure 11 The image shows the spectrum of the laser displacement sensor at 78% of its maximum speed, with a frequency resolution of 0.1 Hz. From this spectrum, the natural frequencies of the blade can be determined to be 47 Hz and 219.3 Hz.

[0132] Next, using different natural frequency identification methods, the natural frequency of the blade was identified by the BTT signal. The block size for both the MUSIC and function beamforming methods was set to 30, and the number of snapshots was set to 180, which means a total of 209 data cycles were required. Figure 12 This paper demonstrates the results of using the block MUSIC algorithm to identify the natural frequencies of a blade under different frequencies. The results show that the frequency values ​​significantly affect the MUSIC method's performance. When the number of frequencies is 1, 2, or 5, the block MUSIC method fails to accurately identify the blade's natural frequencies. Only when 10 frequencies are selected can the block MUSIC method identify the two natural frequencies of the blade, but there is an absolute error of 0.2 Hz for the first natural frequency. Because the blade displacement signal obtained from actual tests has complex components, selecting the number of frequencies for the block MUSIC algorithm is relatively difficult.

[0133] Figure 13 The results of conventional beamforming and function beamforming with different ξ settings for identifying the natural frequencies of rotating blades are shown. Figure 13 As shown, conventional beamforming algorithms suffer from severe frequency aliasing, failing to identify the target's natural frequencies. Increasing the power ξ to 2 improves the aliasing somewhat, but still fails to identify the natural frequencies. When ξ is 5, two natural frequency peaks are clearly observed. Further increasing ξ to 100 and 10000 further suppresses the aliasing, allowing for better identification of the two natural frequencies. At ξ values ​​of 100 and 1000, the absolute error for identifying the two modal frequencies is 0 Hz. Furthermore, the described function beamforming method does not require pre-estimation of the number of frequencies; simply setting ξ to a large value suffices. Simulation results show that the value of ξ has no significant impact on algorithm efficiency. Therefore, in practical applications, a value greater than 100 can be chosen for ξ.

[0134] It is worth noting that the above description is only a preferred embodiment of the present invention and does not limit the scope of patent protection of the present invention. The present invention can also be replaced by equivalent technologies. Therefore, all equivalent changes made based on the description and figures of the present invention, or direct or indirect applications to other related technical fields, are included within the scope of the present invention.

Claims

1. A method for spectral analysis of undersampled harmonic signals, characterized in that: Includes the following steps: S1. Acquire multi-channel leaf tip timing signals; S2. Construct a multi-block snapshot matrix based on the multi-channel leaf tip timing signal: ； in, Let K be the number of snapshots in the matrix. This is the transfer matrix from frequency to tip displacement. Here is the frequency pseudo-amplitude matrix to be solved. It is Gaussian white noise; S3. Calculate the cross-correlation matrix C: ； S4. Calculate the following power expression for the cross-correlation matrix: ； Where U is a unitary matrix, Let ξ be the singular values ​​of the matrix, and ξ be the exponent greater than 1. S5, Calculate Frequency The guide vector : ； in, … This is the delay factor; S6, Function Beamforming Output : ； The blade tip displacement signal obtained by the m-th sensor in step S1 during the l-th revolution. for: ； Where T is the sampling period and P is the number of frequency components of the blade. Let be the coefficient of the p-th frequency component. For the corresponding phase, ； The highest frequency expected to be analyzed. This is the delay factor; Step S2 further includes: dividing matrix A into blocks, with M rows forming one block, which can be divided into... Block, the i-th block for: ； Matrix X is: 。 2. The spectral analysis method for undersampled harmonic signals according to claim 1, characterized in that: In step S4, ξ is greater than or equal to 100.

Images