Method and system for diagnosis prediction and intelligent monitoring of electromechanical equipment

Abstract

This invention relates to the field of electromechanical equipment diagnosis and prediction technology, and particularly to a method and system for electromechanical equipment diagnosis, prediction, and intelligent monitoring. It extracts the time delay and attenuation coefficients of vibration signals at multiple locations using a cross-correlation function to construct a spatial transfer matrix. Through sparse optimization, it separates the source signals of bearings, gears, and motors, and extracts the physical layer feature vectors of each source signal. The physical layer features are organized into feature trajectories over time, and a dynamic time warping algorithm is used to match historical trajectories to obtain a time scaling factor. Based on the corresponding physical degradation equation, this factor is used to adjust the time scale coefficient. The degradation state and equation parameters are estimated using a physically constrained extended Kalman filter to calculate the remaining lifetime distribution. A causal directed acyclic graph containing load wear observation features and lifetime nodes is constructed. Maintenance interventions are performed to calculate the remaining lifetime extension. Combined with cost considerations, the optimal operation and maintenance strategy is selected, achieving high-precision remaining lifetime prediction and optimized maintenance decision-making.

Description

Technical Field

[0001] This invention relates to the field of electromechanical equipment diagnosis and prediction technology, and in particular to electromechanical equipment diagnosis, prediction and intelligent monitoring methods and systems. Background Technology

[0002] In the field of condition monitoring and fault diagnosis of electromechanical equipment, existing technologies typically use vibration signals from a single sensor for time-frequency analysis. These methods rely on empirical selection of sensor placement locations, ignoring the spatial correlation between different measuring points, resulting in poor signal source separation. Conventional practices treat the vibrations of components such as bearings, gears, and motors as a whole signal, making it difficult to effectively distinguish the characteristics of multi-source coupled faults.

[0003] Existing health management methods often set alarm limits through threshold discrimination or simple trend extrapolation, lacking modeling of the physical mechanisms of the degradation process. The construction of degradation trajectories is mostly based on the statistical characteristics of fixed time windows, without considering the nonlinear scaling of time scales under different operating conditions, resulting in large deviations in remaining life prediction.

[0004] The main drawbacks of these conventional practices are: the time delay and attenuation effects of multi-source vibration signals in spatial propagation are ignored; the feature extraction accuracy is limited after the separation of mixed signals; and misjudgment is easily caused by noise or mode aliasing. Degradation modeling lacks a joint filtering mechanism of physical constraints and state observation; statistical fitting alone is difficult to adapt to changes in working conditions; the uncertainty of remaining life prediction is difficult to quantify; and maintenance strategies cannot dynamically balance extending life and execution costs. Summary of the Invention

[0005] The present invention provides a method and system for diagnostic prediction and intelligent monitoring of electromechanical equipment, which can solve the problems in the prior art.

[0006] A first aspect of the present invention provides a method for diagnostic prediction and intelligent monitoring of electromechanical equipment, comprising:

[0007] Vibration signals from multiple spatial locations of electromechanical equipment are acquired, and the time delay and attenuation coefficients between locations are extracted using a cross-correlation function to construct a spatial transfer matrix;

[0008] Based on the spatial transfer matrix, sparse optimization separation is performed to obtain bearing source signals, gear source signals, and motor source signals, and physical layer feature vectors of each source signal are extracted;

[0009] The physical layer feature vectors are organized into feature trajectories according to time. A dynamic time warping algorithm is used to match the historical trajectory database to identify reference degraded trajectories. The time scaling factor of the current trajectory relative to the reference degraded trajectory is calculated.

[0010] The state evolution relationship is established based on the physical degradation equation corresponding to the reference degradation trajectory. The time scale coefficient is adjusted using the time scaling factor. The degradation state and equation parameters are estimated using a physically constrained extended Kalman filter. The remaining lifetime distribution is obtained by sampling the parameter distribution.

[0011] Construct a causal directed acyclic graph containing load, wear, observation features, and remaining lifetime nodes. Traverse the preset set of maintenance actions, perform intervention calculations on the causal directed acyclic graph for each maintenance action to obtain the remaining lifetime extension, and select the maintenance action with the largest net benefit as the operation and maintenance strategy based on the execution cost.

[0012] Vibration signals from multiple spatial locations of electromechanical equipment are acquired, and the time delay and attenuation coefficients between locations are extracted using a cross-correlation function to construct a spatial transfer matrix, including:

[0013] Vibration sensors are arranged in the bearing housing, gearbox housing, and motor end cover of the electromechanical equipment to collect vibration acceleration signals. After time-domain alignment of the vibration acceleration signals collected at each location, the cross-correlation function of the vibration acceleration signals between any two locations is calculated. The signal propagation delay is determined by identifying the peak position of the cross-correlation function, and the signal propagation attenuation coefficient is determined by the ratio of the peak amplitude of the cross-correlation function to the normalization coefficient.

[0014] The propagation delay between each position pair is constructed into a delay matrix, and the signal propagation attenuation coefficient between each position pair is constructed into an attenuation matrix. By applying an exponential decay kernel function transformation to the delay matrix and multiplying it element by element with the attenuation matrix, a spatial transmission matrix describing the propagation characteristics of the vibration signal in space is obtained.

[0015] Based on the spatial transfer matrix, sparse optimization separation is performed to obtain bearing source signals, gear source signals, and motor source signals. Physical layer feature vectors of each source signal are extracted, including:

[0016] The inverse of the spatial transfer matrix is ​​used as the initial estimate of the unmixing matrix. During the sparse optimization iteration process, orthogonality constraints and energy conservation constraints are applied to the unmixing matrix so that the separated bearing source signal, gear source signal and motor source signal satisfy the physical propagation reversibility condition.

[0017] The bearing source signal is subjected to minimum entropy deconvolution to obtain an impact pulse sequence. By calculating the time interval between adjacent impact events in the impact pulse sequence, periodic patterns corresponding to the inner ring fault frequency and outer ring fault frequency of the bearing are identified, and the stability index of the periodic pattern is extracted as the bearing degradation feature. Frequency domain features and phase features are extracted from the gear source signal and the motor source signal, respectively.

[0018] A time consistency verification mechanism for cross-source signals is constructed. By calculating the correlation of the changing trends of degradation features of each source signal, when the changing trend of the bearing degradation features deviates in the opposite direction from the changing trend of the frequency domain features of the gear source signal, it is determined that there is crosstalk in the source separation process and the adaptive calibration of the spatial transfer matrix is ​​triggered. The calibrated spatial transfer matrix is ​​reused for sparse optimization separation until the changing trends of degradation features of each source signal satisfy physical causal consistency, and the physical layer feature vector is obtained.

[0019] Construct a timing consistency verification mechanism for cross-source signals, including:

[0020] Calculate the rate of change of bearing degradation characteristics, gear source signal frequency domain characteristics, and motor source signal phase characteristics within a continuous time window. Construct a rate of change vector sequence from the rate of change of each source signal characteristic. Calculate the covariance matrix of the rate of change vector sequence. Extract the eigenvector corresponding to the largest eigenvalue of the covariance matrix as the main degradation direction.

[0021] Calculate the projection component of the rate of change vector of each source signal feature in the main degradation direction. When the projection component of the bearing degradation feature and the projection component of the frequency domain feature of the gear source signal have opposite signs and the difference in absolute value exceeds a preset difference threshold, it is determined that there is crosstalk in the source separation.

[0022] By calculating the change in the condition number of the spatial transfer matrix before and after the crosstalk occurs, the position pairs that cause crosstalk are identified. Time-domain fine-tuning is applied to the propagation delay corresponding to the position pair, and frequency-domain compensation is applied to the attenuation coefficient corresponding to the position pair to obtain the calibrated spatial transfer matrix.

[0023] The inverse of the calibrated spatial transfer matrix is ​​used as the new unmixing matrix to re-perform sparse optimization separation. The temporal consistency check is repeated until the projected components of each source signal feature have the same sign and the condition number of the covariance matrix is ​​lower than the preset condition number threshold.

[0024] The physical layer feature vectors are organized into feature trajectories according to time. A dynamic time warping algorithm is used to match historical trajectory databases to identify reference degraded trajectories. The time scaling factor of the current trajectory relative to the reference degraded trajectory is calculated, including:

[0025] The physical layer feature vectors are arranged in timestamp order to form the current feature trajectory, and historical trajectories marked with degradation stage divisions are selected from the historical trajectory database as candidate reference trajectories;

[0026] A dynamic time warping algorithm is used to calculate the warping path between the current feature trajectory and each candidate reference trajectory. When constructing the cumulative distance matrix, time monotonicity constraints and slope smoothness constraints are applied to the warping path. The time monotonicity constraint requires that the time index pairs of the warping path be monotonically increasing, and the slope smoothness constraint limits the slope change of adjacent time index pairs of the warping path. The candidate reference trajectory with the smallest cumulative distance that simultaneously satisfies the time monotonicity constraint and the slope smoothness constraint is selected as the reference degenerate trajectory.

[0027] The reference degradation trajectory is divided into an initial degradation segment, a stable degradation segment, and an accelerated degradation segment according to the degradation stage division points marked on it. The mapping relationship of the curved path in the time dimension is extracted, and the time scaling factor of the current feature trajectory relative to the reference degradation trajectory in the corresponding interval of each degradation stage is calculated respectively.

[0028] The confidence index of the time scaling factor is constructed based on the cumulative distance of the curved path. When the cumulative distance exceeds a preset distance threshold, the historical trajectory library is re-searched to update the reference degraded trajectory.

[0029] A state evolution relationship is established based on the physical degradation equation corresponding to the reference degradation trajectory. The time scale coefficient is adjusted using the time scaling factor. The degradation state and equation parameters are estimated using a physically constrained extended Kalman filter. The remaining lifetime distribution is obtained by sampling the parameter distribution, including:

[0030] Based on the type of electromechanical equipment corresponding to the reference degradation trajectory, a physical degradation equation is selected, and the physical degradation equation is expressed in state-space form, where the state variable is the degree of degradation and the observation variable is the physical layer feature vector;

[0031] The time scaling factor is used as the scaling coefficient of the time variable in the physical degradation equation. The corresponding time scaling factor is applied to adjust the time scale for the initial degradation segment, the stable degradation segment, and the accelerated degradation segment, respectively.

[0032] An extended Kalman filter framework is constructed, which uses the physical degradation equation as the state transition equation, imposes degradation monotonicity constraints on the state variables, and imposes physical feasible region constraints on the equation parameters;

[0033] The confidence index of the time scaling factor is mapped to the adjustment coefficient of the observation noise covariance matrix. When the cumulative distance corresponding to the confidence index is lower than a first preset threshold, the observation noise covariance is reduced; when the cumulative distance exceeds the first distance threshold, the observation noise covariance is increased. The posterior distribution of the equation parameters is obtained through extended Kalman filtering iteration.

[0034] A parameter sample set is generated by sampling from the posterior distribution of the equation parameters. For each parameter sample, the physical degradation equation is extrapolated to the failure threshold to obtain the remaining lifetime sample, and the remaining lifetime distribution is statistically obtained.

[0035] A causal directed acyclic graph (DAG) is constructed, containing load, wear, observed features, and remaining lifetime nodes. A pre-defined set of maintenance actions is traversed, and for each maintenance action, intervention calculations are performed on the causal DAG to obtain the remaining lifetime extension. Based on the execution cost, the maintenance action with the largest net benefit is selected as the operational strategy, including:

[0036] Calculate the partial derivatives of the state transition matrix of the physical degradation equation with respect to the load variable, wear variable, and observation variable. Use the variable pairs corresponding to the non-zero elements in the partial derivative matrix as directed edges of a causal directed acyclic graph. Use the posterior distribution of the equation parameters as the weight distribution of the directed edges to construct a causal directed acyclic graph with load as the root node, wear amount as the intermediate node, physical layer feature vector as the observation node, and remaining lifetime distribution as the terminal node.

[0037] Maintenance actions are selected sequentially from a preset set of maintenance actions, including load adjustment actions and component replacement actions.

[0038] For load adjustment actions, the incoming edges pointing to the load nodes are deleted and the load node value is fixed to the post-intervention load level. The weight distribution of the edges is used to propagate forward along the directed edges to the remaining lifetime nodes to obtain the post-intervention remaining lifetime distribution. The monotonicity constraint of the physical degradation equation and the failure threshold constraint are applied to the distribution for verification. The difference between the expected value of the post-intervention remaining lifetime distribution that passes the verification and the expected value of the remaining lifetime distribution is taken as the remaining lifetime extension.

[0039] For component replacement actions, delete the incoming edges pointing to the wear node and reset the wear node to zero wear state. Propagate along the directed edges to calculate and verify the remaining lifetime distribution after intervention, and calculate the remaining lifetime extension.

[0040] Calculate the net benefit of each maintenance action and select the maintenance action with the largest net benefit as the operation and maintenance strategy.

[0041] A second aspect of the present invention provides a diagnostic prediction and intelligent monitoring system for electromechanical equipment, comprising:

[0042] The spatial matrix unit is used to acquire vibration signals from multiple spatial locations of electromechanical equipment. The time delay and attenuation coefficients between locations are extracted through the cross-correlation function to construct the spatial transfer matrix.

[0043] The source separation feature unit is used to perform sparse optimization separation based on the spatial transfer matrix to obtain bearing source signals, gear source signals and motor source signals, and extract the physical layer feature vector of each source signal;

[0044] The trajectory matching unit is used to organize the physical layer feature vectors into feature trajectories according to time, use a dynamic time warping algorithm to match the historical trajectory library to identify reference degraded trajectories, and calculate the time scaling factor of the current trajectory relative to the reference degraded trajectory.

[0045] The degradation estimation unit is used to establish the state evolution relationship based on the physical degradation equation corresponding to the reference degradation trajectory, adjust the time scale coefficient using the time scaling factor, estimate the degradation state and equation parameters using a physically constrained extended Kalman filter, and obtain the remaining lifetime distribution by sampling on the parameter distribution.

[0046] The operation and maintenance decision unit is used to construct a causal directed acyclic graph containing load, wear, observation characteristics and remaining lifetime nodes. It traverses the preset set of maintenance actions, performs intervention calculations on the causal directed acyclic graph for each maintenance action to obtain the remaining lifetime extension, and selects the maintenance action with the largest net benefit as the operation and maintenance strategy based on the execution cost.

[0047] A third aspect of the present invention provides an electronic device, comprising:

[0048] processor;

[0049] Memory used to store processor-executable instructions;

[0050] The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0051] A fourth aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0052] This invention accurately separates the source signals of bearings, gears and motors from the mixed vibration signals of complex electromechanical equipment, extracts physical layer feature vectors, effectively suppresses multi-source interference and noise effects, greatly improves the accuracy of fault location and the ability to identify early weak faults, and provides pure input with clear physical meaning for subsequent degradation analysis.

[0053] The dynamic time warp algorithm matches the historical trajectory library and identifies the reference degradation trajectory. It adaptively adjusts the time scale through the time scaling factor, eliminating the interference of operating condition fluctuations and individual differences on the degradation process. Combined with the physical degradation equation and the extended Kalman filter with physical constraints, the remaining life estimation not only follows the data change law, but also meets the inherent physical evolution constraints of the mechanical system. The prediction results are more accurate, interpretable and robust.

[0054] A causal directed acyclic graph (DAG) is constructed, incorporating load, wear, observed features, and remaining lifetime nodes. Intervention calculations are performed on each candidate maintenance action to quantify its extension of remaining lifetime. The net benefit is then ranked based on the overall execution cost, and the optimal maintenance strategy is automatically output. This method achieves end-to-end intelligent management from condition monitoring and lifetime prediction to maintenance decision-making, significantly reducing the probability of unplanned downtime, minimizing blind and excessive maintenance, effectively extending equipment service life, and reducing maintenance costs. Attached Figure Description

[0055] Figure 1 A flowchart illustrating the diagnostic prediction and intelligent monitoring methods for electromechanical equipment;

[0056] Figure 2 Remaining lifetime distribution map based on Monte Carlo sampling. Detailed Implementation

[0057] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0058] The technical solution of the present invention will be described in detail below with reference to specific embodiments. These specific embodiments can be combined with each other, and the same or similar concepts or processes may not be described again in some embodiments.

[0059] Figure 1 This is a flowchart illustrating the electromechanical equipment diagnosis, prediction, and intelligent monitoring method according to an embodiment of the present invention.

[0060] Methods for diagnostic prediction and intelligent monitoring of electromechanical equipment include:

[0061] Vibration signals from multiple spatial locations of electromechanical equipment are acquired, and the time delay and attenuation coefficients between locations are extracted using a cross-correlation function to construct a spatial transfer matrix;

[0062] Based on the spatial transfer matrix, sparse optimization separation is performed to obtain bearing source signals, gear source signals, and motor source signals, and physical layer feature vectors of each source signal are extracted;

[0063] The physical layer feature vectors are organized into feature trajectories according to time. A dynamic time warping algorithm is used to match the historical trajectory database to identify reference degraded trajectories. The time scaling factor of the current trajectory relative to the reference degraded trajectory is calculated.

[0064] The state evolution relationship is established based on the physical degradation equation corresponding to the reference degradation trajectory. The time scale coefficient is adjusted using the time scaling factor. The degradation state and equation parameters are estimated using a physically constrained extended Kalman filter. The remaining lifetime distribution is obtained by sampling the parameter distribution.

[0065] Construct a causal directed acyclic graph containing load, wear, observation features, and remaining lifetime nodes. Traverse the preset set of maintenance actions, perform intervention calculations on the causal directed acyclic graph for each maintenance action to obtain the remaining lifetime extension, and select the maintenance action with the largest net benefit as the operation and maintenance strategy based on the execution cost.

[0066] Vibration signals from multiple spatial locations of electromechanical equipment are acquired, and the time delay and attenuation coefficients between locations are extracted using a cross-correlation function to construct a spatial transfer matrix, including:

[0067] Vibration sensors are arranged in the bearing housing, gearbox housing, and motor end cover of the electromechanical equipment to collect vibration acceleration signals. After time-domain alignment of the vibration acceleration signals collected at each location, the cross-correlation function of the vibration acceleration signals between any two locations is calculated. The signal propagation delay is determined by identifying the peak position of the cross-correlation function, and the signal propagation attenuation coefficient is determined by the ratio of the peak amplitude of the cross-correlation function to the normalization coefficient.

[0068] The propagation delay between each position pair is constructed into a delay matrix, and the signal propagation attenuation coefficient between each position pair is constructed into an attenuation matrix. By applying an exponential decay kernel function transformation to the delay matrix and multiplying it element by element with the attenuation matrix, a spatial transmission matrix describing the propagation characteristics of the vibration signal in space is obtained.

[0069] In practical engineering scenarios for monitoring the operational status of electromechanical equipment, the placement of vibration signal acquisition points has a decisive impact on the accuracy of subsequent signal separation and degradation analysis. Based on the structural characteristics of typical electromechanical equipment, acceleration-type vibration sensors are installed at three key locations: bearing housings, gearbox housings, and motor end covers. Sensors at the bearing housings are preferentially positioned near the rolling element load area to maximize the capture of impact vibration components from the inner and outer rings of the bearing and the rolling elements. Sensors at the gearbox housings are installed near the rigid ribs of the housing to avoid signal distortion due to localized flexibility. Sensors at the motor end covers are positioned near the center of the end cover to acquire vibration signals generated by electromagnetic excitations such as rotor eccentricity and uneven air gap. The vibration acceleration signals acquired by each sensor are time-domain aligned after being constrained by a synchronous sampling clock. Alignment can be achieved using hardware synchronization based on reference pulse triggering or software alignment based on cross-correlation peak search, ensuring that signals from different locations have a unified reference benchmark on the time axis and eliminating systematic time differences introduced by inherent delays between acquisition card channels.

[0070] After completing temporal alignment, for any two sensor positions and Calculate the cross-correlation function of the vibration acceleration signals between them. Let the position... The vibration acceleration signal is ,Location The vibration acceleration signal is The cross-correlation function between the two is defined as follows: ,in This represents the time delay variable. In discrete implementations, the Fast Fourier Transform (FFT) is used to accelerate cross-correlation calculations for finite-length signal sequences. The time-domain convolution is transformed into a frequency-domain multiplication, followed by an inverse transform to obtain the cross-correlation sequence, significantly reducing computational complexity. Cross-correlation function. The peak position corresponds to the vibration signal from position. Spread to location The time delay experienced is denoted as the time delay corresponding to the peak position. ,Right now Peak location identification employs parabolic interpolation to fit three points near the discrete peak, thereby obtaining a sub-sampling precision time delay estimate and improving time delay resolution.

[0071] The signal propagation attenuation coefficient is determined based on the ratio of the peak amplitude of the cross-correlation function to the normalization coefficient. Specifically, the normalization coefficient is taken as the geometric mean of the amplitudes of the autocorrelation functions of the two signals at zero time delay, i.e. ,in and Positions With position The value of the autocorrelation function of the signal at zero time delay. This ratio. The physical meaning is that vibrational energy travels from position... Spread to location The normalized attenuation level during the process ranges from 0 to 1. A value closer to 1 indicates smaller vibration propagation loss between the two positions and a tighter structural coupling; a value closer to 0 indicates greater energy dissipation along the propagation path, which may originate from factors such as material damping, contact interface friction, or geometric diffusion. After performing cross-correlation calculations on all position pairs sequentially, a complete set of time delays and attenuation coefficients between positions can be obtained.

[0072] The propagation delay between each position pair Arrange by position index and construct the delay matrix. , where the matrix is ​​numbered Line number The elements of the column are For the same position itself, the time delay is defined as zero, that is, the diagonal elements of the matrix are zero. Time Delay Matrix This reflects the temporal relationship of vibration signal propagation between sensor nodes and is one of the core inputs for constructing the spatial transfer matrix. Similarly, the attenuation coefficients between each pair of locations... Arrange the elements according to the same indexing rule to construct the attenuation matrix. The diagonal elements of the matrix are all 1, indicating that the signal has no attenuation at its current location. Attenuation matrix The amplitude attenuation distribution characteristics of vibration energy during its propagation in space are described.

[0073] The spatial transfer matrix is ​​constructed by applying an exponential decay kernel function to the time delay matrix and then multiplying it element-by-element by the decay matrix. The introduction of the exponential decay kernel function has a clear physical basis: when a vibration signal propagates in an elastic medium, its frequency components attenuate to varying degrees with increasing propagation distance, with high-frequency components attenuating faster than low-frequency components. This characteristic can be described by an exponential frequency-distance attenuation relationship. The specific operation of applying the exponential decay kernel function to the time delay matrix is ​​as follows: for each element in the matrix... calculate ,in The attenuation rate parameter can be calibrated based on the damping characteristics of the equipment's structural materials, or determined by fitting historical data under known degradation states. The kernel matrix is ​​obtained after transformation. Its elements This reflects the exponential decay contribution of time delay to the energy propagation of vibration.

[0074] kernel matrix With decay matrix Element-wise multiplication yields the spatial transfer matrix. ,Right now Spatial transfer matrix Each element of the spatial transfer matrix integrates information on both propagation delay and energy attenuation, comprehensively characterizing the amplitude-frequency and temporal characteristics of the vibration signal propagating from one spatial location to another. Diagonal elements are 1, indicating distortion-free transmission of the signal from each sensor location; the magnitudes of off-diagonal elements quantify the strength of vibration coupling between locations. In practical calculations, due to the limited number of sensors, the dimension of the spatial transfer matrix is ​​usually consistent with the number of sensors. For example, in the case where one sensor is placed in each of the bearing housing, gearbox housing, and motor end cover, the spatial transfer matrix is ​​a 3×3 real-valued matrix.

[0075] Completed spatial transfer matrix In subsequent steps, this serves as the core input constraint for sparse optimization separation. By combining the observed multi-location mixed vibration signals with the spatial transfer matrix, a mixed observation model is established, and then the signal components of each source are solved under sparse prior constraints. The accuracy of the spatial transfer matrix directly determines the separation quality of the bearing source signal, gear source signal, and motor source signal. Therefore, when the operating conditions of the equipment change significantly (such as a step change in speed or a sudden change in load), it is necessary to update the spatial transfer matrix online, recalculate the cross-correlation function, and refresh the time delay matrix and attenuation matrix to ensure that the spatial transfer matrix always reflects the actual vibration propagation characteristics under the current operating conditions.

[0076] Based on the spatial transfer matrix, sparse optimization separation is performed to obtain bearing source signals, gear source signals, and motor source signals. Physical layer feature vectors of each source signal are extracted, including:

[0077] The inverse of the spatial transfer matrix is ​​used as the initial estimate of the unmixing matrix. During the sparse optimization iteration process, orthogonality constraints and energy conservation constraints are applied to the unmixing matrix so that the separated bearing source signal, gear source signal and motor source signal satisfy the physical propagation reversibility condition.

[0078] The bearing source signal is subjected to minimum entropy deconvolution to obtain an impact pulse sequence. By calculating the time interval between adjacent impact events in the impact pulse sequence, periodic patterns corresponding to the bearing inner ring fault frequency and outer ring fault frequency are identified, and the stability index of the periodic pattern is extracted as a bearing degradation feature.

[0079] Frequency domain features and phase features are extracted from the gear source signal and the motor source signal, respectively.

[0080] A time consistency verification mechanism for cross-source signals is constructed. By calculating the correlation of the changing trends of degradation features of each source signal, when the changing trend of the bearing degradation features deviates in the opposite direction from the changing trend of the frequency domain features of the gear source signal, it is determined that there is crosstalk in the source separation process and the adaptive calibration of the spatial transfer matrix is ​​triggered. The calibrated spatial transfer matrix is ​​reused for sparse optimization separation until the changing trends of degradation features of each source signal satisfy physical causal consistency, and the physical layer feature vector is obtained.

[0081] In obtaining the spatial transfer matrix Then, inverse its matrix As an initial estimate of the unmixing matrix, this serves as the starting point for sparse optimization separation. Using the inverse process of the physical propagation relationship directly as the initial value allows the iterative process to unfold within a physically reasonable parameter space, avoiding local extremum traps caused by random initialization. During the sparse optimization iteration process, the unmixing matrix... Simultaneously apply orthogonality constraints and energy conservation constraints. Orthogonality constraints require... Approaching the identity matrix ensures that the signals from different sources are linearly independent after separation, preventing signal components from interpenetrating. The energy conservation constraint requires that the total energy of the signals before and after separation remain unchanged; that is, the sum of the energy of each observation channel is statistically equal to the sum of the estimated energy of each source signal. These two constraints together guarantee that the separated bearing source signal, gear source signal, and motor source signal satisfy the physical propagation reversibility condition, meaning that the original observed signal can be reconstructed through the inverse transformation of the unmixing matrix, and the error with the actual acquired data is controlled within a preset threshold. The objective function of sparse optimization introduces a parameter beyond the data fitting term. The regularization term causes each source signal to have sparse excitation characteristics in the time domain, which conforms to the physical excitation mechanisms of bearing impact, gear meshing and motor harmonics.

[0082] The separated bearing source signal is subjected to minimum entropy deconvolution processing. Minimum entropy deconvolution maximizes the kurtosis of the filtered output signal by designing a finite impulse response filter, thereby recovering the impact pulse sequence hidden in the background noise and transmission path response. After obtaining the impact pulse sequence, the time interval sequence between adjacent impact events is calculated, denoted as... ,in This serves as an index for the impact event sequence. Based on the bearing geometry and rotational speed, the characteristic frequency of the bearing inner ring failure is pre-calculated. and outer ring fault characteristic frequency The corresponding theoretical impact periods are respectively and The measured time interval sequence By matching with theoretical periods, periodic patterns corresponding to the failure frequencies of the inner or outer races are identified. Stability indices are calculated by determining the standard deviation of the measured interval sequence relative to the theoretical period. To quantify this, a smaller deviation indicates more stable periodicity, corresponding to a deeper degree of bearing degradation and more regular fault characteristics. The slope of the bearing's degradation characteristics over time is incorporated into the bearing degradation feature vector and used as input for subsequent degradation state estimation.

[0083] For the gear source signal, extract the meshing frequency in the frequency domain. The amplitude and phase of the gear's harmonics, along with the width and amplitude ratio of its sidebands, constitute the gear's frequency domain feature vector. The appearance and expansion of sidebands are usually associated with degradation modes such as tooth surface wear, pitting, or tooth breakage, and their amplitude growth rate reflects the gear's degradation speed. For the motor source signal, the phase stability of the rotational frequency and its harmonics, as well as the amplitude of characteristic frequency components in the stator current spectrum related to rotor eccentricity and winding asymmetry, are extracted to form the motor's phase feature vector. The bearing degradation features, gear frequency domain features, and motor phase features are concatenated to form a complete physical layer feature vector for subsequent degradation trajectory analysis.

[0084] After obtaining the degradation characteristics of each source signal, a time-series consistency verification mechanism for cross-source signals is constructed. Linear trends are fitted to the bearing degradation characteristic sequence and the gear frequency domain characteristic sequence within a sliding window to obtain their respective trend slopes. and Under normal operating conditions, the overall degradation process of electromechanical equipment exhibits physical causal consistency; that is, the degradation characteristics of each component should show a trend of change in the same direction, and the product of their slopes should be positive. When detected... and A reverse deviation occurs, that is (in When the preset deviation judgment threshold is reached, crosstalk occurs in the judgment source separation process, and component aliasing occurs between the bearing source signal and the gear source signal. The current demixing matrix can no longer accurately reflect the actual propagation path.

[0085] After triggering crosstalk detection, an adaptive calibration process for the spatial transfer matrix is ​​initiated. The calibration process utilizes the newly acquired multi-location vibration signals at the current time to recalculate the cross-correlation function between each pair of locations and update the time delay matrix. With decay matrix Thus, the kernel matrix is ​​reconstructed. and spatial transfer matrix The calibrated spatial transfer matrix reflects the actual propagation characteristics of the vibration signal in the current degraded state of the equipment. This is because as bearing wear intensifies or gear meshing stiffness changes, the dynamic transmission path of the structure will drift, causing deviations in the original transfer matrix parameters. (The calibrated spatial transfer matrix is ​​used as a reference.) The unmixing matrix is ​​reinitialized, and the sparse optimization iteration is re-executed under the conditions of applying orthogonality constraints and energy conservation constraints to obtain a new round of source signal separation results.

[0086] After re-separation, the slope of the degradation characteristics of each source signal is calculated again to verify... Whether it recovers to a positive value and the absolute value exceeds If the physical causal consistency condition is still not met, calibration is triggered again, and the above process is repeated until the changing trends of the degradation characteristics of each source signal meet the physical causal consistency condition. In practical engineering scenarios, after 1 to 3 calibration iterations, the source separation results can usually converge to a physically consistent state. The final output bearing degradation characteristics, gear frequency domain characteristics, and motor phase characteristics together constitute the physical layer feature vector. This vector accurately reflects the true degradation state of each component, providing reliable input data for subsequent degradation trajectory matching and remaining life prediction based on the dynamic time bending algorithm.

[0087] Construct a timing consistency verification mechanism for cross-source signals, including:

[0088] Calculate the rate of change of bearing degradation characteristics, gear source signal frequency domain characteristics, and motor source signal phase characteristics within a continuous time window. Construct a rate of change vector sequence from the rate of change of each source signal characteristic. Calculate the covariance matrix of the rate of change vector sequence. Extract the eigenvector corresponding to the largest eigenvalue of the covariance matrix as the main degradation direction.

[0089] Calculate the projection component of the rate of change vector of each source signal feature in the main degradation direction. When the projection component of the bearing degradation feature and the projection component of the frequency domain feature of the gear source signal have opposite signs and the difference in absolute value exceeds a preset difference threshold, it is determined that there is crosstalk in the source separation.

[0090] By calculating the change in the condition number of the spatial transfer matrix before and after the crosstalk occurs, the position pairs that cause crosstalk are identified. Time-domain fine-tuning is applied to the propagation delay corresponding to the position pair, and frequency-domain compensation is applied to the attenuation coefficient corresponding to the position pair to obtain the calibrated spatial transfer matrix.

[0091] The inverse of the calibrated spatial transfer matrix is ​​used as the new unmixing matrix to re-perform sparse optimization separation. The temporal consistency check is repeated until the projected components of each source signal feature have the same sign and the condition number of the covariance matrix is ​​lower than the preset condition number threshold.

[0092] In actual electromechanical equipment operation, the three types of source signals—bearings, gears, and motors—have physical coupling paths. When the estimation accuracy of the spatial transfer matrix deviates, energy crosstalk occurs among the source signals after sparse optimization separation, leading to distortion in subsequent degradation feature extraction. To address this, a temporal consistency verification mechanism for cross-source signals is constructed. This mechanism uses joint analysis of multi-source features to detect and correct this crosstalk, ensuring the quality of source separation.

[0093] Within the continuous acquisition time window, bearing degradation features, gear source signal frequency domain features, and motor source signal phase features are extracted respectively, and the rate of change of each feature between adjacent windows is calculated. Specifically, let the rate of change of bearing degradation features be... The frequency domain characteristic change rate of the gear is The motor phase characteristic change rate is The three in the The time windows constitute the rate of change vector. In accumulation After several consecutive time windows, the rate of change vector sequence is obtained. Calculate the sample covariance matrix for this sequence. ,in The Line number The column elements reflect the first The rate of change of source signal characteristics and the first The linear covariant relationship between the rate of change of characteristics of similar source signals. Perform eigenvalue decomposition and extract the largest eigenvalue. corresponding feature vector As the main degradation direction, this main degradation direction physically represents the comprehensive direction where energy changes are most concentrated during the degradation process of multi-source signals. Under normal degradation conditions, the rate of change of each source signal should maintain a projection component with the same sign in this direction, that is, the degradation trends of various components evolve together macroscopically.

[0094] The rate of change vector for each time window Projected onto the principal degradation direction The scalar projection components of the source signal characteristics are obtained above. Let the bearing degradation characteristics be at the [missing information]. The projection component of the window is The contribution of the bearing component is given by the projection component of the gear frequency domain characteristics. .when and The signs are opposite, and the difference in their absolute values ​​is... Exceeding the preset difference threshold At that moment, it is determined that crosstalk exists in the source separation. Opposite signs indicate a divergence in the evolution trends of bearing degradation characteristics and gear degradation characteristics along the main degradation direction, which is physically unreasonable—components of the same unit should degrade in a coordinated direction under similar operating conditions; a significant reversal indicates energy leakage or aliasing in the separation results. Difference threshold. The standard deviation is determined based on the statistical distribution of historical normal operation data, and is usually three times the standard deviation of the difference in the projected components under normal conditions, in order to avoid misjudgment caused by noise.

[0095] After crosstalk is detected, the sensor pair causing the crosstalk is located. The spatial transfer matrix condition number is calculated at each time step before and after the crosstalk occurrence. The condition number reflects the numerical stability of the transfer matrix; a sudden increase in the condition number indicates that the propagation parameter estimates between corresponding position pairs have become singular or ill-conditioned. This can be analyzed by comparing the conditions before and after the crosstalk time. The change in condition number is used to identify the position pairs with the largest changes as the root cause of crosstalk, denoted as position pairs. The processing of this position pair involves two dimensions: first, applying time-domain fine-tuning to its propagation delay, i.e., searching for the correction amount that minimizes the condition number based on the original estimated delay with small steps, and then substituting the corrected delay into the delay matrix. The update is complete; secondly, frequency domain compensation is applied to its attenuation coefficient. The amplitude-frequency response deviation is analyzed in the frequency domain based on the cross-correlation function of the position pair, compensation weights are applied to each frequency band, and the attenuation matrix is ​​corrected. The corresponding elements in the equation make the frequency domain energy allocation more consistent with the laws of physical propagation. The time domain fine-tuning of the time delay is usually searched within ±5% of the original estimate in steps of 0.1 sampling intervals; the frequency domain compensation is weighted and corrected based on the energy ratio deviation of the frequency band near the fault characteristic frequency to ensure that the compensation does not introduce new frequency band imbalances.

[0096] After correcting for time delay and attenuation coefficients, the kernel matrix is ​​recalculated according to the kernel matrix transformation rules. This leads to an update of the spatial transfer matrix. The calibrated spatial transfer matrix is ​​obtained. .by The inverse matrix is ​​used as the new unmixing matrix. The multi-channel vibration observation signals were re-sparsely optimized and separated to obtain updated bearing source signals, gear source signals, and motor source signals. During the sparse optimization separation process, the sparsity (i.e., the proportion of non-zero coefficients) of each source signal was used as the optimization objective. The source signal estimates that maximize the sparsity of the separation results are obtained under constraints, ensuring the independence of each source signal in the time and frequency domain after separation.

[0097] After separation, extract the features of each source signal again and repeat the above time series consistency check process: recalculate the rate of change vector sequence and update the covariance matrix. The new principal degradation direction is extracted, the projection components of each source signal are calculated, and the sign consistency and difference are judged to meet the requirements. Simultaneously, the condition number of the updated covariance matrix is ​​checked. Is it below the preset threshold number? A low condition number in the covariance matrix indicates a stable correlation structure between the characteristic rates of change of each source signal, uniform energy distribution in all directions, and no numerical ill-conditioning caused by extreme energy concentration in any one direction. If inconsistencies still exist in the signs of the projected components, or Then, the crosstalk position pairs are repositioned, time delay fine-tuning and frequency domain compensation are performed again, the spatial transfer matrix is ​​iteratively updated and the signals are separated again, until the projected component signs of all source signal features are consistent and until.

[0098] The iterative process sets a maximum upper limit on the number of iterations to prevent infinite loops caused by crosstalk that cannot be completely eliminated under extreme operating conditions. If the convergence condition is not met even after reaching the maximum number of iterations, the current optimal unmixing matrix is ​​recorded and an alarm is output, prompting maintenance personnel to check whether there are any abnormal changes in the sensor installation status or propagation path. Through the above-mentioned cross-source signal timing consistency verification mechanism, source separation errors can be dynamically corrected online, ensuring the quality of input features for subsequent degradation state estimation and remaining lifetime prediction, and significantly improving the robustness of the overall diagnostic prediction system.

[0099] The physical layer feature vectors are organized into feature trajectories according to time. A dynamic time warping algorithm is used to match historical trajectory databases to identify reference degraded trajectories. The time scaling factor of the current trajectory relative to the reference degraded trajectory is calculated, including:

[0100] The physical layer feature vectors are arranged in timestamp order to form the current feature trajectory, and historical trajectories marked with degradation stage divisions are selected from the historical trajectory database as candidate reference trajectories;

[0101] A dynamic time warping algorithm is used to calculate the warping path between the current feature trajectory and each candidate reference trajectory. When constructing the cumulative distance matrix, time monotonicity constraints and slope smoothness constraints are applied to the warping path. The time monotonicity constraint requires that the time index pairs of the warping path be monotonically increasing, and the slope smoothness constraint limits the slope change of adjacent time index pairs of the warping path. The candidate reference trajectory with the smallest cumulative distance that simultaneously satisfies the time monotonicity constraint and the slope smoothness constraint is selected as the reference degenerate trajectory.

[0102] The reference degradation trajectory is divided into an initial degradation segment, a stable degradation segment, and an accelerated degradation segment according to the degradation stage division points marked on it. The mapping relationship of the curved path in the time dimension is extracted, and the time scaling factor of the current feature trajectory relative to the reference degradation trajectory in the corresponding interval of each degradation stage is calculated respectively.

[0103] The confidence index of the time scaling factor is constructed based on the cumulative distance of the curved path. When the cumulative distance exceeds a preset distance threshold, the historical trajectory library is re-searched to update the reference degraded trajectory.

[0104] The physical layer feature vectors are arranged in timestamp order to form the current feature trajectory. The physical layer feature vectors output by each time window contain feature components of three types of source signals: bearings, gears, and motors. These are concatenated in strictly ascending order of the acquisition timestamps to form a trajectory sequence that evolves over time in a high-dimensional feature space. The historical trajectory database stores a large number of equipment degradation trajectories that have completed full lifecycle monitoring. Each historical trajectory is accompanied by degradation stage division information marked manually or by algorithms, including the start and end time nodes of the three stages: initial degradation, stable degradation, and accelerated degradation. Historical trajectories that match the current equipment type and operating load range are selected from the historical trajectory database as a candidate reference trajectory set. The size of the candidate set is usually controlled between several dozen and one hundred to ensure the efficiency and coverage of subsequent matching calculations.

[0105] The dynamic time warp algorithm needs to construct a cumulative distance matrix when calculating the similarity between the current feature trajectory and the candidate reference trajectory. Let the length of the current feature trajectory be... The length of a candidate reference trajectory is The dimension of the cumulative distance matrix is Each element in the matrix Indicates the current trajectory before Before each time step and candidate reference trajectory The optimal alignment cumulative distance between time steps, where The time step index for the current trajectory. This provides the time step index for the candidate reference trajectory. The local distance uses Euclidean distance to measure the difference in feature vectors between the two trajectories at the corresponding time step.

[0106] When constructing the cumulative distance matrix, two types of constraints are imposed on the curved path. The time monotonicity constraint requires that the time indexes traversed by the curved path be... It must satisfy the condition that it monotonically increases as the path progresses, that is, in subsequent time index pairs. and All values ​​must be at least equal to the corresponding component of the current index pair. This constraint ensures that the time mapping relationship will not result in physically unreasonable situations such as time regression. The slope smoothness constraint restricts the change in the step direction of the curved path between adjacent time index pairs. Specifically, the slope of the curved path between two adjacent steps is defined as... ,in and These represent the step size of adjacent time index pairs in the current trajectory direction and the reference trajectory direction, respectively. The slope smoothness constraint requires that the slope change within a series of consecutive steps does not exceed a preset slope change upper limit. This avoids drastic folding in the curved path, ensuring the smoothness and physical rationality of the time mapping relationship. The cumulative distance corresponding to the curved path that satisfies the above two types of constraints is the matching distance of the candidate reference trajectory. The trajectory with the smallest matching distance and simultaneously satisfies both types of constraints is selected from all candidate reference trajectories as the reference degenerate trajectory.

[0107] After determining the reference degradation trajectory, it is divided into three segments based on the marked degradation stage division points: initial degradation, stable degradation, and accelerated degradation. The initial degradation segment corresponds to the early stage when the device transitions from a normal state to a detectable degradation signal; the stable degradation segment corresponds to the middle stage where degradation accumulates slowly; and the accelerated degradation segment corresponds to the late stage where the degradation rate accelerates significantly and approaches failure. The time dimension mapping relationship is extracted from the curved path output by the dynamic time warping algorithm, that is, for the time interval in the reference degradation trajectory belonging to a certain degradation stage... Find the time interval corresponding to the current feature trajectory. Time scaling factor Defined as the ratio of the interval length corresponding to the current trajectory to the interval length corresponding to the reference degenerate trajectory, i.e. subscript This indicates the degradation stage number, corresponding to the initial, stable, and accelerated stages. When When this occurs, it indicates that the current device is evolving slower than the reference trajectory during this degradation phase, and the degradation process is relatively delayed; when This indicates that the current equipment degradation process is relatively accelerated. The time scaling factors of the three stages together constitute a set of quantitative indicators describing the current equipment degradation rate, providing a basis for subsequent timescale adjustments to the physical degradation equation.

[0108] A confidence index for constructing a time scaling factor based on the cumulative distance of the curved path. The confidence index reflects the similarity between the current feature trajectory and the reference degenerate trajectory; a higher value indicates better matching quality and stronger reliability of the time scaling factor. Specifically, it calculates the cumulative distance of the optimal match. The mean of the cumulative distances to all candidate reference trajectories in the historical trajectory database For comparison, the confidence index is defined as follows: ,when much smaller hour, Approaching 1 indicates that the reference degradation trajectory is highly representative; when near hour, A value approaching 0 indicates a lack of references in the historical trajectory database that are highly similar to the current trajectory. When the accumulated distance... Exceeding the preset distance threshold When this occurs, a re-retrieval process of the historical trajectory database is triggered. During the re-retrieval, the operating condition range of the candidate reference trajectory is expanded, or historical trajectories of different equipment models are introduced. Dynamic time warping matching is re-executed to obtain a more suitable reference degradation trajectory and update the current time scaling factor estimate. This dynamic update mechanism ensures that the reference degradation trajectory can be adjusted in a timely manner when equipment operating conditions drift or atypical degradation modes occur, avoiding systematic deviations in remaining lifetime prediction due to reference trajectory mismatch.

[0109] In practical engineering applications, maintaining the historical trajectory database is equally important. Whenever a device completes full lifecycle monitoring and records a failure event, its complete characteristic trajectory, along with the degradation stage annotation, is added to the historical trajectory database, continuously expanding the database's coverage. For newly added historical trajectories, clustering methods are used to remove redundancy from existing trajectories, retaining typical representative trajectories for each degradation mode to avoid decreased retrieval efficiency due to excessive database size. (Upper limit of slope variation) With distance threshold All of them can be adaptively adjusted according to the equipment type and historical statistical results to adapt to the differences in degradation modes of different electromechanical equipment.

[0110] A state evolution relationship is established based on the physical degradation equation corresponding to the reference degradation trajectory. The time scale coefficient is adjusted using the time scaling factor. The degradation state and equation parameters are estimated using a physically constrained extended Kalman filter. The remaining lifetime distribution is obtained by sampling the parameter distribution, including:

[0111] Based on the type of electromechanical equipment corresponding to the reference degradation trajectory, a physical degradation equation is selected, and the physical degradation equation is expressed in state-space form, where the state variable is the degree of degradation and the observation variable is the physical layer feature vector;

[0112] The time scaling factor is used as the scaling coefficient of the time variable in the physical degradation equation. The corresponding time scaling factor is applied to adjust the time scale for the initial degradation segment, the stable degradation segment, and the accelerated degradation segment, respectively.

[0113] An extended Kalman filter framework is constructed, which uses the physical degradation equation as the state transition equation, imposes degradation monotonicity constraints on the state variables, and imposes physical feasible region constraints on the equation parameters;

[0114] The confidence index of the time scaling factor is mapped to the adjustment coefficient of the observation noise covariance matrix. When the cumulative distance corresponding to the confidence index is lower than a first preset threshold, the observation noise covariance is reduced; when the cumulative distance exceeds the first distance threshold, the observation noise covariance is increased. The posterior distribution of the equation parameters is obtained through extended Kalman filtering iteration.

[0115] A parameter sample set is generated by sampling from the posterior distribution of the equation parameters. For each parameter sample, the physical degradation equation is extrapolated to the failure threshold to obtain the remaining lifetime sample, and the remaining lifetime distribution is statistically obtained.

[0116] Combination Figure 2 Using a Monte Carlo sampling-based remaining lifetime distribution map as an example, this figure illustrates the remaining lifetime sample distribution obtained by extrapolating the posterior distribution of the physical degradation equation parameters 2000 times to the failure threshold. The histogram exhibits a right-skewed distribution, with peaks concentrated in the 140-160 hour range. The dashed line marking 108 hours corresponds to the 5th percentile, representing a conservative lifetime estimate under high-risk scenarios, providing probabilistic lifetime input for subsequent maintenance decisions.

[0117] Based on the type of electromechanical equipment corresponding to the reference degradation trajectory, a suitable physical degradation equation is selected. For bearing equipment, the Paris crack propagation equation or the Hertz contact fatigue equation is typically used to describe the cumulative contact damage between the rolling elements and the raceway; for gear equipment, the Arcard wear equation based on the relationship between wear volume and load cycle number is used; and for motor equipment, the insulation aging equation or the Arrhenius thermal degradation model is used. All of the above physical degradation equations are uniformly expressed in state-space form: state variables... The observed variable represents the degree of degradation at the current moment (such as crack depth, wear amount, or decrease in insulation resistance), and is the physical layer feature vector extracted from the source signal. The state transition equation describes the evolution of the degradation degree over time, while the observation equation describes the mapping relationship between the physical layer eigenvectors and the degradation degree. Together, they constitute a complete state-space model.

[0118] During the timescale adjustment phase, the piecewise time scaling factor obtained from dynamic time curvature matching is introduced into the time variable of the physical degradation equation. For the initial degradation segment, the stable degradation segment, and the accelerated degradation segment, the corresponding time scaling factors are applied respectively. (in The time variable in the physical degradation equation is scaled (numbered according to the degradation stage). Specifically, the time step in the original physical degradation equation is scaled. Replace with the equivalent time step after scaling adjustment ,in Indicates the first The equivalent time step adjusted for the degradation phase. This is the original sampling time step. This adjustment ensures that the evolution rate of the degradation equation is consistent with the actual degradation rate of the current device: when the device degradation rate is faster than the reference trajectory, The equivalent time step is amplified, accelerating equation evolution; however, when the degradation rate is slower than the reference trajectory, The equation evolution is slowed down. The time scaling factors for the three degradation stages are applied independently to ensure that the equations can accurately capture the degradation rhythm of the current device in different degradation stages.

[0119] When constructing the extended Kalman filter framework, the physical degradation equation adjusted for time scale is used as the state transition equation, and the state variables are... Imposing a degenerate monotonicity constraint, that is, requiring This is to prevent the physical inconsistency of a regression in the degree of degradation during the filtering iteration process. (Regarding the equation parameters...) (e.g., material constants in the Paris equation, wear coefficients in the Archard equation) Physical feasible region constraints are applied to limit the parameters to a reasonable range determined by materials mechanics or tribology theory. In the prediction step of the extended Kalman filter, the degenerate state and equation parameters are jointly predicted using a nonlinear state transition function. The Jacobian matrix is ​​obtained by performing a first-order Taylor expansion on the nonlinear function, and the prediction covariance matrix is ​​approximately calculated. In the update step, the physical layer eigenvectors observed at the current time are used. The predicted state is corrected, the Kalman gain is calculated, and the state estimate and covariance matrix are updated.

[0120] Adjustment mechanism of observation noise covariance matrix and confidence index of time scaling factor Directly hooked. Minimum cumulative distance from the optimal matching reference trajectory Mean cumulative distance from all candidate reference trajectories This is a joint decision reflecting the reliability of the current trajectory's match with the historical database. When Below the first preset threshold When the current trajectory closely matches the reference trajectory, indicating a high degree of confidence in the match, the noise covariance matrix will be observed. Multiply by a scaling factor less than 1 (satisfy This allows the filter to place more trust in the predictions of the physical degradation equations, reducing its reliance on observational data and thus mitigating the impact of noise interference on state estimation. Exceeding the first distance threshold ( When the current trajectory deviates significantly from any reference trajectory in the historical database, the matching confidence is low, and in this case... Multiply by a scaling factor greater than 1 (satisfy This allows the filter to rely more on real-time observation data to correct state estimates, reducing the impact of model prediction errors on the results. Through the aforementioned adaptive adjustment mechanism, the extended Kalman filter can maintain reasonable state estimation accuracy under different matched confidence conditions, and outputs equation parameters after multiple iterations. The posterior distribution of the parameter describes the uncertainty of the parameter in the form of a mean vector and a covariance matrix.

[0121] A parameter sample set is generated from the posterior distribution of the equation parameters using a Monte Carlo sampling method. ,in The total number of samples, For the first Each parameter sample. The corresponding physical degradation equation will be changed from the current degradation state. Start by extrapolating the values ​​until the degree of degradation is reached. Reaching a pre-set failure threshold During extrapolation, the time step is adjusted using the time scaling factor corresponding to each degradation stage to ensure that the extrapolated trajectory is consistent with the actual degradation rate of the current device. The difference between the extrapolated failure time and the current time is the remaining lifetime sample corresponding to that parameter sample. For all Statistical analysis was performed on a sample of remaining lifetimes to calculate the mean, standard deviation, 5th percentile, and 95th percentile, constructing a probability distribution of the remaining lifetimes. This distribution not only provides a point estimate of the remaining lifetimes but also quantifies the prediction uncertainty caused by both parameter uncertainty and model error, providing a probabilistic lifetime prediction result for subsequent maintenance decisions. When the mean of the remaining lifetime distribution is lower than the preset maintenance trigger time limit, or its 5th percentile (corresponding to a conservative estimate in high-risk scenarios) is lower than the safety margin threshold, the subsequent maintenance action optimization process is triggered, and the operation and maintenance strategy is selected based on a causal directed acyclic graph.

[0122] A causal directed acyclic graph (DAG) is constructed, containing load, wear, observed features, and remaining lifetime nodes. A pre-defined set of maintenance actions is traversed, and for each maintenance action, intervention calculations are performed on the causal DAG to obtain the remaining lifetime extension. Based on the execution cost, the maintenance action with the largest net benefit is selected as the operational strategy, including:

[0123] Calculate the partial derivatives of the state transition matrix of the physical degradation equation with respect to the load variable, wear variable, and observation variable. Use the variable pairs corresponding to the non-zero elements in the partial derivative matrix as directed edges of a causal directed acyclic graph. Use the posterior distribution of the equation parameters as the weight distribution of the directed edges to construct a causal directed acyclic graph with load as the root node, wear amount as the intermediate node, physical layer feature vector as the observation node, and remaining lifetime distribution as the terminal node.

[0124] Maintenance actions are selected sequentially from a preset set of maintenance actions, including load adjustment actions and component replacement actions.

[0125] For load adjustment actions, the incoming edges pointing to the load nodes are deleted and the load node value is fixed to the post-intervention load level. The weight distribution of the edges is used to propagate forward along the directed edges to the remaining lifetime nodes to obtain the post-intervention remaining lifetime distribution. The monotonicity constraint of the physical degradation equation and the failure threshold constraint are applied to the distribution for verification. The difference between the expected value of the post-intervention remaining lifetime distribution that passes the verification and the expected value of the remaining lifetime distribution is taken as the remaining lifetime extension.

[0126] For component replacement actions, delete the incoming edges pointing to the wear node and reset the wear node to zero wear state. Propagate along the directed edges to calculate and verify the remaining lifetime distribution after intervention, and calculate the remaining lifetime extension.

[0127] Calculate the net benefit of each maintenance action and select the maintenance action with the largest net benefit as the operation and maintenance strategy.

[0128] After obtaining the degradation state estimate and remaining lifetime distribution, it is necessary to further transform this physical information into actionable operational decisions. To this end, a causal directed acyclic graph (DAG) is constructed to encode the operating physical mechanism of electromechanical equipment in the form of a causal structure. Specifically, the partial derivative matrices of the state transition matrix of the physical degradation equation with respect to load, wear, and observation variables are calculated. Variable pairs corresponding to non-zero elements in the partial derivative matrices are identified as node pairs with causal relationships, and these node pairs are added as directed edges to the causal DAG graph structure. The weights of the directed edges are not fixed scalars, but rather use the posterior distribution of the physical degradation equation parameters as the weight distribution, thus embedding parameter uncertainty into the graph propagation process. The final constructed causal DAG graph has load as the root node, wear as the intermediate node, physical layer feature vectors as observation nodes, and remaining lifetime distribution as the terminal node, forming a complete causal link from external operating conditions to equipment health status and then to lifetime prediction.

[0129] The load node, as the root node, represents the external operating conditions the equipment experiences, and its value directly affects the evolution rate of the wear quantity node. The wear quantity node, as an intermediate node, aggregates the causal influences from the load and further transmits them to the observation node, i.e., the physical layer feature vector node. The physical layer feature vector node comprehensively reflects the current health status of the equipment and ultimately determines the value of the remaining lifetime distribution node. This hierarchical causal structure allows different types of maintenance actions to be accurately modeled on the graph in an "intervention" manner, avoiding the decision-making bias caused by the confusion between causality and correlation in traditional regression models.

[0130] Each candidate maintenance action is selected sequentially from the preset set of maintenance actions. Maintenance actions are divided into two main categories: load adjustment actions and component replacement actions. These two types of actions have different intervention positions on the causal directed acyclic graph, and their corresponding graph operations are also different, requiring separate handling.

[0131] For load adjustment actions, the "do operator" operation in Pearl intervention calculus is executed: all incoming edges pointing to load nodes in the causal directed acyclic graph are deleted, severing the observable association between load nodes and their upstream factors, causing load nodes to decouple from their original joint distribution and become fixed nodes forcibly assigned values ​​by the intervention. The value of the load node is fixed to the target load level after the intervention, such as reducing the rated load to a certain preset percentage. After completing the incoming edge deletion and node fixing, forward propagation is performed along the direction of the directed edges using the edge weight distribution, sequentially updating the distribution of wear nodes, observation nodes, and ultimately remaining lifetime nodes. Since the edge weights themselves are posterior parameter distributions, Monte Carlo integration is used to sample the weight distribution during the forward propagation process. Each sampling yields a complete causal link propagation result, which is finally summarized to obtain the remaining lifetime distribution after the intervention.

[0132] After obtaining the remaining lifetime distribution after intervention, the distribution needs to be physically constrained to ensure that it meets two basic constraints: first, the monotonicity constraint, that is, the degradation state must not show a reversal trend that violates the physical mechanism as time progresses; second, the failure threshold constraint, that is, the degradation trajectory after intervention must eventually reach the failure threshold within a certain finite time. Numerical anomalies that never fail are not permitted. For samples that fail validation, they are removed from the distribution and resampled until all retained samples satisfy the physical constraints. The expected value of the validated post-intervention remaining lifetime distribution is denoted as... Let the expected value of the original remaining life distribution before intervention be denoted as The remaining life extension caused by the load adjustment action for: .

[0133] For component replacement, the intervention location is shifted to the wear quantity node. All incoming edges pointing to the wear quantity node are deleted from the causal directed acyclic graph, detaching the wear quantity node from the causal influence of the load node. The state of the wear quantity node is then reset to zero wear, simulating the physical process of the equipment returning to its initial health level after component replacement. After resetting the wear quantity node, forward propagation is performed along the directed edge direction, sequentially updating the distribution of the observed node and the remaining lifetime node. The propagation results are then subjected to the same monotonicity constraints and failure threshold constraints as those applied to the load adjustment action. The expected value of the remaining lifetime distribution after the intervention, which passes the validation, is denoted as... The remaining life extension caused by component replacement for: .

[0134] After obtaining the remaining life extension for each maintenance action, the net benefit is calculated based on the execution cost of each action. Execution costs include direct economic costs (such as spare parts procurement costs, labor costs, and downtime losses) and indirect costs (such as the discounted loss due to reduced capacity caused by load adjustments). The remaining life extension is converted into equivalent economic benefit by multiplying the extended lifespan by the equipment operating value coefficient per unit time. to gain life extension benefits For maintenance actions Its net income The calculation method is as follows: ;in For the first The remaining life extension of each maintenance action. For the first The cost of performing a maintenance action This represents the equipment operating value coefficient per unit time. Iterate through all candidate actions in the preset maintenance action set, calculate the net benefit for each, and select the maintenance action with the largest net benefit as the current operation and maintenance strategy output.

[0135] When the maintenance action with the greatest net benefit is "no intervention" (i.e., maintaining the current operating state), the net benefit is zero. This indicates that the life extension benefit of any maintenance action is insufficient to cover its execution cost, and continuing to monitor and wait is the optimal strategy. When a certain quantile of the remaining life distribution (such as the remaining life corresponding to the 5th quantile) is lower than the preset emergency response threshold, the highest priority maintenance action is forcibly triggered, regardless of the net benefit ranking, to ensure the bottom line of safe equipment operation.

[0136] The core advantage of this causal intervention framework is that by embedding the posterior distribution of the parameters of the physical degradation equation into directed edge weights, each intervention propagation can reflect the impact of parameter uncertainty on the decision outcome. The output remaining lifetime extension also carries confidence interval information, providing maintenance personnel with a basis for risk perception decision-making rather than a single-point prediction value, thereby enabling the generation of more robust intelligent maintenance strategies under complex operating conditions.

[0137] A second aspect of the present invention provides a diagnostic prediction and intelligent monitoring system for electromechanical equipment, comprising:

[0138] The spatial matrix unit is used to acquire vibration signals from multiple spatial locations of electromechanical equipment. The time delay and attenuation coefficients between locations are extracted through the cross-correlation function to construct the spatial transfer matrix.

[0139] The source separation feature unit is used to perform sparse optimization separation based on the spatial transfer matrix to obtain bearing source signals, gear source signals and motor source signals, and extract the physical layer feature vector of each source signal;

[0140] The trajectory matching unit is used to organize the physical layer feature vectors into feature trajectories according to time, use a dynamic time warping algorithm to match the historical trajectory library to identify reference degraded trajectories, and calculate the time scaling factor of the current trajectory relative to the reference degraded trajectory.

[0141] The degradation estimation unit is used to establish the state evolution relationship based on the physical degradation equation corresponding to the reference degradation trajectory, adjust the time scale coefficient using the time scaling factor, estimate the degradation state and equation parameters using a physically constrained extended Kalman filter, and obtain the remaining lifetime distribution by sampling on the parameter distribution.

[0142] The operation and maintenance decision unit is used to construct a causal directed acyclic graph containing load, wear, observation characteristics and remaining lifetime nodes. It traverses the preset set of maintenance actions, performs intervention calculations on the causal directed acyclic graph for each maintenance action to obtain the remaining lifetime extension, and selects the maintenance action with the largest net benefit as the operation and maintenance strategy based on the execution cost.

[0143] A third aspect of the present invention provides an electronic device, comprising:

[0144] processor;

[0145] Memory used to store processor-executable instructions;

[0146] The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0147] A fourth aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0148] This invention can be a method, apparatus, system, and / or computer program product. The computer program product may include a computer-readable storage medium having computer-readable program instructions loaded thereon for performing various aspects of the invention.

[0149] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for diagnostic prediction and intelligent monitoring of electromechanical equipment, characterized in that, include: Vibration signals from multiple spatial locations of electromechanical equipment are acquired, and the time delay and attenuation coefficients between locations are extracted using a cross-correlation function to construct a spatial transfer matrix; Based on the spatial transfer matrix, sparse optimization separation is performed to obtain bearing source signals, gear source signals, and motor source signals, and physical layer feature vectors of each source signal are extracted; The physical layer feature vectors are organized into feature trajectories according to time. A dynamic time warping algorithm is used to match the historical trajectory database to identify reference degraded trajectories. The time scaling factor of the current trajectory relative to the reference degraded trajectory is calculated. The state evolution relationship is established based on the physical degradation equation corresponding to the reference degradation trajectory. The time scale coefficient is adjusted using the time scaling factor. The degradation state and equation parameters are estimated using a physically constrained extended Kalman filter. The remaining lifetime distribution is obtained by sampling the parameter distribution. Construct a causal directed acyclic graph containing load, wear, observation features, and remaining lifetime nodes. Traverse the preset set of maintenance actions, perform intervention calculations on the causal directed acyclic graph for each maintenance action to obtain the remaining lifetime extension, and select the maintenance action with the largest net benefit as the operation and maintenance strategy based on the execution cost.

2. The method according to claim 1, characterized in that, Vibration signals from multiple spatial locations of electromechanical equipment are acquired, and the time delay and attenuation coefficients between locations are extracted using a cross-correlation function to construct a spatial transfer matrix, including: Vibration sensors are arranged in the bearing housing, gearbox housing, and motor end cover of the electromechanical equipment to collect vibration acceleration signals. After time-domain alignment of the vibration acceleration signals collected at each location, the cross-correlation function of the vibration acceleration signals between any two locations is calculated. The signal propagation delay is determined by identifying the peak position of the cross-correlation function, and the signal propagation attenuation coefficient is determined by the ratio of the peak amplitude of the cross-correlation function to the normalization coefficient. The propagation delay between each position pair is constructed into a delay matrix, and the signal propagation attenuation coefficient between each position pair is constructed into an attenuation matrix. By applying an exponential decay kernel function transformation to the delay matrix and multiplying it element by element with the attenuation matrix, a spatial transmission matrix describing the propagation characteristics of the vibration signal in space is obtained.

3. The method according to claim 1, characterized in that, Based on the spatial transfer matrix, sparse optimization separation is performed to obtain bearing source signals, gear source signals, and motor source signals. Physical layer feature vectors of each source signal are extracted, including: The inverse of the spatial transfer matrix is ​​used as the initial estimate of the unmixing matrix. During the sparse optimization iteration process, orthogonality constraints and energy conservation constraints are applied to the unmixing matrix so that the separated bearing source signal, gear source signal and motor source signal satisfy the physical propagation reversibility condition. The bearing source signal is subjected to minimum entropy deconvolution to obtain an impact pulse sequence. By calculating the time interval between adjacent impact events in the impact pulse sequence, periodic patterns corresponding to the inner ring fault frequency and outer ring fault frequency of the bearing are identified, and the stability index of the periodic pattern is extracted as the bearing degradation feature. Frequency domain features and phase features are extracted from the gear source signal and the motor source signal, respectively. A time consistency verification mechanism for cross-source signals is constructed. By calculating the correlation of the changing trends of degradation features of each source signal, when the changing trend of the bearing degradation features deviates in the opposite direction from the changing trend of the frequency domain features of the gear source signal, it is determined that there is crosstalk in the source separation process and the adaptive calibration of the spatial transfer matrix is ​​triggered. The calibrated spatial transfer matrix is ​​reused for sparse optimization separation until the changing trends of degradation features of each source signal satisfy physical causal consistency, and the physical layer feature vector is obtained.

4. The method according to claim 3, characterized in that, Construct a timing consistency verification mechanism for cross-source signals, including: Calculate the rate of change of bearing degradation characteristics, gear source signal frequency domain characteristics, and motor source signal phase characteristics within a continuous time window. Construct a rate of change vector sequence from the rate of change of each source signal characteristic. Calculate the covariance matrix of the rate of change vector sequence. Extract the eigenvector corresponding to the largest eigenvalue of the covariance matrix as the main degradation direction. Calculate the projection component of the rate of change vector of each source signal feature in the main degradation direction. When the projection component of the bearing degradation feature and the projection component of the frequency domain feature of the gear source signal have opposite signs and the difference in absolute value exceeds a preset difference threshold, it is determined that there is crosstalk in the source separation. By calculating the change in the condition number of the spatial transfer matrix before and after the crosstalk occurs, the position pairs that cause crosstalk are identified. Time-domain fine-tuning is applied to the propagation delay corresponding to the position pair, and frequency-domain compensation is applied to the attenuation coefficient corresponding to the position pair to obtain the calibrated spatial transfer matrix. The inverse of the calibrated spatial transfer matrix is ​​used as the new unmixing matrix to re-perform sparse optimization separation. The temporal consistency check is repeated until the projected components of each source signal feature have the same sign and the condition number of the covariance matrix is ​​lower than the preset condition number threshold.

5. The method according to claim 1, characterized in that, The physical layer feature vectors are organized into feature trajectories according to time. A dynamic time warping algorithm is used to match historical trajectory databases to identify reference degraded trajectories. The time scaling factor of the current trajectory relative to the reference degraded trajectory is calculated, including: The physical layer feature vectors are arranged in timestamp order to form the current feature trajectory, and historical trajectories marked with degradation stage divisions are selected from the historical trajectory database as candidate reference trajectories; A dynamic time warping algorithm is used to calculate the warping path between the current feature trajectory and each candidate reference trajectory. When constructing the cumulative distance matrix, time monotonicity constraints and slope smoothness constraints are applied to the warping path. The time monotonicity constraint requires that the time index pairs of the warping path be monotonically increasing, and the slope smoothness constraint limits the slope change of adjacent time index pairs of the warping path. The candidate reference trajectory with the smallest cumulative distance that simultaneously satisfies the time monotonicity constraint and the slope smoothness constraint is selected as the reference degenerate trajectory. The reference degradation trajectory is divided into an initial degradation segment, a stable degradation segment, and an accelerated degradation segment according to the degradation stage division points marked on it. The mapping relationship of the curved path in the time dimension is extracted, and the time scaling factor of the current feature trajectory relative to the reference degradation trajectory in the corresponding interval of each degradation stage is calculated respectively. The confidence index of the time scaling factor is constructed based on the cumulative distance of the curved path. When the cumulative distance exceeds a preset distance threshold, the historical trajectory library is re-searched to update the reference degraded trajectory.

6. The method according to claim 1, characterized in that, A state evolution relationship is established based on the physical degradation equation corresponding to the reference degradation trajectory. The time scale coefficient is adjusted using the time scaling factor. The degradation state and equation parameters are estimated using a physically constrained extended Kalman filter. The remaining lifetime distribution is obtained by sampling the parameter distribution, including: Based on the type of electromechanical equipment corresponding to the reference degradation trajectory, a physical degradation equation is selected, and the physical degradation equation is expressed in state-space form, where the state variable is the degree of degradation and the observation variable is the physical layer feature vector; The time scaling factor is used as the scaling coefficient of the time variable in the physical degradation equation. The corresponding time scaling factor is applied to adjust the time scale for the initial degradation segment, the stable degradation segment, and the accelerated degradation segment, respectively. An extended Kalman filter framework is constructed, which uses the physical degradation equation as the state transition equation, imposes degradation monotonicity constraints on the state variables, and imposes physical feasible region constraints on the equation parameters; The confidence index of the time scaling factor is mapped to the adjustment coefficient of the observation noise covariance matrix. When the cumulative distance corresponding to the confidence index is lower than a first preset threshold, the observation noise covariance is reduced; when the cumulative distance exceeds the first distance threshold, the observation noise covariance is increased. The posterior distribution of the equation parameters is obtained through extended Kalman filtering iteration. A parameter sample set is generated by sampling from the posterior distribution of the equation parameters. For each parameter sample, the physical degradation equation is extrapolated to the failure threshold to obtain the remaining lifetime sample, and the remaining lifetime distribution is statistically obtained.

7. The method according to claim 6, characterized in that, A causal directed acyclic graph (DAG) is constructed, containing load, wear, observed features, and remaining lifetime nodes. A pre-defined set of maintenance actions is traversed, and for each maintenance action, intervention calculations are performed on the causal DAG to obtain the remaining lifetime extension. Based on the execution cost, the maintenance action with the largest net benefit is selected as the operational strategy, including: Calculate the partial derivatives of the state transition matrix of the physical degradation equation with respect to the load variable, wear variable, and observation variable. Use the variable pairs corresponding to the non-zero elements in the partial derivative matrix as directed edges of a causal directed acyclic graph. Use the posterior distribution of the equation parameters as the weight distribution of the directed edges to construct a causal directed acyclic graph with load as the root node, wear amount as the intermediate node, physical layer feature vector as the observation node, and remaining lifetime distribution as the terminal node. Maintenance actions are selected sequentially from a preset set of maintenance actions, including load adjustment actions and component replacement actions. For load adjustment actions, the incoming edges pointing to the load nodes are deleted and the load node value is fixed to the post-intervention load level. The weight distribution of the edges is used to propagate forward along the directed edges to the remaining lifetime nodes to obtain the post-intervention remaining lifetime distribution. The monotonicity constraint of the physical degradation equation and the failure threshold constraint are applied to the distribution for verification. The difference between the expected value of the post-intervention remaining lifetime distribution that passes the verification and the expected value of the remaining lifetime distribution is taken as the remaining lifetime extension. For component replacement actions, delete the incoming edges pointing to the wear node and reset the wear node to zero wear state. Propagate along the directed edges to calculate and verify the remaining lifetime distribution after intervention, and calculate the remaining lifetime extension. Calculate the net benefit of each maintenance action and select the maintenance action with the largest net benefit as the operation and maintenance strategy.

8. An electromechanical equipment diagnostic prediction and intelligent monitoring system, used to implement the method as described in any one of claims 1-7, characterized in that, include: The spatial matrix unit is used to acquire vibration signals from multiple spatial locations of electromechanical equipment. The time delay and attenuation coefficients between locations are extracted through the cross-correlation function to construct the spatial transfer matrix. The source separation feature unit is used to perform sparse optimization separation based on the spatial transfer matrix to obtain bearing source signals, gear source signals and motor source signals, and extract the physical layer feature vector of each source signal; The trajectory matching unit is used to organize the physical layer feature vectors into feature trajectories according to time, use a dynamic time warping algorithm to match the historical trajectory library to identify reference degraded trajectories, and calculate the time scaling factor of the current trajectory relative to the reference degraded trajectory. The degradation estimation unit is used to establish the state evolution relationship based on the physical degradation equation corresponding to the reference degradation trajectory, adjust the time scale coefficient using the time scaling factor, estimate the degradation state and equation parameters using a physically constrained extended Kalman filter, and obtain the remaining lifetime distribution by sampling on the parameter distribution. The operation and maintenance decision unit is used to construct a causal directed acyclic graph containing load, wear, observation characteristics and remaining lifetime nodes. It traverses the preset set of maintenance actions, performs intervention calculations on the causal directed acyclic graph for each maintenance action to obtain the remaining lifetime extension, and selects the maintenance action with the largest net benefit as the operation and maintenance strategy based on the execution cost.

9. An electronic device, characterized in that, include: processor; Memory used to store processor-executable instructions; The processor is configured to invoke instructions stored in the memory to execute the method according to any one of claims 1 to 7.

10. A computer-readable storage medium having computer program instructions stored thereon, characterized in that, When the computer program instructions are executed by the processor, they implement the method described in any one of claims 1 to 7.

Images