Rotating unbalance fault detection method and system based on non-contact rotor micro-displacement trajectory

Abstract

The application discloses a rotating unbalance fault detection method and system based on non-contact rotor micro-displacement trajectory, and the method is as follows: S1, collecting multi-direction surface vibration data caused by the operation of a rotating device, i.e. an original data frame; S2, based on the original data frame, analyzing vibration target position information by using a vibration positioning algorithm; S3, based on the original data frame and the vibration target position information, extracting a target vibration signal by using a vibration phase extraction algorithm; S4, based on the target vibration signal, correcting and extracting a target phase signal by using a phase correction algorithm; S5, based on the target phase signal, extracting a vibration feature by using a vibration feature extraction algorithm; and S6, based on the extracted vibration feature, diagnosing a rotating unbalance fault by using a feature analysis algorithm. The application proposes to track a rotor trajectory by using two synchronous millimeter wave radars, calculate the eccentricity and flatness of the rotor trajectory, and combine the amplitude and rotating speed to diagnose the rotating unbalance fault, so that the precision is high and the cost is low.

Description

Technical Field

[0001] This invention belongs to the field of millimeter-wave vibration detection technology, and particularly relates to a millimeter-wave-based micro-displacement detection technology for vibration surfaces, specifically a method and system for detecting rotational imbalance faults based on non-contact rotor micro-displacement trajectories. Background Technology

[0002] Rotating equipment is a type of mechanical device capable of rotational motion, typically composed of rotating components, a drive mechanism, a support structure, and a control system. In industry, rotating equipment plays a wide and vital role in numerous fields, including manufacturing, processing, energy, transportation, chemicals, and pharmaceuticals. The rotor, as the core rotating component, converts electrical energy into mechanical kinetic energy. During operation, the rotor bears significant mechanical and thermal stresses, making it prone to failure. Rotational imbalance is a very common rotor fault, referring to the equipment's imbalance due to an imbalance of centrifugal forces during rotation. This fault can cause abnormal vibration, noise, and accelerated wear, and in severe cases, may lead to equipment damage or even accidents. Therefore, timely detection and repair of rotor imbalance faults are crucial for rotating equipment to ensure its normal operation and safety.

[0003] While various sensor devices, such as accelerometers and laser vibrometers, can be used to detect rotational imbalance faults, these devices are currently expensive and require demanding installation conditions. Millimeter-wave vibration measurement methods can conveniently and accurately extract the time-domain information of vibration displacement from multiple targets simultaneously, enabling large-scale, multi-scale, high-precision, and efficient multi-point synchronous deformation and vibration measurements. This provides a new non-contact method and approach for measuring the operating status of rotating equipment. Rotational imbalance faults can be identified by detecting abnormal frequencies, amplitudes, and phases caused by the fault using frequency and time-domain analysis methods. Based on this, this invention innovatively uses two synchronous millimeter-wave radars to track the rotor trajectory, calculate the eccentricity and flatness of the rotor trajectory, and combine this with amplitude and rotational speed to diagnose rotational imbalance faults. Summary of the Invention

[0004] In view of the above situation, the present invention proposes a method and system for detecting rotational imbalance faults based on non-contact rotor micro-displacement trajectory.

[0005] To achieve the above objectives, the present invention adopts the following technical solution:

[0006] A method for detecting rotational imbalance faults based on non-contact rotor micro-displacement trajectories includes the following steps:

[0007] S1. Radar data acquisition: Acquire multi-directional surface vibration data caused by the operation of rotating equipment, i.e., raw data frames.

[0008] S2. Vibration localization: Based on the original data frame, the vibration localization algorithm is used to analyze the position information of the vibration target.

[0009] S3. Vibration signal extraction: Based on the original data frame and the vibration target position information, the vibration phase extraction algorithm is used to extract the target vibration signal.

[0010] S4. Phase Correction: Based on the target vibration signal, a phase correction algorithm is used to correct and extract the target phase signal.

[0011] S5. Vibration Feature Extraction: Based on the target phase signal, vibration features are extracted using a vibration feature extraction algorithm.

[0012] S6. Feature Analysis: Based on the extracted vibration features, feature analysis algorithms are used to diagnose rotational imbalance faults.

[0013] Furthermore, the method involves two synchronized millimeter-wave radars and a radar synchronization system. The two synchronized millimeter-wave radars are placed radially along the rotor of a rotating device, with the radars at a 90-degree angle to each other relative to the rotor. Controlled by the radar synchronization system, the two millimeter-wave radars operate on different frequency bands and acquire data synchronously. The operating frequency bands of the two millimeter-wave radars are 77–79 GHz and 79–81 GHz, respectively.

[0014] Furthermore, the original data frame is a four-dimensional matrix of size [F,T,N,M], where F represents the number of frames in the target data frame (frame dimension), T represents the number of receiving antennas of the millimeter-wave radar (antenna dimension), N represents the number of chirp signals transmitted by the millimeter-wave radar (chirp dimension), and M represents the number of sampling points for each chirp echo (sampling dimension).

[0015] Furthermore, the vibration localization algorithm in step S2 is specifically implemented through the following sub-steps:

[0016] S2.1 For the original data frame, calculate the target distance matrix. Specifically, first extract a single original data frame (of size [T, N, M]) from the original data frame, and then perform a Fast Fourier Transform on the sampling dimension to obtain a target distance matrix of size [T, N, M].

[0017] S2.2. Based on the target range matrix, calculate the optimal antenna weighting coefficients and range-angle power spectrum. Specifically, first calculate the radar antenna's angle scan vector. Assuming the millimeter-wave radar is fixed, angle scanning is achieved by weighting along the antenna dimension. The formula for calculating the radar antenna's angle scan vector is shown below:

[0018] a(θ)=[a(θ1),a(θ2),...a(θc (1)

[0019]

[0020] Where T is the number of radar receiving antennas, the spacing between the antenna arrays is d, λ is the signal wavelength, and θ = {θ1, θ2, ..., θ...} C} represents the angular scanning range.

[0021] Then, the angular distance power spectrum and the optimal antenna weighting coefficient are calculated, and the calculation formulas are shown below:

[0022]

[0023]

[0024] Where p is the angular range power spectrum, w is the optimal antenna weighting coefficient matrix, and x is the target range matrix. H The conjugate transpose, * T Indicates the transpose operation, * -1 This indicates the inverse operation.

[0025] S2.3. Calculate vibration location information based on the angular distance power spectrum and the optimal antenna weighting coefficient matrix.

[0026] First, static background removal technology is used to remove interference from static backgrounds in the environment. Specifically, interference with insignificant power (static background) targets is achieved by subtracting the average power from the angular distance power spectrum. The calculation formula is shown below:

[0027]

[0028] Where p′ represents the calculation result, and p is the angular distance power spectrum. This indicates an average calculation.

[0029] Then, calculate the maximum index (C′,R′) of p′. (C′,R′) is used as the vibration orientation information, denoted as (θ,r).

[0030] Finally, the optimal antenna weighting coefficients for the vibrating target are extracted from the optimal antenna weighting coefficient matrix based on the vibration azimuth information. The extraction formula is shown below:

[0031] α=w(θ,r) (6)

[0032] Here, w is a three-dimensional matrix w[C,R,T]. For a given index (θ,r), the corresponding one-dimensional data w[C=θ,R=r] can be extracted from w, and the length of this data is T.

[0033] Furthermore, the phase extraction algorithm in step S3 is specifically implemented through the following sub-steps:

[0034] S3.1. Reassemble the original data frames to obtain the original data groups. To ensure the continuity of the target vibration signal, the original data frames of size [F,T,N,M] are reassembled, and the F-frame data is reassembled into G groups of data frames.

[0035] First, identify the frame indices contained in each data group. Specifically, for the frame sequence [1,2,3,4,5,6,7,...,F], perform sliding window sampling with a window size of width and a step size of step to divide the frame sequence into G groups. Taking the frame sequence [1,2,3,4,5,6,7], width=3, step=2 as an example, the result after division is [1,2,3], [3,4,5], [5,6,7].

[0036] Then, each set of data is merged. Specifically, multiple frames within a set of data are stacked along the chirp dimension. Taking a set of data [3,T,N,M] as an example, it represents three data frames of size [T,N,M]. After stacking along the chirp dimension, the result is [T,3*N,M]. Since each set has partially overlapping data, the extracted target vibration signal is continuous.

[0037] S3.2 For the original data set, calculate the target distance matrix. Specifically, for the original data set of size [G,T,N,M], perform FFT (Fast Fourier Transform) processing on the sampling dimension to obtain the target distance matrix x of size [G,T,N,M].

[0038] S3.3. Based on the target distance matrix and the vibration target position information, extract the target vibration signal. The extraction formula is as follows:

[0039] z=x(T,N,r)*α (7)

[0040] Where z is the extracted target vibration signal, r is the distance index of the vibration target location, α is the optimal antenna weighting coefficient, and x(T,N,r) is the target vibration signal at distance r extracted based on the distance index r. The vibration signal at distance r is weighted and summed using the optimal antenna weighting coefficient α to extract the vibration signal at angle θ.

[0041] Furthermore, the phase correction algorithm in step S4 is specifically implemented through the following sub-steps:

[0042] S4.1. For a given sub-vibration signal, use Gaussian filtering to remove Gaussian white noise from the radar. Preferably, the key parameter Σ for Gaussian filtering is 4.

[0043] S4.2. Transform the vibration signal into a point set Z in the IQ domain. Specifically, take the real part of the vibration signal as the x-coordinate and the imaginary part as the y-coordinate, and represent the vibration signal as a point set Z = {(x1,y1),(x2,y2),...,(x...} in the IQ domain. n ,y n )}.

[0044] S4.3. Based on the point set Z, use the circle center direction estimation algorithm to estimate the circle center direction vector v.

[0045] S4.4. Based on the estimated center direction vector v, the correction point set Z is obtained, and the correction point set Z′ is obtained.

[0046] S4.5 Extract the phase from the point set Z′.

[0047] S4.6. Improve the signal-to-noise ratio of the phase signal using a denoising algorithm. The denoising algorithm includes Gaussian filtering and high-pass filtering. In a preferred embodiment, the key parameter Σ of the Gaussian filter is 2; the cutoff frequency of the high-pass filter is set to 10Hz.

[0048] Furthermore, step S4.3, the circle center direction estimation algorithm, is specifically implemented through the following sub-steps:

[0049] S4.3.1 Input a point set Z = {(x1,y1),(x2,y2),...,(x...} of length n. n ,y n )}.

[0050] S4.3.2 Initialize the point set index i = 1, and the center direction vector v = 0;

[0051] S4.3.3 Calculate the average center point (x0, y0) of Z, and the center vector v1 = (-x0, -y0).

[0052] S4.3.4 Calculate the current index point (x) i ,y i ) and all points after it {(x i+1 ,y i+1 ),(x i+2 ,y i+2 ),...,(x n ,y n The distance L from )}

[0053]

[0054] S4.3.5 Calculate the first peak index sequence P of L, where peak index j = i + P, and peak index point (x j ,y j ).

[0055] S4.3.6 Calculate the current endpoint index {(x i ,y i ) and peak index point (x j ,y j The perpendicular vector v2 of )}:

[0056] v2=(y i -y j ,x j -x i (9)

[0057] S4.3.7 Update v by checking if the center vector v1, perpendicular vector v2, and center direction vector v are initialized:

[0058]

[0059] S4.3.8 Let i = j, repeat (4.3.4) to (4.3.8) until the peak index sequence P no longer exists.

[0060] S4.3.9 Output the center direction vector v.

[0061] Furthermore, the method for correcting the point set in step S4.4 is as follows: First, calculate the angle β between the center direction vector v and the average center vector v1, as shown in the following formula:

[0062]

[0063] Then, correct all points within the point set Z according to the following formula:

[0064]

[0065] Finally, the set of corrected points Z′={(x′1,y′1),(x′2,y′2),...,(x′...} is obtained. n ,y′ n )}.

[0066] Further, the method for extracting the phase in step S4.5 is as follows: First, the phase of Z′ is extracted using four-quadrant arctangent demodulation. The extraction formula is as follows:

[0067] p = arctan2(y,x) (13)

[0068]

[0069] Then, the extracted phase p is expanded. The phase expansion algorithm is as follows: For a given phase sequence p, starting from the second phase value, calculate the difference d between the current phase value and the previous phase value. If d is greater than 0, subtract d from the current point and all subsequent points; if d is less than 0, add d to the current point and all subsequent points. Repeat this process until the entire phase sequence p has been traversed.

[0070] Furthermore, the vibration feature extraction algorithm in step S5 is specifically implemented through the following sub-steps:

[0071] S5.1, Based on the target phase signal Φ a ,Φ b The amplitude is calculated using a peak search algorithm. Specifically:

[0072] First, the target phase signal is centered. The target phase signal is zero-mean processed using the following formula:

[0073] Φ′=Φ-μ (15)

[0074] Where Φ is the target phase signal and μ is the mean of the target phase signal.

[0075] Then, a peak search algorithm is used to retrieve the peak sequence of the target phase signal. The average value of the peak sequence is the amplitude characteristic.

[0076] S5.2. Based on the target phase signal Φ′, calculate the rotational speed using a spectrum analysis algorithm. Specifically, calculate the dominant frequency component f of Φ′ using a Fast Fourier Transform (FFT). Then, the rotational speed ω = f * 60.

[0077] S5.3. Based on the target phase signal Φ′, the eccentricity of the rotor trajectory is calculated using an ellipse fitting algorithm. The standard ellipse equation can be vectorized as follows:

[0078] F(H,X)=H*X=Ax 2 +Bxy+Dy 2 +Ex+Iy+J=0 (16)

[0079] Where H = [A, B, D, E, I, J], X = [x 2 ,xy,y 2 [x,y,1] T By using a radar target phase signal Φ a As the x-coordinate, another radar's target phase signal Φ b Using the y-coordinate as the reference, the target phase signals from the two radars are converted into a point set Q = {(x1,y1),(x2,y2),...,(x...}. n ,y n)}={X1,X2,...,X n The ellipse fitting problem can be modeled as the following optimization problem and optimized using the least squares method:

[0080]

[0081] Finally, the optimal value of H is obtained. Based on the solved H = [A, B, D, E, I, J], the eccentricity e of the ellipse is calculated:

[0082]

[0083] Finally, the eccentricity e is used as the eccentricity characteristic of the rotor trajectory.

[0084] Among them, A, B, D, I, and J are unknowns set by the standard formula of the ellipse equation from formula (16), which are solved by formula (17). For formula (18), these are known.

[0085] S5.4 Calculate the flatness of the rotor trajectory based on the point set Q using the convex hull algorithm. First, calculate the minimum bounding rectangle of the point set Q using the convex hull algorithm. Then, use the aspect ratio of the rectangle as the flatness feature of the rotor trajectory.

[0086] Furthermore, in step S6, the feature analysis algorithm employs a multi-threshold analysis method. Specifically, different thresholds are set for different rotational speeds.

[0087] This invention also discloses a rotational imbalance fault detection system based on non-contact rotor micro-displacement trajectory, which, based on the above method, includes the following modules:

[0088] Radar data acquisition module: Collects multi-directional surface vibration data caused by the operation of rotating equipment through a rotational imbalance fault detection device, i.e., raw data frames.

[0089] Vibration positioning module: Based on the original data frame, it uses a vibration positioning algorithm to analyze the position information of the vibration target.

[0090] Vibration signal extraction module: Based on the original data frame and the vibration target position information, the vibration phase extraction algorithm is used to extract the target vibration signal.

[0091] Phase correction module: Based on the target vibration signal, a phase correction algorithm is used to correct and extract the target phase signal.

[0092] Vibration feature extraction module: Based on the target phase signal, vibration features are extracted using a vibration feature extraction algorithm.

[0093] Feature Analysis Module: Based on the extracted vibration features, feature analysis algorithms are used to diagnose rotational imbalance faults.

[0094] The beneficial effects of this invention are:

[0095] (1) The technical solution proposed in this invention can simultaneously collect multi-directional surface vibration data of rotating equipment.

[0096] (2) The technical solution of the present invention relates to a phase correction algorithm that can remove DC offset in vibration signals.

[0097] (3) The technical solution of the present invention relates to a rotor trajectory detection algorithm, which can detect the rotor running trajectory of rotating equipment. Attached Figure Description

[0098] Figure 1 This is a flowchart of a preferred embodiment of the present invention, which describes a method for detecting rotational imbalance faults based on non-contact rotor micro-displacement trajectories.

[0099] Figure 2 This is a schematic diagram of the apparatus involved in the rotational imbalance fault detection method of the present invention.

[0100] Figure 3 This is a flowchart of the circle center direction estimation algorithm involved in this invention.

[0101] Figure 4 This is a diagram showing the results of the circle center direction estimation algorithm involved in this invention.

[0102] Figure 5 This is a diagram illustrating the extraction of rotor trajectory eccentricity features.

[0103] Figure 6 This is a block diagram of a rotational imbalance fault detection system based on non-contact rotor micro-displacement trajectory, according to a preferred embodiment of the present invention.

[0104] Figure 7 This is a circuit diagram of a radar synchronization system according to a preferred embodiment of the present invention. Detailed Implementation

[0105] The preferred embodiments of the present invention will now be described in detail with reference to the accompanying drawings.

[0106] like Figure 1 As shown in the figure, this embodiment of a method for detecting rotational imbalance faults based on non-contact rotor micro-displacement trajectories specifically includes the following steps:

[0107] S1. Radar data acquisition: Acquire multi-directional surface vibration data caused by the operation of rotating equipment, i.e., raw data frames.

[0108] like Figure 2As shown, the device involved in the rotational imbalance fault detection method includes two synchronized millimeter-wave radars and a radar synchronization system. The two synchronized millimeter-wave radars are placed radially on the rotor of the rotating equipment, with the radars at a 90-degree angle to each other relative to the rotor. Two of the millimeter-wave radars are existing commercial millimeter-wave radars (TI AWR1843), which can accept synchronization signals to control the radar's on / off state. The radar synchronization system was developed on an Arduino UNO R3, and its pin diagram is shown below. Figure 7 As shown. The radar synchronization system is connected to the two synchronized millimeter-wave radars by connecting the PB5 pin of the radar synchronization system, which serves as the synchronization signal output port, to the synchronization signal input pin of the millimeter-wave radar. The synchronization signal input pin of the millimeter-wave radar uses a high-level signal to enable the radar and a low-level signal to disable it. Therefore, by controlling the high / low level signal output of the PB5 pin of the radar synchronization system, the synchronous on / off state of the two radars can be controlled. Under the control of the radar synchronization system, the two millimeter-wave radars operate on different frequency bands and synchronously acquire data. In a preferred embodiment of the invention, the operating frequency bands of the two millimeter-wave radars are 77-79 GHz and 79-81 GHz, respectively. The millimeter-wave radar uses linear frequency modulated continuous wave (hereinafter referred to as chirp) signals as the transmission signals of the transmitting antenna. In a radar system with T receiving antennas, each antenna transmits N chirp signals, and the receiving antenna of the millimeter-wave radar acquires M sampling points for each chirp echo. Therefore, a frame of echo data with a size of [T, N, M] can be obtained. Finally, the original data frame obtained after accumulating multiple frames is a four-dimensional matrix of size [F,T,N,M], where F represents the number of frames in the target data frame (frame dimension), T represents the number of receiving antennas of the millimeter-wave radar (antenna dimension), N represents the number of chirp signals transmitted by the millimeter-wave radar (chirp dimension), and M represents the number of sampling points for each chirp echo (sampling dimension). In this embodiment, R is used to represent the original data frame of size [F,T,N,M]. This step explains the data representation at the algorithm level. A special note is needed. Figure 1 In this context, subscripts a and b are used to distinguish between two radars, for example, R a and R b These represent the raw data frames from the two radars. The data processing for the two radars is consistent in S2 through S4. For simplicity, subscripts 'a' and 'b' are not used to distinguish between the radars in S2 through S4.

[0109] S2. Vibration Localization: Based on the original data frame, a vibration localization algorithm is used to analyze the location information of the vibration target. The location information of the vibration target includes azimuth information (θ, r) and antenna weighting vector α. For example... Figure 3 As shown, the vibration localization algorithm is implemented through the following sub-steps:

[0110] S2.1 For the original data frame, calculate the target distance matrix. Specifically, first extract a single original data frame (of size [T, N, M]) from the original data frame, and then perform a Fast Fourier Transform (FFT) on the sampling dimension to obtain a target distance matrix of size [T, N, M]. The size of the data matrix remains unchanged before and after the FFT process. To distinguish the target distance matrix from a single original data frame, [T, N, R] is used to represent the target distance matrix. After the FFT process, targets at different distances are mapped to different frequency points, thus realizing the mapping from the sampling dimension (M) to the distance dimension (R).

[0111] S2.2. Based on the target range matrix, calculate the optimal antenna weighting coefficients and range-angle power spectrum. Based on the target range matrix that distinguishes environmental targets by range, use angle scanning to calculate the range-angle power spectrum.

[0112] First, the angle scan vector of the radar antenna is calculated. Assuming the millimeter-wave radar is fixed, angle scanning is achieved by weighting the vectors along the antenna dimension. The formula for calculating the radar antenna angle scan vector is shown below:

[0113] a(θ)=[a(θ1),a(θ2),...a(θ c (1)

[0114]

[0115] Where T is the number of radar receiving antennas, the spacing between the antenna arrays is d, λ is the signal wavelength, and θ = {θ1, θ2, ..., θ...} C} represents the angular scanning range.

[0116] Then, the angular distance power spectrum and the optimal antenna weighting coefficient are calculated, and the calculation formulas are shown below:

[0117]

[0118]

[0119] Where p is the angular range power spectrum, w is the optimal antenna weighting coefficient matrix, and x is the target range matrix. H The conjugate transpose, * T Indicates the transpose operation, * -1 This represents the inverse operation. After calculation, p is a [C,R] matrix of size n, and w is a [C,R,T] matrix of size n. Here, C is the angular dimension related to the angular scanning range θ, and R is the distance dimension. Then, for a given θ and r, w(θ,r) can be represented as the optimal weight vector at angle θ at distance r, with a length of T.

[0120] S2.3. Calculate vibration location information based on the angular distance power spectrum and the optimal antenna weighting coefficient matrix.

[0121] First, static background removal techniques are used to remove interference from static backgrounds in the environment. Specifically, this is achieved by subtracting the average power of the angular distance power spectrum to interfere with targets with insignificant power (static background). The calculation formula is shown below:

[0122]

[0123] Where p′ represents the calculation result, and p is the angular distance power spectrum. This indicates an average calculation.

[0124] Then, calculate the maximum index (C′,R′) of p′. (C′,R′) is used as the vibration orientation information, denoted as (θ,r).

[0125] Finally, the optimal antenna weighting coefficients for the vibrating target are extracted from the optimal antenna weighting coefficient matrix based on the vibration azimuth information. The extraction formula is shown below:

[0126] α=w(θ,r) (6)

[0127] In summary, this step calculates the position information of the vibrating target, including azimuth information (θ, r) and antenna weighting vector α, given a single raw data frame. It should be noted that for a stable environment, step S2 does not need to be repeated; that is, the azimuth information (θ, r) and antenna weighting vector α can be reused.

[0128] S3. Vibration Signal Extraction: Based on the original data frame and the vibration target location information, the target vibration signal is extracted using a vibration phase extraction algorithm. For example... Figure 3 As shown, the phase extraction algorithm is implemented through the following sub-steps:

[0129] S3.1. Reassemble the original data frames to obtain the original data groups. To ensure the continuity of the target vibration signal, the original data frames of size [F,T,N,M] are reassembled, and the F-frame data is reassembled into G groups of data frames.

[0130] First, identify the frame indices contained in each data group. Specifically, for the frame sequence [1,2,3,4,5,6,7,...,F], perform sliding window sampling with a window size of width and a step size of step to divide the frame sequence into G groups. Taking the frame sequence [1,2,3,4,5,6,7], width=3, step=2 as an example, the result after division is [1,2,3], [3,4,5], [5,6,7].

[0131] Then, each set of data is merged. Specifically, multiple frames within a set of data are stacked along the chirp dimension. Taking a set of data [3,T,N,M] as an example, it represents three data frames of size [T,N,M]. After stacking along the chirp dimension, the result is [T,3*N,M]. Since each set has partially overlapping data, the extracted target vibration signal is continuous.

[0132] S3.2. For the original data set, calculate the target distance matrix. Specifically, for the original data set of size [G,T,N,M], perform FFT (Fast Fourier Transform) processing on the sampling dimension to obtain the target distance matrix x of size [G,T,N,M]. This sub-step is processed in the same way as sub-step (2.1) in S2.

[0133] S3.3. Extract the target vibration signal based on the target distance matrix and the vibration target position information. The extraction process involves extracting the vibration signal from the target distance matrix x based on the vibration target position (θ, r) information. The extraction formula is as follows:

[0134] z=x(T,N,r)*α (7)

[0135] Where z is the extracted target vibration signal, r is the distance index of the vibration target location, α is the optimal antenna weighting coefficient, and x(T,N,r) is the target vibration signal at distance r extracted based on the distance index r. The vibration signal at distance r is weighted and summed using the optimal antenna weighting coefficient α to extract the vibration signal at angle θ.

[0136] For data frames G, the target vibration signal z can be represented as {z1, z2, ..., z}. G} Where G represents the number of groups, and each group is a target vibration signal of length N.

[0137] In summary, this step extracts the target vibration signal z under the given conditions of the original data frame and the position information of the vibration target.

[0138] S4. Phase Correction: Based on the target vibration signal, a phase correction algorithm is used to correct and extract the target phase signal. Ideally, the vibration signal in the IQ domain (IQ domain == IQ plane == complex plane) represents an arc with its center at the origin. Due to DC offset, the center of the vibration signal in the IQ domain is not at the origin. Therefore, phase correction is required. The phase correction algorithm corrects the phase of each sub-vibration signal z within the target vibration signal z. i Corrected to Specifically, this is achieved through the following sub-steps:

[0139] S4.1. For a given sub-vibration signal, use Gaussian filtering to remove Gaussian white noise from the radar. Gaussian white noise disrupts the distribution of the target vibration signal in the IQ domain, affecting the estimation of the center direction vector. The key parameter for Gaussian filtering is Σ = 4.

[0140] S4.2. Transform the vibration signal into a point set Z in the IQ domain. Specifically, take the real part of the vibration signal as the x-coordinate and the imaginary part as the y-coordinate, and represent the vibration signal as a point set Z = {(x1,y1),(x2,y2),...,(x...} in the IQ domain. n ,y n )}.

[0141] S4.3. Based on the point set Z, estimate the center direction vector v using the center direction estimation algorithm. The flowchart of the center direction estimation algorithm is as follows: Figure 3 As shown. The process is briefly described below:

[0142] S4.3.1 Input a point set Z = {(x1,y1),(x2,y2),...,(x...} of length n. n ,y n )}.

[0143] S4.3.2 Initialize the point set index i = 1 and the center direction vector v = 0.

[0144] S4.3.3 Calculate the average center point of Z (x0, y0), v1 = (-x0, -y0).

[0145] S4.3.4, Calculate (x) i ,y i The distance L between it and all points thereafter:

[0146]

[0147] S4.3.5 Calculate the peak index sequence P of L, j = i + P[0].

[0148] S4.3.6, Calculate endpoint pairs {(x i ,y i ),(x j ,y j The perpendicular vector v2 of )}:

[0149] v2=(y i -y j ,x j -x i (9)

[0150] S4.3.7 Update v by checking v1, v2, and whether v is initialized:

[0151]

[0152] S4.3.8 Let i = j, repeat (4.3.4) to (4.3.8) until the peak index sequence P no longer exists.

[0153] S4.3.9 Output the center direction vector v.

[0154] This process estimates the center direction by calculating the direction perpendicular to the movement of the midpoint Z. There are two perpendicular directions; the direction vector forming an acute angle with v1 is selected as the currently estimated center direction vector. For multiple estimated direction vectors, their average direction vector is taken as the final center direction estimate. The vibration of the rotating equipment can be viewed as a reciprocating motion, corresponding to the movement of the midpoint Z from one end of the arc to the other. The direction perpendicular to the movement of the midpoint Z is the center direction. The center direction estimation result is as follows: Figure 4 As shown.

[0155] S4.4. Based on the estimated center direction vector v, the correction point set Z is obtained, and the correction point set Z′ is obtained.

[0156] First, calculate the angle β between the center direction vector v and the average center vector v1, as shown in the following formula:

[0157]

[0158] Then, correct all points within the point set Z according to the following formula:

[0159]

[0160] Finally, the set of corrected points Z′={(x′1,y′1),(x′2,y′2),...,(x′...} is obtained. n ,y′ n )}.

[0161] S4.5 Extract the phase from the point set Z′. The specific steps are as follows:

[0162] First, the phase of Z′ is extracted using four-quadrant arctangent demodulation. The extraction formula is shown below:

[0163] p = arctan2(y,x) (13)

[0164]

[0165] Here, p is the extracted phase, and x and y are points within the point set Z′. The range of the arctangent function in the four quadrants is (-π, π), therefore, p will fold within this interval, resulting in a truncated phase. For example, suppose...

[0166] p = [130°, 150°, 170°, -10°], where a phase jump occurs clearly from 170° to -10°.

[0167] Then, the extracted phase p is unfolded. Folded phases cause phase jumps. To obtain continuous phases, the folded phases need to be unfolded. The phase unfolding algorithm is as follows: For a given phase sequence p, starting from the second phase value, calculate the difference d between the current phase value and the previous phase value. If d is greater than 0, subtract d from the current point and all subsequent points; if d is less than 0, add d to the current point and all subsequent points. Repeat this process until the entire phase sequence p has been traversed.

[0168] S4.6. Use denoising algorithms to improve the signal-to-noise ratio (SNR) of the phase signal. Denoising algorithms include Gaussian filtering and high-pass filtering. Gaussian filtering can remove Gaussian white noise from the radar, improving the SNR of the phase signal. High-pass filtering can effectively remove the trend term of the phase signal. In this embodiment, the key parameter Σ for Gaussian filtering is 2; the cutoff frequency for high-pass filtering is set to 10Hz.

[0169] Each set of target vibration signals z is corrected using a phase correction algorithm (steps S4.1 to S4.6). i For target phase signal Yes, the target phase signal can be obtained. Finally, the target phase signal is unfolded into one-dimensional data to obtain a target phase signal of size [G*N], denoted as Φ.

[0170] In summary, the target phase signal Φ is extracted from the target vibration signal z using a phase correction algorithm.

[0171] S5. Vibration Feature Extraction: Based on the target phase signal, vibration features are extracted using a vibration feature extraction algorithm. These vibration features mainly include amplitude, rotational speed, rotor trajectory eccentricity, and rotor trajectory flatness. The extraction of vibration features is specifically achieved through the following sub-steps.

[0172] S5.1, Based on the target phase signal Φ a ,Φ b The amplitude is calculated using a peak search algorithm.

[0173] First, the target phase signal is centered. The target phase signal is processed to have a zero mean using the following formula.

[0174] Φ′=Φ-μ (15)

[0175] Where Φ is the target phase signal and μ is the mean of the target phase signal.

[0176] Then, a peak search algorithm is used to retrieve the peak sequence of the target phase signal. The average value of the peak sequence is the amplitude characteristic.

[0177] S5.2. Based on the target phase signal Φ′, calculate the rotational speed using a spectrum analysis algorithm. Specifically, calculate the dominant frequency component f of Φ′ using a Fast Fourier Transform (FFT). Then, the rotational speed ω = f * 60.

[0178] S5.3. Based on the target phase signal Φ′, the eccentricity of the rotor trajectory is calculated using an ellipse fitting algorithm. The standard ellipse equation can be vectorized as follows:

[0179] F(H,X)=H*X=Ax 2 +Bxy+Dy 2 +Ex+Iy+J=0 (16)

[0180] Where H = [A, B, D, E, I, J], X = [x 2 ,xy,y 2 [x,y,1] T By using a radar target phase signal Φ a As the x-coordinate, another radar's target phase signal Φ b Using the y-coordinate as the reference, the target phase signals from the two radars are converted into a point set Q = {(x1,y1),(x2,y2),...,(x...}. n ,y n )}={X1,X2,...,X n The ellipse fitting problem can be modeled as the following optimization problem and optimized using the least squares method:

[0181]

[0182] Finally, the optimal value of H is obtained. Based on the solved H = [A, B, D, E, I, J], the eccentricity e of the ellipse is calculated:

[0183]

[0184] Finally, the eccentricity e is used as the eccentricity feature of the rotor trajectory. The fitted ellipse and eccentricity are as follows: Figure 5 As shown.

[0185] S5.4 Calculate the flatness of the rotor trajectory based on the point set Q using the convex hull algorithm. First, calculate the minimum bounding rectangle of the point set Q using the convex hull algorithm. Then, use the aspect ratio of the rectangle as the flatness feature of the rotor trajectory.

[0186] S6. Feature Analysis: Based on the extracted vibration features, a feature analysis algorithm is used to diagnose rotational imbalance faults. For normally operating rotating equipment, its vibration parameters (amplitude, rotational speed, eccentricity of the rotor trajectory, and flatness) are within a certain range. When a rotational imbalance fault occurs, these vibration parameters will change abruptly. For example, a common phenomenon of rotational imbalance faults is that the overall amplitude of the equipment increases exponentially. By combining threshold values ​​and detecting the trend of vibration parameter changes, rotational imbalance faults can be identified. This is based on the centrifugal force formula F = m * v. 2 As can be seen from the figure, centrifugal force is directly proportional to rotational speed. The degree of change in vibration parameters varies at different rotational speeds. When the rotational speed is low, even if a rotational imbalance fault exists, the change in vibration parameters is not significant. Therefore, when considering the influence of rotational speed, different threshold values ​​can be set for different rotational speeds.

[0187] like Figure 6 As shown, this preferred embodiment discloses a rotational imbalance fault detection system based on non-contact rotor micro-displacement trajectory, which, based on the above-described method embodiment, includes the following modules:

[0188] Radar data acquisition module: Collects multi-directional surface vibration data caused by the operation of rotating equipment through a rotational imbalance fault detection device, i.e., raw data frames.

[0189] Vibration positioning module: Based on the original data frame, it uses a vibration positioning algorithm to analyze the position information of the vibration target.

[0190] Vibration signal extraction module: Based on the original data frame and the vibration target position information, the vibration phase extraction algorithm is used to extract the target vibration signal.

[0191] Phase correction module: Based on the target vibration signal, a phase correction algorithm is used to correct and extract the target phase signal.

[0192] Vibration feature extraction module: Based on the target phase signal, vibration features are extracted using a vibration feature extraction algorithm.

[0193] Feature Analysis Module: Based on the extracted vibration features, feature analysis algorithms are used to diagnose rotational imbalance faults.

[0194] Other aspects of this embodiment can be found in the above method embodiments.

[0195] This invention innovatively proposes to use two synchronous millimeter-wave radars to track the rotor trajectory, calculate the eccentricity and flatness of the rotor trajectory, and combine the amplitude and rotational speed to diagnose the rotational imbalance fault. This invention has the advantages of high accuracy and low cost.

[0196] Those skilled in the art will recognize that various substitutions and modifications are possible without departing from the spirit and scope of the invention and the appended claims. Therefore, the scope of protection of the invention should not be limited to the content disclosed in the embodiments.

Claims

1. A method for detecting rotational imbalance faults based on non-contact rotor micro-displacement trajectories, characterized in that, Includes the following steps: S1. Place two synchronous millimeter-wave radars in the radial direction of the rotor of the rotating equipment and make a 90-degree angle between the radars relative to the rotor to collect multi-directional surface vibration data caused by the operation of the rotating equipment, i.e., raw data frames. S2. Based on the original data frame, the vibration positioning algorithm is used to analyze the position information of the vibration target; The vibration localization algorithm is implemented through the following sub-steps: S2.1 For the original data frame, calculate the target distance matrix; S2.

2. Based on the target range matrix, calculate the optimal antenna weighting coefficients and range-angle power spectrum; S2.

3. Calculate vibration location information based on the angular distance power spectrum and the optimal antenna weighting coefficient matrix; S3. Based on the original data frame and the vibration target position information, the vibration phase extraction algorithm is used to extract the target vibration signal; The vibration phase extraction algorithm is implemented through the following sub-steps: S3.1 Perform frame reassembly on the original data frames to obtain the original data groups; S3.2 For the original data set, calculate the target distance matrix; S3.3 Extract the target vibration signal based on the target distance matrix and the vibration target position information; S4. Based on the target vibration signal, a phase correction algorithm is used to convert a set of target vibration signals... Correct and extract the target phase signal; The phase correction algorithm described herein is as follows: S4.

1. For a given sub-vibration signal, use Gaussian filtering to remove Gaussian white noise from the radar; S4.2, Converting the vibration signal to a point set in the IQ domain ; S4.3, Point Set Based The center direction vector is estimated using a circle center direction estimation algorithm. ; The specific algorithm for estimating the center direction is as follows: S4.3.1 Input a point set of length n S4.3.2 Initialize the point set index , center direction vector =0; S4.3.3, Calculation average center point , center vector ; S4.3.4 Calculate the current index point All points after it distance : （8） S4.3.5, Calculation The first peak index Peak Index Peak index point ; S4.3.6 Calculate the current index point and peak index point vertical vector : （9） S4.3.7, Through the central vector and vertical vector and the center direction vector To update, initialize or not? : （10） S4.3.8, Order Return to step S4.3.4 until the peak index is reached. It does not exist; S4.3.9 Output the center direction vector ; S4.4, Based on the estimated center direction vector Calibration point set The set of correction points is obtained. ; This step S4.4 is implemented through the following steps: S4.4.1, Calculate the center direction vector. With the average center vector The included angle β between them is calculated using the following formula: S4.4.2, the correction point set is performed according to the following formula. All points inside: S4.4.3, obtain the set of correction points. S4.5, Pair Set Phase extraction; S4.

6. Use denoising algorithms to improve the signal-to-noise ratio of the phase signal; S5. Based on the target phase signal, vibration features are extracted using a vibration feature extraction algorithm; The vibration feature extraction algorithm is implemented through the following sub-steps: S5.1, Based on the target phase signal The amplitude is calculated using a peak search algorithm: First, the target phase signal is centered using the following formula to achieve a zero mean: （15） in, It is the target phase signal. It is the mean of the target phase signal; Then, the peak sequence of the target phase signal is retrieved using the peak search algorithm, and the average value of the peak sequence is the amplitude characteristic. S5.2, Based on the target phase signal The rotational speed is calculated using a spectrum analysis algorithm; S5.3, Based on target phase signal The eccentricity of the rotor trajectory is calculated using an ellipse fitting algorithm; the standard ellipse equation can be vectorized into the following expression: （16） in , By using a radar target phase signal As the x-coordinate, the target phase signal of another radar Using the y-coordinate, the target phase signals from the two radars are converted into a point set. The ellipse fitting problem can be modeled as the following optimization problem and optimized using the least squares method: Based on the solution Calculate the eccentricity of the ellipse : in, ) are the coordinates of the center of the ellipse, a is the major semi-axis of the ellipse, and b is the minor semi-axis of the ellipse; Formula (18) is obtained through Calculate the coordinates of the center of the ellipse Then calculate the major semi-axis 'a' and the minor semi-axis 'b' of the ellipse, and finally calculate the eccentricity of the ellipse. ; Eccentricity of the fitted ellipse Eccentricity characteristics of the rotor trajectory; S5.

4. Based on the point set Q, the convex hull algorithm is used to calculate the rotor trajectory flatness; S6. Based on the extracted vibration features, a feature analysis algorithm is used to diagnose rotational imbalance faults.

2. The method for detecting rotational imbalance faults based on non-contact rotor micro-displacement trajectories according to claim 1, characterized in that, In step S1, the two millimeter-wave radars operate on different frequency bands and collect data synchronously under the control of the radar synchronization system.

3. The method for detecting rotational imbalance faults based on non-contact rotor micro-displacement trajectories according to claim 1, characterized in that, In step S1, the original data frame is a four-dimensional matrix of size [F,T,N,M], where F represents the number of frames in the target data frame, i.e., the frame dimension, T represents the number of receiving antennas of the millimeter-wave radar, i.e., the antenna dimension, N represents the number of chirp signals transmitted by the millimeter-wave radar, i.e., the chirp dimension, and M represents the number of sampling points for each chirp echo, i.e., the sampling dimension.

4. The method for detecting rotational imbalance faults based on non-contact rotor micro-displacement trajectories according to any one of claims 1-3, characterized in that, In step S6, the feature analysis algorithm uses a multi-threshold analysis method.

5. A rotational imbalance fault detection system based on non-contact rotor micro-displacement trajectory, wherein the system is based on the method described in any one of claims 1-4, characterized in that... Includes the following modules: Radar data acquisition module: Collects multi-directional surface vibration data caused by the operation of rotating equipment through a rotational imbalance fault detection device, i.e., raw data frames; Vibration positioning module: Based on the original data frame, it uses a vibration positioning algorithm to analyze the position information of the vibration target; Vibration signal extraction module: Based on the original data frame and vibration target position information, the vibration phase extraction algorithm is used to extract the target vibration signal; Phase correction module: Based on the target vibration signal, it uses a phase correction algorithm to correct and extract the target phase signal; Vibration feature extraction module: Based on the target phase signal, it uses a vibration feature extraction algorithm to extract vibration features; Feature analysis module: Based on the extracted vibration features, Feature analysis algorithms are used to diagnose rotational imbalance faults.

Images