Pump station unit status trend prediction and early warning method

Abstract

This invention discloses a method for predicting and warning the status trend of pumping station units. The method acquires real-time multi-source operational monitoring data of the pumping station units and inputs it into a pre-trained deep learning trend prediction model to obtain the unit status trend prediction result. Based on the prediction result and monitoring data, a multi-factor coupling sensitivity matrix is ​​calculated. The diagonal elements of this matrix represent the independent influence of individual monitoring data, while the off-diagonal elements represent the coupled interactive influence between different monitoring data. This matrix is ​​then input into a preset sensitivity pattern rule base for matching and retrieval to determine the sensitivity pattern to which the current operating state belongs. When the sensitivity pattern meets the conditions, an early warning signal is generated. This invention effectively solves the problems of poor interpretability of deep prediction models and difficulty in calculating multivariate coupling effects by introducing physical causal constraints and coupling sensitivity analysis, achieving accurate early warning with causal analysis capabilities.

Description

Technical Field

[0001] This invention relates to the field of intelligent monitoring and early warning technology for electromechanical equipment, and in particular to a method for predicting and warning the status trend of pump station units. Background Technology

[0002] Large pumping station units are complex systems involving the coupling of multiple physical fields, including hydraulic, mechanical, and electromagnetic fields. Currently, using multi-source sensor data (such as vibration, temperature, water level, and blade angle) for full lifecycle status perception and trend analysis of the units has become a key means of achieving predictive maintenance. This technology aims to identify equipment performance degradation trends in advance by mining the nonlinear characteristics in massive time-series data, thereby avoiding the risk of downtime caused by sudden failures.

[0003] Currently, deep learning algorithms are widely used in pump station unit status prediction. Existing mainstream technical solutions typically employ architectures such as Convolutional Neural Networks (CNN) or Long Short-Term Memory Networks (LSTM) to construct end-to-end prediction models. These methods learn the spatiotemporal dependencies of historical operating data to establish a mapping model between input monitoring indicators and output status indicators, directly outputting future vibration amplitude or temperature trend curves. In engineering applications, the system typically compares the numerical sequence predicted by the model with statically set safety thresholds. When the predicted value exceeds the threshold, an alarm is triggered, and the prediction results are displayed to maintenance personnel in the form of a trend graph.

[0004] However, existing technologies, when dealing with the complex operating conditions of pumping stations, mainly suffer from a lack of interpretability due to the black-box nature of deep learning models, and a technical problem of separating physical causal mechanisms from numerical algorithms. Specifically, while existing models can provide high-precision predictions, they cannot reveal the specific causes leading to these prediction trends. Standard sensitivity analysis methods (such as simple gradient or perturbation methods) typically assume that each input feature is independent, ignoring the objectively existing physical causal chains in the pumping station system (e.g., water level changes precede blade adjustments, thus affecting vibration), making it impossible to accurately isolate the coupling and interaction effects between multiple variables. This deficiency makes it difficult to transform continuous predicted values ​​into discrete decision rules with clear physical meaning, resulting in maintenance personnel struggling to quickly locate the source of the fault upon receiving warnings, and the early warning system's decision support capability being insufficient. Summary of the Invention

[0005] Purpose of the invention: To provide a method for predicting and warning the status trend of pump station units, in order to solve the above-mentioned problems in the prior art.

[0006] Technical solution: A method for predicting and warning the status trend of pump station units, comprising:

[0007] Real-time multi-source operation monitoring data of the pump station unit is acquired, and the real-time multi-source operation monitoring data is input into a pre-trained deep learning trend prediction model to obtain the unit status trend prediction result.

[0008] Based on the unit status trend prediction results and real-time multi-source operation monitoring data, a multi-factor coupling sensitivity matrix is ​​calculated.

[0009] The multi-factor coupled sensitivity matrix is ​​input into a preset sensitivity mode rule base for matching and retrieval to determine the sensitivity mode to which the current operating state belongs.

[0010] When the warning triggering conditions are met in the sensitivity mode, a warning signal is generated according to the matched rules.

[0011] Beneficial effects: By introducing physical causal constraints and coupling sensitivity analysis, this invention effectively solves the problems of poor interpretability of deep prediction models and difficulty in calculating multivariate coupling effects, and realizes accurate early warning with causal analysis capabilities. Attached Figure Description

[0012] Figure 1 This is a flowchart illustrating the accurate prediction and early warning of unit status based on a multi-factor coupled sensitivity matrix in an embodiment of this application.

[0013] Figure 2 This is a flowchart illustrating the steps for obtaining the multi-factor coupling sensitivity matrix in the embodiments of this application.

[0014] Figure 3 This is a flowchart illustrating the steps for generating a warning signal in an embodiment of this application.

[0015] Figure 4 This is a flowchart illustrating the steps involved in generating a rule base in an embodiment of this application. Detailed Implementation

[0016] Example 1 describes the process of achieving accurate prediction and early warning of unit status based on a multi-factor coupled sensitivity matrix, such as... Figure 1 As shown.

[0017] Step 101: Obtain real-time multi-source operation monitoring data of the pump station unit, and input the real-time multi-source operation monitoring data into the pre-trained deep learning trend prediction model to obtain the unit status trend prediction result.

[0018] In this embodiment, real-time multi-source operation monitoring data refers to the collection of physical quantities acquired in real time by sensor arrays deployed at key locations in the pumping station. This type of data includes not only unit body status parameters such as vibration (e.g., X-axis, Y-axis, Z-axis acceleration), temperature (e.g., stator winding temperature, bearing temperature), current (e.g., three-phase current), and voltage, but also hydraulic and control parameters such as inlet pool water level, outlet pool water level, blade angle, and guide vane opening. Specifically, the above data can be continuously acquired via industrial Ethernet, fieldbus (e.g., PROFIBUS, Modbus), or wireless sensor networks at a preset sampling frequency (e.g., 1Hz for status monitoring data and 1kHz for vibration signals).

[0019] The deep learning trend prediction model is a nonlinear regression model pre-trained on a large amount of historical data. Its input is current real-time multi-source operation monitoring data (usually a sequence of data within a historical time window), and its output is the predicted value of unit status indicators for a future period (such as the vibration amplitude trend or temperature rise trend over the next 12 time steps). This model can uncover the spatiotemporal correlation features hidden in the multi-source data, providing a high-precision prediction benchmark for subsequent sensitivity analysis. Preferably, before inputting the data into the model, the real-time multi-source operation monitoring data is standardized (e.g., Z-score standardization) to eliminate the influence of different physical dimensions.

[0020] Step 102: Calculate the multi-factor coupling sensitivity matrix based on the unit status trend prediction results and real-time multi-source operation monitoring data.

[0021] In this embodiment, the multi-factor coupling sensitivity matrix is ​​a structured causal relationship representation tool, not a simple correlation coefficient matrix, but an analytical result obtained through gradient backpropagation or integral calculation of a deep learning trend prediction model. The rows and columns of this matrix correspond to the various input monitoring indicators. Specifically, this matrix is ​​an n×n symmetric matrix (where n is the dimension of the input features), containing two key layers of information: first, the diagonal elements, representing the independent contribution of a single monitoring data point (e.g., water level) to the prediction result (e.g., vibration trend), i.e., single-factor sensitivity; and second, the off-diagonal elements, representing the nonlinear interaction effect generated when two different monitoring data points (e.g., water level and blade angle) change simultaneously, i.e., two-factor coupling sensitivity.

[0022] By calculating this matrix, this embodiment makes the prediction logic of the black-box model transparent and calculates the causal contribution of each input factor to the output result. For example, when the element value corresponding to water level-vibration in the matrix increases, it indicates that the current predicted vibration trend is mainly caused by changes in water level.

[0023] Step 103: Input the multi-factor coupled sensitivity matrix into the preset sensitivity mode rule base for matching and retrieval to determine the sensitivity mode to which the current operating state belongs.

[0024] The preset sensitivity pattern rule base is a knowledge base that stores a large number of historical typical operating conditions (including normal, warning, fault, and other states) and their corresponding sensitivity feature patterns. Each rule describes a specific combination of sensitivity patterns, such as high water level sensitivity and extremely strong water level-vibration coupling sensitivity. In this embodiment, by calculating the matching degree between the multi-factor coupling sensitivity matrix calculated at the current moment and each rule in the rule base, the underlying causal patterns of the current unit operating state can be identified.

[0025] Specifically, this type of matching retrieval can employ similarity calculation or logical condition judgment. For example, the system checks whether the current sensitivity matrix contains a combination of salient features defined in the rule base, determining which known sensitivity mode the current operating state belongs to. This process realizes the transformation from numerical calculation to logical diagnosis.

[0026] Step 104: When the sensitivity mode meets the warning triggering conditions, generate a warning signal according to the matched rules.

[0027] Once the current sensitivity mode is determined, the system checks whether the rules corresponding to that mode are associated with a warning level. If the matched rule is defined as having potential risk (e.g., associated with a Level 1 or Level 2 warning label), the warning mechanism is immediately triggered. The generated warning signal not only includes alarm prompts (such as audible and visual alarms, SMS push notifications) but also provides interpretable auxiliary information based on the content of the matched rules. For example, the warning signal might indicate that the current vibration prediction value is abnormal, mainly due to flow instability caused by fluctuations in the inlet pool water level, and recommends checking the forebay water level control strategy. This type of warning signal, which includes causal analysis, provides direct decision support for maintenance personnel.

[0028] The following example of a typical early warning scenario illustrates the specific execution process of this method.

[0029] Scenario setting: A large axial flow pump station is operating during the flood season. The current water level in the inlet pool is 8.5 meters (7.0 meters higher than the normal water level), the unit vibration value is 4.2 mm / s (close to the warning threshold of 5.0 mm / s), and the bearing temperature is 58℃ (within the normal range).

[0030] Step 1: The system collects the above real-time monitoring data, inputs it into the deep learning model, and predicts the vibration trend in the next 12 time steps. The prediction results show that the vibration value will reach 5.8 mm / s in the 8th time step.

[0031] Step 2: The system calculates the multi-factor coupling sensitivity matrix. The results show that the single-factor sensitivity of water level is 0.42, the single-factor sensitivity of vibration is 0.18, and the coupling sensitivity of water level and vibration is 0.25, indicating that water level change is the main reason for the upward trend of vibration.

[0032] Step 3: Match the sensitivity matrix with the rule base. The rule with the highest soft matching degree is the water level-dominated vibration anomaly mode, with a matching degree of 0.85, corresponding to the level 2 warning.

[0033] Step 4: The system generates an early warning signal, which reads: It is predicted that the unit vibration will exceed the safety threshold in about 8 minutes. The main cause is the unstable flow caused by the abnormal rise in the water level of the inlet pool. It is recommended to check the water level control of the forebay and appropriately reduce the operating power of the unit.

[0034] In actual operation, this method also includes the following exception handling mechanisms:

[0035] Sensor fault handling: When an abnormality is detected in the data of a certain monitoring indicator (such as a sudden change in value exceeding the threshold, data remaining at zero or data loss), the system adopts the strategy of filling with the most recent valid value or replacing redundant data from multiple sensors. At the same time, the indicator is downweighted when calculating the sensitivity matrix, and a sensor abnormality alert is sent to the operation and maintenance personnel.

[0036] No matching rule handling: When the soft matching degree of the sensitivity matrix with all rules in the rule base is lower than the minimum confidence threshold (e.g., 0.3), the system determines that the current state is an unknown mode, triggers the manual review process, and records the sample to the queue to be analyzed for subsequent rule base expansion.

[0037] Handling of identical priority scores: When multiple rules have the same priority score, the rule with the newest creation time is selected first. If the creation time is also the same, the rule with the higher warning level is selected to ensure safety.

[0038] Example 2 describes the specific network structure, parameter configuration, and data processing details of the deep learning trend prediction model.

[0039] Step 201: The deep learning trend prediction model includes a convolutional neural network layer, a bidirectional long short-term memory network layer, and a temporal attention mechanism layer connected in sequence.

[0040] The convolutional neural network layer is used to extract the local correlation features between various monitoring indicators in real-time multi-source operation monitoring data; the bidirectional long short-term memory network layer is used to extract the forward and backward temporal dependency information of the monitoring data based on the local correlation features; the temporal attention mechanism layer is used to assign weights to the temporal dependency information at different historical moments, thereby enhancing the contribution of key historical information to the prediction results of unit status trends.

[0041] In this embodiment, the model adopts a serial hybrid deep neural network architecture. Real-time multi-source operation monitoring data is input in the form of time series, and the dimension of the input data can be represented as (B,T,N), where B is the batch size, T is the time window length (e.g., 60 time steps), and N is the number of features (e.g., 16 monitoring indicators).

[0042] Specifically, the convolutional neural network (CNN) layer is mainly used to extract local spatial correlation features of various monitoring indicators. Preferably, this layer contains two one-dimensional convolutional layers (1D-CNN, 1D Convolutional Neural Network). The first layer has 64 kernels, a kernel size of 3, and a stride of 1, used to capture primary correlations between adjacent time steps and adjacent features; the second layer increases the number of kernels to 128, while maintaining the kernel size of 3, used to further extract deeper combined features. Each convolutional layer is followed by a Rectified Linear Unit (ReLU) activation function, and a max pooling layer (with a pooling kernel size of 2) is used to reduce the data dimensionality and computational cost.

[0043] A Bidirectional Long Short-Term Memory (BiLSTM) network layer is connected after the CNN layer to further explore the long-term dependencies of data over time based on the feature vector sequences extracted by the CNN. Since the operating state of the pumping station is not only affected by the current moment but also closely related to historical states, BiLSTM, by simultaneously training two LSTM networks (forward and backward), can more comprehensively capture the temporal evolution patterns. Preferably, this embodiment uses a two-layer stacked BiLSTM structure, with each layer containing 128 hidden units, and a dropout rate of 0.2 is set between layers to prevent overfitting.

[0044] Step 202: The temporal attention mechanism layer is used to assign weights to temporal dependency information at different historical moments, thereby enhancing the contribution of key historical information to the prediction results of unit status trends.

[0045] The Temporal Attention layer is located after the BiLSTM layer and before the output layer. Since not every data point in a long input sequence (e.g., 60 time steps) is equally important for future predictions, this layer employs a Scaled Dot-Product Attention mechanism to automatically calculate the weight coefficients for each time step. Specifically, the model generates a query vector based on the current hidden state and matches it with the key vectors from each historical time step to obtain an attention score. After normalization using the Softmax function, this score is used as weights to perform a weighted summation of the value vectors from each historical time step, generating the final context vector. This mechanism allows the model to focus on the key historical moments (such as abrupt changes or turning points) that have the greatest impact on the prediction results, improving prediction accuracy.

[0046] The model's output layer is a fully connected layer, mapping the context vector output by the attention layer to the final prediction target dimension. For example, if the prediction target is the root mean square value of vibration over the next 12 time steps, the output layer dimension is 12. During model training, mean squared error (MSE) is preferably used as the loss function, employing the Adam adaptive moment estimation optimizer (with an initial learning rate of 0.001), and incorporating an early stopping strategy. This strategy automatically stops training when the validation set loss no longer decreases within 10 consecutive training epochs, thus obtaining the optimal model parameters.

[0047] Example 3 describes the construction of integral paths under physical constraints, which addresses the shortcomings of traditional interpretable algorithms that ignore physical causal logic.

[0048] Step 301: Obtain the preset physical causal order between various monitoring indicators of the pumping station.

[0049] In this embodiment, the physical causal order refers to an ordered list or directed acyclic graph structure that defines the order of influence transmission among various monitoring indicators of the pumping station. Determining this order is a prerequisite for subsequently constructing the integral path. This embodiment provides three specific implementation methods for obtaining this order.

[0050] The first approach uses a pre-defined physical model. Specifically, based on the energy and signal flow of the pumping station system, the monitoring indicators are divided into five levels: the first level is the external input level, including the inlet water level, ambient temperature, and grid voltage; the second level is the control command level, including start / stop commands and frequency converter settings; the third level is the motion conversion level, including motor speed and motor current; the fourth level is the fluid transport level, including flow rate and outlet pressure; and the fifth level is the state response level, including bearing temperature and unit vibration. In this sequence, indicators with lower level numbers are defined as upstream variables of indicators with higher level numbers. For example, changes in the inlet water level are causally preceded or synchronized with changes in unit vibration.

[0051] The second approach utilizes the PC algorithm to perform conditional independence checks on historical operational data, constructs an undirected graph and orients it as a directed acyclic graph (DAG), and determines the causal order based on topological sorting. This approach is suitable for scenarios with complex physical mechanisms. The specific steps include: constructing an undirected graph containing all monitoring index nodes; performing conditional independence checks (such as the Fisher-Z test) on each pair of adjacent nodes in the undirected graph; if the two nodes are independent given a certain set of conditions, the connecting edge is deleted; next, identifying the V-shaped structure in the graph (i.e., the XZY structure, where X and Y are not adjacent); if Z is not in the separating set of X and Y, orienting X to Z and Y to Z; applying Meek rules (such as rules to prevent loops) to orient the remaining undirected edges; and performing topological sorting on the generated DAG to output a linear causal order.

[0052] The third approach utilizes the Granger causality test to calculate the time-lag cross-correlation between the time series of various monitoring indicators, and determines the causal order based on the lead-lag score. Specifically, for any two variable time series, the cross-correlation coefficient is calculated at different lag orders. If historical data of variable A can improve the prediction accuracy of variable B, and the optimal lag time is positive, then A is determined to lead B. The causal order can be obtained by calculating the cumulative score of each variable leading all other variables and sorting them in descending order of score.

[0053] Step 302: Determine the baseline state of the real-time multi-source operation monitoring data, and construct a causal constraint integral path connecting the baseline state and the real-time multi-source operation monitoring data according to the physical causal order.

[0054] In this embodiment, the baseline state is typically defined as the average or median vector of various indicators in the historical data of the pump station's normal operation, denoted as x'. The real-time multi-source operation monitoring data at the current moment is denoted as x. t The goal of constructing an integral path is to find a path in high-dimensional space pointing from x' to x. t The curve.

[0055] Step 303: The causal constraint integral path is defined by a set of interpolation variables, which satisfy the causal monotonicity constraint. The causal monotonicity constraint means that the rate of change of the interpolation variable corresponding to the monitoring indicator located upstream in the physical causal sequence is always greater than or equal to the rate of change of the interpolation variable corresponding to the monitoring indicator located downstream.

[0056] Specifically, suppose the system has n variables, arranged in causal order as x1, x2, ..., xn. n (Where x1 is the upstream). Any point x(α) on the defined integration path can be represented as a linear combination of the baseline state and the change, but each dimension i has an independent interpolation coefficient α. i The causal monotonicity constraint mathematically requires the interpolation coefficient vector (α1, α2, ..., α...) to... n The expression must fall within a specific simplex region, i.e., satisfy 1 ≥ α1 ≥ α2 ≥ ... ≥ α n ≥0. During the transition from the baseline state to the current state, upstream variables (such as water level) always begin to change before or simultaneously with downstream variables (such as vibration). There will be no path where the vibration changes to the current value while the water level remains at the baseline value.

[0057] To numerically implement this path, this embodiment employs an ordered sampling algorithm to generate path points. The specific steps are as follows: generate n independent random numbers following a uniform distribution from 0 to 1; sort these n random numbers in descending order; and directly assign the sorted random number sequence to interpolation coefficients α1 to α2 arranged in causal order. n Each sample point generated in this way satisfies the aforementioned causal monotonicity constraint, eliminating the need for complex rejection sampling and improving computational efficiency.

[0058] Step 304: Perform gradient integral calculation on the deep learning trend prediction model along the causal constraint integral path to obtain the multi-factor coupling sensitivity matrix, such as... Figure 2 As shown.

[0059] In this embodiment, m sample points (e.g., m=50) are collected along the constructed causal path. The gradient of the model at these points is calculated (for single-factor problems) or the second-order mixed partial derivative (for two-factor problems), and the average value is obtained to obtain an approximate value of the integral gradient. This method, by introducing physical priors, eliminates noise interference caused by non-physical paths (such as straight paths) in traditional integral gradient methods, making the calculated sensitivity more reflective of the real fault propagation mechanism.

[0060] In a preferred embodiment of the present invention, the integral gradient calculation adopts a causal constraint integral domain to ensure that the integral path conforms to the physical causal laws of the pumping station system.

[0061] The causal constraint integration domain ωπ is defined as a subset within an n-dimensional unit hypercube that satisfies the causal order constraint. Specifically, this integration domain is defined as: ωπ = (α1, α2, ..., αn): 0 ≤ αk ≤ 1;

[0062] It holds for all k, and α σ (1)≧α σ (2)≧...≧α σ (n);

[0063] Where σ is the reverse permutation of the causal order π, and σ(k) represents the index of the original variable corresponding to the k-th position in the causal chain.

[0064] The causal constraint integration domain is geometrically an n-dimensional simplex with the following volume:

[0065] Volume(ω π ) = 1 / n!;

[0066] Compared to the full-space integration used in traditional integral gradient methods, the causal constraint integration domain occupies only 1 / n! of the full space, but it includes all integration paths that conform to the laws of physical causality.

[0067] The boundary of the causal constraint integral domain consists of the following parts: the first type of boundary is the set of points where the interpolation coefficients of adjacent causal level variables are equal, representing the boundary case where adjacent level variables change synchronously; the second type of boundary is the set of points where the interpolation coefficient of the downstream variable is zero, representing the case where the causal chain has not yet been passed to the downstream variable; the third type of boundary is the set of points where the interpolation coefficient of the upstream variable is 1, representing the case where the upstream variable has changed to the current state.

[0068] To perform numerical integration within the causal-constrained integration domain, an ordered sampling method is used to generate uniformly distributed sampling points. This method includes the following steps:

[0069] Step S3a: Generate n independent random numbers that follow a uniform distribution from 0 to 1, denoted as u1, u2, ..., u n .

[0070] Step S3b: Sort the n random numbers in descending order of their values ​​to obtain an ordered sequence v1, v2, ..., v n The condition is satisfied that v1 ≥ v2 ≥ ... ≥ v n The sorting operation is represented as:

[0071] (v1,v2,...,v n =Sort descending (u1,u2,...,u n );

[0072] Step S3c: Assign interpolation coefficients according to causal order. For the variable at position k in the causal chain, whose original index is σ(k), the interpolation coefficient of this variable is assigned the following values:

[0073] α σ (k)=v k ;

[0074] The mathematical principle of the ordered sampling method is that the descending statistics of n independent and identically distributed uniform random variables follow a uniform distribution on the corresponding ordered simplex. The sampling points generated by this method are uniformly distributed within the causal constraint integration domain and can be used for subsequent Monte Carlo integration calculations.

[0075] Based on the above sampling method, the Monte Carlo method is used to calculate the causal constraint integral. The specific process is as follows:

[0076] Step S3d: Determine the number of samples M. Choose an appropriate number of samples based on the number of variables n: when n is between 3 and 5, M is 1000; when n is between 6 and 10, M is 5000; when n is between 11 and 20, M is 10000; when n exceeds 20, M is 50000. The above sample numbers ensure that the relative error at a 95% confidence level is less than 3%. The relationship between the number of samples and the number of variables can be expressed as:

[0077] M=1000, when 3≤n≤5;

[0078] M=5000, when 6≤n≤10;

[0079] M=10000, when 11≤n≤20;

[0080] M=50000, when n>20.

[0081] Step S3e: Call the ordered sampling method to generate M sampling points, denoted as α. m , where m ranges from 1 to M.

[0082] Step S3f: For each sampling point α m Construct the corresponding interpolation state vector x m :

[0083] x m =x ref +α m *Δ x ;

[0084] Where * denotes element-wise multiplication, Δ x =x current -x ref This is the state change vector. Calculate the gradient vector Grad of the prediction model at this point. F (x m) and Hessian matrix H F (x m ).

[0085] Step S3g: Average the gradient and Hessian values ​​of all sampling points, and multiply by n! for normalization to obtain the final sensitivity vector and coupled sensitivity matrix:

[0086] S i =n!*Δx i *(1 / M)*Σm=1 to M;

[0087] S ij =n!*Δx i *Δ xj *(1 / M)*Σ[m=1 to M](d2F / dx i *dx j |x m *w(α m )).

[0088] Example 4 describes the specific calculation formula and normalization method of the two-factor coupling sensitivity matrix, and in particular corrects the defect of the conventional method that ignores the collaborative changes of multiple variables.

[0089] Step 401: The multi-factor coupling sensitivity matrix is ​​a symmetric matrix, wherein: the diagonal elements of the multi-factor coupling sensitivity matrix represent the degree of independent influence of a single monitoring data on the prediction result of the unit status trend; the off-diagonal elements of the multi-factor coupling sensitivity matrix represent the degree of two-factor coupling interaction between two different monitoring data.

[0090] In this embodiment, the multi-factor coupling sensitivity matrix is ​​denoted as M. The diagonal elements of the matrix are denoted as M0. ii For the single-factor sensitivity corresponding to the i-th monitoring indicator, the direct contribution of the change in only the indicator itself to the prediction result, considering the causal path, was calculated. The off-diagonal elements of the matrix are denoted as M. ij (where i≠j) corresponds to the two-factor coupling sensitivity of the i-th and j-th indicators.

[0091] Step 402: The degree of two-factor coupling interaction between two different monitoring data consists of a direct interaction term and a causal transit term. The direct interaction term represents the integral contribution of the second-order mixed partial derivative calculated along the causal constraint integral path, and is used to calculate the nonlinear synergistic effect generated by the simultaneous change of the two monitoring data. The causal transit term represents the indirect influence contribution of the upstream monitoring data in the causal chain to the downstream monitoring data through intermediate variables.

[0092] Specifically, to more accurately calculate the direct interaction term, this embodiment employs a generalized interpolation path. The method for calculating the nonlinear synergistic effect between two variables on the generalized interpolation path is as follows:

[0093] ;

[0094] Among them, S ij direct Let x be the direct interaction term in the two-factor coupling sensitivity between the i-th and j-th monitoring indicators; F(·) is the pre-trained deep learning trend prediction model function; x ij (α,β) is the generalized two-parameter interpolation path function, which defines the trajectory of change from the baseline state to the current state; α is the interpolation coefficient (integral variable) corresponding to the i-th monitoring indicator; β is the interpolation coefficient (integral variable) corresponding to the j-th monitoring indicator. The model output is relative to the input variable x. i and x j The second-order mixed partial derivatives (Hessian matrix elements); w(α,β) is a weighting function defined according to the physical causal order. Unlike the simplified path that changes only two variables, this path allows for analysis of variable x... i and x j When the interaction occurs, the remaining variables x k Also based on its relationship with x i ,x j The dependencies undergo corresponding coordinated changes. Direct interaction item S ij_direct The specific calculation formula can be expressed as the weighted integral of the change in the variable and the second-order mixed partial derivative along the path.

[0095] In the specific implementation, the calculation involves the second-order mixture partial derivatives (i.e., Hessian matrix elements) of the model output F with respect to the input x. For deep learning models that support automatic differentiation (such as models built on deep learning frameworks that support automatic differentiation), the automatic differentiation interface provided by the framework can be used for efficient calculation. For black-box models that do not support automatic differentiation, this embodiment provides an alternative based on finite differences. Specifically, the central difference method is used, employing four perturbation points to approximate the calculation of the second-order mixture partial derivatives, with the following formula:

[0096] ;

[0097] in, is the approximation of the second-order mixed partial derivative of the model at the current state point; F(·) is the predicted output value of the deep learning trend prediction model; x ++ The perturbation input vector is the one whose i-th component is increased by hi and whose j-th component is increased by hj, based on the current state vector; x +−To increase h in the i-th component based on the current state vector i The j-th component minus h j The perturbation input vector after; x −+ To subtract h from the i-th component based on the current state vector. i The j-th component increases by h j The perturbation input vector after; x −− To subtract h from the i-th component based on the current state vector. i The j-th component minus h j The perturbation input vector after h; i h is the adaptive differential step size for the i-th monitoring indicator; j This is the adaptive differential step size for the j-th monitoring indicator.

[0098] For example, for variable x i and x j At the current point, a small step size h (e.g., 10 to the power of -4) is applied for both positive and negative perturbations. The output value of the fourth-order model is calculated, and the mixed partial derivative is obtained by calculating the difference formula. To ensure numerical stability, the step size h can be adaptively adjusted according to the magnitude of the variable's value; that is, the step size is equal to the base step size multiplied by the maximum value of the absolute value of the variable and 1.

[0099] In one embodiment of the present invention, the calculation of the second-order mixed partial derivative can be performed using an automatic differentiation method or a finite difference method.

[0100] When the prediction model is implemented based on a deep learning framework that supports automatic differentiation, the automatic differentiation method is preferred. This method calculates the second derivative through two backpropagations: it calculates the first gradient of the predicted output with respect to the input variables, and then differentiates each component of the first gradient to obtain the Hessian matrix. The automatic differentiation method has high computational accuracy and no truncation error.

[0101] When the prediction model does not support automatic differentiation or is a black-box model, the finite difference method is used to calculate the second-order partial derivatives. This method includes the following steps:

[0102] Step S4a: Determine the adaptive step size based on the numerical magnitude of each variable. For variable i, the step size calculation formula is:

[0103] h i =h0*max(|x i |,1);

[0104] Here, h0 is the basic step size, with a value of 1e-4. The adaptive step size can avoid numerical accuracy problems caused by differences in the magnitude of the dependent variable.

[0105] Step S4b: Calculate the second-order mixed partial derivatives using the central difference scheme. For the second-order mixed partial derivatives of variables i and j (i≠j), the formula is:

[0106] d2F / dx i *dx j =(F(x pp )-F(x pm )-F(x mp )+F(x mm )) / (4*h i *h j );

[0107] The perturbation points are defined as follows:

[0108] x pp =x+h i *e i +h j *e j ;

[0109] x pm =x+h i *e i -h j *e j ;

[0110] x mp =xh i *e i +h j *e j ;

[0111] x mm =xh i *e i -h j *e j .

[0112] e i and e j These are the unit vectors in the i-th and j-th directions, respectively.

[0113] Step S4c: For the diagonal elements, i.e., the second-order partial derivatives of variable i with respect to itself, the following formula is used for calculation:

[0114] d2F / dx i 2 =(F(x+h i *e i )-2*F(x)+F(xh i *e i )) / h i 2 ;

[0115] Step S4d: The number of times the prediction model needs to be called to calculate the complete Hessian matrix is:

[0116] N calls =1+2n+4n*(n-1) / 2=2*n 2 +1;

[0117] Here, n represents the number of variables. For a typical pump station monitoring scenario, there are 15 to 20 variables, and a single Hessian matrix calculation requires 450 to 800 model calls. Assuming a single model call takes 1 millisecond, the calculation of the complete coupled sensitivity matrix can be completed within seconds, meeting the requirements for real-time early warning.

[0118] Step 403: Calculate the sum of the absolute values ​​of all elements in the multi-factor coupling sensitivity matrix as a normalization factor; divide each element in the multi-factor coupling sensitivity matrix by the normalization factor to obtain the relative sensitivity matrix, and use the relative sensitivity matrix for subsequent rule matching.

[0119] The formula for calculating the relative sensitivity matrix, which aims to eliminate the influence of different physical dimensions, is as follows:

[0120] ;

[0121] Among them, M' ij M represents the element value in the i-th row and j-th column of the normalized relative sensitivity matrix; ij The value of the element in the i-th row and j-th column of the original multi-factor coupling sensitivity matrix; |M ij | represents the absolute value of the original sensitivity value; n represents the total dimension (number of features) of the monitoring indicators; ∑∑|M pq | represents the sum of the absolute values ​​of all elements in the original sensitivity matrix (normalization factor).

[0122] Because pump station monitoring data contains multiple physical dimensions (such as vibration measured in millimeters per second, temperature in degrees Celsius, and water level in meters), directly calculated sensitivity values ​​may vary significantly in magnitude, hindering the development of unified rules. Therefore, normalization is necessary. The specific operation involves taking the absolute value of each element in the calculated original sensitivity matrix M and summing all the absolute values ​​to obtain a scalar sum, i.e., the normalization factor. Each element in matrix M is then divided by this normalization factor.

[0123] After this processing, the sum of the absolute values ​​of all elements in the resulting relative sensitivity matrix is ​​1, and the value of each element is between -1 and 1, making the sensitivities with different physical meanings comparable. For example, if the normalized relative sensitivity of water level is 0.4 and the relative sensitivity of vibration is 0.1, it can be determined that water level change is the main cause of the current state. This type of relative sensitivity matrix serves as the standard input for the subsequent rule matching module.

[0124] Example 5 describes the automatic extraction of interpretable sensitivity rules from massive historical data, which solves the technical problem that manually defined rules are difficult to cover complex coupled operating conditions.

[0125] Step 501: Collect historical operating data of the pumping station, calculate the multi-factor coupled sensitivity matrix corresponding to each historical moment, and flatten it into a historical sensitivity sample vector.

[0126] In this embodiment, the first step in constructing the rule base is to transform continuous, high-dimensional historical operational data into a structured sensitivity feature space. A representative historical time period (e.g., operational data from the past year) is selected, covering normal operation, various faults, and different water level conditions. For each sampling time t within this time period, the corresponding multi-factor coupled sensitivity matrix M is calculated. t To facilitate subsequent cluster analysis, matrix M will be... t Perform a flattening operation. Because matrix M... t It is a symmetric matrix, requiring only the upper or lower triangular elements and diagonal elements to be retained. If the number of monitoring indicators is n, then the flattened historical sensitivity sample vector v t The dimension of the vector is n plus n multiplied by (n-1) divided by 2. For example, when n=10, the vector dimension is 55. This vector fully preserves the numerical characteristics of all single-factor influences and two-factor coupled influences at that moment.

[0127] Step 502: Perform cluster analysis on the historical sensitivity sample vector to identify several cluster centers that represent typical operating states.

[0128] In this embodiment, an unsupervised learning algorithm is used to perform pattern mining on the collected massive vector set. Preferably, the Ward linking method in hierarchical clustering is used. This method constructs a clustering tree by minimizing the intra-cluster variance, which can effectively discover irregularly shaped clusters. To determine the optimal number of clusters K, this embodiment uses the silhouette coefficient as an evaluation metric. The average silhouette coefficient corresponding to different K values ​​(e.g., K between 8 and 20) is calculated, and the K corresponding to the maximum value of the coefficient is selected as the final number of clusters. After clustering, the arithmetic mean of all sample vectors within each cluster is calculated to obtain the cluster center vector c of that cluster. k Each cluster center mathematically represents a typical pump station operation mechanism.

[0129] Step 503: Extract significant sensitivity features from each cluster center. Significant sensitivity features include single-factor sensitivity components or two-factor coupled sensitivity components whose values ​​exceed a threshold.

[0130] Since the cluster center vector contains information from all dimensions, feature filtering is required to generate concise and readable rules. In this embodiment, an adaptive statistical threshold is used for filtering. Specifically, the median or 75th percentile of the absolute values ​​of all elements in the cluster center vector is calculated and set as the significance threshold τ. The cluster center vector c is then traversed. k For each component, if the absolute value of a component is greater than the threshold τ, it is marked as a significant sensitivity feature. For example, if the value of the 3rd component of the cluster center (corresponding to the single-factor sensitivity of vibration) is 0.5 and the threshold is 0.2, then the component is retained; if the value of the 15th component (corresponding to the coupling sensitivity of temperature and water level) is 0.05, which is less than the threshold, then it is ignored.

[0131] Step 504: The significant sensitivity feature combination is transformed into the corresponding rule antecedent logic, and the risk level corresponding to the cluster center is associated as the rule consequent to generate a rule base, such as... Figure 4 As shown.

[0132] In this embodiment, the selected numerical features are converted into logical semantics. For significant single-factor sensitivity S... i If its value is positive and greater than the threshold, the generation rule condition S is... i If it is high; if it is negative and its absolute value is greater than the threshold, generate S. i The sensitivity is low for significant two-factor coupling. ij If its value is positive and greater than the threshold, the generation rule condition S is... i With S j Positive correlation enhancement; if negative and the absolute value is greater than the threshold, generate S.i With S j Negative correlation enhancement. Connect all significant conditions using a logical AND operation to form the rule antecedents.

[0133] Furthermore, this embodiment provides specific calculation logic for the process of associating the risk level corresponding to the cluster center as a rule consequent. This includes the following sub-steps:

[0134] The percentage of samples labeled as faulty within the sample set belonging to each cluster center is calculated. Historical data typically includes maintenance records or fault labels (e.g., normal, bearing overheating, blade cavitation). For the k-th cluster, the percentage of samples labeled as abnormal within that cluster is calculated, denoted as P. k .

[0135] The failure rate is compared with preset classification thresholds. In this embodiment, two key thresholds are preset: a high-risk threshold (e.g., 0.3) and a medium-risk threshold (e.g., 0.1).

[0136] If the failure rate is higher than the high-risk threshold, the rule corresponding to the cluster center is associated with a Level 1 warning; if the failure rate is between the high-risk and medium-risk thresholds, it is associated with a Level 2 warning. If the failure rate is lower than the medium-risk threshold, the cluster represents a normal or low-risk operating condition, and alarm rules are usually not generated or only a prompt message is generated. Through a probability-based grading method, an automatic mapping from unsupervised clustering to supervised risk management is achieved.

[0137] In one embodiment of the present invention, when multiple rules simultaneously meet the triggering conditions, rule conflict detection and resolution are required.

[0138] The rule conflict detection method is as follows: for any two rules R a and R b If the condition sets of the two rules intersect and the warning actions are inconsistent, they are considered a potentially conflicting rule pair. The formula for calculating the degree of conflict is:

[0139] Conflict(R a ,R b =|Conds a ∩Conds b | / min(|Conds a |,|Conds b |)*I(Action a !=Action b );

[0140] Where |.| represents the number of elements in the set, ∩ represents the intersection operation, and I(.) is an indicator function that takes the value 1 when the condition is true and 0 otherwise.

[0141] The following strategies are used to resolve rule conflicts:

[0142] Strategy 1: Priority resolution. When Conflict(R) a ,R b When the score is greater than 0.5 and both rules are triggered simultaneously, only the rule with the higher priority score will be executed.

[0143] Execute(R a If Priority(R) a Priority(R) b );

[0144] Execute(R b If Priority(R) b Priority(R) a ).

[0145] Strategy Two: Merge and Resolve. When the warning actions of conflicting rules are compatible, generate a comprehensive warning action. For example, when R... a It is recommended to reduce power and R b When checking for vibration, it is recommended to combine two suggestions and perform them simultaneously.

[0146] Strategy 3: Conservative Resolution. When priorities cannot be determined or merging is not possible, a conservative strategy is adopted, which means selecting the rule with the higher warning level for execution.

[0147] The quality assessment of rules includes the following metrics:

[0148] Indicator 1: Rule accuracy.

[0149] Precision(R) = TP / (TP + FP + ε);

[0150] Where TP is the number of true positives (the rule is triggered and the exception actually occurs), and FP is the number of false positives (the rule is triggered but no exception occurs).

[0151] Indicator 2: Rule recall rate.

[0152] Recall(R) = TP / (TP + FN + ε);

[0153] FN represents the number of false negatives (an anomaly occurs but the rule is not triggered).

[0154] Indicator 3: Rule F1 score.

[0155] F1(R)=2*Precision*Recall / (Precision+Recall+ε);

[0156] Indicator 4: Rule stability.

[0157] Stability(R) = 1 - Std(Match scores ) / Mean(Match scores );

[0158] Among them, Match scores This represents the sequence of matching scores for the rule across historical triggering scenarios, where Std represents the standard deviation and Mean represents the mean. Higher stability indicates that the rule performs more consistently across different scenarios.

[0159] Indicator 5: Rule interpretability score.

[0160] Interpretability(R)=1 / (1+log(|Conds|));

[0161] The fewer conditions a rule has, the more interpretable it is.

[0162] The formula for calculating the overall quality score is:

[0163] Quality(R)=0.3*F1+0.2*Stability+0.2*Interpretability+0.3*Score recency ;

[0164] When the overall quality score of the rule is lower than the threshold Q min When the default value is 0.3, the rule will be marked as pending optimization or deletion.

[0165] Example 6 describes an engineering optimization strategy for the rule matching and execution phases, which solves the practical problems of boundary jitter and multi-rule conflicts.

[0166] Step 601: For each rule in the sensitivity pattern rule base, calculate the numerical deviation between the element value in the multi-factor coupling sensitivity matrix and the threshold set by the rule.

[0167] During real-time operation, for each rule r in the rule base, the system extracts all salient feature dimensions defined in its antecedents. Assuming rule r requires vibration sensitivity > 0.3, and the currently calculated vibration sensitivity value is s, then the calculated numerical deviation x = s - 0.3.

[0168] Step 602: Using a nonlinear mapping function that includes a temperature parameter, the numerical deviation is mapped to a soft matching degree with a value between 0 and 1; the temperature parameter is used to control the fuzzy boundary range of the matching determination.

[0169] To avoid frequent start-stop oscillations caused by hard thresholds (i.e., alarming when the threshold is greater than 0.3, but not alarming when the threshold is 0.29), this embodiment employs a soft matching mechanism. The nonlinear mapping function is preferably a Sigmoid function with a temperature parameter. The specific formula can be expressed as: μ = 1 / 1 + exp(-(s-θ) / τ); where μ is the soft matching degree for a certain rule condition (ranging from 0 to 1); s is the relative sensitivity value calculated at the current moment; θ is the preset sensitivity trigger threshold in the rule; s-θ is the sensitivity numerical deviation; and τ is the temperature parameter, used to adjust the steepness of the function. The smaller the τ value, the closer the function is to a step function; the larger the τ value, the more blurred the boundary. In a preferred embodiment, the temperature parameter τ can be set to 0.1 times the absolute value of the rule threshold. Through this calculation, any numerical deviation can be converted into a probability value between 0 and 1, i.e., the soft matching degree.

[0170] Step 603: When the soft matching degree is greater than the preset confidence threshold, it is determined that the current running state belongs to the sensitivity mode corresponding to the rule.

[0171] The system sets a global confidence threshold (e.g., 0.7). A rule is considered activated only when its overall soft match score (usually the weighted geometric mean of the soft match scores of each condition) exceeds this threshold.

[0172] Step 604: When multiple rules belonging to different sensitivity modes are found during matching, calculate the priority score of each matching rule.

[0173] In actual operation, multiple rules may be triggered simultaneously (for example, the conditions for both Level 1 and Level 2 warnings may be met at the same time). In this case, arbitration is required by calculating priority scores.

[0174] Step 605, the priority score is calculated based on the following weighted indicators: the urgency of the warning defined by the rule, the matching confidence of the rule, and the scenario similarity between the current real-time multi-source operation monitoring data and the pre-stored historical trigger scenarios of the rule.

[0175] In this embodiment, the priority score P(r) is calculated as follows: P(r) = w1 × U(r) + w2 × M(r) × conf(r) + w3 × Sim(r, xt); where P(r) is the comprehensive priority score of rule r; w1, w2, w3 are the weighting coefficients of each evaluation index (e.g., w1 = 0.3, w2 = 0.3, w3 = 0.4); U(r) is the urgency of the warning, which is determined by the warning level of the rule's consequent (e.g., 1.0 for red warning and 0.5 for yellow warning); M(r) is the currently calculated soft matching degree; conf(r) is the static confidence degree of the rule itself; and Sim(r, xt) is the scene similarity.

[0176] Specifically, scenario similarity measures the distance between the current real-time monitoring data vector xt and the typical state vector (historical trigger scenario) when the rule was in effect historically. In the calculation, normalized Euclidean distance is used, and a Gaussian kernel function is applied to convert it into a similarity value between 0 and 1. This metric is introduced because if the current operating condition is highly similar to a rule's historical effective scenario, it indicates that the rule has a higher relevance to the current state and should be given higher priority.

[0177] Step 606: Select the rule with the highest priority score and generate an early warning signal based on that rule, such as... Figure 3 As shown.

[0178] The system prioritizes all triggering rules and executes only the warning action corresponding to the rule with the highest score, or executes the top-N rules with the highest scores, thus eliminating semantic conflicts and ensuring the uniqueness and accuracy of control decisions.

[0179] Example 7 describes a system that corrects characteristic changes caused by equipment degradation, achieving adaptive optimization throughout the entire lifecycle.

[0180] Step 701: Establish a sliding time window for each rule in the sensitivity mode rule base, and calculate the trigger accuracy of the rule within the sliding time window based on the subsequent actual running results.

[0181] The system maintains a sliding counter for each online rule, recording its performance over a recent period (e.g., the last 7 days or the last 100 triggers). Statistical metrics include the number of triggers and the number of hits (i.e., the number of times a trigger was manually confirmed or post-event analysis confirmed that a corresponding fault actually occurred). The ratio of these two is the trigger accuracy rate.

[0182] Step 702: Construct a statistical hypothesis test statistic to measure whether the trigger accuracy has decreased relative to the historical baseline accuracy of the rule.

[0183] Historical baseline accuracy refers to the accuracy of the rule during its stable period immediately after implementation or the last update. To scientifically determine whether a decrease in accuracy is caused by random fluctuations, this embodiment employs statistical hypothesis testing methods. Specifically, a Z-test statistic or an exact binomial test statistic is constructed. Taking the Z-test as an example, the calculation of the statistic Z involves the difference between the current accuracy and the baseline accuracy, divided by the standard error.

[0184] Step 703: When the statistical hypothesis test statistic exceeds the significance level threshold, the rule is deemed invalid.

[0185] In this embodiment, a significance level threshold Z needs to be set. α (For example, 1.645 corresponding to a 95% confidence level). If the calculated Z-value is less than a negative Z... α If the null hypothesis is rejected, it is considered that the accuracy has decreased in a statistically significant way, and the judgment rule is invalid.

[0186] Furthermore, this embodiment introduces physical constraints into the threshold setting. Specifically, the significance level threshold is dynamically adjusted based on the equipment degradation model. The equipment degradation model defines the functional relationship between the equipment degradation factor and the equipment service time. Preferably, an exponential decay model is adopted: D(T) = 1 - exp(-λ × T); where D(T) is the degradation factor of the equipment at service time T (ranging from 0 to 1); T is the cumulative service time of the equipment; and λ is the degradation rate constant, estimated from historical maintenance records. As the equipment service time increases, the significance level threshold is reduced using the equipment degradation model, thereby increasing the sensitivity to rule failure determination.

[0187] Adjusted threshold Z adj =Z α ×(1-γ×D(T)); where Z adj Z is the significance level threshold after equipment degradation correction (the critical value for determining rule failure); α Z represents the initial significance level threshold set based on statistical standards (e.g., 1.645); γ is the adjustment coefficient used to control the weight of the degradation factor on the threshold; D(T) is the degradation factor of the equipment at service time T (ranging from 0 to 1). The older the equipment, the larger D(T) is, and the higher the threshold Z is. adj The smaller the value (the smaller the absolute value), the more sensitive the system is to the decline in rule accuracy, prompting the rules to be updated more frequently to adapt to the rapidly changing operating characteristics of older equipment.

[0188] Step 704: Collect historical sensitivity samples that trigger the rule within the sliding time window, perform local re-clustering, and generate new rules based on the new cluster centers to replace the invalid rules.

[0189] When a rule is deemed invalid, the system does not need to retrain the entire model or reconstruct the entire rule base. Instead, it only extracts the sensitivity samples captured by that rule during the sliding window. These local samples are then re-clustered and feature extracted to generate new rule parameters (such as new thresholds and new coupling relationships), which are then used to replace the invalid old rules in situ. This local update mechanism reduces computational costs and enables the system's online adaptive evolution.

[0190] Example 8: Description of the hardware entity architecture for performing the above method.

[0191] This embodiment provides an electronic device, which can specifically be an edge computing server for a pumping station, a local control unit (LCU), or a cloud-based operation and maintenance server. The electronic device includes a processor, memory, communication interface, and display unit connected via a system bus.

[0192] The processor can be a central processing unit (CPU), a graphics processing unit (GPU), a field-programmable gate array (FPGA), or an application-specific integrated circuit (ASIC). Especially when performing deep learning predictive model inference and calculating second derivatives based on automatic differentiation, a GPU with parallel computing capabilities is preferred to meet millisecond-level real-time requirements. The processor is configured to execute computer program instructions stored in memory to implement the pump station unit status trend prediction and early warning method.

[0193] The memory includes non-volatile storage media and internal memory. The non-volatile storage media (such as hard disks or SSDs) stores the operating system, database management system, and program code for the pump station unit status trend prediction and early warning system of this embodiment. Specifically, the physical causal order, sensitivity pattern rule base, and equipment degradation model parameters are persistently stored in this non-volatile storage media in the form of files or database tables. Internal memory (RAM) provides a high-speed cache environment for the processor's operations and is used to temporarily store intermediate variables such as real-time multi-source operation monitoring data and multi-factor coupled sensitivity matrices.

[0194] The communication interface is used to enable data exchange between electronic equipment and field sensors (such as vibration sensors and water level gauges) at the pump station, as well as data acquisition and monitoring control systems (SCADA) or remote monitoring terminals. For example, real-time multi-source operation monitoring data can be read in via this communication interface using TCP / IP or Modbus protocols; generated early warning signals can also be sent to the central control room's large screen or the mobile terminals of maintenance personnel via this interface.

[0195] Furthermore, this embodiment also provides a computer-readable storage medium on which a computer program is stored. When the computer program runs on a computer or processor, it causes the computer to execute a method for predicting and warning the status trend of pump station units. The storage medium includes, but is not limited to, a USB flash drive, a portable hard drive, an optical disc, a read-only memory (ROM), and a random access memory (RAM).

[0196] To address the black-box nature and lack of interpretability of deep learning models, this solution transforms the obscure model prediction logic into a calculable causal contribution by calculating a multi-factor coupling sensitivity matrix, thus intuitively demonstrating the specific impact weight of each monitoring indicator on the unit's state trend.

[0197] To address the shortcomings of existing analytical methods that neglect physical causal chains and fail to accurately isolate multivariate coupling effects, this solution employs an integral gradient algorithm based on physical causal constraints. By constructing an integral path that strictly follows the principle that upstream variables change before downstream variables, and further calculating the two-factor coupling sensitivity including direct interaction terms and causal transmission terms, the complex hydraulic-mechanical-electromagnetic multiphysics interaction mechanism of the pumping station is successfully decoupled at the numerical computation level, thus reconstructing the actual fault propagation path.

[0198] To address the problem that continuous predicted values ​​are difficult to convert into discrete decision rules, this solution automatically extracts and constructs a sensitivity pattern rule base containing significant causal features by performing unsupervised clustering analysis on historical sensitivity samples. This achieves a leap from simple numerical alarms to logical diagnosis with causal analysis capabilities.

[0199] In this application, the coupling strength between variables is calculated to provide a basis for the automatic generation of rule conditions; sensitivity information is transformed into interpretable early warning rules; and a dynamic update mechanism is used to adapt to the long-term evolution of the system.

[0200] Sensitivity analysis alone can only provide a ranking of feature importance, but when combined with rule generation, it can be directly used for early warning decisions. Rule learning alone may generate a large number of redundant or low-quality rules, but introducing sensitivity priors improves the relevance and accuracy of the rules. Drift detection alone requires users to manually process failed rules, but when linked with the rule generation module, it can achieve automatic rule evolution.

[0201] 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 method for predicting and warning the status trend of pump station units, characterized in that, include: Real-time multi-source operation monitoring data of the pump station unit is acquired, and the real-time multi-source operation monitoring data is input into a pre-trained deep learning trend prediction model to obtain the unit status trend prediction result. Based on the unit status trend prediction results and real-time multi-source operation monitoring data, a multi-factor coupling sensitivity matrix is ​​calculated. The multi-factor coupled sensitivity matrix is ​​input into a preset sensitivity mode rule base for matching and retrieval to determine the sensitivity mode to which the current operating state belongs. When the warning triggering conditions are met in the sensitivity mode, a warning signal is generated according to the matched rules; The steps for calculating the multi-factor coupling sensitivity matrix specifically include: Obtain the pre-defined physical causal order among various monitoring indicators of the pumping station; Determine the baseline state of real-time multi-source operation monitoring data, and construct a causal constraint integral path connecting the baseline state and real-time multi-source operation monitoring data according to the physical causal order; Gradient integral calculation is performed on the deep learning trend prediction model along the causal constraint integral path to obtain the multi-factor coupling sensitivity matrix; The causal constraint integral path is defined by a set of interpolation variables that satisfy causal monotonicity constraints; whereby the causal monotonicity constraints refer to: In the physical causal sequence, the rate of change of the interpolated variable corresponding to the upstream monitoring indicator is always greater than or equal to the rate of change of the interpolated variable corresponding to the downstream monitoring indicator. The off-diagonal elements of the multi-factor coupling sensitivity matrix, representing the degree of two-factor coupling interaction between two different monitoring data, consist of direct interaction terms and causal transmission terms. The direct interaction term characterizes the integral contribution of the second-order mixed partial derivative calculated along the causal constraint integral path, and is used to calculate the nonlinear synergistic effect generated by the simultaneous change of two monitoring data. The causal transit term characterizes the indirect impact contribution of upstream monitoring data to downstream monitoring data through intermediate variables in the causal chain.

2. The method according to claim 1, characterized in that, The multi-factor coupling sensitivity matrix is ​​a symmetric matrix, where: The diagonal elements of the multi-factor coupling sensitivity matrix represent the degree of independent influence of a single monitoring data point on the prediction results of the unit's state trend. The off-diagonal elements of the multi-factor coupling sensitivity matrix characterize the degree of two-factor coupling interaction between two different monitoring data.

3. The method according to claim 1, characterized in that, The multi-factor coupled sensitivity matrix is ​​input into a preset sensitivity pattern rule base for matching and retrieval, specifically including: For each rule in the sensitivity pattern rule base, calculate the numerical deviation between the element values ​​in the multi-factor coupled sensitivity matrix and the threshold set by that rule; A nonlinear mapping function incorporating temperature parameters is used to map numerical deviations into soft matching degrees ranging from 0 to 1; the temperature parameter is used to control the fuzzy boundary range of the matching determination. When the soft matching degree is greater than the preset confidence threshold, the current running state is determined to belong to the sensitivity mode corresponding to the rule.

4. The method according to claim 1, characterized in that, Based on the matched rules, an early warning signal is generated, specifically including: When multiple rules belonging to different sensitivity modes are found during the matching retrieval, the priority score of each matching rule is calculated. The priority score is calculated based on the following weighted indicators: the urgency of the warning defined by the rule, the matching confidence of the rule, and the scenario similarity between the current real-time multi-source operation monitoring data and the pre-stored historical trigger scenarios of the rule; Select the rule with the highest priority score and generate an early warning signal based on that rule.

5. The method according to claim 1, characterized in that, The deep learning trend prediction model consists of a convolutional neural network layer, a bidirectional long short-term memory network layer, and a temporal attention mechanism layer connected in sequence, specifically including: Convolutional neural network layers are used to extract local correlation features between various monitoring indicators in real-time multi-source operation monitoring data; Bidirectional long short-term memory network layers are used to extract forward and backward temporal dependency information of monitoring data based on local correlation features; The temporal attention mechanism layer is used to assign weights to temporal dependency information at different historical moments, thereby enhancing the contribution of key historical information to the prediction results of unit status trends.

6. The method according to claim 1, characterized in that, The preset sensitivity pattern rule base is pre-built through the following steps: Collect historical operating data of the pumping station and calculate the multi-factor coupled sensitivity matrix corresponding to each historical moment, and flatten it into a historical sensitivity sample vector; Cluster analysis was performed on the historical sensitivity sample vectors to identify several cluster centers that represent typical operating states; Extract significant sensitivity features from each cluster center. Significant sensitivity features include single-factor sensitivity components or two-factor coupled sensitivity components whose values ​​exceed the threshold. The significant sensitivity features are combined and transformed into corresponding rule antecedent logic, and the risk level corresponding to the cluster center is associated as the rule consequent to generate a rule base.

7. The method according to claim 1, characterized in that, The method also includes the step of dynamically updating the sensitivity pattern rule base: A sliding time window is established for each rule in the sensitivity mode rule base, and the trigger accuracy of the rule within the sliding time window is statistically analyzed based on the subsequent actual running results. Construct a statistical hypothesis test statistic to measure whether the trigger accuracy has decreased relative to the historical baseline accuracy of the rule; The rule is deemed invalid when the statistical hypothesis test statistic exceeds the significance level threshold. Collect historical sensitivity samples that trigger the rule within the sliding time window, perform local re-clustering, and generate new rules based on the new cluster centers to replace the invalid rules.

Images