Pump station unit state comprehensive evaluation method

Abstract

The application discloses a kind of pump station unit state comprehensive evaluation method, belong to pump station unit detection technical field.This method includes: obtaining multi-source monitoring data and generating standardization monitoring data set;Time-varying mutual information between indexes is calculated, and dynamic threshold is determined in combination with current working condition parameter and water level difference correction term, to build dynamic correlation network;Physical mechanism characteristic decoupling is carried out to data, including stripping vibration signal working condition drift component based on reference curve to obtain vibration residual, and calculating equivalent standard working condition temperature based on heat balance principle;First weight based on data statistics and second weight based on network topology are calculated respectively, and comprehensive coupling weight is adaptively generated according to the consistency coefficient thereof;Finally, health state grade is determined using cloud model.The application solves the problem of difficult extraction of fault characteristics and poor model robustness under variable working conditions through deep integration of physical mechanism and data-driven, and realizes accurate evaluation of unit state.

Description

Technical Field

[0001] This invention belongs to the field of pump station unit testing technology, and in particular to a comprehensive evaluation method for the condition of pump station units. Background Technology

[0002] With the construction and development of large-scale industrial equipment, pumping station units are evolving towards larger capacity and higher speeds. These units involve complex coupling of mechanical, electrical, and hydraulic systems, and face operational challenges such as large water level fluctuations and frequent switching of operating conditions. Real-time and accurate comprehensive condition assessment of pumping station units is of significant engineering and technical importance for early identification of equipment degradation trends, prevention of catastrophic failures, and optimization of scheduling operations, and is a crucial link in ensuring the safety of water resource allocation.

[0003] Currently, the condition monitoring of pump station units mainly relies on vibration, temperature, sway, and electrical parameters collected by the SCADA (Supervisory Control and Data Acquisition) system. Existing assessment methods are generally divided into two categories: one is a discrimination method based on fixed thresholds, which sets a single alarm limit according to industry standards; the other is a data-driven intelligent diagnostic method, which uses algorithms such as neural networks and support vector machines to explore the mapping relationship between monitoring data and fault states, or analyzes the linear correlation between multiple parameters through Pearson correlation coefficient analysis. This method has achieved certain application results under steady-state conditions and has been widely deployed in automated monitoring systems.

[0004] However, existing solutions still have significant limitations in complex environments with varying operating conditions. These limitations mainly manifest as incomplete decoupling between physical mechanisms and statistical data, and insufficient adaptability when fusing heterogeneous features. Specifically, the vibration and correlation characteristics of pump station units are affected by the upstream and downstream water level difference and blade angle. Existing methods often employ static thresholds or pure data correlation analysis, neglecting the physical modulation effect of hydraulic conditions on monitoring indicators. This leads to spurious correlations or false alarms under non-design water level conditions. Temperature and vibration signals contain normal baseline components that vary with load. Existing data models often directly process the raw signals without removing the operating condition drift components, causing early, weak fault characteristics to be submerged by the large fluctuations in background noise. In terms of weight allocation, simple statistical weights (such as random forests) and mechanistic weights (such as expert experience) often conflict when dealing with outlier data. Existing technologies lack dynamic fusion mechanisms based on consistency detection, making it difficult to effectively utilize physical topology for robust correction when data patterns fail. Summary of the Invention

[0005] The purpose of this invention is to provide a comprehensive evaluation method for the condition of pump station units, in order to solve one of the aforementioned problems in the existing technology.

[0006] Technical solution: A comprehensive assessment method for the condition of pump station units, including:

[0007] Multi-source monitoring data of pump station units are acquired and time-aligned and standardized to obtain a standardized monitoring dataset.

[0008] Calculate the time-varying mutual information between various monitoring indicators in the standardized monitoring dataset, and determine the dynamic threshold based on the current operating condition parameters extracted from the standardized monitoring dataset to construct a dynamic correlation network;

[0009] Decoupling the physical mechanism features of the standardized monitoring dataset includes: removing the operating condition drift component of the vibration signal based on a preset benchmark curve to obtain the vibration residual features, and converting the temperature signal into an equivalent standard operating condition temperature based on the principle of thermal balance.

[0010] The first weight based on the statistical features of the standardized monitoring dataset and the second weight based on the topological features of the dynamic association network are calculated respectively, and the consistency coefficient between the first weight and the second weight is calculated to generate a comprehensive coupling weight.

[0011] Based on the comprehensive coupling weight, vibration residual characteristics, and equivalent standard operating temperature, the health status level of the unit is determined using a pre-built health level cloud model, and the evaluation results are obtained.

[0012] Beneficial effects: This invention solves the problems of difficult extraction of fault features and poor model robustness under varying operating conditions by deeply integrating physical mechanisms and data-driven approaches, thereby achieving accurate assessment of unit status. Attached Figure Description

[0013] Figure 1 This is a flowchart illustrating the steps of the comprehensive evaluation method for the status of pump station units in this application.

[0014] Figure 2 This is a detailed flowchart of the data preprocessing and correlation analysis steps in the embodiments of this application.

[0015] Figure 3 This is a flowchart illustrating the steps involved in generating a weighted adjacency matrix in an embodiment of this application.

[0016] Figure 4 This is a flowchart illustrating the steps of dual-source weighted adaptive coupling in the embodiments of this application.

[0017] Figure 5 This is a flowchart illustrating the steps of the evaluation process using cloud model theory in the embodiments of this application. Detailed Implementation

[0018] Example 1, such as Figure 1 As shown, a comprehensive evaluation method for the status of pump station units is described, which addresses how to integrate physical mechanisms and multi-source data under varying operating conditions to achieve an accurate evaluation of the operating status of pump station units.

[0019] Step 101: Obtain multi-source monitoring data of the pump station unit, and perform time alignment and standardization processing to obtain a standardized monitoring dataset; wherein, the multi-source monitoring data includes vibration signal, temperature signal, swing signal, electrical parameters and hydraulic parameters.

[0020] This embodiment forms the data foundation for the entire evaluation process. Specifically, multi-source monitoring data can include vibration acceleration signals, swing displacement signals, temperature signals at various locations, electrical parameter signals (such as voltage, current, and power), and hydraulic parameter signals (such as upstream and downstream water levels, blade angles, and flow rates) collected from the pump station's on-site control unit. To address the issue of inconsistent sampling frequencies among different sensors, a unified time reference (such as GPS timing or the main control system clock) is needed. Interpolation or resampling algorithms are used to align the timestamps of all signals to a unified time axis, for example, uniformly resampling to 1000Hz or higher. Furthermore, to eliminate the influence of different physical dimensions, the data needs to be standardized, for example, normalizing the vibration signals to the [0,1] interval and converting electrical parameters to per-unit values ​​to obtain a standardized monitoring dataset.

[0021] Step 102: Extract current operating condition parameters from the standardized monitoring dataset, calculate the time-varying mutual information between each monitoring indicator in the standardized monitoring dataset, determine the dynamic threshold based on the current operating condition parameters extracted from the standardized monitoring dataset, and construct a dynamic correlation network.

[0022] This embodiment aims to construct a dynamic network that reflects the nonlinear coupling relationships within the generating unit. Specifically, unlike the traditional Pearson correlation coefficient, time-varying mutual information can capture the nonlinear correlations between variables. Since changes in unit operating conditions (such as water level and power) affect the strength of the natural correlation between indicators, using a fixed threshold can lead to misjudgments of the network structure. Therefore, this embodiment introduces the concept of a dynamic threshold, which is not fixed but adjusted in real time based on current operating parameters (especially water level difference). Only when the mutual information value between two indicators exceeds this dynamic threshold is a correlation determined and a connection established in the network. In this way, the constructed dynamic correlation network can adaptively filter out false weak correlations caused by changes in operating conditions, preserving the true physical coupling relationships.

[0023] Step 103: Decouple the physical mechanism features of the standardized monitoring dataset to obtain the decoupled feature set. The decoupling includes: removing the working condition drift component of the vibration signal based on the preset working condition-vibration reference curve library to obtain the vibration residual features, and converting the temperature signal into an equivalent standard working condition temperature based on the principle of thermal balance.

[0024] This embodiment aims to eliminate the interference of operating condition changes on the monitoring of a single indicator, belonging to a mechanism-driven feature enhancement process. For vibration signals, the normal vibration level of the unit differs under different loads, and directly setting a uniform alarm value can easily lead to false alarms. This embodiment employs a subtraction strategy, subtracting the preset baseline mean value under the current operating condition from the measured signal, stripping away the normal operating condition drift component to obtain a pure vibration residual characteristic, which more sensitively reflects potential faults. For temperature signals, changes in ambient temperature and cooling water temperature can mask the equipment's own heating trend. This embodiment utilizes the principle of thermal balance to convert the current measured temperature into a temperature value under standard operating conditions (such as rated power and an ambient temperature of 20 degrees Celsius), i.e., an equivalent standard operating condition temperature, making temperature data comparable under different seasons and loads.

[0025] Step 104: Calculate the first weight based on the statistical features of the standardized monitoring dataset and the second weight based on the topological features of the dynamic association network, and calculate the consistency coefficient between the first weight and the second weight to generate the comprehensive coupling weight.

[0026] This embodiment addresses the potential conflict between data-driven and mechanism-driven approaches in feature importance assessment. The first weight is typically mined from historical data using statistical learning algorithms (such as random forests), reflecting correlation at the data distribution level. The second weight, however, is calculated based on the topology of a dynamically connected network (such as node centrality), reflecting correlation at the physical structure level. Under certain complex fault conditions, statistical regularities may fail, while physical topology is often more robust. Therefore, this embodiment calculates the consistency coefficient (such as the rank correlation coefficient) between the two. When the consistency is high, the two are weighted and fused; when the consistency is low, the fusion ratio is adaptively adjusted (e.g., increasing the proportion of the second weight) to generate a more robust comprehensive coupled weight.

[0027] Step 105: Based on the comprehensive coupling weight, vibration residual characteristics and equivalent standard operating temperature, the health status level of the unit is determined using a pre-built health level cloud model, and the evaluation results are obtained.

[0028] This embodiment describes the decision-making stage of the comprehensive assessment. The health level cloud model is an effective tool for converting qualitative concepts into quantitative values. It utilizes three numerical features—expectation, entropy, and hyperentropy—to describe fuzzy concepts such as excellent and good. By inputting the processed feature values ​​into the cloud model, the certainty of the current state belonging to each health level can be calculated. Finally, by combining comprehensive coupling weights to weightedly synthesize the membership degrees, the final health status level of the unit is determined, and an assessment result including the level conclusion and key indicator analysis is generated.

[0029] Example 2, as follows Figure 2As shown, this paper details the data preprocessing and correlation analysis steps, as well as the specific calculation details of time-varying mutual information, and solves the problem of how to accurately extract time-varying correlation features from nonlinear and non-stationary signals.

[0030] In pump station unit condition monitoring scenarios, the correlations between various monitoring indicators are not constant but dynamically change with variations in operating conditions and equipment status. Calculating time-varying mutual information using a sliding window approach can capture the dynamic evolution of these correlations. When the unit is operating normally, the indicators maintain a stable correlation pattern; however, when anomalies or fault symptoms occur, the correlations between some indicators change, manifested as changes in mutual information values. Therefore, the time-varying mutual information sequence can serve as a sensitive indicator of unit status evolution.

[0031] Step 201: Divide the standardized monitoring dataset into sliding window segments to generate a continuous time window sequence;

[0032] This embodiment forms the basis of time-varying analysis. Sliding window segmentation refers to extracting fixed-length data segments on the time axis and moving them forward according to a set step size. Specifically, the window length should cover the typical fluctuation period of the unit's state, for example, set to 300 seconds (corresponding to 30,000 sampling points, if the sampling rate is 100Hz). The sliding step size determines the temporal resolution of the analysis, for example, set to 20% of the window length, i.e., 60 seconds. In this way, the continuous monitoring data stream is segmented into a series of overlapping time window sequences containing local dynamic characteristics, allowing for subsequent window-by-window static analysis to approximate the dynamic process.

[0033] Step 202: For the monitoring index data within each time window, calculate the adaptive bandwidth based on the standard deviation and kurtosis coefficient of the sample, and use the Gaussian kernel function to estimate the probability density based on the adaptive bandwidth.

[0034] This embodiment describes a crucial step in mutual information computation—the estimation of the probability density function. To improve estimation accuracy, this embodiment abandons the fixed bandwidth method and instead employs an adaptive bandwidth strategy. Specifically, for any index sequence x within the window, its standard deviation σ is calculated. x And the number of samples n. Basic bandwidth h base The Silverman empirical formula can be used to calculate h. base =1.06*σ x *n -0.2 Among them, h base The reference bandwidth; σ x Let be the standard deviation of the sample sequence; n be the sample size. Based on this, the kurtosis coefficient k is introduced. x The bandwidth is corrected, and the corrected adaptive bandwidth h is obtained. opt=h base *(3 / k x ) 0.2 Among them, h opt The corrected adaptive bandwidth; h base The reference bandwidth; k x Let be the kurtosis coefficient of the sample sequence. This correction mechanism narrows the bandwidth when the data distribution is angular (large kurtosis) to capture details, and widens the bandwidth when the distribution is flat (small kurtosis) to smooth noise. Based on the adaptive bandwidth, a Gaussian kernel function is used to perform convolution operations on the data to obtain the marginal probability density function p(x) of the index within the window and the joint probability density function p(x,y) between the two indices.

[0035] Considering that the distribution characteristics of pump station unit monitoring data often deviate from a normal distribution, this invention introduces a bandwidth correction mechanism based on kurtosis. The kurtosis coefficient reflects the sharpness or flatness of the distribution; the kurtosis coefficient of a normal distribution is three. When the kurtosis coefficient of a sample sequence is greater than three, it indicates a peaked distribution, and the bandwidth should be appropriately reduced to preserve distribution details near the peak. When the kurtosis coefficient is less than three, it indicates a flat distribution, and the bandwidth should be appropriately increased to avoid generating false multimodal structures. The corrected bandwidth is equal to the initial bandwidth multiplied by three and divided by the square root of the actual kurtosis coefficient.

[0036] Step 203: Based on the estimated marginal probability density and joint probability density, calculate the mutual information value between each monitoring index pair to form a time-varying mutual information matrix sequence corresponding to the time window sequence; the time-varying mutual information matrix sequence is used for subsequent determination of dynamic thresholds and construction of dynamic correlation networks.

[0037] After obtaining the probability density function, numerical calculations are performed using the definition of mutual information. Specifically, mutual information I(X;Y) is defined as the expectation of the logarithm of the ratio of the product of the joint probability density and the marginal probability density. Numerically, a two-dimensional trapezoidal integral method can be used. A uniform grid is divided on the two-dimensional plane, and the integrand value at each grid point is calculated: f(x,y)=p(x,y)*log(p(x,y) / (p(x)*p(y))). Here, f(x,y) is the integrand value at grid point (x,y); p(x,y) is the joint probability density function value between the two indices; p(x) is the marginal probability density function value of the first indices; and p(y) is the marginal probability density function value of the second indices. The values ​​at all grid points are weighted and summed (integrated) to obtain the mutual information value between the two indices within that window. This process is repeated for all possible indice pairs within the window (e.g., vibration and temperature, vibration and oscillation) to construct a symmetric matrix, i.e., the local mutual information matrix. As the window slides, this series of matrices forms a time-varying mutual information matrix sequence, which fully records the evolution trajectory of the relationship between various components of the unit over time.

[0038] When calculating mutual information through numerical integration, it is necessary to handle cases where the denominator of the integrand is close to zero. When the product of marginal probability densities at a grid point approaches zero, the logarithmic term in the integrand tends towards negative infinity, leading to numerical overflow. To avoid this problem, this invention sets a lower bound truncation threshold of 10 to the power of negative 10 for the probability density product. When the actual calculated value is less than this threshold, the threshold substitutes for the actual value in the logarithmic operation. This truncation ensures the stability of the numerical calculation without affecting the accuracy of the mutual information estimation, because the contribution of the extremely low probability region to the integration result is negligible.

[0039] Example 3 describes the working condition-water level correction threshold and network feature extraction process.

[0040] Step 301: The current operating condition parameters include the unit's current calculated water level difference and operating power; determine the dynamic threshold based on the current operating condition parameters extracted from the standardized monitoring dataset; determine the current operating condition category by matching the current operating condition parameters with the preset operating condition cluster centers; and query the mutual information benchmark statistical features corresponding to the operating condition category.

[0041] An adaptive threshold model incorporating a water level difference correction term is used to calculate the dynamic threshold.

[0042] This embodiment clarifies the key physical quantities for calculating the intervention threshold. The calculated water level difference is usually obtained by subtracting the downstream water level measurement from the upstream water level measurement; it is a parameter that determines the hydraulic load and flow stability of the pumping station unit.

[0043] Step 302, the adaptive threshold model is defined as: based on the mutual information baseline statistical characteristics of the corresponding working condition category, a water level difference correction term is superimposed;

[0044] The silhouette coefficient is a comprehensive indicator for evaluating clustering quality, considering both intra-cluster compactness and inter-cluster segregation. For any sample point, the average Euclidean distance between the sample and all other sample points within the same category is calculated; this average distance is denoted as the intra-cluster distance, used to measure the degree of intra-cluster compactness. The distance between the sample and other categories is calculated; the distance to each category is defined as the average Euclidean distance between the sample and all sample points within that category. The minimum value among these inter-cluster distances is taken as the inter-cluster distance, used to measure the degree of separation from the nearest neighbor category.

[0045] The silhouette coefficient of a sample is defined as the difference between the inter-class distance and the intra-class distance, divided by the larger of the intra-class distance and the inter-class distance. The silhouette coefficient ranges from -1 to +1. A silhouette coefficient close to +1 indicates that the sample is correctly clustered and far from the boundaries of other classes; a silhouette coefficient close to zero indicates that the sample is near the class boundary and its classification is unclear; a negative silhouette coefficient indicates that the sample may have been incorrectly assigned to the current class. The overall silhouette coefficient is obtained by averaging the silhouette coefficients of all samples, and is used to evaluate the clustering quality at a specific number of clusters. The optimal number of clusters is selected by iterating through the range of candidate cluster numbers and maximizing the overall silhouette coefficient.

[0046] After determining the optimal number of clusters, the K-means clustering algorithm, specifically the K-means++ initialization strategy, is used to select initial cluster centers, improving the stability and quality of the clustering results. Traditional K-means algorithms randomly select initial cluster centers, which can easily lead to local optima and are significantly affected by initial values. The K-means++ strategy uses probability sampling to distribute the initial cluster centers as widely as possible. Specifically, a sample is randomly and uniformly selected from the sample set as the first cluster center. For each subsequent cluster center selection, the minimum distance from each unselected sample to an existing cluster center is calculated, and the square of this minimum distance is used as the probability weight for selecting that sample. A roulette wheel sampling method is then used to select the next cluster center based on this probability weight. This sampling process is repeated until all initial cluster centers are selected. This initialization strategy ensures that the initial cluster centers are evenly distributed in the sample space, effectively avoiding convergence difficulties and local optima problems caused by overly concentrated cluster centers.

[0047] Step 303: The water level difference correction term is positively correlated with the degree of deviation of the calculated water level difference from the design water level difference of the pumping station, and its correction strength is controlled by a preset water level sensitive factor.

[0048] This embodiment describes the calculation logic of the dynamic threshold. Specifically, for each time window, its operating condition category c is identified, and the mean mutual information μ under that category is queried. c and standard deviation σ c As a baseline, calculate the current dynamic threshold θ. k The specific expression for the adaptive threshold model can be: θ k =μ c +λ*σ c *(1+η*|ΔH-ΔH design | / ΔH design ), where θ k μ is the dynamic threshold for the current time window k; c σ is the benchmark for the mutual information mean under the current operating condition category c; λ is the confidence coefficient (e.g., 1.5 or 2.0); cThe mutual information standard deviation benchmark is used for the current operating condition category c; η is the water level sensitivity factor; ΔH is the current calculated water level difference; ΔH design The design water level difference of the pumping station. The physical meaning of this formula is: the further the actual water level difference ΔH deviates from the design conditions (whether too high or too low), the stronger the hydraulic excitation and nonlinear coupling within the unit will be. In this case, the threshold for determining the correlation (i.e., the threshold θ) should be appropriately increased. k (Increase the level) to avoid misinterpreting normal hydraulic fluctuations as fault associations. The water level sensitivity factor η is used to adjust the sensitivity of this correction, and is usually obtained through calibration using historical data, for example, with a value of 0.2.

[0049] Step 304, as follows Figure 3 As shown, the time-varying mutual information between each monitoring indicator in the standardized monitoring dataset is calculated to obtain the time-varying mutual information matrix. The elements in the time-varying mutual information matrix are compared with the dynamic threshold element by element. When the mutual information value is greater than or equal to the dynamic threshold, a connection edge is established between the corresponding monitoring indicator nodes, and the mutual information value is used as the edge weight to generate a weighted adjacency matrix.

[0050] This embodiment realizes the transformation from a continuous numerical matrix to a sparse network structure. Specifically, for an element M(i,j) in the mutual information matrix, if M(i,j)≧θ k If the adjacency matrix A(i,j) = M(i,j), then A(i,j) = M(i,j); otherwise, A(i,j) = 0. The constructed weighted adjacency matrix retains the strength information of strongly associated edges while eliminating noise interference from weakly associated edges.

[0051] Step 305: Based on the weighted adjacency matrix, calculate the degree centrality, betweenness centrality, and clustering coefficient of each monitoring index node as topological features of the dynamic association network.

[0052] This embodiment defines specific network characteristic metrics. Degree centrality reflects the breadth of a node's connections and is calculated as the sum of the weights of all edges connecting that node. Betweenness centrality reflects a node's bridging role in network information transmission, specifically obtained by calculating the proportion of paths passing through that node in all shortest paths using the Floyd-Warshall algorithm. Clustering coefficient reflects the tightness of interconnections between a node's neighbors. For example, for node i, its weighted clustering coefficient can be calculated as the ratio of the sum of the actual weights of edges connecting neighboring nodes to the theoretically maximum possible sum of weights.

[0053] The calculation of betweenness centrality relies on enumerating the shortest paths between any two points in a network. This invention employs the Floyd-Warshall algorithm to calculate the full-source shortest paths in a weighted network. This algorithm iteratively updates the shortest path length between any two points by progressively introducing intermediate nodes. Initially, the path length between any two points is set to the distance between the edges between them; if there is no direct edge between the two points, it is set to infinity. In the k-th iteration, it considers whether using node k as an intermediate node can shorten the path between any two points i and j. If the path length through node k is less than the currently known shortest path length, the shortest path is updated. After traversing all possible intermediate nodes, the algorithm obtains a matrix of shortest path lengths between any two points.

[0054] Step 306, constructing the dynamic association network also includes the fusion of the time dimension of topological features; specifically, the topological feature matrix calculated for each window in the sliding window sequence is weighted and averaged, and the weight allocation adopts a mechanism that decays exponentially over time, so that the topological features corresponding to the window that is closer to the current time have a greater weight contribution.

[0055] This embodiment further smooths the network features in the temporal domain. Since the topological features of a single window may fluctuate, this embodiment introduces an exponentially decaying time-weighted mechanism. Specifically, let the current window be time t, and the preceding windows be t-1, t-2, ..., tn (n is a positive integer representing the number of preceding windows to be traced back), then the fused feature F... smooth F(t) can be updated recursively using the following formula: smooth (t)=α*F raw (t)+(1-α)*F smooth (t-1). Wherein, F smooth (t) represents the smoothed topological features after fusion at time t; α is the smoothing coefficient (e.g., 0.3), F raw (t) represents the original topological features calculated in the current window; F smooth (t-1) represents the smoothed topological features at time t−1. This processing ensures that the final features can quickly respond to changes in the current state while retaining the inertial information of historical states. The smoothing coefficient α typically ranges from 0.2 to 0.5. When α is smaller, the fused features retain historical states for a longer period; when α is larger, the fused features are more sensitive to the current state. In this embodiment, a value of 0.3 is preferred to balance response speed and stability.

[0056] When fusing topological features from multiple time windows into a comprehensive topological feature, an exponentially decaying time-weighted strategy is adopted. The design principle of this strategy is that the state information of recent windows is more valuable for the current evaluation and should be given a larger fusion weight; the state information of earlier windows is relatively less valuable and should be given a smaller fusion weight.

[0057] The formula for calculating time weight is as follows: the weight of the k-th window is equal to the power function value with the natural constant e as the base, a negative decay coefficient multiplied by the difference between the total number of windows and the current window number as the exponent. A typical value for the decay coefficient is 0.1. This value causes the weight of historical data ten windows ago to decay to about 37% of the weight of the current window, which can both retain some historical information and highlight the influence of recent state.

[0058] Example 4 describes the decoupling of operating conditions and the extraction of sparse features, and addresses how to accurately extract subtle early fault features when the vibration baseline fluctuates significantly due to variable operating conditions.

[0059] Operating condition drift refers to the phenomenon of baseline deviation of monitoring indicators caused by changes in the operating conditions of the unit. During the operation of the pumping station unit, when operating parameters such as operating power, water level difference, and blade angle change, the values ​​of monitoring indicators such as vibration and temperature will change accordingly. This change is a normal physical response of the unit to changes in operating conditions, rather than a symptom of failure.

[0060] Taking vibration signals as an example, when the operating power increases, the vibration amplitude will increase accordingly due to the increased hydraulic load; when the blade angle is adjusted, the vibration spectrum characteristics will also change due to the change in flow regime. If the vibration signal is not decoupled from its operating condition and the increase in vibration amplitude is directly judged as abnormal, a large number of false alarms will occur. Operating condition drift component stripping eliminates the normal response components caused by changes in operating conditions and extracts the residual signal that truly reflects changes in the unit's health status. If there are changing components in the residual signal, it may indicate fault symptoms such as mechanical wear, imbalance, or loosening.

[0061] Step 501: Determine the operating condition category of the unit based on the current operating condition parameters, and retrieve the vibration mean curve corresponding to the operating condition category from the preset operating condition-vibration reference curve library.

[0062] This embodiment establishes a physical benchmark for signal processing. The operating condition-vibration benchmark curve library is a multidimensional lookup table or database, with its index key being the operating condition category label (e.g., the category ID obtained through K-means clustering). For each operating condition, the library stores a standard waveform with a length of one complete vibration cycle, i.e., the vibration mean curve. This curve represents the ideal vibration pattern of the unit under a healthy state at that specific head and power combination. In specific implementation, the system reads the current power P, speed n, and water level difference H, calculates their respective cluster centers, matches the closest operating condition category ID, and loads the corresponding benchmark waveform vector V from memory or the database. mean .

[0063] Step 502: Subtract the real-time vibration signal from the vibration mean curve in the standardized monitoring dataset point by point to remove the vibration baseline determined by the working conditions and obtain the residual signal containing non-stationary abnormal components.

[0064] This embodiment performs decoupling. Traditional vibration monitoring often directly sets an alarm threshold for the total vibration amplitude (e.g., peak-to-peak value not exceeding 500 μm). However, under high load or severe hydraulic conditions, normal vibration may itself be close to this alarm threshold, leading to false alarms. Conversely, under low load, the small vibration increments caused by early faults are easily overwhelmed. This embodiment calculates the difference signal: V residual (t)=V measured (t)-V mean (t). Wherein, V residual (t) represents the vibration residual signal at time t; V measured (t) represents the measured vibration signal at time t; V mean (t) represents the value of the vibration mean curve under the corresponding working condition at time t. This effectively filters out the conventional background vibrations determined by fluid dynamics and mechanical rotation. The resulting residual signal has its energy mainly concentrated on abnormal events such as impacts, friction, or component loosening, thus improving the signal-to-noise ratio.

[0065] Step 503: Perform energy normalization processing on the residual signal to obtain vibration residual characteristics;

[0066] This embodiment is used to standardize data units. Since the residual amplitude may still differ under different operating conditions, in order to facilitate subsequent standardized processing, the residual signal is usually divided by its root mean square (RMS) value to give it unit energy.

[0067] An overcomplete dictionary is a dictionary matrix whose number of atoms exceeds the signal dimension. In traditional signal representation, the number of dictionary atoms equals the signal dimension, and the signal can be uniquely represented by dictionary atoms. In overcomplete dictionary representation, because the number of atoms exceeds the signal dimension, the same signal can be represented by multiple combinations of atoms. This redundancy makes sparse representation of the signal possible. By optimizing the selection of a small number of most relevant atoms for linear combination, an overcomplete dictionary can represent the signal with fewer non-zero coefficients, achieving sparse decomposition of the signal.

[0068] In fault feature extraction applications, an overcomplete dictionary is characterized by the fact that it can contain feature atoms of various fault modes. When the measured vibration signal contains a certain type of fault component, the atoms corresponding to that type of fault will be activated (with non-zero sparsity coefficients). By analyzing which atoms are activated and their activation intensity, the fault type can be identified and the fault degree can be calculated.

[0069] K-SVD is a dictionary learning algorithm based on singular value decomposition. It achieves adaptive dictionary learning by alternately optimizing sparse coefficients and dictionary atoms. The algorithm consists of two alternating phases: sparse encoding and dictionary updating.

[0070] In the sparse coding stage, with the current dictionary matrix fixed, the sparse representation coefficients of each sample signal in the training sample library are calculated. The sparse coefficients are calculated using the orthogonal matching pursuit algorithm, which minimizes the number of non-zero coefficients while ensuring the reconstruction error is less than a given threshold. The sparsity threshold is set to 10% to 15% of the number of dictionary atoms; this range ensures both sufficient signal representation and maintains the sparsity of the coefficients.

[0071] During the dictionary update phase, each atom in the dictionary is updated sequentially. For the j-th atom, a subset of samples using this atom is selected from the training samples, i.e., samples whose j-th element in the sparse coefficient vector is non-zero. The residual of this sample after removing contributions from other atoms is calculated, and all residuals are arranged column-wise to form a residual matrix. Singular value decomposition is performed on the residual matrix, decomposing it into the product of a left singular vector matrix, a singular value diagonal matrix, and a right singular vector matrix. The left singular vector corresponding to the largest singular value is taken as the updated j-th atom, and the corresponding right singular vector is multiplied by the largest singular value to obtain the updated sparse coefficients. This update method is optimal in the least squares sense, and can minimize the reconstruction error between the sample and the dictionary representation.

[0072] One iteration is completed by sequentially performing the above update operation on all atoms in the dictionary. The sparse coding phase and dictionary update phase are repeated until the dictionary matrix converges or the maximum number of iterations is reached. The convergence criterion is that the norm of the difference between the dictionary matrices of two adjacent iterations is less than a preset threshold. The maximum number of iterations is set to fifty.

[0073] Step 601: Load the pre-built component-sensitive atom library, which contains multiple characteristic atoms associated with specific fault types;

[0074] This embodiment introduces sparse representation theory. The component-sensitive atom library is trained on a large number of historical fault samples using dictionary learning algorithms such as K-SVD. Each atom in the library is a short-time waveform fragment, and each atom has a clear label indicating its high sensitivity to a specific fault (such as bearing inner ring wear, blade cracks, or seal wear). For example, atom d1 may correspond to the impact attenuation waveform caused by bearing ball spalling, and atom d5 may correspond to the low-frequency sinusoidal waveform caused by rotor imbalance.

[0075] Step 602: The orthogonal matching pursuit algorithm is used to perform sparse decomposition on the vibration residual features and solve the sparse coefficient vector under the sensitive atom library of the component.

[0076] This embodiment utilizes the Orthogonal Matching Pursuit (OMP) algorithm to decompose complex residual signals into linear combinations of multiple atoms. The specific process is as follows: Initialize the residual r0 to equal the original signal, and set the sparse coefficient vector x to all zeros. In the k-th iteration, calculate the current residual r0. k-1 With all atoms d in the dictionary j Find the atom d with the largest absolute value of the inner product. best Add its index to the support set. Project the original signal onto the subspace spanned by the support set, and update the sparse coefficients using the least squares method. Update the residual r. k The above process is repeated until the residual energy is less than a preset threshold or the set sparsity (number of non-zero coefficients) is reached. The resulting sparse coefficient vector is a vector in which most elements are 0 and only a few positions have non-zero values. The position and magnitude of these non-zero values ​​accurately reveal the fault components contained in the signal.

[0077] Step 603: Based on the atomic tags corresponding to the non-zero elements in the sparse coefficient vector, classify and accumulate the vibration energy contribution values ​​of each component, and generate the component sensitive vibration feature vector.

[0078] This embodiment transforms mathematical sparse coefficients into physical fault diagnosis indicators. Assuming the 1st, 5th, and 10th elements in the sparse coefficient vector are non-zero, and a lookup table shows that the 1st and 10th elements belong to bearing-related fault characteristics, while the 5th element belongs to impeller-related fault characteristics, the system sums the absolute values ​​(or squares) of the 1st and 10th coefficients to obtain the vibration energy contribution value for the bearing component; and uses the absolute value of the 5th coefficient as the contribution value for the impeller component. The final generated component-sensitive vibration feature vector is in the form of C. bearing C impeller Cseal It directly calculates the current abnormality level of each component.

[0079] In some alternative implementations, a weighting factor can be introduced to further improve the robustness of the diagnosis. For critical components (such as main bearings), a weighting coefficient greater than 1 can be preset to amplify the cumulative contribution value, giving priority to warning of potential risks of critical components under the same signal strength.

[0080] The orthogonal matching pursuit algorithm is a greedy iterative algorithm that achieves sparse decomposition of the signal by progressively selecting the atoms most correlated with the residual signal. The algorithm's execution process is as follows:

[0081] During the initialization phase, the residual signal is set as the decoupled vibration signal of the working condition to be decomposed, the selected atom index set is set to an empty set, and the iteration counter is set to zero.

[0082] In each iteration, an atom matching step is performed: the inner product of the current residual signal and each atom in the dictionary is calculated, and the absolute value of the inner product reflects the correlation between the atom and the residual. The atom with the largest absolute value of the inner product is selected as the atom selected in the current iteration, and its index is added to the set of selected atom indices.

[0083] The orthogonal projection step involves extracting the dictionary submatrix corresponding to the selected atom index set, and using the least squares method to solve for the sparse coefficient vector, minimizing the mean square error between the linear combination of the selected atoms and the original signal. The least squares solution is calculated as follows: the sparse coefficient vector equals the transpose of the dictionary submatrix multiplied by the inverse of the dictionary submatrix, then multiplied by the transpose of the dictionary submatrix, and finally multiplied by the original signal vector.

[0084] Perform the residual update step: Subtract the linear combination of selected atoms from the original signal to reconstruct the signal, and obtain the new residual signal.

[0085] Termination criteria: If the energy of the residual signal (the square of the norm) is less than 5% of the energy of the original signal, or if the number of iterations reaches the sparsity limit, then the iteration is terminated; otherwise, the iteration counter is incremented, and the process returns to the atomic matching step to continue. The residual energy threshold is set to 5% of the energy of the original signal. This value ensures reconstruction accuracy while avoiding overfitting.

[0086] After the iteration terminates, the final set of selected atom indices and the corresponding sparse coefficient vector are output. The sparse coefficients are categorized and summarized according to the fault type labels of the selected atoms. For each fault type, the absolute values ​​of the sparse coefficients of all selected atoms belonging to that type are summed to obtain the vibration response contribution value of the component corresponding to that fault type.

[0087] In some optional implementations, the component-sensitive vibration feature vector can replace or supplement the vibration residual feature as input for subsequent health level cloud model evaluation. Specifically, when the residual energy of sparse decomposition accounts for less than a preset threshold (e.g., 10%), it indicates that the signal features have been fully characterized by the atomic library, and the component-sensitive vibration feature vector is preferentially used for evaluation; otherwise, the vibration residual feature is still used for evaluation to ensure the robustness of the evaluation.

[0088] Example 5 describes thermal balance compensation and thermal inertia treatment, addressing how to eliminate the interference of ambient temperature, cooling conditions, and variable load processes on temperature monitoring values, and achieving equivalent evaluation.

[0089] Step 701: Establish the heat balance equation of the pump station unit, which includes electrical loss heat source, mechanical friction heat source, hydraulic loss heat source and cooling heat dissipation. Among them, the electrical loss heat source is proportional to the square of the current, the mechanical friction heat source is proportional to the rotational speed and bearing load, and the hydraulic loss heat source is related to the efficiency loss of flow rate and head.

[0090] This embodiment constructs a physical model of the unit's temperature change. The heat balance equation, based on the law of conservation of energy, describes the relationship between the generation and dissipation of heat within the unit per unit time. Specifically, the electrical loss heat source term Q... elec It is usually proportional to the square of the current; the mechanical friction heat source term Q mech Proportional to rotational speed and bearing load; hydraulic loss heat source term Q hydro It is related to the product of flow rate and head, as well as the efficiency loss factor. Cooling and heat dissipation term Q cool This depends on the cooling water flow rate and the temperature difference between its inlet and outlet. The static heat balance equation can be expressed as: T steady =T env +R th *(Q elec +Q mech +Q hydro -Q cool ), where T steady T is the theoretical steady-state temperature. env R represents the ambient temperature. th Q is the equivalent thermal resistance of the unit; elec For electrical loss heat source items; Q mech For mechanical friction heat source; Q hydro For the hydraulic loss heat source term; Q cool This is for cooling and heat dissipation.

[0091] Step 702: Substitute the current operating parameters in the standardized monitoring dataset into the heat balance equation to calculate the theoretical steady-state temperature under the current operating conditions; at the same time, substitute the preset rated operating parameters into the heat balance equation to calculate the reference steady-state temperature under the standard operating conditions.

[0092] This embodiment calculates two theoretical values ​​for comparison. Current parameters such as power P, speed n, and flow rate Q are collected and substituted into the above equation to calculate the equilibrium temperature that the unit should reach under the current operating conditions, denoted as T. steady_curr Set a set of standard operating parameters (e.g., rated power P). rate Rated speed n rate Design flow Q design Substituting the standard ambient temperature (20 degrees Celsius) into the same equation, the reference temperature under standard operating conditions is calculated and denoted as T. steady_ref .

[0093] The three normalized topology indicators are weighted and fused to form a comprehensive network topology weight. The fusion weight of the three topology indicators is determined based on their physical meaning and contribution to the evaluation.

[0094] The fusion weight for degree centrality is set to 0.4. Degree centrality reflects the breadth of direct correlation between monitoring indicators and other indicators. Indicators with high degree centrality are at the core of the correlation network, and their state changes will directly affect the values ​​of multiple related indicators. Therefore, they should be given a higher weight in the comprehensive evaluation.

[0095] The fusion weight for betweenness centrality is set to 0.4. Betweenness centrality reflects the pivotal role of monitoring indicators in the information transmission path. Indicators with high betweenness centrality are bridge nodes connecting different subsystems, and their anomalies often indicate the propagation of faults between different subsystems, thus having important early warning value and should be given a higher weight.

[0096] The clustering coefficient is calculated using a weighted average after reverse processing, with a weighting factor set to 0.2. The clustering coefficient reflects the degree of interconnectivity between neighboring nodes of a monitoring indicator. Indicators with low clustering coefficients indicate a lack of direct connections between their neighboring nodes; these indicators carry unique information that cannot be replaced by other indicators. Therefore, in the fusion calculation, the value minus the normalized clustering coefficient is used as the contribution term for that indicator, giving higher weights to indicators with low clustering coefficients.

[0097] The formula for calculating the comprehensive weight of the network topology is: the comprehensive weight of the i-th index is equal to 0.4 multiplied by the normalized degree centrality of the index, plus 0.4 multiplied by the normalized betweenness centrality of the index, plus 0.2 multiplied by 1 minus the normalized clustering coefficient of the index.

[0098] Step 1201 calculates the theoretical steady-state temperature under the current operating conditions, taking into account the thermal inertia hysteresis effect of each component of the unit; specifically, it uses the thermal time constants of each temperature measuring point obtained in advance to construct a first-order inertial dynamic response model, and combines the instantaneous operating parameters at the current moment with the temperature state at the previous moment to recursively calculate the theoretical steady-state temperature containing dynamic response characteristics.

[0099] This embodiment is a dynamic correction to the static calculation. Because large-scale units have a large heat capacity, when the operating conditions undergo a step change, the actual temperature will not immediately abruptly change to a new steady-state value, but will slowly approach it exponentially. To simulate this process, this embodiment introduces a thermal time constant τ. The discretized recursive calculation formula is as follows: T dynamic (k)=(τ / (τ+Δt))*T dynamic (k-1)+(Δt / (τ+Δt))*T steady_curr (k). Wherein, T dynamic (k) represents the theoretical temperature at the current moment after considering thermal inertia; τ is the thermal time constant; Δt is the sampling interval; T dynamic (k-1) is the theoretical temperature at the previous time k−1; T steady_curr (k) is the instantaneous theoretical steady-state temperature calculated based on the operating parameters at the current time k. This dynamic correction can eliminate transient false alarms during operating condition switching (e.g., when power suddenly increases, the model predicts that the temperature should also increase slowly, rather than jump immediately).

[0100] The thermal time constant is a characteristic parameter characterizing the dynamic response speed of a thermal system. It is defined as the time required for the system to reach 63.2% of the steady-state change in response to a step thermal input. This definition originates from the step response characteristics of a first-order inertial system: when the input undergoes a step change, the output gradually approaches the new steady-state value according to an exponential law, and after the time constant, the output change reaches approximately 63.2% of the steady-state change.

[0101] The method for identifying the thermal time constant τ is as follows: Select periods from historical operating data where the operating conditions undergo a step change (e.g., a power surge exceeding 20% ​​of the rated power), and extract the response curves of each temperature measuring point within those periods. For each measuring point, fit its temperature response curve to the step response model of a first-order inertial element, and use the nonlinear least squares method to solve for the time constant τ that minimizes the fitting residual. Typically, for large pump station units, the thermal time constant of the upper guide bearing is approximately 15 to 30 minutes, and the thermal time constant of the stator winding is approximately 20 to 40 minutes.

[0102] Due to differences in mass and specific heat capacity, the various components of the generator set have different thermal time constants. Smaller components, such as rolling bearings, have shorter thermal time constants, typically five to ten minutes; larger components, such as the stator core, have longer thermal time constants, typically twenty to forty minutes. These differences in thermal time constants result in varying times for the temperature at different measuring points to reach a new steady state after changes in operating conditions. During the transition phase between operating conditions, directly comparing the temperature with the expected steady-state temperature without considering thermal inertia can lead to misjudgments.

[0103] The continuous form of the first-order inertial dynamic response model is: the time constant multiplied by the derivative of the actual temperature with respect to time, plus the actual temperature, equals the steady-state temperature. The backward difference method is used to discretize this continuous differential equation, yielding a recursive calculation formula. Let the sampling time interval be Δt. The predicted steady-state temperature at the k-th sampling moment equals the actual measured temperature at the previous moment plus the sampling interval divided by the time constant, multiplied by the instantaneous steady-state temperature at the current moment, and then divided by the sum of the two sampling intervals divided by the time constant.

[0104] The physical meaning of this discretization formula is as follows: the expected temperature at the current moment is a weighted average of the historical temperature and the new steady-state temperature driven by the change in operating conditions. The weights are determined by the ratio of the sampling interval to the time constant. When the sampling interval is much smaller than the time constant, the historical temperature has a larger weight, and the temperature change exhibits a gradual characteristic. When the sampling interval is close to the time constant, the weight of the new steady-state temperature driven by the current operating conditions increases, and the temperature response accelerates.

[0105] Step 703: The difference between the theoretical steady-state temperature and the reference steady-state temperature is used to compensate for the real-time temperature monitoring value, and the equivalent standard operating temperature that eliminates the influence of operating condition differences is calculated.

[0106] This embodiment outputs the final standardized temperature index. The calculation formula is: T equiv =T measured -(T dynamic -T steady_ref Or it can be understood as: T equiv =T measured -T dynamic +T steady_ref Among them, T equiv The equivalent standard operating temperature; T measured This is the measured temperature; T dynamic The theoretical temperature under current operating conditions (including thermal inertia correction); T steady_ref This is the reference steady-state temperature under standard operating conditions. The physical meaning of this formula is: to calculate the measured temperature T. measured Compared with the current theoretical temperature T dynamic The deviation (reflecting abnormal temperature rise) is added to the standard reference temperature T. steady_ref After processing, regardless of whether the unit is currently operating at low load in winter or high load in summer, its converted equivalent standard operating condition temperature is unified to the same reference plane. Maintenance personnel only need to focus on T. equiv Whether the alarm threshold under rated operating conditions has been exceeded, there is no need to set a complex variable threshold table for different operating conditions.

[0107] Example 6, as Figure 4As shown, this paper describes the adaptive coupling of dual-source weights to address the inconsistency in feature importance determination between data-driven algorithms (random forests) and mechanism-driven methods (network topology).

[0108] Step 801: Construct a random forest classification model using the standardized monitoring dataset and corresponding historical health status labels, and determine the first weight of each monitoring indicator by calculating the decrease in accuracy of feature permutation of out-of-bag data.

[0109] This embodiment obtains statistical weights. The random forest model consists of multiple decision trees. During construction, approximately one-third of the samples are not used in the training of any particular tree; this is called out-of-bag (OOB) data. To evaluate a feature X... j The importance of this approach is highlighted in this embodiment, which employs a feature permutation strategy: while keeping other features unchanged, feature X in the OOB data is randomly shuffled. j The numerical order of features X disrupts their correspondence with the labels. j This is important; shuffling the model will decrease its prediction accuracy. The decrease in accuracy is used as the importance score for that feature, and after normalization, the first weight W is obtained. RF .

[0110] The Spearman rank correlation coefficient is used to measure the consistency between the two weight vectors. The Spearman rank correlation coefficient is a rank-based correlation measure that evaluates the consistency of the two sets of variables in their order, rather than their numerical correlation.

[0111] The calculation process is as follows: Sort each element in the random forest feature importance vector from largest to smallest value, and assign a rank to each element. The element with the largest value has a rank of one, the second largest has a rank of two, and so on. The same method is used to assign ranks to the network topology weight vector. Calculate the rank difference between the two sets of ranks for each indicator, which is the difference between the rank of the indicator in the random forest importance ranking and its rank in the network topology weight ranking.

[0112] The Spearman rank correlation coefficient is calculated as follows: the coefficient equals one minus six multiplied by the sum of the squares of the rank differences of all indicators, then divided by the total number of indicators multiplied by the square of the total number of indicators minus one. The coefficient ranges from negative one to positive one. A value close to positive one indicates that the two sets of weights are highly consistent in their ranking, a value close to zero indicates that the rankings are uncorrelated, and a value close to negative one indicates that the rankings are reversed.

[0113] Step 802: Extract the degree centrality index, betweenness centrality index, and clustering coefficient index of each monitoring indicator node in the dynamic association network, and perform weighted fusion of the above three indicators according to the preset physical meaning importance weight coefficient to determine the second weight of each monitoring indicator; wherein, the degree centrality index reflects the direct association breadth of the node, the betweenness centrality index reflects the fault propagation hub role of the node, and the clustering coefficient index reflects the local interconnection density of the node.

[0114] This embodiment obtains weights at the topology level. Unlike statistical correlation, network topology reflects the structural vulnerability within the system. Degree centrality measures the direct impact range; betweenness centrality measures the pivotal role of fault propagation; and clustering coefficient measures local interconnect density. Based on physical experience, this embodiment presets fusion coefficients (e.g., degree weight 0.4, betweenness weight 0.4, clustering weight 0.2) and synthesizes the three normalized indicators into a comprehensive topological importance, i.e., the second weight W. NET .

[0115] The fusion ratio is adaptively determined based on the consistency coefficient of the two-source weights. The design of the fusion coefficient follows these principles: when the two weight rankings are highly consistent, it indicates that the judgments driven by statistics and by mechanisms are converging, and the two weight sources should contribute equally; when there is a discrepancy between the two weight rankings, the contribution ratio of the network topology weight should be appropriately increased, because the network topology weight contains the physical correlation mechanism between monitoring indicators and can provide robust prior knowledge constraints when data patterns are not obvious or when new operating conditions are encountered.

[0116] The adaptive fusion coefficient is calculated as follows: the fusion coefficient α of the network topology weights is equal to 0.5 plus 0.3 multiplied by 1 minus the Spearman rank correlation coefficient. When the correlation coefficient is equal to 1, α is equal to 0.5, with each weight contributing 50%; when the correlation coefficient is equal to 0, α is equal to 0.8, with the network topology weights contributing 80%; and when the correlation coefficient is equal to -1, the calculated value is 1.1.

[0117] To avoid the omission of certain weight sources in extreme cases, a boundary constraint is set for the fusion coefficient: the value of α is limited to between 0.3 and 0.9. If the α value calculated by the formula is less than 0.3, then 0.3 is used; if it is greater than 0.9, then 0.9 is used. This boundary constraint ensures that even if the judgments of the two methods are highly consistent or opposite, a moderate contribution of the other method is preserved, enhancing the robustness of weight fusion.

[0118] Step 901: Sort the first weight and the second weight respectively, obtain the rank vector, and calculate the Spearman rank correlation coefficient between the two as the consistency coefficient.

[0119] This embodiment calculates the degree of consistency between the two weight sources. Because W RF and W NETThe numerical distributions of these factors may differ, making a direct comparison of their magnitudes meaningless. Therefore, this embodiment compares their ranking. For example, if both vectors rank first for the vibration index, then they have high consistency. Specifically, the weight vector is converted into a rank vector R. RF and R NET Calculate the sum of squares of the differences in rank between the two, D=Σ(R RF_i -R NET_i ) 2 The Spearman rank correlation coefficient ρ is obtained as 1-6D / (N*(N)). 2 -1), where D is the sum of squares of the rank differences; R RF_i R is the rank of the i-th index in the first weight vector; NET_i Let be the rank of the i-th indicator in the second weight vector; ρ is the Spearman rank correlation coefficient (i.e., the consistency coefficient); N is the total number of indicators. The consistency coefficient ρ ranges from [-1, 1].

[0120] Step 902: Adaptively calculate the network topology fusion coefficient based on the consistency coefficient. The lower the consistency coefficient, the larger the value of the network topology fusion coefficient. When there is a divergence in the weight distribution, the contribution ratio of the weight based on the physical mechanism is increased.

[0121] When ρ approaches 1, it indicates a high degree of consistency between the statistical regularity and the physical structure, allowing for balanced acceptance (e.g., each accounting for 50%). When ρ approaches 0 or even becomes negative, it suggests abnormal disturbances in the statistical data or the emergence of new fault modes not seen in the training set. In this case, data-driven results are unreliable, while the physical topology (such as the connection between the impeller and the motor in the main bearing) represents relatively constant and reliable prior knowledge. Therefore, this embodiment designs a negative correlation mapping function, for example, α = 0.5 + 0.5 * (1 - ρ) / 2, where α is the network topology fusion coefficient and ρ is the Spearman rank correlation coefficient. Alternatively, a simpler linear piecewise function can be used to ensure that the network topology fusion coefficient α monotonically increases as ρ decreases.

[0122] Step 903: Based on the network topology fusion coefficient, the first weight and the second weight are weighted and summed to obtain the comprehensive coupling weight.

[0123] This embodiment outputs the final weight. The calculation formula is: W final =α*W NET +(1-α)*W RF Among them, W final The comprehensive coupling weight is α, which is the network topology fusion coefficient; W is the network topology fusion coefficient. NET The second weight (network topology weight); W RFThe first weight (Random Forest weight) is used. Through this dynamic weighting mechanism, the system can leverage the statistical advantages of big data under normal operating conditions (RF weight dominance), while automatically retreating to a robust physical mechanism defense line (topology weight dominance) when sudden unknown failures or data quality degradation occur, ensuring the robustness of the evaluation results.

[0124] In some alternative implementations, a safety lower limit can also be set, for example, specifying that (1-α) must not be lower than 0.2, to prevent the statistical characteristics of real-time data from being ignored in extreme cases.

[0125] The cloud model is an artificial intelligence approach that addresses the uncertainty of the transformation between qualitative concepts and quantitative values. Traditional fuzzy set theory uses deterministic membership functions to describe fuzzy concepts, with membership values ​​uniquely determined by the element's position. However, in actual cognition, qualitative judgments of the same thing inherently exhibit random fluctuations; elements at the same position may acquire different membership degrees at different evaluation times. The cloud model adds randomness to the fuzziness by introducing a hyperentropy parameter, making the membership calculation exhibit random fluctuations, which better aligns with the uncertain nature of human cognition.

[0126] The cloud model uses three numerical features to describe the computational representation of qualitative concepts. The expected value Ex represents the central position of the cloud droplet distribution in the universe of discourse, signifying the most typical quantitative manifestation of the qualitative concept—that is, the value that best represents the concept. Entropy En measures the fuzziness of the qualitative concept, reflecting the coverage area of ​​the cloud droplet in the universe of discourse; a larger entropy value indicates a more fuzzy boundary and a wider coverage area. Hyperentropy He is a measure of the uncertainty of entropy, reflecting the degree of random fluctuation in the cloud droplet boundary; a larger hyperentropy indicates more significant random fluctuations in membership degrees each time a cloud droplet is generated.

[0127] In unit health status assessment, the boundaries of health levels are inherently ambiguous and random. For example, the boundary between a good state and a satisfactory state is not a clear dividing line, but rather a gradual transition zone. The classification of the same comprehensive score may fluctuate at different assessment times and under different operating conditions. The cloud model naturally characterizes this dual uncertainty through three parameters, making the assessment results more robust and interpretable.

[0128] Example 7 describes the comprehensive assessment and source analysis of health status, addressing how to handle the randomness and ambiguity in monitoring data, achieving quantitative grading and qualitative interpretation of unit status, and quickly locating key factors leading to status deterioration.

[0129] Step 1001: Determine the health status level of the unit using a pre-built health level cloud model, including:

[0130] Multi-source monitoring data of pump station units are acquired and time-aligned and standardized to obtain a standardized monitoring dataset.

[0131] Calculate the time-varying mutual information between various monitoring indicators in the standardized monitoring dataset, and determine the dynamic threshold based on the current operating condition parameters extracted from the standardized monitoring dataset to construct a dynamic correlation network;

[0132] Decoupling the physical mechanism features of the standardized monitoring dataset includes: removing the operating condition drift component of the vibration signal based on a preset benchmark curve to obtain the vibration residual features, and converting the temperature signal into an equivalent standard operating condition temperature based on the principle of thermal balance.

[0133] The first weight based on the statistical features of the standardized monitoring dataset and the second weight based on the topological features of the dynamic association network are calculated respectively, and the consistency coefficient between the first weight and the second weight is calculated to generate a comprehensive coupling weight.

[0134] Based on the comprehensive coupling weight, vibration residual characteristics, and equivalent standard operating temperature, the health status level of the unit is determined using a pre-built health level cloud model, and the evaluation results are obtained.

[0135] Step 1002, as follows Figure 5 As shown, the vibration residual characteristics, equivalent standard operating condition temperature, and swing characteristics, electrical characteristics, and hydraulic characteristics in the standardized monitoring dataset are input into the preset health level cloud model. The positive cloud generator is used to calculate the degree of membership of each indicator relative to different health levels.

[0136] This embodiment employs cloud model theory to address uncertainties in the evaluation process. The health level cloud model consists of multiple sub-clouds, each corresponding to a health level (e.g., Excellent, Good, Acceptable, Attention, Deterioration), defined by three numerical features: expected value Ex, entropy En, and hyperentropy He. Ex represents the central value of the level; En represents the acceptable range of the level; and He represents the ambiguity of the level's boundaries. The forward cloud generator is an algorithm that maps quantitative values ​​to qualitative concepts. Specifically, for a given input index value x... i Given a specific cloud level (Ex, En, He), the algorithm generates normally distributed random numbers En' with En as the expectation and He as the standard deviation; calculate x. i The membership degree of this level is μ = exp(-(x) i -Ex) 2 / (2*En' 2 Where μ is the degree of certainty of membership of a single index; x iThe input index value is denoted as ; Ex is the expected value of the health level cloud model; En′ is a normally distributed random number generated with entropy En as the expected value and hyperentropy He as the standard deviation. To eliminate errors caused by randomness, the above process is usually repeated multiple times (e.g., 100 or 1000 times), and the average of all calculation results is taken as the final single-index membership determination. This method reflects the random fluctuation characteristics of actual monitoring data under critical conditions better than traditional hard threshold determination or fuzzy membership functions.

[0137] When calculating the membership certainty of a single indicator to each health level, the forward cloud generator algorithm is used. Because the hyperentropy parameter of the cloud model introduces randomness into the membership calculation, the results of a single calculation are volatile. To obtain stable membership certainty estimates, the Monte Carlo method is used to perform multiple repeated calculations and then average them.

[0138] The specific process is as follows: To calculate the membership certainty of an indicator's contribution value relative to a certain health level, a positive cloud generator process is executed one hundred times. In each execution, a random temporary entropy value is generated, which follows a normal distribution with mean En and standard deviation He. Based on this temporary entropy value, the membership certainty is calculated. The certainty is equal to the power function value of the difference between the indicator's contribution value and its expected value, with a base of the natural constant e and a negative exponent, divided by the square of twice the temporary entropy value. The arithmetic mean of the membership certainty values ​​obtained from the one hundred calculations is taken as the stable membership certainty estimate of the indicator for that health level.

[0139] Setting the number of repetitions to one hundred strikes a balance between computational efficiency and estimation accuracy. Too few repetitions will lead to large fluctuations in the estimated value, while too many repetitions will increase the computational burden but offer limited improvement in accuracy.

[0140] Step 1003: Based on the comprehensive coupling weight, the membership certainty of the single index is weighted and summed to generate a comprehensive membership certainty vector that reflects the distribution probability of the overall status of the unit at each health level.

[0141] This embodiment performs multi-indicator information fusion. Assume the unit has M monitoring indicators and 5 health levels. After the previous calculation, an M-row, 5-column membership matrix is ​​obtained. Now, the comprehensive coupling weight vector W is used... final (Of length M), perform a weighted summation on each column of this matrix. The formula is: μ combined (j)=Σ i=1 M (W final (i)*μ(i,j)). Where μ combined (j) represents the overall probability that the current state of the unit belongs to the j-th health level; M is the total number of monitored indicators; W final(i) represents the comprehensive coupling weight of the i-th indicator; μ(i,j) represents the single-indicator membership certainty of the i-th indicator to the j-th health level. These five calculated values ​​constitute the comprehensive membership certainty vector. After normalization, this vector can intuitively display the distribution of the unit's status at each health level. For example, a vector of [0.1,0.7,0.2,0,0] indicates that the unit has a 70% probability of being in a good condition, but also a 20% probability of having slipped into an acceptable condition.

[0142] Step 1004: Determine the current health status level based on the maximum value in the comprehensive membership determination vector, calculate the distribution entropy and concentration of the comprehensive membership determination vector, and evaluate the reliability of the determination result.

[0143] This embodiment outputs the final conclusion and performs self-verification. Based on the principle of maximum membership, the level corresponding to the element with the largest value in the vector is selected as the final health status level. Simultaneously, to prevent misjudgment, this embodiment introduces a credibility assessment mechanism. Concentration can be measured by calculating the ratio of the difference between the maximum and second-largest membership; a larger difference indicates a more definitive judgment. Distribution entropy can be calculated using the information entropy formula: H = -Σp j *log(p j ), where H is the membership degree distribution entropy; p j This represents the normalized membership degree. If the distribution entropy is too high or the concentration is too low, it indicates that the unit is in the transition zone (i.e., the fuzzy zone) between two health levels. In this case, the system will automatically mark the assessment result as low confidence, prompting maintenance personnel to manually intervene and verify it, rather than blindly outputting a conclusion.

[0144] The reliability of the membership degree calculation results is evaluated, and the results are graded based on two indicators: the degree of concentration and the degree of dispersion of the membership degree distribution.

[0145] The membership concentration index is calculated as follows: extract the maximum and second largest values ​​from the normalized comprehensive membership determination vector. The concentration is equal to the difference between the maximum and second largest values ​​divided by the maximum value. A concentration close to one indicates that the maximum membership degree far exceeds other levels, and the state determination is highly certain; a concentration close to zero indicates that the membership degrees of the two largest levels are close, and the state attribution is ambiguous.

[0146] The membership degree distribution entropy is calculated as follows: For each element in the normalized comprehensive membership degree vector, calculate its product with its own natural logarithm. Sum all the products and take the negative value to get the distribution entropy. The smaller the distribution entropy, the more concentrated the distribution is in a few levels, and the lower the uncertainty; the larger the distribution entropy, the more dispersed the distribution is in multiple levels, and the higher the uncertainty.

[0147] The credibility level is graded based on the two indicators mentioned above. High credibility is defined as follows: concentration greater than 0.5 and distribution entropy less than 0.8. In this case, the state determination is clear, and the assessment result can be directly accepted. Medium credibility is defined as: concentration between 0.3 and 0.5, or distribution entropy between 0.8 and 1.2. In this case, the state determination has some uncertainty, and it is recommended to combine it with other information for comprehensive judgment. Low credibility is defined as: concentration less than 0.3, or distribution entropy greater than 1.2. In this case, the state assignment is ambiguous, and manual review or reassessment after supplementing data is required.

[0148] Step 1401, after determining the health status level of the unit, also includes the identification of key influencing indicators:

[0149] The overall health score is mapped to a standard scoring range of 0 to 100 points for easier intuitive understanding and horizontal comparison. The mapping formula is: the percentage health score equals 100 multiplied by the difference between the overall health score and the expected value of the deterioration level, then divided by the difference between the expected value of the excellent level and the expected value of the deterioration level.

[0150] The mapping formula is designed based on the following principle: the expected value of an excellent level corresponds to 100 points, the expected value of a poor level corresponds to 0 points, and intermediate scores are calculated using linear interpolation. For example, if the overall health score is exactly equal to the expected value of an excellent level, the percentage score is 100 points; if it is equal to the expected value of a poor level, the percentage score is 0 points; if it is in between, the corresponding percentage score is calculated proportionally.

[0151] The percentage-based scoring system can intuitively display the relative position of the unit's current health status within the excellent to deteriorated range, making it easy for maintenance personnel to quickly judge the degree of condition, and also facilitating horizontal comparison and trend tracking between different units at different times.

[0152] Step 1402: Calculate the degree of deviation of the current value of each monitoring indicator from the expected value of the excellent state, and calculate the contribution of each monitoring indicator to the decline in health score by combining the comprehensive coupling weight.

[0153] This embodiment aims to calculate the responsibility of each indicator for a decline in overall health. Simply knowing that the unit is in a state of alert is insufficient; maintenance personnel are more concerned with where the problem lies. Specifically, it calculates the current value x of each indicator. i The expected value Ex of the excellent cloud model excellent The standardized distance between them, i.e., the degree of deviation D i =|x i -Ex excellent | / Range i Multiply the degree of deviation by the weight W of the indicator. final (i) yields the contribution of this indicator, Impact. i =W final(i)*D i Among them, D i x represents the degree of deviation of the i-th indicator; i Let Ex be the current value of the i-th indicator; excellent The expected value of the excellent-level cloud model; Range i The normal fluctuation range or normalized amplitude of the i-th indicator; Impact i W represents the contribution of the i-th indicator. final (i) represents the comprehensive coupling weight of the i-th indicator. This formula reflects that important indicators with large deviations are the main cause of the decline in health scores. If an indicator has a large deviation but a very low weight, or a high weight but no deviation, its contribution is relatively small.

[0154] Step 1403: Sort the monitoring indicators according to their contribution, extract the indicators with the highest contribution as the key influencing factors leading to the current health status and output them.

[0155] This embodiment outputs diagnostic clues. The system's contribution to all indicators is shown in the figure. i Sort the results in descending order and extract the top N (e.g., top 3) indicators along with their corresponding physical values ​​and deviation directions. For example, the output might show: Key Influencing Factor 1: Water-guided bearing vibration (contribution 0.45, deviation +30%); Key Influencing Factor 2: Stator winding temperature (contribution 0.30, deviation +15%). This provides direct guidance for subsequent inspection and maintenance.

[0156] Example 8 describes the offline preparation work before putting this example into online operation, and solves how to extract prior knowledge from historical data and build various model libraries and parameter sets required to support online evaluation.

[0157] Step 1101, the pre-built models or parameter libraries involved in the method are pre-constructed in the following ways:

[0158] The working condition categories and cluster centers are obtained by training based on the historical working condition feature vector set, using the K-means clustering algorithm and determining the optimal number of clusters with the silhouette coefficient as the optimization index.

[0159] The Working Condition-Vibration Reference Curve Library is obtained by selecting vibration signal segments under historical normal conditions for each type of working condition, aligning them over time, and calculating the mean curve.

[0160] The parameters of the heat balance model are obtained by identifying the coefficients of the heat balance equation based on historical steady-state operating data, with the goal of minimizing temperature reconstruction error;

[0161] The health level cloud model is obtained by using historical sample data with calibrated health levels and extracting the expected value, entropy, and hyperentropy parameters of each level using the reverse cloud generator algorithm.

[0162] Step 1102, the working condition category and cluster center, are obtained by training based on the historical working condition feature vector set, using the K-means clustering algorithm and determining the optimal number of clusters with the silhouette coefficient as the selection index;

[0163] This embodiment describes the construction process of the operating condition knowledge base. Operating records covering the entire operating condition range are extracted from the historical database, and a feature vector set composed of parameters such as upstream water level, downstream water level, power, blade angle, and rotational speed is constructed. To avoid local optima caused by random initialization, the K-means++ algorithm is preferably used to initialize the cluster centers. Regarding the selection of the number of clusters K, this embodiment uses the silhouette coefficient method: possible K values ​​are traversed within a preset range (e.g., 2 to 10), and the average silhouette coefficient of all samples is calculated for each case. The silhouette coefficient combines intra-cluster compactness and inter-cluster separation; a value closer to 1 indicates better clustering results. The K value corresponding to the maximum silhouette coefficient is selected as the optimal number of clusters, and the K centroid vectors obtained from this clustering are saved as a reference benchmark for online operating condition identification.

[0164] The reverse cloud generator is an algorithm that extracts the three parameters of a cloud model from sample data. For a historical sample set of a certain health level, assuming the set contains N samples, the comprehensive score of each sample constitutes a numerical sequence.

[0165] The expected value Ex is calculated using the sample mean formula: Ex equals the sum of all sample scores divided by the sample size N. The expected value represents the central position of the overall health score for that health level.

[0166] The entropy En is calculated using an estimation method based on the mean of the sample absolute deviation. The mean of the sample absolute deviation (MAD) is calculated as the sum of the absolute values ​​of the differences between the overall sample score and the expected value Ex, divided by the sample size N. Based on the theoretical relationship between the mean of the absolute deviation and the standard deviation under a normal distribution, the estimated entropy is equal to 1.2533 multiplied by MAD. This estimation method is more robust than directly using the sample standard deviation and is insensitive to outliers.

[0167] The calculation of hyperentropy He is based on the relationship between variance, entropy, and hyperentropy in cloud model theory. The square of the sample variance S is calculated as the sum of the squares of the differences between the overall sample score and the expected value Ex, divided by the sample size minus one. According to cloud model theory, the sample variance equals the square of the entropy plus the square of the hyperentropy. Therefore, the hyperentropy He is equal to the square root of the difference between the sample variance and the square of the entropy. If the value within this square root is negative, it indicates that the hyperentropy is extremely small; in this case, the hyperentropy is taken as one percent of the entropy as the minimum value.

[0168] Step 1103, the working condition-vibration reference curve library, is obtained by selecting vibration signal segments under historical normal conditions for each type of working condition, aligning them over time, and calculating the mean curve.

[0169] This embodiment describes a method for extracting physical benchmarks. After completing the operating condition clustering, the periods marked as healthy in the historical data are archived according to the operating condition category. For all vibration signal segments belonging to the same operating condition category c, phase alignment is performed (e.g., alignment using key phase signals or maximum cross-correlation values) to ensure that the peaks and troughs of the waveforms correspond. The average value of all aligned segments is calculated point-by-point on the time axis to obtain the vibration mean curve V under that operating condition. mean_c In addition, the standard deviation curve V can be calculated. std_c This is used to construct the envelope range. The baseline curve, obtained based on a large number of historical samples, reflects the actual inherent characteristics of the unit under specific operating conditions better than a single theoretical waveform.

[0170] Step 1104, the thermal balance model parameters, are obtained by identifying the coefficients of the thermal balance equation based on historical steady-state operating data and with the goal of minimizing temperature reconstruction error, using the regularized least squares method.

[0171] This embodiment describes the process of identifying thermodynamic model parameters. The heat balance equation includes unknown parameters such as electrical loss coefficient, mechanical friction coefficient, hydraulic loss coefficient, and heat dissipation coefficient. This embodiment extracts steady-state segments from historical operating data and constructs a linear regression equation system Y=X*β, where Y is the measured temperature rise vector, X is the design matrix composed of physical quantities such as power squared, rotational speed, and flow rate, and β is the parameter vector to be identified. Since multicollinearity may exist among the physical quantities, ordinary least squares methods may lead to unstable parameter estimation (e.g., the appearance of physically meaningless negative values). Therefore, this embodiment uses Ridge Regression or Lasso Regression, adding a regularization term (such as the L2 norm ||β||) to the loss function. 2 By constraining the modulus of the parameters, stable and physically interpretable thermal equilibrium model parameters can be obtained.

[0172] Step 1105, the health level cloud model, is obtained by extracting the expected value, entropy, and hyperentropy parameters of each level based on historical sample data with calibrated health levels using the reverse cloud generator algorithm.

[0173] This embodiment describes the calculation process of the evaluation criteria. Experts categorize historical samples into levels such as excellent and good based on historical maintenance records and industry standards. For each sample set under level j, the cloud parameters are calculated using the inverse cloud generator algorithm. The specific formula is as follows: Expected value Ex j =(1 / N)*Σxi (Sample mean); Entropy En j =sqrt(π / 2)*(1 / N)*Σ|x i -Ex j |(Based on mean absolute deviation estimation); Hyperentropy He j =sqrt(S j 2 -En j 2 (Based on sample variance S) j (Relationship with entropy). Where, Ex j Let x be the expected value of the j-th health level; N be the sample size; x i For sample values; En j The entropy of the j-th health level; He j S represents the hyperentropy of the j-th health level; j 2 Let represent the sample variance. Through the above calculations, abstract expert experience is transformed into computer-processable numerical features (Ex, En, He), completing the transformation from human-based standards to a data-driven model.

[0174] In some embodiments, to further improve the adaptability of the model, an online update mechanism can be established. For example, every month, using newly accumulated manually verified sample data, the above-mentioned clustering and parameter identification process is re-executed to incrementally update the pre-set model or parameter library, adapting to the slow drift of characteristics caused by unit aging.

[0175] Step 1106, the method for constructing the component sensitive atom library is as follows:

[0176] Vibration signal samples corresponding to various typical faults (such as bearing inner ring wear, blade cracks, seal wear, and rotor imbalance) are collected from a historical fault case database. Each type of fault sample undergoes preprocessing, including detrending, filtering, and segmentation. All fault samples are used as a training set, and an iterative optimization algorithm based on K-Singular Value Decomposition (K-SVD) is employed: in the sparse coding stage, the atom library is fixed, and the orthogonal matching pursuit algorithm is used to solve for the sparse representation coefficients of each sample; in the dictionary update stage, the sparse coefficients are fixed, and atoms are updated column by column to minimize the reconstruction error. After iterative convergence, an original dictionary containing a predetermined number of atoms (e.g., 100 to 200) is obtained. Based on the activation frequency of each atom in various types of fault samples, a fault type label is assigned, forming a component-sensitive atom library.

[0177] Step 1107, the calibration method for the water level sensitivity factor is as follows:

[0178] Mutual information statistics under different water level difference conditions are extracted from historical operating data. A linear regression analysis is performed using the water level difference deviation rate (i.e., the relative difference between the current water level difference and the design water level difference) as the independent variable and the increase in the standard deviation of mutual information under the corresponding conditions relative to the design conditions as the dependent variable. The slope of the regression line is the water level sensitivity factor η. Typically, for pumping stations with large head variations, η ranges from 0.15 to 0.30.

[0179] Step 1108, the method for identifying the thermal time constant is as follows:

[0180] Periods of abrupt changes in operating conditions were identified from historical operational data. For each temperature measuring point, the temperature response curve for that period was extracted and fitted to the theoretical model of the step response of a first-order inertial element. The optimal thermal time constant τ was determined using the nonlinear least squares method, with the objective of minimizing the mean square error between the measured and theoretical curves. For measuring points exhibiting multiple step changes in operating conditions, the median of the multiple identification results was taken as the final thermal time constant.

[0181] Example 9: Supplementary description of the evaluation system index construction process, addressing how to unify and integrate the characteristic parameters of heterogeneous and multi-physics fields to form a standardized interface suitable for subsequent weight calculation and cloud model input.

[0182] Step 331: Unify and standardize the features of multi-source signals after specialized processing, and construct a standardized evaluation index set;

[0183] This step occurs after feature decoupling and before weight calculation. To achieve comprehensive coverage of the unit's status, this embodiment constructs a hierarchical index system. Specifically, it includes the following six subsets: Vibration index subset: In addition to component-sensitive vibration characteristics, it also includes the passband amplitude, vibration intensity (effective velocity value), and kurtosis index of the original vibration signal. Temperature index subset: This refers to the equivalent standard operating temperature, covering key measurement points such as the upper guide, lower guide, water guide bearings, and stator windings. Swing index subset: Extracts the peak-to-peak value and phase information of the swing signals from the upper guide, lower guide, and water guide components. Clearance index subset: Extracts the minimum value and rate of change of the runner chamber clearance (air shroud clearance) for monitoring blade scratch risk. Electrical index subset: Extracts the stator current imbalance, power factor deviation, and voltage harmonic distortion (THD) for monitoring motor electrical faults. Hydraulic indicators subset: including operating efficiency deviation (the difference between measured efficiency and theoretical efficiency), net positive suction head (NPSH) margin, and tailrace pressure pulsation coefficient.

[0184] Step 332: Establish a unified numbering system and attribute descriptions for all indicators to form a standardized evaluation indicator set Ind full ;

[0185] To facilitate automated computer processing, this embodiment addresses the above-mentioned n total Each indicator is uniformly coded (e.g., V). 01 Represents upper guide vibration, T 01 (Represents stator temperature). For each indicator, its metadata is recorded, including: indicator name, current measured value (or calculated value), physical unit, design normal range (for normalization reference), and source step identifier. The final output is a standardized evaluation indicator set Ind. full It is a structured data object that eliminates the differences between different physical dimensions and data sources, serving as a unified input interface for subsequent dual-source weight calculation and cloud model evaluation.

[0186] To address the false alarm problem caused by the physical modulation of monitoring indicators by hydraulic operating conditions, this embodiment adopts an adaptive threshold model driven by both operating conditions and water level. By superimposing a water level difference correction term on the mutual information statistical features, the association judgment threshold is dynamically adjusted according to the degree of deviation of the actual water level from the design water level, effectively filtering out false weak associations caused by hydraulic excitation and solving the problem of network construction distortion under non-design operating conditions.

[0187] To address the issue of early, subtle fault characteristics being masked by baseline drift under varying operating conditions, this embodiment employs a physical feature decoupling strategy. Regarding vibration, differential calculations are performed using the baseline mean curve for the corresponding operating condition to eliminate background noise from normal operating conditions. Regarding temperature, the equivalent standard operating temperature is calculated using the thermal balance equation and thermal inertia model, eliminating environmental and load interference. This approach improves the signal-to-noise ratio and enables the sensitive detection of subtle anomalies.

[0188] To address the conflict between statistical weights and mechanistic weights under abnormal operating conditions, this embodiment constructs a dual-source weight adaptive coupling mechanism based on rank correlation. By calculating the consistency coefficient between random forest weights and network topology weights, when the discrepancy between the two is significant (i.e., data patterns fail), the proportion of the more robust physical topology weights is automatically increased, ensuring the model's robustness in the face of unknown or complex faults.

[0189] The preferred embodiments of the present invention have been described in detail above. However, the present invention is not limited to the specific details in the above embodiments. Within the scope of the technical concept of the present invention, various equivalent transformations can be made to the technical solutions of the present invention, and these equivalent transformations all fall within the protection scope of the present invention.

Claims

1. A comprehensive evaluation method for the condition of pump station units, characterized in that, include: Multi-source monitoring data of pump station units are acquired and time-aligned and standardized to obtain a standardized monitoring dataset. Calculate the time-varying mutual information between various monitoring indicators in the standardized monitoring dataset, and determine the dynamic threshold based on the current operating condition parameters extracted from the standardized monitoring dataset to construct a dynamic correlation network; Decoupling the physical mechanism features of the standardized monitoring dataset includes: removing the working condition drift component of the vibration signal based on a pre-set working condition-vibration reference curve library to obtain vibration residual features, and converting the temperature signal into an equivalent standard working condition temperature based on the principle of thermal balance. The first weight based on the statistical features of the standardized monitoring dataset and the second weight based on the topological features of the dynamic association network are calculated respectively, and the consistency coefficient between the first weight and the second weight is calculated to generate a comprehensive coupling weight. Based on the comprehensive coupling weight, vibration residual characteristics and equivalent standard operating temperature, the health status level of the unit is determined by using a pre-built health level cloud model, and the evaluation results are obtained. Current operating parameters include the unit's current calculated water level difference and operating power; dynamic thresholds are determined based on current operating parameters extracted from standardized monitoring datasets, including: The current operating condition parameters are matched with the preset operating condition cluster centers to determine the current operating condition category; the mutual information benchmark statistical features corresponding to the operating condition category are queried. An adaptive threshold model incorporating a water level difference correction term is used to calculate the dynamic threshold. The adaptive threshold model is defined as follows: based on the mutual information baseline statistical characteristics of the corresponding working condition category, a water level difference correction term is superimposed. The water level difference correction term is positively correlated with the degree of deviation of the calculated water level difference from the design water level difference of the pumping station, and its correction strength is controlled by a preset water level sensitivity factor. Constructing dynamic interconnected networks includes: Calculate the time-varying mutual information among the monitoring indicators in the standardized monitoring dataset to obtain the time-varying mutual information matrix; The elements in the time-varying mutual information matrix are compared with the dynamic threshold element by element. When the mutual information value is greater than or equal to the dynamic threshold, a connection edge is established between the corresponding monitoring index nodes, and the mutual information value is used as the edge weight to generate a weighted adjacency matrix. Based on the weighted adjacency matrix, the degree centrality, betweenness centrality, and clustering coefficient of each monitoring index node are calculated as topological features of the dynamic association network. Vibration residual characteristics are obtained by stripping the operating condition drift component of the vibration signal based on a preset reference curve, including: The operating condition category of the unit is determined based on the current operating condition parameters, and the vibration mean curve corresponding to the operating condition category is retrieved from the preset operating condition-vibration reference curve library. The real-time vibration signal in the standardized monitoring dataset is subtracted point by point from the vibration mean curve to remove the vibration baseline determined by the working conditions, and a residual signal containing non-stationary abnormal components is obtained. The residual signal is subjected to energy normalization to obtain the vibration residual characteristics.

2. The method according to claim 1, characterized in that, Calculate the time-varying mutual information among the monitoring indicators in the standardized monitoring dataset, including: The standardized monitoring dataset is divided into sliding window segments to generate a continuous time window sequence; For the monitoring index data within each time window, the adaptive bandwidth is calculated based on the standard deviation and kurtosis coefficient of the sample, and the probability density is estimated using the Gaussian kernel function based on the adaptive bandwidth. Based on the estimated marginal probability density and joint probability density, the mutual information values ​​between each monitoring index pair are calculated to form a time-varying mutual information matrix sequence corresponding to the time window sequence.

3. The method according to claim 1, characterized in that, After obtaining the vibration residual features, the process also includes sparse fault feature extraction from the vibration residual features: Load a pre-built component-sensitive atom library, which contains multiple characteristic atoms associated with specific fault types; The orthogonal matching pursuit algorithm is used to perform sparse decomposition on the vibration residual characteristics and solve the sparse coefficient vector under the sensitive atom library of the component. Based on the atomic labels corresponding to the non-zero elements in the sparse coefficient vector, the vibration energy contribution values ​​of each component are classified and accumulated to generate a component sensitive vibration feature vector.

4. The method according to claim 1, characterized in that, Based on the principle of thermal balance, the temperature signal is converted into an equivalent standard operating temperature, including: Establish the heat balance equations for the pump station units, including heat source terms for electrical losses, heat source terms for mechanical friction, heat source terms for hydraulic losses, and heat dissipation terms. Substitute the current operating parameters from the standardized monitoring dataset into the heat balance equation to calculate the theoretical steady-state temperature under the current operating conditions; at the same time, substitute the preset rated operating parameters into the heat balance equation to calculate the reference steady-state temperature under the standard operating conditions. The difference between the theoretical steady-state temperature and the reference steady-state temperature is used to compensate for the real-time temperature monitoring values, and the equivalent standard operating temperature that eliminates the influence of operating condition differences is calculated.

5. The method according to claim 1, characterized in that, Calculate the first weight of the statistical features based on the standardized monitoring dataset and the second weight of the topological features based on the dynamic relational network, including: A random forest classification model was constructed using a standardized monitoring dataset and corresponding historical health status labels. The first weight of each monitoring indicator was determined by calculating the decrease in accuracy of feature permutation of out-of-bag data. The degree centrality, betweenness centrality, and clustering coefficient of each monitoring indicator node in the dynamic association network are extracted, and the above three indicators are weighted and fused according to the preset physical meaning importance weight coefficient to determine the second weight of each monitoring indicator.

6. The method according to claim 5, characterized in that, Calculate the consistency coefficient between the first weight and the second weight, and generate a comprehensive coupling weight based on this, including: The first and second weights are sorted separately to obtain the rank vector, and the Spearman rank correlation coefficient between the two is calculated as the consistency coefficient. The network topology fusion coefficient is adaptively calculated based on the consistency coefficient. The lower the consistency coefficient, the larger the value of the network topology fusion coefficient. When there is a divergence in the weight distribution, the contribution ratio of the weight based on the physical mechanism is increased. Based on the network topology fusion coefficient, the first weight and the second weight are weighted and summed to obtain the comprehensive coupling weight.

7. The method according to claim 1, characterized in that, The health status level of the unit is determined using a pre-built health level cloud model, including: The vibration residual characteristics, equivalent standard operating temperature, and swing characteristics, electrical characteristics, and hydraulic characteristics in the standardized monitoring dataset are input into the preset health level cloud model. The positive cloud generator is used to calculate the degree of membership of each indicator relative to different health levels. The membership certainty of a single indicator is weighted and summed based on the comprehensive coupling weight to generate a comprehensive membership certainty vector that reflects the distribution probability of the overall status of the unit at each health level. The current health status level is determined based on the maximum value in the comprehensive membership degree vector, and the distribution entropy and concentration of the comprehensive membership degree vector are calculated to assess the reliability of the determination result.

Images