Method and system for detecting composite fault of rotating machine under time-varying working condition

Abstract

The invention relates to a rotary machine composite fault detection method and system under a time-varying working condition, and belongs to the technical field of rotary machine fault diagnosis. The method comprises the steps of collecting a non-stationary vibration signal under a variable speed working condition, fusing a third-order Renyi entropy and kinematics prior to construct an adaptive time-frequency enhancement model, and outputting weighted time-frequency distribution; global optimization is carried out through a priori constraint Viterbi algorithm, and an instantaneous frequency ridge line is extracted under the condition of no rotating speed sensor; performing angular domain resampling by taking the ridge line as a phase reference, and converting a non-stationary signal into an angular domain stationary signal; constructing an extended difference mode decomposition convex optimization model, and iteratively solving an optimal difference spectrum; and determining a frequency spectrum segmentation threshold value based on self-adaptive change point analysis, performing decoupling separation and time domain reconstruction on the optimal difference spectrum, and completing composite fault diagnosis through feature comparison. The method does not need an external tachometer, can adaptively extract and separate the composite fault features under the variable speed working condition, is high in diagnosis accuracy, and is good in engineering application value.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of rotating machinery fault diagnosis, in particular to a rotating machinery compound fault detection method and system under time-varying working conditions. BACKGROUND

[0002] As the core bearing component of the transmission system of large equipment such as wind power generation and heavy machinery, the running state of the planetary gearbox directly determines the running safety, reliability and service life of the whole machine. In the industrial production scene, the failure of the planetary gearbox will cause the equipment to stop and the production to be interrupted, and in severe cases, it will even cause a catastrophic safety accident, causing huge economic losses. Therefore, developing efficient and accurate fault diagnosis technology to realize the early identification and accurate positioning of planetary gearbox faults has important engineering application value for promoting the construction of equipment predictive maintenance system, reducing operation and maintenance cost, and avoiding safety risks.

[0003] However, in actual engineering applications, the planetary gearbox often works in complex working conditions with time-varying speed, and its fault diagnosis faces multiple technical challenges, and the existing related technologies have significant defects, as follows: First, the traditional spectrum analysis method is difficult to adapt to the characteristics of non-stationary signals: speed fluctuations will cause the interval of fault characteristic pulses in the planetary gearbox vibration signal to change dynamically with the speed, making the signal exhibit strong non-stationary characteristics. Traditional frequency spectrum analysis methods such as Fourier transform strictly rely on the stationarity assumption of the signal, and cannot effectively process time-varying frequency components, resulting in fault characteristics being masked by spectral blurring effects, making it difficult to achieve effective fault identification.

[0004] Second, the order tracking technology is limited by the speed measurement conditions: order tracking technology converts time-domain signals to angular-domain signals, which can eliminate the influence of speed fluctuations on fault characteristics, and is an ideal technical path to deal with variable speed conditions. However, the implementation of this technology requires accurate speed reference information, and the planetary gearbox often has compact structure, limited installation space, and strong vibration and electromagnetic interference in the working environment, making it difficult and costly to install a speed sensor, or even impossible to install one, which severely limits the application of order tracking technology in actual scenarios.

[0005] Third, traditional filtering and decomposition methods cannot accurately separate compound faults: under variable speed conditions, the fault signal of the planetary gearbox usually has the characteristics of weak amplitude, strong background noise and mutual modulation of multiple fault components. The traditional filtering method based on fixed frequency band division has insufficient resolution and cannot effectively separate overlapping fault characteristics; although the difference mode decomposition (DMD) and other decomposition techniques have certain advantages in the separation of compound faults of stationary signals, this technology has strict requirements for frequency stability and cannot cope with the frequency drift phenomenon under variable speed conditions, making it difficult to accurately decouple and separate compound fault components, and limiting the accuracy of fault diagnosis.

[0006] In summary, under the time-varying working condition without direct speed measurement, how to robustly extract fault features from strong noise and multi-component modulated non-stationary vibration signals and accurately separate complex faults has become a key technical bottleneck in the field of planetary gearbox fault diagnosis, and the existing technology cannot meet the needs of fault diagnosis accuracy and adaptability in engineering practice. SUMMARY

[0007] To this end, the present application aims to solve the technical problem of how to robustly extract clear fault features from strong noise and multi-component modulated non-stationary vibration signals and accurately separate weak fault enhancement and complex faults under the condition of variable speed operation of planetary gearboxes without direct speed measurement, thereby breaking through the bottleneck of failure of traditional methods under variable speed working conditions and achieving high-precision and high-reliability fault diagnosis.

[0008] To solve the above technical problems, the present application provides a rotating machinery complex fault detection method under time-varying working condition, which comprises the following steps: Step S1: Collecting the original vibration signal under variable speed working condition, performing short-time Fourier transform on it, and obtaining weighted time-frequency distribution based on Renyi entropy and prior meshing order; Step S2: Based on the weighted time-frequency distribution, a global optimization model is constructed by fusing energy statistical prior constraint and physical kinematics prior constraint to extract the instantaneous frequency ridge line which is energy concentrated, frequency continuous and consistent with mechanical motion law; Step S3: Based on the instantaneous frequency ridge line, the angle domain stationary signal is generated by performing equal angle resampling on the original vibration signal; the angle domain stationary signal is defined as a mixed signal; the angle domain stationary signal under the same type and working condition without fault operation of the rotating machinery is selected as a comparison signal, and the normalized Fourier amplitude spectrum of the mixed signal and the comparison signal is calculated by performing fast Fourier transform on them respectively; Step S4: Taking the normalized Fourier amplitude spectrum as input, an extended difference mode decomposition convex optimization model is constructed, and the weight vector is solved by maximum likelihood and sparse regularization iteration to obtain the optimal difference spectrum which removes common mode components and highlights fault features; Step S5: Based on the amplitude percentile sequence of the optimal difference spectrum, a frequency spectrum segmentation threshold for distinguishing fault principal components and noise components in the difference spectrum is calculated by using change point analysis; Step S6: Based on the frequency spectrum segmentation threshold, a frequency domain mask is constructed to extract the target fault frequency spectrum component in the optimal difference spectrum, and the inverse fast Fourier transform is performed on the extracted fault frequency spectrum component to reconstruct each independent time domain fault component, and the time domain fault component is compared with the theoretical fault feature order to output the complex fault diagnosis result.

[0009] In one embodiment of the present invention, step S1, the method for obtaining the weighted time-frequency distribution includes: The original vibration signal is subjected to a short-time Fourier transform to obtain its time-frequency energy distribution. Based on the time-frequency energy distribution, the normalized two-dimensional probability density function is calculated, expressed as: , in, The normalized two-dimensional probability density function represents the signal energy in the time domain. With frequency domain Normalized distribution on; The time-frequency energy distribution is expressed as follows: , The original vibration signal, For the Hanning window function, For integration time variable, As a time-centered variable, The imaginary unit; Introducing the theory of generalized information entropy, based on the aforementioned two-dimensional probability density function Calculate its The Rényi entropy of order Rényi is given by the formula: , in, Let be the order parameter, satisfying and ; Establish a kinematic harmonic frequency model for the planetary gearbox, and define the k-th type characteristic frequency component corresponding to the k-th type characteristic frequency component in the vibration signal of the planetary gearbox. The characteristic frequency of the subharmonic is Its expression is: , in, ω is the instantaneous angular velocity, representing the rotational speed of the gearbox; is the kinematic order relative to the axis; These are the harmonic order coefficients, representing the frequency components. The order of the i-th harmonic; Fusion of the Rényi entropy Calculate the adaptive time-frequency enhancement weight matrix using the kinematic harmonic frequency model. The formula is: , in, This is the noise time-frequency energy distribution threshold, used to distinguish between signal and noise components; Let be the conditional expectation operator, representing that when the condition is satisfied... Statistical average under the given conditions; For the planetary gearbox number Second harmonic, the first A normalized two-dimensional time-frequency probability density function subset of class feature frequency components; For the planetary gearbox number Class characteristic frequency components, the first subharmonics Rényi entropy; The adaptive time-frequency enhancement weight matrix With the time-frequency energy distribution Perform element-wise multiplication to generate a weighted time-frequency distribution. The calculation formula is: , in, This indicates element-wise multiplication.

[0010] In one embodiment of the present invention, step S2, the method for extracting an instantaneous frequency ridge line that is energy-concentrated, frequency-continuous, and conforms to the laws of mechanical motion, includes: The instantaneous frequency ridge extraction is transformed into a dynamic programming optimization problem. The optimal path is determined by minimizing the global cost function containing dual prior constraints, as shown in the following formula: , in, The optimal instantaneous frequency ridge estimation path for the final output; L represents the set of potential ridge paths; For time step; For the first Frequency index corresponding to each time step; This indicates that the weighted time-frequency distribution is at time 10:00. Frequency Index The amplitude at that point; This represents the amplitude penalty function corresponding to the amplitude dominance constraint. This represents a priori constraints on energy statistics; This represents the smoothing penalty function corresponding to the path smoothing constraint. Represents prior constraints in physical kinematics; The initial time step, The termination time step; Set the initial time step For all potential frequency indices Calculate the initial cost That is, it only includes the amplitude penalty value at the start time; For each subsequent time step Iterate through all potential frequency indices The cumulative cost is calculated using the following formula, and the optimal preceding path node is recorded: ; in, For the first Time step, frequency index The corresponding cumulative cost, This indicates selecting the optimal preceding node with the largest cumulative cost among all frequency indices at the previous time step; when At that time, find the terminal node with the highest cumulative cost. By backtracking from this terminal node to the optimal preceding node for all time steps, the complete optimal instantaneous frequency ridge can be obtained. .

[0011] In one embodiment of the present invention, the amplitude penalty function The expression is: , in, Indicates the time step At this point, after sorting the weighted time-frequency distribution amplitude sequence in descending order, the corresponding number is... The frequency range with the maximum amplitude; Sort by amplitude; The smoothing penalty function The expression is: , Where a and b are the adjacent time steps, respectively. and Instantaneous frequency index; The maximum allowable change threshold for instantaneous frequency ridge estimation depends on the physical inertial constraints of the system. The weighting coefficient is used to smooth out the penalty.

[0012] In one embodiment of the present invention, step S3, the method for obtaining the normalized Fourier amplitude spectrum includes: Based on the instantaneous frequency ridge Ridge frequency By combining the characteristic order relationship of the planetary gearbox, the instantaneous rotational speed of the reference shaft is obtained through inversion. The calculation formula is as follows: Where O is the characteristic order of the planetary gearbox; Based on the instantaneous rotational speed A quadratic polynomial fitting model for the rotation phase of the reference axis is constructed. This model is used to derive the non-uniform time-domain nodes for equal-angle resampling, ensuring that the signal is uniformly distributed in the angular domain after resampling. The recursive calculation formula for the k-th resampling time is: , in, For the first Each resampling time, The sampling interval is equal in angle, and k is the sequence number of the resampling point. , , These are the quadratic phase fitting coefficients; A series of non-uniformly distributed time-domain nodes are obtained through the above formula. The amplitudes of the original vibration signal at these nodes are extracted to obtain a stationary angular domain signal uniformly distributed over the angular domain. ; Define the angular domain stationary signal The signal is a mixed signal, comprising fault characteristic components of the rotating machinery and common-mode components of normal operation. A stationary angular domain signal obtained through the same resampling process under the same model and operating conditions as the rotating machinery, without faults, is selected as the comparison signal, denoted as [signature]. The comparison signal includes only the common-mode component of normal operation; For the mixed signal respectively Compared with the aforementioned signal Perform a Fast Fourier Transform to convert the time-domain signal to the frequency domain, and then normalize it using the L1 norm to eliminate the difference in energy magnitude between the two, thus obtaining the normalized Fourier amplitude spectra of the two.

[0013] In one embodiment of the present invention, step S4, the method for obtaining the optimal difference spectrum that removes common-mode components and highlights fault characteristics, includes: The set of spectral samples corresponding to the normalized Fourier amplitude spectrum of the mixed signal is defined as the mixed signal sample set {NFS}. M The set contains fault characteristic components and normal common-mode components, with a sample size of P; the set of spectral samples corresponding to the normalized Fourier amplitude spectrum of the comparison signal is defined as the comparison signal sample set {NFS}. R This set contains only normal common-mode components, and the number of samples is Q; The total spectrum sample set is obtained by fusing the two sample sets, and the total number of samples is And assign a logical label to each sample. When the i-th sample belongs to the set of contrast signals {NFS} R When}, mark When the i-th sample belongs to the mixed signal sample set {NFS} MWhen}, mark ; Define the feature vector of the i-th sample. This vector is composed of the amplitude sequence of the corresponding spectrum sample and is used to characterize the frequency domain features of the sample; A convex optimization model based on a binary logistic regression framework is constructed. The goal is to solve for the optimal parameter vector by combining maximum likelihood estimation with sparse regularization constraints. The formula is as follows: , in, The target loss function consists of a regularization term and a likelihood loss term, and is used to measure the deviation between the model's predicted values ​​and the true labels. Let be the parameter vector to be optimized. , For the weight vector, For bias; is the regularization coefficient, a non-negative constant used to control the sparsity of the weight vector; For the target loss function Perform iterative optimization until the difference in the loss function between two consecutive iterations is less than a preset threshold, at which point the optimal parameter vector is obtained. ,in This is the optimal weight vector; Based on the optimal weight vector The optimal difference spectrum is calculated and expressed as: , in, and These are the optimal weight vectors. The minimum and maximum values ​​in the array.

[0014] In one embodiment of the present invention, step S5, based on the amplitude percentile sequence of the optimal difference spectrum, uses change point analysis to calculate the spectral segmentation threshold for distinguishing the principal component of the fault and the noise component in the difference spectrum, including: Extract all amplitude data of the optimal difference spectrum SR to form an amplitude set. K represents the total number of spectral lines in the optimal difference spectrum. Let the amplitude of the k-th spectral line be denoted by ; sort the set of amplitudes in ascending order to obtain a sorted amplitude sequence; based on the sorted amplitude sequence, calculate its amplitude percentile sequence. , This represents the nth percentile amplitude of the sorted amplitude sequence. The length of the amplitude percentile sequence; Construct based on candidate split points The piecewise linear fitting model is used to fit the amplitude percentile sequence piecewise. The model expression is: , Where T is the piecewise linear fitting model; For sequence index variables; , and , These are the dividing points. The linear fitting coefficients of the two consecutive sequences; For each candidate split point The least squares method was used to solve for the linear fitting coefficients of the two segments separately. and And define a piecewise fitting error function to quantify the fitting accuracy: , , in, Candidate split points The fitting error of the corresponding first half of the sequence, Candidate split points The fitting error of the corresponding second half of the sequence; Define the total fitting error function Iterate through all candidate split points Calculate the total fitting error for each candidate point; find the candidate segmentation point that minimizes the total fitting error, which is the optimal segmentation point. The expression is: , Based on the optimal segmentation point Extract the corresponding index from the percentile sequence of magnitude. elements This element is used as the spectrum segmentation threshold.

[0015] In one embodiment of the present invention, step S6, the method for outputting the composite fault diagnosis result includes: Based on spectrum segmentation threshold A frequency domain separation mask is constructed to filter the spectral components belonging to the target fault in the optimal difference spectrum, retaining only the spectral amplitude of the spectral line index within the fault principal component interval, and suppressing the spectral amplitude of noise and minor components. The frequency domain separation mask is multiplied element-wise with the optimal difference spectrum to extract the spectral components of the target fault. , is represented as: , in, For spectral line indexing; For the spectral components Performing an inverse fast Fourier transform converts the frequency domain signal back to the time domain, reconstructing the independent fault component signal. , is represented as: , in, This represents the operator for taking the real part; This represents the inverse fast Fourier transform operator; For the fault component signal Perform order spectrum analysis to extract the order of its fault characteristics. To obtain the theoretical fault characteristic order of rotating machinery ; The order of extracted features is determined by the relative error criterion. Compared with the theoretical fault characteristic order A comparison is performed to determine if a corresponding fault exists. The criterion formula is as follows: , in, For harmonic order, The allowable relative error threshold; Repeat the above steps to extract, reconstruct, and compare the order of all target fault spectral components in the optimal difference spectrum, and finally output a composite fault diagnosis result containing all fault types.

[0016] Based on the same inventive concept, this invention provides a system for detecting complex faults in rotating machinery under time-varying operating conditions, comprising the following modules: The data acquisition and time-frequency enhancement module is used to acquire the raw vibration signal under variable speed conditions, perform a short-time Fourier transform on it, and based on... The weighted time-frequency distribution is obtained by combining the Rényi entropy of order of meshing with the prior of meshing order. The rotational speed sensorless ridge extraction module is used to extract instantaneous frequency ridges that are energy-concentrated, frequency-continuous, and conform to the laws of mechanical motion by constructing a global optimization model based on the weighted time-frequency distribution and by fusing energy statistical prior constraints and physical kinematic prior constraints. The angular domain resampling and spectrum calculation module is used to perform equal-angle resampling on the original vibration signal based on the instantaneous frequency ridge to generate an angular domain stationary signal; the angular domain stationary signal is defined as a hybrid signal; an angular domain stationary signal of the same model and operating condition as the rotating machinery under fault-free operation is selected as a comparison signal, and fast Fourier transform is performed on the hybrid signal and the comparison signal respectively to calculate their normalized Fourier amplitude spectra; The convex optimization difference spectrum generation module is used to take the normalized Fourier amplitude spectrum as input, construct an extended difference mode decomposition convex optimization model, and obtain the optimal difference spectrum that removes common mode components and highlights fault characteristics by solving the weight vector through maximum likelihood and sparse regularization iteration. An adaptive threshold decision module is used to calculate a spectrum segmentation threshold for distinguishing fault principal components and noise components in the difference spectrum based on the amplitude percentile sequence of the optimal difference spectrum using change point analysis. The signal decoupling and diagnostic output module is used to construct a frequency domain mask based on the spectrum segmentation threshold, extract the target fault spectrum component in the optimal difference spectrum, perform inverse fast Fourier transform on the extracted fault spectrum component to reconstruct each independent time-domain fault component, compare the time-domain fault component with the theoretical fault feature order, and output the composite fault diagnosis result.

[0017] Furthermore, the present invention provides an electronic device including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the steps of the rolling bearing fault diagnosis method described above.

[0018] The technical solution of the present invention has the following advantages compared with the prior art: First, the robustness of time-frequency enhancement and ridge extraction is improved: This invention innovatively constructs an adaptive time-frequency enhancement model that integrates kinematic harmonic priors and α-order Rényi entropy. Addressing the technical shortcomings of traditional algorithms in strong noise environments, such as ridge line breakage and susceptibility to local extrema, a Viterbi global path search mechanism based on dual prior constraints is proposed. Unlike the traditional Viterbi algorithm, which relies solely on pure mathematical deduction of state transition probabilities, this mechanism deeply integrates the physical inertia priors of rotating machinery and the statistical priors of energy distribution, reconstructing the path search cost function. It strengthens the priority of energy-dominant node traversal through an amplitude ranking mechanism and suppresses non-physical frequency jumps by combining physical constraints of frequency continuity. Finally, through dynamic programming search, it obtains the instantaneous frequency ridge with the minimum cumulative cost and optimal frequency change continuity, significantly improving the anti-interference capability and accuracy of ridge extraction.

[0019] Secondly, the sensorless closed-loop architecture overcomes the limitations of variable speed diagnostics: It constructs a sensorless soft-measurement closed-loop system of "ridge extraction—angular domain resampling—differential mode decoupling," using the extracted instantaneous frequency ridge as a virtual tachometer to drive the calculation of order tracking, converting time-varying non-stationary vibration signals into angular domain stationary signals, thus eliminating spectral ambiguity and modulation effects caused by speed fluctuations at their source. This architecture breaks through the strict dependence of traditional differential mode decomposition (DMD) algorithms on spectral stability, solving its technical bottleneck of being difficult to adapt to variable speed signals; at the same time, only a single-channel accelerometer is needed to complete the entire diagnostic process, eliminating the need for additional hardware encoders, effectively reducing equipment deployment costs and space constraints in industrial settings, and improving the engineering applicability of the technology.

[0020] Third, the accuracy optimization of adaptive decoupling for complex faults: Addressing the challenge of feature separation in complex faults, an adaptive threshold decision module based on convex optimization learning and change point detection is constructed, transforming the complex fault decoupling problem into a statistical mutation detection process. By iteratively solving for the optimal weight vector using maximum likelihood estimation combined with sparse regularization, an optimal difference spectrum that highlights fault features and suppresses common-mode components is generated. Then, change point analysis is used to segment and fit the percentile sequence of the difference spectrum amplitude, adaptively determining the spectrum segmentation threshold, completely eliminating the reliance on manual empirical parameters in traditional methods. This design ensures that the diagnostic model can accurately isolate independent fault modes even under non-ideal conditions where fault features are weak and background noise and signal energy are deeply coupled, guaranteeing high accuracy and reliability in complex fault diagnosis. Attached Figure Description

[0021] To make the content of this invention easier to understand, the invention will be further described in detail below with reference to specific embodiments and accompanying drawings.

[0022] Figure 1 This is a flowchart of a method for detecting complex faults in rotating machinery under time-varying operating conditions provided by the present invention; Figure 2 This is the time-domain waveform diagram of the original composite fault signal of the planetary gearbox under variable speed conditions in this invention, corresponding to a mixed fault including cracks in the root of the planetary gear teeth and cracks in the root of the sun gear teeth; Figure 3 This is a waveform of a stationary signal in the angular domain obtained by angular domain resampling using the extracted instantaneous frequency ridge in this invention. Figure 4 The time-domain waveforms of the independent fault components and their corresponding order spectra obtained by differential mode decomposition (DMD) and adaptive threshold decoupling in this invention are shown in (a), where (b) represents the extracted sun gear fault component, (c) represents the extracted planet gear component, and (d) represents the order spectrum of the planet gear fault component. Figure 5This is a schematic diagram of the structure of a composite fault detection system for rotating machinery under time-varying operating conditions provided by the present invention; Explanation of reference numerals on the accompanying drawings: 100. Data Acquisition and Time-Frequency Enhancement Module; 200. Ridge Extraction Module without Rotational Speed ​​Sensor; 300. Angular Domain Resampling and Spectrum Calculation Module; 400. Convex Optimization Difference Spectrum Generation Module; 500. Adaptive Threshold Decision Module; 600. Signal Decoupling and Diagnostic Output Module. Detailed Implementation

[0023] The present invention will be further described below with reference to the accompanying drawings and specific embodiments, so that those skilled in the art can better understand and implement the present invention. However, the embodiments described are not intended to limit the present invention.

[0024] Example 1: Reference Figure 1 As shown, the present invention provides a method for detecting complex faults in rotating machinery under time-varying operating conditions, comprising: Step S1: Acquire the original vibration signal under variable speed conditions, perform a short-time Fourier transform on it, and based on... The weighted time-frequency distribution is obtained by combining the Rényi entropy of order of meshing with the prior of meshing order. Step S2: Based on the weighted time-frequency distribution, a global optimization model is constructed by fusing energy statistical prior constraints and physical kinematic prior constraints to extract instantaneous frequency ridges that are energy-concentrated, frequency-continuous, and conform to the laws of mechanical motion. Step S3: Based on the instantaneous frequency ridge, perform equal-angle resampling on the original vibration signal to generate an angular domain stationary signal; define the angular domain stationary signal as a mixed signal; select the angular domain stationary signal of the same model and operating conditions as the rotating machinery operating without faults as a comparison signal, and perform fast Fourier transform on the mixed signal and the comparison signal respectively to calculate their normalized Fourier amplitude spectra; Step S4: Using the normalized Fourier amplitude spectrum as input, construct an extended differential mode decomposition convex optimization model, and solve the weight vector through maximum likelihood and sparse regularization iterations to obtain the optimal differential spectrum that removes common mode components and highlights fault characteristics. Step S5: Based on the amplitude percentile sequence of the optimal difference spectrum, use change point analysis to calculate the spectrum segmentation threshold used to distinguish the fault principal component and the noise component in the difference spectrum. Step S6: Construct a frequency domain mask based on the spectrum segmentation threshold, extract the target fault spectrum component in the optimal difference spectrum, perform inverse fast Fourier transform on the extracted fault spectrum component to reconstruct each independent time-domain fault component, compare the time-domain fault component with the theoretical fault feature order, and output the composite fault diagnosis result.

[0025] Further, in step S1, this embodiment collects the original composite fault signal of the planetary gearbox under variable speed conditions. The composite fault is a crack in the root of the planetary gear teeth and a crack in the root of the sun gear teeth, such as... Figure 2 As shown. The signal sampling frequency is 512Hz, and the number of teeth on the sun gear is... Number of planetary gear teeth There are three planetary gears in total, and the number of teeth on the gear ring is... The input shaft frequency increased from 0Hz to 10Hz within 6 seconds. This is in response to the aforementioned non-stationary vibration signal. Perform the following specific processing steps: First of all, Figure 2 The original vibration signal shown Perform a Short-Time Fourier Transform (STFT) using the Hanning window function. As a sliding window function, its time-frequency energy distribution is obtained. ,as follows: ; Based on the time-frequency energy distribution Calculate the normalized two-dimensional probability density function. Characterizing signal energy in the time domain With frequency domain The normalized distribution on; its expression is: , in, For integration time variable, As a time-centered variable, It is the imaginary unit.

[0026] To address the issue of signal energy dispersion under variable speed operating conditions, the theory of generalized information entropy is introduced. In this embodiment, an order α=3 is selected to calculate the third-order Rényi entropy of the time-frequency energy distribution, quantifying the degree of signal energy concentration. The calculation formula is as follows: .

[0027] Simultaneously, based on the aforementioned tooth number parameters of the planetary gearbox, and utilizing the kinematic equations of planetary gear transmission, a kinematic harmonic frequency model is constructed, incorporating the doubling relationship between the sun gear rotation frequency and the meshing frequency. The k-th type characteristic frequency component corresponding to the k-th type of characteristic frequency component in the planetary gearbox vibration signal is defined. The characteristic frequency of the subharmonic is Its expression is: , in, ω is the instantaneous angular velocity, representing the rotational speed of the gearbox; The kinematic order relative to the shaft, such as gear ratio or meshing order; These are the harmonic order coefficients, representing the frequency components. The order of the i-th harmonic.

[0028] Finally, the calculated third-order Rényi entropy is fused. An adaptive time-frequency enhancement weight matrix is ​​generated using the kinematic harmonic frequency model. The calculation formula is as follows: , in, This is the noise time-frequency energy distribution threshold, used to distinguish between signal and noise components; Let be the conditional expectation operator, representing that when the condition is satisfied... Statistical average under the given conditions; For the planetary gearbox number Second harmonic, the first A normalized two-dimensional time-frequency probability density function subset of class feature frequency components; For the planetary gearbox number Class characteristic frequency components, the first subharmonics Rényi entropy.

[0029] The adaptive time-frequency enhancement weight matrix With the time-frequency energy distribution Modulation is performed, and a weighted time-frequency distribution is output. The calculation formula is: , in, This indicates element-wise multiplication. This distribution effectively suppresses background noise and highlights the true frequency components related to gear meshing.

[0030] In step S2, the weighted time-frequency distribution output in step S1 is... A global path search model based on the Viterbi algorithm with prior constraints is constructed. This step aims to address the ridge line breakage problem caused by signal non-stationarity under variable speed conditions (acceleration process from 0-10Hz in this embodiment).

[0031] Specifically, this invention models the extraction process of instantaneous frequency ridges as a dynamic programming optimization problem, defining an objective function for global path search. This function requires the final extracted frequency path... To minimize the cumulative cost over the entire time period, a Viterbi path search algorithm based on prior constraints is defined, which finds the optimal path by minimizing the global cost function that includes the prior constraint term: , in, The optimal instantaneous frequency ridge estimation path for the final output; L represents the set of potential ridge paths; For time step; For the first Frequency index corresponding to each time step; This indicates that the weighted time-frequency distribution is at time 10:00. Frequency Index The amplitude at that point; This represents the amplitude penalty function corresponding to the amplitude dominance constraint. This represents a priori constraints on energy statistics; This represents the smoothing penalty function corresponding to the path smoothing constraint. Represents prior constraints in physical kinematics; The initial time step, The termination time step; To achieve the above objective function, we first construct an amplitude penalty function based on sorting. Since the third-order Rényi entropy has already been passed in step S1 Kinematic priors on time-frequency energy distribution Weighted enhancement has been applied, and this function is used here to further drive the search algorithm to prioritize traversing energy-dominant nodes (i.e., the actual fault or meshing frequency components). The calculation formula is as follows: , in, Indicates the time step At this point, after sorting the weighted time-frequency distribution amplitude sequence in descending order, the corresponding number is... The frequency range with the maximum amplitude; The magnitude is sorted by position.

[0032] Meanwhile, considering that the rotational speed of the planetary gearbox in this embodiment changes continuously, there is no physical instantaneous drastic change in frequency. Therefore, a path smoothing penalty function based on Euclidean distance is defined. This is used to suppress non-physical frequency jumps caused by strong background noise. In this embodiment, a reasonable frequency change threshold is set. A penalty is applied when the frequency change between adjacent time points exceeds this threshold. The calculation formula is as follows: , Where a and b are the adjacent time steps, respectively. and Instantaneous frequency index; The maximum allowable change threshold for instantaneous frequency ridge estimation depends on the physical inertial constraints of the system. The weighting coefficient is used to smooth out the penalty.

[0033] Set the initial time step For all potential frequency indices Calculate the initial cost That is, it only includes the amplitude penalty value at the start time; For each subsequent time step Iterate through all potential frequency indices The cumulative cost is calculated using the following formula, and the optimal preceding path node is recorded: ; in, For the first Time step, frequency index The corresponding cumulative cost, This indicates selecting the optimal preceding node with the largest cumulative cost among all frequency indices at the previous time step; when At that time, find the terminal node with the highest cumulative cost. By tracing back from this terminal node to the optimal preceding node for all time steps, the optimal instantaneous frequency ridge that accurately reflects the variation law of the input shaft speed is obtained. This provides an accurate phase reference for subsequent angular domain resampling.

[0034] In step S3, this embodiment utilizes the optimal instantaneous frequency ridge. As a phase reference, angular domain resampling is performed on the original vibration signal. This step aims to eliminate the severe frequency ambiguity and modulation effect caused by the rotational speed increasing from 0Hz to 10Hz in this embodiment.

[0035] First, based on the extracted instantaneous frequency ridge... Ridge frequency By combining the characteristic order relationship of the planetary gearbox, the instantaneous rotational speed of the reference shaft is obtained through inversion. The calculation formula is as follows: , Where O represents the characteristic order of the planetary gearbox. Next, in order to convert the non-stationary time-domain signal into a stationary angular-domain signal, based on the instantaneous rotational speed... A quadratic polynomial fitting model of the rotation phase of the reference axis is constructed to solve for the equal-angle resampling time. Assume the rotation phase of the reference axis... It follows a quadratic polynomial distribution, as shown in the following equation: , in, , , These are the quadratic phase fitting coefficients. Based on equal-angle sampling intervals. Given the constraints, the k-th equal-angle resampling time t is derived. k The recursive calculation formula is as follows: ; A series of non-uniformly distributed time-domain nodes are obtained through the above formula. The amplitudes of the original vibration signal at these nodes are extracted to obtain a stationary angular domain signal uniformly distributed over the angular domain. The reconstructed signal after this step is as follows: Figure 3 As shown, it can be seen that the frequency components that originally varied with time in the time domain were successfully transformed into a stable waveform with constant frequency in the angular domain, laying the foundation for subsequent fault feature separation.

[0036] Define the angular domain stationary signal The signal is a mixed signal, comprising fault characteristic components of the rotating machinery and common-mode components of normal operation. A stationary angular domain signal obtained through the same resampling process under the same model and operating conditions as the rotating machinery, without faults, is selected as the comparison signal, denoted as [signature]. The comparison signal includes only the common-mode component of normal operation; For the mixed signal respectively Compared with the aforementioned signal Perform a Fast Fourier Transform (FFT) to convert the time-domain signal to the frequency domain, and then normalize using the L1 norm to eliminate the difference in energy magnitude between the two, obtaining the normalized Fourier amplitude spectra of both: , , in, Describing the L1 norm, The normalized Fourier amplitude spectrum of the mixed signal. To compare the normalized Fourier amplitude spectrum of the signal.

[0037] In step S4, the set of spectral samples corresponding to the normalized Fourier amplitude spectrum of the mixed signal is defined as the mixed signal sample set {NFS}. M The set contains fault characteristic components and normal common-mode components, with a sample size of P; the set of spectral samples corresponding to the normalized Fourier amplitude spectrum of the comparison signal is defined as the comparison signal sample set {NFS}. R This set contains only normal common-mode components, and the number of samples is Q; The total spectrum sample set is obtained by fusing the two sample sets, and the total number of samples is And assign a logical label to each sample. When the i-th sample belongs to the set of contrast signals {NFS} R When}, mark When the i-th sample belongs to the mixed signal sample set {NFS} M When}, mark Define the feature vector of the i-th sample. This vector is composed of the amplitude sequence of the corresponding spectral sample and is used to characterize the frequency domain features of the sample.

[0038] Next, a convex optimization model based on a binary logistic regression framework is constructed. The goal is to solve for the optimal parameter vector by combining maximum likelihood estimation with sparse regularization constraints. The formula is as follows: , in, The target loss function consists of a regularization term and a likelihood loss term, and is used to measure the deviation between the model's predicted values ​​and the true labels. Let be the parameter vector to be optimized. , For the weight vector, For bias; is the regularization coefficient, a non-negative constant used to control the sparsity of the weight vector.

[0039] For the target loss function Perform iterative optimization until the difference in the loss function between two consecutive iterations is less than a preset threshold, at which point the optimal parameter vector is obtained. ,in This is the optimal weight vector; Based on the optimal weight vector The optimal difference spectrum is calculated and expressed as: ,in, and These are the optimal weight vectors. The minimum and maximum elements in the spectrum. In this embodiment, the generated optimal difference spectrum... The inherent meshing frequency interference in the planetary gearbox was successfully suppressed, highlighting the potential fault impact frequency band.

[0040] In step S5, the aim is to determine how to obtain the optimal difference spectrum from step S4. In the process, a boundary is automatically determined to separate the effective components representing fault characteristics from the background components representing noise.

[0041] Extract all amplitude data of the optimal difference spectrum SR to form an amplitude set. K represents the total number of spectral lines in the optimal difference spectrum. Let the amplitude of the k-th spectral line be denoted by ; sort the set of amplitudes in ascending order to obtain a sorted amplitude sequence; based on the sorted amplitude sequence, calculate its amplitude percentile sequence. , This represents the nth percentile amplitude of the sorted amplitude sequence. The length of the amplitude percentile sequence; Construct based on candidate split points The piecewise linear fitting model (i.e., the change point) is used to fit the amplitude percentile sequence piecewise. The model expression is: , Where T is the piecewise linear fitting model; For sequence index variables; , and , These are the dividing points. The linear fitting coefficients of the two sequences.

[0042] For each candidate split point The least squares method was used to solve for the linear fitting coefficients of the two segments separately. and And define a piecewise fitting error function to quantify the fitting accuracy: , , in, Candidate split points The fitting error of the corresponding first half of the sequence (high amplitude range, usually corresponding to the fault principal component), Candidate split points The fitting error of the corresponding second half of the sequence (low amplitude range, usually corresponding to noise or minor components).

[0043] Define the total fitting error function Iterate through all candidate split points Calculate the total fitting error for each candidate point; find the candidate segmentation point that minimizes the total fitting error, which is the optimal segmentation point. The expression is: , Based on the optimal segmentation point Extract the corresponding index from the percentile sequence of magnitude. elements This element is used as the spectrum segmentation threshold.

[0044] In step S6, the optimal spectral threshold calculated in step S5 is... Frequency domain decoupling is performed on the optimal difference spectrum SR to extract the dominant frequency band of a specific fault, and the independent time-domain fault components are reconstructed through inverse Fourier transform.

[0045] Based on spectrum segmentation threshold A frequency domain separation mask is constructed to filter the spectral components belonging to the target fault in the optimal difference spectrum, retaining only the spectral amplitude of the spectral line index within the fault principal component interval, and suppressing the spectral amplitude of noise and minor components. The frequency domain separation mask is multiplied element-wise with the optimal difference spectrum to extract the spectral components of the target fault. , is represented as: , Where k is the spectral line index. For the spectral components... Performing an inverse fast Fourier transform converts the frequency domain signal back to the time domain, reconstructing the independent fault component signal. , is represented as: , in, This represents the operator for taking the real part; This represents the inverse fast Fourier transform operator. For the fault component signal... Perform order spectrum analysis to extract the order of its fault characteristics. To obtain the theoretical fault characteristic order of rotating machinery ; The order of extracted features is determined by the relative error criterion. Compared with the theoretical fault characteristic order A comparison is performed to determine if a corresponding fault exists. The criterion formula is as follows: , in, For harmonic order, This is the allowable relative error threshold.

[0046] Repeat the above steps to extract, reconstruct, and compare the order of all target fault spectral components in the optimal difference spectrum SR, and finally output a composite fault diagnosis result containing all fault types.

[0047] Figure 4 The results of fault feature decoupling and extraction obtained after the original planetary gearbox variable speed composite fault signal is processed through steps S1-S6 are shown. The original aliased composite fault signal is successfully separated into two independent fault components.

[0048] like Figure 4 (a) and Figure 4 As shown in (b), the first component obtained after separation is the sun gear fault component. Figure 4In the order spectrum of (b), a prominent characteristic order with a value of 2.53 can be clearly observed. This is very close to the theoretical characteristic order of 2.5 for a planetary gearbox sun gear failure, and the relative error is much smaller than the preset threshold. In addition, obvious harmonic components were observed at its octaves, which is consistent with the failure mechanism of local tooth breakage in the gear.

[0049] like Figure 4 (c) and Figure 4 As shown in (d), the separated second component is the planetary gear fault component. Figure 4 In the order spectrum of (d), the extracted feature order is 0.4166, which is in excellent agreement with the theoretical feature order of 0.417 for planetary gear failure. Furthermore, the time-domain waveform of this component clearly demonstrates the unique periodic impact characteristics of planetary gear failure, without being affected by the sun gear failure component.

[0050] In conclusion, through comparison Figure 4 The extraction results and theoretical fault model clearly indicate that the planetary gearbox simultaneously exhibits both sun gear tooth root crack faults and planet gear tooth root crack faults. This result confirms that the present invention can achieve accurate decoupling and identification of compound faults under complex operating conditions with variable speed and without tachometer assistance.

[0051] Example 2: like Figure 5 As shown, this embodiment provides a composite fault detection system for rotating machinery under time-varying operating conditions, including: a data acquisition and time-frequency enhancement module 100, a ridge line extraction module without rotational speed sensing 200, an angular domain resampling and spectrum calculation module 300, a convex optimization difference spectrum generation module 400, an adaptive threshold decision module 500, and a signal decoupling and diagnostic output module 600.

[0052] The data acquisition and time-frequency enhancement module 100 is used to acquire the original vibration signal under variable speed conditions, perform a short-time Fourier transform on it, and based on... The weighted time-frequency distribution is obtained by combining the Rényi entropy of order of meshing with the prior of meshing order. The non-rotational speed sensing ridge extraction module 200 is used to extract instantaneous frequency ridges that are energy-concentrated, frequency-continuous, and conform to the laws of mechanical motion by constructing a global optimization model based on the weighted time-frequency distribution and by fusing energy statistical prior constraints and physical kinematic prior constraints. The angular domain resampling and spectrum calculation module 300 is used to perform equal-angle resampling on the original vibration signal based on the instantaneous frequency ridge to generate an angular domain stationary signal; the angular domain stationary signal is defined as a mixed signal; an angular domain stationary signal that is running without faults under the same model and working conditions as the rotating machinery is selected as a comparison signal, and fast Fourier transform is performed on the mixed signal and the comparison signal respectively to calculate their normalized Fourier amplitude spectra; The convex optimization difference spectrum generation module 400 is used to take the normalized Fourier amplitude spectrum as input, construct an extended difference mode decomposition convex optimization model, and obtain the optimal difference spectrum that removes common mode components and highlights fault characteristics by solving the weight vector through maximum likelihood and sparse regularization iteration. The adaptive threshold decision module 500 is used to calculate the spectrum segmentation threshold for distinguishing fault principal components and noise components in the difference spectrum based on the amplitude percentile sequence of the optimal difference spectrum using change point analysis. The signal decoupling and diagnostic output module 600 is used to construct a frequency domain mask based on the spectrum segmentation threshold, extract the target fault spectrum component in the optimal difference spectrum, perform inverse fast Fourier transform on the extracted fault spectrum component to reconstruct each independent time-domain fault component, compare the time-domain fault component with the theoretical fault feature order, and output the composite fault diagnosis result.

[0053] Example 3: This embodiment provides an electronic device, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, it implements the steps of the rolling bearing fault diagnosis method described in Embodiment 1.

[0054] Those skilled in the art will understand that embodiments of this application can be provided as methods, systems, or computer program products. Therefore, this application can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, this application can take the form of a computer program product embodied on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0055] This application is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of this application. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart... Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0056] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The functions specified in one or more boxes. These computer program instructions may also be loaded onto a computer or other programmable data processing apparatus to cause a series of operational steps to be performed on the computer or other programmable apparatus to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable apparatus for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0057] Obviously, the above embodiments are merely illustrative examples for clear explanation and are not intended to limit the implementation. Those skilled in the art will recognize that other variations or modifications can be made based on the above description. It is neither necessary nor possible to exhaustively list all possible implementations here. However, obvious variations or modifications derived therefrom are still within the scope of protection of this invention.

Claims

1. A method for detecting complex faults in rotating machinery under time-varying operating conditions, characterized in that, Includes the following steps: Step S1: Acquire the original vibration signal under variable speed conditions, perform a short-time Fourier transform on it, and based on... The weighted time-frequency distribution is obtained by combining the Rényi entropy of order of meshing with the prior of meshing order. Step S2: Based on the weighted time-frequency distribution, a global optimization model is constructed by fusing energy statistical prior constraints and physical kinematic prior constraints to extract instantaneous frequency ridges that are energy-concentrated, frequency-continuous, and conform to the laws of mechanical motion. Step S3: Based on the instantaneous frequency ridge, perform equal-angle resampling on the original vibration signal to generate an angular domain stationary signal; define the angular domain stationary signal as a mixed signal; select the angular domain stationary signal of the same model and operating conditions as the rotating machinery operating without faults as a comparison signal, and perform fast Fourier transform on the mixed signal and the comparison signal respectively to calculate their normalized Fourier amplitude spectra; Step S4: Using the normalized Fourier amplitude spectrum as input, construct an extended differential mode decomposition convex optimization model, and solve the weight vector through maximum likelihood and sparse regularization iterations to obtain the optimal differential spectrum that removes common mode components and highlights fault characteristics. Step S5: Based on the amplitude percentile sequence of the optimal difference spectrum, use change point analysis to calculate the spectrum segmentation threshold used to distinguish the fault principal component and the noise component in the difference spectrum. Step S6: Construct a frequency domain mask based on the spectrum segmentation threshold, extract the target fault spectrum component in the optimal difference spectrum, perform inverse fast Fourier transform on the extracted fault spectrum component to reconstruct each independent time-domain fault component, compare the time-domain fault component with the theoretical fault feature order, and output the composite fault diagnosis result.

2. The method for detecting combined faults in rotating machinery under time-varying operating conditions according to claim 1, characterized in that: In step S1, the methods for obtaining the weighted time-frequency distribution include: The original vibration signal is subjected to a short-time Fourier transform to obtain its time-frequency energy distribution. Based on the time-frequency energy distribution, the normalized two-dimensional probability density function is calculated, expressed as: ， in, The normalized two-dimensional probability density function represents the signal energy in the time domain. With frequency domain Normalized distribution on; The time-frequency energy distribution is expressed as follows: , The original vibration signal, For the Hanning window function, For integration time variable, As a time-centered variable, The imaginary unit; Introducing the theory of generalized information entropy, based on the aforementioned two-dimensional probability density function Calculate its The Rényi entropy of order Rényi is given by the formula: ， in, Let be the order parameter, satisfying and ; Establish a kinematic harmonic frequency model for the planetary gearbox, and define the k-th type characteristic frequency component corresponding to the k-th type characteristic frequency component in the vibration signal of the planetary gearbox. The characteristic frequency of the subharmonic is Its expression is: ， Where Ω is the instantaneous angular velocity, representing the rotational speed of the gearbox; is the kinematic order relative to the axis; These are the harmonic order coefficients, representing the frequency components. The order of the i-th harmonic; Fusion of the Rényi entropy Calculate the adaptive time-frequency enhancement weight matrix using the kinematic harmonic frequency model. The formula is: ， in, This is the noise time-frequency energy distribution threshold, used to distinguish between signal and noise components; Let be the conditional expectation operator, representing that when the condition is satisfied... Statistical average under the given conditions; For the planetary gearbox number Second harmonic, the first A normalized two-dimensional time-frequency probability density function subset of class feature frequency components; For the planetary gearbox number Class characteristic frequency components, the first subharmonics Rényi entropy; The adaptive time-frequency enhancement weight matrix With the time-frequency energy distribution Perform element-wise multiplication to generate a weighted time-frequency distribution. The calculation formula is: ， in, This indicates element-wise multiplication.

3. The method for detecting combined faults in rotating machinery under time-varying operating conditions according to claim 1, characterized in that: In step S2, the method for extracting the instantaneous frequency ridge line that is energy-concentrated, frequency-continuous, and conforms to the laws of mechanical motion includes: The instantaneous frequency ridge extraction is transformed into a dynamic programming optimization problem. The optimal path is determined by minimizing the global cost function containing dual prior constraints, as shown in the following formula: ， in, The optimal instantaneous frequency ridge estimation path for the final output; L represents the set of potential ridge paths; For time step; For the first Frequency index corresponding to each time step; This indicates that the weighted time-frequency distribution is at time [time]. Frequency Index The amplitude at that point; This represents the amplitude penalty function corresponding to the amplitude dominance constraint. This represents a priori constraints on energy statistics; This represents the smoothing penalty function corresponding to the path smoothing constraint. Represents prior constraints in physical kinematics; The initial time step, The termination time step; Set the initial time step For all potential frequency indices Calculate the initial cost That is, it only includes the amplitude penalty value at the start time; For each subsequent time step Iterate through all potential frequency indices The cumulative cost is calculated using the following formula, and the optimal preceding path node is recorded: ； in, For the first Time step, frequency index The corresponding cumulative cost, This indicates selecting the optimal preceding node with the largest cumulative cost among all frequency indices at the previous time step; when At that time, find the terminal node with the highest cumulative cost. By backtracking from this terminal node to the optimal preceding node for all time steps, the complete optimal instantaneous frequency ridge can be obtained. .

4. The method for detecting combined faults in rotating machinery under time-varying operating conditions according to claim 3, characterized in that: The amplitude penalty function The expression is: ， in, Indicates the time step At this point, after sorting the weighted time-frequency distribution amplitude sequence in descending order, the corresponding number is... The frequency range with the maximum amplitude; Sort by amplitude; The smoothing penalty function The expression is: ， Where a and b are the adjacent time steps, respectively. and Instantaneous frequency index; The maximum allowable change threshold for instantaneous frequency ridge estimation depends on the physical inertial constraints of the system. The weighting coefficient is used to smooth out the penalty.

5. The method for detecting combined faults in rotating machinery under time-varying operating conditions according to claim 1, characterized in that: In step S3, the method for obtaining the normalized Fourier amplitude spectrum includes: Based on the instantaneous frequency ridge Ridge frequency By combining the characteristic order relationship of the planetary gearbox, the instantaneous rotational speed of the reference shaft is obtained through inversion. The calculation formula is as follows: ， Where O represents the characteristic order of the planetary gearbox; Based on the instantaneous rotational speed A quadratic polynomial fitting model for the rotation phase of the reference axis is constructed. This model is used to derive the non-uniform time-domain nodes for equal-angle resampling, ensuring that the signal is uniformly distributed in the angular domain after resampling. The recursive calculation formula for the k-th resampling time is: ， in, For the first Each resampling time, The sampling interval is equal in angle, and k is the sequence number of the resampling point. , , These are the quadratic phase fitting coefficients; A series of non-uniformly distributed time-domain nodes are obtained through the above formula. The amplitudes of the original vibration signal at these nodes are extracted to obtain a stationary angular domain signal uniformly distributed over the angular domain. ; Define the angular domain stationary signal The signal is a mixed signal, comprising fault characteristic components of the rotating machinery and common-mode components of normal operation. A stationary angular domain signal obtained through the same resampling process under the same fault-free operation conditions as the rotating machinery of the same model is selected as the comparison signal, denoted as... The comparison signal includes only the common-mode component of normal operation; For the mixed signal respectively Compared with the aforementioned signal Perform a Fast Fourier Transform to convert the time-domain signal to the frequency domain, and then normalize it using the L1 norm to eliminate the difference in energy magnitude between the two, thus obtaining the normalized Fourier amplitude spectra of the two.

6. The method for detecting complex faults in rotating machinery under time-varying operating conditions according to claim 1, characterized in that: In step S4, the method for obtaining the optimal difference spectrum that removes common-mode components and highlights fault characteristics includes: The set of spectral samples corresponding to the normalized Fourier amplitude spectrum of the mixed signal is defined as the mixed signal sample set {NFS}. M The set contains fault characteristic components and normal common-mode components, with a sample size of P; the set of spectral samples corresponding to the normalized Fourier amplitude spectrum of the comparison signal is defined as the comparison signal sample set {NFS}. R This set contains only normal common-mode components, and the number of samples is Q; The two sample sets are fused to obtain the total spectrum sample set, with a total number of samples. And assign a logical label to each sample. When the i-th sample belongs to the set of contrast signals {NFS} R When}, mark When the i-th sample belongs to the mixed signal sample set {NFS} M When}, mark ; Define the feature vector of the i-th sample. This vector is composed of the amplitude sequence of the corresponding spectrum sample and is used to characterize the frequency domain features of the sample; A convex optimization model based on a binary logistic regression framework is constructed. The goal is to solve for the optimal parameter vector by combining maximum likelihood estimation with sparse regularization constraints. The formula is as follows: ， in, The target loss function consists of a regularization term and a likelihood loss term, and is used to measure the deviation between the model's predicted values ​​and the true labels. Let be the parameter vector to be optimized. , For the weight vector, For bias; is the regularization coefficient, a non-negative constant used to control the sparsity of the weight vector; For the target loss function Perform iterative optimization until the difference in the loss function between two consecutive iterations is less than a preset threshold, at which point the optimal parameter vector is obtained. ,in This is the optimal weight vector; Based on the optimal weight vector The optimal difference spectrum is calculated and expressed as: ， in, and These are the optimal weight vectors. The minimum and maximum values ​​in the array.

7. The method for detecting combined faults in rotating machinery under time-varying operating conditions according to claim 1, characterized in that: In step S5, the method for calculating the spectral segmentation threshold for distinguishing between the fault principal component and the noise component in the difference spectrum based on the amplitude percentile sequence of the optimal difference spectrum includes: Extract all amplitude data of the optimal difference spectrum SR to form an amplitude set. K represents the total number of spectral lines in the optimal difference spectrum. Let the amplitude of the k-th spectral line be denoted by ; sort the set of amplitudes in ascending order to obtain a sorted amplitude sequence; based on the sorted amplitude sequence, calculate its amplitude percentile sequence. , This represents the nth percentile amplitude of the sorted amplitude sequence. The length of the amplitude percentile sequence; Construct based on candidate split points The piecewise linear fitting model is used to fit the amplitude percentile sequence piecewise. The model expression is: ， Where T is the piecewise linear fitting model; For sequence index variables; , and , These are the dividing points. The linear fitting coefficients of the two consecutive sequences; For each candidate split point The least squares method was used to solve for the linear fitting coefficients of the two segments separately. and And define a piecewise fitting error function to quantify the fitting accuracy: ， ， in, Candidate split points The fitting error of the corresponding first half of the sequence, Candidate split points The fitting error of the corresponding second half of the sequence; Define the total fitting error function Iterate through all candidate split points Calculate the total fitting error for each candidate point; find the candidate segmentation point that minimizes the total fitting error, which is the optimal segmentation point. The expression is: ， Based on the optimal segmentation point Extract the corresponding index from the percentile sequence of magnitude. elements This element is used as the spectrum segmentation threshold.

8. The method for detecting composite faults in rotating machinery under time-varying operating conditions according to claim 7, characterized in that: In step S6, the method for outputting the composite fault diagnosis result includes: Based on spectrum segmentation threshold A frequency domain separation mask is constructed to filter the spectral components belonging to the target fault in the optimal difference spectrum, retaining only the spectral amplitude of the spectral line index within the fault principal component interval, and suppressing the spectral amplitude of noise and minor components. The frequency domain separation mask is multiplied element-wise with the optimal difference spectrum to extract the spectral components of the target fault. , is represented as: ， in, For spectral line indexing; For the spectral components Performing an inverse fast Fourier transform converts the frequency domain signal back to the time domain, reconstructing the independent fault component signal. , is represented as: ， in, This represents the operator for taking the real part; This represents the inverse fast Fourier transform operator; For the fault component signal Perform order spectrum analysis to extract the order of its fault characteristics. To obtain the theoretical fault characteristic order of rotating machinery , The extracted feature order is compared with the theoretical order using a relative error criterion to determine whether a corresponding fault exists. The criterion formula is as follows: ， in, For harmonic order, The allowable relative error threshold; Repeat the above steps to extract, reconstruct, and compare the order of all target fault spectral components in the optimal difference spectrum, and finally output a composite fault diagnosis result containing all fault types.

9. A composite fault detection system for rotating machinery under time-varying operating conditions, characterized in that, Includes the following modules: The data acquisition and time-frequency enhancement module is used to acquire the raw vibration signal under variable speed conditions, perform a short-time Fourier transform on it, and based on... The weighted time-frequency distribution is obtained by combining the Rényi entropy of order of meshing with the prior of meshing order. The rotational speed sensorless ridge extraction module is used to extract instantaneous frequency ridges that are energy-concentrated, frequency-continuous, and conform to the laws of mechanical motion by constructing a global optimization model based on the weighted time-frequency distribution and by fusing energy statistical prior constraints and physical kinematic prior constraints. The angular domain resampling and spectrum calculation module is used to perform equal-angle resampling on the original vibration signal based on the instantaneous frequency ridge to generate an angular domain stationary signal; the angular domain stationary signal is defined as a hybrid signal; an angular domain stationary signal of the same model and operating condition as the rotating machinery under fault-free operation is selected as a comparison signal, and fast Fourier transform is performed on the hybrid signal and the comparison signal respectively to calculate their normalized Fourier amplitude spectra; The convex optimization difference spectrum generation module is used to take the normalized Fourier amplitude spectrum as input, construct an extended difference mode decomposition convex optimization model, and obtain the optimal difference spectrum that removes common mode components and highlights fault characteristics by solving the weight vector through maximum likelihood and sparse regularization iteration. An adaptive threshold decision module is used to calculate a spectrum segmentation threshold for distinguishing fault principal components and noise components in the difference spectrum based on the amplitude percentile sequence of the optimal difference spectrum using change point analysis. The signal decoupling and diagnostic output module is used to construct a frequency domain mask based on the spectrum segmentation threshold, extract the target fault spectrum component in the optimal difference spectrum, perform inverse fast Fourier transform on the extracted fault spectrum component to reconstruct each independent time-domain fault component, compare the time-domain fault component with the theoretical fault feature order, and output the composite fault diagnosis result.

10. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that: When the processor executes the computer program, it implements the steps of the method as described in any one of claims 1 to 8.

Images