SAH and SSET combined harmonic reducer fault diagnosis method and system under time-varying rotating speed

Abstract

The invention provides a harmonic reducer fault diagnosis method and system combining SAH with SSET under a time-varying rotating speed, and belongs to the field of harmonic reducer fault diagnosis. The invention aims to solve the problems of large time-frequency ridge line extraction error and inaccurate fault feature order extraction caused by the non-stable and signal modulation characteristics of a harmonic reducer vibration signal under a time-varying rotating speed. According to the SAH construction method based on sub-band rearrangement dual-tree complex wavelet packet transformation, a band elimination filter is designed according to an optimal frequency band, and a resonance frequency band where a fault impact component is located is filtered out; a time-frequency ridge extraction method combining SAH-SSET with a fast path optimization algorithm is provided, a frequency ambiguity phenomenon is improved, order analysis is carried out based on the frequency ambiguity phenomenon, and then the fault type of the harmonic reducer is judged according to an envelope order spectrum.

Description

Technical Field

[0001] The present invention relates to the technical field of harmonic reducer fault diagnosis, and in particular to a harmonic reducer fault diagnosis method and system combining SAH with SSET under time-varying speed. Background Art

[0002] Industrial robots have been widely used in automated production due to their high precision and low operating costs. [1,2] The harmonic reducer has the characteristics of compact structure, light weight, large reduction ratio and precise transmission. It is the core transmission component of industrial robots. Its normal operation directly determines whether the industrial robot can work safely, accurately and reliably. Once the harmonic reducer fails, it will cause the equipment to stop working at the least, and even endanger personal safety at the worst. [3,4] In actual application scenarios, harmonic reducers often work under continuously changing working conditions. It is of great significance to conduct in-depth research on the fault diagnosis of harmonic reducers under time-varying speed. [5,6] .

[0003] In actual operation, harmonic reducers often accelerate and decelerate frequently, and the speed fluctuations make traditional signal processing methods ineffective. Order analysis is the most direct and effective method to solve the time-varying speed problem. [7] Many scholars use order analysis to diagnose faults under time-varying speeds. At present, there are few literatures on the fault diagnosis of harmonic reducers under time-varying speeds at home and abroad. Bearings, gearboxes and harmonic reducers are all rotating elements. This paper takes the fault diagnosis method of bearings and gearboxes under time-varying speeds as an example to explain the research status of order analysis in the field of time-varying speed fault diagnosis. Reference [8] proposed a gearbox rolling bearing fault diagnosis method combining fast spectral kurtosis with order analysis, and verified the effectiveness of the proposed method using the fault signal of the outer ring of the gearbox rolling bearing. Reference [9] proposed a fault diagnosis method based on local mean decomposition of empirical optimal envelope combined with order analysis, which eliminated the influence of speed fluctuation on fault feature extraction and achieved good results when the speed change was small. Reference

[10] introduced the coherent resonance theory, judged the bearing fault mode by the resonance factor index of the response order spectrum, and realized the rolling bearing fault diagnosis under time-varying speed.

[0004] Traditional order analysis methods provide a way to realize fault diagnosis of rotating components under time-varying speed, but they require additional sensors to provide the instantaneous frequency (IF) of the shaft, such as tachometers or encoders.

[11] However, traditional order analysis is greatly affected by the key phase pulse resolution, has poor performance when the speed changes sharply, and has complex assembly conditions.

[12] , which limits its application in actual industrial scenarios. To address the above problems, tachometer-free order analysis methods have been rapidly developed, and extracting time-frequency ridges from time-frequency spectra as the shaft IF has gradually become a research hotspot.

[0005] To obtain the time-frequency ridge, we first need to generate a time-frequency spectrum with clear time-frequency ridges. The time-frequency spectrum is based on the time-frequency representation (Time Frequency R e It is a visualization form of time-frequency representation. In order to solve the problem that the time-frequency ridges in the time-frequency spectrum of the traditional time-frequency analysis (TFA) method are often affected by noise, many scholars have used the advanced TFA method to highlight these ridges. Reference

[13] extracted the time-frequency ridges as the shaft IF from the time-frequency spectrum obtained by synchronous compression wavelet transform, and realized the extraction of the fault characteristics of the inner ring of the bearing with a speed change of about 300r / min within 2s. Reference

[14] obtained the logarithmically rearranged time-frequency spectrum through Wigner-Ville transform, and used the Crazy climber method to extract the time-frequency ridges as the estimated shaft IF. The effectiveness of the method was verified by the local wear fault signal of the gear. Reference

[15] extracted the time-frequency ridges as the shaft IF based on the time-frequency spectrum of Fourier synchronous compression transform, and realized the diagnosis of bearing faults under time-varying speed.

[0006] Some scholars hope to further reduce the time-frequency ridge error by optimizing algorithms or constructing operators. Reference

[15] introduced the scale space improved generalized linear wavelet transform (GLCT) for the case of monotonic speed variation trend, extracted the time-frequency ridge from its TFR as the shaft IF, and realized the tachometer-free order analysis of the bearing under time-varying speed. Reference

[16] introduced the gray wolf optimizer and the Gini index to search for the optimal parameters of GLCT for the non-monotonic speed problem, realized the non-monotonic time-frequency ridge extraction, and realized the bearing fault diagnosis under time-varying speed by order analysis based on the extracted time-frequency ridge. Reference

[17] proposed a parameterized multi-synchronous compression transformation method based on weighted least squares method, which improved the energy aggregation of TFR and verified the effectiveness of the method using the inner race fault signal and outer race fault signal of the bearing. Reference

[18] combines the power spectrum-based autoregressive variational mode decomposition with the window width optimization algorithm, and proposes an improved adaptive window width local maximum synchronous compression transform, which obtains a TFR with high energy concentration and improves the accuracy of IF. Reference

[19] proposes a generalized synchronous compression transform, constructs a generalized synchronous compression operator through iterative operation, and rearranges the time-frequency energy to obtain a more concentrated TFR. The effectiveness of the method is verified by the bearing inner ring fault signal. Some other scholars focus on signal preprocessing methods to reduce the time-frequency ridge error. Reference

[20] proposes a variational mode decomposition (VMD) combined with SET tachometer-free order tracking method to achieve fault diagnosis of wind turbine gearbox under time-varying conditions. Reference

[21] separates frequency components based on spectrum amplitude modulation, selects appropriate weights to highlight the time-frequency ridge, and uses the peak search algorithm to extract the time-frequency ridge as the shaft IF. The effectiveness of the method is verified by the bearing outer and inner ring fault signals.

[0007] In summary, extracting time-frequency ridges as the shaft interferometer (IF) is the key to solving fault diagnosis problems under time-varying speeds using tachometer-free order analysis. Harmonic reducer vibration signals under time-varying speeds exhibit nonstationary and modulated characteristics, leading to large errors in time-frequency ridge extraction and inaccurate fault feature order extraction. To address these issues, a harmonic reducer fault diagnosis method under time-varying speeds is proposed, combining the subband averaging Hoyergram (SAH) with the second-order synchronous extraction transform (SSET). This method uses the SAH to determine the optimal frequency band for band-stop filtering and combines it with the SSET to generate a time-frequency spectrum with clear time-frequency ridges, thus improving frequency ambiguity. Fast path optimization is used to extract time-frequency ridges from the time-frequency frame (TFR) as the shaft interferometer (IF). A Hilbert transform is performed on the original signal to obtain an envelope signal. This envelope signal is then resampled in the angular domain based on the extracted time-frequency ridges to obtain an envelope order spectrum, ultimately enabling effective diagnosis of harmonic reducers under time-varying speeds. Summary of the Invention

[0008] The technical problems to be solved by the present invention are:

[0009] In order to solve the problem that the vibration signal of the harmonic reducer under time-varying speed has the characteristics of non-stationary and signal modulation, which leads to large errors in time-frequency ridge extraction and inaccurate fault feature order extraction.

[0010] The present invention is to solve the above technical problems using the following technical solutions:

[0011] The present invention provides a harmonic reducer fault diagnosis method combining SAH with SSET under time-varying speed, comprising the following steps:

[0012] S100, collecting time domain vibration signals of typical faults of harmonic reducers under time-varying speed;

[0013] S200, using SAH to determine the optimal filtering parameters of the time domain vibration signal of the harmonic reducer collected in step S100, including the center frequency and the optimal bandwidth, and perform band-stop filtering;

[0014] S300, performing Hilbert transform on the filtered signal and using the envelope time-frequency diagram obtained by SSET, then performing fast path optimization and extracting the time-frequency ridge as the rotation axis IF;

[0015] S400, performing Hilbert transform on the original harmonic reducer vibration signal and obtaining its envelope, then performing angular domain resampling on the envelope signal according to the time-frequency ridge extracted in step S300 to obtain an angular domain envelope signal, and performing fast Fourier transform to obtain an envelope order spectrum;

[0016] S500. Determine the fault type by comparing the theoretical fault characteristic order of the harmonic reducer with the envelope order spectrum.

[0017] Furthermore, in step S200, it includes:

[0018] S210, using the wavelet coefficients to construct the subband spectral kurtosis, as shown in formula (1), the spectral kurtosis is the square of the spectral L2 / L1 norm, and the Hoyer index is the normalized form of the L2 / L1 norm, as shown in formula (2),

[0019]

[0020] Where a is the L2 / L1 norm, N is the signal length, n represents the discrete time variable, K k,r is the spectral kurtosis, H k,r is the Hoyer index, d k,r is the sub-band wavelet coefficient rearranged in descending order of center frequency, k is the number of decomposition layers, and each layer divides the center frequency into 2 k intervals, r is one of the intervals;

[0021] S210, the step of constructing SAH includes:

[0022] S211, using a sliding window, evenly divide the collected harmonic reducer vibration signal x(t) into M sub-signals {x m (t)|m=1,2,...,M}, the window length of the sliding window is N x / M, the step size is equal to the window length, N x is the length of signal x(t);

[0023] S212, using sub-band rearrangement DWCWPT to decompose each sub-signal into a set of sub-band wavelet coefficients

[0024] S213, calculating the Hoyer index of the subband;

[0025] S214, average the Hoyer index of the corresponding sub-band to obtain the sub-band average Hoyer index As shown in formula (3):

[0026]

[0027] S215, draw SAH, select the best frequency band corresponding to the maximum sub-band average Hoyer index (f c ,B ω ), center frequency f c and the optimal bandwidth B ω , which is the optimal filtering center frequency and bandwidth.

[0028] Furthermore, in step S200, it also includes:

[0029] When using SAH to obtain the resonance frequency band (f c ,B ω ), the vibration signal of the harmonic reducer is band-stop filtered to obtain the vibration signal x1(t). The STFT expression of x1(t) based on the window function g(t) is:

[0030]

[0031] Where τ is the time-shift variable, * is the conjugate operation, g(t) is the window function, g(τ-t) is the sliding window function, i is the imaginary unit, ω is the normalized angular frequency of the input signal, and t represents the continuous time variable;

[0032] The Gaussian window function with standard deviation σ is expressed as:

[0033]

[0034] Furthermore, in step S300, it includes:

[0035] Definition Multi-component AM / FM signal is defined as:

[0036]

[0037] Among them, K is the number of modes, A k (t),φ k (t) represents the instantaneous amplitude and instantaneous frequency of each mode respectively;

[0038] The ideal time-frequency representation of x1(t) is shown in formula (7):

[0039]

[0040] Where δ(·) is the delta function, which represents the impulse response; φ′ k (t) is the ideal instantaneous frequency;

[0041] Assume that the single-component signal is locally approximated by a Gaussian modulated linear chirp signal at every moment:

[0042]

[0043] in, represents the instantaneous amplitude of x1(t); c0 is the initial phase, c1 is the initial frequency, c2 is the modulation frequency, and t0 is the envelope center time;

[0044] According to formula (8), the derivative of x1(t) with respect to t is expressed as:

[0045]

[0046] Among them, the imaginary part is the signal IF; q x1 、p x1 is the composite parameter obtained from derivative analysis, A(t) is the amplitude;

[0047] Derived The partial derivative with respect to t is:

[0048]

[0049] Then deduce The partial derivative with respect to t is:

[0050]

[0051] in, represents the STFTs obtained using g″(t),tg′(t),g′(t) as window functions, which can be obtained by using formulas (10) and (11):

[0052]

[0053] φ′(t,ω) is expressed as formula (13):

[0054]

[0055] in, To take the imaginary part;

[0056] The TFR of SSET is defined as formula (14):

[0057]

[0058] Furthermore, in step S300, the method further includes extracting a time-frequency ridge based on fast path optimization, using a support function to optimally extract the time-frequency ridge from the TFR with minimal frequency jump;

[0059] Assume that the TFR of signal x1(t) is SSET X1(t,τ), where τ is the time variable and f is the frequency variable, and the time-frequency ridge f is extracted from X1(τ,f) using fast path optimization p (τ), and the upper boundary of the time-frequency ridge f up (τ) and the lower bound f down (τ);

[0060] The time τ in TFR n The number of peaks is recorded as N p (τ n), the frequency corresponding to the m-th peak is recorded as V m (τ n ), the TFR amplitude of the mth peak is recorded as Q m (τ n ); τ n Peak value at time v m (τ n ) is defined by formula (15):

[0061]

[0062] If X1(τ,f) has a time span [τ1, τ2, τ3, …, τ n ], path optimization is described as:

[0063]

[0064] Among them, m c (τ n ) is to determine the peak to be extracted as τ n The ridge at the position, F[·] is the selected optimization support function, {m1,m2,m3,....,m N} represents the sequence of peak numbers along the time span. Fast path optimization only relies on a limited number of preceding points, and the support function is:

[0065]

[0066] Among them, f d (t n-1 ) is the candidate ridge point at τ n-1 The frequency, f d is the frequency of a series of candidate ridge points in history [τ1, τ2, τ3, ..., τ n-1 ], △f d f d The derivative of the series, m[·] is the median of the series, IQR is the interquartile range, λ1 and λ2 are penalty factors, and the weight functions ω1(·) and ω2(·) are used to suppress the atypical changes of the ridge frequency value and the derivative, respectively.

[0067] Furthermore, in step S400, it includes:

[0068] The fault characteristic order of the inner and outer rings of the flexible bearing is derived through the calculation formula of the bearing fault characteristic frequency:

[0069]

[0070] Among them, O i , O o are the orders of the inner and outer rings of the flexible bearing of the harmonic reducer, f i ,f oare the fault characteristic frequencies of the inner and outer rings of the flexible bearing of the harmonic reducer, f r is the rotation frequency, Z is the number of rolling elements, D is the pitch diameter, d is the ball pitch diameter, and α is the contact angle;

[0071] It is inferred that the flexspline and the rigid wheel have a common fault characteristic frequency of 2 times the rotation frequency, and the flexspline also has a fault characteristic frequency caused by self-rotation (2 / Z F ) times the rotation frequency, Z F is the number of flexspline teeth. According to the relationship between the fault characteristic order and the fault characteristic frequency, formulas (20) and (21) are obtained:

[0072] O CMT =O FMT1 =f CMT,FMT1 / f r =2f r / f r =2 (20)

[0073] O FMT2 =f FMT2 / f r =(2 / Z F )*f r / f r =2 / Z F (twenty one)

[0074] Among them, O CMT and O FMT are the orders of the rigid wheel and flexible wheel of the harmonic reducer, O FMT 1 is one of the fault characteristic orders of the flexible wheel, f CMT,FMT1 is the fault characteristic frequency of the harmonic reducer rigid wheel and flexible wheel, O FMT2 is the flexspline rotation order of the harmonic reducer, f FMT2 It is the fault characteristic frequency caused by the self-rotation of the flexspline of the harmonic reducer.

[0075] Furthermore, in step S500, the flexspline and rigid-wheel faults are distinguished from each other in terms of amplitude. When the flexspline fault occurs, the amplitude is relatively obvious due to the change in the fault point position, and abundant sidebands appear around it. When the rigid-wheel fault occurs, the fault position remains unchanged, the amplitude is relatively small, and the sidebands are relatively few. Accordingly, in the order spectrum, the characteristic order amplitude of the flexspline fault is larger than that of the rigid-wheel fault.

[0076] A harmonic reducer fault diagnosis system combining SAH with SSET under time-varying speed, the system has a program module corresponding to the above steps, and executes the steps in the harmonic reducer fault diagnosis method combining SAH with SSET under time-varying speed when running.

[0077] A computer-readable storage medium stores a computer program, wherein the computer program is configured to implement the steps of a harmonic reducer fault diagnosis method combining SAH with SSET under time-varying speed when called by a processor.

[0078] Compared with the prior art, the present invention has the following beneficial effects:

[0079] (1) This paper proposes a method for constructing a SAH based on DTCWPT to determine the parameters of the band-stop filter and filter out the resonant frequency band of the fault impact component. When processing harmonic reducer fault signals under time-varying speeds, the fast spectral kurtosis diagram cannot obtain the filter parameters, while the SAH can determine the frequency band of the band-stop filter. The Hoyer index in the SAH is between 0.05 and 0.3, while the Hoyer index of the Hoyer diagram based on DTCWPT is between 0.0085 and 0.0115. In comparison, the SAH has a larger Hoyer index, and the differences between different frequency bands are clearer, resulting in higher frequency resolution.

[0080] (2) The present invention proposes a time-frequency ridge extraction method combining SAH-SSET with fast path optimization. SAH-SSET can effectively highlight the time-frequency ridge, improve the frequency ambiguity, and reduce the time-frequency ridge error. Fast path optimization is combined with PSO-VMD-SSET and ICEEMDAN-SSET as comparison methods, and the error of the time-frequency ridge is evaluated using quantitative indicators. In comparison, the time-frequency ridge extracted by SAH-SSET combined with fast path optimization has an RMSE of 0.0072, which is closer to 0, and the determination coefficient R is 0.0072. 2 It is 0.8960, which is closer to 1, indicating that the average deviation between the time-frequency ridge and the actual rotation frequency is small, and the time-frequency ridge can better reflect the changing trend of the actual rotation frequency.

[0081] (3) This paper studies the operating status of harmonic reducers, analyzes the theoretical calculation formula for the fault characteristic order of the flexspline and rigid pulley, and proposes a method for diagnosing harmonic reducer faults under time-varying speeds. This method verifies the theory of the fault characteristic order of the flexspline and rigid pulley of harmonic reducers. The maximum error between the experimental fault characteristic order and the theoretical fault characteristic order is 1.59%. The fault characteristic order in the envelope order spectrum can be used to accurately identify the type of harmonic reducer fault, verifying the effectiveness of the proposed method. BRIEF DESCRIPTION OF THE DRAWINGS

[0082] Figure 1 A flowchart of constructing a SAH in an embodiment of the present invention;

[0083] Figure 2 This is a structural diagram of a harmonic reducer in an embodiment of the present invention;

[0084] Figure 3Flowchart of a harmonic reducer fault diagnosis method combining SAH with SSET under time-varying speed according to an embodiment of the present invention;

[0085] Figure 4 This is a structural diagram of the harmonic reducer vibration signal acquisition experimental platform in an embodiment of the present invention;

[0086] Figure 5 This is a structural diagram of the measurement and control system of the harmonic reducer vibration signal acquisition experimental platform in an embodiment of the present invention;

[0087] Figure 6 This is a diagram showing the deployment of sensors on the harmonic reducer vibration signal acquisition experimental platform in an embodiment of the present invention;

[0088] Figure 7 This is a diagram of the time-varying speed setting interface in an embodiment of the present invention;

[0089] Figure 8 This is a square wave signal diagram of the speed and torque sensor in an embodiment of the present invention;

[0090] Figure 9 This is a graph showing the input shaft rotation frequency in an embodiment of the present invention;

[0091] Figure 10 This is a diagram of an inner race fault signal in an embodiment of the present invention;

[0092] Figure 11 : This is a SAH result diagram of the inner race fault signal in an embodiment of the present invention;

[0093] Figure 12 The fast spectral kurtosis graph and Hoyer graph in the embodiment of the present invention;

[0094] Figure 13 Spectrogram of the envelope of STFT, SSET, SAH-STFT and SAH-SSET in the embodiment of the present invention;

[0095] Figure 14 A comparison diagram of the sensor rotation frequency and time-frequency ridge line in an embodiment of the present invention;

[0096] Figure 15 This is a hyperparameter optimization curve diagram in an embodiment of the present invention;

[0097] Figure 16 This is a diagram of the PSO-VMD decomposition results in an embodiment of the present invention;

[0098] Figure 17 This is the ICEEMDAN decomposition result diagram in an embodiment of the present invention;

[0099] Figure 18 The time-frequency diagram of PSO-VMD-SSET and the comparison diagram of time-frequency ridge and sensor rotation frequency in the embodiment of the present invention are shown;

[0100] Figure 19 The time-frequency diagram of ICEEMDAN-SSET and the comparison diagram of time-frequency ridges and sensor rotation frequency in an embodiment of the present invention are shown;

[0101] Figure 20 This is the envelope order spectrum of the inner race fault in the embodiment of the present invention;

[0102] Figure 21 The fault envelope order spectrum diagrams in the embodiment of the present invention include the outer ring fault envelope order spectrum diagram, the flexspline tooth missing envelope order spectrum diagram, and the rigid wheel tooth missing envelope order spectrum diagram;

[0103] Figure 22 This is the envelope order spectrum of PSO-VMD-SSET in the embodiment of the present invention;

[0104] Figure 23 This is the envelope order spectrum of ICEEMDAN-SSET in an embodiment of the present invention. DETAILED DESCRIPTION

[0105] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, specific embodiments of the present invention are described in detail below with reference to the accompanying drawings.

[0106] 1 Subband Average Hoyer Plot Based on Dual-Tree Complex Wavelet Packet

[0107] 1.1 Dual-tree complex wavelet transform

[0108] Dual-tree Complex Wavelet Transform (DTCWT) has the advantages of approximate shift invariance, reduced aliasing and low computational cost.

[22] The DTCWT operates a pair of real filter bank trees in parallel on the input signal to form the real and imaginary parts of the complex wavelet transform. The wavelets corresponding to the real and imaginary trees are denoted as ψ Re (t) and ψ Im (t), satisfying the scaling relationship of formulas (1) and (2):

[0109]

[0110] in, * represents the real or imaginary part, and are the filter banks of real trees and imaginary trees (including low-pass filters and high-pass filters), t represents the continuous time variable, and n represents the discrete time variable. DTCWT ensures that the wavelet pair satisfies the approximate Hilbert transform relationship, as shown in formula (3):

[0111] ψ Im (t)≈H{ψRe (t)} (3)

[0112] H{·} is the Hilbert transform operation, so that the complex wavelet ψ(t) = ψ Re (t)+iψ Im (t), where i represents the imaginary unit; if the two real discrete wavelet transforms are orthogonal, then the necessary and sufficient condition of formula (3) is:

[0113]

[0114] In the frequency domain, the condition is equivalent to:

[0115]

[0116] Z′(·) is the z-transform operation, ω is the normalized angular frequency of the input signal. In the DTCWT theory, the conditions given by formulas (5)-(7) are called half-sample delay conditions. The key to DTCWT is to design a filter bank pair that meets the following conditions: (1) approximate half-sample delay characteristics; (2) orthogonal or biorthogonal; (3) compact support; (4) moment vanishing; and (5) linear phase.

[0117] 1.2 Dual-tree complex wavelet packet transform with subband rearrangement

[0118] Dual-tree Complex Wavelet Packet Transform (DTCWPT) is a grouping form of DTCWT that can achieve finer frequency resolution in the high-frequency region.

[22] To implement DTCWPT, the filter bank tree consists of three types of filters:

[0119] (1) First-stage filter and

[0120] (2) Double-tree filter and

[0121] (3) General filter {f0(n), f1(n)}.

[0122] Compared to DTCWT, DTCWPT adds only one extra rule: whether in the real tree or the imaginary tree, it is necessary to ensure that the same filter bank {f0(n), f1(n)} is applied to its own output and the high-pass output of the dual-tree filter.

[0123] DTCWPT decomposes the signal into multiple frequency sub-bands. Due to the bandpass characteristics of the wavelet filter bank in some frequency bands, the center frequencies of the sub-bands are disordered, which is called sub-band disorder. Therefore, the sub-bands are rearranged after the sub-bands are directly obtained, and rearranged in descending order according to their center frequencies, which is recorded as {d k,r |r=1,2,...,2 k}, d k,r The sub-band wavelet coefficients are rearranged in descending order of center frequency, k represents the number of decomposition layers, and each layer divides the frequency into 2 k intervals, r is one of them.

[0124] 1.3SAH construction process

[0125] Reference

[22] uses wavelet coefficients to construct the subband spectral kurtosis, as shown in formula (8). In reference

[23] , spectral kurtosis is interpreted as the square of the spectral L2 / L1 norm, and the Hoyer index is the normalized form of the L2 / L1 norm. Specifically, as shown in formula (9):

[0126]

[0127] Among them, a is the L2 / L1 norm, N is the signal length, K k,r is the spectral kurtosis, H k,r is the Hoyer index.

[0128] The specific process of constructing SAH using sub-band rearranged dual-tree complex wavelet packet transform is as follows: Figure 1 shown. Figure 1 The main steps of constructing SAH are as follows:

[0129] (1) The collected harmonic reducer vibration signal x(t) is evenly divided into M sub-signals {x m (t)|m=1,2,...,M}, the window length of the sliding window is N x / M, the step size is equal to the window length, N x is the length of signal x(t);

[0130] (2) Use subband rearrangement DWCWPT to decompose each sub-signal into a set of sub-band wavelet coefficients

[0131] (3) Calculate the Hoyer index of the subband;

[0132] (4) Calculate the average value of the Hoyer index of the corresponding sub-band to obtain the sub-band average Hoyer index As shown in formula (10):

[0133]

[0134] (5) Draw SAH and select the best frequency band (f c ,B ω ), center frequency f c and the optimal bandwidth B ω , which is the optimal filtering center frequency and bandwidth.

[0135] 2 Second-order synchronous extraction transform

[0136] When using SAH to obtain the resonance frequency band (f c ,B ω ), the vibration signal of the harmonic reducer is band-stop filtered to obtain the vibration signal x1(t). The STFT (Short Time Fourier Transform) of x1(t) based on the window function g(t) is expressed as:

[0137]

[0138] Where τ is the time-shift variable, * is the conjugate operation, g(t) is the window function, g(τ-t) is the sliding window function; i is the imaginary unit;

[0139] Among them, the Gaussian window function with standard deviation s is expressed as:

[0140]

[0141] Among them, σ represents the standard deviation. By changing σ, the width and decay speed of the window can be controlled.

[0142] However, since the time-frequency ridge always has a frequency bandwidth in the TFR of STFT, STFT cannot present the ideal time-frequency characteristics. The multi-component AM / FM signal is defined as:

[0143]

[0144] Among them, K is the number of modes, A k (t),φ k (t) represents the instantaneous amplitude and instantaneous frequency of each mode respectively.

[0145] The ideal time-frequency representation (ITFR) of x1(t) is shown in formula (14):

[0146]

[0147] Where δ(·) is the delta function, which represents the impulse response; φ′ k (t) is the ideal instantaneous frequency;

[0148] Considering the universality, it is assumed that the single component signal is locally approximated by a Gaussian modulated linear chirp signal at every moment.

[24] :

[0149]

[0150] in, represents the instantaneous amplitude of x1(t); c0 is the initial phase, c1 is the initial frequency, c2 is the modulation frequency, and t0 is the envelope center time.

[0151] According to formula (15), the derivative of x1(t) with respect to t can be expressed as:

[0152]

[0153] Among them, the imaginary part is the signal IF, q x1 、p x1 is the composite parameter obtained from derivative analysis, A(t) is the amplitude, and p and q are complexes obtained from derivative analysis.

[0154] It can be deduced The partial derivative with respect to t is:

[0155]

[0156] Similarly, it can be deduced that The partial derivative with respect to t is:

[0157]

[0158] in, represents the STFTs obtained using g″(t),tg′(t),g′(t) (second-order Gaussian window function; first-order derivative-weighted Gaussian window function; first-order Gaussian window function can also be called the first-order derivative of the Gaussian window function) as the window function. According to formulas (17) and (18), we can get:

[0159]

[0160] So φ′(t,ω) can be expressed as formula (20):

[0161]

[0162] in, To take the imaginary part;

[0163] The TFR of SSET is defined as formula (21)

[24] :

[0164]

[0165] 3 Time-frequency ridge extraction based on fast path optimization

[0166] The fast path optimization method uses support functions to optimally extract time-frequency ridges from TFR with minimal frequency jumps.

[25] Assume that the TFR of signal x1(t) is SSET X1(t,τ), where τ is the time variable and f is the frequency variable, then the time-frequency ridge f can be extracted from X1(τ,f) using fast path optimization p (τ), and the upper boundary of the time-frequency ridge f up (τ) and the lower bound f down (τ).

[0167] The time τ in TFR n The number of peaks is recorded as N p (τ n ), the frequency corresponding to the m-th peak is recorded as V m (τ n ), the TFR amplitude of the mth peak is recorded as Q m (τ n ). τ n Peak value at time v m (τ n ) can be defined by formula (22):

[0168]

[0169] If X1(τ,f) has a time span [τ1, τ2, τ3, …, τ n ], path optimization can be described as:

[0170]

[0171] Among them, Q m (τ n ) is the TFR amplitude of the mth peak, is τ n Peak value at time, m c (τ n ) is to determine the peak to be extracted as τ n The ridge at the position, F[·] is the selected optimization support function, {m1,m2,m3,....,m N} represents the sequence of peak numbers along the time span. Fast path optimization only relies on a limited number of preceding points, and the support function is:

[0172]

[0173] Among them, f d (t n-1 ) is the candidate ridge point at τ n-1The frequency, f d is the frequency of a series of candidate ridge points in history [τ1, τ2, τ3, ..., τ n-1 ], △f d f d The derivative of the series is m[·], m[·] is the median of the series, IQR is the interquartile range, and λ1 and λ2 are penalty factors. The weight functions ω1(·) and ω2(·) are used to suppress atypical changes in the ridge frequency value and derivative, respectively.

[0174] 4. Harmonic reducer theoretical fault characteristic order

[0175] Typical fault types of harmonic reducers include flexible bearing inner ring, flexible bearing outer ring, rigid wheel missing teeth, and flexible wheel missing teeth. The main components of the harmonic reducer are wave generator, flexible wheel and rigid wheel. The wave generator consists of flexible bearing and cam. The structure is as follows: Figure 2 shown.

[0176] 4.1 Flexible bearing fault characteristic order

[0177] The fault characteristic frequency of the flexible bearing can refer to the calculation formula of the fault characteristic frequency of the ordinary bearing. The calculated fault characteristic frequency fluctuates within the range of 1% to 2% relative to the average value. Therefore, the fault characteristic order of the inner and outer rings of the flexible bearing can be deduced. [26,27] :

[0178]

[0179] Among them, O i , O o are the orders of the inner and outer rings of the flexible bearing of the harmonic reducer, f i ,f o are the fault characteristic frequencies of the inner and outer rings of the flexible bearing of the harmonic reducer, f r is the rotation frequency, Z is the number of rolling elements, D is the pitch diameter, d is the ball pitch diameter, and α is the contact angle.

[0180] 4.2 Fault characteristic order of flexible pulley and rigid pulley

[0181] According to the operating principle of the harmonic reducer, when there is a fault in the rigid wheel, when the wave generator rotates one circle, the passing frequency of the rigid wheel fault point is 2, which is also the passing frequency of the flexspline fault point. At the same time, the flexspline itself has a rotation frequency. When the wave generator rotates one circle, the flexspline rotates 2 teeth relative to the rigid wheel. Therefore, the passing frequency of the second fault point of the flexspline is (2 / Z F ), Z F = is the number of teeth of the flexspline. Therefore, it can be inferred that the flexspline and the rigid wheel have the same fault characteristic frequency of 2 times the rotation frequency. In addition, the flexspline also has a fault characteristic frequency caused by self-rotation of (2 / Z F) times the rotation frequency, according to the relationship between the fault characteristic order and the fault characteristic frequency, we can get formulas (27) and (28):

[0182] O CMT =O FMT1 =f CMT,FMT1 / f r =2f r / f r =2 (27)

[0183] O FMT2 =f FMT2 / f r =(2 / Z F )*f r / f r =2 / Z F (28)

[0184] Among them, O CMT and O FMT are the orders of the rigid wheel and flexible wheel of the harmonic reducer, O FMT 1 is the fault characteristic order of the flexible wheel

[0185] 1,f CMT,FMT1 is the fault characteristic frequency of the harmonic reducer rigid wheel and flexible wheel, O FMT2 is the flexspline rotation order of the harmonic reducer, f FMT2 It is the fault characteristic frequency caused by the self-rotation of the flexspline of the harmonic reducer.

[0186] Since the harmonic transmission ratio is equal to the number of teeth at the output end (flexspline or rigid wheel) divided by the difference in the number of teeth between the rigid wheel and the flexspline (generally 2), and in general, the commonly used transmission ratios of the flexspline output are 30, 50, 80, 100, 120, and 160. This means that in general, the number of teeth of the flexspline should be at least 60, so O FMT2 With O FMT1 There is a difference of at least 66.67 times between them. In this case, it is often difficult to observe O in the order spectrum. FMT2 Therefore, the fault characteristic order diagram of the flexspline should have only one fault characteristic order O FMT1 and its multiples.

[0187] In summary, theoretically, the order of the fault characteristics of both the flexspline and the rigid pulley in the order spectrum is 2, so it is necessary to distinguish between flexspline and rigid pulley faults in terms of amplitude. During the operation of the harmonic reducer, the flexspline is constantly moving and deforming, while the rigid pulley is fixed. When the flexspline fails, the amplitude is relatively obvious due to the change in the fault point location, and there are abundant side frequencies around it; while when the rigid pulley fails, the fault location remains unchanged, the amplitude is relatively small, and there are fewer side frequencies.

[28] Correspondingly, in the order spectrum, the characteristic order amplitude of the flexspline fault should be larger than that of the rigid-wheel fault.

[0188] 5. Fault diagnosis method and process of harmonic reducer under time-varying speed

[0189] Aiming at the problem that the time-frequency ridge error in the time-frequency spectrum of the harmonic reducer vibration signal under time-varying speed is large, which leads to the inaccurate order of fault feature analysis without tachometer order, a harmonic reducer fault diagnosis method combining SAH with SSET under time-varying speed is proposed. The flow chart is as follows: Figure 3 The specific steps are as follows:

[0190] (1) Collect the time domain vibration signals of typical faults of harmonic reducer under time-varying speed.

[0191] (2) Using SAH to determine the optimal filtering parameters of the time domain vibration signal of the harmonic reducer, including the center frequency f c and the optimal bandwidth B ω , perform band-stop filtering.

[0192] (3) Perform Hilbert transform on the filtered signal and use the envelope time-frequency diagram obtained by SSET, combine it with formula (24) to perform fast path optimization, and extract the time-frequency ridge as the rotation axis IF.

[0193] (4) Performing Hilbert transform on the original harmonic reducer vibration signal and obtaining its envelope, then resampling the envelope signal in the angular domain according to the time-frequency ridge extracted in step (3) to obtain the angular domain envelope signal, and performing Fast Fourier Transform (FFT) to obtain the envelope order spectrum.

[0194] (5) The fault type is determined by comparing the theoretical fault characteristic order of the harmonic reducer with the envelope order spectrum.

[0195] 6 Application and Analysis

[0196] 6.1 Experimental Description

[0197] The proposed method is verified by measuring the actual signal of the harmonic reducer under time-varying speed. The relevant parameters are shown in Table 1. The data is collected using the harmonic reducer vibration signal acquisition experimental platform. The platform is shown in the figure. Figure 4 As shown in Figure 2, speed and torque sensors are installed at both the input and output ends to measure speed and torque. The square wave signal sampling frequency of the speed and torque sensors is 100kHz.

[0198] Table 1 Harmonic reducer related parameters

[0199]

[0200] The measurement and control system of the harmonic reducer vibration signal acquisition experimental platform is as follows: Figure 5Install the acceleration vibration sensor on the harmonic reducer holder. The installation position is as shown in Figure 6 As shown, the vibration signal sampling frequency is 100kHz. The vibration signal collected by the 9 o'clock direction sensor 1 is selected for subsequent experiments.

[0201] In actual application scenarios, it is common to use triangular waves or trapezoidal waves to control the operation of industrial robots. This experiment simulates the trapezoidal wave process through the main control system and sets a trapezoidal speed curve. In addition, referring to the high-load working environment and rapidly changing speed of the harmonic reducer in actual industry, the torque load is set to 10N·m, and the speed is set to accelerate from 400 (r / min) to 1000 (r / min) in a step-by-step manner. After staying at 1000 (r / min) for 0.5 seconds, it is decelerated from 1000 (r / min) to 500 (r / min) in a step-by-step manner. The step-by-step acceleration is set to increase the speed by 100 (r / min) every 100ms, and the step-by-step deceleration is set to decrease the speed by 100 (r / min) every 100ms. The time-varying speed setting interface is as follows: Figure 7 shown.

[0202] In order to simulate the failure of various components of the harmonic reducer, the inner and outer rings of the flexible bearing were processed by laser pitting method respectively, and the rigid wheel tooth missing and flexspline tooth missing fault parts were made by wire cutting process. A total of four types of fault types of harmonic reducers were obtained.

[0203] 6.1.1 Input shaft speed

[0204] The input speed torque sensor is located between the drive motor and the harmonic reducer. The three are connected in pairs with couplings. When the drive motor controls the shaft to rotate, the speed code disk rotates at the same frequency and outputs a square wave signal through the photoelectric switch, such as Figure 8 The number of teeth of the encoder used is 1080, that is, when the speed measuring code disk rotates 1 circle, 1080 square waves will appear.

[0205] First, determine the position of each square wave, and then calculate the time of each square wave based on the sampling rate. Then calculate the time difference between adjacent square waves, and calculate the rotation frequency based on the time difference and the number of encoder teeth. Finally, assume that there is a linear relationship between the rotation frequency of the last pulse and the rotation frequency of the first two square waves, and extrapolate to supplement the rotation frequency at the last square wave to obtain the shaft rotation frequency. The result is as follows Figure 9 shown.

[0206] 6.1.2 Theoretical Fault Characteristic Order

[0207] In this experiment, the input shaft rotation frequency order is set to 1. The inner race fault characteristic order and outer race fault characteristic order can be obtained according to equations (26) and (27) in 3.1 and the relevant parameters of the harmonic reducer in Table 1. The flexspline and rigid spline fault orders are determined according to the theory in 3.2. In summary, the theoretical fault characteristic orders of each component of the harmonic reducer used in the experiment are shown in Table 2.

[0208] Table 2 Theoretical fault characteristic order of harmonic reducer

[0209]

[0210] 6.2 SAH Experimental Analysis

[0211] The signal measured by the harmonic reducer with inner race fault is selected for analysis. In order to reduce the computational complexity, the collected vibration signal is downsampled by 8 times. Figure 10 As shown, the sampling frequency f after downsampling s The signal is subjected to SAH, with the number of sub-signals M set to 50 and the decomposition level k set to 5. The result is as follows: Figure 11 As shown, the horizontal axis is the frequency interval, the vertical axis is the level k, and the color represents the sub-band average Hoyer index size.

[0212] In order to verify the effectiveness of SAH based on dual-tree complex wavelet packet transform in fault diagnosis of harmonic reducer under time-varying speed, the fast spectral kurtosis and the Hoyer diagram based on DTCWPT are used as control groups, such as Figure 12 shown.

[0213] Depend on Figure 12 (a) shows that the optimal frequency band in the fast spectral kurtosis graph belongs to level 0, and filtering operation cannot be performed; while in 12 (b), the required optimal frequency band (f c ,B ω ), but the frequency resolution is low and the Hoyer index value is too small, so the effect is poor; Figure 11 The desired optimal frequency band (f c ,B ω ), the rest of the frequency resolution is high, and the center frequency f corresponding to the maximum Hoyer index c The optimal bandwidth is 2830Hz. ω It is 196Hz, paving the way for generating a time-frequency spectrum with clearer time-frequency ridges.

[0214] 6.3 Experimental Analysis of Time-Frequency Ridge Extraction

[0215] 6.3.1SAH-SSET Combined with Fast Path Optimization to Extract Frequency

[0216] Use STFT and SSET to directly perform envelope analysis on the downsampled signal. Figure 13(a) Figure 13 (b) As shown. The band-stop filter is set based on SAH, and the passband cutoff frequency is [0f s / 2], f s is the sampling frequency of the signal, and the stopband cutoff frequency is [(f c -B ω / 2)(f c +B ω / 2)], after removing the resonance component, the spectrum of STFT and SSET envelope is as follows Figure 13 (c) and 13(d).

[0217] Figure 13 (a) Figure 13 (b) In the low-frequency region, the frequency information is locally affected by the cross term, and accurate time-frequency ridges cannot be extracted; Figure 13 In (c), we can preliminarily see the complete frequency conversion and its 2nd frequency curve, but the time-frequency resolution is very low and the energy is divergent; Figure 13 In (d), the frequency conversion and its double frequency curve can be clearly seen, and the time-frequency resolution is relatively high. Figure 13 (c) High, high energy Figure 13 (c) Concentration. By comparing STFT, SSET, SAH-STFT, and SAH-SSET, it is verified that the proposed method can effectively highlight the time-frequency ridges, improve the frequency ambiguity, and reduce the time-frequency ridge error. Figure 13 (d) The time-frequency graph extracts the time-frequency ridge as the rotation axis IF, and compares the sensor rotation frequency with the time-frequency ridge. Figure 14 shown.

[0218] Figure 14 The proposed SAH-SSET method combined with fast path optimization was verified to automatically extract accurate time-frequency ridges from the time-frequency spectrum. These ridges closely matched the sensor's measured speed, demonstrating that the extracted ridges roughly reflect the actual speed and clearly depict its changing trends. The main source of error is the inevitable frequency fluctuations that occur during actual time-varying speeds.

[0219] 6.3.2 Comparative Experiment

[0220] To verify the effectiveness of the proposed method for extracting time-frequency ridges, a comparative experiment was conducted using the signal measured from a harmonic reducer with an inner race fault as an example, using particle swarm optimization (PSO)-based VMD and improved complete ensemble empirical mode decomposition with adaptive noise (ICEEMDAN) as comparison methods. The inner race fault vibration signal was first processed using the comparison method. Time-frequency ridges were then extracted using SSET to generate a time-frequency plot combined with a fast path optimization algorithm. Finally, the accuracy of the time-frequency ridges was evaluated using quantitative metrics.

[0221] The particle swarm optimization algorithm is used to search for the optimal parameter combination penalty factor α and modal decomposition number K of VMD, with envelope entropy as the fitness value and the maximum number of iterations set to 10. The iterative optimization process is as follows: Figure 15 As shown in Figure 2, the iterative process takes about 30 minutes.

[0222] The optimal penalty factor α and modal decomposition number K are 2291 and 14 respectively. The decomposition results of the VMD for the inner ring fault signal of the harmonic reducer under time-varying speed are as follows: Figure 16 As shown in the figure, the correlation between each mode obtained by decomposition and the original signal is calculated, and the modes with a Pearson correlation coefficient greater than 0.6 are selected to form the reconstructed signal

[29] .

[0223] ICEEMDAN is developed based on Empirical Mode Decomposition (EMD). Setting the Gaussian white noise standard deviation to 0.1, the number of average decompositions of the signal to 10, and the maximum number of allowed screening iterations to 10, the decomposition results of the inner race fault signal of the harmonic reducer under time-varying speed by ICEEMDAN are as follows: Figure 17 As shown in the figure, the present invention also calculates the correlation between the IMF components obtained by ICEEMDAN decomposition and the original signal, and selects the modal components with a Pearson correlation coefficient greater than 0.6 to reconstruct the signal

[29] .

[0224] Figure 16 The Pearson correlation coefficient between IMF3 and the original signal in all modes is greater than 0.6, which is 0.6062, so IMF3 is the reconstructed signal. Figure 17 The Pearson correlation coefficients of IMF1, IMF2 and the original signal in all modes are greater than 0.6, which are 0.6718 and 0.7428 respectively. The reconstructed signal is obtained by summing IMF1 and IMF2. The reconstructed signals of PSO-VMD and ICEEMDAN are transformed by SSET, and the time-frequency diagram is made and the time-frequency ridge is extracted. The time-frequency diagram, time-frequency ridge and sensor frequency comparison of PSO-VMD-SSET are shown in the figure below. Figure 18As shown, the time-frequency diagram, time-frequency ridge and sensor frequency comparison diagram of ICEEMDAN-SSET are shown in Figure 19 shown.

[0225] from Figure 18 (a) It can be seen that in the time-frequency diagram, the frequency conversion and its multiples are locally blurred and discontinuous, and the time-frequency feature expression effect is poor. This is reflected in the jump phenomenon of the time-frequency ridge in 18 (b). This is because the local energy of the ridge representing the frequency conversion is too low during the time-frequency ridge extraction process. As a result, although there is a penalty factor to suppress frequency jumps during the extraction process, the fast path optimization method still selects high-energy points. Figure 19 (a) It can be seen that there are local fluctuations in the rotation frequency and its multiples in the time-frequency diagram, which is reflected in the local change trend of the time-frequency ridge in 19 (b) and the problem of large distortion error. The root mean square error (RMSE) and determination coefficient R are calculated using the sensor rotation frequency and time-frequency ridge. 2 , used to measure the error of the time-frequency ridge line. The specific results are shown in Table 3. RMS and determination coefficient R 2 The calculation of is shown in formulas (29) and (30).

[0226]

[0227] Among them, n represents the number of sample points used for comparison, y i Indicates the sensor rotation frequency, represents the time-frequency ridge, In this experiment, the number of sample points of the time-frequency ridge and the sensor rotation frequency is 20,000, so n is set to 20,000.

[0228] Table 3 Evaluation index results of different analysis methods

[0229]

[0230]

[0231] The RSME value of the proposed SAH-SSET combined with fast path optimization is closer to 0 than that of the other two analysis methods, indicating that the average deviation between the time-frequency ridge and the actual rotation frequency is small; and the determination coefficient R 2 The value is closer to 1, which shows that the time-frequency ridge can better reflect the changing trend of the actual rotation frequency, and verifies that the proposed method can extract the time-frequency ridge with smaller error.

[0232] 6.3.3 Order Analysis

[0233] SAH-SSET is used in combination with the fast path optimization method to extract the time-frequency ridge line and perform order analysis on the harmonic reducer inner ring fault envelope signal. First, the rotation angle of the harmonic reducer is calculated based on the time-frequency ridge line. A signal is collected every 1°, and angular domain resampling is performed. The envelope order spectrum is obtained by performing FFT on the angular domain envelope signal. Figure 20 As shown. r is the input shaft rotation frequency order, O i is the characteristic order of the inner race fault.

[0234] Figure 20 The input shaft rotation frequency order and its 2-fold and 3-fold orders as well as the inner ring fault characteristic order and its 2-fold order can be clearly observed, and it can be determined that the inner ring of the harmonic reducer is faulty. Similarly, the method proposed in this invention is used to analyze the flexible bearing outer ring fault signal, the flexible wheel tooth missing signal and the rigid wheel tooth missing signal. The results are as follows: Figure 21 As shown. r is the input shaft rotation frequency order, O o , O CMT and O FMT They are the orders of outer ring, rigid pulley and flexspline respectively.

[0235] Figure 21 In (a), the characteristic order of the outer ring fault and its multiple orders can be clearly observed, and it can be determined that the outer ring of the harmonic reducer has a fault. Figure 21 (b) Figure 21 (c) It can be seen that the fault characteristic order of the flex spline tooth missing fault is the same as that of the rigid wheel tooth missing fault, while the characteristic order generated by the flex spline self-rotation is too small to be reflected in the envelope order spectrum, which verifies the theoretical analysis of the fault characteristic orders of the flex spline and the rigid wheel in Section 4.2 of the present invention. Therefore, the amplitude of the fault characteristic order is needed as the basis for distinguishing the two. In terms of the working principle, during the operation of the harmonic reducer in this experiment, the flex spline is constantly moving and deforming, while the rigid wheel is fixed, so that the flex spline will produce a large amplitude impact at the fault point. It can be concluded that the characteristic order amplitude of the flex spline tooth missing fault should be larger than the characteristic order amplitude of the rigid wheel tooth missing fault. From the experimental results, the order amplitude of the flex spline tooth missing fault is larger, the fault characteristic order amplitude is above 0.45, and the 2-fold order amplitude is above 0.25; the order amplitude of the rigid wheel tooth missing fault is smaller, the fault characteristic order amplitude is around 0.25, and the 2-fold order amplitude is around 0.15, which is in line with the operating principle. Thus, it can be concluded that according to Figure 21 (b) Determine if the harmonic reducer flexible wheel is faulty. Figure 21 (c) Determine if the harmonic reducer pulley is faulty.

[0236] The difference between the experimentally obtained fault feature order and the theoretical fault feature order was divided by the theoretical fault feature order to calculate the error between orders in the envelope order spectrum. The results are shown in Table 4. It can be seen that the fault feature orders obtained by the proposed method are consistent with the theoretical fault feature orders. Except for the inner ring fault feature order, the error is less than 1%. The inner ring fault feature order error is larger than the other fault orders, but it is only 1.59%. This verifies that the proposed method can accurately extract the fault feature order of the harmonic reducer under time-varying speed and realize the fault diagnosis of the harmonic reducer under time-varying speed.

[0237] Table 4 Experimental error analysis

[0238]

[0239] In order to further verify the effectiveness of the method, the order analysis of the envelope signal of the inner ring fault of the wave reducer is carried out based on the time-frequency ridge obtained by PSO-VMD-SSET. The results are as follows: Figure 22 Based on the time-frequency ridge obtained by ICEEMDAN-SSET, the order analysis of the envelope vibration signal of the inner ring of the wave reducer is performed, and the results are shown as follows. Figure 23 shown.

[0240] Figure 22 The input shaft rotation frequency order and its 2 times and 3 times orders, as well as the inner race fault characteristic order can be observed, but the 2 times order of the inner race fault characteristic order is not prominent. Figure 23 The input shaft rotational frequency order and its 2nd and 3rd orders, as well as the inner race fault characteristic order and its 2nd order, can be observed. Table 2 shows that the theoretical fault characteristic order for the inner race fault is 12.60. By dividing the difference between the experimental and theoretical fault characteristic orders obtained by PSO-VMD-SSET by the theoretical fault characteristic order, the inner race fault characteristic order error is calculated to be 8.25%. By dividing the difference between the experimental and theoretical fault characteristic orders obtained by ICEEMDAN-SSET by the theoretical fault characteristic order, the inner race fault characteristic order error is calculated to be 5.16%. The error in the 2nd order of the inner race fault characteristic order is 3.02%. In comparison, the tachometer-free order analysis method based on SAH-SSET provides more accurate fault characteristic orders and smaller errors, enabling more precise determination of harmonic reducer fault types.

[0241] Although the present invention is disclosed as above, the scope of protection disclosed by the present invention is not limited thereto. Those skilled in the art of the present invention may make various changes and modifications without departing from the spirit and scope of the present invention, and these changes and modifications will fall within the scope of protection of the present invention.

[0242] References

[0243] [1] CHEN Zhuofan, ZHOU Kun, QIN Feifei, et al. Inverse Kinematics Solution of Robots Based on IQPSO Algorithm[J]. China Mechanical Engineering, 2024, 35(02): 293-304. doi:10.11999 / JEIT240203.

[0244] [2] ZHI Zhuo, LIU Liansheng, LIU Datong, etal. Fault detection of the harmonic reducer based on CNN-LSTM with a novel denoising algorithm [J]. IEEESensorsJournal, 2021, 22(3):2572-2581.doi:10.1109 / JSEN.2021.3137992.

[0245] [3]LIU Liansheng,ZHI Zhuo,YANG Yufei,etal.Harmonic reducer faultdetection with acoustic emission[J].IEEE Transactions on InstrumentationMeasurement,2023,72:1-12.doi:10.1109 / TIM.2023.3291747.

[0246] [4] KANG Shouqiang, YANG Jianxuan, WANG Yujing, et al. A Fast Classification Method of Rolling Bearing State Under Different Loads Based on Improved BroadModel Transfer Learning[J]. Journal of Electronics & Information Technology, 2023, 45(05): 1824-1832. doi: 10.11999 / JEIT220401.

[0247] [5]JIAYunzhao, LI Yuqing, XU Minqiang, et al. A fault diagnosis scheme for harmonic reducer under practical operating conditions [J]. Measurement, 2024, 227: 1-23. doi: 10.1016 / j.measurement.2024.114234.

[0248] [6]WANG Jing,WAN Zhihua,DONG Zhurong,etal.Research on performancetest system of space harmonic reducer in high vacuum and low temperatureenvironment[J].Machines,2020,9(1):1-11.doi:10.3390 / machines9010001.

[0249] [7] LI Yifan, ZHANG Xin, CHEN Zaigang, et al. Time-frequency ridgeestimation: An effective tool for gear and fault bearing diagnosis at time-varying speeds[J]. Mechanical Systems and Signal Processing, 2023, 189: 1-17. doi: 10.1016 / j.ymssp.2023.110108.

[0250] [8] ZHANG Xuhui, ZHANG Chao, FAN Hongwei, et al. Improved Fault Diagnosis of Rolling Bearing by Fast Kurtogram and Order Analysis[J]. Journal of Vibration, Measurement & Diagnosis, 2021, 41(06): 1090-1095+1235. doi: 10.16450 / j.cnki.issn.1004-6801.2021.06.007. ZHANG Xuhui, ZHANG Chao, FAN Hongwei, et al. Improved Fault Diagnosis of Rolling Bearing by Fast Kurtogram and Order Analysis[J]. Journal of Vibration, Measurement & Diagnosis, 2021, 41(06): 1090-1095+1235. doi: 10.16450 / j.cnki.issn.1004-6801.2021.06.007.

[0251] [9] ZHANG Chao and MAIMAITIREYIMU Abulizi. Fault Diagonosis of Variable Rotating Speed ​​Rolling Bearing Based on EOE_LMD and Order Tracking Analysis[J]. Journal of Vibration and Shock, 2024, 43(07): 308-316. doi: 10.13465 / j.cnki.jvs.2024.07.032.

[0252]

[10] YANGJianhua,YANG Chen,ZHUANG Xuzhu,etal.Unknown bearing faultdiagnosis under time-varying speed conditions and strong noise background.[J]Nonlinear Dyn.,2022,107(3):2177–2193.doi:10.1007 / s11071-021-07078-8.

[0253]

[11] ZHAO Ming and LINJing,Health assessment ofrotating machineryusing a rotary encoder.[J]IEEE Transactions on Industrial Electronics.2018,65(3):2548–2556.doi:10.1109 / TIE.2017.2739689.

[0254]

[12] XU Lei,DING Kang,HE Guolin,etal.A novel tacholess order trackingmethod for gearbox vibration signal based on extremums search ofgearmeshharmonic[J]Mechanical Systems and Signal Processing,2023,189:1-20.doi:10.1016 / j.ymssp.2022.110070.

[0255]

[13] ZHANG Yan, HE Shubei, WANG Ping, et al. Tacholess quantitative characterization of rolling bearing fault feature under varying conditions[J]. Chinese Journal of Scientific Instrument, 2021, 41(08): 104-114. doi: 10.19650 / j.cnki.cjsi.J2107704. ZHANG Yan, HE Shubei, WANG Ping, et al. Tacholess quantitative characterization of rolling bearing fault feature under varying conditions[J]. Chinese Journal of Scientific Instrument, 2021, 41(08): 104-114. doi: 10.19650 / j.cnki.cjsi.J2107704.

[0256]

[14] MENG Lingxia, Xu Xiaoli, and ZUO Yunbo. Fault feature extraction of logarithmic time-frequency ridge order spectrum of planetary gearbox under time-varying conditions[J]. Journal of Vibration and Shock, 2020, 39(07): 163-169+179. doi: 10.13465 / j.cnki, jvs. 2020.07.023.

[0257]

[15] ANIL Kumar, ZHOU Yuqing, and XIANG Jiawei. Optimization of VMD using kernelbased mutual information for the extraction of weak features to detect bearing defects [J] Measurement, 2021, 168: 1-13.doi: 10.1016 / j.measurement.2020.108402.

[0258]

[16] DUAN Rongkai,LIAO Yuhe,YANG Lei,et al.Adaptive tacholess ordertracking method based on generalized linear chirplet transform and itsapplication for bearing fault diagnosis [J].ISA Transactions,2022,127:324-341.doi:10.1016 / j.isatra.2021.08.039.

[0259]

[17] LI Xinyan,ZHAO Huimin,YU Ling,et al.Feature Extraction UsingParameterized Multisynchrosqueezing Transform.[J]IEEE Sensors Journal,2022,22(14):14263-14272.doi:10.1109 / JSEN.2022.3179165.

[0260]

[18] TANG Lei,SHANG Xuqiang,HUANG Tianli,et al.An improved localmaximum synchrosqueezing transform with adaptive window width forinstantaneous frrequency identification of time-varying structures[J].Engineering Structures,2023,292:1-15.doi:10.1016 / j.engstruct.2023.116543.

[0261]

[19] BAO Wenjie,TU Xiaotong,LI Fucai,et al.GeneralizedSynchrosqueezing Transform:Algorithm and Applications[J].IEEE Transactions onInstrumentation and Measurement,2023,72:5293-5302.doi:10.1109 / TIE.2020.2984983.

[0262]

[20] WAN Shuting, WANG Yanjie, ZHANG Xiong, et al. Tacho-less order tracking method of wind turbine gearbox under time-Varying conditions[J]. Journal of Vibration Engineering, 2023, 36(01): 266-279. doi: 16385 / j.cnki.issn.1004-4523.2023.01.028. WAN Shuting, WANG Yanjie, ZHANG Xiong, et al. Tacho-less order tracking method of wind turbine gearbox under time-Varying conditions[J]. Journal of Vibration Engineering, 2023, 36(01): 266-279. doi: 16385 / j.cnki.issn.1004-4523.2023.01.028.

[0263]

[21] JIANG Zuhua, ZHANG Kun, ZHANG Xiangfeng, et al. A tacholess ordertracking method based on spectral amplitude modulation for variable speedbearing fault diagnosis [J]. IEEE Transactions on Instrumentation and Measurement, 2023, 72: 1-8.doi: 10.1109 / TIM.2023.3280512.

[0264]

[22] WANG Lei, LIU Zhiwen, CAO Hongrui, et al. Subband averaging kurtogram with dual-tree complex wavelet packet transform for rotating machinery faultdiagnosis [J]. Mechanical Systems and Signal Processing, 2020, 142: 1-21.doi: 10.1016 / j.ymssp.2020.106755.

[0265]

[23] WANG Dong, Spectral L2 / L1 norm: A new perspective for spectralkurtosis for characterizing non-stationary signals [J]. Mechanical Systems andSignal Processing, 2018, 104: 290-293.doi: 10.1016 / j.ymssp.2017.11.013.

[0266]

[24] TU Xiaotong, HE Zhoujie, HU Yue, et al. The Second OrderSynchroextracting Transform with Application to Bearing Fault Diagnosis underVariable Speed ​​Condition[C]. Asia Pacific Conference of the Prognostics and Health Management, Beijing, China, 2019: 1-5.

[0267]

[25] D.Iatsenko, PVEMcclintock, and A.Stefanovska, Extractionofinstantaneous frequencies from ridges in time-frequency representations ofsignals[J].Signal Process.2016, 125: 290-303.doi: 10.1016 / j.sigpro.2016.01.024.

[0268]

[26] Zhong Bowen. Dynamic simulation and fault characteristic analysis of harmonic drive system[D]. [Master's thesis]. Harbin Institute of Technology, 2023.

[0269] ZHONG Bowen.Research on Harmonic Transmissionsvstem DvnamicsSimulation and Faultcharacteristic Analysis[D].Master dissertation.HarbinInstitute of Technology, 2023.

[0270]

[27] GUO Yingying. Study on Kinematic Characteristics and Fault Diagnosis of Flexible Thin-wall Elliptical Bearing[D]. Ph.D. dissertation. South China University of Technology, 2021.

[0271]

[28] JIA Yunzhao, XU Minqiang, CHENG Yao, et al. Calculation Method for The Harmonic Drive Fault Characteristic Frequency[J]. Journal of Vibration Engineering, 2025, 44(02): 279-291. doi: 10.13465 / j.cnki.jvs.2025.02.028.

[0272]

[29] ZHAO Li, WANG Xiaogang, WANG Ning, et al. Fully Generalized Space Modulation Visible Light Communications System Based on Pearson Correlation Coefricient Selection[J]. Acta Optica Sinica, 2024, 44(04): 116-124. doi: 10.3788 / AOS231379.

Claims

1. A harmonic reducer fault diagnosis method combining SAH with SSET under time-varying speed, characterized in that: The following steps are involved: S100, collecting time domain vibration signals of typical faults of harmonic reducers under time-varying speed; S200, using SAH to determine the optimal filtering parameters of the time domain vibration signal of the harmonic reducer collected in step S100, including the center frequency and the optimal bandwidth, and perform band-stop filtering; S300, performing Hilbert transform on the filtered signal and using the envelope time-frequency diagram obtained by SSET, then performing fast path optimization and extracting the time-frequency ridge as the rotation axis IF; S400, performing Hilbert transform on the original harmonic reducer vibration signal and obtaining its envelope, then performing angular domain resampling on the envelope signal according to the time-frequency ridge extracted in step S300 to obtain an angular domain envelope signal, and performing fast Fourier transform to obtain an envelope order spectrum; S500. Determine the fault type by comparing the theoretical fault characteristic order of the harmonic reducer with the envelope order spectrum.

2. The method for diagnosing harmonic reducer faults using SAH combined with SSET under time-varying speed according to claim 1, characterized in that: In step S200, it includes: S210, construct the subband spectral kurtosis using the wavelet coefficients, as shown in formula (1), the spectral kurtosis is the square of the spectral L2 / L1 norm, The Hoyer index is the normalized form of the L2 / L1 norm, as shown in formula (2), Where a is the L2 / L1 norm, N is the signal length, n represents the discrete time variable, K k,r is the spectral kurtosis, H k,r is the Hoyer index, d k,r is the sub-band wavelet coefficient rearranged in descending order of center frequency, k is the number of decomposition layers, and each layer divides the center frequency into 2 k intervals, r is one of the intervals; S210, the step of constructing SAH includes: S211, using a sliding window to evenly divide the collected harmonic reducer vibration signal x(t) into M sub-signals {x m (t)|m=1,2,...,M}, the window length of the sliding window is N x / M, the step size is equal to the window length, N x is the length of signal x(t); S212, using sub-band rearrangement DWCWPT to decompose each sub-signal into a set of sub-band wavelet coefficients S213, calculating the Hoyer index of the subband; S214. Calculate the average value of the Hoyer index of the corresponding sub-band to obtain the sub-band average Hoyer index As shown in formula (3): S215, draw SAH, select the best frequency band corresponding to the maximum sub-band average Hoyer index (f c ,B ω ), center frequency f c and the optimal bandwidth B ω , which is the optimal filtering center frequency and bandwidth.

3. The method for diagnosing harmonic reducer faults using SAH combined with SSET under time-varying speed according to claim 2, characterized in that: In step S200, it also includes: When using SAH to obtain the resonance frequency band (f c ,B ω ), the vibration signal of the harmonic reducer is band-stop filtered to obtain the vibration signal x1(t). The STFT expression of x1(t) based on the window function g(t) is: Where τ is the time-shift variable, * is the conjugate operation, g(t) is the window function, g(τ-t) is the sliding window function, i is the imaginary unit, ω is the normalized angular frequency of the input signal, and t represents the continuous time variable; The Gaussian window function with standard deviation σ is expressed as:

4. The method for diagnosing harmonic reducer faults using SAH combined with SSET under time-varying speed according to claim 3, characterized in that: In step S300, it includes: Definition Multi-component AM / FM signal is defined as: Among them, K is the number of modes, A k (t),φ k (t) represents the instantaneous amplitude and instantaneous frequency of each mode respectively; The ideal time-frequency representation of x1(t) is shown in formula (7): Where δ(·) is the delta function, which represents the impulse response; φ′ k (t) is the ideal instantaneous frequency; Assume that the single-component signal is locally approximated by a Gaussian modulated linear chirp signal at every moment: in, represents the instantaneous amplitude of x1(t); c0 is the initial phase, c1 is the initial frequency, c2 is the modulation frequency, and t0 is the envelope center time; According to formula (8), the derivative of x1(t) with respect to t is expressed as: Among them, the imaginary part is the signal IF; is the composite parameter obtained from derivative analysis, A(t) is the amplitude; Derived The partial derivative with respect to t is: Then deduce The partial derivative with respect to t is: in, represents the STFTs obtained using g″(t),tg′(t),g′(t) as window functions, which can be obtained by using formulas (10) and (11): φ′(t,ω) is expressed as formula (13): in, To take the imaginary part; The TFR of SSET is defined as formula (14):

5. The method for diagnosing harmonic reducer faults using SAH combined with SSET under time-varying speed according to claim 4, characterized in that: In step S300, the method further includes extracting a time-frequency ridge based on fast path optimization, using a support function to optimally extract the time-frequency ridge from the TFR with minimal frequency jump; Assume that the TFR of signal x1(t) is SSET X1(t,τ), where τ is the time variable and f is the frequency variable, and the time-frequency ridge f is extracted from X1(τ,f) using fast path optimization p (τ), and the upper boundary of the time-frequency ridge f up (τ) and the lower bound f down (τ); The time τ in TFR n The number of peaks is recorded as N p (τ n ), the frequency corresponding to the m-th peak is recorded as V m (τ n ), the TFR amplitude of the mth peak is recorded as Q m (τ n ); τ n Peak value at time v m (τ n ) is defined by formula (15): If X1(τ,f) has a time span [τ1, τ2, τ3, …, τ n ], path optimization is described as: Among them, m c (τ n ) is to determine the peak to be extracted as τ n The ridge at the position, F[·] is the selected optimization support function, {m1,m2,m3,....,m N } represents the sequence of peak numbers along the time span. Fast path optimization only relies on a limited number of preceding points, and the support function is: Among them, f d (t n-1 ) is the candidate ridge point at τ n-1 The frequency, f d is the frequency of a series of candidate ridge points in history [τ1, τ2, τ3, …, τ n-1 ], △f d f d The derivative of the series, m[·] is the median of the series, IQR is the interquartile range, λ1 and λ2 are penalty factors, and the weight functions ω1(·) and ω2(·) are used to suppress the atypical changes of the ridge frequency value and the derivative, respectively.

6. The method for diagnosing harmonic reducer faults using SAH combined with SSET under time-varying speed according to claim 1, characterized in that: In step S400, it includes: The fault characteristic order of the inner and outer rings of the flexible bearing is derived through the calculation formula of the bearing fault characteristic frequency: Among them, O i , O o are the orders of the inner and outer rings of the flexible bearing of the harmonic reducer, f i ,f o are the fault characteristic frequencies of the inner and outer rings of the flexible bearing of the harmonic reducer, f r is the rotation frequency, Z is the number of rolling elements, D is the pitch diameter, d is the ball pitch diameter, and α is the contact angle; It is inferred that the flexspline and the rigid wheel have a common fault characteristic frequency of 2 times the rotation frequency, and the flexspline also has a fault characteristic frequency caused by self-rotation (2 / Z F ) times the rotation frequency, Z F is the number of flexspline teeth. According to the relationship between the fault characteristic order and the fault characteristic frequency, formulas (20) and (21) are obtained: O CMT =O FMT1 =f CMT,FMT1 / f r =2f r / f r =2 (20) O FMT2 =f FMT2 / f r =(2 / Z F )*f r / f r =2 / Z F (21) Among them, O CMT and O FMT are the orders of the rigid wheel and flexible wheel of the harmonic reducer, O FMT1 is one of the fault characteristic orders of the flexspline, f CMT,FMT1 is the fault characteristic frequency of the harmonic reducer rigid wheel and flexible wheel, O FMT2 is the flexspline rotation order of the harmonic reducer, f FMT2 It is the fault characteristic frequency caused by the self-rotation of the flexspline of the harmonic reducer.

7. The method for diagnosing harmonic reducer faults using SAH combined with SSET under time-varying speed according to claim 1, characterized in that: In step S500, the flexspline and rigid-wheel faults are distinguished from each other in terms of amplitude. When the flexspline fault occurs, the amplitude is relatively obvious due to the change in the fault point location, and abundant sidebands appear around it. When the rigid-wheel fault occurs, the fault location remains unchanged, the amplitude is relatively small, and the sidebands are relatively few. Accordingly, in the order spectrum, the characteristic order amplitude of the flexspline fault is larger than that of the rigid-wheel fault.

8. A harmonic reducer fault diagnosis system combining SAH with SSET under time-varying speed, characterized by: The system has a program module corresponding to the steps of any one of claims 1 to 7, and executes the steps of the harmonic reducer fault diagnosis method combining SAH with SSET under time-varying speed during operation.

9. A computer-readable storage medium, characterized in that: The computer-readable storage medium stores a computer program, and the computer program is configured to implement the steps of the harmonic reducer fault diagnosis method combining SAH with SSET under time-varying speed according to any one of claims 1 to 7 when called by a processor.

Images