A method and device for predicting oil equipment failure based on big data

Abstract

The application discloses a kind of based on big data's petroleum equipment fault prediction method and device, and the present application relates to equipment fault prediction technical field.The steps of the method include: collecting the original space-time data of each petroleum equipment and carrying out space-time alignment, obtain pre-processing petroleum equipment space-time data;Using variable time lag transfer entropy algorithm is analyzed to pre-processing petroleum equipment space-time data, obtain causal relationship, and generate space-time behavior association chain;Using sliding window algorithm is analyzed to space-time behavior association chain locally, obtain fault evolution trend value, and generate fault evolution sequence;Using bayesian variable point detection algorithm is identified to fault evolution sequence mutation point, obtain fault key turning point;Construct fault prediction model, and input fault key turning point and fault evolution sequence, output obtains cross-petroleum equipment fault prediction result.The present application generates cross-petroleum equipment fault prediction result by big data, realizes the prediction of petroleum equipment fault.

Description

Technical Field

[0001] This invention relates to the field of equipment failure prediction technology, specifically to a method and apparatus for predicting failures in petroleum equipment based on big data. Background Technology

[0002] In oil extraction, transportation, and processing, oil equipment operates under complex conditions of high temperature, high pressure, and high load for extended periods. The stability of its operation directly impacts production efficiency and operational safety. As the oil industry develops towards large-scale and intelligent operations, the degree of equipment clustering is continuously increasing. Equipment forms close coupling relationships through media flow and energy transfer. If a potential failure in a single piece of equipment is not promptly detected, it can rapidly spread to related upstream and downstream equipment, triggering a chain reaction of failures and causing severe economic losses and safety risks. Therefore, accurately identifying key turning points in the evolution of oil equipment failures and predicting failure trends and their impact range in advance has become a core requirement for ensuring the stable operation of oil production systems.

[0003] Existing methods for predicting faults in oilfield equipment typically begin by collecting physical parameters during equipment operation via sensor networks. After simple preprocessing such as noise reduction and missing value imputation, the core approach uses a fixed threshold method to identify critical inflection points in the fault. This method pre-sets fixed safety thresholds for various monitoring parameters based on historical operating experience. During real-time monitoring, the preprocessed parameter data is continuously compared with the preset thresholds. Once the data exceeds the threshold range, that moment is identified as a critical inflection point in the fault evolution, triggering the corresponding fault prediction process.

[0004] However, during the evolution of oil equipment failures, critical turning points often manifest as signal features that are difficult to capture through conventional numerical comparisons, such as sudden changes in behavior patterns and trend jumps. Judgment methods based on fixed thresholds can only mechanically compare whether parameter values ​​exceed preset ranges and cannot identify these structural change features that reflect the essence of the failure. As a result, it is difficult to accurately capture the critical turning point from latent to deteriorating failures, ultimately leading to misjudgment of critical turning points and failing to provide timely and accurate triggering basis for subsequent failure prediction. Summary of the Invention

[0005] To address the shortcomings of existing technologies, this invention provides a method and apparatus for predicting oil equipment failures based on big data, thereby resolving the problems existing in the background technology.

[0006] To achieve the above objectives, the present invention provides the following technical solution: a method and apparatus for predicting oil equipment failures based on big data, comprising the following steps: Step S1: Collect raw spatiotemporal data of each oil equipment and preprocess it to obtain preprocessed spatiotemporal data of the oil equipment; Step S2: Use the variable time-delay transfer entropy algorithm to perform causal analysis on the preprocessed spatiotemporal data of the oil equipment, obtain the causal relationship of the oil equipment, and generate the spatiotemporal behavior association chain of the oil equipment; Step S3: Use the sliding window algorithm to perform local feature analysis on the spatiotemporal behavior association chain to obtain the failure evolution trend value of the oil equipment. Integrate the failure evolution trend value with the timestamps in the spatiotemporal behavior association chain to generate the failure evolution sequence of the oil equipment. Step S4: Use the Bayesian change point detection algorithm to identify abrupt change points in the fault evolution sequence to obtain the critical turning points of the oil equipment failure; Step S5: Construct a fault prediction model based on a long short-term memory network, input the critical fault inflection point and fault evolution sequence into the fault prediction model, and output the fault prediction results across oil equipment.

[0007] Preferably, the acquisition of raw spatiotemporal data from each oil equipment includes: By deploying a sensor network at key nodes in oil production, raw spatiotemporal data reflecting the operating status of equipment are collected. This data can be formally represented as a raw dataset. :

[0008] in, It is the original dataset; It is an identifier that represents the spatial location and device type of the data source; It is the precise timestamp of data collection; It is a vector of multidimensional physical parameters that is strongly correlated with the device type and is measured at that time and location.

[0009] Preferably, the preprocessing of the original spatiotemporal data to obtain preprocessed spatiotemporal data of oil equipment includes: Preprocessing is performed on the raw data in the original dataset to improve its quality and consistency. The preprocessing operations mainly include denoising, normalization and missing value imputation. In the denoising stage, the Kalman filter algorithm is used to suppress measurement noise and extract the true signal trend; In the normalization process, in order to eliminate the differences between different physical dimensions and numerical ranges, the minimum-maximum normalization method is used to scale the data to a uniform range. For missing value imputation, linear interpolation is used for estimation, taking into account the characteristics of time series data.

[0010] Finally, spatial alignment is performed. The primary task of spatial alignment is to construct a spatial topology model of the entire oil equipment, G=(V,E), where V is the vertex set, representing each equipment node, and E is the edge set, representing the physical connection relationship between equipment and the direction of medium flow. Based on spatial topology, time alignment requires precise compensation for the time delay of signal propagation between devices. To improve the accuracy and efficiency of time delay estimation, a "data-driven approach guided by a physical model" is adopted: Based on the different physical nature of signals, the propagation mechanism and velocity value rules are clearly distinguished, thereby providing accurate prior constraints for data-driven computation and estimating the theoretical time delay range according to rules. , ],in and These are the minimum and maximum time delay estimates, respectively; within the theoretical time delay range [ , Under the physical constraints of By finding the cross-correlation function Time delay estimation to achieve the maximum value The value can be used to obtain the optimal time delay estimate from device i to device j. Ultimately, this extends the optimal time delay estimation between the two devices to global time synchronization across the entire device cluster.

[0011] After the above processing, preprocessed spatiotemporal data of petroleum equipment are obtained. .

[0012] Preferably, the step of using the variable time-delay transfer entropy algorithm to perform causal analysis on the preprocessed spatiotemporal data of the oil equipment to obtain the causal relationships of the oil equipment includes: Before calculating the variable time-delay transfer entropy, a physical connection strength threshold λ is set; only when the connection weight... Only when the value is greater than λ will the device pair be retained and subsequent causal strength quantification and direction determination be performed; Subsequently, the variable time-delay transfer entropy is calculated to quantify the causal strength between devices; the basic formula for variable time-delay transfer entropy is defined as:

[0013] in, This represents the variable time-delay transfer entropy value from device X to device Y; It is conditional entropy; It is the state value of device Y at time t; This represents the set of historical states of device Y over the past k time steps; It is the set of past states of device X after time delay adjustment.

[0014] After calculating the variable time-delay transfer entropy, a normalization process is performed, and the formula is as follows:

[0015] in, It is the ratio of transfer entropy; It is the variable time-delay transfer entropy from X to Y; It is the variable time-delay transfer entropy from Y to X.

[0016] Preferably, the spatiotemporal behavioral correlation chain of the generated oil equipment includes: Based on the calculated causal relationships, the next step is to generate spatiotemporal behavioral correlation chains, and the specific construction steps are as follows: First, from all the calculated causal relationships, identify device pairs with a significant causal direction, i.e., those that satisfy... The "cause-effect" combination is used as a candidate link; Subsequently, to address the issue of unclear priority in association chain construction when multiple causal paths may point to the same device, a two-level priority rule was developed to filter and sort causal links. After determining the main link, based on the spatial topology model G=(V,E) and the variable time delay parameters... The selected events are time-aligned; Finally, the cause-effect events that have completed time alignment are linked sequentially to form a chain sequence; The resulting spatiotemporal behavioral correlation chain serves as a structured framework, clarifying the causal order and propagation path of events within the device group.

[0017] Preferably, the step of using the sliding window algorithm to perform local feature analysis on the spatiotemporal behavioral correlation chain to obtain the failure evolution trend value of the oil equipment includes: Transforming spatiotemporal behavioral correlation chains into a comprehensive evolutionary index sequence representing health status. Incorporating physical connection weights The optimized evolution index value is obtained by weighted summation of all significant causal events at time t. :

[0018] in, C is the evolution index value at time t; C represents the value that satisfies the following at time t. The set of significant causal relationships (X→Y)>1; It is the normalized ratio of transfer entropy.

[0019] Next, the evolution index sequence { The sliding window is used for partitioning, and the window size is selected based on the time-scale characteristics of equipment failure evolution. Its mathematical expression is:

[0020] in, It is the k-th sliding window; k is the window index; It is the evolutionary index sequence at time 10:00. data points It is the time index; w is the window size. In each sliding window Within this process, multi-dimensional feature extraction is performed. The feature extraction process is represented as follows:

[0021] in, It is the feature vector of the k-th window; It is an intensity characteristic; Clustering characteristics It is a diffusion characteristic; Subsequently, the fault evolution trend value is calculated. This fault evolution trend value is used to smooth the feature sequence using a weighted moving average method, and its mathematical expression is as follows:

[0022] in, is the fault evolution trend value at time point t; t is the time series index; j is the summation index; h is the window size for trend smoothing; It is an exponentially decaying weight; It is the weighted Euclidean norm; It is an eigenvector.

[0023] Preferably, the step of integrating the fault evolution trend value with the timestamps in the spatiotemporal behavior association chain to generate the fault evolution sequence of the oil equipment includes: Based on the timestamps corresponding to each node in the spatiotemporal behavior association chain and its implied multidimensional feature vectors And its corresponding fault evolution trend value By performing correlation and integration, a structured fault evolution sequence E is generated, synchronously correlated with a set of key contributing devices. This key contributing device is for the collection The screening follows a dual criterion: only device pairs that make significant contributions and have strong physical correlations are selected. Therefore, the fault evolution sequence sequence The mathematical representation of is optimized as follows:

[0024] in, It's a timestamp; It is the fault evolution trend value; It is the multidimensional feature vector at that moment; It is the selected core equipment set

[0025] Preferably, the step of using the Bayesian change point detection algorithm to identify abrupt change points in the fault evolution sequence to obtain the critical inflection points of the oil equipment failure includes: To achieve collaborative analysis of multi-dimensional features, fault evolution trend values ​​are... With feature vectors splicing together to form a four-dimensional observation sequence }; Subsequently, an improved multidimensional local outlier factor algorithm was used to analyze the four-dimensional observation sequence { Anomaly detection is performed; subsequently, based on the purified four-dimensional observation sequence... Perform Bayesian multivariable detection, assuming the sequence has K variable points. < <⋯< The sequence is divided into K+1 segments, and the data generation model for each segment m is as follows:

[0026] in, They are polynomial basis vectors; It is a 4×3 coefficient matrix; It is observation noise; the optimal number of change points K is determined by minimizing the Bayesian information criterion; Ultimately, each identified critical inflection point of failure The output is enhanced into a structured triple:

[0027] in, It is the location of the critical turning point of the fault; This represents the posterior probability of the critical inflection point of the fault. It is the critical turning point of the failure. A set of key contributing devices associated with each other.

[0028] Preferably, the fault prediction model constructed based on long short-term memory networks, which inputs the critical inflection points and fault evolution sequences into the fault prediction model and outputs fault prediction results across oil equipment, includes: Based on the fault evolution sequence E, a time-series framework is established, identifying key inflection points. The feature vectors are inserted at the corresponding timestamp positions to form the model input; Long Short-Term Memory (LSTM) is a special type of recurrent neural network that effectively captures long-term dependencies in time series. Its core lies in selectively memorizing and forgetting information through a gating mechanism. During the training phase, historical fault data is used to supervise the learning of the LSTM network, and the weighted loss function is dynamically selected according to the prediction target. During the prediction phase, the features of the critical inflection points and failure evolution sequences monitored in real time are input into the trained LSTM model, which outputs the failure probability distribution of each device over a future period:

[0029] Wherein, P( () represents the failure probability distribution over the next T time steps; the softmax function transforms the network output into a probability distribution; and These are the weights and bias parameters of the output layer; It is the hidden state at the current time step; Finally, a cross-equipment fault prediction result is generated. This cross-equipment fault prediction result is a structured and quantifiable decision support output. A cross-equipment association regularization term is introduced to constrain the learning process, and the total loss function of model training is optimized as follows:

[0030] in, The core classification loss is λ = 0.01, the regularization coefficient is C, and C is the set of significant causal relationships. It is the ratio of transfer entropy; It is the physical connection weight; and It predicts the probability of failure; Based on this, the raw single-device failure probability output by the model will be further processed to generate actionable information containing the following dimensions: device location, which clearly indicates the ID of the source device with the highest failure risk and its predicted specific failure type; time window, which predicts the time interval during which the high-risk state will occur or last; and probability intensity, which is the risk level quantified by probability value. Ultimately, all the information was integrated into a structured instruction that the project could use directly.

[0031] An artificial intelligence-based power distribution box regulating device includes a memory, a processor, and a computer program stored in the memory, wherein the processor executes the computer program to implement the steps of the above method.

[0032] Beneficial effects This invention provides a method for predicting oilfield equipment failures based on big data, involving machine learning and deep learning technologies, which has the following beneficial effects: (1) This big data-based method for predicting oil equipment failures employs a variable time-delay transfer entropy algorithm based on piecewise aggregation approximation. This algorithm not only effectively captures the nonlinear and time-varying causal mechanisms between equipment, but also quantifies the causal strength through the transfer entropy ratio, accurately mining potential cross-equipment relationships. Furthermore, this method combines pre-processing physical filtering with a weighted aggregation strategy, ensuring that the generated spatiotemporal behavioral association chain retains the algorithm's advantages in mining complex temporal relationships while strengthening its engineering physical logic foundation. Compared to traditional fixed-delay methods, it significantly improves the accuracy and reliability of causal discovery.

[0033] (2) This big data-based method for predicting oil equipment failures employs a Bayesian change point detection algorithm and constructs a piecewise regression framework based on a probability model. It eliminates the need to pre-determine the number of change points and can adaptively capture abrupt changes such as trend jumps based on sequence characteristics. This method extends the algorithm into a detection framework that integrates multi-dimensional features and equipment-related priors. This significantly enhances the accuracy and anti-interference capability of identifying key turning points in complex failure evolution patterns while retaining the advantages of adaptability and uncertainty quantification, thus realizing a shift from passive alarm to proactive prediction.

[0034] (3) This big data-based petroleum equipment fault prediction method incorporates the Local Outlier Factor (LOF) algorithm before Bayesian change point detection. By calculating the local outlier factor value, it identifies and removes local anomalies in the sequence, effectively suppressing the interference of transient noise. This method extends the LOF algorithm to a multi-dimensional feature space for collaborative anomaly detection and deeply integrates it with subsequent detection processes. This fusion design not only more accurately purifies the sequence but also ensures that Bayesian change point detection can focus on analyzing real fault evolution and abrupt changes, thereby systematically enhancing the robustness of the overall inflection point identification process and the stability of the results. Attached Figure Description

[0035] Figure 1 This invention proposes a method for predicting petroleum equipment failures based on big data.

[0036] Figure 2 The present invention provides a hierarchical graph of spatiotemporal behavioral correlation chains for a big data-based method for predicting oil equipment failures.

[0037] Figure 3 This invention presents a hierarchical diagram of cross-oil equipment fault prediction results obtained through a big data-based oil equipment fault prediction method. Detailed Implementation

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

[0039] Please see Figures 1-3 This invention provides a technical solution: a method for predicting oil equipment failures based on big data. Specifically, the method for predicting oil equipment failures based on big data is provided below; please refer to [link / reference]. Figure 1 The method includes the following steps: Step S1: Collect the raw spatiotemporal data of each oil equipment and preprocess it to obtain preprocessed spatiotemporal data of the oil equipment.

[0040] The goal of this step is to acquire high-quality, spatiotemporally consistent data on oil equipment, laying the foundation for subsequent causal analysis and fault prediction. The process mainly includes three core stages: data acquisition, data preprocessing, and spatiotemporal alignment.

[0041] First, data acquisition is conducted. A sensor network deployed at key nodes in oil production (such as oil wells, pump units, compressors, and pipeline valves) collects raw spatiotemporal data reflecting the equipment's operating status. To overcome the limitation of general parameter lists in accurately characterizing the unique fault mechanisms of various equipment types, this method abandons homogeneous parameter acquisition approaches and instead classifies and collects core monitoring indicators based on the equipment's core operating mechanisms and fault mechanisms. This data can be formally represented as a raw dataset. :

[0042] in, It is the original dataset; It is an identifier (such as device ID) that represents the spatial location and device type of the data source. It is the precise timestamp of data collection; It is a vector of multidimensional physical parameters that are strongly correlated with the equipment type and are measured at that moment and location. Specifically, for oil wells, the core parameters include sucker rod vibration acceleration, wellhead temperature, bottom hole pressure, crude oil flow rate, sucker rod load, and pump efficiency; for pump units, the core parameters include motor vibration velocity, bearing temperature, inlet and outlet pressure, medium flow rate, motor speed, and lubricating oil pressure; for compressors, the core parameters include cylinder vibration, exhaust temperature, intake pressure, exhaust pressure, exhaust volume, and crankcase temperature; for pipeline valves, the core parameters include pressure difference across the valve, medium flow rate, valve opening, valve body vibration, sealing surface temperature, and leakage.

[0043] The parameter vector has a unified format as follows: =( , ,..., The dimension m takes values ​​from 5 to 7 depending on the specific equipment type, and each dimension strictly corresponds to the aforementioned unique core parameter and its standard dimension (e.g., vibration acceleration: mm / s²). 2 Pressure: MPa).

[0044] Then, the raw data in the original dataset is preprocessed to improve its quality and consistency. The preprocessing operations mainly include denoising, normalization, and missing value imputation.

[0045] In the denoising stage, the Kalman filter algorithm is used to suppress measurement noise and extract the true signal trend. The core of this algorithm is to construct a closed-loop process that includes prediction and correction steps, as follows: First, define the core parameters, specifically including: the state vector. It is composed of the core monitoring parameters of the equipment at time t; the state transition matrix is ​​taken as... (m-order identity matrix) to reflect the smooth transition characteristics of parameters during normal operation of oil equipment; control variables The impact of equipment operating load (such as motor power or crude oil production) is modeled by a control input matrix B (an m×1 vector) obtained by fitting historical data; the observation noise covariance matrix R is set according to the sensor accuracy.

[0046] The complete algorithm flow begins with state prediction, as shown in the following formula:

[0047] in, It is the prior estimate of the state at time t.

[0048] Subsequently, the observation equation is introduced, and the specific formula is as follows:

[0049] in, These are sensor observations; H= ; It is observation noise that follows N(0,R).

[0050] The key correction step of the algorithm involves calculating the Kalman gain. (Initial covariance matrix) =0.1 This is implemented and the state is updated accordingly. Covariance Update Through this iterative process of "first predicting the trend, then correcting for deviations with actual observations," the algorithm can effectively separate noise and ultimately output a sequence { This is the real signal after noise reduction.

[0051] In the normalization process, to eliminate differences in physical dimensions and numerical ranges, a minimum-maximum normalization method is used to scale the data to a uniform range. To prevent extreme outliers in the original data (such as sensor fault peaks) from distorting the normalization boundary and causing data distribution distortion, this method adds a rigorous outlier detection and processing step before normalization.

[0052] First, for each monitoring parameter Calculate its mean across the entire dataset. and standard deviation And outlier detection is performed based on the 3σ criterion: if If the value is greater than a certain threshold, the data point is considered an outlier. Then, the "boundary replacement method" is used to handle outliers, specifically those values ​​greater than a certain threshold. +3 Replace outliers with this upper bound, and replace values ​​smaller than this upper bound. -3 Outliers are replaced with this lower bound. After this processing, a new minimum value for each parameter is calculated based on the corrected clean dataset. With the maximum value And perform normalization calculations:

[0053] in, This is the normalized result; d is the original value of a certain feature. To ensure that the normalized data is evenly distributed and can effectively reflect the parameter fluctuation characteristics, we define the coefficient of variation of the final result as ( ). / The value must be no less than 0.1. If this condition is not met, the outlier handling logic needs to be re-examined.

[0054] For missing value imputation, linear interpolation can be used to estimate the value, taking into account the characteristics of time series data.

[0055] in, It represents the estimated value of the missing data at time point t; t is the time point where the missing value is located. and It is the timestamp of the nearest valid data point before and after the missing point, satisfying... <t< ; and Is and The actual measured value at that moment.

[0056] Finally, spatiotemporal alignment is performed. The core purpose is to solve the problem of data asynchrony caused by the physical distribution of devices and signal propagation delays, and to unify multi-source data into a spatiotemporal framework that can truly reflect physical interactions.

[0057] The primary task of spatial alignment is to construct a spatial topology model of the entire oilfield equipment, G=(V,E), where V is the vertex set representing each equipment node, and E is the edge set, quantitatively characterizing the strength of physical connections between equipment. To address the issue that the edge set in the original model only qualitatively describes connectivity and cannot provide physical constraints for subsequent causal analysis, this method assigns a connection weight to each edge. Its calculation comprehensively considers the connection type and key physical parameters.

[0058] Specifically, for pipe connections, their connection weights The value is determined by both the pipe length and diameter, and the calculation formula is as follows:

[0059] in, Pipe length (unit: m); Where α is the pipe diameter (in meters); α=0.01 and β=0.5 are engineering experience coefficients, with weight values ​​constrained to the range of 0.05 to 1.0.

[0060] For drive shaft connections, the weight is proportional to the rotational speed; the connection weight... The calculation formula is:

[0061] in, It is the transmission shaft speed (unit: r / min), with a weight range between 0.1 and 1.0.

[0062] For electrical connections (such as motor-pump assembly), connection weight The power is determined by the rated power of the equipment, and the calculation formula is as follows:

[0063] in, It is the rated power (unit: kW), with a weight range between 0.08 and 1.0.

[0064] All connection weights are stored in the adjacency matrix W, i.e. = If the two devices are not directly connected, the corresponding matrix element is 0. In the subsequent causal analysis in step S2, this weight will be incorporated into the calculation as a physical prior. Specifically, when calculating the variable time-delay transfer entropy, the original value will be used. With connection weights Multiplication, that is This effectively weakens the causal correlation weight between weakly connected devices and avoids misjudgment.

[0065] Based on spatial topology, time alignment requires precise compensation for the time delay of signal propagation between devices. To improve the accuracy and efficiency of time delay estimation, a "data-driven approach guided by a physical model" is adopted. First, this method clearly distinguishes between propagation mechanisms and velocity rules based on the different physical nature of the signals, thus providing precise prior constraints for data-driven calculations. Specifically, based on the spatial topology between devices (such as pipe length L) and the signal type, the theoretical time delay range is estimated according to the following rules: , ]: For transient signals (such as pressure waves and vibration waves), their propagation occurs in the form of waves at extremely high speeds. The propagation speed v=c, which is the speed of sound in the medium, and its value is determined by the type of medium and operating conditions (temperature, pressure). For example, crude oil at 20°C has c≈1300-1500 m / s, and natural gas has c≈300-400 m / s.

[0066] For steady-state signals (such as temperature conduction or changes in medium flow rate), their propagation depends on the macroscopic flow of the medium itself, and the speed is relatively slow. Propagation speed v = The average velocity of the medium is calculated from the pipe flow rate and cross-sectional area.

[0067] Where Q is the volumetric flow rate (unit: m3 / h); D is the pipe inner diameter (unit: m).

[0068] Determine the propagation speed v of the corresponding signal type and its possible fluctuation range. Then, the theoretical time delay boundary can be calculated:

[0069]

[0070] in, and yes of ; and This is the theoretical time delay range. This range will be used to search for the optimal time delay during subsequent cross-correlation function calculations. The physical constraint range.

[0071] Subsequently, within the theoretical time delay range [ , Under the physical constraints of [the original cross-correlation function], the optimal time shift is found by calculating the cross-correlation of the signal sequence. The results of previous methods are subject to interference from the absolute amplitude of the signal (such as temperature signals being much larger than vibration signals) and cannot truly reflect the timing synchronization relationship. This method uses a normalized cross-correlation function for calculation, and its formula is as follows:

[0072] in, It is the normalized cross-correlation function value of the signal of device i and the signal of device j when the time delay is τ. The larger the value, the higher the similarity between the two signals under the time delay τ. It is the signal value collected from device i (usually a suspected fault source or upstream device) at time point t; At a certain point in time The signal value collected from device j (usually the affected device or downstream device); It is a time index; i and j are device indexes; and These are the mean values ​​of the signal sequences from device i and device j, respectively. This function eliminates the DC component by subtracting the mean and normalizes the signal by dividing by the product of the standard deviations of the two signals, making... The value of is strictly limited to [-1, 1]. The closer its value is to 1, the stronger the waveform similarity and time synchronization of the two signals under the time delay τ. Optimal time delay By searching within the constraint interval The maximum value τ is used to determine this. To ensure the effectiveness of the time alignment, a threshold is set for verification: the final optimal time delay must satisfy... ≥0.7. If it is below this threshold, the signal synchronization between the devices is deemed insufficient, and it will not be considered in the causal analysis of the subsequent step S2.

[0073] By finding the cross-correlation function To obtain the maximum value The optimal time delay from device i to device j can be obtained by finding the value of the normalized cross-correlation function. By obtaining the maximum value of τ, the optimal time delay estimate from device i to device j can be obtained. The fundamental principle lies in the fact that the core function of the cross-correlation function is to quantify the similarity between two signals at different time shifts. Its maximum value corresponds to the state where the waveforms of the signals from device i and device j have the highest degree of overlap. Since signal propagation between devices involves a time delay determined by physical distance and medium characteristics, only when the time offset equals this actual propagation time delay can the causal signal (such as pressure fluctuations in device i) and the resulting signal (such as vibration response in device j) achieve precise alignment on the time axis. Therefore, maximizing the cross-correlation function to its maximum value τ essentially marks this most matching and cooperative alignment point, which is the optimal time delay estimate. The signal propagation of oil equipment follows certain physical laws, and this time delay is regular. By traversing and calculating and finding the peak value of the cross-correlation function, this real propagation delay can be effectively captured.

[0074] Finally, to extend the optimal time delay estimation between two devices to global time synchronization of the entire equipment cluster, this method implements a four-step process: "benchmark unification - time delay propagation - compensation alignment - verification and correction". First, a key hub device (such as a long-distance pipeline starting pump station or core controller) is selected from the petroleum equipment cluster as the global time benchmark, and its original timestamp is directly determined as the benchmark time. Next, based on the constructed spatial topology model, starting from the benchmark device, the total time delay of each device in the cluster relative to the benchmark is calculated through path accumulation and propagation. Subsequently, the original timestamp of each device is uniformly subtracted from its corresponding total time delay, achieving mapping and compensation alignment of all data to the benchmark time axis. To ensure the accuracy of synchronization, a consistency verification is finally performed: selecting clearly defined time points for physical events such as equipment start-up and shutdown, and valve switching, the standard deviation of the compensated timestamps of multiple devices is verified (required to be ≤0.1 seconds); if the deviation exceeds the standard, the optimal time delay estimate on the relevant path is recalculated and corrected. Through the above systematic process, it is ensured that multiple device data points representing the same physical event are strictly synchronized on the time axis.

[0075] After the above processing, preprocessed spatiotemporal data of petroleum equipment are obtained. It accurately depicts the collaborative operation status of the equipment group, providing reliable input for the causal analysis in the subsequent step S2.

[0076] Step S2: Use the variable time-delay transfer entropy algorithm to perform causal analysis on the preprocessed petroleum equipment spatiotemporal data, obtain the causal relationship of the petroleum equipment, and generate the spatiotemporal behavior association chain of the petroleum equipment.

[0077] The core of this step lies in using a variable time-delay transfer entropy algorithm based on piecewise aggregation approximation to preprocess the spatiotemporal data of oil equipment. Causal analysis is performed to uncover the causal relationships between equipment, and a spatiotemporal behavioral correlation chain describing the overall evolution trend of the equipment group is generated. This step uses the preprocessed spatiotemporal data of oil equipment output from step S1. Based on this, it is important to clarify that, in order to fairly compare the causal strength between parameter sequences originating from different devices and having different physical meanings, this method does not directly compare the original parameter values, but instead maps them to a unified fault-sensitive feature space.

[0078] Specifically, before calculating the transfer entropy from device X to device Y, the normalized core parameter sequences of devices X and Y are first transformed into a set of statistical feature vectors of the same dimension, for example:

[0079]

[0080] in, It is the intensity, which can be characterized by the root mean square (RMS) value of the sequence within a sliding window, reflecting the average energy level of the signal; It is volatility, which can be characterized by standard deviation and reflects the degree of dispersion of the signal; The trend slope, which can be characterized by the coefficients of the linear fit, reflects the local change trend of the signal. These features collectively characterize the evolution pattern of the equipment state and are strongly correlated with the fault mechanism.

[0081] Based on this, the variable time-delay transfer entropy The calculation is performed between the eigenvector sequences. The sequences processed by the algorithm... and , respectively, the historical sequences of fault feature vectors of equipment X and equipment Y { }and{ This approach ensures that comparisons of causal strength are performed at a level of general and comparable abstract features that reflect the evolution of device states, thereby significantly improving the reliability, robustness, and engineering interpretability of cross-device causal discovery.

[0082] The entire analysis process emphasizes the variability and nonlinearity of causal direction. The final output of the correlation chain is not a simple stacking of events, but a chain structure with a clear causal order, which is reconstructed from related events such as vibration anomaly, temperature rise, pressure fluctuation, and flow disturbance based on engineering physical logic (such as pressure transmission path or vibration propagation mechanism), providing structured input for fault prediction.

[0083] First, core parameters are initialized. To adapt to the physical characteristics and interaction modes of different device types, this method adopts a strategy of "unified underlying algorithm and categorized configuration of core physical parameters." The values ​​of key parameters are not fixed, but are dynamically loaded from a predefined parameter configuration library based on the type of the device pair. This configuration library is built based on historical data and engineering experience and mainly includes: The physical connection strength threshold λ is used for pre-filtering. For example, for pipe connections with significant signal attenuation, a stricter threshold (e.g., λ=0.08) is set to retain more potential causal paths; for drive shafts or electrical connections with direct signal transmission, a relatively lenient threshold (e.g., λ=0.12) is set. Theoretical time delay range [ , ] is used to constrain the time delay estimation in causal discovery, and its value strictly follows the theoretical time delay range determined in step S1 based on the physical nature of the signal (transient wave or steady flow) and the propagation speed (sound speed or medium flow speed).

[0084] Before calculating the variable time-delay transfer entropy, a pre-filter is performed. The logic is as follows: determine the connection weight between any two devices i and j. If the connection weight If the physical connection strength threshold λ is less than or equal to the physical connection strength threshold λ, then it is directly determined that there is no physical basis to support a significant causal relationship between the two, and subsequent calculations are skipped; only if the connection weight is less than or equal to the physical connection strength threshold λ is the subsequent calculations skipped. Only when the physical connection strength threshold λ is reached is the device pair retained and subsequent causal strength quantification and direction determination performed. This process ensures that causal analysis focuses only on device pairs with a strong physical connection foundation, fundamentally eliminating spurious causal associations caused by weak connections or noise.

[0085] Subsequently, variable time-delay transfer entropy is calculated to quantify the causal strength between devices. Variable time-delay transfer entropy is an extension of traditional transfer entropy; it considers the variability of time delay in causal interactions and is more suitable for scenarios where signal transmission delays are not fixed under complex operating conditions in oil equipment. Its basic formula is defined as:

[0086] in, It represents the variable time-delay transfer entropy value from device X to device Y, which measures the reduction in uncertainty of the current state of Y by introducing past information of device X, given that the past state of device Y itself is known. It is the conditional entropy, used to calculate the uncertainty of a random variable under given conditions; It is the state value of device Y at time t; This represents the set of historical states of device Y over the past k time steps, i.e., { , ,…,}, where parameter k is the embedding dimension of the historical state of device Y. To accurately reconstruct the state space of device Y, the pseudo nearest neighbor method is used to independently calculate and determine it: for the state sequence of device Y, starting from the embedding dimension d=1, the pseudo nearest neighbor ratio FNN(d) under each dimension is gradually increased. When FNN(d) is lower than the preset threshold (e.g., 5%) for the first time, the corresponding dimension d is determined to be k. It is the set of past states of device X after time delay adjustment, defined as , where l is the embedding dimension of the historical state of device X, which is also determined independently based on the sequence of device X using the pseudo-nearest neighbor method. It is a variable time delay parameter, and its value is constrained within the theoretical time delay range. , In this study, the dynamic time warping algorithm is used to estimate the time from the data, and the specific rules are as follows: First, the optimal matching path between sequences X and Y is obtained by minimizing the cumulative distance D(p,q). ,in , These are the time indices for X and Y, respectively; subsequently, for each time step t (corresponding to the time index of X)... ), find the time index of the matching Y in path P. Then the time delay at this moment is calculated as follows: = - This ensures the consistency of the time delay and signal matching relationship. The calculation of conditional entropy relies on the estimation of the joint probability distribution. Given that oil equipment data are continuous variables, this method preferentially uses kernel density estimation to better capture nonlinear dependencies, with a Gaussian kernel selected as the kernel function. Its bandwidth h is determined by cross-validation: Define the log-likelihood loss function. ,in The kernel density estimate is used as the final selection of h that minimizes L(h) as the optimal bandwidth.

[0087] After calculating the variable time-delay transfer entropy, it is necessary to further determine the validity of the causal direction. Since the transfer entropy value may be affected by data noise or random fluctuations, directly using the original value may lead to misjudgment. Therefore, a transfer entropy ratio is introduced for normalization to enhance the robustness of the causal direction. The formula is as follows:

[0088] in, It is the transfer entropy ratio, used to determine the causal relationship between two devices. It is a ratio that determines the causal direction from X to Y. The causal direction is determined by comparing the magnitude of the transfer entropy in the two directions from X to Y and Y to X. If the ratio is greater than 1, then X is considered to be the cause variable of Y, that is, the two have a causal relationship. The larger the ratio, the stronger the causal relationship. It is the variable time-delay transfer entropy from X to Y; It is the variable time-delay transfer entropy from Y to X (reverse causality). In the ratio calculation, both the denominator and numerator are based on the same probability estimation method to avoid bias, and the estimation error is reduced by multiple sampling or bootstrapping.

[0089] Based on the calculated causal relationships, the next step is to generate spatiotemporal behavioral correlation chains. This process aims to extract the main fault propagation sequences with the highest confidence from the complex device causal network. The specific construction steps are as follows: First, from all the calculated causal relationships, identify device pairs with a significant causal direction, i.e., those that satisfy... The "cause-effect" combination is used as a candidate link.

[0090] Subsequently, to address the issue of unclear priority in association chain construction when multiple causal paths may point to the same device (path conflict), this method establishes a clear two-level priority rule to filter and sort causal links: First-level priority (causal strength) is based on the ratio of transition entropy. The magnitude of the ratio is used to sort all candidate causal relationships in descending order. The larger the ratio, the stronger the information transmission and causal driving effect from X to Y, and the higher the priority of this link in reflecting the main propagation path of the fault.

[0091] Secondary priority (physical connection strength) refers to the ratio of the transfer entropy of two causal paths pointing to the same resulting device being very close (the difference Δ) When the value is ≤0.2), then the connection weights in the spatial topology model constructed in step S1 are used as the basis. Make a decision. Prioritize connection weights. A path with a larger value is preferred because it represents a stronger physical connection between devices, which better aligns with propagation priorities in practical engineering.

[0092] After determining the main link using the above rules, based on the spatial topology model G=(V,E) already constructed in S1 and the estimated variable time delay parameters... The selected events are time-aligned. For example, pressure fluctuations at upstream pumping stations are aligned with vibration responses at downstream valves by using dynamic time-bending distances, compensating for variable time delays and ensuring they are synchronized on the time axis.

[0093] Finally, the time-aligned cause-effect events are linked sequentially to form a chain sequence such as "abnormal vibration → increased temperature rise → pressure fluctuation," clearly depicting the spatiotemporal propagation trajectory of the fault. For causal relationships that are not included in the main chain due to priority rules but are still significant (such as secondary paths in path conflicts), they are marked as supplementary associations, thereby ensuring the clarity and reliability of the main chain while preserving the integrity of the causal network.

[0094] The resulting spatiotemporal behavioral correlation chain serves as a structured framework, clarifying the causal sequence and propagation path of events within the equipment group, thus providing a foundation for the subsequent fault evolution analysis in step S3. Step S2, by combining variable time delay processing and information theory metrics, effectively improves the accuracy and engineering practicality of causal relationship discovery in non-stationary, nonlinear oilfield equipment.

[0095] Step S3: Use the sliding window algorithm to perform local feature analysis on the spatiotemporal behavior association chain to obtain the failure evolution trend value of the oil equipment. Integrate the failure evolution trend value with the timestamps in the spatiotemporal behavior association chain to generate the failure evolution sequence of the oil equipment.

[0096] This step uses a spatiotemporal behavioral correlation chain as input to construct a fault evolution sequence for oil equipment that reflects the transition from an early latent state to a clearly abnormal state and then to a high-risk stage. The construction of this sequence follows the evolutionary pattern of "from weak to strong, from scattered to concentrated, and from local to global." Its focus is not on identifying anomalies at a certain instant, but on identifying the subtle trend changes that the equipment exhibits over time.

[0097] First, the spatiotemporal behavioral correlation chain is transformed into an evolutionary index sequence that can comprehensively characterize health status. To eliminate interference between different devices To address the differences in the dynamic range of values, the transition entropy ratio of each causal event is first normalized, mapping it to a uniform [0,1] interval. The normalization formula is:

[0098] in, and These are the minimum and maximum values ​​of the transfer entropy ratio, respectively, obtained from fault-free historical data over the past 6 months. If sufficient historical data is unavailable, empirical values ​​can be used. =1.0, =5.0.

[0099] Based on normalization, physical connection weights that characterize the tightness of physical connections between devices are incorporated. (Given by the spatial topology model in step S1, with values ​​ranging from 0.05 to 1.0), the optimized evolutionary index value is obtained by weighted summation of all significant causal events at time t. :

[0100] in, It is the evolution index value at time t, and is a comprehensive index characterizing the intensity of overall abnormal causal activities; This means that at time t, the following condition is met. The set of significant causal relationships; This is the normalized transition entropy ratio calculated in step S2, representing the causal strength from device X to Y. This evolutionary index value... The larger the value, the more intense the anomalous causal activity at time t. Weighted summation ensures stronger physical connections (physical connection weights). Larger causal paths contribute more significantly to evolutionary indices, increasing their value. It is more in line with the actual laws of physical propagation.

[0101] It is important to note that all subsequent feature analysis will be based on the evolutionary index sequence { The numerical variation pattern of}, and no longer directly dependent on the composition A single causal event (X→Y) and its intensity .

[0102] Next, the evolution index sequence { A sliding window approach is used to divide the continuous time series data into multiple overlapping windows for local feature analysis. The window size is chosen based on the time-scale characteristics of equipment failure evolution, and its mathematical expression is:

[0103] in, It is the k-th sliding window; k is the window index; It is the evolutionary index sequence at time 10:00. data points This is the time index; w is the window size, i.e., the number of data points contained in each window, which can be obtained as follows: First, calculations are performed based on the inherent timescale of fault evolution: base window size. The duration of the fault latency phase is determined by both the average duration of the fault latency phase and the data sampling frequency, and the calculation formula is as follows:

[0104] in, It is the average duration of the fault latency period (in seconds) obtained from historical data. This refers to the data sampling frequency (unit: Hertz). For example, the average latency from the onset of a certain type of pump failure to its significant deterioration. =30s, data sampling frequency =1Hz, then the base window size =30 data points. Subsequently, to find the optimal analytical granularity, based on the aforementioned... Based on [0.5], cross-validation was used in the [0.5] range. ,2 Optimization is performed within this range: Different window sizes are iterated within this range, and the correlation coefficient between the "fault evolution trend value" sequence extracted from each window and the corresponding "actual fault stage" label sequence is calculated. Finally, the window size that maximizes the correlation coefficient is selected as the optimal value. Furthermore, to balance computational efficiency and trend continuity, and to avoid information breaks between windows, this method fixes the window overlap rate at 50%, i.e., the window sliding step size is step = 0.5 × Through the above steps, the window size... The determination is based on objective processes of data and quantitative verification, which ensures the standardization of the analysis process and the reproducibility of the results.

[0105] In each sliding window Within this framework, multi-dimensional feature extraction is performed. These features are specifically designed to capture the evolutionary patterns of "from weak to strong, from scattered to clustered, and from local to global." The feature extraction process can be represented as:

[0106] in, This is the feature vector of the k-th window; the feature function is defined as follows, mainly including intensity features. Clustering characteristics and diffusion characteristics :

[0107] in, It is an intensity characteristic (from weak to strong) used to measure the average intensity and growth trend of abnormal activity within a window, ensuring that intensity changes are synchronized with the physical connection basis of critical equipment; Indicates at time All significant causal events In the context, the maximum physical connection weight between corresponding devices. . This ensures that events with stronger physical connections contribute more to the evolutionary indices, thus linking intensity characteristics to the actual impact capabilities between devices. It is the trend progression ratio, reflecting the ratio of the average intensity of the current window to the average intensity of the previous adjacent window (ε= (Used to avoid division by zero). When the ratio is greater than 1, it indicates that the anomaly intensity is increasing, which intuitively matches the evolutionary pattern of "from weak to strong".

[0108]

[0109] in, It is a clustering feature (from scattered to clustered), used to quantify the degree of concentration of anomalous events over time and their dynamic convergence process; and These are windows Internal evolution indicators The ratio of the standard deviation and the mean of . It measures the dispersion of data; the smaller the value, the more concentrated the data. This is the clustering change ratio, which dynamically captures changes in central tendency by comparing the coefficient of variation of the previous window and the current window: if the ratio is greater than 1, it indicates that the dispersion of the current window has decreased compared to the previous window, meaning a convergence process has occurred "from dispersion to clustering". Finally, it is multiplied by the average physical connection weight of all significant causal events within the window. This makes the evolution of aggregation degree also related to the tightness of physical connections.

[0110]

[0111] in, This is a diffusion characteristic (from local to global), where m is the number of backtracking windows (calibrated using historical data, defaulting to m=3). This feature calculates the average change in intensity characteristics between the current window and the previous m windows. This multi-window cumulative difference design smooths out random fluctuations and more robustly reflects trends: when this value is consistently positive and gradually increases, it indicates that not only is the fault intensity increasing, but the rate of increase is also accelerating, clearly demonstrating the dynamic process of anomalies spreading from a local time window to a longer time range (i.e., the "global"). This feature calculation is based on... The information already implies physical weights.

[0112] In summary, the optimized feature function, by introducing the trend progression ratio, the clustering degree change ratio, and the multi-window cumulative difference, allows the formula to directly represent the evolutionary laws of "from weak to strong," "from scattered to clustered," and "from local to global." Simultaneously, it embeds the maximum or average physical connection weights. This ensures that the extracted evolution patterns are consistent with the actual physical correlation strength between devices, enhancing the physical interpretability and engineering consistency of the features.

[0113] Subsequently, the fault evolution trend value is calculated. This fault evolution trend value is smoothed using a weighted moving average method to give higher weight to recent feature windows, thus making them more sensitive to the latest changes in the trend. Its mathematical expression is as follows:

[0114] in, The fault evolution trend value at time point t quantifies the local weighted average change trend of fault risk; t is the time series index; j is the summation index, representing the backtracking distance relative to the current time t; h is the trend smoothing window size, controlling the smoothing degree of the feature sequence, and is determined to be h=0.5w according to empirical rules (w is the sliding window size of the local feature analysis in step S3), so as to achieve a balance between smoothing noise and preserving trend sensitivity; It is an exponentially decaying weight. τ is the attenuation coefficient, which, based on engineering experience, is usually taken between 3 and 5. This weighting makes the values ​​closer to the current time... The greater the contribution to the fault evolution trend value; It is the weighted Euclidean norm; It is an eigenvector.

[0115] Then use the weighted Euclidean norm. This is to amplify the contribution of key features to the overall situation. Weighted Euclidean Norm Defined as:

[0116] in, , and These represent the intensity, aggregation, and diffusion features extracted at time point t, respectively; α, β, and γ are the corresponding feature weight coefficients (α+β+γ=1), whose values ​​are determined through feature importance analysis of historical data (such as feature importance ranking based on random forest or gradient boosting tree models) to ensure that features that are more indicative of fault evolution occupy a higher proportion in norm calculation.

[0117] Finally, based on the timestamps corresponding to each node in the spatiotemporal behavior association chain... and its implied multidimensional feature vectors And its corresponding fault evolution trend value By performing correlation and integration, a structured fault evolution sequence E is generated. To explicitly achieve dual localization of "time + device", a set of key contributing devices is synchronously associated with each node in the generated sequence. This tightly links the time dynamics of failure risk to the specific physical equipment source.

[0118] This key contributing device is for the collection The screening follows a dual criterion to ensure the reliability of its associations: First, based on the current evolutionary indicators for each causal event... Contribution value First, sort the devices in descending order. Second, set a threshold constraint (empirical value > 0.3) to select only device pairs that contribute significantly and have strong physical correlations. Therefore, the fault evolution sequence... The mathematical representation of is optimized as follows:

[0119] in, It's a timestamp; It is the fault evolution trend value; It is the multidimensional feature vector at that moment; It is a set of selected core equipment pairs, in the form of It includes specific device pairs, physical connection weights, and normalized causal strength.

[0120] This sequence fully records the continuous evolution of the equipment's state from early slight deviations, mid-term performance degradation, to late-stage severe anomalies, and establishes the correlation between each evolution step and the core driving equipment. For example, before the oil pump enters the impeller fatigue stage, the periodic slight asymmetry in its vibration mode will be reflected in... The characteristic shows a periodic, slight increase, while simultaneously, the weak causal events leading to this change, from the bearing to the impeller (such as the event with the highest contribution value), will be recorded at the corresponding time. Similarly, for long-distance pipelines, the diffusion characteristics of flow disturbances caused by early corrosion correspond to... The values ​​remained consistently positive, while the core pipeline segment driving this diffusion and its upstream equipment were precisely recorded in the sequence. inside.

[0121] This enhanced structure enables precise localization in subsequent analysis: key inflection points are identified in step S4. At the same time, its associated data can be extracted simultaneously. This directly identifies the core equipment that triggered the sudden change in situation; furthermore, it can be combined with the spatiotemporal behavioral correlation chain of step S2, along... The cause-and-effect relationship is traced backward to locate the source device of the fault. In the prediction of step S5, These can be used as key features input to the LSTM model, enabling the model to not only predict "when the failure will occur" but also simultaneously output "which devices are at high risk." For example, if the failure evolution trend value at t=100s... A sudden increase, corresponding to If the equation is {(A→B,0.8,0.9),(C→B,0.6,0.7)}, then equipment B can be preliminarily identified as the core affected equipment. By tracing back, equipment A can be identified as the potential source of the fault, thus achieving coordinated location and prediction in terms of time and equipment dimensions.

[0122] Through continuous modeling in the above process, this step transforms the equipment's faults from indescribable to a structured sequence that can be identified segment by segment through the evolution path. This sequence will be used in subsequent step S4 to extract key nodes and implement predictive triggers, providing an important basis for fault prediction. The entire method not only considers the temporal characteristics of fault development but also incorporates a feature fusion strategy targeting evolutionary patterns, making the fault evolution sequence both highly interpretable and highly accurate in prediction.

[0123] Step S4: Use the Bayesian change point detection algorithm to identify abrupt change points in the fault evolution sequence to obtain the critical turning points of the oil equipment failure.

[0124] Based on the fault evolution sequence generated in step S3 This step aims to accurately identify key turning points that characterize drastic changes in failure risk. This is to overcome the limitations of processing only one-dimensional trend sequences. To overcome the limitations of traditional methods and fully capture the multi-stage, nonlinear fault evolution characteristics, this method deeply integrates multi-dimensional feature vectors. With core equipment collection This paper constructs a system flow of "multidimensional noise filtering - Bayesian multivariable point detection". The core improvements are: extending the local outlier factor algorithm to multidimensional detection to utilize feature collaborative information and clarify its parameters; upgrading the Bayesian variable point detection model to a framework that supports multivariable points and piecewise multinomial regression; and incorporating the physical connection weights of devices as prior information into the detection process, ultimately outputting the key turning point that is associated with the core driving device set.

[0125] First, multidimensional sequence construction and noise filtering are performed. To achieve collaborative analysis of multidimensional features, fault evolution trend values ​​are... With feature vectors splicing together to form a four-dimensional observation sequence },in Subsequently, an improved multidimensional local outlier factor algorithm was used to analyze the four-dimensional observation sequence { Anomaly detection is performed. The rule for determining the number of neighboring points m, a key parameter, is as follows: , where N is the four-dimensional observation sequence { The total length of the sequence ensures the objectivity of parameter selection and adaptability to sequences of different lengths. The algorithm calculates the multidimensional Euclidean distance between point p and point o:

[0126] in, This is the d-th dimension value of the p-th data point. Based on this distance, the local outlier factor value of each data point is calculated according to the standard procedure. Set the threshold to 1.5, and... Points with a value >1.5 are considered outliers and further categorized: if a point exhibits anomalies in more than two dimensions, it is classified as a "multi-dimensional collaborative anomaly" and directly removed; if it is anomaly only in one dimension, it is labeled as a "single-dimensional anomaly" and handled prudently in conjunction with engineering logic (e.g., only...). abnormal (Under normal circumstances, outliers can be retained but marked as low-risk). After outlier removal, the resulting sequence missing points are filled: if the number of consecutive missing points does not exceed 3, linear interpolation is used; if it exceeds 3, the historical similar sequence library constructed in step S3 is called, and the "interpolation method based on historical similar sequences" is used for filling. After filling, it is necessary to verify that the relative error between the filled data and the adjacent valid data does not exceed 5%; otherwise, the outlier determination logic needs to be re-checked.

[0127] Subsequently, based on the purified four-dimensional observation sequence { Bayesian multivariable point detection is performed. To adapt to the multi-stage and nonlinear evolution characteristics of faults, a piecewise multinomial regression model is adopted. It is assumed that the sequence has K variable points. < <⋯< The sequence is divided into K+1 segments, and the data generation model for each segment m is as follows:

[0128] in, They are polynomial basis vectors. ; It is a 4×3 coefficient matrix; It is observation noise. Furthermore, it allows different segments to have different covariance matrices. For each candidate segment, the linear models (constrained) are compared using the F-test. The last row is zero. The goodness of fit between the model and the quadratic polynomial model is calculated. If the p-value is less than 0.05, the quadratic polynomial model is used to better capture the nonlinear acceleration trend; otherwise, a linear model is used. The optimal number of change points K is determined by minimizing the Bayesian information criterion.

[0129] Where L(K) is the model likelihood; p(K) is the total number of model parameters. The posterior probability distribution at the change point is calculated... At the same time, deeply integrate device association information: Physical connection weights in As a priori constraint, Zhongqiang contribution equipment to ( For the time period corresponding to ≥0.5, a sensitivity enhancement factor η>1 is introduced, which is equivalent to assigning a higher prior probability of the change point occurring during that time period. The posterior probability of the multiple change point locations is calculated using an efficient dynamic programming algorithm.

[0130] Ultimately, each identified critical inflection point of failure The output is enhanced into a structured triple:

[0131] in, It is the location (timestamp) of the critical turning point of the failure. =π( The posterior probability of the critical inflection point of the fault is used to characterize the confidence level. The critical turning point of the failure The set of key contributing devices associated with each other, from This output enables precise location of the sudden change in fault condition in the time dimension and correlation with the source in the device dimension, providing directly usable input for fault location and accurate prediction in step S5.

[0132] Step S5: Construct a fault prediction model based on a long short-term memory network, input the critical fault inflection point and fault evolution sequence into the fault prediction model, and output the fault prediction results across oil equipment.

[0133] The core of this step lies in constructing a time series prediction model based on a Long Short-Term Memory (LSTM) network. To ensure the model input is explicit and computable, the fusion method, specific dimensions, and sequence length of the input features must first be clearly defined. The model input uses the fault evolution sequence E generated in step S3 as the basic time series framework, and incorporates the key inflection points output in step S4. The feature vectors are inserted at the corresponding timestamp positions, thus forming a model input that fuses "continuous evolutionary time-series features" and "discrete key node features." Specifically, at each time step t, the model's input feature vectors... It is a 9-dimensional vector, defined as:

[0134] in, This represents the fault evolution trend value. , , The corresponding three weighted features are intensity, aggregation, and diffusion. ; Key contributing devices to the collection at the current moment Physical connection weights The mean; For the corresponding normalized transfer entropy ratio The mean; Confidence level for inflection points (this value is set to 0 at non-critical inflection point moments); and Encoded values ​​representing the identity of core contributing devices (e.g., the first two bits of one-hot encoding, used to identify the most important source and target devices in the causal relationship).

[0135] The length L of the input sequence to the model is determined by the historical time window required for prediction, and the calculation formula is as follows: .in, The input time window length is typically set to 5 to 10 minutes based on engineering experience, to ensure that the model can capture the complete evolutionary segment from the emergence of the anomaly to the current state. The sampling frequency of the data (in Hz) is used. Therefore, the final input to the model is a three-dimensional tensor of shape (L,9).

[0136] The preprocessed sequence data mentioned above { , ,..., The input is fed into an LSTM network for training. LSTM, through its input gate, forget gate, and output gate mechanisms, can learn and remember long-term dependencies in the fault evolution process. The model ultimately outputs the results for a future period (e.g., future...). The probability distribution of failure of each device node (minutes) P( , ,...∣ , This enables cross-device, forward-looking fault risk prediction.

[0137] Long Short-Term Memory (LSTM) networks are a special type of recurrent neural network capable of effectively capturing long-term dependencies in time series data. Their core lies in selectively memorizing and forgetting information through gating mechanisms. The computation of an LSTM unit at each time step t includes the following gating mechanisms and state updates: The input gate controls the degree to which the current input information is retained, and its calculation formula is as follows:

[0138] in, It is the output value of the input gate at time step t, and its value range is [0,1]. It is the sigmoid activation function, which maps the input to the interval [0,1]. It is the weight matrix of the input gate, obtained through model training; It is the hidden state of the previous time step, carrying historical information; It is the input feature of the current time step, that is, the feature vector extracted from the critical turning point of the fault; This is the bias term of the input gate, learned through training. The input gate uses a sigmoid function to determine which new information needs to be added to the cell state.

[0139] The forget gate controls the degree to which historical information is discarded, and its calculation is expressed as follows:

[0140] in, It is the output of the forget gate, which controls the proportion of cell state retained from the previous time step; It is the weight matrix of the forget gate, which is optimized during training; These are the bias parameters for the forget gate. The forget gate mechanism allows LSTM to selectively retain or discard long-term historical information, avoiding the vanishing or exploding gradient problem.

[0141] After the gating calculation, the cell state needs to be updated. First, the candidate cell state is calculated:

[0142] in, It represents the state of candidate cells, including new information at the current time step; It is the hyperbolic tangent activation function, which constrains the output to the range [-1, 1]. It is the weight matrix of cell states; This is the corresponding bias term. The candidate cell state contains new information that the current input may add to the cell state.

[0143] Update the current cell state based on the outputs of the input gate and the forget gate:

[0144] in, It is the cell state after the current time step update, serving as a carrier of long-term memory; It represents the cell state at the previous time step; dot product represents element-wise multiplication. The updating of the cell state combines selective forgetting of historical information with limited addition of new information, forming the memory mechanism of LSTM.

[0145] The output gate controls the output content in the currently hidden state:

[0146] The final hidden state (output) is calculated as follows:

[0147] in, It is the control signal for the output gate; The hidden state is the current time step and is used as the output of the L unit. The output gate uses a sigmoid function to determine which cell state information should be output to the hidden state.

[0148] During the training phase, historical fault data is used to supervise the learning of the LSTM network. To improve the model's adaptability to different prediction tasks and address the common class imbalance problem in the data, this method abandons the single mean squared error loss and instead dynamically selects a weighted loss function based on the prediction target. Its core logic and specific formulas are as follows: When the prediction objective is to determine "whether the equipment is faulty," and multiple equipment failures are allowed simultaneously, a weighted binary cross-entropy loss is used. This loss function provides an independent binary classification objective for the failure probability of each equipment output by the model, and balances the difference in the number of normal samples and faulty samples in the data through the weight ω, preventing the model from being biased towards the majority class due to sample imbalance. Its formula is defined as:

[0149] Where N is the total number of training samples; D is the total number of devices; ∈{0,1} represents the true state of device d in the i-th sample (0 represents normal, 1 represents fault). The probability of device d failing as predicted by the model; It's weight. , and These represent the total number of normal states and the total number of fault states for all devices in the training set across all samples. When the prediction objective is to accurately determine the specific fault type of the equipment (such as impeller fatigue, bearing wear, pipeline corrosion, etc.), a weighted cross-entropy loss function is used. This loss function assigns a probability distribution across multiple fault types to each piece of equipment and sets different weights for each fault category c. This is used to address the model bias caused by the varying frequencies of different fault types. The formula is defined as:

[0150] Where C represents the total number of fault types; ∈{0,1} indicates whether the true fault type of the i-th sample is c (1 if yes, 0 otherwise); The probability that the model predicts the fault type as c; It is the weight of the c-th type of fault. , To obtain the number of samples of this type of fault in the training set, This represents the total number of samples.

[0151] In summary, if only early warning of fault occurrence is required during model training, then use... Perform optimization; if further identification of the specific fault type is required, then use... Optimization is performed. This dynamic selection mechanism ensures a high degree of match between the loss function and the prediction task, thereby effectively improving the final prediction accuracy of the model.

[0152] During the prediction phase, the features of the critical inflection points and failure evolution sequences monitored in real time are input into the trained LSTM model, and the network outputs the failure probability distribution of each device over a future period:

[0153] Wherein, P( () represents the failure probability distribution over the next T time steps; the softmax function transforms the network output into a probability distribution; and These are the weights and bias parameters of the output layer; It is the hidden state of the current time step.

[0154] Finally, cross-oil equipment failure prediction results are generated.

[0155] The cross-equipment fault prediction result is a structured and quantifiable decision support output. To ensure that the cross-equipment fault prediction result strictly conforms to the physical laws of fault propagation between equipment, this method introduces a cross-equipment association regularization term to constrain the learning process during the model training phase. Specifically, the total loss function of model training is optimized as follows:

[0156] in, The core classification loss is λ = 0.01, the regularization coefficient is C, and C is the set of significant causal relationships identified in step S2. It is the ratio of transfer entropy; It is the physical connection weight; and It predicts the probability of failure. This regularization penalizes those with strong causality. High, physical connection weight Large, but its predicted failure probability and For equipment pairs with excessively large differences, the model is forced to learn the logic that "when the failure probability of the cause equipment increases, the failure probability of the result equipment should also increase synchronously in accordance with physical laws".

[0157] Building upon this, the raw single-device failure probability output by the model will be further processed to generate actionable information encompassing the following dimensions: device location, explicitly indicating the ID of the source device with the highest failure risk and its predicted specific failure type; time window, the predicted time interval for the occurrence or duration of this high-risk state; and probability intensity, the risk level quantified by probability values. For the related impact dimension, refined calculations and classifications will be performed, with the specific rules as follows: First, based on the source device failure probability output by the model. Combined with causal strength Physical connection weights Calculate the probability of a fault being transmitted to related devices. Subsequently, based on The value is used to classify the associated devices into three levels of risk: when When the value is ≥0.3, the equipment is identified as high-risk and must be immediately added to the inspection or maintenance list; when the value is ≤0.1, the equipment is considered high-risk and must be immediately added to the inspection or maintenance list. When the value is less than 0.3, the device is classified as a medium-risk associated device, and the monitoring frequency needs to be increased; when... When the value is less than 0.1, it is considered a low-risk associated device and is only recorded.

[0158] Ultimately, all the information is integrated into a structured instruction that can be directly used by the project, with the following format example: Prediction result = {Fault source device ID:} Fault type: bearing wear; Fault probability: 0.92; Prediction time window: [2023-10-27 14:00, 2023-10-27 16:00]; High-risk associated equipment: [ (Transmission probability = 0.45) (Transmission probability = 0.38)], Medium-risk associated equipment: [ (Transmission probability = 0.15)], Risk transmission path: -> -> Maintenance recommendation: Isolate and inspect immediately. Synchronous inspection and The coordinated execution mechanism constitutes a precise and operable core for preventative maintenance decision-making, capable of automatically triggering tiered early warnings and generating clear inspection work orders, achieving a leap from isolated "point" early warnings to systemic "area" early warnings.

[0159] The results of this cross-oil equipment failure prediction (equipment, time, probability, and related effects) constitute the core of accurate and actionable preventive maintenance decision support, which can automatically trigger early warnings, generate inspection work orders, or provide operational suggestions to operators.

[0160] The innovation of this step lies in combining the time series modeling capabilities of LSTM with the spatiotemporal characteristics of critical inflection points in faults to achieve cross-device correlation prediction. For example, when a pipeline section experiences a propagation trend inflection point, the model not only provides a high-risk warning for crack propagation in that section, but also, based on the spatiotemporal correlation structure of S2, indicates the potential impact of fluctuations on valve chambers and bypass pump groups connected to that section. This provides directly actionable guidance for on-site engineers, significantly enhancing the practical value of fault prediction.

[0161] This technical solution achieves accurate cross-equipment fault prediction in the oilfield industry through a five-step approach. First, raw spatiotemporal data of the equipment is collected and preprocessed with denoising, normalization, and spatiotemporal alignment to lay a reliable data foundation. Next, a variable time-delay transfer entropy algorithm is used to mine causal relationships among the equipment, generating spatiotemporal behavioral correlation chains. A sliding window algorithm is used to extract multi-dimensional features, obtaining a fault evolution sequence. Combining the LOF algorithm with Bayesian change point detection, key fault inflection points are accurately identified. Finally, a model is built based on an LSTM network; inputting the key inflection points and the fault evolution sequence, the model outputs the fault location, type, and impact range. This solution balances data quality and causal correlation, improving the accuracy and foresight of fault prediction and providing strong support for oilfield production safety.

[0162] The present invention also protects a big data-based petroleum equipment fault prediction device, comprising a memory, a processor, and a computer program stored in the memory, wherein the processor executes the computer program to implement the steps of the above method.

[0163] It should be noted that, in this document, relational terms such as "first" and "second" are used only to distinguish one entity or operation from another, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Furthermore, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or apparatus. Without further limitations, the phrase "comprising an element defined as..." does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes said element.

[0164] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their likenesses.

Claims

1. A method for predicting oilfield equipment failures based on big data, characterized in that, Includes the following steps: Step S1: Collect raw spatiotemporal data of each oil equipment and preprocess it to obtain preprocessed spatiotemporal data of the oil equipment; Step S2: Use the variable time-delay transfer entropy algorithm to perform causal analysis on the preprocessed spatiotemporal data of the oil equipment, obtain the causal relationship of the oil equipment, and generate the spatiotemporal behavior association chain of the oil equipment; Step S3: Use the sliding window algorithm to perform local feature analysis on the spatiotemporal behavior association chain to obtain the failure evolution trend value of the oil equipment. Integrate the failure evolution trend value with the timestamps in the spatiotemporal behavior association chain to generate the failure evolution sequence of the oil equipment. Step S4: Use the Bayesian change point detection algorithm to identify abrupt change points in the fault evolution sequence to obtain the critical turning points of the oil equipment failure; Step S5: Construct a fault prediction model based on a long short-term memory network, input the critical fault inflection point and fault evolution sequence into the fault prediction model, and output the fault prediction results across oil equipment.

2. The method for predicting oil equipment failures based on big data according to claim 1, characterized in that, The collection of raw spatiotemporal data from each oil equipment includes: By deploying a sensor network at key nodes in oil production, raw spatiotemporal data reflecting the operating status of equipment are collected. This data can be formally represented as a raw dataset. : ； in, It is the original dataset; It is an identifier that represents the spatial location and device type of the data source; It is the precise timestamp of data collection; It is a vector of multidimensional physical parameters that is strongly correlated with the device type and is measured at that time and location.

3. The method for predicting oil equipment failures based on big data according to claim 2, characterized in that, The preprocessing process, which yields preprocessed spatiotemporal data of the oil equipment, includes: Preprocessing is performed on the raw data in the original dataset to improve its quality and consistency. The preprocessing operations mainly include denoising, normalization and missing value imputation. In the denoising stage, the Kalman filter algorithm is used to suppress measurement noise and extract the true signal trend; In the normalization process, in order to eliminate the differences between different physical dimensions and numerical ranges, the minimum-maximum normalization method is used to scale the data to a uniform range. For missing value imputation, linear interpolation is used to estimate the value, taking into account the characteristics of time series data. Finally, spatial alignment is performed. The primary task of spatial alignment is to construct a spatial topology model of the entire oil equipment, G=(V,E), where V is the vertex set, representing each equipment node, and E is the edge set, representing the physical connection relationship between equipment and the direction of medium flow. Based on spatial topology, time alignment requires precise compensation for the time delay of signal propagation between devices. To improve the accuracy and efficiency of time delay estimation, a "data-driven approach guided by a physical model" is adopted: Based on the different physical nature of signals, the propagation mechanism and velocity value rules are clearly distinguished, thereby providing accurate prior constraints for data-driven computation and estimating the theoretical time delay range according to rules. , ],in and These are the minimum and maximum time delay estimates, respectively; within the theoretical time delay range [ , Under the physical constraints of [ ]; by finding a cross-correlation function Time delay estimation to achieve the maximum value The value can be used to obtain the optimal time delay estimate from device i to device j. Ultimately, the optimal time delay estimate between the two devices is extended to global time synchronization of the entire device cluster; After the above processing, preprocessed spatiotemporal data of petroleum equipment are obtained. .

4. The method for predicting oil equipment failures based on big data according to claim 3, characterized in that, The step of using the variable time-delay transfer entropy algorithm to perform causal analysis on the preprocessed spatiotemporal data of oil equipment to obtain the causal relationships of oil equipment includes: Before calculating the variable time-delay transfer entropy, a physical connection strength threshold λ is set; only when the connection weight... Only when the value is greater than λ will the device pair be retained and subsequent causal strength quantification and direction determination be performed; Subsequently, the variable time-delay transfer entropy is calculated to quantify the causal strength between devices; the basic formula for variable time-delay transfer entropy is defined as: ； in, This represents the variable time-delay transfer entropy value from device X to device Y; It is conditional entropy; It is the state value of device Y at time t; This represents the set of historical states of device Y over the past k time steps; It is the set of past states of device X after time delay adjustment; After calculating the variable time-delay transfer entropy, a normalization process is performed, and the formula is as follows: ； in, It is the ratio of transfer entropy; It is the variable time-delay transfer entropy from X to Y; It is the variable time-delay transfer entropy from Y to X.

5. The method for predicting oil equipment failures based on big data according to claim 4, characterized in that, The spatiotemporal behavioral correlation chain of the generated oil equipment includes: Based on the calculated causal relationships, the next step is to generate spatiotemporal behavioral correlation chains, and the specific construction steps are as follows: First, from all the calculated causal relationships, identify device pairs with a significant causal direction, i.e., those that satisfy... The "cause-effect" combination is used as a candidate link; Subsequently, to address the issue of unclear priority in association chain construction when multiple causal paths may point to the same device, a two-level priority rule was developed to filter and sort causal links. After determining the main link, based on the spatial topology model G=(V,E) and the variable time delay parameters... The selected events are time-aligned; Finally, the cause-effect events that have completed time alignment are linked sequentially to form a chain sequence; The resulting spatiotemporal behavioral correlation chain serves as a structured framework, clarifying the causal order and propagation path of events within the device group.

6. The method for predicting oil equipment failures based on big data according to claim 5, characterized in that, The step of using the sliding window algorithm to perform local feature analysis on the spatiotemporal behavioral correlation chain to obtain the failure evolution trend value of the oil equipment includes: Transforming spatiotemporal behavioral correlation chains into a comprehensive evolutionary index sequence representing health status. Incorporating physical connection weights The optimized evolution index value is obtained by weighted summation of all significant causal events at time t. : ； in, C is the evolution index value at time t; C represents the value that satisfies the following at time t. The set of significant causal relationships (X→Y)>1; It is the normalized ratio of transfer entropy; Next, the evolution index sequence { The sliding window is used for partitioning, and the window size is selected based on the time-scale characteristics of equipment failure evolution. Its mathematical expression is: ； in, It is the k-th sliding window; k is the window index; It is the evolutionary index sequence at time 10:

00. Data points; It is the time index; w is the window size; In each sliding window Within this process, multi-dimensional feature extraction is performed. The feature extraction process is represented as follows: ； in, It is the feature vector of the k-th window; It is an intensity characteristic; It is a clustering characteristic; It is a diffusion characteristic; Subsequently, the fault evolution trend value is calculated. This fault evolution trend value is used to smooth the feature sequence using a weighted moving average method, and its mathematical expression is as follows: ； in, is the fault evolution trend value at time point t; t is the time series index; j is the summation index; h is the window size for trend smoothing; It is an exponentially decaying weight; It is the weighted Euclidean norm; It is an eigenvector.

7. The method for predicting oil equipment failures based on big data according to claim 6, characterized in that, The process of integrating the fault evolution trend value with the timestamps in the spatiotemporal behavior association chain to generate a fault evolution sequence for oil equipment includes: Based on the timestamps corresponding to each node in the spatiotemporal behavior association chain and its implied multidimensional feature vectors And its corresponding fault evolution trend value By performing correlation and integration, a structured fault evolution sequence E is generated, synchronously correlated with a set of key contributing devices. This key contributing device is for the collection The screening follows a dual criterion: only device pairs that make significant contributions and have strong physical correlations are selected. Therefore, the fault evolution sequence sequence The mathematical representation of is optimized as follows: ； in, It's a timestamp; It is the fault evolution trend value; It is the multidimensional feature vector at that moment; It is the set of core equipment pairs selected through screening.

8. The method for predicting oil equipment failures based on big data according to claim 7, characterized in that, The method of using the Bayesian change point detection algorithm to identify abrupt change points in the fault evolution sequence, thereby obtaining the critical inflection points of petroleum equipment failure, includes: To achieve collaborative analysis of multi-dimensional features, fault evolution trend values ​​are... With feature vectors splicing together to form a four-dimensional observation sequence }; Subsequently, an improved multidimensional local outlier factor algorithm was used to analyze the four-dimensional observation sequence { Anomaly detection is performed; subsequently, based on the purified four-dimensional observation sequence... Perform Bayesian multivariable detection, assuming the sequence has K variable points. < <⋯< The sequence is divided into K+1 segments, and the data generation model for each segment m is as follows: ； in, They are polynomial basis vectors; It is a 4×3 coefficient matrix; It is observation noise; the optimal number of change points K is determined by minimizing the Bayesian information criterion; Ultimately, each identified critical inflection point of failure The output is enhanced into a structured triple: ； in, It is the location of the critical turning point of the fault; This represents the posterior probability of the critical inflection point of the fault. It is the critical turning point of the failure. A set of key contributing devices associated with each other.

9. The method for predicting oil equipment failures based on big data according to claim 8, characterized in that, The fault prediction model based on long short-term memory networks inputs the critical inflection points and fault evolution sequences into the model, and outputs fault prediction results across oil equipment, including: Based on the fault evolution sequence E, a time-series framework is established, identifying key inflection points. The feature vectors are inserted at the corresponding timestamp positions to form the model input; Long Short-Term Memory (LSTM) is a special type of recurrent neural network that effectively captures long-term dependencies in time series. Its core lies in selectively memorizing and forgetting information through a gating mechanism. During the training phase, historical fault data is used to supervise the learning of the LSTM network, and the weighted loss function is dynamically selected according to the prediction target. During the prediction phase, the features of the critical inflection points and failure evolution sequences monitored in real time are input into the trained LSTM model, which outputs the failure probability distribution of each device over a future period: ； Wherein, P( () represents the failure probability distribution over the next T time steps; the softmax function transforms the network output into a probability distribution; and These are the weights and bias parameters of the output layer; It is the hidden state at the current time step; Finally, a cross-equipment fault prediction result is generated. This cross-equipment fault prediction result is a structured and quantifiable decision support output. A cross-equipment association regularization term is introduced to constrain the learning process, and the total loss function of model training is optimized as follows: ； in, The core classification loss is λ = 0.01, the regularization coefficient is C, and C is the set of significant causal relationships. It is the ratio of transfer entropy; It is the physical connection weight; and It predicts the probability of failure; Based on this, the raw single-device failure probability output by the model will be further processed to generate actionable information containing the following dimensions: device location, which clearly indicates the ID of the source device with the highest failure risk and its predicted specific failure type; time window, which predicts the time interval during which the high-risk state will occur or last; and probability intensity, which is the risk level quantified by probability value. Ultimately, all the information was integrated into a structured instruction that the project could use directly.

10. A big data-based petroleum equipment fault prediction device, comprising a memory, a processor, and a computer program stored in the memory, characterized in that, The processor executes the computer program to implement the steps of the method according to any one of claims 1-9.

Images