Dam leakage fusion method based on physical constraint collaborative cross-modal attention

Abstract

The application discloses a dam leakage fusion method based on physical constraint collaborative cross-modal attention, and relates to the technical field of image data processing and multi-modal image fusion. The method acquires GPR root mean square amplitude images, GPR instantaneous phase images and ERT apparent resistivity images, and extracts a main direction of apparent resistivity gradient; after being decomposed and filtered by a dual-tree complex wavelet, high-frequency subbands in a seepage direction and a cross-seepage direction are weighted by applying anisotropic weights; a data-driven and physically-guided dual-path channel attention mechanism is constructed to complete low-frequency and high-frequency multi-scale feature fusion; three types of constraints, i.e., depth reliability, seepage-electricity coupling and seepage path continuity, are introduced to construct a total loss function, and after the fusion weight is optimized, a fusion attribute image is output through wavelet inverse transformation. The method can effectively enhance the seepage directionality feature, suppress single-modal artifacts, and improve dam leakage hidden danger detection accuracy and result reliability.

Description

Technical Field

[0001] This invention relates to the field of image data processing and multimodal image fusion technology, and in particular to a dam seepage fusion method based on physical constraint collaborative cross-modal attention. Background Technology

[0002] Accurate identification of seepage hazards in earth-rock dams is a core element in ensuring the safe operation of water conservancy projects. Ground-penetrating radar (GPR), based on the reflection characteristics of high-frequency electromagnetic waves, possesses high shallow detection resolution and can effectively capture the characteristics of fine geological structures and small seepage channels. High-density electrical resistivity tomography (ERT), based on differences in geological resistivity, achieves deep detection, is sensitive to water bodies, and has strong penetrating power. The integration of these two technologies can achieve complementary advantages in shallow and deep detection, representing an important technological development direction in the field of dam seepage detection.

[0003] Existing multimodal geophysical image fusion techniques have been extensively studied: multi-scale transform methods can achieve multi-resolution image decomposition and reconstruction, but they generally apply equal weights to sub-bands in all directions, failing to consider the inherent characteristic of seepage channels extending along the hydraulic gradient direction, which easily leads to the weak features of small seepages being submerged by the background stratigraphic structure; attention-based fusion methods can adaptively allocate feature fusion weights, but weight calculations largely rely on statistical features of the data, failing to embed seepage-related physical property criteria as prior information into the fusion process, and single-modal anomalies can easily cause fusion artifacts, reducing the accuracy of seepage identification; fusion schemes that introduce physical constraints mostly remain at the macroscopic geological structure level, failing to delve into the microscopic mechanism of seepage and electrical coupling, and lacking specificity in constraining seepage areas. Furthermore, existing multimodal fusion frameworks for dam safety are mostly decision-level fusions, lacking deep feature fusion designs at the GPR and ERT attribute image level, making it difficult to fully leverage the complementary advantages of multi-source data. Summary of the Invention

[0004] To address the aforementioned problems, this invention proposes a dam seepage fusion method based on physical constraint-based collaborative cross-modal attention.

[0005] To achieve the above objectives, the present invention employs the following technical solution:

[0006] A fusion method for dam seepage based on physical constraint-based collaborative cross-modal attention includes:

[0007] S1. Acquire geophysical property images, including GPR root mean square amplitude image, GPR instantaneous phase image and ERT apparent resistivity image; calculate the apparent resistivity gradient vector field based on the ERT apparent resistivity image and extract the principal direction of the apparent resistivity gradient.

[0008] S2. The geophysical property image is decomposed using dual-tree complex wavelet transform to obtain low-frequency approximate components and high-frequency detail sub-bands in multiple directions. Based on the main direction of apparent resistivity gradient, sub-bands along the seepage direction and sub-bands across the seepage direction are selected from the high-frequency detail sub-bands, and anisotropic weights are assigned to obtain the high-frequency detail components.

[0009] S3. Based on the low-frequency approximation component, the initial fused low-frequency component is obtained by channel attention weighted summation; the high-frequency detail component is mapped to a multi-channel feature map, the cosine similarity between cross-modal channel features is calculated and normalized to generate data-driven channel weights; the GPR root mean square amplitude normalized contrast and ERT apparent resistivity normalized contrast are calculated, and physical enhancement coefficients are defined; the enhanced feature map is calculated based on the physical enhancement coefficients and the multi-channel feature map to generate physical guided channel weights; based on the data-driven channel weights and the physical guided channel weights, the multi-channel feature map and the enhanced feature map are weighted and mixed to obtain the initial fused high-frequency component.

[0010] S4. Construct depth confidence constraints, seepage-electric coupling constraints, and seepage path continuity constraints, calculate the physical constraint loss, obtain the data fidelity term based on the initial fused low-frequency component, the initial fused high-frequency component, and the geophysical attribute image, and construct the total loss function in combination with the physical constraint loss; iteratively solve the total loss function as the optimization objective, obtain the low-frequency fused component and the high-frequency fused component, and then perform a dual-tree complex wavelet inverse transform to obtain the fused attribute image.

[0011] As a preferred embodiment of the present invention, the step of calculating the apparent resistivity gradient vector field based on the ERT apparent resistivity image and extracting the principal direction of the apparent resistivity gradient includes:

[0012] The expression for calculating the apparent resistivity gradient vector field is:

[0013]

[0014] In the formula, The coordinates are in the horizontal direction. The coordinates are in the depth direction. For the apparent resistivity gradient vector field, For ERT apparent resistivity image in ERT apparent resistivity at the location, The partial derivative in the horizontal direction, This is the partial derivative in the depth direction;

[0015] The horizontal pixel values ​​of the ERT apparent resistivity image are statistically averaged for each depth to obtain the baseline value of apparent resistivity at the corresponding depth. The reduction of ERT apparent resistivity relative to the baseline value of apparent resistivity at the same depth is calculated for each pixel. Pixel regions with a reduction exceeding a preset threshold for abnormal apparent resistivity amplitude are defined as low-resistivity abnormal regions. Within the low-resistivity abnormal regions, the main direction of the apparent resistivity gradient is extracted based on the apparent resistivity gradient vector field.

[0016] As a preferred embodiment of the present invention, the geophysical attribute image is decomposed using dual-tree complex wavelet transform to obtain low-frequency approximate components and high-frequency detail sub-bands in multiple directions; based on the principal direction of the apparent resistivity gradient, sub-bands along the seepage direction and sub-bands across the seepage direction are selected from the high-frequency detail sub-bands, and anisotropic weights are assigned to them respectively to obtain high-frequency detail components, including:

[0017] The GPR root mean square amplitude image, GPR instantaneous phase image and ERT apparent resistivity image are decomposed by dual-tree complex wavelet to obtain low-frequency approximate components and high-frequency detail subbands in multiple directions.

[0018] From the high-frequency detail subbands, the subband with the smallest angle between the center angle and the main direction of the apparent resistivity gradient is selected as the subband along the leakage direction, and the subband with the center angle orthogonal to the subband along the leakage direction is selected as the subband across the leakage direction.

[0019] The magnitude of a local anomaly is defined by the expression:

[0020]

[0021] In the formula, This represents a localized abnormal amplitude. Depth coordinates The apparent resistivity baseline value of the formation, For ERT apparent resistivity image in ERT apparent resistivity at the location;

[0022] Based on the local anomaly amplitude, the anisotropy weights of the sub-zones along and across the leakage direction are calculated separately, as expressed by:

[0023]

[0024] In the formula, The anisotropy weights of the sub-bands along the leakage direction are... For the anisotropy weights of the sub-bands across the leakage direction, These are the anisotropy weighting coefficients;

[0025] The subband along the leakage direction and the subband across the leakage direction are weighted pixel by pixel using corresponding anisotropic weights. The weighted subband along the leakage direction and the subband across the leakage direction together constitute the high-frequency detail component.

[0026] As a preferred embodiment of the present invention, the step of obtaining the initial fused low-frequency component based on the low-frequency approximation component through channel attention weighted summation includes:

[0027] Global average pooling is performed on the low-frequency approximate components to obtain the channel description vector. After nonlinear transformation by the hidden layer and ReLU activation function, the initial weight vector is generated by the output layer and Softmax activation function. The initial weight vector includes: GPR root mean square amplitude mode fusion weights. GPR instantaneous phase mode fusion weights and ERT apparent resistivity mode fusion weights ;

[0028] Based on the initial weight vector and the low-frequency approximation component, the initial fused low-frequency component is obtained by substituting them into the low-frequency component weighted fusion formula; the low-frequency component weighted fusion formula is as follows:

[0029]

[0030] In the formula, For the initial fusion of low-frequency components, The low-frequency approximation component of the root mean square amplitude of the GPR. This is the low-frequency approximation component of the instantaneous phase of the GPR. This is the low-frequency approximation component of the apparent resistivity of ERT.

[0031] As a preferred embodiment of the present invention, the step of mapping high-frequency detail components to multi-channel feature maps, calculating and normalizing the cosine similarity between cross-modal channel features, and generating data-driven channel weights includes:

[0032] High-frequency detail components are passed through 1 1. Convolution maps to a multi-channel feature space with C channels, resulting in a multi-channel feature map; the cosine similarity between the feature vectors of each channel in the multi-channel feature maps of any two modalities is calculated, expressed as:

[0033]

[0034] In the formula, For cosine similarity, For modality No. Channel feature vectors For modality No. Channel feature vectors and These are the L2 norms of the corresponding channel feature vectors;

[0035] All cosine similarities are combined to obtain a cosine similarity matrix. Softmax normalization is then performed row-by-row on the cosine similarity matrix to obtain a cross-modal channel attention map composed of cosine similarity attention values. For each channel, the cosine similarity attention values ​​between that channel and all channels of other modalities are extracted from the cross-modal channel attention map, and the average is taken. After nonlinear transformation by the hidden layer and ReLU activation function, the data-driven channel weights for each channel of the corresponding modality are generated by the output layer and Softmax activation function.

[0036] As a preferred embodiment of the present invention, the calculation of the GPR root mean square amplitude normalized contrast and ERT apparent resistivity normalized contrast, and the definition of physical enhancement coefficients; the calculation of enhancement feature maps based on physical enhancement coefficients and multi-channel feature maps, and the generation of physical guidance channel weights, including:

[0037] The root mean square amplitude baseline value is obtained by depth-wise statistical analysis of the low-frequency approximation component of the GPR amplitude, and the apparent resistivity baseline value is obtained by depth-wise statistical analysis of the low-frequency approximation component of the ERT apparent resistivity. Combining the GPR root mean square amplitude image and the ERT apparent resistivity image, the normalized contrast of the GPR amplitude and the normalized contrast of the ERT apparent resistivity are calculated. The expressions are as follows:

[0038]

[0039] In the formula, The GPR root mean square amplitude normalized contrast. ERT apparent resistivity normalized contrast. The root mean square amplitude of GPR is the measured value. The root mean square amplitude baseline value, This is the baseline value for apparent resistivity;

[0040] The physical enhancement factor is defined based on the GPR root mean square amplitude normalized contrast and the ERT apparent resistivity normalized contrast. The expression is:

[0041]

[0042] The physical enhancement coefficients are multiplied by the multi-channel feature map to obtain the enhanced feature map. The enhanced feature map is then subjected to global average pooling, followed by a nonlinear transformation through the hidden layer and the ReLU activation function. Finally, the physical guidance channel weights for each channel of the corresponding modality are generated through the output layer and the Softmax activation function.

[0043] As a preferred embodiment of the present invention, the step of weighted mixing of multi-channel feature maps and enhanced feature maps based on data-driven channel weights and physically guided channel weights to obtain initial fused high-frequency components includes:

[0044] Based on data-driven channel weights, multi-channel feature maps, physically guided channel weights, and enhanced feature maps, the initial fused high-frequency components are calculated using a dual-path weighted fusion formula for high-frequency components. The dual-path weighted fusion formula for high-frequency components is as follows:

[0045]

[0046] In the formula, For the initial fusion of high-frequency components, For data-driven path weight coefficients, For direction index, For the decomposition level index, For modal indexing, For channel index, Data-driven channel weights For multi-channel feature maps, For physical guidance channel weights, To enhance the feature map, These are the GPR root mean square amplitude image, the GPR instantaneous phase image, and the ERT apparent resistivity image, respectively. This represents the total number of channels.

[0047] As a preferred embodiment of the present invention, the construction of depth credibility constraints, seepage-electric coupling constraints, and seepage path continuity constraints includes:

[0048] The normalized detector sensitivity is calculated based on GPR and ERT. A depth confidence constraint loss is constructed using the deviations between the depth amplitude distributions of the initial fused low-frequency and high-frequency components and their corresponding detector sensitivity benchmarks. The expression for calculating the detector sensitivity benchmark is as follows:

[0049]

[0050] In the formula, For depth The detection sensitivity benchmark value at that location, This represents the GPR-normalized detection sensitivity. ERT normalized detection sensitivity The effective electromagnetic wave attenuation coefficient of the medium;

[0051] The ERT apparent resistivity is substituted into Archie's formula to obtain the water saturation of the ERT side of the formation. The relative permittivity of the soil is obtained by inversion based on the root mean square amplitude image and instantaneous phase image of GPR. The volumetric water content is calculated by the Topp formula and combined with the formation porosity to obtain the water saturation of the GPR side of the formation. The absolute difference between the water saturation of the ERT side and the water saturation of the GPR side of the formation is taken as the seepage-electric coupling constraint loss.

[0052] Seepage path continuity constraint loss for:

[0053]

[0054] In the formula, For a set of valid connected components, For the first The principal axis direction of each connected component. This represents the principal direction of the apparent resistivity gradient. For connected component index, This represents the total area of ​​the abnormal region. For the first The area of ​​each connected region.

[0055] In a preferred embodiment of the present invention, the physical constraint loss is calculated, and a data fidelity term is obtained based on the initial fused low-frequency component, the initial fused high-frequency component, and the geophysical attribute image. A total loss function is constructed by combining the physical constraint loss. The total loss function is iteratively solved with the optimization objective to obtain the low-frequency fused component and the high-frequency fused component, including:

[0056] Physical constraint loss for:

[0057]

[0058] In the formula, , and These are the loss weights for depth confidence constraints, seepage-electric coupling constraints, and seepage path continuity constraints, respectively. For depth-based credibility constraint loss, For seepage-electric coupling constraint loss, Loss due to seepage path continuity constraints;

[0059] The initial fused low-frequency component and the initial fused high-frequency component are subjected to dual-tree complex wavelet inverse transform to obtain the initial fused attribute image. The mean square error between the initial fused attribute image and the geophysical attribute image is used as the data fidelity term.

[0060] The total loss function is obtained by adding the data fidelity term and the physical constraint loss. With the goal of minimizing the total loss function, the parameters of the hidden layer and the output layer are optimized by gradient descent to obtain the optimized low-frequency component fusion weights and high-frequency component fusion weights. The optimized low-frequency component fusion weights are substituted into the low-frequency component weighted fusion formula to obtain the low-frequency fusion components. The optimized high-frequency component fusion weights are substituted into the high-frequency component dual-path weighted fusion formula to obtain the high-frequency fusion components corresponding to each layer and each direction.

[0061] As a preferred embodiment of the present invention, it further includes:

[0062] A comprehensive evaluation index is constructed based on the fused attribute image. The expression is:

[0063]

[0064] In the formula, , , , , and These are the weights of the evaluation indicators. For information entropy, For the average gradient, For structural similarity, The root mean square error, To improve the leakage detection and location rate, Cross-modal consistency retention rate;

[0065] The particle swarm optimization algorithm is adopted to maximize the comprehensive evaluation index and to constrain the leakage detection and localization rate to be no less than the preset localization rate threshold. The weight coefficients of the hidden layer, data-driven path, anisotropy weighting coefficient, physical constraint loss weight, and evaluation index weight are adaptively optimized. Based on the optimal hyperparameter combination obtained by adaptive optimization, steps S2 to S4 are re-executed to calculate the final fused attribute image.

[0066] The beneficial effects of this invention are as follows: First, through multi-scale anisotropic decomposition guided by the leakage direction, high-frequency sub-bands along and across the leakage direction are screened based on the principal direction of the ERT apparent resistivity gradient. Adaptive anisotropic weights are applied in combination with local anomaly amplitudes to specifically enhance the linear extension characteristics of leakage channels and suppress background interference such as stratification of formations across the leakage direction. This effectively avoids the problem of weak features of small leakages being submerged by the background under traditional dual-tree complex wavelet equal-weight decomposition, thus improving the detection capability of shallow and fine leakage channels. Second, a data-driven and physically guided dual-path cross-modal channel attention fusion mechanism is constructed. The data-driven path mines feature associations through cross-modal channel cosine similarity, while the physically guided path embeds the normalized contrast prior of GPR root mean square amplitude and ERT apparent resistivity. The dual-path weighted mixing generates fusion weights, which not only makes full use of the statistical complementarity of multimodal data but also effectively suppresses fusion artifacts caused by single-modal anomalies, improving the physical rationality of weight allocation and leakage identification accuracy. Finally, a multi-dimensional seepage physical constraint loss system was established, integrating depth reliability constraints, seepage-electrical coupling constraints, and seepage path continuity constraints. A total loss function was constructed using data fidelity terms to optimize the fusion parameters, forcing the fusion results to conform to geophysical detection laws and seepage mechanics mechanisms. This approach balances the fidelity of the original data with physical consistency, significantly improving the reliability of the fused attribute image and the identification accuracy of deep seepage areas. This method achieves complementary advantages between ground-penetrating radar and high-density electrical resistivity tomography (EDS), providing a reliable attribute image fusion scheme for the accurate detection of seepage hazards in earth-rock dams. Attached Figure Description

[0067] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the description of the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort. Wherein:

[0068] Figure 1 This is a flowchart of the method in an embodiment of the present invention;

[0069] Figure 2 This is a schematic diagram of the joint constraint process in an embodiment of the present invention. Detailed Implementation

[0070] 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, not all, of the embodiments of the present invention. All other embodiments obtained by those skilled in the art based on the described embodiments of the present invention are within the scope of protection of the present invention.

[0071] like Figures 1 to 2 As shown, this is an embodiment of the present invention, which provides a dam seepage fusion method based on physical constraint collaborative cross-modal attention, including:

[0072] S1. Obtain geophysical property images, including GPR root mean square amplitude image, GPR instantaneous phase image and ERT apparent resistivity image; calculate the apparent resistivity gradient vector field based on the ERT apparent resistivity image and extract the principal direction of the apparent resistivity gradient.

[0073] Specifically, for the target dam detection section, ground-penetrating radar (GPR) data and high-density electrical resistivity tomography (ERT) data were collected. After data preprocessing, attribute extraction, and spatial coordinate registration, three two-dimensional attribute images of consistent size and corresponding coordinate systems were obtained: the GPR root-mean-square amplitude image, the GPR instantaneous phase image, and the ERT apparent resistivity image. The GPR root-mean-square amplitude image is used to characterize the differences in electromagnetic wave reflection intensity of the formation medium; the GPR instantaneous phase image is used to characterize the location information of formation interfaces and structural boundaries; and the ERT apparent resistivity image is used to characterize the spatial distribution characteristics of the apparent resistivity of the formation medium. Based on the ERT apparent resistivity image, the central finite difference method was used to calculate the partial derivatives of the ERT apparent resistivity in the horizontal and depth directions pixel by pixel, constructing the apparent resistivity gradient vector field of the entire detection section.

[0074] In this embodiment, the expression for calculating the apparent resistivity gradient vector field is:

[0075]

[0076] In the formula, The coordinates are in the horizontal direction. The coordinates are in the depth direction. For the apparent resistivity gradient vector field, For ERT apparent resistivity image in The apparent resistivity at the ERT point was obtained by acquiring potential difference and current data between electrodes in the field using high-density electrical resistivity tomography (EPT). This data was then processed using least-squares inversion or the Occam inversion algorithm to obtain a two-dimensional apparent resistivity distribution matrix, which was then mapped into an image format. The coordinates were then... The corresponding pixel value is ; The partial derivative in the horizontal direction reflects the rate of change of the apparent resistivity of ERT in the horizontal direction; This is the partial derivative in the depth direction, reflecting the rate of change of apparent resistivity ERT in the depth direction.

[0077] The horizontal pixel values ​​of the ERT apparent resistivity image are statistically averaged for each depth to obtain the baseline value of apparent resistivity at the corresponding depth. The reduction of ERT apparent resistivity relative to the baseline value of apparent resistivity at the same depth is calculated for each pixel. Pixel regions with a reduction exceeding a preset threshold for abnormal apparent resistivity amplitude are defined as low-resistivity abnormal regions. The main direction of apparent resistivity gradient is extracted based on the apparent resistivity gradient vector field within the low-resistivity abnormal region.

[0078] Specifically, the ERT apparent resistivity image is traversed row by row along the depth direction. The arithmetic mean of the ERT apparent resistivity values ​​of all horizontal pixels at the same depth is calculated to obtain the baseline value of apparent resistivity corresponding to that depth. The decrease in ERT apparent resistivity at the current position relative to the baseline value of apparent resistivity at the same depth is calculated pixel by pixel. Pixel areas where the decrease exceeds a preset threshold for abnormal apparent resistivity amplitude are designated as low-resistivity abnormal areas. The preset threshold for abnormal apparent resistivity amplitude is adapted and adjusted according to the lithology of the dam construction area, the parameters of the detection equipment, and engineering detection experience, and the value ranges from 0.2 to 0.5. In this embodiment, it is set to 0.3.

[0079] Within the low-resistivity anomaly region, extract all coordinates within that region. The corresponding apparent resistivity gradient vector field is combined into a two-dimensional vector sample set. The sample covariance matrix of this two-dimensional vector sample set is calculated, and then eigenvalue decomposition is performed on the sample covariance matrix to obtain all eigenvalues ​​and their corresponding eigenvectors. The direction of the eigenvector matching the largest eigenvalue is selected as the main direction of the apparent resistivity gradient in the low resistivity anomaly region.

[0080] S2. The geophysical attribute image is decomposed using dual-tree complex wavelet transform to obtain low-frequency approximate components and high-frequency detail sub-bands in multiple directions. Based on the main direction of the apparent resistivity gradient, the sub-bands along the leakage direction and the sub-bands across the leakage direction are selected from the high-frequency detail sub-bands and anisotropic weights are assigned to obtain the high-frequency detail components.

[0081] In this embodiment, dual-tree complex wavelets are used to decompose the GPR root mean square amplitude image, GPR instantaneous phase image, and ERT apparent resistivity image to obtain low-frequency approximate components and high-frequency detail subbands in multiple directions.

[0082] Specifically, for the GPR root mean square amplitude image, GPR instantaneous phase image, and ERT apparent resistivity image, a dual-tree complex wavelet transform is used for multi-level decomposition. In this embodiment, the decomposition has four levels. Each level outputs one low-frequency approximation component and six high-frequency detail sub-bands in fixed directions. The center angles of the six fixed directions are ±15°, ±45°, and ±75°, respectively. The low-frequency approximation component represents the large-scale background contour information of the corresponding image, while the high-frequency detail sub-bands represent the image edges, textures, and abrupt structural information in the corresponding directions. The low-frequency approximation component includes the GPR root mean square amplitude low-frequency approximation component, the GPR instantaneous phase low-frequency approximation component, and the ERT apparent resistivity low-frequency approximation component.

[0083] From the high-frequency detail subbands, the subband with the smallest angle between its center angle and the principal direction of the apparent resistivity gradient is selected as the subband along the leakage direction, and the subband with its center angle orthogonal to the subband along the leakage direction is selected as the subband across the leakage direction.

[0084] Specifically, based on the extracted principal direction of the apparent resistivity gradient, the angle between the center angle of each high-frequency detail subband under the current decomposed layer and the principal direction of the apparent resistivity gradient is calculated one by one. The high-frequency detail subband with the smallest angle is selected as the percolation-direction subband, used to characterize the linear structural features along the extension direction of the percolation channel. The high-frequency detail subband with a center angle orthogonal to the percolation-direction subband is selected as the cross-percolation-direction subband, used to characterize the stratigraphic interface features perpendicular to the extension direction of the percolation channel. High-frequency detail subbands that are not selected are not included in the calculation.

[0085] Based on the depth baseline statistics of the ERT apparent resistivity image, the relative deviation of the ERT apparent resistivity at the current position relative to the background at the same depth is calculated pixel by pixel, defined as the local anomaly amplitude, and expressed as:

[0086]

[0087] In the formula, This represents a localized abnormal amplitude. Depth coordinates The apparent resistivity baseline value of the formation.

[0088] The larger the magnitude of the local anomaly, the greater the deviation of the ERT apparent resistivity at that location from the depth background, and the higher the confidence level of leakage occurrence.

[0089] Based on the local anomaly amplitude, the anisotropy weights of the sub-zones along and across the leakage direction are calculated separately, as expressed by:

[0090]

[0091] In the formula, The anisotropy weights of the sub-bands along the leakage direction are... For the anisotropy weights of the sub-bands across the leakage direction, The anisotropy weighting coefficient is a hyperparameter that controls the anisotropy enhancement intensity. Its typical value range is 0.3 to 0.8, and in this embodiment, it is 0.5.

[0092] The subband along the leakage direction and the subband across the leakage direction are weighted pixel by pixel using corresponding anisotropic weights. The weighted subband along the leakage direction and the subband across the leakage direction together constitute the high-frequency detail component.

[0093] Specifically, the anisotropic weights of the subband along the leakage direction and the corresponding subband along the leakage direction are multiplied pixel by pixel, and the anisotropic weights of the subband across the leakage direction and the corresponding subband across the leakage direction are multiplied pixel by pixel; the weighted subband along the leakage direction and the subband across the leakage direction together constitute the high-frequency detail components of this decomposition layer.

[0094] S3. Based on the low-frequency approximation component, the initial fused low-frequency component is obtained by channel attention weighted summation; the high-frequency detail component is mapped to a multi-channel feature map, the cosine similarity between cross-modal channel features is calculated and normalized to generate data-driven channel weights; the GPR root mean square amplitude normalized contrast and ERT apparent resistivity normalized contrast are calculated, and physical enhancement coefficients are defined; the enhanced feature map is calculated based on the physical enhancement coefficients and the multi-channel feature map to generate physical guided channel weights; based on the data-driven channel weights and the physical guided channel weights, the multi-channel feature map and the enhanced feature map are weighted and mixed to obtain the initial fused high-frequency component.

[0095] In this embodiment, global average pooling is performed on the low-frequency approximate components to obtain the channel description vector. After nonlinear transformation by the hidden layer and ReLU activation function, the initial weight vector is generated by the output layer and Softmax activation function. The initial weight vector includes: GPR root mean square amplitude mode fusion weights. GPR instantaneous phase mode fusion weights and ERT apparent resistivity mode fusion weights .

[0096] Specifically, the low-frequency approximation components of the GPR root mean square amplitude, GPR instantaneous phase, and ERT apparent resistivity are each treated as independent feature channels. For each independent feature channel, a global average pooling operation is performed on the two-dimensional low-frequency approximation component matrix to obtain the arithmetic mean of all pixel values ​​in that channel, thus obtaining the channel global statistics corresponding to that mode. The channel global statistics corresponding to the GPR root mean square amplitude mode, GPR instantaneous phase mode, and ERT apparent resistivity mode are sequentially concatenated to generate a channel description vector with a dimension length equal to the total number of modes.

[0097] The channel description vector is input into a two-layer fully connected network of the channel attention module for feature transformation. Specifically, the first layer is a hidden layer with 8-16 neurons, which performs nonlinear mapping using the ReLU activation function to learn the dependencies between different modal channels. The second layer is an output layer, which performs normalization using the Softmax activation function and outputs an initial weight vector whose values ​​sum to 1. The initial weight vector includes GPR root mean square amplitude mode fusion weights, GPR instantaneous phase mode fusion weights, and ERT apparent resistivity mode fusion weights.

[0098] Based on the initial weight vector and the low-frequency approximation component, the initial fused low-frequency component is obtained by substituting them into the low-frequency component weighted fusion formula; the low-frequency component weighted fusion formula is as follows:

[0099]

[0100] In the formula, For the initial fusion of low-frequency components, The low-frequency approximation component of the root mean square amplitude of the GPR. This is the low-frequency approximation component of the instantaneous phase of the GPR. This is the low-frequency approximation component of the apparent resistivity of ERT.

[0101] In this embodiment, the high-frequency detail components are transmitted through 1 1. Convolution maps to a multi-channel feature space with C channels, resulting in a multi-channel feature map; the cosine similarity between the feature vectors of each channel in the multi-channel feature maps of any two modalities is calculated, expressed as:

[0102]

[0103] In the formula, The cosine similarity value ranges from -1 to 1. A higher value indicates a stronger correlation in the feature space distribution of the two channels. For modality No. The channel feature vector is a two-dimensional feature map of the k-th channel in the multi-channel feature map of modality i. The two-dimensional feature matrix of H×W is flattened by rows or columns to obtain a one-dimensional channel feature vector of length H×W; where H is the feature map height and W is the feature map width. For modality No. Channel feature vectors and are the L2 norms of the corresponding channel feature vectors, and is the arithmetic square root of the sum of squares of all elements in the channel feature vector.

[0104] Specifically, the high-frequency detail components are transformed using 1×1 convolution to map the single-channel high-frequency components to a unified feature space with C channels, resulting in multi-channel feature maps for each modality. The 1×1 convolution operation maintains the spatial size of the feature map consistent with the input high-frequency detail components, performing feature dimensionality upscaling and spatial alignment only in the channel dimension. A typical value for the number of channels C is 64. Channel feature vectors are extracted channel by channel from the multi-channel feature map of each modality: the two-dimensional feature map of a single channel is flattened into a one-dimensional channel feature vector, which serves as the global feature representation for that channel.

[0105] All cosine similarities are combined to obtain a cosine similarity matrix. Softmax normalization is then performed row-by-row on the cosine similarity matrix to obtain a cross-modal channel attention map composed of cosine similarity attention values. For each channel, the cosine similarity attention values ​​between that channel and all channels of other modalities are extracted from the cross-modal channel attention map, and the average is taken. After nonlinear transformation by the hidden layer and ReLU activation function, the data-driven channel weights for each channel of the corresponding modality are generated by the output layer and Softmax activation function.

[0106] Specifically, the cosine similarities corresponding to all cross-modal channels are arranged and combined according to the modality index and the channel index to generate a cosine similarity matrix. Softmax normalization is performed on each row of the cosine similarity matrix to map the cosine similarity to the (0,1) interval and the sum of the elements in each row is 1, resulting in a cross-modal channel attention map composed of cosine similarity attention values, which quantitatively represents the feature association strength between channels of different modalities.

[0107] For each channel of each modality, the cosine similarity attention value corresponding to all channels of all other modalities is extracted from the cross-modal channel attention map, and the arithmetic mean is calculated to obtain the cross-modal association statistics of that channel. The cross-modal association statistics of all channels of the same modality are concatenated into a one-dimensional vector in channel order and input into a two-layer fully connected network for feature transformation. Specifically, the hidden layer and ReLU activation function are used to complete the non-linear mapping and learn the non-linear dependency relationship between channels. Then, the output layer and Softmax activation function are used to perform normalization processing and output the data-driven channel weights of each channel of the modality.

[0108] In this embodiment, the root mean square amplitude baseline value is obtained by depth-by-depth statistical analysis based on the low-frequency approximation component of the GPR root mean square amplitude, and the apparent resistivity baseline value is obtained by depth-by-depth statistical analysis based on the low-frequency approximation component of the ERT apparent resistivity. Combined with the GPR root mean square amplitude image and the ERT apparent resistivity image, the normalized contrast of the GPR root mean square amplitude and the normalized contrast of the ERT apparent resistivity are calculated.

[0109] Specifically, for the low-frequency approximation component of the root mean square amplitude of ground penetrating radar (GPR), the arithmetic mean of the amplitude values ​​of all horizontal pixels at the same depth is calculated by traversing the depth line one by one to obtain the baseline value of the root mean square amplitude at the corresponding depth, which represents the root mean square amplitude level under the background of normal strata; for the low-frequency approximation component of the apparent resistivity of high-density electrical resistivity (ERT), the same depth-by-depth mean statistical method is used to obtain the baseline value of the apparent resistivity at the corresponding depth.

[0110] By combining the GPR root mean square amplitude image and the ERT apparent resistivity image, the coordinates of the corresponding image are traversed pixel by pixel. The ratios of the measured GPR root mean square amplitude at the current pixel to the baseline GPR amplitude at the same depth, and the ratio of the apparent resistivity (ERT) at the current pixel to the baseline apparent resistivity at the same depth, are calculated respectively. These ratios are then used to generate pixel-by-pixel GPR root mean square amplitude normalized contrast and ERT apparent resistivity normalized contrast, expressed as follows:

[0111]

[0112] In the formula, The GPR root mean square amplitude normalized contrast. ERT apparent resistivity normalized contrast. The root mean square amplitude of the ground-penetrating radar (GPR) is measured by preprocessing, attribute extraction, and spatial registration of the raw data collected in the field. The resulting image is a two-dimensional image, and the pixel value at the corresponding coordinate is the root mean square amplitude of the GPR. The root mean square amplitude baseline value, This is the baseline value for apparent resistivity.

[0113] The physical enhancement factor is defined based on the GPR root mean square amplitude normalized contrast and the ERT apparent resistivity normalized contrast. The expression is:

[0114]

[0115] The physical enhancement coefficients are multiplied by the multi-channel feature map to obtain the enhanced feature map. The enhanced feature map is then subjected to global average pooling, followed by a nonlinear transformation through the hidden layer and the ReLU activation function. Finally, the physical guidance channel weights for each channel of the corresponding modality are generated through the output layer and the Softmax activation function.

[0116] Specifically, the pixel-wise physical enhancement coefficients are multiplied element-wise with the corresponding modality's multi-channel feature map. The physical enhancement coefficients are broadcast along the channel dimension to all channels. The same spatial weighting is applied to the two-dimensional feature map of each channel to obtain the physically enhanced multi-channel feature map, i.e., the enhanced feature map. Global average pooling is performed on the enhanced feature map channel by channel to obtain the global statistical values ​​of each channel, which are then concatenated into a channel description vector. This vector is then input into a two-layer fully connected network for feature transformation. Specifically, it first undergoes a non-linear mapping through a hidden layer and a ReLU activation function, and then undergoes normalization processing through an output layer and a Softmax activation function to generate the physical guidance channel weights for each channel of the corresponding modality.

[0117] In this embodiment, based on data-driven channel weights, multi-channel feature maps, physically guided channel weights, and enhanced feature maps, the initial fused high-frequency components are calculated using a high-frequency component dual-path weighted fusion formula; the high-frequency component dual-path weighted fusion formula is as follows:

[0118]

[0119] In the formula, For the initial fusion of high-frequency components, This is the data-driven path weight coefficient, with a value range of [0,1]. The larger the value, the more the fusion result is biased towards data statistical characteristics, and the smaller the value, the more it is biased towards physical prior constraints. A typical value is 0.6. For direction index, For the decomposition level index, For modal indexing, For channel index, Data-driven channel weights For multi-channel feature maps, For physical guidance channel weights, To enhance the feature map, These are the GPR root mean square amplitude image, the GPR instantaneous phase image, and the ERT apparent resistivity image, respectively. This represents the total number of channels.

[0120] Specifically, for multi-channel feature maps, a channel-by-channel spatial weighting operation is performed based on the data-driven channel weights of the corresponding modality. Specifically, the two-dimensional feature map corresponding to each channel of the modality is taken, and the data-driven channel weights matching the channel are multiplied pixel by pixel. All two-dimensional feature maps that have completed the weighting operation within the same modality are summed pixel by pixel along the channel dimension to generate a single-channel weighted feature map of the data-driven path corresponding to the modality.

[0121] Reusing the above-mentioned operational logic of channel-by-channel weighting and channel dimension summation, the enhanced feature map of this mode is weighted channel by channel using the weight of the physical guidance channel as the weighting coefficient. After all channels are weighted, the summation is carried out along the channel dimension to generate a single-channel weighted feature map of the physical guidance path corresponding to this mode.

[0122] A data-driven path weight coefficient is introduced as a balance coefficient for the output features of the two types of paths. The single-channel weighted feature map of the data-driven path and the single-channel weighted feature map of the physical guidance path are multiplied by the corresponding allocation coefficient and summed pixel by pixel to complete the linear fusion of dual-path features under a single mode. Dual-path fusion operation is performed sequentially on the GPR root mean square amplitude mode, GPR instantaneous phase mode, and ERT apparent resistivity mode. The single-channel feature map after the fusion of the three modes is accumulated and summed pixel by pixel along the mode dimension, and the initial fused high-frequency component corresponding to the current decomposition layer number and the current direction index is output.

[0123] S4. Construct depth confidence constraints, seepage-electric coupling constraints, and seepage path continuity constraints, calculate the physical constraint loss, obtain the data fidelity term based on the initial fused low-frequency component, the initial fused high-frequency component, and the geophysical attribute image, and construct the total loss function in combination with the physical constraint loss; iteratively solve the total loss function as the optimization objective, obtain the low-frequency fused component and the high-frequency fused component, and then perform a dual-tree complex wavelet inverse transform to obtain the fused attribute image.

[0124] like Figure 2 The diagram shown illustrates the joint constraint process in step S4 of this embodiment; it includes parallel construction of depth credibility constraints, percolation-electric coupling constraints, and percolation path continuity constraints, as well as total loss and weight optimization; wherein the total loss and weight optimization involves constructing a total loss function and optimizing the fusion weights of low-frequency components and high-frequency components. The specific details are as follows:

[0125] (1) Constructing deep credibility constraints

[0126] The depth weight distribution function is constructed as follows: In this embodiment, a detection sensitivity benchmark value is calculated based on the normalized detection sensitivity of GPR and ERT. A depth confidence constraint loss is constructed using the deviation between the depth amplitude distributions of the initially fused low-frequency components and the initially fused high-frequency components and their corresponding detection sensitivity benchmark values. The detection sensitivity benchmark value is then calculated using the depth weight distribution function, the expression of which is:

[0127]

[0128] In the formula, For depth The detection sensitivity benchmark value at that location, To normalize the detection sensitivity of GPR, and to provide the depth sensitivity of ground penetrating radar. The relative detection capability at a given location is calculated using an exponential decay model, where depth is the independent variable. The normalized detection sensitivity of ERT is obtained by calculating the mean of the diagonal elements of the Jacobian matrix obtained by high-density electrical resistivity tomography (ERT) and normalizing it to the maximum value. The effective electromagnetic wave attenuation coefficient of the medium is a constant characterizing the electromagnetic wave attenuation characteristics of the dam soil. It is obtained through exponential fitting based on the measured ground-penetrating radar amplitude-depth curves in the work area; alternatively, empirical values ​​can be selected according to the soil type in the work area, with a typical range of values. In this embodiment, 0.12 is used. .

[0129] The calculation of depth confidence loss and the setting of constraint targets are as follows: the initial fusion low-frequency component and the initial fusion high-frequency components of all decomposition layers and all directions are used to form a fusion multi-scale coefficient set, and the arithmetic mean of the amplitude of all pixels in each depth row is calculated to generate a depth amplitude distribution curve. The difference between the depth amplitude distribution curve and the detection sensitivity benchmark value of the corresponding depth is calculated for each depth. The arithmetic mean of the squared differences of all depths is obtained to obtain the mean square error, which is used as the depth confidence constraint loss.

[0130] (2) Constructing a seepage-electric coupling constraint

[0131] In this embodiment, the apparent resistivity of ERT is substituted into Archie's formula as the formation resistivity to obtain the formation water saturation on the ERT side; the relative permittivity of the soil is obtained by inversion based on the root mean square amplitude image and instantaneous phase image of GPR, and the volumetric water content is calculated by the Topp formula and combined with the formation porosity to obtain the formation water saturation on the GPR side; the absolute difference between the formation water saturation on the ERT side and the formation water saturation on the GPR side is used as the seepage-electric coupling constraint loss;

[0132] The ETR water saturation calculation is as follows: the apparent resistivity of ETR is substituted into Archie's formula as the formation resistivity, and combined with the formation porosity, lithology coefficient, cementation index and saturation index of the work area, the inversion operation is performed pixel by pixel on the entire exploration section to obtain the formation water saturation on the ETR side; in this embodiment, the lithology coefficient is taken as 1.0, the cementation index is taken as 1.8, the saturation index is taken as 1.8, and the formation porosity is taken as 0.4.

[0133] The GPR water saturation estimation is as follows: the relative permittivity of the soil is obtained by inverting the root mean square amplitude image and the instantaneous phase image of the GPR, and the volume water content of each pixel is calculated by substituting it into the Topp formula. The water saturation of the formation on the GPR side is then calculated by combining the formation porosity.

[0134] The seepage-electric coupling loss construction and penalty constraint are specifically as follows: the absolute difference of formation water saturation between the ERT side and the GPR side is calculated pixel by pixel, and the arithmetic mean of the absolute differences of all pixels in the entire detection section is taken as the seepage-electric coupling constraint loss. This constraint loss uses the consistency of the physical properties of the formation water saturation in the dual modes as the penalty target. The greater the seepage-electric coupling deviation of the fusion result, the higher the loss value, thereby forcing the fusion result to conform to the physical laws of the coupling between seepage and electrical properties of the dam, effectively suppressing false anomalies caused by contradictions in physical properties.

[0135] (3) Constructing seepage path continuity constraints

[0136] Specifically, the first step is to perform connected component filtering, which involves: segmenting the initial fused attribute image according to the bimodal leakage anomaly determination threshold, extracting all pixels that meet the bimodal leakage anomaly determination threshold to form a leakage anomaly region; using an eight-neighbor connected component labeling algorithm to group and label the pixels within the leakage anomaly region, resulting in several independent connected components; and removing scattered small connected components whose total pixel area is less than a preset minimum area threshold, leaving all remaining connected components to form a valid connected component set. In this embodiment, the bimodal leakage anomaly determination threshold includes an amplitude threshold and a resistivity threshold, where the amplitude threshold is 0.7 and the resistivity threshold is 0.7; the preset minimum area threshold is 50 pixels.

[0137] Next, weighted connectivity loss calculation is performed. Specifically, for each connected component in the effective connected component set, the two-dimensional coordinate data of all pixels it covers are extracted to construct a coordinate sample set. The sample covariance matrix of the coordinate sample set is calculated, and then eigenvalue decomposition is performed on the sample covariance matrix to obtain the corresponding eigenvalues ​​and eigenvectors. The direction of the eigenvector corresponding to the largest eigenvalue is selected as the principal axis direction of the connected component, which is used to characterize the dominant spatial extension trend of the connected component.

[0138] Finally, the seepage path continuity constraint rules were set, specifically as follows: Using the principal direction of the apparent resistivity gradient as the leakage reference direction, the cosine of the angle between the principal axis direction of each connected region and the leakage reference direction was calculated, serving as the directional consistency weight for that connected region; using the area of ​​each connected region as a weighting coefficient, the directional consistency weights of all connected regions were weighted and summed, and the ratio of the weighted sum to the total area of ​​the leakage anomaly region was used as the directional continuity index; based on this directional continuity index, the seepage path continuity constraint loss was constructed; the seepage path continuity constraint loss... for:

[0139]

[0140] In the formula, A set of valid connected components; For the first The principal axis direction of each connected component; This represents the principal direction of the apparent resistivity gradient. For connected component index, The total area of ​​the abnormal region is denoted as , and the total number of all abnormal pixels is denoted as . For the first The area of ​​a connected region is calculated by counting the number of pixels within the corresponding connected region and converting the area to the actual size of the pixels.

[0141] (4) Total Loss and Weight Optimization

[0142] In this embodiment, the weighted calculation of the joint physical constraint loss specifically involves: weighting and summing the depth confidence constraint loss, the seepage-electric coupling constraint loss, and the seepage path continuity constraint loss to obtain the physical constraint loss. for:

[0143]

[0144] In the formula, , and These are the loss weights for depth confidence constraints, seepage-electric coupling constraints, and seepage path continuity constraints, respectively. For depth-based credibility constraint loss, For seepage-electric coupling constraint loss, This represents the loss due to the continuity constraint of the seepage path.

[0145] in To constrain the loss weight based on the depth of confidence, an empirical value is initially set based on the range of detection depth in the work area, typically ranging from 0.1 to 0.5, and in this embodiment, it is set to 0.3; The weight for the seepage-electric coupling constraint loss is initially set based on empirical values ​​according to the lithology and water-bearing conditions of the work area, with a typical range of 0.2 to 0.6. In this embodiment, 0.4 is used. The loss weight for the continuity constraint of the seepage path is initially set with an empirical value based on the development characteristics of the seepage channel. The typical value range is 0.1 to 0.4, and in this embodiment, it is set to 0.3.

[0146] The initial fused low-frequency components and initial fused high-frequency components are subjected to dual-tree complex wavelet inverse transform to obtain the initial fused attribute image. The mean square error between the initial fused attribute image and the geophysical attribute image is used as the data fidelity term. The data fidelity term and the physical constraint loss are added to obtain the total loss function. With the goal of minimizing the total loss function, the parameters of the hidden layer and the output layer are optimized by gradient descent to obtain the optimized low-frequency component fusion weights and high-frequency component fusion weights. The optimized low-frequency component fusion weights are substituted into the low-frequency component weighted fusion formula to obtain the low-frequency fused components. The optimized high-frequency component fusion weights are substituted into the high-frequency component dual-path weighted fusion formula to obtain the high-frequency fused components corresponding to each layer and each direction.

[0147] Specifically, the data fidelity term is calculated as follows: The initial fused low-frequency component and the initial fused high-frequency components from all decomposition layers and directions are subjected to a dual-tree complex wavelet inverse transform to reconstruct the initial fused attribute image; the mean square error between the initial fused attribute image and the geophysical attribute image is calculated, and the average value is taken as the data fidelity term. The data fidelity term is used to constrain the fusion result from deviating from the basic information of the original detection data, avoiding excessive correction of physical constraints that could lead to the loss of original effective information.

[0148] The total loss function is constructed by linearly adding the data fidelity term and the physical constraint loss. This total loss function balances the fidelity of the original data with the consistency of seepage physics, achieving synergistic guidance from statistical data characteristics and prior physical knowledge.

[0149] In this embodiment, the method for weight optimization is parameter gradient descent optimization, specifically: with minimizing the total loss function as the optimization objective, gradient descent is used to iteratively update the parameters of the hidden and output layers of the trainable network for all channels of the attention module; during the iteration, the gradient value of the total loss function with respect to the parameters of each layer is solved using the backpropagation algorithm, and the parameters of each layer are updated successively according to a preset learning rate, which is 0.005 in this embodiment; the iteration continues until the value of the total loss function converges to a preset convergence threshold, which is 0.005 in this embodiment. The process terminates when the number of iterations reaches the preset maximum number of iterations. In this embodiment, the preset maximum number of iterations is 50 rounds, and the optimized low-frequency component fusion weight and high-frequency component fusion weight are output.

[0150] Substitute the optimized low-frequency component fusion weights into the low-frequency component weighted fusion formula, and perform pixel-by-pixel weighted summation on the low-frequency approximate components to obtain the low-frequency fusion components; substitute the optimized high-frequency component fusion weights into the high-frequency component dual-path weighted fusion formula, and perform layer-by-layer and direction-by-direction calculations according to the decomposition layer index and direction index to obtain the high-frequency fusion components corresponding to the seepage direction and the cross-seepage direction under each decomposition layer.

[0151] To meet the input format requirements of the dual-tree complex wavelet inverse transform, a multi-scale coefficient set with dimension matching is constructed, and the fusion coefficients are organized layer by layer according to the hierarchical order and direction definition corresponding to the dual-tree complex wavelet decomposition. The specific implementation method is as follows:

[0152] The low-frequency fusion components are used as the low-frequency approximations corresponding to the maximum decomposition level and are assigned to the corresponding level of the multi-scale coefficient set. For each decomposition level, the high-frequency fusion components along the seepage direction and across the seepage direction are placed at the high-frequency detail sub-band positions corresponding to the center angles during the decomposition process, corresponding one-to-one with the directional sub-bands in the decomposition stage. For the remaining four high-frequency detail sub-bands in each decomposition level other than those along the seepage direction and across the seepage direction, a two-dimensional coefficient matrix with all zeros and the same pixel size as the high-frequency detail sub-bands of the current decomposition level is used to fill them, ensuring that the number of high-frequency directional sub-bands in each decomposition level matches the inherent number of six directional sub-bands in dual-tree complex wavelet decomposition, which conforms to the input dimension specification of dual-tree complex wavelet inverse transform. Using the constructed multi-scale coefficient set as input, the inverse wavelet transform is performed layer by layer from the highest decomposition level to the lowest decomposition level according to the standard reconstruction algorithm of dual-tree complex wavelet inverse transform, finally outputting a two-dimensional fusion attribute image in the spatial domain.

[0153] The fused attribute image output by the inverse transform is normalized to map the image pixel values ​​to the same numerical range as the original geophysical attribute image. At the same time, the pixel size and coordinate system of the fused attribute image are verified to be completely aligned with the input geophysical attribute image to ensure the spatial positioning accuracy of subsequent leakage anomaly interpretation.

[0154] In this embodiment, a comprehensive evaluation index is constructed based on the fused attribute image. The expression is:

[0155]

[0156] In the formula, , , , , and These are the weights of the evaluation indicators. In this embodiment, 0.1 0.15 0.15 0.1 0.3 It is 0.2; Information entropy is used to quantify the information richness of fused attribute images. It is calculated by statistically analyzing the occurrence probability of each gray level based on the image gray-level histogram and substituting it into the information entropy formula. The higher the entropy value, the more effective information the image contains. This is the average gradient, used to quantify the sharpness of edge details in the fused attribute image. The gradient magnitude in the horizontal and depth directions is calculated pixel-by-pixel, and the average value across the entire image is taken. A higher gradient value indicates richer image texture details and sharper boundaries. Structural similarity is used to quantify the consistency of structural information between the fused attribute image and the geophysical attribute image, with a value range of [0,1]. The structural similarity between the fused attribute image and the geophysical attribute image is calculated separately, and the average value is taken. A higher value indicates better preservation of structural information. The root mean square error is used to quantify the numerical deviation between the fused attribute image and the geophysical attribute image. The root mean square error of the fused attribute image and the geophysical attribute image are calculated separately, and the square root of the average value is obtained. The lower the value, the higher the data fidelity. The leakage detection and localization rate is used to quantify the ability of the fused attribute image to capture the real leakage area. Based on the dual anomaly criteria of the GPR root mean square amplitude normalized contrast being lower than a preset amplitude threshold and the ERT apparent resistivity normalized contrast being lower than a preset resistivity threshold, the detection leakage area is delineated in the fused attribute image. The intersection-union ratio of the detection leakage area and the real leakage marked area is used as the leakage detection and localization rate. The higher the value, the higher the leakage localization accuracy. The cross-modal consistency retention rate is used to quantify the ability of the fusion process to maintain the balance of the advantageous features of the two modes. The structural similarity between the fused attribute image and the ground-penetrating radar attribute image, and between the fused attribute image and the high-density electrical resistivity tomography (HRT) attribute image, is calculated separately. The ratio of the smaller to the larger of these two values ​​is taken as the cross-modal consistency retention rate. The degree of retention of ground-penetrating radar features and high-density HRT features in the fused attribute image is statistically analyzed, and the normalized result of their ratio is calculated. The closer the value is to 1, the better the cross-modal balance.

[0157] The particle swarm optimization algorithm is adopted to maximize the comprehensive evaluation index and to constrain the leakage detection and localization rate to be no less than the preset localization rate threshold. The weight coefficients of the hidden layer, data-driven path, anisotropy weighting coefficient, physical constraint loss weight, and evaluation index weight are adaptively optimized. Based on the optimal hyperparameter combination obtained by adaptive optimization, steps S2 to S4 are re-executed to calculate the final fused attribute image.

[0158] Specifically, the particle population size, maximum number of iterations, inertia weight, and learning factor are set as follows: the particle population size is 20~30, the maximum number of iterations is 20~30, the inertia weight is 0.7, the individual learning factor is 1.5, and the global learning factor is 1.5. The position vector of each particle corresponds to a set of hyperparameters to be optimized. The upper and lower bounds of the values ​​are preset according to the physical meaning and algorithm constraints of each hyperparameter to be optimized. The initial position and initial velocity of each particle are randomly initialized within the corresponding value range.

[0159] For each hyperparameter combination corresponding to a particle, the fusion process from step S2 to step S4 is fully executed to obtain the corresponding fused attribute image; the comprehensive evaluation index value of the fused attribute image is calculated as the fitness of the particle; at the same time, it is verified whether the leakage detection location rate is not lower than the preset location rate threshold. In this embodiment, the preset location rate threshold is 0.85; particles that do not meet the constraints are given extremely low fitness to ensure that the optimization results always meet the leakage detection accuracy requirements.

[0160] The current fitness of each particle is compared with the historical best fitness of other particles to update the individual's optimal position; the optimal fitness of all particles is compared to update the global optimal position. According to the standard update formula of the particle swarm optimization algorithm, the velocity and position of each particle are updated by combining the individual optimal and global optimal positions, and the hyperparameters are constrained within a standard preset range to avoid invalid parameters. The fitness calculation and particle update steps are repeated until the maximum number of iterations is reached or the global optimal fitness converges. The final global optimal position vector is the optimal hyperparameter combination.

[0161] Substitute the optimal hyperparameter combination obtained by optimization into the fusion process, and re-execute the dual-tree complex wavelet decomposition and anisotropic weighting in step S2, the dual-path cross-modal attention fusion in step S3, and the physical constraint loss optimization and dual-tree complex wavelet inverse transform in step S4 to output the final fused attribute image.

[0162] In summary, this invention first employs a dual-tree complex wavelet anisotropic decomposition guided by the seepage direction. Based on the principal direction of the high-density electrical resistivity gradient, it filters high-frequency subbands along and across the seepage direction. Adaptive weights are applied in conjunction with local anomaly amplitudes to specifically enhance the linear characteristics of seepage channels, suppress background interference such as formation bedding, and improve the detection capability of shallow, subtle seepage channels. Secondly, a data-driven and physically-guided dual-path cross-modal channel attention mechanism is constructed. The data-driven path mines statistical correlations of cross-modal channel features, while the physically-guided path embeds normalized contrast priors of ground-penetrating radar amplitude and high-density electrical resistivity, effectively suppressing fusion artifacts caused by single-modal anomalies and improving the physical rationality of weight allocation. Then, three types of physical constraints are established: depth credibility, seepage-electrical coupling, and seepage path continuity. A total loss function is constructed using data fidelity terms to optimize fusion parameters, forcing the fusion results to conform to geophysical exploration laws and seepage mechanics mechanisms, thereby improving the deep identification accuracy and physical property consistency of the fused attribute image. Finally, a two-layer optimization architecture of inner-layer gradient descent and outer-layer particle swarm optimization is adopted. Combined with the exclusive comprehensive evaluation index for leakage detection, key hyperparameters are adaptively tuned to improve the method's adaptability to different dam construction areas and improve the overall detection accuracy and reliability of dam leakage hazards.

[0163] In the above embodiments, implementation can be achieved, in whole or in part, through software, hardware, firmware, or any other combination thereof. When implemented in software, it can be implemented, in whole or in part, as a computer program product, which includes one or more computer instructions. When the computer program instructions are loaded and executed on a computer, all or part of the processes or functions according to this application are generated. The computer can be a general-purpose computer, a special-purpose computer, a computer network, or other programmable device. The computer instructions can be stored in a computer-readable storage medium or transferred from one computer-readable storage medium to another.

[0164] Furthermore, the functional units in the various embodiments of this application can be integrated into a processing module, or each unit can exist physically separately, or two or more units can be integrated into a module. The integrated module can be implemented in hardware or as a software functional module. If the integrated module is implemented as a software functional module and sold or used as an independent product, it can also be stored in a computer-readable storage medium. This storage medium can be a read-only memory, a disk, or an optical disk, etc.

[0165] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any person skilled in the art can easily conceive of various variations or substitutions within the technical scope disclosed in this application, and these should all be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

Claims

1. A fusion method for dam seepage based on physical constraint collaborative cross-modal attention, characterized in that, include: S1. Obtain geophysical property images, including GPR root mean square amplitude image, GPR instantaneous phase image and ERT apparent resistivity image; The apparent resistivity gradient vector field is calculated based on the ERT apparent resistivity image, and the principal direction of the apparent resistivity gradient is extracted. S2. The geophysical property image is decomposed using dual-tree complex wavelet transform to obtain low-frequency approximate components and high-frequency detail sub-bands in multiple directions. Based on the main direction of apparent resistivity gradient, sub-bands along the seepage direction and sub-bands across the seepage direction are selected from the high-frequency detail sub-bands, and anisotropic weights are assigned to obtain the high-frequency detail components. S3. Based on the low-frequency approximation component, the initial fused low-frequency component is obtained by channel attention weighted summation; High-frequency detail components are mapped to multi-channel feature maps, and the cosine similarity between cross-modal channel features is calculated and normalized to generate data-driven channel weights. Calculate the root mean square amplitude normalized contrast ratio (GPR) and the apparent resistivity normalized contrast ratio (ERT), and define the physical enhancement factor. The enhanced feature map is calculated based on the physical enhancement coefficient and the multi-channel feature map, and the physical guidance channel weight is generated. The initial fused high-frequency component is obtained by weighted mixing of the multi-channel feature map and the enhanced feature map based on the data-driven channel weight and the physical guidance channel weight. S4. Construct depth confidence constraints, seepage-electric coupling constraints, and seepage path continuity constraints, calculate the physical constraint loss, obtain the data fidelity term based on the initial fused low-frequency component, the initial fused high-frequency component, and the geophysical attribute image, and construct the total loss function in combination with the physical constraint loss; iteratively solve the total loss function as the optimization objective, obtain the low-frequency fused component and the high-frequency fused component, and then perform a dual-tree complex wavelet inverse transform to obtain the fused attribute image.

2. The dam seepage fusion method based on physical constraint collaborative cross-modal attention according to claim 1, characterized in that, The calculation of the apparent resistivity gradient vector field based on the ERT apparent resistivity image and the extraction of the principal direction of the apparent resistivity gradient include: The expression for calculating the apparent resistivity gradient vector field is: ； In the formula, The coordinates are in the horizontal direction. The coordinates are in the depth direction. For the apparent resistivity gradient vector field, For ERT apparent resistivity image in ERT apparent resistivity at the location, The partial derivative in the horizontal direction, This is the partial derivative in the depth direction; The horizontal pixel values ​​of the ERT apparent resistivity image are statistically averaged for each depth to obtain the baseline value of apparent resistivity at the corresponding depth. The reduction of ERT apparent resistivity relative to the baseline value of apparent resistivity at the same depth is calculated for each pixel. Pixel regions with a reduction exceeding a preset threshold for abnormal apparent resistivity amplitude are defined as low-resistivity abnormal regions. The main direction of apparent resistivity gradient is extracted based on the apparent resistivity gradient vector field within the low-resistivity abnormal region.

3. The dam seepage fusion method based on physical constraint collaborative cross-modal attention according to claim 1, characterized in that, The method employs dual-tree complex wavelet transform to decompose the geophysical attribute image, obtaining low-frequency approximate components and high-frequency detail sub-bands in multiple directions. Based on the principal direction of the apparent resistivity gradient, sub-bands along the seepage direction and sub-bands across the seepage direction are selected from the high-frequency detail sub-bands, and anisotropic weights are assigned to them respectively to obtain the high-frequency detail components, including: The GPR root mean square amplitude image, GPR instantaneous phase image and ERT apparent resistivity image are decomposed by dual-tree complex wavelet to obtain low-frequency approximate components and high-frequency detail subbands in multiple directions. From the high-frequency detail subbands, the subband with the smallest angle between the center angle and the main direction of the apparent resistivity gradient is selected as the subband along the leakage direction, and the subband with the center angle orthogonal to the subband along the leakage direction is selected as the subband across the leakage direction. The magnitude of a local anomaly is defined by the expression: ； In the formula, This represents a localized abnormal amplitude. Depth coordinates The apparent resistivity baseline value of the formation, For ERT apparent resistivity image in ERT apparent resistivity at the location; Based on the local anomaly amplitude, the anisotropy weights of the sub-zones along and across the leakage direction are calculated separately, as expressed by: ； In the formula, The anisotropy weights of the sub-bands along the leakage direction are... For the anisotropy weights of the sub-bands across the leakage direction, These are the anisotropy weighting coefficients; The subband along the leakage direction and the subband across the leakage direction are weighted pixel by pixel using corresponding anisotropic weights. The weighted subband along the leakage direction and the subband across the leakage direction together constitute the high-frequency detail component.

4. The dam seepage fusion method based on physical constraint collaborative cross-modal attention according to claim 3, characterized in that, The initial fused low-frequency component, obtained by channel attention weighted summation based on the low-frequency approximation component, includes: Global average pooling is performed on the low-frequency approximate components to obtain the channel description vector. After nonlinear transformation by the hidden layer and ReLU activation function, the initial weight vector is generated by the output layer and Softmax activation function. The initial weight vector includes: GPR root mean square amplitude mode fusion weights. GPR instantaneous phase mode fusion weights and ERT apparent resistivity mode fusion weights ; Based on the initial weight vector and the low-frequency approximation component, the initial fused low-frequency component is obtained by substituting them into the low-frequency component weighted fusion formula; the low-frequency component weighted fusion formula is as follows: ； In the formula, For the initial fusion of low-frequency components, The low-frequency approximation component of the root mean square amplitude of the GPR. This is the low-frequency approximation component of the instantaneous phase of the GPR. This is the low-frequency approximation component of the apparent resistivity of ERT.

5. The dam seepage fusion method based on physical constraint collaborative cross-modal attention according to claim 4, characterized in that, The process of mapping high-frequency detail components to multi-channel feature maps, calculating and normalizing the cosine similarity between cross-modal channel features, and generating data-driven channel weights includes: High-frequency detail components are passed through 1 1. Convolution maps to a multi-channel feature space with C channels, resulting in a multi-channel feature map; the cosine similarity between the feature vectors of each channel in the multi-channel feature maps of any two modalities is calculated, expressed as: ； In the formula, For cosine similarity, For modality No. Channel feature vectors For modality No. Channel feature vectors and These are the L2 norms of the corresponding channel feature vectors; All cosine similarities are combined to obtain a cosine similarity matrix. Softmax normalization is then performed row-by-row on the cosine similarity matrix to obtain a cross-modal channel attention map composed of cosine similarity attention values. For each channel, the cosine similarity attention values ​​between that channel and all channels of other modalities are extracted from the cross-modal channel attention map, and the average is taken. After nonlinear transformation by the hidden layer and ReLU activation function, the data-driven channel weights for each channel of the corresponding modality are generated by the output layer and Softmax activation function.

6. The dam seepage fusion method based on physical constraint collaborative cross-modal attention according to claim 5, characterized in that, The calculation of GPR root mean square amplitude normalized contrast and ERT apparent resistivity normalized contrast is described, and the physical enhancement coefficient is defined. An enhanced feature map is calculated based on the physical enhancement coefficient and multi-channel feature map, and physical guidance channel weights are generated, including: The root mean square amplitude baseline value is obtained by depth-wise statistical analysis of the low-frequency approximation component of the GPR amplitude, and the apparent resistivity baseline value is obtained by depth-wise statistical analysis of the low-frequency approximation component of the ERT apparent resistivity. Combining the GPR root mean square amplitude image and the ERT apparent resistivity image, the normalized contrast of the GPR amplitude and the normalized contrast of the ERT apparent resistivity are calculated. The expressions are as follows: ； In the formula, The GPR root mean square amplitude normalized contrast. ERT apparent resistivity normalized contrast. The root mean square amplitude of GPR is the measured value. The root mean square amplitude baseline value, This is the baseline value for apparent resistivity; The physical enhancement factor is defined based on the GPR root mean square amplitude normalized contrast and the ERT apparent resistivity normalized contrast. The expression is: ； The physical enhancement coefficients are multiplied by the multi-channel feature map to obtain the enhanced feature map. The enhanced feature map is then subjected to global average pooling, followed by a nonlinear transformation through the hidden layer and the ReLU activation function. Finally, the physical guidance channel weights for each channel of the corresponding modality are generated through the output layer and the Softmax activation function.

7. The dam seepage fusion method based on physical constraint collaborative cross-modal attention according to claim 6, characterized in that, The initial fused high-frequency component is obtained by weighted mixing of multi-channel feature maps and enhanced feature maps based on data-driven channel weights and physically guided channel weights, including: Based on data-driven channel weights, multi-channel feature maps, physically guided channel weights, and enhanced feature maps, the initial fused high-frequency components are calculated using a dual-path weighted fusion formula for high-frequency components. The dual-path weighted fusion formula for high-frequency components is as follows: ； In the formula, For the initial fusion of high-frequency components, For data-driven path weight coefficients, For direction index, For the decomposition level index, For modal indexing, For channel index, Data-driven channel weights For multi-channel feature maps, For physical guidance channel weights, To enhance the feature map, These are the GPR root mean square amplitude image, the GPR instantaneous phase image, and the ERT apparent resistivity image, respectively. This represents the total number of channels.

8. The dam seepage fusion method based on physical constraint collaborative cross-modal attention according to claim 1, characterized in that, The construction of depth credibility constraints, seepage-electric coupling constraints, and seepage path continuity constraints includes: The normalized detector sensitivity is calculated based on GPR and ERT. A depth confidence constraint loss is constructed using the deviations between the depth amplitude distributions of the initial fused low-frequency and high-frequency components and their corresponding detector sensitivity benchmarks. The expression for calculating the detector sensitivity benchmark is as follows: ； In the formula, For depth The detection sensitivity benchmark value at that location, This represents the GPR-normalized detection sensitivity. ERT normalized detection sensitivity The effective electromagnetic wave attenuation coefficient of the medium; The ERT apparent resistivity is substituted into Archie's formula to obtain the water saturation of the ERT side of the formation. The relative permittivity of the soil is obtained by inversion based on the root mean square amplitude image and instantaneous phase image of GPR. The volumetric water content is calculated by the Topp formula and combined with the formation porosity to obtain the water saturation of the GPR side of the formation. The absolute difference between the water saturation of the ERT side and the water saturation of the GPR side of the formation is taken as the seepage-electric coupling constraint loss. Seepage path continuity constraint loss for: ； In the formula, For a set of valid connected components, For the first The principal axis direction of each connected component. This represents the principal direction of the apparent resistivity gradient. For connected component index, This represents the total area of ​​the abnormal region. For the first The area of ​​each connected region.

9. The dam seepage fusion method based on physical constraint collaborative cross-modal attention according to claim 7, characterized in that, The calculation yields the physical constraint loss, and the data fidelity term is obtained based on the initial fused low-frequency component, the initial fused high-frequency component, and the geophysical attribute image. The total loss function is then constructed by combining the physical constraint loss. Iteratively solving the total loss function as the optimization objective yields the low-frequency fusion component and the high-frequency fusion component, including: Physical constraint loss for: ； In the formula, , and These are the loss weights for depth confidence constraints, seepage-electric coupling constraints, and seepage path continuity constraints, respectively. For depth-based credibility constraint loss, For seepage-electric coupling constraint loss, Loss due to seepage path continuity constraints; The initial fused low-frequency component and the initial fused high-frequency component are subjected to dual-tree complex wavelet inverse transform to obtain the initial fused attribute image. The mean square error between the initial fused attribute image and the geophysical attribute image is used as the data fidelity term. The total loss function is obtained by adding the data fidelity term and the physical constraint loss. With the goal of minimizing the total loss function, the parameters of the hidden layer and the output layer are optimized by gradient descent to obtain the optimized low-frequency component fusion weights and high-frequency component fusion weights. The optimized low-frequency component fusion weights are substituted into the low-frequency component weighted fusion formula to obtain the low-frequency fusion components. The optimized high-frequency component fusion weights are substituted into the high-frequency component dual-path weighted fusion formula to obtain the high-frequency fusion components corresponding to each layer and each direction.

10. The dam seepage fusion method based on physical constraint collaborative cross-modal attention according to claim 9, characterized in that, Also includes: A comprehensive evaluation index is constructed based on the fused attribute image. The expression is: ； In the formula, , , , , and These are the weights of the evaluation indicators. For information entropy, For the average gradient, For structural similarity, The root mean square error, To improve the leakage detection and location rate, To maintain cross-modal consistency, a particle swarm optimization algorithm is adopted. With the goal of maximizing the comprehensive evaluation index and the constraint that the leakage detection and localization rate is not lower than the preset localization rate threshold, adaptive optimization is performed on the weight coefficients of the hidden layer, data-driven path, anisotropy weighting coefficient, physical constraint loss weight, and evaluation index weight. Based on the optimal hyperparameter combination obtained by adaptive optimization, steps S2 to S4 are re-executed to calculate the final fused attribute image.

Images