A method and system for monitoring the operation of a machine tool using a three-axis gyroscope

Abstract

The application provides a method and system for monitoring the operation of a machining device by using a three-axis gyroscope, and relates to the technical field of mechanical device state monitoring. The method comprises the following steps: reconstructing phase space by purifying three-axis angular velocity time series data to obtain a three-dimensional phase space trajectory; determining a first phase pole, a second phase pole and a third reference node in the three-dimensional phase space trajectory; extracting a middle axis skeleton from the inside of the three-dimensional phase space trajectory by distance transformation, and counting the stem length, branch number and branch point curvature of the middle axis skeleton as the geometric invariant characteristic of the three-dimensional phase space trajectory form. The application can identify early hidden faults of the device, improve the monitoring sensitivity, and realize safe and stable operation and intelligent adaptive control of the device through phase space reconstruction, multi-dimensional feature fusion extraction, health quantification evaluation, hierarchical decision and closed-loop adaptive regulation.

Description

Technical Field

[0001] This invention relates to the field of mechanical equipment condition monitoring technology, and in particular to a method and system for monitoring the operation of machining equipment using a three-axis gyroscope. Background Technology

[0002] In the field of monitoring the operating status of machining equipment, existing technologies include solutions for data acquisition using single-axis or triaxial vibration sensors. However, some conventional methods mainly focus on the vibration amplitude of a single axis or simple time-domain statistics (such as root mean square, peak value, and kurtosis). They often lack effective means of extracting and analyzing the complex motion patterns of equipment in three-dimensional space, especially the coupling relationships between multi-axis motions and the internal geometric features of the motion trajectory. This may lead to existing monitoring methods being unable to sensitively identify the deep features characterizing the essential changes in motion from the raw data when early anomalies occur in the equipment.

[0003] For example, a five-axis machining center used for machining turbine disks in aero-engines requires its spindle to feed at high speed along a spatial curve when finishing the tenon and groove structure of the turbine disk, with the coordinated operation of three linear axes (X, Y, and Z) and two rotary axes. After long-term continuous operation, the Z-axis lead screw and nut pair may experience slight uneven wear, causing a small amount of periodic wobble in the spindle during vertical feed. At this time, a three-axis gyroscope mounted on the spindle will collect small angular velocity fluctuations in the X and Y axes that are related to the phase of the Z-axis motion. In existing typical monitoring methods, the standard deviation or peak value of the angular velocity in each axis is calculated separately and compared with a preset fixed threshold. Because the wobble amplitude is small, the individual statistics of each axis may still be within the normal range, making it difficult to trigger an alarm. However, this tiny change in phase coupling actually reflects the early performance degradation of the Z-axis drivetrain. If it is ignored, it may cause the positional accuracy of the tenon and groove to gradually exceed the tolerance, eventually resulting in unqualified blade assembly clearance.

[0004] In addition, although some existing systems attempt to fuse and analyze multi-axis data, they usually only calculate the mean or variance of the synthesized vector magnitude, and lack in-depth feature extraction methods for the morphological evolution of motion trajectories in phase space (such as the topological structure of the trajectory skeleton and the coupling degree of phase shifts in each axis). Summary of the Invention

[0005] This invention provides a method and system for monitoring the operation of machining equipment using a three-axis gyroscope, which can effectively identify early latent faults such as wear of transmission mechanisms, abnormal clearances, and guide rail runout.

[0006] To solve the above-mentioned technical problems, the technical solution of the present invention is as follows: In a first aspect, a method for monitoring the operation of machining equipment using a three-axis gyroscope, the method comprising: Step 1: Collect raw triaxial angular velocity time series data and preprocess it to obtain purified triaxial angular velocity time series data; Step 2: Reconstruct the phase space of the three-axis angular velocity time series data to obtain a three-dimensional phase space trajectory. Within the three-dimensional phase space trajectory, determine the first phase pole, the second phase pole, and the third reference node. Extract the central axis skeleton from the interior of the three-dimensional phase space trajectory through distance transformation, and statistically analyze the trunk length, number of branches, and curvature of the branch points of the central axis skeleton as geometric invariant features of the three-dimensional phase space trajectory morphology. Determine the phase offset modulus components and phase deflection angles of each axis based on the first phase pole, the second phase pole, and the third reference node. Calculate the cross-correlation matrix between the phase offset modulus components of different axes to obtain the multidimensional motion coupling index. Combine the phase offset modulus components of each axis, the phase deflection angle, the multidimensional motion coupling index, and the geometric invariant features to obtain a feature index vector. Step 3: Based on the comparison between the feature index vector and the pre-built feature benchmark library, obtain the real-time health score of the device, and monitor the deviation direction and magnitude of each component in the feature index vector relative to the pre-built feature benchmark library to determine the fault diagnosis prompt set. Step 4: Based on the real-time health score and fault diagnosis prompts of the equipment, a hierarchical decision is made to obtain a set of control commands; Step 5: Execute the set of control commands and verify the effect to form a closed-loop adaptive monitoring process.

[0007] Secondly, a system for monitoring the operation of machining equipment using a three-axis gyroscope includes: The preprocessing module is used to collect raw triaxial angular velocity time series data and perform preprocessing to obtain purified triaxial angular velocity time series data; The calculation module is used to reconstruct the phase space of the purification triaxial angular velocity time series data to obtain a three-dimensional phase space trajectory. Within the three-dimensional phase space trajectory, the first phase pole, the second phase pole, and the third reference node are determined. The central axis skeleton is extracted from the interior of the three-dimensional phase space trajectory through distance transformation, and the main trunk length, number of branches, and curvature of the branch points of the central axis skeleton are statistically analyzed as geometric invariant features of the three-dimensional phase space trajectory morphology. Based on the first phase pole, the second phase pole, and the third reference node, the phase offset modulus components and phase deflection angles of each axis are determined. The cross-correlation matrix between the phase offset modulus components of different axes is calculated to obtain the multidimensional motion coupling index. The phase offset modulus components of each axis, the phase deflection angle, the multidimensional motion coupling index, and the geometric invariant features are combined to obtain a feature index vector. The evaluation module is used to obtain the real-time health score of the equipment based on the comparison between the feature index vector and the pre-built feature benchmark library, and to monitor the deviation direction and magnitude of each component in the feature index vector relative to the pre-built feature benchmark library, and determine the set of fault diagnosis prompts. The decision-making module is used to make hierarchical decisions based on the real-time health score and fault diagnosis prompts of the equipment, and to obtain a set of control commands. The monitoring module is used to execute a set of control commands and verify their effects, forming a closed-loop adaptive monitoring process.

[0008] The above-described solution of the present invention has at least the following beneficial effects: This invention overcomes the limitations of relying solely on single statistical quantities in the time and frequency domains by reconstructing phase space and analyzing three-dimensional phase space trajectories. It can comprehensively characterize the complex motion patterns and internal coupling relationships of equipment in three-dimensional space, improving the ability to identify early, weak, and coupling anomalies. Based on the central axis skeleton, geometric invariant features such as trunk length, number of branches, and curvature of branch points are extracted. These features are unaffected by data amplitude fluctuations, temporal offsets, and operating condition disturbances, resulting in higher feature stability and discriminative power, ensuring the consistency and accuracy of equipment condition assessment.

[0009] By using phase poles and reference nodes, the phase offset modulus and phase deflection angle can be accurately characterized, and the multidimensional motion coupling index can be quantified to effectively identify early hidden faults such as wear of transmission mechanism, abnormal clearance, and guide rail runout.

[0010] Quantitative scoring and deviation analysis are performed using feature index vectors and a feature benchmark library. Compared to fixed threshold judgment, this method offers stronger anti-interference capabilities and more comprehensive diagnostic evidence, enabling the identification of fault types and anomaly severity. Based on health scores and fault indications, a tiered control strategy is implemented, achieving a gradient response from parameter maintenance and fine-tuning optimization to intervention and protection, ensuring processing continuity while improving equipment operational safety. Real-time data acquisition and effect verification after control command execution support strategy iteration and statistical parameter updates, forming a self-optimizing and self-correcting intelligent monitoring and control closed loop. Long-term operation continuously improves monitoring accuracy and decision reliability. Non-invasive, multi-dimensional motion monitoring is achieved using a three-axis gyroscope, compatible with mainstream CNC systems and industrial communication protocols. Attached Figure Description

[0011] Figure 1 This is a flowchart illustrating a method for monitoring the operation of machining equipment using a three-axis gyroscope, as provided in an embodiment of the present invention.

[0012] Figure 2 This is a schematic diagram of a system for monitoring the operation of machining equipment using a three-axis gyroscope, provided by an embodiment of the present invention. Detailed Implementation

[0013] Exemplary embodiments of the present disclosure will now be described in more detail with reference to the accompanying drawings. While exemplary embodiments of the present disclosure are shown in the drawings, it should be understood that the present disclosure may be implemented in various forms and should not be limited to the embodiments set forth herein. Rather, these embodiments are provided so that this disclosure will be thorough and complete, and will fully convey the scope of the disclosure to those skilled in the art.

[0014] like Figure 1 As shown, an embodiment of the present invention proposes a method for monitoring the operation of machining equipment using a three-axis gyroscope, the method comprising the following steps: Step 1: Collect raw triaxial angular velocity time series data and preprocess it to obtain purified triaxial angular velocity time series data; Step 2: Reconstruct the phase space of the three-axis angular velocity time series data to obtain a three-dimensional phase space trajectory. Within the three-dimensional phase space trajectory, determine the first phase pole, the second phase pole, and the third reference node. Extract the central axis skeleton from the interior of the three-dimensional phase space trajectory through distance transformation, and statistically analyze the trunk length, number of branches, and curvature of the branch points of the central axis skeleton as geometric invariant features of the three-dimensional phase space trajectory morphology. Determine the phase offset modulus components and phase deflection angles of each axis based on the first phase pole, the second phase pole, and the third reference node. Calculate the cross-correlation matrix between the phase offset modulus components of different axes to obtain the multidimensional motion coupling index. Combine the phase offset modulus components of each axis, the phase deflection angle, the multidimensional motion coupling index, and the geometric invariant features to obtain a feature index vector. Step 3: Based on the comparison between the feature index vector and the pre-built feature benchmark library, obtain the real-time health score of the device, and monitor the deviation direction and magnitude of each component in the feature index vector relative to the pre-built feature benchmark library to determine the fault diagnosis prompt set. Step 4: Based on the real-time health score and fault diagnosis prompts of the equipment, a hierarchical decision is made to obtain a set of control commands; Step 5: Execute the set of control commands and verify the effect to form a closed-loop adaptive monitoring process.

[0015] In this embodiment of the invention, phase space reconstruction and three-dimensional phase space trajectory analysis overcome the limitations of relying solely on single statistics in the time and frequency domains, improving the ability to identify early, weak, and coupled anomalies. Based on the central axis skeleton, geometrically invariant features such as trunk length, number of branches, and branch point curvature are extracted, unaffected by data amplitude fluctuations, temporal offsets, and operational disturbances, ensuring the consistency and accuracy of equipment status assessment. Precise characterization of phase offset modulus and phase deflection angle is achieved through phase poles and reference nodes, and a multi-dimensional motion coupling index is quantified, effectively identifying early latent faults such as transmission mechanism wear, clearance anomalies, and guide rail runout. Quantitative scoring and deviation analysis are performed using feature index vectors and a feature benchmark library, resulting in stronger anti-interference capabilities and more comprehensive diagnostic evidence, enabling the location of fault types and anomaly degrees. A graded control strategy is executed based on health scores and fault prompts, achieving a gradient response from parameter maintenance and fine-tuning optimization to intervention and protection, improving equipment operational safety while ensuring processing continuity. The self-optimizing and self-correcting intelligent monitoring and control closed loop continuously improves monitoring accuracy and decision reliability over long-term operation.

[0016] In a preferred embodiment of the present invention, step 1 includes: Step 100: Rigidly fix the three-axis gyroscope sensor to the reciprocating motion component of the machining equipment, and synchronously collect the instantaneous angular velocity data of the reciprocating motion component in three-dimensional space at a preset sampling frequency to form the original three-axis angular velocity time series data, specifically including: The triaxial gyroscope sensor is rigidly mounted on the reciprocating motion component of the machining equipment using a rigid connection structure. During installation, it is ensured that there is no relative sliding, loosening, or gap between the sensor body and the reciprocating motion component. The sensor's measurement coordinate system is aligned with the standard motion coordinate system of the machining equipment. Specifically, the sensor's X-axis is aligned with the horizontal motion axis of the machining equipment (e.g., the transverse feed direction of the worktable), the sensor's Y-axis is aligned with another horizontal motion axis of the machining equipment (e.g., the longitudinal feed direction of the worktable), and the sensor's Z-axis is aligned with the vertical motion axis of the machining equipment perpendicular to the worktable plane (e.g., the spindle lifting direction). This ensures that the sensor can realistically and without distortion follow the reciprocating motion component synchronously. A fixed sampling frequency of 200Hz is used to continuously and synchronously collect the instantaneous angular velocities generated by the reciprocating motion component during its three-dimensional motion. The X-axis, Y-axis, and Z-axis angular velocity data collected at each moment are arranged sequentially according to the time sequence of data acquisition, forming complete, continuous, and time-aligned raw triaxial angular velocity time-series data.

[0017] Step 101 involves preprocessing the raw triaxial angular velocity time series data to obtain purified triaxial angular velocity time series data. Specifically, this includes performing a complete data preprocessing operation on the raw triaxial angular velocity time series data obtained in Step 100, sequentially removing outliers and mutations, and filtering out high-frequency electrical noise and environmental interference signals. First, a fixed sliding time window with a length of 20 sampling points is used to traverse the raw triaxial angular velocity time series data point by point. Within each sliding window, the arithmetic mean and three times the standard deviation of all data within the current window are calculated to construct a normal value range centered on the window mean and with three times the standard deviation as the interval width. If the actual value of a data point is greater than the sum of the window mean and three times the standard deviation, or less than the difference between the window mean and three times the standard deviation, then the data point is identified as an outlier or mutation and is directly removed from the time series data to ensure the continuity and effectiveness of the time series data. Signals with frequency components greater than 50Hz are identified as high-frequency electrical noise and environmental interference signals. A low-pass digital filter with linear phase characteristics is selected to filter the time series data after outlier removal. The filter cutoff frequency is set to 50Hz. During the filtering process, effective motion signals at or below 50Hz are retained, while interference signal components above 50Hz are completely filtered out. This makes the time series data smoother, more stable, and more in line with the actual motion state of the equipment. Finally, the purification three-axis angular velocity time series data that can accurately reflect the operating characteristics of the equipment is obtained.

[0018] This embodiment, by rigidly fixing a three-axis gyroscope sensor to a reciprocating motion component and simultaneously acquiring three-dimensional angular velocity data, can completely obtain the motion information of the device in three-dimensional space, avoiding the limitations of single-dimensional monitoring and ensuring the comprehensiveness and real-time nature of the raw data. Preprocessing the raw data by removing outliers and filtering noise can improve data quality and eliminate the influence of interference factors on subsequent analysis results.

[0019] In a preferred embodiment of the present invention, step 2 includes: Step 200a: Extract the triaxial angular velocity component sequence within a fixed analysis time window from the purified triaxial angular velocity time series data. Using the X-axis, Y-axis, and Z-axis angular velocity components as the three coordinate axes of a three-dimensional rectangular coordinate system, map the data points at each moment in the triaxial angular velocity component sequence to the original coordinate points in three-dimensional space, forming the original point set. Specifically, this includes: From the purified triaxial angular velocity time-series data obtained after processing in step 101, the triaxial angular velocity component sequence corresponding to a fixed analysis time window of 1 second is extracted to ensure that the extracted sequence can completely reflect the motion state of the equipment in the current time period. Using the X-axis angular velocity component in this sequence as the X-axis of a three-dimensional Cartesian coordinate system, the Y-axis angular velocity component as the Y-axis, and the Z-axis angular velocity component as the Z-axis, the triaxial angular velocity value corresponding to each sampling moment in the sequence is synchronously mapped to a coordinate point in three-dimensional space. All coordinate points corresponding to all sampling moments are summarized according to their temporal relationship to form the original point set for subsequent spatial analysis.

[0020] Step 201a: Based on the principal eigenvalues ​​of the covariance matrix between each data point in the original point set and the remaining data points in each data point's neighborhood, the mutual information in the time neighborhood, and the Euclidean distance gradient, determine the local clustering metric and the local discretization metric for each data point. Specifically, this includes: For each data point in the original point set, a corresponding spatial neighborhood and a temporal neighborhood are determined. The spatial neighborhood is centered on the current data point and includes the 15 nearest spatially adjacent data points within its radius. The temporal neighborhood is centered on the current data point and includes 10 temporally adjacent data points before and after it, forming a local analysis sample set including the current point. For all data points within the neighborhood, the average values ​​of the X, Y, and Z coordinate dimensions are calculated first. Based on these average values, a 3×3 covariance matrix is ​​calculated and constructed. Eigenvalue decomposition is performed on the covariance matrix to obtain eigenvalues ​​sorted from largest to smallest. The largest eigenvalue is taken as the principal eigenvalue of the current data point. Specifically, the decomposition process involves constructing the characteristic equation corresponding to the covariance matrix, i.e. =0, where It is a 3×3 covariance matrix; These are the eigenvalues ​​to be solved; It is the identity matrix of the same order as the covariance matrix. By solving the characteristic equation, all eigenvalues ​​corresponding to the covariance matrix are obtained. All eigenvalues ​​are arranged in descending order of value to form an ordered eigenvalue sequence. The eigenvalue that is the first and largest in the ordered eigenvalue sequence is selected as the principal eigenvalue corresponding to the current data point. The magnitude of this principal eigenvalue directly characterizes the degree of dispersion of the data points in the neighborhood along the main distribution direction in three-dimensional space. The larger the principal eigenvalue, the more dispersed the distribution of the neighborhood data points along the main extension direction and the stronger the spatial extension; the smaller the principal eigenvalue, the more concentrated the distribution of the neighborhood data points along the main direction and the higher the degree of aggregation. This achieves a quantitative characterization of the spatial aggregation and dispersion characteristics of the data points.

[0021] For the preceding and following sequences within the time neighborhood, calculate the information entropy of a single sequence and the joint entropy of the two sequences: ; ; in, , It is a sequence ,sequence Information entropy; It is a sequence and The joint entropy; , It is the probability distribution of a single variable; It is the joint probability distribution of two variables.

[0022] Mutual information is calculated based on information entropy and joint entropy, i.e., the mutual information between the beginning and end of a time series. For the current data point and other data points in its neighborhood, calculate the three-dimensional Euclidean distance, and perform a first-order difference operation on the three-dimensional Euclidean distance along the spatial position direction to obtain the Euclidean distance gradient. The local clustering metric is obtained by weighted summation of the principal eigenvalues, mutual information, and Euclidean distance gradient. ,in , , These are preset weighting coefficients, all positive numbers, with values ​​of 0.4, 0.3, and 0.3 respectively. + + =1, It is the largest principal eigenvalue of the covariance matrix. This refers to mutual information. Local clustering metrics quantify the density of a current data point within its neighborhood in three-dimensional space. A higher value indicates a denser distribution of data points and more concentrated motion trajectory features, effectively reflecting the stability and consistency of the device's motion state and providing a density judgment basis for subsequent adaptive adjustment of data point coordinates. The local dispersion metric is obtained by weighted summing the Euclidean distance gradient with the reciprocal of the principal eigenvalue. ,in , These are preset weighting coefficients, all positive numbers, with values ​​of 0.5 and 0.5 respectively. =1. Local discreteness measure is used to quantify the degree of dispersion and offset between the current data point and its neighboring data points in three-dimensional space. The larger the value, the looser the distribution of data points in space and the more obvious the offset. It can accurately identify distortion, fluctuation and abnormal discrete areas in the motion trajectory of the equipment, and provide a basis for judging the degree of dispersion for subsequent coordinate expansion adjustment and trajectory smoothing.

[0023] Step 202a: Based on the local clustering metric and local discrepancy metric of each data point, adaptively adjust the coordinates of each data point in the original point set. Apply a shrinkage offset along the direction of the corresponding point's neighborhood center for points whose local clustering metric is greater than the neighborhood mean, and apply an expansion offset along the direction of the corresponding point's neighborhood outside for points whose local discrepancy metric is greater than the neighborhood mean. Alternately perform shrinkage and expansion offsets until the coordinate change of all data points is less than a preset convergence threshold to obtain the adjusted point set, specifically including: Based on the local clustering metric and local discrepancy metric calculated for each data point in step 201a, an adaptive coordinate adjustment operation is performed on the coordinates of all data points in the original point set. First, all data points in the original point set are traversed, and the local clustering metric values ​​of all data points are summed. The sum is then divided by the total number of data points to obtain the overall neighborhood mean of the local clustering metric. Next, the local discrepancy metric values ​​of all data points are summed, and the sum is divided by the total number of data points to obtain the overall neighborhood mean of the local discrepancy metric. For data points whose local clustering metric is greater than the overall neighborhood mean, a fixed shrinkage offset of 0.0005 is applied along the geometric center direction of their neighborhood, causing these data points to cluster towards the neighborhood center. For data points whose local discrepancy metric is greater than the overall neighborhood mean, a fixed expansion offset of 0.0005 is applied along the outer normal direction of their neighborhood, causing these data points to expand appropriately outwards from their neighborhood. The contraction and expansion offset operations are performed alternately and cyclically according to the above rules. After each adjustment, the coordinates of each data point after the current iteration are compared with the coordinates of the previous iteration. The absolute value of the difference is taken to obtain the coordinate change. The coordinate change of each data point is checked one by one. When the coordinate change of all data points is less than 0.001, it is determined that the preset convergence threshold has been reached, and the coordinate adjustment iteration process is stopped, resulting in an adjusted point set with a more stable coordinate distribution that better reflects the true motion characteristics.

[0024] Step 203a involves connecting the data points in the adjusted point set sequentially with smooth curves according to the original time order to obtain a first-order continuously differentiable three-dimensional spatial trajectory. This three-dimensional spatial trajectory is defined as a three-dimensional phase space trajectory, specifically including: Arrange all data points in the adjusted point set obtained in step 202a according to the original time sequence of data acquisition. Using time series t as the independent variable and the three-dimensional coordinates (x, y, z) of each data point as the dependent variable, establish a piecewise cubic polynomial. This ensures that the function values ​​and first-order derivatives of adjacent piecewise curves are continuous at the connection points, achieving first-order continuous differentiability. The expression for a single-segment cubic spline function is as follows: ; in, It is the first Segment cubic spline curve, It is a time-series independent variable; It is the first Each segment node time; These are the coefficients of the polynomial to be determined. The continuity constraint for adjacent piecewise segments is: ,in It is the first Segment curves at connecting nodes The function value at that point, i.e., the three-dimensional coordinate value corresponding to that time series; It is the first +1 curve segment at the connecting node The function value at that point, and Equal to ensure that adjacent curve segments are connected continuously without breaks at the nodes; It is the first The first derivative of a segment curve characterizes the curve in the corresponding time series. The rate of change at a point reflects the smoothness of the trajectory; It is the first Segment curves at connecting nodes The first derivative value at; It is the first Segment curves at connecting nodes The first derivative value at that point, and Equal to ensure a smooth transition between adjacent curve segments at nodes, achieving first-order continuous differentiability; It is the first Section and the The time of the common connection node of the segments is the connection time of the two curve segments; To determine the coefficients of each piecewise cubic polynomial (Solve separately for each of the three coordinate dimensions), it is necessary to establish and solve the three bending moment equations. The formulas for the three bending moment equations are as follows (applicable to a single coordinate dimension): ; in It is the first Nodes The second derivative at the point corresponds to independent bending moment values ​​in each of the three coordinate dimensions; It is the first The time interval length of the segment curve, i.e. = - , which is the difference between two adjacent time points, and is always a positive number; It is the first The time interval length of the segment curve, i.e. = - ; It is the first Nodes The dependent variable value at a given location corresponds to any one of the three dimensions (X, Y, Z) in the three-dimensional coordinate system. When substituting these values, they can be replaced with the values ​​for the three dimensions respectively. (X-axis coordinate) (Y-axis coordinate) (Z-axis coordinates), used to solve for the bending moment in the corresponding dimension; It is the first Nodes The dependent variable value at the location (corresponding to the same coordinate dimension); It is the first Nodes The dependent variable value at that location (corresponding to the same coordinate dimension).

[0025] Based on the time series values ​​of all nodes Calculate the length of each time series interval Then, combining the three-dimensional coordinate values ​​corresponding to each node, a system of three bending moment equations is established for the X, Y, and Z dimensions respectively, and the bending moment at each node is obtained by solving the equations. According to bending moment With piecewise polynomial coefficients Based on the correspondence, the coefficient values ​​of all segments are calculated. The specific correspondence is as follows (taking a single coordinate dimension as an example; this applies to all three dimensions X, Y, and Z, only the parameter symbols for the corresponding dimensions need to be replaced). For the th Piecewise cubic spline curves, their coefficients The correspondence between the bending moment value, coordinate value, and time interval length of the corresponding node in this segment is calculated using the following formula: , , , ; in, It is the first The bending moment values ​​at each node (corresponding to the same coordinate dimension).

[0026] Based on the aforementioned clear correspondence, for the X, Y, and Z dimensions, the obtained bending moment values ​​at each node, the length of each time interval, and the corresponding node coordinate values ​​are substituted into the above coefficient calculation formula to calculate the polynomial coefficients for all segments in each of the three dimensions. Using the piecewise cubic polynomials with the obtained coefficients, the data points are smoothly connected sequentially according to the time sequence, forming a continuous and smooth transition curve between adjacent data points. Finally, a three-dimensional spatial trajectory with first-order continuous differentiable properties is generated, and this trajectory is formally defined as a three-dimensional phase space trajectory used to characterize the motion state of the equipment.

[0027] Step 200b: Calculate the local curvature value of each data point on the three-dimensional phase space trajectory, and select the top N data points with the largest local curvature values ​​as a candidate pole set; from the candidate pole set, select the data point with the largest curvature value and the earliest corresponding time, and determine it as the first phase pole representing the abrupt change in the geometric shape of the motion trajectory, specifically including: Traverse all data points on the three-dimensional phase space trajectory generated in step 203a, and calculate the local curvature value corresponding to each data point. For any point to be calculated on the trajectory, select that point and its two temporally adjacent data points, and denote these three spatial points as Q1( ), Q2( ), Q3( Construct two spatial vectors derived from the three points: a vector pointing from point Q1 to point Q2, and a vector pointing from point Q1 to point Q3. After constructing these two vectors, calculate their cross product to obtain the normal vector of the spatial plane containing the three points. (A, B, C), through this normal vector and any one of the three points (here, Q1 is chosen). The equation of a plane is uniquely determined by the fact that three points are coplanar. The formula for the plane equation is as follows: A(u- )+B(v- )+C(w- )=0; After determining the plane equation, the three spatial points are projected onto the plane, transforming the three-dimensional circular arc fitting problem into a three-point concyclic fitting problem in a two-dimensional plane, thus simplifying the calculation process. Let Oc(p, q) be the center of the fitted circle in the projected two-dimensional plane, and let the radius be... Based on the core geometric condition that all three points lie on the circumference and are equidistant from the center of the circle using Euclidean distances, the following system of circle fitting equations is constructed: ; Expanding the above system of equations yields three quadratic equations. Then, subtracting each pair of equations eliminates the quadratic equations. The terms and quadratic terms are transformed into a system of two linear equations in two variables about the center coordinates p and q. Solving this system of linear equations yields the specific coordinates of the center Oc(p, q) of the fitted circle. The center coordinates are then substituted back into any original equation to calculate the radius r of the fitted circle. This radius is the radius of the arc obtained by fitting the three spatial points. The reciprocal of this arc radius is used as the local curvature value of the current point to be calculated, which is used to quantify the curvature of the three-dimensional phase space trajectory at this position. After completing the local curvature calculation of all data points, all data points are sorted in descending order of local curvature value. The top N (10) data points with the largest local curvature value are selected to form a candidate pole set. In this candidate pole set, data points that simultaneously meet the two conditions of the largest local curvature value and the earliest original data acquisition time are further selected and determined as the first phase pole. This first phase pole is used to characterize the key feature position where the geometric shape of the motion trajectory undergoes a significant change.

[0028] Step 201b: Calculate the energy density distribution of the three-dimensional phase space trajectory within the analysis time window, determine the location of the energy density centroid, and identify the trajectory data point closest to the energy density centroid as the second phase pole characterizing the centroid of the motion energy distribution. In the three-dimensional phase space trajectory, select the trajectory data point corresponding to the midpoint of the analysis time window as the third reference node for calibrating the motion state time series reference. Specifically, this includes: Based on the three-dimensional phase space trajectory with first-order continuous differentiability generated in step 203a, the energy density distribution of the three-dimensional phase space trajectory within the fixed analysis time window of 1 second is calculated, and then the second phase pole and the third reference node are determined. Specifically, all data points within the fixed analysis time window of 1 second on the three-dimensional phase space trajectory are selected as calculation samples. These data points all contain the corresponding three-dimensional coordinates of the X-axis, Y-axis and Z-axis, and the timing corresponds to the order of data acquisition.

[0029] The energy density of each data point on the trajectory is calculated. Energy density characterizes the concentration of energy in the device's motion at the corresponding moment. Specifically, for each trajectory data point, the squares of its X-axis, Y-axis, and Z-axis angular velocities are calculated, and the sum of these squares is the energy density of the current data point. The angular velocities of each axis are obtained by calculating the first derivative of the piecewise cubic spline function fitted in step 203a. The moment corresponding to the derivative calculation is the acquisition moment of that data point. Following the above calculations, the energy density of all trajectory data points within a fixed 1-second analysis time window is calculated one by one. The energy densities of all data points together constitute the complete energy density distribution within that time window. The larger the energy density value, the more concentrated the energy of the device's motion at the corresponding moment.

[0030] The three-dimensional coordinates of the energy density centroid are calculated. The energy density centroid is the center of gravity of the device's kinetic energy within the entire time window, and its coordinates represent the concentrated location of the device's kinetic energy within the entire time window. After the energy density centroid is calculated, all data points on the three-dimensional phase space trajectory within the 1-second fixed analysis time window are traversed. The three-dimensional Euclidean distance from each trajectory data point to the energy density centroid is calculated one by one. All distance values ​​are compared, and the trajectory data point with the smallest distance is selected and determined as the second phase pole. The second phase pole is used to represent the center of gravity of the device's kinetic energy distribution and provides a core reference node for subsequent phase shift, phase deflection, and other analyses. The start and end times of the 1-second fixed analysis time window are defined. In the three-dimensional phase space trajectory, the trajectory data point corresponding to this intermediate time is found and determined as the third reference node. The third reference node is used to provide a unified time series reference for phase analysis in subsequent steps (such as phase deflection angle calculation).

[0031] Step 202b: After determining the three poles, the three-dimensional phase space trajectory is discretized into a three-dimensional grid point set. The shortest Euclidean distance from each grid point inside the three-dimensional phase space trajectory to the boundary of the three-dimensional phase space trajectory is calculated to form a range field. Wavefronts propagating simultaneously inward from the boundary of the three-dimensional phase space trajectory are simulated. The positions in the range field with local maxima where the wavefronts first meet are marked as central axis points. All central axis points constitute a central axis point set, specifically including: The three-dimensional phase space trajectory obtained in step 203a is discretized using a uniform grid partitioning method. The spatial range containing the three-dimensional phase space trajectory is divided into a uniformly distributed set of three-dimensional grid points according to a preset uniform grid step size (the grid step size is set to 0.001, consistent with the coordinate adjustment convergence threshold). Each grid point corresponds to a unique three-dimensional spatial coordinate. At the same time, the coordinates of all grid points are normalized to eliminate the influence of differences in the dimensions of different coordinates on subsequent distance calculations and wavefront simulations. The processing formula is as follows: ; in, , , These are the normalized 3D coordinates of the grid points; , , These are the original 3D coordinates of the grid points; , , It is the minimum value of the coordinates corresponding to all data points of the three-dimensional phase space trajectory; , , It is the maximum value of the coordinates corresponding to all data points of the three-dimensional phase space trajectory.

[0032] After discretization and dimensionless transformation, for each grid point within the trajectory, the Euclidean distance to all boundary points (trajectory start point, end point, and edge extreme points) is calculated, and the minimum value is taken as the shortest distance from that point to the trajectory boundary. The shortest distances of all grid points constitute a distance field. The smaller the distance, the closer it is to the boundary; the larger the distance, the closer it is to the center; and a distance of 0 indicates that it falls on the boundary. Through this numerical mapping relationship, the relative positional distribution of each grid point within the entire trajectory space can be clearly and quantitatively presented.

[0033] Next, wavefront simulation is performed to extract the central axis point. The wavefront propagates into the trajectory at a uniform speed of 0.001 unit length / unit time, consistent with the grid step size. All boundary points (trajectory start point, end point, and edge extreme points) are used as wavefront starting points, and synchronous triggering control is implemented for all starting points. By monitoring the propagation distance of each wavefront in real time (propagation distance = speed × time), a uniform propagation state is maintained, with the wavefront advancing into the trajectory by 0.001 unit length per unit time, ensuring that wavefronts starting from different directions and different starting points always advance synchronously without any leading or lagging phenomena.

[0034] The wavefront propagation direction follows the principle of being perpendicular to the trajectory boundary. The propagation direction of each wavefront element is perpendicular to the tangent direction of the trajectory boundary at its corresponding position, and it only propagates into the three-dimensional phase space trajectory, strictly prohibiting diffusion outward from the trajectory. Specifically, the propagation direction is restricted by presetting the normal vector direction of the trajectory boundary. First, the three-dimensional coordinates of each boundary point and its two adjacent boundary points are extracted, and the trajectory tangent vector at that boundary point is constructed (calculated from the difference in coordinates between two adjacent points, consistent with the direction of the first derivative of the piecewise cubic spline curve fitted in step 203a at that point); the perpendicular vector of the trajectory tangent vector is selected (in three-dimensional space, the vector perpendicular to the tangent vector point). The vector with a product of 0 is used as a reference. The four fingers of the right hand are bent along the tangent direction, and the direction pointed by the thumb is the direction inside the trajectory. The vector pointed by the thumb is the preset normal vector direction. The calculated normal vector is normalized (the normal vector magnitude is adjusted to 1. The normalization process is to divide the three components of the original normal vector by their magnitudes to obtain the unit normal vector. At this time, the magnitude of the unit normal vector is always 1). This ensures that the normal vector direction of each boundary point is uniform and the magnitude is consistent. The propagation direction of the wavefront element is consistent with the direction of the normal vector, which restricts the wavefront to propagate only into the trajectory and prevents it from spreading outward.

[0035] During propagation, when the trajectory boundary is a straight segment, the normal vectors of all boundary points in that segment remain parallel, and the corresponding wavefront infinitesimal propagation directions are consistent, thus maintaining a straight wavefront shape. When the trajectory boundary is a curved segment, the normal vectors of the boundary points in that segment change synchronously with the trajectory curvature (the normal vector of each boundary point always points inside the trajectory and is perpendicular to its own tangent direction), and the corresponding wavefront infinitesimal propagation direction adaptively adjusts with the normal vector direction, so that the wavefront synchronously presents an arc shape corresponding to the curvature of the trajectory boundary, always maintaining a perpendicular relationship with the trajectory boundary. At the same time, during the entire wavefront simulation process, the distance field is synchronously sampled and monitored in real time, with the sampling frequency consistent with the wavefront propagation time unit (sampling once per unit time). The wavefront arrival time, wavefront source (i.e., which boundary point the wavefront originates from), and the change of the shortest Euclidean distance for each grid point are recorded in real time.

[0036] During wavefront simulation, the central axis point is simultaneously selected from the range field, which must satisfy two core characteristics. The first characteristic is a local maximum, which means that the adjacent grid points are the eight grid points directly adjacent to the current location in the X, Y, and Z axes (i.e., all grid points in 3D space whose coordinate differences are within one grid step, and whose coordinate differences satisfy the absolute values ​​of the coordinate differences in the X-axis direction |Δx|≤0.001, the Y-axis direction |Δy|≤0.001, and the Z-axis direction |Δz|≤0.001, excluding the current location itself). The shortest Euclidean distance of this location is extracted, as are the shortest Euclidean distances of its eight adjacent grid points. If the shortest Euclidean distance of this location is greater than the shortest Euclidean distances of all its adjacent grid points, then this location is determined to be a local maximum point; if the shortest Euclidean distance of any adjacent grid point is greater than or equal to the shortest Euclidean distance of this location, then this location does not satisfy the local maximum condition.

[0037] The second characteristic is the first encounter of wavefronts. This involves recording the arrival time of wavefronts at each grid point. If a point is simultaneously reached by wavefronts originating from two or more different boundary points (with an arrival time difference ≤ 0.0001), and this is the first time the point has been simultaneously reached by multiple wavefronts (previously it had not been reached by any wavefront or only by a single wavefront), then it is considered the first encounter point of the wavefronts. Only grid points that simultaneously meet both of these conditions are marked as central axis points. After marking, every grid point in the entire distance field is traversed, and each point is judged against the above selection criteria. After collecting all central axis points that meet the conditions, the criterion for determining adjacent central axis points is defined as follows: the three-dimensional Euclidean distance between two central axis points ≤ 2 times the grid step size (i.e., ≤ 0.002 unit length). Central axis points that meet this condition are classified as adjacent nodes. Finally, all central axis points are arranged in an orderly manner according to this proximity relationship, forming a complete set of central axis points. This set of central axis points precisely corresponds to the core centerline position of the three-dimensional phase space trajectory.

[0038] Step 203b: Connect adjacent central axis points and determine the connection path based on the gradient direction of the distance field to construct a skeleton graph with a tree-like topology. For the tree-like skeleton structure, calculate the path length from the root node to each leaf node, and take the length of the longest path as the trunk length. Count the number of all branch nodes in the skeleton graph as the number of branches. At each branch node, fit the direction vectors of the two adjacent skeleton segments of the corresponding node, calculate the angle between the direction vectors of the two skeleton segments, and define the angle as the curvature change at the corresponding branch point. Calculate the arithmetic mean of the curvature changes at all branch points, denoted as the average branch point curvature. Use the trunk length, the number of branches, and the average branch point curvature together as geometric invariant features of the three-dimensional phase space trajectory morphology, specifically including: The spatial positions of all the central axis points obtained in step 202b are analyzed to clarify the three-dimensional coordinate information of each central axis point. Adjacent central axis points are identified according to a preset adjacency criterion: the three-dimensional Euclidean distance between two central axis points is less than or equal to twice the grid step size (the grid step size remains consistent with step 202b, set to 0.001 units, meaning the distance between adjacent central axis points is ≤ 0.002 units). For the identified adjacent central axis points, they are connected in an orderly manner according to their spatial proximity, forming preliminary central axis point connection paths, ensuring that each adjacent central axis point has a unique initial connection path. The gradient direction of the distance field constructed in step 202b is calculated. First, the distance value D of each grid point in the distance field is determined (i.e., the shortest Euclidean distance from the grid point to the boundary of the three-dimensional phase space trajectory, already calculated in step 202b). The partial derivatives of each grid point in the X, Y, and Z axes are calculated using the central difference method. The partial derivatives in the three axes are combined to obtain the gradient vector of the grid point. The direction of this gradient vector is the gradient direction of the distance field. The core characteristic of the distance field gradient direction is that it points in the direction of increasing distance field value. Since the larger the distance field value, the closer the corresponding grid point is to the center region of the 3D phase space trajectory, the gradient direction of the distance field is the direction inside the trajectory. Based on this gradient direction, the paths of the initially connected central axis points are optimized one by one, eliminating unreasonable cross connections and redundant connections. Cross connections refer to connections between different paths that intersect each other and have no practical significance, while redundant connections refer to unnecessary connections that take detours. The optimized connection paths must satisfy the shortest path principle and always conform to the trend of the trajectory's central region. Finally, all central axis points are connected in an orderly manner through the optimized paths to construct a skeleton diagram with a tree-like topology. This skeleton diagram can accurately represent the core morphological features of the 3D phase space trajectory.

[0039] Calculate the 3D Euclidean distance from each central axis point in the central axis point set to the first phase pole (determined in step 200b), and select the central axis point with the smallest distance as the root node of the skeleton graph. Traverse the entire tree-like skeleton structure and filter out central axis points with no subsequent connected nodes, i.e., nodes that are only connected to one preceding central axis point and have no subsequent central axis points connected to them. These nodes are the leaf nodes of the skeleton graph. Calculate the path length from the root node to each leaf node. After all the path lengths from the root node to the leaf node have been calculated, compare all the path lengths one by one and select the path length with the largest value. This path length is determined as the trunk length of the skeleton graph. The trunk length is used to characterize the main extension direction and overall scale of the 3D phase space trajectory.

[0040] The branch count is calculated for all bifurcated nodes in the skeleton diagram. A bifurcated node is defined as a central point connecting two or more subsequent nodes; that is, in addition to connecting to the preceding central point, it also connects to two or more subsequent central points, forming a branch structure. The entire tree-like skeleton structure is traversed, and all bifurcated nodes that meet the definition are identified one by one. The total number of these nodes is the branch count of the skeleton diagram. The branch count characterizes the complexity of the 3D phase space trajectory. A higher branch count indicates more bifurcations during the extension of the 3D phase space trajectory, resulting in a more complex trajectory shape and indicating a more unstable motion state and more complex motion laws for the moving parts. Conversely, a lower branch count indicates fewer bifurcations in the 3D phase space trajectory, with the trajectory shape closer to a straight or simple branch structure, indicating a more stable motion state and simpler motion laws for the moving parts. Its characterization logic directly corresponds to the actual morphological characteristics of the 3D phase space trajectory, accurately reflecting the differences in trajectory complexity.

[0041] The curvature change at each bifurcation node is calculated as follows: A bifurcation node is selected, and its normalized coordinates are denoted as... , , The normalized coordinates of its two adjacent central axis points are respectively , , (Precedence midline point) and , , (Subsequent central axis point), the two skeleton segments are the skeleton segment formed by the bifurcation node and the preceding central axis point, and the skeleton segment formed by the bifurcation node and the subsequent central axis point, respectively. The direction vectors of the two skeleton segments are obtained by fitting, and are denoted as vectors respectively. sum vector ,in = , = Calculate the dot product of the two direction vectors and their magnitudes. Then, calculate the cosine of the angle between the two direction vectors and perform an inverse cosine operation to obtain the spatial angle between them. This spatial angle is defined as the curvature change at the current branch point. The curvature change characterizes the degree of morphological abrupt change at the bifurcation point of the skeleton. The larger the curvature change angle, the more significant the directional difference between the two adjacent skeleton segments at the bifurcation point, and the more drastic the skeleton's turn at that position, i.e., the more obvious the morphological abrupt change of the three-dimensional phase space trajectory at that bifurcation point. The smaller the curvature change angle, the closer the directions of the two adjacent skeleton segments at the bifurcation point, and the gentler the skeleton's turn at that position, i.e., the milder the morphological abrupt change of the three-dimensional phase space trajectory at that bifurcation point. Its characterization logic is directly related to the actual morphological change at the trajectory bifurcation point and can accurately reflect the morphological abrupt change difference at the bifurcation point. Following the same steps, calculate the curvature change at all bifurcation points one by one. The arithmetic mean of the curvature changes (spatial angles) at all branch points is calculated and denoted as the average branch point curvature. The three scalars, trunk length, number of branches, and average branch point curvature, are used together as geometric invariant features of the trajectory morphology in three-dimensional phase space. This feature is scale independent and can stably characterize the core geometric morphology of the trajectory.

[0042] Step 200c: Calculate the absolute values ​​of the coordinate differences between the first and second phase poles along the X, Y, and Z axes, and define these three absolute values ​​as the X-axis phase offset modulus, Y-axis phase offset modulus, and Z-axis phase offset modulus, respectively. Simultaneously, calculate the Euclidean distance between the first and second phase poles in the three-dimensional Cartesian coordinate system, defining it as the total phase offset modulus, used to quantify the attitude divergence of the moving part at the microscale. Specifically, this includes: The three-dimensional coordinate data of the first phase pole and the second phase pole are read. The coordinates of the two poles have been time-aligned and dimensionless to eliminate the influence of dimensional differences on the calculation results. The absolute values ​​of the coordinate differences between the first and second phase poles along the X, Y, and Z axes are calculated respectively. These three absolute values ​​are defined as the X-axis phase offset modulus, Y-axis phase offset modulus, and Z-axis phase offset modulus, respectively. The X-axis phase offset modulus represents the degree of offset between the two poles along the X-axis; the Y-axis phase offset modulus represents the degree of offset between the two poles along the Y-axis; and the Z-axis phase offset modulus represents the degree of offset between the two poles along the Z-axis. The total phase offset modulus is the Euclidean distance between the first and second phase poles in the three-dimensional Cartesian coordinate system. The total phase offset modulus is used to comprehensively quantify the degree of attitude divergence exhibited by the moving part at the microscale. The larger the total phase offset modulus value, the greater the spatial distance between the two phase poles and the more obvious the attitude divergence of the moving part; the smaller the total phase offset modulus value, the closer the two phase poles are and the more stable the attitude of the moving part.

[0043] Step 201c: Construct a first vector pointing from the first phase pole to the second phase pole, and a second vector pointing from the second phase pole to the third reference node; calculate the angle between the first and second vectors, and define the angle as the phase deflection angle, which is used to quantify the degree of deviation of the main vibration direction of the moving part from the timing reference, specifically including: The three-dimensional coordinate information of three key nodes is defined, including the first phase pole and the second phase pole. The coordinates of the two poles have been determined in steps 200b and 201b. The third reference node is the trajectory data point corresponding to the midpoint of the fixed analysis time window determined in step 201b. Its coordinates have been recorded synchronously to provide a unified time series reference. The coordinates of all three nodes have been normalized. Two core vectors are constructed, corresponding to the main vibration direction of the moving component and the timing reference direction, respectively. The first vector, the first vector, is constructed by starting from the first phase pole and ending at the second phase pole, pointing from the first phase pole to the second phase pole. This vector characterizes the core trend of the main vibration direction of the moving component and is a key vector reflecting the main vibration direction. The direction of this vector directly corresponds to the core direction of the main vibration of the moving component, and the magnitude of the vector corresponds to the overall amplitude of the main vibration. The larger the vector magnitude, the larger the amplitude of the main vibration. The second vector, the second vector, is constructed by starting from the second phase pole and ending at the third reference node, pointing from the second phase pole to the third reference node. This vector characterizes the direction of the timing reference and is the core direction of the timing reference. It serves as a reference standard for measuring the degree of deviation of the main vibration direction of the moving component. The vector direction is fixed, and the deviation of the main vibration direction from the timing reference can be quantified by the difference between its direction and that of the first vector.

[0044] After constructing the two vectors, the dot product of the first and second vectors and their magnitudes are calculated. The ratio of the dot product to the product of the magnitudes of the two vectors is then calculated to obtain the cosine of the angle between the two vectors. An inverse cosine operation is performed on the cosine to obtain the spatial angle between the two vectors. This spatial angle is defined as the phase deflection angle, in radians, and is used to quantify the degree of deviation of the main vibration direction of the moving component from the timing reference. The closer the phase deflection angle is to 0, the more consistent the main vibration direction of the moving component is with the timing reference, and the smaller the deviation. This indicates that the main vibration direction of the moving part deviates more severely from the timing reference.

[0045] Step 200d involves extracting the time series sequences of the X-axis phase shift modulus, Y-axis phase shift modulus, and Z-axis phase shift modulus, respectively; calculating the Pearson cross-correlation coefficients between each pair of the X-axis, Y-axis, and Z-axis phase shift modulus time series sequences, and constructing a cross-correlation coefficient matrix, specifically including: Based on the chronological order of data acquisition, the time series sequences corresponding to the X-axis phase offset modulus, Y-axis phase offset modulus, and Z-axis phase offset modulus calculated in step 200c are extracted respectively. The X-axis phase offset modulus time series sequence is denoted as... The Y-axis phase offset modulus time series is denoted as The Z-axis phase offset modulus time series is denoted as ,in, This represents the length of the time series, i.e., the total number of data acquisitions.

[0046] Calculate the Pearson cross-correlation coefficient between any two time series. There are three combinations of the three axial phase shift modulus time series: X and Y, X and Z, and Y and Z. The calculation steps for these three combinations are identical. The Pearson cross-correlation coefficient ranges from [-1, 1]. The closer the coefficient is to 1, the stronger the positive correlation between the two series; the closer it is to -1, the stronger the negative correlation; and the closer it is to 0, the weaker the correlation. Three cross-correlation coefficients are obtained. , , Combining the self-correlation coefficient (which is perfectly correlated with itself and always has a value of 1), a 3×3 three-axis phase offset modulus cross-correlation matrix is ​​constructed, with the matrix form as follows: ; in, , , The diagonal elements of the matrix are the cross-correlation coefficients of the phase offset modulus time series of the X-axis, Y-axis, and Z-axis, respectively, and their values ​​are always 1. , , , , , Let be the off-diagonal elements of the matrix, representing the cross-correlation coefficients of the two corresponding time series sequences. Due to the symmetry of the cross-correlation coefficients, the following condition is satisfied: This matrix comprehensively presents the correlation distribution among the phase shift moduli of the three axes through the magnitude and sign of the matrix elements. The value of each off-diagonal element in the matrix directly corresponds to the correlation degree of the time series of the phase shift moduli of the two axes corresponding to its row and column. The closer the value is to 1, the stronger the correlation between the two axes; the closer the value is to -1, the stronger the negative correlation between the two axes; and the closer the value is to 0, the weaker the correlation between the two axes. Through the symmetrical distribution characteristics of the matrix, the correlation difference between any two axes can be intuitively compared, clearly presenting the distribution of the correlation strength between the X-axis and Y-axis, X-axis and Z-axis, and Y-axis and Z-axis, thereby comprehensively grasping the correlation relationship among the phase shift moduli of the three axes.

[0047] Step 201d: Calculate the mean of the off-diagonal elements in the cross-correlation matrix, using it as a multidimensional motion coupling index characterizing the tightness of topological constraints on multi-source sensor data in the feature space; arrange the X-axis phase offset modulus, Y-axis phase offset modulus, Z-axis phase offset modulus, total phase offset modulus, phase deflection angle, multidimensional motion coupling index, and geometric invariant features in a preset order, and combine them to generate a feature index vector, specifically including: Read the 3×3 cross-correlation matrix constructed in step 200d, and filter out all off-diagonal elements, totaling 6 elements. Cross-correlation coefficients are symmetric. The values ​​of the 6 filtered off-diagonal elements are summed sequentially to obtain the total off-diagonal element sum. This sum is divided by 6 (the number of off-diagonal elements) to obtain the arithmetic mean. This mean is determined as the multidimensional motion coupling index, with a value range of [-1, 1]. The closer the index is to 1, the stronger the correlation between the phase offset moduli of each axis, and the tighter the topological constraint of the multi-source sensor data; the closer the index is to 0, the weaker the correlation between the phase offset moduli of each axis, and the looser the topological constraint of the multi-source sensor data.

[0048] The following characteristic parameters calculated in the previous steps are arranged in a fixed order: X-axis phase offset modulus, Y-axis phase offset modulus, Z-axis phase offset modulus, total phase offset modulus, phase deflection angle, multidimensional motion coupling index, trunk length (scalar), number of branches (scalar), and mean curvature of branch points (arithmetic mean of the angles of curvature change at all branch points, scalar). These nine parameters are combined in this order to generate a characteristic index vector. This vector covers information in four dimensions: equipment micro-attitude divergence (the first four items), deviation of the main vibration direction (the fifth item), multi-axis motion coordination (the sixth item), and motion trajectory geometry (the last three items). The determination of this order follows the core logic of moving from basic to comprehensive and from local to global, and aligns with the representation priority of the equipment's operating status. First, the phase offset moduli of the X, Y, and Z axes are arranged, as they are fundamental parameters for quantifying the divergence of the motion component's attitude, directly reflecting the core offset information of the device's micro-attitude in each axis. Next, the total phase offset modulus is arranged to comprehensively quantify the overall attitude divergence of the motion component in three-dimensional space, supplementing the comprehensive spatial deviation information that cannot be reflected by single-axis offset. Then, the phase deflection angle is arranged to characterize the deviation of the main vibration direction of the motion component from the time-series reference, improving the directional characteristics of the device's motion state. Then, the multidimensional motion coupling index is arranged to comprehensively reflect the tightness of the topological constraints of multi-source sensor data, reflecting the coordination and consistency of the device's multi-axis motion. Finally, the three geometrically invariant features of trunk length, number of branches, and mean curvature of branch points are arranged, as they can respectively characterize the main extension scale, complexity, and morphological change intensity at the branch points of the three-dimensional phase space trajectory, reflecting the regularity and stability of the device's motion from the overall morphology, achieving comprehensive coverage from local parameters to overall morphology.

[0049] This feature index vector covers multiple core dimensions of equipment operating status. It characterizes the divergence of the equipment's micro-attitude in each axis using the phase offset modulus along the X, Y, and Z axes. A larger offset modulus value for a given axis indicates a more significant offset of the moving parts in that axis, resulting in more severe overall attitude divergence; a smaller value indicates more stable attitude. The total phase offset modulus characterizes the comprehensive attitude divergence of the moving parts in three-dimensional space. A larger value indicates a greater spatial distance between the two phase poles, resulting in more significant attitude divergence; a smaller value indicates a more concentrated and stable attitude. The phase deflection angle characterizes the deviation of the principal vibration direction from the time series reference. The closer the phase deflection angle is to 0, the more consistent the principal vibration direction is with the time series reference, and the smaller the deviation; the closer to 0, the smaller the deviation. (180°), the greater the deviation. The degree of collaborative constraint of multi-source data is characterized by the multi-dimensional motion coupling degree index. The closer the index is to 1, the stronger the correlation of the phase offset modulus of each axis, the tighter the topological constraint of the multi-source sensor data, and the better the coordination of equipment motion; the closer it is to 0, the looser the constraint and the worse the coordination. The overall shape of the equipment motion trajectory is characterized by the trunk length, the number of branches, and the mean curvature of the branch points. The longer the trunk length, the wider the trajectory extension range; the more branches, the more complex the trajectory shape; the larger the mean curvature of the branch points, the more obvious the abrupt change in shape at the trajectory bifurcation; conversely, the trajectory is smoother and simpler. The parameters complement each other and progress layer by layer, from the four core dimensions of micro attitude (axial offset and comprehensive offset), directional deviation, data coordination, and overall shape, to achieve a comprehensive, stable, and accurate characterization of the equipment's operating status.

[0050] This embodiment employs phase space reconstruction combined with local aggregation and local discretization metrics to adaptively adjust the coordinates of data points, thereby enhancing the spatial distribution characteristics of the data points. By determining the phase poles and reference nodes, it is possible to locate abrupt changes in trajectory morphology, the center of gravity of energy, and the time series reference point. By calculating the phase offset modulus of each axis, the total phase offset modulus, and the phase deflection angle, early latent anomalies in components such as transmission mechanisms and guide rails can be detected. By constructing a cross-correlation matrix of the three-axis phase offset modulus and calculating the multidimensional motion coupling index, coupling faults that cannot be detected by single-axis monitoring can be effectively identified.

[0051] In a preferred embodiment of the present invention, step 3 includes: Step 300: Call the pre-built feature benchmark library. The feature benchmark library stores the standard state vectors and statistical distribution parameters of the device when it is running in a healthy state, specifically including: The equipment was calibrated to factory standard parameters, confirming that key moving parts were free of abnormalities and that the three-axis gyroscope's acquisition accuracy met standards. The acquisition frequency was set to 200Hz, with each test lasting 30 minutes, and 30 sets of repeated tests were conducted (with 10-minute intervals between each set) to ensure data representativeness. Each set of tests acquired X, Y, and Z-axis coordinate data of the three-dimensional phase space trajectory, as well as process parameters. A fixed sliding time window of 20 sampling points combined with the 3σ criterion was used to remove outlier data points. After time alignment, min-max normalization was used to eliminate dimensional differences. For each set of preprocessed data, the following feature parameters were extracted according to steps 200a to 201d: X-axis phase offset modulus, Y-axis phase offset modulus, Z-axis phase offset modulus, total phase offset modulus, phase deflection angle, multidimensional motion coupling index, trunk length, number of branches, and mean curvature of branch points.

[0052] After extraction, the nine feature parameters of each group of experiments are arranged in the above order to form a single standard state vector. A total of 30 standard state vectors are generated from the 30 groups of experiments, forming a standard state vector set. Based on these 30 standard state vectors, the statistical distribution parameters of each feature component are calculated one by one, i.e., the mean, standard deviation, and variance of each component are calculated; the covariance matrix (9×9) among the nine components is also calculated. The standard state vector set (30 vectors) and the calculated mean, standard deviation, and covariance matrix are stored uniformly in a relational database to form a pre-built feature benchmark library. When called, data is read through an interface, and the number of vector components (must be 9), their arrangement order, and the completeness of the statistical parameters are verified.

[0053] Step 301: Calculate the Mahalanobis distance between the current feature index vector and the standard state vector in the feature benchmark library to obtain the real-time health score of the device, specifically including: The feature index vector generated in step 201d is normalized using the min-max normalization method, which maps the values ​​of each component in the feature index vector to the interval [0, 1]. The specific calculation formula is as follows: ; in, This represents the original value of a component in the feature index vector. This is the minimum value of this component in the standard state vector of the feature benchmark library. This is the maximum value of the component in the standard state vector of the feature benchmark library. This is the normalized value of the component. Following this formula, each component in the feature index vector is normalized one by one to obtain the normalized feature index vector. Simultaneously, all standard state vectors in the feature benchmark library are extracted, and their components are processed using the same min-max normalization method to ensure consistency with the normalization standard of the current feature index vector, resulting in a set of normalized standard state vectors. The Mahalanobis distance between the normalized feature index vector obtained by applying min-max normalization to the current equipment operating state feature index vector in step 301 and the normalized standard state vector set obtained by applying the same min-max normalization method to all standard state vectors in the feature benchmark library is calculated. This distance quantifies the degree of difference between the two. The smaller the Mahalanobis distance value, the smaller the difference between the current feature index vector and the standard state vector in a healthy state, and the closer the equipment operating state is to a healthy state; the larger the Mahalanobis distance value, the greater the difference between the two, and the more abnormal the equipment operating state. The Mahalanobis distance is converted into a real-time device health score. The conversion logic is a reverse mapping, meaning that the smaller the Mahalanobis distance, the higher the health score. The specific conversion formula is as follows: Real-time health score = 100 - ×100, where the real-time health score ranges from [0, 100]. The preset maximum Mahalanobis distance threshold is set to 5.0. (This clarifies the criteria for determining equipment failure status. When key moving parts of the equipment exhibit abnormalities such as wear, jamming, or excessive vibration, or when the three-dimensional phase space trajectory shows significant distortion, the phase offset modulus of each axis continuously exceeds the mean of the healthy state by ±3 times the standard deviation, or the phase deflection angle is greater than 1.5 radians, the equipment is determined to be in a failure state. Multiple sets (no less than 20 sets) of equipment failure state tests are conducted. For each set of tests, characteristic index vectors under failure states are collected. Using the same normalization and Mahalanobis distance calculation method as in step 301, the Mahalanobis distance between each failure state characteristic index vector and the standard state vector in the feature benchmark library is calculated one by one. It is assumed that the maximum Mahalanobis distance under all failure states is 4.545. This maximum value is increased by 10% (4.545 × 1.1 ≈ 5.0) as the threshold.) The final value is determined with a certain safety margin to ensure that the Mahalanobis distance does not exceed 5.0 under all fault conditions. After substituting this value into the health score formula, the health score under fault conditions is 0.

[0054] Step 302: Using the feature index vector as the center point in the feature space, a high-dimensional local neighborhood window is determined based on the center point. The Euclidean distance and cosine similarity between each standard state vector in the neighborhood and the center point are calculated. A two-dimensional joint histogram is constructed, and the histogram is expanded into a high-dimensional local feature descriptor, specifically including: The normalized feature index vector (9-dimensional) from step 301 is directly used as a point in the feature space, which has a dimension of 9. All standard state vectors in the feature benchmark library have also been normalized. First, the Euclidean distance between each pair of standard state vectors in the benchmark library is calculated, and 1.5 times the average of these distances is taken as the local neighborhood radius. Using the center point as the sphere's center and this radius as its length, a spherical neighborhood is constructed in the 9-dimensional feature space, and all standard state vectors falling within this neighborhood are selected. If the number of points in the neighborhood is less than 3, the radius is appropriately increased (by 10% each time) until the number of points is no less than 3, or the radius reaches a preset maximum value (such as twice the maximum pairwise distance in the benchmark library). These standard state vectors within the neighborhood constitute the local neighborhood point set. For each standard state vector in the local neighborhood point set, calculate the straight-line distance between the vector and the center point. The smaller the value, the closer the two points are. Calculate the cosine of the angle between the vector and the center point, and use the cosine value as the cosine similarity. The value ranges from [-1, 1]. The closer the cosine similarity is to 1, the more consistent the directions of the two vectors are; the closer it is to -1, the opposite the directions are; and the closer it is to 0, the orthogonal.

[0055] The Euclidean distance range [0, neighborhood radius] is uniformly divided into 8 intervals, and the cosine similarity range [-1, 1] is uniformly divided into 8 intervals, forming an 8×8 two-dimensional grid. For each point in the local neighborhood point set, the row and column indices are determined according to its Euclidean distance and cosine similarity, respectively, and the count in the corresponding grid cell is incremented by 1. After traversing all neighborhood points, 8×8=64 count values ​​are obtained. These count values ​​are expanded into a 64-dimensional vector in row-first, column-second order. This vector is the high-dimensional local feature descriptor. If the local neighborhood point set is empty (no standard state vector falls into the neighborhood), all 64 components of the descriptor are set to 0. To eliminate the amplitude difference caused by the number of points in the neighborhood, the 64-dimensional descriptor is L2 normalized, that is, the magnitude of the descriptor is calculated, and then each component is divided by the magnitude. If the magnitude is 0 (i.e., a zero descriptor), the zeros are kept unchanged. This descriptor uses Euclidean distance to reflect the overall proximity between neighboring points and the center point, and cosine similarity to reflect directional consistency. The joint distribution of these two factors can effectively characterize the clustering patterns and outlier features of standard state vectors within a local neighborhood. Specifically, if the count value of a cell exceeds 25% of the total number of points in the local neighborhood, it is determined that the cell exhibits significant clustering, indicating that multiple standard state vectors are highly consistent within the distance and direction intervals. If the count values ​​of two or more consecutive adjacent cells both exceed 15% of the total number of points in the neighborhood, it is determined to be a local clustering pattern. If the count values ​​of all cells are less than 10% of the total number of points in the neighborhood, it is determined to be a sparse distribution. If the count values ​​are mainly concentrated (more than 50% of the total number of points) in the cells with the 6th or higher distance sub-interval index (i.e., distance greater than 0.75 times the neighborhood radius) and the 2nd or lower similarity sub-interval index (i.e., similarity less than -0.5), it is determined that there are outliers in the neighborhood or that the current center point is severely deviated from the healthy state distribution area.

[0056] Step 303: Calculate the matching score between the high-dimensional local feature descriptor and the standard descriptors in the pre-built feature benchmark library. Based on the matching score, determine the deviation direction and magnitude of each component in the feature index vector relative to the feature benchmark library. If any deviation magnitude exceeds a preset warning threshold, then based on the component type and deviation direction exceeding the threshold, retrieve the predefined fault mapping table to determine the corresponding fault mode description. All determined fault mode descriptions are then combined into a fault diagnosis prompt set, specifically including: In addition to storing standard state vectors (9 dimensions) and their statistical parameters, the feature benchmark library also pre-stores standard descriptors corresponding to each standard state vector. The method for generating standard descriptors is completely consistent with step 302, that is, taking each standard state vector itself as the center point, constructing a local neighborhood on other standard state vectors in the benchmark library (or all standard state vectors excluding itself), calculating the Euclidean distance and cosine similarity between each point in the neighborhood and the center point, generating an 8×8 two-dimensional joint histogram, expanding it into a 64-dimensional vector and L2 normalizing it, generating 30 standard descriptors for 30 health state experiments. The high-dimensional local feature descriptor generated by the current device operating state (denoted as the current descriptor, 64 dimensions) is compared with the 30 standard descriptors one by one, and the maximum value among all similarities is taken as the matching score, with a value range from -1 to 1. The closer the matching score is to 1, the more similar the current descriptor is to the standard descriptor of a certain health state, that is, the closer the local neighborhood distribution of the current feature index vector is to the health state. Find the standard descriptor with the highest matching score to the current descriptor; its corresponding standard state vector is called the best matching standard vector. Compare the current feature index vector (9 components) with the corresponding components of the best matching standard vector one by one, and calculate the difference (current component - standard component). The sign of the difference indicates the direction of deviation (positive value indicates that the current component is higher than the standard value, negative value indicates that it is lower than the standard value). Divide the absolute value of the difference by the maximum absolute value of the component among all standard state vectors in the benchmark library (pre-statistically counted and stored) to obtain the normalized deviation magnitude, which ranges from [0, +∞). The larger the value, the more severe the deviation.

[0057] A warning threshold of 0.2 is set (i.e., a warning is triggered when the deviation exceeds 20%). The normalized deviation of all nine components is iterated to determine if any component exceeds 20%. If not, it indicates that the deviation of each component in the feature index vector is within the allowable range, and no fault indication is given. If any component has a deviation greater than 0.2, a predefined fault mapping table is retrieved based on the component's type (e.g., X-axis phase offset modulus, total phase offset modulus, phase deflection angle, multidimensional motion coupling index, trunk length, number of branches, and average curvature of branch points) and deviation direction (a positive deviation value indicates the current component is higher than the standard component, and a negative value indicates it is lower than the standard component). This fault mapping table is constructed based on the equipment's full lifecycle test data. The specific construction process is as follows: an equipment fault simulation test platform is built, combining clearly defined equipment design parameters, operating limit parameters, and historical fault data. The equipment design parameters include a spindle speed design value of 2000 rpm, a feed rate design value of 50 mm / min, and a key component fit clearance design value of 0.02 mm. Operating limits include spindle speed limits of 1500 to 2500 rpm, feed rate limits of 30 to 70 mm / min, clearance limits for key components of 0.01 to 0.05 mm, and lubrication system oil supply pressure limits of 0.3 to 0.8 MPa. Equipment failure records from the past three years include 42 cases of component wear failures, 38 cases of main vibration direction deviation failures, 29 cases of attitude divergence failures, and 35 cases of abnormal phase deflection angle failures. Each failure is documented with corresponding deviation data for characteristic parameters and abnormal manifestations.

[0058] Typical faults are simulated in the following four specific ways: adjusting the clearance of key components to the design limit of 0.05 mm to simulate component wear; adjusting the vibration excitation parameters to 1.2 times the design value to simulate deviation of the main vibration direction; changing the oil supply pressure of the lubrication system to the minimum value of 0.3 MPa to simulate attitude divergence caused by insufficient lubrication; and adjusting the phase control parameters to ±0.5 radians to simulate abnormal phase deflection angle. For each type of fault, the positive and negative deviations of the corresponding characteristic components are simulated separately. Real-time data of characteristic parameters, changes in three-dimensional phase space trajectory, and abnormal equipment operation feedback are collected simultaneously. The specific indicators for changes in three-dimensional phase space trajectory are: the maximum distance of trajectory deviation from the standard trajectory is ≥0.08 mm, and the degree of trajectory distortion is ≥30%. The specific indicators for abnormal equipment operation feedback are: the vibration angular velocity amplitude exceeds 1.5 times the average value under healthy conditions (for example, when the angular velocity amplitude under healthy conditions is 0.1 radians / second, the excess value is ≥0.15 radians / second), and the operating noise increases by more than 30 dB compared to the noise value under healthy conditions (when the noise value under healthy conditions is assumed to be 55 dB, the excess value is ≥85 dB).

[0059] Record the deviation magnitude, deviation direction, and corresponding abnormal behavior of characteristic parameters under each fault. Using component type and deviation direction as search keywords, organize the abnormal behaviors into standardized fault mode descriptions, establish a one-to-one correspondence, and store them to form a fault mapping table. For example, if the deviation magnitude of the X-axis phase offset modulus component is positive and exceeds the warning threshold (i.e., deviation magnitude > 20%), the corresponding fault mode description is an abnormal increase in attitude divergence in the X-axis direction; if the deviation magnitude of the phase deflection angle component is positive and exceeds the warning threshold, the corresponding fault mode description is an aggravated deviation between the main vibration direction and the time series reference. Collect and organize all retrieved fault mode descriptions to form a fault diagnosis prompt set, which contains fault mode descriptions corresponding to all components whose deviation magnitude exceeds the warning threshold (20%).

[0060] This embodiment ensures the reliability and consistency of the evaluation benchmark by calling a pre-built feature benchmark library. It converts the difference between the feature index vector and the standard state vector into an intuitive health score, accurately reflecting the degree of difference between the current operating state and the health state of the equipment, avoiding errors from subjective judgment. High-dimensional local feature descriptors can capture local feature patterns in the feature space. Combined with matching calculations with standard descriptors, it can not only determine the deviation magnitude of feature parameters but also clarify the direction of deviation. Then, by retrieving the corresponding fault mode through a fault mapping table, it achieves precise fault location. Simultaneously, by setting warning thresholds, it can promptly issue fault prompts when the equipment exhibits minor anomalies, achieving early warning of faults, preventing further deterioration of the fault, and reducing the risk of equipment damage.

[0061] In a preferred embodiment of the present invention, step 4 includes: Step 400: Preset a first health threshold, a second health threshold, and a third health threshold, wherein the first health threshold is greater than the second health threshold, and the second health threshold is greater than the third health threshold; when the real-time health score of the equipment is greater than or equal to the first health threshold, the equipment is determined to be in a stable operating state, and a control command to maintain the current process parameters is generated; when the real-time health score of the equipment is less than the first health threshold but greater than or equal to the second health threshold, the equipment is determined to be in a slightly abnormal state. If the fault diagnosis prompt set includes vibration-related fault modes, a fine-tuning control command is generated, the fine-tuning control command including... The feed rate is reduced or auxiliary lubrication is triggered according to a preset ratio; when the real-time health score of the equipment is less than the second health threshold and greater than or equal to the third health threshold, the equipment is determined to be in a significantly abnormal state. If the phase offset modulus shows a continuous increasing trend, an intervention control command is generated. The intervention control command includes performing a short pause to attenuate mechanical resonance or switching to a conservative machining parameter group; when the real-time health score of the equipment is less than the third health threshold, the equipment is determined to be in a serious fault state, and a protective control command is generated. The protective control command includes sending an emergency stop or pause command and triggering an alarm, specifically including: Three health thresholds are preset: the first health threshold (T1), the second health threshold (T2), and the third health threshold (T3). These thresholds are set based on equipment lifecycle test data, the boundary characteristics between healthy and abnormal states, and the actual operational safety, processing accuracy requirements, and fault handling costs in production. They were determined through multiple sets of comparative tests, following the order T1 > T2 > T3. The first health threshold, T1 = 85, is based on test data under healthy equipment conditions. Verified through 30 sets of health state tests, the health score of the equipment during healthy operation is consistently 85 or higher, and all characteristic parameters are within the standard range (X-axis, Y-axis, and Z-axis phase offset modulus are all 0.1 to 0.3, phase deflection angle is 0.2 to 0.8 radians, multidimensional motion coupling index is 0.4 to 0.6, and trunk length in geometric invariants is 0.5 to 0.7). The first health threshold has 3 to 5 branches and a branch point curvature change of 0.1 to 0.2. Without any abnormal feedback, setting T1=85 accurately determines the stable operating status of the equipment. The second health threshold T2=60 is based on boundary test data between minor and significant abnormalities. When the health score is between 60 and 84, the equipment only shows slight parameter deviations (deviation ≤ 0.3), with no serious risk of failure, and can be restored to normal through fine-tuning. Setting T2=60 distinguishes between minor and significant abnormalities. The third health threshold T3=30 is based on boundary test data between significant and serious failures. When the health score is between 30 and 59, the equipment shows multiple parameter exceedances, increasing the risk of failure, but this can still be mitigated through intervention. When the score is below 30, the equipment shows signs of serious failure, and continued operation will lead to component damage. Setting T3=30 accurately triggers protective control. All three thresholds range from [0, 100].

[0062] Obtain the real-time health score of the equipment calculated in step 301, and combine it with the fault diagnosis prompt set generated in step 303. Determine the equipment operating status according to the following logic and generate corresponding control commands. When the real-time health score of the equipment is greater than or equal to the first health threshold T1, the equipment is determined to be in a stable operating state. At this time, the difference between the current operating state and the health state of the equipment is minimal, all characteristic parameters are within the standard range, the fault diagnosis prompt set is empty, and there is no need to adjust the equipment operating parameters. Generate control commands to maintain the current process parameters. The core content of this control command is to keep the current process parameters such as feed rate, processing pressure, and lubrication frequency unchanged, continuously monitor the equipment operating status, collect data in real time and update the characteristic index vector to ensure stable equipment operation.

[0063] When the real-time health score of the device T2 ≤ the real-time health score of the device < T1, it is determined that the device is in a slightly abnormal state. At this time, the current operating state of the device shows a slight deviation, and the deviation range of some characteristic parameters is close to or exceeds the warning threshold (Twarn = 0.2). The fault diagnosis prompt set may include relevant fault modes. If the fault diagnosis prompt set includes vibration-related fault modes (such as deviation of the main vibration direction, abnormal vibration amplitude, etc.), a fine-tuning control instruction will be generated. The specific fine-tuning operation is as follows: If it is prompted that the deviation of the main vibration direction from the timing reference is increasing, the feed speed will be reduced by a fixed ratio of 15% (fine-tuned from the designed value of 50 mm / min to 42.5 mm / min), and at the same time, the lubrication frequency will be increased from the original 20 times / minute to 25 times / minute to reduce mechanical vibration and correct the main vibration direction; If it is prompted that the abnormal divergence of the X-axis direction attitude is increasing, the X-axis phase control parameter will be fine-tuned by 0.1 radian (approaching the healthy standard range of 0.2 to 0.8 rad), and at the same time, the processing pressure will be reduced by 10% (fine-tuned from the original designed value of 0.5 MPa to 0.45 MPa) to relieve the attitude divergence; If multiple vibration-related faults occur simultaneously, the above fine-tuning operations can be combined and executed, and the feed speed and lubrication parameters are adjusted first; If there is no vibration-related fault mode in the fault diagnosis prompt set, only a status monitoring instruction will be generated, continuously tracking the parameter changes, shortening the monitoring period from the original 10 s to 5 s, and updating the characteristic index vector in real time to ensure that the abnormality does not deteriorate further.

[0064] When the real-time health score of the device T3 ≤ the real-time health score of the device < T2, it is determined that the device is in a significant abnormal state. At this time, the current operating state of the device deviates significantly from the healthy state, and the deviation amplitudes of multiple characteristic parameters exceed the warning threshold, with a relatively high failure risk. If the phase shift modulus of each axis calculated in step 200c shows a continuous increasing trend (determined by comparing the phase shift modulus data of three consecutive adjacent fixed analysis time windows, and the length of each analysis time window is the same as that in step 200a, i.e., 1 second per window; the preset increase threshold is 0.05, which is set based on the standard range of the phase shift modulus (0.1 to 0.3) and is 1 / 6 of the maximum value of the standard range; if the phase shift modulus of the latter window is greater than that of the previous window, and the difference between adjacent windows exceeds 0.05, it is determined to be a continuous increasing trend), an intervention control instruction is generated; if there is no continuous increasing trend in the phase shift modulus, a strengthened monitoring instruction is generated to shorten the monitoring cycle and track the development trend of the failure in real time. Among them, the intervention control instruction includes two optional operations. One is to execute a short-term pause to attenuate mechanical resonance. The pause time is preset to be 5 to 10 seconds. During the pause, the feed mechanism is closed, and other systems operate normally. The mechanical resonance is attenuated using the pause time to relieve the abnormal vibration of the device. The other is to switch to a conservative processing parameter group, which is preset as a parameter combination with a 30% to 50% reduction in the feed speed and a 20% to 30% reduction in the processing pressure. By reducing the processing load, the loss of device components is reduced, and the abnormal state is prevented from deteriorating further.

[0065] When the real-time health score of the device < T3, it is determined that the device is in a serious failure state. At this time, the current operating state of the device seriously deviates from the healthy state, with a major failure risk. Continuing to operate may cause damage to device components, a significant decline in processing accuracy, and even lead to safety accidents. Therefore, a protective control instruction is generated. This control instruction includes two core operations. One is to send an emergency stop or pause instruction to immediately stop the processing operation of the device and close key mechanisms such as the feed and spindle to prevent the failure from expanding further. The other is to trigger an alarm, send an alarm signal through the device's audible and visual alarm system, and at the same time transmit the failure information to the monitoring terminal to remind the operator to handle the failure in a timely manner and investigate the cause of the device abnormality. When generating various control instructions, it is necessary to ensure the accuracy and executability of the instruction parameters. All control instructions must include specific operation parameters (such as the adjustment ratio of the feed speed, pause time, etc.), be compatible with the device's control system, and can be directly executed by the device. At the same time, record the generation time of the control instruction, the corresponding health score, and the failure prompt for subsequent failure traceability and parameter optimization.

[0066] Step 401, encapsulate the control instructions for maintaining the current process parameters, fine-tuning control instructions, intervention control instructions, and protective control instructions into a control instruction set, specifically including: Collect all control commands generated in step 400, including control commands to maintain current process parameters, fine-tuning control commands, intervention control commands, and protective control commands. Each control command contains complete operation content, execution parameters, and triggering conditions. The collected control commands are categorized and sorted according to the priority of equipment operating status (stable operating status, minor abnormal status, significant abnormal status, and severe fault status) to ensure the logical consistency of the control command set. This facilitates the control system's rapid recall of the corresponding control command based on the current equipment status. Simultaneously, each control command is standardized and encapsulated, unifying the command format and clearly defining the command's identifier, execution target, operation parameters, execution time, and termination conditions.

[0067] For example, the encapsulation format of the control instruction to maintain the current process parameters is as follows: Instruction identifier - maintenance parameter, execution object - equipment feed system, spindle system, operation parameters - current feed speed, machining pressure, lubrication frequency, execution time - continuous execution, termination condition - equipment status changes; the encapsulation format of the fine-tuning control instruction is as follows: Instruction identifier - fine-tuning control, execution object - equipment feed system, lubrication system, operation parameters - feed speed reduction ratio, lubrication start time, execution time - until the equipment status returns to stability, termination condition - health score ≥ T1 or fault indication disappears.

[0068] After encapsulation, all standardized control commands are combined to form a control command set. This set contains control commands corresponding to all possible operating states of the equipment, and an command index is added to facilitate the control system's quick retrieval and execution of corresponding control commands based on the equipment's real-time health score and fault diagnosis prompts. The control command set is then validated to ensure correct command format, complete parameters, and logical coherence. Once validated, it is transmitted to the equipment control system.

[0069] This embodiment divides the equipment operating status into four levels based on the health threshold, and generates corresponding control commands for different states. The stable operating state maintains the current parameters to ensure processing efficiency; the minor abnormal state is fine-tuned to alleviate the abnormality and avoid excessive efficiency decline; the abnormal state is intervened to prevent the fault from worsening; and the severe fault state is protected to avoid equipment damage and safety accidents. This achieves a balance between efficiency and safety and solves the problem of a single control method that cannot be flexibly adjusted according to the state.

[0070] In a preferred embodiment of the present invention, step 5 includes: Step 500: The control command set is sent to the machining equipment control system via industrial communication protocol, and the control system executes the corresponding control commands in the control command set. After executing the control commands, a new round of original three-axis angular velocity timing data is collected, and a new feature index vector and a new real-time equipment health score are obtained after executing the control commands, specifically including: All control commands generated in step 400 based on the four-level health interval classification are summarized and integrated into a complete control command set. After unified format encapsulation and verification, the control command set is sent to the machining equipment control system using a standard industrial communication protocol. Upon receiving the command set, the equipment control system matches and executes the corresponding category of control commands according to the operating status levels defined in step 400. After the control commands are executed, a new round of original three-axis angular velocity time-series data is collected after a stable buffer period. The previously unified min-max normalization method is used to unify the dimensions of all feature parameters, eliminating dimensional differences between different physical indicators. Following the complete feature extraction, vector construction, and Mahalanobis distance calculation process described in steps 200 to 302, the calculations are performed sequentially to obtain a new feature index vector after the control command execution. Finally, the updated real-time health score of the equipment is calculated.

[0071] Step 501: If the new real-time health score of the device evolves towards an increasing trend compared to the real-time health score of the device before executing the control command, and the multidimensional motion coupling index in the new feature index vector shows a converging trend compared to the multidimensional motion coupling index before executing the control command, then the control strategy is confirmed to be effective, and the statistical parameters in the feature benchmark library are updated according to the new feature index vector and the new real-time health score of the device, specifically including: Extract core quantitative data before and after the execution of control commands: the data before execution (health score, feature index vector, multidimensional motion coupling index, Mahalanobis distance) comes from steps 301 and 200; the data after execution comes from the re-collected calculation results in step 500. Calculate the deviation correction difference of each feature component dimensionally (value after execution minus value before execution); simultaneously calculate three core constraint indicators, namely, health improvement = health score after execution - health score before execution (greater than 0 indicates improvement); multidimensional motion coupling index convergence deviation = |Cnew - Ccenter| - |Cold - Ccenter|, where Ccenter is the center value of the health standard interval, and according to the standard interval [0.4, 0.6] defined in step 201d, Ccenter = 0.5; Cold is the coupling index before execution, and Cnew is the coupling index after execution; convergence deviation < 0 indicates convergence towards the standard interval; Mahalanobis distance reduction = Mahalanobis distance before execution - Mahalanobis distance after execution (greater than 0 indicates a reduction in the difference between the feature vector and the benchmark library).

[0072] Based on the above three indicators, the Exponentially Weighted Moving Average (EWMA) method is used to update the statistical parameters of the corresponding working conditions in the feature benchmark library, that is, to update the mean vector. ,in =0.1 is the same smoothing factor as in the mean update. It is the health status mean vector stored in the feature benchmark library before the update. It is a 9-dimensional vector, and each component corresponds to the arithmetic mean of the feature index vector (X-axis phase offset modulus, Y-axis phase offset modulus, Z-axis phase offset modulus, total phase offset modulus, phase deflection angle, multidimensional motion coupling index, trunk length, number of branches, and mean curvature of branch points) in the healthy state. The new feature index vector after executing the instruction; update the covariance matrix. ,in It is the health status covariance matrix stored in the feature benchmark library before the update, which is a 9×9 symmetric positive definite matrix. It is the transpose symbol; the updated one and It is immediately used for subsequent calculations of Mahalanobis distance and health scores.

[0073] Step 502: If the new device real-time health score does not increase or the multidimensional motion coupling index does not converge, the current control level is upgraded by one level according to the preset upgrade rules, an upgraded control command is generated and executed again. If there is still no improvement after a preset number of consecutive attempts, a manual intervention alarm is triggered, forming a closed-loop adaptive monitoring process, specifically including: Comparing the two core judgment indicators before and after the command execution, if the improvement in the real-time health score of the new equipment is less than 5 (this threshold is based on the health score range of 0 to 100, with 5 being the minimum meaningful improvement; values ​​below this are considered no significant improvement), does not show a positive increase, or the multidimensional motion coupling index does not show a convergence correction trend towards the 0.4 to 0.6 standard range defined in step 200, and any one of these conditions is not met, the control strategy for this round is deemed to have failed. The control level is adjusted according to a clear step-by-step upgrade rule: a minor anomaly control level corresponds to a health score of 60 to 85, and after the upgrade, it switches to a significant anomaly control level; a significant anomaly control level corresponds to a health score of 30 to 60, and after the upgrade, it switches to a severe fault control level. After the level upgrade, referring to the higher-level anomaly control logic recorded in step 400, a control command adapted to the current level is generated and reissued for execution. The maximum number of consecutive failures for improvement is fixed at 3 times. This threshold is uniformly set by the system. If, after 3 consecutive level upgrade adjustments, the equipment health score and various characteristic indicators still do not meet the improvement requirements, a manual intervention alarm is immediately triggered. Based on the entire process of issuing and executing instructions, repeatedly collecting operational status data, quantitatively verifying control effects, progressively upgrading control levels, and triggering alarms for extreme anomalies, a complete closed-loop adaptive monitoring and control process is formed.

[0074] This embodiment categorizes and prioritizes control commands based on the actual operating status of the equipment. Combined with a command index-based rapid retrieval mechanism, it can quickly match corresponding control schemes by integrating real-time health scores and fault diagnosis results, improving the speed of equipment anomaly response and the efficiency of control execution. A tiered control level escalation mechanism and manual intervention alarm thresholds are established, forming a graded handling mode of hierarchical control, step-by-step intervention, and extreme alarms. This avoids the continuous deterioration of equipment anomalies due to the failure of a single control method, effectively reducing component wear and the risk of fault escalation.

[0075] like Figure 2 As shown, embodiments of the present invention also provide a system for monitoring the operation of machining equipment using a three-axis gyroscope, comprising: The preprocessing module is used to collect raw triaxial angular velocity time series data and perform preprocessing to obtain purified triaxial angular velocity time series data; The calculation module is used to reconstruct the phase space of the purification triaxial angular velocity time series data to obtain a three-dimensional phase space trajectory. Within the three-dimensional phase space trajectory, the first phase pole, the second phase pole, and the third reference node are determined. The central axis skeleton is extracted from the interior of the three-dimensional phase space trajectory through distance transformation, and the main trunk length, number of branches, and curvature of the branch points of the central axis skeleton are statistically analyzed as geometric invariant features of the three-dimensional phase space trajectory morphology. Based on the first phase pole, the second phase pole, and the third reference node, the phase offset modulus components and phase deflection angles of each axis are determined. The cross-correlation matrix between the phase offset modulus components of different axes is calculated to obtain the multidimensional motion coupling index. The phase offset modulus components of each axis, the phase deflection angle, the multidimensional motion coupling index, and the geometric invariant features are combined to obtain a feature index vector. The evaluation module is used to obtain the real-time health score of the equipment based on the comparison between the feature index vector and the pre-built feature benchmark library, and to monitor the deviation direction and magnitude of each component in the feature index vector relative to the pre-built feature benchmark library, and determine the set of fault diagnosis prompts. The decision-making module is used to make hierarchical decisions based on the real-time health score and fault diagnosis prompts of the equipment, and to obtain a set of control commands. The monitoring module is used to execute a set of control commands and verify their effects, forming a closed-loop adaptive monitoring process.

[0076] It should be noted that this system is a system corresponding to the above method. All implementation methods in the above method embodiments are applicable to this embodiment and can achieve the same technical effect.

[0077] The above description represents the preferred embodiments of the present invention. It should be noted that those skilled in the art can make various improvements and modifications without departing from the principles of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A method for monitoring the operation of machining equipment using a three-axis gyroscope, characterized in that, The method includes: Step 1: Collect raw triaxial angular velocity time series data and preprocess it to obtain purified triaxial angular velocity time series data; Step 2: Reconstruct the phase space of the three-axis angular velocity time series data to obtain a three-dimensional phase space trajectory. Within the three-dimensional phase space trajectory, determine the first phase pole, the second phase pole, and the third reference node. Extract the central axis skeleton from the interior of the three-dimensional phase space trajectory through distance transformation, and statistically analyze the trunk length, number of branches, and curvature of the branch points of the central axis skeleton as geometric invariant features of the three-dimensional phase space trajectory morphology. Determine the phase offset modulus components and phase deflection angles of each axis based on the first phase pole, the second phase pole, and the third reference node. Calculate the cross-correlation matrix between the phase offset modulus components of different axes to obtain the multidimensional motion coupling index. Combine the phase offset modulus components of each axis, the phase deflection angle, the multidimensional motion coupling index, and the geometric invariant features to obtain a feature index vector. Step 3: Based on the comparison between the feature index vector and the pre-built feature benchmark library, obtain the real-time health score of the device, and monitor the deviation direction and magnitude of each component in the feature index vector relative to the pre-built feature benchmark library to determine the fault diagnosis prompt set. Step 4: Based on the real-time health score and fault diagnosis prompts of the equipment, a hierarchical decision is made to obtain a set of control commands; Step 5: Execute the set of control commands and verify the effect to form a closed-loop adaptive monitoring process.

2. The method for monitoring the operation of machining equipment using a three-axis gyroscope according to claim 1, characterized in that, Step 1 includes: A three-axis gyroscope sensor is rigidly fixed to the reciprocating motion component of the machining equipment, and instantaneous angular velocity data of the reciprocating motion component in three-dimensional space are synchronously collected at a preset sampling frequency to form raw three-axis angular velocity time series data. The raw triaxial angular velocity time series data is preprocessed to obtain purified triaxial angular velocity time series data.

3. The method for monitoring the operation of machining equipment using a three-axis gyroscope according to claim 2, characterized in that, Phase space reconstruction was performed on the time-series data of the purification triaxial angular velocity to obtain the three-dimensional phase space trajectory, including: The three-axis angular velocity component sequence within a fixed analysis time window is extracted from the purified three-axis angular velocity time series data. The X-axis, Y-axis and Z-axis angular velocity components are used as the three coordinate axes of the three-dimensional rectangular coordinate system. The data points at each moment in the three-axis angular velocity component sequence are mapped to the original coordinate points in three-dimensional space to form the original point set. Based on the principal eigenvalues ​​of the covariance matrix between each data point in the original point set and the remaining data points in the neighborhood of each data point, the mutual information in the time neighborhood, and the Euclidean distance gradient, the local clustering metric and the local discreteness metric of each data point are determined. Based on the local clustering metric and local discrete metric of each data point, the coordinates of each data point in the original point set are adaptively adjusted. Points with a local clustering metric greater than the neighborhood mean are given a shrinking offset along the direction of the corresponding point's neighborhood center, and points with a local discrete metric greater than the neighborhood mean are given an expanding offset along the direction of the corresponding point's neighborhood outside. The shrinking offset and expanding offset are alternately executed until the coordinate change of all data points is less than the preset convergence threshold, thus obtaining the adjusted point set. The data points in the adjusted point set are connected sequentially with smooth curves according to the original time order to obtain a first-order continuous and differentiable three-dimensional spatial trajectory, which is defined as a three-dimensional phase space trajectory.

4. The method for monitoring the operation of machining equipment using a three-axis gyroscope according to claim 3, characterized in that, In the three-dimensional phase space trajectory, the first phase pole, the second phase pole, and the third reference node are determined respectively. The central skeleton is extracted from the interior of the three-dimensional phase space trajectory through distance transformation. The trunk length, number of branches, and curvature of the branch points of the central skeleton are statistically analyzed as geometric invariant features of the three-dimensional phase space trajectory morphology, including: Calculate the local curvature value of each data point on the three-dimensional phase space trajectory, and select the top N data points with the largest local curvature values ​​as the candidate pole set; in the candidate pole set, select the data point with the largest curvature value and the earliest corresponding time, and determine it as the first phase pole representing the abrupt change in the geometric shape of the motion trajectory. The energy density distribution of the three-dimensional phase space trajectory within the analysis time window is calculated, the position of the energy density centroid is determined, and the trajectory data point closest to the energy density centroid is determined as the second phase pole characterizing the centroid of the motion energy distribution. In the three-dimensional phase space trajectory, the trajectory data point corresponding to the middle moment of the analysis time window is selected as the third reference node for calibrating the motion state time series reference. After the three poles are determined, the three-dimensional phase space trajectory is discretized into a three-dimensional grid point set. The shortest Euclidean distance from each grid point inside the three-dimensional phase space trajectory to the boundary of the three-dimensional phase space trajectory is calculated to form a range field. Wavefronts propagating simultaneously from the boundary of the three-dimensional phase space trajectory are simulated. The positions in the range field that have local maxima and where the wavefronts first meet are marked as the central axis points. All central axis points constitute the central axis point set. Connect adjacent central axis points and determine the connection path based on the gradient direction of the distance field to construct a skeleton graph with a tree-like topology. For the tree-like skeleton structure, calculate the path length from the root node to each leaf node, and take the length of the longest path as the trunk length. Count the number of all branch nodes in the skeleton graph as the number of branches. At each branch node, fit the direction vectors of the two adjacent skeleton segments of the corresponding node, calculate the angle between the direction vectors of the two skeleton segments, and define the angle as the curvature change at the corresponding branch point. Calculate the arithmetic mean of the curvature changes at all branch points, denoted as the average branch point curvature. Use the trunk length, the number of branches, and the average branch point curvature together as geometric invariant features of the trajectory morphology in three-dimensional phase space.

5. The method for monitoring the operation of machining equipment using a three-axis gyroscope according to claim 4, characterized in that, The phase offset modulus components and phase deflection angles for each axis are determined based on the first phase pole, the second phase pole, and the third reference node, including: Calculate the absolute values ​​of the coordinate differences between the first phase pole and the second phase pole along the X-axis, Y-axis, and Z-axis, and define the three absolute values ​​as the X-axis phase offset modulus, Y-axis phase offset modulus, and Z-axis phase offset modulus, respectively. At the same time, calculate the Euclidean distance between the first phase pole and the second phase pole in the three-dimensional rectangular coordinate system, and define it as the total phase offset modulus, which is used to quantify the degree of attitude divergence of the moving parts at the microscale. Construct a first vector pointing from the first phase pole to the second phase pole, and a second vector pointing from the second phase pole to the third reference node; calculate the angle between the first vector and the second vector, and define the angle as the phase deflection angle, which is used to quantify the degree of deviation of the main vibration direction of the moving part from the timing reference.

6. The method for monitoring the operation of machining equipment using a three-axis gyroscope according to claim 5, characterized in that, Calculate the cross-correlation matrix between phase offset modulus components of different axes to obtain the multidimensional motion coupling index; combine the phase offset modulus components of each axis, phase deflection angle, multidimensional motion coupling index, and geometric invariant features to obtain the feature index vector, including: Extract the time series of X-axis phase offset modulus, Y-axis phase offset modulus, and Z-axis phase offset modulus respectively; calculate the Pearson cross-correlation coefficients between each pair of X-axis, Y-axis, and Z-axis phase offset modulus time series, and construct the cross-correlation coefficient matrix; The mean of the off-diagonal elements in the cross-correlation matrix is ​​calculated as a multidimensional motion coupling index that characterizes the tightness of topological constraints of multi-source sensor data in the feature space. The X-axis phase offset modulus, Y-axis phase offset modulus, Z-axis phase offset modulus, total phase offset modulus, phase deflection angle, multidimensional motion coupling index, and geometric invariant features are arranged in a preset order and combined to generate a feature index vector.

7. The method for monitoring the operation of machining equipment using a three-axis gyroscope according to claim 6, characterized in that, Step 3 includes: Call the pre-built feature benchmark library, which stores the standard state vectors and statistical distribution parameters of the device when it is running in a healthy state; Calculate the Mahalanobis distance between the current feature index vector and the standard state vector in the feature benchmark library to obtain the real-time health score of the device. The feature index vector is used as the center point in the feature space. The center point is used to determine the high-dimensional local neighborhood window. The Euclidean distance and cosine similarity between each standard state vector in the neighborhood and the center point are calculated. A two-dimensional joint histogram is constructed and the histogram is expanded into a high-dimensional local feature descriptor. Calculate the matching score between the high-dimensional local feature descriptor and the standard descriptor in the pre-built feature benchmark library. Based on the matching score, determine the deviation direction and magnitude of each component in the feature index vector relative to the feature benchmark library. If any deviation magnitude exceeds a preset warning threshold, then based on the component type and deviation direction that exceeds the threshold, retrieve the predefined fault mapping table, determine the corresponding fault mode description, and assemble all determined fault mode descriptions into a fault diagnosis prompt set.

8. The method for monitoring the operation of machining equipment using a three-axis gyroscope according to claim 7, characterized in that, Step 4 includes: The system presets a first health threshold, a second health threshold, and a third health threshold, wherein the first health threshold is greater than the second health threshold, and the second health threshold is greater than the third health threshold. When the real-time health score of the equipment is greater than or equal to the first health threshold, the equipment is determined to be in a stable operating state, and a control command to maintain the current process parameters is generated. When the real-time health score of the equipment is less than the first health threshold but greater than or equal to the second health threshold, the equipment is determined to be in a slightly abnormal state. If the fault diagnosis prompt set includes vibration-related fault modes, a fine-tuning control command is generated, which includes reducing the feed speed according to a preset ratio or triggering auxiliary lubrication. When the real-time health score of the equipment is less than the second health threshold but greater than or equal to the third health threshold, the equipment is determined to be in a significantly abnormal state. If the phase offset modulus shows a continuous increasing trend, an intervention control command is generated, which includes performing a short pause to attenuate mechanical resonance or switching to a conservative processing parameter group. When the real-time health score of the equipment is less than the third health threshold, the equipment is determined to be in a serious fault state, and a protective control command is generated, which includes sending an emergency stop or pause command and triggering an alarm. The control commands for maintaining current process parameters, fine-tuning control commands, intervention control commands, and protective control commands are encapsulated into a set of control commands.

9. The method for monitoring the operation of machining equipment using a three-axis gyroscope according to claim 8, characterized in that, Step 5 includes: The control command set is sent to the control system of the machining equipment through the industrial communication protocol. The control system executes the corresponding control command in the control command set. After the control command is executed, a new round of original three-axis angular velocity timing data is collected, and a new feature index vector and a new real-time health score of the equipment are obtained after the execution of the control command. If the new real-time health score of the device evolves towards an increasing trend compared to the real-time health score of the device before the execution of the control command, and the multidimensional motion coupling index in the new feature index vector shows a convergent trend compared to the multidimensional motion coupling index before the execution of the control command, then the control strategy is confirmed to be effective, and the statistical parameters in the feature benchmark library are updated according to the new feature index vector and the new real-time health score of the device. If the real-time health score of the new device does not increase or the multidimensional motion coupling index does not converge, the current control level will be upgraded by one level according to the preset upgrade rules, the upgraded control command will be generated and executed again. If there is still no improvement after a preset number of consecutive times, a manual intervention alarm will be triggered, forming a closed-loop adaptive monitoring process.

10. A system for monitoring the operation of machining equipment using a three-axis gyroscope, the system implementing the method as described in any one of claims 1 to 9, characterized in that, include: The preprocessing module is used to collect raw triaxial angular velocity time series data and perform preprocessing to obtain purified triaxial angular velocity time series data; The calculation module is used to reconstruct the phase space of the purification triaxial angular velocity time series data to obtain the three-dimensional phase space trajectory; In the three-dimensional phase space trajectory, the first phase pole, the second phase pole, and the third reference node are determined respectively. The central skeleton is extracted from the interior of the three-dimensional phase space trajectory through distance transformation. The trunk length, branch number, and branch point curvature of the central skeleton are statistically analyzed as geometric invariant features of the three-dimensional phase space trajectory morphology. The phase offset modulus components and phase deflection angles of each axis are determined based on the first phase pole, the second phase pole, and the third reference node; the cross-correlation matrix between the phase offset modulus components of different axes is calculated to obtain the multidimensional motion coupling index; the phase offset modulus components of each axis, the phase deflection angle, the multidimensional motion coupling index, and the geometric invariant features are combined to obtain the feature index vector. The evaluation module is used to obtain the real-time health score of the equipment based on the comparison between the feature index vector and the pre-built feature benchmark library, and to monitor the deviation direction and magnitude of each component in the feature index vector relative to the pre-built feature benchmark library, and determine the set of fault diagnosis prompts. The decision-making module is used to make hierarchical decisions based on the real-time health score and fault diagnosis prompts of the equipment, and to obtain a set of control commands. The monitoring module is used to execute a set of control commands and verify their effects, forming a closed-loop adaptive monitoring process.

Images