Method for monitoring abnormality of key auxiliary equipment of thermal power unit based on vibration data

Abstract

The application relates to the technical field of data processing, and discloses a thermal power unit key auxiliary equipment abnormality monitoring method based on vibration data, which comprises the following steps: collecting original vibration acceleration time sequence signals of key auxiliary equipment of a thermal power unit, classifying device states, and constructing a training sample set; constructing an angle domain impact feature of the vibration signal of the auxiliary equipment of the thermal power unit without a tachometer; classifying auxiliary equipment abnormality based on a order-phase double-channel special-shaped convolution network; training an abnormality classification model and optimizing parameters; and online abnormality identification of the key auxiliary equipment of the thermal power unit. The application has the beneficial effect that an instantaneous rotational frequency estimation method of multi-harmonic energy barycenter weighted fusion is proposed, a plurality of low-order harmonics are synchronously tracked, and the harmonics are fused by taking the frequency band energy as a weight, so that the problem that a single fundamental frequency spectrum peak is prone to jumping and breaking under noise interference is solved.

Description

Technical Field

[0001] This invention relates to the field of data processing technology, specifically to a method for monitoring anomalies in key auxiliary equipment of thermal power units based on vibration data. Background Technology

[0002] Key auxiliary equipment in thermal power units, such as forced draft fans, induced draft fans, feedwater pumps, and condensate pumps, are crucial for ensuring the stable operation of the power generation process. Their health status directly affects the reliability and economy of the entire unit. Anomaly monitoring of these devices, especially online monitoring based on vibration signals, has become an indispensable technical means in industrial operation and maintenance. However, in actual engineering scenarios, existing monitoring methods generally suffer from the following problems: First, instantaneous frequency estimation without a tachometer often relies on single-peak tracking. When the fundamental frequency energy is locally weakened due to fluid disturbances or structural resonance, the frequency trajectory is prone to jumps, breaks, or misjudgments, resulting in poor estimation continuity. Second, fault impact extraction methods typically use fixed waveform templates or single-band filtering, lacking adaptability to differences in impact morphology and failing to simultaneously and effectively retain the characteristics of both broadband transient impacts and narrowband damped oscillations—two types of heterogeneous faults. Third, existing vibration signal image representations often use general time-frequency diagrams, failing to decouple and organize the two types of cyclic stationary information—order distribution and phase distribution—resulting in high feature dimensionality and unclear physical meaning. Fourth, conventional convolutional networks use uniform square convolutional kernels to process feature maps, failing to distinguish between the periodic continuous properties of the phase axis and the discrete physical properties of the order axis. This easily leads to the mixing and modeling of two unrelated types of features, reducing classification interpretability and introducing noise interference. Summary of the Invention

[0003] To address the shortcomings of existing technologies, this invention provides a method for abnormal monitoring of key auxiliary equipment in thermal power units based on vibration data, thereby solving the problems mentioned in the background art.

[0004] To achieve the above objectives, the present invention employs the following technical solution: A method for monitoring anomalies in key auxiliary equipment of thermal power units based on vibration data includes the following steps: S1. Collect the original vibration acceleration time-series signals of key auxiliary equipment of thermal power units, and label the equipment status to construct a training sample set; S2. Perform angular domain transformation and impact feature extraction on the original vibration acceleration time series signal without tachometer to obtain the angular domain impact feature signal; convert the angular domain impact feature signal into an impact energy map with order-phase joint distribution; S3. Construct an order-phase dual-channel irregular convolutional network. Input the impact energy map into the order-phase dual-channel irregular convolutional network for feature extraction and fusion, and output the device status anomaly classification result. The order-phase dual-channel irregular convolutional network includes: order physical location encoding injection module, parallel phase convolution branch and order convolution branch, order gated adaptive fusion module, angle global cumulative pooling layer and classification layer. S4. Using the training sample set, the parameters of the order-phase dual-channel irregular convolutional network are iteratively optimized by minimizing the loss function to obtain the anomaly classification model. S5. Collect the real-time raw vibration acceleration time-series signal of the device to be identified, construct the impact energy map to be identified according to step S2, input the anomaly classification model, and output the device status anomaly identification result.

[0005] Furthermore, the vibration acceleration sensor is installed at any of the following locations on the wind turbine equipment: blower, induced draft fan, feedwater pump, condensate pump, circulating water pump, and other auxiliary equipment with rotating parts; the label includes at least one of the following: normal condition, abnormal rubbing, local bearing damage, abnormal loose connection, and abnormal local impeller impact.

[0006] Furthermore, S2 specifically refers to: S21. Based on the instantaneous frequency shift estimation of the energy centroid of the multi-harmonic frequency band, a smooth instantaneous frequency shift estimation sequence is obtained; S22. Based on the smooth instantaneous frequency conversion estimation sequence, the original vibration acceleration time series signal is resampled at equal angles to obtain the angular domain vibration signal; S23. Construct low-quality factor atomic libraries and high-quality factor atomic libraries, and use a matching pursuit strategy to sparsely reconstruct the angular domain vibration signal to obtain the angular domain impact characteristic signal. S24. Convert the angular domain impact characteristic signal into an impact energy map with a combined order-phase distribution.

[0007] Furthermore, S21 further includes: A short-time Fourier transform is performed on the original vibration acceleration time-series signal to obtain the time spectrum. Within a preset fundamental frequency search range, a coarse fundamental frequency estimate for each analysis frame is determined based on the concentration of spectral energy. Based on the coarse fundamental frequency estimate, harmonic search intervals are established for multiple low-order harmonics. Local band energy and local energy centroid frequency are calculated within each harmonic search interval. The local energy centroid frequencies of each harmonic are converted back to the fundamental frequency, and weighted fusion is performed using the corresponding local band energy as weights to obtain the instantaneous frequency shift estimate. Anomaly frame correction and smoothing processing are performed on the instantaneous frequency shift estimate sequence to obtain a smoothed instantaneous frequency shift estimate sequence.

[0008] Furthermore, S22 further includes: Based on the smooth instantaneous frequency shift estimation sequence, the frame-level instantaneous frequency shift is interpolated to the original sampling point position to obtain the sampling point-level instantaneous frequency shift sequence; Calculate the cumulative rotation angle sequence based on the instantaneous frequency conversion sequence at the sampling point level; Set the number of sampling points per rotation angle and establish an equal-angle sampling grid; use the cumulative rotation angle sequence as the independent variable and the original vibration signal as the dependent variable for interpolation and resampling to obtain the angular vibration signal.

[0009] Furthermore, S23 further includes: The angular domain vibration signal is divided into multiple overlapping angular domain blocks; real-valued decaying oscillation atoms are constructed, and low-quality factor atom libraries and high-quality factor atom libraries are constructed according to the range of decay coefficient values; for each angular domain block, the atom with the greatest correlation to the current residual signal is selected from the low-quality factor atom library and the high-quality factor atom library as the optimal atom, the projection coefficient is calculated and the residual signal is updated. Repeat the above atom selection and residual update process until the stopping condition is met; The selected atomic projection components within each corner domain block are accumulated and fused through overlapping regions to obtain the corner domain impact feature signal.

[0010] Furthermore, S24 further includes: Based on the number of sampling points per rotation angle, the angular domain impact characteristic signal is divided into multiple complete rotation cycles; the phase range of a single rotation cycle is equally divided into multiple phase segments; a local angular domain analysis window is extracted at the center of each phase segment in each rotation cycle, and a discrete Fourier transform is performed to obtain the local order spectrum; the preset order range is divided into multiple order frequency bands, and the average spectral energy in each phase segment and each order frequency band is statistically analyzed to construct an impact energy map with a size equal to the number of phase segments multiplied by the number of order frequency bands.

[0011] Furthermore, the order physical location encoding injection module in step S3 is further used for: The impact energy map is expanded into a multi-channel energy projection tensor by 1×1 convolution; a trainable order position coding matrix is ​​defined, the size of which is the number of order frequency bands multiplied by the number of channels; the energy projection tensor and the order position coding matrix are added element-wise along the order dimension to obtain the input tensor with injected order physical position coding.

[0012] Furthermore, the parallel phase convolution branch and the order convolution branch in step S3 further include: the phase convolution branch uses a convolution kernel along the phase direction and adopts a cyclic padding method to output a phase feature map; the order convolution branch uses a convolution kernel along the order direction and adopts a zero padding method to output an order feature map; the phase feature map and the order feature map have the same size.

[0013] Furthermore, the order-gated adaptive fusion module in step S3 is further used to: perform global average pooling on the order feature map along the phase dimension to obtain the order summary matrix; input the order summary matrix into the gating network to generate gating coefficients; broadcast the gating coefficients along the phase dimension, multiply them element-wise with the phase feature map, and then add them element-wise with the order feature map to obtain the fused feature map.

[0014] Furthermore, the angle-based global cumulative pooling layer in step S3 is further used to: simultaneously perform max pooling and average pooling on the fused feature map along the phase dimension; combine the max pooling result and the average pooling result according to a preset ratio to obtain an order-channel compressed feature matrix; flatten the order-channel compressed feature matrix and input it into the classification layer, and output the device state probability vector through the Softmax function.

[0015] Compared with the prior art, the beneficial effects of the present invention are as follows: 1. This invention proposes an instantaneous frequency estimation method based on the weighted fusion of multiple harmonic energy centroids. By synchronously tracking multiple low-order harmonics and fusing them with their frequency band energy as weights, it solves the problem that a single fundamental frequency peak is prone to jumps and breaks under noise interference.

[0016] 2. This invention constructs a dual quality factor composite dictionary to achieve impact matching tracking sparse reconstruction, allowing two types of atoms, which are good at expressing broadband transient impacts and narrowband damped oscillations respectively, to compete locally, enabling adaptive matching of impact modes of different faults such as rubbing and bearing damage.

[0017] 3. This invention transforms the angular domain impact characteristic signal into an order-phase joint impact energy map, simultaneously encoding two key diagnostic information types—the rotational phase distribution of the fault impact and the energy accumulation in the order frequency band—in a fixed-size two-dimensional image.

[0018] 4. This invention constructs a dual-channel heterogeneous convolutional network of order and phase, and uses a cyclically filled phase convolutional branch and a zero-filled order convolutional branch to model two types of physical characteristics respectively. It also uses an order gating mechanism to adaptively fuse the two types of physical characteristics, thus solving the problem of heterogeneous features being modeled by conventional square convolutional kernels. Attached Figure Description

[0019] Figure 1 This is a flowchart of the present invention; Figure 2 It is the time spectrum obtained by performing a short-time Fourier transform on the original vibration signal; Figure 3 It is a comparison chart of the actual instantaneous frequency conversion, the instantaneous frequency conversion estimated by the center of gravity of the multi-harmonic energy, and the frequency conversion trajectory after abnormal frame correction and smoothing; Figure 4 It is an angular domain vibration signal diagram obtained by resampling the original vibration signal in the angular domain based on the smooth instantaneous frequency shift. Detailed Implementation

[0020] The present invention will be further illustrated below with reference to specific embodiments. It should be understood that these embodiments are for illustrative purposes only and are not intended to limit the scope of the invention. Furthermore, it should be understood that after reading the teachings of this invention, those skilled in the art can make various alterations or modifications to the invention, and these equivalent forms also fall within the scope defined in this application.

[0021] like Figure 1 As shown, a method for abnormal monitoring of key auxiliary equipment in thermal power units based on vibration data is proposed, the main contents of which are as follows: S1. Collect the original vibration acceleration time-series signals of key auxiliary equipment of thermal power units, and label the equipment status to construct a training sample set; This invention first collects vibration monitoring data of key auxiliary equipment of thermal power units. The objects of collection may include blowers, induced draft fans, feedwater pumps, condensate pumps, circulating water pumps and other auxiliary equipment with rotating parts. Since some operating equipment is not equipped with independent speed sensors or key phase signal channels, only acceleration sensors installed on bearing seats, housings or support structures are used to obtain vibration acceleration signals during the collection stage.

[0022] In one implementation, each auxiliary machine is equipped with at least one radial vibration acceleration sensor; when improved spatial orientation diagnostic capabilities are required, sensors can also be installed in the horizontal, vertical, and axial directions. The sampling frequency is denoted as... , This characterizes the number of vibration sampling points collected per unit time, for example, 25600 Hz. The sampling duration of a single training sample is denoted as... For example, if 4 seconds is chosen, then each training sample contains One vibration sampling point.

[0023] The first The original vibration acceleration time-series signal corresponding to each training sample is defined as follows: Original vibration acceleration time sequence signal This refers to discrete vibration signals directly acquired by an accelerometer and obtained through analog-to-digital conversion. It is the original vibration acceleration time sequence signal At the sampling point A single value at, where, Indicates the sample index. This represents the index of the original time sampling point.

[0024] The data collection process covers various operating conditions such as equipment startup, shutdown, speed increase, speed decrease, load adjustment, and stable operation, so that the training samples fully include vibration change patterns under variable speed operating conditions. For abnormal samples, the samples are categorized by combining maintenance records, fault reports, on-site inspection results, and expert interpretation results.

[0025] In a convenient implementation method, the number of device status categories is denoted as... For example, five equipment states can be selected, namely normal state, abnormal contact and rubbing, local damage to bearing, abnormal loose connection, and abnormal local impact of impeller. Each training sample is configured with a unique real category label, which can be recorded in the form of one-hot encoding.

[0026] S2. Construction of angular domain impact characteristics of vibration signals of auxiliary equipment of thermal power units without tachometer: The original vibration acceleration time series signal is transformed and impact characteristics are extracted in the angular domain without tachometer to obtain the angular domain impact characteristic signal; The angular domain impact characteristic signal is converted into an impact energy map with order-phase joint distribution; During the start-up, shutdown, and load change processes of key auxiliary equipment in thermal power units, the speed usually changes continuously. The original vibration acceleration time sequence signal has obvious non-stationary frequency modulation characteristics. The interval of the occurrence of fault impact on the time axis will be stretched or compressed as the speed changes. If fixed window spectrum analysis or conventional convolution modeling is directly performed on the time domain signal, it is easy to cause impact energy diffusion, ambiguity of fault phase relationship and unstable characteristics.

[0027] This invention first estimates the instantaneous rotation frequency using the energy centroid of the multi-harmonic band, then constructs a cumulative rotation angle sequence based on the instantaneous rotation frequency and performs equal-angle resampling to obtain the angular domain vibration signal; next, it reconstructs the fault-dominated impact component through dual-quality-factor composite dictionary matching and tracking; finally, it converts the angular domain impact characteristic signal into an impact energy map with a joint order-phase distribution. The specific steps are as follows: S21. Instantaneous frequency conversion estimation based on the energy centroid of multiple harmonic bands. The fundamental frequency and its low-order harmonics in the vibration signal of auxiliary equipment will move synchronously with the rotational speed of the equipment. However, in the actual operating environment, the fundamental frequency may be weakened locally due to the influence of fluid disturbance, structural resonance and random noise. If the rotational speed is estimated based on only a single frequency peak, it is easy to cause frequency trajectory jumps, breaks or misjudgments.

[0028] This invention utilizes the synchronous variation relationship among multiple low-order harmonics to obtain an instantaneous frequency shift estimation sequence through the fusion of energy centroids in multiple harmonic frequency bands. The specific steps are as follows: 1> The original vibration acceleration time sequence signal Perform a short-time Fourier transform to obtain the time spectrum. , Time spectrum At point A single numerical value, where, Indicates the analysis frame index, Indicates frequency index.

[0029] In one implementation, the analysis window length of the short-time Fourier transform can be 2048 sampling points, the window function can be a Hanning window, and the overlap rate between adjacent analysis windows can be 75%. This setting can achieve a good balance between the speed of rotational change and the frequency resolution capability.

[0030] 2> Locate a coarse base frequency within the preset base frequency search range. The preset base frequency search range can be set according to the auxiliary machine's rated speed and operating range, for example, it can be from 5 Hz to 100 Hz; Within each analysis frame, the location with relatively concentrated spectral energy within that frequency range is selected as the rough fundamental frequency estimate. Rough fundamental frequency estimation It refers to the first Within each analysis frame, within the preset fundamental frequency search range, the estimated value of the rotating fundamental frequency is initially determined based on the degree of spectral energy concentration.

[0031] In one embodiment, for example, the preset baseband search range is 5Hz to 100Hz; in the... Within a single analysis frame, if multiple spectral peaks exist in the range of 5Hz to 100Hz, but the spectral energy is most concentrated near 15Hz, then... The frequency was set to 15Hz.

[0032] In one implementation, to reduce interference from isolated noise peaks, the spectral energy within the search range can be locally averaged first, and then the frequency position with the largest local average energy can be selected.

[0033] 3> Based on rough fundamental frequency estimation Establish multiple harmonic search intervals.

[0034] Specifically, the highest harmonic order involved in the instantaneous frequency conversion estimation is denoted as... For example, it can be 3 or 5; for the first The first harmonic has a theoretical center frequency of . Expanding to the left and right of the center frequency Hertz forms the harmonic search range. For example, 2 Hz can be used.

[0035] In one embodiment, as an example, suppose , , Then according to the first The center frequency of the first harmonic theory and left and right expansion Hertz's rules establish the following search intervals for each harmonic: First harmonic search interval: ; Second harmonic search interval: ; Third harmonic search interval: .

[0036] 4> Calculate the local frequency band energy and local energy centroid frequency within each harmonic search interval. Specifically, for the th harmonic search interval... The analysis frame, the first First harmonic, local frequency band energy denoted as The local energy centroid frequency is denoted as ,in, It refers to the first Local frequency band energy within the first harmonic search interval It refers to the first The local energy centroid frequency corresponding to the first harmonic.

[0037] In practical implementation, By examining the first The spectral energy within the first harmonic search interval is obtained by summing the spectral energies. The local energy centroid frequency is obtained by weighting the corresponding spectral energy as the weights and using each frequency point as the location. It can characterize the true concentrated frequency of the harmonic in the current analysis frame more stably than a single peak location.

[0038] In one embodiment, for example, for the first Frame, second harmonic interval Assuming that there are only [cases] within this interval and The spectral energies at the two frequency points are respectively and Therefore, the local frequency band energy Local energy centroid frequency .

[0039] 5> Convert the local energy centroid frequencies of each harmonic back to the fundamental frequency, and then perform weighted fusion based on the corresponding local frequency band energy to obtain the instantaneous frequency shift estimate. , It refers to the first The instantaneous frequency estimate corresponding to each analysis frame, in Hertz, is calculated as follows: ; in, The characterization is by the first The fundamental frequency estimate obtained by harmonic distortion; This indicates the prevention of extremely small positive numbers with a denominator of 0, for example, we can take... .

[0040] 6> Estimation sequence of instantaneous frequency conversion Abnormal frame correction and smoothing are performed. Specifically, when the energy of all harmonic frequency bands involved in the estimation within a certain analysis frame is lower than a preset threshold, the instantaneous frequency shift estimation result of that analysis frame is deemed unreliable. The estimated value of that analysis frame is not directly used; instead, linear interpolation or spline interpolation is performed based on the preceding and following valid analysis frames. Then, the corrected instantaneous frequency shift estimation sequence is smoothed by a median filter or Savitzky-Golay smoothing filter with a length of 5 frames to obtain a smoothed instantaneous frequency shift estimation sequence. .

[0041] In one embodiment, for example, suppose three consecutive analysis frames The frequencies are 29.8Hz, 31.2Hz, and 30.5Hz respectively; if the energy of all harmonic frequency bands in the second frame... If all values ​​are below a preset threshold, then the frame... Invalid; at this point, linear interpolation is performed on the first frame (29.8Hz) and the third frame (30.5Hz) to obtain the corrected estimated value of 30.15Hz for the second frame.

[0042] It should be noted that this invention estimates the instantaneous frequency shift by using multi-harmonic redundancy fusion rather than single-peak tracking. A single frequency peak is susceptible to local noise interference, causing jumps, while the probability of multiple low-order harmonics being simultaneously submerged or weakened by noise is much lower than that of a single fundamental frequency. By converting the energy centroid frequencies of each harmonic back to the fundamental frequency and using their local frequency band energy as weights for weighted fusion, even if one harmonic is interfered with, other harmonics can still provide accurate frequency information, thus enabling the frequency shift trajectory to have better continuity and noise resistance even without an independent tachometer.

[0043] In one embodiment, such as Figure 2 As shown, the short-time Fourier transform time spectrum (the time spectrum obtained by performing a short-time Fourier transform on the original vibration signal) is analyzed. The horizontal axis represents time (unit: seconds), the vertical axis represents frequency (unit: Hertz), and the color intensity represents the logarithmic amplitude (unit: dB). The experimental results reflect the non-stationary frequency modulation characteristics of the signal and the broadband transient characteristics of the fault impact.

[0044] S22, Angle Domain Resampling Based on Instantaneous Phase Integral Based on the smoothed instantaneous rotational frequency estimation sequence, the original vibration acceleration time series signal is converted from time coordinates to rotational angle coordinates, thereby realigning the periodic impacts that were originally stretched or compressed on the time axis to the angle axis, reducing the impact of rotational speed changes on fault feature extraction. The specific steps are as follows:

[0045] 1> Smooth the instantaneous frequency shift estimation sequence at the frame level Interpolate to the original vibration sampling point location to obtain the instantaneous frequency transfer sequence at the sampling point level. Sampling point level instantaneous frequency conversion sequence At the sampling point The value is denoted as , It refers to the first The instantaneous frequency estimation value corresponding to each original time sampling point.

[0046] In the specific implementation, the frame index is analyzed. The corresponding time is the independent variable, with As the dependent variable, construct a one-dimensional interpolation function; then, interpolate the time corresponding to each sampling point of the original vibration signal. Substituting into the interpolation function, we obtain the instantaneous frequency conversion sequence at the sampling point level. For example, if the smoothed frequency transition of the first frame (time 0.1s) is 15Hz, and the second frame (time 0.2s) is 15.1Hz, then for a sampling point located between the two frames with a timestamp of 0.15s, its frequency can be obtained through linear interpolation. It is 15.05Hz.

[0047] 2> Based on the instantaneous frequency conversion sequence at the sampling point level Calculate the cumulative rotation angle sequence The sequence at the sampling point The value is denoted as This establishes a mapping relationship from "time" to "angle," providing an accurate angle reference for subsequent equal-angle resampling, thereby transforming a non-stationary signal in the time domain into a stationary signal in the angle domain. This refers to the period up to the [number]th The cumulative rotation angle traversed by each original time sampling point, in radians, is calculated as follows: ; in, This represents the index of the sampling point during the integration and summation process.

[0048] 3> Set the number of sampling points per corner region , This represents the number of angle sampling points retained after each complete rotation cycle, for example, 1024. Further, through... Calculate the fixed angle sampling interval This allows for the establishment of a uniform sampling grid on the angle axis, ensuring a fixed number of sampling points per rotation cycle and aligning signal features of the same phase across different cycles. It represents the fixed angle difference between adjacent angle sampling points.

[0049] 4> Constructing an equal-angle sampling position sequence ,in By establishing a uniform, discrete sampling grid on the angle axis, the sequence of target positions to be interpolated is clearly defined; where, Indicates the index of the angular domain sampling point.

[0050] Furthermore, using the cumulative rotation angle sequence The value of As the independent variable, the original vibration acceleration time series signal is used. At the sampling point The value at the location An interpolation function is established for the dependent variable, and sampling positions are taken at each equal angle. The value is taken at the point to obtain the angular domain vibration signal. The signal is indexed at the angular domain sampling point. The value at that location is denoted as Angular domain vibration signal It refers to the vibration signal after equal-angle resampling based on the instantaneous frequency estimation result, and its independent variable is the rotation angle rather than time.

[0051] In the specific implementation, a cumulative rotation angle sequence is used. The value of The independent variable is the original vibration signal. At the sampling point The value at the location For the dependent variable, establish a cubic spline or linear interpolation function; then, query each isoangular sampling point on this interpolation function. The function value at that point yields the angular domain vibration signal. For example, if , Equal angle sampling points If it falls exactly between the two, then it can be obtained through linear interpolation. .

[0052] In one embodiment, such as Figure 3 As shown, the actual instantaneous frequency (black dashed line), the instantaneous frequency estimated by the multi-harmonic energy centroid (red scatter dots), and the frequency trajectory after anomalous frame correction and smoothing (blue solid line) are compared. The horizontal axis represents time (in seconds), and the vertical axis represents frequency (in Hertz). The fluctuations of the red scatter dots reflect the noise interference and local jumps present in the direct estimation, while the smoothed blue curve is continuous and consistent with the trend of the actual frequency, proving that the proposed multi-harmonic redundancy fusion and smoothing processing can obtain stable and continuous tachometer-less frequency estimation under the condition of no independent tachometer.

[0053] S23. Shallow Reconstruction Based on Impact Matching Tracking Using a Dual Quality Factor Composite Dictionary Angular vibration signal The time-scaling effect caused by speed variations has been mitigated, but it still includes periodic harmonics, fluid disturbances, mechanical background vibrations, and random noise. The impact patterns corresponding to different faults are not consistent: rubbing anomalies are closer to transient impacts with short duration and obvious broadband components; local bearing damage is closer to narrowband impacts accompanied by damped oscillations after excitation.

[0054] This invention constructs low-quality factor atom libraries and high-quality factor atom libraries, and employs a local competitive block matching and tracking strategy to obtain corner domain impact feature signals. The specific steps are as follows: 1> Angular domain vibration signal Divided into multiple overlapping corner domain blocks, each corner domain block contains Each corner domain sampling point, For example, 256 or 512 can be used; overlapping areas are set between adjacent corner domain blocks, and the overlap rate can be, for example, 50%.

[0055] In one embodiment, as an example, suppose a corner domain block contains There are 10 sampling points with an overlap rate of 50%. Therefore, the first block covers the index. Up to 255, the second block is covered. Up to 383, the third block is covered. Up to 511, and so on.

[0056] 2> Iterate through multiple sets of structural parameters to generate discretized atoms, constructing real-valued damped oscillating atoms to describe local impact waveforms. This provides a redundant waveform dictionary for the matching pursuit algorithm, capable of sparsely representing different types of impacts. A single atom can be represented as: ; in, Indicated by A candidate atom at the center position specifically refers to a corner sampling point. As the independent variable, by , , The three parameters together control its shape as a decaying cosine oscillation function; This represents the position of the atom's central corner region. Within each corner region block, all possible corner region positions are traversed with a step size of 1 sampling point. The angular domain oscillation frequency is determined based on the order range of the fault characteristics to be analyzed and the number of sampling points per revolution in the angular domain. , which is converted into angular domain oscillation frequency and taken at equal or logarithmic intervals within a preset range; This represents the attenuation coefficient. The larger the atom, the shorter its duration, making it more suitable for describing broadband transient shocks; The smaller the value, the longer the tail of the atomic oscillation, making it more suitable for describing narrowband damped oscillation shocks.

[0057] Furthermore, based on the range of attenuation coefficient values, low-quality factor and high-quality factor atomic libraries are constructed. The low-quality factor atomic library is used to describe broadband transient impacts with short duration and strong time-domain concentration, and the attenuation coefficient... The value can be taken from 0.08 to 0.20; a high-quality factor atom library is used to describe damped resonance impacts with long durations and obvious oscillation characteristics, and the damping coefficient... It can be between 0.01 and 0.05. In the actual implementation, the atoms in both atomic libraries are normalized to unit energy.

[0058] In one embodiment, for example, within a 512-point block, the low-quality-factor atomic library may include: , Take 20 discrete values, The values ​​were 0.08, 0.12, 0.16, and 0.20; therefore, this atomic library contains a total of Candidate atoms normalized to unit energy.

[0059] 3> Initialize the residual signal sequence for each corner domain block. The sequence is at index The value at that location is denoted as Its elements Equals the current corner domain block in the index The original values ​​are then calculated; then, the normalized correlation between all candidate atoms in the low-quality factor atom library and the high-quality factor atom library and the current residual signal is calculated, and the atom with the largest absolute value of correlation is selected as the optimal atom for the current iteration round.

[0060] In one implementation, if the candidate atom sequence is denoted as The current residual signal sequence is denoted as They are in the index The values ​​at each location are denoted as follows: and Then the first The optimal atom selected in the round of iteration satisfies ;in, Indicates the current residual signal and the first The inner product of candidate atoms, Indicates all candidate atom indices In the middle, find the one that makes the absolute value of the inner product equal to the value of the inner product. The index corresponding to the point where the maximum value is reached; Indicates the first The index of the optimal atom selected in the round of iteration.

[0061] In one embodiment, for example, suppose the current residual signal It is a block signal, candidate atoms It is a fast-decaying waveform centered at index 10. It is a slowly decaying waveform centered at index 50; calculate separately. and , The absolute value of the inner product, if , Then this round of selection The optimal atom.

[0062] 4> Calculate the projection coefficient corresponding to the optimal atom, and subtract the projection component of that atom from the current residual signal to obtain the residual signal for the next round; repeat the atom selection and residual update process until the stopping condition is met.

[0063] In one implementation, the stopping condition can be any of the following: a> The ratio of the current residual signal energy to the initial block energy is lower than the preset residual energy ratio threshold; b>The preset maximum number of iterations has been reached.

[0064] In this implementation, the residual energy ratio threshold can be set to 0.05, and the maximum number of iterations can be set to 20.

[0065] 5> The selected atomic projection components within each corner domain block are summed to obtain the local impact reconstruction signal for that block; for overlapping regions of adjacent blocks, windowed averaging or linear weighted averaging is used for fusion to obtain the complete corner domain impact feature signal. The signal is indexed at the angular domain sampling point. The value at that location is denoted as Angular domain impact characteristic signal It refers to the fault-dominated angular domain signal obtained after double quality factor composite dictionary matching and tracking. It mainly retains the transient impact and decaying oscillation structure related to the fault, and suppresses periodic harmonics, stationary background components and random noise.

[0066] In one embodiment, as an example, assume that the 50% overlap region length of blocks A and B is 128 points; within the overlap region, the reconstructed signal amplitude of block A is... Block B is When using linear weighted averaging fusion, the full weight at the starting point of the overlap region comes from... The final weight comes entirely from Midpoint The fusion value ,in For point The distance to the starting point of the overlapping area.

[0067] It should be noted that the dual-quality-factor composite dictionary used in this invention contains two types of atom libraries, each adept at expressing broadband transient impacts and narrowband damped oscillations, respectively. During signal reconstruction, a fixed waveform template is not forced; instead, the two types of atoms compete for selection within each local region of the signal. For broadband short-duration impacts caused by rubbing anomalies, lower-quality-factor atoms are more likely to be selected; while for narrowband damped oscillations induced by localized bearing damage, higher-quality-factor atoms are more advantageous.

[0068] In one embodiment, such as Figure 4 As shown, the angular domain vibration signal (equal-angle resampling) is analyzed. Specifically, it is the angular domain vibration signal obtained by resampling the original vibration signal in the angular domain based on the smoothed instantaneous rotation frequency. The horizontal axis represents the number of rotations (dimensionless, unit is revolutions), and the vertical axis represents the vibration acceleration amplitude. The impacts that were originally stretched or compressed on the time axis are now uniformly distributed, with each rotation cycle occupying a fixed number of sampling points, and the impact positions are aligned to the same rotation phase interval.

[0069] S24, Construction of Order-Phase Joint Impact Energy Map Angular domain impact characteristic signal This invention effectively preserves the fault-induced impact structure by converting the angular domain impact characteristic signal into an impact energy map with a joint order-phase distribution. The specific steps are as follows: 1> Based on the number of sampling points per corner region , angular domain impact characteristic signal Divided into multiple complete rotational cycles, the number of complete rotational cycles is denoted as . For angular sampling points that are less than one complete rotation cycle at the end, they can be discarded directly or filled with zero padding.

[0070] In one embodiment, for example, if the angular domain signal There are a total of 2600 sampling points, and the number of sampling points per revolution cycle is... Then two complete cycles can be divided ( , ), The last 552 sampling points ( If it is less than a complete cycle, discard it directly.

[0071] 2> Phase range of a single rotation cycle Divided into equal parts Each phase segment, For example, 64 can be taken; the... Each phase segment is used to aggregate local impact information located near the rotation phase, where, .

[0072] In one embodiment, for example, let the number of phase segments be... The 0th phase segment covers The first segment coverage And so on, the first Segmented coverage All angular domain signal points falling within these phase intervals will be aggregated.

[0073] 3> Within each rotation cycle, a local angular domain analysis window of fixed length is extracted around the center of each phase segment, and a discrete Fourier transform is performed on the local analysis window to obtain the local order spectrum. The spectrum is indexed at order frequency points. The value at that location is denoted as , It refers to the first One complete rotation cycle, the first Each phase segment corresponds to a local analysis window at the order frequency index. The spectral value at the order of . Where, Indicates the rotation period index; This indicates the order frequency point index.

[0074] In one embodiment, for example, for the 0th rotation cycle ( ), in the 0th phase segment ( )center( Near the signal, a local angular domain analysis window is extracted, and the discrete Fourier transform is performed on the signal within the window to obtain the local order spectrum. ; Corresponding to order 0, For the first order, and so on.

[0075] 4> Divide the preset order range into The first order frequency band, the first The set of order indices corresponding to each order frequency band is denoted as . ,in, For example, when focusing on the fault energy distribution in the range of 0 to 20, this range can be divided into 32 order frequency bands.

[0076] It should be noted that if the interval is directly set to 1 order, only 21 frequency bands can be obtained, resulting in a fixed resolution. By dividing the frequency bands into 32 bands, non-uniform or higher-density frequency band divisions can be introduced across the entire range of interested orders (e.g., using finer granularity in some intervals). This allows for more detailed analysis of fault-sensitive key order regions while compressing the order dimension to a fixed size suitable for the input layer of subsequent convolutional networks. This allows for the retention of diagnostic information while controlling the complexity of the feature maps.

[0077] 5> Calculate the average spectral energy within each phase segment and each order frequency band to obtain the impact energy map. The diagram is in the phase segment index. Order frequency band index The element value at position is denoted as This allows the high-dimensional, redundant local order spectrum information to be compressed into a fixed-size, low-dimensional, physically clear two-dimensional image that can simultaneously encode the energy distribution pattern of the fault impact in the two key physical dimensions of "phase" and "order". It refers to the first The sample at the th The first phase segment, the... The impact energy value at each order frequency band is calculated as follows: ; in, Indicates the first The number of effective rotation cycles actually involved in the statistics at each phase segment; Indicates the first The number of order frequency points contained within each order frequency band.

[0078] 6> All According to phase index and order frequency band index Arranged to form a size of Impact energy diagram Impact energy diagram It refers to the two-dimensional order-phase energy representation extracted from the angular domain impact characteristic signal. Its horizontal dimension represents the distribution change of the fault impact in the rotating phase, and its vertical dimension represents the degree of energy accumulation of the fault impact in different order frequency bands.

[0079] In one embodiment, for example, if a certain type of rubbing abnormality always occurs near... Near the phase, the impact energy map will show a significant local energy enhancement in the corresponding phase segment; if a certain type of bearing local damage causes a continuous impact near a specific order, the impact energy map will show an energy concentration region in the corresponding order frequency band. This two-dimensional representation can simultaneously retain two types of key diagnostic information: "near which phase the impact occurs" and "in which order frequency bands the impact is concentrated".

[0080] S3. Construct an order-phase dual-channel irregular convolutional network, input the impact energy map into the order-phase dual-channel irregular convolutional network for feature extraction and fusion, and output the equipment status anomaly classification results; Impact Energy Diagram It has a distinct two-dimensional heterogeneous structure. The phase axis reflects the location of the impact within the rotation cycle of the auxiliary machine and has periodic continuity. The order axis reflects the characteristic order distribution corresponding to different fault mechanisms and has discrete physical meaning. If a conventional square convolution kernel is used for unified processing, it is easy to mix the phase change mode and the local order structure in the model, which weakens the physical interpretability.

[0081] This invention designs an order-phase dual-channel irregular convolutional network, which is composed of the following modules connected in series: a> Order-based physical location coding injection: by It is formed by adding convolutional channel expansion and trainable positional encoding matrices; b> Dual-channel parallel feature extraction: includes a "phase convolution branch" with cyclic padding and a "order convolution branch" with zero padding and convolution along the order direction. c> Order-gated adaptive fusion: It consists of a global average pooling, a two-layer fully connected gating network, and element-wise multiplication / addition; d> Angle-based global cumulative pooling: Simultaneously performs max pooling and average pooling along the phase dimension, and then combines them with weights; e> Classification layer: Consists of a flattening operation and a fully connected layer connected to the Softmax function.

[0082] The specific steps are as follows: S31, Order of Physical Location Encoding Injection Different order frequency bands in the impact energy map have different physical meanings. If the order position information is not introduced, the network needs to rely entirely on the training data to learn "which vertical position corresponds to which order range", which can easily increase the training difficulty.

[0083] This invention explicitly injects order-based physical information through trainable position encoding, and the specific steps are as follows: 1> Impact energy diagram Expand to Channel energy projection tensor The tensor at position The element value at position is denoted as , specifically adopted Convolution expands the channels, thereby increasing the feature dimension and providing a multi-channel representation space for subsequent addition with positional encoding. This allows the network to learn diverse features related to order position across different channels, as shown below: ; in, Indicates the channel index. ; This indicates the total number of channels after expansion, for example, 16. Indicates the first Trainable scaling parameters corresponding to each projection channel; Indicates the first Trainable bias parameters corresponding to each projection channel.

[0084] 2> Define a trainable order positional encoding matrix Its size is The first in the matrix Line number Column elements Used to indicate the first The first order frequency band in the... Learnable physical location labels on each feature channel, order of location encoding matrix The size is .

[0085] Furthermore, the energy projection tensor With order position encoding matrix Element-wise addition (implemented via broadcasting) yields the input tensor. Furthermore, the absolute physical location information of the order is injected into the network input as prior knowledge, so that the network "knows" the order corresponding to each position in the vertical direction of the feature map from the beginning; input tensor This refers to the network input tensor with injected order physical location encoding, and its size is... ,definition The summation of the energy projection tensor and the position encoding matrix at corresponding positions is expressed as follows: .

[0086] It should be noted that the order physical location coding is set only along the order frequency band dimension, and the order axis has a clear absolute physical layering meaning. The physical order corresponding to each frequency band is fixed and ordered. Introducing position coding can strengthen this vertical discrete physical property. The phase axis represents a complete rotation cycle, which is connected end to end and has natural cyclic continuity. There is no absolute "start" or "end" position concept, so no fixed position coding is set.

[0087] S32. Parallel extraction of phase trend features and order local features The change in the impact energy map along the phase direction reflects the distribution pattern of the fault impact within the rotation period, and the change along the order direction reflects which order frequency bands the fault energy is concentrated in. This invention uses two parallel heterogeneous convolution channels to model these two types of physical characteristics respectively. The specific steps are as follows: 1> Input tensor The data is fed into the phase convolution branch, which uses a convolution kernel that expands along the phase direction. The kernel size is set to... ,in, This indicates the phase axis convolution length, for example, it can be 7.

[0088] In practical implementation, since the phase axis corresponds to a complete rotation cycle and has a head-to-tail relationship, the phase convolution branch preferably adopts a cyclic filling method, so that the 0th phase segment is connected to the 1st phase segment. Each phase segment maintains its adjacency during convolution processing; the phase convolution branch outputs a phase feature map. .

[0089] It should be noted that the phase convolution branch consists of two consecutive convolutional layers. The kernel size of the first convolutional layer is... With a step size of 1, using cyclic filling, the number of output channels is [number missing]. The second layer is followed by batch normalization and ReLU activation; the kernel size, padding method, and number of output channels are the same as the first layer; its output is a phase feature map. Its function is to extract local continuous change patterns and periodic impact trends along the rotation phase direction.

[0090] 2> Input tensor The input is fed into the order convolution branch, which uses a convolution kernel that expands along the order direction. The kernel size is set to... ;in, This indicates the convolution length of the order axis, for example, it can be 5.

[0091] In the actual implementation, the order axis does not have a beginning-end loop relationship, so the order convolution branch uses regular zero padding or boundary copy padding; the order convolution branch outputs the order feature map. .

[0092] It should be noted that the order convolution branch consists of two consecutive convolutional layers. The kernel size of the first convolutional layer is... With a step size of 1, zero-padding is used, and the number of output channels is [number missing]. The second layer is followed by batch normalization and ReLU activation; the kernel size, padding method, and number of output channels are the same as the first layer; its output is a feature map of order. Its function is to extract local energy accumulation and spectral peak abrupt structure along the order axis.

[0093] 3> Phase feature map With order feature map Set to the same size ,in, This represents the number of output channels, specifically determined by uniformly setting the number of output convolution kernels for both convolution branches. To achieve this, for example, both can be set to 32, so that the feature maps output by the two branches have the same channel dimension.

[0094] Among them, the phase feature map The key features are the continuous variation trend of the impact along the rotating phase axis; order characteristic diagram. The focus is on characterizing the local accumulation and abrupt changes of energy in adjacent order frequency bands.

[0095] S33, Order-gated adaptive fusion

[0096] Different order frequency bands do not necessarily need to introduce phase trend information to the same extent. If a certain order frequency band does not show obvious abnormal energy, the phase characteristics corresponding to that order frequency band may mainly come from noise.

[0097] To reduce interference from invalid phase modes, this invention utilizes order features to generate gating coefficients and performs conditional modulation on the phase features. The specific steps are as follows: 1> For order feature maps Along the phase dimension Perform global average pooling on the feature map at location The element value at position is denoted as The order summary matrix is ​​obtained. The matrix in the th The first order frequency band, the first The element values ​​at each channel are denoted as This further compresses the spatial information of the phase direction, resulting in a summary representation that is only related to the order and channel, expressed as: ; in, Indicates the first The sample at the th The first order frequency band, the first The overall order response intensity on each channel.

[0098] 2> Divide the order summary vector Input a lightweight gating network to generate gating vectors The vector at the th The gating coefficient for each channel is denoted as By utilizing the strength of the order features themselves, the inflow of phase trend features is adaptively controlled; gating coefficient The closer the value is to 1, the more phase trend characteristics need to be introduced into that frequency band; gating coefficient The closer it is to 0, the more likely the phase trend characteristics corresponding to that order frequency band should be suppressed.

[0099] In one implementation, the gated network consists of two fully connected layers: the weight matrix of the first fully connected layer of the gated network has a size of... The bias vector length is The first layer will change the channel dimension from Compress to ,in, This represents the compression ratio, for example, 4; then it's followed by the ReLU activation function; the size of the weight matrix of the second fully connected layer is... The bias vector length is The second layer restores the channel dimension to Finally, the Sigmoid activation function is applied to restrict the output to between 0 and 1.

[0100] 3> Gating coefficient Broadcast along the phase dimension and multiply element-wise with the phase feature map, then add it to the order feature map to obtain the fused feature map. , fusion feature map This refers to the comprehensive feature representation after order-gated adaptive fusion, which simultaneously includes phase dynamic trend information and order-local energy accumulation information adaptively filtered by the gated signal. The elements are denoted as... , It refers to the fusion of feature maps In phase index Order frequency band index and channel index The specific element value at a given location is used to uniquely identify a data point in this tensor, and its calculation method is expressed as follows: .

[0101] In one embodiment, if a certain order frequency band corresponds to the local damage characteristic order of the bearing and exhibits a strong abnormal response in the order convolution branch, the gating network will assign a higher gating coefficient to that order frequency band, so that the impact location change information extracted by the phase convolution branch is preserved; if the overall response of a certain order frequency band is weak, the corresponding gating coefficient will be reduced, thereby reducing the noise phase mode from entering the subsequent classification process.

[0102] S34, Angle-based global cumulative pooling and device state classification Since auxiliary machine failures can manifest as either strong localized shocks within a few phase intervals or as a weak, persistent abnormal trend throughout the entire rotation cycle, this invention employs a global cumulative pooling method combining max pooling and average pooling. The specific steps are as follows: 1> For fused feature maps Along the phase dimension Simultaneously performing max pooling and average pooling compresses the phase space dimension, thus... Feature maps are aggregated into The matrix, while preserving order information, integrates the impact performance of different phase intervals to prepare for final classification.

[0103] 2> Combine the results of max pooling and average pooling in proportion to obtain the order-channel compressed feature matrix. The matrix in the th The first order frequency band, the first The element values ​​at each channel are denoted as Then, the local strong impact characteristics (captured by max pooling) and the global persistence characteristics (captured by average pooling) are weighted and synthesized to obtain more robust device state characteristics, represented as: ; in, The combination coefficient between the max pooling component and the average pooling component is a hyperparameter between 0 and 1, used to adjust the contribution ratio between the max pooling result and the average pooling result. For example, it can be set to 0.7. The larger the size, the more it emphasizes the characteristics of strong local impact; The smaller the value, the more emphasis is placed on the average trend characteristics over the entire cycle.

[0104] 3> Compress the order-channel feature matrix Flattening them in a fixed order yields a one-dimensional classification feature vector. If the number of order frequency bands is The number of channels is ,but The length is .

[0105] Furthermore, the one-dimensional classification feature vector Input the fully connected classification layer and output the device state probability vector through the Softmax function. , is represented as: ;

[0106] in, This represents the Softmax function; Indicates the first The device state probability vector corresponding to each sample; Represents the classification layer weight matrix; This represents the classification layer bias vector; if the number of device state categories is... ,but for A probability vector of dimension.

[0107] S4. Anomaly Classification Model Training and Parameter Optimization After obtaining labeled training samples, the network parameters are iteratively optimized by minimizing the loss function. The specific steps are as follows: First, a training sample set is constructed; specifically, each original vibration acceleration time-series signal is used. The instantaneous frequency estimation, angle domain resampling, impact sparse reconstruction, and impact energy map construction are performed sequentially to obtain the corresponding impact energy map. .

[0108] Furthermore, set the training batch size. For example, 32 can be chosen; each training iteration selects from the training sample set. Each impact energy map and its corresponding label is input into the network.

[0109] Furthermore, class weights are set based on the number of samples of different device state categories in the training set. Specifically, the first... The number of samples in each category is denoted as The maximum number of samples in all categories is denoted as , No. The loss weights corresponding to each category are denoted as follows: ,definition ;in, The outlier class with fewer samples receives a larger loss weight, thus reducing the bias towards the majority class during training.

[0110] Further, forward propagation and loss calculation are performed. In one implementation, class-weighted cross-entropy loss is used as the objective function, and the Adam optimizer is used to update the network parameters. The initial learning rate can be set to... If the validation set loss does not decrease for several consecutive training epochs, the learning rate will be reduced by a preset ratio, such as multiplying by 0.5; the number of training epochs can be an integer between 100 and 200.

[0111] After each training cycle, the model performance is evaluated using a validation set. Evaluation metrics can include macro-average F1 score, accuracy, and recall. The network parameters corresponding to the highest macro-average F1 score are saved as the parameters for the final anomaly classification model. The final anomaly classification model refers to the model based on the impact energy map. As input, a device state probability vector This is the output auxiliary equipment anomaly identification model.

[0112] S5. Online anomaly identification of key auxiliary equipment in thermal power units After training the anomaly classification model, this invention can be deployed in the online condition monitoring system for key auxiliary equipment in thermal power units to identify anomalies without a tachometer in real-time vibration data. The specific steps are as follows: The system collects real-time raw vibration acceleration time-series signals of the device to be identified online. The system generates samples to be identified according to a preset sliding window. The length of the sliding window is consistent with that of the training phase, for example, 4 seconds. A 50% window overlap can be set between adjacent samples to be identified to improve the temporal continuity of abnormal responses.

[0113] Furthermore, steps S201 to S204 are executed sequentially on the current sample to be identified to obtain the impact energy map to be identified. Then, the impact energy map to be identified... Input the final anomaly classification model, output the device state probability vector. The current device state is determined based on the device state probability vector. If the maximum predicted probability corresponds to the normal category, the current sample is determined to be in a normal state; if the maximum predicted probability corresponds to an abnormal category, the abnormal category and its prediction confidence are output.

[0114] In one implementation, a continuous window confirmation mechanism can be set up to reduce false alarms caused by sporadic noise. When the same anomaly category is detected in consecutive... All samples to be identified were judged as abnormal, and their corresponding predicted probabilities were all higher than the alarm threshold. When this occurs, a device malfunction alarm is triggered. For example, We can choose 3. A value of 0.7 can be used. In addition, while outputting abnormal alarm results, the system can also retain the phase segment and order frequency band position with the strongest response in the corresponding impact energy map as auxiliary diagnostic information, which is used to indicate which rotating phase interval and which order range the fault impact is mainly concentrated in.

[0115] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for monitoring anomalies in key auxiliary equipment of thermal power units based on vibration data, characterized in that, Includes the following steps: S1. Collect the original vibration acceleration time-series signals of key auxiliary equipment of thermal power units, and label the equipment status to construct a training sample set; S2. Perform angular domain transformation and impact feature extraction on the original vibration acceleration time series signal without tachometer to obtain the angular domain impact feature signal, and convert the angular domain impact feature signal into an impact energy map with order-phase joint distribution; S3. Construct an order-phase dual-channel irregular convolutional network, input the impact energy map into the order-phase dual-channel irregular convolutional network for feature extraction and fusion, and output the device status anomaly classification result; The order-phase dual-channel irregular convolutional network includes: an order physical location encoding injection module, parallel phase convolution branches and order convolution branches, an order gated adaptive fusion module, an angle global cumulative pooling layer, and a classification layer; S4. Using the training sample set, the parameters of the order-phase dual-channel irregular convolutional network are iteratively optimized by minimizing the loss function to obtain the anomaly classification model; S5. Collect the real-time raw vibration acceleration time-series signal of the device to be identified, construct the impact energy map to be identified according to step S2, input the anomaly classification model, and output the device status anomaly identification result.

2. The method for abnormal monitoring of key auxiliary equipment in thermal power units based on vibration data according to claim 1, characterized in that, S2 specifically refers to: S21. Based on the instantaneous frequency shift estimation of the energy centroid of the multi-harmonic frequency band, a smooth instantaneous frequency shift estimation sequence is obtained; S22. Based on the smoothed instantaneous frequency conversion estimation sequence, the original vibration acceleration time series signal is resampled at equal angles to obtain the angular domain vibration signal; S23. Construct a low-quality factor atom library and a high-quality factor atom library, and use a matching pursuit strategy to sparsely reconstruct the angular domain vibration signal to obtain the angular domain impact characteristic signal. S24. Convert the angular domain impact characteristic signal into an impact energy map with a combined order-phase distribution.

3. The method for abnormal monitoring of key auxiliary equipment in thermal power units based on vibration data according to claim 2, characterized in that, S21 further includes: A short-time Fourier transform is performed on the original vibration acceleration time-series signal to obtain the time spectrum. Within a preset fundamental frequency search range, a coarse fundamental frequency estimate for each analysis frame is determined based on the degree of spectral energy concentration. Based on the coarse fundamental frequency estimate, harmonic search intervals are established for multiple low-order harmonics. Local band energy and local energy centroid frequency are calculated within each harmonic search interval. The local energy centroid frequencies of each harmonic are converted back to the fundamental frequency, and weighted fusion is performed using the corresponding local band energy as weight to obtain the instantaneous frequency shift estimate. Anomaly frame correction and smoothing processing are performed on the instantaneous frequency shift estimate sequence to obtain a smoothed instantaneous frequency shift estimate sequence.

4. The method for abnormal monitoring of key auxiliary equipment in thermal power units based on vibration data according to claim 3, characterized in that, S22 further includes: Based on the smooth instantaneous frequency shift estimation sequence, the frame-level instantaneous frequency shift is interpolated to the original sampling point position to obtain the sampling point-level instantaneous frequency shift sequence; Calculate the cumulative rotation angle sequence based on the instantaneous frequency conversion sequence at the sampling point level; Set the number of sampling points per rotation angle and establish an equal-angle sampling grid; use the cumulative rotation angle sequence as the independent variable and the original vibration signal as the dependent variable for interpolation and resampling to obtain the angular vibration signal.

5. The method for abnormal monitoring of key auxiliary equipment in thermal power units based on vibration data according to claim 4, characterized in that, S23 further includes: The angular domain vibration signal is divided into multiple overlapping angular domain blocks; real-valued decaying oscillation atoms are constructed, and the low-quality factor atom library and the high-quality factor atom library are constructed according to the range of decay coefficient values; for each angular domain block, the atom with the greatest correlation to the current residual signal is selected from the low-quality factor atom library and the high-quality factor atom library as the optimal atom, the projection coefficient is calculated and the residual signal is updated. Repeat the above atom selection and residual update process until the stopping condition is met; The selected atomic projection components within each corner domain block are accumulated and fused through overlapping regions to obtain the corner domain impact feature signal.

6. The method for abnormal monitoring of key auxiliary equipment in thermal power units based on vibration data according to claim 5, characterized in that, S24 further includes: Based on the number of sampling points per rotation angle, the angular domain impact characteristic signal is divided into multiple complete rotation cycles; the phase range of a single rotation cycle is equally divided into multiple phase segments; a local angular domain analysis window is extracted at the center of each phase segment in each rotation cycle, and a discrete Fourier transform is performed to obtain the local order spectrum; the preset order range is divided into multiple order frequency bands, and the average spectral energy in each phase segment and each order frequency band is statistically analyzed to construct the impact energy map with a size equal to the number of phase segments multiplied by the number of order frequency bands.

7. The method for abnormal monitoring of key auxiliary equipment in thermal power units based on vibration data according to claim 1, characterized in that, The order physical location encoding injection module in step S3 is further used for: The impact energy map is expanded into a multi-channel energy projection tensor using 1×1 convolution; a trainable order position encoding matrix is ​​defined, the size of which is the number of order frequency bands multiplied by the number of channels; the energy projection tensor and the order position encoding matrix are added element-wise along the order dimension to obtain the input tensor with injected order physical position encoding.

8. The method for abnormal monitoring of key auxiliary equipment in thermal power units based on vibration data according to claim 1, characterized in that, The parallel phase convolution branch and the order convolution branch in step S3 further include: the phase convolution branch uses a convolution kernel along the phase direction and adopts a cyclic filling method to output a phase feature map; the order convolution branch uses a convolution kernel along the order direction and adopts a zero-filling method to output an order feature map; the phase feature map and the order feature map have the same size.

9. The method for abnormal monitoring of key auxiliary equipment in thermal power units based on vibration data according to claim 1, characterized in that, The order-gated adaptive fusion module in step S3 is further used to: perform global average pooling on the order feature map along the phase dimension to obtain an order summary matrix; input the order summary matrix into a gating network to generate gating coefficients; broadcast the gating coefficients along the phase dimension, multiply them element-wise with the phase feature map, and then add them element-wise with the order feature map to obtain a fused feature map.

10. The method for abnormal monitoring of key auxiliary equipment in thermal power units based on vibration data according to claim 1, characterized in that, The angle global cumulative pooling layer in step S3 is further used to: simultaneously perform max pooling and average pooling on the fused feature map along the phase dimension; combine the max pooling result and the average pooling result according to a preset ratio to obtain an order-channel compressed feature matrix; The order-channel compressed feature matrix is ​​flattened and input into the classification layer, and the device state probability vector is output through the Softmax function.

Images