A pipeline risk dynamic perception and visualization method based on multi-dimensional space-time graph convolution

Abstract

The application discloses a kind of based on multi-dimensional space-time graph convolution pipe risk dynamic perception and visualization method, comprising: constructing high-precision pipe internet of things perception matrix, collecting multi-physical field operation data, and constructing high-fidelity space-time tensor;From time-frequency domain and topological space domain dual dimension extraction data's deep feature;Construct double-flow space-time graph convolution attention network, output rich semantic high-dimensional risk characteristics;Introduce cloud model theory to construct adaptive risk entropy field, convert risk value into cloud drop distribution, realize the unified quantification of risk certainty and randomness;Based on force-oriented layout algorithm and particle system design dynamic three-dimensional visualization engine, risk data is mapped into holographic field effect in digital twin interface in real time;Finally, a hierarchical closed-loop response mechanism is established.The application not only significantly improves the identification rate of micro-leakage and early fault, but also greatly improves the operation and maintenance decision efficiency and situational awareness ability through immersive visualization means.

Description

Technical Field

[0001] This invention relates to the fields of industrial Internet of Things and digital twin technology, specifically to a pipeline risk dynamic perception and visualization method based on multidimensional spatiotemporal graph convolution, which is used to perform real-time perception of all elements of risk in industrial pipeline networks, dynamic evolution analysis and immersive holographic display, so as to improve the safety and management efficiency of pipeline network operation. Background Technology

[0002] Industrial pipeline networks (such as oil and gas, water supply, and heating networks) are the "lifeline" of modern industry and urban operations, and their safe and stable operation is crucial. With the expansion of pipeline network scale and the increase in service life, pipelines face various risks such as corrosion, leakage, and overpressure. However, traditional pipeline network monitoring and risk management methods mainly rely on manual inspections or simple threshold alarms based on SCADA systems, which have many limitations. First, traditional methods are mostly based on single-point monitoring, lacking consideration of the overall fluid dynamics and topology of the pipeline network, resulting in delayed perception and a high false alarm rate. Second, existing data analysis methods struggle to effectively integrate multi-source heterogeneous data (such as pressure, vibration, and acoustic emission), failing to deeply explore the spatiotemporal dependencies within the data and making it difficult to detect early risks such as minor leaks. Finally, existing visualization methods are mostly limited to two-dimensional charts, lacking intuitive, immersive risk display and interactive capabilities, making it difficult for managers to accurately assess the risk situation in a timely manner. With the rapid development of the Internet of Things, artificial intelligence, and digital twin technologies, there is an urgent need for a pipeline network risk perception method that can deeply integrate spatiotemporal characteristics, possess high-precision risk reasoning capabilities, and further achieve holographic visualization. Summary of the Invention

[0003] The purpose of this invention is to provide a pipeline risk dynamic perception and visualization method based on multidimensional spatiotemporal graph convolution. By using a global pipeline perception network, a deep spatiotemporal graph convolutional network, and a cloud model evaluation system, the method performs multidimensional fusion analysis and dynamic risk perception on pipeline operation data, and displays the risk situation holographically through a three-dimensional visualization engine.

[0004] The objective of this invention is achieved through the following technical solution: This method and system for dynamic perception and visualization of pipeline risks based on multidimensional spatiotemporal graph convolution includes the following steps: S1: Construct a global pipeline sensing network, collect multi-physics field operation data of pipelines and perform spatiotemporal alignment, and construct a multi-source heterogeneous spatiotemporal data tensor; S2: Perform adaptive variational mode decomposition and graph signal processing on the spatiotemporal data tensor to construct a weighted graph Laplace matrix that reflects the pipeline topology; S3: Construct a dual-stream spatiotemporal graph convolutional attention network to extract the temporal evolution features and spatial diffusion features of risk in parallel, and fuse them to output a risk latent feature map; S4: Construct a risk assessment system based on cloud models, calculate the expected value, entropy, and hyperentropy of risks, and quantify the dynamic risk probability and confidence level; S5: Design a dynamic 3D visualization engine that uses a particle system to render the pipeline risk field and anomaly propagation streamlines in real time based on risk cloud droplet parameters; S6: Establish a hierarchical closed-loop response mechanism to trigger edge control strategies based on risk levels.

[0005] Furthermore, step S1 specifically includes: S101: Based on the physical space and fluid dynamics characteristics, the pipeline network topology graph G=(V, E, A) is constructed by mapping physical entities to mathematical space. S102: Collect data from each node v i Construct the original measurement vector from the multidimensional attribute data at time step t. ; S103: Introduce a weighted covariance alignment algorithm to perform time synchronization correction on the original measurement vectors that are out of sync due to sensor transmission delay; S104: Perform time window slicing on multi-source data to construct a network-level spatiotemporal input tensor X. raw ; Furthermore, step S101 specifically includes: First, determine the set V of sensor nodes, and map the N sensors distributed in the pipeline network as graph nodes, denoted as V={v1, v2, ..., v N}, each node v i It represents the location of a specific physical monitoring point; Then, a set of pipe physical connections and fluid dynamics associated edges E is constructed to capture direct physical connections and indirect hydraulic coupling relationships. Finally, the weighted adjacency matrix W, which reflects the medium transmission resistance or distance attenuation factor between nodes, is calculated.

[0006] Furthermore, step S102 specifically includes: By deploying pressure, flow, temperature, vibration, and acoustic emission sensors, multi-source data are collected simultaneously at the same time. For scalar data such as pressure and flow, Z-Score normalization is used to eliminate dimensional differences, mapping different physical quantities to the same distribution scale. The normalization formula is as follows: Where x is the original physical quantity currently collected, and x' is the standardized feature value. This represents the average value of the sensor over a historical sliding window. The standard deviation is used to eliminate the impact of absolute numerical differences on the convergence speed of subsequent gradient descent. For high-frequency vibration signals, the Fast Fourier Transform (FFT) is used to extract frequency domain features. FFT is an efficient algorithm for Discrete Fourier Transform, and its transformation formula is as follows: Where x[n] is the time-domain discrete vibration signal sequence, N is the number of sampling points, X[k] is the frequency-domain complex sequence, k is the frequency index, and j is the imaginary unit; Based on this, the top k dominant frequency components with the largest amplitude spectrum |X[k]| are selected as features; Finally, the processed data is encapsulated into a high-dimensional vector with feature dimension F, and the original measurement vector is represented as: in, , , These are instantaneous pressure, flow rate, and temperature observations, respectively. These represent the historical mean and standard deviation of the corresponding physical quantity, respectively. k (|X vib [k]|) represents the amplitude values ​​of the first k dominant frequencies of the vibration signal. The logarithmic value representing the energy of the acoustic emission signal is used to amplify the characteristics of minute leaks. Furthermore, step S103 specifically includes: First, calculate the spatiotemporal correlation coefficient between node i and each node j in its neighboring node set N(i), and use the Gaussian kernel function to measure the similarity of the data feature distribution; Then, the normalized weighted alignment factor w is calculated based on the correlation coefficient. ij To ensure that the sum of the contribution weights of neighboring nodes to the current node is 1, the weighted alignment factor is calculated as follows: Among them, the denominator term This is the sum of the correlation coefficients of all neighboring nodes of node i, used for normalization. Finally, a weighted alignment factor is used to aggregate the data of neighboring nodes to correct the observation value of the current node.

[0007] Furthermore, step S104 specifically includes: Set the historical observation time window length to T in Using a sliding window mechanism, all N nodes in the network are viewed in the past T... in Measurement vectors at each time step Stacked along the time dimension, they ultimately form a three-dimensional spacetime tensor: in, This is the state snapshot matrix of the entire network at time t; Furthermore, step S2 specifically includes: S201: Adaptive variational mode decomposition technology is used to automatically decompose the non-stationary time series signal f(t) collected by each node into K eigenmode functions with specific center frequencies in order to separate environmental noise and fault characteristics; S202: Based on the weighted adjacency matrix W constructed in step S101, construct the normalized graph Laplace matrix L, which reflects the topological structure and functional dependence of the pipeline system, from the perspective of graph theory. S203: Based on the graph Laplace matrix L, a graph filter is designed using Chebyshev polynomial approximation theory to project the pipeline spatiotemporal signal into the graph frequency domain. Furthermore, step S201 specifically includes: First, the input object and feature target are defined. f(t) is the raw instantaneous pressure collected in step S102. Or vibration signal within the historical time window T in A continuous signal sequence assembled sequentially within the time frame; specifically, the historical observation vector P at the current time t is defined. window : in, The number of sampling points within the time window. The sampling interval; To meet the requirements of AVMD for continuous domain signal processing, the discrete vector is reconstructed into a continuous-time function f(τ) using the Shannon interpolation formula. The concatenation and reconstruction formula is as follows: Wherein, sinc(x)=sin(πx) / (πx) is the interpolation kernel function, which transforms discrete sensor readings into a smooth signal containing transient pressure waveforms and nonlinear vibration modes; Then, the analytical form of the signal is constructed to obtain its spectral characteristics. The Hilbert transform is then used to convert the real signal into a complex analytic signal. The Hilbert transform formula is: in, Let pv denote the Hilbert transform operator, pv denote the Cauchy principal value integral, and u denote the principal value integral. k (τ) represents the k-th modal component; Based on this, the formula for constructing an analytical signal is: Among them, z k Let δ(t) be an analytic signal with a one-sided spectral characteristic (containing only positive frequencies), and let δ(t) be the Dirac distribution function. j is the imaginary unit; Furthermore, we construct a variational optimization problem, assuming the original signal consists of K frequencies revolving around the center frequency ω. k amplitude modulation and frequency modulation signal u k (t) are superimposed. In order to minimize the bandwidth of each mode, the spectrum of the analytic signal is modulated to the baseband by shifting the exponential term. The objective functional formula for AVMD optimization is: Among them, J AVMD Let ω represent the total bandwidth target to be minimized. k That is, the specific center frequency corresponding to this mode. This is a complex exponential twitch factor used to perform spectrum shifting operations, changing the center frequency ω... k Move to zero frequency, Here, δ(t) is the time gradient operator used to calculate the smoothness (i.e., bandwidth) of the demodulated baseband signal, and δ(t) is the Dirac distribution function. Constraints Ensure that the sum of the K decomposed modes can be used to reconstruct the original signal without loss; To solve this constrained variational problem, we introduce the Lagrange multiplier λ(t) and the quadratic penalty factor α, transforming it into an unconstrained optimization problem. The augmented Lagrange function formula is as follows: Where L is the augmented Lagrangian function; the first term of the formula The second term is a quadratic penalty term, where α is a penalty parameter balancing data fidelity and bandwidth constraints, used to ensure reconstruction accuracy while constraining modal bandwidth; The residual fidelity term represents the energy difference between the original signal and the sum of the decomposed modes, i.e., the reconstruction error; the third term in the formula... For Lagrange multipliers, This indicates the inner product operation, used to strictly enforce the satisfaction of equality constraints; Furthermore, the alternating direction multiplier method is used to alternately update the variables in the frequency domain until convergence is obtained to obtain K eigenmode functions u. k (t); Finally, a signal reconstruction operation is performed, and a subset S of modes containing valid fault information is selected based on the kurtosis criterion. valid Noisy modes are removed, and the clean signal is reconstructed to build a denoised feature tensor X. feat The formula for denoising feature reconstruction is: in, Let be the denoised feature vector of node i at time t. This represents the superposition of effective modes, indicating the denoised dynamic signal (such as a pressure wave). For static attribute features that do not require decomposition (such as temperature, flow rate baseline); Based on this, construct the denoising feature tensor X. feat This tensor is obtained by placing all nodes in the time window T. in Internally reconstructed feature vector Stacked along the time axis, the formula for constructing the denoising feature tensor is: in, Let X be the state snapshot matrix of the entire network after denoising at time t; the tensor X feat This will be used as the input data for the dual-stream network in step S3; Furthermore, step S202 specifically includes: First, calculate the degree matrix D, which is a diagonal matrix with diagonal elements D. ii The sum of all elements in the i-th row of the weighted adjacency matrix W reflects the connection density or importance of node i in the network. The calculation formula is: Among them, D ii W represents the diagonal element value in the i-th row and i-th column of the degree matrix. ij The element in the i-th row and j-th column of the weighted adjacency matrix represents the connection weight between node i and node j. Then, construct the degree matrix D, which is derived from D. ii The resulting diagonal matrix, i.e., the elements excluding the main diagonal, are D. ii Except for the elements at all other positions, all other elements are 0, satisfying D = diag(D 11 ,…,D NN ); Finally, the adjacency matrix is ​​symmetrically normalized using the degree matrix to construct the graph Laplacian matrix, which is calculated as follows: Among them, I N Let W be the N-order identity matrix, and W be the weighted adjacency matrix from step S101. Let L be the inverse square root of the degree matrix. The normalized Laplace matrix L has positive semi-definite properties, and its eigenvalues ​​λ i Distributed within the interval [0,2], this ensures the numerical stability and convergence of subsequent graph convolution operations; Furthermore, step S203 specifically includes: First, calculate the rescaled Laplacian matrix, mapping the eigenvalues ​​to the interval [-1, 1] to satisfy the domain requirement of the Chebyshev polynomial. The rescaling formula is: in, Let λ be the rescaled Laplace matrix. max The largest eigenvalue of the Laplace matrix L; Then, a Chebyshev polynomial sequence is defined using recursion to quickly compute feature extraction in , and the recursive formula is: in, Using the k-th order Chebyshev polynomial operator, the energy distribution characteristics of the pipeline network at different "spatial frequencies" can be extracted without explicitly calculating the eigenvectors, thereby quickly identifying local mutation anomalies (corresponding to high-frequency components) and global evolution trends (corresponding to low-frequency components). Furthermore, step S3 specifically includes: S301: Construct a two-stream feature extraction architecture, and convert the denoised feature tensor X output from step S2 into a single stream. feat Simultaneously inputting into two parallel processing branches: spatial graph convolutional stream and temporal attention stream, the aim is to capture the topological spatial dependency and temporal evolution of risk propagation in pipeline networks, respectively, and improve the feature extraction efficiency of the model by decoupling spatiotemporal features; S302: In spatial graph convolutional flow, design a multi-level dynamic graph convolutional layer, where the layer uses tensor X... feat Taking the spatial slices at each time step as input, convolution operation is performed by combining static physical topology with dynamic data association, and feature aggregation is performed using the Chebyshev polynomial generated in step S203. S303: In temporal attention flow, we design a causal gated temporal convolutional network and a multi-head temporal self-attention mechanism to capture long-short-term temporal dependencies. S304: Design a spatiotemporal fusion gating unit to adaptively fuse spatial flow characteristics H S With time flow characteristics H T ; Furthermore, step S302 specifically includes: First, calculate the spatial self-attention matrix A. att This matrix reflects the similarity weights between the features of each node at the current time. Then, construct the fusion graph Laplacian operator L. mix The physical topology and data-driven correlation are weighted and fused together. Finally, calculate the spatial stream output characteristics. ; Furthermore, step S303 specifically includes: First, local temporal features are extracted using dilated causal convolution. This step converts the tensor X... feat Unfolding along the node dimension, convolving the time series of each node, the output formula for dilated causal convolution is: Among them, h t x is the local temporal feature vector output by the convolution. seq The input time series vector corresponds to X feat The feature sequence of a specific node (i.e., X) feat [:, i, :]), W represents a convolution operation with an expansion rate of d. f b f For filter parameters, W g b g For gating parameters, ⊙ represents the Hadamard product; Then, the convolution output h t Input a multi-head self-attention layer and compute global dependencies through multi-head self-attention; Finally, output the temporal flow latent feature H. T The formula for calculating the latent features of time flow is: Among them, H T The temporal latent features representing the final output are fused with local convolutional features and global attention features, and their dimension is similar to the spatial flow output features H. S Maintaining consistency is essential for future integration; Furthermore, step S304 specifically includes: First, calculate the fusion gating coefficient z. The formula for calculating the fusion gating coefficient is: Where z is the fusion gating coefficient, || represents feature concatenation, and W z b z For gating network parameters; Then, the final risk latent feature tensor Z is generated using the gating coefficients. risk The formula for calculating the risk latent feature tensor is: Among them, Z risk The final output risk latent feature tensor enables adaptive adjustment of the spatiotemporal feature weights according to different risk types; Furthermore, step S4 specifically includes: S401: Introducing cloud model theory to address the fuzziness and randomness in risk assessment, defining three core numerical characteristics of the risk cloud model: expectation Ex, entropy En, and hyperentropy He; S402: Construct a deep cloud mapping network and convert the high-dimensional risk latent features Z output in step S3 into... risk Mapping to the cloud feature space, calculate the forward cloud generator parameters of pipeline node i at the current time; S403: Based on the Forward Normal Cloud Generator algorithm, Monte Carlo simulation is performed to generate M risk cloud droplets to quantify the risk distribution; S404: Calculate the comprehensive risk assessment index R i and assess confidence level C i The centroid method of cloud droplet swarms is used to transform the uncertain distribution of cloud droplets into a single decision value. Step S401 specifically includes: First, we define Ex, which reflects the central location or average intensity of risk occurrence and represents a qualitative concept of risk level. Then, we define En, which reflects the uncertainty and distribution range of risk performance and represents the ambiguity of qualitative concepts; Finally, we define He, which reflects the uncertainty of entropy, that is, the discreteness or random fluctuation of risk generation; In step S402, the cloud parameter mapping formula is: in, Let be the risk latent feature vector of node i; , , These are the weight matrices for the expectation, entropy, and hyperentropy mapping networks, respectively. , , This is the corresponding bias vector; The sigmoid activation function normalizes the expected value to the interval [0, 1] to match the risk probability definition. The exp(⋅) function ensures that the entropy value is positive. softplus(⋅) = ln(1+e x The function ensures that the superentropy is positive and smooth and differentiable. Through this network, abstract deep learning features are transformed into digital features of cloud models with physical meaning. Furthermore, step S403 specifically includes: First, according to entropy En i and superentropy He i Generate with En i For expectations, He i Normal random entropy with standard deviation The formula for generating normal random entropy is: in, The generated mean is Normal random numbers, This reflects the uncertainty bandwidth of risk in a single sampling. Then, according to the expectation Exi and the generated random entropy Generate a quantitative risk value x m The formula for generating the quantitative risk value is: Where, x m Let m be the specific value of the m-th cloud droplet in the number domain space, representing a specific risk prediction value; Finally, calculate the quantitative value x. m Certainty μ of the concept of "high risk" m The formula for calculating the degree of certainty is: Where, μ m x represents m The probability intensity belonging to the current risk level, with a value range of [0, 1], is repeated M times to obtain the cloud droplet set {(x1, μ1), (x2, μ2), … , (x M , μ M This set constitutes a complete cloud map description of the risk; Furthermore, step S404 specifically includes: Among them, R i The comprehensive risk assessment index for node i is a weighted index that takes into account the magnitude and certainty of each cloud droplet. The formula for calculating the confidence level is: Among them, C i This represents the average degree of certainty of the evaluation results based on the cloud droplet swarm generated from the current M Monte Carlo simulations; if R i Higher and C i A lower R value indicates that while cloud droplets are biased towards high risk, their certainty is dispersed, suggesting high uncertainty in the risk and requiring manual verification; if R... i Higher and C i A high level indicates that the cloud droplets are concentrated and definite, thus indicating a definite high risk and automatically triggering an alarm; Furthermore, step S5 specifically includes: S501: Based on pipeline geographic information system data and building information model, a high-precision three-dimensional geometric model of the pipeline network is constructed, and an octree structure is used for spatial index optimization to support real-time rendering in large-scale pipeline network scenarios. S502: Design a four-dimensional risk field color mapping function Establish a risk index R i The non-linear mapping relationship between node color and transparency; S503: The force-directed algorithm is used to simulate the dynamic diffusion effect of risk energy in the pipeline network. For the high-risk nodes identified, a dynamic particle flow is generated. S504: In the WebGL rendering pipeline, write a custom fragment shader to perform physically based halo rendering on the pipeline walls, visually demonstrating the uncertainty boundary of risks; Furthermore, step S6 specifically includes: S601: Establish a risk classification and determination mechanism, compare the comprehensive risk assessment index R_i calculated in step S404 with the preset multi-level risk threshold vector r, and determine the risk level of the current node; S602: Construct an intelligent focusing linkage mechanism based on a visualization engine; S603: Generates edge collaborative control instruction packages and executes closed-loop feedback; Furthermore, step S602 specifically includes: When the risk level of any node i reaches Alarm or Emergency, the viewpoint control algorithm of the digital twin interface is automatically triggered, driving the virtual roaming camera in the 3D scene constructed in step S501 to move smoothly from the current viewpoint to the optimal observation viewpoint of the abnormal node. The formula for calculating the smooth motion trajectory of the camera is: Among them, Cam pos (t) represents the spatial coordinates and attitude quaternion of the camera at time t; Cam current The current roaming viewpoint; The optimal viewing point for the abnormal node i (usually set to the node normal offset position); Slerp(⋅) is a spherical linear interpolation function used to achieve a smooth and seamless camera movement; Smooth(t) is a easing function to ensure that the camera movement conforms to visual inertia; Furthermore, step S603 specifically includes: Based on the determined risk level and type, a corresponding control strategy is automatically generated. If the level is Warning, an adjustment instruction is generated. Among them, u control For the control quantity issued (such as valve opening degree, pump speed), η is the adjustment step size, and ∇J opt (u) represents the optimization gradient based on the pipeline hydraulic balance model; If the level is Emergency, an emergency cut-off signal is issued directly. The instruction is sent to the field edge computing gateway for execution via the low-latency MQTT protocol. The device status and sensor feedback data after execution are used as the input data for constructing the original spatiotemporal tensor Xraw in the next moment. Step S1 is then restarted to achieve real-time closed-loop verification from perception to control and back to perception.

[0008] The beneficial effects of this invention include: (1) This invention constructs a global pipeline sensing network and uses a weighted covariance alignment algorithm to perform spatiotemporal alignment of multi-source heterogeneous data, which solves the asynchronous problem caused by the time delay of data transmission from distributed sensors and provides a high-quality data foundation for subsequent high-precision analysis. (2) The present invention adopts adaptive variational mode decomposition technology, which can adaptively separate environmental noise and fault features in the signal, effectively solve the problem of complex noise interference in industrial sites, and significantly improve the ability to extract early fault features such as micro-leakage; (3) This invention innovatively constructs a dual-stream spatiotemporal graph convolutional attention network. Through parallel processing and adaptive fusion of spatial and temporal streams, it deeply explores the topological spatial dependence and temporal evolution law of risk propagation in pipeline networks, which greatly improves the accuracy and robustness of risk identification. (4) This invention introduces cloud model theory to construct a risk assessment system, and uses three numerical features, expectation, entropy and hyperentropy, to comprehensively describe the fuzziness and randomness of risk, thereby realizing the scientific quantification of risk uncertainty and providing decision support for managers including confidence level; (5) This invention designs a dynamic three-dimensional visualization engine based on force-guided particle system and halo rendering, which maps abstract risk data into intuitive visual effects, realizes the "what you see is what you get" holographic display of the risk field, and greatly reduces cognitive load; (6) This invention establishes a hierarchical closed-loop response mechanism. Through intelligent focusing linkage and edge collaborative control, it realizes closed-loop management of the entire process from risk perception, assessment, display to disposal, which significantly improves the response speed and intelligence level of pipeline safety management. Attached Figure Description

[0009] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will now be described in further detail with reference to the accompanying drawings, wherein: Figure 1 This is a schematic diagram of the process of a pipeline risk dynamic perception and visualization method based on multidimensional spatiotemporal graph convolution according to the present invention; Figure 2 This is a schematic diagram illustrating the specific operation process of weighted covariance alignment according to the present invention; Figure 3This is a schematic diagram illustrating the specific operation process of the adaptive variational mode decomposition of the present invention; Figure 4 This is a schematic diagram of the dual-stream spatiotemporal network architecture of the present invention; Figure 5 This is a schematic diagram illustrating the specific operation process of the hierarchical closed-loop response of the present invention; Figure 6 This is a schematic diagram of the early warning situation of the pipeline risk perception visualization system in an embodiment of the present invention; Figure 7 This is a schematic diagram of the alarm status of the pipeline risk perception visualization system in an embodiment of the present invention; Figure 8 This is a schematic diagram of an emergency situation in the pipeline risk perception and visualization system according to an embodiment of the present invention. Detailed Implementation

[0010] The preferred embodiments of the present invention will now be described in detail. It should be understood that the preferred embodiments are for illustrative purposes only and are not intended to limit the scope of protection of the present invention.

[0011] like Figure 1 As shown, the present invention provides a method for dynamic perception and visualization of pipeline risks based on multidimensional spatiotemporal graph convolution, comprising the following steps: Step S1: Construct a global pipeline sensing network, collect pipeline multi-physics field operation data and perform spatiotemporal alignment, and construct a multi-source heterogeneous spatiotemporal data tensor; Step S2: Perform adaptive variational mode decomposition and graph signal processing on the spatiotemporal data tensor to construct a weighted graph Laplacian matrix that reflects the pipeline topology; Step S3: Construct a dual-stream spatiotemporal graph convolutional attention network to extract the temporal evolution features and spatial diffusion features of risk in parallel, and fuse them to output a risk latent feature map; Step S4: Construct a risk assessment system based on a cloud model, calculate the expected value, entropy, and hyperentropy of risk, and quantify the dynamic risk probability and confidence level; Step S5: Design a dynamic 3D visualization engine to render the pipeline risk field and anomaly propagation streamlines in real time using a particle system based on risk cloud droplet parameters; Step S6: Establish a hierarchical closed-loop response mechanism and trigger edge control strategies based on the risk level.

[0012] like Figures 2 to 8 As shown, the specific steps of the above method will be further explained below through a specific embodiment.

[0013] In this embodiment, step S1 includes the following sub-steps: Step S101: Based on the physical space and fluid dynamics characteristics, map the physical entities to the mathematical space to construct the pipeline network topology. G=( V , E , A ); In this embodiment, the set of sensor nodes is determined. V The 50 sensors distributed in the pipeline network are mapped as graph nodes, denoted as... V ={ v 1, v 2, ..., v 50}, each node vi It represents the location of a specific physical monitoring point; Then, construct the set of physical connections and fluid dynamics associated edges of the pipeline. E The aim is to capture direct physical connections and indirect hydraulic coupling relationships. The specific determination logic is as follows: (1) Calculate the propagation velocity of pressure waves inside the pipe using the Zhukovsky formula. a The formula is: In this embodiment, the fluid is crude oil, with a bulk modulus of K Take 1.5 × 10 9 Pa, density ρ Take 850kg / m 3 The pipe is made of carbon steel with a modulus of elasticity of [missing value]. Epipe Take 206 GPa, pipe diameter (diam) of 0.5 m, and wall thickness... δ The value is 0.01m. Calculations show that the pressure wave propagation speed is a≈1100m / s.

[0014] (2) Calculate any two nodes vi and vj Fluid transport delay between τij This refers to the time required for the pressure wave to propagate along the pipe path, and its formula is: in, Lij This represents the shortest pipeline distance between nodes; (3) Define the radius of influence in fluid dynamics R hydro That is, within the set data observation time window The formula for the maximum effective propagation distance of an internal pressure wave is: In this embodiment, the data observation time window Set to 2 seconds. Calculated R hydro ≈2420m.

[0015] (4) Construct the edge setE traverse the node set V Determine nodes based on pipeline construction drawings v i and v j Are there physical pipe segments connected between them? If the nodes v i and v j If there is a direct physical pipe connection between them, then a physical edge is established; or if the following conditions are met... That is, equivalent to fluid transport delay If a spatiotemporal causal relationship exists between the two, then an edge is established between them. e ij , forming an edge set E Finally, the weighted adjacency matrix reflecting the medium transmission resistance or distance attenuation factor between nodes is calculated. The adjacency matrix weights are calculated using the following formula: In this embodiment, the distance attenuation constant Set the value to 500, the adjustment factor. γ Take 0.8, coefficient of friction f Take 0.02, gravitational acceleration g Take 9.8 m / s 2 .

[0016] Step S102: Collect data from each node. v i At time step t Multidimensional attribute data to construct the original measurement vector By deploying pressure, flow, temperature, vibration, and acoustic emission sensors, multi-source data are collected simultaneously at the same time. For scalar data such as pressure and flow, Z-Score normalization is used to eliminate dimensional differences, mapping different physical quantities to the same distribution scale. The normalization formula is as follows: in, x These are the raw physical quantities currently being collected. x' These are the standardized eigenvalues. μ window This represents the average value of the sensor over a historical sliding window. The historical standard deviation; In this embodiment, the sliding window size is set to 100 sampling points.

[0017] For high-frequency vibration signals, the Fast Fourier Transform (FFT) is used to extract frequency domain features. The transform formula is as follows: In this embodiment, the number of sampling points N =1024, select amplitude spectrum | X [ k The top 5 largest frequency components are used as features; Finally, the processed data is encapsulated into feature dimensions of F The high-dimensional vector, the original measurement vector is represented as: In this embodiment, feature dimension F Take 8; Step S103: Introduce a weighted covariance alignment algorithm to perform time synchronization correction on the original measurement vectors that are out of sync due to sensor transmission delay; In this embodiment, the specific operation process of weighted covariance alignment is as follows: Figure 2 As shown; First, compute nodes i Its neighboring node set N ( i Each node within) j The spatiotemporal correlation coefficient between data points is measured using the Gaussian kernel function to assess the similarity of data feature distributions. The formula for calculating the correlation coefficient is as follows: In this embodiment, the bandwidth parameter of the Gaussian kernel σ k Set the value to 1.0; Then, the normalized weighted alignment factor is calculated based on the correlation coefficient. w ij To ensure that the sum of the contribution weights of neighboring nodes to the current node is 1, the weighted alignment factor is calculated as follows: Among them, the denominator term For nodes i The sum of the correlation coefficients of all neighboring nodes is used for normalization. Finally, a weighted alignment factor is used to aggregate the data of neighboring nodes to correct the observation value of the current node. The weighted covariance alignment correction formula is as follows: In this embodiment, the alignment strength adjustment coefficient α Take 0.6; Step S104: Perform time window slicing on the multi-source data to construct a network-level spatiotemporal input tensor. X raw ; Set the historical observation time window length to Δ T Using a sliding window mechanism, all network nodes are displayed.N A node in the past T in The measurement vectors at each time step are stacked along the time dimension to form a three-dimensional spacetime tensor: in, for t A real-time snapshot matrix of the entire network's status; In this embodiment, the sampling frequency Δ t 50Hz, time step It is 100; Step S2 includes the following sub-steps: Step S201: Adaptive variational mode decomposition technology is used to decompose the non-stationary time series signals collected by each node. f ( t Automatically decomposed into K An eigenmode function with a specific center frequency; In this embodiment, the specific operation process of adaptive variational mode decomposition is as follows: Figure 3 As shown; First, the input object and feature target are clearly defined. In this embodiment, nodes are used as the input object. v 10 Taking the instantaneous pressure signal as an example, f ( t The pressure is the raw instantaneous pressure collected in step S102 within the historical time window T. in A continuous signal sequence assembled in time sequence, specifically, defining the current time. t Historical observation vector P window : in, The number of sampling points within the time window. The sampling interval; To meet the requirements of AVMD for continuous domain signal processing, Shannon interpolation formula is used to reconstruct discrete vectors into continuous-time functions. f ( τ The formula for splicing and reconstructing is: in, sinc ( x )= sin ( πx ) / ( πx ) is the interpolation kernel function, which transforms discrete sensor readings into a smooth signal containing transient pressure waveforms and nonlinear vibration modes; Then, the analytical form of the signal is constructed to obtain its spectral characteristics. The Hilbert transform is then used to convert the real signal into a complex analytic signal. The Hilbert transform formula is: in, Represents the Hilbert transform operator. p . v . represents the Cauchy principal value integral. For the first k One modal component; Based on this, the formula for constructing an analytical signal is: in, z k ( t The signal is an analytic signal with a one-sided spectral characteristic (containing only positive frequencies). δ ( t ) is the Dirac distribution function. , j The imaginary unit; Furthermore, we construct a variational optimization problem, assuming the original signal is given by... K A frequency around the center ω k AM / FM signals uk ( t The signals are superimposed to form a matrix. To minimize the bandwidth of each mode, the spectrum of the analytic signal is modulated to the baseband using an exponential shift. The objective functional formula for AVMD optimization is: in, J AVMD This represents the total bandwidth target to be minimized. ω k This refers to the specific center frequency corresponding to this mode, which is a complex exponential twitch factor used to perform spectrum shifting operations, changing the center frequency... ω k Move to zero frequency, ∂t This is a time gradient operator used to calculate the smoothness (i.e., bandwidth) of the demodulated baseband signal. δ ( t ) is the Dirac distribution function. Constraints Ensure the decomposition K The sum of the modalities can losslessly restore the original signal; In this embodiment, the number of modes is set. K =4, which aims to extract four main components: fluid flow, mechanical vibration, environmental noise, and potential leakage impact.

[0018] To solve this constrained variational problem, Lagrange multipliers are introduced. λ ( t ) and secondary penalty factor α Transforming it into an unconstrained optimization problem, the augmented Lagrange function formula is: in, L To augment the Lagrange function; the first term of the formula This is a secondary penalty item. α The penalty parameter, used to balance data fidelity and bandwidth constraints, is employed to ensure reconstruction accuracy while constraining modal bandwidth; the second term of the formula... The residual fidelity term represents the energy difference between the original signal and the sum of the decomposed modes, i.e., the reconstruction error; the third term in the formula... For Lagrange multipliers, This indicates the inner product operation, used to strictly enforce the satisfaction of equality constraints; In this embodiment, the secondary penalty factor α Take 2000; Furthermore, the alternating direction multiplier method is used to alternately update the variables in the frequency domain until convergence is obtained to obtain the four eigenmode functions. u k ( t The specific iterative steps are as follows: (1) The modal component update formula is: (2) The center frequency update formula is: (3) The Lagrange multiplier update formula is: in, ; In this embodiment, the multiplier update step size ζ Set the value to 0.01 and iterate using the alternating direction multiplier method until the convergence condition is met. satisfy.

[0019] Finally, a signal reconstruction operation is performed, and a subset of modes containing valid fault information is selected based on the kurtosis criterion. S valid ; In this embodiment, the kurtosis threshold is set to 3.0. The kurtosis values ​​of the four decomposed modes are as follows: K u1 =2.8, K u2 =4.5, K u3 =3.2,K u4 =2.9. Therefore, the effective modal subset S valid = { u 2, u 3}.

[0020] Noisy modes are removed, and the clean signal is reconstructed to build a denoised feature tensor. X feat The formula for denoising feature reconstruction is: in, For nodes i At any moment t The denoised feature vector, This represents the superposition of effective modes, indicating the denoised dynamic signal (such as a pressure wave). For static attribute features that do not require decomposition (such as temperature, flow rate baseline); Based on this, construct the denoising feature tensor X feat This tensor is obtained by placing all nodes in a time window. T in Internally reconstructed feature vector Stacked along the time axis, the formula for constructing the denoising feature tensor is: in, for t The state snapshot matrix after full network denoising at any given moment; this tensor Xfeat This will be used as the input data for the dual-stream network in step S3; Step S202: The weighted adjacency matrix constructed based on step S101 W From the perspective of graph theory, a normalized graph Laplace matrix reflecting the topological structure and functional dependencies of a pipeline system is constructed. L ; The specific steps for constructing the graph Laplacian matrix in step S202 are as follows: First, calculate the degree matrix. D It is a diagonal matrix, with diagonal elements Dii Weighted adjacency matrix W No. i The sum of all elements in a row is calculated using the following formula: In this embodiment, node v 10Three neighbors are connected, with pipe distances of 22.0m, 50.0m, and 68.0m respectively. Substituting into the weight formula in step S101, the weights are 0.8, 0.6, and 0.5 respectively. Then the diagonal elements of the degree matrix... D ii =0.80 + 0.60 + 0.50 = 1.90. Then, construct the degree matrix D, which is composed of... D ii The constructed diagonal matrix, that is, the elements excluding the main diagonal are D ii Except for the elements at all other positions, all other elements are 0, satisfying the condition. D =diag ( D 11 ,…, D NN ); Finally, the adjacency matrix is ​​symmetrically normalized using the degree matrix to construct the graph Laplacian matrix, which is calculated as follows: in, I N It is a 50th order identity matrix. W This is the weighted adjacency matrix in step S101. The inverse square root matrix of the degree matrix, and the normalized Laplace matrix. L It possesses the property of being positive semidefinite, and its eigenvalues λ i Distributed within the interval [0,2], this ensures the numerical stability and convergence of subsequent graph convolution operations; Step S203: Based on the graph Laplacian matrix L A graph filter is designed using Chebyshev polynomial approximation theory to project the spatiotemporal signal of the pipeline into the graph frequency domain. The specific steps of the Chebyshev polynomial in step S203 are as follows: First, calculate the rescaled Laplacian matrix, mapping the eigenvalues ​​to the interval [-1, 1] to satisfy the domain requirement of the Chebyshev polynomial. The rescaling formula is: in, This is the rescaled Laplace matrix. λ max Laplace matrix L The largest eigenvalue; Then, a Chebyshev polynomial sequence is defined using recursion to quickly compute feature extraction in , and the recursive formula is: in, fork Chebyshev polynomial operators; In this embodiment, the polynomial order k Taking 3, the generated operator sequence is: T 0, T 1, T 2; Step S3 includes the following sub-steps: In this embodiment, the dual-stream spatiotemporal network architecture is as follows: Figure 4 As shown; Step S301: Construct a two-stream feature extraction architecture and extract the denoised feature tensor output from step S2. Xfeat Simultaneously inputting into two parallel processing branches: spatial graph convolutional stream and temporal attention stream, the aim is to capture the topological spatial dependency and temporal evolution of risk propagation in pipeline networks, respectively, and improve the feature extraction efficiency of the model by decoupling spatiotemporal features; In this embodiment, the input tensor Xf eat The dimension is ( T in , N , F ),in T in =100, N =50, F =8; Step S302: In the spatial graph convolutional flow, design a multi-level dynamic graph convolutional layer, which uses tensors... X feat Taking the spatial slices at each time step as input, convolution operation is performed by combining static physical topology with dynamic data association, and feature aggregation is performed using the Chebyshev polynomial generated in step S203. The specific steps for feature aggregation of the spatial graph convolutional flow in step S302 are as follows: First, calculate the spatial self-attention matrix. Aatt This matrix reflects the similarity weights between the features of each node at the current time. The formula for calculating the spatial self-attention matrix is: in, Aatt For spatial self-attention matrix, X (l) For the first l The input feature matrix, when l When =0, X (0) is the denoised feature vector input in step S301. X feat At the present moment t Spatial feature slices, i.e. X feat[t, :, :], W Q and W K The weight matrix is ​​a learnable linear mapping. d k The feature dimension scaling factor; In this embodiment, the feature dimension scaling factor d k Take 64, weight matrix .

[0021] Then, construct the fusion graph Laplacian operator. L mix The physical topology and data-driven association are weighted and fused together, and the fusion formula is as follows: in, L mix To fuse graph operators, This is the static Chebyshev reference matrix. γ The learnable balance coefficient; In this embodiment, the initial balance coefficient γ Take 0.5; Finally, calculate the spatial stream output characteristics. The formula for calculating the spatial flow output characteristics is: in, For the first l +1 layer output features, T k (⋅) represents the Chebyshev polynomial operator defined in step S203. For the first k The convolution kernel parameter matrix of order 1. It is the ReLU activation function; In this embodiment, the convolution kernel weight matrix After two layers of graph convolution, the output spatial flow feature tensor is obtained. H S The dimensions become (100, 50, 64); Step S303: In the temporal attention stream, design a causal gated temporal convolutional network and a multi-head temporal self-attention mechanism to capture long-short-term temporal dependencies; The specific steps in step S303 for capturing long- and short-term time dependencies using time attention flow are as follows: First, local temporal features are extracted using dilated causal convolution. This step converts the tensor... X feat Unfolding along the node dimension, convolving the time series of each node, the output formula for dilated causal convolution is: in, h t This is the local temporal feature vector output by the convolution. x seq The input time series vector corresponds to X feat The feature sequence of a specific node, i.e. X feat [:, i , :], Indicates the expansion rate d Convolution operation, W f , b f For filter parameters, W g , b g For gating parameters, ⊙ represents the Hadamard product; In this embodiment, the causal-gated temporal convolutional network contains three residual blocks with a dilation rate of [missing information]. d Set them to 1, 2, and 4 respectively; Filter W f and gate controller W g The kernel size is set to 3, and the input and output channels are mapped to 64. Then, output the convolution. h t Input a multi-head self-attention layer, and calculate global dependencies through multi-head self-attention. The multi-head self-attention calculation formula is as follows: in, MHSA (⋅) represents the output of the multi-head self-attention module. head i Indicates the first i The output of each attention head, Q i , K i , V i It is by h t The query, key, and value matrix obtained by linear transformation M mask This is a causal mask matrix to ensure that future information is not leaked. W O To output the projection matrix; In this embodiment, the number of attention heads is set. h=4, and each head has a dimension of 16; Finally, output the temporal flow latent features. H T The formula for calculating the latent features of time flow is: in, H T The temporal latent features representing the final output are fused with local convolutional features and global attention features, and their dimensionality and spatial flow output features are similar. H S Keep them consistent, all being (100, 50, 64); Step S304: Design a spatiotemporal fusion gating unit to adaptively fuse spatial flow features. H S With time flow characteristics H T ; The specific steps for fusing spatial flow features and temporal flow features in step S304 are as follows: First, calculate the fusion gating coefficient. z The formula for calculating the fusion gating coefficient is: in, z To fuse the gating coefficients, || denotes feature concatenation. W z , b z For gating network parameters; In this embodiment, the feature dimension after splicing is 128. .

[0022] Then, the final risk latent feature tensor is generated using the gating coefficients. Z risk The formula for calculating the risk latent feature tensor is: in, Z risk The final output risk latent feature tensor enables adaptive adjustment of the spatiotemporal feature weights according to different risk types; In this embodiment, the final output risk latent feature tensor Z risk The dimensions are (100, 50, 64); Step S4 includes the following sub-steps: Step S401: Introduce cloud model theory to address the fuzziness and randomness in risk assessment, and define three core numerical characteristics of the risk cloud model: expectation. Ex ,entropy En and hyperentropy He ; The specific steps of step S401 are as follows: First, define Ex It reflects the central location or average intensity of the risk occurrence and represents a qualitative concept of risk level; Then, define En It reflects the uncertainty and distribution range of risk performance, representing the ambiguity of qualitative concepts; Finally, define He It reflects the uncertainty of entropy, that is, the degree of dispersion or random fluctuation in risk generation; Step S402: Construct a deep cloud mapping network and convert the high-dimensional risk latent features output in step S3 into... Z risk Mapped to the cloud feature space, calculate the pipeline node at the current time. i Parameters of the positive cloud generator; The specific steps for constructing the deep cloud mapping network in step S402 are as follows: in, For nodes i The risk hidden feature vector; These are the weight matrices for the expectation, entropy, and hyperentropy mapping networks, respectively. This is the corresponding bias vector; The expected value is normalized to the [0, 1] interval to match the definition of risk probability. exp The (⋅) function ensures that the entropy value is positive. The function ensures that the superentropy is positive and smoothly differentiable. Through this network, abstract deep learning features are transformed into digital features of cloud models with physical meaning. In this embodiment, the deep cloud mapping network consists of three parallel fully connected layers. The expected prediction layer weights... bias Entropy prediction layer weights bias ; Hyperentropy prediction layer weights bias To monitor nodes v 10 For example, after network forward propagation calculation, the cloud digital feature parameters obtained are: expected value Ex 10 =0.82 indicates that the node has a relatively high average risk; entropy En 10 =0.15 indicates that the risk has a certain range of uncertainty; hyperentropy He 10 =0.05 indicates that the entropy itself has small random fluctuations, meaning that the cloud droplet thickness is relatively thin.

[0023] Step S403: Perform Monte Carlo simulation based on the Forward Normal Cloud Generator algorithm to generate... M A risk cloud droplet To quantify risk distribution; The specific steps for Monte Carlo simulation using the forward normal cloud algorithm in step S403 are as follows: First, according to entropy En i and hyperentropy He i Generate with En i For expectations, He i Normal random entropy with standard deviation The formula for generating normal random entropy is: in, Norm ( μ , σ 2 ) indicates that the generated mean is The normally distributed random numbers reflect the uncertainty bandwidth of risk in a single sampling. Then, according to expectations Ex i and the generated random entropy Generate quantitative risk values x m The formula for generating the quantitative risk value is: in, x m For the first m The specific value of a cloud droplet in the data domain space represents a specific risk prediction value. Finally, the quantitative value is calculated. x m Certainty regarding the concept of "high risk" μ m The formula for calculating the degree of certainty is: in, μ m express x m The probability strength belonging to the current risk level, with a value range of [0, 1], is determined by repeating the above process. M Next, we obtain the cloud droplet set {( x 1, μ 1), ( x 2, μ 2), … , ( x M , μ M This set constitutes a complete cloud map description of the risk; In this embodiment, the number of Monte Carlo simulations M Take 1000. (Based on node...) v 10 For example, the generation process and calculation results of the first 5 cloud droplets are shown: (1) Sample 1: Generate standard normal random numbers N 1a =0.4, N 1b =-0.2. Then En =0.15 + 0.4 × 0.05 = 0.17, Risk Value x 1 = 0.82 + (-0.2) × 0.17 = 0.786, Determinism μ 1= exp (-(0.786-0.82) 2 / (2×0.17 2 ))≈0.980.

[0024] (2) Sample 2: Generate standard normal random numbers N 2a =-0.6, N 2b =0.5. Then En =0.15 + (-0.6) × 0.05 = 0.12, Risk Value x 2 = 0.82 + 0.5 × 0.12 = 0.880, Determinism μ 2= exp (-(0.880-0.82) 2 / (2×0.12 2 ))≈0.980.

[0025] (3) Sample 3: Generate standard normal random numbers N 3a =1.2, N 3b =-0.8. Then En =0.15 + 1.2 × 0.05 = 0.21, Risk Value x 3 = 0.82 + (-0.8) × 0.21 = 0.652, Determinism μ 3= exp (-(0.652-0.82) 2 / (2×0.21 2 ))≈0.726.

[0026] (4) Sample 4: Generate standard normal random numbers N 4a =-0.2, N 4b =0.1. Then En =0.15 + (-0.2) × 0.05 = 0.14, Risk Value x 4 = 0.82 + 0.1 × 0.14 = 0.834, Determinism μ 4= exp (-(0.834-0.82) 2 / (2×0.14 2 ))≈0.995.

[0027] (5) Sample 5: Generate standard normal random numbers N 5a =0.0, N 5b =-1.5. Then En =0.15 + 0.0 × 0.05 = 0.15, Risk Value x 5 = 0.82 + (-1.5) × 0.15 = 0.595, Determinism μ 5= exp (-(0.595-0.82) 2 / (2×0.15 2 ))≈0.325.

[0028] Step S404: Calculate the comprehensive risk assessment index Ri and assess confidence level Ci The centroid method of cloud droplet swarms is used to transform the uncertain distribution of cloud droplets into a single decision value. In step S404, the risk assessment index is calculated. R i and assess confidence level C i The specific steps are as follows: in, R i For nodes i The comprehensive risk assessment index is a weighted average of the magnitude and certainty of each cloud droplet. The formula for calculating the confidence level is: in, C i Indicates based on the current M The average degree of certainty of the cloud droplet swarms generated by the Monte Carlo simulation for the evaluation results; if Ri higher and C i A lower value indicates that while the cloud droplet is biased towards high risk, its certainty is dispersed, suggesting high uncertainty in the risk and requiring manual verification; if R i higher and C i A high level indicates that the cloud droplets are concentrated and definite, thus indicating a definite high risk, and the system can automatically trigger an alarm. In this embodiment, for M =Accumulate and calculate 1000 cloud droplets. , Then the node v 10 Comprehensive risk assessment index R 10 =750.5 / 920.0≈0.815, assess confidence level C 10 =920.0 / 1000=0.92. This result indicates that the node... v 10 Currently in a high-risk state R 10 >0.8, and this evaluation result C 10 =0.92 has an extremely high confidence level and should trigger a high-level alarm immediately.

[0029] Specifically, step S5 also includes the following sub-steps: Step S501: Based on pipeline geographic information system data and building information model, construct a high-precision three-dimensional geometric model of the pipeline network, and use an octree structure for spatial index optimization to support real-time rendering in large-scale pipeline network scenarios. The specific steps for constructing the three-dimensional geometric model in step S501 are as follows: First, multi-source heterogeneous data is analyzed to obtain the pipeline's geometric and attribute information; specifically, pipeline geographic information system (GIS) data refers to the geospatial coordinate sequence of the pipeline centerline. The data in the Building Information Modeling (BIM) system specifically refers to the cross-sectional properties of the pipes, the material texture mapping coordinates, and the detailed 3D mesh model of the ancillary facilities (valves, pump stations). Then, a unified spatial coordinate system was established, converting the geographic coordinates of the GIS data into Cartesian rectangular coordinates used by the rendering engine, for high-risk nodes. v 10 For example, let's assume its geographical coordinates are... Using Web Mercator projection, the projection reference point (0,0) corresponds to (104.0). ∘ E, 30.0∘ N), scale factor s=111000, approximately 111km per degree; coordinate transformation formula is: in,[ x k , y k , z k ] T These are the transformed local coordinates; M proj The Mercator projection matrix is ​​used to project spherical coordinates onto a plane. T offset This is the scene center offset vector, used to solve the problem of precision loss in large coordinate floating-point numbers; In this embodiment, T offset Set the scene center offset vector to [-6200, -73000, 0]; Substituting the specific values, the calculation is: x = (104.065 - 104.0) × 111000 × cos(30.65) ∘ y = (30.658 - 30.0) × 111000 ≈ 73038.0 m, z = 520 m, the transformed local coordinates are [5.5, 38.0, 520]. T ; Furthermore, a 3D mesh of the pipeline is generated. For each pipeline segment, the mesh is based on its centerline coordinates. x k , y k , z k )and( x k+1 , y k+1 , z k+1 and pipe diameter diam Parametric modeling is used to generate a set of triangular facets for a cylinder. The formula for generating vertices is: in, v ( θ , h () represents the coordinates of the vertex on the surface of the cylinder. c ( h ) is the parameter on the center line h The point at that location, n , b Let these be the normal vector and the binormal vector at that point. θ ∈[0,2π) is an inscribed angle; In this embodiment, the number of circumferential segments is set to 36; Finally, an octree spatial index is constructed, and the generated millions of triangular facet primitives are recursively divided into eight child nodes to establish hierarchical bounding boxes. The virtual camera is then initialized, and the formula for calculating the axis-aligned bounding box of each node is as follows: in, v set The set of all vertices contained in this node; In this embodiment, the maximum depth of the octree is set to 8, and a split is triggered when the number of triangular faces contained in a leaf node exceeds 1000. Meanwhile, a virtual roaming camera is constructed in the 3D scene, and its intrinsic parameter matrix and initial extrinsic parameter matrix are defined. This camera serves as the controlled object for intelligent focusing linkage in the subsequent step S6, and is used to generate the final rendering viewport. Step S502: Design a four-dimensional risk field color mapping function Establish a risk index Ri The non-linear mapping relationship between node color and transparency; The node color mapping formula in step S502 is: in, Color i For nodes i An RGBA color vector containing red, green, blue, and alpha components. Color safe The default safety status color. Color danger The preset danger state color, R i The nodes calculated in step S404 i The comprehensive risk assessment index, r low , r high These are the low and high thresholds for risk visualization, respectively. This is a smoothing step function used to... R i exist[ r low , r high The smooth interpolation coefficients that map to [0, 1] within the interval. It is a linear interpolation function that mixes safety and danger colors based on the interpolation coefficients, intuitively converting risk values ​​into visual color signals; In this embodiment, the preset r low =0.3,r high =0.8. Color safe The value is green (0, 1, 0, 0.5). Color danger The value is red (1,0,0,1.0). For the node... v 10 Its comprehensive risk index R 10 =0.815. Because R 10 > r high ,but SmoothStep (0.3, 0.8, 0.815) = 1.0. Substitute into the interpolation formula: Color 10 = Lerp (Green,Red,1.0)=(1,0,0,1.0), meaning the node is pure red and completely opaque.

[0030] Step S503: Use the force-directed algorithm to simulate the dynamic diffusion effect of risk energy in the pipeline network. For the high-risk nodes identified, generate dynamic particle flow. The specific steps for generating the dynamic particle flow in step S503 are as follows: First, the motion of the particle flow is influenced by the combined effects of risk gradient force, topological constraint force, and random perturbation force. The particle generation rate is defined as follows: for the risk index... R i Nodes with a value >0.8 generate a certain number of particles per second. N p =100× Ri .for v 10 It generates 81 particles per second, and the dynamic equation of motion of the particles is: in, p The position vector of the particle in three-dimensional space ( x , y , z ), where is the particle's acceleration vector, and is the particle's velocity vector. F risk The risk gradient force is the primary driving force behind particle motion. k d This is the motion damping coefficient, used to simulate fluid viscous drag and prevent particle velocity from increasing indefinitely. F topology Topological constraints restrict particles to tangential movement along the pipe axis, preventing them from breaking through the mold and overflowing the pipe wall. F randomTo simulate the random forces of Brownian motion in turbulence, and to increase the naturalness and randomness of the visual effect; Step S504: In the WebGL rendering pipeline, write a custom fragment shader to perform physically based halo rendering on the pipeline walls, visually demonstrating the uncertainty boundary of the risk. The formula for calculating halo intensity in step S504 is: in, I glow The final halo intensity of the current fragment. α Based on the luminous intensity coefficient, the overall brightness is controlled. En i The node risk entropy calculated in step S402 is as follows: the larger the risk entropy, the stronger and wider the halo, which intuitively represents the "uncertainty range of risk". n Let be the normal vector of the pipe wall surface. v Let be the gaze vector from the fragment to the camera. This is the Fresnel effect term. p The edge attenuation index, this factor makes the tube wall edges exhibit a stronger luminous effect, enhancing the three-dimensional volumetric feel. The item introduces a "breathing" effect over time; ω The respiratory rate is set to be correlated with the risk index. R i Proportional ; In this embodiment, the basic luminous intensity coefficient α Take 1.5, edge decay index p Use version 3.0. (For the node) v 10 respiratory rate This manifests as rapid, flickering breathing. Assuming the current moment... t =2.5s, entropy value En 10 =0.15. If the camera's gaze vector... v with pipe wall normal n cosine value of the included angle The instantaneous intensity of the halo. I glow =1.5×0.15×(0.8) 3 ×sin(5.12×2.5)≈0.019.

[0031] In this embodiment, step S6 further includes the following sub-steps: In this embodiment, the specific operational procedures for hierarchical closed-loop response and the pipeline risk perception visualization system are as follows: Figure 5 and Figures 6-8 As shown; Step S601: Establish a risk classification and determination mechanism by combining the comprehensive risk assessment index calculated in step S404 with the preset multi-level risk threshold vector. r Compare the data to determine the risk level of the current node; The risk level determination logic in step S601 is as follows: Among them, Level i For nodes i Risk level status; r =[ r 1, r 2, r [3] is an increasing risk threshold vector, corresponding to the trigger limits for early warning, alarm, and emergency shutdown, respectively; In this embodiment, the threshold vector r Set to [0.4, 0.7, 0.9]; Step S602: Construct an intelligent focusing linkage mechanism based on a visualization engine; The specific steps for constructing the intelligent focusing linkage mechanism based on the visualization engine in step S602 are as follows: When any node is detected i When the risk level reaches Alarm or Emergency, the viewpoint control algorithm of the digital twin interface is automatically triggered, driving the virtual roaming camera in the 3D scene constructed in step S501 to move smoothly from the current viewpoint to the optimal observation viewpoint of the abnormal node. The formula for calculating the smooth motion trajectory of the camera is: Among them, Cam pos ( t ) for t Spatial coordinates and attitude quaternions of the time-lapse camera; Cam current This refers to the current roaming viewpoint; it also refers to the viewpoint for abnormal nodes. i The optimal viewing point (usually set to the node normal offset position); It is a spherical linear interpolation function used to achieve a smooth, stutter-free push-pull effect; Smooth ( t () is a easing function to ensure that the camera movement process conforms to visual inertia; at the same time, in specific implementation, it is also possible to set up an automatically expanding multi-dimensional data holographic dashboard to highlight the real-time pressure waveform and risk cloud droplet distribution map of the node to assist manual review; In this embodiment, v 10When in Alarm state, the camera is triggered to move, and the camera is currently in roaming view. Cam current =[200,100,100], according to S501 v 10 Local coordinates P node =[5.5,38.0,520], set the optimal observation offset. D view =[20,20,20], then the target location Cam target = P node + D view =[25.5,58.0,540]; Smooth ( t The EaseInOutCubic easing function is used, with a motion duration set to 1.5 seconds. At the halfway point of the motion, the easing function value k≈0.5. The camera position will smoothly transition to approximately [112.75,79.0,320], while simultaneously adjusting the quaternion rotation angle to align with the node. v 10 ; Step S603: Generate edge collaborative control instruction package and execute closed-loop feedback; The specific steps in step S603 are as follows: Based on the determined risk level and type, a corresponding control strategy is automatically generated. If the level is Warning, an adjustment instruction is generated. in, u control For the control quantities issued (such as valve opening degree, pump speed), η To adjust the step size, ∇ J opt ( u () represents the optimization gradient based on the pipeline hydraulic balance model; If the emergency level is set to Emergency, an emergency shutdown signal is issued directly. The command is sent to the field edge computing gateway via the low-latency MQTT protocol for execution, and the device status and sensor feedback data after execution are used to construct the original spatiotemporal tensor for the next moment. Xraw The input data is used to re-enter step S1, realizing real-time closed-loop verification from perception to control and back to perception.

[0032] In this embodiment, the learning rate / step size is set. η =0.05.

[0033] In this embodiment, the pipeline multiphysics monitoring data underwent thorough spatiotemporal alignment and noise reduction preprocessing. Figure 2 and Figure 3 The working principles of the weighted covariance alignment algorithm in multi-source data synchronization and adaptive variational mode decomposition in signal denoising and separation are respectively demonstrated. Figure 4 The dual-stream spatiotemporal network structure in the text can effectively utilize spatial topology and temporal dependence to deeply mine pipeline risk characteristics; Figure 5 This is the workflow of the hierarchical closed-loop response mechanism of the present invention; Figure 6 , Figure 7 and Figure 8 This invention demonstrates an example of pipeline risk visualization and classification based on digital twins. The invention improves the ability to identify minute leaks and early failures through multi-dimensional spatiotemporal feature fusion and cloud model evaluation mechanism, which helps to increase the accuracy and real-time perception of pipeline safety operation and includes more intuitive decision data information.

[0034] Specific embodiments have been used to illustrate the principles and implementation methods of this invention. The descriptions of these embodiments are merely illustrative and are intended to aid in understanding the method and core ideas of this invention. Those skilled in the art will recognize that the embodiments described herein are for the purpose of helping readers understand the principles of this invention and should be understood as not limiting the scope of protection of this invention to such specific statements and embodiments. Those skilled in the art can make various other specific modifications and combinations based on the technical teachings disclosed in this invention without departing from the essence of this invention, and these modifications and combinations are still within the scope of protection of this invention.

Claims

1. A method for dynamic perception and visualization of pipeline risks based on multidimensional spatiotemporal graph convolution, characterized in that: Includes the following steps: Step S1: Construct a global pipeline sensing network, collect pipeline multi-physics field operation data and perform spatiotemporal alignment, and construct a multi-source heterogeneous spatiotemporal data tensor; Step S2: Perform adaptive variational mode decomposition and graph signal processing on the spatiotemporal data tensor to construct a weighted graph Laplacian matrix that reflects the pipeline topology; Step S3: Construct a dual-stream spatiotemporal graph convolutional attention network to extract the temporal evolution features and spatial diffusion features of risk in parallel, and fuse them to output a risk latent feature map; Step S4: Construct a risk assessment system based on a cloud model, calculate the expected value, entropy, and hyperentropy of risk, and quantify the dynamic risk probability and confidence level; Step S5: Design a dynamic 3D visualization engine to render the pipeline risk field and anomaly propagation streamlines in real time using a particle system based on risk cloud droplet parameters; Step S6: Establish a hierarchical closed-loop response mechanism and trigger edge control strategies based on the risk level.

2. The pipeline risk dynamic perception and visualization method based on multidimensional spatiotemporal graph convolution as described in claim 1, characterized in that: Step S1 specifically includes: Step S101: Based on the physical space and fluid dynamics characteristics, map the physical entities to the mathematical space to construct the pipeline network topology graph G=(V, E, W); Step S102: Collect data from each node v i Construct the original measurement vector from the multidimensional attribute data at time step t. The specific operation is as follows: Multiple source data are simultaneously collected at the same time using deployed pressure, flow, temperature, vibration, and acoustic emission sensors. For scalar data, including pressure and flow, Z-Score normalization is used to eliminate dimensional differences, mapping different physical quantities to the same distribution scale. The normalization formula is: Where x is the original physical quantity currently acquired, x' is the standardized feature value, and μ window Let σ be the mean value of the sensor within the historical sliding window. window The standard deviation is used to eliminate the impact of absolute numerical differences on the convergence speed of subsequent gradient descent. For high-frequency vibration signals, the Fast Fourier Transform (FFT) is used to extract frequency domain features. FFT is an efficient algorithm for Discrete Fourier Transform, and its transformation formula is as follows: Where x[n] is the time-domain discrete vibration signal sequence, N is the number of sampling points, X[k] is the frequency-domain complex sequence, k is the frequency index, and j is the imaginary unit; Based on this, the top k dominant frequency components with the largest amplitude spectrum |X[k]| are selected as features; Finally, the processed data is encapsulated into a high-dimensional vector with feature dimension F, and the original measurement vector is represented as: in, These are instantaneous pressure, flow rate, and temperature observations, respectively. These are the historical mean and standard deviation of the corresponding physical quantity, respectively. This represents the amplitude values ​​of the first k dominant frequencies of the vibration signal. The logarithmic value representing the energy of the acoustic emission signal is used to amplify the characteristics of minute leaks. Step S103: Introduce a weighted covariance alignment algorithm to perform time synchronization correction on the original measurement vectors that are out of sync due to sensor transmission delay. The specific calculation steps are as follows: First, calculate the spatiotemporal correlation coefficient between node i and each node j in its neighboring node set N(i). Use the Gaussian kernel function to measure the similarity of data feature distributions. The correlation coefficient calculation formula is as follows: Among them, s ij This is the similarity value output by the Gaussian kernel function, with a value range of (0,1]. The squared Euclidean distance between the feature vectors of two nodes reflects the absolute difference in data distribution. The exponential structure reflects the characteristics of the radial basis function: when the difference approaches 0, the similarity approaches 1; as the difference increases, the similarity rapidly decreases. The bandwidth parameter of the Gaussian kernel is used to control the sensitivity of similarity to distance changes; Then, the normalized weighted alignment factor w is calculated based on the correlation coefficient. ij To ensure that the sum of the contribution weights of neighboring nodes to the current node is 1, the weighted alignment factor is calculated as follows: Among them, the denominator term This is the sum of the correlation coefficients of all neighboring nodes of node i, used for normalization. Finally, a weighted alignment factor is used to aggregate the data of neighboring nodes to correct the observation value of the current node. The weighted covariance alignment correction formula is as follows: in, Here is the corrected node feature vector, and α is the alignment strength adjustment coefficient. The fluid transport delay calculated in step S101, This method utilizes physical delay to perform spatiotemporal backtracking of neighborhood data, eliminating asynchronous errors caused by transmission delay. Step S104: Perform time window slicing on the multi-source data to construct the network-level spatiotemporal input tensor X. raw The specific process is as follows: Set the length of the historical observation time window to T. in Using a sliding window mechanism, all N nodes in the network are viewed in the past T... in Measurement vectors at each time step Stacked along the time dimension, they ultimately form a three-dimensional spacetime tensor: in, Let be the state snapshot matrix of the entire network at time t.

3. The pipeline risk dynamic perception and visualization method based on multidimensional spatiotemporal graph convolution according to claim 2, characterized in that: In step S101, the step of constructing the pipeline network topology graph G=(V, E, W) is as follows: First, determine the set V of sensor nodes, and map the N sensors distributed in the pipeline network as graph nodes, denoted as V={v1,v2, ..., v N }, each node v i It represents the location of a specific physical monitoring point; Then, a set E of physical connections and fluid dynamic relationships between pipelines is constructed to capture direct physical connections and indirect hydraulic coupling relationships. The specific determination logic is as follows: (1) The propagation velocity 'a' of the pressure wave inside the pipe is calculated using the Zhukovsky formula, which is: Where K is the fluid bulk modulus, ρ is the fluid density, and E pipe δ represents the elastic modulus of the pipe material, where diam is the pipe diameter and δ is the pipe wall thickness. (2) Calculate v for any two nodes i and v j Fluid transport delay τ between ij This refers to the time required for the pressure wave to propagate along the pipe path, and its formula is: Among them, L ij This represents the shortest pipeline distance between nodes; (3) Define the radius of influence R in fluid dynamics hydro That is, the maximum effective propagation distance of the pressure wave within the set data observation time window ΔT, and its formula is: Wherein, κ is the safety redundancy coefficient, with a value of 1.0~1.2, used to compensate for errors caused by uneven flow velocity; (4) Construct the edge set E, traverse the node set V, and determine the node v based on the pipeline construction drawings. i With v j Are there physical pipe segments connected between them? If node v i With v j If there is a direct physical pipe connection between them, then a physical edge is established; or if L is satisfied... ij ≤R hydro That is, equivalent to fluid transport delay If a spatiotemporal causal relationship exists between the two, then an edge e is established between them. ij This forms the edge set E; Finally, calculate the weighted adjacency matrix that reflects the medium transmission resistance or distance attenuation factor between nodes. The adjacency matrix weights are calculated using the following formula: Among them, W ij The first term represents the weight of node association strength, and is a Gaussian decay term based on Euclidean distance. The first term is the distance attenuation constant, the second term is the hydraulic friction loss attenuation term based on the Darcy-Weisbach formula, γ is the adjustment factor, f is the friction coefficient, and v flow The formula, where is the average fluid velocity and g is the gravitational acceleration, comprehensively measures spatial proximity and hydraulic connectivity.

4. The pipeline risk dynamic perception and visualization method based on multidimensional spatiotemporal graph convolution as described in claim 2, characterized in that: Step S2 specifically includes: Step S201: Adaptive variational mode decomposition (VMD) technology is used to decompose the non-stationary time series signal f(t) collected by each node into K eigenmode functions with specific center frequencies to separate environmental noise and fault characteristics. The specific decomposition steps are as follows: First, the input object and feature target are defined. f(t) is the raw instantaneous pressure collected in step S102. Or vibration signal within the historical time window T in A continuous signal sequence assembled sequentially within the time frame; specifically, the historical observation vector P at the current time t is defined. window : in, The number of sampling points within the time window. The sampling interval; To meet the requirements of AVMD for continuous domain signal processing, the discrete vector is reconstructed into a continuous-time function f(τ) using the Shannon interpolation formula. The concatenation and reconstruction formula is as follows: in, As an interpolation kernel function, this process transforms discrete sensor readings into a smooth signal that includes transient pressure waveforms and nonlinear vibration modes; Then, the analytical form of the signal is constructed to obtain its spectral characteristics. The Hilbert transform is then used to convert the real signal into a complex analytic signal. The Hilbert transform formula is: in, Let pv denote the Hilbert transform operator, pv denote the Cauchy principal value integral, and u denote the principal value integral. k (τ) represents the k-th modal component; Based on this, the formula for constructing an analytical signal is: Among them, z k Let δ(t) be an analytic signal with a one-sided spectral characteristic (containing only positive frequencies), and let δ(t) be the Dirac distribution function. j is the imaginary unit; Furthermore, we construct a variational optimization problem, assuming the original signal consists of K frequencies revolving around the center frequency ω. k amplitude modulation and frequency modulation signal u k (t) are superimposed. In order to minimize the bandwidth of each mode, the spectrum of the analytic signal is modulated to the baseband by shifting the exponential term. The objective functional formula for AVMD optimization is: Among them, J AVMD ω represents the total bandwidth target to be minimized; k That is, the specific center frequency corresponding to this mode; This is a complex exponential twitch factor used to perform spectrum shifting operations, changing the center frequency... Shift to zero frequency; δ(t) is the time gradient operator used to calculate the smoothness of the demodulated baseband signal, i.e., the bandwidth; δ(t) is the Dirac distribution function. Constraints Ensure that the sum of the K decomposed modes can be used to reconstruct the original signal without loss; Finally, a signal reconstruction operation is performed, and a subset S of modes containing valid fault information is selected based on the kurtosis criterion. valid Noisy modes are removed, and the clean signal is reconstructed to build a denoised feature tensor X. feat The formula for denoising feature reconstruction is: in, Let be the denoised feature vector of node i at time t. The superposition of effective modes represents the denoised dynamic signal. These are static attribute features that do not require decomposition. Based on this, construct the denoising feature tensor X. feat This tensor is obtained by placing all nodes in the time window T. in Internally reconstructed feature vector Stacked along the time axis, the formula for constructing the denoising feature tensor is: in, Let X be the state snapshot matrix of the entire network after denoising at time t; the tensor X feat This will be used as the input data for the dual-stream network in step S3; Step S202: Based on the weighted adjacency matrix W constructed in step S101, construct the normalized graph Laplace matrix L, which reflects the topology and functional dependencies of the pipeline system, from the perspective of graph theory. Step S203: Based on the graph Laplace matrix L, a graph filter is designed using Chebyshev polynomial approximation theory to project the pipeline spatiotemporal signal into the graph frequency domain. The specific processing procedure is as follows: First, calculate the rescaled Laplacian matrix, mapping the eigenvalues ​​to the interval [-1, 1] to satisfy the domain requirement of the Chebyshev polynomial. The rescaling formula is: in, Let λ be the rescaled Laplace matrix. max The largest eigenvalue of the Laplace matrix L; Then, a Chebyshev polynomial sequence is defined using recursion to quickly compute feature extraction in , and the recursive formula is: in, It is a k-th order Chebyshev polynomial operator.

5. The pipeline risk dynamic perception and visualization method based on multidimensional spatiotemporal graph convolution as described in claim 4, characterized in that: In step S201, the objective functional formula J for AVMD optimization is solved. AVMD The steps for this constrained variational problem are as follows: First, by introducing the Lagrange multiplier λ(t) and the quadratic penalty factor α, we transform the problem into an unconstrained optimization problem. The augmented Lagrange function formula is as follows: Where L is the augmented Lagrangian function; the first term of the formula The second term is a quadratic penalty term, where α is a penalty parameter balancing data fidelity and bandwidth constraints, used to ensure reconstruction accuracy while constraining modal bandwidth; The residual fidelity term represents the energy difference between the original signal and the sum of the decomposed modes, i.e., the reconstruction error; the third term in the formula... For Lagrange multipliers, This indicates the inner product operation, used to strictly enforce the satisfaction of equality constraints; Then, the alternating direction multiplier method is used to alternately update the variables in the frequency domain until convergence is obtained to obtain K eigenmode functions u. k (t), the specific iteration steps are as follows: (1) The modal component update formula is: (2) The center frequency update formula is: (3) The Lagrange multiplier update formula is: in, The frequency domain representation of the variable is given, where n is the number of iterations. Update the step size for the multiplier.

6. The pipeline risk dynamic perception and visualization method based on multidimensional spatiotemporal graph convolution according to claim 4, characterized in that: The specific calculation steps for step S202 are as follows: First, calculate the degree matrix D, which is a diagonal matrix with diagonal elements D. ii The sum of all elements in the i-th row of the weighted adjacency matrix W reflects the connection density or importance of node i in the network. The calculation formula is: Among them, D ii W represents the diagonal element value in the i-th row and i-th column of the degree matrix. ij The element in the i-th row and j-th column of the weighted adjacency matrix represents the connection weight between node i and node j. Then, construct the degree matrix D, which is derived from D. ii The resulting diagonal matrix, i.e., the elements excluding the main diagonal, are D. ii Except for the elements at all other positions, all other elements are 0, satisfying D = diag(D 11 ,…,D NN ); Finally, the adjacency matrix is ​​symmetrically normalized using the degree matrix to construct the graph Laplacian matrix, which is calculated as follows: Among them, I N Let W be the N-order identity matrix, and W be the weighted adjacency matrix from step S101. Let L be the inverse square root of the degree matrix. The normalized Laplace matrix L has positive semi-definite properties, and its eigenvalues ​​λ i Distributed within the interval [0,2], this ensures the numerical stability and convergence of subsequent graph convolution operations.

7. The pipeline risk dynamic perception and visualization method based on multidimensional spatiotemporal graph convolution according to claim 4, characterized in that: Step S3 specifically includes: Step S301: Construct a two-stream feature extraction architecture, and convert the denoised feature tensor X output from step S2 into a single stream. feat Simultaneously inputting into two parallel processing branches: spatial graph convolutional stream and temporal attention stream, the aim is to capture the topological spatial dependency and temporal evolution of risk propagation in pipeline networks, respectively, and improve the feature extraction efficiency of the model by decoupling spatiotemporal features; Step S302: In the spatial graph convolutional flow, design a multi-level dynamic graph convolutional layer, which uses tensor X. feat Using spatial slices at each time step as input, convolution operations are performed by combining static physical topology with dynamic data association, and feature aggregation is performed using the Chebyshev polynomial generated in step S203. The specific calculation process is as follows: First, calculate the spatial self-attention matrix A. att This matrix reflects the similarity weights between the features of each node at the current time. The formula for calculating the spatial self-attention matrix is: Among them, A att Let X be the spatial self-attention matrix. (l) Let X be the input feature matrix of the l-th generation. When l=0, X(0) is the denoised feature vector X input in step S301. feat The spatial feature slice at the current time t, i.e., X feat [t, :, :]), W Q and W K Let d be a learnable linear mapping weight matrix. k The feature dimension scaling factor; Then, construct the fusion graph Laplacian operator L. mix The physical topology and data-driven association are weighted and fused together, and the fusion formula is as follows: Among them, L mix To fuse graph operators, γ is the static Chebyshev baseline matrix, and γ is the learnable balance coefficient. Finally, calculate the spatial stream output characteristics. The formula for calculating the spatial flow output characteristics is: in, For the output features of the (l+1)th layer, T k (⋅) represents the Chebyshev polynomial operator defined in step S203. Let K be the parameter matrix of the k-th convolution kernel. It is the ReLU activation function; Step S303: In the temporal attention flow, a causal gated temporal convolutional network and a multi-head temporal self-attention mechanism are designed to capture long-short-term temporal dependencies. The specific calculation process is as follows: First, local temporal features are extracted using dilated causal convolution. This step converts the tensor X... feat Unfolding along the node dimension, convolving the time series of each node, the output formula of the dilated causal convolution is: Among them, h t x is the local temporal feature vector output by the convolution. seq The input time series vector corresponds to X feat The feature sequence of a specific node (i.e., X) feat [:, i, :]), W represents a convolution operation with an expansion rate of d. f b f For filter parameters, W g b g For gating parameters, ⊙ represents the Hadamard product; Then, the convolution output h t Input a multi-head self-attention layer, and calculate global dependencies through multi-head self-attention. The multi-head self-attention calculation formula is as follows: Where MHSA(⋅) represents the output of the multi-head self-attention module, head i Q represents the output of the i-th attention head. i K i V i It is h t The query, key, and value matrix M obtained by linear transformation mask W is a causal mask matrix to ensure that future information is not leaked. O To output the projection matrix; Finally, output the temporal flow latent feature H. T The formula for calculating the latent features of time flow is: Among them, H T The temporal latent features representing the final output are fused with local convolutional features and global attention features, and their dimension is similar to the spatial flow output features H. S Maintaining consistency is essential for future integration; Step S304: Design a spatiotemporal fusion gating unit to adaptively fuse spatial flow characteristics H S With time flow characteristics H T The specific integration process is as follows: First, calculate the fusion gating coefficient z. The formula for calculating the fusion gating coefficient is: Where z is the fusion gating coefficient, || represents feature concatenation, and W z b z For gating network parameters; Then, the final risk latent feature tensor Z is generated using the gating coefficients. risk The formula for calculating the risk latent feature tensor is: Among them, Z risk The final output risk latent feature tensor enables adaptive adjustment of the spatiotemporal feature weights according to different risk types.

8. The pipeline risk dynamic perception and visualization method based on multidimensional spatiotemporal graph convolution according to claim 1, characterized in that: Step S4 specifically includes: Step S401: Introduce cloud model theory to address the fuzziness and randomness in risk assessment, and define three core numerical characteristics of the risk cloud model: expectation Ex, entropy En, and hyperentropy He; where Ex reflects the central location or average intensity of risk occurrence and represents the qualitative concept of risk level; En reflects the uncertainty and distribution range of risk performance and represents the fuzziness of the qualitative concept; and He reflects the uncertainty of entropy, that is, the discreteness or random fluctuation of risk generation. Step S402: Construct a deep cloud mapping network and convert the high-dimensional risk latent feature Z output in step S3 into a deep cloud mapping network. risk Mapping to the cloud feature space, calculate the forward cloud generator parameters of pipeline node i at the current time. The cloud parameter mapping formula is: in, Let be the risk latent feature vector of node i; , , These are the weight matrices for the expectation, entropy, and hyperentropy mapping networks, respectively. , , This is the corresponding bias vector; The sigmoid activation function normalizes the expected value to the interval [0, 1] to match the definition of risk probability. The exp(⋅) function ensures that the entropy value is positive. The function ensures that the superentropy is positive and smoothly differentiable. Through this network, abstract deep learning features are transformed into digital features of cloud models with physical meaning. Step S403: Perform Monte Carlo simulation based on the forward normal cloud algorithm to generate M risk cloud droplets. To quantify risk distribution, the specific generation steps are as follows: First, according to entropy En i and superentropy He i Generate with En i For expectations, He i Normal random entropy with standard deviation The formula for generating normal random entropy is: in, This indicates the generation of normally distributed random numbers with mean μ and variance σ. This reflects the uncertainty bandwidth of risk in a single sampling. Then, according to the expectation Ex i and the generated random entropy Generate a quantitative risk value x m The formula for generating the quantitative risk value is: Where, x m Let m be the specific value of the m-th cloud droplet in the number domain space, representing a specific risk prediction value; Finally, calculate the quantitative value x. m Certainty μ of the concept of "high risk" m The formula for calculating the degree of certainty is: Where, μ m x represents m The probability intensity belonging to the current risk level, with a value range of [0, 1], is repeated M times to obtain the cloud droplet set {(x1, μ1), (x2, μ2), … , (x M , μ M This set constitutes a complete cloud map description of the risk; Step S404: Calculate the comprehensive risk assessment index R i and assess confidence level C i By utilizing the centroid method of cloud droplet swarms, the uncertain distribution of cloud droplets is transformed into a single decision value. The risk assessment index is calculated using the following formula: Among them, R i The comprehensive risk assessment index for node i is a weighted index that takes into account the magnitude and certainty of each cloud droplet. The formula for calculating the confidence level is: Among them, C i This represents the average degree of certainty of the evaluation results based on the cloud droplet swarm generated from the current M Monte Carlo simulations; if R i Higher and C i A lower R value indicates that while cloud droplets are biased towards high risk, their certainty is dispersed, suggesting high uncertainty in the risk and requiring manual verification; if R... i Higher and C i A high level indicates that the cloud droplets are concentrated and clearly defined, thus indicating a confirmed high risk and automatically triggering an alarm.

9. The pipeline risk dynamic perception and visualization method based on multidimensional spatiotemporal graph convolution according to claim 8, characterized in that: Step S5 specifically includes: Step S501: Based on pipeline geographic information system data and building information model, construct a high-precision 3D geometric model of the pipeline network, and use an octree structure for spatial index optimization to support real-time rendering in large-scale pipeline network scenarios. The specific construction steps are as follows: First, multi-source heterogeneous data is analyzed to obtain the pipeline's geometric and attribute information; specifically, the pipeline geographic information system data refers to the geospatial coordinate sequence of the pipeline centerline. The building information model data specifically refers to: the cross-sectional properties of the pipes, the material texture mapping coordinates, and the detailed three-dimensional mesh model of the ancillary facilities. Then, a unified spatial coordinate system is established, converting the geographic coordinates of the GIS data into Cartesian rectangular coordinates used by the rendering engine. The coordinate transformation formula is as follows: Among them, [x k , y k , z k ] T M represents the transformed local coordinates. proj T is the Mercator projection matrix used to project spherical coordinates onto a plane. offset This is the scene center offset vector, used to solve the problem of precision loss in large coordinate floating-point numbers; Furthermore, a 3D mesh of the pipeline is generated. For each pipeline segment, the mesh is generated based on its centerline coordinates (x, y, y). k , y k , z k ) and (x k+1 ,y k+1 ,z k+1 And the pipe diameter (diam), using parametric modeling to generate a set of cylindrical triangular facets, the vertex generation formula is: Where v(θ,h) is the coordinate of the vertex of the cylinder surface, c(h) is the point at parameter h on the center line, n and b are the normal vector and binormal vector at that point, and θ∈[0,2π) is the inscribed angle; Finally, an octree spatial index is constructed, and the generated millions of triangular facet primitives are recursively divided into eight child nodes to establish hierarchical bounding boxes. The virtual camera is then initialized, and the formula for calculating the axis-aligned bounding box of each node is as follows: Among them, v set The set of all vertices contained in this node; Meanwhile, a virtual roaming camera is constructed in the 3D scene, and its intrinsic parameter matrix and initial extrinsic parameter matrix are defined. This camera serves as the controlled object for intelligent focusing linkage in the subsequent step S6, and is used to generate the final rendering viewport. Step S502: Design a four-dimensional risk field color mapping function Establish a risk index R i The non-linear mapping relationship between node color and transparency makes high-risk areas more visually prominent. The node color mapping formula is as follows: Among them, Color i Let i be the RGBA color vector containing red, green, blue, and alpha components. safe Color is the preset safety status color. danger R is the preset color for dangerous situations. i r is the comprehensive risk assessment index of node i calculated in step S404. low r high These represent the low and high thresholds for risk visualization, respectively. SmoothStep(⋅) is a smoothing step function used to adjust R... i In [r low , r high The smooth interpolation coefficients are mapped to [0, 1] within the interval. Lerp(⋅) is a linear interpolation function. The risk value is converted into a visual color signal by mixing between safe and dangerous colors based on the interpolation coefficients. Step S503: The force-directed algorithm is used to simulate the dynamic diffusion effect of risk energy in the pipeline network. For the identified high-risk nodes, a dynamic particle flow is generated. The motion of the particle flow is influenced by the risk gradient force, topological constraint force, and random perturbation force. The dynamic equation of motion of the particles is: Where p is the particle's position vector (x, y, z) in three-dimensional space. Let be the particle's acceleration vector. Let F be the velocity vector of the particle. risk The risk gradient force is the main driving force for particle motion, k d F is the motion damping coefficient, used to simulate fluid viscous drag and prevent particle velocity from increasing indefinitely. topology As a topological constraint, particles are restricted to tangential movement along the pipe axis to prevent them from penetrating the pipe wall and overflowing. F random To simulate the random forces of Brownian motion in turbulence, and to increase the naturalness and randomness of the visual effect; Risk gradient force F risk The calculation formula is: Where ∇R(p) is the gradient of the risk field at particle position p, N(p) is the set of neighborhood nodes of the particle's current position, and R p R is the risk value obtained by interpolating the particle's location. j Let be the risk value of the neighboring node j; this force simulates the natural diffusion process of risky substances from high-risk areas to low-risk areas with the risk potential energy gradient; Step S504: In the WebGL rendering pipeline, write a custom fragment shader to perform physically based halo rendering on the pipeline walls, visually representing the uncertainty boundary of the risk. The halo intensity calculation formula is: Among them, I glow The final halo intensity of the current fragment, α is the basic luminous intensity coefficient that controls the overall brightness, En i The node risk entropy is calculated in step S402; n is the normal vector of the pipe wall surface, v is the line-of-sight vector from the fragment to the camera, and max(0, n⋅v) p Here, p represents the Fresnel effect term, and p is the edge decay exponent. The item introduces a "breathing" effect over a time dimension; ω represents the breathing frequency, set to be correlated with the risk index R. i Proportional This represents the phase offset.

10. The pipeline risk dynamic perception and visualization method based on multidimensional spatiotemporal graph convolution according to claim 1, characterized in that: Step S6 specifically includes: Step S601: Establish a risk classification and determination mechanism, and use the comprehensive risk assessment index calculated in step S404. The risk level of the current node is determined by comparing it with the preset multi-level risk threshold vector r; the risk level determination logic is as follows: Among them, Level i Let r represent the risk level of node i; r = [r1, r2, r3] is an increasing risk threshold vector, corresponding to the trigger limits for early warning, alarm, and emergency shutdown, respectively. Step S602: Construct an intelligent focusing linkage mechanism based on a visualization engine. When the risk level of any node i reaches Alarm or Emergency, the viewpoint control algorithm of the digital twin interface is automatically triggered, driving the virtual roaming camera in the 3D scene constructed in step S501 to move smoothly from the current viewpoint to the optimal observation viewpoint of the abnormal node. The formula for calculating the smooth motion trajectory of the camera is: Among them, Cam pos (t) represents the spatial coordinates and attitude quaternion of the camera at time t; Cam current The current roaming viewpoint; The optimal viewing point for anomaly node i; is a spherical linear interpolation function used to achieve a smooth, seamless camera movement effect; Smooth(t) is a easing function to ensure that the camera movement conforms to visual inertia. Step S603: Generate an edge collaborative control instruction package and execute closed-loop feedback. Based on the determined risk level and type, automatically generate the corresponding control strategy. If the level is Warning, generate adjustment instructions. Among them, u control η is the control quantity issued, η is the adjustment step size, and ∇J opt (u) represents the optimization gradient based on the pipeline hydraulic balance model; If the emergency level is set to Emergency, an emergency shutdown signal is issued directly. The command is sent to the field edge computing gateway via the low-latency MQTT protocol for execution, and the device status and sensor feedback data after execution are used to construct the original spatiotemporal tensor X for the next moment. raw The input data is used to re-enter step S1, realizing real-time closed-loop verification from perception to control and back to perception.

Images