Three-dimensional reconstruction method and system for gas leakage based on thermal-optical fusion

Abstract

The application provides a gas leakage three-dimensional reconstruction method and system based on thermo-optical fusion, relates to the technical field of three-dimensional reconstruction, and comprises the following steps: multi-scale features are extracted from a thermal radiation image sequence and a visible light image sequence, and are fused based on mutual information; a space-time causal graph is constructed to perform directed message passing to obtain enhanced node features; the enhanced node features are mapped to a continuous space-time domain and are transformed by a Fourier neural operator to obtain a feature field; the feature field is used as a condition to guide iterative denoising; and finally, a three-dimensional space distribution model of gas leakage is generated. The application can accurately reconstruct the three-dimensional space distribution of gas leakage, and improve the accuracy and reliability of leakage monitoring.

Description

Technical Field

[0001] This invention relates to three-dimensional reconstruction technology, and more particularly to a method and system for three-dimensional reconstruction of gas leaks based on thermo-optical fusion. Background Technology

[0002] Gas leak detection and visualization is an important research direction in the fields of industrial safety and environmental protection. Traditional gas leak detection methods mainly include thermal infrared imaging and visible light imaging. Thermal infrared imaging identifies gas leak areas by detecting the temperature difference between the gas and the environment, enabling non-contact, long-distance detection; while visible light imaging provides high-resolution environmental context by capturing the structural and textural information of the scene. However, existing gas leak detection technologies have the following drawbacks and limitations:

[0003] Most existing methods rely on a single imaging modality, resulting in incomplete information acquisition. While thermal infrared imaging can effectively detect the presence of gases, it lacks environmental context information; visible light imaging can provide rich environmental information, but its ability to identify the gases themselves is limited. Both modalities have information blind spots, making it difficult to achieve accurate leak source location and concentration distribution estimation.

[0004] Existing technologies generally lack sufficient consideration of the spatiotemporal dynamics of gas leaks. Gas leaks are a continuously changing dynamic process, and their diffusion behavior is affected by a variety of environmental factors. However, most methods only focus on static images or simple inter-frame differences, which cannot effectively capture the spatiotemporal variation patterns of gas leaks, resulting in inaccurate reconstruction results.

[0005] Existing 3D reconstruction methods mostly employ discrete representations, making it difficult to depict the continuous characteristics of gas distribution. The distribution of gas in space is essentially a continuous concentration field, while discrete representations lead to loss of detail and unnatural boundaries. Especially with limited resolution, it is impossible to accurately reconstruct the true 3D spatial distribution of gas leaks, affecting subsequent analysis and decision-making. Summary of the Invention

[0006] The present invention provides a method and system for three-dimensional reconstruction of gas leaks based on thermo-optical fusion, which can solve the problems in the prior art.

[0007] A first aspect of the present invention provides a method for three-dimensional reconstruction of gas leaks based on thermo-optical fusion, comprising:

[0008] Multi-scale feature maps are extracted from the thermal radiation image sequence and visible light image sequence of the acquired target area, respectively. The local mutual information distribution of thermal radiation features and visible light features is calculated at each scale, weighted and fused, and then stitched together to generate a multi-scale feature representation.

[0009] The multi-scale feature representations of consecutive frames are spatially gridded to construct a spatiotemporal candidate graph. The causal strength between nodes is calculated. The spatiotemporal causal graph is obtained by pruning according to the causal strength. Causal-guided directed message passing is performed on the spatiotemporal causal graph to obtain enhanced node features.

[0010] The enhanced node features are mapped to the continuous spatiotemporal domain, and after feature transformation in the frequency domain by Fourier neural operators, an inverse transformation is performed to obtain the continuous spatiotemporal feature field.

[0011] Using the continuous spatiotemporal feature field as a condition, the initial noise field is iteratively denoised. Each iteration predicts and removes noise components until convergence is achieved to obtain the depth field distribution map and the concentration field distribution map.

[0012] Based on the depth field distribution map and the concentration field distribution map, a three-dimensional spatial distribution model of the gas leak is generated.

[0013] The steps for generating a multi-scale feature representation by extracting multi-scale feature maps from the acquired thermal radiation image sequence and visible light image sequence of the target area, calculating the local mutual information distribution of thermal radiation features and visible light features at each scale, weighting and fusing them, and then stitching them together include:

[0014] The thermal radiation image sequence and the visible light image sequence are downsampled by a multi-layer convolutional network. After each downsampling, the thermal radiation feature map and the visible light feature map at the current scale are extracted to obtain multiple scale levels from high resolution to low resolution.

[0015] At each scale level, for the corresponding spatial locations of the thermal radiation feature map and the visible light feature map, the KL divergence between the joint probability distribution and the marginal probability distribution of the feature vectors is calculated as the local mutual information value, and the local mutual information value is normalized and used as the fusion weight.

[0016] The thermal radiation feature vector and visible light feature vector at the same spatial location are weighted and summed using the fusion weights to obtain weighted fusion feature maps at each scale. The weighted fusion feature maps at different scales are upsampled to the same spatial size and then stitched together in the channel dimension to generate a multi-scale feature representation.

[0017] The steps of constructing a spatiotemporal candidate graph by spatially meshing the multi-scale feature representations of multiple consecutive frames, calculating the causal strength between nodes, pruning according to the causal strength to obtain a spatiotemporal causal graph, and performing causally guided directed message passing on the spatiotemporal causal graph to obtain enhanced node features include:

[0018] The multi-scale feature representations of multiple consecutive frames are spatially gridded, each grid cell is mapped to a node on the time axis and candidate edges are established to construct a spatiotemporal candidate graph;

[0019] Historical state sequences are extracted from the nodes in the spatiotemporal candidate graph. The baseline error predicted using only the historical state of the target node, the source error predicted by adding the historical state of the source node, and the conditional error predicted by adding the historical state of other nodes are calculated respectively. When the difference between the baseline error and the source error exceeds the error improvement threshold and the difference between the source error and the conditional error is lower than the conditional influence threshold, the difference between the baseline error and the source error is normalized and used as the hysteresis causality strength.

[0020] The conditional mutual information between node pairs is calculated using the kernel density estimation method, and the instantaneous causal strength is obtained by combining spatial distance priors and adaptive threshold filtering.

[0021] The causal intensity and the instantaneous causal intensity are nonlinearly combined to obtain a causal intensity matrix. The pruning threshold is determined based on the statistics of the causal intensity matrix, and low-intensity edges are removed to obtain a spatiotemporal causal graph.

[0022] On the spatiotemporal causal graph, historical source node messages and simultaneous source node messages are aggregated respectively, and after being fused through a gating mechanism, they are combined with the original node features to obtain enhanced node features.

[0023] The steps of mapping the enhanced node features to a continuous spatiotemporal domain, performing feature transformation in the frequency domain using Fourier neural operators, and then performing inverse transformation to obtain a continuous spatiotemporal feature field include:

[0024] Coordinate embedding is performed on the enhanced node features, and the spatial and temporal coordinates of the discrete grid nodes are encoded into positional representations, which are then concatenated with the enhanced node features as input features.

[0025] The input features are transformed in the frequency domain through multiple Fourier layers. Each Fourier layer performs a fast Fourier transform to map the input to the frequency domain space. In the frequency domain space, a learnable frequency domain convolution kernel is used to perform a weighted transformation on different frequency components. Then, an inverse fast Fourier transform is used to map the input back to the spatial domain. The frequency domain convolution kernel captures global spatiotemporal dependencies.

[0026] In the spatial domain, nonlinear activation and residual connection are performed on the frequency domain transformed features to obtain frequency domain enhanced features;

[0027] The dominant causal path is extracted from the spatiotemporal causal graph and represented as a spatiotemporal propagation direction field; a direction-adaptive spatiotemporal kernel function is constructed, the anisotropy of which is determined by the propagation direction field.

[0028] For any query coordinate in the continuous spatiotemporal domain, the contribution weight of each grid node to the query coordinate is calculated based on its spatiotemporal distance from the discrete grid nodes and the spatiotemporal kernel function; the frequency domain enhancement features are weighted and combined according to the contribution weights to obtain the continuous spatiotemporal feature field.

[0029] The steps for calculating frequency domain enhancement features include:

[0030] The input features are subjected to Fast Fourier Transform along the spatial and temporal dimensions to obtain a spectral representation. Based on the energy distribution of the frequency components in the spectral representation, the spectrum is divided into different frequency regions, and feature transformation is performed on different frequency regions using frequency domain convolution kernels with different channel capacities.

[0031] The spectrum representations after transformation of different frequency regions are spliced ​​in the frequency domain space, and different frequency components are adaptively weighted through a frequency domain attention mechanism, which calculates the attention weights based on the spectral energy distribution.

[0032] Perform an inverse fast Fourier transform on the weighted spectral representation to map it back to the spatial domain, and obtain the frequency domain transformed features;

[0033] The frequency-domain transformed features are nonlinearly activated; the activated features are then residually connected with the input features to obtain frequency-domain enhanced features.

[0034] The steps of iteratively denoising the initial noise field using the continuous spatiotemporal feature field as a condition, predicting and removing noise components in each iteration until convergence is achieved to obtain the depth field distribution map and the concentration field distribution map, include:

[0035] An initial noise field is randomly generated, following a zero-mean, unit-variance distribution. The spatiotemporal dimension of the initial noise field is correlated with the target depth field distribution and concentration field distribution. Figure 1 To;

[0036] The initial noise field is iteratively denoised. In each iteration, the current noise field and the continuous spatiotemporal feature field are interacted through a cross-attention mechanism to obtain a conditionally enhanced noise field representation.

[0037] The noise component is predicted by a denoising prediction network based on the conditionally enhanced noise field representation and the current iteration step index; temporal causal constraints are extracted from the spatiotemporal causal graph and encoded as a causal mask matrix to constrain the noise component; the noise component constrained by the causal mask is subtracted from the current noise field to obtain the denoised noise field; the noise removal step size is dynamically adjusted according to the current iteration step index, and the denoised noise field is used as the current noise field for the next iteration;

[0038] The iteration is terminated when the rate of change of the noise field between adjacent iterations is lower than the preset rate of change threshold or the maximum number of iterations is reached; the final denoised noise field is decoupled into a depth field distribution map and a concentration field distribution map.

[0039] The steps for generating a three-dimensional spatial distribution model of gas leakage based on the depth field distribution map and the concentration field distribution map include:

[0040] Spatial geometry reconstruction is performed on the depth field distribution map to convert depth values ​​into three-dimensional spatial coordinates, establishing a three-dimensional geometric structure of the leakage scene. Concentration values ​​in the concentration field distribution map are mapped to corresponding spatial positions in the three-dimensional geometric structure, establishing a correspondence between concentration values ​​and three-dimensional spatial coordinates through voxelization. Spatial interpolation is performed on the voxelization representation, calculating interpolated concentration values ​​at continuous spatial positions between voxel grids to obtain a continuous three-dimensional concentration distribution. Concentration isosurfaces are set based on the three-dimensional concentration distribution, and three-dimensional surfaces corresponding to different concentration thresholds are extracted. These three-dimensional surfaces represent the spatial boundaries of gas diffusion.

[0041] The three-dimensional geometric structure, the three-dimensional concentration distribution, and the three-dimensional surface are fused to generate a three-dimensional spatial distribution model of gas leakage.

[0042] A second aspect of the present invention provides a three-dimensional reconstruction system for gas leaks based on thermo-optical fusion, comprising:

[0043] The first unit is used to extract multi-scale feature maps from the thermal radiation image sequence and visible light image sequence of the acquired target area, respectively. The local mutual information distribution of thermal radiation features and visible light features is calculated at each scale, weighted and fused, and then spliced ​​to generate a multi-scale feature representation.

[0044] The second unit is used to construct a spatiotemporal candidate graph by spatial grid partitioning of the multi-scale feature representations of multiple consecutive frames, calculate the causal strength between nodes, prune according to the causal strength to obtain a spatiotemporal causal graph, and perform causal-guided directed message passing on the spatiotemporal causal graph to obtain enhanced node features.

[0045] The third unit is used to map the enhanced node features to the continuous spatiotemporal domain, and then perform feature transformation and inverse transformation in the frequency domain using Fourier neural operators to obtain the continuous spatiotemporal feature field.

[0046] The fourth unit is used to iteratively denoise the initial noise field based on the continuous spatiotemporal feature field. In each iteration, noise components are predicted and removed until convergence is obtained to obtain the depth field distribution map and the concentration field distribution map.

[0047] The fifth unit is used to generate a three-dimensional spatial distribution model of gas leakage based on the depth field distribution map and the concentration field distribution map.

[0048] A third aspect of the embodiments of the present invention,

[0049] An electronic device is provided, comprising:

[0050] processor;

[0051] Memory used to store processor-executable instructions;

[0052] The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0053] Fourth aspect of the embodiments of the present invention,

[0054] A computer-readable storage medium is provided, having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0055] This invention establishes an efficient three-dimensional reconstruction technology for gas leaks by fusing multi-scale features from thermal radiation and visible light image sequences. It can accurately characterize the distribution and concentration changes of gases in space, achieving precise mapping from two-dimensional images to three-dimensional space, and providing a more intuitive visualization for gas leak monitoring.

[0056] This invention introduces a spatiotemporal causal graph and a causal-guided message passing mechanism to effectively capture the spatiotemporal dynamic evolution of gas leakage, significantly improve the modeling ability of gas diffusion behavior in complex environments, solve the problem of insufficient reconstruction accuracy of traditional methods under non-uniform backgrounds and illumination changes, and enhance the environmental adaptability of the system.

[0057] This invention employs a continuous feature field representation based on Fourier neural operators and a conditional denoising iterative mechanism, which overcomes the resolution limitations caused by discrete sampling and achieves high-precision synchronous reconstruction of depth and concentration fields. Attached Figure Description

[0058] Figure 1 This is a schematic flowchart of the gas leakage three-dimensional reconstruction method based on thermo-optical fusion according to an embodiment of the present invention;

[0059] Figure 2 Flowchart for constructing and computing enhanced node features for spatiotemporal causal graphs. Detailed Implementation

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

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

[0062] Figure 1 This is a schematic flowchart of the gas leakage three-dimensional reconstruction method based on thermo-optical fusion according to an embodiment of the present invention, as shown below. Figure 1 As shown, the method includes:

[0063] Multi-scale feature maps are extracted from the thermal radiation image sequence and visible light image sequence of the acquired target area, respectively. The local mutual information distribution of thermal radiation features and visible light features is calculated at each scale, weighted and fused, and then stitched together to generate a multi-scale feature representation.

[0064] The multi-scale feature representations of consecutive frames are spatially gridded to construct a spatiotemporal candidate graph. The causal strength between nodes is calculated. The spatiotemporal causal graph is obtained by pruning according to the causal strength. Causal-guided directed message passing is performed on the spatiotemporal causal graph to obtain enhanced node features.

[0065] The enhanced node features are mapped to the continuous spatiotemporal domain, and after feature transformation in the frequency domain by Fourier neural operators, an inverse transformation is performed to obtain the continuous spatiotemporal feature field.

[0066] Using the continuous spatiotemporal feature field as a condition, the initial noise field is iteratively denoised. Each iteration predicts and removes noise components until convergence is achieved to obtain the depth field distribution map and the concentration field distribution map.

[0067] Based on the depth field distribution map and the concentration field distribution map, a three-dimensional spatial distribution model of the gas leak is generated.

[0068] In one optional implementation, the steps of extracting multi-scale feature maps from the acquired thermal radiation image sequence and visible light image sequence of the target area, calculating and weighting the local mutual information distribution of thermal radiation features and visible light features at each scale, and then stitching them together to generate a multi-scale feature representation include:

[0069] The thermal radiation image sequence and the visible light image sequence are downsampled by a multi-layer convolutional network. After each downsampling, the thermal radiation feature map and the visible light feature map at the current scale are extracted to obtain multiple scale levels from high resolution to low resolution.

[0070] At each scale level, for the corresponding spatial locations of the thermal radiation feature map and the visible light feature map, the KL divergence between the joint probability distribution and the marginal probability distribution of the feature vectors is calculated as the local mutual information value, and the local mutual information value is normalized and used as the fusion weight.

[0071] The thermal radiation feature vector and visible light feature vector at the same spatial location are weighted and summed using the fusion weights to obtain weighted fusion feature maps at each scale. The weighted fusion feature maps at different scales are upsampled to the same spatial size and then stitched together in the channel dimension to generate a multi-scale feature representation.

[0072] For example, the thermal radiation image sequence and the visible light image sequence are acquired through a synchronous acquisition device, which includes a thermal infrared camera and a visible light camera, both of which are time-stamped and spatially registered. The thermal radiation image sequence has a resolution of 640×512 pixels, a frame rate of 30 frames per second, and a wavelength range of 8 micrometers to 14 micrometers; the visible light image sequence has a resolution of 1920×1080 pixels and a synchronized frame rate of 30 frames per second. The acquired image sequences are arranged in chronological order, with each sequence containing at least 10 consecutive images to ensure the integrity of the temporal information.

[0073] The multi-layer convolutional network employs an encoder structure, comprising four downsampling levels. The first level receives the original resolution input, extracts features using a 3×3 convolutional kernel with a stride of 1, and uses zero-padding to maintain spatial dimensions, resulting in 64 output channels. Downsampling is achieved through 2×2 max pooling with a stride of 2, halving the spatial dimensions. The second level takes the pooled output from the first level as input, also using a 3×3 convolutional kernel, increasing the number of output channels to 128. The third and fourth levels have 256 and 512 output channels, respectively. Thermal radiation image sequences and visible light image sequences are processed through their respective independent convolutional networks. The hierarchical structure and parameter configuration of the two networks are completely symmetrical, but the weight parameters are trained independently. The feature maps extracted after each downsampling step are the thermal radiation feature map and visible light feature map at the current scale. The spatial dimensions corresponding to the four scale levels are the original size, half, quarter, and eighth of the original size, respectively.

[0074] At each scale level, the corresponding spatial locations of the thermal radiation feature map and the visible light feature map refer to the pixel locations with the same coordinate indices after spatial registration. For a given spatial location, the thermal radiation feature vector and the visible light feature vector are extracted, and the dimension of the feature vector is equal to the number of channels at that scale level. When calculating the joint probability distribution, all combinations of the two feature vectors are used as the joint event space, and the probability is estimated by counting the frequency of feature values ​​within a local neighborhood window around the current location. The size of the local neighborhood window is set to 7×7 pixels, and the feature vectors at all locations within the window participate in the probability statistics. The joint probability distribution represents the probability that thermal radiation feature values ​​and visible light feature values ​​occur simultaneously, while the marginal probability distributions represent the probabilities that thermal radiation feature values ​​and visible light feature values ​​occur individually, respectively. The KL divergence is used to measure the difference between the product of the joint probability distribution and the marginal probability distribution; the larger the difference, the stronger the correlation between the two features at that location. In the specific calculation, the feature values ​​are first discretized, mapping continuous feature values ​​to 256 discrete intervals, with the width of each interval evenly divided according to the range of feature values. The joint probability distribution is estimated by counting the number of times each characteristic value falls into each discrete interval combination within a statistical window, and then dividing by the total number of samples within the window. The marginal probability distribution is estimated similarly, by counting the number of times each thermal radiation characteristic value and visible light characteristic value falls into each interval. The KL divergence is calculated as the expected value of the ratio of the logarithm of the product of the joint probability distribution and the marginal probability distribution, accumulated over all discrete interval combinations. The calculated KL divergence value is the local mutual information value at that spatial location.

[0075] Local mutual information values ​​form a mutual information distribution map with the same spatial size as the feature map across the entire feature map. To convert these mutual information values ​​into fusion weights, the mutual information distribution map is normalized. The normalization method uses the softmax function, calculating the exponential function value of the mutual information value at each location and then dividing it by the sum of the exponential function values ​​at all locations. The normalized fusion weights range from 0 to 1, and the sum of the weights at all locations is 1. A higher fusion weight indicates a stronger correlation between the thermal radiation characteristics and visible light characteristics at that location, and should be assigned a greater weight during the fusion process.

[0076] When using fusion weights to weight and sum the thermal radiation feature vector and visible light feature vector at the same spatial location, for each channel component of the feature vector, the sum of the thermal radiation component multiplied by the fusion weight and the visible light component multiplied by the residual weight is calculated separately. The residual weight is defined as 1 minus the fusion weight, ensuring that the sum of the contribution ratios of the two features is 1. The feature vector obtained after weighted summation has the same dimension as the input feature vector; this vector is the weighted fusion feature for that spatial location. Repeating the above weighted summation operation for all spatial locations on the feature map yields the weighted fusion feature map at that scale level.

[0077] Four weighted fused feature maps are generated at four scale levels, with spatial dimensions corresponding to the original size, 1 / 2, 1 / 4, and 1 / 8 of the original size, respectively, from largest to smallest. To achieve the concatenation of feature maps at different scales, all feature maps need to be upsampled to the same spatial size. The target upsampling size is chosen to be the original size, i.e., the largest spatial size. The upsampling method uses bilinear interpolation. For the feature map to be enlarged, the new pixel value is calculated at the new pixel location using the weighted average of its four nearest neighbor pixels. The weights are determined by the inverse ratio of the distance between the new pixel location and its four neighbors. The four upsampled feature maps have completely identical spatial dimensions, with channel dimensions of 64, 128, 256, and 512, respectively. Concatenation along the channel axis refers to stacking the four feature maps along the channel axis. The spatial size of the concatenated feature map remains the original size, and the number of channels is the sum of the number of channels in the four feature maps, 960. This concatenated feature map is the final multi-scale feature representation, containing multi-level semantic information from fine-grained to coarse-grained.

[0078] This invention achieves adaptive fusion of thermal radiation and visible light features by calculating the local mutual information distribution, avoiding the information loss problem caused by fixed weight fusion. The construction of multi-scale feature representation retains hierarchical information from details to the global.

[0079] In one optional implementation, the steps of constructing a spatiotemporal candidate graph by spatially meshing the multi-scale feature representations of multiple consecutive frames, calculating the causal strength between nodes, pruning according to the causal strength to obtain a spatiotemporal causal graph, and performing causally guided directed message passing on the spatiotemporal causal graph to obtain enhanced node features include:

[0080] The multi-scale feature representations of multiple consecutive frames are spatially gridded, each grid cell is mapped to a node on the time axis and candidate edges are established to construct a spatiotemporal candidate graph;

[0081] Historical state sequences are extracted from the nodes in the spatiotemporal candidate graph. The baseline error predicted using only the historical state of the target node, the source error predicted by adding the historical state of the source node, and the conditional error predicted by adding the historical state of other nodes are calculated respectively. When the difference between the baseline error and the source error exceeds the error improvement threshold and the difference between the source error and the conditional error is lower than the conditional influence threshold, the difference between the baseline error and the source error is normalized and used as the hysteresis causality strength.

[0082] The conditional mutual information between node pairs is calculated using the kernel density estimation method, and the instantaneous causal strength is obtained by combining spatial distance priors and adaptive threshold filtering.

[0083] The causal intensity and the instantaneous causal intensity are nonlinearly combined to obtain a causal intensity matrix. The pruning threshold is determined based on the statistics of the causal intensity matrix, and low-intensity edges are removed to obtain a spatiotemporal causal graph.

[0084] On the spatiotemporal causal graph, historical source node messages and simultaneous source node messages are aggregated respectively, and after being fused through a gating mechanism, they are combined with the original node features to obtain enhanced node features.

[0085] Combination Figure 2 The flowchart illustrating the construction and computation of spatiotemporal causal graphs to enhance node features is provided. Multi-scale feature representations across multiple consecutive frames are arranged chronologically, with each frame's feature representation having a spatial size of 640×512 pixels and 960 channels. A uniform block partitioning strategy is used, dividing the spatial dimension of each frame's feature representation into 32×32 grid cells, each with a size of 20×16 pixels. After grid partitioning, the feature vectors of all pixels within each grid cell are aggregated into a single feature vector through average pooling, maintaining a dimension of 960. For a feature representation sequence containing 15 frames, 1024 grid cells are generated per frame. On the time axis, each grid cell maps to a node, and the spatiotemporal coordinates of a node are determined by its frame's time index and the grid cell's spatial location index. Candidate edges are established following the spatiotemporal adjacency principle. For a node with time index t, its candidate edges connect to all nodes with time indices between t-1 and t+1 and whose spatial location is within a 3×3 neighborhood. The direction of the edges is determined by temporal order, pointing from nodes with smaller time indices to nodes with larger time indices. The constructed spatiotemporal candidate graph contains 15,360 nodes and approximately 138,240 candidate edges.

[0086] Historical state sequence extraction is based on node time indices. For a target node with time index t, its historical state sequence contains node feature vectors for all corresponding spatial locations from time index t-5 to t-1, with a historical sequence length of 5. The baseline error is calculated using an autoregressive prediction model, which predicts the target node's feature vector at time index t using only the historical state sequence. The prediction model employs a single-layer Long Short-Term Memory (LSTM) network, containing one LSTM unit layer and a fully connected output layer with 512 hidden units. During input processing, the feature vectors from the five historical time points are concatenated in chronological order to form the input sequence. Each feature vector has a dimension of 960, resulting in a 5×960 dimension input sequence. The LSTM network expands step-by-step, receiving one historical feature vector at each time step. The hidden state output of the last time step is mapped to a 960-dimensional predicted feature vector via a fully connected layer. The weight matrix of the fully connected layer has a dimension of 512×960, the bias vector has a dimension of 960, and linear activation is used. The baseline error is defined as the mean squared error between the predicted and true feature vectors. It is calculated by subtracting the corresponding elements of the predicted and true vectors, squared, and then averaging over all 960 elements. The source error is calculated by adding the historical state sequence of the source nodes as an additional input to the baseline error. This historical state sequence also contains feature vectors with time indices from t-5 to t-1. Candidate edges connect the source and target nodes. The source error is calculated using an extended autoregressive prediction model. This extended model uses the same single-layer long short-term memory (LSTM) network structure as the baseline model, but increases the number of hidden units to 768 to accommodate more input information. The input to the extended model is a concatenation of the historical sequences of the target and source nodes. The concatenation method involves chaining the target and source node feature vectors together at each time step to form a joint feature vector with a dimension of 1920. The input sequence for the five time steps has a dimension of 5×1920. After processing the joint input sequence by the LSM network, the hidden states at the final time step are passed through a fully connected layer to output the predicted feature vector. The weight matrix of the fully connected layer has a dimension of 768×960. The source error is also calculated using the mean squared error between the predicted and true vectors. The conditional error calculation further incorporates the historical state sequences of all neighboring nodes other than the source node. These other nodes refer to all nodes connected to the target node via edges in the candidate graph but not including the current source node. Assuming the target node has k other neighboring nodes, the input to the conditional error prediction model is a concatenation of the historical sequences of the target node, the source node, and the k other nodes. The input dimension at each time step is 960×(2+k). The conditional error prediction model uses the same Long Short-Term Memory (LSTM) network architecture, with the number of hidden units dynamically adjusted according to the input dimension, set to 512+128×k to ensure the model capacity increases with input complexity. The weight matrix of the fully connected output layer has a dimension of (512+128×k)×960.The three prediction models were trained using the same optimizer and loss function: the Adam optimizer with a learning rate of 0.001 and mean squared error loss. The training dataset was sampled from historical frames of the spatiotemporal candidate graph. Each node's training samples included its historical sequence and the true current-time feature vector. During training, the batch size was set to 128, the training epochs were 50, and the validation set comprised 10% of the total data for an early stopping strategy. Training was stopped when the validation set loss did not decrease for five consecutive epochs. The error improvement threshold was set to 0.05, and the conditional influence threshold was set to 0.02. A source node was considered to have a significant lagged causal influence on the target node when the difference between the baseline error and the source error was greater than 0.05 and the difference between the source error and the conditional error was less than 0.02. The lagged causal strength was normalized by subtracting the minimum difference of all candidate edge differences from the difference between the baseline error and the source error, then dividing by the range of the differences. The normalized lagged causal strength ranged from 0 to 1.

[0087] Kernel density estimation is used to calculate the conditional mutual information between node pairs, which are two nodes connected by a candidate edge at the same time index t. Conditional mutual information measures the statistical dependence of the source node's feature vector on the target node's feature vector given the feature vectors of other neighboring nodes. Kernel density estimation uses a Gaussian kernel function with a bandwidth parameter set to 0.3. For each node pair, feature vectors of the source and target nodes, as well as feature vectors of other nodes adjacent to the target node, are extracted. Conditional mutual information is calculated by estimating the log-expected value of the ratio of the joint probability density to the conditional marginal probability density, and the estimation process is numerically integrated over sampling points in the feature space. The spatial distance prior is defined as the reciprocal of the Euclidean distance between node pairs, calculated based on the spatial location index of the nodes; closer node pairs are assigned higher prior weights. Instantaneous causal strength is obtained by multiplying the conditional mutual information by the spatial distance prior. The adaptive threshold is determined based on the distribution of instantaneous causal strength across all node pairs, set to the 75th percentile of the instantaneous causal strength; node pairs below this threshold have their instantaneous causal strength set to zero.

[0088] The causal strength matrix is ​​constructed by nonlinearly combining lagged and instantaneous causal strengths using a weighted geometric mean with weighting coefficients of 0.6 and 0.4, respectively. For edges possessing both lagged and instantaneous causal strengths, their causal strength is the weighted geometric mean of the two; for edges possessing only one type of strength, their causal strength is equal to that strength value multiplied by the corresponding weighting coefficient. The causal strength matrix has a dimension of 15360×15360, and the matrix elements represent the causal strength from row index nodes to column index nodes. The pruning threshold is determined based on the statistics of the causal strength matrix, specifically the mean of the causal strengths plus 0.5 times the standard deviation. All edges with causal strengths below the pruning threshold are removed, and the remaining edges constitute the spatiotemporal causal graph.

[0089] Directed message passing on the spatiotemporal causal graph consists of two phases. First, historical source node messages are aggregated. For each target node at time index t, all edges from source nodes at time index t-1 to that target node are collected. The feature vectors of the source nodes are weighted and summed to obtain the historical message, with the weights being the normalized causal strength values ​​of the corresponding edges. Simultaneous source node messages are also aggregated. All edges from source nodes at time index t to that target node are collected, and the aggregation method is the same as for historical messages. The gating mechanism uses a gated recurrent unit structure. The input is a concatenated vector of historical messages and simultaneous messages, and the output is the fused message. The gated recurrent unit has 960 hidden units. The activation functions for the update and reset gates are sigmoid functions, and the activation function for candidate activations is the tanh function. The fused message is added to the original node feature vector through a residual connection. The resulting vector undergoes layer normalization to obtain enhanced node features, maintaining a dimension of 960.

[0090] This invention accurately identifies temporal dependencies and spatial correlations between nodes through the collaborative calculation of lag causal strength and instantaneous causal strength. The pruning operation effectively removes redundant edges and reduces computational complexity. The causal-guided directed message passing mechanism dynamically aggregates neighbor information based on causal strength, enhancing the ability of node features to represent spatiotemporal evolution patterns.

[0091] In one optional implementation, the steps of mapping the enhanced node features to a continuous spatiotemporal domain, performing feature transformation in the frequency domain using Fourier neural operators, and then performing inverse transformation to obtain a continuous spatiotemporal feature field include:

[0092] Coordinate embedding is performed on the enhanced node features, and the spatial and temporal coordinates of the discrete grid nodes are encoded into positional representations, which are then concatenated with the enhanced node features as input features.

[0093] The input features are transformed in the frequency domain through multiple Fourier layers. Each Fourier layer performs a fast Fourier transform to map the input to the frequency domain space. In the frequency domain space, a learnable frequency domain convolution kernel is used to perform a weighted transformation on different frequency components. Then, an inverse fast Fourier transform is used to map the input back to the spatial domain. The frequency domain convolution kernel captures global spatiotemporal dependencies.

[0094] In the spatial domain, nonlinear activation and residual connection are performed on the frequency domain transformed features to obtain frequency domain enhanced features;

[0095] The dominant causal path is extracted from the spatiotemporal causal graph and represented as a spatiotemporal propagation direction field; a direction-adaptive spatiotemporal kernel function is constructed, the anisotropy of which is determined by the propagation direction field.

[0096] For any query coordinate in the continuous spatiotemporal domain, the contribution weight of each grid node to the query coordinate is calculated based on its spatiotemporal distance from the discrete grid nodes and the spatiotemporal kernel function; the frequency domain enhancement features are weighted and combined according to the contribution weights to obtain the continuous spatiotemporal feature field.

[0097] For example, the dimension of the enhanced node features is 960, and each node corresponds to a grid cell in the spatiotemporal candidate graph. The coordinate embedding is implemented by converting the spatial and temporal coordinates of discrete grid nodes into high-dimensional vector representations through a positional encoding function. The spatial coordinates include the row and column indices of the grid cell in two-dimensional space, ranging from 0 to 31, respectively. The temporal coordinate is the temporal index of the frame in which the node resides, ranging from 0 to 14. Positional encoding uses a sine-cosine encoding method. For the row index i and column index j of the spatial coordinates, a row position vector and a column position vector of dimension 128 are generated. The formula for calculating the k-th element of the row position vector is: sin(i / 10000^(k / 128)) when k is even, and cos(i / 10000^(k / 128)) when k is odd. The column position vector is calculated in the same way as the row position vector, except that i is replaced by j. The time coordinates are encoded to generate a 256-dimensional time position vector. The encoding formula is: sin(t / 5000^(k / 256)) when k is even, and cos(t / 5000^(k / 256)) when k is odd, where t is the time index. The three position vectors are concatenated to obtain a 512-dimensional position representation. This position representation is then concatenated with the augmented node features to form an input feature of 1472 dimensions.

[0098] The input features are fed into multiple Fourier layers for frequency domain transformation, with a maximum of four layers. Each Fourier layer's processing flow includes three stages: frequency domain mapping, frequency domain convolution, and spatial domain mapping. In the frequency domain mapping stage, the input features are first spatially rearranged, recombining the 15360 node features according to their spatial location and time index into a 32×32×15 three-dimensional tensor, with each spatial location corresponding to a feature vector at 15 time steps. The Fast Fourier Transform (FFT) is performed separately in the spatial and temporal dimensions. The spatial FFT uses a two-dimensional FFT, mapping the 32×32 spatial distribution to a 32×32 frequency domain space. The temporal FFT uses a one-dimensional FFT, mapping the 15 time steps to 15 frequency components. The FFT is implemented using the Cooley-Tukey algorithm with single-precision floating-point accuracy. In the frequency domain convolution stage, each frequency component in the frequency domain space corresponds to a complex value. The learnable frequency domain convolution kernel is a complex matrix with a dimension of 1472×1472. Frequency domain convolution is implemented through complex matrix multiplication, multiplying the feature vector at each frequency position with the frequency domain convolution kernel to obtain the transformed frequency domain features. The frequency domain convolution kernel is initialized using the Xavier initialization method, with the real and imaginary parts sampled from a normal distribution with a mean of 0 and a variance of 2 / (input dimension + output dimension). The convolution kernel parameters for different frequency components are learned independently. Low-frequency components in the spatial frequency domain correspond to the global spatial pattern, while high-frequency components correspond to local details. Low-frequency components in the temporal frequency domain correspond to slowly changing trends, while high-frequency components correspond to rapid fluctuations. In the spatial domain mapping stage, the frequency domain features are mapped back to the spatial domain through an inverse fast Fourier transform. The inverse transform is performed sequentially in the temporal and spatial dimensions, outputting a 32×32×15 three-dimensional tensor with the same size as the input, which is then flattened into a feature vector with 15360 nodes.

[0099] The frequency-domain transformed features are nonlinearly activated in the spatial domain using the GELU activation function. The residual connection adds the activated features to the input features of the Fourier layer. Before addition, a linear projection layer adjusts the dimension of the input features to match the activated features. The weight matrix and bias vector of the linear projection layer are 1472×1472. The residual-connected features undergo layer normalization, which calculates the mean and standard deviation along the feature dimension, scaling the features to a distribution with a mean of 0 and a standard deviation of 1. Fourier layers are stacked sequentially, with the output of each layer serving as the input to the next. The output of the final layer is the frequency-domain enhanced feature, maintaining a dimension of 1472.

[0100] The extraction of the dominant causal path is based on the edge weights of the spatiotemporal causal graph, where edge weights represent causal strength. For each node, the edge with the highest causal strength among the incoming edges is selected as the dominant causal edge. The connection direction between the source node and the target node of the dominant causal edge is defined as the propagation direction of the target node. The propagation direction is represented by calculating the spatial position difference and temporal index difference between the source and target nodes. The spatial position difference is the difference between the row and column indices of the source node and the target node, and the temporal index difference is the difference between the time index of the source node and the time index of the target node. The propagation direction vector is a three-dimensional vector containing row, column, and temporal components. The spatiotemporal propagation direction field is a vector field composed of the propagation direction vectors of all nodes, with a dimension of 15360×3.

[0101] The spatiotemporal kernel function employs an anisotropic Gaussian kernel, and its form is determined by the covariance matrix and the distance vector. For node n, its propagation direction vector is denoted as d = (d r ,d c ,d t The components in the row, column, and time directions are represented by ), respectively. The covariance matrix Σ is constructed as follows: Calculate the unit vector u = d / ||d|| of the propagation direction vector, where ||d|| is the Euclidean norm of d. Construct a rotation matrix R such that the first column of R is the unit vector u, and the other two columns are orthogonalized to u using Schmitt orthogonalization. The diagonal elements of the diagonal matrix D are set to the square of the standard deviation of the principal axes σ. parallel 2 =4.0 and the square of the two vertical standard deviations σ perp 2 =0.25. The covariance matrix is ​​calculated using the formula Σ=R×D×R T The calculation yields R, where R is the largest known value. T Let be the transpose of R. During the Schmidt orthogonalization process, the second column is normalized by selecting the coordinate axis vector least related to u as the initial vector and subtracting its projection onto u; the third column is obtained by the cross product of u and the second column and then normalized. The formula for calculating the kernel function value is: K(Δ) = exp(-0.5 × Δ) T ×Σ -1 ×Δ), where Δ is the three-dimensional distance vector between the query coordinates and the node coordinates, Σ -1 Let Δ be the inverse of the covariance matrix. T This is the transpose of the distance vector. Matrix inversion is performed using LU decomposition with a precision of double-precision floating-point numbers.

[0102] In the continuous spatiotemporal domain, the query coordinates are arbitrary real-valued triples, containing spatial row coordinates, spatial column coordinates, and time coordinates. The spatial coordinates of the query coordinates can be any real number between 0 and 31, and the time coordinates can be any real number between 0 and 14. The spatiotemporal distance is calculated using Euclidean distance. For query coordinates q and grid node n, the distance vector Δ = qn is a three-dimensional vector consisting of the differences between the query coordinates and the node coordinates in the row, column, and time dimensions. The contribution weight is calculated using a spatiotemporal kernel function. Substituting the distance vector Δ into the anisotropic Gaussian kernel function of node n, the kernel function value is the unnormalized contribution weight of that node to the query coordinates. After calculating the contribution weights of all grid nodes to the query coordinates, the weights are normalized so that their sum equals 1. The normalization method is: the normalized weight of node n is equal to the kernel function value of that node divided by the sum of the kernel function values ​​of all nodes, where the summation iterates through all nodes. During weighted combination, the frequency domain enhancement feature of each node is multiplied by its normalized contribution weight, and the weighted feature vectors of all nodes are summed to obtain the continuous spatiotemporal feature field at the query coordinates. The feature field has a dimension of 1472.

[0103] For example, the query coordinates are set to spatial row coordinates 15.3, spatial column coordinates 20.7, and time coordinates 9.5. The four grid nodes closest to the query coordinates are located at coordinates (15, 20, 9), (15, 21, 9), (16, 20, 9), and (16, 21, 9), respectively. Additionally, the nodes adjacent in the time dimension are considered: (15, 20, 10), (15, 21, 10), (16, 20, 10), and (16, 21, 10). The propagation direction vector of node (15, 20, 9) is (-0.8, 1.2, -1.0), and the vector norm is 1.75 obtained by calculating the square root of the sum of 0.64, 1.44, and 1.0. The unit vector u = (-0.457, 0.686, -0.571). The principal axis of the constructed covariance matrix is ​​along this direction, and the diagonal elements of the diagonal matrix D are (4.0, 0.25, 0.25). The distance vector between this node and the query coordinates is (0.3, 0.7, 0.5), and the Euclidean norm of the distance vector is 0.91. Substituting into the Gaussian kernel function, the transpose of the quadratic distance vector, multiplied by the inverse of the covariance matrix, and then multiplied by the distance vector, yields 0.58. The kernel function value is then calculated using an exponential function of -0.29, resulting in 0.75. The sum of the kernel function values ​​of the eight neighboring nodes is 4.82, and the contribution weight of node (15, 20, 9) after normalization is 0.156. The mean value of each channel of the frequency domain enhancement feature of this node is 0.42, and the weighted contribution value is 0.066. After summing the weighted features of the eight nodes, the mean value of each channel of the continuous spatiotemporal feature field at the query coordinates is 0.39.

[0104] This invention captures global spatiotemporal dependencies in the frequency domain using Fourier neural operators, avoiding the limitations of local receptive fields in convolutional neural networks and improving the ability to model long-range correlations. The orientation-adaptive spatiotemporal kernel function dynamically adjusts the interpolation weights according to the causal propagation direction, enabling the continuous spatiotemporal feature field to accurately reflect causally driven spatiotemporal evolution patterns, providing high-quality feature representations for queries at arbitrary resolutions.

[0105] In one alternative implementation, the step of calculating the frequency domain enhancement features includes:

[0106] The input features are subjected to Fast Fourier Transform along the spatial and temporal dimensions to obtain a spectral representation. Based on the energy distribution of the frequency components in the spectral representation, the spectrum is divided into different frequency regions, and feature transformation is performed on different frequency regions using frequency domain convolution kernels with different channel capacities.

[0107] The spectrum representations after transformation of different frequency regions are spliced ​​in the frequency domain space, and different frequency components are adaptively weighted through a frequency domain attention mechanism, which calculates the attention weights based on the spectral energy distribution.

[0108] Perform an inverse fast Fourier transform on the weighted spectral representation to map it back to the spatial domain, and obtain the frequency domain transformed features;

[0109] The frequency-domain transformed features are nonlinearly activated; the activated features are then residually connected with the input features to obtain frequency-domain enhanced features.

[0110] For example, the input features are transformed into a spectral representation through Fast Fourier Transform (FFT) in both spatial and temporal dimensions. The spatial dimension of the input features corresponds to a 32×32 grid layout, the temporal dimension corresponds to 15 time steps, and the feature channel dimension is 1472. The spatial FFT uses a two-dimensional Fourier transform, mapping the 32×32 spatial distribution to a 32×32 frequency domain space. Each position in the frequency domain space corresponds to a spatial frequency component, identified by row and column frequency indices, both ranging from 0 to 31. The temporal FFT uses a one-dimensional Fourier transform, mapping the 15 time steps to 15 frequency components, with frequency indices ranging from 0 to 14. The FFT is implemented using the Cooley-Tukey algorithm, with a computational precision of single-precision floating-point numbers. The transformed spectrum is represented as a complex tensor with dimensions of 32×32×15×1472, where the feature vector at each frequency position contains both real and imaginary components.

[0111] The energy distribution of frequency components in the spectral representation is obtained by calculating the square of the complex modulus at each frequency position. For spatial frequency components, the Euclidean distance between the row and column frequency indices is defined as the amplitude of the spatial frequency. Frequency components with amplitudes less than 5 are classified as low-frequency regions, those with amplitudes between 5 and 15 as mid-frequency regions, and those with amplitudes greater than 15 as high-frequency regions. For time frequency components, frequency indices less than 3 are classified as low-frequency regions, those with frequency indices between 3 and 8 as mid-frequency regions, and those with frequency indices greater than 8 as high-frequency regions. After combining the spatial and time frequency regional divisions, the spectrum is divided into 9 frequency regions, corresponding to spatial low-frequency-time low-frequency, spatial low-frequency-time mid-frequency, spatial low-frequency-time high-frequency, spatial mid-frequency-time low-frequency, spatial mid-frequency-time mid-frequency, spatial mid-frequency-time high-frequency, spatial high-frequency-time low-frequency, spatial high-frequency-time mid-frequency, and spatial high-frequency-time high-frequency.

[0112] Feature transformation is performed using frequency-domain convolution kernels with different channel capacities in different frequency regions. The low-frequency region, corresponding to global spatial patterns and slow temporal changes, uses a frequency-domain convolution kernel with a channel capacity of 2048; the kernel is a complex matrix with a matrix dimension of 1472×2048. The mid-frequency region, corresponding to medium-scale spatial structures and medium-speed temporal evolution, uses a frequency-domain convolution kernel with a channel capacity of 1472; the matrix dimension is 1472×1472. The high-frequency region, corresponding to local details and rapid fluctuations, uses a frequency-domain convolution kernel with a channel capacity of 1024; the matrix dimension is 1472×1024. Frequency-domain convolution is achieved through complex matrix multiplication; the feature vector at each frequency position is multiplied by the corresponding frequency-domain convolution kernel to obtain the transformed frequency-domain features. The implementation of complex matrix multiplication is as follows: for a complex matrix A and a complex vector x, the real part of the result y is equal to the product of the real part of A and the real part of x, minus the product of the imaginary part of A and the imaginary part of x. The imaginary part is equal to the product of the real part of A and the imaginary part of x, plus the product of the imaginary part of A and the real part of x. The frequency domain convolution kernel is initialized using the Xavier initialization method, with the real and imaginary parts sampled from a normal distribution with a mean of 0 and a variance of 2 divided by the sum of the input and output dimensions.

[0113] The frequency domain features transformed from different frequency regions are concatenated along the channel dimension. The number of feature channels after transformation is 2048 in the low-frequency region, 1472 in the mid-frequency region, and 1024 in the high-frequency region. The feature vector dimension of each frequency position after concatenation is 4544. The frequency domain attention mechanism calculates attention weights based on the spectral energy distribution and adaptively weights different frequency components. The spectral energy is calculated by taking the square of the complex modulus of the concatenated feature vector for each frequency position and summing it over all feature channels to obtain the total energy at that frequency position. The total energy at all frequency positions constitutes an energy distribution map with a dimension of 32×32×15. The attention weights are calculated by applying the Softmax function to the energy distribution map. The Softmax function normalizes the data at all frequency positions so that the sum of the attention weights at all frequency positions equals 1. The weighting operation multiplies the feature vector at each frequency position by its corresponding attention weight to obtain the weighted spectral representation, maintaining the dimension of 32×32×15×4544.

[0114] The weighted spectrum is then mapped back to the spatial domain using an inverse fast Fourier transform (IFT). The inverse transform is performed sequentially along the time and spatial dimensions. The time-dimension IFT maps the 15 frequency components back to 15 time steps, while the spatial-dimension IFT maps the 32×32 frequency domain space back to a 32×32 spatial distribution. The output of the inverse transform is a real tensor with dimensions 32×32×15×4544. Flattening the tensor by spatial location and time index yields a feature vector with 15360 nodes, each with a dimension of 4544. This feature vector represents the features after the frequency domain transformation.

[0115] The features after frequency domain transformation are nonlinearly activated using the GELU activation function. The GELU activation function is calculated as follows: for an input value x, the output value equals x multiplied by the value of the cumulative distribution function of the standard normal distribution at x. The cumulative distribution function is approximated using an error function, specifically calculated using the hyperbolic tangent function, which is 0.5 multiplied by 1 plus the hyperbolic tangent function. The input to the hyperbolic tangent function is the sum of 0.7978845608 multiplied by x plus 0.044715 multiplied by the cube of x. The activated feature dimension remains 15360×4544.

[0116] The activated features are residually concatenated with the input features. Before the concatenation, a linear projection layer is used to adjust the dimension of the input features from 1472 to 4544. The weight matrix of the linear projection layer has a dimension of 1472×4544, and the bias vector has a dimension of 4544. The weight matrix is ​​initialized using the Kaiming initialization method, sampling from a normal distribution with a mean of 0 and a variance of 2 divided by the input dimension. The bias vector is initialized to 0. The linear projection is calculated by multiplying the input feature matrix by the weight matrix and adding the bias vector to obtain the projected features, which have a dimension of 15360×4544. The projected features are then added element-wise to the activated features to obtain the residually concatenated features, maintaining the dimension of 15360×4544. The residually concatenated features undergo layer normalization, which calculates the mean and standard deviation along the feature dimension, scaling the features to a distribution with a mean of 0 and a standard deviation of 1. Layer normalization is calculated as follows: for each node's feature vector, subtract the mean of all elements in the vector, then divide by the standard deviation of all elements plus a small constant 1e-5 to prevent division by zero. The normalized feature is then multiplied by a learnable scaling parameter gamma and a learnable offset parameter beta. Both gamma and beta have a dimension of 4544, initialized to 1 for gamma and 0 for beta. The normalized feature is the frequency domain enhanced feature, with a dimension of 15360×4544.

[0117] This invention uses a frequency domain attention mechanism to adaptively weight different frequency components according to energy distribution. This frequency-adaptive feature transformation strategy improves the modeling accuracy of multi-scale spatiotemporal patterns. At the same time, it preserves the original feature information through residual connections, thereby enhancing stability and generalization ability.

[0118] In one optional implementation, the step of iteratively denoising the initial noise field using the continuous spatiotemporal feature field as a condition, predicting and removing noise components in each iteration until convergence is achieved to obtain the depth field distribution map and the concentration field distribution map includes:

[0119] An initial noise field is randomly generated, following a zero-mean, unit-variance distribution. The spatiotemporal dimension of the initial noise field is correlated with the target depth field distribution and concentration field distribution. Figure 1 To;

[0120] The initial noise field is iteratively denoised. In each iteration, the current noise field and the continuous spatiotemporal feature field are interacted through a cross-attention mechanism to obtain a conditionally enhanced noise field representation.

[0121] The noise component is predicted by a denoising prediction network based on the conditionally enhanced noise field representation and the current iteration step index; temporal causal constraints are extracted from the spatiotemporal causal graph and encoded as a causal mask matrix to constrain the noise component; the noise component constrained by the causal mask is subtracted from the current noise field to obtain the denoised noise field; the noise removal step size is dynamically adjusted according to the current iteration step index, and the denoised noise field is used as the current noise field for the next iteration;

[0122] The iteration is terminated when the rate of change of the noise field between adjacent iterations is lower than the preset rate of change threshold or the maximum number of iterations is reached; the final denoised noise field is decoupled into a depth field distribution map and a concentration field distribution map.

[0123] For example, the initial noise field is generated using a random number generator, employing a standard normal distribution with zero mean and unit variance. The random number generator uses the Mason twitch algorithm, with a fixed seed value of 42 to ensure reproducibility. When randomness is required, the seed value can be dynamically generated based on timestamps. The spatiotemporal dimensions of the initial noise field are consistent with the target depth and concentration field distribution maps, with a spatial dimension of 256×256 grid resolution, a temporal dimension of 30 time steps, and a physical field dimension of 2, corresponding to the depth and concentration fields respectively. The data structure of the initial noise field is a four-dimensional tensor with dimensions of 256×256×30×2. Each grid cell contains two independently sampled noise values ​​at each time step. The noise values ​​are sampled with single-precision floating-point numbers, theoretically ranging from negative infinity to positive infinity; in actual sampling, 99.7% of the values ​​fall within the range of -3 to +3. The mean of each channel of the generated initial noise field is controlled between -0.01 and 0.01, and the standard deviation is controlled between 0.98 and 1.02. This is achieved by standardizing the generated noise field by subtracting the actual mean and then dividing by the actual standard deviation.

[0124] The iterative denoising process is implemented using a loop structure, processing the current noise field in each iteration. The iteration step index decreases from the maximum iteration step number to 0, with the maximum iteration step number set to 50 and a decrease step size of 1. The current noise field and the continuous spatiotemporal feature field interact through a cross-attention mechanism. The continuous spatiotemporal feature field has a dimension of 256×256×30×1472, while the current noise field has a dimension of 256×256×30×2. The cross-attention mechanism is implemented by using the current noise field as the query and the continuous spatiotemporal feature field as the key and value. The query is mapped from 2 to 512 through a linear projection layer, with the weight matrix of the linear projection layer having a dimension of 2×512 and the bias vector having a dimension of 512. The key and value are mapped from 1472 to 512 through their respective linear projection layers, with the weight matrix of both projection layers having a dimension of 1472×512 and the bias vector having a dimension of 512. The dimensions of the projected query, key, and value are all 256×256×30×512. Attention weights are calculated by performing matrix multiplication of the query and the transpose of the key, dividing the result by the square root of 512 (22.627), and then normalizing using the Softmax function. Normalization is performed along the spatial-temporal dimensions of the key, ensuring that the sum of the attention weights for each query position to all key positions equals 1. Weighted summation of attention is achieved by multiplying the attention weight matrix with the value matrix, yielding the interacting features with dimensions 256×256×30×512. These interacting features are then fused with the current noise field via residual connections. Before fusion, the interacting features are mapped back to dimension 2 using a linear projection layer. The weight matrix of the linear projection layer has dimensions 512×2, and the bias vector has dimensions 2. The mapped interacting features are then added element-wise to the current noise field to obtain a conditionally enhanced noise field representation with dimensions 256×256×30×2.

[0125] The denoising prediction network predicts noise components based on a conditionally enhanced noise field representation and the current iteration step index. The iteration step index is converted into a temporal embedding vector using sine and cosine positional encoding, with an encoding dimension of 128. The encoding method is as follows: for iteration step index t, the k-th element of the temporal embedding vector is a sine function when k is even, with the input to the sine function being t divided by 10000 and k divided by 128; when k is odd, a cosine function is used, with the input calculated in the same way. The temporal embedding vector is mapped to a dimension of 512 through two fully connected layers, with the SiLU activation function used between the two layers. The SiLU activation function is calculated by multiplying the input value by the Sigmoid function at that input value. The conditionally enhanced noise field representation and the temporal embedding vector are added together via a broadcast mechanism, with the temporal embedding vector copied in both spatial and temporal dimensions to match the dimension of the noise field. The added features are then fed into the U-Net architecture denoising prediction network, which contains four downsampling layers and four upsampling layers. Downsampling layers are implemented using convolutions with a stride of 2, halving the spatial resolution at each layer, with channel numbers of 64, 128, 256, and 512 respectively. Upsampling layers are implemented using transposed convolutions, doubling the spatial resolution at each layer, with channel numbers of 256, 128, 64, and 2 respectively. Each convolutional layer is followed by batch normalization and a ReLU activation function. Batch normalization calculates the mean and standard deviation in both the batch and spatial dimensions. The normalized features are multiplied by a learnable scaling parameter gamma and a learnable offset parameter beta. Features are passed between downsampling and upsampling layers via skip connections, which concatenate the downsampling and upsampling features at the corresponding scale along the channel dimension. The noise component of the network output has dimensions of 256×256×30×2, representing the predicted depth field noise and concentration field noise.

[0126] Temporal causality constraints are extracted from a spatiotemporal causal graph, which records the causal relationships between nodes and their temporal order. The causal mask matrix is ​​constructed such that, for a grid cell at the target time step t, only causal influences from grid cells at time steps less than or equal to t are allowed. The causal mask matrix is ​​30×30 in dimension. Matrix elements are 0 when the row index is greater than the column index, indicating that future time steps cannot influence past time steps; and 1 when the row index is less than or equal to the column index, indicating allowed causal directions. Noise components are constrained by element-wise multiplication with the causal mask matrix. Before multiplication, the causal mask matrix is ​​broadcast in both the spatial and physical dimensions to match the dimensions of the noise components. After constraint, noise components that violate causality are set to zero at the locations where causality is violated, while noise predictions that conform to causality are preserved.

[0127] The denoised noise field is obtained by subtracting the noise component constrained by the causality mask from the current noise field. The step size for noise removal is dynamically adjusted based on the current iteration step index. The adjustment method is that the step size equals 1 divided by the current iteration step index plus 1. As the current iteration step index decreases from 50 to 1, the step size gradually increases from 0.0196 to 0.5. The actual amount of noise removed is obtained by multiplying the step size by the constrained noise component. The denoised noise field is then subtracted from the current noise field. The denoised noise field serves as the current noise field for the next iteration, and the next iteration begins after the iteration step index is decremented by 1.

[0128] The rate of change of the noise field between adjacent iterations is obtained by dividing the L2 norm of the difference between the denoised noise field and the current noise field by the L2 norm of the current noise field. The L2 norm is calculated by summing the squares of the values ​​of all grid cells, all time steps, and all physical fields, and then taking the square root. Iteration terminates when the rate of change falls below a preset threshold of 0.001. Iteration terminates again when the maximum number of iterations (50) is reached, regardless of whether the rate of change condition is met. After iteration termination, the final denoised noise field has dimensions of 256×256×30×2.

[0129] The final denoised noise field is decoupled into a depth field distribution map and a concentration field distribution map. The decoupling method involves extracting the first channel from the physical field dimension of the four-dimensional tensor as the depth field distribution map and the second channel as the concentration field distribution map. The depth field distribution map has dimensions of 256×256×30. Its numerical range is linearly mapped from the numerical range of the denoised noise field to a physical depth range of 0 to 10 meters. The mapping formula is: depth value equals noise field value minus the minimum noise field value, multiplied by 10 and divided by the result of the maximum and minimum noise field values. The concentration field distribution map also has dimensions of 256×256×30. Its numerical range is linearly mapped to a physical concentration range of 0 to 100 milligrams per liter. The mapping formula is the same as the depth field, except that the coefficient 10 is replaced with 100.

[0130] This invention utilizes a cross-attention mechanism to deeply fuse continuous spatiotemporal feature fields with noise fields, enabling the denoising process to fully leverage causal-driven spatiotemporal evolution information. The constraint of the noise components by the causal mask matrix ensures that the generated results conform to temporal causal logic, avoiding the non-physical influence of future information on past states. The dynamically adjusted noise removal step size rapidly approaches the target in the early stages of iteration and is finely adjusted in later stages, improving convergence speed and generation quality, resulting in more accurate and stable predictions of the depth and concentration fields.

[0131] In one optional implementation, the step of generating a three-dimensional spatial distribution model of the gas leak based on the depth field distribution map and the concentration field distribution map includes:

[0132] Spatial geometry reconstruction is performed on the depth field distribution map to convert depth values ​​into three-dimensional spatial coordinates, establishing a three-dimensional geometric structure of the leakage scene. Concentration values ​​in the concentration field distribution map are mapped to corresponding spatial positions in the three-dimensional geometric structure, establishing a correspondence between concentration values ​​and three-dimensional spatial coordinates through voxelization. Spatial interpolation is performed on the voxelization representation, calculating interpolated concentration values ​​at continuous spatial positions between voxel grids to obtain a continuous three-dimensional concentration distribution. Concentration isosurfaces are set based on the three-dimensional concentration distribution, and three-dimensional surfaces corresponding to different concentration thresholds are extracted. These three-dimensional surfaces represent the spatial boundaries of gas diffusion.

[0133] The three-dimensional geometric structure, the three-dimensional concentration distribution, and the three-dimensional surface are fused to generate a three-dimensional spatial distribution model of gas leakage.

[0134] For example, the data structure of the depth field distribution map is a three-dimensional array with dimensions of 256×256×30, corresponding to spatial row index, spatial column index, and time index t, respectively. Each grid cell stores a depth value, with a physical range of 0 to 10 meters. The origin of the three-dimensional spatial coordinate system is set at the lower left corner of the ground in the monitoring area, with the X-axis along the spatial column direction, the Y-axis along the spatial row direction, and the Z-axis vertically upward. The spatial row index i and spatial column index j are converted to physical coordinates by multiplying by the spatial resolution, which is set to 0.1 meters per grid. The conversion formula is: X coordinate = j × 0.1, Y coordinate = i × 0.1. The depth value d is directly used as the Z coordinate. For the grid cells with indices i, j, and t in the depth field distribution map, their corresponding three-dimensional spatial coordinates are X = j × 0.1, Y = i × 0.1, and Z equal to the depth value d. The three-dimensional geometric structure is represented by a point cloud data structure. The point cloud contains 256×256×30 points, and each point stores three coordinate components: X, Y, and Z. The data type of the coordinate components is single-precision floating-point numbers. The spatial range of the point cloud is 0 to 25.5 meters in the X direction, 0 to 25.5 meters in the Y direction, and 0 to 10 meters in the Z direction.

[0135] Dimensional and depth field distribution of concentration field distribution map Figure 1The voxel grid is 256×256×30, with each cell storing a concentration value ranging from 0 to 100 mg / L. For cells with indices i, j, and t in the concentration field distribution map, their concentration value c is mapped to a point in the 3D geometry with coordinates X = j × 0.1, Y = i × 0.1, and Z equal to the depth value d. Voxelization is achieved by dividing the 3D space into regular cubic grids with a resolution of 0.1 meters, consistent with the spatial resolution of the depth field. The voxel grid has dimensions of 256×256×100, corresponding to 256 voxels in the X direction, 256 voxels in the Y direction, and 100 voxels in the Z direction. The number of voxels in the Z direction is obtained by dividing the maximum depth of 10 meters by the voxel resolution of 0.1 meters. Each voxel stores the concentration value at its center, and the coordinates of the voxel center are the voxel index multiplied by the voxel resolution plus half the voxel resolution. The concentration value is determined by consulting the depth field distribution map and the concentration field distribution map. For voxels with voxel indices ix, iy, and iz, their center coordinates are X = ix × 0.1 + 0.05, Y = iy × 0.1 + 0.05, and Z = iz × 0.1 + 0.05. The corresponding grid cell indices i and j are determined using the X and Y coordinates. i is equal to the Y coordinate divided by 0.1 and rounded down, and j is equal to the X coordinate divided by 0.1 and rounded down. The depth value d is obtained by consulting the depth field distribution map using the grid cell indices i, j, and time index t. The Z coordinate is compared with the depth value d. If the Z coordinate is less than or equal to the depth value d, the voxel is located inside the leakage scene, and the concentration value is set to the concentration value at indices i, j, and t in the concentration field distribution map. If the Z coordinate is greater than the depth value d, the voxel is located outside the leakage scene, and the concentration value is set to 0.

[0136] Voxelization is used for spatial interpolation, calculating interpolated concentration values ​​at continuous spatial locations between voxel grids. The spatial interpolation employs a trilinear interpolation method. For the queried 3D spatial coordinates X, Y, and Z, the corresponding voxel index floating-point number is calculated by dividing by the voxel resolution of 0.1: X / 0.1 for the X direction, Y / 0.1 for the Y direction, and Z / 0.1 for the Z direction. The floating-point index is decomposed into an integer part and a fractional part. The integer part is obtained by rounding down, and the fractional part is obtained by subtracting the integer part. The integer part determines the eight vertex voxels surrounding the queried location. The indices of these eight vertices are combinations of the integer part and the integer part plus 1 in the X, Y, and Z directions. The concentration values ​​of the eight vertices are looked up in the voxelization representation. If a vertex index exceeds the range of the voxel grid, its concentration value is set to 0. The trilinear interpolation is calculated by performing linear interpolation along the X direction between four pairs of vertices, with the interpolation weight being the fractional part of the X direction, resulting in four interpolation results. Linear interpolation is performed along the Y-axis between two pairs of interpolation results, with the interpolation weight being the decimal part of the Y-axis, resulting in two interpolation results. Linear interpolation is then performed along the Z-axis between the last two interpolation results, with the interpolation weight being the decimal part of the Z-axis, resulting in the final interpolated concentration value. The linear interpolation is calculated as follows: the interpolation result equals the first value × (1 - weight) + the second value × weight. The interpolated concentration value is a single-precision floating-point number, ranging from 0 to 100 mg / L. Continuous three-dimensional concentration distribution is achieved by calculating the interpolated concentration value for any queried three-dimensional spatial coordinate. The precision of the query coordinate is 0.01 meters, and the calculation error of the interpolated concentration value is less than 0.5 mg / L.

[0137] A 3D concentration distribution isosurface was established to extract 3D surfaces corresponding to different concentration thresholds. The concentration thresholds were set to four levels: 10 mg / L, 20 mg / L, 50 mg / L, and 80 mg / L, representing low-concentration, low-to-medium-concentration, medium-to-high-concentration, and high-concentration boundaries, respectively. The 3D surface extraction employed a moving cube algorithm, which traverses each cube cell in a voxel mesh. Each cube cell consists of eight adjacent voxel vertices. For each cube cell, the concentration values ​​of the eight vertices were compared with the concentration threshold. Vertices with concentrations greater than the threshold were marked as 1, and those less than or equal to the threshold were marked as 0. The combinations of these eight vertex markings formed an 8-bit binary number, ranging from 0 to 255. A predefined triangular facet configuration table was then used to look up the binary number. This table stored the vertex positions and connections of triangular faces corresponding to 256 vertex marking combinations. The vertices of the triangular faces were located on the edges of the cube cell, and linear interpolation was used to calculate the precise positions on the edges where the concentration value equaled the threshold. For an edge of a cube element, the concentration values ​​at its two endpoints are c1 and c2, respectively, and the concentration threshold is th. The proportion of the interpolation position on the edge is (th-c1) / (c2-c1). The 3D coordinates of the interpolation position are the coordinates of the first endpoint + the proportion × (the coordinates of the second endpoint - the coordinates of the first endpoint). After the moving cube algorithm traverses all cube elements, the resulting set of triangular patches constitutes a 3D surface of the concentration isosurface. The data structure of the 3D surface is a triangular mesh, storing a list of vertex coordinates and a list of triangular patch indices. The vertex coordinate list contains the X, Y, and Z coordinates of all triangular patch vertices, and the triangular patch index list contains the indices of the three vertices of each triangular patch in the vertex coordinate list.

[0138] The fusion of 3D geometry, 3D concentration distribution, and 3D surfaces is achieved by organizing these three components into a unified data structure. This data structure contains three fields: a geometry field, a concentration distribution field, and an isosurface field. The geometry field stores point cloud data, which is a 256×256×30 point coordinate array. Each point contains X, Y, and Z coordinate components and a time index t. The concentration distribution field stores a voxel representation, which is a 256×256×100×30 four-dimensional array. Each voxel contains a concentration value and a time index t. The isosurface field stores triangular meshes corresponding to four concentration thresholds. Each triangular mesh contains a list of vertex coordinates, a list of triangular facet indices, and a time index t. The 3D spatial distribution model associates the geometry, concentration distribution, and isosurface at different time steps through the time index t. The time index ranges from 0 to 29, corresponding to 30 time steps. The model is stored as a structured binary file. The file header contains the model's metadata, including spatial resolution of 0.1 meters, voxel resolution of 0.1 meters, time steps of 30, and a list of concentration thresholds of 10, 20, 50, and 80 mg / L. The file body contains geometric structure data, concentration distribution data, and isosurface data sorted by time index. Each time step's data block contains the point cloud coordinates, voxel concentration values, and triangular mesh vertices and indices for that time step. Model loading is achieved by reading the file header, parsing the metadata, locating the data blocks according to the time index, and loading the data for the corresponding time step.

[0139] This implementation establishes a precise mapping between concentration values ​​and three-dimensional space through voxelization representation, achieves a smooth concentration distribution in continuous space through trilinear interpolation, and clearly displays the diffusion boundaries of different concentration thresholds through isosurfaces extracted by the moving cube algorithm. The three-dimensional spatial distribution model integrates geometric structure, concentration distribution, and isosurfaces, providing complete spatial information support for the visualization, risk assessment, and emergency decision-making of leakage scenarios.

[0140] A second aspect of the present invention provides a three-dimensional reconstruction system for gas leaks based on thermo-optical fusion, comprising:

[0141] The first unit is used to extract multi-scale feature maps from the thermal radiation image sequence and visible light image sequence of the acquired target area, respectively. The local mutual information distribution of thermal radiation features and visible light features is calculated at each scale, weighted and fused, and then spliced ​​to generate a multi-scale feature representation.

[0142] The second unit is used to construct a spatiotemporal candidate graph by spatial grid partitioning of the multi-scale feature representations of multiple consecutive frames, calculate the causal strength between nodes, prune according to the causal strength to obtain a spatiotemporal causal graph, and perform causal-guided directed message passing on the spatiotemporal causal graph to obtain enhanced node features.

[0143] The third unit is used to map the enhanced node features to the continuous spatiotemporal domain, and then perform feature transformation and inverse transformation in the frequency domain using Fourier neural operators to obtain the continuous spatiotemporal feature field.

[0144] The fourth unit is used to iteratively denoise the initial noise field based on the continuous spatiotemporal feature field. In each iteration, noise components are predicted and removed until convergence is obtained to obtain the depth field distribution map and the concentration field distribution map.

[0145] The fifth unit is used to generate a three-dimensional spatial distribution model of gas leakage based on the depth field distribution map and the concentration field distribution map.

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

[0147] processor;

[0148] Memory used to store processor-executable instructions;

[0149] The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

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

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

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

Claims

1. A method for three-dimensional reconstruction of gas leakage based on thermo-optical fusion, characterized in that, include: Multi-scale feature maps are extracted from the thermal radiation image sequence and visible light image sequence of the acquired target area, respectively. The local mutual information distribution of thermal radiation features and visible light features is calculated at each scale, weighted and fused, and then stitched together to generate a multi-scale feature representation. The multi-scale feature representations of consecutive frames are spatially gridded to construct a spatiotemporal candidate graph. The causal strength between nodes is calculated. The spatiotemporal causal graph is obtained by pruning according to the causal strength. Causal-guided directed message passing is performed on the spatiotemporal causal graph to obtain enhanced node features. The enhanced node features are mapped to the continuous spatiotemporal domain, and after feature transformation in the frequency domain by Fourier neural operators, an inverse transformation is performed to obtain the continuous spatiotemporal feature field. Using the continuous spatiotemporal feature field as a condition, the initial noise field is iteratively denoised. Each iteration predicts and removes noise components until convergence is achieved to obtain the depth field distribution map and the concentration field distribution map. Based on the depth field distribution map and the concentration field distribution map, a three-dimensional spatial distribution model of the gas leak is generated; The steps of constructing a spatiotemporal candidate graph by spatially dividing the multi-scale feature representations of multiple consecutive frames into a grid, calculating the causal strength between nodes, pruning according to the causal strength to obtain a spatiotemporal causal graph, and performing causally guided directed message passing on the spatiotemporal causal graph to obtain enhanced node features include: The multi-scale feature representations of multiple consecutive frames are spatially gridded, each grid cell is mapped to a node on the time axis and candidate edges are established to construct a spatiotemporal candidate graph; Extract historical state sequences from the nodes in the spatiotemporal candidate graph; The baseline error predicted using only the historical state of the target node, the source error predicted by adding the historical state of the source node, and the conditional error predicted by adding the historical state of other nodes are calculated separately. When the difference between the baseline error and the source error exceeds the error improvement threshold and the difference between the source error and the conditional error is lower than the conditional influence threshold, the difference between the baseline error and the source error is normalized and used as the hysteresis causality strength. The conditional mutual information between node pairs is calculated using the kernel density estimation method, and the instantaneous causal strength is obtained by combining spatial distance priors and adaptive threshold filtering. The causal intensity and the instantaneous causal intensity are nonlinearly combined to obtain a causal intensity matrix. The pruning threshold is determined based on the statistics of the causal intensity matrix, and low-intensity edges are removed to obtain a spatiotemporal causal graph. On the spatiotemporal causal graph, historical source node messages and simultaneous source node messages are aggregated respectively, and after being fused through a gating mechanism, they are combined with the original node features to obtain enhanced node features.

2. The method according to claim 1, characterized in that, The steps for generating a multi-scale feature representation by extracting multi-scale feature maps from the acquired thermal radiation image sequence and visible light image sequence of the target area, calculating the local mutual information distribution of thermal radiation features and visible light features at each scale, weighting and fusing them, and then stitching them together include: The thermal radiation image sequence and the visible light image sequence are downsampled by a multi-layer convolutional network. After each downsampling, the thermal radiation feature map and the visible light feature map at the current scale are extracted to obtain multiple scale levels from high resolution to low resolution. At each scale level, for the corresponding spatial locations of the thermal radiation feature map and the visible light feature map, the KL divergence between the joint probability distribution and the marginal probability distribution of the feature vectors is calculated as the local mutual information value, and the local mutual information value is normalized and used as the fusion weight. The thermal radiation feature vector and visible light feature vector at the same spatial location are weighted and summed using the fusion weights to obtain weighted fusion feature maps at each scale. The weighted fusion feature maps of different scales are upsampled to the same spatial size and then stitched together in the channel dimension to generate a multi-scale feature representation.

3. The method according to claim 1, characterized in that, The steps of mapping the enhanced node features to a continuous spatiotemporal domain, performing feature transformation in the frequency domain using Fourier neural operators, and then performing inverse transformation to obtain a continuous spatiotemporal feature field include: Coordinate embedding is performed on the enhanced node features, and the spatial and temporal coordinates of the discrete grid nodes are encoded into positional representations, which are then concatenated with the enhanced node features as input features. The input features are transformed in the frequency domain through multiple Fourier layers. Each Fourier layer performs a fast Fourier transform to map the input to the frequency domain space. In the frequency domain space, a learnable frequency domain convolution kernel is used to perform a weighted transformation on different frequency components. Then, an inverse fast Fourier transform is used to map the input back to the spatial domain. The frequency domain convolution kernel captures global spatiotemporal dependencies. In the spatial domain, nonlinear activation and residual connection are performed on the frequency domain transformed features to obtain frequency domain enhanced features; Extract the dominant causal path from the spatiotemporal causal graph and represent the dominant causal path as a spatiotemporal propagation direction field; Construct a direction-adaptive spatiotemporal kernel function, wherein the anisotropy of the spatiotemporal kernel function is determined by the propagation direction field; For any query coordinate in a continuous spatiotemporal domain, the contribution weight of each grid node to the query coordinate is calculated based on its spatiotemporal distance from discrete grid nodes, combined with the spatiotemporal kernel function. The frequency domain enhancement features are weighted and combined according to the contribution weights to obtain a continuous spatiotemporal feature field.

4. The method according to claim 3, characterized in that, The steps for calculating frequency domain enhancement features include: The input features are subjected to Fast Fourier Transform along both the spatial and temporal dimensions to obtain a spectral representation; Based on the energy distribution of frequency components in the spectral representation, the spectrum is divided into different frequency regions, and feature transformation is performed on different frequency regions using frequency domain convolution kernels with different channel capacities. The spectrum representations after transformation of different frequency regions are spliced ​​in the frequency domain space, and different frequency components are adaptively weighted through a frequency domain attention mechanism, which calculates the attention weights based on the spectral energy distribution. Perform an inverse fast Fourier transform on the weighted spectral representation to map it back to the spatial domain, and obtain the frequency domain transformed features; The frequency-domain transformed features are then subjected to nonlinear activation. The activated features are residually concatenated with the input features to obtain frequency domain enhanced features.

5. The method according to claim 1, characterized in that, The steps of iteratively denoising the initial noise field using the continuous spatiotemporal feature field as a condition, predicting and removing noise components in each iteration until convergence is achieved to obtain the depth field distribution map and the concentration field distribution map, include: An initial noise field is randomly generated, which follows a zero-mean, unit-variance distribution. The spatiotemporal dimensions of the initial noise field are consistent with the target depth field distribution map and concentration field distribution map. The initial noise field is iteratively denoised. In each iteration, the current noise field and the continuous spatiotemporal feature field are interacted through a cross-attention mechanism to obtain a conditionally enhanced noise field representation. The noise components are predicted by a denoising prediction network based on the conditionally enhanced noise field representation and the current iteration step index; Temporal causal constraints are extracted from the spatiotemporal causal graph and encoded into a causal mask matrix to constrain the noise components; Subtract the noise component constrained by the causal mask from the current noise field to obtain the denoised noise field. Dynamically adjust the noise removal step size according to the current iteration step index, and use the denoised noise field as the current noise field for the next iteration. The iteration is terminated when the rate of change of the noise field between adjacent iterations is lower than the preset rate of change threshold or the maximum number of iterations is reached. The final denoised noise field is decoupled into a depth field distribution map and a concentration field distribution map.

6. The method according to claim 1, characterized in that, The steps for generating a three-dimensional spatial distribution model of gas leakage based on the depth field distribution map and the concentration field distribution map include: Spatial geometry reconstruction is performed on the depth field distribution map to convert the depth values ​​into three-dimensional spatial coordinates and establish the three-dimensional geometric structure of the leakage scene; The concentration values ​​in the concentration field distribution map are mapped to the corresponding spatial positions of the three-dimensional geometric structure, and the correspondence between the concentration values ​​and the three-dimensional spatial coordinates is established through voxelization. Spatial interpolation is performed on the voxelized representation, and interpolated concentration values ​​are calculated at continuous spatial locations between voxel grids to obtain a continuous three-dimensional concentration distribution; Based on the three-dimensional concentration distribution, a concentration isosurface is set, and three-dimensional surfaces corresponding to different concentration thresholds are extracted. The three-dimensional surfaces represent the spatial boundaries of gas diffusion. The three-dimensional geometric structure, the three-dimensional concentration distribution, and the three-dimensional surface are fused to generate a three-dimensional spatial distribution model of gas leakage.

7. A three-dimensional reconstruction system for gas leaks based on thermo-optical fusion, used to implement the method of any one of claims 1-6, characterized in that, include: The first unit is used to extract multi-scale feature maps from the thermal radiation image sequence and visible light image sequence of the acquired target area, respectively. The local mutual information distribution of thermal radiation features and visible light features is calculated at each scale, weighted and fused, and then spliced ​​to generate a multi-scale feature representation. The second unit is used to construct a spatiotemporal candidate graph by spatial grid partitioning of the multi-scale feature representations of multiple consecutive frames, calculate the causal strength between nodes, prune according to the causal strength to obtain a spatiotemporal causal graph, and perform causal-guided directed message passing on the spatiotemporal causal graph to obtain enhanced node features. The third unit is used to map the enhanced node features to the continuous spatiotemporal domain, and then perform feature transformation in the frequency domain by Fourier neural operators and inverse transformation to obtain the continuous spatiotemporal feature field. The fourth unit is used to iteratively denoise the initial noise field based on the continuous spatiotemporal feature field. In each iteration, noise components are predicted and removed until convergence is obtained to obtain the depth field distribution map and the concentration field distribution map. The fifth unit is used to generate a three-dimensional spatial distribution model of gas leakage based on the depth field distribution map and the concentration field distribution map.

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

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

Images