Multi-impact-source coupled vibration signal time-frequency domain adaptive decomposition method and system

Abstract

The invention discloses a multi-impact-source coupling vibration signal time-frequency domain adaptive decomposition method and system, and belongs to the field of mechanical signal processing. Firstly, the time-frequency energy gradient of a multi-impact-source coupling vibration signal is calculated; secondly, estimating a noise amplitude mean value of each sub-sequence in an energy gradient plane in a signal time direction, and calculating a weight sequence based on impact morphological characteristics of the sub-sequences so as to obtain a weighted impact recognition threshold value; further, detecting the number of impact sources in an energy gradient plane in the time direction as a target decomposition number, and calculating a time-frequency filtering boundary of sub-impact based on an impact filtering boundary self-adaptive positioning strategy; and finally, based on the time-frequency filtering boundary of the sub-impact, obtaining decomposed sub-signals and time-frequency centers thereof through time-frequency inverse transformation. The method is efficient and accurate, fully pays attention to the energy change rule of the impact signal, realizes time-frequency domain adaptive decomposition of the sub-impact component, and is beneficial to real-time state monitoring and rapid fault identification of reciprocating machinery.

Description

Technical Field

[0001] The present invention relates to the field of mechanical signal processing, and in particular to a method and system for adaptively decomposing a multi-impact source coupled vibration signal in time and frequency domains. Background Art

[0002] The vibration mechanism of reciprocating machinery is closely related to the multi-source excitation generated during the working process. The explosion pressure in the cylinder, the closing of the intake and exhaust valves, the collision between the piston and the cylinder liner, and the occurrence of faults will all generate impact sources with different amplitudes and frequencies. Various impact forces are transmitted to the cylinder head through the mechanical structure, and coupled to form the cylinder head vibration response of the reciprocating machinery. Therefore, the vibration signal of the reciprocating machinery presents the time-frequency coupling characteristics of multiple impact sources in the time-frequency domain. How to decouple the coupled vibration signal of multiple impact sources is a key research difficulty in the field of reciprocating machinery vibration signal processing.

[0003] In recent years, some scholars have begun to conduct signal decomposition research on the signal characteristics of multiple impact source coupling in reciprocating machinery. The time domain decomposition method can separate the impact sources generated by different working cycles in the vibration signal of reciprocating machinery, but due to the limitations of the analysis domain, multi-frequency impacts in the same time domain range cannot be separated. Time-frequency analysis provides a more comprehensive analysis perspective for the time-frequency coupled signals of multiple impact sources in reciprocating machinery. However, the existing time-frequency decomposition methods still do not adequately describe the impact characteristics, and it is difficult to accurately locate the time-frequency filtering boundaries of the impact, resulting in insufficient impact positioning accuracy and computational efficiency that needs to be improved. Therefore, in view of the limitations of the existing methods, the present invention discloses a time-frequency domain adaptive decomposition method and system for multi-impact source coupled vibration signals, which can accurately and efficiently decouple the time-frequency coupled signals of multiple impact sources, greatly shorten the processing time of single working cycle signals, and has engineering application promotion value. Summary of the invention

[0004] Existing reciprocating mechanical signal decomposition methods pay little attention to the energy variation law of impact and it is difficult to remove the noise in the impact area, so there is room for improvement in the description of impact characteristics and the location of impact boundaries. The present invention proposes a time-frequency domain adaptive decomposition method and system for multi-impact source coupled vibration signals. Based on the time-frequency characteristics of the signal, it accurately describes the time-frequency energy variation law of the impact, can adaptively define the weighted impact recognition threshold, accurately locate the time-frequency filtering boundary of the sub-impact, and realize the time-frequency domain adaptive decomposition of multi-impact source coupled signals.

[0005] The present invention provides a method and system for adaptively decomposing a multi-impact source coupled vibration signal in time and frequency domains, which is characterized by comprising the following steps:

[0006] S1: Collect and obtain the mechanical multi-impact source coupled vibration acceleration signal x(t) to be processed, where t is the sampling time;

[0007] S2: Calculate the time-frequency energy gradient of the multi-impact source coupled vibration signal x(t);

[0008] S3: Estimate the noise amplitude mean of each subsequence in the energy gradient plane of the signal time direction, calculate the weight sequence based on the shock morphological characteristics of the subsequence to obtain the weighted shock recognition threshold;

[0009] S4: Detect the number of impact sources in the energy gradient plane in the time direction as the target decomposition number, and then calculate the time-frequency filtering boundary of the sub-impact based on the impact filtering boundary adaptive positioning strategy;

[0010] S5: Based on the time-frequency filtering boundary of the sub-shock, the decomposed sub-signal and its time-frequency center are obtained through inverse time-frequency transform;

[0011] In one embodiment of the present invention, in step S2, the calculating of the time-frequency energy gradient of the multi-impact source coupled vibration signal x(t) comprises the following steps:

[0012] S201: Loading multiple impact source coupled vibration signals x(t);

[0013] S202: Perform a short-time Fourier transform on x(t) to obtain a time-frequency representation G(t,f) of x(t), where t represents the sampling time and f represents the frequency. The short-time Fourier transform is implemented by the package function "stft" in MATLAB.

[0014] S203: Calculate the time-frequency energy gradient of x(t), which is the energy gradient in the time direction and frequency direction energy gradient composition:

[0015] E(t,f)=|G(t,f)| 2

[0016]

[0017] Where E(t,f) represents the time-frequency energy, represents the energy gradient in the time direction, express The sampling time in the plane is t i , frequency is f j The gradient value at time, E(t i+1 ,f j ) indicates that the sampling time in the E(t,f) plane is t i+1 , frequency is f j The energy value at time, E(t i ,f j ) indicates that the sampling time in the E(t,f) plane is t i , frequency is f j The energy value at time ; similarly, represents the energy gradient in the frequency direction, express The sampling time in the plane is t i , frequency is f j The gradient value at time, E(t i ,f j+1 ) indicates that the sampling time in the E(t,f) plane is t i , frequency is f j+1 energy value at the time; i and j represent the index of sampling time and frequency respectively, and |·| represents the absolute value operator;

[0018] In one embodiment of the present invention, in step S3, the noise amplitude mean of each subsequence in the energy gradient plane of the estimated signal time direction is calculated based on the impact morphological characteristics of the subsequence to obtain a weighted impact recognition threshold, including the following steps:

[0019] S301: Energy gradient in time direction The time-frequency noise of each subsequence in the plane is estimated, and the time-frequency noise estimation method is as follows:

[0020] (1) Assuming that the energy gradient in the time direction The size of is P × Q, then The plane contains P subsequences, each subsequence has Q points, record where g p express The pth subsequence of the plane, g p =[g p (1),g p (2),…,g p (q),…,g p (Q)], q = 1, 2, ..., Q, where g p (q) represents the subsequence g p The qth value of , [·] T Represents the transpose of a matrix;

[0021] (2) Define the subsequence g p The noise component is g Noise , g Noise The mean amplitude of is T, and we have:

[0022] ∫ 1 Q g Noise (q) 2 dq=∫ 1 Q T 2 dq

[0023] Among them, g Noise Represents the subsequence gp The noise component, g Noise (q) represents g Noise The qth value of Noise Therefore, in order to estimate the noise amplitude mean T, it is necessary to introduce the denoising sequence g Denoised (T), and use T as the variable to find the optimal solution. Denoised The mathematical formula of (T) is:

[0024] g Denoised (T) = ReLU(|g p |-T)

[0025]

[0026] Among them, g Denoised (T) represents the denoised sequence, T is g Noise The amplitude mean, ReLU(·) represents the activation function, x is the variable of the function ReLU(·), and |·| represents the absolute value operator;

[0027] Define the residual energy function E(T):

[0028] E(T)=-∫ 1 Q g p (q) 2 dq+∫ 1 Q g Denoised (T) 2 dq+∫ 1 Q T 2 dq

[0029] Function E(T) represents the denoised sequence energy g Denoised (T) 2 The square of the noise amplitude mean T 2 If g p The maximum absolute value of (q) is T max , then there must be a zero point T∈[0,T max ], so that E(T) = 0. Use the bisection method to iteratively approximate the equation E(T) = 0 step by step, and you can get the zero point, and this zero point is the estimated noise amplitude mean value T;

[0030] (3) Based on (1)-(2) in step S301, we calculate The time-frequency noise amplitude mean of each subsequence in the plane is obtained to obtain the noise mean sequence T Noise =[T 1 ,T 2 ,…,T p ,…,TP ],p=1,2,…,P,P represents The total number of subsequences of the plane, where T p represents the noise mean sequence T Noise The pth value of ;

[0031] S302: Calculate a weight sequence based on the impact morphological features of the subsequence to obtain a weighted impact recognition threshold. The adaptive delineation method of the weighted impact recognition threshold is as follows:

[0032] (1) Based on the noise mean sequence T obtained in step S301 Noise ,calculate The shock identification threshold of each subsequence in the plane. p The shock recognition threshold T I (p) is calculated as:

[0033] T I (p) = γ·T Noise (p)

[0034] Among them, T I (p) represents the subsequence g p The shock recognition threshold, γ represents the threshold adjustment coefficient, T Noise (p) indicates T Noise The pth value of g. p The sequence follows a normal distribution, and the value of γ is 2.0≤γ≤4.0;

[0035] Calculate separately The shock recognition threshold of each subsequence in the plane is obtained, and the shock recognition threshold sequence T is obtained. I =[T I (1),T I (2),…,T I (p),…,T I (P)], p = 1, 2, ..., P, P represents The total number of subsequences of the plane, where T I (p) represents the shock recognition threshold sequence T I The pth value of , that is, the subsequence g p The shock recognition threshold;

[0036] (2) Based on the shock recognition threshold sequence T I , calculate the weighted shock recognition threshold T WI , the calculation method is as follows:

[0037] First, calculate Shock distortion index of each subsequence in the plane:

[0038]

[0039] Among them, V(g p ) represents the subsequence g p Impact distortion index, g p (i) represents the subsequence g p The i-th value of , i = 1, 2, ..., Q-3, Q is the subsequence g p The number of sample points, max(·) represents the maximum value function, min(·) represents the minimum value function, PE(g p ) represents the subsequence g p The permutation entropy is calculated by the MATLAB package function "getPermEn", where "·" represents the multiplication operator;

[0040] Calculate separately The shock distortion index of each subsequence in the plane is obtained, and the shock distortion index sequence V = [V(1), V(2), ..., V(p), ..., V(P)], p = 1, 2, ..., P, P represents The total number of subsequences of the plane, where V(p) represents the p-th value of the shock distortion index sequence V;

[0041] Secondly, the weight of each subsequence is calculated based on the shock distortion index sequence V:

[0042]

[0043] Among them, w(i) represents the weight corresponding to the i-th subsequence, T I (i) indicates T I The i-th sample value of , V(i) represents the i-th sample of V. The weights corresponding to each subsequence in the plane constitute the sequence w = [w(1), w(2), ..., w(p), ..., w(P)], p = 1, 2, ..., P, where P represents The total number of subsequences in the plane. Probability normalization of the weight sequence w:

[0044]

[0045] Among them, w N (i) represents the normalized weight of the i-th subsequence;

[0046] Finally, calculate Weighted shock recognition threshold for a plane:

[0047] T WI =w N T I T

[0048] Among them, T WI express Weighted shock recognition threshold for the plane, w N express Normalized weight of the plane, T I express The shock recognition threshold of the plane, the superscript T indicates the transpose of the matrix;

[0049] In one embodiment of the present invention, in step S4, the number of impact sources in the energy gradient plane in the detection time direction is used as the target decomposition number, and then the time-frequency filtering boundary of the sub-impact is calculated based on the impact filter boundary adaptive positioning strategy, including the following steps:

[0050] S401: Adaptively obtain the target decomposition number. The specific method is as follows: Based on the weighted impact recognition threshold T obtained in step S3 WI , search In-plane greater than threshold T WI The maximum value of each connected area, the time-frequency position of the maximum value is the peak position of the impact source. If the number of impact source peaks is K, then the target decomposition number is K;

[0051] S402: adaptively calculating the time-frequency filtering boundary, the specific method is as follows:

[0052] First, the time domain filtering boundary is calculated. Based on the peak position of each impact source obtained in step S401, the energy gradient in the time direction is The time-direction energy gradient sequences corresponding to the K impact sources are stripped off on the plane respectively, and the first zero-crossing point immediately before the peak of each impact source is taken as the starting boundary of the time-domain filtering, and the second zero-crossing point after the peak is taken as the ending boundary of the time-domain filtering, thereby obtaining the time-domain filtering boundary of each impact source;

[0053] Secondly, the frequency domain filter boundary is calculated. Similarly, based on the peak position of each impact source, the energy gradient in the frequency direction is calculated. The frequency direction energy gradient sequences corresponding to the K impact sources are stripped off on the plane respectively, and the first zero-crossing point immediately before the peak of each impact source is taken as the starting boundary of the frequency domain filtering, and the second zero-crossing point after the peak is taken as the ending boundary of the frequency domain filtering, so as to obtain the frequency domain filtering boundary of each impact source;

[0054] Finally, the time domain filtering boundary and frequency domain filtering boundary of each impact source calculated together constitute the time-frequency filtering boundary of each impact source;

[0055] In one embodiment of the present invention, in step S5, the sub-shock-based time-frequency filtering boundary obtains the decomposed sub-signal and its time-frequency center through inverse time-frequency transformation, including the following steps:

[0056] S501: Based on the time-frequency filtering boundary of each sub-impact source, the decomposed sub-signal corresponding to each impact source is obtained by inverse short-time Fourier transform to form a decomposed sub-signal set. The inverse short-time Fourier transform is implemented by the package function "istft" in MATLAB;

[0057] S502: Obtain the time-frequency center of the decomposed sub-signals to form a time-frequency center set;

[0058] S503: The decomposed sub-signal set and the time-frequency center set are used as the time-frequency domain adaptive decomposition result of the multi-impact source coupled vibration signal and output.

[0059] Another object of the present invention is to provide a time-frequency domain adaptive decomposition system for multi-impact source coupled vibration signals, characterized in that it includes:

[0060] A signal acquisition module, used to acquire mechanical multi-impact source coupled vibration acceleration signals to be processed;

[0061] A time-frequency transformation module is used to perform short-time Fourier transformation on the multi-impact source coupled vibration acceleration signal to be processed to obtain a time-frequency representation of the signal;

[0062] The signal time-frequency domain decomposition module is used to obtain the time-frequency filtering boundary of each sub-impact source through the time-frequency domain adaptive decomposition method according to the time-frequency representation of the signal. The time-frequency domain adaptive decomposition method includes: calculating the time-frequency energy gradient of the multi-impact source coupled vibration signal; estimating the noise amplitude mean of each sub-sequence in the energy gradient plane of the signal time direction, calculating the weight sequence based on the impact morphological characteristics of the sub-sequence to obtain the weighted impact recognition threshold; detecting the number of impact sources in the energy gradient plane of the time direction as the target decomposition number, and then calculating the time-frequency filtering boundary of the sub-impact source based on the impact filter boundary adaptive positioning strategy;

[0063] The time-frequency inverse transform module is used to perform short-time Fourier inverse transform processing on each sub-shock source according to the time-frequency filtering boundary of the sub-shock source, obtain the decomposed sub-signals corresponding to each shock source, and form a decomposed sub-signal set; and obtain the time-frequency center of the decomposed sub-signal to form a time-frequency center set; the decomposed sub-signal set and the time-frequency center set are used as the time-frequency domain adaptive decomposition results of the multi-shock source coupled vibration signal and are output.

[0064] The present invention aims at the problems that existing time-frequency decomposition methods are difficult to accurately locate the time-frequency filtering boundaries of impacts, resulting in insufficient impact positioning accuracy and the need to improve the calculation efficiency. A time-frequency domain adaptive decomposition method and system for multi-impulse source coupled vibration signals are proposed. Based on the time-frequency characteristics of multi-impulse source coupled signals, the time-frequency energy gradient of the signal is calculated; the mean noise amplitude of each subsequence in the time-direction energy gradient plane of the signal is estimated, and the weight sequence is calculated based on the impact morphology characteristics of the subsequence to obtain the weighted impact recognition threshold; further, the number of impact sources in the time-direction energy gradient plane is detected as the target decomposition number, and then the time-frequency filtering boundaries of sub-impacts are calculated based on the impact filtering boundary adaptive positioning strategy; finally, based on the time-frequency filtering boundaries of sub-impacts, the decomposed sub-signals and their time-frequency centers are obtained through inverse time-frequency transformation, realizing the accurate and efficient decoupling of multi-impulse source time-frequency coupled signals.

[0065] The beneficial technical effects of the present invention are as follows:

[0066] 1. The time-frequency energy gradient defined by the present invention accurately describes the time-frequency domain characteristics of impact components, highlighting the energy change law of impacts, and thus has excellent noise robustness.

[0067] 2. The weighted impact recognition threshold adaptive calculation method proposed by the present invention delimits the threshold based on the time-frequency noise energy, thereby adaptively obtaining the impact decomposition number, which helps to accurately identify impact components.

[0068] 3. The impact filtering boundary adaptive positioning strategy proposed by the present invention utilizes the random energy distribution characteristics of noise in the time-frequency domain, thereby efficiently and accurately obtaining the time-frequency filtering boundaries of sub-impacts.

[0069] 4. The present invention has both high decomposition accuracy and high decomposition efficiency, greatly shortening the processing time of single working cycle vibration signals, and has broad application prospects and popularization significance in the field of reciprocating machinery condition monitoring and fault diagnosis. Description of the Drawings

[0070] Figure 1 is a flowchart of a time-frequency domain adaptive decomposition method for multi-impulse source coupled vibration signals provided by an embodiment of the present invention;

[0071] Figure 2 is a time-domain diagram of a diesel engine cylinder head vibration signal provided by an embodiment of the present invention;

[0072] Figure 3 is a frequency-domain diagram of a diesel engine cylinder head vibration signal provided by an embodiment of the present invention;

[0073] Figure 4 is a short-time Fourier transform time-frequency diagram of a diesel engine cylinder head vibration signal provided by an embodiment of the present invention;

[0074] Figure 5 is a time-frequency diagram of a sub-signal obtained by adaptive decomposition in an embodiment of the present invention;

[0075] Figure 6 is a time domain diagram of a sub-signal obtained by adaptive decomposition in an embodiment of the present invention;

[0076] Figure 7 It is a schematic diagram of a time-frequency domain adaptive decomposition system for multi-impact source coupled vibration signals provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0077] In order to make the objectives, technical solutions and advantages of the present invention more clear, the embodiments of the present invention will be further described in detail below with reference to the accompanying drawings.

[0078] Figure 1 This is a flow chart of a method for adaptively decomposing a multi-impact source coupled vibration signal in the time-frequency domain provided by an embodiment of the present invention. Figure 1 A method for adaptively decomposing vibration signals in time and frequency domains from multiple impact sources includes the following steps:

[0079] S1: Collect and obtain the mechanical multi-impact source coupled vibration acceleration signal x(t) to be processed, where t is the sampling time;

[0080] S2: Calculate the time-frequency energy gradient of the multi-impact source coupled vibration signal x(t);

[0081] S3: Estimate the noise amplitude mean of each subsequence in the energy gradient plane of the signal time direction, calculate the weight sequence based on the shock morphological characteristics of the subsequence to obtain the weighted shock recognition threshold;

[0082] S4: Detect the number of impact sources in the energy gradient plane in the time direction as the target decomposition number, and then calculate the time-frequency filtering boundary of the sub-impact based on the impact filtering boundary adaptive positioning strategy;

[0083] S5: Based on the time-frequency filtering boundary of the sub-shock, the decomposed sub-signal and its time-frequency center are obtained through inverse time-frequency transform.

[0084] In one embodiment of the present invention, in step S2, the calculating of the time-frequency energy gradient of the multi-impact source coupled vibration signal x(t) comprises the following steps:

[0085] S201: Loading multiple impact source coupled vibration signals x(t);

[0086] S202: Perform a short-time Fourier transform on x(t) to obtain a time-frequency representation G(t,f) of x(t), where t represents the sampling time and f represents the frequency. The short-time Fourier transform is implemented by the package function "stft" in MATLAB.

[0087] S203: Calculate the time-frequency energy gradient of x(t), which is the energy gradient in the time direction and frequency direction energy gradient composition:

[0088] E(t,f)=|G(t,f)| 2

[0089]

[0090] Where E(t,f) represents the time-frequency energy, represents the energy gradient in the time direction, express The sampling time in the plane is t i , frequency is f j The gradient value at time, E(t i+1 ,f j ) indicates that the sampling time in the E(t,f) plane is t i+1 , frequency is f j The energy value at time, E(t i ,f j ) indicates that the sampling time in the E(t,f) plane is t i , frequency is f j The energy value at time ; similarly, represents the energy gradient in the frequency direction, express The sampling time in the plane is t i , frequency is f j The gradient value at time, E(t i ,f j+1 ) indicates that the sampling time in the E(t,f) plane is t i , frequency is f j+1 energy value at the time; i and j represent the index of sampling time and frequency respectively, and |·| represents the absolute value operator;

[0091] In one embodiment of the present invention, in step S3, the noise amplitude mean of each subsequence in the energy gradient plane of the estimated signal time direction is calculated based on the impact morphological characteristics of the subsequence to obtain a weighted impact recognition threshold, including the following steps:

[0092] S301: Energy gradient in time direction The time-frequency noise of each subsequence in the plane is estimated, and the time-frequency noise estimation method is as follows:

[0093] (1) Assuming that the energy gradient in the time direction The size of is P × Q, then The plane contains P subsequences, each subsequence has Q points, record where gp express The pth subsequence of the plane, g p =[g p (1),g p (2),…,g p (q),…,g p (Q)], q = 1, 2, ..., Q, where g p (q) represents the subsequence g p The qth value of , [·] T Represents the transpose of a matrix;

[0094] (2) Define the subsequence g p The noise component is g Noise , g Noise The mean amplitude of is T, and we have:

[0095] ∫ 1 Q g Noise (q) 2 dq=∫ 1 Q T 2 dq

[0096] Among them, g Noise Represents the subsequence g p The noise component, g Noise (q) represents g Noise The qth value of Noise Therefore, in order to estimate the noise amplitude mean T, it is necessary to introduce the denoising sequence g Denoised (T), and use T as the variable to find the optimal solution. Denoised The mathematical formula of (T) is:

[0097] g Denoised (T) = ReLU(|g p |-T)

[0098]

[0099] Among them, g Denoised (T) represents the denoised sequence, T is g Noise The amplitude mean, ReLU(·) represents the activation function, x is the variable of the function ReLU(·), and |·| represents the absolute value operator;

[0100] Define the residual energy function E(T):

[0101] E(T)=-∫ 1 Q g p (q) 2 dq+∫1 Q g Denoised (T) 2 dq+∫ 1 Q T 2 dq

[0102] Function E(T) represents the denoised sequence energy g Denoised (T) 2 The square of the noise amplitude mean T 2 If g p The maximum absolute value of (q) is T max , then there must be a zero point T∈[0,T max ], so that E(T) = 0. Use the bisection method to iteratively approximate the equation E(T) = 0 step by step, and you can get the zero point, and this zero point is the estimated noise amplitude mean value T;

[0103] (3) Based on (1)-(2) in step S301, we calculate The time-frequency noise amplitude mean of each subsequence in the plane is obtained to obtain the noise mean sequence T Noise =[T 1 ,T 2 ,…,T p ,…,T P ],p=1,2,…,P,P represents The total number of subsequences of the plane, where T p represents the noise mean sequence T Noise The pth value of ;

[0104] S302: Calculate a weight sequence based on the impact morphological features of the subsequence to obtain a weighted impact recognition threshold. The adaptive delineation method of the weighted impact recognition threshold is as follows:

[0105] (1) Based on the noise mean sequence T obtained in step S301 Noise ,calculate The shock identification threshold of each subsequence in the plane. p The shock recognition threshold T I (p) is calculated as:

[0106] T I (p) = γ·T Noise (p)

[0107] Among them, T I (p) represents the subsequence g p The shock recognition threshold, γ represents the threshold adjustment coefficient, T Noise (p) indicates T Noise The gth value of . pThe sequence follows a normal distribution, and the value of γ is 2.0≤γ≤4.0;

[0108] It should be noted that the value of the threshold adjustment coefficient γ can be determined according to the actual decomposition task requirements. If it is necessary to extract the weak impact component in the signal, γ takes a small value; if the weak impact component in the signal is not concerned, γ takes a large value. Preferably, γ is set to 3.0;

[0109] Calculate separately The shock recognition threshold of each subsequence in the plane is obtained, and the shock recognition threshold sequence T is obtained. I =[T I (1),T I (2),…,T I (p),…,T I (P)], p = 1, 2, ..., P, P represents The total number of subsequences of the plane, where T I (p) represents the shock recognition threshold sequence T I The pth value of , that is, the subsequence g p The shock recognition threshold;

[0110] (2) Based on the shock recognition threshold sequence T I , calculate the weighted shock recognition threshold T WI , the calculation method is as follows:

[0111] First, calculate Shock distortion index of each subsequence in the plane:

[0112]

[0113] Among them, V(g p ) represents the subsequence g p Impact distortion index, g p (i) represents the subsequence g p The i-th value of , i = 1, 2, ..., Q-3, Q is the subsequence g p The number of sample points, max(·) represents the maximum value function, min(·) represents the minimum value function, PE(g p ) represents the subsequence g p The permutation entropy is calculated by the MATLAB package function "getPermEn", where "·" represents the multiplication operator;

[0114] Calculate separately The shock distortion index of each subsequence in the plane is obtained, and the shock distortion index sequence V = [V(1), V(2), ..., V(p), ..., V(P)], p = 1, 2, ..., P, P represents The total number of subsequences of the plane, where V(p) represents the p-th value of the shock distortion index sequence V;

[0115] Secondly, the weight of each subsequence is calculated based on the shock distortion index sequence V:

[0116]

[0117] Among them, w(i) represents the weight corresponding to the i-th subsequence, T I (i) indicates T I The i-th sample value of , V(i) represents the i-th sample of V. The weights corresponding to each subsequence in the plane constitute the sequence w = [w(1), w(2), ..., w(p), ..., w(P)], p = 1, 2, ..., P, where P represents The total number of subsequences in the plane. Probability normalization of the weight sequence w:

[0118]

[0119] Among them, w N (i) represents the normalized weight of the i-th subsequence;

[0120] Finally, calculate Weighted shock recognition threshold for a plane:

[0121] T WI =w N T I T

[0122] Among them, T WI express Weighted shock recognition threshold for the plane, w N express Normalized weight of the plane, T I express The shock recognition threshold of the plane, the superscript T indicates the transpose of the matrix;

[0123] In one embodiment of the present invention, in step S4, the number of impact sources in the energy gradient plane in the detection time direction is used as the target decomposition number, and then the time-frequency filtering boundary of the sub-impact is calculated based on the impact filter boundary adaptive positioning strategy, including the following steps:

[0124] S401: Adaptively obtain the target decomposition number. The specific method is as follows: Based on the weighted impact recognition threshold T obtained in step S3 WI , search In-plane greater than threshold T WIThe maximum value of each connected area, the time-frequency position of the maximum value is the peak position of the impact source. If the number of impact source peaks is K, then the target decomposition number is K;

[0125] S402: adaptively calculating the time-frequency filtering boundary, the specific method is as follows:

[0126] First, the time domain filtering boundary is calculated. Based on the peak position of each impact source obtained in step S401, the energy gradient in the time direction is The time-direction energy gradient sequences corresponding to the K impact sources are stripped off on the plane respectively, and the first zero-crossing point immediately before the peak of each impact source is taken as the starting boundary of the time-domain filtering, and the second zero-crossing point after the peak is taken as the ending boundary of the time-domain filtering, thereby obtaining the time-domain filtering boundary of each impact source;

[0127] Secondly, the frequency domain filter boundary is calculated. Similarly, based on the peak position of each impact source, the energy gradient in the frequency direction is calculated. The frequency direction energy gradient sequences corresponding to the K impact sources are stripped off on the plane respectively, and the first zero-crossing point immediately before the peak of each impact source is taken as the starting boundary of the frequency domain filtering, and the second zero-crossing point after the peak is taken as the ending boundary of the frequency domain filtering, so as to obtain the frequency domain filtering boundary of each impact source;

[0128] Finally, the time domain filtering boundary and frequency domain filtering boundary of each impact source calculated together constitute the time-frequency filtering boundary of each impact source;

[0129] In one embodiment of the present invention, in step S5, the sub-shock-based time-frequency filtering boundary obtains the decomposed sub-signal and its time-frequency center through inverse time-frequency transformation, including the following steps:

[0130] S501: Based on the time-frequency filtering boundary of each sub-impact source, the decomposed sub-signal corresponding to each impact source is obtained by inverse short-time Fourier transform to form a decomposed sub-signal set. The inverse short-time Fourier transform is implemented by the package function "istft" in MATLAB;

[0131] S502: Obtain the time-frequency center of the decomposed sub-signals to form a time-frequency center set;

[0132] S503: The decomposed sub-signal set and the time-frequency center set are used as the time-frequency domain adaptive decomposition result of the multi-impact source coupled vibration signal and output.

[0133] According to a second aspect of the present invention, a time-frequency domain adaptive decomposition system for multi-impact source coupled vibration signals is characterized by comprising:

[0134] A signal acquisition module, used to acquire mechanical multi-impact source coupled vibration acceleration signals to be processed;

[0135] A time-frequency transformation module is used to perform short-time Fourier transformation on the multi-impact source coupled vibration acceleration signal to be processed to obtain a time-frequency representation of the signal;

[0136] The signal time-frequency domain decomposition module is used to obtain the time-frequency filtering boundary of each sub-impact source through the time-frequency domain adaptive decomposition method according to the time-frequency representation of the signal. The time-frequency domain adaptive decomposition method includes: calculating the time-frequency energy gradient of the multi-impact source coupled vibration signal; estimating the noise amplitude mean of each sub-sequence in the energy gradient plane of the signal time direction, calculating the weight sequence based on the impact morphological characteristics of the sub-sequence to obtain the weighted impact recognition threshold; detecting the number of impact sources in the energy gradient plane of the time direction as the target decomposition number, and then calculating the time-frequency filtering boundary of the sub-impact source based on the impact filter boundary adaptive positioning strategy;

[0137] The time-frequency inverse transform module is used to perform short-time Fourier inverse transform processing on each sub-shock source according to the time-frequency filtering boundary of the sub-shock source, obtain the decomposed sub-signals corresponding to each shock source, and form a decomposed sub-signal set; and obtain the time-frequency center of the decomposed sub-signal to form a time-frequency center set; the decomposed sub-signal set and the time-frequency center set are used as the time-frequency domain adaptive decomposition results of the multi-shock source coupled vibration signal and are output.

[0138] The specific embodiments are as follows:

[0139] Taking a certain model of diesel engine as an example, in order to fully demonstrate the complete content of the present invention, the vibration acceleration signal of the diesel engine cylinder head is selected as an example. The diesel engine cylinder head vibration signal is collected under the condition of a speed of 1000rpm, and the sampling frequency is 51.2kHz. The single working cycle signal x(t) is intercepted as the vibration signal to be processed, and its time domain is as follows Figure 2 shown. Figure 3 and Figure 4 The frequency domain diagram and short-time Fourier transform time-frequency diagram of the diesel engine single working cycle cylinder head vibration signal x(t) to be processed are shown respectively. Figure 2-Figure 4 It can be seen that the signal to be processed exhibits the characteristics of time-frequency coupling of multiple impact sources.

[0140] The method of the present invention is used to analyze and process the collected cylinder head vibration signal. According to step S2, first, load the multi-impact source coupled vibration signal x(t); then, perform short-time Fourier transform on x(t) to obtain the time-frequency representation of x(t), such as Figure 4 As shown; finally, the time-frequency energy gradient of x(t) is calculated. According to step S3, first, the time-frequency noise of each subsequence in the energy gradient plane in the time direction is estimated; then, the weight sequence is calculated based on the impact morphological characteristics of the subsequence, so as to obtain the weighted impact recognition threshold T WI=0.0072. According to step S4, first, the number of impact sources in the energy gradient plane in the time direction is detected as the target decomposition number, and the adaptively obtained target decomposition number K=10; then, the time-frequency filtering boundaries of 10 sub-impacts are calculated based on the impact filter boundary adaptive positioning strategy. According to step S5, first, based on the time-frequency filtering boundaries of each sub-impact source, the decomposed sub-signals corresponding to each impact source are obtained by short-time Fourier inverse transform to form a decomposed sub-signal set. The time-frequency diagram of the decomposed sub-signals is shown as follows: Figure 5 As shown, the time domain diagram of the decomposed sub-signal is as follows Figure 6 As shown. Figure 5 and Figure 6 It can be known that the time-frequency domain adaptive decomposition method of multi-impact source coupled vibration signals provided by the present invention can effectively decouple multi-impact source time-frequency coupled signals. Then, the time-frequency center of the decomposed sub-signal is obtained to form a time-frequency center set; finally, the decomposed sub-signal set and the time-frequency center set are used as the time-frequency domain adaptive decomposition results of the multi-impact source coupled vibration signal and output. The time-frequency center set of 10 sub-signals is shown in Table 1. The unit "s" of the time center in Table 1 represents "seconds", and the unit "Hz" of the frequency center represents "Hertz". The total running time of this embodiment is 3.52 seconds, which shows that the present invention can quickly process single working cycle signals and has engineering application promotion value.

[0141] Table 1 Time-frequency center set of decomposed sub-signals

[0142]

[0143] Based on the above embodiments, it can be seen that the present invention accurately describes the time-frequency domain characteristics of the signal impact component based on the time-frequency energy gradient, and can adaptively define the impact recognition threshold through time-frequency noise estimation, and adaptively obtain the target decomposition number. In addition, the impact filter boundary adaptive positioning strategy provided by the present invention can accurately calculate the time-frequency filter boundary of the sub-impact, and realize the time-frequency domain adaptive decomposition of the multi-impact source coupling signal.

[0144] Compared with the above-mentioned method for adaptively decomposing a vibration signal of multiple impact sources in time and frequency domain, another embodiment of the present invention provides an adaptive decomposition system of a vibration signal of multiple impact sources in time and frequency domain, such as Figure 7As shown in the figure, it mainly includes: a signal acquisition module, a time-frequency transformation module, a signal time-frequency decomposition module, and a time-frequency inverse transformation module. Specifically, the signal acquisition module is used to acquire the coupled vibration acceleration signals of mechanical multi-impulse sources to be processed; the time-frequency transformation module is used to perform short-time Fourier transform processing on the coupled vibration acceleration signals of the multi-impulse sources to be processed to obtain the time-frequency representation of the signals; the signal time-frequency domain decomposition module is used to obtain the time-frequency filtering boundaries of each sub-impulse source according to the time-frequency representation of the signals through a time-frequency domain adaptive decomposition method. The time-frequency domain adaptive decomposition method includes: calculating the time-frequency energy gradient of the coupled vibration signals of multi-impulse sources; estimating the mean noise amplitude of each subsequence in the plane of the signal time-direction energy gradient, and calculating a weight sequence based on the impact morphological characteristics of the subsequence to obtain a weighted impact recognition threshold; detecting the number of impulse sources in the plane of the time-direction energy gradient as the target decomposition number, and then calculating the time-frequency filtering boundaries of the sub-impulse sources based on the impact filtering boundary adaptive positioning strategy; the time-frequency inverse transformation module is used to perform short-time Fourier inverse transform processing on each sub-impulse source respectively according to the time-frequency filtering boundaries of the sub-impulse sources to obtain the decomposed sub-signals corresponding to each impulse source, forming a set of decomposed sub-signals; and obtaining the time-frequency centers of the decomposed sub-signals, forming a set of time-frequency centers; and outputting the set of decomposed sub-signals and the set of time-frequency centers as the time-frequency domain adaptive decomposition result of the coupled vibration signals of multi-impulse sources.

[0145] It should be noted that the above embodiment of the time-frequency domain adaptive decomposition system for coupled vibration signals of multi-impulse sources corresponds to the embodiment of the time-frequency domain adaptive decomposition method for coupled vibration signals of multi-impulse sources and has the same technical effects. For specific descriptions, please refer to the embodiment of the time-frequency domain adaptive decomposition method for coupled vibration signals of multi-impulse sources. In addition, the embodiment of the time-frequency domain adaptive decomposition system for coupled vibration signals of multi-impulse sources is obtained based on the embodiment of the time-frequency domain adaptive decomposition method for coupled vibration signals of multi-impulse sources, so it is briefly described here and will not be elaborated further.

[0146] The above is only the preferred embodiment of the present application and the description of the applied technical principles. Those skilled in the art should understand that the scope involved in the present application is not limited to the technical solution formed by the specific combination of the above technical features, but should also cover other technical solutions formed by any combination of the above technical features or their equivalent features without departing from the inventive concept. For example, the technical solutions formed by the mutual replacement of the above features and the technical features (but not limited to) disclosed in the present application with similar functions.

[0147] Except for the technical features described in the specification, the remaining technical features are known to those skilled in the art. To highlight the innovative features of the present invention, the remaining technical features will not be elaborated further here.

Claims

1. A method for adaptively decomposing vibration signals in time and frequency domains from multiple impact sources, characterized in that: The following steps are involved: S1: Collect and obtain the mechanical multi-impact source coupled vibration acceleration signal x(t) to be processed, where t is the sampling time; S2: Calculate the time-frequency energy gradient of the multi-impact source coupled vibration signal x(t); S3: Estimate the noise amplitude mean of each subsequence in the energy gradient plane of the signal time direction, calculate the weight sequence based on the shock morphological characteristics of the subsequence to obtain the weighted shock recognition threshold; S4: Detect the number of impact sources in the energy gradient plane in the time direction as the target decomposition number, and then calculate the time-frequency filtering boundary of the sub-impact based on the impact filtering boundary adaptive positioning strategy; S5: Based on the time-frequency filtering boundary of the sub-shock, the decomposed sub-signal and its time-frequency center are obtained through inverse time-frequency transform; Wherein, in step S2, the time-frequency energy gradient of the multi-impact source coupled vibration signal x(t) is calculated, comprising the following steps: S201: Loading multiple impact source coupled vibration signals x(t); S202: Perform short-time Fourier transform on x(t) to obtain the time-frequency representation G(t, f) of x(t), where t represents the sampling time and f represents the frequency. The short-time Fourier transform is implemented by the package function "stft" in MATLAB. S203: Calculate the time-frequency energy gradient of x(t), which is the energy gradient in the time direction and frequency direction energy gradient composition: E(t,f)=|G(t,f)| 2 Where E(t, f) represents the time-frequency energy, represents the energy gradient in the time direction, express The sampling time in the plane is t i , frequency is f j The gradient value at time, E(t i+1 , f j ) indicates that the sampling time in the E(t, f) plane is t i+1 , frequency is f j The energy value at time, E(t i , f j ) indicates that the sampling time in the E(t, f) plane is t i , frequency is f j The energy value at time; similarly, represents the energy gradient in the frequency direction, express The sampling time in the plane is t i , frequency is f j The gradient value at time, E(t i , f j+1 ) indicates that the sampling time in the E(t, f) plane is t i , frequency is f j+1 energy value at the time; i and j represent the index of sampling time and frequency respectively, and |·| represents the absolute value operator; Among them, in step S3, the noise amplitude mean of each subsequence in the energy gradient plane of the estimated signal time direction is calculated based on the impact morphological characteristics of the subsequence to obtain the weighted impact recognition threshold, including the following steps: S301: Energy gradient in time direction The time-frequency noise of each subsequence in the plane is estimated, and the time-frequency noise estimation method is as follows: (1) Assuming that the energy gradient in the time direction The size of is P × Q, then The plane contains P subsequences, each subsequence has Q points, record p=1,2,...,P,where g p express The pth subsequence of the plane, g p =[g p (1), g p (2), ..., g p (q), ..., g p (Q)], q = 1, 2, ..., Q, where g p (q) represents the subsequence g p The qth value of , [·] T Represents the transpose of a matrix; (2) Define the subsequence g p The noise component is g Noise , g Noise The mean amplitude of is T, and we have: ∫1 Q g Noise (q) 2 dq=∫1 Q T 2 dq Among them, g Noise Represents the subsequence g p The noise component, g Noise (q) represents g Noise The qth value of Noise The amplitude mean; therefore, in order to estimate the amplitude mean T of the noise, it is necessary to introduce the denoising sequence g Denoised (T), take T as the variable to find the optimal solution; g Denoised The mathematical formula of (T) is: g Denoised (T)0ReLU(|g p |-T) Among them, g Denoised (T) represents the denoised sequence, T is g Noise The amplitude mean, ReLU(·) represents the activation function, x is the variable of the function ReLU(·), and |·| represents the absolute value operator; Define the residual energy function E(T): E(T)=-∫1 Q g p (q) 2 dq+∫1 Q g Denoised (T) 2 dq+∫1 Q T 2 dq Function E(T) represents the denoised sequence energy g Denoised (T) 2 The square of the noise amplitude mean T 2 The difference between p The maximum absolute value of (q) is T max , then there must be a zero point T∈[0, T max ], so that E(T) = 0; the equation E(T) = 0 is solved by iterative approximation using the bisection method, and the zero point can be obtained, and the zero point is the estimated noise amplitude mean value T; (3) Based on (1)-(2) in step S301, we calculate The time-frequency noise amplitude mean of each subsequence in the plane is obtained to obtain the noise mean sequence T Noise =[T1, T2, ..., T p , ..., T P ], p = 1, 2, ..., P, P represents The total number of subsequences of the plane, where T p represents the noise mean sequence T Noise The pth value of ; S302: Calculate a weight sequence based on the impact morphological features of the subsequence to obtain a weighted impact recognition threshold. The adaptive delimitation method of the weighted impact recognition threshold is as follows: (1) Based on the noise mean sequence T obtained in step S301 Noise ,calculate The impact recognition threshold of each subsequence in the plane; subsequence g p The shock recognition threshold T I (p) is calculated as: T I (p)=γ·T Noise (p) Among them, T I (p) represents the subsequence g p The shock recognition threshold, γ represents the threshold adjustment coefficient, T Noise (p) indicates T Noise The pth value of g p The sequence follows a normal distribution, and the value of γ is 2.0≤γ≤4.0; Calculate separately The shock recognition threshold of each subsequence in the plane is obtained, and the shock recognition threshold sequence T is obtained. I =[T I (1), T I (2), ..., T I (p), ..., T I (P)], p = 1, 2, ..., P, P represents The total number of subsequences of the plane, where T I (p) represents the shock recognition threshold sequence T I The pth value of p The shock recognition threshold; (2) Based on the shock recognition threshold sequence T I , calculate the weighted shock recognition threshold T WI , the calculation method is as follows: First, calculate Shock distortion index of each subsequence in the plane: Among them, V(g p ) represents the subsequence g p Impact distortion index, g p (i) represents the subsequence g p The i-th value of , i = 1, 2, ..., Q-3, Q is the subsequence g p The number of sample points, max(·) represents the maximum value function, min(·) represents the minimum value function, PE(g p ) represents the subsequence g p The permutation entropy is calculated by the MATLAB package function "getPermEn", where "·" represents the multiplication operator; Calculate separately The shock distortion index sequence V = [V(1), V(2), ..., V(p), ..., V(P)], p = 1, 2, ..., P, P represents The total number of subsequences of the plane, where V(p) represents the p-th value of the shock distortion index sequence V; Secondly, the weight of each subsequence is calculated based on the shock distortion index sequence V: Among them, w(i) represents the weight corresponding to the i-th subsequence, T I (i) indicates T I The i-th sample value of , V(i) represents the i-th sample of V; The weights corresponding to each subsequence in the plane constitute the sequence w = [w(1), w(2), ..., w(p), ..., w(P)], p = 1, 2, ..., P, P represents The total number of subsequences of the plane; the probability normalization of the weight sequence w: Among them, w N (i) represents the normalized weight of the i-th subsequence; Finally, calculate Weighted shock recognition threshold for a plane: T WI =w N T I T Among them, T WI express Weighted impact recognition threshold for the plane, w N express Normalized weight of the plane, T I express The shock recognition threshold of the plane, the superscript T indicates the transpose of the matrix; Among them, in step S4, the number of impact sources in the energy gradient plane in the detection time direction is used as the target decomposition number, and then the time-frequency filtering boundary of the sub-impact is calculated based on the impact filter boundary adaptive positioning strategy, including the following steps: S401: Adaptively obtain the target decomposition number. The specific method is as follows: Based on the weighted impact recognition threshold T obtained in step S3 WI , search In-plane greater than threshold T WI The maximum value of each connected area, the time-frequency position of the maximum value is the peak position of the impact source. If the number of impact source peaks is K, then the target decomposition number is K; S402: adaptively calculating the time-frequency filtering boundary, the specific method is as follows: First, the time domain filtering boundary is calculated; based on the peak position of each impact source obtained in step S401, the energy gradient in the time direction is The time-direction energy gradient sequences corresponding to the K impact sources are stripped off on the plane respectively, and the first zero-crossing point immediately before the peak of each impact source is taken as the starting boundary of the time-domain filtering, and the second zero-crossing point after the peak is taken as the ending boundary of the time-domain filtering, thereby obtaining the time-domain filtering boundary of each impact source; Secondly, the frequency domain filtering boundary is calculated; similarly, based on the peak position of each impact source, the energy gradient in the frequency direction is calculated. The frequency direction energy gradient sequences corresponding to the K impact sources are stripped off on the plane respectively, and the first zero-crossing point immediately before the peak of each impact source is taken as the starting boundary of the frequency domain filtering, and the second zero-crossing point after the peak is taken as the ending boundary of the frequency domain filtering, so as to obtain the frequency domain filtering boundary of each impact source; Finally, the time domain filtering boundary and frequency domain filtering boundary of each impact source calculated together constitute the time-frequency filtering boundary of each impact source; Wherein, in step S5, the sub-shock-based time-frequency filtering boundary obtains the decomposed sub-signal and its time-frequency center through inverse time-frequency transformation, including the following steps: S501: Based on the time-frequency filtering boundary of each sub-impact source, the decomposed sub-signal corresponding to each impact source is obtained by inverse short-time Fourier transform to form a decomposed sub-signal set. The inverse short-time Fourier transform is implemented by the package function "istft" in MATLAB; S502: Obtain the time-frequency center of the decomposed sub-signals to form a time-frequency center set; S503: The decomposed sub-signal set and the time-frequency center set are used as the time-frequency domain adaptive decomposition result of the multi-impact source coupled vibration signal and output.

2. A system for implementing the method of claim 1, characterized in that: include: A signal acquisition module, used to acquire mechanical multi-impact source coupled vibration acceleration signals to be processed; A time-frequency transformation module is used to perform short-time Fourier transformation processing on the multi-impact source coupled vibration acceleration signal to be processed to obtain a time-frequency representation of the signal; A signal time-frequency domain decomposition module, used to obtain the time-frequency filtering boundary of each sub-impact source through a time-frequency domain adaptive decomposition method according to the time-frequency representation of the signal; The time-frequency domain adaptive decomposition method comprises: calculating the time-frequency energy gradient of the multi-impact source coupled vibration signal; Estimate the noise amplitude mean of each subsequence in the energy gradient plane of the signal time direction, calculate the weight sequence based on the shock morphological characteristics of the subsequence to obtain the weighted shock recognition threshold; detect the number of shock sources in the energy gradient plane of the time direction as the target decomposition number, and then calculate the time-frequency filtering boundary of the sub-shock source based on the shock filter boundary adaptive positioning strategy; The time-frequency inverse transform module is used to perform short-time Fourier inverse transform processing on each sub-shock source according to the time-frequency filtering boundary of the sub-shock source, obtain the decomposed sub-signals corresponding to each shock source, and form a decomposed sub-signal set; and obtain the time-frequency center of the decomposed sub-signal to form a time-frequency center set; the decomposed sub-signal set and the time-frequency center set are used as the time-frequency domain adaptive decomposition results of the multi-shock source coupled vibration signal and are output.

Images