Underground equipment fault real-time diagnosis method and system based on edge calculation

Abstract

The invention provides an underground equipment fault real-time diagnosis method and system based on edge calculation, and relates to the technical field of coal mine safety production, and the method comprises the steps: collecting multi-modal data through a distributed sensor network, extracting multi-scale time sequence features, projecting the features to a Lie group manifold space, constructing a coupling mapping relation matrix, obtaining fusion features, and carrying out the real-time diagnosis of an underground equipment fault; and constructing a causal directed acyclic graph based on a topological connection relationship and a Granger causal coefficient, executing Bayesian probabilistic reasoning, determining an execution strategy in combination with entropy similarity matching, and performing deep time-frequency analysis and causal chain verification. High-precision real-time diagnosis of equipment faults in an underground complex environment is realized, and the fault early warning accuracy is improved.

Description

Technical Field

[0001] This invention relates to the field of coal mine safety production technology, and in particular to a method and system for real-time fault diagnosis of underground equipment based on edge computing. Background Technology

[0002] With the deepening of intelligent construction in coal mines, the safe and stable operation of underground equipment is crucial to ensuring mine production. The underground environment is complex and ever-changing, posing numerous challenges to equipment operation status monitoring and fault diagnosis. Traditional underground equipment fault diagnosis mainly relies on manual experience or single sensor data analysis;

[0003] With the development of IoT and edge computing technologies, fault diagnosis methods based on multi-source data have been widely used. By deploying distributed sensor networks to collect multimodal data and using data analysis techniques to extract fault features and identify patterns, the introduction of edge computing technology allows data processing to be performed close to the data source, effectively reducing data transmission latency and improving real-time performance.

[0004] However, existing downhole equipment fault diagnosis technologies still suffer from several problems, including a lack of in-depth exploration of the intrinsic correlations between different modal data, resulting in limited expressive power of fused features, inability to fully reflect the complex characteristics of equipment faults, neglect of the physical connections and causal effects between equipment components, lack of interpretability of diagnostic results, difficulty in tracing fault propagation paths, lack of adaptive diagnostic strategy selection mechanisms, inability to dynamically adjust the analysis depth according to different fault types and the uncertainty of diagnostic results, and difficulty in achieving efficient and accurate real-time diagnosis in resource-constrained edge computing environments. Summary of the Invention

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

[0006] A first aspect of this invention provides a real-time fault diagnosis method for downhole equipment based on edge computing, comprising:

[0007] 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.

[0008] 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.

[0009] 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.

[0010] 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.

[0011] In one alternative implementation,

[0012] 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:

[0013] 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.

[0014] 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.

[0015] 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.

[0016] 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.

[0017] In one alternative implementation,

[0018] 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:

[0019] 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.

[0020] 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.

[0021] 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.

[0022] In one alternative implementation,

[0023] 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:

[0024] 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.

[0025] 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.

[0026] 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.

[0027] 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.

[0028] In one alternative implementation,

[0029] 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:

[0030] 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.

[0031] 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.

[0032] 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.

[0033] In one alternative implementation,

[0034] 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:

[0035] 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.

[0036] 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.

[0037] 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.

[0038] In one alternative implementation,

[0039] 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:

[0040] 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.

[0041] 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.

[0042] 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.

[0043] A second aspect of this invention provides a real-time fault diagnosis system for downhole equipment based on edge computing, comprising:

[0044] 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.

[0045] 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.

[0046] 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.

[0047] 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.

[0048] A third aspect of the present invention provides an electronic device, comprising:

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

[0050] A fourth aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0051] In this invention, multi-scale temporal features are projected onto the Lie group manifold space, and geodesic distance and curvature change rate are calculated. A coupling mapping relationship matrix is ​​constructed for feature fusion, which effectively captures the nonlinear relationship between different modal signals and improves the feature representation capability. A causal directed acyclic graph is constructed through temporal Granger causality test and topological connectivity relationship, avoiding spurious associations that may be caused by traditional correlation analysis, and providing an interpretable causal link for fault diagnosis. The adaptive selection mechanism of diagnostic strategy based on Bayesian probabilistic inference and entropy calculation realizes intelligent response to complex fault scenarios, which significantly improves the accuracy and efficiency of diagnosis. The optimized diagnostic process through deep time-frequency analysis and causal chain verification improves the accuracy of fault location, reduces false alarms and false negatives, and provides a reliable basis for downhole equipment maintenance decisions. Attached Figure Description

[0052] Figure 1 This is a flowchart illustrating the real-time fault diagnosis method for downhole equipment based on edge computing, according to an embodiment of the present invention.

[0053] Figure 2 This is a flowchart illustrating the intelligent fault source diagnosis process of the real-time fault diagnosis method for downhole equipment based on edge computing, as described in an embodiment of the present invention. Detailed Implementation

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

[0055] The technical solution of the present invention will be described in detail below with reference to specific embodiments. These specific embodiments can be combined with each other, and the same or similar concepts or processes may not be described again in some embodiments.

[0056] Figure 1 This is a flowchart illustrating the real-time fault diagnosis method for downhole equipment based on edge computing, as described in an embodiment of the present invention. Figure 1As shown, the method includes:

[0057] 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.

[0058] 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.

[0059] 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.

[0060] 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.

[0061] In one alternative implementation,

[0062] 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:

[0063] 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.

[0064] 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.

[0065] 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.

[0066] 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.

[0067] A distributed sensor network is used to collect current, voltage, temperature, and pressure signal sequences corresponding to downhole equipment. This network comprises multiple current, voltage, temperature, and pressure sensors, distributed across key components of the downhole equipment, such as motor windings, bearings, and the pump body. The sampling frequency for the current and voltage sensors is set to 10kHz, and the sampling frequency for the temperature and pressure sensors is set to 1kHz. Each collected data point includes a precise timestamp for subsequent time-series alignment. For example, for a downhole electric pump, the current sensor collects three-phase current values; the voltage sensor collects three-phase voltage values; the temperature sensor collects motor and pump body temperatures; and the pressure sensor collects inlet and outlet pressures.

[0068] Time alignment is performed based on timestamps to align data collected from different sensors to a unified timeline. During processing, interpolation is applied to data with different sampling frequencies to ensure that all data are precisely consistent in the time dimension. Temperature and pressure data are improved to the same time accuracy as current and voltage data through linear interpolation, generating an aligned multimodal dataset. The time alignment accuracy is controlled within 0.1ms to ensure the temporal correlation of the data.

[0069] Multimodal data is transmitted to edge computing nodes at the wellhead or well site via industrial Ethernet or wireless communication networks. These edge computing nodes are industrial-grade devices with sufficient computing power, featuring processors with a clock speed of at least 2.5 GHz and memory capacity of at least 8 GB to meet real-time computing requirements. Multiple time windows of varying lengths (0.5 seconds, 1 second, 2 seconds, 5 seconds, and 10 seconds) are set on the edge computing nodes for data analysis, covering different timescales from transient to steady-state characteristics.

[0070] Fourier transforms are performed on the current and voltage signal sequences within each time window. For the current or voltage signal within each window, its complex spectrum in the frequency domain is calculated. Taking the first phase current within a 0.5-second window as an example, the time-domain signal is converted to a frequency-domain representation using a Fast Fourier Transform algorithm, obtaining a complex spectrum containing real and imaginary parts. The modulus of the complex spectrum is calculated, and the square root of the sum of the squares of the real and imaginary parts is calculated to obtain the spectral energy distribution. For electrical equipment with a power frequency of 50Hz, the energy values ​​of the fundamental frequency (50Hz), second harmonic (100Hz), third harmonic (150Hz), fifth harmonic (250Hz), and seventh harmonic (350Hz) are extracted and recorded as fundamental energy, second harmonic energy, third harmonic energy, etc., respectively.

[0071] The phase difference of harmonic components is extracted from the spectral energy distribution. The phase value is calculated by taking the arctangent of the ratio of the imaginary part to the real part of the complex spectrum. Characteristics such as the phase difference between the second harmonic and the fundamental wave, the phase difference between the third harmonic and the fundamental wave, and the phase difference between the fifth harmonic and the fundamental wave are calculated, and these phase difference characteristics constitute a set of phase shift characteristics. Normally operating equipment exhibits a stable phase difference value, while the phase difference changes significantly under fault conditions. For example, a bearing failure may cause the phase difference between the second harmonic and the fundamental wave to increase from the normal 20 degrees to more than 35 degrees.

[0072] The temperature and pressure signal sequences are subjected to Hilbert-Huang transform within each time window. An empirical mode decomposition algorithm is then used to decompose the temperature or pressure signal into several intrinsic mode functions (EMFs) and a residual term. For pump body temperature and pressure signals, typically 4 to 6 EMF components can be decomposed, each representing the temperature and pressure variation characteristics at different frequency scales. A Hilbert transform is then performed on each EMF to construct an analytic signal containing the original EMFs and their Hilbert transform results.

[0073] The instantaneous frequency distribution is calculated based on the analytic signal obtained from the Hilbert transform. The instantaneous frequency is defined as the derivative of the instantaneous phase with respect to time, where the instantaneous phase is the arctangent of the ratio of the Hilbert transform result to the original intrinsic mode function. In practical calculations, the derivative is approximated by dividing the phase difference between adjacent time points by the time interval. For each intrinsic mode function component, the corresponding instantaneous frequency sequence is calculated, and the average frequency, standard deviation of frequency, peak frequency, and other characteristics within the window are statistically analyzed.

[0074] The time derivative of the instantaneous phase is calculated based on the instantaneous frequency distribution to obtain the phase evolution characteristics. The first-order and second-order difference sequences of the instantaneous phase are calculated, and the statistical characteristics such as the mean, variance, skewness, and kurtosis of the difference sequences are statistically analyzed. The phase evolution of normally operating equipment usually exhibits a stationary change, while faulty equipment may show abrupt or nonlinear change patterns.

[0075] The spectral energy distribution and phase shift characteristics are arranged in ascending order of time window length to form an electrical mode characteristic sequence. For example, for the first phase current, the corresponding characteristic sequence includes the spectral energy distribution and phase shift characteristics for a 0.5-second window, the spectral energy distribution and phase shift characteristics for a 1-second window, the characteristics for a 2-second window, and so on up to the characteristics for a 10-second window.

[0076] The instantaneous frequency distribution and phase evolution characteristics are arranged in ascending order of time window length to form a thermodynamic mode characteristic sequence. Taking the pump body temperature and pressure signals as an example, the characteristic sequence includes the instantaneous frequency distribution and phase evolution characteristics of the 0.5-second window, the characteristics of the 1-second window, the characteristics of the 2-second window, and so on up to the characteristics of the 10-second window.

[0077] The electrical and thermodynamic modal feature sequences are organized according to modal type to obtain multi-scale time series features. Specifically, all current-related features are grouped into a current mode feature set, all voltage-related features are grouped into a voltage mode feature set, all temperature-related features are grouped into a temperature mode feature set, and all pressure-related features are grouped into a pressure mode feature set.

[0078] In this embodiment, by uniformly acquiring heterogeneous signals such as current, voltage, temperature, and pressure at the distributed sensor network level and performing strict time-series alignment based on timestamps, joint characterization of different physical quantities under the same time reference is achieved, improving the accuracy and consistency of multi-source information fusion. By setting time windows of various lengths and performing frequency domain analysis on electrical signals, it is possible not only to capture the main energy distribution characteristics under steady-state operation, but also to identify harmonic changes caused by short-term disturbances and abnormal operating conditions, thereby improving the sensitivity of early fault identification. 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, more accurately depicting the dynamic evolution process of the equipment's thermodynamic state and enhancing the ability to characterize gradual anomalies under complex operating conditions.

[0079] In one alternative implementation,

[0080] 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:

[0081] 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.

[0082] 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.

[0083] 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.

[0084] Nonlinear mapping of multi-scale temporal features yields manifold points in the Lie group manifold space. This nonlinear mapping is achieved through an exponential function, mapping the multi-scale temporal features from the original feature space to the manifold space. First, the multi-scale temporal features are normalized to ensure that all features are within the same order of magnitude. For current mode features, normalization involves subtracting the mean from the original feature value and then dividing by the standard deviation. For temperature and pressure mode features, a maximum-minimum normalization method is used: subtracting the minimum value from the original value and then dividing by the difference between the maximum and minimum values. The normalized features are then transformed to the Lie group manifold space through an exponential mapping function. The base of the exponent in the mapping function is set to the natural constant e, and the exponent 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, it is set to 0.75; for temperature features, it is set to 0.9; and for pressure features, it is set to 0.85.

[0085] A tangent space is constructed at each point of the manifold, and its basis vectors are calculated. During the tangent space construction, the dimension of the manifold is first determined; in this embodiment, the manifold dimension is set to 6. The basis vectors in the tangent space are calculated using the orthogonal basis method, obtained through the Gram-Schmidt orthogonalization process. During the basis vector calculation, a random vector set is generated, and orthogonalization is achieved by eliminating components through iterative projection. Finally, the orthogonalized vectors are normalized to obtain a set of normalized orthogonal basis vectors. Specifically, the random seed is set to 42, initial random vectors are generated, the number of orthogonalization iterations is set to 10, and the vector normalization precision is controlled within 0.0001.

[0086] Two manifold points corresponding to different modal features are selected as the start and end points. A direction vector is set in the tangent space of the start point. Stepping movement is performed until the end point is reached, and the geodesic distance is obtained by integrating the Riemannian metric along the stepping path. The direction vector is set based on the basis vectors of the tangent space of the start point. The vector from the start point to the end point is projected onto each basis vector in the tangent space of the start point, obtaining the projection coefficients of each basis vector. The basis vectors are then weighted and summed using these projection coefficients to obtain the initial direction vector. The stepping movement uses the fourth-order Runge-Kutta method with a step size of 0.05 and an upper limit of 500 iterations. In each iteration, the tangent vector is transmitted along the path via parallel transmission to ensure that the movement direction always lies within the tangent space of the manifold. The Riemannian metric integration uses the Simpson integral method, with the integration interval evenly divided into 100 sub-intervals. For example, the geodesic distance between the corresponding manifold points for the current mode feature and the temperature mode feature is calculated to be 3.28.

[0087] Neighborhood sampling points are obtained by sampling the local neighborhood of the manifold point. A spherical uniform sampling method is used, generating uniformly distributed sampling points on a sphere centered on the manifold point. The sampling radius is set to 0.2, and the number of sampling points is 12, distributed across three orthogonal planes, with 4 points evenly distributed on each plane. For manifold points corresponding to current mode characteristics, after generating the corresponding neighborhood sampling points, the actual distance from each sampling point to the center point is calculated to ensure sampling quality; the distance deviation is required to be no more than 0.01.

[0088] The tangent vector corresponding to the neighborhood sampling point is calculated based on the basis vectors. The method for calculating the tangent vector is to project the vector difference between the sampling point and the center point onto the basis vectors of the tangent space of the center point to obtain the projection coefficients of each basis vector. The basis vectors are then weighted and summed using the projection coefficients of each basis vector as weights to obtain the tangent vector. For example, for the manifold point corresponding to the current mode feature and its neighborhood sampling points, the calculated tangent vector length is between 0.18 and 0.22.

[0089] The deviation vector is obtained by transmitting the tangent vector of the neighboring sampling points along a closed path and comparing it with the initial tangent vector corresponding to the starting point. The closed path is set to start from the starting point, pass through two adjacent neighboring sampling points, and return to the starting point, forming a triangular closed path. Parallel transmission along the closed path adopts an discrete stepping method with a step size of 0.01, and a local tangent space coordinate transformation is performed at each step point. For the manifold point corresponding to the current mode characteristic, after its tangent vector is transmitted around the closed path, the magnitude of the deviation vector from the initial tangent vector is approximately 0.032.

[0090] The curvature tensor components are calculated based on the deviation vector, and the scalar curvature is obtained through tensor contraction. The curvature tensor is calculated by dividing the deviation vector by the area enclosed by the closed path, yielding an approximate value. The tensor contraction operation uses a trace-finding method, which calculates the sum of the diagonal elements of the curvature tensor to obtain the scalar curvature. For the manifold points corresponding to the current mode characteristics, the calculated scalar curvature is approximately 0.57.

[0091] The rate of change of curvature is obtained by performing finite difference on the scalar curvature. The finite difference method employs a central difference scheme. Within a local region of a manifold point, the scalar curvature values ​​of two points on the left and right are taken along a specific direction, and the difference is calculated and divided by the distance to obtain the rate of change of curvature. The difference step size is set to 0.1, and the difference direction is chosen to be the direction of the most significant change in scalar curvature. For the manifold point corresponding to the current mode characteristic, the calculated rate of change of curvature is approximately 1.25.

[0092] Geometric similarity is obtained by mapping geodesic distance using a preset exponential kernel function. The exponential kernel function is in the form of a negative power of the exponent base e, where the power is the ratio of geodesic distance to the kernel width parameter. The kernel width parameter is set to 2.0. For current mode features and temperature mode features, the geometric similarity calculated based on a geodesic distance of 3.28 is approximately 0.19.

[0093] The difference in the rate of curvature change between different modal features is calculated, and the evolutionary consistency is obtained by mapping using a reciprocal function. The difference in the rate of curvature change is calculated as the absolute difference between the rates of curvature change at corresponding manifold points of the two modal features. The reciprocal function mapping uses the reciprocal plus one form to ensure the mapping result is between 0 and 1. For the current modal feature and the temperature modal feature, the rates of curvature change are 1.25 and 2.15, respectively, with a difference of 0.9, and the resulting evolutionary consistency is approximately 0.53.

[0094] The coupling strength is obtained by performing an element-wise product operation on geometric similarity and evolutionary consistency. The element-wise product operation directly multiplies the geometric similarity and evolutionary consistency to obtain a joint metric as the coupling strength. For current mode characteristics and temperature mode characteristics, the coupling strength calculated based on a geometric similarity of 0.19 and an evolutionary consistency of 0.53 is approximately 0.10.

[0095] The coupling strength between all modal feature pairs is organized into a coupling mapping matrix. The matrix rows and columns are labeled with the respective modal features, and the matrix element values ​​are the coupling strength between the corresponding modal features. The diagonal elements of the matrix are set to 1, indicating that the mode is fully coupled with itself. In this embodiment, the four modal features (current, voltage, temperature, and pressure) constitute a 4×4 coupling mapping matrix.

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

[0097] The modal features are weighted and fused to obtain the fused features. During the weighted fusion process, each modal feature is summed according to its corresponding fusion weight to obtain the fused feature representation. In this embodiment, the dimension of the fused feature is the same as the dimension 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.45×1.2+0.23×0.8+0.19×0.6+0.13×0.5=0.895.

[0098] In this embodiment, by performing nonlinear mapping on multi-scale temporal features and characterizing 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, thus improving the accuracy of characterizing 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 characterized, enhancing the sensitivity to gradual and hidden state changes. By coupling geometric similarity with evolutionary consistency and adaptively generating fusion weights based on the coupling strength, dynamic adjustment of the contribution degree of different modes is realized, improving the discriminative ability and stability of fused features.

[0099] In one alternative implementation,

[0100] 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:

[0101] 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.

[0102] 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.

[0103] 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.

[0104] 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.

[0105] The process involves acquiring physical connection information between the components of downhole equipment and establishing topological connections based on the connection directions. Downhole equipment typically consists of multiple components such as motors, pump bodies, bearings, and sealing devices, with clearly defined physical connections between them. Physical connection information, including the relative positions and connection methods of each component, is obtained through equipment structure diagrams and on-site installation records. The connection directions between components are determined based on energy flow and signal transmission directions. For example, if a motor drives a bearing, the connection direction is from the motor to the bearing; if a bearing supports a pump body, the connection direction is from the bearing to the pump body. Taking a downhole electric pump as an example, its components include nine key components: a power control unit, a motor stator, a motor rotor, a main bearing, an auxiliary bearing, a sealing device, a pump inlet, a pump impeller, and a pump outlet. The power control unit is connected to the motor stator, which drives the motor rotor to rotate. The motor rotor is connected to the sealing device through the main and auxiliary bearings. The sealing device is connected to the pump inlet, which is connected to the pump impeller, and the pump impeller is connected to the pump outlet, forming a complete topological connection.

[0106] Monitoring data of each component before and after historical failures were acquired as historical failure sample data. Monitoring data was collected from various sensors installed on downhole equipment, including current sensors, voltage sensors, temperature sensors, and vibration sensors. For historical failure events, monitoring data from 7 days before the failure to 1 day after the failure were extracted, with sampling frequencies 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. Significant outliers, such as values ​​exceeding three standard deviations of the normal operating range, were removed during data preprocessing. Missing data points were filled in using linear interpolation. Taking bearing failure as an example, temperature and vibration data of the main bearing, current data of the motor stator, and pressure data at the pump inlet were collected before and after the failure. In the 4 days before the failure, the main bearing temperature gradually increased from the normal 65°C to 80°C, the vibration amplitude increased from 0.5 mm / s to 2.8 mm / s, and the motor current fluctuation increased by 20%.

[0107] The stationarity of time-series monitoring data for each component in the historical fault sample data was tested. The augmented Dickey-Fowler test was used to test each type of monitoring data sequence for each component separately. The significance level was set to 0.05, and the lag order was automatically selected based on the Schwarz information criterion, with a maximum lag order of 12. After performing the stationarity test on each time-series data, probability values ​​were obtained. These probability values ​​were compared with the significance level; if the probability value was less than 0.05, the stationarity test was passed; otherwise, differencing was required. The stationarity test on the main bearing temperature data yielded a probability value of 0.32, which is greater than the significance level of 0.05, therefore the time series was judged to be non-stationary. The test on the motor stator current data yielded a probability value of 0.02, which is less than the significance level, therefore the time series was judged to be stationary.

[0108] Differential data are obtained by differencing the time-series monitoring data that failed the stationarity test. Differentiation is achieved by calculating the difference between data points at adjacent time points, i.e., the value at the current time point minus the value at the previous time point. The order of differencing is determined based on the data characteristics. Generally, first-order differencing is performed first. If the stationarity test is still not passed after first-order differencing, second-order differencing is performed, up to a maximum of third-order differencing. After first-order differencing of the main bearing temperature data, a stationarity test was performed again, yielding a probability value of 0.03, which is less than the significance level of 0.05. Therefore, the stationarity test was passed, and the first-order differencing result was adopted. The differrated main bearing temperature data reflects the rate of temperature change, with a data range from -0.5°C / hour to 1.8°C / hour.

[0109] The differential data is merged with the time-series monitoring data that has passed the stationarity test to obtain stationary time-series data. During the merging process, the temporal correspondence between the data sequences is maintained; for differentially processed data, the processed sequence is used to replace the original sequence. The merged stationary time-series data includes both the stationary original data and the differentially processed data, together forming the stationary dataset. The main bearing temperature data after first-order differencing, together with the originally stationary motor stator current data, constitutes the stationary time-series dataset for the corresponding component.

[0110] Stationary time series data corresponding to any two components are selected as the target sequence and reference sequence, respectively. Taking the main bearing and motor stator as an example, the first-order difference 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 over 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 lag order ranges from 1 to 10, generating sequences with lags of 1 time unit, 2 time units, up to 10 time units. The motor stator current data is shifted backward by 1 to 10 time units to obtain 10 lag sequences, which are then combined with the original sequence to form a lag sequence group. A lag of 1 time unit means shifting the entire motor stator current data one sampling point backward, i.e., the current time point uses the value from the previous time point.

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

[0113] The causal significance statistic is obtained by solving for the regression coefficients and lag orders. The causal significance statistic is calculated using the variance ratio statistic, which is the sum of squares of regression coefficients for all lag orders divided by the sum of squares of the residuals, then multiplied by the ratio of degrees of freedom. In the calculation of the variance ratio statistic, the numerator (degrees of freedom) is the lag order, and the denominator (degrees of freedom) is the total number of samples minus the number of regression coefficients minus 1. In the causal test of the effect of motor stator current on the temperature change of the main bearing, the calculated variance ratio statistic is 8.76.

[0114] The Granger causality test results are obtained by comparing the causal significance statistic with a preset critical value. The critical value is determined based on the quantiles of the variance ratio distribution, with a significance level set to 0.01. In this embodiment, the critical value for the variance ratio distribution is 4.25. Since the calculated variance ratio statistic of 8.76 is greater than the critical value of 4.25, it is determined that there is a significant Granger causal relationship between the motor stator current and the main bearing temperature change, i.e., the change in motor stator current can predict the change in 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 that passed the significance test in the Granger causality test results. The Granger causality coefficient is calculated by weighting the regression coefficients, with the weights being the reciprocals of the lag orders. Only regression coefficients that passed the significance test are included in the calculation. The criterion for significance testing is that the t-statistic obtained by dividing the regression coefficient by its standard error is greater than the critical value of the t-distribution for the corresponding degrees of freedom. The Granger causality coefficient of the motor stator current on the main bearing temperature is calculated to be 0.37. Similarly, the Granger causality coefficients between other component pairs are calculated.

[0116] Treating each component in the topological connection as a node, directed edges are established between components with a Granger causality coefficient greater than a preset threshold. The preset threshold 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 dependent variable component, representing the direction of the causal relationship. Based on the calculation results, 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 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 directed edges, and directed edges forming loops are deleted to obtain a causal directed acyclic graph. Loop detection uses a depth-first search algorithm, starting from each node and traversing along the directed edges. If the traversal returns to the starting node, a loop is identified. After a loop is detected, the directed edge with the smallest Granger causality coefficient in the loop is deleted. Loop detection is repeated until no loops exist. In this embodiment, a loop is detected between the motor rotor and the auxiliary bearing. The Granger causality coefficient from the auxiliary bearing to the motor rotor is 0.22, which is less than the 0.31 from the motor rotor to the auxiliary bearing. Therefore, the directed edge from the auxiliary bearing to the motor rotor is deleted, ultimately resulting in a loop-free causal directed graph.

[0118] In this embodiment, by introducing the real physical connection relationship between components as a priori constraint, the potential action path is limited in the causal modeling process. This effectively avoids non-physical associations and spurious causal relationships that occur in traditional pure data-driven methods, and improves the ability to specifically describe fault triggering and propagation behavior. By performing stationarity tests on component monitoring data and differentiating non-stationary sequences, the risk of spurious regression and misjudgment of causal relationships is significantly reduced, and the stability and credibility of Granger causality test results are improved. By introducing regression analysis with multiple time lag orders and screening causal relationships based on significance statistics, the fault impact path and key transmission nodes can be more precisely characterized, which helps to identify the core components in fault propagation.

[0119] In one alternative implementation,

[0120] 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:

[0121] 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.

[0122] 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.

[0123] 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.

[0124] Based on component type, fused features are assigned to a causal directed acyclic graph and multi-dimensional feature decoupling is performed to obtain state degradation feature vectors and dynamic evolution feature vectors. During feature assignment, different mapping rules are used for different component types to map fused features to corresponding nodes. Fusion features for motor components include dimensions such as current, voltage, and temperature, and all dimensions are preserved during mapping; fused features for bearing components include dimensions such as temperature, vibration, and noise, and all dimensions are preserved during mapping; fused features for pump components include dimensions such as pressure, flow rate, and vibration, and all dimensions are preserved during mapping. Feature decoupling uses independent component analysis (ICA) to decompose the fused features into mutually independent components. The number of decoupling components is set to 2, and a fast ICA algorithm is used to solve the problem. The number of iterations is set to 1000, and the convergence threshold is set to 0.0001, yielding state degradation feature vectors and dynamic evolution feature vectors. For the main bearing component, the fused feature is a 16-dimensional vector, which, after feature decoupling, yields an 8-dimensional state degradation feature vector and an 8-dimensional dynamic evolution feature vector. The state degradation feature vector mainly characterizes the steady-state fault characteristics of the component and reflects the inherent performance degradation degree of the component; the dynamic evolution feature vector characterizes the dynamic fault characteristics of the component and reflects the changing trend of the component's operating state.

[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, taking the state degradation feature vector as the observed data, and establishing a Gaussian mixture model as the prior distribution. The mixture component is set to 3, corresponding to the normal state, the minor fault state, and the severe fault state, respectively. The variational inference algorithm iteratively solves for the posterior distribution parameters, with the number of iterations set to 500 and the convergence threshold set to 0.0001. The posterior distribution parameters include the mean vector and covariance matrix. For the main bearing assembly, the mean vector for the normal state is [0.12, 0.08, 0.15, 0.10, 0.07, 0.09, 0.11, 0.13], the mean vector for the minor fault state is [0.45, 0.38, 0.42, 0.40, 0.35, 0.39, 0.43, 0.41], and the mean vector for the severe fault state is [0.82, 0.78, 0.85, 0.80, 0.75, 0.79, 0.83, 0.81]. By sampling 1000 times from the posterior distribution and statistically analyzing the distribution proportions of the samples in each state, the probability value of the fault state is calculated. For the main bearing assembly, the sampling results show that the probability of normal state is 0.25, the probability of minor fault state is 0.55, and the probability of serious fault state is 0.20. The initial fault probability value is the sum of the probabilities of minor fault and serious fault, which is 0.75.

[0126] Starting from the root node of the causal directed acyclic graph, obtain the initial fault probability value of the current node, the fault probability value of the corresponding parent node, and the dynamic evolution feature vector. The root node is a node with an in-degree of 0, i.e., a node with no directed edges pointing to it. In the causal directed acyclic graph of downhole equipment, the power control unit is usually the root node. Taking the motor stator node as an example, obtain the corresponding initial fault probability value of 0.35, the fault probability value of the parent node power control unit of 0.28, and the dynamic evolution feature vector passed by the power control unit as [0.23, 0.25, 0.20, 0.22, 0.24, 0.21, 0.26, 0.19].

[0127] A time-varying conditional probability transfer kernel is constructed based on the Granger causality coefficient of directed edges and the dynamic evolution feature vector. This kernel is then convolved with the fault probability value using a tensor to obtain the causal propagation probability distribution. The time-varying conditional probability transfer kernel represents the degree of influence of the parent node's fault state on the child node's fault state. Its construction is based on the Granger causality coefficient and the dynamic evolution feature vector. The transfer kernel construction method uses the Granger causality coefficient as the basic weight, and then adjusts it according to the values ​​of each dimension of the dynamic evolution feature vector. The adjustment formula is the basic weight multiplied by the weighted sum of the dynamic evolution features, where 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. Therefore, the time-varying conditional probability transfer kernel is 0.525. The parent node's fault probability value is then convolved with the time-varying conditional probability transfer kernel using a tensor, which multiplies the parent node's fault probability value by the time-varying conditional probability transfer kernel, to obtain the causal propagation probability distribution. The causal propagation probability distribution of the power control unit to the motor stator is 0.28 multiplied by 0.525, which equals 0.147.

[0128] The aggregated prior probability is obtained by performing a marginalization integral on the causal transitivity probability distribution. Marginalization integration refers to the process of integrating or summing the conditional probabilities 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 transitivity probability distributions of each parent node and then divide by the number of parent nodes to obtain the aggregated prior probability. For the motor stator node, which has only one parent node (the power control unit), the aggregated prior probability is equal to the corresponding causal transitivity probability distribution of 0.147.

[0129] The updated fault probability value is obtained by optimizing the solution based on the aggregated prior probability and the initial fault probability value, and then passed to the child nodes along with the dynamic evolution feature vector. The optimization solution employs a Bayesian update method, using the aggregated prior probability as the prior and the initial fault probability value as the likelihood, and calculating the posterior probability as the updated fault probability value. The calculation formula is: prior probability multiplied by likelihood probability divided by normalization constant. The normalization constant ensures that the sum of probabilities is 1 through integration or summation. For the motor stator node, the aggregated prior probability is 0.147, the initial fault probability value is 0.35, and the calculated updated fault probability value is 0.41. The updated fault probability value, along 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 node, is passed to the motor rotor, a child node of the motor stator. This process is repeated, traversing all nodes in the causal directed acyclic graph, to complete one iteration of the update.

[0130] The Kourbak-Leibler divergence of the fault probability values ​​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 terminates. The Kourbak-Leibler divergence measures the degree of difference between two probability distributions. It is calculated by multiplying the logarithm of the probability value after iteration by the probability value after iteration, then subtracting the logarithm of the probability value before iteration multiplied by the probability value after iteration, and summing the results over all states to obtain the divergence value. The preset convergence threshold is 0.01. Before the first iteration, the fault probability value of the motor stator node is 0.35, and after the first iteration it is 0.41. The calculated Kourbak-Leibler divergence is 0.024, which is greater than the convergence threshold of 0.01, so the iteration continues. After the second iteration, the fault probability value of the motor stator node is 0.43, and the calculated Kourbak-Leibler divergence is 0.008, which is less than the convergence threshold of 0.01. At this point, the motor stator node is considered to have converged. The overall iteration process terminates when the Kourbach-Leibler divergence of all nodes is less than the convergence threshold.

[0131] A comprehensive fault score is calculated by extracting the fault probability value and the dynamic evolution feature vector of each node. The comprehensive fault score is calculated by weighting the fault probability value and the norm of the dynamic evolution feature vector, with weighting coefficients of 0.7 and 0.3 respectively. The norm of the dynamic evolution feature vector is calculated as the square root of the sum of squares of each dimension. For the main bearing node, the final converged fault probability value is 0.82, the dynamic evolution feature vector norm is 0.65, and the calculated comprehensive fault score is 0.77. For the motor stator node, the final converged fault probability value is 0.43, the dynamic evolution feature vector norm is 0.42, and the calculated comprehensive fault score is 0.43. After calculating the comprehensive fault score for all nodes, the component corresponding to the node with the highest comprehensive fault score is identified as the fault source. In this embodiment, the main bearing node has the highest comprehensive fault score of 0.77, therefore, the main bearing is determined to be the fault source.

[0132] The diagnostic reasoning results are obtained by integrating the fault information of the fault source component and its related components. The comprehensive fault score of the fault source component is used as the main indicator, and the comprehensive fault scores of the child nodes of the fault source component are combined as related fault indicators to form a complete diagnostic report. The diagnostic report includes the name of the fault source component, the fault probability, the fault score, and the possible fault type. The fault type is determined by pattern recognition of the state degradation feature vector. By comparing it with a pre-labeled fault pattern library, the fault type with the highest matching degree is found. For the main bearing assembly, its state degradation feature vector has a matching degree of 92% with the bearing inner ring wear fault mode, so the fault type is judged to be inner ring wear. The final diagnostic reasoning result is: the main bearing assembly has an inner ring wear fault, with a fault probability of 0.82 and a comprehensive fault score of 0.77; the associated affected components include the sealing device, with a fault probability of 0.45 and a comprehensive fault score of 0.38, which may lead to poor sealing; it is recommended to replace the main bearing and check the sealing device.

[0133] In this embodiment, by mapping fused features to corresponding nodes in a causal directed acyclic graph according to component type, and decoupling the state degradation and dynamic evolution of features, the fault representation no longer mixes information with different time scales and physical meanings. This helps to distinguish the different effects of long-term degradation and short-term disturbances on the fault probability, improving the precision and reliability of fault state characterization. By propagating probability along the causal direction in the causal directed acyclic graph and introducing a time-varying conditional probability transfer mechanism jointly determined by Granger causality coefficient and dynamic evolution features, the transmission process of fault impact can adaptively adjust with changes in operating state, improving the physical consistency and temporal rationality of fault propagation modeling. By aggregating the causal transmission probability distribution and jointly optimizing it with the initial fault probability of the node itself, misjudgment and probability oscillation are significantly reduced.

[0134] Figure 2 This is a flowchart illustrating the intelligent fault source diagnosis process of the real-time fault diagnosis method for downhole equipment based on edge computing, as described in an embodiment of the present invention.

[0135] In one alternative implementation,

[0136] 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:

[0137] 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.

[0138] 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.

[0139] 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.

[0140] The failure probability values ​​of each component are extracted from the diagnostic inference results and normalized to obtain the failure probability distribution. The diagnostic inference results contain failure probability values ​​for each component. After extracting the failure probability values, normalization is required to ensure data comparability. The normalization process uses a summation normalization method, which divides the failure probability value of each component by the sum of the failure probability values ​​of all components, ensuring that the sum of the normalized failure probability values ​​is 1. Taking a downhole electric pump system as an example, the diagnostic inference results contain nine components: power control unit, motor stator, motor rotor, main bearing, auxiliary bearing, sealing device, pump inlet, pump impeller, and pump outlet. Their original failure probability values ​​are 0.28, 0.43, 0.39, 0.82, 0.35, 0.45, 0.30, 0.25, and 0.20, respectively, with a total failure probability value of 3.47. The fault probability distributions obtained after normalization are 0.081, 0.124, 0.112, 0.236, 0.101, 0.130, 0.086, 0.072, and 0.058, respectively.

[0141] The information entropy of the fault probability distribution is calculated to obtain the current task entropy value. Information entropy measures the degree of uncertainty of the system. It is calculated by multiplying the normalized fault probability value of each component by the negative of its logarithm, and then summing all the results. The base of the logarithm is set to the base of the natural logarithm, e. For the fault probability distribution of the downhole electric pump system, the calculated information entropy value is 2.14. The information entropy value reflects the complexity and uncertainty of the current fault diagnosis task. A higher current task entropy value indicates greater uncertainty in fault diagnosis, requiring a more complex diagnostic strategy; a lower current task entropy value indicates a more certain fault location, allowing for a more targeted diagnostic strategy.

[0142] The causal propagation probability distribution between components is extracted from the diagnostic inference results, and the conditional entropy is calculated to obtain the causal association entropy value. The diagnostic inference results contain the causal propagation probability distribution between components, which reflects the possible paths and intensity of fault propagation between components. The conditional entropy calculation method is as follows: for each component pair, the uncertainty of the target component's fault probability distribution is calculated under the condition of known source component fault probability. The source component's fault probability value is multiplied by the corresponding conditional probability to obtain the joint probability. The negative of the logarithm of the joint probability is calculated and then multiplied by the joint probability. The summation is then performed over all possible component pairs. In the downhole electric pump system, the causal propagation probability from the power control unit to the motor stator is 0.525, the causal propagation probability from the motor stator to the motor rotor is 0.480, the causal propagation probability from the motor rotor to the main bearing is 0.420, and so on, forming a complete causal propagation probability matrix. The conditional entropy calculated based on these causal propagation probabilities is 1.75. This value reflects the complexity and uncertainty of the fault propagation relationship between components.

[0143] The comprehensive entropy feature vector is obtained by weighted summing of the current task entropy value and the causal correlation entropy value. The weighting coefficients are set empirically, with the current task entropy value weighted at 0.6 and the causal correlation entropy value weighted at 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 comprehensive entropy feature vector obtained by weighted summation is 2.14 × 0.6 + 1.75 × 0.4 = 1.984. The comprehensive entropy feature vector comprehensively considers the uncertainty of fault distribution and the complexity of fault propagation relationships, providing a basis for subsequent diagnostic strategy selection.

[0144] The historical task entropy value and historical causal association entropy value corresponding to each diagnostic strategy in the preset diagnostic strategy library are obtained and weighted according to preset weights to obtain the historical comprehensive entropy feature vector. The diagnostic strategy library contains a variety of diagnostic strategies for different fault scenarios, and each strategy has corresponding applicable scenario characteristics. The historical task entropy value and historical causal association entropy value are obtained based on the statistics of historical diagnostic cases, reflecting the characteristics of the applicable scenarios of various diagnostic strategies. The diagnostic strategy library for the downhole electric pump system includes 5 strategies: single component deep diagnosis strategy, multi-component parallel diagnosis strategy, causal chain tracing diagnosis strategy, full system scan diagnosis strategy, and regional stepwise diagnosis strategy. The historical task entropy value of the single component deep diagnosis strategy is 1.2, and the historical causal association entropy value is 0.8; the historical task entropy value of the multi-component parallel diagnosis strategy is 1.8, and the historical causal association entropy value is 1.5; the historical task entropy value of the causal chain tracing diagnosis strategy is 2.2, and the historical causal association entropy value is 1.9; the historical task entropy value of the full system scan diagnosis strategy is 2.5, and the historical causal association entropy value is 2.0; the historical task entropy value of the regional stepwise diagnosis strategy is 1.9, and the historical causal association entropy value is 1.6. Using the same weighting coefficients as the current comprehensive entropy feature vector, the historical entropy values ​​of each strategy are weighted and summed to obtain historical comprehensive entropy feature vectors of 1.04, 1.68, 2.08, 2.30 and 1.78 for each strategy.

[0145] The Euclidean distance between the current comprehensive entropy feature vector and historical comprehensive entropy feature vectors is calculated, and a reciprocal transformation is performed to obtain the distance similarity. The Euclidean distance is calculated by subtracting the current comprehensive entropy feature vector from each historical comprehensive entropy feature vector and taking the absolute value. The reciprocal transformation uses 1 divided by the Euclidean distance plus 1 to ensure that the distance similarity value is between 0 and 1. The current comprehensive entropy feature vector of the downhole electric pump system is 1.984, and its Euclidean distances with each historical comprehensive entropy feature vector are 0.944, 0.304, 0.096, 0.316, and 0.204, respectively. After the reciprocal transformation, the distance similarities obtained are 0.514, 0.767, 0.912, 0.760, and 0.830, respectively. A higher distance similarity indicates that the features of the current diagnostic task are closer to the applicable scenario of the corresponding diagnostic strategy.

[0146] The entropy pattern similarity is obtained by calculating the cosine similarity between the current task entropy value and the historical task entropy values. The cosine similarity is calculated by multiplying the current task entropy value by the historical task entropy value and dividing by the product of their moduli. Since entropy is a scalar, the cosine similarity simplifies to the smaller of the ratio of the two values ​​and 1. The current task entropy value of the downhole electric pump system is 2.14, and its cosine similarities with the entropy values ​​of each historical task are 0.561, 0.841, 0.973, 0.856, and 0.888, respectively. Entropy pattern similarity reflects the degree of similarity between the current diagnostic task and historical diagnostic cases in terms of complexity and uncertainty.

[0147] A weighted fusion of distance similarity and entropy pattern similarity was performed to obtain a comprehensive similarity score. The weighting coefficients were empirically set, with a weighting coefficient of 0.7 for distance similarity and 0.3 for entropy pattern similarity. The comprehensive similarities obtained after weighted fusion of distance similarity and entropy pattern similarity for each diagnostic strategy of the downhole electric pump system were 0.528, 0.789, 0.930, 0.789, and 0.847, respectively. This comprehensive similarity score fully considers the matching degree between the current diagnostic task and the applicable scenarios of each historical diagnostic strategy, providing a reliable basis for the selection of diagnostic strategies.

[0148] The diagnostic strategies were sorted in descending order based on comprehensive similarity, with the top-ranked strategy being the current execution strategy. The ranking of the diagnostic strategies from highest to lowest comprehensive similarity was as follows: causal chain tracing diagnostic strategy (0.930), regional stepwise diagnostic strategy (0.847), multi-component parallel diagnostic strategy (0.789), full system scan diagnostic strategy (0.789), and single-component deep diagnostic strategy (0.528). Therefore, the causal chain tracing diagnostic strategy was selected as the current execution strategy. The causal chain tracing diagnostic strategy is a method of step-by-step diagnosis along the fault propagation path, suitable for situations where the fault has a clear propagation link and there is a strong causal relationship between components. For the downhole electric pump system, the diagnosis starts from the main bearing, tracing upwards along the causal chain to the motor rotor, motor stator, and finally the power control unit, while simultaneously tracing downwards to components such as the sealing device and pump inlet, gradually confirming the fault state and impact of each component.

[0149] In this embodiment, by modeling the failure probability distribution of each component using information entropy, the degree of uncertainty contained in the diagnostic results can be quantified, which can more comprehensively reflect the complexity and information completeness of the current diagnostic task. By weightedly fusing task entropy and causal association entropy to form a comprehensive entropy feature vector, a unified representation of the overall state of the diagnostic task is achieved. By measuring the similarity between the comprehensive entropy feature of the current diagnostic task and the comprehensive entropy feature corresponding to historical strategies, adaptive strategy matching based on historical experience is achieved, improving the accuracy and stability of strategy recommendation. By ranking diagnostic strategies through comprehensive similarity and selecting the optimal strategy, a diagnostic strategy with controllable risk and higher information benefit can be automatically matched under different failure complexity and causal clarity conditions. This significantly improves the adaptability and decision rationality of the diagnostic process under complex working conditions, helps to reduce the risk of misdiagnosis and improve the overall diagnostic efficiency.

[0150] In one alternative implementation,

[0151] 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:

[0152] 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.

[0153] 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.

[0154] 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.

[0155] 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 acquired and subjected to adaptive wavelet decomposition to obtain frequency band components. The diagnostic inference results contain the fault probability values ​​of each component, and the component with the highest fault probability value is selected as the target diagnostic component by comparison. Taking the downhole electric pump system as an example, the diagnostic inference results show 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. Real-time monitoring time-series data of the main bearing is acquired, including vibration signals, temperature signals, etc., with a sampling frequency of 1000Hz for vibration signals and 1Hz for temperature signals, and a data duration of 30 minutes. Adaptive wavelet decomposition is performed on the acquired real-time monitoring time-series data, using an improved empirical mode decomposition method to decompose the signal into multiple frequency band components. The parameter settings for adaptive wavelet decomposition include a decomposition level of 5 levels, a db4 wavelet as the wavelet basis function, a soft threshold as the threshold function, and a maximum-minimum threshold criterion for threshold selection. After performing five-level wavelet decomposition on the vibration signal of the main bearing, five detail components and one approximate component are obtained, corresponding to different frequency bands: detail component d1 corresponds to the frequency band 500-1000Hz, detail component d2 corresponds to the frequency band 250-500Hz, detail component d3 corresponds to the frequency band 125-250Hz, detail component d4 corresponds to the frequency band 62.5-125Hz, detail component d5 corresponds to the frequency band 31.25-62.5Hz, and approximate component a5 corresponds to the frequency band 0-31.25Hz.

[0156] The instantaneous frequency and instantaneous amplitude of the frequency band components are calculated to construct a time-frequency joint characterization matrix. The instantaneous frequency is calculated using the Hilbert transform method; a Hilbert transform is performed on each frequency band component to obtain an analytic signal, and the phase derivative of the analytic signal is calculated to obtain the instantaneous frequency. The instantaneous amplitude is calculated as the magnitude of the analytic signal. For the six frequency band components obtained from the decomposition of the main bearing vibration signal, their instantaneous frequency and instantaneous amplitude are calculated respectively, resulting in a time-frequency joint characterization matrix with dimensions of 6×2×18000, where 6 represents the number of frequency band components, 2 represents the characteristic dimension (frequency and amplitude), and 18000 represents the number of time points (30 minutes × 60 seconds × 10 sampling points / second). In the d3 frequency band component (125-250Hz), the instantaneous frequency is stable at around 185Hz under normal conditions, and the instantaneous amplitude fluctuates within the range of 0.5-0.8mm / s; however, at abnormal times, the instantaneous frequency shows a significant drift, reaching a maximum of 220Hz, and the instantaneous amplitude suddenly increases to 2.5mm / s.

[0157] The gradient rate of change corresponding to the joint time-frequency characterization matrix is ​​calculated to obtain the time-frequency evolution trajectory, and amplitude abrupt change points and frequency drift points are extracted as anomalous time-frequency feature points. The gradient rate of change is calculated by dividing the difference between the instantaneous frequency and instantaneous amplitude at adjacent time points by the time interval. For each frequency band component, the gradient rate of change of its instantaneous frequency and instantaneous amplitude is calculated to form the time-frequency evolution trajectory. The extraction of anomalous time-frequency feature points is based on an anomaly detection algorithm, with a threshold set at three times the standard deviation of the gradient rate of change under normal conditions. For the d3 frequency band component, the standard deviation of the gradient rate of change of instantaneous frequency under normal conditions is 0.5 Hz / s, and the standard deviation of the gradient rate of change of instantaneous amplitude is 0.1 mm / s. 2 Therefore, the threshold for determining the frequency drift point is 1.5 Hz / s, and the threshold for determining the amplitude abrupt change point is 0.3 mm / s. 2 Based on these thresholds, 15 frequency drift points and 12 amplitude abrupt change points were detected in the d3 frequency band component and marked as anomalous time-frequency feature points.

[0158] Density distribution calculation and peak detection were performed on abnormal time-frequency feature points to determine the abnormal clustering times. Density distribution calculation employed a kernel density estimation method, projecting each abnormal time-frequency feature point onto the time axis and weighted summing to form a density curve. A Gaussian kernel was selected as the kernel function, with a bandwidth parameter set to 30 seconds. Peak detection used a local maximum detection algorithm, setting the minimum peak height to twice the average density and the minimum peak distance to 60 seconds. After calculating the density distribution of the abnormal time-frequency feature points of the main bearing, three density peaks were detected, corresponding to the 8th, 15th, and 22nd minutes after the start of recording. These time points were identified as the abnormal clustering times.

[0159] Data segments for abnormal time periods are extracted based on the anomaly clustering moments. The extraction method involves extending 30 seconds forward and backward from the anomaly clustering moment to form a 1-minute data segment. For the three detected anomaly clustering moments, corresponding data segments are extracted to obtain three abnormal time period data segments. The first abnormal time period data segment corresponds to the time from 7 minutes 30 seconds to 8 minutes 30 seconds, the second corresponds to the time from 14 minutes 30 seconds to 15 minutes 30 seconds, and the third corresponds to the time from 21 minutes 30 seconds to 22 minutes 30 seconds. These abnormal time period data segments contain the most significant characteristic information of the fault, providing crucial data for subsequent causal verification.

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

[0161] Cross-correlation analysis was performed on historical monitoring time-series data and data segments from abnormal periods to obtain time-delay correlation coefficients. The cross-correlation analysis method calculates the correlation coefficient between two signals at different time delays, identifying the maximum correlation coefficient and its corresponding time delay. For each pair of adjacent components, the cross-correlation function between its historical data and data segments from abnormal periods was calculated, with a time delay range of -5 seconds to 5 seconds and a step size of 0.1 seconds. For the motor rotor and main bearing, the cross-correlation function between its current and vibration signal was calculated, yielding a maximum correlation coefficient of 0.65 and a corresponding time delay of 0.3 seconds; for the main bearing and auxiliary bearing, the cross-correlation function between their vibration signals was calculated, yielding a maximum correlation coefficient of 0.72 and a corresponding time delay of 0.2 seconds; for the main bearing and sealing device, the cross-correlation function between their vibration and pressure signals was calculated, yielding a maximum correlation coefficient of 0.48 and a corresponding time delay of 0.5 seconds. These maximum correlation coefficients are recorded as time-delay correlation coefficients.

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

[0163] The corrected causal verification link is determined based on the causal deviation and a preset verification threshold. The verification threshold is set to 0.25, meaning that when the causal deviation is greater than 0.25, the causal relationship is considered to have changed significantly and needs correction; when the causal deviation is less than 0.25, the causal relationship is considered relatively stable and the original link is retained. For the link from the motor rotor to the main bearing, the causal deviation is 0.23, which is less than the threshold of 0.25, so the link is retained. For the link from the main bearing to the auxiliary bearing, the causal deviation is 0.37, which is greater than the threshold of 0.25, so the link is marked as needing correction. For the link from the main bearing to the sealing device, the causal deviation is 0.09, which is less than the threshold of 0.25, so the link is retained. The corrected causal verification link is: motor rotor → main bearing → sealing device, removing the link from the main bearing to the auxiliary bearing.

[0164] The causal verification confidence level 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 initial causal verification link contains 4 components: motor rotor, main bearing, auxiliary bearing, and sealing device; the corrected causal verification link contains 3 components: motor rotor, main bearing, and sealing device. The ratio of the number of remaining components to the number of initial components is 3 / 4 = 0.75, meaning the causal verification confidence level is 0.75.

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

[0166] The optimized diagnostic results were obtained by combining the corrected fault probability value and the corrected causal verification link. The optimized diagnostic results include information such as the faulty component, the corrected fault probability value, the fault type, and the fault propagation path. For the downhole electric pump system, the optimized diagnostic results are: wear fault exists in the inner ring of the main bearing, with a corrected fault probability value of 0.72; the fault propagation path is motor rotor → main bearing → sealing device, indicating that the main bearing fault is affected by the operating state of the motor rotor, and also affects the sealing device; recommended measures include replacing the main bearing, checking the operating 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 characterization on the component with the highest failure probability, the anomaly analysis is upgraded from the original low-dimensional statistical features to a high-resolution time-frequency evolution level, which improves the sensitivity of identifying sudden and hidden faults. The abnormal time period is determined based on the density clustering of abnormal feature points, which effectively avoids false triggering caused by noise fluctuations, making the anomaly location more concentrated and accurate on the time axis. By comparing the time delay correlation coefficient with the Granger causality coefficient and correcting the causal verification link accordingly, the causal structure can be dynamically verified and corrected according to real-time anomaly features, which significantly improves the adaptability of the diagnostic system under conditions of changing operating conditions and structural uncertainty. By calculating the causal verification confidence based on the scale of the corrected causal link and using this confidence to make a second correction to the failure probability, the final failure probability not only reflects the model inference result, but also integrates the anomaly time-frequency features and causal consistency verification information, which significantly improves the reliability and interpretability of the diagnostic results.

[0168] A second aspect of this invention provides a real-time fault diagnosis system for downhole equipment based on edge computing, comprising:

[0169] 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.

[0170] 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.

[0171] 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.

[0172] 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.

[0173] A third aspect of the present invention provides an electronic device, comprising:

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

[0175] A fourth aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0176] This invention can be a method, apparatus, system, and / or computer program product. The computer program product may include a computer-readable storage medium having computer-readable program instructions loaded thereon for performing various aspects of the invention.

[0177] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features therein. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

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