Bearing residual service life prediction method based on dynamic graph neural network

Abstract

The invention discloses a bearing residual service life prediction method based on a dynamic graph neural network, and the method comprises the following steps: building a bearing dynamical model based on a Hertz contact theory, which is used for describing an expansion process of a defect of a bearing outer ring and calculating contact deformation; taking the defect size as a variable, and using multi-objective optimization to match the dynamic response of the digital twin of the bearing and the real bearing to obtain a preliminary full life cycle defect size curve; extracting a high-fidelity defect size through a bearing time domain vibration signal, calculating a calibration defect size, and performing two-stage updating of a life cycle defect size; a bidirectional long-short-term memory network is used as a training model for training; and taking the trained feature real-time mapping network model as a feature supplement tool, taking the mapped defect size as an additional feature, and inputting the supplemented feature set data into the double-correlation dynamic graph neural network to obtain a final residual service life prediction result, thereby improving the model prediction performance.

Description

Technical Field

[0001] The present invention relates to the technical field of bearing remaining service life prediction, and in particular to a bearing remaining service life prediction method based on a dynamic graph neural network. Background Art

[0002] Bearings play a vital role in industrial equipment and are widely used in industries such as rail transportation, aerospace, and wind power generation. However, most bearings in industrial equipment operate under harsh conditions, making them susceptible to corrosion, wear, and other failures. Bearing failures can affect the normal operation of the system and, in more serious cases, lead to major safety incidents. Accurate and efficient remaining service life prediction can guide industrial equipment maintenance, develop reasonable maintenance plans, and avoid the economic losses and casualties caused by unexpected bearing failure. Therefore, to improve the accuracy of bearing remaining service life prediction, both industry and academia are actively exploring and researching better prediction methods.

[0003] Existing deep learning-based bearing remaining service life prediction methods do not observe the motion behavior of the bearing during degradation and ignore the real-time interaction between the feature space and the bearing physical space. Summary of the Invention

[0004] The purpose of the present invention is to provide a method for predicting the remaining useful life of bearings based on a dynamic graph neural network. The real-time correlation between physical space and digital space is captured through a multi-objective optimized digital twin, and the model parameters are updated to realize the interaction between the twin and the actual bearing. The defect curve is further updated using high-fidelity calibration defects in two stages. Real-time feature mapping is realized by establishing a model based on a bidirectional long short-term memory network, and a dynamic graph neural network with dual correlation is used to improve the model's ability to extract correlation, thereby improving the model prediction performance.

[0005] To achieve the above object, the present invention provides a method for predicting the remaining service life of a bearing based on a dynamic graph neural network, comprising the following steps:

[0006] Step S1: Establish a bearing dynamics model based on Hertz contact theory to describe the expansion process of the bearing outer ring defect and calculate the contact deformation;

[0007] Step S2: Using the defect size as a variable, multi-objective optimization is used to match the dynamic response of the bearing digital twin and the real bearing to obtain a preliminary full life cycle defect size curve;

[0008] Step S3: extracting high-fidelity defect size from the bearing time-domain vibration signal to calculate the calibrated defect size, and then perform a two-stage update of the lifecycle defect size;

[0009] Step S4: Using the obtained full life cycle defect size as a label dataset and the root mean square value of the real bearing vibration signal as input, a bidirectional long short-term memory network is used as a training model to obtain a feature mapping network model capable of real-time feature mapping;

[0010] Step S5: Collect vibration data during the operation of the bearing and convert its samples into feature channels. Use the trained feature real-time mapping network model as a feature supplement tool, use the mapped defect size as an additional feature to supplement the feature set data of the real bearing, and input the supplemented feature set data into the dual-correlation dynamic graph neural network to obtain the final remaining service life prediction result.

[0011] The above solution is further preferred, wherein in step 1, establishing a bearing dynamics model based on Hertz contact theory to describe the expansion process of the bearing outer ring defect comprises the following steps:

[0012] Step S0101: Introduce Hertz contact theory to establish a dynamic model of the bearing. The established dynamic model satisfies:

[0013]

[0014] Where M is the mass of the system, C is the damping coefficient, K is the stiffness, and F is the excitation load. x represents the acceleration, velocity and displacement of the system respectively;

[0015] Step S0102: Newton's second law is used to derive a two-degree-of-freedom dynamic equation for a rolling bearing that takes into account damping and a constant external load, which is expressed as the following equation:

[0016]

[0017] Among them, M b is the equivalent mass of the bearing, C b is the equivalent damping of the bearing, and is the vibration acceleration in the x and y directions, and is the vibration speed, F Cx and F Cy is the contact force component between the rolling element and the outer ring in the x and y directions, Q x and Q y is the bearing load component;

[0018] Step S0103: All contact force components are considered as nonlinear contact forces, and the following equation is calculated using Hertz contact theory:

[0019]

[0020] Where K represents the total contact stiffness of the bearing, N represents the number of rolling elements, μ is the factor that determines whether the roller is in contact with the outer ring, and θ i Indicates the angular position of rolling element i, for cylindrical roller bearings For ball bearings

[0021] Step S0104: Substitute the equation in step 1.3 into the equation in step 1.2 to calculate and obtain the final defect expansion model of the rolling bearing twin, which is expressed by the following equation:

[0022]

[0023] Where K represents the total contact stiffness of the bearing, N represents the number of rolling elements, and μ is a factor that determines whether the roller is in contact with the outer ring. μ satisfies the following expression:

[0024]

[0025] Among them, Δ i Indicates the contact deformation between rolling element i and raceway, Δ i Satisfies the following equation:

[0026] Δ i =xcosθ i +ysinθ i -cH;

[0027] Where H represents the time-varying displacement excitation, θ i represents the angular position of rolling element i at time t;

[0028] Step S0105: Calculate the angle position θ i First, convert the angular velocity of the shaft into radians, then calculate the speed of the bearing cage, and finally calculate the position of the rolling element i. The specific calculation process is as follows:

[0029]

[0030] Among them, n s Indicates the shaft speed, ω s represents the angular velocity of the shaft, d represents the diameter of the rolling element, D p Indicates the pitch diameter of the bearing, ω c represents the speed of the bearing cage, and Z represents the number of rolling elements;

[0031] Step S0106: Calculate the time-varying displacement excitation ΔH of the bearing over its entire life cycle in two stages.

[0032] The above solution is further preferred, wherein the two-stage calculation of the time-varying displacement excitation ΔH of the bearing over its entire life cycle includes the following steps:

[0033] Step 3.1: First, simulate the displacement excitation generated by the evolution of a real bearing defect. To simplify the dynamic model and ignore the contact deformation between the bearing rolling elements and the edge of the bearing outer ring defect, the bearing outer ring defect is considered as a rectangle with a length of L and a depth of H, and a two-stage displacement excitation model is established.

[0034] At the initial stage of defect generation, the rolling element cannot completely contact the bottom of the defect, and its maximum displacement excitation H max Always smaller than the defect depth H, H max Calculated using the following expression:

[0035] H max =0.5d-((0.5d) 2 -(0.5L) 2 );

[0036] At this time, the time-varying displacement excitation function satisfies the following expression:

[0037]

[0038] Among them, b i is the residual angular position of rolling element i at time t, b i =mod(θ i ,2π),θ d is the angle corresponding to the defect center, and δ1 represents the span angle corresponding to the defect;

[0039] Step 3.2: In the defect expansion stage, as the defect size increases, the maximum displacement excitation gradually increases until it is equal to the defect depth H. At this time, the rolling element will completely fall into the defect. At this time, the time-varying displacement excitation function in the defect expansion stage satisfies the following expression:

[0040]

[0041] Among them, mod(·) represents the remainder operation, b i is the residual angular position of rolling element i at time t,

[0042] b i =mod(θ i ,2π),θ d is the angle corresponding to the defect center, δ1 represents the span angle corresponding to the defect, and δ2 represents the span angle of the ball center when the rolling element contacts the bottom of the defect.

[0043] The above solution is further preferred. In step 2, the defect size is used as a variable, and multi-objective optimization is used to match the dynamic response of the twin and the real bearing to obtain a preliminary full life cycle defect size curve, including the following steps:

[0044] Step 2.1: Use multi-objective optimization matching to iterate the parameters of the twins. The optimization target takes into account the similarity of time domain and frequency domain to ensure the best matching defect parameters.

[0045] The matching expression satisfies the following:

[0046]

[0047] Among them, λ op represents the optimized defect parameters, including defect length L and defect depth H; T and P represent the dynamic responses of the twin and the real bearing, respectively; F(·) represents fast Fourier transform; dtw(·) represents DTW distance; pearson(·) represents Pearson correlation coefficient;

[0048] Step 2.2: Use the particle swarm algorithm to iteratively update the position of the group, guide it to converge towards the optimal solution, and finally find the global optimal solution. The specific update process is expressed as follows:

[0049]

[0050] in, represents the velocity vector of particle k after the tth iteration, represents the position vector of particle k after the tth iteration, τ represents the inertia coefficient, η is a random number in [0,1]; α1 represents the individual learning factor, α2 represents the social learning factor, represents the optimal position of particle i in t iterations, gbest t Indicates the optimal position in the group;

[0051] Step 2.3: After calculating the velocity of the t+1th iteration, the position of each particle is updated according to the following formula. The position of the particle represents the size of the defect:

[0052]

[0053] The above solution is further preferred, in step S3, the two-stage update of the defect size of the bearing throughout its life cycle includes the following steps:

[0054] Step S301: Use the bearing time-domain vibration signal to describe the physical process of the rolling element passing through a local defect at a certain moment, which is expressed by the following equation:

[0055] (OF·cosα1-Rcosα2) 2 +(OF·sinα1-Rsinα2-R·ω c ·t p ) 2 =r b2 ;

[0056] Among them, α1 represents the span angle between the rolling body center point A and point E when the rolling body just contacts the defect entrance E, α2 represents the span angle between the rolling body center point C and point E when the rolling body just leaves the defect entrance E, t p The time it takes for the rolling element to rotate from point C to just contact the defect outlet point F is the time it takes for the bearing to rotate, r b is the radius of the rolling element;

[0057] Step S302: Extract the position of each point of the bearing time domain vibration signal to obtain the span angle α1 of the defect, and calculate the bearing calibration defect size U according to the span angle t , bearing calibration defect size U t Use the following expression for calculation:

[0058]

[0059] Step S303: Apply steps S301-S302 to the vibration signal at the selected moment to obtain the calibration defect size at that moment, and apply the following expression to the full life cycle bearing defect curve generated in one stage:

[0060] The final defect size of the two-stage update model is calculated using the formula:

[0061]

[0062] Among them, L t represents the first stage optimization defect corresponding to the calculated first calibration defect, n represents the total time of bearing failure, γ is the number of calibration defects obtained, U t represents the calibration defect size at the tth moment.

[0063] The above solution is further preferred, in step 4, a bidirectional long short-term memory network is trained,

[0064] The steps to obtain a network model that can perform real-time feature mapping include:

[0065] Step S41, using a bidirectional long short-term memory network to establish a connection between the virtual and real spaces, and calculating the root mean square value of the bearing's real-time vibration signal;

[0066] Step S42: Obtain the defect size of the bearing throughout its life cycle and use it as a training data set.

[0067] Step S43: Using the root mean square value of the bearing throughout its life cycle as input data and the defect size of the bearing throughout its life cycle as output data, a bidirectional long short-term memory network is trained, and a well-performing feature mapping model is obtained through verification;

[0068] Step S44: Input the real-time root mean square value of the bearing into the trained feature mapping network model, and then use the trained model to output the real-time defect size of the twin, and quickly feed it back to the real bearing to guide its remaining service life prediction.

[0069] The above solution is further preferred, in step 5, inputting the supplemented feature set into the dual-correlation dynamic graph neural network to obtain the final remaining useful life prediction result, including the following steps:

[0070] Step S501: Extract the characteristics of the vibration data during the operation of the bearing, use the root mean square value of the entire life cycle of the bearing as the input feature, and input it into the trained bidirectional long short-term memory network to obtain the twin feature, and then fuse the twin feature with the characteristics of the real bearing;

[0071] Step S502: Using a moving window with a stride of 1 and a length of L, sample the feature set data that concentrates on the defect features of the bearing digital twin;

[0072] Step S503: Using the columns of the window data as graph node vectors, a spatial correlation graph is established; using the rows of the window data as graph node vectors, a temporal correlation graph is established;

[0073] Step S504: Use cosine similarity to define the similarity matrix between graph node vectors. The cosine similarity matrix C (i,n) Elements in It is calculated by the following expression:

[0074]

[0075] in represents the f-th row feature of the n-th sampling data of the i-th bearing, ||·|| is the second norm of the vector;

[0076] Step S505: Apply function o(·) to the cosine similarity matrix C (i,n) The elements in can be used to get the adjacency matrix Elements in δ is a threshold defined artificially, and the function o(·) is expressed as:

[0077]

[0078] Then, the normalized graph Laplacian matrix is defined by the obtained adjacency matrix A, which is expressed as follows:

[0079]

[0080] Where D represents the degree matrix, and the expression is I Nrepresents the identity matrix;

[0081] Step S506: In order to obtain the correlation information of the feature set data through the adjacency matrix, the information transfer process of the graph convolutional neural network model is defined. This process can update the feature set data, which is expressed as follows:

[0082]

[0083] in, Is an adjacency matrix with self-loop, that is, a self-connection is added to the original adjacency matrix, W (L) represents the trainable weight matrix of the Lth layer, σ(·) represents the ReLU activation function, and H (L) Represents the input of the Lth layer. When L=0, H (0) =0;

[0084] Step S507: Input the fused feature channel data into the dual-correlation dynamic graph neural network to extract spatiotemporal correlation. The dual-correlation dynamic graph neural network update mechanism is as follows:

[0085]

[0086] Where σ represents the Sigmiod activation function, F(·) represents the graph convolution process, W and b represent the trainable weight matrix and bias, i t 、f t 、o t are the input gate, forget gate, and output gate at the tth position, respectively, c t is the cell state at the tth position, h t is the hidden layer information at the tth position;

[0087] In step S508, after the spatiotemporal correlation is extracted using the dual-correlation dynamic graph neural network, the updated feature set data is flattened and input into the fully connected layer to obtain the final remaining useful life prediction result, which is expressed as follows:

[0088] y t =σ(W t ·Flatten(g t ));

[0089] Among them, y t Represents the prediction result, σ represents the ReLU activation function, W t Represents the trainable weight matrix, g t Indicates the updated feature set data.

[0090] In summary, since the present invention adopts the above technical solution, the present invention has the following beneficial technical effects:

[0091] (1) This paper proposes a two-stage updated digital twin model. Unlike traditional dynamic modeling, it can capture the real-time correlation between two types of spaces through multi-objective optimization, update the model parameters to realize the interaction between the twin and the actual bearing, and further update the defect curve using high-fidelity calibration defects in two stages.

[0092] (2) In the present invention, a real-time mapping model based on a bidirectional long short-term memory network is established, and the defect size obtained by the two-stage updating digital twin model is used as a training data set.

[0093] (3) In the present invention, a dynamic graph neural network based on dual correlation is proposed, which can extract the dynamic spatiotemporal correlation of feature channel data, thereby better extracting the correlation of features between physical space and digital space.

[0094] (4) In the present invention, a dynamic gated graph convolution layer is proposed, which can capture the current spatiotemporal correlation while retaining effective historical information, and also includes the temporal and spatial correlation information of the fused feature channels, thereby improving the stability of the model.

[0095] (5) The prediction model designed by the present invention based on two-stage updated digital twins and dual-correlation dynamic graph neural network has the significant advantages of high prediction accuracy, strong prediction robustness and good generalization performance. It can meet the practical application needs of bearing life prediction in predictive maintenance and has broad application potential in other fields. BRIEF DESCRIPTION OF THE DRAWINGS

[0096] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0097] Figure 1 It is a general flow chart of the digital twin architecture of the present invention.

[0098] Figure 2 This is a schematic flow chart of the steps of a bearing remaining service life prediction method based on a two-stage updated digital twin and a dual-correlation dynamic graph neural network of the present invention.

[0099] Figure 3 It is a schematic diagram of the two-degree-of-freedom twin bearing dynamics model of the present invention.

[0100] Figure 4 It is a schematic diagram of the defect expansion of the twin bearing outer ring of the present invention.

[0101] Figure 5 It is a schematic diagram describing the physical process of the defect motion of the outer ring of the rolling element of a twin bearing provided by the present invention.

[0102] Figure 6 It is a flow chart of the dual-correlation dynamic graph neural network model of the present invention.

[0103] Figure 7 It is a schematic diagram of the dynamic graph neural network module of the present invention.

[0104] Figure 8 It is a dynamic response effect diagram of the twin model in an embodiment of the present invention.

[0105] Figure 9 This is an envelope spectrum analysis comparison diagram of a twin bearing in an embodiment of the present invention and a real bearing.

[0106] Figure 10 1 is a first-stage bearing defect size curve obtained based on multi-objective optimization in an embodiment of the present invention.

[0107] Figure 11 1 is a bearing defect size curve diagram after averaging processing in an embodiment of the present invention.

[0108] Figure 12 1 is a graph of bearing defect size after two-stage updating using calibration defects in an embodiment of the present invention.

[0109] Figure 13 3 is a schematic diagram comparing the full life cycle defect sizes obtained by the embodiment of the present invention and other methods in the literature.

[0110] Figure 14 This is an effect diagram of feature mapping using a bidirectional long short-term memory network in an embodiment of the present invention.

[0111] Figure 15 This is a comparison chart of the score index effects of various models of the XJTU-SY dataset in the S2 scenario in an embodiment of the present invention.

[0112] Figure 16 This is a comparison chart of the score index effects of various models of the XJTU-SY dataset in the S1 scenario in an embodiment of the present invention.

[0113] Figure 17 3. It is a schematic diagram comparing the prediction effects of various models of the XJTU-SY dataset under the S1 scenario in an embodiment of the present invention.

[0114] Figure 18 3. It is a schematic diagram comparing the prediction effects of various models of the XJTU-SY dataset under the S2 scenario in an embodiment of the present invention. DETAILED DESCRIPTION

[0115] The following describes embodiments of the present invention in detail, examples of which are shown in the accompanying drawings, wherein the same or similar reference numerals throughout represent the same or similar elements or elements having the same or similar functions. The embodiments described below with reference to the accompanying drawings are exemplary and are intended to be used to explain the present invention, and are not to be construed as limiting the present invention.

[0116] The following are annotations for the corresponding English abbreviations to facilitate subsequent descriptions:

[0117] RUL: Remaining useful life;

[0118] DGCN: Dynamic Graph Neural Network module;

[0119] DT: Digital Twin;

[0120] Combined diagram Figure 1 and Figure 2 The present invention provides a bearing remaining service life prediction method based on a two-stage updated digital twin and a dual-correlation dynamic graph neural network, comprising the following steps:

[0121] Step S1: Establish a bearing dynamics model based on Hertz contact theory to describe the expansion process of the bearing outer ring defect and calculate the contact deformation of the bearing;

[0122] Step S2: Using the defect size as a variable, multi-objective optimization is used to match the dynamic response of the bearing digital twin and the real bearing to obtain a preliminary full-life cycle defect size curve. Using the defect size as the variable x, multi-objective optimization continuously changes the value of x to match the optimal fitness.

[0123] Step S3: Calculate the calibration defect size using theoretical methods from existing literature. Extract high-fidelity defect size from the bearing's time-domain vibration signal to calculate the calibration defect size. Perform a two-stage update of the full life cycle defect size based on this. Update the curve obtained in step 2 again to obtain a full life cycle defect size curve, thereby obtaining a higher-precision curve.

[0124] Step S4: Using the obtained full life cycle defect size as a label dataset and the root mean square value of the real bearing vibration signal as input, a bidirectional long short-term memory network is used as a training model to obtain a feature mapping network model capable of real-time feature mapping;

[0125] Step S5: Collect vibration data during the operation of the bearing and convert it into feature channels by sampling. Use the trained feature real-time mapping network model ((the feature mapping model is the bidirectional long short-term memory network trained in step 4));) as a feature supplement tool, use the mapped defect size as an additional feature to supplement the feature set data of the real bearing, and input the supplemented feature set data into the dual-correlation dynamic graph neural network to obtain the final remaining service life prediction result; the feature channel specifically includes time domain features and frequency domain features, and the features in the time domain include variance, standard deviation, maximum and minimum values, root mean square, skewness, kurtosis, peak-to-peak value, and entropy. The frequency domain features are extracted by fast Fourier transforming the original signal, including median frequency, average frequency, maximum amplitude, maximum power spectrum density, average power spectrum density, variance of power spectrum density, skewness, kurtosis, standard deviation frequency, asymmetry coefficient of frequency distribution, and kurtosis of frequency distribution.

[0126] The overall flow diagram of the present invention is shown in Figure 2 In step S1, the present invention constructs a bearing dynamics model process based on Hertz contact theory and is used to describe the expansion process of the bearing outer ring defect, which mainly includes the following steps:

[0127] Step S0101: Introduce Hertz contact theory to establish a dynamic model of the bearing, so as to more accurately reflect the expansion process of the local defects of the bearing. The established dynamic model is derived from the differential equation:

[0128]

[0129] Where M is the mass of the system, C represents the damping coefficient, K is the stiffness, and F represents the excitation load, which represent the acceleration, velocity, and displacement of the system respectively;

[0130] Step S0102: Based on Newton's second law, the two-degree-of-freedom dynamic equation of the rolling bearing considering damping and constant external load can be obtained, which can be expressed as the following equation:

[0131]

[0132] Where M is the equivalent mass of the bearing, C b is the equivalent damping of the bearing, C is the equivalent damping of the bearing, and is the vibration acceleration in the x and y directions, and is the vibration speed, F Cx and F Cy is the contact force component between the rolling element and the outer ring in the x and y directions, Q x and Q y is the load component of the bearing in the x and y directions;

[0133] Step S0103: The contact force component is regarded as a nonlinear contact force and is calculated using the Hertz contact theory. The following equation can be obtained using the Hertz contact theory:

[0134]

[0135] Where K represents the total contact stiffness of the bearing, N represents the number of rolling elements, μ is the factor that determines whether the roller is in contact with the outer ring, and θ i Indicates the angular position of rolling element i; for cylindrical roller bearings For ball bearings

[0136] Step S0104: Based on Newton's second law, the two-degree-of-freedom dynamic equation of the rolling bearing considering damping and constant external load can be obtained. Substitute the equation of step 1.3 into the equation of step 1.2 to calculate the final defect expansion model of the rolling bearing twin, which is expressed by the following equation:

[0137]

[0138] Where K represents the total contact stiffness of the bearing, N represents the number of rolling elements, and μ is a factor that determines whether the roller is in contact with the outer ring. μ satisfies the following expression:

[0139]

[0140] Among them, Δ i Indicates the contact deformation between rolling element i and raceway, Δ i Satisfies the following equation:

[0141] Δ i =xcosθ i +ysinθ i -cH;

[0142] Where H represents the time-varying displacement excitation, θ i represents the angular position of rolling element i at time t;

[0143] Step S0105: Calculate θ using the following expression: i , first convert the angular velocity of the shaft into radians, then calculate the speed of the bearing cage, and finally calculate the position of the rolling elements:

[0144]

[0145] Among them, n s represents the rotational speed of the shaft in the rotating system, ω s represents the angular velocity of the shaft, d represents the diameter of the rolling element, D pIndicates the pitch diameter of the bearing, ω c represents the speed of the bearing cage, and Z represents the number of rolling elements;

[0146] Combine Figure 3 The figure shows a schematic diagram of the two-degree-of-freedom twin bearing dynamics model provided by the present invention. i represents the angular position of the i-th rolling element, ω s represents the angular velocity of the shaft. To establish a real-time interactive digital twin framework, a two-degree-of-freedom bearing dynamics model based on Hertzian contact theory was built. In this example, M is set to 0.3 kg, C is set to 200 N·s / m, and n is set to 10 / 9 to obtain a dynamic response that is most similar to that of a real bearing.

[0147] Step S0106: Calculate the time-varying displacement excitation ΔH over the entire life cycle of the bearing in two stages. Regarding the time-varying displacement excitation H, the entire life cycle of the bearing is usually divided into two stages for calculation to simulate the displacement excitation generated by the evolution of a real bearing defect. To simplify the dynamic model, the contact deformation between the rolling element and the defect edge is ignored, and the defect is regarded as a rectangle with a length of L and a depth of H. An expression for the two-stage displacement excitation function is established. The time-varying displacement excitation is calculated to simulate the operation process of a real bearing. When a defect occurs in a real bearing, the displacement excitation at the defect changes constantly, so a mathematical model needs to be designed to represent this change process. The two stages are the initial generation stage and the defect expansion stage. In each stage, the maximum displacement excitation and the time-varying displacement excitation are calculated respectively.

[0148] In the present invention, the two-stage calculation of the time-varying displacement excitation of the bearing over its entire life cycle includes the following steps:

[0149] Step 3.1: First, simulate the displacement excitation generated by the evolution of the real bearing defect. In order to simplify the dynamic model and ignore the contact deformation between the bearing rolling element and the edge of the bearing outer ring defect, the bearing outer ring defect is regarded as a rectangle with a length of L and a depth of H, and a two-stage displacement excitation model is established. At the initial stage of the defect, due to the small length of the defect, the rolling element cannot completely contact the bottom of the defect, and its maximum displacement excitation H is max Always smaller than the defect depth H, H max It can be calculated using the following expression:

[0150] H max =0.5d-((0.5d) 2 -(0.5L) 2 )

[0151] At this time, the time-varying displacement excitation function can be expressed as the following expression:

[0152]

[0153] Among them, b i is the residual angular position of rolling element i at time t, b i =mod(θ i ,2π),θ d is the angle corresponding to the defect center, and δ represents the span angle corresponding to the defect;

[0154] Step 3.2: During the defect expansion phase, as the defect size increases, the maximum displacement excitation gradually increases until it reaches the defect depth H. At this point, the rolling element completely falls into the defect. This phenomenon can usually be observed in the middle and late stages of failure. At this time, the time-varying displacement excitation function is expressed as follows:

[0155]

[0156] Among them, mod(·) represents the remainder operation, b i is the residual angular position of rolling element i at time t,

[0157] b i =mod(θ i ,2π),θ d is the angle corresponding to the defect center, δ1 represents the span angle corresponding to the defect, and δ2 represents the span angle of the ball center when the rolling element contacts the bottom of the defect.

[0158] In the present invention, Figure 4 The figure shows a schematic diagram of the defect expansion process of the outer ring of the twin bearing of the present invention. It shows how to divide the entire life cycle of the bearing into two stages for calculation to simulate the displacement excitation generated by the evolution of real bearing defects. In order to simplify the dynamic model, the contact deformation between the rolling element and the edge of the defect is ignored, and the defect is regarded as a rectangle with a length of L and a depth of H. It shows the expansion process of the local defect of the bearing considered by the present invention. On this basis, the expression of the displacement excitation function of the two stages is established; wherein δ1 represents the span angle of the defect, and δ2 represents the span angle of the moving distance of the center of the rolling element in the defect. At the initial stage of defect generation, Figure 4 (a) It can be seen that at this time, due to the small length of the defect, the rolling element cannot completely contact the bottom of the defect. Figure 4 (b) As the defect size increases, the maximum displacement excitation gradually increases until it reaches the defect depth H, at which point the rolling element is completely engulfed in the defect. This phenomenon is typically observed in the middle and late stages of failure. To obtain a time-varying displacement excitation that closely reflects the actual bearing degradation process, the entire bearing lifecycle is typically divided into two stages for calculation. To simplify the dynamic model, the contact deformation between the rolling element and the defect edge is ignored, and the defect is considered a rectangle with a length of L and a depth of H. Based on this localized bearing defect expansion process, an expression for the two-stage displacement excitation function is established.

[0159] In step 2 of the present invention, the defect size is used as a variable, and multi-objective optimization is used to match the dynamic responses of the twin and the real bearing, achieving real-time interaction between the twin and the real bearing, and obtaining a preliminary full-life cycle defect size curve. The defect size is used as the variable x, and the multi-objective optimization continuously changes the value of x to match the optimal fitness. After matching the dynamic responses at all times, the defect size x at all times can be obtained, that is, the size curve is obtained, including the following steps:

[0160] Step 2.1: Use multi-objective optimization matching to iterate the parameters of the twin. The optimization target considers the similarity of time domain and frequency domain at the same time to ensure the best matching defect parameters. The multi-objective optimization matching satisfies the following expression:

[0161]

[0162] Among them, λ op represents the optimized defect parameters, including defect length L and defect depth H. T and P represent the dynamic responses of the twin and the real bearing, respectively; F(·) represents fast Fourier transform; dtw(·) represents DTW distance, and pearson(·) represents Pearson correlation coefficient;

[0163] Step 2.2: Use the particle swarm algorithm to iteratively update the position of the flock, guiding it to converge toward the optimal solution, ultimately finding the global optimal solution. The particle swarm algorithm is an algorithm that simulates the foraging behavior of birds. As a swarm intelligence evolutionary algorithm, its principle is to compare the location of food to the optimal solution to the optimization problem, replacing the movement direction and speed of particles with the flight direction and speed of birds. During the swarm update process, each bird will simultaneously consider its own memory of the food location and the influence of collective memory, thereby adjusting its flight direction. In this way, the position of the bird flock is continuously iteratively updated, guiding it to converge toward the optimal solution, and ultimately finding the global optimal solution. The specific update process of the particle swarm algorithm can be expressed as follows:

[0164]

[0165] in, represents the velocity vector of particle k after the tth iteration, represents the position vector of particle k after the tth iteration, τ represents the inertia coefficient, η is a random number in [0,1]; α1 represents the individual learning factor, α2 represents the social learning factor, represents the optimal position of particle i in t iterations, gbest t Indicates the optimal position in the group;

[0166] Step 2.3: After calculating the velocity of the t+1th iteration, the position of each particle is updated according to the following formula. The position of the particle represents the size of the defect:

[0167] After iteration, the velocity at the next moment is obtained, and then the particle updates its position again according to this formula; the defect size is used as the variable x, and the multi-objective optimization continuously changes the value of x to match the optimal fitness. After matching the dynamic responses at all moments, the defect size x at all moments can be obtained, that is, the size curve is obtained.

[0168] In step 3 of the present invention, the two-stage update of the defect size of the bearing throughout its entire life cycle includes the following steps:

[0169] Step S301: Use the bearing time domain vibration signal to describe the rolling element passing through a local defect. The physical process at a certain moment is expressed by the following equation:

[0170] (OF·cosα1-Rcosα2) 2 +(OF·sinα1-Rsinα2-R·ω c ·t p ) 2 =r b 2 ;

[0171] Among them, α1 represents the span angle between the rolling body center point A and point E when the rolling body just contacts the defect entrance E, α2 represents the span angle between the rolling body center point C and point E when the rolling body just leaves the defect entrance E, t p The time it takes for the rolling element to rotate from point C to just contact the defect outlet point F is the time it takes for the bearing to rotate, r b is the radius of the rolling element;

[0172] The present invention uses the most realistic defect size possible to readjust the obtained full life cycle defect curve, introduces a signal processing technology from the literature (Signal processing techniques for rolling element bearing spallsize estimation), and solves the more accurate bearing defect size in a certain time period by analyzing the bearing time domain vibration signal.

[0173] S0302: Extract the position of each point of the bearing time domain vibration signal to obtain the span angle α1 of the defect, and calculate the bearing calibration defect size U based on the span angle t , bearing calibration defect size U t Use the following expression for calculation:

[0174]

[0175] In the present invention, after obtaining the span angle, the bearing defect size can be calculated with high accuracy:

[0176] Step S303: Apply steps 301-302 to the vibration signal at the selected moment to obtain the calibration defect size at that moment, and apply the following expression to the full life cycle bearing defect curve generated in one stage to calculate

[0177] Calculate the final defect size of the two-stage update model:

[0178]

[0179] Among them, L t represents the first stage optimization defect corresponding to the calculated first calibration defect, n represents the total time of bearing failure, γ is the number of calibration defects obtained, U t represents the calibration defect size at the tth moment.

[0180] In the present embodiment, Figure 5 Figure 2 shows a schematic diagram depicting the physical process of a defect in the outer race of a twin bearing rolling element, as provided by the present invention. In this figure, when the rolling element reaches point A, it begins to enter the defect. Due to the reduced contact area, the resulting dynamic response gradually decreases until it reaches a local minimum, point B. Since point B is not required for defect size calculation, it is not shown in the figure. Subsequently, the dynamic response gradually increases, reaching its peak at point C, just as the rolling element is about to lose contact with the raceway. After point C, the rolling element passes through the defect and reaches point D, where it collides with the defect edge and generates a high-frequency response.

[0181] The following is further described with reference to specific embodiments:

[0182] Taking the XJTU-SY dataset as an example, eight time-domain features, 11 frequency-domain features, and eight time-frequency-domain features were extracted from the vertical direction of the original vibration signal, for a total of 28 features. In the XJTU-SY dataset, 15 bearings were selected under three operating conditions: Operating Condition 1 (2100 rpm / 12000 N), bearings 1-1 through 1-5; Operating Condition 2 (2250 rpm / 11000 N), bearings 2-1 through 2-5; and Operating Condition 3 (2400 rpm / 10000 N), bearings 3-1 through 3-5. The eight time-domain features in this example include variance, kurtosis, entropy, maximum and minimum values, peak-to-peak value, skewness, root mean square (RMS), and standard deviation. Frequency-domain features are extracted from the raw signal after fast Fourier transform (FFT). These features include median frequency, mean frequency, maximum amplitude, maximum power spectral density, mean power spectral density, power spectral density variance, skewness, kurtosis, standard deviation frequency, frequency distribution asymmetry coefficient, and frequency distribution kurtosis. Eight time-frequency domain features are obtained through three-layer wavelet packet decomposition using Daubechies wavelet basis functions. These features have been shown to be effective in predicting bearing RUL. A trained bidirectional long-short-term memory network is then used to map real-time defect size features and add them to the feature space. Two graph data sets with different adjacency matrices are constructed, and the DGCN module is used for spatiotemporal feature extraction. When the rolling element reaches point A, it begins to enter the defect. Due to the reduction in contact area, the dynamic response gradually decreases until it reaches a local minimum point B. Since point B is not required for defect size calculation, it is not shown in the figure. Subsequently, the dynamic response gradually increases and reaches its peak at point C, indicating that the rolling element is about to lose contact with the raceway. After point C, the rolling element passes through the defect and arrives at point D, where it collides with the edge of the defect and generates a high-frequency response. In order to obtain the high-fidelity bearing defect size at some moments, the calculation theory in the literature is introduced. It should be noted that the present invention only selects specific moments to calculate the calibration defect. There are two reasons for this: (1) The calculation process of the above-mentioned calibration defect is based on the double pulse phenomenon of the bearing. In the early stage of bearing failure, it is often difficult to capture this phenomenon. Therefore, this method is only accurate after the defect expands to a certain size. (2) The calculation cost of the calibration defect is relatively high. It is unrealistic and unreasonable to calculate the calibration defect of the bearing at every moment.

[0183] In step S4 of the present invention, the execution process is to use a bidirectional long short-term memory network to establish a connection between the virtual and real spaces. This connection is based on the superior learning ability of the neural network, can receive the real-time vibration signal of the bearing, and then use the trained neural network to output the real-time defect size of the twin, and quickly feed it back to the real bearing to guide its remaining service life prediction. Training the bidirectional long short-term memory network to obtain a network model that can perform real-time feature mapping includes the following steps:

[0184] Step S41: Using a bidirectional long short-term memory network to establish a connection between the virtual and real spaces, the root mean square value is calculated based on the real-time vibration signal of the bearing. The bidirectional long short-term memory network is a feature mapping model, and the root mean square feature is used to map the defect size feature.

[0185] Step S42: Obtain the full life cycle defect size of the bearing through the two-stage update model of step 3 and use it as a training data set; the two-stage update of the present invention is to use the calibration defect to update the defect curve obtained in the first stage again;

[0186] Step S43: Using the root mean square value of the bearing throughout its life cycle as input data and the defect size of the bearing throughout its life cycle as output data, a bidirectional long short-term memory network is trained, and a well-performing feature mapping model is obtained through verification;

[0187] Step S44: Input the real-time root mean square value of the bearing into the trained feature mapping network model, and then use the trained model to output the real-time defect size of the twin, and quickly feed it back to the real bearing to guide its remaining service life prediction;

[0188] In step 5 of the present invention, the supplemented feature set data is input into the dual-correlation dynamic graph neural network to obtain the final remaining service life prediction result. The execution process mainly obtains the twin features through the trained bidirectional long short-term memory network, fuses them with the features of the real bearing, and finally inputs them into the dual-correlation dynamic graph neural network to obtain the final prediction result, including the following steps:

[0189] Step S501: Extract the features of the vibration data during the operation of the bearing, wherein the features of the vibration data mainly include time domain features and frequency domain features. The features in the time domain include variance, standard deviation, maximum and minimum values, root mean square, skewness, kurtosis, peak-to-peak value, and entropy. The frequency domain features are extracted by subjecting the original signal to fast Fourier transform. The features in the frequency domain include median frequency, average frequency, maximum amplitude, maximum power spectrum density, average power spectrum density, variance of power spectrum density, skewness, kurtosis, standard deviation frequency, asymmetry coefficient of frequency distribution, and kurtosis of frequency distribution. The root mean square value of the entire life cycle of the bearing is used as the input feature and input into the trained bidirectional long short-term memory network to obtain twin features. The twin features are then fused with the features of the real bearing to obtain the fused feature channel data.

[0190] Step S502: using a moving window with a stride of 1 and a length of L to sample the feature set that integrates the twin defect features;

[0191] Step S503: Using the columns of the window data as graph node vectors, a spatial correlation graph is established; using the rows of the window data as graph node vectors, a temporal correlation graph is established;

[0192] Step S504: Use cosine similarity to define the similarity matrix between graph node vectors. First, build an adjacency matrix based on the feature set data, and then update the feature set data using the built adjacency matrix. The cosine similarity matrix C (i,n) Elements in It is calculated by the following expression:

[0193]

[0194] in represents the f-th row feature of the n-th sampling data of the i-th bearing, ||·|| is the second norm of the vector;

[0195] S0505: Apply function o(·) to the cosine similarity matrix C (i,n) The elements in can get the adjacency matrix Elements in δ is a threshold defined artificially, and the function o(·) is expressed as:

[0196]

[0197] S0506: Define the normalized graph Laplacian matrix, expressed as follows:

[0198]

[0199] Where D represents the degree matrix, and the expression is IN represents the identity matrix;

[0200] By building two different types of time correlation graphs and then merging them, the purpose of building the graph is to use graph neural networks to extract correlation information;

[0201] S0506: In order to obtain the correlation information of the feature set data through the adjacency matrix, the information transfer process of the graph convolutional neural network model is defined. This process can update the feature set data, which is expressed as follows:

[0202]

[0203] in, Is an adjacency matrix with self-loop, that is, a self-connection is added to the original adjacency matrix, W (L) represents the trainable weight matrix of the Lth layer, σ(·) represents the ReLU activation function, and H (L) Represents the input of the Lth layer. When L=0, H (0) =0;

[0204] Step S508: Input the fused feature channel data into the dual-correlation dynamic graph neural network to extract spatiotemporal correlation. The dual-correlation dynamic graph neural network update mechanism is as follows:

[0205]

[0206] Among them, σ represents the Sigmiod activation function, F(·) represents the graph convolution process (i.e., the convolution process of the graph neural network model GCN), W and b represent the trainable weight matrix and bias, i t 、f t 、o t are the input gate, forget gate, and output gate at the tth position, respectively, c t is the cell state at the tth position, h t is the hidden layer information at the tth position;

[0207] Step S508: After the spatiotemporal correlation is extracted using the dual-correlation dynamic graph neural network, the updated feature set data is flattened and input into the fully connected layer to obtain the final remaining useful life prediction result, which is expressed as follows: t =σ(W t ·Flatten(g t ));

[0208] Among them, y t Represents the prediction result, σ represents the ReLU activation function, W t Represents the trainable weight matrix, g t Indicates the updated feature set data.

[0209] In the present invention, Figure 6 As shown, Figure 6 This is a schematic diagram of the DGCN module provided by the present invention. Graph convolutional neural networks can capture the spatial correlation between features well by aggregating neighboring nodes. However, using only graph convolutional neural networks cannot obtain the temporal dependency between features, that is, the temporal correlation. Therefore, this paper designs a DGCN module to simultaneously obtain the spatiotemporal correlation of fused feature channel data and use it for RUL prediction. This module can retain valid historical information while capturing the current spatiotemporal domain correlation, so its hidden layer simultaneously contains the temporal and spatial correlation information of the fused feature channel. In this example, the dimension of the hidden layer of the DGCN module is taken as 6 to obtain the best overall performance. Among them, Figure 7 This is a rendering of the dynamic response of the twin model provided by the present invention. It can be seen that the dynamic response of the twin exhibits significant periodicity due to the presence of a localized defect. Since the location of the outer ring defect remains constant, the impact signals generated are uniformly spaced. The impact period caused by the rolling element passing through the defect is the inverse of the characteristic frequency period of the outer ring fault.

[0210] Figure 8 This is a schematic diagram comparing the envelope spectrum analysis of the twin bearing provided by the present invention and the real bearing. It can be observed that the envelope signals of the twin signal and the real signal are basically similar. Due to the presence of noise in the operating environment of the real bearing, its dynamic response will show irregularity, and the amplitude is slightly larger than the twin signal, while the twin signal is more regular. The envelope spectrum of the simulation data can clearly show the fault characteristic frequency of the bearing and its multiple frequencies. When compared with the envelope spectrum of the real bearing, the difference is very small. Therefore, it can be considered that this twin model can represent the operating status of the real bearing to a certain extent.

[0211] Figure 9 This is a first-stage bearing defect size curve obtained based on multi-objective optimization, as provided by the present invention. It shows the first-stage bearing defect curve obtained after multi-objective optimization. The different defect sizes corresponding to each minute in the figure are the Pareto optimal solutions solved by the multi-objective optimization at that time.

[0212] Figure 10 This is a first-stage bearing defect size curve after averaging provided by the present invention. The defect sizes under Pareto dominance obtained every minute are averaged to obtain a first-stage bearing life cycle defect size curve based on multi-objective optimization.

[0213] Figure 11This is a graph of the bearing defect size curve after a two-stage update using calibration defects, as provided by the present invention. It can be seen that the bearing defect curve derived in the first stage can relatively accurately track the overall trend. However, in the later stages of bearing failure, the defect size curve obtained by the multi-objective optimization algorithm differs from the curve after the two-stage update. This is because during periods of larger bearing defects, the dynamic response is more unstable, leading to errors in defect matching. The two-stage calibration defects, however, can guide the further update of the first-stage defect curve, resulting in a high-fidelity bearing defect curve for the entire life cycle.

[0214] Figure 12 This figure shows a comparison of the lifecycle defect size curves obtained by this invention and other methods reported in the literature. The two-stage updated lifecycle defect size curves show similar trends to those in other relevant literature, validating the effectiveness of the proposed DT model. Because this paper employs multi-objective optimization to fit the single-stage curves, the defect size curves are more stable, facilitating subsequent RUL prediction.

[0215] Figure 13 This is a schematic diagram of the effect of feature mapping using a bidirectional long short-term memory network, as provided by the present invention. It can be seen that the mapping curve differs little from the defect size curve obtained by the two-stage update, verifying the effect of the bidirectional long short-term memory network trained in step 4 of the present invention and proving the effectiveness of the training. The effectiveness of the mapping model has been verified, and it can effectively map the dynamic response of the real bearing to the defect size of the twin. This real-time mapping reveals the relationship between the bearing vibration signal and the operating defect size, and it has excellent immediacy and accuracy. This article will subsequently use this trained model to obtain the real-time operating status of the bearing and incorporate it into the bearing's RUL prediction.

[0216] Figure 14 and Figure 15 The following are comparisons of the score indicators of the various models provided by the present invention for the XJTU-SY dataset in the S2 and S1 scenarios. It can be seen that the proposed model has the lowest score, which indicates that it has the highest prediction accuracy and its predictions are more conducive to guiding practical work. Although the convolutional gated recurrent unit has a good root mean square error, it performs poorly in terms of the score indicator, which means that its prediction curve cannot be well applied, thus verifying the superiority of the invention proposed in the remaining useful life prediction.

[0217] Figure 16 and Figure 17Figure 2 shows a comparison of the prediction performance of various models for the XJTU-SY dataset under the S2 and S1 scenarios, respectively. It can be seen that most models accurately reflect bearing degradation in their overall trends, but TCN (Transistor-Converted Genetic Network) does not achieve satisfactory results in terms of accuracy. Compared to the two models that extract temporal correlations separately, the model that simultaneously extracts spatiotemporal correlations demonstrates better accuracy. Furthermore, after adding the full-life defect size as a feature, the model's prediction curves more closely approximate the actual life curves. In both scenarios, DC-DGCN tracks bearing degradation trends better than other models and demonstrates superior RUL prediction capabilities. A comparison of the two scenarios shows that the inclusion of the real-time defect dimension significantly improves the DC-DGCN's prediction curves. The results demonstrate that the proposed method effectively captures the general degradation characteristics of bearings by incorporating defect size features, improving the accuracy of bearing RUL prediction, especially in the later stages of bearing life. This is crucial for RUL prediction, as the initial RUL typically requires less attention. In contrast, accurate late-stage RUL predictions help improve equipment reliability and reduce wear caused by unexpected failures. Furthermore, it's noteworthy that incorporating defect size features into the proposed method's predictions did not significantly improve the prediction results for some bearings. This phenomenon is attributed to the lack of visible outer ring defects during operation. The twin model constructed in this paper only describes the defect growth process in the outer ring.

[0218] Figure 18 Figure 3 compares the prediction performance of the dual-correlation dynamic graph neural network proposed in this paper for a portion of the XJTU-SY dataset under two different scenarios. The model in scenario S1 demonstrates higher prediction accuracy, with smaller errors in the later stages of the prediction compared to scenario S2. This advantage arises because the DC-DGCN model effectively utilizes defect sizes across the entire lifecycle, and the correlations captured by the model significantly enhance RUL predictions.

[0219] The above disclosure is only a preferred embodiment of the present invention, and certainly cannot be used to limit the scope of the rights of the present invention. Ordinary technicians in this field can understand that all or part of the processes of the above embodiment and equivalent changes made in accordance with the claims of the present invention are still within the scope of the invention.

Claims

1. A method for predicting the remaining service life of a bearing based on a dynamic graph neural network, characterized in that: The following steps are involved: Step S1: Establish a bearing dynamics model based on Hertz contact theory to describe the expansion process of the bearing outer ring defect and calculate the contact deformation; Step S2: Using the defect size as a variable, multi-objective optimization is used to match the dynamic response of the bearing digital twin and the real bearing to obtain a preliminary full life cycle defect size curve; Step S3: extracting high-fidelity defect size from the bearing time-domain vibration signal to calculate the calibrated defect size, and then perform a two-stage update of the lifecycle defect size; Step S4: Using the obtained full life cycle defect size as a label dataset and the root mean square value of the real bearing vibration signal as input, a bidirectional long short-term memory network is used as a training model to obtain a feature mapping network model capable of real-time feature mapping; Step S5: Collect vibration data during the operation of the bearing and convert its samples into feature channels. Use the trained feature real-time mapping network model as a feature supplement tool, use the mapped defect size as an additional feature to supplement the feature set data of the real bearing, and input the supplemented feature set data into the dual-correlation dynamic graph neural network to obtain the final remaining service life prediction result.

2. A method for predicting the remaining useful life of a bearing based on a dynamic graph neural network according to claim 1, characterized in that: In step 1, establishing a bearing dynamics model based on Hertz contact theory to describe the expansion process of the bearing outer ring defect includes the following steps: Step S0101: Introduce Hertz contact theory to establish a dynamic model of the bearing. The established dynamic model satisfies: Where M is the mass of the system, C is the damping coefficient, K is the stiffness, and F is the excitation load. x represents the acceleration, velocity and displacement of the system respectively; Step S0102: Newton's second law is used to derive a two-degree-of-freedom dynamic equation for a rolling bearing that takes into account damping and a constant external load, which is expressed as the following equation: Among them, M b is the equivalent mass of the bearing, C b is the equivalent damping of the bearing, and is the vibration acceleration in the x and y directions, and is the vibration speed, F Cx and F Cy is the contact force component between the rolling element and the outer ring in the x and y directions, Q x and Q y is the bearing load component; Step S0103: All contact force components are considered as nonlinear contact forces, and the following equation is calculated using Hertz contact theory: Where K represents the total contact stiffness of the bearing, N represents the number of rolling elements, μ is the factor that determines whether the roller is in contact with the outer ring, and θ i Indicates the angular position of rolling element i. For cylindrical roller bearings For ball bearings Step S0104: Substitute the equation in step 1.3 into the equation in step 1.2 to calculate and obtain the final defect expansion model of the rolling bearing twin, which is expressed by the following equation: Where K represents the total contact stiffness of the bearing, N represents the number of rolling elements, and μ is a factor that determines whether the roller is in contact with the outer ring. μ satisfies the following expression: Among them, Δ i Indicates the contact deformation between rolling element i and raceway, Δ i Satisfies the following equation: D i =xcosθ i +ysinθ i -cH; Where H represents the time-varying displacement excitation, θ i represents the angular position of rolling element i at time t; Step S0105: Calculate the angle position θ i First, convert the angular velocity of the shaft into radians, then calculate the speed of the bearing cage, and finally calculate the position of the rolling element i. The specific calculation process is as follows: Among them, n s Indicates the shaft speed, ω s represents the angular velocity of the shaft, d represents the diameter of the rolling element, D p Indicates the pitch diameter of the bearing, ω c represents the speed of the bearing cage, and Z represents the number of rolling elements; Step S0106: Calculate the time-varying displacement excitation ΔH of the bearing over its entire life cycle in two stages.

3. The method for predicting the remaining useful life of a bearing based on a dynamic graph neural network according to claim 2, wherein: The two-stage calculation of the time-varying displacement excitation ΔH over the entire life cycle of a bearing includes the following steps: Step 3.1: First, simulate the displacement excitation generated by the evolution of a real bearing defect. To simplify the dynamic model and ignore the contact deformation between the bearing rolling elements and the edge of the bearing outer ring defect, the bearing outer ring defect is considered as a rectangle with a length of L and a depth of H, and a two-stage displacement excitation model is established. At the initial stage of defect generation, the rolling element cannot completely contact the bottom of the defect, and its maximum displacement excitation H max Always smaller than the defect depth H, H max Calculated using the following expression: <h2 style=";text-align:left;direction:ltr">H<h2 style=";text-align:left;direction:ltr"> max <h2 style=";text-align:left;direction:ltr"> =0.5d-((0.5d)<h2 style=";text-align:left;direction:ltr"> 2 <h2 style=";text-align:left;direction:ltr"> -(0.5L)<h2 style=";text-align:left;direction:ltr"> 2 <h2 style=";text-align:left;direction:ltr"> ); At this time, the time-varying displacement excitation function satisfies the following expression: Among them, b i is the residual angular position of rolling element i at time t, b i =mod(θ i ,2π),θ d is the angle corresponding to the defect center, and δ1 represents the span angle corresponding to the defect; Step 3.2: In the defect expansion stage, as the defect size increases, the maximum displacement excitation gradually increases until it is equal to the defect depth H. At this time, the rolling element will completely fall into the defect. At this time, the time-varying displacement excitation function in the defect expansion stage satisfies the following expression: Among them, mod(·) represents the remainder operation, b i is the residual angular position of rolling element i at time t, b i =mod(θ i ,2π),θ d is the angle corresponding to the defect center, δ1 represents the span angle corresponding to the defect, and δ2 represents the span angle of the ball center when the rolling element contacts the bottom of the defect.

4. The method for predicting the remaining useful life of a bearing based on a dynamic graph neural network according to claim 1, wherein: In step 2, the defect size is used as a variable, and multi-objective optimization is used to match the dynamic response of the twin and the real bearing to obtain a preliminary full life cycle defect size curve. The steps include the following: Step 2.1: Use multi-objective optimization matching to iterate the parameters of the twin. The optimization target considers the similarity of time domain and frequency domain at the same time to ensure the best matching defect parameters. The multi-objective optimization matching satisfies the following expression: Among them, λ op represents the optimized defect parameters, including defect length L and defect depth H; T and P represent the dynamic responses of the twin and the real bearing, respectively; F(·) represents fast Fourier transform; dtw(·) represents DTW distance; pearson(·) represents Pearson correlation coefficient; Step 2.2: Use the particle swarm algorithm to iteratively update the position of the group, guide it to converge towards the optimal solution, and finally find the global optimal solution. The specific update process is expressed as follows: in, represents the velocity vector of particle k after the tth iteration, represents the position vector of particle k after the tth iteration, τ represents the inertia coefficient, η is a random number in [0,1]; α1 represents the individual learning factor, α2 represents the social learning factor, represents the optimal position of particle i in t iterations, gbest t Indicates the optimal position in the group; Step 2.3: After calculating the velocity of the t+1th iteration, the position of each particle is updated according to the following formula. The position of the particle represents the size of the defect:

5. The method for predicting the remaining useful life of a bearing based on a dynamic graph neural network according to claim 1, wherein: In step S3, the two-stage update of the defect size of the bearing throughout its life cycle includes the following steps: Step S301: Use the bearing time-domain vibration signal to describe the physical process of the rolling element passing through a local defect at a certain moment, which is expressed by the following equation: (OF·cosα1-Rcosα2) 2 +(OF·sinα1-Rsinα2-R·ω c ·t p ) 2 =r b 2 ; Among them, α1 represents the span angle between the rolling body center point A and point E when the rolling body just contacts the defect entrance E, α2 represents the span angle between the rolling body center point C and point E when the rolling body just leaves the defect entrance E, t p The time it takes for the rolling element to rotate from point C to the moment it contacts the defect outlet point F is the time it takes for the bearing to rotate, r b is the radius of the rolling element; Step S302: Extract the position of each point of the bearing time domain vibration signal to obtain the span angle α1 of the defect, and calculate the bearing calibration defect size U according to the span angle t , bearing calibration defect size U t Use the following expression for calculation: Step S303: Apply steps S301-S302 to the vibration signal at the selected moment to obtain the calibrated defect size at that moment. Apply the following expression to the full life cycle bearing defect curve generated in the first stage to calculate the final defect size of the two-stage updated model: Among them, L t represents the first stage optimization defect corresponding to the calculated first calibration defect, n represents the total time of bearing failure, γ is the number of calibration defects obtained, U t represents the calibration defect size at the tth moment.

6. The method for predicting the remaining useful life of a bearing based on a dynamic graph neural network according to claim 4, wherein: In step 4, training a bidirectional long short-term memory network to obtain a network model capable of performing real-time feature mapping includes the following steps: Step S41, using a bidirectional long short-term memory network to establish a connection between the virtual and real spaces, and calculating the root mean square value of the bearing's real-time vibration signal; Step S42: Obtain the defect size of the bearing throughout its life cycle and use it as a training data set. Step S43: Using the root mean square value of the bearing throughout its life cycle as input data and the defect size of the bearing throughout its life cycle as output data, a bidirectional long short-term memory network is trained, and a well-performing feature mapping model is obtained through verification; Step S44: Input the real-time root mean square value of the bearing into the trained feature mapping network model, and then use the trained model to output the real-time defect size of the twin, and quickly feed it back to the real bearing to guide its remaining service life prediction.

7. The method for predicting the remaining useful life of a bearing based on a dynamic graph neural network according to claim 5, wherein: In step 5, the supplemented feature set is input into the dual-correlation dynamic graph neural network to obtain the final remaining useful life prediction result, which includes the following steps: Step S501: Extract the characteristics of the vibration data during the operation of the bearing, use the root mean square value of the entire life cycle of the bearing as the input feature, and input it into the trained bidirectional long short-term memory network to obtain the twin feature, and then fuse the twin feature with the characteristics of the real bearing; Step S502: using a moving window with a stride of 1 and a length of L to sample the feature set data that concentrates on the defect features of the bearing digital twin; Step S503: Using the columns of the window data as graph node vectors, a spatial correlation graph is established; using the rows of the window data as graph node vectors, a temporal correlation graph is established; Step S504: Use cosine similarity to define the similarity matrix between graph node vectors. The cosine similarity matrix C (i ,n) Elements in It is calculated by the following expression: in represents the f-th row feature of the n-th sampling data of the i-th bearing, ||·|| is the second norm of the vector; Step S505: Apply function o(·) to the cosine similarity matrix C (i,n) The elements in can be used to get the adjacency matrix Elements in δ is a threshold defined artificially, and the function o(·) is expressed as: Then, the normalized graph Laplacian matrix is defined by the obtained adjacency matrix A, which is expressed as follows: Where D represents the degree matrix, and the expression is I N represents the identity matrix; Step S506: In order to obtain the correlation information of the feature set data through the adjacency matrix, the information transfer process of the graph convolutional neural network model is defined. This process can update the feature set data, which is expressed as follows: in, Is an adjacency matrix with self-loop, that is, a self-connection is added to the original adjacency matrix, W (L) represents the trainable weight matrix of the Lth layer, σ(·) represents the ReLU activation function, and H (L) Represents the input of the Lth layer. When L=0, H (0) =0; Step S507: Input the fused feature channel data into the dual-correlation dynamic graph neural network to extract spatiotemporal correlation. The dual-correlation dynamic graph neural network update mechanism is as follows: Where σ represents the Sigmiod activation function, F(·) represents the graph convolution process, W and b represent the trainable weight matrix and bias, i t 、f t 、o t are the input gate, forget gate, and output gate at the tth position, respectively, c t is the cell state at the tth position, h t is the hidden layer information at the tth position; Step S508: After the spatiotemporal correlation is extracted using the dual-correlation dynamic graph neural network, the updated feature set data is flattened and input into the fully connected layer to obtain the final remaining useful life prediction result, which is expressed as follows: t =σ(W t ·Flatten(g t )); Among them, y t Represents the prediction result, σ represents the ReLU activation function, W t Represents the trainable weight matrix, g t Indicates the updated feature set data.

Images