Water supply pump shaft center locus dynamic prediction system based on eddy current displacement signal

Abstract

This invention relates to the field of industrial monitoring technology and discloses a dynamic prediction system for the shaft center trajectory of a water pump based on eddy current displacement signals. The system addresses the problem that conventional methods often mistakenly omit drift and underestimate runout amplitude due to the initial attenuation of the recentering force caused by the internal throttling seal of the water pump, followed by the expansion of the loop due to the thermal effect of the sliding bearing. This is because conventional methods often incorrectly omit drift amounts. The system utilizes dual probes to acquire the original displacement sequence, extracts the dominant period to divide the trajectory window, and performs smoothing. It extracts the shaft center and loop area separately, calculates the maximum time-delay correlation between the center migration intensity and the area growth rate to construct a hysteresis index, extrapolates the future center trajectory using this index, and maps the delayed center migration intensity to the shape parameters of the future loop. Finally, it generates the future boundary envelope and calculates the radial distance difference between the farthest and nearest points as the maximum radial runout amplitude output. This invention effectively preserves the time-staggered evolution information and improves the accuracy of early warning.

Description

Technical Field

[0001] This invention relates to the field of industrial monitoring technology, and more specifically, to a dynamic prediction system for the shaft trajectory of a water pump based on eddy current displacement signals. Background Technology

[0002] The feedwater pump rotor operates in a multi-stage controlled clearance seal, balancing drum, and sliding bearing complex environment, resulting in static deflection typically greater than that of a typical single-stage rigid pump. During actual operation, the recentering effect of throttling seals such as interstage bushings and balancing devices changes with pressure differential, clearance, and leakage, initially altering the shaft's average equilibrium position within the bearing clearance. Long-term cumulative factors such as sliding bearing wear, oil film thermal state, and oil supply variations then significantly impact the local loop area, major and minor axes, and stability of the trajectory after a certain period. Existing dynamic monitoring schemes for shaft trajectory often rely on general image recognition or direct sequence prediction, failing to adequately consider the staggered evolutionary relationship between the initial average position drift and subsequent trajectory loop expansion. If the slow average position migration is treated as ordinary baseline drift and removed during preprocessing, subsequent prediction extrapolation utilizes only local loop information, ultimately leading to a systematic underestimation of the future trajectory boundary and maximum radial runout prediction within the preset time window. Summary of the Invention

[0003] This invention provides a dynamic prediction system for the shaft trajectory of a water pump based on eddy current displacement signals, which solves the technical problems mentioned in the background art.

[0004] This invention provides a dynamic prediction system for the shaft center trajectory of a water pump based on eddy current displacement signals, comprising: The signal preprocessing module is used to acquire the original displacement sequence through dual probes, extract the dominant period using the synthetic autocorrelation function to divide multiple trajectory windows, and smooth the original displacement sequence of each trajectory window. The parameter extraction module is used to extract the axis center and loop area of ​​the smoothed displacement sequence within the trajectory window, calculate the center migration intensity based on the axis center of the adjacent trajectory windows, calculate the area growth rate based on the loop area of ​​the adjacent trajectory windows, and calculate the maximum time lag correlation between the center migration intensity and the area growth rate to construct a lag index. The prediction extrapolation module is used to extrapolate the future center trajectory using the lag exponent, and after delaying the center migration intensity according to the maximum time lag correlation, it is mapped to the shape parameters of the future hysteresis loop. The result output module is used to generate a future boundary envelope using the future center trajectory and the shape parameters. Taking the bearing clearance center as the origin, it calculates the radial distance difference between the farthest point and the nearest point of the future boundary envelope and outputs it as the maximum radial runout amplitude.

[0005] The beneficial effects of this invention are as follows: Addressing the complex time-varying evolution phenomenon caused by the different environments of the multi-stage controlled clearance seals and sliding bearings inside the water pump, this invention decouples the slow migration of the shaft center position from the expansion of the local loop. By constructing a hysteresis index, the sequential driving relationship between the center shift and the loop expansion is determined. During denoising and prediction extrapolation, the true center offset, which is usually mistakenly deleted as baseline drift, is fully preserved. The delayed center migration intensity is directly mapped to the shape parameters of the future loop. This overcomes the technical defect of traditional single-sequence prediction, which easily and systematically underestimates the maximum radial runout boundary, significantly improving the accuracy and reliability of dynamic prediction for complex rotating machinery under transitional operating conditions or thermal changes. Attached Figure Description

[0006] Figure 1 This is a block diagram of the dynamic prediction system for the shaft trajectory of a water pump based on eddy current displacement signals according to the present invention. Detailed Implementation

[0007] The subject matter described herein will now be discussed with reference to exemplary embodiments. It should be understood that these embodiments are discussed only to enable those skilled in the art to better understand and implement the subject matter described herein, and changes may be made to the function and arrangement of the elements discussed without departing from the scope of this specification. Various processes or components may be omitted, substituted, or added as needed in the examples. Furthermore, features described in some examples may be combined in other examples.

[0008] like Figure 1 As shown, the dynamic prediction system for the shaft center trajectory of a water pump based on eddy current displacement signals includes: The signal preprocessing module is used to acquire the original displacement sequence through dual probes, extract the dominant period using the synthetic autocorrelation function to divide multiple trajectory windows, and smooth the original displacement sequence of each trajectory window. The parameter extraction module is used to extract the axis center and loop area of ​​the smoothed displacement sequence within the trajectory window, calculate the center migration intensity based on the axis center of the adjacent trajectory windows, calculate the area growth rate based on the loop area of ​​the adjacent trajectory windows, and calculate the maximum time lag correlation between the center migration intensity and the area growth rate to construct a lag index. The prediction extrapolation module is used to extrapolate the future center trajectory using the lag exponent, and after delaying the center migration intensity according to the maximum time lag correlation, it is mapped to the shape parameters of the future hysteresis loop. The result output module is used to generate a future boundary envelope using the future center trajectory and the shape parameters. Taking the bearing clearance center as the origin, it calculates the radial distance difference between the farthest point and the nearest point of the future boundary envelope and outputs it as the maximum radial runout amplitude.

[0009] Preferably, the signal preprocessing module acquires the original displacement sequence through dual probes, extracts the dominant period using a synthetic autocorrelation function to divide the trajectory into multiple trajectory windows, and performs smoothing processing on the original displacement sequence of each trajectory window, including: Obtain the original displacement sequence of the two channels. and Calculate the synthesis autocorrelation function of the original displacement sequence. :

[0010] When calculating the synthetic autocorrelation function, the original displacement sequence of the two channels must first undergo dimensionless preprocessing with a fixed reference. The fixed reference used in the preprocessing is the full-scale displacement calibration value of the eddy current displacement probe of the corresponding channel, which remains constant throughout the entire operating cycle of the equipment. Each sample value of the original displacement sequence is divided by the full-scale calibration value of its respective channel to obtain a dimensionless displacement relative ratio sequence. Then, the synthetic autocorrelation function is calculated so that the function output is a dimensionless normalized correlation value.

[0011] The sample length corresponding to the first positive local maximum of the synthetic autocorrelation function is taken as the dominant period. :

[0012] According to the dominant cycle The original displacement sequence is sliced ​​over time to obtain the multiple trajectory windows. And according to the preset prediction time With sampling interval Calculate the number of future windows that need to be extrapolated :

[0013] For the original displacement sequence within the trajectory window Construct a smooth optimization model with second-order difference constraints:

[0014] Adaptive calculation of optimal smoothing intensity using generalized cross-validation method :

[0015] According to the optimal smoothing strength The original displacement sequence is solved to obtain the smoothed displacement sequence. ; in, The total number of samples, This is the sample lag. For sample serial number, This is the starting sample number of the trajectory window. To smooth the sequence during the optimization process. For smooth strength, This is the optimal smoothing sequence for the current channel. The corresponding smoothing matrix, and These are the optimal smoothing sequences for two channels, For the trace operation of a matrix, It is a 2-norm.

[0016] The original displacement sequence xn of the dual-channel probe is a discrete sequence of axial displacement relative to the probe output by the first channel at each sampling time. It can be continuously acquired by the first eddy current displacement probe installed at the bearing of the water pump.

[0017] The original displacement sequence yn of the dual-channel probe is a discrete sequence of axis-relative displacements output by the second channel at each sampling time. It can be continuously acquired by a second eddy current displacement probe installed at the bearing of the water pump and orthogonal to the first probe.

[0018] The synthetic autocorrelation function is a periodic evaluation function obtained by jointly normalizing the similarity of the original displacement sequences of the two channels under different sample hysteresis.

[0019] The sample lag is the number of samples that the current comparison sample is shifted backward relative to the reference sample when calculating the synthetic autocorrelation function.

[0020] The dominant period is the length of a locally repeating trajectory corresponding to the first positive local maximum of the synthetic autocorrelation function.

[0021] The trajectory window is a set of local trajectory samples obtained by segmenting the original displacement sequence according to the dominant period.

[0022] The starting sample number of a trajectory window is the first sample number of a trajectory window in the original displacement sequence.

[0023] The preset prediction time is the period for observing future risks set in advance by the system. It is preferably the time length corresponding to 2 to 5 dominant cycles, because it needs to cover the successive evolution process of the axis center migration and the lap area expansion, while avoiding the cumulative error caused by too many prediction steps.

[0024] The sampling interval is the time interval between two adjacent sampling points. It is preferable to include at least 128 sampling points in each dominant period, because it is necessary to simultaneously resolve the slowly varying axis center and local loop details.

[0025] The number of future windows is the number of future extrapolated windows calculated from the preset prediction time, dominant period, and sampling interval.

[0026] The original displacement sequence zn is one of the original displacement sequences that participates in the current channel smoothing optimization calculation.

[0027] The total number of samples is the total number of sample points currently participating in the synthetic autocorrelation calculation and smoothing optimization calculation.

[0028] The sample sequence number is a discrete number used to identify the position of a single sample in the current sequence.

[0029] A smooth sequence is the target displacement sequence to be solved in a smooth optimization model.

[0030] Smoothing strength is a weighting parameter that balances the fitting error of the original displacement sequence with the strength of the second-order difference smoothing constraint.

[0031] The optimal smoothing sequence for the current channel is the optimal displacement sequence obtained by the current channel after satisfying the smoothing optimization objective.

[0032] The optimal smoothing strength is the smoothing strength that is adaptively selected by the generalized cross-validation method to achieve the best trade-off between prediction residuals and model complexity.

[0033] A smoothing matrix is ​​a linear operator matrix that maps an original displacement sequence to a smoothed sequence for a given smoothing strength.

[0034] The smoothed displacement sequence is a two-dimensional displacement sequence composed of two channels that have been smoothed separately.

[0035] The optimal smoothed sequence xn of the dual channels is the optimal displacement sequence obtained after the smoothing optimization of the first channel.

[0036] The optimal smoothed sequence yn for the two channels is the optimal displacement sequence obtained after the second channel has been smoothed.

[0037] In detail, because the trajectory of the water pump shaft may experience speed fluctuations, thermal changes, and pseudo-period drift under transitional operating conditions, using a fixed-length window or relying on external key phase signal windowing will fragment a locally complete trajectory loop. Therefore, using the sample length corresponding to the first positive local maximum of the synthesized autocorrelation function as the dominant period can make the trajectory window more stably correspond to a locally complete loop.

[0038] In detail, because the slow migration of the axis center itself is a real evolutionary information that needs to be preserved, if a high-pass detrending or simple baseline subtraction is used directly, the average position change that must be used for subsequent prediction will be mistakenly deleted. Therefore, a smooth optimization model with second-order difference constraints is adopted to retain the slow-changing center information while suppressing spikes and high-frequency noise.

[0039] In detail, because the required smoothing intensity varies under different working conditions, different measuring points, and different noise levels, fixing empirical smoothing parameters can lead to undersmoothing or oversmoothing. Therefore, using the generalized cross-validation method to adaptively calculate the optimal smoothing intensity can automatically achieve a balance between data fitting and smoothing constraints in the smoothing results.

[0040] In detail, the method for determining the first positive local maximum of the synthetic autocorrelation function is as follows: First, calculate the synthetic autocorrelation function value point by point according to the sample lag from small to large. Then, calculate the function difference under the adjacent sample lag. When the current difference is greater than 0, the next difference is less than or equal to 0, and the current synthetic autocorrelation function value is greater than 0, the sample lag is determined to be a positive local maximum, and the smallest sample lag that meets the conditions is taken as the dominant period.

[0041] In detail, the backoff method for the synthetic autocorrelation function when there is no positive local maximum is as follows: within a search interval not exceeding half of the total number of samples, the sample lag with the largest synthetic autocorrelation function value is selected as the dominant period to ensure that each original displacement sequence can obtain an effective time slice length.

[0042] In detail, the integerization method for the number of future windows is as follows: first, divide the preset prediction time by the product of the dominant period and the sampling interval to obtain the theoretical window number, then round up the theoretical window number, and limit the result to a positive integer not less than 1.

[0043] In detail, the smoothing matrix is ​​formed as follows: First, a second-order difference discretization operator is constructed based on the total number of samples. Then, the identity matrix and the smoothing intensity are multiplied by the transpose of the second-order difference discretization operator and itself, and the result is added together. Finally, the inverse of the resulting matrix is ​​obtained, thus forming a linear mapping matrix from the original displacement sequence to the smoothed sequence.

[0044] In detail, the solution method of the smoothing optimization model is as follows: establish smoothing optimization objectives independently for the first channel and the second channel respectively, solve the optimal smoothing sequence of the current channel using a system of linear equations after giving the smoothing intensity, and combine the solution results of the two channels into a smoothed displacement sequence.

[0045] In detail, the numerical method for obtaining the optimal smoothing strength is as follows: perform a one-dimensional search for the smoothing strength within the range of positive real numbers, calculate the generalized cross-validation target value point by point, and take the smoothing strength that minimizes the target value as the optimal smoothing strength.

[0046] Preferably, the parameter extraction module extracts the axis center and hysteresis area of ​​the smoothed trajectory window, including: Calculate the trajectory window The smoothed displacement sequence described above The arithmetic mean of the axes is used to obtain the center of the axis. :

[0047] The smoothed displacement sequence With the axis center By subtracting, we obtain the centered trajectory sequence. :

[0048] Calculate the centered trajectory sequence The covariance is used to obtain the local loop shape matrix. :

[0049] set up and The local loop shape matrix eigenvalues, and Extract the half length of the spindle and the second shaft half length And calculate the principal axis direction angle :

[0050]

[0051] The centered trajectory sequence is closed-loop connected at the beginning and end of time, let and Calculate the geometric area enclosed by the closed and connected sequence of centered trajectories, and use it as the area of ​​the lap. :

[0052] in, and All of these represent the number of samples in the trajectory window. The sample number. , and These are the elements of the local loop shape matrix. and These are the components of the centered trajectory sequence in the dual-channel direction. and These are the components of the first sample within the trajectory window in the dual-channel direction. and These are the components of the newly added last sample in the dual-channel direction after the closed-connection process.

[0053] The axis center is the arithmetic mean position vector of the smoothed displacement sequence within the trajectory window in the dual-channel direction.

[0054] A centered trajectory sequence is a local trajectory sequence obtained by subtracting the axis center from each smoothed displacement sample within the trajectory window.

[0055] The local loop shape matrix is ​​a two-dimensional symmetric matrix composed of the covariance of the centered trajectory sequence. It is used to characterize the degree of dispersion and coupling relationship of the local loop in the two channel directions.

[0056] Elements of the local loop shape matrix It is the variance term of the centered trajectory sequence in the direction of the first channel.

[0057] Elements of the local loop shape matrix It is the covariance term of the centered trajectory sequence in the two channel directions.

[0058] Elements of the local loop shape matrix It is the variance term of the centered trajectory sequence in the direction of the second channel.

[0059] Eigenvalues It is a large eigenvalue of the local loop shape matrix, used to characterize the expansion intensity in the main direction.

[0060] Eigenvalues It is a small eigenvalue of the local loop shape matrix, used to characterize the expansion intensity in the secondary direction.

[0061] The principal axis half-length is the length of the principal direction half-axis of the local loop obtained by taking the square root of the larger eigenvalue of the local loop shape matrix.

[0062] The secondary axis semi-length is the length of the secondary axis of the local loop obtained by taking the square root of the smaller eigenvalue of the local loop shape matrix.

[0063] The main axis direction angle is the direction angle of the local loop main axis relative to the reference direction of the first channel.

[0064] The lap area is the geometric area enclosed by the closed loop of the centered trajectory sequence.

[0065] The number of samples in the trajectory window is the total number of sample points within the trajectory window.

[0066] Number of samples in the trajectory window It is a representation of the number of window samples used for calculating the area of ​​the closed loop.

[0067] Components of the centered trajectory sequence in the dual-channel direction It is the component value of the nth centered trajectory sample in the direction of the first channel.

[0068] Components of the centered trajectory sequence in the dual-channel direction It is the component value of the nth centered trajectory sample in the direction of the second channel.

[0069] The components of the first sample in the trajectory window in the dual-channel direction It is the component value of the first sample in the first channel direction of the centered trajectory sequence.

[0070] The components of the first sample in the trajectory window in the dual-channel direction It is the component value of the first sample in the second channel direction of the centered trajectory sequence.

[0071] The components of the newly added end sample in the dual-channel direction after closed-loop processing. It is the first channel component formed after the first sample is copied to the end of the sequence.

[0072] The components of the newly added end sample in the dual-channel direction after closed-loop processing. It is the second channel component formed after the first sample is copied to the end of the sequence.

[0073] In detail, because the smoothed displacement sequence within the trajectory window contains both average position changes and local loop changes, directly using the original two-dimensional displacement points to describe the loop shape would mix rigid body translation with shape changes. Therefore, the arithmetic mean is used to extract the axis center first, and then a centered trajectory sequence is constructed around the axis center.

[0074] In detail, since the expansion scale and directionality of the local loop in both channel directions can be stably characterized by the covariance structure, the shape matrix of the local loop can be calculated and eigenvalue decomposed on the centered trajectory sequence to obtain the principal axis half-length, secondary axis half-length and principal axis direction angle.

[0075] In detail, based on the expansion and contraction of the local loop itself, rather than the apparent displacement caused by the average position translation, closing the beginning and end of the centered trajectory sequence and calculating the loop area according to the polygon area principle can effectively separate positional changes from shape changes.

[0076] In detail, the extraction method of the principal axis half length and the secondary axis half length is as follows: first, perform eigenvalue decomposition on the local loop shape matrix, then sort the eigenvalues ​​from largest to smallest, take the square root of the larger eigenvalue to obtain the principal axis half length, and take the square root of the smaller eigenvalue to obtain the secondary axis half length.

[0077] In detail, the method for making the main axis direction angle continuous is as follows: First, calculate the main axis direction angle of the current trajectory window, and then subtract it from the main axis direction angle of the previous trajectory window. When the difference is greater than 90 degrees, subtract 180 degrees from the current main axis direction angle. When the difference is less than -90 degrees, add 180 degrees to the current main axis direction angle, thereby maintaining the continuity of the main axis direction angle of adjacent trajectory windows.

[0078] In detail, the closed connection method of the centralized trajectory sequence at the beginning and end of time is as follows: the successor point of the last centralized trajectory point in the trajectory window is designated as the first centralized trajectory point in the trajectory window, and the two channel components of the first sample are copied as the two channel components of the newly added last sample.

[0079] In detail, the area sign of the loop area is handled as follows: first, calculate the sum of the directed area terms segment by segment according to the centered trajectory sequence after the beginning and end are closed, and then take half of its absolute value as the loop area to ensure that the loop area is a non-negative geometric quantity.

[0080] Preferably, the parameter extraction module calculates the center migration intensity based on the axis center of adjacent trajectory windows and calculates the area growth rate based on the loop area of ​​adjacent trajectory windows, including: The center migration intensity is obtained by calculating the L2 norm of the vector difference between the axis center of the current trajectory window and the axis center of the previous trajectory window. :

[0081] The area growth rate is obtained by calculating the difference between the natural logarithm of the lap area of ​​the current trajectory window and the natural logarithm of the lap area of ​​the previous trajectory window. :

[0082] in, The center migration intensity, The axis center of the current trajectory window, The axis center of the previous trajectory window, Represents the L2 norm operation; This refers to the area growth rate. The area of ​​the loop in the current trajectory window. The area of ​​the loop in the previous trajectory window. Represents the natural logarithm operation.

[0083] The center migration strength is the L2 norm of the difference vector between the axis center of the current trajectory window and the axis center of the previous trajectory window.

[0084] The area growth rate is the difference between the natural logarithm of the current trajectory window loop area and the natural logarithm of the previous trajectory window loop area.

[0085] In detail, since the actual migration of the pump shaft center may occur along the first channel direction or the second channel direction, comparing only the single channel increment will lose two-dimensional displacement information. Therefore, the L2 norm of the difference vector between adjacent shaft centers is used as the center migration intensity.

[0086] In detail, since the area of ​​the loop may expand or shrink exponentially under different operating conditions, directly using the area difference is easily affected by the absolute scale. Therefore, the difference of the natural logarithm of the areas of adjacent loops is used as the area growth rate, which can transform the relative expansion process into a comparable increment.

[0087] In detail, the window alignment method for center migration intensity is as follows: First, number all axis centers according to the time sequence of the trajectory window, then pair the axis center of the current trajectory window with the axis center of the previous trajectory window one by one, subtract them, and calculate the L2 norm to obtain the center migration intensity.

[0088] In detail, the logarithmic safety handling method for the area growth rate is as follows: the loop area should be the loop area after absolute value processing. When the loop area of ​​the current trajectory window or the previous trajectory window is 0, the area growth rate is not updated and the previous valid area growth rate is used. Normal calculation is resumed after the loop area becomes positive again.

[0089] Preferably, the parameter extraction module calculates the maximum time-lag correlation between the center migration intensity and the area growth rate to construct a lag index, including: Calculate the normalized time-delay correlation coefficient between the center migration intensity and the area growth rate under different lags. :

[0090] Find the optimal lag that maximizes the normalized time-delay correlation coefficient. and the corresponding maximum time-delay correlation. :

[0091]

[0092] The optimal hysteresis Correlation with the maximum time delay Multiply to obtain the lag index. :

[0093] in, It is a lag quantity. This refers to the sequence number of the current trajectory window. The center migration intensity, This is the average value of the center migration intensity corresponding to multiple consecutive trajectory windows. The area growth rate is given when the lag exists. This is the average of the area growth rate corresponding to multiple consecutive trajectory windows. The normalized time-delay correlation coefficient is... This is the optimal hysteresis. The maximum time-delay correlation, The lag index is... To perform the maximum value operation, The operation is performed to find the value of the independent variable that maximizes the objective function.

[0094] The first difference sequence is obtained by subtracting the average center migration intensity from the center migration intensity corresponding to multiple consecutive trajectory windows.

[0095] The second difference sequence is obtained by subtracting the average area growth rate from the area growth rate corresponding to multiple consecutive trajectory windows with lag.

[0096] The average value of the center migration intensity corresponding to multiple consecutive trajectory windows is the arithmetic mean of the center migration intensity within the consecutive window segments currently involved in the time delay correlation calculation.

[0097] The average of the area growth rate corresponding to multiple consecutive trajectory windows is the arithmetic mean of the area growth rate within the current consecutive window segment participating in the time delay correlation calculation.

[0098] The normalized time-lag correlation coefficient is the normalized correlation between the first difference sequence and the second difference sequence under a given lag.

[0099] The hysteresis is the number of trajectory windows in which the area growth rate shifts backward relative to the intensity of center migration.

[0100] The optimal lag is the lag that maximizes the normalized time-delay correlation coefficient.

[0101] The maximum time-delay correlation coefficient is the largest value among the normalized time-delay correlation coefficients corresponding to all candidate time delays.

[0102] The lag index is the product of the optimal lag and the maximum time-delay correlation.

[0103] The current trajectory window number is a discrete number used to indicate the temporal order of the trajectory windows.

[0104] In detail, because the throttling and sealing effect inside the water pump and the thermal effect of the sliding bearing are not synchronized, the intensity of center migration and the area growth rate often do not reach their peak values ​​simultaneously in the same window. Therefore, by searching for the maximum time-delay correlation coefficient using the normalized time-delay correlation coefficient, the time-staggered relationship between the center's initial movement and the hysteresis loop's subsequent expansion can be directly quantified.

[0105] In detail, since the optimal lag alone cannot reflect the strength of the correlation, and the maximum time-delay correlation alone cannot reflect the time interval between events, the lag exponent is formed by multiplying the optimal lag and the maximum time-delay correlation. This exponent can simultaneously compress and express both the time interval between events and the coupling strength.

[0106] In detail, the selection method for multiple consecutive trajectory windows is as follows: taking the current trajectory window as the endpoint, select a window segment that is continuous in time and whose center migration intensity and area growth rate are defined as the calculation interval, and ensure that the window segment contains no less than 5 valid trajectory windows.

[0107] In detail, the lag search method is as follows: on the premise that there are overlapping samples between the center migration intensity sequence and the area growth rate sequence, the lag is increased one by one starting from 0 until the two sequences no longer meet the overlap calculation conditions after further lag is added.

[0108] In detail, the parallel processing method for optimal lags is as follows: when there are multiple lags with the same maximum normalized time lag correlation coefficient, the positive lag with the smallest value is selected as the optimal lag to ensure that the earliest occurrence relationship is retained first.

[0109] In detail, the maximum time-delay correlation coefficient is determined as follows: only the largest positive value among the normalized time-delay correlation coefficients is taken as the maximum time-delay correlation coefficient. When all normalized time-delay correlation coefficients are not positive, the maximum time-delay correlation coefficient is set to 0, and the optimal lag is set to 0.

[0110] Preferably, the prediction extrapolation module uses the lag exponent to extrapolate the future center trajectory, including: Calculate the difference between the axis center of the current trajectory window and the axis center of the previous trajectory window to obtain the first-order change vector of the center. :

[0111] Subtracting twice the center of the previous trajectory window from the center of the current trajectory window, and then adding the center of the second trajectory window from the previous one, yields the second-order center change vector. :

[0112] For the set number of extrapolation steps The current trajectory window's axis center, the first predicted component obtained by multiplying the extrapolation step number by the first-order change vector of the center, and the extrapolation step number and the maximum time-delay correlation are used to determine the predicted component. The second predicted component, formed by multiplying the second-order change vector of the center, is added to the first predicted component, and the result is taken as the future center trajectory. :

[0113] in, Let be the first-order change vector of the center. The axis center of the current trajectory window, The axis center of the previous trajectory window, The central second-order transformation vector, The axis center is the second trajectory window that is pushed forward. The extrapolation step number is... The maximum time-delay correlation, This refers to the future center trajectory.

[0114] The first-order change vector of the center is the difference between the axis center of the current trajectory window and the axis center of the previous trajectory window.

[0115] The central second-order change vector is a discrete second-order difference vector formed by the axis centers of the current trajectory window, the previous trajectory window, and the second-to-last trajectory window.

[0116] The extrapolation step number is an integer sequence number that progressively numbers the future windows. It is preferably from 1 to the number of future windows, because the extrapolation step number is used to enumerate every future trajectory window within the preset prediction time. Exceeding the number of future windows will cause the trajectory to fall outside the prediction range.

[0117] The first predictive component is the linear trend component obtained by multiplying the extrapolation steps by the central first-order change vector.

[0118] The second prediction component is the curvature correction component, which is formed by the quadratic combination of extrapolation steps, the maximum time-delay correlation, and the central second-order change vector.

[0119] The future center trajectory is the future axis center position vector predicted under the corresponding extrapolation step number.

[0120] In detail, since the future center trajectory is affected by both the continuous effect of the current migration trend and the curvature change reflected by the time-staggered coupling of center migration and loop outward expansion, incorporating the first-order center change vector and the second-order center change vector weighted by the maximum time-delay correlation into the extrapolation of the future center trajectory can simultaneously preserve trend information and curvature information.

[0121] In detail, because the prediction object is the axis center rather than the entire original waveform, the first and second prediction components are directly superimposed on the axis center of the current trajectory window.

[0122] In detail, the initial conditions for determining the future center trajectory are as follows: the future center trajectory extrapolation is only started when the axis centers of the current trajectory window, the previous trajectory window, and the second trajectory window before it have all been calculated; otherwise, only parameter accumulation is performed and the future center trajectory is not output.

[0123] In detail, the channel-by-channel calculation method for the first and second predicted components is as follows: First, calculate the first-order change vector and the second-order change vector of the center for the axis center components on the two channels respectively. Then, calculate the first and second predicted components of the two channels respectively with the same extrapolation steps and the same maximum time delay correlation. Finally, synthesize the future center trajectory.

[0124] In detail, the extrapolation step number is determined as follows: the extrapolation step number starts from 1 and increases by integers until the number of future windows ends. Each extrapolation step number corresponds to the future center trajectory of a future trajectory window.

[0125] Preferably, the prediction extrapolation module, after delaying the center migration intensity based on the maximum time-delay correlation, maps it to the shape parameters of the future hysteresis loop, including: Calculate the half length of the spindle respectively With the half length of the secondary shaft The natural logarithm value is used to obtain the half-length of the logarithmic principal axis. With the half length of the logarithmic secondary axis :

[0126] When performing natural logarithmic calculations, the semi-lengths of the main shaft and the secondary shaft must first undergo dimensionless preprocessing with a fixed reference. The fixed reference used in the preprocessing is the radial design clearance value of the bearing calibrated during the equipment installation and commissioning phase, which remains constant throughout the entire operating cycle of the equipment. The semi-lengths of the main shaft and the secondary shaft to be calculated are divided by this fixed reference to obtain the dimensionless relative ratios of the semi-lengths of the main shaft and the secondary shaft, respectively. Then, the natural logarithmic calculation in this step is performed.

[0127] Using the optimal hysteresis Determine the effective index set of historical sequences And calculate the average value of the center migration intensity after the delay. :

[0128]

[0129] The center migration strength involved in the calculation needs to undergo dimensionless preprocessing based on a fixed reference. The fixed reference used in the preprocessing is completely consistent with the fixed reference for bearing radial design clearance used in the natural logarithmic calculation of the main shaft half-length and the secondary shaft half-length. This reference remains constant throughout the entire operating cycle of the equipment. The center migration strength to be calculated is divided by this fixed reference to obtain the dimensionless relative ratio of the center migration strength. Then, it is used to solve the mapping gain and calculate the recursive formula to ensure that all physical quantities involved in the addition and subtraction operations are dimensionless values.

[0130] Calculate the first mapping gain respectively Second mapping gain With the third mapping gain :

[0131]

[0132]

[0133] Based on the maximum time delay correlation By recursively working backwards, we can obtain the half-length of the principal axis of the future logarithm for each step. Future logarithmic subaxis half length With future principal axis direction angle :

[0134]

[0135]

[0136] With the extrapolation steps Continuous recursion, for the th The future logarithmic principal axis half length of the step With the future logarithmic secondary axis half length Perform natural index reduction calculations to obtain the future principal axis half-length. With the future secondary shaft half length :

[0137] The extrapolation steps are as described above. The following is the future spindle half length The future sub-shaft half-length and future main axis direction angle As the shape parameter of the future lap; in, This refers to the sequence number of the current trajectory window. The traversal sequence number within the valid index set. The delayed center migration strength, The main axis direction angle is... The number of elements in the valid index set. For natural logarithm operations, It is a natural constant.

[0138] The logarithmic principal axis half-length is the length state quantity obtained by performing a natural logarithmic transformation on the principal axis half-length.

[0139] The logarithmic subaxis half-length is the length state quantity obtained by performing a natural logarithmic transformation on the subaxis half-length.

[0140] The effective index set of historical sequences is the set of historical indexes that ensures both the difference between adjacent sequences and the center migration strength after delay when calculating the mapping gain.

[0141] The average value of delayed center migration strength is the arithmetic mean of delayed center migration strength within the effective index set of the historical sequence.

[0142] The first mapping gain is the mapping coefficient of the delayed center migration strength to the difference between adjacent sequences at half the logarithmic principal axis.

[0143] The second mapping gain is the mapping coefficient of the delayed center migration strength to the difference between adjacent sequences at the half length of the logarithmic subaxis.

[0144] The third mapping gain is the mapping coefficient of the delayed center migration intensity to the difference between adjacent sequences of the principal axis direction angle.

[0145] The single-step future logarithmic principal axis half-length is a single-step recursive prediction of the logarithmic principal axis half-length of the next trajectory window.

[0146] The single-step future logarithmic subaxis half-length is a single-step recursive prediction of the logarithmic subaxis half-length of the next trajectory window.

[0147] The future principal axis direction angle of a single step is a single-step recursive prediction of the principal axis direction angle of the next trajectory window.

[0148] The half-length of the future logarithmic principal axis at step h is the recursive result of the half-length of the logarithmic principal axis corresponding to the h-th future trajectory window.

[0149] The future logarithmic subaxis semi-length at step h is the recursive result of the logarithmic subaxis semi-length corresponding to the h-th future trajectory window.

[0150] The future principal axis half-length is obtained by restoring the future logarithmic principal axis half-length at step h to its natural exponential form.

[0151] The future subaxis semi-length is obtained by restoring the future logarithmic subaxis semi-length at step h to its natural exponential form.

[0152] The future principal axis direction angle is the principal axis direction angle obtained recursively from the h-th future trajectory window.

[0153] The traversal sequence number within the valid index set is a discrete number used when traversing each index in the valid index set of the historical sequence.

[0154] The number of elements in the effective index set is the total number of valid historical samples within the effective index set of the historical sequence.

[0155] In detail, since the semi-lengths of the principal axis and the secondary axis must remain positive, and their changes are closer to relative proportional changes than absolute difference changes, taking the natural logarithm of the semi-lengths of the principal axis and the secondary axis first, and then recursively extrapolating within the logarithmic domain, can improve the stability of future lap shape parameter predictions.

[0156] In detail, since the delayed center migration intensity is regarded as the leading driver of future changes in the shape parameters of the loop, the correspondence between the delayed center migration intensity and the logarithmic principal axis half-length, logarithmic secondary axis half-length, and principal axis direction angle differences between adjacent sequences can be obtained by statistically analyzing the effective index set of historical sequences.

[0157] In detail, since the future principal axis half-length and the future secondary axis half-length will eventually directly participate in the generation of future boundary points, performing natural exponential restoration after completing the logarithmic domain recursion can restore the intermediate prediction results to the future hysteresis shape parameters.

[0158] In detail, the effective index set of the historical sequence is constructed as follows: only historical indexes that simultaneously satisfy the conditions that the difference between adjacent sequences is defined, the center migration strength after delay is defined, and the index order does not exceed the boundary are retained, and the effective index set of the historical sequence is formed in chronological order.

[0159] In detail, the processing method when the denominators of the first mapping gain, the second mapping gain, and the third mapping gain are 0 is as follows: when the delayed center migration intensity does not have discrete fluctuations in the effective index set of the historical sequence, the first mapping gain, the second mapping gain, and the third mapping gain are uniformly set to 0, so that the prediction degenerates into an inertial recursion that only depends on the difference between adjacent sequences.

[0160] In detail, the multi-step recursive method for the future logarithmic principal axis half-length, the future logarithmic secondary axis half-length, and the future principal axis direction angle is as follows: when the required delayed center migration intensity is still within the historical interval, the historical value is directly called; when the required delayed center migration intensity has crossed the historical interval, the corresponding future center migration intensity is first calculated based on the future center trajectory, and then the future center migration intensity is substituted into the subsequent recursion.

[0161] In detail, the continuous processing method of the future principal axis direction angle is as follows: after each step of recursion is completed, the current future principal axis direction angle is compared with the previous future principal axis direction angle. When the difference is greater than 90 degrees, 180 degrees are subtracted from the current future principal axis direction angle. When the difference is less than -90 degrees, 180 degrees are added to the current future principal axis direction angle to maintain the continuity of the angle sequence.

[0162] Preferably, the result output module generates a future boundary envelope using the future center trajectory and the shape parameters, and calculates the radial distance difference between the farthest point and the nearest point of the future boundary envelope with the bearing clearance center as the origin, and outputs it as the maximum radial runout amplitude, including: Using the future principal axis direction angle Construct rotation matrix :

[0163] Combined with parameter angle The future spindle half-length With the future secondary shaft half length Generate boundary coordinates using the rotation matrix. Rotate the boundary coordinates and compare the rotation result with the future center trajectory. Add them together to get the number of extrapolation steps. The following future boundary point :

[0164] All set future window numbers that need to be extrapolated The set of all the future boundary points within the specified range is merged to obtain the set of future boundary points. :

[0165] Calculate the set of future boundary points The convex hull is used as the future boundary envelope. :

[0166] Using the calibrated bearing clearance center as the origin, calculate the future boundary envelope respectively. All points on radial distance to the origin of the coordinate system :

[0167] The maximum and minimum values ​​of the radial distances are selected, and the point corresponding to the maximum value is taken as the farthest point. The point corresponding to the minimum value is taken as the nearest point. :

[0168]

[0169] Subtracting the minimum value from the maximum value yields the maximum radial runout amplitude. And output it:

[0170] in, This is the sequence number of the current trajectory window. The future principal axis direction angle, For the parameter angle, For the future spindle half-length, For the future sub-axis half-length, For the future center trajectory, For the future boundary point, The extrapolation step number is... The number of future windows, Let be the set of future boundary points. For the union operation, For convex hull operations, Let the future boundary envelope be... For points on the future boundary envelope, The radial distance is... Represents the L2 norm operation. For the farthest point, For the nearest point, The operation is to find the value of the independent variable that maximizes the objective function. The operation is to find the value of the independent variable that minimizes the objective function. The maximum radial runout amplitude is given.

[0171] A rotation matrix is ​​a two-dimensional orthogonal rotation matrix constructed from the future principal axis direction angle, used to rotate standard boundary coordinates to the future principal axis direction.

[0172] The parameter angle is an angle variable used to iterate through the future lap boundary in various directions.

[0173] The future boundary point is the boundary coordinate point determined by the future center trajectory, the future principal axis half length, the future secondary axis half length, and the rotation matrix under a certain extrapolation step and a certain parameter angle.

[0174] The set of future boundary points is the set of points formed by merging all future boundary points within a preset prediction time.

[0175] The future boundary envelope is the convex hull boundary obtained from the set of future boundary points.

[0176] The bearing clearance center is the geometric center reference point of the bearing, which is used as the origin of the coordinate system after installation and calibration.

[0177] A point on the future boundary envelope is any point located on the future boundary envelope and is used to calculate the radial distance.

[0178] The radial distance is the L2 norm distance from a point on the future boundary envelope to the center of the bearing clearance.

[0179] The farthest point is the point with the largest radial distance on the future boundary envelope.

[0180] The closest point is the point with the smallest radial distance on the future boundary envelope.

[0181] The maximum radial runout amplitude is the difference between the radial distance of the farthest point and the radial distance of the nearest point.

[0182] In detail, since the future loop boundary is essentially a local closed boundary that unfolds around the future center trajectory and is oriented along the future principal axis, the future boundary points can be generated by constructing a rotation matrix using the future principal axis direction angle and combining the future principal axis half-length, the future secondary axis half-length, and the parameter angle to transform the future loop shape parameters into geometric boundaries.

[0183] In detail, since the upper bound of risk within the preset prediction time is not determined by a single future trajectory window, but by all future trajectory windows, merging the future boundary points under all extrapolation steps into a set of future boundary points and obtaining the convex hull can yield a future boundary envelope covering the entire prediction time window.

[0184] In detail, since the engineering significance of the maximum radial runout amplitude is the maximum radial expansion range of the future boundary relative to the center of the bearing clearance, the maximum radial runout amplitude is output by taking the center of the bearing clearance as the origin and using the radial distance difference between the farthest point and the nearest point on the future boundary envelope. This can correspond to the upper limit of equipment operation risk.

[0185] In detail, the parameter angle is discretized by discretizing the parameter angle in the range of 0 degrees to 360 degrees with equal angular step size, and ensuring that the number of discrete parameter angles in a future trajectory window is not less than the number of samples corresponding to the dominant period, so as to form a sufficiently dense future boundary point.

[0186] In detail, the merging of the future boundary point set and the calculation of the convex hull of the future boundary envelope are as follows: First, write all future boundary points into a unified point set in the order of extrapolation steps, then execute the convex hull algorithm according to the planar coordinates of the points, and finally output the convex hull vertices as the future boundary envelope in a counterclockwise order.

[0187] In detail, the calibration method for the bearing clearance center is as follows: During the equipment installation and commissioning phase, the bearing clearance center is determined by combining the bearing geometric center and the zero-position calibration results of the dual probes, and then the bearing clearance center is fixed as the coordinate origin in subsequent online calculations.

[0188] In detail, the search method for the farthest and nearest points is as follows: First, calculate the radial distance from all vertices of the future boundary envelope to the center of the bearing clearance. Then, for each boundary line segment, determine whether the foot of the perpendicular from the center of the bearing clearance to the line segment is inside the line segment. If the foot of the perpendicular is inside the line segment, compare the radial distance of the foot of the perpendicular. Finally, determine the point corresponding to the maximum value among all candidate distances as the farthest point and the point corresponding to the minimum value as the nearest point.

[0189] In detail, the output method of the maximum radial runout amplitude is as follows: at least the maximum radial runout amplitude, the coordinates of the farthest point, the coordinates of the nearest point, and the sequence of future boundary envelope vertices are output synchronously so that the host computer can display and trigger alarms.

[0190] It should be noted that the interval and threshold sizes are set for ease of comparison. The size of the threshold depends on the amount of sample data and the base number set by those skilled in the art for each set of sample data, as long as it does not affect the proportional relationship between the parameter and the quantized value. Furthermore, the above formulas are all dimensionless calculations, and the formulas are derived from software simulations using a large amount of collected data to obtain the most recent real-world results. The preset parameters in the formulas are set by those skilled in the art according to the actual situation.

[0191] The embodiments of this example have been described above. However, this example is not limited to the specific implementation methods described above. The specific implementation methods described above are merely illustrative and not restrictive. Those skilled in the art can make many other forms based on the guidance of this example, and all of them are within the protection scope of this example.

Claims

1. A dynamic prediction system for the shaft trajectory of a water pump based on eddy current displacement signals, characterized in that, include: The signal preprocessing module is used to acquire the original displacement sequence through dual probes, extract the dominant period using the synthetic autocorrelation function to divide multiple trajectory windows, and smooth the original displacement sequence of each trajectory window. The parameter extraction module is used to extract the axis center and loop area of ​​the smoothed displacement sequence within the trajectory window, calculate the center migration intensity based on the axis center of the adjacent trajectory windows, calculate the area growth rate based on the loop area of ​​the adjacent trajectory windows, and calculate the maximum time lag correlation between the center migration intensity and the area growth rate to construct a lag index. The prediction extrapolation module is used to extrapolate the future center trajectory using the lag exponent, and after delaying the center migration intensity according to the maximum time lag correlation, it is mapped to the shape parameters of the future hysteresis loop. The result output module is used to generate a future boundary envelope using the future center trajectory and the shape parameters. Taking the bearing clearance center as the origin, it calculates the radial distance difference between the farthest point and the nearest point of the future boundary envelope and outputs it as the maximum radial runout amplitude.

2. The dynamic prediction system for the shaft center trajectory of a water pump based on eddy current displacement signals according to claim 1, characterized in that, The signal preprocessing module acquires the original displacement sequence through dual probes, extracts the dominant period using a synthetic autocorrelation function to divide the trajectory into multiple trajectory windows, and performs smoothing processing on the original displacement sequence of each trajectory window, including: The original displacement sequence acquired by the dual probes is obtained, the synthetic autocorrelation function of the original displacement sequence is calculated, and the sample length corresponding to the first positive local maximum of the synthetic autocorrelation function is taken as the dominant period. The original displacement sequence is sliced ​​according to the dominant period to obtain the multiple trajectory windows, and the number of future windows to be extrapolated is calculated according to the preset prediction time and the dominant period. For the original displacement sequence within the trajectory window, a smoothing optimization model with second-order difference constraints is constructed. The optimal smoothing intensity is adaptively calculated using the generalized cross-validation method. The original displacement sequence is then solved based on the optimal smoothing intensity to obtain the smoothed displacement sequence.

3. The dynamic prediction system for the shaft trajectory of a water pump based on eddy current displacement signals according to claim 2, characterized in that, The parameter extraction module extracts the axis center and hysteresis area of ​​the smoothed displacement sequence within the trajectory window, including: Calculate the arithmetic mean of the smoothed displacement sequence within the trajectory window to obtain the axis center; The smoothed displacement sequence is subtracted from the axis center to obtain the centered trajectory sequence; The covariance of the centered trajectory sequence is calculated to obtain the local loop shape matrix. The local loop shape matrix is ​​then decomposed into eigenvalues ​​to extract the principal axis half-length and the secondary axis half-length, and the principal axis direction angle is calculated. The centered trajectory sequence is closed and connected at the beginning and end of time, and the geometric area enclosed by the closed and connected centered trajectory sequence is calculated based on the principle of polygon area calculation, and it is used as the loop area.

4. The dynamic prediction system for the shaft trajectory of a water pump based on eddy current displacement signals according to claim 3, characterized in that, The parameter extraction module calculates the center migration intensity based on the axis center of adjacent trajectory windows and calculates the area growth rate based on the loop area of ​​adjacent trajectory windows, including: Calculate the L2 norm of the vector difference between the axis center of the current trajectory window and the axis center of the previous trajectory window, and use the result as the center migration intensity; Extract the lap area of ​​the current trajectory window and the lap area of ​​the previous trajectory window, calculate the natural logarithm of each, subtract the natural logarithm of the previous trajectory window from the natural logarithm of the current trajectory window, and use the difference as the area growth rate.

5. The dynamic prediction system for the shaft trajectory of a water pump based on eddy current displacement signals according to claim 4, characterized in that, The parameter extraction module calculates the maximum time-lag correlation between the center migration intensity and the area growth rate to construct a lag index, including: The first difference sequence is obtained by subtracting the center migration intensity corresponding to multiple consecutive trajectory windows from their own arithmetic mean. The second difference sequence is obtained by subtracting the area growth rate corresponding to the multiple consecutive trajectory windows with lag from their own arithmetic mean. Calculate the sum of the products of the first difference sequence and the second difference sequence, and divide it by the product of the square root of the sum of squares of the first difference sequence and the square root of the sum of squares of the second difference sequence to obtain the normalized time-delay correlation coefficient. Find the lag that maximizes the normalized time-delay correlation coefficient, take it as the optimal lag, and take the corresponding normalized time-delay correlation coefficient as the maximum time-delay correlation degree; The optimal lag is multiplied by the maximum time-delay correlation, and the product is used as the lag index.

6. The dynamic prediction system for the shaft center trajectory of a water pump based on eddy current displacement signals according to claim 5, characterized in that, The prediction extrapolation module uses the lag exponent to extrapolate the future center trajectory, including: Calculate the difference between the axis center of the current trajectory window and the axis center of the previous trajectory window to obtain the first-order change vector of the center; Subtract twice the center of the previous trajectory window from the center of the current trajectory window, and add the center of the second trajectory window to the left to obtain the second-order change vector of the center. For a given number of extrapolation steps, the extrapolation steps are multiplied by the first-order change vector of the center to obtain the first prediction component; Calculate the product of the extrapolation step number and the difference between the extrapolation step number and one, divide it by two, and then multiply it by the maximum time-delay correlation degree and the central second-order change vector to obtain the second prediction component; The axis center, the first predicted component, and the second predicted component of the current trajectory window are added together, and the result of the addition is taken as the future center trajectory.

7. The dynamic prediction system for the shaft trajectory of a water pump based on eddy current displacement signals according to claim 6, characterized in that, The prediction extrapolation module, after delaying the center migration intensity based on the maximum time-delay correlation, maps it to the shape parameters of the future hysteresis loop, including: Calculate the natural logarithm values ​​of the principal axis half-length and the secondary axis half-length respectively to obtain the logarithmic principal axis half-length and the logarithmic secondary axis half-length; The effective index set of the historical sequence is determined using the optimal lag, and the average value of the center migration intensity after the delay is calculated. The adjacent sequence differences of the logarithmic principal axis half length, the logarithmic secondary axis half length and the principal axis direction angle are calculated respectively. Each adjacent sequence difference is multiplied by the deviation between the delayed center migration intensity and its average value and then summed. The sum of these values ​​is then divided by the sum of the squares of the deviations to obtain the first mapping gain, the second mapping gain and the third mapping gain in sequence. The future logarithmic principal axis half-length is obtained by adding the current logarithmic principal axis half-length, the difference between its adjacent sequences, the maximum time-delay correlation, the first mapping gain, and the delayed center migration intensity. Similarly, the future logarithmic secondary axis half-length and the future principal axis direction angle are obtained by recursively using the second mapping gain and the third mapping gain. The future logarithmic principal axis half-length and the future logarithmic secondary axis half-length are calculated using natural exponential reduction for each extrapolation step to obtain the future principal axis half-length and the future secondary axis half-length; the future principal axis half-length, the future secondary axis half-length, and the future principal axis direction angle set are used as the shape parameters of the future hysteresis loop.

8. The dynamic prediction system for the shaft trajectory of a water pump based on eddy current displacement signal according to claim 7, characterized in that, The result output module generates a future boundary envelope using the future center trajectory and the shape parameters. Taking the bearing clearance center as the origin, it calculates the radial distance difference between the farthest and nearest points of the future boundary envelope and outputs this difference as the maximum radial runout amplitude. This includes: Construct a rotation matrix using the future principal axis direction angle; The boundary coordinates are generated by combining the parameter angle, the future principal axis half length, and the future secondary axis half length. The boundary coordinates are rotated using the rotation matrix, and the rotation result is added to the future center trajectory to obtain the future boundary points under each extrapolation step. The future boundary points under all extrapolation steps are set and merged to obtain the future boundary point set, and the convex hull of the future boundary point set is calculated, and the convex hull is used as the future boundary envelope. Using the calibrated bearing clearance center as the origin, calculate the L2 norm of all points on the future boundary envelope to the origin, and use it as the radial distance; The maximum and minimum values ​​of the radial distances are selected, and the point corresponding to the maximum value is taken as the farthest point, and the point corresponding to the minimum value is taken as the nearest point. Subtract the minimum value from the maximum value to obtain the maximum radial runout amplitude and output it.

Images