Edge computing-based real-time diagnosis method and system for downhole equipment failure

Abstract

The application provides a kind of based on edge computing's downhole equipment fault real-time diagnosis method and system, it is related to coal mine safety production technical field, including: through distributed sensing network acquisition multimodal data and extract multiscale time series feature, feature is projected to Lie group manifold space and constructs coupling mapping relationship matrix to obtain fusion characteristics, based on topological connection relationship and granger causality coefficient constructs causal directed acyclic graph, executes bayesian probability inference and combines entropy value similarity matching to determine execution strategy, carries out deep time-frequency analysis and causal chain verification.Realize the high-precision real-time diagnosis of equipment fault in the complex environment of downhole, improve the fault early warning accuracy.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of coal mine safety production, and particularly relates to a real-time fault diagnosis method and system for underground equipment based on edge computing. BACKGROUND

[0002] With the deepening of the intelligent construction of coal mines, the safe and stable operation of underground equipment is crucial to the protection of mine production. The underground environment is complex and changeable, and the equipment operation state monitoring and fault diagnosis face many challenges. Traditional underground equipment fault diagnosis mainly relies on artificial experience judgment or single sensor data analysis;

[0003] With the development of Internet of Things and edge computing technology, the fault diagnosis method based on multi-source data has been widely applied. Through the deployment of distributed sensor networks to collect multi-modal data and the use of data analysis technology to realize fault feature extraction and pattern recognition, the introduction of edge computing technology enables data processing to be carried out near the data source, effectively reducing data transmission delay and improving real-time performance;

[0004] However, the existing underground equipment fault diagnosis technology still has the problems of lacking in-depth mining of the internal correlation between different modal data, resulting in limited expression ability of the fused features, failing to fully reflect the complex characteristics of equipment faults, ignoring the physical connection relationship and causal influence between equipment components, leading to lack of explainability of the diagnosis results, difficulty in tracing the fault propagation path, and lack of adaptive diagnosis strategy selection mechanism, which cannot dynamically adjust the analysis depth according to different fault types and the uncertainty of the diagnosis results, and is difficult to realize efficient and accurate real-time diagnosis in the resource-limited edge computing environment. SUMMARY

[0005] The present application provides a real-time fault diagnosis method and system for underground equipment based on edge computing, which can at least solve some of the problems in the prior art.

[0006] In a first aspect, the present application provides a real-time fault diagnosis method for underground equipment based on edge computing, comprising:

[0007] The multi-modal data corresponding to the underground equipment is collected by a distributed sensor network and is subjected to time-frequency domain transformation by an edge computing node. The frequency spectrum features and phase features of each modal signal at different time scales are extracted and organized to obtain multi-scale time sequence features;

[0008] The multi-scale time sequence features are projected into a Lie group manifold space, the geodesic distance and curvature change rate of each modal feature in the manifold space are calculated, and a coupling mapping relationship matrix between the modal features is constructed. The fusion weight of each modal feature is calculated according to the coupling mapping relationship matrix, and the fused features are obtained by weighted summation;

[0009] Obtaining a topological connection relationship of downhole equipment and historical fault sample data, performing time series Granger causality test on the historical fault sample data and calculating Granger causality coefficients between time series data of each component, constructing a causal directed acyclic graph based on the topological connection relationship and the Granger causality coefficients;

[0010] Inputting the fusion features into each node of the causal directed acyclic graph and performing Bayesian probability reasoning along the direction of the directed edges, iteratively updating the fault probability values of each node and solving to obtain a diagnostic reasoning result, calculating the entropy values of the fault probability distribution of each component in the diagnostic reasoning result and performing similarity matching with historical task entropy values corresponding to a preset diagnostic strategy to determine a current execution strategy, and performing deep time-frequency analysis and causal chain verification on the component with the highest fault probability value based on the current execution strategy to obtain an optimized diagnostic result.

[0011] In an optional implementation,

[0012] Collecting multi-modal data corresponding to downhole equipment through a distributed sensing network and performing time-frequency domain transformation through an edge computing node, extracting frequency spectrum features and phase features of each modal signal at different time scales and organizing to obtain a multi-scale time series feature tensor, including:

[0013] Collecting current-voltage signal sequences and temperature-pressure signal sequences corresponding to downhole equipment through a distributed sensing network and performing time series alignment according to timestamps to obtain multi-modal data;

[0014] Sending the multi-modal data to an edge computing node and setting multiple time windows with different lengths at the edge computing node, performing Fourier transform on the current-voltage signal sequences in each time window to obtain frequency domain complex spectra and performing modulus calculation to obtain frequency spectrum energy distributions, and extracting phase difference of harmonic components in the frequency spectrum energy distributions to obtain phase shift features;

[0015] Performing Hilbert-Huang transform on the temperature-pressure signal sequences in each time window, obtaining intrinsic mode functions through empirical mode decomposition and obtaining instantaneous frequency distributions through Hilbert transform, and calculating time derivatives of instantaneous phases based on the instantaneous frequency distributions to obtain phase evolution features;

[0016] Arranging the frequency spectrum energy distributions and the phase shift features according to time window lengths to obtain an electrical modal feature sequence, arranging the instantaneous frequency distributions and the phase evolution features according to time window lengths to obtain a thermodynamic modal feature sequence, and organizing the electrical modal feature sequence and the thermodynamic modal feature sequence according to modal types to obtain a multi-scale time series feature.

[0017] In an optional implementation,

[0018] projecting the multi-scale time sequence features to a Lie group manifold space, calculating geodesic distances and curvature change rates of the modal features in the manifold space, and constructing a coupling mapping relationship matrix between the modal features, calculating fusion weights of the modal features according to the coupling mapping relationship matrix, and performing weighted summation to obtain a fusion feature, including:

[0019] performing nonlinear mapping on the multi-scale time sequence features to obtain manifold points in the Lie group manifold space, constructing a tangent space at the manifold points and calculating base vectors of the tangent space, selecting two modal feature corresponding manifold points as a starting point and an end point, setting a direction vector in the tangent space at the starting point based on the base vectors, performing step movement until the end point is reached, and performing Riemannian metric integration on the step path to obtain a geodesic distance;

[0020] sampling a local neighborhood of the manifold point to obtain neighborhood sampling points and calculating tangent vectors corresponding to the neighborhood sampling points based on the base vectors, transmitting the tangent vectors of the neighborhood sampling points along a closed path and comparing them with an initial tangent vector corresponding to the starting point to obtain a deviation vector, calculating a curvature tensor component based on the deviation vector, and obtaining a scalar curvature through tensor contraction operation, and performing finite difference on the scalar curvature to obtain a curvature change rate;

[0021] mapping the geodesic distance through a preset exponential kernel function to obtain a geometric similarity, calculating a difference between the curvature change rates of different modal features and mapping the difference through an inverse function to obtain an evolution consistency degree, performing element-by-element multiplication operation on the geometric similarity and the evolution consistency degree to obtain a coupling strength, and organizing to obtain a coupling mapping relationship matrix, normalizing row elements in the coupling mapping relationship matrix to obtain fusion weights, and performing weighted fusion on the modal features to obtain a fusion feature.

[0022] In an optional implementation,

[0023] obtaining a topological connection relationship of downhole equipment and historical fault sample data, performing time sequence Granger causality test on the historical fault sample data and calculating Granger causality coefficients between time sequence data of each component, and constructing a causal directed acyclic graph based on the topological connection relationship and the Granger causality coefficients, including:

[0024] obtaining physical connection information between components of downhole equipment and establishing a topological connection relationship according to the connection direction, and obtaining monitoring data of each component before and after a historical fault as historical fault sample data;

[0025] performing stationarity test on time sequence monitoring data of each component in the historical fault sample data, performing difference processing on time sequence monitoring data that does not pass the stationarity test to obtain difference data, and merging the difference data and time sequence monitoring data that passes the stationarity test to obtain stationary time sequence data;

[0026] selecting stationary time series data corresponding to any two components as a target sequence and a reference sequence, shifting the reference sequence by different time lag orders to obtain a lag sequence group, performing multiple linear regression on the target sequence and the lag sequence group to obtain regression coefficients corresponding to each lag order, solving a causal significance statistic based on the regression coefficients and the lag orders, comparing the causal significance statistic with a preset critical value to obtain a Granger causality test result, and calculating a Granger causality coefficient based on regression coefficients that pass the significance test in the Granger causality test result;

[0027] treating each component in the topological connection relationship as a node, establishing a directed edge between components with a Granger causality coefficient greater than a preset threshold value and performing loop detection on the directed edge, and deleting the directed edge forming a loop to obtain a causal directed acyclic graph.

[0028] In an optional implementation,

[0029] inputting the fusion features into each node of the causal directed acyclic graph and performing Bayesian probability reasoning along the direction of the directed edge, iteratively updating the fault probability value of each node, and solving to obtain a diagnostic reasoning result, including:

[0030] allocating the fusion features to the causal directed acyclic graph according to the component type and performing multi-dimensional feature decoupling to obtain a state degradation feature vector and a dynamic evolution feature vector, calculating the posterior distribution parameter of the fault state based on a variational inference algorithm and the state degradation feature vector, and sampling to obtain an initial fault probability value;

[0031] starting from the root node of the causal directed acyclic graph, obtaining the initial fault probability value of the current node and the fault probability value and the dynamic evolution feature vector of the parent node corresponding to the current node, constructing a time-varying conditional probability transition kernel based on the Granger causality coefficient of the directed edge and the dynamic evolution feature vector, and performing tensor convolution operation with the fault probability value to obtain a causal transmission probability distribution, marginalizing the causal transmission probability distribution to obtain an aggregated prior probability, and based on the aggregated prior probability and the initial fault probability value, performing optimization solving to obtain an updated fault probability value and transmitting the updated fault probability value and the dynamic evolution feature vector to the child node, repeating traversal and updating until a single iteration is formed;

[0032] calculating the KL divergence of the fault probability value of each node before and after a single iteration and comparing it with a preset convergence threshold value, if it is less than the convergence threshold value, terminating the iteration, extracting the fault probability value and the dynamic evolution feature vector of each node to calculate a comprehensive fault score, and taking the component corresponding to the node with the maximum comprehensive fault score as the fault source and integrating to obtain a diagnostic reasoning result.

[0033] In an optional implementation,

[0034] The entropy value of each component failure probability distribution in the diagnostic reasoning result is calculated and similarity matching is performed with the historical task entropy value corresponding to the preset diagnostic strategy to determine the current execution strategy, including:

[0035] The failure probability value of each component is extracted from the diagnostic reasoning result and normalized to obtain a failure probability distribution. The information entropy of the failure probability distribution is calculated to obtain a current task entropy value. The causal transmission probability distribution between components in the diagnostic reasoning result is extracted and the conditional entropy is calculated to obtain a causal correlation entropy value. The current task entropy value and the causal correlation entropy value are weighted and summed to obtain a comprehensive entropy feature vector.

[0036] The historical task entropy value and the historical causal correlation entropy value corresponding to each diagnostic strategy in the preset diagnostic strategy library are obtained, and the historical comprehensive entropy feature vector is obtained by weighted summation according to the preset weight.

[0037] The Euclidean distance between the comprehensive entropy feature vector and the historical comprehensive entropy feature vector is calculated and inverse transformed to obtain a distance similarity. The cosine similarity between the current task entropy value and the historical task entropy value is calculated to obtain an entropy value mode similarity. The distance similarity and the entropy value mode similarity are weighted and fused to obtain a comprehensive similarity. The diagnostic strategies are ranked in descending order based on the comprehensive similarity, and the diagnostic strategy ranked first is taken as the current execution strategy.

[0038] In an optional implementation,

[0039] Based on the current execution strategy, the component with the highest failure probability value is subjected to deep time-frequency analysis and causal chain verification to obtain an optimized diagnostic result, including:

[0040] The component with the highest failure probability value is extracted from the diagnostic reasoning result as a target diagnostic component, and the corresponding real-time monitoring time series data is obtained for adaptive wavelet decomposition to obtain frequency band components. The instantaneous frequency and instantaneous amplitude of the frequency band components are calculated to construct a time-frequency joint representation matrix. The gradient change rate corresponding to the time-frequency joint representation matrix is calculated to obtain a time-frequency evolution trajectory, and the amplitude mutation point and the frequency drift point are extracted as abnormal time-frequency feature points. The abnormal time-frequency feature points are subjected to density distribution calculation and peak value detection to obtain an abnormal aggregation time, and the abnormal time period data segment is intercepted based on the abnormal aggregation time.

[0041] extracting parent node components and child node components of the target diagnostic component from the causal directed acyclic graph, constructing a causal verification link and obtaining corresponding historical monitoring time series data, performing cross-correlation analysis on the historical monitoring time series data and the abnormal period data segment to obtain a time lag correlation coefficient, calculating a difference value between the time lag correlation coefficient and a Granger causality coefficient corresponding to the directed edge to obtain a causal bias amount, determining a modified causal verification link according to the causal bias amount and a preset verification threshold;

[0042] calculating a proportion of a remaining component quantity in the modified causal verification link to an initial component quantity to obtain a causal verification confidence, modifying a fault probability value based on the causal verification confidence to obtain a modified fault probability value, and combining the modified fault probability value and the modified causal verification link to obtain an optimized diagnostic result.

[0043] In a second aspect of the embodiment of the present application, a downhole equipment fault real-time diagnosis system based on edge computing is provided, comprising:

[0044] A feature extraction unit is configured to collect multi-modal data corresponding to downhole equipment through a distributed sensing network, perform time-frequency domain transformation through an edge computing node, extract frequency spectrum features and phase features of each modal signal at different time scales, and organize to obtain multi-scale time series features.

[0045] A feature fusion unit is configured to project the multi-scale time series features to a Lie group manifold space, calculate the geodesic distance and curvature change rate of each modal feature in the manifold space, and construct a coupling mapping relationship matrix between the modal features, calculate the fusion weight of each modal feature according to the coupling mapping relationship matrix, and weightedly sum to obtain a fusion feature.

[0046] A graph construction unit is configured to obtain a topological connection relationship of downhole equipment and historical fault sample data, perform time series Granger causality test on the historical fault sample data, calculate Granger causality coefficients between time series data of each component, and construct a causal directed acyclic graph based on the topological connection relationship and the Granger causality coefficients.

[0047] A fault diagnosis unit is configured to input the fusion feature into each node of the causal directed acyclic graph and perform Bayesian probability reasoning along the direction of the directed edge, iteratively update the fault probability value of each node, and solve to obtain a diagnostic reasoning result, calculate the entropy value of the fault probability distribution of each component in the diagnostic reasoning result, and perform similarity matching with a historical task entropy value corresponding to a preset diagnostic strategy to determine a current execution strategy, perform deep time-frequency analysis and causal chain verification on a component with the highest fault probability value based on the current execution strategy to obtain an optimized diagnostic result.

[0048] In a third aspect of the embodiment of the present application, an electronic device is provided, comprising:

[0049] A processor and a memory for storing processor-executable instructions, wherein the processor is configured to invoke the instructions stored by the memory to execute the aforementioned method.

[0050] In a fourth aspect of the embodiments of the present application, a computer readable storage medium is provided, which stores computer program instructions, and the computer program instructions are executed by a processor to implement the aforementioned method.

[0051] In the present application, the multi-scale time sequence features are projected to the Lie group manifold space, the geodesic distance and the curvature change rate are calculated, the coupling mapping relationship matrix is constructed for feature fusion, the nonlinear relationship between different modal signals is effectively captured, the feature expression ability is improved, the Granger causality test and the topological connection relationship are constructed to form a causal directed acyclic graph, the false correlation caused by the traditional correlation analysis is avoided, the explainable causal link is provided for fault diagnosis, the adaptive selection mechanism of the diagnosis strategy based on Bayesian probability reasoning and entropy calculation is realized, the intelligent response to the complex fault situation is realized, the accuracy and efficiency of the diagnosis are significantly improved, the optimization diagnosis process of the deep time-frequency analysis and the causal chain verification is realized, the accuracy of fault positioning is improved, the false positive and false negative rates are reduced, and a reliable basis is provided for the maintenance decision of the downhole equipment. BRIEF DESCRIPTION OF DRAWINGS

[0052] Figure 1 A flowchart of the downhole equipment fault real-time diagnosis method based on edge computing of the embodiments of the present application is shown in the figure.

[0053] Figure 2 A fault source intelligent diagnosis flowchart of the downhole equipment fault real-time diagnosis method based on edge computing of the embodiments of the present application is shown in the figure. DETAILED DESCRIPTION

[0054] To make the objectives, technical solutions and advantages of the embodiments of the present application clearer, the technical solutions in the embodiments of the present application will be described clearly and completely below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, rather than all the embodiments of the present application. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor fall within the scope of protection of the present application.

[0055] The technical solutions of the present application will be described in detail below with specific embodiments. The following specific embodiments can be combined with each other, and the same or similar concepts or processes can not be described in some embodiments.

[0056] Figure 1 A flowchart of the downhole equipment fault real-time diagnosis method based on edge computing of the embodiments of the present application is shown in the figure. Figure 1As shown, the method comprises:

[0057] The multi-modal data corresponding to the downhole equipment is collected through the distributed sensing network and time-frequency domain transformation is performed through the edge computing node, the frequency spectrum features and phase features of each modal signal at different time scales are extracted and multi-scale time sequence features are organized;

[0058] The multi-scale time sequence features are projected into the Lie group manifold space, the geodesic distance and curvature change rate of each modal feature in the manifold space are calculated, and a coupling mapping relationship matrix between the modal features is constructed, the fusion weight of each modal feature is calculated according to the coupling mapping relationship matrix, and the fusion feature is obtained by weighted summation;

[0059] The topological connection relationship and historical fault sample data of the downhole equipment are acquired, the time sequence Granger causality test is performed on the historical fault sample data, and the Granger causality coefficient between the time sequence data of each component is calculated, and a causal directed acyclic graph is constructed based on the topological connection relationship and the Granger causality coefficient;

[0060] The fusion feature is input into each node of the causal directed acyclic graph, and Bayesian probability reasoning is performed along the direction of the directed edge, the fault probability value of each node is iteratively updated, and a diagnostic reasoning result is obtained by solving, the entropy value of the fault probability distribution of each component in the diagnostic reasoning result is calculated, and similarity matching is performed with the historical task entropy value corresponding to the preset diagnostic strategy to determine the current execution strategy, and the component with the highest fault probability value is subjected to deep time-frequency analysis and causal chain verification based on the current execution strategy to obtain an optimized diagnostic result.

[0061] In an alternative embodiment,

[0062] The multi-modal data corresponding to the downhole equipment is collected through the distributed sensing network and time-frequency domain transformation is performed through the edge computing node, the frequency spectrum features and phase features of each modal signal at different time scales are extracted and multi-scale time sequence features are organized, including:

[0063] The current-voltage signal sequence and temperature-pressure signal sequence corresponding to the downhole equipment are collected through the distributed sensing network, and the multi-modal data is obtained by time sequence alignment according to the time stamp;

[0064] The multi-modal data is sent to the edge computing node, a plurality of time windows with different lengths are set in the edge computing node, the frequency domain complex spectrum obtained by performing Fourier transform on the current-voltage signal sequence in each time window is calculated, and the frequency spectrum energy distribution is obtained by modulus calculation, and the phase difference of the harmonic component in the frequency spectrum energy distribution is extracted to obtain the phase shift feature;

[0065] performing Hilbert-Huang transform on the temperature and pressure signal sequence in each time window, obtaining intrinsic mode functions through empirical mode decomposition and obtaining instantaneous frequency distribution through Hilbert transform, and calculating time derivative of instantaneous phase based on the instantaneous frequency distribution to obtain phase evolution characteristics;

[0066] arranging the frequency spectrum energy distribution and the phase shift characteristics according to the time window length to obtain an electrical modal feature sequence, arranging the instantaneous frequency distribution and the phase evolution characteristics according to the time window length to obtain a thermodynamic modal feature sequence, and organizing the electrical modal feature sequence and the thermodynamic modal feature sequence according to modal types to obtain a multi-scale time sequence feature.

[0067] The current and voltage signal sequence and the temperature and pressure signal sequence corresponding to the downhole equipment are collected through a distributed sensing network. The sensing network includes multiple current sensors, voltage sensors, temperature sensors and pressure sensors, which are distributed at key positions of the downhole equipment, such as motor windings, bearings, pump bodies and the like. The sampling frequency of the current and voltage sensor is set to 10 kHz, and the sampling frequency of the temperature and pressure sensor is set to 1 kHz. Each collected data point includes accurate timestamp information for subsequent time sequence alignment processing. For example, for a certain downhole electric pump equipment, the current sensor collects three-phase current values; the voltage sensor collects three-phase voltage values; the temperature sensor collects motor temperature and pump body temperature; and the pressure sensor collects inlet pressure and outlet pressure.

[0068] The data collected by different sensors are aligned to a unified time axis according to the timestamp. During processing, interpolation is performed on data with different sampling frequencies to ensure that all data are accurate and consistent in the time dimension. The temperature and pressure data are improved to the same time precision as the current and voltage data through linear interpolation method to generate an aligned multi-modal data set. The time sequence alignment accuracy is controlled within 0.1 ms to ensure the time correlation of the data.

[0069] The multi-modal data is transmitted to the edge computing node at the wellhead or well site through an industrial Ethernet or wireless communication network. The edge computing node is an industrial-grade device with sufficient computing power, with a processor frequency of not less than 2.5 GHz and a memory capacity of not less than 8 GB, which can meet the real-time computing demand. Multiple time windows of different lengths are set on the edge computing node for data analysis, and the time window lengths are 0.5 seconds, 1 second, 2 seconds, 5 seconds and 10 seconds, covering different time scale features from transient to steady state.

[0070] The Fourier transform is performed on the current and voltage signal sequences in each time window. For each window, the frequency domain complex spectrum of the current or voltage signal is calculated. Taking the first phase current in a 0.5 second window as an example, the time domain signal is converted into a frequency domain representation by the fast Fourier transform algorithm to obtain a complex spectrum containing a real part and an imaginary part. The modulus of the complex spectrum is calculated to obtain the frequency spectrum energy distribution by calculating the square root of the sum of the squares of the real part and the imaginary part. For electrical equipment with a power frequency of 50 Hz, the energy values of the fundamental wave 50 Hz, the second harmonic 100 Hz, the third harmonic 150 Hz, the fifth harmonic 250 Hz, and the seventh harmonic 350 Hz are extracted, and are recorded as the fundamental wave energy, the second harmonic energy, the third harmonic energy, and the like.

[0071] The phase difference of the harmonic components is extracted from the frequency spectrum energy distribution. The phase value is calculated by the inverse tangent of the ratio of the imaginary part to the real part of the complex spectrum. The phase difference between the second harmonic and the fundamental wave, the phase difference between the third harmonic and the fundamental wave, the phase difference between the fifth harmonic and the fundamental wave, and the like are calculated, and the phase difference characteristics constitute a phase shift characteristic set. The equipment in normal operation shows a stable phase difference value, while the phase difference changes significantly in the fault state. Exemplarily, bearing failure can cause the phase difference between the second harmonic and the fundamental wave to increase from 20 degrees in the normal state to more than 35 degrees.

[0072] The Hilbert-Huang transform is performed on the temperature and pressure signal sequences in each time window, and the temperature or pressure signal is decomposed into a plurality of intrinsic mode functions and a residual term by the empirical mode decomposition algorithm. For the pump body temperature signal and the pressure signal, 4 to 6 intrinsic mode function components can be generally decomposed, and the intrinsic mode function components represent temperature and pressure change characteristics of different frequency scales, respectively. The Hilbert transform is performed on each intrinsic mode function to construct an analytical signal, which includes the original intrinsic mode function and the Hilbert transform result thereof.

[0073] The instantaneous frequency distribution is calculated based on the analytical signal obtained by the Hilbert transform, and the instantaneous frequency is defined as the derivative of the instantaneous phase with respect to time, wherein the instantaneous phase is the inverse tangent value of the ratio of the Hilbert transform result to the original intrinsic mode function. In actual calculation, the phase difference between adjacent time points is divided by the time interval to approximately calculate the derivative. For each intrinsic mode function component, the corresponding instantaneous frequency sequence is calculated, and the average frequency, the frequency standard deviation, the frequency peak value, and the like in the window are counted.

[0074] The time derivative of the instantaneous phase is calculated based on the instantaneous frequency distribution to obtain the phase evolution characteristics, and the first order difference sequence and the second order difference sequence of the instantaneous phase are calculated, and the mean, the variance, the skewness, the kurtosis, and the like of the difference sequence are counted. The phase evolution of the equipment in normal operation generally shows a stable change, and the fault equipment can show a sudden change or a nonlinear change mode.

[0075] The spectral energy distribution and phase shift features are arranged in time window length from short to long to form an electrical modal feature sequence. For example, for the first phase current, the corresponding feature sequence includes spectral energy distribution and phase shift features of 0.5 second window, spectral energy distribution and phase shift features of 1 second window, 2 second window features, and so on until 10 second window features.

[0076] The instantaneous frequency distribution and phase evolution features are arranged in time window length from short to long to form a thermodynamic modal feature sequence. Taking the pump body temperature and pressure signals as an example, the feature sequence includes instantaneous frequency distribution and phase evolution features of 0.5 second window, 1 second window features, 2 second window features, and so on until 10 second window features.

[0077] The electrical modal feature sequence and the thermodynamic modal feature sequence are organized by modal type to obtain multi-scale time sequence features. The specific organization method is to group all current-related features into a current modal feature set, group all voltage-related features into a voltage modal feature set, group all temperature-related features into a temperature modal feature set, and group all pressure-related features into a pressure modal feature set.

[0078] In this embodiment, by uniformly collecting heterogeneous signals such as current, voltage, temperature and pressure at the distributed sensing network level and completing strict time sequence alignment based on timestamps, joint representation of different physical quantities under the same time reference is realized, the accuracy and consistency of multi-source information fusion are improved, by setting time windows of multiple lengths and performing frequency domain analysis on electrical signals, not only the main energy distribution features under steady-state operation can be captured, but also the harmonic changes caused by short-time disturbances and abnormal working conditions can be identified, the sensitivity of early fault identification is improved, by performing Hilbert-Huang transform on temperature and pressure signals, nonlinear and non-stationary signals can be adaptively decomposed to obtain instantaneous frequency and phase evolution information, the dynamic evolution process of the device thermodynamic state is more accurately described, and the representation ability of gradual abnormalities under complex working conditions is enhanced.

[0079] In an alternative embodiment,

[0080] The multi-scale time sequence features are projected into the Lie group manifold space, the geodesic distance and curvature change rate of each modal feature in the manifold space are calculated, and a coupling mapping relationship matrix between the modal features is constructed, the fusion weight of each modal feature is calculated according to the coupling mapping relationship matrix, and the fusion feature is obtained by weighted summation, including:

[0081] The manifold points in the Lie group manifold space are obtained by performing nonlinear mapping on the multi-scale time sequence features, a tangent space is constructed at the manifold points and a basis vector of the tangent space is calculated, the manifold points corresponding to two modal features are selected as a starting point and an end point respectively, a direction vector is set in the tangent space at the starting point based on the basis vector, step moving is performed until the end point is reached, and a geodesic distance is obtained by performing Riemannian metric integration on the step path;

[0082] The local neighborhood of the manifold points is sampled to obtain neighborhood sampling points, and tangent vectors corresponding to the neighborhood sampling points are calculated based on the basis vector. After transmitting the tangent vectors of the neighborhood sampling points along a closed path, the tangent vectors are compared with the initial tangent vector corresponding to the starting point to obtain a deviation vector. The curvature tensor components are calculated based on the deviation vector, and a scalar curvature is obtained through tensor contraction operation. The scalar curvature is subjected to finite difference to obtain a curvature change rate;

[0083] The geodesic distance is mapped by a preset exponential kernel function to obtain a geometric similarity. The difference between the curvature change rates of different modal features is calculated and mapped by an inverse function to obtain an evolution consistency degree. The geometric similarity and the evolution consistency degree are subjected to element-by-element multiplication operation to obtain a coupling strength, and a coupling mapping relationship matrix is organized. The row elements in the coupling mapping relationship matrix are normalized to obtain a fusion weight, and the modal features are weighted and fused to obtain a fused feature.

[0084] The manifold points in the Lie group manifold space are obtained by performing nonlinear mapping on the multi-scale time sequence features. The nonlinear mapping is realized by an exponential function, which maps the multi-scale time sequence features in the original feature space to the manifold space. The multi-scale time sequence features are first normalized to ensure that each feature is within the same order of magnitude. For current modal features, the normalization process is to subtract the mean value from the original feature value and divide it by the standard deviation. For temperature modal features and pressure modal features, the maximum and minimum value normalization method is used, that is, the original value is subtracted from the minimum value and divided by the difference between the maximum value and the minimum value. The normalized features are converted to the Lie group manifold space through an exponential mapping function. The base of the mapping function is set to the natural constant e, and the exponential power is the product of the normalized feature value and the mapping coefficient. For current features, the mapping coefficient is set to 0.8. For voltage features, the mapping coefficient is set to 0.75. For temperature features, the mapping coefficient is set to 0.9. For pressure features, the mapping coefficient is set to 0.85.

[0085] A tangent space is constructed at the manifold point and a basis vector of the tangent space is calculated. In the tangent space construction process, the dimension of the manifold is first determined. In this embodiment, the dimension of the manifold is set to 6. The basis vector in the tangent space is calculated using an orthogonal basis method, which is obtained through a Gram-Schmidt orthogonalization process. In the basis vector calculation process, a random vector group is generated, and orthogonalization is achieved by iterative projection to eliminate components. Finally, the normalized orthogonal basis vector group is obtained by unitizing the orthogonalized vectors. In specific implementation, the random seed is set to 42 to generate the initial random vector, the orthogonalization iteration number is set to 10, and the vector unitization accuracy is controlled within 0.0001.

[0086] Two manifold points corresponding to two modal features are selected as the starting point and the ending point. A direction vector is set in the tangent space of the starting point, and step movement is performed until the ending point is reached. The geodesic distance is obtained by integrating the Riemann metric of the step path. The direction vector is set based on the basis vector of the tangent space of the starting point. The vector from the starting point to the ending point is projected onto each basis vector of the tangent space of the starting point to obtain the projection coefficients of each basis vector. The initial direction vector is obtained by weighting and summing the basis vectors with the projection coefficients as weights. The step movement adopts the fourth-order Runge-Kutta method, and the step size is set to 0.05. The upper limit of the iteration step number is 500. In each iteration, the tangent vector is transmitted along the path through parallel transmission to ensure that the moving direction is always within the tangent space of the manifold. The Riemann metric integration adopts the Simpson integration method, and the integration interval is uniformly divided into 100 subintervals. For example, the geodesic distance between the manifold points corresponding to the current modal feature and the temperature modal feature is 3.28.

[0087] The local neighborhood of the manifold point is sampled to obtain neighborhood sampling points. The sampling method adopts spherical uniform sampling to generate uniformly distributed sampling points on the sphere centered at the manifold point. The sampling radius is set to 0.2, and the number of sampling points is 12. The sampling points are distributed on three orthogonal planes, with 4 points uniformly distributed on each plane. For the manifold point corresponding to the current modal feature, the actual distances of the corresponding neighborhood sampling points to the center point are calculated after the neighborhood sampling points are generated to ensure the sampling quality, with a distance deviation requirement of not more than 0.01.

[0088] The tangent vector corresponding to the neighborhood sampling point is calculated based on the basis vector. The tangent vector calculation method is to project the vector difference between the sampling point and the center point onto the basis vector of the center point tangent space to obtain the projection coefficients of each basis vector. The tangent vector is obtained by weighting and summing the basis vectors with the projection coefficients of each basis vector as weights. Exemplarily, for the manifold point and its neighborhood sampling points corresponding to the current modal feature, the calculated tangent vector length is between 0.18 and 0.22.

[0089] The tangent vector of the neighborhood sampling point is transmitted along the closed path and compared with the initial tangent vector corresponding to the starting point to obtain the deviation vector. The closed path is set to start from the starting point, pass through two adjacent neighborhood sampling points, and return to the starting point to form a triangular closed path. The parallel transmission along the closed path adopts a discrete step method, and the step length is set to 0.01. The coordinate transformation of the local tangent space is performed at each step point. For the manifold point corresponding to the current modal feature, the modulus of the deviation vector of the tangent vector after being transmitted along the closed path for one cycle is about 0.032.

[0090] The curvature tensor component is calculated based on the deviation vector, and the scalar curvature is obtained through the tensor contraction operation. The curvature tensor calculation method is to divide the deviation vector by the area enclosed by the closed path to obtain an approximate value of the curvature tensor. The tensor contraction operation adopts the trace method, that is, the sum of the diagonal elements of the curvature tensor is calculated to obtain the scalar curvature. For the manifold point corresponding to the current modal feature, the calculated scalar curvature is about 0.57.

[0091] The curvature rate of change is obtained by finite difference of the scalar curvature. The finite difference method adopts the central difference format. In the local area of the manifold point, the scalar curvature values of the left and right points are taken along a specific direction, and the difference value is calculated by dividing the distance to obtain the curvature rate of change. The difference step length is set to 0.1, and the difference direction is selected as the direction with the most significant change in the scalar curvature. For the manifold point corresponding to the current modal feature, the calculated curvature rate of change is about 1.25.

[0092] The geometric similarity is obtained by mapping the geodesic distance through a preset exponential kernel function. The exponential kernel function is in the form of the negative power of the base number e, and the power is the ratio of the geodesic distance to the kernel width parameter. The kernel width parameter is set to 2.0. For the current modal feature and the temperature modal feature, the geometric similarity calculated based on the geodesic distance 3.28 is about 0.19.

[0093] The evolution consistency degree is obtained by mapping the difference value of the curvature rates of change between different modal features through an inverse function. The calculation method of the difference value of the curvature rates of change is to take the absolute difference value of the curvature rates of change of the manifold points corresponding to the two modal features. The inverse function mapping adopts the form of inverse plus one to ensure that the mapping result is between 0 and 1. For the current modal feature and the temperature modal feature, the curvature rates of change are 1.25 and 2.15 respectively, the difference value is 0.9, and the mapped evolution consistency degree is about 0.53.

[0094] The coupling strength is obtained by performing element-by-element product operation on the geometric similarity and the evolution consistency degree. The element-by-element product operation is to directly multiply the geometric similarity and the evolution consistency degree to obtain the joint measurement value of the two as the coupling strength. For the current modal feature and the temperature modal feature, the coupling strength calculated based on the geometric similarity 0.19 and the evolution consistency degree 0.53 is about 0.10.

[0095] The coupling strength between all pairs of modal features is organized into a coupling mapping relationship matrix. The matrix row and column identifiers are respectively the individual modal features, and the matrix element values are the coupling strength between the corresponding modal features. The diagonal element of the matrix is set to 1, indicating complete coupling of the modal with itself. In this embodiment, the four modal features (current, voltage, temperature, and pressure) form a 4x4 coupling mapping relationship matrix.

[0096] The row elements in the coupling mapping relationship matrix are normalized to obtain the fusion weight. The normalization method is to divide each row of the matrix by the sum of the elements in the row, ensuring that the sum of the elements in each row is 1. For the row corresponding to the current modal feature, the fusion weights obtained after normalization are 0.45, 0.23, 0.19, and 0.13, respectively, corresponding to the fusion weights of the current itself, voltage, temperature, and pressure modal features.

[0097] The weighted fusion of each modal feature obtains the fused feature. In the weighted fusion process, each modal feature is weighted and summed according to the corresponding fusion weight to obtain the fused feature representation. In this embodiment, the dimension of the fused feature is the same as that of each modal feature, and the element value is the weighted sum of the corresponding elements of each modal feature. For the row corresponding to the current modal feature, the first element value of the fused feature is 0.45x1.2+0.23x0.8+0.19x0.6+0.13x0.5=0.895.

[0098] In this embodiment, by performing nonlinear mapping on multi-scale time series features and describing the geodesic distance between feature points in the manifold space, the difference measurement between modal features is no longer limited by linear assumptions and coordinate selection, improving the accuracy of describing modal similarity under complex operating conditions. By introducing the calculation of curvature and its rate of change in the local neighborhood of the manifold, the structural change trend of different modal features in the time evolution process is described, enhancing the sensitivity to gradual and hidden state changes. By coupling the geometric similarity and evolution consistency, and adaptively generating fusion weights based on coupling strength, dynamic adjustment of the contribution degree of different modalities is realized, improving the discrimination ability and stability of the fused feature.

[0099] In an alternative embodiment,

[0100] The topological connection relationship of the downhole equipment and the historical fault sample data are obtained, the time series Granger causality test is performed on the historical fault sample data, the Granger causality coefficients between the time series data of each component are calculated, and a causal directed acyclic graph is constructed based on the topological connection relationship and the Granger causality coefficients, including:

[0101] Obtain physical connection information between components of downhole equipment and establish a topological connection relationship according to the connection direction. Obtain monitoring data of each component before and after the historical failure as historical failure sample data.

[0102] Perform stationarity test on the time series monitoring data of each component in the historical failure sample data. Perform difference processing on the time series monitoring data that does not pass the stationarity test to obtain difference data. Merge the difference data and the time series monitoring data that passes the stationarity test to obtain stationarized time series data.

[0103] Select the stationarized time series data corresponding to any two components as a target sequence and a reference sequence, respectively. Shift the reference sequence according to different time lag orders to obtain a lag sequence group. Perform multiple linear regression on the target sequence and the lag sequence group to obtain regression coefficients corresponding to each lag order. Solve to obtain a causality significance statistic based on the regression coefficients and the lag orders. Compare the causality significance statistic with a preset critical value to obtain a Granger causality test result. Calculate a Granger causality coefficient based on the regression coefficients that pass the significance test in the Granger causality test result.

[0104] Take each component in the topological connection relationship as a node. Establish a directed edge between the components whose Granger causality coefficients are greater than a preset threshold value and perform loop detection on the directed edge. Delete the directed edge that forms a loop to obtain a causal directed acyclic graph.

[0105] Obtain physical connection information between components of downhole equipment and establish a topological connection relationship. Downhole equipment is usually composed of multiple components such as motors, pump bodies, bearings, and sealing devices. There is a clear physical connection relationship between the components. The physical connection information is obtained through device structure diagrams and field installation records, including the relative positions and connection methods of the components. Determine the connection direction between the components according to the energy flow and signal transmission direction, for example, the motor drives the bearing, so the connection direction points from the motor to the bearing; the bearing supports the pump body, so the connection direction points from the bearing to the pump body. Take a downhole electric pump as an example. Its components include a power control unit, a motor stator, a motor rotor, a main bearing, an auxiliary bearing, a sealing device, a pump body inlet, a pump body impeller, and a pump body outlet. The power control unit is connected to the motor stator, the motor stator drives the motor rotor to rotate, the motor rotor is connected to the sealing device through the main bearing and the auxiliary bearing, the sealing device is connected to the pump body inlet, the pump body inlet is connected to the pump body impeller, and the pump body impeller is connected to the pump body outlet, forming a complete topological connection relationship.

[0106] The monitoring data of each component before and after the historical failure is obtained as historical failure sample data. The monitoring data is collected by various sensors installed on the downhole equipment, including current sensors, voltage sensors, temperature sensors, vibration sensors, etc. For historical failure events, monitoring data from 7 days before failure to 1 day after failure is extracted, with a sampling frequency of 10 sampling points per second for current and voltage, 1 sampling point per minute for temperature, and 50 sampling points per second for vibration. In data preprocessing, obvious outliers are removed, such as values exceeding 3 times the standard deviation of the normal working range, and missing data points are completed by linear interpolation. Taking bearing failure as an example, temperature data and vibration data of the main bearing before and after failure, current data of the motor stator, pressure data of the pump body inlet, etc. are collected. The temperature of the main bearing gradually rises from 65°C to 80°C within 4 days before failure, and the vibration amplitude rises from 0.5 mm / s to 2.8 mm / s, while the motor current fluctuation increases by 20%.

[0107] The time series monitoring data of each component in the historical failure sample data is subjected to stationarity test. The augmented Dickey-Fuller test method is used to test each type of monitoring data sequence of each component. The significance level is set to 0.05 during the test, and the lag order is automatically selected based on the Schwarz information criterion, with a maximum lag order of 12. After stationarity test of each time series data, the probability value is obtained, and the probability value is compared with the significance level. If the probability value is less than 0.05, the stationarity test is passed, otherwise the difference processing is needed. The temperature data of the main bearing is subjected to stationarity test, and the probability value is 0.32, which is greater than the significance level 0.05, so it is judged that the time series is not stationary; the current data of the motor stator is tested to obtain a probability value of 0.02, which is less than the significance level, and it is judged that the time series is stationary.

[0108] The time series monitoring data that does not pass the stationarity test is subjected to difference processing to obtain difference data. Difference processing is achieved by calculating the difference between adjacent time points, i.e. the value of the current time point minus the value of the previous time point. The difference order is determined according to the characteristics of the data, generally first-order difference is performed first, if the first-order difference still does not pass the stationarity test, second-order difference is performed, and at most third-order difference is performed. The temperature data of the main bearing is subjected to first-order difference, and the probability value obtained by stationarity test again is 0.03, which is less than the significance level 0.05, so the first-order difference result is passed. The difference of the main bearing temperature data reflects the temperature change rate, and the data range is from -0.5°C / hour to 1.8°C / hour.

[0109] The differential data is combined with the time series monitoring data that passes the stationarity test to obtain stationary time series data. The time correspondence of each data sequence is maintained during the combination process, and for the differential processed data, the processed sequence is used instead of the original sequence. The combined stationary time series data includes stationary original data and differential data, which together constitute a stationary data set. The data of the main bearing temperature after first-order difference and the originally stationary motor stator current data together constitute the stationary time series data set of the corresponding components.

[0110] The stationary time series data corresponding to any two components are selected as the target sequence and the reference sequence, respectively. Taking the main bearing and the motor stator as an example, the first-order differential data of the main bearing temperature is selected as the target sequence, and the motor stator current data is selected as the reference sequence. The length of the target sequence is the total number of sampling points for 8 days, and the length of the motor stator current data as the reference sequence is the same.

[0111] The reference sequence is shifted according to different time lag orders to obtain a lag sequence group. The selection range of the lag order is 1 to 10, i.e. a sequence lagging 1 time unit, lagging 2 time units, and so on to lagging 10 time units is generated. The motor stator current data is shifted by 1 to 10 time units to obtain 10 lag sequences, which are combined with the original sequence to form a lag sequence group. Laging 1 time unit means shifting the motor stator current data as a whole by one sampling point, i.e. using the value of the previous time point at the current time point.

[0112] Multiple linear regression is performed on the target sequence and the lag sequence group to obtain the regression coefficients corresponding to each lag order. In the multiple linear regression model, the target sequence is the dependent variable, and each sequence in the lag sequence group is the independent variable. The regression calculation uses the least squares method to calculate the regression coefficients and standard errors corresponding to each lag order. In the regression analysis of the main bearing temperature differential data on the motor stator current lag sequence, the regression coefficient of lag 1 is 0.42, the standard error is 0.11; the regression coefficient of lag 2 is 0.35, the standard error is 0.13; the regression coefficient of lag 3 is 0.28, the standard error is 0.12; the regression coefficient of lag 4 is 0.15, the standard error is 0.10; and the regression coefficients of the remaining orders are all less than 0.1.

[0113] The causal significance statistic is solved based on the regression coefficients and the lag orders. The causal significance statistic is calculated by the variance ratio statistic, i.e. the sum of the squares of all lag order regression coefficients divided by the sum of squares of residuals multiplied by the degree of freedom ratio. In the calculation of the variance ratio statistic, the numerator degree of freedom is the lag order, and the denominator degree of freedom is the total number of samples minus the number of regression coefficients minus 1. In the causal test of the motor stator current on the main bearing temperature change, the calculated variance ratio statistic is 8.76.

[0114] The Granger causality test result is obtained by comparing the causality significance statistics with a preset threshold value. The threshold value is determined based on a quantile of a variance ratio distribution, and a significance level is set to 0.01. In this embodiment, the threshold value of the variance ratio distribution is 4.25, and because the calculated variance ratio statistics 8.76 is greater than the threshold value 4.25, it is determined that there is a significant Granger causality relationship between the motor stator current and the main bearing temperature change, that is, the change of the motor stator current can predict the change of the main bearing temperature. Similarly, the Granger causality test results between other component pairs are calculated.

[0115] The Granger causality coefficient is calculated based on the regression coefficients of the Granger causality test results that pass the significance test. The Granger causality coefficient is calculated by a weighted sum of the regression coefficients, and the weight is the reciprocal of each lag order. Only the regression coefficients that pass the significance test are involved in the calculation, and the standard is that the t-statistic obtained by dividing the regression coefficient by its standard error is greater than the t-distribution threshold value corresponding to the degree of freedom. The Granger causality coefficient of the motor stator current to the main bearing temperature is calculated as 0.37. Similarly, the Granger causality coefficients between other component pairs are calculated.

[0116] Each component in the topological connection relationship is taken as a node, and a directed edge is established between components with a Granger causality coefficient greater than a preset threshold value. The preset threshold value is set to 0.2. For component pairs with a Granger causality coefficient greater than 0.2, a directed edge is established from the dependent variable component to the independent variable component, indicating the direction of the causal relationship. According to the calculation result, a directed edge is established between the motor stator and the main bearing with a Granger causality coefficient of 0.37; a directed edge is established between the main bearing and the pump body inlet with a Granger causality coefficient of 0.25; and a directed edge is established between the motor rotor and the auxiliary bearing with a Granger causality coefficient of 0.31.

[0117] Loop detection is performed on the directed edges, and the directed edges forming a loop are deleted to obtain a causal directed acyclic graph. The loop detection uses a depth-first search algorithm. Starting from each node, traversal is performed along the directed edges. If the starting node is returned during traversal, it is judged that a loop exists. After detecting a loop, the directed edge with the smallest Granger causality coefficient in the loop is deleted, and loop detection is repeated until there is no loop. In this embodiment, a loop is detected between the motor rotor and the auxiliary bearing, and the Granger causality coefficient from the auxiliary bearing to the motor rotor is 0.22, which is smaller than 0.31 from the motor rotor to the auxiliary bearing. Therefore, the directed edge from the auxiliary bearing to the motor rotor is deleted, and finally a loop-free causal directed graph is obtained.

[0118] In this embodiment, by introducing the real physical connection relationship between components as a prior constraint, the potential action path is limited in the causal modeling process, effectively avoiding the non-physical association and false causal relationship in the traditional pure data-driven method, improving the targeted description ability of the fault triggering and propagation behavior, through the stationarity test on the component monitoring data and the difference processing on the non-stationary sequence, the risk of false regression and false causal relationship is significantly reduced, the stability and credibility of the Granger causality test result is improved, by introducing the regression analysis of multiple time lag order and selecting the causal relationship based on the significance statistic, the fault influence path and key conduction node can be described more finely, which is helpful to identify the core components in fault propagation.

[0119] In an alternative embodiment,

[0120] Inputting the fusion features into each node of the causal directed acyclic graph and performing Bayesian probability reasoning along the direction of the directed edge, iteratively updating the fault probability value of each node and solving to obtain a diagnostic reasoning result includes:

[0121] According to the component type, the fusion features are allocated to the causal directed acyclic graph and decoupled in multiple dimensions to obtain a state degradation feature vector and a dynamic evolution feature vector, the posterior distribution parameters of the fault state are calculated based on the variational inference algorithm and the state degradation feature vector, and the initial fault probability value is sampled;

[0122] Starting from the root node of the causal directed acyclic graph, the initial fault probability value of the current node and the fault probability value and dynamic evolution feature vector of the parent node corresponding to the current node are obtained, the time-varying conditional probability transition kernel is constructed based on the Granger causal coefficient of the directed edge and the dynamic evolution feature vector, and the tensor convolution operation is performed on the fault probability value to obtain the causal transmission probability distribution, the marginal integral is performed on the causal transmission probability distribution to obtain the aggregated prior probability, and the updated fault probability value is obtained by optimization solving based on the aggregated prior probability and the initial fault probability value, and is transmitted to the child node with the dynamic evolution feature vector. Repeat traversal and update until a single iteration is formed;

[0123] The KL divergence of the fault probability value of each node before and after a single iteration is calculated and compared with a preset convergence threshold, if it is less than the convergence threshold, the iteration is terminated, the comprehensive fault score is calculated based on the fault probability value and the dynamic evolution feature vector of each node, and the component corresponding to the node with the maximum comprehensive fault score is taken as the fault source, and the diagnostic reasoning result is integrated.

[0124] The fusion feature is assigned to the causal directed acyclic graph according to the component type, and the state degradation feature vector and the dynamic evolution feature vector are obtained by multi-dimensional feature decoupling. In the process of assigning the fusion feature, different mapping rules are adopted for different component types to map the fusion feature to the corresponding node. The fusion features of motor components include current, voltage, temperature and other dimensions, all dimensions are kept during mapping; the fusion features of bearing components include temperature, vibration, noise and other dimensions, all dimensions are kept during mapping; the fusion features of pump body components include pressure, flow, vibration and other dimensions, all dimensions are kept during mapping. The feature decoupling adopts independent component analysis method to decompose the fusion feature into independent components. The number of decoupled components is set to 2, the fast independent component analysis algorithm is used for solving, the iteration number is set to 1000, and the convergence threshold is set to 0.0001, to obtain the state degradation feature vector and the dynamic evolution feature vector. For the main bearing component, the fusion feature is a 16-dimensional vector, and after feature decoupling, an 8-dimensional state degradation feature vector and an 8-dimensional dynamic evolution feature vector are obtained. The state degradation feature vector mainly represents the steady-state fault feature of the component, reflecting the inherent performance degradation of the component; the dynamic evolution feature vector represents the dynamic fault feature of the component, reflecting the trend of the running state of the component.

[0125] The posterior distribution parameters of the fault state are calculated based on the variational inference algorithm and the state degradation feature vector, and the initial fault probability value is obtained by sampling. The variational inference algorithm uses the mean field approximation method, takes the state degradation feature vector as the observation data, establishes a Gaussian mixture model as the prior distribution, and sets the number of mixture components to 3, corresponding to the normal state, the slight fault state and the severe fault state. The variational inference algorithm iteratively solves the posterior distribution parameters, and the iteration number is set to 500 and the convergence threshold is set to 0.0001. The posterior distribution parameters include the mean vector and the covariance matrix. For the main bearing component, the mean vector of the normal state is [0.12, 0.08, 0.15, 0.10, 0.07, 0.09, 0.11, 0.13], the mean vector of the slight fault state is [0.45, 0.38, 0.42, 0.40, 0.35, 0.39, 0.43, 0.41], and the mean vector of the severe fault state is [0.82, 0.78, 0.85, 0.80, 0.75, 0.79, 0.83, 0.81]. 1000 samples are sampled from the posterior distribution, the distribution proportion of each state is calculated, and the probability value of the fault state is calculated. For the main bearing component, the sampling result shows that the probability of the normal state is 0.25, the probability of the slight fault state is 0.55, and the probability of the severe fault state is 0.20, and the initial fault probability value is the sum of the probabilities of the slight fault and the severe fault, i.e. 0.75.

[0126] From the root node of the causal directed acyclic graph, the initial failure probability value of the current node and the failure probability value and dynamic evolution feature vector of the parent node corresponding to the current node are obtained. The root node refers to the node with an in-degree of 0, that is, the node without a directed edge pointing to it. In the causal directed acyclic graph of the downhole device, the power control unit is usually the root node. Taking the motor stator node as an example, the corresponding initial failure probability value is 0.35, the failure probability value of the parent node power control unit is 0.28, and the dynamic evolution feature vector transmitted by the power control unit is [0.23, 0.25, 0.20, 0.22, 0.24, 0.21, 0.26, 0.19].

[0127] The time-varying conditional probability transition kernel is constructed based on the Granger causality coefficient of the directed edge and the dynamic evolution feature vector, and the tensor convolution operation is performed with the failure probability value to obtain the causal transmission probability distribution. The time-varying conditional probability transition kernel represents the influence degree of the parent node failure state on the child node failure state, which is constructed based on the Granger causality coefficient and the dynamic evolution feature vector. The transition kernel construction method is to take the Granger causality coefficient as the basic weight, and then adjust it according to the dimension value of the dynamic evolution feature vector. The adjustment formula is the basic weight multiplied by the weighted sum of the dynamic evolution feature, wherein the weighting coefficient is a preset influence factor. For the directed edge from the power control unit to the motor stator, the Granger causality coefficient is 0.42, and the adjustment coefficient calculated based on the dynamic evolution feature vector is 1.25, so the time-varying conditional probability transition kernel is 0.525. The tensor convolution operation is performed between the parent node failure probability value and the time-varying conditional probability transition kernel, that is, the parent node failure probability value is multiplied by the time-varying conditional probability transition kernel, to obtain the causal transmission probability distribution. The causal transmission probability distribution of the power control unit to the motor stator is 0.28 multiplied by 0.525, equal to 0.147.

[0128] The aggregated prior probability is obtained by marginalizing the integral of the causal transmission probability distribution. Marginalizing integration refers to the process of integrating or summing the conditional probability over all possible parent node states. When a node has multiple parent nodes, the influence of all parent nodes needs to be considered comprehensively. The marginalization calculation method is to add the causal transmission probability distributions of all parent nodes, and then divide by the number of parent nodes to obtain the aggregated prior probability. For the motor stator node, it has only one parent node, the power control unit, so the aggregated prior probability is equal to the corresponding causal transmission probability distribution 0.147.

[0129] The updated failure probability value is obtained by optimization solving based on the aggregated prior probability and the initial failure probability value, and is transmitted to the child node together with the dynamic evolution feature vector. The optimization solving adopts the Bayesian updating method, taking the aggregated prior probability as the prior and the initial failure probability value as the likelihood, and calculating the posterior probability as the updated failure probability value. The calculation formula is the prior probability multiplied by the likelihood probability divided by the normalization constant. The normalization constant is ensured to be 1 by integration or summation. For the motor stator node, the aggregated prior probability is 0.147, the initial failure probability value is 0.35, and the updated failure probability value is calculated as 0.41. The updated failure probability value is transmitted to the child node of the motor stator, i.e. the motor rotor, together with the dynamic evolution feature vector [0.32, 0.30, 0.35, 0.33, 0.29, 0.31, 0.34, 0.28] of the motor stator. The above process is repeated to traverse all nodes in the causal directed acyclic graph, and one iteration update is completed.

[0130] The Kullback-Leibler divergence of the failure probability values of each node before and after a single iteration is calculated and compared with the preset convergence threshold. If it is less than the convergence threshold, the iteration is terminated. The Kullback-Leibler divergence is used to measure the difference between two probability distributions. The calculation method is to multiply the logarithm of the iteration probability value by the iteration probability value, and then subtract the logarithm of the iteration probability value multiplied by the iteration probability value. The divergence value is obtained by summing all states. The preset convergence threshold is 0.01. The failure probability value of the motor stator node before the first iteration is 0.35, and after the first iteration, it is 0.41. The calculated Kullback-Leibler divergence is 0.024, which is greater than the convergence threshold 0.01, so the iteration continues. After the second iteration, the failure probability value of the motor stator node is 0.43, and the calculated Kullback-Leibler divergence is 0.008, which is less than the convergence threshold 0.01. At this time, it is judged that the motor stator node converges. When the Kullback-Leibler divergence of all nodes is less than the convergence threshold, the overall iteration process is terminated.

[0131] The comprehensive failure score is calculated by extracting the failure probability value and the dynamic evolution feature vector of each node. The comprehensive failure score calculation method is the weighted sum of the failure probability value and the norm of the dynamic evolution feature vector, and the weight coefficients are 0.7 and 0.3 respectively. The norm calculation method of the dynamic evolution feature vector is the square root of the sum of squares of each dimension. For the main bearing node, the final converged failure probability value is 0.82, the norm of the dynamic evolution feature vector is 0.65, and the calculated comprehensive failure score is 0.77. For the motor stator node, the final converged failure probability value is 0.43, the norm of the dynamic evolution feature vector is 0.42, and the calculated comprehensive failure score is 0.43. After calculating the comprehensive failure score of all nodes, the component corresponding to the node with the maximum comprehensive failure score is determined as the fault source. In this embodiment, the comprehensive failure score of the main bearing node is 0.77, which is the highest, so it is judged that the main bearing is the fault source.

[0132] The fault information of the fault source component and its associated components is integrated to obtain a diagnosis reasoning result, the comprehensive fault score of the fault source component is taken as a main index, and the comprehensive fault score of the child node of the fault source component is taken as an associated fault index to form a complete diagnosis report. The diagnosis report includes the fault source component name, the fault probability, the fault score and the possible fault type. The fault type is determined according to the pattern recognition of the state degradation feature vector, and the fault type with the highest matching degree is found by comparing with the pre-labeled fault pattern library. For the main bearing component, the matching degree of the state degradation feature vector and the bearing inner ring wear fault pattern reaches 92%, and therefore the fault type is judged as inner ring wear. The final diagnosis reasoning result is that the main bearing component has an inner ring wear fault, the fault probability is 0.82, the comprehensive fault score is 0.77, the associated affected components include the sealing device, the fault probability is 0.45, the comprehensive fault score is 0.38, and the sealing may be not tight, and it is suggested to replace the main bearing and check the sealing device.

[0133] In the embodiment, the fusion features are mapped to corresponding nodes in the causal directed acyclic graph according to the component type, the features are decoupled and modeled for state degradation and dynamic evolution, the fault representation is no longer mixed with information of different time scales and physical meanings, the different influences of long-term degradation and short-term disturbance on the fault probability are distinguished, the precision and reliability of the fault state description are improved, the probability propagation is performed along the causal direction in the causal directed acyclic graph, the time-varying conditional probability transfer mechanism determined by the Granger causality coefficient and the dynamic evolution feature is introduced, the transmission process of the fault influence can be adaptively adjusted with the running state, the physical consistency and time sequence rationality of the fault propagation modeling are improved, and the false judgment and probability shock phenomenon are significantly reduced by aggregating the causal transmission probability distribution and jointly optimizing the initial fault probability of the node itself.

[0134] Figure 2 The fault source intelligent diagnosis flowchart of the downhole equipment fault real-time diagnosis method based on edge computing in the embodiment of the application.

[0135] In an alternative embodiment,

[0136] The entropy value of the fault probability distribution of each component in the diagnosis reasoning result is calculated, and similarity matching is performed with the historical task entropy value corresponding to the preset diagnosis strategy to determine the current execution strategy.

[0137] extracting a fault probability value of each component from the diagnostic reasoning result and performing normalization processing to obtain a fault probability distribution, calculating an information entropy of the fault probability distribution to obtain a current task entropy value, extracting a causal transmission probability distribution between each component in the diagnostic reasoning result and calculating a conditional entropy to obtain a causal correlation entropy value, and performing weighted summation on the current task entropy value and the causal correlation entropy value to obtain a comprehensive entropy feature vector;

[0138] obtaining a historical comprehensive entropy feature vector by obtaining a historical task entropy value and a historical causal correlation entropy value corresponding to each diagnostic strategy in a preset diagnostic strategy library and performing weighted summation according to a preset weight;

[0139] calculating an Euclidean distance between the comprehensive entropy feature vector and the historical comprehensive entropy feature vector and performing inverse transformation to obtain a distance similarity, calculating a cosine similarity between the current task entropy value and the historical task entropy value to obtain an entropy value mode similarity, performing weighted fusion on the distance similarity and the entropy value mode similarity to obtain a comprehensive similarity, and arranging the diagnostic strategies in descending order based on the comprehensive similarity and taking a diagnostic strategy ranked first as a current execution strategy.

[0140] A fault probability value of each component is extracted from the diagnostic reasoning result and normalized to obtain a fault probability distribution. The diagnostic reasoning result contains fault probability values of each component. After extracting the fault probability values, normalization processing is needed to ensure data comparability. The normalization processing adopts a summation normalization method, that is, the fault probability values of each component are divided by the sum of all component fault probability values, so that the sum of the normalized fault probability values is 1. Taking a downhole electric pump system as an example, the diagnostic reasoning result contains 9 components: power supply control unit, motor stator, motor rotor, main bearing, auxiliary bearing, sealing device, pump body inlet, pump body impeller and pump body outlet. The original fault probability values are 0.28, 0.43, 0.39, 0.82, 0.35, 0.45, 0.30, 0.25 and 0.20, and the sum of the fault probability values is 3.47. The normalized fault probability distribution is 0.081, 0.124, 0.112, 0.236, 0.101, 0.130, 0.086, 0.072 and 0.058, respectively.

[0141] The information entropy is calculated by the fault probability distribution of the current task. The information entropy is used to measure the degree of uncertainty of the system. The calculation method is to multiply the normalized fault probability value of each component by the negative of its logarithmic value, and then add all the results. The base of the logarithm is set to the natural logarithm base e. For the fault probability distribution of the downhole electric pump system, the calculated information entropy value is 2.14, and the information entropy value reflects the complexity and uncertainty of the current fault diagnosis task. The higher the current task entropy value, the greater the uncertainty of fault diagnosis, and more complex diagnosis strategy is needed; the lower the current task entropy value, the more determined the fault location, and a targeted diagnosis strategy can be used.

[0142] The causal correlation entropy value is obtained by extracting the causal transmission probability distribution between components in the diagnosis reasoning result and calculating the conditional entropy. The diagnosis reasoning result contains the causal transmission probability distribution between components, which reflects the possible path and strength of fault propagation between components. The conditional entropy calculation method is to calculate the uncertainty of the target component fault probability distribution under the condition that the source component fault probability is known for each component pair. The fault probability value of the source component is multiplied by the corresponding conditional probability to obtain the joint probability. The negative of the logarithmic value of the joint probability is calculated, and then multiplied by the joint probability. Sum all possible component pairs. In the downhole electric pump system, the causal transmission probability from the power supply control unit to the motor stator is 0.525, the causal transmission probability from the motor stator to the motor rotor is 0.480, and the causal transmission probability from the motor rotor to the main bearing is 0.420. The complete causal transmission probability matrix is constructed in this way, and the conditional entropy calculated based on these causal transmission probabilities is 1.75, which reflects the complexity and uncertainty of the fault propagation relationship between components.

[0143] The comprehensive entropy feature vector is obtained by weighted sum of the current task entropy value and the causal correlation entropy value. The weight coefficient of weighted sum is set by experience. The weight coefficient of the current task entropy value is set to 0.6, and the weight coefficient of the causal correlation entropy value is set to 0.4. The current task entropy value of the downhole electric pump system is 2.14, and the causal correlation entropy value is 1.75. The weighted sum of the comprehensive entropy feature vector is 2.14x0.6+1.75x0.4=1.984. The comprehensive entropy feature vector considers the uncertainty of fault distribution and the complexity of fault propagation relationship, which provides a basis for subsequent diagnosis strategy selection.

[0144] The historical task entropy value and the historical causal correlation entropy value corresponding to each diagnostic strategy in the preset diagnostic strategy library are obtained, and a historical comprehensive entropy feature vector is obtained by weighting and summing according to a preset weight. The diagnostic strategy library contains various diagnostic strategies for different fault scenarios, and each strategy has a corresponding applicable scene feature. The historical task entropy value and the historical causal correlation entropy value are obtained based on historical diagnostic cases, and reflect the characteristics of the applicable scene of various diagnostic strategies. The diagnostic strategy library of the downhole electric pump system contains five strategies: single-component deep diagnostic strategy, multi-component parallel diagnostic strategy, causal chain tracking diagnostic strategy, full-system scanning diagnostic strategy, and regional step-by-step diagnostic strategy. The historical task entropy value of the single-component deep diagnostic strategy is 1.2, and the historical causal correlation entropy value is 0.8; the historical task entropy value of the multi-component parallel diagnostic strategy is 1.8, and the historical causal correlation entropy value is 1.5; the historical task entropy value of the causal chain tracking diagnostic strategy is 2.2, and the historical causal correlation entropy value is 1.9; the historical task entropy value of the full-system scanning diagnostic strategy is 2.5, and the historical causal correlation entropy value is 2.0; the historical task entropy value of the regional step-by-step diagnostic strategy is 1.9, and the historical causal correlation entropy value is 1.6. The historical entropy values of each strategy are weighted and summed using the same weight coefficient as the current comprehensive entropy feature vector to obtain the historical comprehensive entropy feature vectors of each strategy, which are 1.04, 1.68, 2.08, 2.30 and 1.78, respectively.

[0145] The Euclidean distance between the comprehensive entropy feature vector and the historical comprehensive entropy feature vector is calculated and inversely transformed to obtain the distance similarity. The Euclidean distance calculation method is to subtract the absolute value of each historical comprehensive entropy feature vector from the current comprehensive entropy feature vector. The inverse transformation adopts the method of 1 divided by the Euclidean distance plus 1, which ensures that the value range of the distance similarity is between 0 and 1. The current comprehensive entropy feature vector of the downhole electric pump system is 1.984, and the Euclidean distances of each historical comprehensive entropy feature vector are 0.944, 0.304, 0.096, 0.316 and 0.204, respectively. After inverse transformation, the distance similarities are 0.514, 0.767, 0.912, 0.760 and 0.830, respectively. The higher the distance similarity, the closer the characteristics of the current diagnostic task to the applicable scene of the corresponding diagnostic strategy.

[0146] The cosine similarity between the current task entropy value and the historical task entropy value is calculated to obtain the entropy value mode similarity. The cosine similarity is obtained by calculating the product of the current task entropy value and the historical task entropy value divided by the modulus of the two values. Since the entropy value is a scalar, the cosine similarity is simplified to the smaller value of the ratio of the two values to 1. The cosine similarities of the current task entropy value of the downhole electric pump system with each historical task entropy value are 0.561, 0.841, 0.973, 0.856 and 0.888, respectively. The entropy value mode similarity reflects the similarity of the complexity and uncertainty of the current diagnostic task and the historical diagnostic cases.

[0147] The distance similarity and the entropy value mode similarity are fused to obtain a comprehensive similarity. The weight coefficients of the weighted fusion are set by experience, and the weight coefficient of the distance similarity is set to 0.7, and the weight coefficient of the entropy value mode similarity is set to 0.3. After the distance similarity and the entropy value mode similarity of each diagnostic strategy of the downhole electric pump system are fused by weighting, the comprehensive similarity obtained is 0.528, 0.789, 0.930, 0.789 and 0.847 respectively. The comprehensive similarity comprehensively considers the matching degree of the current diagnostic task and the applicable scene of each historical diagnostic strategy, and provides a reliable basis for the selection of the diagnostic strategy.

[0148] The diagnostic strategies are ranked in descending order based on the comprehensive similarity, and the diagnostic strategy ranked first is selected as the current execution strategy. The diagnostic strategies are sorted according to the comprehensive similarity from large to small, and the sorting result is: the causal chain tracking diagnostic strategy (0.930), the regional step-by-step diagnostic strategy (0.847), the multi-component parallel diagnostic strategy (0.789), the full-system scanning diagnostic strategy (0.789), and the single-component deep diagnostic strategy (0.528). Therefore, the causal chain tracking diagnostic strategy is selected as the current execution strategy. The causal chain tracking diagnostic strategy is a method of diagnosing step by step along the fault propagation path, and is suitable for the case where there is a clear propagation link of the fault and there is a strong causal relationship between components. For the downhole electric pump system, the diagnosis starts from the main bearing, traces upwards to the motor rotor, the motor stator and the power supply control unit, and traces downwards to the sealing device, the pump body inlet and other components, and gradually confirms the fault state and the influence degree of each component.

[0149] In this embodiment, the information entropy of the fault probability distribution of each component is modeled, so that the degree of uncertainty contained in the diagnostic result can be quantified and described, and the complexity and information completeness of the current diagnostic task can be more comprehensively reflected. The task entropy and the causal correlation entropy are fused by weighting to form a comprehensive entropy feature vector, which realizes the unified representation of the overall state of the diagnostic task. The similarity between the comprehensive entropy feature of the current diagnostic task and the comprehensive entropy feature corresponding to the historical strategy is measured, which realizes the adaptive strategy matching based on historical experience, improves the accuracy and stability of the strategy recommendation, and sorts the diagnostic strategies by the comprehensive similarity and selects the optimal strategy. Under different fault complexity and causal clarity, the diagnostic strategy with controllable risk and higher information benefit can be automatically matched, which significantly improves the adaptability and decision rationality of the diagnostic process under complex working conditions, and helps to reduce the risk of misdiagnosis and improve the overall diagnostic efficiency.

[0150] In an optional implementation manner,

[0151] The optimized diagnostic result obtained by performing deep time-frequency analysis and causal chain verification on the component with the highest fault probability value based on the current execution strategy includes:

[0152] extracting a component with the highest fault probability value from the diagnostic reasoning result as a target diagnostic component and obtaining corresponding real-time monitoring time series data to perform adaptive wavelet decomposition to obtain frequency band components, calculating instantaneous frequency and instantaneous amplitude of the frequency band components to construct a time-frequency joint representation matrix, calculating a gradient change rate corresponding to the time-frequency joint representation matrix to obtain a time-frequency evolution trajectory and extracting amplitude mutation points and frequency drift points as abnormal time-frequency feature points, performing density distribution calculation and peak value detection on the abnormal time-frequency feature points to obtain an abnormal aggregation time, and intercepting abnormal period data segments based on the abnormal aggregation time;

[0153] extracting parent node components and child node components of the target diagnostic component from the causal directed acyclic graph to construct a causal verification link and obtain corresponding historical monitoring time series data, performing cross-correlation analysis on the historical monitoring time series data and the abnormal period data segments to obtain a time lag correlation coefficient, calculating a difference between the time lag correlation coefficient and a Granger causality coefficient corresponding to the directed edge to obtain a causal bias, and determining a modified causal verification link according to the causal bias and a preset verification threshold;

[0154] calculating a proportion of the number of remaining components in the modified causal verification link to the number of initial components to obtain a causal verification confidence, modifying a fault probability value based on the causal verification confidence to obtain a modified fault probability value, and combining the modified fault probability value and the modified causal verification link to obtain an optimized diagnostic result.

[0155] The component with the highest fault probability value is extracted from the diagnostic reasoning result as the target diagnostic component, and the corresponding real-time monitoring time series data is obtained for adaptive wavelet decomposition to obtain frequency band components. The diagnostic reasoning result contains the fault probability value of each component. By comparing, the component with the highest fault probability value is found as the target diagnostic component. Taking the downhole electric pump system as an example, the diagnostic reasoning result shows that the fault probability value of the main bearing is 0.82, which is higher than that of other components, so the main bearing is selected as the target diagnostic component. The real-time monitoring time series data of the main bearing is obtained, including vibration signals, temperature signals, etc. The sampling frequency of the vibration signal is 1000Hz, the sampling frequency of the temperature signal is 1Hz, and the data duration is 30 minutes. The adaptive wavelet decomposition is performed on the obtained real-time monitoring time series data, and the improved empirical mode decomposition method is used to decompose the signal into multiple frequency band components. The parameter settings of adaptive wavelet decomposition include 5 layers of decomposition, db4 wavelet as the wavelet basis function, soft threshold value as the threshold value function, and maximum minimum threshold criterion as the threshold value selection. After 5-layer wavelet decomposition of the vibration signal of the main bearing, 5 detail components and 1 approximation component are obtained, which correspond to different frequency bands respectively: detail component d1 corresponds to frequency band 500-1000Hz, detail component d2 corresponds to frequency band 250-500Hz, detail component d3 corresponds to frequency band 125-250Hz, detail component d4 corresponds to frequency band 62.5-125Hz, detail component d5 corresponds to frequency band 31.25-62.5Hz, and approximation component a5 corresponds to frequency band 0-31.25Hz.

[0156] The instantaneous frequency and instantaneous amplitude of the frequency band component are calculated to construct a time-frequency joint representation matrix. The instantaneous frequency is calculated using the Hilbert transform method, and the analytic signal is obtained by performing Hilbert transform on each frequency band component. The instantaneous frequency is calculated by calculating the phase derivative of the analytic signal. The instantaneous amplitude is calculated as the modulus of the analytic signal. The instantaneous frequency and instantaneous amplitude of the six frequency band components obtained by decomposing the vibration signal of the main bearing are calculated, and a time-frequency joint representation matrix with dimensions of 6x2x18000 is obtained, where 6 represents the number of frequency band components, 2 represents the feature dimension (frequency and amplitude), and 18000 represents the number of time points (30 minutesx60 secondsx10 sampling points / second). In the d3 frequency band component (125-250Hz), the instantaneous frequency is stable at about 185Hz under normal conditions, and the instantaneous amplitude fluctuates in the range of 0.5-0.8mm / s; while in abnormal time, the instantaneous frequency drifts obviously, reaching a maximum of 220Hz, and the instantaneous amplitude suddenly increases to 2.5mm / s.

[0157] The gradient change rate corresponding to the time-frequency joint representation matrix is calculated to obtain the time-frequency evolution trajectory, and the amplitude mutation point and the frequency drift point are extracted as the abnormal time-frequency feature points. The gradient change rate calculation method is the difference between the instantaneous frequency and the instantaneous amplitude of adjacent time points divided by the time interval. For each frequency band component, the gradient change rates of the instantaneous frequency and the instantaneous amplitude are calculated to form the time-frequency evolution trajectory. The abnormal time-frequency feature point extraction is based on the anomaly detection algorithm, and the threshold is set to be 3 times the standard deviation of the gradient change rate in the normal state. For the d3 frequency band component, the standard deviation of the gradient change rate of the instantaneous frequency in the normal state is 0.5 Hz / s, and the standard deviation of the gradient change rate of the instantaneous amplitude is 0.1 mm / s 2 ; therefore, the judgment threshold of the frequency drift point is 1.5 Hz / s, and the judgment threshold of the amplitude mutation point is 0.3 mm / s 2 . Based on these thresholds, 15 frequency drift points and 12 amplitude mutation points are detected on the d3 frequency band component and marked as abnormal time-frequency feature points.

[0158] The abnormal time-frequency feature points are subjected to density distribution calculation and peak value detection to obtain the abnormal aggregation time. The density distribution calculation adopts the kernel density estimation method, projects each abnormal time-frequency feature point on the time axis and accumulates the weighted sum to form a density curve with time as the horizontal axis. The kernel function is selected as the Gaussian kernel, and the bandwidth parameter is set to 30 seconds. The peak value detection adopts the local maximum value detection algorithm, and the minimum peak value height is set to 2 times the average density, and the minimum peak value distance is 60 seconds. After the density distribution calculation of the abnormal time-frequency feature points of the main bearing, 3 density peaks are detected, and the corresponding time points are the 8th minute, the 15th minute and the 22nd minute after the start of recording. These time points are determined as the abnormal aggregation time.

[0159] Based on the abnormal aggregation time, the abnormal period data segment is intercepted. The interception method is to extend 30 seconds forward and backward from the abnormal aggregation time to form a 1-minute data segment. For the 3 abnormal aggregation times detected, the corresponding data segments are intercepted to obtain 3 abnormal period data segments. The first abnormal period data segment corresponds to the time from 7 minutes and 30 seconds to 8 minutes and 30 seconds, the second abnormal period data segment corresponds to the time from 14 minutes and 30 seconds to 15 minutes and 30 seconds, and the third abnormal period data segment corresponds to the time from 21 minutes and 30 seconds to 22 minutes and 30 seconds. The abnormal period data segment contains the most prominent characteristic information of the fault, and provides key data for subsequent causal verification.

[0160] The parent node components and child node components of the target diagnostic component are extracted from the causal directed acyclic graph to build a causal verification link and obtain the corresponding historical monitoring time series data. The directed edges pointing to the main bearing (target diagnostic component) and the directed edges starting from the main bearing are found in the causal directed acyclic graph to determine the parent node components and child node components of the main bearing. The parent node component of the main bearing is the motor rotor, and the Granger causality coefficient is 0.42; the child node components of the main bearing are the auxiliary bearing (Granger causality coefficient is 0.35) and the sealing device (Granger causality coefficient is 0.39). These components are connected to form a causal verification link: motor rotor → main bearing → auxiliary bearing / sealing device. The corresponding historical monitoring time series data of these components are obtained, including the current and rotational speed data of the motor rotor, the vibration and temperature data of the main bearing, the vibration data of the auxiliary bearing, and the pressure data of the sealing device. The historical data come from normal operation and fault cases in the past 30 days, the data length is 24 hours per day, and the sampling frequency is the same as that of the real-time monitoring data.

[0161] The time-lag correlation coefficient is obtained by cross-correlation analysis of the historical monitoring time series data and the abnormal period data segment. The cross-correlation analysis method is to calculate the correlation coefficient of two signals at different time delays, and find the maximum correlation coefficient and its corresponding time delay. For each pair of adjacent components, the cross-correlation function of the historical data and the abnormal period data segment is calculated, the time delay range is set to -5 seconds to 5 seconds, and the step is 0.1 second. For the motor rotor and the main bearing, the cross-correlation function of the current and vibration signals is calculated, and the maximum correlation coefficient is 0.65, and the corresponding time delay is 0.3 seconds; for the main bearing and the auxiliary bearing, the cross-correlation function of the vibration signals is calculated, and the maximum correlation coefficient is 0.72, and the corresponding time delay is 0.2 seconds; for the main bearing and the sealing device, the cross-correlation function of the vibration and pressure signals is calculated, and the maximum correlation coefficient is 0.48, and the corresponding time delay is 0.5 seconds. These maximum correlation coefficients are recorded as time-lag correlation coefficients.

[0162] The causal deviation is obtained by calculating the difference between the time-lag correlation coefficient and the Granger causality coefficient corresponding to the directed edge. The causal deviation calculation method is to subtract the absolute value of the Granger causality coefficient from the time-lag correlation coefficient. For the link from the motor rotor to the main bearing, the time-lag correlation coefficient is 0.65, the Granger causality coefficient is 0.42, and the causal deviation is 0.23; for the link from the main bearing to the auxiliary bearing, the time-lag correlation coefficient is 0.72, the Granger causality coefficient is 0.35, and the causal deviation is 0.37; for the link from the main bearing to the sealing device, the time-lag correlation coefficient is 0.48, the Granger causality coefficient is 0.39, and the causal deviation is 0.09.

[0163] The modified causal verification link is determined according to the causal deviation quantity and a preset verification threshold. The verification threshold is set to 0.25, that is, when the causal deviation quantity is greater than 0.25, it is considered that the causal relationship has a significant change and needs to be modified; when the causal deviation quantity is less than 0.25, it is considered that the causal relationship is relatively stable and the original link is retained. For the link from the motor rotor to the main bearing, the causal deviation quantity is 0.23, which is less than the threshold 0.25, and the link is retained; for the link from the main bearing to the auxiliary bearing, the causal deviation quantity is 0.37, which is greater than the threshold 0.25, and the link is marked as needing to be modified; for the link from the main bearing to the sealing device, the causal deviation quantity is 0.09, which is less than the threshold 0.25, and the link is retained. The modified causal verification link is: motor rotor→main bearing→sealing device, and the link from the main bearing to the auxiliary bearing is removed.

[0164] The causal verification confidence is obtained by calculating the proportion of the number of remaining components in the modified causal verification link to the number of initial components. The initial causal verification link contains 4 components: motor rotor, main bearing, auxiliary bearing and sealing device; the modified causal verification link contains 3 components: motor rotor, main bearing and sealing device. The proportion of the number of remaining components to the number of initial components is 3 / 4=0.75, that is, the causal verification confidence is 0.75.

[0165] The modified failure probability value is obtained by modifying the failure probability value based on the causal verification confidence. The modification method is to multiply the original failure probability value by the causal verification confidence and add a compensation term. The compensation term is designed as the original failure probability value multiplied by the remaining causal verification confidence (1 minus the causal verification confidence) and multiplied by an adjustment coefficient 0.5. The original failure probability value of the main bearing is 0.82, and the causal verification confidence is 0.75. The calculated modified failure probability value is 0.82×0.75+0.82×(1-0.75)×0.5=0.615+0.1025=0.7175, which is about 0.72.

[0166] The optimized diagnosis result is obtained by combining the modified failure probability value and the modified causal verification link. The optimized diagnosis result includes information such as fault component, modified failure probability value, fault type, fault propagation path, etc. For the downhole electric pump system, the optimized diagnosis result is: there is a wear fault in the inner ring of the main bearing, and the modified failure probability value is 0.72; the fault propagation path is motor rotor→main bearing→sealing device, indicating that the main bearing fault is affected by the running state of the motor rotor, and also affects the sealing device; it is suggested to take measures including replacing the main bearing, checking the running state of the motor rotor, and closely monitoring whether the sealing device is damaged.

[0167] In this embodiment, by performing adaptive wavelet decomposition and time-frequency joint representation on the component with the highest failure probability, the abnormality analysis is improved from the original low-dimensional statistical features to the high-resolution time-frequency evolution level, the recognition sensitivity of sudden and hidden faults is improved, the abnormal period is determined based on the density aggregation of abnormal feature points, the false triggering caused by noise fluctuation is effectively avoided, the abnormal positioning on the time axis is more concentrated and accurate, the time delay correlation coefficient is compared with the Granger causality coefficient, and the causality verification link is corrected accordingly, so that the causality structure can be dynamically checked and corrected according to the real-time abnormal features, and the adaptability of the diagnosis system under the conditions of working condition change and structural uncertainty is significantly improved. The causality verification confidence is calculated based on the size of the corrected causality link, and the confidence is used to modify the failure probability again, so that the final failure probability not only reflects the model inference result, but also integrates the abnormal time-frequency features and the causality consistency verification information, and the reliability and explainability of the diagnosis result are significantly improved.

[0168] In a second aspect of the embodiment of the present application, a downhole equipment failure real-time diagnosis system based on edge computing is provided, comprising:

[0169] A feature extraction unit is configured to collect multi-modal data corresponding to downhole equipment through a distributed sensing network and perform time-frequency domain transformation through an edge computing node, extract frequency spectrum features and phase features of each modal signal under different time scales, and organize multi-scale time sequence features.

[0170] A feature fusion unit is configured to project the multi-scale time sequence features to a Lie group manifold space, calculate the geodesic distance and curvature change rate of each modal feature in the manifold space, and construct a coupling mapping relationship matrix between the modal features, calculate the fusion weight of each modal feature according to the coupling mapping relationship matrix, and weightedly sum to obtain a fusion feature.

[0171] A graph construction unit is configured to obtain the topological connection relationship of downhole equipment and historical failure sample data, perform time sequence Granger causality test on the historical failure sample data, calculate the Granger causality coefficient between the time sequence data of each component, and construct a causality directed acyclic graph based on the topological connection relationship and the Granger causality coefficient.

[0172] A fault diagnosis unit is configured to input the fusion feature into each node of the causality directed acyclic graph and perform Bayesian probability reasoning along the direction of the directed edge, iteratively update the failure probability value of each node, and solve to obtain a diagnosis reasoning result, calculate the entropy value of the failure probability distribution of each component in the diagnosis reasoning result, and perform similarity matching with the historical task entropy value corresponding to the preset diagnosis strategy to determine the current execution strategy, and perform deep time-frequency analysis and causality chain verification on the component with the highest failure probability based on the current execution strategy to obtain an optimized diagnosis result.

[0173] In a third aspect, the present application provides an electronic device, comprising:

[0174] A processor and a memory for storing processor-executable instructions, wherein the processor is configured to invoke the instructions stored by the memory to perform the method described above.

[0175] In a fourth aspect, the present application provides a computer-readable storage medium having stored thereon computer program instructions, which when executed by a processor, implement the method described above.

[0176] The present application can be a method, apparatus, system and / or computer program product. The computer program product can include a computer-readable storage medium having computer-readable program instructions stored therein, which are used to perform various aspects of the present application.

[0177] Finally, it should be noted that: the above embodiments are only used to illustrate the technical solutions of the present application, and not to limit them; although the present application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand: it can still modify the technical solutions recorded in the foregoing embodiments, or make equivalent replacement for part or all of the technical features; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the scope of the technical solutions of the embodiments of the present application.

Claims

1. A real-time fault diagnosis method for downhole equipment based on edge computing, characterized in that, include: Multimodal data corresponding to downhole equipment is collected through a distributed sensor network and time-frequency domain transformation is performed through edge computing nodes. The spectral and phase features of each modal signal at different time scales are extracted and organized to obtain multi-scale time series features. The multi-scale temporal features are projected onto the Lie group manifold space. The geodesic distance and rate of curvature change of each modal feature in the manifold space are calculated, and the coupling mapping relationship matrix between each modal feature is constructed. The fusion weight of each modal feature is calculated based on the coupling mapping relationship matrix, and the weighted sum is obtained to obtain the fused feature. Obtain the topological connection relationship and historical fault sample data of the downhole equipment, perform time-series Granger causality test on the historical fault sample data and calculate the Granger causality coefficient between the time-series data of each component, and construct a causal directed acyclic graph based on the topological connection relationship and the Granger causality coefficient. The fused features are input into each node of the causal directed acyclic graph and Bayesian probabilistic inference is performed along the directed edge direction. The fault probability value of each node is iteratively updated and the diagnostic inference result is obtained. The entropy value of the fault probability distribution of each component in the diagnostic inference result is calculated and matched with the entropy value of the historical task corresponding to the preset diagnostic strategy to determine the current execution strategy. Based on the current execution strategy, deep time-frequency analysis and causal chain verification are performed on the component with the highest fault probability value to obtain the optimized diagnostic result.

2. The method according to claim 1, characterized in that, Multimodal data corresponding to downhole equipment is collected through a distributed sensor network and time-frequency domain transformation is performed through edge computing nodes. The spectral and phase features of each modal signal at different time scales are extracted and organized to obtain a multi-scale time-series feature tensor, including: Multimodal data is obtained by collecting current and voltage signal sequences and temperature and pressure signal sequences corresponding to downhole equipment through a distributed sensor network and aligning them according to timestamps. The multimodal data is sent to the edge computing node, and multiple time windows of different lengths are set on the edge computing node. The frequency domain complex spectrum of the current and voltage signal sequence is obtained by performing Fourier transform on each time window, and the modulus is calculated to obtain the spectral energy distribution. The phase difference of the harmonic components in the spectral energy distribution is extracted to obtain the phase shift feature. The temperature and pressure signal sequence is subjected to Hilbert-Huang transform within each time window. The intrinsic mode function is obtained through empirical mode decomposition and Hilbert transform is performed to obtain the instantaneous frequency distribution. The time derivative of the instantaneous phase is calculated based on the instantaneous frequency distribution to obtain the phase evolution characteristics. The electrical mode feature sequence is obtained by arranging the spectral energy distribution and the phase shift feature according to the time window length. The thermodynamic mode feature sequence is obtained by arranging the instantaneous frequency distribution and the phase evolution feature according to the time window length. The multi-scale time series feature is obtained by organizing the electrical mode feature sequence and the thermodynamic mode feature sequence according to the mode type.

3. The method according to claim 1, characterized in that, The multi-scale temporal features are projected onto the Lie group manifold space. The geodesic distance and rate of curvature change of each modal feature in the manifold space are calculated, and a coupling mapping matrix between each modal feature is constructed. Based on the coupling mapping matrix, the fusion weights of each modal feature are calculated and weighted summation is performed to obtain the fused features, including: The multi-scale temporal features are nonlinearly mapped to obtain manifold points in the Lie group manifold space. A tangent space is constructed at the manifold points and the basis vectors of the tangent space are calculated. The manifold points corresponding to two modal features are selected as the starting point and the ending point, respectively. A direction vector is set in the tangent space of the starting point based on the basis vector. Step movement is performed until the ending point is reached, and the step path is integrated by Riemann metric to obtain the geodesic distance. The local neighborhood of the manifold point is sampled to obtain neighborhood sampling points, and the tangent vector corresponding to the neighborhood sampling points is calculated based on the basis vector. The tangent vector of the neighborhood sampling points is transmitted along a closed path and compared with the initial tangent vector corresponding to the starting point to obtain the deviation vector. The curvature tensor components are calculated based on the deviation vector and the scalar curvature is obtained through tensor contraction operation. The scalar curvature is then subjected to finite difference to obtain the rate of change of curvature. Geometric similarity is obtained by mapping the geodesic distance using a preset exponential kernel function. The difference in the rate of curvature change between different modal features is calculated and mapped using a reciprocal function to obtain the evolutionary consistency. The geometric similarity and the evolutionary consistency are multiplied element-wise to obtain the coupling strength and then organized into a coupling mapping relationship matrix. The row elements in the coupling mapping relationship matrix are normalized to obtain the fusion weights, and the modal features are weighted and fused to obtain the fusion features.

4. The method according to claim 1, characterized in that, Acquire the topological connectivity of downhole equipment and historical fault sample data; perform a time-series Granger causality test on the historical fault sample data and calculate the Granger causality coefficients between the time-series data of each component; construct a causal directed acyclic graph based on the topological connectivity and the Granger causality coefficients, including: Obtain the physical connection information between the components of the downhole equipment and establish the topological connection relationship according to the connection direction. Obtain the monitoring data of each component before and after the occurrence of historical failures as historical failure sample data. The time-series monitoring data of each component in the historical fault sample data are subjected to stationarity test. The time-series monitoring data that fails the stationarity test are differentially processed to obtain differential data. The differential data is then merged with the time-series monitoring data that passes the stationarity test to obtain stationary time-series data. Stationary time series data corresponding to any two components are selected as the target sequence and reference sequence, respectively. The reference sequence is shifted according to different time lag orders to obtain a lag sequence group. Multiple linear regression is performed on the target sequence and the lag sequence group to obtain the regression coefficients corresponding to each lag order. The causal significance statistic is obtained based on the regression coefficients and the lag order. The causal significance statistic is compared with a preset critical value to obtain the Granger causality test result. The Granger causality coefficient is calculated based on the regression coefficients that pass the significance test in the Granger causality test result. Each component in the topological connection relationship is taken as a node. Directed edges are established between components whose Granger causality coefficient is greater than a preset threshold, and loop detection is performed on the directed edges. The directed edges that form loops are deleted to obtain a causal directed acyclic graph.

5. The method according to claim 1, characterized in that, The fused features are input into each node of the causal directed acyclic graph, and Bayesian probabilistic inference is performed along the directed edge direction. The fault probability value of each node is iteratively updated, and the diagnostic inference result is obtained by solving the problem, including: According to the component type, the fused features are assigned to the causal directed acyclic graph and multi-dimensional feature decoupling is performed to obtain state degradation feature vector and dynamic evolution feature vector. Based on the variational inference algorithm and the state degradation feature vector, the posterior distribution parameters of the fault state are calculated and the initial fault probability value is sampled. Starting from the root node of the causal directed acyclic graph, the initial failure probability value of the current node, the failure probability value of the corresponding parent node, and the dynamic evolution feature vector are obtained. A time-varying conditional probability transfer kernel is constructed based on the Granger causality coefficient of the directed edge and the dynamic evolution feature vector, and tensor convolution is performed with the failure probability value to obtain the causal transmission probability distribution. The causal transmission probability distribution is marginalized and integrated to obtain the aggregated prior probability. Based on the aggregated prior probability and the initial failure probability value, the updated failure probability value is obtained by optimization and passed to the child node with the dynamic evolution feature vector. The traversal and update are repeated until a single iteration is formed. Calculate the KL divergence of the fault probability values ​​of each node before and after a single iteration and compare it with a preset convergence threshold. If it is less than the convergence threshold, terminate the iteration. Extract the fault probability values ​​of each node and the dynamic evolution feature vector to calculate the comprehensive fault score. Take the component corresponding to the node with the largest comprehensive fault score as the fault source and integrate it to obtain the diagnostic reasoning result.

6. The method according to claim 1, characterized in that, Calculating the entropy value of the failure probability distribution of each component in the diagnostic inference result and performing similarity matching with the entropy value of the historical task corresponding to the preset diagnostic strategy to determine the current execution strategy includes: The fault probability values ​​of each component are extracted from the diagnostic reasoning results and normalized to obtain the fault probability distribution. The information entropy of the fault probability distribution is calculated to obtain the current task entropy value. The causal transmission probability distribution between each component in the diagnostic reasoning results is extracted and the conditional entropy is calculated to obtain the causal association entropy value. The current task entropy value and the causal association entropy value are weighted and summed to obtain the comprehensive entropy feature vector. Obtain the historical task entropy value and historical causal association entropy value corresponding to each diagnostic strategy in the preset diagnostic strategy library, and sum them according to the preset weights to obtain the historical comprehensive entropy feature vector. The distance similarity is obtained by calculating the Euclidean distance between the comprehensive entropy feature vector and the historical comprehensive entropy feature vector and performing a reciprocal transformation. The entropy pattern similarity is obtained by calculating the cosine similarity between the current task entropy value and the historical task entropy value. The comprehensive similarity is obtained by weighted fusion of the distance similarity and the entropy pattern similarity. The diagnostic strategies are sorted in descending order based on the comprehensive similarity, and the diagnostic strategy ranked first is taken as the current execution strategy.

7. The method according to claim 1, characterized in that, Based on the current execution strategy, in-depth time-frequency analysis and causal chain verification are performed on the component with the highest failure probability value to obtain optimized diagnostic results, including: The component with the highest fault probability value is extracted from the diagnostic inference results as the target diagnostic component, and the corresponding real-time monitoring time series data is obtained and adaptive wavelet decomposition is performed to obtain the frequency band component. The instantaneous frequency and instantaneous amplitude of the frequency band component are calculated to construct the time-frequency joint characterization matrix. The gradient change rate corresponding to the time-frequency joint characterization matrix is ​​calculated to obtain the time-frequency evolution trajectory, and amplitude mutation points and frequency drift points are extracted as abnormal time-frequency feature points. The density distribution of the abnormal time-frequency feature points is calculated and peak detection is performed to obtain the abnormal clustering time. Based on the abnormal clustering time, abnormal time period data segments are extracted. The parent and child nodes of the target diagnostic component are extracted from the causal directed acyclic graph to construct a causal verification link and obtain the corresponding historical monitoring time series data. Cross-correlation analysis is performed on the historical monitoring time series data and the abnormal period data segments to obtain the time-delay correlation coefficient. The difference between the time-delay correlation coefficient and the Granger causality coefficient corresponding to the directed edge is calculated to obtain the causal deviation. The corrected causal verification link is determined based on the causal deviation and the preset verification threshold. The causal verification confidence score is obtained by calculating the ratio of the number of remaining components in the corrected causal verification link to the number of initial components. The fault probability value is corrected based on the causal verification confidence score to obtain the corrected fault probability value. The corrected fault probability value is combined with the corrected causal verification link to obtain the optimized diagnostic result.

8. A real-time fault diagnosis system for downhole equipment based on edge computing, used to implement the method of any one of claims 1-7, characterized in that, include: The feature extraction unit is used to collect multimodal data corresponding to downhole equipment through a distributed sensor network and perform time-frequency domain transformation through edge computing nodes to extract the spectral and phase features of each modal signal at different time scales and organize them to obtain multi-scale time-series features. The feature fusion unit is used to project the multi-scale temporal features onto the Lie group manifold space, calculate the geodesic distance and curvature change rate of each modal feature in the manifold space, construct the coupling mapping relationship matrix between each modal feature, calculate the fusion weight of each modal feature according to the coupling mapping relationship matrix, and obtain the fused feature by weighted summation. The graph construction unit is used to acquire the topological connection relationship and historical fault sample data of downhole equipment, perform time-series Granger causality test on the historical fault sample data and calculate the Granger causality coefficient between the time-series data of each component, and construct a causal directed acyclic graph based on the topological connection relationship and the Granger causality coefficient. The fault diagnosis unit is used to input the fused features into each node of the causal directed acyclic graph and perform Bayesian probabilistic inference along the directed edge direction, iteratively update the fault probability value of each node and solve to obtain the diagnosis inference result, calculate the entropy value of the fault probability distribution of each component in the diagnosis inference result and perform similarity matching with the historical task entropy value corresponding to the preset diagnosis strategy to determine the current execution strategy, and perform deep time-frequency analysis and causal chain verification on the component with the highest fault probability value based on the current execution strategy to obtain the optimized diagnosis result.

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

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

Images