Online testing system for full life cycle performance decline of reciprocating pump

Abstract

The invention discloses an online testing system for full life cycle performance decline of a reciprocating pump, and particularly relates to the technical field of reciprocating pump state monitoring. A dynamic feature matching module identifies working conditions and outputs a fault feature set; the weak signal enhancement module enhances and extracts nonlinear weak fault features from the original signal and outputs an enhanced feature signal; the coupling and decoupling analysis module performs multi-source decoupling on the enhanced signal and outputs an independent fault characteristic value sequence and a fault propagation path; the full life cycle mapping module generates a health state quantitative index according to the sequence and the path; the reference alignment module carries out dynamic space-time reference alignment and outputs an aligned track geometric data set; the dynamic feature matching module adapts to working condition changes, early weak faults are extracted through the weak signal enhancement module, multi-cylinder interference is separated through the coupling decoupling and mapping module, the health state is quantified, and online testing of full-life periodic performance decline of the reciprocating pump is achieved.

Description

Technical Field

[0001] This invention relates to the field of reciprocating pump condition monitoring technology, and more specifically, to an online testing system for the performance degradation of a reciprocating pump throughout its entire life cycle. Background Technology

[0002] A reciprocating pump is a positive displacement pump that uses the reciprocating motion of a piston or plunger to transport fluids. It is widely used in industrial fields such as oil drilling, chemical processes, mine drainage, and power transmission. During long-term operation, key components of a reciprocating pump, such as pump valves, plungers, seals, and bearings, will gradually experience wear, fatigue, corrosion, and other performance degradation. Timely detection of these degradation conditions is crucial for ensuring safe equipment operation and reducing maintenance costs.

[0003] Currently, the performance monitoring of reciprocating pumps mainly adopts a combination of regular inspections and online parameter threshold alarms. Common online monitoring methods include vibration monitoring, temperature monitoring, pressure pulsation analysis, and flow monitoring. These methods determine whether the equipment is abnormal by setting fixed thresholds or simple statistical indicators (such as vibration amplitude and kurtosis factor). Some advanced testing systems introduce time-domain statistical features, frequency-domain features (such as sideband amplitude), and machine learning classification methods (such as support vector machines and neural networks) for fault identification. These methods can achieve good diagnostic results under steady-state conditions.

[0004] In actual industrial settings, the operating conditions of reciprocating pumps are often not constant. Factors such as speed regulation, load changes, and alterations in the characteristics of the transported medium lead to strong time-varying characteristics in the operating conditions. Most existing monitoring methods rely on steady-state assumptions or diagnostic models trained for specific operating conditions, resulting in fluctuations in detection accuracy as operating conditions change. Furthermore, fault signals generated in the early stages of performance degradation have low amplitudes and are easily affected by vibrations from normal equipment operation and environmental noise. Conventional signal processing methods struggle to reliably extract these weak features. For multi-cylinder reciprocating pumps, vibration signals between cylinders are coupled through mechanical structures, and the signals collected by sensors are typically a superposition of vibrations from multiple cylinders. This adds further difficulty to accurately locating the fault source and tracing the fault propagation path between cylinders. Therefore, this invention proposes an online testing system for the entire lifecycle performance degradation of reciprocating pumps to address the aforementioned problems. Summary of the Invention

[0005] To achieve the above objectives, the present invention provides the following technical solution: An online testing system for the performance degradation of a reciprocating pump throughout its entire life cycle includes: The dynamic feature matching module is used to identify changes in the operating conditions of reciprocating pumps in real time and dynamically output the corresponding fault feature set. The weak signal enhancement module is used to receive the fault feature set, enhance and extract nonlinear weak fault features from the original sensing signal, and output the enhanced feature signal. The coupling and decoupling analysis module is used to perform multi-source information coupling and decoupling analysis on the enhanced feature signal, and outputs the decoupled independent fault feature value sequence and fault propagation path; The full lifecycle mapping module is used to generate quantitative indicators of the current health status of reciprocating pumps based on the independent fault characteristic value sequence and fault propagation path mapping.

[0006] In a preferred embodiment, the operating conditions of the reciprocating pump obtained by the dynamic feature matching module include at least the rotational speed, output pressure, instantaneous flow rate, and characteristics of the conveying medium; the characteristics of the conveying medium include medium density, viscosity, and sand content.

[0007] The specific steps of the dynamic feature matching module to identify changes in the operating conditions of the reciprocating pump in real time and dynamically output the corresponding fault feature set are as follows: The real-time speed, output pressure, instantaneous flow rate, and medium density, viscosity, and sand content of the reciprocating pump are continuously collected by speed sensors, pressure sensors, flow sensors, and medium characteristic sensors to form a sequence of operating parameters. A sliding time window is used to determine the steady state of the operating condition parameter sequence. The coefficient of variation of each parameter within the time window is calculated. When the coefficients of variation of all parameters are lower than the corresponding preset variation threshold, it is determined to be a steady state operating condition; otherwise, it is determined to be a dynamic operating condition. Based on the steady-state or dynamic discrimination results, the pre-trained working condition classification model is invoked to output the working condition category label; Based on the operating condition category label, the corresponding fault feature set is retrieved from the preset fault feature knowledge base and dynamically output; the fault feature set includes the time domain feature set and the wavelet packet energy feature set.

[0008] In a preferred embodiment, the operating condition classification model maps the current operating condition parameter vector to the operating condition category cluster based on a self-organizing map neural network.

[0009] In a preferred embodiment, the weak signal enhancement module receives the fault feature set in the following manner: The fault feature set is analyzed to extract the sensitive frequency band and sensitive time scale parameters related to weak faults; among them, the original sensing signal is the vibration acceleration signal of multiple measuring points of each cylinder of the reciprocating pump collected by the vibration sensor. The sensitive frequency band is extracted as follows: taking the vibration acceleration signal of the reciprocating pump when it is running in a healthy state as a reference, the original sensing signal is decomposed into three layers of wavelet packets to obtain eight frequency band subspaces. The rate of change of wavelet packet energy of each subspace relative to the healthy reference is calculated, and the frequency ranges corresponding to the first three subspaces with kurtosis index greater than three are determined as the sensitive frequency bands. The sensitive time scale parameter is extracted as follows: the wavelet packet decomposition level is determined based on the ratio of the lowest frequency in the sensitive frequency band to the sampling frequency, and the time resolution corresponding to the level is used as the sensitive time scale parameter.

[0010] In a preferred embodiment, the specific steps of the weak signal enhancement module in enhancing and extracting nonlinear weak fault features from the original sensing signal and outputting the enhanced feature signal are as follows: Based on the sensitive time scale parameter, the original sensing signal is segmented according to the time scale to obtain multiple signal segments; Each signal segment is normalized by its maximum and minimum values ​​so that all signal amplitudes are between negative one and positive one, thus obtaining a normalized signal segment. Each normalized signal segment is input into a bistable stochastic resonance system. The system parameters of the bistable stochastic resonance system include potential well parameters and noise intensity parameters. The signal-to-noise ratio gain maximization criterion is used as the objective function for parameter optimization. The center frequency of the sensitive frequency band is used as the characteristic frequency of the system's expected response. The search range of the potential well parameters is set from a minimum value of 0.1 to a maximum value of 10, and the search range of the noise intensity parameters is set from a minimum value of 0.01 to a maximum value of 100. The optimal potential well parameters and noise intensity parameters are determined through iterative search, so that the system is in a stochastic resonance state. The output signals of each normalized signal segment after random resonance processing are spliced ​​together in the original time sequence to obtain the random resonance enhanced signal. The random resonance enhancement signal is subjected to three-layer wavelet soft thresholding denoising to eliminate residual noise and obtain the final enhanced feature signal.

[0011] The enhanced feature signals received by the coupling and decoupling analysis module are the enhanced feature signals corresponding to multiple measuring points in each cylinder of the reciprocating pump; Multi-source information coupling and decoupling analysis refers to separating the coupling components in the enhanced characteristic signals of each cylinder block measuring point, eliminating the mutual interference between vibration signals of adjacent cylinder blocks, and outputting a decoupled independent fault characteristic value sequence corresponding to each cylinder block.

[0012] In a preferred embodiment, the specific steps for the coupling-decoupling analysis module to output the decoupled independent fault feature value sequence and fault propagation path are as follows: The enhanced characteristic signals of each cylinder block measuring point are constructed into an observation signal matrix. Using the independent component analysis algorithm with kurtosis as the comparison function, the unmixing matrix is ​​calculated iteratively to decompose the observation signal matrix into mutually statistically independent source signal components. Each source signal component corresponds to the vibration contribution of a cylinder block, thus obtaining the independent fault characteristic value sequence of each cylinder block. Using the obtained independent fault feature value sequence of each cylinder as input, the cross-correlation function of the independent fault feature value sequence between adjacent cylinders is calculated in time order. When the peak value of the cross-correlation function exceeds the preset cross-correlation threshold and there is a time delay, the propagation direction is determined according to the sign of the time delay corresponding to the peak value of the cross-correlation function: if the delay is positive, it is determined that the fault propagates from the previous cylinder to the next cylinder; if the delay is negative, it is determined that the fault propagates from the next cylinder to the previous cylinder. The propagation directions are connected in sequence to form the fault propagation path. The fault propagation path is output in the form of cylinder block number sequence, indicating the initial fault source cylinder block and the subsequent cylinder block sequence affected by the propagation.

[0013] In a preferred embodiment, the specific steps of the full lifecycle mapping module in generating quantitative indicators of the current health status of the reciprocating pump based on the independent fault feature value sequence and fault propagation path mapping are as follows: For each cylinder, its current independent fault feature value sequence is compared with the health baseline feature value sequence, and the health score of each cylinder is calculated. The health baseline feature value sequence is the independent fault feature value sequence of each cylinder of the reciprocating pump in a healthy state. The health score is equal to the ratio of the Euclidean distance between the current independent fault feature value sequence and the health baseline feature value sequence to the sum of the modulus of the health baseline feature value sequence and the modulus of the current independent fault feature value sequence, and the value is zero when the calculation result is less than zero. Based on the fault propagation path, determine the sequential position of each cylinder in the propagation path, and correct the health score of the cylinder that is later in the propagation sequence by multiplying it by the propagation attenuation coefficient. The propagation attenuation coefficient ranges from 0.8 to 1, and the coefficient decreases as the propagation distance increases. The health scores of all cylinders after correction are weighted and averaged. The weight of each cylinder is pre-assigned according to its propagation order in the fault propagation path. The cylinders that are earlier in the propagation order have lower weights. This yields a quantitative health status index between zero and one, which represents the current degree of degradation of the reciprocating pump throughout its entire life cycle.

[0014] The technical effects and advantages of this invention are as follows: This invention identifies changes in the operating conditions of a reciprocating pump in real time through a dynamic feature matching module and dynamically outputs the corresponding fault feature set. This enables the system to adaptively adjust the types of fault features of interest according to changes in speed, pressure, flow rate and medium characteristics, avoiding misjudgment or omission of fixed feature sets under changing operating conditions. This significantly improves the dynamic adaptability and operating condition robustness of performance degradation testing throughout the entire life cycle.

[0015] This invention enhances and extracts nonlinear weak fault features from the original sensing signal through a weak signal enhancement module and outputs enhanced feature signals. Combined with the extraction of sensitive frequency band and sensitive time scale parameters, it can effectively enhance weak fault signals such as early wear and micro-leakage from strong background noise, solving the problem that early performance degradation signals are difficult to capture by conventional methods, and realizing early online identification of hidden faults in reciprocating pumps.

[0016] This invention performs multi-source information coupling and decoupling analysis on enhanced feature signals through a coupling and decoupling analysis module, outputs the decoupled independent fault feature value sequence and fault propagation path, and uses a full life cycle mapping module to generate quantitative health status indicators. It can eliminate mutual interference of multi-cylinder vibration, trace the propagation sequence of faults between cylinders, and provide maintenance personnel with quantitative values ​​of the overall pump health, thereby supporting predictive maintenance decisions and extending the effective service life of equipment. Attached Figure Description

[0017] To facilitate understanding by those skilled in the art, the present invention will be further described below with reference to the accompanying drawings; Figure 1 This is a schematic diagram of an online testing system for the performance degradation of a reciprocating pump throughout its entire life cycle, as described in this invention. Detailed Implementation

[0018] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative effort are within the scope of protection of the present invention.

[0019] Reference Figure 1 The following examples were obtained: Example 1: An online testing system for the performance degradation of a reciprocating pump throughout its entire life cycle, comprising: The dynamic feature matching module is used to identify changes in the operating conditions of reciprocating pumps in real time and dynamically output the corresponding fault feature set. The weak signal enhancement module is used to receive the fault feature set, enhance and extract nonlinear weak fault features from the original sensing signal, and output the enhanced feature signal. The coupling and decoupling analysis module is used to perform multi-source information coupling and decoupling analysis on the enhanced feature signal, and outputs the decoupled independent fault feature value sequence and fault propagation path; The full lifecycle mapping module is used to generate quantitative indicators of the current health status of reciprocating pumps based on the independent fault characteristic value sequence and fault propagation path mapping.

[0020] The dynamic feature matching module acquires the operating conditions of the reciprocating pump, including at least rotational speed, output pressure, instantaneous flow rate, and characteristics of the transported medium. The transported medium characteristics include density, viscosity, and sand content. A speed sensor is installed at the crankshaft end of the reciprocating pump to measure the number of crankshaft rotations per minute (rpm). A pressure sensor is installed in the pump outlet pipe to measure the fluid pressure at the pump output end (MPa). A flow sensor is installed in the outlet pipe to measure the volume of fluid discharged by the pump per unit time (m³ / h). Medium characteristic sensors include a densitometer, a viscometer, and a sand content detector, used to measure the density (kg / m³), dynamic viscosity (Pascal-second), and volume fraction of solid particles (percentage), respectively. These sensors continuously collect data at a sampling frequency of ten times per second, and the data is arranged in chronological order to form a sequence of operating condition parameters. For example, during normal drilling operations, a three-cylinder reciprocating pump can collect the following data within 30 seconds: rotational speed series fluctuating between 120 and 135 revolutions per minute; pressure series fluctuating between 15 and 18 MPa; flow rate series fluctuating between 30 and 33 cubic meters per hour; medium density series fluctuating between 1,200 and 1,250 kilograms per cubic meter; viscosity series fluctuating between 0.05 and 0.06 Pascals per second; and sand content series fluctuating between 3% and 5%.

[0021] The specific steps of the dynamic feature matching module to identify the operating condition changes of the reciprocating pump in real time and dynamically output the corresponding fault feature set are as follows: continuously collect the real-time speed, output pressure, instantaneous flow rate, and medium density, viscosity, and sand content of the reciprocating pump through speed sensor, pressure sensor, flow sensor and medium characteristic sensor to form a sequence of operating condition parameters. A sliding time window is used to determine the steady-state condition of the operating condition parameter sequence. The coefficient of variation (COP) of each parameter within the time window is calculated. A steady-state condition is defined as when the COPs of all parameters are below their corresponding preset COP thresholds; otherwise, a dynamic condition is defined. For example, the sliding time window length is set to ten seconds, meaning each time window contains one hundred sampling points. The time window slides forward in steps of one step per second. For each operating condition parameter within each time window, the arithmetic mean of the parameter within that time window is first calculated, followed by the standard deviation. The COP is then obtained by dividing the standard deviation by the arithmetic mean. The COP is dimensionless.

[0022] Each operating condition parameter corresponds to a preset variation threshold. These thresholds are set based on the following: the speed variation threshold is set to 5%, because the speed fluctuation of a reciprocating pump during normal operation usually does not exceed 5% of the rated speed; the output pressure variation threshold is set to 8%, because the outlet pressure of a reciprocating pump is affected by pipeline impedance and medium characteristics, and the allowable fluctuation range is slightly greater than that of the speed; the instantaneous flow variation threshold is set to 8%, corresponding to the pressure fluctuation range; the medium density variation threshold is set to 3%, because the drilling fluid density is relatively stable during circulation; the viscosity variation threshold is set to 5%, because viscosity changes relatively slowly; and the sand content variation threshold is set to 10%, because sand content fluctuates greatly due to the influence of formation sand production, so the threshold setting is relatively lenient.

[0023] When the coefficients of variation of all six parameters within the current time window are lower than their respective preset thresholds, the condition is considered steady-state. If the coefficient of variation of any single parameter exceeds its corresponding threshold, the condition is considered dynamic. For example, in the aforementioned 30 seconds of drilling data, the coefficients of variation for the rotational speed (4%), pressure (6%), flow rate (7%), density (2%), viscosity (3%), and sand content (8%) were all lower than their corresponding thresholds in the first 10 seconds, thus indicating a steady-state condition. Between the 20th and 30th seconds, the sand content jumped from 3% to 8%, with its coefficient of variation reaching 15%, exceeding the 10% threshold, therefore indicating a dynamic condition.

[0024] Based on the steady-state or dynamic discrimination result, the pre-trained operating condition classification model is invoked to output operating condition category labels. When the discrimination result is a steady-state operating condition, the operating condition classification model outputs the corresponding operating condition category label based on the average operating condition parameter vector within the current time window. When the discrimination result is a dynamic operating condition, the operating condition classification model outputs the corresponding transition operating condition category label based on the operating condition parameter change trajectory within a sliding time window. The operating condition category labels are predefined, such as "low speed steady state," "high speed steady state," "pressure jump transition," and "sudden increase in sand content transition."

[0025] Based on the operating condition category label, the system retrieves and dynamically outputs the corresponding fault feature set from a pre-set fault feature knowledge base. The fault feature set includes a time-domain feature set and a wavelet packet energy feature set. The pre-set fault feature knowledge base is a data structure that stores the mapping relationship between different operating condition category labels and their corresponding fault feature sets. Each operating condition category label corresponds to a set of fault feature sets specifically designed for that operating condition.

[0026] The time-domain feature set in the fault feature set refers to the set of dimensionless statistical indices directly calculated from the original time series of the vibration acceleration signal, including but not limited to kurtosis, peak factor, impulse factor, margin factor, and waveform factor. The kurtosis reflects the intensity of the impact component in the signal; the peak factor reflects the ratio of the signal's peak value to its effective value; the impulse factor reflects the ratio of the signal's peak value to its rectified average value; the margin factor reflects the ratio of the signal's peak value to its square root amplitude; and the waveform factor reflects the ratio of the signal's effective value to its rectified average value. These time-domain features are sensitive to mechanical faults in reciprocating pumps, such as pump valve impact and bearing wear. The wavelet packet energy feature set refers to the distribution vector of signal energy in each frequency band subspace after wavelet packet decomposition of the vibration acceleration signal. Specifically, the signal is decomposed into eight frequency band subspaces through three levels of wavelet packet decomposition. The sum of squares of the signal in each subspace is calculated as the energy of that subspace. Then, the eight energy values ​​are normalized (i.e., the energy of each subspace is divided by the total energy of the eight subspaces), resulting in a wavelet packet energy feature set composed of eight normalized energy values. Wavelet packet energy characteristics can reflect the changes in energy distribution of fault signals at different frequency scales and are sensitive to early, weak faults.

[0027] For example, when the operating condition category label is "low-speed steady state," the time-domain feature set retrieved from the fault feature knowledge base includes two features: kurtosis index and peak factor. This is because vibration energy is low at low speeds, and kurtosis and peak factor are most sensitive to weak impacts. The wavelet packet energy feature set selects the normalized energy of the first, second, and third subspaces in the third-level decomposition, because fault feature frequencies are concentrated in the low-frequency band at low speeds. When the operating condition category label is "high-speed steady state," the time-domain feature set adds the impulse index and margin index, and the wavelet packet energy feature set selects the normalized energy of the fourth to eighth subspaces, because fault feature frequencies shift to the higher frequency band at high speeds. When the operating condition category label is "pressure jump transition," the time-domain feature set only retains the waveform index, and the wavelet packet energy feature set uses the normalized energy of all eight subspaces, because a wider bandwidth of energy distribution information is needed to distinguish normal pressure adjustment from real faults under transitional operating conditions. The dynamically output fault feature set will be passed to the weak signal enhancement module as a guide for subsequent enhancement extraction.

[0028] The operating condition classification model is trained using a self-organizing map neural network, which consists of an input layer and a competition layer. The input layer has six nodes, corresponding to the six dimensions of the operating condition parameter vector. The competition layer is a two-dimensional grid structure, initially set to 100 nodes (10 rows x 10 columns), with each node corresponding to a weight vector. This weight vector has the same dimensions as the input vector, i.e., six dimensions. The training process includes the following steps: Throughout the entire lifecycle of the reciprocating pump, including its healthy state and varying degrees of performance degradation, operating parameters are continuously collected using speed sensors, pressure sensors, flow sensors, and media characteristic sensors, forming a large number of operating parameter vector samples. Each sample is also manually labeled with its corresponding operating condition category, such as "low-speed steady state," "high-speed steady state," "pressure jump transition," "sudden increase in sand content transition," "viscosity increase transition," and "density decrease transition." A total of no fewer than ten thousand samples are collected, covering all typical operating conditions that a reciprocating pump may encounter.

[0029] A small random number is assigned to the six-dimensional weight vector of each node in the competition layer, with the random number ranging from zero to one. At the same time, the initial learning rate is set to 0.5, and the initial neighborhood radius is 5 (that is, covering all nodes in the competition layer within a range of 5 nodes from the current winning node).

[0030] For each input working condition parameter vector sample, calculate the Euclidean distance between this vector and the weight vector of each node in the competition layer, and find the node with the smallest distance as the winning node. Then, update the weight vectors of the winning node and all nodes in its neighborhood. The update formula is: the new weight vector equals the old weight vector plus the learning rate multiplied by the neighborhood function multiplied by the difference between the input vector and the old weight vector. The value of the neighborhood function decreases as the distance between the node and the winning node increases, and it adopts the form of a Gaussian function. The learning rate and neighborhood radius gradually decrease with the increase of the number of training iterations: after training 500 samples, the learning rate is multiplied by 0.95, and the neighborhood radius is halved, until the learning rate drops below 0.01 and the neighborhood radius drops to 1, at which point the process stops.

[0031] After iterative training on all training samples, adjacent nodes with similar weight vectors in the competition layer automatically aggregate to form different regions. The input samples corresponding to the nodes within each region are statistically analyzed, and the operating condition category with the largest number of samples in that region is taken as the category label for that region. This region is thus a operating condition category cluster. For example, the upper left region of the competition layer gathers all low-speed steady-state samples, so this region is labeled as the "low-speed steady-state cluster"; the upper right region gathers all high-speed steady-state samples, labeled as the "high-speed steady-state cluster"; and the middle region gathers samples during the pressure jump process, labeled as the "pressure jump transition cluster".

[0032] A working condition cluster refers to a set of multiple operating states with similar operating parameters. Each working condition cluster corresponds to a range of parameters for a set of stable or transitional states of a reciprocating pump under specific operating conditions. The working condition parameter vector consists of six parameters: rotational speed, output pressure, instantaneous flow rate, medium density, medium viscosity, and medium sand content, forming a six-dimensional vector. A self-organizing mapping neural network automatically organizes these six-dimensional vectors into different regions on a two-dimensional plane according to their similarity; each region is a working condition cluster. For example, the "low-speed steady-state cluster" contains multiple working condition parameter vectors with rotational speeds between 100 and 130 rpm and pressures between 10 and 15 MPa; the "sand content surge transition cluster" contains intermediate state vectors during the process of sand content rapidly increasing from 3% to 8%. Vectors within the same working condition cluster are close in Euclidean distance, representing similar operating conditions; vectors between different clusters are far apart, representing different operating conditions.

[0033] After training, for any working condition parameter vector collected in real time, it is input into the self-organizing map neural network. The Euclidean distance between the vector and the weight vector of each node in the competition layer is calculated. The nearest node is found, and the working condition category cluster to which the node belongs is the mapping result of the current working condition parameter. The corresponding working condition category label is output. This label is used to retrieve the corresponding fault feature set from the preset fault feature knowledge base.

[0034] The raw sensor signal refers to the vibration acceleration signal from multiple measuring points in each cylinder of the reciprocating pump, directly acquired by vibration sensors. Specifically, a piezoelectric accelerometer is installed at the valve box and cylinder head of each cylinder of the reciprocating pump. The sensor sensitivity is 100 millivolts per gram, the measurement range is -50 grams to +50 grams, and the frequency response range is 0.5 Hz to 5 kHz. Two measuring points are set for each cylinder, one at the suction valve cover and the other at the discharge valve cover. For a three-cylinder reciprocating pump, a total of six vibration sensors are installed. All sensors synchronously acquire vibration acceleration waveforms during the operation of each cylinder at a sampling frequency of 20,000 times per second. The raw sensor signal contains various components such as pump valve opening and closing impacts, piston reciprocating motion, fluid pulsation, and mechanical transmission. Among these, the signal component generated by weak faults has a very small amplitude and is usually submerged in background noise.

[0035] The health baseline refers to the vibration acceleration signal collected under standard operating conditions (rated speed, rated pressure, clean water medium) when the reciprocating pump is first put into operation after a new or major overhaul. The specific acquisition steps are as follows: After the reciprocating pump has been installed, commissioned, and run continuously and stably for two hours, confirm that all monitoring parameters (vibration amplitude, temperature, pressure pulsation) are within the factory specifications; at this point, the equipment is considered to be in a healthy state. Vibration acceleration signals are continuously collected for 60 seconds at a sampling rate of 20,000 times per second as the health baseline signal. This baseline signal is then subjected to three-level wavelet packet decomposition to obtain eight frequency band subspaces. The wavelet packet energy value of each subspace is calculated and stored as the health baseline energy vector. This baseline remains unchanged throughout the entire lifespan of the equipment unless a new baseline is established after a major overhaul.

[0036] Three-level wavelet packet decomposition is a process for multi-resolution analysis of the original sensor signal. The difference between wavelet packet decomposition and ordinary wavelet decomposition is that it decomposes not only the low-frequency component but also the high-frequency component simultaneously, thus obtaining a uniform frequency band division. The Daubechies fourth-order wavelet is selected as the mother wavelet. The first level of decomposition decomposes the signal into low-frequency approximation coefficients and high-frequency detail coefficients, corresponding to two frequency band subspaces: 0 to 5 kHz and 5 kHz to 10 kHz. The second level of decomposition further decomposes each of the two subspaces obtained in the first level into low-frequency and high-frequency components, resulting in four subspaces: 0 to 2,500 Hz, 2,500 to 5 kHz, 5,000 to 7,500 Hz, and 7,500 to 10 kHz. The third level of decomposition further decomposes each of the four subspaces in the second level, ultimately obtaining eight frequency band subspaces, each with a frequency width of 10 kHz divided by eight, which equals 1,250 Hz. The frequency ranges of the eight subspaces are: 0 to 1250 Hz, 1250 to 2500 Hz, 2500 to 3750 Hz, 3750 to 5000 Hz, 5000 to 6250 Hz, 6250 to 7500 Hz, 7500 to 8750 Hz, and 8750 to 10000 Hz.

[0037] For the currently acquired vibration acceleration signal segment, a three-layer wavelet packet decomposition, identical to that performed on the healthy baseline, is first performed to obtain eight frequency band subspaces. The wavelet packet energy of each subspace is then calculated. The energy of each subspace is defined as the sum of the squares of all wavelet packet coefficients within that subspace. Let the energy of the current i-th subspace be Ei, and the energy of the i-th subspace under the healthy baseline be E0i. The formula for calculating the energy change rate of the i-th subspace is: the absolute value of Ei minus E0i (within parentheses) divided by E0i, then multiplied by 100%. The energy change rate reflects the degree of increase or decrease in energy of the current signal relative to the healthy state in that frequency band. A larger change rate indicates a significant change in the vibration energy of that frequency band, potentially corresponding to the appearance of fault characteristic frequencies. For example, the energy change rate of a healthy reciprocating pump in the 0 to 1250 Hz frequency band is typically within 5%; however, when a slight leak occurs in the discharge valve of a cylinder, the energy in that frequency band may increase by more than 30%.

[0038] Kurtosis is a dimensionless time-domain statistical parameter that reflects the sharpness of a signal waveform and is extremely sensitive to impulsive components. The kurtosis is calculated as the fourth-order central moment of the signal divided by the fourth power of the signal's standard deviation. For pure Gaussian noise, the kurtosis value is three. When periodic impulsive components are present in the signal, the kurtosis value is greater than three; the stronger the impulsive component, the larger the kurtosis value. In reciprocating pump fault diagnosis, a kurtosis value greater than three is often used as a threshold value to determine the presence of abnormal impulsive components. This threshold of three is based on statistics: the fourth-order cumulant of a Gaussian distribution is zero, corresponding to a kurtosis value of three; any signal deviating from the Gaussian distribution will cause the kurtosis value to deviate from three. In engineering practice, a kurtosis value exceeding 3.5 is considered to indicate a significant impulsive component, and exceeding four is considered to indicate a serious fault. Here, a value greater than three is chosen as a preliminary condition for screening sensitive frequency bands to include those frequency bands that have deviated from the Gaussian noise background in the candidate range, avoiding the omission of early, weak faults. For example, the kurtosis value of the wavelet packet reconstructed signal in each frequency band subspace of a healthy reciprocating pump is usually between 2.8 and 3.2; when there is slight wear on the pump valve sealing surface, the kurtosis value of the corresponding frequency band may rise to 3.5 to 4.0.

[0039] The determination of sensitive frequency bands involves two sub-steps. The first sub-step calculates the rate of change of wavelet packet energy relative to a healthy baseline for each frequency band subspace, and simultaneously calculates the kurtosis index of the reconstructed signal for that subspace. The second sub-step combines the energy change rate and kurtosis index to select the three subspaces with the largest energy change rates, requiring that the kurtosis index of these three subspaces is greater than three. The frequency ranges corresponding to these subspaces are then determined as sensitive frequency bands. If any of the top three subspaces in terms of energy change rate has a kurtosis index less than three, the process continues to the next subspace with a relatively large energy change rate, until three subspaces satisfying the kurtosis index greater than three are selected. For example, in a certain monitoring, the energy change rates of the eight subspaces were ranked as follows: 3rd subspace, 45%; 1st subspace, 38%; 5th subspace, 32%; and 2nd subspace, 15%. The kurtosis index of the third subspace is 4.2, the kurtosis index of the first subspace is 3.8, and the kurtosis index of the fifth subspace is 2.9 (less than 3). Therefore, the fifth subspace is excluded, and the second subspace (kurtosis index 3.5) is selected. The final sensitive frequency bands are determined to be 2,500 to 3,750 Hz for the third subspace, 0 to 1,250 Hz for the first subspace, and 1,250 to 2,500 Hz for the second subspace. These three frequency ranges will serve as the target frequency bands for subsequent weak signal enhancement processing.

[0040] The sensitive timescale parameter is a time resolution parameter derived from the sensitive frequency band, used to guide subsequent signal segmentation. The extraction method is as follows: First, find the lowest frequency from the three sensitive frequency bands; for example, in the above example, the lowest frequency of the three sensitive frequency bands is 0 Hz (from the first subspace). Then calculate the ratio of this lowest frequency to the sampling frequency. When the sampling frequency is 20,000 Hz and the lowest frequency is 0 Hz, it cannot be directly calculated; in actual processing, the minimum center frequency of the sensitive frequency band is used instead. More specifically: take the minimum center frequency in the sensitive frequency band, i.e., the center frequency of the first subspace is 625 Hz, then the ratio is 625 divided by 20,000 equals 0.03125. Based on this ratio, determine the wavelet packet decomposition level, which is equal to the integer part of the logarithm (base 2, sampling frequency divided by the lowest frequency). Alternatively, a simplified engineering method can be used: directly use the period corresponding to the lowest frequency in the sensitive frequency band as the timescale reference. It is clarified here that the sensitive timescale parameter is equal to half the reciprocal of the lowest frequency in the sensitive frequency band, in seconds. For example, the lowest frequency is 625 Hz, and its period is 0.0016 seconds. Taking half of that, we get 0.0008 seconds, or 0.8 milliseconds. This time scale is used as the window length for subsequent segmentation, ensuring that each segment contains at least one complete lowest frequency vibration cycle.

[0041] The specific steps of the weak signal enhancement module in extracting nonlinear weak fault features from the original sensing signal and outputting the enhanced feature signal are as follows: Step 1: Based on the sensitive time scale parameter, the original sensing signal is segmented according to the time scale to obtain multiple signal segments; based on the sensitive time scale parameter, the original sensing signal is segmented according to the time scale to obtain multiple signal segments.

[0042] Step 2: Perform maximum and minimum value normalization on each signal segment to ensure all signal amplitudes fall between -1 and +1, resulting in a normalized signal segment. For each signal segment obtained in Step 1, perform the following operations: Find the maximum and minimum values ​​among all sampling points within the segment. Let the maximum value be Vmax and the minimum value be Vmin. For each sampling point X in the segment, calculate the normalized value using the formula: the normalized value equals (X minus Vmin) divided by (Vmax minus Vmin), multiplied by two, and then subtracted by one. This formula linearly maps the original amplitude range to the interval between -1 and +1. For example, in an original signal segment, the maximum value Vmax is 5 millivolts and the minimum value Vmin is -3 millivolts, so the amplitude range is 8 millivolts. A sampling point in the segment is 2 millivolts. First, subtract -3 to get 5, divide by 8 to get 0.625, multiply by 2 to get 1.25, and subtract 1 to get 0.25. After normalization, this sampling point becomes 0.25. After normalization of all segments, the waveform shape of the signal remains unchanged, but the amplitude is unified to a standard range. This is a typical requirement of bistable stochastic resonant systems for the amplitude range of the input signal, because the design of the potential well parameters and noise intensity parameters is usually based on a unit amplitude signal. The normalized signal segment is called the normalized signal segment.

[0043] Step 3: Input each normalized signal segment into a bistable stochastic resonance system. The system parameters of this bistable stochastic resonance system include potential well parameters and noise intensity parameters. The signal-to-noise ratio gain maximization criterion is used as the objective function for parameter optimization. The center frequency of the sensitive frequency band is taken as the characteristic frequency of the system's desired response. The search range for the potential well parameters is set to a minimum value of 0.1 to a maximum value of 10, and the search range for the noise intensity parameters is set to a minimum value of 0.01 to a maximum value of 100. The optimal potential well parameters and noise intensity parameters are determined through iterative search, bringing the system to a stochastic resonance state. This step is described in the following steps: Mathematical Model of a Bistable Stochastic Resonance System: A bistable stochastic resonance system is described by the Langevin equation, which consists of three parts: the potential well force, the input signal, and noise. The potential well force is determined by the potential well parameters, which control the potential barrier height and the potential well width. When the potential well parameters are small, the potential barrier is low, and particles can easily transition between the two potential wells; when the potential well parameters are large, the potential barrier is high, and particles are confined to a single potential well. The noise intensity parameter represents the intensity of the internal thermal noise or external noise of the system. The input of the system is a normalized signal segment, and the output is the signal enhanced by stochastic resonance. The stochastic resonance phenomenon occurs when the signal, noise, and nonlinear system achieve optimal matching. At this point, the weak periodic signal can be synergistically amplified by the energy of the noise, and the output signal-to-noise ratio reaches its maximum value.

[0044] The signal-to-noise ratio (SNR) gain maximization criterion is followed: SNR gain is defined as the ratio of the output SNR to the input SNR. The SNR is calculated by performing a Fast Fourier Transform (FFT) on the signal to obtain the power spectrum, extracting the signal power near the center frequency of the sensitive frequency band, and simultaneously extracting the average noise power in the surrounding neighborhood of that frequency; the ratio of these two is the SNR. The SNR maximization criterion refers to adjusting the potential well parameters and noise intensity parameters to maximize the increase in output SNR relative to input SNR. This criterion can automatically find the optimal combination of system parameters without prior knowledge of the specific noise intensity.

[0045] Obtaining the center frequency of the sensitive frequency band: Sensitive frequency bands typically have three frequency ranges, each with a center frequency. For example, when the sensitive frequency band is 0 to 1250 Hz, the center frequency is 625 Hz; when the sensitive frequency band is 1250 to 2500 Hz, the center frequency is 1875 Hz. The characteristic frequency of the system's desired response refers to the fault characteristic frequency that we want to amplify from the noise. This frequency is taken as the center frequency of the sensitive frequency band. Because the periodic impacts generated by the fault often fall within this frequency band, using the center frequency as the target allows the stochastic resonant system to selectively amplify that frequency component.

[0046] Search ranges for the potential well parameter and noise intensity parameter: The search range for the potential well parameter is set from a minimum value of 0.1 to a maximum value of 10. The lower limit of 0.1 ensures that the system has a significant bistable potential barrier (the barrier height is positive), avoiding degradation into a monostable system; the upper limit of 10 ensures that the barrier is not so high that the signal cannot induce a transition. The search range for the noise intensity parameter is set from a minimum value of 0.01 to a maximum value of 100. The lower limit of 0.01 corresponds to extremely low noise, at which point the system is close to a deterministic system, and random resonance phenomena are not obvious; the upper limit of 100 corresponds to high-intensity noise, at which point the signal is completely submerged in noise, and the output signal-to-noise ratio is extremely low. This range covers the noise levels that may be encountered in most industrial vibration signals.

[0047] Iterative search to determine optimal parameters: For each normalized signal segment, the following search process is executed: First, the potential well parameters and noise intensity parameters are sampled in a grid within their respective ranges, with step sizes set to 0.1 for the potential well parameters and 1 for the noise intensity parameters. For each parameter combination, the normalized signal segment is input into a bistable stochastic resonance system, and the Langevin equation is numerically solved (using the fourth-order Runge-Kutta method, with the time step taken as one-tenth of the sampling interval) to obtain the output signal. The signal-to-noise ratio (SNR) gain of the output signal is calculated. After traversing all parameter combinations, the set of parameters with the largest SNR gain is selected as the optimal potential well parameters and optimal noise intensity parameters for the current signal segment. Then, using this set of parameters as the center, the search range is reduced to 20% of the original range, and a fine grid search is performed again (step size halved) to further improve the accuracy of the optimal parameters. The iteration terminates when the rate of change of the SNR gain corresponding to the optimal parameters obtained in two consecutive iterations is less than 1%. The final potential well parameters and noise intensity parameters bring the system into a stochastic resonance state, i.e., the weak periodic signal and noise work together to maximize the output of signal energy.

[0048] For example, suppose the current normalized signal segment originates from the vibration signal of a cylinder block, with a sensitive frequency band center frequency of 625 Hz. After the first coarse search, it was found that the signal-to-noise ratio (SNR) gain reaches its maximum of 12 dB when the potential well parameter is 1.5 and the noise intensity parameter is 8. The second fine search is performed within the range of potential well parameters from 1.2 to 1.8 and noise intensity parameters from 6.4 to 9.6, with the step size halved. Finally, the optimal potential well parameter is determined to be 1.35, and the optimal noise intensity parameter is 7.8, achieving an SNR gain of 13.2 dB. At this point, the amplitude of the impact component near 625 Hz in the system output signal is approximately four times greater than the amplitude of the same component in the input signal.

[0049] Step 4: The output signals of each normalized signal segment after random resonance processing are spliced ​​together in their original time order to obtain the random resonance enhanced signal. After the processing in Step 3, each normalized signal segment corresponds to an output signal segment enhanced by random resonance. The length of this output signal segment is the same as that of the input segment (sixteen sampling points). These output signal segments are spliced ​​together sequentially according to their time order in the original signal to form a complete output signal sequence. During splicing, it is ensured that adjacent segments are continuous and seamless, without overlap or gaps. The spliced ​​signal is called the random resonance enhanced signal. The amplitude of the periodic weak impact components (such as weak impacts caused by micro-leakage in pumps and valves, and periodic pulses caused by early wear of bearings) that were originally masked by noise in this signal is significantly improved, increasing the signal-to-noise ratio by more than ten decibels.

[0050] A three-layer wavelet soft thresholding denoising process is applied to the stochastic resonance enhanced signal to eliminate residual noise and obtain the final enhanced feature signal. Although the stochastic resonance enhanced signal has already enhanced the target signal components, it may still contain some broadband noise, especially random noise in the high-frequency band. To further improve signal quality, a three-layer wavelet soft thresholding denoising process is adopted. The specific process is as follows: Daubechies fourth-order wavelet is selected as the mother wavelet, and the stochastic resonance enhanced signal is decomposed into three layers of detail coefficients and one layer of approximation coefficients. An unbiased risk estimation thresholding method is used to calculate an independent threshold for each layer of detail coefficients. For each layer of detail coefficients, the threshold amount is reduced to zero for coefficients whose absolute value is greater than the threshold, and the threshold amount is set to zero for coefficients whose absolute value is less than or equal to the threshold. This is the soft thresholding operation. The processed detail coefficients and the unprocessed approximation coefficients are then reconstructed using wavelet denoising to obtain the denoised signal. Soft thresholding can effectively suppress Gaussian white noise while preserving the signal's peak characteristics (such as impulse pulses). The denoised signal is the final enhanced feature signal, which is output to the coupling-decoupling analysis module. The signal-to-noise ratio of the weak fault features in this signal is increased by more than 15 dB compared to the original sensing signal, and the nonlinear weak fault components are significantly enhanced.

[0051] The enhanced feature signals received by the coupling and decoupling analysis module are the enhanced feature signals corresponding to multiple measuring points in each cylinder of the reciprocating pump; Multi-source information coupling and decoupling analysis refers to separating the coupling components in the enhanced characteristic signals of each cylinder block measuring point, eliminating the mutual interference between vibration signals of adjacent cylinder blocks, and outputting a decoupled independent fault characteristic value sequence corresponding to each cylinder block.

[0052] The specific steps for the coupling / decoupling analysis module to output the decoupled independent fault characteristic value sequence and fault propagation path are as follows: The enhanced characteristic signals of each cylinder block measuring point are constructed into an observation signal matrix. Using the independent component analysis algorithm with kurtosis as the comparison function, the unmixing matrix is ​​calculated iteratively to decompose the observation signal matrix into mutually statistically independent source signal components. Each source signal component corresponds to the vibration contribution of a cylinder block, thus obtaining the independent fault characteristic value sequence of each cylinder block. The observation signal matrix is ​​constructed as follows: Vibration sensors are arranged in the valve box and cylinder head of each cylinder of the reciprocating pump, and each measuring point corresponds to one enhanced characteristic signal. Assuming that the reciprocating pump is a three-cylinder structure, two measuring points are set in each cylinder (located in the suction valve cover and the discharge valve cover respectively), then there are a total of six enhanced characteristic signals. Each enhanced characteristic signal is the final enhanced characteristic signal output after being processed by the weak signal enhancement module, and the signal length is L sampling points (for example, L equals 32,000 points, corresponding to 1.6 seconds). These six signals are denoted as X1(t), X2(t), X3(t), X4(t), X5(t), and X6(t), respectively, where t=1,2,...,L are the sampling times. These six signals are arranged in rows to form an observation signal matrix X with six rows and L columns, where the first row is X1(1) to X1(L), the second row is X2(1) to X2(L), and so on. Each row in the observation signal matrix X represents the vibration signal acquired by one sensor channel, and each column represents the sampled values ​​of all sensor channels at the same time. This matrix contains mixed information between the vibration signals of each cylinder block, because the signal acquired at each measuring point is actually a linear superposition of the vibration contributions of multiple cylinder blocks.

[0053] The basic principle of Independent Component Analysis (ICA) is that the observed signal is formed by mixing several unknown source signals through a linear mixing matrix, and the source signals are statistically independent of each other. In the vibration signal of a reciprocating pump with multiple cylinders, the vibration sources of each cylinder (such as pump valve impact and piston motion) are driven by different crankshaft phases and are statistically independent of each other. The goal of ICA is to find a demixing matrix W such that the output signal Y obtained after linearly transforming the observed signal matrix X is equal to W multiplied by X, and the components of each row of Y are statistically independent. The output signal Y is the estimate of the source signal components, and each source signal component corresponds to the vibration contribution of one cylinder.

[0054] This algorithm uses the absolute value of kurtosis as the comparison function, aiming to maximize the absolute kurtosis value of each component of the output signal. Kurtosis is a statistical measure of the non-Gaussianity of a signal. For a discrete signal sequence y, the kurtosis K is calculated as follows: K equals the average of the fourth power of all sample points divided by the square of the average of the squares of all sample points, minus three. This definition makes the kurtosis value of Gaussian noise zero. The larger the absolute value of kurtosis, the stronger the non-Gaussianity of the signal, and thus the stronger the statistical independence between signals.

[0055] The iterative calculation of the unmixing matrix W employs a fixed-point iterative algorithm. First, each row of the observed signal matrix X is centered by subtracting its mean, making the mean of each row zero. Then, whitening is performed by linearly transforming the centered observed signal matrix into an identity matrix, resulting in a whitened signal matrix Z. The unmixing matrix W is initialized as a six-row, six-column identity matrix. For each row vector w of the unmixing matrix W (each w corresponds to a source signal component to be extracted), iterative updates are performed according to the following formula: the new value of w equals the average of the cube of the whitened signal matrix Z multiplied by the cube of (w's transpose multiplied by Z), minus three times w. After each update, w is orthogonalized: w equals w minus the transpose of the updated rows in W multiplied by the current w, then multiplied by the cumulative sum of the updated rows. Orthogonalization ensures that the extracted source signal components are uncorrelated. After orthogonalization, w is normalized by dividing it by its Euclidean norm, so that the unit length of w is 1. Following the above update steps, the unmixed vector w1 corresponding to the first source signal component, the unmixed vector w2 corresponding to the second source signal component, and so on, up to the unmixed vector w6 corresponding to the sixth source signal component. After each w is extracted, it is used as an updated row in the orthogonalization process of subsequent vectors.

[0056] The iteration termination condition adopts the following two criteria; the iteration terminates if either one is satisfied. Criterion 1: The rate of change threshold of the unmixed matrix. The rate of change between the unmixed matrices W_old and W_new obtained from two consecutive iterations is defined as follows: calculate the Frobenius norm of the matrix W_new minus W_old, divide by the Frobenius norm of W_old, and then multiply by 100%. When this rate of change is less than 0.1%, the unmixed matrix is ​​considered to have converged, and the iteration terminates. The Frobenius norm is the square root of the sum of the squares of all elements of the matrix, and can comprehensively reflect the overall change of the matrix.

[0057] Rule 2: Maximum number of iterations, set to 1000. When the number of iterations reaches 1000, the iteration is forcibly terminated regardless of whether the rate of change of the unmixed matrix meets the threshold, and the current unmixed matrix is ​​output as the final result. In typical three-cylinder reciprocating pump vibration signal processing, the rate of change threshold condition is usually met between 150 and 300 iterations.

[0058] After independent component analysis, the unmixed matrix W is obtained. W is multiplied by the observed signal matrix X (in actual calculations, it is multiplied by the whitened signal matrix Z) to obtain the output signal matrix Y, which has the same dimensions as X, six rows by L columns. Each row of the output signal matrix Y represents an estimate of a source signal component. The correspondence between the source signal components and the cylinder blocks is established based on the synchronization relationship between the impact phase of each source signal component and the crankshaft angle: cross-correlation analysis is performed on each source signal component and the crankshaft speed signal to calculate the time interval between the impact pulses in the source signal components. For each revolution of the crankshaft of the three-cylinder reciprocating pump (360 degrees), the three cylinder blocks sequentially complete one intake and exhaust stroke, with the crankshaft phase difference between the pump valve impacts of each cylinder block being 120 degrees. By extracting the phase information of the impact pulses in the source signal components, they are matched to the corresponding cylinder block numbers. After matching, each cylinder block corresponds to one or more source signal components. The source signal components belonging to the same cylinder block are superimposed to obtain the independent fault characteristic value sequence of that cylinder block. This sequence is a time series of length L, reflecting the vibration intensity change of that cylinder block on the time axis.

[0059] Using the obtained independent fault characteristic value sequences of each cylinder block as input, the cross-correlation function of the independent fault characteristic value sequences between adjacent cylinder blocks is calculated in time order. For two independent fault characteristic value sequences of adjacent cylinder blocks, the cylinder blocks are numbered as cylinder block 1, cylinder block 2, and cylinder block 3 according to the crankshaft rotation order. Adjacent pairs include (cylinder block 1, cylinder block 2) and (cylinder block 2, cylinder block 3). For each pair of adjacent cylinder blocks, two sequences are extracted, denoted as sequence A (corresponding to the previous numbered cylinder block) and sequence B (corresponding to the next numbered cylinder block), both sequences having a length of L sampling points. The formula for calculating the cross-correlation function is as follows: for each time delay τ (the value of τ ranges from negative L plus one to L minus one), the cross-correlation function value R(τ) is equal to the sum of the value of sequence A at time t multiplied by the value of sequence B at time t plus τ. When τ is positive, it indicates that sequence A precedes sequence B (i.e., the change of A occurs before B); when τ is negative, it indicates that sequence B precedes sequence A (i.e., the change of B occurs before A). To eliminate the influence of the magnitude of the sequence's own amplitude on the correlation value, a normalized cross-correlation function is used, which is R(τ) divided by the product of the standard deviation of sequence A and sequence B and the sequence length. The normalized cross-correlation function value ranges from negative one to positive one.

[0060] The preset cross-correlation threshold is based on both statistical hypothesis testing and engineering experience. Regarding statistical significance testing, under the null hypothesis (no correlation between the two sequences), for independent and identically distributed Gaussian sequences of length L, the distribution of their cross-correlation function values ​​approximates a normal distribution with a mean of zero and a standard deviation of 1 / 1 square root of L. Taking a significance level of 5%, the corresponding critical value is approximately 1.96 divided by the square root of L. For example, when L equals 32,000, the square root of L is approximately 179, and 1.96 divided by 179 is approximately 0.011. That is, when the peak value of the normalized cross-correlation function exceeds 0.011, there is a 95% confidence that a significant correlation exists between the two sequences. Considering the robustness of engineering applications, the actual threshold is taken as five to ten times this theoretical value.

[0061] Regarding engineering experience thresholds, in practical applications of reciprocating pump fault diagnosis, the peak value of the normalized cross-correlation function between independent fault characteristic value sequences of adjacent cylinders under healthy conditions is typically between 0.02 and 0.05. When a fault occurs in a cylinder, the peak value of the cross-correlation function rises above 0.1. Based on the above, the preset cross-correlation threshold is set to 0.08. When the peak value of the normalized cross-correlation function exceeds 0.08, a significant fault propagation relationship is considered to exist between the two cylinders; when the peak value is below or equal to 0.08, a significant fault propagation relationship is considered not to exist.

[0062] For each pair of adjacent cylinders, calculate the normalized cross-correlation function R(τ), and find the maximum value of R(τ) and its corresponding delay time τ_max. If R(τ_max) is greater than the preset cross-correlation threshold of 0.08, a fault propagation relationship is determined. Then, the propagation direction is determined based on the sign of τ_max: if τ_max is positive, it means that the independent fault feature value sequence of the preceding cylinder number precedes the sequence of the following cylinder number, that is, the fault feature change of the preceding cylinder appears first, and the fault feature change of the following cylinder appears later. Therefore, it is determined that the fault propagates from the preceding cylinder to the following cylinder. If τ_max is negative, it means that the sequence of the following cylinder number precedes the sequence of the preceding cylinder number, that is, the fault feature change of the following cylinder appears first, and the fault feature change of the preceding cylinder appears later. Therefore, it is determined that the fault propagates from the following cylinder to the preceding cylinder. Connect the propagation directions of each adjacent pair in sequence to form a complete fault propagation path. For example: For adjacent cylinder block 1 and cylinder block 2, the calculated peak value of the normalized cross-correlation function R_max is 0.35, which exceeds the threshold of 0.08. The corresponding time delay τ_max is positive 0.002 seconds, indicating that the fault propagates from cylinder block 1 to cylinder block 2. For adjacent cylinder block 2 and cylinder block 3, the calculated R_max is 0.28, and τ_max is negative 0.0015 seconds, indicating that the fault propagates from cylinder block 3 to cylinder block 2. Therefore, the fault propagation path is cylinder block 1 to cylinder block 2 and cylinder block 3 to cylinder block 2.

[0063] The fault propagation path is output as a sequence of cylinder numbers, indicating the initial fault source cylinder and the subsequent cylinders affected by the propagation. Based on the propagation directions of each adjacent pair obtained above, a directed graph is constructed with cylinders as nodes and propagation directions as edges. Nodes without incoming edges (i.e., no other cylinders from which the fault propagation reaches) are identified in this directed graph; these nodes are the initial fault source cylinders. Starting from each initial fault source cylinder, all reachable nodes are traversed along the propagation direction, resulting in a sequence of cylinder numbers originating from the initial fault source cylinder. Subsequent nodes in this sequence are the cylinders affected by the propagation. The fault propagation path is output as a sequence of cylinder numbers, starting from the initial fault source cylinder and listing the subsequent cylinders affected by the propagation in the order of the fault propagation, until no further propagation occurs.

[0064] It should be noted that when there are multiple propagation directions, such as simultaneously detecting propagation from cylinder 1 to cylinder 2 and from cylinder 3 to cylinder 2, the peak values ​​of the cross-correlation functions of each pair of adjacent cylinders are compared, and the propagation direction corresponding to the pair with the largest peak value is selected as the main propagation path. Only the cylinder number sequence corresponding to the main propagation path is output.

[0065] The specific steps of the full lifecycle mapping module in generating quantitative indicators of the current health status of the reciprocating pump based on the independent fault feature value sequence and fault propagation path are as follows: For each cylinder block, its current independent fault feature value sequence is compared with the health baseline feature value sequence to calculate the health score of each cylinder block. The health baseline feature value sequence is the independent fault feature value sequence of each cylinder block in a healthy state for the reciprocating pump. The health score is equal to the ratio of the Euclidean distance between the current independent fault feature value sequence and the health baseline feature value sequence to the sum of the magnitudes of the health baseline feature value sequence and the current independent fault feature value sequence, and is zero when the result is less than zero. The independent fault feature value sequence of each cylinder block is obtained from the coupling-decoupling analysis module and is a time series of fixed length, for example, 1,000 points. The health baseline feature value sequence is pre-acquired and stored through the same coupling-decoupling analysis process when the reciprocating pump is first put into normal operation and confirmed to be without any faults; it is also a sequence of 1,000 points. The Euclidean distance between the two sequences is defined as: squaring the numerical difference at each corresponding position, adding all the squared values, and then taking the square root. The smaller the Euclidean distance, the closer the vibration characteristics of the current cylinder block are to a healthy state; the larger the Euclidean distance, the more obvious the degradation. The denominator is the sum of the modulus of the health baseline characteristic value sequence (i.e., the square root of the sum of the squares of all points in that sequence) and the modulus of the current independent fault characteristic value sequence. This design ensures that when the current sequence is exactly the same as the baseline sequence, the Euclidean distance is zero, and the health score is equal to one; when the current sequence differs greatly from the baseline sequence, the Euclidean distance approaches the larger of twice the modulus, and the health score approaches zero. If the calculation result is negative, it is set to zero, ensuring that the health score always remains within the range of zero to one.

[0066] For example, consider a three-cylinder reciprocating pump. The health baseline characteristic value sequence of cylinder one has a magnitude of 10.0, and the current independent fault characteristic value sequence has a magnitude of 10.5. The Euclidean distance between the two sequences is 2.0. The denominator is 10.0 plus 10.5, which equals 20.5. The ratio 2.0 divided by 20.5 is approximately 0.098. The health score is 1 minus 0.098, which equals 0.902, indicating that the cylinder block is in good health. If the Euclidean distance of the other cylinder is 15.0, and the magnitudes are 10.0 and 18.0 respectively, the sum of the denominators is 28.0, the ratio is 0.536, and the health score is 0.464, indicating that the performance of this cylinder block has significantly deteriorated.

[0067] Based on the fault propagation path, the sequential position of each cylinder block within the propagation path is determined. The health score of cylinders later in the propagation sequence is corrected by multiplying it by a propagation attenuation coefficient. The propagation attenuation coefficient ranges from 0.8 to 1, decreasing with increasing propagation distance. The fault propagation path is output by the coupling / decoupling analysis module and is presented as a sequence of cylinder block numbers, such as "Cylinder Block 1, Cylinder Block 2, Cylinder Block 3," indicating that the fault propagates from Cylinder Block 1 to Cylinder Block 2, and then to Cylinder Block 3. Based on this sequence, the propagation sequence position of each cylinder block is determined: Cylinder Block 1 is position 1, Cylinder Block 2 is position 2, and Cylinder Block 3 is position 3.

[0068] The propagation attenuation coefficient is used to correct the health score of cylinders affected by propagation. Its physical meaning is that the additional damage to subsequent cylinders caused by a fault gradually weakens during propagation. Therefore, the health score of subsequent cylinders should not be directly included in the overall pump index based on the original calculated value, but should be multiplied by a coefficient less than one to reflect that the health of that cylinder is partially affected by the upstream fault. The specific setting method for the propagation attenuation coefficient is as follows: the propagation distance is defined as the current cylinder's sequential position minus the sequential position of the initial fault source cylinder. For example, if the initial fault source cylinder is cylinder one (sequential position one), then the propagation distance of cylinder one is zero, and it is not multiplied by the attenuation coefficient (or multiplied by one); the propagation distance of cylinder two is one, and the attenuation coefficient is 0.9; the propagation distance of cylinder three is two, and the attenuation coefficient is 0.8. The attenuation coefficient ranges from 0.8 to 1. For every increase of one in the propagation distance, the attenuation coefficient decreases by 0.1, but it is never lower than 0.8. If there are only two cylinders in the propagation path, only cylinder two is multiplied by 0.9. If cylinder 1 propagates to cylinder 2 in the propagation path, and cylinder 3 also propagates to cylinder 2 at the same time (i.e., there are multiple propagation directions), then the sequential position is determined according to the main propagation path (the one with the largest peak value of the cross-correlation function), and the cylinders on the non-main path have their attenuation coefficients calculated separately according to their respective distances relative to the initial source.

[0069] For example, the fault propagation path is "Cylinder Block 2, Cylinder Block 3, Cylinder Block 1", meaning the initial fault source is Cylinder Block 2, propagates to Cylinder Block 3, and then to Cylinder Block 1. The propagation distance of Cylinder Block 2 is zero, and the health score remains unchanged at the original calculated value of 0.7. The propagation distance of Cylinder Block 3 is one, multiplied by an attenuation coefficient of 0.9, resulting in a corrected health score of 0.63. The propagation distance of Cylinder Block 1 is two, multiplied by an attenuation coefficient of 0.8, resulting in a corrected health score of 0.56 (assuming the original score is 0.7).

[0070] The corrected health scores of all cylinders are weighted and averaged. The weight of each cylinder is pre-assigned according to its propagation order in the fault propagation path, with cylinders appearing earlier in the propagation order having lower weights. This yields a quantitative health status index between zero and one, which represents the current degree of degradation throughout the reciprocating pump's entire lifespan. The weighting principle for the weighted average is: cylinders appearing earlier in the propagation order (i.e., closer to the initial fault source) contribute less to the overall pump health index; cylinders appearing later in the propagation order contribute more.

[0071] The rationale behind this design is that the initial health score of the fault-causing cylinder accurately reflects the root cause of the fault. However, the health scores of subsequent cylinders may contain inaccuracies due to propagation effects. Therefore, higher weights are assigned to subsequent cylinders to highlight the cumulative effect of overall pump performance degradation. The weighting method is as follows: Assuming there are N cylinders in the fault propagation path, number them from one to N according to the propagation order. Define the base weight as the square of the sequential position number. For example, the cylinder at position one has a weight of one squared, equal to one; the cylinder at position two has a weight of two squared, equal to four; and the cylinder at position three has a weight of three squared, equal to nine. Then, the weights of each cylinder are normalized, meaning the final weight of each cylinder is equal to its base weight divided by the sum of the base weights of all cylinders. After normalization, the sum of the weights of all cylinders is one.

[0072] For example, the fault propagation path is "Cylinder Block 1, Cylinder Block 2, Cylinder Block 3", where N equals 3. The base weights are 1, 4, and 9, totaling 14. After normalization, the weights are: Cylinder Block 1 = 1 divided by 14, approximately 0.0714; Cylinder Block 2 = 4 divided by 14, approximately 0.2857; Cylinder Block 3 = 9 divided by 14, approximately 0.6429. Assume the corrected health scores are: Cylinder Block 1 0.90, Cylinder Block 2 0.80, Cylinder Block 3 0.70. The weighted average is calculated as: 0.90 multiplied by 0.0714 plus 0.80 multiplied by 0.2857 plus 0.70 multiplied by 0.6429, equaling 0.0643 plus 0.2286 plus 0.4500, which equals 0.7429. The health status quantification index is 0.74, indicating that the overall health of the reciprocating pump is currently approximately 74%, in a mild to moderate decline stage. This indicator gradually decreases as the equipment operates for longer periods. When the indicator falls below a preset alarm threshold (e.g., 0.4), the system can issue a maintenance warning.

[0073] The above-mentioned models or function formulas are all dimensionless and numerical calculations. The models or function formulas are obtained by software simulation based on a large amount of collected data to obtain the most recent real situation. The preset parameters in the models or function formulas are set by those skilled in the art according to the actual situation.

[0074] It should be understood that in the various embodiments of this application, the order of the above-mentioned processes does not imply the order of execution. The execution order of each process should be determined by its function and internal logic, and should not constitute any limitation on the implementation process of the embodiments of this application.

[0075] Those skilled in the art will recognize that the units and algorithm steps of the various examples described in conjunction with the embodiments disclosed herein can be implemented in electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are implemented in hardware or software depends on the specific application and design constraints of the technical solution. Those skilled in the art can use different methods to implement the described functions for each specific application, but such implementation should not be considered beyond the scope of this application.

[0076] Those skilled in the art will understand that, for the sake of convenience and brevity, the specific working processes of the systems, devices, and units described above can be referred to the corresponding processes in the foregoing method embodiments, and will not be repeated here.

[0077] The above are merely specific embodiments of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

Claims

1. An online testing system for the performance degradation of a reciprocating pump throughout its entire life cycle, characterized in that, include: The dynamic feature matching module is used to identify changes in the operating conditions of reciprocating pumps in real time and dynamically output the corresponding fault feature set. The weak signal enhancement module is used to receive the fault feature set, enhance and extract nonlinear weak fault features from the original sensing signal, and output the enhanced feature signal. The coupling and decoupling analysis module is used to perform multi-source information coupling and decoupling analysis on the enhanced feature signal, and outputs the decoupled independent fault feature value sequence and fault propagation path; The full lifecycle mapping module is used to generate quantitative indicators of the current health status of reciprocating pumps based on the independent fault characteristic value sequence and fault propagation path mapping.

2. The online testing system for the performance degradation of a reciprocating pump throughout its entire life cycle, as described in claim 1, is characterized in that... The operating conditions of the reciprocating pump obtained by the dynamic feature matching module include at least the speed, output pressure, instantaneous flow rate, and characteristics of the conveyed medium; the characteristics of the conveyed medium include the medium density, viscosity, and sand content.

3. The online testing system for the performance degradation of a reciprocating pump throughout its entire life cycle, as described in claim 2, is characterized in that... The specific steps of the dynamic feature matching module to identify changes in the operating conditions of the reciprocating pump in real time and dynamically output the corresponding fault feature set are as follows: The real-time speed, output pressure, instantaneous flow rate, and medium density, viscosity, and sand content of the reciprocating pump are continuously collected by speed sensors, pressure sensors, flow sensors, and medium characteristic sensors to form a sequence of operating parameters. A sliding time window is used to determine the steady state of the operating condition parameter sequence. The coefficient of variation of each parameter within the time window is calculated. When the coefficients of variation of all parameters are lower than the corresponding preset variation threshold, it is determined to be a steady state operating condition; otherwise, it is determined to be a dynamic operating condition. Based on the steady-state or dynamic discrimination results, the pre-trained working condition classification model is invoked to output the working condition category label; Based on the operating condition category label, the corresponding fault feature set is retrieved from the preset fault feature knowledge base and dynamically output; the fault feature set includes the time domain feature set and the wavelet packet energy feature set.

4. The online testing system for the performance degradation of a reciprocating pump throughout its entire life cycle, as described in claim 3, is characterized in that... The operating condition classification model uses a self-organizing map neural network to map the current operating condition parameter vector to the operating condition category cluster.

5. The online testing system for the performance degradation of a reciprocating pump throughout its entire life cycle, as described in claim 1, is characterized in that... The specific method by which the weak signal enhancement module receives the fault feature set is as follows: The fault feature set is analyzed to extract the sensitive frequency band and sensitive time scale parameters related to weak faults; among them, the original sensing signal is the vibration acceleration signal of multiple measuring points of each cylinder of the reciprocating pump collected by the vibration sensor. The sensitive frequency band is extracted as follows: taking the vibration acceleration signal of the reciprocating pump when it is running in a healthy state as a reference, the original sensing signal is decomposed into three layers of wavelet packets to obtain eight frequency band subspaces. The rate of change of wavelet packet energy of each subspace relative to the healthy reference is calculated, and the frequency ranges corresponding to the first three subspaces with kurtosis index greater than three are determined as the sensitive frequency bands. The sensitive time scale parameter is extracted as follows: the wavelet packet decomposition level is determined based on the ratio of the lowest frequency in the sensitive frequency band to the sampling frequency, and the time resolution corresponding to the level is used as the sensitive time scale parameter.

6. The online testing system for the performance degradation of a reciprocating pump throughout its entire life cycle, as described in claim 5, is characterized in that... The specific steps of the weak signal enhancement module in extracting nonlinear weak fault features from the original sensing signal and outputting the enhanced feature signal are as follows: Based on the sensitive time scale parameter, the original sensing signal is segmented according to the time scale to obtain multiple signal segments; Each signal segment is normalized by its maximum and minimum values ​​so that all signal amplitudes are between negative one and positive one, thus obtaining a normalized signal segment. Each normalized signal segment is input into a bistable stochastic resonance system. The system parameters of the bistable stochastic resonance system include potential well parameters and noise intensity parameters. The signal-to-noise ratio gain maximization criterion is used as the objective function for parameter optimization. The center frequency of the sensitive frequency band is used as the characteristic frequency of the system's expected response. The search range of the potential well parameters is set from a minimum value of 0.1 to a maximum value of 10, and the search range of the noise intensity parameters is set from a minimum value of 0.01 to a maximum value of 100. The optimal potential well parameters and noise intensity parameters are determined through iterative search, so that the system is in a stochastic resonance state. The output signals of each normalized signal segment after random resonance processing are spliced ​​together in the original time sequence to obtain the random resonance enhanced signal. The random resonance enhancement signal is subjected to three-layer wavelet soft thresholding denoising to obtain the final enhanced feature signal.

7. The online testing system for the performance degradation of a reciprocating pump throughout its entire life cycle, as described in claim 6, is characterized in that... The enhanced feature signals received by the coupling and decoupling analysis module are the enhanced feature signals corresponding to multiple measuring points in each cylinder of the reciprocating pump; Multi-source information coupling and decoupling analysis refers to separating the coupling components in the enhanced characteristic signals of each cylinder block measuring point, eliminating the mutual interference between vibration signals of adjacent cylinder blocks, and outputting a decoupled independent fault characteristic value sequence corresponding to each cylinder block.

8. The online testing system for the performance degradation of a reciprocating pump throughout its entire life cycle, as described in claim 7, is characterized in that... The specific steps for the coupling / decoupling analysis module to output the decoupled independent fault characteristic value sequence and fault propagation path are as follows: The enhanced characteristic signals of each cylinder block measuring point are constructed into an observation signal matrix. Using the independent component analysis algorithm with kurtosis as the comparison function, the unmixing matrix is ​​calculated iteratively to decompose the observation signal matrix into mutually statistically independent source signal components. Each source signal component corresponds to the vibration contribution of a cylinder block, thus obtaining the independent fault characteristic value sequence of each cylinder block. Using the obtained independent fault feature value sequence of each cylinder as input, the cross-correlation function of the independent fault feature value sequence between adjacent cylinders is calculated in time order. When the peak value of the cross-correlation function exceeds the preset cross-correlation threshold and there is a time delay, the propagation direction is determined according to the sign of the time delay corresponding to the peak value of the cross-correlation function: if the delay is positive, it is determined that the fault propagates from the previous cylinder to the next cylinder; if the delay is negative, it is determined that the fault propagates from the next cylinder to the previous cylinder. The propagation directions are connected in sequence to form the fault propagation path. The fault propagation path is output in the form of cylinder block number sequence, indicating the initial fault source cylinder block and the subsequent cylinder block sequence affected by the propagation.

9. The online testing system for the performance degradation of a reciprocating pump throughout its entire life cycle, as described in claim 8, is characterized in that... The specific steps of the full lifecycle mapping module in generating quantitative indicators of the current health status of the reciprocating pump based on the independent fault feature value sequence and fault propagation path are as follows: For each cylinder, its current independent fault feature value sequence is compared with the health baseline feature value sequence, and the health score of each cylinder is calculated. The health baseline feature value sequence is the independent fault feature value sequence of each cylinder of the reciprocating pump in a healthy state. The health score is equal to one minus the ratio of the Euclidean distance between the current independent fault feature value sequence and the health baseline feature value sequence to the sum of the magnitude of the health baseline feature value sequence and the magnitude of the current independent fault feature value sequence, and the value is zero when the calculation result is less than zero; Based on the fault propagation path, determine the sequential position of each cylinder in the propagation path, and correct the health score of the cylinder that is later in the propagation sequence by multiplying it by the propagation attenuation coefficient. The propagation attenuation coefficient ranges from 0.8 to 1. The health scores of all cylinders after correction are weighted and averaged. The weight of each cylinder is pre-assigned according to its propagation order in the fault propagation path to obtain a quantitative health status index between zero and one.

Images